scripts/notebooks/reconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstruction of 1, 2, 1+2 Neutron Responses for Constant P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H:\\GitHub\\pytracer\n"
     ]
    }
   ],
   "source": [
    "%cd H:\\GitHub\\pytracer\\"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scripts.utils import nice_double_plot\n",
    "import pytracer.algorithms as algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "response_single = np.load(r'scripts\\data\\fission_response_single.npy')\n",
    "response_double = np.load(r'scripts\\data\\fission_response_double.npy')\n",
    "single_probs = np.load(r'scripts\\data\\single_probs.npy')\n",
    "double_probs = np.load(r'scripts\\data\\double_probs.npy')\n",
    "\n",
    "d1 = single_probs.reshape((-1))\n",
    "G1 = response_single.reshape((response_single.shape[0], -1)).T\n",
    "d2 = single_probs.reshape((-1))\n",
    "G2 = response_single.reshape((response_double.shape[0], -1)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single Neutron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.logspace(-5, -1, 100)\n",
    "norms, residuals = algorithms.trace_lcurve(d1, G1, alphas)\n",
    "plt.loglog(residuals, norms)\n",
    "\n",
    "plt.figure()\n",
    "c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)\n",
    "plt.semilogx(c_alphas, curve)\n",
    "\n",
    "m1_alpha = algorithms.solve_tikhonov(d1, G1, 0.007)\n",
    "m1_alpha = m1_alpha.reshape((30, 50))\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(m1_alpha, interpolation='none', vmax=1.0, vmin=0)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Double Neutron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.logspace(-5, -1, 100)\n",
    "norms, residuals = algorithms.trace_lcurve(d2, G2, alphas)\n",
    "plt.loglog(residuals, norms)\n",
    "\n",
    "plt.figure()\n",
    "c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)\n",
    "plt.semilogx(c_alphas, curve)\n",
    "\n",
    "m2_alpha = algorithms.solve_tikhonov(d2, G2, 0.007)\n",
    "m2_alpha = m2_alpha.reshape((30, 50))\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(m2_alpha, interpolation='none', vmax=1.0, vmin=0)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single + Double Neutron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0109854114199\n",
      "[  1.00000000e-06   1.32571137e-06   1.75751062e-06   2.32995181e-06\n",
      "   3.08884360e-06   4.09491506e-06   5.42867544e-06   7.19685673e-06\n",
      "   9.54095476e-06   1.26485522e-05   1.67683294e-05   2.22299648e-05\n",
      "   2.94705170e-05   3.90693994e-05   5.17947468e-05   6.86648845e-05\n",
      "   9.10298178e-05   1.20679264e-04   1.59985872e-04   2.12095089e-04\n",
      "   2.81176870e-04   3.72759372e-04   4.94171336e-04   6.55128557e-04\n",
      "   8.68511374e-04   1.15139540e-03   1.52641797e-03   2.02358965e-03\n",
      "   2.68269580e-03   3.55648031e-03   4.71486636e-03   6.25055193e-03\n",
      "   8.28642773e-03   1.09854114e-02   1.45634848e-02   1.93069773e-02\n",
      "   2.55954792e-02   3.39322177e-02   4.49843267e-02   5.96362332e-02\n",
      "   7.90604321e-02   1.04811313e-01   1.38949549e-01   1.84206997e-01\n",
      "   2.44205309e-01   3.23745754e-01   4.29193426e-01   5.68986603e-01\n",
      "   7.54312006e-01   1.00000000e+00]\n",
      "[ 11.81618787   1.18161878   1.18161877   1.18161875   1.18161872\n",
      "   1.18161866   1.18161856   1.18161839   1.18161808   1.18161753\n",
      "   1.18161658   1.1816149    1.18161195   1.18160677   1.18159766\n",
      "   1.18158166   1.18155355   1.18150419   1.18141759   1.18126585\n",
      "   1.18100053   1.18053846   1.17973912   1.1783723    1.17608056\n",
      "   1.17236011   1.1666196    1.1583995    1.14774239   1.13544879\n",
      "   1.1228788    1.11140271   1.10192007   1.09457507   1.08865002\n",
      "   1.08274445   1.07513446   1.06404513   1.04784332   1.02528647\n",
      "   0.99562345   0.95802197   0.91042472   0.84929328   0.77134904\n",
      "   0.67606697   0.56678826   0.45031636   0.33624686   0.23517299]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xt0pHWd5/H3tyqpyrXTufX9kqa5\nCKI0EBDBYdABl3Udb+NlcFZhlhF11z26617QOWfXmbOz6+4RnZ3L6mkPjJw9KnqQGdRhR0BhEGSg\n09hA001D3+1rLp17UrlU/faPqkqHdDqpSj11eX75vA45lXrqqdT3RzqfPPnV7/k+5pxDRET8Eyl3\nASIiUhwKeBERTyngRUQ8pYAXEfGUAl5ExFMKeBERTyngRUQ8pYAXEfGUAl5ExFMKeBERT1WV8sXa\n2tpcR0dHKV9SRCT0du7c2euca8/3eSUN+I6ODrq6ukr5kiIioWdmR5byPE3RiIh4SgEvIuIpBbyI\niKcU8CIinlLAi4h4SgEvIuIpBbzIMnG0b4zB8alylyElpIAXWQYmp1O876+f5s8ff63cpUgJKeBF\nloHnDvUxMDbFsf7xcpciJaSAF1kGHt9zGoDekYkyVyKlpIAX8ZxzjscyAd8zrIBfThTwIp7bc3KI\nE4MJ2hpi9I5M4Jwrd0lSIgp4Ec89vqcbM/jQVRtITKUYnUyWuyQpEQW8iOce23uKKzeu5E1rGgHo\n1TTNsqGAF/HYycFxdh8f4pbL1tDWEAegR2+0LhuLBryZ1ZjZ82b2opm9YmZ/ktm+xcyeM7PXzewH\nZhYrfrkiko/H93YDcMtlq2YCXkfwy0cuR/ATwLucc1cA24Bbzew64H8C33DOXQT0A3cWr0wRWYrH\n9pxmS1s9W9sbaG/MBLyO4JeNRQPepY1k7lZnPhzwLuDBzPb7gQ8UpUIRWZLhxBTPHujl5ktXYWa0\n1MeIGPSMTJa7NCmRnObgzSxqZruAbuAx4AAw4JybzuxyDFhfnBJFZCmeeq2XqaTj5ktXAxCNpENe\na+GXj5wC3jmXdM5tAzYA1wKXzrfbfM81s7vMrMvMunp6epZeqYjk5fG9p2muq+bqzc0z29oa4pqi\nWUbyWkXjnBsAngSuA1aaWfai3RuAE+d5znbnXKdzrrO9Pe+LgovIEkwnU/zi1W7e+aZVVEXP/pgr\n4JeXXFbRtJvZyszntcDNwF7gCeDDmd1uBx4uVpEikp8dh/sZHJ/ilsz0TFZ7owJ+OalafBfWAveb\nWZT0L4QfOud+amZ7gAfM7L8BvwbuLWKdIpKHx/eeJhaNcOPFb/yrua0hPQfvnMPMylSdlMqiAe+c\newm4cp7tB0nPx4tIBXHO8fje01x/YSv18Tf+iLc1xGfaFTTEczm+kzDTmawinnm9e4QjfWMzq2dm\n08lOy4sCXsQz2dbA8wW8TnZaXhTwIp55bM9p3rqhiTVNNec8NtOPRkfwy4ICXsQj3cMJdv1m4JzV\nM1ltjemWUTqCXx4U8CIe+UWmudjNl80f8C11MUztCpYNBbyIR365v5d1TTUzvd/nqopGaK2P6Qh+\nmVDAi3jkWP84W1c1LLjGva0hrjn4ZUIBL+KRU4PjrFlx7purs6ldwfKhgBfxxFQyRffwBGtX1i64\nX/bi2+I/BbyIJ7qHJ3AO1s6zPHK2toY4vcOTODdvA1jxiAJexBOnBscB5l3/Plt7Y5zxqSSjk8lS\nlCVlpIAX8cTJwQSQ2xE8qF3BcqCAF/HEqZmAX2QOXu0Klg0FvIgnTgwkqItFWVGzcJfItgadzbpc\nKOBFPHFqaJw1TTWL9nnPNhzTWnj/KeBFPHFyMLHo/DuoXcFyooAX8cSpwQRrViw8/w7pdgUtdVoL\nvxwo4EU8MJ05yWndysWP4CG7Fl4B7zsFvIgHekYmSKbcomvgs9oaY/ToCN57CngRD+S6Bj6rXf1o\nlgUFvIgHsmvgc5mDh7PtCsRviwa8mW00syfMbK+ZvWJmn89s/4qZHTezXZmP9xS/XBGZT/YIPuc5\n+Gy7gonpYpYlZbbwGRFp08AXnXMvmFkjsNPMHss89g3n3NeKV56I5OLU4Dg11RGaaqtz2n/2tVnr\n47nEgITRokfwzrmTzrkXMp8PA3uB9cUuTERyd2Iwwdqm2kVPcspqV7uCZSGvOXgz6wCuBJ7LbPqc\nmb1kZveZWXPAtYlIjtJr4HObngG1K1gucg54M2sAfgR8wTk3BHwT2ApsA04C95zneXeZWZeZdfX0\n9ARQsojMdSrHs1iz2rNTNDqb1Ws5BbyZVZMO9+865x4CcM6dds4lnXMp4NvAtfM91zm33TnX6Zzr\nbG9vD6puEclIphynhxKszfENVoCW+ky7Ap3s5LVcVtEYcC+w1zn39Vnb187a7YPA7uDLE5HF9I5M\nMJ1yrFmkTfBsalewPOTy9vkNwCeAl81sV2bbl4HbzGwb4IDDwKeLUqGILGjmJKc85uBB7QqWg0UD\n3jn3NDDfW/OPBF+OiOQr10v1zdXWqCN43+lMVpGQy7dNQVZbQ1z9aDyngBcJuVODCWJVEVrqY3k9\nT+0K/KeAFwm5E5klkrme5JTVrnYF3lPAi4TcqcHxvE5yysq2K9A8vL8U8CIhl+ul+ubKns2qtfD+\nUsCLhFhq5iSn3NfAZ+kI3n8KeJEQ6x2dYCrplnQEv6pR7Qp8p4AXCbGzF/rIP+Cz7Qp0spO/FPAi\nIXZ2DXz+UzRV0QjNdbo2q88U8CIhNnMEv4QpGki/0aojeH8p4EVC7ORgglg0QmueJzlltTfq4ts+\nU8CLhNjJwXFWN8WJRPI7ySmrrSFOr95k9ZYCXiTETg4mWLsi//n3rLaGuNbBe0wBLxJipwYTS55/\nh3TAq12BvxTwIiHlnEtfqi+PKznNpWuz+k0BLxJSZ0YnmUym8r7Qx2ztjTqb1WcKeJGQOjmzRLKw\nOXhQPxpfKeBFQmqpF/qYrV3tCrymgBcJqeyl+goJ+OxFQnSyk58U8CIhdXIwQVXEZqZZlqI6mr4S\nlObg/aSAFwmpk4MJVq+oWfJJTlltDTHNwXtKAS8SUicHxwuanslKn82qgPfRogFvZhvN7Akz22tm\nr5jZ5zPbW8zsMTN7PXPbXPxyRSSr0JOcstSuwF+5HMFPA190zl0KXAf8GzO7DLgb+Llz7iLg55n7\nIlICzjlODiZYt4QrOc2lI3h/LRrwzrmTzrkXMp8PA3uB9cD7gfszu90PfKBYRYrIGw2MTTExnVrS\nhT7mam+MMzapdgU+ymsO3sw6gCuB54DVzrmTkP4lAKwKujgRmd+JAJZIZqldgb9yDngzawB+BHzB\nOTeUx/PuMrMuM+vq6elZSo0iMkehF/qYrU3tCryVU8CbWTXpcP+uc+6hzObTZrY28/haoHu+5zrn\ntjvnOp1zne3t7UHULLLsFXKpvrnaZ9oV6I1W3+SyisaAe4G9zrmvz3rox8Dtmc9vBx4OvjwRmc+p\nwQTRiM20GihEa2aKpm9UR/C+qcphnxuATwAvm9muzLYvA18FfmhmdwJHgY8Up0QRmevkYILVjXGi\nBZ7kBGfbFZzRUknvLBrwzrmngfP9K/qdYMsRkVycHBwPZP4dIF4VpTFeRd+oAt43OpNVJIRODSYC\nmX/PammIcUYB7x0FvEjIZE9yCmKJZFZLvQLeRwp4kZAZGp9mfCoZ2BQNQGt9TFM0HlLAi4TM2ZOc\nApyiqY/Rp3Xw3lHAi4RMkCc5ZbXUx+kfm8Q5F9jXlPJTwIuETBCX6purtT7GVNIxlFA/Gp8o4EVC\n5siZUWLRCKsDaDSWNbMWXvPwXlHAi4TMoZ5RNrXWBXKSU1ZLQzbgNQ/vEwW8SMgc7htlS1t9oF+z\nNXME36ezWb2igBcJkWTKcbhvLPCA1xSNnxTwIiFyYmCcyelUEY7g003LtBbeLwp4kRA53DcKEHjA\n18ai1FZHdQTvGQW8SIgc6k0H/AUBBzyoXYGPFPAiIXKwZ5T6WDSQPvBztTWoXYFvFPAiIXKod5SO\ntnrS1+EJVvoIXsskfaKAFwmRYiyRzGqpj+uiH55RwIuExOR0it+cCX6JZFZrZopG/Wj8oYAXCYmj\nZ8ZIueBX0GS11MeYmE4xNpksyteX0lPAi4TE4d7iLJHM0slO/lHAi4TEoSIH/Ey7AgW8NxTwIiFx\nsHeU5rpqVtbFivL1zx7BayWNLxYNeDO7z8y6zWz3rG1fMbPjZrYr8/Ge4pYpIod7i7eCBma1K9BK\nGm/kcgT/HeDWebZ/wzm3LfPxSLBlichc2TXwxXK2ZbAC3heLBrxz7ingTAlqEZHzGJ2Y5tRQoigt\nCrLqY1FiVRHNwXukkDn4z5nZS5kpnObAKhKRc5xtMtZQtNcwM1rrY5qi8chSA/6bwFZgG3ASuOd8\nO5rZXWbWZWZdPT09S3w5keXtcO8YULwVNFlqV+CXJQW8c+60cy7pnEsB3wauXWDf7c65TudcZ3t7\n+1LrFFnWDvWOANDRVlfU11FHSb8sKeDNbO2sux8Edp9vXxEp3MHeUdasqKEuVlXU12mtV0dJnyz6\nr8XMvg/cBLSZ2THgvwI3mdk2wAGHgU8XsUaRZa/YSySzWurjOoL3yKIB75y7bZ7N9xahFhE5j0O9\no9x6+drFdyxQa0OMsckkiakkNdXRor+eFJfOZBWpcANjk/SPTRV1iWSW2hX4RQEvUuGK3YNmtpl2\nBVoq6QUFvEiFywZ8Mc9izWptyB7Ba6mkDxTwIhXuUO8oEYNNLcVdIgnpN1lB7Qp8oYAXqXAHe0fZ\n2FJHrKr4P67qCe8XBbxIhTvcO0pHa/GnZwBW1FRRHTW9yeoJBbxIBXPOcahEa+Ah3Y+muS6mN1k9\noYAXqWDdwxOMTSa5oL00AQ/paRodwftBAS9SwUq5RDKrtUENx3yhgBepYDNLJEs0Bw9qV+ATBbxI\nBTvUO0qsKsK6lbUle001HPOHAl6kgh3sGaWjtY5oxEr2mi31MYYT00xOp0r2mlIcCniRCna4r3RL\nJLO0Ft4fCniRCpVMOY70jbKlhCtoYHbDMb3RGnYKeJEKdbx/nKmkK0kXydl0BO8PBbxIhTrUV/oV\nNACtDepH4wsFvEiFOtSTvg5r2aZodDZr6CngRSrUod5RGuJVtGeOqEulqbaaaMR0BO8BBbxIhTrU\nN0ZHWx1mpVsiCRCJGM111VoL7wEFvEiFOtQ7wpa2hrK8dku92hX4QAEvUoEmppMc6x8vaQ+a2dIB\nryP4sFPAi1SgI31jOAdb2op/Faf5tNbHNUXjgUUD3szuM7NuM9s9a1uLmT1mZq9nbpuLW6bI8vL8\noTMAXLFhZVleX0fwfsjlCP47wK1ztt0N/Nw5dxHw88x9EQnIrw70sq6ppqxTNANjU0wn1Y8mzBYN\neOfcU8CZOZvfD9yf+fx+4AMB1yWybCVTjl8d6OOGC9tKvoImq7UhvRa+f2yqLK8vwVjqHPxq59xJ\ngMztqvPtaGZ3mVmXmXX19PQs8eVElo89J4YYGJviHRe1la0GtSvwQ9HfZHXObXfOdTrnOtvb24v9\nciKh9/T+XgDevrW1bDW0qOGYF5Ya8KfNbC1A5rY7uJJElrdfHejlktWNrGqsKVsNrfXqR+ODpQb8\nj4HbM5/fDjwcTDkiy1tiKsnzh85ww4Xlm54BTdH4Ipdlkt8HngUuMbNjZnYn8FXgFjN7Hbglc19E\nCvTCkX4mplO846LyTc8ANNdVA2o4FnZVi+3gnLvtPA/9TsC1iCx7T+/vpSpiXLulvAFfFY2wsq5a\nR/AhpzNZRSrIMwf62LZxJQ3xRY+9ii598W29yRpmCniRCjE4NsXLxwbKPv+e1Vof1xRNyCngRSrE\nswf7SDkqJuDVriD8FPAiFeKZ/b3UxaJs21ie/jNztTQo4MNOAS9SIZ450MvbtrQQq6qMH8vW+hj9\nY5OkUq7cpcgSVca/JJFl7sTAOAd7RitmegbSUzQpBwPj6kcTVgp4kQrwTKY9QaUFPKArO4WYAl6k\nAjyzv5e2hhiXrG4sdykzsu0KtJImvBTwImXmnOOZA31cv7WNSKQ87YHno3YF4aeAFymz17tH6Bme\n4B0VND0DZ3vC69J94aWAFymzp19Pz79ff2F52xPM1VynI/iwU8CLlNkz+3vpaK1jQ3N5LrB9PrGq\nCI01VQr4EFPAi5TRVDLFPx3sq6jVM7Ol+9Eo4MNKAS9SRi8dG2B0Mllx8+9Z6XYFWiYZVgp4kTJ6\n+vU+zMp7eb6FtKjhWKgp4EXK6Jn9vVy+romVmTc0K02b+tGEmgJepEwO946y82g/N15cmdMzkJ6i\n6R+bxDn1owkjBbxImfz5469RHTVuv76j3KWcV0t9jKmkYygxXe5SZAkU8CJlsO/UMA+/eILbr+9g\nVWNNucs5rw3NtQAc7BkpcyWyFAp4kTK459F9NMSq+MyNW8tdyoKu2tQMwM4j/WWuRJZCAS9SYi/+\nZoBH95zmj37rAprrK/PN1axVK2rY3FrHjsNnyl2KLEFBV/Y1s8PAMJAEpp1znUEUJeKzrz26j+a6\nav7VOzrKXUpOOje38OS+bpxzmFVOMzRZXBBH8O90zm1TuIss7rmDffzy9V4+e9NWGmuqy11OTq7p\naKZvdJJDvaPlLkXypCkakRJxzvG1R/exekWcT769o9zl5KyzowWArsOahw+bQgPeAY+a2U4zuyuI\ngkR89eRrPew43M/n3nURNdXRcpeTs63t9TTXVWsePoQKmoMHbnDOnTCzVcBjZvaqc+6p2Ttkgv8u\ngE2bNhX4ciLh5Jzjnkf3saG5lo91bix3OXkxMzo7WujSSprQKegI3jl3InPbDfwtcO08+2x3znU6\n5zrb29sLeTmR0PqH3afYfXyIL9x8MbGq8M2MXtPRzKHeUXqG1XgsTJb8L83M6s2sMfs58G5gd1CF\nifgimXLc89hrbG2v54NXri93OUuSnYffeUTTNGFSyKHEauBpM3sReB74e+fcPwRTlog//u7Xx9nf\nPcIX330J0Qq65mo+Ll/XRLwqwg690RoqS56Dd84dBK4IsBYR7zx3sI//8vBu3rK+iVvfvKbc5SxZ\nrCrCto0r6dIbraESvslAkZD4x9d6uP1vnmdNUw3f/mQnkZAevWdd09HC7hNDjE2q8VhYKOBFiuBn\nr5ziU/d3cUFbAz/89NtZ01S5DcVy1dnRTDLl2HV0oNylSI4U8CIBe3jXcf71d1/gsnUr+P6nrqO1\nIV7ukgJx1eZmzNA8fIgUug5eRGZ54PmjfOlvX+bajhbuveMaGuL+/IitqKnmktWNdGklTWjoCF4k\nIPc9fYi7H3qZGy9q5zt/eK1X4Z51TUcLLxzpZzqZKncpkgMFvEiBXj01xOcf+DV/+tM93PrmNWz/\n5NXUxsLTiiAfnR3NjE4mefXUcLlLkRz4d4ghUiI7Dp/hm08e4BevdlMXi/KZ397Kf3j3xVRF/T1u\numam8dgZLl/fVOZqZDEKeJE8pFKOJ/Z1880nD9B1pJ/mumr+/S0X88m3b2ZlXWVfvCMI61bWsn5l\nLTuO9HPHDVvKXY4sQgEvsoju4QQvHBnghaP9PPFqN693j7B+ZS1f+d3L+Og1G6mLLa8fo86OZv7p\nYJ8uABICy+tfpsgCUilH3+gkJwbGeen4IC8c6WfnkX6OnhkDIBaNcMXGJr7xsSt471vXUe3xVMxC\nOjtaeHjXCY71j7Oxpa7c5cgCFPDiBecc0ynH5HSKqWSKyekUk8kUiakUoxPTjGQ+RjMfIxNJ+kYm\nODmU4PRggpODCbqHE0wl3czXbG+Mc/WmZj5x3Wau2tzM5etXEK/y883TfFzTkb4Q947DZxTwFc7b\ngHcuPVd6uHeMwfEpBsenGMrcDo5PMTIxTTLlcEDKOUj/h3PpH/BoxKiKRIhEjKqIZe4bZuBm7Zu+\nXbye2X/JluuP2hzKPPc5CzzpnIdm7exmbXKZe9mH3Zz/19l93vD/1aW/L6nMrXPproxJ50il0mGe\nvU06x1QyldP3Ybaa6ghrm2pZs6KGt21pYXVTDWubali9oobL1q5gQ3OtpiDmcfGqRhprqthxuJ8P\nXbWh3OXIArwM+IGxSf7Tgy/x6J7TM9sa41WsqK2mKfOxqaWOaCawzQzj7C1A0jmSyUyQZI4OZ6/9\nNQPDZoJ7oSBwCySPc28M/0q00NjmPjLfL7LZ/1/PPp75f88b/19mvx8RMyIGkexzjZlftNGIETUj\nGokQjUAkYsSjEWJVEapn30YjxKsjNMSrqI9X0ZD5yH5eUx1RgC9BJGJ0bm5W47EQ8C7guw6f4fMP\n7OL0UIIvv+dNfOTqjTTWVHm9dE2k1Do7Wnhi3z76Rydprvd/9VBYeRPwyZTjW/94gK8/9hrrV9by\no89ezxUbV5a7LBEvXTNzAZB+br5sdZmrkfPxIuC7hxL8ux/u4pn9ffzuFev47x+8nMaa6nKXJeKt\nt25oIhaNsOPIGQV8BQtlwE8lU5waTHB8YJwDPSN8/dHXGJ2c5n/93lv5SOcGzauKFFlNdZS3bGji\n6dd7mX53SlOgFSoUAf/wruP84tVujvePc2JgnFNDCVKz3rd805pGfvDx67hwVWP5ihRZZj589Qa+\n9NDLfO57v+Z/37ZNS0grUCgC/tVTw7xwtJ91TbVct7WVDStrWd9cO3Pa9ObW+tBe61IkrG67dhOJ\nqSR/8pM9/NH9XXzrX15NvYcdNMPMFlrCF7TOzk7X1dVVstcTkeL70c5j/McHX+SKjSv5zh3X0lSn\n97+CZmY7nXOd+T5PE2ciUpDfu3oD/+cPruaV40N8bPuzdA8nyl2SZBQU8GZ2q5ntM7P9ZnZ3UEWJ\nSLjcevka7rvjGo6eGeOj33qWY/1j5S5JKCDgzSwK/DXwz4HLgNvM7LKgChORcHnHRW383zvfxpnR\nST78zWd5cl83g+NT5S5rWSvkHZFrgf3OuYMAZvYA8H5gTxCFiUj4XL25mR98+u188r7nueNvdgCw\npa2ey9c38db1TbxlQxNvXrdC56mUSCEBvx74zaz7x4C3FVaOiITdpWtX8Isv/ja/PjrAy8cHeflY\nuvXyT148MbNPXSxKY00VjTXVb7yNV1FTHSVeFSGevc18XjN3W1WUeHWEmqooNdURamNRaqujM8/X\n+TCFBfx8//fOWZJjZncBdwFs2rSpgJcTkbBorKnmxovbufHi9pltfSMTvHx8kFdODNE/OslwYprh\niSmGE9MMjk9xrH+MkcQ0E9MpJqaTJKaWfmFvM6itTgd+fbyKxpp0g7nsL5OGeBUraqtoqY/T1hCj\nrSFOa+a2uS7mzbLrQgL+GLBx1v0NwIm5OznntgPbIb1MsoDXE5EQa22Ic9Mlq7jpklU57e+cYzKZ\nSgf+VIrEVDLT4z85sy37iyB9m2R8MsnYVJLEZJLxqSRjk8mZ6wEMJaY5MZBgeGKYkUT6fjJ1biSZ\nQWt9nA3NtZmPOjY0p8+92Zi5X1MdjpO6Cgn4HcBFZrYFOA78PvDxQKoSkWXPzNLTMFVRqAn+66dS\njqHEFL0jk/SOTNA3czvB6aEJjg+Ms/v4ID975dQbLgQDsHpFnM0t9WxqraOjtY5NrfVsbqlj1Yr0\nXwCV8gtgyQHvnJs2s88BPwOiwH3OuVcCq0xEpIgiEWNlXYyVdTEuXNVw3v1SKUfPyATH+sc41j/O\n0b4xjpwZ42jfGE+91sODwxPnPKcuFqW5LkZLfYzm+hgtddV86sYLePO6pmIO6RwFnVfsnHsEeCSg\nWkREKk4kYqxekb7S19Wbz318fDLJ0TNjHOkbpXdkkv6xSc6Mpm/7Ryc5MzbF4d5RxiaTJa9djSNE\nRApQG4tyyZpGLllTec0O1apARMRTCngREU8p4EVEPKWAFxHxlAJeRMRTCngREU8p4EVEPKWAFxHx\nVEmvyWpmPcARoA3oLeBLNQGDBe4332OLbZv7ePb+7O2Fju18deSzz1LGNvf+fJ9X6tjm216pY1ts\nv3y/T3Pvn+9zH37mfBvb3PsL/bvc7JxrJ1/OuZJ/AF0FPn97ofvN99hi2+Y+nr0/Z5+Cxpbr+IIe\n20JjnTXOihzbYmOppLEF8b3Ld2xBjK8SfuZ8G1su4yt0bGGdovlJAPvN99hi2+Y+/pPzbC9ULl8v\n6LHNvb/QuAtRjLHNt71Sx7bYfvl+n+beD/vY5m7zeWxz7wc+vpJO0cy8qFmXc66z5C9cAhpbOPk8\nNvB7fBrb+ZXrCH57mV63FDS2cPJ5bOD3+DS28yjLEbyIiBRfWOfgRURkEQp4ERFPKeBFRDxVUQFv\nZhEz+zMz+0szu73c9QTNzG4ys1+a2bfM7KZy1xM0M6s3s51m9t5y1xIkM7s08z170Mw+W+56gmRm\nHzCzb5vZw2b27nLXEzQzu8DM7jWzB8tdSxAyP2P3Z75nf7DY/oEFvJndZ2bdZrZ7zvZbzWyfme03\ns7sX+TLvB9YDU8CxoGoLQkDjc8AI6WvEV8z4AhobwH8GflicKpcmiLE55/Y65z4DfBSomOV4AY3t\n75xznwLuAD5WxHLzFtD4Djrn7ixupYXJc5wfAh7MfM/et+gXL+QsqTlnZN0IXAXsnrUtChwALgBi\nwIvAZcBbgJ/O+VgF3A18OvPcB4OqrYLGF8k8bzXw3XKPKeCx3Qz8PumgeG+5xxTk2DLPeR/wK+Dj\n5R5T0GPLPO8e4Kpyj6mI46uoPClgnF8CtmX2+d5iXzuwi247554ys445m68F9jvnDgKY2QPA+51z\n/wM45894MzsGTGbulv4S5AsIYnyz9APxYtS5FAF9794J1JP+RzhuZo8451JFLTwHQX3fnHM/Bn5s\nZn8PfK94FecuoO+bAV8F/p9z7oXiVpyfgH/mKlY+4yT9l/8GYBc5zMAEFvDnsR74zaz7x4C3LbD/\nQ8BfmtlvAU8Vs7CA5DU+M/sQ8M+AlcBfFbe0guU1NufcHwOY2R1AbyWE+wLy/b7dRPpP4zjwSFEr\nK1y+P3P/lvRfX01mdqFz7lvFLC4A+X7vWoE/A640sy9lfhGEwfnG+RfAX5nZvyCHdgbFDnibZ9t5\nz6xyzo0BFT1fNke+43uI9C+xMMhrbDM7OPed4EsJXL7ftyeBJ4tVTMDyHdtfkA6NsMh3fH3AZ4pX\nTtHMO07n3Cjwh7l+kWKvojneyRfDAAAAo0lEQVQGbJx1fwNwosivWUo+j09jCyefxwb+jy8rkHEW\nO+B3ABeZ2RYzi5F+E+7HRX7NUvJ5fBpbOPk8NvB/fFnBjDPAd4K/D5zk7BLHOzPb3wO8Rvod4T8u\n9zvWGp/GprGF48P38ZVinGo2JiLiqYo6k1VERIKjgBcR8ZQCXkTEUwp4ERFPKeBFRDylgBcR8ZQC\nXkTEUwp4ERFPKeBFRDz1/wGFGWJ3TZNABQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbd9af8a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmYFOW59/HvzQCCMAwCwyKIqCzu\nRhkNMTFucUEGjcYjak6OnhCXvDExiBrNoqjxxETU6IkRMRrimxPU5KgRUBODC0qMr6AYURMXQGRf\nRJYxw3q/fzzVQ9PTM9Pd0z1dM/37XFdfTVc9VXV3M92/fqqeqjZ3R0REJG7aFbsAERGRdBRQIiIS\nSwooERGJJQWUiIjEkgJKRERiSQElIiKxpIASEZFYUkCJiEgsKaAKyMwmmNmaHJc9y8yeNbNPzGyz\nmb1rZj82s175rrNQmvMczGyKmc1JenyOmV2Y5/p22UbcZPoa5ON5mNmFZvaKmdWY2QYze8HMTs9x\nXQWrM8s6Yv3/m2BmB5rZTDP71MyWmdmNZlaWj+XMbLCZ3Wtmb5jZdjN7Ps16mmxTLO2LXYDUZ2a3\nAd8Ffg3cAWwADgQuBQ4CzixedZnJw3O4Ceic9PgcoBcwJY9lpm4jblriNcDM7gEuAn4J/JDwuXAu\n8Eczu8bdf5rlKhuqM+6vd4szsz2AvwBvA2cA+wG3EToPP8zDcgcBpwF/Azo2sLpM2hSHu+tWoBsw\nAViT5TKjAQe+nmZeGTCymTWVAR0L/Lxzfg4N1Qf8AXg+T/UV/DUo0Oua9jUgBMGcHNf55ej/6tI0\n834KbAeOyEedRXi9cn5dWrDGa4F1QLekaVcDnyZPy3U5oF0Gfz9Ntina61PsAtryLceAehaYm2Hb\n54E/pEw7LvrAOTh6PAWYE30QvQVsBX4EbAa6pyx7ULTsiUnTvgC8EP3hrwXuA8rz+BzS1XdM8odL\n9G9PuU3ItMZMtpFS0znAm9Fr9BFwM9A+zfpOAv4O1AAvAQc18jxPiOreM2nay4QA6J407U3g5uTt\nNPUa5FJP0vaeA94DytLMqwA+AR5I8zr+A6iNtnNgymvTaJ1pXsdRhJ7Ap8AMoAcwOKqtJmpzaEpt\nnwOeAJZFbeYBX033t5XBa9ADuBdYA3wM/DSa/hLw7QJ/RswCHkqZNjB63UbnczkyCJ9M2rTkTceg\nYsTMOgBHA0/nedWDgJ8BPyF05R+MpqfuZhsDrCIEH2b2eWAmsAI4m7DL7jTCbru0cnwOqfUtTJl/\nE+HD6nXCB9PngF9lWWNT20jUfzLwMPAaYdfJfwNXAr9IaToQuJUQXucBvYFHzMwaeI4vszMYMbPd\ngeHAFuDz0bQehC8JL6ZZvsHXIMd6MLP20Xqmufv21Pnuvj7a5heTJu8N3B7Vcz4hxP5kZp0yrDPV\nQOBGwm6piwl/O5OBh6Lb2YRdjg+lPJe9gdnANwg99v8Ffm1m5zWyrXrMrDvwV8L/xdeBicDVZvZN\nwm6z+xpYzsysfVO3DErYnxD2ddx9MSGs9y/Acq2KjkHFS09gN2BxAdb7JXefl5hgZk8TAin5g3wM\n8PukD6tbgL+6+5ik5ZYCM83sYHefn6fnkK6+upnu/oGZfUzYFfG3lGUzrbHRbSS5kfAN8oLo8dNR\nu5+Y2Y/dfUk0vQfweXd/L1pXO+AxYBgpHxzRc/iXmc0lBNTDwAjCcbm/RNNmEHqCTvjATF2+sdcg\n63oivQj/Vx82MJ9o3qkpy5zh7n+NtjMX+AC4EJiUQZ3p6v6cu38Qre9Q4CrgAnd/MJpmhNdnf+Ad\nAHd/KLGCaP4sYADhWNrUDLab8AOgLzDY3dcAT5jZOMIXmZvcvbaB5S6gkS9qSRr8ghDZg9BLTbUu\nmpfv5VoV9aCKJM03sOTRN/n+DZSlyR/MkYeBExMj6szsM8DQaHriG/7nCN/Ck78RvkToCQxvYpvZ\nPId09TUpyxqb3Eb0f3AE8PuUWQ8T3iufS5q2KBEGkbej+wGNbOJFoh4UoVfyImHXZPK0N9x9Q2N1\nNiCXenKxKhFOAO7+ITAXOCrH9S1KhFPk/ej+2TTT+icmmNkeZnaXmX1I+L/eSuiBDc10w1GwXQDc\nF4VTwvpofZMaWXwacGQGt0yke69YA9PzsVyroYAqnmPZ+cbaSthNtZZw3GNgnre1Ms20J6LtnhU9\nHgMsJXy4Q/gWVkYY2ZVc52agA7BXA9vK5Tmkqy8T2dSYyTZ6Rcultk087pE0LfXb65bovhMNmwUc\nHO1WOoYQUC8CVdEussS0XORSzxrCa7V3I232JvxdJKxK02YV0K+pAhvQUN2fpJmW/FymEP5mbwVO\nJoTBAzT+fFMNAyqBZ1KmlwF3untNI8t+TDju1dStKeuA7mmmJ47/5Xu5VkW7+IpnLrt+w9ro7lvN\nbDZwCo0MMU1SS/1hoT3StKv3jcrdN5nZDMKbfDJhYMAjHh0pJfyRO2Ggx5Np1rksXUE5PIe09WUo\nmxoz2cYaQsD1TpneJ7r/OMv6Us0mfMM9jrCL73uEQRubgBMJvbdbm7mNjLn7NjN7GRhlZle6+47k\n+WbWLar1saTJqa9NYtpbBSs0RRTmo4DL3H1S0vRsv3DvE93X7eI0s6MJxyubCpd87eL7BynHjMxs\nL6ALDe+abc5yrYoCqkjcfSNhdFKqnxP2g1/g7r9JnhG9AU9298QAhCXsegAbwkiuTD0EPGxmo4F9\no8eJ+mrM7G/AMHe/MYt1ZvscMrWFlG/HzayxHnffHh1T+TfgnqRZ5wA7CAMdmrP+dWY2HxhHGL33\nuru7mb1EGCKc2D3ZkHqvQR7cSQigbxC+qCS7BujGrgNEepvZ0UnHoAYSgjX5w7oQdSbbjdDL2ZyY\nYGblwOlk92Uncay1R7QOIwykgaaDJbGLr7meAq4ys/LoMwHCl8Z/EXb/5nu5VkUBVXgdzezsNNNf\ncPfVqRPdfZqZ3Q7cH41Q+yPhG/b+hJNcF7FzhNxjwFgzu4NwEPl4Qs8lUzMIo37uBRa6+/9LmX81\nYbDBDsLw042EXXejgB+4+7vpVprlc8jUP4AzzOzLhGBe5u7Lcq2xEdcTRqX9mhDYhxBGpt2XNECi\nOWYB3wL+lDQY5UVCz+k9d1/RyLINvQY5c/fHzWwScLeZHQhMJ3wujCEMfLjW3V9LWmQN8H/N7EeE\nD8MbCbv4phSyzpSa15vZq8B1ZraB8OXhGsKxo25ZrGouYS/Ez8zsZsJz7kEYdn+umb3p7osaqGEt\nYXd2c00CvgM8amY/JXxRnADcnjgWaWb/Qdh9uV90zC+j5aJldyeMWoVwDK9b0ufRk+7+aSZt8vA8\nc1Psce5t+Ub4g0k9JyRxO66JZb9CGK67nvCN9F3CENi+Ke2uJZyrsxH4LTu/Re5yHlQj2/lt1P4n\nDcz/LCFMNhDON3mbMMy4IoPn3+RzaKi+1OmE40OPEXazpZ4H1WiNmW4jafoYwvlIWwgfsGnPg0pZ\nZlBUV3UTr8mYqN33U+p3ovONsn0NmlNP1NYIYfRK9PptJHwLPz1dPYTjlu8SejCzE39r2dbZQN0X\nRst0bey5EM6TejaqdzHhi8oEks47bOj/N2V75xBOOdhCCKx9CLvvNgFXF/LzIamGA6Pn8i9gOeEL\nUVnS/MRrMiib5VJeu3S3QZm2KdbNogJFRBplZlMIYVRV7FqkNGgUn4iIxJICSkREYkm7+EREJJbU\ngxIRkVgq+WHmvXr18kGDBhW7DBGRkjF37tw17l7ZVLuSD6hBgwYxZ07sf3RTRKTNiK6h2CTt4hMR\nkVhSQImISCwpoEREJJYUUCIiEksKKBERiSUFlIiIxJICKo/mz5/P0UcfzUUXXVTsUkREWr2SPw8q\nnzZt2sTLL7/Mjh07mm4sIiKNUg9KRERiSQElIiKxpIASEZFYanMBZWb7mtn9ZvaHYtciIiK5axUB\nZWYPmNkqM5ufMv1UM/unmb1vZtcAuPsCdx9bnEpFRCRfWkVAAVOAU5MnmFkZcDcwEjgQOM/MDmz5\n0kREpBBaRUC5+yzg45TJRwHvRz2mLcBDwBktXpyIiBREqwioBvQHPkp6vATob2Y9zWwScLiZXZtu\nQTO72MzmmNmc1atXt0StIiKSpdZ8oq6lmebuvha4tLEF3X0yMBmgqqrKC1CbiIg0U2vuQS0B9kp6\nPABYVqRaREQkz1pzQL0KDDGzfcysI3Au8ESRaxIRkTxpFQFlZlOBl4FhZrbEzMa6+zbgMuBPwDvA\nI+7+VhbrHG1mk9evX1+YokVEpFlaxTEodz+vgelPAk/muM5pwLSqqipdelxEJIZaRQ9KRERKjwJK\nRERiSQElIiKxpIASEZFYKtmA0ig+EZF4K9mAcvdp7n5xRUVFsUsREZE0SjagREQk3hRQIiISSwoo\nERGJJQWUiIjEUskGlEbxiYjEW8kGlEbxiYjEW8kGlIiIxJsCSkREYkkBJSIisaSAEhGRWFJAiYhI\nLJVsQGmYuYhIvJVsQGmYuYhIvJVsQImISLwpoEREJJYUUCIiEksKKBERiSUFlIiIxJICSkREYqlk\nA0rnQYmIxFvJBpTOgxIRibeSDSgREYk3BZSIiMSSAkpERGJJASUiIrGkgBIRkVhSQImISCwpoERE\nJJYUUCIiEkslG1C6koSISLyVbEDpShIiIvGWUUCZWXczG2VmR5uZpczrYmbX5WMZERGRhCYDyswO\nAt4B/gi8BLxqZnsnNekKXN/cZURERJJl0oP6CfAyUAH0BxYAs81sSJ6XERERqdM+gzYjgOPdvQao\nAc4xs9uB583seCDdKINclhEREamTSUDtBnjyBHe/Ijqu9AJwXp6WERERqZNJQP0TqALeTp7o7uPM\nrB3hOFM+lhEREamTyTGox2igx+PulwO/BSxlVi7LiIiI1GkyoNz9J+4+spH533L3ds1dRkREJJlC\nQkREYkkBJSIisZTJIAkAzGxgHrf7ibtvyOP6RESkjck4oIBFpAwdb4YbgBvztK6cmNloYPTgwYOL\nWYaIiDQgm4C6kfwF1At5Wk/O3H0aMK2qquqiYtciIiL1ZRxQ7j6hgHWIiIjsQoMkREQklhRQIiIS\nSwooERGJJQ0zFxGRWCrZYeYiIhJvJTvMXERE4k3DzEVEJJY0SEJERGIpq4Ays8nRr+KKiIgUVLY9\nqG8Aj5lZp0IUIyIikpBtQD0KnA78xcz2aKqxme1mZt/JqTIRESlp2QbUvwG/BI4GXjKzvdI1MrNO\nZvZdYCFwR/NKFBGRUpRVQHlwGfAj4ADgr2Z2cGK+mXU2sysIwXQb0A24PY/1iohIicjmPKg67n6z\nmS0D7gVmmdm/AwcCVwKVwL8IAXWru6/OV7EiIlI6cgooAHf/tZltAh4GpkWTPyUE08/cfU0e6hMR\nkRKV03lQZra7mV0F/CIxKbq/xd2vbg3hZGajzWzy+vXri12KiIikke15UF3M7HuE6/L9FOgM3AJ8\nAfgIuNHMrst3kYXg7tPc/eKKiopilyIiImlku4tvEdAD2AjcDNzu7usAzOxo4GngejPrC3zL3fN1\n7T4RESkx2e7i6wD8GBjk7j9KhBOAuy8DjgFeAi4Bfm9mHfNWqYiIlJRse1B7u3uDB23cfb2ZnQxM\nBc4C/mRmZ+i3n0REJFvZngfV5IgCd98MnA1MBo4FZuVWmoiIlLKch5k3xt13AJea2Urgh4XYhoiI\ntG0FCagEd78+OqFXREQkKxkHlJkNzHEbT6VZ9hMdlxIRkcZk04NaRP5+8v0Gwk/Ii4iIpJVNQN1I\n/gLqhTytR0RE2qiMA8rdJxSwDhERkV3kdC0+ERGRQlNAiYhILCmgREQkllpimHk6GmYuIiKN0jBz\nERGJJQ0zFxGRWNIwcxERiSUNkhARkVhSQImISCwpoAqhpgaWLCl2FSIirVrOP7dhZgsyaLYD2AC8\nAzzq7v+b6/Zi7/nnYdy48O/582HgQDjxRLjpJhgxoqiliYi0Rs35Pah20fJ7Ro+3AWuBnknrXQb0\nBj4DnGtmTwJfdvftzdhu/EybBmedBdu27ZzmDn/5C7z4Ijz5JJxwQvHqExFphZqzi+9QYCnwIvAF\noJO79wM6AcdE05cA/YFhwNPAacDlzSk4drZuhYsugm3bKEtMSp6/eXOY7/kaoS8iUhqaE1A3AxXA\nie7+1+hn3nH3He4+GzgJ6A7c7O7vAf9GCLSvNrPmeHn8cVi5EghJDPA+sCW5zYIF8Oc/t3BhIiKt\nW3MC6kzgCXfflm6mu28BpgFnRY8/BWYCQ5uxzfiZP7/un3sSupUbgK8ATwEric5uTmonIiJNa84x\nqJ5AxybadIjaJaxo5jbjp3PnXR5OAk4Gpkc3gM5A/5/+lP7TptGrVy969OhR77bHHntQXl6+y61L\nly60a6eBliJSmsxzPDZiZm8B3YAD3X1jmvndgLeAje5+YDRtCnCSu/dPbZ8vZtYF+CVhL9vz7v4/\njbWvqqryOXPm5L7Bf/wDDjhgl0mLgXsJ13OaD6zPcdVmRteuXesFV3l5Od26daO8vJyuXbvucuvS\npUu9acnTd9ttN8ws9+crItJMZjbX3auaatec3sxk4A7gFTO7GZhN2KPVhzBo4geEvV5XRAUZcBww\nL9sNmdkDQDWwyt0PTpp+KnAnUAb8yt1vIexS/IO7TzOzh4FGA6rZ9t8fRo8OI/kiAwkH6BI2nHMO\nSydMYOnSpXz88cf1buvWrWPdunVs3LiRDRs2sHHjRjZu3Minn35a9+98KSsrazLMkm+JEEy9T/53\np06dFHoiknc5B5S732lmw4BLgQfTNDFgsrvfGT3uDUwFnslhc1OAXyRvx8zKgLsJgzGWAK+a2RPA\nAODNqFnLDGd/8EGorobZs+vPO+kkuj3wAN26dOGAlJ5WU7Zv386mTZt2Ca3ELTFt06ZN1NTUsGnT\nprpbY4+3bNnC+vXrWb8+135dfYnQSxdejQVbcm+woqKi7taxY1N7jkWkFDTreJC7/x8z+x1wIeFc\npwrCGIHXgQfdfVZS25XAtTluZ5aZDUqZfBTwvrsvADCzh4AzCGE1gNBTa5kDON27w6xZ4XynBx8M\no/r694cLL4STToIcexdlZWV1H9r5smXLFmpqahoMsY0bN1JTU1MXfqn36aZt3rw5r6HXqVOneqFV\nUVFRb1pjbXbffXf16kRauWYPWHD3l4CX8lBLtvoDHyU9XgJ8FrgL+IWZjSKMIqzHzC4GLgYYODBP\nv8PYrl3oRVVX52d9BdKxY0c6duzIHnvskbd1JkKvoVBrLOQ2bNhQF26Jf9fW1lJbW8uqVatyrik5\n3NMFWUVFBT179qRnz5706NGj7t89e/ake/futG/ftsbyiLRGeXsXRoMiKoD1LfRruem+Hru71wD/\n2diC7j6ZcAyNqqoqnUHbTPkMPXfn008/3SWwUgOssWnJIZc4xpeL7t27pw2v5Mep88rLy9VrE8mj\nZgVUdBzoKuAbwD5J0xcCvwImNnSeVB4sAfZKejyAcGklacXMjC5dutClSxf23HPPphdoQOJYW0Oh\ntm7dOj7++GPWrl3L2rVrd/n3J598Unf74IMPMt5m+/bt60KrsWDr2bMnffr0oV+/flRUVCjURBrQ\nnIvFdiRcvuhYwrmoHwHLgX7AIMJAtlPN7OTopN18exUYYmb7EK5QcS5wfgG2I61Qx44dqayspLKy\nMutlt2/fzieffFIvuFIfp86rqalh1apVWe2a7NSpE/369Wvwtueee9KvXz969uypc+Kk5DSnB3UF\nYdj4dGB8dDkjAMxsP+A2YHTU7pZmbAczmxptq5eZLQGud/f7zewy4E+EYeYPuPtbWaxzNDB68ODB\nzSlN2qCysrK6nk42Nm/eXBdajQXZ2rVrWbFiBcuXL2fTpk0sXLiQhQsXNrru9u3b07dv30bDrF+/\nfvTp00fHz6TNaM6Jun+P/vmZxHX4Uua3I4ykM3c/JPcSC6vZJ+qKNMOmTZtYvnw5y5YtY/ny5Q3e\n1q1bl9H6zIzKysome2X9+/enQ4cOBX52Ium1xIm6g4H/ThdOEC4aa2ZPAd9uxjZE2rSuXbsyZMgQ\nhgwZ0mi72traul5Xulsi4FavXl23m/GNN95ocH1lZWXsvffe7LfffmlvXbp0yfdTFclacwJqC9C1\niTZdSPn1CRHJXqdOnRg0aBCDBg1qtN22bdtYuXJlo72xZcuWsXTpUhYsWMCCBQt45pn658737du3\nwfDq1auXBnZIi2jOLr5ZhN95OtjdV6eZ34twKbp33f2LzaqygLSLT0pRbW0tCxcu5IMPPtjl9v77\n77Nw4UK2bm34e2V5eXldWA0ePHiX8BowYABlZWUNLisCme/ia05AnQM8BHwI/Bh4jjCKry9hQMMP\nCaP5znP3R3LaSAtQQInsavv27SxZsqReeCVuGzY0fJpjx44dGTRoUNqe17777kunTp1a8JlIXBU8\noKKN/BdwDdFPHqXOBn7m7tfkvIECShrFd9F7773XZHsRCSdSr1mzpsHwWrFiRYPLmhn9+/dnyJAh\nHHHEEQwfPpzhw4czePBgDaEvMS0SUNGGRgBjgcOJriRBuBbfA+7+crNW3gLUgxLJn5qaGhYsWMD7\n779fL7w+/PBDtm+vf/3mbt26cfjhh9cF1vDhwxkyZIhCqw1rsYBq7RRQIi1j69atLF68mLfffpvX\nXnuNuXPnMnfuXJYtq38BmPLy8nqhNXToUIVWG6GAypACSqS4li9fXhdWjYVW165d60KrqqpKodWK\n5T2goh8NzIW7+9gcly04BZRI/KxYsaJeaC1durReu+TQSu5paSRhvBUioNKekJsBd/fY/rUooERa\nh5UrV9YLrSVLltRr16VLl3qhNWzYMIVWjBQioPbOtRh3/zDXZQtFo/hEWr9Vq1bVC62PPvqoXrsu\nXbpw1FFHMWrUKKqrqxk6dKhONi4iHYPKkHpQIm3LqlWrdhmEMXfuXBYvXrxLm8GDB1NdXU11dTXH\nHHMMHTt2LFK1pUkBlSEFlEjbt3r1ap599lmmT5/Ok08+ucsPWZaXl3PKKadQXV3NyJEj6d27dxEr\nLQ0KqAwpoERKy7Zt23jllVeYPn0606dPZ/78+XXzzIzPfvazdb2rQw89VLsCC0ABlSEFlEhpW7Ro\nETNmzGD69Ok8++yzbNmy8/dVBwwYUHfc6oQTTmD33XcvYqVthwIqQwooEUnYtGkTM2fOrOtdJV+6\nqVOnTpx44olUV1czatQo9tprryJW2ropoDKkgBKRdHbs2MHrr79eF1apnxOHHXZY3a7AI488UsPY\ns6CAaoKGmYtINpYvX86TTz7JjBkz+POf/0xNTU3dvMrKSk477TSqq6s5+eST6datWxErjT8FVIbU\ngxKRbG3evJkXXniB6dOnM23aNBYtWlQ3r3379px22mmMHz+eY445RoMs0lBAZUgBJSLN4e688847\ndbsCZ8+ezY4d4cI7VVVVjB8/nrPPPpv27ZvzA+ZtiwIqQwooEcmnlStXcs8993D33XezZs0aAAYO\nHMh3v/tdxo4dq91/ZB5QugywiEge9enThwkTJrB48WImTZrE0KFDWbx4MVdccQV77bUXV111VdrL\nMUl9CigRkQLo3Lkzl1xyCe+88w5PPPEExx57LBs2bGDixInsu+++fPWrX+W1114rdpmxpoASESmg\ndu3aMXr0aJ5//nleffVVzj33XNyd3/3udwwfPpwTTjiBGTNm1B23kp0UUCIiLaSqqoqpU6fywQcf\nMG7cOMrLy3nuueeorq7moIMO4r777qO2trbYZcZGyQ6S0HlQIlJs69ev57777uPOO++s+22ryspK\nLrvsMr75zW9SWVlZ5AoLQ6P4MqRRfCJSbFu3buX3v/89EydO5PXXXwfCpZUuuOACxo0bx7Bhw4pc\nYX5pFJ+ISCvRoUMHzj//fObOncuzzz7LqFGjqK2t5d5772X//ffnjDPOYNasWZRah0IBJSISE2bG\n8ccfz/Tp03n77be56KKL2G233epGAR511FFMmzat2GW2GAWUiEgMHXDAAUyePJnFixdz3XXX0atX\nL+bMmcPpp5/OhRdeyIYNG4pdYsEpoEREYqx3797ccMMNLF68mNtuu41OnTrxm9/8hsMOO4wXX3yx\n2OUVlAJKRKQV6Ny5M1dccQWvvfYaRxxxBIsWLeLYY4/l2muv3eVHFtsSBZSISCtywAEH8PLLL/P9\n738fM+OWW25hxIgRvP3228UuLe8UUCIirUzHjh25+eabmTVrFvvssw+vv/46w4cP56677mpTV6RQ\nQImItFKf//zneeONN/j6179ObW0tl19+OaeeeipLly4tdml5UbIBZWajzWzy+vXri12KiEjOysvL\nuf/++3n00Ufp1asXzzzzDIcccgiPPPJIsUtrtpINKHef5u4XV1RUFLsUEZFmO/PMM3nzzTcZOXIk\n69atY8yYMXzta1+jNX8JL9mAEhFpa/r27cuMGTO455576Ny5M7/97W859NBDeeGFF4pdWk4UUCIi\nbYiZcemllzJv3jyOPPJIFi9ezPHHH89VV13F5s2bi11eVhRQIiJt0NChQ5k9ezbXXXcd7dq1Y+LE\niRx11FG8+eabxS4tYwooEZE2qkOHDtxwww289NJL7Lfffvz973+nqqqK22+/vVUMR1dAiYi0cSNG\njGDevHlcfPHFbNmyhfHjx/OlL32p7jeo4koBJSJSArp27cq9997LE088QWVlJc899xzHHHMMK1as\nKHZpDVJAiYiUkNGjRzN//nyOPPJIFi1axMiRI2N7ZXQFlIhIienduzczZsxgyJAhzJs3j7POOiuW\nF5xVQImIlKDKykqefvpp+vTpw8yZM7nwwgtjN3BCASUiUqL23XdfnnrqKbp27crUqVO5+uqri13S\nLhRQIiIl7PDDD+fRRx+lQ4cO3Hbbbdx+++3FLqlOyQaULhYrIhKcdNJJTJkyBYDx48czderU4hYU\nKdmA0sViRUR2Ov/885k4cSIAF1xwATNnzixyRSUcUCIisqvx48czbtw4tm7dyplnnsm8efOKWo8C\nSkRE6kycOJFzzz2XjRs3MnLkSBYuXFi0WhRQIiJSp127dkyZMoUTTjiBFStWcMopp7BmzZri1FKU\nrYqISGzttttuPPbYYxx22GEKL/2NAAAKGUlEQVS89957VFdXU1NT0+J1KKBERKSebt268dRTTzFo\n0CBeeeUVxowZw7Zt21q0BgWUiIik1a9fP55++ml69uzJjBkzuOSSS3D3Ftu+AkpERBo0bNgwpk+f\nTufOnXnggQe4/vrrW2zbCigREWnUiBEjeOSRRygrK+Omm25i0qRJLbJdBZSIiDSpurqae++9F4Bv\nfetbzJgxo+DbbF/wLYiISJswduxYli1bxuOPP05VVVXBt2ctecArjqqqqnzOnDnFLkNEpFVwd2pr\na+ncuXPO6zCzue7eZMJpF5+IiGTMzJoVTtlQQImISCwpoEREJJYUUCIiEksKKBERiaWSDSj9oq6I\nSLyVbEDpF3VFROKt5M+DMrPVwIeNNKkAmupm9QKK84MpxZHJa9JSWqKWfG6juevKdflslsu0rd4b\n9cXpvQGFryfX9e/t7pVNtnJ33Rq5AZMzaDOn2HXG7TVpS7XkcxvNXVeuy2ezXKZt9d4o7N9Ka6in\n0Osv2V18WZhW7AJiKE6vSUvUks9tNHdduS6fzXKZto3T30FcxO01KXQ9BV1/ye/iywczm+MZXLZD\npNTovSHNoR5UfkwudgEiMaX3huRMPSgREYkl9aBERCSWFFAiIhJLCigREYklBVSBmdmXzew+M/uj\nmZ1c7HpE4sLM9jWz+83sD8WuReJJAdUIM3vAzFaZ2fyU6aea2T/N7H0zu6axdbj74+5+EXAhMKaA\n5Yq0mDy9Nxa4+9jCViqtmUbxNcLMvghsAh5094OjaWXAu8BJwBLgVeA8oAz4Scoqvu7uq6LlbgP+\nx91fa6HyRQomz++NP7j72S1Vu7Qe7YtdQJy5+ywzG5Qy+SjgfXdfAGBmDwFnuPtPgOrUdZiZAbcA\nTymcpK3Ix3tDpCnaxZe9/sBHSY+XRNMa8m3gS8DZZnZpIQsTKbKs3htm1tPMJgGHm9m1hS5OWh/1\noLJnaaY1uJ/U3e8C7ipcOSKxke17Yy2gL23SIPWgsrcE2Cvp8QBgWZFqEYkTvTckrxRQ2XsVGGJm\n+5hZR+Bc4Iki1yQSB3pvSF4poBphZlOBl4FhZrbEzMa6+zbgMuBPwDvAI+7+VjHrFGlpem9IS9Aw\ncxERiSX1oEREJJYUUCIiEksKKBERiSUFlIiIxJICSkREYkkBJSIisaSAEhGRWFJAiYhILCmgRCTW\nzOxaM3vVzDaY2Wozm2ZmBxe7Lik8BZSIxN1xwC+Bo4ETgG3AX8ysRzGLksJTQEmrZmaDzMzNbEoR\ntj0h2vagbOaVuoZem4amu/sp7v5rd5/v7m8CXwMqgc+3UMlSJAookSIys67RxVbdzB5vpF13M/vE\nzNaaWXma+e3M7A/Ren5V2KrrtvmjprZnZsPNbKOZrTKzIXnadDnhs2tdntYnMaWAEimuCYRfnXXg\niIYaufsnhB++7EG4Yniqu4CvANOBS/JeZXp3EkLia2Y2IHWmmQ0GniT8kOEod38vj9udR7iaurRh\nCiiRIjGzA4HvAE8Dfwf2MrNejSxyB7ABGG9mXZPW8wPgW8DfgDHuvr1wVe/k7huimjoCVybPM7M+\nhJ/d6AGc7e6v5mObZnY78AXgKy31PKV4FFDSZpnZOWY2y8zWm9m/zOzNaETYbmnampldbmZvm1mt\nmS01s1+YWYWZLTKzRQUo8ReE3sUVhIACOLyhxu6+DvhvoCchkDCz/wR+DPwTqHb3TwtQZ2MSvaiL\nEuEa7YJ8CtgXGOvuT+djQ2Z2B3AecIK7L8jHOiXeFFDSJpnZfwEPAwcAv2NnGPwX8Ccz65CyyN3A\nz4EKYDIwFTgZeAZIbZuP+s4Hjgd+6e7vAG9GsxrczRe5HdgEXGlm50S1LgdOdfe1+a6zKVEv6ufA\n7sDl0S/pPkYI2u+5+4P52I6Z3QmcTwinf+RjnRJ/Cihpc8zsc8C1wEfAIe7+TXe/CvgM4RjNscBV\nSe2PAb4JvAsc5O7fcfcrgYOBWmDPPNdXDtwKrCUcg4KdPahGA8rdPyaEbS9CAH8KjHT3RfmsMUs/\nJ/SiLiME+4nAne7+s3ys3MzuBv6T0HtaZ2Z9o1vXJhaVVk4BJW3R16P7H7v7isTE6CfJxwM7gG8k\ntb8gur85GoyQaL+FEHT5dgMh9K6PdtvBzh5Ug7v4kkxP+vdX3f2NfBaXraReVHfgLEJwjsvjJv4P\nYeTeTEJvMXG7srGFpPVTQElblOiFPJs6w93fBZYA+5hZ92hyIhReSrOuvxFODM0LMzsI+DbwFjAp\nqa5lhB7VYDPr1sjyewL/kzTpwCy3vygaGp7p7bcZrvqh6H4J8B/u7tnU1Rh3twZuE/K1DYmn9sUu\nQKQAKqL75Q3MXw4MjNp9ktR+ZWpDd99uZvk8tnM34X03Ls0otDcJV034DDArdcEoUJ8G9gauA75H\nOBZ1t7vXZLj9Dwi7LTO1LMN2B0X3r0Q9T5FmU0BJW7Q+uu9L+EBO1S+l3Ybovg+wy+gwMysjjJpb\n2tyizOzfCce/AP5sZg01PYKUgDKzTsAfgUOAG939pqindSXh+NnETGpw9xNzKD0Tw6P7OQVav5Qg\n7eKTtuj16P641BnRyaMDgIVJx5sS7b+QZl0jyMMXuShMbgW2Ag8A96e5zYyaH56ybBlhJOIXgcnu\nfn0062eEQRJXmdnuza2xmaqi+7wHlC4WW7oUUNIWPRDd/9DMKhMTow/6iYS/+/uT2ieGQv/AzCqS\n2nckDEvPhxsJPbo73H2su38j9QZcE7VNHcl3N3Am8DhhwAAA7r6acBHV3sCleaozV4ke1NwCrPs4\ndLHYkqSAkjbH3f9K6F0MAuab2d1m9jPC5XHOIAyGuDWp/QuE84mGAm+Z2Z1mNpFwTKgz4TjMjlzr\nMbNDCCfWfkQIqobMB7YDB5hZ52jZGwiXLnoROC/NcatbCb2oqxPLtDQz25sw7H1B0qjEvNHFYkuX\nAkraJHf/HuG8mfeA/yBcUqgd8EPgpDQH8r9JuKLDJkJv5HzgL8BJQDd2HqfKRWJgxOWNDWZw99qo\n3jLgUDO7lDAYYj5wejQ/dZlVwD2E42ctdQ2+VC19/EkXiy0RGiQhrVp0gmra0Qbu/hA7hz83tZ4d\nhOvK3ZE8PboCd1fgnWbU+MUs2h6Q9PAVkoaiN7LMlRTxnCB3f5QG/g8KRBeLLREKKBHAzPoCq6Kg\nSkzbnXACKoTL90iRJV0s9gu6WGzbp4ASCb4LnGdmzxPOk+pLuGTPAMKFT39fvNIE6i4Wey5wvC4W\nWxoUUCLBM8BhhAvE9iCMFHuX8DtLP8/nlREke9HFYs8FjtPFYkuHAkoEcPeZ7DwPKVPPR/efZDmv\n1D0f3ae+NmmnRxeL/RrwZaKLxUazNrn7pgLVKDFg+mIoInFmZg19SN2g6/G1bQooERGJJZ0HJSIi\nsaSAEhGRWFJAiYhILCmgREQklhRQIiISSwooERGJJQWUiIjEkgJKRERi6f8DfPz2nFJxdJcAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbdd98b400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAEHCAYAAADxpHypAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xu4HEWd//H3J4EEcUXuiAQFJaLg\nBSUCuoAIgqgroIKCKEH5ibLiendBV0TQFXRdb4to1Ai4uMCCSFZhI4hBcQETIAIBAyFECOEidxRI\ncnK+vz+qhkwm3XNm5sycuZzP63n6mZm+VNf0XL5d1dVVigjMzMyseRO6nQEzM7N+5SBqZmbWIgdR\nMzOzFjmImpmZtchB1MzMrEUOomZmZi0a10FU0pGSQtKR42nfg07SHElN3buVP4s5HcqSDThJe+Xv\n0IndzouNrYEKopImSvqApCskPSRppaT7Jd0g6YeSDuh2HseapEmSPirpj5Iez8dlnqSTJK3fYppR\nMy2X9BdJ1+Xj/CZJE9v9XvqZpDMKjtsTkm6W9HVJm3U7j90kaZt8TM7odl6K9Hr+rHvW6XYG2iX/\naf8C2B94BPglsBTYGHgh8G7gxcCsqs0uBK4G7hnTzI6tmcDhwL3AmcBKYEfguLxsySjS/mJ+nAhs\nmNN9L3AUME/S4RFx6yjSH0QXAfPz8y2ANwOfAN4haeeIeLBrObPR+APwEuCBbmfExtbABFHgMFIA\n/SPwuoh4tHphLnXtWj0vr7PGeoNE0rNJJw+PAq+IiPurlm0M/G006UfEiQX73AL4DnAIcJmkadX7\nNX4eEWdUXkhaj3Qi9wrgWFafmFgfiYgngD91Ox829gapOve1+fGM2gAK6UseEb+pnld2XVLSkjyt\nL+lrku7MVZaLJP2zJNWmr+SjuXruKUl3S/oPSc+upNfoG5E0JW+7OO/3QUmzJL260TSyFaSS5921\ngSwiHoqI5U2mN6KIuA84FJgDbA18tnYdSVMlnZWP0QpJy/LrqQXrVqpBtylYVvc6lKTJkr4k6Y58\nHG+X9AVJkxp9P5LWkfSPkq6W9Fiugr1e0rGSRv37iYingLPzy8LPt9nvQ76s8SFJv5f0qKQn83f3\nh7XHOH8/vyJpYf7ePixptqQ3FKT79PGWtJOkX0p6JB+TKyS9tmCbZ0n6vKSb8vF7PH8O50raOa9z\nInBH3mR6TZX3kQX73iXv+6Hq74bqXNce4Xu0S87P3fn43iPpV5Le2Ur+CtJv5vt+Yk5nL0kHS/pD\nPr4PSTpH0lZF78+6Z5BKopVqsBe1Kb11gV8BzwUuAYaAg4BTgPVYu8RwGnAMsAyYQQpgBwC75LRW\nNrJTSa/K+90YmA38DNg07/tKSW+LiIsbSSsinpR0GvBxparVs0fcqA0iYljSl4C9gMMkfTxyJ835\nj/8y4FmkqvWbSdXshwMHStonIua1KSvnkQLT+aTjfyBwIjBN0gExQsfRktYF/gd4I7AQ+CnwFPB6\nUml7V1L19WhVTsrW+o40+33IJwi/BN4A3JXz/BiwDfA24ErgtrzuhsDvgR2AucA3c9rvBH4l6ZiI\n+H5BfqcBnwGuAn4IPA94B/BrSTtFxMKcvoD/JZ3gVtYdIp1c7QX8DriWdMK1IfBRUk3Sz6v2NZ81\nvQY4Pr+PmTm/Kwry2BBJHwBOB1aRvo+3AZvn9/iPpO9QM/mrTb/V7/s/kv4/ZgFXkL5r7wJekY9x\n20+ArUURMRAT8ErSj2kY+AnwduD5I2xzJBDAkTXzl+T5FwPPqJq/Oel66yPAulXz98jrLwQ2rJo/\nCfhtXrZkpH2TTmoWkf6oX1ez/nOBu0nXbyc3eEwmAt/N+1kJHN6mYx3pq1N3ncl5nwFsW4mhwC15\n3uE1678rz/8TMKFq/hl5/jYF+9grLzuxZv6cPP9WYKOq+euR/swDeG/Be5pTM+/EPP87wMSa4/qj\nvOzABo9Z5X3UfteeAdyQl32yZlnT3wfgX3Nas2q/J/kz2azq9ffzut8HVDV/KukSwPLq4151vIve\nxwfz/O9WzXtZnndhwfGYUPPZbJPXPaPk+FXv+4N1vpdzSpat9T0inTysBB4CdizYZkoL+Tuxal4r\n3/fKd+4x4GU12/w0L3tnO37HntozDUx1bkRcD7wHuC8/XgAsyVVfF0p6awvJ/lNEPFm1j/tJDUOe\nDWxftd70/PjliHikav0VpLPmRr2F1AjqOxFxRfWCiFgGfBV4DrBPg+l9n1Q6/iSp0dVPJP1T7UqS\nZuQqpFc1kde6Ip0pV2oHKi1PX0s6C78qakrFEXEuqXSxPbB7m7JxckQ8XLWPp1j9eby/3oa5qvZY\nUoOsj0fEqqp0VpGOaZBKFM04KFfZnSjpu6QTr5eRTrZOr1m3qe+DUuO6fwSeBD4UNaWViFgeEX/J\n665L+p38FTg+8r90Xu824Nukk8AjCt7D76Pqum42k1TK3KVg/SdrZ0TEcPVn04T5UVw6bsUxpBOV\nkyNiQe3CiFg6yvRH833/dkTcWDPvB/mx6BhblwxSdS4RcZ6kC0nVbbuTSqe7k6q+DpJ0FukMupF7\nCB+NiEUF8+/KjxtVzXtlfryyYP2rSX8ujXhNfnx+0bUVUgkBUivAulW6+ZrWUcA5EfHvkr4D/Dfw\nLUmbRcTna9J9glTV1E6VasrK8a4E6ctL1r+c1Z/bb9uw/ysK5v2O9Hm8smBZtRcBm5Cq9/5Fa18G\nhxQcXtJkng7MU7VLgbdERG11brPfhxeTTvCuyUG2nhcD65MC4kMFyy8H/oXi47RW9WNErJR0H2v+\nLm4mVXceJun5pBPQK4F5+QSzFX9ocbsiu+XHS9qYZrXRfN+LqniL/nusywYqiEL6MZOuIf0Knj47\nfwfpTPkI0m0tPy9NYLVHSuZXAmL1fZDPzo/3FeRnlaRGb1vYJD8eMsJ6f9dAWpUS0vdzPlbmhhLn\nk4LCpsCH8z5fC/wsl9TaQqnV6cb55V/yY+U4ld1SVJm/YZuyUe/z2HyEbSufxVTgC3XWa+SzqPa+\niDgjfy9fAJxMqto7Hfh/JXlo9PtQOW53N5CP0XwW9X4bT/8u8rHeGzgBOBg4NS96XNKZpBLwXxvI\na7V7m1y/nmaOVyvafYyL/nusywamOrdMRKyKiPOAb+RZe3dgN4/lxy1qF+Q/y01q55eotCo+MCJU\nZ2rkNojn5sen71vLZ/8Hk6p2PwScA3ye1PDpqw3msVG7k07S7ouIJXle5f09p2SbLWvWg3SNG4pP\n+EYKtvU+j8fWXn0NlTxcOMJnse0I6RTK38vbSLcgXQMcpbU7A2n2+1D5422kBWcrn0XTIuLhiPh4\nRGxNOiH5f6TrgMeydvV1Q0mOsKysYFAvUHWqxeuYHGPrroEPolUez4+F9XKjdH1+LLq2sRuNl/iv\nzo97jDpHqztR2Kt6Zg6k7yC14DwE+AjwvXxNuS3y9cTP5Zc/rVpU2ccaeapSmX9d1bzKdbOtC9af\nNkJWXlcwbw/S5zHS+/0T6U92t3z9sCMiYpjU6hPgq1qzp6dmvw+VPL9c0nNHWHchqQp/J0lF1YOv\nz4/XFSxrSUQsiogfkT6Xv7JmtXblmvNoSlkPU/A9ycd0p4L1K8f3TQ2k3Ur+Wvm+W58ZmCAq6TBJ\n+6rg3j1JzwE+kF+241pbrbPy4+eUOjio7HcSqbVkoy4Cbgc+LOnNRStIeo0a665vJunM/GRJe1Uv\nyIH0/KpZWxQdt1ZI2pxUwt0LuJM13//vSX/eu0s6uGa7g4E9SS1qq68tV66BfaBm/ZexOviU+Xx1\ngMhVzF/JL39cb8OIGCK1yt0S+LakZ9SuI2lLSTuMkIcRRcQ1pNqB7VmzIU9T34fc4Om7pBa/35M0\nuWbdScrdC+bvwNmkquCTatZ7IfBPpJarP2n1fUnaVtKOBYs2IrUUrm5w9DDp+/q8VvdH+q48T9J+\nNfP/BXh+wfqnk6pIP1/0OUqaMsr8tfJ9tz4zSNdEdyX9qd4r6UpW3xy9LamV4zNIf0rnF2/euoi4\nQtIM4GhggaQLSH9AbyVV1SxjdbVkvXRWSno76X7AX0r6P1LDjCdIZ9ivJl1H2zLPq5fWNZI+BXwN\n+I2ky0j3uK1LKmW8jPQjHybdDnQ66TaFhlU1dpnA6m7/die16vwDqVl/dXVySJpOakhzrqSLSKWn\n7UmNvx4Hjsils4qLSI17Dst/ateQ/sgOzMveWSeLt5A+j+r7RF9IKoU3EhxOJvUk9CHgrZIuJ10/\n25xUNfn3pBJ3OxpknUD6nn5B0tkRsaLF78MXSb+FtwK3SvoF6bhuDewHfJp0uwekrh/3AI7N9zP+\nhtX3iT4LODYiKr+jVrwCuFDStcBNpN/BZqTPYV1WXyMlIv4q6RpgD0lnk4LLKmBWRNzQ4P7+jXRP\n70WSziXduvJa0n/AHNaulblZ0j8C3wOuz9/H20jV/dNIx+31reavxe+79Zt23CfTCxPpT+LDpIZD\nC0nXvFaQLt5fTGrOP6FmmyMpv090Scl+Tszb7FUzfwLwcdKPZDnpD+M0UuOCx0lN80fcd162OalT\nh5tIf45/Jf24z8/vY50mjstrSa1y7yUFkgdIN8AflvO8MekPIUhN/RtJM2qm5Tnda0nN8PevPdY1\n229PCmL35DzdA/wnsH2dz7byp/gkqWOAtzPyfaKTgS+RTqiWA4tJjYTWus+WknsMSdX/7wV+nfe/\nghRIryT1xrR1g8fsjLLPu2qdC/I6HxnN94F0cnws6UTmr6TuHW8jdQKyXc26G5KC2W35GD1C+tPf\nryB/hce77HcDTCHVRPw+f/+Wk/qzvgR4U8H225E6t3iQdHL39PEaad9VaRxAatn6VE7nHFIptHL8\ntynY5jX52N+fP99lpN/Iwe3IH0183yn5f8nLtqHOvaqeujMpfzjWIUpde91KutXksG7nx8zM2mdg\nrol2m6Tn1F5XzNeqvplfXjj2uTIzs05yEG2fjwF3SDpT0ilK4w4uJA11dQmpStXMrO9I2l9pkIJF\nko4rWP4JpcE3bpD069y5RmXZdEm35Wl61fydJd2Y0/y2Sno06XWuzm0TSfsAnyI1pd+Y1OrvVtIt\nHt+MtXujMTPrefkWoVuBfUnXtOcCh0XEzVXrvJ7UU9YTko4hXdN9l9KQi/NIDbWC1G5i54h4WNIf\nSI1Brya1W/l2RHSq96iOGaTWuV0VEb8mNT4xMxskuwCLImIxgKRzSC2snw6iseYwk1eTGrxBai19\naeSuJSVdCuyvNGTdBhFxVZ5/FqnFct8FUVfnmplZPVuxut9eSKXRer08HcXqYFi27Vb5eaNp9qyB\nKolO0uRYj2d2OxtmZh31FH9jRSx/+hriG1//zHjwoVX1Nqnr2huWLyDdFlQxIyJm5OdF1yoLrwNK\neg+p6rbSW1jZtg2n2esGKoiuxzPZVY2OEmZm1p+uiTWvHD3w0CqumT2lZO2Rrbvl7U9FRFk3mktZ\nszvFKaR7adeQR476HGns2+VV2+5Vs+2cPH9KzfyRRh7qSa7ONTPrc0GwMla1PI1gLjA1d+M4CTiU\nNOj70yS9kjRi1AGRxl2umA3sJ2mj3AXnfsDsiLiHNJrPbrlV7hGkHsj6zkCVRM3MxqvhkXsWbUlE\nDEk6lhQQJwIzI2KBpJNIY8POInUv+nfAf+c7Ve6MiAMi4iFJJ5MCMcBJsXr82mNIPUk9g3QNte8a\nFYGDqJlZ30sl0c51wRsRF5NuQ6med0LV8zfU2XYmaUCM2vnzgJe2MZtd4SBqZtbnAljVn+1y+p6D\nqJnZABh2EO0KB1Ezsz4XwEr3PtcVDqJmZn0uCFfndomDqJlZvwtY5RjaFQ6iZmZ9LhArCzsBsk5z\nEDUz63MBDLsk2hUOomZmA2CVS6Jd4SBqZtbnUutc9+LaDQ6iZmZ9LnW24JJoN/R0EJW0IfBDUtdQ\nAby/MoirmZklgVjl8US6oqeDKPAt4H8j4uA8esD63c6QmVmvcXVu9/RsEJW0AbAncCRARKwAVnQz\nT2ZmvUmschDtip4NosALgL8AP5b0CuBa4KMR8bfqlSQdDRwNsJ4LqmY2DgUw7Orcrujlo74O8Crg\n9Ih4JfA34LjalSJiRkRMi4hp6zJ5rPNoZtZ1EWJFTGx5stb1chBdCiyNiGvy6/NJQdXMzGoMo5an\nkUjaX9JCSYskrVWYkbSnpOskDUk6uGr+6yXNr5qeknRQXnaGpDuqlu3U1gMyRnq2Ojci7pV0l6Tt\nI2IhsA9wc7fzZWbWawKxIjrzdy5pInAasC+pcDNX0qyIqP4/vpPUfuVTa+Qr4jfATjmdjYFFwK+q\nVvl0RJzfkYyPkZ4NotlHgLNzy9zFwPu6nB8zs57T4WuiuwCLImIxgKRzgAOpKtRExJK8bLhOOgcD\nl0TEE53KaDf0dBCNiPnAtG7nw8ys162KUXW2sKmkeVWvZ0TEjPx8K+CuqmVLgV1b2MehwL/XzPuy\npBOAXwPHRcTyFtLtqp4OomZmNrJArBxdde4DEVFWYCmKzk11dy9pS+BlwOyq2ccD9wKTgBnAPwMn\nNZNuL3AQNTPrc6nbv45V5y4Ftq56PQVY1mQa7wQujIiVlRkRcU9+ulzSj6m5ntoverl1rpmZNSAQ\nq6L1aQRzgamSts3tUw4FZjWZxcOA/6qekUunSBJwEHBTk2n2BJdEzcz6XASjrc6tk3YMSTqWVBU7\nEZgZEQsknQTMi4hZkl4NXAhsBLxV0hcjYkcASduQSrJX1CR9tqTNSNXF84EPdeQNdJiDqJlZ32vs\nfs9WRcTFwMU1806oej6XVM1btO0SUuOk2vl7tzeX3eEgambW5wLcd26XOIiamfW51DrX3fd1g4Po\nWFBJNUs01Urc2smfiQ0YjyfaHQ6iZmZ9LoBhV+d2hYOomVmfc3Vu9ziImpkNgFUdbJ1r5RxEzcz6\nXIRcndslDqJmZn0uwNW5XeIgamsra7laT1mr1l5tBdvt/Zu1lXyfaJc4iJqZ9TmXRLvHQdTMrM8F\nYnh044lai3o+iEqaCMwD7o6If+h2fszMetGwO1voip4PosBHgVuADbqdETOzXhQBK4cdRLuhp4+6\npCnAW4AfdjsvZma9KlXnTmh5stb1ekn0m8BngGd1OyP9bOV+0wrnL3lL8cc/Yag4nVZ+axoqvk4z\nPKm8dey2F60snL/O5dc2nwGzccKdLXRHz56CSPoH4P6IqPvPKeloSfMkzVvJ8jHKnZlZ7wjE0PDE\nlqeRSNpf0kJJiyQdV7B8T0nXSRqSdHDNslWS5udpVtX8bSVdI+k2SedKmtSWgzHGejaIAn8PHCBp\nCXAOsLek/6xdKSJmRMS0iJi2LpPHOo9mZj1hOA/M3cpUT27ceRrwJmAH4DBJO9SsdidwJPDTgiSe\njIid8nRA1fxTgW9ExFTgYeCo1t55d/VsEI2I4yNiSkRsAxwKXB4R7+lytszMek4ErAq1PI1gF2BR\nRCyOiBWkQs2Ba+4/lkTEDcBwI/mVJGBv4Pw860zgoGbec6/o9WuiZmY2gkp17ihsKmle1esZETEj\nP98KuKtq2VJg1ybSXi+nPQScEhE/BzYBHomISguMpXk/facvgmhEzAHmdDkbZmY9a6Rq2RE8EBHF\nLRApTLiZfjOfFxHLJL0AuFzSjcBjo0yzZ/Rsda6ZmTUmDcqtlqcRLAW2rno9BVjWcN4iluXHxaTC\n0CuBB4ANJVUKck2l2Uv6oiTa98ags/Pbv/aa8oVqz/5jYnk6Gi7+Ica6xduozpWTP7953eIFb9mt\ncPYLP3l1eWJm40GMujq3nrnAVEnbAneT2qi8u5ENJW0EPBERyyVtSmow+tWICEm/AQ4mXWOdDlzU\nkdx3mEuiZmZ9Luhc69x83fJYYDap97jzImKBpJMkHQAg6dWSlgKHAN+XtCBv/hJgnqQ/Ar8hXRO9\nOS/7Z+ATkhaRrpH+qL1HZWy4JGpmNgA62QF9RFwMXFwz74Sq53NJVbK12/0f8LKSNBeTWv72NQdR\nM7M+F8CQ+87tCgdRM7M+56HQusdB1Mys3wUMuSP5rnAQbSeVnAm2sXXu7V8vbqEaE8r3Ufbb0qrm\n9j3pkeZ/pCs2LG6GO2Fl+VnzcEnj3LJGxrf/W/ExAXjhp9xy1wZf5RYXG3sOomZmA8BBtDscRM3M\n+lwgVrlhUVc4iJqZDYBRdvtnLXIQNTPrcxGuzu0WB1Ezs77n6tyRSJoNzAeuz48LI0bf6tNBtJ3G\noBVuKyasKD5DfdGM+wrnr7ptcdv2XWadrZ5buuxPn3xe4fyy/najTpehi08t7lP4Bf98VflGZn0o\nXBIdyfXATsARwBbAE3lEmfmsDq43RMRTzSTqIGpm1ud8i8vIIuK4ynNJW5BGk9kpTx8HtgNC0m0R\nsUOj6TqImpn1u4BVDqINi4j7gP/NEwCSnkEKqC9vJi0HUTOzPhe4One0IuJJ4Ko8Naxng6ikrYGz\ngOcAw8CMiPhWd3NlZtaL3HduMyRNBj4GvJg06Ph8YH5E3N5sWr3cnGsI+GREvATYDfiwpIbrqc3M\nxpPhYbU8jUTS/pIWSlok6biC5XtKuk7SkKSDq+bvJOkqSQsk3SDpXVXLzpB0h6T5edqpbQdjZN8D\nPgKsII1r+mPgVkmPSPptMwk1XBKVtD6wO7AjsDmpBuEvwE3A7yPiiWZ2PJKIuAe4Jz9/XNItwFbA\nzXU3NDMbZyI6V50raSJwGrAvqdQ2V9KsqsG1Ae4EjgQ+VbP5E8AREXGbpOcC10qaHRGP5OWfjojz\nO5Lx+t4CHBoRl0t6N/BqUmHtK0BTHW6PGEQlvQn4ELB/Xr/2kwpgSNIlwPci4n9pM0nbkFpSXVOw\n7GjgaID1WL/du27OhJJ7LYaLe3qf8PIXlyZV9nuIdYpvo5mwvLxS4YWfLq7ib7L/+bYauntZ6bLt\nPlG87LZv71o4X0N1OrOfVHy8hvbZuXD+Or++tjQts17WwercXYBFeRBtJJ0DHEhVgSYiluRla9yI\nFhG3Vj1fJul+YDPgEbrrGcBt+fkKYEJEnCnpWaRWug0r/eeVtIekPwC/JI1MPhN4H7AHsAOpRLoH\n8H5SUXgn4GJJcyXt3tz7KSfp74ALgI9FxGO1yyNiRkRMi4hp6zK5Xbs1M+srHazO3Qq4q+r10jyv\nKZJ2ASYB1dcdv5yreb+Rr1OOlcWsfg93A1Py80uAQ5tJqF5JdA7wc9J1yd/VWe/3wJkAkl5Hulg7\nZ4S0GyJpXVIAPTsifjba9MzMBlGg0VbnbippXtXrGRExIz8vSripnmUkbQn8BJgeEZXS6vHAvaTA\nOoN0bfKkpnLduvOA/UhVt3OAo4BLSQXG9ZpJqF6ge1VE/LGZxCLiCuAKSa9oZrsikgT8CLglIv59\ntOmZmQ2s0fed+0BETCtZthTYuur1FKD8ekwNSRuQajT/JSKevt6Y270ALJf0Y9a+ntoxEfHlqpdf\nJV3nfQj4O+D7zaRVGkSbDaDt2rbK3wPvBW6UND/P+2xEXNyGtM3MBkv7eh2tNReYKmlbUtXnocC7\nG9lQ0iTgQuCsiPjvmmVbRsQ9ucB0EKmR6piLiKWSdiQ1Nnqw2RjTs/eJRsSVFFcjmJlZjU61zo2I\nIUnHArOBicDMiFgg6SRgXkTMkvRqUrDcCHirpC9GxI7AO4E9gU0kHZmTPDIi5gNnS9qM9D8/n9SA\ndUzkPnN3j4hH83t8iFTd3LSmgqikZ5LOQKYCm1DQUjcijmolIwOhpBUuKv5y3/q+DcvTKvk9lHUm\nv925j5cm1bkT1LH14m/fXzj/1mO2KN0mSprO3XHAuoXzp/666WyZdV1AQ/d7tpx+Kp1dXDPvhKrn\nc1ndOKd6nf8E/rMkzb3bnM1m7Ahrt0SV9GzgXyPiw40m1Mx9oq8FZgEb11ktSBdozcxsrATl98XZ\n0yRdDPyBdMS2BmrPzNcHPgi0P4gC3yF1v3cg8Luqm2XNzKzL2jgS4yBbAOxFquv7g6THgT+Sh0Ej\ndQN4T+nWBZoJojsAJ0TE/zSzAzMz6zQRHazOHRQR8WkAScuB1wDPZfVwaG8hxcTPNJNmM0H0HmBl\nM4mbmdkYcUm0Gc+MiCHgOuAXo0momQ7ofwi8O/ejaGZmvSL3ndvqNA5NlHSqpFskLZZ0kaRDWkmo\nmZLoV0hF36sknQ4soaD71Yhoqgf88axen69MKD6tHJ5cPD/mdeUWqzG1atEdhfOH1ylvnauyBtPt\nyJBZLxmfwbBV/wb8A/BdYDnpcuVMSYcC78ql1IY0E0SfQbqtZWdSqbSWSBUKLqmamY01V+c24xDg\n7RHxf5UZkr5Auo3nOOBLjSbUTBA9jXTj7M+B3wEPN7GtmZl1koNoM9aj5vaWiLhP0sdJA6p0JIge\nSOqp4gNNbGNmZp0WuHVuc64g9WlwfM38pUD59aECzQRRkfpQNDOzXuOSaDOOA/5P0ibAN4E/kUaT\n+SjpXtKGNdM6dw5p5G8zM+s1odancSYibgFeB7yc1PH9U8BfSZ3rf7SZtJopiX4MuFzSJ4D/iIgV\nzezIzMw6JEDDI69mq0XEDcBukl4EvBR4HLgmIh5rJp1mguhvgGcCXwNOkXQPa9/iEhHxwmYyMC6U\n9Mc1PKn5+pd1H2um8mB8mLi8/Ey67BiX3fpi1p/GZ4myVdWjuETErcCtrabVTBC9E9e6m5n1JpdE\nmzH2o7hExF6NrmtmZmPMRZwRdWIUl56uG5S0v6SFkhZJOq7b+TEz60mVodA61LBopP9iSXtKuk7S\nkKSDa5ZNl3RbnqZXzd9Z0o05zW9LJQMvt1ftKC6PSLpC0jclvR/4OE2O4tJwEJX0BklfqbP8K5Je\n38zOR9jfRFIHD28idcl0mKQd2pW+mdkg0XDrU910G/svvhM4EvhpzbYbA18AdgV2Ab4gaaO8+HTg\naGBqnvYfxdtvSER8OteqrgReDbwHuJQ0oPjxpJ6MOjaKy2eAR+ss3xb4Z1IDpHbYBVgUEYsBJJ1D\n6vDh5jalb2ZmIxvxvzgiluRltSH5jcClEfFQXn4psL+kOcAGEXFVnn8WcBBwSUffyWpdGcXlFcDV\ndZZfk9dpl62Au6peL83z1iDpaEnzJM1byfI27r7z6p4ZKgqnmEDhNK4Nl08xIYon3y5nA6bkL6Oh\nCdi08j+ap6Orkm7ov7hE2bbADMNSAAAYoElEQVRb5eetpNkwSZdJOkXSuyRtV5lf3cG8pF0kPafV\nfTRTEn028Lc6y58ENqqzvFlFf2lrXTqPiBnADIANtLEvrZvZ+BPA6Lr9eyAippUsa+i/uMltR5Nm\nM/bOUwBIehy4nlQCvRa4kVTLujWpyrlpzQTRu0kjuJTZGbi3lUyUWEp6YxVTgGVtTN/MbHB0rggx\nmv/ipaSGPNXbzsnzp7SYZjM2BF5Fik87A9OAPUm9FVWOmKhfQKyrmYrAXwLTJb2hdoGkfYDppGFk\n2mUuMFXStpImkbpjmtXG9M3MBsYoq3PrGc1/8WxgP0kb5QZF+wGzI+Ie4HFJu+VWuUcAF7X0xuuI\niMciYk5EfD0i3h0RLyLVmL4FODuvdi9wQqv7aKYk+mXgHcBsSZcA80mR/JWkVlv3Aie3mpFaETEk\n6VjShzCRNIJMUx0Dm5mNGx3qbKHsv1jSScC8iJgl6dXAhaQA9VZJX4yIHSPiIUkns3rwkpMqjYyA\nY4AzSGNVX8IYNSrK3fpdAlwi6RfAD0ijurSkmc4W7pP0WlKz5DcBb64syhk6Np9dtE1EXEx7S7dm\nZgOnwRJly4r+iyPihKrnc1mzerZ6vZnAzIL580h91nZNRJwn6cPAqcC+raTRTEmUiPgz8OZcLN+O\nVJd8W0R4gO56Su4h1lB5Q4Cyvl2Hnum+vWqtWr/8mExYWXKM3QTNBo2bl7fqj8D7Wt24qSBakYOm\nxxY1M+sRHsVlbZKOJLXCXRARZUdoE9LdJS0pDaKSNmq1hDmabc3MrAWuXSkyk3RknpJ0AymgVqa7\ngT2AtwPfaXUH9UqiSyR9HTgtIh5sJDFJmwH/BHyE1LTYzMw6rcPXRPvYx0iNXyu3uezK2qcb9wGL\nJU0DboyIpnrtqRdEjwNOBD6bW+NWer+/PSL+CiDpWaQ+D3cjNTTaD3iI1P2fmZmNFVfnriUivl15\nLmky8HJSQK1MLwWeQ+obGGCVpD8B10fEdBpQGkQj4nRJZ5OGhPkAqa/ESq8PQ6RGRRMr+QMWA58D\nvhcRjzf4Hs3MrA1cEq0vlzDnUtWeR9I6pEBaHVhfThpvdHRBNO/0MeArkk4hdUL8OlIv/puRAupf\ngJuAORFxbXNvyczM2sZBFABJ3wR+BvwuIuoeldyH7vw8zczbTwBe3Oj+GmqdmzNyTZ6sWWWf44Ty\nz3d43eJlE58q7mRqnS3L+08euqedvTF2z9A+xb1Olt7GQnmr/yg5vmZ9yddEqx1LapfzoKRZpIB6\nWUSsaGTj3Iq34dHCxvv4H2ZmgyFGMQ2W55J6Q7qWNF7o/wAPSDonj+byrHburKX7RM3MrHcI3yda\nERH3k0b2miFpA+AfgLeRGr++E1gu6dekbgpnRcRfRrM/l0TNzAaBS6JryR3Q/zQiDiG15TkIOId0\nq8sPgGWSrpD0UUnPb2UfDqJmZv1uFCO4jJdrqRGxPCJmRcT7gC2AfYDvAdsC3yDdK9p0A1kHUTOz\nQTA8imkckHS6pPUhNR6KiN9ExEci4nmku09OJY0o05S610QlnQd80F34Naiko/my1rnb/KK8sdgd\nB0wqnL9qveJv/MJPblOa1gs/NRitc//85nUL59ftd7tk2Sbzff5og2W8lChH4ShSB0JPSLoAmF7p\nOCiPKDMP+GyziY70T/J2YIGktzabsJmZjaEOXhOVtL+khZIWSTquYPlkSefm5ddI2ibPP1zS/Kpp\nWNJOedmcnGZl2eajPQQjuIt0LRRSQ6P125HoSEF0L+BvwM8lzWx302AzM2uDSK1zW53qkTSR1C3e\nm0id7RwmaYea1Y4CHo6I7UjXF08FiIizI2KniNgJeC+wJCLmV213eGV5blXbSacAF0i6jnTq8H5J\ne+YWvC2rG0Qj4kpSF0jfJXWBdKOkfUazw0ZI+pqkP0m6QdKFktyZvZlZPZ0rie4CLIqIxbnDgnNI\n3cBWOxA4Mz8/H9hHWuv61mHAfzX7ttolIn5A6uLvHNKFniOBy4CHJS2W9DNJJ0g6oJl0R7xPNCKe\nBD6S65BnAr+S9H3g6pL1z2omAyUuBY6PiCFJpwLH407tzcxKdfCa6FakqtCKpayuFl1rnfy//Shp\nnM4HqtZ5F2sH3x9LWgVcAHxppG76RisiFgJflXQUsDvwOPAy4BXATsC+wCdoYhSyhjtbiIg5knYn\n9TH4wTxVE+mcZtRBNCJ+VfXyauDg0aZpZjawgtG2st1U0ryq1zMiYkZ+XtQ8rzbY1V1H0q7AExFx\nU9XywyPi7nyZ8AJSdW87CmEjiojtq15WGhW1pOEgKmlvUkl0U1JvEFe1utMmvR84t06+jgaOBliv\nPdeJW9fkSdTE31xXukxvfk3JPso2KN/P7V/frXD+i06/r3D+qkV3lCfWJhNfMrV02cIPblK8oOxP\not5FiZJtNp45Vl9fs84Tdf8CGvFAREwrWbYU2Lrq9RRgWck6S/PIKM8mDYtZcSg1VbkRcXd+fFzS\nT0nVxm0NopIuIwXI64FrI2JRwTq7AHdGREu3MYwYRCU9A/gqqS/CZcAbI+LSVnZWk+5lpHHcan0u\nIi7K63wOGALOLksnny3NANhAG7uRt5mNT53795sLTJW0LXA3KSC+u2adWaR2M1eRag4vr1TN5lFR\nDgH2rKycA+2GEfGApHVJXfNd1oG8752nSl4eJwXU60h9694IfIZ0AlBbRd2Qke4TfS1wBrAd8BPg\nnyLi0VZ2VCsi3jDCvqeTDuw+na4nNzPrd53qOzdf4zwWmE0aQ3pmRCyQdBIwLyJmAT8CfiJpEakE\nemhVEnsCSyNicdW8ycDsHEAnkgLoDzqQ/Q1JY4TunKdpOT+vY/Vph0h3obRkpJLob0kXht9WKR2O\nBUn7kxoSvS4inhir/ZqZ9a0OFjUi4mLg4pp5J1Q9f4pU2izadg6wW828v5GCWkflMbHn5AmAfEvL\n35NaCx8O3At8rdV9jHSf6IXAS8cygGb/ATwLuDTfhPu9Md6/mVn/cN+5Dcud0l8SEUeQAukzgSta\nTa9uSTT3fD/m8g27ZmbWIA+F1ryIOE/Sh0mdQ+zbShruQNTMbBB4KLRW/ZGa6uZmeFDudirrgL5M\nnfZSL/x08S0Yi75Z/FlHnbPQWKd4P7d+sKhxNDBhi8LZE58of38TVhUvW/l3LZweN/mj1sryfJUd\nR7OBEi6JFpF0JKkV7oKI0n/JTYAnW92Hg6iZ2SBwibLITNKReUrSDaSAWpnuBvYgDbTynVZ34CBq\nZtbnxPhrINSgjwGvZPVtLruy9unGfaQBuacBN0bE8mZ24CBqZjYANOwoWisivl15LmkyaUCVV1VN\nLyV1+nNaXm2VpD8B10fE9Eb24SBqZtbv3EBoRLmEOTdPwNM9J72UNQPry4EdST0wjchB1MxsALg6\nt3kRMUQaVGU+6fpppZvCFzeahoNoO41B74TbfaxwBDoWn1LSYT3l/bYPTy5eouGSlrYbNt/8T0Ml\nLWfrNWQu2c2EkrRe8Bm3wDVz69z2yK14b250fQdRM7NB4JJoVziImpn1u3HYfV+vcBA1M+tzwtW5\n3eIgamY2CDxiZFc4iJqZDQBX53aHg+iAeMFx5S1Un3hb8YDty/Ysbu06YXnx/JhYp0nthJJfcAs/\n7O2/e2/h/KHFS5pPzGw8CNCqziWfx3j+FmkA7R9GxCk1yycDZ5F6BXoQeFdELJG0DXALsDCvenVE\nfChvszNwBvAM0lilH43ov+K0R3ExMxsEHRrFRdJEUo8+bwJ2AA6TtEPNakcBD+dhLL9BGlqs4vaI\n2ClPH6qafzpwNDA1T/s394Z7g4OomdkA6OCg3LsAiyJicUSsAM4BDqxZ50DgzPz8fGAfqXxYK0lb\nAhtExFW59HkWcFALb7vrejqISvqUpJC0abfzYmbWsyL1ndvqBGwqaV7VdHRV6lsBd1W9XprnUbRO\n7gXoUdIQYwDbSrpe0hWS9qhaf+kIafaFnr0mKmlr0kjjd3Y7L2ZmPW90VxMfiIhpJcuKSpS1eytb\n5x7geRHxYL4G+nNJOzaYZl/o5ZLoN4DP0KcH1sxsrChaL4U2MPrLUmDrqtdTgGVl6+RO3Z8NPBQR\nyyPiQYCIuBa4HXhRXn/KCGn2hZ4MopIOAO6OiD92Oy9mZv2gg9dE5wJTJW0raRJwKDCrZp1ZrB71\n5GDg8ogISZvlhklIegGpAdHiiLgHeFzSbvna6RHARW05EGOsa9W5ki4jjeNW63PAZ4H9GkznaFIL\nL9Zj/bblr+9MmFi6aP0LrymcP/V/ij/+GBpqeh8Ml7Svr7dNiaGytMqUt18AlZwnRkn3Lv3Xwn7w\n1ft8/Xmt1qFDERFDko4FZpNucZkZEQsknQTMi4hZwI+An0haBDxECrQAewInSRoCVgEfioiH8rJj\nWH2LyyV56jtdC6IR8Yai+ZJeBmwL/DE37poCXCdpl4hY6wbCiJgBzADYQBv7F2Vm40+AVnXu7y8i\nLibdy1k974Sq508BhxRsdwFwQUma80hjefa1nmtYFBE3AptXXktaAkyLiAe6likzs17nIkRX9FwQ\nNTOz5rnbv+7o+SAaEdt0Ow9mZr2ugVa21gE9H0TNzGwEDXTfZ53hIDoomm3RSp1WuGWtIVvYR0vb\nNKteC80o2X+9Fp/WW9wCd0Qi3StqY89B1MxsAHSyda6VcxA1M+t3rs7tGgdRM7O+F6727hIHUTOz\nAeDWud3hIGpm1u8CVNKTpXWWg6itbTxUC42H92jji7/TXeEgamY2AFyd2x0OomZmg8Al0a5wEDUz\n63OK8H2iXeIgamY2CFwS7YqSEYvNzKyvRLQ+jUDS/pIWSlok6biC5ZMlnZuXXyNpmzx/X0nXSrox\nP+5dtc2cnOb8PG1em24/cEnUzKzfdXBQbkkTgdOAfYGlwFxJsyLi5qrVjgIejojtJB0KnAq8C3gA\neGtELJP0UmA2sFXVdofnwbn7lkuiZmaDoHMl0V2ARRGxOCJWAOcAB9ascyBwZn5+PrCPJEXE9RGx\nLM9fAKwnaXKb3nFPcBA1M+t7owigKYhuKmle1XR0VeJbAXdVvV7KmqXJNdaJiCHgUWCTmnXeAVwf\nEcur5v04V+V+XurPoZVcnWtm1u8CGF117gMRMa1kWVFwq91Z3XUk7Uiq4t2vavnhEXG3pGcBFwDv\nBc5qPMu9oWdLopI+ki86L5D01W7nx8yslymi5WkES4Gtq15PAZaVrSNpHeDZwEP59RTgQuCIiLi9\nskFE3J0fHwd+Sqo27js9WRKV9HpSHfvLI2J5v7baMjMbM527xWUuMFXStsDdwKHAu2vWmQVMB64C\nDgYuj4iQtCHwS+D4iPh9ZeUcaDeMiAckrQv8A3BZp95AJ/VkEAWOAU6p1J1HxP1dzo+ZWe+KgFWd\n6YE+IoYkHUtqWTsRmBkRCySdBMyLiFnAj4CfSFpEKoEemjc/FtgO+Lykz+d5+wF/A2bnADqRFEB/\n0JE30GG9GkRfBOwh6cvAU8CnImJu0Yr5AvjRAOux/tjl0Mysl3Sws4WIuBi4uGbeCVXPnwIOKdju\nS8CXSpLduZ157JauBVFJlwHPKVj0OVK+NgJ2A14NnCfpBRFrf0siYgYwA2ADbewuO8xsfHKPRV3R\ntSAaEW8oWybpGOBnOWj+QdIwsCnwl7HKn5lZ34iAVau6nYtxqVdb5/4c2BtA0ouASaSeL8zMrEgH\nu/2zcr16TXQmMFPSTcAKYHpRVa6ZmZHuyPR4ol3Rk0E0dy31nm7nw8ysbwx3pnWu1deTQdTMzJrh\natlucRA1M+t3gRsWdYmDqJnZIHBJtCscRM3M+l64YVGXOIiamfW7gHB1blc4iJqZDQJX53aFg6iZ\nWb+L8C0uXeIgamY2AFyd2x0OomZmfc/3iXZLr/ada2Zmjap0+9fqNAJJ+0taKGmRpOMKlk+WdG5e\nfo2kbaqWHZ/nL5T0xkbT7BcuiZqZ9bmgc9W5kiYCpwH7AkuBuZJmRcTNVasdBTwcEdtJOhQ4FXiX\npB1IA3TvCDwXuCwPKkIDafYFl0TNzPpdBMRw61N9uwCLImJx7tf8HODAmnUOBM7Mz88H9pGkPP+c\niFgeEXcAi3J6jaTZFxxEzcwGQAxHy9MItgLuqnq9NM8rXCcihoBHgU3qbNtImn1hoKpzH+fhBy6L\n8//cpd1vyvge83Q8v3+/9/Gpm+/9+dUvHufh2ZcNn7fpKNJbT9K8qtczImJGfq6C9Wsjb9k6ZfOL\nCnB92TJqoIJoRGzWrX1LmhcR07q1/24bz+/f793vvdsiYv8OJr8U2Lrq9RRgWck6SyWtAzwbeGiE\nbUdKsy+4OtfMzOqZC0yVtK2kSaSGQrNq1pkFTM/PDwYuj4jI8w/NrXe3BaYCf2gwzb4wUCVRMzNr\nr4gYknQsMBuYCMyMiAWSTgLmRcQs4EfATyQtIpVAD83bLpB0HnAzMAR8OCJWARSlOdbvrR0UvkG3\nLSQdXXUNYdwZz+/f793v3cYvB1EzM7MW+ZqomZlZixxER0nSIZIWSBqWNK1mWWF3V4NI0omS7pY0\nP09v7naeOm1Qui1rhaQlkm7Mn/W8kbfob5JmSrpf0k1V8zaWdKmk2/LjRt3Mo3WHg+jo3QS8Hfht\n9cya7q72B76bu88aZN+IiJ3ydHG3M9NJVV2hvQnYATgsf+bjyevzZ90Tt3l02Bmk33G144BfR8RU\n4Nf5tY0zDqKjFBG3RMTCgkVl3V3ZYBiYbstsZBHxW1Kr02rVXd2dCRw0ppmynuAg2jkD061VE46V\ndEOu+hr0qq3x+PlWC+BXkq6VdHS3M9MlW0TEPQD5cfMu58e6wPeJNkDSZcBzChZ9LiIuKtusYF5f\nN4WudxyA04GTSe/xZODrwPvHLndjbuA+3yb9fUQsk7Q5cKmkP+XSmtm44iDagIh4QwubNdJVVl9p\n9DhI+gHwiw5np9sG7vNtRkQsy4/3S7qQVL093oLofZK2jIh7JG0J3N/tDNnYc3Vu55R1dzWQ8p9I\nxdtIDa4G2cB0W9YsSc+U9KzKc2A/Bv/zLlLd1d10oKxWygaYS6KjJOltwHeAzYBfSpofEW+s193V\ngPqqpJ1IVZpLgA92NzudVdYVWpezNVa2AC5Mw0WyDvDTiPjf7mapsyT9F7AXsKmkpcAXgFOA8yQd\nBdwJHNK9HFq3uMciMzOzFrk618zMrEUOomZmZi1yEDUzM2uRg6iZmVmLHETNzMxa5CBqVoekzSU9\nJukDXdq/JF0n6cfd2L+Z1ecgagNF0r9KCknvK1gmSVdIWi7ppQ0m+SVSTzRdCWKR7kE7ETgi34dr\nZj3E94naQMm9B11L6pLvZRFxV9WyjwP/DhwfEac0kNYUUscRn4yIb3Umx42RdDtwXUT4hn6zHuKS\nqA2UPCzZdOCZwA8r8yVtD3wZuAb4WoPJfZDUA9N/tTmbrfhP4EBJRQMAmFmXOIjawImI64CvAPtJ\nOjoPoH0WaeSV6U10v3gIMC8i1upYPFcNf0DSNZL+mqcbJZ1Utc6RuWp5H0knSPqzpCfzNrvldV4n\n6UpJf5N0j6TPl+TlEmBdPGalWU9xELVBdTIwH/g3Ut/GuwCfLRlAfS2StgC2p3zQgJ8AM0gl1S8D\nnwYuBw4uWPcUUvD7FvBF4AXAbEkHAT8Dfgd8CvgTcJKk9xSkcR2wnNR/q5n1CHdAbwMpIlZKmg7M\nA44BriQFsUbtkB9vr10g6Z3A4aQq1ukRMVy1rOjEdCKwW65qRtLNpBE/zgdeExFz8/wfAX8GPpzT\nrn4/K3LH5zs28R7MrMNcErVB9hip9AZwcXWwa8Bm+fGhgmWH58dP1aZZso/TKwE0+11+vLoSQPO2\nK0gl36kleXoQ2HykjJvZ2HEQtYGkNE7Xj4FJwC3Av0h6YRNJVJqtq2DZVOCeiLivwbQWr5FwxMP5\n6R0F6z4MbFKSjqryZWY9wEHUBtVHSNcPv0hqILQOMDMH10b8JT9uXLCs2WBW1pCp2fFlN2Z1vsys\nBziI2sCRNJXUOncucGoeLPuLwJ6k4NqIBaRAWVS1uhB4bm58NCYkTSbd+3rTWO3TzEbmIGoDJTfs\nOYPUmKf6dpZTSY2MvtJItW5E/AW4GditYPHZ+fGrtQ2JmijpNuuVpKrpKzqUvpm1wEHUBs0ngdcC\nJ0TELZWZOZgeSXPVuv8N7CRpy+qZEfHfwLnAEcCVko6T9EFJXwdubM/bWMtbgJXAzzuUvpm1wEHU\nBoaklwAnAVcDX69dnqt1T6Txat0fkKp0312w7N3AscB6wAmk+1H3JQXeTjgcuCgi7u1Q+mbWAved\na1aHpO8B+wHbR8TKLuXhQFKnDNMi4vpu5MHMijmImtUhaXPgNuDTETGjS3m4DvhjRKw1Mo2ZdZeD\nqJmZWYt8TdTMzKxFDqJmZmYtchA1MzNrkYOomZlZixxEzczMWuQgamZm1iIHUTMzsxY5iJqZmbXo\n/wOn5hpoZmWx2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cbddba4400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = np.concatenate((d1, d2)) * 0.2\n",
    "G = np.concatenate((G1, G2))\n",
    "\n",
    "alphas = np.logspace(-6, 0, 50)\n",
    "norms, residuals = algorithms.trace_lcurve(d, G, alphas)\n",
    "\n",
    "plt.figure()\n",
    "c_alphas, curve = algorithms.lcurve_curvature(alphas, norms, residuals)\n",
    "plt.semilogx(c_alphas, curve)\n",
    "\n",
    "best_alpha_arg = np.argmax(curve)\n",
    "best_alpha = c_alphas[best_alpha_arg]\n",
    "\n",
    "norms[0] *= 10\n",
    "plt.figure()\n",
    "plt.loglog(residuals, norms, 'k', lw=2)\n",
    "plt.xlabel(r'$\\log{||A X - Y||_2^2}$', size=20)\n",
    "plt.ylabel(r'$\\log{||X||_2^2}$', size=20)\n",
    "plt.title(r'L-Curve Criterion with Optimal $\\alpha$ = {:.3f}'.format(best_alpha), size=15)\n",
    "plt.scatter(residuals[best_alpha_arg], norms[best_alpha_arg], c='red', s=100, edgecolors='none')\n",
    "plt.tight_layout()\n",
    "\n",
    "m_alpha = algorithms.solve_tikhonov(d, G, 0.01)\n",
    "m_alpha = m_alpha.reshape((30, 50))\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "ax = plt.gca()\n",
    "im = plt.imshow(m_alpha, interpolation='none', vmax=0.2, vmin=0, extent=[-25. / 2, 25. / 2, -15. / 2, 15. / 2],\n",
    "                cmap='viridis')\n",
    "plt.title('Single & Double Reconstruction', size=20)\n",
    "plt.xlabel('X (cm)', size=18)\n",
    "plt.ylabel('Y (cm)', size=18)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "cb = plt.colorbar(im, cax=cax)\n",
    "cb.set_label(r'$\\mu_f / \\mu_{total}$', size=20, labelpad=15)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}