Scalable Epistatic Net (EN-S).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script implements the EN-S regularization algorithm. \n",
    "\n",
    "Aghazadeh A, Nisonoff H, Ocal O, Brookes DH, Huang Y, Koyluoglu OO, Listgarten J, Ramchandran K. Sparse epistatic regularization of deep neural networks for inferring fitness functions. bioRxiv. 2021 Jan 1:2020-11. doi: https://doi.org/10.1101/2020.11.24.396994\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from itertools import product\n",
    "import random\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "from torch.utils.data import TensorDataset, DataLoader, random_split\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import seaborn as sns\n",
    "from utils import *\n",
    "\n",
    "sns.set_style('white')\n",
    "sns.set_context('talk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def walsh_hadamard_matrix(L=13, normalize=False):\n",
    "    '''\n",
    "    Compute the WHT matrix for domain of dimension L\n",
    "    '''\n",
    "    H1 = np.asarray([[1.,1.], [1.,-1.]])\n",
    "    H = np.asarray([1.])\n",
    "    for i in range(L):\n",
    "        H = np.kron(H, H1)\n",
    "    if normalize:\n",
    "        H = (1 / np.sqrt(2**L)) * H\n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>binary</th>\n",
       "      <th>amino acid</th>\n",
       "      <th>counts</th>\n",
       "      <th>counts.1</th>\n",
       "      <th>counts.2</th>\n",
       "      <th>Unnamed: 5</th>\n",
       "      <th>brightness</th>\n",
       "      <th>brightness.1</th>\n",
       "      <th>Unnamed: 8</th>\n",
       "      <th>brightness.2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>genotype</td>\n",
       "      <td>sequence</td>\n",
       "      <td>input</td>\n",
       "      <td>red</td>\n",
       "      <td>blue</td>\n",
       "      <td>NaN</td>\n",
       "      <td>red</td>\n",
       "      <td>blue</td>\n",
       "      <td>NaN</td>\n",
       "      <td>combined</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>'0000000000000'</td>\n",
       "      <td>DVLTFNSAAYNNK</td>\n",
       "      <td>5431</td>\n",
       "      <td>12</td>\n",
       "      <td>7846</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.085317</td>\n",
       "      <td>1.574628</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.574635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>'0000000000001'</td>\n",
       "      <td>DVLTFNSAAYNNR</td>\n",
       "      <td>6574</td>\n",
       "      <td>11</td>\n",
       "      <td>9047</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.075637</td>\n",
       "      <td>1.544269</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.544276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>'0000000000010'</td>\n",
       "      <td>DVLTFNSAAYNKK</td>\n",
       "      <td>10493</td>\n",
       "      <td>33</td>\n",
       "      <td>13352</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.103633</td>\n",
       "      <td>1.49045</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.490458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>'0000000000011'</td>\n",
       "      <td>DVLTFNSAAYNKR</td>\n",
       "      <td>9545</td>\n",
       "      <td>35</td>\n",
       "      <td>12513</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.106846</td>\n",
       "      <td>1.511977</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.511984</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            binary     amino acid counts counts.1 counts.2  Unnamed: 5  \\\n",
       "0         genotype       sequence  input      red     blue         NaN   \n",
       "1  '0000000000000'  DVLTFNSAAYNNK   5431       12     7846         NaN   \n",
       "2  '0000000000001'  DVLTFNSAAYNNR   6574       11     9047         NaN   \n",
       "3  '0000000000010'  DVLTFNSAAYNKK  10493       33    13352         NaN   \n",
       "4  '0000000000011'  DVLTFNSAAYNKR   9545       35    12513         NaN   \n",
       "\n",
       "  brightness brightness.1  Unnamed: 8 brightness.2  \n",
       "0        red         blue         NaN     combined  \n",
       "1   0.085317     1.574628         NaN     1.574635  \n",
       "2   0.075637     1.544269         NaN     1.544276  \n",
       "3   0.103633      1.49045         NaN     1.490458  \n",
       "4   0.106846     1.511977         NaN     1.511984  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_excel(open('41467_2019_12130_MOESM7_ESM.xlsx', 'rb')) \n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/0m/znv2x4f15l370c2fvtf4t8kh0000gn/T/ipykernel_75286/1444983188.py:6: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.\n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  brightness = data['brightness.2'][1:].to_numpy(dtype=np.float)\n"
     ]
    }
   ],
   "source": [
    "# get genotypes as strings. the excell file has some odd formating\n",
    "# so we remove the first row and get rid of extra quotes at beginning\n",
    "# and end\n",
    "genotypes = data['binary'].to_numpy()[1:]\n",
    "genotypes = [x[1:-1] for x in genotypes]\n",
    "brightness = data['brightness.2'][1:].to_numpy(dtype=np.float)\n",
    "\n",
    "genotype_to_brightness = {g:b for g,b in zip(genotypes, brightness)}\n",
    "\n",
    "# get dimension of input\n",
    "L = len(genotypes[0])\n",
    "# size of the landscape\n",
    "N = 2**L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all 2^L inputs\n",
    "all_inputs = np.asarray(list((product((0,1), repeat=L))))\n",
    "# convert to strings\n",
    "all_inputs_as_str = [''.join(map(str, x)) for x in all_inputs]\n",
    "# get y-values for all possible x-values\n",
    "landscape = np.asarray([genotype_to_brightness[g] for g in all_inputs_as_str])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1399e5c70>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAEDCAYAAAALG2EMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABYF0lEQVR4nO3dd3QU1dsH8O/uppEEQhJ67y2hN+lSQ1MQC2AB64sFFREEFEUsFAUUOwjID0FUQEHpRVpCTYBQA4QkJAHSe99y3z822exs383szOzm+ZzjMTN7d+ayM3PnmdtGxhhjIIQQQgghgpCLnQFCCCGEkOqEgi9CCCGEEAFR8EUIIYQQIiAKvgghhBBCBETBFyGEEEKIgDzEzoAtOnXqBI1GA39/f7GzQgghhBBiUUFBAeRyOa5fv27yc5eo+dJoNKAZMQghhBDiChhj0Gg0Zj93iZqvihqvyMhIkXNCCCGEEGJZr169LH7uEjVfhBBCCCHugoIvQgghhBABUfBFCCGEECIgCr4IIYQQQgREwRchhBBCiIAo+CKEEEIIERAFX4QQQlza/ZxinI3LFDsbhNiMgi9CCCEurf+y/zB57RmcuJUudlYIsQkFX4QQQtzC7sv3xc4CITah4IsQQohboLfQEVdBwRchhBBCiIAo+CKEEEIIERAFX4QQQgghAqLgixBCCCFEQBR8EUIIIYQIiIIvQgghboEGOxJXQcEXIYQQQoiAKPgihBBCXER8RiH+vpgMlVojdlZIFXiInQFCCCGE2GboimMAgMyCMrw8qJW4mSEOo5ovQggRUalKLXYWiMjUGvt7q+2+/MAJOSFCoeCLEDtlFpSixfw9aDF/j9hZIVXEGMOCv67g74vJoux/e1QyQhcdwLbIJM56jYYhOikHZSpqWrKHK75e6KNdV9Hj00OISy8AAETdzcbVe7ki54o4GwVfelJyS1CipKdQYtmQL4/p/s4rUfK+/cTMIkz4Lhy/nU3kfduE66Nd17D1XCLe+SMaxWXCX/tztkVDqWaYu/0yfj2dgKzCMgDAlwdvYsL3EXhr60XB81RdZReW4dfTCUjPLxV0v5tO30VusRIf7rqKpKwiPP7jKYz/Nhy5RfyXLUQ6qM9Xudi0Aoz66ji6NwvEjtf6i50dImEFpSrd3xoHmguseXfbJUQn5yI6+Qqe7tuM9+2TSr+euav7u0ytgadahm/+i0W3pgFIzy9Fi2A/9G0VLEhePtx1Db+fT8LMoW3w47E7AID911IE2TcBZmyOwrn4LPx65i4OvjPE6fvbEB6PyLtZuuWI2Ew89sMp3XJyThECfAPMft8FK/mIHgq+ysWlF0DDgJsp+Zz1jDHsuHAPbev5o2vT2uJkjrg1pVqD6/fzENo4AAq5DPdzSnSfaTQMcrlMxNxVHx5yGTafuYtvjtzmrE9YNk6wPFy7n4fXtlwQbH+k0rl4bSB0K7XA6ft6kFuMT3ZfN1qfUVBZ65ZfojL6nMMV21iJDjU7llOU3+D0Oz5uCI9HywV7MWdbNCZ8H+H0PBSXqXU1KWn5JRjy5VEs/vea0/dLHHcvp7jKHaYX/n0VE76PwJcHbhptq9snB/H14VtV2j6xjVwmw2Xqa1NtXL2XiylrT+PErXSb0j/ILcbH/1zDjQd5Vd63LU2bqw7Rde/OKPgqV1G7UKxUo6hMBZVaY/LJxFlO3EpHx4/2o9X7ewEAKw7cxN3MIvwSkSBYHtxRdmEZ/oxMQm6xc/pPjPsmHI9+W7XA/I/yztY/Hb+Dzh8fxL2cYt1neSUqfH34trmvkiqq5WNQ+U+VCdXGU2tO40xcFqZtOGdT+lc2RWLjqQSMWX2yyvuWy6zXZidmFlV5P0S6KPgyYcuZRKfdrM1563dux1r9KueJ30dgyd4bgubHXTy/8Tze234Zb1hoyknLL8H8HZdxNi7ToX3cTM23nshGNLpNWCM61tf9Hfb1Cfx18Z6IuSFVxeyInovsHGBx9V7Va7zsobbSrEjPCa7N7uDrxo0bCAkJQUqK5Y6gu3btQvv27Y3+++STTxzOrDPpd5zOKCxFQmahUZq/LyZj1yXnFM51/b11f687GYd9Vyt/30tJOVh7Is4p+3V30Uk5AIDw2Ayzad79Mxq/n0/C5LVnAGibI1rM34NXNkUKkUUiIv0bWGIW1TQQYdhS8+XI3F/EddjV4f7OnTuYMWMGVCorHQEBxMTEoHnz5vjiiy846+vUqWNfDgWi0jvR1xyPw5rjxsHOO39EAwB6NAtE0yBfXvfvqaiMgz/b4161XNsikxDaOAAdG9YSOysmRd3N5iyP/zYcAHDoeioKSlXw9pBj7rZotK1fE28MbSNGFjmik3Lwye7reHNYGzzcvp7Y2XFpzM5Oy9ujkvHziTh8NbkbOjWS5vlc3anUGrz9+yW0rOOHOWHtxc6OSXIbqj2svT6I+tu7NpuCL5VKhT/++AMrV66Ep6enTRu+efMmQkJC0K1bt6rkTzD2PGWk5pWgaZAv4tILwKCtNSsqU1dpNKTChhFt2yKT0KdlEJoH+zm8H6G9tz0af0ZqJ7CM/XwMPBTSa+m29Mvfyy5GfEYhdl66DwCY3r+FIHmy5ImfTkGpZnj+l/NIWDYOaXklqOPvTaMiHWBv5cKcbdoHsBc3nseZ94c7IUekqnZeuo89V7Szv78xtA1qeClEzpExmcVSR4tqvtybTXfCqKgorFixAi+++CLmzJlj04ZjYmLQvr00nzpMUdlxoqs1DKl5JRi28jiGrzyOkV+dwITvI7AhPN7h/dsSfM3dfpkzwacrqAi8AECplmZhYqkJ4Fx8JlLzKqd+KCqzXuvrDN8fjdVNg6L/O64Pj0efJUfwzp+XRMmXq3P0jEzROyfEcDezEL9ExHPmnCNa2eUT1QKASiPNPpS2PCcpKfhyazYFX61bt8bhw4cxc+ZMKBTWnyLS0tKQmZmJ69evY/To0QgJCUFYWBh27txZ1fw6jdqOi1TDgDMmOmd/svs6HuQWm/iGdbYEX67Ohm4OotDP15S1pzmfZRaW4a8LlQGkWB3ivzxwE2FfnzBa/2n5iNxd5TVzxD57rzj+frwSpRpX7+Xa3XRpr0K9AOu97dFoMX8Phnx5DIv/vY6Pdl116r5dDuNez5kFZebTSpy1ssaewQVEemwKvurUqYPgYNtneY6JiQEAJCcnY+7cuVizZg06d+6MefPmYceOHY7l1MmsNK9zaCwUtvoTZNpD4UBkcp+HOabcTZlKgwPXUjiTFVoj0/vtz8RlcT77+vBtRCdXzv0k0Qdp4iBbmnaOxqRxaj8rTFl7BuO/DcdmvVnynaHbJweRVD4YQL8mGQD+uiDd0ZkPcoudHpiaon89z9tx2abvCN3EZ20kI3F/TumAExoaip9++gmbNm3C8OHDMXDgQKxYsQL9+/fH6tWrnbHLKguoYVtfNkAbfPF97dgbe11IzEb/Zf9h3Dfh/GbExX11+BZm/BqF8Xb8LlKtkSPS8MLG8xj0xVGj9ZfKR9J+eeCmU/evVDN891+sU/fBt81n7qLf0v/wuQiDh/QvZ8PBNOZ0/vgAZ3ljRDzumhjxzhd77jfmlCjpSdCVOSX4CgoKwtChQ+Hv789ZP2TIEKSmpiIrK8vMN8UzvIPto8ZyipQWapwci8rsbXasePdbbJrzX4XhStaf1Pa7s6dPDsVexBqx51/bfiHZeiIJWbhT2xy6rgr9YB2l/zBl64OV4ZxfH/97HcNXHucxV/y4qvcGhti0Artq+Im0OCX4unjxIrZt22a0vrS0FB4eHqhZs6YzdlslcrkMDWr52JR28b/X8M0Rfp9Eq0OfL6nWtNsy5w5xT0Yz3DtAiNOaRr7ZTv9qrkqZo9IwPLf+rK6GUwpO3ubOV1iVQV5EXE4Jvi5duoSFCxfq+n4BgEajwYEDB9CjRw+bp6sQmq3xT0ZBGecVMPzs2/YAIC2/BCVK6utlSpk9nffK2RN7mevkmpYv7ug34phAPy+b097m8U0GjnhqzWnriQinz5c9o9hNOXk7AxMFeK+vOZeTc1BcptaV94ZlFfUdc128BF9ZWVm4dOkSCgq0TWCTJk1Co0aNMHPmTOzevRtHjx7FjBkzcOvWLZunqhCDjJcaEMe2Yc8l1OfzI0ZPQK5AuqNzqn7chep756w3LFRX9ty7Rn5lPNpUuxF+8mLNuXjpddeoUKJUY//VB5xpHsTA4F59OB/9LgIdP9qPvkuOoFSlNiqp1BKdvodYx0vwdezYMUyePBnXrl0DAAQEBGDz5s3o0qULli5dilmzZqGoqAgbN25E165d+dglcUFSfUizp7A2N89aer4wfS/e/v2SIPshxB4Tv4/Aq5svYPRqMwGqgNwo9tLJLVZiz+UH2HoukbO+qjV7RDx2d3iYNGkSJk2aZHVd48aNsWrVqqrljhBCJI5uf0BM+QTAqXnUAdxZZv8ZbbROjKk8CD+k964Xlyf8xWDtHWBSIdViwh2flIltpNsUThymV5XdjOd38EoNnb2ui4IvN9Dmg32Yt922yQSdLauwDCsO3MSNB3lGn0n1KS1NoCZD4p7oFT/iSsws0v2dUVDKeZjyleB7HQkBKPiSjKoGJn9EJvGUk6p5b3s0vjsaizGrT4qdFVLNJWYW4c/IJHoLhJvTnwPt5O0MTllaHabwIa6p6pPcEKLn8I003d95JUoE+noiu0gJgKrIibAGf6mdlf5uZiHmhnUQOTfEWQzjK/0+6HzN4TffxtcU2YqvRgCJNiYQG1DNF+/oSatCYakKHRrU0i27e0GRWx5kEmn5W8LvPyRVZxhg6b97l6+Kr9/PS6NlgVimVGvw7p/R+CXC/OSzSrUGGgmMEqXgy42cinW9ub+cacFfwvaDyyuh4EuKLBWz7v5AIIaxAnY5+PN8ElYdumU+gZMm/dp9+b4kmrNpwAjXb2cTseNCMhb/e93k5/klSgxc/h8mrz0teh9kCr54J94BfXrdWdH2bYrMsBZQ4J9m6zlhn1Y9FXQ5EXLdxGAbZ3nPRHPgT8fv6P52VjvEzN8uotdnh520ddvRwwOXtTfP7IhKRmpeKc4nZCO3WNyHZerzpcedZkaWApVGw/lNnfmUlpxdhP+dSnDa9m1BnXuliW5Q1YtQc43ll4g/ypVObS792iylWmP0QKzf2mhUOSAwelQnTiVUQDtt/Tn8fFLcl8x6UPDlcigwc2+hjWtZT0Tcxp30Qt3f/8WkWUgpPgq+iGCceaOLyyi0nogQUq3Ur+kjdhacih4euJR6E46rJP7eSwq+iGCkfSkQQoiroVJVn34waqrVRUq/FvX50kN9vvjFmPjt6oTQiDBC3FtRmQqrDt5CuAuN+KeaLyIYsYf2kurJ3U+7nKIy0fYdm5aPN7ZcEGXfUninbZlK3DyUqdz85LbR6iO3sS5c3D6/9qLgi/DKy4N7SlFtIiHOdeh6qmj7fvrns9hz5YEo+84qFC/orLDy0E1R90+1ulrRSTlG66R+66Hgi/BqXOeGnGV3r3XQV43+qS7F3Y/LehGf+KX+UnpnH/tfwhMc+p67n5NCc8XuLRR8EcFQgUMI/2JS8sXOQvUl9j2fClUADrawiHzsKPjS44rRs5QUlRlPOkjNjkRs1an2lXA5u/gRu3jT0MkNwDXvMzTakWcSeF+nKPZfTcFrW6LoRkckyPxJSYNASFWIfdOns1fLVMWJ2MfGGgq+CC9e3RxltM7wvkb3OUIIX6g4oTK1gtQDLVOo2ZE4zcJdV3HytuvMu1JVVIsiTZYOCx0x1yWFJjexu6pI4TcgjqHgizjNiVvpYmeBEAqwCE7FZmBbZBLv2xW7xoXObS2ZjQdCSg/I1OyoR+wLyd3RnDREDHRZV29qDcPT684CABoE+GBQ27q8bVv0c4uKVAAOHgeRfzuq+SKEuDVLZayEHoSJEzAAKk3lLPSXk3PFy4wT0AOt66LgSyLoJkAIIbaztcx0Zr8sW5u7nIXuG1omX6Jt5bcRO3Cl4IsIx80LCjf/57kl6mrg/px5k3W0DxFffY8o+NJyxcuYgi89rngACSHGSpRq3d+W3gFINy/35uwyvbBMbT2RE8VlFIi6fylQqjVQO3Adi33tU/DFM7EPKLHNfzGpGLriGI7TiEy3lCmBly4T56KiFkjOLhY7C6JSqTUI+/qEQyPrxT5/KPiSiOrQ/MGgfQXRkz+dwqe7r4ualxc3RiI+oxDTN5wTNR/EOarB5URsRA/E7isuoxBx6YUmP5P6PZWCL4moLgXELxEJOJ+QjfXh8YLve8ra07icnCP4fonw5LbO+yP682/VdWtaW+wsuAwpzfNEqk4hN3+dW+1wL/K5QMEXEVReiVK0fZ+Jy8Kj30Xg6M00p2yfynXpsFAmux1XOO1Sckvw1aFbSMoqEnS/hr/NgWupgu7f2WiGe9dFwZeehExhC4bqSOzXcQDAC7+cFzsLxNnMnGY5RWU4dD0VpSptR2m6dwlj2oazWH3kNh774RRv27S15kI/2ZV7udBoGH49c5e3fBDxVOVuIvalT8EXEQxj0m+HJ+7BXLPjlLVn8MqmSCzbFyNwjpzHFS6pW6naUXkZBaUi5wTYFX0PH+68KnY2SDVHwRcRjdht7sR9mQq+7mYWIiYlH4C27yGpHgz79e27kiJSTvhHRajjxP7t6N2OEuEOHX9toX9LpJow4iymTqshXx4TOhuCqB4lh2NWHbqFVYducdbdz3Wf6Rmq+7G39IYBqf82VPNFBMPAOMGWu3UWrS4BtK1yisrw7ZHbuJMu/ESQFNS7P0eLDyn0O3Wzos8liV1eU/DFM2pKsx39Uu7th2N3sPLQLQxfeVzwfUvhBkvMu5tZiP1XH0CjEb4UuHLPfV6uTWe566LgiwhK/6ZIcap7E2MuN3vRKSiOIV8ew6ubL+DPyCSxs0JcWJUqO+j1QgSoHk/qhtdJQakKag2DSq2hGkM35P5ntERI/NqxdG3TwAdSFf9GPxA7Cw6jDvcSIXb7s1D0++L0+PQQ2tevieyiMrRvUBO/vtRXvIwRQpyix6eHzH52MzVfwJy4n8aBNcTOgpGKYNtSZ3i+HL5hftJcqzPc85wXe1HNF+HFgDbBNqUzvBxvpuYjLb8UJ29nQKXW8J8xIYl9NRMA2pftFivVNqWVeKWRW8guEu+tFu7OUy6tWzhjDFN/PoOJP5wSpDyXu/CrLKjmi/AioIan1TTW7nNCPCkR4YhxOBljGLP6JG6nWR9heSe9QBKTflY3haUqsbNAnOTGg3ycicsCAJyMzcDQ9vWcur8qzXBPfb5IdZGWVyJ2FoiAlGrhS7eswjKTgde9HOO5ncQYhUmAbp8cFDsLbkNq3VXUeqNXVXrX//9OJWDruUQxsiRZVPNFBGPtvW7avgJU+0X4l11YJnYWSDm+gnKxay6kQGq/gana7hsP8rDon2sAgIFt6qBpkC9v+6tKq6PYgSvVfPFMYtcCIdWKqabrlnX8sHy/+7zLkZAKUr7fVHS8f6D3RoEsnh+C7O2qIqVglWq+JEJKJ4UjqsNUGda4+CF0eSm5JbjxIM9ofZPAGjh5O0OEHBFSfZmaP1dKZaTY91wKvohkSOnCJK7noaVHTK4vKrNt5CNxLWI3G0mB1OZH1K+I2nQ6AaNDGzj1wdzeLevnT+xfjpodCSFuzd3rZMW+iRDxSPnYn7qTabROyGAx2N9LsH05goIviXD5WRZcPf/Ebbn8tUWIGR0b1hI7CxaVKNUoVTmv5rkq17bYtYbU7EgIIYS4oC6NA8TOgkUDl/+HjILKTvZChjumYisptdJSzZdESOmkcAQflQuu/hu4ev6FUFSmwq5L9wSd+oGOi3ui4yq9ZkfD/l36gRcAZBbwPNqxCncesc8fCr4IIYKZv+MK3v79Ep5edxaMMeSX0KtnCHEX1poBF+26KkxGIP0BGXYHXzdu3EBISAhSUlIspissLMTixYsxYMAAdO/eHa+88goSEhIczSchxA38E30fgHbixU92X0fnjw/i5O10p+5T2kUwIdXH/dwSnLjl3Ovd0KnYDFxKyhF0n7awK/i6c+cOZsyYAZXK+ru53nnnHezfvx9z5szB8uXLkZqaimnTpiE/373fYi92VaYrk/qTCuHXLxEJAIDXN19w6n7E7lhLSHVhSwf4aRvOIfx2BnL5eOG6pf0xIDYtH0+vO4uJ30fgQmJ21ffHI5uCL5VKhS1btuDJJ59Eaan1F9FGRkbi+PHjWL58OR577DGMGjUKGzduRH5+PrZu3VrlTBNCiK0o9HJPdFxd92H/2fVn0fWTg4g18R5WPl1OztX9PemHU7iZWln5I/ZvZ1PwFRUVhRUrVuDFF1/EnDlzrKaPiIiAn58fBgwYoFsXFBSE3r1748SJE47n1o2JfSJUlb2veXBHVHMnTa5+bVnj7v8+4r5+OBbrtG0zAKUqDWfdgWuWu0sJyabgq3Xr1jh8+DBmzpwJhUJhNX1cXByaN29ulLZZs2aIj493LKeEEJdVolRj/9UHouybYhMiNe4aMEvtNXM/Hb/DWdbPndgPyzbN81WnTh27NlpQUAB/f3+j9X5+figocG41IxEHTTVBLPn4n2v4/XySyc/osBPiGLEDCLFZ7PLFgLuZRdx1zs2OXZwy1YSlDq5yOc1uYQq12hF3Zi7wEkJGvvV+qkQctXwcn+ebBlJIjyvdx8Q+fZwSCfn7+6OwsNBofWFhockaMUJI9VVQan30dFXcyyl26vYJIY6578Rrs1hp4rVGEorXnRJ8tWzZEklJSUZPJnfv3kXLli2dsUuXJ3YUTqqOjqHjxOoPRsRFA3WqtzNxWTanLTEIpm48yMPZePPf/8NEbXu+3oOe2MW1U4KvgQMHIi8vD6dOndKty8rKQmRkJPr37++MXUpGdW+DJ8QRrzp5ri8iTdR0WDVS+/mcFUofvJaCzh8fwM8n4gAASrUGY1aftPidwzdSnZQbfvASfGVlZeHSpUu6zvS9e/dGnz59MHv2bGzbtg2HDh3C888/j5o1a2Lq1Kl87JJIDD3AEkKEJLG4g8B594H/+zUKSjXD53tvAACKykw0Kdrpyr1c64mciJfg69ixY5g8eTKuXbumW/fdd99h2LBh+OKLLzB//nw0aNAAGzduRECAtN/CTsQjtac4Qggh7umtrRdF3b/dQ00mTZqESZMmWV0XEBCApUuXYunSpVXLIZEktYbh6r1cdGpUC54KGsFKiFiqS1eHEqUaa0/EoV/rYPRuESR2dohJwjSBpOWVCLIfZ6K7pkS4WgH6xYEYTPg+ArP+uARAqEtO2lzrCBIivrwSFfZesW2wxfdHY7Hq0C08+dNpAFRTXp39fDJO7CxUGQVfxCFrjmtP/j2XaZQaIcRxr2+xbbDF6TuZTs6JcFztYVtq3CHwpuCL8IKPIeN3s4znhiOEEIBqll2BvbeBZx9q5tB+3OFcoOBLIqT2Tix78TFk/IkfT/OQE0KIO6JpKdyPRzV+4031/ZcTyTGc6fyz3dcxbMUxpJt4PczaE3fw/dFYobJGCBGZcehFwZjUAlLXrkIQFgVffHPwWnD1PgDOmKl6XXg84jIKseLATc76+IxCLNkbgy8P3MS1++LO1aJPagWhu9p75QGOxqSJnY1q6fSdTNzNFKd7gIYuL8mz9z7gaJnpDkUtBV/Equ1RyVi48wrKVBrB9nlVbwK8nOIyzmdZhZU1YaZeIUHcV0puCV7fcgEvbDyPxMwisbNTrZyJy8TUn89gyJfHxMmAO9xxCS9cvbICcGCeL1L9zNkWDQBoEeyHlwe1MpmG73qv8NgM3d+WZjO+kJjN856JlCXo1bok5xShWbCviLmRBqFikiNOfl3LvZxi+Ht7IKCGp8nPXf92SwxV53d7UvBFbJYgYHODQu+iVKq5NW76Nxt5Nb54papEqcaNB3m8bY8xBg0DFHKDY013Y0E580aZkFGIh1ccAwDELx1rcl+GQSZVhEnvN6DS2HYUfBGzPvn3OsJj03XLFi90nq46xhhkMhnkejfavGJuR3z9vh90sUvPc+vP4nwCfzWSr2yKxNV7eTgwazAdbxE587f/9cxd3d/dPz2Ens0Csf753pw0+u/iyy1SOjE3xFH2xucO95OVWNDpCOrzRczaEBGPW6kFumUhznd1eWSVX1JZuMakcGtRNPoXrIRqvqT2FOosjDFkF5aZ/ZzPwAsADt9IQ0peCX45Fc/NB697IQCg0TAUm2vmd+Kltj688tjmFClxJCYNhQajn/UVK6v+YmXCPyGmTHpu/VncSMl3+n6cjYIviXCFG/eJW+nWE1XRo99FgDGG1LzKTvUapm3KqhCZkKX727AlijhfywV70f3TQzguwPmg/2SsoeFuJvFZdjy/8Tx6fHoI93OK8cHfVzBnW7TuGOjfWIUY2Xudx6ZrqXGF8l6qTt7O4LVbg1go+CI2S84uRvjtDKP1fI6CvP4gDwevp+KiQUf6Dh/uR9Rd7bptUcm69dTnS1jn9QLf6RvOOXVfjDFcSspx6j4I14lb6ShWqvHW1ovYcjYR26OSEVl+3elfamLHwTIZ1XwC0vsN7C2Oq3OHewq+eOboxSC1c9Cwqa/CkRjjEU/tFu7D2bgsE6kdM+PXKMSYqFZ+/MdTAAw73PO2W2KDB7klgu1rx4V7eOyHU5UrZDJOYa1/HlTnOdacUXboX38VzX/615oQv/eTP53Gp7uvm/yMLnv3deJWOiavOc2ZbsgdUfBFAAAHrqXgw51Xdf09Rn990mQ6c2XuvZxiZ2XNSG3fyqHog9vWFWy/RNhg98dj3DcYyCC9hxR3Zfi2CYDb7ChUzZd+XzAibcVlapPnjSX6QbyqfFT7tA3ncDY+C0+tce/XzdFoRwJAW9sEAP4+Hpg3uoPZdBrGRK9l8FRUPjPQPE/CErOZ13DXFRMtHr6eKug0KNUVt9mx+tY0EmNKtQYPLT2C3GLHRqGuOxmHlQdvYc1zPXXrLM3v6A4o+JIIqZRl1joyasrnXBKT2MFfdSZEzZdaw/BnZBLupHMDKpmJsVQaDcPLmyKdnykJq8rloFRr8OZvF9Guvj9i0wtMpqlo6hUr7NZoGGfqmQpUDEjnN4jPKHQ48AKAz/bcAKCt9aouKPgiHMdupiPFQr+evGKVqP0tDEe8VecOm+Jw/u/9x/kkvP/3FeM9m9i1RO49LuvN3y5i/7UU7L9mPs2qQ7cwpF1dzgEQsuZLpWHwMgy+6LJ3C4VuXrtlCfX5IkY+3WO6kysA/BN9X8CcGNtxIZmzTGWwMEx1ugaAd/64hAe5/Pb3O3nbtikspPLU76qKylTYfy3Farro8hGn3A73TsqUCZN+jDBaJ8R8UsR2jh6N7VHJ+P1cIq95cRUUfLkZPprk9lx+wENOnOPUnUyq7RDY8v0x6LL4IE7eTjfq8/X3xXuY+dtFXvdnrjLz9J1MXEzM4abldc/VS2GpfbUO3A73wl2FV+/lIa/EuEnLHV6uXFVS+Q2q0gAx/y/jWu4Kt1Pz8dZWfssXqaDgi9hNzMu9RKlG49o1dMupecJNfWCNu9bE/HjsDtQahpf+Fwm5iRKjYv41vpj7HU/HZeLzvTcq00E7Sre6E2p6G/30pqaCcSY+5xIkzuCcx6DHfzwlemuLs1DwJRFuet/mnYYxjO/SSLe8RO9mTJxLITPV5Z1/N228se+/moLXtlxwcm5IBf1mx61nhW0qUqndq4R0xX+NSq0RfLBTXol9U1e4Egq+eOautR9S4SGXc24CYo+8rE4Ucpkg82zFZdg2bcTRmDQn58Q1OHpI7P2eTKQO96b8fTEZ8ek0vYhQisvUeHjFMTzx02kabc4TGu0oEXzd0xhz74koDYecT+3TDJkFpWZSEz7JZdoATCoKy9z3qVgIqio8uQh9+zUs05bsjRE4B9IkVBy0+/J9JGcXIzm7GA9yS9BIr+sH4N73HGehmi9iNzGffDwMbv6dGwfg1c1RIuWmelHIZVBIqJTNd+MmCXs4ejXuunTP5rSFpSrODfa/mDRsi0xycM/E1Vir6ZROqeA6qOaLuBS5TMa52ZiaD0osUhl55CwKuczkZJd8OnqTmhKFkldse/BqOP1HfokKc7df5jtLZtHNXdpovkX7Uc2XRLjSbXvU1ydE2zdd4+KRy2RObXZMyS3BC7+cd9r2CZc9x7JUpUFEbIYTc2OZK5WPhNiCar6I3eJE7OhKfT3Fo5DLnPp6odtpwk5fUN3ZE3y9/fsl52WEOEyo4lB/Tjh3KoJLlGr4eCpE2TfVfLkZd7owTKGRNuKRy2RObV4Q86XdrszRa0JCYyesosteXD8ci7X4uQudShy3UsV74KPgi7gUtYRLYQlnzS4aDTN6hyYAeCic2+E+3sYpJgg/qJ8OsVVGQZnFz131VBLzgY+aHXnm7p2uxaaQy9wmyJEijYah1ft7AQB73hqIujW9dZ8p5M7t87Vw51WnbZsYc6WaRipXpc1V37Up5iVAwZdUUNlik04Na4mdBbd2U68aftw34ZzP4tILMf7bcMOvEJE5WoOloHYP0fDWfUKEJ9GcojIE+3mJ1leKT2LOW0jBl5vRXtSu+RRii2M309G6rp/Y2XBbUppElTiXr5frFP+n72SKnQVSbtw34Qjy88L5D0boygsXqkTlEHPeQte5+tydi568QguPzUC4iEPe3Z2fNxUJrsbRWpR29WvynBPnmf1ntNhZcHv3coqRnl+KnKIydGpYC/Vq+ZhNm1VYhrxiJQL9vATMIf/E7PdIJS0hPHGHlmMpzWBPCBFGbrESA5b9p1uWy4C4peNEzJEwqNmRECIqxhhm/BqFmBSaa4sQV8HXA1/MgzzOsi2v/dR/TnPVZzYxe1lQ8EUIQUJmEQ5eTxU7G4QQEThSA+SqIxz1iTnil8a7SAVfg1/42QypRjadTsCKgzfFzgYRGE1YTCoo1cbnQov5e/Da5ijzX3L92Mvp76q1uG/R9kwIEV1WYRk+2nUNey4/EDsrhBAB5ZcoUViqfbn6zyfjTKbZdzUF93KK8UjXRha35aoT9lKzoxuhh8nqyxVrEoqVauuJCCGS5GiRU1CqQtfFB+GpkCN60ShcTMw2mzY6KQeeJqIU/fLONUMvanYkhIiERjcSUv2ci8+EhgGlKg1iUvLhaWHGXXPvdHXBZ00jYpZ+VPPlZtzhgiDCkdPjV7VFRQUBgMd+iLB43/DxlJsczegO54+YzaVU9BJSjfFV7X42jmYgJ8RWfAUujr7z8l5OSeU2rGzC3GuEOM2OLlqBLma+KfgipBrjq9lx8tozvGyHEOJ8H+2y/SX2U9aewfaoZKP1+jGbq047IWauKfgihCfuUA1PCHF/fHRP0bhBHxdqdiSEEEKI62D6f7pmIEY1X4Q3rnoREEIc42gFhBtUXFR7Yh5Ddzh9aKoJQkETEQWddYQQR7hF8E4d7t2HO5yPhBBCiCXuUGFAox0JcQNu8SRIXI6rDvMnrk2/vHPVso8mWbUF0wDZd81/7lcH8PLT/q1RA7nGQ2M5/OsDnj7av9VKIO8+msjSzSZPZYFQlv9c3ihDXVmuyXTeBclAdhFQqzGgKP95y4qAQvPbBsA5e31RgiBZvsWk91BXt1wTRQiQFWoXshMBD4OYWu4BBDSuXC7OBkryOEn0/+1lzANpCNQtByIPfrJSs/kpZZ5IR23dcjByUUNWZjZ9EfNGFmrplusiB94ypdn0BcwHOaipW66PLHjKzL8WJ4/5Ig9+uuWGyIRCpjGbPof5oQC+umVL5wEAZLGaKIL23JFBg8Yy7RxXHnlJgKev8Rf86wGeNbR/q1VA3j3Ox4b7S2O1UQZPAIAXlKgny7GYnwcsCGqUz8WjLAYK0iymR0BT3eyqrDTf6r83mVWea/4oQu2Kc81Q9l1Azp0TqBYKUUtWZHbbSqZAKoJ0y7WRD39Zidn0hudaEPLga+HcLGZeyESAbtnauVbIvJGtd27WQza8ZCqz6fNYDeTBv3JF7j1AYz49atQGfCrz0xjpFoMn/XMNYGgiyzBKU09dVFk2+tUFvMrPQRPnmj6vgmx4owyl8AIAeEKF+jLzr5kBgBQWCJUN5WCF+ywYmvJn/BooQbCFcg3gnmt+KEagrMBsWg2T4T7q6JatnWsqJkcKgnXLAShATVmx2fSlzAPptpSDFb+9Zw3ttV6hIE17PZrgkVeEQORxzrW6yIa3hXMtn9VArt651gCZqFlSA8iuYfoLPgHa861CTpL2PgrTZVw280chKrdlrVwIv34XT/XvoF3QqK2mT2cBunPNAyo0sHKu2XrPraBfDvqgFHVkeRbTJ7M6qFExh1lpAVBkYa5CmQyo3axyuSQXKM4xn15uPbSyOfjavXs3fvzxRyQlJaFx48aYMWMGJk6caDb9rl278N577xmtf+aZZ/DRRx/ZuttKpfnA6i7mP39qE9BpgvbvklzLaQFg+m6g5SDt39l3ge96ItzbfPKw0mW4ybQ/fogsAX95f2w64T/l/58dA9RqqP07/jiwdYrF7HjU/0f391jFWazwXGM2bQHzQWjpBt3yZMVRLPTcol343sQXApoB71ypXI74BghfxUmi/2+/ommBR8qW6JbnevyBpz2Oms3PCXVnTFMu0C1/5rkBYxTnzabfpe6Pt5Uzdcvfen2Lh+Q3zKbfqBqFj1XP65b/57UcHeRJZtN/rZqEr1VP6JZ3en+I+hYCmEXK6fifOky3/J/XbHhZCO7eKnsD/2gGAAB8UYpw77e1H2w284Vn/wLaDNf+nXfP6Nw0PO/Gl36Gq6wVAKCdLAm7vReazQsADCxdXXnTunsK2DzJYnosSAa8tcGs150Dlfk3oYwp0K70V93y44qTWOz5P9OJVwPwbwCg8tx60WMfZnn8ZXb7MZqmGF22XLc8y2MHnvc4aDb9GU1HTCn7ULe8yHMTJihOmU2/T90brynf0S2v9PwRgxVXzKb/TTUU76te0S2v9/oSneUJZtN/r3oUX6r0ru0No4HcRLPpMeozoP+busUD3vMsBptzlDOwXT0EAOAJteljlQ/tbw8AU38H2o/R/l2YbrEc7AkgRPYxLrB2AIBWsvs44D3ffN4BDC1diXimLde6y2Pxu9dnFtN3K1mje3AaIr+Mn7y+tpi+Rclvur8fVZzCUs/1ZtNmM390L12rW35WcQjvef5pNn2cpgGGlVWemzM9duIVj71m00dp2uLxssW65Q88f8MTihPGCSt++3ajgaf/qFz/z5vArf0mt90cwAeegzFH+apu3U9eX6On/LbZ/PysGovPVc/qln/z+hytLqUAl8x8YdiHwOA5lctrBmkfvGFc5gDAAuVL2Koerlu2VC4AwKu7Z+Gp/trfR1aaazX9lLKFOKPpBABoKkvHUe93Laa3+Z5brk/J97pKgwHyq1jvtdJi+rYlm+BR8Vql67uAXa+bT+zlD7yv9yBzYRNw0EK5HNAMgIWAAjYGX3v37sWcOXMwffp0DBw4EIcPH8a8efPg4+OD0aNHm/xOTEwMmjdvji+++IKzvk6dOibTV3euWm1LCCGEuBo/L9Mz9wtFxpj12/7IkSMRGhqKr776Srdu1qxZuHnzJvbt22fyOy+++CICAgI433FUr169AKZB5OG/zSfiodlx4Bfma3dsrQJd8WRXPNQyyO5mxye2JiMyMQdA1Zodj8weAm8Hmh31/+3U7OhYs+Pml/qiRbD9zY6G552jzY4Jy8bZ3eyYkZWJiV/+azG5rc2O4e8NBeQKtFgarVtXXZodE5aN066ws9lxwPyNVW52bFXHD5te7KNdsKPZMSoxG09vjadmRzNsbXYMf2+o9g87mh0TMovw2LpLVW52nP5QE7w2pLXpL1hodjR1r7O32TGT1cSNZY8DAO5nFeCpL7dbTC+1Zsdsz/q49kl5LbETmh17Dde2xEVGRppMYrXmKykpCYmJiZg9ezZnfVhYGPbt24ekpCQ0bdrU6HsxMTGYNm2atc3bTiYHApvbllausD0tACg8gcDmnAvfklJ4mU1b6t8ECDT4zMsX8LKcH5m8spAsgg+KmI+F1Fz58EU+0xa4rHYzwMy7uHRqBGr/02Pp356NWsi2o2YuEwF2DftMR2270qciyK70DxBsV3pbzwMAYJDr0qtqNQECa1r+gsLD6Ny0tL8yeNqVH3jWsO/c9/K3a/sF8EUBMxFgAnr7rQy+8uCHPOZnOr0JOaiJHGblN9SThVrIcuK5loZA+4Yw6z/k2OAe6tqxfZnJY+Ur9zd9zE2ca/pKs/xQispyRwkPu84FS+WgKcXwQbId5VohaqCQmenPZIK951ou/JHL/K0nLGe2HDT3G+sHYgZUygJkI46zLt3Ocy0Fwcj1bmT79V678j5ty3GzqxyUK+xKr3LyuVYCb6vp/WV6lRTe/tr/bOUTwHmIcoTV0Y5xcdoTpGXLlpz1zZtrD3h8fLzRd9LS0pCZmYnr169j9OjRCAkJQVhYGHbu3FmlzBJCCOEJdXUgRDRWa77y87XVxP7+3KjQz0/7hFFQYFwtHBMTAwBITk7G3Llz4e3tjZ07d2LevHlQq9V4/PHHq5xxqbKhFZcQyaDTlZDqx8dTjhKl+a4YxPmsBl/Wggm53LjyLDQ0FD/99BN69+6tC9oGDhyIzMxMrF692q2DL0fRTZAQQogQBreti4PXU6u8HZVaUzli0MWUKM33GxaC1V+tZk1t/4vCQm4n24oar4rP9QUFBWHo0KFGtWVDhgxBamoqsrKyHM4wIVJFT5KEECE5Osv8wLb8zDqg0rhurYHYebcafFX09UpM5M5dc/fuXc7n+i5evIht27YZrS8tLYWHh4fJgI3wg2rQxLNw51Wxs0CIzaioqL74uk+oywMY6m5jP6vBV/PmzdGkSRPs38+dLO7gwYNo0aIFGjVqZPSdS5cuYeHChbq+XwCg0Whw4MAB9OjRA56enjxkXZoyC8xPsUDc26WkHLGzYDd3eD8bIUQcagq6HGbTJKtvvPEGFixYgICAADz88MM4cuQI9u3bp5vDKysrC4mJiWjTpg38/f0xadIkbNq0CTNnzsSsWbPg5+eH3377Dbdu3cKWLVuc+g8S29J9N/B4zyZiZ4MQQggxia+aKkY9LRxmU0+5SZMmYfHixQgPD8cbb7yB8+fPY/ny5Rg7diwA4NixY5g8eTKuXbsGAAgICMDmzZvRpUsXLF26FLNmzUJRURE2btyIrl27Ou9fIwEZVPNFCCHEIp5qjESueKKaL8fZ/G7HKVOmYMoU0+8nnDRpEiZN4r5PrnHjxli1apXJ9IQQiaCy0+XR/Y+IRe3CHe7F5ppjRIlZ1IeHEGILCtqqL74OfcX9hs4l+1HwRQghhLg4jRi1UBR0OYyCL0KqMSo7qye1hmHejstiZ4NUUcX1m5RVhD5LDuODv6+Isn9iPwq+CHEjafkleGvrRRyNSRM7K0TC9lx5gHs5xWJng/Bk6b4byCgow5azidYTg79mwortFJap+NlgNULBl0Tw1gZPjyLV2vwdV/BP9H28sPG82FkhAtsRlYyX/xeJrELrI67T80sFyBERSplKnDkfNOU3nG+O3BZl/66Mgi9C3Mj1+3l2padg3X28uy0ah2+k4rPd162mlQmQH+J8js7XxV+He2DzmbvYeyWFpy1WHxR8EeLGcouUSMgotJ6QuI27WUVW08go+nKa1LwSHLqeCpVanNqoglIV9l55gIJS5zcFMsbotWoOsnmeL0KI6+n26UEwBuyfNQgdGtRy6r40GoY0as4SnGEthrw8sGKMgTFALjeOtCj24k9xmRo1vBS65cFfHEWpSoMPxnbEK4Nb6dan55cioIYnvDz4rfPIL1HiVmqBbvm59WdxMTEHABC3ZKzu+J+Lz8LdzEI8weMbWKjm3HFU80WIm7iYmG1Uo1FROO68eN/kd/icF+61LVF4aOkR3rZHHCOTycAYw+M/nkLY1ydQqlKbTEP4sXTfDc5yaXn/q63nKju/30zJR+/PD2P8tycB8Bu0DF1xHIl6tZ0VgRcA9Pr8MO5mamu+n1pzGnO3X8aRG2n41IamaeJcFHy5GXoQqb4e++EUZ1l/3h8h7rUHrqU6fyfEKhmAhMwiXEjMwe20Ahy6bnxcKPbiz/aoZKtpwr4+AQCcGio+MAZkFJivbc4qLMMn/3IDrX+iTT+IObp/4hgKvggREGMMJUo13tx6EZtOJ/C+/Qe5Jbq/9d+7di+bphVwpjvp/N5Uq0Iuk3GatopKTdR8CZkhNycD8NWhW5jwfQSy9UaaZhWV4WZKvngZK1em1iAtv7JcOHk7nbdtn43P5G1b1Q0FXwRn4+gC4ktsWgH+ib7Paeqp+Pv6/Ty0XLAXHT7cj3+j7+OjXdecmhf99679E30f93OKEZuWj7Un7iC/RAmAnlz5EpcunUENcjngqdfPq8jUHExU9cWbwjI1Vh+5jeikHKw8dFO3PqdIibCvTyA6KYeTvqhMBaVauAvv5O0M9Pm8sjtAdpGSt23P3U4T9TqKOtwTTF57RuwsuI0Rq44DAF4a2BLzRndAz08PIV+AUUemqAxeN9J/2X+6v2+lFuDNYW3w8v8ihc6WW2KM4U56ARIyCjGsQz1dn6r8EiWik3LxUKsgeCiEedaVQcbpfnAz1bj2hUIv59h8xniS078ucJsll+2LwabTd3nZHz07uS4KvtyMo/O+EH6tD49Hen6paIEXAPT67JDZz7ZHJePGgzzcTpNOc5krYwCGr9QG3muf64lRIQ0AAFN/PoOr9/Iwc2gbzAlrL0heZDJujebWc0lYOqkLYtPyEZtWgLCQBlTxJaD/GQRafAVegLacIa6Jmh0lgoIm98Nnx1ZHlCgtzzN0zc4JWYl5+pevfgf3q/e0v/F3R2MFy4tcJtPNPK5vxKoTeHXzBey9kgIZ1X0RIioKvpwkPb8Ub/9+EYdNjDQypNYwXNAbHuwMFX18CCH8e3VzlO7virBHrMlt5TLLzVEHr6fAxNRfvBrQJhit6vo5dyeEuDBqdnSCNcfv4EJiNg5cS8WuS/eRsGwc1BrtKDc/bw+k5ZXgw11XMbFbY4zp3BBfHbrl1Pwcvp6KlzdFYmK3RmgQUAP9WwfjxK10PN23GVrV9Xfqvollw1ceEzsLhGcVlU4PrzhmMV2K3sjUqu2PcWrOLyTmWKxJV2uYIM2OcmrbJMQsqvlygqX7YozmPJr68xl0//QQ7ucUY/5fV3DgWipe23IBAP9NEpEJWXj794tIzNROvPfyJm2n6p2X7uOn43cwbcM5rAuPx4TvInjdL7HfHQmNkiP82HEh2eSrZf6Nvo+526JRUKrCv9H3eZ2QVn9sRW6x0mgUq/5UGKaaJPnGGFCiNJ7ighCiRTVfAjkXnwWAO+LMGRiAJ346DQC4ei8XR9592GxaMTuDE+LOvj96x2jdm1svAgBq1fDErzx2us4tVnKmFQGM5x2rGAwAAHuvpCAtj9/XQL03uj2+2F85zULflsE4dYemsCHEHKr5EsB/Meb7fd02MQycL7bUqkTEZjht/2KJXDhC7CyQai7ybpbZz9aHx6OMx5cuZxSU4XJyDmfdAytNmpF3s3nbPwC8/nAbeOtN7DqiUz1et0+Iu6HgSwAvbjQ/l9LIr07wuq9By49ylk9ZCa6eWXeW1/2LbcbgVgjy9RI7G6SaK7Uy0pRvFbXdFcQYPK0fUAo9mvKXF3pzlsd1aSjo/gmxFwVfbia3mDuq8WkRgqsPx3fiLI/oWF+wfc8a0Y7mMCKiO2+h5ksIQvTrMtSwlo/ub0+FsBfh0PbcmraF4zoKun9C7EXBF+Ho1rR2lbcxrV9zzvLqKd2qvE1b1fBSQCaTIdDXEwDw87ReNn83Ydk4zvKO1/ohLES4wJGIZ+9bg3D+gxG49NFIXrbHR+xTlekgzsYLH/ztnDlA93ebeuKOom4YUEPU/RNiDQVfelY82RUA8PrDrVHTWxpjEXo2D8TJ94YKtr/a5UELnxTOnlTIhFPzh+PcB8MxspNtwdPLA1sCAL6d2h1ymfbJuWfzIKx5rjJ4e39sB0zp3RSeCpnT/02HZw/G3rcGOXUfhpoH+wq6PylpFuyLujW9UdvXC0/1amL0+bkPhguanxcGtICmCgHcvyJM8Fuvpg8Slo1DwrJxutcruYJuTWtjx2v98O7IdmJnRXSeChmuLg7Djtf6i50Vp2tXX9wHBGlEGBLxRM8mGNahHoL8vDA3rD1WHLyJn0/Go0xV2Zehce0auJdTrFt+7eHW8PVUYKWT5ura8Vp/pObxMx+QLY7drHzjfePaNeDv7cF5N9y4zg2RXVRmdiTTwnEd4an3DrsGtXw4y0Kp4aVADS+FzenHd20EAHikayMM71gPvl6Vl8YvL/TGhbvZeL5/S3gqZFg4vhN8PRX49r9Y3EkvMDmT/ZTeTfH7+aQq/AtkkMmEbTr6+/UB8PaQ4+N/riEtvxQdGtTEmhNxnDRvDG2N2LQCo6lUTPlgbEd8vveGs7Jr0YiO9eDv7YGdl+wPQr54oise7doYz67XNtn/+lIf1Kvpg0Ft6+DkbWEGqCx6JAT3c4p1v7NhuUOM1fb1RI4DL43e+Ya2xq6Wj6fJcvzTiaHYHpmE6OTcKudRiiIXjoBCJsOd9AL0bB4ImUyGns0D8UjXRg4F8cF+XsgsLHNCTvnVLEjch02q+TIQ5KftrC2TyTA3rANufDKa87nc4BebN7oD3hze1il5qeXjwcmTOWNCG/C2z25Na6NVXT8M71AP4fOG4sA7gzmff/5YKHq3CNIt31kyFhO6NdItP2zQ98LXS8FbLVFXHppEzdHPoX7gBWj7k7w7qj28POSQyWTw9/aAXC7D2yPa4gMzfUuWTuqMiXq/i0N5MvjZejUP5CyPsrFWz1YKuQx+3h748smu+N+LfbBgbEfc/nwMzr4/HH1aBOGTCSGYG9YBq6d0t2l7rwxuxVmuqF0UQkANL/RvXcfm9IZnqH6fqRqe2iD+0a6Vx/P43IfhxeNDxZvD2hitWzapC2aNaIuD7wxGxPxhvO3LlD1vDcSiRzqZ/Izv88xZqvqQ17Z+Tfz4TA/4eFZu563hbfFs32ZYPaW7rjx+pGvVrmupqePvjUA/L/RqEcSpsfxmSjdELxqlWx7WwfQI1q8nd8PGF3qjprcH6tX0xpn3hyOkUS2n59vVUfBlhWHg8OIA4W4gnZsEANAWKpELR2BYh3poYdA01Lh2DUzu3ZS3fdb08cCR2UOw/vneugvx26mVN1sPhRz9WwfrluUy4MiNNN2yudaGioLLHjte68dZ3jajn1Eaw35aptjS56ymA/kDjPv2zBjcCkfeHQKZTIavp3RH1MIRaBKo7X9irRMyt78dg8Lgx5QbnIsdGtR0KM/mmAqSPRVy1K/lgz9f7Ydp/VoAAHw8FXhhQAuL2zo8Wxu092hWW7cuzOAh4c6SsWa/v/etQbpuAI5oW8UmBf3DWnEYJnZvjAndGmHe6A5oHuyHOWH8NVP1aRmE9vW5xzPQzwuzRrRDu/r8HmdTQhoF4Pn+LTjrvniiC94Z0Q6rJndz+v6t8bWhFpuPfnZjOjfE4LZ1dcuzR7aDTCZDizp+uPxxGBKWjeOUh9YY/qauRCaTIaCGJ9ZN64XxXRpi2eOdMTrE+EF/fJeGeLh9PVz+eBTOvj8cngo5Vk/pDl8vBSb1aGxz1w8A6Nw4oMr5XvtcTxv7GYvbNE7Bl52m9WthFADxqWUd0+9Dq+PvjQ3P98axudz+X4F+/PbRkslkRv01xnZuiCd7NsGcUe3g7+2Bvq2CseH5Xtj39iDIZDIU6E3Wam6I/cF3hljd96tDWnOWezYP4ix7ecix+NEQW/8pOhO6NcZbFmonGwX4OPyaJaZ3m54/pgMWjO2I1nrbCvb3xp63BmHtcz1x8aNRGNpeW7AH+Xlh/pgOaFPPHx+M7Yjp/ZpzRmgF+XlztgOAE4w1CvDBDIPfCwB+/7+HHPp3AICHHTWU74V14OTXsP9Em3ragGHHa/3x3dPdsXpKN06NKaAN3F8yUxvWqVEtPNHTuO+VrWYMbsU5NtYY1pro13xVXA8VN5XXHtb+7hXBKB/6tAzCgXcG45fne+sCV0sufzzKahp7yWQyvD+2AwDtCOWnejXF2yPawt+g/2uXJlW/QepbN60XJvVobPbzbk1r4/BsbvmhH6jWrekNwPgBb2of6w+lFd/VZ0tfu9DGttXsBNow7Y09D6ZjQhsgfqn5hxZnGNGpPr57ugfq1fTByqe64tOJoZzPPcqvHf17R5t6/rjw4UisfLIrehrU2OvrY1AmWHvQ6N3C/LYqjAppgAndzJ9PFcTulkjBlw3WPNcT9Wp646dne0Ihl+HAO4Px54x+uP35GN731dqGl9HufWsQ6tb0hq+XAssf72J0k+abQi7Dl092xcxhlQHMsA710bGhcQHUONBglFH5Cd4gwAfjOluee0cu0z5FAeYLzmf6NrMj55Vmj2yHhGXjjGoXEpaNw6kFjnem1n/aDjQzWCGghidGhTSAv7cHfnmhD+KXjsWZBcPx6pDWODx7CF4Z3AqLJ4SiV4sgfDO1O9ZP74UgPy/I5TLdjR7gFhanFgyHn8FN8eR7Q/FQq2A4yp538dXwUuDlQa1wav4wvDuyHTa92BdvDNXm9adne+jlWYbxXRrpCsO5Ye0BQNe3xHBaEkOvDKoMzurW9MaSxzrblD+ZTGbzXFMNavnAy4NbFLbSewhqXNv0yDkfTwVvI3k9y/szDO1QTxe4GtK/3mr58D8wBgBeGdQK+2cNwo96xxAAuuoFXP97oQ++mmy9VtLcyOklj3XGhucrB7L4+3hg1VPd8NvLfTGld1PEfMrt6vHZxFA0ql0DG57vhffHdkD80rHY93blYJTg8m4Z3z3dAw1q+WDBGG0AaercerJnE+x4rR8e6doI307tjv/eNX4oVGmsz9G2+NFQs5810J9yw0OG8x9YnvT57zcG4NmHmlnsnnHug+G4s2Qsfny2p6iDGfy8PTC5V2XZbGk0uY+nduT5SwNb4s1hbfDLC70RMX8YNr/UF4C2GfPPV/shVu8+WtPHAy8NbIm+LYOMtjeuc0Nse7U/dr85kJd/y+sPGz+8Cok63NsgLKQBwvSqW709FOhjcHK0quuHOB7e0/faw21wWK8Zz5ROjWrh/AcjUKbS6G4aQ9vXxVG9zvKOamIYPNkpoIb5m8L3z/SA+tco7L+Wgu7NaqNfq2D8cKzyNSwyGbDqqW6Y3r+FruBe9VRXzP4zGu+M0DbxeJjo17F/1iC8/9cV/J9BHyNTxnRuwBlAUFWcGhIbb/YymQxeHqbTPmrQn0R/+z6exk0vC8d1xGd7tJ3aTT3F28ORvnmNatfQ9XmcG9YBM4a0thgYvDG0DV4e1NJsf6lAX08sGFtZozZvdAf8fDIegDY4n9SjMf44n4gWdfywekp3tFu4D2UqDRoG+GBavxZYvj9GV4syvmtDfH34FurW8sHCcR3xpMFEpIse6YQWwX4mA9bmwX5YN60XPBQy1Ne7mRqa0K0x3v79ktnPbWXL/XTtcz2x+sht3QPIsw81w+YziVXeNzcfMnRoYPxQ9UTPJroO54F+XnisexM81CoYc7ddNtlfDeAGIfqm9G6K8wmVU2Foyqua+repg/5tuP30wkLqI7S8KWpYh/oY1qF+eT4r07QI1gbKnRrVwukFw3TBia+XB94a1gbf/Ff57twlkzrDUyE3qlXXp7Th7QM9mwfio/Gd4Okhx2e7r6NUb1DWaw+3xqJ/rmnzCRnq1vSGt4eck0Zf67r++GxiZzzUKhgzf7to9PmCMR1Qr6b5c9Cark1rIzopx+HvG/LykGP/rEHILlSiX2vrD3ueCjneHdVet9y4dg3Efj5GV5YblukVQXOL+Xt064L9vPDd09rm3lAzTZMPtQqyGBTrOzV/GBqZeagSCgVfPNn95kDEZxSiuExtNNu0vtufj8GpO5mYvuGcyc/1q2it9WHQf1r/Zmp3DPnyGLKqOMrkvbD21hMZ2PvWIIz95qTZDrv6fnimBwrKVKjl44nk7CJO8NW7RRC8POSc5qlJPZpgeIf6CNCrVXqoVRDOxGVh1gjtTb9Dg1r46/XKOYYs6aVX6Jpr4rWHn17nfKNaPx5o9NpAFj3SCdFJObraQQB4YUBL3EkvRMs6viaDM3NaBPsis7AM+SXaJuMRHevxMjDClhoZbw/z+bzw4UjOk71+wdyodg34eCqwa2blk++tz7i1z0/1aoJgf20Q6uvlgePvDYVCJoNcLsOuNwbgyTWndaOXX7DSf3OEjX1VJnZr5NCoygon3xtqU21G0yBfTj+4zyZ25j34Mufpvs1Rv5YPQvRufA0DamDzy33NfqdBgOmAQS6Xcfovqi0UdBXH0pTlj3fGv9EP8MmEyq4Ihr/j7FHt8dKgVjgVm4HuzQJt6pT/3EMtEBGbadTcaujF8ibz7/+LRUr5iPRnH2qGx3s20QVfFTa92Acv/y8SM4e1wdJ9Mbr1+sdzXOeGaD7TDwE1PHE6LgPZRUr0bhFoMlC88clodPxoPwDrD/4/PtMDz64/y0vlQAVTAbo9TD1EA6YfQoa2r4tfXuhjMv30fs3xUKtg9GwRaBSgymSm76FvDmsjeuAFUPDFG18vD4Q0CkCpSs1Zv2BMB93Ftua5nvBUyDGkXV1djU6v5oG8vGetpo8nIuYN012QAHB49hCMWHXcwreM1Xbg1TydGtUy2/HdsIZDLpfpbtBNAn2x562BuH4/T/e7mBJg0Jy38YU+uHY/F92aWm//NzSwbR2sntINZSoNL6OWAv288PljoUjPL+UMROCLfiHVPNgP5z8YwblxKeQyLJ1kW1Ocvv/KX7i+dN8NtKzjj6cdbM7lw0sDW2J9eDyefaiZySBk+eOdsfVcEr62oeO34c1a/2bbtWltHJvzMN747QKm8DhIZcmkzhaDr/ljOuDX03fNThXR1MlD3n97uW+V33ShkMswykRna0tmjWiLqLvZ6Nc6GCdupSMmpbLG2dLoYn2WQtLJvZthcm/r521ADU+MsdLlQV9YSH38/Xp/XY2aNeum98KMX6Pwf4NbYXr/Figqq+wDW3E6920VjEuLRkEhl6Fjw1pY8NcVzB/TgVMGyWQy3SCrZsGW/101vBRIWDYOjDEwBlxKzsGkH06ZTCt236aqMhWaT+3TFKfvZGLWiHYINDMbgEImg8og+hrUto5R32KxUPDFM28PBYa0q4vjt7RNgK8MaoUezQPRoJYPp5Cd1KMJJvXQdijWr16tCsN5reydZXquA7Ve5swY3Aq/nU3Eyqcs9w0JaRSAkEb2deD18VRYbDawxpbOmPZ4pm9z64kc9H+DWmHXxXvoWV4baDji0VEV2/lgnPXaSmf7YGxHPNGziVF/vAq23mRt0ah2DfxtYy2prXy9PHQBpCmvDmmN7k1rY/LaMwCAqIUjEFDDE39EJqFrk9q85sWUTiaG/fdsHoio8oe+XW/w+3tUqO3rhX/L++fEpRdwgq/uzQLRr1UwvDzk6O7EKWQcIZPJ0L2Z7Q92oY0DOFOBeOjNR6Rfu15Rszy4XV3epg7RdnI3fsjlpIEMyyZ1wVNrtC0yK5/sCrkc+HDnNSx/vAsv+eCLqa4bPU0ci6WTuoAxZrHGWNuHtTL48lTI8OtL5mtqhUbBlxN8M6U71kfEY0xoA8jlMqNRXtZUBG/zyzuO2mPjC73x0/E7WPSItir+6uIwhC46AED7BNy/TR0MWPafyafwHnYUONYsGNsRc8Pam61eJrYJ9PNC+LxhdgVdW195CCsP3jRbo2rYuVxs8vLaAFe2YEwHHLuZhjvphejfOthoEuK+rYLx2cRQ1PH30tXO8R20vze6PZoE+uKtrdx+QxqmHVH48qZI3bo/Z/RDTlEZCkvVaCbAmw0+f6wzipWXdA8+CrkMWy2MzB3Utg6i7mbjjaGm+5NJmZeHHF9N7oq7mUWizJG2dFJnBPp64tXNFxDo64k6/l5oEOCD6I9GoVYND13AMqFrY94e5qrKz0uBwjI1nupdOcL5t1f64lx8Fv5viOm+vNaa6r+Z2h2vbo7CuC4N8VDLIKM5KMUmY0yEN7DaqVcv7YiKyMhIKyldU/jtDKw6dBMfjOuEns0DodEw5BYrzVan2qtUpUZGQZlu1FZ2YRku38uFRsPw+/lE3Szav//fQ1UaLUekxVyN6ofjO5md4oE4rlSlxuk7mejdIgjRSTl4et1Z/PBMD4y1o8nLXhXHuHVdPxwpb0pmjKHlgr26NLc+GwMvDznnfIhfOlbSrwDSaBhKVGqLTZKk0pXkXDzyXTgA7QjEkZ3qI7dICS8PuV1v+hBLQakK2YVlvDfBZxSUItjPS5Rz3VrcQme2BAxsWwcD21aO8pHLZbwFXoC2KVR/uHygn5euf5WXh1wXfNkz1QCRPv3Xg2x4vhdKlBpcSsoxevE54Ye3h0L3dN2/TR2bJgDmi/7NxfBGU1HTufmlvrpXJkk58AK0ZSAFXrbTf1VOxezyhn1lpczf28PqAAdH1LEwYENsdHZXc9zJJEXMCOHdN1O64bOJofBSVD79OrMWhghvXJeG2HP5AT6faH2I/cC2dXB1cRh8JNbsTKouwNcT+94eBKVaI4mRfMQ6ugqrueZB1ieTJK6p4vUgrtDsQBzz3dTuuPjhSPQ16C5QMVeS4Yhef28P6ofppjo2rIUuAgzgIPygmq9qrlmwr27mfnpiIsS1yGSmuyi8NLAlpvRuatM7EQkhwqPgi2B0qH3z9xBCpM/w9VOEEOmg+mdCCCGEEAFR8EUIIYQQIiAKvgghhBBCBETBFyGEEEKIgCj4IoQQQggREAVfhBBCCCECcomxyAUFBWCM6d6VRAghhBAiVfn5+RZf4+USNV9yuVzy7yIjhBBCCAG0EyDL5eZDLBljei/3I4QQQgghTuUSNV+EEEIIIe6Cgi9CCCGEEAFR8EUIIYQQIiAKvgghhBBCBETBFyGEEEKIgCj4IoQQQggREAVfhBBCCCECouCLEEIIIURAFHwRQgghhAiIgi9CCCGEEAFR8EUIIYQQIiAKvsrt3r0b48aNQ5cuXTBmzBjs3LlT7CwRKzQaDbZu3YpHHnkE3bt3x4gRI7B06VIUFBTo0ly5cgXPPfccunfvjoEDB2LVqlVQKpWc7SQkJODVV19Fr1690LdvXyxatIizDQDIyMjAu+++i759+6Jnz56YPXs20tPTBfl3EtNmzpyJkSNHctaFh4fj8ccfR9euXTFs2DBs2LDB6Ht8nRPE+c6fP4+pU6eia9euGDhwID799FMUFhbqPqfj7V62bt2KMWPGoFu3bnjkkUfwzz//cD53q+PNCNuzZw9r3749W7JkCTtx4gT76KOPWLt27di+ffvEzhqxYM2aNaxjx45sxYoVLCIigm3evJn16dOHvfjii4wxxhISEliPHj3YSy+9xI4dO8bWr1/PQkND2eLFi3XbyMnJYYMHD2aPP/44O3z4MPvjjz9Yr1692P/93//p0iiVSvboo4+yESNGsL1797J//vmHDRkyhD366KNMqVQK/u8mjO3cuZO1a9eOjRgxQrcuKiqKhYSEsDlz5rDjx4+zVatWsfbt27N169bp0vB1ThDnu3jxIgsJCWFvvvkmi4iIYL/99hvr1asXmzVrFmOMjre7+f3331m7du3YsmXLWEREBFuyZAlr164d27t3L2PM/Y43BV+MsREjRugu6Apvv/02Gz16tEg5ItZoNBrWu3dv9vHHH3PW79mzh7Vr145dv36dvf/++2zIkCGstLRU9/mWLVtYx44dWUpKCmOMse+//55169aNZWVl6dIcO3aMtWvXjl26dIkxVnmjj42N1aW5ffs2a9++PduzZ48z/5nEhJSUFNa7d282ePBgTvA1ffp09uSTT3LSfvHFF6xXr166c4Cvc4I43zPPPMOeeeYZptFodOs2b97Mhg8fzoqKiuh4u5nJkyez5557jrPu6aefZs8++yxjzP2u72rf7JiUlITExESMGjWKsz4sLAxxcXFISkoSKWfEksLCQjz66KMYP348Z32rVq0AAImJiYiIiMDQoUPh5eWl+3z06NFQq9UIDw8HAERERKB3794IDAzUpRk4cCD8/Pxw/PhxXZo2bdqgdevWujQVyxVpiHAWLlyIAQMGoF+/frp1paWliIyMNHkd5+Xl4cKFCwDA2zlBnCsrKwuRkZGYOnUqZDKZbv0zzzyDw4cPQy6X0/F2M6WlpfDz8+Osq127NnJyctzy+q72wVdcXBwAoGXLlpz1zZs3BwDEx8cLnidinb+/PxYuXIiePXty1h8+fBgA0Lp1azx48MDouAYFBcHf3193XOPi4ozSKBQKNGnSxGIaAGjWrBmdHwLbtm0brl27hg8//JCzPikpCUql0uJ1XFxczNs5QZzr1q1bYIwhICAAs2bNQrdu3dCzZ08sWrQIJSUldLzd0LRp03Dy5Ens27cPBQUF2L9/P44dO4YJEya45fH2EGxPEpWfnw9AezPXVxGBU6dL1xEdHY21a9dixIgRqFWrFgDj4wpoj23Fcc3Pz7cpTZs2bUymuXv3Lp//BGLBvXv3sHTpUixduhRBQUGcz2y5js2lqUhnzzlBnCsrKwsAMH/+fIwcORI//vgjbt68ia+//hqlpaWYPHkyADre7mTcuHE4c+YMZs2apVv32GOP4eWXX8bFixcBuNfxrvbBF2PM4udyebWvHHQJUVFRePXVV9GkSRN89tlnKCsrs5jeluPKVxpSdYwxvP/++xgyZAjCwsJMfm6JXC7n5Vqn4y2MitFpPXr0wKJFiwAA/fr1A2MMy5cvx1NPPWXx+3S8Xc9rr72GixcvYsGCBejUqROio6Pxww8/wN/fH2PHjrX4XVc83tX+zKpZsyYAcIYvA5U1XhWfE+nau3cvXnjhBTRs2BAbN25EYGCg7snG8LgC2mNbcVz9/f3NpqnYhi1piHNt2bIFN2/exPvvvw+VSgWVSqUrbFUqlU3XMZ/nBHGuihqNwYMHc9YPHDgQjDFcuXIFAB1vd3HhwgWEh4dj4cKFeP7559GnTx+88sormD9/Pn799Vf4+voCcK/jXe2Dr4q238TERM76iuYkU319iHT88ssvmD17Nrp164YtW7agXr16ALSFd/369Y2aBTMzM1FYWKg7ri1btjRKo1arkZycbDENoD1n6PwQxoEDB5CdnY2BAwciJCQEISEh2LlzJxITExESEoLIyEgoFAqj67hiuWXLlryeE8S5WrRoAQBGNdgVNWJNmjSh4+1G7t+/D0Bb06mvV69eAIAbN2643fGu9sFX8+bN0aRJE+zfv5+z/uDBg2jRogUaNWokUs6INdu2bcOyZcswZswYrFu3zqiWcsCAATh69CinAD9w4AAUCgX69OmjS3P27Fnk5OTo0oSHh6OoqAj9+/cHoH3avn37tm5wBgDExsbizp07ujTEuRYvXozt27dz/hs6dCgaNGiA7du3Y/To0ejVqxcOHjzIaX44cOAAatasidDQUAD8nRPEuVq3bo3GjRtj7969nPVHjx6Fh4cHunfvTsfbjVQEPVFRUZz1ly5dAqAdxe52x1uwSS0kbMeOHaxdu3Zs8eLF7Pjx47pJVmkOJ+nKyMhgXbt2ZUOHDmXnz59nFy9e5PyXmZnJYmNjWefOndn06dPZf//9xzZs2MBCQ0PZokWLdNvJzMxkffv2ZRMmTGAHDx5kf/75J+vduzd7+eWXdWlKS0tZWFgYGzp0KPv333/Zv//+yx5++GE2fvx4mmRVRPPmzePM83Xq1CnWvn179vbbb7Njx46xr776irVv356tXbtWl4avc4I4X8Xk1++++y6LiIhga9asYSEhIWzJkiWMMTre7ub1119n3bt3Z7/88gs7c+YMW7duHevevbvuOLjb8abgq9zWrVvZyJEjWWhoKBszZgz7+++/xc4SseDvv/9m7dq1M/vfzp07GWOMnT9/nj355JMsNDSUDRo0iK1cuZKVlZVxtnXz5k02ffp01qVLF9avXz/24Ycfsvz8fE6a+/fvszfeeIN169aN9enTh73zzjssNTVVsH8vMWYYfDHG2MGDB9n48eNZSEgIGzZsGFu/fr3R9/g6J4jzHTp0iE2cOJGFhoayIUOGsB9++IGp1Wrd53S83UdpaSlbuXIlGzJkCAsNDWVhYWHs22+/5UyY6k7HW8aYlSEChBBCCCGEN9W+zxchhBBCiJAo+CKEEEIIERAFX4QQQgghAqLgixBCCCFEQBR8EUIIIYQIiIIvQgghhBABUfBFCCGEECIgCr4IIYQQQgREwRchhBBCiID+H+PN8jZm2XXPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,4))\n",
    "ax.plot(landscape)\n",
    "ax.axhline(np.mean(landscape), color='C1', linestyle='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create WHT matrix \n",
    "H = walsh_hadamard_matrix(L=L)\n",
    "# get WH coefficients\n",
    "spectrum = np.dot(H, landscape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'WHT')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAEXCAYAAAD/WnOEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9nElEQVR4nO3deXxM5/4H8M9kkkmYDJIiaCyhTRBEmCSWNEotoS1t/aq9XBellquttaWtVt0qxbVcWhclXKXa0lZra9AGTVTIUtRWsgdBEpFNklme3x+RI2NmkpD1JJ/36+XVzjnfOfOcOSczn3mesyiEEAJEREREVKPZVHcDiIiIiKh0DG1EREREMsDQRkRERCQDDG1EREREMsDQRkRERCQDDG1EREREMsDQRkR1wvTp0+Hh4YHbt2+bzZszZw48PDwwZcoUs3k5OTno2LEjZs6ciTVr1sDDwwPh4eEWXyM5ORkeHh6YO3cuAEj1pf0bPXp0xa4sEdVKttXdACKiquDj44MDBw7gzJkz6NOnj8m88PBw2NnZ4dSpUzAYDFAqldK8M2fOwGAwwM/PDzdv3nyo1xwwYABatWolPY6NjcW6deswYMAADBgwQJreuHHjR1wrIqpLGNqIqE7w8fEBALPQFh8fj+vXr+PFF1/EDz/8gHPnzqFLly7S/KioKACAr68v9u7d+1Cv2b59e7Rv3156HB4ejnXr1sHDwwPDhg0rz+oQUR3E4VEiqhOefPJJNGrUCKdPnzaZfuLECdjY2GDy5MlQKBT4/fffTeZHR0ejadOmcHNzq8rmEhGZYWgjojpBoVBAq9Xi7NmzJtNPnDiB9u3bo02bNvDw8MCJEyekeUIInDlzBr6+vlXdXCIiMwxtRFRn+Pj4ICMjA/Hx8QAKQ9nJkyfh5+cHAOjRoweioqJQUFAAAIiJicGdO3fQo0cPk+VkZWUhPT3d7F9mZmaVrg8R1S08po2I6oyiHrPTp0+jTZs2+Ouvv5CWliaFMj8/P2zZsgVRUVHo0aMHoqOjTZ5XZOrUqVXbcCIiMLQRUR3Svn17aDQanDlzBsOGDcOJEyegVCqh1WoBFPbEKZVKnDx5Uup1a9asGVq3bm2ynDlz5picYFAkNTUVb7/9dpWsCxHVPQxtRFRn2NjYoHv37tLJCCdOnECnTp3g6OgIANBoNOjQoQMiIyMBAH/88Yc0dFqcp6enxenJycmV2Hoiqut4TBsR1SlarRYXL15Efn4+IiIizI5X8/Pzw5kzZ5Ceno64uDiehEBENQZDGxHVKb6+vtDpdDhw4AAyMzPNesx69OiB3Nxc7Nu3D0IIiz1qRETVgaGNiOoUT09P1K9fH19//TXs7OzQvXt3k/ndu3eHra0tfvjhBzz++ONo2bJlNbWUiMgUQxsR1Sm2trbw9vZGdHQ0vLy84ODgYDJfrVajc+fOOHfuHIdGiahGYWgjojqn6JZWDx7PVqRoSJShjYhqEoUQQlR3I4iIiIioZOxpIyIiIpIBhjYiIiIiGWBoIyIiIpIBhjYiIiIiGajVt7Hq2LEjjEajdIsaIiIiopoqOzsbNjY2OH/+vMX5tbqnzWg0gifHEhERkRwIIWA0Gq3Or9U9bUU9bBEREdXcEiIiIqKSabXaEufX6p42IiIiotqCoY2IiIhIBhjaiIiIiGSAoY2IiIhIBhjaiIiIiGSAoY2IiIhIBhjayslgFAi9nIr0nILqbgoRERHVYgxt5XTwXAr+vikcb+88Xd1NISIiolqMoa2crt3JM/kvERERUWVgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhl4pND2xhtvYMCAASbTQkNDMXz4cHh5eaFfv34ICgoye97Zs2cxevRoeHt7w9/fHytWrIBOpzOpiY+Px+TJk6HVauHn54f58+cjOzv7UZpZpYQQ1d0EIiIiqsUeOrT9+OOPOHTokMm0qKgoTJ48GW3btsWaNWvw/PPPY+nSpdi0aZNUk5CQgLFjx8Le3h6rVq3Ca6+9hs2bN2Px4sVSzZ07dzBmzBikpqZiyZIlmDVrFvbv349Zs2aVYxUrl6K6G0BERER1gu3DFN+4cQOffPIJmjVrZjJ99erV6NixI5YtWwYACAgIgF6vx7p16zB69GioVCps2LABGo0Ga9euhUqlQp8+feDg4ICFCxdi0qRJcHFxwfbt25GZmYndu3fDyckJAODi4oKJEyfi9OnT8PLyqqDVJiIiIpKXh+ppmzdvHnr37o2ePXtK0/Lz8xEREYGBAwea1A4aNAiZmZmIiooCAISFhaFv375QqVRSTWBgIAwGA0JDQ6UaHx8fKbABgL+/P9RqNY4ePfrwa0dERERUS5Q5tO3cuRPnzp3DBx98YDI9KSkJOp0Obm5uJtNbt24NAIiLi8Pdu3dx/fp1sxpnZ2c4OjoiLi4OABAbG2tWo1Qq4erqKtUQERER1UVlGh69evUqFi9ejMWLF8PZ2dlkXlZWFgDA0dHRZLparQYAZGdnW60pqis60SArK6vUGiIiIqK6qNSeNiEE3nvvPfTp0weDBg2yOL/EF7CxKVNNacpSQ0RERFRbldrTtn37dly6dAl79uyBXq8HcD+o6fV6aDQaAEBOTo7J84p6xjQajdR79mBNUV3RMhwdHa3WtGjRoswrRURERFTblBragoODcfv2bfj7+5vN8/T0xEcffQSlUonExESTeUWP3dzcoFar4eLigoSEBJOatLQ05OTkSMexubm5mdUYDAYkJydb7OUjIiIiqitKHXNcsGABdu3aZfKvb9++aNasGXbt2oXAwEBotVocPHjQZBg0ODgYGo0GnTp1AgD07t0bISEhKCgoMKlRKpXw9fWVasLDw5GRkSHVhIaGIjc3F7169aqodSYiIiKSnVJ72tq2bWs2rVGjRlCpVOjcuTMAYMqUKRg3bhxmzJiBF198EdHR0di0aRNmzZqFevXqAQAmTJiAffv2YeLEiRgzZgzi4+OxYsUKjBgxQhr6HDlyJLZt24axY8di6tSpyMjIwLJlyxAQEIBu3bpV5HoTERERyUqFHN3fs2dPrFmzBjExMZg6dSr27NmDd955B6+//rpU065dOwQFBSE3NxdvvfUWNm/ejHHjxuH999+XapydnbF161Y0atQIs2fPxsqVKxEYGIiVK1dWRDOJiIiIZEshavFNM7VaLQAgIiKi0l4jKDQO/9p7Hu2bafDz9IBKex0iIiKq3UrLLbyORjkpePNRIiIiqgIMbUREREQywNBGREREJAMMbUREREQywNBGREREJAMMbUREREQywNBGREREJAMMbUREREQywNBGREREJAMMbUREREQywNBWQWrvzcCIiIioJmBoKyfexYqIiIiqAkMbERERkQwwtBERERHJAEMbERERkQwwtBERERHJAEMbERERkQwwtBERERHJAEMbERERkQwwtBERERHJAEMbERERkQwwtBERERHJAENbBRHgzUeJiIio8jC0lZNCwbuPEhERUeVjaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhkoU2gTQmDLli0YNGgQunTpgqFDh2LPnj0mNaGhoRg+fDi8vLzQr18/BAUFmS3n7NmzGD16NLy9veHv748VK1ZAp9OZ1MTHx2Py5MnQarXw8/PD/PnzkZ2dXY5VJCIiIpI/27IUrV+/HqtXr8abb76Jrl274tixY5g9ezaUSiWGDBmCqKgoTJ48GYMHD8a0adMQGRmJpUuXQgiB8ePHAwASEhIwduxYeHt7Y9WqVYiJicHKlSuRnZ2NDz/8EABw584djBkzBk2aNMGSJUuQlpaGZcuWISUlBevXr6+8d4GIiIiohis1tOl0OgQFBeFvf/sbpkyZAgDo2bMn/vzzT2zbtg1DhgzB6tWr0bFjRyxbtgwAEBAQAL1ej3Xr1mH06NFQqVTYsGEDNBoN1q5dC5VKhT59+sDBwQELFy7EpEmT4OLigu3btyMzMxO7d++Gk5MTAMDFxQUTJ07E6dOn4eXlVYlvBREREVHNVerwqFKpxJdffomJEyeaTLezs0N+fj7y8/MRERGBgQMHmswfNGgQMjMzERUVBQAICwtD3759oVKppJrAwEAYDAaEhoZKNT4+PlJgAwB/f3+o1WocPXr00deSiIiISOZKDW02Njbw8PCAi4sLhBBITU3Fhg0bcPz4cbzyyitISkqCTqeDm5ubyfNat24NAIiLi8Pdu3dx/fp1sxpnZ2c4OjoiLi4OABAbG2tWo1Qq4erqKtXUVIK3HiUiIqJKVKZj2oocPHgQb731FgDg6aefxtChQ3HhwgUAgKOjo0mtWq0GAGRnZyMrK8tiTVFd0YkGWVlZpdbUNLz1KBEREVWFh7rkR8eOHbFt2zZ88MEHiIqKwsSJE2E0Gkt+ARsbiFK6oWxsSm9GWWqIiIiIaquH6mlr2bIlWrZsCR8fHzg6OmLOnDnSvJycHJPaop4xjUYj9Z49WFNUp9FoABT2xFmradGixcM0lYiIiKhWKbX7KiMjA7t378aNGzdMpnfs2BEAkJycDKVSicTERJP5RY/d3NygVqvh4uKChIQEk5q0tDTk5ORIx7G5ubmZ1RgMBiQnJ5sd60ZERERUl5Qa2oxGI+bOnYtvvvnGZHpYWBgAoHPnztBqtTh48KDJMGhwcDA0Gg06deoEAOjduzdCQkJQUFBgUqNUKuHr6yvVhIeHIyMjQ6oJDQ1Fbm4uevXq9ehrSURERCRzpQ6POjs7Y+TIkdiwYQMcHBzQuXNnREZGYv369Xj55ZfRtm1bTJkyBePGjcOMGTPw4osvIjo6Gps2bcKsWbNQr149AMCECROwb98+TJw4EWPGjEF8fDxWrFiBESNGSEOfI0eOxLZt2zB27FhMnToVGRkZWLZsGQICAtCtW7fKfSeIiIiIajCFKO0sARReYHfLli3YtWsXrl27hmbNmmHEiBEYP368dILAoUOHsHr1asTFxcHFxQWjRo3Ca6+9ZrKciIgILF26FBcuXICTkxNeeOEFvPnmm7Czs5Nq/vrrLyxatAjR0dFQq9Xo378/3nnnHYtnlZZGq9VKr1tZtv4ejw9/PIcnmzri0Mw+lfY6REREVLuVllvKFNrkiqGNiIiI5KK03MLraBARERHJAEMbERERkQwwtBERERHJAENbBam1BwYSERFRjcDQVk689SgRERFVBYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaCMiIiKSAYY2IiIiIhlgaKsgQvBGVkRERFR5GNrKS8EbWREREVHlY2gjIiIikgGGNiIiIiIZYGgjIiIikgGGNiIiIiIZYGgjIiIikgGGNiIiIiIZYGgjIiIikgGGNiIiIiIZYGgjIiIikgGGNiIiIiIZYGirILzzKBEREVUmhrZy4p1HiYiIqCowtBERERHJQJlCm9FoxI4dO/D888/D29sb/fv3x+LFi5GdnS3VnD17FqNHj4a3tzf8/f2xYsUK6HQ6k+XEx8dj8uTJ0Gq18PPzw/z5802WAQCpqamYNWsW/Pz80L17d8ycORO3bt2qgFUlIiIiki/bshRt3LgRq1atwvjx49GzZ0/ExcVh9erVuHLlCjZt2oSEhASMHTsW3t7eWLVqFWJiYrBy5UpkZ2fjww8/BADcuXMHY8aMQZMmTbBkyRKkpaVh2bJlSElJwfr16wEAer0e48ePR25uLj766CPo9XosX74cEyZMwHfffQdb2zI1l4iIiKjWKTUFCSGwceNGvPLKK5g1axYAoFevXnBycsKMGTNw4cIFbNu2DRqNBmvXroVKpUKfPn3g4OCAhQsXYtKkSXBxccH27duRmZmJ3bt3w8nJCQDg4uKCiRMn4vTp0/Dy8sK+fftw8eJF7N+/H+3atQMAdOjQAc899xwOHjyIIUOGVOJbQURERFRzlTo8mpOTg6FDh+K5554zmd62bVsAQGJiIsLCwtC3b1+oVCppfmBgIAwGA0JDQwEAYWFh8PHxkQIbAPj7+0OtVuPo0aNSzRNPPCEFNgDS46IaIiIiorqo1J42R0dHzJs3z2z64cOHAQDt2rXD9evX4ebmZjLf2dkZjo6OiIuLAwDExsZi6NChJjVKpRKurq4mNQ8uBwBatWol1RARERHVRY909ujp06exYcMG9O/fHw0aNABQGO4epFarpRMNsrKyKqSGiIiIqC566NAWGRmJCRMmwNXVFQsXLoQQJV9W1sam9JeoqBq5CruSil8u3KjuZhAREVEN9lBJaP/+/Rg3bhyaN2+OLVu2wMnJSeoZy8nJMavPzs6GRqMBUNgTZ62maBllqalt0nMKMGpjOMb/LwJXbrI3kYiIiCwrc2jbvHkzZs6cia5du2L79u1o2rQpgMKhSxcXFyQkJJjUp6WlIScnRzpGzc3NzazGYDAgOTm5xBqg8GQHS8e61QY3s/Kk/798I6saW0JEREQ1WZlC286dO/Hpp59i8ODB2Lhxo9R7VqR3794ICQlBQUGBNC04OBhKpRK+vr5STXh4ODIyMqSa0NBQ5ObmolevXgAKzya9fPkyYmNjpZorV64gJiZGqqmxePNRIiIiqkSlnj2alpaGTz75BI8//jhGjRqF8+fPm8xv1aoVJkyYgH379mHixIkYM2YM4uPjsWLFCowYMQItWrQAAIwcORLbtm3D2LFjMXXqVGRkZGDZsmUICAhAt27dAABDhgzBunXrMGHCBMycORMAsHz5cri7u2Pw4MEVve4VQsGbjxIREVEVKDW0/fbbb7h79y6uXr2KUaNGmc1funQphg0bhqCgICxduhRvvfUWnJycMG7cOLz55ptSnbOzM7Zu3YpFixZh9uzZUKvVCAwMxDvvvCPVqFQqbN68GZ988gk++OADqFQq9O7dG3PnzuXdEIiIiKhOKzUJvfDCC3jhhRdKXZBWq8W3335bYo27uzu2bNlSYk3z5s3x2Weflfp6RERERHVJ7b2OBhEREVEtwtBGREREJAMMbUREREQywNBGREREJAMMbUREREQywNBGREREJAMMbUREREQywNBWQXgXKyIiIqpMDG3lpADvY0VERESVj6GNiIiISAYY2oiIiIhkgKGNiIiISAYY2oiIiIhkgKGNiIiISAYY2oiIiIhkgKGNqo3eYMTw/x7Ha1tOQQhe6Y6IiKgkttXdAKq7Qi7dQmTCbQBAzK0cPNHUsZpbREREVHOxp42qTYHeKP2/kT1tREREJWJoIyIiIpIBhrYKwmOyyodvHxERUckY2spJwVuPEhERURVgaKMageGXiIioZAxtVC4ZuQX4POQK4lJzqrspREREtRpDG5XLzG9PY1nwJfRbfqS6m0JERFSrMbRRufx68SYAnkhARERU2RjaiIiIiGSAoY2IiIhIBhjaiIiIiGSAoY2IiIhIBhjaiIiIiGTgoUPbhQsX4OnpiZSUFJPpoaGhGD58OLy8vNCvXz8EBQWZPffs2bMYPXo0vL294e/vjxUrVkCn05nUxMfHY/LkydBqtfDz88P8+fORnZ39sM2UDZ51SURERGVh+zDFMTExmDRpEvR6vcn0qKgoTJ48GYMHD8a0adMQGRmJpUuXQgiB8ePHAwASEhIwduxYeHt7Y9WqVYiJicHKlSuRnZ2NDz/8EABw584djBkzBk2aNMGSJUuQlpaGZcuWISUlBevXr6+gVa4czF5ERERUmcoU2vR6Pb755hssX74cdnZ2ZvNXr16Njh07YtmyZQCAgIAA6PV6rFu3DqNHj4ZKpcKGDRug0Wiwdu1aqFQq9OnTBw4ODli4cCEmTZoEFxcXbN++HZmZmdi9ezecnJwAAC4uLpg4cSJOnz4NLy+vClz1ilHeuy/x9k1ERERUFmUaHo2MjMS///1vvPbaa5g9e7bJvPz8fERERGDgwIEm0wcNGoTMzExERUUBAMLCwtC3b1+oVCqpJjAwEAaDAaGhoVKNj4+PFNgAwN/fH2q1GkePHn20NSQiIiKqBcoU2tq1a4fDhw/jjTfegFKpNJmXlJQEnU4HNzc3k+mtW7cGAMTFxeHu3bu4fv26WY2zszMcHR0RFxcHAIiNjTWrUSqVcHV1lWqIiIiI6qIyDY82btzY6rysrCwAgKOjo8l0tVoNAMjOzrZaU1RXdKJBVlZWqTVEREREdVG5L/khSjn90cbGpkw1pSlLDREREVFtVe4kpNFoAAA5OTkm04t6xjQajdR79mBNUV3RMhwdHa3WWOqBIyIiIqoryh3aWrVqBaVSicTERJPpRY/d3NygVqvh4uKChIQEk5q0tDTk5ORIx7G5ubmZ1RgMBiQnJ5sd60ZERERUl5Q7tNnb20Or1eLgwYMmw6DBwcHQaDTo1KkTAKB3794ICQlBQUGBSY1SqYSvr69UEx4ejoyMDKkmNDQUubm56NWrV3mbSkRERCRbFXKg2JQpUxAVFYUZM2bg6NGjWLVqFTZt2oRJkyahXr16AIAJEybg1q1bmDhxIkJCQrB582YsXrwYI0aMQIsWLQAAI0eOhEqlwtixY3Ho0CHs3LkTb7/9NgICAtCtW7eKaCoRERGRLFVIaOvZsyfWrFmDmJgYTJ06FXv27ME777yD119/Xapp164dgoKCkJubi7feegubN2/GuHHj8P7770s1zs7O2Lp1Kxo1aoTZs2dj5cqVCAwMxMqVKyuimZUqIS0XBiPvi0BERESV46FuYwUAL730El566SWz6QMGDMCAAQNKfK5Wq8W3335bYo27uzu2bNnysM2qEb45lYSRfq2quxmyxHuwEhERlYzX0Sin4rehik68XX0NISIiolqNoa2aXb+TV91NqBF4D1YiIqKSMbRVoEcZ4Ru3+VSFt4OIiIhqH4Y2qhF4TBsREVHJGNqIiIiIZIChrQKxt+jR8Zg2IiKikjG0EREREckAQ1sFEo90KgIRERFR6RjaiIiIiGSAoa0isaONiIiIKglDGxEREZEMMLRVIHa0ERERUWVhaCsnBXitCiIiIqp8DG0VSPBCbURERFRJGNqIiIiIZIChjYiIiEgGGNqIiIiIZIChrQLxiDYiIiKqLAxtFYjnIRAREVFlYWgjKge9wVjdTSAiojqCoa0ClbejLTY1p0LaQVVjz+lr8JwfjJ0RSdXdFCIiqgMY2mqQZcGXcDMrr7qbQWX05o5o5OuNeHvXmUd6vtEosOaXy9h/9noFt6zmE0IgMiGd+zsR0UNgaKtAFXFx3TNJdyqgJSQH+85ex/JDf+Gf26OquylV7silWxj+39/h+8kv1d0UIiLZYGijSpeanY+g0DjczKxbvSpCiBLX+fKNLJPaumRnJIeUiYgeFkNbeRW79aixjn3xltX4/0XgX3vPY+TG8OpuSpVa8vMl+C76Bf87Hm+5QHF/5zHWsV2H9+wlInp4DG0VaP/ZlHIvQ87f3Tn5eovTTydlAACu3MyuwtZUv3VHYwAA8386V2qtXAP/I/cQlpLZhBCISryNO3d1j7Z8IqJaiKGNKsy4zace+bklffen5xQgT2d45GXXVMVzi6XQZjQKfH0yEWeSM6qsTQ/j65OJ8P74EI7+davCl/1D9FW8tPY4Bq08VuHLJiKSK4Y2mdh75hp+vXijuptRopPx6RW+zGsZd9Ht40N4etmRWnfcV7HRUYuhdVdUMuZ+fxZDPwurukY9hLnfn0VGrg5jgk4+9HNLGxz94rc4AEBKHTsOkoioJAxtMnDheibe+Coar22JQHpOQXU3p1L8kXTb4vSdEckACr+8a/NxX5Z62kIvp1ZDS6qGQlFybLPhIW9ERGZqZGjbu3cvnn32WXTp0gWDBw/G7t27q7tJVcbSd9XlYseCpWXnW3xens4AYw1MNYYytmnOd2eRmWd+/FLxL2+9sXbdfaD4wfiW3qZScg2CQuMw8osTuG0lyGfm6ZCQZv2CzVduZmH8llM4HlM54TAyIR3zdp/FrSzzfba0TGZT2sqXkxACKXfYi0dE8lLjQtv+/fsxe/Zs+Pv74/PPP4evry/mzJmDn3/+ubqb9lDSsvNx5WZW6YUPsBRxin99Wfouu51TgB6Lf8ErG36vcUOIF65nlrk2tdiXe9F6FF/fsgbA6nIt4y4+D7lS5kubKB4487hAb8S735/BN6cSC+eX8vx/7T2P4zFpWPLzRYvzeyz6BX2WHbG6DcYEncIvF29i5BeWz+o1GgUupWQ98o+B4f/9HdtOJGLmt3+YzbO0H+uK3RKsLJlN3HvPHsXaIzHosfgXfB5y5ZGeT0RUHWpcaFu5ciUGDx6Md999F0899RQWLFiAwYMH4z//+U91N82izAfObvvyRALe/f4sui88jP4rjuH7qGQ8s/wIlh+8VOZAJYRA8u3cMtdvO5GAjFwdTsXfht7CF2zIpZt4fk0oYm5ZPntz9S+XEbjqWIk9D7dzCsrUnne/P4M3voqSzvp78Mv3+JWSe3WEEPjy93h4LTiIsCupJsNoD65bgd6IRfsvYPSmcEQnWh5eLYt3dp1Gv+VHcCe3fGcq/n1jOJYFX8IYKydk5OkMiEq8bTF8CmPhdtxxMglzvjsLoPQhxCJXM+5anJ5bUHjyRtEQc1mfV+TTny9i0Kpj+Hjf+VLbIITAidg0TP0qColpuSbzjsekmdU/uGYxt7LR7V+HMPe7wrtLlGXdx2w+Bb9Fh5Fqpfe5JMuCL5n890F5OgPCY9N4b1l6aAajwN2C2nfiFNUMNSq0JSUlITExEQMHDjSZPmjQIMTGxiIpqeZdkLPosg5FPtj9J3acTJQez/z2NGJu5WDNr1fwxo5obP09Hu9+fxYvfB6GxQcumC3v9a0RmPhlJPyXhOCzX6+YvUb/Fcfw+wNfgtfu3P/yLQoExQPWuM2ncPbqHTyz/KjFdVhx6C9cTMnCx3sLv5yNRiF9ERbojXCfdwDeHx/CmzuipWXn6y1/KO04mYS9Z67Da8FB3Lmrg/KBg5NGbgzH9XvtTc81Hdb7+lQS3N7djw9+PIfMPD1GbQzHidj767rkwEWMCTophYLt4QnYcCwWv11OxYtrj0t1l29k4UYJvV3FQ1O+3oBvI5IReysHG36LsViflp0vvZ86K1/ibebuk+4da6lnKyO3AOM2n8JLa4/jk32F2714Bp72TTQ+LdZjJoQwCTa5BXpcv3MXA1YcxYpDf5ksu7TLhYh7/be5BfqH6v3dcCwWALA5LB5A4UWST8SmWXwPPvv1Cl7dcAL7zlzHhK2mobW0HtLk27lYtO8CsvL1+PpU4d+4pcgWl5oj9awJIXDsr1u4navDF7/FmtUKITB1exSGfR6GAr0Rf93IwuQvI032p5LM+OYPvLLhBP59sPC9ztMZEHo51eQsZiEEvgpPxMk4yyfg3M4pwG+Xb1n9sWM0CmRbuUwOULivhcemWf1bC7l0E3tOXyvT+ly5mY2kdNMwHXYlFRP+F4FYKz/mSnMt467VtpWFzmAscf3lSAiB/1t3HH6LDls9lKU0V25m4aOfzpltr4dtR3U8Fyj8zrC2XwghSvxsDruSiinbIs1++D2M8v7QKtAbkVtQc/dLhahB42lHjx7FxIkT8eOPP6J9+/bS9PPnz+PFF1/EF198gYCAgDIvT6vVAsKIiMM/WC9SNwZU6sL/NxqAO5Z7JSSOLoCdQ+H/G3Twn7e9xPIbwgk62AIA7FGAJoqSb1N1XTjDACUAwAH5aKwoeXgxWTRG0VdcfeTBWZGFFg0dcO1er1nAk41x7N4B7UIAV9FEeq4GuWiosH7Mk17YIAWPSY8bIBsNFNZ7ZwqELW7CSXrc0iEXIt/y8ts1VuP8rQLcQiNp2mO4g3oK6yda5Ap7pKOB9LgJMmCvsN471qtja3x77v4fvwvSYaew/iWTKeojE2rp8dA2RkQlWP6SV6uUuJrvgGzUl6a5Kkq+9EW60CAXhfuOAkY8rig5QNwSDZEPFQBACQOaK6yfnWujUCDF2BAFsAMANLAzooHe+vJHaFviP6dypH2tLPvmVfEYxL3feUX7WkmSxf19zRG5aFTCvmYQNrh+b18b5OmC38/FooHC8ge3xt4W6XkCN+AsTWuELDgqrH8Z5As7k33NGZmoryj8Um39WH0kFPuSUCoUyDbaIQ0NpWml7Ws6m3q4YXCUHjfFbagU1j/4M0U9ZMJRen2H3BQIox65BQY0cbRHY0cVLqTcf3/vCDWyiu1rj+NWiUPIxfc1QMBVUfgZoFIqUGAw/8hXOzXDpfTCv422zvZQ59/E7Xs/quxtbZD/wDB08X3TDnq4KAp7ur0ebwhHBzuEFTtO0lmtwvlsNR5roIbGwQ5JN9PRrl4OMvPuvz9+bs4ILxZ+r4nHYLy3r9VDHtrVz8MTTR1x+Wa22bX7XDQOyFU/jkspmTAKQI27cFJYDqIqpQLPd3XFfyLu7ysNkGNxX2vgYIvMPL3J52ADB1so8jKgsfA52LyhA7RtnLDrj1u4Vexz0AmZUCssBzi3xvWhcWyAk7ds0aCeLWJv5aAx7sChhM/BHGGP2/c+Bwd5ukB/JwWXrlofwcgS9XAH9/fNZkiDrcKIdo3ViEvLNfvR9+DnYAukwkZhOSZ4t2yEEykCtwpU0rTSPgfThAZ37+2bNjCihSINTTT2aFTPzuT47S6PN8SZq3dM9jVHOyMaFftc83Vzxsm4dLRrrMYL3o9j+aG/rH7nOtVXSft0kQEdmmLreT0MUMJZrUJuTpb0ndvl8YZwqq/CqYR05BYYoG3thDy9EbNffgZPuNz7HsrPBnJL+BxXKIBGre4/zrsD3M2wXm9jC+0zwwAAERERFktsrT+76mVlFX5IOTo6mkxXqwt3oOzsR/hFmJ8F/KeL9fkjtgIdC98k5N0puRYAxuwF3J4q/P/bCQi1n1Zi+aD8T3FJFG40T0U8vrf/qMR637zPpeDT2+ZPbFItL7H+ybyt0g46RBmOf9utB/IA2N8rSLz//9nCAZ3yg6TnvqIMwTw766EzWTSGf/5q6fEk272YavuT1fqzxjZ4vmCR9HiKfjtG2odYLs4Cjtl1xj9070qTFtoFYbDS+rXefjT0wjTdG9LjNao16GFj3ltZZMulgQDGSo//p1qC9jbWe2tX6V/CKv3/SY/fv/4GXOwzrNbPV47B/wyDpMe/qmZCVUIofKtgKn4y9gYA1Ed+qfvO6IK5+M1YuD82V6SXWv9c/kL8KdoCAFrp47HXfp714rPAt4r/SMHK1+YivlR9WuLyPfM2IQf1AAD9bSKxWvW51doCoYR7/pfS4+HK37DA7n9W62+IRvDLXwsACD53A9NtD2C67fdW6y+qWiKwYIn0eLrtdxhre9Bq/QljB7xa8IH0eL7dVgxT3uudzcb9v5d7Dhh8MEU3Q3q83O6/CFCetbr8r/R98R5elx5vUi1DZ5t4q/Wf64dimf5VAEBCWi5C7ecVBit7ADoAt03btFA3ChsNz0qPg+3nlBhSZ+smYZehDwDADgbTfcfCp/74jFm4hO4AgOz0FPzqMNX0PVGa1r+U/xGihDsAoK3iGoLt5xbOKMoOxZ+rB/oqliMu0xY3MvPRw+YKvhYLTWuumT6na956ZEADAOhjcwbrjKuAomuXP7CtUAC0yfxKejhUeRyL7TaZr+Q9t8864j/YID3+u/IQ3rH71rxQFL5WrLEZ+hWsAABk5unxvu1uvG6737w+D8BFoL/qSQwvWCBNft/uK/yf0sr1BrOAwxne2K97Wxrd+NRuA/oro622f5chALN1kwEU/q18p/oI3e0vW63/Qj8En+j/Lj3+SvUJ2tqkAFkAVOb1S3UjsNbwgvR4n/17VkMwbgLvGsZjB56RJpX2OTW5YDp+NvoCKAzMofbTgAIU/iu+bVMLH79aMA8njB0BAE30NxBiP+t+TdF+kwXgGPCifQnfuQaY7zuxwD4Ufuem5xTgmeLfuUX7suLe8+7tf0+u3IoLnzwPW6UNcP5H4Md/Wl9ZlSPw3tX7j6O2AgdL+Fxu2MpCI03VqOHR0jr9bGxqVHOJiIiojnnwsJ+qVKOGR48cOYJJkybhp59+goeHhzT93LlzeOmllxAUFITevXuXeXlVMTyKzGuY8e0fOBVv+UD4qhwebazSwbOR3uoJB/VUtricf7/bflh7NSIvxVtddnmHR0saFgCA9o83xuHk+zt/eYdH7WwU0BU7fipbOEi/1oGSh0dtbRSo18AZF27f/2HQHGlQKqwfH5Eh1DVmeBQAbopG0vCoCjo0VWSUWF98X3tw37Q0lFBVw6OA9SGrIjqhfOTh0aYae9TT3YbBytA9ANwVqhKHR5s1cEDf9k2w42Rhz23xIavZA92x9WC4NDz61BONoTMacSL2/vYrPjwKAK1sb8No0MPBzgZ9PZriwJ+mt8TLs9UgVecgPX5wePRpjyY4cun+/vfg8KifUzauZtx/fx4c8kwVDZB37xf+g/taU409bmblo3lDB1y/d9hF8X0zoF0DxMYWHg9qo1CYDbc9plZB28UTm44XfraW5XPwweHRx8q4r/Vr3xThFxOkniFtaydEJJh+NhuFAtfQWHrcubFATkaaxWFjoPBzsFOHDjgek4aebR9DxMVYaXj0yXtDtsU1dFTjXNb9z4VpPZ3x3QnLJ7wAQJ5QIbXYvlZ8eNTWRmF2Albxfc2pvh1sc2/CvoSheGvDo891bo7GGnscOp9ism88zPAoANwWjshBPfi6OcPPzRk/hPwuzVPA/IoIloZHi4bt1fZK5OTf/4x+b0gHfPjLTaTmFe4LrRrawZhZeBxnxxYNkHInz+zapUXfuY3q28HPVY1zl633QgLWv3N92jiZfa+P7tEKr/TvjUbqe71h1TA8WqNCW0JCAgYOHIjPPvsMAwYMkKbv378fM2bMQEhICFq0aFHm5Wm1WgDWV74iTf86Grv/MD0oOP7TwuGMooP5H1Q0/1rGXfxrz3m8rHXF6aQMrL53AsLHwzzx0Z7z0oHcH7/QCalZ+Th84Qbe7PcEAjs1R5u5+8yWV9yD841GgT1nrqFj8wZ40kWDURtPIOxKmsn88Lh0eDTTwFmtMnt+vt6AO3d1aKop/KP78kQCtv2egH/2bQefNs7YcCwWT7o4okBvROvH6uO1Lfff+wv/CkR4XBo6tmggPR8Apn4VhX1nrpu0+9yCQTh3LRMj1hd+AMQsGoJdkUlYtP8ixvZqgxkD3E3Wr3lDB/z+7jN40NojVxB87ga++Ed3NNU4QG8wFnZrA/D+10HcvnfGqKX3rstHwdJxN8W3pVEIONgp0Xl+MLKKHUi9eawPjKLw/XumfVMk3b6L2TtPS/OPvv00Wj92/8Ow+HsLAFc+GYzE9FyobG3g6lQfm8PisGDP/TM34z99Fnk6Axzs7o9VrT1yBUt/vgQ/N2d8M6kngMID3G3u/RL0X/Irkm/fRcjsp9G8oQNy8vVwVqugUCgghMBrW04hIS0X+6c9ZbJcoOiyNdm4mJKFZ7s0R2NHe+nyH3vPXsdbO+4P4UwMaCuduHDlk8H4PCQGKw8XHsT/XJfm+PD5jtAbBOyUNmiiscf28AS8/8OfJq/30fMdceDPFHw2shuaaOxN3p/lL3uhX/umcFIXBoVfLtzA+P9F4NnOzfHZSG+zs03zdAYs3HceY3u1wRNNNWbvCwAc/esWfv4zBU97NMEgz2Z4UIHeiIspmejUoiEEgMT0XLR5rL7VM1uFENK8ora7OtVD6Jx+FutLMm7zSYTcC2GW9s0iD65TRUnPKYDaXgl7W6XF+UII5BQY4Ghv+Qib4zGpOJN8B68/1dasV+L6nbtY9vMlvKxtiZ7tHjN77i8XbmD6N39g3rMd8IpPK7P5Ra+ffPsuXJ3qlfks6wddzbgL5/oq1FNZXseSCCFQYDBafX+y8nTIyTegWUMHs3n3L2f0aO3OyC2Ag53S7O8VKDzpZ8OxWHi2aIAA98Igm5CWgyYae9RXle1oKKNRIF9vtPq+FN/P5UYIAb2x8HPIkh//uIqrGXcxpU+7alvH0nJLjQptAPDMM8+ga9euWL78/rFc06dPx4ULFxAcHPxQy6rK0GY0CmTc1SH2VjbGBJ3ExIB2mNb/SWn+ybh0pGTmoatrI/znl8v4e49W8G7lZLacm5l58F30CwAg/L1ncDwmFTO+KfziPzDtKXRo3sCkvqyhTa1S4ty/As3mv7rhd6kHwNLzey7+Rfp1XdKXhyV/Xr2D59aEltg+oOiMonw01dhjz5lreExtD/8nC38J//xnCpo1dEDXlo0sPrdo/YqHlrLq9vEh6VeapbadSc7A3O/OYmrfJ/Bsl+Zm87+NSMI7u85IjzeP80Ffj6YW2wcAJ99/xiSsFgWPIpba8F1kMmbdC37W3r/k27lwaeBg8YMoT2dAZp7O5HUf9Cgfwj/+cRXTvv5Denzx40BcTMlCo3p2aNNYjRWH/sLqXwp/4f61cDBUtqZtu5GZB797+zkAuLs44uCMPiY1/9weif1nC3ucLK37zaw8NHG0r5FfIM+u/g3nrmXiv6O6YXBn832nNLey8vF5yBUM69rC4udEbVdZYZSopistt9SoExEAYOrUqXj33XfRsGFDPP300/jll19w4MABrFy5srqbViIbGwWc1So4q51xev5AqTeniK/b/aGc5SO8rC6naQMH7JxcGD5cGjign4eLNM+tsdra00r11JNNLE4v7SYDn430xsgvwjG2d5tHfu3SKBQK6RfpsK6Pm8wL7GTeA1Lc0uFd8G1EEv79svX39FF1cW2E/dOesjp/hLYlIIB37l1brLSr+D8YnJ7p4IL2zTS4mJKFdwI9LD7nRe/HYatUoH2zBhbnA4CrU32r86z9Ii+uvKFnw+jucLBTmgRrlfL+Mu2U5st3aWD6Xlh67z4a6ol8nRFDu1ruXS8piFa3XZN7ITE9Fx7NNKUXW9BEY4+PhnpWcKvkg4GNyLIaF9peeuklFBQUICgoCDt37kTLli2xZMkSDBkypLqbVmYPBraH5dPmfsBrWN8OkfP6w1ZpY/HLt+hyBX/vYXkYYfsEP+w9cw2zB1oOBdP6P4lRG8PRsbnlUNC9tTPOfDTQ6jBASRzsKv/EkRE+LTHCp+UjPXeoVwtsOR5frtdv2+R+kH6Ur5nvpvTC+euZ6GalN8XGRmEWZGuC4v3zT7qYB5PRPdvgm4gkdHFtVKZQWN/CUExTjQM2jfUpVzurSz2V8pEDGxGRNTUutAHAq6++ildffbW6m1FjPOZo/RTg76b0wu8xaRjQ0cXi/N5PNEbvJxpbnFc0/9jbfeHS0PprPEpgA4B2TRxLL6pGcwLbo20TNfxLeH9KU/zYAkvZ5POR3bAs+CKWDLd8KRm1va1JSJeL4gebW4pkDevZ4djbfUsMbB2bN8D5exciXvp/pVxqh4iIamZoo7Jr7GiP573KfnKGJa0esz68Vh418Vij4uqplPhHzzblWkbzYgcat7QwTPlsl+YWj4eTu+I9bdY2c2nb//t/9kJSei6eaOpY4/cVIqKagKGNKtWkgLZYfywWrz7iEGZN5+pUH5+P7AaDEGhTjmMO5aZ4D2Npx/JZ42CntDi0SkREljG0UaWaE9geQ7u2KPFAermrjT1ppalhJ50TEdUJDG1UqWxsFPBs0bD0QpKV4iG8YX27amwJEVHdwdBGRA+ts2tDrP6bNxrWs0MDB4Y2IqKqwNBGRI9kaDlPgCEioofDO7ATERERyQBDGxEREZEMMLQRERERyQBDGxEREZEMMLQRERERyQBDGxEREZEM1OpLfmRnZ0MIAa1WW91NISIiIipRVlZWifdirtU9bTY2NrwRNREREcmCQqGAjY31aKYQvIkgERERUY1Xq3vaiIiIiGoLhjYiIiIiGWBoIyIiIpIBhjYiIiIiGWBoIyIiIpIBhjYiIiIiGWBoIyIiIpIBhjYiIiIiGWBoIyIiIpIBhjYiIiIiGWBoIyIiIpIBhrZy2rt3L5599ll06dIFgwcPxu7du6u7SVQKo9GIHTt24Pnnn4e3tzf69++PxYsXIzs7W6o5e/YsRo8eDW9vb/j7+2PFihXQ6XQmy4mPj8fkyZOh1Wrh5+eH+fPnmywDAFJTUzFr1iz4+fmhe/fumDlzJm7dulUl60nm3njjDQwYMMBkWmhoKIYPHw4vLy/069cPQUFBZs+rqP2BqsapU6fwt7/9DV5eXvD398fHH3+MnJwcaT63ee2yY8cODB48GF27dsXzzz+Pn376yWR+rdregh7Zvn37hIeHh1i0aJE4duyY+PDDD4W7u7s4cOBAdTeNSrB+/XrRoUMH8e9//1uEhYWJbdu2CV9fX/Haa68JIYSIj48X3bp1E+PHjxdHjhwRmzZtEp06dRILFiyQlpGRkSECAgLE8OHDxeHDh8U333wjtFqtmDhxolSj0+nE0KFDRf/+/cX+/fvFTz/9JPr06SOGDh0qdDpdla93Xbd7927h7u4u+vfvL02LjIwUnp6eYvbs2eLo0aNixYoVwsPDQ2zcuFGqqaj9gapGdHS08PT0FG+++aYICwsTX331ldBqtWL69OlCCG7z2ubrr78W7u7u4tNPPxVhYWFi0aJFwt3dXezfv18IUfu2N0NbOfTv31/6ICgybdo0ERgYWE0totIYjUbh4+MjPvroI5Pp+/btE+7u7uL8+fPivffeE3369BH5+fnS/O3bt4sOHTqIlJQUIYQQn3/+uejatatIT0+Xao4cOSLc3d3FH3/8IYS4HxKuXLki1Vy+fFl4eHiIffv2VeZq0gNSUlKEj4+PCAgIMAltY8aMES+//LJJ7dKlS4VWq5W2f0XtD1Q1Ro0aJUaNGiWMRqM0bdu2beKZZ54Rubm53Oa1zCuvvCJGjx5tMm3kyJHi73//uxCi9v2Nc3j0ESUlJSExMREDBw40mT5o0CDExsYiKSmpmlpGJcnJycHQoUPx3HPPmUxv27YtACAxMRFhYWHo27cvVCqVND8wMBAGgwGhoaEAgLCwMPj4+MDJyUmq8ff3h1qtxtGjR6WaJ554Au3atZNqih4X1VDVmDdvHnr37o2ePXtK0/Lz8xEREWHxbzgzMxNRUVEAUGH7A1W+9PR0RERE4G9/+xsUCoU0fdSoUTh8+DBsbGy4zWuZ/Px8qNVqk2mNGjVCRkZGrfwbZ2h7RLGxsQAANzc3k+mtW7cGAMTFxVV5m6h0jo6OmDdvHrp3724y/fDhwwCAdu3a4fr162bb1dnZGY6OjtJ2jY2NNatRKpVwdXUtsQYAWrVqxf2jCu3cuRPnzp3DBx98YDI9KSkJOp2uxL/hu3fvVtj+QJXvr7/+ghACDRs2xPTp09G1a1d0794d8+fPR15eHrd5LfSPf/wDv/32Gw4cOIDs7Gz8/PPPOHLkCIYNG1Yrt7dtlb1SLZOVlQWgMAQUV5T4eTCqfJw+fRobNmxA//790aBBAwDm2xUo3LZF2zUrK6tMNU888YTFmoSEhIpcBbLi6tWrWLx4MRYvXgxnZ2eTeWX5G7ZWU1T3MPsDVb709HQAwNy5czFgwAD897//xaVLl7Bq1Srk5+fjlVdeAcBtXps8++yzOHHiBKZPny5Ne/HFFzFhwgRER0cDqF3bm6HtEQkhSpxvY8NOTDmIjIzE5MmT4erqioULF6KgoKDE+rJs14qqofIRQuC9995Dnz59MGjQIIvzS2JjY1Mhf+fc1lWn6Gy/bt26Yf78+QCAnj17QgiBJUuWYMSIESU+n9tcfqZMmYLo6Gi8++676NixI06fPo21a9fC0dERQ4YMKfG5ctze3LMekUajAQCT08iB+z1sRfOp5tq/fz/GjRuH5s2bY8uWLXBycpJ+ST24XYHCbVu0XR0dHa3WFC2jLDVUebZv345Lly7hvffeg16vh16vlz6g9Xp9mf6GK3J/oMpX1IMSEBBgMt3f3x9CCJw9exYAt3ltERUVhdDQUMybNw9jx46Fr68vXn/9dcydOxdffvkl6tevD6B2bW+GtkdUNLadmJhoMr1o2MvSsUxUc2zevBkzZ85E165dsX37djRt2hRA4Ye+i4uL2fBlWloacnJypO3q5uZmVmMwGJCcnFxiDVC4z3D/qHzBwcG4ffs2/P394enpCU9PT+zevRuJiYnw9PREREQElEql2d9w0WM3N7cK3R+o8rVp0wYAzHrMi3rgXF1duc1rkWvXrgEo7FktTqvVAgAuXLhQ67Y3Q9sjat26NVxdXfHzzz+bTD948CDatGmDFi1aVFPLqDQ7d+7Ep59+isGDB2Pjxo1mvaK9e/dGSEiIyQd/cHAwlEolfH19pZrw8HBkZGRINaGhocjNzUWvXr0AFP66v3z5snTSCgBcuXIFMTExUg1VngULFmDXrl0m//r27YtmzZph165dCAwMhFarxcGDB02GSIKDg6HRaNCpUycAFbc/UOVr164dHn/8cezfv99kekhICGxtbeHt7c1tXosUhaXIyEiT6X/88QeAwqsC1LrtXWUXF6mFvvvuO+Hu7i4WLFggjh49Kl1cl9fgqrlSU1OFl5eX6Nu3rzh16pSIjo42+ZeWliauXLkiOnfuLMaMGSN+/fVXERQUJDp16iTmz58vLSctLU34+fmJYcOGiYMHD4pvv/1W+Pj4iAkTJkg1+fn5YtCgQaJv375iz549Ys+ePeLpp58Wzz33HC+uW03mzJljcp2248ePCw8PDzFt2jRx5MgRsXLlSuHh4SE2bNgg1VTU/kBVo+ii57NmzRJhYWFi/fr1wtPTUyxatEgIwW1e2/zzn/8U3t7eYvPmzeLEiRNi48aNwtvbW9oOtW17M7SV044dO8SAAQNEp06dxODBg8UPP/xQ3U2iEvzwww/C3d3d6r/du3cLIYQ4deqUePnll0WnTp3EU089JZYvXy4KCgpMlnXp0iUxZswY0aVLF9GzZ0/xwQcfiKysLJOaa9euialTp4quXbsKX19fMWPGDHHjxo0qW18y9WBoE0KIgwcPiueee054enqKfv36iU2bNpk9r6L2B6oahw4dEi+88ILo1KmT6NOnj1i7dq0wGAzSfG7z2iM/P18sX75c9OnTR3Tq1EkMGjRIrFmzxuRCubVpeyuEKOXUCSIiIiKqdjymjYiIiEgGGNqIiIiIZIChjYiIiEgGGNqIiIiIZIChjYiIiEgGGNqIiIiIZIChjYiIiEgGGNqIiIiIZIChjYiIiEgG/h9kx3x/DJA2yAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,4))\n",
    "ax.plot(spectrum)\n",
    "ax.axhline(np.mean(spectrum), color='C1', linestyle='--')\n",
    "ax.set_title('WHT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_sorted = np.sort(np.abs(spectrum))[::-1]\n",
    "# get normalized cumulative energy\n",
    "cumulative_energy = np.zeros(N)\n",
    "total_energy = np.linalg.norm(spectrum_sorted)\n",
    "for ind in range(N):\n",
    "    cumulative_energy[ind] = np.linalg.norm(spectrum_sorted[0:ind]) / total_energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Cumulative Energy (Top 100 Coefficients)')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEXCAYAAABS7RbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA71UlEQVR4nO3deVxU5eIG8GcGBBUIcy0TkdTBK4iiiCmQu4CK11DDJbDSXKJbLpSOv27ptUzMJLfSMk0zTdGrqSwSKq6ZYmVpajdRGXcJF8AEhnl/f+AZZpiBGRSE0zzfz4cPzHvec847h5l55j3LexRCCAEiIiKSFWV1N4CIiIgqjgFOREQkQwxwIiIiGWKAExERyRADnIiISIYY4ERERDLEAK/BcnNzsXLlSoSHh6NTp07o0KEDhg4dig0bNkCn01V388rk6emJ6dOnP9C8Go3G6HFkZCR69epVGc2yyuLFi+Hp6Vnuzz//+c9H1p5HTafTYciQIdi2bRv++9//WtwWnp6ej/T/I9m+fTs8PT3NTsvLy8OcOXPQvXt3dOjQAVFRUTh58uQD1yuLEAIJCQl48cUXERAQAG9vbwQHB2PBggXIzc194Odmre3bt6NXr15o164dpk6dioKCAqjVanTs2BEdO3bE7t270atXL0RGRlZouY/iPZebm4vs7OwKzZOeno4ePXrg7t27VdQq+bGv7gaQeRkZGZg4cSIuXbqEsLAwDBkyBAUFBUhNTcU777yDo0eP4sMPP4RCoajuplaazZs3Y9asWfjll1/0ZRMmTMBff/31yNsyYcIEPP3002an1atX79E25hFav349CgoKEBYWhosXL2LevHn6aTdv3sQHH3wAPz8/PP/88/pyJyenR9rGU6dOYebMmWVOnzJlCg4ePIjIyEg89dRTWLt2LSIjI7Flyxa4u7tXuJ45ubm5iImJwZ49exAYGIiXXnoJdevWxbFjx/D5559j9+7dWLt2bZW9Vm7evAm1Wo1mzZrh7bffhru7OzZu3Ij//ve/+Oc//4nOnTvD29sbM2bMQJ06dSq07Kp+z504cQITJ07E/Pnz0aVLF6vn8/PzQ6tWrbBkyRK89dZbVdY+WRFU49y7d08EBwcLf39/cerUKZPpM2fOFCqVSqxevboaWmeZSqUS06ZNq/B806ZNEyqVqgpaZL1FixYJlUolDh8+XK3tqA45OTmiU6dOYvv27WanazSaB/7fVpbvvvtOdOzYUahUKrOvlQMHDgiVSiU2bNigL8vKyhKdO3cWU6ZMqXC9srz11luiTZs2YsuWLSbTEhIShKenp4iOjq7gs7PesWPHhEqlEmvWrNGXvf3220KlUomcnJwqW29l2Lx58wO/x44ePSratm0rMjMzq6Bl8sNd6DXQunXrcO7cOajVarRp08Zk+rRp0+Dq6opvvvmmGlpHf1ebN2+GVqtFnz59qrspZr3zzjuIjo6Gu7s7AgMDzdZJSEhAnTp1MHjwYH1ZgwYNEBISgl27diE/P79C9cw5efIktm7disGDBxvNL+nfvz8CAgKQlpaGq1evPtBztaSwsBCA8d4PqczZ2blK1lkT+Pn5oXnz5li7dm11N6VGYIDXQAkJCahbty4GDBhgdnrt2rWxceNGbN26VV9W1rGu0uW9evXCf/7zH8THxyM4OBg+Pj4YMmQIfvnlF9y4cQNvvPEGfH19ERQUhAULFhgday/r2LalY96FhYVYvnw5Bg0ahPbt28PHxweDBg3Cpk2b9HWkXZell2d4PO6zzz6Dp6en2eOUvXr1QlRUlP7xH3/8gejoaPj5+aF9+/YYPnw49u/fX2YbH1RkZCTGjBmDffv2ITw8HO3atUP37t2xePFik/MUrGmTtLy4uDj4+vqia9euOHPmDADg+PHjiIqK0v9/Fi9ejCVLluiPBe/fvx+enp74+uuvTdo5efJkBAYGoqioqMznsm7dOgQGBqJ27doV3g6XLl3Cm2++iWeeeQbt2rXDoEGDsHHjRqM606dPR9++ffHTTz8hPDwcPj4+CAkJwfr1661ax9mzZ/H666/jm2++QaNGjczWOXHiBFq3bg0HBwejci8vL/z11184e/ZsheqZk5CQAACIiIgos86cOXNw6NAhPPHEE/oya7YRAFy9ehVvvfWWvt7gwYOxbds2/fTp06frX+tqtVp/LoLh+0d6z5v7XDh+/DheeeUV+Pn5oUuXLhg3bpz+NQaYPwZekdduee+FxYsXQ61WAwCioqL06xFCYMmSJQgODka7du3QrVs3vPnmm7hy5YrJ9unXrx82b96Me/fulbn9bQWPgdcwQgicOnUKHTt2RK1atcqs16JFiwdeR2pqKlJSUjB69GgIIfDpp5/iX//6F1xcXNC6dWtMnz4dKSkpWL58OTw8PPDcc8898LqA4g+ZpKQkjBgxApGRkbh58yY2btyI//u//0OjRo3QvXt3TJgwATqdDunp6Zg3bx6aN29uspyBAwdiwYIFSEpKgpeXl778+PHjuHTpEiZOnAgAOHPmDEaOHImGDRti/PjxqFWrFnbs2IFx48bho48+Qv/+/S22OScnp8yTbJydnY0++H///XdMmjQJERERiIiIwI4dO7BkyRLUr18fo0aNqnCbfvzxR2g0Grz55pu4ePEiWrVqhRMnTiAqKgoNGzZEdHQ0/vrrL6xZswZKZcl38G7duqFBgwZITk7WrxcA7t69iz179mDo0KGws7Mz+5zOnz+P8+fPY+zYsRa3TWkajQbPP/888vPz8cILL6BRo0ZISUnBv//9b5w/f97oeOWtW7cwduxYdO/eHeHh4UhJScHMmTNx584djB8/vtz1rFq1yiRwS7t27Ro6d+5sUi4F/uXLl9G2bVur65lz8uRJ1KpVy+g1WFqTJk2MHlu7ja5du4Zhw4ZBCIHIyEi4urpi165dePPNN3H9+nWMHTsWERERaNKkCZYtW4aIiAh06tQJOp0OmzZt0r9/GjZsaLZd6enpePHFF9G4cWOMHTsWtWvXxpo1axAVFYXNmzejWbNmJvNU5LVr6b3Qt29f3LhxAxs2bMCECRPQrl07AMCyZcuwdOlSjBo1Cp6enrh48SLWrFmDEydOYMeOHUav2y5dumDZsmX48ccf0a1btzL/BzahevfgU2l//vmnUKlUYvLkyRWar2fPnuKFF16wWN6zZ0/h6ekpTp8+rS+LjY0VKpVKTJo0SV+Wl5cnvLy8jI4HlnX8s3S54ePr168LT09PMX/+fKN5zp49K1QqlZg9e7a+zNwx8BdeeEH07NlT/3jUqFGid+/eRnXmzJkjvL29xe3bt/Xz9OnTR+Tl5enrFBYWipEjR4pu3bqJ/Px8k+cgkY6Bl/fz3XffGbVPpVKJXbt26cvu3bsnOnfuLCIiIozqWdMmaXk///yzUbuioqJE586dxZ9//qkvO3nypGjTpo3RNps9e7Zo06aNuH79ur5s+/btZpdpaNOmTUKlUonjx4+XWaesY+CTJk0Sbdq0ESdOnNCXFRUVifHjxwtPT0/x+++/CyFK/r/vvfeevp5WqxWjRo0SPj4+4tatW2Wuu7Syzpfw8vISU6dONSk/ePCgUKlUYtu2bRWqZ07//v1FQECA1W0VomLbyN/fX1y7dk1fT6fTiSlTpghvb2+RlZUlhBDi8OHDQqVSic2bN+vrmdsmpd//Q4cOFQEBASI7O1tflpGRIdq0aSNiY2OFEKbvuYq+di29F8wdAw8NDRXjxo0zavv69evFoEGDxIULF4zKb9y4IVQqlVi4cKGwddyFXsNIParydnU+rObNmxtdguPh4QEA6Nu3r76sbt26aNCgAW7cuPFQ62rUqBGOHTuGV199VV8mhIBWqwVQfClPRYSFhUGj0eDEiRP6ZSUlJaFHjx547LHHcPPmTRw5cgTdu3fHvXv3kJ2djezsbNy5cwd9+/ZFVlYWfv31V4vrmTZtGlatWmX2p2PHjkZ169Spgx49eugfOzo6wsPDA1lZWQBQ4TbVrl1b3zMBgNu3b+PIkSMYNGgQ6tevry9v27YtAgICjNoycOBA6HQ67Ny5U1+WkJAANzc3tG/fvsznK12+Z64HVp6ioiKkpaUhMDDQqEeqVCoxYcIECCGwe/duo3kMe9p2dnaIiorCvXv3cOjQoQqtuyzlXZlhOM3aeqUplcoKvT+t3UY6nQ6pqanw8/ODvb29/nVy8+ZN9OvXDwUFBTh48KDV6y3tzz//xC+//IKwsDA8/vjj+nIPDw9s3rwZr7zyisk8FX3tWnovlOWJJ57ADz/8gNWrV+vrDh8+HN9++63J3riGDRuiTp06uHjx4oNshr8V7kKvYVxdXVGrVq0KXyNZEQ0aNDB6LO2eMgwHqVxUwt1mHRwcsG3bNhw4cADnz5/HhQsX9MFd0eWHhIRg9uzZSE5Ohre3N44dO4Zr165h4MCBAEqC6KuvvsJXX31ldhnmjquV5uXlZfUlLvXq1TPalQ0UP2fpuF9F21R6eRqNBjqdzuylTU8//bTRscgOHTqgefPmSE5OxgsvvICcnBzs378fY8aMKfc53Lp1C0DFT4C6efMm7t69q/8SaKhly5YAio/9Gj630rt3pedlWO9B1a1b1+wJaNLxUun5WVvPnEaNGiEjIwOFhYXlHuaSWLuNbt68iZycHKSmpiI1NdXssqx57ZZF2r7mXkdlHS542NcuYPxeKMtbb72FiRMnYs6cOfjggw/g5eWFXr164fnnnzd7voOzszNu3rxZ7jJtAQO8hlEoFPD19cWJEyeg1Wphb2/+XxQXFweNRgO1Wl3mCT2A+Z58Wct8kGvKLfVE8vPzMXLkSJw6dQpdunRB165d8eKLL8Lf39/om7q1XF1dERQUhOTkZMTExCAxMREuLi7o2bOnUXtGjRpV5tnUrVq1qvB6y1P6A6u0irap9HFqaW+FueO/jo6OJmUDBgzA8uXLcf36dRw4cACFhYX6LziWnkNFv1CVV1/60DZst7nAk+qVdXy+Ip588klcv37dpFwqk45NW1vPHF9fXxw8eBAnT55Ehw4dzNZJTU3F+vXr8eqrr5o9n0NiuI2k10lwcDCGDx9utr6bm1uZy7JEWldF3ucVfe1aei+UpU2bNti5cyf279+PPXv2YP/+/Vi0aBFWrVqFDRs26L/oSHQ6XaW8XuSOAV4D9e3bF0eOHEFCQoLZUb/u3buHTZs2oaioSD9QhFKpREFBgVE9rVaLmzdvlvsBUhHm1mFp11hSUhJOnDiB999/H0OHDtWXX7t27YHbERYWhsmTJ+PUqVNISUlBv3799CHx1FNPASgOg9InuPzxxx+4ePFihQe2eFgP2ybpQ/v8+fMm0y5cuGBSFhYWhk8//RRpaWnYu3cvPD090bp163LbKO2VuXXrVrnhVVr9+vVRt25dZGRkmEw7d+4cABidiZ2VlYW8vDyjy5+k52Vp8BRrtG3bFsnJySZffn/77Tc4Ojrqt4O19czp27cvlixZgvj4+DIDfNOmTThw4ACmTp1q9TaqX78+6tSpA61Wa/I6uXz5Mn777beHeu0++eSTAIDMzEyTaR9++CFcXV0xbtw4o/JH8X4qKirC6dOn4ezsjN69e6N3794AgMTEREyePBnx8fEmV7ncvn3bZE+iLeIx8BooIiICTz31FObNm4fff//daFpRURFmzpyJrKwsvPLKK/oeTcOGDXHu3DmjSyt2795d7vWsFdWwYUOcPn3aqNeVmJhY7jzSrtnSvd41a9YAKOldAiXf3i3tbuvVqxecnJywcOFC3LhxA2FhYfppjRs3hre3N7Zs2WL0JaGwsBAzZszA66+/brTOR+Fh29SgQQP4+vpix44duH37tr5co9Fg3759JvVbtmyJtm3bIjU1Fd9//73F3jdQ8kFd0V20dnZ2CAoK0vdIJUIIfP7551AoFEZ7WoQQRpe5abVarF69Gi4uLujatWuF1m1OcHAw7t69q7+kCgCys7ORnJyM4OBgfVhbW8+cNm3aoH///tiyZYv+kjJD8fHx2LNnD3r06IG2bdtavY3s7e3x7LPPYu/evTh9+rTRMufOnYvo6OiH2m3cpEkTtGnTBgkJCUZDvWo0GqxZs8bsl/GqeD+Vfp8XFRUhKioKc+bMMaonnbNRuld/48YNaLVa/RcSW8YeeA3k6OiIJUuW4OWXX8bQoUMRFhaGdu3a4datW0hOTsapU6cQEhKCl156ST/PwIEDMXv2bIwdOxaDBg3ChQsXsHHjRv0Hc2UYOHAgVq5ciddeew09evTAyZMnkZSUZHLs3FC3bt1gb2+Pt956C6NGjYK9vT327NmDAwcOoFatWkYnsUnLWbRokX53uzm1a9dGv379sGXLFjRu3NjkWPXbb7+N0aNHY8iQIRgxYgTq1auHhIQEHD9+HFOnTjU6gacshw4dKncQjr59+6Ju3boWl1NZbZo2bRoiIyMxdOhQDB8+HAUFBfjqq6/K/LIzcOBAzJs3DwqFoszxBAw988wzAIovySurV1mWmJgY/PDDD4iMjERkZCQaNWqE7777DocPH8ZLL71k8uXtk08+waVLl9C6dWskJSXhp59+wvvvv18pe0Z69OiBLl26YNasWdBoNGjSpAnWrl0LIQRee+21Ctcry7vvvguNRoMpU6bg22+/1Z9M+P3332PPnj1o2bIl3n///QpvI6neqFGjMGrUKDRt2hRpaWnYs2cPIiIiLO5JsUStVmPs2LEYMmQIhg0bBqVSibVr1+Kxxx4zexIbUDnvJ0PS+3z9+vXIyspCWFgYIiMj8emnnyI6OhpBQUG4d+8eNmzYgDp16mDIkCFG8x8/fhwAKuULn9wxwGuotm3b4ttvv8WXX36Jffv2ITExEUIIeHp6Ys6cOQgPDzc6ljVy5EjcunULmzZtwuzZs9GmTRssWbIEK1eurLTB/9944w1otVokJCTgwIEDaN++PVavXo2YmJgy51GpVFi0aBGWLFmCBQsWwMnJCa1bt8aqVauwbt06HDlyRH8y0IgRI3D48GGsWLECv/76a7lv0LCwMGzZsgUDBgww+Ybu6+uL9evXY/HixVi1ahW0Wi08PDwwd+5cq69pX7ZsWbnTd+3aVaEAf9g2+fr6YsWKFYiLi8PHH3+MevXqITIyEmfPnjU641wycOBAzJ8/H+3bt7fqS9wTTzwBlUqFY8eOYfTo0VY/L6D4qoaNGzfi448/xjfffIN79+7pA8zwsInkiy++wMyZM7Flyxb92NaGV0A8rKVLl2L+/PnYsGEDCgsL4ePjgwULFpjsore2njn16tXDV199hfj4eGzbtg2ffPIJ8vLy4ObmhujoaIwZM8boMIG120iqt2jRImzcuBF3796Fm5sb1Gp1hW9KYs4zzzyD1atXY9GiRVi6dCkcHR3RuXNnvPnmm2WeS1MZ7ydDXbt2RWhoKPbs2YPDhw+jX79+eP3111GvXj1s3rwZsbGxsLOzQ8eOHfHhhx+aHP8+duwYXF1dK/xF8+9IISrjNGMiqlJZWVlmB+eYMGECTp8+jbS0NKPy69evo3v37vj3v/+NkSNHWrWO1atXY8GCBTh48GCVDMc5ffp0bNmyxWjUL6KK0Ol06NmzJ0JCQvQjutkyHgMnkoFhw4aZXAqWlZWFH374AT4+Pib1N27cCAcHB6t2n0uGDh0KR0dHJCUlPXR7iarCDz/8gKysrArvJfq74i50IhkYNGgQli1bhqlTp6JLly64c+cONm7cCJ1Oh+joaH29jz76CP/73/+wd+9ejBo1Cq6urlavw8nJCRMnTsQXX3yB8PBwXqZDNc7y5csxYsQING3atLqbUiMwwIlk4I033kDDhg2xceNG7Nq1C46OjujYsSMWLVpkNKre3bt3cfjwYfTp0wdTpkyp8HqioqKwfft2bN++3eydtoiqy5EjR3Du3DksWbKkuptSY/AYOBERkQzJpgfetm1b6HS6v/W9bomIiCS5ublQKpX47bffzE6XzUlsOp2uUsblJiIikgMhRLkDW8mmBy71vNPT06u5JURERFXPz8+v3Omy6YETERFRiQoH+KlTp+Dl5VXuMJNA8X2eZ82ahYCAAPj6+uKVV14xezMGIiIiqrgKBfjZs2cxfvx4qwavnzx5sv6Wj7Gxsbh27RqioqKQk5PzwI0lIiKiYlYFuFarxddff41hw4ZZdXer9PR07N27F7GxsXjuuefQr18/fPnll8jJycH69esfutFERES2zqoAP3bsGObPn4+XX3653BtXSA4ePAgnJyf9HXqA4jvQdO7c2eztD4mIiKhirArwli1bIjU1Fa+99ppVwytmZGTA3d3dpG7z5s31N7AnIiKiB2fVZWTm7oJUntzcXLMDrjg5ORndSJ6IqCYRQkAIQEh/A/cfi/vTof8tYL4uStUva1nQT5MeCoO/S9pSen1G01EyHSbrNF6uufUYLrt0e0vWY7zskkf365eqa/j8SteV1m+4Lv121zfItF7pdZs+97LbWNa6jdtmPMaI2f+JMK1fus32dgr09GyMx50c8ChUyXXg5Q24UvrezUSVRQiBIp1AkRDQ6QCtTqf/LZUV/y6pVzwPUKQT0BnNL+6XQV+uE1Kd4jKdwfSypknr0OkEiu5/2BXXLa4nBIyWI4ThMo3ri3LqiNKPYbB8fXiYqWcwXSor3pYl8wsU10fpZd2fB6XCSaczDoGS5Rh88JZ6rK8jSv6XOpOQLJ7RbCiWsdzSQQuTZRoEL1ElCPZqguWR5V+/XVmqJMCdnZ1x8eJFk/K8vDwOhVoDCSFQUKRDvlaHAq3h7yIUaHUovD+tsEigUKtDQVFxWYH0t1YHrU7c/1ugsEiHQl3x31pd8XzaIh2KdAKFOoGi+2VFuuK6RToB7f3A1N6frr0/3XiazqRMClatTvBDmKiGUSju/9Y/VpR6DChQUklhMp+izOUoDApMll+qvKx1G9YyXUdJGwynKRSGSzWopwBqKZXo1/YJPCpVEuAeHh74/vvvIYQwerIXLlyAh4dHVazyb0sIgb8Ki5Cbr0VefhHy8rXFPwXFj+8a/i4owl/3f+5ppd865BcWIf9+MOdri5BfKIV0STmVTaEA7BQKKJUKKKW/7z+2U97/WwH93woFoFQUT5P+Vup/my9XGExXKos/NKT1KVB+neLpxR8sJfWk6QBwfzkG80sfnIblUJR8WCkN/pY+LKV60t8Kg+Xcn71kWSiZT2lUt7iSfhkGbVBAakPJ8zVcvuEyza0XRo8N22m87LKWAZNlmi4f+vnMLwcWHpe3LKnUbHCU+j8YTjfMIYWF5RqFpcG8+kA0036Umteo3Ny8pQKOqk6VBHhgYCCWLVuGQ4cO6c9Ez87ORnp6OsaPH18Vq6zxhBDIKyhCdm4B/szLR3ZeAbLzCnDrbiFu/SX9LsRtg8e3/ypEbr62xvQsa9kpUMtOCQd7JRzslEZ/29+fJtWxt1PCwU4Be2XJNDulArXsigPMXqm8/7cS9koF7O0UsFcWP7ZTQj+fnVIBO4UC9nbF5XZKJewU0jLuB+j9v6VwLF1PKpOmKw3KpXCWlmNXqo6SH0hEVENVSoBnZ2cjMzMTrVq1grOzMzp37gx/f39MmTIFMTExqFevHhYvXgwXFxeMGDGiMlZZoxXpBHb8chkJv1zBtZx8ZOXkIys3v1J7uo72Sjg72qOOgx2cHe1R18EOTo72qF3LDnUd7FCnlh1q17KDYy0latsX/3awU8LRXgnHWnbFv+2l30o41ip+7HD/sYN9yWOH+8HMICMiqjkqJcDT0tKgVquxZs0adOnSBQCwZMkSzJ07F/PmzYNOp0OnTp3w8ccfw9XVtTJWWSNpi3T49ufLWLrnD2Rk5ZVb18FeiQZODqhX1wGP162FenVrwbWO4d/FPy61a8HJ0R7OjsUBXdfBHk4OdrC348mARES2TCFkco9O6a4sNfFuZEIIfPvzZSz47ndkZt/Vl/dt2wQd3OqhobMDGjo7ooGzIxo4OaC+kwPqOtixR0tERGWylHuyuZ1oTXXi0m3M3HYS6RduAig+iWOgT1O81rMVPJ9wqebWERHR3xUD/AHduluAD3eewbojmfqTzEK9n0BMsCdaNuKlckREVLUY4A/gXmERhn92GKevFt9ZrXVjZ8wc5IWAVhUbsY6IiOhBMcAfwNyk0zh9NQd2SgVm9P8Horq6oxZPKiMiokeIAV5Be3+/gS8PnQcATOrdGmMCOTANERE9euw2VkB2XgFi4o8DAPzcH8erPVtVc4uIiMhWMcCtJITA9M2/4EZOPpwd7REX0QF2Sl4GRkRE1YMBbqWN6Rqk/HYNADBrkBfc6tet5hYREZEtY4Bb4WZeAd7bcQoAMKDdkwjv+FQ1t4iIiGwdA9wKn+49i5x8LVxq2+O9wd4cQY2IiKodA9yCq7fvYfX9s84ndG+Jx50cqrdBREREYIBbtGj3/5Cv1aGhsyNeCmhR3c0hIiICwAAv1/msPGw8qgEA/KtXK9R14GXzRERUMzDAyxGX+ju0OoGn6tXBcH+36m4OERGRHgO8DKeu3MG245cBAJP7quBob1fNLSIiIirBAC/DRylnIATQqrEznvPlZWNERFSzMMDNuHzrL6Seug4AmNpXxRHXiIioxmGAm3Hgf1kAAJfa9ujbtkk1t4aIiMgUA9yMff+7AQDo1rIB7HmbUCIiqoGYTqXodAKHzv4JAAhs3aiaW0NERGQeA7yU367cQXZeAQDg2dYNq7k1RERE5jHAS5F2n7vVrwP3Bk7V3BoiIiLzGOClSCewBXH3ORER1WAMcAN/FRQh/fxNAEBQK+4+JyKimosBbuCHc3+ioEgHpQLo1pIBTkRENRcD3IC0+9ynWT241q1Vza0hIiIqGwPcwIE/pOPf7H0TEVHNZnWA79ixAwMGDICPjw9CQ0OxdevWcutnZ2dDrVYjMDAQ/v7+GD9+PM6fP/+Qza061+/cw+mrOQCAQB7/JiKiGs6qAE9MTERMTAwCAwOxdOlS+Pv7Y9q0aUhOTjZbXwiB6Oho7Nu3DzExMZg3bx5u3LiBqKgo3L59u1KfQGWRet9ODnbwbf54NbeGiIiofPbWVIqLi0NoaCjUajUAICgoCLdv38bChQsREhJiUv/8+fP48ccfERsbi8GDBwMAWrZsiT59+mD37t147rnnKu8ZVBLp+PczTzeAgz2PLBARUc1mMak0Gg0yMzPRr18/o/Lg4GBkZGRAo9GYzJOfnw8AcHIqGQjF1dUVAHDr1q2HaW+VEEJg//0eeCCPfxMRkQxYDPCMjAwAgIeHh1G5u7s7AODcuXMm87Rp0wZdunTB0qVLcfbsWWRnZ+O9995D3bp10adPn8pod6XSZP+FGznFXzoCePybiIhkwOIu9Jyc4hO7nJ2djcql3nVubq7Z+WbOnImxY8eif//+AAAHBwcsXboUbm5uD9XgqpCTX6j/u4lL7WpsCRERkXUs9sCFEOUvQGm6iLNnzyIiIgKPP/44li5dii+++AI9e/bE66+/jvT09AdvbRUpLCp5jjz+TUREcmCxB+7i4gIAyMvLMyqXet7SdENffvklAGDlypX6Y98BAQEYOXIk5syZg//+978P1ejKVqDV6f+uZaeoxpYQERFZx2J3Uzr2nZmZaVR+4cIFo+mGLl++jJYtW+rDGwAUCgU6deqEP/7446EaXBWkAFcqAHs79sCJiKjms5hW7u7uaNasmck13ykpKWjRogWaNm1qMo+Hhwf+97//4c6dO0blx48fx1NPPfWQTa58hUXFAV6L4U1ERDJh1XXg0dHRUKvVcHV1RY8ePbBr1y4kJSUhLi4OQPGoa5mZmWjVqhWcnZ3x4osvYtu2bXj55Zcxbtw41K5dG99++y2OHDmin6cmyb/fA+fxbyIikgurAjw8PBwFBQVYuXIl4uPj4ebmhtjYWP0Z5mlpaVCr1VizZg26dOmCZs2aYf369Zg/fz7UajUUCgVUKhVWrVqFbt26VekTehBSD9yBPXAiIpIJhbB0mnkN4efnBwBVchb75mMXMTX+OJ50rY3v1b0rfflEREQVZSn32OWEQQ+cu9CJiEgmmFgACngSGxERyQwTCyWXkfEYOBERyQUTCwY9cO5CJyIimWBioaQH7sgeOBERyQQTCwYDudhzGFUiIpIHBjh4DJyIiOSHiYWSu5HxLHQiIpILJhY4lCoREckPEwscyIWIiOSHiQUeAyciIvlhYsEgwNkDJyIimWBigfcDJyIi+WFioWQkNvbAiYhILphYKNmFzh44ERHJBRMLJT1wR/bAiYhIJphYMLiMjD1wIiKSCSYWDHehcyx0IiKSBwY4DC8js6vmlhAREVmHAQ7DsdDZAyciInlggINjoRMRkfwwscCT2IiISH6YWOBQqkREJD9MLHAoVSIikh+bTyydTkCrKz6JjT1wIiKSC5tPLGkUNoABTkRE8mF1Yu3YsQMDBgyAj48PQkNDsXXr1nLr63Q6fPrpp+jduzd8fHwQFhaGhISEh21vpTMKcO5CJyIimbC3plJiYiJiYmIwevRoBAYGIjU1FdOmTUPt2rUREhJidp45c+Zgw4YNmDJlCtq0aYOEhARMnToVzs7O6N69e6U+iYchncAGsAdORETyYVWAx8XFITQ0FGq1GgAQFBSE27dvY+HChWYDPDMzE19//TX+85//YNiwYQCArl274vz589i/f3+NCvBCgx44T2IjIiK5sBjgGo0GmZmZmDJlilF5cHAwkpKSoNFo4ObmZjQtNTUVtWvXxuDBg43K165d+/AtrmTsgRMRkRxZTKyMjAwAgIeHh1G5u7s7AODcuXMm85w5cwYeHh44dOgQBg0ahLZt26Jfv35ITEysjDZXKuMeOIdSJSIiebAY4Dk5OQAAZ2dno3InJycAQG5ursk82dnZuHLlCmbMmIEXXngBK1asgJeXFyZPnozDhw9XRrsrTb5BD9zRjjczISIiebC4C10IUe50pdL0O0BhYSGys7OxbNky9OzZEwDwzDPPICMjA0uWLMEzzzzzgM2tfNKNTADuQiciIvmwmFguLi4AgLy8PKNyqectTTfk5OQEOzs7BAQElKxIqUS3bt1w5syZh2pwZTM8Bs5d6EREJBcWA1w69p2ZmWlUfuHCBaPphtzd3aHT6aDVao3KCwsLoVDUrJCUAlypAOx5FjoREcmExcRyd3dHs2bNkJycbFSekpKCFi1aoGnTpibzBAUFQQiBpKQkfZlWq8X+/fvRqVOnSmh25eE46EREJEdWXQceHR0NtVoNV1dX9OjRA7t27UJSUhLi4uIAFJ+0lpmZiVatWsHZ2Rldu3ZF9+7d8d577+Hu3bto0aIF1q1bh0uXLuGjjz6q0idUUbwXOBERyZFVAR4eHo6CggKsXLkS8fHxcHNzQ2xsLPr37w8ASEtLg1qtxpo1a9ClSxcAwKJFi7Bw4UJ89tlnuH37Ntq2bYuVK1fC29u76p7NA+C9wImISI4UwtJp5jWEn58fACA9Pb1Sl7v52EVMjT+OJ11r43t170pdNhER0YOylHs23+3U98C5C52IiGTE5lOrgCexERGRDNl8akmXkfEYOBERyYnNp5a+B85d6EREJCM2n1pSD9yRPXAiIpIRm08t/UAu9jVrhDgiIqLy2HyA8xg4ERHJkc2nlnQ3Mp6FTkREcmLzqcWhVImISI5sPrU4kAsREcmRzacWj4ETEZEc2XxqFXAXOhERyZDNpxbvB05ERHJk86lVwGPgREQkQzafWtIudPbAiYhITmw+taQeuCN74EREJCM2n1r6y8jYAyciIhmx+dQq2YXOsdCJiEg+GOD6y8jsqrklRERE1rP5AC8ZC509cCIikg+bD3COhU5ERHJk86nFk9iIiEiObD61OJQqERHJkc2nFodSJSIiObLp1NLpBLS64pPY2AMnIiI5senUkkZhAxjgREQkLzadWkYBzl3oREQkI1an1o4dOzBgwAD4+PggNDQUW7dutXolV65cQadOnfDJJ588SBurjHQCG8AeOBERyYtVqZWYmIiYmBgEBgZi6dKl8Pf3x7Rp05CcnGxxXiEEZsyYgdzc3IdubGUrNOiB8yQ2IiKSE3trKsXFxSE0NBRqtRoAEBQUhNu3b2PhwoUICQkpd95169YhIyPj4VtaBdgDJyIiubKYWhqNBpmZmejXr59ReXBwMDIyMqDRaMqdd/78+Zg9e/bDt7QKGPfAOZQqERHJh8UAl3rPHh4eRuXu7u4AgHPnzpmdT6fTYfr06QgNDcWzzz77sO2sEvkGPXBHO97MhIiI5MPiLvScnBwAgLOzs1G5k5MTAJR5bHv16tW4ePEili1b9rBtrDLSjUwA7kInIiJ5sRjgQohypyuVpsF39uxZfPzxx1i0aBFcXFwevHVVzPAYOHehExGRnFjsdkoBnJeXZ1Qu9bxLB3RRURHUajVCQkIQEBAArVYLrVYLoHi3uvR3TSAFuFIB2PMsdCIikhGLqSUd+87MzDQqv3DhgtF0yZUrV3D8+HFs3boVXl5e+h8AWLx4sf7vmoDjoBMRkVxZ3IXu7u6OZs2aITk5GX379tWXp6SkoEWLFmjatKlR/caNG2PTpk0myxk6dChGjBiBIUOGVEKzKwfvBU5ERHJl1XXg0dHRUKvVcHV1RY8ePbBr1y4kJSUhLi4OAJCdnY3MzEy0atUKzs7OaNeundnlNG7cuMxp1YH3AiciIrmyKrnCw8Mxa9YsHDhwANHR0Th69ChiY2PRv39/AEBaWhoiIiJw8uTJKm1sZeO9wImISK6s6oEDwPDhwzF8+HCz08LDwxEeHl7u/GfOnKlYyx4BHgMnIiK5sunkku5Gxh44ERHJjU0nl34XOnvgREQkMzadXFIPvBZ74EREJDM2nVxSD9yRPXAiIpIZm04u/Uls9hxGlYiI5MWmA5zHwImISK5sOrmku5HxMjIiIpIbm04uDqVKRERyZdPJVcjrwImISKZsOrl4DJyIiOTKppOLY6ETEZFc2XRycSx0IiKSK5tOLo6FTkREcmXTySXtQmcPnIiI5Mamk0vqgTuyB05ERDJj08lVcgycQ6kSEZG82HSA8zIyIiKSK5tOrpLLyOyquSVEREQVY9MBXjIWOnehExGRvNh0gHMsdCIikiubTi79WOg8Bk5ERDJj08nFoVSJiEiubDq5OJQqERHJlc0ml04noNUVn8TGHjgREcmNzSaXNAobwAAnIiL5sdnkMgpw7kInIiKZsdnkkk5gA9gDJyIi+bE6uXbs2IEBAwbAx8cHoaGh2Lp1a7n1b9y4gbfffhs9e/aEr68vwsPDkZSU9LDtrTSFBj1wnsRGRERyY29NpcTERMTExGD06NEIDAxEamoqpk2bhtq1ayMkJMSkfkFBAcaOHYucnBy8/vrraNy4MXbu3IlJkyahqKgIAwcOrPQnUlHsgRMRkZxZFeBxcXEIDQ2FWq0GAAQFBeH27dtYuHCh2QDft28fTp8+jfj4ePj4+AAAAgICcPnyZXz++ec1IsCNe+AcSpWIiOTFYtdTo9EgMzMT/fr1MyoPDg5GRkYGNBqNyTxOTk6IiIhAu3btjMqffvppZGZmPmSTK0e+QQ/c0Y43MyEiInmx2APPyMgAAHh4eBiVu7u7AwDOnTsHNzc3o2ldu3ZF165djcoKCwuxd+9etG7d+qEaXFmkG5kAQC179sCJiEheLPbAc3JyAADOzs5G5U5OTgCA3Nxcq1b04Ycf4vz58xg3blxF21gljI6B8yQ2IiKSGYs9cCFEudOVyvLDTwiBDz/8EKtXr8aYMWPQp0+firWwikgBrlQA9gxwIiKSGYsB7uLiAgDIy8szKpd63tJ0cwoKCjB9+nQkJCRgzJgxeOuttx6mrZWK46ATEZGcWQxw6dh3ZmYmPD099eUXLlwwml5abm4uxo8fjx9//BEzZszA6NGjK6O9lYb3AiciIjmzmF7u7u5o1qwZkpOTjcpTUlLQokULNG3a1GSeoqIiTJw4EcePH0dcXFyNC2+A9wInIiJ5s+o68OjoaKjVari6uqJHjx7YtWsXkpKSEBcXBwDIzs5GZmYmWrVqBWdnZ3zzzTc4cuQIIiIi8MQTT+Dnn3/WL0uhUKB9+/ZV8mQqgvcCJyIiObMqwMPDw1FQUICVK1ciPj4ebm5uiI2NRf/+/QEAaWlpUKvVWLNmDbp06YKdO3cCADZs2IANGzYYLcvOzg6//fZbJT+NiuMxcCIikjOrAhwAhg8fjuHDh5udFh4ejvDwcP3jNWvWPHzLqph0NzL2wImISI5sNr30u9DZAyciIhmy2fSSeuC12AMnIiIZstn0knrgjuyBExGRDNlseulPYuM46EREJEM2G+A8Bk5ERHJms+kl3Y2Ml5EREZEc2Wx6cShVIiKSM5tNLw6lSkREcmaz6cWhVImISM5sNr0Y4EREJGc2m14cC52IiOTMZtOLY6ETEZGc2Wx6SbvQ2QMnIiI5stn0knrgjuyBExGRDNlsepUcA+dQqkREJD82G+AcSpWIiOTMZtOr5DIyu2puCRERUcXZbICXjIXOXehERCQ/NhvgHAudiIjkzGbTi2OhExGRnNlsenEoVSIikjObTS8OpUpERHJmk+ml0wlodcUnsbEHTkREcmST6SWNwgYwwImISJ5sMr2MApy70ImISIZsMr2kE9gA9sCJiEiebDK9Cg164DyJjYiI5Mjq9NqxYwcGDBgAHx8fhIaGYuvWreXWz8vLw6xZsxAQEABfX1+88sorOH/+/EM2t3KwB05ERHJnVXolJiYiJiYGgYGBWLp0Kfz9/TFt2jQkJyeXOc/kyZORnJyMmJgYxMbG4tq1a4iKikJOTk6lNf5BGffAOZQqERHJj701leLi4hAaGgq1Wg0ACAoKwu3bt7Fw4UKEhISY1E9PT8fevXvx+eef49lnnwUA+Pn5oXfv3li/fj3GjRtXiU+h4vINeuCOdryZCRERyY/FHrhGo0FmZib69etnVB4cHIyMjAxoNBqTeQ4ePAgnJycEBAToy+rXr4/OnTtj3759ldDshyPdyAQAatmzB05ERPJjMcAzMjIAAB4eHkbl7u7uAIBz586Zncfd3R12pXq3zZs3N1v/UTM6Bs6T2IiISIYsppd0zNrZ2dmo3MnJCQCQm5trMk9ubq5JfWkec/UfNSfH4i8W9erWgj0DnIiIZMjiMXAhRLnTlUrTACxvHnP1H7W2Tz6GeUN90KKBU3U3hYiI6IFYDHAXFxcAxZeFGZJ60tJ0Q87Ozrh48aJJeV5entme+aOmUCjwvJ9bdTeDiIjogVnsDkvHvjMzM43KL1y4YDS99DwajcakJ37hwgWz9YmIiKhiLAa4u7s7mjVrZnLNd0pKClq0aIGmTZuazBMYGIg7d+7g0KFD+rLs7Gykp6ejW7duldBsIiIi22bVdeDR0dFQq9VwdXVFjx49sGvXLiQlJSEuLg5AcThnZmaiVatWcHZ2RufOneHv748pU6YgJiYG9erVw+LFi+Hi4oIRI0ZU6RMiIiKyBVYFeHh4OAoKCrBy5UrEx8fDzc0NsbGx6N+/PwAgLS0NarUaa9asQZcuXQAAS5Yswdy5czFv3jzodDp06tQJH3/8MVxdXavu2RAREdkIhbB0mnkN4efnB6B4lDciIqK/O0u5V/3XdBEREVGFWbULvSbIzc2FEEL/jYSIiOjvLCcnBwpF2cN9y6YHrlQqy30iREREfycKhaLcwc9kcwyciIiISsimB05EREQlGOBEREQyxAAnIiKSIQY4ERGRDDHAiYiIZIgBTkREJEMMcCIiIhligBMREckQA5yIiEiGGOBEREQyxAAnIiKSIZsM8B07dmDAgAHw8fFBaGgotm7dWt1NqtF0Oh3Wr1+PsLAw+Pr6ok+fPvjggw+Qm5urr/Prr78iMjISvr6+CAwMxIIFC1BYWFiNra75XnvtNfTt29eo7MCBAxgyZAjat2+PXr16YeXKldXUuprt6NGjGDFiBNq3b4/AwEDMnj0beXl5+uncjtZZv349QkND0aFDB4SFhWHbtm1G07kdy3fq1Cl4eXnh6tWrRuXWbLdK+cwUNiYhIUF4enqKOXPmiH379ol33nlHqFQqkZSUVN1Nq7GWL18u/vGPf4j58+eLgwcPirVr1wp/f3/x8ssvCyGEOH/+vOjYsaMYM2aMSEtLE1988YXw9vYWs2bNquaW11xbt24VKpVK9OnTR1927Ngx4eXlJWJiYsTevXvFggULhKenp1ixYkU1trTm+emnn4SXl5f417/+JQ4ePCjWrVsn/Pz8xKRJk4QQ3I7W+uabb4RKpRJz584VBw8eFHPmzBEqlUokJiYKIbgdLfnjjz9EUFCQUKlU4sqVK/pya7ZbZX1m2lyA9+nTR/9Gl7zxxhsiJCSkmlpUs+l0OtG5c2cxc+ZMo/KEhAShUqnEb7/9JmbMmCG6d+8u8vPz9dO//vpr8Y9//ENcvXr1UTe5xrt69aro3LmzePbZZ40CfPTo0WLYsGFGdefNmyf8/PyMtq2tGzVqlBg1apTQ6XT6srVr14revXuLu3fvcjtaKSIiQkRGRhqVjRw5UrzwwgtCCL4ey1JYWCjWrl0rfH19hb+/v0mAW7PdKusz06Z2oWs0GmRmZqJfv35G5cHBwcjIyIBGo6mmltVceXl5GDRoEAYOHGhU/vTTTwMAMjMzcfDgQfTs2RMODg766SEhISgqKsKBAwceaXvl4O2330ZAQAC6du2qL8vPz0d6errZ1+adO3fw448/Pupm1kjZ2dlIT0/HiBEjoFAo9OWjRo1CamoqlEolt6OV8vPz4eTkZFRWr1493Lp1i6/Hchw7dgzz58/Hyy+/jJiYGKNp1m63yvrMtKkAz8jIAAB4eHgYlbu7uwMAzp0798jbVNM5Ozvj7bffRqdOnYzKU1NTAQAtW7bElStXTLZp/fr14ezszG1aSnx8PE6ePIl///vfRuUajQaFhYV8bVrw+++/QwgBV1dXTJo0CR06dECnTp3w7rvv4t69e9yOFRAVFYX9+/cjKSkJubm5SE5ORlpaGv75z39yO5ajZcuWSE1NxWuvvQY7OzujadZst7/++qvSPjPtH/A5yFJOTg6A4lAyJH0LNTwpi8p2/PhxfPbZZ+jTpw8ee+wxAKbbFCjertymJS5duoQPPvgAH3zwAerXr280ja9N62RnZwMApk+fjr59++LTTz/FmTNn8PHHHyM/Px8REREAuB2tMWDAABw+fBiTJk3Slz333HMYO3YsfvrpJwDcjuY0bNiwzGnWvI/LqiPVq8i2takAF0KUO12ptKkdEg/k2LFjmDBhApo1a4b33nsPBQUF5dbnNi0mhMCMGTPQvXt3BAcHm51eHm7HYtJZuh07dsS7774LAOjatSuEEIiNjcXzzz9f7vzcjiUmTpyIn376CWq1Gm3btsXx48fxySefwNnZGf379y93Xm5H86x5H1fme92mAtzFxQUAjC43AUq+TUrTybzExERMnz4dLVq0wIoVK/D444/rt2XpbQoUb1du02Jff/01zpw5g+3bt0Or1QIoebNrtVq+Nq0k9WSeffZZo/LAwEDMnTsXv/76KwBuR0t+/PFHHDhwAB988AHCw8MBAP7+/njsscfwzjvvYOjQoQC4HSvKmvex1POujM9Mmwpw6ZhDZmYmPD099eUXLlwwmk6mVq1ahdjYWPj7+2Pp0qX6F5mTkxOaNGmi34aSP//8E3l5edym9+3cuRM3b95EYGCgyTQvLy/MnDkTdnZ2yMzMNJomPeZ2LNaiRQsAMNnzI/XMmzVrxu1ohcuXLwMo3pNhyM/PD0Dx9c3cjhXXvHlzi9utMj8zbWo/iLu7O5o1a4bk5GSj8pSUFLRo0QJNmzatppbVbPHx8Zg7dy5CQ0OxYsUKk2+IAQEB2LNnj9GH6s6dO2FnZwd/f/9H3dwaadasWdi0aZPRT8+ePfHEE09g06ZNCAkJgZ+fH1JSUox2se3cuRMuLi7w9vauxtbXHC1btsRTTz2FxMREo/I9e/bA3t4evr6+3I5WkELi2LFjRuU///wzgOKrTLgdK87R0dGq7VZpn5kPcBmcrG3evFmoVCoxa9YssXfvXv1ALgkJCdXdtBopKytLtG/fXvTs2VMcPXpU/PTTT0Y/f/75p/jjjz9Eu3btxOjRo8Xu3bvFypUrhbe3t3j33Xeru/k12rRp04yuAz906JDw9PQUb7zxhkhLSxNxcXHC09NTfPbZZ9XYyppHGoxp6tSp4uDBg2L58uXCy8tLzJkzRwjB7WitV199Vfj6+opVq1aJw4cPixUrVghfX18xduxYIQS3ozWkPDG8Dtya7VZZn5k2F+BCCLF+/XrRt29f4e3tLUJDQ8WWLVuqu0k11pYtW4RKpSrzZ+vWrUIIIY4ePSqGDRsmvL29RVBQkPjoo49EQUFBNbe+Zisd4EIIkZKSIgYOHCi8vLxEr169xBdffFFNravZvvvuOzF48GDh7e0tunfvLj755BNRVFSkn87taFl+fr746KOPRPfu3YW3t7cIDg4WixcvNhpchNuxfOYCXAjrtltlfGYqhLBwShwRERHVODZ1DJyIiOjvggFOREQkQwxwIiIiGWKAExERyRADnIiISIYY4ERERDLEACciIpIhBjgREZEMMcCJiIhk6P8BQ/miVWn8hhEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "top_k = 100\n",
    "ax.plot(cumulative_energy[0:top_k])\n",
    "ax.set_title(f'Cumulative Energy (Top {top_k} Coefficients)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's train a network using EN-regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Net(nn.Module):\n",
    "    def __init__(self, n, multiplier):\n",
    "        super(Net, self).__init__()\n",
    "        self.fc1 = nn.Linear(n, multiplier*n)\n",
    "        self.bn1 = nn.BatchNorm1d(multiplier*n)\n",
    "        self.fc2 = nn.Linear(multiplier*n, multiplier*n)\n",
    "        self.bn2 = nn.BatchNorm1d(multiplier*n)\n",
    "        self.fc3 = nn.Linear(multiplier*n, n)\n",
    "        self.bn3 = nn.BatchNorm1d(n)\n",
    "        self.fc4 = nn.Linear(n, 1)\n",
    "        torch.nn.init.xavier_uniform_(self.fc1.weight)\n",
    "        torch.nn.init.xavier_uniform_(self.fc2.weight)\n",
    "        torch.nn.init.xavier_uniform_(self.fc3.weight)\n",
    "        torch.nn.init.xavier_uniform_(self.fc4.weight)\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.bn1(F.leaky_relu(self.fc1(x)))\n",
    "        x = self.bn2(F.leaky_relu(self.fc2(x)))\n",
    "        x = self.bn3(F.leaky_relu(self.fc3(x)))\n",
    "        x = self.fc4(x)\n",
    "        return x.reshape(-1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_wht_spectrum(model, H, all_inputs, device=\"cuda\"):\n",
    "    '''\n",
    "    Compute the WH spectrum of the DNN model\n",
    "    Used for EN-regularization\n",
    "    '''\n",
    "    all_inputs = all_inputs.to(device)\n",
    "    H = H.to(device)\n",
    "    # get all outputs of model (aka landscape)\n",
    "    landscape = model(all_inputs).reshape(-1)\n",
    "    # compute spectrum by matrix-vector multiplication with WHT basis\n",
    "    spectrum = torch.matmul(H, landscape)\n",
    "    # clear up gpu memory\n",
    "    H.to(\"cpu\")\n",
    "    all_inputs.to(\"cpu\")\n",
    "    return spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(model, train_dl, val_dl, H, all_inputs, lr=1e-3, weight_decay=0.0, reg_lambda=1e-2, rho=0.01, num_epochs=20, device=\"cuda\"):\n",
    "    model.to(device)\n",
    "    optim = torch.optim.Adam(model.parameters(), lr=lr, weight_decay=weight_decay)\n",
    "    \n",
    "    l2_loss = torch.nn.MSELoss()\n",
    "    \n",
    "    # these are the uniform sampling protein pattern\n",
    "    # parameters are set as N=13 mutations, m=4 anticipated sparsity level, and d=3 number of offsets\n",
    "    sampling_matrix1  = pickle.load(open('N13-m4-d3/sampling-matrix-1.p','rb'))\n",
    "    sampling_matrix2  = pickle.load(open('N13-m4-d3/sampling-matrix-2.p','rb'))\n",
    "    sampling_matrix3  = pickle.load(open('N13-m4-d3/sampling-matrix-3.p','rb'))\n",
    "    delays_matrix1 = pickle.load(open('N13-m4-d3/delays-1.p','rb'))\n",
    "    delays_matrix2 = pickle.load(open('N13-m4-d3/delays-2.p','rb'))\n",
    "    delays_matrix3 = pickle.load(open('N13-m4-d3/delays-3.p','rb'))\n",
    "    all_sampling_locations1 = pickle.load(open('N13-m4-d3/sampling-locations-1.p','rb'))\n",
    "    all_sampling_locations2 = pickle.load(open('N13-m4-d3/sampling-locations-2.p','rb'))\n",
    "    all_sampling_locations3 = pickle.load(open('N13-m4-d3/sampling-locations-3.p','rb'))\n",
    "    \n",
    "    X_all = np.concatenate((np.vstack(all_sampling_locations1),np.vstack(all_sampling_locations2),np.vstack(all_sampling_locations3)))\n",
    "    X_all,X_all_inverse_ind = np.unique(X_all, axis=0, return_inverse='True')\n",
    "    \n",
    "    # initialzie ADMM \n",
    "    Hu = np.zeros(len(X_all))\n",
    "    lam = np.zeros(len(X_all))\n",
    "    \n",
    "    for epoch in range(num_epochs):\n",
    "        ###############################\n",
    "        ########## TRAIN LOOP #########\n",
    "        ###############################\n",
    "        model.train()\n",
    "        for X, y in train_dl:\n",
    "            optim.zero_grad()\n",
    "            X, y = X.to(device), y.to(device)\n",
    "            y_hat = model(X)\n",
    "            loss = F.mse_loss(y, y_hat)\n",
    "\n",
    "            # Apply EN-S regularization\n",
    "            #spectrum = compute_wht_spectrum(model, H, all_inputs, device=device)\n",
    "            #reg_loss = F.l1_loss(spectrum, torch.zeros_like(spectrum))\n",
    "            #loss += reg_lambda * reg_loss\n",
    "            \n",
    "            #spectrum2 = compute_wht_spectrum(model, H, torch.tensor(X_all, dtype=torch.float), device=device)\n",
    "            wht_out = model(torch.tensor(X_all, dtype=torch.float)).reshape(-1)\n",
    "            wht_diff = wht_out - torch.tensor(Hu, dtype=torch.float) + torch.tensor(lam, dtype=torch.float)\n",
    "            reg_loss = rho/2 * l2_loss(wht_diff,torch.zeros_like(wht_diff))\n",
    "            loss += reg_lambda * reg_loss\n",
    "\n",
    "            loss.backward()\n",
    "            optim.step()\n",
    "        \n",
    "        with torch.no_grad():\n",
    "            y_hat_all = model(torch.from_numpy(X_all).float())\n",
    "        y_hat_all = y_hat_all.numpy().flatten()\n",
    "    \n",
    "    \n",
    "        # run spright to do fast sparse WHT\n",
    "        spright = SPRIGHT('frame', [1,2,3],sampling_matrix1,sampling_matrix2,sampling_matrix3,delays_matrix1,delays_matrix2,delays_matrix3,all_sampling_locations1,all_sampling_locations2,all_sampling_locations3)\n",
    "        spright.set_train_data(X_all,  y_hat_all + lam, X_all_inverse_ind)\n",
    "        spright.model_to_remove = model\n",
    "        flag = spright.initial_run()\n",
    "        if not flag:\n",
    "            continue\n",
    "        spright.peel_rest()\n",
    "        \n",
    "        M = make_system_simple(np.vstack(spright.model.support), X_all)\n",
    "        Hu = np.dot(M,spright.model.coef_)\n",
    "        \n",
    "        # update the dual \n",
    "        lam = lam + y_hat_all - Hu\n",
    "\n",
    "        ###############################\n",
    "        ########## VAL LOOP #########\n",
    "        ###############################\n",
    "        model.eval()\n",
    "        with torch.no_grad():\n",
    "            X, y = next(iter(val_dl))\n",
    "            X, y = X.to(device), y.to(device)\n",
    "            y_hat = model(X)\n",
    "            val_loss = F.mse_loss(y, y_hat).item()\n",
    "            print(f\"Validation Loss: {val_loss:.3f}\")\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################\n",
    "######## Get WHT Matrix and All Inputs #########\n",
    "################################################\n",
    "\n",
    "H = walsh_hadamard_matrix(L=L)\n",
    "all_inputs = np.asarray(list((product((0,1), repeat=L))))\n",
    "H = torch.tensor(H, dtype=torch.float)\n",
    "all_inputs = torch.tensor(all_inputs, dtype=torch.float)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################\n",
    "###### Get Dataset and Split Train/Test ########\n",
    "################################################\n",
    "\n",
    "torch.manual_seed(1)\n",
    "random.seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "n_train = 60\n",
    "bs = 16\n",
    "\n",
    "X_all = torch.tensor([[float(x) for x in g] for g in genotypes], dtype=torch.float)\n",
    "y_all = torch.tensor(brightness, dtype=torch.float)\n",
    "n = X_all.shape[0]\n",
    "ds = TensorDataset(X_all,y_all)\n",
    "\n",
    "# Create Datasets / DataLoaders\n",
    "train_ds, val_ds, test_ds = random_split(ds, lengths=[n_train, n_train, n - 2*n_train])\n",
    "train_dl = DataLoader(train_ds, batch_size=bs, shuffle=True)\n",
    "val_dl = DataLoader(val_ds, batch_size=len(val_ds), shuffle=False)\n",
    "test_dl = DataLoader(test_ds, batch_size=len(test_ds), shuffle=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Validation Loss: 0.263\n",
      "Validation Loss: 0.320\n",
      "Validation Loss: 0.486\n",
      "Validation Loss: 0.691\n",
      "Validation Loss: 0.881\n",
      "Validation Loss: 1.018\n",
      "Validation Loss: 1.056\n",
      "Validation Loss: 1.061\n",
      "Validation Loss: 0.910\n",
      "Validation Loss: 0.788\n",
      "Validation Loss: 0.751\n",
      "Test Loss: 0.653\n"
     ]
    }
   ],
   "source": [
    "lr = 1e-3\n",
    "weight_decay = 0.0\n",
    "reg_lambda = 1e-2  ## with EN regularization\n",
    "num_epochs = 20\n",
    "rho = 0.01\n",
    "device = \"cpu\"\n",
    "\n",
    "torch.manual_seed(1)\n",
    "random.seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "model = Net(L, 2)\n",
    "model = train_model(model, train_dl, val_dl, H, all_inputs, \n",
    "                    lr=lr, weight_decay=weight_decay, reg_lambda=reg_lambda, rho=rho, \n",
    "                    num_epochs=num_epochs, device=device)\n",
    "\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    X, y = next(iter(test_dl))\n",
    "    X, y = X.to(device), y.to(device)\n",
    "    y_hat = model(X)\n",
    "    test_loss = F.mse_loss(y_hat, y).item()\n",
    "    print(f\"Test Loss: {test_loss:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Validation Loss: 1.248\n",
      "Validation Loss: 2.106\n",
      "Validation Loss: 2.118\n",
      "Validation Loss: 2.028\n",
      "Validation Loss: 1.830\n",
      "Validation Loss: 1.336\n",
      "Validation Loss: 1.295\n",
      "Validation Loss: 1.293\n",
      "Validation Loss: 1.265\n",
      "Validation Loss: 1.254\n",
      "Test Loss: 1.091\n"
     ]
    }
   ],
   "source": [
    "lr = 1e-3\n",
    "weight_decay = 0.0\n",
    "reg_lambda = 0     ## without EN regularization\n",
    "num_epochs = 20\n",
    "device = \"cpu\"\n",
    "model = Net(L, 2)\n",
    "model = train_model(model, train_dl, val_dl, H, all_inputs, \n",
    "                    lr=lr, weight_decay=weight_decay, reg_lambda=reg_lambda, \n",
    "                    num_epochs=num_epochs, device=device)\n",
    "\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    X, y = next(iter(test_dl))\n",
    "    X, y = X.to(device), y.to(device)\n",
    "    y_hat = model(X)\n",
    "    test_loss = F.mse_loss(y_hat, y).item()\n",
    "    print(f\"Test Loss: {test_loss:.3f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  },
  "name": "EN_V2.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}