TomographyReconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "51dec861-35e3-4f08-91a9-191ba91e9219",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "from cil.framework import  AcquisitionGeometry\n",
    "from cil.optimisation.functions import L2NormSquared, BlockFunction, MixedL21Norm, IndicatorBox\n",
    "from cil.optimisation.operators import GradientOperator, BlockOperator\n",
    "from cil.optimisation.algorithms import PDHG\n",
    "from cil.plugins.astra.operators import ProjectionOperator\n",
    "from cil.utilities.display import show2D\n",
    "from cil.utilities import dataexample\n",
    "from cil.plugins.astra.utilities import convert_geometry_to_astra\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from cvxpy import *\n",
    "import astra\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4af233b2-a324-41c1-a82a-8b77854d7287",
   "metadata": {},
   "source": [
    "## Create AcquisitionGeometry, ImageGeometry and ProjectionOperator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "36a3c313-e8ec-49d6-9223-039b15ded289",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Detectors\n",
    "N = 32\n",
    "detectors =  N\n",
    "\n",
    "# Angles\n",
    "angles = np.linspace(0,180,90, dtype='float32')\n",
    "\n",
    "# Setup acquisition geometry\n",
    "ag = AcquisitionGeometry.create_Parallel2D()\\\n",
    "                        .set_angles(angles)\\\n",
    "                        .set_panel(detectors)\n",
    "# Get image geometry\n",
    "ig = ag.get_ImageGeometry()\n",
    "\n",
    "# Create projection operator using Astra-Toolbox. Available CPU/CPU\n",
    "A = ProjectionOperator(ig, ag, device = 'cpu')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fae01aa-b7e5-4bc3-b93f-85a0880d5af1",
   "metadata": {},
   "source": [
    "## Create phantom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7722a2c9-2f98-4a34-bb8b-2d7b2f16f477",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get phantom\n",
    "phantom = dataexample.SIMPLE_PHANTOM_2D.get(size=(N, N))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f182be1a-443d-420f-91e3-53607b60d07a",
   "metadata": {},
   "source": [
    "## Create sinogram "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0c4a9704-1676-4159-9233-aeea17014199",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an acqusition data\n",
    "sino = A.direct(phantom)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e00b3f2-40c0-41ab-8882-e9fd13af8751",
   "metadata": {},
   "source": [
    "## Get sparse projection matrix from the ProjectionOperator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "446aaabf-fbe5-4453-a398-1be4c81a89af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to astra geometries\n",
    "ig_astra, ag_astra = convert_geometry_to_astra(ig, ag)\n",
    "\n",
    "# projection id\n",
    "proj_id = astra.create_projector('line', ag_astra, ig_astra)\n",
    "\n",
    "# matrix id\n",
    "matrix_id = astra.projector.matrix(proj_id)\n",
    "\n",
    "# Get the projection matrix as a Scipy sparse matrix.\n",
    "W = astra.matrix.get(matrix_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3f9cd8e-f9c5-4b8d-b0ed-667e6cd9ca29",
   "metadata": {},
   "source": [
    "## Compare matrix-free with matrix*vector implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "08a1c5fc-c51d-4937-8bf2-9ac9b78b5f97",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJgAAAEDCAYAAADJFM0VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABKhklEQVR4nO29e4xtW3bW9xtzrrX3rtd53Uff9u222223DZIjg4NwoBHhjUMsbKEQ2QmIV0T+wYEElNhIiSUkpI4UIVpKImHZBkwaCLRtYSGwaRssIIKWuxuDsdvtR7t9+z7PufecOqeq9mutNUf+GHOux37U4579Ot37k0pVtffaa62qPfaYY47xjW+IqrLHHuuC2/YN7PGljb2B7bFW7A1sj7Vib2B7rBV7A9tjrdgb2B5rxd7AnhGIyA+KyH0R+Y8rOl8lIj8bv35sFedceJ19HuzZgIj8TuAc+CFV/YYVnO9cVY+f/s4ux96DPSNQ1X8JPGw/JiJfIyI/LiKfFpF/JSK/YUu3txR7A3u28X3Ad6nqfwr8JeD/vsFrByLyKRH5tyLy7Wu5OyBb14n3WC9E5Bj47cA/FJH0cD8+90eAv7LgZa+p6h+MP3+lqr4uIh8E/rmI/Jyq/uqq73NvYM8uHHCqqr9p9glV/RHgRy57saq+Hr9/XkR+GvjNwMoNbL9EPqNQ1SfAr4nIHwUQwzde57UicldEkrd7Hvgw8AvruM+9gT0jEJG/B/wb4OtF5FUR+TPAfwv8GRH598DPA992zdP9RuBT8XX/AviIqq7FwJ4qTSEi3wJ8FPDA96vqR1Z1Y3t8aeBdG5iIeOCXgN8PvAr8DPCd6/ok7PFs4mmC/N8K/Iqqfh5ARP4+5qKXGlhP+jrg6Cku+eWHMRdMdSJXH7kcf/B3H+k7D6trHfvp/zD5CVX9lqe5XhtPY2AvA19s/f4q8M2XvWDAEd8sv/cpLvnlh0/qTz31Od5+WPHJn3jftY7N3/urzz/1BVt4GgNb9KmaW29F5M8CfxZgwOFTXG7uxCBX71HcoI9kq8nGaFkSxpNrHBhgp0pwSqVhK1d+mv/8q8D7W7+/D3h99iBV/T4s48wtubfS/7q4K1YO75GDATIYrOaC4zFSVVBdvtzo9VajjUGBMP/Z3wiexsB+BviQiHw18BrwHcB/s5K7ug7ELfdgThARcPa8liXiHPTy5eerFnzCvWueqyoQMW8oglbBPNWiexDdOSsLPGMeTFVLEflzwE9gaYofVNWfX9mdXQfJkKA2JsAModezn4sSnRbI0SHa79nSOgtVZFp0lzXn0Nz+PTKZmoFFj4gqXAzRMtg9eA8hUO/Id8u2UJTiGVwiUdV/AvyTFd3Lu70HZJHRtBHiP1cEnECIhpCW2BCfU20MsH3OoI3xiZv3TiF0j90xKFA9g0vkdqEBLUr7ceYpcYJOi8a7eN88GdQ8UllBL0f7uRlPUZqX6vfQXg4hIKOJeabJBKZFc2lVtApojMd0dpnckre4DM9iDLZd6PI4p35/RZB+v2tg0Zh0WiDeAdHAqgotK6SPxV4hQFGgZYlOplAUtgxHT6dV1Xi1HYu3ZqFAtaVd7bNrYDAfT9XLmKUw2rtMLUskpRjSLrCskPHUjCm9djJFqmCeq2oeV1UkBDOyEMxL6pKleadSFIZV+lQRuQN8P/ANmP3+aVX9N4uOfXYNrJUHS4ZUexVxSJ5ZbBYDcKYFqhetGMvZMlraMptiJx2N0Co0r4U6xtLo6dLrpeUYAbSOv3YrD6boqmOwjwI/rqr/lYj0YHmC89kzsEUBfYqBxDI+HbQ9UOv3+riwIM0QAupck0lWNQOczbuJWxxviUO8mMGF7S+fqlCsyL5E5BbwO4E/aefWKTBddvyzZWAiuH7fUghFiZZFfDOdLV1OGg9TdT+z0u83xlmW5qWyDLIZA8sy2xjE6816ovScUpnRzRqZBtzhIXJ4CGVJ9fjJDhiZUC0svCzE8yLyqdbv3xeT5QkfBB4AfzPyzz4N/HlVvVh0smfLwMByUd63lkPMuLxHQkC9N+MK2rzx4pC2oQSNcddi71O/FzOGA742UhFBHd20RDxe8gzp91AnFqtteVOp3Ch78raq/pZLns+Ab8J6AT4pIh8Fvhv4X5cdvPsQQbK8XqK0FcyLk26MlP6T6c3OcvNwIpaRdwK9HCFfmIVPx2pZ2e4xesQ6cRuDfJwzg3Y0idV4Po2bCS1K+yDsAG7gwa7Cq8CrqvrJ+PvHMQNbiGfCwMR7pJd3s/XQvKEp+F7wMZU8a4rdVQV45HAAmYeyskC/8wKBLLPnRiO0qnD9fm1gqtoxstpbthAmExiPV/K3rwKWaF2NganqmyLyRRH5elX9HPB7uYSi9WwYWK+HOzmOmXh7c7Uc2A6wl8PJkXmUswv0Ymi5q5SSaJdwAGmvV1Vl53CC5HmT45pO7XFxceMQd6guxnvE5G4IFhMeHkJVEUbjnfFYbShQ6ErZ8d8FfCzuID8P/KllBz4TBuZu36L8qhfR1i5OgiJloLjd5/FX96h6wt1fmnDw+XcsU//oFK0q2z0W0Uv5JoayZbAknJ0hvR5ycmy1x9MzwunjaHQZkmdWMShLpN9HjuKOfDxGQ0Bu36J68TZuXCKvvI6en+9UigJAEaoVtl+o6s8Cl8VpNZ4JA8N7QubAC5qC7EoR76gGjuJYCD2oBh7NM6PUpKBftalXpp/TeatgCdjIkEAEQoWWRVNikibHJu1EbkLmCX37Nzrvds64EsKypPCa8UwYWHj8hN4XXfeNjW9k9qBH/50j1DuyhxfIE/MgMhhAvweTqS150GwGzs7BiS2jkXqjZ+fgfSwheVseUxE7LqtaVeh4HJ9TO04VVwSkWBADOo+LsaNOp2hK6m4Yq4zBbopnw8DOz9HRqPtgewf4K/bv014OvR7kmeWhVGMs1dp9VpUteRrq4FyrCj2/MEOoKju3k07sBthr46ZAssx2pKpIUSFlmDtenCADq4Vqive2AqFabQx2bTwTBgbM7dTEzSeXNBmUKogZpBZlN33RnCCWemJSNRaybWmkm/qIxzPLzACYTHFPhtFwuztSDWrGnAiLW4LVN/YGthwLmBOq85x8LcqawiMzHq9toHWdMrEjWukPy4PF9ENaWsUhLjTpEqjTIuHsHH10CkGtstBGqAjDYfM3bAmqwlT91QeuAVea9SLhMxG5JyKfEJFfjt/vrvc2l9zbIk6+hnr5S1/LT2AeS7xDpMWOvSz1HloeLV2/qlpVgwUMjx0I/ANyra9V4zp+828Bs31y3w38lKp+CPgpLsnkrgxpl+c8kmW4fh93eIg7GNSGJpFgaNl7X38ldoXr5bheXhuS9HrI0WHcEPQhpiTCxchyWjOGoUVpydfxxKoBx8e4k2Pcndu44yOrNohr7nVHYEG+u9bXqnHlGRcJn2ENtn87/vy3gW9f7W0tQGzyEBeXxjy3N7nX6y6V3kcmayqCS8PdTzGUi8VxHxtB8lhOcrYr1GI6v9wBaCBMi24TyaBvhfR+v07CNuyOXTEyC/Kv87VqvNsY7D2q+gaAqr4hIi8uO3AVfZGS93AHrdYzJ8jtW4STI6Qo8f1ek89qJ1PB6okhphBi4rUOuKcFKrGkk3oZvcelNrcU0CdP5mPaoc0zK0qjVFeRoEhcuiW3pbrNfN0SvqSD/FX0Rcqgj9y+Zb+EACJUL9xm/OIBfhzo9TOkqAi5R/MW1aZS3MXYOoaKEh2PkaBoZLDqdIpOJhZ/9XpmuFnWtLfNMCVkMECPD+3c08LqlZMJOhxFBmyMxXzPzldV8+mVLaF6xhKtb4nIe6P3ei9wf5U3tRDRsOovVcvmBzXuYLDfoeocY3z7MJc07XQCkSjRgAtGQnTS9EWqQhWPSf2TZWzATfRqVfNu0BTDwyWbhQ1CEQrdTsLg3V71x4A/AXwkfv9HK7ujNpyPHUJTwuljSy8cHYF3uOGU/tsgRYW7GNkbXhRQBYuF8vhmx8aN1NhRL5VVsMA/y2KeLPbdlPaz9HrIrRPL9I8n9trJxJZDsOUwhNgsYsuze+6edSo9OSM8fhI5advfQaYgfxu40sCi8NnvwpiOrwLfixnWP4giaK8Af3QdN1cH9EEtg66KHBFpMgE3LqNhWcuZprJQnjVvbFnOJzmdq+Mt8syWzem06RYKsV7ZZsqm3FjN4Q91AleD1kG/DnowzHaKVaHI7i6RqvqdS57amEyOeFd3VYe7t9CDRgJAJgUyjEvmoI8cHphRZDEWS7u+tGz2e1Qv3aE8zGvmqh9X5K+8jQ6HiPNNLOa9sVIHffOINcHRDFKi8WpRxt2sPSe9Hv7eXXtuPKmDf8AC/y2UjL5kg/yngQa1co5ztmQNBhR3BpRHGa6wIrOHJuE56Fs90gnqPRJZFJK8kPeEkwNOP3TE+J4gASRA74ly7/ExLgTzQnnW4fNrP7JZ84zQy+y85yNkEruSsqIhNaqigx6SZ+ZZofmOtc9tuiapyr4WuRRpq1+WMJmSPRnjyh5SBvuaFMaYKIwFIbH9v87Kp9Y0EXAV4h0HD0tc5WsDy88r65msKiikyWenc4RgwXxR4krznhI5YlpaoK9aIBOPFK7eXGgZS1ft5XIL0gIW5G+nVLTbBhYq2/SlZWg0Qi4ucHVJRy3Hk3jzF8OmrpgozdAwVasK8Z7DJ+ccZZnFUSHYmz6ZEFTr5Rjn68RpuBjakuh9zaJIcgFalk3N8jw21rT6AuZisS0F/Tsb5O8ENCpcVaCpeJzQFqJbVJv0viuOEoLlrXxqfSubLL9bLAkl0RglbgpSotVe51oKP8njtXoxq6rZECT+/obb2BTZEw4XIhqPOxgYvytUhCfnTRkndXbHjHvdH5loN1DvAOX2CdVzJ4TcUxznhNzRf2dMdv8xgNGlnUP7mSVrU1zuheH7Dhnf9qiDFMpkE0UqQEAFQiZM7hizFjUqv5vC4FHAT5TBOyX5ozH+bEj49dfQYmmv6lqw92CLEGuPMugjJ0dQVUirsaJmU7QF5yI7grzX8Wjh1iGjrziiHAiT246qJxwPHMdTO5fmHvVCNcgIPUuNuEqpeo7Tr8m4eDmYMTmL27Kh4KaCeiVkEPoK7x1xcDBFVQgKF8MeF68NyEaOo9ccR295Bg887lWHLih1rgvWF7k3sHloQIOzRGcWSX2pHCMu5p9i2sBh4iRgzRgzrf7ubMTgQY+Qe3pnGSET+g8nyLmVciR6PpdndblJQiDLHLe+4MnPXO2tRMGPwVUBjQ24IRcmbx4S8lhvVTiaQv+hkk0CB2+X9B5NcKcXhEVqimvFjTq7V4odNzAjGobhEEatojTUcYyqIFlMxhLLQtqoYSU+mBtP8I9O8eLI63NoS5UwiZrEXaSLHDEn3PqC55abUzpp7vEyJA2z0ppwQztZuyFY29p+F7kcs4zWOVLfvFaqtnZyEEs7kU+fsvVpaV1IqYamFtnamdZJ03ZA3+L3t++lVv1JMgZb2kHakr1fIq+PttQlyZiqlmhJW1qpFexDtwAdBVKWsV61qhAnuNu3jEkxHqNPzmMLW9d4xHvrm0w7zbbH04CeXxC22O29T7Q+DZLX0DDflBEhIixiN1xVM9QQiYi9vJOB73is6Jkkz0ydJ5Iha21XVXS0PeOyJM8+Brs50rJY8/AdQtN2BjQ5KXRuWYMmxdFpCnHS9WqTiXWVTyZL+fqaCvJVMEOcTrsJ11kNjI1itW1rIvIF4AyTfSkvU+N5dg0spjA6hqCWjO20o8XHgVpsrpPmqFUMy85jIo0RhtHYjCtea+668RpzU0A68k9b7CpiLZ3dv1tV377qoGfXwCLESRP/twP9dpqi7em6L+6eZ8ZotDaoSEJsB/Od87S6kdqCdDvABYN9LfKpIe24K77BJjEeeyBjZ7UWjQdaeJ7ZJTTGV4pHXOiWeRIZsi3BmV4DO2NcCTeg61ylcAjmFP+ZiCjwNxY8X+M6hMP3Az8EvIQVUL5PVT8qIveA/xf4APAF4L9W1UfX/SuuhSVTOZrnZzxTir2CqcI1bWxtfbDW6xYZWvscM5qvbWXpRIasd6sxBuzcT/vv2Grj7Y04+VcpHAJ8OA6UfxH4hIj8Yuw+m8N1PFgJ/EVV/YyInACfFpFPYCKwP6WqHxGR78Z6I/+X6/4VV8J5/AvP2egWqLn1Sf+rlh4PFnfV3iPE2uW9O2iW2tQEKStcbOPX0agmIup0aoYSUwxyeACDfre3sSgbdkaeWcDfYnRIkuScTI2dEc+9EEVhAnUbNrhVxmCtgfL3ReRHsdmh787AYntaalE7E5HPYrMivw2jUoP1Rv40KzQwyTO4d5vy9oE9EBQpAt47GE+MizUtqKUx6xsOcDCgfPE2oRdb2BxIqbhpZRz+R4L6CQxHkcJjTbj0cvTWMeHWAeqkLhm5cYEbl4ReRjjMOzplycCkDPiLKRSl6ZPFgH92zI1OprZh2CCMTbGaXaSIkdajLRwBfwD4K8uOv1EMJiIfAH4z8Emu2Rv5rvsiqwoZTfBt2cyysmx8GqmXvFha0lLROyhuWOCmrRxXGYw6U1Ywmc5z9dUaOKQokXFpq+gk6VxMkaLEFZnVO0XQzBmdugrWzRR7JEn8su4/wYiJYfNlImDVCofvAX40fnAy4O+q6o8vO/jaBiYix8APA39BVZ9cOYAq4t32RWpVEe6/3cgkeW9iIzPLY7w38JE7n2WWWX/zAfXYmGhIdX2yXlrjrk+cdRlRoudD3HhSL3VpJpImUbvYHe6OjyxjX5TWzQR2bTDjTUYWZdJ1PDbZga10Gq3Og8UR2t943eOvZWAikmPG9TFV/ZH48Hp7IzVKH6XYx8dlMHmdlvZqezZk/dppEbVcy0uDe2mXkERs6dVQx2ptoqC40rj+3jL7EuWiatkmte4iVV1eH92SpvnOZvLFXNUPAJ9V1b/Wemq9vZESeWDtIaJB0aKoh4FqUhkc9KxoPS3Q0bgzBqaWH2/NI6qL3ekcbsYQoSsxQKuFLp1rNLYeyaBNuUnHKNg9R0+qcUWULMMfH5ku7Gi0US92w13kSnEdD/Zh4I8DPyciPxsf+8usuzdSHBLrejWlpjQxuZoTBoCvNb00Nsaq87hBYra6enmd06aI0GAyuXXhu+1l2oIr7ddMu4xUTR5XLeNUp0ZS3sw7q1Ou5r9zY+wsm0JV/zUs/b+srzdS45xGbYLmRKupE5ypBlmUsV2/YVnYY/G2O9IB2qlFdhCF5jQsWNpiXo3grK65ALbcxum3k64ubOKXzRrmJrDn5C+Cqi0lo1G3pJP0vqIhaOwIqks0dbdP5H4V5bwHahmXzJaT4rntsK4hWklqhn3RHhsYG0FsxuSSVMQWEq4KlLvqwbaKmZmP9eArjEQntIJ3J5fPym4H3cuGybeuY6gaw53hoHVPrQ25MBlQOn72WsLGu4pgh5fIrUEEf3zUSCnNlIUkCeu2vJOORi0V6Lx7vqSCmKSa0tCslLKIqYU6ZgOkcs0uMqpPywzfLMWDWmH3U3tSK1O5w8MofBdb4iYTqvOLzRqZ7pfIeYiVZSSfmVGUtL0kRONq90W2qMpLCtqp84iYX5OZsX82LsbHVIOJoHToPTPnFRW0bFOGqsbTJXZrlllTifeotqblbmi53BMOF0GDsUCL+cz3nB59u8UfzNMsWM6St+rMlExVgZTrymNit9J5Pf2gxhODhs8vkfG6aNMQxU9wsaQVh9VvQ/Vw78FmoWrdRAumewCxtT8NB03Ph1avpOu+Nua+bPxeHNyQtCPaUpfJSEOwjULHEAJE+fQ6uasa+fnN0pk2B1pVDVV6m00f7A2shr9zG977ItrLqI77VH1P3Y8YwI9LpIy7PMBfTNFXXp+XqkwzHY+OLPEZdShExKSYnMDF0EbIJGqO0CzHaals58zqXsxgm4wwky+LsgISBe7aEJ/Xyoch0ac3FIcpQrlodPQGsHMGVv7GD/DKtxwxvVfxnq99m99w9z49V3LsJzwsjvjZt17mYtinKh2h8Ay+cMIH/86Y8td+vSYBooEwKXD9Pvre55g+1zTDqhNCLqgTDl47Q351bB4t7k6TLpiAKUfrzFIZqriDxZbHmtvvcLdOIM9NKCV5rsQbOz5CD/rIZIo8fmL1zclkY8XvfQwWEQae6XMV2XNjfsd7Ps+HT36ZEzfixI15UJ2QSeC14W2GRY+LaY+Hj+8Z76uFVFBWVcg95UGb8UptYJp7C747S29r41APHm1IjHYBrT1Z8zrL1JN5pHSRMxbTJ86CfR3kJGVFCcvStWuA7pfIGv3PP+Ar/+lLFEcH/OQLv40fP/ptjeiIQjYGKUEqcKXyFW8HeHLeeK+UhI0xkX/njIOiqompkoSBVZEnFyb9FEKjrgPx+dBQgubKR2KCvynVUcV468lZ67FErXaIlNDvIWOPTIulfZjrwj4Ga6H89S/Sf+VV+sBxu0M6Teu4fctyY1WoWRN6cTHHy09UZj19jCTdrlgmCsNhk+GflXdKCNopijfnjp4pz5EsiwLD1qJWxet0upXSa4oCmdoYZ7uV/S5ye6j/+antH+rRxtNp16tEzYcOLz8hTjuT5DGcRHXo1ptbVWgrgapVaKajXeZpQqJoN+erO5Nm78GBTCZRxTrStBfM+l4XFKHaB/kLMJMi0ArC+UXtcWqV55ROcJ7ZMX+a3thYZqo5ZMkYQmWvq1zNQdOymGNRdPUlQi07UOfUUrA/k4hN3U1VcQ6cd47dZNpiH+S30WqoSL/XS2DKZS1o3b8KnWVJYtC+SEglGZHMlJvSvYkZqsyeL9U7YzdT2iDMCrFsOh+m+yC/gWQZcmCNHiYBXuEOBriTYwumx2MbrpA8fmglVJN8gPdNTirPOsE4UAfvks//+XWaokVaNI3+ohY4sQe1YdXWXis2iRweIoO+0a6HI8uXuehBi9issuE5RrqrBiYiA6wlqR+P/7iqfu/a+iK9R5JseGFlFen14PDA3qhiilauHvhuN9nupg6Ab5QO+/06GJf0xtY07HhMCuhVzRDsDzcDDIoSl2BxNUNWJ5OFgbqIWLpi0EfKqmkuSYYb53w3G4hNJFt3u9g9AX6Pqp5Hbv6/FpF/CvwR1tEX2SqvJI+j0ymcXZBERNJU21nmaQeJYx+ZrHVw3fI8QiwzpZEwLYOpdcPAlry8Zx4oNp3UNcyW0Eq9G41xn3hn2rJpAH1iehweIKqmN7upbP6uGpjafzlGp+TxS1lTX6SW5Vx2OwyHtgubTRks41zF57SqkFBZ0jNlzts7t2R80Bhtenn8Ls7oPdK38TAaVXPMA8UNAr6rBzsao+cX+FvH6PN30dzj7z+iOn1sTcHHJ3bu8WQjYsDWs7wdA7vW3lVEfOTj3wc+oapzfZHA0r5IEfmUiHyq4N03nHZ5+PXJl78gNXkk7xWz77PnqKffdpppL+n8aetWdOTRW3KcGqykVJTWZ1lUNd9Mq2C5sHKzjIptjVS+VpCvqhXwm0TkDtZ0+Q3XvcBTz4tcsKNsnlvkuVpSAjC/nKacV2sKbloqk2RTCovScYmQqFXVnYQ7q41RWX1SY2+ATibw+luICGEysdeOAiG2uYUNaYYpO7xEtqGqpyLy09gM783NjFz0SV9iXOl7I0Qy8ya2X5dGKxNjrlRuXCSOUid3U1wWWa+zxMYOh79Cz887978oBFg/thfkX7lEisgL0XMhIgfA7wN+kaYvEtY5MzLFWS6Nd/FztOX28Hd7IOrrt5s9oudKbItZsmFzstZrExFxhoeW6p5zSdU42i+dp/5y/vLlfANolWAv/boOYsj070TkH1917HU82HuBvy1ic8+Af6Cq/1hE/g0bmBmZIC12gqQO7NYuTkQ6y1vNq6cyFnOrPpjoN/Uuso269axpqBXROu8214XUgcWJHWOvVa5X8V9491jxEvnngc8Ct6468Dq7yP+ACZ7MPv4OK+yLlH6/mXNdv7HSHWgF8x6luZ/oMbSjCZZ6GoEraotNSmL+5KEWV+lIc6bXdXJwNFQeaTWQQD2Mq25vS3/HmrXzbRe5MnWd9wH/JfBXgf/pquN3I5PvPP6F5wl3TpDJFHcW2Q95boNFU/xTFITTx03zaqsuCJgoic/NUFPGvTpDyxlDW4QkIzCTcmjiuQpxGS6eNyVmlQptq/jEe9IKJHN2vPd2z2WJ6/eRk+Oabg2gF0Oq09O17ipXeOq/DvzPwMl1Dt4NAwMj5PU9ErKG6pJ5NE0zkwAhcXBS/isWvWNfYjM1TZrCtpvZ6S1CO3aaRbtKAI3xqesyLtrvYLuHMnUixc5u62jyHQOrBz6sETdYIpdKaIrItwL3VfXTIvK7rnOy3TGw8QR3liGjCeHJGfX4vNZsxrqeKFIH451YKNUgAR2aMk6dvW8nZJP2aqtLqI6t0nJIi3EBJKZFCKM5D2cHL15atYpzk9LSnkpW7bmWa95VKnITA7tMQvPDwB8WkT8EDIBbIvL/qOofW3ay3TAwDZHHnqGjEWFma9/pIFokQY7Fa5L6HcvSpMdTQTkdk7qv00Mau4QW3E/7NfXxoWpKO7UI8BVrTwhNOSoZWlE0giw0JbF1YhUrpKp+D/A9ANGD/aXLjAt2xcAgymGOGtGQGdRLXk2/0kY6AEtwUpZ1O1rdpNHOQS1Yzt4t6a9Dz14S2mmIHqsd0CdPmrhoqVy1zqy+gm6pVLQbBqZqXuuipQPRgs3/ya0drChj0qaaC6sSURCYbwnTBcyFRY9dF5HCbfezZInTCp2GmjFbD0xNiBN3Z+ug68CqM/mq+tNY/flS7IaBwdVv9oI5QwvPwfqXmxqJNn0Z0j2pNL2UwYiO2p5/uWZsqed3hwzsEhglOi1513lDuTRznhKhT/Pm1qIn11li04YC4+HXOa9299Ea8czUIjeG2QK36uJE5OxxC881owtWl28EwaPlJQa2RFMs/T7ncJfJHKTjU29luWEZAWXxZmYD2E0Dm/3nt3aRHa8xmyKY49yz1DDSvMe5ks8ij3YZsRGoWS7tzqYFjzW1y/i3tAXr1p2q2C+RyyHem0Q5mGR40uwKbrFBiHRTEnPeJG4Q2s0kzVN13mzuHDJf4LbrucZDJQEWn+qZ6ZDGoFOwL6l5d7xuCQH5Mt9FXoXUf7jo8TZuwlhI6YI5BcIFmw1xWJfRNRUSF6BWQUzGVydab0BjeBrsPdhyaFU1OqyLMNMBPhd8zxheh+2wSFBOW8bnJGb7/cLjgW4LXX3JFr8sFcFLbbyxiOX8ig0IoOg+yL8cqQsnMiVgSYJ0Ztra4uWzMZwrUWuNXePY1Ge5DHOcssiUra65E31a7D3YFWixFOz3bjnn2kaTXntZ9v0GsPKUNB+C5D2r2Y1B3DnG7H5dAE+107V3F+0oozVhlsUoIvdE5BMi8svx+9313WaEpvJKq7RyE48EdX1y2deN82Kxf1KyrFMwT+eaY75qFBwuyg77Y+0I1/xaMW4SrSYWY8J3Y32RHwJ+Kv6+OSwI6EUiZadNl36XX8um4i5DQ3hsUbOX3XpUsq4p2+teIlMe7DpfK8Z1h2EtYjF+G2ucF3ktJENo67J6XytHi7/01auBas2GsPRDbINLBe7Z6kNifABU1caGk+56HuyvM89iXM+8yFZuaqF2w3VoMtGTqWpXiWfF6OiKJf2JJEsAJm+wKLpu59429c7vqoG9GxZjGzfti/T37qIvvwgB3GtvUT1+0txLnuFu3UJ6OTocEZKwXN3J02qGTRoTYOTFVXf1LNDur+dX9nKjbJelCQ2XRp7UoEaZPoziLsMRjMerva9l2OE0xUIWI2vqi5SDA8YvHUNQDk4PkbOz5rksQ06O0H6OhICkeY5tb7bII7ibhJrvArPXzDLLdU1dI2NOQ5LUXl7TujcF2ZIHu/IvVNXvUdX3qeoHgO8A/nlkMa6lL1JHI/r3h/QfDE2qafZ5J6j3cDDA3bmNOzluhiesAiE0Xzc4TrIMd/cO7oXn4YV7VC/cRu/e6i6HYCo/FyPkbFiPvVk7VCBc82vFeJo82FrmRVaPHuPipz7MUlnSOJbcE24dwq1Dm6ddhZUNWq89opNLM0dzx/VyqpeeozrKKY9ziiNH/1Gf/v13Ohr+Op5QTd+2nzdANGwuvLlLtXFT6YCfJrIYV90X2VwkLJ+pGGWW3LRERcz/bmt7tAASrMHDFQE/EVwRuglX6CZ5Nzle+VkwsI0gcb8WzRoajXFv3G8UDLNs5U2rVqe8+o1Pg+hTfKeTKfLmO+RZRhLe1GlBOL9YTCrc9Adjb2AzWPQGaEW4iFLhWdbM4141rrMpmD2mqtCzc8tpTqdbEDi5BHvC4c2hQW0S2ybKLDfBZU28W8S2dpHPrIHVlOUti4osxCZjq+tib2A3xKY7iG6CHdp4JOw92B7rxT4GuyGWif/uBDbcNXQVlP0SeWOIWzhIYRdg3d47FhzuDexmEDffEbQzqOZlDbaNa6T2rneeJYM5lh2/mwa2jPkgDnd0aEXirDUiZjLt6qxuGnHepBwMTAQvapdpEsxbdm+bbr5dDRYO5lDVf7vo4N0zsAW9ivVT/T7y0gvoQa82QhlNkQcPN9KCvxCxcUPyHO7eRg96lMd9yuOc/MkU/4tlQyuawaZmFYmubhd5yWCOhdi5KFmyHHf7Fu72rVpntX7Oe7TfIwxyNI/KzVtWb+4gea7cEXIhZG6OTYH3Njmk19vsJuX6lOnn0+CM+PVnZ0+1ZDDHQuycB3Mf/Ere+eYXkKDc++R99NdeaZ7MM6rbA6a3cgZvDnFvPUSLwkbEVFVXNHgWreFZ7zp2a58jXSsRG0OAh6e4LCOr7gIH+FExt3S74yN4z/Mggn/tTarTx+/uXm6K63uwyxQO7VQLBnOo6n9cdOzOGVj53BGnXweudNz53FH3SWcD3stDZ8OkHp12CYf+Cm/WZr6+WxLiomtFvVU9v0DBxg8e5shkvqFDej2m94w63ntnAGzGwNaRaJ0ZzPFsGFj24Ix7P39gtJdH53WuXoOi4zH9+xdkF33ck+FcHl9VjRM/i0XG1D7OuWapbVGfr0Q6zrk6FgPQ8yHeO2Q8JUynndqkjkbk989sKR2OFp119dCV7iJfAIpoXGkwx/++7PjrdhV9ATjDKn+lqv6Wdc2LDL/2CnfefGA/R15YGn0cLoa4X/0imfc1GbFu0Y9YOMNxZnDo7DESd4E4t1SHv8aCc6Tzp17I6p2HyKNH9X23UT1+Ug+p3+jGZHUebOFgjmUH38SD/W5Vfbv1e+qLXOm8SC1L9GIIMN9bGOdF4lvarNfYhV11TD0e+TpU6eucf5F2WBszE3XXPYjB7mlFp1kymGMZnmaJXF9fZEiztGeC8ZaX6KClBfZuqDKmJN0y2tY5OrIE7R7HBd3k1k53iWgd4I6PkZdegMzHzYEiT84p33xrrSmLXS92K/DPxGRm/kZsRVtPX2R9RW0kj2afmg2cXTPC5d2k0OeVChvdMZsBPn/tazX1ips7uWQZethHe1ltUH5aLDz2SwHXNbAPq+rr0Yg+ISK/eN0LPNW8yJaoG0AtANeSaGrr5du8xxgjJUNYZqBRTnNOOKWqrqm7Oi+gstDAZ9QXdTLBnZ5Hgbo4EOJiuH4O2S57MFV9PX6/LyI/CvxWNjEvcpHydJKdrACqWgsCkRigx6UsvfuLhFFCGZ+SbjzkXByZfEU81J5HGX/v3DN0E8At76TTKeHhqT1eFFH7f/06+avaRd4UVyaDRORIRE7Sz8AfwHIea+mLxHncYIAbDJo5iyLGwU+SAp25QDP/ufbI40XPtec5xsfqNzkJkXQMJp5nNv5qo/17qi60xVTS48RegsMDq6keHVqnd7aBbJFe82vFuM5f9h4sW5uO/7uq+uMi8jOsoS/SHQxw96IS1Olj6ymM5RWC7SLrgFwD7Y2AViFKnrfyTlVlatLtzUFcGm1Js+bX+hWLRgW2upzmNF3bBquhOSYaVn2v4hAvyMkJ4X0voN7hzyfIeApnF9bXuSYvJuxwkK+qnwe+ccHj6+mLdM6mkRFzXPG75bvcYs9kN8SiKSH2XGiOATodtdd9U9sT1K6DuOtMhpxiQfEOzW2KnGYOvNuMhMCuGtimoeMJ+ugUghJG42aq7XRqycyZAFyrCqbTZilyjdx5fUztZaT7mIbLjWaBTtichPqiv6E9R0kckks9dCE8OSP79ahRMbEZkmuXcFohm+Km2D0DK6ZUp93Oblv25g2n3pnFJaye3d3efabhC7NxVcIynfxFjNm0eVjkKUVaj8VEqzhcLy6V8X7DcEgYDq/zr1gtthTk74aBiVgCsr+AwqKx/T7JXI7G84Ou0vzIXs/yTGWJjifdcyzDsufi4PjmGt3luTOyOczsAttLcdTTcAcHNng0bSjS5gLbWeqKtDWW4cvbg4lD3vM81fMLpvQGxQ0L3NkFMi3MeCYzycsYVMvhAXJwAOOJfT3lsrO0fCPSGPN0ShhfkiD13vTCnrtLuH0IZcBNS9MsK80w5eyC6v6DtS+T28BOGJh4T/X8CefvP5h/TiG7CPQe9/EXU9zFkKr9aW+NLpZ+Hz05tF2Tk/XmLlu5s6WHpPgt2AZBvTQjbNI5NoEv964id+uYX/nDR3zTf/45+z3+N0Lc7r16dofXH50QXr/Nh34I+A9nnddrUFwmlO9/nrMPHHL41jG9x0+ozi+a5bQ9NXcmlZGevxKzsdfMrKQ5IqPExG1R4M6HeFXzwmfntkymvoLR+lUOv6yXSOn1yL7ujI994CcB8K04rNLAz00L/r/R1/Lx57+J4u5zsyXw+o2e3u4xfNHhpzm93GTFVS34Fmdy48bbWlAjvIbGhW0cUg6u9Y61R8y0z5NGKAM6HluqZTymevQYNOAODupldu3c/C9nA9OLIb1/eYuvG/2ZxQck6v0bfT706PH8hkhs3PLgjXPuygm9h2PLQWVZTEdUjVgKcflkScK0c1q59HmrDLTycipWoWp9QFIOrN50FGX0dDH5G8JGlHi2VSraCQOrnjzhvT/w7/mKj/XnnxQhfNVLnH3wmHxY4d56aAY2W/OrKvilLzD4Vd8wVPMcqQIad3/15DPvLTeVEBbTgGont2g0YIgDHdozLeus/UyuLihMJrZTTNpmxLrkJmSevtxjMMB0vy4uLFbKs04B2o0K8mHAj1qf9k7eyaCVxUgigub5/EVSTqxz4QVL3WU1x2WUms4w0hZTIxXTE9r0bJmhaq8JAlsaJLNDBpbgjg5xLzyH5hnhsG8tYD2PHwf8OA6En0mcJspNmqCBqs2VhIV5LpuwQef1nYGncZmbC9o11FPLahasnaC5VqukJL0e7vYtC+ZTzqsoYDS21+dZPVheiyWyoavCl7sHS5AsIxwO0NxTHfcIPWeljkqR2fHHbU8AdRtZoj3rokB8FjHrXy+ZybhaE0TqYQ4V9fLX9kqdwRAdmo7AoI/m0cCSd51MoNIOZXrd+LLeRc5CQoAC674JAT3sU570kZC4U2Gxhwmt2mKex0kbcXloTZq9tIMoy7rLife4O7fRgz4uDQ8NAX1y3sm+mxeNy3LK9h8fMfnAc5SHHj8K+EllubyytNrqaGzeaxPNH3sDi4jLjBQlvPWAcH6Bf+k9uIPnzIPFN1l8O+0wU5/0Pi6XsV+RGNhnUU8/y8yjhYYDVrfxp9fGx8R79M4JxV1LAouCFBVZWaGTiXmkNBA1GnAaciWHB5y/3GN6IvTOlPwio3/q6D/uI5HJuvalEdgm4XD3DAwa/fmjI5z36KAl9psMRdzyT34IZlxp2QpqtDHXDqhDNKJ24Tsthe2eyZhHSx5Atfak9fVdmvLm6qJ7OlwqcKV9F1UbJDHomZc8O2dj2GUPFlvEvx/4BuxW/zTwOdbQF1nPsHae6sW74AVNvYgiuH4PKQ+sQBxpOq6Xz5Vs6ibYtAQ5QULWdAehtRq0pJSGiNU6kwxBnscNhSJVQCpbHmViij46ndY6E7OzkOxDIOTDADiyiSIlqBfKu4e4SQ93dm475w1gZ0fJRHwU+HFV/Q0Y+fCzrGtepBPwDvVCOMgoj3uEfoZ6Z+36eY708uvpSyzrcWyTFiOrIc2ZTI8BzWNBG+Oq4u6xvYtNY/zaQ9/zDHWCKxVfKK5I5zAvZt5uk+In1/xaMa4zbe0W8DuBPwmgqlNgKiLfxjr6Ivt9yjuHhIFnejuj7DuLeyrFT3sMMocbneAensE7Dxefox28u6bAbHGWa4aWeg894u4zpRbMgGsERcaT7icxBMh7RjFKqRFo0hPHR+jRAWGQ4ccBKZX8vMSNCzPQorKqwpopOm3s8i7yg8AD4G+KyDcCn8am316rL/KmkDynOswoDz2TW55yEAPrSnAFSNknG2T0piVylncD9KUnTXkqrcXigDgZd+Zf4L0VvjXYOD6AydTirrQLVTVa92Awfy0n6OHANgUiuDLgpuDPxrjzcXMPZbWZLD6YZ1qdNsX7gR8CXopn/T5V/eiy469jYBnwTcB3qeonReSj3GA5vGnjrRYF2YV90tVBNhJCLlQ9W27yYYk/nyKTwnZvLlhTx2zKIRlcokZX1Px+6jE0LcariwXPKkAxaZa6NNi0LCE0BqZl1epAj0iTb8+G5DMcfhmOYTxp1IDSXO8NYMVNHyXwF1X1M7Hb7NMi8glV/YVFB1/HwF4FXm2JjH0cM7Br9UXeuPF2NMa/fYb3nvy0h+aeyb0+oxcyfKHkbw+R0zMrGbkoL+BpdnYtjYlakCRRdsrS8lz+EMlzVCqIaQKJpaUwGRLOz3H9Pu7ObTOwskKnk8YI07nibteaaG3TQFUZJfrNq/9ULTc0zg9WFl/FVSutXGci8lngZeDdGZiqvikiXxSRr1fVz2GdRL8Qv/4ENtbvT7CqeZFVhZtEgRNVJPNkfU927MlGAaaFCaCkHFb3Xptu6SVLpsbeSA3BzhEUXCMs1369ViHmxOK1Khp1nSoKpairNSaSgWtRLqV1bwtyWQjRxfMi8qnW798XncT8OUU+gAmhPLXC4XcBHxORHvB54E8RpXtW3RcZRmN48HYz4F2E/EGf3hcHUMaBU7HDyF6Q2AyNPkUtB9B6HKi9Tzi/6PC2RATicFBLaUSPdGGCcpoacmfPF1GnyK7qOJp9fFMiwDfbIV6pcAggIsfADwN/QVWfLDvuutIBPwssuujq+yJDdTnHPXmCRe369a8yL07S8iB19nwmKVof6sRipKRPVjf6XvEuXUaBTjHZlgY0rHIXGdWlfxj4mKr+yGXH7mYm/yq0a5FtyaZZSrM4Fk7LSm90u7U/ianEUpFWC4LwtoEsotlcZTxbnP6xQoVDAX4A+Kyq/rWrjt85lekrkVgJkbgnmVFe2u36XU2I5k8UJw1LNSn1pJgq7TYTh8tHL5liqVprYlZ7YltMqxtidYnWDwN/HPg9IvKz8esPLTt49zyYCJLZjm6uOSNp6MemWMkySyFUlX1EW3ytjpdrZ/0XkQ5nkXJd4oBg9+NSd3brduIutjbOdK/O14bcloqaw+xGYF3Q1S2Rqvqv4fr8xZ0zsDo9AITTx9ZWHyHe1wQ+cc6+pybbqjLDKMrWslnhj4+g34diGnefLUNIO7+2XlhaIqGmNrt+H3q51R+jllftNaHRaI3XdoM+cnhgu9HUn1m31yVVn2Dy65tMtm4BO2dgVsfLm5/bn3zvzajyvMmqA/iCWmfCxcaL5DW8RzKPBg9Smq59IiSyfDqttIc8RC/YFimZfV5i91D9e5bFNIdvap3Q7F5bNJ91Y8WJ1hth5wxMJxP08RMj+p0c196MEBrGKlg+LGbTay8Rgr2RWda09k8mTVtYW3ac6OdjMVvL0Owckyp1CvInE0usFkW9FIZpgfim3gnNkklVoeOxbRSmU/OUeeyBLKomX7bBRKtsaczz7hlYWVI9eYLkPfydW4ST2O0dbCmS8TSqMlfoxUW37SuJn7TmHYXhMGb9m24eY7jGYF6NjaclsQRUzqUl0u6yU/MMFRpLRfUGo6UJpqNIZExxZOKMVdX61XRmsSamxHWwcwZWIxLzqqN+lBg3op9ThUKQQb8OvCWWbdpep1ObjDvGepmaHcLgXL2xSNfuGAdY3NXup6xFUCL7wkVKtkQmbSwf1QXzXo4kic5pwaYFf/eM1hmI95R3Dhm/2Dc2aFBcqfRErNB9aD2UUimU1sAqT87rtnxNtcJ4rhRHoZFoWFVmHKk1rpd3i97Tguq8Ymbb2MgxxUK5DPq4o0PIMmPeeo8OclORho6RI4I7vbD73NQOMmHvweahXgiZ2OqT+IG5w4WWbKaPybwS8xCJuhxZrEliMxnSIhKipLyWEysTOQEfaumBtLNM5SsAcdGok75E8lpRuTDkMc4L0Rt6BwIu284Q1X2QvwDtOYfqhJBDeZwjBxluGnBFZeS9aWFSSM4hxzZAqw5qZ2qOqmp6r1VV70q7OSyp6c6Ok+ilBrFRRGoPp7n9LkVpS15ZIRcje101sJixqpDhBAmKHvQJvQyZbJBBkaBsrYqw0wYGdFy7OqE8cEhlNy7B2KFSRoaoMxG6Dtr6p6lhQ6RhSoiAOmuIDYk6HZs3wJ4/GKD9GKNFrVgdZITckz0ZI9OiMw5GfEwETwvk3DYZQvS0iyaVbAD7GGwBVAABPwm4if2HQt+WQTcN+LNx7cGoIvkv9T0uqhWK2LscPZGqNilpEeP819x8gaxvBpX5pg5ZVohXKF1cmqNhtdSkKStkPLHX3jZRPQ0BGU1gMl2ae1sX9nmwZRBBBfw40Ht7iGaOyfMHaC74UYF754mlFSLVxsiEWjNRbZcZl7/MNyWgRd3fqZyUWKyZh34PFWm4VFVUJSyjMZcOmVgDLQA9KylpUUBRIEeHTL/iNiFz9F9/gj48tXvdtAeL9dttYGcNTLKM6e2c0XOOkOXAISETJncyQg6uPGBQ3oEy2BKllswMKVOfpJpS+5p3ZmRVMO+S1Kvb8gCRTq2DnjXQ5lGwLgSLqVJ9sa0BlkRaEtreKQRkGnABM8oo17QN7D3YDOTkmDd/m+foP3nIg/MDyicxtvIViCJDTza8jRSQnwtS2Y5SHfQfKXc+X+AnlTW6Ri8kQXGTCv/4jHBmtGgOBh0ln/KFW0xeGCClkl2USKW4MrJXq6aSoJmzWOxwYFr3kQzJdNLsJs8vyH4lDh0tS0ttbGt4/a4amIh8PdZgm/BB4H/DOktW33ibkGdMXyz51pd/mdfGd3j17A6VSs2WGE5zJpOccuopHveQUtBMUa+EXsbRW952nplYKqNU2xwE8JE/r7M6FM56MacnHj9V/DQAxtAQHEqAkLyYA2eGJppbegSrRJinc6aBf35hObcoGmwHbf7d3lkPFnn4vwkgTjl9DfhR1jSQtL7u+QV3P53xj06/mfxc6J0CAuWBear8Ao4vFAnEtnzLOUlQemclh69c2AYgBe6qUJlGKmVpmXtVk1KCmH5w5G/2OY4pEHcxsYA9JnLrrnMR9KBvS+60sPJVWRIiW0OL0mK/dp4MupuBTUIxNZ8t4KZL5O8FflVVf31tjbcR4fEZL/2LB7zwmQP8w3P09beQXo6+773oQY5/8Jhw/21rlL19y2Km8yG6aMhBy2MoWAoiz2zyxrjb/OomE7L7/VY1IBAiu7UWlEu0oX4PHY0JZ+dNgRw6HUfSMx0KHU8Wiw9vCDvrwWbwHcDfiz+veSBpQC5GeFU4uyCMJ0hV4c+HaNFDL0Y184F8bLu38dgklWLwnhKrrXvpXmOh7mrosC/S6+tGkhCL40Vhy2t7OFf7vBqsj7I9uGuLlOmd30XGjqI/DHzPTS7wbgeSallSvfWg5scTKnQaCG/et3pffGM1KHL62F5Ta4cpohK9VZNo1RTptjqQ6r/PW4lIq9jR3aZaRy0Mrap43WABPXRUdhKrotF7rRpKzjaNi2fDg/0XwGdU9a34+9oHkmoxTdP24gNqVJd6rItiCtKzcU1ULYzdQTVm1J+bh6VTUrLWtxmdVmcZ/3gAWuicV+q8bssG1cEzQtf5TprlEZqBpCttvL02Il++RpuW3JqS1s2aWxwliXLdOV98bSqStwebJrZq8lazvZata89Jae6AoQmRdbIFXFcf7BD4/cB/33r4I6yh8fbqm3HzP7fH8rXncAP1dLQaRsmu6c+z4wFToraqCNPoGVM1oC1DUN9Dqyw1e59bqDkuww06u1eK6zbeDoHnZh57h3U03l55M5d0BbXnO9bDEcKch5G2lms67YLm2tR3maSeFC7t8N5ZPCNL5G4gxl1dNefIvUqeq61cUzNStbNstUXnbBRzK7WQdL7SvMjEB6siJTvStOfa1YRGcXqnDHBfi7wcLY594sXXLNWkDwZNbTE42zalpXOZEIpql8YyMwO8g9mUhjjExcGnhK733EE8C7vIrcHfuwvP30OKkvDmfcJwiHvuHty7jYynhDfesi6fXm58MA21+k07e14Lo6S2srJshFSCmpGkuqS2DLco48TaUJ/DMvQ5Lhl6VTUjoGdnhO8C9h5sOWQwoLh3hBuXuLcjk/TAHvPnGdz3QBElMU31UHyjiJNiJ3GpubZhUSRuvV3IdYc5tIXi2ukIjRn6KGQneW4tdJPJQimMrUNXt4sUkR8EvhW4r6rfcNXxO21gkveiPlcge2Lyk3L3NtnJMYiQPTgzFsPhIX7QN2pznpnHGUbPNfvJdQL9vnk6N2n0wCK9p9a49x4Z9LvK094j/b6xJw4O7FrTAh3bvbl+H4107CT5VCdaZ8fNbBqru+TfAv5PjOxwJXbXwERwBwPo940N+s4p9HtUL9wm9DPyN04Jv/4qrt9HnruL5lGJ2jlkOEafnHX5VzE2E++RwwP0oN8Vl4sDHVL7v/T7uFsntvxNJkZ/PjxAbx1D5imP+2juyd4+Rx+dWmfS0aEtvUWBTCbmHcfUxlm3w21i+MLsv3NFRq2q/zIKz10Lu2tg0MiCB0WD6dnLtMTFBtYOUhyFeQgR6X5o2zu7EJrc1iXqOLYJaCQImtcKUsRzJS+ZUh8Sa6DiEDHJ8s4Y5W1hH4NdgmJKuBiBE9zF0HJSTnC3bpkBTguj4bTlM9PE28ReJWbZixI9vzCqczWzXDlB8JAMYloYezXLkH7Pjn/n1A5NY2kSXXu2MuAENHV9S7OR2Eb6In7urolrS2heB8+GgbVimSqqDvqTE+TIWK61dEBVmWdKzbOtoe1NzkrMyC5rvEjxUlXZ8YnTVZboaNQVMwFblmc9Ya0h1urTTPexYQh6kyXyWhKa18VOG5iOJ3U7mjs+biTIRSzY7uXGdR8Ou0tm+lFcPVihTVdOgbz9MkPnyVtNvdGYRZwZRhUaPYs2/6som1jPe6Tfs01ATJHUMZgkFscWtppb6gXYXQNTJYytW8cdHeFOjs2L9Huojx3cqXexsHxW22MBiGjcMeb2ZicNi9GoTk908lpJeNg7m3EUs/XOeyRpUfgZTf6Y6NUqNgA7wfXuweHAlu1xS+jEe4QCLZcnf9fzv+QmS+SlEJG/hxFNnxeRV4HvVdUfWHb87hpYG1VlrWBVMGZAGv0iYn2GZdkoHbagziHaSra21aJn0ZojWU9bE9eICc++LhbFlVZclTzdtLDOpdmeSWzjsB1O/sp2kd95k+OfCQMLk0ntHdqtZvZka6C793UMlqSSFFv6tNV5DTRpCxc6AnOJyCgpTYIZRc1sTbvE1HvZMZjYHHx2Buc2RW3RsPmtYL+LvARqPt6EbqouJSfKKtnEjwUD2EOwDqC2oFw8zviDvjusVGLR26dNgph6YRXmA/kY3OvM4xuTxbw29sXuSyH9Pu7wsOu57t4m3DmCSnHjaT0nqBYkGY2tYTZhMqk9k/R6dT1SvIM8R28doV5wD88Ij06tFHR8ZNoWVbB8WAh2zqpCkqcb9PEnJ0bpvhguNa4rB3atE7veVSQi/yPw32G3+nPYpI9D1tkX2b5+r4ccHTZJUScUL95i+N6BaYY96TeJT8CPS9ypWFd3TKpqiC1rYGJwWWYbhl4OvZzi7gGaO3qTEjk/N9Hfg35tYDYrMpjxOhclABTyHnJ8GJtFmq6hepB8W3B4w6JzbWyLcHilCq2IvAz8D8BvicVNj3UXrWcg6fIb6XowERBQJ1Q9R+hH6cxKkTI0ijvJs6WidVsQLn1XxY8KsvOYrO33zQBTctQ7M7ik7eXsd+n3LQHby9F+D3frBHfnNu74yLxknKK7E0ix4lVfK8Z1l8gMOBCRAvNcr2PdRb8rPr/yvsg2jEKTdnlWktFoXOqgPPTWcPtYbQBoUZkgSVU1ja9JkDflsSCmJcxD+YfnTeno8CAaa5Qn7/fRXoYbxcx9FSzPBab7ddgHB5JUFy/GNr5vWtiSuu2YTFm+e14zrvRgqvoa8H9gvPs3gMeq+s+Y6YsElvZFisinRORTBSua8Jp0Jkqt50pqnO2tQuN5qlYyFGqSovR6MOg3skxAGhKKRjXC+Fg9qzt9wmfLSy3Zp9lgf3dwTe+1DQ8mIneBbwO+GjgF/qGI/LHrXuDd9kV2kGZcq1reqwpkb59zVASqw4zhiz2qTPAHDldluElpbIbp1LqqAbzDnZwg/R7lV71IcZzjxxV+VJiI3cUo9kNKVOGp4tAFRYoSP4w0oIvYOZ52nYOevSYE3PnI4rzJ1Cg8UWpgrvC+DezwLvL3Ab+mqg8ARORHgN/OBvoi21Bn8VVqjJXh2BSn9RCkT8gkftlSmvJeUlUQMktp9DwcDJje7jG95emdOVwZajKeVgFJYtMxYVo31ZZ5zYAFawXDu1oeigAynsby1rTRLNtSiaYDpYlFN4zrGNgrwH8WW9dGWCfRp4ALNtQXqVWwIaVOkH4P+uY1ZDzFiXD0qke9w49Lo/OcjSzlWVNsUt1RYDxh8MY5vceN1KbGGqcAOuih8fyMcqSKM7xj/VMGA1tGU5NIVdks7hTjFZGlkZKwS4SHN4vtNaFcR13nkyLyceAz2Lzmf4ctecdsqi+yKNDxxEbCHB9C5o1UeDGE4ZDs/tt2XKLGpNnZaYxLUdbqOVpN4Ndew4WAe/F5yhdugRdLV3iHHvYJBznOC240aJZN56DfIxybwK87uzBvNZlaTiwEdDhqSWlKx8is22ht/6GrscNLJKr6vcD3zjw8YVN9kXHeEN5b+qGKATlEhmhp/0DfyrancTGdzL50d1NlhZsmaabIdHCOkHuTftK4JDupGRR10Via9raUfO1gpxo+2Nou8tnI5A/66J0TM6qHp+hobAyJvAeEpn7oIpO0lyPHRzWfvk4/iINMcP24PBYF7q2H9nNSJcwdxUmGVAE3Hlu6IwbpUgWcc0ZCFIFBHyZT+9JmTpJOp3GWURxfswu7y132YFuH86i3XBdFSRhPjELTl5pW3aElO1dr3dc1SO8jCyOmICSK9Y7HVtvs9wDbGaaUB9Di9Df8fVFXpzikjMXs2LTbwS4UuRP2BrYcenGBPBBrD3v5PbjcI2VomBMvvwgihIPclrcqoKPSdFmnls0Px32KOwPcNND7wgPC6eOGDChSszGytx5zeD5Fc0/59e+35pNpBWVcBscFOsg5+7rbjO84jl8rOfrcA9txTqYQKsvm5zkUBeHJWe1ht/cP1PklfEN4JgwsDC2Y93duU3zNi0zu5ORnFdn5lNDzTO/0CD1heuwoDyC/UA7vlzYNZFohQZnc63P+FRn5ULn3ihilJkJFLOPsnAXugHv5JR7/xhcoB0I+VPxUyZ+U9N+8oDrqc/q1nuH7KqpBzuD+EW44NdZFEdDjQ8KdI5PgPH28lS6iOew92NXQoiR7PIHKaoduOMUlqfHc4SeecuDIRoHsrGhGzVRKzzsO+g4/Ns8neWsiiOvy9wGkKOk/qsgGjmxY4acBf1EgkyneweEbh0jwHN4PuPOJNZEUUf91PMWdeWQ8IWwp/zSHvYFdjTAc4n7x82SptT8OKc3jpLRepN90GAwxbsre6nHrC3EW+LTA3TrunnyW0zUccfQf3zDjqzuQonE64YWHZ9YfEDVa2yOViR1Q2+qBnIfud5HXQkpwtj+NRdm8sWkSmve2y4Sm2aMK1trvxCg2eb7oCh3jtEm6qYeyaeogKOHRaT25tkOJDgq6PbHfhVDQXU207hLc0RF87VdSnvRNLUYVfz5BvvC6xWlpPmQ/Q46P4y6vtM5wSe1nDk2F7llo5NMXpenA9m3QaHXUI/Si8Sq4aYV75S3Cw9PmteKaQfFt6YCd8GDsdKloZyC9HqP3HjO+560DTGHwKGfweg5DGi8Tp+XiHFLE5KyPqQXvCbHtfw5BccMYh/VywlGfkDumd3tMjy1NIgGyceDofr+bhnBiXtN7KGKerKo230G0CO0uqA3jmTIwNNiUtamjOLAdoys9B3nenRfU7xPuHNkUjjLYAIaqQiZxftGkssfFcl41QmR+pipAUeFUyYaVKRJ4oerPNti6esakHB9br6bGKW7jCfr2w93wYvsg/2poFfDjkmzoGN91DF+y7PztQb8WqAOQo0PGLx5QDRyoGU3+pKL/1nk9ak+qgGbNoPgOvFGDZFoihZBPS3LnKE/6jF/odXlfUV5TBn3Cc7eojmLfZiZkTybIk7OdMDDde7BrIATcqCTrebKhJxsKfsy8+w8BVyo6jZ9aAZeIgqo1V18qQd1l//gY0wUBCbiJx48yXBEWxzSRrm0vFft528sjwL6r6HoIozH+C6/T6+U898ohd48HyGhKeHha63sB6MNTDn6hqgN6Mg9F2Yz9qzlc8826CxE9VnYxInsLKCvCk7P6HBpKCEPkjft411qqiynVeEUs3qfBvth9TYSK6lFsXJKW8vPsYTHzj/fWgDEYLD6f6o1a6nU8JpxfdMsuaWpIGdDHT+znbbaoLYDC1kpVz5aBzeBKpZpQGhuiWE3ThZZlx1M+M9DVEQ5F5FuAj2LdZd+vqh+57Phn18D0en2GYdH0tS9DrEI2Ko5z/L+woRyvAj8jIj+mqr+w7DVXdhXt8SWCNBbnqq/L8VuBX1HVz6vqFPj7WEPQUmzUg53x6O2f1I9fAG9v8rorxvNs9v6/6mlPcMajn/hJ/fjz1zx8cInC4cvAF1vPvQp882Un26iBqeoLIvKpVSrobRrP4v2r6res6FSLqLmXrr37JXKPm+BV4P2t39+Hdfkvxd7A9rgJfgb4kIh8dRxQ+x3YWMel2MYu8l0rFu8InvX7f9dQ1VJE/hzwE1ia4gdV9ecve43oDiUE9/jSw36J3GOt2BvYHmvFRg1MRL5FRD4nIr8iIusVrFsBROT9IvIvROSzIvLzIvLn4+P3ROQTIvLL8fvdbd/rrmJjMVgsM/wSrTID8J2XlRm2jaga9F5V/YyInACfBr4d+JPAQ1X9SPyg3FXVtYjvPevYpAe7cZlh21DVN1T1M/HnM+CzWDb72zBVR+L3b9/KDT4D2KSBLSozvLzB6z8V4gi73wx8kmuqO+6xWQO7cZlhVyAix8APA39BVZ9s+36eJWzSwG5cZtgFiEiOGdfHVPVH4sNvxfgsxWlrVXd8lrFJA7txmWHbEBMA+wHgs6r611pP/Rim6ghrVnd81rHRTL6I/CHgr9OUGf7qxi7+LiAivwP4V9jwiUSW+stYHPYPgK8kqjuq6sOt3OSOY18q2mOt2Gfy91gr9ga2x1qxN7A91oq9ge2xVuwNbI+1Ym9ge6wVewPbY634/wHbRzZ/puLGMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sino_astra = np.reshape(W*phantom.array.ravel(),ag.shape)\n",
    "plt.figure()\n",
    "plt.imshow(np.abs(sino.array-sino_astra))\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7865cce2-35f8-4da0-bcba-777e9980494a",
   "metadata": {},
   "source": [
    "## Regularisation parameter for TV regularisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9ff5f899-4f98-4f55-9d8f-d4690a10125b",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e1c6c0f-aa9f-4bb5-8e63-b4e8a1ec1bc1",
   "metadata": {},
   "source": [
    "## Solve with cvxpy\n",
    "\n",
    "To solve with cvxpy, we need a matrix representation of the ProjectionOperator.\n",
    "\n",
    "An alternative is to use [ProxImal](https://github.com/comp-imaging/ProxImaL) library as a matrix-free option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0c61b022-17d0-4d23-b9f8-a69beae9b29b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                    v1.1.15                                    \n",
      "===============================================================================\n",
      "(CVXPY) Dec 15 08:21:24 AM: Your problem has 1024 variables, 1 constraints, and 0 parameters.\n",
      "(CVXPY) Dec 15 08:21:24 AM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Dec 15 08:21:24 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Dec 15 08:21:24 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 15 08:21:24 AM: Compiling problem (target solver=SCS).\n",
      "(CVXPY) Dec 15 08:21:24 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS\n",
      "(CVXPY) Dec 15 08:21:24 AM: Applying reduction Dcp2Cone\n",
      "(CVXPY) Dec 15 08:21:24 AM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Dec 15 08:21:24 AM: Applying reduction ConeMatrixStuffing\n",
      "(CVXPY) Dec 15 08:21:24 AM: Applying reduction SCS\n",
      "(CVXPY) Dec 15 08:21:24 AM: Finished problem compilation (took 1.504e-01 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 15 08:21:24 AM: Invoking solver SCS  to obtain a solution.\n",
      "----------------------------------------------------------------------------\n",
      "\tSCS v2.1.4 - Splitting Conic Solver\n",
      "\t(c) Brendan O'Donoghue, Stanford University, 2012\n",
      "----------------------------------------------------------------------------\n",
      "Lin-sys: sparse-direct, nnz in A = 115831\n",
      "eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00\n",
      "acceleration_lookback = 10, rho_x = 1.00e-03\n",
      "Variables n = 1986, constraints m = 6789\n",
      "Cones:\tlinear vars: 1024\n",
      "\tsoc vars: 5765, soc blks: 962\n",
      "Setup time: 2.34e-01s\n",
      "----------------------------------------------------------------------------\n",
      " Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)\n",
      "----------------------------------------------------------------------------\n",
      "     0| 2.26e+20  3.30e+20  1.00e+00 -1.92e+26  4.81e+23  1.16e+26  5.37e-03 \n",
      "   100| 3.31e-04  4.17e-03  2.98e-04  6.63e+03  6.63e+03  2.60e-12  6.95e-01 \n",
      "   200| 1.23e-04  1.26e-03  2.56e-04  6.61e+03  6.61e+03  3.25e-12  1.34e+00 \n",
      "   300| 4.13e-05  6.12e-04  1.30e-04  6.62e+03  6.62e+03  7.00e-12  1.97e+00 \n",
      "   400| 1.48e-05  1.89e-04  1.23e-04  6.63e+03  6.63e+03  3.57e-11  2.61e+00 \n",
      "   460| 9.37e-06  9.84e-05  2.26e-05  6.63e+03  6.63e+03  1.68e-12  3.10e+00 \n",
      "----------------------------------------------------------------------------\n",
      "Status: Solved\n",
      "Timing: Solve time: 3.10e+00s\n",
      "\tLin-sys: nnz in L factor: 648247, avg solve time: 3.19e-03s\n",
      "\tCones: avg projection time: 3.22e-05s\n",
      "\tAcceleration: avg step time: 3.21e-03s\n",
      "----------------------------------------------------------------------------\n",
      "Error metrics:\n",
      "dist(s, K) = 2.2737e-13, dist(y, K*) = 8.5265e-14, s'y/|s||y| = -1.0523e-17\n",
      "primal res: |Ax + s - b|_2 / (1 + |b|_2) = 9.3725e-06\n",
      "dual res:   |A'y + c|_2 / (1 + |c|_2) = 9.8413e-05\n",
      "rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 2.2632e-05\n",
      "----------------------------------------------------------------------------\n",
      "c'x = 6634.5068, -b'y = 6634.2065\n",
      "============================================================================\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 15 08:21:28 AM: Problem status: optimal\n",
      "(CVXPY) Dec 15 08:21:28 AM: Optimal value: 6.642e+03\n",
      "(CVXPY) Dec 15 08:21:28 AM: Compilation took 1.504e-01 seconds\n",
      "(CVXPY) Dec 15 08:21:28 AM: Solver (including time spent in interface) took 3.352e+00 seconds\n"
     ]
    }
   ],
   "source": [
    "N, M = ig.shape\n",
    "u_cvx = Variable(N*M)\n",
    "fidelity = sum_squares(W@u_cvx.flatten() - sino.array.flatten())\n",
    "constraints = [u_cvx>=0]\n",
    "\n",
    "regulariser = (alpha) * tv(reshape(u_cvx, (N,M))) \n",
    "obj =  Minimize( regulariser +  fidelity)\n",
    "prob = Problem(obj, constraints)\n",
    "\n",
    "res = prob.solve(verbose = True, solver = SCS, eps=1e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cef55a5-eb40-432a-b45e-42f809c0d7a9",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Solve with CIL using the PDHG algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9329f8c7-9462-439e-b864-e8a2f0dc75df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialised GradientOperator with C backend running with  20  threads\n",
      "PDHG setting up\n",
      "PDHG configured\n",
      "     Iter   Max Iter     Time/Iter        Primal          Dual     Primal-Dual\n",
      "                               [s]     Objective     Objective             Gap\n",
      "        0       3000         0.000    3.25987e+05  -0.00000e+00    3.25987e+05\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/cil_devel_epaps/lib/python3.7/site-packages/cil/framework/framework.py:2421: RuntimeWarning: invalid value encountered in true_divide\n",
      "  pwop(self.as_array(), x2.as_array(), *args, **kwargs )\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      500       3000         0.007    6.88734e+03   5.73411e+03    1.15323e+03\n",
      "     1000       3000         0.008    6.71964e+03   6.36772e+03    3.51920e+02\n",
      "     1500       3000         0.008    6.68013e+03   6.53405e+03    1.46077e+02\n",
      "     2000       3000         0.008    6.65576e+03   6.58625e+03    6.95139e+01\n",
      "     2500       3000         0.008    6.64719e+03   6.61190e+03    3.52965e+01\n",
      "     3000       3000         0.008    6.64371e+03   6.62034e+03    2.33727e+01\n",
      "----------------------------------------------------------------------------\n",
      "     3000       3000         0.008    6.64371e+03   6.62034e+03    2.33727e+01\n",
      "Stop criterion has been reached.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Setup PDHG\n",
    "\n",
    "# separable sum\n",
    "f1 = L2NormSquared(b=sino)\n",
    "f2 = alpha * MixedL21Norm()\n",
    "f = BlockFunction(f1, f2)\n",
    "\n",
    "# non-negativity constraint\n",
    "g = IndicatorBox(lower=0.0)\n",
    "Grad = GradientOperator(ig)\n",
    "\n",
    "# BlockOperator\n",
    "K = BlockOperator(A, Grad)\n",
    "\n",
    "# stepsize\n",
    "sigma = 10\n",
    "\n",
    "# Run PDHG\n",
    "pdhg = PDHG(f = f, g = g, operator = K, update_objective_interval=500, max_iteration=3000, sigma=sigma)\n",
    "pdhg.run(verbose=2) \n",
    "u_cil = pdhg.solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2fe5efd2-c4f9-436e-87c9-c36797665812",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CVX objective = 6641.616740187816\n",
      "CIL objective = 6643.71435546875\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmGElEQVR4nO3deXCk9X3n8fe3Wy21bqm7Jc0FaJgMx4A5zHDYG5sxdghHxRNwCExcxFDrYqlgb8C1a7NZeyFx2ZUQHzgxNou9VJyU8SyLIRBnHK9hjbFjCDM43Ic94AHE6Gzdt9T92z+elkbSaEY9M63+9fF5VQl1P/2o9W16+qOnn36+z9ecc4iISPEL+S5ARERyQ4EuIlIiFOgiIiVCgS4iUiIU6CIiJaLC1y9OJBKuvb3d168XESlKzzzzTJ9zrmW527wFent7O3v27PH160VEipKZvXmo21bc5WJm95pZj5m9eIjbzcz+xsz2mtnzZvbuYylWRESOTjb70P8OuOQwt18KbM583QB889jLEhGRI7VioDvnngD6D7PKduDvXeApoMnM1uaqQBERyU4u9qGvB95ecL0js6xz6YpmdgPBVjzHH398Dn61iJSKmZkZOjo6mJyc9F1KQYhGo2zYsIFIJJL1z+Qi0G2ZZcueIMY5dw9wD8DWrVt1EhkRmdfR0UF9fT3t7e2YLRcr5cM5RzKZpKOjg40bN2b9c7k4Dr0DOG7B9Q3A/hzcr4iUkcnJSeLxeNmHOYCZEY/Hj/jdSi4C/RHgjzNHu1wADDnnDtrdIiKyEoX5AUfz/2LFXS5m9j1gG5Awsw7gNiAC4Jy7G9gFXAbsBcaB64+4Cil7ycEhnnvqMZidIJyewVLThNPThNLThNIziy6bS+Eye/rShMAM5wxnFizPfA/WMZyFliw/cD2NYXN7CJ0DHJb5vujygh2EbtELLbjs5vc8Hrwncbl9izZ/2uqFd5z5nQuWtZ15KZvPvCDb/41S5lYMdOfcjhVud8BNOatIytLL3/8CF739P32XUXBeeOencObjvssoK11dXdx8883s3r2bqqoq2tvbufPOO7nyyit58cUXefzxx/nSl77ED37wA9+lHsRbp6jIQlVD++ijmdA1/wDhSlxFJYSrIFyFC1diFQe+Y6Hg7aib2w5PZy4DLk2wRZ3GzAVbvc4Bacw5LHO7ZdbBBdvoWLA1j819hQDDMt8xC1ZZNBBm4Zb9wmUHv1Ve7t3z/H1nVli84W+8dNdHiY2/cVT/P+XoOOe44oor+NjHPsbOnTsBePbZZ+nu7vZcWXYU6FIQqie76Yus5ZRT3ue7lIIxWbeBxOhuXDqNhXQevXz4yU9+QiQS4cYbb5xfdtZZZ7Fv3z5/RR0BBboUhMaZHvrqTvZdRkGx+rXUdE8xNNhHY6zVdzl59ef/9BIv7x/O6X1uWdfAbb932mHXefHFFznnnHNy+nvzSX/2xTuXTpNIJ5mpXeO7lIJSEdsAQHL/bzxXIsVCW+ji3UB/LzGbwhrW+S6loNQmgm7qkd63gfP9FpNnK21Jr5bTTjuNBx54wMvvzgVtoYt3yc59AERixx1+xTLTtKYdgMn+tw+/ouTMRRddxNTUFN/61rfml+3evZs33zzkGWsLigJdvBvtCV4sdQmd32eh+JrjSTsjPfiO71LKhpnx0EMP8eMf/5hNmzZx2mmncfvtt7Nu3eJ3j4899hgbNmyY/3ryySc9VbyYdrmId3NboM1r2/0WUmAilVH6rJHwqM6kkU/r1q3j/vvvP2j5iy8GIyG2bdvGxMREvsvKirbQxTs39A4pZzS3apfLUgMVCaITxXEMtPinQBfvwqOd9Fsz4Uil71IKzlhlKw3TPb7LkCKhQBfvqie6GahI+C6jIE3VrCGWTvouQ4qEAl28a5jpYayqzXcZBcnVr6eBMUZHBn2XIkVAgS5eOeeIp/qYrtXUwuVUNK8HILl/n99CpCgo0MWr4cEB6m0C1FS0rOp4cCjncE9xHActfinQxatkV9DWHmne4LmSwtTYFgT6RJ+ai/Klq6uLa665hk2bNrFlyxYuu+wyzIzXXntt0Xo333wzd9xxBw8++CAf/OAH55f//Oc/56yzzmJ2djbfpSvQxa/h7mDLszahQxaXk1jXDkBqqMNvIWVi7vS527Zt4/XXX+fll1/mi1/8Itu2bZs/nS5AOp3mgQce4Oqrr+bKK68kGo1y3333MTs7y5/8yZ/wjW98g4qK/Lf5KNDFq8lksOXZmGlzl8WiNfUMUkdoRFMd8+FQp8/92te+tijQn3jiCdrb2znhhBMA+Nu//Vs++9nPctttt3Huuefy3ve+F4Dbb7+da6+9losuuojNmzfPn1Lg2muv5eGHH56/v49+9KM88sgjx1y/OkXFq/RgsOUZX3OC50oKV384QdV4mTUX/fBW6Hoht/e55l1w6V8edpVDnT73jDPOIBQK8dxzz3HmmWeyc+dOduw4MMztxBNP5Oqrr+brX/86r7/++qKfff7553nqqacYGxvj7LPP5vLLL+fjH/84X/3qV9m+fTtDQ0P84he/4Dvf+c4xP0RtoYtX4dFOkjQSqar2XUrBGqlspU7NRd7t2LGDnTt3Mjs7y8MPP8xVV101f1s6nebRRx+lrq7uoBN5bd++nerqahKJBB/4wAd4+umnufDCC9m7dy89PT1873vf4yMf+UhOdtFoC128ik50MRBOEPddSAGbrF7DcROv+i4jv1bYkl4thzt97o4dO7j44ou58MILOeOMM2htPTB05K677uL000/n85//PDfddBNPPvlkMCYR5r/Pmbt+7bXX8t3vfpedO3dy77335qR+baGLV/XTPYyqqeiwXN1aYgwzOTHuu5SSd6jT5/70pz9l06ZNxONxbr311kW7W7q6uvjKV77CHXfcwSWXXML69ev59re/PX/7ww8/zOTkJMlkkscff5xzzz0XgOuuu44777wTCP6Q5IICXbxSU9HKwk2Z5qLOtzxXUvpWOn3ujh07ePXVV7niiivmf+ZTn/oUn/70p2lpaQHgzjvv5Atf+AL9/f0AnHfeeVx++eVccMEFfO5zn5u/r7a2Nk499VSuv/76nNWvXS7izcjwII02RrpeTUWHE80c0jnU/RvWn3iK52pK36FOnwtwyy23cMsttyxadt999y26ftxxxy0aKn3SSSdxzz33HHRf4+Pj/PrXv160tX+stIUu3vR3BR8eRTLt7bK8htbgCKDxpJqLSsWjjz7KKaecwic/+UkaGxtzdr/aQhdvhjJNRXPt7bK82NqNAMwMqLmo2Nx+++3LLv/Qhz7EW2/lfheattDFm4lk8A+6qa3dbyEFrr4xxqirxoZLf3KRc853CQXjaP5fKNDFm9RAMCsztlZNRStJhuNUjnf5LmNVRaNRksmkQp0gzJPJJNFo9Ih+TrtcxJvQ6H4Gqaepps53KQVvJNJC7VRpNxdt2LCBjo4Oent7fZdSEKLRKBs2HNlJ6xTo4k3VeBf94QRNvgspAhPVa2gd/DffZayqSCTCxo0bfZdR1LTLRbypn+5hpLJ15RWFVN1a4m6AmZlp36VIAVOgizexVB9TNWt8l1EUrHE9YXMku3ToohyaAl28mBgfI8awmoqyFI0H+1IHu/f5LUQKWlaBbmaXmNlrZrbXzG5d5vZGM/snM3vOzF4ys9z1skpJ6usMjkGvaFJTUTbqEpnmol61/8uhrRjoZhYG7gIuBbYAO8xsy5LVbgJeds6dCWwDvmxmlTmuVUrIUGZLU01F2YmvbQdgOnOop8hystlCPw/Y65x7wzk3DewEti9ZxwH1FpwXsg7oB/I/UE+KxnhfsKXZ0KZj0LPRGG9jykVwwwp0ObRsAn09sPCTmI7MsoW+DpwK7AdeAP7UOZdeekdmdoOZ7TGzPTrWtLzNZrY042oqyoqFQvSG4kTGSru5SI5NNoFuyyxb2sr1u8CzwDrgLODrZtZw0A85d49zbqtzbuvcqSalPNnIfoapoaa+2XcpRWMo0kLtZJmNopMjkk2gdwALR7JvINgSX+h64EEX2Av8BtB5PuWQqsY76Q8lfJdRVCaibTTO6p2tHFo2gb4b2GxmGzMfdF4DLB1P/RbwQQAzawNOBt7IZaFSWuqm1FR0pGZr15JIJ0mnUr5LkQK1YqA752aBTwA/Al4B7nfOvWRmN5rZjZnVPg+818xeAB4DPuOc61utoqX4xVK9TFZr9NyRsMZ1VFqK/r5O36VIgcrqXC7OuV3AriXL7l5weT9wcW5Lk1I1NTVBzA2xt05NRUeiMhbs+Rzo3Eei7chO2iTlQZ2iknfJzrcImZuflSnZqU0Ex+yP9r7puRIpVAp0ybu5SUXR+HErrCkLxda0AzDdr8lFsjwFuuTd2FxTUauOQT8Sza3rmXFh0kNqLpLlKdAl7+ZmY8bWnei5kuISrqggac1UjOlDUVmeAl3yzob3M+ai1DfGfJdSdAYjLVRPlPbkIjl6CnTJu8hYJ33hBNhyTchyOONVrTTMKNBleQp0ybvaqR5GIjr1w9GYzjQXufRBp0oSUaBL/jXP9jJRrUlFR6VhHTU2xfBg0nclUoAU6JJXszPTxN0Aqbq1vkspSpHmoKGov1Nn1pCDKdAlr5I9HVRYmlCjmoqORm0iOHZ/uEezReVgCnTJq4HOfYCaio5WU6a5aLJfgS4HU6BLXs1NKqpTU9FRia85nrQz0oNqLpKDKdAlr+ba1udmZMqRiVRG6bdGwqNLRxKIKNAlz9zwO0y6CA3NOhf60RqoaCE6oclFcjAFuuRVZKyLvlACC+mf3tEarWylYVrNRXIwvaokr2onuxlSU9Exma5pI5bWcehyMAW65FXTTA8TUU0qOhaufj0NjDE2Mui7FCkwCnTJm3QqRcL1M6umomNS0Rwcw5/MHAIqMkeBLnnT3/sOEUthaio6JtXxYHLRcLcmF8liCnTJm7mmoqqYmoqORWNbcAz/RFLNRbKYAl3yZrQ301TUcrznSopbYl07ALODGkUniynQJW+mM+3qzWoqOibRmjoGqSM0oslFspgCXfImPbSfaRemObHOdylFrz+UoErNRbKEAl3yJjK2n75QnFA47LuUojdS1UrdlJqLZDEFuuRN9WQ3QxVqKsqFyeo1xFK9vsuQAqNAl7xpnOllXE1FOZGuX0eMYSYnxnyXIgVEgS554dJpWtJJZmrVVJQLFY3B5xDJzrc8VyKFRIEueTGY7KbKZkBNRTkRzUwuGure57cQKSgKdMmL/s7fAFCZmYkpx6ahJWguGldzkSygQJe8GMk0FdWqqSgnYus2AjA7oOYiOUCBLnkxldmSbF6z0XMlpaG+Mcaoq4ZhTS6SAxTokhfpoXeYdSHibdrlkivJcILK8S7fZUgBySrQzewSM3vNzPaa2a2HWGebmT1rZi+Z2U9zW6YUu4rRTvosRriiwncpJWMkkqBWzUWywIqBbmZh4C7gUmALsMPMtixZpwn4BvBh59xpwFW5L1WKWfVkF4NqKsqpieo1NM+quUgOyGYL/Txgr3PuDefcNLAT2L5knT8CHnTOvQXgnNNmgyzSMNPLWFSDoXNptm4tcTfA7My071KkQGQT6OuBhcdGdWSWLXQS0Gxmj5vZM2b2x8vdkZndYGZ7zGxPb6+2LMqFS6dJpPqYqVFTUS6FGtcTNkeyW4cuSiCbQLdllrkl1yuAc4DLgd8FPmdmJx30Q87d45zb6pzb2tKit9/lYngoSY1NQYPOsphL0XjwAfNA1z6/hUjByOYTqg5g4YiZDcDSY6U6gD7n3BgwZmZPAGcCv8pJlVLU+vf/hkagQk1FOVXX0g7AeK/a/yWQzRb6bmCzmW00s0rgGuCRJes8DLzPzCrMrAY4H3glt6VKsRrpCWZfqqkot+Jrgm7R6YF3PFcihWLFLXTn3KyZfQL4ERAG7nXOvWRmN2Zuv9s594qZ/QvwPJAGvu2ce3E1C5fiMZlpKmrKzMKU3GiMtzHlImouknlZHRTsnNsF7Fqy7O4l1/8a+OvclSalIj30Dmln81uUkhsWCtEbihMZ0yg6CahTVFZdeKSTpDVRWVXlu5SSMxRpoWZSo+gkoECXVRed7GJATUWrYiLaRqOaiyRDgS6rrn66h9FKNRWthpnadSTSSdKplO9SpAAo0GXVxdN9TNeu8V1GSbKGtVRaioE+7UcXBbqsstHhfuqZwNWrqWg1VMaCFhE1Fwko0GWVJfcHx6BH1FS0Kuoyx/aP9qi5SBTosoqmJkbo3v19AKoTaipaDbE17QBM9et8LpLlcegi2ZqemuTVn/8jM8/9H04d+hnn2RT7aWHtb53lu7SS1Ny6nnFXRdtr/8Bre87h5K0f8l2SeKRAl2M2OzPDy0/9kIlf3s/JAz/hDEYZopbnmi8m+u6rOf09lxCJRHyXWZLCFRX86v1/w7qf/Tfaf/ARdv/r5Wz+oy/T1KIzW5Yjc27piRPzY+vWrW7Pnj1efrccu1Qqzct7Hmdo905O6vsxrfQz7qp4ufF9VJxxFae+bztVVdW+yywbo8MDvHDfn3Fu505GrYZfn/FfOGf7fyYUDvsuTXLMzJ5xzm1d9jYFumRrcnKCV578IeMvPMKm/idYQ5JpF+bVuvNJn/4HnPL+q4jWNvgus6y98fLTTD50C1tmXuS1yClEPvxVTnzXe32XJTmkQJejNjzUz2s/fxBe+WdOHnmKBhtnwlXyq/rz4OTL2Pz+q6lpTPguUxZIp9Ls/qdv8lvP/hVNbphn2v6ALR+9g7rGmO/SJAcU6HJE+rre5vWf/W+ir/8Lp078O5U2ywANvN7820Tf9WFOes/vUVld57tMWcFQfy8vf/fTnN/3EAPWyJvv/gwtJ18AGNiCuTVmOGx+eXDZgXMYbv7yPBcsa0qspyHeltfHJAp0ydLo8AAv7LyNd79zH1U2wzvWRkfbRTSd/ftsPudDhCr0GXoxevWXT2D//ClOTv06p/c7SB0Nn3tL++nz7HCBrleokEql2PPwXWx6/su8h0H2NP4OsYv/Kxu3nMv6kFoVit0p734/qTOe4t+feIjZ8aGDtroPXF9wmcyWurHgss1fTu99jPMHdzHQ30OzjqgpGAr0MvfCk/9C1aP/nfNTe3ktcgr9l/wdW8/5gO+yJMfCFRWcfdFVObu/Z3YZPL2LwZ63FOgFRIFept5+41W6vv8Zzh17nG7i/HLrX3P2ZR/HtEUuWaiOrQdgtE/j7wqJAr3MDA0N8MLO2zl3/3dJYDzdfgNn/OHnaNPhhnIE6hPBuXkmBzT+rpAo0MvIM4/dz/E/+zS/zQD/3vQ7HPeHf8V56zf5LkuKUKwtCPTUUJfnSmQhBXoZafrXLzJjVbzx4Yc4++yLfJcjRay2vokxF8XGFOiFRDtMy0gs3cfbsfdwosJccqA/FCMy3uO7DFlAgV4mJifGaWYEV6dGEMmNkUic6FTSdxmygAK9TCQ7gwEI4UZNDpLcGK9K0DCrQC8kCvQyMdwbBHo0pslBkhuz1a3E0v2+y5AFFOhlYqwvmGhT13Kc50qkVLi6NmpsirHhAd+lSIYCvUzMDAbHC8fWaBSc5EZFwxoA+rs1/q5QKNDLxUgnUy5CQ3Or70qkRERjwecxI70dniuROQr0MhEZ6yYZiqm1X3KmLtMtOqFu0YKhV3eZqJ7qYahCgygkd5pag89jZoY6PVcicxToZaJ+po/xqhbfZUgJaYq1Mu0qYKTbdymSoUAvAy6dJp5OMlOjpiLJHQuF6LcmwuMK9EKhQC8DIyOD1NoUNOi81ZJbQxVxqid7fZchGVkFupldYmavmdleM7v1MOuda2YpM/uD3JUox2qg600AKtQlKjk2XhmnbkbdooVixUA3szBwF3ApsAXYYWZbDrHeXwE/ynWRcmyGe4Iu0Zq4ukQlt6arW2lMq7GoUGSzhX4esNc594ZzbhrYCWxfZr1PAt8HdPq1AjOZDKbK1LeqqUhyK13bSjMjTE9N+i5FyC7Q1wMLW8E6Msvmmdl64Arg7sPdkZndYGZ7zGxPb6/2u+XL7FAQ6PE1avuX3AplPpcZULdoQcgm0G2ZZW7J9TuBzzjnUoe7I+fcPc65rc65rS0tOoQuX2y0ixGqqalr8l2KlJiqpuBzmSF1ixaEbCYWdQALN+02AEtbw7YCO80MIAFcZmazzrl/zEWRcmwqx3voDyWo912IlJzaRPBmfSypYdGFIJtA3w1sNrONwDvANcAfLVzBObdx7rKZ/R3wA4V54aid6mEkoi5Ryb3GzNk7pwfV/l8IVtzl4pybBT5BcPTKK8D9zrmXzOxGM7txtQuUY9c428dEVCflktxrbl1H2hlO3aIFIash0c65XcCuJcuW/QDUOXfdsZcluZJOpYi5AfbVqktUci8SqSRpDYTGFOiFQJ2iJW6gr5NKS2HqEpVVMhiKUTWho9YKgQK9xA10B01Flc3rV1hT5OiMViaomenzXYagQC95c6Pn1CUqq2Uq2kLjrGaLFgIFeomb6g8OJ2tSl6isktmaVmJukPTsrO9Syp4CvcSlMsMHNEtUVkuovo0KSzOY7PJdStlToJe48FgnSRqJVFb5LkVKVKQp+MB9UN2i3inQS1zVRA8DYTUVyeqpiQUfuI/2KdB9U6CXuNrpPsYqFeiyehpagg/cpzQs2jsFeolrTvUxVa0uUVk9zW1B+39qWPvQfVOgl7CZ6SniDJGqU1ORrJ6a2nqGXQ02qm5R3xToJSyZOUd1WF2issoGQjEqJzTbxjcFegkb7A5miUZj6hKV1TUaiVE9pW5R3xToJWwiGRx1UJvQpCJZXRNVLdTPali0bwr0EjbVHxx10NSmpiJZXTM1rcTSA7h02ncpZU2BXsLcyH6mXZjmhPahyyqra6PaphkbGfRdSVlToJewirFu+i1GKBz2XYqUuIrGuWHRb3mupLwp0EtYdLKHoYq47zKkDESbg2HRI+oW9UqBXsIaZvoYq1JTkay++ky36ER/p+dKypsCvYQ1p5PM1CjQZfU1twaBPjOk9n+fFOglanx0iAbGSdet8V2KlIGGpgRTLgLqFvVKgV6ikl3Bh1MVTWoqktVnoRDJUDORMXWL+qRAL1HDPUHbf7W6RCVPhsNxolMaFu2TAr1ETSSDQK9vUZeo5Md4ZYK6GXWL+qRAL1GzmQ+nYmvb/RYiZWO6poWmtIZF+6RAL1UjXYy7Kurqm3xXImXC1bbSyBhTk2O+SylbCvQSFRnrIhmKYyE9xZIfc6dpHuhWc5EverWXqJqpXoYjGj0n+VOZGRY91Pu250rKlwK9RDXM9jERVVOR5E9tPGguGk+qucgXBXoJcuk0iXQ/MzVtvkuRMtKUmS06Paj2f18U6CVoeKCXKpvBNHpO8qg5sZaUM9yIhkX7okAvQf1dwei5SNM6z5VIOamIROi3JsJjav/3JatAN7NLzOw1M9trZrcuc/tHzez5zNcvzOzM3Jcq2RrJfChVE1dTkeTXUDhG1aS6RX1ZMdDNLAzcBVwKbAF2mNmWJav9BrjQOXcG8HngnlwXKtmbHHgHgMZWjZ6T/BqNxKmdVreoL9lsoZ8H7HXOveGcmwZ2AtsXruCc+4VzbiBz9SlgQ27LlCORHgwCPbZGW+iSX1PRFhpTCnRfsgn09cDCA0s7MssO5T8CPzyWouTY2GgXA9QTra71XYqUmVRtK81uiNTsrO9SylI2gW7LLHPLrmj2AYJA/8whbr/BzPaY2Z7eXu1nWy2VEz0MhjR6TvIvVL+WsDkG+3Qsug/ZBHoHsPC9+wbgoGfLzM4Avg1sd84t+57LOXePc26rc25rS0vL0dQrWaib7mWkUl2ikn/z3aI9av/3IZtA3w1sNrONZlYJXAM8snAFMzseeBC41jn3q9yXKUeicTbJpLpExYO58++PJhXoPlSstIJzbtbMPgH8CAgD9zrnXjKzGzO33w38DyAOfMPMAGadc1tXr2w5lNTsLHE3wOsaPSceNGSGRU8NaJeLDysGOoBzbhewa8myuxdc/jjw8dyWJkdjoKeDhDmsQU1Fkn/xNUGgp4bVLeqDOkVLzEB3MEu0qlmBLvkXra5liFpCGhbthQK9xIz1BUeY1iV0DLr4MRCKUTmhYdE+KNBLzFSmS7Sp7QTPlUi5GqmIU61uUS8U6CUmPdxJyhnNLdrlIn5MRltomFWg+6BALzEVo10krZmKSMR3KVKmZqtbiKUHcOm071LKjgK9xFRN9jBYoaYi8ah+DVU2w+iwttLzTYFeYuqn+xhTl6h4VNE4Nyxas0XzTYFeYprTfUxXq0tU/KnOHDI70qtu0XxToJeQqckxmhglXafRc+JPXSJoLppUt2jeKdBLSLIzeIsbbtQRLuJPc1swWGV2SN2i+aZALyHDPcEs0Whc80XEn4bGJsZdFYwq0PNNgV5CxpNBU1FdiwJd/LFQiIFQExXj6hbNNwV6CZmeGz3Xqi5R8Wu4Ik71lIbY5JsCvZQMdzLpIjTGNDxE/BqvTFA3o+PQ802BXkIqxrtJhmJYSE+r+DVd3UpTemDlFSWn9MovITWT3QyrS1QKgKtro4FxJsdHfZdSVhToJaR+JsmYRs9JAQg3BBOzBnrULZpPCvRS4RzxdJLZ6jbflYgQbQ5miw6rWzSvFOglYmSonxqbwtWrS1T8q40HzW1zh9JKfijQS0R/VzB6LqLRc1IAGluDXoiZQbX/55MCvUSM9gaBXh1TU5H4F2tZx4wL40Y0WzSfFOglYqI/2FdZ36JZouJfOBym3xoJq1s0rxToJWJ2sBOA+JrjPVciEhgKx6maVLdoPinQS4SNdjLsaqitb/RdiggAY5Vxaqf7fJdRVhToJaJyvJv+cNx3GSLzpqMtNKb6fZdRVhToJaJmqpeRiLpEpXCkattodsOkZmd8l1I2FOglonG2j0l1iUoBCTW0ETLHQI+ai/JFgV4C0qkUcTfATO0a36WIzKtsCnoihtQtmjcK9BIw2NdJxFKE1CUqBaQ6FrT/j6lbNG8U6CVgsDtoKqrMnD9DpBA0ZiZnTQ90eq6kfCjQS8BoX3BGu+qEukSlcMTagn+PqWEFer4o0EvAZH/wlnZu2rpIIYhGqxmgHhtT+3++VGSzkpldAnwNCAPfds795ZLbLXP7ZcA4cJ1z7pc5rhWA1OQIs0PZThN3wX/Tbv66w4FzwS3OLVgLMGN+ceYyZJZlrptzuMzPsdz9zt/ZgvudX+AW32aGc3bg987VsHTZ/M8u+N2Z6ziH9b5M2hmxVrX9S2EZDMWonOiFdJq5f6+BJZcPknn9zb8OD7VsBUt/xyGvL6nDLa3pULctV/tCS2qdqz0UgYrKFX72yK0Y6GYWBu4CfgfoAHab2SPOuZcXrHYpsDnzdT7wzcz3nHv2/z3AOU/fvBp3XbRagG6L01ZV5bsUkUWGIwnOGvtX+Itm36UUlOdOuI4zr/9azu83my3084C9zrk3AMxsJ7AdWBjo24G/d8Hm41Nm1mRma51zOd95Fj/5Ah4d/nz2P2BL/6of+G524O/r/NLMX99gdbdgmWPR1vrcFjuGwxbd36K/yUu2JubLcQ5j8V/64F4yW91zv985zEKZW4Pf6MwwO/B7AZrb34VGW0ihmdn2WX70zMPBu9H510pg6b/pOeaCd7yLXzkHtqiz2T63zD3Mv+vlwIa1W/g6dku3sZerjyX3syALlvzMYkvv+cD1xIb3ZvEojlw2gb4eWDhHqoODt76XW2c9kPNAb990Ku2bTs313YrIKtj6novgPRf5LqNsZPOh6HJ/fpbuOMpmHczsBjPbY2Z7ent1FjYRkVzKJtA7gIWftm0Alo4hyWYdnHP3OOe2Oue2trS0HGmtIiJyGNkE+m5gs5ltNLNK4BrgkSXrPAL8sQUuAIZWY/+5iIgc2or70J1zs2b2CeBHBIct3uuce8nMbszcfjewi+CQxb0Ehy1ev3oli4jIcrI6Dt05t4sgtBcuu3vBZQfclNvSRETkSKhTVESkRCjQRURKhAJdRKREmDvonAV5+sVmvcCbR/njCaAUps+WwuPQYygMegyFIR+P4QTn3LLHfXsL9GNhZnucc1t913GsSuFx6DEUBj2GwuD7MWiXi4hIiVCgi4iUiGIN9Ht8F5AjpfA49BgKgx5DYfD6GIpyH7qIiBysWLfQRURkCQW6iEiJKLpAN7NLzOw1M9trZrf6rudomNk+M3vBzJ41sz2+68mGmd1rZj1m9uKCZTEz+7GZ/TrzvaDnjB3iMdxuZu9knotnzewynzWuxMyOM7OfmNkrZvaSmf1pZnnRPBeHeQxF81yYWdTMnjaz5zKP4c8zy70+D0W1Dz0z3/RXLJhvCuxYMt+04JnZPmCrc65omijM7P3AKMGowdMzy+4A+p1zf5n549rsnPuMzzoP5xCP4XZg1Dn3JZ+1ZcvM1gJrnXO/NLN64Bng94HrKJLn4jCP4Q8pkufCzAyodc6NmlkE+Dnwp8CVeHweim0LfX6+qXNuGpibbyqrzDn3BNC/ZPF24DuZy98heFEWrEM8hqLinOt0zv0yc3kEeIVg3GPRPBeHeQxFwwVGM1cjmS+H5+eh2AL9ULNLi40D/q+ZPWNmN/gu5hi0zQ0yyXxv9VzP0fqEmT2f2SVTsLsqljKzduBs4N8o0udiyWOAInouzCxsZs8CPcCPnXPen4diC/SsZpcWgf/gnHs3cClwU2ZXgPjxTWATcBbBUPMve60mS2ZWB3wfuNk5N+y7nqOxzGMoqufCOZdyzp1FMHLzPDM73XNJRRfoWc0uLXTOuf2Z7z3AQwS7kopRd2Z/6Nx+0R7P9Rwx51x35oWZBr5FETwXmX223we+65x7MLO4qJ6L5R5DMT4XAM65QeBx4BI8Pw/FFujZzDctaGZWm/kgCDOrBS4GXjz8TxWsR4CPZS5/DHjYYy1HZe7Fl3EFBf5cZD6M+1/AK865ryy4qWiei0M9hmJ6LsysxcyaMpergQ8Br+L5eSiqo1wAMocy3cmB+aZf8FvRkTGzEwm2yiEYAXhfMTwGM/sesI3g9KDdwG3APwL3A8cDbwFXOecK9kPHQzyGbQRv8R2wD/hPhTzg3Mx+G/gZ8AKQziz+M4J90EXxXBzmMeygSJ4LMzuD4EPPMMGG8f3Oub8wszgen4eiC3QREVlese1yERGRQ1Cgi4iUCAW6iEiJUKCLiJQIBbqISIlQoIuIlAgFuohIifj/DdejMHQs/+cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAE6CAYAAADgN6LXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABGKElEQVR4nO3dfXCV9Z3//9fJObkh3MpNMBAoskFuAiEKCHb6tVqHYjNtrLJr6borDtos1I5ub3T9x7Hd2ba43c7aKU5pWkuxW6E3s5puC6motWvtSkRFp2VbUwtbiBG5CwRyc07OOb8/+BGJJNfnk5zr5LquD8/HDDMkn/f5XJ/rOtd5nSufc53rimWz2awAAAAAAABCrCDoAQAAAAAAAJgwgQEAAAAAAEKPCQwAAAAAABB6TGAAAAAAAIDQYwIDAAAAAACEHhMYAAAAAAAg9HKawGhqatLcuXNVWVmpjRs3+jUmAIAB+QsAwSB/ASA4sWw2mx3OA9PptC6//HLt2rVLFRUVWrZsmbZt26YFCxb4PUYAwHnIXwAIBvkLAMEa9hkYzc3Nqqys1OzZs1VUVKQ1a9aosbHRz7EBAAZA/gJAMMhfAAjWsCcwWltbNWPGjL6fKyoq1Nra6sugAACDI38BIBjkLwAEKzHcBw70zZNYLHbB7xoaGtTQ0CBJevnll4e7OAAY1KRJk3T06NELfn/DDTcM+HuTyZMnq6mpyY+h5QX5i8EMtB9E2TC/5YoRRP7a5W8i4X3I7ce+nk6nc+5DkuLxeE7tkozra9NPUVGRsY/CwkJjjU0/p06dyrmPVCplrDl9+rRneyaTMfZhs/0LCsyfUff29uY8ltGjRxtrbPYF01hs9m0/avx6DdkwPUc27+fjx4+/qPL3nGFPYFRUVOjgwYN9Px86dEjTpk27oK6+vl719fWS3DuwshGmdR6psfixnJHcbqaDBpsAHykczA9s1qxZA/7+6NGjeumll4bc37Jly3IcUX4NN39tDmgQHD9yz7Xn2CZ/ycVgkb/m/C0sLNSkSZM8+00mk57tNvt5R0eHscbmNTVu3DjP9vHjxxv7mDhxorHGtE0qKiqMfVx66aXGmpkzZxprfvnLX+bcx+HDh401zz//vGd7Z2ensY8JEyYYa8aMGWOsMf2Ba7M/LV++3FhjM9729nbPdtMEkySdPHnSWHP8+HHPdtMEk2Q3yWHzfj5q1CjP9uLiYmMfg71GXM3fc4Z9pLNs2TK1tLRo//79SiaT2r59u+rq6vwcGwDkLJvNDvlf2JG/AKKA/AWAYLiYv+cM+wyMRCKhTZs2adWqVUqn01q3bp2qqqr8HBsA5CxKgWyL/AUQBeQvAATDxfw9Z9gTGJJUW1ur2tpav8YCAL6K2ozyUJC/AMKM/AWAYLicv1KOExgAEHYuBzgAhBn5CwDBcDl/3braFwAAAAAAcBJnYABwmssz0AAQZuQvAATD5fxlAgOA01wOcAAIM/IXAILhcv4ygYFA2Nwf+WLkctgEhW0KAMEgf8+KxWKKx+OeNaZ2G4mE+bC+p6fHWJPJZDzb0+m0sQ+bGtNyUqmUL8tJJpM597Nv3z5jH5dddpmx5u233/Zsr6ioMPZRWlrqS824ceM827u7u419tLe3G2tKSkqMNaZl+bXPmf7+sMksm/3JZp392P+9uJy/TGAAcJbrV2EGgLAifwEgGK7nLxMYAJzmcoADQJiRvwAQDJfzlwkMAE5zOcABIMzIXwAIhsv5ywQGAKe5HOAAEGbkLwAEw+X8ZQIDgNNcDnAACDPyFwCC4XL+MoEBwFmuX8QIAMKK/AWAYLiev0xgAHCaywEOAGFG/gJAMFzOXyYwLiI2O7Lp/sg2/OjDLzZjMdWY7tOMcHM5wHFxCVO22ojaeOE/8vddpm1RUFDg2Z5Op43LMPUh2b0uu7q6PNuLi4uNfXR3dxtrenp6PNuTyaSxj87OTmPNiRMnjDVHjhzxbP/f//1fYx/PPvussWbs2LGe7RMnTjT2MWPGDGONzXNk2r4TJkww9mGz/c+cOWOsMe0vvb29xj5sakzrbLPPFRYWGmtstr/ptZjre6jL+WtOOgAAAAAAgIBxBgYAp7k8Aw0AYUb+AkAwXM5fJjAAOMv1ixgBQFiRvwAQDNfzlwkMAE5zOcABIMzIXwAIhsv5ywQGAKe5HOAAEGbkLwAEw+X8ZQIDgNNcDnAACDPyFwCC4XL+MoEBwGkuBzgAhBn5CwDBcDl/mcAA4CzXL2IEAGFF/gJAMFzPXyYwMCRFRUW+9JPJZDzbY7GYsY9Ewrz72tQUFBR4tp86dcrYh2l9JLdnQsOM7Q5XFBcXG2v8yCJTJkpSPB431tjkrynrbfLXBjkQDLb7u0yvGdPrzuZ1OWrUKGONTUaYanp6eox92NR0d3d7tp85c8bYx8mTJ401paWlxpp9+/Z5ts+cOdPYx5/+9CdjjWnbHjhwwNiHzevKJqP37t3r2V5dXW3so7Oz01izf/9+Y43p/aKkpMTYh80+19XV5dnu1/tfb2+vsca0LJvXqheX85cJDABOcznAASDMyF8ACIbL+csEBgCnuRzgABBm5C8ABMPl/GUCA4DTXA5wAAgz8hcAguFy/pq/6AMAAAAAABAwzsAA4CzXr8IMAGFF/gJAMFzPXyYwADjN5QAHgDAjfwEgGC7nLxMYAJzmcoADQJiRvwAQDJfzlwkMAE5zOcABIMzIXwAIhsv5ywRGntnsPLFYbARGYqepqcmzffny5cY+CgpyvzZsOp021nR3dxtrjh8/bqx59dVXPdv/4R/+wdhHV1eXsSaTyRhrXA6boLBNERVPPfWUZ/uyZcuMffjxfmKTVSOVv+vXrzf20dnZaawhB4LBdn+X6bXpx7ayOf4qLi421pw5c8az3eaYx2Y5p06d8mwvKSkx9hGPx401vb29xpqWlhbP9r//+7839nH06FFjTVFRkWd7ImH+0yyVShlrbN4LrrvuOs/2N954w9jHqFGjjDU27xemHLd5X7LZL01s9ieb15nN82hi8zx7cTl/mcAA4CzXL2IEAGFF/gJAMFzPXyYwADjN5QAHgDAjfwEgGC7nLxMYAJzmcoADQJiRvwAQDJfzlwkMAE5zOcABIMzIXwAIhsv5m9MExqxZszR27FjF43ElEgnt2bPHr3EBgC9cDnAyGECYkb8AEAyX8zfnMzB+9atfafLkyX6MBQB85fpFjCQyGEA4kb8AEAzX85evkABwmssBDgBhRv4CQDBczl/zjWw9xGIxffjDH9aSJUvU0NDg15gAABbIYAAIBvkLAMHI6QyMF154QdOmTdM777yjlStXat68ebrmmmv61TQ0NBDsBqYZslgs5styJkyYYKy5+uqrPdvj8bixj+LiYtsh5WTMmDHGmksuucRYYzr9c+zYscY+UqmULzUmLs+m5ovL28yUweRveNhk0VVXXeXZXlBg/syhpKTEWOPHe4pN/tq850yaNCnn5dhkazKZNNa4nBVBcXmbDiV/M5mM8fjJtK1sjr8ymYyxxuYYrbe317O9s7PT2Ed3d7ex5syZM57t7e3txj5s2FyfZPXq1Z7tbW1txj5OnTplrLE5pjSx2f42+8KJEyc8223eK2zel6ZOnWqsaW1t9Ww3jVWyG6/pdVRYWJhzH7ZjSafTnu255qfL+ZvTGRjTpk2TJJWVlemmm25Sc3PzBTX19fXas2cPFzcCEIhz3wMcyr+oMGUw+QsgSOTv2fy1+YMHAPzkcv4OewLjzJkz6ujo6Pv/U089pYULF/o2MADwQz4CvKmpSXPnzlVlZaU2btx4QfvJkyf1sY99TIsXL1ZVVZW2bNni+3qRwQDCjvwFgGC4PIEx7K+QHD58WDfddJOks6ea/e3f/q1uuOEG3wYGALnKRyCn02nddddd2rVrlyoqKrRs2TLV1dVpwYIFfTWPPPKIFixYoP/6r//SkSNHNHfuXN16660qKirybRxkMIAwI38BIBhRm5AYqmFPYMyePVuvvfaan2MBAN/5HeDNzc2qrKzU7NmzJUlr1qxRY2NjvwPoWCymjo4OZbNZnT59WhMnTlQi4e9Nn8hgAGFH/gJAMJjAAICI8jvAW1tbNWPGjL6fKyoqtHv37n41n/nMZ1RXV6dp06apo6NDP/rRj6wudAUALiF/ASAYLk9gkOgAnDac7wAeOXJES5cu7ft3/p08BnpDeO/Vpn/5y1+qpqZGb731lvbu3avPfOYzVlcnBwCXkL8AEIwgrkGUzWZ19913q7KyUtXV1XrllVeMj927d69WrFihmpoaLV26dMCbgrwXZ2AAcNpwZqCnTJky6J07KioqdPDgwb6fDx061Hc1+nO2bNmi+++/X7FYTJWVlbrsssv0hz/8wXibTABwCfkLAMEI4hpEO3fuVEtLi1paWrR7925t2LBBu3fv9nzsfffdpwcffFAf+chHtGPHDt1333167rnnPMfCBIYjxowZY6yxuQ+z6TTLUaNGWY8pDGxuXTZlyhTP9okTJxr78ON+55L53usunw6WD/m4iNGyZcvU0tKi/fv3a/r06dq+fbsef/zxfjUzZ87UM888o//3//6fDh8+rD/+8Y9939mGe0Yqf0tLS63HlG8297i3uWjiSOVvJpMx1qTT6Zz7wLvI36Exvf5tXnM229t0nCFJxcXFOS/H5pjHtM42x3B/+MMfjDU24z19+rRn++jRo419mPJMMufISH7dqaSkxLPdlImS1NXVZayxeS+YMGGCZ3tnZ6exj8LCQmONaZ38uuWxzbZLpVKe7abnx0s+8tfmGkSNjY267bbbFIvFtGLFCrW3t6utrU0HDhwY9LGxWKzvLLmTJ09eMCk9ECYwADjN7wBPJBLatGmTVq1apXQ6rXXr1qmqqkqbN2+WJK1fv14PPPCAbr/9di1atEjZbFYPPfSQJk+e7Os4ACDsyF8ACEYQ1yAaqKa1tdXzsQ8//LBWrVqlL3zhC8pkMvrtb39rHAsTGACclo+zVmpra1VbW9vvd+vXr+/7/7Rp0/TUU0/5vlwAiBLyFwCCMZz8PXcNonPq6+tVX18/aH/vPTNssBqvx37rW9/Sv//7v2v16tX68Y9/rDvuuENPP/205ziZwAAAAAAA4CKW6zWIBqtJJpODPnbr1q36xje+IUn6m7/5G915553GcXIXEgBOy8dVmAEAZuQvAATD7/w9/xpEyWRS27dvV11dXb+auro6PfbYY8pms3rxxRc1fvx4lZeXez522rRp+vWvfy1JevbZZzVnzhzjunEGBgCncUAMAMEgfwEgGEFcg6i2tlY7duxQZWWlSktLtWXLFs/HStJ3vvMd3XPPPert7VVJSUm/W2cPOhZf1wwAQoRP9AAgGOQvAAQjX/lrugZRLBbTI488Yv1YSfrABz6gl19+eUjjYAIDgNM4gAaAYJC/ABAMl/OXCQwATnM5wAEgzMhfAAiGy/nLBEYEjB8/3lhz6aWXGmvmzZtnrCkuLrYak0tML/Dz71s8mHQ6baw5cuSIsaazs9Ozvaury9gH+nM5wJF/EyZMMNaUl5cbay6//HJjzcWYvyY2+dvb22usOXbsmLHmzJkznu3d3d3GPtAf+Wsvk8l4thcUmK+7X1RUZKzx4zmxec3ZjPf06dOe7TbHVqbtJkmTJk0y1piOr+LxuLEPG6blFBYWGvtIJMx/vtm8n1x99dWe7YcOHTL2kUqljDWm51mSxowZ49lus2/b7C+m7fve24IOxGafs+nH9Dwmk0ljH15czl8mMAA4zeUAB4AwI38BIBgu5y8TGACcxUXkACAY5C8ABMP1/GUCA4DTXA5wAAgz8hcAguFy/jKBAcBpLgc4AIQZ+QsAwXA5f5nAAOA0lwMcAMKM/AWAYLicv+ZLBAMAAAAAAASMMzAAOM3lGWgACDPyFwCC4XL+MoEBwFmuX4UZAMKK/AWAYLiev0xghMDy5cs922tqaox9jB071lgzevRoY008HjfWuCYWi3m2X3HFFcY+KisrjTXHjh0z1vzlL3/xbH/hhReMfaA/lwMcubvqqqs828nfYNlsf5v8PX78uLHm//7v/zzbf/vb3xr7QH/k71nZbFa9vb3GGi+ZTMaXsdjkjKnGdNwkybi+knTZZZd5tr/99tvGPk6fPm2sOXXqlLFmwoQJnu1dXV3GPmy2i6kmmUwa++js7DTWnDhxwljz5z//2bPdZp0LCsxXIxg/fryxxrQsm21rkzfpdDrn5fjFj9eZF5fzlwkMAE5zOcABIMzIXwAIhsv5ywQGAKe5HOAAEGbkLwAEw+X8ZQIDgNNcDnAACDPyFwCC4XL+MoEBwFmuX8QIAMKK/AWAYLiev0xgAHCaywEOAGFG/gJAMFzOXyYwADjN5QAHgDAjfwEgGC7nLxMYAJzmcoADQJiRvwAQDJfzlwmMPLO537bpPveXXHKJsY/CwkJjTWlpqbEGF7LZbn5t/+LiYs/2F154wdgH+nM5wOHN5v70V1xxhWf7hAkTjH2Qv/kzkvlbVFTk2f7b3/7W2Af6I3/t9fb2erbb7OemPiQpFosZa0zLGj16tLGPo0ePGmtOnjzp2d7Z2Wnsw4bNe4HpeL2rq8vYh83+bjrOs1lOImH+8y2VShlrTOtssxybde7o6DDWmLaLzXNo8xoxrZPN+vj1OjP1Y9OHF5fz17w3AAAAAAAABIwzMAA4y/WrMANAWJG/ABAM1/OXCQwATnM5wAEgzMhfAAiGy/nLBAYAp7kc4AAQZuQvAATD5fw1XgNj3bp1Kisr08KFC/t+d/z4ca1cuVJz5szRypUrdeLEibwOEgCG69xpdEP5FxbkL4Aoi3L+SmQwgOiKev56MU5g3H777Wpqaur3u40bN+r6669XS0uLrr/+em3cuDFvAwSAXEQ5wMlfAFEW5fyVyGAA0RX1/PVinMC45pprNHHixH6/a2xs1Nq1ayVJa9eu1ZNPPpmXwQFALoYT3mEKcPIXQFRFPX8lMhhANLmQv16GdQ2Mw4cPq7y8XJJUXl6ud955x9dBAYBfohTINshfAFHhWv5KZDCAaHAxf8/J+0U8Gxoa1NDQkO/FhJbNzjN27FjP9kTC/DQVFhYaa+LxuLEGF7LZbjbPc2lpqbHm3EER/ONygJtc7PlrY8yYMZ7t5G+wbLa/DZvtf+mll/qyLLyL/D2bv+l0WqlUyrM+Fot5tqfTaV/GlclkjDVdXV0595FMJo013d3dnu2dnZ3GPkzbTbLbD9vb2z3bi4qKjH309PQYa0zPo8362Gx/m7GMHj3as90mNwsKjCfzWz2Ppn3hkksuMfZheo1J5u1rs21t9ieb44J8czl/zXvdAKZOnaq2tjZJUltbm8rKygatra+v1549e7Rnz57hjRAAcuDaKXTkL4CocC1/JfsMPj9/bf7IAwA/uZi/5wwrUevq6rR161ZJ0tatW3XjjTf6OigA8ItrAU7+AogK1/JXIoMBRIOL+XuOcQLjk5/8pK6++mr98Y9/VEVFhR599FHdf//92rVrl+bMmaNdu3bp/vvvH4mxAsCQRP0iRuQvgKiKev5KZDCAaHIhf70Yv1y6bdu2AX//zDPP+D4YAMC7yF8ACA4ZDADhk/eLeAJAkKI0owwALiF/ASAYLucvExgAnOZygANAmJG/ABAMl/OXCQwATnM5wAEgzMhfAAiGy/nLBAYAp7kc4AAQZuQvAATD5fxlAiPPMpmMsSaR8H4a4vG4sY9YLGas4T7kw2Oz3fza/sXFxVZjgp2oXVUZ/hqp/LV5bZO/w2OTrTZstn9JSYkvy8JZ5O+7stmsent7PWtM+6jNPtzT02OsKSwsNNak02nPdpvn1WY5HR0dnu1FRUXGPkzbVbI7tkomk57tXV1dxj5s8so0Xpv3LZv3JZsa0zrb7HM2z5HNvjB+/HhjjYnNfml6jmz68Os931RjOj7x4nr+MoEBwGkuBzgAhBn5CwDBcDl/mcAA4DSXAxwAwoz8BYBguJy/TGAAcJrLAQ4AYUb+AkAwXM5fJjAAOM3lAAeAMCN/ASAYLucvExgAnOX6RYwAIKzIXwAIhuv5ywQGAKe5HOAAEGbkLwAEw+X85b5uAAAAiISmpibNnTtXlZWV2rhx44A1zz33nGpqalRVVaUPfvCDIzxCAEA+MYEBwGnnTqMbyj8TDqABwMzv/E2n07rrrru0c+dO7du3T9u2bdO+ffv61bS3t+vTn/60fvazn+n3v/+9fvKTn+RzFQEglPJx/BsWfIUkz2KxWM41BQXmeSabmijtmGFi8xz61U88HvdlWXiX3/v9uQPoXbt2qaKiQsuWLVNdXZ0WLFjQV3PuALqpqUkzZ87UO++84+sYYMfmNWfKTptstVkO+Ts8fuWvDZvnGkPj937f3NysyspKzZ49W5K0Zs0aNTY29svfxx9/XDfffLNmzpwpSSorK/N1DMNl2hbJZNKz3eb4wK99OJVKebaXlJQY+zhz5oyxZsyYMZ7tnZ2dxj5s1rm3t9dYk06nPdttssi03STzfmCznEwmk/NyJPM62yzHZr+0GYtp29k8z0VFRcYaP57nRMKfP59N65Trclw+7uDdGoCzhjP7bAr88w+gi4qK+g6gzxfWA2gAGCn5yN/W1lbNmDGj7+eKigq1trb2q3njjTd04sQJXXvttVqyZIkee+yxvKwfAIRVPvJXMp+BnM1mdffdd6uyslLV1dV65ZVXrB77zW9+U3PnzlVVVZXuu+8+4zg4AwOA0/yegR7oAHr37t39at544w2lUilde+216ujo0D333KPbbrvN13EAQNgNJ3+PHDmipUuX9v1cX1+v+vr6Qft77yemvb29evnll/XMM8+oq6tLV199tVasWKHLL798yGMBgKgK4gzknTt3qqWlRS0tLdq9e7c2bNig3bt3ez72V7/6lRobG/X666+ruLjY6qxlJjAAOI0DaAAIxnDyd8qUKdqzZ8+AbRUVFTp48GDfz4cOHdK0adMuqJk8ebJGjx6t0aNH65prrtFrr71G/gK4qATxFb7GxkbddtttisViWrFihdrb29XW1qYDBw4M+thvfetbuv/++1VcXCzJ7qxlvkICwGnDOYXu3AH0uX/nJi8k+wPoG264QaNHj9bkyZP7DqAB4GLi9ynMy5YtU0tLi/bv369kMqnt27errq6uX82NN96o559/Xr29vers7NTu3bs1f/78fK4mAIROEF/hG6zG67FvvPGGnn/+eS1fvlwf/OAH9dJLLxnXjQkMAE7jABoAguF3/iYSCW3atEmrVq3S/Pnzdcstt6iqqkqbN2/W5s2bJUnz58/XDTfcoOrqal111VW68847tXDhwpFYXQAIjeHk77kzkM/9a2ho6Nffe733DOTBarwe29vbqxMnTujFF1/U1772Nd1yyy3m9wKrLQAAEWR7UaKhOP8AOp1Oa926dX0H0JK0fv36fgfQBQUFHEADuOjkI38lqba2VrW1tf1+t379+n4/33vvvbr33nt9XzYARMFw89ePr/ANVJNMJgd9bEVFhW6++WbFYjFdddVVKigo0NGjRzVlypRBx8kEBgCncQANAMHIR/4CAMz8zt/zz0CePn26tm/frscff7xfTV1dnTZt2qQ1a9Zo9+7dGj9+vMrLyzVlypRBH/vxj39czz77rK699lq98cYbSiaTmjx5sudYmMDIM5udx3QfYJt7Evt132jANRxAX7xsnntTdvqVv6Z7z2NgNtsW4UX+nhWLxYzHeiapVMpYU1hY6Es/Jl1dXTn3IUmnT5/2bE8kzH+m+JWtyWTSs72oqGhExjKSmWd6fdrsTz09PcYam35M+9To0aONfdj8nROPxz3b/dr+Ntln2r9zzYwgzkCura3Vjh07VFlZqdLSUm3ZssXzsZK0bt06rVu3TgsXLlRRUZG2bt1qfB6YwADgNA6gASAY5C8ABCOIM5BjsZgeeeQR68dKZycH/+M//mNI42ACA4DTOIAGgGCQvwAQDJfzl7uQAAAAAACA0OMMDADOytdV8AEA3shfAAiG6/nLBAYAp7kc4AAQZuQvAATD5fxlAgOA01wOcAAIM/IXAILhcv4ygQHAaS4HOACEGfkLAMFwOX+ZwADgNJcDHADCjPwFgGC4nL9MYIRALBYbkeUUFHDTmeEYqecH/nP9IkbIHfkbbuRvdJG/78pms+rt7fWsSafTnu2JhPmQPZPJGGtM45Ck4uJiz/bOzk5jH368dm3Wx2Y5PT09OfdjM5bCwkJjTSqVymkckt0fpjbvOaZ1sllnm7GY9m3JvF/abBebGtM62bzO/HpfMj1HuSzH9fxlAgOA01wOcAAIM/IXAILhcv4ygQHAaS4HOACEGfkLAMFwOX+ZwADgNJcDHADCjPwFgGC4nL9MYABwmssBDgBhRv4CQDBczl/jFV7WrVunsrIyLVy4sO93X/ziFzV9+nTV1NSopqZGO3bsyOsgAWA4zl3EaKj/woL8BRBVUc9fiQwGEE0u5K8X4wTG7bffrqampgt+/9nPflZ79+7V3r17VVtbm5fBAcDFjPwFgOCQwQAQPsYJjGuuuUYTJ04cibEAgO+iPANN/gKIsijnr0QGA4iuqOevl2HfmH7Tpk2qrq7WunXrdOLECT/HBAC+cTHAyV8AUeBi/kpkMIDwczV/pWFexHPDhg164IEHFIvF9MADD+jzn/+8vve97w1Y29DQoIaGhpwGebGLxWK+1GQyGT+GA0RKlALZBvkbTa7th4ANF/d72ww+P38zmcyIHIOl02ljjc04ent7PduTyaSxj3HjxhlrTMeuPT09xj5s1icejxtrCgq8P9O16cNmvLmOQ/LvmD+R8P4z0LQf2I4llUoZa8aPH2+sMbHZ/03rbLP9bXLNtBzJvO1yzU8X8/ecYZ2BMXXqVMXjcRUUFOhTn/qUmpubB62tr6/Xnj17tGfPnmEPEgCGy7UZaPIXQFS4lr+SfQafn782f+QBgJ9czN9zhjWB0dbW1vf/J554ot/VmQEgLIYT3mEPcPIXQBS4mL8SGQwg/FzN33OM57d88pOf1HPPPaejR4+qoqJCX/rSl/Tcc89p7969isVimjVrlr797W+PxFgBYMiiFMjvRf4CiLIo569EBgOIrqjnrxfjBMa2bdsu+N0dd9yRl8EAgN+iHODkL4Aoi3L+SmQwgOiKev56GdZFPAEgKlwOcAAIM/IXAILhcv4ygQHAaS4HOACEGfkLAMFwOX+ZwADgrKhdlAgAXEH+AkAwXM9fJjAAA25/Fm0uBzjCgYwABkb+vivXbZFOp32pscmr3t5ez/Z4PG7so7u721iTyWRyapekwsJCY42NcePGebbbrE8qlTLWJBLef3qZtr1kt/0LCsw3mjTtkzb7rM3+ZDOWU6dOebabnh/Jbl8wrZPN9jc9hzbLsamx2f9zHUNUDes2qgAAAAAAACOJMzAAOM3lGWgACDPyFwCC4XL+MoEBwGkuBzgAhBn5CwDBcDl/mcAA4DSXAxwAwoz8BYBguJy/TGAAcJbrV2EGgLAifwEgGK7nLxMYAJzmcoADQJiRvwAQDJfzlwkMAE5zOcABIMzIXwAIhsv5ywQGAKe5HOAAEGbkLwAEw+X8ZQIjBNLp9IgsJxaLjchygDBxOcCRu0wmMyLLIX+Hh+0WbeTvWbFYTPF43LPGdCxoc6xo83rp7e011piWZbOcVCplrPFj/7BZH9O2l8zvBTbrY8O0zomE+U8zmxq/9hcTm22bTCZzHotfz7Np+9v04ddYTM+RTR9eXM5fJjAAOMv1ixgBQFiRvwAQDNfzlwkMAE5zOcABIMzIXwAIhsv5ywQGAKe5HOAAEGbkLwAEw+X8LQh6AAAAAAAAACacgQHAaS7PQANAmJG/ABAMl/OXCQwATnM5wAEgzMhfAAiGy/nLBAYAZ7l+FWYACCvyFwCC4Xr+MoEBwGkuBzgAhBn5CwDBcDl/mcAIAdMO5tcOmMlkfOnnYuNyAFwMeP7gZaRykfwdHl6/0cbzd1Y2m1UymfSsSSS8D8ltMiQWixlrent7jTUFBd7X+E+lUsY+ioqKjDUmNsuxWeeSkpKcl+XX9vejD7/eT+LxuGe7zevXpsZmXzCt06hRo4x92OzbJul02pflmLatzbJyXR+X85cJDABOcznAASDMyF8ACIbL+csEBgCnuRzgABBm5C8ABMPl/PU+RwwAIuzcRYyG+g8AkJt85W9TU5Pmzp2ryspKbdy4cdC6l156SfF4XD/96U/9XC0ACD3Xj385AwOA06IUyADgEr/zN51O66677tKuXbtUUVGhZcuWqa6uTgsWLLig7p/+6Z+0atUqX5cPAFHh8vEvZ2AAcBqfAAJAMPzO3+bmZlVWVmr27NkqKirSmjVr1NjYeEHdN7/5Ta1evVplZWX5WjUACDWXz8BgAgOA0/wO8HOfAO7cuVP79u3Ttm3btG/fvgHr+AQQwMXM7/xtbW3VjBkz+n6uqKhQa2vrBTVPPPGE1q9fn5d1AoAoCOIDvGw2q7vvvluVlZWqrq7WK6+8Yv3Yf/u3f1MsFtPRo0eN42ACAwCGgE8AASB/jhw5oqVLl/b9a2ho6Gsb6AD7vbed/Md//Ec99NBDVrcxBADYsfkAb+fOnWppaVFLS4saGhq0YcMGq8cePHhQu3bt0syZM63GwjUwADjN71PiBvoEcPfu3RfUPPHEE3r22Wf10ksv+bp8AIiK4eTvlClTtGfPngHbKioqdPDgwb6fDx06pGnTpvWr2bNnj9asWSNJOnr0qHbs2KFEIqGPf/zjQx4LAESV38e/53+AJ6nvA7zzr0HU2Nio2267TbFYTCtWrFB7e7va2tp04MABz8d+9rOf1b/+67/qxhtvtBoLExh5ZvMJQDqdznk57/0EYiCJBE83Li7D/U7fuU8Az6mvr1d9fX1fn+/FJ4Dh5Ef++nUAQP7iYpOP71QvW7ZMLS0t2r9/v6ZPn67t27fr8ccf71ezf//+vv/ffvvt+uhHPxr45EU2mzVmjek4LpVKGZdjkzM2x4t+HJfa9DFq1CjPdpt1ttnHkslkzjUlJSXGPmzGa3pf6u3tNfZRVFRkrLHZLplMJuex2OxzNu/FpmWdPn3a2IfNdvHjdWbDZv/v7u72bLd5rQ4mH8e/th/gDfQ1P6/H/uxnP9P06dO1ePFi63FyRAXAaXwCCADB8HsCI5FIaNOmTVq1apXS6bTWrVunqqoqbd68WZK47gUA/P/8Pv61+QBvsJrBft/Z2akvf/nLeuqpp4Y0TiYwADiNTwABIBj5uKp9bW2tamtr+/1usImL73//+74vHwCiwO/8tfkAb7CaZDI54O/ffPNN7d+/v+/si0OHDunKK69Uc3OzLr300kHHwkU8ATjN76swn/8J4Pz583XLLbf0fQJ47lNAAIDbt/EDgDDzO3/P/wAvmUxq+/btqqur61dTV1enxx57TNlsVi+++KLGjx+v8vLyQR+7aNEivfPOOzpw4IAOHDigiooKvfLKK56TF5LFGRgHDx7UbbfdprffflsFBQWqr6/XPffco+PHj+sTn/iEDhw4oFmzZunHP/6xLrnkEovNCQAjJ8qfAJK/AKIsyhMS5C+AKAviK3y1tbXasWOHKisrVVpaqi1btng+dthjsRns17/+dV155ZXq6OjQkiVLtHLlSn3/+9/X9ddfr/vvv18bN27Uxo0b9dBDDw17IADgt6h/okf+Aogq8hcAgpGv/DV9gBeLxfTII49YP/a9Dhw4YDUO41dIysvLdeWVV0qSxo4dq/nz56u1tVWNjY1au3atJGnt2rV68sknrRYIACMpyqcwk78Aooz8BYBgRDl/TYZ0Ec8DBw7o1Vdf1fLly3X48GGVl5dLOhvy77zzTl4GCAC5iFIgeyF/AUQN+QsAwXAlfwdiPYFx+vRprV69Wg8//LDGjRtnvYCGhgY1NDQMa3AuKCgwXyfVdM9hP+7HLdndhxkXGskAsLlXOYbGhQAnf4fHj/zNZDLGPmzu1U7+Ds9Ivn5N+wKGjvx9N39NeWRqt8kQm+1tU2OTaSZFRUXGGtM6JxLmP1NsXre9vb3GGhOb4zOb7Wba/jbPj83fBTb7i6nG5j3U5nnu7OzMeSw2z6HNWLq7uz3bbfY5m+3S1dVlrPHr77vBuJC/g7G6C0kqldLq1at166236uabb5YkTZ06VW1tbZKktrY2lZWVDfjY+vp67dmzZ9B7ygJAPkX9FDryF0BUkb/kL4BgRD1/vRgnMLLZrO644w7Nnz9fn/vc5/p+X1dXp61bt0qStm7dqhtvvDF/owSAYRhOeIcpwMlfAFFF/gJAMKKevybG82ReeOEF/eAHP9CiRYtUU1MjSfrKV76i+++/X7fccoseffRRzZw5Uz/5yU/yPVYAuKiQvwAQDPIXAMLJOIHxgQ98YNAZmWeeecb3AQGAn6I0o/xe5C+AKCN/ASAYUc5fkyHdhQQAosblAAeAMCN/ASAYLucvExgAnOZygANAmJG/ABAMl/OXCQwATnM5wAEgzMhfAAiGy/nLBAYAZ0XtqsoA4AryFwCC4Xr+MoGRZwUFxjvVqqury7N99OjRxj5sdlKbsZj6icVixj5c41cAZDIZY01nZ6cvy8K7XA5weBup/LV5bcfjcWMN+Xshv16/Nv2Qv/4jf9/V29vr2Z5Opz3bbbalXzUmNnnW3d1trDHlq81YU6mUsWbUqFHGGtP2t1mOqQ/J/L5ks21t9PT0GGsSCe8/A23e2/za5yZPnuzZfuzYMWMfyWTSWGNaJ5t1ttkXbBQWFnq2nzlzJqf+Xc5fJjAAOM3lAAeAMCN/ASAYLucvExgAnOZygANAmJG/ABAMl/OXCQwATnM5wAEgzMhfAAiGy/nLBAYAZ7l+ESMACCvyFwCC4Xr+MoEBwGkuBzgAhBn5CwDBcDl/zZdoBwAAAAAACBhnYABwmssz0AAQZuQvAATD5fxlAgOA01wOcAAIM/IXAILhcv4ygZFnsVjMWHPixAnP9tLSUmMfRUVFxppEwvx0m3Z2m/VxjU0A2NR0dXUZa9566y2rMcGO6xcxgjebvGpvb/dsHzVqlLEPm/yNx+PGGvL3Qn69fm3y9+233/ZlWTiL/H2XzbYwvb79OhaxkclkPNsLCszfQLfJqzFjxni222Rrd3e3LzWm8dqsczqdNtb09vbmNA5bNv0kk8mcl1NSUmKssdl206ZN82z3K59N78U9PT3GPmxeZzbrbOrH5rjBq2+X85cJDABOcznAASDMyF8ACIbL+csEBgCnuRzgABBm5C8ABMPl/GUCA4DTXA5wAAgz8hcAguFy/jKBAcBpLgc4AIQZ+QsAwXA5f5nAAOAs1y9iBABhRf4CQDBcz18mMAA4zeUAB4AwI38BIBgu5y8TGACc5nKAA0CYkb8AEAyX85cJDABOcznAASDMyF8ACIbL+csERp6l02ljzaFDhzzbCwoKjH2kUiljzahRo4w1pvHajMU1mUzGWHPq1CljTWtrq7HmzTfftBoTADM/8jcWixn76O3tNdbY5K8pa8jfgdnkb1tbm7HmT3/6k9WYgKGKxWLGLDG127wWbGpsJBLefx6MHz/e2Mfp06eNNfF43LN9woQJxj5sjn9PnjxprDGxyXkbpj8qbf7oND0/fjE9P5JUUlJirLF5/zNtX5s+urq6jDW5vg4lu+1i81o0PdcX43u+LSYwADjN5RloAAgz8hcAguFy/jKBAcBZrl+FGQDCivwFgGC4nr9MYABwmssBDgBhRv4CQDBczl8mMAA4zeUAB4AwI38BIBgu5y8TGACc5nKAA0CYkb8AEAyX85cJDABOcznAASDMyF8ACIbL+csEBgBnuX4RIwAIK/IXAILhev4ygZFnNveNfv311z3b33zzTWMfl1xyibFm1qxZxpq77rrLs72wsNDYh2v27NljrHnppZeMNWfOnDHWnDp1ympMsOdygMNbKpUy1rz22mue7X7l7/ve9z5jzac//WnP9kTi4nvLtsnWl19+2VjT2dlprDl58qTVmGCP/D0rm80qnU4ba3Jpl6RYLDakcQ2moKDAs90mW8vLy401b731Vs7LsVln07aXpHg8bqwxsclom7GYZDKZnPuQpL/+67/2bD98+LCxj7fffttYY/MeeezYMc/20tJSYx82TH+X2ewHfu2XptdZrq9nl/P34jsaAnBRcTnAASDMyF8ACIbL+csEBgCnuRzgABBm5C8ABMPl/GUCA4DTXA5wAAgz8hcAguFy/np/+QYAAAAAACAEjBMYBw8e1HXXXaf58+erqqpK3/jGNyRJX/ziFzV9+nTV1NSopqZGO3bsyPtgAWAozl2Feaj/woL8BRBV5C8ABCPq+Wti/ApJIpHQ17/+dV155ZXq6OjQkiVLtHLlSknSZz/7WX3hC1/I+yABYLiiFMjvRf4CiDLyFwCCEeX8NTFOYJSXl/fdCmns2LGaP3++Wltb8z4wAPBDlAOc/AUQZeQvAAQjyvlrMqRrYBw4cECvvvqqli9fLknatGmTqqurtW7dOp04cSIvAwSAXLhyCh35CyBqyF8ACIYr+TsQ67uQnD59WqtXr9bDDz+scePGacOGDXrggQcUi8X0wAMP6POf/7y+973vXfC4hoYGNTQ0+DroKLHZGbq6ujzbk8lkzn34VVNaWmrsI2oymYxn+8svv2zso6Ojw1hz6tQp6zHBP1EK5MGQv/ljyrxUKpVzH5LU2dmZcz8XY/6++uqrxj7I3/DKR/42NTXpnnvuUTqd1p133qn777+/X/sPf/hDPfTQQ5KkMWPG6Fvf+pYWL1487OX5lb+xWMxzOabXgunxNn1IUlFRkbGmsLDQsz0ejxv7SKfTxpqpU6d6tp88edLYx+nTp40106ZNM9YcOXLEs91mfWy2S29vr2e7adtLdq+rSZMmGWuam5s92ydMmGDso6qqylhjk+Om/XLcuHHGPhIJ85+1PT09nu027+c2f5fZvF5NNTZ9eHHh+HcwVmdgpFIprV69WrfeeqtuvvlmSWdDJx6Pq6CgQJ/61KcGfRHU19drz5492rNnj3+jBgAL+bqIUVNTk+bOnavKykpt3LjxgvYf/vCHqq6uVnV1td7//vfrtddeG/Y6kL8Aoigf+ZtOp3XXXXdp586d2rdvn7Zt26Z9+/b1q7nsssv061//Wq+//roeeOAB1dfXD3sdyF8AURTU8W82m9Xdd9+tyspKVVdX65VXXjE+9t5779W8efNUXV2tm266Se3t7cZxGCcwstms7rjjDs2fP1+f+9zn+n7f1tbW9/8nnnhCCxcuNC4MAEZalA+gyV8AUeZ3/jY3N6uyslKzZ89WUVGR1qxZo8bGxn4173//+3XJJZdIklasWKFDhw4Ne+zkL4CoCuL4d+fOnWppaVFLS4saGhq0YcMG42NXrlyp3/3ud3r99dd1+eWX66tf/apx3Yzn2rzwwgv6wQ9+oEWLFqmmpkaS9JWvfEXbtm3T3r17FYvFNGvWLH372982LgwARprfp9CdfwAtqe8AesGCBX0173//+/v+n8sBNPkLIMr8zt/W1lbNmDGj7+eKigrt3r170PpHH31UH/nIR4a1LPIXQJQFcfzb2Nio2267TbFYTCtWrFB7e7va2tp04MCBQR/74Q9/uO/xK1as0E9/+lPjWIwTGB/4wAcG3AC1tbXmNQWAgEX5AJr8BRBlw8nfI0eOaOnSpX0/19fX953FNlB/g31P/Fe/+pUeffRR/eY3vxnyGCTyF0C0BXH8O1BNa2ur9bHz9773PX3iE58wjsX6Ip4AEEXDCfCjR4+G4gAaAKJsOPk7ZcqUQa8bUVFRoYMHD/b9fOjQoQEv0Pj666/rzjvv1M6dO60uZggArgni+HewGpvHfvnLX1YikdCtt95qHCcTGACcNdzbQnEADQC5GW7+el15f9myZWppadH+/fs1ffp0bd++XY8//ni/mr/85S+6+eab9YMf/ECXX375kJcPAFEX1PHvYDXJZNLzsVu3btXPf/5zPfPMM1Z3X7G6CwkA4KzzD6CTyaS2b9+uurq6fjUcQAOA/xKJhDZt2qRVq1Zp/vz5uuWWW1RVVaXNmzdr8+bNkqR//ud/1rFjx/TpT39aNTU1/T5NBAAMj83xb11dnR577DFls1m9+OKLGj9+vMrLyz0f29TUpIceekg/+9nPrG8XzxkYAJzm93cAzz+ATqfTWrduXd8BtCStX7++3wH0ucdwKz0AFxu/81c6ew2K916HYv369X3//+53v6vvfve7vi8XAKIkiOPf2tpa7dixQ5WVlSotLdWWLVs8HytJn/nMZ9TT06OVK1dKOnshz3N9DiaWzce7y2ALszgl5GJk2i422y0ejxtrSkpKjDWmP7IqKyuNfRQUhOfEHpvd++233/ZsX7RokbGP48ePW48p30bwJR0aS5YsGXDfLS0tHdYZEC5OOMRisVC9NqPCJn8TCfNnAaNGjTLWNDc3e7b/1V/9lbGPMD3HNll0/i0pB3LFFVcY+zh27Jj1mLxcjNnphyuuuIL8NfDjOC6dThv7sHn9FxcX59yPzaekRUVFxpqxY8d6tttk65EjR4w1c+fONda8+uqrnu1nzpwx9pHJZIw1JjbPj80xfyqVMtaYjm9tnkPTdpNk9TXaiRMnerbbbBeb56i7u9uz3WZ/SiaTxhqb7W96Hm1ez9XV1Rdl/nIGBgCn8UcJAASD/AWAYLicv0xgAHCaywEOAGFG/gJAMFzOXyYwADhruFdhBgDkhvwFgGC4nr9MYABwmssBDgBhRv4CQDBczl8mMAA4zeUAB4AwI38BIBgu5y8TGACc5nKAA0CYkb8AEAyX85cJDABOcznAASDMyF8ACIbL+csERgTY7IA29we3uSfx3r17PdsrKiqMfdjcH9wPNtvl8OHDxpqf//znnu3t7e22Q0LIuH4RI+SfX/lrc994U/5Onz7d2EfU8nfHjh2e7SdPnrQekxdyYOSRv0OTyWQ82+PxuLGPRMJ8WG9TU1xc7NleUlKScx82/ZSXlxv7mDJlirHmL3/5i7GmrKzMs/3EiRPGPjo6Oow1pmNx034gSePGjTPWzJgxw1jzpz/9ybP9Rz/6kbGPtWvX+jKWwsJCz3abfc7mvdi0fW3eQwsKCow1sVjMWGMar00fg3E9f5nAAOA0lwMcAMKM/AWAYLicv+YpJAAAAAAAgIBxBgYAp7k8Aw0AYUb+AkAwXM5fJjAAOM3lAAeAMCN/ASAYLucvExgAnOZygANAmJG/ABAMl/OXCQwAznL9KswAEFbkLwAEw/X8ZQIDgNNcDnAACDPyFwCC4XL+MoEBwGkuBzgAhBn5CwDBcDl/mcDIs1gsNiLLsdlJe3p6jDWf/OQn/RgOAmSzz7kcau91Ma0rgpHJZIw13d3dxppbb73Vj+FcdHiNhxfPzbtM783xeNyzvaCgwLiMwsJCY41pOZJUXFyc83LGjBmTc01JSYmxj1mzZhlrampqjDWmffXkyZPGPmzeC0zPo00fpufHtmb58uWe7b/5zW+MfSxZssRYc+bMGWPNpEmTPNs7OjqMfdisc29vr2f7qFGjjH3YPEc2Nclk0liTC5fzlwkMAE5zOcABIMzIXwAIhsv5ywQGAGe5fhEjAAgr8hcAguF6/jKBAcBpLgc4AIQZ+QsAwXA5f5nAAOA0lwMcAMKM/AWAYLicv0xgAHCaywEOAGFG/gJAMFzOX/MljQEAAAAAAALGGRgAnObyDDQAhBn5CwDBcDl/mcAA4CzXr8IMAGFF/gJAMFzPXyYwgItQLBbzbHcp9FxaF+Biw+s32nj+3pXr+248HjcuI5EwH9aXlJTkXFNYWGjso6ioyFgzevRoz/bS0tKc+5CkUaNGGWtMz4+pXZJOnz5trDHJZDLGGpt9wY9t953vfMfYxzXXXGOsmTx5srEmmUx6ttvsc37slzbbzeY5ssk+0zqn02ljH7mOIaqYwADgNJcDHADCjPwFgGC4nL9MYABwmssBDgBhRv4CQDBczl8mMAA4zeUAB4AwI38BIBgu56/xNqrd3d266qqrtHjxYlVVVenBBx+UJB0/flwrV67UnDlztHLlSp04cSLvgwWAoTh3EaOh/gsL8hdAVJG/ABCMqOeviXECo7i4WM8++6xee+017d27V01NTXrxxRe1ceNGXX/99WppadH111+vjRs3jsR4AWBIohzg5C+AKCN/ASAYUc5fE+MERiwW05gxYyRJqVRKqVRKsVhMjY2NWrt2rSRp7dq1evLJJ/M6UAAYjigHOPkLIMrIXwAIRpTz18Q4gSGdvY1LTU2NysrKtHLlSi1fvlyHDx9WeXm5JKm8vFzvvPNOXgcKAMMR9QAnfwFEFfkLAMGIev56sZrAiMfj2rt3rw4dOqTm5mb97ne/s15AQ0ODli5dqqVLlw57kAAwXFEPcPIXQFSRv+QvgGBEPX+9DOkuJBMmTNC1116rpqYmTZ06VW1tbSovL1dbW5vKysoGfEx9fb3q6+slnT0dD8Fh++ePzbaNUjAgfMhfAAhGrvlbUFCgRCK3G//F4/ERqzG9XxQXFxv7KCoqMtaUlpbm1C5JJSUlxpqJEycaa0x6enqMNTbPcSqVynksNseTmUzGWPPnP//Zs/3YsWPGPvzYtpJ06tQpz/auri5jHzbbv7Cw0LO9oMD82b5NjWk5knnfTSaTxj4uVsZn4MiRI2pvb5d0dud5+umnNW/ePNXV1Wnr1q2SpK1bt+rGG2/M60ABYKiGM/scpokm8hdAVJG/ABCMqOeviXGqqq2tTWvXrlU6nVYmk9Ett9yij370o7r66qt1yy236NFHH9XMmTP1k5/8ZCTGCwBDEqVAfi/yF0CUkb8AEIwo56+JcQKjurpar7766gW/nzRpkp555pm8DAoA/BLlACd/AUQZ+QsAwYhy/prk9oU8AAg5lwMcAMKM/AWAYLicv0xgAHCaywEOAGFG/gJAMFzOXyYwADgrahclAgBXkL8AEAzX85cJDABOcznAASDMyF8ACIbL+csERg5M98kGEDyXAxxwHa/fYGUymZwez/N3VjabNW7LeDzu2W7zXNhs74KCAmONSWFhobEmkTD/iWFaZ5v1GTVqlLFm/vz5xpqWlhbP9uLiYmMfNtvW9LeDzTrbbNtkMmmsOXz4sGd7WVmZsY/S0lJjTW9vb841Nn9zmfYnybzvFhUVGfuwGYtNjek56u7uNvaxdOnSQdtczl8mMAA4zeUAB4AwI38BIBgu5y8TGACc5nKAA0CYkb8AEAyX85cJDADOcv0iRgAQVuQvAATD9fzN/YtwAAAAAAAAecYEBgCnnZuFHso/AEDu8pG/TU1Nmjt3riorK7Vx48YBl3n33XersrJS1dXVeuWVV/KxagAQai4f//IVEgBOi1IgA4BL/M7fdDqtu+66S7t27VJFRYWWLVumuro6LViwoK9m586damlpUUtLi3bv3q0NGzZo9+7dvo4DAMLO5eNfzsAA4DSXZ6ABIMz8zt/m5mZVVlZq9uzZKioq0po1a9TY2NivprGxUbfddptisZhWrFih9vZ2tbW15XM1ASB0XD7+ZQIDgNNcDnAACDO/87e1tVUzZszo+7miokKtra1DrgEA17l8/DuiXyEZPXq05s2bN5KLHLYjR45oypQpQQ/DGuPNnyiNVYrWeP0a64EDBwb8/apVq3T06NEh9zd58uQcRxQ+Ucpf6eLcj0dKlMYbpbFK0RpvWPO3q6tLS5cu7fu5vr5e9fX1kgY+JToWi/X72aZmpLmWvzZ/6HR0dPhSY/LUU08Z9+P//M//zHk5folSRmQyGeNYn3766REajZkf2zYejxtr/DpGLCoq8hzv+Tk4mIv1+HdEJzDmzZunPXv2jOQih23p0qWRGavEePMpSmOVojXefI+1qakpb31HTZTyV2I/zqcojTdKY5WiNd4o5m9FRYUOHjzY9/OhQ4c0bdq0IdeMNPI3f6I0Vila443SWCXGez7Xj3/5CgkAAABCb9myZWppadH+/fuVTCa1fft21dXV9aupq6vTY489pmw2qxdffFHjx49XeXl5QCMGAPiNu5AAAAAg9BKJhDZt2qRVq1YpnU5r3bp1qqqq0ubNmyVJ69evV21trXbs2KHKykqVlpZqy5YtAY8aAOCnEZ3AOPcdxiiI0lglxptPURqrFK3xRmmsURe1bR2l8UZprFK0xhulsUrRGm+Uxnq+2tpa1dbW9vvd+vXr+/4fi8X0yCOPjPSwPEVtW0dpvFEaqxSt8UZprBLjvZjEslG65CgAAAAAALgocQ0MAAAAAAAQeiMygdHU1KS5c+eqsrJSGzduHIlF5mTWrFlatGiRampqrG5hM9LWrVunsrIyLVy4sO93x48f18qVKzVnzhytXLlSJ06cCHCE7xporF/84hc1ffp01dTUqKamRjt27AhwhP0dPHhQ1113nebPn6+qqip94xvfkBTO7TvYWMO6fbu7u3XVVVdp8eLFqqqq0oMPPigpnNvWNVHKYPLXP+Rv/pC/sEX++idK+StFK4PJ3/whf/Mgm2e9vb3Z2bNnZ998881sT09Ptrq6Ovv73/8+34vNyfve977skSNHgh7GoH79619nX3755WxVVVXf7+69997sV7/61Ww2m81+9atfzd53331BDa+fgcb64IMPZr/2ta8FOKrBvfXWW9mXX345m81ms6dOncrOmTMn+/vf/z6U23ewsYZ1+2YymWxHR0c2m81mk8lk9qqrrsr+z//8Tyi3rUuilsHkr3/I3/whf2GD/PVXlPI3m41WBpO/+UP++i/vZ2A0NzersrJSs2fPVlFRkdasWaPGxsZ8L9Zp11xzjSZOnNjvd42NjVq7dq0kae3atXryyScDGNmFBhprmJWXl+vKK6+UJI0dO1bz589Xa2trKLfvYGMNq1gspjFjxkiSUqmUUqmUYrFYKLetS8hgf5G/+UP+5g/5Gwzy119Ryl8pWhlM/uYP+eu/vE9gtLa2asaMGX0/V1RUhHonk87uaB/+8Ie1ZMkSNTQ0BD0cK4cPH+67z3l5ebneeeedgEfkbdOmTaqurta6detCe8rUgQMH9Oqrr2r58uWh377nj1UK7/ZNp9OqqalRWVmZVq5cGYltG3VRy2DyN//Cmg/nI3/9R/6OPPI3/6K4D4c1I84hf/1H/vor7xMY2QFuchKLxfK92Jy88MILeuWVV7Rz50498sgj+u///u+gh+SUDRs26M0339TevXtVXl6uz3/+80EP6QKnT5/W6tWr9fDDD2vcuHFBD8fTe8ca5u0bj8e1d+9eHTp0SM3Nzfrd734X9JCcF7UMJn/zK8z5cA75mx/k78gjf/FeYc4IifzNF/LXX3mfwKioqNDBgwf7fj506JCmTZuW78Xm5Nz4ysrKdNNNN6m5uTngEZlNnTpVbW1tkqS2tjaVlZUFPKLBTZ06VfF4XAUFBfrUpz4Vuu2bSqW0evVq3Xrrrbr55pslhXf7DjbWMG9fSZowYYKuvfZaNTU1hXbbuiJqGUz+5lfY84H8zT/yd+SQv/kXtX04zBlB/uYf+euPvE9gLFu2TC0tLdq/f7+SyaS2b9+uurq6fC922M6cOaOOjo6+/z/11FP9rh4cVnV1ddq6daskaevWrbrxxhsDHtHgzr1YJemJJ54I1fbNZrO64447NH/+fH3uc5/r+30Yt+9gYw3r9j1y5Ija29slSV1dXXr66ac1b968UG5bl0Qpg8nf/AtrPkjkbz6Rv8Egf/MvavtwWDOC/M0f8jcPRuJKob/4xS+yc+bMyc6ePTv7L//yLyOxyGF78803s9XV1dnq6ursggULQjneNWvWZC+99NJsIpHITp8+Pfvd7343e/To0eyHPvShbGVlZfZDH/pQ9tixY0EPM5vNDjzWv/u7v8suXLgwu2jRouzHPvax7FtvvRX0MPs8//zzWUnZRYsWZRcvXpxdvHhx9he/+EUot+9gYw3r9n3ttdeyNTU12UWLFmWrqqqyX/rSl7LZbDaU29Y1Uclg8tdf5O/IjzWs25f8DQ75658o5W82G60MJn/zh/z1XyybHeALegAAAAAAACGS96+QAAAAAAAA5IoJDAAAAAAAEHpMYAAAAAAAgNBjAgMAAAAAAIQeExgAAAAAACD0mMAAAAAAAAChxwQGAAAAAAAIPSYwAAAAAABA6P1/LH02/I9DQgoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x1080 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# reshape cvxpy solution\n",
    "u_cvx_arr = u_cvx.value.reshape(ig.shape)\n",
    "\n",
    "# print objectives\n",
    "print(\"CVX objective = {}\".format(obj.value))\n",
    "\n",
    "print(\"CIL objective = {}\".format(pdhg.objective[-1]))\n",
    "\n",
    "# show middle line profiles\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(u_cil.array[int(N/2)], label=\"CIL\")\n",
    "plt.plot(u_cvx_arr[int(N/2)], label=\"CVXpy\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "show2D([u_cil.array, u_cvx_arr, np.abs(u_cil.array - u_cvx_arr)], num_cols = 3, origin=\"upper\")\n",
    "\n",
    "# compare solution\n",
    "np.testing.assert_almost_equal(u_cil.array, u_cvx_arr, decimal=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9450173-cbdf-4a51-8df8-2ee08fae4b8d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:cil_devel_epaps] *",
   "language": "python",
   "name": "conda-env-cil_devel_epaps-py"
  },
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}