08a-CDMS_NNVAE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "954ee364",
   "metadata": {},
   "source": [
    "## Introduction to Variational Auto Encoders (VAEs)\n",
    "\n",
    "This notebook will introduce generative models, VAEs in particular. Earlier, we used an MLP to predict the impact location on the CDMS prototype detector based on timing measurements. However, in many applications in HEP, we need to simulate the behavior of the detector and obtain simulated dataset to be matched against the data. In this notebook, our focus is on learning the joint distribution of the timing measurements and the impact location via a NN. \n",
    "\n",
    "A VAE consists of two networks. The first network is called an **Encoder**, that maps the input features $(\\mathbf{x})$ to a lower dimensional *latent* feature space $(\\mathbf{z})$. The idea of the Encoder network is to map the unknown, perhaps complicated and untractable distribution of the input feature space to a known and tractable distribution in the latent space. Generally, the preferred choice of distribution for the latent features space is a Gaussian distribution whose mean $(\\mu)$ and variances $(\\sigma^2)$ are predicted by the Encoder network.\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu, \\log\\sigma^2 = \\mathrm{EncoderNN}_\\theta(x)\n",
    "\\end{equation}\n",
    "\n",
    "here $\\theta$ represents the trainable parameters of the Encoder network.\n",
    "\n",
    "**Can you think of a reason why we calculate the logarithm of the variance from the Encoder network instead of the variance directly?**\n",
    "\n",
    "The second part of the VAE is, as you may intuitively guess, is called a **Decoder**. It is trained to predict a sample point in the input feature space from a sample in the latent space. The input dimension of the Decoder network is the same as the dimension of the latent space and the output is of the same dimension as the original input feature space.\n",
    "\n",
    "\\begin{equation}\n",
    "z \\sim \\mathcal{N}(\\mu, \\sigma^2) \\\\\n",
    "\\tilde{x} = \\mathrm{DecoderNN}_\\phi(z)\n",
    "\\end{equation}\n",
    "\n",
    "here $\\phi$ represents the set of trainable parameters of the Decoder network and $\\tilde{x}$ is the prediction of the Decoder network for a sample, $z$, taken from the latent space. \n",
    "\n",
    "To train a VAE means to have the distribution of the Decoder predicted sample, $\\tilde{x}$, be as close as to the input feature distribution, $x$. Simultaneously, we want the latent space distribution to be tractable and impervious to overfitting the inherent statistical fluctuations in the input feature distributions. The first part can be accomplised by minimizing the MSE loss between $x$ and $\\tilde{x}$.\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathcal{L}_{MSE} = \\frac{1}{N}\\sum_{i=0}^{N-1}\\left(x_i - \\tilde{x}_i\\right)^2\n",
    "\\end{equation}\n",
    "\n",
    "For the second part, we want to make the parameterized, marginal distribution of the latent space $p_\\theta(z|x)$ to be as close as possible to the prior distribution of $z$, given by $q(z) \\sim \\mathcal{N}(0,1)$. This is done by minimizing the so-called **KL Divergence** between these two distributions. The KL-divergence between two distributions is obtained by-\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{KL}(p_1(y)||p_2(y)) = \\mathbf{E}_{p_1}\\left(\\log{\\frac{p_1(y)}{p_2(y)}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "KL divergence is non-negative and is minimum when $p_1 = p2$. Also, it is not symmetric, i.e. $\\mathrm{KL}(p_1(y)||p_2(y)) \\neq \\mathrm{KL}(p_2(y)||p_1(y))$. Going back to the loss function for VAE training- it can be written as:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L}_{VAE} &= \\mathcal{L}_{MSE} + \\mathcal{L}_{DL} \\\\\n",
    "& = \\frac{1}{N}\\sum_{i=0}^{N-1}\\left(x_i - \\tilde{x}_i\\right)^2 + \\mathrm{KL}(p_\\theta(z|x)||q(z))\n",
    "\\end{align}\n",
    "\n",
    "It can be shown that for $p_\\theta(z|x) \\sim \\mathcal{N}(\\mu, \\sigma^2)$ and $q(z) \\sim \\mathcal{N}(0,1)$, the KL-diveregence takes the following form:\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{KL}(p_\\theta(z|x)||q(z)) = -\\frac{1}{2N} \\sum_{i=0}^{N-1} \\left[ 1 + \\log\\sigma_i^2 -\\sigma_i^2 -\\mu_i^2 \\right]\n",
    "\\end{equation}\n",
    "\n",
    "A full derivation of the mathematical derivation of the VAE loss function can be found [in this paper](https://doi.org/10.48550/arXiv.1907.08956).\n",
    "\n",
    "In the following sections, we setup the code for training the VAE. Let's start by loading the necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "id": "31da1c8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.model_selection import *\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch import optim\n",
    "from torch.utils.data import DataLoader\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "297fa8ab",
   "metadata": {},
   "source": [
    "Now we will get the data. Unlike our previous exercises, we will load the all the data features into the input vector for the Enoder network. So our input dimension for the Encoder network will be 20. We also want to make sure that the dataset has equal representation of all impact locations, so a uniform distribution can be assumed for the impact locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "id": "2a8cf63a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Get the dataset as numpy array\n",
    "\n",
    "data = pd.read_csv('data/CDMS_Dataset.csv', sep=',')\n",
    "\n",
    "features = list(data.columns)[1:]\n",
    "xy_data = data.values[:,1:]\n",
    "# n_features = len(features) - 1\n",
    "\n",
    "# unique_ys = np.unique(xy_data[:,-1])\n",
    "# counts = []\n",
    "\n",
    "# for y in unique_ys:\n",
    "#     counts.append(int(np.sum(xy_data[:,-1] == y)))\n",
    "\n",
    "# min_count = min(counts)\n",
    "\n",
    "# xy_data_uniformized = np.zeros((min_count*len(unique_ys), n_features+1))\n",
    "# for ii,y in enumerate(unique_ys):\n",
    "#     indices = xy_data[:,-1] == y\n",
    "#     selected = xy_data[indices][:min_count, :]\n",
    "#     xy_data_uniformized[ii*min_count:(ii+1)*min_count, :] = selected\n",
    "\n",
    "# print(xy_data_uniformized)\n",
    "# xy_data = xy_data_uniformized\n",
    "y_norm = -41.9 # largest y dimension to be used to normalize y data\n",
    "x_mean = np.mean(xy_data[:,:-1], axis = 0).reshape(1, -1)\n",
    "x_std  = np.std(xy_data[:,:-1], axis = 0).reshape(1, -1)\n",
    "\n",
    "\n",
    "def split_data(split_type='random', labels = []):\n",
    "    if split_type == 'random':\n",
    "        XY_train, XY_test  = train_test_split(xy_data, test_size=0.25, stratify=xy_data[:, -1])\n",
    "    if split_type == 'label':\n",
    "        XY_train = data.loc[~data[features[-1]].isin(labels)].values[:, 1:]\n",
    "        XY_test = data.loc[data[features[-1]].isin(labels)].values[:, 1:]\n",
    "    XY_train[:,:-1] = (XY_train[:,:-1] - x_mean)/x_std\n",
    "    XY_test[:,:-1]  = (XY_test[:,:-1] - x_mean)/x_std\n",
    "    XY_train[:,-1]  = 1*(2*XY_train[:,-1]/y_norm - 1)\n",
    "    XY_test[:,-1]   = 1*(2*XY_test[:,-1]/y_norm -1)\n",
    "    \n",
    "    return XY_train, XY_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "id": "e30e54e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Defining a class for a generic MLP\n",
    "\n",
    "class Network(nn.Module):\n",
    "    def __init__(self, Layers, device, activation, \n",
    "                 end_in_sigmoid = False, \n",
    "                 add_dropout = False, \n",
    "                 add_batchnorm = False):\n",
    "        super(Network, self).__init__()\n",
    "        self.Layers = Layers\n",
    "        self.device = device\n",
    "        self.activation = activation\n",
    "        self.end_in_sigmoid = end_in_sigmoid\n",
    "        self.add_dropout = add_dropout,\n",
    "        self.add_batchnorm = add_batchnorm,\n",
    "        self.NN = self.build_model().to(device)\n",
    "        \n",
    "    def build_model(self):\n",
    "        Seq = nn.Sequential()\n",
    "        for ii in range(len(self.Layers)-1):\n",
    "            this_module = nn.Linear(self.Layers[ii], self.Layers[ii+1])\n",
    "            Seq.add_module(\"Linear\" + str(ii), this_module)\n",
    "            if not (ii == len(self.Layers)-2):\n",
    "                if self.add_batchnorm[0]:\n",
    "                    Seq.add_module(\"BatchNorm\" + str(ii), nn.BatchNorm1d(self.Layers[ii+1]))\n",
    "                Seq.add_module(\"Activation\" + str(ii), self.activation)\n",
    "                if self.add_dropout[0]:\n",
    "                    Seq.add_module(\"Dropout\" +str(ii), nn.Dropout(p=0.5, inplace=False))\n",
    "            if self.end_in_sigmoid:\n",
    "                Seq.add_module(\"Sigmoid\" + str(ii), nn.Sigmoid())\n",
    "        return Seq\n",
    "    \n",
    "    def forward(self, X):\n",
    "        X = X.to(self.device)\n",
    "        return self.NN(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "id": "ecee55b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "class CDMS_NNSuite:\n",
    "    def __init__(self, XY_train, XY_test, device):\n",
    "        self.XY_train = torch.tensor(XY_train, dtype = torch.float32, device = device)\n",
    "        self.XY_test  = torch.tensor(XY_test,  dtype = torch.float32, device = device)\n",
    "        self.device = device\n",
    "\n",
    "    \n",
    "class CDMS_VAE(CDMS_NNSuite):\n",
    "    def __init__(self, XY_train, XY_test, \n",
    "                 enc_Layers, dec_Layers,\n",
    "                 device, activation=nn.Sigmoid()):\n",
    "        super().__init__(XY_train, XY_test, device)\n",
    "        self.enc_Layers = enc_Layers\n",
    "        self.dec_Layers = dec_Layers\n",
    "        self.Encoder = Network(enc_Layers, device, activation,\n",
    "                              end_in_sigmoid = False, \n",
    "                              add_dropout = False, \n",
    "                              add_batchnorm = False)\n",
    "        self.Decoder = Network(dec_Layers, device, activation)\n",
    "        \n",
    "    \n",
    "    def Evaluator(self, xy):\n",
    "        enc_output = self.Encoder.forward(xy).to(self.device)\n",
    "        mu, logvar = enc_output[:, :int(self.enc_Layers[-1]/2)], enc_output[:, int(self.enc_Layers[-1]/2):]\n",
    "        #print(mu.shape, logvar.shape)\n",
    "        z = mu + torch.exp(logvar/2)*torch.randn_like(mu) # reparameterization\n",
    "        z = z.reshape((-1, int(self.enc_Layers[-1]/2)))\n",
    "        #print(z.shape)\n",
    "        dec_output = self.Decoder(z).to(self.device)\n",
    "                \n",
    "        mse_loss = 0.5* torch.sum((dec_output.reshape(-1) - xy.reshape(-1))**2)/(xy.shape[0])\n",
    "        KL_div = -0.5 * torch.sum(1 + logvar - mu**2 - torch.exp(logvar))/(xy.shape[0])\n",
    "        \n",
    "        return mu, logvar, z, dec_output, mse_loss, KL_div\n",
    "    \n",
    "    \n",
    "    def Train(self, n_epochs, batch_size):\n",
    "        Train_MSE_Losses = []\n",
    "        Train_KL_Losses = []\n",
    "        Test_MSE_Losses = []\n",
    "        Test_KL_Losses = []\n",
    "        params = list(self.Encoder.NN.parameters()) + list(self.Decoder.NN.parameters())\n",
    "        optimizer = optim.Adam(params, lr=1e-3)\n",
    "        for ii in range(n_epochs):\n",
    "            self.Encoder.train()\n",
    "            self.Decoder.train()\n",
    "            DL = DataLoader(self.XY_train, batch_size=batch_size, shuffle=True, drop_last=True)\n",
    "            for xy in DL:\n",
    "                mu, logvar, z, dec_output, mse_loss, KL_div = self.Evaluator(xy)\n",
    "                optimizer.zero_grad()\n",
    "                total_loss = mse_loss + KL_div\n",
    "                \n",
    "                total_loss.backward()\n",
    "                optimizer.step()\n",
    "                \n",
    "            self.Encoder.eval()    \n",
    "            self.Decoder.eval()\n",
    "            with torch.no_grad():\n",
    "                mu, logvar, z, dec_output, mse_loss, KL_div = self.Evaluator(self.XY_train)\n",
    "                Train_MSE_Losses.append(mse_loss.item())\n",
    "                Train_KL_Losses.append(KL_div.item())\n",
    "                \n",
    "                mu, logvar, z, dec_output, mse_loss, KL_div = self.Evaluator(self.XY_test)\n",
    "                Test_MSE_Losses.append(mse_loss.item())\n",
    "                Test_KL_Losses.append(KL_div.item())\n",
    "                \n",
    "                if ii % int(n_epochs/10) == 0 or ii == n_epochs - 1:\n",
    "                    print('''Epoch {0}/{1}, \n",
    "                        Training MSE Loss = {2:0.4f} \n",
    "                        Training KL Div   = {3:0.4f}\n",
    "                        Test MSE Loss     = {4:0.4f}\n",
    "                        Test KL div       = {5:0.4f}\n",
    "                        '''.format(ii, n_epochs, \n",
    "                                   Train_MSE_Losses[-1], Train_KL_Losses[-1],\n",
    "                                   Test_MSE_Losses[-1],  Test_KL_Losses[-1]\n",
    "                                  )\n",
    "                         )\n",
    "            \n",
    "        \n",
    "        return Train_MSE_Losses, Train_KL_Losses, Test_MSE_Losses, Test_KL_Losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "id": "6ad83289",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0/4000, \n",
      "                        Training MSE Loss = 9.7111 \n",
      "                        Training KL Div   = 0.0027\n",
      "                        Test MSE Loss     = 9.8097\n",
      "                        Test KL div       = 0.0028\n",
      "                        \n",
      "Epoch 400/4000, \n",
      "                        Training MSE Loss = 1.9707 \n",
      "                        Training KL Div   = 1.9067\n",
      "                        Test MSE Loss     = 2.0060\n",
      "                        Test KL div       = 1.9093\n",
      "                        \n",
      "Epoch 800/4000, \n",
      "                        Training MSE Loss = 1.8418 \n",
      "                        Training KL Div   = 1.9700\n",
      "                        Test MSE Loss     = 1.8860\n",
      "                        Test KL div       = 1.9730\n",
      "                        \n",
      "Epoch 1200/4000, \n",
      "                        Training MSE Loss = 1.7817 \n",
      "                        Training KL Div   = 2.0245\n",
      "                        Test MSE Loss     = 1.8147\n",
      "                        Test KL div       = 2.0259\n",
      "                        \n",
      "Epoch 1600/4000, \n",
      "                        Training MSE Loss = 1.7918 \n",
      "                        Training KL Div   = 2.0305\n",
      "                        Test MSE Loss     = 1.9061\n",
      "                        Test KL div       = 2.0301\n",
      "                        \n",
      "Epoch 2000/4000, \n",
      "                        Training MSE Loss = 1.8048 \n",
      "                        Training KL Div   = 1.9494\n",
      "                        Test MSE Loss     = 1.7866\n",
      "                        Test KL div       = 1.9527\n",
      "                        \n",
      "Epoch 2400/4000, \n",
      "                        Training MSE Loss = 1.7348 \n",
      "                        Training KL Div   = 1.9933\n",
      "                        Test MSE Loss     = 1.8110\n",
      "                        Test KL div       = 1.9973\n",
      "                        \n",
      "Epoch 2800/4000, \n",
      "                        Training MSE Loss = 1.6835 \n",
      "                        Training KL Div   = 2.0045\n",
      "                        Test MSE Loss     = 1.7692\n",
      "                        Test KL div       = 2.0081\n",
      "                        \n",
      "Epoch 3200/4000, \n",
      "                        Training MSE Loss = 1.6915 \n",
      "                        Training KL Div   = 1.9901\n",
      "                        Test MSE Loss     = 1.7720\n",
      "                        Test KL div       = 1.9930\n",
      "                        \n",
      "Epoch 3600/4000, \n",
      "                        Training MSE Loss = 1.7012 \n",
      "                        Training KL Div   = 1.9714\n",
      "                        Test MSE Loss     = 1.7845\n",
      "                        Test KL div       = 1.9740\n",
      "                        \n",
      "Epoch 3999/4000, \n",
      "                        Training MSE Loss = 1.6840 \n",
      "                        Training KL Div   = 1.9943\n",
      "                        Test MSE Loss     = 1.7865\n",
      "                        Test KL div       = 1.9952\n",
      "                        \n"
     ]
    }
   ],
   "source": [
    "NLatent = 4\n",
    "NNodes_enc = 32\n",
    "NNodes_dec = 128\n",
    "NhiddenLayers_enc = 2\n",
    "NhiddenLayers_dec = 2\n",
    "enc_Layers =  [n_features+1] + NhiddenLayers_enc*[NNodes_enc] + [2*NLatent]\n",
    "dec_Layers =  [NLatent]  + NhiddenLayers_dec*[NNodes_dec] + [n_features+1]\n",
    "# gen_Layers = [3, 32, 32, 20]\n",
    "# disc_Layers = [20, 128, 64, 1]\n",
    "XY_train, XY_test = split_data('random')\n",
    "lossfn = nn.MSELoss(reduction='mean')\n",
    "device = torch.device('cpu')\n",
    "#activation = nn.LeakyReLU(negative_slope=0.01)\n",
    "activation = nn.Sigmoid()\n",
    "n_epochs = 4000\n",
    "batch_size = 128\n",
    "\n",
    "\n",
    "NNvae =  CDMS_VAE(XY_train = XY_train, \n",
    "                  XY_test = XY_test, \n",
    "                  enc_Layers = enc_Layers, \n",
    "                  dec_Layers = dec_Layers,\n",
    "                  device = device, \n",
    "                  activation = activation)\n",
    "\n",
    "Train_MSE_Losses, Train_KL_Losses, Test_MSE_Losses, Test_KL_Losses = NNvae.Train(n_epochs = n_epochs,\n",
    "                                                                                 batch_size = batch_size\n",
    "                                                                                 )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "id": "6ba8c712",
   "metadata": {},
   "outputs": [],
   "source": [
    "with torch.no_grad():\n",
    "    mu, logvar, z, dec_output, mse_loss, KL_div = NNvae.Evaluator(NNvae.XY_train)\n",
    "    mu2, logvar2, z2, dec_output2, mse_loss2, KL_div2 = NNvae.Evaluator(NNvae.XY_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44e5c0f9",
   "metadata": {},
   "source": [
    "Now let's explore the behavior of the model. First, we want to see how well the model behaves in reproducing the impact location from the decoder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "id": "8a1f61a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa233523090>"
      ]
     },
     "execution_count": 278,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABthklEQVR4nO2dd3xUVfbAv+fNpJKELiJVWFACShEQKTbsYmMtKLsruAu6Yl8rKAjYVndtuNh2BWwruj8rimBFqgoIUqJBIPTekkCSKe/+/njzJjOTSTJJZjKTcL+fz3wy995X7tzM3PPuOeeeI0opNBqNRqMJxIh3BzQajUaTeGjhoNFoNJoyaOGg0Wg0mjJo4aDRaDSaMmjhoNFoNJoyOOPdgWjQrFkz1b59+3h3Q6PRaOoUy5Yt26uUah6urV4Ih/bt27N06dJ4d0Oj0WjqFCKyqbw2rVbSaDQaTRm0cNBoNBpNGbRw0Gg0Gk0Z6oXNIRxut5utW7dSXFwc764cVaSmptK6dWuSkpLi3RWNRlMD6q1w2Lp1K5mZmbRv3x4RiXd3jgqUUuzbt4+tW7dy/PHHx7s7Gk29xjRNDMMot1xT6q1aqbi4mKZNm2rBUIuICE2bNtWrNY0mxlzz8mIunrIA0zQBSzBcPGUB17y8OGr3qLfCAdCCIQ7oMddoYotpmqzfU0jOjgK/gLh4ygJydhSwfk+hX2DUlHotHDQaTd0kdIKL1oQXK9xud4XlaCIi7C10AZC2cxn/eOgm0nYuA2BvoStqD2haOMSIffv20aNHD3r06MGxxx5Lq1at/GWXy1XhuUuXLuW2226r9B79+/ePSl+PHDnC8OHDOemkk+jWrRsDBw6ksLCwwnMee+yxqNxbownlmpcXMyREZTIkyiqTaNJ53Gec8NBcv0Bwu92c8NBcOo/7LCb3s+/TS3J5J3ky9zhn8k7yZHpJblB7TdHCwUdo0qOaJkFq2rQpK1asYMWKFdx0003ceeed/nJycjIej6fcc3v37s3zzz9f6T0WLVpUoz7aPPfcc7Ro0YJVq1axevVq/vOf/1TqbaSFgyYWmKZJQbGbtTsK/AJiyJQFrN1RQEGxO+FWEG63G5dXYQLjJ9yHOeNy6y/g8qqYrCBuGf8FAKMds0jCiwgk4WW0Y1ZQe01JWOEgIheIyK8i8puI3B/Lez3zRS6TZq31CwSlFJNmreWZL3Kjep8RI0Zw1113cdZZZ3Hffffxww8/0L9/f3r27En//v359ddfAfj2228ZMmQIAA8//DA33HADZ555Jh06dAgSGhkZGf7jzzzzTK688kpOPPFEhg8f7v8sn332GSeeeCIDBw7ktttu8183kB07dtCqVSt/+YQTTiAlJQWAN998k759+9KjRw9uvPFGvF4v999/P0VFRfTo0YPhw4dHdYw0RzeGYZCR4iTFabB2RwEdxs5m7Y4CUpxWfTS9caKB1+sFYJjxFY8l/QfZ8A2PJf2HYcZXQe3RpMD391wjOGSQXS4gOiTWSPsQEQfwL+BCIBu4VkSyY3EvpRT5xW6mLczzC4hJs9YybWEe+cXuGq8gQsnNzeXLL7/kn//8JyeeeCLfffcdP/30E5MmTWLs2LFhz/nll1+YM2cOP/zwAxMnTgz7NPLTTz/x7LPPsnbtWjZs2MDChQspLi7mxhtvZPbs2SxYsIA9e/aEvf4NN9zA3//+d0477TQefPBB1q1bB0BOTg4zZ85k4cKFrFixAofDwVtvvcUTTzxBWloaK1as4K233ore4GiOekzTZPX2fEo8wSuEEo9Vn2grB5t7HDMBsNX9djkWDP6d9TfUsiAh7TUlUfc59AV+U0ptABCRd4DLgLXRvpGIMH6IJXemLcxj2sI8AEYOaM/4IdlR97656qqrcDgcABw6dIjrr7+edevWISLlLkEvvvhiUlJSSElJ4ZhjjmHXrl20bt066Ji+ffv663r06EFeXh4ZGRl06NDBv+fg2muv5ZVXXilz/R49erBhwwbmzp3Ll19+SZ8+fVi8eDFfffUVy5Yto0+fPgAUFRVxzDHHRG0sNJpwtGmUwq+7j9BLculn5LDE7MJy1Zk2jVLi3bUy2CuZLDkcVG+XY7HS6datKfy2jyK34rY5xfy1dzKntHQEt0eBRBUOrYAtAeWtwKmBB4jIaGA0QNu2bWt0M1tA2IIBiIlgAGjQoIH//UMPPcRZZ53FBx98QF5eHmeeeWbYc2wVD4DD4Qhrrwh3TFVWPRkZGQwdOpShQ4diGAafffYZycnJXH/99Tz++OMRX0ejqSmGw0EvyeWt5MdIwoMbJ8NdYylynBLvrpXB/o25lQOnlP4u3coR1B5NnvhwHwXLPyXjC0uB1CRVgoTDEx/u48N+Nb9PQqqVKLtiAggaZaXUK0qp3kqp3s2bhw1HHjG2KimQQBtErDh06JBf1z99+vSoX//EE09kw4YN5OXlATBzZvil7sKFCzlw4AAALpeLtWvX0q5dOwYPHsz//vc/du/eDcD+/fvZtMmK8JuUlBRTdz3N0YlhGBw64uI0Rw5JeHCKSRIeTnPkcOiIK+FsDrbn4WFJD6q3y5V5JlaVHTt28NHfh7D/ixcB+NPJSTx1XmrQMRX7GUZOYo10KVuBNgHl1sD2WNwo0MYwckB7Nj5+ESMHtA+yQcSKe++9lwceeIABAwbExHCVlpbG1KlTueCCCxg4cCAtWrSgYcOGZY5bv349Z5xxBieddBI9e/akd+/e/P73vyc7O5tHHnmE8847j5NPPplzzz2XHTt2ADB69GhOPvlkbZDWRBWv18uughIWe7vgxolHGbhxstjbhV0FJTH5ndQEW0X8nXkSAPZ0YZft9mhw0003cdxxx/nLebdlMP3ytDLHXRQlza/E+um4OoiIE8gFBgPbgB+B65RSa8Id37t3bxWa7CcnJ4cuXbpEdL9nvsglv9jtVyXZAiMrNYk7z+1co88SbwoLC8nIyEApxZgxY+jUqRN33nlnTO9ZlbHXaALxeDx0enAOCsrYHARY98j5OJ2Jow0vKCjgpEe/4yXn05zvWIqIJSDmeHtzk+cuVo07nczMzBrd46effqJXr17+8pgxY5iVcSE5SX8k1fD671lsOujifoM3hzRk4MCBEV1bRJYppXqHa0ucUQ5AKeURkVuAOYADeK08wRAN7jy3M0opv43BtkHUh1AQr776KjNmzMDlctGzZ09uvPHGeHdJoykXp9NJn/aN+GnzQZabnVnutR7Okgzo2bZRQgkGgPR0S33UQXYE1dtlu706eL1eTjvtNH788UfA0gTs2bOHAwcOMOuFlWXmJ7scrZTJiapWQin1mVKqs1Kqo1Lq0Vjfr7yBruvYm+/Wrl3LW2+9VaMvq0ZTG2w9UARYK4ebHR/5d/7a9YmEbXfboFoG1dvl6trlPvzwQ5xOp18wzJo1iyNHjtCgQQO/ai3PtPRHtvLHLkdL9ZZYYlij0RzVeDweduaX0INc/pv8iN9b6VrXg6zI74zH40mo1YNtIP/W7M75jqX+ifpbs3tQe6QUFBTQsGFDv61z0KBBfPvtt0HXsaMXNDKOAPjVSnY5WrlUEnbloNFojj5EBAGGOuaTjAdDIBkPQx3zERJvRW+vDM40VgKlm+DsclVWDo8//jhZWVl+wbBq1Sq+++67MgImOTkZgCbkB9XbZbu9piSOCNZoNEc9hmHQLCMFwmiQmmWkJJwrq/2U3sf4JajeLkfyFJ+XlxeUHOu2227jueeeK/d4e0+TCvH4t8uBe55qQqUjLSIDROQLEckVkQ0islFENkTl7hqNRhOAiNC2SRrvewfhwolXCS6cvO8dRNsmaQm3crA3pGZyJKjeLlcUYFMpxbBhw4IEw65duyoUDFCqqgoNJGKGtNeUSFYO/wHuBJYBieVknMDs27ePwYMHA7Bz504cDgf2Zr0ffvih0qXft99+S3Jyctiw3Lt27eLPf/4zW7Zswe120759ez77rPzwwAcPHuTtt9/m5ptvrsEn0mhij9vtZummgyg6c63rwWBX1k0HcbvdCZWfvKLJv6L2RYsWMWDAAH/51Vdf5S9/+UtE9ywqspZVHpIInJKtstUeGImhukQiHA4ppWbX+E5HGXbIbrAiq2ZkZHD33XdHfP63335LRkZGWOEwfvx4zj33XG6//XYAfv755wqvdfDgQaZOnaqFgybhcTgc/lAIy1WpKytYIRKiuaksGth5T7wYJAU8y3t9SpnCwkKysrL89S6Xi+zsbNavXw9Aq1atWL9+fZVUQfaxm8wWdHVs8tdvMlsEtdeUSNYf34jIUyJymoj0sl9RuftRxrJlyzjjjDM45ZRTOP/88/27jZ9//nmys7M5+eSTGTZsGHl5ebz00ks888wz9OjRg/nz5wddZ8eOHUGB904++WT/+6eeeoo+ffpw8sknM2HCBADuv/9+1q9fT48ePbjnnntq4ZNqNNVDqVJNeqgrqxCbWEU1wVZz2cFi7e7Z5UA12BtvvEFKSopfMHzzzTds3bq12pP5VppXWK4pkawc7IB3gbvoFHB2VHuSCGz5AfLmQ/tB0KZvVC+tlOLWW2/lo48+onnz5sycOZNx48bx2muv8cQTT7Bx40ZSUlI4ePAgjRo14qabbip3tTFmzBiuueYaXnjhBc455xxGjhzJcccdx9y5c1m3bh0//PADSikuvfRSvvvuO5544glWr17tX8loNImKiGAIdKds4L2VdE44m8Npzy0HIM2w1Ee2W6ldPu255Sy/9zSaNi2NlHrppZfy4YcfVvuzVLZXKVp7mSoVDkqps6Jyp0Rnyw8w41LwusCRDNd/HFUBUVJSwurVqzn33HMBa6NKy5bWRhk7RtHll1/O5ZdfXum1zj//fDZs2MDnn3/O7Nmz6dmzJ6tXr2bu3LnMnTuXnj17AtaSdt26dTWOWqvR1BamaeIwoB+lgfdQHvoZOaymM6ZpJpTH0trxZ5E96Zty23/v/pqmTUsTbK1bt47f/a5mCRdKSkoqbY+GgKhUOIhIQ2ACcLqvah4wSSl1qMZ3TyTy5luCQXmtv3nzoyoclFJ07dqVxYvL5sH99NNP+e677/j444+ZPHkya9ZUHimkSZMmXHfddVx33XUMGTKE7777DqUUDzzwQJkQGXZUVo0m0XE4HDRKS2LJYSvwHspaOSwxu9CoQVLC2RwKCsLnXftlr5dNLw7haV95/PjxTJw4MSr3rGwMojVGkYjg17Ayz13te+UD06Jy90Si/SBrxSAO62/7QVG9fEpKCnv27PELB7fbzZo1azBNky1btnDWWWfx5JNPcvDgQQoLC8nMzCz3i/f1119z5IjlKldQUMD69etp27Yt559/Pq+99prfSLZt2zZ2795d4bU0sSHaOcmPFqzMjB6Wq85MdP+RhWZXJrr/yHLVmfziquUoqQ3sidi2MZim4vw3D9P1xdLkPwcOHIiaYIDKw2PUZviMjkqp3weUJ4rIiqjcPZFo09dSJcXI5mAYBv/73/+47bbbOHToEB6PhzvuuIPOnTvzhz/8gUOHDqGU4s4776RRo0ZccsklXHnllXz00UdMmTKFQYNKhdWyZcu45ZZbcDqdmKbJX/7yF3+2tpycHE477TTASuDz5ptv0rFjRwYMGEC3bt248MILeeqpp6L62TTB1Ocov7FGKUWJR9FLcpmQ9AZJeOhr/Equqw0/eTonnHDo/Y/vAfid521e2nA1F/63dL9Ds0vvpUGX02nUqFFU71mZATta3kqRCIciERmolFoA1qY4wu5frAe06Rt1oQCWK6vNd999V6Z9wYIFZeo6d+5crovqPffcU67X0e233+53cQ3k7bffjrC3mpoQmJMcrIyCgflCAqP/asKjgH5GWZtDoFtropD3xMW0/dv77Jz6By4stqbFbsc46DjiCVZIF/KeuDgm9zWARiFpfRpRGNV4SJEIh78CM3y2BwH2AyOi2AeNpt5Q2znJ6xt2bKUlZhe8GIhSeDFYYnZJyNhKL7zwAluevtVf/v4vDTilVRL/dOfSs9+FMXkY8Hq9mEATCVYVN5ECTF97NOwOkXgrrQC6i0iWr5xf8RkazdFNbeYkr480TndCEQjK//LXJwg7duwIysp22qBBLDwvF/G6UI5kmp5wNvmpSTH5n9uqtf0qOImQXY6W6q3c0RaRPyil3hSRu0Lq7Q48HfbEBEIv4WufRNMJx4PycpJrAREZhghDHfNJwspylqS8DHXM51npGu+uAVa6zpdfftlf3rJlC61atUK2/gh585H2g7ihdZ+Y/69/oxWn8ktQOZpUJIrt4Bzhctwl/AyQmprKvn37aNq0qf5B1hJKKfbt20dqamrlB9dTQnOSB9ocQK8gImHvYTfNnMGe8s04xN7D1UucEy1C03U+/fTTwSl3A2yWsfwP2yqj1WZ7cJTuyl5ttg9qrynlCgellC0av1RKLQxs8xmlE5rWrVuzdetW9uzZE++uHFWkpqYGhfY42hARslKTgmwMtg0iK0ZqhvqEHahuLw2D6u2yx+OJWr6CSPF6vZx66qksW7YMsHYg7969OyrB7aqDYRg4JTiHhFJW+X9qcK1GZZ0ChMZSCleXUCQlJQWFwtVoaov6nJM81thPvfkqDSh9KrbLtb0J7oMPPmDo0KH+8qxZs7j44th4IEWKaZp4FHSRvKD6LpKHRxG1XeQV2RxOA/oDzUPsDllAYm1T1GgSjPqakzzW2ON0mpHjK1sC4jQjB7y1N46RpOuMF/YmNwlRXtllr9cblVSqFX3SZCADS4BkBrzygStrfGeNRqMJwVYrJRGcB8EuV5Y/IRpEmq4zXtj9WqvaBdXb5Zh7Kyml5gHzRGS6UmpTecdpNBpNtLAn4EwOB9Xb5VhO0KHpOm+//XaeffbZmN2vutg7oL81u3O+Y6lf9fat2T2ovaZEsvY4IiJPAV0BvxuKUqr+hezWaDQJQW2q5ZRSXHPNNbz33nv+ut27d/szNyYa9uqpiRRiAg4Bj7LKdns0suVFIobfAn4BjgcmAnnAjzW+s0aj0ZSD23RUWI4WixYtwjAMv2B49dVXUUpVWTDUZqBF0xflL4PDGFg2GYevHNheUyJZOTRVSv1HRG4PUDXNi8rdNRqNJgA7P3Ka4Qqqt8tFRUVkZobbelU1XC4XXbp0YcOGDYDl+v7bb79VSyVT24EW7cm/PKN9tIRDJCsHe+fJDhG5WER6AkevI7tGo4kZ9uS82OwClLqy2uVo6NNff/11UlJS/ILhm2++YcuWLdW6dmCgxUmz1gZtgswvdsdkBWF7Iu1SjYPq7XI0PJUgspXDI76ge3/D2t+QBdxZ8SkajUZTdezJtMhn3rSfiu1yTSbbffv20axZM3/5sssu44MPPqiRPSMegRbta77iHcJgx3IcysSLwSveIUHtNaXSlYNSapZS6pBSarVS6iyl1ClKqY+jcneNRqMJw+/YVmG5qtx3331BgmHdunU1yuMcSKCAsInlpsdAAal8extUwJ6HaK1WKhUOIjJDRBoFlBuLyGtRubtGo9EEUFxcDEAXCfaet8t2e6T88ssviAhPPvkkABMmTEApVeM8zoGUF2gxVkZp+7p2cEJDIAkrOGFge02JRK10slLqYEDHDvjsDhqNRhNV7H0M6RIsBOxypPsclFKce+65fPXVV/66AwcORD0rWzwCLdpj0IyywQkD22t8n0iOERG/5UNEmhCZUNFoNJoqYU+kXhU8NdnlSCbauXPnYhiGXzDMnDkTpVTUBYPdn3CBFkcOaB+zQIv2NcsLThite0Yyyf8TWCQi//OVrwIejcrdNRqNJoDCQmsjlxEyv9nlwsJCMjIywp5bVFREq1atOHDgAAAnnXQSy5cvj5r3TnnUdqBF+7rvewdxlWMeTuXFg4P3vYOC2mtKJJngXheRpYC9I3qoUmptRedoNJrYEJrAqr4ltCr9LKF6cxXSHsy//vUvbrnlFn/5hx9+oE+fPjHoYXhqc0e3223tLliuOnOt60H6GTksMbuwXHX2t0fD5TdSkZqElb9C+d5rNJpaprY3W8UDO1eDQfBGLrscmsshNF3niBEjmDZtWox7WZbaFNqBNoVzjKVc4PiRDA6z3Nu5THuN7lPZASJyO1YIjWbAMcCbInJrxWdpNJpoEo/NVvHg0CHLqBo6rUpIO8Do0aODBMOWLVviIhie+SI3yDvJ/t8880VuTO97r+Nt/uqcxfGyi786Z3Gv4+2oXj+SlcOfgVOVUocBROTvwGKsDXEajaYWiMdmq3iwfft2oHzhsH37dg4cOFBxus5aJFBoA0HeSiMHtI/JCsJeGVzusBJ02hsFL3cs5EnvdbWaCU4Ab0DZS2xTpGo0dZ5YqBlsAWFPRFD/clLbmd4UwZOMApTp5c9//jO//PILAA0aNGDXrl1xS9cJ8RHaduyk/SqLlnLAX79fZfnbo5ExLxIRMw34XkQeFpGHgSXAf2pyUxG5SkTWiIgpIr1D2h4Qkd9E5FcROb8m99EkBrUZsTIReOaLXCZ+siZIzTDxkzVh1QxVGZva3mwVD8pTK32Q42bzU5f5BcOnn35KYWFhXAWDTW3vkLYzwa1QwRv57LLdXlMiCZ/xNDAS2A8c8L1/t4b3XQ0MBb4LrBSRbGAYVu6IC4CpIqJTktZh4qWPjRdKKebl7mb6ok1+ATHxkzVMX7SJebm7gybyqoxN6GarjY9fxMgB7YNsEIlATR8EQlUi+SUKY1I+V/3PitZ60kkn4fV6ueiii2rW0ShS20LbXhW87x2ECydeJbhw+l1Zo5VnOyJvJaXUcmC5XRaRzUDb6t5UKZXju05o02XAO0qpEmCjiPwG9MWycWjqGPHQxyYCPds0ZsWWQ0xftInpizYF1dtUdWzK22wFxGyzVVWJhjdVoHB4bH4JD35T4i+3vOEFnh3eJWHSdUJ8dkjb1yvPlbU2N8GF7V9U7l6WVlhqK5utvrqyHRAZDYwGaNu22nJKE0PqghE12rYBEeG9ZVtwCHgDHhwdAu8t28KES7v6j6vq2NT2ZquqEK0HgYYNG+I+uBHHy/n+uttPTebp81Lp4GpPw4YNKzi79rGF9oj+7YKEtlIqZkLb4XCQnmxwxGWyXHX2u7ACpCcbtbtyCEOl6yUR+RI4NkzTOKXUR+WdFum9lFKvAK8A9O7dOzHW1JoyJLIRNRb7BjweD26vChIMYAkKt1fh8Xj8O3arMza1udmqKviFnQoRdv0jfxBQSnH77bezfeFCf93uuzNolm748zrs3LkzFt2vMRIydYWWo4nX66XEHT6hT4nbxOv1RmVXeLnrMxGZIiLPh3lNARpVdmGl1DlKqW5hXuUJBrBWCm0Cyq2B7ZF+GE3ikahG1FjtG3A4HHRqnh62rVPz9KCnukQdm+ry7JfrUCHPcgrFs1+uq/TchQsXYhgGC32C4ZUhqZjjs2jeIHiKOvbYcM+b8cP/PVoU8j1aFLv9JyJS5uHDxqtqR620tJptNeFj4G0ReRo4DugE/BCje2liTDz0sZESK5WXUoqdBa6wbTsLXH71SiKPTXVQSnGoyBVkYwGYvmgTI/q3K1et5HK5OPHEE9m4cSNgTf7Jw1/khvTrg6/v+5uVlRWT/leXeCX7cRqCxywrIZyGxF44KKVmROUOYRCRK7A20TUHPhWRFUqp85VSa0TkXWAt4AHGKKWi45elqXUS3YgaC5WXYRi0bZTC3sKyAqJtoxS/MTXRx6Y6lKdKKa/+9ddf5/rrS4XAt99+S/fu3enxxEL2qUyaS4G/bZ+y8kY3b948ij2ODrWtOjUMg15tG7FudyEHjrj99Y3Tk+h0TEatboKLOkqpD4APyml7FB31td6Q6EbUcGqdmvTP7XazcltB2LaV2wpwu90kJVnhyRJ5bKqKiJCVZhlmA1cPI/q3IystWNiFpuu8/PLLef/99xERDh48CEAxKUDpOFrlxCQW36PK7te1VUN+yDtQZtXZtVXDqHkBJo5PmKbekohG1FjtG0hKSiIt2RE2/ENassMvGPz1CTg21Ufx/YZ9QTVWuXQsQ9N1/vbbb0F5nG1DajKeoOvY5ViH364q8dh/ErTq7F6ILHia8d0Lo55DotKRFpEBSqmFldVpNHWJWKl1lFJc1btNkIoBrOnxqt5tavxUl6ghu03T5Iu1u8jZWUh2y0xm3TqQIVMWsHZHASK7OO84D926dfMf//DDDzNhwoQy17F39x6iAS046K8/RIOg9kQhXurBO8/tjNr8PfL6ZeB1IY5kxv/pI6TtqVG7RyRieArQK4I6jaZOEQu1joiwZtshmqQnsT9AH9wkPYk12w7V6NqJHLLbMAwOHrHsLGt3FNBh7GzAEl6LptxBtzt+8h9bUbrO1NRUAKZ5L+Bx4z9+F9Zp3guC2hOJeKkHZdMC8LpAeS0BsWkBRFE4VOTKepqI/A1oLiJ3BbweBnRIC029INpqHdM0KSzxsP+Im+yWmWx47EKyW2ay/4ibwhKPP2haVUn0kN2mabL9kLWbuZfkcrPjI1pt/JjNT17CofWWYHj33XcrTddpq41yVRvcOFCAGwe5qk1Qe6IRF/Vg+0HgSAZxWH/bD4rq5Ssa6WQgw3dMZkB9PnBlVHuh0dQTDMPg3OwWQPATdHbLTM7NblFtT5LAFcnKRXN5aslzrDS70CS9a41XJNFgz549gCUYpssE2v6zAJ+sIKl5ezb+tIBWrcIGOwjC47FsC/2MHAxMDAFDmfQzclju7YzH44naDuA6T5u+cP3HkDffEgxt+kb18hW5ss4D5onIdKXUpvKO02g0wdx57gncPriTXzAAzLp1YI1cDO0VSfuiNfw3+RGS8ODGybVFD1JYcgqmacY15pC9c3nkzqdoOL3Uy+iTkS255ZgX2Lt3b0TCwWa/ysCBQilwoNivwueNPupp0zfqQsEmkm/Tv0WkkV0QkcYiMicmvdFo6gFV2fkcaRRTwzDYX1jMUMd8kvFgiOXBM9Qxn/2FxXEPRpeens7+L1/msuk7ABjRIwlzfBYXtT4MwDHHHBPRdexVQTcjD7AS2QSW9aqh9ojkG9VMKXXQLiilDmClC9VoNCEopbhi6kL/zuCNj1/k9/2/YurCaofs9ng87Cxwl6kH2Fng9qtj4sEnn3zCmWeeScGyWdzcJ4mD92Uy7bK0oGMiTXhfmXos3uqzo4lIhIMpIv6wpyLSjggC72k0Nkdbsh87fqSE/A2MK1ldA3Maxb7zg8vxYMeOHVx11VVceumlpKenc+wfnuL5C9LISik7gW/evDmia9oG+/e9g3DjwFSWQdrOVVBdg76m6kRi+h8HLBCReb7y6fhCZWs0lZHI7pexQET44Ob+TPxkDdMW5TFtUR5g7RSecEnXIHfHrFQnXY7NCIrJ0+XYDLJSneU+Ifc0fvOdbwkIu1ybmKbJq6++yn333UdxcTGPPvooPXv25K/zyp+49+7dG9G1A9VjyidMVYBQjbf67Ggikkxwn2PtaZiJlQHuFKWUtjloKiXR3S9jxbBXloTdKTzsldJUJVbGuD3k7CwMOi5nZyHzcveUGRt7Uvzc28d3PkHl2po0165dy+mnn85NN91Er169WLVqFWPHjuV3v/tdhee1a9cu4nsYWN5KDp+3kgPT572kqU0idRr2AruBVCDb9wT4XSXnaI5y6kKyn2hjmib5Ra6wk36XYzOCvIp6tGnEii2HylyjR5tGZepcLmuD2ZPe6wC4wPEjn3v7+MsulyumG8RKSkp47LHHePzxx8nMzGTatGlcf/31/v9hRkbF3kSVtduICCawxOyCGycoyytridkFE21zqE0iCZ/xF+B2rNwKK4B+WGk7z45pzzT1gkRO9hMLRIQUZ3iPmhSnI+hz/7T5QNjjwtUHrgye9F7nFwrh2qPNd999x+jRo/n1118ZPnw4Tz/9dBnvoyNHjgBgErxD1gxprwzbsL5cdWa4a2yZFJgej4fk5OSafBxNhETyjbod6ANsUkqdBfQE9sS0V5p6Q31LaBMJO/LDG4lD63eG7CjuJblB9YEETv6hx4e2R4sDBw4wevRozjjjDEpKSpg9ezZvvvlmhW6p4QIOVoVAV9XlqjNTvZf5BUNouya2RKJWKlZKFYsIIpKilPpFRE6Iec80dZ66kNAmFoHsLux2bJmkN3Z94H28pkkvyeWt5Mf8m9qGu8ay2exaph9ut+XG2ktyeTd5Ig4UXoSrXRNYrjrjdrujFlpCKcV7773Hbbfdxt69e7n77rt5+OGHadCgQbnnNGnSBAAvBgalhmmv7/nTbq8MK5ENeMLYtp2GVivVJpE8bmz1bYL7EPhCRD5Cp+7UREB5ESujHVq4ulRln0GkKKWY+eMWoOwT/swftwStmA4Weehn5JCEB6eYJGGVDxaVv2dhkmMaDhQi1s7hSY5p1e5rODZv3swll1zCNddcQ+vWrfnxxx956qmnKhQMAEVFRdbfkLwLdtlu19QdKn3UUEpd4Xv7sIh8AzQEPo9przT1hkRNaBPoSQUErWpGDmhf4xVEL8kNDnPhepAcTgw6plFaEkuOlDW8NkpPKnM9uy8djeDnMrtc0/H0er288MILjBs3DqUUTz/9NLfeemvEq5H0dCtvtpPgkNp22W6vDNM0KU/jqBRxDxNyNBHRf15EegEDsTa/LVRKhU+Sq9GEIRET2sTKk0pEuPqU1nRa+irJeBCBZGWFuVh3yjlBQvKik1rw+hJXGcPrn05qUe79XcpJqriDyjVlxYoVjBo1iqVLl3LhhRcydepU2rdvX6VrlO7QDp3ZVUh7xTgcDlKcBkfcZfVKKU5D2xxqkUpFsIiMB2YATYFmwDQReTDWHdNoYk2ggLCJxqrmf8u3RVQ/Z80uoKzh1a4PxN4ZfIhgl1C7XJ2dw0eOHOG+++6jd+/ebN68mf/+9798+umnVRYMUJpnITSdp12O1M1WRBh1egdOaBGsxjqhRQNGnd4hIR4sjhYiWZ9dC/RRSk1QSk3AcmUdHttuaTSxJxaeVEop0pIcYcM/pCU5guwbHjP8fTymKrcP21XTCsuRMnfuXLp168aTTz7JiBEjyMnJYdiwYdWefIuLLU+sQhUsBOyy3V4ZSikOFbn5ddfhIJvNr7sOc6io/m6cTEQiWZPmYW1+s/+7KcD6WHVIo6kNAj2pHupewA2ttvHatlZMjoInlQpQrUg59VC+ei1cvZ0es4MErz7scqTpM/fs2cNdd93Fm2++SefOnfnmm28488wzIzq3Imybwg6a0ZbSUBk7aBbUHgkrthyil+QyM/VxnHjw4OSa4gdYsaVRjfupiZxIhEMJsEZEvsBSIJ6LFWvpeQCl1G0x7J9GExNsT6qHuhdww/o7kFwXNziSofuz5NfAk0pEyC+ybAxJeBGBJOVlqGM+E4tO8F/XMAw6Ns+geYabtF3L/DaHohankJWWVMbompRkGambSPDOa7tst5eHUoo33niDu+66i/z8fB566CHGjh0btV3VtnA6GJJ3wS5HKrxEhDM6N+fkpJ04t7sRTJwoxhy/k5+Pb67VSrVIJMLhA9/L5tvYdEWjqV3uPLcz6rtPkFwrD694XdzQahty+rBqX1MphaMcZa3DCN5H8c7ofrz61tv88UDpPoc3mj3LqOHXlTnXtinkqzQaS+lu43yVFtQejt9++42bbrqJr776iv79+/PKK6/QtWvX6n7EsNjCqZtsCKq3y5UJr0DuPLczaumJyHbrMwkmZ/c6kcG961+gxkQmElfWGbXREU39JRYbzaJG+4FW/l2vy5eHd2CNLmeaJi6P4n0ZxFWOeTiVF4/P5uBSyu+KqZSi3+Nf8fvD80hOsrIlC24O5syj3+MtWPLA4KAxsp+899KYxpQKh700DmoPxO12889//pOJEyeSnJzM1KlTufHGG2PiCmpfs7kEx4qyy1W9pxTtBzFAmSCGVdbUKpHEVhoCTAba+Y4XQCmlsmLcN009IJFDdlt9y2D8nz5CNi1AtRvIpJUZZP2aW+2+ORwOmmcms7ygMxPc13Oh4wdme/uyXHWmRWay3xVTKcW+Qhf7ycDwWSMMXzrMfYWucgXofpVZYdnm+++/Z9SoUaxatYqhQ4cyZcoUjjvuuGp9pkiw+3pEpZAcsLI5olKC2iOm/SBwpAQI7UFR66smMiIR588C1wNNlVJZSqlMLRg0kZDIIbuD+rYyAzXwLiatzIhK347NSqWX5DIxaQYDjdVMTJpBL8nl2Kxg/X7DNCdNpBBTWemAvEpoIoU0TCv7zGYHmzsY4spql+32goICbrvtNk477TT279/PBx98wP/93//FVDBAaXiPIxLsymqX7faIadMXrv8Yzh5n/Y1RnmRN+URic9gCrFbah0xTReyENtktM4M2mmW3zKwwoU1t9W38kGyUqYL6NuK0djXyVDJNk5Xb8nnEOd+/Cc7wbYJ7cFtnTNPE4bCis7ZunMaSbV3w4ADlxYuDJWYXWjdOK3N/e3JtHRLz0i673W7mzJnDmDFj2LZtG2PGjOHRRx8lK6t2nuPsFdESM5uhjoX+Xc5LzOyg9irRpq8WCnEkEuFwL/CZLxOcP1ykUurpmPVKUy+wns49rN1REFS/dkcBp3ZoGnfbw7BXlpBfHPxE+33efoa9soSZN55WrWtW5pXj9Xr9E+XOQyW0AgxMBOUPWBcuKqs9Tk2M4LFsYhTgKdzPH/7wBz766CO6devGe++9R79+/arV/+pi2xSKSPX11wp3YZeranNIaDvVUUIk/7FHgSNYex0yA14aTYWICA9d3IXslsFfl+yWmTx0cZe4/thN0yS/2E1OiODK2VFAfrG72rmKk5KSSHVYOZBdOPEqwYWT972DSHWUeu14vV52FZQw2jELpy/jmROT0Y5Z7CooKSNk7PIOsxHgizOkFC8s9bD933/l888/57HHHmP58uW1Lhig8sm/KsIhFgERNVUnkpVDE6XUeTHviabeoZRi8qc5YVcOkz/NiWsAPhEhuRyf02SHUaN+lXhhOZ251vVgUMwkCbOoaCEHKizb2MKhsXEYgJy9XkZ/UszCLQWktjuZ+f/7D7179652n2uKLUzf95b10rLbI1EtxTogoiZyIhEOX4rIeUqpuTHvjaZeISJkplg2h0ABkd0yk8yU+NocoFRV00tygyfxGvTL4/H490EvV51Z7i31elK+9qSkJBwOB+lJBjO9Z9LDWO/X0c/0nkl6UtkAcykplmF3nasZ0xdv4okFLjJThHuH/I53sh8lOzs4RlRtYwuv5aoz//Fc4E9jaseL8nq9Ee11SGQ71dFGJMJhDHCviJQAbrQrqyZClFIUlFg2h9BkP/G2OYgISQb0NnJ5I6l0E9of3WMRo2/1Q2dU4rcRqCop8Zi8owbT1rPLP5m+Yw7GocwyY2OaJsVbVjNmzo+s3+dh+ElJ/PO8FMy0Ema6pdpqsGhhh/YeZnzFX52zAPircxabVQveMQdHHPo70e1URxOVKgJ9rquGUipNu7JqqkIiJ/vxeDz8vC2fvhKcbKev5PDztvyIQ0yHUpnqxG43TROvslYtI51zaSt7GOmcSy/JxauCdzwfOHCAW265hV1v34/HazL7unTeHJrGMQ0MMsUKeRZvZ0LbpnCh4wfAMkgHliO1OSSynepoo1xx7svhUC5KqeXR746mvpGoyX6cTicuj8kSKZtsx6XMaqfcDJzUQ9VVdrvD4fCrYQIzwaGsTHDLvZ39Xk2B6Tqz+g5l6hnbOT9trV8N9aNpXffw4cNkZsbPT8T+3GvMdpxurPL3b43Zzt8eiYBIZDvV0UZFv4B/VtCmgLOj3BdNPSURk/243W4Ulo48NNmO3V6VeEA29qTfS3KZmTwJJyYeDK5xjWe56uzXvdtP+kvMssIJIC8vj7/97W98+umnnHLKKcycOZPrZxdiOC0PcttVtMSXL8G2ScQLwzBomp5EYUkDTASHKLxKKKQBTdPLBhIsj0S3Ux1NlCsclFJn1WZHNJraxOFw4AC8lDUcO6jmpi1K1Tv3Ot7BiYkIOJXJvY53GOYZX0b9s1x1ZqL7j/4wG8u8HSlY9hGnvPB2ULrOgoICmL2IPvJL0Pmh5XhyqNjDElW6qc/j29R3qDhyFV0i26mONmqeY1CjqYOYpkl529W8VD9XsS1UOhrBeRfsst0e6Ck1IekNkvCQsnsNL3/8Pw7u3ML555/Pyy+/TLt2llrGVnNlSFHQde1yddVg0UJEUKYCwbehTyE+vy1lRj6h23aqEf3bMb57IbLgacZ3H4hS7eJupzraiEumbhF5SkR+EZGfReQDEWkU0PaAiPwmIr+KyPnx6J+m/hNqG7AzjoVrrwoul5VePTS3s1222+3r9zNycLvdjPvyCANfPYgnfy/NLr2X999/3y8YoFRd5Q35ydrlSPMlxAqPx4MJwXkssPJYmESeQ9qmTeFqeP0y+PpReP0yq6ypVeL1uPEF8IBSyiMifwceAO4TkWxgGNAVOA5rj0VnpVR8v/maeoe9KuglubyVXOrKOtw1luWqc7XDWtuTfrIET4Z22W63J8uP1hu8OKeAvIMmN/RMZtfpY1md2rPMZF9SYoXUSCY43IddttvjhdPpJD3ZAWFkanqyo4qurG7ca75GJZUgmJieEvat+ZqifidptVItUu4vQER6VfSqyU2VUnOVUvavZwnQ2vf+MuAdpVSJUmoj8BugI29poo4dyK6fkUMKLpxikoKLfkZOUHtVsSfpfLNBUL1dttt37NjB3ln/ZNHMl0l2wDd/Suffl6TSLX1/2Pv7vb1C7ich7fFCROh6XBY5HA/g91bK4Xi6HpdVJbXS+CHZNO16NiXKiUcZuJSTpl3P1p5KtUwk3kqpQG9gJdZ38WTge6BmWVFKuQGY6XvfCktY2Gz11ZVBREYDowHatm0bpa5oYkEiBlGzJ9+OstW3q9P6cneUrf72tLS0Kl/XNji3lL1B9XbZNE1mzJjBnXfeyeFDBdw8sBn/OL2EtCRBKbjB8TnvmIPLrFzskNxuHKQEWEvcOILa44XX62Xd7kJ6qwK/t5JHCQ1VAet2FwYFHKyMZ79cx8EGXYO8yE5q0JVnv1wX9xwgRxOVeiuJyDvAaKXUKl+5G3B3ZRcWkS+BY8M0jVNKfeQ7ZhzgAd6yTwvXlXL69wrwCkDv3r11OPEEJVGT/djqnbONn4BS11C7XF2bgy0cUkPUSqniwX1gO1deeSULFiygT58+bD1pBI8dN4lUo/Rr3wJr5RCqo7ddVUO/6CqkPV6ICEUuj7VPhCSSAlxzi1yeiB8GlFJ8++suVm7NB0q9yJYv3kz31lnccU6nuD9YHC1Eogg80RYMAEqp1SLSo7KTlFLnVNQuItcDQ4DBAbkitgJtAg5rDWyPoI/1ny0/QN58KyNWHYlxn8hB1Ow9DOkUB9Xb5erscQAoKioqU+f2Kv6xyMWO+bdQkJbCiy++yMUXX8ygf/1MugTbCuxyqA2huNjqV3KIj5VdLi4urtZKJ5oYhsFyT7Br7nLVmfQq2G+UUri94Z/13F6VEKvOo4VIhEOOiPwbeBPrQeUPQE5NbioiFwD3AWcopY4ENH0MvC0iT2MZpDsBP9TkXvWCLT/AjEtLUybWkcxYtv7Y9JpBQdSu79c27vrj/Px8oKzRzQhob9CgAVXlyJEjQeXvt3oY9Ukxq3abpHfuz8cvP8qZZ57Jli1bAJAQC65dDl252MKhPCprjzWGYfDnAe35Yf4cJjpm4MRLPyOHjd529B1wfsQGfsMwODe7BUCZTXDnZreISf5rTXgiGemRwBrgduAOYK2vria8gJUT4gsRWSEiLwEopdYA7/ru8TkwRnsqYa0YvC5QXutv3vx49yhiBv79az5auT3IXfSjldsZ+Pev49ovW20TqjwyQ9qryqFDhwAoKFHc9nkxp/3nCPuLFO9fnUbzK8b61T/25F+ecAoVDunp6UD5Bmm7PV6YpslXv+zmEuaRjAeHKJLxcAnz+OqX3VVS091xTmf6tG8cVNenfWPuOEfbG2qTSlcOSqli3+T9mVLq12jcVCn1uwraHsVKMKSxaT/IWjHUsWTrXq+XfYUusr2/8N/kR/zuotcWP8ha94lVMlJGm4IC66nUESIe7LLdXlVcLhdH1n1Pty8L2ZavGNMniUfOTiUzWbjTVbrPwXZVNYHAEbB7E+rKaj8xu0IM0i7f2fF+ojYMg4wUJylOI8gwkuK06iPtn1KKiZ+sYcbizUH1MxZvRkSYcElXrVaqJSr9j4nIpcAKrCd5RKSHiHwc435pAqmjydZFhOObNWCow8qnbAgkY+VTPr5Zg7j+yG0zV3lP7tWJcrp9+3buvXcse96fTKMUYeEN6Uy5KI2slNLPuXNnIVC6wihvJWC3h1JAeoXleKGUomurhvy3ZABeIwmF4DWS+G/JALq2alil8VyxxfrsI/u3Z+PjFzGyf/ugek3tEInNYQLWXoNvAZRSK0SkfQz7pAlHHUy2bhgGh4pcGCEzoCFY9XF82o0070IkmKbJK6+8wn333UdRURGNTv8TSwd+QLKzrPBr1Mj6yaWm+nIth7TbZbs9tD+rzfac6SiNerrabF/l/sYCO+xF9/7n4ejeH9m0AEe7gXRfmVGlsBciwhmdm9OzTSPGX+IL836JlcgoK02Hz6hNIhEOHqXUIf1PqYQYehMl4j6BSDBNk4ZpSfxc0J5rjdKNUT9729MwLana8YuiQWFhYY3abdauXcvo0aNZuHAhZ599NsOHD2dSbguSHB9WeF5l+xJC2+0dxragtV1v7XK8YytBSHj2tqciwPg2Vf+uhg3zfoneAFfbRPLLXC0i1wEOEekkIlOARTHuV93C9ib6+lHr75boOVjV5WTrSil+3VVIEynEizWheYEmUsivuwrj+rRbWbiJytqLi4sZP348PXr0ICcnh+nTp/Pll1/StGmLSs6zDN2VuZ2W127nRwjNl5AoRCs8eyKGeT/aiEQ43IoV66gEeBs4hOW5pLGJkTdR4D4BW0DY+wTyi91xVyVEgqlgv8rAgTWhObDKZpy7npGRUe32efPm0aNHDyZPnszVV19NTk4O119/PSLC4cMV39duP1zJgaHttuG+kAZBgraQBkHtGk20iGQterFSahwwzq4QkauA92LWq7pGjLyJ7H0CQNA+gcC0m4mOYK0UAkMqNJHCsFvha5OqqnXAStd577338u9//5vjjz+ezz//nPPPDw4c3KZNFvycX+5127SxMuxWVTjZaiMrOVBymeRAiaBW0tQvIlk5PBBh3dFLDL2JAgWETV0RDJbeGH9IBY8ycJPEErOLT2cev+WDvQkuknalFDNnzqRLly5MmzaNe+65h1WrVpURDADNmjUDyt8/YbfbBufyjgs1SAM4pDRz3dOeq/wRZB2J/1WoMqHfjbqwSq5vVJRD+kLgIqCViDwf0JSFFQ9JE0iMvIlsVVIgk2atrRMCwuFwkGzAcm9npnnO4wLHj3zu7cNy1ZlUR3xVIY0bNwYOVtIOmzZt4uabb+azzz7jlFNOYfbs2fTs2bPc8+ywG24cOMIEyLPb7ZVJCcmk4/IfV0JyULuNiCACqLKZ60Tql04+UeNxHW1UtHLYDiwFioFlAa+PAZ2EpxYItDGMHODz+R7QPsgGkciYpolCGGZ8xV+dszhedvFX5yyGGV+hkGoHt4sGlYXGSElJ4ZlnniE7O5t58+bxzDPPsGTJkgoFg30egIT8a+xyaIC8A2ZGhWUbwzA4JjOV5JBlQrJDOCYzNe6b4KJFfbCz1Rcqisq6ElgpIh8Ah+0wFiLiAOIbAvIoIShlou8pavyQbJRSdSJlolKKEq/imqRvgVL3y2sc3/KOe3Bcf+j22IXboezatZ4rrniAVatWcdFFFzF16tSgrGwVYU/SSRIs+Oyy3W4LidA53S6HChHTNGmUnsT2Q8Vkt8xk1q0DGTJlAWt3FNAoPb5uwdGkPtjZ6guRfJvmAoF+dWnAl7HpjiYcEmK+DS0nKiJWT3ep4Dg5u1RjhPiqQmy1jQr4CRxxK+77soQdM+5k586dzJw5k1mzZkUsGACysiyD8yEz2BXVLtvt9mffYTay+uGTk3Y5dGwMw+CcLi3IbpnJ2h0FdBg7m7U7Cshumck5XRInIF00bAV12c5Wn4jkG5WqlPLvCPK9T4w9+/Uc/xJ7UcgSe1HdWGIbhkFassG3ZnegdAL81uxOWrIR1wnN9u7x+ATtnN88nPRiIU8vLiHjpHNYtGgRV199dZUnJNuOsklaBtXbZbvdnxTIOAhYq6rAcrj/7V3nncCsW4NzbM26dSB3nXdClfoYK6K1J6c8O1uif9/rG5H8Og8HpgUVkVOAskHrNVHHfoKy7QzHP/CZ3/5QF56klFIUuUyaSCEKn1oJy7W1yGXG9cdurxx2Hhb++EERF7x1hCSHMPtPDWl64W0cc8wx1bquPflvUFaeK/sj2mW73TAMBNhsBt9ns3kMQvhAekopJn8aHC1/8qc5CTFpRstWUNftbPWJSJyj7wDeExE76U5L4JqY9UgThC0gbN0r1J0ltmWQhgwOB6XizOAwytcer9WDaZoUrvqS3t8coKBE8dDpyTwwMAWnQ8BVfddJe2d1B9kJlNpZ7HJJSQkpKSkopUhJMnjSM4x3HRMxlMJEeNI7jJQkI6x6JnDSDEycBPH/TkTLVmDb2QLPs69bF+xs9YlIQnb/KCInAidg/bZ/UUpVL/u6psrUZVdWe4IbHJKKc7DxE096r4vbU+Bvv/3GqFGj2Pftt/Rv4+CVIal0PcaBUuDxdam6OZnt85qrA0H1dtluNwzDHxfJiwMDD147/LaUXTnUhUkzWg8yYWMr1YHve32jon0OZyulvhaRoSFNnXy+x+/HuG9HPYn+tBgpaQF+/OHKtYXb7eYf//gHkyZNIikpiSbn3cy7fT7mOEfphre9KhOwnvBrknaziVFQYVkpRXqyk37eHByYGAIOZdLPyGFzcrewwRUTfdKM5oOMjq0Ufypa05/h+3tJmNeQGPdLQ/lPiyMHtE+Yp8WKcDqddDu2AWtVsLfPWtWObsc2qNWQD0uWLKFXr16MHTuWiy++mNWrV5PZ8yI+UqcDpbaBD0wr9Ell4S3Kw37iz1fB+yjsst1ueXIpXzgMp2/3uBUOQyg/kmmiTpraVlD/qGifwwTf35qmBNXUgER/WqwIr9fL6p2H+dbozvmOpUHeSqt3Hsbr9cZcQOTn5zNu3Dj+9a9/cdxxx/Hhhx9y2WWX4Xa7gVUU0gCFpcoxVWkgu+pOZvZ575sD+asxy/+Z3zcHBrUrpVCIPxxGPyOHJWYXlqvONEPqTFh2m7qg9qpvxDqUf0Vqpbsq6djTUeuFpkJEJChfhNSRpD92HuZuRh5QanPoZuSBabXHUjh89NFHjBkzhu3bt3PLLbfwyCOP+PcZ2OxXGUHG8v2qeisGG9sbqTyDdKC3UrJD/PGS7HAYDrF2PSfKvoWqUJcfZOoatRFipKJvYKbv1Rv4K9DK97oJyK7gvDpDtIJ7xTxIWAzzRcQSe1I4jTVB9XY5VpPG9u3b+f3vf8/ll19OkyZNWLx4Mc8//3yQYLBzOVsRY4NzTQS2V5cWsr/Cstfrpdht4g35qngVVn1IDum6QqKqveoTtRVipCK10kQAEZkL9FJKFfjKD1MPwnVHS/L6r9O9ENm0ANVuIJN8qRGjFiQsXL6IOrB6sCeGFkaw545djvbEYZomL7/8Mvfffz8lJSU89thj3H333f5gd4HY4SkyOIxBaa6JDA4HtVeXmd6z6GFs8KuVZnrPCmp3OBwMP7Utr8zfSImnNNRGitNg+KltdX4GTbnUVoiRSNb0bSHIvcQFtI/K3eNEoOQFgryARg5oH7Huzr7OykVz8Sx/HCcePDhZWfwA3fufFz0dYIzyRcQa+7MfNBvQwFGaWe2g2SCoPRqsWbOG0aNHs2jRIs4++2xeeuklOnXqVGHfHAZ0NTb5ypaA6GpswqGq3zfDMHAa8I45GNxwoeMHZnv78o45GKdRapBWSlHo8gYJBoASj0mhy1vnbA6a2qU29j9FIhzeAH7wBeBTwBXA61HrQRyI5oad8UOy+XrH68g2NyImohRjjt/J2dH8R9n5ImKUozpW2KoRZ0gQOrvs9XrDPtVXheLiYh577DGeeOIJMjMzmT59On/6058qHXuHw0HP1lns3W6pmuwn/L0qi56ts6r95B5433fMwZaQCNMuIqzeepAUp1Fm5bB660EtGCCmednrOrWx/6lSq5dS6lFgJHAAKwD+SKXUY1G5exyJVnAvEeHsC4YGuSOefcHQ6P+42/SFQX+rkz8Sh3grLFeXefPm0b17dyZPnsw111zDL7/84k/XWRlut5ulm/PpJnlAaWyjbpLH0s35Pm+mqmOaJp5yIpF7TPxhyr1eL+v3HKbEY9IkPYn1j15Ak/QkSjwm6/ccrrM2h6hRR+1stUFtuQ1H6hKRDuQrpZ4DtorI8VG5exyJVnAvpRSTVmYw0f1HFppdmej+I5NWZmi/bkpXDqvN9kDp07ldru4EuH//fv7yl79w5pln4na7mTNnDm+88QbNmzeP+Br2/6e8DXrRck4or93hcNCpRQaN05PYf8RNx3Gfs/+Im8bpSXRqkaFtDjHKy14fqK39T5WqlURkApbH0gnANCAJeBMYEJUexIFo7Ty2r7Ny0Vxmpr6JEw/9nblcs6gNk6pwnfqKrTJqJEeAUr2+Xa6qSkkpxbvvvsttt93Gvn37uOeee5gwYUKliXvCYU++21VT2rDHX79dNQ1qr85105MMjrjLLh/Sk4yg6/br0IwTW2QyY8lmf92lJ7ekUQOdLqWu2tlqi9pwG47E5nAF0BNYDqCU2i4imVHrQRwIkrzdC5EFTzO++0CgapLXvs6Y43fi3O5GMHFi2Rx+jqIEDw1QV1cSu9gqlCQVnFXWLlclE1xous7PP/+80qxsFWHnmjhI8L6Gg2TUKNeEUoq0coRDmi+gnu0dl1/kDhIMADOWbGZk/8idIuotddTOVpvE2m04EuHgUkopESvRoYhU/TEtAbnz3M6ozd8jr18GXhfiSGb8nz5C2p5a9essPRHZbk0GgsnZvU5kcO/ouLFe8/Ji8otcfHrbIAzDwDRNLn5+Pllpycy88bSo3CNW2JN/sgQLB7sciXDweDxMmTKFBx98EBHhmWee4dZbb62x2sVW7+ylYVC9Xa6JWumwK7y67LDPC8nm3aVbAOgluUE7pN9duoXxl9SLrUQ1I0Z52TWREcnj57si8jLQSERGYWWBezW23aodZNOCIL2mbFpQvesU7QfxDaUYVjkKmKbJ+j2F5Ows5OLn5/sFQ87OQtbvKYxrDuZIsJ9k9qvghaZdruxJ56effqJfv37cddddnHnmmaxZs4Y77rgjKvp4r9eLAt73DsKFE68SXDh53zsIRfXtIUopij3hBUuxR/mFg2maeJWil+TyVvJj3OV8j7eSH6OX5OJVKuH/t5r6T4XCQaxf70zgf8D/YdkdxiulptRC32KPrdcUR830mu0HgSPFd52UqOlHRYQhJ1sZxHJ2FtJh7Gxydlo7eIec3DLh1Q52aIxwqpvA9lAOHz7MPffcQ58+fdi6dWu10nVG2rflqjMT3NezwOzGBPf1LFedK+xbZRiGgaOcf4sjIBS3w+Fg9YRzGeD8hSQ8OMUkCQ8DnL+wesK52iCtiTsV/gJ86qQPlVKnAF/UUp9qj2jpNWOkHxURJlzSFYDpizb560f0b8eES7omvHCo7Ok73D6HOXPm8Ne//pWNGzcyatQo/v73v9O4ceNyrlB9TNNEgJ6Sy4SkN0jCQ1/jV3JdbfhJda6RXadxejJ7D7vKqIsap5fmiFBK8ejsX1npOZGbk52gPLhxstBzIoWzfz3qnRk08SeSx6MlItJHKfVjzHsTD6Kl14yRflQpxfcb9gXVfb9hX50wWNqqkdYB3kCB5UDVye7du7nrrrt46623OOGEE5g3bx6nn356zPomIiign5FDEi6cAigX/Ywclns7V3tsRYQit8evLkrCmvSHu8byq/vEIO+SjGSDNY4Tg6KyrnGcyIBkI+H/t5r6TySPRmdhCYj1IvKziKwSkZ9j3TENQTaGQAJtEImMvXL4nbEtqN4ue72WgXb69Ol06dKFd999l/Hjx7NixYqYCgYoNThncBgHZWMr1cQgbSpb6JSqi/oZOZiq9Lper5e3vt9CicckL60rf5v0InlpXSnxmLz1/Ra9CU4TdyIRDhcCHYCzKU30c0ksO6WxEBEOHLF26o7o346Nj1/EiP6W3v3AEXfCP13aE6GTYG8lu/zbb78xePBgRo4cSZcuXVixYgUTJ04kNTU15n2zxy4wtlJwufpj2yDFGTaJT4OU0oW6w+EgLdlBk/Qkfhw3GIfDwY/jBtMkPYm0ZIe2OWjiTkX5HI4BxgK/A1YBjyul8ss7/mgmVkk3RIRr+rTlUJHLb2OwbRAN05ITXjjY+RxCcXsVhxa/y2nPziQlJYWXXnqJUaNG1ereDVtwrTHbcbqxyr97e43ZLqi9qhiGQcfmDcg1uzKx+I/+wHsb07rSuXmDoM+48P7BeL1evyCwBYQWDJpEoCKbw+vAMmAK1mrheWBELfSpThHrpBvhdkLWBWM0WHmYQ1my1cOoT4o5uOd1hgwZwssvv8xxxx1X632zx6+QBngBp4AnIBNcTTbBlXhMji9aw4TkAEN3URtKPH3KPDiECgItGDSJQkWPascqpcYppeYopW4FTq6tTtUVaivpRl1NoFJUVARYGdbySxS3fFbEgNeOcLBY0XzoQzz33HNxEQxg/e8MwtscDGpmc9iy/0hYm8OW/UfKXDfUbpTodiTN0UNFKwcRkcZYv20AR2BZKVXtnV4iMhm4DDCB3cAIpdR2X9sDwJ+xEnPdppSaU937xJraSrpRV1m/fj0AH//q5tbZxWwvUIzpm8zkM1PoIaeyfv16OnToEJe+GYaBUU4+B0NRYxXXErMLJoKpwERYYnYpc8w1Ly+moNjNrFsH+ne/D5mygMzUpITf/a6p/1T0C2iIpVayX1lY8ZWWAUtreN+nlFInK6V6ALOA8QAikg0MA7oCFwBTRSSh19kiwkMXB//wH7q4y1EvGACe+nAHez54jKHvFtEkTVj853SevyCVrBRrbKZ8VbNUnDVBKYXHhNnevr6yVT/b2xePWTObQ4dmDeji2EISXgRIwksXxxY6NCu1OZimSUGxm7U7ChgyZYFfMKzdUUBBsVuvIDRxp6I0oe1jddMQw3YDrCRCYK0m3lFKlQAbReQ3oC+wOFZ9qSnPfPErX6zdFVQ3ZMoCzs1uwZ3nnhCnXsUXO13nohn3UFTi5rGzU/jbackkO4XAOXfs5U3j10nAYUCuaoMbB0l4ceMgV7XBUYNFg2mabNx3hFvEyj9gr0jOkx+Ys+9C/+Y6wzCYdetAv0DoMHY2ANktM/0rCY0mnsTtGygij4rIFmA4vpUD0ArYEnDYVl9duPNHi8hSEVm6Z8+ecIfEHNM0+WLtLtbuKCC7ZSYbHruQ7JaZrN1RwBdrdx2VT39r1qxh0KBB3HzzzXTv3p2WN7zAfQNSSAoTU6Iq+ReijYiQ7DDoZ+RgoBABA0U/I4dkR803oYVbkYRiC4heksvNjo/oJblaMGgShuoFkIkAEfkSODZM0zil1EdKqXHAOJ+N4RZgAqX2jUDCru+VUq8ArwD07t07Lpl1DMPgnC4tAMo8/Z3TpcVR9SMvLi7m0Ucf5e9//7s/XWfHjh35w6xD5Z6zc+dOOnbsWIu9LMUwDE5q1ZBVO7vhVh/4w1esSurGScc2rPb/TkRwSPgc0i0k2JnANE3ufubfQTup736mAf+48y9H1XdHk5jETDgopc6J8NC3gU+xhMNWoE1AW2tge5S7FlXuOu8E7jink18wAEfd09+8efMYPXo0ubm5/OEPf+Dpp5+mefPmDLv/U8DyOgg0HNnrqSmfHGRAnFJGKaXIPi6L6XkdGS5jA2IgdWTEcVk12qvSsmEauwpcZXJIt2yY5n9v2xjO3L+U5CQPDkzAw7H7lzJkyolH3XdIk3hUFpXVEJHV0b6piHQKKF4K/OJ7/zEwTERSfKlIOwExSx4banSsjhFSKcXkT3OC6iZ/mnNUpAmtLF3npBusf/NSdWLQeXbZbo8XK7ZYq5rup53HPY+8TPfTzguqry7d2zSstN4wDDJTk9jZpDeG04ro63CmsLNJbzJTk7Rg0MSdyqKymiKyUkTaKqU2V3RsFXlCRE7AeojcBNzku98aEXkXWAt4gDFKqZgEmYnG5rVopRuN5D6x2IFdk/7MnDmT22+/nX379nHvvfcyYcIE0tPTg46z9jCs40nPMGYmT8KpTDwYPOkZFtAeH0SEMzo3p2ebRoy/xJeH15dgJyutZln8Vm4NH0ggtH7mjadhmqci23pA3nyk/SD+0aq3FgyahCAStVJLYI2I/AC+qGSAUurS6t5UKfX7CtoeBR6t7rUjvL9/8xoQNKmPHBB5ikY7TeiI/u2CEn0rpaKW6DvWO7CrSl5eHjfffDOzZ8+md+/ezJkzhx49eoQ91o6RtFx15hrX+KDw1YHt8SJsHt5LaibQRYTTOzXD5fGydkeBvz67ZSand2pW5tqGYQRF9NViQZMoRCIcJsa8F7VMtDevCQJbfoBNC6DdQCQkuU11iZYQiwah6TqfffZZbrnllgrDPQSq1parziz3di63PV5Ee/e5UoqCEg9rdxSUWU2e2qFp3Fd9Gk2kVCoclFLzRKQd0Ekp9aWIpBNsX6yT2ALCnnih6moge/JeuXgunp8ex4kHD05WFj9A9/7n1XgisFYmTrJbZgYJseyWmWSlOmttkvnpp58YNWoUy5Yt4+KLL2bq1Km0bdu20vMqixNUH+MI2avJkf3bB60mUURtNanR1AaVrmJ9eaP/B7zsq2oFfBjDPtUKtnomEDs+UqTYP/wxx+9ETDeivIjpZszxO6Nib7CEjydIPQGW22x+sSfmT97h0nV+8sknEQkGiCwTXH1FhXhgh5Y1mkQnEhXnGGAAkA+glFoHHBPLTsWaQEPyQ90L2HjRLzzUvcAfQC9081pFk/CzX65jY0bPoNj9GzN68uyX62rcTzs0R3bLzKD67JaZMQ/RMWfOHLp168Y//vEPbrjhBnJycrj66qurdM/Kjq2PT9H2anL6ok1BwRinL9oU1WCMGk2sicTmUKKUcgUY7ZyUszGtrmAv/R/qXsAN6+9Acl3c4EiG7s8yd3s+kz/NicgArJTiUJGL537O5LNAX/mfMxnR31VjtZJSiqEvLg67chj64mI+uLl/1CfY3bt3c+edd/L222/XOF2nYRg4BLxhvi0OqXlwu0REB2PU1Bci+XXOE5GxQJqInAu8B3wS227FnjvP7czIVlsRrwuUF/G6GHncVrKPy6pSCO7Zq3YClsF1qvcyvyeOXV9TlAofgqO8+urfRzFt2jS6dOnCe++9F5V0nSJCt1ZZAEEhIgC6tcqqtxNloICw0YJBU9eIZOVwP1YI7VXAjcBnwL9j2ana4Jkvcsna25obHMmI14VyJDNte2uymjkZ0b9d0FNfoKtqIEopvOWoCbxKRcUzpVfbJmH95nu1bVKj6waybt06brzxRr755hsGDBjAK6+8QnZ2duUnVoJpmvyys4BekhsUImK4ayxrdp6AaZr10ihdnj1LCwhNXaLSlYNSylRKvaqUukopdaXvfZ1WK9l64ckrM3mt47Oos8bxWsdnmbwyk3m5e8sozQQJa0MwDIPvHzibRqnBE1yjVAffP3B2jdUmIkJWmuX5EsjI/u1rvFELwOVy8dhjj3HSSSexbNkyXnrpJb777ruoCAawxsdUhE18Y0YhZ0IiEroxcuPjFzFyQPug1ahGUxeoKIf0KiqwLSil6mxmuMBl/+SFeUxeaYVzGNnf2jswffGmoOOnLcrzt4XuVP79S4s5WBzsdXOw2MvvX1rMBzcPqPEEfsc5nXj44+AIJqYyueOcmoWeWLJkCaNGjWL16tVceeWVMcnKppTCaQhLPF0wxYFCYYqDJWYXnE6plz7/IsLa7flBTgMPXdyF7zfsY+32/Hr3eTX1l4oe3YYAlwCf+17Dfa/PsFxb6zTh9MIPDekSPi4sFbgilic+o/CAqJTiiqkLmbF4MyP6t2Pj4xcxon87ZizezBVTF1brKTQ/P59bbrmF/v37c/DgQT766CPee++9mISysA3/ACjT6q/PVlJfff7tgH5rdxT4Y2xN/jTHCuvuC+in0dQFKkr2swlARAYopQJjZ94vIguBSbHuXCwJpxee/GmOPxzG9EWlq4cR/dvRMC05rM1hy4EjgGVwDQwPseXAkSg9Gfu8xEL+livFKuDDDz/klltuYfv27dx666088sgjZGZmVn5iDWjZMJWhRfOtrGgCScrLUMd8/tewT0zvGy+0t5KmvhCJ0reBiAy0CyLSHyt7W52lon0OBeVsLgunxhERnA7Db3C9y/kebyU/Ri/JxRmFhDEiwgc397cM5IvyOP6Bz5i2KI8R/dtVyY1127ZtDB06lCuuuIImTZqwePFinnvuuZgLhqMV7a2kqQ9EIhz+DPxLRPJEJA+YCtwQ017FGFsvfEXzbdY+h28e5Yb1d3BF8228u3QLMxZvZmR/nzGxf3umL9rEFVMXlREaIsLi+8/mkoYbggyulzTcwOL7z47KZFDeZrpINtmZpsmLL75IdnY2s2fP5vHHH2fZsmWceuqpNe5XJCilKPZ4ed87CBdOvEpw4eR97yCKPd56q2KJxu57jSbeROKttEwp1R04GeiulOqhlFoe+67FDqUUm/YdpuX+pShPCSgvpqeElvuXUuKxdOK2jaHU1lD2h23rkz851CFoh/QnhzpEJaeDvckuUMUFMH3RJg4VuSq8fmC6zj59+rBq1Sruv/9+kpKSatSnqiAi7D/sYrnqzLWuB/mn52qudT3IctWZ/YddifMkveUHmP9P628N0d5KmvpCpfscRKQhVpa2033lecAkpVTNMqLEEaUUjdOTWFLQhRLl9PvfLzG70KlFA/q0b8z0RZv8k/KI/u2YcEnXsBE8124/xMa0rgwvKt0hvTGtK0nbD0Vl8pNybAvl1Qem68zKymLGjBn88Y9/jMtErJTCa1qTYWhUVq8ZnX0gNWbLDzDjUvC6wJEM13/sD59dHWojjLtGUxtEolZ6DSgArva98oFpsexUrDEMg3Ozj7UmdddYnvZcxXDXWDamdeXc7GOZeNlJQceHEwxgqW3W7znMgSNuilr04u7JL1HUohcHjrhZv+dwmRhNVaWq+xy+/fZbTj75ZB555BGGDRtGTk4Of/rTn+I2IYkIrRunh21r3Tg9MSbKvPmWYFBe62/e/Hj3SKNJCCIRDh2VUhOUUht8r4lAh1h3LJYopSgo9nDgiDso7MWBI27yi91l9hVM/GRNWHWAlUzemuBydhbSYexscnYWAuAQicrkd8c5ncJG+Aw0kNvpOs866yw8Hg9z5szh9ddf96frjCc784urVF/rtB9krRjEYf1tP6hGl1NKMS93jxV47xNfCJZPrMB783L3aLWSps4QSfiMIhEZqJRaAJZrK1AU227FHrOc2ESfrNjOviNuvypp4idrmL5oEyu2HAy7qc2euENdWaMRojkwomdo4hh7c9W7775babrOeHJsViq78kvC1icEbfpaqqS8+ZZgqIFKycaOezVtUR7TFuWVqddo6gKRCIe/AjN8tgeAA8CImPWoFlBK8e7SrWHb9h9xA5HtK1BK0bRBCq0LV/Pf5Ef8totrXQ9S1OCUKCX7SQrykbddJD0HdzFkyD0RpeuMJ73aNmbl1rLmqV5tG8ehN2VRSiEBaTqjGQ8r9IEhmvGwNJpYE0kmuBVAdxHJ8pXDZ0+vQ4iUZ86FVKdwTd+2QU99FRmk9x0u4TbHfJLxIALJysNQx3yeP9wtKmolf57jrT9C3ny8rU8jI3cx48ePjzhdZ7wItJkEPkFHKzZUTYlFfm4R4aEhXTi0biGP5z+CEy8eHDyQ9TgPDbkw7p9Zo4mUSDLBPSYijZRS+UqpfBFpLCKP1EbnYoVSivQUSy6GhpIWMVBmWZVQuH0Fpmmyp8AV9h57Clw1Nkj78XnULH9zIqcOPIO7776bs846i7Vr13L77bcnpGCwicRmEg8C83NHGp490utOnpVDrwOfk4wHhyiS8dDrwOdMnlXWvbmyskYTLyJRK12olBprF5RSB0TkIuDB2HUrtogIrRul0fbw6jKhpDcnd2XGks1Bx09ftIkR/duFVTkI8L53EFc7vsWprKfE972DqhHcIjzPfJFLp1/e4ZvPDvLckhKapQsPjrqUxlc+GXG6znhRmc0knruGA1V0KxfN5aklz7HS7MLIAefVuF/LNx+gYzn1gcRi5aLRRItIvJUcIpJiF0QkDUip4PiER0R4/+b+XNJwfdDO5iFZ6xnSPXwAuvIUUU7fQ7uBGfTXGYWHeaUU7374Cdc8OIOnF5cwsmcyP9/ShDVNz+S7dXsT/imzPJvJyAHtE8LnX0QY370wKPTJ+O6FUQl78r53EF4jCYXgNZKsB4aQiL6xWLloNNEikpXDm8BXIjINa5vwDcCMmPYqxtg/xJ8PdWRYshOUtXKYld8R7+aDZQLvlacjdzgcJDkM7pV3cKAQAVGKex3v8Gfj4Rqpe+x0nUvefhtnk9YMuu46OrV3cZPPuJm+q26Yfvw2k9I0swkTZ0gpxdefv88ZvgcElIevP3+fs0f1rXb/RIQzOjcnv815OHr0RzYtwNFuIN1XZAR9h3SAPk2iE4lB+kkR+Rk4B0uLMlkpNSfmPYsxs1ftZJfqzHDX2CCPkmMOFdFTNQo6tjwducfj4YjLpGPy9qD6jsZ2jrhMPB4PTmck8jfgXkoxffp0/va3v1FYWEjTQdfRoO9VbHYmMTUgbYSpqDOZ1MIZ8uON/YCwcuOxnJ6ahMKDEif/2ngsC2qYtS1IILY9FQHGtwmjkvQJCFswgA7Qp0kcIk3FlQN8rpT6GzBfROp0OM/A9J7Fx57C3ZNfovjYUwDIL/YwffGmoLg40xdtqjAuznrzuArLkZKbm8vgwYO54YYb6Nq1Kz/99BMdzrsecZaNh5SV4qiXmdRqC1vl1b3/eThHfoKcPQ7nyE/o3v+8qKi8IhGIOkCfJpGJJLbSKGA00AToCLQCXgIGx7ZrscMwDK7r25Yvc3axdkcBHcbOBiC7ZSYZKU66tmpYZl9BuAnDfmr/0BzAqY5fsH/TH5oDgtorw+Vy8dRTTzF58mRSU1N5+eWX+ctf/oJSit2FeWHP2V3orjMrh0Ql0if8WBAaoC/QWA96BaGJP5HoPMYAfYHvAZRS60TkmJj2qha467wTuO3sjvzuwVIN2cdj+uN0OjFNM0g3/NDFXcI+pZumiQK6GXm+Y0Epq/yOGZnaZ/HixYwePdqfrvP555+nZcuWALjdbv9xoRuq7Ptr4VAzYqXyCtybQvtBqNZ9gq5d0QbHRDDWazSRCIcSpZQrYLJ0EpUkmPHl6pcWsW53QVBdn8e+Ii3JwfldWzL+klL3wsmzcshKK+teaE/ezQjeAWyX3W53uZN3fn4+Y8eOZerUqbRq1YqPPvqISy+9NOgYe8ztZEKBLrfLVWc9gSQoz3yRS9be5VauEK8L5UjmtY7Pkt+sV9B3KHSDo7QfxPghffT/VZMQRKK0niciY4E0ETkXeA/4JLbdii1er5cf8g5w4IiHxulO1j96AY3TnRw44mH7oRKmLcoLCpo2bVFe2KBpleVGKK/9ww8/JDs7m6lTp3Lrrbeydu3aMoIhkH5GDsm4cYpJMm76GTlV/9BxJHQzYNQ2ByYgtovqvtVfB+UK2bf667AuqrL1Rytk+NePwoxLrbJGkwBEIhzuB/YAq4Abgc+owxvgIHhyOnDEQ8dxn3PgiCfomMC0nAA92jQklMoMh6Htgek6mzZtypIlSyJK15nBYQwUSoGBIoPDFR6fSFzz8mKGTFngH3PTNBkyZQHXvLw4zj2LDSJCZoqTnU16U6KsBFAuZZUzU5xlVwU6ZLgmQYkkE5wJfAjcrJS6Uin1qqrj7hSBqp7Q8BkAQzoFp8i+5ISMsLGVbLXSXoIFh122203TZOrUqXTp0oXZs2fzxBNPsHTpUvr2rTgCqNdr+a52Naw9F/bt7bLdnqiYpklBsZu1Owr8AmLIlAWs3VFAQbE7IVYQ0V7VKKUoKPHw/p5WQblC3t/TioKSMPnJoxwyXKOJFuXaHMSaCScAt2DtbxAR8QJTlFKTaql/McHlsuIh9ZJc3k2eiAOFF+Fq1wSWq87MWhf8ZP7Jr4V88sBn5D1xcVC9PZHkqzQAv7eSXTZNkzVr1jBq1CgWL17M4MGDeemll/jd734XUT/tyX+vygq6vl1OdOFgGAbndGnBtgNHgrzCGqY6OKdLi7i74l7z8mIOHi5m9h1nYBgGpmly4bPzaNQglZk3nlata9oODN9v2MfyHaXZ77JbZvLQxV3KrhxiEDJco4kGFf067wAGAH2UUk2VUk2AU4EBInJnbXQuVuTl5QEwyTHNv7PZgWKSozTB3VLnKDYkX8dS5yh/3d69e4Ouc+CAFSvnHOMnoPTJ/hzjJ5THxdixY+nZsye5ubnMmDGDL774ImLBALBhwwYAukle0PXtst2eqCileOHr3zhUHCzEDhV7eeHr3+Lqz2+aJss27efX3Ue48Nl5fsHw6+4jLNu0v9orCKUUQ19czNodwc4Oa3cUMPTFxeE/c5u+MOhvWjBoEoqKhMOfgGuVUhvtCqXUBuAPvrY6y5jp6wHoYmwKqrfLS52jaOo4jAg0dRz2C4hXXvk+6Ph//WslAC3YH1Sfm7eT7a/dwvPPP1+jdJ3/eNPqTxt2BdXbZbs9UfF4POHSYFiIrz1OKKXITLUcBn7dfYQOY2fz6+4jAGSmJtVIcNnnhqos67g2VnOUUZFwSFJK7Q2tVErtASp200lwbrI2Q5eZt+xyE8NSK9lzuV0eMCAj6Phevay/mWIlxttfpPjzx0Vc8uYBUCYPPji5Ruk6T21n/U2V4EnULtvtiYrD4aBhanhX3oapjrju0XA4HCwdN5jG6cFf5cbpSSwdN7hGfevVtrHf/dgO6NdLchMmwZFGEwkVCYfwiQoqb0t4zjijT1Suc+ml5wDWE+F/V7vp8q9CXl/p5t7+ybS84QUefPCeGl3/nHMqDttcWXu8ERGSyglPm+R0xN2fX0Q4JiNYOByTUbMNaHaCozHH7wyK+Dvm+J0JkeBIo4mUioRDdxHJD/MqAE6Kxs1F5G4RUSLSLKDuARH5TUR+FZHzo3GfULZt21at8zZtKgwq7969G8+hXVz83yMMf7+Idg2FH0c14PHBqRhJqezevbtG/azMppDoNgeAC7LDb6Yvr762CLQxBBJog6gud5zTiY0ZPXFjubK6cbIxo2fcExxpNFWhXOGglHIopbLCvDKVUjVWK4lIG+BcYHNAXTYwDOgKXABMFZGo6x7unbm98oPC8GpAjDSPx8PIkQ+z/T83s2Czl2fOT2HxnxvQ41iHf/v4yy//XKN+vvpNxTr5ytrjjcfj4Y3vw+fqfuP7rXG3Oazfa6kDG6cn+TZCWl/r9XuLapQJbtKstTzycyZvnzgFx+AHefvEKTzyc6YOqqepU8TTl/AZ4F6CQ3FcBryjlCrxGcJ/w4rrFFVG+WwFLhUsd0LLoYzuav1dvnw5p556Kl999RqpbU/mT6P+zO2npmD4VAbjPH8GoE8NtVdntKpZe7xJSkrC6SgvSZJUusM8ljgcDlpkpQTZGGwbRIuslGrbHAJjJt0w7Brk9L9xw7BrEibBkUYTKVVLNhAlRORSYJtSamXIj6UVsCSgvNVXF+4ao7GixVY5XeYFFwzi3uXzKZAGpFCaNKdAGlRwFjRtKtx9990888wzNG/enHfeeYd7lzfgMxEecGdxoeMHZnv78o5pBay95JILq9SvUK64ojv/fGElimDjuQpoT2SUUvyhX7uwkUf/0C982tXaZOH9g/F6vX5BYAuImhrKEznBkUYTKTETDiLyJXBsmKZxwFjgvHCnhakLuw5XSr0CvALQu3fvKq3V7c1X73lP56/OWf7NZe95Ty/3nKL1Sxk5bSq7d+9m9OjRPPHEEzRs2JCHV8+hyG3yjjnYLxQA0pJqviiz9d4eHCRTulfAgyOoPVGpC5FHQwVBtDyoEjHBkUZTFWImHJRS54SrF5GTgOMBe9XQGlguIn2xVgptAg5vDVTPQFABaWnWDuYnvdcBcIHjRz739vGX95sNaOo4jFKwq9BkzBwPu9c8TOvWrfnuu+8YNGiQ/RnJSkuiyF3C9ae15eFLu/Hwx6uZsXhzVDxTGjVqBMAssx9DHQv9QmyW2S+oPZHRT9EaTd2k1tVKSqlVgN9VRUTygN5Kqb0i8jHwtog8DRwHdAJ+iEU/BGtJ8qT3Or9QsOnteZUf1V/4aNVB7vmihIMuBw0HXMvX0ybSqVOpx4mIMKxPWw4eKeHhS7shIjx8aTcAGqWn1HgCTElJAeBvnjEAXJGxlg8Ks/1luz3R0U/RGk3dIy42h/JQSq0RkXeBtYAHGKOUinoAIYfDQcuGqewtLMYVcPVkA1wmuPdvo+Pc5hRs2klGm2yOO+9Wkpq1oWPHjmWuFe7J2BYUNSVw8h/y0AcYKSkMKSnhbxO+LNOu0Wg00STuSYiVUu0Dd2IrpR5VSnVUSp2glJodo3tyfrdjcXkJyhVd4nbjXfY/drx2C8bBzbz88ssc2LCS445tw6nHNyk3UFwsn4zznriYXyee4xcEKSkp/DrxnDJBADUajSaaJNTKobYIZyg9t+kh/vXuPezZ/FuZdJ3fT7owrhFEQ1cIesWg0WhizVEpHKBUHWSn63zxxRdp1aoVH3/8MZdccknQsfEOLa3RaDS1zVErHACWLVvGZZddxo4dO7j11lt55JFHKs3KptFoNEcDR7Vw6NChA127duWDDz6oNCubRqPRHE0c1cKhSZMmzJ07N97d0Gg0moRDK9M1Go1GUwYtHDQajUZTBi0cNBqNRlMGLRw0Go1GUwYtHDQajUZTBi0cNBqNRlMGLRw0Go1GUwYtHDQajUZTBqkPCc9FZA+wqZZv2wzYW+lRRzd6jCpHj1HF6PGpnJqMUTulVPNwDfVCOMQDEVmqlOod734kMnqMKkePUcXo8amcWI2RVitpNBqNpgxaOGg0Go2mDFo4VJ9X4t2BOoAeo8rRY1QxenwqJyZjpG0OGo1GoymDXjloNBqNpgxaOGg0Go2mDFo4VBMRuVtElIg0C6h7QER+E5FfReT8ePYvnojIZBH5WURWiMhcETkuoO2oHyMReUpEfvGN0Qci0iig7agfHwARuUpE1oiIKSK9Q9r0GAEicoFvDH4TkfujfgOllH5V8QW0AeZgbbxr5qvLBlYCKcDxwHrAEe++xml8sgLe3wa8pMcoaHzOA5y+938H/q7Hp8wYdQFOAL4FegfU6zGyxsHh++wdgGTfmGRH8x565VA9ngHuBQKt+ZcB7yilSpRSG4HfgKMyMbVSKj+g2IDScdJjBCil5iqlPL7iEqC1770eHx9KqRyl1K9hmvQYWfQFflNKbVBKuYB3sMYmamjhUEVE5FJgm1JqZUhTK2BLQHmrr+6oREQeFZEtwHBgvK9aj1FZbgBm+97r8akcPUYWMR8HZzQvVl8QkS+BY8M0jQPGYqkFypwWpq7e+glXNEZKqY+UUuOAcSLyAHALMIGjaIwqGx/fMeMAD/CWfVqY4+vl+EBkYxTutDB19XaMKiDm46CFQxiUUueEqxeRk7D0nCtFBCx1wHIR6YsludsEHN4a2B7jrsaN8sYoDG8Dn2IJh6NmjCobHxG5HhgCDFY+JTJH0fhAlb5DgRxVY1QBMR8HrVaqAkqpVUqpY5RS7ZVS7bH+Qb2UUjuBj4FhIpIiIscDnYAf4tjduCEinQKKlwK/+N7rMcLyMgHuAy5VSh0JaNLjUzl6jCx+BDqJyPEikgwMwxqbqKFXDlFCKbVGRN4F1mKpCsYopbxx7la8eEJETgBMLI+um0CPUQAvYHnbfOFbgS5RSt2kx6cUEbkCmAI0Bz4VkRVKqfP1GFkopTwicguW16QDeE0ptSaa99DhMzQajUZTBq1W0mg0Gk0ZtHDQaDQaTRm0cNBoNBpNGbRw0Gg0Gk0ZtHDQaDQaTRm0cNAkBCJyhS/K7YkRHHuHiKTX4F4jROSFSOtjjYg0EpGbK2gvrI37ishxIvK/WNxLU/fQwkGTKFwLLMDazFMZdwDVFg4JSCOgXOFQW/dVSm1XSl0Zh35oEhAtHDRxR0QygAHAnwkQDiLiEJF/iMgqX+6DW0XkNuA44BsR+cZ3XGHAOVeKyHTf+0tE5HsR+UlEvhSRFlXo03QReVFEvhGRDSJyhoi8JiI59vXte4vIP0VkuYh8JSLNffWjRORHEVkpIv9nr3REpIUvh8NK36s/8ATQ0Zf/4qkI+9dDRJYE5IRo7Kv/ne+zrvT1qaOIZPj6ttw3lnb0zqD7ikh7EVntu06qiEzzHf+TiJzlqx8hIu+LyOcisk5Enox0TDV1jHjHJdcv/QL+APzH934RVkgSgL8C/0dp7oMmvr95+PJo+MqFAe+vBKb73jemdKPnX4B/+t6PAF4I0w9/PTAdKwyyYIVCzgdOwnqgWgb08B2ngOG+9+MDzm8acN1HgFt972cCd/jeO4CGQHtgdQXjUxim7mfgDN/7ScCzvvffA1f43qdirbCc+HJsAM2wwlxL6H0Dy8DfgGm+9ycCm33XGwFs8PU7FWsHfJt4f4f0K/ovvXLQJALXYk3E+P5e63t/DlaiIA+AUmp/Fa/bGpgjIquAe4CuVTz/E2XNjquAXcqKrWUCa7AmUrBChMz0vX8TGOh7301E5vvuPTzg3mcDL/o+j1cpdaiKfUJEGgKNlFLzfFUzgNNFJBNopZT6wHf9YmXFbhLgMRH5GfgSK7RzZauogcAbvuv8giUEOvvavlJKHVJKFWOFsWhX1c+gSXx0bCVNXBGRplgTZjcRUVhP00pE7sWa1CKJ7xJ4TGrA+ynA00qpj0XkTODhKnavxPfXDHhvl8v77dh9mQ5crpRaKSIjgDOreO/qEC6MM1jCqTlwilLKLSJ5BI9TVa4FwWPhRc8j9RK9ctDEmyuB15VS7ZQV7bYNsBHryXUucJOIOAFEpInvnAIgM+Aau0Ski4gYwBUB9Q2Bbb7318eo/4bvMwBch2VUx9e/HSKShDU523yFpS6zbSpZlP08FeJbbRwQkUG+qj8C85SVgW+riFzuu36Kz9bRENjtEwxnUfqkX9F9v7P7LSKdgbZAuMxsmnqKFg6aeHMt8EFI3f9hTbT/xtJ1/ywiK311AK8As22DNHA/MAv4GtgRcJ2HgfdEZD6wNya9h8NAVxFZhrUCmuSrfwhL//8FpSHLAW4HzvKpm5YBXZVS+4CFIrK6HIN0uohsDXjdhSXsnvKpinoE3PePwG2++kVYyXTeAnqLyFKsCf8XgEruOxVw+Po5ExihlCpBc9Sgo7JqNDVARAqVUhnx7odGE230ykGj0Wg0ZdArB41Go9GUQa8cNBqNRlMGLRw0Go1GUwYtHDQajUZTBi0cNBqNRlMGLRw0Go1GU4b/BzDta1ZRVCkOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f1_org = NNvae.XY_train[:,-1].numpy()\n",
    "f1_pred = dec_output[:,-1].detach().numpy()\n",
    "f1_org2 = NNvae.XY_test[:,-1].numpy()\n",
    "f1_pred2 = dec_output2[:,-1].detach().numpy()\n",
    "\n",
    "f1_org = (f1_org + 1)*y_norm/2\n",
    "f1_pred = (f1_pred + 1)*y_norm/2\n",
    "f1_org2 = (f1_org2 + 1)*y_norm/2\n",
    "f1_pred2 = (f1_pred2 + 1)*y_norm/2\n",
    "plt.scatter(f1_org, f1_pred, marker='x', label='Training Set')\n",
    "plt.scatter(f1_org2, f1_pred2, marker='.', label='Test Set')\n",
    "plt.plot(f1_org, f1_org, color='k')\n",
    "plt.xlabel('Actual Impact Location')\n",
    "plt.ylabel('Decoder Predicted Impact Location')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cef94e6",
   "metadata": {},
   "source": [
    "We can see that the model is somewhat poorly replicating the impact location at the extreme points while the performance is more reasonable for the intermediate locations. We can more closely inspect the distribution of the predicted impact locations from the decoder. Remember, in the normalized dataset, the location $y = -41.9$ represents at normalized value of 1 while $y = 0$ represents a normalized value of -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "id": "4e9953a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "XY = NNvae.XY_train\n",
    "output = dec_output\n",
    "\n",
    "indices = XY[:,-1]==-1\n",
    "data0 = XY[indices]\n",
    "sim0 = output[indices]\n",
    "\n",
    "indices = XY[:,-1]==1\n",
    "data1 = XY[indices]\n",
    "sim1 = output[indices]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "id": "fb1a8bbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'PCstart')"
      ]
     },
     "execution_count": 280,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAV3UlEQVR4nO3df5BV5Z3n8fdX7dC4EEQhDtIqjIUORAbiEHU0UUfdwEYLoykTZscZU2NiZYrND5NoIGOyZFOkmEoqMVatThnNak2cOIxRJOo6449YMT8mClEjSBhxZLGFVSQLA/6KMN/94x6sFvrH7e57+3Y/vF9VXffcc55zzveh4XMPzz33uZGZSJLKclCrC5AkNZ7hLkkFMtwlqUCGuyQVyHCXpAId0uoCACZMmJBTpkxpdRmSNKKsXr365cyc2N22YRHuU6ZMYdWqVa0uQ5JGlIj4Pz1tc1hGkgpkuEtSgQx3SSrQsBhzl6SBevPNN+ns7OT1119vdSlN097eTkdHB21tbXXvY7hLGtE6OzsZO3YsU6ZMISJaXU7DZSbbtm2js7OTqVOn1r2fwzKSRrTXX3+dI444oshgB4gIjjjiiH7/z8RwlzTilRrsew2kf4a7JBXIMXdJRTl92UO8sP21hh1v8mGj+dmis/u1z8c//nE+97nPMWPGjIbV0V+Gu9SHRodFvQYSKoIXtr/GxmXnNex4Uxbd0+99brzxxoadf6AMd6kPjQ6Leg0kVDT0XnnlFT7ykY/Q2dnJnj17+PKXv8z111/PN7/5TebMmcOYMWNYuHAhDzzwAOPHj+frX/86V111FZs2beKaa65h/vz5TanLMXdJGoT77ruPo446iieffJI1a9Ywb968t21/5ZVXOOuss1i9ejVjx47l6quv5v777+fOO+/kK1/5StPqMtwlaRBmzpzJAw88wBe/+EUeeeQRxo0b97bt73jHO94K/JkzZ3LmmWfS1tbGzJkz2bhxY9PqclhGkgbh+OOPZ/Xq1dx7770sXryYD3zgA2/b3tbW9tatjAcddBCjRo16a3n37t1Nq8twl6RB2Lx5M4cffjiXXHIJY8aM4eabb251SYDhLqkwkw8b3dA3oycfNrrX7U899RRXXnklBx10EG1tbVx//fV84QtfaNj5B8pwl1SUob59dO7cucydO/dt6x5++OG3lnft2vXW8pIlS97Wruu2RvMNVUkqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgb4WUVJZvz4Qdmxp3vHHHwBVP1d18yZIljBkzpsd73VesWMHxxx/f9OmADXdJZdmxCZbsaNzxlozru00/rFixgvPPP7/p4e6wjCQN0tKlSznhhBM499xzWb9+PQDf/e53ee9738usWbP48Ic/zKuvvsrPf/5zVq5cyZVXXsns2bN59tlnu23XCIa7JA3C6tWrue2223j88ce54447eOyxxwC46KKLeOyxx3jyySeZPn06N910E6eddhrz58/nG9/4Bk888QTHHXdct+0awWEZSRqERx55hAsvvJBDDz0U4K0v31izZg1XX30127dvZ9euXftNUbBXve36y3CXpEHaO6VvVx/72MdYsWIFs2bN4uabb37bfDMDaddfDstI0iCcccYZ3Hnnnbz22mvs3LmTH/3oRwDs3LmTSZMm8eabb3Lrrbe+1X7s2LHs3Lnzrec9tRssr9wllWXcMY29w2XcMb1uPumkk/joRz/K7NmzOfbYY3n/+98PwNe+9jVOOeUUjj32WGbOnPlWoC9YsIBPfOITXHvttdx+++09thusyMyGHGgw5syZk6tWrWp1GVK3piy6p2VfkN2K844069atY/r06a0uo+m662dErM7MOd21d1hGkgpUd7hHxMER8XhE3F09Pzwi7o+IZ6rH8V3aLo6IDRGxPiIa89avJKlu/bly/wywrsvzRcCDmTkNeLB6TkTMABYA7wbmAddFxMGNKVeS9jcchpebaSD9qyvcI6IDOA+4scvqC4BbquVbgA91WX9bZr6Rmc8BG4CT+12ZJNWhvb2dbdu2FRvwmcm2bdtob2/v13713i1zDXAVMLbLuiMzc0t18i0R8a5q/WTgX7q066zWvU1EXA5cDnDMMb2/Gy1JPeno6KCzs5OtW7e2upSmaW9vp6Ojo1/79BnuEXE+8FJmro6Is+o45v5388N+L6mZeQNwA9TulqnjuJK0n7a2NqZOndrqMoadeq7cTwfmR8QHgXbgnRHxfeDFiJhUXbVPAl6q2ncCR3fZvwPY3MiiJUm963PMPTMXZ2ZHZk6h9kbpQ5l5CbASuLRqdilwV7W8ElgQEaMiYiowDXi04ZVLkno0mE+oLgOWR8RlwCbgYoDMXBsRy4Gngd3AwszcM+hKJUl161e4Z+bDwMPV8jbgnB7aLQWWDrI2SdIA+QlVSSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBWoz3CPiPaIeDQinoyItRHx1Wr94RFxf0Q8Uz2O77LP4ojYEBHrI2JuMzsgSdpfPVfubwBnZ+YsYDYwLyJOBRYBD2bmNODB6jkRMQNYALwbmAdcFxEHN6F2SVIP+gz3rNlVPW2rfhK4ALilWn8L8KFq+QLgtsx8IzOfAzYAJzeyaElS7+oac4+IgyPiCeAl4P7M/CVwZGZuAage31U1nww832X3zmrdvse8PCJWRcSqrVu3DqILkqR91RXumbknM2cDHcDJEXFiL82ju0N0c8wbMnNOZs6ZOHFiXcVKkurTr7tlMnM78DC1sfQXI2ISQPX4UtWsEzi6y24dwObBFipJql89d8tMjIjDquXRwLnAb4CVwKVVs0uBu6rllcCCiBgVEVOBacCjDa5bktSLQ+poMwm4pbrj5SBgeWbeHRG/AJZHxGXAJuBigMxcGxHLgaeB3cDCzNzTnPIlSd3pM9wz89fAe7pZvw04p4d9lgJLB12dJGlA/ISqJBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDml1AdJw99NRn4Yl/7UF550AnDfk51UZDHepDx3xMizZMfTnXTJuyM+pcjgsI0kFMtwlqUCGuyQVyHCXpAIZ7pJUoD7DPSKOjogfR8S6iFgbEZ+p1h8eEfdHxDPV4/gu+yyOiA0RsT4i5jazA5Kk/dVz5b4b+HxmTgdOBRZGxAxgEfBgZk4DHqyeU21bALwbmAdcFxEHN6N4SVL3+gz3zNySmb+qlncC64DJwAXALVWzW4APVcsXALdl5huZ+RywATi5wXVLknrRrzH3iJgCvAf4JXBkZm6B2gsA8K6q2WTg+S67dVbr9j3W5RGxKiJWbd26dQClS5J6Une4R8QY4IfAZzPz33tr2s263G9F5g2ZOScz50ycOLHeMiRJdagr3COijVqw35qZd1SrX4yISdX2ScBL1fpO4Oguu3cAmxtTriSpHvXcLRPATcC6zPxWl00rgUur5UuBu7qsXxARoyJiKjANeLRxJUuS+lLPxGGnA38OPBURT1TrvgQsA5ZHxGXAJuBigMxcGxHLgaep3WmzMDP3NLpwSVLP+gz3zPwp3Y+jA5zTwz5LgaWDqEuSNAh+QlWSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFaie+dwltUBnTqBjybjWnHzcMXDFU605txrCcJeGqfe9cS0bl53XmpO36kVFDeOwjCQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDnltHI8O2ZsGNTS07dmRPoaMmZpYEz3DUy7NgES3a05NTvW3QPG1tyZmngHJaRpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAfYZ7RHwvIl6KiDVd1h0eEfdHxDPV4/gu2xZHxIaIWB8Rc5tVuCSpZ/Vcud8MzNtn3SLgwcycBjxYPSciZgALgHdX+1wXEQc3rFpJUl36DPfM/Anw231WXwDcUi3fAnyoy/rbMvONzHwO2ACc3JhSJUn1GuiY+5GZuQWgenxXtX4y8HyXdp3Vuv1ExOURsSoiVm3dunWAZUiSutPoN1Sjm3XZXcPMvCEz52TmnIkTJza4DEk6sA104rAXI2JSZm6JiEnAS9X6TuDoLu06gM2DKVBSC4w7BpaMa815r3hq6M9boIGG+0rgUmBZ9XhXl/V/HxHfAo4CpgGPDrZISUOsVQHbiheUQvUZ7hHxA+AsYEJEdAL/nVqoL4+Iy4BNwMUAmbk2IpYDTwO7gYWZuadJtUuSetBnuGfmn/aw6Zwe2i8Flg6mKEnS4PgJVUkqkOEuSQUy3CWpQIa7JBXIcJekAg30PndJTTb5sNFMWXRPy879s0Vnt+TcagzDXRqmWhmurXpRUeM4LCNJBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIuWXUP9+eCTs2Df15xx0z9OeURjDDXf2zYxMs2dHqKiT1wWEZSSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoH8EJOk/Uw+bDRTFt0z5Ofd2D7kpyyW4S5pPz9bdHZrTrykNactkeEu6YB3+rKHeGH7ay059+TDRjflxdRwH4laNXkXOIGXivTC9tfYuOy8lpy7WcNfhvtI5ORdKlRnTqBjybghP+9PR00AWhPuzWK4Sxo23vfGtS25gm7FC0qzeSukJBXIK3dJw0arbsH8RftEJrXo6r1ZQ0JNC/eImAd8BzgYuDEzlzXrXJLK0LJbMFs43t6sIaGmhHtEHAz8T+A/A53AYxGxMjOfbsb5WsavnBsyrb5VTRppmnXlfjKwITP/DSAibgMuAJoT7q0MWe9aGRKtvFVNGomaFe6Tgee7PO8ETunaICIuBy6vnu6KiPXV8gTg5SbV1WBr4HPRiAONoD43TL/7HH/TpEqGjr/nA0P/+/zVAefIsT1taFa4d1dpvu1J5g3ADfvtGLEqM+c0qa5hyT4fGOzzgWG49LlZt0J2Akd3ed4BbG7SuSRJ+2hWuD8GTIuIqRHxDmABsLJJ55Ik7aMpwzKZuTsi/hvwT9RuhfxeZq6tc/f9hmoOAPb5wGCfDwzDos+RmX23kiSNKE4/IEkFMtwlqUAtD/eI+EZE/CYifh0Rd0bEYd20OToifhwR6yJibUR8pgWlNkw9fa7azYuI9RGxISIWDXGZDRURF1e/u/+IiB5vE4uIK6p2ayLiBxExYr94rR99Piwibq/+TqyLiD8eyjobqd4+V20PjojHI+LuoaqvGerpcysyrOXhDtwPnJiZfwj8K7C4mza7gc9n5nTgVGBhRMwYwhobrc8+d5nC4b8AM4A/HeF9XgNcBPykpwYRMRn4NDAnM0+k9mb8gqEpryn67HPlO8B9mfkHwCxgXbMLa6J6+wzwGUZ2X/eqp89DnmEtD/fM/OfM3F09/Rdq98Tv22ZLZv6qWt5J7S/E5KGrsrHq6TNdpnDIzN8Be6dwGJEyc11mru+7JYcAoyPiEOBQRvDnI+rpc0S8EzgDuKna53eZuX0IymuKen/PEdFBbbauG5tfVXPV0+dWZFjLw30ffwn8794aRMQU4D3AL4eioCHQU5+7m8JhxL6g1SMzXwC+CWwCtgA7MvOfW1tV0/0+sBX4X9UQxY0R8Z9aXdQQuAa4CviPFtcx5IYqw4ZkPveIeAD4vW42/XVm3lW1+Wtq/3W5tZfjjAF+CHw2M/+9GbU2SgP63OcUDsNNPX3uY//x1P53MhXYDvxjRFySmd9vaKENNNg+U/s3eBLwqcz8ZUR8B1gEfLmBZTZUA37P5wMvZebqiDirweU1RQN+z3uPM2QZNiThnpnn9rY9Ii4FzgfOyR5uvI+INmp/KLdm5h2Nr7KxGtDnETeFQ199rsO5wHOZuRUgIu4ATgOGbbg3oM+dQGdm7r2Ku51auA9bDejz6cD8iPgg0A68MyK+n5mXDL665mhAn4c8w1o+LFN9qccXgfmZ+WoPbYLamOS6zPzWUNbXDPX0mQNzCodNwKkRcWj1Oz+HMt5w61Fm/l/g+Yg4oVp1Ds2aGnuYyMzFmdmRmVOo/b1+aDgHeyO0JMMys6U/wAZqY8tPVD9/W60/Cri3Wn4ftSGJX3dp98FW197MPlfPP0jtbppnqf33r+W1D6LPF1K7Sn0DeBH4px76/FXgN9TuQPg7YFSrax+CPs8GVlV/v1cA41tde7P73KX9WcDdra672X1uRYY5/YAkFajlwzKSpMYz3CWpQIa7JBXIcJekAhnuktRgEfG9iHgpItY06Hh7IuKJ6qeuW6INdxWtyz+KNRHxjxFxaLX+9yLitoh4NiKejoh7I+L4Xo7zpQGe/7N7z6kDys3AvAYe77XMnF39zK9nB8Ndpdv7j+JE4HfAJ6sPlNwJPJyZx2XmDOBLwJG9HKff4V7N7PlZahOg6QCSmT8Bftt1XUQcFxH3RcTqiHgkIv6gmTUMyfQD0jDxCPCHwJ8Ab2bm3+7dkJlPAETEJOAfgHdS+/fxV9RmLxwdEU8AazPzzyJiBbXpIdqB72TmDdX+u4BvAXOBe6h9kOXHEfFyZv7JEPRRw9cNwCcz85mIOAW4Dji7zn3bI2IVtbmolmXmir528ENMKlpE7MrMMdUUwj8E7gPagKmZeUU37T8PtGfm0urK+9DM3Ln3OF3aHZ6Zv42I0dSmijgzM7dFRAIfzczlVbuN1Oanf7npndWwUs3+eHdmnlhNGLYV6Do18KjMnB4RFwH/o5tDvJCZc6tjHZWZmyPi94GHqM1J9Wxv5/fKXaXbe8UNtSv3m4BP9tL+MeB71SRPK/Ze0Xfj0xFxYbV8NDAN2AbsofYiInV1ELA9M2fvuyFrk4j1OpFYZm6uHv8tIh6mNmVwr+HumLtK1/WNqE9l7YtP1gJ/1F3jaqz0DOAF4O8i4i/2bVNNU3su8MeZOQt4nNrwDMDrmbmn8d3QSJa16X2fi4iLoTaRWETMqmffiBgfEaOq5QnUZtXsc3I5w10HooeAURHxib0rIuK9EXFmRBxLba7x71K7yj+pavJmdTUPMA74f5n5avWm2Km9nGsnMLbxXdBwFhE/AH4BnBARnRFxGfBnwGUR8SS1C4x6v1ltOrCq2u/H1Mbc+wx3x9xVtH3HyrusP4ratwH9EfA6sJHanS2nAVcCbwK7gL/IzOci4m+A+cCvqH171gpq34y1HpgILMnMh7sZm/8UsBDY4huqGkqGuyQVyGEZSSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK9P8By7U4xwZRJhYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "feat_to_plot = 1\n",
    "_data = data0[:,feat_to_plot].numpy()*x_std[0,feat_to_plot] + x_mean[0,feat_to_plot]\n",
    "_sim = sim0[:,feat_to_plot].numpy()*x_std[0,feat_to_plot] + x_mean[0,feat_to_plot]\n",
    "\n",
    "bins = np.arange(min(_data), max(_data) + (max(_data) - min(_data))/10, (max(_data) - min(_data))/10)\n",
    "plt.hist(_sim, bins = bins, histtype='step', label='sim') # The simulated distribution of feature 1 for y = -1\n",
    "plt.hist(_data, bins = bins, histtype='step', label='data') # The actual distribution of feature 1 for y = -1\n",
    "plt.legend()\n",
    "plt.xlabel(features[feat_to_plot])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "id": "a66e9a9f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'PCstart')"
      ]
     },
     "execution_count": 281,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWT0lEQVR4nO3df7CcVZ3n8fcXEkgwIQkksCEx3EgBA2sGZOKAMPwQGGHE4oeKMiMjuiily6LiiARRyY6FFUtrRMqVqQgO2V1WlmUgxIFlB4KsCCNLIsEkRJZf2XghQ0IcMuFHJMHv/tFPUjfhktu3u+/tvn3er6pUdz99nud8T93cTz/3dD+nIzORJJVht3YXIEkaPoa+JBXE0Jekghj6klQQQ1+SCjKq3QUATJ48OXt6etpdhiSNKEuXLn0xM6cMZp+OCP2enh6WLFnS7jIkaUSJiP832H2c3pGkghj6klQQQ1+SCtIRc/qS1KgtW7bQ29vL5s2b213KkBkzZgzTp09n9OjRTR/L0Jc0ovX29jJ+/Hh6enqIiHaX03KZyYYNG+jt7WXmzJlNH8/pHUkj2ubNm9l33327MvABIoJ99923ZX/JGPqSRrxuDfxtWjk+Q1+SCuKcvqSucty8+3jupddadrxpE8fy4JyTB7XPpz71Kb74xS9y+OGHt6yOVjH0R6BW/6cejEZ+AaTh9NxLr7F63hktO17PnDsHvc/111/fsv5bzdAfgVr9n3owGvkFkLrZK6+8wkc+8hF6e3t54403+NrXvsZ1113Hd77zHWbPns24ceO4+OKLuffee5k0aRLf/OY3+fKXv8yaNWu45pprOPPMM4e1Xuf0JakJd999NwcccACPPfYYK1as4PTTT9/h+VdeeYWTTjqJpUuXMn78eL761a9yzz33cPvtt/P1r3992Os19CWpCbNmzeLee+/l8ssv54EHHmDChAk7PL/HHntsfyGYNWsWJ554IqNHj2bWrFmsXr162Ot1ekeSmnDIIYewdOlS7rrrLq644gre97737fD86NGjt3/kcrfddmPPPffcfn/r1q3DXq+hL0lNeP7559lnn304//zzGTduHDfeeGO7S9olQ19SV5k2cWxLP3AwbeLYXT6/fPlyLrvsMnbbbTdGjx7Nddddx5e+9KWW9d9qhr6krjLcHyk+7bTTOO2003bYdv/992+///LLL2+/P3fu3B3a9X1uuPhGriQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIH9mU1F2+Ows2rmnd8SbMgEuX19187ty5jBs37i0/q79w4UIOOeSQti27bOhL6i4b18Dcja073twJA7cZhIULF/KBD3ygbaHv9I4kNenqq6/m0EMP5dRTT+WJJ54A4Ic//CHvfve7OeKII/jQhz7Eq6++ykMPPcSiRYu47LLLOPLII3n66af7bTeUBgz9iPhRRKyLiBV9tu0TEfdExJPV7aQ+z10REU9FxBMRcVr/R5Wk7rB06VJuvvlmHn30UW677TYeeeQRAD74wQ/yyCOP8Nhjj3HYYYdxww03cOyxx3LmmWfy7W9/m2XLlnHQQQf1224o1XOmfyNw+k7b5gCLM/NgYHH1mIg4HDgP+LfVPj+IiN1bVq0kdZgHHniAc845h7322ou99957+5eirFixguOPP55Zs2Zx0003sXLlyn73r7ddqwwY+pn5M+C3O20+C1hQ3V8AnN1n+82Z+bvMfBZ4Cvjj1pQqSZ1p29LJfX3iE5/g+9//PsuXL+eqq65i8+bN/e5bb7tWaXROf//MXAtQ3e5XbZ8G/KZPu95q25tExEURsSQilqxfv77BMiSpvU444QRuv/12XnvtNTZt2sRPfvITADZt2sTUqVPZsmULN9100/b248ePZ9OmTdsfv1W7odLqT++8+eUOsr+GmTkfmA8we/bsfttI0qBNmNHaT9xMmLHLp4866ig++tGPcuSRR3LggQdy/PHHA/CNb3yDo48+mgMPPJBZs2ZtD/rzzjuPT3/601x77bXceuutb9luqDQa+i9ExNTMXBsRU4F11fZe4O192k0Hnm+mQEkalEF8pr5VrrzySq688so3bf/sZz/7pm3HHXccjz/++A5t+ms3VBqd3lkEXFDdvwC4o8/28yJiz4iYCRwM/J/mSpQktcqAZ/oR8WPgJGByRPQCVwHzgFsi4kJgDXAuQGaujIhbgMeBrcDFmfnGENUuSRqkAUM/M//8LZ465S3aXw1c3UxRkjQYmdnvJ2i6RWbr3vb0ilxJI9qYMWPYsGFDS4Oxk2QmGzZsYMyYMS05nmvvSBrRpk+fTm9vL9380e8xY8Ywffr0lhzL0Jc0oo0ePZqZM2e2u4wRw+kdSSqIoS9JBXF6RxpIq7+Uo16D/PIOqR6GvjSQVn8pR71a/OUdEji9I0lFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcS1d6RONWFG+9bfcbG3rmXoS52qnaHrYm9dy+kdSSqIoS9JBXF6R4MybeJYeubc2ZZ+H5xz8rD3K3UbQ1+D0q7gbccLjdSNnN6RpIIY+pJUEENfkgpi6EtSQZoK/Yi4NCJWRsSKiPhxRIyJiH0i4p6IeLK6ndSqYiVJzWk49CNiGvA5YHZmvhPYHTgPmAMszsyDgcXVY0lSB2h2emcUMDYiRgF7Ac8DZwELqucXAGc32YckqUUaDv3MfA74DrAGWAtszMx/BPbPzLVVm7XAfq0oVJLUvGamdyZRO6ufCRwAvC0izh/E/hdFxJKIWLJ+/fpGy5AkDUIz0zunAs9m5vrM3ALcBhwLvBARUwGq23X97ZyZ8zNzdmbOnjJlShNlSJLq1UzorwGOiYi9IiKAU4BVwCLggqrNBcAdzZUoSWqVhtfeycyHI+JW4JfAVuBRYD4wDrglIi6k9sJwbisKlSQ1r6kF1zLzKuCqnTb/jtpZvySpw3hFriQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVJCm1t6RhsvP9/wczP2L9nQ+YUZ7+pWGgKGvEWF6vAhzN7a7DGnEc3pHkgpi6EtSQQx9SSqIoS9JBfGNXI0YPXPubEu/0yaO5cE5J7elb6nVDH2NGKvnndGWftv1YiMNBad3JKkghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBWkqdCPiIkRcWtE/DoiVkXEeyJin4i4JyKerG4ntapYSVJzmj3T/x5wd2b+AXAEsAqYAyzOzIOBxdVjSVIHaDj0I2Jv4ATgBoDMfD0zXwLOAhZUzRYAZzdXoiSpVZo5038HsB74u4h4NCKuj4i3Aftn5lqA6na//naOiIsiYklELFm/fn0TZUiS6tVM6I8CjgKuy8x3Aa8wiKmczJyfmbMzc/aUKVOaKEOSVK9mQr8X6M3Mh6vHt1J7EXghIqYCVLfrmitRktQqDYd+Zv4z8JuIOLTadArwOLAIuKDadgFwR1MVSpJaptlvzroEuCki9gCeAT5J7YXkloi4EFgDnNtkH5KkFmkq9DNzGTC7n6dOaea4kqSh4RW5klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVpNmLs6SuN23iWHrm3NmWfh+cc/Kw96vuZuhLA2hX8LbjhUbdz9CX9GYTZsDcCe3p99Llw99vQQx9SW/WruBtxwtNYXwjV5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgrSdOhHxO4R8WhE/EP1eJ+IuCcinqxuJzVfpiSpFVpxpv95YFWfx3OAxZl5MLC4eixJ6gBNhX5ETAfOAK7vs/ksYEF1fwFwdjN9SJJap9kz/WuALwO/77Nt/8xcC1Dd7tffjhFxUUQsiYgl69evb7IMSVI9Gg79iPgAsC4zlzayf2bOz8zZmTl7ypQpjZYhSRqEUU3sexxwZkS8HxgD7B0R/xV4ISKmZubaiJgKrGtFoZKk5jV8pp+ZV2Tm9MzsAc4D7svM84FFwAVVswuAO5quUpLUEkPxOf15wJ9GxJPAn1aPJUkdoJnpne0y837g/ur+BuCUVhxXktRaXpErSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUkJassqmCfHcWbFwz/P1OmDH8fUpdyNDX4GxcA3M3trsKSQ1yekeSCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0JakgDYd+RLw9In4aEasiYmVEfL7avk9E3BMRT1a3k1pXriSpGc2c6W8F/iozDwOOAS6OiMOBOcDizDwYWFw9liR1gIZDPzPXZuYvq/ubgFXANOAsYEHVbAFwdpM1SpJapCVz+hHRA7wLeBjYPzPXQu2FAdjvLfa5KCKWRMSS9evXt6IMSdIAmg79iBgH/D3whcz813r3y8z5mTk7M2dPmTKl2TIkSXVo6usSI2I0tcC/KTNvqza/EBFTM3NtREwF1jVbZKc6bt59PPfSa8Pe77SJY4e9T0ndoeHQj4gAbgBWZebf9HlqEXABMK+6vaOpCjvYcy+9xup5Z7S7DEmqWzNn+scBfwksj4hl1bavUAv7WyLiQmANcG5TFUqSWqbh0M/MnwPxFk+f0uhxJUlDxytyJakghr4kFaSpT+9IGjrTJo6lZ86dbev7wTknt6VvDS1DX+pQ7Qzddr3YaOg5vSNJBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXEtXckdY4JM2DuhPb0e+ny4e+3DQx9SZ2jXcHbjheaNnF6R5IKYuhLUkGc3hmJvjsLNq5pT98TZrSnX0ktYeiPRBvXwNyN7a5C0gjUHaHfpjPfn+85GThj2PuVpEZ1R+i36cx3ekHv+EvqDr6RK0kF6Y4z/TbpzcntOdv3zVSptdp1Udi2vofx+gRDvwl/8rtrWT3POX1pxGvn1bjD/GLj9I4kFcTQl6SCdM30Ts+cO4e9z2kTxw57n5LUjCEL/Yg4HfgesDtwfWbOG6q+AOfWJakOQzK9ExG7A/8J+DPgcODPI+LwoehLklS/oZrT/2Pgqcx8JjNfB24GzhqiviRJdRqq6Z1pwG/6PO4Fju7bICIuAi6qHr4cEU8M4viTgRd32PIfY/BVdp43j2vk68YxQXeOa4cxxbfaWElrdf7PqrH8mgwcONidhir0+xtB7vAgcz4wv6GDRyzJzNmN7NvJunFc3Tgm6M5xdeOYoOvH1TPY/YZqeqcXeHufx9OB54eoL0lSnYYq9B8BDo6ImRGxB3AesGiI+pIk1WlIpncyc2tE/Afgf1H7yOaPMnNlC7toaFpoBOjGcXXjmKA7x9WNYwLHtYPIzIFbSZK6gsswSFJBDH1JKkhHh35ErI6I5RGxLCKW9PN8RMS1EfFURPwqIo5qR52DUceYPlaN5VcR8VBEHNGOOgdroHH1affuiHgjIj48nPU1qp5xRcRJ1fMrI+J/D3eNg1XH/8EJEfGTiHisGtMn21HnYEXExIi4NSJ+HRGrIuI9Oz0/EvNioDENPi8ys2P/AauBybt4/v3A/6R2XcAxwMPtrrkFYzoWmFTd/7ORMKZ6xlW12R24D7gL+HC7a27Rz2si8Dgwo3q8X7trbsGYvgJ8q7o/BfgtsEe7665jXAuAT1X39wAm7vT8SMyLgcY06Lzo6DP9OpwF/Oes+QUwMSKmtruoZmTmQ5n5L9XDX1C7xqFbXAL8PbCu3YW00F8At2XmGoDM7IaxJTA+IgIYRy30t7a3pF2LiL2BE4AbADLz9cx8aadmIyov6hlTI3nR6aGfwD9GxNJq2Yad9bfcw7RhqaxxA42prwupnZmMBLscV0RMA84B/nbYK2vOQD+vQ4BJEXF/1ebjw1xfIwYa0/eBw6hdULkc+Hxm/n44C2zAO4D1wN9FxKMRcX1EvG2nNiMtL+oZU1915UWnr6d/XGY+HxH7AfdExK8z82d9nh9wuYcONNCYAIiI91L7If7JsFfYmIHGdQ1weWa+UTuBHDEGGtco4I+AU4CxwD9FxC8y8/+2o9g6DTSm04BlwMnAQVWbBzLzX9tQa71GAUcBl2TmwxHxPWAO8LU+bUZaXtQzJmBwedHRZ/qZ+Xx1uw64ndrqnX2NuOUe6hgTEfGHwPXAWZm5YXgrbEwd45oN3BwRq4EPAz+IiLOHs8ZG1Pl/8O7MfCUzXwR+BnT0m+91jOmT1KasMjOfAp4F/mB4qxy0XqA3Mx+uHt9KLTB3bjOS8qKeMQ06Lzo29CPibRExftt94H3Aip2aLQI+Xr0rfwywMTPXDnOpdatnTBExA7gN+MsOP1vcrp5xZebMzOzJ2gJRtwL/PjMXDnetg1Hn/8E7gOMjYlRE7EVtNdlVw1tp/eoc0xpqf7kQEfsDhwLPDGedg5WZ/wz8JiIOrTadQu0N9r5GVF7UM6ZG8qKTp3f2B26vpgJGAf8tM++OiM8AZObfUvsUyPuBp4BXqZ2hdLJ6xvR1YF9qZ8IAW7PzVwisZ1wj0YDjysxVEXE38Cvg99S+JW7nEO0k9fysvgHcGBHLqU2JXF79FdPpLgFuitp6X88AnxzheQEDj2nQeeEyDJJUkI6d3pEktZ6hL0kFMfQlqSCGviQVxNCXpBaLiB9FxLqIaMknuaK2SOGy6l9T30Jo6Kur9fllWRER/6P6LD0R8W8i4uaIeDoiHo+IuyLikF0c5ysN9v+FbX2qKDcCp7fweK9l5pHVvzObOZChr2637ZflncDrwGeqhcRuB+7PzIMy83BqK0vuv4vjDDr0I2J34AuAoV+YalmL3/bdFhEHRcTd1ZpHD0REW65y7uSLs6RWewD4Q+C9wJa+F41l5jKAatXF/w7sTe3347PAGcDYiFgGrMzMj0XEQmqX9I8BvpeZ86v9Xwb+htr6NXcCBwA/jYgXM/O9wzBGda75wGcy88mIOBr4AbX1jeoxJmrffbAVmNfM1exenKWuFhEvZ+a4iBhFbVnnu4HRwMzMvLSf9n8FjMnMq6sz9b0yc9O24/Rpt09m/jYixgKPACdm5oaISOCjmXlL1W41MHuEXNGqFoqIHuAfMvOdETGO2oqZT/RpsmdmHhYRHwT+up9DPJeZp1XHOqBaJO8d1L6T4pTMfLqRujzTV7fbdoYOtTP9G4DP7KL9I8CPImI0sHDbXwD9+FxEnFPdfztwMLABeIPai4vU127AS5l55M5PZOZt1NbPeUt9Fsl7JiLuB94FNBT6zumr2/V9A+ySzHwdWEltOeQ3qeZiTwCeA/5L9LM+fkScBJwKvCczjwAepTbNA7A5M99o/TA0klXLUj8bEefC9q9urGs11oiYFBF7VvcnA8fx5sXk6mboq0T3AXtGxKe3bYjad/eeGBEHAusy84fU/irYtpTtlursH2AC8C+Z+Wr1Ztwxu+hrEzC+9UNQJ4uIHwP/BBwaEb0RcSHwMeDCiHiM2onHWXUe7jBgSbXfT6nN6Tcc+s7pq6vtPBffZ/sB1L7Y5Y+AzdS+N/YL1L5z9DJgC/Ay8PHMfDYivgWcCfwS+HfAQmrfuvQEte+RnZuZ9/cz938JcDGw1jdy1QkMfUkqiNM7klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQV5P8D5WK9a28lhDMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "feat_to_plot = 1\n",
    "_data = data1[:,feat_to_plot].numpy()*x_std[0,feat_to_plot] + x_mean[0,feat_to_plot]\n",
    "_sim = sim1[:,feat_to_plot].numpy()*x_std[0,feat_to_plot] + x_mean[0,feat_to_plot]\n",
    "\n",
    "bins = np.arange(min(_data), max(_data) + (max(_data) - min(_data))/10, (max(_data) - min(_data))/10)\n",
    "plt.hist(_sim, bins = bins, histtype='step', label='sim') # The simulated distribution of feature 1 for y = -1\n",
    "plt.hist(_data, bins = bins, histtype='step', label='data') # The actual distribution of feature 1 for y = -1\n",
    "plt.legend()\n",
    "plt.xlabel(features[feat_to_plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2775f4b0",
   "metadata": {},
   "source": [
    "We can see that the model clearly has some bias in predicting the individual feature distributions both for $y=-1$ and $y=1$. The decoder-constructed distribution of \\texttt{PCstart} for $y = -1 (1)$, corresponding to impact location of 0 (-41.9), is rather right (left)-shifted compared to what we see in data. This trend also reflects the model's poorer reconstruction of the impact parameter for these extreme values.\n",
    "\n",
    "*Need to explore some additional models to get it right. Let's see*\n",
    "\n",
    "Next, let's try to look at the latent space. The latent space is expected to learn about the distribution of the input features. Hence, it can show strong correlation with the input features. Let's see how those distributions look like for the different impact locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "id": "c9dbc288",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa233542410>"
      ]
     },
     "execution_count": 282,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5yklEQVR4nO3deXhcV3038O9P+2pLtmVb3mRb8hrbMbZIyG6bkB1SeEvJQghP82JCSAtPobzwwgNtoTxtaSlPsxmH0BDCC6WEkIJJszmy4yQOkZN43+RFXiRbY1nWvo3mvH/oDJ2fPJp779xz75wjnc/zKF9HGl0dzdw5uufcs5AQApZlWZa5sjJdAMuyLMsfW5FblmUZzlbklmVZhrMVuWVZluFsRW5ZlmW4nEz80ClTpoi5c+dm4kdblmUZa8eOHeeEEBUjP5+Rinzu3Lmor6/PxI+2LMsyFhE1Jvu87VqxLMsynK3ILcuyDGcrcsuyLMPZityyLMtwtiK3LMsynK3ILcuyDGcrcsuyLMMpq8iJKJuI3iWi36k6pmWNFQPRGF7edxYD0Vimi2KNQSqvyL8AYL/C41nWmLH1UAT3P70DWw9FMl0UawxSUpET0SwAtwL4kYrjWdZYc+3CCmz45Gpcu/Ci2dWW5ZuqK/IfAPgKgFHbjUS0nojqiag+ErFXJdb4kpeTheuXTkNejr0tZann+6wiotsAtAghdqR6nBBioxCiVghRW1Fhr0osSxe2/958Ki4PrgLwESI6DuAXANYR0dMKjmtZVghs/735fFfkQoivCSFmCSHmArgDwGYhxCd9l8yyrFDY/nvzZWQZW8uy9BHvv7fMpbQiF0LUAahTeUzLsiwrNXsL3bIsy3C2IrcsyzKcrcgty7IMZytyy7Isw9mK3LIsy3C2IrcsyzKcrcgtKwR2GrwVJFuRW0YyrWK00+CtIBlVkZv25rWCY1rFaKfBjx+ZqKeMqshNe/NawTGtYrTL2I4fmainjDqrTHvzWpY1/mSinjKqIresONs6s3SVidaXURX55v1nsf6n9di8/2ymi2JlmG2dWYnG+/0zoypyEEDD/7ECpvsbw/Y5W4nGewvNqHfBusXT8MN7VmPdYrt2ctDG6xtD9z9gVnLjvYWmYs/OAiL6AxHtJKK9RPS3KgqWjL0KC894fWOMx+67sfDHa7zXDSp+634A64QQlwJYCeAmIvqAguNaGTRe3xjRmGA5HozX1tdYomLPTiGE6JL/mys/xsy7YCxcrVju5WQRy/FgvLa+xhIll1tElE1E7wFoAfCSEOKtJI9ZT0T1RFQfiZjzl99erYwv65ZMw8Z7arFuyfi5DzMWWl/j/YJLySsnhBgSQqwEMAvAZUS0LMljNgohaoUQtRUV6f3lz8SLZa9WxpegKrXxXtEETacLLuOn6AshLmB48+WbVB43LhMv1li4WknHeK14gvq9dapoxiKdLriMnKJPRBVEVCb/XQjgegAH/B43GZ1erLFuvFY8Qf3e9twNlk4XXJl4rUkIf/cliWgFgJ8AyMbwH4ZfCiH+LtX31NbWivr6el8/1wrWQDSGrYciuHZhhRZvjrB09UWxYcsR3H9dNUoKcjJdHMtiiGiHEKJ25Od9n6lCiF0A3uf3OJZe4lc44832o63YsOUIVs4uG5e/v2Wm8XOpZVkufGD+ZNx/XTU+MH9ypotieRDkPR0T7hfZinyMMeGk01n8inz70Valx7WvS7CCvKfj9djGj1oZi0x7A47Xm5SqBHWjavOBs/jsT3dg84HxM/U/TEHeYPR6bCNHrYx1plWMdnSEP4GNfhCAGP6PNcbZjSU0ZFrFqNMwLOt/XL2gAg+sqcHVC8w4j0wTZIvH67FTvQeDauHbd7tlhWBbQwSP1R3BtgYzWnbGCbLFo/DYQbXwjarIM9FfbVrXiml9+uky7ve0XSuBCnKNHJXHDqqFb1RFnolK1bSulfGynrbXc6GrL4p/fuEguvqiAZcsufG4GNd4lYlzzaiKPBNjfAeiMbx38oIxV359gzEIMZxjmddzYcOWI3j41QZs2HIk5eOCutLX+d6Fca2bJILsI/d60ZDqXLNdKwDqDp7Fo682oO5geFebbisAXRw82wkhcyzzOt77/uuq8eDaGtx/XXXKxwVVIehcWZrWfZhUgF1XXlvlqc4127UCYF9TB2Iyw+K2AtDF+mvm49bllVh/zXxfx9G54gG8vyHycrKwcnaZ4xVxNBpDTAhEFf/eOo8jN637MJkgu668tqZKCnLw5RsXhbpWj1EV+WeuqcatyyvxmWvCq1Qz8aL4seNEG17YewY7TrT5Oo7uV2le31xuf5++wRgEAuiasjc7A2VKF6jtWoG6SsoL3a9MR1o9pxw3XjIdq+eU+zrOWLhKS+T293lh3xmWqug8jlz3P9puBNkFqrIOCOo+n1EVuapKygvTTvLHtx3Bpt3NeHybvxNa55tzgPc3l9vfZ1Z5IUtVth2O4NG6Bmw7rN95NBb+aAfZBeq1DmjrGsDnf/YO2roGLvpaUPMJ9HyXjmLDlgZs2t2MDVsaQvuZpq2Gt7CihGW6dG+JBDX88NJZ5SCZKkVjMUDItJQbjMZw7Fw3BgM4X73WAd94bg827W7GN57bc9HXgroHo2KHoNlE9CoR7SeivUT0BRUFSybebxnm0DrTZuQV5OcgmwgF+f769HVviXi9inTb9C7IzQLRcKqUk50FIkJOtn7XTmNh7kGqytMvryOkvnP7Mty6vBLfuf2irYuRk5OFLCLkKG7pqriDFwXwJSHEO0RUCmAHEb0khNin4NjMmfY+lmHo6o1iSAh09WZmIolXl8+djJuWTcflc/21IHRviXjd+CLe5HZqel8+bzJuXlaJy+ep/b1VvS6BIICG/2Osb922FGc7+vCt25YqP7bX90JuThbmTSlGbpLK+uqaCnxuTTWurtFs+KEQolkI8Y78dyeA/QBm+j1uMt/88FLUVBTjmx9W/2KN5pUDLSx198Trx7BpdzOeeP2Yr+PUHWzBo3UNqDtoxu/txO3oo42vHcWm3c3Y+NpRpT9f1esShHWLp+GH96zGusXmzjp950TbHz9UqztwFo+82oA6l0NHH3rlEB5+tQEPvXLooq8Z0UdORHMxvO3bW0m+tp6I6omoPhJJ75f49zeOoyHSjX9/47ivcnrx3Y8ux63LK/Hdjy4P7Wf6oeqmz+5TFxATw6mjoPrwayqKQTJVuu+qebh1eSXuu2qe0uOq0N0XxbPvnkZ3hpYvUKFvcEjOaB5Sfuw/HGuDkOnGkODJBDQMVVlFTkQlAJ4B8EUhxEUzdoQQG4UQtUKI2oqK9JoVJ1t7WIYhVTNJR24nvjhZUlkKkqkjrxNs3N7s3NPUDiFTpdePRPD87ma8fkS/ew5ffWYXNu1uxlef2ZXpoqRtf3MHhEzV3jvVxtLJylkTQTJHCmoYqpLaiYhyMVyJ/0wI8WsVx0zmyzcuRHlhDr5848KgfsRFTJuir+omZU52NoiGU0ser2zcvo5NF/pYqpKJWcluVZTmszRRkBceq6vKWTq5YVklHv9ULW5YVnnR14Iahqpi1AoBeALAfiHE9/0XaXR//cwutPVG8dchXjnc+f7ZqKkoxp3vnx3az/SjalIhSvOzUTXJ3zjonCxiqRuvNyXddjl98zZ5H0bxTbNPXzkPtVXl+PSV+nWtRLr6WZooyAuPL92wGA+urcGXbljs6vGpZplGY4KlKiquyK8CcA+AdUT0nvy4RcFxL3KpbKpcmqTJEpRvPrcbDZFufPO53aH9TD/uenw72nqjuOvx7b6Oo/NMRAB4/DU58ek1dy0lt11OP3/7JBoi3fj52ydVFPOPnth2FPWNbXhim9qbqCr8w8dW4NbllfiHj63IdFHSJwSEGE7VvHZXPlp3GA+/2oBH6w5f/MWAyqli1Mo2IQQJIVYIIVbKj9+rKNxI75tVxjIMzR39LHWXl0ss06XzTEQAWDpjArJkuvH7XU34zFP1+P2uppSPC6oF1itvwvUGcDPOL9PuAyUTv/gNYv6a1+7KpdMngGSOFNQ4crNeOcoaHupK4RV7Skk+S91dVT2FZdo0H1u8ZtE0PLC2BmsWuRsy98LesxAyU/m73+1DQ6Qbf/c7tdMgTp7vYamThzbL4XKbLx4uZ4qdp4ZHlux0eUPSC6+Tzy6bNxmrq8pxWZJuP23HkYdpz+kLwyMKTl8I7WdWy6nu1T6nvIelrSfKMl1aT2CB9/G437h1CWoqivGNW5ekfNyMsgKWqmTJew1ZOt5zECPSQE1tfSxV8rqy4td/swf1jW34+m8unmXqdZaoW0ZV5F19gyzDcOnMCSx1d6Oc7Xijh1mPyXjtgw6b1zUrfvTaUTREuvEjh4k+NVOLWapy8yXTWOrkvqvno7aqHPdd7W8N+0wKcrKg15FrA0NRlons6ocAXj96nmUYVK1dEhYiAsn0Y+G0EpBMLRGBhn9RVw/fKSc27XSY4PTMO00s1aERqY8n3zyG+sY2PPmmfrNO3frJm41oiHTjJ282Kj+210l2eXLkTF6SETRGzOwMWnEesQxDUH1aQXlR9gW/6NAX7CxrROrF6/DIh+5chZqKYjx056qUj/vHjy1HeWEO/vFjamfy/m53E0udqFoxM5PiMzqDmNnp1d98+BLUVBTjbz58yUVf6+sfXrupr1/tLFo936WjEPJqRoR4VZOJfUL9uHbRFJbp2idnyu0LYKacCl639nr3RBuORLrxrsNaHP/0wkG09UbxTy8cVFHMPzoi91A9ouNeqrIV57Z1o6Pjrd0sVXr0VTmc8NUkwwmTePqtE2iIdOPpt05c9LWdpztYqmJURd7VP8gyDLtPtSMm0wR1ByIs07VwaglL0+2WU+93O0y9j9/QUr2GS0lhHkudvNM4POLjncbwdt5S7fi5TpYqee1mvGTGBGTRcI50qq2XpSpGVeQQxDMESyonyqm/4U1C8qNqSjHLdB0408FSN17H9s6WO/7Mdtj5Z0hO1BhSPGFjubxZvlzDm+a7TrezNFFt1SSWannrZryqegpuXlaZdAiwqsEIIxlVkWeij7yrvx9CpgnWXzM8AmH9Nf5GIAzJnWyGNN3Rxuvd/1+/28RyNOe7+1mqsq+5k6VO7pCTn+4wZBmKZOpla6I+gFbF/jPDrbn9Z9z9oXukbngns0fqLt7JrKQgB1kE5Zu5G1WR9wzyDMOGrcdY6u6xLQ2ob2zDYz63w8vLzWapG693/zfcvRq1VeXYcPfqlI+LiRhLVYbk8YYUH1eFAjnLsMDgmZ1F+dksVVo+Y7hVvnyGu1b50UgHy0Rar34YljmTC1mGoVuOWe8Ocey6H/WN51mm6/5ra3Dr8krcf22NimIp53UceVF+Dj4wfzKKHIaRXugaZKnK6jmTWeqkXt4Arg9gU4awZMvZ3tkBzPpes3gaPr+2BmtcbrwxIBciH0iyILkdfgigrXuAZRimluSx1F2FLGeFz/JuOXgWm3Y3Y4umo3X6ZQXe77Iidzvy4IJcr/yC4k0WdB7M+drhcyxNtHJ2GUuVvO6Wdep8L8tErV39GBICrYpXmtTxvBpVu/zl20NcbjM/L5el7q6QN1iu8LnWyvP7zrLUzX/vOQMh042qyUUsR3Pp7IksVTnR1s1SJ1dWT2ZpovjKhH43VEnmzYZWxMRwuvG/r5nLMtGPtx1jqYpRFXlfjGcYHrt7NWoqivGYQ9+qLjbWHWGZrivml7PUzdQJ+SydvCJ3iH/FYaf4qkklLFU529HLUie5clJVro7rwLiULcueHcDv8MaxcyydFORks0z0/T+7FPk5hO//2aXqCgjDKvJcirEMw9vHhyeSvH3cjP7DBdNLWabr1+80s9TNqbYelk7KivNYjmb7sVaWqnT2DbHUyfHz3SxNtFiO8V4cwJISGz+5GuWFOdj4SXcXc8/ubGKZ6G9/uw/9UYG//a3a1TWNqshbu2Msw/DbXachZJogS87Oy/I5S2+SrPAmOVR8meJ1V5u3j51nOZqcrCGWqsyV4/rn+hzfH4QjkS6WJnr90DmWKv1H/Sm09UbxH/WnXD1+16kOlonik4pUr2Gkas/OHxNRCxFdvG6jQkMjMgwDQ0Msddd4rotluq5fUsHSdHddXsVyNCdaB1iqkpudxVInpfl5LE20Sd7L2RTAPZ0T8r10wuV76oMLprJM1NYbZamKqrPqSQA3KTrWqDKxbHJb9yBL3RXm57JM12vyxs5rLm/whO3SWeUsnTz77imWo5k3pZClKrPLi1jqJN54M3ipFXxo4RSWKsXX83G7rk9zZy/LRF9YV4PivCx8YZ3aYb1KKnIhxFYAga8tm5fNMwzF+cRSdz1yHZoen+vRXFMziaVumi70snTSLB/X7PD4VXMnsVSlbzDKUifxsfhux+TrqLG9j6XSY8thhI1JhhMm866cXfpuklmmn316B7oHYvjs0zvUFRAh9pET0Xoiqiei+kgkvcHw5QVZLMNwuKWHpe6K5NTfIp9TgDcfOMdSN61yCn2ry6n0QnaNCYcusl1yvfJdDuuWe1V38BxLnTS197I00cTCHJYqfaJ2Fion5OMTtbNcPb4oL4tlotaOPpaqhFYjCiE2CiFqhRC1FRXp9bsOyeVrh0JcxjY2FGOpu245kaXb54SWZjlMrlnD4XKA95tGUXmqRx1O+aBmCN5zxRyWOinJz2Jpoguyz/mC4r5nAPjiL99Dc0c/vvjL91w9flDWFYNJ6ozSohyWqhj1yp3rHmIZhv7BGEvdne8cYJmumLxyjWl6k3fr4VaWTm6+ZAbL0SyZPpGlKs/sOM1SJz0DgqWJJsoW6ETFi1EB3i8aJsuRXpOTjPgaiAqWqhhVkU+U/dQTQ+yvXjZrIkvdTSjKZZmu3kGeulk4rZilk+sWVyCLhjOVfc3tLFUpyM1iqZMbl09naaKTcn3vk4rX+QaAvXIN+70Oa9nHnWgfYJmI5FANUjxkQ9Xww58DeBPAIiI6RUT3qTjuxT+HZxiOn+thqbs2uTRkm88lIqvkwmRVIS5Q5oXXmXxrFk3FA2tqsGbRxUPCEkVlCySquCUyUW4oMVHDjSVI8DRRkZyaXxTAFP0h2UUy5LJ7taI4m2WiRbKlt0hxi0/VqJU7hRCVQohcIcQsIcQTKo47Uq9s+vWG2ASsmVbEUncDMZ7pmjelhKVuWjr7WTr5zbsn8PCrDfjNuxdvv5XodGsPS1UOnW1nqZOX5bIFLzssX6CzroEoS5Vy5CbKOUk2U04mOyuLZaLcHGKpin7tvBT6YzzDkJ+Tw1J32SMybZkYtO+B13sBD20+ynI0XUM8VdH5Xktb7xBLE3X2DLBUacm0CSydnOkcZJkoIs/XiM97WCMZVZHnjsgwDMjm1IAho1amlOSyTNdpORTttKZD0u65Yi5LJxvvWYX8HMLGe1alfFypHNpaqniIa7xrXMMuchRl8TRRxyBPlfbJVtQ+l62pVEsWt/f2s1TFqJcuRjzDcKCpg6XucmSTLcdn021SUR5L3Tz55jGWTv715YPojwr868sHUz6uVD5vpYqbvoNyn9nBEPebdSsTq4qqVib/8JYFMMckKkeYRF2ONJlclM0yUX+UpypGVeQzynJYhiEmN+GNKd6MNyhDcnbekM9Zeq83tLDUTZ/sC+1z2Se693Qny9GUFBewVGXmhDyWOskmniaaUpzLUqVznX0snfTI915PkvfgGrnF25rxvNXb6bYoyzC0y8uUdkMuVzrlRKBOnxOCzsmtzs4p3vJMldmTClk6uXZRBcvRRNp7WKrSM0QsdTIoeJqouWOApUqTS/NZOsmRw39ykgwDapGrdbaM5x2CYiMyDBPkbLcJhsx6qywrYpmuVVXlLHWzas4klk7crl/eK2uzXsW1WkdXH0udFOUSSxNNKMhmqVJ8hKvbPSuKCvJYJjrUfIGlKmbUThlUImeKlQQwYywIPQODLNPVH42y1E2LvOvf4vLu/265NvTuJGtEJxoY4qlKJkZcuTUkNwkeSrJZsClI1rIUwA5B8afF7dPT1TXAMlGkZ4ilKkZV5HnEMwzXLJzKUndD8mbakM+bartOdbLUzVTZzJ3qsrn7xL2rkUXDmUpQrb74LRYdb7VoPtLUlVhMsFSpMdLD0snAiEw0tzyfpSpGVeTx1yiA12pUx1o7WepO1TjV6ikFLHWz61Q7SycbthxDTAxnJpTKJn9pAE1/v2bISmWG4solTJHOKEuV5k4uZukkOsQzUZZsMWQpbjkYVZFnYoegXSc6WOouZ0Sm60/fX8VSN15HE+2R62TscblehmoPyI0EHlC8oYAKEXmDMBLAjcKwBNmqKJHDCEuSDCdMpmJCLstEx871sVTFqIo8E03AeXKtkXmarjlykVSzETz4921HWepmcCjK0kmNvJqqcXlVpdrBM10sdRLUQk5hmlWex1KlnfIibqfLizkakYmCGiZpVEWeidlx+872sNSdqrVWInKKXCSIqXIKNMormkaXVzZn5XCvsw7DvuKDk1QPUnrundMsdRIfYarpSFNXOuRw2w6fw26TWbdoMksnJFf1oySr+/XLC49+lxcgbhlVkWdi4kImbrDqYP7UIpa68TqS4M7LZ7McTXxvZNV7JGdi6Kxb1ZPyWJqoT1659Pm9gkli+/E2lk4G5IixgSQjx3r6BEtVjKrI+4Z4hmG6bKpND6DJFoRUzTov7rhsDkvdXD5/Eksnv95xkuVo4vciVd+TLJLNyCINF1tZNL2MpYkWzyxhqdKC6SUsnZzrFSwTDY5IVfQ7qzRz6vwAS92puo/wH39oZKmbdXKDiHUOG0XENbT0shzN+X6eqvTKK8XeAK4Y/fr9vhaWJjrR2sdSpeNy2OFxl8MPi+XEquIkE6yCus+namOJm4joIBE1ENFXVRxTF/k5PMeL4629LHXzWN0Rlk4yMeIpUV4OT0utqBzrF0025s+nqSUFLJ3MmFDAMtHkAmKpiu+KnIiyATwC4GYASwHcSURL/R5XF5NK8lnqLr4mk9+1mQrlDuCFSXYC10GfXD6uz+UycvG9bhXveetaT5SnpVZHv2CpUnlxDksn3YODLBMFtdKkinfpZQAahBBHhRADAH4B4HYFx9XCqlllLHXXPcAzXZUTC1nqZqJcXneiy2V2b1hWyTJsY2HN7/Hq+Plelk7a5V/r9iR/tWeVFbBURcVpNRNA4h2kU/JzDBGtJ6J6IqqPRCIKfmw46g5HWOpOVReC7jc7vY5gOnimg2XYosRTJ/Gl1xUvwR6qwmyeKkXa+1k66Y7yTFQuW/blilv4KiryZC//Re0bIcRGIUStEKK2okLtWrxBKpZnd7HJZ3kaNu+PsNRNldxLtMrlnqLn5JIF5xyWLlA16mekmXKW38wks/0yLb5fgst9E7QU5DID1dNKWPpxsLmDpSoqKvJTABIH584C0KTguFqIrz2l4cYuSd0kR3Hc5HI0x2gONrex1E1WNk8nH3vfLJajmSh3mJmoeKeZ+MWcy4s6y6M2uZpgm+JVBQHgY6tmsHSSepmMYC4VVJytbwNYQETziCgPwB0A/kvBcbVwQW5Ie8GQjWlPyD02T/jca1P3DTUiHf0snfzo9WMsRxOTjcmY4gFipbJzvFTDTnJlG3ZnULec0dkdwMzOp95sZOkkPqEzycTOwEbX+D6rhBBRAA8CeAHAfgC/FELs9XtcXcyRO9DMcbkTTaYJuaGp8Lmx6dJZJSx10yt3o+91uSt9XjbP0XTKGXedimfeNbX2s9RJZVkeSyMF1ScG4PSFfpZOUi2a5XVGsltKLg+EEL8XQiwUQlQLIf5exTGTUbQelDep/rzqKCvGM019Azx18/HaWSydLJs1keVogpqwUSrvwpUGcTfOp9MXBlia6LoFU1iqNH1iLksn8+R9m3lJ7t8EtQGGfu28FDKx+mHThV6WusuW+wRmJ9kv0IuyglyWuvnJG40snfQNCJajCaqboVX23bYG0Ifr11jYWGLr4XMsVZo3uZSlk5i8oRZLcmPt2uopLFWxFbmDiYW5LHW3p6mHZbrO9wyw1E1ZYRZLJ/MnF7EcTdWkApaqXD2/nKVOglpfJkwFcuJaQQAT2M5397N0svtkG8tErx09z1IVoyryTKiUw5kqDdk9RVX3U7ecMdntcuZk2HKys1k62XywheVoWuXmyK2KN0l+82gbS53EG28+G3EZVS731C0PYG/dlo4+lk6mlOSyTDRNbk04zeUWhW4ZVZEHOeh/NIeauljqbkpxNst0XT5vMkvdEGWxdFJZXshyNO0DPFUpySeWOjHtNlAyJ9sGWKo0NDTE0kmTPHmakpxEff2DLFUxqiKPjwAMcyRgsVwtq9iQVbNUjac9K68+zrq8Cgmb1/UvKicWsBzNlJIclqrky9Wy8jVcNSsTy0OrNkm+XpMUv24A0D/E00mqjbZFVhZLVYyqyDNhYDDKUnfxrny/XfpHznWx1E1b9yBLJ26btCR3biHFO7ic7xxkqZNCuWtKocG7p/QPDrFUqUj2uxe57H9PtUl8UNvq2YrcwZTSQpa6i8+H8DsvQsizUCQ7GzWQ6qonmbfkzaW3HG4yRXp5qhIdkTqZIVspMxxaKzoTsotNuOxq8+La6qksnUyVS49OTbIE6YzyIpaqGFWRx9dpT7Jee2BOn+9hqbtZE/NZpisnJ5ulbrz2686Wm2fPdthEO6h5JZUT8lnqpKtvkKWJumS/UFcA/UPvnLrA0klRLs9Ee+Ux9ro8lltGVeRer8JU6B7iqbuznf0s0/XQHSuRn0N46I6VCkql3oqZZSwdiSyeoz1sRKpSlBNjqZMsuWBNltuFazTkduZuOgZkN9uAy+62VDde44PAVA8GM6oi93qDS4UbllSw1F1JYR7LdH3lmV3ojwp85ZldKoql3EnZQjrpsqV0orWL5WiCWje88cIgS53ky31E8zXcT9Stwvxsliqd74qydBLfzS/Zrn6FOTxVMeqVi8gFfiPJFvoNyKTifJa6a5HLtLY4LNfqpOFsN0vdHDzbwdLJNDkPYJrDfICgdnDJEzx1cuxcL0sT3bi0kqVKD6ytZulkclE2SyY+WmU8j1qJPy/Jnp+gvHm0laXupsnhV9N8DsOaW1HEUjc3yZ1+bnK548/Rlh6Wo4kvp+FyWQ3X4mtwKV6LS4mKkmyWJponz9N5AZyvWw5GWDrpkEN/O5IMAa6UN0Ar/e7FOIJRFXn8eQlzuYorqyez1F2W3DIny+3WOaOYLm+WTvd50zQof7luIWqryvGX6xa6evyi6aUsR9Me5alKJm7Uu1VeXMDSRI++epSlSm4XXIubM7mAZaLOgSGWqhhVkWdCrmwC5SpuCgUlWw6/yvY5DCsiu2YiPrtogvLE68dQ39iGJxzWF4/Ll6Nv8h1G4eQLnqqsrCpnqZP27j6WJpor18aZq3iNHMB7HdDSNcAyUaucR9CqeD6BGbWTlIkF8BvlzbRGQ4YfXi5bDpf7bEHovupjZ2+UpZPpctjfdIfhf/0jUpUiOTO4SMMZwudkE/echiszurVXLqGxN4ClNIbkMLkhl8PlVsyeyDLRTUunsVTFV0VORB8nor1EFCOiWlWF0klLVy9L3e1v6mCZrg8tmcZSN3ub2lk62d/cxXI0QQ0/nCH/gMzQcBz5lKIcliaKr53mcg01Txpbe1g6uXL+FJbsWPLCqFHxBZLfK/I9AD4GYKuCsjhat3AyyzB09A6y1N1QLMYyXR1yamhHAFtnqbBQboS70OWGuDcum8ZyNEG1+k609bHUycSiPJYmijd0gmjwtMiVMFtcroj5g5cPs0y0VN6jWepwr8YrXxW5EGK/EOKgqsI42d7YxjIMEwryWOqupz/GMl0VJQUsdXO4pYulkx9vPcpyNJfOnsBSlVNtPSx1kiX7frMMuQ+UTEleLkuVOnoGWToZivFM1Cwn6jX7nLA3UmivHBGtJ6J6IqqPRNwN4xmpJD+LZRhK5YlRGsAJEoTVs8tYpuvAmQ6Wuvl47UyWTirkOiIVDuuJdA9EWaryyF3vQ3FeFh65631Kj6tCmbwSLzP4irxMXnCUBXDhIUSMpZNUrYPz3QMsVXGsEYnoZSLak+Tjdi8/SAixUQhRK4SorahIb5bk/IpSlmHY13yBpe7q5a4k9Ul2J/HithWVLHVz7FwPSyf3XjGX5Wg65FodHYrX7PjeiwfRPRDD914MrQHr2t4T51ma6H2zy1mq1CdXVOxzubLiRLnV0sQkWy4tnjqBpSqOFbkQ4nohxLIkH88pLYkLmVgAf8WcMpa6U7W569NvHmepmwH5phpw+eYalJsCDDpsDnCuo5+lKkcj3Sx1Ep957nIGupYK5fIChQEsM1Ak15AvcrmWfGlRActEe89cYKmKUZ1iR1q6WYahalIRS9196ooqlumaLfe2nO2wx2WmnJR3/U+6vPv/xLZGlqMJaiJUrbxSrA3gitGvqRNyWZpolhw/PiuAceRrFk9l6aTpfDfLRB1yuGyHy2GzbvkdfvhRIjoF4AoAm4joBTXFSq65Y4BlGPY0d7LU3esN51mmq2ZqKUvdxOSY3pjLsb2XzytnOZqm9n6Wquw9085SJ1fI7fyu0HRbPzeeevMES5VqqyaDZLpx47IZLBNVTy1hqYrfUSvPCiFmCSHyhRDThBA3qipYMgumFLMMw9LKUpa6a+3qZ5muObIFMkfTlkiWXDE8y+XK4YWyWVzo0DwOat2RiJzlF0ky2y/T6k9cYGmiOy+bw1KlHNlNmeOyu7Ih0sUy0QcXV7BUxaiulaysIZZh2HnyAkvd3XXZLJbpeqzuCEvdzCgrZOkkWw6ty3YYYvf5dYtYqjJVbjE3VfHu6SqUFuSyNNFvdzWxVCknJwtZRMjJcVddrpTdZyuTdKP9545TLFUxqiLvGeAZhla51kirpmuOjFSQmweS6cd1CytY6qZFjsNtcTket6aikOVocuXUwFzFUwSXz5zIUie1c8pZmsjrBDEvLp87GTctm47L57rrWlkxs5Rlovg8PZ/z9S5iVEVeUpjNMgyfXTOfpe52nroAIdOPEnl1VqLpVdqqqjKWTl7aH2E5miffOMpSlT1NF1jqZMWsCSxN9IdjbSxV2rC1AZt2N2PD1gZXjz8g1/A/kGQt/yXTJ7JUxaiKfFJRPsswnO8aZKm7Y+e6WaZrMBpjqZsfbjnKUpWgVn3slOPSOwPYU9KvvU2dLE00b0oRS5UG5HtgwOV7IdUyGYdaOlmqYlRFfl72qZwPsW/l01fOQ21VOT595bzQfqYf2fKGTLbPceTvnmpjqZvPXDOPpZObL5nGcjTXyIWOrkmy4JEf8ytKWOrkhFw24ISGywe49aElU1mqlC0nrmS7nMDS0TvAMtFK2epZqbj1Y1RF3i0XcOoOcSGnJ98cXvf6yTfdrXudaW4rLCcxOR055nJactheORBh6aQgN0feO0g9auWIHPt7JMkYYD9ukMuW3qB4+VIVJsn9XSf53Oc1k7YcamWp0so55cii4XTjpX0tLBOtkJtTrHC5SYVbRlXkV1ZXsAzDQnkFtVDDK6lkSgpzkU2EkkJ/fdvTS4tY6mbD3atRW1WODXevdvcNRMMzgh2uqipKC1iqUie3CatzuV1YmN5uPM/SRBWleSxVuqp6Cm5eVomrqt210pbLBdeWJ1l47XCkh6UqRlXkJQU5LENBNDxSOcx1AXzweod9NDdeUsFSN1MnFuBXn7sSUx0WwYqLRocQE8OZypc+tADFeVn40ocWqCjmH02VfximKv4DocJl8yaxNNHe050sVXri9WPYtLvZ9W5UH14+k2Wi2XLz79kOm4B7ZVRF7nYImUr7z3RAyDTBv71yCJt2N+PfXjnk6zg6r58NAF19UfzzCwfR5bKb7ZCcnHEoySSNRJ/72Q50D8TwuZ/t8F3GRGfae1jqpFWuxNeqeEW+MN20fBpLle6/rhoPrq3B/ddVu3r8rlPtLBP9akcTS1WMqsg3y6Fjmx2GkKm0QHapLDCka2X78VaW6XpgzQI8uLYGD6xRe2WqyoYtR/Dwqw3YsMXdhKW7L6tCTUUx7r4s9Ro0Xqf+uzWzvJilTmbKSVUzXU6u0tGTrx9nqVJeThZWzi5DnssJQS/ta2aZaIWcR7BC8XwCoyryKbJZOiXE5ukRuUzqEZfLpWZapKOPZbq8nrxh83qV9N3nD6Ah0o3vPn8g5eM+e+18lqoU5Gaz1En8Ndb1tXZjnZzyvk7x1HcA2HzgLD770x3YfOCsq8cvrpzAMtGqqnKQTJWMeuXqT7SyDMMnamehckI+PlHrb8p7WD5zbTXLdG3efxbrf1qPzfvdnbxh8/qH5lu3LUVtVTm+ddvSlI8ryMlmqcq9Vwy3CO71uSplEOZOLmRpoqtqKpAlU7VoNIaYEIi6HEd+07JKlokK8rKRRYSCPLXnl1EV+YlzvSzD8Ff/uRPNHf34q//cGdrP9OMPx86zTBsBNPwfLW09FMH9T+/A1kPuutl2nW7HeycvYNfp1KsP/vfeFpaquG0RZMKz7zazNNGaRdPwwNoarFmkvo/c61ormw+0sEykajDCSEZV5BPyeIbhrvfPZqm7WEywTNfVNRX43JpqXB3AFY4KH5g/GfdfV40PzHf3hrh2YQU2fHI1rnVYO0bV5tUjfef2Zbh1eSW+c/sypcdV4V8+fikqJ+TjXz5+aaaLkra6A2fxyKsNqHPZ/eHFusXT8MN7VmPdYv9/JLxO93fLqIo8mpXNMgwfXDIdty6vxAeXTA/tZ/pRLW/KVvu8ObvtcASP1jVg22H9xj0DwPajrdiw5Qi2H1XbzVYt11+vVrwOe3vvAN5pPI/2JLP9Mu3JNxrR3NGPJ99IvemGznY3DY8u292U+dFlKSfliRGpiFEV+SUzJrAMg9cxpJn2+bULcOvySnx+rb/RJlF5RR/1eWUfFK9X5G67YtZfMx+1VeVYf43am513Pr4dzR39uPPx7UqPq8J7conm9wxZqjmZJdNLQDJV89qNV1KYJyflXdx18BcfXIgH19bgLz64UGkZ/e4Q9D0iOkBEu4joWSIqU1SupP7k0lkgmWG576p5uHV5Je67yoy1Vt461orn9zTjrWM+r1SFgBDDqaO6gy14tK4BdQfd9WW7rfif2t6I+sY2PLVd7dXpE/fWojgvC0/cW6v0uCrc8f4ZLE1UkJsDIuclGNLhtlsuLlW3ZElBDr584yLlkxr9XpG/BGCZEGIFgEMAvua/SKMrzh9+sYrzw5vZueNEG17YewY7Tui5eNRI3f1RCDGcvric0p4pe5s6EBPD6ca2hggeqzuCbQ2pr6q8Dmt067e7zqB7IIbf7jqj9LgqlBUXIJsIZcX6zTp16+oFFXhgTQ2uXqD+nk5eThauXzrN9QipVOfaQDSGl/eddb2Solt+t3p7UQgRrzG2Awj0UjnIF2s0Xv8aZ9pL+85AyPTD6/ZWYVt/zXzcurzSdReI2yFkQV0xXTJjArIo3G5Bt4IaSREmt3+oQyEAMfyfi3gdk+6Wyj7yPwfw/GhfJKL1RFRPRPWRSHpPdiZuwHX3RfHsu6dDXXHRj+vlMp7X+1zO8/J5k3HzskpcrumGvF5bSm6HkHmd+u/WmkVT8cCaGqxZpH6ZVb8eqTuMTbub8Ujd4UwXJX0pKs+wrVsyDRvvqcW6JclvdgZRTseKnIheJqI9ST5uT3jM1wFEAfxstOMIITYKIWqFELUVFeld3UZjMUDIDMk3ntuDTbub8Y3n9oT2M/1IdaPFC91v8nq92el2CNkPXj6Eh19twA9e9rdWzUg6jwI6eb6XpYl0v/CIC6qcjhW5EOJ6IcSyJB/PAQAR3QvgNgB3CxHsnbHokICQGRadx/8mo2rM65+umonKCfn401UXr+CmA6/DD932c55o7WKpjMYTrL59+zLUVpXj24ac48kEeeHhtV871SiXx7cdwabdzXh8m9pNzf2OWrkJwP8B8BEhROCLkRw62wkhMyzlJXl45O5VKC8xd9H9dHz5V7vQ3NGPL/9qV6aLkpTXexdu34y3yGnVtySZXu2HzhOs3M561VlQN6kB78MPV88px42XTMfqJBtRBLW/gd8+8ocBlAJ4iYjeI6INCso0qgfWyhX5fI6RHste3NuM9U/V48W9/qZbP3rXKtRWlePRu1YpKplaXkcSuH0z3nLpTPzoU7W45VK1LZGgJjCpkKriMUVQN6kB7xcNqVoHBfk5yCZCgeKRd76OJoSoUVUQN+IvVpgGojFsPRTBtQsrjFgdbvepdsRk3uajMopv3DBWZHr0UaZ/fiobXzuKTbubMWdSEf7PzYszXRztxC8a3Iq3CpK1DuJdn6rPA/1rpgRBjShIxWuzKtN03uRXJa/9lm6v4IN6vb22IMKk89BIt4Ian52OVK2DoM4D/c6qFB6tO4yHX23AoyEOk/I6OiLTXt53luVY5bXCdftG1/nKOSg6D410S6cLrkz8UTGqIp8/uYhlGLxOBc+02ZOKWI5VXitcnd7outG5/96tIC+4VI5aCYpRFfnhlm6WYfA6FTzTVs0pG96BZE5ZposSKK9NVLcV/4t7mvGZp+rx4h5z1+b2aiy0QoIcp++1Ys7E82lURb581kRkyQzL59fW4MG1Nfj82lDv66ZtzeJp+PzaGqzxOY5cpz7HZILqI3/vVDuEzPFC5/571wIcp++1Yk71fGq51krYbrikEhs/VYsbLlE7xjeVIIc1BUFVM1n3roigypdNPC0zqNz8IUhBnbcU8GTMpGpra0V9fX3oP3c8UDVcUvdhl0GVr6svig1bjuD+66qN+eNtBevlfWdx/9M7sOGTqz0NQ0zG73lLRDuEEBethWxURa575aID+xxZllo6vadGq8iNeqdnYmd33fuKRwpqmUxr7DLtHA+bCfcQ9C1ZMhlYeEj3vuKLaLScZ5CCqnzGY6Vm3DmeRJCvm9dj23HkDjKx8JBpQ7NSroU8hgQ1IWgsVGpemXaOJxPk6+b12Bk5h4QQoX+sXr1apOOlvWdE9dc2iZf2nknr+62xo39wSLy094zoHxxy9Xi3505n76D43n8fEJ29gyqKaYXE6/kQ5LFTPd5vOQHUiyR1qlFX5GPhyiFo46VrIKgJQVptGWa5FmQ/ttdjp3p8UFfrRlXkJtx0yLTx2DXghutzJ6B7DOPlD6yVWlAXo7ZGHGNsq8WfoO4x2D+wFqDp6odE9G0i2iU3lXiRiGaoKphlZUJQbzSd/8Da1oL5/J6t3xNCrBBCrATwOwDf9F8kyw975acnnbsF7TljPl9nlRAicUnAYoz50cv60/nKz9KTPWfCo+2iWUT090R0EsDdSHFFTkTriaieiOojEfuXPyg6X/mZYDx2M9hzJjwZG7VCRC8T0Z4kH7cDgBDi60KI2QB+BuDB0Y4jhNgohKgVQtRWVNi//LobjxUaYLsZrGBlbNSKEOJ6IcSyJB/PjXjo/wPwv5SWzsqY8Vqhmba1n2UWXUetLEj4348AOOCvOJYuxmu/6VjY9swaf/wuuPwPRLQIQAxAI4D7/RfJ0kH8ymG8Ga9/wCyz+arIhRC2K8UaU8brHzDLbPY2tWVZluFsRW5ZlmU4W5FblmWlYMJQXFuRW5ZlpWDCUFxbkVuWZaVgwkgmv8MPLcuyxjQTRjLZK3LLsizD2YrcsizLcLYityzLMpytyC3LsgxnK3LLsizD2YrcsizLcLYityzLMpytyC3LsgxnK3LLsizDKanIiejLRCSIaIqK41mWZVnu+a7IiWg2gA8BOOG/OJZlWZZXKq7I/xXAVwAIBceyLMuyPPK7+fJHAJwWQux08dj1RFRPRPWRiL7LQVqWZZnGcfVDInoZwPQkX/o6gP8L4AY3P0gIsRHARgCora015up9IBrD1kMRXLuwAnk59t6wZVn6cayZhBDXCyGWjfwAcBTAPAA7ieg4gFkA3iGiZJW+sUxYVN6yrPEt7fXIhRC7AUyN/7+szGuFEOcUlEsbJiwqb1nW+GY3lnBgwqLylmWNb8oqciHEXFXHsizLstyzd+8sy7IMZytyy7Isw9mK3LIsy3C2IrcsyzKcrcgty7IMZytyy7Isw5EQ4c+WJ6IIgMY0v30KAB0nHdlyeWPL5Y0tlze6lgvwV7YqIcRFsxMzUpH7QUT1QojaTJdjJFsub2y5vLHl8kbXcgHBlM12rViWZRnOVuSWZVmGM7Ei35jpAozClssbWy5vbLm80bVcQABlM66P3LIsy+JMvCK3LMuyEtiK3LIsy3BaVuRE9HEi2ktEMSIadZgOEd1ERAeJqIGIvprw+UlE9BIRHZZZrqhcjsclokVE9F7CRwcRfVF+7W+I6HTC124Jq1zycceJaLf82fVevz+IchHRbCJ6lYj2y9f8CwlfU/p8jXa+JHydiOjf5Nd3EdEqt98bcLnuluXZRURvENGlCV9L+pqGVK41RNSe8Pp80+33Blyuv04o0x4iGiKiSfJrgTxfRPRjImohoj2jfD3Yc0sIod0HgCUAFgGow/CuQ8kekw3gCID5APIA7ASwVH7tnwB8Vf77qwD+UVG5PB1XlvEMhgfxA8DfAPhyAM+Xq3IBOA5git/fS2W5AFQCWCX/XQrgUMLrqOz5SnW+JDzmFgDPAyAAHwDwltvvDbhcVwIol/++OV6uVK9pSOVaA+B36XxvkOUa8fgPA9gcwvN1LYBVAPaM8vVAzy0tr8iFEPuFEAcdHnYZgAYhxFEhxACAXwC4XX7tdgA/kf/+CYA/UVQ0r8f9IIAjQoh0Z7G65ff3zdjzJYRoFkK8I//dCWA/gJmKfn6iVOdLYnmfEsO2AygjokqX3xtYuYQQbwgh2uT/bsfw/rhB8/M7Z/T5GuFOAD9X9LNHJYTYCuB8iocEem5pWZG7NBPAyYT/P4X/qQCmCSGageGKAgl7i/rk9bh34OKT6EHZtPqxqi4MD+USAF4koh1EtD6N7w+qXAAAIpoL4H0A3kr4tKrnK9X54vQYN98bZLkS3YfhK7u40V7TsMp1BRHtJKLniegSj98bZLlAREUAbgLwTMKng3q+nAR6bmVsz04iehnA9CRf+roQ4jk3h0jyOd9jKVOVy+Nx8gB8BMDXEj79GIBvY7ic3wbwLwD+PMRyXSWEaCKiqQBeIqID8koibQqfrxIMv+G+KITokJ9O+/lK9iOSfG7k+TLaYwI51xx+5sUPJFqL4Yr86oRPK39NPZTrHQx3G3bJ+xe/AbDA5fcGWa64DwN4XQiReKUc1PPlJNBzK2MVuRDiep+HOAVgdsL/zwLQJP99logqhRDNsvnSoqJcROTluDcDeEcIcTbh2H/8NxE9DuB3YZZLCNEks4WInsVws24rMvx8EVEuhivxnwkhfp1w7LSfryRSnS9Oj8lz8b1BlgtEtALAjwDcLIRojX8+xWsaeLkS/uBCCPF7InqUiKa4+d4gy5XgohZxgM+Xk0DPLZO7Vt4GsICI5smr3zsA/Jf82n8BuFf++14Abq7w3fBy3Iv65mRlFvdRAEnvcAdRLiIqJqLS+L8B3JDw8zP2fBERAXgCwH4hxPdHfE3l85XqfEks76fkCIMPAGiXXUJuvjewchHRHAC/BnCPEOJQwudTvaZhlGu6fP1ARJdhuD5pdfO9QZZLlmcigOuQcM4F/Hw5CfbcUn33VsUHht+0pwD0AzgL4AX5+RkAfp/wuFswPMrhCIa7ZOKfnwzgFQCHZU5SVK6kx01SriIMn9ATR3z/TwHsBrBLvliVYZULw3fFd8qPvbo8XxjuJhDyOXlPftwSxPOV7HwBcD+A++W/CcAj8uu7kTBiarRzTdHz5FSuHwFoS3h+6p1e05DK9aD8uTsxfBP2Sh2eL/n/nwbwixHfF9jzheGLtmYAgxiuu+4L89yyU/Qty7IMZ3LXimVZlgVbkVuWZRnPVuSWZVmGsxW5ZVmW4WxFblmWZThbkVuWZRnOVuSWZVmG+/95p5CCuvT0RwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = XY[:,-1].numpy()\n",
    "_z0 = z[:,0].detach().flatten().numpy()\n",
    "_z1 = z[:,1].detach().flatten().numpy()\n",
    "_z2 = z[:,2].detach().flatten().numpy()\n",
    "plt.scatter(y,_z0,label='z0',marker='.', s=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "id": "75c31ebb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa23116a9d0>"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYp0lEQVR4nO3dfZAcdZ3H8c93N9lsNiEBsouEJLi5Yl0QDJCslAJ1B+jlgBI4L8QnfCpTpqhS0SoiArEkV2qVwMlV3WnOC4Q6uPKJ0nCkNJ6sd4IFCLrJkWfFaEBDnjaYLFnI7mZ3v/fHzMr0ZnYedn4z3Z19v6qmvrPTPb/+pqfznZ5fd//a3F0AgPSqizsBAEBlKOQAkHIUcgBIOQo5AKQchRwAUm5SHAttbm721tbWOBYNAKm1cePGQ+7eMvr1WAp5a2ururq64lg0AKSWmb2U73W6VgAg5SjkAJByFHIASDkKOQCkHIUcAFKOQg4AKUchB4CUq7iQm1mjmf3KzDab2XYz+8cQiSXFod5+ffLhLh3q7Y87FQDIK8Qeeb+kq9z9QkkXSbrazN4RoN1EWLV+uzp3HNCq9dvjTgUA8qr4yk7P3JmiN/vn5OzjpLlbxYrF7Zo6uV6fuvKcuFMBgLyC9JGbWb2ZPS/poKROd38uzzzLzazLzLq6u7tDLLYmWpun6d6lF6q1eVrcqQBAXkEKubsPuftFkuZKusTMLsgzzxp373D3jpaWE8Z8AQCMU9CzVtz9iKQnJF0dst04cbATQNKFOGulxcxOzT6fKundkn5TabtJwcFOAEkXYhjb2ZIeMrN6Zb4YHnH3HwVoNxE42Akg6UKctbJF0sUBckmkkYOdAJBUXNkJAClHIS+Cg50Aki6WW72lyR3rtqhzx0FJrvs/+va40wGAE7BHXsTgkEciACQNhbyIW97VptkzG3XLu9riTgUA8qKQF7H2qd3a19OntU/tjjsVAMiLQl7EBWedEokAkDQU8iK+9t8vRCIAJA2FvIimyXWRCABJQ3Uq4vTpDZEIAElDIS9iz+G+SASApKGQA0DKUcgBIOUo5ACQchRyAEg5CjkApByFHABSjkIOAClHIQeAlKu4kJvZPDP7uZntNLPtZvbZEIkBAEoT4g5Bg5JudfdNZnaKpI1m1unuOwK0DQAoouI9cnff5+6bss+PStopaU6l7QIAShO0j9zMWiVdLOm5PNOWm1mXmXV1d3eHXCwATGjBCrmZTZf0Q0mfc/dXR0939zXu3uHuHS0tLaEWCwATXpBCbmaTlSni33b3dSHaBACUJsRZKyZpraSd7n5f5SkBAMoRYo/8MkkfkXSVmT2ffVwboF0AQAkqPv3Q3Z+SZAFyAQCMA1d2AkDKUcgBIOVSVcgP9fbrkw936VBvf9ypAEBihLhEv2ZufWSznnyhWwODw3roE5fEnQ4AJEKq9si3v3wkEgEAKSvkZ5wyJRIBACkr5Dv290YiACBlhRwAcCIKOQCkHIUcAFKOQg4AKUchB4CUo5ADQMpRyAEg5SjkAJByFPIi6i0aASBpKORF3H3jAk1rqNfdNy6IO5WSMEIkMPFQyIu4ov0MXXpOs65oPyPuVEpy6yOb1bnjgG59ZHPcqQCokVQNYxuHz3zn//TLP7yi3r5BfXf5O+JOp6ide3siEcDJL8geuZk9aGYHzWxbiPaSZONLf47EpLvv/RepdVaT7nv/RXGnAqBGQnWt/IekqwO1lSitzdMiMenOnT1DbW86RefOnhF3KgBqJEghd/dfSErHLmuZjrw+EIlJRx85MPFwsLOIWxe/RXWWiWmwdc+RSARw8qtZITez5WbWZWZd3d3dtVpsxe7r/J2GPRPToLd/MBIBnPxqVsjdfY27d7h7R0tLS60WW7HWWdMiMekGhjwSUR7Ow0ca0bVSxO3XnKvZMxt1+zXnxp0KamDV+u3q3HFAq9ZvjzsVoGShTj/8rqRfSmo3sz1mtixEu0mw9qnd2tfTp7VP7Y47FdTAisXtWrporlYsbo87FaBkQS4IcvcPhmgniZZdNl+bXjqsZZfNjzsV1EBr8zTdu/TCuNMAykLXShFrn96tvT19Wvs0e+Q4OXFcIP0o5EXwUzuZKD7hcFwg/SjkRYz81E7LlZ0TRbWKz0T8gmBnJf1SVcibGuojERNXucWn1AL9me9sUueOA/rMdzaFSDMVpjdO0pFjxzW9kTH00ipVhXxK9u4OU2p4l4eJuIeWBuX+Uip1D/7XLx6OxFCSvB3RtZJ+qSrk7WfOiMRamKgbeZILz3iUugc/OOyRGEqSt6MlF8/RlEl1WnLxnLhTwTilqpB/bckCLV00V19bUru79UzU/sMkF57xOHC0T489/7IOHO0rOF/dqBjKssvm66yZjYk8jfW2H25R/+CwbvvhlrhTwTilqpDHceBxoh7sPNm+wD7ywHMaGHJ95IHnCs7X2FAXiaGsfnKX9vb0afWTu4K2G0L/0HAkIn1SVchPtp/71TB75pRIHK+T7QtsMDv2zGCRMWj6jg9HYih7jxyLxCQpdd0guVJVyE+2n/vVsPqmRZo9s1Grb1pUUTtJ/9J8bvcresvKDXpu9yslzT88Ko45n0djKDv2Ho3EJHl9YCgSkT6pKuQn28/9agg1NkzSb1Dx4fuf1cCQ68P3Pxt3KiWZlu2qmRa4ywaQUnbzZcbBKG7F4nZNnVyvT115TkXt/Hr3oUhMmpHu3LR06zY1TFLvwICaGpL3X25aQ51eGxjmSybFkrdVoSKhvuxeP+6RmDSldpUkxcHegUhMktOaGvTaQJ9Oa2qIOxWME1/BSKVL558eiRi/PUf6IhHpQyFHKj2z+8+RmHQzGusjEQiJQg7UwKt9Q5EIhEQhB4CUo5ADQMpRyAEg5SjkAJByQQq5mV1tZr81s11mdnuINgEApam4kJtZvaRvSrpG0lslfdDM3lppuwCA0oTYI79E0i53/4O7D0j6nqQbArQLAChBiEI+R9Kfcv7ek30twsyWm1mXmXV1d3cHWCwAQApTyPPdQPOEATrcfY27d7h7R0tLS4DFAgCkMIV8j6R5OX/PlbQ3QLsAgBKEKOS/ltRmZvPNrEHSByStD9AuAKAEFQ9j6+6DZvZpST+VVC/pQXfnFj4AUCNBxiN39w2SNoRoCwBQHq7sBICUo5ADQMpRyAEg5SjkyKuhPhqTZsbU+kgEJjIKOfIaGIrGpHm9fzgSi2mZNjkSgZMJhRx5TZ9SH4lJMzjskViMm0UiTi6Hevv1yYe7dKi3P+5UYkEhR153L1mgqZPrdPeSBXGnklfbGdMjsZizZk6NxLHUWTQiHVat367OHQe0av3EvISFQo68Hnt+r44dH9ZjzydztIXTmiZHYjH7e45F4lhGdvBL3NFHQixZOEdTJtVpycITxuurWMi9/Wr9cqCQI6/jQ8ORmDQvvvJaJBaz6vrzNXVynVZdf37B+eotGpEOt/1gi/oHh3XbD7YEbzvk3v6d67aqc8cB3blua4DM3kAhR153XXe+li6aq7uuK1z44nLf+y5S66wm3fe+i0qaf8O2/Tp2fFgbtu0vON/8WU2RGMqMxkmRiLAGBocjMaT5s6ZGYiWqtYNEIUderc3TdO/SC9XaPC3uVPK6vK1FT3z+Sl3eVtqQyCsWt2vporlasbi94Hx/PHwsEkM5b/aMSERYQ9m+sKEq9In968//EInFPL5jv+bf/mM9vuPEnYZbrmrTWTMbdctVbUFzpJBjQpjeOElHjh3X9CJ7xPNOmxqJody9ZIGWLpqb2IPHadebPU+2NwHny9788EZ5No628rGt2tvTp5WP0bUClK3Ufs49R45FYiilfpHEoTXbZdAaoOsgLlMmWSTG6S1nTo/EXDv2Ho3EUCjkmBBK7Vr56nvfpmkN9frqe98WdPlJPj3u9exe7OsJ2Jsdr/5Bj8Q4fevDHVq6aK6+9eGOmi2TQo4JodQ+/yvaz9Cl5zTrivYzgi5/yaK5mjLJtGTR3KDthrDquuwZPQk9sJ02cRxfopADOap1etjKdVvVP+haGbjdEL7ftUfHjg/r+1174k7lpBDHVaYUciBHtU4Pu+fGBWqd1aR7bkzewc7tL/dEIipz6yOb1bnjgG59ZHPNlpm8Iy9AjO667nw1T9+lT115TtB2z509Q21vOkXnJvD0w7NnNenQawM6O/C58xPV5j8djsRaoJADOUb6N0O7Y90Wde44KMl1/0ffHrz9SjRlxypuSuqYxSnT2zcYibVQUdeKmS01s+1mNmxmtTtEC6TM4JBHIk5ecYzXU+ke+TZJ/yDp3wPkApy0qtVlE4J7NKIyw6NiLVRUyN19pyQZYzwDBVWryyaEhkl1kYj0qdknZ2bLzazLzLq6u7trtVigLNU6dSzJNz5I+gBpKK5oITezn5nZtjyPG8pZkLuvcfcOd+9oaSltoCOg1qp1BWaSr+xM8vABKE3RT87d312LRIAkWLG4XVMn1wfvy65WuyGMfMlMmVSnb3xoYdzpYBz4CgZyVKsvO8l95En+kkFpKj398L1mtkfSOyX92Mx+GiYtALWS9LHnUVylZ608KunRQLkAAMaB840AIOUo5ACQchRyAEg5CjkAFPChS+ZFYhJRyAGggHWbXo7EJKKQA0ABt19zruotE5OKQg4ABTz9+0Ma8kxMKgo5ABQQcix5GxVD4RJ9ACgg5FjyPiqGQiEHgAJCjpMzc+ok9Rwb1MypYUsvXSsAUCOv9Q9GYigUcgCokcHhaAyFQg4ABYS8u9O8UxsjMRQKOQAUEPLuTgeO9kdiKBzsBIACQt5447zZM7R5T4/Omz0jQGZvYI8cAAoIeeONVdefr9kzG7Xq+rA3uqaQA0CNrH5il/b19Gn1E7uCtkshB5B6dRaNSRXyKtFc9JEDSL3pUybp1b5BTZ+S7JIW8irRXOyRA0i91llNkRhSyNMPq3Wj64oKuZnda2a/MbMtZvaomZ0aKC8AKNmew69HYkghTz984cBRvf0rnXrhwNEAmb2h0j3yTkkXuPsCSS9IuqPylACgPGYWiSGtWNyupYvmasXi9orb+viDv1J374A+/uCvAmT2hooKubs/7u4jgwY8K2lu5SkBQHk+f/W5qrNMDC1kd8g9Ny5Q66wm3XPjggCZvSFkH/knJP1krIlmttzMusysq7u7O+BiAUx0X//pbzXsmZhkl7e16InPX6nL21qCtlv0EK+Z/UzSmXkmrXT3x7LzrJQ0KOnbY7Xj7mskrZGkjo6O0MPxApjAhtwjcaIpukfu7u929wvyPEaK+MckvUfSTe4TdC0CiNVXbrhAUyfX6Ss3XBB3KgWFPAMmV6VnrVwt6QuSrnf38IeLAaAEG7bt17Hjw9qwbX/cqRQU8gyYXJX2kX9D0imSOs3seTP7VoCcAKAsSxbO0ZRJdVqycE7cqRS07LL5Omtmo5ZdNj9ou5WetXKOu89z94uyj5tDJQYApVr56Db1Dw5r5aPb4k6l4HABq5/cpb09fVr9JGOtAEDEF99znqZOrtMX33Ne3Klo2KMxV7XGWqGQA0i9DVuzfeRbk91HfstVbTprZqNuuaotaLsUcgCpF/Lqy0o11EVjrrVP79benj6tfXp30GUme6gwACjByNWXSeBmkjwbo0LebSgXe+QAUq9a52ePx5tPmxqJuRI5+iEAJEG1zs8ejz8dPhaJtUAhB5B6Seoj/8y7zonEXIm8shMAkqBaXRbj8dAzL0Virmr9cuBgJwAEVGgALw52AkAKfOHv2lVnmTgaBzsBIAX+6fEXNOyZWCsUcgCpV617YY5HHGOjU8gBpF617oU5HnGMjU4hB5B61boX5njEMTY6Z60ASL2Re2EmwbLL52vjS4e17PKwY44Xwh45AAS0+old2tfTp9VPhB1zvBAKOYDUS9JYK9Uac7wQulYApN7IFZNTJtXpGx9aGGsud113vpqn7wp+0U8hFHIAqVetKybHI44hdSsq5Gb2ZUk3SBqWdFDSx919b4jEAKBUSRqPPA6V9pHf6+4L3P0iST+S9KXKUwKA8iSpjzwOFRVyd381589pkmrXuw8AWXeu26rOHQd057qtcacSi4r7yM3sq5I+KqlH0pgncprZcknLJenss8+udLEA8BfHh4YjcaIxLzIegJn9TNKZeSatdPfHcua7Q1Kju99VbKEdHR3e1dVVbq4AkNeLh17TN3+eOVMkCWOSV4uZbXT3jhNeL1bIy1jAmyX92N2LDjBAIQeA8o1VyCvqIzeztpw/r5f0m0raA4DxSNLBzjhyqbSP/Gtm1q7M6YcvSbq58pQAoDxJuiAojlwqKuTuviRUIgAwXkm6ICiOXIL1kZeDPnIAKF9V+sgBAFFx9JFTyAEgoJE+8lXrt9dsmQyaBQABxdFHTiEHgIDiGMCLrhUACIg+cgBIOfrIASDl6CMHgJSjjxwAUDYKOQCkHIUcAFKOQg4AKUchB4CUo5ADQMpRyAEg5SjkAJByFHIASDkKOQCkHIUcQOrFMeJgkgQp5Ga2wszczJpDtAcA5YhjxMEkqXjQLDObJ+lvJf2x8nQAoHxxjDiYJCH2yP9Z0m2SPEBbAFC2kREHW5unxZ1KLCoq5GZ2vaSX3X1zCfMuN7MuM+vq7u6uZLEAgBxFu1bM7GeSzswzaaWkOyUtLmVB7r5G0hpJ6ujoYO8dAAIpWsjd/d35Xjezt0maL2mzmUnSXEmbzOwSd98fNEsAwJjGfbDT3bdKOmPkbzN7UVKHux8KkBcAoEScRw4AKRfsnp3u3hqqLQBA6dgjB4CUM/fan0BiZt2SXhrn25slJbEfnrzKQ17lIa/yJDUvqbLc3uzuLaNfjKWQV8LMuty9I+48RiOv8pBXecirPEnNS6pObnStAEDKUcgBIOXSWMjXxJ3AGMirPORVHvIqT1LzkqqQW+r6yAEAUWncIwcA5KCQA0DKJbKQm9lSM9tuZsNmNuZpOmZ2tZn91sx2mdntOa+fbmadZva7bDwtUF5F2zWzdjN7Pufxqpl9LjttlZm9nDPt2lrllZ3vRTPbml12V7nvr0ZeZjbPzH5uZjuzn/lnc6YFXV9jbS85083M/iU7fYuZLSz1vVXO66ZsPlvM7BkzuzBnWt7PtEZ5XWFmPTmfz5dKfW+V8/p8Tk7bzGzIzE7PTqvK+jKzB83soJltG2N6dbctd0/cQ9J5ktolPaHMQFz55qmX9HtJfyWpQdJmSW/NTrtH0u3Z57dLujtQXmW1m81xvzIn8UvSKkkrqrC+SspL0ouSmiv9d4XMS9JsSQuzz0+R9ELO5xhsfRXaXnLmuVbSTySZpHdIeq7U91Y5r0slnZZ9fs1IXoU+0xrldYWkH43nvdXMa9T810n63xqsr7+WtFDStjGmV3XbSuQeubvvdPffFpntEkm73P0P7j4g6XuSbshOu0HSQ9nnD0n6+0CpldvuuyT93t3HexVrqSr998a2vtx9n7tvyj4/KmmnpDmBlp+r0PaSm+/DnvGspFPNbHaJ761aXu7+jLsfzv75rDJDRldbJf/mWNfXKB+U9N1Ayx6Tu/9C0p8LzFLVbSuRhbxEcyT9KefvPXqjALzJ3fdJmUKhnOF2K1Ruux/QiRvRp7M/rR4M1YVRRl4u6XEz22hmy8fx/mrlJUkys1ZJF0t6LuflUOur0PZSbJ5S3lvNvHItU2bPbsRYn2mt8nqnmW02s5+Y2fllvreaecnMmiRdLemHOS9Xa30VU9VtK9joh+WyAncecvfHSmkiz2sVn0tZKK8y22mQdL2kO3Je/jdJX1Ymzy9L+rqkT9Qwr8vcfa+ZnSGp08x+k92TGLeA62u6Mv/hPufur2ZfHvf6yreIPK+N3l7Gmqcq21qRZZ44o9mVyhTyy3NeDv6ZlpHXJmW6DXuzxy/+S1Jbie+tZl4jrpP0tLvn7ilXa30VU9VtK7ZC7mPceagMeyTNy/l7rqS92ecHzGy2u+/L/nw5GCIvMyun3WskbXL3Azlt/+W5md0v6Ue1zMvd92bjQTN7VJmfdb9QzOvLzCYrU8S/7e7rctoe9/rKo9D2UmyehhLeW828ZGYLJD0g6Rp3f2Xk9QKfadXzyvnClbtvMLPVZtZcynurmVeOE34RV3F9FVPVbSvNXSu/ltRmZvOze78fkLQ+O229pI9ln39MUil7+KUop90T+uayxWzEeyXlPcJdjbzMbJqZnTLyXJl7rW4r9f1VzMskrZW0093vGzUt5PoqtL3k5vvR7BkG75DUk+0SKuW9VcvLzM6WtE7SR9z9hZzXC32mtcjrzOznJzO7RJl68kop761mXtl8Zkr6G+Vsc1VeX8VUd9sKffQ2xEOZ/7R7JPVLOiDpp9nXz5K0IWe+a5U5y+H3ynTJjLw+S9L/SPpdNp4eKK+87ebJq0mZDXrmqPf/p6StkrZkP6zZtcpLmaPim7OP7UlZX8p0E3h2nTyffVxbjfWVb3uRdLOkm7PPTdI3s9O3KueMqbG2tUDrqVheD0g6nLN+uop9pjXK69PZ5W5W5iDspUlYX9m/Py7pe6PeV7X1pcxO2z5Jx5WpXctquW1xiT4ApFyau1YAAKKQA0DqUcgBIOUo5ACQchRyAEg5CjkApByFHABS7v8BxYfk7CXuK50AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(y,_z1,label='z1',marker='*',s=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "id": "d7f69a35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa233227ad0>"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOrUlEQVR4nO3daYwkdRnH8d9PlBcqUXFXQRYcjcSIRiPpEBXjHQOLiL4wwROPZOMLEki8lpAYE994RKPGKyuSeBDxhXgEMboexKiB2IucriCYJSysMJ5oTETk8UXXOMVMTx/T/+rqp+v7SSbbU9VV9Wx19a+f/ld1jyNCAIC8HtF2AQCA2RDkAJAcQQ4AyRHkAJAcQQ4AyT2yjY3u2LEjVlZW2tg0AKR14MCBP0XEzo3TWwnylZUV9fv9NjYNAGnZvnPYdIZWACA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghyoWdn7fa3s/X7bZQBTIcgBILlWPtkJLKpDHzmr7RKAqdGRA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDmG4jtHgDwIcgBIju9awVB85wiQBx05ACRHkANAcgQ5ACRHkANAcgQ5UuLySGAdQQ4AyXH5IVLi8khgXbGO3PZRtn9j+8pS6wQAjFdyaOUCSQcLrg8AMIEiQW57l6SzJF1SYn0AgMmV6sg/Jen9kh7a6g6299ju2+6vrq4W2iwAYOYgt/0aSfdFxIFR94uIfRHRi4jezp07Z90sAKBSoiM/XdJrbR+SdLmkV9j+eoH1AgAmMHOQR8RFEbErIlYknSvppxHxlpkrAwBMhA8EAUByRT8QFBFXS7q65DoBAKPRkQNAcgQ5ACRHkGOoZft2wWX7/wB1BDkAJOeImPtGe71e9Pv9qZdb66j45jsAXWT7QET0Nk6nIweA5FJ9HzmdOABsRkeOoTg5CORBkANAcqmGVjA/DGMBedCRA8CcNDVkSZADQHIMrSAlPlOAjJo6XunIMRRXrQB50JEjJTpxYB0dOQAkR0eOoeh4gTzoyDEUY+RAHgQ5ACRHkCMl3jEA6whyAOl1/YU91clOPgSCNRwDWFRt5FSqIAeAYbr+wp4qyLv+YM0T+xrII9UYeRvjYF0de+vq/xvIKFWQA9vFCxOWWaqhlTYwxAB027QnL9vIDIJ8yZQ6Y75sL2DL9v8B6ghyDMWlnsDAtM+BNp47jJGPwdgqgEVHR75kutJB844BWEeQj0FQAIuv6y/sBPkYXT9AFhWPB7COIAeQXpMv7BmaOYIcAAoaFfhNvShw1QqA9Lp+dRkd+RiL/HYKQPMyXEdOkCOlDOOWmJ+uHwczB7ntEyV9VdJxkh6StC8iPj3rehcFgQEsviafpyW/a6WpHCnRkT8o6T0RcZ3tYyQdsL0/In5bYN3AULywYlGlHFqJiCOSjlS3/2H7oKQTJBHkANIrGcgprlqxvSLp+ZKuHTJvj+2+7f7q6mrJzQJApxU72Wn7sZK+JenCiLh/4/yI2CdpnyT1er0otV08HGP6QFkln1NNPS+LdOS2H6VBiF8WEVeUWCcALJumrncvcdWKJX1Z0sGI+OTsJWEWXenEeecBrCvRkZ8u6a2SXmH7+upnd4H1Yhu6/gk3oItKXLXyC0kuUAsAYBv4ZCeA9Lo+xMaXZgE1DE3l1PXHjSAHgOQYWgFquv4WHTkR5AAwQoYXd4ZWAGCEDOPvBDk6IcOTEdgughyoIfCREWPkADAni/yHJYClkeHEFvJq6juCCHJ0AgG9Nb6ALD/GyNEJk459M0aOjOjIgY6jE8+PjhwAkqMjRyfQdWKZ0ZEDQHIEOTqBk5hYZgQ5MAe8kKBJBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJMcflsBQ/CEGIA86cgBIjiDHUHx/NpAHQQ4AyTFGjqEYIwfyoCMHgOQIcgBIrkiQ2z7D9q22b7e9t8Q6AQCTmTnIbR8l6XOSzpR0iqQ32j5l1vUCmA+uUMqvREd+mqTbI+IPEfGApMslnVNgvQCACZQI8hMk3VX7/XA17WFs77Hdt91fXV0tsFkAgFQmyD1kWmyaELEvInoR0du5c2eBzQIApDJBfljSibXfd0m6p8B6AQATKPGBoF9LOtn20yTdLelcSW8qsF4Ac8CHv/KbOcgj4kHb50v6oaSjJF0aEbfMXBmAuVi7YoVAz6vIR/Qj4ipJV5VYFwBgOnyyEwCSI8gBIDmCHACS42tsgY7jJGd+dORAx/FdK/nRkSMlukhgHR05ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAcgQ5ACRHkANAcnyNLVDD1+MiI4IcnTBpQK/9gQUCHZkwtAIAyRHkAJAcQQ4AyTFGDtQwNo6M6MiBGv6iPDKiI0cntH01Cp0+mkSQAzUELjIiyNEJBDSWGWPkQA1j5MiIIEcnENBYZgytADUMwSAjghydQEBjuzIcOwT5kslw0AGZtH3p6iQYI18yjAUD3UNHDtRk6L6w2No4dgjyMXhCA4tvkZ6nbTQDMwW57Y9LOlvSA5LukPSOiPhbgbqAVixSIGByXX8nNesY+X5Jz4mI50q6TdJFs5e0WLKNOR/6yFmdPZiBrpopyCPiRxHxYPXrNZJ2zV4SZpHthQfA7EqOkb9T0je3mml7j6Q9knTSSScV3CwALI6FPNlp+8eSjhsy6+KI+G51n4slPSjpsq3WExH7JO2TpF6vF9uqFgCwydggj4hXjZpv+zxJr5H0yoggoAEslWk77FEnXpvq1me9auUMSR+Q9NKI+FeZkjALTnQCZWW4ImbWq1Y+K+kYSfttX2/7iwVqwgw42Ql0z0wdeUQ8o1QhALBdi9wtz0OqT3a28WB1/QABMmhy+CNDBvClWQCQXKqOvA0ZTnQAXdfk87NkBjSVJ6mCnFAdj30DdE+qIAeAeSvZHC3kdeTzRrcJYNGNyqmmRhU42QkAyaXqyNsYI+ddAIBSmsoTOvIlwyc7ge5J1ZEDKI93nfPD5YeYCE9KdFHXL01OFeRdfZCAJnU9BEtL9zW2XcBBDiy+rj8/UwU5oQqUx/OprDauI08V5G3gIAew6NzGX2fr9XrR7/fnvl0AyMz2gYjobZye6jpyrpEGgM1SBTkAYLNUY+SMVwPIjC/NEkMrADBMqo4cABYdHwgag6EVANgsVZBjPD40BbSrjedeqjFyAMBmdORLhk4c6B46cgBILlWQc/khAGyWKsgBAJulGiNn/BcANqMjB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASI4gB4DkCHIASM4RMf+N2quS7tzm4jsk/algOaVQ13SoazrUNZ1FrUuarbanRsTOjRNbCfJZ2O5HRK/tOjairulQ13SoazqLWpfUTG0MrQBAcgQ5ACSXMcj3tV3AFqhrOtQ1HeqazqLWJTVQW7oxcgDAw2XsyAEANQQ5ACS3kEFu+w22b7H9kO0tL9OxfYbtW23fbntvbfqxtvfb/n317xMK1TV2vbafafv62s/9ti+s5n3I9t21ebvnVVd1v0O2b6q23Z92+Sbqsn2i7Z/ZPlg95hfU5hXdX1sdL7X5tv2Zav6Ntk+ddNmG63pzVc+Ntn9l+3m1eUMf0znV9TLbf689Ph+cdNmG63pfraabbf/X9rHVvEb2l+1Lbd9n++Yt5jd7bEXEwv1IepakZ0q6WlJvi/scJekOSU+XdLSkGySdUs37mKS91e29kj5aqK6p1lvV+EcNLuKXpA9Jem8D+2uiuiQdkrRj1v9XybokHS/p1Or2MZJuqz2OxfbXqOOldp/dkn4gyZJeIOnaSZdtuK4XSXpCdfvMtbpGPaZzqutlkq7czrJN1rXh/mdL+ukc9tdLJJ0q6eYt5jd6bC1kRx4RByPi1jF3O03S7RHxh4h4QNLlks6p5p0j6SvV7a9Iel2h0qZd7ysl3RER2/0U66Rm/f+2tr8i4khEXFfd/oekg5JOKLT9ulHHS73er8bANZIeb/v4CZdtrK6I+FVE/LX69RpJuwpte6a6Glq29LrfKOkbhba9pYj4uaS/jLhLo8fWQgb5hE6QdFft98NaD4AnR8QRaRAUkp5UaJvTrvdcbT6Izq/eWl1aaghjirpC0o9sH7C9ZxvLN1WXJMn2iqTnS7q2NrnU/hp1vIy7zyTLNllX3bs06OzWbPWYzquuF9q+wfYPbD97ymWbrEu2Hy3pDEnfqk1uan+N0+ix1dofX7b9Y0nHDZl1cUR8d5JVDJk287WUo+qacj1HS3qtpItqk78g6cMa1PlhSZ+Q9M451nV6RNxj+0mS9tv+XdVJbFvB/fVYDZ5wF0bE/dXkbe+vYZsYMm3j8bLVfRo51sZsc/Md7ZdrEOQvrk0u/phOUdd1Ggwb/rM6f/EdSSdPuGyTda05W9IvI6LeKTe1v8Zp9NhqLcgj4lUzruKwpBNrv++SdE91+17bx0fEkerty30l6rI9zXrPlHRdRNxbW/f/b9v+kqQr51lXRNxT/Xuf7W9r8Lbu52p5f9l+lAYhfllEXFFb97b31xCjjpdx9zl6gmWbrEu2nyvpEklnRsSf16aPeEwbr6v2gquIuMr2523vmGTZJuuq2fSOuMH9NU6jx1bmoZVfSzrZ9tOq7vdcSd+r5n1P0nnV7fMkTdLhT2Ka9W4am6vCbM3rJQ09w91EXbYfY/uYtduSXl3bfmv7y7YlfVnSwYj45IZ5JffXqOOlXu/bqisMXiDp79WQ0CTLNlaX7ZMkXSHprRFxW236qMd0HnUdVz1+sn2aBnny50mWbbKuqp7HSXqpasdcw/trnGaPrdJnb0v8aPCkPSzp35LulfTDavpTJF1Vu99uDa5yuEODIZm16U+U9BNJv6/+PbZQXUPXO6SuR2twQD9uw/Jfk3STpBurB+v4edWlwVnxG6qfWxZlf2kwTBDVPrm++tndxP4adrxIerekd1e3Lelz1fybVLtiaqtjrdB+GlfXJZL+Wts//XGP6ZzqOr/a7g0anIR90SLsr+r3t0u6fMNyje0vDZq2I5L+o0F2vWuexxYf0QeA5DIPrQAARJADQHoEOQAkR5ADQHIEOQAkR5ADQHIEOQAk9z86kjt367PXqwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(y,_z2,label='z2',marker='+', s=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c8db3fe",
   "metadata": {},
   "source": [
    "As you can see, some of these latent spaces actually show a strong correlation with impact location while others are rather randomly distributed, showing almost no correlation. These patterns will vary every time you run the code and retrain the model, since there is no *a priori* way of telling how these correlations will build up during model training.\n",
    "\n",
    "You may have noted that the model's performance is not so great. Can you think of ways to improve it? Try changing the hyperparameters (number of layers, number of nodes per layer, number of epochs etc.) Also, there are modified architectures for VAE, like the so-called $\\beta$-VAE [(link to paper)](https://openreview.net/forum?id=Sy2fzU9gl). The $\\beta$-VAE architecture is expected to disentangle the dimensions of the latent space representation. Hence, it may allow a single latent space dimension to more strongly correlate with the impact location compared what you are seeing with the plain VAE architecture"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}