{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6946d418-e2b3-4500-b6d7-d1ffcd0bcd95",
   "metadata": {},
   "source": [
    "# 1-D prior toy example (RBF covariance)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10d5c70-b5dc-4de6-a6d2-4d10e3b1d833",
   "metadata": {},
   "source": [
    "> Demonstration of prior error bounds for a 1-D toy example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80eebe9b-ec85-49a8-b895-e7724eddecc0",
   "metadata": {},
   "source": [
    "In this notebook we will demonstrate the error bounds for the statFEM prior for the toy example introduced in <a href=\"/statFEM/oneDim.html#1-dimensional-case-(ODE)\"><code>oneDim</code></a>. We first import some required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0bf09fed-be73-42d8-b397-fe5d23f0b67b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "set_log_level(LogLevel.ERROR)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (10,6)\n",
    "# import required functions from oneDim\n",
    "from statFEM_analysis.oneDim import mean_assembler, kernMat, cov_assembler\n",
    "from scipy.stats import linregress\n",
    "from scipy import integrate\n",
    "from scipy.linalg import sqrtm\n",
    "from tqdm 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": "d6cac7ef-3ef0-4dd2-b835-268b1ce4bcea",
   "metadata": {},
   "source": [
    "The random forcing term has the following GP distribution:\n",
    "\n",
    "$$f\\sim\\mathcal{GP}(\\bar{f},k_{f})$$\n",
    "\n",
    "where,\n",
    "\n",
    "$$\\bar{f}(x)=1$$\n",
    "\n",
    "$$k_{f}(x,y)=\\sigma_{f}^{2}\\exp\\left(-\\frac{|x-y|^{2}}{2l_{f}^{2}}\\right)$$\n",
    "\n",
    "$$\\sigma_{f}=0.1$$\n",
    "\n",
    "$$l_{f}=0.4$$\n",
    "\n",
    "This forcing gives rise to the true solution mean and covariance functions:\n",
    "\n",
    "$$\\mu_{\\star}(x)=\\frac{1}{2}x(1-x)$$\n",
    "\n",
    "$$c_{u}(x,y)=\\int_{0}^{1}\\int_{0}^{1}G(x,w)k_{f}(w,t)G(t,y)dtdw$$\n",
    "\n",
    "where $G(x,y)$ is the Green's function for our problem:\n",
    "\n",
    "$$G(x,y)=x(1-y)\\Theta(y-x)+(1-x)y\\Theta(x-y) \\quad \\forall x,y\\in[0,1]$$\n",
    "\n",
    "(note: $\\Theta(x)$ is the Heaviside Step function). We now set up all these functions, noting that we use quadrature to compute $c_u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8ce7fd7c-4452-485a-b543-bae276c3e81a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up true mean and f_bar\n",
    "μ_star = Expression('0.5*x[0]*(1-x[0])',degree=2)\n",
    "f_bar = Constant(1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4a7520dc-9375-4412-9649-ebdcc94e9b93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up kernel functions for forcing f\n",
    "l_f = 0.4\n",
    "σ_f = 0.1\n",
    "\n",
    "def c_f(x,y):\n",
    "    return (σ_f**2)*np.exp(-(x-y)**2/(2*(l_f**2)))\n",
    "\n",
    "# translation invariant form of c_f\n",
    "def k_f(x):\n",
    "    return (σ_f**2)*np.exp(-(x**2)/(2*(l_f**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d7878225-c01a-4f36-a5ea-9ab35df5eb7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up true cov function for solution\n",
    "# compute inner integral over t\n",
    "def η(w,y):\n",
    "    I_1 = integrate.quad(lambda t: t*c_f(w,t),0.0,y)[0]\n",
    "    I_2 = integrate.quad(lambda t: (1-t)*c_f(w,t),y,1.0)[0]\n",
    "    return (1-y)*I_1 + y*I_2\n",
    "\n",
    "# use this function eta and compute the outer integral over w\n",
    "def c_u(x,y):\n",
    "    I_1 = integrate.quad(lambda w: (1-w)*η(w,y),x,1.0)[0]\n",
    "    I_2 = integrate.quad(lambda w: w*η(w,y),0.0,x)[0]\n",
    "    return x*I_1 + (1-x)*I_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2c69805-4662-4429-ae6c-8e904e27d4e6",
   "metadata": {},
   "source": [
    "We now set up a reference grid on which we will compare the true and statFEM covariance functions. We take a grid of length $N=51$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "956ccb91-e41a-49e4-aabe-f9cb807cf6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up reference grid\n",
    "N = 51\n",
    "grid = np.linspace(0,1,N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f064ffbc-10bb-47b2-9e33-901f102c469a",
   "metadata": {},
   "source": [
    "We now compute the true prior covariance on this reference grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6074fbde-0a03-4aa6-b7f1-e6e04bdccf27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the true cov matrix and its square root\n",
    "C_true = kernMat(c_u,grid,True,False)\n",
    "tol = 1e-15\n",
    "C_true_sqrt = np.real(sqrtm(C_true))\n",
    "rel_error = np.linalg.norm(C_true_sqrt @ C_true_sqrt - C_true)/np.linalg.norm(C_true)\n",
    "assert rel_error <= tol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ae6dcb-3fe0-4856-9929-ce519e23b1b7",
   "metadata": {},
   "source": [
    "We now set up a function to get the statFEM prior, for a FE mesh size $h$, using functions from <a href=\"/statFEM/oneDim.html#\"><code>oneDim</code></a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1b53cd72-ae20-4e5a-b568-f965ce7e737a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up function to compute fem_prior\n",
    "def fem_prior(h,f_bar,k_f,grid):\n",
    "    J = int(np.round(1/h))\n",
    "    μ = mean_assembler(h,f_bar)\n",
    "    Σ = cov_assembler(J,k_f,grid,False,True)\n",
    "    return μ,Σ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b13d0567-190d-4b5e-a64d-dc9578d47943",
   "metadata": {},
   "source": [
    "We now set up a function to compare the covariance functions on the reference grid. This function computes an approximation of the covariance contribution to the 2-Wasserstein distance discussed in <a href=\"/statFEM/oneDim.html#Difference-between-true-prior-covariance-and-statFEM-prior-covariance\"><code>oneDim</code></a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "00c476b6-ec3b-4caa-9185-f884a3dd01f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to compute cov error needed for the approximate wasserstein\n",
    "def compute_cov_diff(C1,C_true,C_true_sqrt,tol=1e-10):\n",
    "    N = C_true.shape[0]\n",
    "    C12 = C_true_sqrt @ C1 @ C_true_sqrt\n",
    "    C12_sqrt = np.real(sqrtm(C12))\n",
    "    rel_error = np.linalg.norm(C12_sqrt @ C12_sqrt - C12)/np.linalg.norm(C12)\n",
    "    assert rel_error < tol\n",
    "    h = 1/(N-1)\n",
    "    return h*(np.trace(C_true) + np.trace(C1) - 2*np.trace(C12_sqrt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "331c810b-5779-4503-9032-4810e4c356e5",
   "metadata": {},
   "source": [
    "With all of this in place we can now set up a function which computes an approximation of the 2-Wasserstein distance between the true and statFEM priors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c99edfd1-9bed-40c5-b2ae-97bbe2fa4987",
   "metadata": {},
   "outputs": [],
   "source": [
    "def W(μ_fem,μ_true,Σ_fem,Σ_true,Σ_true_sqrt):\n",
    "    mean_error = errornorm(μ_true,μ_fem,'L2')\n",
    "    cov_error = compute_cov_diff(Σ_fem,Σ_true,Σ_true_sqrt)\n",
    "    cov_error = np.sqrt(np.abs(cov_error))\n",
    "    error = mean_error + cov_error\n",
    "    # error = np.sqrt(mean_error**2 + cov_error**2)\n",
    "    return error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "430432e8-fb90-4bb3-9643-d1da2c98d2a9",
   "metadata": {},
   "source": [
    "We set up a range of $h$ values on which to compute this error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7a73360d-5abd-4a64-ba25-cab1c770027e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.25, 0.2 , 0.17, 0.14, 0.12, 0.11, 0.1 , 0.09, 0.08, 0.08, 0.07,\n",
       "       0.07, 0.06, 0.06, 0.06, 0.05, 0.05, 0.05, 0.05, 0.04, 0.04, 0.04,\n",
       "       0.04, 0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set up range of h values to use\n",
    "h_range_tmp = np.linspace(0.25,0.02,100)\n",
    "h_range = 1/np.unique(np.round(1/h_range_tmp))\n",
    "np.round(h_range,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73a37e69-b0a6-4ee6-ac90-3ad679328736",
   "metadata": {},
   "source": [
    "We now compute the errors for this range of $h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "77abe78e-40df-4ff8-969c-a6d66b17ae8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "set_log_level(LogLevel.ERROR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5b92f86d-15c4-46c9-aa6c-c89df817c56e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 30/30 [00:00<00:00, 46.23it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 657 ms, sys: 3.44 ms, total: 660 ms\n",
      "Wall time: 653 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "errors = []\n",
    "for h in tqdm(h_range):\n",
    "    # obtain the fem prior for this value of h\n",
    "    μ, Σ = fem_prior(h,f_bar,k_f,grid)\n",
    "    # compute the error between the true and fem prior\n",
    "    error = W(μ,μ_star,Σ,C_true,C_true_sqrt)\n",
    "    # append to errors\n",
    "    errors.append(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "659a14fc-c982-4d61-b758-46ca276a3959",
   "metadata": {},
   "source": [
    "We now analyse the results by plotting the errors on a log-log scale and using ordinary least squares regression to estimate the slope of the line of best fit. This gives us an estimate of the rate of convergence as if we assume:\n",
    "\n",
    "$$W(\\nu_{\\star},\\nu_{h})\\leq\\mathcal{O}(h^{p})\\text{ as }h\\rightarrow 0$$\n",
    "\n",
    "then we have \n",
    "\n",
    "$$W(\\nu_{\\star},\\nu_{h})=Ch^{p}+\\mathcal{O}(h^{p+1})$$ \n",
    "\n",
    "from which it follows that:\n",
    "\n",
    "$$\\log W(\\nu_{\\star},\\nu_{h})=p\\log h + \\log|C| + \\mathcal{O}(h)$$\n",
    "\n",
    "and so the slope is the rate of convergence $p$. Our theory gives us $p=2$ and so we should expect to obtain a slope of $2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "cd8ac4ab-fb39-49cc-8086-6ce980eedd5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "errors = np.array(errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "740d82e2-04a4-4665-8181-bb662a3e52f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope: 2.000523    intercept: -2.184977    r value: 0.999881     p value: 0.000000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAGJCAYAAAA66h/OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3gU1R7G8e8hhB66IL2LIC0QmlwFpIkgIqJSRBALzYpiVLyK6FWSKE00CIKIhiIgTTpCLChiIPSidAhIUQIJJKSd+8eGSCeQhNlN3s/z5CE7Ozv72wzB1zNzfsdYaxERERER95PN6QJERERE5PIU1ERERETclIKaiIiIiJtSUBMRERFxUwpqIiIiIm5KQU1ERETETSmoiYiIiLgpBTURERERN6WgJiI3zBiz1xjT8ia8T1VjzHpjTJQx5vmMrsMY84Ex5sUbfO0aY8wd6VXLVd4n3T7z5T6vMeaAMabuZfa9KZ9PRFwU1ETkqpIDQYwxJtoYc8QYM8kYk+8GjpGWUPEqsNJa62OtHZ2G41yTMeYW4HHgs+TH640xD533fG1jjDXG9Dxv2x3JIbIg8CEwNCNrTE8Xf97kbYWAksC2y7zEoz6fiKdTUBOR1LjfWpsPqAv4AW/e5PcvB2y5Se/VC1horY1JfnwSyH/e8y9cZttzwGRrbSQwD2hujLn1JtSaHnpx4ecFqAnsvWjbOZ72+UQ8moKaiKSatTYCWATUuPg5Y0w1Y0yoMSbSGLPFGNMheftXQFlgfvKo3KuXO/ZVXr8CaA6MSX79bdeq80rHSn6urjEmPHkEbIYxZrox5r3zXt4W+OG8xymhLHn0qQPwBVAgeVtBoDvwcfLPKBZYC7S5Qm3+xpiZF20bZYwZfd7zEcn17TDGtLjJnxegFrDTGDPaGHPMGHPIGNMqNZ9PRNKXgpqIpJoxpgxwHxB+0XZvYD6wFCiGa4QpxBhT1VrbA9hP8qictTbwMse92uvvAX4Cnk1+/R/XqPGKxzLG5ABmA5OAwsBU4MGLDlET2HHe4/NHz/oA05I/z7ltvYFV1trt571mG1D7CiVOA+4zxvgk1+sFPAJMMcZUBZ4F6ltrfXCFob03+fOe2+YHLACK47os6p/Kzyci6UhBTURSY44xJhL4Gdfoy/sXPd8IyAcMs9bGWWtXAN8BXVN5/LS+PrXHagRkB0Zba+Ottd8Cay56fUEg6rzHJ4H8yYGoLzAaOAUUMMZkAwYAoy46RlTycS5hrd0HrOPfwHQPcMZauxpIBHIC1Y0x3tbavdbaXTf584JrRO1/1tol1tokYGtqP5+IpC8FNRFJjY7W2oLW2nLW2v6XuXepJHAg+T/q5+wDSqXy+Gl9fWqPVRKIsNba8547cNHrTwA+5z0+N6LWGdiQPKJ3Ctelz3ZAArD4omP4AJFXqXEK/4bQbsmPsdbuBF4EhgBHjTHTjDElr3IcSOfPa4wxuC5tzz9vnxpcGNau9flEJJ0oqIlIejgElEkeYTqnLBCR/L299CXX9fr0quUwUCo5jJxT5qLXbwTOvw/uXFB7gX9Hzk4lb3sO+PiiIARQDdhwlRpnAM2MMaVxjaxNOfeEtXaKtfY/uCZQWCDgKseB9P+8FZL/3HneNl9g/XmPr/X5RCSdKKiJSHr4DTgDvGqM8TbGNAPux3U/FsARoGIaXp9etfyK6/Lis8aY7MaYB4AGF71+IdD0vMcngbsBH2vt0uRtp4A6QH1c93+lMMbkAuoBy65UoLX2GBCKa1LCHmvttuTXVjXG3GOMyQnEAjFA0pWOk0Gftxaw6aLw6UtyMEvN5xOR9KOgJiJpZq2NwxUO2gLHgU+Bx8+7wf4D4M3kWYmv3MDr06WW5Oc6AU/iunT3GK77uc6ed4jJuG72z538+CSuS4jn9287BZQAJllroy8q4X4g1Fp76BqlTgFact5oGq7704Yl1/0XrskBr9/kz1uT80bLjDFFgVuBzdf5+UQkHZhLR+xFRLIOY8xvwFhr7RfnbXsfOGqtHXmDx3vSWrv5mjs7IK2f190/n0hmo6AmIlmKMaYprnYUx3H1PxsLVLTWHna0sAyS1T6vSGaT3ekCRERusqrAN0BeYDfQOZOHlqz2eUUyFY2oiYiIiLgpTSYQERERcVMKaiIiIiJuKlPeo1a0aFFbvnz5q+5z+vRp8ubNe3MKkmvS+XAvOh/uRefDveh8uJfMcD7Wrl173Fp7y+Wey5RBrXz58oSFhV11n9DQUJo1a3ZzCpJr0vlwLzof7kXnw73ofLiXzHA+jDH7rvScLn2KiIiIuCkFNRERERE3paAmIiIi4qYy5T1qlxMfH8/BgweJjY0FoECBAmzbts3hquScAgUKsGfPHkqXLo23t7fT5YiIiLiFLBPUDh48iI+PD+XLl8cYQ1RUFD4+Pk6XJclOnTpFXFwcBw8epEKFCk6XIyIi4hayzKXP2NhYihQpgjHG6VLkMowxFClSJGXEU0RERLJQUAMU0tyczo+IiMiFslRQc9ro0aOpVq0a3bt3Z968eQwbNgyAOXPmsHXr1pT9Jk2axKFDh67r2Hv37qVGjRpX3R4WFsbzzz+fhk8gIiIiN1OWuUfNHXz66acsX76c0qVLA9ChQwfAFdTat29P9erVAVdQq1GjBiVLlkzX9/fz88PPzy9djykiIiIZRyNqN0nfvn3ZvXs3bdu2ZcSIEUyaNIlnn32WX375hXnz5jFo0CDq1KlDQEAAYWFhdO/enTp16hATE8PatWtp2rQp9erVo02bNhw+fBiAtWvXUrt2bWrXrs0nn3xyzRpCQ0Np3749AEOGDKF37940a9aMihUrMnr06JT9vv76axo0aECdOnXo06cPiYmJGfNDERERkavKmiNqL75I7rVrwcsr/Y5Zpw6MHHnFp8eOHcvixYtZuXIlRYsWZdKkSQDceeeddOjQgfbt29O5c2cAFi1axIcffoifnx/x8fE899xzzJ07l1tuuYXp06czePBgJk6cyBNPPMGYMWO4++67GTRo0HWXvH37dlauXElUVBRVq1alX79+7Ny5k+nTp7Nq1Sq8vb3p378/ISEhPP744zf0YxEREZEblzWDmgfZsWMHmzdvplWrVgAkJiZSokQJIiMjiYyM5O677wagR48eLFq06LqO3a5dO3LmzEnOnDkpVqwYR44c4fvvv2ft2rXUr18fgJiYGIoVK5a+H0pERMSNWWs5efIkBQsWdLqULBrURo4kxkP6qFlrueOOO/j1118v2B4ZGZnmY+fMmTPley8vLxISErDW0rNnTz744IM0H19ERMQTLV68mKFDh/LLL7843pFA96i5AR8fH6Kioi77uGrVqhw7diwlqMXHx7NlyxYKFixIwYIF+fnnnwEICQlJl1patGjBzJkzOXr0KAD//PMP+/btS5dji4iIuKO4uDgmTZrEF198AUCrVq14+umnsdY6XJmCmlvo0qULQUFB+Pr6smvXLnr16kXfvn2pU6cOiYmJzJw5E39/f2rXrk2dOnX45ZdfAPjiiy8YMGAAderUSbe/TNWrV+e9996jdevW1KpVi1atWqVMXhAREclMoqOjGTFiBJUqVeKJJ55g2rRpAGTPnp3evXuTLZvzMcm4Q1pMb35+fjYsLOyCbdu2baNatWopj7WElHs5dz4uPk/ijNDQUJo1a+Z0GZJM58O96Hy4lxs9HyEhITz33HOcOHGCpk2b8tprr9GmTRtHLnUaY9Zaay/bPytr3qMmIiIiWc6+ffvImTMnt956KyVLlqRp06b4+/vTqFEjp0u7IufH9EREREQy0KZNm+jRoweVKlVKmSzXvHlzZs+e7dYhDRTUREREJJNatWoV7du3p1atWsyePZvnn3+eV155xemyrosufYqIiEimYa1Nuc9s/Pjx/PbbbwwdOpQBAwZQuHBhh6u7fh4zomaMedkYY40xRZ2uRURERNxLfHw8kydPpmbNmqxbtw6AwMBA9u3bx3//+1+PDGngIUHNGFMGaA3sd7oWERERcR8xMTGMGjWKSpUq0bNnT7Jly0Z0dDQAxYoVI0+ePA5XmDYeEdSAEcCrQObrJXIZP/30E3fccUfKouzuqFevXsycOROAp556iq1btzpckYiIZDWJiYk8+eSTvPjii5QvX54FCxawYcOGlOUVMwO3v0fNGPMAEGGt3XC13ibGmGeAZwCKFy9OaGjoBc8XKFDggu7/iYmJFzy+may1WGuv2Ehv0qRJvPjii3Tp0oWEhIRr1nmt42WE+Ph4YmJiiIqKYsSIEQBp+nmeOx+xsbGXnDu5+aKjo3Ue3IjOh3vR+XDWX3/9xfLly+nevTvGGLp27UqFChWoUaMGAD/88IPDFaazc/+Rd/ILWA5svszXA8BvQIHk/fYCRa91vHr16tmLbd269YLHp06dumSfjLRnzx5722232R49etjq1avbvXv32iVLlthGjRpZX19f27lzZxsVFWXHjx9vCxUqZMuXL2+7detmrbU2MDDQ+vn52Zo1a9q33nrruo5nrbXlypWzb731lvX19bU1atSw27Zts9ZaGxUVZXv16mVr1Khha9asaWfOnGmttVc8zvl69uxpZ8yYYa21tmnTpvb333+31lqbN29e+8Ybb9hatWrZhg0b2r/++staa+3Ro0dtp06drJ+fn/Xz87M///zzBcc7dz4uPk/ijJUrVzpdgpxH58O96Hw4Y/Pmzfbxxx+32bNnt9mzZ7cbNmyw1maO8wGE2StkGre49GmtbWmtrXHxF7AbqABsMMbsBUoD64wxt6b1Pe+77z6aNWt2wdenn34KwJkzZy55rlmzZkyaNAmA48ePX/Jcavz555/079+fLVu2kDdvXt577z2WL1/OunXr8PPzY/jw4Tz11FN06NCBoKAgQkJCWLp0KX/++Sdr1qxh/fr1rF27lh9//DHVxzunaNGirFu3jn79+vHhhx8C8O6771KgQAE2bdrExo0bueeeezh+/PhVj3Mtp0+fplGjRilDz+PHjwfghRde4KWXXuL3339n1qxZPPXUU6k+poiIZF1HjhyhQ4cO1KhRg5kzZ/Lss8+ya9cuatWq5XRpN4VbX/q01m4Cip17nBzW/Ky1xx0rKg3KlSuX0lhv9erVbN26lSZNmgCuBWEbN258yWuWLl3K0qVL8fX1BVxD7n/++Sdly5a9ruN16tQJgHr16vHtt98CsHz58pR1zQAKFSrEd999l6q6riRHjhy0b98+5b2WLVuW8l7n38d26tQpoqOjyZcvX6qPLSIiWUNSUhL79++nfPnyFCpUiAMHDvDOO+8wYMAAihQp4nR5N5VbB7WMtHDhwiuu9ZknT56r3n9QtGjRG7o/IW/evCnfW2tp1aoVU6dOveprrLW8/vrr9OnT54Lte/fuva7j5cyZEwAvLy8SEhKu+n6pqetKvL29U/rXnP9eSUlJrF69mly5ct3QcUVEJPOLj49n2rRpBAYGcuLECXbv3k2OHDlYt26dI2twugO3uPSZWtba8p46mnaxRo0asWrVKnbu3Am4Lhn+8ccfl+zXpk0bJk6cmDLVOCIigqNHj97w8c7XqlUrPvnkk5THJ06cuKHjpEbr1q35+OOPUx6vX78+zccUEZHM4fTp04wePZrKlSvz+OOPY63lgw8+SJkkl1VDGnhYUMtMbrnlFiZNmkTXrl2pVasWjRs3Zvv27Zfs17p1a7p160bjxo2pWbMmnTt3vuzsytQe73xvvvkmJ06coEaNGtSuXZuVK1fe0HFSY/To0YSFhVGrVi2qV6/O2LFj03xMERHJHH7++WdeeOEFypQpw/z589m4cSM9evQge/Yse+EvhXFNNshc/Pz8bFhY2AXbtm3bRrVq1VIeR0VFXfHSp9x8587HxedJnBEaGprqSTKS8XQ+3IvOR9odOHCA4cOHU6BAAYYMGYK1lrCwMOrXr3/dx8oM58MYs9Za63e55zSiJiIiIjfF1q1b6dWrFxUrVmTMmDEcO3YMcF3avJGQlhUoqImIiEiG++ijj7jjjjuYMWMG/fv3Z+fOnRfcJy2Xp4u/IiIiku6stSxatIhKlSpRtWpVWrZsydtvv82zzz5L0aJFnS7PY2hETURERNJNQkICISEh1K5dm3bt2jFmzBgAateuzZAhQxTSrpOCmoiIiKSLCRMmUKVKFR577DESExP58ssvr2t1G7mULn2KiIjIDYuMjKRAgQIYY9i8eTMlSpRg9OjRtGvXLqUPmtw4/QQdNGTIkJR1Ny9nzpw5Fyy7JCIi4i4OHjzIyy+/TJkyZVixYgUAAQEBrFq1ivvvv18hLZ3op+jGFNRERMTdbNu2jd69e1OxYkVGjRrFAw88QOnSpQHXes9ZeRWBjKCgdgVzwiNoMmwFFV5bQJNhK5gTHpEux/3f//7Hbbfdxn/+8x927NgBwPjx46lfvz61a9fmoYce4syZM/zyyy/MmzePQYMGUadOHXbt2nXZ/URERG6WhIQEWrRowbRp0+jTpw87d+7k66+/pmrVqk6XlmkpqF3GnPAIXv92ExGRMVggIjKG17/dlOawtnbtWqZNm8b69etZuHAhv//+OwCdOnXi999/Z8OGDVSrVo0JEyZw55130qFDB4KCgli/fj2VKlW67H4iIiIZxVrL4sWL6datG/Hx8WTPnp3p06ezb98+Pv74Y8qXL+90iZmegtplBC3ZQUx84gXbYuITCVqyI03H/emnn3jwwQfJkycP+fPnp0OHDgBs3ryZu+66i5o1axISEsKWLVsu+/rU7iciIpIWCQkJTJ06FV9fX9q2bcuPP/7I7t27Abjrrru45ZZbHK4w61BQu4xDkTHXtT2tevXqxZgxY9i0aRNvv/02sbGxadpPRETkRu3fv5/bbruNbt26cfbsWb744gt2796ty5sOUVC7jJIFc1/X9tS6++67mTNnDjExMURFRTF//nzAtSB5iRIliI+PJyQkJGV/Hx8foqKiUh5faT8REZG0OHHiBD/88AMApUuXpkmTJsyZM4ctW7bQq1cvcuTI4XCFWZeC2mUMalOV3N5eF2zL7e3FoDZp+7+JunXr8uijj1K7dm3atm2bsgDtu+++S8OGDWnSpAm33357yv5dunQhKCgIX19fdu3adcX9REREbkRERASvvPIKZcuW5cEHHyQmJoZs2bLx1Vdf8cADD6jFhhtQw9vL6OhbCnDdq3YoMoaSBXMzqE3VlO1pMXjwYAYPHnzJ9n79+l2yrUmTJhe05+jXr99l9xMREbkee/bs4X//+x+TJ08mMTGRLl268Oqrr5I7d9quHEn6U1C7go6+pdIlmImIiLiL+Ph4vL29OXr0KCEhITz99NO88sorVKhQwenS5AoU1ERERDIxay1Lly4lICCAihUr8vnnn9OwYUMOHTpEoUKFnC5PrkEXn0VERDKhhIQEpk+fTt26dbn33nvZsWMHtWvXTnleIc0zZKmgZq11ugS5Cp0fEZH0M2TIELp06UJMTAwTJ05kz549PPfcc06XJdcpy1z6zJUrF3///TdFihTROmRuyFrL33//Ta5cuZwuRUTEI0VGRhIcHEyzZs1o3LgxTz31FPXq1dPsTQ+XZYJa6dKlOXjwIMeOHQMgNjZWocCNxMbGUrBgwZSFfUVEJHUOHTrEyJEjGTt2LFFRUbz99ts0btyY8uXLa4mnTCDLBDVvb+8LZrWEhobi6+vrYEVyPp0PEZHr9/rrrzN8+HASEhJ45JFH8Pf3p06dOk6XJelIY6EiIiIeJDw8nMRE13rUhQoVonfv3vzxxx9MnTpVIS0TUlATERFxc9Zali1bRsuWLalbty6zZ88G4NVXXyU4OJhKlSo5XKFkFAU1ERERN5WUlMSMGTPw8/OjdevWbN26laCgIFq3bu10aXKTZJl71ERERDyFtRZjDNZaBg8ejDGG8ePH06NHD3LmzOl0eXITKaiJiIi4iZMnTzJ27FhCQkL49ddfyZs3L8uWLaN06dJ4eXk5XZ44QJc+RUREHHb48GFee+01ypYty2uvvcatt97K8ePHAShXrpxCWhamETUREREH7d69m+rVqxMfH0/nzp3x9/enbt26TpclbkJBTURE5CZbu3YtGzZsoHfv3lSsWJF33nmHhx56iMqVKztdmrgZXfoUERG5Cay1fP/997Rq1Qo/Pz/eeOMNzp49C4C/v79CmlyWgpqIiEgG+/3332nQoAEtW7Zk8+bNBAQEsGPHDs3glGvSpU8REZEMcPbsWSIjIylevDj58+fn1KlTfPbZZzz++ONaa1pSTUFNREQkHZ0LZCNGjKBx48bMmjWLqlWrsn37dowxTpcnHkZBTUREJB0cOXKEUaNG8emnn3Ly5ElatGhBv379Up5XSJMb4fZBzRgzBHgaOJa86Q1r7ULnKhIREbnUmDFjGDZsWEqLjXr16jldktygOeERfLRoK3+dOEOxIj4MalOVjr6lHKnFUyYTjLDW1kn+UkgTERHHhYeH06VLF7777jsAXnzxRXbs2ME333yjkObB5v22m7X/DWLSR0/Qbf0iIiJjeP3bTcwJj3CkHk8JaiIiIo6z1hIeHk6bNm2oW7cuCxcu5NChQwAUKVKEKlWqOFyh3LBTpyAoiDtbNeDdBaM4nSM3+wqWACAmPpGgJTscKcvtL30me9YY8zgQBrxsrT3hdEEiIpL1PProo8yYMYPixYvzwQcf0LdvXwoWLOh0WZIWR47AqFHw6adw8iTby9Xh0/te4pdyteG8+woPRcY4Up6x1jryxhcUYcxy4NbLPDUYWA0cByzwLlDCWtv7Msd4BngGoHjx4vWmTZt21feMjo4mX758aaxc0ovOh3vR+XAvOh/OiYuLY8WKFdxzzz3kyJGD77//nn/++YcOHTqoB5qbuNHfj1wREZT55htKLFqESUjg2N13c6BrV8IKlCQuMemS/XN4ZaPqrT7pUfIlmjdvvtZa63e559wiqKWWMaY88J21tsbV9vPz87NhYWFXPVZoaCjNmjVLt9okbXQ+3IvOh3vR+bj5oqKiUlpsHDp0iKlTp9KlSxdA58PdXPf5CA+HgACYMQOyZ4eePeGVV+C22wDXRILXv91ETHxiyktye3vxQaeaGTahwBhzxaDm9pc+jTElrLWHkx8+CGx2sh4REcm84uLiGDp0KJ988gmRkZE0b96cSZMm0bJlS6dLk7SwFkJDYdgwWLoUfHxc4ezFF6FEiQt2PRfGgpbs4FBkDCUL5nZ01qfbBzUg0BhTB9elz71AH2fLERGRzObUqVPkz58fb29vli1bxj333IO/vz8NGjRwujRJi6QkmDvXFdDWrIFixeD996FfP7jKvYUdfUs5Fswu5vZBzVrbw+kaREQkc9qwYQMBAQEsWLCA3bt3U6RIEX788Ufdf+bp4uLg668hKAi2b4eKFSE42HWZM3dup6u7Lm4f1ERERNKTtZYffviBgIAAFi9eTL58+ejbt2/K8wppHiwqCsaNgxEjICIC6tSBqVOhc2fX/WgeyDOrFhERuUE7d+6kefPmFCtWjP/973/069ePQoUKOV2WpMXRo/DxxzBmDERGQvPmMGECtG59QYsNT6SgJiIimVpcXBwhISH8+eefvP/++1SpUoX58+fTokULcnvYZTC5yJ49VBk1ChYvhrNn4cEHwd8fMtG9hVqZQEREMqWoqCiGDx9OxYoV6d27N0uXLiUuLg6A9u3bK6R5so0boXt3qFKFEt995/p+2zaYNStThTTQiJqIiGRCS5YsoUuXLkRGRtKsWTMmTJhA69atMR5+GSxLsxZ++sk1g3PRIsiXD156idUNGnDnww87XV2G0YiaiIhkCnv37mXzZlerzZo1a9KyZUtWr17NypUradOmjUKapzrXYuPOO6FpUwgLg/feg/37ISiIuFtucbrCDKWgJiIiHm3Tpk089thjVK5cmZdeegmAkiVLMmPGDBo2bOhwdXLD4uLgyy+hRg3o2BH++gs++QT27YPBgyGLTABRUBMREY+0evVq2rVrR61atZgzZw4vvPACX3zxhdNlSVpFR8PIkVCpEvTqBd7eMGUK/Pkn9O/vcX3Q0kr3qImIiMdISkrCWouXlxerVq1izZo1vPvuu/Tv35/ChQs7XZ6kxfHjrhYbH38MJ064LnOOGwf33uvxLTbSQiNqIiLi9uLi4vjyyy+pWbMmX331FQD9+/dn3759vPnmmwppnmzfPnj+eShbFoYOdQW0X391rc3Ztm2WDmmgoCYiIm4sOjqakSNHUrlyZXr16oWXlxfFihUDIHfu3OTJk8fhCuWGbdoEPXq4LnEGB8Ojj8KWLTB7NjRq5HR1bkOXPkVExG3df//9hIaGcvfddzN27Fjatm2r2ZseZk54BEFLdnAoMoaSBXPzQbGT3P3tBFiwAPLmdY2mvfQSlCnjdKluSUFNRETcxr59+xg9ejRvvfUWBQoUYMiQIeTIkYPGjRs7XZrcgDnhEbz+7SZi4+JpviuMfqtnUj9iK2cLFibnO+/As8+CLltflYKaiIg4bvPmzQQGBjJlyhSMMdxzzz20a9eOpk2bOl2apMHwhVu4N3wZfX+bSdXj+zmYvxhvtezDz3d1YMVb9zldnkdQUBMREcfExsby8MMP891335E3b16ee+45Bg4cSBldBvNsp0/DhAlMDXyPUqeOsb1oOV5o/zILbr+LBK/smDPW6Qo9hoKaiIjcVElJSWzcuJE6deqQK1cufHx8GDp0KAMGDNDsTU/3998wZoyrxcbff3OsfE0Gt+5PaEW/C2ZvliyYtXqhpYWCmoiI3BTx8fFMmzaNgIAAduzYwa5duyhbtixTpkxxujRJq/37YfhwGD8ezpyBDh3A35+9ucvx27ebID4xZdfc3l4MalPVwWI9i4KaiIhkqNOnTzNhwgQ++ugj9u/fT40aNZg0aRIlSpRwujRJq61bITAQQkJcj7t1g0GDXMs+AR2Tdzt/1uegNlXp6FvKmXo9kIKaiIhkCGstxhiOHz/OwIEDufPOOwkODlaLjczgl18gIADmzYM8eWDAABg40NW09iIdfUspmKWBgpqIiKSr/fv3M3z4cA4cOMCsWbMoV64c27dvp3Llyk6XJmlhLSxc6ApoP/3kaqsxZIirxUaRIk5Xl2lpZe2bVaIAACAASURBVAIREUkXW7ZsoWfPnlSqVIlPPvkEHx8f4uPjARTSPFlCguvSZu3a0L497N0Lo0a57kt7+22FtAymETUREUmzb775hkcffZQ8efIwYMAABg4cSNnLXAYTD3LmDEycCB995Apnd9wBX34JXbuCt7fT1WUZCmoiInLdrLUsXLiQXLly0aJFC1q3bs3QoUPp378/RTTC4tn++Qc++QRGj4bjx+HOO13ft2sH2XQh7mbTT1xERFItISGBkJAQateuTfv27RkxYgQABQsW5L///a9Cmic7ePDfCQFvvQUNG7ruRVu1Cu6/XyHNIfqpi4hIqkyfPp0qVarw2GOPkZiYyJdffsns2bOdLkvSats2eOIJqFjRNXLWsSNs2ADffQf/+Y/T1WV5uvQpIiJX9M8//5ArVy7y5MnDmTNnKFmyJKNHj6Zdu3Zk0wiLZ1u92jWDc84cyJ0b+vZ1jaiVL+90ZXIe/ZaJiMglDh48mDIhYNy4cQD07NmTVatWcf/99yukeSprYdEiaNYMGjeGH36A//4X9u1zjaYppLkdjaiJiEiKbdu2ERgYSEhICElJSXTt2pVWrVoBKJx5soQEmDHDNYK2YQOULu1a8unppyFfPqerk6tQUBMRkRT9+/fnt99+o2/fvgwcOJDyGmHxOHPCI1KWbCqfNxsjz4RTe/rnsGcP3H47fPGFa6mnHDmcLlVSQUFNRCSLstayePFihg8fzuTJkylRogTBwcEULVqUokWLOl2e3IA54RG8/u0mvKNO0i98Ib3D5lL0zEn+qVmXwrOHuxZL18ioR1FQExHJYhISEvjmm28ICAhg48aNlC5dml27dlGiRAluv/12p8uTNPhixipe/H463dYvwicuhpUV6xHc6GEiavixqmMLp8uTG6CgJiKShcTExFCrVi127txJtWrVmDRpEl27diWHLoN5th07ICiIbyZ9SfakJL67/S7GNnqIbcUqAmBOxjpcoNwoBTURkUzuxIkTLFmyhC5dupA7d2569OhBnTp1aN++vSYIeLrff4dhw2D2bMiZk+/q38fI2h04UPDWC3YrWTC3QwVKWimoiYhkUgcPHmTEiBGMGzeO06dP07hxY8qVK8dbb73ldGmSFtbCsmWugLZyJRQsCG+8Ac8/j1dEPMe/3QTxiSm75/b2YlCbqg4WLGmhoCYikskcPnyYN998k6+++oqkpCS6dOnCq6++Srly5ZwuTdIiIQFmzoTAQAgPh5Il4cMP4ZlnwMcHgI7FXLuem/VZsmBuBrWpSkffUg4WLmmhoCYikklERUXh4+ODt7c3c+fOpU+fPrz88stqseHpYmNh0iQICoLdu6FqVZgwAbp3h5w5L9m9o28pBbNMREFNRMSDWWtZunQpAQEBnD59mtWrV1O0aFEOHjxIrly5nC5P0uLkSQgOhpEj4cgRaNDAFdY6dlSLjSxEZ1pExAMlJCQwbdo06taty7333ssff/zBI488QmKi694khTQPdvgw+PtD2bLw+utQpw6sWOFam7NTJ4W0LMYjRtSMMc8BA4BEYIG19lWHSxIRcdTXX3/NE088we23387EiRPp3r27Wmx4uj//dI2Yffml6360zp3htdfA19fpysRBbh/UjDHNgQeA2tbas8aYYk7XJCJys0VGRvLpp59SpkwZevToQZcuXShUqJAWSM8M1q51zeCcNcu1rFPv3vDKK1CpktOViRvwhN/ufsAwa+1ZAGvtUYfrERG5aQ4dOsTYsWMpW7YsgwcP5qeffgJclzYfeOABhTRPZS0sXw4tW4KfHyxd6ho927vXdV+aQpok84Tf8NuAu4wxvxljfjDG1He6IBGRm+HDDz+kQoUKzJgxg/bt2xMeHs64ceOcLkvSIjERZsyA+vWhVSvYsgUCAmD/fnj/fbj11msfQ7IUY611ugaMMcuBy/3tHAz8D1gJPA/UB6YDFe1FhRtjngGeAShevHi9adOmXfU9o6OjyZcvX9qLl3Sh8+FedD6cs337dkqUKEGBAgX49ddf+e2332jfvj2VK1d2ujRJlprfj8iYeI6cjCUuMYkcXtkokScbVX9cSZnp08lz8CBnSpfmwKOPcqR1a5J0b2GaZIZ/r5o3b77WWut3uefcIqhdjTFmMRBgrV2Z/HgX0Mhae+xKr/Hz87NhYWFXPW5oaCjNmjVLz1IlDXQ+3IvOx81lrWX58uUMGzaMFStWMGTIEN5+++2U53U+3Mu1zsec8Ahe/3YTMfGJ5Dt7hu7rF/JU2FxuiT4B9eq5ZnJ27AheXjev6EwsM/x+GGOuGNTcfjIBMAdoDqw0xtwG5ACOO1uSiEj6mDlzJh988AHr1q2jZMmSBAUF8cwzzzhdlqRB0JId5I08zrNh8+gRvpD8Z0/zY3lf3n30DUaPfxmMcbpE8SCeENQmAhONMZuBOKDnxZc9RUQ8SUJCAtmzu/75/frrr4mOjubzzz/nscceI+dlOs2LB9m1i37TP+ThTcvxTkxgYdUmjG3Umc23VsYAoxXS5Dq5fVCz1sYBjzldh4hIWp08eZLg4GA+/vhjQkNDqVKlChMnTqRAgQJ46TKYZwsPd00KmDGDR4wXM2u0YFyDB9lb+N+lnEoWzO1ggeKp3D6oiYh4usOHDzNy5EjGjh3LqVOnaN26NWfPngWgcOHCDlcnN8xaWLnSFdCWLnUtjP7KK3zf4lHe/ekoMfGJKbvm9vZiUJuqDhYrnkpBTUQkA505c4Zq1aoRFRVF586d8ff3p27duk6XJWmRmAjffutqUvv771C8OHzwAfTtCwUL0hY4e0sEQUt2cCgyhpIFczOoTVUtlC43REFNRCSdrV27lrlz5zJ06FDy5MlDcHAwDRo0oJKamHq2s2fh669p8M47cOCAqynt2LHQsydctLZqR99SCmaSLjyh4a2IiNs712KjVatW+Pn5MWrUKA4ePAhA165dFdI8WVQUfPghVKwITz1FYq5cMH067NgBffpcEtJE0pOCmohIGv355580aNCAVq1asWXLFgIDAzlw4AClS5d2ujRJi6NH4c03oWxZGDQIqlWDpUtZ+9ln8Mgj6oMmN4UufYqI3ICzZ8+yd+9eqlatSsmSJfH29mbcuHE8/vjjarHh6Xbvho8+gokTXZc7O3UCf3/Xsk8AoaGOlidZi4KaiMh1OHXqFGPHjmXkyJH4+PiwdetW8ubNyy+//OJ0aZJWGza4ZnBOn+4aLevZE155BapqtqY4R0FNRCQV/vrrL0aNGkVwcDAnT56kVatW+Pv7ky2b7iDxaNbCjz+6ZnAuXgz58sHAgfDii1BKkwHEeQpqIiJXYa3FGMOPP/5IQEBASouNevXqOV2aXKc54f+2zCiVPycf5dpHw2/Gw2+/wS23wHvvQf/+UKiQ06WKpFBQExG5jHXr1hEQEICvry+vvfYaDz30EH/88QeVK1d2ujS5AecWSk+IjaXzllD6/DaLyv8c5HSpsuT99FPo1Qtya+UAcT8KaiIiyay1rFy5kmHDhrFs2TLy58+Pr68vAF5eXgppHmzMvPV0+2UOT62ZTYnov9larALP3T+IDQ1b8WO/Vk6XJ3JFCmoiIskGDhzIyJEjKV68OMOGDaNv374UKFDA6bIkLY4dg48/ZmbQCArGRvNr2Zr4t32eHyvUBWMwUXFOVyhyVQpqIpJlnT17lq+++oqWLVtSvnx5unTpwu23307Pnj3JpSamnm3fPleLjc8/h9hYNlRvwoi6D7K+5IUzOLVQurg7TVcSkSzn1KlTfPjhh1SoUIGnn36aKVOmANCwYUP69OmjkObJNm2CHj3+Xd6pa1fYupUTX01jR7nqF+yqhdLFE2hETUSylHfeeYeRI0cSGRlJixYt+PLLL2nZsqXTZUlaWAs//+zqgbZgAeTN62qv8eKLkLw6RMfkXbVQungaBTURyfQOHTpEyZIlAThw4AAtWrTA39+f+uc6zYtnSkqC775zBbRffoGiRWHoUBgwAAoXvmR3LZQunkhBTUQyrfXr1xMQEMA333zD6tWrqV+/PuPGjVOTWk8XFwdTp0JgIGzdCuXLw8cfQ+/ekCeP09WJpCsFNRHJVKy1/PDDDwwbNowlS5bg4+PDyy+/TJkyZQAU0jxZdLRrcsDw4XDgANSsCSEhrgXSs+s/Z5I56W+2iGQqp0+fpmPHjuTMmZP333+ffv36UbBgQafLkrQ4ftw1YjZmDPzzD9x9t2uiQNu2YIzT1YlkKAU1EfFocXFxhISEMH/+fGbOnEm+fPlYunQpNWvWJLc6zXuc85d58rWnCDqwnEpzp8GZM/DAA+DvD40bO12myE2joCYiHikqKopx48YxYsQIIiIiqFOnDkePHuXWW2+lQYMGTpcnN+DcMk9lDu3iw99m0WHbjwDsb9eJssOGQPXqVz+ASCakoCYiHmfDhg00a9aMyMhImjdvzoQJE2jdujVGl8E82uLPZjH6+xBa7VzDGe+cTK7bngn1H8CULccqhTTJohTURMQj7Nmzh127dtGyZUuqV6/Oww8/zJNPPknDhg2dLk3SwlpYuBCGDWPszz9zIpcPI5p048t67YnMnR8AExnjcJEizrnuoGaMyQvEWmsTM6AeEZELbNy4kYCAAKZPn07p0qXZvXs33t7ejBs3zunSJC3i42H6dFcPtM2boWxZRrQfwLgqzYnJceHKEFrmSbKya85TN8ZkM8Z0M8YsMMYcBbYDh40xW40xQcaYyhlfpohkNeHh4dx3333Url2befPm8cILL7Bq1Sq11/B0Z864ZnBWruxa6slamDwZdu6kwtDXXasKnEfLPElWl5oRtZXAcuB1YLO1NgnAGFMYaA4EGGNmW2u/zrgyRSQrSEpKIjY2ljx58nD8+HHCwsJ477336N+/P4UKFXK6PEmLv/+GTz5xhbTjx6FJE1e7jXbtIDl8n1s1QMs8ifwrNUGtvbU26uKN1tp/gFnALGOMd7pXJiJZRlxcHFOmTCEwMJC2bdvy0Ucf0bJlS/bt26cWG57uwAFXg9rx4+H0aWjf3tVi4z//uezuWuZJ5EKpuYawzRjzqTGm1pV2sNbGp2NNIpJFREdHM2LECCpVqsQTTzyBt7c3d955JwDGGIU0T7Z1K/TqBRUrukbROnWCTZtg/vwrhjQRuVRqRtSqAt2Az40x8cBnwHRr7dkMrUxEMr0XX3yRCRMm0LRpU8aPH0+bNm3UYsODnN+cNuUy5dkDMGwYzJ0LuXND//4wcCCUK+d0uSIe6Zojatba09ba8dbaBsAAoDGuUbaPjDFVMrxCEck09u3bx/PPP8+GDRsA8Pf359dffyU0NJR7771XIc2DnGtOGxEZg7WWKmt/pFTHe12rBvz4I7z1FuzfD6NGKaSJpME1R9SSZ3XmB3yS//wZ2AU8BbwIeGVkgSLi+TZt2kRgYCBTp04lW7Zs1KxZk9q1a1OlShWqVNH/73mioCU7iDsbR4ftP9Fv9UyqHdvLIZ+ijGrXnxemBUC+fE6XKJIppObS5x9ABDAbOAFEA1HAu8l/iohcUdeuXZk2bRp58+bl+eefZ+DAgZQuXdrpsiQtzpyhxYoZPL1mNmVOHuHPImV4+b6XmFf9bhK8vHlBIU0k3aQmqNUF+uBqxTENmGytPZKhVYmIx0pKSmLFihW0aNECYww1atTg3XffpX///hQuXNjp8iQtTpxwtdgYNYqhx4+ztuTtDG3xNMsrN8Aa1500pdScViRdXTOoWWvXA/2SVyToBswxxhwAxlprV2R0gSLiGeLj45k6dSoBAQFs3bqV5cuX06JFCwYPHux0aZJWEREwYgR89hlER8N99/FTpyd5Zk9uYhKSUnZTc1qR9Hc9Lb6TgLlAD1wNcD81xmzLkKpExGOcPXuW0aNHU7lyZXr27ImXlxdff/01TZs2dbo0Savt2+HJJ6FCBRg5Ejp0gA0bYMEC7nqyEx88VItSBXNjcI2kfdCppnqgiaSz1EwmOAFY4DRwKvkrCtic/L2IZEGJiYl4ebnmEgUEBFCpUiWCg4Np27atZm96ujVrXC025syBnDmhTx94+WUoX/6C3dScViTjpeYetcLWWpvhlYiIR9i/fz/Dhw9n6dKlbNiwgZw5c7Ju3TqKFy/udGmSFtbC0qWuRdJXroSCBWHwYHjuOShWzOnqRLKs1AS1azLGGIU5kcxty5YtBAYGMmXKFAC6d+9OVFQUhQsXVkjzEJdtUFuzOMyc6Qpo69dDqVLw0Ufw9NPg4+N0ySJZXqoWZTfGzALmWmv3n9tojMkB/AfoiWvh9kkZUaAxZjqu1REACgKR1to6GfFeInJ5a9asoWHDhuTJk4dnn32Wl156ibJlyzpdllyHcw1qY+ITATh+LJL1g4fRatN88h7cB1WrwsSJ0L075MjhcLUick5qgtq9QG9gqjGmAhAJ5MLV6HYpMNJaG55RBVprHz33vTHmI+BkRr2XiLhYa1mwYAF//fUXTz75JH5+fowePZpu3bpRpEgRp8uTGxC0ZAcx8Ynkj43msfCFPBE2j1vORLKl9O3c8e238MADkO165peJyM2QmvYcscCnuGZ5egNFgRhrbWRGF3c+47o7+RHgnpv5viJZSXx8PNOnT+ett95iz5491KpViyeeeIJs2bLx3HPPOV2epEHCgYO8FjaX7usX4RMXQ2iFeoxt9BC/lanJngfbO12eiFyB8ZRby4wxdwPDrbV+V3j+GeAZgOLFi9ebNm3aVY8XHR1NPnXPdhs6H85bs2YNw4cP58iRI5QtW5bu3btzzz33kD17utzKKmmQlt+P3AcOUGb6dIovWYJJTOLPRk1Y1/5BjperCEAOr2xUvVX3ol0P/XvlXjLD+WjevPnaK+WbVP8LbIwZeJnNJ4G1yU1xb5gxZjlw62WeGmytnZv8fVdg6pWOYa0dB4wD8PPzs82aNbvqe4aGhnKtfeTm0flwxt9//83Zs2cpWbIk+fLlY/78+Xz++efkyZOHe+7R4LW7uKHfj7Aw1wSBWbMgZ052P9idPsWa8me+Yq7GSptcDWo/6FSTZmqxcV3075V7yezn43puSPAD+gKlkr/64Lp/bbwx5tW0FGGtbWmtrXGZr7kAxpjsQCdgelreR0RcDhw4kDIh4NzKAX5+fvz000+0b9+ebLpXyTNZC8uWQcuWUL8+LF8Ob7wB+/ZR8ZtJDHiytRrUiniY67mmURqoa62NBjDGvA0sAO4G1gKB6V9eipbAdmvtwQx8D5FMb+vWrQQGBhISEoK1lm7duvHyyy87XZakVWKia+QsIADWrYMSJSAw0NWoNn/+lN3UoFbE81xPUCsGnD3vcTxQ3FobY4w5e4XXpJcuXOWyp4ikTnBwMN988w39+/dn4MCBlCtXzumSJJUu2wOtWhGYPBmCgmDnTrjtNvj8c3jsMdeKAiLi8a4nqIUAvxlj5gIGaA9MSV6sfWtGFHeOtbZXRh5fJDOy1rJo0SICAgIYOnQoTZs25b///S9vv/02RYsWdbo8uQ6RMfG8/v2/PdBOHTnOrle+JnbDd+T6+xj4+bma1nbsCMnLeolI5pDqoGatfdcYswhokrypr7U2LPn77ulemYjckISEBKZPn05AQACbNm2iTJkyREa6uukU01JAHunIyVhi4rNxS/Q/9A6bR/fwheSPO8OayvVo8M00aN4ctL6qSKZ0vfPu44EkXIu0x6d/OSKSFtZamjRpwpo1a6hevTqTJk2ia9eu5FCneY+WOyKC9xfP5aHNy8melMSiqk0IbvgQW2+tzB7NzhXJ1K6nPccLwNPALFyXPr82xoyz1n6cUcWJyLX9888/TJ48mWeffZbs2bPz3HPPkT9/fs3ezAzWrYOAAB6bMZP4bF7MrNmScQ06sa9QScA1c1NEMrfrGVF7EmhorT0NYIwJAH4FFNREHHDw4EFGjBjBZ599xunTp/H19aVp06Y89thjTpcmaWEtrFwJw4a5Wm3kz8/Ozg/zZKkOHMhZIGW33N5eDGpT9SoHEpHM4Hr+d9sAiec9TkzeJiI30alTp+jduzcVK1Zk1KhRPPjgg2zcuJGmTZs6XZqkxbkWGw0aQIsWsGmTK6zt38+h/n15+fGm6oEmkgVdz4jaF7hmfc7GFdA6AhMzpCoRucSRI0coXrw4efPmZf369fTt25eBAwdSvnx5p0uTtDh7Fr76ytVi448/oHJlGDcOevSAXLlSdlMPNJGs6XpmfQ43xoTy76zPnmldOkpErs5ay5IlSxg2bBibNm1i//795M2bl99//x0vtWHwKBf3QXu9SUnar54PI0bA4cNQty588w106qQWGyKS4ppBzRgThWuWZ8qm856z1tr8l75KRNIiISGBGTNmEBAQwIYNGyhVqhRvvvlmyuQAhTTPMic8gte/dfVBK3r6BF1/nMzdQxbA2dOuy5yTJ7v+VIsNEbnINYOatdbnZhQiIv8KCwujW7du3H777XzxxRd069ZNLTY8WNCSHRQ9FsEza77lkY3L8E5MYFHVO/m2ZXcmjOnndHki4saut4+aiGSAEydOEBwczJkzZ3jvvfdo1KgRK1asoGnTpmqx4enCw3ntyyHct2MVCdmyMatGC8Y36MSewqU0G0tErklBTcRBERERjBw5krFjxxIdHU3Hjh2x1mKMoXnz5k6XJzfKWggNdS2SvmQJ9+TMw/gGDzLB7wGO5SucsltJ9UETkWtQUBNxyFdffcWTTz5JYmIijz76KP7+/tSuXdvpsiQtkpJg7lxXW401a6BYMXj/fUKbPMCoZftT1uoE9UETkdRRUBO5idasWUPevHm54447aNy4MU8//TQvv/wyFStWdLo0SYu4OPj6awgMhB07oGJFCA6Gnj0hd27aAfE+BS6Y9TmoTVW12xCRa1JQE8lg1lqWLl1KQEAAK1eupGvXrkyZMoXKlSvzySefOF2eXIdLWmz8pxTtf/vO1WIjIgLq1IFp0+ChhyD7hf+8qg+aiNwIBTWRDDRv3jyGDBlCeHg4JUuW5MMPP+SZZ55xuiy5Aee32Ch85iSP/vgVd72zAGKjoXlzmDABWrdWiw0RSVcKaiLpLDY2lhw5cpAtWzbWrl3LmTNnmDBhAt27dydnzpxOlyc3KGjJDooci+Dp32fz6MZl5EiIZ8ltjZnVsjuff9Lf6fJEJJNSUBNJJ5GRkQQHBzNq1CjGjh1Lx44def3113n77bfVYsPTbdzIoK+G0n7bTySZbMy+oznjGnZiV5EyarEhIhlKQU0kjQ4dOpTSYiMqKoo2bdpQqpTrXqRc563VKB7GWvjpJ9cMzkWLaJ0jNxP9HmBC/Qc44lM0ZTe12BCRjKSgJpIG1lpatmzJjh07eOSRR3j11Vfx9fV1uixJi6QkmD/fFdBWr4ZbboH33mPlfzoyYvkBtdgQkZtKQU3kOoWFhREcHMyYMWPInTs3wcHBlC5dmkqVKjldmqTSxbM3B7WpSsc7boGpU11Nardtg/Ll4ZNP4Ikn/m2xkb+gWmyIyE2loCaSCtZavv/+e4YNG8b3339PgQIFeOaZZ2jYsCFNmzZ1ujy5DufP3gQ4cfQftvkPpfXG+eQ5chhq1YIpU+Dhh9ViQ0Qcp6Amcg3//PMPrVu3Zu3atZQoUYLAwED69OlD/vz5nS5NbkDQkh3ExCdS6MxJeq39jp7r5lMwNprwCrXwXfA5tG2rFhsi4jYU1EQuIzY2lnXr1nHnnXdSqFAhbrvtNvr27UuPHj3UYsPDmX37ePv32XTZsJTcCWdZWqURYxs+RHipauy57z6nyxMRuYCCmsh5Tp48ydixYxk5ciQnT57k4MGDFC5cmClTpjhdmqTVpk0QGEjolClYDHOqN+ezhp3YWbQsAKU0e1NE3JCCmghw7Ngxhg8fzqeffsqpU6do1aoVr732GoUKFXK6NEmrn392zeBcsADy5mVv1yd5uvBd7MlTOGUXzd4UEXeloCZZWlJSEtmyZePo0aMEBQXRqVMnXnvtNerWret0aXIdLpnF2aoKHQ9vcAW0VaugSBEYOhQGDKBy4cK8cLlZn5okICJuSEFNsqR169YREBCAl5cXU6ZM4Y477uDAgQOUKFHC6dLkOp0/izN7YgINfl7AHcNnwbF9UK4cfPwx9O4NefKkvEazN0XEU2hdG8kyzrXYaN26NfXq1WPx4sVUqFABay2AQpqHClqyA06fplfYPELHPc2IBcNJtPDOI6/Bn3/Cs89eENJERDyJRtQkyxg1ahQvvfQSt956K8OGDaNv374UKFDA6bIkLf7+m87fTaDnuu8oHHOKNaWr81arfqyoVB9jDG97eztdoYhImiioSaZ19uxZJk+ezG233UbTpk159NFHyZMnD48//rjW4PR0+/fD8OEwfjwvnTnDssoNGNuwM2tLV0/ZRWtwikhmoKAmmc6pU6f47LPPGDFiBIcPH6ZPnz40bdqUEiVK8MwzzzhdnqTF1q0QGAghIa7HXbvy/f09eX5jvNbgFJFMSUFNMpURI0bwzjvvcPLkSVq2bMnkyZNp0aKF02XJdbjsOpyx+10zOOfNc91vNmAADBwIZcvSAvigsmZxikjmpKAmHm/Xrl2ULl2anDlz4uXlRatWrfD398fP7//t3Xd4VFX+x/H3N9RggIihKAFEERFQRAKKK5qgUsSVooKILgqKIAqLCrbddXf9rYSu2GiLIougIE0QiKIRdaVEkV5ViqGJECQIJGTO748JGNhQU+bOzOf1PDxkZs4MX57DTT7ce8/3xAW6NDlLx+3D6RyXfbOA2Dd6w9aVUK4c/P3v/sUBF1xw3Pu0ilNEQpVWfUrQWrp0Kffccw81a9ZkQvalsF69ejF58mSFtCA1aN46Mg5ncMfqZOa89ThvT/kHF6bt5OXbH/Xfl/bCC/8T0kREQpnOqElQcc6RnJxMYmIiSUlJlC5dmieffJKWLVsGujTJq99+45b57/PQkulU2beTdTFVeaJVH2ZecRNZRYry5/POCHiZ2gAAHw5JREFUC3SFIiKFTkFNgk6fPn3YsWMH/fv3p3v37kRHRwe6JMmLPXvgjTfglVf4x+7dpFS+gr/f0o1PL22IM/9Jf+3DKSLhSkFNPO3w4cNMmDCBESNGkJSURHR0NFOmTKFy5cpERuqHd7DIdYFAeQfDhsHIkXDgALRqxYJ2XXnkx0it4BQRyeb5e9TM7GozW2hm35lZipk1CnRNUvD279/P4MGDueSSS+jatStHjhwhNTUVgBo1aiikBZGjCwRS0w7igJIb15P14IP4ql8Cr7wCbdrA8uUwaxY3dmlL/3ZXUjk6EsN/Jq1/uyu1UEBEwlYwnFEbCPzDOTfHzG7Lfhwf2JKkIO3evZuaNWuyd+9emjZtyltvvcWtt96KmQW6NDkHg+at42BmFvVT19J90RRu3bCIw0WL80HD27n73aFw8cXHjdcKThGR3wVDUHNAmeyvywLbAliLFJDt27czbtw4OnfuTExMDE8++SS33norjRrpBGpQc44a337JkIWTuW7rStJKRvHq9R0Y1+CP7C1VlrtPCGkiInK8YAhqfwbmmdlg/Jdqrw9wPZKPli1bxoABA3jvvfcoWbIkrVu3Jjo6mueffz7QpUleHDkCkyfDgAGMW7aMbaVjeLHpQ0ys15zfivsvW2uBgIjI6ZlzLtA1YGafAJVyeel54Gbgc+fcB2bWHujmnLsll8/oBnQDqFixYoNJkyad8s9MT08nKioqz7XLudmyZQuvv/46ixcvJjIykhYtWtCxY0fKly8f6NKEMzs+0g5msnPfITKyfBQvEkHFsiUpF+Gj0ty5VHnvPSK3b+dA1aqsu7M9i66+jiNFfv9/YYQZlc+PJDpSm6afCX2/8hbNh7eEwnwkJCR845zLtQGoJ4LaqZjZPiDaOefMf5PSPudcmVO9Jy4uzqWkpJzyc5OTk4mPj8+/QuW0fD4fe/bsISYmhq1bt9K4cWN69OjBo48+yrJlyzQfHnK64+O4HQSAMofSeXDZHHosn03JPbvh2mvh2Wfhj3+EiIjcV33qPrQzpu9X3qL58JZQmA8zO2lQC4ZLn9uAm4BkoCmwIaDVyFnLyMhgwoQJDBw4kNjYWJKSkoiNjWXz5s0UKVIk0OXJOTi6QKDC/l/omjKDe7+bQ+mMg3xdsxGNp06BG2+EHIs/tEBAROTcBENQexh4xcyKAofIvrwp3peens7o0aMZOnQoP/30E/Xq1aNr164AmJlCWhAr+f0GEhdPpe2qTynq8zGrVhNGXHcnaytcwo833RTo8kREQobng5pz7kugQaDrkLM3YsQI+vbtS3x8PGPGjKFZs2ZqsRHsliyBxEQ+njaNjCLFeO+q5oxu1Jat0f5bTLVAQEQkf3k+qEnw2LRpE0OGDKFJkya0b9+ehx9+mBtuuIHrrrsu0KXJWZi+NJWdO/bz4DOz/feTNatJm92rYcAA+PRTiI5mQ5fH6VKmManFSx97n3YQEBHJf57fmUC8b/ny5dx3333UqFGDkSNHsn79egDKli2rkBZkji4SyMjyYb4s6n89j5qt4qF5c1i7FgYNgi1buHzMK/S9v4l2EBARKWA6oyZ50rt3b4YPH05UVBS9e/emT58+xMbGBrosOUeD5q3Dd/AgdeZ/zKfTpnNx2na+LxfLS3c+xXMT/g9KlDg2VgsEREQKnoKanBWfz8esWbNISEigdOnSNGnShPLly/Poo49Srly5QJcnebFvH3fMfYcuKTMofyCN7y68jEcSniPpsuvAInguR0gTEZHCoaAmZyQjI4N3332XgQMHsmbNGl577TV69uzJXXfdFejS5AydtJfZ9u3w8svw5ps8vX8/n1e/hl86tuWJI1cfa7GhRQIiIoGhoCan5PP5GD58OEOGDOGnn37iqquuYsKECbRv3z7QpclZOLFBbWraQUaMmcvV2z/l4tlT/Fs+tW/PZ3c8wKOrHY/WOgwr/CFNiwRERAJHQU1ydfDgQSIjI4mIiGDmzJlceumljBo1ihYtWqjFRhA62qAWoO6OjXRfOIXb1n1FZtGi8FBXeOopuPRSEoD+S1PZue5bDLSLgIhIgCmoyXE2b97M0KFDeeedd1i5ciWVK1dm5syZQb+PWrjbtvc3/rB5Gd0XTqHJ5u/4tXgp3rzuLt6Ou4Mlr9533Ng29SuTvG8DPybGB6ZYERE5RkFNAFi5ciUDBw5k4sSJANx33334fD4AhbRglpUFU6cyZ8JfqJW6nl3nnU//+Ad49+qW7C9xnu49ExHxOAU1YefOndSvX58SJUrw+OOP06dPH6pUqRLosuQsnLhQ4On4atyxfL6/79nGjcRWrc7fbuvFe1fEc7hocUD3nomIBAMFtTDk8/mYPXs2X3/9NS+99BIVK1Zk4sSJJCQkcMEFFwS6PDlLORcKRB3+jdvnTaHxv2ZA+l5o0AAmTyaqbVuuWb6D+bmt+hQREc9SUAsjmZmZTJw4kYEDB7Jq1SqqVavGM888Q5kyZdRmI4gNmreO89J281jKTO5f+hFlDh/gi2pX82KH5xg++sljLTbUoFZEJPgoqIWJxYsXc/fdd7Nlyxbq1q3L+PHj6dChA8WKFQt0aZIX339Pj/cGc/eKTyiWdYSPLv8DI667i5WVamDAcK3QFREJagpqIeyXX35hx44d1KlThxo1alCrVi1ef/11WrVqpRYbQSTXRrXs8m+SPnky7a0IU+o2ZVSjdmwq9/sZs4u0UEBEJOgpqIWgLVu2MHToUEaPHk3t2rVZsmQJ5cqVY968eYEuTc7ScY1qnaPqskXEjO4LP3wLpUvDU08x/+YOvPjFrmN90kALBUREQoWCWghZu3YtiYmJTJgwAYBOnTrRt2/fAFcleTFo3joOZWTSfP1CeiyazNXbN/DzedG82fwhekwaBNHRtAQOlz/J9lAiIhLUFNRCgM/nIyIiggULFjB58mR69uzJE088QdWqVQNdmuTF4cPc8PkMui2eyqV7fmJT9IU817wnH9S9mYyixekRHX1sqBYKiIiEJgW1IOWc46OPPmLAgAG0b9+exx57jM6dO3PnnXeqxUYQyfX+sxplYORIGDaMAdu2saLipfS842nmXH49vogigDZJFxEJFwpqQSYzM5P33nuPgQMHsmLFCqpUqULZsmUBKFGiBCVKlAhwhXKmTtwo/VDqdnY8PpqM5XMovn8fNG3KV38ZzEOp0Rw84jv2Pt1/JiISPhTUgkynTp2YPHkytWvX5p133uGee+5Ri40gdXSj9CppO3h48TTar/iY4kcy+bzuDSSMHQING/IH/Juk6/4zEZHwpKDmcXv27OGNN96gW7duVKhQgV69enH//ffTqlUrIiIiAl2e5EHZ9avot/ADbl/7BVkWwdTsFhs/XhDLjw0bHhun+89ERMKXgppHbd26lWHDhjFq1CgOHDhAlSpV6Ny5MzfccEOgS5O8cA4WLIDERD6aO5f04pGMadiGsXF3sLN0DKD7z0RE5HcKah7j8/l4+OGHeeedd3DO0bFjR/r168eVV14Z6NIkL3w+mDkTEhNh0SIoX57VPZ/mwZIN2Fm01LFhuv9MRERy0rUzj1i/fj0AERERZGZm0qNHDzZu3Mj48eMV0oJZRga89RbUqQNt28KuXfDGG7B5M7VfS+TZTtdTOToSw38mrX+7K3WZU0REjtEZtQByzjF37lwSExP54osvWL16NbVq1WLcuHHa4ilI5Npeo35lSE+H0aNhyBBITYV69WDiRLjrLij6+2Gn+89ERORUFNQC4MiRI0yePJnExESWL19ObGwsQ4cOJTY2FkAhLUic2F4jNe0gA8d/Qa03F1Jryjuwdy/Ex8OYMdC8OWheRUTkLCmoBcCePXvo0qUL1atX56233uLee++lePHigS5LztLR9hoAlfft4qEl07hnWRIlsjKgTRt4+mm49toAVykiIsFMQa0Q7N27lzfeeIMlS5Ywbdo0KlSowKJFi6hbt65abASxbWkHufznTTyy6APuWP05Potgep14Rl57J/NHdw90eSIiEgIU1ApQamoqw4YNY+TIkaSnp9OyZUsOHDhAVFQUV111VaDLk7z48kv+M/1F/rBuEQeKlWRsXGvGxrVmR5kYtdcQEZF8o6BWQObPn0/Lli3x+Xx06NCBfv36Ua9evUCXJXnh88GsWTBgAPz3v8RFl+OVm/7E2Hot2RdZGlB7DRERyV+67paPFi1aRFJSEgCNGzemV69ebNiwgQkTJiikBbPMTBg3Dq68Elq39q/ifPVVSqRupdqwl4i6sILaa4iISIHQGbU8cs6RlJREYmIiycnJNGzYkGbNmlGqVCkGDx4c6PLkDOXaZqNmtH/F5pAhsHWrP6j95z/Qvj1k76/apn4pBTMRESkwOqOWBx9//DHXXHMNLVq0YMOGDQwZMoT58+cHuiw5S0fbbKSmHcQBv23bwU+9nyajciz8+c9QvTrMng3LlkGnTsdCmoiISEHTGbU82LdvH4cOHWLs2LF06tRJLTaC1NE2Gxf9uouHFk/nnuXzKJV5mAW1r+fGOYOhceNAlygiImFKQS0P2rVrR7t27dRiI8idt2EtQxb7W2wAzKgdz4hr7+T7mKr8qJAmIiIBpKCWBwpoweXE+9D6V9rPjVP/TdKHH/JbsRK8c83tjGnYhu1lygOozYaIiAScgpqEhWPbPWUcIeGHFHosnEyjn1ZzOPp8fuj+JF1KNWJ7sfOOjVebDRER8QIFNQkLQz9aRYulH/PIog+otXszqaXL8/ebu/HFjXcw/2+tePpkm6uLiIgEkOeDmpnVA0YAUcAmoJNz7teAFiXB48ABGDuWiQNfpPKvP7Mupip9Wj3Bh1fcyJEiRbHf/MPa1K+sYCYiIp7j+aAGjAGecs59bmZdgL7AXwNck3jdL7/A66/D8OHwyy/srlaXv97ag88ujcPZ7/cWXqT70ERExMOC4W74msCC7K8/Bu4MYC3idVu3Qp8+ULUqvPCCv7XGF1/w47S5fH1F4+NCmu5DExERrwuGM2qrgNbAdOBuoEpgyxFPWrMGBg707xzgHNx7L/TrB3XrAtAme5juQxMRkWBizrlA14CZfQJUyuWl54F1wHDgAmAm0Ms5d0Eun9EN6AZQsWLFBpMmTTrln5menk5UVFQeK5f8cq7zUWb1aqq++y4xX31FVokSbG/Viq13383hSrn9c5IzpePDWzQf3qL58JZQmI+EhIRvnHNxub3miaB2psysJvAf51yjU42Li4tzKSkpp/ys5ORk4uPj87E6yYuzmg/nYO5cGDAAPv8cypWDxx6Dxx+HmJgCrTNc6PjwFs2Ht2g+vCUU5sPMThrUPH/p08wqOOd2mVkE8Bf8K0AlDJzYoLbfzZfSesNX/oC2fDnExsLQofDwwxDk/5sSERHJjeeDGtDRzHpmfz0VeCuQxUjhONagNjOLkpmHaPrpLBokToN9O+GKK+Dtt6FjR9D+qiIiEsI8H9Scc68ArwS6Dilcg+ato9j+fXT9dhYPfPMhMb/t49uLLufVOx5jwNvPgbbvEhGRMOD5oCZhKDWVB6a+Rsdlc4nKOMinl8Tx5nV3sSS2DmbGAIU0EREJEwpq4h1r18KgQTB+PA9mZTGrVhNGXHsXaytUPzZEDWpFRCScKKhJoTlxccCxPmaLFlHnb3+DL7+EEiXgkUf4tEUnnl2YxsHMrGPvV4NaEREJNwpqUihyLg4ASN37G7MGvc0fNsymfMrXRJcuDc8/72+xUaECzYD+F2mjdBERCW8KalIoBs1bx8HMLIr4srht7Zd0X/QBdXb9wK4yMTBkCAtr1aLJbbcd9x5tlC4iIuFOQU0KxS8/p9Fp5Xy6LZ5KtbQdfF8ulr4tezOjTjzrn2hDVnJyoEsUERHxHAU1KVhpafDmm/x31EDKpafx3YU1eSmhC0mXXYezCCprcYCIiMhJKahJwdi2DV5+GUaMgP37ybw+nj9d0pIFF9UGM0CLA0RERE5HDakkf61f79/SqXp1GDIEbr8dli6l4lef0e6JTlQ+vxQGVI6OpH+7K3UPmoiIyCnojJrkj5QU/x6cH3zgb7HRtSs89RRccsmxIVocICIicnYU1OSM5NoD7eqLYP58SEz0/162LDz7LPTqBRUrBrpkERGRoKegJqd1Yg+07XvSSf6/N7hpzSzOX7McLrzQfzate3coUybA1YqIiIQOBTU5raM90EocyaDdyk/ptvgDqu/dzpYLKnP+6NFw//3+y50iIiKSrxTU5LT279xN96Vz6JIygwoH9rKs0mX0aP0MSTUb8/1DdwS6PBERkZCloCYnt2MHvPwy/33zNaIOH2DBxfXp/cen+LrqVWCmHmgiIiIFTEFN/tfGjTB4MLz9NmRmknZzKzpXbcY3MdWPDVEPNBERkYKnPmryu2+/hQ4d4PLL/SHtgQdg3Tpik2Zyf8+2VI6OVA80ERGRQqQzauHOOfjsM/+qzaQk/6rNfv2gd2+oVOnYMPVAExERKXwKauEqKwumT/f3QEtJ8YeyAQPgkUf8/dBEREQk4BTUws3hwzB+PAwa5N/uqUYNGDkS/vQnKFky0NWJiIhIDgpqISjXXQQuLe0PZMOGwfbtcM018P770K4dFCkS6JJFREQkFwpqIebEXQQOp25jx+OjyFw2h2Lpv8LNN8O4cXDLLWAW4GpFRETkVBTUQszRXQSqpO2g2+Kp3L3iE4ofySS5bhOavjUE4uICXaKIiIicIQW1EFN23SqeXjSFVmu/5EhEBFPrNGXUtXeyqVxlflRIExERCSoKaqHAOfj8cxgwgI/mzmV/8UhGN2zD2LjW7Cp9AYB2ERAREQlCCmrBzOeDGTP8bTUWLYIKFVj12DM8WLIBu4r8Hsy0i4CIiEhw0s4EwSgjA8aOhTp1/Ks2f/4Z3nwTNm2izqv9ee7extpFQEREJATojFow2b8fRo3yt9hITYWrr4aJE+Guu6Do71OpXQRERERCg4JaMNi1C4YPh9dfh7Q0SEiAf/8bmjVTiw0REZEQpqDmISc2qn2hbiTN5k7wh7LDh6FtW3j6aWjUKNClioiISCFQUPOInI1qa+36ke4fTqHpmi/wFSlCxJ/uh759oVatQJcpIiIihUhBzSMGzV3LlT8so/vCKTT9IYX04pGMjWvN7Js7MKN/h0CXJyIiIgGgoBZoPh98+CHDX3uGBtvW8ktkGQY3uY/x9VuxL7I05gJdoIiIiASKglqgZGT4V2wOGABr1nDh+ZX4663dmXzlLRwqVvLYsIvUqFZERCRsqY9aYUtPh5dfhho14IEHoFgxmDCBJfMWMuXa1seFNDWqFRERCW86o1ZYdu+GV1/1/9q7F268EUaMgJYtwYzWgCta9LhVn32bX65+aCIiImFMQa2gbd4MQ4bAmDFw8CC0bu1vsdG48f8MVaNaERERyUlBraCsWAEDB/rvQzODTp2gXz+oXTvQlYmIiEiQ8MQ9amZ2t5mtMjOfmcWd8NqzZrbRzNaZWfNA1XjGvvwSbr8drroKpk2DXr3ghx/g7bcV0kREROSseOWM2kqgHTAy55NmVhu4B6gDXAR8YmY1nXNZhV/iKfh8MHu2fwXnV19BTAz885/QsyeUKxfo6kRERCRIeSKoOefWANj/7lvZGpjknDsM/GhmG4FGwNeFW+FJZGbCpEn+gLZqFVSr5t+Ts2tXKFUq0NWJiIhIkDPnvNNR1cySgaeccynZj18DFjrn/pP9+N/AHOfclFze2w3oBlCxYsUGkyZNOuWflZ6eTlRU1DnXGrllC/X69aPkzp2kV6/O1o4d2ZWQgCvqiewbdPI6H5K/NB/eovnwFs2Ht4TCfCQkJHzjnIvL7bVCSxVm9glQKZeXnnfOzcjr5zvnRgGjAOLi4lx8fPwpxycnJ3O6MaeUmQmzZkHnzkTddhtXmHHFuX9a2MvzfEi+0nx4i+bDWzQf3hLq81FoQc05d8s5vC0VqJLjcWz2c4FXrBi8/36gqxAREZEQ5vXrdDOBd81sKP7FBJcBiwNbkt/0palqTisiIiIFyhNBzczaAq8C5YHZZvadc665c26Vmb0PrAaOAD29sOJz+tJUnp26goOZ/lJS0w7y7NQVAAprIiIikm880UfNOTfNORfrnCvhnKvonGue47V/Oecudc5d7pybE8g6jxo0b92xkHbUwcwsBs1bF6CKREREJBR5IqgFm21pB8/qeREREZFzoaB2Di6Kjjyr50VERETOhYLaOejb/HIiixU57rnIYkXo2/zyAFUkIiIiocgTiwmCzdEFA1r1KSIiIgVJQe0ctalfWcFMRERECpQufYqIiIh4lIKaiIiIiEcpqImIiIh4lIKaiIiIiEcpqImIiIh4lIKaiIiIiEcpqImIiIh4lIKaiIiIiEcpqImIiIh4lIKaiIiIiEeZcy7QNeQ7M/sZ2HyaYTHA7kIoR86M5sNbNB/eovnwFs2Ht4TCfFRzzpXP7YWQDGpnwsxSnHNxga5D/DQf3qL58BbNh7doPrwl1OdDlz5FREREPEpBTURERMSjwjmojQp0AXIczYe3aD68RfPhLZoPbwnp+Qjbe9REREREvC6cz6iJiIiIeFrYBTUze9LMnJnFnOT1zma2IftX58KuL1yY2YtmttzMvjOzJDO76CTjsrLHfGdmMwu7znBxFvOh46MQmNkgM1ubPSfTzCz6JOM2mdmK7HlLKew6w8VZzEcLM1tnZhvN7JnCrjNcmNndZrbKzHxmdtLVnqFyfITVpU8zqwKMAWoBDZxzu094vRyQAsQBDvgme9zewq411JlZGefcr9lf9wJqO+e65zIu3TkXVegFhpkzmQ8dH4XHzJoBnzrnjpjZAADn3NO5jNsExJ34vUzy15nMh5kVAdYDtwI/AUuAjs651YVdb6gzsysAHzASeMo5l2sIC5XjI9zOqA0D+uH/IZOb5sDHzrk92T98PgZaFFZx4eRoKMh2HiefEykEZzgfOj4KiXMuyTl3JPvhQiA2kPWEuzOcj0bARufcD865DGAS0Lqwagwnzrk1zrl1ga6jsIRNUDOz1kCqc27ZKYZVBrbmePxT9nNSAMzsX2a2FegE/O0kw0qaWYqZLTSzNoVYXtg5g/nQ8REYXYA5J3nNAUlm9o2ZdSvEmsLZyeZDx4f3hMTxUTTQBeQnM/sEqJTLS88DzwHNCrei8Haq+XDOzXDOPQ88b2bPAo8BL+QytppzLtXMLgE+NbMVzrnvC7DskJVP8yH55HTzkT3meeAIMOEkH3ND9vFRAfjYzNY65xYUTMWhLZ/mQ/LJmczHGQiJ4yOkgppz7pbcnjezK4HqwDIzA/9p62/NrJFzbkeOoalAfI7HsUBygRQbBk42H7mYAHxELsHAOZea/fsPZpYM1AcU1M5BPsyHjo98dLr5MLMHgNuBm91JbibOcXzsMrNp+C+/Bd0PIi/Ih/lIBarkeByb/Zycg7P4fnWqzwiJ4yMsLn0651Y45yo45y52zl2M/5T0NSeENIB5QDMzO9/Mzsd/Bm5eIZcbFszsshwPWwNrcxlzvpmVyP46BvgDoBtzC8CZzAc6PgqNmbXAfz/tHc65304y5jwzK330a/zzsbLwqgwfZzIf+BcPXGZm1c2sOHAPoJXqARJKx0dYBLVTMbM4MxsD4JzbA7yI/4BbAvwz+znJf4lmttLMluM/gHrD8fMBXAGkmNky4DMgUSuoCsxp50PHR6F6DSiN/3LNd2Y2AsDMLjKzj7LHVAS+zD4+FgOznXNzA1NuyDvtfGQvNngM/39e1gDvO+dWBargUGZmbc3sJ6AxMNvM5mU/H5LHR1i15xAREREJJmF/Rk1ERETEqxTURERERDxKQU1ERETEoxTURERERDxKQU1ERETEoxTURERERDxKQU1ERETEoxTURCQsmVl6Ht4baWafm1mR7Mc3m9n4E8YUN7MFZhZSW/WJSOFSUBMROXtdgKnOuazsx/WApTkHOOcygPlAh0KuTURCiIKaiIQ1M3sie/uslWb25xzP/9XM1pnZl2Y20cyeyvG2TsCMHI/rAZWyz6BtMbOjG0pPzx4rInJOdEpeRMKWmTUAHgSuBQxYZGaf4//eeCf+AFYM+Bb4Jvs9xYFLnHObcnxUPWCCc+5GM2uLP5x9gn8T6IaF87cRkVCkoCYi4ewGYJpz7gCAmU0FmuC/2jDDOXcIOGRmH+Z4TwyQdvSBmRUDLgCGZD9V7OjrzrksM8sws9LOuf0F/rcRkZCjS58iImfnIFAyx+MrgGXOOV/246vwn0k7qgRwqJBqE5EQo6AmIuHsC6CNmZUys/OAttnPfQX80cxKmlkUcPvRNzjn9gJFzOxoWKsHLMvxmVcBywHM7AJgt3Mus+D/KiISinTpU0TClnPuWzN7G1ic/dQY59xSADObiT9w7QRWAPtyvDUJ/2XTT/AHtSU5XqvL72fUEoDZBVW/iIQ+c84FugYREc8xsyjnXLqZlQIWAN2cc99mv3YN0Mc5d/9pPmMq8Ixzbn3BVywioUhn1EREcjfKzGrjvx9t3NGQBsfOxH1mZkVy9FI7Tvbq0OkKaSKSFzqjJiIiIuJRWkwgIiIi4lEKaiIiIiIepaAmIiIi4lEKaiIiIiIepaAmIiIi4lEKaiIiIiIepaAmIiIi4lEKaiIiIiIe9f8Cc1RAvSprywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "log_h_range = np.log(h_range)\n",
    "log_errors = np.log(errors)\n",
    "\n",
    "# perform linear regression to get line of best fit for the log scales:\n",
    "lm = linregress(log_h_range,log_errors)\n",
    "print(\"slope: %f    intercept: %f    r value: %f     p value: %f\" % (lm.slope, lm.intercept, lm.rvalue, lm.pvalue))\n",
    "\n",
    "# plot line of best fit with the data:\n",
    "x = np.linspace(np.min(log_h_range),np.max(log_h_range),100)\n",
    "plt.scatter(log_h_range,log_errors,label='data')\n",
    "plt.plot(x,lm.intercept + lm.slope * x, 'r', label='fitted line')\n",
    "plt.plot(x,-1+2*x,'--',c='black',label='reference line')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.xlabel('$\\log(h)$')\n",
    "plt.ylabel('$\\log(W)$')\n",
    "plt.title('Plot of $\\log(W)$ vs $\\log(h)$')\n",
    "# plt.savefig('1D_prior_results.png',dpi=300,bbox_inches='tight',facecolor=\"w\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9181463-ec9e-4114-84b0-4bb3fea7ad3f",
   "metadata": {},
   "source": [
    "We can see from the results above that we indeed obtain a slope of $p=2$ for this example."
   ]
  }
 ],
 "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
}