{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2f9f2f21-569b-4429-8f5c-081e6089bc78",
   "metadata": {},
   "source": [
    "# Building up tools to compute an approximation of the 2-Wasserstein distance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b42cdee-d0de-4562-8669-06916312af00",
   "metadata": {},
   "source": [
    "> In this section we create a function to compute an approximation of the 2-Wasserstein distance between two univariate data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb9a10f4-12da-4456-a50c-ed52bf6b2e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "import ot\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from statFEM_analysis.oneDim import mean_assembler, kernMat, cov_assembler, sample_gp\n",
    "from scipy.stats import linregress\n",
    "from scipy import integrate\n",
    "from scipy.linalg import sqrtm\n",
    "from tqdm.notebook import tqdm\n",
    "import sympy; sympy.init_printing()\n",
    "# code for displaying matrices nicely\n",
    "def display_matrix(m):\n",
    "    display(sympy.Matrix(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d30ae0ca-9076-4be4-be54-63eaf72656c0",
   "metadata": {},
   "source": [
    "## Computing the 2-Wasserstein distance between two data-sets\n",
    "\n",
    "We start by creating a function [wass()](statFEM_analysis.rst#statFEM_analysis.maxDist.wass) to estimate the 2-Wasserstein distance between two data-sets `a` and `b`, using the Python package [POT](https://github.com/PythonOT/POT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "17b61beb-dd7d-4ec6-ba3a-d2326075eaf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statFEM_analysis.maxDist import wass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "781816a8-b1fd-4a1c-a56f-f5b9f1f25953",
   "metadata": {},
   "source": [
    "`wass` takes in the two datasets `a` and `b` as well as an argument `n_bin` which controls how many bins are used to create the histograms for the datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e6cb9d6-dd22-48ef-acd5-90edf9645c7a",
   "metadata": {},
   "source": [
    "Let's test this function out. First we make sure it gives $\\operatorname{wass}(a,a) = 0$ for any dataset $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f5409c91-138f-42ac-a0ea-c5fbeeee6aa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# standard normal\n",
    "N = 1000 # number of samples\n",
    "n_bins = 10 # number of bins\n",
    "np.random.seed(134)\n",
    "a = np.random.normal(size=N)\n",
    "assert wass(a,a,n_bins) == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd46535-f54c-4332-9187-e9323d32de47",
   "metadata": {},
   "source": [
    "We also test it on samples from 2 different Gaussians, $a\\sim\\mathcal{N}(m_0,s_0^{2})$ and $b\\sim\\mathcal{N}(m_1,s_1^{2})$. We expect, theoretically, that $\\operatorname{wass}(a,b)=\\sqrt{|m_0-m_1|^{2}+|s_0-s_1|^{2}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1a09e1d6-5fe0-49ee-b4e7-2f11d4fdcb79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up means and standard deviations\n",
    "m_0 = 7\n",
    "m_1 = 58\n",
    "s_0 = 1.63\n",
    "s_1 = 0.7\n",
    "\n",
    "# draw the samples\n",
    "N = 1000\n",
    "#####################################\n",
    "n_bins = 50 # number of bins\n",
    "#####################################\n",
    "np.random.seed(2321)\n",
    "a = np.random.normal(loc = m_0, scale = s_0,size=N)\n",
    "b = np.random.normal(loc = m_1, scale = s_1,size=N)\n",
    "\n",
    "# tolerance for the comparison\n",
    "tol = 1e-1\n",
    "\n",
    "# compute the 2-wasserstein with our function and also the true theoretical value\n",
    "W = wass(a,b,n_bins)\n",
    "W_true = np.sqrt(np.abs(m_0-m_1)**2 + np.abs(s_0-s_1)**2)\n",
    "# compare\n",
    "assert np.abs(W - W_true) <= tol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a20642-7942-4e9e-9e58-63e678a1073a",
   "metadata": {},
   "source": [
    "Let's take the previous example and compute the distance for a range of different means and standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e09dbc00-94f1-4bad-bc11-24c4111a8022",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up range for means and standard deviations\n",
    "n = 40\n",
    "m_range = np.linspace(m_0 - 2, m_0 + 2, n)\n",
    "s_range = np.linspace(s_0/4, 2*s_0, n)\n",
    "\n",
    "# set up arrays to hold results with our function, the theoretical results, \n",
    "# and theoretical results using estimated means and standard deviations\n",
    "W = np.zeros((n, n))\n",
    "W_0 = np.zeros((n, n))\n",
    "W_est = np.zeros((n,n))\n",
    "\n",
    "N = 10000 # number of samples\n",
    "################################################\n",
    "n_bins = 100 # number of bins\n",
    "################################################\n",
    "np.random.seed(2321)\n",
    "a = np.random.normal(loc = m_0, scale = s_0,size=N)\n",
    "m_a_est = np.mean(a)\n",
    "s_a_est = np.std(a)\n",
    "\n",
    "# sample for each m,s in the ranges and compute the results\n",
    "for i, m in enumerate(m_range):\n",
    "    for j, s in enumerate(s_range):\n",
    "        b = np.random.normal(loc = m, scale = s, size = N)\n",
    "        m_est = np.mean(b)\n",
    "        s_est = np.std(b)\n",
    "        \n",
    "        W[i,j] = wass(a,b,n_bins)\n",
    "        W_0[i,j] = np.sqrt(np.abs(m - m_0)**2 + np.abs(s - s_0)**2)\n",
    "        W_est[i,j] = np.sqrt(np.abs(m_est - m_a_est)**2 + np.abs(s_est - s_a_est)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e20cccc-a9b7-440c-adde-aec590edb288",
   "metadata": {},
   "source": [
    "Let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8d86ba25-6a18-4bfb-bb82-a19a07b7e818",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAFgCAYAAADAT84SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df/xkV13n+fen0+QHdmwaogMkaWgh9hiC/DBGRMfNDhgTFpPdFTVxRPGh5uGsjDK6ziLOI2pG3dWZcWYRViYSjLDKD1F5BDZI68qI7MiPiEnID1s7NHYSgkDSNmkT0qb7s39UVbq6un7cunXvPZ9zzuuZx33k+/1WddWpe895n3vO/VHm7gIAAAAAAEC+tqQuAAAAAAAAADbDBA8AAAAAAEDmmOABAAAAAADIHBM8AAAAAAAAmWOCBwAAAAAAIHNM8AAAAAAAAGSOCR4AACAz+2dmtjd1OeYxs4vN7N4E7/tpM3vp0O8LAADQBhM8KM54h/wRMztsZn9nZjeY2bbxYy83s4+Z2T+Y2QNm9ttmds74sdeN/81hM/uSmR2d+v2OtJ8KAOabybzJ8oYG/87N7NmT3939z9x9d09lvMHMfqGP1x6//gmfpeVr9FpGAOjaTO4fm+kL/kXq8gEYHhM8KNW3u/s2SS+UdKGkf2tmr5D0O5L+s6SzJD1H0qOSPmxmO9z9l9x92/jf/YikP5/87u7PSfQ5AKCJb5/Kq23u/urUBYrEzLamLgMAdG069yUd0Il9wW/PPp8sBMrHBA+K5u73SXq/pOdK+o+SfsHdf8fdH3H3z0r6IUmHJf3rhMUEgF6Y2bPN7E/N7JCZfcHM3jn++4fGT7l1fKT3u2cvgxqfGfRTZnbb+KzH683sn5jZ+83sITP7YzPbMfX83zWzz47f60Nm9pzx36+W9C8k/Zvxe713/Penm9nvmdnnzWy/mf3Y1GudMT6j5qCZ3Snp65d8xoWfxcz+NzP7rKTfNLNXmdmHZ/6tj9fR3DKOPX+8Dg6Z2TvN7PS1NwQAJGBmvzDOrbeb2UOSvtfM/m8z+7mp57zUzD499fs5ZvYHU9n8owmKDqAlJnhQNDM7V9LLJD0saaek351+3N2PSfo9Sd86fOkAoHf/TtIeSTsknSPp1yTJ3b9l/Pjzxkd637ng33+HRvn41ZK+XaMJ89dJ+gqN9iF+bOq575d0nqSvlPQJSb89fq/rxj//yvi9vt3Mtkh6r6RbJZ0t6SWSXmNm3zZ+rZ+V9Kzx8m2Svn/RB1zyWZ4q6cmSniHp6kX/flEZpx7+LkmXStol6WslvWrZawFAMP+TRmewb5e0KOslSeNsfp+kj2uUzd8q6afM7CV9FxJAN5jgQaneY2Z/L+nDkv5U0uvHf79/znPv1+iSLQDI1XvM7O+nlh8e//0fNZrgeLq7f8ndP7zkNeb5NXf/u/HZkH8m6aPu/pfu/iVJfyDpBZMnuvtb3P0hd39U0s9Jep6ZbV/wul8v6Svc/Vp3P+Lun5L0G5KuHD/+XZJ+0d0fdPd7dDzD13FM0s+6+6Pu/kiLfz/xenf/jLs/qNGk1PM3eC0AGNqH3f297n6sQRZ+o6QvH9+24Ii775N0vY5nM4DguA4Tpfof3f2PJ7+Y2T8d//g0Sftnnvs0SV8YqmAA0IMTMm/Kv9HoLJ6PmdlBSf/R3d+yxuv+3dTPj8z5fXID+1Mk/aKk79To7J5j4+ecJenQnNd9hqSnjyfiJ07RaBJJkp4u6Z6px/52jTJPfH48EbWpz079/LBGZQOAXNyz+imPe4aknXOy+b92WiIAvWGCB7XYK+lejQYfvzL54/hU1O+Q9J5E5QKA3ozvNfbDkmRm3yzpj83sQ+Ojsl36HklXSHqppE9rdCnAQUk2KcrM8++RtN/dz1vwevdLOlfS5BsMd7Yo0+x7/oOkJ05+MbOnrng+AJRgaRZqdDnrxD2S/sbdv6b3UgHoBZdooQru7pL+V42+Tet7zOz08c79myV9uaT/lLSAANADM/tOMztn/OtBjXb0J2fX/J2kr+rorc7U6FsJH9Bo4PBLM4/PvtfHJD00vgnyGWZ2ipldYGaTmym/S9JPm9mOcfn/1Yr3b/JZbpX0HDN7/vhGyT/X4jUAIHe3SPofxvn6NJ14L7U/l3TEzH5yvK98ipk918y+Lk1RAayLCR5UY3zjzVdq9I1ZD0i6U9IZkr7J3R9IWTYA2NB7x9/+NFn+YPz3r5f0UTM7LOlGST8+vt+NNJrg+K3xPXu+a8P3f6tGl1Hdp1G2fmTm8eslnT9+r/e4+1FJL9fofjb7NbpM9s0anfkjST8/fr39Gt0k+m0r3n/lZ3H3v5Z0raQ/lvQ3Gt2jbWEZV7wfAOTqBkl3aZSxfyjpHZMH3P0xjb6c5CKNzsb8gqT/otHBUAAZsNGJDQAAAAAAAMgVZ/AAAAAAAABkjgkeAAAAAACAAZnZW8zsc2Z2+4LHzcxeb2b7zOw2M3vhqtdkggcAAAAAAGBYN0i6dMnjl0k6b7xcLenXV70gEzwAAAAAAAADcvcPSXpwyVOukPRWH/mIpCeNv/1uoa1dFnAIO568xZ9+zubFfuDotg5Kk84Xj5yeugjAQl9+6pdSF2Gpp5xyuJfXvfOT//gFd/+KXl48APJ3hPxFZNHzV+ong8nfZnLMXzIXOcghe2d1lcWl5++3/fdf5g88eLTVv/2L2x69Q9J05bjO3a9b4yXOlnTP1O/3jv92/6J/kN0Ez9PP2arfed8/2fh13nbwxR2UJp09B3anLgKw0CU796Yuwkqv3PHfOn/N5z/j3r/t/EUDIX9H2fuVqQsBrBA9g8nf9dWav2QuchI9e2d1lcWl5+8DDx7Vxz6ws9W/PeVpf/Mld7+w4yItVe0lWn3sXAyFyR1El0MdzW0nFwCaip7Bbzv4YjI4gdzWefR6DMzKrc7mlgmpuKRjLf/rwH2Szp36/Zzx3xaqdoIHQL9y6OTo2LCOHOo0MJFDfSWDh5Pbus6h/gLz5FZ3c8uGNFxH/VirpQM3Svq+8bdpvUjSIXdfeHmWlOElWrXLLTRycnj/dknStl2HEpekHHsO7A5/uuqkY8v5rL7cvHLHf2OHAo+bZO8EGVyXtx18MfmLE7Cvu7nZXG2DLG4vh/1fNDc6g8d7eW0ze7ukiyWdZWb3SvpZSU+QJHd/k6SbJL1M0j5JD0v6gVWvyQQPinJ4//a1O6TZTpCJnm7l0skxyMAyDDiWm87Rptm5aABCBneH/IWU1xF6sra5LiZx2rw+2dxMLvkrkcFNdHS51Unc/aoVj7ukH13nNau+RCu3ikynt9jh/dsf74gmP6/q+FY9p8lroJlc6m5OO8FAaouydvrv8zK0abaSwd3IKX/J4LrlUldTmM3VlNkYqSzR5VSnyd9ycAYPstV0gDAxOeKwbkfU5qwgnCyXIxnTHVxuk8DoR047aH1rsyO/6c5/m7ODcKJc8lfistmu5TJoI2dPlNukCZfaLpZT/mI+l+uo93OJVh+Y4MlE7R1fFx3dJq/BaardyK2TY6CBmkUcYMwrEzncDPlbn1wmdzASMXPbYmL+RLnkL5dqLdbXPXj6UPUEDx1fXLl0chyxWF8unRz6k1P21ja5nkv2TttkIFHb/X7IX0RUW85OyzFz18Vkzwj5my+XdJQJnvhyGWDU1unl3tGt24nR6aE2uWSvVFf+5p69E00nbEr5vMAiuWRtTTk7rdYMYr83Ps7imY8zeIKj04ul1E5u0UBj2TfH1NLZcRQDkZG9+SN/FyN/EUUtWTtRcua2UdsZlFI++cskz4lc4h48kTG5E0ctHV0tnxNYJYf8JXvLss5nrWmSB2Uja2OpKXPbqG2iJ5dJHpyony9J70fVX5Me0Z4Du4vt9PhKxdVqWi+l1nPMx4AjHbK3uVrWT6l1HWRtJGTuempaXzm0gRyyZCgu19GWSwpVncETvaLm0NjXUUtId42jyMCwSsteifzFahxFRgol5u0s8ncztZzRQwajL5zBE0RJHV5NM/B9qWX9lVTvsVjkyfXS6iD5uznWH3IVOWtrQP52q4b1GX0fhEwZc+loyyWFaiZ4IlfQ6I17HaUH8ZBqWZcl1X/kpbS6V0tmDKGGgYVU9mXhiKfUulZLXqTC+k0r8hh6KK7RPXjaLClUMcETuWKW1NkRvt2r5d4Z0QcZkTMkuqjrLnJ9W1fp+ZBSTRkcVdQMiSb6eopcxzZRci5EU2oOl9o2ymI62nJJofgJnsgdXikNutTAjaj0wUbkiZ7IWRJV1HUWtY61UWIORFZyBpO/+Yq+fqLWq02UmAG5KHG9R28j0TOmby7pmLdbUih6gidyZYzekJsqMWRzUuoORtT2ETlToom6rqLWrTZKbPu5KTGDo7aRqJmC5aLWp7ZKbPM5KnE7lNZWSsMZPAFE3hEopQGXFqw5K7GjiypytkQRdR2VlL2091hK2x5R20rUbEmJdTKc0tp5Cdgmw6k5a1xM8CQXtQJGPv15XQRqTCVtl8htJWrGRBB13USuT+soqY2Xhom3YUTNGJyspNylbcdV0vYppc2U6JhbqyWF4iZ4onb8pTTYkkK0VCVto8jtJmrWpBR1nUSuR+sopV2XrpTtFLndvO3gi8PmzZAir4PI9WcdpbTnGpSyrSK3nciZg+OKmuCJWukiN9R1lBKctShle0VuP1EzJ4Wo6yJy/VlHKe25FqVsr1LaD7Cukg6W1aSU7Ub2xsIlWmNmdrqZfczMbjWzO8zs5+c85zQze6eZ7TOzj5rZM/sqDzZTQljWqJTtFrmjizixQf6WpZR2XJtSthv5u56h8jfiZ5+IXGeaKKXt1oxt2J/I2dMXl+motrRaUujzXR+V9M/d/XmSni/pUjN70cxzflDSQXd/tqT/JOmX275ZjZUNaKKUoxlYC/mr/AcZEjupuSv1K9Wx1KD5G03uuUtbLUfu2zL3tlQa7sEjyUcOj399wniZ/Tb4KyT91vjnd0t6iZmlWRM9KaFx5h6QGMl9oFFCWxrKkPkbdXKnBLm2VcyXcwaTv80Nkb9Rczf3epJj28RyuW/TqG0qagb1hUu0ppjZKWZ2i6TPSfojd//ozFPOlnSPJLn7Y5IOSXrKnNe52sxuNrObDz547KT3qa2SDSn3YMR8uQ40onZ0EQ2Vv1FRVxAZ+Vu22vM3R7m1RzSXY94iGtNR39JqSaHXd3X3o+7+fEnnSLrIzC5o+TrXufuF7n7hjiefWOTIkzu57wwRhnXIbTvn3q6GUnv+5i63dol2ctvOEfM3Yg71mb8RP68Us240weC/Hrlu56htK2oW9cElHdOWVksKg7yru/+9pA9KunTmofsknStJZrZV0nZJDwxRJiyXawiindy2d7TOLnInR/7mJ7f2iM3kNsCMlr+R1ZK/udaJnNodupHrNs+1jZWES7QkmdlXmNmTxj+fIelbJf3VzNNulPT9459fIelP3H32OuWFIg+qcpZr+GEzuW13OrvFas/fnOtGbu0Q3clp2+fcxvrWZ/5Gzt3c5NTe0C22fXdqySR3LtGaeJqkD5rZbZI+rtE1yO8zs2vN7PLxc66X9BQz2yfpJyS9tsfyDCrXnR9Cr265bf9c29kAqs7fXOXW/lA38nehqvI3x3pA1iLHOpBjW0MaW/t6YXe/TdIL5vz9mqmfvyTpO9u8fi0zhkPKMewAnIz8zQ/5C2lUD7btOpS6GNkZZdK7UhdDUn/5+8DRbZsXrmM5DjjJWkzkmLd7DuzWJTv3pi7GCSLlb5+OJbrcqo005w0hHDo8YDO1THrU8jmHRP4iVzkO8FEvshazqBNoYvQ16VtaLSlkOcHDAAPoT26dHQMMTMutPuTW3gBAImtRjtzqRm5trwzcg6dXEU9RnZVbw8st2ACkkUP+5oTsBYD+kbVYJbc6kttYM3d8TTqyklugAfPQ0SE3ZC9KQf7WJ6dtTtaiKeoKljnq1mpJgQmeSh3ev50gA4AEyF4sQ/1AZEzuoGQ51Zmc2mLuXMY9eBBbTuGFNHKsI3R0iF4HmFhHqaK3PdSHrEVb1B3Mc8y3tFpSYIKnQ3sO7A69k8PgAuugriAnkbNXoj1hPdQXRBQ9ZydoP9hULnUolzaJYTHB05HoDSyXoEIs1BvkIHL+MrGOtqg3wPpoN+hKLnUp8j5QKfia9ApFb1i5BBRiov4gssj5S9vBpqhDiCJy1k7QXtA16hSk8T14uMlyPaJ3eAQTukA9QkSR85c2g65Ql5Ba5KydoJ2gLznUrRzaaO74mvRKRG9MOQQS8pFDfYreJtGdyNs6h7aCvORQpyK3SZQth/aBvFHH6uYuHfUtrZYUspvgecoph1MXQVL8HRmCCH2gXtWN/F2NNoK+ULeQQuS8BYYUPYNpq30yHWu5pLA1ybtmLnoDih5AyNvh/du1bdeh1MVApSLnL9mLvpG/wInI3fnOvHuzY/gPPetYRyUB8udSsrNx2shugueBo9uSvj+DC4BBRq3I38XI3/YWDUQYYMwXOX/3HNitS3buTV0MdCRy5krkrrT5RM66r1t7LkfOX4kM7lOqb8RqI7sJHpystg7uzLu3dNbBNOkYa+/MFonayU12SOngMISa8nc6LzfJxaYDkmXPqz2Xo+avxAADw6gpe6f1NaHT9v1rzOLI+QtITPCsJeKRjNo6uEnHMvn/uh1Lm46x6b+hk4uFQUZZyN+0ZnNwnQzuY0DCICN+/kpMtOcsYuZO1JS9UvpJnWW6mvjPTfT8JXu75TIdS/SV520wwZOxmjq4RZ3bqo5lyE6xyzOLchK9k5MYZKB7teTvqgydl8EpBiNtJ/1zFzl/JQYawCYiT+zMU9tkT/T8RbdyukQrn5ImFvlIRsnOvHvLWqfzzy5DS/W+qUUf7NJ+8xZt+0Wv711ZN8si5F+EMgwten2M1n6xWuRtFr2+byrlPmyXSvgMOYvchnPkko75llZLCrS8TJXewUn5HbmYVmPHFr1O7jmwmw4PG4tez7tQQn7lXv51Ra+X5C+6EL2eb6KE3J2n1M81UXKdxDTT0ZZLCuW2uA5F2ympIUxK6QxK+RwlidaesVyk7UX25qX0gUWOIrVnzBd1G5Wav7XkVMmfM2rdjNqWc8QZPOhV1BDpUmkdQMmd2qwa6ieGw87JsErNqVoymPwF8lFLLs0q9XNHzV/2o7qT0xk83GR5CRrFsEoM/Gm13XwuMm78GR/5O6zS81cig6Mgf+OKmrtRB89t1JC1TdR6Y3zkyd2SnY3TRj4lHVjUTq5UtXV4pdxEb56SdsSQBvk7nFJzaJVSM5j8RWlKqtOl5U0XSsrhkuoqhmNml5rZXjPbZ2avnfP4TjP7oJn9pZndZmYvW/WaZbSojjG4GFYpwd5WSZ0bsInIN2ItcceN3BkhgwH2fftExqxWyvqJuK9A2+7GUd/SalnGzE6R9EZJl0k6X9JVZnb+zNP+raR3ufsLJF0p6f9aVdYyWlNHIg8uSkSHdyLWx7Bo67GwPYZF1pyMdTIc2juaiDhYXgf7dethfSEql3RM1mpZ4SJJ+9z9U+5+RNI7JF0x5+2/fPzzdkmfWfWitKKxHHY2cu/ophHgi7FuUJsc8rckZMxiDDBQo4gZnPs+LznSXu7rLmLdjdjG82K9nMEj6WxJ90z9fu/4b9N+TtL3mtm9km6S9K9WvWjeLagjVPph5R7cQ8h9kBGxc0NM5O9wcs+VIeW8nshfIB1ythu5r0dyuCyjr0m3Vouks8zs5qnl6jXf/ipJN7j7OZJeJultZra0cfAtWhhMzkGdypl3byniGwbmdXTbdh1KUJIT8W0uaTG5Mxzyd325f8vLsgEG+YuJiDmc6+CYnO1eKfvBEZC5mzna/ryYL7j7hQseu0/SuVO/nzP+27QflHSpJLn7n5vZ6ZLOkvS5RW9YfRJF7NhKRKfXXq5HMQ7v3/74AiCNXPMjkhzXH7mLXOVad3PMiVzkum5zrcs4mavd2TvjM3iW+bik88xsl5mdqtFNlG+cec4BSS+RJDP7GkmnS/r8shfNs8VUKNeQYHDRHdYlkEbO+YtukL8AFiEb+sc67gYnNrR3TFtaLcu4+2OSXi3pA5Lu0ujbsu4ws2vN7PLx035S0g+b2a2S3i7pVe7uy16XS7TQC4K4PyWcrhrh8gCgVORvf8hflCTaYC/HyXTydjg5XjZ7eP92MhdLuftNGt08efpv10z9fKekb1rnNZngyUBuHR6dXf9y7OSAHJG/mEX+AiBr08ltop1Jnvy5S0dXX24VBukUHIMLLMP6BvpD/mIZ1jdyxtk77dH202MbtBet7eeip3vw9ILWgc4QtmlwbwigezkNNiTyNxXyF9hcTnlLe48jp22RUx3HyUY3Wd7Sakkhn5bRk8hfF5dTGOQUsqViGwDdyCl7Jdp+BEz0ICccwW+HNh4P26QdMmB9R2WtlhRoFUHlNMAgXONgkLE+Orl0Ik+w54L2HgvbYz3kL3LZ36Vtx5XLtsmlruNkLi7Ryg6DjPZyCdXasF2QC/K3Pdp5TGwXoCy06fjYRutjgn0dXKIFIAA6u+bo5DCRyxE22ndsnE3ZHPk7rEjrO4e8pR3nI4dtlUOdx3zHZK2WFOK3hIFwFHl9OQRp7aJuIzo4TCN/1xO1XSMP5C/QDFmbH7YZ+jD5mvQ2Swq0gilRBhnRd744MpkXtlUzkY5qAovQnvPC9mqG/B1GpPWcw74u8hR920Wq+5EyAd2J3QIQTvTQxHwRt1ukDg7pRZlgjyxiO8ZqEbcb+QugZBFzdxoZnB/uwRMYM5XtRQ9LYFPkQ3/2HNjN+m2Jsybzx/ZbjXyoR/TBLe21DGxHdMXV7hu0ivsWLTM718w+aGZ3mtkdZvbjc55zsZkdMrNbxss1fZVHOr7zwEBjfYRk/iJuw+g7ebmKmL8T5O96IrZblIH87UfU/CV3myFzMZQoGUw2NJPTTZa39vjaj0n6SXf/hJmdKekvzOyP3P3Omef9mbu/vMdyLBSxQkdp7NPo7Mpx5t1b9NCzjqUuRmh7Duwu4XIh8reFaPlL9paF/F2N/C1ftJydRuaWh9xFF1xKdjZOG70lmbvf7+6fGP/8kKS7JJ3d1/s1EX2nIWKnR2eHvkWs97mLmL/RRauHZC+GEK3el4D8zROZWy627WoRD7pFwz14ZpjZMyW9QNJH5zz8jWZ2q5m938yes+DfX21mN5vZzQ8ffLTHkmIagVgmtmtdyN/80EbLxbatS5T8jTJ4Y0IRqUTNXtpEJlref6e4e/BMmNk2Sb8n6TXu/sWZhz8h6Rnu/jxJvybpPfNew92vc/cL3f3CJ+44rd8CJxKtgUcNQpQpWv0vRaT8jXwGZaT6R/ZiaJHqf0ki5S+WI3cBLOPK6x48vSaamT1Bo87tt93992cfd/cvuvvh8c83SXqCmZ3VZ5miObx/e7idKzq68rGNl4tytHMT5O9qEfMX5SN/lyN/yxQ1a2mP9Yi6raO0jRKyFyN9fouWSbpe0l3u/qsLnvPU8fNkZheNy/NAX2WKJOrAImr4oXwR20OuyN/lyF/gRBHbQ66i5S+DtsXI3PqwzdFWTpdo9fktWt8k6ZWSPmlmt4z/9jpJOyXJ3d8k6RWS/qWZPSbpEUlXurv3WKbkIu9EEXp14ZsFikb+zkH+Igryt2jk74zI2Yv6RMzfw/u3a9uuQ6mLgQVy+xat3iZ43P3D0vILz9z9DZLe0FcZoqGDA5ajg+sG+XuyyPnL5A4iIH+7Qf7mgdwFTrbnwO7Q90xMKacJHtINkujoahVxu0cZiHNaO4A+kb+Lkb/dYD3OF7HtYVgR60CU/MXJXHyLFoCMROzkomDnuDusy/lof3Vj+y9GZpSBQSuiIn+xjpy+RavPe/CEw87CfAQcol2PzKUCQHl27D2ig7tPTV2McMhfYDg17fPu2Huk1b8jp+vGZVpzeF6XaFU1wQPMdnZ0YscxyJiPjg59qWWgMZ278wYc5HDM/JWUPIPJ3/YiHNSMdvZOyZnbdjKn6WuVmtMRszd17uJk3GQZWSm5s5u2qONb1iGW2pktE7Gjk9IPMgCsr+mAo8nzasjjaPkrMdgAoutyYqfJ+9SQxUDu6hjdK8aRDHaShrdj75HWnd/k304vSCP1UcAI+ZEz1l9d+sjLWvI34kEX8hdtpK43syK2rbZS7peWuE8crW5EazsY4SbLQGJ9dT4ldWiLROvoJujwUJKo7WwTfe/015C/UZG/eWFS7ESl5G20iZVo5dlEKXWkK2TIiXL7Fi0u0RpQtB2k0sKM01S7E/FSAc6AyxM7CScje7t5r1IzOGL+Akgn+iRKDfvFgHMPnlgYYJyslAFGyk6v9IEGgwyge2Rvf2UoKYej5S8T7FhHpAOaOWduhJxdR+4TPZFyl/ufxZPqK8/bqGKCJwI6u25F6/S4SShwsiiT61Hyl+zt16qykcEAchE5a1fJfaIHI3yL4XGe2dek57+3uUKUAUYUuQ8wcr7eN8dy515fACZ3upFz9k7kVv7c6wzSYL/3uBzbUAlZO5Hj54hUZ6Lsv2DE3VotKRR9Bk+UTi5KA40UWuvKsZOYh6Ma7XCaan6i5G8EZG8cO/YeySp/I1wyQP5iHVH2eXNUWt5Kee73RshdRJPuhslt5LvXibUwwIglp8+Uc91B3RhobCannFpHSUfIAcyX075LDZlU+ucDIskn/dYU5ehxhAFGTp3ctNI7vJw+W651CGlEyd8Icmw7pWfvRC6fMcc6hDTI3vzkkkNdyOmzRsnd1ONIMuW4nC7RilF7C5W6UUpxAmodtQwupLo+KzAk8red2vKots8L9ClC7kr5ZG+N+VPjZ0b+XKObLLdZUsgjAdFabteQEvxx5bLD1BeOYqBkNU825/C5yV/yF2XJIXf6kstnrz13McVH36TVZkmBmosQah5cSPl0dqlEOSrIICMPUepLLsgfLBOlPZG/aCKHQTmZy37/OqJkcO2OyVotKcRPwUzRGJsj4EdyWA857Dj1jUEGSsEO9nE5rAfyl/xdJvW6Yb+3mRyyZkisj/hSZ0sELu7BAzTC4ALriLTzSGcXV6R6EhXZOx/rZLFI7Yr8xSLRJ0HJmPxEr1MYSrv773APnoJE2hGKik5uPtZLPhhkYJWI90UpdikAACAASURBVEAjY/LGYGNkz4HdZDCyQvYuxrpZjbFletyDB1iCIF8u+vphgHEcAwzkJHq2RMA6WiziAIMMjiFC3Yi8b0KurBZ5HUWuWxgOl2hVLEInF1nkAEd8EdsXA4w4ItaPKMheoFz0Q3GRvc2xruIiY/JS5AQPlfC4SLPOBHdzrKt8XLJzb+oihEL+Hhcpf9Ec+QugC2TJ+qKuswj9OQex0hldbsUZPNXatutQ6iI8LkIYoUwp6xYdHHJA/uYt6iBDIn+nMcGeXoQ6ETFvI2dIdKw7RMRNlgFsjA4uPgYXAJAG+QuUi31gRMNNlhPi8oCRiEczsD46uPkiHDHEycjfkWj5S460w3qbj/xFJNHyViI7ShWxrg2p9n08LtFCUrUHEDAEjh5jHvIX6B/5e1ztg65omNzpDuvyZEywp+FqN7nDBA+KRki3x7qLh8EFUIcde4+QwcGQv3Ew2DwRWdG9aOuUgzj18pZLCkXVUo5iEDylitbBRcCOZSzkL/lbMjL4ROQvIiBzAQyCb9EC0IdoAwx2rIDjoraHaLkBACUia/vDugXWE3OPNGMpj2oxwCgf6zIGzlZBDsiLbkVan1H7+yGQv5DqbgM1In+PSznWrDp/e7pGy8wuNbO9ZrbPzF674DnfZWZ3mtkdZvY7q16TdESvIgUyykMnBwBpkL/ppVwPXKZ3HPu6QPn6uETLzE6R9EZJl0k6X9JVZnb+zHPOk/TTkr7J3Z8j6TWrysoET4fo7E5Eh9ePSOs19VGM1BhkxEH+nihSTpQk0nolf8lfpBcpE0rHukZK7u2WFS6StM/dP+XuRyS9Q9IVM8/5YUlvdPeDo3L451a9aN17B+gNIQwAw+Nbn/rH+gXSqn1yE+lRB+vi2ugMnrPM7Oap5eqplz5b0j1Tv987/tu0r5b01Wb2/5nZR8zs0lXlpXaiUwwuhhFpHafu5FKfucFR5PTrIHUdSN0GJiLlAoaRuu6lbnupswd1I3OHxzofSZ291XFJbu0W6QvufuHUct2a775V0nmSLpZ0laTfMLMnLfsHMfZKO5C6k0/d0FLv5EmELuqVOn8A8ndYrO84yN/hpd7nBZBOrZnb0yVa90k6d+r3c8Z/m3avpBvd/R/dfb+kv9Zowmeh9LMCKAI7u8OLtM4jTDACNeKsSZC/dap1kBUFuZtOlHVP9lamn2/R+rik88xsl5mdKulKSTfOPOc9Gp29IzM7S6NLtj617EWpmdhYlKBFvTiaCGBo9H0j5C+GxKAaEvmLMrj7Y5JeLekDku6S9C53v8PMrjWzy8dP+4CkB8zsTkkflPRT7v7Astfd2mehAfRrx94jOrj71NTFAJJJPbhksAEAdWFyAajN6q88b8vdb5J008zfrpn62SX9xHhphD3TAqQcYNDJpcc2GEk90AeGRttPj20wQv4CwLDI3YH1c4lWL5jgAQrAIGOEzg7A0CLkL2eSYSgp+9kI9TxCe8dIhG0RoU5iAL7R16QPjloJoBN0cqgNZ09igu3BBDsA1KLKm7xzBo9kZuea2QfN7E4zu8PMfnzOc8zMXm9m+8zsNjN7YV/l6VOtRzPYoY2F7THCIKOu/AUQB/lL/paM/ax42CYYjrVchtfnTZYfk/ST7v4JMztT0l+Y2R+5+51Tz7lMo+9xP0/SN0j69fH/AbTATZcxRv4WjB3amMhfjPWev1UePUdrZ9x+X6PnPXLB2T2XpExn3r1FDz3rWOpioG+JzsZpo7cJHne/X9L9458fMrO7JJ0tabqDu0LSW8d3h/6ImT3JzJ42/rdYgbN3ujfpBHPu5BhkgPztH5ckdmt2AJJrBpO/IH/7Qeaup+mkzqJ/k1MG15y7h/dv17Zdh1IXow5M8JzIzJ4p6QWSPjrz0NmS7pn6/d7x307o4MzsaklXS9KZT3tiX8VExWY7wlw7udQ4ihEP+ds9Jte7s2gQUsqEz5DI33hKy9+aL8HLJXvbTOwsex2yF8hP73upZrZN0u9Jeo27f7HNa7j7de5+obtf+MQdp3VbwEwxwOjOqs7wjNvve3wBcjJE/nKpANpaJ1PJX+SG/V8Mqa/9VPZ/AY1vmGztlgR6nSUwsydo1Ln9trv//pyn3Cfp3Knfzxn/bS0pBxgpjmZwmmo32nRauXRyKSfhUtdPTlUdGSp/U6otf3OYXG+SkW0HDOTvauRvDDXkb02iZ+8Q2Rg9f6NvI+TPvd2SQp/fomWSrpd0l7v/6oKn3Sjp+8bfJvAiSYe4/ni51DtvOQRokzNuNumoondyAPnbj9T5G90kG6czeDaPN81P8hfRkb/dI3sXGzITyd/FqKMVyOhr0vu8B883SXqlpE+a2S3jv71O0k5Jcvc3SbpJ0ssk7ZP0sKQf6LE82SM8Vlt2L52u34frkhfjXhDJkb8dI3+Xa3rmTpfvRQbPR/4mR/5iECkmXMhfTOw5sFuX7NybuhjDSXS5VRt9fovWh7Xiy9/H3x7wo32VAd2KfvbO0B1d9Ememr9VoHbkb7ciTO5Ezt9UR3UjZzD5W6++85f7ng0rYvZGOJMmYv7Wmrt8k9YwLKNv0Uq/15qxIe//EGGAEVnKAQbmo86iTzV/m0skqTMw9ftHRf6ia2RuDJEyL1JZgF61vTyrtHvwDIUjGcOJOiueuoPhGwZQK/K3blFyL0o5AKBPEbMuYplSYWK9ZC2/QavEb9EC+hapY4lUlomIpxb3jaOMQP+i5V208kjkL9AFBs3xRcrfGnMXA+EMHnSNGyaeLFKHMhGxTADyF+kMyqg5x9mUALoUabKAbAPQFBM8yFLkjo5BRnocRQb6kUO25VDGkpG/QLdyyLQcyghshDN4UKoIR5Fz6URyKSfQRk333+EMyvyQvwAwLHIXRWOCZxg1DTCQpwidXaRTjIfEUeRysW3TiJBn64hwNiX5C+QpSttNnWE5irLthkTm9szFTZZrQENKI8eOLscyY31MOKNvKc+gzDnHci47miF/UaIcsyvHMneptpuC15S95u2WFOqqhUAitXd4ADZX41FJAOmlOqhZ22AZQGBcooVSpRxgMEkCAABqVdPR8trlvM+bc9mBEqyc4DGzD5nZl49//hEze42Zpb/TrqRLdu5NXQQgC6km5jj6tpnI+YthMbmeL/I3T+QvkC/OeEXXSrtEa7u7f9HMvk7SD0vaIek3+i0WIiIsgcGFzV8m2AEULmz+ol+p93dLmFQv4TMAudra4Dn/aGZbJX2fpF9293eZ2c09lws4AR0FKkX+Ah044/b79MgFZ6cuBvJC/gIARhJ9I1YbTSZ4Xi/pVkmnS3rt+G/beisRAGCC/AWXZwFpkL8YXEmZy8Q6ipHwhsltrLxEy93fKukbJF3g7o+Y2bMl/XnvJUMoqU9XLUHqTpttmB/yFygD+Zsf8hcA8LjSvkXL3Q+7+yPjn/e5+w/0WyzguNQTI9gMN/rcDPlbNyYGsAnydzPkL7CZlPvw3OAeXcrpJstNLtECAAAAAPSEA5pAYCVdogVwBBlAFIf3b09dhGow2ACQSqqzINjnRa7YP+pZRpdoZX0Gz54Du5O8Lw0IQO1S5S8AAAAwlJSXW7XBGTwIjSPIAAAAQJ7YlweGlfUZPAAAAAD6wVnrw2ASBAjOLXUJGmOCBwCAgFLeC4LBBgAAwFhGl2gxwYOluNkcgFncfwcAAAC17BPmdA8eJngAAFgi1be5AAAAIICMJnjYa81EjQMMLhEAYrpk597URQDWRp8CAGmkyl+uREAn/Pg3aa27pFDfrAGQEAMMAAAAAMiIt1wS4BItAAAAAEiAg39ABrhEC12q8fIsAIigxvwtebBR8mcDAACob8+1A9t2HUpdBGSMAcYwaKcAkAb5ixxxvxYAi+R0Dx4u0UJIfU6CbPWjuuaB90mSfunJl+l1D75fknTtU16ux+yU3t4XwOa27Tqkw/u3py4GWiJ/h3fm3Vv00LOOpS4GgMTIX6AOnMGD6lzzwPv03Ec/o+c++hm97bO/+fjPk04P3avxMhf0Y8jJHept98hfAEiD/AU2kNFNltl7bYGjx2U4XY9pmx/R6Xps8PfmMi0gPs566A/5CwBppMxfIEt8Tfow9hzYnboIyNQvPfkyPWYnVv3HbIt+8ckvS1QiIC/kL9oif+vAgTAgHvIX2ABn8KAUB3efmroInXvdg+/XVj/xyPxWP6afefCmRCVCHxhglIdt2r++z24hfwEgDfIX2EBGEzzcZHlNDDDK8SVt1WO25aTODt3jUhcA08hfAEiD/AXWY0p3uVUbnMGDcPo+gnztU16uT572dH3ytKfrlU/9gcd/vvYpL+/1fYEScHnWcEo8gzJK/tZ0Hx4m2AFIcfIXQL+yPIOHAQY28ZidomvOuuLx36d/HtIZt9+nRy44O8l7A0AKUfIXAGpD/gIb4AyeMnF5FrrE5A7QHPkLAOhTiWdNAuhAj9+iZWaXmtleM9tnZq9d8rzvMDM3swtXvWaWEzyX7NybugiD4vRqYH3bdh1KXQQAQTHB3i/yF2iOPAIy0MNNls3sFElvlHSZpPMlXWVm58953pmSflzSR5sUNcsJHiB3dObIVW0T7Cnt2HskdRGKRP4CwLBS5S5nZaEz/XyL1kWS9rn7p9z9iKR3SJp37eS/k/TLkr7UpKhM8GApBhjdY3ABAGmQvwAAYF0bXKJ1lpndPLVcPfWyZ0u6Z+r3e8d/O/6+Zi+UdK67/z9Ny5rtBA9HkZEjBhcoQU35yyWyAAAA81WzT9j+DJ4vuPuFU8t1Td/SzLZI+lVJP7lOUXub4DGzt5jZ58zs9gWPX2xmh8zslvFyTV9lQV6YBMGmar//A/kLnIy+ZRjkL/kLAEVpO7mz+hKt+ySdO/X7OeO/TZwp6QJJ/9XMPi3pRZJuXHWj5T7P4LlB0qUrnvNn7v788XLtum9QzYwhisDgAgO6QeQv8DjyFwO6QT3nLwBgWD19i9bHJZ1nZrvM7FRJV0q6cfKgux9y97Pc/Znu/kxJH5F0ubvfvOxFe5vgcfcPSXqwr9cHcsLgAkMif/PH/c+APJWWv7WfkTUU9hOB+rj7Y5JeLekDku6S9C53v8PMrjWzy9u+bup78Hyjmd1qZu83s+csepKZXT25MdHDBx8dsnwAUCryFydhkAEMgvxFNehXUIR+LtGSu9/k7l/t7s9y918c/+0ad79xznMvXnX2jpR2gucTkp7h7s+T9GuS3rPoie5+3eTGRE/ccdpgBYyEG33m7Yzb71v9JGA45C8ApEH+AkBmerpEqxfJJnjc/Yvufnj8802SnmBmZ6UqD2Jhtr97B3efmroIgzm8f3vqIoRG/q6HCXZsivzFBPnbTKrcramtoixcTtmzns7g6UOyCR4ze6qZ2fjni8ZleSBVeQCgFuQvlmGCHegP+YtFSszelJ+JyTp0pr9v0erF1r5e2MzeLuliSWeZ2b2SflbSEyTJ3d8k6RWS/qWZPSbpEUlXuvvaq+GSnXu158DuzsqNOB654GwubQJaGCp/0Z+Du0/lRstAhshfoG6c9VseGy+56G2Cx92vWvH4GyS9oa/3B4Bakb8AkAb5CwAFymgaPvW3aGENNc4Il3i6KlCKS3buTV0E9ITsBYDhlZS9JX0WgJssAwAAAAAAYDBM8KyhxruTc4OybnAUA0BuSsmtUj4HgGGk3vctIbNK+AzACTK6yTITPAiPTmJzKXdWary0sBZD3uA+9QQ79RhtpR4sAkBtyF10jgme4Qz9DVqpBxnIDxNUw6OdogTsoG6u9vxNMTFJ/vantvue1TyxnnN25Vx2YK6W99/hHjwAACzBwHF47KgDIHvTIH+BQDiDByXhKHJ7dM4AkAb5C6At9n3bqT13azvrrKazCjmDp3AcyRhe7R0GmqN9lq227csgAzmprX0CfWP/d3019ptk7wA4gwcAnTKAEuSYZTmWGQByRu6iZJzBg+JEmA2n42iHb9ACsCnytx3yF8CmcsjfHMoI1IIJHmQllw4kl3ICiC/CBHtOyF+gHCknKiNlb+Rci1w2oBNtL8/iDJ68cK1jOnQkeeDoMVAWsjcf5C+AGnHWJHrDBA9KxJGMZiKVLdI2K11N3ySA4UVpy49ccHaojJsVqWxRtlkNyF/UIFK+TUQsE9A1E/fgQY+YHT4u4kAjWnlSoZ6iT5xBmR5ZFxf5iz7UmrvRJmoj7ftGKQcwCM7gQamidXRSrM4ukojbCsgZA+cTRcveSGUhf9EXzlaClDbvomX/RK25W+vk69DMvdWSAhM8KEbqDidiZ5dCjYNgdrgxhKg7r6mzF8eRv0A9UuQuWY+JqrKXmywDadHhAcDwmGAH0IfUE5dRJ9elYbOPnF0sdR1F/3K6B8/WNG8L9GvSCZ1x+329v0dEkXdGAJRriOydfa9oyF8AQ3rkgrN7z9yoeTtB7gLHZX8GT8rTw1Jd88gscXNdd0iTSxGid3QpnXl3uljhOuRh1Zi/aKavjCSDlyN/gTr1lYvkLTCW0SVanMGDtR3cfap27D2SuhiNbXpEObeOjaMYQH8eetaxpIPonHR1Nk9OGUz+omTbdh3S4f3bUxcjiVz2fafzspb9XmAIqS63aoMJHlRj3VNY6eAAYHNtJ3rIYAATTK6vb93Jnlwzt+aJdc6cHBATPMO6ZOde7TmwO3UxkIHZzmteh5drByfV3ckBNcjlKPI8TSd6cs5gYAjs92Jd5Gp/uHVGBRLeMLmNIiZ4Uqr5dNUS0OEBm2GggXWVmrtMsANly3mCvVTkbhpVfUX6REYTPJzrmClmiwHUjlOTAQAA0CdTXl+TzgQPWmHGPB62CTCM1BPstHUAqaScWE+dvRL5G0mEbRGhTmIg7u2WBJjgQWsRghXx1NbZVXmaKoATROkPyV8AqANnMQ+LM3gqk6qB1bYjByxDR4eaRJlQACTyF3Uhf9NjGwCLMcGDjRCwAFKqeYKd/AXqVevZSxGyF0BlfIMlASZ4sDEGGelF2QbseAHDitL2kR75CwyH7E0nyrqvNXNrnVy2Y+2WFJjg6UjNR5GlOGELALUhf9Ni/aM2XJIHoDqcwYMasZMLIIXaJ9gBoCZRspf93uGxzkeYZB0eN1lOoNbTxYAonV2UHS4Mj/xNL0oO1CbKeid/AWA4ZG5lXHxNOoYVKWSi7OyiLhzJQCqR8hdIgfxFzdjvHQ7rGilxBg+AQUTp7BjkAulFyYNaRFnfteYvZw6mXQepJ/Yi1fsoWVAy1jGS4x48adTc2dPR1Yf1nF7NmQPUjPwFgDqlHnOlnFxlvzcPRU3wABhe6o4OMdTc6UdqA0w81CVS3QNqRvb2h3WL1ExcogWgZ3R2AJAG+QvEEG2Ck2wACtX2BsvcZBmbitTR0cn1h3V7otT3AQAk8rcWrNsTkb91Y/ujb9EyN1Jfj2FxBk9CNV8mEE20UEb36OgAIA3yFxL7vdGw79sd1uXJmFRNiJss1ylCo2OHr2x0doiMgUYcZEX3WKdAPBH3e8mKzUVchxHr2pBq38fjDB5gLGJA5yrauozQ0fFNAphggv1k0TIjZ9HWZYS6Rv4Ci0XLDAAbcEnHvN2SABM8BYqw4zeNTm5zrEMAbZAdAPrE5Dq6Rr81X4S2VjUu0ZLM7C1m9jkzu33B42ZmrzezfWZ2m5m9sKv3Tnlkh8Y3H2FdFnamYkuZv7WL2DbI381EW38R6xiOS5G/nNEUU7TsyEHUdUbuIid9nsFzg6RLlzx+maTzxsvVkn69x7JUJ2IQRQ3t6Fhv8zGZutQNSpS/TLCjJOTvfLS1pW4Q+78YI0OaY13FxSQy9+CRJLn7hyQ9uOQpV0h6q498RNKTzOxpfZUHyFHEzi7i5CFORP6mFbGNRMyS6CKus4h1CyeqMX8jTPhFbhsRsySayOsoct3CgNzbLSuY2aVmtnd8Vudr5zz+E2Z25/iMz//XzJ6x6jVT3oPnbEn3TP1+7/hvJzGzq83sZjO7+eGDjw5SuE1E6Oiiihzg0bCu4irgSEax+YvFyJTmWFdxkb/IEZmyGOtmNcaW6fVxBo+ZnSLpjRqd2Xm+pKvM7PyZp/2lpAvd/WslvVvSr6wqaxY3WXb369z9Qne/8Ik7Tmv0bwrYAdhY1Blngny1qOsoSp2ioxtOm/xNKUrdiNJWZkXNlkiirqModSpKG6tBbvmbWpQ2skjUbEkp+jqJXqcwkLY3WF59As9Fkva5+6fc/Yikd2h0lufxt3b/oLs/PP71I5LOWfWiKSd47pN07tTv54z/hgpED/SUWDexFTJ53Gv+FrKOikXGLMa6ia2QbOklf1OvGyb+miFjjmNd5CF1tkRgksy91bJC4zM6x35Q0vtXvWjKCZ4bJX3f+NsEXiTpkLvfn7A8nYrS0UWeeT64+1TCfQbrY7UobStz5O8AoucvTsQ6WS1K28pc0fmbWuTcnSBr8lgHUeoSuRvEsZaLdNbkUtvxcnWbtzez75V0oaR/v+q5W9u8QcNCvF3SxRp9qHsl/aykJ0iSu79J0k2SXiZpn6SHJf1AX2VBbAd3n6ode4+kLkZy0Tu7KB0dViN/0QTZmw/yN5+jyOQvmqg5f6Pv7wLzNDgbZ5EvuPuFCx5rdEanmb1U0s9I+u/cfeUN2Xqb4HH3q1Y87pJ+tK/3R15q7ujQXOqjGLkMMMhfNDXZ0a49fxlwrJY6f3ORMn8v2blXew7s7uOlG9m265AO79+e7P0nHnrWMZ15d/zbjNa475tL1jKpjoF8XNJ5ZrZLo4mdKyV9z/QTzOwFkv6LpEvd/XNNXjR++m0g9WAsys5QLiGVS+j3oebPjjKRvyPkb3w1f/ZcpM4ToC813a6gls/ZpdT7MmTvWE83WXb3xyS9WtIHJN0l6V3ufoeZXWtml4+f9u8lbZP0u2Z2i5nduKq4vZ3BA7RR29HkXDq7CINUOjmgX7UdTSZ/m0udv8C6cjmLZ6Lk/M0layciZC6ican9JVrLX9n9Jo0u3Z3+2zVTP7903dfMJ/kyFWWnKLewyq0zaKOGzwggL7UcTa7hMwIpRNnvzVGJuVTiZ0KdzNstKRQ/wcNR93yVOtDI7XPlNjkITEQZaOTYhnLKqHWQvygd+73H5dh+csuoRXL9HJHqTJR9GIy5t1sSKH6CJ4IoDTRSaK0jxw5ikZI+y5BStyF2mNthvR2XY/7muoO+SEmfZUjkLzC8nPM313LjRGTvFJfsWLslBSZ4kIWcOzop3/LnOCgFpqUenJYgx+yaRv4Cw4qUu7m3o5yyK9esnYhUVyK1IYxxBg9mRWmokcKrjZw6jklHl1OZp+VeV4Bocm5TOWZZjmWeyLmuANHk3p6iZ1n08gGd6OFbtPpSxQQPp5iVJfrESeSyAUOLkL9RJthLED3fovcPwBDI3TJFy7Zo5dlE7pOAXYuQIdGYe6slBb4mvUK5fXXkMlG+Vr2UDm6Cjg7oRyn5GyV7J0rKYPIX6F4p2SudmHdDZ3BJWTsRLXOZHMWmmOAZ0LZdh3R4//bUxZBUVkcnpevs6Oj6R0eHLkTK35Iw0OgW+QtgHUNkcIlZC6wt0dk4bVQzwXPJzr3ac2B36mJgAH13dnR0ADZR2gT7xFCTPWRwPbhMoL0I+73RJtZLzd6Jedm4bhbXlq9Mqp+M3J3DJcWqKktVM8GDk5Xe0UndDjhq6PSidXQR0NFtLsJAA8PqerKH/AXQhRr2fafVkJ1A30zp7qfTBhM8A4t2NKMmdHLLRRxcRDiSAfShpkEG2bsa+YtSsd+LqCLmLgLLaIKnjr1LLES4AahFtAEr+QugT5yBOh/Zi4h1INo+Cma4t1sSYIIHQMiOLgJ2jgH0jfydj/xFn2h3wMnI3QUm9+BpsyRQ1QQPlXY+Orm6Rd3+HMlADaK2Pwwj6vYnf9El6hMiiZi7tJH4zL3VkkJVEzxRRGzEEcMOQDmiTLBHzF8A6EuU7I2Ifd/6sM1RAyZ4gIrR0S3GTjGGQjusE9t9MfK3PFEn12mH9Yi6raO0DXJ3Be7BE1eUyhulMU+LGnzoR+TtHbF9oBwR61fk9oi6RGwfQJ/IXwDLtZzcYYIHEdDJAehLlAn2qMjfekTd1kzulClK9lK/kAqZu1yUjAjLxQRPdFEqcZRGPStqCKI7kbdxhHYRJSPQnwj1bJ7IbRPdYBsvR/4iBdpludi26ATfooWmGGRgaJG3bdT2gO4weFstchvFZiJvW/IXQ4hczyK3T7QTeZtGbgs4Gd+ilQEGGatFDkUA2FTknSvyF0CX2O9thuwtB9uyGbKhIS7RwjoYZGAokbdn5HaAbkXamYhc7yK3V6wv8vaM0g4iZQP6E6W+LRK5raIM0dsAZrikY95uSaDqCZ5IOxKRGzodXRnYjs1EygWAdlsGtiMioH9rjjabN7Yfalb1BA+aIyjzFn37RZ7gRD8iDTSof+gT+QucLId6F73tYr7o2y1S3Y+0LxYbX5OOQkUPTMzHdmuOjg4R0YbzxbZrjvxFRLThvLC90BsmeNBGpBldlCGHjo56jwii18Mc2jJOlMM2i17v0b1IE2nUP3SJzF1PpCzIAhM8+aByryeH8MRIDtsqUkeH4ZG/68mhTWMkh20VKX/JAkSWQ3uuHdsIveImy9hEpB2uRQjR+NhG62OAAfIXXWAbIbpI/V0OuSvRriPLZdvkUtcxj0t+rN2SABM8aCWXMK1RLtuGjg5oJ5c2XqNctg35i0hyqY+5tO+asE3aiTTJmw0u0cIm6OjQFtukHTo6TOSSv4iH/G2H/EVOaOdx5LQt2LfIHJdo5Yedi/ZyCtfS5bQt6OgwQf6289CzjmXV5kuX07YgfyHFy96c6iX5m15O6z9a3Y7W9tE9JniCihYGy+QUsqViG7RHR4dZ5C/WwTZoj/xFzmj7w2NyDclwiVZ+2MnYDGGbTm7rPqfBM4ZB/m4mtwwoSW7rnvzFtGjZm2P9zC0Dcpbjus6xTmMBJnjQhdxCIcfgzVmORzGi1eloO7eIyJRWpAAADaBJREFUI1pdXSW3LMgd+bs58hfzRKunTeSWBTliHXeD3G2r5eQOEzzpRaz0uXV0BPAwWM8oTcT8zQ25MAzWM0pC9naDXOhPrus2tzEclnBJx461WxJgggedyzWIc5Hr+o3W0bFTi1Wi1dkmcs2HXOS6fqPVZfIXy0Srr03leGZfZDmvz4h1mNzdEGfw5Cti5Y8YEqvkHMqRsU6BYeWav+ge6xVAE2TF5liHCIcJHmCEgO5G7hNm0QbJESdywXbpUu6ZEUnu65L8xSoRt0m0eruu3HMjpdzXW8S6G7GN58WlYy2XBJjgmSNiI4gYFk3RyW0m93WXc90FpLzrcO75kVru6y/nuguUUH9zz5AhlTBeKKHOIn9M8GQk99DIPbSHRkfXj4gTuDgu6vaJWJebKiFLhlbCOotYZ6O2b7Bt+lRCnvSJ9dMv2nYHXHI/1mpJodcJHjO71Mz2mtk+M3vtnMdfZWafN7NbxssP9VmedURtDBF32NZBiDdTwjrKva7mLuf8RT/I32ZKWEfkb1rkb3dKqsslZEvXSlonJdVVzMElWpKZnSLpjZIuk3S+pKvM7Pw5T32nuz9/vLy5r/IgFgYa87Fe0IXc85cJ9n6RM/OxXtCFnPOX7O0fOTNS2nqIWkejtukscZNlSdJFkva5+6fc/Yikd0i6osf361zURhE1RNooLeDbKm09RK2jUdt0D7LP36ii1u02SsudtkpbD1HrKPmLTUSt122VljtN1fq5kTl36dixdksCfU7wnC3pnqnf7x3/bdZ3mNltZvZuMzt33guZ2dVmdrOZ3fzwwUf7KGt2Su3oagv9Ej9zaXUzU9nnb+TBYGl1vMQcaqLEz11a3cxU1vkbOXtLVGIOzVPy54yau7TljnEGT2PvlfRMd/9aSX8k6bfmPcndr3P3C939wifuOG3QAkYWNVA2VXInMFHDZ4yGju4k5O8GSszfWibaa/iM0ZC/JyF/WygxdydKzaVSP9dEyXUSJ/Jjx1otKfQ5wXOfpOkjEueM//Y4d3/A3SeHJN4s6et6LE8r0XdKSg6WEjuFEj/TtJLrY2bIX2ykxKwq8TNNI3/DyD5/I2dv6fW8hIn2Ej5DE5HrYuQ2nKeWZ+8kOoNna4+v/XFJ55nZLo06tislfc/0E8zsae5+//jXyyXd1WN5kKlJB3Hm3alPOGuv9E5OoqMLhvwdwLZdh3R4//bUxegV+ZsH8jcU8hedmM6u6BlcQ86iYq5k34jVRm9p4e6PSXq1pA9o1HG9y93vMLNrzezy8dN+zMzuMLNbJf2YpFf1VZ5NRN85ibxj16XcjgbUcgRDqqcO5oL8HU4tdT+3LCN/kUop+Rs5e2us8xEzLWKZhhK5DkZuuxhGn2fwyN1vknTTzN+umfr5pyX9dJ9lqEUNR5InIh/RqLGTi67Wjo78RR/IX6yD/D3hb+Rvh2ra7501m3VDZTEZOxJ5cgc98nzqf68TPCW5ZOde7TmwO3UxMCNVJ7fo/WtER4e+Rc/fWgcaqSd7yF/yF/0ie/OwLAvbZDPZulj0zK11Ur1vLskzukSLCZ6C0NEt7pS6GnzQ6Z2Ijg4YqT1/+87eZe9RK/IXIHtXITe7Ez1z0SP33s7gMbNLJf2fkk6R9GZ3/z9mHj9N0ls1uhn/A5K+290/vew1meBZQ/QjGRId3SJ0cN2jo8OQcshfnIzs7Qf5i6HkkL3s+wJMqvetjzN4zOwUSW+U9K2S7pX0cTO70d3vnHraD0o66O7PNrMrJf2ypO9e9rqxLqAHkIUcBhd0dBhaDu0C+cuhnpG/AEqSQ+6iZ36s3bLcRZL2ufun3P2IpHdIumLmOVdI+q3xz++W9BIzs2Uvap7o+9nbMrOHJEXbczhL0hdSF2IGZWqGMjVDmZrZ7e5npi5EX8jfxihTM5SpGcrUDPk7vIj1gDI1Q5maoUzNlJ6/f6jRem/jdElfmvr9One/bvy6r5B0qbv/0Pj3V0r6Bnd/9dR73z5+zr3j3+8eP2dhHcjxEq297n5h6kJMM7ObKdNqlKkZytRM1DKlLkPPyN8GKFMzlKkZytQM+Tu8qPWAMq1GmZqhTM2Unr/ufmnqMqyDS7QAAAAAAACGc5+kc6d+P2f8t7nPMbOtkrZrdLPlhZjgAQAAAAAAGM7HJZ1nZrvM7FRJV0q6ceY5N0r6/vHPr5D0J77iHjs5XqJ1XeoCzEGZmqFMzVCmZijT8CJ+PsrUDGVqhjI1Q5mGF/HzUaZmKFMzlKkZylQId3/MzF4t6QMafU36W9z9DjO7VtLN7n6jpOslvc3M9kl6UKNJoKWyu8kyAAAAAAAATsQlWgAAAAAAAJljggcAAAAAACBzYSd4zOzTZvZJM7tl3lev2cjrzWyfmd1mZi8MUKaLzezQ+PFbzOyaAcr0JDN7t5n9lZndZWbfOPN4ivW0qkyDricz2z31XreY2RfN7DUzzxl0PTUsU4r69K/N7A4zu93M3m5mp888fpqZvXO8nj5qZs8MUKZXmdnnp9bTDw1Qph8fl+eO2e02fnzwdtcl8rdxmcjf1eUhf5uXi/xtVibyl/wlf5uVh/xtXi7yt1mZis7fYrh7yEXSpyWdteTxl0l6vyST9CJJHw1QposlvW/g9fRbkn5o/POpkp4UYD2tKtPg62nqvU+R9FlJz0i9nhqUadD1JOlsSfslnTH+/V2SXjXznP9F0pvGP18p6Z0ByvQqSW8YcD1dIOl2SU/U6Eb1fyzp2VHqU0efkfxtVibyd72ykb+Ly0H+NisT+Uv+Tt6T/F2vbOTv4nKQv83KVHz+lrKEPYOngSskvdVHPiLpSWb2tNSFGpKZbZf0LRrdXVvufsTd/37maYOup4ZlSuklku5297+d+XvK+rSoTClslXSGmW3VKMA/M/P4FRrtwEjSuyW9xMwscZmG9jUadVgPu/tjkv5U0v8885zS86n0z7cS+dsK+bsc+bsa+Vv+51uJ/G2F/F2O/F2N/M1E5Akel7THzP7CzK6e8/jZku6Z+v3e8d9SlkmSvtHMbjWz95vZc3ouzy5Jn5f0m2b2l2b2ZjP7spnnDL2empRJGnY9TbtS0tvn/D1FfZpYVCZpwPXk7vdJ+g+SDki6X9Ihd98z87TH19M43A9JekriMknSd4xPBX23mZ3bV3nGbpf0z8zsKWb2RI2OVsy+Z8r61AXydzXyd33k7wLkb2PkL/krkb9tkL8LkL+N1ZC/RYg8wfPN7v5CSZdJ+lEz+5bUBdLqMn1Co9MMnyfp1yS9p+fybJX0Qkm/7u4vkPQPkl7b83uu0qRMQ68nSZKZnSrpckm/O8T7NbGiTIOuJzPbodHM+y5JT5f0ZWb2vX2+Z0dleq+kZ7r710r6Ix0/wtILd79L0i9L2iPpDyXdIulon++ZAPm7Gvm7BvJ3ZVnI3wbI32TI39XI3zWQv52UifzFXGEneMYzl3L3z0n6A0kXzTzlPp04a3jO+G/JyuTuX3T3w+Ofb5L0BDM7q8ci3SvpXnf/6Pj3d2vUuUwbej2tLFOC9TRxmaRPuPvfzXls8Pq0qkwJ1tNLJe1398+7+z9K+n1JL555zuPraXzK6HZJD6Qsk7s/4O6Pjn99s6Sv67E8k/e83t2/zt2/RdJBSX8985RU9akT5G8j5O96yN/lyN+GyF/yV+Tvusjf5cjfhkrP31KEnOAxsy8zszMnP0u6RKPTwqbdKOn7bORFGp26dn/KMpnZUyfXY5rZRRqt394av7t/VtI9ZrZ7/KeXSLpz5mmDrqcmZRp6PU25SotPBR10PTUpU4L1dEDSi8zsieP3fYmku2aec6Ok7x///ApJf+I+uqtaqjLZidf2Xj77eB/M7CvH/9+p0fXHvzPzlFT1aWPkbzPk79rI3+XI34bIX/KX/F0b+bsc+dtQyflbFA9wp+fZRdJXSbp1vNwh6WfGf/8RST8y/tkkvVHS3ZI+KenCAGV69fixWyV9RNKLB1hXz5d0s6TbNDqFcUfK9dSwTCnW05dp1Dlsn/pb6vW0qkwp1tPPS/orjXbe3ibpNEnXSrp8/PjpGp1Ou0/SxyR9VYAy/e9T6+mDkv7pAGX6M4123G6V9JII9anDz0b+Ni8X+dusTORvszKRv83KRP6Sv+Rv8zKRv83KRP42K1Ox+VvSYuONAQAAAAAAgEyFvEQLAAAAAAAAzTHBAwAAAAAAkDkmeAAAAAAAADLHBA8AAAAAAEDmmOABAAAAAADIHBM8AAAAAAAAmWOCBwAAAAAAIHNM8KA6Zva7ZvYGM/uwmf2tmX2zmb3NzP7azK5PXT4AKBX5CwBpkL9AHbamLgCQwHMl/bm7v9rMXifpekkXS/q8pHvN7DR3fzRlAQGgUOQvAKRB/gIVYIIHVTGz0yU9SdJ/Hv/JJV3v7vePHz8q6Uii4gFAschfAEiD/AXqwSVaqM1zJH3C3Y+Nf3+epI9KkpmdI+kz7u6pCgcABSN/ASAN8heoBBM8qM1zJd069fvXSrpt/PPzpn4GAHSL/AWANMhfoBJM8KA2z5V0i/T46apnuPvB8WPTnR0AoFvkLwCkQf4ClTDOxgMAAAAAAMgbZ/AAAAAAAABkjgkeAAAAAACAzDHBAwAAAAAAkDkmeAAAAAAAADLHBA8AAAAAAEDmmOABAAAAAADIHBM8AAAAAAAAmfv/Ad7NvrSxRrYwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x360 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M, S = np.meshgrid(m_range, s_range,indexing='ij')\n",
    "plt.rcParams['figure.figsize'] = (16,5)\n",
    "fig, axs = plt.subplots(ncols=4, gridspec_kw=dict(width_ratios=[4,4,4,0.2]))\n",
    "axs[0].contourf(M, S, W)\n",
    "axs[0].scatter([m_0],[s_0],marker='X',c='red')\n",
    "axs[0].set_xlabel('$m$')\n",
    "axs[0].set_ylabel('$s$')\n",
    "axs[0].set_title('POT')\n",
    "\n",
    "axs[1].contourf(M, S, W_est)\n",
    "axs[1].scatter([m_0],[s_0],marker='X',c='red')\n",
    "axs[1].set_xlabel('$m$')\n",
    "axs[1].set_ylabel('$s$')\n",
    "axs[1].set_title('Estimated truth')\n",
    "\n",
    "axs[2].contourf(M, S, W_0)\n",
    "axs[2].scatter([m_0],[s_0],marker='X',c='red')\n",
    "axs[2].set_xlabel('$m$')\n",
    "axs[2].set_ylabel('$s$')\n",
    "axs[2].set_title('True')\n",
    "fig.colorbar(axs[np.argmax([W.max(), W_est.max(),W_0.max()])].collections[0], cax=axs[3])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}