TVM_Ellipses_LoDoPab.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "TVM-Ellipses-LoDoPab.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "uco9yMnAmbIm",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "86adc580-d1d9-443e-80ba-0d7ab1740dbe"
      },
      "source": [
        "import pickle\n",
        "import zipfile\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "# update path to import from Drive\n",
        "import sys\n",
        "sys.path.append('content/drive/MyDrive')\n",
        "\n",
        "import numpy as np\n",
        "import torch\n",
        "import matplotlib.pyplot as plt\n",
        "from torch import nn\n",
        "from torch.utils.data import Dataset, TensorDataset, DataLoader\n",
        "import torch.optim as optim\n",
        "import time\n",
        "from time import sleep\n",
        "from tqdm import tqdm\n",
        "import os\n",
        "\n",
        "seed = 42\n",
        "torch.manual_seed(seed)\n",
        "\n",
        "device = 'cuda:0' # \"cpu\" \n",
        "\n",
        "use_ellipses = False"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qUdNaX-9LRQq",
        "outputId": "56cf1f4a-11a8-4d86-da16-720e190098af"
      },
      "source": [
        "folder_path = './'\n",
        "if use_ellipses:\n",
        "    data_path = folder_path + 'FFPN-Ellipse-TrainingData-0.015IndividualNoise.pkl'\n",
        "else:\n",
        "    data_path = folder_path + 'FFPN-Lodopab-TrainingData-0.015IndividualNoise.pkl'\n",
        "\n",
        "if os.path.isfile(data_path):\n",
        "    print(\"FFPN data .pkl file already exists.\")\n",
        "else:\n",
        "    print(\"Extracting data from .pkl file.\")\n",
        "    if use_ellipses:\n",
        "        with zipfile.ZipFile('/content/drive/MyDrive/FixedPointNetworks/FFPN-Ellipse-TrainingData-IndividualNoise.zip', 'r') as zip_ref:\n",
        "            zip_ref.extractall('./')\n",
        "    else:\n",
        "        with zipfile.ZipFile('/content/drive/MyDrive/FixedPointNetworks/FFPN-Lodopab-TrainingData-IndividualNoise.zip', 'r') as zip_ref:\n",
        "            zip_ref.extractall('./')        \n",
        "    print(\"Extraction complete.\")\n",
        "\n",
        "    sys.path.append('content/drive/MyDrive/FixedPointNetworks')\n",
        "    sys.path.insert(0,'/content/drive/MyDrive/FixedPointNetworks')  \n",
        "\n",
        "state = torch.load(data_path)\n",
        "A = state['A'].to(device)\n",
        "u_train = state['u_true_train']\n",
        "u_test = state['u_true_test']\n",
        "data_obs_train = state['data_obs_train']\n",
        "data_obs_test = state['data_obs_test']"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Extracting data from .pkl file.\n",
            "Extraction complete.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7KSnNLuuUyNk",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 273
        },
        "outputId": "f652ef69-ad9a-432b-9fba-27a806c4dadc"
      },
      "source": [
        "plt.figure()\n",
        "plt.subplot(1,2,1)\n",
        "plt.imshow(u_train[0,0,:,:])\n",
        "plt.subplot(1,2,2)\n",
        "plt.imshow(data_obs_train[0,0,:,:])\n",
        "plt.colorbar()"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.colorbar.Colorbar at 0x7f1f116e7790>"
            ]
          },
          "metadata": {},
          "execution_count": 3
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAADvCAYAAAD4ic/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9a6wtW3bf9RtzVtWqtdZ+nHPuedxHu91P23FMbIJlJwQhB5NgIIq/RFEIRBGK1HyAKOEhkiChIARSkBDgDxHQIkAQISYxsWIhEwdFiiAScWzHzsPv7k63+97uvvee136sV1XNOfgw5qyqtc4+957bfU/v43vrf7S1z16rVq2q9Rhz1H/8x3+IqjJhwoQJEz54cNd9ABMmTJgw4flgCvATJkyY8AHFFOAnTJgw4QOKKcBPmDBhwgcUU4CfMGHChA8opgA/YcKECR9QTAF+woQJE34TQUT+XRH5RRH5xyLyl0Wkfuq2kw5+woQJE54f/qXfvdQHD8Mzbftz/3D3U6r6Q0+7X0ReA/4O8J2quhGRvwL8pKr+z1dtX3wdxzthwoQJE54RDx4G/t5PffSZtvWv/PrtZ9isAOYi0gIL4CtP23CiaCZMmDDhOUKB+Iz/3nVfqm8A/yXwG8BXgTNV/ZtP234K8BMmTJjwHKEorYZn+gFui8jPjn4+M96XiNwEfhj4OPAqsBSRf+Npzz1RNBMmTJjwnPEs2XnCfVX93ne4/18E/omqvg0gIn8N+GeB//WqjacAP2HChAnPEYoS3j8xy28Av0NEFsAG+EHgZ5+28RTgJ0yYMOE5I/L+BHhV/WkR+THg7wMd8PPAZ5+2/RTgJ0yYMOE5QoHwPgV4AFX9s8CffZZtpwA/YcKECc8Z71cG/14xBfgJEyZMeI5QoL2mhtIpwE+YMGHCc4Si7ytF814wBfgJEyZMeJ5QCNfkCDMF+AkTJkx4jrBO1uvBFOAnTJgw4blCCMi1PPMU4CdMmDDhOcKKrFOAnzBhwoQPHEwHPwX4CRMmTPhAIk4Z/IQJEyZ88DBl8BMmTJjwAYUihGtyZp8C/IQJEyY8Z0wUzYQJEyZ8AKEIjfpree4pwE+YMGHCc4Q1Ok0UzYQJEyZ8IDEVWSdMmDDhAwhVIeiUwU+YMGHCBxJxyuAnTJgw4YMHK7JeT6i9nuuGCRMmTPiQIBdZn+Xn3SAi3y4ivzD6OReRP/m07acMfsKECROeM8L7pINX1V8FvgdARDzwBvDjT9t+CvATJkyY8BzxHDtZfxD4vKp+6WkbTAF+woQJE54z4vNR0fwh4C+/0wbPjYMXkR8SkV8Vkc+JyJ9+Xs8zYcKECS8yzGzMPdMPcFtEfnb085mr9ikiFfD7gb/6Ts/9XDL4xA39eeD3AK8DPyMiP6Gqv/Q8nm/ChAkTXlQoQvvsVgX3VfV7n2G7fxn4+6r65jtt9Lwomu8DPqeqXwAQkR8FfhiYAvyECRM+VFDleTQ6/Wu8Cz0Dzy/AvwZ8efT368D3P23jSmZas9y7TZwjLmvaY6Ff/J51MvlhwVqfct/h/gQQ7bcRr3gfcShddKgKs6IjqNCtS4oV+E0LIQ471NG+nnoMYj/99qPfAOLs8XrwQBntVJ/2HHr1efWPEXAOYrz6eK96jifOSYZf4oiVJ8yEWGHvlRye8AgOfBEpfaCQQFBHFx1t5yErDTQ/fPT3IfKxiSJBKNbgL3borERChBD692UTL2l0ez2dJhMmIO9ro5OILDF25N96t22vrciauKXPANQs+H7/e4f7nOCOlpz/4Hfw5vc54kwhWuy9EgcBXCUFGqcQBYnsB6kUtCTFZRXAgXpFS4Ui4ueBm6crvItsmpK283zspYfcXy+5/Lt3uPPzHcvPP0IenaNdZ0EZIIb9Y3P+ydtEkLq2525b6DoLSHbyyKyy+1QhpgPXiJQlGmPazI0Cdd4kQtehTYuIgPegcTg2sH3XM3S1sefV0bx3cVAU/fPTdWgY3e8kHb5AUVhqIoLWFeH2MauPzHn8Sc/mlYiWijp9chHyiqs7ZvOWk8WWygdWTcnjx0ti64f3Mwp0gkTpx9JLPFh8AK0U2QnL1x0f+b/eRmcl6gX/5mPio8cA/N31/8mECdcF5f3N4FV1Bbz0LNs+rwD/BvAto78/km7roaqfBT4LcCK3ngzdsxnNsUOLZ03b9yE5OZUh1ZWD5LaXpgqoU3sup+BgVjd4F9m1BQqcLLas24qzn7/Naz/dUH/hAVysUjY+evNELPCl/4t3aA6i/e2uD45XQUNEvENEUIcF+fyYqH2g3cvG8+948HqJG4K4OFQ1JdjRfvL2Tobb8r7BFgEnaaE6eO4QwDnkck2xazm+PMJ1JzzSgs0dJc7j8P6NVue482yj0DYFZdURgiM23t4cR3oP0uLR2NWMqKCS3sP0O1Wv0FLZ3la2r50we2tFnM/wZYF4j4aDxXXChGvAB23gx88AnxaRj2OB/Q8Bf/iZHqkRpETqGc2JWBaYs7XxF3tMOxygD+Sj7UQP7mO4Tx2QM36nuDIwKzvWTYmqUJcdlQ+8/g9f5hM/uab80tvobpcyXpeCeNprdEAKkt4PC4DG/YAegt1/1fnHAN7ZAiGCpiuAvWw6w+Urh/37VBXJ2Xt+/hzQx8fhrlhk0rmoKtp1iHcW4MfBXTUtOOm82hb34DFH2x1+e5v731WyveMIy7i/SCfqRVtH1zm6TfoIdsMbJY6eKtNgV2E9fdZ/CNKiHQUtlO4ocv7xijv3N2jh0LKw13cK8BOuGYp8sAZ+qGonIv8O8FNY6PwfVfUXn+nBkgLmck47ouX7BPAq7pkrbtOn/H+0L810TaFoEcErUihFFVhvK5xTThZbztc1lz9zm0//+EPkzYfoiDZBHXiPhDAEeXFDMFS1gJ2DarovZ9KWpaeFIAwZtKoieYFI1AsAzhtF4txon+NMPi0KISTKXaEsoJjZtjnA58Dfv0bRrhxItIx3SFEYFQO23/z/EIxaGmX/mTKSbcPiH3yZVx/e4f73HHH5rY7mVrDMPL8PkbSqMgTtQzonB/RSobGFXg6/JJL2FQQK5fzjcONzNSqCViVSlbDbMWHCdUKB9pq8aJ7bs6rqTwI/+XU9uCyJyxlxlnd2uHNGRcGrnvyK+55S/DTePUKhiI+IVzQK5SxQ+sBbn3+JW//Acff/fQsePLJgXng0CkS1QCxiQTFlzFJ4qEp7gqY9eMIInVEtCkiRqIQcwHF95qkhjPZXod0hv2/HYGX6YP9PQVe8twXHCeOMXbxL24bhPnEQA5Iz/ZT1Kp1dSagaH58oGslXFuM3xjkofF/YLL7ykHsXW45fv8Gjbyu5+ERA6ziqgaTMPGfnXpFiVBTJa6WPqE8PimlRVqOuJJctImiA9lbg8admHH2lIy5KinmNXq6u+IBMmPDNhEx+8GOId3R1QZgNASQnb4eF1v6K/WkYFVPz330h1pGKqim4O7VYKErTeNrPnfDKP4ic/uIDeHg2ZLyqls06ZwFw14wKqRHtQFRhVtk2LhU7x0VTsKBJKlqOccCjaxdGVwQRVXnygiUH86qEYEF5b5tUyFVxdmUAdjVQlcZ1b7fpxXJAsCCvI6opZsJb0HyFkheY8XEnXp6ysGz+1x9Qv7Xg6CtH3P/ugualAF6HjD1Y4CbK8GYevqFebTsBnKlmNHH0WXEjwa6KLj4G5crjtxV+NUceXs+otAkTMpTn1sn6rnghAzxArBzv1hvwVFprpJR56kYpWOAVGdEDGqELBcVXZtz5hcjxr50hD8+sUJoCnZKCZ6ZJqtLokZRxQ6JYQuKnnQzcvJch0I8za0mBMx/E3olGIHHggScLqVlZk+SXStomLypdN3pMsO2LAqlntv1ma9tkSihz9bmomh8brQCNDIuSOoyeEbHFjsIoIYAQkG3ANy2n64Zic4O3v6dgey+gsySL8unNCIKKgMfej/Fb5dSutDRx7uOrNx3qKtIK7Ulkc9vjdwXlWYnLV1ITJlwjpgx+DHHE0qH+Kan5O7xWe/FcDrL3g32otwCfI4RGh0ZB1p7TX4XjXz/HPRjJIHPQxCgBka6nOcS5PQpjjx8fnddw/E9KJ5/IiHOmHN0QUMf8PlxZJO3VN3kxGFFJppZRpKrQWWULwAbj2nOmLg7coLy58nmKArrOgntRGD0z3i4remKELiK7hmWItIvbPI6e3R0hLkKvWjJZZFoEM//en1B6j0SG976/pBuxb0GI88juJpQXjvqoYjafP/H6TJjwzYSqTBk80HPYAOpTHe4pMf7ditKHgf5wP5p5X58CduaDW0f9pufmr6xw98/Qph1RFOwpUjQK0rSWwXtvtEzWn4NltIfZdk+1SL9oaOprEhFUBz18rzv3Q0F1/P+xbLGnedIx5EJr/xLk4J5f31lpmXYXhiAdgvHq+Yojc/Qj2Wd/DkVSqFQlmrJkaTvj4LPaZ6xg8R7Z7jj9xce4cMrj4Fm/qmgO8mCvf5BUm4A91UyvrEnnO35Z+6APRGhPIs2Joz0uqI4XTJhwnbAi6/VQhS9WgM9wQixl/xL8sGB6yMlfVXhN/OyVPH2iZyAF92AZpFs7bvx6pPyNkRQyHRNdlwIhQ3COKUCWQ+PPXjY+li9mJUpSvGhucIqOvQU+BXcZXwVklcxYQePS//P+MlWTm5DycXOQXZel1Qe6MCxQXRgKvflKJTqjlJJE8onj8R49WqCSFrqQ6wwj3X/etrTzdustJ/+oo1jd5GFbcvEphcqK3ATS+5BURuVo8QVrRpNUJ8mF2NF7LQrSOsIi0h47dqeO2el8/ypqwoRvOqaZrD0kNdWEyqFOccG+yELK6K/6rh4E/H4xGKkx9jZPXasIaOcGxcbOsXzdcePn37KiY0zFQFcYBTOWOo4RI1p4pAv7wd3lAqIOwb3waFkgbWeHfBh8ZjPTnmdKJunN+6uD8WN6TtyuGiRr4g+1394kkX327605iaY1uSNYcC4KUwjRQpclm+lxo2On8HYFAEgbkO0O2TborhmeYxzc04/6QbO/+OJjyssj/G7O4+/qoIyp8KpDJu8E8ZpqqaNzHa/2+b1PN7sOwgzaE2Vzx1E/qq7uN5gw4ZsEK7JOHHwP8Y5Qsc+zZpolB+yr6JtxRjfeX5/lWypvjU26FxxoHbOHjtv/uDHFjPNACthRUZIEsp6h6206hhEtAgMl0zckWYCTUWCjY/BKgUESmTnyOOjge0S1TLo/Tx0WkvFtOYsf0UlSuCHAhYDG5IqRg3uWd+YFqZjZYpqkk3TdsMh0tijFeUV3VFG+eY5sdklxk66GQjSrgyeUQel587GHSPnmOfd+uiP6Y85/SzfUQ3JzW+P6Ymwf15/yPcmMFwGkE0IdaU4925seLa4ne5owIeOD1sn69UMsIIVSUqHvad/o8WOe8vuqh/Wt8PvbFeee5VeU+kuPk0bcjzLlYJl4UaD1DNk1Q0BOwVOiWmaeMlmwhUpPjpBMx4wz8TGdkxql9qiYcRNRVunA0LkKFpzHBd2MLNksCjvOTJHEiGx3sGtMepmDu/cwK9JC4qBjXyLp7Bw1H38ELZxdsaTgryEiVUm8fWpF1fV2WPiyZ026ctFRHcGdrbnzD0va45rNqx1U6XLNp0uxVuw2GIK/U6NyRp8F0eEKTzpBS6WbK7tTNwX4CdeKD1wn6zcEJ1AWhBn7/Cr7V+ZXat+fksHbfUOxTnOmmK8KWsfsoXD0RoucXexfGIy061IUxHmFNLUFN+/QrB7JwRmG4CjeAuthANaRquWJAqZ1kGY+XA4D+LiY+k7Ix7xnSzA0UQ22A6bTl7LslTAaE5deFAPtM95H4ZCgQ01ABFnU6NGC3Z0Fs7dWthA4h7rhvHQ2dMXmBi3pAtXrj7nxubuE2tO8hHHyeSEOYpRNpmry4pxZqAOxzfjKLpZqdhf+er5cEyZkPMtA7eeBFyvASwpus4owz4FsfD9DUD4osB12rcr4MSPGRK/I3v2lY/E1Zf7GpdEWVxXlknJEC2eqkSJl+GWBekG6aFltU0DrB6+ZlOWOC5Q6LkIeInWO9oHVuVHB9IoPSaZqukELr4le2nOcbNq9KwQBtCiGxUgVXc6h7ewKRVLNoBx9RNSuUrRw+G1n5wZWbJ1VhOMZsXTItt0/5sKjmYcvi3QVEHpaSZqWk1+9YHvzhFh6uptqRVc/BHijYLR/D/tWr/x5SO+7uhToI8mjRrmm+taECYB9zNv4IQ/wkoqQMq+JRzO6OUNgzshf5EM+ntHth+gfo70tcN9FmQp68zcdx1/emSwym2plPXhuTvLeVCCbFkkUDM7Z7n1JXFSEeUkBuC6g67XpxM8uLBvPRdAcYHNm3jcIJSlmm5Qsh0VV7+ktCTKymuZgQZLCgyvJjo89L589Y/J23qFdhN2OuN0SX3uJWM0pywI5Xw30Sz2z16RN1gVB8WeXqbHLnlvaDn+2YXG2GeiZYN44WpSo91aQHats8nF7hz9bcfNXa8JsxsXcmXzSj4J8//7pflAfv+/jz0cUcEq30HdtmJsw4XnCKJoPcYAXEdzpif1/NkNLTzcfy2EOkDP28ReafTnkXtd73jbTM7mIFwVZFZx+ITD7yrn5qBcWyAeKpBgCrAiy3lrnZ9Te0lcuQE+O6I4qupOaIoJr2z33x7FBWKZgpCwHeWXbjTj/A+fJHOwLb9LF7OSYi7IpaPcySe/RRY00rdUDMicf00uWrROSd7wmH5uwLNm+VOKaQHG+MlloB3psrm/SWtbudi16sbL9Zhqnae1+GIquZYkua9qbc6SLlF9b2zZR7conXwU5o7pmrz/mZn2LWJZcfirRNGWExvWeNeJTV+soUx9e5NF/01uspRKnAD/hmvHh7mTN3O9qDSKEuiDOSIM6RlE744rgDjxpJcz+dlkaCfSX/ss3HEdfOjflTAx2jR9tuMZ44MVecTTpxLVLjVlOkO2O6o3HSRIo6OmxZfq7BvBDt2fKqiXGgToBekOvTOdkJUqIgw1Cpn3GWXuIvQxRMp0TI3K5RptmMCHzfnCGjKG/XaoSKZdI4cl2Z7h0VeEqtJ6hM4/srKkpLyomnfS236weyueUFjJd1Gjp8esW//DSpJlgOvx0zlp4249z0HbMv7LiRn3M7qWC9k5rfWXpQorsQZPllCqpC9Y+Kzno9x+ZxN1PGfyE64TyoZdJCrrdEZsW7z1hlmwKlCeCuwpXmoeNdjUk/qPtdC97t+AunXD8pYh/eGlZbFLOaNvua9SzwVfOnFPzkHWepuPbNYP3TOFNKdIN8kH8oAOXNlEfB4XUfoJT4UYFWx05TaZTHQfVjNz0lLLnfl/irM/Jj1wum3af3nG2kBSPd1Slw22zNbEdg7vYWuDuPWbi/oKlmmgsN9AyVZlURR2u7ZCd1QC0C9Y9nIqwvcVB7sDdNMzf3LF8fcHjmzIYk+VO4950bFSfHloZYGiGHj4DEwc/4VpxfRTNi/HRF/pvqziHFvLeCmMjPn7sOjleG0wemZ8LiFBcOpZvbGEz6Nrz2Lv9/adgORqtt6dOCakLNP+0XaIzumEByPRMl8y+yvIJhYvUMwvCkDTz4cmmpX7Ax+FVzbBI9ME9Oz6OLAoQsSBflWhdoXXVLwxutcXvwv7zhIDsGjufRNFILq5mPn1MGXnXN0NJjNbQ1Xa2OMyqgQJLmbsWbuDr0/785Y7FmwHZmTfQnld8anzrM3l3xWjAHPg1fyCuKGZPmPBNRESe6edZICI3ROTHRORXROSXReR3Pm3brzuDF5FvAf4X4B72Dfqsqv6IiNwC/nfgY8AXgT+oqo/edX9FgbgGvCcWMiqk7lM072gNDMbTHhZghaFglwqrEoXZQ6F8sLKgns3EsvNir+HOmvRgk5XGk43GOJQ+dqMBIN6bXDDdbrc5iCP9e9aJ52NoRkqU0WvEiIYBhkDs3HBpc7goiI0OHDxy0hVGfqxYg5R0wYZljGoBEpO3zNgDJwf0fC59wVhHUsyYHm/noIWZC0lIVwKFN/li3pdqf1UiXaR+0FGeVTSLMEzbGvvH5yxebBCIKuYPf/DWvOvnZcKE5wxT0byvPOGPAH9DVf+AiFTAUw2XvpEMvgP+fVX9TuB3AP+2iHwn8KeBv6Wqnwb+Vvr73TGbIVVlGWA20HovX85R1+uewiZRNlfRM4uvKXKxHnmvJA/0NIAjNycBA5dNCrTeD8qYXIAdWwUEkymKcz13LV2ahBTCqNuVPsBJm+efun7B2BvTVxaDrj6EgfMOw6LRH0O6kiAEO955jS7nJnOclYMiJ+8HLPMWcE0Hmy263fU0k5YFuqzR2h6rTbung7fzHrJwybRM0/ZXNJIeo0mlJG1A2sE/Pi4q4mKGFo7yomHxpgwNTVn9lLl1oZe8qleTVXLw3ve4Hv5zwgTIKppn+3k3iMgp8M8DfwFAVRtVffy07b/uDF5Vvwp8Nf3/QkR+GXgN+GHgB9JmfxH428CfesedxWiX7ou5aczHqpjD7D1z7k8J/pmjf8JtsBgu4yUKbissv9L0Q6VF7KXQ8XCOkb97j+ykmGWLTWut+alIOKZ38qAOUUU2u2HMXeH2gqPmYmtW8WTevZdK+v2GIzc8TlTRkzTbsO0QLdMC5AafdxiolRjt+FNxV2O0+oP3yK5h9qWH0HZ2TPnqYlbZD1hA3jX0HvWHWX0UZDtupEo1icIj69YWhqZN3cE6+N8s57jRVYmWyskXO84/5YneAjk+ZfHY+y9RbBqXV1s0gsO10h/OUI+Z0vgJ14tnpV+eAR8H3gb+JxH5buDngD+hqleOLntfiqwi8jHgnwZ+GriXgj/A1zAK56rHfAb4DEAtS1N8jLtB4Z0Tr8Mgn3XSV0ALhSLSm48FKFfC7P7mScsBAJxlziP/d1OipN8jGoYQ0KOFZc+JVuk5cI22cORMWgSqYrDZ9c68XLrQLyq6HfH/zo5LqsoC81hK6Zwda1GgVYF0g3qFtoEuXVGUpZ1j29q2J0vLnp0VgSXG4WXcNf3rL1njn56LXYNbhREt46wYna9wCuPYtfB2TiJWeM2drl1AVyu7KhjvN0+1apoU7AuorD4xf2NF9fCUbe2QOo0RlGj2BVHoJ0EVlsFrpWiQpL4af4SmDH7C9UF5Tyqa2yLys6O/P6uqnx39XQC/HfjjqvrTIvIjGEvyH1+1s284wIvIEfB/AH9SVc/HhUNVVZGriZZ00J8FOHEvqWnQk19KUPZ9CQ4fzMDDHh7PwbOpWIAXp2g3XBqUlyCXm3wsw/CKcSv/2BfdYy399WzoDk3t/BR+CL6Seiy9gCT5Z+L1BzfHtJh0Rk/I0ixt9XINbQqAY5vfcTZ/2M2qirvY9P8Hem/3PljnLHo5Z/2xE+ZvrGzBydSQc5axj9U5zpmzJwzeMWE4Z+lfp9Spmrp7ZZdqB1kTP6ZwxA2LXzre/i0M6TmKolcbuYs19f1TdncTx+7UdPDq+mJrLpjnGosWCp3sL/5TAj/hmvEeVDT3VfV73+H+14HXVfWn098/xjvQ4N+QikZESiy4/yVV/Wvp5jdF5JV0/yvAW++6oxEfLCEiIdUL87dfD36PHyo8NXO3g2R/IRArxpXnimybXsnSNyKVxaBwyV2mIhbcs8IlhH0JZGrv74O+iGm965ll4JJpg1zETZTMuCPWu6EQmjJ305zrfnAfd8Dmvw+7Q7OuvD9nQZNixW+jcd/Z9TLJKnsNfbZfmJV9927/POP3y6faQlkYd59eB9pu/3jy9nlhOLSB6KkqNzwfmAInKvXjiDTpTRbMPriIwye3V8pgippiZE0wBfYJLwBUhU7dM/28+770a8CXReTb000/CPzS07b/RlQ0ghH9v6yq/9Xorp8A/ijw59Lvv/4s+9PE90qfwT9tQ3mCjnlaMbbfJjsTim3sWmF2Hvs5pBKsI5OyQOvKWu2BfnpTP2Q6BcEuDgHMe6NmRvSMpExYC99TGIIVbTXvl8TPJ+pFclD3Sc85HjQy6qTNXjW9k+SBoke9s+A9DvwidlsbqN5eDQtLDrzpPHprgqSwUYmDx71Pi18u6PYKnLRAji2Q8+Iz0vID5opZVdDsu3FKVVkdoyjSa2LFZS081UXEN54YBLJvm1e0iNBlDX3+IGB0TZsv8XhK0XXChG8u3udGpz8O/KWkoPkC8G8+bcNvhKL5XcAfAf6RiPxCuu0/wgL7XxGRPwZ8CfiDz7S3nKmqWjBS6SVwOB3icySNvE5BY1R4lZzsHhZYRW28nrPfroXq/EDiV3h0bjYJ8jjZDPRGXNGomjSs2hp6ZAiQWWmSg2WSFkpoLJseB0ZVCObpIlLY1KjUHCR1PWT8I659fCWhVZmKpfm50v5g5FyZMuwu+9BEJARTsjiTY6oTpCz25Y/JQlhCNHMwEbstvS7q6kE2GYairbSjFzyfaw7yIdiCmJuaZjMkL6pO7O96NkhE+45dW1CKVcBvCtqQr7SSq2ShezN7RRgNBaH/bFxTA+GECT3eIwf/7vtT/QXgnWicHt+Iiubv8PTq1Q++9x0OxUwJmqxftS+mCvZb3ehLfWhGlrbZ+zNTNGKX9tp5/E6ozhrzgwGjU9QoG0iKmG5H75UuAhrRxgJafPUObtciF2t0k/hvP+psHf9kHh3ozcsc5keTKRsAZ9YE4sxDXusS2bbW3p88YSjLRAWNMu/s+1IUg1LmUKcfApoLq2U1ZOUwBHHVfnTfmGLRuiIczfCXg2QSVePa25bebz5PfKqr/gqjb4hyzigfEathaJWsENLVyK4Ztk3HJNsGmpbqvqM6n7G75WAB4rLfjtgiNOpsk0zLF4oEsY/UFOAnvAD4kFsVJKQuUt9G5KCZ9DBwH05pylax+xvRyyNdGXFOiVFwDbhNO+y6rgZKI0T09MjiwphyUEV3Dbre0N6scaGiBFivh2DprXvTdqr7QTQP5U7k8d4VQmpuUidwcpS8bdqB18+SyK7b58PzMY+PFQbPGTAFTqKX1LveMGzfZsANzU+5DlB4Yl0QFpVZAzetBd1cWC0L9OgE9T45aqYCcxuGhSb3DiSNv1almaAlfb12u8GLJ1FRUpYj50zFna0oz2/iGiGSLwyU0KaXMow+G2KFWDMjUyRLKqcgP+EaMQ38GEMVaSMuB/iBjUn3y0hZjAIAACAASURBVD4Hc6CT7zvThd5/RpzivfnJkCgaaQMym5kkrzf0SgE3DYju4VzycinMQVIV6RIPnqYdjRua+kPNU5TScBAo+uxYqtIy63TOiA21VrAaQNasp+cX1eTjokNH6ngRyRjz/1EHBYzIcI4xUU6O/r6sXNHa6Jg4L+gWBX4X8ffP9/3kc5DvkkfeKFPvz5e073TRYldmsVfRaNNC2/SToKSuh4wehgWoCxQ7xXUyDLly6b3MVqG50p6HuIw9/yf+fcILgPdRB/+e8MIFeA0R6WLfsNJj/GV9msMkw82ainEU0cyt8mZBbPHIAXTsqzJWf4wte5N5FlWJNC3lg7UtELsmBW6G4JQDVOo27RcPGMy4eoVMNGno+DTazkbqZeojB/OiAG2HYzt0lczyxlyArcrkwjhcUYwpk/640/6kS6+Bc8Sjiu6oJHph9vbaulpVh2NNSh/pwqCayeecZZDZUz+9fkBy0bRtxDuUCpkJMq+Ng8+F2pG9gYjsqaqc077MMIYqw1fIGZUnmbOZMOEaYbqMb0iw+HXjhQnwmi/l1Yp8PUWjB7+ftXKWW9h91k5bdkzEPEuypDE1APUBr8kZbhyCKwyc8XqLe3D+pN9LjHuyR2CQBI6z7LHyZbz/zKXvrDZgyh2GoFmkwuTYnyYfc7YMzjbAo8Vmj28fc+Ipw+8DakhGavMZ3bIkzBzFOuDO11efS74i2I0a1MbnGXLFOxVBso4e7WsGMi/QqiSmWa007egYbV/mGU8/0akoxlz96PUfmYpJKspflQdMmHAdmCgaSIHGMnjf8PTs6ynZ+95VeyrQilOci8Tg97eXwepXkoIGsCwyYxxMQ+ytCHr+PMsnR94yklr6CeHJLByGrHsceNNiozn7DtHsDHIR1qdW/1AMi086Ni0LpEsF0LYdLIgP6Y58Plknf1Ua7BzNvWM2d0okKPXb6YphVu01J/U/o6KyqA469tzXkPbZm5VlFDbiLy5r4rzArxrkbGsF6xz8U0avhZnPqbP3svTBsqFcPB9/DEQH6iZRNToR8BOuGRMHn6FZfhfxO/a/wPn/esVt413k2xxGzzhNxVWTSKpXQkVSjHTGoXdhaK8fyxOz3C8mG4Lkgy6l7FE6moNZVaK3Tm0263qHrDZ9O74drwxdsMUoGObTiSmwJz8aUtu+lgN3T4yD9FHEzsOblYC2Zuxl9EiiO4rUiOQcortery/dk1cD3Ss3uf/batRD/UDpFiX66k38xdaKp8FMxNju0NUOOT6CkyN7bdKCpTMb5iFJcrk3rATQ5Zz25VPaE9tudn+Le3hhjykKO97agr80HYRIuY64zhGBqgiExtl8VqfoVVd0mZ7LUtYJE64ZOgV4o2k06bVd1lbnRDNl5rah7MsldfidZZHj4drOKcfLLW3wrBpPN4ewrPBtBy4VOnfQzyzNmeiYpskBMRdgM7VReBuRt93abQ8e22FGRZPMsp/FCoxH7PUFzhTs+5F3eYYq2IIHxJNF78ku0Gf3str0Vgg4DzM/NGWFsL+QNK2pcQ66VuNyzuqTJ7z12wt2dzuqBx6/gVCbqyPRzl9CtOB7+4RYONQL5Zvn5i8Tdbh6ydbFhTevmA6TpIrYfFdg9nCHf7iCtx/2r5OdrxmoSdR+pF/0OYOPeHeVNpZ+4e+z+KSmkTBl8BOuH1ORNSMEpGktg8+4gky9KnHr/xZScE80hgqzsqPwkV1dEhYl3VGFLwv6+aTiBj13liMeShKBeLo0CeNYQilpzuq4mxOMwkmcuM5nQ6CtSuK8xK2bwbMl7ysH9z7zlCRvHO13vFik1wxx+57veX9juqQq7Ryzmkas6am9s+DRtxXsXmnxZwVHX4aTL3XM7m9wq10ve9SyYPPRUx5+R0Us4eWf2Zh0UsSuPADZtejxgjiv7HVqWnBWoNX5DApHcdmYD9Djc2v0AmTmenmnNm1qXbCrl1iCFuCdDraqOry/vYVweq9FdC8BmAqtE64TpjuYAjxAz0P73dVSyP7vsWri4K6YZnaKywF+dL/YjM5YOfOXSXJFstFVXSGbgfsGBnsA59EkWTRt+35xUaSzCURJRSJlkk161882tYzdNOlu3YyUJ0/JTLG6gI69bmBYgMyFi35GrGrSkjtwqbko1xXGRdB0TO3dY84+PmNzL+IuPTd+DU4/31DdXyFtQL235/eOeDJnfbdkc0+ZPRKKx3mUX7lf8M0vem/I5tFyZrRLvq0q4cYJbrkgD+EGkp9Nu1driIVRa1andcQolqXnL40C8YrL4H7i+oQJ1wkhfNhVNBmSsthiq4P393veyeiHIauzxNjcHlUYfF7ygOtDp8bMledGp8JbB+vIOVJdsqdNvLmABRafujqrNFyjcP1+TGIZBm15piScIHnYaEiBPCrEES2TKZ68yGiSJHY6KHHywnKokR8XR6uS9u4x5x+rufwWC6DHX3Dc/JUN5VsX1phUeCgHOWacFaiD8lI4ej3izlZ90VZ9mtaU5rtKbnZKFJTWZWpMMimmVs723YbhXLowXBllC4hZRaggloqIEqIQo+vFOKj0E536/ofJn2DCC4aJgwcz6UrccbGJSPRPl7plbv5pr9uIu89Jr6plfhLBBU2qnWCBONnfmj97tx/wR3y1rDY9967JGGtvslF+TPKgiYsK6ZKD42pr+x8rTLImPf1fvex3w+agno5vr8A67mQdz2jNgT0H82xAlvXldUl3+4izT9ZcfKvQ3IjUbznu/Pya6iuPbTFLPjWytWxaouKawPJrLfMHjuXnH1m2Xg72vvm5Jfnf97NXnUN2ne2rLMxjqDNHTZtilRaszXZvwHg+/lALWg0SyRhT9j5OAFKwf+ql3YQJ1wTlQ07RiHO4ed0XNnXXUKw6pCsGQ6nxlzdna4f7ybcdXg2psN5VhOgsZkSG4dIkaWSwpifdNUavjIO1c8RFTTitKc/T4JQsUyy8ZePjhiYxmka2Da4LxlN3Bw1BY6tfSGPs0gzTQ+WH98MA77w4JMVJtnewq5ECmc2MZmq7XvqJpDmrdQVdINxa8uC7Flx8DOIsMn/T8cr/t6H68oOhgSq/NrvGtO6zCvd4xfzxyoKyau8v078e40UnF4/BePjtzsYGquLOduhmg262xBDsvS8KdL2xIrv3yLzuqa5uDlQR39dU0ksUxArA48/D4fdoCvYTrhtXfKW/WXghAnwuBGrTkicb+csG19WM5etPIMeWgyxfRUfSGvu1Xs+MpmkdrhPcLkAMaJu6SfNQbDcKWK1RKHpc0d5eoE4otlvbXq313pQio8Cdg7Aq2jTDMI18yIdXBvl3ukLQEJPsMF1FJF8YXc4t2G53+3NaZ5WteV3XT0MyCqTsjysuKrS0LL65MePhb5mxfkXxW+Hk83Drl1aUv/G28eCJ95ZMBYkY5RNt7ODeFcVsKPZmionK5I/j+ayIpONv4eEFMc9mzYtb0yLRCrFSleQB6LpdE268Snui+Noy+y54YnS95wxRIJDMxaSvu/SdTlOBdcILgA+1ikZV7UufIDHi1jv8Vgizd/iGHlA0exr4sfglSW5UBTqH24FLo/Gkt/89UMCk4xJAumDKj63pstWrDf/wiabJlMRome49V4hD0TUH7Gydmxub0nMosP34LTYv3eb4ixuKt88Hvt9ZRixFgUiiarY7G8UHpgaqZymwWzcwqoTbJ2zuzWmOHephe8vRHMPia8LJb3Qsv3iJe/vx0AE7LpZmU7Kx7YCdnEk0U9epusHLRrbNfrqSufTtlrhregdKoH9dNA/+ds4CfLZEaBqaWxXtkVJXHc5FQvBoEBDMHTSaa6REW+V7meSECS8I9DdzkVVEPPCzwBuq+vtE5OPAjwIvYQNh/4iqNu+6o+ysmMa6ue0OvwVuMNAy/bY5eutgB5yS9t5RMlsN58v6/NjeiyYOGXIwmkNzsBQHXgb9eoy4840F4TyVKNsDDJVcy77HfHT2kYf94A4jI7Hh8dKanv384w5kzokqxdeME5ddOzymnlm2rQoXRhnJck68cYSWHne2RuuK5t4RF98yY3tTCHV6bQRmZ3D6+Zb5Gxe4R5eDimh4UxPNJUM/QP57rMjps30hu1X2Reem7RdO3W5tAHjq7r3y7c8Zf5Ne0bJCFgvWdwq07kb8uzOlUn7f+4am0edib8dTsJ9w/fjNTNH8CeCXgZP0938B/Neq+qMi8t8Bfwz4b9/THqPCrqHYfB1HI5juWkwSaeqT0Zc8gusY1BpgdEHuvNSYMshRoO5CXwzs56rmCVTwROYuI2dIgCf82fNQjNzC74eCbbHqKDYzqyGW3qwI2q7PcPWkprsxpzmt0EKYvz5H2kB7e8HuVkkshPKipls4Ll/17G5BLHRgrILgWqU6a3CXSeaY58pC31xEPi/nBiOyMc/uRpl+GiJuReUZbrVF1huQJEFt2ie9e656XaKiBAvy3sOtG2zuCm7R4VNmHoLbq8f09sCHdZnJhGbCC4TflCoaEfkI8K8C/znw76Uxfv8C8IfTJn8R+E94lgAv+5mvti3FRq1AIVfUyq76Aqdt+8z9itdUoljAyi6IOfBko650LD3PnTsnQ7IXhiGgjwunsB+wRrNY+8eMm552jXnIYEFVk/tjcbbh7s9Ir1rRLOEsPO2tBZt7M9Z3HLtbQrdQju+eUG6U9V1HcwPCTJFY0C2VuGiRnaN87PAbO5ZQw/aWsLlX45qAe7wynXvWoZcFWubiaJcOXXvKJ59bLqKKKrrd2ku9nNPdmFHtbCKWlKQrm9AvfO+KqGjskCrQ3T1hczdS1dak1gVHyJ2peVeidrWGqWr6RqepujrhBYHlRr8JAzzw3wD/IXCc/n4JeKyqWev2OvDaM+0pUzQAWOGy2Oh+x+rIMfCdLIOv3n/OSk1FI1Htqj4H8j2OOe4Py07BrbcEGKNN3i5lacM1YLAKgMHiIHP0eej2rLJFJurAsxceOV9RrbdJOz4jLkq6u0suPlrx8LcCr26Z1Ynx6jwPXq7AQbloqOuW47KjLjra6Hh4viS8vqC8FMpLe8i2hO29yNtzx+70mNPPl5RfTdLIWYXOy2GGihfcttuftpQz/GgSye7WEf5hgV6s4MFjqrwQpKL1mG83+2J90kr48DVN22/uzgi3Oo4ro86CJg38mGvvF3QGNc0U2ye8YHg/ZZIi8kXgAht306nqU8f3fSNDt38f8Jaq/pyI/MDX8fjPAJ8BqFnYJTmApEJbVVI/DEhboFUm2EfyPbXihejIGjavESNRi1yxCNggkKxkSQFJU5OOt0LmnlTQu30xRn7MOFjNa+LpktzI1N9feBu7lw/IOWS9RY8W5svSB8AIO9PkS0xzaduOWNU8+K0z2h8445M3zgDo1LHrCrrouCwDqlBXLfOyo/RDMC6KQIhCsYbF25FiEykvPY8Xwu5bGt5+xfH422pu/NrL3PwlWwHaGzP8LpiKadMOmXqyKshUjmx2yHpL98oJF5+4R7mKLL58iWwawq0lXvMAcjNb06ZNtJXuB/TDhRXsPTg+4vxjnvp4jXORpvM0TUHszElSxwqZnMX3oiSZuPcJLxSeAwf/u1X1/rtt9I0O3f79IvKvADXGwf8IcENEipTFfwR446oHq+pngc8CnJZ31J2epIhsPLd2HfXbW1x7hL7TUaasbS+O9wU3xeVJTr3eXG1/2YcGEm/vBw9y3b/A17F8EhjPVkUEWdTE5dwGU3tBqgJpOqSLqDw5yALnLHsvPJo7XIPaKDsRo0hUicc1lx+tOfttLR87uSQidNERosOJUvnAUW00T+VDb8SlKgQV2+1RoD0qaM8F38DirQ7feM4+NWN7N9C+1PHguz3re8fU95VqpSy/GnAbG9EXbi7p7i4pH5qVr3Qxje6z55q9/pjm9Da7E4e8umTx61v8RfKXicbZS12b9BMQ2l5h9ASS5YI7WtJ86h6r1yKLsiNGR9MWdG1hssgshRwppXrJ5IQJLxiUdOV5DfhGhm7/GeDPAKQM/j9Q1X9dRP4q8AcwJc0fBf76u++M3jALIroLaNPgzjf43TGhZojgh1l8ZmpSUBdN1IvmGqFSVw1tsAwwOIge1DlcHlOX+WFv3adsB6ezJ7zf+0Dvh8cUHrfa2GS6kxntSUUsa/w2Uj40ukUgZZmjn7ZDgkNnBTrzaLnAbZvUAerY3p1z8VHHrZcfUrlAp+4JLq9KGbsbyQOjSt/URRVpbhpvEUvP8Zdbjr+4YXZecflywe6WJ5ZQrGF2oczfaigebVJHrwX69tUlYXZEed5QnCXXzCxr3OyYv7XD3Rw5Se6afcdM5/peAwWkbQe6BobXF6zOcfclHn26JpyaRUMXHaFzxLEzZJ+9pw+CG9021VcnvGB4Dx/J2yLys6O/P5uS4cPd/U0xeuK/v+L+Hs9DB/+ngB8Vkf8M+HngL7zrI2I0y9kEDeaxXqw2+MaoFy3Sb3cQ6Me8/N4+rfFFVTipd3gX+drZMU0OBl6eLPwld0VJskkgeauQ6AO/b7WbF4R0DhSO3a0Zm9sWNOf3I+VDzNgrFypHz9V7pgM4R3eUGpLaSJwXrO8UbO5GXlmYnChE9wSX59LCN749JMfFGB04pTsKqHPE0lGuC+qvXrI42zB7ULO7OaM9tiuZ+n5L9dbloGUvzCvfbwLtSYF0BW5X4NfJDjiZsPnLHbP+rUi+Ol72/WnSwijeocEhMQ5BPi+egLt1k/VHT1m9Jrh5hxOlC+m8o3nOPFFfT8KePTfRnORPwX7CdeO9FVnvvxOnnvDPqeobInIX+L9F5FdU9f+5asP3JcCr6t8G/nb6/xeA73uPjyfudk/evtlQrEFuPcMKOH4BMyUThRiFumh5bXHG/cul3e2woRzB3BJzAVDAgnxdIZepEUoEHSt8xpmpN+13HnS9uz03aeJNo0OWX1Vck7xWiqxGST7nhbe/O7NHxjtiNSfOPMVFQ5gXNKdCvNHi0J6eyTVEJ4qTbMBlxxRV9l4nVaBN982VpohcfsRx+msF7uyM4uyS4s2SeDynO52bP0yyL9Yke0SV2ddWxOqYWDnCokJ2c9x6a9slTt5f7oyW2jXDgBKXjjfr98fWyM6l1zsO0tKqonvtFucfK9ndDpSlXZ2E4NJizUDRDJX3kWrq4Es0BfcJLwrex8+iqr6Rfr8lIj+OxdvnF+Dfd2T729Wa6lxZd6AzbJYqgOzH86cV1DQKXet5vJ3zqeP7zMoOSiUWEBYVLgfuPOgjGWVpavcfNyP1g0DGXjCpESoe14R5ye6Wtekv3lSOvtJSf/kM2bXE4zmIdb6SO3ZTo5RW1j0aZyUuRPx5g4SALEvUgysDTfQUiV9vg8eJUrjY0zJRpc/kMzIf71ae6ix173b22sm2Q1dru0pZO+RyTdneoLl3bDJJurR4CbJpcSEyfyO5SZaecFoTFyX+Yms00zpNwxpx69J2tvj1hmOeuEtNWSJIPUuDxE1FpKrw2j3OPmUeOboIVFVHFx3triA2PnnPSBqoPX6j8+diiugTXky8XzJJEVkCTlUv0v9/L/CfPm37Fy/Aj7Ll2LTUDwOuKWB5qIZ5kpPfQ6JuROBsNcfdjZzWWx7VS7plyeZuRfWFwa99eNI4WPKmAN/TOIUf5I4xQpcy/MJRXOw4fbBieXMBqhSP1si2Mf59tU3HLEOW6Z15y6duULdrkQfBqIvCfGNcC3HnaUaGPG3ncWmyUaZlxsE9Z/delPmspfyK4+7P7aj/yX20Kjn77tu0txfMHtcW5NM4Pdns6JY3cDcX+LOtNYJFjGKJZg3s0jlo4Qk3l8Sj2jztk7dML/XMxdUQzeSsKuHWKbK1uasK0LheOaVdhzs94dF3v8TZJxztzY7Z8Y7SBx4/XqI7D53YdCYFLQWKA//8VFTevw0zI5sw4RqhJAfU9wf3gB9PMakA/jdV/RtP2/jFC/C54Jay+NmDHX5X0DcxpXF9AuxN9RkjyebER8qq43ixZdXNaIPHF4FmGdncdpxWZa9FP+zU1JOlZadtN9Kyp2gx6vyk7SjePLNtvMPNS8tcc7QZz3aFPgjGRWV2CQEkDOegzqGzkjDzxlEH6YumhYvU1eDZE1Twok/w8lGFwkUzVyvAhYieX0BUjv/JnPVrc6q6gu1gHyBtR1c74u2Kxa5DzprBCjmdW+9d3wXctu3PX2eVDT8JCrpDFzOzB267FOQDlAXu9ISYC7RgdRfnkNmMi+/7KGefcOxeisiiY1E3bHYVurb3XrJz5Mgd4on3fI+me/eP2oQJ3xT0qr73YVdGgX/3s27/4gT4FND3oErxeE2xWeK6FNijPPnlHTU9DZ2sIM48xG/WG5bFDhEbwK11ZHvLE5dz3KPzoaAKlnnuGmJdQjmDxlvRMTs8pm2A3sMG6D3dpYtD1phMwnDFoBgpPHFeEktnuykcmvzi8YI6R1gUrO8W7G6CzCJdomUAvGhfRM1/j3t79rJ5UbqFEmaeUgTdbfFfe4R+ZIEua1zboX5ngbzw+EbZnTrKkxmzTZssBkbeM6n+oE7M2z29RxQOLb2dg3P9gqWFN+OzNEuWskAWc5sNm66EpCrpvuOjPPp0QXND0UVgVreIKLtV1WftCGiZ9O7PSMXIpIWf8ILgN7MXzXOFnK8oVreRDmKdmpnGEQ3obQl0/zZVoes8IsrNct1ntZSR9lQJpzXu/iOTQma3x+R+KJsGqtJ06vOZyS+TGZl0HejIx9iNvc874+tzYAwRnVemc9ekLnFiCwEYz53WCC0ccVbQnJSs7zmamxFXRrZtQdN5FChcxDsrrh4OoD4M7qWPhLkSqnRFFCN6do7fvUJ3o6ZsgwV3AOfwu4hvUiNR6fc7WBP2rk7SFU+sTM/vOhuY3Y/q81ZfUDpbJMsCqWtUHNI2UFboa3d4+J1zdreVsIi4eUdVdezaEt36vnlNvZqxWDHi368K4H1H6xTcJ7xA+FAH+Ce400zTCHpxSbVSfOMI80MenqtfuEzlBGibgqCO28UlM5+VMdAtIrtbFeUX0pNnJU0yEpPLtVExJ0ujU6oCv92hbYvGRBHloRv5aUO0AmX6f1byaJVUJSmoS7SGIS1s9Jyo9vr99rhkddezvaPE447SRTPY8pH1aoZzSlEGqrJjVpqMcPzyjYN8XXSEZSTUghS2QMTdjnLVsbtR4jYzim2yPVDF7wL1faVYt0OAdokaGev3cy0h001VWtySTbGO6CnNhmzZDbNMbpzLOeH2CY+/44jzT5iHjtaB2cy075t1hbRD9p6dI8Xtv+F9ZqTJC36iZia8cJBr86K5nvaqA4g4pCz2CqwZcbOlOg/4LYMV8HjY8ghmh5u+5ElSF4PwUr3iE7M3KSTSth7tBF0EVi8nZYumZqZolr+9f0rXwa5BvbB9dYHePNnTbAMWBLtgo/xS8JOoA72RaI1YemJdoD5x/H7Yh3or1AJ0S8/6FaFbRHCKLyK3j1d87ytfZr5o6M4rtg/mXF7UrLYVTVcMp5v17+n1KX3An7Rsbzr0eNn765Rvr+jmjrAo7cojZ90RyvMd/uEKd7ExWmWsjolp1mp2t0y1BFHwm9aKybljuCrtfPNM12zpECPx5hGX/9TLfPV3nfDgu4RQK3EWKZctVdnRdp54WaaKMXvvZy+XTLTcYcDvPx/5MVMiP+FFgD7jz/uMFyODd4KbzSxRbLs9Ll67ltnDHeWlSRD3FsKRFHov08uZfTQ65MsXN7i4Pafy3VDrrCLdXPYDGNa8I2WNLud9xuofrZkntYgUHu0CGoNRMWWBJg2/5AHb0Ge3eIfsWsRXxMIhzoKPy7SMYDLJumR3q+LRt3naI6s1FHXHR289ovSBoMLvfO2L/FzxEc7OFsTOsduVzMqul0nmIA9AdLTBU88bzj5ZU17e4dauId5/CG89QOItuqWnWNb4RxdGtZQOCR5XFmkAduw1/Djzet+TqjQtbtcgy5pe9y6CLmaDEsmBagFhBk5YfeomD7+9MEqmSjRVofjTluVix64t2F7OjHvP63iR6Bmnqb4hUEa7SjukYyR9fMa1minIT7hOKDbD4BrwYgR4gHmNc4LKzpqeMk0TA/7hivJiaTp4seCn4+wsB+0DIQUKdI63Hx3zE8ffw9ubI3uhO4c2jmKlQ5G0KGxK06yybtb1dm/6kHQ2tUmzKZpLQ7WPatxma37yTWtBHwY6o7OFQILiUvVVReiOSprTgu2ppz0W2iNoTpQwj0hntEi7LvnaxTEvH1/w5vqEedFyNGvY1BW7tS14ZkkgT6gBQ5pgFaPQ3W158NtKuvoj3Pl7C+TxBa41bj4uSvz92HfkqkvzW12SOXoHfmhc2pv4lCdHbdvhKiB3+karR2hhHbqbTx1z8RGjnkJtXazqQMuIO2k5OV7TdBbcdZ24dwGt7EqGvdm84/f4oOh+RcCfMOH68WEO8DFZBtQ1iLW3azca4bfaUK0irnHE2T7Puvdd3rvUyZkshM7xj956hRCcte8XEdl5fB7gVJZm3+sHp0RtR3NRxRnXPi4uqnV9ukujJfLiYM8Zh0Jj4S1gqhLLgva4YPNSwfqe0NwwlUusA5SKlDaCjsZZBts4Ls7mzKuW24sV21BQeuPfu+TZHqKjDQ7vhs7W7IkO5ijpqkB7w3H2bR71tzj66gnrO87UpF1F8aZ1rbo2IkF7Lh2wYSRlYUZoWSqaF72sd8/e+oB6z+7Ogs3dklAJoYJuLrQn0C6VOFO0sB8KxdUd84XVATarKhVWZRi2njL3YdZqfqMPIEN5YMKEFw4f5iKrJq91mc2gKhPVEdEUQHS9oVjFVGjlSgvg/tL8kK6JggbH6qIeinB5ap7HmnB6T/jYSx/3JjLBYFkQsyukmgUu7I/q68LeIBAtPWE5o7lRsX3Js7nt2NxT2pda/LKjLDu8j7ikjAnB0bae0Hm08WgUVruK24sVQc1J0ruILwKh830DhY6kopJsDCJJdeMjcRZpWrXMSwAAIABJREFUb8LZpz3bl0rCzAzGqgtnlArgt91w5UFSzPh3LtMMtsuW4YfTOWefrLj8KDZ0pI5QJeloNgtzQGEKoWrWURUdm11lVsBlREt738izV99NMQVJIz/i3vPxvc8BX0R+CHNN9cD/oKp/7v19hgkfSHyYAzxgMztTAVOqin5IRjBnyfKyw29KuBWhO/jWHlA1hwZk2qUo4cfffIgV6KIeGpp2DdoFC9Z5ViuYDDJEC+7QF2TBKA0AbVu7EmjbPuBrVdLdnLO+W7F6xbN+WWnvtNQ3thzXO0pvHaljDXtITU0hOjZNSdtaEF+1FfOipQ0m+/TJBtmCuVkT5IUvK2ucKIWPOB8Rb4Ez3AusTj3SCfKmJ5RCrEvcrkM2LaRib7YZsI7U0P+/LzCLoM4KrepNIx8WJavXas4/Cfraljs3L7i7vKSQwOsXNzhf1bZwpbfK+Zg8621RK+uOatYhoqwuanQ7kqKO5vD2N/UKGoa5rD4+Sdu8T0jzh/888HuwYTY/IyI/oaq/9P4/24QPDA5pw28iXogALyJo06LdubWvJ6dGVxTWEIOpNMrVnI2wN2h77Etzlctg5uEBtGXE5UKYCc3Lx0inlA/XyMOzdBlgfvQy0nxrSHNFZzPj3zN3n60LoOcIpPAwr2k+coOH3zHj4mPQ3d2xPN3y0nxLmawGDqVTIopTofKBwkWOZjua4AnRcbap+cTdBzzyc3s9nBKi0AVP4QNl0sSPu1r7hSNRHCrgqoCbdbTbAnVp5F7hkAujo7TDznnskgnovBq861WJlblumnbf05wUrG97Lj8qdCfWifry8oJXF2cUEjgqd7xZn/DW5RGbXdlfecQodNFzcrThYzceclTueGtzzOfWM8KoqamvXXMQvGN6f2N6UzNSJi/vr5Lm+4DPpW5CRORHgR8GpgA/4R3x4W50ciaT1KYhtp0pTcT13aF4jzQd1bkOX1rdD+xXjvVLd7hdzkoVgqCiEIV2CeffOmPxdkfxyLzkpSzt+crEuTctukv2ud7b/SGgGgc+Og0L0bZF5jXh9gmX37rk8ac9q081LG9teKneURdd35x01QivnJFnYzEn2gfuzW7Og92S7zz5Km9WJ3zx8hbbrmB7xcs5VtO41L2bZ5RoEIp5RwgO12JjEaP2xdF+GHm2byg86hzdSU2YDeZs3cKxuenZviQ0p0p3pMR5gCLi54G6anl7s+SsqZkXLffmF3z76ZvcmV/y9uaIty6PuEgZ/bfee8CnTt7mpNiyCRWXxYzlcsv5dmnrqyhJTT/w7Ok29QoSjdLJK3yU59XF+hrw5dHfrwPff7jReFqZx/8zi34e/TNAQHzRf/a7uVlP90Xnw9MaK8lGv/c2e1pwOaxXg6mePIMkefz49LrjwBeWpJQuUEnAS6SQQCERwQQF+V2zr6sQMEFA/juq3RbU0UZHp542mJ10viI7tJzK330J6f9jB+7DU73qIzB6vSTauRKh2MT0+TeKth/l+Yy44NF9Vb3z1A0+3Coa/f/Ze9cY27Ltru835pxrrf2oXVXn0aeft++9+F4TG4PNQwEJkjiQECAozgMhEwVIQkI+YEICUWL4AhJCIhJ5EIkgmeDEKCRgOSAs5EAAOUEIYXyvbXxt7Gv7dve93ae7z7uqdu3Hesw58mHMtfbadep0n+4+fc+x3EMqVdXea++11lxrjTnmf/zHf1gJO6Db2rB3IhqjOXvv8edbyvOUHbRdHLGP7jXllgh0PX5uAlUSIQX28HktEtsbENaOxZv5Ys6m1qxCFZ2U+fv6Rtw26ehms/s/66GjpinDlSO2r17h9DMl55+G+iVz7gfZuctIN6bHyC+zfpuUHb5gu/na6TFXqjUOZVHUFC7S5gh/rAvfrwz6/YnsnpTU2UPmQ7TCrgTSRoviu/ER7bD4dFCiXkiVY3M1sLkhbK8r7ZWIzDqcN8zfYCIoisi2KdjUJc4pZeg4bypuzJZ8w8FdmuR5Z7kgtY5QRr7p6BbPlUsKidQpEFzihcWS5cksH8poVdInU8dPr9P9h/ljxOAfx8bdyg7lqv56+S2P3lgEN53iXrhBvDKnuTZlfSOwec7Rmro1xdImYl+D6xQXzQEnL+aMBTRALGTPOfeTwng8xpXg1hrSXo+loMGek3ahdNPxMfZwJKRZQqYdk3nDwbTmxvyc5ybnTH3L1WLFUVhTSOS5sCTmwS8lsko5zyOJbSqIOJI67ndz7jQL7jZzTuopJ9sp66Zgs66IrbOke2MFgVoo0gl+5fC1taIsVvkQ+3HJvsHX4Fs7P4mKr3Mgk3NwsXSkQugqIU6gWdhvSVCeKIubkcmdLeHWKenWHdK23kG0l9jf1x/46nvdE09L6PQjOXgROQb+F+BbsKH7j4EvA38d+AzwBvC7VfXBe37RSNNFvB9gGSBXjUaoG4rzhGv9nkrAQ8ekFqXiLNLQSSKW7DuAJLk6E4pzxa9Nx7zvVKRti5xhf6ed0xPvLKLtIuO+rFIVpBtXOP+GI85f9qyfV5rrkeKgYVY1g3N32eH2ztjJw0JhwyFmRy2iVKGjKCKrTcUby6vMi2Zo2acqJmOgsiddEJPLWvFGlRw8XRTa1hQpu5nSzkaYenCDpEIvM5CmBc1Ryfr5gvXz5ti7qy1h3lKFaNo+mY/e/+46N+xTVagpOF9XrJqCZTPh7nrGel2hrSM6ZRMLJq5lmwpWXcU2FkxDS6g6ujrslRYM98VlSVaw1OeFrl9PEI+/CXxq9P8jW1K+l0kIuNkMuXaF9oVjNsclzZGnnQrdXOhm0CystWSxFPr5WUdCa67DWE8pT9JDS0ohBSF5SMEcfgq7z/bOXzE6bP8sSY5TfA0SBddCe5AbpauQpgkOWw4WWwofqYqOyke65HhrdUybFU8LH03sThJH1YaXpqfcKJZcDecUEonqWPgt21RwNy540M04aaecNROcKIuqZl42bCc1TRdYbctBQ6rrPNvTipgEFyEVFuwVK0WiTX6hNmfu2txesl8B+BwYZp0kVxu8GAohloLf2oTaHgqrT8Hys56wOqB6cMD03otM73ZU754j79xDVyvSZvP4uMuTvQc/kH3UCP7PA39HVX+XiJTADPgTwD9Q1T8rIt8NfDfW5emRpqoWuffR8MX3kzFWwqbDtYEY9kfrIv9dclJDi0h50Bg9snMDhZAkaOtw58LkNOGXtWmldLk5hzgrXorJ+PHeAW7IDQzqkoBMKuIrz3H6uTmrFx31FaU9MkXEyaTdg2VgX0pg7NzH1MZxFA4QJBnTZF1xsp6yLYPJARctwVtjbbnke825u8Hx9sU/qfFQCDoxcbE4KwmnG3PyuR8s2EMRD0rWNwqWnxY2L3aE44bFzAq72jbQdY6x2Bsqw2vD+SRBu8DKVWzqknpb2DGokFrPa8trXK/O8STqZM3ESxc5Wmy4Vy9GY/TwvdNH9P2+P+aQ/UeBz4vIZzHH/p3Av//YnxbBHRzgrl8lXj+kuVLRzj3tTOhmMtBKUwBcjqYPFdcKrgEfzYH3Xc0kw1LWGSuft5oTCy27RiuQnbs5fvVCCoo6IRbkSUCMXJZsP5LvlzhT2kXi6FOnvHh4RnAGxTQpsGpLlk2FsKu7qKOny87+wXbK186uUPrINLRcmaz59Ow+r5QPWKeS282CVVdRuY6Jb+nURPXqGKzHcL7Hk4pJR6+n1GUBW3vgiyWUZ4pvLHrvdZViKbY6DT5rIzFQeum7iDm7cVxLpgaD62wyjaXQPh9JpZEH2gPH9mpJdeOY6tVDJne3hHcekO6fkM7PH8PRf+z35SPtQzt4ETkC/mXgPwRQ1QZoROQ7gG/Pm30f1unpPR08mEbKZVIFg7UNbtPhGqGb6eXw2nhp7hRXRg7nWyvp1wJE8RbuQCXo/UB1v8Ott2hjXGzxHqlKw93J9MfCiorIxU79Uk3KkvTcMSe/fM7ZZx3dXImloqXii0QVIoWPl7bVe+jYh2Si7L3Wf7bw9rBum4KmswdoXRR4pzxK56JvBpJyiX/v4LXJJVeV0hzD9nrJwalBT1q47ESEVHnqKwXnrwibT7XMrq2ZVS1JYVOXtE0gxQwjudE55P24oBkr313X7apEG7+LaKLw7skhr0+vcaM6x0nKzUwSLy9OeXA2I9b+oQdEnF76XA35brnQjvEJmKp2IvJdwN/F1grfq6o//bifd1WFu3aFeG1Bd1CQCrGf0DteBkjF1ZYrSqWyeUHppo7yzKJVFcG3mqugBdfpHiwzJJZH8IxmCQ1r4avZ5+QIPzt8MEefQoY6EboF6DQNUOBhuaWJntJ1tM4a0bTRo8mBaGaACVXR4UWpu8B557m/mvHVe1f4QvcqziXm04bj2YZFWeNQgouU0rHuDBptoud0Pc2ieZGqEIoQWRxuaKYt22mFa0tcm483gagMkJNExbcZhop5haO2OuknAsmEK8Ted3nlkgrQaYStJ04UiWKQVXKkIHTTKdNpoJxU+FuOeHZuEO173jyPe5c8WfsoEfxngTvA/yoi3wp8EfgjwPOq+k7e5l1MoP4jm7Ydbtvit2TMdRT9XkwWCeBNxyX4aBHjJoBXi+onHUfzNXeKqQlr1Y1BLT1rpghI16FkMbG+BV1KaNPatrMp8YUrnH7jgpNvFLqpjh4sozFWobsUhrmsA9N7WSLTIV2iaz2qAU2wpSQ1HldGJrOGqugIfn8FlFRI0aGZaSJRLI+Ro7/mSDl/2TO9VeHXLVp4UuGIlaNdeM4+7Vl/pmV+fc20bGk6z2Zb0NV260hO4PYQzSAF4fJ5gjUoKZSy6Gg2xYWLBc028Av3r9MdO16YLjkuNiy7CXPfcPPwiPsP5qTWD9uDTSJ7OjQXseaPKWJS1R8CfuhDfbbr0NUGN5/iqoCUucI5gosmjicOXMTOp7b8kRbQXEm0h+BrobonlGdY9NlDNzCCa8yxJS8XkpQ5N6O7aL6HbWxSt1VDKiyK7SYgLfiTwIoZr9eBd6cLqqLjaLLFocxDw2zSsO5Klm1F3QVqDQNs2CVn0Mq6RM+DFbE1wlJnnIZjdBaRMlJO22FlWHeeGB3BJ67MNkyCVSRGdaRK2HaBB6JsXnTEyhPWQtiAy/2bbUyF2PVReZYZxyYuFxXXmeO3nIUQC6FZCMvPQPNiY5NBa+Nv8I99DvIKyH/A++uD5WyfmH0UBx+AXwP8YVX9ERH58xgcM5iqqjxUlWQ2ZhpMmD28QS9V0FuMSN3hGwxbl4ehGehnZB0cbdN5gwM6gVZIjaN1UB5GtNJBAIuUmTGqxuiJJkCm261NAD3/u+uQowWrb3qek88XrF6y6kzJevW9QJb3acAjL9oYg4eHIRm58H6Xk6hJxc6lP9dOkI3x5GtvkXwR4h4kZHIFzjoidbZUFMWOVyBVyuoVqE6mLN7ID0GXiEeBs097lt/YMrm6pQwdy3VFsy3QziE+IVkM7TJLPVwiSghWbOVEh+h+R4ZXUus5X1eczqb8+itv8K/Mf5YGz/fd/k38qutv8+PpZU7PZqTa70elY+YM7FFg7USe2sr4UtOuI969i1suKa5dJRwv8M/NaRcBiQ6ZQuoMWqNkiDg1YnLZwYrH2sOERKE4dRTnUJybA5Iu3/oRg9tGz8iAwYs58D6C76GgFMzZDUQkn7cj6yadBLqtY7kIrMrI+aYa7vHeRJT1tqJrPbHxaO2QrTfiTWu4uUSG6286fZ40hzo66vtT/LlDktAtIuGo4WBSczU0HBQ1TTJ3dW87p53WuBvKelbRngf8yuMb8Fu7z12XjztDL+RJz9c2mUq0ZzWWQnsgbG5Y3kxmHbQOWXv8VggrS+b2Se7qLDG501C+cYd0/wFpvX5/iCZDxk/DPoqDfwt4S1V/JP//A5iDvyUiL6rqOyLyInD7sg9fZBo8tMEFuEZVcU1r+jGPYwopOcMDxz4oV4xu2gJ/7nB1t3eBxOVG2j0cM7ayQF99gbvfdsz5K0J7pKRCM1WvX9b2+370Bb2s1d6w/0dg9DE5NMpD95JOElJG3Ih+qckaWqeMh6dOBgevQQ0KiDKwLGKlLF91TO4XTO5sqa9WLF8OnP+yiF+YwuOmLmnrYAlsUcTrjl/f5zwyBW5w5KNzKnzkfDUZahIG1EYw3Z1VybvFgi9OXuVKWPFrJ2+wiiVvnx/x0uEZIsrp2ZyY9erF9ZNFTpgngWDe7NKmMM+KqZK2W9LNt5HbJeXNKdXBAfGFK9TXpzQLRztzdDNzPqkwR9uzwVzG4VNhjr65ahCCq80Z+WYXcQIDw6Z34OrZTa7AQDnM3z3UD8gu4aoepFD81pE2JXGS2EyKHXuprzrGot7+mroo+H4Vkp+PIT8WlFSAawS/DbhW8LVp0sWpgldi43n3zhG3/aHdQ6VVfXtRZlXD1fmaqmhZlhPiwtOq0HQO2XikFlyTHX0rGWuXIWpPQWyyPI7ILOLLiItCWplfKJaOsIFwDsVamd3pmLy7xr99j3R6RvdBkqzsxyFfT/vQDl5V3xWRN0Xkl6vql4HfghV8/HPg9wN/Nv/+W0/kSHOkHbbYAz1+hnNUtyuAMrw59c4kVzaKCopFuts27B6CXk0yeOPkb2qL5nvzHjmYU3/uBnd/5YTtdTWJ25CxTwc6JuZHoWsDbfRDkmh3Gg/DNWNaox3PvnNUFdo4OpfxeTtztr31ydqkQozOKkc7i4okO8IUUlZV3i3hu7ly9mqgPZizet6xekVh0VJVreGqeSxdSEPzLRk58t65S38MI/gkZWXLHtbZozX2EXkS1mcTfkpfBOCF50+4Xq74yeVLNNFzOKkpQ+T+2Zw2t3C0L8cmjQtO6xeDadsQTxo4OYWbb1N5z3Q2gxvX6J4/oj0I1MeB5lBoZ+bsNfRQSqZEBiVWSncY8Z/dsFisKXJD9nVbsGkK6m1BzPRY3fqhEXtYQbFUWxUrw4Sf+kVijnBTZUWB6iC1EM6twUsqIRWWJyDfXwMNM1sqFK166CTfKy7ha6E4FcI2Q03RnHsq7PP+zNvMJKaSGovcEdJbQHUCFiQ4xR+0TOcNziVmuaXluilomkCKjlbFJoZcOV36yOGkHgoIV+uK7qzErxzT+47pXWV+q6M4i1RvPkDfvkXabNEUGfV++4AX+8N+8GHLFdVfAG6q6u98r20/KovmDwN/NTNoXgP+I2zYv19E/gDwVeB3P9Y3vVeCFQY9cdfp5YM19nuKRRTpAtadHYkm6Dqfl64W7e6aWySjQYLp4hwu6K4vWL88Mzz6hXxD95OJ6LAE22GhQqw9mzZQhd0t8V7O/VHWJ0rbaJjmftHJaP/0E4PuVCQ7R4q76KqPrERtaY6oOepozmP9EmyedzTHiXTUZZ2YyHpbkpLbu0TidQ8jkwt42fBXHv8Y3RD97zl36fMpClHYrkpunh/R3gh8y/wtflg+x+lqihwoVyYbZkXL3fM5q3VlyVevuKollBFNQnte7r5bd/t/5i3rMcWzMzg7Q74ilOKYXjlCjg7RSUl3PGX18oR27mjnQjeH1SuJG5+7x2eO7vPp2X1eLE94IZxy7Ne8FE5ZppJ/vP48P/Lgs/zcveeov3bM7F2hPDX2iXHFdXDwfXQbM31QvZJS7jrmbW62+gl7zXkhVZrx/HwPOM1KoYpWWQ00CbJ1xmFvZNhXf8uo5WgJK8sz9PkDk6juqZ75XtG8/7zaUB/YzifEeaK7vubF4zNeOjjluck514oVM9/wanmXG2HJVgtutld4pznmjfU1fvz8ZdzXprz0TyKzry5x5xs4WRLv3gXV982dPiX7I8DPwPtX0H0kB6+qPwH8ukveeo/Kjg9nMkjYcgnwzs7h9fhuMoxWRO0GG80Amrnh6kELZ7BMry2TnbssFqRrh6xfmnP+UmD1stBcSRf06HUPE95zXq2w3lZMi24Pp9w77Etw996hj/H5Phrfq4YbLV80wzLqoj2MYlFzSkYPHdgzubBJEztYCYN9UqHEuaJlwk2Md+99InjrKKWJYaWgSfDu4dXGxcth58Ugorb/xu597f9X++6mCxTS8fnqXaqi42w5Y1WXXJuu+WWLe0xDy1fa66TOEYrI1aMV33jlNjdXx7y2vmG9e58l8P3DmibS6RnStLijQ0IXOVy3bJ+fsXylYP0CuGsNnzm6zy8/uMWR33Ds1wDMXM0ylSzTlJlrWBQ1MboMmWhOPOpQMDVg9d4ce/IZLpnklcPIGctocduzVyjVGEFVQkuFkGxl1dpzizdn3xWKX3qKlcEmfeWsa7EG9ALS7vbRryaKmr0gqs8nGPNIkBZQq8FYNSWFM5rlr5q9yUvFAybS4lHuxTkvFw+oU8HPpRvU25LyXAjnEWkjst7aJPuEtQWeVJwhIq8A/ybwZ4A/+n7bPyOVrPBQUnVk4v0g5esu+Mq9wFFHv3X0/x7bYocb2/LS77jtYA09jhesPnPM6WcLtjegOUqkKhoNq7Xk1R7MoKP/+9eSUG8KlkXFrGqGZfPAxb9sCB6BzTfR0za7gp+9pI1YdJSSINHR69en6Eit22HefaJNe1xVc7LaepwmMdwTtWpXc+bJEn5DmNVTIjVH/26nYpm59t4nUnS76Lwf9vE5DU7+wuoqY+vBW+n7KlUEbyJqXec5qycclxtmobEEX9VxON/yrddv8u9d/QI/VH4rr735nDmcXySB+2AiWYfJIWWBu3JMOl4Qjya0i4L62AqhYmUwTQo5ufrzU/7p/c/zY9c/xaeee8BnDu5zrzaJiPO6Yl0XuQ5CiJ2j/WxDfS0QlkKxcsavrxnGPwVIuTCwd64Ddj+6Fex1He6rvpuaawSNQprY9uRtcGoy2K1F+CmAdLKj7wu7iSQ7/THck9g/Fst95ePNRVjVpKUqbcV8Wk/4se0r/Pjdlwftp6NyS+mNc397s+C8qagmDfU3wesvFZQPrjC5c5Xp3VeZv9NSvXWKvvUOWtdDk/gPZRk5eEy7LiJfGP3/PTlf2dv/CPzXwILHsGfHwY+d+wVnL70mTRcNohlH45eNm+x+JCcEtcfJ/S5BmDLKQPAwXRCvHbB6ecrZp3NP1CINEYOr3XCzSX+jDw4y69v0++6PrXFsNiUxOaqipQw7UbAxfbJ36BcJR6pCmxzruqTdhh3UcsE0J7kizsTD+tZ2gwjX6LgcQ4Srmh/QvltS30Up5oe0XwFl/meKeVDVrlco4m7lkcdFnHGaY2dSx5oZRSEoNZXto9d4vxDNh2nHteNzfu31t5hIyw+ffzOrbUlZtVyZb5gVDV++f8Nok53DV5Gm87yxvMZfT7+eH7/98g6P/8Xg4EWQUOCmE2Q2Ra8c0l2d0y4K2gNHO3V008xNLwySTkVmuJQMUGFYOfzdObe/OOdu+4rJFpTmCLvZKPBQCDnqjVMlzoxL7xsZVlGp2K1Kh0g9T+KDky/U4JL+lu8nVGF4vU/oD0nvzSjBWSpporSANDkZmiPwobJ2D7PXwbHHaUILhSrhSgvKtPEQhfqsom6nnNc2keyxdoCbefKI04QuOqq51b44H9FDaLzSHng2N4Szz1SEzQ2K8+dYvNUxvblCbt4hnZ1ZzcwHdfaPv/ldVb0MFUFEfidwW1W/KCLf/jhf9uw4+EdE8K4srNAoc9Ct6mwMuPeff8TX5qhTISdlMhSdIZq7v3KK++aXaedCu8ga5j2+DGOI+6Fy8X7/w81+0f8moWu8YeitpygiZeiYlq2JgI1gGMe+g1c1vu/5tmJ9XkHjjCWScfThnPPDpWpkdHVqWPfYuY/HI0fue1TF8fj1idjO0eVVk/eJLssU26RpD0VsvfG2nUkSO2f8fzt+SLkApus806ohzDq682KIZnQMbznl1Rv3+f2v/GP+3YO3iCh/6d3nKIMJl51uJrx96xhdWz2DlInUOlarCa+tK34hPUfaemuaEk2D6MM6ehH5FPBXsBoOxaKoPy8ifwr4T7H6D4A/kXnxH3QHuKpCFgs4XpCO5zRHJd3c000c3URI5YhF0zv2YueAw1ooz4TiXAkbkJh2sgxisIs6Rp+XITLuud/q+klDaRdqkfAYZuyfs7w6ddPO6h5guM52OppXCJ7UOrR1+5F+SDBLSJEoCpO3AIhRaDcFugyWC+vZPzmap8xEhioHEf0qE6BxpK0zbZpGcJ0xifyWPYrkZYlf9Z448cSyojlW4lGHTCIEg5Gi9uKEVojWzQqWLx8xvb9g/rU1/qu30OW56dP0kiXvd8mfTMDxG4F/S0R+BzABDkXkf1fV/+BRH3h2HPwlJkVAytIKRNp8ES6DY/Y+NPp9EQIYb6tCKhPnnzEa5S5K0Ev59QNrp3foQ5Qjw1KVUVQD9r9GqxrVaA9A4wJ1GwbHNZb59fl42+hNcKkuqTeFObXLTnbsINX2MejsXNAN2+eG646mOHzeTljUeNeIQ0Ni2xTMJg2rVO247UBsLUIXn4xnnzB9+H43spuwYhS2TcHBfMtpdgD7yWKlmLYclVZN+1oHX6pf4ayZcOPgnLfPDlmvTLuGIiFBh7FOnaDRZ6wXcMnwXt1FgB/COuCPqeqPicgC+KKI/L383v+gqn/uw34xzuMPD+C5a6SjGd2ipJt5uqkz4atSMpvEnLtBMjbsroPyVCjWit9aghS1AhwZler3sgQ9pm44uw4BSu/8+wkgTu3/ThypSsMzZM+bEspIWZnkRlV0TILllXr9GSdK2zOloqNuQ84BCbNJw2JSD3mlXia7S44melZFZJ2mVsnrdJiYtbQyU7d1SN1XqO7EAyWxo0E25ElOrR4g0zGlLxhTdpWsPeTjbbUab0N7ULB5LtBejQbBessn7BITQBI2Vx2xnDOfv0L12h24f2ISK137/oHEE3DwqvrHgT8OkCP4/+q9nDs8Sw5+j6LhTDIg68FfbMT9vnS47ITF6ZAM3PscBiVIVqakV+Bz+uiZ9sLkMDA08uf3mon0++oflExlFDLs0gTa1tN0npCbfvi+o1MS6jZT2xqPNjnxicfpAAAgAElEQVQS97rvFMe0zPEx9cea2Qx7r/Und3Gh1E8ALkf4CkTTianrwiCWImYRMUvcahSjTPYfV0vsaitIdvQ9vCNiya+qgGrWUq+tWKrftziYTlqCS7xe3+Bud8iXzl+mdJFlDNR1sO2H3IfuVgHjlUj/t9cdN/5DWK7Efif/vRSRn8Gkgj+yiffIZEKalsRJIBVuYK3soBjJFET7TFhjkWlrlEbfmJhWz1vvy/H39pMvvGbNGQ+j6F6xhi0G98QaK3oSaL2YauMkMj/acGW2wbvENBjdt/SR0lmj9049y6Zi25kbmRUtKQhXZhuCWCV33e1cjIgO0fFxZbmUdVdyejhh1ZRsMrWxWZWG1+fexK4zxo5EY9j4hpw70CHq902vHHlhHLJi5rC6scGxgM7ZBOFrY++sm0B9Ne6goQxZubYPfuz7kxe0LKypz+Ne+CcTwX9ge3Yc/MjEiemxO4c2zU5d8lHaIhccfs+HF29StXvb5ERnUUS6ViiWMix91clQJYt/2Nn3kMLu9SxqFNnh2KNlrYSEC4oPxkgZlrLREaOjqR3tyAmCJUcHxx5lgKN2jmwXpe1POqPffeJz/F4Pn/cO8sIyXDR3aPIKfYTdOKJ4al9kVo0aDCRKrwak/Zhi9FNTnnFD9D6cV3I0nTe5Ah9slTBMLMpiYmXqb22vUCfP15ZXuT4958F6aqsFZcf3113Bl/Tn1J+L2qpCo7dJ/iPqcIvIZ4BfDfwItkT+LhH5fRgP+Y+9r1LqBdOuJZ0tcSEQRIzFVTpTR+xPr6/UjlCse4fOoJg4yOD28rfwUKTavyYNIGq8ebcbC3UmLxwzxz2VDNTG6fU1n7pywjS0OEkkdcPvJnk23ZQ2+iEqP6q2HJZbU5F0aY8OPPEtpYu4DOgfFxs69dTRXM82But74CO1BMqyM7prKlCcrahXzoqmMlYvGYKRhDn7NuvOxN6Z63DP9xj8oLY5KGgKqOLyb++E8gHEyplTz9XDrjWufnmmVGeJ6e2G4vYS7twjrTY5en9v7z2ksZ6gqer/i+l8vac9cw5enCAhICGQ6tqi9950FKlc5tz2vshwwlnR2uLAmQMWp4QicjzbsNUF1QOGm1yzrGoqlVjBmDUwXKSx35SM77Vibet6Jw9ISJTTlmnVDkyQPrEaMxWwbT1dGwy37PnqKS9DxxFo75Aln5i7cML9cV2WcdbxQ635GNktP/OkpJlOKl5RkpVrR0vUtltj8BRFxJUJcZ7oINYeSTaJmRCoFTgNOvT5moExcxoKZNIOLBwyddKHxPWpCXt3OcLzLrGNBWfn06wEyvCU6AWnLXlCHs27O72ij/BgicgB8H8B/4WqnonIXwT+dB7tPw38d5hE9sXPPVqGQ5W0XqNNg5yeUZ1eITx3RHOlQhYeFZeLeezYXZude8pQzIW8Sv//QB+Mpr3S48exh3ry751ssDFxgAEKSqXdHzE6TusJnTqOyg1OEst2Qhs9s6JhGlrToAkNwUUOQkPpOrrkCS4ycw1VriJs1TNxLafdlFY9m1iQ1LGJBeddxbKpOK+r3J3MWje2mwIah2scfiNZfmD/XLspMDWHb5WqmrH4XNnbT4TxwiTYP7Oi4HLzEb8rsupVJV0tFCuo7ivT+4nprZpw6xTunZDOV4/l2Pfsl3bDD+iTrBIMd7/o3CVHH0PUMopg+9lx8GXZkTmfOCo3LA42rMNOV/pgWvPLju5y68VD9KszK/poGYqomkNhc8PRzhWKnoGzu8EGx99HTRGkduZAi4SvIrNZzWJS48X6rA6MGcB5pQw1sXSDpG8bPV1nCUlNxmLR6HaVmhetd8YD5n6JM+vfg12eoI+Ce5iHPtq3JfrgPEPaTVpqybCuCRRVR1W1pCKySaXlGPrLlMdk4Lz7fT66Yu85UTvcvBQuq5ZZMEZDmR1D4SI3T4+I27AXAo1VMfdgsXFeYTi/0e8PaCJSYM79r6rq3wBQ1Vuj9/8S8Lcv++z7ynCA5ZWWS9JyCV8VqumU6Y3rxGsLmisT1i8UNAsZonmXHdX4Xk+FoBPZ4epZPyYWFqz02/dObhzBA8P2ffWoq4XY2QrzbD3hfFtxx82ZlS1Xp2temZ9wo1qy8CY0NnEtrXqcJGauYSL9/0qdCmauZqsFp92MynVs2pJlO9lz7KttSYzOyAit5VKkdbitENaC38gwwQ0O3ufIPEMvBmdlxcd8wY1Bo4MMQw9tSTKc3hLOQjcxUbV2bhNgcSZM7ikHb7dM7mzxt09Id+6RtjXxPRp+vJ89rXq7Z8fBkxkz3hss0+4XBeuwfuUhR/ZQ4NpH7MC6K1lMarpkT39wiUVVc1xseP7qGc3plOOfX+PPa+topEqalWxfmHHyDQXb60Y1S4V970WMb7fPDM9Ei1xnOYHaR7IDt72HLjAmSul377WZVx6To4tuUOLrWk/KmDgZqjCxL0hRdknLAcbp99DvchehX+rwxq/17JMx3t2P88bTbDxNUVLMGuaLrR1f54wTj62aNI2KtfJKwYXMbimV6bRhI0rbmghajI779YxPH9ynkMRZZ91/1ttyqFnoJ54hWu9xUr2A3I3H2V188/FMTGf4LwM/o6r//ej1F0dKqf8O8FMf+MsvsxzVpze+Bm9AWZRMr12B2ZR4NOf0X1iQCpO53VwVYuWGrk49b7zHnyWzUFy7mxCGwCRqLmCy7kXdeHEhWRN944i+oFvYdTqY1Dw3XXFQ1BwVG14sT7jqz5m7hkI6Wg2DkwdoNbBNBQu/YZUqTvNO6rxcWHcld9Zzk5zYBrQ2IUBR00oStXMIK0exMqwd7Bz2zjMzY0QV2l0OAXbMIc268DtxQhmYNYjJNJhEgnD1Z1qmb5/jlhv09Ix4/wGq+uGlCR66xk/qiz6YPTMOflzMpN0jhjXtHNQexDx6bTyOKTreOTukS7b0A4jOcfd8zo91n+JsWzEB/FmNW66Gz/n1ltl5jd8esXylZP2Co75q+jNDwAt7Sz8NtmqQtSduPbdbx/zQVBiLDNHAw0VMdlojZyyWN/DOEZIj+khXOJom0PYrBmfbOW+rntjz04cBufDvGOLZc4Z547GUQSs7ZlCh0Ln9lagKbIW2mdCWJb6KhKLDF+3gSy92YEopV9sGmJQtB5VJG58A7XlJsyl45+yQw3LLItR0yQ+5iv57NO1feOmdfEaa+n4sT8h+I/B7gS+JyE/k1/4E8HtE5Nts0HgD+M/e74uWPDj/+/oDX/5Ae2/IKd5sP/aBPn2ZXQfufuRv+ej2LBzHx3UMn37kOx8DBv+49sw4+IFFo7qL1i+ae7wnuC86Ssmx3pZDWzkwyp6q8CBadFyGHAXEZB4iWEcjWW2o3hZcPafYTFi2js0NS0LBDsK2/WFwQUjAjn+ektBFjxN2BUEXiprGDn+QBM7cfetxmouNyuwocys+cblYK5d360XnfekYj+Ct8VpfybK0QnHi8bU1VemuJYt4WkdPr+yjICJocsTGEX2AkAYhMh8iVWWsmP48C2+NPA7KmlloSOqYFS1v6xHtWclSZ3zFX+PFxZKDot7T8FG14x3j8DZg7LDNlOVx3W7bD6tWoKr/6BGj+WF04L/8qMKVr5eJyBee9jE8K8fx1I7hl7yDB/qmGpdZj8H34kP5mb/cesebIHYewqjHq5qTj12wyFlGbxh52yhQ7Ro5WVJuG/x2gevmiHq2VyFVGdPO+xoolt66OSFKNbXKVcjsCBXTMsu76+Hi4f9LTkZEB8xeVQghDQrGbpxozQnkgaj/iPHYs977JQxW6gTtTHK2OoWmEbojU4DqteOlX+6Ozl06gcbe7PnWaZrhlKplWraUPjIrGq5PVsxDncfEEVzkTjmn1Qrdeh7cP6DpAtfmawofLcWR++v2zn2QCRZyspadxnzKS3HY4fZPJ7f1iX1ie3YxOf71smfLwWcTJxbFX6xulZFi4wiLF0bRWh/JDmJYO3bH2FStWrNYq+nQ5IbbOEEPJkgd0NUGluf4tuNg21GcH/DgcyWb5yFOMyzTR+89gjDrCGXH8cEm89pdFhDbP8ceox7L+/bVrTH10IRF8UNE79OgM993ULLxUqvCk/6NPJMM2FU/YP3J78aOzuG2Dr81JpAGcI1SrARprG/rsFLxth/6Vnx9orS10nBpBBEhEWhqTzpw+JzzOChqXp3ex4lyHivuNzNOm6lBZ1U0KGgVOF8fcF7OKOaNnV8vn+DS4NAvOu1Bk37QzcnbZVmGT+wT+6Vqz4yD1xgHzRkJwVrjjZz7oPbo5dKobHD4Y/y3j/Jgh+fCIF2gq8Ds3RbJDb8FoG5wDzSvJrLWxXaLu9sxXW3wm6ucbCesXhKaI0gh7VW2hjIynzbMyyZTIh1tcoPTvhipX2zw0QuSjbe1CkAdRKO0x7Tze9I7srxSGJz8GIvvZ8ExHh+tUrC6Zwmt808L7SKxSS4nqdSgjzLBJBreHiIhGPTSq0Q2tRVlpa3PXaNsLEzH3rNuC+ah4ayb8rX1FW6eHbJcTeiaACkXTFWRuAnWrOE0EJeBNEnDOfRJtL6A5yItdEi+9vfA+Lyfrn3P+2/ysduzcAzwbBzH0zmGTyAa0GiVkXiPKyE17fCeOAHv6SpBL/PxFxKL9hnrY6rJ0ugycp6pc0xuBap3HkDTmhyxmHPUbW3OvQ+7fX5vW1O8c8KxHuHaKeevCNvnrPLP4AFzoCnryIztYhu98d+P6Gq4l3ztNd5TdEMUP6xK+kg+R+pDNAu7FUa/bYY7VEEaT3XPsfiqaYOvXrIy8FQocQIUCaIQDlpms5ppabh6P9z913fTmk1TsNmUtKsCqb1RIreejZSoGu//a6fHnJzM0XWwiUBsVaB4ZBJxZbQg3DmkdUjtLG3QGa9bvVqV5ZjdM8olaGuQ0nCuTyuzNbILSoC/ZI8Bno3jeCrH8Is1ySoi/yXwn2DP+5ewhh8vAn8NuIY14v69qto81hdqsv6TYE5+kitZ+y7ozhHLHCr3EdyosnOcaOzbyPWOPWXcdmgtt/VMboOstwPtQ1PaOfX+dy68oiwgJmRTU9xbsSgcSIV6x/Y5c1QSLPHZO3MhR9hcYPc84vT7JOtYhMzOxYpA+gYe40kh9VWdumOv2IfyT6/bMoIqNAk0jvKeY35TOXxjS3sQcDHQzRNxKsPn/CRyMN8yq5pBK6ffb3+cwSUWk9q02/2E+qyyCaYT4rJgdV6wyisGyU0/NOwicYkCy2BCWiFZ7YEkY/Ok0TZqCphaReTCbK4Zg3/cRPwn9ol9Xe0pOfjLBdgfw0TkZeA/B36dqn4LpsjyncB/iwkyfQ54APyBx/1OTWraMzEZJp71aMTlDF8IWav6kghtHFL2L+l+FNw799RZd5nJSTL2jJ2QRfGaO2KEQN9FSrvOGm9Ho6zItqF4sGV2q2X2rlKciVWg+kQIaZAFDmPxrdFPj7Vflli9yKoB68dat4GuCUOhj/bCXxmKGfqT9j/DTnWY2AZL4LaOyV3h4O2O4u45YRttFRIUnUZ0EnEhUVYtk7Ldc+6XHWdwiUnomE9rXBWH8ZbO4c49s9cLrvyE5+hnPLO3AsWJR5o8gfUSs61NPD3uPm5ivufPx7zY8f9+fyLb0+d5CiYiv01EviwivyAi3/3+n3hi+31DRL4kIj/Ra4uLyFUR+Xsi8vP595UnvM/vFZHbIvJTo9cu3aeY/U95XH5SRH7Nx3wcf0pEbubx+Imsxti/98fzcXxZRP6NJ3UcD5k+5s8Ttg/t4LMFYCoiAZhh7N3fjDXgBvg+4N/+oF86OPmUjBvvveHzRSAWl0Ro42d6vHwf4bRj+EI7h2uEYpX2G2v3RGpxSPDgsjh1TObkXSZbp4SstlT3thy83TG9Lbi1DWXhreejH1EEL9o4T/gQJn/J9k3naeoCbRza89K1T7Syq+7sZVGHyU5HWPTIGXaO8oFj/m5i8vYSOVni6rhT3wumoyNeKUYa9r1dbD3YmxNlWnSEMg45CTqhWDoOX09c/8k11356y5Wfiyxeh+k7nrD0phfjsUkyNwcfoHanez92AJeM6cWJ4GN6YB7Xct/MvwD8duCbMQ79N38dD+FfVdVvG1ECvxv4B6r6eeAf5P+fpP1vwG+78Nqj9vnbgc/nnz8I/MWP+TjAgs5vyz8/BJCvx3cCvyJ/5n/O1+2JmsBQift+P0/aPrSDV9WbwJ8DvoY59lMMkjlR1Z7E/BYfUoVPYxycal8EpcGTKvajVNgPj/dgmnESLkMlfbVpJ4SNVa72XZAMmO6de/4dgvVr9fYaAEmRLuKWWybvrFncjJSnDo1u0Jx5HBs7+R7KGb8OFr03bSCug0W3nR3/0Ow69ZIGsi9XMHZ4mfHSY/Zu45i9q8zf3CC37qObDX7dDqp5/eedWHeli9r1vV3m6INLlGWXawJM/a86gcOvrCjevEf1+l0Of/o+1//Ziutf6li8BsXJWDHzgpN37PRzIE8aj7htnz7kPrZ/EfgFVX0tQ5R/DfiOp3g834EFXPAhA6/3MlX9h8D9x9zndwB/Rc3+CXAsIi9+jMfxKPsO4K+paq2qrwO/gF23J2s55nycnydtHwWiuYIN0GeBl4A5l8+cj/r8HxSRL4jIF1rqS7fRttvBImWJlgWxujAQvTO/4ODHeuQ6wrKHzwj4OppzT4p2edWAYfHamJiQlAUynSKTyiaDxtIJ6q3BhlvXTG7XTO6Cbv0giQC8JwxzMcC8CN0ULtFGz3I9YftggqxzpKsyOHPtLKIf5AVg5+izHo+EhC+t+43k/6UTyjMlnG1hs0WbFne2tq46QK+4iezkli92n4JHQ03zqsFXFsUXS2F2K+HqDp1N0GkFXSTcOuXgS+/w4t99m0/9/ZqD1z2uFmPCjFU5RUdOfrQiuXiv9Of97NjLwJuj/z90sPMhTIH/R0S+mIXPAJ4fySy8izUz+bjtUft8GmPzXRkO+t4RPPX1O46nBNF8lCTrvwa8rqp3AETkb2Al3sciEnIU/wpw87IPP44gE5iTl8KSnGlWEEvYq3JS9pflMJoS7W8FUDHVxgxlSAS/7qDrLKHqZAfRdHliSNHe9x6ZVOhihizX5vi3OW8cPGFZM707YXXi2RxWbMuWSej2WDtjJ9g7y3EE3GPcKdMfl9uK5WpCPCtxWzcU8Qw3QsrJ1Zx83Jv1+gmslxHwCRGjirZNoLovzG41yNmK1DRoTOhyRbHCuO/ehNP6hPHFSP1RcgvjpKsTRdeO6j5M7nVo4SEqsq2RmKz2oG3RuqG4e59PfeWI+//SK5x8o6e+Hq0Bc8xt3/r9OcmrLPY1aDJVdC+Z/GxF819v+02qelNEbgB/T0R+dvymqqo8Cj/8mOxp7HNkj6UC+rHaUzrzj4LBfw34DSIyy+JMvwX458APA78rb/P7gb/1Qb6014J3ZWGOXZNx5EOgW5SkSh8OfXvrnf0l6I0qJnblEzKJJlxUt+YleuZFD2r3lElnTZBJaklWEXQxI3tLc/51gztbM7+5Zfa2EB9ULDfVpZH82DH2EXtvUa3/6qouufNgwcndA+JpaQ73YmDaFwB1vXg9+7BVSPhJx2y+ZTYx3Rfnksk03K64+rORyRv30NMzW7VogrYhbNSw8FGicnzsl/2MrZ8I6i4QO0dxLpRLJaytY71bbZC2s9zKAIvljlb3H3DlJ+5z+BU1yMbn7j49+0kYViX4XNk6HpgLwmh7F//p2E3gU6P/HxnsPGnL8Cmqehv4mxjscKuHQfLv21+HQ3nUPr+uY6Oqt1Q1qmoC/hI7GObrdhy/6CAaVf0RLJn6YxhF0mER+X8D/FER+QWMKvmXH+8Lk3W7KUvD3HPB01DsVAS6mR+63DzkyMeDM4rqdfwaGbOuPWEtSBvJyl2jcLD/QpcLr8Qmma6DtqM7noF3BuN0HdpF48ffX7N4MzJ9y7O+N2O5rYYq1rFddI4xN8JY1yWnqymnpzPa08ogmWaQwbPR1dwA5L1uhKCEKjKf1cyzqFcZ4sDoKZbC9N0tenpG2mzpq4U1JlMijEBn8E9KVqj1XoHXxRVKmxybpiDVnuJcqM4i4XSDrLYWtSeDxAgeqtJWR2KTuqw2LN6sOfia5QkMqmF/4pYLF/tRgjNPH635UeDzIvJZESmxZN4Pftw7FZF5bjOIiMyB34qpXv4gFnDBhwi8PqQ9ap8/CPy+zKb5DcDpCMp54nYB3x+rgP4g8J0iUonIZ7Gk7z/9WA7iFyFEg6r+SeBPXnj5NT5MokKcRezeMxJc2e2rCLQzRwrv0VYPdoM09gGZHglYgrV2lEuQbbM/piMnLzmxKslWED0so17QIiDbXFOZYQPZ1ExvbTlYzEhVwVmYw1U4mNQPRevAICrWRUfTBZrG020LdO1xtRuc2tCg4VHskLEjKxJ+EplmLfrCx6H/ZXQm6VsshXC2tWKuvg2iGmNJHRRnDtSRKui80sw8k/Deoql9pyrAzqUOyCpQnijlSYesa1stqVrC2jl0UjJQU+vGktkilHdWLKae+mrB+iVM6njApsbXlIfFx/rxefrOHVXtROS7gL+LrYu+V1V/+uuw6+eBv2mLagLwf6jq3xGRHwW+X0T+APBV4Hc/yZ2KyP8JfDtwXUTewvzCn33EPn8I+B1YUnON1c98nMfx7ZepgKrqT4vI92PIQwf8IVX98KLvjzL9eBgyj2PPTCVrX6kKObnqBBfC7rEuC7qpoCEZX/qyB3lwgqM3NBfLjBJw0grFmcLYyY01b3oe/Ig2SezQpiGcN0P0Kd7von8R/Kphdqsg+QJSyZlCuipMy3Y43L6Yqe4CMQmbuqSpA2kTkLXHN6aJ3TcB7xUgRR8drNoxKq4y5344qZkW1kPTOev5WneBbhuo7iuMKoQ1qY29c8SpMLkDkwdKN4EzLTmfVEwzD34so9Bbfz5gMNOmKehyhezsTqQ42SJdfma8rZa0LNCqIJUeVwZk09h7XYS6YfrOisPDQ7bXhTQZQXL9IPQDMaZvjlZo4i84/adkmY73YRQoP8o+XwO+9ZLX72Ew6se139/ziLce2qcabe0PfR2P45Eogqr+GeDPfBzHsr+jJ/M1IjIB/iFQYf77B3Kgfak9Mw7esO4RNjt+y3vStKCdy6NZFA8Hefvv58QkUfC1MDlNVrnaQzG9TIEqYKuIvcrWEBDnkJNzNHi0KmEi5phykxCcI6w6Dt6GchU4Pys5/3Rgfb2mnHSEEAcGTV0X1shjFfDnnmplbcKsy46ihZAkO3ns3IScR8iqxH2pP6JQJsqqM1gmR9xOFIdSx8Dp+YTyZsHxV2pzuGW5yyOIw81nNIe5RdmdjmLdMTmpuC1zHoTE4WxL6aPt9gL+rirU0WCm7baAs4LFV5XprS2yrncXIAHBxkvWNb52aBnQKiB1XiU4h6y2HP1MYvXCVZbfEDP9U/ev/VDslq+7jpLOl0X2n9gn9hTtCd6ONfCbVfU8dx37RyLyf2e66UP27Dh4TcO6uxcd61+XEOjmBe3BaPsxbMH+62OGzVh1EYAEfitU93fsmUFoDAwyKAp0WllCsKdqOgfBm5RwFyF2Fo32yo91BBfxXcJvO4ozT3lWMr3rqY+ndDOIJXTejjE0WZhxA9WJMnkQ8dtM0wxCO3PUx47tNWs2MibK9F18epliPPgyURYdpe8liu2MOnWcbCZ0t2bc+Fmlev2OjfGksryCKm46of6ml6mvKNUD8HUiPNiwWHe4bsattODOKwWzucE+PvddjVnvftxLs7hV8tw/U45+7gx3uja2UQ+5eW9FYnWT8fiETCdo8EjT2hirQheR9ZbnvzChvjaluZZyC7r3wV6G6F4/cfKf2LNlT+hWzKuf8/xvkX8e+e3PjoOHvcrVwcQhsyntPJgAVp9kfJRzf9SpZg8pyeEaKJZGc5TQdx7Oy/0QDIJJaScj3FsPyaSdxIG0mF5NzEVTuUjKOcGvasqTkjgNqMtyvM5YMa61nqe+tgnBJgiyngpUhWd6J9DcDqyf82xuCM3hTsOlr3xLONRFnDOJhJ7WqBnjX24rlvfmTG87JvdbO6cikJ47Rl+4hhaO7XHFvW8qifOO+jiQSofUHVJ3HLwWce0B979pyuZGxflBgjLhiqy0GR3aONzKM7vjOP75yOHPL3GnK9jWaLKmJcO4xtEqLTt0aTs7Lkk7KmRKFO+ccPDVKacTR+fVxM/2WFMXLnhuXK4h+3dhaLDyiX1iT80+WAL1ei8vke17Lgqk5WrbLwKfA/5CJrxcas+Mg9ekJuYFUBTD61IWMJvSzRyp0KE/9KXB3JCNlke/HgVfg2wyDp01b4ZJJSVo8nGI7E026p1t3+vX9LRKHDoJEOOO462KnCckzlER3LZD+pUAmaIJ0HbD5DBm82hV4DYBV5e4rkRdIAWhOVa0MrjDNQbXqJre/UXp4VVTcHJrQfVOweyWEjYRnU1Ix3NOPz+nWQipEFKAzQ0Fr9TXEpvrnupOhX+wAidMbm04KmeUJ47mKNDNbPXgWnAdhDWUp8r0Xsfk7hZ3th6cO2BjmCe+fmzox7e1WgTNiVYJfnc96obj12o2NybEyhkcdRFf74vYsrjZ7vqPIJxP7BN7ivYBF5N336/jVE4Ef5uIHGMJ9W9R1Uv7Az8zDh4ATeasklEmewkBrQq6yhzR+86Ee06+xzT6qN+a7PoagwlgJ1HcSxJ0JpEgIjCx5s9D1Km6Sxj27ZmGaFSMYthH99Y6Cg0O10bcujYYok/cbrYD5r/XotCJYf1Ni3iHrAJuOwOdmZMvhLZUdJIGXR6ZmE67HapNINsu8ODBAdOvFczfVma3oyWInWP98pSTzznahaX3JcrQbzYuIqsXCyb3Z8zONvb9SQmbRNgadbFYQnmuFCvFtTZx+E2HX7c2cbXd4Nx7GWYNHkm6O/++/qCvINYEKebmLDvNn/LmKbN3K7oDoanEJBB6KEZG0shZfgL3xFbDn9gn9sTs41hIqtYeAgIAACAASURBVOqJiPwwpiDw7Dt4TTkKVcPGaawoRgtPrLLuenbg4zqXnlwhPd6uD+fkepMouNYctYI5FU07Pfjgh+hdq3IHH/TOvYsWjcqww+zck70/LuQpCrQMuLONyRJ3u0SiNo1tW4QhwdvDQVqUO4aPOPz5mtn2EF8f4rqSpTiaGwl30OKcUlatFTOJ0kRPTI6z1YTizYqj1xLTux1h2SLbDoJn+Upg+1JnBV+tg16DRoAisXk+cbYMVPcPkDrSHFdsr3rqI9OLn92JLF5f4c42u1WMc+DdQIeUsSPvnbzDIvQME6l3xoFPCdWQNx9dNeeQs3MWbx6zvV7RHmdq6phSlJ17D93pxYqRTyL4T+xZsCfHonkOaLNznwL/Oqbge6k9Uw4eYKhcTcH45zGS5hXdTEhBcRco2Q8VL15Iqu4lWLNPdh3mZGDHntnWECIyn6IHs5xYlREMw84Jh1GOwGfIZb2193JULmIMG39via7WO6bfMJGEXPBj7BxUjR5aFkhZmExCMqepqnD3hOrOfW68fZX5rSs8+MaC5eeFgxeWHE23gNEUmy5wvqlIb8y5/lPKwZtbXBNJVWDz6oLtVc/yM2mAOqSwXACdDBWs8UrH2TcUqD9gejdx9hlHKqA4h4O3I4c/e4Kcr/ej8aRZfiBfoB7ayuMrXUSLkBdSik5K0rxC6hZ3PydcndutpGKyya9zTN5eMn25ZPuco5v2EytobhrS69RocvuiZB9XeeAn9ol9UHtyt+GLwPdlHN4B36+qf/tRGz9zDl6TMhbslLKkPSpp5zAA8A99iB0urxh0koSLVaSQk5M9lAJGz+xrG5JVpTKtSAczq3TtfI7cO4sUndtF6c4kjOn163t8HXaR62abvzvDTShChxSFRfJdRMnOzZs0gtY1OI+UWc0SLOKvO3jnNvP7p8xev8rpzWPe+a1zrn5qbb1cO6NDup+b8+r/11C9fQZJSYdTNs+X3PsWz/bVxnTri0RRdrjczzRGR+w8sXMoEI87Tn6FcCIwub6ifmfGwZuOg6+ukJOlnTfYxLTJk5s4y5lcUhms3hkDaZtpkzEhmxa32qB1Y5/rmVN1Y8lZV6FVgaw2HLzTsr1esjwWKPLT0rlBusCurez3YO1bf30SxX9iT9OeYJyhqj8J/OrH3f6Zc/CD5WYbMqtoFs4YNLAHqw+2x6zAqh97dcEhirOHXSKEOkfiY2pmVVhyN3hk2+BE0DIM0APJmZMfJYBtfxmicSNqZ0/zHEX/vQTxMLF0HUwqROudXHEIQx9YyZWeZEkEYoSqQsoS2hY5X1OcH+LKHS3ybDMhvj3jxR+NVP/sDYucD+Z0Lx2yuuHZ3ugoZg3TiUE63qWBbWNt9TybbUFsvUmuiyBBiZ2z6lqBVAXcbALBE4+mkBS/3CJnK3S7RTvZX6X017IISNsNOQxZbZCusxVMVf7/7Z1rjGXZVd9/a59z7qte/Zyenoc9Y3uwwHZsE8e2wAgHApjEiSNFQkZKYgiSQYEPkSIlJvkAX5CQIkUiCkSMhGWQCAYlQTiSg3GQIoSCg81DfuHH2B57pufZU93VVfd1Hnvlw9r7nHOrqrtrxt1dNT37L5Wq6ta5t865j7XX/q//+i9bMNQK7a2523gYFk9PcXXJ6HLB3qtDsbX3erffGwkbhiiVjAlByuITjhnH9BY8cQFe+iPX1JtEcuLwA2W/OuZQekbo/EugI+fjZ74WkygG73eyMJIv3yeBLHJbN/w+hUukD7x1u+pkZEoY6Lj0iDpYHcdsNgb3EOi1CM1TQftN4OVlMrbHqetWq47XjvYYDmhOb3DtoZzxZErVZDQqTHdHjC87hpfnSFGgm+vMH9zi6usGTF+lyKRhNKqsMzV0gXrMyVJVGBZhnqx3aBMC5TSjuVqw9rRQTBt7fgcFmgl7r5qw8xqH1BusPX2Wra/NyJ+5ajuc0DdAfF7A6hl9uisidrhORmjg5KmbYA1hf3OziuE1TzZz1JOmkya4m3xyUvaecALwircqaCFhwIc3dYmOh9Rj8FkI6PsC9oG7x6RZu7jcbtP701NC0a9FlD42XUFVvbegFP82KIy6Cd4zmmf49RFub2EBWHtf3ncBHVpunnAqkd5pTplSR5YVLnO2YOzOjJIpq06HX+QtzeNPrTN9eJ29B5WRKLOyYFEWcGXA4CrghOrhC+w9OGbvAcf8glKfrRiMKya94A6d1UDmPDlQZ469eU62XVBMBbc07n182dsELDDrX0wqadp8u55iPmZt2ZDtTOn3BLTPb/D2RwQKWVkANHP4UdgxLSvcfAmLpTVCZQ5pGvJZQz7LqDXUC1YmOAW5ZOLdE04gjusteaICvHmiiE1vChSKjoamoMl9F9BvlpWtzCUNQdVpq7BZQZy52vtdgt5dqtVAraMBzeaQZpSTzUcgUE9yRtN+O36UTbIS4LWvkgEkz1HnaEY5mgsut0x1ft86k69WMJvZeQXJKGFerA4L5vevc/W1GfXZ0jpIq4zl7pDhVUe2VOb3DJmfc+y8FuqzFW5cMygaJqOSQdbQHOJwGb1mMqdk2wVbj5laRgVG25VJHL1JJjWz+4+2a9afMI+gfK40hdBsDnCLMlAtDg27DlmW9lwP846/rxtbQEW6cYi1RxYVOptDXdvuKhuAV/JFQzGFRSyg9iZWpVQ94cQiJJzHgRMV4AGzDBbBB48Uvz7AD+k8WUSv27G+f+4F0I6pU4cVSTEVB9AG40ihaOTTxfTZK7x607T88exCgS8GiFeKmTL+2rIdYNGqQGLQisXY1nArs4Iq4GYLBnuz9m86HpLPx0FJMwi7gmh85pDhgNnrznHlOwpm93mk8GbpW2e4aznSwPKMsDifMbuvga2KrPA45ymKhtHAmqv6mx+hc7rMnKfxwugF4dRXl4y+cRkdDWzISebQ0dACtHNo7hg8P+feS7udeVmeWUCvahgO8OMCqb25STa+NWnzg9x6HRrfZfRNY9YGswU6nVp/Qp6jZWnnOxpAo+RTtdcxly647x/fuPIGSNl8wglACvCEbM2yd61qJMuo1gqagpZisQO57hO2EvyDisKGQ5jawhdKMwhyPLXA1FInfZlecLQ0dDRDcWmbU9fWqM5PqNZzo4LGQ2Q6b/3NNVAPogrLpivoRjhn04wWS6jCZKgsQ+qa4bfCbmFtbMXHZYnWDTIomL3hIs+9tWBx3qMDRRcZg40l070Bk+cd+cwWMl/AcDuj9II/XZKPG4qsG57dd4DcP7jj6pU1zj2tDJ/ZNYoE0KpCGtOsx12IDAqqs2swF7LpHN2bWsadZchkwvT151iczsiWymh7TDavcfMacodUDVJ2BVdUkUbRqKTJsjATICyG3noMXNlQzKyOokV4ffs0zcqbIAX2hJOBdpN5DLhpgBeRDwPvBZ5T1TeG284Avws8hPkr/5iqXgmTnX4F83qeAT+hqn95lBORLDOFSpiHKs5Mx+pJhhbhs3yDnfgBO92+uiIG+czjh1CuBVVHoGbaWax4uy3rqTogBPqgimkUtztluDdjUOTmKglokXdNS33qIc87O4Q+hgMLcHkWFp8MnYwsg10uw7CTsONYn7D75gs8/+aceqwMtx2uhvm9QlVljJ4ccPaLNcMrVetlI7WyPD3ghTeOmD6cMbhg19ofv7c/uF/a3mLtcyO2vj63rB1W9emEQilA48lmBfXmkHrjHFKdwdUeFWF+ccTOwxnNwHoO5mcda8/lTJ6ckU3LFY8fzZzRP2WFbozx61vIch23s2eLoPemm99dkgPDi+OOgjOpj70xfK/+kpBwwiD+eN6YR8ngPwL8Z+C3erd9CPhjVf1lEflQ+P3fAj+KTUV5BHgHNgvxHUc5Ecmie6Qa/44F/Xrs8Af8Rw7evx26HWpsGo5TFSTcKE7wA6UZWes8TTS3ilYEYjp8ytVGnf5Iv0gp1A2yLK3BKR7X+C5gxQAfLQ5aFU4I6GXVBvZWcSLSDvvWxsPWOvWFLXZeM2b31XYOm1+H4Y6nmgjzC1BOB5x/XJk80XWWSuOhrCieKij2zvLC3oidapPi1Q2jwhQ/eTAliza/uztjxl8acf6zJcUzO0Y3xQYmWFUTxeauZYNbNNQbBfWmGaoNdmrymWfzcchKRWrFVZ583uCWXZevcfrevHTyzCY61b2axfoYmqEVsxdL08qPalMjxJSo7SYOb4CWd9qnhU9IOE6cZA5eVf9ERB7ad/P7sKkpAL8J/B8swL8P+K1gafkpETklIhePNI4rWMlqz71Rsox6KN3gC6TNyA/gwDadlQ+8AOI81chTrYVsO+iw96d9GgqrbeYaFTKRT4/BWrXzrYnZdh2lll3W3na2toXE2gK4k6CoCfctw2CQjXWa0xssLk7Yu5izOC9kJUyeUdafLBFVdh8Yog6yKzlrT1dmhzC3Yq9WlS0g0ynDsuJ0fpF6PGTn9Bp+fU6eeXbrjMWyoJwOwmM4zn6xYvTNXpdqnttrEp+L9rm26xXvyWYlmrvgkqnkswq3XTPy3hrFoGf9rAfVRuG1j9y9i5bJQTuvw8KoLm+7HV/IdSiZQ5CCe8IJwYmlaK6DC72g/Qw2JgzgfuCJ3nFPhtsOBHgR+SDwQYAREwsCTdMVPoOaphmw2n4e77+PgrEbaQO7qPmWxOEU4hTnlHpcU09y/LDA9XlxERDXNljhw1CMPvrabmcBsCvY+gNDQgTamaOqavpuAuevHsjahYF4v1MblPdtMb1oSphqDbI5nPp6zeRbu8hsSX1uw4zXRBk9nzHcDkM8MhsQLs4CLg3o7h6jp3ZZP18wfXDANY81ge0W5Ncc69vC5Fll8/E5xdNXkencAvpoiI6GtkvxflXTHs9ZFZmX5CJki+AtP10iuzN0sbBrLwrTzbv+Ith011vk6HgA3nculKpGUTXehquMBkbTDAqaQlDnu4w9Turqf4D62by2r0RCwvHhZRbgW6iqyo2mMl//fo9iQ7rZdGeCwNu3xUgNPi2+iDWzm3xI92f2cUyfQDGoyTKP9w5XeJqR4se51edym57RFlqjKVZWdAVXO2FrxCkyZF6ahLIX6FZmyWKUUxwkoiFzVYxPFpFWKtnOKC1y/OaY2f0Trr0qpxnCYFc586WGtcd3cc9daXcA0qyZdXolbD7ukbJu6R8dOKSyIi5FWKAaz2DXM346p7kyZnQZNi41FNdK8mltQ7F39roCalGgoyH1mTWQNfLLex09Eymn6D0TpjNJ/7nIs66GECZeSaOd2qbVwuc0W2vUW0OK7ZkVlBcL2kHo8/D4k3HbN6AuLN5eV2VT+xb2eCoJCScBL7cM/tlIvYSJ5c+F2y8BD/aOeyDcdjQ4Z3LAoPaQPEPzQ7SPfRyQxdFldsFhsGkcVZlTVxlaO9ssjHJysIYir+aEmOeWlR9SFNUip7x3A3XC8MnQ/FRVFsQjlRO7UrOYlZv23f7eSS61DDrxQUFzbpNqc0C1nlFuOJabjtEVz8bjCwZPXUV3p/ZYcQhJCK6uVoaXHRtf3zXuOi5QddO1+weaRb1nuL3k9FfsvMbPLs06uPZmSzxf2mNHWiYM4ZDGs/PIGqfKBre76II8gJqyRTPXdaiGhiUtejub1ikz7JaiBj4saq6sGT6+i+7s4suyo7YGBTIemX6+VwdZ6QjsBXc5LJMHk1SmQJ9w3HiZBfiPAR/ApqZ/APiD3u0/JyIfxYqrO0fi3/sX7wStOm7WF6waSMXP8Q3iviiWwdeCrzKazOagamODMXyGcblegSYMmsjbLFyryjh4keCy6KwwuVvi5pUZZgUFjuQ5BIMtHQ87BU1Qf7SuiFFiuLQgJutrVPeeYnFhiM+EfO7ZeGLJ1rwmuxr04MvSdjRZBuNRmKHa2LDqBjYf92Tbe13QjdRHoLoUhXph5/6CY31Z44cm7ZSqZ20cFUB9GwHvyXbmTJ4fmqwxBlqxqVRAsDsedS6SkauP19vP6nuLQysjXZboC1Oa3b12fKJdRw2zudFE65N2J+AnA6q1QNHE94LTjqbpUXTtm6Q3bD0h4VignFyrAhH5Haygek5EngR+AQvsvyciPwV8E/ixcPjHMYnkY5hM8iePfCYxg+wjc52dQFBNHGBqetRM6wkfP+iNQGVTi7Lc451anTZXmpFbndakijRh+EQcvwcW/EdD/Naaabf3ZibpixxyuB9FYTJH7y3YqaKba1BWpmppGrtf45GiwG+t27Sk50vcssbNzJogSgZpmqAscvbYwwJZKOQmLXS1MtyuWzmmDdTwEJQokjnrhA2qHVmWOO9xg8KCbZwkhT1eeBI6S4Rg0zC6tGsa/8WyVfxIeL20rFZrEPF7HHgSTdecs3d4VN9E/XtVd3YMEs4hum6WJbo3RbY2TYIqQjPMusHrKpD1ayj995K0nkTSpACfcLyIoq/jwFFUND9+nT/94CHHKvCzL+lMoqujdBzqSkB/CYoIUcvatTEHwrarNVPqses1MtFlvLHwGqgCGY/wG2P8KCe7vGvNSeF41FshMwsyv9CQ0wa7yFXHwd2q5vV+aoNmY4ibVd2kpzjKri3iuvY5aINuLKQCrlKKnSDRbHzHgcdidXwO4uPEgB6lnDG4x+Kn92ZDAN20Ku9tQVuWphaCliaKU5i0rKyOEREz+NaO2TL+9hi/L5Vx0gX3/vxb79DFEpnUds15hhaOZtx/gbsfVYJSKQ4AgTAIhISE48cxFYROTier99bodJMnoj+9aT+vtTLIJ8RfqQWtBR2AZHZnzUKzk8tae97VTlMB75DxGH9qnWZtgFtUpjDpZ5vxfxWFtdJHygNo7Q6WZRcEJ2P85oTyzBhXNrjdWcd/9wN6f5BGnllALqtghTAytUqpuFlpfPciNA9F2+JWsx6zYt8NvwZwPcvj/q6pN0y8pVXqaBAWCtGt1LHH80f5Z4BmbmUxEN/bJcTFRTWokqKvvF8ZXSguqJLKMtRiMpphRjPsP/G994AL7wuPDd4OvHzK4BNOAk5sBn+nEBUnKwFHtX1iVgJ6/zMrdE1O9I6N372gjR2UFzW+cTSZUq2LZdOzOrS7SxewAIqC5p4t6s0hUnmbzFSF5iRicAvWCmtj497ny45rjha4dRgqfWqD5atOM71YkM+VjS9P7fiYtfYpjrDI6WhgYwObBhZlOz9WGjPekmWFn4zCHNSeD06W2bk6se9AtEYGuq666JqZZ90uIe6kfAjkzlnBs+9smTtbJNy+XRAYzTRbdMXV8JpKnqHDgUkiA/8vZRUopNXgvvKeCOekg5x64qjHijSy6gkf98Bx9a8l2FOkAmvCCcAhyeidwskI8IJ5mGQ2TUicrDQ87fcWkX2BfoV336eeI2R0vhZkAHnRUI8aylPBy302Dw8SFC5RFZM53KIm9+CWlXmtRDS+C5giFqj6ARCA3OgaEdhcZ+fN59l5yGiK+/50D3n2BRvL19ffQ8fXO4cOB9Rnx0jlKa5cg+USNtZwVdMeJ96bTjwMDG8zbjC6S1y3O4lBt7eIAEilqHO2WJRVRzPFRWpQ2MJX1d1c2eidr7JaQC1L85CJlEz8P01jXj2q+HEBk4FZJC+WQfLog4IqvKjRj8abKsiPC5abjmYcHs8prmiMftNwPvH90PTfHAffbgkJdxq3qsgqIg9irgIXsIj3qKr+yvWOPxEBXuKot+EQqetOuuiVOFszTl+7rheNXy2ytrbxaoU2bRzeC4NBQzGuKE8VNGfXya7sBDrC+HQLtliGOV/iQmYtmxtWJF0sghVwExqperNds6y1HpYwhUmXJZLnjF6oyOeeYqcif3bHqJey6XYEhIDYNGYulpvbZOG9dYpGv5pgaRB16LKz10kPYwCPNrvt0xUKynG8YKSBgqyxky6C3wj+64rp6xu1JiPvkVLROL4QOoVPnCvrXBg1WBxYRHAOvzHCzUpkujBv/XGBP3earPHo7m7rCSSZMzfN4dCy+8WSZpSz3BJ04O3pyhWXKU1rHxGegrjYx59f9LsxIeHW4xaqaGrgX6vqX4rIBvAXIvJJVf3iYQefiAAPGP0QBkK44RAfFSoeUOnUGhzCwfc+0PF7P8jHwlvTOFSFwaCm3KhZnh0y+ZojjC86yEeXVauJ12FhTEC0GYhZcfR5j66HoWtVGx+yzwaKnGJ7TvG84qbzjo5xQS4YaJVWolnXKAPrIm2arvgZm4fi/Yu8y9qhO4/R0Lj53WC72xYpc4R6VbEUeXbvkdpZ8BzZQiFlbQtdNANblp37JdhiGCWmga6RPENnC/s5jjcMfLvUpjCSqoZlRTYNC82ggK3NsLibRFVVbQBKAzIasjxTUG12/9oNGvKiQVVoenb+/Vr8fuouIeFYoNyyImuQnT8dft4Vkb/B3AJOdoBvC3zBokACPSD+OpnYYc9XTzLZh3hT0/ja0TSOQVFTrFVML0xYi5k0tNm0Ru2091BZEbVtNGof1AKatk08zUE+2iuMR0ahLCrTki+WXTNR3ze+abrxfOJsh1BiCwzhHDKzIJDaW7CMxcpBYfRJRCh6av9N1fhOnti3GqgbaEJQnYxCg5PRIsa7h9ek5dNdS6e01+i9+bxPRtapW9ui1ip6SjMuk3nZGqFpr6lJBgPL2HuKH1VtLZ11Y8Jy01GtRXoGsrxBVoot/SJM772QAnzCCcDtKLIGj7C3Av/vesecnAAPQUoYOOk4Cajh5vvsQz7MK/bBkaapHU0jUMB4XDK9f80aiNriKS3FsaKwEWk7RCON0Abz6FsTOzBd1v0tc8hkjM8y3HJhihoNlEfjjaeOhl5Rf58Fbj8sGu0mZRKUOtBm1lRVS9u0wTs0OrUF61AniNr3dhhJRN30CrFBLrmouqapQLsgVpQmc90OpjUBG9JsrVGdHlHsLMgG1hNA3VhvQWiAkpidlyU6X6BNYzw7dHJOl0FBd/5NQ3VmwnJL8CNvI/lypSgavO8boMVgv/o+SNbwCScCR38PnhORz/R+fzTYuqxARNaB/w78K1W9dr0HOxEBXsMW3pQ0qzSJRB6+hwOrYT9LO0RpIwpam7qiqTO81qwNS669qqI5u0E2m3dBrt9RGX+PiNbCAGRBm17bbT4Q/nnMjq2hScfDtvlJ+0E4mo1Fzbn6VU044TG9DR5RJ/jJwDjsZdVl7HFubLx8722nEOmcKGGMC5dmXadp26SkRq30aaq+qoes8+TJM3tqG9922vozG8wfWKPccJx+ds/OowiLzj67YZ3OLLhX9rwp2K4mPtd5joxG3a6k8cwuDilPgR94cOAmNYO8Yb50K+rW674nEhKOEVHkdURcVtW33fDxRAosuP+2qv6PGx17IgI86OpQjJgdNg3S7PusXoeGORRBJRcHbWst+DJjmeeMBxXnLu5w9TvOcPaFMewE//NB3kolVT1kwQmyrCwwxnPsZ+4QqCVnVEoMnnm20sAkw0EIstr5wtdNJ6VsPCJ1NyM2UlVksLOHy7I22MddAIDosvPCAQvqsTAa/1cM3E3Tes/3LXs1WCxY5p1Zl2i8tn7zVh4sfLPQwTscU54eoxkMrzTWGFU3UHv8+hAZFci8hNnCFrnZ3II7BCoqyCDLqvXOIc9sutO1XeT0Fnv3Z5SnPOQKThkM7f51naFxm9ZW4Xtvlfh7CvYJx4ko970FCEOVfgP4G1X9jzc7/mQEeMVG9LX6cdcqUlxFJ3280UPIIatkbzGQ2PRUCU1lQy4ubl3j8TcJW1/bJJ/NQ4ANAT0O5whqHg0FzzYjjVl71qlfAAvg8ba6m+OqE8tKpYqFS9+2+1PXrTEZeW76/NhpGqgbyRpksURHw6BcyaAyCkPVOoBbdUxcHPu3ZQ6IDo9hcYn8eft8mVxSJyMkz9pB2fYEdguChMzaOHbP4NJVBk9Jj9JR3N6MZrTJ8sIa9WSDYrdh9OVn7LF6OyTp00WZLYIa5JQyGnHl7fcyvV/xaw0UnmzYsDFZMFsOQn9D77XuX0pMm5InfMJJwK2jCb8X+GfA50Tkr8Nt/05VP37YwScjwINt1+vastBByFzrGldz5CdHXeDssfu09Gsb5AX1ZkBWlhlVkyEPTdl7cMzW5Qm8cNUKrF5X/VmiP4r6Vn6oDaZsbAuZWbAB0M6rZSV4Di0rzRwyN6907dsjx8cpAkURA3/T2IxagNnCMoFQwJQ87+Jb2VO30AucvSx9xVAswjkkb++Emy26+0S6Kkor4+P1+fm6Rua9xyw6i2V3bc5wUTMozGZAl6X1N/R0763zJtj/6hmf1fefZfs7Hc2pChk2ZLlnMlniRCmXub1W7fg+WQ3mYbVPM7cTTgJu1ftQVf+UF7EnPRkBPp5uEygPH6iQqsbV2lnB3ughYrLWK6weQCy2NsbFz6qC86f22H1wnfUnNsheuNoN247BImawUfa3z+ulRW4TolYCa98Tpqpt8EURBknve4zY5GW/hC5bJ1DaY7WPG+mb4CEjznXcfmyW6lsdxMamRg5m4xH9aU1R9dJ25Gbdz70OXZyz4+J9e772UtW2QHqP25221sgadyqtxNSbDDLaK8TnLXSuXn39OsuLFdlaRZ43DAY1G6MlZZPRNL1zDtn7dT9EKconHCeU1VreHYS7+SG3HyKuLS5qDPKBKnH1YRrJfbjB31c+23EbHxQ187JgUlTM7vcszo+MIwfagRPRhzxMbyI08Kzo4OP/yfqyyeBK2Q+iy7Kb0xooiBZRGjoctNr4dmhG5lqdvGQODSoUlmXg9rs1WkRaT/sDzpxtVu67JqS+E2RbVDXuvlULxfuEGkN0h9TCuHLyrLNHCLbJdk3hrVU36GKBTmfdc9Z7vS3IN10xOCxW5ZkxV18P4zNzhsOK4cAK42tFSVVnnZQVwIslAZE926+iSUg4bugRv24xTkSAJ8uQ0bAtummrBfe4JsqcezK4/vc+9j9Rh6hvpJHWgGy5KKi9Qy8s2buYoVsbdmBfRRMKf7F4qmXZFUohyCl7mXEb3H2XAQdvqVj+qQAAEiVJREFUGJktkLLGb5rH+UqmnefIcBgmTAV6pB+wo199bDoKnvU6Gqyeb3Rz7D23LWefudZTvtXQxzGDMeCH3YLkNnTDnkftfGqi33zTtDa+feVN64VfW4MWzgrP7TCVIj+0X6CteYjgxwU7rx3gXr/HmfUZo0HFZFgyKUzpNFt0oxSNUmPVNVJ6kT5l7wknAG1X/U2+bjVOBEWjgxy5eA/y/DZ+b2r67uHALHibXqEskOo2tu3g4xzpCVKMrqgddZlxbTFkfX3B9IER84dOMb68HY7rPVhdW7027CzaZTF6yKi3gBu7Or1H65DVx38rIbhWNTou0I015Npe9zgQCqYKfT+X3FQ9WjfoMpiZRYtfF5wZiwKtqrCwlNaRCm2wj9x/W/SN/HnoGI3/V4aDNlibNUHVdtayLLsCOFgxdGNi6pjwvEifzokZeSzu5lm7w9DZHK3r1m9IssxsKkTQ9THXXrfB9ruWfP+Dj/PUdAsFtoYLcue5dG2Tatl72yodD690r018z5yMFCbhFY5bpaJ5sbjp219EPiwiz4nI53u3/QcR+ZKIfFZEfl9ETvX+9vMi8piIfFlEfuQoJyFlyPqGA8vkIXQ/OpNJNqxswQ/FSuZ+SHrfu680IKXRNLvTEU6U5bmG3VcVyMZ6598eeXIJhdfQWNQqXkIRtE+TUORGwThzklTVTvYYArLbmYVuz6It2uLDcbEI6QNtEYP8cICMRqsyzDhcRMTG2w2KbjcRs+/apjHJbGG2wtBl9Pu598bbkJLwGrTumWGkIUGmibPpVTosukUC7Ji4mHiPOrEsv++VMx4haxNkYx13aovszGncmdO2GI6GXH3zWZ7+u55/8qa/olFh2eScH08ZZRXXliN2die9nZoVV6WR1utj5aV3aq6TaapTwnHiqPTMMVE0HwHes++2TwJvVNW/BXwF+HkAEfku4P3AG8J9fk1EMm4CrWv02cv46cxuKIouu2u0+/A6bq6VhIN0Tvv33p8agcpRlzmLsoBxw+weobrvdEdr9PXjYLRGa++rbVAEOqomyCDt/+lqcSUEZFmUXVdpq6s3n3mJwXwQaIiYcfc58ywL9glBFRMDcZEja2NbYPq69yDn1PFwNaj3IaaIiefe0kfhvpLnZikQArTfmqBFZjuRwWB1hu1+/j9ee/wediV9ekhV2XnbRZ7+ezX/8G1/xdDVXJqeYpRXbA7m7JRjnru2btLI+HLG3VikZxyr0siw4zuucWkJCUBodNIjfd1q3DTAq+qfANv7bvsjVY3tk5/ChmsDvA/4qKouVfUb2Oi+t9/0LFTR+bzzJ4k8r3pc5YPhGC9CHBQf98Z/My7eUS5zRKDaVKYPjLtdRESkYgIFI8MBMhm3X/0GLV0ug2GWdoEu/F29tyDafvUaiKAttsampW6+ql8phHZDvnv8d9Cla5Gb42TjO28bMEVOpIwOK8CG21v/mT6cBPO0IOkcDWhGObKs0EnI5ONjex8mTwV7gzivdaVZKjMKbji0axFh9t2v5tKPeN79xi9z7+AaX9m7h0Wdc2405fJinef31lguCmtsapVSUbEUvoWJXW3NRkLN5RDRU0LCHYU/4tctxq3g4P8F8Lvh5/uxgB/xZLjtAETkg8AHAUZMUK8IjXmDO6NEBHC9DB4g6tvtQei26jci4PvUTf+4RswrvrQA3Ew8e/dlbJ3Zgue2DyhlWjVN0alU1Akyx4J6qx3XTg3j4hhC+5uGV1Ggy2Db82m6wmx77tr70R5DMtfdXnSBXvsBH7pGJ1Wk8aiEkX/R9yYi7jKy0HkbawErz6EtGBLtihVkUeE3x8bdx91LlDnGYnG8ht7UKs1cr8Epo3rgLE/8UMY/eMtf8eBom8dm9/DU3hYbgyW1Op6+tslsNuyy9yiH7W9r9wd9BypKtnQpwCccO25Hdn4UfFsBXkT+PeZP/Nsv9r7BQOdRgE05o0AY/NCbQOQVqXzg4OkGcLMvVt8guNsYN7Ggu0LEGxffdkMKaK7MLyizhzZZ29kzrxU7MfOKD7JEjYOr7UkIHH1Pi34YRbHf3yaT3uDu+D00NKm25lttoI28PJiFQjhG1wbdrNZAG4lXc7As8m5Ih6oNBynybmRe+yTFpiPXqWv2F4UkBOXgQ59Ny5Dtj/CjHNf3xBGxLL1t/Op530Qv++A/5DcnXPr+CT/zw59gJDV/dvU1PHb1HEXW8Or1bT77wn1c2x2jsXDctxZt5ZHd7F7x4bV01j+RTx0uBfiE48Rt4tePgpcc4EXkJ4D3Aj+oXbp3CXiwd9gD4bYX8cDBjlaDJ3jZ4Erbasct+E3RTvvo3xZu720BREEqQQcdZ19veC6/qWD0zCmyZ650yhVoC6EynVvwD8VTzfMuUMfuzJD9R+uDFuE4rUK23tfUB314HJ6hZWX+61EiSXC4bB0vQ7PRoKDZtGEabmdpE5AuWN07297risEi5gcfNPtmMtZ1w7a0UpymJNIuchI17iEzl6XtWNzeHB0P0CyDYWE7hbjDCb7zmmeWwbhgdxz+T31ug2e+Z4N/+U//Jxtuzoe/9S4u762xPlry5rNP8dVr53nm6dPtAt5fNzUO1+5tRMRbsdUXiqiQ7TnGz4HUx/TpSkgA4NZ50bxYvCQRmYi8B/g3wD9S1VnvTx8D3i8iQxF5GHgE+PMX/w86OZ6bVWSltpkZ7M/eD7u/rv7Yz/BVVoq2UoVAAZApWnimr6557u2b+POnOpVLK/ljtQtVvQXdPDep32QclC6h+WdzHRkNu+w3jKCzbLsX3EOjk/T/R1wwQvZsgzs6XbxWFexO0dyxODei3hxZEK0qshd2cXvLrkAbERuuIhWWZ91gjkgXNSbJ1MXCfnddsI7HydwcIKWqcbtzszjA6JfWdrgOowVDkNdRV/ydv+4c33zvOu/5wP8lw/NrX383V2ZjXnP2Bd52/gmenJ3isa/dC6WD2oGX7jLUZK5S95RVIbP3OeCg2HasPyFsPFnj6lRlTThm9PpMbvh1i3HTDF5Efgd4N+ZT/CTwC5hqZgh8MnCpn1LVn1HVL4jI72HTRWrgZ1X129ogy7IMGTyt4Vgbr/vPx35p5E3UNuIVzeiKrQ7IvfnLOLj2OiVbbnHm8xn5t54LjUUaAm74x06ADBmN0LWxURStzjw0Mu0vuLbn4FdsEFpEZYz3yHBo94uTm4KlceudHughW8O0y5bzvFP8OAdFcfANFOephvOSA+cXeXnpFrlwu01kijsL1/L5EorEWkh3bd4jNZbVh8Ef2++4l2e+3/N9b/4C03rIr3/r+xBRvvf+bzBwNV+5dg9ffvKCLbyRa++fWiNmtLYS3MPvueJmjtFlYXzZky16XboJCceBY1Ry3TTAq+qPH3Lzb9zg+F8CfunbOakWjTXmZCHAHwntYFZuEuRpA4Q0Ak1UmtgfmjXPtddmSLPOaSD7xjPQGF2hqj1axYJkG9wj3x25+ZbH7xd3V71jOg/4HjLrJNXRAD8ucPMK2dkzF8oIkWCfAG7pccu6G/zh4iLmLPjHYdwxWNd1V9CNNFLw0ml94ftWBsGCmD5/3+9wDdehfR4/HtM0Zm88GXH575zh8g8seetDT1CI58+fexUiyjsvfpOxK/nq3j1c2tkiHzTI2ZpqmbdFcFRCcHf2mvUy+pZ7z5R8DsVUyZferC6OaXuckNDi5Vhkve0IkrusVFzluGGlImZx/e83Qszi48M2GHeexcxVKc82XHtNjvh1zizOIpeeox3Usa/TteXjey6MAgcnQPWxX2d/3WsT1IUiZ1/HHhcXVbJFg1tUXZOWc5Y25EWXPvQCd2uJ3D+3sgoUS6TI9i0m+7T47F9seqqdft0CJ/jNCdtvPsXlH1jyrkceoxDP43tnqOqM7zz3LGeKKZ/fuY/t+YSN0ZKHT28zrwu+tX2aZeNsYEuwJZAm/LzPMTIW4fOpkC0VqcGVB32DEhLuOF5uRdY7gZhV5kvftd/HLfv+uNgP6kfwLBANhbqonW4EqtUsnlxZXqi5mme45hRn9hbo7vRARqhV1WXDrSrGhcUjNvi47ucwuSly3ytF2Phz00AJ4j2u3zgFXSCNt3lPtqiRZbQXDpYFed41PfU6amW/D074XdXTlmWc2DlD528TFw9oRw220s4+jRMHjgBkDr8x4eobNth57x4//Z2f4nO79/PlnXPU3vGWC5d40/ol/tezb+DZ3XVedeoq7zzzDc7lu3z62sNc2tkKL5aDRkJw3/e6x6cvV6QRir0gr1Ulm1fH9uFKSIiQY6IJT1aA76tJAketdUO2UBv88WLQ6qLj99UVoS+ri1m81JYha9Z0B+We8pzy/FtzsvICW59+Cp0tVjNVEePLmyZkkea7oo1fUQV1JmLOtOOBG5eoE4/fBwVIaFgKmbLszawRzEXXyaDAqW2XI2Xd6errYHkcaRXv26lRuAwGPdlib1KTuVC63jBwaOWMG2NktuzsDuLzqNq6XraPmTnwoKMhi/vXufymAaN3X+YXH/kkH99+E1/bOcvZ8Yy3nHqSc8UuH/3m23j2mVPcd982f/v0t3jd8Fn+ZnEfAGWZ2WvpNAT2ngQqUjOK6d4zJdtz5HPF1eAqxe0ujm17nJAAsH8YzZ2EHGhoOY6TEHkemAKXj/tc7hDO8cq5Vjj+6321qp4/xv+f8ArG1tp9+s7v+ukjHftHn/nFv7jZTNYXgxORwavqeRH5zK28sJOMV9K1wivvehMSDuCYEulkppqQkJBwu3GLdPCHufveCCnAJyQkJNxORA7+1piNfYSD7r7XxYmgaAIePe4TuIN4JV0rvPKuNyFhBbdKRaOqfyIiDx31+BMT4IP52CsCr6RrhVfe9SYkrOL22BAcBScmwCckJCTclVBeTIA/JyKf6f3+6LeTIB07By8i7wnj/R4TkQ8d9/ncDojI4yLyORH56/jiicgZEfmkiHw1fD993Of5UnCdkY6HXpsY/lN4rT8rIt99fGeekHAHcXQO/rKqvq339W3tfo81wIdxfr8K/CjwXcCPh7F/dyP+rqq+pScX/BDwx6r6CPDH4feXIz7CwaLP9a7tRzGH0UewYS//5Q6dY0LCseLEjuy7zXg78Jiqfl1VS+Cj2Ni/VwLeB/xm+Pk3gX98jOfyknHYSEeuf23vA35LDZ8CTonIxTtzpgkJx4hbJ5P8HeDPgNeLyJMi8lM3Ov64Ofj7gSd6vz8JvOOYzuV2QoE/EhEFfj1suy6o6tPh788AF47t7G49rndth73e9wNPk5BwtyK6sN6ShzrU3fe6OO4A/0rBu1T1kojcg3nof6n/R1XVEPzvOtzN15aQcGS8QjtZv/0Rfy8DqOql8P054PcxaurZSE+E788d3xneclzv2l4Rr3dCwgEc00Sn4w7wnwYeEZGHRWQAvB8b+3fXQETWRGQj/gz8MPB57Do/EA77APAHx3OGtwXXu7aPAf88qGneCez0qJyEhLsTShjXeYSvW4xjpWhUtRaRnwM+gQ3L+7CqfuE4z+k24ALw+8F/PQf+q6r+oYh8Gvi9UCT5JvBjx3iOLxnXGen4yxx+bR8H/j7wGDADfvKOn3BCwh2HHtvQmRNhF5yQkJBwt2JrcEG/596j1Ub/8IlfufvsghMSEhLuaiSrgoSEhIS7FCnAJyQkJNyNSGZjCQkJCXcnlG728R1GCvAJCQkJtxspg09ISEi4G3HrrApeLFKAT0hISLidUNBj0sGnAJ+QkJBwu3EbulSPghTgExISEm43EgefkJCQcBdCNaloEhISEu5apAw+ISEh4W6Eok1zLP85BfiEhISE24loF3wMOG4/+ISEhIS7H+qP9nUEiMh7ROTLIvKYiHzoRsemDD4hISHhNkIBvUUZvIhkwK8CP4TNNP60iHxMVb942PEpg09ISEi4nVC9lRn824HHVPXrqloCHwXed72DUwafkJCQcJtxC4us9wNP9H5/EnjH9Q5OAT4hISHhNmKXK5/43/rfzh3x8JGIfKb3+6Oq+uhL/d8pwCckJCTcRqjqe27hw10CHuz9/kC47VAkDj4hISHh5YNPA4+IyMMiMgDeD3zsegenDD4hISHhZQJVrUXk54BPABnwYVX9wvWOFz2mFtqEhISEhNuLRNEkJCQk3KVIAT4hISHhLkUK8AkJCQl3KVKAT0hISLhLkQJ8QkJCwl2KFOATEhIS7lKkAJ+QkJBwlyIF+ISEhIS7FP8fvEismQKvhp4AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 3 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bzUyStSFmu0p",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "81e6eeec-2f20-454b-e237-2458c23e8014"
      },
      "source": [
        "# Create training datasets\n",
        "batch_size  = 15\n",
        "data_train  = TensorDataset(u_train, data_obs_train)\n",
        "data_loader = DataLoader(dataset=data_train, batch_size=batch_size, shuffle=True)\n",
        "n_batches   = int(len(data_loader.dataset)/batch_size)\n",
        "# Info about the dataset\n",
        "print()\n",
        "print(f'u_train.min(): {u_train.min()}')\n",
        "print(f'u_train.max(): {u_train.max()}')\n",
        "print(\"data_obs_train.shape = \", data_obs_train.shape)\n",
        "print('n_batches = ', n_batches)\n",
        "#print()*7410"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "u_train.min(): 0.0\n",
            "u_train.max(): 1.0\n",
            "data_obs_train.shape =  torch.Size([20000, 1, 30, 183])\n",
            "n_batches =  1333\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dMpjvT8vbI5K"
      },
      "source": [
        "class TVM_Net(nn.Module):\n",
        "    def __init__(self, A, lambd=0.1, alpha=0.1, beta=0.1, eps=1.0):\n",
        "        super().__init__()\n",
        "        self.A = A\n",
        "        self.At = A.t() \n",
        "        self.lambd = lambd\n",
        "        self.alpha = alpha\n",
        "        self.beta = beta\n",
        "        self.shrink = torch.nn.Softshrink(lambd=lambd)\n",
        "        self.eps = eps\n",
        "\n",
        "    def name(self) -> str:\n",
        "        return \"TVM_Net\"\n",
        "\n",
        "    def device(self):\n",
        "        return \"cuda:0\" #next(self.parameters()).data.device\n",
        "\n",
        "    def box_proj(self, u):\n",
        "        return torch.clamp(u, min=0.0, max=1.0)\n",
        "    \n",
        "    def D(self, u):  \n",
        "        u = u.view(128, 128, u.shape[-1])\n",
        "        Dux = torch.roll(u, 1, 0) - u\n",
        "        Dux = Dux.view(128 ** 2, u.shape[-1])\n",
        "        Duy = torch.roll(u, 1, 1) - u\n",
        "        Duy = Duy.view(128 ** 2, u.shape[-1])\n",
        "        Du  = torch.cat((Dux, Duy), 0)\n",
        "        return Du\n",
        "\n",
        "    def Dt(self, p): \n",
        "        p    = p.view(256, 128, p.shape[1])\n",
        "        px   = p[0:128, :, :]\n",
        "        Dtpx = torch.roll(px, -1, 0) - px\n",
        "        Dtpx = Dtpx.view(128 ** 2, p.shape[2])\n",
        "\n",
        "        py   = p[128:256, :, :]\n",
        "        Dtpy = torch.roll(py, -1, 1) - py\n",
        "        Dtpy = Dtpy.view(128 ** 2, p.shape[2])\n",
        "        Dtp  = Dtpx + Dtpy\n",
        "        return Dtp        \n",
        "\n",
        "    def ball_proj(self, w, d, eps):\n",
        "        ''' Project w onto the ball B(d, eps)\n",
        "        ''' \n",
        "        dist = torch.norm(w - d, dim=0)\n",
        "        scale = torch.min(torch.ones(dist.shape, device=d.device), self.eps/dist)\n",
        "        proj = d + scale * (w - d)\n",
        "\n",
        "\n",
        "\n",
        "        #dist = torch.norm(w - d, dim=0) \n",
        "        #out_ball = dist > eps  \n",
        "        #proj = w.clone() \n",
        "        #proj[:, out_ball] = d[:, out_ball] + (w[:, out_ball] - d[:, out_ball]) / dist[out_ball]\n",
        "        return proj\n",
        "\n",
        "    def forward(self, d, tol=1.0e-3, max_depth=12, \n",
        "                depth_warning=False):\n",
        "      \n",
        "        self.depth = 0.0\n",
        "\n",
        "        # Initialize sequences \n",
        "        d    = d.view(d.size()[0],-1).to(self.device()) \n",
        "        d    = d.permute(1,0)         \n",
        "        uk   = torch.zeros((128 ** 2, d.size()[1]), device=self.device())\n",
        "        pk   = self.D(uk) \n",
        "        wk   = torch.matmul(self.A, uk)  \n",
        "        nuk1 = torch.zeros(pk.size(), device=self.device())\n",
        "        nuk2 = torch.zeros(d.size(), device=self.device())   \n",
        "\n",
        "        for _ in range(max_depth):\n",
        "\n",
        "            # TVM updates\n",
        "            res1 = self.Dt(nuk1 + self.alpha * (self.D(uk) - pk))\n",
        "            Auk  = torch.matmul(self.A, uk)\n",
        "            res2 = torch.matmul(self.At, nuk2 + self.alpha * (Auk - wk))\n",
        "            rk   = self.beta * (res1 + res2)\n",
        "            uk   = self.box_proj(uk - rk)\n",
        "            \n",
        "            res  = self.lambd * (nuk1 + self.alpha * (self.D(uk)-pk))\n",
        "            pk   = self.shrink(pk + res)\n",
        "\n",
        "            Auk  = torch.matmul(self.A, uk)\n",
        "            res  = self.lambd * (nuk2 + self.alpha * (Auk - wk))\n",
        "            wk   = self.ball_proj(wk + res, d, self.eps)\n",
        " \n",
        "            nuk1 = nuk1 + self.alpha * (self.D(uk) - pk)\n",
        "            nuk2 = nuk2 + self.alpha * (torch.matmul(self.A, uk) - wk)\n",
        "\n",
        "        self.depth = max_depth\n",
        "        if self.depth >= max_depth and depth_warning:\n",
        "            print(\"\\nWarning: Max Depth Reached - Break Forward Loop\\n\")\n",
        "\n",
        "        uk = uk.permute(1,0)        \n",
        "        return uk.view(uk.shape[0], 1, 128, 128)"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IZuPOkV4dokn"
      },
      "source": [
        "Setup Network"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mMcRhjnqm-Uy",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c59949f7-35e4-4aa2-bf39-d1681cb5ec26"
      },
      "source": [
        "Phi = TVM_Net(A.to(device))\n",
        "Phi = Phi.to(device)\n",
        "print(Phi)\n",
        "pytorch_total_params = sum(p.numel() for p in Phi.parameters() if p.requires_grad)\n",
        "print(f'Number of trainable parameters: {pytorch_total_params}')\n",
        "\n",
        "max_epochs = 2000 \n",
        "max_depth  = 250 \n",
        "criterion  = torch.nn.MSELoss()  \n",
        "fmt        = '[{:2d}/{:2d}]: train_loss = {:7.3e} | ' \n",
        "fmt       += 'depth = {:5.1f} | lr = {:5.1e} | time = {:4.1f} sec'"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "TVM_Net(\n",
            "  (shrink): Softshrink(0.1)\n",
            ")\n",
            "Number of trainable parameters: 0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_gvGr543dimO"
      },
      "source": [
        "Train the Network"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LEr2wb91nCoa"
      },
      "source": [
        "best_loss = 1.0e10 \n",
        "\n",
        "for epoch in range(max_epochs): \n",
        "  sleep(0.5)  # slows progress bar so it won't print on multiple lines\n",
        "  tot = len(data_loader)\n",
        "  loss_ave = 0.0\n",
        "  start_time_epoch = time.time() \n",
        "  with tqdm(total=tot, unit=\" batch\", leave=False, ascii=True) as tepoch:      \n",
        "    for idx, (u_batch, d) in enumerate(data_loader): \n",
        "        u_batch    = u_batch.to(device) \n",
        "        batch_size = u_batch.shape[0]\n",
        "        train_batch_size = d.shape[0] # re-define if batch size changes\n",
        "        u = Phi(d.to(device), max_depth=max_depth) # add snippet for hiding\n",
        "        output = criterion(u, u_batch)\n",
        "        train_loss = output.detach().cpu().numpy()\n",
        "        loss_ave += train_loss * train_batch_size\n",
        " \n",
        "        if idx%2 == 0:\n",
        "            # compute test image \n",
        "            Phi.eval()\n",
        "            u_test_approx = Phi(data_obs_test[0:1,:,:,:], max_depth=max_depth)\n",
        "\n",
        "            plt.figure()\n",
        "            plt.subplot(2,2,1)\n",
        "            plt.imshow(u_batch[0,0,:,:].cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u true train')\n",
        "            plt.subplot(2,2,2)\n",
        "            plt.imshow(u[0,0,:,:].detach().cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u approx train')\n",
        "            plt.subplot(2,2,3)\n",
        "            plt.imshow(u_test[0,0,:,:].cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u true test')\n",
        "            plt.subplot(2,2,4)\n",
        "            plt.imshow(u_test_approx[0,0,:,:].detach().cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u approx test')\n",
        "            plt.show() \n",
        "            Phi.train() \n",
        "\n",
        "        tepoch.update(1)\n",
        "        tepoch.set_postfix(train_loss=\"{:5.2e}\".format(train_loss),\n",
        "                            depth=\"{:5.1f}\".format(Phi.depth))\n",
        "    \n",
        "\n",
        "  loss_ave = loss_ave/len(data_loader.dataset)\n",
        "  end_time_epoch = time.time()\n",
        "  time_epoch = end_time_epoch - start_time_epoch\n",
        "  #lr_scheduler.step()\n",
        "  print(fmt.format(epoch+1, max_epochs, loss_ave, Phi.depth, \n",
        "                   0.0,\n",
        "                   time_epoch))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iyxCl1Id_Fwb"
      },
      "source": [
        "from skimage.metrics import structural_similarity as ssim\n",
        "from skimage.metrics import peak_signal_noise_ratio as psnr\n",
        "\n",
        "n_samples = u_test.shape[0]\n",
        "\n",
        "data_test     = TensorDataset(u_test, data_obs_test)\n",
        "test_data_loader     = DataLoader(dataset=data_test, batch_size=batch_size, shuffle=True)\n",
        "\n",
        "\n",
        "def compute_avg_SSIM_PSNR(u_true, u_gen, n_mesh, data_range):\n",
        "    # assumes images are size n_samples x n_features**2 and are detached\n",
        "    n_samples = u_true.shape[0]\n",
        "    u_true = u_true.reshape(n_samples, n_mesh, n_mesh).cpu().numpy()\n",
        "    u_gen  = u_gen.reshape(n_samples, n_mesh, n_mesh).cpu().numpy()\n",
        "    ssim_val = 0\n",
        "    psnr_val = 0\n",
        "    for j in range(n_samples):\n",
        "        ssim_val = ssim_val + ssim(u_true[j,:,:], u_gen[j,:,:], data_range=data_range)\n",
        "        psnr_val = psnr_val + psnr(u_true[j,:,:], u_gen[j,:,:], data_range=data_range)\n",
        "    return ssim_val/n_samples, psnr_val/n_samples"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xrt0leh7qg59"
      },
      "source": [
        "\n",
        "test_loss_ave = 0\n",
        "test_PSNR_ave = 0\n",
        "test_SSIM_ave = 0\n",
        "with torch.no_grad():\n",
        "    for idx, (u_batch, d) in enumerate(test_data_loader): \n",
        "\n",
        "        u_batch    = u_batch.to(device) \n",
        "        batch_size = u_batch.shape[0]\n",
        "        temp       = u_batch.view(batch_size, -1)\n",
        "        temp       = temp.permute(1,0)        \n",
        "        test_batch_size = d.shape[0] # re-define if batch size changes\n",
        "        Phi.eval()\n",
        "        u = Phi(d, max_depth=max_depth) # add snippet for hiding\n",
        "        output = criterion(u, u_batch)\n",
        "        test_loss = output.detach().cpu().numpy()\n",
        "        test_SSIM, test_PSNR = compute_avg_SSIM_PSNR(u_batch, u, 128, 1)\n",
        "        test_PSNR_ave += test_PSNR * test_batch_size\n",
        "        test_loss_ave += test_loss * test_batch_size\n",
        "        test_SSIM_ave += test_SSIM * test_batch_size\n",
        "\n",
        "        print('test_PSNR = {:7.3e}'.format(test_PSNR))\n",
        "        print('test_SSIM = {:7.3e}'.format(test_SSIM))\n",
        "        print('test_loss = {:7.3e}'.format(test_loss))\n",
        "        if idx%1 == 0:\n",
        "            # compute test image \n",
        "            plt.figure()\n",
        "            plt.subplot(1,2,1)\n",
        "            plt.imshow(u_batch[0,0,:,:].cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u true')\n",
        "            plt.subplot(1,2,2)\n",
        "            plt.imshow(u[0,0,:,:].detach().cpu(), vmin=0, vmax=1)\n",
        "            plt.title('u approx')\n",
        "            plt.show()  \n",
        "\n",
        "print('\\n\\nSUMMARY')\n",
        "print('test_loss_ave =  {:7.3e}'.format(test_loss_ave / len(data_loader.dataset)))\n",
        "print('test_PSNR_ave =  {:7.3e}'.format(test_PSNR_ave / len(data_loader.dataset)))\n",
        "print('test_SSIM_ave =  {:7.3e}'.format(test_SSIM_ave / len(data_loader.dataset)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qDWZLuTs_Dpl"
      },
      "source": [
        "if use_ellipses:\n",
        "    ind_val = 0\n",
        "else:\n",
        "    ind_val = 1000\n",
        "\n",
        "u = Phi(data_obs_test[ind_val,:,:,:], max_depth=max_depth).view(128,128)\n",
        "u_true = u_test[ind_val,0,:,:]\n",
        "def string_ind(index):\n",
        "    if index < 10:\n",
        "        return '000' + str(index)\n",
        "    elif index < 100:\n",
        "        return '00' + str(index)\n",
        "    elif index < 1000:\n",
        "        return '0' + str(index)\n",
        "    else:\n",
        "        return str(index)\n",
        "\n",
        "cmap = 'gray'\n",
        "fig = plt.figure()\n",
        "plt.imshow(np.rot90(u.detach().cpu().numpy()),cmap=cmap, vmin=0, vmax=1)\n",
        "plt.axis('off')\n",
        "\n",
        "data_type = 'Ellipse' if use_ellipses else 'Lodopab'\n",
        "\n",
        "save_loc = './drive/MyDrive/FixedPointNetworks/Learned_Feasibility_' + data_type + '_TVM_ind_' + string_ind(ind_val) + '.pdf'\n",
        "plt.savefig(save_loc,bbox_inches='tight')\n",
        "plt.show()\n",
        "\n",
        "print(\"SSIM: \", compute_avg_SSIM_PSNR(u_true.view(1,128,128), u.view(1,128,128).detach(), 128, 1))\n",
        "\n",
        "############\n",
        "# TRUE\n",
        "###########\n",
        "\n",
        "cmap = 'gray'\n",
        "fig = plt.figure()\n",
        "plt.imshow(np.rot90(u_true.detach().cpu().numpy()),cmap=cmap, vmin=0, vmax=1)\n",
        "plt.axis('off')\n",
        "\n",
        "\n",
        "save_loc = './drive/MyDrive/FixedPointNetworks/Learned_Feasibility_' + data_type + '_GT_ind_' + string_ind(ind_val) + '.pdf'\n",
        "plt.savefig(save_loc,bbox_inches='tight')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}