Neural Network/Neural Network Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural Network Method for Black Hole Imaging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries and modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: ehtim in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (1.1.1)\n",
      "Requirement already satisfied: scikit-image in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (0.18.1)\n",
      "Requirement already satisfied: matplotlib in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (3.3.4)\n",
      "Requirement already satisfied: astropy in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (4.2.1)\n",
      "Requirement already satisfied: future in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (0.18.2)\n",
      "Requirement already satisfied: h5py in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (2.10.0)\n",
      "Requirement already satisfied: networkx in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (2.5)\n",
      "Requirement already satisfied: scipy in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (1.6.2)\n",
      "Requirement already satisfied: ephem in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (3.7.7.1)\n",
      "Requirement already satisfied: numpy in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from ehtim) (1.19.0)\n",
      "Requirement already satisfied: pyerfa in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from astropy->ehtim) (1.7.2)\n",
      "Requirement already satisfied: six in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from h5py->ehtim) (1.15.0)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from matplotlib->ehtim) (8.2.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from matplotlib->ehtim) (2.4.7)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from matplotlib->ehtim) (2.8.1)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from matplotlib->ehtim) (1.3.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from matplotlib->ehtim) (0.10.0)\n",
      "Requirement already satisfied: decorator>=4.3.0 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from networkx->ehtim) (5.0.6)\n",
      "Requirement already satisfied: imageio>=2.3.0 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from scikit-image->ehtim) (2.9.0)\n",
      "Requirement already satisfied: tifffile>=2019.7.26 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from scikit-image->ehtim) (2021.3.31)\n",
      "Requirement already satisfied: PyWavelets>=1.1.1 in /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages (from scikit-image->ehtim) (1.1.1)\n",
      "Welcome to eht-imaging! v  1.1.1\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import csv\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "from numpy.random import randint\n",
    "from scipy.ndimage import gaussian_filter\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "\n",
    "from keras.callbacks import ModelCheckpoint\n",
    "from keras.initializers import RandomUniform, Constant\n",
    "import keras.models\n",
    "from keras.models import Sequential\n",
    "import keras.layers\n",
    "from keras.layers import Layer, Activation, LeakyReLU\n",
    "from keras.layers import Input, InputLayer, AveragePooling2D, Conv2D, MaxPooling2D, UpSampling2D\n",
    "from keras.models import Sequential, Model, load_model\n",
    "from keras.layers.normalization import BatchNormalization\n",
    "from keras.layers.merge import Concatenate, Add\n",
    "from keras.layers import Dense, Lambda, Reshape\n",
    "import keras.initializers\n",
    "import keras.regularizers\n",
    "import keras.callbacks\n",
    "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
    "from keras import backend as K\n",
    "from keras import losses\n",
    "from keras.datasets import fashion_mnist\n",
    "from keras.datasets import mnist\n",
    "\n",
    "# Import eht imaging package\n",
    "!pip3 install ehtim\n",
    "import ehtim as eh \n",
    "import ehtim.const_def as ehc\n",
    "import ehtim.observing.obs_helpers as obsh\n",
    "from ehtim.observing.obs_helpers import *\n",
    "\n",
    "# Import helpers from cosense\n",
    "import helpers_posci as hp\n",
    "\n",
    "# Import utilities for computing data terms, losses, and gradients\n",
    "import helpers_posci as hp\n",
    "from models_posci import IsingVisNet\n",
    "from layers_posci import _unet_from_tensor, Ising_sampling, Ising_sampling2, Site_mask_prob\n",
    "from layers_posci import realnvp_encoder, realnvp_decoder, Lambda_Gaussian, Lambda_logdet_sigmoid, STE_layer, Gaussian_sampling\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Imaging Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Define observation parameters.'''\n",
    "eht_array='EHT2019'\n",
    "target='sgrA'\n",
    "\n",
    "nsamp = 10000\n",
    "npix = 32 \n",
    "fov_param = 100.0\n",
    "flux_label = 1\n",
    "sefd_param = 1\n",
    "\n",
    "tint_sec = 5    # integration time in seconds\n",
    "tadv_sec = 600  # advance time between scans\n",
    "tstart_hr = 0   # GMST time of the start of the observation\n",
    "tstop_hr = 24   # GMST time of the end\n",
    "bw_hz = 4e9     # bandwidth in Hz\n",
    "\n",
    "stabilize_scan_phase = False # if true then add a single phase error for each scan to act similar to adhoc phasing\n",
    "stabilize_scan_amp = False # if true then add a single gain error at each scan\n",
    "jones = False # apply jones matrix for including noise in the measurements (including leakage)\n",
    "inv_jones = False # no not invert the jones matrix\n",
    "frcal = True # True if you do not include effects of field rotation\n",
    "dcal = True # True if you do not include the effects of leakage\n",
    "dterm_offset = 0 # a random offset of the D terms is given at each site with this standard deviation away from 1\n",
    "dtermp = 0\n",
    "\n",
    "array = '/Users/Johanna/Desktop/Proximal Gradient Descent/arrays/' + eht_array + '.txt'\n",
    "eht = eh.array.load_txt(array)\n",
    "\n",
    "# Define observation field of view\n",
    "fov = fov_param * eh.RADPERUAS\n",
    "\n",
    "# define scientific target\n",
    "if target == 'm87':\n",
    "    ra = 12.513728717168174\n",
    "    dec = 12.39112323919932\n",
    "elif target == 'sgrA':\n",
    "    ra = 19.414182210498385\n",
    "    dec = -29.24170032236311\n",
    "\n",
    "rf = 230e9\n",
    "mjd = 57853 # day of observation\n",
    "fwhm = 1.117609542559987e-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Helper Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' \n",
    "    Prepare and return dataset. Fashion MNIST and Digits MNIST images are blurred by 0.1*fwhm.\n",
    "    \n",
    "    -----------------------------------------------------------------------------------------------------\n",
    "    Parameters:\n",
    "        -dataset: 'fashion' (fashion MNIST, 'mnist' (MNIST digits), or 'bh_data' (simulated black hole images)\n",
    "        -flux: sum of pixels per image\n",
    "    ----------------------------------------------------------------------------------------------------------\n",
    "    \n",
    "'''\n",
    "def get_data(dataset='fashion', flux=1):\n",
    "    xdata = []\n",
    "    pad_width = 2\n",
    "    if (dataset == 'fashion' or dataset == 'all'):\n",
    "        npix = 32\n",
    "        (x_train, y_train), (x_test, y_test) = fashion_mnist.load_data()\n",
    "\n",
    "        xdata = 1.0*x_train[[k%60000 for k in range(int(nsamp))]]\n",
    "        xdata = np.pad(xdata, ((0,0), (pad_width,pad_width), (pad_width,pad_width)), 'constant')\n",
    "        xdata = xdata[..., np.newaxis]/255\n",
    "        \n",
    "        xdata = xdata.reshape((-1, npix*npix))\n",
    "        \n",
    "        # Blur images by 0.1*fwhm\n",
    "        xdata_blur = []\n",
    "        for X in xdata:\n",
    "            pred_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "            pred_blur.imvec = X.flatten()\n",
    "            pred_blur = pred_blur.blur_circ(fwhm_i=0.1*fwhm, fwhm_pol=0.1*fwhm)\n",
    "            xdata_blur.append(pred_blur.imvec)\n",
    "        xdata = xdata_blur\n",
    "            \n",
    "    if (dataset == 'mnist' or dataset == 'all'):\n",
    "        npix = 32\n",
    "        (x_train_mnist, y_train_mnist), (x_test_mnist, y_test_mnist) = mnist.load_data()\n",
    "        \n",
    "        xdata_train = 1.0*x_train_mnist[[k%60000 for k in range(int(nsamp))]]\n",
    "        xdata_train = np.pad(xdata_train, ((0,0), (pad_width,pad_width), (pad_width,pad_width)), 'constant')  # get to 160x160\n",
    "        xdata_train = xdata_train[..., np.newaxis]/255\n",
    "        \n",
    "        xdata = xdata_train.reshape((-1, npix*npix))\n",
    "        \n",
    "        # Blur images by 0.1*fwhm\n",
    "        xdata_blur = []\n",
    "        for X in xdata:\n",
    "            pred_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "            pred_blur.imvec = X.flatten()\n",
    "            pred_blur = pred_blur.blur_circ(fwhm_i=0.1*fwhm, fwhm_pol=0.1*fwhm)\n",
    "            xdata_blur.append(pred_blur.imvec)\n",
    "        xdata = xdata_blur\n",
    "        \n",
    "    if (dataset == 'bh_data'):\n",
    "        # Download at https://www.dropbox.com/s/kv0x5eolg10w52g/bh_sim_data.npy?dl=0\n",
    "        bh_sim_data = np.load('bh_sim_data.npy', allow_pickle=True).item()\n",
    "        bh_data = bh_sim_data['image']\n",
    "        \n",
    "        # resize images to 32 x 32 and fov = 100\n",
    "        bh_data = np.array(bh_data)\n",
    "        bh_data_reshape = []\n",
    "        for i in range(len(bh_data)):\n",
    "            bh_img = eh.image.make_empty(160, 160, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "            bh_img.imvec = bh_data[i].flatten()\n",
    "            bh_img_reshape = bh_img.regrid_image(100, 32)\n",
    "            bh_data_reshape.append(bh_img_reshape.imvec)\n",
    "        xdata = np.array(bh_data_reshape).reshape((-1, 32*32))\n",
    "    \n",
    "    return xdata\n",
    "\n",
    "'''\n",
    "    Generate measurement data for image X using specified parameters. \n",
    "    \n",
    "    -----------------------------------------------------------------------------------------------------\n",
    "    Parameters:\n",
    "        -X: target image\n",
    "        -th_noise: If True, include thermal noise in measurements\n",
    "        -amp_error: If True, include amplitude error in measurements\n",
    "        -phase_error: If True, include phase error in measurements\n",
    "        -gainp: Amount of site-wise standard deviation in gain error to include in measurements\n",
    "        -gain_offset: Amount of fixed gain error to include in measurements\n",
    "    -----------------------------------------------------------------------------------------------------\n",
    "'''\n",
    "def get_measurements(X, th_noise=False, amp_error=False, phase_error=False, gainp=0.1, gain_offset=0.1):\n",
    "    # Define noise parameters\n",
    "    add_th_noise = th_noise # False if you *don't* want to add thermal error. If there are no sefds in obs_orig it will use the sigma for each data point\n",
    "    phasecal = not phase_error # True if you don't want to add atmospheric phase error. if False then it adds random phases to simulate atmosphere\n",
    "    ampcal = not amp_error # True if you don't want to add atmospheric amplitude error. if False then add random gain errors \n",
    "\n",
    "    simim = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "    simim.imvec = X\n",
    "    \n",
    "    # generate the discrete Fourier transform matrices for complex visibilities\n",
    "    obs = simim.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, add_th_noise=add_th_noise, ampcal=ampcal, phasecal=phasecal, \n",
    "                    stabilize_scan_phase=stabilize_scan_phase, stabilize_scan_amp=stabilize_scan_amp,\n",
    "                    jones=jones,inv_jones=inv_jones,dcal=dcal, frcal=frcal, dterm_offset=dterm_offset, \n",
    "                    gainp=gainp, gain_offset=gain_offset)\n",
    "    obs_data = obs.unpack(['u', 'v', 'vis', 'sigma'])\n",
    "    \n",
    "    uv = np.hstack((obs_data['u'].reshape(-1,1), obs_data['v'].reshape(-1,1)))\n",
    "    \n",
    "    # Extract forward model (Discrete Fourier Transform matrix)\n",
    "    F_vis = ftmatrix(simim.psize, simim.xdim, simim.ydim, uv, pulse=simim.pulse)\n",
    "    vis = obs_data['vis']\n",
    "    sigma = obs_data['sigma']\n",
    "    \n",
    "    if th_noise:\n",
    "        sigma = 224.46 * np.concatenate([np.expand_dims(obs.data['sigma'], -1), np.expand_dims(obs.data['sigma'], -1)], -1)\n",
    "    else:\n",
    "        sigma = None\n",
    "        \n",
    "    t1 = obs.data['t1']\n",
    "    t2 = obs.data['t2']\n",
    "    \n",
    "    print(\"Finished computing visibilities...\")\n",
    "    \n",
    "    # generate the discrete Fourier transform matrices for closure phases\n",
    "    obs.add_cphase(count='max')\n",
    "    # Extract forward models for telescopes 1, 2, and 3\n",
    "    tc1 = obs.cphase['t1']\n",
    "    tc2 = obs.cphase['t2']\n",
    "    tc3 = obs.cphase['t3']\n",
    "    \n",
    "    cphase = obs.cphase['cphase']\n",
    "    sigma_cphase = obs.cphase['sigmacp']\n",
    "    cphase_map = np.zeros((len(obs.cphase['time']), 3))\n",
    "\n",
    "    zero_symbol = 10000\n",
    "    for k1 in range(cphase_map.shape[0]):\n",
    "        for k2 in list(np.where(obs.data['time']==obs.cphase['time'][k1])[0]):\n",
    "            if obs.data['t1'][k2] == obs.cphase['t1'][k1] and obs.data['t2'][k2] == obs.cphase['t2'][k1]:\n",
    "                cphase_map[k1, 0] = k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 0] = zero_symbol\n",
    "            elif obs.data['t2'][k2] == obs.cphase['t1'][k1] and obs.data['t1'][k2] == obs.cphase['t2'][k1]:\n",
    "                cphase_map[k1, 0] = -k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 0] = -zero_symbol\n",
    "            elif obs.data['t1'][k2] == obs.cphase['t2'][k1] and obs.data['t2'][k2] == obs.cphase['t3'][k1]:\n",
    "                cphase_map[k1, 1] = k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 1] = zero_symbol\n",
    "            elif obs.data['t2'][k2] == obs.cphase['t2'][k1] and obs.data['t1'][k2] == obs.cphase['t3'][k1]:\n",
    "                cphase_map[k1, 1] = -k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 1] = -zero_symbol\n",
    "            elif obs.data['t1'][k2] == obs.cphase['t3'][k1] and obs.data['t2'][k2] == obs.cphase['t1'][k1]:\n",
    "                cphase_map[k1, 2] = k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 2] = zero_symbol\n",
    "            elif obs.data['t2'][k2] == obs.cphase['t3'][k1] and obs.data['t1'][k2] == obs.cphase['t1'][k1]:\n",
    "                cphase_map[k1, 2] = -k2\n",
    "                if k2 == 0:\n",
    "                    cphase_map[k1, 2] = -zero_symbol\n",
    "\n",
    "    F_cphase = np.zeros((cphase_map.shape[0], npix*npix, 3), dtype=np.complex64)\n",
    "    cphase_proj = np.zeros((cphase_map.shape[0], F_vis.shape[0]), dtype=np.float32)\n",
    "    for k in range(cphase_map.shape[0]):\n",
    "        for j in range(cphase_map.shape[1]):\n",
    "            if cphase_map[k][j] > 0:\n",
    "                if int(cphase_map[k][j]) == zero_symbol:\n",
    "                    cphase_map[k][j] = 0\n",
    "                F_cphase[k, :, j] = F_vis[int(cphase_map[k][j]), :]\n",
    "                cphase_proj[k, int(cphase_map[k][j])] = 1\n",
    "            else:\n",
    "                if np.abs(int(cphase_map[k][j])) == zero_symbol:\n",
    "                    cphase_map[k][j] = 0\n",
    "                F_cphase[k, :, j] = np.conj(F_vis[int(-cphase_map[k][j]), :])\n",
    "                cphase_proj[k, int(-cphase_map[k][j])] = -1\n",
    "                \n",
    "    clamparr = obs.c_amplitudes(mode='all', count='max',\n",
    "                                        vtype='vis', ctype='camp', debias=True, snrcut=0.0)\n",
    "\n",
    "    uv1 = np.hstack((clamparr['u1'].reshape(-1, 1), clamparr['v1'].reshape(-1, 1)))\n",
    "    uv2 = np.hstack((clamparr['u2'].reshape(-1, 1), clamparr['v2'].reshape(-1, 1)))\n",
    "    uv3 = np.hstack((clamparr['u3'].reshape(-1, 1), clamparr['v3'].reshape(-1, 1)))\n",
    "    uv4 = np.hstack((clamparr['u4'].reshape(-1, 1), clamparr['v4'].reshape(-1, 1)))\n",
    "    camp = clamparr['camp']\n",
    "    sigma_camp = clamparr['sigmaca']\n",
    "    \n",
    "    mask = []\n",
    "    # shape: (4, 2022, npix**2)\n",
    "    F_camp = (ftmatrix(simim.psize, simim.xdim, simim.ydim, uv1, pulse=simim.pulse, mask=mask),\n",
    "          ftmatrix(simim.psize, simim.xdim, simim.ydim, uv2, pulse=simim.pulse, mask=mask),\n",
    "          ftmatrix(simim.psize, simim.xdim, simim.ydim, uv3, pulse=simim.pulse, mask=mask),\n",
    "          ftmatrix(simim.psize, simim.xdim, simim.ydim, uv4, pulse=simim.pulse, mask=mask)\n",
    "          )\n",
    "    \n",
    "    return obs, vis, F_vis, sigma, t1, t2\n",
    "\n",
    "# Post process predicted image by normalizing flux to that of the target image.\n",
    "def post_process(Z, img):\n",
    "    flux, X_max = np.sum(img), np.max(img)\n",
    "\n",
    "    # Normalize flux to target image\n",
    "    Z_flux = np.sum(np.abs(Z))\n",
    "    Z = (flux/Z_flux)*np.abs(Z)\n",
    "\n",
    "    # Normalize Z between 0 and X_max\n",
    "    Z = np.maximum(np.zeros(np.shape(Z)), Z)\n",
    "    Z = np.minimum(X_max*np.ones(np.shape(Z)), Z)\n",
    "    \n",
    "    return Z\n",
    "\n",
    "# Compute chi-squared loss between true and predicted image\n",
    "def chisq_loss(x_true, pred_vis):\n",
    "    print('HERE')\n",
    "    # compute true visibility\n",
    "    true_vis = keras.layers.Lambda(hp.Lambda_dft(global_F))(x_true)\n",
    "    \n",
    "    # compute chisq loss\n",
    "    num = tf.reduce_sum(tf.square(tf.divide(tf.abs(tf.subtract(pred_vis, true_vis)), global_S)))\n",
    "    chisq = tf.divide(num, tf.cast(tf.multiply(2, tf.size(true_vis[0])), tf.float32))\n",
    "    \n",
    "    return chisq\n",
    "\n",
    "# Compute normalized cross-correlation between images X and Z\n",
    "def compute_xcorr(X, Z):\n",
    "    target_img = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "    target_img.imvec = X.flatten()\n",
    "    recon_img = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "    recon_img.imvec = Z.flatten()\n",
    "    xc = target_img.compare_images(recon_img)[0][0] \n",
    "    return xc\n",
    "    \n",
    "# Compute MAE loss \n",
    "from sklearn.metrics import mean_absolute_error\n",
    "def compute_mae(X, Z):\n",
    "    mae = mean_absolute_error(X.flatten(), Z.flatten())\n",
    "    return mae\n",
    "\n",
    "''' Show results '''\n",
    "def visualize(target, pred, obs):\n",
    "    max_color = 0.01\n",
    "    if dataset == 'fashion':\n",
    "        max_color = 0.005\n",
    "    \n",
    "    # Compute nominally blurred target image\n",
    "    fwhm = obs.res()\n",
    "    target_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "    target_blur.imvec = target.flatten()\n",
    "    target_blur1 = target_blur.blur_circ(fwhm_i=0.3*fwhm, fwhm_pol=0.3*fwhm)\n",
    "    target_blur2 = target_blur.blur_circ(fwhm_i=0.7*fwhm, fwhm_pol=0.7*fwhm)\n",
    "    \n",
    "    # Show nominally blurred target image\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(figsize=(13, 3), ncols=3)\n",
    "\n",
    "    ground_truth = ax1.imshow((target).reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax1.title.set_text('Ground Truth')\n",
    "    fig.colorbar(ground_truth, ax=ax1)\n",
    "    blur1 = ax2.imshow(target_blur1.imvec.reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax2.title.set_text('0.3 * fwhm Blurred Truth')\n",
    "    fig.colorbar(blur1, ax=ax2)\n",
    "    blur2 = ax3.imshow(target_blur2.imvec.reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax3.title.set_text('0.7 * fwhm Blurred Truth')\n",
    "    cbar = fig.colorbar(blur2, ax=ax3)\n",
    "    cbar.minorticks_on()\n",
    "    plt.show()\n",
    "    \n",
    "    # Compute nominally blurred predicted image\n",
    "    fwhm = obs.res()\n",
    "    pred_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "    pred_blur.imvec = pred.flatten()\n",
    "    pred_blur1 = pred_blur.blur_circ(fwhm_i=0.3*fwhm, fwhm_pol=0.3*fwhm)\n",
    "    pred_blur2 = pred_blur.blur_circ(fwhm_i=0.7*fwhm, fwhm_pol=0.7*fwhm)\n",
    "    \n",
    "    # Visualize blurred results\n",
    "    fig3, (ax7, ax8, ax9) = plt.subplots(figsize=(13, 3), ncols=3)\n",
    "\n",
    "    pred_img = ax7.imshow(pred.reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax7.title.set_text('Predicted Image')\n",
    "    fig3.colorbar(pred_img, ax=ax7)\n",
    "    blur1_img = ax8.imshow(pred_blur1.imvec.reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax8.title.set_text('0.3 * fwhm blur')\n",
    "    fig3.colorbar(blur1_img, ax=ax8)\n",
    "    blur2_img = ax9.imshow(pred_blur2.imvec.reshape(npix, npix), vmin=0, vmax=max_color)\n",
    "    ax9.title.set_text('0.7 * fwhm blur')\n",
    "    fig3.colorbar(blur2_img, ax=ax9)\n",
    "    cbar.minorticks_on()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading dataset...\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "# Select dataset: 'fashion' for Fashion MNIST, 'digits' for MNIST digits, or 'bh_data' for black hole simulated images\n",
    "print('Loading dataset...')\n",
    "dataset = 'fashion'\n",
    "xdata = get_data(dataset)\n",
    "print('Done.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select target image and generate observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating empty observation file . . . \n",
      "Producing clean visibilities from image with nfft FT . . . \n",
      "Adding gain + phase errors to data and applying a priori calibration . . . \n",
      "Finished computing visibilities...\n",
      "Getting bispectra:: type vis, count max, scan 49/106 "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Johanna/opt/anaconda3/lib/python3.7/site-packages/ehtim/obsdata.py:486: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray\n",
      "  return np.array(datalist)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Getting bispectra:: type vis, count max, scan 106/106 \n",
      "\n",
      "Updated self.cphase: no averaging\n",
      "updated self.cphase: avg_time 0.000000 s\n",
      "\n",
      "Getting closure amps:: type vis camp , count max, scan 106/106\n",
      "\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "# Select target image at index 0\n",
    "target = xdata[0]\n",
    "\n",
    "# Normalize target images\n",
    "target /= np.sum(target)\n",
    "\n",
    "obs, visibility, F_vis, sigma, t1, t2 = get_measurements(target.flatten())  # data terms\n",
    "fwhm = obs.res()\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Pre-Trained Model\n",
    "\n",
    "Pre-trained model can be downloaded at https://www.dropbox.com/s/fszl68xkdr6hdt6/pretrained_model?dl=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading pre-trainend model...\n",
      "WARNING:tensorflow:From /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "If using Keras pass *_constraint arguments to layers.\n",
      "WARNING:tensorflow:From /Users/Johanna/opt/anaconda3/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:4074: The name tf.nn.avg_pool is deprecated. Please use tf.nn.avg_pool2d instead.\n",
      "\n",
      "Replacing input layer...\n",
      "Copying shared layer weights to vis_model...\n",
      "Model: \"model_2\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "vis (InputLayer)                (None, 1691, 2)      0                                            \n",
      "__________________________________________________________________________________________________\n",
      "reshape_3 (Reshape)             (None, 3382)         0           vis[0][0]                        \n",
      "__________________________________________________________________________________________________\n",
      "dense_dirtyimage (Dense)        (None, 1024)         3464192     reshape_3[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "reshape_4 (Reshape)             (None, 32, 32, 1)    0           dense_dirtyimage[0][0]           \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_20 (Conv2D)              (None, 32, 32, 64)   640         reshape_4[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_19 (LeakyReLU)      (None, 32, 32, 64)   0           conv2d_20[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_19 (BatchNo (None, 32, 32, 64)   256         leaky_re_lu_19[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_21 (Conv2D)              (None, 32, 32, 64)   36928       batch_normalization_19[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_20 (LeakyReLU)      (None, 32, 32, 64)   0           conv2d_21[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_20 (BatchNo (None, 32, 32, 64)   256         leaky_re_lu_20[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling2d_5 (AveragePoo (None, 16, 16, 64)   0           batch_normalization_20[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_22 (Conv2D)              (None, 16, 16, 128)  73856       average_pooling2d_5[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_21 (LeakyReLU)      (None, 16, 16, 128)  0           conv2d_22[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_21 (BatchNo (None, 16, 16, 128)  512         leaky_re_lu_21[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_23 (Conv2D)              (None, 16, 16, 128)  147584      batch_normalization_21[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_22 (LeakyReLU)      (None, 16, 16, 128)  0           conv2d_23[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_22 (BatchNo (None, 16, 16, 128)  512         leaky_re_lu_22[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling2d_6 (AveragePoo (None, 8, 8, 128)    0           batch_normalization_22[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_24 (Conv2D)              (None, 8, 8, 256)    295168      average_pooling2d_6[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_23 (LeakyReLU)      (None, 8, 8, 256)    0           conv2d_24[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_23 (BatchNo (None, 8, 8, 256)    1024        leaky_re_lu_23[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_25 (Conv2D)              (None, 8, 8, 256)    590080      batch_normalization_23[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_24 (LeakyReLU)      (None, 8, 8, 256)    0           conv2d_25[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_24 (BatchNo (None, 8, 8, 256)    1024        leaky_re_lu_24[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling2d_7 (AveragePoo (None, 4, 4, 256)    0           batch_normalization_24[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_26 (Conv2D)              (None, 4, 4, 512)    1180160     average_pooling2d_7[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_25 (LeakyReLU)      (None, 4, 4, 512)    0           conv2d_26[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_25 (BatchNo (None, 4, 4, 512)    2048        leaky_re_lu_25[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_27 (Conv2D)              (None, 4, 4, 512)    2359808     batch_normalization_25[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_26 (LeakyReLU)      (None, 4, 4, 512)    0           conv2d_27[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_26 (BatchNo (None, 4, 4, 512)    2048        leaky_re_lu_26[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "average_pooling2d_8 (AveragePoo (None, 2, 2, 512)    0           batch_normalization_26[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_28 (Conv2D)              (None, 2, 2, 1024)   4719616     average_pooling2d_8[0][0]        \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_27 (LeakyReLU)      (None, 2, 2, 1024)   0           conv2d_28[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_27 (BatchNo (None, 2, 2, 1024)   4096        leaky_re_lu_27[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_29 (Conv2D)              (None, 2, 2, 1024)   9438208     batch_normalization_27[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_28 (LeakyReLU)      (None, 2, 2, 1024)   0           conv2d_29[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_28 (BatchNo (None, 2, 2, 1024)   4096        leaky_re_lu_28[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "up_sampling2d_5 (UpSampling2D)  (None, 4, 4, 1024)   0           batch_normalization_28[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_5 (Concatenate)     (None, 4, 4, 1536)   0           batch_normalization_26[0][0]     \n",
      "                                                                 up_sampling2d_5[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_30 (Conv2D)              (None, 4, 4, 512)    7078400     concatenate_5[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_29 (LeakyReLU)      (None, 4, 4, 512)    0           conv2d_30[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_29 (BatchNo (None, 4, 4, 512)    2048        leaky_re_lu_29[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_31 (Conv2D)              (None, 4, 4, 512)    2359808     batch_normalization_29[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_30 (LeakyReLU)      (None, 4, 4, 512)    0           conv2d_31[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_30 (BatchNo (None, 4, 4, 512)    2048        leaky_re_lu_30[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "up_sampling2d_6 (UpSampling2D)  (None, 8, 8, 512)    0           batch_normalization_30[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_6 (Concatenate)     (None, 8, 8, 768)    0           batch_normalization_24[0][0]     \n",
      "                                                                 up_sampling2d_6[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_32 (Conv2D)              (None, 8, 8, 256)    1769728     concatenate_6[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_31 (LeakyReLU)      (None, 8, 8, 256)    0           conv2d_32[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_31 (BatchNo (None, 8, 8, 256)    1024        leaky_re_lu_31[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_33 (Conv2D)              (None, 8, 8, 256)    590080      batch_normalization_31[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_32 (LeakyReLU)      (None, 8, 8, 256)    0           conv2d_33[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_32 (BatchNo (None, 8, 8, 256)    1024        leaky_re_lu_32[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "up_sampling2d_7 (UpSampling2D)  (None, 16, 16, 256)  0           batch_normalization_32[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_7 (Concatenate)     (None, 16, 16, 384)  0           batch_normalization_22[0][0]     \n",
      "                                                                 up_sampling2d_7[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_34 (Conv2D)              (None, 16, 16, 128)  442496      concatenate_7[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_33 (LeakyReLU)      (None, 16, 16, 128)  0           conv2d_34[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_33 (BatchNo (None, 16, 16, 128)  512         leaky_re_lu_33[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_35 (Conv2D)              (None, 16, 16, 128)  147584      batch_normalization_33[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_34 (LeakyReLU)      (None, 16, 16, 128)  0           conv2d_35[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_34 (BatchNo (None, 16, 16, 128)  512         leaky_re_lu_34[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "up_sampling2d_8 (UpSampling2D)  (None, 32, 32, 128)  0           batch_normalization_34[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_8 (Concatenate)     (None, 32, 32, 192)  0           batch_normalization_20[0][0]     \n",
      "                                                                 up_sampling2d_8[0][0]            \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_36 (Conv2D)              (None, 32, 32, 64)   110656      concatenate_8[0][0]              \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_35 (LeakyReLU)      (None, 32, 32, 64)   0           conv2d_36[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_35 (BatchNo (None, 32, 32, 64)   256         leaky_re_lu_35[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_37 (Conv2D)              (None, 32, 32, 64)   36928       batch_normalization_35[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "leaky_re_lu_36 (LeakyReLU)      (None, 32, 32, 64)   0           conv2d_37[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "batch_normalization_36 (BatchNo (None, 32, 32, 64)   256         leaky_re_lu_36[0][0]             \n",
      "__________________________________________________________________________________________________\n",
      "conv2d_38 (Conv2D)              (None, 32, 32, 1)    65          batch_normalization_36[0][0]     \n",
      "__________________________________________________________________________________________________\n",
      "recon (ReLU)                    (None, 32, 32, 1)    0           conv2d_38[0][0]                  \n",
      "__________________________________________________________________________________________________\n",
      "pred_vis (Lambda)               (None, 1691, 2)      0           recon[0][0]                      \n",
      "==================================================================================================\n",
      "Total params: 34,865,537\n",
      "Trainable params: 34,853,761\n",
      "Non-trainable params: 11,776\n",
      "__________________________________________________________________________________________________\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "print(\"Loading pre-trainend model...\")\n",
    "model = IsingVisNet(t1, t2, F_vis, n_ising_layers=3, slope_const=3, obs_prob=1)\n",
    "model.load_weights('pretrained_model')\n",
    "\n",
    "# Remove input image layer, replace with direct complex visibility input\n",
    "print(\"Replacing input layer...\")\n",
    "vis_input = keras.layers.Input(shape=(1691,2), name=\"vis\")\n",
    "\n",
    "filt = 64\n",
    "kern = 3\n",
    "acti = None\n",
    "F = F_vis\n",
    "\n",
    "vis_reshape = keras.layers.Reshape((2*F.shape[0], ))(vis_input)\n",
    "dirty_im = keras.layers.Dense(32*32*1, activation=acti, use_bias=True, kernel_initializer=RandomUniform(minval=-5e-4, maxval=5e-4, seed=None), name='dense_dirtyimage')(vis_reshape)\n",
    "dirty_im_reshape = keras.layers.Reshape((32, 32, 1))(dirty_im)\n",
    "pred_img = _unet_from_tensor(dirty_im_reshape, filt, kern, acti)\n",
    "pred_img = keras.layers.ReLU(name='recon')(pred_img)\n",
    "pred_vis = keras.layers.Lambda(hp.Lambda_dft(F), name='pred_vis')(pred_img)\n",
    "\n",
    "# Model returns predicted image and predicted visibilities\n",
    "vis_model = keras.models.Model(inputs=vis_input, outputs=[pred_img, pred_vis])\n",
    "\n",
    "print(\"Copying shared layer weights to vis_model...\")\n",
    "for i in range(1, len(vis_model.layers)-1):\n",
    "    vis_model.layers[i].set_weights(model.layers[i+1].get_weights())\n",
    "    \n",
    "opt = keras.optimizers.Adam(learning_rate=0.001)\n",
    "vis_model.compile(optimizer=opt, loss='mae')\n",
    "\n",
    "print(vis_model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Noise Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "npix = 32\n",
    "flux = 1            # image flux\n",
    "th_noise = False    # add thermal noise\n",
    "phase_err = False   # add phase noise\n",
    "amp_err = False     # add amplitude noise\n",
    "gainp=0.0           # std of site-specific gain error \n",
    "gain_offset=0.1     # fixed offset of site-specific gain error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Prediction with Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwMAAADWCAYAAACNOV4iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABD/0lEQVR4nO29e7hsV1nm+/vW2pe19wZ3THYSFJEIhksIfYIShG5bIjbYKAhe6OamHQOdg3hAjx3UhhZ3sLWx2xx57IeYs7l0bDuCRwUhSjSA2ZCowSQGhMR0NBIMbYckJDsgyd57Xb7zx5wrVKreWWvUmrVWrTXn+3ueetaqUWOOOWZVvWPWmGO+3xeZiTHGGGOMMaZ/zM26A8YYY4wxxpjZ4MmAMcYYY4wxPcWTAWOMMcYYY3qKJwPGGGOMMcb0FE8GjDHGGGOM6SmeDBhjjDHGGNNTej0ZiIidEfGTEXFNRFxV/70uIt4eEc+Zdf/WIiLOjYhPRkRjfNiIuLSuc7h+3BkRRwaeXxsRh6fUn3Mi4tyhsosi4vZp7cNsbSLi6bWOPh4RN0TEDxds87KI+OOI+EhEXB0R10fEDxXu79yIOFhY9zkR8Rf19/6TEfGCku0Gtv/htfQ2TSLiiiG9Xh0RfxsRnxgcnyLibVtdY3X/jzR9VhFxmhiXUoxd50ypPwcj4rSB50+t95HT2kcfWaf+jw58xquPmyPicwXbWv/W/3r6Y/0Pk5m9fAC7gauA9wMnDJQ/BrgO+Pys+1h4HOdUH2Pj65cC5ww9Pzzw/LTB5y37clC11VTuR7cewKOB+4Dvr58/FjgCfO8a210P/MjA8xcCy8CZDfX3Av8R2AecW3+/dgBvBB49Zj9/B/xo/f8zgeet4xjH6m0D3tNhve4Efg/4EnDKQPmW1xhwGDjY8NppjI5LKcauc6bUl1RtNZX7UfSerlf/t4uytwO/3FDf+rf+2/bF+h969Hll4CBwBvDKzDyyWpiZdwDnzahPG8HvA7ePef0IlciMacvrgbsy8/0Amfk54LeBn19ju9cBvzXw/DDVquXpDfV3APPANVQn9ScCHwe+iWowb+KbgM/Wfbs2M69co19bjsxcBP478Ejg22bcnWlyhLXHod9n/FhmZst69f+8wScRsRt4KfCuhvrWv/VvpkwvJwMRMQ+8Bvj9zPzK8OuZ+WngTXXdk+vlo6MR8e8i4r/Vy3Urq8tJEfGD9ZLoJyLi0xHxHyJirn7twuElvIj4wOCS2dAS1asj4nfr9q6OiMcO9f3lEXFrvYz2fuDx4441M38/M28f8/qRzLw0It4QEbfUfX1FRFweEf9Q92u4v980vKQWEW+iukpz1sCy3llDfX9DRHw4Iv4mIi4Y12+zLXke1araIJ8Azo6IE5s2ysw/z8wlqG7dA94A3Ax8uKH+lzLz31P9YHgR8ALgpzPz32bmPwzXX9VX/fRt9XfzFRHxqfo7fE1EPLWu+9aIeN/AtldExJcj4m1DbX5vrZG/jojfjIhddflzB5a4XxIR76t1dVlE7I6IN0Z1S+Jfx/pvRdxZ/11RL0bEswb6cFpd9mOr+h6od1m99P4bEfELEXFlRHylXkL/aK35t0bEr0XEn0TEYtS3AQ4c559GxJ9FxJsjYsdA2wfqceyzdbtvGHdAq+PQGnV+H/jKmPH4+2PoNo7h8TciHi2+C28b2tU3132/Xo3BppH16v/WoaIXAzeJ8tX61n+F9W/9T49ZL03M4gE8heoKwusn2OZ24K+Ak+vnbwe+HfhO4Dhwdl1+APhb4D8MbHuQoSU8xJJZ3ac/BHbWzz8A/LeB159ONQA8v37+SOBPmWDZkqFlx6HXzgUeAH6yfv4Y4A/W6O85445zoPzLwHfXz/8PqttAvmnW3wU/pvegWrr+laGy59ffk7MLtn871W0GfwacNqbeI4FfAm4ELgHeU2/z68Cjxmw3/H09tdbTcwbKbgb+EVion38D1UWD1dfPqdt5Y/18L/C/gH8zUOe0us7b6+cL9fhxBfCkuuzHgdsK3pOH6RU4AfgT4KbVPtblB9HL7KcNlJ3L0C0ZdftHgG+rn78M+Nn6/8P1sT2ufv4zwMuBJwFHgX828Hl8CviFgXYvBz4G7K6f/0j9vh4s/C6t9v+chtdvR4/H5zA0Hg6/N+q7MFQ+OAa/n4Ex2I+N0//ANn/EwG2D4nXr3/q/Het/qo9ergxQCQpgZFVgDd6fmXcDZOaPZ+Y1wJuBD2fmdXX5PcAh4I1RLXdOym9ntRQIleifNvDa/011xeSKel9fBn5jHfsYxw7g4rr9OzJzIpPVGO7OzD+u2/0UcD/VpMB0h33AsaGyYwOvjSUzfxw4iWoA/ouIOKOh6jLVoP3twLXA/wS+A/g81fe3iMz8AtWVy+8DiIjTgVupTt7fVVf7PuCDYvP/UbfxQN3G00Sd99R1jlJdMd2ZmbfUrx0GHhcRJxR0dXW17Wqq+54PAK+q250Gn8rMT9R9fU9mvnXgtY9k5t/Vr/1yZv4W8LPA9Zn5p3X5l6nej9fDQ+/jC4Bfy8xjdZ3/TqX5aaLG42kwOAZfhf5szSit9A8QEd8APAP43THVrH/rH6z/qdLXycCR+u/DBqiIOLEW3TX10tIzh7b7e9HWPwH+ZqjsVmAPzfc8j+N/Dfz/JWD/wPMzqFYdBvnsOvYxji9k5vEptwkPPy4YPTaz/fkKlTF/kN0Dr61JZq5k5n8B7gTe0lDngcx8Uw7c4peZS5n5i5n5+Qn7/EHqHwNUtxxcSrXatlr2QuAPxHbjdLrK/x74/yviOQ3bDfPJzDwnM/95Zp4I/DRwVUT8nwXblqDGtXGvPRU4PQaivwD/GrgvIr6WapyC0bHq9rYdLejbNCj5bM0orfUP/BuqH2MPNFWw/q3/MX2bBr3Uf/EsumPcQjUhOHOwMDPvBc6p77P7LNUVgkGWG9rLoecxVD78OjS/94P7yIG2NouiYxy8P7BFu5t9bGZjuQ141FDZ1w28NkJEBNUVs+EJ6F9TsHKUa9xnWsAHgV+KiDOp7nm+GHgccEFE/DTVUvxdYr/D32f1XR6uMxUNZOYfRcRVVKuS/29TNVFWMuaUvnZtZr5IvVB9pJuC6tskx13arsepMibWv+BcqltVirD+rf8hrP910suVgVrIvw68OCIe0bK5TwFPGCo7nere+9WZ8ZeArxmq8w3r2NdngG8eKvumdbSzHoaP4TGizkOGpojYFRF7NrxXZitxJZWvZZBnUC0p39uwzWOBvxTlX8/oatLUycybqHT6KuCB+orkB6h+xPwc1f3LW5FFvmokVHyp/ruWZtfDp6juG36I2jD4jvrpTfXf4bHqsVPa/zi+VPdnreN+6EdDRHxNbOIvmA6zHv0/RET8c+ArmXn9RnROYf2vC+u/g/RyMlBzIdWX9j0RcdJqYR0V4Pn1U+nWH+ItwHMj4ux6+wPA+cB/Wr1fjurHzpMj4pS6zoup7o2elF8FzoiI59ftPBJ49TraWQ9/CXx7VJGYQIdf/QKwGjXip6jeY9Mffg04NSJeBBAR30i1fHxwtUJEfE9E/O+IGLzqf8bqNnWdFwD/jMp7sxlcDvxfVMYxMvM2qpWJn0DfLzxTai/Fd1KF2pNk5n1Uq5vn1NucSBWlZRq8FTgtIl5etx1UWl+9f/dv6769ftU3FRGvpDJsbjR/QxWs4Jx6v4/lq/d/D/IF4MSoolf9PYX3tJuxrFf/q/wozeFENxLrfzKs/y6yGS7lrfoAdlH9aL2Wyvl+LdWs97eoE6VQ3ft/mMo9fwvwO6KdHwJuoDISfYbqisLcUJ3/Um//YeCCus3bgf9MdXX/MNVs9ZPAv6RaLr2l3u9hYFfdzkupPAmfoLqX8XX1doeBp65xvB+luhf7SF3/5IHXfmxof987tO3XUl0x+TTVoLkaJeKTwMvrOidS3W/5p/V7+XiqQeL2ep/vresNvp/nz/p74MdUNXU2Vfzvq2tN/MjQ6y+uvwvfWj9foIpScV39vbmm/m6/dIp9euqQvv5g6PVnU038HzVQ9p+Avxmq9/319qt6exTwtlpTd1LdYvCs+ruf9d9vHarzX4Gzhuo8q6HfVwzp9XD9Pv018CvAnrre2wY0duXA9s+hGo+uoTL4XTCg7xOpbjFY7ddhBhI21Vo/Urd7GJgf6ttz6s/phrr9twI7Bl4/kcoE+lngI1QXTT5Wt3fRGp/XDw+8P58E3jDwWsl4/NL6tY9RrQD/snhvzqeKHHMd1Y/V4jHYj+npf6B8H9WPya/dgD5Z/9b/8Htj/Q89on5jjDHGGGOMMT2jz7cJGbPpRMRCRFwaVcKW6yPieWPqviKq5HPXR8RFg/c1rtVOVMlkDg88/uNGHpcxZm2sf2P6y1bWf1+jCRkzKw4CkZnPjIgnANdGxJOzinn9EHV0i4uoIl7dS3WL12upkquUtPNHmXnuhh+NMWYSDmL9G9NXDrJF9e+VAWM2iYiYozJ8vwsgM2+lyqL5SlH9VcCHMvOezFwB3g28Zh3tGGO2ANa/Mf1lq+t/KpOBSZY+jOkxj6OKInXLQNnNjIbjg8qIN1zvKVGFay1p58kRcUVUCfQORRXlakOw/o0pwvo3pr9saf1Pa2XgIPWSBfBy4L0RsRlhpIzZTqxq4shA2f3AKQ11h+sFVSr6knZuogrr98/r1/5o8J7DKXMQ69+YtbD+jekvW1r/rT0DA0sWPwDVkkVErC5ZXKS22RW7c8FhXTeGjRjuOx5w6svcd09mnrxWve/+zn35xXubkzXe8FfHbqIKQ7bKocwsiZU/yac2ru5Dr2XmQ3kgIuIgVQjdZ1CFg5sa1v8MKf7WlFZsKfRtOk5Y/+vH+m+J/DTLvg7l8m/5o6Aw4mS5/AtrbtJ4Yv1XTMNAPMnSBwAL7OPbQuWB6CGlQo2yRZyYm/5sIFdKxVuSo61p29n9kvhI/u7nSurdc+8Sf/ZHj258feHrP3s0Mxu/98BqWvsTBv7fP/D/cN0TBp7vpxoe7wZ2T9AOmfmViLiXKgPkVH8MYP23Q+l/2lovbK+Vftm+44T13wrrf5gJNC01rOqKesULPXMtbwBZKdOrDFOvxoRC/RePJ80NFFX7yMrvWP9M5zahoqWPiDi/vp/w+kWOYcx2I4EVsvFRwG1UkQGeOFB2BlXik2GuE/VuyswH12onIp4TEQ8NSlFl1T4B+IeSTk6I9W96gfUvsf5NL+i6/jcymtDDpq2ZeSgzn56ZT9/50MTGmO1DkizmcuNjze2rqADvAM4DiIjTqbJRXhYRByLi6qjSxgO8E/ieiDipXoo/F7hkrXbqbb8R+LGBewRfB/wd078qOA7r33QK638irH/TKbqu/2ncJjTJ0kd/meS+PbFMWLycqCi9naBhWS7mRaFcghutWL7U1+4Whc26zajwCsA4DgKXRMS1VPp7WWbeGRGPAZ4E7AXuzczPRMQFwJVUb87HqdLdj22nfu1PgO8Arq4HhC8DL8zMxbadF1j/JTTpv1Trhdu2uZ0gc4IxSug65kqX/1vctjTJOLoBY4L1P0K/9b8RnuxSDavbf0rrTYJqU2gr1O1E86PbynGmxXgCDb8zSn8fTUCX9T+NycDgksXqAHAG8KEptG3MliGBxZaTlsw8SjXLHy6/Azh5qOwyvjrbL2qnfu3vqa8abALWv+kF1r/E+je9oOv6bz11KliyMKYTJLCc2fjoI9a/6QvW/yjWv+kLBfrfX8f0f+GMu7ouprEyAOOXLIzpBEmyuF3jJ24sB7H+Tcex/hs5iPVvOk6B/u/PzPM3qz/TZiqTgXFLFsZ0hUxY9G+BEax/0wesf431b/pA1/U/rZUBM0hbs/C8cOwqU5Gq18ZA1BRPWMYPHq2r4gxHlMYZbmM+hlYG5OLdBMsbktXNdIoNiDMe88osLMrUmFBI4zdbjQtdGhOsf1NCi5xAjUEBhF5Lz+ulvxNKAw00BRCRKL2K7XN5NMpOqLEj1NhR3sfyICeCtQMBre6l0/r3ZMCYQhJYnCTiijGmM1j/xvSXruvfkwFjCkno9JUBY0wz1r8x/aXr+vdkwJhCqisDG5mnzxizVbH+jekvXde/JwPGFJIEyxuatNsYs1Wx/o3pLwX63x8Rh4DLM/PyTerW1PBkYJOQZh/QZqOd4mNRRqMdop4yBs0VmgqbDDfK8KOMQYUGolLzoa7XYHIS9/JNZjZcm65fGTDroNAs3NZAKLWuxhRlNG5rIFTjwrLSa2E9pXUxThSbj5sovb93Ap+x9W8eRqnWG87/rbQu6slAA22zJAsdKm3KLMIr4vyvfhOosaf0twM0BC8oPO7CvNwF+ndoUWP6QbDsHwPG9BTr35j+0m39ezJgTCFVOvL1h240xmxfrH9j+kvX9e/JgDGFZAaL2d3BwBjTjPVvTH/puv49GTCmkCq0WHeXCY0xzVj/xvSXruvfk4G2lGYbbcpAKszCsWvnaMWdu0br7SgzFWlToeh3kzmn0ASYS0uju1EmINVeqVlI1WvoT7EBuTADYRIspiVjxjORgVCZAFUAAaV/NU6oMaE0A3kT0vAryoT+9dhRaCosNR83bD99A6H13xvaBAYozSpMg1lY6FrW2yn0X5ipPJsCGghUYBCZRbjFOJGiXrH5GPR5vbSPhRTo39GEjOkLyx3OQGiMGY/1b0x/WUP/jiZkTB/wlUFj+ov1b0x/6br+u3tkxkyZrt8zaIxpxvo3pr90Xf+eDBhTSHVloLvRBIwxzVj/xvSXruvfkwFjCul6aDFjTDPWvzH9pev692RgEkojDKhoAipCCA0RQXbvHq23MFqWMurQ6H5yR7toQrEsypfKnP6h6skyEYmkNDpJw76LIxQd102ObAqdzkBo1qBF5DAZDQSIXSJK0O7RMkSZ1L8oy3nV7wmMsEL/MkqY0HU7/Yt6KsIItI9GVID1b0YoPf836F9GDhLnf6lrGXWsxfm/gZRR/cSYoLRVOCaUjidN0YRkJMPSaIKFdF3/ngwYU0jXlwmNMc1Y/8b0lwL9O7SoMX2hywYiY8x4rH9j+ssa+t/WoUU9shlTyOqVgaZHCRGxEBGXRsS1EXF9RDxvTN1XRMQNdb2LIr66tlvaTkRcGBG3T3qsxpiHY/0b01+6rn+vDBhTSAIr7e8ZPAhEZj4zIp4AXBsRT87MLwxWiogzgYuAM4F7gY8CrwXeXtpORJwCvJriHKvGmCasf2P6S9f178lAS1Q68hApwRsNRDuFgVCZhfcujJStLIxuu7JQZiBSpkJpuAFpFppbEunDF4UxqLCMRWEAUgYiUQ8gj49+34vTmX9FNjm6j5bRBCJijkqcP1C1l7dGxI3AK6mEP8irgA9l5j31tu8Gfhp4+wTtvBm4pG7LTEJpsIBS/TcEEFDGQJT+94yWtdH/RGvC0kAsxoRp619pXY0JQC4K/Zeaigux/jtKm8AgQuuIelLnNJiFVbAApfXdwlS8e3TfKfqYYtxqRPwuCGHEDfWbQOiteExQRuOG838sifG1NIBIIV3Xv28TMqaQhLbLhI8DTgJuGSi7GXi6qHu2qPeUiNhT0k5EPB74euDqko4ZY8Zj/RvTX7quf08GjCkmWM65xgdwoL5/b/UxbCY6tf57ZKDsfuAUsbNTRb0ADhS28wvAz5cfmzFmPNa/Mf2l2/r3bULGFLJ6ZWAM92SmmuWvxQRrtmPrBkBEfCuwmJmfjohz1tEfY8wQ1r8x/aXr+vdkwJhCkmAlJ9HtCHfVf08Y+H//wP/DdU8YeL6fajy6G1i90bSpnbcAP96mo8aYh2P9G9Nfuq5/TwaaaJFtVBqImjIQl5qF943WW3rEqKloeWF038u7R49FGYiaDMQhfDixNFp3/vioWUgajVW948JAeEwYiI7pdMGxQ7znE2QwLCFzzSsDa3EbVWSAJ/JV4Z4BfEjUva6ux0C9mzLzwYhobCciHgl8M3BpHYnsBOBREXEY+O3M/PU2B2CGkFoXZsEJMhBLs7DQ//K+0f0sCf2v7BQGwglG/lL9zy2O6npelLXSvwgUUJVb/1j/a1OYeVcFBkCVFWYbDjUmgM4srMzCe9W5fnQ/K8JAvLyrLIBANrw3OpOvaFIGFVBaV2VK68IsvNgwcJUGIGmTgbzj+vdkwJgJaHNlIDNXIuIdwHnA1RFxOnAW8IqIOAC8H3hRZt4LvBP4cEScBNwHnEsVGWBsO5n5ZQYGkXqZ8NLMPGfdHTfGANa/MX2my/q3gdiYQqaRdIQ6PnBEXAu8B3hZZt4J7AGeBOwFyMzPABcAVwKfAD4JXFzQzkNExJuBt1FfGYiIZ6zjsI0xWP/G9Jmu698rA8YUkgRLK62WCcnMo1Sz/OHyO4CTh8ouAy6bpJ2hOm+hun/QGNMS69+Y/tJ1/XtlwJgJWCEaH8aYbmP9G9Nf1tD//og4FBEvnHU/14NXBiag1FQkMxU2GIhSZBFUmUWVWXjxEaMf3+K+0fnd8i5hIJ5gghvKLCQNhKP7nl8U9Y4J8/Gx0WOZPzpqCpp7UHd87qh4z48Js2FbA2HLKwNm+yKzDQvTnc5A3KB/mUVYmAWFWVjqf2+Z/ld2lP94lfqXZsEZ6l8FEFBZyVtlILb+O0lhYBCtdfF9UMECGjIQS/3vGa27LMqW9o3uZ3lhtI/LOwv13zQkCP9w8W+CJTEmCP3PiQAA88dGy+aONWQgPi4+Q5XVeGP1f39mDucW2DZ4MmBMIVMILWaM2aZY/8b0l67r35MBYwpJYCl9Z50xfcT6N6a/dF3/ngwYMwErHR4MjDHjsf6N6S9d1r8nA8YUkhmdvjJgjGnG+jemv3Rd/54MQKtsw9JApLINN2Ug3K0MRKPbL+0d3c/iI0b7c3yfMBAJU9HKJBlIlYFoabRNaRYUCYPnj4/W23F09Fh2iMyJO4Q5s9qRKJpTpiJtQCohgaWV7g4GvWXa2cYnMBCq8pW9wiy8d7TN448QY8Jeof/R5MWkMBA2neeUWVBlJZ5bHG1zTnj4dwgD4Y6jYowSxsdG/c8Lc6fSf5sAAlj/W5LCrMJVXXEOb5FtWJ3rQ2UVVr8JgNwlsggvCLOwOP8viWABS+Jcv7RbaF10Z5Lb4bWBWOlfBAvYVaZ/NUbNq88FmCsN5LKx+t8fEYeAyzPz8nXvaEZ4MmDMBHTZQGSMGY/1b0x/WUP/jiZkTB9Iur1MaIxpxvo3pr90Xf+eDBhTSvrKoDG9xfo3pr90XP/F05yIODsibouIc4fKz4qIP4+IayLi8og4aeq9NGYLsHrPYNOjy1j/pu9Y/9a/6S9d13/RykBEfD/wEuDIUPku4APAuZl5VURcCFxS192aFJqNWpmKlIGwwUC0IgxEywvCQLxn9MumzIJL0kAs9qsMRE3fZ2EgnhM+nCVlKlYG4mOiP8IsqLMkN3xlU3RSEMJoWErXk4400Xn9F5oKdWZhYSAU+m8yEK4IrS8vFJoF95Tpf0npv8HPXEqslJmFlf5XRmMmsLxz9Jh3Kq02jd9C/3L4ECbnUqz/Tdb/JMZguX2hWVgZTVWZ0nDpub4hgMDK7jL9L6vzv9K/KFsWBmKl/0nugNEGYlEmxoSVYyqgyah+xZDQaHLeoeI9zI92UmVELqXr+i+9Tei6zHx/RBweKn8+sJyZV9XP3wncHhEnZ+bd0+qkMVuF5Q7fMzgG698YrP+hcuvf9Iou67/oyDLz8w0vnQ3cMlDvDuAB4Fvad82YrUUmLK/MNT66ivVvjPUvsP5NbyjQ//6IOBQRL5x1X9dDWwPxqQwtHQL3A6cMV4yI84HzARbY23K3xsyCbi8TrgPr3/QI638I69/0iDX179CigpF3LDMPAYcAviZOLLvB25gtREKnrwBOEevfdA7rvxjr33SOruu/7WTgLuBpQ2X76/LNYwOMRqXZhpVZUBoIdzcYiKSBsDCzoCrbM7oPaSDeOToeT2QgUhlIhYFYGZWkWXhOnB9Eh9R+AeYWVbZBkRFZb15GFvuU+8Lm63+GpsLyYAHiSz+JgXD3aB+VCVBpXZmFpf5FFtDJ9C8MxMJAqMzCpfoPcRUuGk7GsST0vyL03+b7Y/0Ps379l3wO6rwsqzW0pbYvzVSrdK3O9ULrKljAStP5XxqI1fm/8DeBNBCL/bY1EAsdqHPz/HHxm0AZg+W+2/3wnhffsWwRQKTr+m87zbkOeOLqk4h4DLAXuKFlu8ZsOZLKQNT06CHWv+kN1v8I1r/pDV3Xf9sjuALYERHPrp+fB7zPkQRMN6nuGWx69BDr3/QI638I69/0iG7rv2gyEBHfWocVOwv42Yh4H0BmHgNeDLw1Iq4Bng68ZkN6aswWYGUlGh9dxfo3psL6t/5Nf1lD/92PJpSZNwDnNLx2I/CsKfbJmC3JamixvmH9G2P9N7xm/ZteUKB/RxOaCm2MXYVGo6pqiwyEylSkjIGiLEWmYYCV3SLboDALySyCyhiksnvuHnW9SAPRfLk7RmUgXFEZCGWmY/EZiGW2uaXR/iwLQxLA8q7R92zu+GhZLrc7mbc1EEXEAlWWzidR6e+NmXllQ91XAD9Fdbvix4ALMqsejGsnIh4FXAh8M7AHWAZel5mfbNf7KbJevbfVugwMINpU+t81Ki5pIBRmwWYDYZn+l4T+1ZigzILLC4X6b3prhVtQ6l8FFRB6lSvqyiys2hMZTQHmhf5zcbRMeIonwvqfEkM6bDQBr7EdIE3BAKHGmMJzvTqHx24hLnWu3z06TqhAIQDLe0S28VL9CwOxChawrEz8oqzt+V8FFVBjivbwqmABo7XmGs7fIX57yAAiE5w/FF3W/9aZDBizxUmClfZXBg8CkZnPjIgnANdGxJMz8wuDlSLiTOAi4EzgXuCjwGuBtxe0cyawNzO/q27rQuC3gDPadt6YvmL9G9Nfuq7//q15GtOCHPNYi4iYA14NvAsgM28FbgReKaq/CvhQZt6TmSvAu6nvxy1o50bgZwfa+gvgscUHaYyRWP/G9Jcu698rA8aUkpDtjIKPA04Cbhkou5nKeDfM2cAHh+o9JSL2AI8e105mfnG1sK7/I8A72nTcmN5j/RvTXzquf68MGDMBmdH4KODU+u+RgbL7gVMa6g7XC+BAaTsR8dPAHcBO4I0lHTTGNGP9G9Nfuqx/TwaMKSRZM7TYgYi4fuBRGllgkssN4+o+7LXM/M/AycCngasiQjjljDElWP/G9Jeu6392twkNO/3bpB5v2lalHi+NMCBSirNTRBMREUZyQUQT2KOjiSwviGgCpVFCVNmussghK6Iecw13vqngPyqagPqqifc7hCV/RUQiWN4pIijt1H3MHSKSgSprE01o7WXCezJTLfmtclf994SB//cP/D9c94SB5/urHnA3sPrJr9lOZmZEvIXqHsLvBD4y7gA2jQK9F2t9I6KJCP2ryEEsCBGqaCJ79FC71EL/KiKILCvU/yTRRNRNsnMioocaO0JcRVORiGS/G85YK2qsEFGi5iY5xmGs/+kQMaq5wmhCUqsNkclkXVUmogRGaZRAca5PETlIRQ0DWN5dFjmwWP9KM6XRBBujiYkydf5XoboKo8bJ6GQictiy/hnFnNB/yHO92FEpHde/VwaMmYQ2DiK4jSoywBMHys4ArhN1rxP1bsrMB9dqJyJeEBFf/1CXM5eBo8C+ol4aYzTWvzH9pcP692TAmGKCXGl+rEUdFeAdwHkAEXE6VVbPyyLiQERcHREn1tXfCXxPRJxURw84lyqu8Nh26m2fXtenfv1fU2n9Ey0O3pieY/0b01+6rX9HEzKmlKTUKDSOg8AlEXEtlf5elpl3RsRjqBKI7AXuzczPRMQFwJVUa5sfBy5eq536td8Ffj4iPk41CMwBLxp43RgzKda/Mf2l4/r3ZMCYSWg5GGTmUQZm7QPld1CZfQbLLuOrs/2idurXPgO8pFVHjTGjWP/G9JcO6382k4FoMPcMo8x+haZAmMBAJIxm0iwojIEy9fje0bLlJgPhnkID0S5hkBNNpjIG7RA3tAkjXaOBUH3/VVVx05n0FIkltRXxsaQ4vhWdy5xU26vvSqGhqZGW6cgN2kCoUAEASjUNhNB1sYFQ6b/QQLi8r1z/y0L/SwvKQKgM9qPtKYO90r8sm8BcG6KqXClX+pfGYDW+lY15ACm+K9IQ2TB+FGP9tyYidICOYUrP/01tzYntdwj9K12rNmVggNFtl0XZ0h49Rimt6wAC6jfBaHvFpmL126Hl+X9uSehVBAsJESxE9UcFEJlrCiAgxrM5ofVsdEkXMl7/+yPiEHB5Zl7ebkebj1cGjCklafjFY4zpPNa/Mf1lbf3fn5ml4US3HJ4MGDMBKqSqMaYfWP/G9Jcu69+TAWMmob2ByBizXbH+jekvHda/JwPGlJI6OYoxpgdY/8b0l47rfyaTgSCIYdPenDB2qOy1ytjTZEYsNQYKU5EyECmzcCqz0L7RssV9TQZiYSBaGK23IjMQisyCykhbahZuykCsUFmJVT1l4hEftTQAq7IG/0+pgbDdxD46fWVgs2hlIBSaVjqv6paZhZUBUep/YVTXykC4tFfo/xF6jFoU+l8u1L8MFiB2I023SkcTeOtSqD3EoNBG15PpX20v+iNMjuVY/1MhhI5VBuJSs7/67QDlxmCl9V0qi7AoE4EBFveOli3t1X3UAURG6xUbg1UAARmco20GcmXOFb9HxD32IT4CZQxuff5Xv3ta/Zjvtv69MmDMJHT4yoAxZg2sf2P6S4f17wzExkxCu3TkxpjtjPVvTH8Zr//9EXEoIl44s/61wCsDxpSSOkeCMaYHWP/G9Je19e/Qosb0Bl8BNKa/WP/G9JcO639GGYgDhoy8MluoMgtJA6Bw0gC5u8xApMxCqmx596ijZXlhtGxp32jZYoOBaHGfMBAqU6EyEJVmEVS7VhPcpklvaQZiVU+YkpXZR/VRenUa+ijrqv20zECsMq+aCYkghs19Uar/sgyiAKnKS82Cu0b3owyEyyKzqDILNul/aa8yEE8327DWoEoN3NJAqNpUxj4VfGAjxiiZgbxh+0Ks/ykQc7B7yCWrso3LYAFlwT6AYq2X6l+d65cXRAZxYQpeFGWgg4VI/RcaiLVhX+lfdGaSm8aFWVid19WwrsaJFFqVQQEmCSBQmpV8Arqsf68MGFOKM5Aa01+sf2P6S8f178mAMZPQ4SsDxpg1sP6N6S8d1r8nA8ZMQJeTjhhjxmP9G9Nfuqx/hxY1ZhIcWtCY/mL9G9NfHFp0yswFMZTNt9QYrEzBiMzAACsLyhikMgsKs+AulRlQmIWE2WdJmopkF6WBUGYgFtmGlVmo2Cy8WRTuu9Qs3JgAsLRui+lvOLTgdIgYMRDKLKJiTGitf5VFVJkFhdbVmNBa/6UZyJWBsDRD5yRGXIWsK379KsEpx12hsXeSDOLSlCzL1q9f639KxOj5HxVARGQglgFAhKahPIuwDAyi9L9baFXUW1b1hKahwSysMhArs7DKIlxqupUBBES9SSjdvnQ8ah1AoLA/hRTo36FFjekLXV4mNMaMx/o3pr90Wf+eDBgzCb4dwJj+Yv0b0186rH9PBowpJbt9ZcAYMwbr35j+0nH9ezJgzCR0+MqAMWYNrH9j+kuH9T+jDMRzxMKQo0ZlFhQZRHPPqJNuea82EC7vEWYhYe5TZqHlXSoLYJnZRxmIVL2m8hVVJg5RG4hUdr/Cb/AkX/TSui3E0zoDscpq6AzEs2cuiF1DX+g2+t+jMxCvqOzghVlEi/UvdKn0rzRd9We0TOu/LNtwsYFY0ejOLzQLK200tjlEG/MhlGcbdgbi2TMX5MKQcGQAATUmCFOwCBQAsCzK5bm+jf7Vub5wTIAGrZdmG5cBBDbAGKxoo4NZamhj9b8/Ig4Bl2fm5e32tPl4ZcCYSfCPAWP6i/VvTH8Zr/9tHU3IeQaMKaW+Z7DpUUJELETEpRFxbURcHxHPG1P3FRFxQ13vooivXv4c105EPDYi3hURf1K//p6IOLnNoRvTe6x/Y/pLx/XvyYAxhQTtBwPgIBCZ+Uzg5cB7I+LUkX1FnAlcBHw38AzgW4DXFrbzo1RXKZ4DPAtYBg5NdLDGmIdh/RvTX7quf08GjJmEFhlII2IOeDXwLoDMvBW4EXilqP4q4EOZeU9mrgDvBl5T2M5fAm+vX0vgvcBzJz1UY8wQ1r8x/aXD+p9ZBuIRA5HKNqjMgvtGy5b26cNY2qsyhpZlEVUmoBVhIFJmH7mt9jhqs5AwCxZnGxWEMPFV35GRmrqB0vtk1exY+Zm2q4GofWixxwEnAbcMlN0MPF3UPRv44FC9p0TEHuDR49rJzMHtABaAe1r1fJqE0L8yC4psw8t7R8uW9jbof58yEJdlDJWBAXYK/ctxQm0ru0jK8aPULKjbHK0oytT3uKm9QrOwGmek1gvHiS2H9T8dImAoOEAKA3HuHP3Sy6ziIigAwNKesmABS9IEXKb18nqyi6yIoUtmG1b6b2OQL9VlEyoTr9pelMkxobRslnRc/14ZMGYSWlwZAFaX8Y4MlN0PnNJQd7heAAcmbAfgBdRXCowxLbD+jekvHda/owkZMwFrXBk4EBHXDzw/lJkl9+pOsl4xru7IaxHxT4HHA/92gn0YYwTWvzH9pcv692TAmFLWvgJwT2aqJb9V7qr/njDw//6B/4frnjDwfH+997uB1ZtYxrYTEY8F3gr8UGYuju25MWY81r8x/aXj+l/zNqGI2BERPxERhyPiYxHx5xHxXQOvn1WXXRMRl0fESWu1acx2pWU0gduAe4EnDpSdAVwn6l4n6t2UmQ+WtBMRB4DfBM7NzC9ExDdGREPqq2asf2O+ivVv/Zv+0mX9l6wMfAPwk8BZmXl/RDwX+EBEPJFqlvKBeodXRcSFwCXAS8a2qAxEwiykMosqU5AyCgMc3yfMQnuEWajYQCzKlLFHmYJFtlCAFJ+ANAaqzIKlKBOPWG0qzlTcQAhTUSyXGY2UmFqadTaENm9RZq5ExDuA84CrI+J04CzgFbV43w+8KDPvBd4JfLg+ud4HnEulrbHtAETEI4DfA94M3FU/fw1VeLHbJ+z25uh/hzD2Cv2rrOLKKAywqAIIKP0XZhGXZv9SA2CD/tX2xZlFS1FmX2UAbJuBXGl4WZTJ/pSVzdpUaP1PR//D53tpIN4tDMSiTGUVhwazsDjXKwNxG7OwHCd0F+VvhZxvkUW4VJdtU/GKNueWxPlfjgmF9Sb4TVA8VrQcP7qs/5LJwJeBN2fm/XVHPhwRR4F/ChwHljPzqrruO4HbI+LkzLy7oG1jtg/JZBEXNAeBSyLiWir9vSwz74yIxwBPAvYC92bmZyLiAuDKeq8fBy5eq536tX8PfAdwFQ9nPbHGrX9jwPrH+jc9puP6X3MykJlfpFpuACAiAthFdVXgXzAQ3igz74iIB6gSJPzxWm0bs50I2odFzcyjVLP84fI7gJOHyi4DLpuknfq1NwFvatfTh9qy/o3B+gfr3/SXrut/PaFFnw18jmqmMhz+CBpCHEXE+VGlTb7++PID69itMbMnsvnRE6x/01us/ynof8n6N9uTNfS/PyIORcQLZ9zNdTFRNKGIWAB+CfjR+r6lxqrDBXWIpUMA+/d8XX+GTtMttqCPYbOw/k3vsf7b63/v11v/ZnsyXv/3Z+b5m9STqVO8MlAvDx4CfjUzV2OpDoc/guZQScZsb8ZcFej6lUHr3/Qe69/6N/2l4/qfZGXgV4C/yMzfqUMUnUoVyuilqxVqE8Re4IaxLUWMRA9R0URU2coulU5cz2lURBAVOWhpYXRbGRGgMEqIihDSNO1SkYNkmnG58WiRjN6jAocop/4kAQZEJ6X7X0QYmFMRRpZEmaonypr2TYo3qG00gf5eGdxY/YtoYjLCkND/yk79xVXRP2TkoNJoYiryj4wcMvolUxGCoCFyWKEO5QlI6UNFDlJjwgT6lxGBpK7FOKG0Xqr/xmgihVq3/tfL1PSfIaIHqfP/vND6DhEFT5Q11V0WY4U6r8vIQYXRxOQ5vUn/bYL6lJ7/S8eJCdBRgkbL5hbF+V9EvJ8T+p9bGu34JOd/NSbIcWICuqz/opWBiPgZYCdwaR2m6PFUYY2uAHZExLPrqucB73MkAdNZcsyjo1j/xtRY/9a/6S8d1v+aKwMR8QSqLGYArxt46cLMPBYRLwYujohlvhoP1Zjukd2+MqCw/o2psf6tf9NfOq7/ktCitzJm4TozbwSeNc1OGbMVCbo9GCisf2MqrH/5uvVvekHX9b+e0KLG9JcOLxMaY9bA+jemv4zXf39Ci24oyr1aaq5ruG4hjTyqzVLDT2l7hR6ephfk4RQag0sJ6T6coEFlIFYGIlV2fHTbeWUqKjQaVftRZqOysmKyvQHJNFDqXlfamtPbav0LA2KhrtWYIA2AkxhxRZk09xc6DScKAjC834Zt5VChxi0ZQGC0bP5YoalwcXQn86IMtNlQjglt9Gv9T48hfanvnvw+Fuq3sVyOH2X1SstUvxujzahyGQxDtNmiP5Mg+y76OFcYLEBq/XhZvYn0L+qqesWsrf9tHVp060wGjNkGdCGEmDFmfVj/xvSXLuvfkwFjJqDL9wwaY8Zj/RvTX7qsf08GjJmEDl8ZMMasgfVvTH/psP49GTCmlI6HFjPGjMH6N6a/dFz/s5kMZI4YOVMZO6UBdLRakylUZeMrzWorM90pY5By8YjMt43GvpbmnvUzictZbC0z/pUZiOaFWWj+qCg7JgyExxsMRKJ8fnG0k7G0fjVXocU6fGlgMxnW7Ir4rFQGyQlM4XNS/2r7smzayryuvk3K7Jtzuo/SbKgqtshK3hZtIC4MIKAMhEL/O5T+VZnYtmpTmAWPj346c2JMKMX6nyJDH4MKaKGzXLfNSjtapoJcSBN/adZd0RfVHiAHEPmbYgOMwSO7mMDkXPzeFpqFdxwV5+9jo/Wazv/zx4TWlf5bBBAp0P/+iDgEXJ6Zl697RzPCKwPGTECXDUTGmPFY/8b0lzX072hCxvSCbL4KZYzpONa/Mf2l4/r3ZMCYSfCVQWP6i/VvTH/psP49GTCmFCcdMqa/WP/G9JeO639GBmJg8eHOMmWamT82asXJnaNl87u1k0aZUnS2wjID4cpOsaXKiCuOpSlLYiuUsaf0u9ryO11sICo0EM6Lz2rnA6pMO7Hmjwqz0LHRNb1YbLfO53uGp0AmsTz0OQijWcwLrR8fLVsR4wTAyo7RD2tlvjA9qHDnSf0r86HMaLpJkQJafD8nMhBKc6cokwEEhFlQmIV3Pjhab8eD7fQ/Z/3PnEgxDgvRhHD2zu0Q2auXtP5lBlqVRVxlNVYftAqQofQv9yG72DqD+bqZ5LdDqYFYmanVZyDG+h0iWIg0FYsygDlhIFam4lhuFw6oy/r3yoAxhVTRBGbdC2PMLLD+jekvXdf/RlyvNqabZBIrzQ9jTIex/o3pL2vrf39EHIqIF866q+vBKwPGTILP+cb0F+vfmP7S4dCiXhkwppSskt40PUqIiIWIuDQiro2I6yPieWPqviIibqjrXRQDxpq12omKN0REQ5omY8xEWP/G9JeO6392GYiPL46UlTAvDIA7hQGoiVgZnf8sqcyiu0bLlncKo1GpWaitKajwipTOFrr+es3bl2WClAbiUgORMAvveEAbAHc8MNro3IOjZW0NxFO4MngQiMx8ZkQ8Abg2Ip6cmV8YrBQRZwIXAWcC9wIfBV4LvL2wnd8A/h4QttcZk0kcG/pshLFLXqmQZr8GcUmvsNK/MCXuEuZj8U6mGI8m0X9pFtGJ9LrebRvak1ovzdCqDMTCVKiyje9QBuIG/c8/OLqjeaF/WmQgB6z/aZBJDAUQQZyXS/U/36B/rS1hVJZmWBFAQAQkkOf/CQzEisbxrAClVb0Tte0EdVWwgMIAIqX6VwZgFSigqivM5sdFJ0Wm+4nosP69MmDMBLS5Zzgi5oBXA+8CyMxbgRuBV4rqrwI+lJn3ZOYK8G7gNRO080bgnes8TGOMwPo3pr90Wf+eDBgzAZHNjwIeB5wE3DJQdjPwdFH3bFHvKRGxp6SdzPx8UY+MMcVY/8b0ly7r35MBYwqJXPPKwIH6/r3Vx7CZ6NT675GBsvuBU8TuThX1AjgwYTvGmClg/RvTX7quf0cTMmYC1jAK3ZOZapa/ZrNTqrtJma2M6SfWvzH9pcv6n5GBeAWOPdzkHEsqlZ/IXiea29Fwv1Ysj3on5o6PuntkVlNhFi7NaNrWQNgmi3C5gVhlXdW7kdlGizMQqwyEouy4MAsfHf1OzD+gUj5DKLPwsVEjvfyelZK0NRDdVf89YeD//QP/D9c9YeD5/nrvdwO7J2hn67GSo/pfVq478T0R37FG/QvH79yi0P/CqAhVsIAVkf0051WWc2Vyll0sH76nrfXCfUCD1qWpUHxeKqiAyEA8t1im/zlhFAaYO6r0L+oOZ76eBOt/OqyMBhBJkW1cfW/n1He5Sf9qezHOzC2KTOeLZfpfaRlARJuF1/8lm+S8PrptQ3lpsADxORSf/5VZWIwJKqs4QIhyFSwk2hiIO65/3yZkTDGtkw7dRhUZ4IkDZWcA14m614l6N2XmgxO2Y4yZCta/Mf2l2/r3ZMCYSchsfqy5aa4A7wDOA4iI04GzgMsi4kBEXB0RJ9bV3wl8T0ScVEcPOBe4ZK12pnegxpgRrH9j+kuH9e/JgDGlZPukI9TxgSPiWuA9wMsy805gD/AkYC9AZn4GuAC4EvgE8Eng4oJ2oHrhTcB76/8PCzOTMWYSrH9j+kvH9W8DsTGT0DLpSGYepZrlD5ffAZw8VHYZDbP9pnYGXv9F4BdbdNUYM4z1b0x/6bD+ZzMZWEny6LGHlykDkTB7hchUquoBhDSgjB7yyu5RF5A0C+0c7aPMQCpMRRP5vEsNP8oL08YsPJGBeLRQzY5lPfW5iGyBqmwkc+1q+VGRdXs4yyW0MxDS0oBkKnKFPDak/6VR0Sizt9R/w2dSqv/5o5ug/yZa+AdnqX+1LK4NxFPW//GGAAKl+m+pX+t/CmTCkIE4lJFWBXtQmcobrsoWn3t2rV//6t4KGUCg6R6MNvFfNmOcaKpbaCAOaSAWY7jIDD4ntK6MwsBoRmuQ53p1/piENfS/PyIOAZdn5uWtdjQDvDJgTCGREy0HGmM6hPVvTH8p0P/9mbltb8fzZMCYSSgN02iM6R7WvzH9pcP692TAmFIS8JVBY/qJ9W9Mf+m4/j0ZMGYCihM4GWM6h/VvTH/psv5nlIE4yeNDhq854bBRZg9h4GjKKjsnTCVzR0cPOXeJsp2jpqLcIQyE0iw0/azwpSY+SalZqDGTY+G+pYG4zCykzGLSFKTKABZHjcW5pExFbQxE2dqAaKgCCAxnIBYGwtw5qksdVGAC/QsDcan+V4T+Kc02LDONTkCh1ou1OomBUBkDVVWldXUlTehffobi85NjAowYUgFSnRdW2gQQsP6nggogECqAiAgqoAIINJz/c3k0A3kcFxmI1blelYkgJ/L8L4IKtDIKT0LxuV5t2+78rw3Eped/UW+CACC6rthPqx/z3da/8wwYU0rSKumIMWYbY/0b01/W1v/+iDgUES+ccU/XhW8TMmYCHE3EmP5i/RvTXxxNyBhTG4i6u0xojBmD9W9Mf+m4/j0ZMKYY3w5gTH+x/o3pL93WvycDxkxChw1Expg1sP6N6S8d1v9MJgOZORLpRUYTkdEpRt3kMmoMEIuj0QTYISKUiKgluWM0moCMeCQiDEwSOSRF3eLwVTJKSHmUoKJtx5UPIVN9y/6URROSn2vDZ62iDGRhNKpikvL30jSTORL9KZW2ZMQqEU2iYek2RIQZdokII0dV5DARdURpXfRbabp1qIbCr60cO0rLGve9/jblmFCof6V1qemGunJMaHNlz/qfDpmj0Z9UVJ4VcQ5WvwkaxnRZvigiFKnz/7zYt4gSJMctFU2oidLfCm1+EyjEW9P4u6NY64Xn+tKIUGqsb4gcVXz+T5//m/DKgDHFZMvQhMaY7Yv1b0x/6bb+HVrUmFJWrww0PYwx3cX6N6a/rK3/7ocWjYjXAS+mSp1xMnAoM/9r/dpZwK8Dy8B9wLmZ+cWN6KwxM6fD9ww2Yf0bU2P9W/+mv4zXfy9Ci74aeG5m3hURpwN/HRGfAD4JfIBqALgqIi4ELgFesiG9NWaWZDZmQOw41r8x1r/1b/pLx/VfOhn44cy8CyAz/yYi7gNOA74OWM7Mq+p67wRuj4iTM/Pu5uZG31S5yFpqIGwyvqi6yhi0KEzF0hgktpVmQZWOXBuFZGmxCaiFMXAiA2Hh1TDRpjTslZpAhdmn0UAoDUQbINwOhxYbw1T1n4wGECCE3gr1Gw23aOSS0KswoIUwBocINFAaLEAFQ2hE1Z22htsEFWjavnTfpSbe0jGhQdNS69b/tJju+T9FcAilI/UdFWXZ8F1WBuJYFrpWptTCACJq7JCanmRMKGXa+p/ku10aLKBFAJFiUzDoe/mV/tveztdh/Rd5BjLzr1b/j4gfBP4RuBI4G7hloN4dwAPAt0y3m8ZsBbL6gdL06CjWvzFg/Vv/ps90W//F0YQi4qnAbwP7gH+VmUci4lTgyFDV+4FTxPbnA+cDLLB3vf01ZnYkZJvQZNsY69/0Hut/OvqPfRvbWWM2go7rv3gykJmfBs6IiKcBV0TE942pPrImlpmHgEMAXzN3YnfXWky36cAVgPVg/RuD9T8F/e+fP2D9m+3JeP3vj4hDwOWZefkm9WhqTJxnIDNvjIg/BF4PfA542lCV/cBdU+ibMVuLzF5GExnE+je9xfq3/k1/WVv/3Y4mFBEnAedk5u8NFH8FOBG4DnjpQN3HAHuBG8Y2ms2mn4ftG2EgEUbDQGQaRWcwLDYlKWOQMrkIs/BEBsIWFGfTLM5APMGJrs2+S02FMnthg4Gw2CTV7mTe1pQcEQtUETeeRKW/N2bmlQ11XwH8FJW//mPABVkf6FrtTCvk32bpP+ZUplqR3VeMCc37ESY3Wa0s+3EsFRoDVfCBScaEYm2t39ivq7XNSlxqICzbVpoFm/Qn3oviMWECrH+gtf7FjyqVyVu815OcWdUnrQN2iN8Jql6I75gYozbEQDzloAJNQReK2ywtK9Rwlp7rm67MTxJspAVd1n+JgfiRwM9FxN56JydTxRz+KHAFsCMinl3XPQ943/hIQsZsU3IqBqKDQGTmM4GXA++t7719GBFxJnAR8N3AM6hMea8taScidlGF/HtjZn478JdUA8d6sP6NAesf69/0mI7rv2QycCdwOfDRiPgY8BHg3cDFmXmMamB4a0RcAzwdeE1Bm8ZsT3Kl+bEGETFHFbP7XQCZeStwI/BKUf1VwIcy856sXEvvptZWQTvPZzTk3w/UJ/JJsf6NWcX6t/5Nf+mw/te8TSgzjwI/Vz/U6zcCz1qrHWO2O5nZdpnwccBJDITjA26mOokOczbwwaF6T4mIPcCj12hnJORfRKyG/PvjSTps/RtTYf3L161/0wu6rv+iPAPGmIpcycZHAavLgUcGymQovrrucL0ADhS0UxzyzxhTjvVvTH/psv4njiY0Db7Mffd8ZPm3P1c/PQDcIytuv8zPzcey/ejTsTy2pJEvc98ff2Tl/zswpspCRFw/8PxQHVJvLSZxl42ru1Y7m+NsXwOp/+2ndUWfNLOdsP4n38+G8aW8954rH/wfa5//tx8+lq2J9V/w+kwmA5n50L1LEXF9Zqplkm2Hj2VrMq1jycx/2bKJ1ZB7Jwz83xSK7666HgP1Ergb2L1GO3exhUP+Wf9bHx/LKNb/dLD+tz4+llG6rn/fJmTM5nEbcC/wxIGyM6hC9A1znah3U2Y+WNDOw7YtDvlnjNlIrH9j+suW1r8nA8ZsEnVUgHdQheAjIk4HzgIui4gDEXF1RJxYV38n8D0RcVIdPeBc6vBg49qpt3XIP2O2GNa/Mf1lq+t/JrcJDVFyT9V2wceyNdlKx3IQuCQirqXS38sy88569v4kqhn8vZn5mYi4ALgSWAE+Dly8VjsAmXksIl4MXBwRDyUd2YRjWw9b6bNpi49la7KVjuUg1v8gW+mzaYuPZWuylY7lIFtU/zFR1kljjDHGGGNMZ/BtQsYYY4wxxvSUmU0GImIhIi6NiGsj4vqIeN6s+rIeIuLsiLgtIs4dKj8rIv48Iq6JiMsj4qQZdXFNImJHRPxERByOiI/V/f6ugde3zbEARMTrIuKjEfEnEfHpiHjdwGvb6li6jvU/e6x/Myus/9lj/ZuHkZkzeQBvBX6j/v8JVO7oU2fVnwn7/v3Ab1G5s88dKN8FfA74zvr5hcDvzLq/Y47jNOCzwP76+XOBf6TKcLetjqXu46eAU+r/TweWgGdsx2Pp+sP6n/3D+p99n/v6sP5n/7D+Z9/nrfSY1Yc2R5UE4jsGyj4K/LtZvyGF/f+G+u/hocHgRcDfDTx/DFXqtJNn3eeG4zgJ+OGhsnuAl2y3Y6n7+E+Gnt8N/KvteCxdflj/W+Nh/W/dY+nyw/rfGg/rf+seyywes7pN6HFUX8RbBspuBrZFkovM/HzDS2czcEyZeQfwAPAtm9GvScnML2bmb64+j4igmkXfzTY7FoDM/KvV/yPiB6muclzJNjyWjmP9bwGs/617LB3H+t8CWP9b91hmwawmA6fWf48MlN0PnLL5XZkqp/LwY4LtdVzPplpO+zjb9Fgi4qkRcTPw/wAvzcwjbNNj6TDW/9bE+jebgfW/NbH+e8xWyDMwSMy6AxvElj+uiFgAfgn40cxcqS4S6Kqb16vJycxPA2dExNOAKyLi+8ZU39LH0kO6+nls+eOy/s0WoKufx5Y/LuvfzGpl4K767wkDZfsHyrcrd/HwY4JtcFz18uAh4Fcz8/q6eFseyyqZeSPwh8Dr2ebH0kGs/y2E9W82Get/C2H9G5jdZOA2qugBTxwoOwO4bjbdmRrXMXBMUWWV20sVdWAr8yvAX2Tm70TE7oj4RrbZsUSVtvsHh4q/Auxjmx1LD7D+txbWv9lMrP+thfVvZjMZyMwV4B3AeQARcTpwFnDZLPozRa4AdkTEs+vn5wHvy8y7Z9insUTEzwA7gUsj4hHA46n6vd2O5ZHAz0XEXoCIOBl4MVWUiu12LJ3G+t86WP9ms7H+tw7Wv1kl6lBLm7/j6h61S4AnUXkX3piZV86kMxMSEd8KXEQ1gN0J3JyZP1C/9jTgYqrQVfdRhR774oy6OpaIeALwP8VLF2bmwW12LAvAm4B/ARynWhZ8P/CW+h7IbXMsfcD6nz3W/9Y8lj5g/c8e639rHsusmNlkwBhjjDHGGDNbZuUZMMYYY4wxxswYTwaMMcYYY4zpKZ4MGGOMMcYY01M8GTDGGGOMMaaneDJgjDHGGGNMT/FkwBhjjDHGmJ7iyYAxxhhjjDE9xZMBY4wxxhhjeoonA8YYY4wxxvSU/x9Lrgg6xaUNEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 936x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwMAAADWCAYAAACNOV4iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABBXUlEQVR4nO29fbhkVXnm/XvqfPTp0/QHdPMlKgaDIGIGFRzNJNFoNMYxasyYaDCZDhovx4xJXnXU6JvYmonjZMI7ziQSBtBBHaK+8SNqFEURBCeioGgUQsyQQDCK0DQNTX+cj6pn/tj7QHWtZ9dZdapOnzp737/r2tc5tWrttdeu2vfatfbe9/OYuyOEEEIIIYRoHq217oAQQgghhBBibdBkQAghhBBCiIaiyYAQQgghhBANRZMBIYQQQgghGoomA0IIIYQQQjQUTQaEEEIIIYRoKLWaDJjZaWZ2tZkdMrPbyv//2sz+1sz+zMyOXqXtfsjM7jSzS7vKnmRmd5nZCQO080Ize+EI+vOupf3vU6f3s3r/sNsVYtwws7PN7Mtmdo2Zfd3MfjVjnZea2efM7Atmdq2Z3WBm/yZzezvNbFdm3WeY2ddKHX7TzJ6Xs17X+r9arrfq8aHN7Dlmdp2ZuZk9qqLOCV1jys7V7pMQy7FC/R8qj+Pu5WYzuz1jXelf+l+XTK51B0aJu/8d8HQzuw241N13QXGQAtcDjwAGElzmdl/SPREo2Qd8F5gboKkXln//csj+/I6Z7QWe3qdO+FkJURfM7CTg88B57v5xMzsZ+JaZ7XH3T/dZ9XXAf3f395ft/Dzwl2b2L9z9O8F2ZoE3A/+pq2wSeAPwPnf/54rtXAL8gbv/TzN7CrBlkP1z9w+Y2R3AVYOstxLc/bNmdgvwj33q3MlDY4oQa8oQ+r/T3Z/e09a7gQcqtiP9I/2vd2o1GajC3e80s48Cv2VmR7l7KOoRb/O7wE+s9naEEJX8FnCXu38cwN1vN7MPA28F+v0YeA3FxYMlrqa4i3oqkEwGKMbRCeDLwFeBrcA1wE1Av6t2P0J5cnX365bfHSHEAKxU/8/ufmFmG4CXAE+tqC/9i3VPrR4TWoZpoA0smtll5WM97zOzPzCzK8xs/9LtvfIRn6vL22J/bWb/zcw2LTVkZhvN7GIzu93Mvmhm/5muz9LMzirXdzN7elf5MeU2v12+/xUze0353p8AzwGe03Vr8oTyvUeb2V+VtxS/XLaxo6vdlpm9s+zPl8zsQmB2JR9S12fzfjN7R9neLWb2bDN7nJl9uLxl+gEzm+5a7yfN7Mqy/v82s4+UV2a6236imX213P8vmNkvl5/RN628fWtms1Y80vXN8jP4vJn92Er2RTSeZ3P4j3ooTtbnmNkxVSu5+1fcfRHAzKaA/wDcTHGVMap/v7v/LsUPhhdQ3H18g7v/hrt/v7e+mT3eHnqE713lcX6umX2r1MOXzezxZd13mtnHuta93Mz2mdm7etr812b2KSseiXxQm2b2LHvo9v6LzexjpZ4vM7MNZvZmM7uqXO8ZlZ/k4ZxtZp+24rGLb/dbz8zeZj2PLJrZJ8xsb9d4+/iu8fLXyzHma9bnkQQhMlip/r/bU/RC4KagfKm+9F+B9L+OcPfaLcBtwK6u12cCe4D/2VV2KbAX+Jfl65cCbwJ2lHVfWpZPA58DPtC17p8AtwDbytdPo3gs6NKefjjw9K7XVwEfBSbK1z8L3NvTp942poG/B363fG3Ae4Bruuq8DrgbeGT5+nTgLuDqQT+rrn7sAU4vX78K+D7wpvL1RuB24Fe71vmPS++Xr98KfKHr9SxwJ/Cfu/bjfwSf0QeAvwKmytf/FrgX2L7Wx5WW9bUA9wN/3FP2c+Uxd07G+u8uj72/Bh7Vp95m4B3AjcCFwAfLdf4MOKHPer3H/vFAB3hGV9nNFI8nzJSvHw78Zdf7Ty/beXP5ehb4Z+DfdtV5VFnn3eXrmVL3l3dp/DeBW5f5PJba+WiXPv8dsB84savebcDOrte7esciirstu4LP40pgtnz9KcoxTYuWQZdh9d+1zmeBX+vzvvQv/a/7pc53BnaWs80lUb6T4kdtN99y968CuPsH3f2dwL8HDrj7B8vyeYpn+8614sr+UcArgYvdfW9Z50vAt/t1xsyeRiHcd7p7u1zvc8BFy+zHr1CI8F3lOk4x4PykmT2hrPNa4EPu/k9lnVuALyzT7nJ8o2wH4FrgRIrboLj7QeBrwJO66v8p8N+6Xn8IeIaZbSxfn0sx2P1R136c371BM/uRst673H2hrPc+YAr4tSH3RzSPTaSenbmu9/ri7r8JbAc+DnzNzM6oqNqmOJH9BHAd8HfATwHfY4BHMd39hxRXLp8PYGanUviOZoBnltWeD3wyWP1/lW0cKNt4QlBnaUw7RHHFdKpL41cDp5jZtoyuXrikT+BiYB54dcZ6Ofyvch9w959fGtOEWAFD6R/AzB4OPBn4SJ9q0r/0v+6ps2fgUl/eFBsdaI8HNtvhkXhmyronUTwONA38n571Ko01JUuPuvx9d6G7v3GZ9R5PcbXgcjNbKpukuDJ/gpltAR5W0Z+HLdN2P37Q9f/+irKtXa+nKG55Ph5YADZQXP0/ruzrGcBed7+np4/dnFmu84dm9v92ld8JbFvZbogGs5/iOOxmQ9d7y+LuHeC/WPEY29uBJKpQefJ6C8CSRr14zOgPV9DnTwK/AfwOxSMHlwJHU/wI+DTw8xR3y3rpNinez+HaXKJXv4d6XlOut3eZPt629I+7L1oRZeVxy6yTi07+YlQMrX8KrX146QdqhPQv/deBOk8GcmhXlN/mPdEEljCzfzHkNlcSBuxQn/4MFIFgAKLPprfMuv6/nOJRop9x9wPls37/2FMnl1e5+40rWE+Ibm4FekP7ntj1XoIVZ/Op8o5gN38LLKt9d790wD728kngHWZ2JsUzzxcApwCvN7M3UDwucFew3X7aXKK3TqTxlegV+o9r0XtV556qMVmIQRlY/wE7KR4hzkL6z35P+h8z6vyY0Er5FvAjVkQQAIowYaUpZ4biCvwc8KM96/1IRrtQRCR5EDN7kz1kTu50lc+UffgWsMXMHtaz3rvN7GHufj/F7chB+zMyzGw7xZWBT3VdQZnuqfYdYFtZd4nePn6bYuA4vaf9V5vZs0bYZdEMrgDO7il7MnCDu++pWOdk4BtB+cM4/OrbquDuN1GMMS+neFzxAPAJih8xv0fx/PJa86ilf6wIoXgyxbPNVdxPGjbx4aPvlhCHsRL9P4iZ/SSw391vWI3ORUj/Yq3QZCDlTykec3ldV9nrKEy/h9x9P4UH4TeWnq8rB43eQecw3P0aiufy3mhmrXK9XwKeV7YJ8ENgKcrBuygGhA9RPDf41vKqJWb2y8AT/KFIBX8MvMTMHlm+fxrw3JXs/ArZQ3EL8plL+wb8Yk+dP6e4c/CGrrLDnjN099uA9wGvM7PN8OC+vJ5lPBlCBPx34HgzewFAqY9fpjC0UZY918x+0HPH74yldco6zwP+Fcv7e0bFpyi8S58GcPdbKe5M/Dbx88JHmtdYEWUJikcapimuYFbxDeCxZnYcFMkVKbwYQqwmK9X/Er9OEazjSCP9iyPPWjuYR7kAp1H84D5E8Vzb1QRufoooNneWy9XAST3vnwV8Efgb4EsUPwKO6np/piy7nSJC0J9S/Ni9E7isXP9qiqvc3wR+pVxvO/B+iqvkXwI+Bjysq90fpYhIcE35/jFl+aMorg78bdmvD3XvF8Wk7h1lf66h+EH9Lopn/67M/KzeH3w2f1Luy3XlvlxHYRp+V1ed/1Gu9+PADWUf/5LCsL20zpO7PtevUvyw/yzwr8s6/6qrXxvL7f5d2b/PL62vRcugC3AOhfH9WuDr9EQFoQgbuBd4Uvl6BngjhcHuf/NQ7PCXjLBPj+8ZH/6q5/2nUdwl7Nb4fwL+vqfeL5Tre9neCT3avIAiNno//UYaf2rQ5+d01dlJEWHtGxRj2TPLOid0jSm3AK/tWv+/lGWfp5jcX00x7vwRxR3C7s/jT9b6uNFSj2VQ/XeVb6KI0Hf0KvRJ+pf+x26x8osSYtUxs+O863lHK/IQfI9iMpbEYxZCCCGEEKuLHhMSR5K/NrNHdL3+f4DPN2kiUHpBLi0TwdxgZs/uU/fcMrHLDWZ2/tJjYjntWJE47uqu5T+u5n4JIZZH+heiuYyz/pseTUgcWT4MfMrM7qUI8fZdirwCTWIXYO7+FDN7DHCdmT3WixjTD1JGkzifhxLmXUnhsXh3Zjufdfedq743QohB2IX0L0RT2cWY6l93BsQRw93f4u5nuftPu/uPu/tOd797rft1pCjN1a+gNKV5kd7+RuBlQfWXA59x991exLp/L2XSvAHbEUKMAdK/EM1l3PU/ksnAILc+hGgwp1CYyG/pKruZOBLVOUG9x1mR0Tmnncea2eVm9mUzu8jMdoxiByKkfyGykP6FaC5jrf9R3RnYRXnLAvgV4ENmdvyI2haiLixpYm9X2X0UWZqjur31DNiR2c5NFGH0frJ877PdzxyOmF1I/0Ish/QvRHMZa/0P7RnoumXxIihuWZjZ0i2L86N1pm2Dz7ApekuII84+7t3t7scuV+9nf3qT37OnOkHi1/9m7iYOT/F+kbvnxKYf5CTdr+6D77n7eQ8Wmu0CXkuRcOerA2xr+c5I/2KdI/2vHOlfrHek/4JRGIgHufUBwAyb+Jf2zBFsWojh+YJ/5Pacerv3LPLXnz2p8v2Zh/3jIXfvl3xuKazqtq7/t3b931t3W9frrRRxmO+mMF/ntoO77zezPRSZIkf6YwDpX6xzpP+hkP7Fukb6LxjFY0JZtz7M7JXl84Q3LDA3gs0KcWRxoINXLhncShEZ4LSusjMoElz1cn1Q7yZ3P7hcO2b2DDN7cFAys2mKgWM1QrhK/6IRSP8h0r9oBHXX/2pGEzrsdoa7X+TuZ7v72VMPTmyEWD84zoK3K5dl1y+iAlwMnAdgZqdSZH+8zMx2mNm1ZnZMWf0S4Llmtr28Fb8TuHC5dsp1Hwn8u65nBF8D/AOjvyrYD+lf1ArpfyCkf1Er6q7/UTwmNMitDyHWNZlXAPqxC7jQzK6j0N9L3f1OK5KxnQ7MAnvc/Ttm9nrgCorU9NdQpJfv20753heBnwKuLQeEfcDPu/vCsJ0PkP5FY5D+E6T/tWLV/OAZ+NA6WJfUWf+jmAx037JYGgDOAD4zgraFGBscWKAzXBvuhyhm+b3ldwDH9pRdxkOz/ax2yvf+ifKqwRFA+heNQPoPkf5FI6i7/od+TCjjloUQtcCBtnvl0kSkf9EUpP8U6V80hbrrfxR3BqD/LQshaoHjLAx/m7CO7EL6FzVH+q9kF9K/qDl11/9IJgP9blkIURfcYaG+Y8GKkf5FE5D+Y6R/0QTqrv9R3RkQogEY7YHygwgh6oP0L44AkTHY0ie6rZV5LAbrhnj8PLx3ol/Amc/O1+DxmYeot/41GRAiEwcWvL6DgRCiGulfiOZSd/1rMiBEJg61vjIghKhG+heiudRd/5oMCJFJcWVgNfP0CSHGFelfiOZSd/1rMiBEJo7RXtWk3UKIcUX6F6K51F3/mgwIkUndrwyIdUZuBtLIQFhhFgyplQlw5Uj/YuQEGraJibReUBbWi8aEVuYx226HxdZJxwpvB+NHMKaE5uNBxp5k3bUbi+quf00GhMjGaNd4MBBC9EP6F6K51Fv/mgwIkUmRjjy4GiOEqD3SvxDNpe7612RAiEzcjQWv72AghKhG+heiudRd/5oMCJFJEVqsvrcJhRDVSP9CNJe661+TgfVMZD6anErLJtIDODIA+eJCvB0ZCIEimsCCSzJiDcjV+nRaFhoII1PgfKx/j4yFndhsWGekf7FiKsz+kQnYJoNjbCpP6xbUyw004MGYAEA0LiwupmXROBGOHZm/J3INyRV143or/y1Td/3Xd8+EWAXaNc5AKIToj/QvRHOps/41GRAik7pfGRBCVCP9C9Fc6q7/+u6ZECOm7s8MCiGqkf6FaC51178mA0JkUlwZqG80ASFENdK/EM2l7vrXZECITOoeWkwIUY30L0Rzqbv+NRlYL7TSg7C1aTYt27olXXcyOIAPzSVFnfv3hZv2+fm0LHL11zzCiEOtMxCKMSGI/tHasCGttnlzWnZUOib4VDrM20IaDcQfOBB2xw+k5T6Xjh/hmHAEonwcKaR/kUUUvccqjpsomtD0dFo2uzFdd+NMUuRRNLHo/B8RjAkAFkQTsrngN8Fiev63zKhDHuk/qGcVEY/C6IhRJKPI/5v5s6Xu+tdkQIhM6n6bUAhRjfQvRHOpu/41GRBiAOpsIBJC9Ef6F6K51Fn/9d0zIUbM0pWBqiUHM5sxs0vN7Dozu8HMnt2n7rlm9vWy3vlmD91/zm3HzN5mZrcNuq9CiMOR/oVoLnXXv+4MCJGJA53hnxncBZi7P8XMHgNcZ2aPdfcfdlcyszOB84EzgT3AlcCrgXfntmNmxwGvACpSSwshcpH+hWgudde/JgPrhFZgFrITj0vKDj7q6KRsfnM6a516IHXNbLzj/nDb9r07k7JOYCCkE8yOo5TikbEnYsxMhcNGEzCzFoU4X1S05981sxuBl1EIv5uXA59x993luu8F3gC8e4B2fh+4sGxLrBNsMjUBRmZhP3F7UnbouE1JWXtjesxOzKW6nL7nYNifid3puOD3P5CWRabiXK1HBIbkyvYio/KIxw/pX6wUC4zCEJuFQ2NwEBigszkta29Kxw6fyvsB25qPzbmtudQE3DoY/L6MjMaRKTkw+0ZGY4/MwsE2ivLU0GztdL8jo3G2gbjm+tdjQkJk4jDsbcJTgO3ALV1lNwNnB3XPCeo9zsw25rRjZo8GHgZcm9MxIUR/pH8hmkvd9a/JgBDZGG1vVS7AjvL5vaXllT0NHF/+3dtVdh+Q3uIp6vbWM2BHZjt/ALw1f9+EEP2R/oVoLvXWvx4TEiKTpSsDfdjt7tEsfzmi6McrqWsAZvYkYMHdv21mT19Bf4QQPUj/QjSXuutfkwEhMnGMjg+i24S7yr/buv7f2vV/b91tXa+3UoxHdwNLGaiq2nk78JvDdFQIcTjSvxDNpe7612RgDLHJ9GtpHZMagx94TGog3HNGuu78ltRIN7k/neFuOeaYsD9bN6SmpNb3guM3MhUPQWQW9CijIXGW5NEbCJe9MrAct1JEBjiNh4R7BvCZoO71ZT266t3k7gfNrLIdM9sM/ChwaRmJbBtwgpldDXzY3f9smB0QIybILG4zabZhjtmaFB18RGoqvv8Rqf4XjkpPYBOBVGfvjk8Hm2ZTk+PU7tTk2DpwKF05yEoaZmiNiEyFh4JtAJ2DafmoxwTpX2QRZBu2ifiJbAsyBkfZhtuBWXjhmFSDC0GwkMUNmQbixVgbkwfTPk4eTMeE1lyq9dZCYNhdDAzEC0G24czMxwAcCga0YEyIshpT0WQvdde/JgNCDMAwVwbcvWNmFwPnAdea2anAWcC5ZrYD+DjwAnffA1wCfN7MtgP3AjspIgP0bcfd99E1iJS3CS9196evuONCCED6F6LJ1Fn/MhALkckoko5Qxgc2s+uADwIvdfc7gY3A6cAsgLt/B3g9cAXwVeCbwAUZ7TyImf0+8C7KKwNm9uQV7LYQAulfiCZTd/3rzoAQmTjGYpRLYZA23A9RzPJ7y+8Aju0puwy4bJB2euq8neL5QSHEkEj/QjSXuutfkwEhBqAzkPFfCFEnpH8hmkud9a/JwDgSZCvsbN+SlN3/qPTre+CUwGA7lRp25gOv0OJRqVEIYHEm3faWo1Pz0tS+PCdOJ8iIGGVJtCAj4tQP74vbvDM1NHcOHMjqTy7usDDklQEhurFWenKxmVRbi9vTzMIPnJjqf/9J6TYWN6U6ai2k2104Kn5qdGE27c/stsBUuD8wAUZZhIN9js6xE4fS9ibvjTXduvvepKxzf5A5OTQVh02m1aR/0Utgho80HZ3TAWwqyBi8MQ0g0N6SGnbntqX6n9uSanhxY57eqgzEE3Npm5OHgqzmwY+K1kJaZsF2oozoE4fS3zIT++PfGK19wecbBSqoMiBnUHf9azIgRCYjCC0mhFinSP9CNJe661+TASEycWDR5bkXoolI/0I0l7rrX5MBIQagU+PBQAjRH+lfiOZSZ/1rMiBEJu5W6ysDQohqpH8hmkvd9a/JwBhi06lZaP6YNCvhgRMCw09gFo43khbNHxdn973fA6PS0anJKcpK2JlMN9ROd4V24F2ePJiWbfuHIDsrkOZiBb/j+2nZwhAGImCxU9/BQKwBgbHQAsPu/NbAQHhMqq3Fzan+OxtTI2707Gtk7K8qn9+SjgkTc3nmus5U3nO3kwfS8e2oH8RBDtL8rNAKtN6JMpBmDpnSv8gizEBcoY3AQNzZmJYtzKZ6mw8M//NbgjEhEEf0m9Y6sS4nglNmaz7Iah4YiC34SdEKyqJ1pw6kn9n0dPw5RqNCK8g2brnZzwPqrn9NBoQYgDobiIQQ/ZH+hWgudda/JgNCZOLU+zahEKIa6V+I5lJ3/WsyIEQuXu8rA0KIPkj/QjSXmus/e5pjZueY2a1mtrOn/Cwz+4qZfdnMPmVm20feSyHGgKVnBquWOiP9i6Yj/Uv/ornUXf9ZdwbM7BeAFwN7e8qngU8AO939KjN7G3BhWVfkEGUwDDKQzgfZBhe2ZjrfcmnFGQgXtgSmxOkgi3Awae5Mp+v6hqDfgfHZDgbmyk58yG68M82SPLF7T1LWHspAXO+kI1VI/yOgwrgWGQt9OrXDLW6KMoum7XUibU1Euk7LOkF7AJFi2jPBuBVsxgO/X2cyb9yaeiA4wQbBDACm7007P7UnHUc5eCgtW8jqjvQv/acEZmGiDMRTFT+1gvLOdFrW3pin/8gsvLgpyAIedadClu3IBBxkMLfFoCxosxXobWIuCmiQr7XWQhpYxObTjltn5b+Z6q7/3MeErnf3j5vZ1T3lPwe03f2q8vUlwG1mdqy73z2qTgoxLrRr/MxgH6R/IZD+e8qlf9Eo6qz/rD1z9+9VvHUOcEtXvTuAA8ATh++aEOOFO7Q7rcqlrkj/Qkj/AdK/aAx11/+wBuLj6bl1CNwHHNdb0cxeCbwSYCaMCi3EuFPv24QrQPoXDUL670H6Fw2i3vpfrWhCySfm7hcBFwFssWPih9OFGGMcanEF4Agg/YvaIf1nI/2L2lF3/Q87GbgLeEJP2dayXOQQZitMy9rT6YzUp4NsmsMwwBDdmQpMSWFZYNiZzNuQz6T7t7gxzkDYmUrLJ6Ksj5GRM3e/vbhVKB5E+h+W6HgMTIWLGwKD3XRwMEYXrqKy6DiuuOgV6bodySjy+0b6zz2fBlfhFmfjTrZn0s9sKsruPEQGUuk/oVn6j4J9BGbhMNtwKz5veWQg3pDWjfTfjsaE1EdLJ01eHp6/q86DkefWgp8ekVk4NhAH/Q58/ZH+J+bjwaMdfGYTgRHbMrOkh9Rc/8NOc64HTlt6YWaPoMgM//Uh2xVi7HAKA1HV0kCkf9EYpP8E6V80hrrrf9g9uByYNLOnla/PAz6mSAKinhTPDFYtDUT6Fw1C+u9B+hcNot76z80z8CTgfOAs4E1m9nx3f5G7z5nZC4ELzKwN3AvsXKW+CrHmdDrrX/SDIv0LUSD9S/+iudRZ/1mTAXf/OvD0ivduBJ46wj4JMZYshRZrGtK/ENJ/xXvSv2gEddf/akUTEsMwGZiKIt9LmFl0CCrMOVFmwchs6FEG42Em0qGBKK7aWkgdTb4YpE4ckmENRGY2Q5Gl83QK/b3Z3a+oqHsu8FqKxxW/BLzevehBv3bM7ATgbcCPAhuBNvAad//mcL0XoybMQDyZ6rAzGQQQiMaEKA1wRHRbe4Bj24Oxx6PAANEYFZVFV9yCfbGK2/HWidoMsp8PKWDpvyFERvMo23BUFhnXg6ziAJ0NUQbiQP9BNt5O0GRkDA6DfURarQogEB3zkVk4s14r/NUZnOvn0lrtiqzEHo2PrfRz9CA4yyDUWf+aDAiRiWN0hr8ysAswd3+KmT0GuM7MHuvuP+yuZGZnUtyaPxPYA1wJvBp4d0Y7ZwKz7v7Msq23AX8OnDFs54VoKtK/EM2l7vqv7z0PIVYB77Msh5m1gFcA7wFw9+8CNwIvC6q/HPiMu+929w7wXuBVme3cCLypq62vASdn76QQIkT6F6K51Fn/ujMgRC4OPpyB6BRgO3BLV9nNwNlB3XOAT/bUe5yZbQRO6teOu9+zVFjW/zXg4mE6LkTjkf6FaC4117/uDAgxAO5WuWRwfPl3b1fZfcBxFXV76xmwI7cdM3sDcAcwBbw5p4NCiGqkfyGaS531r8mAEJk4RWixqgXYYWY3dC2vzGx6kMsN/eoe9p67/xFwLPBt4CozGyL9ohDNRvoXornUXf96TGiNCVOXBynK29N5EX2G6kvVLbAo8MgwkUNyaacbntwft9d6IAgztLCQlg0TDmD524S73T265bfEXeXfbV3/b+36v7futq7XW4secDewlHR+2Xbc3c3s7RTPEP408IV+OyBWh1DnEGrdN6RhQjrT6aphhLHVCIMdtRldRsoti9qLgpsE+m/NV+h/LoomlpbRDspykf7XF1FEoOx10wM30rAF+g0jB1VEE4q03t6Qbrsd6T+KJhSMCWHkr+i8PMjHFZ1Ho6vjQTShoIhOcKruTOT/5ok2bWEfdf6v6pzuDAgxCMM4iOBWisgAp3WVnQFcH9S9Pqh3k7sfXK4dM3uemT3swS67t4FDwKasXgohYqR/IZpLjfWvyYAQ2RjeqV6Wo4wKcDFwHoCZnUqR1fMyM9thZtea2TFl9UuA55rZ9jJ6wE6KuMJ92ynXPZuuTKBm9ssUWv/qEDsvRMOR/oVoLvXWvx4TEiIXJ9co1I9dwIVmdh2F/l7q7nea2SMoEojMAnvc/Ttm9nrgCoo7q9cAFyzXTvneR4C3mtk1FINAC3hB1/tCiEGR/oVoLjXXvyYDQgzCkIOBux+ia9beVX4Hhdmnu+wyHprtZ7VTvvcd4MVDdVQIkSL9C9Fcaqx/TQaOFBVmJgvSY7e3HZWUHdq+Gs7AwwlNwRAbiCMDUtX6K+5QWjS1P67a2pe+sTgfuJKGZcS7KGpIoHWbjIdam51Nyha2pG7BhU1pmz4VWPFyhwkLDuSKh0azz3/RmNDKFEywEQt2b2KuYtMHA60vpEEFvDOkgKX/I0euATgw+xbFK18/Oi8TGYhnNqRlGzcmZT47E256cVPqAl6YDX4TzKT7EhmIw3N4rrF/EKLfBGEUgKAoCAwQjlsD/OSJhrPILGzSfyWaDAiRiwPDJR0RQqxXpH8hmkvN9a/JgBAD4FFcNCFEI5D+hWgudda/JgNCDMLwBiIhxHpF+heiudRY/5oMCJGLx88xCyEagPQvRHOpuf41GVgNIgPhdJBCEGhtPyYp23dymhvi4IlH4CiMDImAR5kAc42BQxAZjSbn4u36ocBZOPJ7elbrKwNiGSJdR6bCQOu2OQ0KANA+bmtSduD41Bk4n1ajMx0ZiDN1mZtVGMh2zY14TLAgWXBrMd6GLQQZiBcW04pDjQnS/8jIMAeHWbtzzb4ArUwTcK5ZeCpw7AYGYt+UGogXt6VlAPNb059gUbCAxWD1znSUWTit57m6rBo7cg/5aPXge44kNLSsco3Bw2Qgrrn+NRkQYhBqfGVACLEM0r8QzaXG+tdkQIhBqHFoMSHEMkj/QjSXGutfkwEhcnGwGocWE0L0QfoXornUXP+aDAgxCDW+MiCEWAbpX4jmUmP9azIwCGFm0dRU1NqYZhu0owMHILD/jBOSst1npl9Le1uaTfOIcQTMwrlU+neGzSyYSa4/U6wSkfGwIgNpuHqUlTRaP6gXZREOM5AelQYAaB+3LezPAyenGYgfOCntz/zW4GHVKNvoMBeuxkjnwHD7skpI/yPAAnNwrjE4MPHadGDsBYj0GmUCn0xdtz6V1vNgO52NaVn7qDSAwPyWuI9zW9Ntz2/ONBBHGYijLOBhdt8Rjx2V28nd9pAbz802PSR11r8mA0LkUvMMhEKIPkj/QjSXmutfkwEhBqHGVwaEEMsg/QvRXGqsf00GhBiAOicdEUL0R/oXornUWf+aDAgxCDW+MiCEWAbpX4jmUmP9r7/JQKaJF4hNgJFRJTdT4Ya8zKKLx25JyvadkpoKAfY8Nt323MPXyCxcdaCv0WNyHhgkF2dis2hk5AyNoR6kNc3Eah5abOxopRpsBSa+quzeBCbAqM1co6JvTI+xzuY0WMD89tTttz/IKgyw/2Hp8XTo2PTyU2c2uCQVGX4jDQ+bNXM1zIYZ24iyqXam4g37ZGRADRoYAul/VFhyfs0250daj8Z+wGfSup3pwBg8lR4nnemgbENatrgxKJsNzumb4/PWQmQWDn4qtGdSfXSmAl3mxlIIx4QhMxCPmgEyFYflYbCJle9M3fW//iYDQqwhdb5NKIToj/QvRHOps/41GRBiEGp8m1AIsQzSvxDNpcb612RAiFy83lcGhBB9kP6FaC41178mA0IMQo2vDAghlkH6F6K51Fj/Yz0ZiExFrSC7px29LVzfZ1Nzn0+kBhAPjGadjem257elhqSDO9J6+09Mt3Hw+HhK2dk6YrNwdLAuBOa6wAgTZi8EiMqPRLbSqfQzWzgqPmQjc2eUbdaHnNnXOQPhmhKZhTel2XlbW1Nzfmfb5rDJTpAJtBMYTX0qLWsHZsGFTYExcEuQLXhbetzNbQu7yMKWwCy8MdMsHJnZcg3EUb0qb1xkSoz6M4y3Lli3E/h/27FXlM6GdFyYiAzEA2SrjpD+h8daRmtDzxcZBQbYmBrxfXM6JnQ2xQdFezY9JhZnIrNwoP/p9IAMyzbkZQtenI3FEdVtbwzMwtN5BvuhNFgZaGD1tR4OUVH8jwr5enSuD8zCHgWLGIA663+sJwNCjB01HgyEEMsg/QvRXGqs/+GmSUI0ifKZwaolBzObMbNLzew6M7vBzJ7dp+65Zvb1st751hUXt187Znaymb3HzL5Yvv9BMzt2mF0XovFI/0I0l5rrX5MBITIxhh8MgF2AuftTgF8BPmRmxyfbMjsTOB/4WeDJwBOBV2e28+vAfe7+DOCpQBu4aKCdFUIchvQvRHOpu/41GRBiELzPsgxm1gJeAbwHwN2/C9wIvCyo/nLgM+6+2907wHuBV2W28w3g3eV7DnwIeNaguyqE6EH6F6K51Fj/4+MZCMwekVm4ffrJSdnex8TZfQ9tT9vsBHscGXEio1pk9mtvCjLaBsbXoWkHDpvFIMNyUC8qC6k4oMPiI2EgDlhIEz4D4LPBFxYZCBcXV77x4UOLnQJsB27pKrsZODuoew7wyZ56jzOzjcBJ/dpx9+71AGaA3UP1fNT06D3KLNzatjUpmz95R1K2/6Q0UADA3NY8w1+Y6TZIdBqNCYuzqQ7aG9KDxDfEeqk07SeNBrqODMTR8ZlrFq4y50WFYcbPIcaEyEAYZCBvb4g7GRqIo6y2gdGQ3KTk0v9osBb0GIijbMORWXhxW+q4nd8aZyBf2BxkB56JjMHpupFZuBMkEY/Gidx6AO1gXMj9jRLqLfL6RsFCBsoqPkC24hGSayoG4rErCBYzTAbiuutfdwaEGIQhrgwAS7fx9naV3QccV1G3t54BOwZsB+B5lFcKhBBDIP0L0VxqrP/xuTMgxDpgmSsDO8zshq7XF7l7zrO6g1yu6Fc3ec/Mfhx4NPAbA2xDCBEg/QvRXOqsf00GhMhl+SsAu909uuW3xF3l321d/2/t+r+37rau11vLrd8NLN1T79uOmZ0MvBP4N+6+0LfnQoj+SP9CNJea63/Zx4TMbNLMftvMrjazL5nZV8zsmV3vn1WWfdnMPmVm25drU4j1ypDRBG4F9gCndZWdAVwf1L0+qHeTux/MacfMdgAfAHa6+w/N7JFmVpGyqRrpX4iHkP6lf9Fc6qz/nDsDDwd+BzjL3e8zs2cBnzCz0yhmKZ8oN3iVmb0NuBB4cUa7h2HTgcPmxPQRqN0/lpqF7318/E34bK4zbA2ZC7IDBxmDW/ORgTBtLjTYBGah0JBUcUBHBuQoux+B4W/UVGYgnQmyVUcZCIfc/jAeSXfvmNnFwHnAtWZ2KnAWcG4p3o8DL3D3PcAlwOfLk+u9wE4KbfVtB8DMjgI+Cvw+cFf5+lUU4cVuG7Dbo9e/gfWYuy3INtw+LjUQ7zs5NQvvOzm+pjG/NTDnReb+MNNl8EUHmsk2AFcQmvuDYSusl7vp0OwbNOdxgxZUzt7rIQINRGNZZK6EOIusTY5XBmLpv8QM23D4+d5ngmzhQWbhyCw8d3R0MoP5zUF24CATcHROCU28QVkn0H9UrzJzbrR+NM7kHraZx6dF7tyKB1DCsXCoFMRDNDdIBuIoWMAaZiAed/3nTAb2Ab/v7veVHfm8mR0CfhyYB9ruflVZ9xLgNjM71t3vzmhbiPWDUzlhGoBdwIVmdh2F/l7q7nea2SOA04FZYI+7f8fMXg9cUW71GuCC5dop3/td4KeAqziclcQal/6FAOkf6V80mJrrf9nJgLvfQ3G7AQAzM2Ca4qrAz9AV3sjd7zCzAxQJEj63XNtCrCeM4a4MALj7IYpZfm/5HcCxPWWXAZcN0k753luAtwzX0wfbkv6FQPoH6V80l7rrfyX3TJ4G3E4xU+kNfwQVIY7M7JVWpE2+YYG5FWxWiLXHvHppCMPr36V/sT6R/ofX/3zn4Kp3UojVoM76H2gyYGYzwDuAXy+zolVW7S1w94vc/Wx3P3uKgX1MQowHnT5LzRmZ/gf3MQoxHkj/Q+t/upUmDhNiXVBj/WdPBsrbgxcB/9Xdl2Kp9oY/gupQSUKsb/pcFajDlYF+SP+i8Uj/0r9oLjXX/yB5Bv4Y+Jq7/0UZouh4ilBGL1mqUJogZoGvL9taT6QXC1LHLx6TRg46eFzgGh+3qEGLQeSfQ3HEgyhKUGsxqBgdbFEUlGjVzClflLYcIAoyYsHGPZoejzrCUMW+dKbSz7c1EX/mwzBkOvL1zAj1b9AbTSiIJra4Ob2DMHd0egDMb4uPscUtwbgQRcaIDvvcwzaqF+jIgjEBwBbyooQRRv8IIpGE0cTCTedtA6CT+Znlbij8DvLKqqOyRJFDAv1HEUYGQPofgf7NoDfS01R6/m9vSL+/9sb0AFgIIgQBLBwVRRNK67U3rDwiUBhpJ4rUVXUOHnFQnqHD5QVEvws81Gtmf6LfGavxYzrSehQFcQDqrP+sn4lm9kZgCri0DFP0aIqwRpcDk2b2tLLqecDHFElA1Bbvs9QU6V+IEulf+hfNpcb6X/bOgJk9hiKLGcBrut56m7vPmdkLgQvMrM1D8VCFqB9e7ysDEdK/ECXSv/QvmkvN9Z8TWvS79LmZ5e43Ak8dZaeEGEeMeg8GEdK/EAXSf/i+9C8aQd31P4hnQAhRg9uBQogVIv0L0VxqrP+1mwz0pIW3wOzVmQpSzEde4XbFhYsg1Xc2UZuRMTioF5kCW0EZVLjQg6qdqbRiZzpYOUplHmwk7HeVyTHYjEefRWQqjqbSQ3wv4fdfRWu41OMJDhaZKcVAGGC9Rq5A/2GK+VEb7qDCdR+V5emf4Bit1FbmlSYPjPge6SjXQByalKs2HjSZayoc4rMd5MteleMi2Yj0PzJ69R8ZOwNTeGzijTfhwbmwExmDg7JOpLdcba3CsRhLK3NDwx6y0e+HXCNuNHZEcR2iMTMao6rGy9CoHH1oQ3wYNde/7gwIMQB1CCEmhFgZ0r8QzaXO+tdkQIgBqPMzg0KI/kj/QjSXOutfkwEhBqHGVwaEEMsg/QvRXGqsf00GhMil5qHFhBB9kP6FaC411/9YTwYm9y8kZRv2pllJ5/fGu7G4OdNtGia/jFy8Qb3ISJtpCgboBCbAzlRazzcEG5/KdR8GRQuZ6YupyEyc+/BcuG5QL8rkmGs+rCA0OUVlmW0WocVqfGngCOGA9xq5FtO02xMH07LpfWmm4vl9FZmzJ9JjvBNlG801tA9zPFaZHCMZBv3JNwuvXJehGRrifczOahyV5WkwOvFGRkMAa6++LqX/UeHQ6flye18DtpCWtRbSz7+V/kyoLI/KIqNxbJDNzDY8gIF4mGfQh/IpD7TdzHN4ptZb6bBOay6tNzGX1puYD7YBtIJjxdrBABKVZVJ3/Y/1ZECIcaPOBiIhRH+kfyGaS531r8mAELn4gKFNhRD1QfoXornUXP+aDAgxCDW+MiCEWAbpX4jmUmP9azIgRC41TzoihOiD9C9Ec6m5/tduMtA5/H5L5+ChpMrE7vuTsi23z6RNTQWOW+DgcakzaPGoKJNvairxycgtnBZ56meMqXIVhcbAER9w0aajfQ4yvkLFZDg38+Iw+xKYHKsMRKFZKGKYDITU+5nBI4dD+3D9+6FU/5N79idlsz9MBdeZiocxW0wdtgtb0i+wPRMYdoOM3x4cy7nGXq905+eW5ZkXQ3KP2WEzp0bbyc3QHNSbOBiVxTszMReZBUd/T1/6HwHuMHf4QG5BtviJA+l5feqBVOvt6fwAAtFBHh17nTROSZy9eNisxMMcT2HgkyG2UTVEDWNyDj5bCwzEk+nwz9T+dMNTB+Pz/OShVOs2F5QtBBsfgDrrX3cGhMikiCaw1r0QQqwF0r8QzaXu+tdkQIhc3Gt9m1AI0QfpX4jmUnP9azIgxCDUdywQQiyH9C9Ec6mx/itS4QghErxIblS15GBmM2Z2qZldZ2Y3mNmz+9Q918y+XtY737oy4SzXjhX8BzOrcFkIIQZC+heiudRc/2NzZ8AX0j53fvDDpGw2MAVNHjw6bHPfI1IX0P4T0/Xnjk5NLu3NgcE2MN2GBuCQMZtSzqefQ5hpmNg4GRqiIs91J5hvBqaiKDPk9H3pukfdUWEguvdAUubzq3AeHP5r3AWYuz/FzB4DXGdmj3X3ww52MzsTOB84E9gDXAm8Gnh3ZjvvA/4JiN31a4mD95g7O/sPJtVa99yblEX6n5jfHG5man+66we3p+vPb02PxyjQQHs6MhAHGx42AEBoAoz0mlYLyyITb1RvkG5n6z+vP61AqhvuS8s23RWbgifvS48fPxSkMB3WVCz9D0/H8V4DcVDNJlNxTU2mOqiKzWGB47e1EIwfs+m67Q2BqTgwKnci/a+GET/aTKZZOHzOPdd8XEXm+lFc/iiL9ERgIJ5+IO349P2xfif2pz8gLNL/fEW66lxqrH/dGRBiAKzjlcuy65q1gFcA7wFw9+8CNwIvC6q/HPiMu+929w7wXuBVA7TzZuCSFe6mECJA+heiudRZ/5oMCDEA5tVLBqcA24FbuspuBs4O6p4T1HucmW3Macfdv5fVIyFENtK/EM2lzvrXZECITMyXvTKwo3x+b2l5ZU8Tx5d/93aV3QccF2zu+KCeATsGbEcIMQKkfyGaS931PzaeASHWA8sYhXa7ezTLX7bZEdUd9mlVIUQfpH8hmkud9T/Wk4FOkJXUbv/npGz63sBpBmz/59RYvPmkrUnZvoenRuODx6bOoIWj0rIwe2l0v6XqHkxg7mlFppu5IBtn4I+JjDgTc2kfWwMk4osyL0ZEbUYZgyfm0/5MHgrMQnvTlad/kGalBvDvp2bzzlzwAQ2DM6yB6K7y77au/7d2/d9bd1vX663l1u8Glg7YnHbGk55M0L6YGrs6+x5IylpBpumZg/H3PLX3qLTujo1J2aHt6QE+tyUwFW9KRdwOMpCH+q8g23QXaCsy3UfZPVuL6UaibVTd6q4yaCbr545lC2nFyWCMmro/3ZkNd6eBAgBad+9NyjoHA1PxMAZi6X80uKfBHTw9JrqCpzxIJK3pih9oEwupOCcPplpfnE1bXZwJgopMpWWemZU4V0NQMSZEJuBcw/4QpuIqzIMxJbM/0Xg0cSjQ//5U/5P74qAgrftTrVtwXvBhfhPUXP96TEiIbKpvEWYmI7mVIjLAaV1lZwDXB3WvD+rd5O4HB2xHCDESpH8hmku99a/JgBCD4F69LLuqd4CLgfMAzOxU4CzgMjPbYWbXmtkxZfVLgOea2fYyesBO4MLl2hndjgohEqR/IZpLjfWvyYAQuYwg6QhlfGAzuw74IPBSd78T2AicDswCuPt3gNcDVwBfBb4JXJDRDhRvvAX4UPn/1YGZSQgxCNK/EM2l5vofa8+AEGPHsPmk3A9RzPJ7y+8Aju0pu4yK2X5VO13v/yHwh0N0VQjRi/QvRHOpsf7X3WQgylTcvmdPWNfuTw2I03endXfcnpqK28ek5sPFzakhaXEmNRX7RL5bKDTTHUpNbq25oOxgkHVvPnAQLqRlFhgxCQxbAB5kfQ3rdgITWLBtFtN98SgzYGD2qTIFh9mGM27dDYoF+yhGQPBdRd9pJzCAWsUx0bp/X1I2e09gKt66KSlb3JwGFVjclA6XnekgI2qQldQrtBWaeyOD3Xx63LWCMlsM6gVlkVYrqeh7QnB1LNKLLQQm3miMigyAgSkYoLM/NRZ3ojFlyDFB+h8ed0/H+8jYHT2HHXz+E9E5BmjNzaR196fn8PZMqut2cF6PtR6YigfQf0T0/Hlszo0+n8z2IhlUPPeem5k47nfe+NaaD8b14PdN60AQIQXgYFoeZSAPfycMQJ31v+4mA0KsFeYD3Q4UQtQI6V+I5lJ3/WsyIMQgrMLdBiHEOkH6F6K51Fj/mgwIkYsTPgohhGgA0r8QzaXm+tdkQIgBiJKtCCGagfQvRHOps/5rPRkIzcZ7A9PtvtRoaN9PP5qpqamkbDoy104EDqIKc44vBuan0EyVGlc8MAF7kMnR8xJiYK0KA3FkAqqom7XtoI/ZrKkYfTDjpRiOyFQc6KUqq6wFZjELjKa2NzULT21Iy6YD/TMVDKEDmAXD4ykw2Edjgkf1OplGzEE0aJkRqKM2g22H31fwXUdm8XBsHKDN4ZD+R4J7+n1F560oy230PUdGccACA+nEgVTXrelU15NBmU8F5/XJwFQcyb/qfJlt7o20lRaFBtdIB2Gm4gH0EtUNtxOURUEgorEsMIZXZhAOjoHILBz9Zsqn3vqv9WRAiJHi1PqZQSFEH6R/IZpLzfWvyYAQA1DnaAJCiP5I/0I0lzrrX5MBIXJxYKjbjEKIdYv0L0Rzqbn+NRkQIhuv9W1CIUQ/pH8hmku99a/JgBCDUGMDkRBiGaR/IZpLjfXfvMlAboSSKMrPoSAVdm7kkHUwoxwkwMgwAYHWLU5lVCixhlRoKzvyUBCJwlr703pRVJ0gSogNEE0oipgyTPSf3MhhR4x1GzksQPofHT1RrzwKwRN91pnRaSDWv82lEWaYTH8GtSajyEHBz6UoSlCk/6oxYZioPNnrZkYYGoBw3IoIv8MgMloUdSw30iIV43oUgW0Y/dZc/82bDAixYjwO3SiEaADSvxDNpd7612RAiFxqfmVACNEH6V+I5lJz/WdlkzGz15jZlWb2RTP7tpm9puu9s8zsK2b2ZTP7lJltX73uCrHGdDrVS02R/oUokf6lf9Fcaqz/3DsDrwCe5e53mdmpwN+a2VeBbwKfAHa6+1Vm9jbgQuDFq9JbIdYS98pnFmuO9C+E9C/9i+ZSc/3nTgZ+1d3vAnD3vzeze4FHAScCbXe/qqx3CXCbmR3r7nePvLfjyLiZ3MTq0szvuz76D013kdFsiE2sfFUx7kj/q6P/XF0GwvR2/ICDBTHhvZUaiG0iNQt7UBYGBggNxFkPXFQTDT6Zj6fEAQlW4ar1MDoI+jNsIIXQGLwaUU5qrP+so9bd/2bpfzP7ReAB4ArgHOCWrnp3AAeAJ462m0KMA14kHalaaor0LwRI/9K/aDL11n+2gdjMHg98GNgE/JK77zWz44G9PVXvA44L1n8l8EqAGWZX2l8h1g6Pr0w1AelfNB7pX/oXzaXm+s+eDLj7t4EzzOwJwOVm9vw+1ZN7Z+5+EXARwBY7pr73WkS9qcEVgJUg/QuB9C/9iyZTY/0P/HCbu98IfBr4LeAuYFtPla1luRD1wr3W0QRykP5FY5H+pX/RXGqu/2XvDJShwp7u7h/tKt4PHANcD7ykq+4jgFng6yPupxBjQZjpcADMbIYi4sbpFPp7s7tfUVH3XOC1FJ7ULwGv99JptVw7ZnYW8GdAG7iXIuLHPSvor/QvRIn0D6yl/kMDZ4WpNPqqwrKgMDAB+5DZxodhqIy/q8GIH5fJzgw8yHZXwexbZ/3n3BnYDPyemc2WGzkWeCFwJXA5MGlmTyvrngd8bGwjiQgxDD4SA9EuwNz9KcCvAB8qn709DDM7Ezgf+FngyRSmvFfntGNm0xQh/97s7j8BfINi4FgJ0r8QIP0j/YsGU3P950wG7gQ+BVxpZl8CvgC8F7jA3ecoBoZ3mtmXgbOBV2W0KcT6xDvVyzKYWYsiZvd7ANz9u8CNwMuC6i8HPuPuu71wLb2XUlsZ7fwcaci/F5Un8kGR/oVYQvqX/kVzqbH+l31MyN0PAb9XLtH7NwJPXa4dIdY77j7sbcJTgO10heMDbqY4ifZyDvDJnnqPM7ONwEnLtJOE/DOzpZB/nxukw9K/EAXSf/i+9C8aQd31P2R2DCGahXe8cslg6Xbg3q6yMBRfWbe3ngE7MtrJDvknhMhH+heiudRZ/9mhRUfJPu7d/QX/yO3lyx3A7rXoxyqgfRlPltuXk3Ma2ce9n/tC5//f0afKjJnd0PX6ojKk3nIM4kLrV3e5do6M220ZpP91QZP2Rfo/goxc/8P6RIfzhHbTJM2sJ6T/jPfXZDLg7g8+u2RmN7h7dJtk3aF9GU9GtS/u/pwhm1gKubet6/+qUHy9Yfu2Upz27gY2LNPOXcATetobm5B/0v/4o31Jkf5Hg/Q//mhfUuqufz0mJMSR41ZgD3BaV9kZFCH6erk+qHeTux/MaOewdRXyU4ixQPoXormMtf41GRDiCFFGBbiYIgQfZnYqcBZwmZntMLNrzeyYsvolwHPNbHsZPWAnZXiwfu2U6yrknxBjhvQvRHMZd/2vyWNCPeQ8U7Ve0L6MJ+O0L7uAC83sOgr9vdTd7yxn76dTzOD3uPt3zOz1wBUUWXWuAS5Yrh0Ad58zsxcCF5jZg0lHjsC+rYRx+m6GRfsynozTvuxC+u9mnL6bYdG+jCfjtC+7GFP9W3ZmOyGEEEIIIUSt0GNCQgghhBBCNJQ1mwyY2YyZXWpm15nZDWb27LXqy0ows3PM7FYz29lTfpaZfcXMvmxmnzKz7WvUxWUxs0kz+20zu9rMvlT2+5ld76+bfQEws9eY2ZVm9kUz+7aZvabrvXW1L3VH+l97pH+xVkj/a4/0Lw7D3ddkAd4JvK/8/zEU7ujj16o/A/b9F4A/p3Bn7+wqnwZuB366fP024C/Wur999uNRwD8CW8vXzwIeoMhwt672pezjt4Djyv9PBRaBJ6/Hfan7Iv2v/SL9r32fm7pI/2u/SP9r3+dxWtbqS2tRJIH4qa6yK4HXrfUHktn/h5d/r+4ZDF4A/EPX60dQpDQ5dq37XLEf24Ff7SnbDbx4ve1L2ccf63l9N/BL63Ff6rxI/+OxSP/juy91XqT/8Vik//Hdl7VY1uoxoVMoDsRbuspuBtZFkgt3/17FW+fQtU/ufgdwAHjikejXoLj7Pe7+gaXXZmYUs+i7WWf7AuDuf7P0v5n9IsVVjitYh/tSc6T/MUD6H999qTnS/xgg/Y/vvqwFazUZOL78u7er7D7guCPflZFyPIfvE6yv/Xoaxe20a1in+2Jmjzezm4H/D3iJu+9lne5LjZH+xxPpXxwJpP/xRPpvMOOQZ6AbW+sOrBJjv19mNgO8A/h1d+8UFwniqkeuV4Pj7t8GzjCzJwCXm9nz+1Qf631pIHX9PsZ+v6R/MQbU9fsY+/2S/sVa3Rm4q/y7ratsa1f5euUuDt8nWAf7Vd4evAj4r+5+Q1m8LvdlCXe/Efg08Fus832pIdL/GCH9iyOM9D9GSP8C1m4ycCtF9IDTusrOAK5fm+6MjOvp2icrssrNUkQdGGf+GPiau/+FmW0ws0eyzvbFirTdv9hTvB/YxDrblwYg/Y8X0r84kkj/44X0L9ZmMuDuHeBi4DwAMzsVOAu4bC36M0IuBybN7Gnl6/OAj7n73WvYp76Y2RuBKeBSMzsKeDRFv9fbvmwGfs/MZgHM7FjghRRRKtbbvtQa6X98kP7FkUb6Hx+kf7GElaGWjvyGi2fULgROp/AuvNndr1iTzgyImT0JOJ9iALsTuNndX1S+9wTgAorQVfdShB67Z4262hczewzwd8Fbb3P3XetsX2aAtwA/A8xT3Bb8OPD28hnIdbMvTUD6X3uk//HclyYg/a890v947stasWaTASGEEEIIIcTaslaeASGEEEIIIcQao8mAEEIIIYQQDUWTASGEEEIIIRqKJgNCCCGEEEI0FE0GhBBCCCGEaCiaDAghhBBCCNFQNBkQQgghhBCioWgyIIQQQgghREPRZEAIIYQQQoiG8n8BGXlfKgR7BZAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 936x216 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE =  0.00021775\n",
      "Chi^2 =  34.5336289\n"
     ]
    }
   ],
   "source": [
    "vis_split = np.column_stack((np.real(visibility), np.imag(visibility)))\n",
    "vis_split = np.float32(224.46*vis_split.reshape((1, 1691,2)))\n",
    "[pred_img, pred_vis] = vis_model.predict(vis_split, steps=1)\n",
    "\n",
    "# Normalize predicted image by constant \n",
    "# (used for training neural network to prevent dying weights during back propagation)\n",
    "ALPHA = 500\n",
    "pred_img /= ALPHA\n",
    "\n",
    "# Normalize predicted image flux to target image flux\n",
    "pred_img = post_process(pred_img, target)\n",
    "\n",
    "# Apply blur by 0.5*fwhm to target image\n",
    "target_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "target_blur.imvec = target.flatten()\n",
    "target_blur = target_blur.blur_circ(fwhm_i=0.5*fwhm, fwhm_pol=0.5*fwhm)\n",
    "target = target_blur.imvec\n",
    "\n",
    "# Apply blur by 0.5*fwhm to predicted image\n",
    "fwhm = obs.res()\n",
    "pred_blur = eh.image.make_empty(npix, fov, ra, dec, rf=rf, source='random', mjd=mjd)\n",
    "pred_blur.imvec = pred_img.flatten()\n",
    "pred_blur_img = pred_blur.blur_circ(fwhm_i=0.5*fwhm, fwhm_pol=0.5*fwhm)\n",
    "pred_blur = pred_blur_img.imvec\n",
    "\n",
    "mae = np.mean(np.abs((np.subtract(pred_blur.flatten(), target.flatten()))))\n",
    "chisq = obs.chisq(pred_blur_img,dtype='vis') \n",
    "\n",
    "visualize(target, pred_img, obs)\n",
    "\n",
    "print(\"MAE = \", np.round(mae, 8))\n",
    "print(\"Chi^2 = \", np.round(chisq, 8))"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}