TotalGeneralisedVariation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d62f4331-efe2-4bcb-a55c-5e540721b716",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cil.optimisation.operators import GradientOperator, SymmetrisedGradientOperator, IdentityOperator, ZeroOperator, BlockOperator\n",
    "from cil.optimisation.functions import MixedL21Norm, BlockFunction, L2NormSquared, ZeroFunction\n",
    "from cil.optimisation.algorithms import PDHG\n",
    "from cil.utilities import dataexample\n",
    "from cvxpy import *\n",
    "from regularisers import tgv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from cil.utilities.display import show2D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "569d88e9-ec1d-4d6f-9bef-a6a40f5b7f0b",
   "metadata": {},
   "source": [
    "# Load Data and resize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "af88e27d-d47d-4c33-9db6-fe87332ed95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = dataexample.CAMERA.get(size=(32, 32))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95ff0e41-ccbc-49e7-9737-98ad84296cbc",
   "metadata": {},
   "source": [
    "# Setup and run tgv denoising with cvxpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c63d3590-e4c6-47fc-89fd-04ec31b9fe5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                    v1.1.17                                    \n",
      "===============================================================================\n",
      "(CVXPY) Dec 06 07:53:11 PM: Your problem has 3072 variables, 0 constraints, and 0 parameters.\n",
      "(CVXPY) Dec 06 07:53:11 PM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Dec 06 07:53:11 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Dec 06 07:53:11 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:53:11 PM: Compiling problem (target solver=SCS).\n",
      "(CVXPY) Dec 06 07:53:11 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS\n",
      "(CVXPY) Dec 06 07:53:11 PM: Applying reduction Dcp2Cone\n",
      "(CVXPY) Dec 06 07:53:11 PM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Dec 06 07:53:11 PM: Applying reduction ConeMatrixStuffing\n",
      "(CVXPY) Dec 06 07:53:11 PM: Applying reduction SCS\n",
      "(CVXPY) Dec 06 07:53:11 PM: Finished problem compilation (took 2.995e-02 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 06 07:53:11 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: 5121, constraints m: 9218\n",
      "cones: \t  q: soc vars: 9218, qsize: 2049\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): 20994, nnz(P): 0\n",
      "------------------------------------------------------------------\n",
      " iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)\n",
      "------------------------------------------------------------------\n",
      "     0| 1.36e+01  5.11e-01  2.58e+03 -1.29e+03  1.00e-01  7.93e-03 \n",
      "   250| 2.58e-03  1.79e-04  7.20e-03  8.46e+00  1.00e-01  4.30e-01 \n",
      "   500| 1.06e-03  2.19e-05  1.66e-03  8.47e+00  1.00e-01  1.16e+00 \n",
      "   750| 8.23e-04  2.36e-05  9.50e-04  8.47e+00  1.00e-01  1.91e+00 \n",
      "  1000| 5.55e-04  4.07e-06  6.86e-04  8.47e+00  1.00e-01  2.65e+00 \n",
      "  1250| 4.14e-04  6.05e-06  4.16e-04  8.47e+00  1.00e-01  3.19e+00 \n",
      "  1500| 4.23e-04  5.41e-07  3.15e-04  8.47e+00  1.00e-01  3.69e+00 \n",
      "  1750| 2.66e-04  6.43e-06  1.95e-04  8.47e+00  1.00e-01  4.37e+00 \n",
      "  2000| 1.80e-04  2.44e-07  1.37e-04  8.47e+00  1.00e-01  5.18e+00 \n",
      "  2250| 1.42e-04  1.95e-07  1.15e-04  8.47e+00  1.00e-01  5.83e+00 \n",
      "  2500| 1.32e-04  1.65e-07  1.02e-04  8.47e+00  1.00e-01  6.60e+00 \n",
      "  2750| 2.92e-01  2.51e-02  3.18e-06  8.47e+00  1.00e-01  7.31e+00 \n",
      "  2875| 1.01e-04  4.41e-06  6.91e-05  8.47e+00  1.00e-01  7.66e+00 \n",
      "------------------------------------------------------------------\n",
      "status:  solved\n",
      "timings: total: 7.69e+00s = setup: 3.15e-02s + solve: 7.66e+00s\n",
      "\t lin-sys: 5.91e+00s, cones: 5.44e-01s, accel: 4.64e-01s\n",
      "------------------------------------------------------------------\n",
      "objective = 8.473567\n",
      "------------------------------------------------------------------\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Dec 06 07:53:18 PM: Problem status: optimal\n",
      "(CVXPY) Dec 06 07:53:18 PM: Optimal value: 8.476e+00\n",
      "(CVXPY) Dec 06 07:53:18 PM: Compilation took 2.995e-02 seconds\n",
      "(CVXPY) Dec 06 07:53:18 PM: Solver (including time spent in interface) took 7.711e+00 seconds\n"
     ]
    }
   ],
   "source": [
    "# solution\n",
    "u_cvx = Variable(data.shape)\n",
    "w1_cvx = Variable(data.shape)\n",
    "w2_cvx = Variable(data.shape)\n",
    "\n",
    "# regularisation parameters\n",
    "alpha0 = 0.1\n",
    "alpha1 = 0.3\n",
    "\n",
    "# fidelity term\n",
    "fidelity = 0.5 * sum_squares(u_cvx - data.array)   \n",
    "regulariser = tgv(u_cvx, w1_cvx, w2_cvx, alpha1, alpha0)\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": "c3948c92-bf1c-4900-bfc1-c1ff1e4e0538",
   "metadata": {},
   "source": [
    "# Setup TGV denoising using CIL and the PDHG algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d37c9e67-0f4a-4cbb-b4f4-ceb0f2eef4a7",
   "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       1000         0.000    2.15422e+02  -0.00000e+00    2.15422e+02\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/cil_cvxpy/lib/python3.9/site-packages/cil/optimisation/algorithms/PDHG.py:477: UserWarning: Convergence criterion of PDHG for scalar step-sizes is not satisfied.\n",
      "  warnings.warn(\"Convergence criterion of PDHG for scalar step-sizes is not satisfied.\")\n",
      "/opt/anaconda3/envs/cil_cvxpy/lib/python3.9/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       1000         0.008    8.48440e+00   8.46977e+00    1.46298e-02\n",
      "     1000       1000         0.008    8.47692e+00   8.47331e+00    3.61393e-03\n",
      "----------------------------------------------------------------------------\n",
      "     1000       1000         0.008    8.47692e+00   8.47331e+00    3.61393e-03\n",
      "Stop criterion has been reached.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ig = data.geometry\n",
    "\n",
    "K11 = GradientOperator(ig)\n",
    "K22 = SymmetrisedGradientOperator(K11.range)\n",
    "K12 = IdentityOperator(K11.range)\n",
    "K21 = ZeroOperator(ig, K22.range)\n",
    "K = BlockOperator(K11, -K12, K21, K22, shape=(2,2) )\n",
    "\n",
    "f1 = alpha1 * MixedL21Norm()\n",
    "f2 = alpha0 * MixedL21Norm()\n",
    "F = BlockFunction(f1, f2)\n",
    "G = BlockFunction(0.5 * L2NormSquared(b=data), ZeroFunction())\n",
    "\n",
    "sigma = 1./np.sqrt(12)\n",
    "tau = 1./np.sqrt(12)\n",
    "\n",
    "# Setup and run the PDHG algorithm\n",
    "pdhg_tgv = PDHG(f=F,g=G,operator=K,\n",
    "            max_iteration = 1000, sigma=sigma, tau=tau,\n",
    "            update_objective_interval = 500)\n",
    "pdhg_tgv.run(verbose = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7a3dd7e-6501-4f1f-b179-7dfe8cf1b00a",
   "metadata": {},
   "source": [
    "# Compare solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "04b98efe-e3b2-461a-a10a-c72e116f3d40",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CVX objective = 8.476068484702772\n",
      "CIL objective = 8.476922034844755\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwM0lEQVR4nO3deXyU5b3//9cnk30PWUlCFiBA2JEQcQWtIrgh1VbQaje1nOo55fhr3VpbTj09X2t7rJ5qy1Fra904agGpohataFXQBAxbwhKWQICsZCN7MtfvjwQaY4AhTHLNTD7Px4MHmXvuTN7DJG+u3HPd1y3GGJRSSnk/P9sBlFJKuYcWulJK+QgtdKWU8hFa6Eop5SO00JVSykf42/rCcXFxJiMjw9aXV0opr7Rx48YqY0x8X/dZK/SMjAzy8/NtfXmllPJKIlJysvv0kItSSvkILXSllPIRWuhKKeUjtNCVUspHaKErpZSP0EJXSikfoYWulFI+wto8dKWUd2lpaqC26ggN1UdorimjvbEGZ2cHdHZgnN1/OtvB2QmdHeDs+tgvMpH0c68jMS3L9lPweVroSg1xzs5OKo/so3J/IU1lu3Ee3Y+jqYrA1mpC2msI76wl2llHqLSSBCT154ts/0/2+mVQnjSLmKnXkHXOJTj8tX7czaV/URGZCzwOOIBnjDEP97o/CngBSOt+zF8bY/7o5qwAFH32Lo0fPA5I1x/p+tsIX9rm9A/BGZ1GYNxIIpPHkJA+jqiYuIGIpZTHqy4vpax4M41HduKsLiaovoTo5oMkdR4hUdpI7N6vzTiolSgaHNE0BsRQF5rOgeBYTFgcjvB4AqISCYlOJCRiGP4BwTgC/HE4AvDz98ffPxCHvz/+/gE4/ANwOAI4WLyFI3kriTjwd2Yceh7/w89RsyaCPZEzYexcss67jqhh+nPpDnK6KxaJiAPYBVwOlAJ5wCJjTGGPfR4Aoowx94pIPLATSDLGtJ3scXNyckx/Tv3fuu41Ij9cimAQurILBswXbwuGYNNCDPVf+Pw6wqjwH059cCrtkWk4YjOJzZpJ5sSZiJ++paB8T3tbK/nPP0DugT/gkK6fkTbjzxFHEkeDR9ASkYHEjiJs+BhiR4wjMXXUgI2e62qqKF7/Os6dbzO6bgMx1NNh/NgZNIGGjLmMmv0N4pMzBuRr+woR2WiMyenzPhcK/TxgqTHmiu7b9wMYY/5fj33uB0YAdwIZwFpgjDHGebLH7W+hn6nG+qOUH9hF7aHdtFXuQWr2EdJYSkzrYRKd5QRKJwAH/FI4nHo1abO+SfKoCQOeS6nBUFq8jcbl32Fsx07youYQPP1mYtOyB7S0XdXZ0cHuz9dRU/AGiUfeZ6RzP04jFAZNpjFrPmMvuZnouH4d4PFpZ1voNwBzjTG3dd++BTjXGHNXj30igNXAOCACuNEY82Yfj3UHcAdAWlra9JKSk64xMyg62tupOLSXg/lrCNu1gvGtW/ETwy7/MdSMuo7Rl95KbOIIqxmV6g/jdJK38nEmbvl/tIs/xbkPMf3K79qOdUolOws4/NELpJa+yQhzmHbjoDB0Ou3ZX2Xc7IWER8bYjugRzrbQvwZc0avQc40x/9pjnxuAC4C7gVF0jdCnGGPq+3hIYPBG6GfiyME97Hv/zySUrGZ05146jVAYcg6t2deTfckiwiKH2Y6o1GnVVB5h/x+/y7Smj9kWNJX4W54lMXWU7VguM04ne7aup3LDS2QeeYckKmkxARRGnIdMuoHsi64nODTcdkxrBuOQy5vAw8aYf3Tf/jtwnzHms5M9ricWek/7Cjdy5OPnST+8hhRTTosJYOuk+5hxww9tR1PqpLas+wvJ6+4m0hxjU9a/krvoQfwcDtux+s3Z2cmu/Peoy1vO6Kp3iaWOYyaEouiLCZz6NcZfcC0BgUG2Yw6qsy10f7reFP0KcIiuN0VvMsZs77HP74FyY8xSEUkENtE1Qq862eN6eqEf5+x0smPj32HtT0lt34f5wVadKaM8TkvTMQr+uISZla+y328EnQueYdSkmbZjuVVHextF69+kedMrjKtdRyRN1BDBrthLCZ++kOxzr/Dq/7xcdVaF3v0AVwKP0TVt8VljzC9EZDGAMWaZiCQDfwKG0zV38GFjzAunekxvKfTjiresZ/SKuWwYcTszv/tr23GUOmHvtk/xW3EbGc4DbIj/GlO//ZjPH5JobWmi6B+r6NzyKtn1HxMqrVQwjL2Jc4g59ybGTL3IZ2etnXWhDwRvK3SATY9cxeimTTpKVx6jubGBtl+NpZ0ADs9+lMmzr7cdadA1Hauj6INXcGxfwfjGTwmUTkoliYMpV5My61bSsqbYjuhWpyp03/wvbIBEzv0xkTRRtOLh0++s1CAo/OAVomik7PInhmSZA4SGRzH9qtuZes9bNC/ZSd6Uh6gNHE7uwT+Q9uLF7PrPGWx46RdUlR20HXXA6Qj9DOkoXXmSzx+5ktSmQob9pNj6vHJPU3V4P8XvP0fc3tcZ3bmHDuPH9pDptE24gfGzFxEWEWU7Yr/oCN2NTozSV/7SdhQ1xNXXVjOh8VP2JMzRMu9DXHIGM2/+GaMf3ETJjX8nL/VWElr2M2PjvfDrLPIevYHN6/5CR3u77ahuo98FZ2j05PP5/O0LGH/gBepr7yUyWkfpyo4d779ErnQwbOZNtqN4vPTs6aRnT8fZ+ShF+Wtp+PRFxh59j6h1a6lcF8PuxCtJvOhbjJqYazvqWdERej9EXNE1Si9c8YjtKGoIC965kkOSSNbUi21H8Rp+DgfZ584l99+eJ/i+Yj4/77ccCs1mRtlyRr12OcUPncP6Fx+istw7j7droffD6CkX8Hno8VH6SafaKzVgqstLGd/8OQeS5/ns9LyBFhQcyrQrbmXqPW9x7K7tfDbuXoz4cd7uXxPzu8kU/PIK8tc8S0tzo+2oLtPvhH7qGqU36ihdWVG87kX8xUnSBd+wHcUnxMQPJ3fhA2T9JJ+ShX8nP+UbJDfvIuezf6ft4Sw+ffwWtn30Bp2dnbajnpLOcjkLnz8yj1FNm2HJViKjY23HUUNI0S/OJ6SzgYyfbrUdxWd1dnRQtP4NWvNeILvuwxMnL+1JmEPMzJsYa+nkJZ3lMkAirvhJ1yh9pY7S1eApO7Cb7PbtHEm7ynYUn+bw92fiRdcx/e7XkB8Vs2nGrzkcOo7p5a8ybvW1HHpoPJ88/e/sK9pkO+oJOkI/SzpKV4Ntwws/Y2bxYxy69RNSRura/YOtrqaSXe+/SMjOVWS3FOAQwx6/TMrTryZl5g2kZU0e0JG7nvo/gHYXfETWqqvYkL6Ymd/Wuelq4BU/dA5OcTDmJ3m2owx5VWUH2PP+80TvXc3Y9h0AVBJDScQ0OtIuZPiUy0kbPdGtBa+FPsA+f2QuI5u24LdkKxE6SlcD6ODuzYx48WI2ZP1/zLz5p7bjqB6O7N9Baf4bOEo+Jq1hE3HUAlDBMEoiptGZdoFbCv5Uha4nFrlB+JyfELXqKjasfERH6WpAlf7jRVKMMHL2LbajqF6GZ4xjeMY44IcYp5PSPVs5snktjpKPyWzYSNz292D7z7tWhRxzGzNv+rHbM2ihu0HW1Asp+Nt5ZJc8T0PtPTpKVwPCOJ0kl75JUdBEJqRk2o6jTkH8/EjNmkJq1hT6KviA8IG5+pkWupuEzfkxUauu1lG6GjB7t21glLOUT7O+bTuKOkO9C36g6LRFN8maehEFId2j9LqjtuMoH1Sx/mXajYMxs2+2HUV5KJcKXUTmishOESkWkfv6uP9HIlLQ/WebiHSKyJC7onLo5Q8QRSOF7zxjO4ryMcbpJOPIWxSGnENM/HDbcZSHOm2hi4gDeBKYB4wHFonI+J77GGN+ZYyZaoyZCtwPfGCMGXLD1KypF3LMhEDlTttRlI/Zmf8ew6mkLfurtqMoD+bKCD0XKDbG7DXGtAHLgfmn2H8R8LI7wnkb8fOj3H84wccO2I6ifExd3nJaTADjZi+0HUV5MFcKPQXouZZkafe2LxGRUGAu8JeT3H+HiOSLSH5lZeWZZvUKdcEpRLceth1D+ZCO9jZGVb5LYfh5REQNuSOZ6gy4UujSx7aTnY10DfDxyQ63GGOeMsbkGGNy4uPjXc3oVdoi0kjqLMfp4auyKe9RtP6trpNUJg3Na4Yq17lS6KXAiB63U4GTDUEXMkQPtxwnsSMJknYqDu+zHUX5iObP/49jJoTxF99gO4rycK4Ueh6QJSKZIhJIV2mv7r2TiEQBs4DX3RvRu4QmZgFQfUDfGFVnr7WliXE16yiKvojg0HDbcZSHO22hG2M6gLuAd4Ai4BVjzHYRWSwii3vsugD4mzHGey7vMQCGpY4BoKlst+UkyhcU/WMVkTQSOPXrtqMoL+DSmaLGmDXAml7blvW6/SfgT+4K5q0SR4yi3TjoOKqHXNTZ69zyKjVEMP6Ca21HUV5AzxR1M/+AQMr94gmsL7EdRXm5pmN1ZNd/zK7YSwkIDLIdR3kBLfQBcDQohcjmUtsxlJcr/OAVQqWV8Ok691y5Rgt9ADSHjSChQ+eiq7Pj2L6CCoYxLneO7SjKS2ihDwATk0EUjdQd9c2Tp9TAqztayYTGT9mbOAeHvy6KqlyjhT4AAuNHA1BxYIflJMpb7Xz/BQKlk9jzvmE7ivIiWugDIDqla+pi/eFdlpMobxW2cyUHJZnRky+wHUV5ES30AZCYPhaAtso9lpMob1RxaB/ZrVsoHXH1gF49Xvke/W4ZAGER0VQRjaN2v+0oygvtXfdn/MSQepEeblFnRgt9gFT5DyesUacuqjMXu3c1ux2jGZE1xXYU5WW00AfIsdBUYtt06qI6Mwd2FZDVWUz1SD0zVJ05LfQB0h6VQYKporWlyXYU5UUO/eMFnEYYOftW21GUF9JCHyD+cSPxE0P5AZ3polxjnE5SS9+kKGgyCSmZtuMoL6SFPkAikrrmoteUaqEr1xRv+ZgR5jCNYxfYjqK8lBb6AIlLHwdAS4VOXVSuqV7/Im3GwdhLdHaL6h8t9AESm5BKkwnC6DK6ygWdHR2MLH+H7WHnEjXMNy/PqAaeFvoAET8/yh3DCT52wHYU5QWKPn2LBI7inKiXmVP9p4U+gGqDU4hu0bno6vSa8pfTaIIZP0uvTKT6z6VCF5G5IrJTRIpF5L6T7DNbRApEZLuIfODemN6pNSKNpM4yjLPTdhTlwbquG/o+hdEXExIWYTuO8mKnLXQRcQBPAvOA8cAiERnfa59o4HfAtcaYCcDX3B/V+8iwTIKlnaqyg7ajKA9W+OEKImkkaOqNtqMoL+fKCD0XKDbG7DXGtAHLgfm99rkJWGGMOQBgjKlwb0zvFJLYNXWxUpfRVafg3PIqR4kk+4JrbEdRXs6VQk8Beg4xS7u39TQGiBGRdSKyUUT6PM1NRO4QkXwRya+s9P2LPwxL7Vp18diR3ZaTKE91rL6GCQ0fszvuMr1uqDprrhS69LHN9LrtD0wHrgKuAB4UkTFf+iRjnjLG5BhjcuLjfX9qVsKILDqN0Fm913YU5aGK3n+ZYGknKvcm21GUD3Dl2lalwIget1OB3qtOlQJVxphGoFFEPgSmAEP6NMnAoCAO+8UTWF9iO4ryUEFFKzgsCYzN+YrtKMoHuDJCzwOyRCRTRAKBhcDqXvu8DlwkIv4iEgqcCxS5N6p3OhqYTHiTTl1UX1ZdXsr45o2UDJ+nF7JQbnHa7yJjTAdwF/AOXSX9ijFmu4gsFpHF3fsUAW8DW4DPgGeMMdsGLrb3aAxLI6HjiO0YygMVr3sBf3GSdOEttqMoH+HS5cSNMWuANb22Let1+1fAr9wXzTeY6AxijtbTUHeUiKhhtuMoDxJV/Dr7/DLIHD/DdhTlI/T3vAEWGD8SgPL9egRK/dPh/TsZ115IWfrVtqMoH6KFPsAik7sm+zQcHtLvD6teSj74MwDps/RCFsp9tNAHWGJGNgCtlTp1Uf1TUskb7AgYT3LGWNtRlA/RQh9gEVHDqCESqd1vO4ryEPu2f0qmcz91o3ufcK3U2dFCHwQV/sMJa9RldFWXiveeoMP4MXq2XshCuZcW+iBoCEllWFvvc7HUUPTp/z3MuUdXk59wPbGJqbbjKB+jhT4I2qPSSXRW0t7WYjuKsqhg7UvkFD5MQchMcu74ne04ygdpoQ8CR+xIHGIoP1BsO4qyZNemdYz96AfsCchizJ2v4B8QaDuS8kFa6IMgfHgWAEdLd1pOomw4tHc7catvocYvhmG3rSA0PMp2JOWjtNAHQdyIcQA0l+sIfaipqTyC8/kb8MNJx6JXiEsacfpPUqqftNAHQVzSCFpMAOboPttR1CBqaTpG+VNfJcFZSdm8P5I2ZqrtSMrHaaEPAj+HgzJHEkENOnVxqOjs6KDwyUWMaSti+8xfMe7cObYjqSFAC32Q1AanEt1yyHYMNUjynvo+5zR+yGdj/p1z5n3bdhw1RGihD5KW8DQSO49gnE7bUdQA2/DSQ8ys+D82xH+Ncxc9aDuOGkK00AdLTAah0kp1hV7swpdtevtP5O78bz4Pu5AZ31umF65Qg0q/2wZJSOJoAKpKdlhOogZK8eaPGL/+h+wOGEv2nf+Hw9+lyw0o5TZa6IMkJrV7Gd0ynbroq+rXPkKrBBJ/xwqCQ8Ntx1FDkEuFLiJzRWSniBSLyH193D9bROpEpKD7z0/dH9W7JaaNwWmEjipdRtcXVZUdZFLDRxQlXsOwhBTbcdQQddrfCUXEATwJXA6UAnkistoYU9hr138YY/TyKycRFBxKmcQRULffdhQ1AHa/s4zzpJPhl/6L7ShqCHNlhJ4LFBtj9hpj2oDlgC7k3A/VgcMJb9I3RX2Ns7OTtP2vsj1wEuljp9qOo4YwVwo9BTjY43Zp97bezhORzSLylohM6OuBROQOEckXkfzKysp+xPVujWEjiO/QZXR9zfaPXifFlNM8WS8np+xypdClj22m1+1NQLoxZgrwW2BVXw9kjHnKGJNjjMmJj48/o6C+oDM6g1jqaGyotR1FuVH7Z3+khkgmXaYXrFB2uVLopUDPFYVSgS8MM40x9caYY90frwECRCTObSl9RGD8KADKS3TVRV9RVXaAScc+ZmfSNQQFh9qOo4Y4Vwo9D8gSkUwRCQQWAqt77iAiSSIi3R/ndj9utbvDervI7mV0aw/tspxEuUvx28sIkE5SvvI921GUOv0sF2NMh4jcBbwDOIBnjTHbRWRx9/3LgBuAfxGRDqAZWGiM6X1YZshLSM8GoK1S56L7AmdnJ2klr7I9cAoTsqbYjqPU6QsdThxGWdNr27IeHz8BPOHeaL4nalg8dYQhNfttR1FusO3DlUw2FRyZcq/tKEoBeqbooKvwTyak8eDpd1QerzPvWY7qm6HKg2ihD7L6kFSGteoyut6u8vB+JjWuZ+fwawkMCrYdRylAC33QtUWkkeispKO9zXYUdRaK3/k9/uJkxFf0zFDlObTQB5kjdiQB0kn5wT22o6h+6uzoILPkNbYGTSN19ETbcZQ6QQt9kIUldU1drCnVuejeatuHr5FEFe1Tv2k7ilJfoIU+yGLTxgLQWK5TF72Vyf8TVUQz6Ss32Y6i1BdooQ+yhORM2ow/plqX0fVGZQeLmdS4gd3J8wkIDLIdR6kv0EIfZH4OB2WORAIbDtiOovph39+WIUD65d+3HUWpL9FCt6AmKIWoFl1G19t0tLcx6uBf2BYyneTMcbbjKPUlWugWtISnkdhRhnE6bUdRZ2DbutdI4Cgd075lO4pSfdJCt8DEZBIuzdRUHbEdRZ0B2fRHKolh0iVftx1FqT5poVsQnNC1jG7FgR2WkyhXHSnZyaSmPIpTF+ibocpjaaFbEJ3SNRf9WLnOdPEWJWv/F4CMyxdbTqLUyWmhWxCfOhqA9qr9doMol3S0tzGydCVbQ2cwPH2s7ThKnZQWugVhEdHUEIlfva666A22vv8KCRzFOU3PDFWeTQvdkipHAiGNuuqiN3Bs+hMVDNM3Q5XHc6nQRWSuiOwUkWIRue8U+80QkU4RucF9EX1TQ3Ay0W06y8XTHd6/k4nN+ewZ8VX8AwJtx1HqlE5b6CLiAJ4E5gHjgUUiMv4k+/2SrkvVqdNoC08hobNC56J7uJK1v8MAmXN0mVzl+VwZoecCxcaYvcaYNmA5ML+P/f4V+AtQ4cZ8Pkti0gmWdqor9LCLp2pvayXr0Cq2hp5L0ojRtuModVquFHoK0PPdu9LubSeISAqwAFjGKYjIHSKSLyL5lZWVZ5rVpwTFZwJQXbrbchJ1Mlv//jJx1ELOt21HUcolrhS69LHN9Lr9GHCvMabzVA9kjHnKGJNjjMmJj493MaJvikrqKnSdi+65Aj5/jjLimDRL3xJS3sHfhX1KgRE9bqcCh3vtkwMsFxGAOOBKEekwxqxyR0hfFJfadXJRe3WJ5SSqL4f2bmdS6ybWp32PJH9XfkyUss+V79Q8IEtEMoFDwELgCyv7G2Myj38sIn8C3tAyP7WIqGHUEo7U6TK6nujA2mUkGj9GXaFvhirvcdpDLsaYDuAuumavFAGvGGO2i8hiEdHzoM9ClSORYJ2L7nHaWlsYc+R1tobNJCEl8/SfoJSHcOl3SWPMGmBNr219vgFqjPnW2ccaGuqDk4lt3m87hupl63svM506Smd8x3YUpc6InilqUddc9HKdi+5hgjb/iTLimXjRAttRlDojWug2RacRIm26LroHKS3exsTWAvalX49D3wxVXkYL3aLguK7js1WlxZaTqONK3/0dHcaP0fpmqPJCWugWRQ7vutBFQ5nORfcErS1NjC37K1vDzyc+OcN2HKXOmBa6RXHH10Wv3mc5iQLY+t6LxFCPQ98MVV5KC92iyOhY6glD6nRddE8QsvnPHJZEJl50ne0oSvWLFrpllToX3SMc2FXAhLYtlGTcgJ/DYTuOUv2ihW5ZfdBwolp1lotth99bRrtxkHWFniunvJcWumWt4ak6F92yluZGxpa/wdaIC4hLSrMdR6l+00K3LTqNUGmltrrcdpIha9u7LxBDAwG5t9mOotRZ0UK3LCguA4AqXRfdmtCtz3NIEplwwdW2oyh1VrTQLYtMGglAg66LbkXJjk2Mb9vKgcwb9c1Q5fW00C2LGzEGgLaq/XaDDFFH/r6MNuNgzBXfsx1FqbOmhW5ZVEwc9YQitbou+mBrbKglu+INtkZcRGxiqu04Sp01LXQPUOlIJEjnog+6La/+gigaCZv9A9tRlHILLXQPoHPRB19V2QGmlDzHxvBZjMu51HYcpdxCC90DtIalEN9ZoXPRB9GeV39CAB0kXPdftqMo5TYuFbqIzBWRnSJSLCL39XH/fBHZIiIFIpIvIhe6P6oPi04jXJqpr6m0nWRIKNmxielVf2VTwgJGjJ5oO45SbnPaQhcRB/AkMA8YDywSkfG9dnsPmGKMmQp8B3jGzTl9WmD3uuiVOhd9UBxd/WOaCSLraz+3HUUpt3JlhJ4LFBtj9hpj2oDlwPyeOxhjjhljTPfNMMCgXHZiLnrZHstJfF/hhreZ1vQJ20Z+h2EJKbbjKOVWrhR6CtBzfdfS7m1fICILRGQH8CZdo/QvEZE7ug/J5FdW6uGF4+K710Vv1bnoA8o4nfi/+yAVDGPqDffbjqOU27lS6NLHti+NwI0xK40x44DrgIf6eiBjzFPGmBxjTE58fPwZBfVlkTHxHDMhUKvrog+kTW8/x5iOXeyf8u+EhEXYjqOU27lS6KXAiB63U4HDJ9vZGPMhMEpE4s4y25Ahfn46F32AtbW2kJj3MPv8Mph+zfdtx1FqQLhS6HlAlohkikggsBBY3XMHERktItL98TlAIFDt7rC+rC54OFEtJ/1/Up2lTSv+m1RTRsNFD+Lw97cdR6kBcdrvbGNMh4jcBbwDOIBnjTHbRWRx9/3LgOuBW0WkHWgGbuzxJqlyQWtYCvGNBRinE/HT0wPcqb62mrE7f8+2oKlMmvVV23GUGjAuDVWMMWuANb22Levx8S+BX7o32tBiokYQUdlMXW01UcP0/QV32v7KUs6jgeorf6H/WSqfpt/dHiJQ10UfEGUHi5l26GXyIi9n9BQ93035Ni10DxGRNAqA+jJdF92dDv7lJ/hhSL3+F7ajKDXgtNA9RHxqFgCtVfssJ/Ede7ZuYHrN22wafiPD08fajqPUgNNC9xBRwxJoNMGg66K7TeObP6ZBQsn++n/YjqLUoNBC9xBdc9ETCDpWajuKT9j64Uomt+RTlPU9fZNZDRla6B6kLmg4ka1ltmN4PeN0EvzBf3JYEph2/Y9sx1Fq0Gihe5CWsBTiO8ttx/B6O/PWktVZzMHxiwkKDrUdR6lBo4XuQUx0GpE0UldTZTuKV2v86PfUE8akebfZjqLUoNJC9yCBsRkAVJUW2w3ixaoO72dy/YcUJl5LaHiU7ThKDSotdA8Snti1Lnr9ES30/tr91hM4cJJ6+V22oyg16LTQPUjc8XXRq0ssJ/FOba0tZB18la2hM0jVS8upIUgL3YPExA2nyQTpXPR+2vLu88RRC7l32I6ilBVa6B5E/Pyo0Lno/RZW8CylksSki3VFRTU0aaF7mLqg4US0HLEdw+sUb/6Y7PZCSkffjJ/DYTuOUlZooXuYltBk4jv15KIzdXTdkzSZILLn/YvtKEpZo4XuYUxUGlE00lB31HYUr1FXXc7ko39ja+wVepq/GtJcKnQRmSsiO0WkWETu6+P+m0VkS/efT0RkivujDg0Bui76GSt663cESzvxl+pURTW0nbbQRcQBPAnMA8YDi0RkfK/d9gGzjDGTgYeAp9wddKg4vi563RFdF90VnR0dpO15mcLASYyceK7tOEpZ5coIPRcoNsbsNca0AcuB+T13MMZ8Yoyp6b65AUh1b8yhIzalay56S6Wui+6KrR+8RrIpp3nqd21HUco6Vwo9BTjY43Zp97aT+S7w1tmEGsqGxSfTbAJ1LrqL/PKepoJhTL7sJttRlLLOlUKXPraZPncUuYSuQr/3JPffISL5IpJfWVnpesoh5Phc9MBjh2xH8XgHd29mcks+e9K/TkBgkO04SlnnSqGXAiN63E4FDvfeSUQmA88A840x1X09kDHmKWNMjjEmJz5eZyOcTF3gcCJavvRPrHo5tPYJ2oyDrHl32o6ilEdwpdDzgCwRyRSRQGAhsLrnDiKSBqwAbjHG7HJ/zKGlOSyFOF0X/ZQaG2qZUP5XtkRdQlxSmu04SnmE0xa6MaYDuAt4BygCXjHGbBeRxSKyuHu3nwKxwO9EpEBE8gcs8RDgjEojhgaO1decfuchattbTxMhzYRfpCcSKXWcvys7GWPWAGt6bVvW4+PbAL2agJsExKbBXqgsLSZ8/AzbcTyOcTpJ3PFnih2jGDv9UttxlPIYeqaoBzqxLnrZHstJPFPh+rfIcB7g6MRvIX76LazUcfrT4IHiUrMAaKncbzeIh2r9ZBm1hDP5iu/YjqKUR9FC90CxCam0mABMjc5F763sYDGTj33EjqT5BIeG246jlEfRQvdA/5yLfvD0Ow8x+97+LX4Y0q74N9tRlPI4WugeqjYwSddF76W+tpoJh15lc9j5JGeOsx1HKY+jhe6hmkNTiO3Queg9bV/5CJE0EnHFA7ajKOWRtNA9lDNqBMOop+lYne0oHqGh7ijjS56nIPQ8Rk+50HYcpTySS/PQ1eALiM2AfVB5sJj07Om241i3beWvOI9Gwi7X0bmvam9vp7S0lJaWFttRPEJwcDCpqakEBAS4/Dla6B7q+Fz02iN7hnyhH6uvIXv/n9kcksuUaRfbjqMGSGlpKREREWRkZCDS15qAQ4cxhurqakpLS8nMzHT58/SQi4fSddH/aevKXxPNMUIu/7HtKGoAtbS0EBsbO+TLHEBEiI2NPePfVrTQPVRs4gjajD/OIb4uemNDLeP2PceW4BmMOWe27ThqgGmZ/1N//i200D2Un8NBuV8CgQ2ltqNYtWXVo8TQQNBlX7qUrVKqFy10D9Y1F33orovedKyOsXv+yNagcxibc5ntOGqIKCsrY+HChYwaNYrx48dz5ZVXsmvXLiZOnAjAunXruPrqqy2n7JsWugdrDhvac9E3r/oNw6gn4Cv3246ihghjDAsWLGD27Nns2bOHwsJC/uu//ovycu/4OdRZLh7MGZlG7NE66mqriYqOtR1nUDU3NjCm+Fm2BU1lYu4c23HUIPuPv26n8HC9Wx9zfHIkP7tmwin3ef/99wkICGDx4sUntk2dOpX9+/e7NctA0RG6B4udMheAwtW/sZxk8BWseoxY6vC7REfnavBs27aN6dO9d5qwjtA9WNa0i9nyTg7j9j5H07EfERoeZTvSoGhpOkbW7mfYHjiZCTPn2o6jLDjdSFr1zaURuojMFZGdIlIsIl+abiAi40RkvYi0isgP3R9z6Aq89D5iqGfLqsdsRxk0BaseJ45aZLbObFGDa8KECWzcuNF2jH47baGLiAN4EpgHjAcWicj4XrsdBf4N+LXbEw5x43IvZ1vQVEYXP0tL0zHbcQZcS3MjI3c9Q2HgJMaff5XtOGqIufTSS2ltbeXpp58+sS0vL4+SkhKLqVznygg9Fyg2xuw1xrQBy4H5PXcwxlQYY/KA9gHIOOTJrHuJo5aC1//HdpQBV7Dqf0jgKM6L77EdRQ1BIsLKlStZu3Yto0aNYsKECSxdupTk5OQv7Pfee++Rmpp64s/69estJf4iV46hpwA9r7RQCpzbny8mIncAdwCkpaX15yGGpAnnX0nhuomM3Pk0Lc1LCA4JtR1pQLQ0NzFy51MUBUxgwvmeOc9X+b7k5GReeeWVL23ftm0bALNnz6a5uXmwY7nElRF6X+efmv58MWPMU8aYHGNMTnx8fH8eYsjqvPAeEjhKwV+ftB1lwGxe/QQJHKXzonv04s9K9YMrPzWlwIget1OBoXv6oiUTL7yGnf7jyCj8X9pafW950daWJtKL/pcdAdlMuPBa23GU8kquFHoekCUimSISCCwEVg9sLNWb+PnRcsEPSaKSgjeW2Y7jNi3NjWz62wsUPv5Vkqii/YIf6ehcqX467TF0Y0yHiNwFvAM4gGeNMdtFZHH3/ctEJAnIByIBp4gsAcYbY9x7qtcQN3nW9ez++JekbPsd7dd+n4CAQNuR+qW58RhF/1hB5/ZVZNd/wjnSTC3hbBh+C+devMB2PKW8lksnFhlj1gBrem1b1uPjMroOxagBJH5+NJ57N1kf/wufvfk0udfdaTuSy5oa6yn6cCVm+yqyGz7hHGmhhggKYy8jdOr1jJt5JTMDg2zHVMqr6ZmiXmbKVxayZ8MjJG1+gs6rv4fD33NfwpqqMoo/WYXfrrcY17Ce6dJKDZFsj51D6LQbGDdzHrle+luGUp7Ic9tA9Un8/KjP+QHTPl1C/tvPknP1HbYjnWCcTvYV5lG+8XWiS99nTFsRM8RQTRTb4uYRNu16smfOI9ff9WskKjXYysrKWLJkCXl5eQQFBZGRkcFbb73Fjh07GDt27In9lixZQnJyMqNHj+bJJ5/kvffeA+Cjjz7irrvuIj8/H/9BHnBpoXuhKXNuZX/efxO36X9wzvsufg6HtSzNjQ3s3PAmrYVvkV79ESOpYiSw2380eWnfZdi0axg95SLOtZhRKVcdXz73m9/8JsuXLwegoKCA5uZmli9fzs9+9jMAnE4nr732Gh9//DHp6en84Q9/4KWXXuLrX/863//+91m2bNmglzlooXslP4eD6nN+wPT8H7Lxb39m+rxvD9rXdnZ2snfbBqq2riWk9CPGNhcwVdppNMHsCs/h4MjLyTxvAVnJ6WQNWirlc966D8q2uvcxkybBvIdPucvJls99/PHHufHGG08U+ocffkhGRgbp6ekA/Pa3v+Wyyy5j+/btzJgxg/PPPx+ApUuXsmfPHg4dOsTBgwe55557uP3227nlllu44YYbmD+/66T7m2++mRtvvJFrrz27Kbta6F5q6txvc3Djo8TkPYZzzq0DNko3TiclOzdRvvlvBB78mJGNnzOaRkYDByWZgsQFhE28mjG5lzMt2DfPYFVDx8mWz508eTJ+fn5s3ryZKVOmsHz5chYtWnTi/pEjR3LjjTfyxBNPsGfPni987pYtW9iwYQONjY1MmzaNq666ittuu43f/OY3zJ8/n7q6Oj755BOee+65s86vhe6lHP7+lE+9i5zPH2DTey9zzpxvuOVxO9paObR3G+Vb38dx4CMyGjaRQR0ZwGFJZGf0LPxGziJt+lxGpGR84YwzpdzmNCNpGxYtWsTy5cuZMGECr7/+Oj//+c9P3Od0Onn33XcJDw+npKSEuLi4E/fNnz+fkJAQQkJCuOSSS/jss8+47rrruPPOO6moqGDFihVcf/31bjlEo4XuxaZeeTuHCh4n4tPfYC676YxOyGlva+Xw3m0c3beF1iOFBBzdxbCmfaR0lpIunaQDFQxjX2QuezIuJvWcK0jOGEvyaR9ZKe81YcIEXnvttT7vW7RoEXPmzGHWrFlMnjyZhISEE/c9+eSTTJw4kYceeog777yT9evXI9K1asrxv487fvuWW27hxRdfZPny5Tz77LNuya+F7sX8AwI5POn7zNjyMzave5Upl95IW2srdTXlNNZU0lRXSWt9Fe0NVTgbq6G5hsCGAwxr2ktK5+ETxe00wmG/RKqCMymLvhhHYjZJEy5kxKhJJOhZm2oIufTSS3nggQd4+umnuf3224Gu5XObmpqYNWsWsbGx3HfffSxZsuTE55SVlfHoo4/y2WefER8fz9NPP80zzzxz4vNff/117r//fhobG1m3bh0PP9z128e3vvUtcnNzSUpKYsIE91zQQwvdy029ejFHtj7B6A/+lYYPfkCENBMP9LX0WZvxp8IvjsrgTA7HXIp/YjYx6ZNJyZpEaliknhmmhrzjy+cuWbKEhx9+mODgYDIyMnjssceArlH6/fffz4IF/zyj+e677+aee+7h+IKDjz32GBdddBHXX389ALm5uVx11VUcOHCABx988MRSvImJiWRnZ3Pddde5L78x/Vo48azl5OSY/Px8K1/b12z7cCUtG1+iMygaExyDhA3DPzyWoIh4gqPiCIuOJ2JYImFhkbpOivJYRUVFZGdn247hVkuXLiU8PJwf/vDLF3Jrampi0qRJbNq0iaiovi8v2de/iYhsNMbk9LW/jtB9wMSLF4CugaKU13j33Xf5zne+w913333SMu8PLXSllBogS5cu7XP7ZZddxoEDB9z+9fT3b6WUx7B1CNgT9effQgtdKeURgoODqa6u1lKnq8yrq6sJDg4+o8/TQy5KKY+QmppKaWkplZWVtqN4hODgYFJTz2zumRa6UsojBAQEkJmZaTuGV9NDLkop5SO00JVSykdooSullI+wdqaoiFQCJf389Digyo1xbPGF56HPwTPoc/AMg/Ec0o0xfa3uYa/Qz4aI5J/s1Fdv4gvPQ5+DZ9Dn4BlsPwc95KKUUj5CC10ppXyEtxb6U7YDuIkvPA99Dp5Bn4NnsPocvPIYulJKqS/z1hG6UkqpXrTQlVLKR3hdoYvIXBHZKSLFInKf7Tz9ISL7RWSriBSIiFdctklEnhWRChHZ1mPbMBFZKyK7u/+OsZnxdE7yHJaKyKHu16JARK60mfF0RGSEiLwvIkUisl1EftC93Wtei1M8B695LUQkWEQ+E5HN3c/hP7q3W30dvOoYuog4gF3A5UApkAcsMsYUWg12hkRkP5BjjPGakyhE5GLgGPBnY8zE7m2PAEeNMQ93/+caY4y512bOUznJc1gKHDPG/NpmNleJyHBguDFmk4hEABuB64Bv4SWvxSmew9fxktdCRAQIM8YcE5EA4CPgB8BXsfg6eNsIPRcoNsbsNca0AcuB+ZYzDQnGmA+Bo702zwee6/74Obp+KD3WSZ6DVzHGHDHGbOr+uAEoAlLwotfiFM/Ba5gux7pvBnT/MVh+Hbyt0FOAgz1ul+Jl3wjdDPA3EdkoInfYDnMWEo0xR6DrhxRIsJynv+4SkS3dh2Q89lBFbyKSAUwDPsVLX4tezwG86LUQEYeIFAAVwFpjjPXXwdsKXfrY5j3HjP7pAmPMOcA84M7uQwHKjt8Do4CpwBHgv62mcZGIhAN/AZYYY+pt5+mPPp6DV70WxphOY8xUIBXIFZGJliN5XaGXAiN63E4FDlvK0m/GmMPdf1cAK+k6lOSNyruPhx4/LlphOc8ZM8aUd/9gOoGn8YLXovuY7V+AF40xK7o3e9Vr0ddz8MbXAsAYUwusA+Zi+XXwtkLPA7JEJFNEAoGFwGrLmc6IiIR1vxGEiIQBc4Btp/4sj7Ua+Gb3x98EXreYpV+O//B1W4CHvxbdb8b9ASgyxjza4y6veS1O9hy86bUQkXgRie7+OAS4DNiB5dfBq2a5AHRPZXoMcADPGmN+YTfRmRGRkXSNyqHrEoAvecNzEJGXgdl0LQ9aDvwMWAW8AqQBB4CvGWM89k3HkzyH2XT9im+A/cD3jh8D9UQiciHwD2Ar4Oze/ABdx6C94rU4xXNYhJe8FiIyma43PR10DYxfMcb8XERisfg6eF2hK6WU6pu3HXJRSil1ElroSinlI7TQlVLKR2ihK6WUj9BCV0opH6GFrpRSPkILXSmlfMT/D2khnHrsoXA7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAE2CAYAAACX9WKsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABeCklEQVR4nO3dfVRV953v8c/hUcBnBQXRIj2KihpqJQ/TTB7KUFsmxUnSOqS9V7NsysKYSTpJOrUryWpnctuQO51Z6b246rBqreltpJOZRKa3QuKYsUlzqwQT8yAxIQkmggTxAfEBOXA4948saY7A/u0DZx/O2b5fa7GWnN9vf/fv7LPP52x/7LO3JxAIBAQAAAAAABDF4sZ7AAAAAAAAACZMYAAAAAAAgKjHBAYAAAAAAIh6TGAAAAAAAICoxwQGAAAAAACIekxgAAAAAACAqDemCYz6+nrl5eXJ6/WqsrIyXGMCAFyGvAWAyCBvASB6eQKBQGA0C/r9fi1cuFC7d+9Wdna2CgsLtWPHDi1ZsiTcYwSAKxp5CwCRQd4CQHQb9RkYDQ0N8nq9ys3NVVJSksrKylRbWxvOsQEARN4CQKSQtwAQ3UY9gdHW1qa5c+cO/p6dna22trawDAoAoonpdOLTp0/r1ltv1fLly3X11VfrrbfeCuv6yVsAVwryFgBgJWG0Cw73zROPxzPkserqalVXV0uS3n77bX32s58d7SptG24cobSHaz1XWo1YMspvToVcI1rWE45xRHI9Jh9//LFOnDgx5PEvf/nLwz5uZebMmaqvrx+x3e/3a+PGjUGnE5eWlgadTvzjH/9YBQUFevbZZ3X48GFt3LhRe/bsCWkcVsjb6Mm5aKkRS6IlByO1nishb0eTtZJ78/b111/X1KlTLesmJFgfcpvaJSkuzvrvjgMDA8Ya/f39xj59fX2W7RcvXjTWuHDhgmX7lClTjDWmT59u7BMfH2/ZbmebmN5Lfr/fWMO0XU3bVJJ6e3uNfcKxD5jWY2cfSUpKMvYx7dN2PgtNz8fOa9PT02PsM2HCBMv25ORkYw3T+8K0DwwMDAz7fJ3K21g26gmM7OxsHT16dPD31tZWZWVlDelXXl6u8vJySdLy5cv1H//xH5Z1TW9MU7tkDjNTezStx06NcKzHzvON1H9UwiEcH0amPuH4QLNTJ1pq2KkTju1q56B87dq1wz5+4sQJvfLKK8blP62wsNCy/dOnE0saPJ340wfUTU1N+v73vy9JWrRokY4cOaKOjg7NmjUrpLGMhLwlb0NtjyTyNvw17NQZz7wdTdZK7s3b9PR0rV692rJuenq6Zbud8Zv+E2XnP2qnTp0y9mltbbVsf+edd4w1Xn/9dcv2m2++2VijrKzM2GfatGmW7WfPnjXWML3Xzp07Z6zR0dFh2f7xxx8ba7S0tBj7hGMfMK2ns7PTWGP+/PnGPjNmzLBstzNpZ5oUOH36tLHGm2++aeyTl5dn2Z6Tk2OscfjwYct20z7S3d097ONO5W0sG/VXSAoLC9Xc3KyWlhb5fD7V1NSotLQ0nGMDgJAFAoGQfkzsnE581VVX6ZlnnpH0yQH4hx9+aDwADAV5CyDahJq15C0AjI4TeRvLRn0GRkJCgqqqqrRq1Sr5/X6tX79e+fn54RwbAIQs1NDu7OzUypUrB3//9F/VRqp3+V+/N23apPvuu08FBQVatmyZPve5z9n6y4Jd5C2AaDOaA2TyFgBC5/YJiVCNKfFLSkpUUlISrrEAwJiMZtY5PT1djY2NI7bbOZ148uTJ2rZt2+AY5s+fb+v0ylCQtwCixWj/wkfeAkBoroQzKkI16q+QAEA0CvcpdnZOJ+7q6pLP55Mk/fznP9cNN9ygyZMnO/L8ACAaOHFKM3kLAEPxFZJg4TvnDgCiQLhDe6TTibds2SJJqqio0Ntvv621a9cqPj5eS5Ys0datW8M6BgCINk4cIJO3ADCU2yckQsUEBgBXcSLkhzuduKKiYvDf1113nZqbm8O+XgCIVk4dUJO3ABCMCYxgfIUEAAAAAABEPc7AGMHlV71GbAnHTKWpxsDAQFjGEYmxxlINO9vV6TEgssjb2Ebejk+NcBhL3kbLc4gGfr9fZ8+etezzzjvvWLZfvHjRuJ4VK1ZYtmdkZBhrJCYmGvtcfuHUy02bNs1Yw7R/7Nu3b8zjkKR58+ZZtnd1dRlrmK6hkpKSYqzR19dn2d7b22usYdqHJKmjo8Oy/ZZbbjHWMI3l448/Nta4dB0aKxcuXLBst7Mvmvaj1NRUY43Fixcb+5w6dcqy3c5+1N3dbdluutjw4cOHR2wjb4MxgQHANa6ECxcBwHgjawEgMsjboZjAAOAqhDwAOI+sBYDIIG+DMYEBwFUIeQBwHlkLAJFB3gZjAgOAqxDyAOA8shYAIoO8DcYEBgBXIeQBwHlkLQBEBnkbjAkMAK7BhY4AwHlkLQBEBnk7FBMYAFyFkAcA55G1ABAZ5G2wqJvA8Hg8jtcIxzrCWWe81yFFzxvDzvMNx1ij5flK5rHYGWs4atgRC9s+ml7baOemvI3EcwmXaNlHydvQ28NVw45IjHUsoul1HW8XL17U4cOHLfu0trZatk+aNMm4nkOHDlm2JyYmGmtkZWUZ+0yZMsWyfe7cucYaH374oWV7Z2ensUZtba2xz+c+9znL9v7+fmON1NRUy/akpCRjDZMLFy4Y+5w4ccLY54MPPrBs/9KXvmSs8eabb1q2nz592ljDzv5qyogJEyYYa5j2aTv7/Jw5c4x95s2bZ9k+efJkYw3Te9y0PeLi4ka97JUm6iYwAGAsCHkAcB5ZCwCRQd4GYwIDgGvwPUEAcB5ZCwCRQd4OxQQGAFch5AHAeWQtAEQGeRuMCQwArkLIA4DzyFoAiAzyNtjIVwsBAAAAAACIEpyBAcBVmKUGAOeRtQAQGeRtMCYwALgKIQ8AziNrASAyyNtgfIUEgGtculJzKD8AgNCMJmvJWwAInZN5W19fr7y8PHm9XlVWVg677nvvvVder1fLly/Xq6++alz21KlTKi4u1oIFC1RcXKzTp09Lkk6ePKmbb75ZEydO1D333DPseEpLS7V06VLjuCN+BobH4xnT8nFxkZlzGes47daI1PMx7ch2xmqqYefNEo7takc4xhqOg61YqhGOPnZqDAwMjLmGk8u7yZWUt3aQt84gb0Ov4Ya8JWv/pLe3Vy0tLZZ9enp6LNtvvPFG43peeOEFy/bMzExjjZycHGOfOXPmWLYnJiYaa1x//fWW7QkJ5v+CHDhwwNjnt7/9rWX7okWLjDXOnz9v2W5nX4+Pj7dsT0lJMda49B89K319fZbtzzzzjLGGaaxpaWnGGufOnTP2MUlOTjb2Me0nSUlJxhqmHJSkSZMmWba3tbUZa5g+c02vr9/vH7HNibz1+/3auHGjdu/erezsbBUWFqq0tFRLliwZ7FNXV6fm5mY1Nzdr//792rBhg/bv32+5bGVlpYqKirRp0yZVVlaqsrJSjz/+uCZMmKBHH31Ub731lt56660h43nmmWc0ceJEW2PnDAwArsJfBAHAeZyBAQCR4UTeNjQ0yOv1Kjc3V0lJSSorK1NtbW1Qn9raWq1du1Yej0fXXnuturq61N7ebrlsbW2t1q1bJ0lat26ddu7cKemTibHrr79eEyZMGDKWc+fO6Z//+Z/18MMP29oeTGAAcBUOqAHAeUxgAEBkOJG3bW1tmjt37uDv2dnZQ840GamP1bIdHR2DZ4NlZmbq+PHjxrE88sgjeuCBB5SammrsKzGBAcBlOKAGAOcxgQEAkTGavO3s7NTKlSsHf6qrq4fUvNzlX4MZqY+dZe06ePCg3nvvPd166622l2ECA4BrOHVAbbrI0ZkzZ/TVr35VV111lfLz87Vt27ZwPzUAiBqjyVryFgBCN9q8TU9PV2Nj4+BPeXl5UN3s7GwdPXp08PfW1lZlZWXZ6mO17KxZs9Te3i5Jam9vV0ZGhuXz++Mf/6gDBw4oJydH119/vd59913ddNNNlsswgQHAVcJ9QH3pQkV1dXVqamrSjh071NTUFNRn8+bNWrJkiV5//XXt3btXDzzwgHw+n1NPEQDGnRMTGOQtAAzlRN4WFhaqublZLS0t8vl8qqmpUWlpaVCf0tJSPfnkkwoEAtq3b5+mTJmizMxMy2VLS0u1fft2SdL27du1evVqy3Fs2LBBx44d05EjR/SHP/xBCxcu1N69ey2XifhdSADASeE+TfnTFyqSNHihok9fpdnj8ejs2bMKBAI6d+6cpk+fbusK6wAQq5z4Sgh5CwBDOZG3CQkJqqqq0qpVq+T3+7V+/Xrl5+dry5YtkqSKigqVlJRo165d8nq9Sk1NHTzjbaRlJWnTpk1as2aNtm7dqnnz5unpp58eXGdOTo66u7vl8/m0c+dOPf/880H5bnvsYXj+ABA1wh3yw12oaP/+/UF97rnnHpWWliorK0tnz57Vb37zm4jdshMAxoMTB9TkLQAM5dQ1hEpKSlRSUhL0WEVFxeC/PR6PNm/ebHtZSZoxY4b27Nkz7DJHjhyxHE9OTs6wt1i9XMQnMEwX+AjHfevDUSOa1hMJH374obHP5adxXq6/v99Yw3Rv8U8fuIxk2rRpxj6mgxk7QWDqE44advuY2LnHdSTGEQ01Ql3+0kWOLikvLw/6nqCdCxU999xzKigo0AsvvKD3339fxcXF+vM//3NNnjw5xNGHF3kbncjb0PuQt9FXYzTLujVvZ86cqW9+85uWfZ555pkxr6evr8+y/eTJk8YadrJ0+vTplu1nzpwx1liwYIFl+7vvvmussXDhwjH3ef7554017N75wIrptbHD7/cb++Tk5Fi228mnzs5Oy/aJEycaa/T29hr7mLaJnYlH0/Ox83wTExONfUzXaXjnnXeMNS7deWMkPT09lu1WX4XjIsjBOAMDgGuM5kr3ly5yNBI7Fznatm2bNm3aJI/HI6/Xq/nz5+vw4cO6+uqrQ3sCABADRpO1EnkLAKEabd66GefcAXCV8bjI0bx58wZPl+vo6NA777wz+B1uAHCj8bqoHHkL4Erj1F2fYhVnYACABTsXOXrkkUd05513atmyZQoEAnr88cc1c+bMcR45AMQW8hYAYMIEBgBXcWLW2XSRo6ysLFvfswUAtxivi8qRtwCuNG4/oyJUY5rAyMnJ0aRJkxQfH6+EhATL7zUCQCS4NeTJWwDRxK1ZK5G3AKKLm/N2NMZ8BsZ//dd/ceoegKjh5pAnbwFECzdnrUTeAogebs/bUPEVEgCucSVcuAgAxhtZCwCRQd4ONaa7kHg8Hn3pS1/S5z//eVVXV4drTAAwam69SjN5CyCauPmq+OQtgGji5rwdjTGdgfHyyy8rKytLx48fV3FxsRYtWqQbbrghqE91dfVg+J86dWosq5P0yYdKJGpEaj2RqCFJXV1dlu2///3vjTX+4z/+w7Ldzuubn59v2X7zzTcba9i51/v06dMt2xMSzLv+wMCAZXu4wsFUJxzrsVMjXH3GWmOs63BraJO3sVFDIm8vR94612esNcayDrdmrRR63no8Hn3ta1+zrNnf32/Z3t3dbRxXcnKyZfvUqVONNdLS0ox9ioqKLNufeuopY40XXnjBsj01NdVY47Of/ayxz9NPP23ZnpiYaKzR09Nj2W4nw/x+/5ja7a5nwoQJlu0FBQXGGjNmzLBsP3TokLHGuXPnjH1M4uPjjX1M+7xpe0hSUlKSsc9HH31k2f6FL3zBWOP06dPGPlbOnz8/Ypub83Y0xnQGRlZWliQpIyNDt956qxoaGob0KS8vV2NjoxobG40HOQAwVm6doSZvAUQTN/9FMNS8tTNxAACj5ea8HY1RT2CcP39eZ8+eHfz3888/r6VLl4ZtYAAwGm4MePIWQLRx6wE1eQsg2rg1b0dr1F8h6ejo0K233irpk9PivvGNb+jLX/5y2AYGAKFya2iTtwCiiVuzViJvAUQXN+ftaI16AiM3N1evv/56OMcCAGPmxpAnbwFEGzdmrUTeAog+bs3b0eI2qgBchZAHAOeRtQAQGeRtMCYwALgKIQ8AziNrASAyyNtgTGAAcBVCHgCcR9YCQGSQt8HGdBtVAAAAAACASIj4GRgej8cV64jUeuzMuPX29hr7NDU1Wba/8cYbY66RkGDend59913L9uTkZGONtLQ0Y5+FCxdatk+bNs1YIykpybLdzmsTjhnTaKlhp044tslYxsqVmoORt6HVIG+HIm/Hp4adOuM5VrI2WGJiombPnm3Z5wtf+IJl+4EDB4zrmTFjhmX7Sy+9ZKyRmJho7POb3/zGsv3GG2801li0aJFle1tbm7HGc889Z+xz7tw5y3Y7mR0XZ/333L6+PmONcLwfBgYGjH1M223OnDnGGrm5uZbtEydONNb48MMPjX1OnTpl2d7f32+skZqaatk+depUYw07n5der9eyvaCgYMzref/99y3bX3zxxWEfJ2+H4iskAFyFkAcA55G1ABAZ5G0wJjAAuAohDwDOI2sBIDLI22BMYABwFUIeAJxH1gJAZJC3wZjAAOAqhDwAOI+sBYDIIG+DMYEBwDW40BEAOI+sBYDIIG+HYgIDgKsQ8gDgPLIWACKDvA3GBAYAVyHkAcB5ZC0ARAZ5GyziExgej8ey3XQvZtPy4RhDJPX09Fi2nz592lijvb3d2OfgwYOW7W+++aaxhune0lOmTDHWmDlzpmV7fHy8scYHH3xg7GO6z7mde09Pnz59TOuQ7AWOqU84atgRqXuYOz0OQv5PyNtgFy5csGwnb4cib0fXJxI17ORtOJ6vE8u6jd/v16lTpyz7TJo0ybI9JSXFuJ6MjAzLdjvv148++sjY5+abb7ZsT0pKMtbYsmWLZXtLS4uxht/vN/bx+XyW7XbGaqphh+mzLlzZYvoce+edd4w1ent7LdvtfL7MmjXL2Cc5Odmy3fSekcyfHwkJ5v/KTps2zdjH9N4yfZ5K0owZMyzbTe/xV155ZcQ2p/K2vr5e9913n/x+v+666y5t2rRpyHrvu+8+7dq1S6mpqfrlL3+pFStWWC576tQp/fVf/7WOHDminJwc/eu//qumTZumkydP6mtf+5peeeUV3XnnnaqqqpL0yT799a9/Xe+//77i4+P11a9+VZWVlZbjtj56BYAYc+m7gnZ/AAChCzVryVsAGB0n8tbv92vjxo2qq6tTU1OTduzYoaampqA+dXV1am5uVnNzs6qrq7VhwwbjspWVlSoqKlJzc7OKiooGJyMmTJigRx99VD/5yU+GjOXBBx/U4cOH9dprr+nll19WXV2d5diZwADgGk4dUNfX1ysvL09er3fYWeF//Md/VEFBgQoKCrR06VLFx8fb+ssCAMSi0WQteQsAoXMqbxsaGuT1epWbm6ukpCSVlZWptrY2qE9tba3Wrl0rj8eja6+9Vl1dXWpvb7dctra2VuvWrZMkrVu3Tjt37pQkpaWl6frrr9eECROC1pGamjp41ldSUpJWrFih1tZWy7EzgQHAVcZjhvq73/2uDh48qIMHD+qxxx7TjTfeaDwlHgBi2Xj9RZC8BXClcSJv29raNHfu3MHfs7Oz1dbWZquP1bIdHR3KzMyUJGVmZur48eO2n2dXV5d++9vfqqioyLIfExgAXGU8Zqg/bceOHbrjjjvC+ZQAIOqM118EP428BXAlGE3ednZ2auXKlYM/1dXVQ2pe7vJruYzUx86yoerv79cdd9yhe++9V7m5uZZ9uQsJAFexe5qyXcPNMu/fv3/YvhcuXFB9ff3ghYkAwK3CnbUSeQsAwxlN3qanp6uxsXHE9uzsbB09enTw99bWVmVlZdnq4/P5Rlx21qxZam9vV2Zmptrb240XSL2kvLxcCxYs0He+8x1jX87AAHBFC8cM9SW//e1v9YUvfIHTmQFgGOQtAESHwsJCNTc3q6WlRT6fTzU1NSotLQ3qU1paqieffFKBQED79u3TlClTlJmZablsaWmptm/fLknavn27Vq9ebRzLww8/rDNnzuiJJ56wNXbOwADgKqHOUodjhvqSmpoaTmcGcEUYr78IXkLeArhSOHHGW0JCgqqqqrRq1Sr5/X6tX79e+fn5g7dArqioUElJiXbt2iWv16vU1FRt27bNcllJ2rRpk9asWaOtW7dq3rx5evrppwfXmZOTo+7ubvl8Pu3cuVPPP/+8Jk+erB/96EdatGjR4C1a77nnHt11110jjz3sWwMAxkkoV7q369OzzHPmzFFNTY2eeuqpIf3OnDmj3//+9/o//+f/hHX9ABBtnMhaibwFgMs5lbeSVFJSopKSkqDHKioqBv/t8Xi0efNm28tK0owZM7Rnz55hlzly5Miwj4f6/JjAGAPTxh7pRfq0ffv2WbZ/+OGHxhoff/yxsc+n/6IxnIGBAWONKVOmWLbPmTPHWGP27NmW7RcuXDDWOH36tLGP6fY7CQnmXT8+Pt6y3bQ9JHsXtAlHKJlq2FlHuPqMtcZY1xHukLczQy1Jzz77rL70pS8pLS0trOvHJ8jbYOTtUORt6DXGso7x+ougFH1529vba8yP8+fPW7bbeT/29vZatk+cONFYIyUlxdinvb3dsv3Xv/61sUZycrJl+8WLF4014uLM33RPSkqybLeTt6bsMOWTnRp+v99Yw05W+nw+y3Y7d3wwjWWks54+LTU11dinr6/Pst3OPnDy5EnLdjvvm66uLmOfnJwcy3Y77y3T/mraV632IacmMGIVExgAXMWJkDfNUEvSnXfeqTvvvDPs6waAaDRefxGUyFsAVxYmMIIxgQHAVQh5AHAeWQsAkUHeBmMCA4CrEPIA4DyyFgAig7wNxgQGANdw8kJHAIBPkLUAEBnk7VBMYABwFUIeAJxH1gJAZJC3wZjAAOAqhDwAOI+sBYDIIG+DMYEBwFUIeQBwHlkLAJFB3gaLuQmMcNw73U4NO/eNfu+99yzbX3jhBWONvXv3WrYfPnzYWMPOPZBN92ueMWOGsUZiYqJle3d3t7HGtGnTLNvt3PPbzno6Ozst26dMmWKsYepj517qpm0mmfdHO/viWNcRyTrheH86uTz+hLwdirwdirwNbR2RrONk3pK1f+Lz+XT06FHLPn6/37L9/PnzxvV4PB7Ldjvv19OnTxv7vPjii5btdt4nprHaeb8mJSUZ+/T391u2m7a7FJ73SThq2MnkhATr/7rZ+YwyvX49PT3GGrNnzzb2MX022NkXTfuJnc+GCRMmGPucOHHCsv348ePGGjNnzhzTOKxef/I2WMxNYADASLjQEQA4j6wFgMggb4diAgOAqxDyAOA8shYAIoO8DWY+VwkAAAAAAGCccQYGAFdhlhoAnEfWAkBkkLfBjGdgrF+/XhkZGVq6dOngY6dOnVJxcbEWLFig4uJiWxdhAYBIuPRdQbs/0YS8BRArQs1a8hYARifW8zbcjBMYd955p+rr64Meq6ysVFFRkZqbm1VUVKTKykrHBggAoYjlgCdvAcSKWD+gJm8BxIpYz9twM05g3HDDDZo+fXrQY7W1tVq3bp0kad26ddq5c6cjgwOAUMR6wJO3AGLBaLKWvAWA0Lkhb8NtVNfA6OjoUGZmpiQpMzPT1r1xASAS3Bba5C2AaOS2rJXIWwDRyY15OxaOX8Szurpa1dXVkj75buFY2XkBTX38fr+xxttvv23ss3fv3jG1S9Jrr71m2d7b22usER8fb+zT398/5hqmPomJicYal/+143IzZsww1khOTjb28fl8Y2qXpIGBAct2O/uiqUa4hCPYoqWGG8YwXshba+RtMPJ2dKIlK8c768Z7/ePt03nb3d2tY8eOjaleT0+PsU9CgvVh+8SJE4013nvvPWMfU/7Yee1nz55t2T516lRjjcmTJxv7mJ6PneuY2PmcGiuPx2Ps09fXZ+xjyn47r43ps85OltoZqyn7vV6vscbZs2ct28+cOWOs0dnZaexz/vx5y3Y7x1QLFy60bM/JybFst3rtrvS8vdyobqM6a9Ystbe3S5La29uVkZExYt/y8nI1NjaqsbHReDAFAGPltlPsyFsA0ciNpzSPNm9TUlIiNUQAVyA35u1YjGoCo7S0VNu3b5ckbd++XatXrw7roABgtNwW8OQtgGjkxgNq8hZANHJj3o6FcQLjjjvu0HXXXad33nlH2dnZ2rp1qzZt2qTdu3drwYIF2r17tzZt2hSJsQKApVgPePIWQCxww0XlyFsAscANeRtuxmtg7NixY9jH9+zZE/bBAMBYxXJok7cAYkUsZ61E3gKIHbGet+Hm+EU8ASCSCHkAcB5ZCwCRQd4GYwIDgKsQ8gDgPLIWACKDvA02qot4AkC0cuI7gvX19crLy5PX61VlZeWwffbu3auCggLl5+frxhtvDOdTAoCo49R3sslbAAjGNTCCcQYGANdwIrT9fr82btyo3bt3Kzs7W4WFhSotLdWSJUsG+3R1denuu+9WfX295s2bp+PHj4d1DAAQTZw6QCZvASDYlTAhEaqom8AYGBhwvMa7775rrPH8888b+5gu9NTU1GSs4fP5LNvt7LB+v9/Yx1QnPj7eWMPj8Vi2x8WZT+iZOnWqZfvMmTONNazuy37J+fPnLdsTExONNUzP1w63BU44no+pRrRts4aGBnm9XuXm5kqSysrKVFtbG3RA/dRTT+m2227TvHnzJNnbR6MBeRuMvB2KvB0/kXg+0bbNYjVve3t79d5771n2SU5Otmy/ePGicT2TJk2ybO/v7zfWMOWgHcuWLTP2MWV/fn6+scal/cBKT0+PZfvZs2eNNfr6+izb7XxWmnLdTu7bWY/pNbbznjbl7YQJE4w17GzXtLQ0y/aJEycaayQkWP9X1c7ni2kfkaSTJ09atp8+fXrMNbq6uizbL1y4YFxHuNXX1+u+++6T3+/XXXfdNeTOS4FAQPfdd5927dql1NRU/fKXv9SKFSsslz116pT++q//WkeOHFFOTo7+9V//VdOmTdPJkyf1ta99Ta+88oruvPNOVVVVDa7nwIEDuvPOO9XT06OSkhL99Kc/tdxP+QoJAFcJ9yl2bW1tmjt37uDv2dnZamtrC+rz7rvv6vTp07rpppv0+c9/Xk8++WTYnxcARBMnTmkmbwFgKCfy9tIZb3V1dWpqatKOHTuG/DGorq5Ozc3Nam5uVnV1tTZs2GBctrKyUkVFRWpublZRUdHgVwEnTJigRx99VD/5yU+GjGXDhg2qrq4eXFd9fb3l2JnAAOAqoQZ8Z2enVq5cOfhTXV09pN7lLp8V7u/v14EDB/S73/1Ozz33nB599FFbZx4AQKwazQE1eQsAoXNiAuPTZ7wlJSUNnvH2abW1tVq7dq08Ho+uvfZadXV1qb293XLZ2tparVu3TpK0bt067dy5U9InZ+Rcf/31Q87waW9vV3d3t6677jp5PB6tXbt2cJmRRN1XSABgLEI9JTo9PV2NjY0jtmdnZ+vo0aODv7e2tiorK2tIn5kzZyotLU1paWm64YYb9Prrr2vhwoWhDR4AYsRovn5C3gJA6Jz4ut9wZ7zt37/f2Ketrc1y2Y6ODmVmZkqSMjMzjdcpamtrU3Z29pB1WOEMDACuEu4Z6sLCQjU3N6ulpUU+n081NTUqLS0N6rN69Wq99NJL6u/v14ULF7R//34tXrzYqacIAOPOib8IkrcAMNR4nfE2Uh87y4by3EKtxRkYAFzD7kFyKBISElRVVaVVq1bJ7/dr/fr1ys/P15YtWyRJFRUVWrx4sb785S9r+fLliouL01133aWlS5eGdRwAEC2cyFqJvAWAy402b8N1xttwfXw+34jLzpo1S+3t7crMzFR7e7vxQsvZ2dlqbW21HMflmMAA4CpOHFSXlJSopKQk6LGKioqg37/73e/qu9/9btjXDQDRyImslchbALicE3n76TPe5syZo5qaGj311FNBfUpLS1VVVaWysjLt379fU6ZMUWZmptLT00dctrS0VNu3b9emTZu0fft2rV692nIcmZmZmjRpkvbt26drrrlGTz75pP7mb/7GchkmMAC4ilMH1QCAPyFrASAyxuuMt5KSEu3atUter1epqanatm2b5bKStGnTJq1Zs0Zbt27VvHnz9PTTTw+uMycnR93d3fL5fNq5c6eef/55LVmyRD/72c8Gb6P6la98RV/5ylesxx72rREFTPfZ3bdvn7HG3r17jX1MV73u7e011vD7/ZbtcXHmy5TY+c6R6d7Sdu4/bmcsJklJSZbtdu49bUdKSoplu537RpvuPW2HnXt6h2O7mlxJB5pX0nONBuTtUORtMPJ2/Dg5lmh6nuOtp6dHb7zxhmWfadOmWbbbyYWPP/7Yst3Oey05OdnYp7u727LdlKWSVFRUZNl++e0ah5OTk2PsM3XqVMv2+Ph4Yw1TLtjJJ9PrZ8pjyd52NT0fO+9LUx+fz2esYSdve3p6LNvPnTtnrGHa9nZeG9NnlGT+PLTz/jQdp/T391u2W22v8TrjzePxaPPmzbaXlaQZM2Zoz549wy5z5MiRYR9fuXKl3nrrLZujdukEBoArk1PfywYA/AlZCwCRQd4OxQQGAFch5AHAeWQtAEQGeRuMCQwArkLIA4DzyFoAiAzyNhgTGABchZAHAOeRtQAQGeRtMCYwALgKIQ8AziNrASAyyNtgzl+OGwAAAAAAYIw4AwOAa3ClZgBwHlkLAJFB3g7FBAYAVyHkAcB5ZC0ARAZ5GyzmJjDsvIA+n8+y/fjx48YaJ0+eNPa5cOGCZfvAwICxRjh4PJ4x17Az1osXL465Rlyc9beW7DyXc+fOGfvk5eWNeT3x8fHGPuEQqf3kSkHIhw95OxR5OxR5e2Uia/+kv79fJ06csOxjyrmHH37YuJ4tW7ZYth87dsxY4+zZs8Y+puww5ZMkzZ0717K9q6vLWGPbtm3GPpMmTbJsX7hwobHGxx9/bNluJwdTU1Mt28Ox3SXz+87O+9LUx/S5LkkJCeb/Qvb391u2nz592lgjJSXFsj0cn4WSlJiYOKZx2OnT0tJi2e73+0dsI2+DxdwEBgBYIeQBwHlkLQBEBnkbjAkMAK5CyAOA88haAIgM8jYYExgAXIMLHQGA88haAIgM8nYoJjAAuAohDwDOI2sBIDLI22BMYABwFUIeAJxH1gJAZJC3wZjAAOAqhDwAOI+sBYDIIG+DMYEBwFUIeQBwHlkLAJFB3gaL+ARGJF6AiRMnWrbbuSf0a6+9Zuzz0UcfWbab7n8sSR6Px9gnHMKxHtNrZ+dezKZ7S/f29o65hmR9L2XJ3vawc49rBBvvgOVCR8HI22DkbTDy1hmR2s/GM+vI2mBpaWm65pprLPskJSVZtj/66KPG9Zw5c8ayPS4uzlhjypQpxj7p6eljXs+bb75p2X7ixAljjalTpxr7mLI/LS3NWOPaa6+1bG9paTHWMH2OmV5/yd572pT9dmqYMsrOWO1kdk9Pz5jaJfNnkJ2xJicnG/uYno/pWEcy768ZGRmW7W+99dawj5O3Q7nniAEANP6TKABwJSBrASAyyNtgTGAAcBVCHgCcR9YCQGSQt8GYwADgKoQ8ADiPrAWAyCBvg5m/xAYAAAAAADDOjBMY69evV0ZGhpYuXTr42A9/+EPNmTNHBQUFKigo0K5duxwdJADYdeliR3Z/7Kivr1deXp68Xq8qKyuHtO/du1dTpkwZzMR/+Id/GNXYyVsAsSLUrCVvAWB0nMrbWGX8Csmdd96pe+65R2vXrg16/G//9m/14IMPOjYwAAiVE6Ht9/u1ceNG7d69W9nZ2SosLFRpaamWLFkS1O/P//zP9X//7/8d07rIWwCxwKkDZPIWAIJdCRMSoTKegXHDDTdo+vTpkRgLAIxZuGeoGxoa5PV6lZubq6SkJJWVlam2ttaRsZO3AGKFE38RJG8BYCjOwAg26mtgVFVVafny5Vq/fr1Onz4dzjEBwKiFO+Db2to0d+7cwd+zs7PV1tY2pN8f//hHXXXVVfrKV76iQ4cOhfU5kbcAoo0TB9TkLQAMxQRGsFHdhWTDhg165JFH5PF49Mgjj+iBBx7QL37xi2H7VldXq7q6WpJ06tSp0Y80BMnJyZbtixYtMtaw0+fNN9+0bPf5fMYaJh6PJyJ97NSIi7Oe70pKSjLWMP21IyMjw1jDTp+JEydatpv2Ecn8fO2Eg53tGgnRMo5ICDW0Ozs7tXLlysHfy8vLVV5eblnv8u25YsUKffjhh5o4caJ27dqlv/qrv1Jzc3OIIx8eefsJ8jYYeTvUlZRz0WA0B8huzVufz6eBgQHL2gUFBZbtn564GUlDQ4Nl+4EDB4w1zp49a+xz5swZy/ZJkyYZa0yePNmy3ZQbksLyutqp0dHRYdluZ1+fNm2aZXt3d7exRn9/v7FPOJiej50sTUlJGfM4/H6/sc/58+ct2+18Fk6ZMsXYJzU1dczrSUtLs2w3bbN33nlnxDa3T0iEalQTGLNmzRr897e//W3dcsstI/b99IfT8uXLR7M6ALAt1JBPT09XY2PjiO3Z2dk6evTo4O+tra3KysoK6vPpA7WSkhLdfffdOnHihGbOnBnSWIZD3gKIRqM5oHZr3pr+sw4AY8EERrBRfYWkvb198N/PPvts0BWcAWC8OHGKXWFhoZqbm9XS0iKfz6eamhqVlpYG9fn4448HazU0NGhgYEAzZswIy3MibwFEm9FkLXkLAKFzKm8l812fAoGA7r33Xnm9Xi1fvlyvvvqqcdlTp06puLhYCxYsUHFxcdBX8R577DF5vV7l5eXpueeeG3x8x44dWrZsmZYvX64vf/nLOnHihOW4jWdg3HHHHdq7d69OnDih7Oxs/f3f/7327t2rgwcPyuPxKCcnR//yL/9iKgMAERHuWeqEhARVVVVp1apV8vv9Wr9+vfLz87VlyxZJUkVFhf7t3/5NP/vZz5SQkKCUlBTV1NSM6nR28hZArHDiL4LkLQAM5UTe2rnrU11dnZqbm9Xc3Kz9+/drw4YN2r9/v+WylZWVKioq0qZNm1RZWanKyko9/vjjampqUk1NjQ4dOqRjx47pL/7iL/Tuu+8qEAjovvvuU1NTk2bOnKm/+7u/U1VVlX74wx+OOHbjBMaOHTuGPPatb31rdFsKABzmRMiXlJSopKQk6LGKiorBf99zzz265557xrwe8hZArHDqlGbyFgCCOZG3n77rk6TBuz59egKjtrZWa9eulcfj0bXXXquuri61t7fryJEjIy5bW1urvXv3SpLWrVunm266SY8//rhqa2tVVlam5ORkzZ8/X16vVw0NDVq5cqUCgYDOnz+vGTNmqLu7W16v13Lso74LCQBEI67SDADO46r4ABAZo8nbSxdNvvRz6aLDl9i569NIfayW7ejoUGZmpiQpMzNTx48ft6yVmJion/3sZ1q2bJmysrLU1NRknEwe1UU8ASBacZAMAM4jawEgMpy4aLKduz6N1MfOsnbX19fXp5/97Gd67bXXlJubq7/5m7/RY489pocffnjEWpyBAcA1+IsgADjPyYvKAQD+xKm8tXPXp5H6WC07a9aswQsit7e3D96OfaRlDh48KEn67Gc/K4/HozVr1uj//b//Zzn2mDsDw3Sfbcl8H2U79y62c49r073t7VxUytQnHDXs9ImLM89lmW5RNmfOHGON9PR0y/bs7GxjDTt9TK8fB1Ljx7Ttx/ra8NqGD3nrTB/yNraN5oKR48XJvHXb6zoWqamp+tznPmfZp76+3rLdTt6a3vfvvPOOsUZPT4+xT3x8vGX79ddfb6zx9ttvW7abMl365KKuJrNnz7Zsnz9/vrHGiy++aNk+ceJEYw3TZ6Hf7zfWsPN8TXXsvC9Nn0F29kU7TJ8N3d3dY16Hne3a29tr7GPa5+18bjvJibz99F2f5syZo5qaGj311FNBfUpLS1VVVaWysjLt379fU6ZMUWZmptLT00dctrS0VNu3b9emTZu0fft2rV69evDxb3zjG7r//vt17NgxNTc36+qrr1ZHR4eamprU2dmp9PR07d69W4sXL7Yce8xNYAAAAAAAgNGxc9enkpIS7dq1S16vV6mpqdq2bZvlspK0adMmrVmzRlu3btW8efP09NNPS5Ly8/O1Zs0aLVmyRAkJCdq8ebPi4+OVlZWlH/zgB7rhhhuUmJioz3zmM/rlL39pPXbnNgsARB5/FQQA55G1ABAZ43XXJ4/Ho82bN9teVpJmzJihPXv2DLvMQw89pIceemjI4xUVFUHrNWECA4CrcFANAM4jawEgMsjbYExgAHAVQh4AnEfWAkBkkLfBmMAA4Bpc6R4AnEfWAkBkkLdDMYEBwFUIeQBwHlkLAJFB3gZjAgOAqxDyAOA8shYAIoO8DcYEBgBXIeQBwHlkLQBEBnkbLKITGAMDAzp79qxln/j4+DG1S9LFixct2zs6Oow1Tp48aezj8/ks2z0ej7FGpJjGMmnSJGON+fPnW7YvXLjQWGPy5MmW7bNnzzbWsNPH7/dbtvf29hprmF7fpKQkYw07IrGf2Ak+O30GBgbG1C6ZXxtTuwkh/wnydvyQt8FiKW/t7EemjLlS8pas/ZNTp06ppqbGsk9KSoplu50cNL2X7LyedvrExcWNqV2SLly4YNl+4sQJY42EBPN/U0x9jhw5YqyRlpZm2W7KJ8m8XcPxng8X02e7nXGcP3/e2Me0z5s+oyTzPm9nrHZy3bSec+fOGWv09/dbtvf09Fi2W+1D5G0wzsAA4Bpc6AgAnEfWAkBkkLdDMYEBwFUIeQBwHlkLAJFB3gZjAgOAqxDyAOA8shYAIoO8DcYEBgBXIeQBwHlkLQBEBnkbjAkMAK5CyAOA88haAIgM8jYYExgAXIMLHQGA88haAIgM8nYo832QAAAAAAAAxhlnYABwFWapAcB5ZC0ARAZ5GyyiExh9fX06fvy4ZZ/4+HjL9oQE85D9fr9l++HDh401PvjgA2Ofnp4eY59o4fF4LNvT09ONNebNm2fZ/pnPfMZYIy0tbczjSElJMfYxvTZnz5411ujo6DD2MUlKSjL2iYuzPhHK9J6QzME2MDBgrGF630hSf3+/ZbvP5zPWMPXp6+sz1rBCyH+CvB0/5G0w8nYoN+QtWfsnAwMDOn/+vGUf0/vA9HpL0kcffWQch4npPSBJkydPtmw/efKksYYpB8+dO2esMW3aNGMf0+fUiRMnjDVM28TO+9X0XgrH56lk3q52Xl/TWO3koJ3PZFOfSZMmGWskJycb+5icOXPG2Me0P06cONFYIzU11bLdtF2t8pq8DcZXSAC4yqXvCtr9saO+vl55eXnyer2qrKwcsd8rr7yi+Ph4/du//Vu4ng4ARKVQs5a8BYDRcSpvYxUTGABcJdwB7/f7tXHjRtXV1ampqUk7duxQU1PTsP2+973vadWqVU48LQCIKk4cUJO3ADAUExjBmMAA4BpOBHxDQ4O8Xq9yc3OVlJSksrIy1dbWDun3v//3/9btt9+ujIwMJ54aAESN0WQteQsAoXMqb2MZExgAXCXcAd/W1qa5c+cO/p6dna22trYhfZ599llVVFSE/fkAQDRy4oCavAWAoZjACMZdSAC4Sqih3dnZqZUrVw7+Xl5ervLycst6l19E6zvf+Y4ef/xxWxe+AgA3GM0BMnkLAKFz+4REqJjAAOAqoYZ8enq6GhsbR2zPzs7W0aNHB39vbW1VVlZWUJ/GxkaVlZVJ+uSK57t27VJCQoL+6q/+KqSxAECsGM0BNXkLAKFjAiMYExgAXCXcIV9YWKjm5ma1tLRozpw5qqmp0VNPPRXUp6WlZfDfd955p2655RYOpgG4mhMH1OQtAAzFBEawiE5g9Pf3q7Oz07KPnfsXmxw/ftyy/Y033jDWeO+994x9TPdXt7Ozme7nHC6m+09PnTrVWGP27NmW7Xaer+n1Ddcb1FTHzv2rL//e7eXs3MN8woQJxj6m1yYpKclYIyUlZcw17Nx/3LTd7GyT8+fPW7b39vYaa4zEie/9JSQkqKqqSqtWrZLf79f69euVn5+vLVu2SFLUfg+bvB3KTXlrB3k7VCzl7cWLFy3bz549a6zhVN469R3rWM3bQCBgzCjTPpyWlmZcj+m9ZOc1sdPH9F768MMPx7weO58/dr4mdPr0acv2sRxTXDJ58mRjn/7+fst2O9v9woULxj6mzzE7n3OmbW8nn+y8fqax2Mmw1NRUy/Y5c+YYa9j5/DDtR2fOnDHWMOVtYmKiZftI+9CVcE2LUHEGBgBXcSLkS0pKVFJSEvTYSAfSv/zlL8O+fgCINk4dUJO3ABCMCYxgTGAAcBVCHgCcR9YCQGSQt8G4jSoAV+E2UwDgPG7rBwCR4VTe1tfXKy8vT16vV5WVlcOu995775XX69Xy5cv16quvGpc9deqUiouLtWDBAhUXFwd9Peexxx6T1+tVXl6ennvuucHHfT6fysvLtXDhQi1atEj//u//bjluJjAAAAAAALhC+P1+bdy4UXV1dWpqatKOHTvU1NQU1Keurk7Nzc1qbm5WdXW1NmzYYFy2srJSRUVFam5uVlFR0eDkRlNTk2pqanTo0CHV19fr7rvvHrzeyo9+9CNlZGTo3XffVVNTk2688UbLsRsnMI4ePaqbb75ZixcvVn5+vn76059Ksp5dAYDxEst/ESRvAcSKWD8Dg7wFECucyNuGhgZ5vV7l5uYqKSlJZWVlqq2tDepTW1urtWvXyuPx6Nprr1VXV5fa29stl62trdW6deskSevWrdPOnTsHHy8rK1NycrLmz58vr9erhoYGSdIvfvELff/735f0yQViZ86caTl24wRGQkKC/umf/klvv/229u3bp82bN6upqWnE2RUAGC+xfkBN3gKIBaPJWvIWAELnVN62tbVp7ty5g79nZ2cPuSPYSH2slu3o6FBmZqYkKTMzc/BudSMt09XVJUl65JFHtGLFCn39619XR0eH5diNExiZmZlasWKFJGnSpElavHix2traRpxdAYDxFMsH1OQtgFgR6xMY5C2AWDGavO3s7NTKlSsHf6qrq4fUvNzlt74dqY+dZYd7DsMt09/fr9bWVn3hC1/Qq6++quuuu04PPvigZa2Q7kJy5MgRvfbaa7rmmmtGnF0BgPEUbQfJo0XeAohmbslaibwFEN1Gk7fp6elqbGwcsT07O1tHjx4d/L21tVVZWVm2+vh8vhGXnTVrltrb25WZman29nZlZGRY1poxY4ZSU1N16623SpK+/vWva+vWrZbPzfYExrlz53T77bfriSee0OTJk+0upurq6sEZnzNnzujcuXOW/S9dzGMkFy5cMK7z3XfftWx/++23jTUunc5ixbQzRdOH+6RJkyzb09LSjDVMM2sDAwPGGv39/Zbtvb29xhpxceZrz/p8Psv2np4eY40TJ05Ytp88edJYY8KECcY+CQnWb8PExERjjSlTpli223nPml5fScb3b3d3t7HG2bNnx7QOk2h6340WeTsUeRvM9NpJ5O1wyNvQ1mElmt5zYxGOvB0YGDC+D0zvJTvvadM2t3McZuc9bXqftLe3G2tMnz7dst3OWE3b1E4fUw7aGUt8fLyxhinX7Wz3cGwTO9lien3tjNUOU53U1FRjjYkTJ1q2m66VIElz5swx9uns7LRst3OcYprsNB1TWb3+TuRtYWGhmpub1dLSojlz5qimpkZPPfVUUJ/S0lJVVVWprKxM+/fv15QpU5SZman09PQRly0tLdX27du1adMmbd++XatXrx58/Bvf+Ibuv/9+HTt2TM3Nzbr66qvl8Xj01a9+VXv37tUXv/hF7dmzR0uWLLEcu60JjL6+Pt1+++365je/qdtuu03SyLMrlysvL1d5ebkkKScnx87qAGDUYv2gmrwFEAtiPWul8OVtuP7DBwDDcSJvExISVFVVpVWrVsnv92v9+vXKz8/Xli1bJEkVFRUqKSnRrl275PV6lZqaqm3btlkuK0mbNm3SmjVrtHXrVs2bN09PP/20JCk/P19r1qzRkiVLlJCQoM2bNw9ODj7++OP67//9v+s73/mO0tPTB9cz4thNTy4QCOhb3/qWFi9erPvvv3/w8ZFmVwBgvETj96xDQd4CiAWxnrUSeQsgNjiZtyUlJSopKQl6rKKiYvDfHo9Hmzdvtr2sJM2YMUN79uwZdpmHHnpIDz300JDHP/OZz+jFF1+0PW7jBMbLL7+sX/3qV1q2bJkKCgokST/+8Y9HnF0BgPEUywfV5C2AWBHLWSuRtwBiR6znbbgZJzCuv/76ETfaSLMrADBeYjnkyVsAsSKWs1YibwHEjljP23AL6S4kABDtCHkAcB5ZCwCRQd4GYwIDgKsQ8gDgPLIWACKDvA3GBAYA13DDheUAINqRtQAQGeTtUBGdwIiLizPepz0c9zc23Rvdzr187dx7OpZ2JtO9ze3c+zwlJcWy3c59sk3s3K/74sWLYaljYroft2k/k6Senh5jH9N2s7NdTWOx8/qa7gsuSb29vZbt58+fN9Yw9bHz+lqJpfelk8jb8UPeho68HSra8zaW3pNOCwQC8vv9xj5WTMvbqWHndq529j3TWKZPn26skZiYaOxjMmnSJGOfzs5Oy3Y7+2k4boM71tdfsvd8z5w5Y9luJwdNz9f0+SNJSUlJxj5paWmW7cnJycYaplvCL1iwwFjDzr5o2o9Onz5trPH2229btp84ccKy3eq1I2+DcQYGAFch5AHAeWQtAEQGeRts7FOOAAAAAAAADuMMDACuwiw1ADiPrAWAyCBvgzGBAcBVCHkAcB5ZCwCRQd4GYwIDgGtwpWYAcB5ZCwCRQd4OxQQGAFch5AHAeWQtAEQGeRuMCQwArkLIA4DzyFoAiAzyNhh3IQHgKpdOtbP7Y0d9fb3y8vLk9XpVWVk5pL22tlbLly9XQUGBVq5cqT/84Q/hfloAEFVCzVryFgBGx6m8jVURPQMjLi5OKSkpln08Ho9lu50XxOfzWbYPDAwYa9hZT7TsHAkJ5pdx0qRJlu1Tpkwx1jBtt/j4eGMN01jt1EhMTDT2MUlOTjb26e/vt2y38/rbGavpOdupYdqufr/fWCMc7ws7Y01NTbVsT0pKMtYYiROh7ff7tXHjRu3evVvZ2dkqLCxUaWmplixZMtinqKhIpaWl8ng8euONN7RmzRodPnw4rOMIFXnrDPI2dNGUt6ZtYuf1JW+d+052rOatZO81tWJnvzFtc1OmS598Npj09PSMaRySed+z4/Tp08Y+pu0Wjm0ya9YsY40LFy5Ytq9cudJYo6ury9inqanJsv3ixYvGGqZtZspjKTyfDXbyx/T5kZWVZaxh5zPXtM9/8MEHxhqmfcC03Ud6X10JExKh4gwMAK4S7hnqhoYGeb1e5ebmKikpSWVlZaqtrQ3qM3HixMGDpPPnz9s6YAKAWObEXwTJWwAYijMwgjGBAcBVwh3wbW1tmjt37uDv2dnZamtrG9Lv2Wef1aJFi/SXf/mX+sUvfhHW5wQA0caJA2ryFgCGYgIjGBMYAFwl1IDv7OzUypUrB3+qq6uH1LvccH/xu/XWW3X48GHt3LlTjzzyiGPPDwCiwWgOqMlbAAgdExjBuAsJAFcJNbTT09PV2Ng4Ynt2draOHj06+Htra6vldy5vuOEGvf/++zpx4oRmzpwZ0lgAIFaM5gCZvAWA0Ll9QiJUnIEBwDWcmKEuLCxUc3OzWlpa5PP5VFNTo9LS0qA+77333mCtV199VT6fTzNmzHDkOQLAeBtN1pK3ABA6p/I2lnEGBgBXCXdoJyQkqKqqSqtWrZLf79f69euVn5+vLVu2SJIqKir07//+73ryySeVmJiolJQU/eY3v+HCcgBczYkDZPIWAIZy+4REqJjAAACDkpISlZSUBD1WUVEx+O/vfe97+t73vhfpYQGA65C3AAArTGAAcBVmqQHAeWQtAEQGeRssohMYHo9HCQnWq0xOTrZsj4szX7ZjYGDAst00Brvr8fv9lu12TmkMx2mPaWlpxj6TJk0aU7skxcfHj6ndTh87293O62dieu3s9DHtZ5K919f0fBITE4017Gx7k3CEo51xJCUljbmGFUL+E+Tt6PqYhCNvJ0+ebKxB3gYjb0c3DifzlqwNFontYdrH7bxPfD6fsY/pufT19RlrTJgwwbI9NTXVWOPChQvGPqYcs/O6mN73H330kbGG6bV56aWXjDXsXIjWNFY772nTNrOzj9j5/DC9xikpKcYaZ86csWx/9913jTVyc3ONfT744IMxr+f06dOW7abPH6t9lbwNxhkYAFyFkAcA55G1ABAZ5G0wJjAAuMaVcOVlABhvZC0ARAZ5OxS3UQXgKtxmCgCcx239ACAynMrb+vp65eXlyev1qrKyctj13nvvvfJ6vVq+fLleffVV47KnTp1ScXGxFixYoOLi4qCv1jz22GPyer3Ky8vTc889N2R9paWlWrp0qXHcTGAAcBUOqAHAeUxgAEBkOJG3fr9fGzduVF1dnZqamrRjxw41NTUF9amrq1Nzc7Oam5tVXV2tDRs2GJetrKxUUVGRmpubVVRUNDi50dTUpJqaGh06dEj19fW6++67g65/9cwzz2jixIm2tgcTGABchQNqAHAeExgAEBlO5G1DQ4O8Xq9yc3OVlJSksrIy1dbWBvWpra3V2rVr5fF4dO2116qrq0vt7e2Wy9bW1mrdunWSpHXr1mnnzp2Dj5eVlSk5OVnz58+X1+tVQ0ODJOncuXP653/+Zz388MO2tgcTGABchQNqAHAeExgAEBlO5G1bW5vmzp07+Ht2drba2tps9bFatqOjQ5mZmZKkzMxMHT9+3Li+Rx55RA888ICtuxNJTGAAcBEOqAHAeaPJWvIWAEI32rzt7OzUypUrB3+qq6uH1L3c5bcDHqmPnWWHex7DLXPw4EG99957uvXWWy2X/7SI34XEdN9g0z1yJ02aZFxHRkaGZfu0adOMNTo6Oox97Nzb3sS0PezcZ9nONjF9p8jOjJfp3tJ27j1t6mN6/SV728RUJykpyVjDdLBl9/tlJqaxhmubREI4xmFnP7LCQfKfjDVvJ0+ebFwHeTuUKW9TUlKMNcjb0Nol8nY0xpK3ZG3khWOb26nR399v2T4wMGCscfLkScv2xMREYw07+7jp+HXq1KnGGl1dXZbt4cjsBQsWGGvs37/f2Gf69OmW7Xa2a19fn2W7nc9tO30u/2v+5ey8NgcOHLBsv/yaDcOxcyxz9OhRy/acnBxjDdNnv+l9deLEiRHbRvPeT09PV2Nj44jt2dnZQc+7tbVVWVlZtvr4fL4Rl501a5ba29uVmZmp9vb2wePEkWr98Y9/1IEDB5STk6P+/n4dP35cN910k/bu3Tvi2KPjUxgAwoS/CAKA8zgDAwAiw4m8LSwsVHNzs1paWuTz+VRTU6PS0tKgPqWlpXryyScVCAS0b98+TZkyRZmZmZbLlpaWavv27ZKk7du3a/Xq1YOP19TUqLe3Vy0tLWpubtbVV1+tDRs26NixYzpy5Ij+8Ic/aOHChZaTF9I4nIEBAE7iIBkAnEfWAkBkOJG3CQkJqqqq0qpVq+T3+7V+/Xrl5+dry5YtkqSKigqVlJRo165d8nq9Sk1N1bZt2yyXlaRNmzZpzZo12rp1q+bNm6enn35akpSfn681a9ZoyZIlSkhI0ObNm0d9FiATGABchYNqAHAeWQsAkeFU3paUlKikpCTosYqKisF/ezwebd682faykjRjxgzt2bNn2GUeeughPfTQQyOOJycnR2+99ZZx3ExgAHANTlMGAOeRtQAQGeTtUMZrYBw9elQ333yzFi9erPz8fP30pz+VJP3whz/UnDlzVFBQoIKCAu3atcvxwQKAm5G3ABAZ5C0AxCbjGRgJCQn6p3/6J61YsUJnz57V5z//eRUXF0uS/vZv/1YPPvig44MEALtieZaavAUQK2I5ayXyFkDsiPW8DTfjBEZmZqYyMzMlfXL7uMWLFxtviwMA4yWWQ568BRArYjlrJfIWQOyI9bwNt5Buo3rkyBG99tpruuaaayRJVVVVWr58udavX6/Tp087MkAACIVbbutH3gKIZm66jSp5CyCauSlvw8H2RTzPnTun22+/XU888YQmT56sDRs26JFHHpHH49EjjzyiBx54QL/4xS+GLFddXa3q6mpJ0pkzZzQwMGC5Ho/HY9melpZmHGtGRsaY2iXp2LFjxj79/f2W7X6/31gjLs56DikhwfwSTZo0ydhn4sSJlu12bmNj6mOnhun52Klh2kck83ZNTEw01jC9+e2Eg2kfkczPx/Rc7NSwIxxhF47XZqzjcENoR0vepqamGsdK3g5F3gYjb4dyQ966IWul8ORtpJheT1PmS/ZeN9N67OStKSvPnz9vrJGcnGzsY8q5v/iLvzDWeOmllyzb7Yw1KSnJsv3EiRPGGlOnTjX2uXDhgmW7nW128eLFMa1Dsvd5adrX7Hz2z5w507L9lltuMdZ47733jH1mz55t2W56fSXz++Ls2bOW7V1dXSO2uSVvw8XWGRh9fX26/fbb9c1vflO33XabJGnWrFmKj49XXFycvv3tb6uhoWHYZcvLy9XY2KjGxkZNmTIlfCMHgMu4YYaavAUQ7UaTtW7OWwBwilvyNpyMExiBQEDf+ta3tHjxYt1///2Dj7e3tw/++9lnn9XSpUudGSEAhCCWA568BRArYv2AmrwFECtiPW/DzXj+z8svv6xf/epXWrZsmQoKCiRJP/7xj7Vjxw4dPHhQHo9HOTk5+pd/+RenxwoARrEc2uQtgFgRy1krkbcAYkes5224GScwrr/++mE3WklJiSMDAoCxiOWQJ28BxIpYzlqJvAUQO2I9b8MtpLuQAEC0c+IUu/r6euXl5cnr9aqysnJI+69//WstX75cy5cv15/92Z/p9ddfD/fTAoCo4tQpzeQtAATjKyTBbN+FBACi3WhC23Qlf7/fr40bN2r37t3Kzs5WYWGhSktLtWTJksE+8+fP1+9//3tNmzZNdXV1Ki8v1/79+0f1HAAg2o32AJm8BYDQOJW3sYwzMAC4SrhnqBsaGuT1epWbm6ukpCSVlZWptrY2qM+f/dmfadq0aZKka6+9Vq2trY48NwCIFk78RZC8BYChOAMjWMTPwDBtUFO7ndmktLQ0y3Y7txe006e3t9eyvaenx1jD9HxSU1ONNeyMdfLkyZbt8fHxxhqmPnbuCW26t7id1zccb0o76zE938TExDGPQzI/H9M2C8c6pPDM1EZqPWMdw6d1dnZq5cqVg7+Xl5ervLx88Pe2tjbNnTt38Pfs7GzLv/Zt3bpVX/nKV0Iag1PI22DkbTDydijyNrzrv5yb83asIvUfDjv7eH9/v2W7nQzr6+uzbE9OTjbWsLOeM2fOWLbv27fPWMP0fDs6Oow10tPTjX1MkpKSjH1Onjxp2T5hwgRjDdO2N33eStKFCxeMfUz7mp19YPbs2Zbtb7/9trHG2bNnjX1mzpxp2e7z+Yw1TMZybOD2CYlQ8RUSAK4Sasinp6ersbExpHoj/afgv/7rv7R161b94Q9/CGkMABBrRnNATd4CQOiYwAjGBAYAWMjOztbRo0cHf29tbVVWVtaQfm+88Ybuuusu1dXVacaMGZEcIgC4AnkLADDhGhgAXCXc3xEsLCxUc3OzWlpa5PP5VFNTo9LS0qA+H330kW677Tb96le/0sKFC516agAQNZz4TjZ5CwBDcQ2MYJyBAcA1nAjthIQEVVVVadWqVfL7/Vq/fr3y8/O1ZcsWSVJFRYX+4R/+QSdPntTdd989uIzVadIAEMucOkAmbwEg2JUwIREqJjAAuIoTIV9SUqKSkpKgxyoqKgb//fOf/1w///nPw75eAIhWTh1Qk7cAEIwJjGBMYABwFUIeAJxH1gJAZJC3wZjAAOAqhDwAOI+sBYDIIG+DMYEBwFUIeQBwHlkLAJFB3gaL6ARGIBDQwMCAZZ+xtksj3zP8kuTkZGONtLQ0Y5+pU6daticmJhpr9Pf3W7ZPnjzZWGPKlCnGPqY6CQnmXcHUJz4+3ljDTh+TSL2J4+Ksb9Jj57nYGaudfToSwrFdTe89O+uxU8OqNiH/CfJ2KPI2dOStM2I9b8na2GXndTP1sfM+unDhgmV7amqqsYad9Zgy+fjx48YapnyZPn26sca8efMs2+1s946ODmOflJQUy/aLFy8aa5je93bG2tPTY+xj+sw15b4knT592rLdToY9/PDDxj6//vWvLdsPHTpkrJGXl2fZbjo2GGl7kLdDcQYGAFch5AHAeWQtAEQGeRuMCQwArkLIA4DzyFoAiAzyNhgTGABchZAHAOeRtQAQGeRtMPOXjwAghlz6rqDdHwBA6ELNWvIWAEbHqbytr69XXl6evF6vKisrh13vvffeK6/Xq+XLl+vVV181Lnvq1CkVFxdrwYIFKi4uDrqOyWOPPSav16u8vDw999xzkj65Xs5f/uVfatGiRcrPz9emTZuM42YCA4BrcEANAM4bTdaStwAQOqfy1u/3a+PGjaqrq1NTU5N27NihpqamoD51dXVqbm5Wc3OzqqurtWHDBuOylZWVKioqUnNzs4qKigYnN5qamlRTU6NDhw6pvr5ed999t/x+vyTpwQcf1OHDh/Xaa6/p5ZdfVl1dneXYmcAA4CocUAOA85jAAIDIcCJvGxoa5PV6lZubq6SkJJWVlam2tjaoT21trdauXSuPx6Nrr71WXV1dam9vt1y2trZW69atkyStW7dOO3fuHHy8rKxMycnJmj9/vrxerxoaGpSamqqbb75ZkpSUlKQVK1aotbXVcuxMYABwFQ6oAcB5TGAAQGQ4kbdtbW2aO3fu4O/Z2dlqa2uz1cdq2Y6ODmVmZkqSMjMzB29lbGd9XV1d+u1vf6uioiLLsXMRTwAAAAAAXKKzs1MrV64c/L28vFzl5eWDvw83yeHxeIJ+H6mPnWUvZ1qmv79fd9xxh+69917l5uZa1or4BMbAwIBl+6Xvwoy23U6fuDjziScTJkww9pk8efKYa1y8eNGyferUqWMeh50+qampxhppaWmW7cnJycYaCQnWu5xp55ei50q8dsYaHx8/5vXYeb7h2CZ2nk8kxuGGMUQL8jYYeRuMvB2KvI2d9V+JwrHNI1UjHOvp7e019unr67Nsnz17trFGd3e3ZXtOTo6xxqxZsyzbfT6fscbJkyeNfZKSkizbz507Z6xhOjaww87xgcmFCxeMfUyf23Y+C7///e8b+0ybNs2yfc6cOcYaps9c03vCqn0076f09HQ1NjaO2J6dna2jR48O/t7a2qqsrCxbfXw+34jLzpo1S+3t7crMzFR7e7syMjJsra+8vFwLFizQd77zHeNz4yskAFyFU5oBwHl8hQQAIsOJvC0sLFRzc7NaWlrk8/lUU1Oj0tLSoD6lpaV68sknFQgEtG/fPk2ZMkWZmZmWy5aWlmr79u2SpO3bt2v16tWDj9fU1Ki3t1ctLS1qbm7W1VdfLUl6+OGHdebMGT3xxBO2tgdfIQHgGhwkA4DzyFoAiAyn8jYhIUFVVVVatWqV/H6/1q9fr/z8fG3ZskWSVFFRoZKSEu3atUter1epqanatm2b5bKStGnTJq1Zs0Zbt27VvHnz9PTTT0uS8vPztWbNGi1ZskQJCQnavHmz4uPj1draqh/96EdatGiRVqxYIUm65557dNddd4089rBvDQAYRxxUA4DzyFoAiAyn8rakpEQlJSVBj1VUVAz+2+PxaPPmzbaXlaQZM2Zoz549wy7z0EMP6aGHHgp6LDs7O+TnxwQGAFfhoBoAnEfWAkBkkLfBmMAA4CqEPAA4j6wFgMggb4MxgQHAVQh5AHAeWQsAkUHeBmMCA4BrcGE5AHAeWQsAkUHeDhXRCYxAIGC897DpBbJz72JTn7g4891jTfdZlsz3++3v7zfWMN2/Oi0tzVhjypQpxj4pKSljXs/EiRMt2+1s12gRjiDweDwR6xMJkQpH0/Md6zgI+U+Qt0ORt+ODvB3KDXlL1l657LyPTPtHT0+PsYbf7zf2iY+Pt2zv7u421pg0aZJle05OjrFGamqqZXtnZ6exhukzSjI/XzufL6b1+Hw+Yw0773/TfjJ16lRjDdPzaW9vN9aYMGGCsc+ZM2fGXMO0D5j2eattSt4G4wwMAK5CyAOA88haAIgM8jZY7PwJBwBsuHSqnd0fO+rr65WXlyev16vKysoh7YcPH9Z1112n5ORk/eQnPwn3UwKAqBNq1pK3ADA6TuVtrOIMDACuEu7Q9vv92rhxo3bv3q3s7GwVFhaqtLRUS5YsGewzffp0/a//9b+0c+fOsK4bAKKVEwfI5C0ADOX2CYlQcQYGANdwYoa6oaFBXq9Xubm5SkpKUllZmWpra4P6ZGRkqLCwUImJiU49NQCIGqPJWvIWAELnVN7GMuMExsWLF3X11VfrqquuUn5+vn7wgx9Ikk6dOqXi4mItWLBAxcXFOn36tOODBQCTcAd8W1ub5s6dO/h7dna22traHBk7eQsgVjhxQE3eAsBQTGAEM05gJCcn64UXXtDrr7+ugwcPqr6+Xvv27VNlZaWKiorU3NysoqKiYb+nCADRrrOzUytXrhz8qa6uDmof7kPAqTsZkLcA3Iy8BQCMlfEaGB6PZ/B2bn19ferr65PH41Ftba327t0rSVq3bp1uuukmPf74444OFgBMQp11Tk9PV2Nj44jt2dnZOnr06ODvra2tysrKGvX4rJC3AGLFaP7CR94CQOjcfkZFqGxdA8Pv96ugoEAZGRkqLi7WNddco46ODmVmZkqSMjMzdfz4cUcHCgB2hPsUu8LCQjU3N6ulpUU+n081NTUqLS11bPzkLYBY4MQpzeQtAAzFV0iC2boLSXx8vA4ePKiuri7deuuteuutt2yvoLq6evAUwe7ubg0MDFj2H2u73T4mqampEVmP3++3bE9KSjLWsNPHdLErOxfDio+Pt2yPpjdLNI1lrOycPnulPV8r4d4WCQkJqqqq0qpVq+T3+7V+/Xrl5+dry5YtkqSKigp9/PHHWrlypbq7uxUXF6cnnnhCTU1Nmjx5csjrI2+dWw9564xoGstYkbf2ObEdYjVvrzTheO1NeSxJcXHmv7OastLOPn727FnL9rq6OmON6dOnW7YnJycba9j5nDNt+/7+fmMNkwkTJhj79Pb2GvuYxmLn+Zpe37S0NGONcFzHxk6+mLabaX+22h5u+twJh5Buozp16lTddNNNqq+v16xZs9Te3q7MzEy1t7crIyNj2GXKy8tVXl4uSfrsZz879hEDwAicmnUuKSlRSUlJ0GMVFRWD/549e7ZaW1vDuk7yFkC0cvIvfLGYt05dpwMAroQzKkJlnNrs7OxUV1eXJKmnp0f/+Z//qUWLFqm0tFTbt2+XJG3fvl2rV692dKAAYEcsn2JH3gKIFbF+SjN5CyBWxHrehpvxDIz29natW7dOfr9fAwMDWrNmjW655RZdd911WrNmjbZu3ap58+bp6aefjsR4AcBSLIc2eQsgVsRy1krkLYDYEet5G27GCYzly5frtddeG/L4jBkztGfPHkcGBQCjFcshT94CiBWxnLUSeQsgdsR63oZbSNfAAIBoR8gDgPPIWgCIDPI2GBMYAFzjSvjeHwCMN7IWACKDvB2KCQwArkLIA4DzyFoAiAzyNljEJzBML0A4XiDT7azs3CPZdN9h6ZP7lY+V6R7IdtaRlJRk7GPaJnZqjHUdEm/A4YzlvtCXhOMWbpF474VrPeNZP5aQt8HIW5C34cP+hbEwvRcle58NKSkplu3Tpk0z1pgzZ45l+8mTJ401PvzwQ8v2CxcuGGvYeb49PT2W7XYyzPQZZOczyu/3G/uY2Dk+mD9//pjHcfHiRWMf0/54/vz5Ma/HtK9aZSp5G4wzMAC4CiEPAM4jawEgMsjbYExgAHAVQh4AnEfWAkBkkLfBzOdvAQAAAAAAjDPOwADgGlypGQCcR9YCQGSQt0MxgQHAVQh5AHAeWQsAkUHeBmMCA4CrEPIA4DyyFgAig7wNxgQGAFch5AHAeWQtAEQGeRuMCQwArkLIA4DzyFoAiAzyNlhEJzA6Ojr0P//n/4zkKkets7NT6enp4z0MWw4cOBAzY42l7cpYwy9c4zxy5Miwj69atUonTpwIqdbMmTPHPJ5oRN46g7x1RqyMNVbGKTmbt6PJWsm9eZuWlqZFixaN9zBsuRL3YZOLFy9atqelpRlrpKenx8R2jbXXv7u727JPRkaGsY6dPmPV2dk55nxLTU0d9nHydihPIIJTOitXrlRjY2OkVjcmjNUZjNUZsTLWWBmnG8TStmaszmCs4Rcr45Ria6yxLpa2NWN1RqyMNVbGKTFWjCxuvAcAAAAAAABgwgQGAAAAAACIehGdwCgvL4/k6saEsTqDsTojVsYaK+N0g1ja1ozVGYw1/GJlnFJsjTXWxdK2ZqzOiJWxxso4JcaKkUX0GhgAAAAAAACjwVdIAAAAAABA1IvIBEZ9fb3y8vLk9XpVWVkZiVWOWk5OjpYtW6aCggKtXLlyvIczxPr165WRkaGlS5cOPnbq1CkVFxdrwYIFKi4u1unTp8dxhJ8Ybpw//OEPNWfOHBUUFKigoEC7du0axxH+ydGjR3XzzTdr8eLFys/P109/+lNJ0bldRxprNG7bixcv6uqrr9ZVV12l/Px8/eAHP5AUndvVTcjb8IiVrJXIW6eQtzAhb8ODvHUGeesM8jYKBBzW398fyM3NDbz//vuB3t7ewPLlywOHDh1yerWj9pnPfCbQ2dk53sMY0e9///vAgQMHAvn5+YOPffe73w089thjgUAgEHjssccCf/d3fzdewxs03Dh/8IMfBP7xH/9xHEc1vGPHjgUOHDgQCAQCge7u7sCCBQsChw4disrtOtJYo3HbDgwMBM6ePRsIBAIBn88XuPrqqwN//OMfo3K7ugV5Gz6xkrWBAHnrFPIWVsjb8CFvnUHeOoO8HX+On4HR0NAgr9er3NxcJSUlqaysTLW1tU6v1rVuuOEGTZ8+Peix2tparVu3TpK0bt067dy5cxxGFmy4cUarzMxMrVixQpI0adIkLV68WG1tbVG5XUcaazTyeDyaOHGiJKmvr099fX3yeDxRuV3dgrwNn1jJWom8dQp5CyvkbfiQt84gb51B3o4/xycw2traNHfu3MHfs7Ozo3aHlD7ZKb/0pS/p85//vKqrq8d7OLZ0dHQoMzNT0icBcPz48XEe0ciqqqq0fPlyrV+/PipPrTpy5Ihee+01XXPNNVG/XT89Vik6t63f71dBQYEyMjJUXFwcE9s1lpG3zoq1fTcaM+HTyNvwIm8ji7x1Vqztu9GYCZ9G3oYXeTu+HJ/ACAxzkxOPx+P0akft5Zdf1quvvqq6ujpt3rxZL7744ngPyTU2bNig999/XwcPHlRmZqYeeOCB8R5SkHPnzun222/XE088ocmTJ4/3cCxdPtZo3bbx8fE6ePCgWltb1dDQoLfeemu8h+Rq5C0uidZMuIS8DT/yNrLIW1wSrZlwCXkbfuTt+HJ8AiM7O1tHjx4d/L21tVVZWVlOr3bULo0tIyNDt956qxoaGsZ5RGazZs1Se3u7JKm9vV0ZGRnjPKLhzZo1S/Hx8YqLi9O3v/3tqNq2fX19uv322/XNb35Tt912m6To3a4jjTVat60kTZ06VTfddJPq6+ujdru6AXnrrFjad6M5E8hbZ5G3kUHeOiuW9t1ozgTy1lnk7fhwfAKjsLBQzc3Namlpkc/nU01NjUpLS51e7aicP39eZ8+eHfz3888/H3SV4WhVWlqq7du3S5K2b9+u1atXj/OIhnfpTS1Jzz77bNRs20AgoG9961tavHix7r///sHHo3G7jjTWaNy2nZ2d6urqkiT19PToP//zP7Vo0aKo3K5uQd46K5b23WjMBIm8dQp5G3nkrbNiad+NxkyQyFunkLdRIBJXCv3d734XWLBgQSA3NzfwP/7H/4jEKkfl/fffDyxfvjywfPnywJIlS6JyrGVlZYHZs2cHEhISAnPmzAn8/Oc/D5w4cSLwxS9+MeD1egNf/OIXAydPnhzvYQ47zv/23/5bYOnSpYFly5YFvvrVrwaOHTs23sMMBAKBwEsvvRSQFFi2bFngqquuClx11VWB3/3ud1G5XUcaazRu29dffz1QUFAQWLZsWSA/Pz/w93//94FAIBCV29VNyNvwiJWsDQTI20iPNRq3LXk7Psjb8CBvnUHeOoO8HX+eQGCYL/EBAAAAAABEEce/QgIAAAAAADBWTGAAAAAAAICoxwQGAAAAAACIekxgAAAAAACAqMcEBgAAAAAAiHpMYAAAAAAAgKjHBAYAAAAAAIh6TGAAAAAAAICo9/8BvQ5/Cvg1kxAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x1080 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<cil.utilities.display.show2D at 0x7ff0d3876910>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(pdhg_tgv.solution[0].array, u_cvx.value, decimal=3)\n",
    "\n",
    "# print objectives\n",
    "print(\"CVX objective = {}\".format(obj.value))\n",
    "print(\"CIL objective = {}\".format(pdhg_tgv.objective[-1]))\n",
    "\n",
    "# show middle line profiles\n",
    "N, M = data.shape\n",
    "plt.figure()\n",
    "plt.plot(pdhg_tgv.solution[0].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([pdhg_tgv.solution[0].array, u_cvx.value, np.abs(pdhg_tgv.solution[0].array - u_cvx.value)], num_cols = 3, origin=\"upper\")"
   ]
  }
 ],
 "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
}