{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intrusive Galerkin\n",
    "\n",
    "When talking about polynomial chaos expansions, there are typically two\n",
    "categories methods that are used: non-intrusive and intrusive methods. The\n",
    "distinction between the two categories lies in how one tries to solve the\n",
    "problem at hand. In the intrusive methods, the core problem formulation,\n",
    "often in the form of some governing equations to solve is reformulated to\n",
    "target a polynomial chaos expansion. In the case of the non-intrusive methods\n",
    "a solver for deterministic case is used in combination of some form of\n",
    "collocation method to fit to the expansion.\n",
    "\n",
    "The ``chaospy`` toolbox caters for the most part to the non-intrusive\n",
    "methods. However it is still possible to use the toolbox to solve intrusive\n",
    "formulation. It just requires that the user to do more of the mathematics\n",
    "them selves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem revisited\n",
    "\n",
    "This section uses the same example as the [problem\n",
    "formulation](./problem_formulation.ipynb). To reiterate the problem\n",
    "formulation:\n",
    "\n",
    "$$\n",
    "   \\frac{d}{dt} u(t) = -\\beta\\ u(t) \\qquad u(0) = \\alpha \\qquad t \\in [0, 10]\n",
    "$$\n",
    "\n",
    "Here $\\alpha$ is initial condition and $\\beta$ is the exponential growth\n",
    "rate. They are both unknown hyper parameters which can be described through a\n",
    "joint probability distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:36.621214Z",
     "iopub.status.busy": "2021-05-18T10:56:36.620878Z",
     "iopub.status.idle": "2021-05-18T10:56:38.588823Z",
     "shell.execute_reply": "2021-05-18T10:56:38.589079Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "J(Normal(mu=1.5, sigma=0.2), Uniform(lower=0.1, upper=0.2))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from problem_formulation import joint\n",
    "\n",
    "joint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the parameters are positional defined as $\\alpha$ and $\\beta$\n",
    "respectively.\n",
    "\n",
    "First step of intrusive Galerkin's method, we will first assume that the\n",
    "solution $u(t)$ can be expressed as the sum:\n",
    "\n",
    "$$\n",
    "   u(t; \\alpha, \\beta) = \\sum_{n=0}^N c_n(t)\\ \\Phi_n(\\alpha, \\beta)\n",
    "$$\n",
    "\n",
    "Here $\\Phi_n$ are orthogonal polynomials and $c_n$ Fourier coefficients. We\n",
    "do not know what the latter is yet, but the former we can construct from\n",
    "distribution alone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.591144Z",
     "iopub.status.busy": "2021-05-18T10:56:38.590881Z",
     "iopub.status.idle": "2021-05-18T10:56:38.635427Z",
     "shell.execute_reply": "2021-05-18T10:56:38.635655Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q1-0.15, q0-1.5, q1**2-0.3*q1+0.0216666667])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import chaospy\n",
    "\n",
    "polynomial_expansion = chaospy.generate_expansion(3, joint)\n",
    "polynomial_expansion[:4].round(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note again, that the variables are here defined positional. $\\alpha$ and\n",
    "$\\beta$ corresponds to positions 0 and 1, which again corresponds to the\n",
    "polynomial variables `q0` and `q1` respectively.\n",
    "\n",
    "The second step of the method is to fill in the assumed solution into the\n",
    "equations we are trying to solve the following two equations:\n",
    "\n",
    "$$\n",
    "   \\frac{d}{dt} \\sum_{n=0}^N c_n\\ \\Phi_n = -\\beta \\sum_{n=0}^N c_n \\qquad\n",
    "   \\sum_{n=0}^N c_n(0)\\ \\Phi_n = \\alpha\n",
    "$$\n",
    "\n",
    "The third step is to take the inner product of each side of both equations\n",
    "against the polynomial $\\Phi_k$ for $k=0,\\cdots,N$. For the first equation,\n",
    "this will have the following form:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "   \\left\\langle \\frac{d}{dt} \\sum_{n=0}^N c_n \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle -\\beta \\sum_{n=0}^N c_n\\Phi_n, \\Phi_k \\right\\rangle \\\\\n",
    "   \\left\\langle \\sum_{n=0}^N c_n(0)\\ \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle \\alpha, \\Phi_k \\right\\rangle \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Let us define the first equation as the main equation, and the latter as the\n",
    "initial condition equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main equation\n",
    "\n",
    "We start by simplifying the equation. A lot of collapsing of the sums is\n",
    "possible because of the orthogonality property of the polynomials $\\langle\n",
    "\\Phi_i, \\Phi_j\\rangle$ for $i \\neq j$.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "   \\left\\langle \\frac{d}{dt} \\sum_{n=0}^N c_n \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle -\\beta \\sum_{n=0}^N c_n\\Phi_n, \\Phi_k \\right\\rangle \\\\\n",
    "   \\sum_{n=0}^N \\frac{d}{dt} c_n \\left\\langle \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   -\\sum_{n=0}^N c_n \\left\\langle \\beta\\ \\Phi_n, \\Phi_n \\right\\rangle \\\\\n",
    "      \\frac{d}{dt} c_k \\left\\langle \\Phi_k, \\Phi_k \\right\\rangle &=\n",
    "   -\\sum_{n=0}^N c_n \\left\\langle \\beta\\ \\Phi_n, \\Phi_k \\right\\rangle \\\\\n",
    "   \\frac{d}{dt} c_k &=\n",
    "   -\\sum_{n=0}^N c_n\n",
    "   \\frac{\n",
    "      \\left\\langle \\beta\\ \\Phi_n, \\Phi_k \\right\\rangle\n",
    "   }{\n",
    "      \\left\\langle \\Phi_k, \\Phi_k \\right\\rangle\n",
    "   }\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Or equivalent, using probability notation:\n",
    "\n",
    "$$\n",
    "   \\frac{d}{dt} c_k =\n",
    "   -\\sum_{n=0}^N c_n\n",
    "   \\frac{\n",
    "      \\mbox E\\left( \\beta\\ \\Phi_n \\Phi_k \\right)\n",
    "   }{\n",
    "      \\mbox E\\left( \\Phi_k \\Phi_k \\right)\n",
    "   }\n",
    "$$\n",
    "\n",
    "This is a set of linear equations. To solve them in practice, we need to\n",
    "formulate the right-hand-side as a function. To start we create variables to\n",
    "deal with the fact that $\\alpha$ and $\\beta$ are part of the equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.637580Z",
     "iopub.status.busy": "2021-05-18T10:56:38.637307Z",
     "iopub.status.idle": "2021-05-18T10:56:38.644573Z",
     "shell.execute_reply": "2021-05-18T10:56:38.644325Z"
    }
   },
   "outputs": [],
   "source": [
    "alpha, beta = chaospy.variable(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As above, these two variables are defined positional to correspond to both\n",
    "the distribution and polynomial.\n",
    "\n",
    "From the simplified equation above, it can be observed that the fraction of\n",
    "expected values doesn't depend on neither $c$ nor $t$, and can therefore be\n",
    "pre-computed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the denominator $\\mathbb E[\\beta\\Phi_n\\Phi_k]$, since there are both $\\Phi_k$\n",
    "and $\\Phi_n$ terms, the full expression can be defined as a two-dimensional\n",
    "tensor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.646490Z",
     "iopub.status.busy": "2021-05-18T10:56:38.646231Z",
     "iopub.status.idle": "2021-05-18T10:56:38.658420Z",
     "shell.execute_reply": "2021-05-18T10:56:38.658142Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(10,), (10, 10)]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi_phi = chaospy.outer(\n",
    "    polynomial_expansion, polynomial_expansion)\n",
    "[polynomial_expansion.shape, phi_phi.shape]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This allows us to calculate the full expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.660422Z",
     "iopub.status.busy": "2021-05-18T10:56:38.660161Z",
     "iopub.status.idle": "2021-05-18T10:56:38.835901Z",
     "shell.execute_reply": "2021-05-18T10:56:38.835661Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.50e-01, 8.33e-04, 0.00e+00],\n",
       "       [8.33e-04, 1.25e-04, 0.00e+00],\n",
       "       [0.00e+00, 0.00e+00, 6.00e-03]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_beta_phi_phi = chaospy.E(beta*phi_phi, joint)\n",
    "e_beta_phi_phi[:3, :3].round(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the numerator $\\mbox E(\\Phi_k\\Phi_k)$, it is worth noting that these are\n",
    "the square of the norms $\\|\\Phi_k\\|^2$. We could calculate them the same way,\n",
    "but choose not to. Calculating the norms is often numerically unstable, and\n",
    "it is better to retrieve them from three-terms-recursion process. In\n",
    "``chaospy`` this can be extracted during the creation of the orthogonal\n",
    "polynomials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.837936Z",
     "iopub.status.busy": "2021-05-18T10:56:38.837677Z",
     "iopub.status.idle": "2021-05-18T10:56:38.880238Z",
     "shell.execute_reply": "2021-05-18T10:56:38.880452Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.00e+00, 8.33e-04, 4.00e-02, 1.00e-06])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, norms = chaospy.generate_expansion(3, joint, retall=True)\n",
    "norms[:4].round(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having all terms in place, we can create a function for the right-hand-side\n",
    "of the equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.882537Z",
     "iopub.status.busy": "2021-05-18T10:56:38.882272Z",
     "iopub.status.idle": "2021-05-18T10:56:38.889040Z",
     "shell.execute_reply": "2021-05-18T10:56:38.888753Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "\n",
    "def right_hand_side(c, t):\n",
    "    return -numpy.sum(c*e_beta_phi_phi, -1)/norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial conditions\n",
    "\n",
    "The equation associated with the initial condition can be reformulated as\n",
    "follows:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "   \\left\\langle \\sum_{n=0}^N c_n(0)\\ \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle \\alpha, \\Phi_k \\right\\rangle \\\\\n",
    "   \\sum_{n=0}^N c_n(0) \\left\\langle \\Phi_n, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle \\alpha, \\Phi_k \\right\\rangle \\\\\n",
    "   c_k(0) \\left\\langle \\Phi_k, \\Phi_k \\right\\rangle &=\n",
    "   \\left\\langle \\alpha, \\Phi_k \\right\\rangle \\\\\n",
    "   c_k(0) &=\n",
    "   \\frac{\n",
    "      \\left\\langle \\alpha, \\Phi_k \\right\\rangle\n",
    "   }{\n",
    "      \\left\\langle \\Phi_k, \\Phi_k \\right\\rangle\n",
    "   }\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Or equivalently:\n",
    "\n",
    "$$\n",
    "   c_k(0) =\n",
    "   \\frac{\n",
    "      \\mbox E\\left( \\alpha\\ \\Phi_k \\right)\n",
    "   }{\n",
    "      \\mbox E\\left( \\Phi_k \\Phi_k \\right)\n",
    "   }\n",
    "$$\n",
    "\n",
    "Using the same logic as for the first equation we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.891057Z",
     "iopub.status.busy": "2021-05-18T10:56:38.890749Z",
     "iopub.status.idle": "2021-05-18T10:56:38.900930Z",
     "shell.execute_reply": "2021-05-18T10:56:38.901150Z"
    }
   },
   "outputs": [],
   "source": [
    "e_alpha_phi = chaospy.E(alpha*polynomial_expansion, joint)\n",
    "initial_condition = e_alpha_phi/norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equation solving\n",
    "\n",
    "With the right-hand-side for both the main set of equations and the initial\n",
    "conditions, it should be straight forward to solve the equations numerically.\n",
    "For example using `scipy.integrate.odeint`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.903384Z",
     "iopub.status.busy": "2021-05-18T10:56:38.903124Z",
     "iopub.status.idle": "2021-05-18T10:56:38.911000Z",
     "shell.execute_reply": "2021-05-18T10:56:38.910740Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 10)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "coordinates = numpy.linspace(0, 10, 1000)\n",
    "coefficients = odeint(func=right_hand_side,\n",
    "                      y0=initial_condition, t=coordinates)\n",
    "\n",
    "coefficients.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These coefficients can then be used to construct the approximation for $u$\n",
    "using the assumption about the solutions form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.913220Z",
     "iopub.status.busy": "2021-05-18T10:56:38.912894Z",
     "iopub.status.idle": "2021-05-18T10:56:38.936483Z",
     "shell.execute_reply": "2021-05-18T10:56:38.936234Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([q0, -0.01*q0*q1+q0, -0.02*q0*q1+q0, -0.03*q0*q1+q0])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_approx = chaospy.sum(polynomial_expansion*coefficients, -1)\n",
    "u_approx[:4].round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, this can be used to calculate statistical properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.938754Z",
     "iopub.status.busy": "2021-05-18T10:56:38.938439Z",
     "iopub.status.idle": "2021-05-18T10:56:38.953843Z",
     "shell.execute_reply": "2021-05-18T10:56:38.954059Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1.5     , 1.49775 , 1.495503, 1.493259, 1.491019]),\n",
       " array([0.04    , 0.03988 , 0.039761, 0.039643, 0.039525]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean = chaospy.E(u_approx, joint)\n",
    "variance = chaospy.Var(u_approx, joint)\n",
    "\n",
    "mean[:5].round(6), variance[:5].round(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:38.956492Z",
     "iopub.status.busy": "2021-05-18T10:56:38.956225Z",
     "iopub.status.idle": "2021-05-18T10:56:39.024700Z",
     "shell.execute_reply": "2021-05-18T10:56:39.024917Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEECAYAAAA2xHO4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/DElEQVR4nO3dd3xkZ33o/8/0Php1bdP2fWzvuq8BY4MxNiFgyE1IciGElsC9AVIgBEJyQxKbUBxaAuRCSEJ+3OSGNODGAcc0Y1NcsL2212W9j7evpFWXZkbT6++Pc0ZnVqsdjXY1RdL3/XrtSzpnzsx5/Fgz33na97GVy2WEEEKI87G3ugBCCCHamwQKIYQQNUmgEEIIUZMECiGEEDVJoBBCCFGTBAohhBA1ORvxokqpncBHgMeBzcC01vrDC67xAp8CRoDdwJ1a6+fNx94EXA0UgWNa6y81opxCCCGW1qgWRRfwL1rrT2qt3wO8QSl17YJr3guc1lp/HPgL4MsASqnNwPuB92utfx94h1Jqd4PKKYQQYgkNCRRa60e11nctuE9ywWW3AQ+Z1z8NXKmUCgOvBA5orSsrAR8CXtWIcgohhFhaQ7qeqimlfgH4jtb68IKH+oC5quO4ee58589SKpXKxaKsKgdwOOwUi6VWF6MtSF1YpC4sUhcWl8thW+5zGhoolFI3AzdjdDMtNAGEqo7D5rkJYNeC80cXPrlYLBONplasrKtZJOKXujBJXVikLixSF5be3tDSFy3QsFlPSqnbMLqR3gMMKKWuV0p1md1LAHcD15vXXg4c1FrHge8A1yqlKlHveuCeRpVTCCFEbQ0JFObA9b8CLwLuA+4CFPAHwLvNyz4LbFVKfQj4PeDtAFrrYYzZUH+hlPo08Hda6yONKKcQQoil2VZr9th8vliWpqRBmtUWqQuL1IVF6sLS2xta9hiFLLgTQghRkwQKIYQQNUmgEEIIUZMECiGEEDVJoBBCCFGTBAohhBA1SaAQQghRkwQKIYQQNUmgEEIIUZMECiGEEDVJoBBCCFGTBAohhBA1SaAQQghRkwQKIYQQNUmgEEIIUZMECiGEEDVJoBBCCFGTBAohhBA1SaAQQghRk7MRL6qUGgA+Alyptb5ukcdvB14B5M1Te4A3aq3vV0o9DGTM80Wt9S2NKKMQQoj6NCRQADcCdwFXnefxx4BPaq2TSik7cDfwQ/Oxb2utb29QuYQQQixTQwKF1vprSqmX1Xj8W1WHPwfcrbUum8eXK6U+CPiAR7XWdzeijEIIIerTqBbFcrwVeHPV8Z9rrR9RSjmAHyml5rTWP1r4JIfDTiTib1oh25nUhUXqwiJ1YZG6uDgtDRRKqauAo1rrROWc1voR82dRKfVj4GbgnEBRLJaIRlPNKmpbi0T8UhcmqQuL1IVF6sLS2xta9nOaNutJKRVQSvUuOP1bwF9VXXOJUurtVY/vBo41o3xCCCEW16hZTzdhdCdtUEp9CPg08DbgcuCd5jX9gFdrfarqqXHgNqXURiAMDAFfbUQZhRBC1MdWLpeXvqoN5fPFsjQlDdKstkhdWKQuLFIXlt7ekG25z1m1C+4S2UKriyCEEOvCqg0UTw5FmZjLtroYQgix5q3aQFEqlzk8keDEdIrV2n0mhBCrwaoNFBVD0TSHxuYoFEutLooQQqxJqz5QAEyn8jw5EieVK7a6KEIIseasiUABkMoXeWIkxlQy1+qiCCHEmrJmAgVAsVTm0NgcJ2dkGpwQQqyUNRUoKk7PpnlmNC7jFkIIsQLWZKAAmEnleWIkJusthBDiIq3ZQAGQzpc4OBJnXNZbCCHEBVvTgQKgWC6jJxIcmUxQkvUWQgixbGs+UFSMxrMcHImTycsUWiGEWI51EygA5rIFHh+WKbRCCLEc6ypQABTMKbTHp5OS+kMIIeqw7gJFxXA0I11RQghRh3UbKADi0hUlhBBLWteBAqyuqGNTSZkVJYQQi1j3gaJiJJbhyZEYaemKEkKIs0igqJLIFnl8OMZYPNPqogghRNuQQLFAsVTm+ckkh8fnKJSkK0oIIZyNeFGl1ADwEeBKrfV1izz+MuAvgah56m6t9SfNx24FXgdMAGWt9R2L3aPRU1snEjnimSiqL0iHz9XQewkhRDtrSKAAbgTuAq6qcc17tdb3V59QSvmBvwb2aq2zSqmvK6Vu0Vrfu/DJn/vBMV53xQAhT6P+EyBTKPHUmThbOn1s7fRhs9kadi8hhGhXDfmU1Vp/zWw11PJmpdR+IAz8rdZ6CLgeOKW1rmTxewC4DTgnUDw1EuPkdJK3Xr+VvRs7VrD055rOlSjEc1y2IYTf3bjAdKEcDjuRiL/VxWgLUhcWqQuL1MXFadWn3iHgz7TWJ5VSe4HvKaUuA/qAuarr4ua5c6j+IHo8wefvO8ZLdnRx295+XI7GDbkkk1nGp5Ps6PGzIext2H0uRCTiJxqVzZpA6qKa1IVF6sLS2xta9nOW/cmqlLp52XdZQGs9obU+af7+LBABtmCMS1T/V4TNc+d4/yv28OrL+rDb4MfHZ/jsD483fLZSsVzmyGSSZ0fj5AqyKZIQYn1YskWhlHoV8C4gCNiAQWDncm+klAoAfq31pFLqD4C/0VrPKKW6ADcwDkwCW5VSHrP76QbgC4u9nt1u45Y9vezuDfJPjw0zGs/yF/cf5+f2DfDi7Z0NHU+YTuWJD0fZ3ROgJ+hp2H2EEKId1NP19EfAezE+xG3AW5d6glLqJuDNwAal1IeATwNvAy4H3gmcAD6rlDoEXAa8RWudMZ/7LuBzSqlJ4KnFBrKrDXb6eN/NO/iPp8Z45HSUbzw1yuGJBK+/eiPBBg5054tlDo0n6Evm2NUTwNnAbi8hhGgl21LTTJVSn9da/3bV8Q6t9fGGl2wJ9x4aK0cXdDUdHInx70+eIZ0vEfI4+ZVrN6H6gg0vi8dhZ3dfgC6/u+H3Woz0v1qkLixSFxapC0tvb2jZ3S31fOUeUEr9E/C8efwS4Nbl3qgZrtzUwWCnn68eGOb4dIq/efAUN2w3Bro9zsZ9488WSzwzOsdAyMOObr+0LoQQa0o9n2gDwHeBk+a/aOOKc/E6/S7edeM2Xn1ZHw6bjQdOzPCZ+45xaqbx3ybG5rIcGIoxk5JstEKItaOeFsWvaa2PVg6UUvc0sDwrwm4zBrov6Qvy1cdHGItn+fyPTnDLnh5ecUkvTnvjWxf9IQ87pXUhhFgD6vkUSyil/kkp9bRS6h8xBrRXhU0RH++9aQcv29UNwPefn+KzPzzBaKzxSf/GzdaF7HUhhFjt6gkUH8VIx/FW4FvAnQ0t0QpzOey8dt8A775xG11+F2diGf7ih8e578hUw/efyBZLHBqb47nxOVl3IYRYterpejqstf438/fHlVI7GlmgRtnRE+D3bt7JN58d5+GTs3zr2XGeHZvjDddsoifQ2NlKk4kc0XSeHd0B+kOy7kIIsbrU06LYZS6KQynVA6zKQAHgdTn45as28o4XDRL2OjkxneLTPzjGT45PN7x1kS+W0RMJnj4j+3QLIVaXegLF/wEOKqWiwAHg7xtaoia4dCDE+1++k6s2hckVS/y/p8b44k9OMpnILv3kizSbznNgKMZwNN3wVOlCCLESlgwUWusHtdZbgF1a663Akw0vVRME3E7efN0W3vqCLYQ8To5Pp/jUD45xfxPGLorlMsenUzw+HGMuU2jovYQQ4mKdN1CY3UwopV6qlHopcJn583PNKlwzXLExzAdu2cm1WzoolMp889lxPv+jE03ZDjWZK/LkSIyjk0kKRRnsFkK0p1qD2f8IvAr4LGe3Ii5vZIFaIeB28sZrN3PVpg6+9uQZTs+m+cz9x/kZ1cvNu3tw2Bs3I7gMnIlnmEpm2dEdoE8Gu4UQbea8gUJr/Srz19/RWv8YQCllx9hcaE26bCDEB27ZxbeeGefhU7Pc89wET52J8/qrN7Ip4mvovXPFMocnEozPZdnVG8DncjT0fkIIUa96BrMvXfD7rzWoLG3B53Lwy1dv5DdevJVOn4uRWIa//OFx/uvQOPkmdA8Zg91RTk6nKJZksFsI0XrnbVEopcIYGwpdopQaNE8ngXUx+rqnL8gHbtnJ3c9O8MCJGe59foqDI3F+8aoN7OltbEbaUhlOR9OMJ7Ls7Ak0fJ2HEELUUmuM4hcw9pDYBlyFkbqjAHy70YVqFx6ng9dduYGrN3fw70+eYXwuy5ceOMX+LR28dt9AQ/e7AMgWjJXdXX4XO3ukO0oI0Rr17EfxCq3195pUnrr94NBYebYJM5MqCqUS9x+Z5nt6kkKpjN/t4Of2DbB/S0dDd9OrsNtgU4ePLZ0+nAsG1yXXvkXqwiJ1YZG6sFzIfhT1rKM4K0gopd623Js0wuWbOvA38Ru2027nVtXL+1++k109AVK5Iv/y+Ah//cCppizUK5VhKJrmwOko43ONv58QQlTU06J4I/CnQA+Qxtj3uqsJZaspny+WZ2aTDM2mGYqmaea4b7lc5sBQjLueGSOVK+K027hV9XLz7u6GpjCvFvY42dkTIOR1yrelKlIXFqkLi9SFpSEtCuDFGLOdPqG13kwbZY+122xs7fJz7ZYIEZ+rafe12WzsH4zwwVt2sX9LhEKpzLefm+Az9x3n6GSyKWWIZws8MRJDTyTISu4oIUQD1TMaO6y1LimlvObxpqWeoJQaAD4CXKm1vm6Rx98GvAg4BlwDfF5r/aD52MNAZfChqLW+Zan7+VwOrtgYZnwuy/HpJPlic5oXQXNf7v2DHXztyVHG57J88YGTXL25g9fu7aejCcFrfC7LQydm6HLZ2dThbejiQCHE+lRPoHiBUuq1QFYpdS/1bVx0I8YeFled5/FNwHu11hml1AuBv8Na8f1trfXtddzjHP0hD91+F8enU4w1sR9/d2+QD7x8J/cdneb7epInhmMcGpvjlZf0cuOO7oZ/eBdLZU7OpBiNZ9je5ZfV3UKIFbXkGEWFUsoB3AY8pLWerOP6lwGf0lrvX+K664HPaK2vN4+/DjwC+IBHtdZ3L/a8UqlcLtZYABdN5dDjcySyze2WmUpk+fcDwxwcjgGwscPLG67bwp7+UMPu6bDbzlqcF/Y62d0XJOJff+svHA47tf4u1hOpC4vUhcXlciz7m2vdCwG01kXgP5VSv84KpRpXStmA9wDvqzr951rrR8zA9COl1JzW+kcLn1sslpYcnNoT8TISy3BqJk2xSSm9fTZ4y/7NHNrcwX88NcqZWIbPfP8I12zu4LX7+gl7V747KhDwkExaLahkMsvodJLugJsd3f51tf5CBi0tUhcWqQtLb+/yv7QuGSiUUndgLLwrYnQ7hVmBQGEGiU8CX9FaP1Q5r7V+xPxZVEr9GLgZOCdQ1MNms7E54qMn4Ob4dKqp+1dfNhBid2+AHxyZ4gfPT/F4pTvq0j5u2N7VlLGE6WSOmWSODWEvg50+3M7mzMgSQqwt9XxyXANs01rv0FpvB95+ITdSSgWUUr3m7w6MrLTf1Fp/Wyn1i+b5S5RS1a+/G2PA+6J4XQ4uGwixb0MIn6t5H5Yuh51XXtLH79+yi0v7g2QKJe56eozP3HeM5ycSTSlDJTvto0NRTs1I/ighxPLV0/X0BODFWEMBxmdPTUqpm4A3AxuUUh8CPo3RKrkceCdGS+LngSuUUgA7ga8DceA2pdRGjJbLEPDVuv9rltDldxPZ4mJoNs1wNNO07qjugJt3XL+VZ0fn+I+nRxmby/KlB0+xdyDEz+3rpyfY+MHnYqnMqdk0o/EMg51+BsIe7E1YUS6EWP3qWXD32xhrJ8Yxu5601t1NKFtN+XyxfDF9jpl8kWPTKaab2B0FkC+W+NGxab6vp8gVSzhsNl66q4tb9/TivcCxhIVjFPXwuexs7Vx7M6SkL9oidWGRurBcyIK7eloU/w3YqLWOASil3rrcm7Qjr8vB3oEQM6kcx6aSpPPNmRHhcti5ZU8v1w1G+K9DEzx6Osp9R6Z57HSMV13Wx3WDkaZ800/nSxyeSDAUTbO920/XOpwhJYSoTz0d9g9VgoTpZIPK0hJdfjfXbomwrcuPo4ldMWGvizdcs4n33LSDrV0+5rIF/u2JM3z2/uMcn27O6m4wtmN9ZnSOgyMxoul80+4rhFg96ul6ehjoB05gdD0Naq13NqFsNV1s19NiMvkiJ2ZSTCaa2x1VLpd5YjjGtw6NE0sb231ctSnMbZf101XHXhQX0vV0Pp0+F9u6/IS8jU2h3ijSxWCRurBIXVga1fV0Eni9+buNC5z1tBp4XQ4u7Q+xIZzn2FSSZK45i/VsNhvXbImwd0OY+45Mcd+RKZ4cifP06Bw37jDGL/zu5qyFmE3nmR2J0R1ws7XT1/A9N4QQ7a/uldkVSqlLtNaHG1SeujWiRVGtXC4zGs9yciZFoclTSmdTOe45NMEBc3W3z+XgFaqHG7Z34XSc21u4ki2KhXoCbrZ2+Qi4V0fAkG+OFqkLi9SF5UJaFOcNFEqp12itv6WU+pMFD71Ua33rhRRwJTU6UFQUiiVOzqYZjWWWnhe8woajab75zDhHp4wxiy6/i1df1s+Vm8JnDXg3MlBU9AbdDHa2f8CQDwSL1IVF6sKy0mnGK1lfrwZOVf2LLrtkq5jTYWdXT4BrtnTQ2cRU5gCbIz7eecNW3vGiQQZCHmZSef7vY8N87ocnODbVvAFvgMlEjseHYjw3Pkcyty62TRdCmOoZzN6itR6qOr5Ba/1Aw0u2hGa1KBaaTuY4Pt286bQVxVKZR09H+c7hCeIZ44P6soEQr9nbz46BcMNbFAv1BIwWRruNYcg3R4vUhUXqwtKowewPAr8FoJRyAr8HtDxQtEp3wE2n38WZWIbTs+mmjV847DZetK2Tqzd38MOjU9x3ZJpDY3McHp/j+h3d3Lyzm05/81o8U8kcU8kc3QE3gxHfqp0lJYRYWj0tiruAp4F/wkgG6NdaX9mEstXUqhbFWWUoloy0GC0Yv4hn8nz38CQ/PTVLqWwEkhu2d3LLnt6WfMvv8rvYEvE1ZbOmWuSbo0XqwiJ1YVnRwewKpVQ/8EbgzzHyNX1Taz13IQVcSe0QKCpSuSInppNMp5q/YG0ykeX7R6Z57NQsAB6nnZfu7OamXd0tSS/e4XUy2Omjs0UrveUDwSJ1YZG6sDRqz+xHgFcCOzCyuf6f5d5krfO7HezdEOaKDWGCnuZ+OPcGPbzjxu2872U7uLQ/SLZQ4nt6ko999wj3HZki3+TNWmKZAk+PzvH4cJSpRJblTr8WQrSfegLFT4BXaa2HtdZ3YMx8EouI+F1cvamDPb0BPIusd2ikTREf77h+K7/5km3s6PaTyhf51rPjfOx7R3joxEzT04snskUOjSd4bCjKWDxDSQKGEKtWPV1PHeavu4AjQEJr3fI9Bdup62kxxVKZkViGoWi64R/SC9dRlMtl9ESC/zo0wUgsAxiD8D+jerl6c0dTNk1ayOOwsyniZSDsxdnA+0sXg0XqwiJ1YWnUrKebgP8NzAIR4N3At5Z7o/XGYbcx2OljQ9jDqdm0+a26Ofe22Wxc0h9iT1+Qp8/Euee5CSYTOf758RG+pyd5RQsCRrZY4vh0itOzaTaEvWzs8OKRHfeEWBXqeaf+DLBTa30FsAd4TWOLtLa4zAV7126J0Bts7gCv3Wbjyk0dfODlu3jDNZvoDriZShoB4xP3HuWx09Gmd0kVSmWGomkePT2LnkjI4j0hVoF6AsUprXUOQGudAU43tkhrk89MOHj1pg46mrzmwGG3cd1ghA/ecm7A+OQPjnJgKNr0MYRSGcbnshwYivH0mTizqeZm7BVC1K+eT6ydSqn3Accxtizd2tgirW0hr5MrN3Uwk8pxYjrVtAy1YAWMazZ3cGAoyvf1JJOJHF89cHaXVLO3SJ1N55lN5wm4HWzq8NIXkm1ahWgn9QxmB4H/BVwBPAncqbVONL5otbX7YHa9JuaMDLWZwoXPD7jQpIDFUnk+YFTWgPQG3dy6p3WD3gBuh40NYS8bwl7cyxzHkEFLi9SFRerC0qgFd18HPqq1frzeF1VKDQAfAa7UWl+3yON24GPAHLAN+LLW+mHzsVuB1wETQNmcknuOtRIoAErlMmPxLEOzabIXsO7hYrPHFktlHjMDxowZMLr8Ll6+u4f9gxFcTZ7qW2G3GetENnV4615tLh8IFqkLi9SFpVEL7krVQUIpVc+KshuBuzA2OlrMfwfCWuuPYuSS+gellEMp5Qf+GvhdrfXtwBVKqVvquN+qZrfZ2NjhZf9ghO1d/oZOH12Mw27jhVs7+YNbd/Pfr95IT8DNTCrP1w6O8rHvHeH+I1NkC83rIquojGM8Phzj4EhMFvAJ0SL1BIoHlFKXVB1/cKknaK2/htFaOJ/bgIfMa2eADLAXuB5j8Lzy9fgB89p1wWG3saXTxwsGIwxGfE3v+qkEjA/euos37d/MhrCHeKbAN58d5yPfOcJ3D0+QatEspVimwKHxBI+cjjI0m276inMh1rN62vN/BvyOUgqMFkIYo9voYvRxdiCJm+d6z3P+HA6HnUjEf5HFaF893UEuK5Q4NZ1kJJapOY3VYbcRCHhW9P43Ki837OnlmTNx7nlmjONTSb5zeJL7j07z0t093Hppf8sSAE5ki0xPp+kLedjS6SPktcqx1v8ulkPqwiJ1cXHqCRR3aK0/VTlQSv3SCtx3AghVHYfNc+XznD9HsVhaF32OvR4H4S4fQ7NpxuYWX7TXyB3udkS8vPuGrRyfTvH95yd5fiLJ956b4D49yQu2RnjZrh66A61JABify3D0TIywx8mGDi+9QTddnYF18XdRD+mXt0hdWHp7Q0tftMCSgUJr/Sml1GXAZcCzZrfSsimlAhgpyieBu4GXAv+olOoCvMCzgAfYqpTymN1PNwBfuJD7rSUep51dvQE2R7ycnk0zkcg2bZU3GCu9d/YE2NkT4PRsmh88P8nTo3M8eGKWh07McsWmMC/b1cNgp695haoSzxaITyQ4Pm1jd75M0FbG24LMuUKsVfXMevog8BbgGEa+p3/QWt+5xHNuMp/zs8AXgU8Dvw5crrV+pznr6eNAChgE/rZq1tMrgF8CJoH8epj1tFyZfJFTs2km5rKUac6e2QuNxTPcd2SKx4dj80FrR7efm3f3cEl/sGXrIAIBD6lklk6/i4Gwl26/C9s6XZMh36ItUheWRk2P/Qet9Vuqjr+qtX7jBZRvRa3nQFGRzhc5PZsmWYZEormBoiKazvOT49M8dGJ2fi1If8jDTTu7uWZLR9On1i4Mmh6nnYGQh4Hw+sstJR+OFqkLS6Omxx5ecHwQwOyOEi3kczlQfUFetL2L/pDnvHORGynic/GavQP88Sv38Nq9/XT4nIzPZfm3J8/w0e8e4ft6smUzpQCyBWMXwkdOzfLs2BwzqZxMsRVimeppUdwNTGOk8NgBuIFDwEu11rc2vITnIS0KS+XbUjpfZKgFYxjViqUyT47EuP/INGfiRopzt8POC7ZGeMmOLnqCKzs7a6F6uuE8Tjv9IQ8DIc+aHsuQb9EWqQtLo9KM54B7zd9PVJ2PLvdmorF8Lgd7+oIMdvoYiqYZn2t+wHDYbVy7xcgndWQyyX1Hp3h+IslPjs/wwPEZLu0P8pKd3ezuDbRs7CBbKHF6Ns3QbNoYywh56Aq4Jb+UEOdRT6D4Ha310MKTSql7F7tYtJ7X5WB3b5DBTj9D0ebuhVFhs9nY0xdkT1+QM7EMPzo2zRPDMQ6NJzg0nmAg5OHGnV1cuzmy7HxOK6UMzKTyzKTyuBw2+oIeBsIeAu7mZvcVot3V0/V0FfA3GNNjDwHvXE7ep0aRrifLUs3qXKHEcCzNaDzb9P0nqs1lCzx8cpYHT8wQzxjjFj6Xgxdti3DD9i46/Re/HmMlZoCFPU76wx56A26cLcpztRKku8UidWFp1Kynu4CPAkcBBXxIa93ytBoSKCz1vgkKxRLDsQxnYhkKLQwYhVKJp0bi/Pj4DKdn04Cx5P/yjWFesqOL7d3+C+6WWsmpwg6bje6Ai/6Qh4hv9U2zlQ9Hi9SFpVFjFAe01o+Yvz+klHp0uTcR7cHpsLOty8/miI+xeIaRaOaCstVedDnsdq7ZEuGaLRFOzaT48fEZDo7EeOpMnKfOxNnU4eWG7V1cvbmjZd1SAMVymYlEjolEzhgAD3roC3nwu9fuALgQi6mnRfFRjMHsyqyn64B/Bt6ttf6DhpfwPKRFYbnQb0ulcpmJuSxD0TTpfGuT7MXSeR48McNDJ2fnN3PyOu1cNxjhenP6bz2asfgw7HHSF/LQF2zvrin5Fm2RurA0quvpBHBykYcGtdY7l3vDlSKBwnKxb4JyucxUMsdQNE0i2/x04tXyxRIHR+I8eGKGU2a3FMCungAv3t7Jvg3hmll1m7lK3W6DTr+b/pCHLr+r7WZNyYejRerC0qiup9/RWn9z4Uml1KuXezPRnmw2G71BD71BD9FUnqFomtl0viVlcTns7B+MsH8wwnA0zUMnZnl8OMrRqSRHp5KEvU5euLWTF23rJNKi7LUVpTJMJ3NMJ3M47TZ6g276gp6WZdUVolGWbFEspJR6m9b6K40pTv2kRWFpxLelRLbAcDTNZCJHq9cxp/NFDgxFefDELONzRmvBBuzdEOLF27vY3RuY/zbfirxXC1XGM3pD7pZOtZVv0RapC0ujup7eCPwp0AOkMTLAdl1QCVeQBApLI98EmXyRkViGsbnWTq0Fo4vs+HSKB0/M8NSZ+PzakG6/ixds7eS6wQgbe4ItDxTVAm4HfUEPvUF301eBy4ejRerC0qiupxcDlwIf0Fr/uVLq95ddMrFqeV0OdvYE2NrlN2ZKxTJkC60Z+K5Odx7P5HnkVJSHT84yncpzz3MTfPu5CfZtCrN/cweX9oeavkPgYpK5IidmUpyYSdHhddITdNMb8LR0NpcQy1VPoBjWWpeUUl7zeFMjCyTak9NuY3PEx6YOL5OJHCOxDHPZ1iX7C3td3Kp6efmeHp6fSPLTU7M8OzrH0yNxnh6JE/I42T8Y4YVbI/Q2OL9UvWKZArFMgeNTKTp8LvqCbnpW+aI+sT7U0/X0DeD/w9jT+hWATWv98iaUrSbperK0qlkdS+cZiWWYTrZ+HAOMld9Pjyf48fOTTCRy8+d3dvt54bZOrtgYbnra86XYgE6/i57AygcN6W6xSF1YGjJGUaGUcgC3AQ+Zu9S1lAQKS6vfBJVxjPG5bEtXfIMxmJ1IZDg5k+anp2Y5OBIjVzTK5HXauWZLB9cNdrIl4m27ldZ2G3T6XPQEPXT7XRcdNFr9d9FOpC4sDQ0U7UYChaVd3gSFkrGA70wsQyrfmvUYC2c9ZfJFnhiJ8dOTUYai1rqM/pCH/Vs6uHZLpC2ns9ptxl4fF9PSaJe/i3YgdWGRQLFOteObYCaV40wsw0yquesxak2PPRPL8MipWZ4YjpEwV3/bgN19AfZviXD5hnBbDjLbMIJGb9BNl99ddxnb8e+iVaQuLBIo1ql2fhOk80XONLFbqp51FMVSmcPjCR4bivLs2Nz8tF+P086VG8PsH4ywvdvfdiutwQgaYXP2VLe/9pTbdv67aDapC0tTAoVS6vVa63+t47pbgdcBE0BZa33Hgse/DFSnALkcuFZrfVIpdRIrbciI1vpXF76+BArLangTFEtlJhJGt1Qll1MjLHfBXSpX4MmROI+ejs5nsgXo8ru4dkuE/Vs6Gr4r38UIeZx0m91TC5MVroa/i2aRurCsaKAwczwtfNAGhLXW3bVeVCnlB54C9mqts0qprwNf0FrfW3XNfMBRSoWBr2itX2ce3661vr3WPSRQWFbbmyCWzjMazzCVzK34hkoXszJ7fC7LgaEoB4ZiRKtSmAx2+rh6cwdXbQoT9rbfeEaF3+WgO+CmO+Ai7HWtur+LRpK6sKz0grtPaK2/uPCkUupddbzu9cAprXXlHfsAxoyp+UCxoFXy68DfVx2/xFzYFwLu0Vo/WMc9xSrR4XPR4XORK5QYm8syGm/dIr5q/SEPr76sn5+9tI+jU0keOx3l6TNznJ5Nc3o2zX8+Pcau3gBXb+7gig1hfG2WbjyVL5KKphmKpnE7bGwbKOIuFYn42i9hoVhd6up6Ukr1AAPASa11oo7rfwV4vdb6583jdwAv01q/aZFr7cDdwKu11mXz3Au01o+YLZPHgddorY9WP69UKpeLLdhLoR05HHZWc12Uy2WmkzlGomljTcZFtDIcdtuKphrJFUo8NRzl0VOzPHMmPv/aTruNfZvCXLeti8s3tnbfjPOp1IXDbjO6p4JuutfpqvDV/h5ZSS6XY+VTeCil3gT8EfA08FWl1F6t9UeXeNoERmugImyeW8zPAXdXggRAZaMkrXVKKfUkcAPGDnvzisWSNCVNa6FZ7QS2Bt30exyMzWUZj2cvaFOlRiQFvLQ3wKW9AVK5Ik+fifP4cIxjU0meHIrx5FAMj9PO5RtCXL05wu7eQFukDoGz6yI+l+EE1mB4l99NV8C1bvYHXwvvkZXS2xta+qIF6vkruVJrfalS6oNa6/9QSu2r4zkPAVuVUh6z++kG4AtKqS6goLWOV137VuDNlQOl1C2AS2v9bfPULuBYXf81YtXzuhxs6/KztdPHdCrPWDzDbCrfFiu//W4HL9zWyQu3dRJL5zk4YgSNoWiax4ZiPDYUI+h2sG9jmCs3htnZ0z5Bo6KMlUrkxIyxCLEr4Kbbb3QHSheVWEw9gSJm/qy8V5fsmDVbAu8CPqeUmgSe0lrfq5T6BDAD3AmglLoKOLqgO2sCuF0pdQ2wEfiG1vondf3XiDXDZrPNLzbLFkqMxY0Mtu0wlgHGOMtLd3Xz0l3dTCWyPDEc4/HhGBOJHA+fnOXhk7ME3A4u3xDmik1hdrVh0ADIFEqcMfdRd9hsRPwuuvwuuvxuPOuwi0osrp5cT58HXMAG4AjGVNcPNKFsNcmsJ8t6aVaXy2Wi6Txjc1mmzzNjqpX7UZTLZUbjWQ6OxDh4Js5kVb4pv9vB5RtCXLmpo2lB42LrIuB20OV30+l30eF1tl3Kk+VYL++RejRqPwoH8HbgCuAg8GWtdcu/1kmgsKzHN0G+WGJ8Lsv4XPasdRntsHERGEFjLJ7l4Jk4B0diZyUp9Lsc7NsQ4spNYXb3BhsWNFayLpx2GxGf0droXIWtjfX4HjmfZi24e5XW+p7l3milSaCwrPc3QSJbYCyeZSKRxeNzt0WgqFYulxmby3JwJM5TZ+Lzu/QB+FwO9g6E2LchhOoLruiMpEYGzYDbQaffRZfPTdjnbPuxjfX+Hqm20gvu7uPcBXcAW7XWOxc531QSKCzyJjCUymXyDgdHRqJtMwC+mLF4hoNn4jw1EmesKmi4HDZUX5B9G8JcNhC86BlJzWpdOWw2OnxOOv0uOn3nrhBvB/Iesaz0grtHgC8AbwAeA44DO4BbL6h0QjSY3WajP+zFUwqTK5SYSJzbNdUOBsJeBsJeXnlJH+NzWZ4ZjfPMqLGw75nROZ4ZncNug+3dfvZtCLNvQ4guv7vVxT6vYrnMTCpvJoBM4XHa6fS56PS7iPhcbbcHiFi+esYo/qh63YRS6g6t9Z82vGRLkBaFRb4tWRari0S2wPhclslEdn5vinYUS+d5dmyOp0fjHJ1MnjVYv6nDy74NIfZtCLMh7KlrYLldxmuCHgedPjedPlfLuqnkPWJp1J7ZL1RK7cdY8LYH2L/cmwjRSkGPk6DHyY5uP7PpPBNzWaaTeYptljm5w+fixdu7ePH2LtK5Is9NJHhmNM7h8QQjMWO/8u8cnqTL72LvQIhLB0Ls7Pa3/VaqiWyRRNZILWK3QYfXaGlEfC6CHseqnk21XtQTKP4U+BLGVqjPAvXkehKi7dhsNmNFst9NsVRmKpljMpFty/EMn9vBNZs7uGZzB/liiaOTSaNbaizOTCrPj4/P8OPjM7gddlRfgEv7Q1w6EGzrpIUApTLMpvPMmkkXnXYbHT4XEZ+TiG/9rBRfbWQ/ijVAmtWWC6mLXKHEZDLL5FyOeLbQoJKtjFK5zOnZNM+NzXFoLMGZeOasxzdHvFzWH+KygRB7NnWQTuXO80rtye2wzbc2OnwufDX221gOeY9YGrWO4hLgyxjrKJ4E3qG11hdSwJUkgcIibwLLxdZFJl9kMpFjItF+g+CLmU3lOTw+x6GxOY5MJclXjcGEvU4u6Qty6UCIPb2BmpsctSuP007E66LDbHFc6H+DvEcsjQoUX8boeqqMUfyG1vrXLqiEK0gChUXeBJaVrItUrshEIstUIteyPcCXo9JFdWjcaG1U76nhsNnY3u1H9Qe5pC9Y94B4u6kOHMtpcch7xNKowezDlWyuwMNKqZuXexMhViO/20hQuK3LTyJbYCqRYzKZJZ1veWKCRbkcdi41B7lfd0WZWL7MgZPTHBpLcGomxdGpJEenktz97DghjxPVF0T1B9nTGyDoWR1jA9lCifFElvGEMZvL47DPB40Or6st13CsBfX8dexRSl0NnMDYunRXY4skRPupzJza1m0EjclEjqk2Dho2m41NnV4i7l5u2dNLKlfg+ckkh8cT6IkE8UyBx4aiPDYUxQZsjvhQfQFUf5Ctnf62TGC4mGyxxEQiN58ixeWw0eF1EfZWBsdlVtVKqHeM4u+xxijeLmMU7UWa1ZZm10U7tzTOt46iklKkEjSOT6fO2uzJ67Szu9cIGqo3SFegfRf7LcVhtxH2ONnSH4JcgZDHuWqCYKM0K9dTr9Z6crk3WmkSKCwSKCytrItkzggaU8lcWwyE17vgLlsocXwqiZ5IcHgicVbWW4Buv4vdvUF29QbYvYq6qapV6sJug4DbSdhrdFeFPc51t+Nfowazg8ArsHase63W+peXX7yVJYHCIoHC0i51kc4XmUrkmE62bsrtha7Mnknm0BNGa+PoVPKcltKGsIfdvUF29wbY0e1fFbOpatWFz2Un7DFWjYe9Tvyutd1d1ajB7G9hdDnNmsddy72JEOuNz+VgS6ePLZ0+soUS00kjaMQy+UX30WgnXQE312/v4vrtXZTKZYajGY5MJjgymeTEdIrReJbReJYfHZvGboPBTt984Nja6Wv7leILpfMl0nlrgNxptxHyGEEj7HUS8jhX3X/TSqsnUBzVWr+3cqCU2t244gix9nicdjZ2eNnY4aVQKjOTzDGTyjGTylNo86hht9kY7PQx2Onjlj295IslTs2k5wPHUDTNyRnj3/f0JC6HMQ13d0+AnT0BNkd8q25MoFAqn7V63IYxA64SPELroNWxUD1dT78MBLH2rX6z1vp/NLpgS5GuJ0u7dLe0g9VUF6VymVg6z3Qqz3Qyt+LbvDYjKWAmX+TYVMoIHFNJxuJn38/tsLOty8eOngA7u/0MtqjFsdJ1URkkD3mMwBH2OldNltxGdT39GpAFoubx5cu9iRDiXHabjU6/m06/m109ARLZgpGuO5ljLltou/xTi/G6HOzdEGLvBmMIcy5T4MhUkmPmv8lEjucnkzw/mQSMbp2tXT52dhstjq1dvlXzAVutuKDVAcZssZDZVRUyp1OvttbU+dQTKKa01m+pHJhrKpaklLoVeB0wgbHP9h0LHn8b8E6gkqzmy1rrfzQfexNwNVAEjmmtv1TPPYVYzSprNQY7feSLJXOPhxyzq6CLqiLkdc4nMwSIZ/Icn05xbCrJ8ekUY/Esx6ZSHJtKgZ7EYbOxpdPHzh4/O3sCbOvy4XG2/+D4YjKFEplEbn7W2MIuq6DHuWrXddTT9fR7wONYXU9v01p/eInn+IGngL1a66xS6uvAF7TW91Zd8zbgfq31yQXP3YwxgH611rqslHoUeKPW+kj1ddL1ZFlN3S2NthbrolwuE8sUmDUDR71Tb9tlP4pqiWyBE2bgODadYjSWOavlZLfBxg4v27v8bOv2s73LT4fv4jPitktdOGw2gh7H/JeCkMfZ9NXkjep6ej9wuOp4EKgZKIDrgVNa68r/mQeA24B7F1z3W0qpMcAP/JXWegZ4JXBAa135+3kIeBVwBCHWIZvNyqi6vdtPJl9kNm3sKBdN589aLNfugh4nl28Mc/nGMGDk0zoxneLYdJLjU0mGo5n5fz8+PgNAp881HzS2d/sZCHvafo/u8ymaQT+WsaZMO+w2gm7HfOAIeZ0rljV3pdQTKP5Qa/2VyoHZpbSUPmCu6jhunqv2Q+BurfWkUurVwL8Dt9T5XBwOO5GIv46irH1SF5b1UhcD5s9SqUwsk2faXLORyJ79ARQIeFpTwDoFAtDb6ecFu3oAY3D85HSSY5NJjk0mOD6ZNMYChmM8MRwDwOuyG4PjvUF29gbY1r10Ztx2r4sCEC2UiRbykMwbU3Qr4x0+13zLo1XdVksGiuogYR5/v47XncBaoAcQNs9Vv86JqsMfAP+plHKY1+1a8NyjC29QLJbWXBfDhVqL3S0Xaj3WhQ3ocdvpcXvJFUrMpvNEU3lydhvxucySz283W0IetoQ8vGyHsY5jLJ7lxHSKEzMpTk6nmE3nOTQ6x6FR4/ukDaO7aluXn61dxlTenoD7rA/Vdul6Wo7Y3NnHlW6rgNs5333ldzuW3brq7Q0tfdECjVqL/xCwVSnlMbufbgC+oJTqAgpa67hS6uPAH2utC8Bu4KTWuqiU+g7w20opm9n9dD3w+QaVU4g1xe200x/y0B/yEIn4GR6PE03nmU3liWcKbbf961LsNtv8GpQbdhhrfWPpPCdnUpyYTnFyJjW/TexILMMD5tdPv8sxv/5jsMvHJZscrM7OKsti3VZ2m7G40xjzcBB0GwPmKz0FuWE73CmlXgH8EjAJ5LXWdyilPgHMaK3vVEq9B9iHkZX2cuCzWuuHzee+CWNv7iLw/GKznmQw27Iev0Wfj9SFZWFdlMpl4umCETjSeRKrZAruUrKFEkOzaU7Npjg1k+b0bJq5RdKm9ATcbDUDx2Cnn40dHpz21Tc1tx5ep31+llXA4yTodsx3zzUlKWC7kEBhkQ9Hi9SFZam6KBRLRDNG4Iil822RyHAllMvGGofTs2lOz6Q5NZtmJJY+a/c/MNZ0bOrwmq0OP4MRL90LuqzWEqfdmBRx0+UbGzLrSQixBjkddnoCbnrMNOK5Qolo2phJFcvk2y5ter1sNhtdfjddfjdXbTLWc3h9bo6Oxjg9mzZbHSkmEjlOzRqBBHOGlddpZ3PEx+aIly2dPjZHfHT7XWsieBRK5UVbWvWQQCGEAIzxjb6Qh76QMTsoawaO2CoPHGDMejICgI8XbzfOpXJFhqJpTs+kOD2bZiiaYS5bmN8JsMLnqgoe5mt0rZHgUS8JFEKIRXmqBsbBCBwxs8URzxRWxT7itfjdDmM72L7g/LlYOs9wNMNQNM1wNM2wGTyOTCY5MlkdPBxsjnjZHPGxxfy5loOHBAohRF08C1ocuUKJeCZPNFMgls6TyhVX/eB4h89Fh881n7uqXC4TzxTMwJFhOJpmKJomkS2eEzy8ZpbgTeYsrY0dXgZCnjWRolwChRDigriddnqCHnqCRuAoFEvEswXi6QKxTJ65bKHt995Yis1mmw8e+zYYq8krKVWGo2mGZtMMx4wAksgWOT6d4vi0NYHAboP+kOecABJwr66P3tVVWiFE23I67PODyGBMx01kC8TNuf9zmTy54iqPHJydUqUSPMBIgHgmluGMuabjTCzDZCI3v9HTgaHY/LUdPiebwl4zgPjYGPHS5Xe1bWoSCRRCiIaw22yEvS7CXhebzXPpfJF4pkA8Y45zrIHuqorKf+sl/dbK52yhxFj87OAxGs8QSxeIpRMcGk/MX+t22BkIexgIedgQ9jIQNn4GPa3POCuBQgjRND6XA5/LMT9AXiiVSWQKxLNGiyOeLZyz3mE18zjtbO3ys7XLyj9WKpeZTubmA0flZzxTMNZ+zKbPeo2A2zEfOCrBYyDkaepe5RIohBAt47TbiPhdRPwuwAcYiQHjZvBIZAokcqt/rKOa3WajN+ihN+iZX+cBkMwVGItnGY1nzJ9ZxuIZkrniOVN2wciqOx84zJ99QXdDBs8lUAgh2orXZaSbqMyuKpXLJLNF5rIF498amJq7mIDbyc4eJzt7AvPnyuUy0XSBsXiG0bksYzHj58Rcdn6Hveequq9sQE/QPT+tufKvL+jB7bzwACKBQgjR1uw2M+W21/q4KpSMgfJEVfDIrPCe4+3AZrPR6XfR6Xdx6YA19lEslZlK5owAYrY8RuNZppPGDnuTiRzPjFrpZ21Al9/Fhg4vr9k/uOxySKAQQqw6lbxFkard7wrFktnqKM4HkbUYPMBYaV5pLVy5yTqfL5aYTOQYn8ue9W8ykWU6lWc6lT//i9YggUIIsSY4HXY6/W46q/atqgQPu9fN6KSR6yiTL62ZmVYLuRz2+bUa1QqlElOJHDMSKIQQ4myV4BGJ+Okwu+gLpTLJrDFInswWSeSMabpracB8IafdzkDYe9bsq2U9f4XLI4QQbc1pt1ZbV5TLZZK5ovnPCiBraaruxZBAIYRY92w2m7lLnBOw9tbOFUpG4MgV54NHOr+2Wx+LkUAhhBDn4XbacTvPHvcol8uk8kbgSFVaILnimh04BwkUQgixLDabjYDbeU5iv0KpTLoqcKTMf9ni6g8gEiiEEGIFOO3nrvcAI4CkzAHzZK5IOm/8zK6iFkjDAoVS6lbgdcAEUNZa37Hg8Q8CA8AosB/4E631YfOxk8BJ89IRrfWvNqqcQgjRSE67lRyxWrFUng8a6VyRVN5ogWQK7TcG0pBAoZTyA38N7NVaZ5VSX1dK3aK1vrfqsiDwPq11WSn1euCTwGvNx76itb69EWUTQoh24LBXD6BbyuUymULJGPcwg0c6bwSSVs3CalSL4nrglNY6ax4/ANwGzAcKrfUfV11vBxJVxy9RSv0+EALu0Vo/2KByCiFEW7HZbPNZdhcqFEuk8kXS+RLpXJF0wWyF5EsUy40LIo0KFH3AXNVx3Dx3DqWUG3gr8JtVp/9Qa/2I2TJ5XCn1Gq310ernORx2IpELWzyy1khdWKQuLFIXlvVQF1mz1ZHOWd1Z6XyBdL5E0ezL8rouLDFgowLFBEZroCJsnjuLGSS+CPyR1vpY5bzW+hHzZ0op9SRwA3BWoCgWS0SjKQREIn6pC5PUhUXqwrJe6sIG+AG/2w5uO2CMi2QLJdL5IvkLnIHVqF2/HwK2KqUqK1duAO5WSnUppcIwP47xJeAzWusDSqlfNM/fopT62arX2gUcQwghxAXxOO1EfC56g56lL15EQ1oUZkvgXcDnlFKTwFNa63uVUp8AZoA7gf8L7AO2K6UAAsDXMVoetyulrgE2At/QWv+kEeUUQgixNFu5gQMgjZTPF8vroSlZj/XSrK6H1IVF6sIidWHp7Q0tewPuRnU9CSGEWCMkUAghhKhJAoUQQoiaJFAIIYSoSQKFEEKImiRQCCGEqEkChRBCiJokUAghhKhJAoUQQoiaJFAIIYSoSQKFEEKImiRQCCGEqEkChRBCiJokUAghhKhJAoUQQoiaJFAIIYSoSQKFEEKImiRQCCGEqEkChRBCiJokUAghhKjJ2agXVkrdCrwOmADKWus7FjzuBT4FjAC7gTu11s+bj70JuBooAse01l9qVDmFEELU1pAWhVLKD/w18Lta69uBK5RStyy47L3Aaa31x4G/AL5sPncz8H7g/Vrr3wfeoZTa3YhyCiGEWFqjWhTXA6e01lnz+AHgNuDeqmtuA/4XgNb6aaXUlUqpMPBK4IDWumxe9xDwKuBI9Q1cLoettzfUoOKvPlIXFqkLi9SFReriwjVqjKIPmKs6jpvn6rmmnucKIYRokkYFigmgOnyHzXP1XFPPc4UQQjRJowLFQ8BWpZTHPL4BuFsp1WV2LwHcjdFFhVLqcuCg1joOfAe4VillM6+7HrinQeUUQgixBFu5XF76qguglHoF8EvAJJDXWt+hlPoEMKO1vlMp5cOY9TQK7AI+tmDW036MWU/Py6wnIYRonYYFikZaaurteqGU2gl8BHgc2AxMa60/3NpStY755eOnwHe11u9vdXlaRSmlgF8B0sBNwO1a60daW6rWUEp9ANgGTGFMw3+71jrd0kI1kVJqAOMz4kqt9XXmufMuTTifVbfgrs6pt+tFF/AvWutPaq3fA7xBKXVtqwvVQh8Bnmh1IVpJKeUAPgN8WGv958DbgROtLVVrmB+Sfwj8ttb6T4EAxhfM9eRG4C7AVnXuvSyyNKGWVRcoOP/U23VHa/2o1vquqlN2INmq8rSSUurNGH8L6/JDscp1GB8Kv62U+kPgtRjfptejFJDDmBADEASebV1xmk9r/TXOnkUKxuflQ+bjTwNXVo0dL2o1BgqZPrsIpdQvAN/RWh9udVmaTSl1GXCp1vobrS5LG9iK8WXqK+Y3xpcCb21tkVrDnBzzAeBflVJfAYaBoy0tVHtY9mfoagwUMn12AaXUzcDNwO+2uiwt8gtARin1BxhN7Rcopd7b2iK1TBw4rLWOmcc/AV7WuuK0jlLqKoxAcZvW+m0YLas/aWWZ2sSyP0NXY6BYdOptC8vTUkqp2zBWs78HGFBKXd/iIjWd1vqjWusPa63vxPhgfERr/ZctLlar/BToNscqwGhh1ByoXMM2YcyyLJjHo4C3heVpF+dbmnBeq3XW0zlTb1tcpJYwB65/CDxmngoA/1tr/ZWWFaqFlFK/CPwm4Maoh39ucZFawuyGfDnG+2MQYzB33cz0qTCD5eeADBAF9gHv1VqPtrJczaSUugl4C/CzwBeBT5sPLbo04XxWZaAQQgjRPKux60kIIUQTSaAQQghRkwQKIYQQNUmgEEIIUZMECiGEEDVJoBDiIimlPqGUut/8PayU+uFFvNY2pdTPr1TZhFgJEiiEuHhfqPxiLlx62UW81jbg5y+uOEKsrEbtmS1E21BK/QnGIrwscCXwO8CHMVYs78bIi/SAUmrjec7/K7ADY3Hji4F/A34E3AE8AuSr7vUWjEVeEaXUazGyc34TI2XCZcAbtdYnzVbDfwM0cDnwLoyFYW8DrlJK3Q78C0ZOno8Dz2AsjvqS1vqAUurdGKnl54CtWut3rmytCWGRFoVY05RSrwRepLX+kNb6zzB2S/w08G2t9SeAD2EkjbPVOP9BYCNGyupXAv8FfAn4iPmalZXxaK3/AWMVMFrrb2KkFDmstX4H8P+AXzQvncVYJXwncAB4s9Y6B3wFeFJrfbuZ4PFTVWX6GPB35vP/J0YSyI8D/7CSdSbEQhIoxFp3BVUZQ7XWXzbPHTePx4EOoKfGeYCjWuu81nrOTHewFzhiPnZ8iTJU0iNMYiVjSwB/YqYCvx7orVH+F5oJD38FmFBK2TFaHr+hlHoEWM97kIgmkEAh1rqDwM7KgVLq16vPmZvbRDEyi57vPMDCXDeHgD3m7zuWKMNieXL+DrjLbBF8r+p8EbAppbzmTnUHgXvNlsedwD9prUvAFq31GzGyBv+eUqpriTIIccEk15NY88wxCh/GGMA08A3goxgtgl3Al6vGKBY7/xHgV4E/01r/vfma+4FKt5MLeD3wW0AE+CuMbqrHMXZjfAKjK+uvgE6MbqOfxRijuA+jRVA5PwP8O3AY+BZGoPgwcAxjR8OHtNbfUEr9Lcb+CmUgorV+30rXmxAVEiiEEELUJF1PQgghapJAIYQQoiYJFEIIIWqSQCGEEKImCRRCCCFqkkAhhBCiJgkUQgghavr/AXkj/0JzJ8rJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "\n",
    "pyplot.rc(\"figure\", figsize=[6, 4])\n",
    "\n",
    "pyplot.xlabel(\"coordinates\")\n",
    "pyplot.ylabel(\"model approximation\")\n",
    "pyplot.axis([0, 10, 0, 2])\n",
    "sigma = numpy.sqrt(variance)\n",
    "pyplot.fill_between(coordinates, mean-sigma, mean+sigma, alpha=0.3)\n",
    "pyplot.plot(coordinates, mean)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the true mean and variance as reference, we can also calculate the mean\n",
    "absolute error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:39.027149Z",
     "iopub.status.busy": "2021-05-18T10:56:39.026841Z",
     "iopub.status.idle": "2021-05-18T10:56:39.034361Z",
     "shell.execute_reply": "2021-05-18T10:56:39.034125Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.54252e-11, 1.6888e-08)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from problem_formulation import error_in_mean, error_in_variance\n",
    "\n",
    "error_in_mean(mean).round(16), error_in_variance(variance).round(12)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}