02-BayesianBRDF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Towards Bayesian BRDF linear model understanding\n",
    "### J Gómez-Dans (NCEO & UCL)\n",
    "\n",
    "<div>\n",
    "<img style=\"float: right;\" src=\"https://designabetterbusiness.com/wp-content/uploads/2016/06/uncertainty.jpg\" width=\"250\"/>\n",
    "</div>\n",
    "\n",
    "## Uncertainty\n",
    "\n",
    "We haven't said anything about how uncertain data is. In this case, the solution will not change if we assume that all observations have the same uncertainty associated (e.g. $\\sigma^2\\sim 0.005$, say). But depending on things like SZA, VZA and other stuff, this will change. It's only that we don't know how. But let's see how this would work.\n",
    "\n",
    "Assuming no correlation between uncertainty in observations, we have that \n",
    "$$\n",
    "\\hat{\\rho}^{\\lambda_i}(\\Omega, Omega') = \\sum_{i=1}^{N_{kern}}f_{i}(\\lambda_i)K_{i}(\\Omega, \\Omega') + \\mathcal{N}(0, \\sigma_{obs}),\n",
    "$$\n",
    "\n",
    "which is an additive, zero mean, Gaussian noise model. We can calculate the minus log-likelihood derived from a set of observations that we stack in a vector $\\vec{\\rho}$, where we also stack the kernels into a matrix $\\mathbf{K}$, and the uncertainties into a covariance matrix $\\mathbf{C}_{obs}$ (in this case, diagonal), we have:\n",
    "\n",
    "$$\n",
    "\\log\\left[p(\\rho|\\vec{f}\\right] \\propto \\frac{1}{2}\\left[\\vec{\\rho^{\\lambda_i}} - \\vec{f}\\cdot\\mathbf{K}\\right]^{\\top}\\mathbf{C}_{obs}^{-1}\\left[\\vec{\\rho^{\\lambda_i}} - \\vec{f}\\cdot\\mathbf{K}\\right]\n",
    "$$\n",
    "\n",
    "We can find the maximum likelihood estimate by finding the minimum minus log-likelihood estimate: we equate the derivative to 0 and we have your bog-standard linear least squares solution with uncertainty. The optical weights are given by\n",
    "\n",
    "$$\n",
    "\\mathbf{K}^{\\top}\\mathbf{C}^{1}_{obs}\\mathbf{K}\\cdot \\vec{f}_{opt} = \\mathbf{K}^{\\top}\\mathbf{C}^{1}_{obs}\\cdot\\vec{\\rho}\n",
    "$$\n",
    "\n",
    "And the uncertainty by the radius of curvature around the optimal point, or the inverse Hessian\n",
    "$$\n",
    "\\mathbf{C}_{f}^{-1} = \\mathbf{K}^{\\top}\\mathbf{C}^{1}_{obs}\\mathbf{K}\n",
    "$$\n",
    "\n",
    "Implicitly, the solution is assumed to be multivariate Gaussian with mean $\\vec{\\rho}$ and covariance matrix $\\mathbf{C}$\n",
    "\n",
    "You may want to check [this document](https://doi.org/10.6084/m9.figshare.957577.v4) for a more thorough introduction to the maths behind this.\n",
    "\n",
    "Let's see how this works in practice by using the `fit_period_prior` function:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, so let's use this again, and now plot the kernel weight uncertainties as vertical streak. We are not using any prior information here, although `fit_period_prior` can deal with that. It's the same as before, but now showing the uncertainty I have calculated as detailed above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuYAAADcCAYAAAA852SxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdy0lEQVR4nO3de5xcdXn48c8mBKjgpUICxp8BBXxQEOpP8OWtsohajZSLBMT8FBDw2ggFLwheoqGIrVW8hBSrKKiFWtQAEhUpuFCggKho5PIYBSEJGCkFEX6kYWH6xzkLk2V2dvYye2Z2Pu/Xa1+ZOd9z5jz7sLM8893nfE9frVZDkiRJUrVmVB2AJEmSJAtzSZIkqSNYmEuSJEkdwMJckiRJ6gAW5pIkSVIHsDCXJEmSOsAmVZw0IvqAZcB84C5gQWauqRtfAHwa2ACcl5kfKbdfD2xZ7nZDZh46pYFLkiRJbVJJYQ7sB8wGtgcOAU4BDgeIiE2AzwCvBNYCl0XEnsDPgPszc4/JDiYiapk52S8rSZIk1etrNlhVK8t84OzMrAHLgf66sTnAFZl5e2YOAtcBOwNzKWbXJUmSpGmnqsJ8HrAGIDM3ADMjYkb5/M7MfCtARGwHHAxcC2wH7B4Rv4iIayPi5dWELkmSJE2+qlpZasBg3fPBzHy0foeIOJSipeXDmfnriNgC+ApwOrArcH5E7JiZD490kojoZ+PZeEmSJKkjVVWYr6VoTVkZEbOA9fWDEXESsADYJzNvKTffAqws21t+HhHrgG0oZ94bycwBYGC0YCJi8Ti+B0mSJGnSVNXKsgJYWD5eCFwyNBARs4F3AP11RTnAscDfl/vsADwVuHNKopUkSZLarK9Wq035ScvlEr8MvApYDRwEDC19mMB32LjoPpGieD8HeD7wR+DYzLxykuJxVRZJkiS1W9NVWSopzDuNhbna4ZHzv7DR85kHHFNRJJIkqUM0Lcyr6jFXj7FIlSRJaq6qHnNJkiRJdSzMJUmSpA5gK4ukSWPLkiRJ42dhLkkT4IcRSdJksZVFkiRJ6gDOmFegF2fYeuF7HK4Xv2dJkjR+zphLkiRJHcDCXJIkSeoAtrJImjS270iSNH4W5pI0AX4YkaTu0enX+dnKIkmSJHUAC3NJkiSpA9jKUoFO+7OJJEmSqmdhLkmSpJ7Q6ZOjFuaSJKnndfpFge3Qi99zp7PHXJIkSeoAFuaSJElSB6iklSUi+oBlwHzgLmBBZq6pG18AfBrYAJyXmR8Z7ZhusvTGgY2eL9qlv5I4JEmS1Dmq6jHfD5gNbA8cApwCHA4QEZsAnwFeCawFLouIPYG5Ix0jSZI0Eb3YX92L33Onq6qVZT5wdmbWgOVAf93YHOCKzLw9MweB64CdRzlGkiRJ6mpVzZjPA9YAZOaGiJgZETMy89HMvBN4K0BEbAccDPwzsHCkY0Y6SUT0YwEvSZKkLlBVYV4DBuueDw4vsCPiUIqWlg9n5q8jYtRjhsvMAWBgtGAiYnGLcUuSJEltUVUry1qKnnEiYhawvn4wIk4CPgjsk5lntXKMJEmS1M2qmjFfQdGacnH57yVDAxExG3gHsFtm3t/KMZKkqeNNSaY//xtL1aiqML8A2DcibgVWAwdFxKJyLIGnA9dFxND+JwLnDz9makOWJEmS2qeSwrxcWeXoYZuX1j1+ygiHDj9GkiRJmha886ckSZLUAapqZZEkSR3KnnKpGhbmkqQxsWiTpPawlUWSJEnqABbmkiRJUgcYsZUlIr46htepZeZRkxCPJEmS1JOa9ZgfBhwJ9I3yGn3AlwELc0mSJGmcmhXmZ2Xm11t5kYh4xSTFI0lSR/EumJKmyog95pn5hJv5RMSXW91XkiRJUuvGevHn4W2JQpIkSepxYy3MR+s3lyRJkjQOYy3Mj2xLFJIkSVKPa7Zc4iWZ+Zr6bZn5jVb3lSRpOvBiT0lTpdmqLK+KiI+18Bp9wN6TFI8kSZLUk5oV5ktovaf85EmIRZK6ztIbBzZ6vmiX/krikCR1vxEL88z8xFQGIkmSJPWysV78KUmSJKkNmrWytE1E9AHLgPnAXcCCzFwzbJ9ZwPWZuXvdtuuBLcunN2TmoVMUsiRJktRWo86YR8TsEbZvNYHz7gfMBrYHTgNOGfbaBwMDwDZ122YC92fmzuWXRbkkSZKmjWbLJc4pH66JiGey8YWgWwM/Bf5snOedD5ydmbWIWA78w7DxVcCpwFfqts2lmF2XJFXIC14lqT2azZj/nqIQ3gRYVz4f+voF8P0JnHcesAYgMzcAMyPisVgy84bMvGjYMdsBu0fELyLi2oh4+QTOL0mSJHWUZquyzACIiJWZ+YJJPm8NGKx7PpiZj45yzIMUM+inA7sC50fEjpn58EgHREQ/0D+xUCVJkqT2G/Xiz8x8QUTMpegH32TY2BXjPO9aitaUleVFnutbOOYWYGVmDgI/j4h1FD3oa0Y6IDMHKHrVm4qIxS2cX5LUg2zdkTRVRi3MI2IJcAJwO/BI3VANeP44z7sCWAhcXP57SQvHHEtxwej7ImIH4KnAneM8vyRJktRRWlku8Vhgj8xcOYnnvQDYNyJuBVYDB0XEIoDMXDrCMUuBcyLiN8AfgaNaaH+RJEmSukIrhfm9TPJqKJlZA44etvkJBXlmblv3+AGKZRYlSZKkaafZcokvLh9+A/huRHyWokCvDe2Tmde1NzxJkiSpNzSbMf/WsOenDXteA54zueFIkiRJvanZconPnspAJEmSpF7WyqosHxthaBC4B7gsM1dNalSSJElSj2nl4s/tgEOA71Dc9XMn4HXARRTF+aci4rjMPKtdQU43roErSZKk4VopzJ8HvCYzrxnaEBF7AZ/OzBdHxCuBs4Gz2hOiJEmSNP3NaGGf5wPD1zC/jqJgH3o8ZzKDkiRJknpNKzPm36W4sc+pwB0Ud9/8ADAQEVsAnwP+s30hSpIkSdNfKzPm7wZ+CXwd+C1wGbAZxQ2C5lAU929tV4CSJElSLxh1xjwz/wf4aPnVyNsmNSJJkiSpBzW78+cvM3O3iLiZurt91svM57ctMkldZ+mNAxs9dwUiSZJa12zGfFH577umIhBJkiSpl43YY56ZV5T/Xk6x8spTgOcCPwduLrdLkiRJmgSjXvwZEXsDa4DjgS8CzwRuioi/anNskiRJUs9oZVWWzwNHZubeQC0zbwYOBT7b1sgkSZKkHtJKYf4s4OLy8dBFoFcA27UlIkmSJKkHtVKY/ztwUkT01W17G95USJIkSZo0rdz5850UNxe6D9gsIu4CfkfRziJJkiRpErRSmD+SmftGxByK9pU/ZObtEzlpOfu+DJgP3AUsyMw1w/aZBVyfmbu3eowkSZLUrVppZbkjIq6jWJXlacAfJuG8+wGzge2B04BT6gcj4mBgANim1WMkSZKkbtZKYf50ipsNraO42dAtEXF5RHxsAuedD5ydmTVgOdA/bHwVcOoYj5EkSZK61qitLJn5SET8BHio/HqQYvZ6K2DJOM87j2JtdDJzQ0TMjIgZmfloue0G4IaIaPmYRiKiHwt4SZIkdYFRC/OIuAT4C4pZ7CuB84C/zcz/nsB5a8Bg3fPBZgX2eI/JzAGKlpimImLxaPtIkiRJ7dTKxZ/3A38CZpZfmwKbTfC8a4G5wMryIs/1bTpGkjTJFu3SX3UIkjQtjdpjnpkHZeZzgAXANRStITdFxG0TOO8KYGH5eCFwSZuOkSRJkrpCK60szwJeUn69FNiJokC/dALnvQDYNyJuBVYDB0XEIoDMXNrqMRM4vyRNCmePJXWrpTcObPTc32fVa6WVZRXFXT5/DHwAuDYzB5sf0ly5ssrRwzY/oSDPzG1HOUZd4rTzr9no+XEHvKSiSCRJkjpTK4X5n2fmQ22PRJIkSephrSyXaFEuSVIPscVBqkYrNxiSJEmS1GattLJIGgdnnCRJ0liMWJhHxLLRDs7M90xuONL08fCqzTfesEs1cUiSpO7QbMZ83ZRFIUmSJPW4EQvzzPxE/fOIeCqwDbCqXLpQkqRpzzY0SVNl1Is/I+JZEXEpcBfwS2CPiPhJRDyn7dFJkiRJPaKVVVm+QnGDoacBtcz8CfAd4Mx2BiZJkiT1klYK85cBn8zMDcBQC8vngD3bFpUkSZLUY1pZLvFXwN7AirptuwO3tSUiSZKkKeYSt+oErRTm7wYuioifA7Mi4lzglcBb2hqZJEmS1ENGLcwz84aIeC6wL3AF8Afg2Mz8Q7uD0/Qxa6f1VYegKeAMkyRJ49dKjznAiyiK+HUUfeavi4jD2haVJEmS1GNGnTGPiDOAQ4BfAPXTnjXg622KS5IkSeoprfSYHwK8MDNvb3cwkiRJUq9qpTC/lceXSZQk1Tn/0lUbPT9gn50qikSS1O1aKcwvAC6OiC8C/1U/kJn/1paoJEkdyw8jktQerRTmrwJ+Dxw8bHsNGFdhHhF9wDJgPnAXsCAz19SNvwH4AtAHfCQzzym3Xw9sWe52Q2YeOp7zS5IkSZ2maWEeETOAZZl53iSfdz9gNrA9RQ/7KcDh5Tk3BU4D9gIeAK6PiAsoLjy9PzP3mORYppyzTb3BJSIlSdJYNF0uMTMfBT4VEXMn+bzzgbMzswYsB/rrxl4ErMzMNZl5H3Al8ApgLsXsuqQOdf6lqzb6kiRJrWulleUc4PsRsQy4r35gAj3m84A15WtsiIiZETGj/CDw2FjpTmBb4EFg94gYWrbx+My8qtlJIqKfjYt+SZIkqSO1Upi/ArgXePOw7ePuMS+PHax7PlgW5Y3GasAjFIX5V4DTgV2B8yNix8x8eKSTZOYAMDBaMBGxeCzBS5IkSZNt1MI8M/duw3nXUrSmrIyIWWx846KhsSFzgYuBWyhaXAaBn0fEOmAbNp5dlyRJkhrq9Ov8Wrnz56bACcABFBdsvg44AlicmQ+N87wrgIUUBfdC4JK6sWuBMyNia4oe+D2BdwLHl+d/X0TsADyVos1FkiRJ6npNL/4sfZ6ineV9wNbAHcAOFMsdjtcFwMMRcStwJPCJiFgUEYvKGfEPAlcBVwMnZuYGYCmwU0T8hqKF5qi69hdJkiSpq7XSY/5mYIfMvCciyMwHIuJoijuCvm08Jy1XYzl62OaldeMXAhcOO+YBimUWJUmSpGmnlRnzdRQz5fW2oFhjXJIkSdIkaGXG/OPADyNiKTAjIt4JvAv4TDsDkyRJUvss2qW/6hCm3G1/umfYli67+DMzz42IVRS94JcCLwFOyswftDs4TR+9+OaXpqtOW8VAkqaLEQvziHj9UPGdmdcD19eNbRIRn8zMk6YgRkmSNIWcTJGq0WzG/NsRcUhmrqjfGBG7At+kWLrQwlxST3P2WJI0WZoV5kcA50bEWzLzwojoo1jG8BPAecAxUxCfpC5ikSpND51+ExZpuhqxMM/M8yLi/wPnRMRJwKEU65cfUi5nKEmSJHWNWTutH32nCjVdLrFsYzkQ+BSwKbCLRbkkSZI0+UZdxzwzLwNeCzwb2KftEUmSJE2xh1dtvtGXVIVmq7I8BNTqNm0KfCsi/gfoA2qZ+aQ2xydJkiT1hGYXf+48ZVH0GC+ikSRJ0nDNLv68fSoDkSRJknrZqD3mkiRJktqvWSuLpAnwznnS9OCa3pKmioW5JElSD/JDZ+exMJckSRuxQJOqYY+5JEmS1AEqmTGPiD5gGTAfuAtYkJlr6sbfAHyBYr30j2TmOaMdI0mSJHWzqmbM9wNmA9sDpwGnDA1ExKbltr2A/wssiYgtmh0jSZIkdbuqCvP5wNmZWQOWA/11Yy8CVmbmmsy8D7gSeMUox0iSJI3brJ3Wb/QlVaGqwnwesAYgMzcAMyNixvCx0p3AtqMcI0mSJHW1qlZlqQGDdc8HM/PREcZqwCOjHNNQRPTjzLokSRqF955QJ6iqMF8LzAVWRsQsYH2DsSFzgYtHOaahzBwABkbbLyIWtxq4JEnSdOCymJ2nqlaQFcDC8vFC4JK6sWuBF0bE1hExB9gTuGaUYyRJkqSuVtWM+QXAvhFxK7AaOCgiFgFk5tKI+CBwFTATOC4zN0TEE46pKHZJkiRp0lVSmJcrqxw9bPPSuvELgQtbOEaSJEmaFqqaMZckqSvYhytpqrjcoCRJktQBLMwlSZKkDmBhLkmSJHUAC3NJkiSpA1iYS5IkSR3AwlySJEnqABbmkiRJUgewMJckSZI6gIW5JEmS1AEszCVJkqQOYGEuSZIkdQALc0mSJKkDWJhLkiRJHWCTqgOQJEmSpsKiXfqrDqEpZ8wlSZKkDmBhLkmSJHUAC3NJkiSpA1TSYx4RmwHnAnsANwOHZOYfh+1zFPBRYAPw7sy8NCI2AW4DHix3W5GZ75u6yCVJkqT2qOriz3cBt2XmGyPiBOB4YPHQYETMAT4A7AY8HfgB8DzgmcCPM/OwqQ9ZkiRJap+qWlnmA2eVj78FvHbY+D4Us+H3Z+bvgHUREcA8YPVUBSlJkiRNlapmzOcBa8rHdwLbNhmv3+dZwOsiYn/gbuA9mXnzSCeJiH6gf3JCliRJktqnqsK8BgzWPX6kyXj9Pr8HvgB8g2KW/RzghSOdJDMHgIHRgomIxcWEvCRJktQ2tczsG2mw7YV5RHwUeNOwzTsBc4Es/71j2PhaYNe650P73AM8lJmPAj+MiDMjoi8zaxOJsVmCulVEfDwzP151HN3CfI2N+Rob8zU25mtszNfYmbOxMV9jM5F8tb0wz8yTgZPrt0XE3wILKS74PAJYMeywHwEnRMQSiraWLTPzjog4A7gR+GJEvIziAtIJFeWSJElSJ6iqleVLwL9GxG+BlcCbASLiVOC6zFweEaeXYxuAt5XHLQHOjYj3UvSYHz3lkUuSJEltUElhnpkPAfs32H5i3eMvURTw9eN3Anu1PUBJkiRpinnnT0mSJKkDWJhLkiRJHcDCXJIkSeoAFuaSJElSB7AwlyRJkjqAhfn0NVB1AF1moOoAusxA1QF0mYGqA+gyA1UH0GUGqg6gCw1UHUCXGag6gC4zMN4D+2o1788jSZIkVc0Zc0mSJKkDWJhLkiRJHcDCXJIkSeoAFuaSJElSB9ik6gA0cRGxP/DSzPxQROwGfBV4CnAV8PbMHIyIo4CPAhuAd2fmpdVFPPUiYgZwJrAPcC/wHuAu4NvA1sC3M/P4ct+ezhU0zldmXlWOLQdOzMxbyufmq/HP138D5wJbAL8C3pqZD5ivEfM1C/gisBnwM+CwzNxgvkZ9P+4OXJSZzyqfm6/GP199wDeB9eVuizPzW+ZrxHz9DDgH2B24HViQmfeYr8IIOfsSj9fVM4H7MnPPsebMwryLRUQf8FlgIfC1cvMZwDGZeXVEfBp4S0R8H/gAsBvwdOAHwPMqCLlKB1J879tRfO//CvwWWAxcBFwUEXsDN2KuoEG+IuII4ATgAOBEgIiYg/mCxj9ftwNLMvO7EXEq8O6IOBvzBY3ztRnw+sy8NSK+BewXEVdgvqBxvnaLiJnAP1J8qPH9+LhG+foMxfvxq0M7ma/HNMrXd4BfZeaBEbEY+JuIOAPzNeQJOcvMXYcGI2IJcNt4fsZsZel+l1J8qh2yfWZeXT6+DHgNxSe6FZl5f2b+DlgXETG1YVZuW+DszKxl5k3AbGAPipmmGnAe8FrM1ZBG+bob+AbF7O8Q81VolK/NgQvK8auAnTFfQxrl60NlUT4L2BL4I+ZryBPyVU7MvJ/irzJDzFeh0c/XPGD1sP3MV6FRvg4APl+On0bx1wbz9biR3pNExHOAV2fm1xhHzpwx72JlQXlRRGxN8T99gLUR8XLgP4H9efwX0pq6Q++k+KHKKQy3Upl5+tDjiHg7cB+wocwhFDnZq9ze07mChvlalZmrgdUR8f66XXv+ZwtGzNdryud/BhxLUUCZL0bM1/KI+EuKv2DdBVwNLMJ8NcwXsCOwV2bOj4hPlsP+fDFivrYD9omILwA/pWg9MF80/fk6LiIOAH4NvB3z9ZgRfocN1RMnUfwlC8aRMwvz6edoYBnFjNPVwINADRis26cGPDL1oVUrIrYAPge8FDgc+Ke64aGcmKvSsHztN8Ju5qvUKF8R8UKKmabLgbMoZjjNF43zlZn/ERFbAadT5OohzBfwhHztT5GjY4bt5vux1ODn63UU78GrgFOAJRRFkvmiYb5+DWRm7hIRJ1Hk63eYr8eM8Dv/ycDLKD7IwDjek7ayTD9PzsyXZ+buwJXATcBaYG7dPnOBO6oIrioR8STgP4AHKFpYfgpsVbfLUE56PlfwxHxl5q0j7Gq+aJyviNgHuBD4QGa+JzMfxXwBDd+P68uZTDJzkKIFaB7mC2iYr7spLsq7MCJuAbaOiOswX8CIv7++mZlXlrOa3wR2xXwBI+brXooWT4DlQGC+HtPk/5FvAi6umz0fc86cMZ9+Ph0Rx1NcUf0uij+hrwZOKC9GmAdsmZm99mZ6F3BlZh43tCEibigv+LwCeAvwIYo3TK/nChrkawQ/wnxB43x9FjggM39at818FTbKV0T8F8XFnv9Y5uOvgWsxX0OG/3ytB54xNBgRv8/MF0fEbMwXNH4/XhsRh2bmLyj+4uDP1+Ma5esy4A0UK5e9GvgJ5qveSP+PfBXFX2aGjDlnFubTz/solvB5ErAsM38GEBGnAysplut5W3XhVeZFwCsj4rV12/6a4urzrYCzhgoocwU0yFdm7jx8p8y823wBT8zXLIqe1n+pu85neWaeaL6Axu/HY4BLIqJGMRP1tcx82HwBvh/HqtHP1zsp3o8zKS5gPyoz7zdfQON87Q2cFREnU/RDH56ZfzRfjxnpPflS4L1128b8nuyr1Wqj7SNJkiSpzewxlyRJkjqAhbkkSZLUASzMJUmSpA5gYS5JkiR1AAtzSZIkqQO4XKIkdahy6cChu8TVgN8Ap2bm16fg3M8Hvgeszsz+uu3zKG5ctn9mXlq3/b0Uy7U+LzMfmuRY+oEfA49k5hP+vxURA8AZwGuAI4BvZOYRkxmDJE0FZ8wlqbPtWBajTwf+DviniNh/Cs67J7CuvigHKG+O8XfAZyNiBkBEPA1YDBw7GUV5udb0cJc3KsqHxXYUcNREzy9JVXHGXJK6QGb+ieIGKbsCHwUuKAvjzwH/j+L3+b9T3MBiLnADsFVmPggQEWuAwzLzsvrXjYjDytfbGlhBcXOM3YCvAn0RMTC8OAc+AxwOHAl8pTz+msy8oHzN7YEvAy+hmF1/R3nHRSLiw+U5tgSuAY7IzDURcRZwH7APxS3T/36kXETEc4F/obit+o8obqgmSV3PGXNJ6i7fB14YEbOA1wP9wM4UdxqdTXGHvluA2yhupU1E/AWwOXBF/QtFxEuAfwAOBp4DDAJLM/NyipnnHzUoysnMh4G/AU6OiN3LfY8pX3MGcCHwHWBb4GxgeTm2C7CIomCfA9wN1N/SegHwpswcsSgvfZniQ8jW5WvvOcr+ktQVLMwlqbuso/jdvRVwHbAvcD+wDbC+3A5FwfqG8vEbgO9l5uCw1zoC+OfMvCEz7wU+DLwxIvpGC6Kceb+8/DotM28th/YEZmbmGZn5YGYuA2aUBfxq4C+BNRRF+4a6eAG+mZk3NTtvRDwT2ANYXL7+WRSz8pLU9WxlkaTuMpvigtB7gHnA14BnAjcDf16333Lgu+Xj+cCnGrzWdsDVdc/vpphZf3KLsZwMvJGN2062A3aOiPV122ZSFOJrKVpvXgAkMAu4vW6/+1o45/+h6H3fULftzhbjlaSO5oy5JHWXvwKuLNtJlgBXZ+YOmbkv8Nu6/a4HiIi9KXqxL2nwWndTFLpDdgLuy8z7W4zlQeDRzKwvwv8AXJ+Zmw99AS+maKM5DvgTsH1mvha4tsXz1FsHbBMRm9VtmzeO15GkjuOMuSR1gYjYnKJt5f3AgeXmTYDNyrFXUxTtv42IvsysRcQFwGnAD4cVz0POA06PiPMpZp0XU/SET8Q1wJyI2Jfiw8ACij72Z5fxzgI2j4gXUVy0ekUrrTNDMvN3EXET8JGI+BRwKLDjBGOWpI7gjLkkdbbfRMQgRR/5EuDIzLy4HDuZYu3ue4CFFBdWHkexqgoU7Sy7A//W6IUz83vAFylWNlkDPEqxwsq4lR8ADixf516Ktc0PLFtPPk/RdnMPRT/7e4D9KD5wjMWbKVZvuZviA8nFzXeXpO7QV6vVqo5BktQGEfEMigsjnzHCjHlXKG8w9PFGK8Q02PcIoN8bDEnqRs6YS9I0FBFbUCxjeG43F+WS1EsszCVpeloOvIUmN+rpMnuVLT0jiogzgTOnKB5JmnS2skiSJEkdwBlzSZIkqQNYmEuSJEkdwMJckiRJ6gAW5pIkSVIHsDCXJEmSOoCFuSRJktQB/hfGy60FKNSTRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "doys, qa, vza, sza, raa, rho, kern, n_obs = prepare_data()\n",
    "colors = [\"#FC8D62\", \"#66C2A5\", \"#8DA0CB\"]\n",
    "plt.figure(figsize=(12,3))\n",
    "doy_start = 183\n",
    "doy_end = doy_start + 8\n",
    "while doy_end <= 273:\n",
    "    f, fwd, obs, rmse, r2, the_unc, f_upper, f_lower, sigma_f = fit_period_prior(\n",
    "        doys, qa, rho, doy_start, doy_end, 0.005, kern)\n",
    "    t0 = 0.5*(doy_start + doy_end)\n",
    "    for i in range(3):\n",
    "        plt.vlines(t0, f[i]-1.96*sigma_f[i], f[i]+1.96*sigma_f[i], \n",
    "                   color=colors[i], lw=4, alpha=0.7)\n",
    "    \n",
    "    doy_end = doy_end + 8\n",
    "    doy_start = doy_start + 8\n",
    "_ = plt.xlabel(\"Day of Year [d]\")\n",
    "_ = plt.ylabel(\"Kernel weight [-]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we se that the isotropic kernel weights have very small uncertainty, whereas the geometric and particularly the volumetric have pretty large uncertainties. All this means is that there's a large range of e.g. volumetric kernel weights that allow one to fit the observations (within the given uncertainties) equally well. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above isn't very satisfying: things are moving all over the place. It is also clear that in areas of low dynamics, we could possibly extend the temporal window to get more observations, improve angular sampling (and perhaps downweight the contribution of the odd outlier), and get a more stable result. However, it's hard to figure out how to chooose the optimal window size. \n",
    "\n",
    "\n",
    "## Enter the Reverend...\n",
    "\n",
    "Moving on forward, we might consider that we know stuff about this problem that we haven't used. I mean, all we've done is just fitted the data to some ropey observations. If we can bring in any extra (perhaps vague) information to our problem, we may get \"better results\". A neat way of doing this is to use Bayesian concepts.\n",
    "\n",
    "![](http://jim-stone.staff.shef.ac.uk/BookBayes2012/HTML_BayesRulev5EbookHTMLFiles/ops/images/f0001-01.jpg)\n",
    "\n",
    "It's Bayes that came up with the idea, and Laplace did a lot of the legwork in making it useful. It is now the remain of the latte-sipping, podcast-listening classes. I don't feel like giving you the lowdown in Bayesian stats, but [this book](http://jim-stone.staff.shef.ac.uk/BookBayes2012/HTML_BayesRulev5EbookHTMLFiles/ops/xhtml/ch01BayesJVSone.html) (where I nicked the photos above from) might be useful.\n",
    "\n",
    "Anyhow, back to our thing. A clever way to improve on this might be to consider whether we can use the solution from the previous time step to constrain the solution from the current time step. Clever hey... In general, this is called \"using a prior pdf\", and is well known in the realm of Bayesian statistics (yeah, that weird bit that you never paid attention to because you were farting about with random forests or somesuch). Again, you may want to consult [this fine document](https://doi.org/10.6084/m9.figshare.957577.v4), and I will assume you understand most of it.\n",
    "\n",
    "Basically, to the likelihood model we introduced above, Bayes' Rule states that the *a posteriori* probability of the model parameters ($\\vec{f}$) given the observations ($\\vec{\\rho}$, $p(\\vec{f}|\\vec{\\rho})$ is proportional to the product of the likelihood times the *a priori* pdf. \n",
    "\n",
    "$$\n",
    "p(\\vec{f}|\\vec{\\rho})\\propto p(\\vec{f}|\\vec{\\rho})\\cdot p(\\vec{rho})\n",
    "$$\n",
    "\n",
    "It turns out that the solution to this assuming a multivariate Gaussian pdf for the prior ($\\mathcal{N}(\\vec{f}_{prior}, \\mathbf{C}_{prior}^{-1})$ is quite simple:\n",
    "$$\n",
    "\\left[\\mathbf{K}^{\\top}\\mathbf{C}^{1}_{obs}\\mathbf{K} + \\mathbf{C}_{prior}^{-1}\\right]\\cdot \\vec{f}_{opt} = \\mathbf{K}^{\\top}\\mathbf{C}^{1}_{obs}\\cdot\\vec{\\rho} + \\mathbf{C}^{1}_{prior}\\cdot\\vec{f}_{prior}\n",
    "$$\n",
    "\n",
    "The clever bit here is not to use the uncertainty from the posterior from the previous time step directly, as that is only pertinent to the previous timestep, but to add a bit of extra uncertainty to account for changes. How much uncertainty you add is entirely up to you, but if you add too much you'll be back to the original problem, and if you add too little, you'll over costrain the solution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAEiCAYAAABKo6fiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhlVXWw8Xd309AKDpFBggoI6EIbQSP6aWKwkWgUCYK0BIkKAs6tBGcB0wbnCaemQxQVHHBAZRCiSCAFAgqiooxLEEEaEMHIKG3TcL4/zilzu7hVdav6nrpdVe/veeqpe8/e+9xVq25X16q9zz6lqiokSZIkSbPbnEEHIEmSJEkaPItDSZIkSZLFoSRJkiTJ4lCSJEmShMWhJEmSJAmLQ0mSJEkSFoeSNCuVUoZKKdWIj5tKKctKKRuMMe7YUspJUxnr2qzJ4yen4HX2L6XcNoH+ryil3FBKOat5XpVS9pjka19bSvnXyYzt8fz/XEr5VSnlT6WUc0spT+5H/1LKKW3GLUkzkcWhJM1e3wYe23xsA7wW2Bv42Bhj3gq8uv3QBqOU8p5SysUTGLIPcESP555QgTfCt4DtJ9D/XcBZwMsm+kJd/gDwLOALEz1Pj6/1t8CXgSOBpwI/AU4vpfzVZPuXUuaXUl4GvKCNmCVpJrM4lKTZ666qqq5tPn5dVdXJwKeArjNMpZR5VVXdWlXV7yfzYqWUeWsSbNtKKetMdExVVb+rqup/24hnxOvcVVXVbycw5EHAxVVV3diH115eVdUda3qeUbwF+GZVVUdXVXUF9R8f7gf2mkz/pni8i7qAnPD3U5JmO4tDSVKnu2l+qW5muq4tpbyulHIL8KKRs0rNEr/LSyn3lFJ+WUrZp6PtAeNHvli3Zaqdyxib9i+XUj5USrm1lHJHKWVpKaV09H9LKeX6UsrdzTLDHTvadiql/KSUsqKUcnUp5eCOtvc0y0KPKKXcAZwELAF2KKVUTZ9SSnlvc/57SilXlFJe03GOvywrHSvWUsr+wBeBhzVLPLdpPj+/41zzSyl3lVIeMDPbOetYSlnYjH1eKeXS5mu7oJSy9XD+gC2Aj5VShrqca34p5T9KKTc3SzN/Xkp58fDXAOzXfK+v7fL9mNvk64bm6zurlPI3I/LxkVLK55r2W0sph42MocOzge8PP6mq6j7gXOA5k+z/C+DJwJOANS6MJWm2sTiUJA0XQU8CXgd8r6NpM+AlwK50/FLejFkIfAU4hnqJ30eBY4YLjfHGT8DewCOAhdQzRW8Ant/E8ErgcOCNwP8DrgFOKaWsU0rZCjgVOLaJ7wjgiFLKAR3nfiYQwE7AgdQzp1dQL7UFOIB6Ge1rgadQL688upSy3QRj/RbwNuDO5tzXAhcBe3aM/UdgXerlvr34KLCYumDaEPhAc/xZwA3Ae6mXvY50GLBLE+uOwP8A3yylPLyJ+dvAGc15RjqcOheLgWcAPwLOKaVs3tHnYOrC7JnAZ4D3lVKeMPJEpZSHNnH/ZkTTcuCRk+lfVdXdVVVdWlXVpcC9XeKXJI3B4lCSZq+XN7NOK4CVwC+B66kLhGHzgFdUVfWTqqruGjH+MOCrVVUdWVXV5VVVfRlYCryzx/G9uhF4XfNL/2epC8AnNW2HAh+pquqkpiB4A5DAlsA7mviOqqrqsqqqvgR8nLoIHPZn4JVVVV1cVdXNwG3Ayqqqrm3ab25e+7Sqqq6kLh7vB7aeSKzN134rcH+zjHcV8HXqGbrh/4tfDJxRVdUfeszLm6uqGqqq6gLguOGcVFW1HFgF/G9VVb/rMu5XwGuqqjq7qqrLgWXAXGCLqqpupV6W+afmPH9RSlmP+r3x9qqqTmy+54cBl1IXi8MurKpqSVVVlwHvo35vPYkHekjz+Z4Rx+/qaFuT/pKkCbI4lKTZ6xTqJXhPBrYDHlZV1XOaImnYHSOLhA7bA2ePOHYxsG2P43t1abN8cNidwINLKfOpi7Tzhhuqqrqzqqqdq6q6mnqm76DhArgpgg8HHtdxrt9UVfWn0V64qqpTgXtLKUc2y18vof6/s4wypGuso/T9BrAJ8MxSX++4G/C10WLponPjnLFeZ6SvAZuWUj5TSjkVOLM5PtrXNGwrYAPG/57/Ja4mF/eMEtvtzeeRbfOBP/ahvyRpgrxYW5Jmr9ub2bCxVGO0PYh6VqjTg4HOYmus8aOZP+L5fV17wfrUBc1oywfXBY4Cjh7jfGPGV0pZCiwCPg38kHp3zOvHGDJarA9QVdXyUsoPqZeWPqj5OLnX8RN5rRFOol5KezRwOvA76q9rPA9qPo/3Pe8prqqq7iql3Ak8ekTTo6mX3a5Rf0nSxDlzKEmarCuprzvrtBOrz2iNZyX/V3RQ6lsSbNLLwGb55a3UM4TD4x9VSrml2ZzlCuDRVVVdOfxBvRPrgd3P2NXLgUOrqvpAVVUn0v8/qn6dujh8MXBqVVV39vn8q2ny+0LqZaUfb2ZGR72v5QhXUxd+f/meNxsD/T0T+553+gH19Y/D55tHfQ3lGX3qL0maAGcOJUmTdSRwbCnlMuqlnQuBf6HeWKVXlwGLmmJuOfAhJjbb+Ang3aWU66iv9/sIcE1VVb8upXwMOL/ZafO/qYuYJdSF2GhWABuVUhY018zdQH1d4AXAY5rxfwa2LBO/NccK4EHN7p6/bK47PIF6VnI/4KUTPN9k3AncAfxzKeVGYAHw7qbt8dRF3grg0aWUrauq+vXwwKqq7iilHAMc2SzRvRF4DfBw4D8nGc8y4PullPOb134X9Y65pwCUUjYANgKWN/kas78kac04cyhJmpSqqr5OvfnMu4CfAq8HDqiq6qwJnOaz1IXbL6lnpm4EfjaB8R8Gvgp8iXrXzbto7nlXVdVPqQuu1zTxvZF6xux73U8F1LubAlzYfN6fenObi4B/p94A51jqIvYBO2qOY4h6+eP5wKZNjLdSf/13sfousa1oCqyXUt/64WfUG/jsR72T7BebbidQ76jareB6cxPnl6m/jgXA86qqmtQ1f8175VXUBep51LuRPr+qquGlq4uodyd9dI/9JUlroFTVZC4HkSRJ/VBKOQ24vKqqtw06FknS7GZxKEnSADT3YdyO+h6I21ZVdc2AQ5IkzXJecyhJ0mC8g/pG9IdZGEqS1gbOHEqSJEmS3JBGkiRJkjTLlpVGRJWZgw5DkiRJkgaljNbgzKEkSZIkyeJQkiRJkmRxKEmSJEmi5WsOI6IAy4BdgZuARZm5vKN9EfBRYCVwQmYe3hy/CNig6XZxZu7TZpySJEmSNNu1vSHN7sDGwJbU93J6P7AfQESsA3wc2Am4ATgrIp4G/Ay4IzN3bDk2SZIkSVKj7WWluwLHZWYFnAgs7GjbBDgnM6/LzFXAhcC2wGbUs4ySJEmSpCnSdnG4ObAcIDNXAnMjYk7z/MbMfDlARGwBvAS4ANgC2CEifhERF0TE37UcoyRJkiTNem0vK62AVR3PV2Xm/Z0dImIf6uWlh2XmryJifeAY4ChgO+CkiNgmM+8d7UUiYiGrz0pKkiRJkiag7eLwBuplopdExDxgRWdjRBwKLAJ2ycwrm8NXApc0S01/HhE3A4+kmYHsJjOHgKHxgomIJZP4GiRJkiRpxmu7ODwN2Bc4vfl8xnBDRGwMvBrYPjPv6BhzMPUmNm+JiK2BhwE3thznQCy9bGi154sXLBxIHJIkSZLUdnF4MrBbRFwDXA/sFRGLm7YEHgFcGBHD/d8FLAWOj4irgduBA0cuRZUkSZIk9VerxWGzS+lBIw4v7Xj80FGG7t5ORJIkSZKkbtrerVSSJEmSNA1YHEqSJEmSLA4lSZIkSe1vSCNpQO476dOrPZ+7x5sGFMn0YL4kSdJs58yhJEmSJMniUJIkSZJkcShJkiRJwuJQkiRJkoQb0kiSJsENfCRJmnmcOZQkSZIkWRxKkiRJklxWKs1YLvObGPMlSZJmO4tDTRte4yRJkiS1x2WlkiRJkiSLQ0mSJEmSy0olSZPgsm5JkmYei8MBqn537eoHFgwkDElSy7xmWpI0HbisVJIkSZJkcShJkiRJclmpJEmSJPXNdL6UwJlDSZIkSZIzh5o+ptNfXSRJkmaK6TwTpolx5lCSJEmSZHEoSZIkSXJZqSRJrXMJliRpOrA4HKCy6ZaDDkGSJElSH03nPwi6rFSSJEmS5MyhJEmSpNFN55kwTYwzh5IkSZIki0NJkiRJUsvLSiOiAMuAXYGbgEWZubyjfRHwUWAlcEJmHj7eGEmSJGlNeFN3qbu2Zw53BzYGtgQ+Abx/uCEi1gE+DiwEFgA7RcTTxhojSZIkSWpH28XhrsBxmVkBJ1IXgsM2Ac7JzOsycxVwIbDtOGMkSZIkSS1oe7fSzYHlAJm5MiLmRsSczLw/M28EXg4QEVsALwE+C+w72pjRXiQiFmIRKUmSJEmT1nZxWAGrOp6vGlnkRcQ+1MtLD8vMX0XEuGNGyswhYGi8YCJiSY9xS5IkSdKs0vay0huAzQAiYh6worMxIg4F3g7skpnH9jJGkiRJktR/bReHp1EvE6X5fMZwQ0RsDLwaWJiZV/YyRpIkSZLUjraXlZ4M7BYR1wDXA3tFxOKmLYFHABdGxHD/dwEnjRzTcoySJEmaRbx1hdRdq8Vhs+PoQSMOL+14/NBRho4cI0mSJElqUdvLSiVJkiRJ04DFoSRJkiTJ4lCSJEmSZHEoSZIkSaL93UolSZLUsvtO+vRqz92NU9JkOHMoSZIkSbI4lCRJkiRZHEqSJEmS8JrDgVq8YOGgQ5AkSZIkwJlDSZIkSRLOHEqSJE177k4qqR+cOZQkSZIkWRxKkiRJkiwOJUmSJElYHEqSJEmSsDiUJEmSJDHGbqUR8YUJnKfKzAP7EI8kSZIkaQDGupXFK4ADgDLOOQrwOcDiUJIkrbH7Tvr0as+9TYMkTY2xisNjM/NLvZwkIp7Vp3gkSZIkSQMw6jWHmXnQyGMR8ble+0qSJEmSpo+JbkizXytRSJIkSZIGaqLF4XjXH0qSJEmSpqGJFocHtBKFJEmSJGmgSlVVXRsi4ozMfG4vJ5lI30GKiCozBx2GJEmSJA3KqKtBx9qt9DkR8W89nnznCYckSZIkSVprjFUcHkHv1xi+tw+xSJIkSZIGZNRlpTORy0olSZIkzXKjTgBOdEMaSZIkSdIMZHEoSZIkSRq/OIyIjUc5vmH/w5EkSZIkDcKoG9JExCbNw+UR8ShWX5u6EfBT4EFjnTwiCrAM2BW4CViUmctH9JkHXJSZO3QcuwjYoHl6cWbu09uXI0mSJEmajLF2K/0dUFEXhTePaLsPOKWH8+8ObAxsCewNvB/Yb7gxIl4C/CvwyI5jc4E7MnPHHs4vSZIkSeqDUYvDzJwDEBGXZOaTJnn+XYHjMrOKiBOBj4xovwr4IHBMx7HNqGcZJUmSJElTZKyZQwAy80kRsRn17N86I9rOGWf45sDypu/KiJgbEXMy8/7m2MXAxRHROWYLYIeI+AWwAnhzZp431otExEJg4XhfiyRJkiSpu3GLw4g4AngHcB31ctJhFfDEcYZXwKqO56uGC8Mx3E09k3gUsB1wUkRsk5n3jjYgM4eAoXHOS0QsGa+PJEmSJM1G4xaHwMHAjpl5ySTOfwP1MtFLmo1nVvQw5krgksxcBfw8Im6mviZx+djDJEmSJEmT1ct9Dv/I5K8BPA3Yt3m8L3BGD2MOBj4MEBFbAw8Dbpzk60uSJEmSejDWrSye3jz8MvCdiDiSukishvtk5oXjnP9kYLeIuAa4HtgrIhY3Y5eOMmYpcHxEXA3cDhzYw1JUSZIkSdIaKFVVdW2IiN+MM7bKzK36H1J7IqLKzEGHoUlaetnQas8XL1g4kDgkSZKkaayM1jDWrSwe204skiRJkqS1TS+7lf7bKE2rgD8AZ2XmVX2NSpIkSZI0pXrZkGYL4G3AVsCDgR2ob22xA7ATcGFE7N9WgJIkSZKk9vVyK4snAM/NzB8PH4iIZwMfzcynR8ROwHHAse2EKEmSJElqWy8zh08ERt7j8ELqonH48Sb9DEqSJEmSNLV6mTn8DvWtJT4I/BbYmHqZ6VBErA98EvhReyFKkiRJktrWy8zh64BfAl8Cfg2cBawHHEQ9Y7gO8PK2ApQkSZIktW/cmcPM/DPw7uajm1f2NSJJ0lrP+45KkjTzjFocRsQvM3P7iLgCqLr1ycwnthaZJEmSJGnKjDVzuLj5/NqpCESSJEmSNDijXnOYmec0n8+m3pH0ocDjgZ8DVzTHJUmSJEkzwLgb0kTEzsBy4M3AZ4BHAZdHxD+2HJskSZIkaYr0slvpp4ADMnNnoMrMK4B9gCNbjUySJEmSNGV6KQ4fA5zePB7emOYcYItWIpIkSZIkTbleisP/Bg6NiNJx7JV443tJkiRJmjHGvc8h8BrgS8BtwHoRcRNwLfXSUkmSJEnSDNBLcXhfZu4WEZtQLyX9fWZe13JckiRJkqQp1Muy0t9GxIXUu5U+HPh9uyFJkiRJkqZaL8XhI4DFwM3Aa4ErI+LsiPi3ViOTJEmSJE2ZcYvDzLwP+An1xjRnAGcDOwB7txuaJEmSJGmqjFscRsQZ1EtJ/xPYCjgB2Cozt2s5NkmSJEnSFOllWekdwJ3A3OZjXWC9NoOSJEmSJE2tXpaV7pWZWwGLgB8DC4HLI+I3LccmSZIkSZoi497KIiIeAzyj+Xgm8DjqIvHMdkOTJEmSJE2VXu5zeBXwI+B/gLcBF2TmqlajkiRJkiRNqV6Kw7/KzHtaj0SSJEmSNDC9XHNoYShJkiRJM1wvu5VKkiRJkma4XpaVSpIkaS229LKh1Z4vXrBwIHFImt5GLQ4jYtl4gzPz9f0NR5IkSZI0CGPNHN68piePiAIsA3YFbgIWZebyEX3mARdl5g69jpEkSZIk9deoxWFm/nvn84h4GPBI4KrMrHo8/+7AxsCWwN7A+4H9Os75EuBfm/P2NEaSJEmS1H/jbkgTEY+JiDOpZ/F+CewYET+JiK16OP+uwHFNMXkisHBE+1XAByc4RpIkzWBLLxta7UOSNDV62ZDmGOBHwAuA2zPzJxHxbeDzwM7jjN0cWA6QmSsjYm5EzMnM+5tjFwMXR0TPY7qJiIVYREqSJEnSpPVSHP4tsGdTqA0vJ/0kcHgPYytgVcfzVWMVeZMdk5lDwNB4wUTEkvH6SJIkSdJs1EtxeCn1DOFpHcd2AH7Tw9gbgM2AS5qNZ1a0NEazgNtyS5quvM2AJGk66KU4fB1wakT8HJgXEV8DdgJe1sPY04B9gdObz2e0NEaS1oi/vEuSpNlu3OIwMy+OiMcDuwHnAL8HDs7M3/dw/pOB3SLiGuB6YK+IWNycd2mvY3p4HUmSJEnSGuhl5hDgqU3f4XsfPj8iyMwvjTWo2XH0oBGHH1AUZuam44yRNEHOhEmSJGkixi0OI+Jo6vsN/oLVr/+rgDGLQ0mSJEnS9NDLzOHewFMy87q2g5EkSZIkDcacHvpcQz1LKEmSJEmaoXqZOTwZOD0iPgPc2tmQmd9sJSpJkiRJ0pTqpTh8DvA74CUjjleAxaEkSZIkzQBjFocRMQdYlpknTFE8kiRJkqQBGLM4zMz7I+JDEXFeZt44VUFJWnP3XjV/9QMLBhOHZiZvjSJJ0szTy7LS44H/iohlwG2dDV5zKEmSJEkzQy/F4bOAPwIvHXHcaw4lSZIkaYYYtzjMzJ2nIhBJkiRJ0uCMWxxGxLrAO4A9gI2B5wP7A0sy855Wo5MkSZIkTYk5PfT5FPXS0rcAGwG/BbYGlrUYlyRJkiRpCvVSHL4U2DczhwAy8y7gIOqZREmSNI57r5q/2ockSWujXorDm6lnDDutD9zV/3AkSZIkSYPQS3H4HuD7EfEWYE5EvAb4LvDxNgOTJEmSJE2dXnYr/VpEXAUcAJwJPAM4NDO/13ZwkiRJkgZr6WVDqz1fvGDhQOJQ+0YtDiPiBcMFYGZeBFzU0bZORHwgMw+dghglSZIkSS0ba1nptyLihSMPRsR21IXifq1FJUmSJEmaUmMtK90f+FpEvCwzT4mIArwd+HfgBOBNUxCfJEmS1Fcuk5S6G7U4zMwTIuJPwPERcSiwD/X9DffOzFOmKkBJkiRJUvvG3K00M08D9gQ+BKwLLLAwlCRJkqSZZ9xbWWTmWcDzgMcCu7QekSRJkiRpyo21W+k9QNVxaF3gGxHxZ6AAVWY+uOX4JEmSJElTYKwNabadsigk9d28x60YdAiSJEmaRsbakOa6qQxEkgbp3qvmr35gwWDikKTJcLdNae0xnXfDHfeaQ0mSJEnSzDfWslJprXLSmVet9nyPXR43oEgkaWJc5i1Jmg6cOZQkSZIkWRxKkiRJkiwOJUmSJEm0fM1hRBRgGbArcBOwKDOXd7S/EPg09X0TD8/M45vjFwEbNN0uzsx92oxTkjQxXgMsSdLM0/aGNLsDGwNbAnsD7wf2A4iIdYFPAM8G7gIuioiTgRXAHZm5Y8uxSZIkSZIabReHuwLHZWYVEScCH+loeypwyfBMYkScCzwLuJx6llGSpsxjH7LhoEPQDDad7nG1NjBfkjQYbReHmwPLATJzZUTMjYg5mXl/Z1vjRmBT4G5gh4j4BfUs4psz87yxXiQiFgIL+x++JEmSJM0ObReHFbCq4/mqpjDs1lYB91EXh8cARwHbASdFxDaZee9oL5KZQ8DQeMFExJKJBC9JkiRJs0Xbu5XeAGwGEBHzqGcCH9DW2Az4LXAlsDQz783MnwM3A49sOU5JkiRJmtXaLg5PA/ZtHu8LnNHRdgHwlIjYKCI2AZ4G/Bg4GPgwQERsDTyMesmpJEmSJKklbS8rPRnYLSKuAa4H9oqIxQCZuTQi3g6cB8wFDmmuS1wKHB8RVwO3Awd2LEWVJEmSJLWg1eIwMyvgoBGHl3a0nwKcMmLMXdS3wJC0BtztT5IkSRPR9syhJE0L3sRdkiTNdm1fcyhJkiRJmgacOdS04cyOJEmS1B6LQ0mSJM0q9141f/UDCwYTh7S2cVmpJEmSJMniUJIkSZLkslJJ0iR4DbAkSTOPxaEkSVqrnHTmVas9948RkjQ1XFYqSZIkSbI4lCRJkiRZHEqSJEmSsDiUJEmSJGFxKEmSJEnC3UolSZIkjWHxgoWDDkFTxJlDSZIkSZIzh5IkSdOd94aU1A8Wh5IkSZpV5j1uxaBD0Ax271XzVz+wYDBxTIbLSiVJkiRJzhwOkktAJEmSJK0tLA4lSZI0q7j7ptSdy0olSZIkSc4cSpKktYuXWUjSYFgcSpIkSRqV+2TMHi4rlSRJkiRZHEqSJEmSLA4lSZIkSVgcSpIkSZKwOJQkSZIk4W6lkiRJ0567R0prj0P2eMagQ5g0Zw4lSZIkSe3OHEZEAZYBuwI3AYsyc3lH+wuBTwMFODwzjx9vjCRJkiSp/9qeOdwd2BjYEvgE8P7hhohYtzn2bOBvgCMiYv2xxkiSJEmS2tF2cbgrcFxmVsCJwMKOtqcCl2Tm8sy8DTgXeNY4YyRJkiRJLWh7Q5rNgeUAmbkyIuZGxJzMvL+zrXEjsOk4Y7qKiIVYREqSJEnSpLVdHFbAqo7nqzqKvJFtFXDfOGO6yswhYGi8YCJiyfghS5IkSRrmbrizR9vLSm8ANgOIiHnAim5tjc2A344zRpIkSZLUgraLw9OAfZvH+wJndLRdADwlIjaKiE2ApwE/HmeMJEmSJKkFbS8rPRnYLSKuAa4H9oqIxQCZuTQi3g6cB8wFDmmuMXzAmJZjlCRJkqRZr1RVNegYpkxEVJk56DD+4qQzr1rtueu5JUmSJLWsjNbQ9syhxmAxKEmSJGlt0fY1h5IkSZKkacDiUJIkSZJkcShJkiRJsjiUJEmSJGFxKEmSJEnC4lCSJEmShMWhJEmSJAmLQ0mSJEkSFoeSJEmSJGCdQQcw1SJi0CFIkiRJ0qBUmVm6NZSqqqY6GM0gEfGezHzPoOOYrcz/4Jj7wTH3g2PuB8fcD465HxxzP/VcVipJkiRJsjiUJEmSJFkcSpIkSZKwOJQkSZIkYXEoSZIkScLiUJIkSZKExaEkSZIkCYtDSZIkSRIWh1pzQ4MOYJYbGnQAs9jQoAOYxYYGHcAsNjToAGaxoUEHMIsNDTqAWWxo0AHMNqWqqkHHIEmSJEkaMGcOJUmSJEkWh5IkSZIki0NJkiRJEhaHkiRJkiQsDiVJkiRJwDqDDkDTR0S8CHhmZr4zIrYHvgA8FDgPeFVmroqIA4F3AyuB12XmmYOLeHqLiDnA54FdgD8CrwduAr4FbAR8KzPf3PQ1733WLf+ZeV7TdiLwrsy8snlu/vtolPf+/wJfA9YHLgVenpl3mfv+GiX384DPAOsBPwNekZkrzX1/jfMzZwfg1Mx8TPPc3PfRKO/7AnwFWNF0W5KZ3zD3/TVK7n8GHA/sAFwHLMrMP5j7qWFxqHFFRAGOBPYFvtgcPhp4U2aeHxEfBV4WEf8FvA3YHngE8D3gCQMIeabYkzqPW1Dn8evAr4ElwKnAqRGxM3AZ5r0ND8h/ROwPvAPYA3gXQERsgvnvt27v/euAIzLzOxHxQeB1EXEc5r7fuuV+PeAFmXlNRHwD2D0izsHc91u33G8fEXOBj1EX6f7MaUe33H+c+mfOF4Y7mftWdMv9t4FLM3PPiFgCvCEijsbcTwmXlapXZ1L/FWfYlpl5fvP4LOC51H/1OS0z78jMa4GbIyKmNswZZVPguMysMvNyYGNgR+q/HlfACcDzMO9t6Zb/W4AvU89cDTP//dct9/OBk5v284BtMfdt6Jb7dzaF4TxgA+B2zH0bHpD75o+zb6WeNR9m7vuv2/t+c+D6Ef3Mff91y/0ewKea9k9Qz+Ca+ynizKHG1RQip0bERtS/kAHcEBF/B/wIeBH/94N0ecfQG6n/0ecUhjtjZOZRw48j4lXAbcDK5vsBdX6f3Rw3733WJf9XZeb1wPUR8daOrr7v+2yU3D+3ef4g4GDqX5bNfZ+NkvsTI+LvqVcs3AScDyzG3PdVt9wD2wDPzsxdI+IDTbPv+z4bJfdbALtExKeBn1IvdzT3fTbG+/6QiNgD+BXwKsz9lLE41GQdBCyj/ivy+cDdQAWs6uhTAfdNfWgzR0SsD0AfSSsAAAa4SURBVHwSeCawH/AfHc3D+TXvLRmR/91H6Wb+W9At9xHxFOq/IJ8NHEs9o2Lu+6xb7jPzhxGxIXAUdd7vwdz33Yjcv4g6328a0c2fOS3o8r5/PvXPmfOA9wNHUBck5r7PuuT+V0Bm5oKIOJQ699di7qeEy0o1WQ/JzL/LzB2Ac4HLgRuAzTr6bAb8dhDBzQQR8WDgh8Bd1MtJfwps2NFlOL/mvQUj85+Z14zS1fz3WbfcR8QuwCnA2zLz9Zl5P+a+77r83FnRzJyQmauol/Zujrnvuy65v4V6Q45TIuJKYKOIuBBz33ej/Lz/Smae26zW+QqwHea+70bJ/R+pL50BOBEIzP2UceZQk/XRiHgz9Y5Sr6Ve5nU98I6IOIL6l4cNMtN/uJP3WuDczDxk+EBEXNxsQnMO8DLgndQ/HM17/z0g/6P4Aea/37rl/khgj8z8accxc99/q+U+Im6l3oDmY01u/wm4AHPfhpHv+xXAXw83RsTvMvPpEbEx5r7fuv3MuSAi9snMX1DP4vq+b0e33J8FvJB6d/Z/AH6CuZ8yFoearLdQbz38YGBZZv4MICKOAi6h3mb4lYMLb0Z4KrBTRDyv49g/Ue/ktSFw7PAvyua9FQ/If2ZuO7JTZt5i/vtuZO7nUV//89WO/QdOzMx3mfu+6/Zz503AGRFRUf+F/4uZea+57zt/5gxOt/f9a6h/5syl3oTswMy8w9z3Xbfc7wwcGxHvpb6mcL/MvN3cT41SVdX4vSRJkiRJM5rXHEqSJEmSLA4lSZIkSRaHkiRJkiQsDiVJkiRJWBxKkiRJkvBWFpKkGaa55cJ9zdMKuBr4YGZ+aQpe+4nAd4HrM3Nhx/HNgcuBF2XmmR3H30h9a6AnZOY9fY5lIfA/wH2Z+YD/7yNiCDgaeC6wP/DlzNy/nzFIkqYXZw4lSTPRNk1B9AjgfcB/RMSLpuB1nwbc3FkYAjQ3a34fcGREzAGIiIcDS4CD+1EYNvdjG+nsboXhiNgOBA5c09eXJE1/zhxKkmaszLyT+kbW2wHvBk5uirNPAv9C/f/gf1PfUHkz4GJgw8y8GyAilgOvyMyzOs8bEa9ozrcRcBrwRmB74AtAiYihkQUi8HFgP+AA4Jhm/I8z8+TmnFsCnwOeQT3L+OrM/EXTdljzGhsAPwb2z8zlEXEscBuwC/AV4MOj5SIiHg98FdgO+AHw4F5yKEmaPZw5lCTNBv8FPCUi5gEvABYC2wJbABsD+2XmlcBvgH8AiIgnA/OBczpPFBHPAD4CvATYClgFLM3Ms6ln4H7QpTAkM+8F3gC8NyJ2aPq+qTnnHOAU4NvApsBxwIlN2wJgMXXRuAlwC3BIx6kXAf+cmaMWho3PURfCGzXnfto4/SVJs4zFoSRpNriZ+v+8DYELgd2AO4BHAiua41AXTS9sHr8Q+G5mrhpxrv2Bz2bmxZn5R+Aw4MURUcYLopmBPLv5+ERmXtM0PQ2Ym5lHZ+bdmbkMmNMUkdcDfw8spy4cV3bEC/CVzLx8rNeNiEcBOwJLmvMfSz07KUnSX7isVJI0G2xMvUnNH4DNgS8CjwKuAP6qo9+JwHeax7sCH+pyri2A8zue30I9w/iQHmN5L/BiVl8CugWwbUSs6Dg2l7oYvIF6GeyTgATmAdd19Luth9d8NPW1kCs7jt3YY7ySpFnCmUNJ0mzwj8C5zdLOI4DzM3PrzNwN+HVHv4sAImJn6mvzzuhyrluoi61hjwNuy8w7eozlbuD+zOwsBH8PXJSZ84c/gKdTL2k9BLgT2DIznwdc0OPrdLoZeGRErNdxbPNJnEeSNIM5cyhJmrEiYj71EtK3Ans2h9cB1mva/oG6cPx1RJTMrCLiZOATwPdHFHDDTgCOioiTqGffllBfI7gmfgxsEhG7UReki6iva3xsE+88YH5EPJV6I51zelnGOiwzr42Iy4HDI+JDwD7ANmsYsyRphnHmUJI0E10dEauorys8AjggM09v2t5LfW+/PwD7Um/2cgj1bqNQLy3dAfhmtxNn5neBz1Dv+LkcuJ9659FJa4rQPZvz/JH63od7NstAP0W9BPYP1Nc3vh7YnbronYiXUu9qegt1UXz62N0lSbNNqapq0DFIkrTWiIi/pt6s5a9HmTmcFiJiIfCebjundum7P7AwM/dvNypJ0trMmUNJkhoRsT71LSa+Np0LQ0mSJsPiUJKk/3Mi8DLGuJn8NPPsZnntqCLi88DnpygeSdJazGWlkiRJkiRnDiVJkiRJFoeSJEmSJCwOJUmSJElYHEqSJEmSsDiUJEmSJAH/H5CuyRjnJB0UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAEiCAYAAABKo6fiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZglVXn48e9hGEDBJbKIqICAvugghAj+NDEwhGgUCYIsQaKCgGtGibssZgzuGygOE6KoIIoLKgNCFAmkQUBBVJT1FUWQYRMIu4zDQP3+qGpzp7ndfbu51Zfu/n6ep5++VafOqbdPX5p571mqVFWFJEmSJGl2W2XQAUiSJEmSBs/kUJIkSZJkcihJkiRJMjmUJEmSJGFyKEmSJEnC5FCSJEmShMmhJM0KpZShUko14uumUsriUspaY9Q7rpSyZCpjfTRr+vEzU3Cf/Uopd07g+teWUm4opZzdHFellF0nee9rSyn/Opm6Pbb/T6WUX5dS/lhKOa+U8peP5PqJtidJGp3JoSTNHt8BntF8bQa8CdgL+NQYdd4FvKH90AajlPKBUsolE6iyN3B4j21PKMEb4dvAlhO4/mDgbODVE71Rlw8AXgR8aaLt9HivvwZOAI4Angf8FDijlPIXk7m+h/JnlFKWdfnapI2fT5KmO5NDSZo97q2q6trm67dVVZ0CfBboOsJUSplbVdVtVVX9YTI3K6XMfSTBtq2UsupE61RVdXNVVf/bRjwj7nNvVVW/n0CVxwCXVFV1Yx/uvbSqqrsfaTujeCfwraqqjqmq6krqDx8eAnaf5PXjlW8OXAv85Yiv6/v9g0nSTGByKEmz233AqvDnka5rSylvLqXcCrxi5KhSM4XvilLK/aWUX5VS9u4oe1j9kTfrNk21cxpjU35CKeVjpZTbSil3l1IWlVJKx/XvLKVcX0q5r5lGuE1H2XallJ82o0O/KaUc1FH2gWZa6OGllLuBJcBCYKtSStVcU0opH2zav7+UcmUp5Y0dbfx5WulYsZZS9gO+DDyhmeK5WfP9pR1trVFKubeU8rCR2c5Rx1LK/KbuS0oplzU/24WllE2H+w/YCPhUKWWoS1trlFL+o5RySzP18hellFcO/wzAvs3v+touv485TX/d0Px8Z5dS/mpEf3yilPKFpvy2UsqhI2PosD3wg+GDqqoeBM4D/m6S149X/izgl1VVXTXi64ExYpSkWcvkUJJmoSaBeS7wZuD7HUUbAHsCO9Hxj+6mznzgq8Cx1FP4PgkcO5xojFd/AvYCngTMpx4J+hfgpU0MrwMOA94K/D/gGuDUUsqqzVTB04DjmvgOBw4vpezf0fYLgQC2Aw6gHjm9knqqLcD+1NNo3wRsTT298phSyhYTjPXbwLuBe5q2rwUuBnbrqPsPwGrU03178UlgAXVCtDbwkeb8i4AbgA9ST3sd6VBgxybWbYD/Ab5VSnliE/N3gDObdkY6jLovFgAvAH4MnFtK2bDjmoOAG6n79nPAh0opzx7ZUCnl8U3cvxtRtBR48kSv77G9ZwFPbxLi25rkdtsuP6ckCZNDSZpNXtOMOi0DlgO/op5e966Oa+YCr62q6qdVVd07ov6hwNeqqjqiqqorqqo6AVgEvK/H+r26EXhzVVWXVVX1eeoE8LlN2SHAJ6qqWlJV1WXUyVgCGwPvbeI7uqqqy6uq+grwaeokcNifgNdVVXVJVVW3AHcCy6uqurYpv6W59+lVVV1FnTw+BGw6kVibn/024KFmGu8K4BvUI3TD/+99JXBmVVW399gv76iqaqiqqguB44f7pKqqpcAK4H+rqrq5S71fA2+squqcqqquABYDc4CNqqq6DbgX+GPTzp+VUlanfm+8p6qqk5vf+aHAZdTJ4rCLqqpaWFXV5cCHqN9bz+XhHtd8v3/E+Xs7yiZyfS/tBXUC+T5gZ+rEcaiU8swu95OkWW/C6y0kSdPWqdQblwBUwE1d1pbdPTJJ6LAlcOKIc5cAb+mxfq8ua6YHDrsHeGwpZQ3qJO384YKqqu4BdgAopWwNbN2MLg5bhToBHPa7qqr+ONqNq6o6rZTyj6WUI4BNgGc3bZRRqnSNdZRrv0k9+vfCUsqF1MnKQaNc203nxjlj3WekrwN7llI+Rz2KOZy4jfYzDdsEWAs4p0scm3eLq6qqB0sp948S213N95FlawB3TOL6Xtp7HXDH8AcVpZSLqEeVXw+8p8s9JWlWMzmUpNnjrmY0bCzVGGWPoR4V6vRYoDPZGqv+aNYYcfxg16tgTeqEZrT1YqsBRwPHjNHemPGVUhYBewBHAT+i3v1yrM1LRov1YaqqWlpK+RH11NLHNF+n9Fp/IvcaYQn1CNoxwBnAzdQ/13ge03wf73feU1xVVd1bSrkHeNqIoqdRT7ud0PW9tFdV1Uq/u6qqHiqlXA48pZeYJWm2cVqpJKlXV1GvO+u0HSuPaI1nOf+XdFDqRw6s10vFZvrlbdRrAYfrP7WUcmuzOcuVwNM6Nx6h3on1gO4tdvUa4JCqqj5SVdXJ9P9D1G9QJ4evBE5rRj5b0/Tvy6mnlX66qqrTqEcDe/Eb6sTvz7/zZmOgv2Viv/NOP6Re/zjc3lzqNZRnTvL6UctLKU8q9bM8t+8oX4V6t9IrJhm/JM1ojhxKknp1BHBcM/JyPvUmLP9MvbFKry4H9miSuaXAx5jYaOORwPtLKddRr/f7BHBNVVW/LaV8Crig2Wnzv6mTmIXUidholgHrlFLmNWvmbqBeF3gh8PSm/p+AjcvEH82xDHhMs7vnr5p1hydRj0ruC7xqgu1Nxj3A3cA/lVJuBOYB72/KnkWd5C0DnlZK2bSqqt8OV6yq6u5SyrHAEc061RuBNwJPBP5zkvEsBn5QSrmguffB1DvmngpQSlkLWAdY2vTXmNePVV5V1fJSymXA4uY9cSv1GtV1gc9PMn5JmtEcOZQk9aSqqm9Qb+xxMPAz6rWG+1dVdfYEmvk8deL2K+qRqRuBn0+g/seBrwFfod51816aZ9pVVfUz6oTrjU18b6UeMft+96aAendTgIua7/tRb25zMfDv1BvgHEedxD5sR81xDFFPb7wAWL+J8Tbqn/9eVt4lthVNgvUq6kc7/Jw6OdqXeifZLzeXnUS9FvHULk28o4nzBOqfYx7wkqqquq0R7CWes6nX+72f+gOGtYGXVlU1PHV1D+rdR5/Wy/U9tLc39e/26038z2zi73UTIEmaVUpVTWZ5iCRJmoxSyunAFVVVvXvQsUiS1MnkUJKkKdA8h3EL6mcgbl5V1TUDDkmSpJW45lCSpKnxXuoH0R9qYihJejRy5FCSJEmS5IY0kiRJkqRZNq00IqrMHHQYkiRJkjQoZbQCRw4lSZIkSSaHkiRJkiSTQ0mSJEkSLa85jIgCLAZ2Am4C9sjMpR3lewCfBJYDJ2XmYc35i4G1mssuycy924xTkiRJkma7tjek2QVYF9iY+tlOHwb2BYiIVYFPA9sBNwBnR8S2wM+BuzNzm5ZjkyRJkiQ12p5WuhNwfGZWwMnA/I6y9YBzM/O6zFwBXARsDmxAPcooSZIkSZoibSeHGwJLATJzOTAnIlZpjm/MzNcARMRGwJ7AhcBGwFYR8cuIuDAi/qblGCVJkiRp1mt7WmkFrOg4XpGZD3VeEBF7U08vPTQzfx0RawLHAkcDWwBLImKzzHxgtJtExHxWHpWUJEmSJE1A28nhDdTTRC+NiLnAss7CiDgE2APYMTOvak5fBVzaTDX9RUTcAjyZZgSym8wcAobGCyYiFk7iZ5AkSZKkGa/t5PB0YB/gjOb7mcMFEbEu8AZgy8y8u6POQdSb2LwzIjYFngDc2HKcA7Ho8qGVjhfMmz+QOCRJkiSp7eTwFGDniLgGuB7YPSIWNGUJPAm4KCKGrz8YWAScGBG/Ae4CDhg5FVWSJEmS1F+tJofNLqUHjji9qOP140epuks7EUmSJEmSuml7t1JJkiRJ0jRgcihJkiRJMjmUJEmSJJkcSpIkSZIwOZQkSZIk0f6jLCRpWnhwyVErHc/Z9W0DikSSJGkwHDmUJEmSJDlyKM1UjoRJkiRpIkwOJUkT5ocPkiTNPE4rlSRJkiSZHEqSJEmSnFYqSYDTIiVJkhw5lCRJkiSZHEqSJEmSnFYqzVhOk1SbfH9JkjTzmBxq2nDrfEmSJKk9JocDVN187con5g0kDEmSJEkyOZQkqW3OfJAkTQduSCNJkiRJcuRQkiRJkvplOs8WceRQkiRJkmRyKEmSJElyWqmmkek0JC9JkjRTTOdpkpoYRw4lSZIkSY4cSpLUNj9llyRNByaHA1TW33jQIUiSJEnqo+n8gaDTSiVJkiRJJoeSJEmSJKeVSpIkSRrDdJ4mqYlx5FCSJEmS1O7IYUQUYDGwE3ATsEdmLu0o3wP4JLAcOCkzDxuvjiRJkiSp/9oeOdwFWBfYGDgS+PBwQUSsCnwamA/MA7aLiG3HqiNJkiRJakfbaw53Ao7PzCoiTgY+0VG2HnBuZl4HEBEXAZsDLxqjjiRJkvSIPLjkqJWOXVMn1dpODjcElgJk5vKImBMRq2TmQ5l5I/AagIjYCNgT+Dywz2h1RrtJRMynHoGUJEmSJE1C28lhBazoOF4xMsmLiL2pp5cempm/johx64yUmUPA0HjBRMTCHuOWJEmSpFml7TWHNwAbAETEXGBZZ2FEHAK8B9gxM4/rpY4kSZIkqf/aTg5Pp54mSvP9zOGCiFgXeAMwPzOv6qWOJEmSJKkdbU8rPQXYOSKuAa4Hdo+IBU1ZAk8CLoqI4esPBpaMrNNyjJIkSZI067WaHGZmBRw44vSijtePH6XqyDqSJElSX7g7qdRd29NKJUmSJEnTgMmhJEmSJMnkUJIkSZJkcihJkiRJov3dSiVJktSyB5cctdKxG65ImgxHDiVJkiRJJoeSJEmSJJNDSZIkSRKuORyoBfPmDzoESZIkSQIcOZQkSZIk4cihJEnStOfupJL6wZFDSZIkSZLJoSRJkiTJ5FCSJEmShMmhJEmSJAk3pJEkSY8yDy45aqVjN1uRpKkxanIYEV+aQDtVZh7Qh3gkSZIkSQMw1sjha4H9gTJOGwX4AmByKEmSJEnT1FjJ4XGZ+ZVeGomIF/UpHkmSJEnSAIy6IU1mHjjyXER8oddrJUmSJEnTx0R3K923lSgkSZIkSQM10d1Kx1t/KEmS9Ii4O6kkDcZERw73byUKSZIkSdJAjZocRsSZI89l5gm9XitJkiRJmj7Gmlb6dxHxbz20UYAd+hSPJEmSJGkAxkoOD6f3NYYf7EMskiRJkqQBKVVVDTqGKRMRVWYOOgxJkiRJGpRRBwAnuiGNJEmSJGkGMjmUJEmSJI2fHEbEuqOcX7v/4UiSJEmSBmHUDWkiYr3m5dKIeCorz01dB/gZ8JixGo+IAiwGdgJuAvbIzKUjrpkLXJyZW3WcuxhYqzm8JDP37u3H0Uy26PKhlY4XzJs/kDgkSZKkmWis3UpvBirqpPCWEWUPAqf20P4uwLrAxsBewIeBfYcLI2JP4F+BJ3ecmwPcnZnb9NC+JEmSJKkPRk0OM3MVgIi4NDOfO8n2dwKOz8wqIk4GPjGi/Grgo8CxHec2oB5llCRJkiRNkbFGDgHIzOdGxAbUo3+rjig7d5zqGwJLm2uXR8SciFglMx9qzl0CXBIRnXU2AraKiF8Cy4B3ZOb5Y90kIuYD88f7WSRJkiRJ3Y2bHEbE4cB7geuop5MOq4DnjFO9AlZ0HK8YTgzHcB/1SOLRwBbAkojYLDMfGK1CZg4BQ+O0S0QsHO8aSZIkSZqNxk0OgYOAbTLz0km0fwP1NNFLm41nlvVQ5yrg0sxcAfwiIm6hXpO4dOxqkiRJkqTJ6uU5h3cw+TWApwP7NK/3Ac7soc5BwMcBImJT4AnAjZO8vyRJkiSpB2M9yuL5zcsTgO9GxBHUSWI1fE1mXjRO+6cAO0fENcD1wO4RsaCpu2iUOouAEyPiN8BdwAE9TEWVJEmSJD0CY00r/eaI4yNHHFfAJmM1npkVcOCI0w9LCjNz/Y7X91I/AkOSJEmSNEXGepTFM6YyEEmSJEnS4PSyW+m/jVK0ArgdODszr+5rVJIkSZKkKdXLhjQbAe+mnkL6WGAr6kdbbAVsB1wUEfu1FaAkSZIkqX29PMri2cCLM/MnwyciYnvgk5n5/IjYDjgeOK6dECVJkiRJbetl5PA5wMhnHF5EnTQOv16vn0FJkiRJkqZWLyOH36V+tMRHgd8D61JPMx2KiDWBzwA/bi9ESZIkSVLbehk5fDPwK+ArwG+Bs4HVqR9RsR51gvmatgKUJEmSJLVv3JHDzPwT8P7mq5vX9TUiSZIkSdKUGzU5jIhfZeaWEXEl9QPvHyYzn9NaZJIkSZKkKTPWyOGC5vubpiIQSZIkSdLgjLrmMDPPbb6fQ70j6eOBZwG/AK5szkuSJEmSZoBxN6SJiB2ApcA7gM8BTwWuiIh/aDk2SZIkSdIU6WW30s8C+2fmDkCVmVcCewNHtBqZJEmSJGnK9JIcPh04o3k9vDHNucBGrUQkSZIkSZpy4z7KAvhv4JCIWNhx7nX44HtJmrUWXT600vGCefMHEockSeqfXpLDNwJfAe4EVo+Im4BrqaeWSpIkSZJmgF6Swwczc+eIWI96KukfMvO6luOSJEmSJE2hXtYc/j4iLqLerfSJwB/aDUmSJEmSNNV6SQ6fBCwAbgHeBFwVEedExL+1GpkkSZIkacqMmxxm5oPAT6k3pjkTOAfYCtir3dAkSZIkSVNl3OQwIs6knkr6n8AmwEnAJpm5RcuxSZIkSZKmSC/TSu8G7gHmNF+rAau3GZQkSZIkaWr1Mq1098zcBNgD+AkwH7giIn7XcmySJEmSpCky7qMsIuLpwAuarxcCz6ROEs9qNzRJj4QPKZckSdJE9PKcw6uBHwP/A7wbuDAzV7QaldSFyY0kSZLUnl6Sw7/IzPtbj0SSJEmSNDC9rDk0MZQkSZKkGa6X3UolSZIkSTOcyaEkSZIkafQ1hxGxeLzKmfmW/oYjSZIkSRqEsTakueWRNh4RBVgM7ATcBOyRmUtHXDMXuDgzt+q1jiRJkiSpv0ZNDjPz3zuPI+IJwJOBqzOz6rH9XYB1gY2BvYAPA/t2tLkn8K9Nuz3VkSRJkiT137hrDiPi6RFxFvUo3q+AbSLipxGxSQ/t7wQc3ySTJwPzR5RfDXx0gnUkqe8WXT600pckSdJs08tzDo8Ffgy8DLgrM38aEd8BvgjsME7dDYGlAJm5PCLmRMQqmflQc+4S4JKI6LlONxExH5NISZIkSZq0XpLDvwZ2axK14emknwEO66FuBazoOF4xVpI32TqZOQQMjRdMRCwc7xpJkvpt5Gj0gnnzBxKHZi7fY5L6oZdHWVzGw0cItwJ+10PdG4AN4M8bzyxrqY4kSZIk6RHoJTl8M/CfEfE9YG5EfB34LvC2HuqeDuzTvN4HOLOlOpIkSZKkR2DcaaWZeUlEPAvYGTgX+ANwUGb+oYf2TwF2johrgOuB3SNiQdPuol7r9HAfSZIkSdIj0MuaQ4DnNdcOP/vwpRFBZn5lrErNjqMHjjj9sKQwM9cfp44kSZIkqUXjJocRcQz18wZ/ycrr/ypgzORQkiRJkjQ99DJyuBewdWZe13YwkiRJ7rwpSYPRy4Y011CPEkqSJEmSZqheRg5PAc6IiM8Bt3UWZOa3WolKkiRJkjSlekkO/w64GdhzxPkKMDmUJEmSpBlgzOQwIlYBFmfmSVMUjyRJkiRpAMZcc5iZDwEfi4gNpigeSZIkSdIA9DKt9ETgvyJiMXBnZ4FrDiVJkiRpZuglOXwRcAfwqhHnXXMoSbOUjxaQJGnmGTc5zMwdpiIQSZIkSdLgjJscRsRqwHuBXYF1gZcC+wELM/P+VqOTJEmSJE2JMTekaXyWemrpO4F1gN8DmwKLW4xLkiRJkjSFellz+Cpg08y8PSLIzHsj4kDgGuB17YYnSZIk9deiy4dWOnYdtVTrJTm8hXrE8PaOc2sC97YSkaS+eODqNVY+MW8wcUiSJGl66CU5/ADwg4hYBKwSEW8E3gR8us3AJEmaKfywRpI0HfSyW+nXI+JqYH/gLOAFwCGZ+f22g5MkSZI0WE7DnT1GTQ4j4mXDCWBmXgxc3FG2akR8JDMPmYIYJUmSJEktG2u30m9HxMtHnoyILagTxX1bi0qSJEmSNKXGSg73A74eEbsARESJiPdSJ4aXAlu0H56kyXrG49Ze6UuSJEkay6jTSjPzpIj4I3BiRBwC7E39fMO9MvPUqQpQkiRJktS+sUYOyczTgd2AjwGrAfNMDCVJkiRp5hkzOQTIzLOBlwDPAHZsPSJJkiRJ0pQba7fS+4Gq49RqwDcj4k9AAarMfGzL8UmSJEnStDGdH/0x1nMON5+yKCRpwKbTH25JkqQ2jLUhzXVTGYgkSZIkaXDGGjmUpFljyVlXr3S8647PHFAkkiRJg2FyKElSy3zWqNrm1HhJ/TDubqWSJEmSpJnPkUNphnJapPTo4X+PkqTpwORQkjRhrtGUJGnmaTU5jIgCLAZ2Am4C9sjMpR3lLweOon5u4mGZeWJz/mJgreaySzJz7zbj1PTgP0YlSZKk9rQ9crgLsC6wMbAX8GFgX4CIWA04EtgeuBe4OCJOAZYBd2fmNi3HJkmSJElqtL0hzU7A8ZlZAScD8zvKngdcmplLM/NO4DzgRcAG1KOMkiRJkqQp0vbI4YbAUoDMXB4RcyJilcx8qLOscSOwPnAfsFVE/JJ6FPEdmXn+WDeJiPmsnHhKkiRJkiag7eSwAlZ0HK9oEsNuZRXwIHVyeCxwNLAFsCQiNsvMB0a7SWYOAUPjBRMRCycSvCRJkiTNFm0nhzdQTxO9NCLmUo8EjiwbtgFwBnAV9XTTFcAvIuIW4MmsPMooSZJmKB/oLkmD0faaw9OBfZrX+wBndpRdCGwdEetExHrAtsBPgIOAjwNExKbAE6innEqSJEmSWtL2yOEpwM4RcQ1wPbB7RCwAyMxFEfEe4HxgDvD2Zl3iIuDEiPgNcBdwQMdUVEmSJElSC1pNDptdSg8ccXpRR/mpwKkj6txL/QgMSZIkSdIUaXvkUJKmhV13fOagQ5AkSRqottccSpIkSZKmAZNDSZIkSZLTSiVJ0qPLkrOuXunYad+SNDVMDiVJE+Y/1iVJmnmcVipJkiRJcuRQ04cjFZIkqR8euHqNlU/MG0wc0qONyaEkSZJmlWc8bu1BhyA9KpkcSpIkSRrVgnnzBx2CpohrDiVJkiRJJoeSJEmSJJNDSZIkSRImh5IkSZIkTA4lSZIkSbhbqSRJ0rS35KyrVzr22cCSJsORQ0mSJEmSI4eSJOnRxVEvSRoMk0NJkiTNKn4AIXVncjhArg+QJEmSZpan3fzUlU/MG0wck+GaQ0mSJEmSyaEkSZIkyeRQkiRJkoRrDiVJkiSNwX0yZg9HDiVJkiRJJoeSJEmSJJNDSZIkSRKuOZQkSZr2XAMmqR8cOZQkSZIkOXIoSZIkSf0ynUfyHTmUJEmSJLU7chgRBVgM7ATcBOyRmUs7yl8OHAUU4LDMPHG8OpIkSZKk/mt75HAXYF1gY+BI4MPDBRGxWnNue+CvgMMjYs2x6kiSJEmS2tF2crgTcHxmVsDJwPyOsucBl2bm0sy8EzgPeNE4dSRJkiRJLWh7Q5oNgaUAmbk8IuZExCqZ+VBnWeNGYP1x6nQVEfMxiZQkSZL6bjpvsKKJaTs5rIAVHccrOpK8kWUV8OA4dbrKzCFgaLxgImLh+CFLkiRJ0uzT9rTSG4ANACJiLrCsW1ljA+D349SRJEmSJLWg7eTwdGCf5vU+wJkdZRcCW0fEOhGxHrAt8JNx6kiSJEmSWtD2tNJTgJ0j4hrgemD3iFgAkJmLIuI9wPnAHODtzRrDh9VpOUZJkiRJmvVKVVWDjmHKRESVmYMO48+WnHX1Sscu9pUkSZLUsjJaQdsjhxqDyaAkSZKkR4u21xxKkiRJkqYBk0NJkiRJksmhJEmSJMnkUJIkSZKEyaEkSZIkCZNDSZIkSRImh5IkSZIkTA4lSZIkSZgcSpIkSZKAVQcdwFSLiEGHIEmSJEmDUmVm6VZQqqqa6mA0g0TEBzLzA4OOY7ay/wfHvh8c+35w7PvBse8Hx74fHPt+6jmtVJIkSZJkcihJkiRJMjmUJEmSJGFyKEmSJEnC5FCSJEmShMmhJEmSJAmTQ0mSJEkSJoeSJEmSJEwO9cgNDTqAWW5o0AHMYkODDmAWGxp0ALPY0KADmMWGBh3ALDY06ABmsaFBBzDblKqqBh2DJEmSJGnAHDmUJEmSJJkcSpIkSZJMDiVJkiRJmBxKkiRJkjA5lCRJkiQBqw46AE0fEfEK4IWZ+b6I2BL4EvB44Hzg9Zm5IiIOAN4PLAfenJlnDS7i6S0iVgG+COwI3AG8BbgJ+DawDvDtzHxHc6393mfd+j8zz2/KTgYOzsyrmmP7v49Gee//L/B1YE3gMuA1mXmvfd9fo/T9XOBzwOrAz4HXZuZy+76/xvmbsxVwWmY+vTm27/tolPd9Ab4KLGsuW5iZ37Tv+2uUvv85cCKwFXAdsEdm3m7fTw2TQ40rIgpwBLAP8OXm9DHA2zLzgoj4JPDqiPgv4N3AlsCTgO8Dzx5AyDPFbtT9uBF1P34D+C2wEDgNOC0idgAux35vw8P6PyL2A94L7AocDBAR62H/91u39/51wOGZ+d2I+Cjw5og4Hvu+37r1/erAyzLzmoj4JrBLRJyLfd9v3fp+y4iYA3yKOkn3b047uvX9p6n/5nxp+CL7vhXd+v47wGWZuVtELAT+JSKOwb6fEk4rVa/Oov4UZ9jGmXlB8/ps4MXUn/qcnpl3Z+a1wC0REVMb5oyyPnB8ZlaZeQWwLrAN9afHFXAS8BLs97Z06/9bgROoR66G2f/9163v1wBOacrPBzbHvm9Dt75/X5MYzgXWAu7Cvm/Dw/q++XD2XdSj5sPs+/7r9r7fELh+xHX2ff916/tdgc825UdSj+Da91PEkUONq0lETouIdaj/QQZwQ0T8DfBj4BX83x/SpR1Vb6T+jz6nMNwZIzOPHn4dEa8H7gSWN78PqPt3++a8/d5nXSV2cZEAAAc/SURBVPr/6sy8Hrg+It7Vcanv+z4bpe9f3Bw/BjiI+h/L9n2fjdL3J0fE31LPWLgJuABYgH3fV936HtgM2D4zd4qIjzTFvu/7bJS+3wjYMSKOAn5GPd3Rvu+zMd73b4+IXYFfA6/Hvp8yJoearAOBxdSfIl8A3AdUwIqOayrgwakPbeaIiDWBzwAvBPYF/qOjeLh/7feWjOj/XUa5zP5vQbe+j4itqT9BPgc4jnpExb7vs259n5k/ioi1gaOp+/1+7Pu+G9H3r6Du77eNuMy/OS3o8r5/KfXfmfOBDwOHUyck9n2fden7XwOZmfMi4hDqvr8W+35KOK1Uk/W4zPybzNwKOA+4ArgB2KDjmg2A3w8iuJkgIh4L/Ai4l3o66c+AtTsuGe5f+70FI/s/M68Z5VL7v8+69X1E7AicCrw7M9+SmQ9h3/ddl787y5qREzJzBfXU3g2x7/uuS9/fSr0hx6kRcRWwTkRchH3fd6P8vf9qZp7XzNb5KrAF9n3fjdL3d1AvnQE4GQjs+ynjyKEm65MR8Q7qHaXeRD3N63rgvRFxOPU/HtbKTP/Dnbw3Aedl5tuHT0TEJc0mNOcCrwbeR/3H0X7vv4f1/yh+iP3fb936/ghg18z8Wcc5+77/Vur7iLiNegOaTzV9+4/Ahdj3bRj5vl8GPGW4MCJuzsznR8S62Pf91u1vzoURsXdm/pJ6FNf3fTu69f3ZwMupd2f/e+Cn2PdTxuRQk/VO6q2HHwsszsyfA0TE0cCl1NsMv25w4c0IzwO2i4iXdJz7R+qdvNYGjhv+h7L93oqH9X9mbj7yosy81f7vu5F9P5d6/c/XOvYfODkzD7bv+67b3523AWdGREX9Cf+XM/MB+77v/JszON3e92+k/pszh3oTsgMy8277vu+69f0OwHER8UHqNYX7ZuZd9v3UKFVVjX+VJEmSJGlGc82hJEmSJMnkUJIkSZJkcihJkiRJwuRQkiRJkoTJoSRJkiQJH2UhSZphmkcuPNgcVsBvgI9m5lem4N7PAb4HXJ+Z8zvObwhcAbwiM8/qOP9W6kcDPTsz7+9zLPOB/wEezMyH/f8+IoaAY4AXA/sBJ2Tmfv2MQZI0vThyKEmaiTZrEqInAR8C/iMiXjEF990WuKUzMQRoHtb8IeCIiFgFICKeCCwEDupHYtg8j22kc7olhiNiOwA44JHeX5I0/TlyKEmasTLzHuoHWW8BvB84pUnOPgP8M/X/B/+b+oHKGwCXAGtn5n0AEbEUeG1mnt3ZbkS8tmlvHeB04K3AlsCXgBIRQyMTRODTwL7A/sCxTf2fZOYpTZsbA18AXkA9yviGzPxlU3Zoc4+1gJ8A+2Xm0og4DrgT2BH4KvDx0foiIp4FfA3YAvgh8Nhe+lCSNHs4cihJmg3+C9g6IuYCLwPmA5sDGwHrAvtm5lXA74C/B4iIvwTWAM7tbCgiXgB8AtgT2ARYASzKzHOoR+B+2CUxJDMfAP4F+GBEbNVc+7amzVWAU4HvAOsDxwMnN2XzgAXUSeN6wK3A2zua3gP4p8wcNTFsfIE6EV6naXvbca6XJM0yJoeSpNngFur/560NXATsDNwNPBlY1pyHOml6efP65cD3MnPFiLb2Az6fmZdk5h3AocArI6KMF0QzAnlO83VkZl7TFG0LzMnMYzLzvsxcDKzSJJHXA38LLKVOHJd3xAvw1cy8Yqz7RsRTgW2AhU37x1GPTkqS9GdOK5UkzQbrUm9SczuwIfBl4KnAlcBfdFx3MvDd5vVOwMe6tLURcEHH8a3UI4yP6zGWDwKvZOUpoBsBm0fEso5zc6iTwRuop8E+F0hgLnBdx3V39nDPp1GvhVzece7GHuOVJM0SjhxKkmaDfwDOa6Z2Hg5ckJmbZubOwG87rrsYICJ2oF6bd2aXtm6lTraGPRO4MzPv7jGW+4CHMrMzEfwDcHFmrjH8BTyfekrr24F7gI0z8yXAhT3ep9MtwJMjYvWOcxtOoh1J0gzmyKEkacaKiDWop5C+C9itOb0qsHpT9vfUieNvI6JkZhURpwBHAj8YkcANOwk4OiKWUI++LaReI/hI/ARYLyJ2pk5I96Be1/iMJt65wBoR8TzqjXTO7WUa67DMvDYirgAOi4iPAXsDmz3CmCVJM4wjh5Kkmeg3EbGCel3h4cD+mXlGU/ZB6mf73Q7sQ73Zy9updxuFemrpVsC3ujWcmd8DPke94+dS4CHqnUcnrUlCd2vauYP62Ye7NdNAP0s9BfZ26vWNbwF2oU56J+JV1Lua3kqdFJ8x9uWSpNmmVFU16BgkSXrUiIinUG/W8pRRRg6nhYiYD3yg286pXa7dD5ifmfu1G5Uk6dHMkUNJkhoRsSb1Iya+Pp0TQ0mSJsPkUJKk/3My8GrGeJj8NLN9M712VBHxReCLUxSPJOlRzGmlkiRJkiRHDiVJkiRJJoeSJEmSJEwOJUmSJEmYHEqSJEmSMDmUJEmSJAH/HwmpCjW0dQelAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAEiCAYAAABKo6fiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhlVXWw8XfTNIMQNQpIUAEBs1AQNKKfJgqNqFEkCNISJCoIKA6tRBwZTBsc44BT0yGKCqI4oDIIKhJIgYCCqCiDLBkEaUAEIzSNtE3D+f44p8zt4lbVrep76lJV7+956qm6Z5+977qrblfXqr3PPqWqKiRJkiRJs9sagw5AkiRJkjR4FoeSJEmSJItDSZIkSZLFoSRJkiQJi0NJkiRJEhaHkiRJkiQsDiVpViilDJVSqhEft5ZSFpdS1h+j3/GllFOnMtaHsiaPn5yC59m/lHLnBM5/dSnl5lLKuc3jqpSyxySf+4ZSyr9Opm+P4/9zKeXXpZQ/lVIuKKU8dXXO73W8Usrpbb4uSZoJLA4lafb4FvCE5mMr4PXA3sDHxujzduB17Yc2GKWU95ZSLptAl32Ao3oce0IF3gjfBLabwPmHAecCr5zoE3X5A8BzgC9MdJwen+vvgROBo4GnAz8Bziql/PVkzu9lvFLKOqWUVwIvbuM1SdJMYnEoSbPHsqqqbmg+rquq6jTgU0DXGaZSytyqqu6oqur3k3myUsrc1Qm2baWUNSfap6qq31VV9b9txDPieZZVVfXbCXRZF7isqqpb+vDcS6qqWrq644zibcA3qqo6tqqqX1H/8eEBYK9Jnj9me1M8LqMuICf8/Zak2cbiUJJmt3tofmluZrpuKKW8oZRyO/DSkbNKzRK+q0op95ZSfllK2aej7UH9Rz5Zt2WqncsYm/YTSykfLqXcUUpZWkpZVEopHee/rZRyUynlnmYZ4Q4dbTuWUn5SSlleSrm2lHJIR9t7m2WhR5VSlgKnAguB7UspVXNOKaW8rxn/3lLKr0opB3eM8ZdlpWPFWkrZH/gi8IhmiedWzecXdYy1TillWSnlQTOznbOOpZR5Td8XllKuaF7bxaWULYfzB2wGfKyUMtRlrHVKKf9ZSrmtWXr581LKy4ZfA7Bf872+ocv3Y06Tr5ub13duKeXvRuTjI6WUzzXtd5RSjhgZQ4edgO8PP6iq6n7gAuB5kzx/vPZfAE8FngKsduEsSTOdxaEkzUJNAfMU4A3A9zqaNgFeDuxKxy/dTZ95wJeB46iX8H0UOG640Biv/wTsDTwKmEc9E/Qm4EVNDK8BjgTeDPw/4Hrg9FLKmqWULYAzgOOb+I4CjiqlHNAx9rOBAHYEDqSeOf0V9VJbgAOol9G+Hnga9fLKY0sp204w1m8C7wDubsa+AbgU2LOj7z8Ca1Ev9+3FR4EF1AXRo4EPNsefA9wMvI962etIRwC7NLHuAPwP8I1SyiObmL8FnN2MM9KR1LlYADwL+BFwfill045zDqEuvJ4NfAZ4fynlSSMHKqU8vIn7NyOalgCPmej5vYxXVdU9VVVdUVXVFcB9XV6fJKmDxaEkzR6vamadlgMrgF8CN1EXCMPmAq+uquonVVUtG9H/COArVVUdXVXVVVVVnQgsAt7dY/9e3QK8ofml/rPUBeBTmrbDgY9UVXVq8wv/m4AENgfe1cR3TFVVV1ZV9SXg49RF4LA/A6+pquqyqqpuA+4EVlRVdUPTflvz3GdWVXU1dfH4ALDlRGJtXvsdwAPNMt6VwNeoZ+iG/+99GXB2VVV/6DEvh1ZVNVRV1cXACcM5qapqCbAS+N+qqn7Xpd+vgYOrqjqvqqqrgMXAHGCzqqruoF52+admnL8opaxN/d54Z1VVpzTf8yOAK6iLxWGXVFW1sKqqK4H3U7+3nsKD/VXz+d4Rx5d1tE3k/ImOJ0kah8WhJM0ep1MvsXsqsC3wiKqqntcUScOWjiwSOmwHnDfi2GXA1j3279UVzfLAYXcDDyulrENdpF043FBV1d1VVe1cVdW11DN9Bw0XwE0RfCTwxI6xflNV1Z9Ge+Kqqs4A7iulHN0sf72c+v/KMkqXrrGOcu7XgY2AZ5f6esfdgK+OFksXnRvnjPU8I30V2LiU8plSyhnAOc3x0V7TsC2A9Rn/e/6XuJpc3DtKbHc1n0e2rQP8cRLnT3Q8SdI4vDhbkmaPu5rZsLFUY7StSz0r1OlhQGexNVb/0awz4vH9Xc+C9agLmtGWB64FHAMcO8Z4Y8ZXSlkEzAc+DfyQevfLm8boMlqsD1JV1ZJSyg+pl5au23yc1mv/iTzXCKdSL6U9FjgL+B316xrPus3n8b7nPcVVVdWyUsrdwONGND2OetnthM6f6HiSpPE5cyhJ6tXV1NedddqRVWe0xrOC/ys6KPUtBzbqpWOz/PIO6hnC4f6PLaXc3mzO8ivgcVVVXT38Qb0T64HdR+zqVcDhVVV9sKqqU+j/H1G/Rl0cvgw4o6qqu/s8/iqa/L6Eelnpx5uZ0VHvaznCtdSF31++583GQM9lYt/zTj+gvv5xeLy51NdQnj3J8yc6niRpDM4cSpJ6dTRwfCnlSuqlnfOAf6HeWKVXVwLzm2JuCfBhJjbb+AngPaWUG6mv9/sIcH1VVdeVUj4GXNTstPnf1EXMQupCbDTLgQ1KKds018zdTH1d4MXA45v+fwY2LxO/NcdyYN1md89fNtcdnkw9K7kf8IoJjjcZdwNLgX8updwCbAO8p2n7W+oibznwuFLKllVVXTfcsaqqpaWU44CjmyW6twAHA48E/muS8SwGvl9Kuah57sOod8w9HaCUsj6wAbCkydeY5/fQLkmaAGcOJUk9qarqa9SbzxwG/BR4I3BAVVXnTmCYz1IXbr+knpm6BfjZBPr/B/AV4EvUu24uo7mnXVVVP6UuuA5u4nsz9YzZ97oPBdS7mwJc0nzen3pzm0uBf6feAOd46iL2QTtqjmOIennjRcDGTYx3UL/+Zay6S2wrmgLrFdS3dvgZ9QY++1HvJPvF5rSTqXdU7VZQHdrEeSL169gGeGFVVZO6pq95r7yWukC9kHq30RdVVTW8dHU+9e6jj+vl/B7GkyRNQKmqyVweIkmSJqOUciZwVVVV7xh0LJIkdbI4lCRpCjT3YdyW+h6IW1dVdf2AQ5IkaRVecyhJ0tR4F/WN6I+wMJQkPRQ5cyhJkiRJckMaSZIkSdIsW1YaEVVmDjoMSZIkSRqUMlqDM4eSJEmSJItDSZIkSZLFoSRJkiSJlq85jIgCLAZ2BW4F5mfmko72+cBHgRXAyZl5ZHP8UmD95rTLMnOfNuOUJEmSpNmu7Q1pdgc2BDanvrfTB4D9ACJiTeDjwI7AzcC5EfEM4GfA0szcoeXYJEmSJEmNtpeV7gqckJkVcAowr6NtI+D8zLwxM1cClwBbA5tQzzJKkiRJkqZI28XhpsASgMxcAcyJiDWax7dk5qsAImIz4OXAxcBmwPYR8YuIuDgi/qHlGCVJkiRp1mt7WWkFrOx4vDIzH+g8ISL2oV5eekRm/joi1gOOA44BtgVOjYitMvO+0Z4kIuax6qykJEmSJGkC2i4Ob6ZeJnp5RMwFlnc2RsThwHxgl8y8ujl8NXB5s9T05xFxG/AYmhnIbjJzCBgaL5iIWDiJ1yBJkiRJM17bxeGZwL7AWc3ns4cbImJD4HXAdpm5tKPPIdSb2LwtIrYEHgHc0nKcA7HoyqFVHi/YZt5A4pAkSZKktovD04DdIuJ64CZgr4hY0LQl8CjgkogYPv8wYBFwUkRcC9wFHDhyKaokSZIkqb9aLQ6bXUoPGnF4UcfXDx+l6+7tRCRJkiRJ6qbt3UolSZIkSdOAxaEkSZIkyeJQkiRJkmRxKEmSJEnC4lCSJEmShMWhJEmSJAmLQ0mSJEkSFoeSJEmSJCwOJUmSJElYHEqSJEmSgDUHHYAkPRTcf+qnV3k8Z4+3DCiS6cF8SZI08zhzKEmSJEmyOJQkSZIkWRxKkiRJkvCaQ2nG8powSZIkTYTFoSRh8TxR5kuSpJnH4nCAqt/dsOqBbQYShiRJkiRZHEqS1DaXeUuSpgOLQ00b/nIlSZIktcfiUJqhLJ4lSZI0Ed7KQpIkSZJkcShJkiRJclmpJEmtc5m3JGk6sDiUJEmSpD6ZzpsoWhxq2phO/7AkSZKk6cbicIDKxpsPOgRJkiRJAtyQRpIkSZKEM4eSJEmSxjCdr6EbhOmcH2cOJUmSJEkWh5IkSZKklpeVRkQBFgO7ArcC8zNzSUf7fOCjwArg5Mw8crw+kiRJkqT+a/uaw92BDYHNgb2BDwD7AUTEmsDHgR2Bm4FzI+IZwCaj9ZEkSZI0tabzNXSamLaXle4KnJCZFXAKMK+jbSPg/My8MTNXApcAW4/TR5IkSZLUgrZnDjcFlgBk5oqImBMRa2TmA5l5C/AqgIjYDHg58Flg39H6jPYkETEPi0hJkiRJmrS2i8MKWNnxeOXIIi8i9qFeXnpEZv46IsbtM1JmDgFD4wUTEQt7jFuSJEmSZpW2l5XeTH0NIRExF1je2RgRhwPvBHbJzON76SNJkiRJ6r+2Zw7PpF4melbz+ezhhojYEHgdsF1mLu2ljyRJkrS6vKm71F3bxeFpwG4RcT1wE7BXRCxo2hJ4FHBJRAyffxhw6sg+LccoSZIkSbNeq8Vhs+PoQSMOL+r4+uGjdB3ZR5IkSZLUoravOZQkSZIkTQNtLyuVJEmSHlK8xlDqzplDSZIkSZIzh5IkSdOdu29K6gdnDiVJkiRJzhxKkqSHFmfBJGkwnDmUJEmSJFkcSpIkSZJcVjpQC7aZN+gQJEnSDODSW0n9YHEoSZIeUix0JGkwXFYqSZIkSbI4lCRJkiRZHEqSJEmSsDiUJEmSJGFxKEmSJEnC4lCSJEmSxBi3soiIL0xgnCozD+xDPJIkSZKkARjrPoevBg4AyjhjFOBzgMWhJEmSJE1TYxWHx2fml3oZJCKe06d4JEmSJEkDMOo1h5l50MhjEfG5Xs+VJEmSJE0fE92QZr9WopAkSZIkDdRYy0q7Ge/6Q6k1i64cWuXxgm3mDSQOSZIkaSaa6MzhAa1EIUmSJEkaqFGLw4g4e+SxzDyx13MlSZIkSdPHWMtKnxcR/9bDGAXYuU/xSJIkSZIGYKzi8Ch6v8bwfX2IRZIkSZI0IKMWh5n571MZiCRJkiRpcCa6IY0kSZIkaQayOJQkSZIkjV8cRsSGoxx/dP/DkSRJkiQNwqjXHEbERs2XSyLisay6Oc0GwE+BdccaPCIKsBjYFbgVmJ+ZS0acMxe4NDO37zh2KbB+8/CyzNynt5cjSZIkSZqMsXYr/R1QUReFt41oux84vYfxdwc2BDYH9gY+AOw33BgRLwf+FXhMx7E5wNLM3KGH8SVJkiRJfTDWbqVrAETE5Zn5lEmOvytwQmZWEXEK8JER7dcAHwKO6zi2CfUsoyRJkiRpiow1cwhAZj4lIjahnv1bc0Tb+eN03xRY0py7IiLmRMQamflAc+wy4LKI6OyzGbB9RPwCWA4cmpkXjvUkETEPmDfea5EkSZIkdTducRgRRwHvAm6kXk46rAKePE73CljZ8XjlcGE4hnuoZxKPAbYFTo2IrTLzvtE6ZOYQMDTOuETEwvHOkSRJkqTZaNziEDgE2CEzL5/E+DdTLxO9vNl4ZnkPfa4GLs/MlcDPI+I26msSl4zdTZIkSZI0Wb3c5/CPTP4awDOBfZuv9wXO7qHPIcB/AETElsAjgFsm+fySJEmSpB6MdSuLZzZfngh8OyKOpi4Sq+FzMvOSccY/DdgtIq4HbgL2iogFTd9Fo/RZBJwUEdcCdwEH9rAUVZIkSZK0GsZaVvr1EY8/MeJxBWwx1uCZWQEHjTj8oKIwMzfu+HoZ9S0wJEmSJElTZKxbWTxhKgORJEmSJA1OL7uV/tsoTSuBPwDnZuY1fY1K0mpbdOXQKo8XbDNvIHFIkiRpeuhlQ5rNgHdQLyF9GLA99a0ttgd2BC6JiP3bClCSJEmS1L5ebmXxJOAFmfnj4QMRsRPw0cx8ZkTsCJwAHN9OiJIkSZKktvUyc/hkYOQ9Di+hLhqHv96on0FJkiRJkqZWLzOH36a+tcSHgN8CG1IvMx2KiPWATwI/ai9ESZIkSVLbepk5fAPwS+BLwHXAucDa1Leo2Ii6wHxVWwFKkiRJkto37sxhZv4ZeE/z0c1r+hqRJEmSJGnKjVocRsQvM3O7iPgV9Q3vHyQzn9xaZJIkSZKkKTPWzOGC5vPrpyIQSZIkSdLgjFocZub5zefzImJd4PnAxsDXgXUy8/dTE6Ik6aFm0ZVDqzxesM28gcQhSZL6Z9xrDiNiZ+Cb1JvSPBu4APhhRPxLZp7VcnzSX/jLpyRJktSeXnYr/RRwQGbuDFSZ+StgH+DoViOTJEmSJE2ZXu5z+HhgeIZweGOa84HNWolIkgbAZZKSJGm262Xm8L+BwyOidBx7Dd74XpIkSZJmjF5mDg8GvgTcCawdEbcCN1AvLZUkSZIkzQC9FIf3Z+ZuEbER9VLS32fmjS3HJUmSJEmaQr0sK/1tRFwCHAo8EvAWFpIkSZI0w/RSHD4KWADcBrweuDoizouIf2s1MkmSJEnSlBm3OMzM+4GfUG9MczZwHrA9sHe7oUmSJEmSpsq4xWFEnE29lPS/gC2Ak4EtMnPblmOTJEmSJE2RXpaVLgXuBuY0H2sBa7cZlCRJkiRpavWyrHSvzNwCmA/8GJgHXBURv2k5NkmSJEnSFBn3VhYR8XjgWc3Hs4EnUheJ57QbmiRJkiRpqvRyn8NrgB8B/wO8A7g4M1e2GpUkSZIkaUr1Uhz+dWbe23okkiTNUIuuHFrl8YJt5g0kDkmSxtLLNYcWhpIkSZI0w/UycyhpGnJmQpIkSRPRy60sJEmSJEkz3KgzhxGxeLzOmfnG/oYjSZIkSRqEsZaV3ra6g0dEARYDuwK3AvMzc8mIc+YCl2bm9r32kSRJkiT116jFYWb+e+fjiHgE8Bjgmsysehx/d2BDYHNgb+ADwH4dY74c+Ndm3J76SJIGz2taJUmaeca95jAiHh8R51DP4v0S2CEifhIRW/Qw/q7ACU0xeQowb0T7NcCHJthHkiRJktRnvexWehzwI+DFwF2Z+ZOI+BbweWDncfpuCiwByMwVETEnItbIzAeaY5cBl0VEz326iYh5WERKkiRJ0qT1Uhz+PbBnU6gNLyf9JHBkD30rYGXH45VjFXmT7ZOZQ8DQeMFExMLxzpEkSZKk2aiX4vAK6hnCMzuObQ/8poe+NwObAJc3G88sb6mPJK0Wr6GTJEmzXS/3OXwD8F8R8R1gbkR8Ffg28JYe+p4J7Nt8vS9wdkt9JEmSJEmrYdyZw8y8LCL+FtgNOB/4PXBIZv6+h/FPA3aLiOuBm4C9ImJBM+6iXvv08DySJEmSpNVQqmr8u1JExHOBxzOimMzML7UUVysiosrMQYchSZplFl05tMpjlzFLkgaojNYw7sxhRBxLfb/BX7Dq9X8VMK2KQ0mSJElSd71sSLM38LTMvLHtYCRJkiRJg9HLhjTXU88SSpIkSZJmqF5mDk8DzoqIzwB3dDZk5jdaiUqSJEmSNKV6KQ6fB/wOePmI4xVgcShJkiRJM8CYxWFErAEszsyTpygeSZIkSdIAjHsri4i4DnhuZt4yNSG1x1tZSJIkSZrlJn8rC+Ak4LsRsRi4s7PBaw4lSVK/eV9ISRqMXorD5wB/BF4x4rjXHEqSJD0EWFBL6odxi8PM3HkqApEkSZIkDc64xWFErAW8C9gD2BB4EbA/sDAz7201OkmSJEnSlOhlWemngC2AtwHfBX4LbAksBl7TXmiSJElS/7kMV+qul+LwFcCWmfmHiCAzl0XEQcD1WBxKkiRJM5rF9OyxRg/n3AZsMOLYesCy/ocjSZIkSRqEXorD9wLfj4i3AWtExMHAd4CPtxmYJEmSJGnq9LJb6Vcj4hrgAOAc4FnA4Zn5vbaDkyRJkqTpZDovwx21OIyIFw8XgJl5KXBpR9uaEfHBzDx8CmKUJEmSJLVsrGWl34yIl4w8GBHbUheK+7UWlSRJkiRpSo1VHO4PfDUidgeIiBIR76IuDC8Htm0/PEmSJEnSVBh1WWlmnhwRfwJOiojDgX2o72+4d2aePlUBSpIkSZLaN+ZupZl5JrAn8GFgLWAbC0NJkiRJmnnGvZVFZp4LvBB4ArBL6xFJkiRJkqbcWLuV3gtUHYfWAr4eEX8GClBl5sNajk+SJEmSNAXGus/h1lMWhSRJM9ip51yzyuM9dnnigCKRJGl0Y21Ic+NUBiJJkiRJGpxxrzmUJEmSJM18FoeSJEmSpDGvOZQ0jXmNk6TpasE28wYdwrRjziT1gzOHkiRJkiSLQ0mSJElSy8tKI6IAi4FdgVuB+Zm5pKP9JcCnqe+beGRmntQcvxRYvzntsszcp804JUmSJGm2a/uaw92BDYHNgb2BDwD7AUTEWsAngJ2AZcClEXEasBxYmpk7tBybJEmSZiGv0ZS6a7s43BU4ITOriDgF+EhH29OBy4dnEiPiAuA5wFXUs4ySNGXcwEeSpO4spmePtovDTYElAJm5IiLmRMQamflAZ1vjFmBj4B5g+4j4BfUs4qGZeeFYTxIR84B5/Q9fktSNxfTEmJ+J8f0lSYPRdnFYASs7Hq9sCsNubRVwP3VxeBxwDLAtcGpEbJWZ9432JJk5BAyNF0xELJxI8JIkSZI0W7S9W+nNwCYAETGXeibwQW2NTYDfAlcDizLzvsz8OXAb8JiW45QkSZKkWa3tmcMzgX2Bs5rPZ3e0XQx8PiI2oC5SnwEcDBxKvYnN2yJiS+AR1EtOJUmSJOkhbTpfo9l2cXgasFtEXA/cBOwVEQsAMnNRRLwTuBCYA7y1uS5xEXBSRFwL3AUc2LEUVbOY16BIkiRJ7Wm1OMzMCjhoxOFFHe2nA6eP6LOM+hYYklaDxbMkSZImou1rDiVJkiRJ04DFoSRJkiTJ4lCSJEmS1P6GNJI0LXiNpvTQ4b9HSRoMi0NJ0oT5y7skSTOPxaEkSdI05+2eJPWDxaEkSZJmFYtpqTs3pJEkSZIkWRxKkiRJkiwOJUmSJElYHEqSJEmScEMaTSNeLC5JkiS1x5lDSZIkSZLFoSRJkiTJZaWSJEnTnpdeSOoHi0NJkiTNKhbTUncuK5UkSZIkOXMoSZIkaXSnnnPNKo+deZ25LA4HyH9okiRJ0swynX/Hd1mpJEmSJMniUJIkSZJkcShJkiRJwuJQkiRJkoTFoSRJkiQJi0NJkiRJEt7KQpIkSdIYptOtGLR6LA4lSZIkqU+mczHtslJJkiRJksWhJEmSJMniUJIkSZJEy9ccRkQBFgO7ArcC8zNzSUf7S4BPAwU4MjNPGq+PJEmSJKn/2p453B3YENgc+ATwgeGGiFirObYT8HfAURGx3lh9JEmSJEntaLs43BU4ITMr4BRgXkfb04HLM3NJZt4JXAA8Z5w+kiRJkqQWtH0ri02BJQCZuSIi5kTEGpn5QGdb4xZg43H6dBUR87CIlCRJkqRJa7s4rICVHY9XdhR5I9sq4P5x+nSVmUPA0HjBRMTC8UOWJEmSpNmn7WWlNwObAETEXGB5t7bGJsBvx+kjSZIkSWpB28XhmcC+zdf7Amd3tF0MPC0iNoiIjYBnAD8ep48kSZIkqQVtLys9DdgtIq4HbgL2iogFAJm5KCLeCVwIzAHe2lxj+KA+LccoSZIkSbNeqapq0DFMmYioMnPQYfzFqedcs8rjPXZ54oAikSRJkjRLlNEa2p451BgsBiVJkiQ9VLR9zaEkSZIkaRqwOJQkSZIkWRxKkiRJkiwOJUmSJElYHEqSJEmSsDiUJEmSJGFxKEmSJEnC4lCSJEmShMWhJEmSJAlYc9ABTLWIGHQIkiRJkjQoVWaWbg2lqqqpDkYzSES8NzPfO+g4ZivzPzjmfnDM/eCY+8Ex94Nj7gfH3E89l5VKkiRJkiwOJUmSJEkWh5IkSZIkLA4lSZIkSVgcSpIkSZKwOJQkSZIkYXEoSZIkScLiUJIkSZKExaFW39CgA5jlhgYdwCw2NOgAZrGhQQcwiw0NOoBZbGjQAcxiQ4MOYBYbGnQAs02pqmrQMUiSJEmSBsyZQ0mSJEmSxaEkSZIkyeJQkiRJkoTFoSRJkiQJi0NJkiRJErDmoAPQ9BERLwWenZnvjojtgC8ADwcuBF6bmSsj4kDgPcAK4A2Zec7gIp7eImIN4PPALsAfgTcCtwLfBDYAvpmZhzbnmvc+65b/zLywaTsFOCwzr24em/8+GuW9/7/AV4H1gCuAV2XmMnPfX6Pkfi7wGWBt4GfAqzNzhbnvr3F+5mwPnJGZj28em/s+GuV9X4AvA8ub0xZm5tfNfX+NkvufAScB2wM3AvMz8w/mfmpYHGpcEVGAo4F9gS82h48F3pKZF0XER4FXRsR3gXcA2wGPAr4HPGkAIc8Ue1LncTPqPH4NuA5YCJwBnBEROwNXYt7b8KD8R8T+wLuAPYDDACJiI8x/v3V7798IHJWZ346IDwFviIgTMPf91i33awMvzszrI+LrwO4RcT7mvt+65X67iJgDfIy6SPdnTju65f7j1D9zvjB8krlvRbfcfwu4IjP3jIiFwJsi4ljM/ZRwWal6dQ71X3GGbZ6ZFzVfnwu8gPqvPmdm5tLMvAG4LSJiasOcUTYGTsjMKjOvAjYEdqD+63EFnAy8EPPelm75vx04kXrmapj5779uuV8HOK1pvxDYGnPfhm65f3dTGM4F1gfuwty34UG5b/44+3bqWfNh5r7/ur3vNwVuGnGeue+/brnfA/hU0/4J6hlccz9FnDnUuJpC5IyI2ID6FzKAmyPiH4AfAS/l/36QLunoegv1P/qcwnBnjMw8ZvjriHgtcCewovl+QJ3fnZrj5r3PuuT/msy8CbgpIt7ecarv+z4bJfcvaB6vCxxC/cuyue+zUXJ/SkQ8l3rFwq3ARcACzH1fdcs9sBWwU2buGhEfbJp93/fZKLnfDNglIj4N/JR6uaO577Mx3vdvjYg9gF8Dr8XcTxmLQ03WQcBi6r8iXwTcA1TAyo5zKuD+qQ9t5oiI9YBPAs8G9gP+s6N5OG9l4VwAAAauSURBVL/mvSUj8r/7KKeZ/xZ0y31EPI36L8jnAcdTz6iY+z7rlvvM/GFEPBo4hjrv92Lu+25E7l9Kne+3jDjNnzkt6PK+fxH1z5kLgQ8AR1EXJOa+z7rk/tdAZuY2EXE4de5vwNxPCZeVarL+KjP/ITO3By4ArgJuBjbpOGcT4LeDCG4miIiHAT8EllEvJ/0p8OiOU4bza95bMDL/mXn9KKea/z7rlvuI2AU4HXhHZr4xMx/A3Pddl587y5uZEzJzJfXS3k0x933XJfe3U2/IcXpEXA1sEBGXYO77bpSf91/OzAua1TpfBrbF3PfdKLn/I/WlMwCnAIG5nzLOHGqyPhoRh1LvKPV66mVeNwHvioijqH95WD8z/Yc7ea8HLsjMtw4fiIjLmk1ozgdeCbyb+oejee+/B+V/FD/A/Pdbt9wfDeyRmT/tOGbu+2+V3EfEHdQb0Hysye0/ARdj7tsw8n2/HPib4caI+F1mPjMiNsTc91u3nzkXR8Q+mfkL6llc3/ft6Jb7c4GXUO/O/nzgJ5j7KWNxqMl6G/XWww8DFmfmzwAi4hjgcupthl8zuPBmhKcDO0bECzuO/RP1Tl6PBo4f/kXZvLfiQfnPzK1HnpSZt5v/vhuZ+7nU1/98pWP/gVMy8zBz33fdfu68BTg7Iirqv/B/MTPvM/d958+cwen2vj+Y+mfOHOpNyA7MzKXmvu+65X5n4PiIeB/1NYX7ZeZd5n5qlKqqxj9LkiRJkjSjec2hJEmSJMniUJIkSZJkcShJkiRJwuJQkiRJkoTFoSRJkiQJb2UhSZphmlsu3N88rIBrgQ9l5pem4LmfDHwHuCkz53Uc3xS4CnhpZp7TcfzN1LcGelJm3tvnWOYB/wPcn5kP+v8+IoaAY4EXAPsDJ2bm/v2MQZI0vThzKEmaibZqCqJHAe8H/jMiXjoFz/sM4LbOwhCguVnz+4GjI2INgIh4JLAQOKQfhWFzP7aRzutWGI6I7UDgwNV9fknS9OfMoSRpxsrMu6lvZL0t8B7gtKY4+yTwL9T/D/439Q2VNwEuAx6dmfcARMQS4NWZeW7nuBHx6ma8DYAzgTcD2wFfAEpEDI0sEIGPA/sBBwDHNf1/nJmnNWNuDnwOeBb1LOPrMvMXTdsRzXOsD/wY2D8zl0TE8cCdwC7Al4H/GC0XEfG3wFeAbYEfAA/rJYeSpNnDmUNJ0mzwXeBpETEXeDEwD9ga2AzYENgvM68GfgM8HyAingqsA5zfOVBEPAv4CPByYAtgJbAoM8+jnoH7QZfCkMy8D3gT8L6I2L459y3NmGsApwPfAjYGTgBOadq2ARZQF40bAbcDb+0Yej7wz5k5amHY+Bx1IbxBM/YzxjlfkjTLWBxKkmaD26j/z3s0cAmwG7AUeAywvDkOddH0kubrlwDfycyVI8baH/hsZl6WmX8EjgBeFhFlvCCaGcjzmo9PZOb1TdMzgDmZeWxm3pOZi4E1miLyJuC5wBLqwnFFR7wAX87Mq8Z63oh4LLADsLAZ/3jq2UlJkv7CZaWSpNlgQ+pNav4AbAp8EXgs8CvgrzvOOwX4dvP1rsCHu4y1GXBRx+PbqWcY/6rHWN4HvIxVl4BuBmwdEcs7js2hLgZvpl4G+xQggbnAjR3n3dnDcz6O+lrIFR3HbukxXknSLOHMoSRpNvhH4IJmaedRwEWZuWVm7gZc13HepQARsTP1tXlndxnrdupia9gTgTszc2mPsdwDPJCZnYXg74FLM3Od4Q/gmdRLWt8K3A1snpkvBC7u8Xk63QY8JiLW7ji26STGkSTNYM4cSpJmrIhYh3oJ6duBPZvDawJrN23Ppy4cr4uIkplVRJwGfAL4/ogCbtjJwDERcSr17NtC6msEV8ePgY0iYjfqgnQ+9XWNT2jinQusExFPp95I5/xelrEOy8wbIuIq4MiI+DCwD7DVasYsSZphnDmUJM1E10bESurrCo8CDsjMs5q291Hf2+8PwL7Um728lXq3UaiXlm4PfKPbwJn5HeAz1Dt+LgEeoN55dNKaInTPZpw/Ut/7cM9mGeinqJfA/oH6+sY3ArtTF70T8QrqXU1vpy6Kzxr7dEnSbFOqqhp0DJIkPWRExN9Qb9byN6PMHE4LETEPeG+3nVO7nLs/MC8z9283KknSQ5kzh5IkNSJiPepbTHx1OheGkiRNhsWhJEn/5xTglYxxM/lpZqdmee2oIuLzwOenKB5J0kOYy0olSZIkSc4cSpIkSZIsDiVJkiRJWBxKkiRJkrA4lCRJkiRhcShJkiRJAv4/eWjvbouwpvsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAEiCAYAAABKo6fiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhkVXn48e+ZYRCERCMMIsqioC8Kgkbwp4nCENQoEgRZgkQFQUXNKIori8HgGhdwGSZEUUEUF1QWIYoE0iCgICrKIq8oggwgggLDIOM4cH9/3Numpqnuru6p20V3fz/P00933XPPrbfOFE2/9Z5zbqmqCkmSJEnS7DZn0AFIkiRJkgbP5FCSJEmSZHIoSZIkSTI5lCRJkiRhcihJkiRJwuRQkiRJkoTJoSTNSKWUoVJKNeLr1lLK4lLKumP0O7GUcvpUxvpQ1ozjx6fgeQ4opdw1gfNfWUq5uZRyfvO4KqXsPsnnvqGU8ubJ9O3x+v9cSvlFKeWPpZSLSilPW53zV7ddkjQ6k0NJmrm+ATy++doCeB2wD/DRMfq8DXht+6ENRinlPaWUKybQZV/g6B6vPaEEb4SvA9tM4PzDgPOBl0/0ibp8APAc4HMTvU6Pz/V3wMnAMcAzgB8C55RS/mYy5/eh/fGllOVdvp7QxuuXpOnG5FCSZq5lVVXd0Hz9qqqqM4BPAF0rTKWUeVVV3VFV1e8m82SllHmrE2zbSilrTLRPVVW/rarqD23EM+J5llVV9ZsJdFkbuKKqqlv68NxLqqpaurrXGcVbga9VVXV8VVU/p/7w4QFgz0mev7rtWwI3AE8b8XVTn16vJE1rJoeSNLvcC6wBf6l03VBKeX0p5XbgJSOrSs0UvWtKKfeVUn5WStm3o+1B/Uc+Wbdpqp3TGJv2k0spHyql3FFKWVpKWVRKKR3nv7WUclMp5d5mmuB2HW07lFJ+2FR/fllKOaSj7T3NtNCjSylLgdOBo4BtSylVc04ppby3uf59pZSfl1IO7rjGX6aVjhVrKeUA4PPAI5opnls031/Yca21SinLSikPqsx2Vh1LKQuavi8opVzVvLZLSymbD48fsCnw0VLKUJdrrVVK+c9Sym3N1MqflFJeOvwagP2bf+sbuvx7zG3G6+bm9Z1fSvnbEePx4VLKZ5r2O0opR4yMocOOwHeGH1RVdT9wEfAPkzx/ddufBPy0qqprR3z9eYzXIEmzhsmhJM0CTQLzVOD1wLc7mjYC9gZ2oeOP6qbPAuCLwAnUU/Q+ApwwnGiM138C9gEeBSygrvT8K/DCJoZXAUcCbwT+H3A9cGYpZY1mKuBZwIlNfEcDR5dSDuy49rOBAHYADqKunP6ceqotwIHU02hfBzydenrl8aWUrScY69eBtwP3NNe+Abgc2KOj7z8Ca1JP9+3FR4CF1AnPesAHmuPPAW4G3ks97XWkI4Cdm1i3A/4X+Fop5ZFNzN8Azm2uM9KR1GOxEHgW8H3gwlLKJh3nHALcQj22nwLeV0p58sgLlVL+uon71yOalgCPnuj5q9ve/PwkYOMmYb6jSX63HxmLJM1WJoeSNHO9oqk6LQdWAD+jnj73to5z5gGvrKrqh1VVLRvR/wjgS1VVHVNV1TVVVZ0MLALe1WP/Xt0CvL6qqquqqvo0dQL41KbtcODDVVWdXlXVVdTJWAKbAe9s4juuqqqrq6r6AvAx6iRw2J+AV1VVdUVVVbcBdwErqqq6oWm/rXnus6uqupY6eXwA2HwisTav/Q7ggWYa70rgK9QVuuH/174UOLeqqt/3OC6HVlU1VFXVpcBJw2NSVdUSYCXwh6qqftul3y+Ag6uquqCqqmuAxcBcYNOqqu4AlgF/bK7zF6WUh1G/N95RVdVpzb/5EcBV1MnisMuqqjqqqqqrgfdRv7eeyoP9VfP9vhHHl3W0TeT81W2H+oOC9ajfw7tSJ45DpZQndolHkmadCa+/kCRNG2dSb1wCUAG3dllbtnRkktBhG+CUEceuAN7QY/9eXdVM/xt2D/DwUspa1EnaxcMNVVXdA+wEUEp5OvD0pro4bA51Ajjs11VV/XG0J66q6qxSyj+VUo4BngA8ublGGaVL11hHOfer1NW/Z5dSLqVORg4Z5dxuOjfOGet5RvoysHcp5VPUVczhxG201zTsCcC6wAVd4tiyW1xVVd1fSrlvlNjubr6PbFsLuHMS569uO8CrgDuHP8gopVxGXXV+DfCOLjFJ0qxicihJM9fdTTVsLNUYbWtTV4U6PRzoTLbG6j+atUY8vr/rWbAOdUIz2nqwNYHjgOPHuN6Y8ZVSFgF7AZ8Evke9u+VYm5OMFuuDVFW1pJTyPeqppWs3X2f02n8izzXC6dQVsuOBc4DfUr+u8azdfB/v37ynuKqqWlZKuQd43Iimx1FPu53Q+avb3jzHKv+2VVU9UEq5GnhML69JkmY6p5VKkkZzLfW6s047sGpFazwr+L+kg1LfUmCDXjo20y/voF4LONz/saWU25vNWX4OPK5zYxHqnVgP6n7Frl4BHF5V1QeqqjqN/n9o+hXq5PClwFlN5bM1zfi+mHpa6ceqqjqLuhrYi19SJ35/+TdvNgZ6LhP7N+/0Xer1j8PXm0e9hvLcSZ4/6fZSyqNKfa/PHTva51DvVnrNJF+fJM0oVg4lSaM5BjixqaxcTL0Jy79Qb6zSq6uBvZpkbgnwISZWbTwWeHcp5Ubq9X4fBq6vqupXpZSPApc0O23+D3UScxR1Ijaa5cD6pZStmjVzN1OvC7wU2Ljp/ydgszLxW3MsB9Zudvf8WbPu8FTqquT+wMsmeL3JuAdYCvxzKeUWYCvg3U3bk6iTvOXA40opm1dV9avhjlVVLS2lnAAc06xTvQU4GHgk8F+TjGcx8J1SyiXNcx9GvWPumQCllHWB9YElzXiNef7qtFdVtaKUchWwuHnP3E69hnU+8OlJvj5JmlGsHEqSuqqq6ivUG3ccBvyIeq3hgVVVnT+By3yaOnH7GXVl6hbgxxPo/x/Al4AvUO+6uYzmnnVVVf2IOuE6uInvjdQVs293vxRQ724KcFnz/QDqzW0uB/6degOcE6mT2AftqDmOIerpi5cAGzYx3kH9+pex6i6xrWgSrJdR37rhx9TJz/7UO8l+vjntVOq1iGd2ucShTZwnU7+OrYAXVFXVbY1gL/GcT72e793UHzCsB7ywqqrhqat7Ue8u+rhezl/ddurdXS+jXpd5CfDE5vX1ukmQJM1opaoms1xEkiT1opRyNnBNVVVvH3QskiSNxeRQkqQWNPdh3Jr6HohbVlV1/YBDkiRpTK45lCSpHe+kvhH9ESaGkqTpwMqhJEmSJMkNaSRJkiRJs2xaaURUmTnoMCRJkiRpUMpoDVYOJUmSJEkmh5IkSZIkk0NJkiRJEi2vOYyIAiwGdgFuBfbKzCUd7XsBHwFWAKdm5pHN8cuBdZvTrsjMfduMU5IkSZJmu7Y3pNkNmA9sRn2vp/cD+wNExBrAx4AdgJuB8yNie+DHwNLM3K7l2CRJkiRJjbanle4CnJSZFXAasKCjbQPgwsy8MTNXApcBWwIbUVcZJUmSJElTpO3kcBNgCUBmrgDmRsSc5vEtmfkKgIjYFNgbuBTYFNg2In4aEZdGxN+3HKMkSZIkzXptTyutgJUdj1dm5gOdJ0TEvtTTS4/IzF9ExDrACcBxwNbA6RGxRWb+ebQniYgFrFqVlCRJkiRNQNvJ4c3U00SvjIh5wPLOxog4HNgL2Dkzr20OXwtc2Uw1/UlE3AY8mqYC2U1mDgFD4wUTEUdN4jVIkiRJ0ozXdnJ4NrAfcE7z/dzhhoiYD7wW2CYzl3b0OYR6E5u3RsTmwCOAW1qOcyAWXT20yuOFWy0YSBySJEmS1HZyeAawa0RcD9wE7BkRC5u2BB4FXBYRw+cfBiwCTomIXwJ3AweNnIoqSZIkSeqvVpPDZpfSV484vKjj578epetu7UQkSZIkSeqm7d1KJUmSJEnTgMmhJEmSJMnkUJIkSZJkcihJkiRJwuRQkiRJkoTJoSRJkiQJk0NJkiRJEiaHkiRJkiRMDiVJkiRJmBxKkiRJkjA5lCRJkiRhcihJkiRJAtYYdACS9FBw/+mfXOXx3N3fNKBIpgfHS5KkmcfKoSRJkiTJ5FCSJEmS5LRSacZy2p8kSZImwuRwgKrf3rDqga0GEoYkTZgfNkiSNPOYHEoSJjtql5V8SdJ0YHKoacM/riRJkqT2uCGNJEmSJMnKoTRTWVmVJEnSRJgcSpLUMj+skSRNB04rlSRJkiSZHEqSJEmSTA4lSZIkSbjmcKDKhpsNOoRpxTU7kiRJeqibzrdfs3IoSZIkSTI5lCRJkiSZHEqSJEmScM2hJEmSpDFM5zV0gzCdx6fV5DAiCrAY2AW4FdgrM5d0tO8FfARYAZyamUeO10eSJEmS1H9tTyvdDZgPbAYcC7x/uCEi1gA+BiwAtgJ2iIjtx+ojSZIkSWpH28nhLsBJmVkBp1EngsM2AC7MzBszcyVwGbDlOH0kSZIkSS1oe83hJsASgMxcERFzI2JOZj6QmbcArwCIiE2BvYFPA/uN1me0J4mIBZhESpIkSX03ndfQaWLaTg4rYGXH45Ujk7yI2Jd6eukRmfmLiBi3z0iZOQQMjRdMRBzVY9ySJEmSNKu0Pa30ZmAjgIiYByzvbIyIw4F3ADtn5om99JEkSZIk9V/blcOzqaeJntN8P3e4ISLmA68FtsnMpb30kSRJklaXt2aQums7OTwD2DUirgduAvaMiIVNWwKPAi6LiOHzDwNOH9mn5RglSZIkadZrNTlsdhx99YjDizp+/utRuo7sI0mSZgmrOpI0GG2vOZQkSZIkTQNtTyuVJElSy6y2TozjI3Vn5VCSJEmSZOVQkiQ9tFjVkaTBsHIoSZIkSbJyKEmSNN1ZbZXUD1YOJUmSJEkmh5IkSZIkp5UO1MKtFgw6BEmSJEkCrBxKkiRJkjA5lCRJkiRhcihJkiRJwuRQkiRJkoTJoSRJkiQJk0NJkiRJEiaHkiRJkiTGuM9hRHxuAtepMvOgPsQjSZIkSRqAUZND4JXAgUAZ5xoF+AxgcihJkiRJ09RYyeGJmfmFXi4SEc/pUzySJEmSpAEoVVX1fHJEfCYzX9NiPK2KiCozBx2GJmnR1UOrPF641YKBxCFJkiRNY6PODJ3ohjT7r2YgkiRJkqSHoIkmh+OtP5QkSZIkTUMTTQ4PbCUKSZIkSdJAjZocRsS5I49l5sm9nitJkiRJmj7G2q30HyLi33q4RgF26lM8kiRJkqQBGCs5PJre1xi+tw+xSJIkSZIGZNTkMDP/fSoDkSRJkiQNzkQ3pJEkSZIkzUAmh5IkSZKk8ZPDiJg/yvH1+h+OJEmSJGkQRl1zGBEbND8uiYjHsurmNOsDPwLWHuviEVGAxcAuwK3AXpm5ZMQ584DLM3PbjmOXA+s2D6/IzH17ezmSJEmSpMkYa7fS3wIVdVJ424i2+4Eze7j+bsB8YDNgH+D9wP7DjRGxN/Bm4NEdx+YCSzNzux6uL0mSJEnqg7F2K50DEBFXZuZTJ3n9XYCTMrOKiNOAD49ovw74IHBCx7GNqKuMkiRJkqQpMlblEIDMfGpEbERd/VtjRNuF43TfBFjSnLsiIuZGxJzMfKA5dgVwRUR09tkU2DYifgosBw7NzIvHepKIWAAsGO+1SJIkSZK6Gzc5jIijgXcCN1JPJx1WAU8Zp3sFrOx4vHI4MRzDvdSVxOOArYHTI2KLzPzzaB0ycwgYGue6RMRR450jSZIkSbPRuMkhcAiwXWZeOYnr30w9TfTKZuOZ5T30uRa4MjNXAj+JiNuo1yQuGbubJEmSJGmyernP4Z1Mfg3g2cB+zc/7Aef20OcQ4D8AImJz4BHALZN8fkmSJElSD8a6lcUzmx9PBr4ZEcdQJ4nV8DmZedk41z8D2DUirgduAvaMiIVN30Wj9FkEnBIRvwTuBg7qYSqqpBEWXT20yuOFWy0YSBySJEmaHsaaVvrVEY+PHfG4Ap4w1sUzswJePeLwg5LCzNyw4+dl1LfAkCRJkiRNkbFuZfH4qQxEkiRJkjQ4vexW+m+jNK0Efg+cn5nX9TUqSZIkSdKU6mVDmk2Bt1NPIX04sC31rS22BXYALouIA9oKUJIkSZLUvl5uZfFk4PmZ+YPhAxGxI/CRzHxmROwAnASc2E6IkiRJkqS29VI5fAow8h6Hl1EnjcM/b9DPoCRJkiRJU6uXyuE3qW8t8UHgN8B86mmmQxGxDvBx4PvthShJkiRJalsvyeHrgSOBLwAbA38Ezqe+RcUGzTVe0VaAkqSHHu+jKUnSzDNucpiZfwLe3Xx186q+RiRJkiRJmnKjJocR8bPM3CYifk59w/sHycyntBaZJE0hK2GSJGm2G6tyuLD5/rqpCESSJEmSNDijJoeZeWHz/YKIWBt4HrAh8FVgrcz83dSEKNWs5EiSJEntGfdWFhGxE7AEOBT4FPBY4JqI+MeWY5MkSZIkTZFe7nP4CeDAzNwJqDLz58C+wDGtRiZJkiRJmjK9JIcbA+c0Pw9vTHMhsGkrEUmSJEmSplwvyeH/AIdHROk49iq88b0kSZIkzRjj3ucQOBj4AnAX8LCIuBW4gXpqqSRJkiRpBuglObw/M3eNiA2op5L+LjNvbDkuSZIkSdIU6mVa6W8i4jLq3UofCXgLC0mSJEmaYXpJDh8FLARuA14HXBsRF0TEv7UamSRJkiRpyoybHGbm/cAPqTemORe4ANgW2Kfd0CRJkiRJU2XcNYcRcS7wNOA64CLgVODNmfmHlmOTJGlGWHT10CqPF261YCBxSJI0ll6mlS4F7gHmNl9rAg9rMyhJkiRJ0tTqZVrpnpn5BGAv4AfAAuCaiPh1y7FJkiRJkqZIL9NKNwae1Xw9G3gidZJ4XruhSVodTluTJEnSRPRyn8PrgO8D/wu8Hbg0M1e2GpUkSZIkaUr1khz+TWbe13okkqRpw8q0JEkzTy9rDk0MJUmSJGmG62W3UkmSJEnSDNfLtFJJmvGcJilJkma7UZPDiFg8XufMfEN/w5EkSZIkDcJYlcPbVvfiEVGAxcAuwK3AXpm5ZMQ584DLM3PbXvtIkiRJkvpr1OQwM/+983FEPAJ4NHBdZlY9Xn83YD6wGbAP8H5g/45r7g28ubluT30kSZIkSf037oY0EbFxRJxHXcX7GbBdRPwwIp7Qw/V3AU5qksnTgAUj2q8DPjjBPpIkSZKkPutlQ5oTgO8DLwLuzswfRsQ3gM8CO43TdxNgCUBmroiIuRExJzMfaI5dAVwRET336SYiFmASKUmSJEmT1kty+HfAHk2iNjyd9OPAkT30rYCVHY9XjpXkTbZPZg4BQ+MFExFHjXeOJEmSJM1Gvdzn8CoeXCHcFvh1D31vBjaCv2w8s7ylPpIkSZKk1dBLcvh64L8i4lvAvIj4MvBN4E099D0b2K/5eT/g3Jb6SJIkSZJWw7jTSjPzioh4ErArcCHwO+CQzPxdD9c/A9g1Iq4HbgL2jIiFzXUX9dqnh+eRJEmSJK2GUlXj35UiIp4LbMyIZDIzv9BSXK2IiCozBx2GJEmSJA1KGa1h3MphRBxPfb/Bn7Lq+r8KmFbJoSRJkiSpu152K90HeHpm3th2MJIkSZKkwehlQ5rrqauEkiRJkqQZqpfK4RnAORHxKeCOzobM/ForUUmSJEmSplQvyeE/AL8F9h5xvAJMDiVJUl8tunpolccLt1owkDgkabYZMzmMiDnA4sw8dYrikSRJkiQNwJhrDjPzAeBDEbHRFMUjSZIkSRqAXqaVngL8d0QsBu7qbHDNoSRJkiTNDL0kh88B7gReNuK4aw4lSZIkaYYYNznMzJ2mIhBJkiRNjpv4SOqHcZPDiFgTeCewOzAfeCFwAHBUZt7XanSSJEmSpCkx5oY0jU9QTy19K7A+8Btgc2Bxi3FJkiRJkqZQL2sOXwZsnpm/jwgyc1lEvBq4HnhVu+FJkiRJ/eU03IlxvGaPXiqHt1FXDDutAyzrfziSJEmSpEHopXL4HuA7EbEImBMRBwOvAz7WZmCSJEmSNN1M50rruJXDzPwysDf1OsPzgGcBh2fmx1uOTZIkSZI0RUatHEbEizLz2wCZeTlweUfbGhHxgcw8fApilCRJkiS1bKzK4dcj4sUjD0bE1tSJ4v6tRSVJkiRJmlJjJYcHAF+OiN0AIqJExDupE8Mrga3bD0+SJEmSNBVGnVaamadGxB+BUyLicGBf6nWH+2TmmVMVoCRJkiSpfWNuSJOZZwN7AB8C1gS2MjGUJEmSpJln3FtZZOb5EfEC4AxgZ+DrrUclSZJmrem07bskzSRj7VZ6H1B1HFoT+GpE/AkoQJWZD285PkmSJEnSFBircrjllEUhSZIkSRqosTakuXEqA5EkSZIkDc6YG9JIkiRJkmaHcTekkSRJ0kObm/hI6geTQ0mSWnb6edet8nj3nZ84oEgkSRqdyaEkSZJmFSutE+N4zR6uOZQkSZIktVs5jIgCLAZ2AW4F9srMJR3tLwY+SX3fxCMz85Tm+OXAus1pV2Tmvm3GKc1ETmOTJEnSRLQ9rXQ3YD6wGbAP8H5gf4CIWBM4FtgRWAZcHhFnAMuBpZm5XcuxSZIkSVJfTedpuG0nh7sAJ2VmFRGnAR/uaHsGcOVwJTEiLgKeA1xDXWWUJEmzkDMfJGkw2k4ONwGWAGTmioiYGxFzMvOBzrbGLcCGwL3AthHxU+oq4qGZefFYTxIRC4AF/Q9fkiRJkmaHtpPDCljZ8Xhlkxh2a6uA+6mTwxOA44CtgdMjYovM/PNoT5KZQ8DQeMFExFETCV7S7GGlQpIkzXZtJ4c3AxsBV0bEPOpK4Mi2YRsB5wDXUk83XQn8JCJuAx7NqlVGSdIAmUxLkjTztJ0cng3sR5307Qec29F2KfDZiFif+pYa2wMHA4dSb2Lz1ojYHHgE9ZRTSZKmJZNnSdJ00HZyeAawa0RcD9wE7BkRCwEyc1FEvAO4GJgLvKVZl7gIOCUifgncDRzUMRVVkiRJktSCVpPDzKyAV484vKij/UzgzBF9llHfAkNahdPYJEmSpPbMGXQAkiRJkqTBa3taqSRJ0oQ4M0SSBsPkUJqh/ONKkmYPl15I6genlUqSJEmSrBxKEvgpuyTNJlZape5MDiVJE+YfUpIkzTxOK5UkSZIkmRxKkiRJkkwOJUmSJEm45lCSJGnacx2wpH6wcihJkiRJsnIoSZKk2cVKq9SdyaGmDX+RS5IkSe0xOZQkSZI0qtPPu26Vx35gP3O55lCSJEmSZHIoSZIkSXJaqSRJkiT1zXSehmvlUJIkSZJk5XCQpvOnCpIkSZJmFiuHkiRJkiQrh5IkSZJG5+y22cPKoSRJkiTJyqEkSZIk9ct0rrRaOZQkSZIkmRxKkiRJkkwOJUmSJEmYHEqSJEmSMDmUJEmSJGFyKEmSJEnC5FCSJEmSRMv3OYyIAiwGdgFuBfbKzCUd7S8GPgkU4MjMPGW8PpIkSZKk/mu7crgbMB/YDDgWeP9wQ0Ss2RzbEfhb4OiIWGesPpIkSZKkdrSdHO4CnJSZFXAasKCj7RnAlZm5JDPvAi4CnjNOH0mSJElSC1qdVgpsAiwByMwVETE3IuZk5gOdbY1bgA3H6dNVRCzAJFKSJEmSJq3t5LACVnY8XtmR5I1sq4D7x+nTVWYOAUPjBRMRR40fsiRJkiTNPm1PK70Z2AggIuYBy7u1NTYCfjNOH0mSJElSC9pODs8G9mt+3g84t6PtUuDpEbF+RGwAbA/8YJw+kiRJkqQWtD2t9Axg14i4HrgJ2DMiFgJk5qKIeAdwMTAXeEuzxvBBfVqOUZIkSZJmvVJV1aBjmDIRUWXmoMP4i9PPu26Vx7vv/MQBRSJJkiRpliijNbRdOdQYTAYlSZIkPVS0veZQkiRJkjQNmBxKkiRJkkwOJUmSJEkmh5IkSZIkTA4lSZIkSZgcSpIkSZIwOZQkSZIkYXIoSZIkScLkUJIkSZIErDHoAKZaRAw6BEmSJEkalCozS7eGUlXVVAejGSQi3pOZ7xl0HLOV4z84jv3gOPaD49gPjmM/OI794Dj2U89ppZIkSZIkk0NJkiRJksmhJEmSJAmTQ0mSJEkSJoeSJEmSJEwOJUmSJEmYHEqSJEmSMDmUJEmSJGFyqNU3NOgAZrmhQQcwiw0NOoBZbGjQAcxiQ4MOYBYbGnQAs9jQoAOYxYYGHcBsU6qqGnQMkiRJkqQBs3IoSZIkSTI5lCRJkiSZHEqSJEmSMDmUJEmSJGFyKEmSJEkC1hh0AJo+IuIlwLMz810RsQ3wOeCvgYuB12Tmyog4CHg3sAJ4fWaeN7iIp7eImAN8FtgZuBN4A3Ar8HVgfeDrmXloc67j3mfdxj8zL27aTgMOy8xrm8eOfx+N8t7/A/BlYB3gKuAVmbnMse+vUcZ+HvAp4GHAj4FXZuYKx76/xvmdsy1wVmZu3Dx27PtolPd9Ab4ILG9OOyozv+rY99coY/9j4BRgW+BGYK/M/L1jPzVMDjWuiCjAMcB+wOebw8cDb8rMSyLiI8DLI+K/gbcD2wCPAr4NPHkAIc8Ue1CP46bU4/gV4FfAUcBZwFkRsRNwNY57Gx40/hFxAPBOYHfgMICI2ADHv9+6vfdvBI7OzG9GxAeB10fESTj2/dZt7B8GvCgzr4+IrwK7RcSFOPb91m3st4mIucBHqZN0f+e0o9vYf4z6d87nhk9y7FvRbey/AVyVmXtExFHAv0bE8Tj2U8JpperVedSf4gzbLDMvaX4+H3g+9ac+Z2fm0sy8AbgtImJqw5xRNgROyswqM68B5gPbUX96XAGnAi/AcW9Lt/G/HTiZunI1zPHvv25jvxZwRtN+MbAljn0buo39u5rEcB6wLnA3jn0bHjT2zYezb6Oumg9z7Puv2/t+E+CmEec59v3Xbex3Bz7RtB9LXcF17KeIlUONq0lEzoqI9an/IAO4OSL+Hvg+8BL+7xfpko6ut1D/R59TGO6MkZnHDf8cEa8B7gJWNP8eUI/vjs1xx73Puoz/dZl5E3BTRLyt41Tf9302ytg/v3m8NnAI9R/Ljn2fjTL2p0XEc6lnLNwKXAIsxLHvq25jDxC0vaIAAAcKSURBVGwB7JiZu0TEB5pm3/d9NsrYbwrsHBGfBH5EPd3Rse+zMd73b4mI3YFfAK/BsZ8yJoearFcDi6k/Rb4EuBeogJUd51TA/VMf2swREesAHweeDewP/GdH8/D4Ou4tGTH+u41ymuPfgm5jHxFPp/4E+QLgROqKimPfZ93GPjO/FxHrAcdRj/t9OPZ9N2LsX0I93m8acZq/c1rQ5X3/QurfMxcD7weOpk5IHPs+6zL2vwAyM7eKiMOpx/4GHPsp4bRSTdZfZebfZ+a2wEXANcDNwEYd52wE/GYQwc0EEfFw4HvAMurppD8C1us4ZXh8HfcWjBz/zLx+lFMd/z7rNvYRsTNwJvD2zHxDZj6AY993XX7vLG8qJ2TmSuqpvZvg2Pddl7G/nXpDjjMj4lpg/Yi4DMe+70b5ff/FzLyoma3zRWBrHPu+G2Xs76ReOgNwGhA49lPGyqEm6yMRcSj1jlKvo57mdRPwzog4mvqPh3Uz0/9wJ+91wEWZ+ZbhAxFxRbMJzYXAy4F3Uf9ydNz770HjP4rv4vj3W7exPwbYPTN/1HHMse+/VcY+Iu6g3oDmo83Y/hNwKY59G0a+75cDjxlujIjfZuYzI2I+jn2/dfudc2lE7JuZP6Wu4vq+b0e3sT8feDH17uzPA36IYz9lTA41WW+l3nr44cDizPwxQEQcB1xJvc3wqwYX3ozwDGCHiHhBx7F/ot7Jaz3gxOE/lB33Vjxo/DNzy5EnZebtjn/fjRz7edTrf77Usf/AaZl5mGPfd91+77wJODciKupP+D+fmX927PvO3zmD0+19fzD175y51JuQHZSZSx37vus29jsBJ0bEe6nXFO6fmXc79lOjVFU1/lmSJEmSpBnNNYeSJEmSJJNDSZIkSZLJoSRJkiQJk0NJkiRJEiaHkiRJkiS8lYUkaYZpbrlwf/OwAn4JfDAzvzAFz/0U4FvATZm5oOP4JsA1wEsy87yO42+kvjXQkzPzvj7HsgD4X+D+zHzQ/+8jYgg4Hng+cABwcmYe0M8YJEnTi5VDSdJMtEWTED0KeB/wnxHxkil43u2B2zoTQ4DmZs3vA46JiDkAEfFI4CjgkH4khs392Ea6oFtiOCK2g4CDVvf5JUnTn5VDSdKMlZn3UN/Iemvg3cAZTXL2ceBfqP8/+D/UN1TeCLgCWC8z7wWIiCXAKzPz/M7rRsQrm+utD5wNvBHYBvgcUCJiaGSCCHwM2B84EDih6f+DzDyjueZmwGeAZ1FXGV+bmT9t2o5onmNd4AfAAZm5JCJOBO4Cdga+CPzHaGMREU8CvgRsDXwXeHgvYyhJmj2sHEqSZoP/Bp4eEfOAFwELgC2BTYH5wP6ZeS3wa+B5ABHxNGAt4MLOC0XEs4APA3sDTwBWAosy8wLqCtx3uySGZOafgX8F3hsR2zbnvqm55hzgTOAbwIbAScBpTdtWwELqpHED4HbgLR2X3gv458wcNTFsfIY6EV6/ufb245wvSZplTA4lSbPBbdT/z1sPuAzYFVgKPBpY3hyHOml6cfPzi4FvZebKEdc6APh0Zl6RmXcCRwAvjYgyXhBNBfKC5uvYzLy+adoemJuZx2fmvZm5GJjTJJE3Ac8FllAnjis64gX4YmZeM9bzRsRjge2Ao5rrn0hdnZQk6S+cVipJmg3mU29S83tgE+DzwGOBnwN/03HeacA3m593AT7U5VqbApd0PL6dusL4Vz3G8l7gpaw6BXRTYMuIWN5xbC51Mngz9TTYpwIJzANu7Djvrh6e83HUayFXdBy7pcd4JUmzhJVDSdJs8I/ARc3UzqOBSzJz88zcFfhVx3mXA0TETtRr887tcq3bqZOtYU8E7srMpT3Gci/wQGZ2JoK/Ay7PzLWGv4BnUk9pfQtwD7BZZr4AuLTH5+l0G/DoiHhYx7FNJnEdSdIMZuVQkjRjRcRa1FNI3wbs0RxeA3hY0/Y86sTxVxFRMrOKiDOAY4HvjEjghp0KHBcRp1NX346iXiO4On4AbBARu1InpHtRr2t8fBPvPGCtiHgG9UY6F/YyjXVYZt4QEdcAR0bEh4B9gS1WM2ZJ0gxj5VCSNBP9MiJWUq8rPBo4MDPPadreS31vv98D+1Fv9vIW6t1GoZ5aui3wtW4XzsxvAZ+i3vFzCfAA9c6jk9YkoXs017mT+t6HezTTQD9BPQX299TrG98A7Ead9E7Ey6h3Nb2dOik+Z+zTJUmzTamqatAxSJL0kBERj6HerOUxo1QOp4WIWAC8p9vOqV3OPQBYkJkHtBuVJOmhzMqhJEmNiFiH+hYTX57OiaEkSZNhcihJ0v85DXg5Y9xMfprZsZleO6qI+Czw2SmKR5L0EOa0UkmSJEmSlUNJkiRJksmhJEmSJAmTQ0mSJEkSJoeSJEmSJEwOJUmSJEnA/wc1WTlLC6cwKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "for add_unc in [0.01, 0.005, 0.001, 0.0005]:\n",
    "    plt.figure(figsize=(15,4))\n",
    "    colors = [\"#FC8D62\", \"#66C2A5\", \"#8DA0CB\"]\n",
    "    prior_mean = np.array([0.2, 0., 0.])\n",
    "    prior_std = np.array([0.5, 0.5, 0.5])\n",
    "    doy_start = 181\n",
    "    doy_end = doy_start + 8\n",
    "    while doy_end <= 273:\n",
    "        f, fwd, obs, rmse, r2, the_unc, f_upper, f_lower, sigma_f = fit_period_prior(\n",
    "            doys, qa, rho, doy_start, doy_end, 0.005, kern, prior_mean=prior_mean,\n",
    "            prior_std=prior_std)\n",
    "        t0 = 0.5*(doy_start + doy_end)\n",
    "        for i in range(3):\n",
    "            plt.vlines(t0, f[i]-1.96*sigma_f[i], f[i]+1.96*sigma_f[i], \n",
    "                       color=colors[i], lw=4, alpha=0.7)\n",
    "        prior_mean = f\n",
    "        # Add some uncertainty!\n",
    "        prior_std = sigma_f + np.ones(3)*add_unc\n",
    "        doy_end = doy_end + 8\n",
    "        doy_start = doy_start + 8\n",
    "    _ = plt.xlabel(\"Day of Year [d]\")\n",
    "    _ = plt.ylabel(\"Kernel weight [-]\")\n",
    "    _ = plt.title(f\"Prior uncertainty inflation {add_unc:g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the kernel weight uncertainties shrink as we decrease the extra inflation passing information around. An effect of this is that we see how the change of reflectance due to the fire is stretched over and appears smoothed: as we trust the prior retrieval more and more, the system tries to balance this with the goodness of fit. So trusting the previous time step too much results in changes taking longer (e.g. requiring extra evidence)  to appear in the data.\n",
    "\n",
    "In the next notebook, you'll see ways to exploit this for fame and wealth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}