{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phase-Only Compressive Sensing\n",
    "\n",
    "This notebook aims at testing the possibility to recovery a signal from phase-only observations, as described in the papers:\n",
    "\n",
    "* [1] L. Jacques, T. Feuillen, \"The importance of phase in complex compressive sensing\", 2020 ([arXiv](https://arxiv.org/abs/2001.02529))\n",
    "* [2] L. Jacques, T. Feuillen, \"Keep the phase! Signal recovery in phase-only compressive sensing\", 2020 ([arXiv](https://arxiv.org/abs/2011.06499))\n",
    "\n",
    "The code relies on the \"basis pursuit denoising\" (BPDN) solver implemented by [\"SPGL1: A solver for sparse least squares\"](https://friedlander.io/spgl1/), and in particular its [port to Python](https://github.com/drrelyea/spgl1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in [1,2], the phase-only compressive sensing is defined from the following model:\n",
    "\n",
    "$$\n",
    "\\boldsymbol z = {\\rm sign}_{\\mathbb C}(\\boldsymbol \\Phi \\boldsymbol x)\n",
    "\\tag{PO-CS}\n",
    "$$\n",
    "\n",
    "with ${\\rm sign}_{\\mathbb C}(\\lambda) = \\lambda/|\\lambda|$ the complex sign operator, a complex Gaussian random matrix $\\boldsymbol \\Phi \\sim \\mathbb C \\mathcal N^{m \\times n}(0,2)$, and $\\boldsymbol x$ an $k$-sparse vector in $\\mathbb R^n$.\n",
    "\n",
    "The papers above develop a method to recover, up to a global amplitude, $\\boldsymbol x$ from $\\boldsymbol z$.\n",
    "\n",
    "**Remark:** All the codes below are based on a complex Gaussian matrix $\\boldsymbol \\Phi$. However, complex Bernoulli and random partial Fourier sensing are also implemented through the function ``generate_sensing_mtx(m,n,model='gaussian')``, with ``model`` in ``{'gaussian','bernouilli','fourier'}``. Feel free to play with them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Requirements**:\n",
    "* Numpy, SciPy, MatplotLib\n",
    "* [SPGL1](https://github.com/drrelyea/spgl1)\n",
    "* [tqdm](https://github.com/tqdm/tqdm) (just for progress bars ;-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Preliminaries**: Loading important modules. If you cannot run or install the progress bar of ```tqdm```, just comment the lines where it's used. See below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "from scipy.sparse.linalg import LinearOperator\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "#from scipy.sparse.linalg import LinearOperator\n",
    "#from scipy.sparse import spdiags\n",
    "import spgl1\n",
    "\n",
    "# for nice progress bars (comment it if it doesn't work)\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "\n",
    "# Initialize random number generators\n",
    "np.random.seed(43273289)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Noiseless PO-CS reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a few important functions, namely the complex sign operator ${\\rm sign}_{\\mathbb C}$, the constant $\\kappa = \\sqrt{\\pi/2}$, the signal-to-noise ratio metric, as well as a shortcut to call the BPDN solver from SPGL1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# numerical error constant: 10*np.finfo(np.float32).eps\n",
    "epsilon = 10*np.finfo(np.float32).eps \n",
    "\n",
    "# complex sign operator\n",
    "signc = lambda x: x / (np.absolute(x) + (x == 0)*epsilon) \n",
    "\n",
    "# the kappa constant (see [1])\n",
    "kappa = np.sqrt(np.pi/2)\n",
    "\n",
    "# the SNR metric between two signals (in dB)\n",
    "snr = lambda a,b: 20*np.log10(np.linalg.norm(a)/np.linalg.norm(a-b)) \n",
    "\n",
    "# simplifying the call to BP from SPGL1 and storing parameters in the same time\n",
    "#bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-7, bp_tol=1e-7)[0]\n",
    "#bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-8, bp_tol=1e-8, iscomplex=False)[0]\n",
    "bpalg = lambda A, b: spgl1.spg_bp(A, b, opt_tol=1e-7, bp_tol=1e-7)[0]\n",
    "\n",
    "# simplifying the call to BPDN from SPGL1 and storing parameters in the same time\n",
    "bpdnalg = lambda A, b, sigma: spgl1.spg_bpdn(A, b, sigma)[0]\n",
    "\n",
    "# a function to generate the j^th canonical vector in d dimensions\n",
    "def generate_canonic_vec(d,j):\n",
    "    e_j = np.zeros(d)\n",
    "    e_j[j] = 1.0\n",
    "    return e_j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function to generate real $k$-sparse signals in $\\mathbb R^n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_sparse_sig(k,n):\n",
    "    idx = np.random.permutation(n)\n",
    "    x = np.zeros(n)\n",
    "    x[idx[0:k]] = np.random.randn(k)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define several functions generating $m \\times n$ random sensing matrix constructions, namely complex Gaussian, complex Bernoulli, and partial Fourier random matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_complex_Gaussian_mtx(m,n):\n",
    "    Phi = np.random.randn(m, n) + np.random.randn(m, n) *1j;\n",
    "    A = Phi/np.sqrt(m)\n",
    "    return A\n",
    "\n",
    "def generate_complex_Bernoulli_mtx(m,n):\n",
    "    Phi = np.sign(np.random.randn(m, n)) + np.sign(np.random.randn(m, n)) *1j;\n",
    "    A = Phi/np.sqrt(m)\n",
    "    return A\n",
    "\n",
    "class generate_partial_Fourier_mtx(LinearOperator):\n",
    "    def __init__(self, m, n):\n",
    "        self.m = m\n",
    "        self.n = n\n",
    "        self.shape = (m, n)\n",
    "        self.dtype = np.complex128\n",
    "        self.idx = np.random.permutation(n)[0:m]\n",
    "        \n",
    "    def _matvec(self, x):\n",
    "        z = np.fft.fft(x) / np.sqrt(self.m) \n",
    "        return z[self.idx]\n",
    "    \n",
    "    def _rmatvec(self, y):\n",
    "        z = np.zeros(self.n, dtype=complex)\n",
    "        z[self.idx] = y\n",
    "        return np.fft.ifft(z) * self.n / np.sqrt(self.m)\n",
    "\n",
    "def generate_sensing_mtx(m,n,model='gaussian'):\n",
    "    if model == 'gaussian':\n",
    "        return generate_complex_Gaussian_mtx(m,n)\n",
    "    elif model == 'bernoulli':\n",
    "        return generate_complex_Bernoulli_mtx(m,n)\n",
    "    elif model == 'fourier':\n",
    "        return generate_partial_Fourier_mtx(m,n)\n",
    "    else:\n",
    "        print(\"unknown model\")\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function to turn complex linear operators (applied to real vectors) into real linear operators with two times as many rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Op2Real(LinearOperator):\n",
    "    def __init__(self, A):\n",
    "        self.A = A\n",
    "        m, n = A.shape\n",
    "        self.shape = [2*m, n]\n",
    "        self.dtype = np.float64\n",
    "        \n",
    "    def _matvec(self, x):\n",
    "        return np.concatenate(((self.A @ x).real,\n",
    "                               (self.A @ x).imag))\n",
    "    \n",
    "    def _rmatvec(self, y):\n",
    "        m = self.A.shape[0]\n",
    "        return (self.A.T @ y[0:m]).real + (self.A.T @ y[m:]).imag\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following [1], this function constructs the equivalent matrix $\\boldsymbol A_{\\boldsymbol z}$ from both the phase measurements $\\boldsymbol z$ and the sensing matrix $\\boldsymbol A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class generate_Az(LinearOperator):\n",
    "    def __init__(self, A, z):\n",
    "        self.A = A\n",
    "        self.ismat = (type(A) == np.ndarray)\n",
    "        self.z = z\n",
    "        self.D = np.diag(z)\n",
    "        self.shape = (A.shape[0]+2, A.shape[1])\n",
    "        self.dtype = np.float64     \n",
    "        \n",
    "    def _matvec(self, x):\n",
    "        \n",
    "        m = self.A.shape[0]\n",
    "        n = self.A.shape[1]\n",
    "        \n",
    "        H_z_x = self.D.conj().T @ (self.A @ x)\n",
    "        \n",
    "        alpha_adj_x = np.array([np.ones(m) @ H_z_x / (kappa * np.sqrt(m))])\n",
    "                \n",
    "        return np.concatenate([alpha_adj_x.real, alpha_adj_x.imag, H_z_x.imag])\n",
    "    \n",
    "    def _rmatvec(self, y):\n",
    "        # passed the adjoint test\n",
    "        m = self.A.shape[0]\n",
    "        n = self.A.shape[1]\n",
    "        \n",
    "        if self.ismat:\n",
    "            alpha = (self.A.conj().T @ (self.D @ np.ones(m))) / (kappa * np.sqrt(m))\n",
    "            return (alpha.real * y[0]) - (alpha.imag * y[1]) - (self.A.conj().T @ (self.D @ y[2:])).imag\n",
    "\n",
    "        else:\n",
    "            alpha = (self.A.H @ (self.D @ np.ones(m))) / (kappa * np.sqrt(m))\n",
    "            return (alpha.real * y[0]) - (alpha.imag * y[1]) - (self.A.H @ (self.D @ y[2:])).imag\n",
    "            \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Single reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define the sensing model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dimensions\n",
    "n = 128 # space dimension\n",
    "m = 50  # number of measurements\n",
    "k = 10  # sparsity level\n",
    "\n",
    "# Generating the complex Gaussian matrix\n",
    "#A = generate_complex_Gaussian_mtx(m,n)\n",
    "sensing_model = 'gaussian'\n",
    "\n",
    "A = generate_sensing_mtx(m, n, model = sensing_model)\n",
    "\n",
    "# Generating a k-sparse signal, random k-length support, and Gaussian entries\n",
    "x = generate_sparse_sig(k,n)\n",
    "\n",
    "# Linear measurements\n",
    "y = A @ x\n",
    "\n",
    "# Phase only measurements\n",
    "z = signc(y)\n",
    "\n",
    "# Definition of the matrix A_z for sparse signal recovery in PO-CS, such that A_z x0 = (1, 0, ..., 0)  \n",
    "A_z = generate_Az(A,z)\n",
    "\n",
    "#A_z = build_Az(A,z)\n",
    "\n",
    "# Let x0 such that ||Ax0||_1 = kappa sqrt(m) (see [1])\n",
    "# This is the signal we hope to estimate, as the amplitude of x is lost (see [1])\n",
    "x0 = x * kappa * np.sqrt(m) / np.linalg.norm(y,1)\n",
    "\n",
    "# The measurement vector e_1 such that A_z x0 = e_1\n",
    "e_1 = generate_canonic_vec(m+2,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute that $\\|\\boldsymbol A \\boldsymbol x_0 - \\boldsymbol e_1\\|$ is 0, up to numerical errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "||Ax0||_1 = 1.0000000000000002\n",
      "||Ax0 - e_1|| = 2.056073113926681e-16\n"
     ]
    }
   ],
   "source": [
    "print(f\"||Ax0||_1 = {np.linalg.norm(A @ x0,1)/(kappa*np.sqrt(m))}\")\n",
    "print(f\"||Ax0 - e_1|| = {np.linalg.norm((A_z @ x0) - e_1)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve the linear and phase-only basis pursuit (BP) problems, defined by, respectively\n",
    "\n",
    "$$\n",
    "\\boldsymbol x_{\\rm lin} \\in \\arg\\min_{\\boldsymbol u} \\|\\boldsymbol u\\|_1\\ {\\rm s.t.}\\ \\boldsymbol A \\boldsymbol u = \\boldsymbol y,\n",
    "\\tag{BP/Lin}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\boldsymbol x_{\\rm PO-CS} \\in \\arg\\min_{\\boldsymbol u} \\|\\boldsymbol u\\|_1\\ {\\rm s.t.}\\ \\boldsymbol A_z \\boldsymbol u = \\boldsymbol e_1,\n",
    "\\tag{BP/PO-CS}\n",
    "$$\n",
    "\n",
    "rembering that $\\boldsymbol x_{\\rm lin}$ should estimate $\\boldsymbol x$, and $\\boldsymbol x_{\\rm PO-CS}$ should estimate $\\boldsymbol x_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The SNR between x_lin and x is 138.6 dB\n",
      "The SNR between x_pocs and x0 is 110.8458148 dB\n"
     ]
    }
   ],
   "source": [
    "# Linear model (note: running SPGL1 on complex mtx seems wrong). \n",
    "A_equiv = Op2Real(A)\n",
    "y_equiv = A_equiv @ x\n",
    "x_lin = bpalg(A_equiv, y_equiv)\n",
    "\n",
    "# Phase-only model\n",
    "x_pocs = bpalg(A_z, e_1)\n",
    "\n",
    "# SNRs\n",
    "print(f\"The SNR between x_lin and x is {snr(x_lin, x):.4} dB\")\n",
    "print(f\"The SNR between x_pocs and x0 is {snr(x_pocs,x0):.10} dB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x12d4540d0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApw0lEQVR4nO3dfXQV9Z348fcnyQ0xgKiASA0k0AKK8mhwfS4VqqAVZNEfcLIVjnVTtXZxd6tLT7aurebUVVdba2ubWp/TqsUWo9VVGq26PhIUUXkQRB6CKCEoxQRIQj6/P2YmuQk3T8xNZm7m8zrnnjsz95uZb+4lHz738535jqgqxhhjer+0oDtgjDGmZ1jAN8aYiLCAb4wxEWEB3xhjIsICvjHGRERG0B1oy6BBgzQvLy/obhhjTEpZuXLlLlUdnOi10Ab8vLw8Kioqgu6GMcakFBHZ0tZrVtIxxpiIsIBvjDERYQHfGGMiIrQ1/ETq6+uprKxk//79QXcl8rKyssjJySEWiwXdFWNMJ6VUwK+srKR///7k5eUhIkF3J7JUlerqaiorKxkxYkTQ3THGdFJKlXT279/PwIEDLdgHTEQYOHCgfdMKWmkp5OVBWprzXFoadI9MyKVUhg9YsA8J+xwCVloKhYVQW+usb9nirAMUFATXLxNqKZXhG2NcRUVobS0P8W1qOcLZVlsLRUXB9suEmgX8FHPjjTdy++23H7J92bJlrFmzpsv727x5M7///e+b1h944AGuueYaX300PWDrVjYwioU8xJPMbrHdmLb07oAfUI2zoaGhR44Tr72A315/Wgd8kyKGD+fvHAnAXvq32G5MW3pvwPdqnFu2gGpzjdNn0L/pppsYM2YMZ511FgsWLGjKtqdOncq1115Lfn4+P//5zykvL2fSpEmMGzeOyy+/nAMHDgDOlBG7du0CoKKigqlTpwJO5n755ZczdepURo4cyV133dV0zOLiYkaPHs1ZZ53F+vXrD+nTa6+9RllZGddddx0TJ07ko48+OqQ/ixYtYunSpU0/069fPwCWLFnCK6+8wsSJE7nzzjsB+OSTT5gxYwajRo3i+uuv9/V+mW5SXExNn4EA1NDX2ZadDcXFAXbKhF3KDdp2WlFR84CWx6txHuag1ooVK3jiiSd49913qa+vZ/LkyZxyyilNr9fV1VFRUcH+/fsZNWoU5eXljB49mssuu4x77rmHa6+9tt39r1u3jhdffJG9e/cyZswYrrrqKlavXs2jjz7KqlWraGhoOOSYAGeccQazZs3iW9/6Fpdccskh/QFYtGhRwmPecsst3H777Tz99NOAU9JZtWoV77zzDn369GHMmDF8//vfZ9iwYYfxjpluU1BA7btfgduglr6Qm+sEexuwNe3ovRl+W7VMHzXOV199ldmzZ5OVlUX//v256KKLWrw+b948ANavX8+IESMYPXo0AAsXLuTll1/ucP8XXnghffr0YdCgQRx77LF89tlnvPLKK8yZM4fs7GyOPPJIZs2a1en+ev3pqmnTpjFgwACysrIYO3YsW7a0OReTCVDNqd9wnn94M2zebMHedCgpAV9EZojIehHZKCJL2mk3V0RURPKTcdx2tVXL7MYaZ9++fTtsk5GRQWNjI8Ah57H36dOnaTk9Pd33WEB8f+KP29jYSF1dXZs/l+x+mO7hfYFt/UXWmLb4Dvgikg78EpgJjAUWiMjYBO36A4uBN/0es1OKi52aZjyfNc4zzzyTp556iv379/Pll182lUFaGzNmDJs3b2bjxo0APPzww3z9618HnBr+ypUrAXjiiSc6POY555zDsmXL2LdvH3v37uWpp55K2K5///7s3bu3zf3EH7esrIz6+vpO/ZwJr5qals+hZheJhUIyMvxTgY2quklV64BHIf48sSY3Af8N9MzlmQUFUFLi1DZFnOeSEl9fe6dMmcKsWbMYP348M2fOZNy4cQwYMOCQdllZWdx///1ceumljBs3jrS0NK688koA/uu//ovFixeTn59Penp6h8ecPHky8+bNY8KECcycOZMpU6YkbDd//nxuu+02Jk2axEcffXTI6//8z//MSy+9xIQJE3j99debsv/x48eTnp7OhAkTmgZtTWpImQy/m06gMIdBVX09gEuAe+PWvw3c3arNZOAJd/lvQH4b+yoEKoCK4cOHa2tr1qw5ZFtP27t3r6qq1tTU6CmnnKIrV64MuEfBCcPnEWU//rEqqM6eHXRPOpCbqwp6Kz/QNZzgdBqc7SbpgAptI153+6CtiKQBdwD/3lFbVS1R1XxVzR88OOEdugJXWFjIxIkTmTx5MnPnzmXy5MlBd8lEVMpk+Fu3UkeM67mNP7CgxXbTs5JxWuZ2IP6cvRx3m6c/cDLwN3f+leOAMhGZpaopdw9Du0jJhEXK1PCHD6dui3PtiV0kFqxkZPgrgFEiMkJEMoH5QJn3oqruUdVBqpqnqnnAG0BKBntjwiRlMvziYuqOOAqAL3Eu+LOLxILhO8NX1QYRuQZ4DkgH7lPVD0TkJzi1pLL292CMORwpE/ALCqjfkwXfcwO+XSQWmKRcaauqzwDPtNp2QxttpybjmMZEXcqUdIC6i+bC92DvhQvg6QUd/4DpFr33SltjermUyfAB7zq/L78Mth9RZwG/Gz3wwAN88sknTetXXHHFYU1h3JrNcGkgxTJ8C/ihYAG/G7UO+Pfeey9jxx5yEXKXWcA30JzZ19VB2Ge/8AK+XdQdLAv4h+GRRx7h1FNPZeLEiXz3u9/l4MGDLFq0iJNPPplx48Zx5513snTpUioqKigoKGDixIns27ePqVOnNs1e2a9fP6677jpOOukkpk+fzltvvdU0NXJZmTPOvXnzZs4++2wmT57M5MmTee2114BDpzQ+ePAg1113HVOmTGH8+PH85je/Cey9MT0nPrMPe1nHncnDMvyApez0yNdeC6tWJXefEyfCz37Wfpu1a9fy2GOP8eqrrxKLxbj66qu5+eab2b59O++//z4AX3zxBUcddRR33303t99+O/n5h84VV1NTw7nnnsttt93GnDlz+M///E+WL1/OmjVrWLhwIbNmzeLYY49l+fLlZGVlsWHDBhYsWEBFRcUhUxqXlJQwYMAAVqxYwYEDBzjzzDM577zzGDFiRHLfIBMq8UG+thaOPDK4vnTESjrhkLIBPyjl5eWsXLmyaU6bffv2MWPGDDZt2sT3v/99LrzwQs4777wO95OZmcmMGTMAGDduHH369CEWizFu3Dg2b94MQH19Pddccw2rVq0iPT2dDz/8MOG+nn/+eVavXt10g5M9e/awYcMGC/i9XE0NHHMM7N4d/jp+fMBXdaa3Mj0vZQN+R5l4d1FVFi5cyE9/+tMW24uLi3nuuef49a9/zeOPP859993X7n5isRjulcekpaU1TUmclpbWNB3xnXfeyZAhQ3j33XdpbGwkKyurzT794he/4Pzzz/f765kUUlsLX/2qE/DDXtLxAn5jI+zbd+hEtqZnWA2/i6ZNm8bSpUvZuXMnALt372bLli00NjYyd+5cbr75Zt5++23A/9TDe/bsYejQoaSlpfHwww9z8ODBhPs9//zzueeee5qmPP7www+pCXvKZ3zxBmq9KafC/nHH337ByjrBSdkMPyhjx47l5ptv5rzzzqOxsZFYLMYdd9zBnDlzmm4w4mX/ixYt4sorr+SII47g9ddf7/Kxrr76aubOnctDDz3EjBkzEk5pvGjRIhYvXszmzZuZPHkyqsrgwYNZtmxZ0n5nEz5eRu8F/LBn+N6gLTgB/9hjg+tLlIkzm2b45Ofnq3dGi2ft2rWceOKJAfXItGafR3C2b4ecHGda+ZISKCuDVnfcDJXHHoP5853lVatgwoRAu9OrichKVU14V0Er6RiTglItw7eSTjhYwDcmBXk1e6vhm65IuRq+qjad3RKI6mrn+3RdHWRmwvHHw8CBwfUnIGEtBUZFqmX48TV8u9o2OCmV4WdlZVFdXR1csKmudu7H6aUrdXXOenV1MP0JiKpSXV3d5mmipvtZhm8OR0pl+Dk5OVRWVlJVVRVMByor0YMH2cMA+vMl6TinSVJd7YygRUhWVhY5Efudw8TL6I85puV6WFnAD4eUCvixWCzYq0dPOon1OoozWM/dfI/v8Stnu4hzRYkxPcTL6Pv1cy5iSqUM30o6wUlKSUdEZojIehHZKCJLErz+byKyRkRWi0i5iOQm47g9bvhwqnC+Q3vP3nZjepKX0WdnQ9++qZPhp6dbhh8k3wFfRNKBXwIzgbHAAhFpPQfwO0C+qo4HlgK3+j1uIIqLqcp0yhhNAd/uzWkC4GX0ffumRoZfXw8ZGc43Egv4wUlGhn8qsFFVN6lqHfAoMDu+gaq+qKpeDvIGkJrF34ICdn37WgB2Mci5N2dJid2b0/S4VMzwMzOhf38r6QQpGTX844FtceuVwD+00/47wLNJOG4gqkae5jx/Yx68MC/g3pioqqmBtDTo0yc1Mnwv4FuGH6weHbQVkX8C8oGvt/F6IVAIMDykdXHvBKGgThQyBpyMPjvbOV8glTJ8C/jBSkZJZzswLG49x93WgohMB4qAWap6INGOVLVEVfNVNX/w4MGJmgRu166Wz8YEoabGCfSQGhl+fb2VdMIgGQF/BTBKREaISCYwHyiLbyAik4Df4AT7nUk4ZmC8zH7XLudGDsYEwcvwIXUy/FjMMvyg+Q74qtoAXAM8B6wFHlfVD0TkJyIyy212G9AP+KOIrBKRsjZ2F3pewG9ogC++CLQrJsJqa1tm+KkQ8L0M3wJ+cJJSw1fVZ4BnWm27IW55ejKOEwa7djmZSn29s3z00UH3yERRTU3LDD/sJZ34Gr6VdIKTUnPphEFVFYwe3bxsTBBSNcO3kk6wLOB3QU2Ncz/Ose5lZRbwTVBaZ/i1teEeU4oftK2psZlIgmIBvwu8M3O8mzzZmTomKK0zfFXYvz/YPrUnftAWwl+C6q0s4HeBl9F7Ad8yfBOU1hm+ty2s4ks6YGWdoFjA7wIvo8/NhSOOsAzfBKd1hu9tC6v4s3TAAn5QLOB3gZfRDx7sPCzDTyGlpZCX58xHkJfnrKew+Azfew5zhu/V8L0M387UCUZKzYcfNC/ADxrkPCzgp4jSUigsbE6Bt2xx1iElJ75rbHTq9V6G7z2HPcOPxSzDD5pl+F3gnYM/YICT4VtJJ0UUFUFtLT9lCWVc5GyrrXW2p6D4mTLjn8Oc4VsNPxws4HdBVZWT2YtYSSelbN1KHTFu5EZ+zuIW21ORF/BTLcO3kk7wLOB3gRfwwUo6KWX4cN7nZOrow0pOoRFp2p6KvEw+FTN8K+kEywJ+F+za5WT24Dx7F2KZkCsupiLzDAD2cBQf8dWUvlNZKmb4rQdtLeAHwwJ+F1RVtQz4YHX8lFBQQMWZixGcyztXDJqZ0ncqS9UMP/7CKyvpBMMCfhe0Lul420z4VXwxim+cm0ZWFlR8+66UDfZw6KBtKmT4XkknFnPu0mUZfjAs4HdSQwN8/rll+Klo/3547z047TSYNAkqKoLukT/xNzCH8Gf4qs0BH2wCtSBZwO+k6mrnuXXAT7kMv5ddgNQZq1c7/2Hn5zuPt9+GgweD7tXha53hZ2Q4wTSsGX5Dg/McH/CtpBMMC/idFH/RVfxzSgV87wKkLVuctMu7AKmXB30vo/cCfk0NrF8fbJ/8aJ3hQ7hvc1hf7zx7Ad9ughIcC/id5JVuvMz+6KMhPT3FSjpFRWytHcg8HmUVE5xtKXwBUmdVVMCxx0JOjhPwvW2pqnWGD+G+zWFdnfMciznPVtIJTlICvojMEJH1IrJRRJYkeL2PiDzmvv6miOQl47gJeSULEee7bvyzNydCR9sSvFY1bR4AgwvOg9JS0v5QykCtoqr4N0nZf0+0//uW3XyLp3mceVzIX9jOV5z3bMuWlOj/4e6r4v7V5O99Afl9KWNW/p6+UkPFwl+E73frZPuaa/4DgL5Txjr/3ktLyf5sEzW/+0PX9++V9RL93SSj/3l51P1+KeBm+KWl9F/1Mnuffy05/WmvfRg+S690evXVzaXUrn42yaSqvh5AOvARMBLIBN4FxrZqczXwa3d5PvBYR/s95ZRTtMseeUQ1O1vVKVgk9fFLrlJQ3cEQ1VhMNTNTx/K+zuGJbjlesh/1pOtM/qLp1OvP+Bftx991Eit1L30D71t3Pmo4QtNo0Bu4selzO5uX9HReDbxvh/u4kRsUVBtIa/qdJvCOzmLZ4e3T3Ud39Xdb1tcUVH/7nddVs7N1Dk/oON5NXn+6uf+BPrKznbjWBUBFW3E1Iwn/Z5wKbFTVTQAi8igwG1gT12Y2cKO7vBS4W0TE7VzyFBVRVZtNLskvrNcTQ2hkINVQ74xCHctOlnEx2YS0eBqnkTQOkMVvKKSQ3zKKDVzEUwykmnRSeASzA42k0Ug6p7CyqZicTwV38m8p8bklUkcmR1BLOo1Q71xb0J+9PMVFXfqdruFubuU/mt6X03mNd71SXxL8O//DTdxA/X7n7yVz2eNQW8uR/J33GH9IX6/nVm7kx1BfTwPpTGK1c5FcG37IT/kRN0N9PXXEmMgHbCYvaf33K0Y9TzCX6ZQDsIoJTKOcfRzR4c9OYQUvMbW55Jqk04iTEfCPB7bFrVcC/9BWG1VtEJE9wECgRQVcRAqBQoDhh3PZ+9atHEE213B313+2E05gHTEamtZ/wg085U3GlQLGs5p/wvmKeAHP8hQX8TemBtupHtCPLzmf55rW/4W76MMBDpIeYK/8Gcd7LdZ/wg08y8xO//wfWEAF+S22vck/cBpvcBb/57t/pRTwphsG6nBGa2PVnwLwr9zJsexs0f6PXMoLnOsEfGAHQ3mfcVzAXziJDw7Z/x+5lOc5zwn4wGbyWMtYLqKME1jnu//JcCf/yguc2xTwX+YcdjOQxfyMTOra/dnhxM3zlMw5n/yUc9wE/RLg3rj1bwN3t2rzPpATt/4RMKi9/R5WSSc3N/ivYPawRwo8vkG5ns1LTesHEQXVH/OjpOx/Gsv1DP5PFXQ1JyuoLh303Tbb/xMPaR6bmtZf5XQF1Wc5P2H7hdyvx7Otaf1/OU9B9WXOCvy99R4nsEYv5k9N61fyKz2K3drY1X3l5nYpDNJOSScZg7bbgWFx6znutoRtRCQDGABUJ+HYLRUXtzx1obvEYs3nmKWiVO//4eqNv/dh/k4ZNFBPrGkf9bG+TduTIZtaanD2WZd1JACZl81v8+8zh0q2c7wzsV0sxrb0EU3bExnJJraTw376QCzGpvTRTdvD4kTWso4TmtbXcQInstabuq9zkjznUzIC/gpglIiMEJFMnEHZslZtyoCF7vIlwAvu/0TJVVDgzJGSm+usp6e3fB440Hl0tK2913Jz4f774b77Dj1OMvbf3e1Tvf+Hu6/Wv7dI+H63HvwsMzhIAxlN+2i457fu9oak9L+vuAE/N5e6/7gBgMzzpjb/fbZ6/4fJdurJZGfOKXD//VTO+zdnO9sS7n+EbAFgy1fOgPvv5+MZV9GHAwxlRzg+m4EDOeGILWzka9QP/ypcdRXr0k/iBNZ3fv+5ucmf8+lwyjitH8AFwIc4pZoid9tPgFnuchbwR2Aj8BYwsqN9HlZJxxjTKbNmqU6Y0Lz++edO9eCOO5Kz/+98R3XoUGf5Jbdy9MILbbd/8kmnzVtvOeuLF6v266fa2Ji4/f851SJ95hlnfe5c1TFjktP3ZHnoIaePa9c2v7+33tr9x6Wbz9JBVZ8Bnmm17Ya45f3Apck4ljHGv4yM5ikPoHk5IykRwbkQzLvyt/WFV4kMc4vClZUwZYrzPGyY80UgkRFOxYePP25+9raFxQluNWftWvjii5bbgpKkj9cYk0piscQBv72g3BXxV/56Ab+9oYacHOd527bmZ29bIscd58y6GR/wTzvNX5+TzQvu69bBnj3O8oknBtcfsIBvTCRlZDTPcQPNy8nM8BsanGDfmYA/aJATwOMD/sknt90+Lc3J6DdtcrLnzz8PX4bfvz8cf3xzhp+Z6Vw8GyQL+MZEUE+UdMAp67SePC0RESejr6x02n/6afsZPjgB/uOPm7P8kSP99zvZTjyxOcMfPTp57+/hssnTjImgWCxxhp/Mkg44Ab8zGT44Nftt2+CTT5wT0IcNa7/9yJFOhr/JPRMzbBk+OGWddetgzZrg6/dgAd+YSOrJDL8zg7bQnOFXVjavt2fECCdzfvttZz2sGf7evbBxY/D1e7CAb0wk9cSgLXQ9w9++3Zm81Vtvj5fRl5c705UPGHD4/e0u8Vm9ZfjGmEB096Bt/G0XuxLwGxpg5crm9fZ4Gf2KFeEs5wCcuOaJ5uUfXBj4zYYs4BsTQW2VdLojw+/MoC00l3Bef905w+XII9tv7wX5xsZwlnMoLeW46y/jSJxzMkfv+Fvgd5izgG9MBLU1aBtkDd/L6Feu7Di7B6eEc8wxznIoM/yiImRfLSewjuFsoS+1gd9hzk7LNCaCMjKcG7mrOqdEdtegbW1t1wZtwWnf0YAtAKWljPhyHLsZz8jfFcGEscmdd8Yvd1rjJdzCHgYcsj0IluEbE0Fe8D3o3vumu0/LTE9vnhOsLQMHQlaWs9xhhl9aCoWFjKhz7kY/YndF4OWSQ7j39JjDMhbx4CHbg2AB35gI8jJ5L9B392mZnZnB2bv4CjqR4RcVQW1t03TII9kUeLnkEImma0/ydMddZQHfmAjyArsX6Lt70LZTU/aXljJs22sADPvF9e1n625Z5EL+wgyeJY/NLbaHQvx07SLdM91xF1nANyaCvMDuZfjJHrRNT3fmxvEy/A7/I3FLNDkHNgKQs/vd9ks0blnkHF7hWS5ovvVogOWShAoKYPNm51SizZsDH2OwgG9MBLWV4SdzrhdviuROlXTcEs0w9/bYw9jWfokmhOWSVGAB35gI8jJuL9Ane9AWnPjb6YDvlmLO5FW+ysaOSzQhLJekAl8BX0SOEZHlIrLBfT46QZuJIvK6iHwgIqtFZJ6fYxpj/OvuQVtozvA7VcN3SzEX8CwbGUU2+1psTyhk5ZJU4DfDXwKUq+oooNxdb60WuExVTwJmAD8TkaN8HtcY40PrDD/Zg7bQxZKOlWh6hN+APxuaTjB9ELi4dQNV/VBVN7jLnwA7gcE+j2uM8aF1hp/sQVtovutVpwZtrUTTI/x+vENUdYe7/CkwpL3GInIqkIlzs/NErxcChQDDwzbabkwv0lODttXVztk6nTots6DAAnw36/DjFZG/AscleKnF8LmqqohoO/sZCjwMLFTVxkRtVLUEKAHIz89vc1/GGH96YtC2b19nzLVfv04GfNPtOgz4qjq9rddE5DMRGaqqO9yAvrONdkcCfwGKVPWNw+6tMSYpenrQtk+f5O3XHD6/NfwyYKG7vBB4snUDEckE/gw8pKpLfR7PGJMEPZXhd2VqBdP9/Ab8W4BvisgGYLq7jojki8i9bpv/B5wDLBKRVe5jos/jGmN86MkMv1ODtqZH+Pp4VbUamJZgewVwhbv8CPCIn+MYY5Ir0aCtSMczWnZF376wbx/s328ZfljYlbbGRFCikk4ys3toPq1+zx4L+GFhAd+YCEpU0kl2wPdmzPz8cwv4YWEB35gISpThJ7vO7gX8/futhh8WFvCNiaCezPDBMvywsIBvTAQlGrTtrgwfLOCHhQV8YyKoJwZtLeCHjwV8YyLISjrRZAHfmAjqyUHb+OOZYFnANyaCLMOPJgv4xkSQDdpGkwV8YyKoJ6+0BQv4YWEB35gISlTSSXaGn5XlzM8DFvDDwgK+MRHUExm+SHNZxwZtw8ECvjER5M2K2Z2DttAc8C3DDwcL+MZEUFqa8+jO0zLBAn7YWMA3JqJisZZn6ViG3/v5CvgicoyILBeRDe7z0e20PVJEKkXkbj/HNMYkR0ZG9w7agtXww8Zvhr8EKFfVUUC5u96Wm4CXfR7PGJMk8Rl+dwzagmX4YeM34M8GHnSXHwQuTtRIRE4BhgDP+zyeMSZJWmf4FvB7P78Bf4iq7nCXP8UJ6i2ISBrwP8APOtqZiBSKSIWIVFRVVfnsmjGmPa0zfBu07f06/D9dRP4KHJfgpaL4FVVVEdEE7a4GnlHVSvGuwmiDqpYAJQD5+fmJ9mWMSZKMjO4ftPWutrWAHw4dfsSqOr2t10TkMxEZqqo7RGQosDNBs9OBs0XkaqAfkCkiX6pqe/V+Y0w3s0Hb6PFb0ikDFrrLC4EnWzdQ1QJVHa6qeThlnYcs2BsTvG4ftC0tpe8DvwQg8/xvQGlpkg9guspvwL8F+KaIbACmu+uISL6I3Ou3c8aY7tOtGX5pKRQW0vfvnwCQ+cnHUFhoQT9gvgK+qlar6jRVHaWq01V1t7u9QlWvSND+AVW9xs8xjTHJ0a0ZflER1NbSlxoAMqmD2lpnuwmMXWlrTER162mZW7cCcA4vcxFlDKaqxXYTDAv4xkRU/Fk6ST8tc/hwACbzDmXMJkZDi+0mGBbwjYmobp1Lp7i45R1QwFkvLk7iQUxXWcA3JqK8kk5jI6gmOcMvKICSEsjNdSbGz8111gsKkngQ01XdcKmFMSYVxGJw4EBzHT/pp2UWFFiADxnL8I2JKC/D98o63XGlrQkXC/jGRJQ3aOtl+HY1bO9nAd+YiPIGbS3Djw4L+MZEVOuSjmX4vZ8FfGMiysvwu23Q1oSOBXxjIsoy/OixgG9MRLUetLUMv/ezgG9MRNmgbfRYwDcmoqykEz0W8I2JKBu0jR4L+MZElGX40eMr4IvIMSKyXEQ2uM9Ht9FuuIg8LyJrRWSNiOT5Oa4xxj8btI0evxn+EqBcVUcB5e56Ig8Bt6nqicCpJL7ZuTGmB1lJJ3r8BvzZwIPu8oPAxa0biMhYIENVlwOo6peqWuvzuMYYn7wAf+CA82wlnd7Pb8Afoqo73OVPgSEJ2owGvhCRP4nIOyJym4ikJ9qZiBSKSIWIVFRVVfnsmjGmPV6A37fPebYMv/fr8CMWkb8CxyV4qcXdiFVVRUTbOMbZwCRgK/AYsAj4XeuGqloClADk5+cn2pcxJkm8AO8FfMvwe78OA76qTm/rNRH5TESGquoOERlK4tp8JbBKVTe5P7MMOI0EAd8Y03Msw48evyWdMmChu7wQeDJBmxXAUSIy2F0/F1jj87jGGJ+8AL9/f8t103v5Dfi3AN8UkQ3AdHcdEckXkXsBVPUg8AOgXETeAwT4rc/jGmN8spJO9Pj6P11Vq4FpCbZXAFfErS8Hxvs5ljEmuaykEz12pa0xEWUZfvRYwDcmoizDjx4L+MZEVOtBW8vwez8L+MZEVOuSjmX4vZ8FfGMiyko60WMB35iIskHb6LGAb0xEWYYfPRbwjYmo+EFbEUhPOKWh6U0s4BsTUfElHcvuo8ECvjERFV/SsYAfDRbwjYmo+AzfBmyjwQK+MRFlGX70WMA3JqLiB20tw48GC/jGRJQN2kaPBXxjIsrL6g8csAw/KizgGxNR8Vm9ZfjR4Cvgi8gxIrJcRDa4z0e30e5WEflARNaKyF0iIn6Oa4zxLz6rt4AfDX4z/CVAuaqOAsrd9RZE5AzgTJw7Xp0MTAG+7vO4xhif4oO8lXSiwW/Anw086C4/CFycoI0CWUAm0AeIAZ/5PK4xxicr6USP34A/RFV3uMufAkNaN1DV14EXgR3u4zlVXZtoZyJSKCIVIlJRVVXls2vGmPbEZ/WW4UdDh/+vi8hfgeMSvFQUv6KqKiKa4Oe/BpwI5LiblovI2ar6Suu2qloClADk5+cfsi9jTPJYhh89HX7Mqjq9rddE5DMRGaqqO0RkKLAzQbM5wBuq+qX7M88CpwOHBHxjTM+xgB89fks6ZcBCd3kh8GSCNluBr4tIhojEcAZsE5Z0jDE9Jy3NeYCVdKLCb8C/BfimiGwAprvriEi+iNzrtlkKfAS8B7wLvKuqT/k8rjEmCbxAbxl+NPj6mFW1GpiWYHsFcIW7fBD4rp/jGGO6R0aGXWkbJXalrTER5mX2luFHgwV8YyLMy+wtw48GC/jGRJhl+NFiAd+YCLNB22ixgG9MhHmB3ko60WAB35gIs5JOtFjANybCbNA2WizgGxNhluFHiwV8YyLMBm2jxQK+MRFmg7bRYgHfmAizkk60WMA3JsJs0DZaLOAbE2GW4UeLBXxjIswy/GixgG9MhFmGHy0W8I2JMAv40eIr4IvIpSLygYg0ikh+O+1miMh6EdkoIkv8HNMYkzxW0okWvxn++8A/Ai+31UBE0oFfAjOBscACERnr87jGmCSwDD9a/N7icC2AiLTX7FRgo6pucts+CswG1vg5tjHGP8vwo6UnavjHA9vi1ivdbcaYgFmGHy0dfswi8lfguAQvFanqk8nsjIgUAoUAw4cPT+aujTEJ2Fw60dLhx6yq030eYzswLG49x92W6FglQAlAfn6++jyuMaYDNpdOtPRESWcFMEpERohIJjAfKOuB4xpjOmAlnWjxe1rmHBGpBE4H/iIiz7nbvyIizwCoagNwDfAcsBZ4XFU/8NdtY0wy2KBttPg9S+fPwJ8TbP8EuCBu/RngGT/HMsYkn2X40WJX2hoTYTZoGy0W8I2JMBu0jRYL+MZEmJV0osUCvjERZoO20WIB35gIsww/WizgGxNVpaXEbvoRALGZ06G0NOAOme5mAd+YKCothcJCMnbvBCBjx1YoLLSg38tZwDcmioqKoLaWGPUAZNAAtbXOdtNrWcA3Joq2bgVgJs9SxM3ksbnFdtM7WcA3Jorc2WiH8ik38yPS0BbbTe9kAd+YKCouhuzsltuys53tpteygG9MFBUUQEkJ5OaCiPNcUuJsN72WnX1rTFQVFFiAjxjL8I0xJiIs4BtjTERYwDfGmIiwgG+MMRFhAd8YYyJCVDXoPiQkIlXAFh+7GATsSlJ3gmD9D04q9x2s/0ELuv+5qjo40QuhDfh+iUiFquYH3Y/DZf0PTir3Haz/QQtz/62kY4wxEWEB3xhjIqI3B/ySoDvgk/U/OKncd7D+By20/e+1NXxjjDEt9eYM3xhjTBwL+MYYExG9LuCLyAwRWS8iG0VkSdD96YiIDBORF0VkjYh8ICKL3e3HiMhyEdngPh8ddF/bIyLpIvKOiDztro8QkTfdz+ExEckMuo9tEZGjRGSpiKwTkbUicnoqvf8i8q/uv533ReQPIpIV5vdfRO4TkZ0i8n7ctoTvtzjucn+P1SIyObieN/U1Uf9vc//9rBaRP4vIUXGv/dDt/3oROT+QTrt6VcAXkXTgl8BMYCywQETGBturDjUA/66qY4HTgO+5fV4ClKvqKKDcXQ+zxcDauPX/Bu5U1a8BnwPfCaRXnfNz4H9V9QRgAs7vkRLvv4gcD/wLkK+qJwPpwHzC/f4/AMxota2t93smMMp9FAL39FAf2/MAh/Z/OXCyqo4HPgR+COD+Lc8HTnJ/5ldunApErwr4wKnARlXdpKp1wKPA7ID71C5V3aGqb7vLe3GCzfE4/X7QbfYgcHEgHewEEckBLgTuddcFOBdY6jYJbf9FZABwDvA7AFWtU9UvSKH3H+e+FkeISAaQDewgxO+/qr4M7G61ua33ezbwkDreAI4SkaE90tE2JOq/qj6vqg3u6htAjrs8G3hUVQ+o6sfARpw4FYjeFvCPB7bFrVe621KCiOQBk4A3gSGqusN96VNgSFD96oSfAdcDje76QOCLuD+AMH8OI4Aq4H63JHWviPQlRd5/Vd0O3A5sxQn0e4CVpM7772nr/U7Fv+nLgWfd5VD1v7cF/JQlIv2AJ4BrVfXv8a+pc+5sKM+fFZFvATtVdWXQfTlMGcBk4B5VnQTU0Kp8E/L3/2icLHIE8BWgL4eWG1JKmN/vjohIEU6ZtjToviTS2wL+dmBY3HqOuy3URCSGE+xLVfVP7ubPvK+u7vPOoPrXgTOBWSKyGaeEdi5OTfwot8QA4f4cKoFKVX3TXV+K8x9Aqrz/04GPVbVKVeuBP+F8Jqny/nvaer9T5m9aRBYB3wIKtPkCp1D1v7cF/BXAKPcMhUycwZKygPvULrfe/TtgrareEfdSGbDQXV4IPNnTfesMVf2hquaoah7O+/2CqhYALwKXuM3C3P9PgW0iMsbdNA1YQ4q8/zilnNNEJNv9t+T1PyXe/zhtvd9lwGXu2TqnAXviSj+hISIzcMqas1S1Nu6lMmC+iPQRkRE4g89vBdFHAFS1Vz2AC3BGyT8CioLuTyf6exbO19fVwCr3cQFOHbwc2AD8FTgm6L524neZCjztLo/E+Ye9Efgj0Cfo/rXT74lAhfsZLAOOTqX3H/gxsA54H3gY6BPm9x/4A854Qz3ON6zvtPV+A4Jz5t1HwHs4ZyOFsf8bcWr13t/wr+PaF7n9Xw/MDLLvNrWCMcZERG8r6RhjjGmDBXxjjIkIC/jGGBMRFvCNMSYiLOAbY0xEWMA3xpiIsIBvjDER8f8BCRccTauvimMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x0, 'ro', label=\"ground truth\")\n",
    "plt.plot(np.real(x_pocs), 'b', label=\"estimate\")\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Systematic analysis over a range of measurements\n",
    "Let's now try a range of different measurements over many trials to see when the linear and the phase-only methods succeed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "77efb594cd64493da5c7400fc9fc760e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/26 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sensing model\n",
    "sensing_model = 'gaussian'\n",
    "\n",
    "# Dimensions\n",
    "n = 256\n",
    "k = 20\n",
    "m_vec = np.arange(1,256,10)\n",
    "nb_trials = 100\n",
    "\n",
    "# Inits\n",
    "scores_lin = np.zeros((m_vec.size,nb_trials))\n",
    "scores_po = np.zeros((m_vec.size,nb_trials))\n",
    "\n",
    "# Main Monte Carlo Loop\n",
    "pbar = tqdm(total=m_vec.size) #comment this if it doesn't work for you\n",
    "for m_ind, m in enumerate(m_vec):\n",
    "    \n",
    "    # The measurement vector e_1 such that A_z x0 = e_1\n",
    "    e_1 = generate_canonic_vec(m+2,0)\n",
    "    \n",
    "    for s in np.arange(nb_trials):\n",
    "    \n",
    "        # Generating the sensing matrix according to the selected model\n",
    "        A = generate_sensing_mtx(m, n, model = sensing_model)\n",
    "\n",
    "        # Generating a k-sparse signal, random k-length support, and Gaussian entries\n",
    "        x = generate_sparse_sig(k,n)\n",
    "\n",
    "        # Linear measurements\n",
    "        y = A @ x\n",
    "\n",
    "        # Phase only measurements\n",
    "        z = signc(y)\n",
    "\n",
    "        # Definition of the matrix A_z for sparse signal recovery in PO-CS, such that A_z x0 = (1, 0, ..., 0)  \n",
    "        A_z = generate_Az(A,z)\n",
    "\n",
    "        # Let x0 such that ||Ax0||_1 = kappa sqrt(m) (see [1])\n",
    "        # This is the signal we hope to estimate, as the amplitude of x is lost (see [1])\n",
    "        x0 = x * kappa * np.sqrt(m) / np.linalg.norm(y,1)\n",
    "        \n",
    "        # Linear model (note: running SPGL1 on complex mtx seems wrong). \n",
    "        A_equiv = Op2Real(A)\n",
    "        y_equiv = A_equiv @ x\n",
    "        x_lin = bpalg(A_equiv, y_equiv)\n",
    "\n",
    "        # Phase-only model\n",
    "        x_pocs = bpalg(A_z, e_1)\n",
    "        \n",
    "        # Recording SNRs\n",
    "        scores_lin[m_ind,s] = snr(x_lin,x)\n",
    "        scores_po[m_ind,s] = snr(x_pocs,x0)\n",
    "        \n",
    "    pbar.update(1) #comment this if it doesn't work for you\n",
    "            \n",
    "pbar.close() #comment this if it doesn't work for you\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Displaying now the success rates, by counting how many times over the number of trials, the reconstruction SNR was above a certain threshold ```success_score``` (e.g., 60 dB) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m/k: CS 2.05; PO-CS 3.55\n",
      "PO-CS/CS: 1.7317073170731707\n",
      "mCS = 41, mPO = 71, 2*mCS - 2 = 80\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Recovery Rate')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEOCAYAAACXX1DeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5wUlEQVR4nO3deXiU5dX48e/JBiQhC/sWZJEdKyqiLBVesRWtir7+QNTWtVLrUrVq61Jtq1ZxqS/21arYitb1tVCXqiiCoqiIQF3YBBEEwk4gQEhCljm/P+5JMkkmk0xI8jyTnM91zTUzz3omy3Pmfu5NVBVjjDGmJnFeB2CMMcbfLFEYY4yJyBKFMcaYiCxRGGOMicgShTHGmIgSvA6goXXo0EF79erldRjGGBNTli1btltVO4Zb1+wSRa9evVi6dKnXYRhjTEwRkY01rbNbT8YYYyKyRGGMMSYiSxTGGGMiskRhjDEmIksUxhhjIvIsUYjI0yKyU0RW1LBeROQvIrJORL4WkWObOkZjjDHeliieASZEWH8a0C/4mAo83gQxmQgWLYL77nPPjbWP37ZfPmMRC069j+Uz6rZDtNs3xTlifXs/xuS37eu7T12Jl8OMi0gv4E1VHRpm3ZPAAlV9Kfh+DTBOVbdFOubw4cO1vv0onnnmmWrLhgwZwvHHH09xcTEvvPBCtfXDhg1j2LBh5Ofn88orr4SLh6FDh7Jv3z5effXVautHjhzJgAED2L17N2+++Wa19SeddBJ9+vRh+/btvPPOO9XWjx8/nqysLDZv3sz8+fOrrZ8wYQJdunRh/fr1fPTRR9XWn3HGGXTo0IE1a9awKMzV85xzziEtLZ277srmrru6EQgIcXHKrbfO4cgjd3HSSZO55JJkWrXK5cILXyvfb/r08ezZk8yWLe0IBAQRpUePPbRpU1zp+J07dwEgIWE3o0Yt5sEHT6W42H1/ycrK5auv2pGRAR9++CH33deeDRs6lO976FArNm/OJBAg7PHj4xP49NMOZGTAO++8wx/+0I/PP+8DKImJAe655zPS08fw/POQk5NDSUnIvlrCpD0v87ucW0mghABxfJM4iPyEVDp2zKN162JatWpFYUEGO3dCSsoBUkr30WnnRuIIECCONUlDKG6dAUBRcTFQ+X8tLi6OlNJ8+hUuL99nZ6cjKExKZndOKmlpMKB/EgFVtm3dzdZt6SSX5DGweHX59t+2PorCpLYUF1f+uQKkpwvdUwtova7i+GuTBpLZqwSAtLQ01q9vQ3Gxln/2cMcvSGxbvj4+PkC3bvsAyMjIYNOKQ/TNrxz/7uJMCgsTSUhIJE6EQCBASWlJjcfPT0ilS+fdxMW7n8+uXakUFibS78h4MjLi2bNxN+mb1pTvU/Z7SExMRBBKS0spDZSWf+6q58jt0Y/dRa3Zf6A16ekFpKUVAtC6dWc2fg+JhbkMKFpZ5fhtSUpMBKCktIRAIBD8mRbQNnCg0u/5+9SBpHRzn6+0JI6t29IRhMTg/uGOX5DYlsQEtz4pqYDMdvvL98/dHF8p/l1derPlQHsCAan425E4EhJcN7ji4mLalBxgYPFqBKWQ1nz35HyOmjqy2t9EJCKyTFWHh1vn5w533YHNIe+zg8uqJQoRmYorddCzZ88mCa65CwSEb7/txG9/25q33oLs7B4h6+Cbb7pw5JG7KCoSFi6ELl2SK+3/3XcdyctrXf5eVdi8uX2186xd6547d25Lu3ZdKC6OQzUOUDZtakdJScW22dmZrF3bpdL+J7KIcSxggY7js83V/zFC99+woQOqAgglJbB8eXt69ICFCwHak0gRJ/M+5zKbs3mNjuwu3zeeAEOKV0IxsKnyOboC7K+8LJ4Ag4uWQ1G1kGoUT4CuOzcA0BugANjhiv3dg4+q2w8s/AoKazhgQfXtBxWtgrUVy6p9Q6vL8UP27x8m/q5siHDUMMcHyKtYXv4bXu6e2lXZp/z3UOXz1XSO9tlrKP/LKwC2l71ZS0aY7SMeP8zPtG9e5Z9p9xq2rXb8UBV/apV+z/EE6LL9O7rwXZhgwkukiJzZCyDKRBGJn0sUbwLTVPXj4Pv5wG9VNWJx4XBKFC1dcTF8+CHMng2vvgo7dlSs69gR9u51SSIpCd5/H0aOhMJC+PxzaNUKTjihYvtPP4Uvv4Rf/9odNzERHn4YhtZwZWrVyh17/HgoKoKEBLf9FVe4fQGWL3cxlNkzaz5n/O+pxFFKKQl8edqtxJ99Fvk9BxJITgVcjGX7P/88/PznLnkkJcH8+dA1o4CDr86lw4ezaffpv0nMyy0//v7MLNrs3YYQoIREFk9+mIwxQ+ndG9LS3DbbtsHOndClC+QtXkH3P/+aBIorbR9J7scrOOGVin223Pgw7U8aysaN7hy9e7vtDhXBN6urbx/pHGlpEFheOaalFzzMmF9UbL9iBZRWfBmv9fiJiTB4cMX2C59YwfEvVY6fIUPJC7nwR/q8ZccfMsT9zgG++w7y8ij/OX/3RnQ/13A/0+Tjh5b/njp3Dm6XCxs3RvczDfd7/vqShxlx6dBKv6e6fObQ31Po7/mzp6rHX9BnaKXfU82fuYQikhq8ROHnRNHkt55amkWLYMECSE6Gr76C11+HPXsq1vfpA+ee6x7HHw+LF7vtx41zF+BozlHXfeq8fUkJDBlSUSSpqmdPGDTIXdUGD3avBw1i+aw17P2/dzniyESO2Lcc3nqLSle1oUPdB/5//w+GDGH5U5+RM3sB7c8dV6d/vOUzFkW1fX32aWnb+zEmv21f331CxWqi+AlwDXA6cALwF1UdUdsxLVHUzaJFFd/eoeJb5cCBFclh2DAQqfEQ3lGFK6+EGTPc+7g49zjpJPf1fu3aig9WF8cdV/Gh+/evfXtjmiFf1lGIyEvAOKCDiGQDvwcSAVT1CeBtXJJYB+QDl3oTafO0YIG7lpaWumvsKafAI49Uvq3gW/fd55JE69Yu6JycykWQkhJYvx5Wr4ZVqyqev/7a3Qcrc9pp8NhjFeV+Y0xYnpYoGoOVKOpm4UI49VSXLMru19f1dpKnnnsOLrrIFXVmzYL//u+67/vJJy4jFhfH2Ic2pvH5skRhvCUCBQXu1v3f/x4j18v58+Gyy9zr6dOjSxIAo0e7WvhoK1qMaeEsUbRQ2dmulcnxx8fI9XL5cpcYSkpcU6pf/ap+xxk5MkY+sDH+YYmihZoyxV13a2rG6CvZ2a4+Yf9+mDQJHnzQ64iMaVFsUMAWLCkJ2rWrfTtP7dsHp58OW7bAmDHwj3+42ndjTJOx/zjjX0VFrsnq8uWu3e7rr7uWTsaYJmWJogXauxc6dHB3c3xL1XWjnj/fdaWdMycGij/GNE+WKFqglStd14Ndu7yOJII77nBNYVNSXO/pXr28jsiYFssSRQu0apV7HjLE2zhqNGMG/OlPEB8Pr7ziek4bYzxjiaIFWrnSPfuyF/ZDD7nhOQAef9xVZBtjPGXNY1ugskThuxLF7Nlw883udUJCzUPNGmOalJUoWiDf3nqaObPitarrQW2M8ZwlihZm7143h0JyMhxxhNfRVJGf757j4lwnj3HjPA3HGONYomhhykoTgwb5rN+aKqxZ415fd50N2GeMj1gdRQvj24rsjRth61bXV+Khh3yWxYxp2ey/sYXxbUX2p5+655EjLUkY4zP2H9nC+LYiuyxRjB7tbRzGmGosUbQwvr319Mkn7nnUKG/jMMZUY3UULUhRERx1FKSm+mxEjAMH3DSlZRNkGGN8xRJFC5KUBO++63UUYSxeDIGAG6ojOdnraIwxVditJ+O9svoJu+1kjC9ZomhBsrPh4EGvowijrH7CKrKN8SVLFC3IpZe6+om5c72OJERpKXz2mXttJQpjfMkSRQuiComJ0Lev15GEWLXKzYV9xBHQvbvX0RhjwrDK7BZk3jwoLnaNi3zDmsUa43t+umSYJpCY6HUEVVhFtjG+Z7eeWohDh9ytJ9+ximxjfM8SRQvx29+68faee87rSEJs3w7r17t5sY86yutojDE1sETRQqxcCbm5kJHhdSQhFi1yzyee6LOKE2NMKEsULYQvBwO0imxjYoIlihYgN9dN9dCmjc/GeLKKbGNiQlSJQkTiReQiEXleRN4TkWOCyzODy60hvA/5cla7wkJYtgxE3K0nY4xv1fmyISLJwIfAM8BE4GQgM7h6PzAN+GU0JxeRCSKyRkTWicgtYdb3FJEPROQLEflaRE6P5vjG8eXQ4suWueFshwzxWcWJMaaqaL5f/gEYDpwD9AGkbIWqlgL/Ak6t68FEJB54DDgNGAycLyJVL2W/A15R1WOAKcBfo4jXBPlyVjubqMiYmBFNopgEzFDV14FAmPXrgF5RHG8EsE5V16tqEfAyrqQSSoG04Ot0YGsUxzdBVpFtjDkc0SSKbsBXEdbnA22jOF53YHPI++zgslB/AH4qItnA28C14Q4kIlNFZKmILN21a1cUIbQMvrv1pGolCmNiSDSJIofqF/JQQ2j4b/znA8+oag/gdOA5EakWs6rOUNXhqjq8Y8eODRxCbPNli6d162DXLujUCfr08ToaY0wtokkU84FLg5XalYhIb+Ay4J0ojrcFyAp53yO4LNTlwCsAqroIaA10iOIcLV7ZbaeBAyE+3ttYyoU2ixWJvK0xxnPRJIo/4lo5LcG1blJggojcB/wHOATcF8XxlgD9RKS3iCThKqvfqLLNJmA8gIgMwiUKu7cUBavINsYcrjonClVdh7tolwB34Vo93QT8FlfXMF5VN9d8hGrHKwGuAd4FVuNaN60UkbtE5KzgZjcCV4jIV8BLwCWqvhzazrc6dYIJE2DMGK8jCWEV2cbEFKnPdVdEhgKDcMniW1X9oqEDq6/hw4fr0qVLvQ7D1CQ3FzIzISnJTVjUqpXXERljABFZpqrDw62r80hsInISsFpVd6nqCmBFlfUdgMGq+tFhRWuat7KBAIcPtyRhTIyIpo7iA+BHEdaPD25jfKKgAD7/HPLyvI4khI3vZEzMiSZR1NY8JZ7wHfGMR774Ak44AU46yetIQlhFtjExJ9oh4iJVaIwCdh9GLKaBFRa6+YCOPdbrSIJKSmDxYvd65EhvYzHG1FnEOgoRuQ64LmTRdBH5U5hNM3FDbTzdgLGZw3TyyfD1115HEeLrr+HgQTjySOjc2etojDF1VFtldi6wMfi6F6539o4q2yiuYvsz4H8aMDbT3FizWGNiUsREoarPAs8CiMgG4BZVrdopzvjUtm3QpYuPOj9bRbYxManOzWNVtXdjBmIaVm4udOsGXbvCli0+SRZWkW1MTPLLfGemga1e7Z59U6LIzoZNmyA93UfD2Bpj6iLaqVBHi8ibIrJLREpEpLTKo6SxAjXR8d3Q4mWliZEjfTQfqzGmLqKZCvUkXIe6E4DFwX0/wA3uJ7gK7ecaIUZTD76brMgqso2JWdF8tbsd2IabtvSS4LJ7VfVEYALQG/hbg0Zn6s23JQpLFMbEnGgSxQjgb6q6i4oe2HEAqjoXV5q4u2HDM/Xlq+HFDx503cTj4lxXcWNMTIkmUbSiYmKhQ8Hn0KlPvwSOa4CYzGHat8+1dGrdGnr7oa3akiVQWgpHHw2pqV5HY4yJUjSJYhtuFjpU9SCuM97QkPU9cHNVGI/5blY7axZrTEyrcz8KXKV16H/6XOAGEdmISzjX4Cq5jcd8ddsJrCLbmBgXTYni78BuEWkTfH8bUAA8gxvj6RDwmwaNztSLr1o8BQIVc1BYicKYmBRNz+z3gPdC3q8Xkf64eShKgY9VdV/Dh2ii5asWT2vWwN690L07ZGV5HY0xph6iufVUTbCuonzsJxHpraobDjsqc1h8VaIIve3kiy7ixphoNUgXWRE5QkSeAr5piOOZ+gsE4IIL4KyzfNLiySqyjYl5tZYoRCQT18GuH7AHeDk4ZzYi0hn4Y3B9Em6oceOhuDi4/36vowhhFdnGxLzaJi7KAhYBXamYCvVmETkLVy/xf7hJiz4C7lbV+Y0Yq4k1u3fD2rWQkOCm2zPGxKTabj39HpckpgNnANcDecBfgNnAJuBkVR1nScIfFi+Gzz5znaE9N3Omey4pgVNPrWj9ZIyJKbUlilOAF1X1RlV9W1X/gpsatR/wNXCiqi5o5BhNFO680w3QOm+e15EA77xT8bqoCBYs8CwUY0z91VZH0RVYWGVZ2fvHVfUQxlcGDYKdO2Ho0Nq3bXStWrnnuDhISoJx4zwNxxhTP7UlikTcraZQZe+3N3w45nBNn+51BCEKCtzzxRfDFVe4oo4xJubUpR+FRrncGGfjRvd8yy3Qv7+3sRhj6q0uiWKaiNwa8j4elyT+JiJVq0xVVY9usOhMVLZscf0oevTwQd+20lLYvNm9th7ZxsS02iqzN+Hmnmgb8kgOLo+rsrwtkNZokZpaPfww9OwJDzzgdSTA9u2utVOnTtCmTe3bG2N8K2KJQlV7NVEcpgGU9W1TP9wULLvtdMQR3sZhjDlsns5yLyITRGSNiKwTkVtq2GayiKwSkZUi8mJTxxgrFi2Czz93r//4Rx90WbBEYUyz4VmiEJF44DHgNNw83OeLyOAq2/QDbgVGq+oQXIc/E8Z771WUJIqLfdBlwRKFMc2GlyWKEcA6VV2vqkXAy8DEKttcATymqnsBVHVnE8cYM0KHFPdFlwVLFMY0G14miu7A5pD32cFlofoD/UXkExH5TEQmhDuQiEwVkaUisnTXrl2NFK6/deninrOyYP58H3RZsERhTLPhaR1FHSTghgsZB5wPPCUiGVU3UtUZqjpcVYd37NixaSP0ia1b3fPxx/sgSYAlCmOaES8TxRYgtIF9j+CyUNnAG6paHJwQaS0ucZgqyhJFt27exgG4yhJLFMY0G14miiVAPxHpLSJJwBRCZssLeg1XmkBEOuBuRa1vwhhjhq8SxZ49bvjatDTIyPA6GmPMYapzohCRtSLyWxHp0hAnVtUS4BrgXWA18IqqrhSRu4LzXRBclyMiq4APgJtVNachzt/clCWKrl29jQOw0oQxzUw0c2YXA/cBd4vI28DfgLdVNVDfk6vq28DbVZbdGfJagV8HHyYCX5UoLFEY06zUuUQR7McwCngW+C/gdWCziPxJRPo2UnymjixRGGMaS1R1FKr6mapegZun4ufABlyHuLUi8r6IXCAirRohTlOLyZPh/PPdgICe27TJPffs6W0cxpgGUa/KbFXNV9WZqjoGGIjrLDcOeA7YKiLTRcSuEk3orrvgxRd9UndsJQpjmpV6t3oSkXgROQd4GDgPN/T4B8BnuErq1SJStae1aQksURjTrESdKERkoIg8iOvzMBsYDjwE9FfVU1T1J7hSxhrADwNeN3s7d8LChZCd7XUkQZYojGlWomkee7mIfAKsxLVC+gqYBPRQ1VtU9buybVV1HfAXwCq5m8D778NJJ8ENN3gdCa7/xO7dbsCpzp29jsYY0wCiaR77FG6e7GnAU6r6fS3br8LVWZhGlpLihu34wQ+8joTKFdlxfh8hxhhTF9Ekiv8G/q2qpXXZWFU/Bz6vV1QmKmee6R6+YLedjGl26vSVT0RSgT/jKqmNqZklCmOanTolClXNA9oDBxs3HFMfu3e7yYp8wRKFMc1ONDeRP8O1cDI+M2oUtGoFa9Z4HQmWKIxphqJJFLcAk0XkUhGRxgrIREfVDd+hWjF5kacsURjT7ERTmf0wsBc3GOADIvIdkF9lG1XV8Q0VnKndgQOuRWqbNm5Ub89ZojCm2YkmUfTB9b4Otn/EGsn7QOhggJ6X84qLXUAi0L3qrLbGmFhV50Shqr0aMQ5TT9u2uWdfjBq7ZQsEAi5JJCV5HY0xpoFYj6gYZ8OLG2MaWzS3ngAQkV7AKbhbTy+o6vfBqUy7ANtVtahhQzSRWKIwxjS2qEoUInI/8C0wA7gLV28B0Bo3ZMdVDRqdqZUlCmNMY4tmUMBfADcDjwE/BsqrTlV1P/AG4JeBJFoMmyvbGNPYoilRXAW8qqrXA1+EWf81MKAhgjJ1ZyUKY0xjiyZR9Afei7B+F9Dh8MIx0fJVqydLFMY0S9EkikIgJcL6I4Dcw4rGRKWsVzb4IFEEAhVDjFuiMKZZiSZRfA6cE26FiLQGfgZ80hBBmbopKoKJE+FHP4K2bT0OZudOOHQI2rd3E2QYY5qNaJrHPgi8KyLPAU8Hl3URkVOBPwI9gAsaOD4TQatW8NJLXkcRZLedjGm2oumZPU9Efgk8QkVCKJvBrgi4QlUXNXB8JlaUJYqePb2NwxjT4KLqcKeqM0TkDdxc2QNxTWS/BV5R1S2NEJ+JYNcuyMtzTWNbt/Y4GCtRGNNsRT2Eh6puV9X/VdWrVfUqVf0fSxLeePpp6NMH7rjD60iwimxjmrFoOtz9S0TOEpGoh/0wjSMxEbKyfHK3x0oUxjRb0ZQoTgNeBbaKyHQROa6RYjJ19Otfuy/y117rdSRYojCmGYsmUXQGrgTWANcCn4vIChG5SUT8MICE8ZIlCmOarTonClXdr6pPqeoPgb64QQFbAQ8Am0RkjohMaaQ4TRilpV5HELRvn3skJ7t+FMaYZqVe81Go6veq+kdV7QeMAf4OjAaej+Y4IjJBRNaIyDoRuSXCdueKiIrI8PrE2xypuk52XbtCYaHHwYSWJjyfZs8Y09AOa+IiEUnBjQHVHze8R52vEiISjxuJ9jRgMHC+iAwOs11b4Dpg8eHE2tzs3QsFBZCfb01jjTGNK+pEIc6PReR5YAeul/ZQ4FHg+CgONQJYp6rrg5MdvQxMDLPd3cD9uLGmTJBvxngCSxTGNHN1buoqIkOBi3C9srsCxcDbwLPAW6paEuW5uwObQ95nAydUOeexQJaqviUiN0eIbSowFaCnL9qKNj5LFMaYphJNn4ivg89LgfuAl1R1T8OH5IhIHPAwcElt26rqDNysewwfPlwbKyY/sURhjGkq0SSKB4BnVXV1A517C5AV8r5HcFmZtrhbWgvEVZB2Ad4QkbNUdWkDxRCzfDmzXQspzRnT0kQzKGCNrZLqaQnQT0R64xLEFEJGn1XVfYRMhCQiC4CbLEk4VqIwxjSVqIfjEJH/ws1L0Se4aD3wL1VdEM1xVLVERK4B3gXigadVdaWI3AUsVdU3oo2tJfHNzHaFhbBjByQk+CAYY0xjiKYyOw5XcX0BrhlsILgqDrhaRF4ALlbVOtcRqOrbuArx0GV31rDtuLoetyXwTYlic7A9Qo8eEB/vbSzGmEYRTfPYG4ELgVnAMKBN8DEMeCW47tcNG56piW8Shd12MqbZi+bW0yXAXFU9r8ryr3Gd5TKBy4A/N1BspgaBQMWtJ88rsy1RGNPsRVOi6AP8O8L6f1NRb2EaUU4OFBdDRga0aeNxMJYojGn2oilRHMSNIFuTLsFtTCMrKIAxYyA11etIsERhTAsQTaJYCFwjIv+nqitDVwTHaLoaWNCAsZka9OwJCxd6HUWQJQpjmr1oEsWdwGfAFyLyOrAquHwIcCZQBPy+YcMzvmeJwphmL5oOd8tFZCzwCHBu8FHmU+A6VV3ewPGZMPbvh1at3MNTpaWQne1eZ2VF3tYYE7OiGj1WVZeq6mhcXcWJwUdnVR2jqssaI0BT3S23uKHFn3zS40C2boWSEujSxQdjnRtjGkvUPbMBVHUXsKuBYzF1dOiQ6wjdoUPt2zYqG+PJmBahziUKETlPRP4RYf2zIvL/GiYsE8nf/+6SxcRws3c0JaufMKZFiObW0zVUDNsRTilw7eGFY+oqLs6VKjy1aZN7tkRhTLMWTaIYBHwRYf0XuClNTUthJQpjWoRoEkUKrtRQE8XNIWEa0Y4dbnynH//Y60iwRGFMCxFNotgAjImwfgyw6fDCMbXZutWN87Rjh9eRYInCmBYimkTxKjBJRC6vukJELgMmAf9qqMBMeL4ZNVbVEoUxLUQ01aHTgInADBG5AfgyuPxoXN3EGuDeBo3OVOObKVBzciA/H9LT3cMY02zVuUShqgeA0cCTQFfcBEYXAN2Ax4FRqrq/MYI0FXxTorDShDEtRlQNLIPzWF8lIldTMZ/17mhmtTOHxzdToFqiMKbFqG/PbMV6ZnvCShTGmKYW1VhPItJWRO4UkY9F5FsRGRlc3iG4fGDjhGnK+C5R2PAdxjR7dS5RiEhH4GPcLHbrgs9tAFR1t4hcDGRg82Y3Kt8lCitRGNPsRXPr6R7cLHYn4PpL7Kyy/nVgfAPFZcIoKXH9J0Sgc6S5BpuCJQpjWoxobj2dAfxVVf+D64Vd1XrAJiVoRDt3QiAAHTtCYqLHwViiMKbFiCZRdMDdcqpJALBJCRqRb1o85eXBnj1u5qROnTwOxhjT2KJJFNuBvhHWH4MN4dGoCgqgd2/38FTZqLE9e7phbI0xzVo0dRRvA5eLyP/i5scuJyInABcB0xsuNFPVmDGwfr3XUWC3nYxpYaL5OvhHoAQ3nPh9uHqKi0XkJeAjYCtwf4NHaPzHEoUxLUo0Q3hsx82RvRi4DBDgZ8BkYC7wQ1Xd0xhBGicQadqopmSJwpgWJaobzKq6WVUnAu1wzWRPBDqq6pmqmt0YAZoKEye6Fk8ffOBxIJYojGlR6lUTqar7VXWJqn5eVooQkdEiMr9hwzOhtm2D3bshNdXjQCxRGNOi1ClRiEh7ERkhIkeGWXeiiMzF1VOcFM3JRWSCiKwRkXUickuY9b8WkVUi8rWIzBeRFn1lWrzYJYujj/Y4EEsUxrQoEROFiMSLyBPADmARsEZEPhWRTiKSJiIvAp8A/wW8CBxV1xOLSDzwGHAabj6L80Wk6pzbXwDDVfUHwCzggboevzmKj4cuXSApycMgiorcOCJxcdC9u4eBGGOaSm3NY68FpgLZwGfAkbh6iceAHsAI4DngblX9LspzjwDWqep6ABF5GTcx0qqyDVQ19G78Z8BPozyHaWjZ2W52u+7dfdA93BjTFGq79fQzYDkwUFUnq+qxuEmKzsUljTGqekk9kgRAd2BzyPvs4LKaXA7MCbdCRKaKyFIRWbprV/Mc/fzLL+H44+HGGz0OxG47GdPi1JYo+gP/UNX8kGWPB5/vV9VFjRNWZSLyU2A48GC49ao6Q1WHq+rwjh07NkVITe7772HpUli71uNALFEY0+LUduspBTd0R6iy98sP89xbqDyIYI/gskpE5BTgdmCsqh46zHPGLN8ML142fEcLThSBQIDs7GwOHjzodSjG1FliYiKdOnUiLS0t6n3rMoRH1ZFiy94XR322ypYA/USkNy5BTMHNwV1ORI7BzdE9QVWrDmveovhmQEArUbB7925EhAEDBhBnY12ZGKCqFBQUsGWL+y4ebbKoS6I4XUS6hLxPxiWLSSIyrHo8+j91ObGqlojINcC7QDzwtKquFJG7gKWq+gbuVlMq8E8RAdikqmfV5fjNjW9KFJYoyM3NpVevXpYkTMwQEZKTk+nevTtbt25tlERxAVW+6Qf9IswyBeqUKABU9W3cYIOhy+4MeX1KXY/V3Fmi8I/S0lISrcWXiUFt2rShuDj6m0G1JYr/ql84pqGVJYquXT0MIhCoPMR4CxYs4RoTU+r7dxsxUajqh/U6qmlwvihR7NjhOtx16AApKR4GYoxpSnaTNQYUFbkxnuLj3aCAnrHbTqaF++STT+jXrx+pqam89tprXofTZCxRxIDtwQbJXbq4ZOEZSxS+9/HHHzNq1CjS09Np164do0ePZsmSJV6H1WQWLFhAjx49Gu34d955J9dccw15eXmcffbZYbd58cUXGT58OKmpqXTt2pXTTjuNjz/+GHANIS677DK6dOlC27Zt6d+/P9OmTWu0eBtKNDPcGY/44rYTwEcfuecE+7Pxo/3793PGGWfw+OOPM3nyZIqKili4cCGtWrXyOrQGoaqoqqetzTZu3MiQIUNqXP/www8zbdo0nnjiCU499VSSkpJ45513eP311xkzZgw33HADBw8eZPXq1aSnp7N27VpWrFjRhJ+gnsp++M3lcdxxx2lzM3u2KqhOnOhhEJ9+qpqQ4AJJTHTvW6hVq1bVe99PP1W9997G+fEtWbJE09PTa1z/+9//Xi+88MLy9xs2bFBAi4uLVVU1JydHL7nkEu3atatmZGToxJA/uNdee02PPvpobdu2rfbp00fnzJmjqqq5ubl62WWXaZcuXbRbt256++23a0lJiaqqfvvtt3rSSSdpWlqatm/fXidPnqyqqoFAQK+//nrt2LGjtm3bVocOHarLly8PG/PYsWP1tttu01GjRmnr1q3122+/1aeffloHDhyoqamp2rt3b33iiSdUVTUvL09bt26tIqIpKSmakpKiW7Zs0dLSUr3vvvu0T58+2q5dO500aZLm5OTU+HOaMWOG9u3bVzMzM/XMM8/ULVu2qKpqnz59VES0devWmpKSooWFhZX2y83N1ZSUFH3llVdqPPaQIUP01VdfrXF9U6jp7xfXLSHsddXzC3tDP5pjopg5U7VVK9Urr/QwiLvvdn8uoBoX5652LVTVf7SyH0s0jzZtKpJF2bJQxx4bfnkk+/bt03bt2ulFF12kb7/9tu7Zs6fS+toSxemnn66TJ0/WPXv2aFFRkS5YsEBVVRcvXqxpaWk6d+5cLS0t1ezsbF29erWqqp599tk6depUzcvL0x07dujxxx9ffuGeMmWK3nPPPVpaWqoFBQW6cOFCVVV955139Nhjj9W9e/dqIBDQVatW6datW8N+prFjx2pWVpauWLFCi4uLtaioSN98801dt26dBgIBXbBggbZp00aXLVumqqoffPCBdu/evdIxpk+frieccIJu3rxZCwsLderUqTplypSw55s/f762b99ely1bpoWFhXrNNdfoD3/4w/L1RxxxhL733nth950zZ47Gx8eX/zzDufzyy3Xw4MH69NNP69q1a2vcrjFZomimiUJVNRBQrfIFpmndcUf4q1wL1BCJIj6+Itc2VKIoi+3iiy/W7t27a3x8vJ555pm6fft2VY2cKLZu3aoiUi25qKpOnTpVr7/++mrLt2/frklJSZqfn1++7MUXX9Rx48apqurPfvYzveKKK3Tz5s2V9ps/f77269dPFy1apKWlpRE/z9ixY/WOO+6IuM3EiRN1+vTpqho+UQwcOFDnzZtX/n7r1q2akJAQ9oJ+2WWX6c0331z+/sCBA5qQkKAbNmxQ1ciJ4vnnn9fOnTtHjDU/P1//9Kc/6bHHHqsJCQnat29fffvttyPu09DqkyisMjtGiIBnt5oDAZg1y70+91yYPx9GjvQoGP+pa3r49FNo08Y1SEhKgnHjKu8fatmy8MtrM2jQIJ555hmys7NZsWIFW7du5frrr691v82bN9OuXTsyMzPDruvbt2+15Rs3bqS4uJiuXbuSkZFBRkYGv/jFL9i5042288ADD6CqjBgxgiFDhvD0008DcPLJJ3PNNddw9dVX06lTJ6ZOncr+/ftrjC0rK6vS+zlz5nDiiSfSrl07MjIyePvtt9m9e3eN+2/cuJFzzjmnPMZBgwYRHx/Pjh07qm27detWjghprJGamkr79u3Lh76IpH379uzevZuSkpIat2nTpg233XYby5YtIycnh8mTJzNp0iT27NlT6/G9ZInC1O6NN2D1aujRA1580ZJEPY0c6XLs3Xc3Ta4dOHAgl1xySXllaUpKCvn5FQNBb99eMd5nVlYWe/bsITc3t9pxsrKy+O676jMJZGVl0apVK3bv3k1ubi65ubns37+flStXAtClSxeeeuoptm7dypNPPslVV13FunXrAPjVr37FsmXLWLVqFWvXruXBB8MODA1U7iR26NAhzj33XG666SZ27NhBbm4up59+urs9QvgOZVlZWcyZM6c8xtzcXAoLC+keZuKtbt26sbGsdR9w8OBBcnJywm5b1ciRI2nVqlWdm82mpaVx2223cfDgQTZs2FCnfbxiiSIGjB0Lw4bB+vUenFwV7rvPvb75Zo+n14t9I0fCrbc2TpL45ptv+POf/0x2djbgSgIvvfQSJ554IgDDhg3jo48+YtOmTezbt4/7yn6vUN6M86qrrmLv3r0UFxfzUbCV2+WXX87MmTOZP38+gUCALVu28M0339C1a1d+/OMfc+ONN7J//34CgQDfffcdH37o+un+85//LI8lMzMTESEuLo4lS5awePFiiouLSUlJoXXr1nVuyVRUVMShQ4fo2LEjCQkJzJkzh7lz55av79y5Mzk5Oezbt6982ZVXXsntt99engB27drF66+/Hvb4559/PjNnzuTLL7/k0KFD3HbbbZxwwgn06tWr1tjS09O56667uPrqq3nttdfIz8+nuLiYOXPm8Jvf/AaAu+++myVLllBUVERhYSGPPPIIGRkZDBgwoE6f3zM13ZOK1UdzrKPIzHQ3IXbu9ODk8+a5k3fooHrwoAcB+M/htHpqTNnZ2Tpp0iTt1q2bJicna7du3XTq1Km6b9++8m2uuuoqTU9P1759++qMGTOqtXq66KKLtFOnTpqRkaHnnHNO+X7/+te/9KijjtLU1FTt27evvvPOO6rqWvpceeWV2r17d01LS9Nhw4bpSy+9pKqqN998s3br1k1TUlK0T58++uSTT6qq6rx58/Soo47SlJQUbd++vV5wwQV64MCBsJ9p7Nix+tRTT1Va9uijj2qnTp00PT1df/rTn+p5552nt99+e/n6Sy+9VNu1a6fp6enlrZ7+/Oc/a//+/TU1NVX79Omjt956a40/x8cff1z79OmjmZmZ+pOf/KRSHUukOooyzz//vB533HGanJysnTt31tNPP10/+eQTVVW9++67dciQIdq2bVvNzMzUsWPHlq9rKvWpoxCN9iaozw0fPlyXLl3qdRgNats215fimGPcVNVN6pRTKu6X/O53TXxyf1q9ejWDBg3yOgxj6qWmv18RWaaqw8PtYz2nYkDXrh4NBrhkiUsSqalw9dUeBGCM8QOrozA1K7uHfdVVEKY1jDGmZbBE4XPvvw9nnw0zZjTxiVetgldfdW1yb7ihiU9ujPETu/XkcytWwOuvu5apTer++93zpZe60QiNMS2WlSh8zpMBAb//Hl54wfUMu/nmJjyxMcaPLFH4nCeJ4qGHoLQUpkyBPn2a8MTGGD+yROFzTT4F6o4d8Pe/u9e33NJEJzXG+JklCp9r8hLFI49AYSGcdRYMHdpEJzXG+JklCp/bts09N0mi2LcPHnvMvb711iY4oTGxxaZCNb6Tnw+5uW54pXbtmuCEf/0r7N/vhjUNjg9kYotNhertVKi9evWiTZs2pKam0rlzZy655BLy8vLK17/55puMGDGClJQU2rdvz4UXXlg+HlYkn3/+OaeffjoZGRm0a9eOESNGMHPmzPL19957L7179yY1NZUePXpw3nnnNcjnLWOJwsdCSxNhBsVsWAUFMH26e33bbY18MtMYyqZCvfbaa9mzZw9btmzh97//fbOaCjUQCHgaQ21ToQL8+9//Ji8vj//85z8sXbqUe+65B4BZs2ZxwQUXcP3117N7925WrlxJq1atGDNmDHv37q3xeIsWLeLkk09m7NixrFu3jpycHB5//HHmzJkDwLPPPstzzz3HvHnzyMvLY+nSpYwfP77hPjTYoIB+9tFHbjy+UaOa4GSPPupOdtxxbpYkU6PDGhSwEedCtalQvZ0KVbX6oIE33XST/uQnP9FAIKA9e/bU+++/v9L2paWlOmTIkIiTM40ePVqvuuqqGtdfffXVet1119W4viqb4a6ZJYqXX3a/oXPPbeQTFRWp9uzpTjZrViOfLPZV+kerz/R29XnUgU2F6u1UqFXXb9q0SQcPHqy/+93vdPXq1Qro+vXrq+1z55136oknnhj2eAcPHtS4uDh9//33azznc889p5mZmfrAAw/okiVLyhN1TSxRNLNE8fDD7jd07bWNfKJnn3UnGjBAtZapKY1/E0VZbDYVqjdToZatT0lJ0fT0dO3Zs6f+8pe/1Pz8fF24cKECWlBQUG2fxx9/XI888siwx8vOzlagPDHX5Pnnn9fx48drcnKytmvXTqdNm1bjtjYVajNTNgFZo7Z4CgRg2jT3+pZbPBjHPMZFc7kPnQu1TRv3vq771pFNhVpdU02FWua1114jNzeXjRs38te//pU2bdrQoUMHALaVVTyG2LZtW/n6e++9l9TUVFJTU7nyyivJzMwkLi4u7H6hLrzwQubNm0dubi5PPPEEd9xxB++++26dY66NXRV8bNo02LULrryyEU9SNs1pVhZccEEjnsg09VyoNhVqRZxNMRVqJAMGDKBHjx7885//rLQ8EAgwe/bs8srn2267jby8PPLy8njiiSdITk5m5MiRzJ49u07nSUxMZNKkSfzgBz8o/703BEsUPiYCHTpARkYjnUAV7r3Xvb7pJpvmtCk04lyoNhWqt1OhRiIiPPTQQ9xzzz28+OKLFBYWsn37dn7+85+zf/9+bogwQvMDDzzAM888w4MPPkhOTg4AX331FVOmTAHgmWee4a233uLAgQMEAgHmzJnDypUrOeGEEw4r5kpquicVq4/mVEfR6P7yF3djIz3dpjmNgk2FalOh1qS29a+99poOHz5ck5OTNTMzU6dMmaKbNm2qcfsyixcv1gkTJmhaWppmZmbqiBEj9Nlnn1VV1dmzZ+uoUaM0IyOjvBXZzJkzazxWzE2FKiITgEeAeOBvqjqtyvpWwD+A44Ac4DxV/T7SMZvTVKijR7vSxKxZ7pZ2vRUXw3ffuTkmVq92z0uWwLffuvWJifDhh41+K6S5sKlQTSyLqalQRSQeeAz4EZANLBGRN1R1VchmlwN7VfVIEZkC3A80bJfDoFWPvMeuWR+SfMooek8+vtbtN7yyhLy5n9LpnFEMudhtv38/FBVB27Zuvh+AgwddX7bQfVJ/HPkcbdu6a3vap0sYLZ+y4bFRDL649pj4/HOYOxfat4eSkoqk8O237oA1CQRgwQJLFMaYsLycuGgEsE5V1wOIyMvARCA0UUwE/hB8PQt4VEREG7gYtHzGIgZfP4HBBODjkDNG0KHsxadAcMqGtDDbpQQf1fap5RytgDkAGjz+4U4L0asXDB4Mgwa555ISuP56l9mSktywHcYYE4aXiaI7sDnkfTZQtfalfBtVLRGRfUB7oFJbOBGZCkwF6NmzZ9SB5MxegOCGBlDgIMkUSnKN27fWfFLIR4LbFycmk5SezP4D7rqblgZJiW7b/AI3ZlPVfSKdIy0NyM8nsbj6OWqUn+8eZcaNg5//3CWGAQMgJaX6Pkcd5UoS48ZZacIYU6NmMRWqqs4AZoCro4h2//bnjqNwbhsSKaKYJDY8OY+jptZ84Vw+YxF9fzG+fPvvHnXbhytRJAcfVfep7zlqtGgRjB9fUUK4997aL/4jR1qCqCdVDdsU0xg/q+/NGC8TxRYgtCdNj+CycNtki0gCkI6r1G5QR00dyXLmkzN7Ae3PHRf5glyP7ZvkHGVt9K2E0Oji4+MpLi4myZoTmxhTUFBAYmJi1Pt51uopeOFfC4zHJYQlwAWqujJkm6uBo1T1ymBl9n+r6uRIx21OrZ6MP+3cuZNDhw7RvXv3Orf/N8ZLqkpBQQFbtmyhc+fOpKVVv//hy1ZPwTqHa4B3cc1jn1bVlSJyF6497xvA34HnRGQdsAeY4lW8xpTp0KED2dnZrFmzxutQjKmzxMTEGpNEbTztR9EYrERhjDHRi1SisHKzMcaYiCxRGGOMicgShTHGmIgsURhjjImo2VVmi8guYGOtG1bXgSo9vmNQrH+GWI8f7DP4Rax/Bi/iP0JVO4Zb0ewSRX2JyNKaavxjRax/hliPH+wz+EWsfwa/xW+3nowxxkRkicIYY0xEligqzPA6gAYQ658h1uMH+wx+EeufwVfxWx2FMcaYiKxEYYwxJiJLFMYYYyKyRAGIyAQRWSMi60TkFq/jiYaIZInIByKySkRWish1XsdUXyISLyJfiMibXsdSHyKSISKzROQbEVktIjE1KYiI3BD8G1ohIi+JSGuvY6qNiDwtIjtFZEXIsnYi8p6IfBt8zvQyxtrU8BkeDP4dfS0ir4pIhochWqIQkXjgMeA0YDBwvogM9jaqqJQAN6rqYOBE4OoYiz/UdcBqr4M4DI8A76jqQOBoYuiziEh34FfAcFUdihv6PxaG9X8GmFBl2S3AfFXtB8wPvvezZ6j+Gd4DhqrqD3Dz9tza1EGFavGJAhgBrFPV9apaBLwMTPQ4pjpT1W2q+p/g6wO4i1N3b6OKnoj0AH4C/M3rWOpDRNKBk3BzqKCqRaqa62lQ0UsA2gQnFUsGtnocT61U9SPcXDWhJgLPBl8/C5zdlDFFK9xnUNW5qloSfPsZbgZQz1iicBfVzSHvs4nBCy2AiPQCjgEWexxKfUwHfgMEPI6jvnoDu4CZwdtnfxORFK+DqitV3QI8BGwCtgH7VHWut1HVW2dV3RZ8vR3o7GUwDeAyYI6XAViiaCZEJBWYDVyvqvu9jicaInIGsFNVl3kdy2FIAI4FHlfVY4CD+P+WR7ngffyJuITXDUgRkZ96G9XhU9f+P2b7AIjI7bjbyy94GYclCjdfd1bI+x7BZTFDRBJxSeIFVf2X1/HUw2jgLBH5Hnfr72QRed7bkKKWDWSrallpbhYuccSKU4ANqrpLVYuBfwGjPI6pvnaISFeA4PNOj+OpFxG5BDgDuFA97vBmiQKWAP1EpLeIJOEq8N7wOKY6ExHB3RdfraoPex1PfajqraraQ1V74X7+76tqTH2bVdXtwGYRGRBcNB5Y5WFI0doEnCgiycG/qfHEUGV8FW8AFwdfXwy87mEs9SIiE3C3Ys9S1Xyv42nxiSJYYXQN8C7uH+MVVV3pbVRRGQ38DPct/Mvg43Svg2qhrgVeEJGvgWHAvd6GU3fBktAs4D/Acty1wVfDSIQjIi8Bi4ABIpItIpcD04Afici3uJLSNC9jrE0Nn+FRoC3wXvB/+glPY7QhPIwxxkTS4ksUxhhjIrNEYYwxJiJLFMYYYyKyRGGMMSYiSxTGGGMiskRhjDEmIksUxnhIRB4VkS3BTm7G+JIlCmM8EkwOZwOvez1EgzGRWKIwxjvH40Yqfs3jOIyJyBKFMVESkftFREWkv4j8JXjrKD84m1pWcJuficiy4PI1InJ2mEOdA+QCH4QcO01Ebg/ObLZPRPYHZy98rEk+nDFh2BAexkRJROYCY4D1wFfAQuA44OfAv3ED7P0QeB43xPUtQBugl6ruCjnON8AyVb0w+L4V8DlwBDATN6hgMnAU0FNVf9wUn8+YqhK8DsCYGHQM7sJ/r6q+WLZQRIYDZwJv4aYULQ4uL8ZNzPQD3NSciMggYADwu5DjnhXc5tQYnjTINEN268mYKASnbO0AvBGaJIL2AsXAZWVJIqhsIqmSkGVnA4XAOyHLMoPPI0TE/jeNb9gfozHROSb4/H9h1g0FFqhq1YlyBgWf14QsOweYp6p5Ictm4W5l3Q1sFZGnROQMSxrGa/YHaEx0ymat+yx0YbASu2PV5SH7bA1OblRWKhlOldZOqroHV9cxAZeIfoSr8/g4OKmWMZ6wRGFMdI4BclV1fZXlZQnkPzXsE7r8bFwld7WZFFW1VFXfVdXrgL64CvGRwNGHGbcx9WaJwpjoVL3olzku+LwsdKGI9ALaVdnnHODjKi2gOlbtna2qpUApLqnE1DzupnmxVk/G1JGItAN6Aq+EWX0ssEtVN4dZDsFEETzGSbj5kEM9BIwRkdeBdbgvcacCZwAPqurWBvkQxtSDJQpj6q6sIrumEsUXddjnDNz/3atVtpuPa001GVfXsQfXj+JsVX39MGI25rBZhztjmpCIvAr0VtVhXsdiTF1ZicKYprUIeMrrIIyJhpUojDHGRGStnowxxkRkicIYY0xEliiMMcZEZInCGGNMRJYojDHGRGSJwhhjTESWKIwxxkT0/wGbUAn9nwTnPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "success_score = 60\n",
    "cslin_success = np.sum(scores_lin>success_score,axis=1)/nb_trials\n",
    "cspo_success = np.sum(scores_po>success_score,axis=1)/nb_trials\n",
    "\n",
    "mlin_tr = np.where(cslin_success > 0.5)[0][0]\n",
    "mpo_tr = np.where(cspo_success > 0.5)[0][0]\n",
    "\n",
    "print(f\"m/k: CS {m_vec[mlin_tr]/k}; PO-CS {m_vec[mpo_tr]/k}\")\n",
    "print(f\"PO-CS/CS: {m_vec[mpo_tr]/m_vec[mlin_tr]}\")\n",
    "print(f\"mCS = {m_vec[mlin_tr]}, mPO = {m_vec[mpo_tr]}, 2*mCS - 2 = {2*m_vec[mlin_tr]-2}\")\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(m_vec/k, cspo_success*0 + 1.0, '--', color='gray')\n",
    "plt.plot(m_vec/k, cslin_success, 'b.-.',linewidth=2, label='Success rate of CS')\n",
    "plt.plot(m_vec/k, cspo_success, 'r.-',linewidth=2, label='Success rate of PO-CS')\n",
    "\n",
    "plt.legend(fontsize=12)\n",
    "#plt.title('Basis Pursuit Recovery Rates')\n",
    "plt.xlabel('$m/s$',fontsize=18)\n",
    "#plt.xlim((0,70/k))\n",
    "plt.ylabel('Recovery Rate',fontsize=18)\n",
    "#plt.savefig('BP_rec_rate_Gaussian.pdf', bbox_inches='tight', pad_inches=0.025, transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Noisy PO-CS reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now consider a noisy PO-CS model described by \n",
    "\n",
    "$$\n",
    "\\boldsymbol z = {\\rm sign}_{\\mathbb C}(\\boldsymbol \\Phi \\boldsymbol x) + \\boldsymbol \\xi,\n",
    "\\tag{PO-CS}\n",
    "$$\n",
    "\n",
    "with $\\boldsymbol \\xi \\in \\mathbb C^m$ an additive noise with bounded amplitude, i.e., such that $\\|\\boldsymbol \\xi\\|_\\infty \\leq \\tau$.\n",
    "\n",
    "Note that $\\tau < 2$ since there exists a counterexample of additive noise with $\\|\\boldsymbol \\xi\\|_\\infty \\leq 2$ that totally remove any information about the signal in the phase (that is, the random dephasing, see the paper)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cef0731e9f96487a9e6d16e48c23e810",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/23 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sensing model\n",
    "sensing_model = 'gaussian'\n",
    "\n",
    "# Problem dimensions\n",
    "n = 100\n",
    "k = 10\n",
    "m_vec = np.arange(1,70,3)\n",
    "\n",
    "\n",
    "# Levels of noise\n",
    "nb_tau = 16\n",
    "tau_vec = 2*np.logspace(-3, -1, nb_tau)\n",
    "\n",
    "# Number of trials\n",
    "nb_trials = 100\n",
    "\n",
    "# Inits\n",
    "scores_po = np.zeros((m_vec.size,nb_tau,nb_trials))\n",
    "\n",
    "# Main Monte Carlo Loops\n",
    "pbar = tqdm(total=m_vec.size) #comment this if it doesn't work for you\n",
    "for m_ind, m in enumerate(m_vec):\n",
    "    \n",
    "    # The measurement vector e_1 such that A_z x0 = e_1\n",
    "    e_1 = generate_canonic_vec(m+2,0)\n",
    "    \n",
    "    for s in np.arange(nb_trials):\n",
    "        \n",
    "        for t, tau in enumerate(tau_vec):\n",
    "                        \n",
    "            # Generating the complex uniform noise with tau variations\n",
    "            rho = np.sqrt(np.random.uniform(0, 1, m))\n",
    "            phi = np.random.uniform(0, 2*np.pi, m)\n",
    "            xi = tau * rho * (np.cos(phi) + np.sin(phi)*1j)\n",
    "            \n",
    "            # Generating the sensing matrix according to the selected model\n",
    "            A = generate_sensing_mtx(m, n, model = sensing_model)\n",
    "\n",
    "            # Generating a k-sparse signal, random k-length support, and Gaussian entries\n",
    "            x = generate_sparse_sig(k,n)\n",
    "\n",
    "            # linear measurements\n",
    "            y = A @ x\n",
    "            \n",
    "            # phase only measurements\n",
    "            z = signc(y) + xi\n",
    "\n",
    "            # Let x0 such that ||Ax0||_1 = kappa sqrt(m) (see [1])\n",
    "            # This is the signal we hope to estimate, as the amplitude of x is lost (see [1])\n",
    "            x0 = x * kappa * np.sqrt(m) / np.linalg.norm(y,1)\n",
    "\n",
    "            # Definition of the matrix A_z for sparse signal recovery in PO-CS, such that A_z x0 = (1, 0, ..., 0)  \n",
    "            A_z = generate_Az(A,z)\n",
    "\n",
    "            # Defining A_xi to compute varepsilon, the level of noise\n",
    "            A_xi = generate_Az(A,xi)\n",
    "            varepsilon = np.linalg.norm(A_xi @ x0)\n",
    "            \n",
    "            # Phase-only model\n",
    "            x_pocs = bpdnalg(A_z, e_1, varepsilon)\n",
    "\n",
    "            # Recording SNRs\n",
    "            #scores_lin[m_ind,t,s] = snr(x_lin,x0)\n",
    "            scores_po[m_ind,t,s] = snr(x_pocs,x0)\n",
    "            \n",
    "    pbar.update(1) #comment this if it doesn't work for you\n",
    "            \n",
    "pbar.close() #comment this if it doesn't work for you"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-42-bd2bab08cf7d>:23: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  cbar.ax.set_yticklabels(['3.0','','','','2.0','','','','1.0'])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAESCAYAAADJ+2ORAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACX40lEQVR4nOydd1xUx/rGv7OF3qsgoiB2ERXUqIkxaupNYq4x/abf9N7bvb/Ue5Pc9N57M4npPbHEqEmMgiKKKCCI9M5Sdpct8/vjLL0tsMAx2YfPfOacs2fnPGfOcJ6dmXfeV0gpccMNN9xww42hhmakCbjhhhtuuPHXgFtw3HDDDTfcGBa4BccNN9xww41hgVtw3HDDDTfcGBa4BccNN9xww41hgVtw3HDDDTfcGBa4BccNN9xwQ6UQQngJIf4QQqQLIXYLIe7t5hxPIcSHQogcIcQWIcS4EaDqFNyC44YbbrihXpiBJVLKJGAmcJwQ4rBO51wM1EgpE4AngIeHl6LzcAuOG2644YZKIRU0OHb1jtR5tf5y4C3H9mpgqRBCDBPFfkE30gSGA2FhYXLcuHEjTcMNN9w4BJCamloppQwfTBlCJEhocvLskt2Aqd2Bl6WUL7eVJbRAKpAAPCel3NKpgNHAQQAppVUIUQeEApUD5T9U+EsIzrhx49i2bVuv5+Tm5jJ+/PhhYuQc1MZJbXzAzclZuDk5h9zcXBISEg4MviQjcJWT5/7LJKVM6elTKaUNmCmECAI+E0JMl1LuGjzH4cdfQnCcQUhIyEhT6AK1cVIbH3BzchZuTs7BdZw0gLeLylIgpawVQqwHjgPaC04RMAYoFELogECgyqUXdxHcczgONDU52/0dPqiNk9r4gJuTs3Bzcg6u4yRom3LpK/VSihDhjp4NQghv4Gggq9NpXwLnO7ZXAuukSr0yu3s4Dmg06tNetXFSGx9wc3IWbk7OwXWcBC56vUYBbznmcTTAR1LKr4UQ9wHbpJRfAq8B7wghcoBq4ExXXHgo4BYcB/T63n9pjATUxkltfMDNyVm4OTkH13Fq6eEMDlLKncCsbo7/X7ttE3DaoC82DFDfT4wRQkNDQ98nDTPUxkltfMDNyVm4OTkH13Fq6eE4k/46+GvdbS8ICwsbaQpdoDZOauMDbk7Ows3JObiOk2t6OH82uHs4DhQWFo40hS5QGye18QE3J2fh5uQcXMepxUrNmfTXgbuH40BCQsJIU+gCtXFSGx9wc3IWfwZOUkJVNeTkg6Eemi1gsUJzs7LdmlvAYmk7NnManHL80HDqGe4eTndwC44Du3fvJikpaaRpdIDaOKmND7g5OYtDiVN1DWTnKcKSvV/Zbkm1df2/zkVnOi84u3fv7v8FeoT79doZQqXm2i5FSkqK7MvTgBtuuDH0MBqhpByKS6G4zJFKoagUcg8oolJd03a+EBA7GibEwYR4JU8YByFB4OEBHvq2XK/r5pheKaM/EEKk9rby37kyJkp43smzjx709Q4VuCXYgdTUVJKTk0eaRgeojZPa+ICbk7MYak4tw10Hi6GwpC1vEZQWcamp7fpdT0+IjoS4WDjtRIe4OAQmPlb5fLiQmprqopJctg7nTwV3D8cNN9xwCiVlsH0XHCjsKiyFJWAydTxfp4OoCIgepQhKdGT328FB/e+FDCVc08OZIuFNJ88+zN3D+avhr/irtL9QGx9wc3IW/eVUUwvb0mFrOmzdoeRFJW2f63QwehTEREHKDDjlWBgTDTHRyrEx0RARBlpt75xCgtVXT66B22igO7h7OG648RdHU5PSc9maDn9sV/KcvLbPJ8TDnCSYM1MRl/Hj+haTQxmu6eFMk/CBk2cnuXs4fzVkZGSQmJg40jQ6QG2c1MYH3JycRdr2DPwCEtssvtpZfx0oVOZgQOmdzJmpWHbNSYLkGcqQ11BgoPUkJdQ1QWUDVNY7UgNUGLoeq6yH5bPhf056F8vIyOg3n+7h7uF0B7fgODBx4sSRptAFauOkNj7g5tQd8grg+/WQldMmKnkF07HZ2s4JClQm5hekwAWnw+xERWhGRQwfz5Z6khKqGqCsDirqlVTpyCsMbcda9isbwGrrvkwPHYT5Q7g/hPnBrLEwKar/nAYPt9FAd3DXiAMFBQVMmDBhpGl0gNo4qY0PuDm1IP8gfPwVfPy1MucC4OsDE+Nh9nQ4+vAa5iWHMCFOORYaMnwT9c1WKKyGgqq2dKASsg5aqTR6UlAFTc3dfzfIRxGP8ACID4e58Y59/3bC0pL8wM9rcPdVUFAw8C93gLuH0x3cguNAZGTkSFPoArVxUhsf+GtzOlCoiMxHX7WJzJyZ8L9/wYoTIH5s28vXYNAREDA0PFp6KDllkFuu5C3bB6qgpLZtyK4FEQEQE+zFtBg4fgaMCYWooDYxCQ9QBEQ/zG8o1z07wV/NbY0zcAuOA7W1tQQM1X/kAKE2TmrjA389TgcKYfXXisj8sV05lpIED98Fp52krGUZKk4NJthxAPaVthOWcmW7rl3cMiEgJhjGR8Kx0yE2VEljw5Q8JgS8PaCgoIjY2B4IjxBqa2tdVJK7h9Md3ILjgJeX10hT6AK1cVIbHzi0OZlMsCcb6huhoXNq6rqff7CtJ5M8Ax66UxGZ+LGu49QCswV2HoSt+2FrnpJnFrf1VLQaiAuH8REwPwESIhSBSYhQjnt5uJ5Tf9Fgga1V4K+DFCedQLuOk3sOpzu4a8QNN4YJUkLmPvhxA/zwM2z4vetiyfbQasHfD/x8lPmYsBB48A5FZMaPcx0vmx2yituEZWsepBcocy+gDHHNiYeVcyAlDqZEKz2V4R7u6g1SQn4D/FoBv1Uo+c4asEk4axy8v2i4Gbl7ON1BRU1mZGHq7T9/hKA2TmrjA+rnVFkFazYqIvPjL22LJycnwKXnwBHzFL9gfr5K8vVRBMbPV/EH1p8J8MJq+HoHfL0dNmV3tOSSMgwhui+v2domLv5ekBwH1x2jTNDPiVfEZSgMDAbz7Ew2SKtShKVFZEqNymd+OpgXBncmwoJwZXs4OHWEu4fTHdw14kBQUNBIU+gCtXFSGx9QHyebDXZnh/LyB0ovJnWn8us7KBCWHQHHHglHHwljYwZ/LbsdtuUpIvPVdtjhMLCKC4cz5ini0QKz2YKnp2eXyXtQhsemx8CcOMWEWDNMUbL6enZmGxxohLx6yGtoS/sbIKMGmu3KeeP9YVmUIi4LwmF6kHJPQ8HJebh7ON1B1YIjhMgH6gEbYJVSpgghQoAPgXFAPnC6lLKmpzKcRVlZmeomn9XGSW18QF2cdmbCxTfBtnRftFo4bDbccxMcc6RiPeaKlfmNZvhplyIw36Qra1c0AhZMgIfPgBNnKkNenXsk2dnDZ6ptaIbceshvVHoiUoIE7I68Zbu0zEhERAASZeiruKmjsBQ3Kee2wEMDY30hzg+unwLzw5UU6UJjsLKyMtcV5kYXqFpwHDhKSlnZbv92YK2U8iEhxO2O/dsGexG1WcuA+jipjQ+og5PJBA88BQ8/B8GB8NL/LJx+kp6gwP6VY7NDdQOUGzoudGzZzi6FDXuVCf1AHzguURGY42dAqH/vZbuynqSEEqMiKu3T/gYlrzQ7W1Ik5LbtCSDGB+IcPZY4v44p2kcR16GEGtrTnxmHguB0xnJgsWP7LeBnXCA4+/btU507ErVxUhsfGHlOm7bAP2+Gvblw3mnw+N1QXJRFUGDvnOx2eGczvLmpbYV9VUPX9SotCPaF0cFw5RI4cRYcMbF/k/YDrScpFSHZXA6bymFLJWQbwNhufkgjINYXxvvBilhliGu8vyISPjpFSDRCyYVo29+7N4spkyYrnwkI9wSPEfbPtm/fPtcV5qw4/vndWbZimEZrBwwJ/CiESBVCXOo4FimlbPFbWwp0u1JLCHGpEGKbEGJbSUkJlZWVlJSUUFRURE1NDbm5uRiNRjIzM7Hb7VgsFqDNW2xaWhp2u53MzEyMRiO5ubnU1NRQVFRES3n5+fk0NDSQlZWF1WolPT29QxkteUZGBmazmezsbAwGAwUFBZSXl1NeXk5BQQEGg4Hs7GzMZnOrL6fU1FQSExNby0hPT8dqtZKVlUVDQwP5+fl93lNaWppL76k9n4HeU/vcFfdksVhG5Dn99vturroDjvg7mMzwzL3ZvPUUFB5MZ8qUKb3e0+oNhcy9x84Fr0CFQRLtU8PKOXDxnGKe/gf89/g8frrVzmcX5bL/f0b2/CuX3P/W8P1VRdyypIQZEZUUFfbvnry8vJx6ThY7vL1pD09mwtLPa4heDQmfwfmbYdV+G5EeVs4eVcfD0xp4f1YVG+eXUnxCDWuScvnqCCPXBWZy81Q7cVVpzA6FxpxUJgVCfU4a8X52zAczidIZsZXlMjc+Em1dEVpDCZ5NlRQfHNr/J2faXlhYPywM+oLWyQRhLe8qR7q0+wIPfajaW7QQYrSUskgIEQH8BFwDfCmlDGp3To2UMri3cpzxFv1ncCk/1FAbHxgZTt+sgctvU6JUXnsxPHCbYlXWF6eiarj9I3j3V4gOhodPh7Pnu36SvrAR0mvAaqe195Cbk82EhAkdehgCJZns8Eel0ov5vRKaHBZr4/zg8AhYGK7kU4NcO6Sl1vaUkpIyeG/RmhSJl5Me6o2D9059qEDVgtMeQoh7gAbgEmCxlLJECBEF/CylnNTbd93hCdxwBcor4bp/w6ovYNokePVROMyJ96WpGR7/Hv77lSICNx0Hd5yk+P0aLGx2yKhVxGJzOWyugILG/pejETAzGBZGOEQmAkb7DJ6fK2CTkmyrZLtFSTutdhqkMjzTIQklF4AG0XpsmYeGa/ycH6tzSXgCXYokwMl3Ts1fR3BUO4cjhPAFNFLKesf2McB9wJfA+cBDjvwLV1xPrb+21MRJbXxgeDhJCe9+AjfcDYYGuPdmuP1qZZ1Mb5ykhM9T4aYPIK8CVqTAI2dC/CA8MjdYlHmUFoH5rRLqldFgorwVobhhCswJAy9tm4VY5p49TJ48paO1mCPXaRRTYv9htuLt7tmZpCTDItlhkWy32B0CI2ly/C7WA9N0ghAN2FH4W1G27XZHDtiRyucSqnTO/6h2WQA2QctwmRvtoNoejhAiHvjMsasD3pdS/kcIEQp8BMQCB1DMoqt7K8vdw3FjMLj6TnjuTZifDK8+BlOd8GCfcRCufw/WZSprXJ48B5ZO6991Dc3Kavl0R0qtUnKbVN5n04MUgVnoGPYa5zfyoZptUpJvg71WOzlWSYMEC2B15JbWfdm6bwXMErKskj1WSYs9QoCAmXrBTL2GWXrBLL1gik7gMcQ36ZIejkeKJMzJd06Ju4cz4pBS7geSujleBSx19fXS09NJSupyuRGF2jipjQ8MPacX3lLE5vpL4NH/63stTU0jXPFSJR+nhxHoA8+eB5cdBbpevielsu4kvQbSq9sEJq+h7ZwQD5gZAndMVwTmsHAIcsJfWQtcXU/Vdslea9eUY5V0F2mgZRmkXigvHb0AYbHg7aF37AvGawXLvTTM1Atm6TXEaUEzzAraYqjgEqj27TpyUG0Px5VwpodjtVrR6dTVQtTGSW18YGg5/fwrHH0mHLsYvnijb7H5ejtc+gaU1UmuWCq49++9r4/5rAAez1TEpWVYTAATAyApGJJCHHmwMp8ymHfvYOqp2i75tdnOpmbJb812Mq2SSnvb53pgvE4wqX3SCibqBAEa5fPuhEOt7Umv1w++h+OZIhntZA8nz93D+cshJyeHyZMnjzSNDlAbJ7XxgaHjlFcAKy+BhHHw3rO9i01tI9zwPry5ERLHwDMr8jl1cVyP51eb4do/4L08mBwI58W3Ccz0IGXtiqvhbD1JKSmwwaZmOxsdIrPbqvwo1QPJesHfvTQdxCVOK9ANQA3V2p5cAvccTrdwC44DMTEucG7lYqiNk9r4wNBwamiE5RcqFmBfvgmBvXjO+WEn/PN1JcjYXSfDv5eDxRze4/nfFMIlv0GFCe5NgjsSQT8Mq+F6qqcmu2SfTfJrs1QExmyn0NF7CRCwwENwlreWIzw0zPEQeLtwiGso21MzdnJoZA/1jMKLhYQMLycN7vhr3cAtOA5UVlbi5+c30jQ6QG2c1MYHXM/JbofzroXde+G7d2FCfPfnGYxw8wfwys8wdTR8eq3iWRmgpKgrp7pmuGErvJELiUHwzRKYFeoy2j1CSkmZHTZX1dGk8Wa/TZJrla15abuhsWgNHOGh4XBPDUd4CKbrBNohnENx1bNrxk6uQ1yyaCCLBnJoxOZYwn88EU4LTmVlZd8nOQv327UL3FXigNpepKA+TmrjA67ndN/j8Nl38Pg9cMzi7s9Zswsufk0JB3Db3+Cev3cMONaZ00/FcPGvUGRUXOb/3wzwdMFwS7NDTEptklK7pMQGpXZJqU1SZIf9DmFpkoA+DGqtir8yLcRrBSd4aYjXChJ0gnl6DWO1IIZxkr6/z86KpIpmKjCTQyNZNLCHenJoxOoQlwB0TMaPfxDDZPyYgj+j8BwyTj3CPaTWLdyC40CLaxs1QW2c1MYHXMvpk2/g3sfh/NMVq7TOaDDBravghXWKG//N/4bDEnrm1GCBW1PhhX3KXM1vi2HuADynNEvJlyY7X5nsFNskJXZFWKrs3Z8fqoEojTK3ssxTEZVgQw1zIkIYqxV4jbTttAMt9SSR1GGlHDOVDkGpoJlKmqnETLlju5rmDm7H/B3icnaruPgRhRfCaSdmPXMaNNzRCbqFW3AcsNt7+O8dQaiNk9r4gOs4pe9WhtIOS4YXH+pqEfbzHrjwFThQBTceBw+sBO8ezJLtdju/lMGFmxXT5pumwv0zwbuf/237rHZebbTzltFGuR3CNTBeK5igEyzSCEZpBaM0gigtjHLsR2rodp1KUZ2Z0brhd51oQ1JDM+U0U465NVXQzMEQA3UcpBxzJylREIyeMDwIw4PJ+LVuh+NJHD6MHqS4dAeXtXF3/LVu4a4SB3x8VOLHox3UxkltfMA1nCqqFCOB4CD49FVoH9be1Ax3fAxP/gAJkfDLnXB4L46UmqzwwP5QXsqDeH/45Vg4vFv3st3DJCWfGO280mRjQ7NEC5zkpeFSHw3HeGoGPKfiinoyYiOHRmqx0ICVBmw0Ym3dbnBsNzq267FSjaV1LqUFOgQReBKi0TAVXxYTRjieRDjEJAwPQvHAYwR8C7usjQvox0jeXwZuwXGgurqa4OBefYAOO9TGSW18YPCcmpsV8+eySvjlU4hqJw67CuHsFxSvAVctUwKc+fbwErFLeHc/3Lkdipq8uHoSPDQbfJ0cVtllsfNKk513mmzUSIjXwn/9tVzgoyVKO/hf8f2tJzuSAozswsBu6tmFgVwaW70AtIcegR86fNHihw4/dMTijS86wvAgAk9HUraD0KNBkFuQy/jx4wd9b65EdXWvTkuch3tIrVu4BceB6OjokabQBWrjpDY+MHhO1/0f/PK7stZmzkzlmJTw7E9wy4cQ6A3f3AQn9LJI/+dSuGkbpFXDnFB4a56ZpWP6/nlrsEtWm+y80mjjd4tED6zw0nCJr5ajPIRLV9n3VU91WFqFZTf17KaeehS30b5omYY/5xPLVPwJwwNfdPg5BGagPZE/Y3tqhXtIrVsMqEqEEJ6AlFJ258XikEReXh5Tp04daRodoDZOauMDg+P04ttKuvVKOHuFcqysTpmr+W6nIjKv/xMie4jcuc+gGAV8cRDG+MC7h8NZcZC1JxfonlOtXfKVyc7HRjs/mO00A5N0gkcDtJznrSXcBb2Z7tBST3YkJZjIobE1ZdPIQYyAsnwkHl+WEsZ0AphOAGPxRuPiuZL2nNSEvLw81xTkFpxu4ZRrGyHETOB0lEib04AW28EGYBdK1M3VUsrtQ0FysHDGtY3dbkfj6sAkg4TaOKmNDwyc0y+/w9LT4ZgjlcWdWi18s0MRm3oTPHoWXLm0e3cylSa4bye8sFcxBLhjOlw/pc0ooDOnarvkC5Od1UY7P5ntWIAYDaz01nKat4b5ejEk5sh1WMilkVwayZYN5Iom9tNEk2NgTACj8SIBX6bizzQCmIIfvsP0pnRVe7JgxUADdTRQR6MjV5KBBiYylmM4zGlOWq128K5tQlIkS510bbPa7doGACHEicC/gRSU9pkP/AFUOfZDgATgDuB2IcQ24D4p5TdDyHlIsGPHDmbPnj3SNDpAbZzUxgcGxslshotuhLhYeP85aLbBLe/Bc2tgxhhYfwVM62bBudkGz2TBAzuh3gqXToB7kiCy04ryHTt2EJs0i89NdlabbKw1S6zAWC1c66uIzBz94IbMbEiqaaas1fKr/baZEkxUtHOj6WuDSbpATiSSBHwZjy/x+OIzTItF7Ngx0YwRM0ZMGDGzO28PUeNHY8KMETM27NgcgQWUXDr27O1yZbsRY6uoNGHqcj0tGgLxIwA/tP0Y8tuxY4drbtg9h9MteuzhCCHWAUei9F7eAb6TUpb1cG4k8DfgH47vrJdSLhsKwgOBOzyBG+3x0DNwx4Pw4wcQEacYBmQWwQ3Hwn9P67iIE5Q5ndUH4LY0xcz5+NHwSDJMC+p4nlUqczKvNdpY36y42Y/XwmneWlZ6aUjuR09GIqnBwkGMrakQUwez4s7WX55oWifnI/EkHl8SHCkMD5ebEHeGETNFlFNIOYWUUUwljTTRhBlztz6kO0IgHIHTNK25QKB15JrWP4Ev3gTgSyB+rcKibCvHfAYxDOiS8ARhKZKTnXznvOHu4QDUArOklDv7KsQhRK8DrzuG3+52CbthxF81uFh/oDY+0H9OxaXwwFNw8rGwywi33wMhvvD9zXDsjI7nWu3wUT48sht21CguaX5cBkd3mldusEtea7LxRKONAzaIsZq5LciHld4aZup6FxkDFgraiUpLKsDYwSZMiyAKTyLxZDaBrZZfke3yAHQ9ioqrn52BBoewKOJSSDmV1LZ+7o8Po4lgNOF44Yl3Nyk/M49ZU5PwwQsvPNCqYGm+OwDb0MIdnsCNvxTOuxY+/BLOvwVe2QInz4JXL4bwdg46GyzwWg48kQkHGhUvAbdOU7w6a9uNzpTYJM802nih0UathMM9BLf4aTnRU9PrcJkVyUaq+Igi0qhrPa4BRuHFGLy7pGg80Q3zuhQrNqoxUEkNFdRSSQ3l1FBEOQba4liHEkgMEcQQ6UgRBKI+N0jOwiU9nIgUyUon3zkvuHs4fzmkpaWpbn5CbZzUxgf6x+n3VHhnNdx8BbySAafNhQ+vajMMKDMqczTP74WaZjg8Ap6eCyfGgKadfuyx2Hm00ca7TYoBwAovDTf7aTnMQ9MrpzosfEkpqymmFDNReHIpY5mAH2PwZjReQ7bYsTtOduw0YqKeRqqoo4IaKqltzasxINsN23miJ4wgJjG2VVhGE4EPXp0vN2BOI420tDTXFOS2UusWA+7hCCGCgXmAP5Ampcx1JTFXwm2l5hqojQ84z8luh8NOhMISuO9RuOQtWH8HLJ6imDc/thveyoVmO5wSC7dMg/ntogxIKfmlWfJog42vzXa8BVzoreEGPx0JOtHpWh055dDIRxTxPeWYsZNMIGcwmsMJRevieRUbNqqoo5Ja6mmiESMNGGmQjTQKU9s+TRgxdXEo440n4QQTTjBhBBFGEOEEEUYw/vi4dB5Ire3JJVZqkSmSc5zs4Tzh7uG0QghxHXAB0Ay8KKV8QwhxPvAM4NvuvJellFcMFdGhRlZWlurWBKiNk9r4gPOc3v4Ytu6At5+GN/6A+AjwCIEVP8PnBeChgfPHw03TlIib7bHZbOcGg5WtFkmYBu7113Klj5awHtbMZGVlMWnqFDZSxYeOYTNPNBxHBKczmoS2f5sBQSKpp4lyqimnxpErqZI67HT0B6ZFg6dVT5DeH198GE04vnjjhzd++OCHD6EEEkYQvsMYxGUg7cmKGQPFWDC22rD1lfyIJJxe/BF14uQSuOdwukVfZtFnAE8AFsAEvCqE0AIvApnAWhTjv+OBS4UQ26WUL7uSoON624AiKeWJQog4YBUQCqQC57piAWpcXM8RGkcKauOkNj7gHCdDPdz+X8Ux5/zD4bwvYeF8WPgDBHsoIQOumdzVvBngI6ONc2usRGnhhUAd5/toeg1C1oiVLRO8uJM/KMHMKDy5mjhOZhSBA7STtSPJIo80sihzCIsRc+vnOrSEE0wUYSQxkQhCCCcIf3zxwwcvPDBZTXjr1RURrLdnZ6MZA8XUUkAdB6njILUU0kAZdOPoszfEc5TTguOyNu4OwNYt+urhXAHsBRYCNSiWaE8BG4BjpJQ2ACGEF8r6nIsBlwoOcB2wB2j53fkw8ISUcpUQ4kXHNV8Y7EWKi4tV59dJbZzUxgec4/TAk1BWAV+9BW9tUuZjNmvh4gR4cg749aADTzVYucFgY4GH4MsQPSGa3oeT8mjiNnZzQG9kNoFcz3gOJxTdAIehmrGwlUx+ZhtlVOOLN6MJJ5kpRBBCBCFEEkwwAWj6mPtR67MbOz6Gekqoo9CRChzCUop09NQEWgKIIoQ44lhEEGPwwK/ViFo4zKY77mtbtz1w3iFncXGxa27OPYfTLfqqkknAU1LKagAhxJPA+cCbLWIDIKU0CSHeA/7lSnJCiBiU9T3/AW4Uin3pEuBsxylvAffgAsEJCXEuIuBwQm2c1MYH+ua0LxeefBUuPANmz4C/vwGeETAmAp6Z233IALuU3Gaw8WijjRVeGt4N1vUZWnkdFdzPPjzR8HB9HIv9xwz4ngw0sJEdbGIHjRiJIZJzOYFZTEY3wHGakX52JgwYKMJAIXUUY6CQmnEH+Z0qWnosAg1+jCKIMYxlAYGMIYgx+BOFdphWUbq0npwXnDDHovkWvOzqkSK1oK8qiQCK2u2XOPLufgYUQT9+SjiHJ4FbUQwTQBlGq5VSWh37hcBoV1yoqalJdZ6Q1cZJbXygb0433gtenvDfO+CnXVBUAyTC6wu6FxuzlFxYa+UDo52rfDQ8FajrNSSAFcmL5PEOhUzHnweZisVQ2dZi+4EiylnPNlLJwo6N6SSwmBQSiBn0ZP1wPrsmqigjkwr2UEsBBooxY2j9XIsHAUTjb44lwWcJgcQQQDQBRKOlhyBDw4SmpibXFNS/OZzKv4rRQF8mIgI6eCRv2e5uENWlC3ocbnXKpZQDWoklhLhUCLFNCLGtpKSEyspKSkpKKCoqoqamhtzcXIxGI5mZmdjtdvLz84G2hV9paWnY7XYyMzMxGo3k5uZSU1NDUVERLeXl5+fT0NBAVlYWVquV9PT0DmW05BkZGZjNZrKzszEYDBQUFFBeXk55eTkFBQUYDAays7Mxm81kZGS0flej0bSWkZ6ejtVqJSsri4aGBvLz8/u8pxYTT1fdU3s+A72n9rkr7ik/P7/He3pvdTXfrIGrL6jGz6eB//usEfRwxuQKFkZ05fPrrt0cV9HMB0Y7d3tYucVQQlVFRY/3VEMzFxk28w6FHFGu4wWSKEnPQkrp9D1Z7Ta+yPmRZ/mQh3mLHexjQnkUd9ovYmHmJGKMYezP3T/otldWVjZkz6mkdj9bylaz2fosq82X8RmX8StPkc8mGuobGcNcIgqXcaT9DqZm38DJxleYlHsF44pOI7RoIR4lcdgq/TiYXzyk/0/O3FNdXdu6qEGhxbWNM+kvhF7NooUQduAcKeUHjv1QoAJYKqVc3+ncc4C3pZQusc0QQjwInAtYAS+UOZzPgGOBUVJKqxBiPnCPlPLY3spyxiy6srKSsLABxP8dQqiNk9r4QM+cmpshcYnyK2jXOihpgHE3QGA8FN3ZNU5NsU1yfJWFTKvktSAd5/n03owzqed2MqmhmduYwImM6pMTKAYAFdRwkFIOUsZu9lNONYH4sYjZLCRpwOtaeoMrn53Sg9lNGbspZzf1lAKgx4dIphLBNCKZRhBj0fTyM1+t7Sk8PHzwZtFjUyS3OWkWfZXbLLo9HhJC3OHY1qL8D78qhGjsdF4PTtwHBinlHShOQRFCLAZullKeI4T4GFiJYql2PvCFK67X0NCgusavNk5q4wM9c3rmddi3H75+Gzw84LxPAAlPLO8qNnssdo6rtlBth29C9Bzj1XvH/wtKeIQcwvDgZWYypdP4WQsnRVyqKaCsVWAKHQGVAfToiGUUxzGfWUwaUtcuDQ0NhIQFY8aAkRqaqMZEDU3UYKIOG83YsXZJtg77FiyYaKLSwV8RmAkc65TAdMdpqNqTRCIxYqcegQ4toU5zchncRgNd0FeVFKAIjH+nYxq6jlLbHZ8NNW4DVgkhHgC2A6+5olC1vUhBfZzUxge651RWAfc9Accvgb8tg01l8EsaRETAhTM7nrvJbOfkagseAjaE6pnt0bPYmLHzGDl8QSlzCeIBpnQxdS6mgi3Re/mC3yikDDMWQBGX0YQzl2mMIZIxRDKKUJeLjAUjFWRRRQ5NVGOkGiM1NMZW8yuGVsuv9vDEHy2eaNChQYfWkbckHZ7tjnkQQtyABKYz+tuerFTSzB6aycJKFZJ6bNRj7zY1gCOAXABnEcH/hoRTj3BbqXWLXqtESjlumHj0Cinlzyheq5FS7gfmuvoahYWFTJ482dXFDgpq46Q2PtA9pzsfhCYjPHEPGK1wzjdAI9x1asfvfma0cXaNlVit4PtQPXG6nifmyzBxO3vIpJ7zGcNljOvgJcCClR/4jTX8gdDBGEYxj+mMYRRjiCSS0H65yXcWzTRSzh7KyaScTKrZ7xAVgReBeBOENyEIQzBRQXF4E9wheRE4bBZgndFTe7JjwkIOZvbQzB7MZNHMHmyOnpUCPRr80eCHlgAEfugZ7TimHNcQgAY/PHC+zRYWFrrgznALTg9wV4kDCQkJI02hC9TGSW18oCunbenwxodw46UwKQFu2QYF2eCph/MXtJ33TpON82utzNMLvgrR9+g1oBQTm6jmFQ7QjJ2HmcpiOv4KzuEgq/iBcmqYyzROsh9BoGYAZmpOwEy9Q2B2U0YmNeQDEg06QpnANP5OBNMIZyK6dnNBVj8ruhH+d5dYsVOHjTrs1BI9oZp69jj2a2gmh2ayaCYXWoPEeeLBRHxYgidT8GAKnkxGQ8iQhFtwWRt3C063cFeJA7t37yYpqZfA9SMAtXFSGx/oyElKuPZfEB4K/74Bfq+Ax3aBvgLOmAuBDqP9EpvkqjorR3gIvgvR49NuQacVSQYGfqWazVSRi2ImOxFf7mcK49pZ/hsx8yUb2Ew6oQRyJacxmXGk7053up5sNGOmgWYaaaaeZhod+y3H2rYbqaSOg4BiWhzGRBI5jUimEsoEdHg6VU+ugMSCjep2qap12966X4ONWuzUYqMOSaf5kU6jcTrG4MlkfDnWISxT0DMOMYyvqd27d7umILdrm27R45N0BGDrL6SUcukg+IwY1PYiBfVxUhsf6Mjp/c/gt1R47THw9IUL10NwDVQ3w0VHtn3nFoMVs4RXg3T4aAS1WPiNajZTzRZqMGBFi2AmAVxLPAsJYSzeHX5Rp7OPj1lDPU0sIYXjWYinYw1Jb/VkwkA5uyllF2XsxkBvQzgCD3wdyQ9fwhnH4UQwlVAS+jUUNtBnZ8eEmd2Y2YGJHZjZhY1y7O3W1XSGhiC0hKAlBB1RaJmChkC0BDk+C2yXtxwPQKjARthlbVzgdm3TDXr76RBP17U1vtA6nlDryIMceSV0/glz6ODPEFxsqKE2PtDGSUp46FlImgoXnAF37YCsOkishZBIWORwpbXBbOc9o50b/WCzrpD/Uc0uDNiBYPQcQSgLCWEewfh18+9RRwOrWUM62YwmnEtZQWw7k+j2nADMNFBOpsOMeBe1HABAhxfhTGYcC/EiEA/8WoWlbdsH4aJ5H2eencSOhRxMreKyAzN7aJl81zIKL2ag44hWQdESgobQdvvBTvdIFE7qcrfj0gBs7vGjLnA6PIEQIh5YD3wKPCylLHUcHwXcDpwCHCWlzBsaqgOHOwDbnx/bM2D2sfDCQ5B8Ahz2HawIg9UfKGGj7zgJLFIyq8JCvZQkh2dQpGliCn4sJISFhDIZvx7DEtuR/M5OvmADVmwcxwKWkNLFysyGhVJ2UubowVSTB0i0eBDOJCKZTiTTCWU8mmF+I0mkY3irEiuV2KhwbJdgJgMT6a3DXhr88SQJL5LwZCZeJKEjalj5jhRcEoBtYorkWSffOce61+F0hyeAX6WUN7Q/6BCe6x3C8wSK8BxyUPOvd7VAbXygjdPbq5X1NstPgKN/hVFeEFunOOo8b6Fy7tONNnZbJaeHlFKgaeJZEplD3+5eKqnlfb4nh4MkMIYzOZaIbr5Xzh628CIGitCgazfHMp0wJgy5NZjETDM5mMnCQg5WyrE5hMXYXIzGwwAOM+2O0OPJFAJY4RCXmegZ77LeVU8YUHuSVrCWgb0e7GaQJpBmJdkdecuxln3PKeDf69rwDpxcAncPp1v0p0oWo6yB6Qk/Aw8NhsxIQm0vUlAfJ7XxAYWTxaLM35x0NDx/EHbXwhdHwhWPw/EzYHQIFNkk99TbSPZsJs/zAJcx1imxOUAJL/EpNmycxbEcRmIX66hmGtnBe2TzIz6EcQQ3E82sXifxBwOJHSsFDnPhLJrZi5m9WNhPm/cpHVrC0BGOljACPSajJRwdYWgJQ0tE6+cagoZcXLpDh/YkJdgNYCkCa5GSt9+2FjvyMuhmLVGvCLrAacFxWRt3C0636E+VSGBKL59PGySXEUVGRgaJiYkjTaMD1MZJbXxA4VRQlkh5JSw8Fm7ZBefFg64WimvgmX8o591UZ6VZSkTgHuaJIC4gts+yM9nP63yJHz5cyVlE0NWT8EH+YCuvYqKWyfyNGZxJVkYOsYmuExsbNTSxASO/Otam7EPS5mRSRyyeTMaP4/FgMp5MRk9ch7kUVT07KcFSyIG9qxkbXg7GrWBMBXtt13O1waAbDfrR4DlDyfXRoAkE4QkaTyUXXp32Hcc0niCc9ync4nvNJXALThf0p0p+BK4QQqQC70jH5I8jZMB5wGXA5y5nOEyYOHHiSFPoArVxUhsfUDg98CyEhcBn3hBmhyfmwCUvQbg/nDgL1prtfGiyM8W/jEidnXuZ3Gdo5y3s4gO+J5pwLudUAvDr8LmRGrbyGgf5nSBiWcQthDGhldNgIJE0s4dG1tHEOkykAnY0BOLJNAI4s1VYPJiIxokIoiP67KyVDlFpl6xljBVApQ68ZkDg6eA5oU1cdNEOYXG1A/re4bJ60sAQuMQ75NEfwbkRmAO8geJfLdtxfAIQCRx0nHNIoqCggAkTJow0jQ5QGye18QHYtbuQL34cz6kr4f0qeHg22Mzw5Xa47hhAK7m62kKQ1kKAXwH3M4OQXlzgSyQ/sYWv2cgkxnIRy/FuNzQmsZPLWtJ4BxsWkjibqZzcwQBgIPVkp5EmNtHEeppYh9URCcSTRIK5Bl+W4skMxAAXd7j02UkJsglsNWCrBmu1kndO1gow7QBLvuOLAjwng9+x4D2HgspRxE44ETTqeTMXFLjIO5d7SK1bOF0lUspCIcRMlHmc5bS5l9kPvAn8T0pZ62J+w4bIyMiRptAFauOkNj4Am1KjMZvBnASewEUJ8PZ6sNrgokXwRIONLCtMDMnlCjGWWb34mLVjZzVr2cQOkpnCORzfIeCZgSK28BLlZBLJNOZyGQFEdymnt3qyY2q3ULISC7k0sh4jvwPNCPzw4QhCuBEfjkKHa+q8W062emVuxFrqEI86sNcpeet2bcfj9jpFTHqL6i70oA0FbQh4z4GQK5XcezZoA1pPC9IbVCU24MI27hacbtGvKpFS1gF3OtKfCrW1tQQEBPR94jBCbZzUxgfg3U8EkxLgewFnjoNQT3jtFzhsPPhFSu6tsBLiVcMJXlrOpeconBasvM3XpJPNEuZwMke2mkjbsZLJF2SwGh0ezOMKxrOkW9cqVsqosr2DBatDVKqwUoXdkXdZbQ/oGU8Q5+PDUryZgxhsEDJph+b9YCl0CEoxVO+FugawlLQds3d2+N4Owge0gUrSBClzKfpxjmMh7VJwx31diPLdPiKkgjrbU21tresKcwtOF7irxAEvL3X90gL1cVIbn9x82LbTi5Mvgb02uHIS/LEfMovg5QvhWkMzzUgWBJRwD4k9rrFpwsQrfEYuhfydozgKZUmERHKQLezgPeopIZb5pHAR3j1Yt5nZQzHnYQsupQadYyFkKFrC0BOLD6GO/bakIxp9N72kfkFKaN4HjeuhYZ2S2yo7nOKPNzAadFFKT0N3ojJHootSkjbYIS6BSi9EDP2q/2FpTzYL2C2gd24uyGWcBKBzaUzKPwV6c22zVEq5diCFCiGWSSnXDJyWG270jXdWgxCSzDhBSijMDYNLvwQfDwiZaeMLI4zxL+IJXUKXMAItqMHAC6ymglou4CRmOzwLV5BFGu9QyV4CiGExdzCank1mm9hICZeiwQ/f6ncZFbJoSJxLtqL5ADSucwjMOqXHAqCLAf8TwGcReMS1ikpFpYkIFQ6J9htN5VCVDk2lYKoBcw001yq5qd22uQbMtWBpgKn/hKNeGV6eGglevQw7/kXRWw/neyHERuBx4Dsppa2XcxFC6IETgeuB+TDCwcn7CZPJNNIUukBtnNTER0p4ezVMm9nMLo0nb0yCRjOs+h1OnSu5utmEl9bCfX5eJNL9sE0xFbzAasw0cwUrmUgsBorZwXscZAteBDGPy4hnSa9xXwysppxb8CCBaN6muMGCCHGx2FhKlJ5Li8hYHA49tOHgtwR8l4DvUeCR0O1wlsnsotDJLkSv7clug7ocqNzRLqVDU0nXc/X+4BkEnsFKCohv2/YMhgjnF/G7ro1LNFqrU2f2c1XRIY3eBGcWith8CVQIIdYAfwC5QDVKpzEExUrtMGApil+1H4GZQ8Z4iBAUFDTSFLpAbZzUxGfzH5BXALOO1xLiAWeMg49+g3oTNM4xUWrTsTKkjPNFXLff381+3uZr9Oi4jrMIxZOtvEo2P6FFzwzOYAondXDx3xkSSQ1PU82jeHM4o3gJLQEEBfXs2NJpWKuh8WdFYBrXgXmPclwTBL6LIewGRWA8pzk1XzLsz85qhIaDYDUpc0rYHdZtdse+JLS5HooLlM/sVqjd1yYuVRlgdaw10uggeBqMORrCZkJYEviNcQhKkPK5i+CqehJCotP3+hu9FX+lflCPT0pKuQs4RggxH7gSxTLtLLo69BSAAcXH2gtSyq1DxHVIUVZWproJTLVxUhOft1eDjw+kR2u5MUFpyI99D2MjJZ9HQrRXLa95jekyrFVEOV+wgSzyiSSUSziJcn5hE59jw0wCR5PIaXi3+qTtHhILFdyFgQ/w51Qi+F/rZP+A6snWAE0b23owpu2ABI0v+BwBQReC31HgNQtE/02jXf7s7FZoKARDniPtV/J6x35TaZ9FdLt6yDNIEZVplyp5aBKETAHt0Hht6IyysjLXFCRAq3NOcP5K6POngZTyN+A3IYQWSAamAuEowlMB7AK2SykP6Z5hbGzfK8+HG2rjpBY+RiN89BUkzIcMT7hiEjz9E2QchFEX1YPQ8U6gDwHt5m1qqecbNvEHu/DGi+UcyWjq2MRdGKlmDHOZyTkEMLrP69tppJTLaeJngrmGEG7pIGx91pOUissW4w4w/qGITNMWwArCA3wWQMS9yjCZ9xzQDH50ekDPrrke6nLBkKsMb9XlKqk+D+oLoP0ou9AqvY6AOBh7gpL7jwWdr6MHpgHhSAgQGpqtVjw8vFr3CRyvlOFEj22o4Ko2LoR0C0436M86HBvKkNofQ0dn5LBv3z71uP5wQG2c1MLnq5+gzgAHJ8HhAfXomwO4+1MYP8NM7gRPzvZrZIlWcUNjopm1/ME6tmJHspgUljGHbTxNKmmEMZHDuZEIJ8MQWymjhAsxk0k4DxHIOV3O6VBP0grmfUqPxbRDERnTjnZWZBpFVMJvVQTGZwFoXB9IpdtnJyWYqqAuu6Og1OUoImOs6Hi+dzgEjIfI+TDhbEVUWpJvDGj7Z9m2NyODxPEj357aY9++fS4pRxEc5+Zw/kpwm0U7oIYXaWeojZNa+Ly9GkIjoCoGbksJ4Pr3wC4lB/9mJkQPr/sFYcPOb+zkOzZTTxOzmcyJHEEoAfzKMxSTRgoXM5HjnLYmayZbMXummihex5clXU+yVpAYvQmKnnWITIbivRgU/15eiRBwCnjNBO9Zin8wrV/XclwJcx2JkWbY94EyT1KX3Zaba9udKJQeRmACxJ2i9DgCxiv7gfHg4drhVLW0p/ZwFSchJJ4ef6XZGeegWsERQngBv6AsINcBq6WUdwsh4oBVQCiQCpwrZW/Lnp2Dml3vqwVq4FNWAd+vh+jjICAA9v2ezafbJjD55EayAj15PNBKtsjjCzZQRhXxjOYS/s44opFI0niTfDaSxNlM4ninr2vkd0r4JwIPRvMxXszoeIKUYPgIiq8CW5WyrsVrlmOV/SxFYDwnDf36FlMVHPgeitZB7V6ozQZjebsTBPjHQuAEmHAWBE1UtgMTIGDcsM2VgAvbk80K5iZoblLylm1To5IHR0PcLKc5uQICiRanh9TChBDtg+e8LKV82SVEVAbVCg5gBpZIKRscJtebhBDfofhre0JKuUoI8SJwMfDCYC820i/S7qA2Tmrg8/5nYLMpw2kPxMEz704gPsrO3nkaJnoZqPNcy8scJIJg/skpJJLQ2oPJ5HOy+IZJ/I1p/N2p69lppI53qOIR9IwhmrfRd/Y0bS1ThMbwiTI8Fr0GvJKGZy5CSsVc+MC3cOAbKPtdsQLzCoWQaRB3siIorcIyHnTqWMDrdHtqNsGBdMjdBvu3wf5UqK9sExdbdzF+2mHxhXDl667l5AT6ITiV7gBsIwyHN+oWPyB6R5LAEuBsx/G3gHtwgeCo4dd7Z6iNkxr4vP0xRCRA3Sio2gt5FTD+qkaETsvSgA0UU8FKlrKQpA7ROHNYyw7eYxyHk8z5fQ6j2aimljeo4w3s1OHDYiJ5Gm17LwPtezX2eoh8CMJuIjUtneTkoVz02QCFaxWBOfAtNBYpxyNSIOXfMPZvEJHsmKBXkJqaSnK8uiKIdNuerBYo3A05WxVxyd0KBRlKLwYgIBzGp8CEw8DTBzx8lLzztqdv27GgUV0v3gsnV0Dp4Rz6czhCCJuUcmAeY7vB8Edd6geEEFohxA6gHPgJZQ1QrZSy5UkWQvdmRUKIS4UQ24QQ20pKSqisrKSkpISioiJqamrIzc3FaDSSmZmJ3W5HOH6NtjS4tLQ07HY7mZmZGI1GcnNzqampoaioiJby8vPzaWhoICsrC6vVSnp6eocyWvKMjAzMZjPZ2dkYDAYKCgooLy+nvLycgoICDAYD2dnZmM3m1ngcLf+MLWWkp6djtVrJysqioaGB/Pz8Pu8pLS3NpffUns9A76l93t972rHLzo7dUDsd5lPDsz/A7Ol15I7xYLZHITpdHotLk0isiaO0qLT1nlLLvmaLfJGAhgmkWC9nZ3r3fDIyMmg055Fdcx158jBqeBKNeRY+1W9iLfgPjQZt6z3tyfgZDq6Eg2eCx3h2m9+F8NtI37mbpKQklz6n/Xt3YNj3PbW/3Iv54yORr4XCd6cgs1dh8JmG7chX2D13DZy2lVTtSRA5h9S07R2eU0BAwLA9J7vdTlrqNmg2sWPjWijbz57vVmHP/IWCL16i+ed3KV/1IFP2f4/h9ZtpePEKTE+di+nm2cjzA+DWWfDypVg3fQB+oZTM+Qfc9Ak7r/gaXikj48SHMV/wDNkLr8Jw4u0UzDmP8rn/oHzmCgriFmOYegzZfgmYxx9GhtELwsY4fU/R0YN0M+SAQKLD5lRSOVz6y0k4wtqoGkKIIOAz4N/Am1LKBMfxMSheEKb39v2UlBS5bVvv8cXT09NJSkpyDWEXQW2cRprPLffB46+A/Q5IroTcYgk31mMKlFwb8RWBopnbubCDz7QyMlnPAwQRy1LuRk/3FmDNZFPDC9TzGQD+nEIwV+BBp/goUkLdh1ByNdgbIOI+CLsRRNtgwYDryWaGmiyo2gXVGUpelQEN7VzmB01SejDj/gZRh4PWOZPpfnGy28FQDlWFUF2kpJbtmmJoNoLVDBZHsjZ33e9rmKs99F7gEwDRk5Xey/g5Sh45fthNpNPT05k5c2bqYIe4vFOmyrht7/V5ng0t+0TSoK83VGjp4Qgh7pdS/nuw5blsSE0I4QfcIKW831VltkBKWSuEWI/iMidICKFz9HJigCJXXGPaNHUNN4D6OI0kH6sV3v0UfKZDuITUHDjhTCvf+nny98A8jKKE0zmxg9jUkM8GHsKXMBZzZ7diY2I7NTxPIz8g8CSQ8wjiEvTEdEOiDIqvAMNn4D0PRr8BXl2D4PZaTy2myIb9SqrNhmqHsNTua1vbotFD8GSIWgghl0FoorLC3n9g60S65VRdBFmbIOcPqCxoE5bakrYhrBZotMrke0i0MlzlE6CInd5TSbqW3KPjvpcfePuDl7+Se/uDdwB4+2PVeaPzDwatekb2XdfGnRtSG4lhNyHESuBwYCfwdrsRI4QQ30gp/9bN165D+cE/KAz6STuE5lqUyfxgwCWCI4QIBywOsfEGjgYeBtYDK1Es1c4HvnDF9XJycpg82bm1GMMFtXEaST5rN0FpObAENLtgVrzkpyQrwZ51TPfaigchzGJS6/n1lLGOB9DhxRL+jVcnf2omdlLFfzGyGQ2BBHMtQVyIltCuF5cS6lY5ejWNEPk/R6+m+6Ht3L27mRTtCXX721bgtwiMYb/iULI9AuIhZDrEr4DQ6RCSCEETnO69OIOcffuY7GdXBCZrE+zdDBX5yod6LwiLhdAYmLYYQmIgZLSSQh3bgRGK6LgQOVlZTA4Kd2mZg0VOTo5LylHC4ahvuEwIcTXwLxSXZbcAlwkhjpdSVjtOOaKnr/ZQ3uVSyhedvX6fgiOEOBO4A8VnWjXwDnCXlNIuhLgUeAAIA/JxbZycKOAth4cDDfCRlPJrIUQmsEoI8QCwHXjNFReLienmF+0IQ22cRpLP2x+D3ldZytJQBUGnNmMVcKpvCbWUci4noHFMSRqpZT33Y8fKMdyPLx1fao2so5TL0BBAKP8ikHPQdAoh3XbyJii/W/EG4H0YxLyhRK3sDsWbIO1BJh74jg4eoHTeiqgExMHoxW3bLbm+7xDR/UazUZl437sZsjYxae+v0FSrfBY0CiYfDidcp+Rjk0A39OEIOkNt7Rtcx6mfZtHDiauBY6WU6UIIHfAcsE4IscQhOj2NYXo7rITTHSkDsAH/A1wjOEKIk4D3HbuVwCjgVkAKIYKBy4Acx7F3+vIo3R9IKXeiOBDtfHw/bdFGXYbKykr8/IZ4AV4/oTZOI8XHUA+ffgeWqUARnLpU8kkojPErIVJsx5MgZqMMbVloYj3/oYlqlnI3gZ2Crhn4hHJuwpOpRPE2OsK6v2jjRii/FxrXgi4Sop6CkKu69mqkhIIfIPW/ULIRvMOpG385QXELlcWS/nHgE+nauYiWOZbKg1Bd2DWvOqgMj7UMi42eQsP04/CfczxMWgiR8SPqPqYFA25PTU3Q1KjUg3Q4BbXbOybZbjsgECKjnObkKqhUcKKklOkAjqG0y4QQjwPrhRBL6OorswUWFMfMM1B6RlNQLIff6M/F++rhXIdiIXaMlHKnQ2Q+RQlBoEfp+TzWfgzwUIWaXuwtUBunkeLzyTdgMgH+EB4gST/KjIe2mXP86jBoqjiH49GiwUYzG/gftRzgSG4nvN0QG0Atr1DJfXizgCheRYN/14s1boTye5QejS4SRj2uzKFoOgXwstsg9xNIewgqtysr9I94GqZcjLWuCcJ6ELL+wGaFwkxl3cn+VMU8uGWupfOkvN5TGQoLGwOTj4DwsZAwDyYtAP9QzJWV+LuCkwvhVHsqL4WM7bBrh5J274D92YrIOIuzLoInnRsIcVUbV3EPp1IIESelzGs5IKW8UQjxJMp0RU+aYJVSPtH+gBDCS8oWNxrOoS/BmQU86+htIKWsEUL8C9iIsvjy4f5cTM2wWPphVTNMUBunkeLz1scggkHq4ISVkrd0gqkBB/HTZKKz+pKiU3o3W3mNMnYxn2sYzezW70skVTxMLc/hywlE8hSazmEHGn9x9Gj6EBpbM+x9B9IeVlzDBE2CJW/AxLNb51wslgHEnrFZoWiPIiy5jsWNB9KVoTFQJt/HJik9lNAYCB2jiEuLyPiH9dprUVtbgk6cbDZFSNoLy64dUNHOe3NsHEyfCSvOhhDH/WocDkE13SQcn4+NHxinQUAg8cTskrJcjLXABcDd7Q9KKa8XQjwN9GTxu6vzgf6KDfQtOEEoa1/ao2VWbX1/L6Zm2O3qc3atNk4jwSf/IGz4DUiA6QmSTyY2E+Bp4EwvC2UUc2z9XLTBWqrZTy5rmcLJxHNk6/clVsq5nXo+JIB/EM4DiPbB1Bp/cfRo1juE5gkIubSr0FgaIfNV2PGo4pY/PBmOW634HOs0me5UPTXWQsZayNygLHDM39FRXOJmwdGXQ3yykqImOl6iA8OItSUpobYGSouhrFjJS4uhrISg/Fyoq1GOl5VAy8ter4fJ02HpCYrATJ8JU2dAYNCQ03VVPam4h3M1Pbz3pZTXCiEe7XRYOD6b74qL9yU4ArrY7bXsN7mCgFrg4+NczPPhhNo4jQSfdz9xbETBhDOtZAJzAg6iF3sJxp+5THP4SHsLTwKYzqmt37VjpIyraeRHgrmeEG5s8zBg3AaltzqEZlTPQmO3QcZzsO0+xZw5+kg46jUlGFgPPYpu68luU3ou6T8oKXuLcszTF+Jnw7LLFGEZnwJRE1xuETaoZ2cyQV4O5O6F/FyoN4DJqCRjkxIvomW75bjJqMyzVJaDuZtf+kHBeIRFQPQYWLAYIqJg4hRFXCZMAY+RCRjsqjauVsFx+J3s0feklLKg075LnQM4YxY9Tggxu91+oCOfIISo7XyylDLNFcSGG9XV1QQHB/d94jBCbZyGm4/dDk+/DQTDnCMln/nbGeVbwt/1Oooo4jSWYaiuwxycTxm7SeFiPBxhvWzUUcJFmNhKGPcTxAVKodIGFQ8qvRpdmENoLus+JEBVBqz7J5T/oQjMnHsgakGfvFvrqboIdjgEJmMNNFQrIhWfAqfcATOPVeZZhsFCrM9nJ6XS88jZqwhL+/xgfsc5EyHAy1uJgOflDd6d8qAQZdvbG8IjFTEZFa2kyGhlAt/bm/zcXMaPHz/k994fVFdX932Sk/gzuLZxNZwRnPvpfm3N8z2c79qfZsMEV7m0cCXUxmm4+fz8K1SUArMkphMs+GqsTPYvQ0s+gfhxGIk0RxtZx1MEEM0EjgaUmDXFnEszOUTyDP4sVwpszofCf0DTZgg8C6KfB21Q1wtbTbDtAdj+sBLG+Oj3YcKZzll2lWQTu+E5eH4NHNytHAuOgpSTIelYSFwGAUM8ed/cDPsylaGshnpoMDCmuhK+MSm9kwaDktcblM9rq5UeTGO79UHePjB+IsyaC6edC+MnQcIkiJ8Afv4usXJTW/sG13FqcW3jRkf0JTj3DgsLFSAvL4+pU6eONI0OUBun4eZz/4uAHsYeK8nwkMQH5HOSxpd8DnIqS9Gj49faTzB4F3Mkt6FBh4V8ijgHG5VE8wY+LfM5te9B8ZXKdsy7ENQ1cBoAxb/A+kuUVf+Tz4eFjymel/uC1QJfPgKf3IfWboepR8KRFyi9mDHTh9YMuaEetv0Gv2+ELRshbYvDrK8NrQNUOh34B4BfgJL7B8Co0TDvCEVQWoQlavSg5oycgdraNyicXAG1Dqm5AkKIQCnlACxj+hAcKeVfRnDUtKK/BWrjNJx8Kqpgw89AjKRmiYUojyYSvA0IDhKAL/NJpJlGKkZtIJJpjCYFM7sp5lwkVkazCi9mga1WEZq6D8BnoSI2HuO6XtBcC7/eBpkvK4sxT/5RGUZzBtlb4KVLFLPlw1bC+U9CaN+hqgeMygr4Y1ObwGRsV6y8NBqYMRvOvwKSD4OwiFZRsfv6oQkIAk9PVazBgYG3J5vZjL2pCWmzIa1WJe9lWx8Sgs+4cUPKqTMEEg91WqkNCEKIacDTKHHIGoQQ/kAFijuzDGfLUY8ToxHGjh07mD17dt8nDiPUxmk4+Tz4prJ2zycZGkIksYE5nCQC2McB/s5ReKBnOx9hFvXM4jzM7KSYc9Dgy2g+woMEZU1N4T/AUgQR90P47R2cbLZi/2ew4SowlsHMm2Duvc6t/jfWw6p/wffPKH7Gbv0CUk5mR1oas10hODYbHDwA+/dB7j7I2qUIzb49yueenoqwXHsHHHYEpMxXhru6wY60NGbPdt5N/3Cgp/Yk7XbMJSU05eVh3L+/LXdsm4v65z4x5qKLSHzNuXU4O3bs6FfZPeFPOKT2AvBPKWVrDG4hRALwOrDI2ULcguOAml7sLVAbp+HiIyW88DYQBOZjrUzwqWWs3oYkE398WEgSDVSQxdfEsQhf6h1iE8hoPkIvI6H8Lqh4CDziIH4z+MzreqHGYvjlGtj/KYQmwd++VGLKOIO0b+CVK5SV/cdcCWf9V3FoST/rqWWyPnefsg6lRVzyshWLsPbrQgICYc4COO08RWCSUhTRcQLD2ZbsFgvGggKM+fnY6uuxm83YTCbsZnNbMpnwM5vJWrVK2TcaMRUWKuKSl4e9vWWbEHjFxOAdF0fY0UfjHReHLiAAjU4HWm1rLrRahE7XMddq8R471mnurqynP9mQmpa2JTEtyKOfc/Z9ubapp2dXB91BSikD+z5NfVBDcLHOUBun4eLz7LdgqgAxQ2KbYcPDr4CTCWQXGzmZRXigZyvvIxDos30pnnA2GoIVsTE3QeFCMG6FoIsg6knQdvrVb7cpQ2e/3aGEBJj/ECTdCFonrMVqy+DN6+DXDyFmKty/GSZ2XKLQZz3Z7bBpPbzzMqz5RjEfboGXF4xLgIlT4djlysR9/ASInwjhEQMeDnPls5NSYqmpaet15ObStH9/676xoEC5x74gBBpPTyV5eeEVHY3ftGlEnHgi3vHx+MTH4xMXh9fYsWidFNbBYoRCTB8KeAPYIoT4GahBcdR8BPBqfwrpq4eTinOCE4HiW0f9wXV6gJpe7C1QG6fh4CMl3PsKoJPojrIzKqCeeJ3ESia+eHMEs6gih3w2MokFhEx4zCE2H6Ov2wSFF4LGE8ashsBTu16geBNsvAYqd8Doo2DxS4pnZmeIrX8D3rkZzI1w+n1wym2KO/5O6LGeykth1Zvw7itwYD8EBSsWYFNnKIISPwGiY4Zksn6gz85mNtOQkUFdaip127Zh2L6dppwcrHUd54w9IiLwiY8naMECov/xD3zGj8d73Dh0gYGtoqL18moTGE9PhF7fGvhQLXBVG/+zCY6U8lUhxCcofixDgDTgv1LKmv6U05fRwOLePneEJrgZJTQBKC6vD0mkpaWpbghLbZyGg897u6EqFYgWWObZ8PUt5O8EkcovnMgReKAnlbfwxAd/XsNq8iXeazW6pgJlvsZ7Doz5EPSdvP42FsOvt8K+98AvBo5ZBQmn991jMFQqngC+fAR2rYMpR8ClL8PonieXO9STzQYbflJ6Mz9+pQT2WXAk3HY//G2F0qMZBjjz7OzNzdQ7xMWwbRt1qanUZ2QgHcN6+pAQAmbPVgSlpQcyfrwyxDUAH2Rqa99Aa/TVwULFrm0GDIe4/DCYMgY0h+MIGXA5SlyFCOA34DYp5ebBkBlJzJw5c6QpdIHaOA01H5sdbnkFsIFmgiRwciPTPKyYyMQbT45gFoVspYI9xFGEF6FEeaxCZ9VCwUrQjYbYr0AX0q7QZkh/Erbdr2wn36mk7owCmurafJnlblXylpgx3gFwyYuw9JI+eyAzZ86E4kJY9Qa89yoUFkBoGFx2A5zzT2WYbJjR8uysDQ2YioowFRZiKizEXFSEMT+furQ0RVyalUXouqAgApOTibvxRgJSUghMTsZ73DiX9kjU1r7BdZz+bD2cniCEeExKeZOz5/dbcIQQp6PEwEkAsoDLpZQuCYI2ksjKylLdmgC1cRpqPm/nQunPQKDEvshGSGARK0UQv/Izx7MQDzSk8QremIjGg9GsYt+eCqb6XA+2Soj/taPYHPgeNl2nrKkZdyIsfAKCEpTPLGZFVHK2Kj2Y3G1Qsq/tuxFxkDAXjr1S8QwwPkWJVtkTLBbYuxvS/qDx0/fx37JRmcdYtAzufhSOW+4ydy12iwVbYyO2pqYec6vBgLm4WBGWoiLqcnORFRVdhsIA9KGhBCQlMe666whMSSEgORmf+PghH+5SW/sGhZOr8GcSHMd7v8thYFl/ynFacIQQR6FE3EwGSoBLgdellOryMDlAxMXFjTSFLlAbp6HkY7LB7Z+iWPZPFXgmm5njZaaBXLzw4Ehms4dXaKCGqZiJ4UN0jGJC4KNQux5GvwnejvBJdfth842Q9wUEJsDfvoFxJyif2ayw4W34+B4lbgwo3pfjU+DI82D8HMWnmX8viz2lhAN5sP2PtpSRpvgUA/wiRsE1t8PZF8M45z0V94TarVspeO45yr/6CqvBgLQ66TJFCDxHjcIrJobAqVPxjY3FMyYGr5gYvEaPbs213t249RkGDEV7knY7xupqGsrKaCwro7G8nICYGGIPP3xYOQknQ0z3WY4QY4C3gUiUOfKXpZRPdTpHAE8BJ6D4uLxgCFyMvQg8SdcAbf2KHuhMxM8ZKEJzDGBAGUZ7Ukpp7M+F1I7i4mLV+XVSG6eh5PP8Xij/BTR6iX26JHxiEaeLYNazjmOZj4Xf2M2PBGJnOu+gYxTUfYK+9nEIuRyCzwdLkxKfZvv/lPU2hz0IM28AraciEn98BqvugqIspfdy/hMweaESAbM3GI3w688dBabKEajLywsSZ8O5lyluYGbPZb9VMj4hYVD1YTOZKP3oIw48+yx1W7ei9fVl1MqVeEZHo/X1Revj02uu8/PDIyICjV6xvMtVod+yg3l5+DU1UZKWRml6OjazGY1ej9bDA61er2w79lu2W3KzwUBDWRlN5eUdxKWxogJp69izSDrvPKcFp7i42CX35sJ1OFbgJillmmOxZaoQ4icpZWa7c45Hicg8AZiHsmamm3UAg0IW8KKUsrz9QSFEYn8K6css+h3gLBTvok8A/+mvVcKhgpCQkL5PGmaojdNQ8alrhvt/B80usEcKNHMsLPI1Uk0+nngwDxup3I2VcObxb/SMAtMeKLoAqz4Z3agnwVQNXyxTgqFNOAsWPAJ+jsWXGWvh/TuUIbTRU+DmT2HOKX0bDEgJ330B/3e9sgBTCJg0DY45GWbPVQRm8nTFnX47hNQM/F/EeOAABS++yMFXX8VSWYnvpElMfeYZos87D31AwIDLHem2ZDEaKdu5k5LUVErS0ihJTaV81y7sjt6aZ0AAHn5+2CwWZcjQYsHW3Iy9l/g0eh8ffCMj8Y2IIGjcOKLnzsXPsd9y3C8yEv9++EdzZT25YkhNSlmCMqKElLJeCLEHGA20F5zlwNtSSgn8LoQIEkJEOb7rKiyFrlYQUsqV/Smkrx7OOSjduH3AROCtPsZ1pZRyeX8IqAVNTU2q8swM6uM0VHwe2Q21W1CC2I6GkPnlnKMJ4Hv2sYQYyriSUiYzjgWEkwy2eji4AoQ3FT4vENXcBF8eA9W74YQvIe4kpeCcrfDBnYqn5rBYuPINWHSuc67/83PhrmthzbeKqLz7Ncxf1ONK/vbobz1Ju52qtWs58OyzlH/9NQCRJ59M7NVXE7pkiUvmUobq2VnNZsx1dZjq6rrkxupqyjMyKElLoyIzs7Xn4R0aSnRyMomXXcbExYuJmj2boLi4bu9TSom02bA1N3cQI8+AADx8+zWa4xSamlwTdUVxbdNjFIDOCBNCbGu3/7KU8uUuZQoxDiUo5pZOH40GDrbbL3Qcc5nguGpEy5k5HIESx3qGE+cesutwNEPsqHAgUBunoeBT0gSPZ0LgDmgIldhmSI6JraOUXDzRM5HHOch4NHgwiwuUXkfRhWDeB+PWIKoD4KtjoWonHP+5MldTlKW4nNnyiRIJ8/wnlGBmHk6YIBuN8MxD8OzDoPeA+56Ai67q0ovpDX3Vk5SS5ooKTAUF1GzeTMHzz9O4bx8e4eGMv/12xlx2Gd6xsU5fzxWcekJLz6R0+3ZKtm+nPCMDY1VVq7BYOzkJ7QzfyEiik5OZtHw50cnJRM2eTcCYMQghKCkpISoqqtfvCyEQOh0anY6hD+LgujbezzmcSillry4uHEtQPgGul1IaBsuvm/LX9XGKBIxAAfAj8IWjR9Uv9LUOZ8TeeD1NlgkhQoAPgXFAPnC6K4b59P14oQwX1MZpKPjcvxOa88B4EJgi8Jlby4U6P75iL/Opx4yRSjyYznJ8CIXKR8DwCYx6BDySCf11KdSkw3GfQEASvPhPZYGmpw+cdg/87YZWlzN94oev4F/XQUGeEsb4nkeV2C39hNZqpXHfPsW9S0EBps75wYPY272oA+fNY8Y77zDqtNNcsqLeWF1N0R9/oPPywjMwEK/AQKwWC7awMLS9PENTbS0l27dT6kgl27dTuWcP0uE1wCsoiMikJEbNnIlnUBBegYGt5XfIHZ95BQXhGRjYYw9Nbe0bXMfJlb7UhBB6FLF5T0r5aTenFAFj2u3HOI71B/GANxDu2K915EGOvALQoBgmXAZsFkIcL6Vs5yajb6jZl1q3k2Uo8bjXSikfEkLcDtwO3DbYizU0NBAWNsRxSvoJtXFyNZ8cA7ySDROyYJ+HxB4Nxy2oopA96NEwlQ8oYAFe6JnKcmhYB2W3Q8BK8L8UvjoeXfUOOO5j8JwEd8yF+ko4/lpYcScEhPfJAYD8/YrQ/PQ1TJoKn66HhYv7dS8Ne/ZQ9tlnlH76KYbO7lGEwDMqCu/YWAJmzSJy+XK8xozBOzYW30mT8JsypV/X6gwpJZV79rDv66/Z9/XXHNy8uVUkOkPn7d1FILQeHlRkZlLbzjW/f3Q0o2bNYvLf/07UrFlEzZ5N4NixLjWVVlv7BoWTq+CKORyHBdprwB4p5eM9nPYlcLUQYhWKsUDdAOZvFgPrgUeAR6WUFY7rhwO3ACuBo4B64A7gJuD/6Oe7d1CCI4SIRhkrzJZS1g6mrM7oZbJsOUrlALwF/IwLBEdtDR/Ux8nVfP61A/TNkP+rIja6Gc1cH+zJarKYwz7qiaeORhZwHfrmSjh4JnhOhIin4Zu/QdkWTEe+hbfnZLh3sWLy/NA2iHXScMZkguf+B08/qMzr3P0IXHKdU8NnUkoMqamUfvopZZ99RqNj/UbQYYcRe+edBE2ZgndsLF6xsXhFR6Nxcchkq9lM/s8/k/3NN+z7+utWsRg1axZH3HUXcUuWIKVsnU8xlJcjTabW/WaDoXVYzGI0Ep2SwuxLLiFq1ixGzZqFX2SkS/l2B7W1b3AdJxcu/FwInAtkCCF2OI7dCcQCSClfBL5F6XnkoJhFXziA6zwBbJZSdniXOoTnViHEaOAJKeUK4BYhxGTgVFwpOEKImcAS4C0pZVW742HAOyim0gBWIcR/pJT39efizqLTZFlkO/UuRRly6+47l6KsFSI6OprKykosFgt2ux0fHx+qq6uJjo4mLy+PyZMns3XrVo466qhWJ4dpaWnMnDmTrKws4uLiKC4uJiQkhKamJjQaDXq9vvUXWmFhIQkJCezevZukpKTWMlryjIwMJk6cSEFBAZGRkdTW1uLlcGliMpkICgqirKyM2NhY9u3bR2JiIqmpqfj6+tLY2EhycjLp6elMmzaNnJwcYmJiqKysxM/Pr9d7anH/7qp7as9noPfUUi+MS+bDfDimyM6PJg1EwIy5FWSat6DzkEwSv5EjpxHcOJlIzXSM+4/ES2Mkq+7fTNl1FrLkN8Qx77Mzo4m53x+FzWLBeucPFFl8CKmp6f2e0neQdDAX89034VlaRPWiowl56nUyqmqYaLdTkJ3d7T0F+vmR/+23aH7/nZJPPsFWWgpaLaGLFyOWL2fONdewt7ISoddTr9EQGBNDeWUlfgaDS55TZV4ehevXc3DdOg7+/DPWpia0Xl7EL1vGuHPO4ajLLiO7rKxL25sycSLbtm0jMTGxz+eU0PKcIiNbyxiqtmcwGIiIiOiz7W3bupWxfn5s+uADtCUl5G7YgGhuxiYleg8P7FKi0enQaLVIIdB7eGC12fDy8cFkNjPzjDPwXrasQ730dE/5+fmueWe5aB2OlHITXde+dD5HAlcN8lJLgFt7+Xwj8FC7/TWAkwGj2iB6m/cRQrwInCCljO10/HPgZGA/sAM4HGXs71Qp5ef9JdErQWWybAOKSfanQohaKWVQu89rpJS9mt+kpKTIbdu29XYKVqsVnU5dI4xq4+RKPsf8BKlVEP2KZHc1yLmSn58p4VOP95lNGqEIjHhwIk/gXXwnVL8I0e/BL69C8QZY9i74JCHvXaL8N969TvHe3BuqKhWfZm8+DyVFMGEy/OcZOLL3xdK1W7dy8MUXKfvySyyVlWg8PQk79lgiV6wg4sQT8QjtuEh0sPXUWFFBRWYmlXv2ULFnD5WZmVTs2UO9Iw5MwJgxTDzxRCaeeCLjjjoKvRMLN9XWlqBnTmaDgaI//uDgb79R+OuvFG7Zgslhau4VHEzMvHn4hIUh7XbsLQHX2m3bHfst2wnHHcfCW3t7l3bkpNfrU/uaxO8LY1PC5B3bTnLq3CvEm4O+3mAhhKhDWch/Qw+fPwlc2BINQAhxDXBfX+/ezuirBc4Hvut04bEoYpMOHCalNDvG+VKBS4DP+0OgN/QwWVbWYmMuhIgCynsuwXm0/JpSE9TGyVV81pbATyVwnR88tUfANMmYlCbyPFLRYSOGSkoJZAGX413zhSI2ITfCxjeg6GdY9jZ4JcI9i7HaQX//L70602R3OrzyNHz6HpjNiruZh1+AZSeAtmcT6cbsbPbdeSelq1ej9fcn4sQTiVyxgvDjjuvVWaWz9WSsqaF461YqHIJSuWcPFZmZGKtaBxPQ+/oSPmUKcUuWEDF9OgnHHUdEYmK/51JGqi1Jux2L0YilqalL2rtrF2NGjcLS1ITZYKB0xw4Kf/uN8l27FGtEIQifOpUpp57KmPnziZk/n7BJkxBDaL25e/dul5RzCPpSWwNcIYTYIqVc1f4DIcRZKL4zv253eDaK0Va/0JfgRKOswWmPJY78eSmlGZRxPiHEuwxs7LBb9DJZ9iVwPkr37nzAJX7c1PRib4HaOLmCj5RwRxrE+oJhMwgPiQwX3LnIxHb2kkgG5YQQw2zGmcOh+DjwWQSpGVC4Fpa+Ad5JcO9RoPNAf996iO7GGabVCj98Ca8+Db9uAG8fOPNCuPgaxTCgF5hLS8m57z4OvvwyGi8vEv7v/4i7+WZ0/n2vwYGe66mhtJQDGzdy4JdfKPjlF8oyMpQKAbxDQpSX64oVhE2ZQvjUqYRPmUJATIxLXrBD3ZZszc1UZGa2WrmVpKVRnpGB2dC7Be9v7bY9AwKIOewwpqxYwZgFCxg9dy5eQUFDyrszXFlPh5jg3IgSeuA9IcSjtAVbSwCiUObTbwIQQngBY1GsiPuFvgTHjzbzuBbMRTFTXt/peC5KnARXoafJsoeAj4QQFwMHgO6cyvUbagt2Burj5Ao+Xx6ErVXwbBLccptEjgevSBv6CevRYiEYO3a8mGs/H3HwBND4wF4dFPwER70GXjMVsdF7wd3rSS02kNx+IXlNNbz/GrzxnOIdYMxYxRjg7IuV+DO9wFpfT96jj5L32GPYzWbGXHYZCf/+N56j+heauaWeag8c4MAvv7QKTNU+5beb3seHMQsWsPjeexmzYAGRiYn4hIcPqbNMV7al5sZGytLTO5hQl+/ahc3haVrv68uomTNJ/Mc/8I2IQO/j023aX1DA9NmzW/f9o6KGtPfiDFwbgG3wczjDBSnlASFEEorV74m0ucbJB94HHm6Zx5dSmmjrePQLfQlOIYrCtccCoFZK2TncqA5wmU1hH5NlS111nRao6cXeArVxGiwfu4R/74AJ/qBJB2OTAH/JqfON7NTkMZE8TFhYwBV4lz0Kpu1QOR/2r4PFL4PPbLhvibLG5u71MCqB5JZlMgX58PwjSkgAoxEWLIb7noRjT+p12AyUODAFL71E7v3301xRwajTTmPif/6D7wQnArM5IKWkMiuLgo0bKdi0iY0bNlBXUACAZ2AgY484gln//CdjFy0iavbsXtfDDAU6P7v6kpLWnkhFZiZWoxG71dpnsjQ1UZOX19oz8wkLY9SsWcy7/vpWC7fQCROcEo7BuzV1PVwZgM1V63CGC1LKahTDAecmvAaAvgRnG3CeEOJpx5zJfCAR+Libc6cCrvF8NwJQW28C1MdpsHw+zoeMWnjvcHjgCrsyYOuvYfLCb6nGij+NxJDCuAYTVD0Gtrmw+zdY+Bh4J8O9S8DLzyE2iiPK3Z99zLS1X8Gn7ytxalaeq5g2T+vbMYa02yn56CP23XUXxv37CVm8mEkPP0zQ3Ll9ftfW3ExJWhoFmzYpIrN5c+vci0dICAlLljD/5psZu2gREdOno+lD9IYKUkpq8/PZvHo1PnV1iqPM7dtpKC1tPScoLg5Pf3/F0qtd0np4oPfx6XJsxnnntYpLQEzMgHtmamvf4NoeTj9c2/xl0JfgPISy4CdLCLEXmAbYUVxhd8aJdB1mO2SgtoYP6uM0GD5WO9ydDtOCIMEAezI0MM/O+HFmKkbnEUs5XuiZa12JKDwCNLGwaaviiNN/Edy/FHwCFbGJiIPtW+HpB5n27WfK/Mwl18HlN0LU6N55NDbSsGsXhvR0Dr78MobUVPwTE0n59lvCjjuux5en2WCg8PffObBxIwc3baJwyxasjnAEIRMmMOnkk4k9/HBiDz+ckAkThj10spSSxvJyqvbto2rvXir27Gkd7jLV1gIgtFrCp05l/LHHMsqxmHNUUhKeg3AKOhiorX2Dq0NMHzpDagBCCF+U3s3faeuA7gc+BR7pr1eB7tCXa5t0IcTfgQdReja5wN1Syl87ET0WJfLnd11LOTSQkZFBYmK/PG0POdTGaTB83s+DvQb45Eh47DkLeOrAS0Pywg1osBNKKSnycryLbgRrNWR4Q+A0iL0cHjgafILg/9bBvny46lL4ZQ0EBlF27mVE3vGAElGzHaTdjjE/H0N6OvU7d7amptzcton6sWOZ8fbbRJ99NqKbHkhjRQWZq1ez+8MPKdi4EWm3I7RaombNIvmyyxSBWbgQv05zPEP53JobG6nOzqZq3z4q9+6l2pFX7duHuV1wNa2nJ5EzZjD19NOJmj2bxoAAFpxyilMm1MOF4WjfVrMZq8mEV2Cg05xchUNpSM3hMmwjMAXFjc12x0cTUTwKnCaEOMIx7DZg9GmYL6X8mo7mcN2d8wPgnAmPSjFx4vCH/e0LauM0UD4WO9ybDrNCYL4XnLZaC9MlGk8ImL+XCCoYx0ziqndD/VdQGgMNBlj2HPxnJfgGw+K74MKzIfV3CI+E//sfnH85QXoPcPgfq9m8maJ33lHEJSMDW4ubEiHwnTCBgJkzGX3eefjPmIH/jBl4jx3bZa7BVFvLns8+Y/eqVexfuxZpsxE2ZQqH33kn4448ktHz5uHZh7WaK56blJK6goLWIbCStDTKdu7EcPBgh/MCxowhbNIkEs85h7BJkwidOJHQSZMIjI3tMIxnNpvRu8BPmyvhinqyGI3UFRRQm59P3YED1Obnd9iuLylh9iWXcNJLLw0bJzgkzaLvAyYDVwMvSSltAEIILcoC+meAe4BrB3MRda0EG0EUFBQwoR+TxMMBtXEaKJ83c2B/A3y9BG58thm7RY8mQjJ58h58AxsZQy1zTachSpdAcwzsLYRjPoaXbwZjA+QFw8f/hDHj4OHnFfNmx2r5guxsJkyYQMWPP5J28slovLwImDmTmAsvbBUWv2nT0PXiyr65oYG9X33F7lWryPn+e2zNzQTHx7PwttuYfuaZREyf3q8hsv7Wk7TbqcrOVuLEtBOYlsWOQqMhbMoUxi5aRNjkyYS2CMuECeh9fIaE02AgpcRYVUXdwYPYrVY0Wi1Co1FSu+2CgweJi49v3bcYjZgNBprr6zHX13e73ezYNhQVUZufT2NZWYdrC62WwNhYgsaOZfwxxxA4dixjFi50mnuBw9BjsDgEBedk4FUp5fPtDzqE5wUhxCzgFNyC4xpEDoPvqP5CbZwGwsdkg/t2wmFhkIjkw3f0sMiO3aIlZmEuEZSzwH42PoWXg9TDtkJIvhPWfq8ETMvUQKCA596B5Wd08XMWGRlJ1c8/k7Z8Ob6TJzN33To8nAiiZTWZyP72W3atWsW+r7/GajTiP3o0c66+mulnnkl0SsqA52EiIyOxW60Ya2owVlVhrK6myZG33zdVV2MoLKQ0PR1LozI8rvXwICIxkakrV7bOs0QmJjotLL1xchVa5ova9yhqDxygrt12y/24CkKjUWLg+Pvj6e+PX1QUE088kaBx4wgcO5agceMIGjsW/+hoNIPwqODKejrE5nAiaRtG6w5pKOseBwW34DhQW1tLwAhNnvYEtXEaCJ9X9kFhE7y5EE553IK06YlJNFGeq2V8cjazCSeu7Csw7YRdOhh1NNSPgXX/hUINxKTA6jU9Bj4r/uknCs4/H5/x45nz00+9io2UkuJt29jxxhvs+uADTLW1+EZEMPPCC5l+5pnELlw44HUgUkoO/PILW599lpwff6S5lwWPQqvFOyQE75AQ/CIjmXXRRcoE/qxZhE+dOiQm0wNtS2aDgZLt25VInamplO7YQU1eXqvBRAu8goMJGjeO0IkTiT/mGILGjiUwNhath4fiZqbF9YxjW9psVJaXExIS0npc7+3dQVQ8/P3xDAjA098fnbf3sBhi1DoMLAYLDXY8Dy0rtTIUf5U9YZbjnEHBLTgOtDg0VBPUxqm/fJqs8J8MODISmppsbP9Qj26BhcJMX6YdmU6sZymLGuYiqu6Ecj8wh8D0G+GBk6BWAx5T4YPvehSb2i1bOHDBBXiPHs3cNWvwDO8+HEFDWRk7332XHW+8QcXu3ei8vJiyYgVJ559P3JIlg/pF3NzQwM5332Xrc89RvmsXXsHBJCxfTkRCgiIqoaF4h4Tg48i9Q0Px9Pcf9gWOzjw7s8FASVoaxQ5xKUlNbV2sChAQE0PU7NkkHH98l57FQCzdysvLiYiI6Pf3hhKu/J87xIbUvgIuE0KkAa9IKe0AQggN8E/gIsC5ibBe4BYcN4YMz++FMhO8sUiy8n4JEqzSg6mzdnLkuT/xN9th+By8Asy+kN0My16G/54LRjs0joPP1kBw9z2WurQ0th57LPqwMOauW9fFG4DNYiH7m2/Y8cYbZH/7LXarldHz5vG3F19k+hlnDNplStW+fWx9/nl2vPEGZoOBUbNmcfJrrzH9zDOpaWgY8ReptNuVuQ9HOILS/HxqNZq2MNAGQ+t2U0UFpTt2UJ2d3fr9gDFjiE5OZsa55xKVnEx0cjK+KhMHNeMQnMP5PxTvz88D9zqWwQBMQnHMnAPcPdiLuAXHAVMfoXJHAmrj1B8+9RZ4aBccGy15qtRC0xd6iBYsWLqBWf/cymSNjekHXgVrFey0wcKX4MW7oK4CykfBx+shovvxdEN6OluPPhp9UBAx77+P1+i2tTdlGRnseOMNdr77Lk0VFfiNGsVhN9zAzAsuIHxqH96k+4DdZiP722/Z+uyz5P74Ixq9nmmnncacq68m5rDDWod8TJWVg7pOb7CazTSUlFBfUtJt3rLdWF7eav7dEzQ6HZ6BgXgHBxMxfTpJ551HdEoKUbNnD4u4qK19gys5HVrrcKSUVUKIFJT4NqcAcxwf7QdeBf7nitDWbsFxIGiYnQQ6A7Vx6g+fp/ZAlRmmJdh5/DIdIEg+/3dmX/IHEZomLq72RTR8p/xuGncJrP0VDqRCcQC89wvExHZbbn1mJn8sW4bWx4e569ZhdQTMMtXV8fl557H3yy/R6PVMOukkZl54IQnHHTeoITNQhuTS33qLbS+8QG1+Pv6jR3PU/fcz+5JLug1U5ornZjYYKN+1i7KdOynbuZPyjAwq9uzp4Em6BUKjwTcyEv+oKPxHjyZ6zhx8IyPxDg5ujexp0+kIjorqEO1zuOZFeoLa2je4jpPg0FqHA+AQlLscaUjQr/9Ex0rUs4EJQChdfZ1JKeXFLuI2rCgrK1PVBD2oj5OzfGrM8OhuOHyU5PEXgQwNQUdVMv/KDQQKI1eaj8C7+ASo1oJ1JhgmweabocITXtsM8d2b7zbu28fWpUvR6HTMXbcOn/h4srOzaS4pYdXy5dTk5rLkP/8h+dJL8Rlk5Mb6khL2fPopmR9/zIFffgEpGbd4MUc/+iiTTj6514n9/jw3u81GdU5Om7Ds3ElZRkaHcM+eAQFEzpjBlFNPJSAmBv/oaPyjovCLisI/Kgqf8PA+XedkZ2cToyITe1Bf+waFkytwCA6pDQucFhwhxFyUBaC9/SdL4JAUnNjY7n9RjyTUxslZPo9nQp0Rtv4q4UctwsPOsU9+jY8wc7FtGaEHLgaLDfICYNJN8Og50KiFR9bDlOndltmYm8uWJUuQNhtzf/651bGmJTOTV887D62nJ+etW8fYI44Y8P0ZCgtbRaZg82aQkvCpU1n0738z/YwznB6S66meWvyaFW3ZQuGWLRRt2ULpjh2tFl9CoyF00iRGz53L7H/+k8gZM4hITCQwNnbQPRG1tSX4c3NSfKmZXVLWUEAIcd5Aviel7HdIgvboTw/nccADJRzAusG6OFAbWsIgqwlq4+QMnwoTPLFT4rELzAcElElir9hPdHgRZ8rZxBX8B8zZkAkkPQX3XwA2CTethpT53ZZpPHCAP5YswW40Mnf9evynTkVKyeb//Y+1d9zBqJkzOfPzzwkcwMuirqCAzE8+IfPjjyn8TYnOEpGYyOJ77mHqypUDmvdpqSdTbS1FW7dS5BCXoj/+UOZWAJ2XF1HJySRfdhmjZs4kMjGR8KlT0Q2RZaLa2hKol5MrcAh4i34TpYPQn18ykgHEwGmP/ghOMvBfKeXqwVxQrVBbwwf1cXKGz91boHELYIRgaqnz8WfhFT9zApHMKl4Ljd9DNjDhQbjvGtA0wxlPwrJTui3PWFjIlqOOwmowMHftWgJmzMDS1MSXF1/MrlWrmH7mmZz82mtOLYw019dTmZXVGr45f/16iv74A4BRM2dy1AMPMHXlSsImTXK+UtqhrqCA/WvXcmDDBn7ZsoXKrKzWz8ImTybh+OMZPW8eMfPmEZGYOKwhCtTWluDPz0nlQ2pHjcRF+yM4BqDrjOWfBGp1la4mTn3x2VIAL7wvwQITjskn+6Y44q7byymhFSwpD4CaJ5UAFuOfhwceBl0dzL8czrquS1lSSkpWrSLrxhuxNjYyd80aAmfPpq6ggFWnnELpjh0sfeghvJYu7SI2jRUVraLSErq5cs8eDIWFredo9HpGJSWx9MEHmXLqqYQOYH7DWFND/vr17F+zhv1r1rSaFXsEBTHu8MNJPOccRs+bx+g5c4Y9cmVnqK0tgXo5uQJqn8ORUm4YiesK2YfpZOuJQrwAREsplw8tJdcjJSVFbtu2baRp/KmRWQRzHpA0mSH4gjp83qyjdPNo7tryDP+yxaAvOkP5yRL2Htz3GFhTIe5IePjnLmU17N1L5lVXUbV2LYEpKUx/9VUCkpI4sHEjH516KjazmRXvv8/Ev/0NUNac7F+7ltSXXiL/5587WHLpfX0JmzyZ8ClTWkM3h02ZQnB8fL97GFaTiYLNm9m/Zg15a9ZQnJoKUqL39WXc4sXEL1tG3NKl/fa95oa6IIRIlVKmDKaMpBSd/HZbkFPnxoiqA0B7W/qXpZQvD+b6akV/eji3AT8IIZ4BngT2S2fV6hCAWn9tqYlTT3zKDbD4QUmTBTjexjF+n/Ph9xcw9aad3GKLR1+wAuXH3v1w97UgKyFoHDzwQ4dybEYjuQ8+yP6HH0br7c3U554j9rLLEFot2156ie+uvprg+HjO/OILwiZPprGigq/vv5+yb7+lJjcXn7AwJp9yCuHTphE+dSrhU6YoAcKcWNVvt1ppqqqiqaKCpspKGh15y37lnj0UbNqE1WRCo9MRc9hhHHn33cQvXcrouXPRenj0WU8jCTcn5+DKHo6n80YDlYMVuEMF/enh2FEmjXqDlFKqbm2Pu4czdJASTnwcvtslkXNg+coP+f3mRVRuC2f9d19zRMlp4GWDDdNh8y5oBmasgDvfUqJ3OlDx/ffsvuoqjPv3E33OOUx+9FE8R43C1tzMd9ddR+qLL5Jw/PGseO89yjMy2Pbii+z55BNszc2MXbSI5MsvZ8qKFej6cMEvpaTw99/ZtWoVJdu2tYpLi2fm7uAZGEjQ2LGMW7KE+GXLGLtoUZ8hCtw4dOGKHs6sFI1ct805A5AQYRz09Q4V9Ecc3qZvwXEphBCvo0QSLZdSTnccCwE+BMYB+cDpUsqe3xZOIj09naSkpMEW41KojVN3fJ79Cb5NByYJIhOLoBDK1kRz5E27OWL/SgizwydA5j5lzfK1z8N5V7R+31RUxJ7rr6d09Wp8J01i7tq1hC5ZAkDxtm18c+WVFG/dyrzrriNw7FheX7iQyj178AoKIvnyy/E78kiOWLGiV95SSkp37GDXqlXs/vBD6g4cQOvpScxhhzFq5kx8wsPxCQtrzX3b74eGdui9DLSeRhpuTs4hPT3dZWVpbeqdwxkpON3DGQkIIRYBDcDb7QTnf0C1lPIhIcTtQLCU8rbeynGmh2O1WtENckW6q6E2Tp35ZByE5LslliDQLzZz/vGv8OX5K6neGcre/04kPv4A/OoDP0bDrznwxJtwhuLh3G61cuCZZ8j+v/9DWq2M/9e/iLv5ZrSenjRVVrL2rrtIe+UVvENCiExMpPD337GaTIyeN4+Uyy9n2umno/fx6bWOKjIz2bVqFbtWraI6OxuNTsf4Y45h2plnMnn58iELray25wZuTs7CarWi1+sH3eOYnSzkxl+dc9Dq52V393DUACnlL0KIcZ0OLwcWO7bfAn5GmV8aFHJycpg8efJgi3Ep1MapPR9jM6x4VmLRAMl2lh/1ETXbQyj/OYqTL/1UERtDEnznB7//Bs+8Ayv/AUDNb7+x+4orqE9PJ/yEE5j6zDP4xMdjt9nY+sILrLvrLswGA8FxcdTs30/xtm0kXXABKY41Kz1xAqjOzWX3hx+ya9UqyjMyEBoN4446igW33MKUFSvwCQ0d1npSC9ycnENOTo5LyhESdDa7S8r6M2FAPtKFEH5CiBghRGzn5GqC3SBSSlni2C5FCRzUHcdLhRDbhBDbSkpKqKyspKSkhKKiImpqasjNzcVoNJKZmYndbsfgiF/SMmmYlpaG3W4nMzMTo9FIbm4uNTU1FBUV0VJefn4+DQ0NZGVlYbVaW7vjLWW05BkZGZjNZrKzszEYDBQUFFBeXk55eTkFBQUYDAays7Mxm82tMdVTU1OJiYlpLSM9PR2r1UpWVhYNDQ3k5+f3eU9paWkuvaf2fM55tpqcEgEzYNHRa5ngl8PmhxejDzPz/IlXY9UkYbxfg9zyOwV3PAgr/8EfP/zAzvPP5/cFC7BUVRH8+OMkff45Bc3N7F2zhueSkvj2yiuRgLTZaKquZuaNN3JFTg7x11xDxIwZXe6prq6O4rQ0Pr76al6YOZNnEhJYd9ddCC8vFjzwAOfv3MmiV19l0llnUVBRMSTPqX2enp7OqFGjRvQ5dXdPNpttUPc0FG3P399/WP+fnLknT1eF4ZagtTqX/lKQUjqdgDOBXSg2R92m/pTn5DXHAbva7dd2+rymrzKSk5NlX8jLy+vznOGG2ji18Pk0VUrOlZK77HLi/p3yX/Y75Q3fPCiJkvKMRz+QMjNKyjOSpYzSSvnFR9JmNsvc//1P/uDnJ7/z8JBZd9whLfX1Ukop60tL5Wfnny/vAfmAl5e8B+SDQUFywwMPSFNdXbc8rM3NMnfNGvntNdfIR0aPlveAvEcI+dqCBXLzI4/I2gMHhqtKuoXanpuUbk7OIi8vTwLb5CDfW8lJSFnuXHLF9Q6V1B9faqcA7wP7UALxXO7Y16G4s94JfDNw6XMaZUKIKClliRAiCih3RaF+fn59nzTMUBsnPz8/imvgjBck+AmCTy3jjHEfMC1jN1c//gqekUaeXXoTPBYKG9PhlY8oF15kJSbSuG8fESedxOQnnsB3/HhsFgu/P/kk6/71LyxNTQBovbw44l//Yt4113SZXzEbDOR8/z17v/iC7G+/xVRbi87Li5gjj2TJvfcy8cQTu/XcPBJQ23MDNydn4TJOEg6h6ATDhv7M4dwM7EFxceOHIjivSynXCSGmA5uB/7ieYhd8iRJb+yFH/oUrCrVYLK4oxqVQGyez2cKCFyQWi0B/ShOXzniB4J9rueTut6nP8+fO+/9D2Ade8P1eTPc/za4XXqfim2/wnTiRlG+/Jfz44zEbDGR9/jk/3nwzNbm5AHj4+7Pw1luZe801eAUGtl7PajYrk/4ffEDeunXYLRZlrc3f/86k5cuJX7aMKoOBqKiokaqSbqG25wZuTs7CZZzcgtMt+iM4M4AHpJQmIUSLLxEtgJRylxDiZeAOXCQAAEKID1AMBMKEEIUoEeceAj4SQlwMHEBxJjpo2O3qm+BTG6dLvwmgoEDAYTb+uew5ar4M5bF/34WwwSuP/ZN/7v8e+V4VRUeezK7Lr0Pj5UXCAw/AjBmkb9hA/j33UJyainSYi+p9fVl4223Mu/baDkJjqq1l24svsuXpp2koKSEkIYF5113H5OXLiZk/v4MrfruLYtC7Emp7buDm5CxcxsktON2iP4Kjpc2XmtGRB7b7fC9wBS6ElPKsHj5a6srrAPg44fxxuKEmTk9ts/DtZh8YLTn5Hx+S+d5Mfrl3KUGja/jpviuZ/XYqcm0l6TY/St7/BM2MGRTp9Wy55x7sVisavZ7QiRPReXlhM5tZeNttLLjllg5CU3vgAFueeoq0V16huaGB+KOP5pS33iJ+2bIeXcWoqY5a4ObkHP70nNzLcLqgP4JTCIwFkFIahRDlKMNrLd6jJwGNrqU3fKiuriY4OHikaXSAWjitK2vi+te8wROmXJxGxhMp5L00kelzMthw65cEXroeW2UlaTWSAmst+YAxM5PoOXNYcOutjJk/n71ff03aSy8RPm0aK957j1HtFvyVbN/Ob48+yq4PP0QIwfQzz2T+TTd1MYHuDmqpo/Zwc3IOauXkErh7ON2iP4LzK7AM+D/H/pfA9UIII4p59VXAV66lN3yIjo4eaQpdoAZOu5rqOPZZH2gU+J1RSckDY6n9MYwzj/uYV/9Wiudp92BptvJ7HRTFjCP8tNM4ZckSYg8/HA8/P0q2b+fTc86hcs8e5l1/PcsefBCdlxdSSnJ//JFfH3mEvLVr8fDzY95113HY9dcTOGaM0/zUUEed4ebkHP7UnCSoOP7aiKE/gvM88HchhLeU0ogS93oucI/j890ohgWHJPLy8pg6gGBbQ4mR5pRtqeDwtz2wHtAjJjdjfdybxr3ePHLqvzivYQveV6yhwQZrrL5Mff8NFp92Wut37TYbmx5+mPX//je+4eH848cfGX/00Ui7nYwPPmDTgw9SnpGBf3Q0yx5+mORLLx2QC/+RrqPu4ObkHNTKySVw93C6xaBd2wghZqCMVu6RUqpvFhDnXNvY7XY0TngVHk6MJKe99iKO+cZKweexYJNoU214GptZvfwMpq9Zy5hSE6XNkHH8mSx66eUOzizrCgr47LzzOLBhA1NOPZUTX3oJn9BQDmzcyI833kjxtm2ET5vGgptvJvHss/vtq6w93M/NObg5OQe73Y5Wqx20q5mUKUJue9O5c8Vh/GVc2wz6aUspd0opd6tVbJzFjh07RppCF4wUp532/Zy4roGCNWOhCMR6iBYl/HH2HJLXfsOYUhM5Gh8sX2zg6Pc/6CA2Ge+/zwszZlCSmsryN97gtI8/xlRby0crV/LmokXUl5Rwyttvc8XOncy84IJBiQ24n5uzcHNyDi7j1NLDcSb9hdDvHo7DoeYxKC5lHpNSZgkh/IDZwE4pZa3LWQ4S7vAEzuMPmcnZP2jJ/WIS7JSQJ1gQt4kvTl6O/qtqArMhf/YRjPn8J7Tt3ICYamv55sor2fXBB8TMn8+Kd9/FOySEXx54gC1PP41Wr2fh7bez4KabnAoH7YYbIwVXhCdImSzkNidDqIkj3T2cLhBCaIUQHwLrgTuBi4CWGTYr8DlwpasJDhdcFXjJlRhuThtkGqd/7kHuu5NgoyI2/zjsTdafsQSv9dX4Z0P9zfcy7rtfWsWmbOdOvr/hBp5OSGD3Rx+x+L77OG/tWrK//ZanExL47fHHmXHuuVyTnc2R//63y8XG/dycg5uTc3AZJztgcjL9hdCfAGx3AvcBNwHfo3gdWCalXOf4/DVgkpTy8CHiOmC4ezi9QyL5Vm7m4nfGUfZRDOxW4rE/f9RlXJr8OuafQfebJ5o3PkEc/TeMNTXs+uADtr/+OiWpqWj0eiadfDILbr2VpvJyfrz5Zqr27mXcUUdx7OOPO2Xe7IYbaoFLejgThNz2tJPXO8Hdw+kO56HEpXmKjvG3W7AHGO8SViOAFu+2asJwcLIj+ci+hrOemELZizGQBhGRpaSdNYNLk1/HlirQ7w5DfLWRH/bl88lZZ/FYVBTfXnUVduv/t3fn8VHU9x/HXx+ScIlc5RDkVPFEEMVQ631iFRV+eIAHVOtRRQreilrtD61Si1qPeoAHHki9qPjDAxBRaT0I4YiEKARCIAQCkhhCwpLNfn5/7IRG2MCETHYm2c/z8djHsrOzs+/ZBT77nfnO9xtm4JNPctv69Zxy//3MHTeOty64AFQZNmMGIz77rM6LTaJ+bzVlmdzxNJOdw9lNTbpF9wAm7uH5IiBYV3HVwDEB/BVe15nCVPDctk+4/f4z2PFxM/gZBp42k/dTL6HRT2XoVNBwT7495zy+vWgoxWvX0rRNG4697jr6XXMNnfr14+fcXGbddhtLXn+dZm3acO5TT9H/D38gKSWlTrNXSsTvbV9YJnc8y2TdomOqSQtnK9B2D88fAmyqXRz/ZGVl+R1hN3WZKcQO7s//iDFXn8uO95pBmfL4yNHMPGYQyR9tp8kkIauoK39bsIrZTz1Lh6OO4tePP85t+fmc9/TTtO7Rg9l33cXThx7K9//8J7+5/XZGr1zJgNGj41ZsIPG+t31lmdzxLJOyh0lcdrklkJq0cOYDVzpTPP+CiLQh2ongE6+CxVvPnj39jrCbusq0jTJGpn/Je38cBKuERp0q+ObS4+m7bBFJU6CwZQfe+WkjJZHt/Gb8QxwzciQtu3ShrKwMVPnPxIl89fDDbC8qou9VV3H6+PG06haPufd2l0jfW21YJnc8y2QtnJhqUnAeJlp05gKvOsv6ikgv4G5gP6IjOddL69ev5+CDg3UKqi4yFbKJ099Yx5KHz4GfhaYnbSP3iK60f6eQ7Zv241/byvn+py0MuPVOTh43bufgmhqJMP8f/2DZ00/z85o1HDxwIGdNmPCLMdH8kCjfW21ZJnfWr1/vzYaUhOuB5obrgqOqaSIyFJgMvOIs/hsgRCdBG6Kqmd5HjI+2bfd0tNAfXmdaVPY5Z97bm8J3+4Eoba8pYM2WHjR/q4wFW1swZ2MJhwy9mFETJtDmoIN2vm7VnDnMvvNONixaxAH9+nHh5MkcdNZZnmbbV4nwvXnBMrnjWabKQ2rmF2rSwkFVZ4pID+Bs4AiixWYF8KmqlnofL35KS0sDN3KtV5nWVJRybfpS5tx6KmQ3gnZK15tz+P67vmz9PMKra6DxMYdx+dtP0P3kk3e+Lj89nc/uuYfsWbNo1b07pz31FKeMGoUEaDiShvy9eckyuVNa6tF/Y3ZILaYaFRwAVQ0B/+fcGoygjekEtcukqszboUzcup6ZkzvCpAHRWYxSlYOGreSbTwbw2XtbWb1/J858+VH6XHkl0qgR2zZtImPqVJa+9hr56ek0bdOGcyZO5PhRo9i0ZUugig00vO+trlgmdzzLVLOC005Eql4o+KKquhynoH6pccFpqFLi2LPKrX3JtDWivF4W4ZltIZZvTIK7O8OiRrCfwu+Vq/u8wrj372LKh6X0v/0Bzr/jDholJ5P1r3+xZMoUVnz0EZFwmAP69WPgk0/Sd8QImjm/QhvKZ1TXLJM7DT6T+0NqmxPlws89FhwRmVvD7amqej4bZzyUlJTQrl07v2P8Qk0y/RCO8Oy2CK+WhtmqQrM55fCXZlAscCSc/btPeTTzNrJvyiRj4BVct/wRSvLzmX3nnSybNo2yLVtoccABDBg7lr4jRtDx6KNrlSdeLJM7lsmdkpISbzZUObSN+YW9tXBOA8qBHS63V7u5DnwUtL/44C5Tdli5+edyPgkpKRrmQM0hdPdBlM1rBcnQ6ap1vNHpCspu/ZIFBx7DqZOms2n5ct445xw2Z2WR3LQphw8eTN+RIznorLNolFz9X4n6+hnFm2Vyp8Fnsk4Du9nbAcsw0Y4Bc4ArgFaquv8ebi3rPHEdWbdund8RdrOnTKrK5G0V9C3YzvzSrQzY72s6Ls0n5/wj2TG3GY26hXni/tHMntuVRX9eQqPhf6RV9+68PXQoc8eNo3n79lwwaRK3bdjA0Lfe4pBzz91jsdlbHr9YJncskzueZbLpCWLa4+CdItKe6BhqvwOOJNr9+TXgZVX9IR4BqyMi5wJ/B5KAyapa7TVAbgbvDIfDJO/lP9x4qy5TQYVyXVE5M0JKt+Q1tEouJOuePpR/3AxU6XPpYj7+6QSmvxmiWepv2LZpE1tWrKB5u3b0u/Zajr32Wtruw/UP9ekz8pNlcieomVJSUmo/eGdn0bTr3K0r/2uDdwKgqptUdaKqHg2cAHwAXA9kisjXInKtMxdOXIlIEvAs8FuihXC4iNRqrtply5Z5Ec1TsTLN3F5B7w0lfBQqp1vSarZM7UDGwFTKP2xGcrtyXrxrJA+/fixTP25JUePmrP3Pf2jaujWDp0zhlrVrOeuRR/ap2FSXx2+WyR3L5I5nmayFE9O+TMDWFBgKXA2cDpQCN6rqG97HqzbDCcCDqjrQeXwPgKo+Emv9hjA9wbaIckthGZNCSTQpKEVfa8KOz5rBJoFGyqGDM3kh/yy+nrGJHeEKGqWkcPTw4Rw/ahQHpqb6Hd+YesOT6Qk6iaaNdPl+E6yFUy1V3a6qbwIPAJ8RHdLmoD2/ynMHAmurPF7nLNtJRK4XkTQRScvPz2fz5s3k5+eTl5dHYWEh2dnZlJWVkZmZSSQSYfbs2cB/J2BKT08nEomQmZlJWVkZ2dnZFBYWkpeXR+X2cnJyKCkpISsri3A4zJIlS36xjcr7jIwMQqEQK1asoLi4mNzcXAoKCigoKCA3N5fi4mJWrFhBKBQiIyNj52srb//ZEabHmk1M+ioZuTWJ0GVt2DGtObJV6T14Ea/dMZhRU4/ji/c30LhDR3qPHcvY3Fy6jRnDgampnu1T5a02+1T1fsmSJYTDYbKysigpKSEnJ2ev31Pl8PGV25g9e7bv39Ou+7RgwYJa7VNd/N2bN2+er99TrH36+uuvff2eYu3T/Pnz8YRNwBZTjVo4ItIJGEn0nE4vYD3RczrPq+raPbzUUyJyMXCuql7rPL4KGKCqN8davz62cMoIsUxzuDET0mYdDB8mwWqBECS3Lefy817nocwxzJtVwqpy6Ny/PyeNG8dhF1yw15P/xpjqedLC6SCadpnL93vGWjg7iUiKiFwsIjOBXOBBYCkwCOiuqvfGs9g48oCuVR53cZbts6BMdxtiB+/u+IoT5ywhdexBpP3+MHg8CbKE5p1LmDDiDor2b8LFT/+eV2aWoCefwXVpaVy3YAFHDBlSp8UmKJ9RVZbJHcvkjqeZ7BzObvbWS+0p4HKiE6tlAC8Db6jqlvjEqzZXMvAjcCbRQrMAuFxVY57xqw8tHEX5piKTu778mX+/ejyRpSlQoIDQ8pgixh93HzfOeJZvV8JXpdD9vPM5ffx4OvXr53d0YxoUT1o47UXTBrt8v8mJ08LZ28/hm4mOwPUWkO6s/zsRqW59VdUnvItX7ZuEReRm4FOi3aJfrq7YuJWRkcHRMa6uj4fVkVxGf72Kj5//DZHFjeEnoDG0OG0rfzn0LkZNe4Fl85VnSqDzRUO4+oEHfJkWwM/PqDqWyR3L5E7lOZ9as8E7Y9pbCydSw+2pqibVLpL33LRwQqEQTZo0iVOiqC26mlvSMnnzibOpSG8MxUBz2P+MIh7sei9/nPIC+Vsq+HQrtBp8Mafcfz8d+/SJa8aq/PiM9sYyuWOZ3AmFQjRt2rT2LZxfiaad725ded1aOJVOj0uKAMjNzaVXr151/j6KUsIC7luYyT8eG044rWe0Y3kraDusgInltzDinbfYXK5M3wZy4SVc+Kc/0aF37zrPtjfx+oxqwjK5Y5ncyc3N9WZDNgFbTHssOKr6RbyC+K1jx451un2lgq0VbzBh8Xr++sgthBemQghoBwcMz+MfG25k4OsfkhGCSduF9pdcTv8xY+gZoGto6voz2heWyR3L5I5nmWwCtpis/6yjqKiIli3rZii4kGby0vIXGHvvBMoXNo0e2+2sHHhpLs9vuI5ez85h8Tbl78nN6Hf9DVw2diytu3f37teWR+ryM9pXlskdy+ROUVGRNxuyczgxWcFxNG3a1PNtRtjO9wV/4Zw/38DGD5+EsEAPpcuI1Tz90w3s9+g8FheFWd6uPQP+NIYLb7yRZlWmuK2LTLURtDxgmdyyTO54lskKTkxWcOrITyXzGPp8hC+e+zOUCXRSOtywjolJfyA8bg5LCkK0OfhgTn/kDvqOGEFKs2Z+RzbGeMkOqe3GCo5j+3ZvzvDtKC/mltc+5rnHL0YLk6ANpAzdzt2n38cBNzxL9trtHJiayqXP3cVhF11Eo6TqO/V5lckrQcsDlskty+SOZ5lsAraYrOA4WrduXavXRyLw9LuzuWviiYTWXgYtgZMinH7DJwy7YyT5j22mxZlnct7r99P9lFPYw7VMnmXyWtDygGVyyzK541kmO6QWU40H72yoNm7cuE+vU4X3Ps6hwxmbGTv2bEJFzaGP0nbsBiZ2OJPTzjofSerBiLlzGTFnDj1OPdVVsalNproStDxgmdyyTO54mqnC5S2BWAvH0a1btxq/Zt36ck67Jo/spT2ghcJRIIeHufKM5zn8qjtI6dyVy2bM4NBBg1wXmdpmqktBywOWyS3L5I5nmayFE5O1cBw//vhjjV9z6ti1ZGd2hz4KqUL3S7N5YlsqfW58iEEvTObmrCwOu+CCfSo2+5qpLgUtD1gmtyyTO55lsgnYYqrxBGz1UV0M3jlp6gKuv/14OERJ7lvO9ef8nV63/y+/vukhjr/pJpJSUjx9P2NMfHgyeGdz0bRDXL5fRuIMbWMtHEdNhiXfuq2cmyb0hf0guX+IqYf0Z8i/8xm1uIBfjxnjWbEJ2vDtQcsDlskty+SOZ5lsAraYrIWzD84Y9R2fT0+F45VHh45gzEXP0LRVK8+2b4zxjyctnCaiaV1cvt8qa+EkHLe/bL74ajGff3A8HADHXvc1w7pfVWfFJmi/AIOWByyTW5bJHZuArW5ZC6cGKirCtDqjlG2rW5IyopTXdwznsr9+4EFCY0xQeNLCaSya1t7l+61nDbC5yqIXVfXF2rx/UFkLx7FkyZK9rjPijulsW9ESjlLuOGUMQ8e/7XumeApaHrBMblkmdzzLVLNeaptVtX+VW4MsNmAtnJ3C4TDJydVflvTv6dM46bbLoAkc+epC3i1YxxEXDPY4ac0yxVvQ8oBlcssyuRMOh0lJSal9CydZNK21u3XlJzuHk3BWrlxZ7XOFa3IYMuU02C4kX7KDWxeOr/Nis7dMfghaHrBMblkmd7zKFFEo2+7ulkiC9fPCR126xO5SEtq6lXsem8am7+6CQ5UbRz/KZcSnxVtdJr8ELQ9YJrcskzteZVKF8gQbtsaNQLZwROQSEVkmIhER6b/Lc/eIyEoR+UFEBnr1nps3b95tWaSign/+YTgvzL0TmkCPv/7AkG8raNE+PjMVxsrkp6DlAcvklmVyx6tMqhAOu7slkqC2cL4H/gd4oepCETkSGAYcBXQG5ojIoapa698SLVq02G3ZrNtv5255DooaIVeUc0/H+zm1/7TavlWtMvkpaHnAMrllmdzxKpMC5QlWTNwIZAtHVZer6g8xnroImKaqIVVdDawEUr14z/Ly8l88XjhpEjPkV+R/0RW6KMP+NIVTy0fSqFH189d4bddMfgtaHrBMblkmd7zMZJfh7C6QBWcPDgTWVnm8zlm2GxG5XkTSRCQtPz+fzZs3k5+fT15eHoWFhWRnZ1NWVkZmZiaRSITs7GwgeuHX6rlz+fK5x3n+y3Gg0HZ8PoMLZtLhVyeSl5dH5fZycnIoKSkhKyuLcDi8s0tl5cVjlfcZGRmEQiFWrFhBcXExubm5FBQUUFBQQG5uLsXFxaxYsYJQKERGRsbO10YikZ3bWLJkCeFwmKysLEpKSsjJydnrPqWnp/8iR3p6OpFIhMzMTMrKysjOzqawsND1PlXNs6/7VPXei33Kzs6u1T558T3tuk/l5eW+fk+x9ikvL8/X7ynWPhUXF/v6PcXapy1btsT436TmIkCpy1si8a1btIjMAQ6I8dS9qvqBs8484HZVTXMePwN8o6pvOI9fAj5W1Xf39F5uukUXFhbSpk0bfvrxR15KHcDjQ5ax8dPO8NsKHnvhj4zYMZoOzQ+v8X7WRmWmoAhaHrBMblkmdwoLC2nbtm2tuyn3FtF3XK57JNYtus6p6lmq2jvGbU+X7ucBXas87uIsq7XKXzYbly5l2c2PsPHLztBGOe2hWfTeHIl7samaKSiClgcsk1uWyR2vMilQ7vKWSOrbIbUZwDARaSIiPYFewHdebLhz584ArN6wltfm3QBl0OS2Eq7pMJkTOz7oxVvsc6agCFoesExuWSZ3vMpk0+HEFsiCIyJDRGQdcAIwU0Q+BVDVZcDbQCbwCTDKix5qAKtXrwbg4x1tIAc4LsLYyx+jR+kx7C/x6QZdXaagCFoesExuWSZ3vMpkLZzYbGgbRyQSoVGjRvR5cBUZU3ty5NuLuO/ohxmS9DJN8WfqgcpMQRG0PGCZ3LJM7kQiEZKSkmp9TuVwEX3J5bon2TmcxLN48WIA7r54MW1nrGd0n79zGAN9KzZVMwVF0PKAZXLLMrnjVaYIUObylkishbOLBW/eQ97l6YQiLfmfpCmk0LyO0xljgsSL6Ql6iehTLtc9z1o4iaeyf/6vhh3CNtmfo5Mu873YBG2CqqDlAcvklmVyx8tM1mlgd9bC2cVGMlnG+5zKnSTRuI6TGWOCxosWziEi+leX6w61Fk7iqbwyuiNHcgb3BaLYVGYKiqDlAcvklmVyx6tM1kstNmvhOILaYyZImYKWByyTW5bJHa96qR0kouNdrnultXAST1ZWlt8RdhO0TEHLA5bJLcvkjleZrJdabEGdniDuevbs6XeE3QQtU9DygGVyyzK541WmykNq5pesheNYv3693xF2E7RMQcsDlskty+SOV5lsaJvYrIXjaNu2rd8RdhO0TEHLA5bJLcvkjleZrIUTm7VwHKWlwZuZImiZgpYHLJNblskdrzJ52cIRkZdFpEBEvq/meRGRp0RkpYgsFZFjPdmJOmAFxxG03jIQvExBywOWyS3L5I5XmTzuNPAqcO4env8t0ZHzewHXA8/tU+g4CN437pOUlBS/I+wmaJmClgcsk1uWyZ0gZlLVL4E9TdRzEfCaRn0DtBaRTvFJVzMJcQ5n4cKFm0VkzV5WawdsjkeeGghapqDlAcvklmVypx3QvbYbyYdPH4xuy42mIlL1QsEXVfXFGrzdgcDaKo/XOcvya7CNuEiIgqOq7fe2joikBe3iq6BlCloesExuWSZ3nEw9arsdVd3TIbCEZYfUjDGmfssDulZ53MVZFjhWcIwxpn6bAYxweqv9GvhZVQN3OA0S5JCaSzU5ZhovQcsUtDxgmdyyTO4ELpOIvAWcBrQTkXXAA0AKgKo+D3wEnAesBEqBq/1JuncJMXinMcYY/9khNWOMMXFhBccYY0xcJHzBEZFzReQHZ1iIu/3OA3sfyiLeRKSriHwuIpkiskxExgQgU1MR+U5EljiZ/ux3pkoikiQii0Tk//zOAiAiOSKSISKLd7newzci0lpE3hWRLBFZLiIn+JznMOfzqbwVi8hYPzM1RAl9DkdEkoAfgbOJXiy1ABiuqpk+5zoFKCF69XBvP7M4eToBnVQ1XUT2BxYCg/38nEREgP1UtUREUoD5wBjnSmtficitQH+gpaoOCkCeHKC/qgbmIksRmQJ8paqTRaQx0FxVi3yOBez8fyEPGKCqe7tg3NRAordwUoGVqrpKVXcA04gOE+ErF0NZxJWq5qtquvPnrcByolcy+5lJVbXEeZji3Hz/9SQiXYDzgcl+ZwkqEWkFnAK8BKCqO4JSbBxnAtlWbLyX6AWnuiEhTDVEpAfQD/jW5yiVh64WAwXAbFX1PRPwJHAn0fEbg0KBWSKyUESu9zsM0BPYBLziHHqcLCL7+R2qimHAW36HaIgSveCYGhCRFsB7wFhVLfY7j6pWqOoxRK+sThURXw8/isggoEBVF/qZI4aTVPVYoqMKj3IO2fopGTgWeE5V+wHbgKCcP20MXAi843eWhijRC069GRLCb855kveAN1X1fb/zVOUcjvmcPQ/hHg8nAhc650ymAWeIyBv+RgJVzXPuC4DpRA8l+2kdsK5Ki/RdogUoCH4LpKvqRr+DNESJXnAWAL1EpKfzy2YY0WEiTBXOCfqXgOWq+rjfeQBEpL2ItHb+3Ixox48sPzOp6j2q2sUZ/HEYMFdVr/Qzk4js53T0wDlsdQ7ga+9HVd0ArBWRw5xFZwK+dtSpYjh2OK3OJPTQNqoaFpGbgU+BJOBlVV3mc6yYQ1mo6ks+RjoRuArIcM6ZAIxT1Y/8i0QnYIrTo6gR8LaqBqIbcsB0BKZHfzOQDExV1U/8jQTAaOBN54feKgIwHItTkM8GbvA7S0OV0N2ijTHGxE+iH1IzxhgTJ1ZwjDHGxIUVHGOMMXFhBccYY0xcWMExxhgTF1ZwjDHGxIUVHGOqEJFnRCTPudjVGOMhKzjGOJwiMxj4QO0CNWM8ZwXHmP86nuho4f/yOYcxDZIVHFPvicgEEVEROVREnnIOiZWKyGwR6eqsc5UzPH+pM8Pr4BibGgIUER0ItHLbLUXkXhFZKiI/OzNBZorIs3HZOWMaEBvaxtR7IjILOInomFxLgK+A44BrgQ+BXOBk4A2ic8PcDTQDeqjqpirbyQIWquoVzuMmwHdAd+AVogNMNgeOBrqp6jnx2D9jGoqEHrzTNBj9iBaQv6jq1MqFItIfuACYSXSK5XJneTnRidL6AJ85y44ADgPuq7LdC511BqrqrLrfDWMaNjukZuo1Z0rndsCMqsXGUQiUA9dUFhtH5eRx4SrLBgPbgaojKbdx7lNFxP6tGFNL9o/I1Hf9nPt/xniuNzDPmXisqiOc+x+qLBsCzFHVkirL3iV6iG48sF5EJonIICs+xuwb+4dj6rvKmSK/qbrQ6SzQftflVV6z3pkIrLKV1J9deqep6hai54LOJVrQziZ6Tmi+M4+LMaYGrOCY+q4fUKSqq3ZZXlmI0qt5TdXlg4l2JthttldVrVDVT1V1DHAw0Y4HJwB9a5nbmIRjBcfUd7sWj0rHOfcLqy4UkR5A211eMwSYv0uPtfa7jjagqhVABdHilFfr5MYkGOulZuotEWkLdAPejvH0scAmVV0bYzk4BcfZxinAnbus9zfgJBH5AFhJ9MfZQGAQ8JiqrvdkJ4xJIFZwTH1W2WGguhbOIhevGUT038H0Xdb7jGjvt0uJngvaQvQ6nMGq+kEtMhuTsOzCT5PQRGQ60FNVj/E7izENnbVwTKL7GpjkdwhjEoG1cIwxxsSF9VIzxhgTF1ZwjDHGxIUVHGOMMXFhBccYY0xcWMExxhgTF1ZwjDHGxIUVHGOMMXHx/wpq6WmGReIdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean_snr = scores_po.mean(axis=2)\n",
    "\n",
    "fig, ax = plt.subplots(1)\n",
    "\n",
    "# setup the normalization and the colormap\n",
    "normalize = mcolors.Normalize(vmin=np.log10(tau_vec.min()/2), vmax=np.log10(tau_vec.max()/2))\n",
    "colormap = cm.jet\n",
    "\n",
    "#colors = plt.cm.magma(np.linspace(0,1,2*nb_tau))\n",
    "\n",
    "for t in np.flip(np.arange(nb_tau)):\n",
    "    ax.plot(m_vec/k, mean_snr[:,t], color=colormap(normalize(np.log10(tau_vec[t]))))\n",
    "\n",
    "# Displaying a grid\n",
    "ax.grid(linestyle=':')\n",
    "    \n",
    "# setup the colorbar\n",
    "scalarmappaple = cm.ScalarMappable(norm=normalize, cmap=colormap)\n",
    "scalarmappaple.set_array(20)\n",
    "cbar = plt.colorbar(scalarmappaple)\n",
    "cbar.ax.invert_yaxis() \n",
    "cbar.set_label(r'$\\log_{_{10}}\\, \\frac{2}{\\tau}$', fontsize=18)\n",
    "cbar.ax.set_yticklabels(['3.0','','','','2.0','','','','1.0'])\n",
    "    \n",
    "#plt.legend(('$\\pi$/1000','$\\pi$/100','$\\pi$/10'))\n",
    "#plt.title('Basis Pursuit Recovery Rates')\n",
    "ax.set_xlabel(r'$m/s$', fontsize=18)\n",
    "ax.set_ylabel('Mean SNR (dB)', fontsize=18)\n",
    "\n",
    "\n",
    "#fig.colorbar\n",
    " \n",
    "plt.savefig('BPDN_mean_snr.pdf', bbox_inches='tight', pad_inches=0, transparent=True)"
   ]
  }
 ],
 "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.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": false,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": true,
   "report_style_numbering": false,
   "user_envs_cfg": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}