TotalVariation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "379f0a10-b77a-47a1-8858-64190d67c48a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cil.optimisation.functions import L2NormSquared\n",
    "from cil.utilities import dataexample\n",
    "from cvxpy import *\n",
    "from regularisers import tv\n",
    "from cil.optimisation.functions import TotalVariation\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from cil.utilities.display import show2D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a8cd78-a70f-48ba-8eff-6f6506eeb051",
   "metadata": {},
   "source": [
    "# Load data and resize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eec3a2c8-0734-44c1-b362-83e2f13d6ab8",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = dataexample.CAMERA.get(size=(32, 32))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcad0d4d-7db5-48aa-83ff-cb88b30cb2eb",
   "metadata": {},
   "source": [
    "# Setup and run tv denoising with cvxpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df8d2513-528b-4052-b5d8-4ce28eac4bc6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                    v1.1.17                                    \n",
      "===============================================================================\n",
      "(CVXPY) Dec 06 07:52:48 PM: Your problem has 1024 variables, 0 constraints, and 0 parameters.\n",
      "(CVXPY) Dec 06 07:52:48 PM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Dec 06 07:52:48 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Dec 06 07:52:48 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 06 07:52:48 PM: Compiling problem (target solver=SCS).\n",
      "(CVXPY) Dec 06 07:52:48 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS\n",
      "(CVXPY) Dec 06 07:52:48 PM: Applying reduction Dcp2Cone\n",
      "(CVXPY) Dec 06 07:52:48 PM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Dec 06 07:52:48 PM: Applying reduction ConeMatrixStuffing\n",
      "(CVXPY) Dec 06 07:52:48 PM: Applying reduction SCS\n",
      "(CVXPY) Dec 06 07:52:48 PM: Finished problem compilation (took 2.362e-02 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 06 07:52:48 PM: Invoking solver SCS  to obtain a solution.\n",
      "------------------------------------------------------------------\n",
      "\t       SCS v3.0.0 - Splitting Conic Solver\n",
      "\t(c) Brendan O'Donoghue, Stanford University, 2012\n",
      "------------------------------------------------------------------\n",
      "problem:  variables n: 2049, constraints m: 4098\n",
      "cones: \t  q: soc vars: 4098, qsize: 1025\n",
      "settings: eps_abs: 1.0e-05, eps_rel: 1.0e-05, eps_infeas: 1.0e-07\n",
      "\t  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1\n",
      "\t  max_iters: 100000, normalize: 1, warm_start: 0\n",
      "\t  acceleration_lookback: 10, acceleration_interval: 10\n",
      "lin-sys:  sparse-direct\n",
      "\t  nnz(A): 6018, nnz(P): 0\n",
      "------------------------------------------------------------------\n",
      " iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)\n",
      "------------------------------------------------------------------\n",
      "     0| 1.42e+01  5.23e-01  2.10e+02 -1.02e+02  1.00e-01  7.76e-03 \n",
      "   250| 1.47e-03  1.87e-04  3.05e-03  7.89e+00  1.00e-01  1.08e-01 \n",
      "   500| 5.33e-04  1.32e-06  5.43e-04  7.89e+00  1.00e-01  2.26e-01 \n",
      "   750| 2.65e-04  1.07e-06  2.53e-04  7.89e+00  1.00e-01  3.47e-01 \n",
      "  1000| 1.75e-04  1.48e-06  1.32e-04  7.89e+00  1.00e-01  4.77e-01 \n",
      "  1250| 1.11e-04  3.42e-07  7.90e-05  7.89e+00  1.00e-01  5.90e-01 \n",
      "  1475| 6.55e-05  1.24e-06  5.43e-05  7.89e+00  3.18e-01  6.90e-01 \n",
      "------------------------------------------------------------------\n",
      "status:  solved\n",
      "timings: total: 7.09e-01s = setup: 1.89e-02s + solve: 6.90e-01s\n",
      "\t lin-sys: 3.32e-01s, cones: 1.20e-01s, accel: 7.15e-02s\n",
      "------------------------------------------------------------------\n",
      "objective = 7.892653\n",
      "------------------------------------------------------------------\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 06 07:52:49 PM: Problem status: optimal\n",
      "(CVXPY) Dec 06 07:52:49 PM: Optimal value: 7.893e+00\n",
      "(CVXPY) Dec 06 07:52:49 PM: Compilation took 2.362e-02 seconds\n",
      "(CVXPY) Dec 06 07:52:49 PM: Solver (including time spent in interface) took 7.260e-01 seconds\n"
     ]
    }
   ],
   "source": [
    "# solution\n",
    "u_cvx = Variable(data.shape)\n",
    "\n",
    "# regularisation parameter\n",
    "alpha = 0.1\n",
    "\n",
    "# fidelity term\n",
    "fidelity = 0.5 * sum_squares(u_cvx - data.array)   \n",
    "regulariser = alpha * tv(u_cvx)\n",
    "\n",
    "# objective\n",
    "obj =  Minimize( regulariser +  fidelity)\n",
    "prob = Problem(obj, constraints = [])\n",
    "\n",
    "# Choose solver ( SCS, MOSEK(license needed) )\n",
    "tv_cvxpy = prob.solve(verbose = True, solver = SCS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b65944b-0c61-436b-b691-729e7ee6a414",
   "metadata": {},
   "source": [
    "# Use TotalVariation function from CIL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "58ed46f9-1894-4dfc-90e8-95be1b04c8a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialised GradientOperator with C backend running with  20  threads\n"
     ]
    }
   ],
   "source": [
    "# use TotalVariation from CIL (with Fast Gradient Projection algorithm)\n",
    "TV = alpha * TotalVariation(max_iteration=100)\n",
    "u_cil = TV.proximal(data, tau=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba2e34e0-7bed-4bd8-a91d-84cf323704e6",
   "metadata": {},
   "source": [
    "# Compare solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2a82667f-cc2c-4ab9-bf87-e52b4e89ccc3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CVX objective = 7.89286512743248\n",
      "CIL objective = 7.8956286907196045\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsuUlEQVR4nO3deXhc9X3v8fd3Nmm075Jly5ZsDNgGG4IwZCEYQhogiyFJi52UhLQJpYHe0DxtCn1uG0raPKTtzaW3kMsD2dukLk1IIDckaUhMCGGzTIzBGBt5lSxrs/ZlZjQz3/vHjB1Z1jKSRzpzRt/Xgx40Z34z+hyN9PXRmfP7/kRVMcYY434epwMYY4xJDyvoxhiTJaygG2NMlrCCbowxWcIKujHGZAmfU1+4oqJC6+vrnfryxhjjSjt37uxW1crJ7nOsoNfX19PU1OTUlzfGGFcSkSNT3WenXIwxJktYQTfGmCxhBd0YY7JESgVdRK4VkX0i0iwid01yf7GI/EhEXhGRPSLyifRHNcYYM50ZC7qIeIEHgeuAtcBWEVk7YdjtwOuqugHYBPwvEQmkOasxxphppHKEvhFoVtWDqhoBtgGbJ4xRoFBEBCgAeoBoWpMaY4yZVioFfSnQMu52a3LbeA8Aa4A24FXgM6oan/hEInKriDSJSFNXV9ccIxtjjJlMKtehyyTbJvbcfQ+wC7gaWAX8XER+raoDpz1I9WHgYYDGxkbr22tMGrz840eIdOxf+C88rvX2tL/MKbbozq2/lPVX/T5er/fsci1iqRT0VqBu3O1lJI7Ex/sEcJ8mmqs3i8gh4HzgpbSkHOfI/ldoefGHKY0V8SCBIB5/Ht6cPLyBIL7cPPw5+fhyCwgE88gJFhDICeILBAgEcgn4c/B47eIf4w7tLc1c9NJf4pEsOD5q/Sqtv7mHltV/yNrrP01xSZnTiVwnlYK+A1gtIg3AMWAL8JEJY44C7wJ+LSLVwHnAwXQGPam7+WXeceDL8/HUp4yplyhexsRHDB9RfETFSwwfcfESw0tcxn0kb6t4iXlyKHrfvTSsu2xeMxoDcHj7t6gRpeXm56itXzPtWFUl8TbXzFIbBeL53chUn3sy0UiYV5/6d/J++whv3f9PDO17gOcq30fte+6kfvUFc37exUZSWbFIRK4H7ge8wNdV9R9E5DYAVX1IRGqBbwJLSPws3Keq/z7dczY2Nupcpv7HxyJEw8OnbZtqD2KxKOHREcKjQ0RGhxkLDRMNjzAWHiYWHiEWHkEjI8THQsRjYxAbQ6MRiCc+P/l/iUeRWAQ0hkdjiEaR+O8+9yS3e4hxXmQvL9V9gss/+b9nvW/GzNahezcQ9gQ5/3++4HSUtDmw6xn6tv8f1vf9Ei9xduVdhlz+aTa84/321zMgIjtVtXGy+1Lq5aKqTwJPTtj20LjP24DfO5uQqfL4AwT8qV8RmVc8j2Em0fZ3q/ENHF3YL2oWpUN7XqQhfpgXz7vb6Shpteqid8JF76Sn/QjNP/4XVrf8F6XbP8bBX63g+HkfY8lbrmdp/Wpy/H6no2Ycx5pzZaueQC2FIy0zDzTmLHU8+28sUy+rr7rZ6SjzoqxmBRv/+MtEQn/Pb3/2NUpe/Tpv3/sF2PsFRjSH/Z6lnAjWEyo5B0/VeRTVraV25QVUlRSe1ekfN7OCnmbD+cs5t+dpp2OYLBePxWg4/iSv5zWyoWriVcTZJZCbx8Wb/ww+cAdHXn2WngNNxDvfILf/ICtDr1Hd9svEZRq7IKoeDlNNR2A54WA1FFThK6ohWLaE4oqlVC6po7C8FgnkOb1b88IKepppST2lPQMM9vdQWGzv0pv58caLP2UtJ2hZl12nW6Ylwor1V7Bi/RWnbY6Hhug+uoeew68SOf4Gvp79LB05QuHgXkoGBs68Jg8YJo8BXxmhnArCedXEC2rwFNWSU15HQUUdRVXLySmpBZ+7JrxbQU+zQNUqOAgdh/dSuOHtTscxWWpox3cZ1lzWXbXF6SiO8+QWUHXuZVSde+aVZRqN0NfVRndHC/2dxxjpbWOsvwOGOvCHuikY7qZ88LfUdPaSI2NnPL5XiunzVRLyFZH6tT+TmPDQ0LmbufiGz8z9+aZgBT3NimrPBWCgbT9YQTfzIDQ6zPm9v+T1kiu5NL/Q6TgZTXwBSpfUU7qkfsoxY7E4JwbDnOhuZ6jrKKETrUT725DB4/iH28kLd5IbGZzxa015veAkdwyPjKaUf7asoKdZ9YrzAQh3zctl+Maw95nvcTEj5L7Fjs7Twe/1UFMSpKakAc5pcDrOWbGLOtOssLiMHorw9h1yOorJUrr7v+imhDVve5/TUUyGsYI+D7p8S8gbtksXTfr193ZxwdDzNFe9B98s5mOYxcEK+jwYzKujPHLM6RgmC+375b8TkCjlb/1Dp6OYDGQFfR5Ei1ZQFe8mEg45HcVkmYJ93+eoLOWcDe9wOorJQFbQ54G3YiVeUTpa3nQ6iskiHUffZG3kVY4tfz/isV9dcyb7qZgHBUtWA9Dbss/hJCabHHr62wAsv/LjDicxmcoK+jyorEtcujja2exwEpNNqg8/wRu+NSxdOXFJX2MSrKDPg/KaOkY0B+2xSxdNehx8LdFZsX/1DU5HMRnMCvo8EI+HDm8NOYPWRtekR8dvvs1YFndWNOlhBX2e9OUupTTU6nQMkwXisRgrj/+E1/MaKcvyzorm7FhBnyehwhVUx9rReNzpKMbl9r7wU6o5wdgFv+90FJPhrKDPE09ZPUGJ0N1up13M2RluSnZW3GS9W8z0rKDPk2DVOQB0HX3D4STGzcZ3VgxaZ0Uzg5QKuohcKyL7RKRZRO6a5P6/FJFdyY/XRCQmIot6dYfSuvMAGG63yUVm7vY+8z2KrLOiSdGMBV1EvMCDwHXAWmCriJx2Iayq/pOqXqSqFwF3A79S1Z55yOsa1XWriakQ67Y2uuYs7H7UOiualKVyhL4RaFbVg6oaAbYBm6cZvxX4j3SEc7NATi4dnkp8A3YO3cxNf08X64ZesM6KJmWpFPSlwPhesK3JbWcQkTzgWuD7U9x/q4g0iUhTV1fXbLO6Tk+glsIRa6Nr5mbfduusaGYnlYI+2UJ6U6229H7gN1OdblHVh1W1UVUbKysrU83oWsP5y6mKHnc6hnGpvP0/pEVqrbOiSVkqBb0VqBt3exmTrqMNwBbsdMspWlJPKQMM9i/qtxPMHIRDI6wO7eFY1Tuts6JJWSo/KTuA1SLSICIBEkX7iYmDRKQYuBJ4PL0R3StQtQqAjsN7HU5i3ObQ7t+QI2MEVtrRuUndjAVdVaPAHcDPgL3Ao6q6R0RuE5Hbxg29EfhvVR2en6juU1R7LgADbfsdTmLcpveNZwCov/hqh5MYN/GlMkhVnwSenLDtoQm3vwl8M13BskH1ikQb3XCXXbpoZifv+Isc9SxlufVuMbNgJ+fmUWFxGb0U4e2zNromdfFYjIbR12gvvtjpKMZlrKDPs07fEvKG7dJFk7ojbzRRxDCy4m1ORzEuYwV9ng3m1VEeOeZ0DOMina89DUDtejt/bmbHCvo8ixatoCreTSQccjqKcQlv64t0UkZt/XlORzEuYwV9nnkrVuIVpaPFmnSZmWk8Tt3gLloKN9j152bW7CdmnhUsWQ1Ab8s+h5MYNzh+9M3EYhZLL3c6inEhK+jzrLIuceniaGezw0mMG7Tt/iUAles2ORvEuJIV9HlWXlPHiOagPXbpoplZ7PBzDJBH/ZpGp6MYF7KCPs/E46HDW0POoLXRNTOr6fsth4IX4PWlNOfPmNNYQV8AfblLKQ21Oh3DZLjeruOsiLcwUnOp01GMS1lBXwDhwuVUx9rReNzpKCaDHf7tLwAoOf9Kh5MYt7KCvgCkrIGgROhut9MuZmrhA88SVj8rN1zhdBTjUlbQF0Cw6hwAuo6+4XASk8nKTuzkYOBccnLznI5iXMoK+gIorUvM+Btut8lFZnIjQ/00jB2gr9KubjFzZwV9AVTXrSamQqzbLl00kzu461f4JUbeajvdYubOCvoCCOTk0uGpxDdwxOkoJkMN7v81cRXqL7rK6SjGxaygL5CeQC2FI9ZG10yuoGMHh3z1FJdWOB3FuJgV9AUynL+cquhxp2OYDBQdi7Aq9DrdpW9xOopxOSvoC0RL6illgMH+HqejmAxz8NXnyZMwvgZb0MKcnZQKuohcKyL7RKRZRO6aYswmEdklIntE5Ffpjel+gapVAHQc3utwEpNpevYmfl2WX3yNw0mM281Y0EXECzwIXAesBbaKyNoJY0qArwAfUNV1wO+nP6q7FSXb6A607Xc4ick0OW0vckyqqaytdzqKcblUjtA3As2qelBVI8A2YPOEMR8BHlPVowCq2pnemO5XXb8GgHDXQYeTmPk02/YOGo+zYng3x4psQWhz9lIp6EuB8ZdntCa3jXcuUCoiT4vIThH52GRPJCK3ikiTiDR1dXXNLbFLFRaX0UsR3j67Fj1bNf34EXrvXUFr82spP6aleTdlDMByW9DCnL1UCrpMsk0n3PYBlwDvBd4D/I2InHvGg1QfVtVGVW2srKycdVi36/QtIW/YLl3MVnmvfIsyBhj8z1uJRaMpPab91e0A1Fxo15+bs5dKQW8F6sbdXga0TTLmp6o6rKrdwDPAhvREzB6DeXWURSZ+60w2OH5kH2sjr7Lfdy5rxvawY9vfp/Q4Ofo8PRRRd876eU5oFoNUCvoOYLWINIhIANgCPDFhzOPAFSLiE5E84DLALueYIFq0gup4F5FwyOkoJs0OP/0tAAo++m/8Nu9tXPzmAxzZu3PGxy0d2MWR/PW2ILRJixl/ilQ1CtwB/IxEkX5UVfeIyG0icltyzF7gp8Bu4CXgq6qa+onERcJbsRKvKB0t1qQrm2g8Tu2RJ9jrX0dtw/nUfexhRiSXyPf/hLFIeMrHdR07RK12EK69bAHTmmyW0mGBqj6pqueq6ipV/YfktodU9aFxY/5JVdeq6gWqev885XW1guSli70t+xxOYtLpwKvPsyLewsC5HwSgoqaOQ5d9gdXRN2n6zt9O+bijuxILWpSv3bQQMc0iYH/nLaDKuvMBGO1sdjiJSafu5/6NiHo5/+qbT217y3WfoKnwXTQefoTmV34z6eOih37DiObQcIFd4WLSwwr6AiqvqWNUA2iPXbqYLWLRKOd0/IQ9+ZdTXF592n2rb/m/9Esh3sf/lHBo5IzHVvT+lgO5a/H5AwsV12Q5K+gLSDweOrw15AzaUnTZ4vXnfkQFfej6PzjjvuLyao6980s0xI/w8rf/6rT7+nu7aYgeZqh640JFNYuAFfQF1pu7jNJQq9MxTJqEdv4HA+Sx9srJu11suHoLL5W+l43H/o03djx1avvh3/4SjyiF575joaKaRcAK+gILFy6nOtY+6yniJvOMDPWztu9XvFF6FbnB/CnHrbnlAbqkgvwn/4zR4cHEY5ufZUy9rLrYJhSZ9LGCvsCkrIGgROhut9Mubvf60/9JvoTIb/zotOMKi8vovubL1Gkbr3zzzwEo6WrioP8cgvmFCxHVLBJW0BdYsOocALqOvuFwEnO2/Hu+RwflrLn82hnHXvCOD/Bi5Ye5vOu/2L39v1gV2UdPxSULkNIsJlbQF1hp3XkADLfbpYtu1tN5jHUjOzi45L14vN6UHrP+lvtpkVrOefp2AhIld+Xb5zmlWWysoC+w6rrVxFSIdVsbXTd785ffxidxaq6YtLHopIL5hQxf93/IIQJA/VtsQQuTXj6nAyw2gZxc2jyV+AaOOB3FnIWS5h9wwNvAqrWXzupx5298N8/vvZ1A9x4uqaiZp3RmsbKC7oCeQC2FI9ZG161aml/lvOg+XjjnTlbN4fFv/fg/pD2TMWCnXBwxnL+cquhxp2OYOWr91beIq7DyqlucjmLMaaygOyBeUk8pAwz29zgdxcySxuPUtf6I13M3ULW0wek4xpzGCroDcipXAtBx2FrGu82+l7ezTNsZOf/DTkcx5gxW0B1QUJ04shvqtvPobtP/wr8TUj9rrp5+MpExTrCC7oD8kkRXvsjA4loo2+3GImHO7f45e4reQWFxmdNxjDmDFXQHFCcvV4sPdzucxMzGnl//gFIG8W64yekoxkzKCroD8guKiagPHT7hdBQzC7Fd2+ilkLVX3Oh0FGMmZQXdAeLx0CdFeEetoLvFYH8P6waeZX/Fuwnk5Dodx5hJpVTQReRaEdknIs0ictck928SkX4R2ZX8mHohRQPAoLcYf6TP6RgmRW9s/y65MkbxZfZmqMlcM84UFREv8CDwbqAV2CEiT6jq6xOG/lpV3zcPGbPSqK+Y4Fif0zFMinL3fo9jUs15l1ztdBRjppTKEfpGoFlVD6pqBNgGbJ7fWNkvHCglP9rndAyTgq62w6wL7eLosvcjHjtLaTJXKj+dS4HxF0y3JrdN9FYReUVEfiIi6yZ7IhG5VUSaRKSpq2txX7IXzSmlSPudjmFScOz15/GIUrr+OqejGDOtVAq6TLJNJ9x+GVihqhuAfwV+ONkTqerDqtqoqo2VlZWzCppt4nkVFDNMdCzidBQzg/CJxOpSFUtXO5zEmOmlUtBbgbpxt5cBbeMHqOqAqg4lP38S8ItIRdpSZiFPfjkAfSc6HE5iZhLvbyWiXsqqlzkdxZhppVLQdwCrRaRBRALAFuCJ8QNEpEZEJPn5xuTz2jV50/AVJP69G+q1gp7p/ENtdHvKU16ZyBinzHiVi6pGReQO4GeAF/i6qu4RkduS9z8EfBj4UxGJAqPAFlWdeFrGjJNTXAXAcG+nw0nMTPJG2+nzV1PrdBBjZpDSAhfJ0yhPTtj20LjPHwAeSG+07JaX7OcS7rcj9ExXEu2krXC90zGMmZFdg+WQwrJEQR8btH4umSwWjVIZP8FYgR2fm8xnBd0hxeXWoMsNejpb8UsMT0ndzIONcZgVdIcEcnIZ1CAyaqsWZbKetoMA5JQvdziJMTOzgu6gAU8RvpAV9Ew21HkEgMKqemeDGJMCK+gOGvYWE4j0Oh3DTGOs9+SkolUOJzFmZlbQHTTqLyXPGnRltv5jjGgORSXlTicxZkZW0B0UCZRSELN+LpksMNxGl7fSmnIZV7CfUgfFcksp1gGnY5hpFIQ76A9UOx3DmJRYQXeQ5lcSlAijw4NORzFTKIt2EgrWOB3DmJRYQXeQ92SDru62GUYaJ4RDI1TQR6xwsm7RxmQeK+gO8heebNBl/VwyUXdb4pJFb6lNKjLuYAXdQcFkP5fRPivomaiv/RAAwfIVDicxJjVW0B2UX5rouBgZWNyrN2Wqka7DABTXWEE37mAF3UFFZYk326JDVtAzUay3FYBKm1RkXMIKuoMKSyqIqgcdtrVAMpEMHqOXQoL5hU5HMSYlVtAd5PF66ZdCPNbPJSPljrTR413ca98ad7GC7rBBTxF+K+gZqSjcyUCOXYNu3MMKusNGfCXkWj+XjFQe7yKSZwXduIcVdIeF/KXkRfucjmEmGBropYhh4kXLnI5iTMpSKugicq2I7BORZhG5a5pxl4pITEQ+nL6I2W0st5SiuDXoyjQnjiUWtvCX2aQi4x4zFnQR8QIPAtcBa4GtIrJ2inFfAn6W7pDZLJ5bRrEOEo/FnI5ixunvOAxAfqVdg27cI5Uj9I1As6oeVNUIsA3YPMm4PwO+D9i0x1mQ/HK8ogz22dqimSTUnZj2X1JT72wQY2YhlYK+FGgZd7s1ue0UEVkK3Ag8NN0TicitItIkIk1dXTaZBsBXkLgsbqCn3eEkZrx4XytxFSqW1DsdxZiUpVLQZZJtOuH2/cBfqeq05w1U9WFVbVTVxspKu74XIFCU+D4M93Q4nMSM5x1qo1tK8QdynI5iTMp8KYxpBca/M7QMmNjvtRHYJiIAFcD1IhJV1R+mI2Q2C5Yk+rmM9tuZqkwSHD1Or6+KKqeDGDMLqRT0HcBqEWkAjgFbgI+MH6CqDSc/F5FvAv/PinlqCpP9XMYG7Rx6JimOdNKVv9rpGMbMyoynXFQ1CtxB4uqVvcCjqrpHRG4TkdvmO2C2Ky5PFPSYNejKGBqPUxnvIpJf63QUY2YllSN0VPVJ4MkJ2yZ9A1RVbzn7WItHML+QUQ0gI9agK1P0dh+nTMag2CYVGXexmaIZoF+K8YZ6nY5hknqOJxa2yCm3SUXGXaygZ4AhbxGBiDXoyhSDyUlFBVUN0w80JsNYQc8AI/4SgtagK2OETxwFoHSJFXTjLlbQM0AkUEp+zPq5ZIz+ViLqo6zS3hQ17mIFPQNEc8soig84HcMk+YeP0+mpwOP1Oh3FmFmxgp4BNFhGoYwSDo04HcUA+aF2+v02pci4jxX0DODJrwBgoMdmi2aCkrFORoK2sIVxHyvoGcCf7OcyaP1cHBcdi1ChPUQLls482JgMYwU9A+QUJf68H+m1jotO624/ik/ieEpsUpFxHyvoGSC/NFHQwwM2/d9pfclJRbnlyx1OYszsWUHPACcbdEWHrEGX04a7EgtbFNXYNejGfaygZ4DissQRetwKuuPGehKTisprVzqcxJjZs4KeAXz+AP3k4xm1Bl1Ok4FjDGmQopJyp6MYM2tW0DPEgBTjC1uDLqcFho/T7bXVtIw7WUHPEMO+EnIiVtCdVhhuZyBQ7XQMY+bECnqGGPWXkGcNuhxXFusilGeTiow7WUHPEGOBEgqsn4ujQiNDlDFArNAmFRl3soKeIWLBCkp0AI3HnY6yaHW1Ja5B95XawhbGnVIq6CJyrYjsE5FmEblrkvs3i8huEdklIk0i8o70R81ukl9GQKIMDfY5HWXR6m9PFPRghU0qMu40Y0EXES/wIHAdsBbYKiJrJwz7BbBBVS8C/gj4appzZr1TDbpOWD8Xp4x2J65BL6mxa9CNO6VyhL4RaFbVg6oaAbYBm8cPUNUhVdXkzXxAMbOSU5y4VG64zwq6U6K9LQBULLVZosadUinoS4GWcbdbk9tOIyI3isgbwI9JHKWbWchNNuga7bMGXU7xDh6jhyJyg/lORzFmTlIp6DLJtjOOwFX1B6p6PnAD8IVJn0jk1uQ59qauLmtENV5BWeLa58iATf93Su5IOye8trCFca9UCnorMP5t/2VA21SDVfUZYJWIVExy38Oq2qiqjZWVNhtvvKJkg66Y9XNxTFGknaEcm1Rk3CuVgr4DWC0iDSISALYAT4wfICLniIgkP38LEACsMcksFBSVElEvOmzfNqdUxLqJ5C9xOoYxc+abaYCqRkXkDuBngBf4uqruEZHbkvc/BHwI+JiIjAGjwE3j3iQ1KRCPh34pwhuygu6Egb4TFMkoWmwLWxj3mrGgA6jqk8CTE7Y9NO7zLwFfSm+0xWfIU4zfGnQ54sSxZooAv00qMi5mM0UzyLC/hFzr5+KIgY7Ewhb5lSscTmLM3FlBzyBhfwkF0T6nYyxKoROJSUWltrCFcTEr6BkkmltGoVqDLifE+1qIqoeKGjtCN+5lBT2DxIPlFOkw0bGI01EWHd9QG91ShteX0ttKxmQkK+gZxJNfjkeU/p5Op6MsOnmj7fT6bVKRcTcr6BnEW5iYizXUY/1cFlrJWAcjubawhXE3K+gZ5GQ/l+E+O0JfSPFYjMr4CcZsUpFxOSvoGSRYkph2Huq3gr6QerraCEgUKbFr0I27WUHPIIXJBl1jg9a4bCH1tB0EIFBmBd24mxX0DFJcnijo8WFr0LWQhrsOA1BYXe9oDmPOlhX0DJKTm8eQBpGRHqejLCrhE4l2/+W1qxxOYszZsYKeYQY8RfisQdfC6m8lpH5Kyq11rnE3K+gZZshbTCDS53SMRSUw3EaXpxLx2K+DcTf7Cc4wo/4Sgtaga0HlhzroD9ikIuN+VtAzTCRQSmGs3+kYi0pptJPRoF2DbtzPCnqGieWWUWQNus7a/p3badm/a8ZxY5EwFdpLrKB2/kMZM8+soGeavDLyJMzo8KDTSVwrOhah5kcfofw7v8euX2ybdmz38cN4RPHYwhYmC1hBzzCegsTi2f3Wz2XOml/5NUWMEJYcLnzmNl74jy9OObb3+CEAcsuXL1Q8Y+aNFfQM4y9MFPShnnaHk7hX72s/ByD0iafYnf9WLt/3JV74yqeIRaNnjB3pSqxUVGyTikwWSKmgi8i1IrJPRJpF5K5J7v+oiOxOfjwnIhvSH3VxCBYnCvqoNeias8K25zjgXcmSFeex/rM/4oWqm7i881F2f/n9jAyd/obzWG8rAOVLbaUi434zFnQR8QIPAtcBa4GtIrJ2wrBDwJWquh74AvBwuoMuFnmlickt4QHr5zIXo8ODnBveQ1fFZQB4fT4u//TDvHj+Xawffp62/30VXW1HTo33DLQyQD4FRaVORTYmbVI5Qt8INKvqQVWNANuAzeMHqOpzqnpyufoXgGXpjbl4FJUlenJHrUHXnDTvfIqARMk7/12nbb9sy928duVD1EZbiT18NYf2vAhAzshxTngqnYhqTNqlUtCXAi3jbrcmt03lj4GfnE2oxayotJKYCjpi0//nYmjvLxlTL+dc+u4z7ttw9RbaPvgDPMSpfHQzu5/+PoXhDgZybMq/yQ6pFHSZZJtOOlDkKhIF/a+muP9WEWkSkaauLjsCnYzH66VfCvGMWoOuuajofJ43A2vIKyiZ9P5zNrwd/eRTdPhqWLv9k9RFjxDKs0lFJjukUtBbgfEX6S4D2iYOEpH1wFeBzao66eGlqj6sqo2q2lhZaX/mTmXQU4w/ZAV9tvp7OlkVbaZ/ydumHVe9bBXVn9nOnrxGAhIjXjTdH5zGuEcqBX0HsFpEGkQkAGwBnhg/QESWA48BN6vq/vTHXFyGvcXkjPXOPNCc5uCOn+ARpXTdNTOOLSgqZd1nf8xLF97D6vfctgDpjJl/vpkGqGpURO4AfgZ4ga+r6h4RuS15/0PA3wLlwFdEBCCqqo3zFzu7hQMllI4edTqG64Tf3M6I5rDq4itTGu/zB9j4oT+f51TGLJwZCzqAqj4JPDlh20PjPv8k8Mn0Rlu8IjllFAy/5nQM16k98RLNwfWsD+Q6HcXMwdjYGK2trYRCIaejZITc3FyWLVuG3+9P+TEpFXSzsOLBMkp0gHgshsfrdTqOK7S3HGC5HqOt7iano5g5am1tpbCwkPr6epJ/6S9aqsqJEydobW2loaEh5cfZ1P8MJHnl+CTOYL+9MZqqo02JK2WrNrzH4SRmrkKhEOXl5Yu+mAOICOXl5bP+a8UKegbyJfu5DFg/l5TJ4V/RSxH1ay51Ooo5C1bMf2cu3wsr6BkokCzow1bQU6LxOCv6mzhY8BY7RWUWNSvoGShYklgOLdRvk69ScfTNV6iih3j9O52OYrJAe3s7W7ZsYdWqVaxdu5brr7+e/fv3c8EFFwDw9NNP8773vc/hlJOzN0UzUEFZYir6mPVzScnxl3/KCmDpJdc7HcW4nKpy44038vGPf5xt2xKLo+zatYuODnesT2AFPQMVlycbdA11O5zEHXJafk2bVFHbsMbpKCZN/u5He3i9Lb1LMa6tLeLz71837Zjt27fj9/u57bbfTTa76KKLOHz4cFqzzBc75ZKBgnmFhNSPWIOuGUXHxlg1sotjJfZmqDl7r732GpdcconTMebMjtAzkHg89EsRXuvnMqMDu3/DeQzjWXWV01FMGs10JG0mZ0foGWrQW4I/bP1cZtKTXG6u4dJrHU5issG6devYuXOn0zHmzAp6hhr1FRO0Bl0zKmx7lkOeesqq62YebMwMrr76asLhMI888sipbTt27ODIkSPTPCpzWEHPUJGcUvJj/TMPXMRCo8OcE9pDZ3K5OWPOlojwgx/8gJ///OesWrWKdevWcc8991BbW3vauF/84hcsW7bs1Mfzzz/vUOLT2Tn0DDWWU0ZRf3rf5c82bzb9ggtljOB5VzsdxWSR2tpaHn300TO2v/ZaomHepk2bGB0dXehYKbEj9AyleeUUyQiRsHWem8rg3l8QVQ+rLv09p6MYkxGsoGcoT345AAMn3DGhwQkVnc/RHDif/KIyp6MYkxGsoGcof7Kfy2Cv9XOZTH9PN6vG3qSv+q1ORzEmY1hBz1A5RYl+LsO9doQ+mQNNP8UrSnEKy80Zs1hYQc9QeaWJgh4esH4uk4m8uZ1RDbDqYptQZMxJVtAzVGGyQVd00Pq5TKbmxIs0By8kkBt0OooxGSOlgi4i14rIPhFpFpG7Jrn/fBF5XkTCIvIX6Y+5+BQnC3p82Ar6RJ3HDlMfb2Fk6RVORzFZaLL2uSLCvn37Tht355138o//+I889thjvOtd7zq1/dlnn+Wiiy4iGo0udPSZC7qIeIEHgeuAtcBWEVk7YVgP8D+Af057wkXKH8hhgHw81qDrDIeTy81VbHi3w0lMtjnZPnfTpk0cOHCA119/nS9+8Yts2rTpVDtdgHg8zve+9z1uuukmPvjBD5Kbm8t3v/tdotEon/70p/nKV76Cz7fw03xS+YobgWZVPQggItuAzcDrJweoaifQKSLvnZeUi9SAFOGzfi5nOvg0fRTQsM6ucMlaP7kL2l9N73PWXAjX3TftkKna5/7Lv/wLN910E5///OcBeOaZZ6ivr2fFihUA/Ou//ivXXHMNe/bs4dJLL+Vtb3sbAPfccw8HDhzg2LFjtLS08LnPfY5PfepT3HzzzXz4wx9m8+bNAHz0ox/lpptu4gMf+MBZ7WIqp1yWAi3jbrcmt5l5NuQtIRCxgj6exuMs72/ikC03Z+bBVO1z169fj8fj4ZVXXgFg27ZtbN269dT9K1eu5KabbuKBBx7gS1/60mmP3b17Nz/+8Y95/vnnuffee2lra+OTn/wk3/jGNwDo7+/nueee4/rrz36BllSO0CdbqVTn8sVE5FbgVoDly5fP5SkWlVCghMKwXYc+3tHmV1lBN0eX23JzWW2GI2knbN26lW3btrFu3Toef/xx7r333lP3xeNxnnrqKQoKCjhy5AgVFRWn7tu8eTPBYJBgMMhVV13FSy+9xA033MDtt99OZ2cnjz32GB/60IfScoomlSP0VmB8K7tlQNtcvpiqPqyqjaraWFlZOZenWFQigVIKrUHXKdFojNadifPnSy+xdrkm/aZrn7t161YeffRRnnrqKdavX09VVdWp+x588EEuuOACvva1r3H77bej+rtjXpHTj4lP3r755pv5zne+wze+8Q0+8YlPpCV/Kv8k7ABWi0gDcAzYAnwkLV/dTCseLKW4dxCNxxFP5l5hqqpExiKEhocJjw4RDo0QHh1iLDTCWGiEaHiYWCREfGwEHQuhY6MQDcPYKBINQTSExMJ4oqN4oiN4oyP4oyME4qME4iFydJSghsgjxNslRodUUNtgCyCY9Lv66qv567/+ax555BE+9alPAYn2uSMjI1x55ZWUl5dz1113ceedd556THt7O1/+8pd56aWXqKys5JFHHuGrX/3qqcc//vjj3H333QwPD/P0009z332Jvz5uueUWNm7cSE1NDevWpefnecaCrqpREbkD+BngBb6uqntE5Lbk/Q+JSA3QBBQBcRG5E1irqtYu8GzkVZAjY3Teuwpl8nNf47eqTNwy1djTz5jpJI9QBUGTHyS/fuL2755NCWiEXCLkSIycVPZpEqMaICwBwgQIS5CwJ8iYN4/RQDlDvjxi/jzUnw+BfCSQT+maq6jO4H/gjHudbJ975513ct9995Gbm0t9fT33338/kDhKv/vuu7nxxhtPPeazn/0sn/vc5zh51uH+++/niiuu4EMf+hAAGzdu5L3vfS9Hjx7lb/7mb0614q2urmbNmjXccMMN6cs//k+DhdTY2KhNTU2OfG23OH5kHy0/+iISGzu17eSrdebLpif/m2Ts1K+xJD+Q35V7QVER5FT5FhBJjkuU9VOfe3NQfx7iz0X8eUggiNcfxJOTjy8niC8n79SHPycPf24QXyCPQDCfnNw8fP6cjP7rwyycvXv3smZNdi30fc8991BQUMBf/MWZ03NGRka48MILefnllykuLp708ZN9T0Rkp6o2Tjbe+qFnsCUrzmPJHd9yOoYxJs2eeuop/uiP/ojPfvazUxbzubCCbowx8+See+6ZdPs111zD0aNH0/717G9dY0zGcOoUcCaay/fCCroxJiPk5uZy4sQJK+okivmJEyfIzc2d1ePslIsxJiMsW7aM1tZWurqsZTQk/oFbtmzZrB5jBd0YkxH8fj8NDQ1Ox3A1O+VijDFZwgq6McZkCSvoxhiTJRybKSoiXcCROT68AsiGpXyyYT9sHzKD7UNmWIh9WKGqk3Y3dKygnw0RaZpq6qubZMN+2D5kBtuHzOD0PtgpF2OMyRJW0I0xJku4taA/7HSANMmG/bB9yAy2D5nB0X1w5Tl0Y4wxZ3LrEboxxpgJrKAbY0yWcF1BF5FrRWSfiDSLyF1O55kLETksIq+KyC4RccWyTSLydRHpFJHXxm0rE5Gfi8ibyf+XOplxJlPswz0iciz5WuwSkeudzDgTEakTke0isldE9ojIZ5LbXfNaTLMPrnktRCRXRF4SkVeS+/B3ye2Ovg6uOocuIl5gP/BuoJXEAtZbVfV1R4PNkogcBhpV1TWTKETkncAQ8G1VvSC57R+BHlW9L/mPa6mq/pWTOaczxT7cAwyp6j87mS1VIrIEWKKqL4tIIbATuAG4BZe8FtPswx/gktdCRATIV9UhEfEDzwKfAT6Ig6+D247QNwLNqnpQVSPANmCzw5kWBVV9BuiZsHkzcHKNvG+R+KXMWFPsg6uo6nFVfTn5+SCwF1iKi16LafbBNTRhKHnTn/xQHH4d3FbQlwIt42634rIfhCQF/ltEdorIrU6HOQvVqnocEr+kQJXDeebqDhHZnTwlk7GnKiYSkXrgYuBFXPpaTNgHcNFrISJeEdkFdAI/V1XHXwe3FXSZZJt7zhn9zttV9S3AdcDtyVMBxhn/F1gFXAQcB/6Xo2lSJCIFwPeBO1V1wOk8czHJPrjqtVDVmKpeBCwDNorIBQ5Hcl1BbwXqxt1eBrQ5lGXOVLUt+f9O4AckTiW5UUfyfOjJ86KdDueZNVXtSP5ixoFHcMFrkTxn+33gO6r6WHKzq16LyfbBja8FgKr2AU8D1+Lw6+C2gr4DWC0iDSISALYATzicaVZEJD/5RhAikg/8HvDa9I/KWE8AH09+/nHgcQezzMnJX76kG8nw1yL5ZtzXgL2q+uVxd7nmtZhqH9z0WohIpYiUJD8PAtcAb+Dw6+Cqq1wAkpcy3Q94ga+r6j84m2h2RGQliaNySCwB+F037IOI/AewiUR70A7g88APgUeB5cBR4PdVNWPfdJxiHzaR+BNfgcPAn5w8B5qJROQdwK+BV4F4cvNfkzgH7YrXYpp92IpLXgsRWU/iTU8viQPjR1X1XhEpx8HXwXUF3RhjzOTcdsrFGGPMFKygG2NMlrCCbowxWcIKujHGZAkr6MYYkyWsoBtjTJawgm6MMVni/wNuaaaroFpyIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAE4CAYAAACt/wPcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABVsUlEQVR4nO3de3QU553n/0/rLiEuBiTcIAhWGnMRCNlGMbmMb6yCoyQiNnsYvE4sL0k0Is4662Q8Yc8Zjz2T2URMJuc4uzAhmhBHzolRlmRAJxtQTLDxJTsWljH4QmwrDthCFiAB4iIktdTq3x/80CAj1VNSd1d3V71f53AO0vOtp54qVX+69aguvnA4HBYAAAAAAEACS4n3AAAAAAAAAEyYwAAAAAAAAAmPCQwAAAAAAJDwmMAAAAAAAAAJjwkMAAAAAACQ8JjAAAAAAAAACS+iCYzGxkbNnz9fgUBANTU10RoTAGAU5C4AOIvcBYDE4QuHw+HxLBgKhXT99ddrz549KigoUGlpqbZt26ZFixZFe4wAAJG7AOA0chcAEsu4z8DYv3+/AoGACgsLlZGRobVr16qhoSGaYwMAXIHcBQBnkbsAkFjGPYHR1tam2bNnD31dUFCgtra2qAwKAHA1chcAnEXuAkBiSRvvgiNdeeLz+a76Xm1trWprayVJb731lubNmzfeVTpqpG0ZT40TfSC+TFdh2blKa5xXco1ZNMbi1Fjt+OCDD9TZ2XnV9++8884Rv29l+vTpamxsjNbQYoLcJXdxCbkbP+SuOXdfeeWVmI/LrrQ064/6AwMDxj6ysrIiXs/g4KCxj2AwGPF6UlNTjX309PQYa9LT0y3b7byPmF63dvaJndd+Sor136NN2yJJ3d3dxhpTP3aOJdPPT4rOe7Sd/Wba/5MnTzb2cfHiRdtjisSECRM8lbsm457AKCgoUGtr69DXx44d08yZM6+qq6qqUlVVlSTphhtu0LPPPjveVUaV6cVuapfsvQhN/dgJ2mh8qOcDe+yYAru/v9/Yh52aaLwRhkIhY43pA4SdNyinPmyvXr16xO93dnbq5ZdfHlNfpaWl0RhSTJG75C4uSZTctZN1djLTNJZojDVayF1z7vp8PmPOROM93U7GXHPNNZbtp06dMvYRCASMNdOnT7dsv3DhgrGP9957z1gzbdo0y3Y7v3QePnzYWJOfn2/ZbmdSwPTa7+vrM/Zh57VvmmCaMWOGsY+mpiZjjakfO79Im35+kr19a2Jnv5kmslauXGns4/XXX7c9pkhkZmaO+H235q7JuC8hKS0tVUtLi44cOaJgMKj6+npVVFREc2wAMG7hcHhM/5IBuQsgkZG7AOAsN+auybjPwEhLS9OmTZu0cuVKhUIhrVu3TkVFRdEcGwCMm1tC+krkLoBERu4CgLPcmLsm457AkKTy8nKVl5dHaywAEBVummX+MHIXQCIidwHAWW7OXSsRTWAAQKLyYqADQDyRuwDgLC/mLhMYAFzJi4EOAPFE7gKAs7yYu+O+iScAAAAAAIBTOAMDgCt5cUYaAOKJ3AUAZ3kxdz07gWHn2dlO9BEt0Th4E2l7EoUXQ8Et+NklHqdylyxLbsn02k2msTqB/XFJamqqJk6caFlj2leDg4PG9djZ3ykpkZ9snZqaaqzp7++3bM/MzDT20dvba6zp6emxbJ88ebKxj0mTJhlrzpw5E3Efpu0xbYskDQwMGGv6+vos26+77jpjH3aYfsahUMjYh2mskvk93M7xaIdpPa+99pqxj3Pnzhlr7PwMTWbOnDlqmxdz17MTGADcy6t3ZQaAeCF3AcBZXs1dJjAAuJIXAx0A4oncBQBneTF3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7kxUAHgHgidwHAWV7MXSYwALiOV29qBADxQu4CgLO8mruRP1sJABLQ5VC3+w8AEBlyFwCcFYvcbWxs1Pz58xUIBFRTUzPiOh988EEFAgEVFxfrwIEDxmVPnz6tsrIyzZs3T2VlZUOPKv7FL36hkpKSoX8pKSk6ePCg5fg4AyPGTM8YtsPOwRaN9eBqdva9qSaRPqQl0lhizUvbiuhz6vghu6/mxdxNpPFGwi3bEamUlBRNnDjRssb02g+FQsb1dHV1GWvS0qw/6mdnZxv7sPNz7e3ttWzPzc019tHd3W2sSU9Pt2yfMmWKsQ/Tz0aSOjo6LNunTp0acR9ZWVnGPuzo6emJeD12tsd0HAwMDBj7uHDhgrEmJcX67+umY1qy995qWs+pU6eMfdh5Dfb19UU0DkmaOXPmqG3Rzt1QKKQHHnhAe/bsUUFBgUpLS1VRUaFFixYN1ezevVstLS1qaWlRU1OT1q9fr6amJstla2pqtGLFCm3YsEE1NTWqqanRxo0bde+99+ree++VJL3++utatWqVSkpKLMfIGRgAXIm/BAKAs8hdAHBWtHN3//79CgQCKiwsVEZGhtauXauGhoZhNQ0NDbrvvvvk8/m0fPlydXV1qb293XLZhoYGVVZWSpIqKyu1c+fOq9a9bds23XPPPcYxMoEBwHXGGuZ8kAaAyJC7AOCsWORuW1ubZs+ePfR1QUGB2trabNVYLXvixAn5/X5Jkt/v18mTJ69a9y9/+UtbExhcQgIAAAAAgMt1dHRo2bJlQ19XVVWpqqpq6OuRJjk+fEnOaDV2lh1NU1OTcnJytHjxYmMtExgAXIm/7gGAs8hdAHDWWHM3Ly9Pzc3No7YXFBSotbV16Otjx45ddQ+O0WqCweCoy86YMUPt7e3y+/1qb29Xfn7+sD7r6+ttnX0hcQkJAJeKx12Zz549q89//vNaunSpioqK9MQTT0R7swAgYZG7AOCsaOduaWmpWlpadOTIEQWDQdXX16uiomJYTUVFhZ588kmFw2G99NJLmjx5svx+v+WyFRUVqqurkyTV1dVp1apVQ/0NDg5q+/btWrt2ra1t5gwMAK4Uj7syb968WYsWLdJvfvMbdXR0aP78+br33nuVkZER1bEAQCIidwHAWdHO3bS0NG3atEkrV65UKBTSunXrVFRUpC1btkiSqqurVV5erl27dikQCCgnJ2do4ni0ZSVpw4YNWrNmjbZu3ao5c+Zo+/btQ+t8/vnnVVBQoMLCQntjjOoWA0CCiHagX3lnZUlDd1a+8oO0z+fT+fPnFQ6HdeHCBU2dOtXW474AwA3IXQBwViwu3SsvL1d5efmw71VXVw/93+fzafPmzbaXlaRp06Zp7969Iy5z22236aWXXrI9PhIegOvE4g73I91ZuampaVjN17/+dVVUVGjmzJk6f/68fvnLX9p6vjcAJDtyFwCc5dUnOjGBMYpke/Nz6uC1eyfZZGBnn0WrxglOjTVRttdkrOOMxl2Zf/e736mkpETPPPOM3n33XZWVlekv/uIvNGnSpDGOPvlEIxuSLXdNovVaIXfHV+OEZBqrE8jdS3w+n1JTUy1rTGeJ2LkEpqenx1hjGkdOTo6xj4sXL0a8nszMTGMfds6cOXv2rGW7ne2ZMmWKsebKs3xGcvvttxv72LFjR8TjOH78uLHmxIkTlu379u0z9mHn5xMKhSzb7byHDwwMGGv6+/sjGockZWdnG2uCwaBl+x133GHsY9euXcaavLw8y/bBwUFjH1a89B5zGRMYAFwpHndlfuKJJ7Rhwwb5fD4FAgFdd911euutt/Sxj31sbIMHgCRE7gKAs7w4geGuP3cBwP8vHndlnjNnztD1fSdOnNDbb79t+4ZEAJDsyF0AcFYsnv6U6DgDA4ArxeOuzI888ojuv/9+LVmyROFwWBs3btT06dOjOg4ASFTkLgA4yy2TEmPBBAYA14nVLLPprswzZ87U008/HfX1AkCiI3cBwFluOqtiLJjAAOBKXgx0AIgnchcAnOXF3GUCA4AreTHQASCeyF0AcJYXc5ebeAIAAAAAgITHGRgAXMmLM9IAEE/kLgA4y4u568oJDJ/PF5UaJ9g56KIxVjvrsVPzxz/+0bK9r6/P2MeSJUss2zMyMox92GHansHBwYj7SKT1ROtnbGe8ycCLgZ7oEiV3oyUaxxi5O/Y+Emk9TuVusuRZsowz1sLhsPF12d/fb9lu5zWZnp5urAmFQpbtqampxj66urqMNZGOw+5YTK+VY8eOGfs4fvy4sWZgYMCy/dChQ8Y+srOzLduvv/56Yx+TJk0y1pw6dcqy/ZprrjH2ce7cOWON6Zi18/PLzMw01pheO2lp5l9fc3JyjDXnz5+3bJ82bZqxDzvb3NnZadluJzc/8pGPRLS827hyAgOAt3n1rswAEC/kLgA4y6u5ywQGAFfyYqADQDyRuwDgLC/mLhMYAFzJi4EOAPFE7gKAs7yYuxFNYMydO1cTJ05Uamqq0tLS1NzcHK1xAUBE3Bro5C6AREXuAoCz3Jq7ViI+A+PZZ5/V9OnTozEWAIgaNwc6uQsgEZG7AOAsN+fuaLiEBIDrePWmRgAQL+QuADjLq7mbEsnCPp9Pn/70p3XTTTeptrY2WmMCgIhdDnW7/5IFuQsgUZG7AOAst+aulYjOwPjDH/6gmTNn6uTJkyorK9OCBQt0yy23DKupra0dCnvTc3ARX7/4xS+MNX/3d38X8Xry8vIs2zds2GDs4wtf+IKxJhovUtPzxiXzs82jFRamfuysJ5mCK9KxJtO2jkWi5q7P53NkPdFg59hwanvI3auRu/FD7o5srLnb39+vU6dOWfZpOoZzc3ON48rMzDTWBINBy/be3l5jH11dXcaa7OzsiPuwc/ykpUV+8vjAwEDEfdgZ68WLFy3bn3vuOWMfc+fONdYsWbLEsv306dPGPrq7u4016enplu05OTnGPuzku2m/2VmPnZ+x3++3bG9razP2MXXqVGON6TV25swZYx9W3Jq7ViI6A2PmzJmSpPz8fN11113av3//VTVVVVVqbm5Wc3Mz1w4CcIxbZ6TJXQCJity9lLvJNKELILm5NXetjHsCo7u7W+fPnx/6/9NPP63FixdHbWAAEAk3Bjq5CyCRkbsA4Cw35q7JuM/DOnHihO666y5Jl07T+S//5b/ozjvvjNrAAADDkbsA4CxyFwASy7gnMAoLC3Xo0KFojgUAosJNs8xXIncBJCpyFwCc5dbcNeExqgBcyYuBDgDxRO4CgLO8mLtMYABwJS8GOgDEE7kLAM7yYu4ygQHAlbwY6AAQT+QuADjLi7nLBAYAV/JioANAPJG7AOAsL+auKycwovH87cHBQWNNSor5KbSmgyqRnhX+ox/9KOI+QqGQsebMmTOW7T/5yU+MfSxatMhYM2/ePMt2Oz/jYDBorDGxEyx2xmKqsdOHnbE4FYR2xjteXr2pUawkUk5F4xhOpO0hd682MDBgrDH1E62sI3ftI3eHM31GNO2rlpYW4zquv/56Y01vb69le09Pj7GPKVOmGGtmzpxp2W7n2Hv//feNNXPmzLFsLyoqMvbxb//2b8YaJ47lvr4+Y01HR4exxu/3W7YfPXrU2Edubq6xxvQzzMzMNPbR399vrDH1Y6cP03EvSTNmzLBs37dvn7EPO/vN9DpNSzP/Ot7d3T3i92OVu42NjfrGN76hUCikr3zlK9qwYcNV6/3GN76hXbt2KScnRz/72c904403Wi57+vRp/eVf/qWOHj2quXPn6v/8n/+ja665RpL02muv6a/+6q907tw5paSk6OWXX1ZWVtao4zP/Bg4ASciLz8UGgHgidwHAWdHO3VAopAceeEC7d+/W4cOHtW3bNh0+fHhYze7du9XS0qKWlhbV1tZq/fr1xmVramq0YsUKtbS0aMWKFaqpqZF06Y8WX/ziF7Vlyxa9+eab2rdvn9LT0y3HyAQGAFfigzQAOIvcBQBnRTt39+/fr0AgoMLCQmVkZGjt2rVqaGgYVtPQ0KD77rtPPp9Py5cvV1dXl9rb2y2XbWhoUGVlpSSpsrJSO3fulCQ9/fTTKi4u1tKlSyVJ06ZNU2pqquUYmcAA4Ep8kAYAZ5G7AOCsaOduW1ubZs+ePfR1QUGB2trabNVYLXvixImhy538fr9OnjwpSXrnnXfk8/m0cuVK3Xjjjfqnf/on4xhdeQ8MAODDMQA4i9wFAGeNNXc7Ojq0bNmyoa+rqqpUVVVl2d+H7x02Wo2dZT9sYGBAL774ol5++WXl5ORoxYoVuummm7RixYpRl+EMDACuM9bZaLvh39jYqPnz5ysQCAxdu3el73//+yopKVFJSYkWL16s1NRUnT59OtqbBwAJh9wFAGeNJ3fz8vLU3Nw89O/KyQvp0lkTra2tQ18fO3bsqpv0jlZjteyMGTPU3t4uSWpvb1d+fv5QX7feequmT5+unJwclZeX68CBA5bbzQQGAFeK9gdpOzc1evjhh3Xw4EEdPHhQ3/ve93Trrbdq6tSpsdpEAEgo5C4AOCvauVtaWqqWlhYdOXJEwWBQ9fX1qqioGFZTUVGhJ598UuFwWC+99JImT54sv99vuWxFRYXq6uokSXV1dVq1apUkaeXKlXrttdd08eJFDQwM6LnnnjM+9YxLSAC4UrRPZb7yxkSShm5MNFrIbtu2Tffcc09UxwAAiYzcBQBnRTt309LStGnTJq1cuVKhUEjr1q1TUVGRtmzZIkmqrq5WeXm5du3apUAgoJycHD3xxBOWy0rShg0btGbNGm3dulVz5szR9u3bJUnXXHONvvnNb6q0tFQ+n0/l5eX67Gc/az3GqG4xALjUSDcmampqGrH24sWLamxs1KZNm5waHgC4DrkLAM4rLy9XeXn5sO9VV1cP/d/n82nz5s22l5UuPV1k7969Iy7zxS9+UV/84hdtjy/pJjBMNwKxyzRbZWc9g4ODEY/Dznqisc0vvviisWZgYMBYk5eXZ9k+ZcoUYx+33nqrZbud/bpjxw5jzfLlyy3bFy5caOxjwoQJxhrTsWRne6JRY6cPp26wFq1tjkQ8bmp02W9+8xt98pOf5DTmD4nGayUayN2rkbvjqyF3hyN3L0lPT9eMGTMsa/r6+izbe3p6jOuxs79NP3O7l/KYXL6efTQ333yzsY8f/ehHxpqsrCzL9t27dxv7sMP0HmBnv5keBWlql6QzZ84Ya1544QXL9tzcXGMfpuPVDjvvZ+np6cYa0/Fm57Vx8eJFY82xY8cs2+1sT2dnp7Gmt7fXsv0v//IvjX28+uqro7Z58ebJSTeBAQB2jDXQL9/UaDR2bmp0WX19PacxA/AcchcAnOXFCQxu4gnAdWJxN3w7NzWSpLNnz+q5554bujkRAHgBuQsAzorV058SHWdgAHCleNzUSLp0mv2nP/1pW6fBA4CbkLsA4Cy3TEqMBRMYAFwpFoFuuqmRJN1///26//77o75uAEh05C4AOIsJDABwCS8GOgDEE7kLAM7yYu4ygQHAlbwY6AAQT+QuADjLi7nLBAYA13HTjYoAIBmQuwDgLK/mrisnMOz8IKPxfOdoPHP95z//ubGP3/3ud8aa8+fPW7a///77xj7sPCfa9PztoqIiYx8ZGRmW7abnJUv2ns38+uuvW7bbef72vHnzjDWm7bEjGseSnT7sMB370XptxJoXAz2eyN2rkbtXS6TctXO8hUIhy3Zydzhy95L09HRde+21ljXd3d2W7Z2dncb1nD171lhjet2a8kWSpk6daqxpb2+3bP/+979v7MPOMXzw4EHL9vT0dGMf0Xi/ssO0HjuZmpJifmCkab+dO3fO2Ed2draxZtq0aZbtdo6lixcvGmtM+83OPjFlt2Qeb0FBgbGP06dPG2seeughy/ZTp04Z+7Dixdx15QQGAHgx0AEgnshdAHCWF3OXCQwAruTFQAeAeCJ3AcBZXsxdJjAAuJIXAx0A4oncBQBneTF3zRcRAQAAAAAAxBlnYABwHa/elRkA4oXcBQBneTV3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7kxUAHgHgidwHAWV7M3aSbwLDzQ/L5fBH3Y6cPO7Zv327ZXl1dHZX1mMablmb+Ud9www3GmvT0dMv21157zdjHtGnTjDUm586di7imsLDQ2MfZs2eNNZMnTzbWuEmyBGWyjDMZkLsjI3fHXkPujk+y5FmyjDPWUlJSlJOTY1mTmppq2X7hwgXjelpbW4012dnZlu3d3d3GPkxjlaSuri7L9v/8n/+zsQ/TPpOktrY2y/Z9+/YZ+0hJMT/DIBrvV6Y+ovXeaqqx08fp06eNNT09PRGvZ3BwMOKaYDBo7MPO++/ixYst203HtCT9j//xP4w177zzjrEmEl7M3aSbwAAAE6/e1AgA4oXcBQBneTV3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7kxUAHgHgidwHAWV7MXeNdbNatW6f8/PxhNzo5ffq0ysrKNG/ePJWVlenMmTMxHSQAjNXl6wLt/ksk5C6AZETuAoCzkjl3x8s4gXH//fersbFx2Pdqamq0YsUKtbS0aMWKFaqpqYnZAAFgrMYa5okW6OQugGRD7gKAs5I9d8fLOIFxyy23aOrUqcO+19DQoMrKSklSZWWldu7cGZPBAcB4JXOgk7sAkhG5CwDOSubcHS/zg5BHcOLECfn9fkmS3+/XyZMnozooAMBw5C4AOIvcBYDEE/ObeNbW1qq2tlaS1NnZaZz58fl8Ea8zGrNLdvp4/PHHjTWPPfZYxGOJhlAoZKx5//33jTU5OTmW7f39/cY+br/9dsv2jIwMYx929Pb2Wrbb+SBy0003RWUsycR07CfL7G2yjDMWyN3HIh5LNDiVu8Fg0NgHuZvYyN3kd2Xu+nw+LV++3LLedD+NlBTz3xh7enqMNV1dXRGvx869P/bs2WPZfvbsWWMfHR0dxprTp09btpv2uyR997vfNdY48b5pR3p6urHGznuAyZIlS4w1Fy5csGw35b8kXbx40VjT3d1t2W4nZ9LSzL/iXnfddZbtpu2VpN///vfGmkmTJlm2Z2VlGfuw4sXcHdcZGDNmzFB7e7skqb29Xfn5+aPWVlVVqbm5Wc3NzZo+ffr4RgkAY+S2U+rIXQCJjty9lLumSUcAiBa35a4d45rAqKioUF1dnSSprq5Oq1atiuqgACBSbgt0chdAoiN3AcBZscjdxsZGzZ8/X4FAYMSbF4fDYT344IMKBAIqLi7WgQMHjMuO9lSno0ePKjs7WyUlJSopKVF1dbVxfMYJjHvuuUcf//jH9fbbb6ugoEBbt27Vhg0btGfPHs2bN0979uzRhg0bbO0MAHDCWMM80T5Ik7sAkg25CwDOikXuhkIhPfDAA9q9e7cOHz6sbdu26fDhw8Nqdu/erZaWFrW0tKi2tlbr1683Lmv1VKePfvSjOnjwoA4ePKgtW7YYx2i8QGjbtm0jfn/v3r3GzgEgXhLtw/FYkLsAkhG5CwDOinbu7t+/X4FAQIWFhZKktWvXqqGhQYsWLRqqaWho0H333Td0v5+uri61t7fr6NGjoy7b0NCgffv2Sbr0VKfbbrtNGzduHNcYx3UJCQAkunicUidJ+/btU0lJiYqKinTrrbdGc5MAIKGRuwDgrGjnbltbm2bPnj30dUFBgdra2mzVWC1r9VSnI0eO6IYbbtCtt96qF154wTjGmD+FBADiIdoz0pdPi9uzZ48KCgpUWlqqioqKYTPSXV1d+trXvqbGxkbNmTOHR+4B8BRyFwCcNdbc7ejo0LJly4a+rqqqUlVVlWV/H37qzmg1dpb9ML/fr/fff1/Tpk3TK6+8oi984Qt68803LZ/ewgQGAFeKxyl1Tz31lO6++27NmTNHkizvWA8AbkPuAoCzxpq7eXl5am5uHrW9oKBAra2tQ18fO3ZMM2fOtFUTDAZHXfbyU538fv+wpzplZmYqMzNT0qXHqH/0ox/VO++8M2yS5cMcncAIh8MaHByMqA+nnrts50ZN//qv/+rASJxz6tQpY43pmeO5ubnGPlpaWizbP/OZzxj7uOWWW4w1L774omW7nWefu02kr79kEYsbxI10WlxTU9OwmnfeeUf9/f267bbbdP78eX3jG9/QfffdF9VxjBW5m9jIXfcjd8cvWXO3s7NTW7dutawx5Wq09uXcuXMt29955x1jH7W1tcaa48ePW7Z3d3cb++jo6DDWHDt2zLL9lVdeMfaRnp5urDHt/4GBAWMf0XjtB4PBiPuw8x4eje2xM9aLFy8aa/r7+y3bo7U91157rWX7b3/7W2Mfdt7TTGN5/vnnjX1cOWl7pVjkbmlpqVpaWnTkyBHNmjVL9fX1euqpp4bVVFRUaNOmTVq7dq2ampo0efJk+f1+5eXljbrs5ac6bdiwYdhTnTo6OjR16lSlpqbqz3/+s1paWoYmrUfDGRgAXCkep9QNDAzolVde0d69e9XT06OPf/zjWr58ua6//voxjh4Akg+5CwDOivYERlpamjZt2qSVK1cqFApp3bp1KioqGno6SHV1tcrLy7Vr1y4FAgHl5OToiSeesFxWuvRHqjVr1mjr1q2aM2eOtm/fLunSBM7f/d3fKS0tTampqdqyZYumTp1qPcaobjEAJIh4nVI3ffp0TZgwQRMmTNAtt9yiQ4cO8UEagCeQuwDgrFg8/am8vFzl5eXDvlddXT30f5/Pp82bN9teVpKmTZs24lOdVq9erdWrV49pfN47lxOAJ0T7rsxXnlIXDAZVX1+vioqKYTWrVq3SCy+8oIGBAV28eFFNTU1auHBhrDYRABIKuQsAzorF058SHWdgAIANdk6pW7hwoe68804VFxcrJSVFX/nKV7R48eI4jxwAkhO5CwD4MCYwALhSPE6pk6SHH35YDz/8cNTXDQCJjtwFAGe55ayKsWACA4DruOk0OQBIBuQuADjLq7nLBAYAV/JioANAPJG7AOAsL+YuExgAXMmLgQ4A8UTuAoCzvJi7jk9gDA4ORrR8Sor5wSkffkb4SP7pn/7Jsv1nP/uZ3SG5RjAYjLiPvr4+Y80zzzxj2f6pT33K2EdPT4+xpqOjw7J95cqVxj68yC1B6JbtiAZyN3GRu5Dck1du2Y5IDQ4O6uLFi5Y1vb29lu0DAwPG9WRmZhprTp48adk+f/58Yx+nTp0y1pjeZ86ePWvsw5QfdtazYMECYx9vv/12xOtJT0839mHKTDvvzXbeW001EyZMMPZx/vx5Y43p/cp0TEv23q9M+3bSpEnGPvr7+401ptfP+++/b+zDzn47c+aMZXt+fr6xDytezF3OwADgOl69JhAA4oXcBQBneTV3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7kxUAHgHgidwHAWV7MXSYwALiSFwMdAOKJ3AUAZ3kxd5nAAOA6Xr2pEQDEC7kLAM7yau4ygQHAlbwY6AAQT+QuADjLi7nLBAYAV/JioANAPJG7AOAsL+auoxMYg4OD6uvrc3KVo3rttdcs230+n0MjSS7BYDDiPnp7ey3b6+rqIl6HHdXV1Y6sB4gncjf5kbtA8hkcHLRsz8jIsGxPT083ruM//af/ZKzZs2ePZbud94dXX33VWGMa77Fjx4x9vP7668aaW2+91bK9q6vL2Ed3d7exZtKkSZbtOTk5xj78fr9le2trq7GPtDTzr2oDAwOW7bm5ucY+7LzPmI6Vnp4eYx8TJkww1syaNcuyfeHChcY+fve73xlrTNtj5/V18OBBY41p/5uyQJJCoZCxxks4AwOAK3lxRhoA4oncBQBneTF3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgOt49a7MABAv5C4AOMurucsEBgBX8mKgA0A8kbsA4Cwv5i4TGABcyYuBDgDxRO4CgLO8mLtMYABwJS8GOgDEE7kLAM7yYu4ygQHAlbwY6AAQT+QuADjLi7nr2QmMtWvXWrY3Nzcb++jo6IjWcDwlGi80n88XlRqv8UrIefWmRonOlLsvv/yysY/Ozs5oDcdTnHo9pKSkOLKeZOKVLCJ3/8Pg4KCCwaBlTXZ2tmX7xYsXba3HZGBgwLK9tbXV2Me0adOMNVlZWRGvZ9GiRcaaXbt2WbZfd911xj5SU1ONNRMmTLBst5N1pn3S19dn7MPOekyvu+7ubmMfOTk5xppQKBRRu2Rv35s+v2dmZhr7MB33kpSWZv1r8Mc//nFjH3Z+FzQdB6ZxSNLx48dH/L5Xc9ezExgA3M2LgQ4A8UTuAoCzvJi7TGAAcCUvBjoAxBO5CwDO8mLuMoEBwJW8GOgAEE/kLgA4y4u5ywQGAFfyYqADQDyRuwDgLC/mrvHOMOvWrVN+fr4WL1489L3HHntMs2bNUklJiUpKSow31AEA2EfuAoCzyF0AuKSxsVHz589XIBBQTU3NVe3hcFgPPvigAoGAiouLdeDAAeOyp0+fVllZmebNm6eysjKdOXNmWJ/vv/++cnNz9c///M/G8RknMO6//341NjZe9f2HHnpIBw8e1MGDB1VeXm5cEQA45fJdmcfyzw5ToO/bt0+TJ08e+rD7D//wD+MaP7kLINmQuwDgrFjkbigU0gMPPKDdu3fr8OHD2rZtmw4fPjysZvfu3WppaVFLS4tqa2u1fv1647I1NTVasWKFWlpatGLFiqvy/KGHHtJnPvMZW9ttvITklltu0dGjR211BgCJItqn1F0O5T179qigoEClpaWqqKi46rFvf/EXf6H/+3//b0TrIncBJCNyFwCcFe3c3b9/vwKBgAoLCyVJa9euVUNDw7DcbWho0H333Sefz6fly5erq6tL7e3tOnr06KjLNjQ0aN++fZKkyspK3Xbbbdq4caMkaefOnSosLDQ+vviycT+wfdOmTSouLta6deuuOgUEAOIt2jPSVwZ6RkbGUCg7idwFkMjIXQBwVrRzt62tTbNnzx76uqCgQG1tbbZqrJY9ceKE/H6/JMnv9+vkyZOSpO7ubm3cuFGPPvqo7W0e1008169fr0ceeUQ+n0+PPPKIvvWtb+mnP/3piLW1tbWqra2VdOnal0Rxxx13WLZ//etfN/Yxlh2N/zB58mTL9pSUcc+rDXPjjTdatodCoaisJ5nYCS633Awo2tsxUig3NTVdVffv//7vWrp0qWbOnKl//ud/VlFRUVTW74Xc/W//7b8Z+yB3x4fcjR9yd/ySNXdTUlI0ffp0y74vXrxo2T5t2jTj+P7f//t/xpq5c+datr/77rvGPi7/omFl2bJllu2mbJCkf/mXfzHWDA4OWrb39PQY+7Bj4sSJlu19fX3GPt5//33Ldju5a+c1FQwGLdvt5G5mZmbEY7GzPb29vcaajo4Oy/bLf8W3kp6ebqwZGBiwbF+wYIGxj+effz7isaSlmX8dP378+KhtY83djo6OYa/XqqoqVVVVWfbn8/mM6/T5fLaW/bBHH31UDz30kHJzc41jv2xcExgzZswY+v9Xv/pVfe5znxu19sqdUlxcPJ7VAcCYxSPQb7zxRr333nvKzc3Vrl279IUvfEEtLS1jHPnIyF0AiY7cvTR2O788AUA0jDV38/Ly1NzcPGp7QUGBWltbh74+duyYZs6caasmGAyOuuyMGTPU3t4uv9+v9vZ25efnS5Kampr0q1/9Sn/zN3+jrq4upaSkKCsry/JkgnFNYFxeuSTt2LFj2B2bASDexnKDuMuiEeiTJk0a+n95ebm+9rWvqbOz0/iXODvIXQCJjNwFAGeNJ3dNSktL1dLSoiNHjmjWrFmqr6/XU089NaymoqJCmzZt0tq1a9XU1KTJkyfL7/crLy9v1GUrKipUV1enDRs2qK6uTqtWrZIkvfDCC0P9PvbYY8rNzTVeCWGcwLjnnnu0b98+dXZ2qqCgQH//93+vffv26eDBg/L5fJo7d65+/OMfj3nnAEAsxSPQjx8/rhkzZsjn82n//v0aHBy0dRrwh5G7AJIRuQsAzop27qalpWnTpk1auXKlQqGQ1q1bp6KiIm3ZskWSVF1drfLycu3atUuBQEA5OTl64oknLJeVpA0bNmjNmjXaunWr5syZo+3bt49/jKaCbdu2XfW9L3/5y+NeIQA4IR6B/qtf/Uo/+tGPlJaWpuzsbNXX1xuv/RsJuQsgGZG7AOCsWNxDqby8/KrHRldXVw/93+fzafPmzbaXlS7d12fv3r2W633sscdsjW9cl5AAQKKLR6B//etft3UDYABwI3IXAJzllptAjwUTGABcyYuBDgDxRO4CgLO8mLtMYABwnVjc1AgAMDpyFwCc5dXcZQJjFFfe1dor7LwATNeV2nn2uukZ3SNdN/VhpueaS9KZM2cs2zs7O419jOWZxMnASyHnpW11C3J3ZNHI3YsXL1q2k7ux46Us8tK2WklLS9OUKVMsa7Kysizbu7u7jetJSUkx1rS3t1u22/mZdXV1GWuufJLASOzklJ3X/rlz5yzb7TzCNiMjw1gzceJEy/a3337b2Icpu+3s+76+PmON6TgYHBw09nHhwgVjTWpqqmV7Wpr510o723P27NmI1xMIBIw1puPNzs2ABwYGjDWmY9K0X028mLvm5AMAAAAAAIgzzsAA4EpenJEGgHgidwHAWV7MXSYwALiSFwMdAOKJ3AUAZ3kxd5nAAOBKXgx0AIgnchcAnOXF3GUCA4DrePWuzAAQL+QuADjLq7nLBAYAV/JioANAPJG7AOAsL+YuExgAXMmLgQ4A8UTuAoCzvJi7TGAAcCUvBjoAxBO5CwDO8mLuunICIxo/yGeeeSYKI3GfoqIiy/aqqipjH62trZbtCxYsMPaRmppqrDEdB0ePHjX2MXfuXGNNMkmmkIt0rMm0rbiE3B0ZuZvckimLyN3oSElJ0aRJkyxrSkpKLNv37t1rXM/p06eNNZmZmZbtPp/P2Ed+fr6xpqenx7K9qanJ2EdOTo6xxu/3W7bPnj3b2EdWVpax5pVXXjHWRCo3N9dYEwwGjTUDAwOW7Xa2d3Bw0FiTnZ1t2W7nWEpLM//qacqR+fPnG/v49Kc/baz513/9V8v2wsJCYx923tOuvfZay/b9+/cb+7j++utHbfNi7rpyAgOAt3n1pkYAEC/kLgA4y6u5ywQGAFfyYqADQDyRuwDgLC/mLhMYAFzJi4EOAPFE7gKAs7yYu0xgAHAlLwY6AMQTuQsAzvJi7qbEewAAAAAAAAAmnIEBwJW8OCMNAPFE7gKAs7yYu0xgAHAdr96VGQDihdwFAGd5NXeZwADgSl4MdACIJ3IXAJzlxdx1dAKjt7dXb731lmVNVlaWZXt+fr5xPdnZ2caaN99807L9N7/5jbEPL7r//vst29PT0419zJo1y7I9NTV1LEMat46ODmPNqVOnjDXTpk2LeCyDg4PGmpQU61vWhEIhYx92Qs6pIDStx872RNK/VziVuzk5OcaaN954w7Kd3B0ZuXs1cnd8yF1nXLx4Uc3NzZY1pnY7r8kJEyYYay5cuGDZbue1ZOe4mDRpkmX7NddcY+zDTpZNmTLFsj0QCBj7KC8vN9Y8/PDDlu0+n8/Yh4mdDLKzT0x6e3uNNRkZGcYa0+s7Lc38a2V3d7exJhgMWrZ3dnYa+zB9rpHMx1JTU5OxDzvvaXPnzrVsv/HGG419WL2OvZi7nIEBwJW8GOgAEE/kLgA4y4u5ywQGAFfyYqADQDyRuwDgLC/mLhMYAFzHqzc1AoB4IXcBwFlezV0mMAC4khcDHQDiidwFAGd5MXeZwADgSl4MdACIJ3IXAJzlxdy1vr02ACSpy6fV2f1nR2Njo+bPn69AIKCamppR615++WWlpqbqV7/6VbQ2BwASHrkLAM6KRe4mOiYwALhStAM9FArpgQce0O7du3X48GFt27ZNhw8fHrHu29/+tlauXBmLzQKAhEXuAoCz4jFxHA6H9eCDDyoQCKi4uFgHDhwwLnv69GmVlZVp3rx5Kisr05kzZyRJ+/fvV0lJiUpKSrR06VLt2LHDOD5HLyHp7+9Xe3u7ZY3pGcJ2nrdr53nVP/7xj401XpOZmWmsMT0vvKenx9iHnWebm9h5LrbpmeQXL1409mHnGdA333yzZXtOTo6xDzv7pK+vz7Ldzs/PqZlX0zO8JenUqVOW7QMDA+Nefyxmmffv369AIKDCwkJJ0tq1a9XQ0KBFixYNq/vf//t/a/Xq1Xr55Zejuv7xIncTG7l7Nbflrh0+ny/iPkxjlchdJ6WkWP+NMJJ9fZmd11M02Ml30/ZOmDDB2MfEiRONNaaM6erqMvbxgx/8wFhjYuc1a3pvtbO9nZ2dxpqsrCzLdjvHmp3X7tmzZy3b7ewTOzXXX3+9ZfuaNWuMfdTX1xtrTMdsd3e3sY+ZM2caa/bt22fZ3tvba+yjpKRkxO/HIncvTxzv2bNHBQUFKi0tVUVFxbDc3b17t1paWtTS0qKmpiatX79eTU1NlsvW1NRoxYoV2rBhg2pqalRTU6ONGzdq8eLFam5uVlpamtrb27V06VJ9/vOft3z9cAYGAFeK9ox0W1ubZs+ePfR1QUGB2trarqrZsWOHqquro749AJDoyF0AcFa0c/fKieOMjIyhieMrNTQ06L777pPP59Py5cvV1dWl9vZ2y2UbGhpUWVkpSaqsrNTOnTslXfqDw+XJit7eXluTXExgAIAunWWwbNmyoX+1tbXD2kcK/Q+H7H//7/9dGzdujMpfuwHA7chdAEgsdieOR6qxWvbEiRPy+/2SJL/fr5MnTw7VNTU1qaioSEuWLNGWLVuMZy/xFBIArjTWU+ry8vLU3Nw8antBQYFaW1uHvj527NhVpw42Nzdr7dq1ki6d9rlr1y6lpaXpC1/4wpjGAgDJiNwFAGeNNXcvTxxfVlVVpaqqKsv+PjxxPFqNnWVHcvPNN+vNN9/UH//4R1VWVuozn/mM5aVRTGAAcKVoXxNYWlqqlpYWHTlyRLNmzVJ9fb2eeuqpYTVHjhwZ+v/999+vz33uc3yIBuAZ5C4AOCseE8ej1QSDwVGXnTFjhtrb2+X3+9Xe3q78/Pyr1r1w4UJNmDBBb7zxxrBJlg8zXkLS2tqq22+/XQsXLlRRUZF++MMfShr9TqIAkAiifU1gWlqaNm3apJUrV2rhwoVas2aNioqKtGXLFm3ZsiWqYyd3ASQjchcAnBXt3L1y4jgYDKq+vl4VFRXDaioqKvTkk08qHA7rpZde0uTJk+X3+y2XraioUF1dnSSprq5Oq1atknRpEvryjWbfe+89vf3225o7d67lGI1nYKSlpekHP/iBbrzxRp0/f1433XSTysrK9LOf/WzEO4kCQLzF4q7MklReXq7y8vJh3xvtxnE/+9nPxr0echdAsiF3AcBZscjdKyeOQ6GQ1q1bNzRxLF3K3/Lycu3atUuBQEA5OTl64oknLJeVpA0bNmjNmjXaunWr5syZo+3bt0uSXnzxRdXU1Cg9PV0pKSn6l3/5F02fPt16jKaN8Pv9QzfcmDhxohYuXKi2tjY1NDQMPRamsrJSt912G4EOIGE49cjYWCB3ASQjchcAnBWPiWOfz6fNmzfbXla69Ej4vXv3XvX9L33pS/rSl740pvGN6R4YR48e1auvvqqbb77Z8k6iABBvyfxB+krkLoBkQe4CgLPckrtjYXsC48KFC1q9erUef/xxTZo0yfYKamtrhx6Lde7cOQWDQcv63t5ey/bOzk7jOt99911jzbPPPmus8ZqRbqbyYadOnbJs7+7ujngcU6dOjbiPaLnyRjSR1DghEAgYaz75yU9GvJ4rb5g2GjuvwbfeeivisVhxQ6CTu+6Xl5dnrCF3x1fjhHnz5hlropG7f/7zn4015G50RCN3fT6f0tPTLetTUqxvQ3f5mvBI3XPPPZbtTz/9tLEPO2MJhUKW7abtlWT51IHLTFm1aNEiYx//9m//ZqwxPcZxcHAw4j4WL15s7MPOExxM63nxxReNfUybNi3isdiZ3CssLDTWTJkyxbK9pqbG2EdBQYGx5oMPPrBsf++994x92Mk802cw02snGmNwG3OaSOrv79fq1at177336u6775b0H3cSlTTqnUSlS49maW5uVnNzsyZOnBilYQOAtWjf1Mhp5C6AZEPuNlve3R8Aoi3Zc3c8jBMY4XBYX/7yl7Vw4UJ985vfHPr+aHcSBYB4G2uYJ1qgk7sAkg25CwDOSvbcHS/jJSR/+MMf9POf/1xLlixRSUmJJOm73/3uqHcSBYBEkMwhTe4CSEbkLgA4K5lzd7yMExif+tSnRt0xI91JFAASQTIHOrkLIBmRuwDgrGTO3fEa01NIACBZeDHQASCeyF0AcJYXc9fWTTwBAAAAAADiiTMwALiSF2ekASCeyF0AcJYXc5cJDACu46Y7LQNAMiB3AcBZXs1dRycwUlJSlJOTE1Ef2dnZttZjMjg4GNE43KiqqspYY/r5dXd3R2s4ERsYGIioPdH09/dbtre0tBj7CIVCEY+jr6/PWHP+/PmI1xMpLwb6SMjdxPZXf/VXxhpyN35MufvOO+8Y+4jGNpO7ySU9PV0zZsywrAkGg5btJ0+eNK7HTqY+99xzlu12jpvc3Fxjjek4t3Ns2NmeCxcuWLa/+uqrxj7svF+lpVn/imTn81Rvb29E7ZK9/dbZ2WnZbsoxSfL5fMaa6dOnW7bb+awRCASMNcePH7dsnzVrlrGPhQsXGmtuuOEGy/YdO3YY+8jIyDDWmPZtZmamsQ8rXsxdzsAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7kxUAHgHgidwHAWV7MXSYwALiOV29qBADxQu4CgLO8mrtMYABwJS8GOgDEE7kLAM7yYu4ygQHAlbwY6AAQT+QuADjLi7nLBAYAV/JioANAPJG7AOAsL+auoxMYPp9PqampEfVhZ/mDBw9GtA6vysrKMtZ0d3dHvJ6ZM2datpuevW2XqR8722t6rrkkdXR0WLbn5eUZ+5gwYYKxxs4zup1g53nVdmpMzxO34/nnnx/x+169JnAk5G5iI3ev5lTu5ubmGmsSBbmbXDIyMlRQUGBZk5KSYtk+f/5843qeffZZY83p06ct20OhkLGPYDBorDG9bu28ru2sp6enJ6J26dLPx8SUZXaO9f7+fsv2kpISYx85OTnGGtOxZGe/Pvnkk8aahQsXWra3t7cb+9ixY4expqyszLL94sWLxj4OHTpkrCksLLRsnzFjhrGPEydOGGtMx5Kdn89ovJq7nIEBwJW8GOgAEE/kLgA4y4u5az1lBwAAAAAAkAA4AwOAK3lxRhoA4oncBQBneTF3mcAA4EpeDHQAiCdyFwCc5cXcZQIDgCt5MdABIJ7IXQBwlhdzlwkMAK7j1bsyA0C8kLsA4Cyv5i438QTgSpdD3e4/OxobGzV//nwFAgHV1NRc1d7Q0KDi4mKVlJRo2bJlevHFF6O9WQCQsMhdAHBWLHI30XEGBgBXinZIh0IhPfDAA9qzZ48KCgpUWlqqiooKLVq0aKhmxYoVqqiokM/n02uvvaY1a9borbfeiuo4ACBRkbsA4Cy3TEqMhSsnMM6cOWOsSUmxPvkkFApFazgJYe7cuVHpZ3Bw0LLdzosoNzc3KmNxwtSpU401p0+ftmw37TPERrQDff/+/QoEAiosLJQkrV27Vg0NDcM+SF95bHd3d8vn80V1DImM3L0auTs+5G7yIncvCYfDGhgYiPl6TJkqybg/7OyvYDBorOnv74+o3W6NaSzPP/+8sY++vj5jjelYzszMNPZx/fXXW7afPHnS2Ied7M7IyLBsnzhxorGP3t5eY017e7tl+6FDh4x91NfXG2suXLhg2Z6VlWXswzRWScbXqGkcknOvQStenMDgEhIArhTtU+ra2to0e/bsoa8LCgrU1tZ2Vd2OHTu0YMECffazn9VPf/rTqG4TACQychcAnBWPS/fC4bAefPBBBQIBFRcX68CBA8ZlT58+rbKyMs2bN09lZWVDf/jas2ePbrrpJi1ZskQ33XSTnnnmGeP4mMAA4DpjDfNwOKyOjg4tW7Zs6F9tbe1VfX7YSLPmd911l9566y3t3LlTjzzySMy2EQASCbkLAM4aT+6aXL50b/fu3Tp8+LC2bdumw4cPD6vZvXu3Wlpa1NLSotraWq1fv964bE1NjVasWKGWlhatWLFiaHJj+vTp+s1vfqPXX39ddXV1+tKXvmQcoysvIQGAsZ5Sl5eXp+bm5lHbCwoK1NraOvT1sWPHNHPmzFHrb7nlFr377rvq7OzU9OnTxzQWAEhG5C4AOCsel+41NDTovvvuk8/n0/Lly9XV1aX29nYdPXp01GUbGhq0b98+SVJlZaVuu+02bdy4UTfccMNQv0VFRert7VVfX5/lZVqcgQHAlaI9I11aWqqWlhYdOXJEwWBQ9fX1qqioGFbzpz/9aaivAwcOKBgMatq0aTHZPgBINOQuADgr2me+2bl0b7Qaq2VPnDghv98vSfL7/SPeA+bXv/61brjhBuM9ZjgDA4ArRXtGOi0tTZs2bdLKlSsVCoW0bt06FRUVacuWLZKk6upq/frXv9aTTz6p9PR0ZWdn65e//GVC3FAOAJxA7gKAs6J95pudS/dGq7F72d9I3nzzTX3729/W008/baxlAgMAbCovL1d5efmw71VXVw/9/9vf/ra+/e1vOz0sAHAtchcAnGPn0r3RaoLB4KjLzpgxQ+3t7fL7/Wpvb1d+fv6wurvuuktPPvmkPvrRjxrHyCUkAFwpFndlBgCMjtwFAGfF49K9iooKPfnkkwqHw3rppZc0efJk+f1+y2UrKipUV1cnSaqrq9OqVaskSV1dXfrsZz+r733ve/rkJz9pa5s5AwOA6/DhGACcRe4CgLNikbt2Lt0rLy/Xrl27FAgElJOToyeeeMJyWUnasGGD1qxZo61bt2rOnDnavn27JGnTpk3605/+pO985zv6zne+I0l6+umnh52hcdUYo7rFCeLmm2821vz5z3+2bD937ly0hpMQli1bZqyxc43S4OCgZbvppiuSlJKSPCf+pKWZXyJ5eXmW7b29vdEaDsaAD9LOInevRu6OD7mbvMjdSwYHB43HYGpqqmX7hAkTjOtJT0831phe+3Zeb6YMkqT+/n7L9lAoZOxjYGDAWGM6xnp6eiLuQzJns53tufKpDSPp7u429tHX12esCQaDlu0XLlww9pGdnW2s+eCDDyzbP/KRjxj72Lp1q7HmpZdesmy389qw85535swZy/asrCxjH5MmTTLWmI4V02vHJBa5a7p0z+fzafPmzbaXlaRp06Zp7969V33/b//2b/W3f/u3YxqfKycwAIAP0gDgLHIXAJzlxdxlAgOAK3kx0AEgnshdAHCWF3OXCQwAruTFQAeAeCJ3AcBZXsxd4wVCra2tuv3227Vw4UIVFRXphz/8oSTpscce06xZs1RSUqKSkhLt2rUr5oMFADvGekfmRAt/chdAsiF3AcBZyZ6742U8AyMtLU0/+MEPdOONN+r8+fO66aabVFZWJkl66KGH9Nd//dcxHyQAjFUyhzS5CyAZkbsA4Kxkzt3xMk5g+P1++f1+SdLEiRO1cOFCtbW1xXxgABCJZA50chdAMiJ3AcBZyZy74zWm56odPXpUr7766tDj8jZt2qTi4mKtW7fO+CgaAHCSW06pI3cBJAtyFwCc5ZbcHQvbN/G8cOGCVq9erccff1yTJk3S+vXr9cgjj8jn8+mRRx7Rt771Lf30pz+9arna2lrV1tZKks6dOxe9kVu44447jDVvvvmmZXtTU5Oxj0if2+ukG264wZH15OfnO7KeRDJ58mTL9szMTIdGgiu5IaTJ3auRu1cjd69mJ3ftZITP57M9JpC7l3O3v7/fmL2Dg4OW7XbO/rh8xoiVo0ePWrZPmjTJ2EdWVpaxxvRasfNa6uvrM9aY9mtqaqqxj4GBAWNNSor133jtZExLS4uxxsTO9pj2rWlbJCkQCBhr/vjHP1q2d3d3G/vIzc011ixevNiy/eTJk8Y+zp49a6xJS7P+NXjixInGPuz8fILBoGV7pJ9r3JC7Y2XrDIz+/n6tXr1a9957r+6++25J0owZM5SamqqUlBR99atf1f79+0dctqqqSs3NzWpubrYVkgAQKTfc1IjcBZBMyN3/yF07vzACQKTckLvjYUzYcDisL3/5y1q4cKG++c1vDn2/vb196P87duwwzpYBgJOSOdDJXQDJiNwFAGclc+6Ol/ESkj/84Q/6+c9/riVLlqikpESS9N3vflfbtm3TwYMH5fP5NHfuXP34xz+O9VgBwBPIXQBwFrkLAMnBOIHxqU99asTZmvLy8pgMCACiIZlnmcldAMmI3AUAZyVz7o6X7Zt4AkAy8WKgA0A8kbsA4Cwv5i4TGABcaayBztMGACAy5C4AOMuLucsEBgDXGc+NitwQ6AAQL+QuADjLq7nLBAYAV/LiKXUAEE/kLgA4y4u569kJjKysLMv2SZMmGfvo6uoy1oRCIbtDiohpvOnp6cY+BgcHjTU9PT2W7W6Y1RurtLTIX0b9/f3GmoyMjIjX4yVeDPRE57bcnThxomU7uRs75G5iIncvCYVCOnPmjLHGip3j8+zZs2Ma10guXLhgrDl//ryxxvSzX7BggbGPvr4+Y81rr71m2Z6SkmLsIzs721hjymY72X348GHLdjv5kpOTE3GNnfdWOzXTp0+3bD9x4oSxj+PHjxtrbrjhBmONiZ3j4OTJk5btdvIsNTXVWGP67BPp+4wXc9ezExgA3M2LgQ4A8UTuAoCzvJi7TGAAcCUvBjoAxBO5CwDO8mLuMoEBwHXGc1MjAMD4kbsA4Cyv5i4TGABcyYuBDgDxRO4CgLO8mLtMYABwJS8GOgDEE7kLAM7yYu4ygQHAlbwY6AAQT+QuADjLi7lrfsYMAECS1NjYqPnz5ysQCKimpuaq9l/84hcqLi5WcXGxPvGJT+jQoUNxGCUAAADgTkxgAHClyzc2svvPJBQK6YEHHtDu3bt1+PBhbdu27apnvF933XV67rnn9Nprr+mRRx5RVVVVrDYPABJOtHNXYuIYAKzEIncTnWcvIent7bVsnzZtmrGP7OxsY01nZ2dE45Ck9PR0Y831119vrDEZGBgw1vT09Fi2p6QwJ/ZhGRkZxho7x4GdfnBJLEJ6//79CgQCKiwslCStXbtWDQ0NWrRo0VDNJz7xiaH/L1++XMeOHYvqGJKd23J3/vz5xhoTO7lrGi+5ezVy13mxyN3LE8d79uxRQUGBSktLVVFRMSx3L08cX3PNNdq9e7eqqqrU1NQU1XGMVSgUUldXl2VNamqqZXtmZqZxPXbyMBr5Yed1MH36dMv2rKwsYx/vvfeescaUmdddd52xj/b2dmON6ecTCoWMfXzkIx+xbD937pyxD9P7mWTet3Zelzk5Ocaaz33uc5btv/71r6OyHtN+uXDhgrEP089PMh/7fX19xj7svE5NNXY+b4zGTZMSY8GnHgCuFO0Z6ba2Ns2ePXvo64KCArW1tY1av3XrVn3mM5+JyrYAQDKIdu5eOXGckZExNHF8pU984hO65pprJDFxDMB7vHgGBhMYAFxprIHe0dGhZcuWDf2rra29qr8P8/l8I6772Wef1datW7Vx48aYbBsAJKJo5y4TxwBgLR6X7oXDYT344IMKBAIqLi7WgQMHjMuePn1aZWVlmjdvnsrKynTmzBlJ0qlTp3T77bcrNzdXX//6122Nz7OXkABwt7HOMufl5am5uXnU9oKCArW2tg59fezYMc2cOfOqutdee01f+cpXtHv3bluXRACAW0Q7d8czcfziiy+OaQwAkMzicene7t271dLSopaWFjU1NWn9+vVqamqyXLampkYrVqzQhg0bVFNTo5qaGm3cuFFZWVn6zne+ozfeeENvvPGGrTFyBgYAV4r2jHRpaalaWlp05MgRBYNB1dfXq6KiYljN+++/r7vvvls///nPo3JfGgBIJtHO3bFOHDc0NDBxDMBT4nHpXkNDg+677z75fD4tX75cXV1dam9vt1y2oaFBlZWVkqTKykrt3LlTkjRhwgR96lOfsnWfnMs4AwOA68TiOr+0tDRt2rRJK1euVCgU0rp161RUVKQtW7ZIkqqrq/UP//APOnXqlL72ta8NLWP110UAcItY5O6VE8ezZs1SfX29nnrqqWE1TBwD8KpY5O5Il+59+MbIo13eZ7XsiRMn5Pf7JUl+v18nT54c9xiZwADgSrG4UVF5ebnKy8uHfa+6unro/z/5yU/0k5/8JOrrBYBkwMQxADhrrLl7+d5Dl1VVVamqqsqyvw9fujdazVgu+4sEExgAXMktd1oGgGTBxDEAOCse93wbrSYYDI667IwZM9Te3i6/36/29nbl5+ePadxX8uwExuDgoGX7xIkTjX3k5uYaa6ZMmWLZbuf0GTvXBH3iE5+IaBx2ReM5315j5/npdo4lEzsB5qVf6r20rcnCqdydPHmyZXtHR4exD3I3uZG78eGlbbWSlpZmfP2npqZatts5hqNhYGDAWBMKhYw1H3zwgWV7e3u7sY/+/n5jjWm/vffee8Y+7BynpvXceeedxj5M2X3lL3qR1JgeHWznPc/Oz/jdd9+1bLfzHv7Rj37UWGMa7+nTp4195OXlGWsmTJhg2d7b22vsY9KkScaanJyciMZhGks8Lt2rqKjQpk2btHbtWjU1NWny5Mny+/3Ky8sbddmKigrV1dVpw4YNqqur06pVq8Y9Rs9OYABwNz5IA4CzyF0AcFY8Lt0rLy/Xrl27FAgElJOToyeeeMJyWUnasGGD1qxZo61bt2rOnDnavn370Drnzp2rc+fOKRgMaufOnXr66aeHPfXkqjFGdYsBIAHE4qZGAIDRkbsA4KxY5a7p0j2fz6fNmzfbXlaSpk2bpr179464zNGjR8c0Ph6jCgAAAAAAEh5nYABwJf4SCADOIncBwFlezF0mMAC4khcDHQDiidwFAGd5MXeZwADgSl4MdACIJ3IXAJzlxdxlAgOAK3kx0AEgnshdAHCWF3OXCQwArsPd8AHAWeQuADjLq7nr2QkMn89n2Z6VlRWV9eTk5ESlH5OpU6c6sp7Jkyc7sh6MnemYlrw1S+ulbU0WiZK7dl4rdpC7IHeH89K2WklNTTXmQ1qa9UdwO3mYnp5urGltbTXWREMoFLJsT0kxP/jQTk001jM4OGiseeeddyJql6Q//vGPlu0ZGRnGPkzbK0m9vb2W7e+9956xDzvHW25urmX7lClTjH0cOnTIWGPanmi9NjIzMy3be3p6jH2kpqYaa0z77ZprrjH28cEHH4za5sXc9ewEBgB382KgA0A8kbsA4Cwv5i4TGABcyYuBDgDxRO4CgLO8mLtMYABwJS8GOgDEE7kLAM7yYu4ygQHAdbx6UyMAiBdyFwCc5dXcNd7ppre3Vx/72Me0dOlSFRUV6dFHH5UknT59WmVlZZo3b57Kysp05syZmA8WAOy6HOp2/yUSchdAMiJ3AcBZyZy742WcwMjMzNQzzzyjQ4cO6eDBg2psbNRLL72kmpoarVixQi0tLVqxYoVqamqcGC8A2JLMgU7uAkhG5C4AOCuZc3e8jBMYPp9v6PEv/f396u/vl8/nU0NDgyorKyVJlZWV2rlzZ0wHCgBeQe4CgLPIXQBIDuaHJevSM4hLSkqUn5+vsrIy3XzzzTpx4oT8fr8kye/36+TJkzEdKACMRbLPSJO7AJINuQsAzkr23B0PWzfxTE1N1cGDB9XV1aW77rpLb7zxhu0V1NbWqra2VpJ07ty58Y1yjJz64eTk5DiyHjtCoVDEfWRkZERhJEBiSPaQJndHRu4CiYvcvZS74XBYkyZNsqxPS7P+CD558mTjOt9++21jTUqK9d8q7fzMiouLjTV/+tOfLNu7u7uNfUSDab9K0qJFi4w16enplu1Tp0419nHttddatp89e9bYx6lTp4w1HR0dlu1ZWVnGPgYHB4012dnZlu1Tpkwx9vHqq68aa0zH/pIlS4x9tLW1GWt8Pp9lu53Xhp3jbcKECZbtl8/8Gq9kz93xsHUGxmVTpkzRbbfdpsbGRs2YMUPt7e2SpPb2duXn54+4TFVVlZqbm9Xc3GwMcwCIhrHORidy+JO7AJIBufsfuWvnlxoAiJSbcncsjBMYHR0d6urqkiT19PTo97//vRYsWKCKigrV1dVJkurq6rRq1aqYDhQAxiKZA53cBZCMyF0AcFYy5+54GaeI29vbVVlZqVAopMHBQa1Zs0af+9zn9PGPf1xr1qzR1q1bNWfOHG3fvt2J8QKALckc0uQugGRE7gKAs5I5d8fLOIFRXFw84vVK06ZN0969e2MyKACIVDIHOrkLIBmRuwDgrGTO3fHiIj0AruTFQAeAeCJ3AcBZXsxdJjAAuI6brvMDgGRA7gKAs7yau0xgAHAlLwY6AMQTuQsAzvJi7rpyAsP0XF8n12M6qEzP57ZrYGAg4j4yMjKiMBIgMcQi0BsbG/WNb3xDoVBIX/nKV7Rhw4Zh7W+99Zb+63/9rzpw4ID+5//8n/rrv/7rqI8hUTmVu9FA7gKxQe5e4vP5jI9SNeXQmTNnjOuZN2+eseb48eOW7XYy6OjRo8aanp4eY40T7OT75z//+Yj7yc7ONvYxYcKEiNrtrsf0M7TzWF87r91QKGTZfuLECWMfU6ZMMdaYHv9+6tQpYx92mD632Nkndj77mPZ/VlaWsQ8rTGAAgEtEO9BDoZAeeOAB7dmzRwUFBSotLVVFRYUWLVo0VDN16lT9r//1v7Rz586orhsAkgG5CwDO8uIERnT+DAUACSbaz8Xev3+/AoGACgsLlZGRobVr16qhoWFYTX5+vkpLS5Wenh6rzQKAhEXuAoCzop27yYAJDACuM9YwD4fD6ujo0LJly4b+1dbWDuuzra1Ns2fPHvq6oKBAbW1tTm8aACSk8eSuCbkLAKOLRe4mAy4hAQBJeXl5am5uHrV9pNBPpvs+AECiuTxxfFlVVZWqqqqGviZ3AQAfxgQGAFeK9ixzQUGBWltbh74+duyYZs6cGdV1AEAyG2vumiaOyV0AsOaWsyrGgktIALhStE+pKy0tVUtLi44cOaJgMKj6+npVVFQ4sCUAkBzIXQBwFpeQAIBLRDuk09LStGnTJq1cuVKhUEjr1q1TUVGRtmzZIkmqrq7W8ePHtWzZMp07d04pKSl6/PHHdfjwYeMjwQDADchdAHCWWyYlxoIJDACuFItALy8vV3l5+bDvVVdXD/3/2muv1bFjx6K+XgBIBuQuADjLixMYvrCDW52bm6sFCxY4tbqIdHR0KC8vL97DsCWZxiol13gZa2xEa6xHjx5VZ2fnVd+/8847R/y+lenTp6uxsTHiMSUacjc2kmmsUnKNl7HGBrnrHHI3NpJprFJyjZexxga5GxuOTmAsW7bM8mZNiYSxxk4yjZexxkYyjTXZJdO+Zqyxk0zjZayxkUxjTXbJtK8Za+wk03gZa2wk01iTCTfxBAAAAAAACY8JDAAAAAAAkPAcncCoqqpycnURYayxk0zjZayxkUxjTXbJtK8Za+wk03gZa2wk01iTXTLta8YaO8k0XsYaG8k01mTi6D0wAAAAAAAAxoNLSAAAAAAAQMJzZAKjsbFR8+fPVyAQUE1NjROrjMjcuXO1ZMkSlZSUaNmyZfEezjDr1q1Tfn6+Fi9ePPS906dPq6ysTPPmzVNZWZnOnDkTxxH+h5HG+thjj2nWrFkqKSlRSUmJdu3aFccR/ofW1lbdfvvtWrhwoYqKivTDH/5QUmLu29HGmqj7tre3Vx/72Me0dOlSFRUV6dFHH5WUmPvWTcjd6CF3YyOZcldKruwld+OD3I0ecjc2yN3YIXcdFI6xgYGBcGFhYfjdd98N9/X1hYuLi8NvvvlmrFcbkY985CPhjo6OeA9jRM8991z4lVdeCRcVFQ197+GHHw5/73vfC4fD4fD3vve98N/8zd/Ea3jDjDTWRx99NPz9738/jqMa2QcffBB+5ZVXwuFwOHzu3LnwvHnzwm+++WZC7tvRxpqo+3ZwcDB8/vz5cDgcDgeDwfDHPvax8L//+78n5L51C3I3usjd2Eim3A2Hkyt7yV3nkbvRRe7GBrkbO+Suc2J+Bsb+/fsVCARUWFiojIwMrV27Vg0NDbFerWvdcsstmjp16rDvNTQ0qLKyUpJUWVmpnTt3xmFkVxtprInK7/frxhtvlCRNnDhRCxcuVFtbW0Lu29HGmqh8Pp9yc3MlSf39/erv75fP50vIfesW5G50kbuxkUy5KyVX9pK7ziN3o4vcjQ1yN3bIXefEfAKjra1Ns2fPHvq6oKAgYQ+8y3w+nz796U/rpptuUm1tbbyHY3TixAn5/X5Jl17oJ0+ejPOIrG3atEnFxcVat25dQp5GdfToUb366qu6+eabE37fXjlWKXH3bSgUUklJifLz81VWVpYU+zaZkbuxl2zHb6Jmw2XJlLtScmQvuesscjf2ku34TcRcuBK5G33krjNiPoERHuEhJz6fL9arjcgf/vAHHThwQLt379bmzZv1/PPPx3tIrrF+/Xq9++67OnjwoPx+v771rW/Fe0jDXLhwQatXr9bjjz+uSZMmxXs4lj481kTet6mpqTp48KCOHTum/fv364033oj3kFyN3MWVEjkbpOTKXSl5spfcdRa5iyslai5cRu7GBrnrjJhPYBQUFKi1tXXo62PHjmnmzJmxXm1ELo8vPz9fd911l/bv3x/nEVmbMWOG2tvbJUnt7e3Kz8+P84hGN2PGDKWmpiolJUVf/epXE2rf9vf3a/Xq1br33nt19913S0rcfTvaWBN13142ZcoU3XbbbWpsbEzYfesG5G7sJdPxm8jZkEy5KyVn9pK7ziB3Yy+Zjt9EzgVyN/bI3diK+QRGaWmpWlpadOTIEQWDQdXX16uioiLWqx237u5unT9/fuj/Tz/99LC7CieiiooK1dXVSZLq6uq0atWqOI9odJdfwJK0Y8eOhNm34XBYX/7yl7Vw4UJ985vfHPp+Iu7b0caaqPu2o6NDXV1dkqSenh79/ve/14IFCxJy37oFuRt7yXT8Jmo2JFPuSsmVveSu88jd2Eum4zcRc0Eid2OJ3HWQE3cK/e1vfxueN29euLCwMPyP//iPTqxy3N59991wcXFxuLi4OLxo0aKEG+/atWvD1157bTgtLS08a9as8E9+8pNwZ2dn+I477ggHAoHwHXfcET516lS8hxkOh0ce6xe/+MXw4sWLw0uWLAl//vOfD3/wwQfxHmY4HA6HX3jhhbCk8JIlS8JLly4NL126NPzb3/42IfftaGNN1H176NChcElJSXjJkiXhoqKi8N///d+Hw+FwQu5bNyF3o4fcjY1kyt1wOLmyl9yND3I3esjd2CB3Y4fcdY4vHB7hoj0AAAAAAIAEEvNLSAAAAAAAACLFBAYAAAAAAEh4TGAAAAAAAICExwQGAAAAAABIeExgAAAAAACAhMcEBgAAAAAASHhMYAAAAAAAgITHBAYAAAAAAEh4/x+b+01RTYtqewAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x1080 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<cil.utilities.display.show2D at 0x7f1156fb08b0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compare solution\n",
    "np.testing.assert_almost_equal(u_cil.array, u_cvx.value, decimal=3)\n",
    "\n",
    "# print objectives\n",
    "print(\"CVX objective = {}\".format(obj.value))\n",
    "\n",
    "f = 0.5*L2NormSquared(b=data)\n",
    "cil_objective = TV(u_cil) + f(u_cil)\n",
    "print(\"CIL objective = {}\".format(cil_objective))\n",
    "\n",
    "# show middle line profiles\n",
    "N, M = data.shape\n",
    "plt.figure()\n",
    "plt.plot(u_cil.array[int(N/2)], label=\"CIL\")\n",
    "plt.plot(u_cvx.value[int(N/2)], label=\"CVXpy\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "show2D([u_cil.array, u_cvx.value, np.abs(u_cil.array - u_cvx.value)], num_cols = 3, origin=\"upper\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27ead7da-2416-4ab9-8ea0-a803be230341",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:cil_cvxpy]",
   "language": "python",
   "name": "conda-env-cil_cvxpy-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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}