gplvm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Process Latent Variable Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Gaussian Process Latent Variable Model](https://en.wikipedia.org/wiki/Nonlinear_dimensionality_reduction#Gaussian_process_latent_variable_models) (GPLVM) is a dimensionality reduction method that uses a Gaussian process to learn a low-dimensional representation of (potentially) high-dimensional data. In the typical setting of Gaussian process regression, where we are given inputs $X$ and outputs $y$, we choose a kernel and learn hyperparameters that best describe the mapping from $X$ to $y$. In the GPLVM, we are not given $X$: we are only given $y$. So we need to learn $X$ along with the kernel hyperparameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not do maximum likelihood inference on $X$. Instead, we set a Gaussian prior for $X$ and learn the mean and variance of the approximate (gaussian) posterior $q(X|y)$. In this notebook, we show how this can be done using the `pyro.contrib.gp` module. In particular we reproduce a result described in [2]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import torch\n",
    "from torch.nn import Parameter\n",
    "\n",
    "import pyro\n",
    "import pyro.contrib.gp as gp\n",
    "import pyro.distributions as dist\n",
    "import pyro.ops.stats as stats\n",
    "\n",
    "smoke_test = ('CI' in os.environ)  # ignore; used to check code integrity in the Pyro repo\n",
    "# assert pyro.__version__.startswith('1.5.0')\n",
    "pyro.enable_validation(True)       # can help with debugging\n",
    "pyro.set_rng_seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data we are going to use consists of [single-cell](https://en.wikipedia.org/wiki/Single-cell_analysis) [qPCR](https://en.wikipedia.org/wiki/Real-time_polymerase_chain_reaction) data for 48 genes obtained from mice (Guo *et al.*, [1]). This data is available at the [Open Data Science repository](https://github.com/sods/ods). The data contains 48 columns, with each column corresponding to (normalized) measurements of each gene. Cells differentiate during their development and these data were obtained at various stages of development. The various stages are labelled from the 1-cell stage to the 64-cell stage. For the 32-cell stage, the data is further differentiated into 'trophectoderm' (TE) and 'inner cell mass' (ICM). ICM further differentiates into 'epiblast' (EPI) and 'primitive endoderm' (PE) at the 64-cell stage. Each of the rows in the dataset is labelled with one of these stages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data shape: (437, 48)\n",
      "---------------------\n",
      "\n",
      "Data labels: ['1', '2', '4', '8', '16', '32 TE', '32 ICM', '64 PE', '64 TE', '64 EPI']\n",
      "--------------------------------------------------------------------------------------\n",
      "\n",
      "Show a small subset of the data:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Actb</th>\n",
       "      <th>Ahcy</th>\n",
       "      <th>Aqp3</th>\n",
       "      <th>Atp12a</th>\n",
       "      <th>Bmp4</th>\n",
       "      <th>Cdx2</th>\n",
       "      <th>Creb312</th>\n",
       "      <th>Cebpa</th>\n",
       "      <th>Dab2</th>\n",
       "      <th>DppaI</th>\n",
       "      <th>...</th>\n",
       "      <th>Sox2</th>\n",
       "      <th>Sall4</th>\n",
       "      <th>Sox17</th>\n",
       "      <th>Snail</th>\n",
       "      <th>Sox13</th>\n",
       "      <th>Tcfap2a</th>\n",
       "      <th>Tcfap2c</th>\n",
       "      <th>Tcf23</th>\n",
       "      <th>Utf1</th>\n",
       "      <th>Tspan8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.541050</td>\n",
       "      <td>-1.203007</td>\n",
       "      <td>1.030746</td>\n",
       "      <td>1.064808</td>\n",
       "      <td>0.494782</td>\n",
       "      <td>-0.167143</td>\n",
       "      <td>-1.369092</td>\n",
       "      <td>1.083061</td>\n",
       "      <td>0.668057</td>\n",
       "      <td>-1.553758</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.351757</td>\n",
       "      <td>-1.793476</td>\n",
       "      <td>0.783185</td>\n",
       "      <td>-1.408063</td>\n",
       "      <td>-0.031991</td>\n",
       "      <td>-0.351257</td>\n",
       "      <td>-1.078982</td>\n",
       "      <td>0.942981</td>\n",
       "      <td>1.348892</td>\n",
       "      <td>-1.051999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.680832</td>\n",
       "      <td>-1.355306</td>\n",
       "      <td>2.456375</td>\n",
       "      <td>1.234350</td>\n",
       "      <td>0.645494</td>\n",
       "      <td>1.003868</td>\n",
       "      <td>-1.207595</td>\n",
       "      <td>1.208023</td>\n",
       "      <td>0.800388</td>\n",
       "      <td>-1.435306</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.363533</td>\n",
       "      <td>-1.782172</td>\n",
       "      <td>1.532477</td>\n",
       "      <td>-1.361172</td>\n",
       "      <td>-0.501715</td>\n",
       "      <td>1.082362</td>\n",
       "      <td>-0.930112</td>\n",
       "      <td>1.064399</td>\n",
       "      <td>1.469397</td>\n",
       "      <td>-0.996275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.056038</td>\n",
       "      <td>-1.280447</td>\n",
       "      <td>2.046133</td>\n",
       "      <td>1.439795</td>\n",
       "      <td>0.828121</td>\n",
       "      <td>0.983404</td>\n",
       "      <td>-1.460032</td>\n",
       "      <td>1.359447</td>\n",
       "      <td>0.530701</td>\n",
       "      <td>-1.340283</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.296802</td>\n",
       "      <td>-1.567402</td>\n",
       "      <td>3.194157</td>\n",
       "      <td>-1.301777</td>\n",
       "      <td>-0.445219</td>\n",
       "      <td>0.031284</td>\n",
       "      <td>-1.005767</td>\n",
       "      <td>1.211529</td>\n",
       "      <td>1.615421</td>\n",
       "      <td>-0.651393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.732331</td>\n",
       "      <td>-1.326911</td>\n",
       "      <td>2.464234</td>\n",
       "      <td>1.244323</td>\n",
       "      <td>0.654359</td>\n",
       "      <td>0.947023</td>\n",
       "      <td>-1.265609</td>\n",
       "      <td>1.215373</td>\n",
       "      <td>0.765212</td>\n",
       "      <td>-1.431401</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.684100</td>\n",
       "      <td>-1.915556</td>\n",
       "      <td>2.962515</td>\n",
       "      <td>-1.349710</td>\n",
       "      <td>1.875957</td>\n",
       "      <td>1.699892</td>\n",
       "      <td>-1.059458</td>\n",
       "      <td>1.071541</td>\n",
       "      <td>1.476485</td>\n",
       "      <td>-0.699586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.629333</td>\n",
       "      <td>-1.244308</td>\n",
       "      <td>1.316815</td>\n",
       "      <td>1.304162</td>\n",
       "      <td>0.707552</td>\n",
       "      <td>1.429070</td>\n",
       "      <td>-0.895578</td>\n",
       "      <td>-0.007785</td>\n",
       "      <td>0.644606</td>\n",
       "      <td>-1.381937</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.304653</td>\n",
       "      <td>-1.761825</td>\n",
       "      <td>1.265379</td>\n",
       "      <td>-1.320533</td>\n",
       "      <td>-0.609864</td>\n",
       "      <td>0.413826</td>\n",
       "      <td>-0.888624</td>\n",
       "      <td>1.114394</td>\n",
       "      <td>1.519017</td>\n",
       "      <td>-0.798985</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 48 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Actb      Ahcy      Aqp3    Atp12a      Bmp4      Cdx2   Creb312  \\\n",
       "1  0.541050 -1.203007  1.030746  1.064808  0.494782 -0.167143 -1.369092   \n",
       "1  0.680832 -1.355306  2.456375  1.234350  0.645494  1.003868 -1.207595   \n",
       "1  1.056038 -1.280447  2.046133  1.439795  0.828121  0.983404 -1.460032   \n",
       "1  0.732331 -1.326911  2.464234  1.244323  0.654359  0.947023 -1.265609   \n",
       "1  0.629333 -1.244308  1.316815  1.304162  0.707552  1.429070 -0.895578   \n",
       "\n",
       "      Cebpa      Dab2     DppaI  ...      Sox2     Sall4     Sox17     Snail  \\\n",
       "1  1.083061  0.668057 -1.553758  ... -1.351757 -1.793476  0.783185 -1.408063   \n",
       "1  1.208023  0.800388 -1.435306  ... -1.363533 -1.782172  1.532477 -1.361172   \n",
       "1  1.359447  0.530701 -1.340283  ... -1.296802 -1.567402  3.194157 -1.301777   \n",
       "1  1.215373  0.765212 -1.431401  ... -1.684100 -1.915556  2.962515 -1.349710   \n",
       "1 -0.007785  0.644606 -1.381937  ... -1.304653 -1.761825  1.265379 -1.320533   \n",
       "\n",
       "      Sox13   Tcfap2a   Tcfap2c     Tcf23      Utf1    Tspan8  \n",
       "1 -0.031991 -0.351257 -1.078982  0.942981  1.348892 -1.051999  \n",
       "1 -0.501715  1.082362 -0.930112  1.064399  1.469397 -0.996275  \n",
       "1 -0.445219  0.031284 -1.005767  1.211529  1.615421 -0.651393  \n",
       "1  1.875957  1.699892 -1.059458  1.071541  1.476485 -0.699586  \n",
       "1 -0.609864  0.413826 -0.888624  1.114394  1.519017 -0.798985  \n",
       "\n",
       "[5 rows x 48 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# license: Copyright (c) 2014, the Open Data Science Initiative\n",
    "# license: https://www.elsevier.com/legal/elsevier-website-terms-and-conditions\n",
    "URL = \"https://raw.githubusercontent.com/sods/ods/master/datasets/guo_qpcr.csv\"\n",
    "\n",
    "df = pd.read_csv(URL, index_col=0)\n",
    "print(\"Data shape: {}\\n{}\\n\".format(df.shape, \"-\" * 21))\n",
    "print(\"Data labels: {}\\n{}\\n\".format(df.index.unique().tolist(), \"-\" * 86))\n",
    "print(\"Show a small subset of the data:\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to define the output tensor $y$. To predict values for all $48$ genes, we need $48$ Gaussian processes. So the required shape for $y$ is `num_GPs x num_data = 48 x 437`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = torch.tensor(df.values, dtype=torch.get_default_dtype())\n",
    "# we need to transpose data to correct its shape\n",
    "y = data.t()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now comes the most interesting part. We know that the observed data $y$ has latent structure: in particular different datapoints correspond to different cell stages. We would like our GPLVM to learn this structure in an unsupervised manner. In principle, if we do a good job of inference then we should be able to discover this structure---at least if we choose reasonable priors. First, we have to choose the dimension of our latent space $X$. We choose $dim(X)=2$, since we would like our model to disentangle 'capture time' ($1$, $2$, $4$, $8$, $16$, $32$, and $64$) from cell branching types (TE, ICM, PE, EPI). Next, when we set the mean of our prior over $X$, we set the first dimension to be equal to the observed capture time. This will help the GPLVM discover the structure we are interested in and will make it more likely that that structure will be axis-aligned in a way that is easier for us to interpret."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "capture_time = y.new_tensor([int(cell_name.split(\" \")[0]) for cell_name in df.index.values])\n",
    "# we scale the time into the interval [0, 1]\n",
    "time = capture_time.log2() / 6\n",
    "\n",
    "# we setup the mean of our prior over X\n",
    "X_prior_mean = torch.zeros(y.size(1), 2)  # shape: 437 x 2\n",
    "X_prior_mean[:, 0] = time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a sparse version of Gaussian process inference to make training faster. Remember that we also need to define $X$ as a `Parameter` so that we can set a prior and guide (variational distribution) for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = gp.kernels.RBF(input_dim=2, lengthscale=torch.ones(2))\n",
    "\n",
    "# we clone here so that we don't change our prior during the course of training\n",
    "X = Parameter(X_prior_mean.clone())\n",
    "\n",
    "# we will use SparseGPRegression model with num_inducing=32;\n",
    "# initial values for Xu are sampled randomly from X_prior_mean\n",
    "Xu = stats.resample(X_prior_mean.clone(), 32)\n",
    "gplvm = gp.models.SparseGPRegression(X, y, kernel, Xu, noise=torch.tensor(0.01), jitter=1e-5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the [autoguide()](http://docs.pyro.ai/en/dev/contrib.gp.html#pyro.contrib.gp.parameterized.Parameterized.autoguide) method from the [Parameterized](http://docs.pyro.ai/en/dev/contrib.gp.html#module-pyro.contrib.gp.parameterized) class to set an auto Normal guide for $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use `.to_event()` to tell Pyro that the prior distribution for X has no batch_shape\n",
    "gplvm.X = pyro.nn.PyroSample(dist.Normal(X_prior_mean, 0.1).to_event())\n",
    "gplvm.autoguide(\"X\", dist.Normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned in the [Gaussian Processes tutorial](gp.ipynb), we can use the helper function [gp.util.train](http://docs.pyro.ai/en/dev/contrib.gp.html#pyro.contrib.gp.util.train) to train a Pyro GP module. By default, this helper function uses the Adam optimizer with a learning rate of `0.01`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbkklEQVR4nO3df5gcBZ3n8fe3f0zPTCY/gExiSKIJJPwInuFwFhBWjNytJrCS3WfPe2B99JbTYzmB8/YeT+PD7cqdu6su63P4A+VYHp5bPJXzFF1Og6igBy7gZsAQEmJgCIH8gGRCyO/52f29P6p6pqfTM9OT9HR1VX9ezzNPqqurqz9T0/l0dVV1lbk7IiISf6moA4iISG2o0EVEEkKFLiKSECp0EZGEUKGLiCSECl1EJCEiLXQzu9fM9pnZ5iqn/9dm9ryZbTGzb093PhGROLEoj0M3syuAo8B97v72SaZdDnwXuNLd3zSzee6+rx45RUTiINI1dHd/DDhQOs7Mzjazn5jZ02b2uJmdF97174A73f3N8LEqcxGREo24Df1u4BZ3fyfwSeDr4fhzgHPM7B/N7CkzWx1ZQhGRBpSJOkApM+sALgP+j5kVR+fCfzPAcmAVsAh43Mze7u4H6xxTRKQhNVShE3xiOOjuF1a4bxfwlLsPAS+b2TaCgt9Qx3wiIg2roTa5uPthgrL+IIAFVoZ3/xB4bzh+LsEmmO1R5BQRaURRH7b4HeBJ4Fwz22VmHwU+BHzUzJ4FtgBrw8kfBt4ws+eBXwD/2d3fiCK3iEgjivSwRRERqZ2G2uQiIiInL7KdonPnzvUlS5ZE9fQiIrH09NNP73f3zkr3RVboS5Ysobu7O6qnFxGJJTN7Zbz7tMlFRCQhVOgiIgmhQhcRSQgVuohIQqjQRUQSQoUuIpIQKnQRkYSIXaFve/0IX/rpNt44OhB1FBGRhhK7Qn+p9yhffbSHXhW6iMgYsSv01mwQeWCoEHESEZHGErtCz2XSAPQP5SNOIiLSWGJY6OEa+rDW0EVESsWw0IM1dBW6iMhYsSv0kW3ow9rkIiJSKnaFProNXWvoIiKlJi10M7vXzPaZ2eZx7jcz+4qZ9ZjZJjO7qPYxR+W0hi4iUlE1a+j/E1g9wf1rgOXhzw3AN0491vhGdopqDV1EZIxJC93dHwMOTDDJWuA+DzwFzDGzBbUKWK41q52iIiKV1GIb+kJgZ8ntXeG4adGSDiLrOHQRkbFqUehWYZxXnNDsBjPrNrPu3t7ek3qyVMpoSae0hi4iUqYWhb4LWFxyexGwp9KE7n63u3e5e1dnZ8WLVlcll0lpp6iISJlaFPqDwEfCo10uBQ65+2s1mO+4ctm01tBFRMpkJpvAzL4DrALmmtku4LNAFsDd7wLWA1cBPcBx4PrpCluUy6S0DV1EpMykhe7u101yvwM31SxRFXJZbUMXESkXu2+KQvBtUR2HLiIyViwLvTWrnaIiIuViWei5TEpr6CIiZWJa6Gn6tYYuIjJGLAu9NaujXEREysWy0NuyaZ0+V0SkTCwLvTWb1hq6iEgZFbqISELEstBz2ZQ2uYiIlIllobdl0wzmC+QLFU/qKCLSlGJZ6KMXudBmFxGRongWeqZ4kQttdhERKYpnoYdr6H3aMSoiMiKWhd7WEhS6jnQRERkVy0LPZVToIiLlYlnorVltQxcRKRfTQg+PctEauojIiFgXus64KCIyKpaFntNhiyIiJ4hloeuLRSIiJ4ppoWsNXUSkXDwLXYctioicIJaFntMauojICWJZ6MU19J89/3rESUREGkcsCz2VMgCeefVgtEFERBpILAu96MLFc6KOICLSMGJb6CsXzWZ2WzbqGCIiDSO2hd7WkqZvUEe5iIgUxbfQs2mdD11EpERsC729JcPxweGoY4iINIyqCt3MVpvZNjPrMbN1Fe6fbWb/18yeNbMtZnZ97aOO1ZpN6zh0EZESkxa6maWBO4E1wArgOjNbUTbZTcDz7r4SWAV8ycxaapx1jPaWtNbQRURKVLOGfjHQ4+7b3X0QuB9YWzaNAzPNzIAO4AAwrW3b1pLmuHaKioiMqKbQFwI7S27vCseV+hpwPrAHeA74hLufsD3EzG4ws24z6+7t7T3JyIG2bJqB4QKFgp/SfEREkqKaQrcK48pb9P3ARuBM4ELga2Y264QHud/t7l3u3tXZ2TnFqGMVLxStI11ERALVFPouYHHJ7UUEa+Klrgce8EAP8DJwXm0iVtauQhcRGaOaQt8ALDezpeGOzmuBB8umeRX4FwBmNh84F9hey6Dl2sKLXOjLRSIigUkL3d2HgZuBh4GtwHfdfYuZ3WhmN4aTfQ64zMyeAx4BPu3u+6crNMDAcLCJ/qXeo9P5NCIisZGpZiJ3Xw+sLxt3V8nwHuB9tY02sQ07DgBwz+Mvs+rcefV8ahGRhhTbb4recuUyAK5ZeWbESUREGkNsC312W/C9JV0oWkQkENtCn5ELdooeHVChi4hAjAu9eBm6L/7ktxEnERFpDLEt9OJl6M7qnBFxEhGRxhDbQgdYOKdNl6ETEQnFutB3H+zjgWd2Rx1DRKQhxLrQRURklApdRCQhVOgiIgmhQhcRSYhYF/rvLDkt6ggiIg0j1oX+7uXBRTKG87pYtIhIrAu9eJGLYzonuohIvAu9Ixec/ffYwLRej1pEJBZiXegzVOgiIiNiXujFMy6q0EVE4l3oLcEa+nFtQxcRiXmhh5tctIYuIpKQQn95/7GIk4iIRC/Whd6aDeJ/4SFd5EJEJNaF3tmRA+C8t8yMOImISPRiXeiZdIq3zGrlHYtmRx1FRCRysS50CL4tqqNcREQgE3WAU7V9/zH2HOqLOoaISORiv4YO0D+kk3OJiMS+0N9zTmfUEUREGkLsC33FmbPIpg13jzqKiEikYl/os9uyDOWdviHtGBWR5lZVoZvZajPbZmY9ZrZunGlWmdlGM9tiZv+vtjHHN7stC8Cbx4fq9ZQiIg1p0kI3szRwJ7AGWAFcZ2YryqaZA3wduMbdLwA+WPuolb1lVisAew/31+spRUQaUjVr6BcDPe6+3d0HgfuBtWXT/DHwgLu/CuDu+2obc3wzW3VOdBERqK7QFwI7S27vCseVOgc4zcx+aWZPm9lHKs3IzG4ws24z6+7t7T25xGU6wkI/2q9CF5HmVk2hW4Vx5YeUZIB3AlcD7wf+3MzOOeFB7ne7e5e7d3V21uZww+I50Y9oDV1Emlw13xTdBSwuub0I2FNhmv3ufgw4ZmaPASuBF2qScgLa5CIiEqhmDX0DsNzMlppZC3At8GDZNP8AvNvMMmbWDlwCbK1t1MqK50Q/ok0uItLkJl1Dd/dhM7sZeBhIA/e6+xYzuzG8/y5332pmPwE2AQXgHnffPJ3Bi7LpFLlMSmvoItL0qjo5l7uvB9aXjbur7PbtwO21i1a9ma0ZbUMXkaYX+2+KAnTkMjrKRUSaXjIKvTWjC0WLSNNLRKG3ZzMcH1Shi0hzS0ShZzPGUF5nWxSR5hb7KxYBHOobYvPuw1HHEBGJVCLW0FXmIiIJKfQiXeRCRJpZogp9YFjXFhWR5pWoQn/w2fJTzIiINI9EFfqnvrcp6ggiIpFJRKHf928vBuDTq8+LOImISHQSUegXnDkLgPaWdMRJRESik4hCz2WDIu/ZdzTiJCIi0UlEobekg1/jm0+9EnESEZHoJKLQs+ngKnmz27IRJxERiU4ivvpvZlx61ukUdBi6iDSxRKyhQ3BOdF3kQkSaWWIKfWC4wNbXDuvr/yLStBJT6I+/uB+AHW8cjziJiEg0ElPoly87A4CC1tBFpEklptA/fOkSAPqH8tEGERGJSGIKvTUb/Cr9QzrURUSaU2IKfUYuOALzcN9QxElERKKRmEIvfqnoP9z/m4iTiIhEIzGFXtwZeqRfx6KLSHNKTKGfM29m1BFERCKVmEJPpSzqCCIikUpMoZcaGNahiyLSfBJV6LlMeOjioA5dFJHmU1Whm9lqM9tmZj1mtm6C6X7HzPJm9q9qF7F6t11zAQB9+nKRiDShSQvdzNLAncAaYAVwnZmtGGe6LwIP1zpktQaHgzXzHW8ciyqCiEhkqllDvxjocfft7j4I3A+srTDdLcD3gX01zDcl82flAHjz2GBUEUREIlNNoS8Edpbc3hWOG2FmC4E/BO6aaEZmdoOZdZtZd29v71SzTur8BcHFoo8NapOLiDSfagq90vGA5ac0vAP4tLtP2KTufre7d7l7V2dnZ5URqzenvQWAg8e1hi4izaeaS9DtAhaX3F4E7Cmbpgu438wA5gJXmdmwu/+wFiGrNas1+HX2Hu6v59OKiDSEatbQNwDLzWypmbUA1wIPlk7g7kvdfYm7LwG+B3y83mUOwbVFAf7u8Zfr/dQiIpGbdA3d3YfN7GaCo1fSwL3uvsXMbgzvn3C7uYiI1Ec1m1xw9/XA+rJxFYvc3f/k1GOdvPaWNMcH8/QP5WnNpqOMIiJSV4n6piiMHuny7V+/GnESEZH6Slyh3/zeZQCcPa8j4iQiIvWVuELPpIMdoz9/fm/ESURE6itxhV7c5GI6m66INJnEFfoZM4IvF9335CsRJxERqa/EFbpp1VxEmlTiCr3UAZ2kS0SaSCIL/ab3ng3Aa4f6Ik4iIlI/iSz0K8+bD8DGnQejDSIiUkeJLPR5M4Pzot/6g80RJxERqZ9EFvqZc9qijiAiUneJLPR0Ske6iEjzSWShA1yy9PSoI4iI1FViC33ra4cB+OW2yC5xKiJSV4kt9FuuXA7A0YHhiJOIiNRHYgv9inOCa5Z+7dGeiJOIiNRHYgt9Vltw7Y7fvn4k4iQiIvWR2EJfMFuHLopIc0lsoZc6pu3oItIEEl3oq84NtqN/4Ku/ijiJiMj0S3Sh/5erzwdg+/5jEScREZl+iS70ZfNmArBy0eyIk4iITL9EF3rRs7sO0TeYjzqGiMi0aopCB7jjkReijiAiMq0SX+inh9cY3X9EVy8SkWRLfKF/90/fBcD3n9kVcRIRkemV+EJfNq8j6ggiInWR+EIv9eRLb0QdQURk2jRVof/V+uejjiAiMm2qKnQzW21m28ysx8zWVbj/Q2a2Kfx5wsxW1j7qydv631YDsHn34YiTiIhMn0kL3czSwJ3AGmAFcJ2ZrSib7GXgPe7+DuBzwN21Dnoq2lrSI8OFgkeYRERk+lSzhn4x0OPu2919ELgfWFs6gbs/4e5vhjefAhbVNmbtfOaB56KOICIyLaop9IXAzpLbu8Jx4/ko8NCphJpO/7t75+QTiYjEUDWFbhXGVdxuYWbvJSj0T49z/w1m1m1m3b29vdWnrIF7PtJV1+cTEam3agp9F7C45PYiYE/5RGb2DuAeYK27Vzw+0N3vdvcud+/q7Ow8mbwn7V+umD8yPDCs87qISPJUU+gbgOVmttTMWoBrgQdLJzCztwIPAB9294Y/acqnvrcp6ggiIjU3aaG7+zBwM/AwsBX4rrtvMbMbzezGcLK/AM4Avm5mG82se9oS18A/bDzhA4aISOxlqpnI3dcD68vG3VUy/DHgY7WNVnsPfeLdrPny41HHEBGZFk31TdHzF8waGb7xm09HmEREpPaaqtABvvGhiwD4yZbX9SUjEUmUpiv0Nf9swcjwdza8GmESEZHaarpCB7jlymUA3PqDzREnERGpnaYs9Bvfc/bI8Ld+/UqESUREaqcpC31GLsPy8MIXt/5gM+7ali4i8deUhQ7w0z+7YmR46WfWTzCliEg8NG2hmxmf/cDoWYD/4M5/jDCNiMipa9pCB7j+8qUjwxt3HuS//6zhz1ogIjKupi50gBf+cs3I8JcfeZGHnnstwjQiIiev6Qu9JZPi5c9fNXL733/rGW5/+LcRJhIROTlNX+gQbE8vXncU4M5fvMSSdT/mcP9QhKlERKZGhR5qa0nzzJ//3phxf3qfzvciIvGhQi9x+oyWMZtfntz+Bh/46q8iTCQiUj0VehkzY8cXrh65/dzuQyxZ92P2Hu6PMJWIyORU6OPY8YWrueitc0ZuX/LXj7Bk3Y85NjAcXSgRkQmo0CfwwMcvZ9Nt7xsz7oLPPsy6729iOF+IKJWISGUW1XlMurq6vLu7oa9UN8bew/1c8tePnDD+bz+4kj+6aCFmFkEqEWk2Zva0u3dVvE+FPjX7jw7Q9Zc/r3jfi3+1hmxaH3pEZPqo0GvM3Xn9cD/v+vyjJ9zXkk7xwMcv4/wFs0intNYuIrWlQp9GhYLzlUdf5I6fvzjuNP/jw+/k986fT0oFLyKnSIVeJ+7Oxp0H+cOvPzHuNBcunsNt11zAykWztd1dRKZMhR6BQsF5cd9R3n/HY5NOe+HiOXzyfeeycvFsZrZm65BOROJKhd4AhvIFdh44zg9/s5uvPNoz6fS/u2wuv+rZz2Vnn8F7zunk31y2hNZsug5JRaSRqdAbUL7gHB8cZueBPj7/0FYef3H/lB5/WnuW8xfM4sw5bSydO4M/umgR82fltBlHJOFU6DHj7uw+2Meq23/JcGHqf5+WTIrB4bFffDq7cwaZVIo/vuStzJ+VY9m8mcxqy3B6ewspM8zQm4FIDKjQE+aNowPsOdjPKweO8bkfPc/ewwOc95aZvHaon0N9p3bK37kdOVIGs9uyvO2MGSw5o50Nr7zJufM7WDJ3BnM7chhwpH+Yt53Rzqy2LEcHhpk3M0dHLkNbSxoccpk0uWyKXCY4Ll9vFiK1oUJvUoWC0z+cZ8/BfrbsOcSO/cd549gA214/wuLT2zl9RgtPvLSfzbsPjzzmtPYsbx4P3hRmtWboHyowWIPTHHTkMrS3pMmmU8zIpUmZkU4Z2XSK3Qf76Mhl2He4n4vedhpt2TQFh5mtmeCTA0Z7SxozcIfWbIpUymjNpGnJpMgXnGw6Rb5QIJdJk00brdlg+kxqdP5ndLSQSaXIpA0jeJMJ/oVU8EQj49NmtLWkyBegLZtmuFAgk0qFn2SCTMXhVMl8KI4Hsulg+nTKRqah9LGMvtEVH29hiOL7X8qMVPh7p8xIpSBtFg4bQ+HfJh1+ynJnzOGx+YLr+xAJM1GhZ+odRuonlTLaWzIsm9fBsnkdJz2fI/1DpMzoG8rTN5gnX3CODgwzmC9wuG+IjTsPsmB2Kykz9h7up70lw94j/cxqzdJ7ZIDuVw5w1twO2lvS9A3lGRgqMFxwBvMFBofzzGnLcnwwz7HBPDsPHCeXSVNwHzkRWsFhYDgf7nfIn9RmqGaUsmDZlX44SpvhBJv1iksxm0qBBfcVwvFG8EaUL3jwphS+iTjBm0QmFb4jOSPTmzEyT3dK3sDCPOEbuLtTvh5ZzOjh/Nw9eJMtztdLpyt9Q6z8Bkk4j+Lv7x78PoVwvif7RhcsqtFcY36Hkd8lXI4e/A2K40qXx3UXv5WPvfusKT//ZFToMqnioZQzcpVfLqvOnVfPOCMK4ZuCOwwVCmTCAjo+mGcoHF9882nNBmvbA8P5oJgcHB9TIMG/AD7yySSbStE3lCedgnyBkf+ojHksI0VY/MTrTvDJxiHvPvI4D+8cfa6xz106Dwjmmy+MlnPBnULBg3mGn0xGH+vhY4IRhXBtvXR+xYJOlbTRUN5HlkWxdPDRtXsPl3XeHcNIp2C44GNKtnw4KP+xxV1wDz9RBJ86Sku8VDFbwYOsBR8ty9K/kTsV/o6j8wl+z/DNJZxPcb7F18B4WwLH23BRKPvblCyu0b8nkLbR4fLXGA6dM3OVn+AUVVXoZrYa+DKQBu5x9y+U3W/h/VcBx4E/cfdnapxVZIxUymhNBYdytjF6SKeO5ZdmNemZpMwsDdwJrAFWANeZ2YqyydYAy8OfG4Bv1DiniIhMoppTA14M9Lj7dncfBO4H1pZNsxa4zwNPAXPMbEGNs4qIyASqKfSFwM6S27vCcVOdBjO7wcy6zay7t7d3qllFRGQC1RR6pd0G5bsMqpkGd7/b3bvcvauzs7OafCIiUqVqCn0XsLjk9iJgz0lMIyIi06iaQt8ALDezpWbWAlwLPFg2zYPARyxwKXDI3V+rcVYREZnApIctuvuwmd0MPExw2OK97r7FzG4M778LWE9wyGIPwWGL109fZBERqaSq49DdfT1BaZeOu6tk2IGbahtNRESmIrJzuZhZL/DKST58LjC1883WR6PmgsbNplxTo1xTk8Rcb3P3ikeVRFbop8LMusc7OU2UGjUXNG425Zoa5ZqaZstVzU5RERGJARW6iEhCxLXQ7446wDgaNRc0bjblmhrlmpqmyhXLbegiInKiuK6hi4hIGRW6iEhCxK7QzWy1mW0zsx4zWxfB8+8ws+fMbKOZdYfjTjezn5nZi+G/p5VM/5kw6zYze38Nc9xrZvvMbHPJuCnnMLN3hr9Pj5l9xU7xas7j5LrNzHaHy2yjmV0VQa7FZvYLM9tqZlvM7BPh+EiX2QS5Il1mZtZqZv9kZs+Guf5rOD7q5TVershfY+E802b2GzP7UXi7vssruL5fPH4ITj3wEnAW0AI8C6yoc4YdwNyycX8DrAuH1wFfDIdXhBlzwNIwe7pGOa4ALgI2n0oO4J+AdxGcMfMhYM005LoN+GSFaeuZawFwUTg8E3ghfP5Il9kEuSJdZuE8OsLhLPBr4NIGWF7j5Yr8NRbO8z8B3wZ+FMX/ybitoVdzsY0orAX+Phz+e+APSsbf7+4D7v4ywbluLq7FE7r7Y8CBU8lhwUVIZrn7kx68ku4reUwtc42nnrle8/CyiO5+BNhKcM7+SJfZBLnGU69c7u5Hw5vZ8MeJfnmNl2s8dXuNmdki4GrgnrLnr9vyiluhV3UhjWnmwE/N7GkzuyEcN9/Ds0uG/xavmlzvvFPNsTAcrke+m81skwWbZIofOyPJZWZLgH9OsHbXMMusLBdEvMzCzQcbgX3Az9y9IZbXOLkg+tfYHcCngELJuLour7gVelUX0phml7v7RQTXUb3JzK6YYNpGyAvj56hXvm8AZwMXAq8BX4oql5l1AN8H/qO7H55o0npmq5Ar8mXm7nl3v5Dg+gYXm9nbJ5g86lyRLi8z+31gn7s/Xe1DpiNX3Ao98gtpuPue8N99wA8INqHsDT8qEf67L5y83nmnmmNXODyt+dx9b/ifsAD8HaObneqay8yyBKX5LXd/IBwd+TKrlKtRllmY5SDwS2A1DbC8KuVqgOV1OXCNme0g2BR8pZn9L+q9vE51J0A9fwhO97udYCdCcafoBXV8/hnAzJLhJwhe5LczdsfH34TDFzB2x8d2arRTNJz/EsbufJxyDoILmFzK6A6Yq6Yh14KS4T8j2HZY11zhfO4D7igbH+kymyBXpMsM6ATmhMNtwOPA7zfA8hovV+SvsZLnX8XoTtG6Lq+aFEs9fwgupPECwV7hW+v83GeFf4RngS3F5wfOAB4BXgz/Pb3kMbeGWbdRg73oJfP9DsFHyyGCd/WPnkwOoAvYHN73NcJvD9c41zeB54BNBFe3WhBBrt8l+Oi6CdgY/lwV9TKbIFekywx4B/Cb8Pk3A39xsq/1OuWK/DVWMt9VjBZ6XZeXvvovIpIQcduGLiIi41Chi4gkhApdRCQhVOgiIgmhQhcRSQgVuohIQqjQRUQS4v8D8RyHA+CMy+oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# note that training is expected to take a minute or so\n",
    "losses = gp.util.train(gplvm, num_steps=4000)\n",
    "\n",
    "# let's plot the loss curve after 4000 steps of training\n",
    "plt.plot(losses)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After inference, the mean and standard deviation of the approximated posterior $q(X) \\sim p(X | y)$ will be stored in the parameters `X_loc` and `X_scale`. To get a sample from $q(X)$, we need to set the `mode` of `gplvm` to `\"guide\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "gplvm.mode = \"guide\"\n",
    "X = gplvm.X  # draw a sample from the guide of the variable X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s see what we got by applying GPLVM to our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAGJCAYAAABrSFFcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB3b0lEQVR4nO29eZgV1bW//66ekaFbBKQBvQhBwAFFUSQGI7ZiHAgGvUZj4hS//LxXbxyuRjSJQTOIiQYxMfFyjYo3RmLUiK3GqUElKCqKEhQQQaJAIyp2C03PvX9/1KnmnNNVdarOPKz3ec7TfXbt2rWrh7Nq7b3WZ4kxBkVRFEVRCoOiTE9AURRFUZT0oYZfURRFUQoINfyKoiiKUkCo4VcURVGUAkINv6IoiqIUEGr4FUVRFKWAUMOvpAURmSQiC0Vks4i0iciXIvKGiPxMRKqj+pqwV4eIbBSRe0VkWFifC0PHvxJ17rzQORFjhh0XEfmXiLwYen982LWmOvQfLiJdoeOXJOWHkSAicriIPCoiH4lIq4jUi8gSEflBWJ/hoTlfmMJ5pOQaItJbRK4XkbdEZKeItIjIOhH5XfTv2+d4RkRmh72fLSJJy2O2xwt7NYjI6yLyHYe+vu5NRF6MGnOniCwTkW8mONdNInJ/HOedISJXJ3JtJXtQw6+kHBH5b2AZMBD4MXAicA7wLDATuNfhtPuBScDxwG+A6UCdiPSKcbkFQDHQ40M3xPHA/qF+4ewEvufQ/3xgV4xrpg0ROQpYDgwAfgicDFwLrAO+Fda1Huvn91S655gIoQe217Hu7SngLOAU4E6s+/lr5mYXk69hzfE7wBbgQRG52D4Yx72tCrVPAr4P9AYeE5GJqb0NR84A1PDnCSWZnoCS34jIFODXwDxjzFVRh58WkVuAf3c4dYsxZnno+3+IyE6sh4FTgMfcrmeMeUtEVmMZ7NsdupwP7AYeiWp/DDhLRHobY5rC2r8HPApc6HbNNPNfQAMw1RjTGtb+JxHpfpAPHVtO7vF/QDVwtDFmfVj7EhH5PdYDYLbymjGmA0BEngPWAFey58E26L3tDPsfWC4irwAfARcBr6XmFpRCQD1+JdVcB3wW+toDY0yTMeZ+H+O8EfrqZ6l3ATBORMaFN4rIXsCZwN+MMTujznkMMMCMsP5fBUZifWD7QkSOFpEXRGSXiDSJSJ2IHB3V5/7Qlsd4EVkqIrtFZL2IXOrjEv2BL6KMPgDGmK6wa/RYhg9yXRE5UURWhpaiPxCRS0Lnb/LxM/h66L53hn4Gz4rIIT7OOxqoAX4ZZRjt+zPGmMejzvl/IvJOaJ6ficgfRaR/rGv5QUT2EpHfi8jnod/nEyLyNT/bG6EHgJWE/l7juTeHPpuBT7FWrPzM/4rQ0n6LiKwQkckOfQaKyP+IyPuhv4ePReTPIjI0rM/9wAXA0LCth02hYxUiMldEVod+RttEpFZExviZo5IZ1PArKUNESoCvA88bY9oSHO6A0NcGH33/BHTSc+n+DKAvPZf5wVoFeDTqnPOxtig2+plg6EHjJWBvrBWC84F+wEsiclhU937An0NznY71YPOH0AqJF68DY0Tk7tBDRtBVu5jXFZGDsJaid2FtydwAXAGcEGtwETkNqAud+12sZe++wFIR2S/G6SeGvj7h50ZEZA7we+AF4JtYWx7fAP4uIsV+xojB/wCXYG01zcDaTvlzgPMPYM/fa6B7c0JE+gL7ABt89P0+cAewBOvv/n7gIay/zXD6Ay3A9Vg/u2uBUcAyEakI9fkZ8DTWQ4e99WBvK5Vj/X5/DpwG/AdQgbVCMTjwTSrpwRijL32l5AXsi+VF3+JwrCT8FXXMAL8IHasAjsFaNm0ChoT6XBjq9xWXaz8NbAWKw9qeATYDRWFtx4fGORHLsHUCQ7E+0HYA/w8YHupzSYz7fQTrg74qrK1faJzHwtruD403JaytHGtlZH6Ma/QC/hY632A9sDwXmmf4vdpzvjDodbGM26fAXmFt1VgGYlOMa3wA1EXNuV/oGnfEuLc/hMYr9/G3NTz0u7oxqv3Y0BhnRP09zQ57P9v66PMcf3Ro/Fkuc7wwerzQz7IEGAT8NNR2R9B7C/V/EfgHe/5HDsCKAdgOjIhxbhHwMfBMVPu3Q3O43+PcYmC/UL9vRf3tbPYx72JgL6yYmav83Ku+0v9Sj19JJeLYaHkC7eEvB8/1htCxZuDV0PenGmO2+rz2AixjVRO6ZjWWcf+TCVsSj2IJ1oPBd4BpWEb2YZ/XAzgOeNIY02A3GGO+xPLyvh7Vd7cxZklYv1ZgPTGWcY0xzcaYbwEHY3lnfwcmAPOxYiYcf+YBr3sM8LQxZndYv3rgFa+BRWQU1tbIgyJSYr+wHk5exfr5ICJF4cfj9M5PwjJw0dd6DfjSvlYCTAyNH/37X+hxTgvW3+knWH+/dwCzEpjDsez5H9mI9Td5pjEm1grUsNAreu6PAh3RnUXkP0LbJbtCxz8KHRrtZ5IicraIvCYiDaHzm4A+fs9X0o8afiWVfIb1YRhtzD4Djgq9/tfl3HtDx8cDA4wx44wxLwW49iIs7/v80PvzsLyRB9xOMJbL8iDWcv8FwBPGmMYA1+yPFU0fzTZ6LrF+4dCvFWuFIybGmPeMMbcZY84EhmAt3U/FWm71ws91q7E8y2g+iTH2oNDXPxL1YAecjrVMDXBj1LG6UPvHoa//FuM64df6wOFa/cKuFS92Omj0PXv9DI7B+pv9CtDHGHOVMaYldCzIvdm8ExrvGKyo/p3AX0VkYIzzHOdurLiDz8PbROS/2LNdMgM4OnQ98PG3KCLTgL9grch9B+uB6SisFSNff8tK+tGofiVlGGM6RORl4CQRKTOhff7QB9AKABE53eX0emPMigSu3SIiDwPfFZE+WMb8DWPMezFOfQBrv/NgrH3jIOwAnPY1B4eOpYTQvf4aa0/9IODJBIesZ49hDWffGOfZRuV6LEMSjR3nMZ/IOdqBli9gbfFMwzkjw+laU3F+mPncoS0I9gPcvkTGeHj9DN4M/W07EeTebHaF/Q+8JiIfAouxthYu8zgvfO7dhFZEoh+IzsHamvnvsH4H4J9zgA+MMReGnV+K9RCsZCnq8Sup5ldYOee3ZuDaC7D2G28CxuHh7dsYY9YCd2Ht1z8b8HovAaeFgrCA7oCsaaFjCSNhIkZR2FHUTisOQVkOnCpWFoR93WqspWcv1gGbgIONMSscXqsAjDFbo9rXhdpfx/L+bxAXoR4RsVPenge6gP1drvVhAvcP1pZBF3B2VPs58QwW8N7cxliCFd9xicffAVjbVR/Tc+5n0tPZ2wtrlSScixzGbMXa+opmL3puH3wPa3VNyVLU41dSijGmTkRmAXNCUe8PAB9iLQMeiPVB2oQVTBQP3xCRbVFtjcaY540xr4jIeuAqrA+3h3zO+fI45/IzrCXtOhG5FeuersP6cLw5zjGjuVtE9sVKMVyN9QF7FJYozAYsw5AoP8cSl3lWRG7DClr7CdbSsVt8BMYYIyKXAYtEpAxrj/kzLM/zq8BHxpjfxLj297C84zdE5LdYAW5tWA82FwOlwCJjzIbQz/h3IjIa68GqBSsw7STgnvBYhqAYY9aJyJ+Bm8XSR3gjNO6p8Y6Jz3uLMcaNWNkY12FpOjjNvUtEbgLuEZH7sOISvoK1EvNlVPdngOtE5AasjJETsH730bwH9BeR/8BarWsxxvwzdP4ZIjIXaxXnSOAH+Mu+UTJFpqML9VUYLyxv8WEsRbM2rA+gN7C88eqovgb4eYzxLmRPZHv0a3VYvx+H2v7mMs7xoeMnelxrOD6i+kN9J2J9uO/CeqCpwxJsCe9zPw4R0liR3C/GGP9krJWMdVhL5K1YBv93wL4Oc74wnutiGbm3Q+NvBP4/rIeKlV7XCLVPwjICXxDKBMAyPpN8/q30wQqOWxn6GbaG7nceURHtWMZ0eajfLqy95t8Bw6L+nmaHvZ9NjKj+UL+9sKLxd4TGfoI9WQMXRo9HVHZKIvcW+p38w2WMP2MFvVbHuNYVwL9Cv4MVWMqCmwiL6sfy4v+AtSe/M/R7O8DhZ9Yb68H5i9CxTaH2IqwHxa1YQZwvYcXlRFxHX9n1ktAvT1EUxZVQnMQHwFPGmO9nej6ZQkSGY61YXWT8CU8pStahS/2KovQgtBT9CpYnNwTLe9wbyzNVFCWHUcOvKIoTFVgBmftibc28jrUdsiqjs1IUJWF0qV9RFEVRCghN51MURVGUAkINv6IoiqIUEAWxxz9gwAAzfPjwTE9DURRFUdLCm2+++ZkxxlHeuSAM//Dhw1mxIm71V0VRFEXJKUTkX27HdKlfURRFUQoINfyKoiiKUkCo4VcURVGUAqIg9vidaG9vZ/PmzbS0tMTuXCBUVFQwbNgwSktLMz0VRVEUJUUUrOHfvHkzffv2Zfjw4YhIpqeTcYwxfP7552zevJkDDghSjltRFEXJJQp2qb+lpYV99tlHjX4IEWGfffbRFRBFUZQ8p2ANP6BGPwr9eSiKouQ/BW34M83FF1/MoEGDOOSQQzI9FUVRFKVAUMOfQS688EKeeeaZTE9DURRFKSDU8Pvk0W07mPDKu1QveZsJr7zLo9t2JDzmcccdR//+/ZMwO0VRMs6qh2HuITC7yvq66uFMz0hRHCnYqP4gPLptB9es+5jmLquE8ebWdq5Z9zEAZw5Ww60oBc+qh6H2B9DebL1v/Nh6DzDu7MzNS1EcUI/fB7dsrO82+jbNXYZbNtZnaEaKomQVdTfvMfo27c1Wu6JkGWr4fbCltT1Qu6IoBUbj5mDtipJB1PD7YGi5s5KdW7uiKAVG5bBg7YqSQdTw++D6EdX0KorMce9VJFw/ojqhcc8991wmTZrEunXrGDZsGH/84x8TGk9RlAxRcyOU9opsK+1ltStKlqHBfT6wA/hu2VjPltZ2hpaXcv2I6oQD+x566KFkTE9RlExjB/DV3Wwt71cOs4y+BvYpWYgafp+cObi/RvAriuLOuLPV0Cs5gS71K4qiKEoBoYZfURRFUQoINfyKoiiKUkCo4VcURVGUAkINv6IoiqIUEGr4M8THH3/MlClTGDt2LAcffDDz5s3L9JQURVGUAkDT+TJESUkJt99+O0cccQQ7d+7kyCOP5KSTTuKggw7K9NQURVGUPEYNv0/qty1i44bbaGmtp6K8mhEjr6F68PS4x6uurqa62lL+69u3L2PHjmXLli1q+BVFUZSUoobfB/XbFrF27Y/o6rKqb7W0bmXt2h8BJGT8bTZt2sTKlSuZOHFiwmMpiqIoihe6x++DjRtu6zb6Nl1dzWzccFvCY+/atYszzzyTO+64g379+iU8nqIoiqJ4oYbfBy2t9YHa/dLe3s6ZZ57Jeeedx4wZMxIaS1EURVH8kFWGX0S+ISLrROQDEZnlcHy6iKwSkbdFZIWIfC0d86ood67C59buB2MM3//+9xk7dixXX3113OMoiqIoShCyxvCLSDFwF3AKcBBwrohER7rVAYcZYw4HLgbuScfcRoy8hqKiyJKbRUW9GDHymrjHXLZsGf/3f//H4sWLOfzwwzn88MN5+umnE52qoiiKoniSTcF9RwMfGGM2AojIQmA68J7dwRizK6x/b8CkY2J2AF8yo/q/9rWvYUxapq8oiqIo3WST4R8KfBz2fjPQI8xdRL4F3AIMAk5Lz9Qs45+MCH5FURRFySRZs9QPiENbD5fYGPM3Y8wY4AzgZ66DicwMxQGs+PTTT5M3S0VRFEXJYbLJ8G8G9gt7PwzY6tbZGPMyMFJEBrgcn2+MmWCMmTBw4MDkzlRRFEVRcpRsMvxvAKNE5AARKQPOAZ4I7yAiXxERCX1/BFAGfJ72mSqKomQrqx6GuYfA7Crr66qHE+un5B1Zs8dvjOkQkcuBZ4Fi4F5jzLsicmno+N3AmcD5ItIONAPfNhohpyiKYrHqYaj9AbSHBMcaP7beA4w7O3g/JS/JGsMPYIx5Gng6qu3usO9vBW5N97wURVFygrqb9xhzm/Zmqz3coPvtp+Ql2bTUX5B0dnYyfvx4Tj/99ExPRVGUXKdxs792v/2UvCSrPP5CZN68eYwdO5Yvv/wy01NRFCXXqRxmLdtH02tvax+/cbPVp9fe0LzD+Xwl71GP3yerVq1i7ty5zJ49m7lz57Jq1aqEx9y8eTNPPfUUl1xySRJmqChKwVNzI5RGqoxSXAatO0MPBMb62rYLikoj+5X2ss5X8h41/D5YtWoVtbW1NDY2AtDY2EhtbW3Cxv/KK6/kV7/6FUVF+mtQFCUJjDsbpt0JlfsBYn0t6wNd7ZH9OtugvG9kv2l36v5+gaBL/T6oq6ujvT3yH6e9vZ26ujrGjRsX15hPPvkkgwYN4sgjj+TFF19MwiwVRVGwjHe4AZ9d5dyv+Qu47sO0TEnJLtTV9IHt6ftt98OyZct44oknGD58OOeccw6LFy/mu9/9btzjKYqiOOK2b6/7+QWLGn4fVFZWBmr3wy233MLmzZvZtGkTCxcu5IQTTuBPf/pT3OMpiqI44rTvr/v5BY0afh/U1NRQWhoZCFNaWkpNTU2GZqQoiuITp31/3c8vaHSP3wf2Pn5dXR2NjY1UVlZSU1MT9/5+NMcffzzHH398UsZSFEXpQfS+fzirHraEe+xUv5ob9aEgz1HD75Nx48YlzdAriqJkBSrdW5DoUr+iKEqh4iXdq+QtavgVRVEKFZXuLUjU8CuKohQqmupXkKjhVxRFKVQ01a8gUcOvKIpSqGiqX0GiUf0ZZO7cudxzzz2ICIceeij33XcfFRUVmZ6WoiiFhFeqn5KXqMefIbZs2cKdd97JihUrWL16NZ2dnSxcuDDT01IURVHyHPX4fdK0cjtfPruJzoZWiqvK6XfycHqPH5TQmB0dHTQ3N1NaWsru3bsZMmRIkmarKKllzdIlLF34ADs//4y++wxg8jnnM3bylExPS1EUH6jH74OmldtpeGw9nQ2tAHQ2tNLw2HqaVm6Pe8yhQ4dyzTXXsP/++1NdXU1lZSVTp05N1pQVJWWsWbqE5+b/jp2ffQrGsPOzT3lu/u9Ys3RJpqem5AOrHoa5h1hVBeceYr1Xkooafh98+ewmTHtXRJtp7+LLZzfFPeYXX3zBokWL+PDDD9m6dStNTU1apEfJCZYufICOttaIto62VpYufCBDM1LyBltJsPFjwOxRElTjn1TU8PvA9vT9tvvhhRde4IADDmDgwIGUlpYyY8YMXnnllbjHU5R0sfPzzwK1K4pvVEkwLajh90FxVXmgdj/sv//+LF++nN27d2OMoa6ujrFjx8Y9nqKki777DAjUrii+USXBtKCG3wf9Th6OlEb+qKS0iH4nD497zIkTJ3LWWWdxxBFHcOihh9LV1cXMmTMTnKmipJ7J55xPSVnkQ29JWTmTzzk/QzNS8gZVEkwLGtXvAzt6P9lR/TfddBM33XRTMqaoKGnDjt7XqH4l6dTcGFktEFRJMAWo4fdJ7/GDEjb0ipIvjJ08RQ19vrHqYWsvvXGz5WHX3Jh+YR/7epmeR56jhl9RFKXQsaPpbU/bjqaHzBh/NfQpRff4FUVRCh2Npi8o1PAriqIUOhpNX1Co4VcURSl0NJq+oFDDryiKUujU3GhFz4ej0fR5ixr+DHLxxRczaNAgDjnkkIj23/72t4wePZqDDz6YH/7whxmanaK4s2bpEuZfdhG3nzON+ZddpDr9uc64s2HanVC5HyDW12l3RgbZqYZ+3pBVUf0i8g1gHlAM3GOMmRN1/DzgutDbXcB/GGPeSe8sk8eFF17I5Zdfzvnn7xE+WbJkCYsWLWLVqlWUl5ezfXv8hYAUJRXYRXpsvX67SA+gKX65jFc0vVfUP2j6XY6RNYZfRIqBu4CTgM3AGyLyhDHmvbBuHwJfN8Z8ISKnAPOBiemYXyrKkB533HFs2rQpou0Pf/gDs2bNorzcUkYbNEi1A5TswqtIjxr+PMUt6v/v10FHc3akASq+yaal/qOBD4wxG40xbcBCYHp4B2PMK8aYL0JvlwNpiTxJZxnS999/n6VLlzJx4kS+/vWv88YbbyT9GoqSCFqkpwBxi+5v3qFpgDlINhn+ocDHYe83h9rc+D7w95TOKEQ6y5B2dHTwxRdfsHz5cn79619z9tlnY4xJ+nUUJV60SE8BEjS6X9MAs5psMvzi0OZo8URkCpbhv87peKjPTBFZISIrPv3004Qmlk4PZ9iwYcyYMQMR4eijj6aoqIjPPlNPSsketEhPAeIW9d+rv3N/TQPMarLJ8G8G9gt7PwzYGt1JRMYB9wDTjTGfuw1mjJlvjJlgjJkwcODAhCaWTg/njDPOYPHixYC17N/W1saAAepJKdnD2MlTOPjrNUiR9fEhRUUc/PWamPv7mgmQIzhF77tF/Z9yq6YB5iBZE9wHvAGMEpEDgC3AOcB3wjuIyP7AY8D3jDHvp2tik885PyKKGZLj4Zx77rm8+OKLfPbZZwwbNoybbrqJiy++mIsvvphDDjmEsrIyFixYgIjTYoiipA6vYNY1S5fw7kt1mK4uAExXF+++VMfQ0WMZO3mK47mAZgJkG05FecBbs98tYE+j+nMKyab9YxE5FbgDK53vXmPML0TkUgBjzN0icg9wJvCv0CkdxpgJscadMGGCWbFiRUTbmjVrGDt2rO+5pSKqPxsJ+nNR8o/odD2wHnSnzrycsZOnMP+yi6xA1yj6Dhjo+pBcUl5Gy86djufMvOu+1NyI4k50eh5YnnpJLytgL5rK/eCq1cmfgz4wpAwRedPNPmaTx48x5mng6ai2u8O+vwS4JN3zAi1DmsusWrWKuro6GhsbqayspKamhnHjxmV6WllLrHQ9r5gXt3Oj28LPUTKAW3pedJtNsoP1sqkaYAGSVYZfUZLNqlWrqK2tpb29HYDGxkZqa2sB1Pi7ECuYte8+Axw9fhFxbPciOk6mUFbWMk5QQ57sYD2vaoBq+FNONgX3KUrSqaur6zb6Nu3t7dTV1WVoRtlPrGBWp6h+oHvP34nyPn0dz9n52afdgX7p1MsoeNwMea/+6QnW02qAGUUNv5LXNDY2BmpXYqfrjZ08hakzL6fvgIEgYr08KCkrp+bCmXvOiWLnZ5/y9O9u5+nf3Z42vYyCxY7Yb/zY+fjB34qt2e81rl8df60GmFF0qV/JayorKx2NfGVlZQZmkxvYS+teS+52zMuapUt4+ne3u45lB/zZ53oFB7qhcQAeBAmQcwroi2b9c3D6b4Itt8ezX19zo3NwoaYBpgU1/EpeU1NTE7HHD1BaWkpNTU0GZ5X9+A1m9fLG3SL2gxpyVQR0IajBddpXj6bxY8tjDxJhH89+vd2uUf0ZQQ1/hmhpaeG4446jtbWVjo4OzjrrLG666SYArr32WmpraykrK2PkyJHcd999VFVVdZ/7z3/+k+9973sAfPTRR1RWVlJZWcmAAQO45557GDt2LKNHj+7uf/XVV0dUACwk7AA+jepPDtHBd16eu5POxZqlSxAR3zLUqgjoQVCD63f/3I/HHr7S4CywGvt6XroASkpRw58hysvLWbx4MX369KG9vZ2vfe1rnHLKKRxzzDGcdNJJ3HLLLZSUlHDddddxyy23cOutt3afe+ihh/L2228DVmnf008/nbPOOguATZs2MXLkyO7jimX8k23oCzFF0Kkcrxvlffr2WDGwz/cKAgwneptAiSJogFzlMPe9/Wi8HiD8bBnY13NC8/czjhp+n7z/2jZeXbSBXTta6dO/nEnTR3LgxMFxjyci9OnTB7CizNvb27sV+qZOndrd75hjjuGRRx5JbPJKUinUFEGnHH0n7GC+RM63xYIUF1Y9DFIEprPnsV57O58zaiqsuJdID10I7LH72TJw26/X/P2sQKP6ffD+a9tY8uBadu2wPrR27WhlyYNref+1bQmN29nZyeGHH86gQYM46aSTmDhxYo8+9957L6ecckqgcTds2MDhhx/e/Vq6dGlC81QiKdQUQa+9eTvCv++Aga5GO9HzlRC28XQy+gAtjXDrAZER9k9e7Wz0J1wciuB3wM1j91zCj5EJ4LU9oaQN9fh98OqiDXS0RS5PdrR18eqiDQl5/cXFxbz99ts0NDTwrW99i9WrV3PIIYd0H//FL35BSUkJ5513XqBxdak/tRRqiqDbnr5f2d1Ez1dCxPK4Tece2d3Gj+Hx/4SudqeOVhR/0Ah7ty0DP7K+rtsToaDCQlz+z8DWh3r8PrA9fb/tQamqquL444/nmWee6W5bsGABTz75JA8++KAW6cky3FIB8z1FMNFyvFrO11o9XHDDMu66dDELblgW36phUJEbR6MfNpZb5b1o4xOhARD1meQ3Fc81T19C45o9y/+xtACCElRrIB3Yqzepvvco1PD7oE//nopjXu1++PTTT2loaACgubmZF154gTFjxgDwzDPPcOutt/LEE0+w1157xX0NJTXU1NRQWloa0VYIKYLRwj1Bl+Xt88v79O1uKykvS9V0s46kbRkmU+TGHmvc2Za3PrvB+upk9LsNFFhbBiHj71fkB6yHg2hlQKc4g2Qv/2fIwMYkQ1sfutTvg0nTR7LkwbURy/0lZUVMmj4y7jHr6+u54IIL6OzspKuri7PPPpvTTz8dgMsvv5zW1lZOOukkwArwu/vuu72Gi8De47e5+OKL+cEPfhD3XPORRKLyCzVFMFk6+p1tbd3ft+zcWTDleZO2Zei0NB8X4l8wx3F7wQSv2ueUv++WaZBM+d5srQ2QIeliNfw+sP8pkxnVP27cOFauXOl47IMPPvA9zv333x/xfvjw4TQ3J/qBkN8kIyo/FSmC2YxTKl88Bnvxgvmelf/ymaRtGUYbTzuK36mcrk1xGXS2hTWEAvv8Gr1kGqjo/H1XCWETXEzIjWytDeAaL5Fa6WI1/D45cOLghAy9kj14ReUXkjEPglu53b//fi7gz/ivWbqElp07HY8Vgixvn/7ljkY+ri1DJ/Gb2VW4puZNvyuxALJ4DZSfwDWvFYzGj+GxmfDRcktKOF4yZGBjkiHpYt3jVwqOQo3KTwQ3w2y6unxX0POU9y0AWd5J00dSUhb5kZvolmEEroVv9ou9hx8Lp735WAbK7756RHChE8ZKRUxkPz6e+acDv4GVSUYNv1JwFGpUfiJ4GWa/FfS8vPpCiOw/cOJgppw3ptvD79O/nCnnjUl8JTEZ0fZO44VHv8djoIIErtkPJtHz78YkFvCWIQPri0QfyuJAl/qVgkML9wRn8jnnR+zxR+Nnqd4tj99J3jeaZAUWZhqvLcO41EF7yOfa0fahwLugS/qxlPWCjBXPvnoqg/20NkA36vErBce4ceOYNm1at4dfWVnJtGnTdH/fAzsVT4qcPzL8LNW75fE7yfuGYwcW7vzsUzCmO7DQz/ZCrvDSn9fy/H3vBU/1ixVt72XonDz7ZKaXuW49eOyr19yIq9ef6f34PEI9fqUgKbSo/GRge9jRnr9fER77/KCeu1tgYTZnAgTx3t9/bRurX97ao91Xqp+XV+0VWOfm2bulCMbjbccTuDbubCuQL1peOBv24/MINfwZIpGyvGBV4Tv99NNZvdrKoX399de55ppr+OSTTxARvva1r3HnnXfy8MMPc9FFF/HCCy90L2X/7W9/Y8aMGfz1r3/truqnKH6I13iHnx/ed83SJcy/7KKIsaLHd6sCmK2ZALZQj52zb3vv9Rsa2LT6c3btaKWidwkGQ2tTJ+Kx7hoz1c9labyjrYjiv12GmFAKX/SSvZtnL8XuhX+CSuo65ez7Oe/038D+x2gFvxSihj9DJFKWN5pPPvmEf//3f2fhwoVMmjQJYwyPPvooO0OpU4ceeigPPfRQt+FfuHAhhx12WFruMx8pxJK84UQb73hx0gZ45u55GGMwnZ3dbW5kayaAm1BPuFff0tTR/b3xqFIcM9XPwavu6hDo6kKiDXh7Mzz2/0IG1WUf3XRa3nX4Q0FxGbTujNT/91tRL959dd2PTylq+H3SWFvL9rl30FFfT0l1NYOuupLKadPiHi+ZZXnvuusuLrjgAiZNmtQ9drgnP3nyZJYuXUp7ezutra188MEHEcp+in8KpSRvsoLpvMZxWsLv6uhwGqYHUlyctZkAyarhAcRO9QsZx/Y/X0ZJeRudrdZnSHG5Sz4/uBt9m5Je1qv5C8vbbmvqKQ6UDap3StxocJ8PGmtrqf/JjXRs3QrG0LF1K/U/uZHG0Ad+vCSrLO/q1as58sgjXY+LCCeeeCLPPvssixYt4pvf/GZC8y5kCqEkb7KC6WKNk8hSfflee2Xt/n4iNTzCGTa6yl+q37iz+WDRQLa+WkVRCZRUGBKq69W8AzqaYcZ8K0Cw+QvnfplWvVPiRg2/D7bPvQPT0hLRZlpa2D73joTGtcvybt68mddff717v94m3rK8TpxzzjksXLiQhQsXcu655yY8XqFSCOI/XsF0yRwnkaX6ll274j431TgJ9cRDw6f+pbdLqqsZeNhOiko8PP0ghEfyu0XTS1Hmi9zEIhsr8mUBavh90FFfH6g9KImW5T344IN58803PfscffTRrF69ms8++4wDDzwwKfMuROIR/2lauZ36Oa+zedZS6ue8TtPK7amaXlJw88SDeuixxnFK7ysqKcGPu5qt+/uwR6invHdxQuME2TIYdNWVlO7lEJTnSYyfs70l4FhRDyseoPYH8OTV2Wlcs7UiXxaght8HJdXVgdr9kMyyvJdffjkLFizgtdde627705/+xLZtkTnAt9xyC7/85S/jnrMSvCRv08rtNDy2ns4G60O8s6GVhsfWZ7XxdzOqQY1trHHGTp7CwV+v6dYGkKIiDp0ylfLefTzH9Zs+mEkOnDiY0nKPECofS/FBtgwqp03DlO/juz+lvawiPRLj4WTuIdbXaXc6921vtlLvstG4ZqjkbS6gwX0+GHTVldT/5MaI5X6pqGDQVVfGPWYyy/Luu+++LFy4kGuuuYbt27dTVFTEcccdx4wZMyL6xYoVUGITtCTvl89uwrRHhm2b9i6+fHYTvccPSvl83fAKunNT6WtraWHN0iXd/WIFADqNE2601yxdwrsv1WG6rJ+P6eri3ZfqXNUBAfoOGJi1qn3RufueHnvYinxRsWAwEVl0JWVFDD9kHxbcsMy3kl/R6XP8l+s97DtW2tyKe7372YZ82p0e6QdR2wvZEviXrRX5sgAxJkl7QlnMhAkTzIoVKyLa1qxZw9ixY32Pkeyo/mwl6M9F8WbzrKWux4bNmZzGmewhOo0OLIM8deblEUZ98YL5Parp2f3AWcgnfAx7HLeHg/mXXeSYridFRd0PA+H0HTCQmXfdl8Cdp47o3P2gVPQuoaS8uNvIDz9kH9Yu3xYxXklZEVPOs1YFX354Ha1Nnd3nTj77QOuhIFy0x61SH1ge/7Q7vVP7wrEL6PjpC4BY2vOZxK3cr61qmOeIyJvGmAlOx9Tj90nltGl5aegVf8Sbu19cVd69zB/dnin8KOGNnTyFpQsf6GH4w4Pz/KjpeeX8e1X8Kykrj0sdMFM45e4HoaWpg8tuP677/YIbljlqASx9+H1aWzoiVgdamjqoe2ANAAdODMt/d61zzx6vvOZGWHQZdLZ5T7BxsxXl73dFIRvkdTNU8jYX0D1+RYmBnbtvR+7bufurVq2KeW6/k4cjpZH/ZlJaRL+Th6diqr7wG7zn1S8ZAYCuMQADBjJ15uX0HTAQRLrfOyn+3X7ONOZfdlHGdfu9lvXDq/FV9Hb2taSICF1+t/FamjochfW6Og2vLtoQ2egWlGfTuNl6SJh+F/Tq794PLOW+mOVzQ2SLcc3minwZJqs8fhH5BjAPKAbuMcbMiTo+BrgPOAL4kTHmtvTPUslnnDx7r9z9WF6/vY//5bOb6GxopbiqnH4nD8/o/r6bDG60IXatpte7D2UVFb7G8MIrBsBrpcBJ8e+5+b8DyNjev9uefp/+5Vzwy2O737ttCZguWPLgWsAKDIwZI+BAj/62gfvbpc4yvLZXHq6St+ph5xWA5h1w6wFwyq2hErJVuG4lpMK4etUd8EIVAB3JGo9fRIqBu4BTgIOAc0XkoKhuO4AfAGrwlaTj5tknmrvfe/wgqmcdzbA5k6medXRGjT64V8mLXkqffM75SHHPSO72lmZGjD/K1xixKCkv6/6+vE/fHp69E8nSGQjK+69tY8ENy7jr0sUsuGFZhIfulLtfUlbUQ3nPTvVz0ufvaOvi+fveY8ENyxh+yD6O43mlCDpmAYw7G751d0/Pv7QXjJraMw1v3NlQ5pJV0bwDHpsJsytxLTAgxVafZKb1aVpe0skaww8cDXxgjNlojGkDFgLTwzsYY7YbY94A2p0GUJREcPPs3fDK3c9m7BK7Xkvpdr9yh3TSro4ONq58w9cYbthee3gMQWdbjH3mEMnSGQiC7am7lc21DXr4sv6U88Y4RuEfOHGwpz7/rh2trF2+jTHHDI7YGiguFUYdua9jVl1RsbjL+zoteR/2HXjnz87G1E2pD+j28p1WELrbo8ZzE9HxK66jaXlJJ5uW+ocC4ZEom4GeGrZKVlO/bREbN9xGS2s9FeXVjBh5DdWDp8c+MQsIor7nlbufapKho++0lO40rptC3s7PP4sYwz736bt+42tOiZTa9btVkUzcCu+El809cOLg7u/t1L7n73vPMRUv1lJ+R1sXH7y5nY6wVNDWpk7WLt/GwccOYf2bnzhH9bsRveQ99xB3Y+pS8c8RKbb2KaSo58NAezP8/TpL/je6/O9Hy60Hj+h2e67haFpe0skmj99J0iLuXEMRmSkiK0Rkxaefulf4yiQNDQ2cddZZjBkzhrFjx/Lqq69GHL/tttsQET77rKcns2nTJnr16sXhhx/OQQcdxKWXXkpXV1dEu/164IHULoHa1G9bxNq1P6KldStgaGndytq1P6J+26K0XD9R/HrwIsK0adMyUpQnWTr6fsd1E9MJN7LxzCkRr93vVkUycTPSTu2xVgfAR/EdrEA+p4eNTas/55Lbv85l//kZl43+Ad/v+00OfOXEYEvfXsY0VlBgOKbLSttzW8Jo3uH8gLHiXv9efK+9ncfOROZAnkgAZ5Ph3wyEh4sOA7a69I2JMWa+MWaCMWbCwIEDE55cKrjiiiv4xje+wdq1a3nnnXci8uc//vhjnn/+efbff3/X80eOHMnbb7/NqlWreO+993j88ccj2u3X+ecn9wOxftsili2bTN3ir7Bs2eRuw75xw210dUX+M3d1NbNxQ26EZDip8jlhjMlYJb5U7W+7jStCTCMbz5wSUQf0u1WRTNxU9Jza3VYHXljwXsTWQLzs2tGa+L63m9GsHLZnayBWpH/4OIGNsItP1/hx5D2setgqCRxNcVn6MwfyKNYgm5b63wBGicgBwBbgHOA7mZ3SHp7a+BTz3prHtqZtDO49mCuOuILTRpwW93hffvklL7/8Mvfffz8AZWVllJXtCXS66qqr+NWvfsX06bGXyUtKSvjqV7/KBx98wBFHHBH3nPxge/W2gbe9eut759oFbu3ZRrQqn4hwQPsgjuoYSR8q2EULb5Rs4PN9WmKMlDpStb/tdn7Lrl2cetnVnlsL8cwplqpfLLyi/lPBpOkje0TjOwXvgfvqQHTkvhclZUUUlwqtTZ3s+8nrjNj4BBWtX9BSvjdbDj2T9qfvodTNY/YTxe6U4w5WCV7bkJX1tjx2KQ4t4wsRBjs8bc8tZ76kV8+SvrEIX/Kvuxm6HOJsyvr0vM94I//94hVrkGOZA1lj+I0xHSJyOfAsVjrfvcaYd0Xk0tDxu0VkMLAC6Ad0iciVwEHGmC9TObenNj7F7Fdm09JpfeDXN9Uz+5XZAHEb/40bNzJw4EAuuugi3nnnHY488kjmzZtH7969eeKJJxg6dCiHHXaYr7F2795NXV0dN99sLZNt2LCBww8/vPv4b3/7WyZPTo5KnJdXX1FeHVrmj6SiPP6aBulm3Lhx3Q8AaxetoOLVXZRgRVP1pRfHdYylZYS3lnwqSdX+tte4sYys27ki0r3c7/bgkGisQrqwDXW4JK+bhK7X/n14XIBbPymiW6Fv9a8f4MB1f6Y4ZPx6tX7ByHf+j5KRW503R/3ue9uG6u/XRRrm5h3w2P+L7Gs6LSN+2Hdg/XPWNezl98dmWoZv1FTLyNuGsVd/K/UPeo4Xi3Bj6nY/0QGItjfuJ2YgXvIo1iBrDD+AMeZp4OmotrvDvt+GtQWQVua9Na/b6Nu0dLYw7615cRv+jo4O3nrrLX77298yceJErrjiCubMmcP111/PL37xC5577rmYY9gGXkSYPn06p5xyCps2bepe6k8FXl79QQfdHrEaAFBU1IsRI69JyVxSTeWaLjqJDKEuoZjKNV1R+SbpI1FPOei4I8YfZUnrehhnN21/09XF07+7PaItOuc+3uDATBAevOeF0+pAOLaxd1tFCM8G6Kx/uofHK+2ttDWXUr6XgyfstOTu5gnbHrUfj7y9Gd6830oNhJ5GdsUfI/s377AeKk651XoICOr1N35spQ26EX2f6fDG3YIes0GlMCBZZfizlW1N2wK1+2HYsGEMGzaMiROtxIWzzjqLOXPmsGHDBj788MNub3/z5s0cccQRvP766wweHPmhk0oD74aXV29H7+dqVH80TlK7Xu3pIFWestO4I8YfFVE0x00ox/7+77+f66ixH0109H42CvIkgm20X1jwnmPM2/67VrL+hF/QWV/Pcf0H8cEB0/ioz3jnVYQdzlUcP327N/se00xp0Z6/xfauckqj971jecJBvFW7DG+4Z+9F8w6rv5066OccPzgpA6bDG88jCWA1/D4Y3Hsw9U09Pd3BveMP0Bk8eDD77bcf69atY/To0dTV1XHQQQdx6KGHsn37nn/24cOHs2LFCgYMyI764yNGXtPDq4dSOjqbqVv8lZw39uFko84+pG5/O3rc+Zdd5DvlbuzkKTx91298Xyt8/z+R1L50EV55T4qs/Xqv5X67Ldqjr/5sBV95/8977vfzTxjV9CeO+9mBVE47occ4JdXVdGzt+aC9/ZOhvPflSRzT50H6Fn/Gzs4BvM2FHBft3cbyhIOk7tnnBjHg7c3W9kB3QaDNzql/fqncz3nvPh3euH3NVMYRpAk1/D644ogrIvb4ASqKK7jiiCsSGve3v/0t5513Hm1tbYwYMYL77ktO5bHoPf6LL76YH/zgB0kZO9qrLymupLOriY4Oa88tPNgv141/v5OH0/DY+oiyupnW2Q9CvPn+3ec57NuDe9Ce216/W99Y46VSkCcI0TK7thdvp+mBc7DegRMHU7+hgXf/sbU71f3Aj5+E6C2Rlha2zrmdx1/t3yN+wKkkOGXlbBp1BvUtE1jf8nUgsnJfBK6e8MdWSlqvva0I+VhFehLBrgkQLgvcI7AwKnDQEXGvqpcubzxPJIDV8PvA3sdPZlQ/wOGHH050ueBoNm3a5Ng+fPhwVq+2/gna2htoa91GV1c7AwaU0vhlPWWlVQnNzYvqwdO7jfqyZZPpaG2IOG4H++W64c9GnX2/OC2fP/2723n6d7e71rRfs3QJdffPp3WXQ/pUGG6BhCPGH8U7zz/teCyc6JiETAjyBMGr8l60iE8477+2jbXLt3U/KJguKGp0fjAyn3/SI+8f4MBQRdDokuDFg45ip49AQ2+P3ljL8UWlYfvwfgxwQKK9bi/P2auioJf3nkfeeDpQw++T00aclrChTwVt7Q20tGzpdkO6utppaf6Yzo4mevUamvLrpzuFL97yuPHSe/ygnDD00Tgtn9s47aFHPyi44RZIuGbpEt59qc7X3IrD0lbBOThQiotpb23h9nOmZTzYL1axHLfjTg8MLeV706u1pyRuS3mkSE34A4VTSfBKfGoBuKXthdPVbqXuXfdhZCBgr72hbVfs1QB778MJN6/b9pzt64VnB6z8v57XLCp1996jgxdnzFeDH4NsEvBR4qCtdZvjP117+w7a2htSfn23VL1UpPAlUh630Ii1TB4tsOP1oGDjJZTjdr4U9fyIad21M0LZL1qQp7xPX0TE0vFPojphvLiJ98Q67vRAsHHEN+ksihSJ6iwqZeOIb/boW737OdrnjE1MJS5ap98Ne0tg3Nmh6nsN1oPA9Ltil+GVYmu7IJpe/b0r9TkJ4rzzZ9h/UmQRoNLecMbvncfJI1GddKKGP8fpchK3CNHWGn/WgV9GjLyGoqJIec9UpfB5lcdVIvGzTB7+cBDrQaHvgIHMvOs+V6/b7XxjjGXQo+hoa+Xpu37D7d8+nd+c+83utL9TL7uasooKujo6evRPdfU9N5wq79m4ifiA8wPBJ/sezYdHXEDJkCEgQsmQIawd/R0+2ffoiH6jKl5iSr8/UNpiyV8nZNDCjbmbEXdbRrfPdaoMZNPVbgnqhBcBmvG/1oODm7Gee4iV3+8UePjhy1HOjEemiBbwiQs1/DlOUZG7xKzXQ4Ff3OR5baoHT2fMmF9QUT4EECrKhzBmzC9Ssr+faHncQsJPXn/4w4HXg4IfnQAvCV7Xhwpj7SXbKYC2Zx80qDDVRFfes51Rrwp84PzAIMWwfd+jeO7AH/H6t/6XrZfc1cPoAxzT58GIVD0gOQbNSYffTxDckRd6H2/+Ys/DxVWrfXr5bkTFGDjdt/3w4DZODorqpBPd489xysoH09Ls/Mfv9VDgBy953nDDHh7sl0oqKysdjXyulsdNJWMnT2HLujWewXbhxtxNhKeib19OuGBmzP11LwGgVXXPYIy/gLGOtlakqMhRDyCTwX5+xXuiz4E9an/lvYtpb+mipclazdi1o5XVLzuXI+lb7PKQE2bQwlMMIwL8vKRr4w2COz2Uqhkt1GPjN23OyUP3Q7ghf/Jqq8iPVxBiDorqpBM1/DlOWWkVnR1NtLdHKWNJEWXl8esMgLc8byYi9mtqaqitrY1Y7s9kedxs58RL/pO1ry51jNKv6NvXUYQnXmGg6PPLe/ehq6PdV5R/NKari5Ky8qSrE2aC8AeGBTcso7XJn/jTzs4B9CvpufLRUrIvf7lhWY/4ATsToO+WWqrX/MRbujbelLTTfwP7H5NY2ly8nrhtyFc9HNvo56ioTjpRw59BGhoauOSSS1i9ejUiwr333sukSZO6j992221ce+21fPrppxECPv/85z/53ve+B8BHH31EZWUl/fr1oX//Ptx55085+ugzOPDAUUhoX+7qq6+Oq0JfthXdiS6ik46o/lyn5sKZjp74CRfM7NE3qDCQk07AzLvu850h4Iadbuj1EBKvRkG6cPLGY2UHhLN813nWHn+UMt/Ln53Drhb3OgB93/kVSAqlaxNNmwsqGGQzauqe63p6+i4CP0oEavgziF2W95FHHqGtrY3du3d3H/Mqy3vooYd2S/VeeOGFnH766Zx11lkAIa3+r/DOO/9MeH7ZWHQnvIhOLtO0cnu3PoCdOp0KnQAvTz58DkGv7SWz6ydDwA3bs/d6CMl2id9owR/bG6/oXdK9zO/FqIqXOKbPg5RIK12mCKGLnZ0DWb7rvG7BHjd64yKglMw9b7cVAz/V8fykFzqxPlS7xOs+KvdzF/hRIvBt+EXErTC8AVqMMf4ku3KVJJd8TGZZ3lThJM+by0V3soWmldsjFQFDDkxnQysNj60HSLrxjzaI0XMIem03mV2/ev1OuAkL+b12tkj8OuXvd7R1UVwqlJQVuYoBQVg0f8jTF7po7yr3ZfQBmhhIHxz0/YPsecfzWee3Ol7EikEAz982+K4rBqLL+wEIEtW/CfjQ4bUJ2CYiX4jIb0Qk/1YRUpArGl6Wd/z48VxyySU0NTUBBC7LG40t2Wu/li5dGtc46YzYz1WaVm6nfs7rbJ61lPo5r9O00rmoSjhfPrspQgY4HNPexZfPbkryLP3NIci1XdP34jD6JWXlnHr5f3umC/q5drZI/Lot6bc2dUZkB1T0LqG4LDK33imav7SolWP6PBjzuiVlRew87IfxRe3bxPtZFyStzk4RjKUPEI794OKUlYDAhIt1eT8AQYz0ucCvgLuB10JtE4GZwGygCvgxsBP4adJmmA2koORjMsryuhG0al+45G9r62fUb3u/27inK2I/XSSyvO00ll+vOWJpPwbpqP4XtPJg9J56RZ8+lsBOnEhRkZXjH8f+fCISv421tT3kb6NV8RKlT/9yR+Pfp395d7CfvR3Q2Ra5X+0Wze8a5R9CimDMMYOpnnE8rKqMf3Uy3s+6eKrjOS37F5WCSKRyX/iDi0rzJoUghv8/gKuMMY+FtS0WkXXAFcaYr4vIduAm8s3wp6DkYzLK8noRbsyLikopKx/sqN8fLflrTGfeFNmJJtHl7Wi8vObw8Xos7ccgHdX/glQedNpTl+JiikpKegjthOOWlldSVu6qAOgHt9TBWFH/jbW1EQVvOrZupf4nlkFJpvGfNH1kj6p80UI/bvr/btH8OzsHhMYRSkqLe8QKmC5Yu3wb1SOrOHBiAoVk4v2si6c6npsRd2qL3i7IpKFP8rZvJghi+CcCThFjq4GjQt+/CuRfAmUKSj6msiyvMZ099ftbtnQfD38gMKarh+RvvhTZsbH1/b/xycH0JXKZ0MlQ+8Wv1+y1tB9Nuqr/9Tt5OF/8dV2kKFoRjtd22lM3nZ2U9elL2+4m1+V9Oy3v4K/XsHHlG92rBSPGH8XShQ/w9F2/icvjjzf1cPvcOyKr3GFVxts+946kGv7o/H2nIjpu2wFu0fzLd50HQEebAbocAwW9Cgb5Jt7Punir47kZ8Ww1pH5jGbKcIIb/X1jL+tdGtf8/4KPQ9wOBqITyPCBFJR9TVZZ348YP+dqxZ3a//+53z+DSS8+jtWUrBhPxQOBGplL2ko2t79/e3k4fKhz7xLu07tdrjjl+CqP6va8bVYlNnLXc3fbOW5t2ceplV3um7nW0tbJx5RvMvMv6205WRH7Q1EOAjnrnv2m39kSIJfjjth1gB/Ad0+dB+hZ/xs7OAT0C+zraulwDBIOkDDqSiAGHnPeEY5KCbd9MEMTw/zfwqIicCryB9YlxFDASsK3MUUD+VUdI0R91ImV5beysAJvhw4fzySdvOPY1ptP33DKZspdMwvX9d9HSw+OH+JfW+508vMcSvpPH7vaAIL2KGfrTr8Z17UT58tlN0BmVD91pHFc/vPbUI7xvH1K78UbkJyNvv6S6mo6tPdNTS6rT/7futB1gs77l674i+J2IVVAoJol81mV6CT4dpGDbNxP4NvzGmKdEZBTwn8BoLD/lCeBuY8xHoT6/T8kss4Ec+qMuKipNSKc/n1L2wiV+3yjZwOSOsZSyp+BIIkvrtoGMFSzouKwOmLYumlZuT4mHH8tQBgnuc5PzHTHe2uGzve/5l10UM+gunoj8ZK0SDLrqyog9fgCpqGDQVVf6HiNZOG0HDD9kH9Yu3+aZ7mdT3ruYznbjGUcQNzn0WZd2UrDtmwkCpd4ZYz4Grk/RXJQkUVY+OGKPHwApQhBHr1+kGJEiurraESlmzJhfALBs2WRaWuupKK9mxMhr0rLnX79tERs33Ja064br+28s+QSAozpG0ocKSqoqEl5a7z1+UMzze48fRGPtBrp2RwXCuXjYieLHUAYJ7nPT/X/3pTqGjh7bPWasoLs1S5cgIo66/V4R+cnK27f38VMd1e+XaOO/afXnjDlmMJtWfx5zyf64s0dHnGs/OLy6aAPP3/eeY1yBkgRStO2bbgIZfhHZCzgcGESUBkBUtL+SQezo/eiofsDxgaC8Ykj3OeXla4AvfRXnSTZ+iwIFIVrff2PJJ3zcawfTpk1LqwJgD6Mfwmv/P97UQz+G0u82hc3GlT23j6LH9Aq6sx9G3KL8vSLyk5m3XzltWkYMvZOEL9BD4W/t8m3dFf8WOGjyg+Xt2wbd/uqmFhjeR0kCeRLLEES570TgIWAfh8MGwtZPlYxTVlrlmL4HPR8IovtlqjhPKq6bKn3/oEbZy8NuWrk9YkVAehXT67CBNL+5Pa7UQz+G0u82RZAxwT3ozk3GV4qKYqb2JZK3nw24GeWS0p4qfuGR+W5pgba3H46bWmDCUf5KT/JgKySIxz8PeAq4wRjjXEtSyXq8HghsEinO42ep3q1PqooCJVvfPx49ADcPu3zM3nzxyPsRgXamuZPdy7f1GMNv6qFfQ+lnmyLomG64Kv0ZE3fJ31yp1udmlGNF5juV9RWE5+97j1cXbYhYynfbGkg4yl/JS4JI9g4HfqZGP/9xi+iPFelvL9VbhX1M91J9/bZFvvrEe910E4/cbe/xg6iaMap7D71orxIoEcvAR0fXe+An9XDyOedTUha5V5+ooUxkTHtv3wk/Dw5jJ09h6szL6TtgIIjQd8DAhASA0k08xvf916wHvwMnDuaCXx7LIccNobWpszt33141sPu5RfMnHOWv5CVBPP5lWNH8G1I0l4Ij3rK8YKX5jR07ltGj9yz72eV3hw8fTt++fSkqKmLfffflgQceYPDgwb7FgOItzuNnqd6rT7YUBbIFf9y2BoLK3drYHnZQJb9w/KQexitwE3RMP0I8ieztR18/Vwx9NG45++W9i2ltck6xXfrw+xGevlM/P9sCSYnyV/KOIIb/buA2ERmCpeAXkS9mjHkrmRMrBOIty2vjpcm/ZMkSBgwYwA033MAvf/lL7rzzTt/zCjfSLa31lBRXggjvvfff3Qbaac/dz1K9V5/o66Yzm8AmXPAHrHTA2tpaYE+8QJCIeCeCKPmFEyT10DaU3bEIT7VSv+z1hLIYwo2v3xS7RPb2U0UyNAGC4LVX//x97zme09LUASHv3u3hACzP/65LF9Onf3lERoBrVH8eyM0qiRPE8D8S+jrf4VjeB/c9vnILv352HVsbmhlS1YtrTx7NGeOHxj1eusryHnfccYGMvo1dnKc70r4zdqR9RXl1aAmfHu1++2S6KFC44I9Ne3s7dXV13YY/aER8NL6UAotgr6MH07r2C1/Bd07BhkBSaxOE4zfFLpG9/VSQLE2AIHhJ+NptiRKdEeBInsjNKokTxPAfkLJZZDmPr9zC9Y/9k+Z268l7S0Mz1z9mlS2I1/iHl+V95513OPLII5k3bx69e/f2XZbXLr9r89vf/pbJkydH9HnyySc59NBD45ojBIu097NUny3L+W6EC/64tQeNiI/GbcUgnL2OHkz/M0b5Gs8t2FBKi3wVEYoHv1H+2RaRnyxNgKC4Sfh6KfgFJWYUf57IzSqJE0S571+pnEg28+tn13UbfZvm9k5+/ey6uA1/Msryei31T5kyheLiYsaNG8fPf/7zuOYIwSL8/SzVZ8Nyvhfhgj/R7eEEiYiPxmnFIJrmN7fT9G+Vvq7hFmzoNP6/dr7Lqo9eZvc5OxNa5vZr0LMtIj+ZmgDJwGk1IJEVAM9z80RuVkkcT8MvIjOAWmNMe+h7V/JZwGdrQ3Ogdj+kuiyvvcefCFY0fhHQc4/RLdLez1J9ppfzvYgW/AEoLS2lpqYmadeIXjGwi/SEE8Qz91tk6F873+WNz5+h01h7x4ksc/s16KkINEyEbFuBgJ6rAff890ue+/peeEbx54ncrJI4sTz+R4DBwHb27PE7kZQ9fhH5BpZeQDFwjzFmTtRxCR0/FdgNXJiOoMIhVb3Y4mDkh1T1LPjil1SW5U0G9t6+k9HPpqX5ZJMqwZ9owlcMNs9a6tjHr0H3KgJEh+n2/Fd98XK30beJd5k7iEHPpoj8bFuBcOK4s0fzwgPvEaCmFuAjij9P5GaVxPE0/MaYIqfvU4GIFAN3AScBm4E3ROQJY0x42OspwKjQayLwh9DXlHLtyaMj9vgBepUWc+3JPRW0gpBoWd7oPf6LL76YH/zgBwnNycZpb9/C0vLPVo89GSRD8CeIsl+iWQJuwYZV3/wKsGdlYXfnl47nx7vMnU0G3S+ZXIFwku112o93Wv7vaN2Twx+OFFkK3L60+fNEblZJnEBa/SnmaOADY8xGABFZCEwHwg3/dOABY1X5WC4iVSJSbYxJafF4ex8/mVH9kFhZ3uHDh9Pc7LzV4HZOrBK/4bir5XXljdGPla8fL17KftAzKDDRLIFYwYb2176X3RtzmTtWqlv38c8+RYqKMF1d9B0wMKPL90HJxANLUC396OX/6PPB8vA9o/idyAO5WSVxghbp2Q+YjHORnt8kOJehQPgG1GZ6evNOfYYCKTX8YBn/RA19LuEnNS8b8WvM/eTrx4tbsF3DEx9ELL3bDwRVM0ZRNWNU3FkC4C/Y0GuZe83SJSxeMJ+WnTu7j0XHAESnwtmiPOlIict1EtXS90oJVJSgBCnScx5wL9ABfEpkOJIBEjX8Tpqe0VqmfvpYHUVmAjMBTxEcxZlsT7tzIogx95OvHy9ue/OmueemrR3EVz3r6KSX543GbZkb6PFAYBMeA+AmxhPdT+lJMrT03VICFSUoQTz+m4HbgZ8Yp6LuibMZ2C/s/TAg2uX00wcAY8x8QmJDEyZM8C+GrgDJT7vzU7zHDb9efBBj7idfP1785OmHE6RvIjSt3E7Vst6c2vdiivfbs6ow/7KLXA067IkBiBULkKmUuFzALU0vm7X0/cYkKLlHEMO/L1akfSqMPsAbwCgROQDYApwDfCeqzxPA5aH9/4lAYyL7+8YY1+IhhYgVOrGHZKXddav/dcVW/4smiBcfxJj7zdePh34nD+eLv64Dn5osbkF8QUv/euEVdxDLYNsxAOW9+9C6a2fMfkpPck1LP2hMgpJbBInUf5oURtAbYzqAy4FngTXAw8aYd0XkUhG5NGwOG4EPgP8F/jPe61VUVPD555/3MHaFRlt7A7t2raWxcRWbN78D0pD0a3ip/8XCy4uPxs1oO7XX1NRQWloa0ZasfP3e4wdRVOH/mdopiM821PZqgG2om1Zu79HXD14VBb0MdngMQHuLu25FtqXEZRsHThzMlPPGdHv4ffqXBw/MC+P917ax4IZl3HXpYhbcsKy7Sl+y8IpJUHIfPwI+Ns8Dt4rIwTgX6UlYwMcY8zSWcQ9vuzvsewNcluh1wBLQ2bx5M59+2jPKuVDo7NxNe3sDYDB00dnxEU27/0hV1aykRu4HUf+LJogXH0R8J9X5+l27e6ZeBcHLUMfj9XtVFHQK+gMo79OXmgtnMnbyFOZfdhFdHc73VN6nLyLw9F2/YenCB3Iqwj+dJGuPPh3eeDJiEpTsxY+ATzQ3OLTlXJGe0tJSDjigYMsPALBs2WTHyH0nHf5ESCRDIMiSfFBjnox8fTeC7PM7GXMvQ73j8fXsfm2b9V8nsNfE2Lr+XloBY0P1HbzS+Ly2Azrb2tJa9KbQSTRDwA+5GJOg+Me3gI+SfyTiiQchkQyBoBK6qTTmQfCjx2/jZpAdjX+psHt52LKuofu9l/GPpRUQK7fdTepWiooyUvSmkEmHN55rMQlKMLJJwEdJM+nK1U8kQyBVS/KpEu+x6aHHX4yT+jEQktaNws1Qmw7nB4ndr23zNPyJVhR00wBwywbQCP/UkQ5vXHUD8psgefz3Au8aY26Par8aOMgYc0myJ6eklnTm6mdTYZ5UivfYREfkd7V2OObxA5i2LppWbo8wwm6G+ou/rHO+oIH6Oa97GvVEKgq6aQDYKn7RaIR/6kiGN+4nVU91A/KXIB7/qcBvHdoXA9mr6qK4ku0lciE1RjqV4j3gnDrnSadx3Od3MtRfPLzORbKKHhkA9hjJwm07INuL3uQbiXrjmqqnBDH8VcAuh/YmoH9SZqOknWzyxJ1IhZFOpXgPOEfkx8JvIOBeEwdH7vG7YEsEZ0oNUPf3U0si3ng6ggOV7CaI4X8fy+ufF9V+GlZevaIknVQY6VSK90B8Snx+K/HZ+/jhUf1uKwCmubPHFkIqyMUqfYWMpuopQQz/7cDdIjIIa3kfoAa4kiTl1iu5SSJyvLFIhZEOmikQFLeI/KK9SjDG9NjrD1KJDyzjHx7IZ+/tOxFv3r+Sv2iqnuI7Xc8YswDLyJ+PJebzPPA94GpjTLBC8kreYMvxWtkBpluOt37boqSMnwqFvXHjxjFt2rTuh4fKykqmTZuWtMC+ficPR0oj/7WktIjKaSMZ+tOvsve3R3d7+MVV5VTNGJWQcfZ6aEhXHQAld5g0fSQlZZF/n5qqV1hIPJK1IjIwdG58+qFpZsKECSZW3XslPtxEgCrKh3DssUuTco1Up96lgnhEdhJh682vuqoF7v3t0er1KxFoAZ78R0TeNMZMcDoWVx6/MaZwdW6VCNIhApQtojx+aVq5neY3t+/ZezfQ/OZ2mv6tMmUGuHLaSNdUP13uV6LRVL3CxvdSv4j0F5E/iMj7ItIgIl+Gv1I5SSV7cRP7SbYIUC7hpbOfKrwMuy73K4oSThCP/4/AeKwa99aGrlLwpFMEKFfw0tlPJV56/IqiKDZBDH8NcJIx5rVUTUbJPXJBBCjdxDLA0ap+QaRzvYilx68oigLBDP92nAV8lAIn20WA0o2XAXZS9UuWyl6ievz5wFMbn2LeW/PY1rSNwb0Hc8URV3DaiNMyPS1FySqCGP4fATeLyAXGGH0AUBQXvAxw/ZzXXff/k2GgE9Hjz3We2vgUs1+ZTUtnCwD1TfXMfmU2gBp/RQkjiOH/MTAc2C4i/wIidFSNMbkTdq0oKcbNAGdq/78QmPfWvG6jb9PS2cK8t+ap4VeUMIIY/kdSNgtFKRA0AC91bGtyrmHg1q4ohYpvw2+MuSmVE1GUQkAD8FLH4N6DqW/qqR8xuLfmqytKOL7z+BVFSZze4wdRNWNUUiV7FYsrjriCiuKKiLaK4gquOOKK7vdPbXyKqY9MZdyCcUx9ZCpPbXwqJXNJ13UUJR58e/wiUoYV4HcusD8QIaBujClO7tQUJT9JdgCeV3pgqlIHsxF7H98tqj9dwX/puo7K7irx4lurX0RuBb4N3ALMZU+w3znAT4wx/5OiOSaMavUrYGn+//3vf6e52RIb6tWrF6ecckpOyQFHE50eCNbWQdUMqy6A27F8Nf5eTH1kquNWQHXvap4767mcus77r21jyYNr6Wjb87stKStiynlj1PgrgLdWf5Cl/rOBS0MGvhNYZIz5AfBT4KTEp6lkK/XbFrFs2WTqFn+FZcsmJ63yXjpZtWoVjz/+eLfRB2hubmbRokWsWrUqgzNLDC954ExIB2cz6Qr+S8d1Xl20IcLoA3S0dfHqog1Ju4aSvwSJ6t8XeC/0/S6gKvT9M8CtSZyTkkXYZXdtSV677C6QdaI9XlX86urq6Orq6nFOZ2cndXV1Oev1x5MeWKipg+kK/kvHdXbtcP4durUrSjhBPP6PgCGh7z8ATg59PwlodjxDyXk2brgtQocfoKurmY0bbsvQjJxZtWoVtbW1NDY2AtDY2EhtbW23N2+3O+F1LNtxSwMsrir3PFaI+An+y5Xr9Onv/Dt0a1eUcIIY/r9h6fUDzANuEpEPgfuBe5I8LyVLSEfZ3WRQV1dHe3uEphTt7e3U1dUBUFlZ6Xqu17Fsp9/Jw5HSyH9jOz3Q61ghctqI05j91dlU965GEKp7VzP7q7OTHthnCwkVifWzT8V1Jk0fSUlZ5O+2pKyISdNHJu0aSv4SJI//+rDvHxGRj4FjgfeNMU+mYnJK5qkor6aldatjezbh5rXb7TU1NTz++OM9lvuLi4upqalxOjUn8KPPXyhR/X44bcRpKVPxi47m7zJd3Z5+sq9pB/BpVL8SD76i+kWkFPgTcIMxJueiRzSqP36i9/jBKrs7ZswvsmqPf+7cuY7GX0QwxlBZWcmoUaN499138yqqX0kNtude31RPkRTRZbqo7l3tacTTlTWgKH7wiur35fEbY9pFZCpwfczOSl6RK2V3a2pqqK2t7bHcbz/YNjY28s477zBt2jQ19IonTp47xM7HV8lgJVcIEtX/GDADyK6oLiXl5ELZ3fDo/cbGxm5PPxx7z18Nf2ERtFSvU7EfG6+iPyoZrOQKQQz/R8CPRWQysAJoCj9ojPlNMiemKEEZN25ct1GfPXu2Y59cjuBXghOPil4sD93t+BVHXBFxLUhN1oCiJEqQqP4LgS+AccDFwH+FvS5PZBIi0l9EnheR9aGve7v0u1dEtovI6kSup+Q/IhKoXclPvEr1uhHLQ3c7nuqsgcbaWtafUMOasQex/oQaGmtrkzKuUngEieo/wP5eRPqE2nYlaR6zgDpjzBwRmRV6f51Dv/uB3wEPJOm6Sp7iFrTqV6I6VygkLf54iGff/YojrmDW0lmex91IVdZAY20t9T+5EdNiPcR0bN1K/U9uBKBy2rSkX0/JbwJV5xORK0XkI6ARaBSRj0XkKkncjZoOLAh9vwA4w6mTMeZlYEeC11IKALfc/FzO2Y/G1um3lfg6G1ppeGw9TSu3Z3hm2YObd+7l1Z824jS+PfrbjseEzKwYbZ97R7fRtzEtLWyfe0dG5qPkNr4Nv4j8CpgN/A+WNv9JwN3AjSQu2buvMaYeIPRVXRYlIWpqaigtjSggSWlpaU7n7EejWvyxiVdF78fH/JjKsp4PiQbjuU2QKjrqnQWz3NoVxYsgwX2XAJcYYx4Ja1ssIuuwHgZ+6HWyiLwAOD1m/yjAHHwjIjOBmQD7779/Ki6hZDHRUf7R2v35QDw6/YVGrFK9XnzZ9qVjeybS80qqq+nY2lNIq6Q6PiGtxtpats+9g476ekqqqxl01ZW6ZVBABDH8AE5lzFbhY+XAGHOi2zER+UREqo0x9SJSDSS8VmmMmQ/MB0vAJ9HxlNwjPMo/32hauR0EcPjLLlQtfjf87Ls7pfz1K+tHY1vPLJB+Zf1SNVVXBl11ZcQeP4BUVDDoqisDj6XxAkoQw/8AcBkQvUb2H8D/JTiPJ4ALgDmhr7lX91XJWbyq+mUj9t6+k9EvZC1+J/zk8Lul/Nla+9FkIjPENsjb595hef7FxRF7/EEMtle8gBr+wsDT8IvInVF9vysiJwPLQ20TsSr2PZjgPOYAD4vI97H0Av49dP0hwD3GmFND7x8CjgcGiMhm4KfGmD8meG2lgLGr+tmKf3ZVPyBrjb/T3j4AAlUzRhV0VH+4oa8sr2RX2y46TAfgnMP/1ManuOEfN3Sr89m4CfgANLZmRgvCNsrR3vrWH17H7rfeovqnP/U1jsYLKLE8/kOj3r8Z+vpvoa/bQq8xiUzCGPM5eyr/hbdvBU4Ne39uItdRlGi8qvplq+F33cM3FLzRD/fcG1obevQJV96z+0cb/VhkUonPyVvHGBoeWsheRxzhy2NPdryAknt4Gn5jzJR0TURRMkGsqn7ZSHFVuaPxL9S9/fCCOn6wg/O8pHndyLQSn5dX7nepPpnxAkpuEiiPX1HyjVzM9+938nCkNPJft1D39n++/OfMWjrLt9GHPR570Oj8Iili+lemOwYKPrXxKaY+MpVxC8Yx9ZGpPLXxqUBj+8XLK/e7VF85bRrVP7uZkiFDQISSIUOo/tnNur9fQASN6leUvMKpql+25/vby/mFrtj31Man+Mu6vwQ6J9xjdyuq40aX6WLRB4sYP2h8hPH3qgcA8aUSujHoqivZ+sPrwEGBMshSfeW0aWroCxjJNwlTJyZMmGBWrFiR6WkoWUquRfUrFlMfmRrIcFf3ro4wvNEGO8g4z531XFzzqCiuSFi/v/6mm2h4aGFEm1RUBPbaNZc/vxGRN40xExyPqeFXFCUXGbdgHMYpp9GByrJK/nHuP3q0h2cB+B1LEG6ZfEvg88KJfggJSqJGOzqXH6yHh8pvncGul17Wh4E8QA2/Gn5FyTuCevxzJs/xNLR+x+tV3Ivmzmbf13UjGd5/vA8A60+ocYzsjyaelQQlO/Ay/BrcpyhKzhAeRLe7fTelRaWxTwoRS2PfSdc/mhIpSYrRh9jlgaOJLstbf9NN1P/kRsuAG9OtwOenXK/fQEAtBJSfqOFXFCUnsPfk65vqMRga2xoxxlBVXoUgVPeudq2qB1bQnVfk/WkjTmP2V2dT3bs6Yrzw933K+iT1nvxmFthL8+FGvuGhhXFX7AsSCKjCPvmHRvUripITOOXdd5gOepX0Yuk5S7vbnt30rKN4D9C9lO+k4md/77X0Pm6Be9CnvW/fI2DQ3k51kPr1KwbkKNzjgh9D7ZTLj0jC2QJKbqAev6IoOYGbdxzdPuvoWTGX7KHnUrufXHwvQ20H63WvGhgY0Gj4rye6+K8nuihrizSqQcSAgnjdfgy1Uy5/1TnfRioif24q7JOfqMevFDSaypc7uOXdD+49uEcxnulfmc7Lm1/ufu8WtGc/NHjl4oevADh69MC3R3+7u5+9atBYW0v9HeFedRcPTSnm834wOGBUv5vMbjRSWurbUDvl8u91xBGa4lcAqOFXCpZcLNBTyDgZ3YriCo4bdlwPo73og0UREfNuEfu2B++0jRCu629jf+9HlCeiol59PVMaqvn28PgMqaPMbmkppqMjYnk+0SwtFfYpDDSdLwCPbtvBLRvr2dLaztDyUq4fUc2Zg/snYYZKOgn38p2orKzkqquuSvOsFD84ldl10+kPF9pxEusJT6dz0wQQhFUXrErdDQUgOnWvc/duTENDj36tvfqzbOLP6NO/nEnTR3LgxMwVFVIyh1c6n3r8Pnl02w6uWfcxzV3Wh8Pm1nauWfcxgBr/HCLay3cimwv0FDpOwXfXL73esW/43n8sT91tO6CyvJKpj0xNmuRuIoR74421tWy99oeO/cqadwCwa0crSx5cC+Bp/FXBr/DQ4D6f3LKxvtvo2zR3GW7ZqKkuuYRTGd5osrlAj9ITt4C78PbwCn4iQn1TPfPemtcdwOeUw19aVMqutl3d6YP2vn+qCvD4xU7tc6OlfO/u7zvaunh10YaYY8WjBaDkLmr4fbKl1dlYbHZpV7KTWN58thfoAWhauZ36Oa+zedZS6ue8TtPK7ZmeUkZxMtrhEfPh+f9gFdsBIgy5Uw7/XiV70WE6IsYNKrqTCrxS+zqLStk44psRbbt29Czh7DWWaWmh/he/jDmPaEEhfVjIHdTw+2RoubNCmGBtAyi5gZc3X1lZybRp07I6sK9p5XYaHltPZ4P1Yd7Z0ErDY+sL2vg7Ge3wwD6nwD2bcEN+2ojTeO6s51h1wSqeO+s5vmz70vGcoOV8Iblle91S+wywdvR3+GTfoyPa+/QvDz5WQ4OnIdeVgtxGDb9Prh9RTU/5DeufTZf7c4eamhpKSyMf4kpLS5kxYwZXXXVVVht9sErxmvauiDbT3sWXz27KzISyhGijHb4PH8tQ1zfVOxpkP1sIfohWHEx0y8AtT1/22ZfP9zsmsm9ZEZOmjww8FuCpAOi2UqDyvrmBGn6fnDm4v2sNrs2t7VQveZsJr7yr3n+WM27cOKZNm9bt+eeClx+O7en7bVf8GWongxxrC8EvXqmC8TDoqisdhXaGzPpvppw3ptvD79O/nCnnjfEM7PPK+fcSDXI7pvK+uYFG9QdgWHmp656+QSP9c4Vx48bljKGPpriq3NHIF1e5L+cWOm6iO06E5+4Hydn3wq/ioF+i9QHCI/Er8Y7gdxqr/kc/xrS19TjmtRrgJigURN738ZVb+PWz69ja0MyQql5ce/Jozhg/1Pf5Svyo4Q/A9SOqI1L6nAiP9Nec/8KmaeV2vnx2E50NrRRXldPv5OH0Hj8ooTH7nTychsfWRyz3S2kR/U4enuBs85dwA17fVE+RFHUH+DkRnQaYaPqel+JgvCRLaKf+ppscjT7FxZ6rAY6CQgHkfR9fuYXrH/snze2dAGxpaOb6x/4JoMY/DajhD4BtuG/ZWO8Zzb+5tZ3L13zUvTWgKwGFhx2EZxtoOwgPSMj42+cm+4Ei33Ey4LHU/JKFm+Jg0C2DVNDw8F+dD3R1eT5YeK06+OHXz67rNvo2ze2d/PrZdWr404Aa/oCcObg/rzfu4oGtO1z3/IEex+yVADX8hYFXEF6iRrr3+EFq6JNAugxysrYMIAViO52dzu0+FF0TWXXY2tAcqF1JLmr4A/Loth0xjb4bbloASv6hQXjZTzINsp9rJTqunUJnL6/bKXRA/Ma/uNjZ+BcXxztNXwyp6sUWByM/pKpXUq+jcQTOqOEPyC0b6+My+uCuBaDkHxqElxskwyCnC68UungNf9XZ/07DQwsd21PJtSePjtjjB+hVWsy1J49O2jU0jsAdTecLSLxee68i4foR/iNeldym38nDkdLIfy8NwstOkimuk0qSkUIXrba31xFHUHXuOXs8/OJiqs49h+qf/jQZU3bljPFDuWXGoQyt6oUAQ6t6ccuMQ5NqkL3iCAod9fgDMtQjpS8awdrrH6ZR/QWHBuFlB07V/MI9/OiqfXYuP9BjJSDTxWwSTaFz2yqo/tnNKTf0TpwxfmggQx902V7jCNzRsrwBia7SF40ae0XJDmKV4gX3yP7wkr7Q02iClb5W/bOb02b8veYAsSPs159Q4/zgMGQIoxbXpXbyCRK9bA/W1oDXKsGxcxY7xhEMrerFslknpGyu2YJXWV5d6g/ImYP7c9vo/Rjmsl9vG/0VXz1Yjb6iZBA/inl+xXWyQaK2cto0qn92MyVDhoAIJUOGdBv9aN38rdf+kLXjj+D9YyZ1L+s7GX3IDbW9eJbtrz15NL1KI4MUkx1HkKvoUn8cnDm4P2cO7k/1krcdA/00el9RMo8fo+5XXCdbJGqdUujWHjPJsVqfaW6ms9nyeDu2bgURxzS9IGp7mSKeZXt7JUCj+nuihj8B3Pb7NXpfCUoqVP4KHT9G3W8ufzIkalNBY20tpqHBX2cHox9EbS+TxJv+FzSOoFDIiqV+EekvIs+LyPrQ170d+uwnIktEZI2IvCsiGZe9un5ENb2KImv2afS+EhQttZsa/BTZiVXS18atME6mjWY8Ww3RWwXpDFD0w+Mrt3DsnMUcMOspjp2zmMdXbtFl+ySTFcF9IvIrYIcxZo6IzAL2NsZcF9WnGqg2xrwlIn2BN4EzjDHvxRo/mcF90Ty6bYdq8isJUT/nddec/+pZRzucofglVlR/EDId1e/EmrEH+VLZs8n2QD6vID7QZfsgeAX3ZctS/3Tg+ND3C4AXgQjDb4ypB+pD3+8UkTXAUCCm4U8l4fr9W1rbuwv0qPFX/KIqf6kjmQI9ySqMk0zctiCcyIYVilh4BfEtm3VCD0OvynzxkRVL/cC+IcNuG3jPzU0RGQ6MB15L/dS8sdP7Nre2R5TmfXTbjkxPTckR3NT8VOVPiUWfrx9nBe05IFVVFFdVOS7rRwv5NNbWpnHW7gQJ4rNXB7Y0NGPYo8z3+MotKZ5l7pM2j19EXgCcyl79KOA4fYBHgSuNMV969JsJzATYf//9g1wiELdsrO+R068FeZQgaKnd/CN8W6C4spIuwDQ2JnWLoLG2lsa/Pd5jqV/22gvT3EzxXns5Xislmv9JIkgQn1b4i5+0efzGmBONMYc4vBYBn4T28O29fMeoJhEpxTL6DxpjHotxvfnGmAnGmAkDBw5M9u1045a6pyl9il96jx9E1YxR3R5+cVU5VTNGaVR/jmIbVjuvvrOhwYq8D+XY1//kxqR42E7aAgBm927Pa2WDJoEbTkF8AkwZ0/MzXJX54idblvqfAC4IfX8BsCi6g4gI8EdgjTHmN2mcmyduqXua0qcEoff4QVTPOpphcyZTPetoNfo5jJtBtkmWkfWjIeB0rWzRJHDijPFDOfPIoYRvXhjg0Te39FjCd0vlS3aFv3wkWwz/HOAkEVkPnBR6j4gMEZGnQ32OBb4HnCAib4dep2ZmunvQlD5FUcLxY0CTYWT9aghEX8vtvExrEtgsWftpD2E0J5W+RFP8nNIGC4WsMPzGmM+NMTXGmFGhrztC7VuNMaeGvv+HMUaMMeOMMYeHXk97j5x6wiV8BUuu97bR++n+vqIUKH4MaDKMrJO2gJ9rZasmgY3fJfxEKvwVemBgtqTz5TS2hK+iKNlBMvP3gzLoqit7FNMJJ1lG1g7Es4MIpbISmpow7Xvii5yuFX1etmgS2AQJ8ItXma/QAwPV8CuKklf4LbWbqoeDaMOaqqh++1rhY/kVGco2TYLwfPyqvUopLRLaw7Klkq3SV+iBgVmh3JdqUqncpyhKduFWardIijDGMLj3YI4bdhyLPljkWbJXSQ9Oan2lxULvshIam9tTIsxTCCV7tSyvoigFg1tVvi7ThcFQ31TPX9b9JWbJXiU9OC27t3caepeX8OGc0xwV+xKl0LX/1fAripJXRJfUDYLbQ0OqyVYlvXSQiWX3RAID8wHd41cUJa9wKrXrl0QeGuIlm5X00kG8JXdj8ePH/8lDr31MpzEUi3DuxP34+RmHdh8v5JK96vEripJXRJfaLRL/H3PhJXudSIVnns1KeukgFcvuP378n/xp+Ud0hmLYOo3hT8s/4seP/zOhueYLavgVRck7ThtxGs+d9RyrLljFL7/2SyqKY+e7V5ZVegb2RUvxJkt+N5uV9NJBKpbdH3rt40DthYYu9SuKktfYxtxO3assr2RX2y46TEd3n4riCq6feL3nOLE883jz4t1K6yYi8uM3rS9bSPaye6dLtppbe6Ghhl9RlLzntBGnJZzD7+qZhzz/ePfonQR/EhH5ydWYgfBc/qApfNHnFgl0Odj4YpcSxoWG5vEriqL4YP0JNY6euRslQ4YwanGdr77J9NDd5hlkPunGKZe/V2mxryV/p3OLgC6Hvt89Zv+IAL98xiuPXz1+RVEUH8SS4o0myB59MpX0siFmIKj3noiErtO5XUCv0iLaOoxrVH8ho4ZfUZS8JxnyvNFSvBQVQWen+wnGsP6EmrTvr6ciZiAI0R64XQAHcDXiieTyu/Vpae/iwzmqwuiERvUripLX2Nr99U313cp9s1+ZzVMbnwo8VuW0aYxaXMfYNe9Bl9NiciTJivwPQqar73l572645ez7yeVP5NxCRQ2/oih5zby35qVEntevB53unPzKadOo/tnNlAwZAiKUDBlC9c9uTtuqQzzeeyK5/IUuvxsPavgVRclr3GR4E5XndfKs3Uh3Tn74ysSgq65k+9w70iYHHI8Hnkguf6HL78aD7vEripLXDO492LFaX6LyvE517c3u3XQ2NPTom6799Wgykdp37cmjHSP0Y3ngieTyF7L8bjyox68oSl5zxRFX9FDuqyiuiCnP64dwz3rU4jr2/dENGd1fjyYTcsDJ8sAfX7mFY+cs5oBZT3HsnMU8vnJLaiZcgKjHryhKXhOt3BdvVL8fnFYBMqmal6nUvkQ98HgyAxT/qICPoihKlhOvwE8uivkAHDtnsWPFvqFVvVg264QMzCj38BLw0aV+RVEUB1JRiS+eayVSHCjTqX3x4mT0vdqVYKjhVxRFiSJVlfjiuVYi+/SZTu2LFy9Nfd3rTxzd41cUJSUkUnQl03gZ22QbzVjXimefPteq80XjVUXPj4yv4o16/IqiJB07OGtLQzOGPcFZueKtpTMoLta13FIB3drTuVqRKoZ65Pz7kfFVvFHDryhK0olHtjWbCGpsU3mtoPv0mUjhSzbXnjwat8V+leJNHDX8iqIknUSKrmQD6QyKi3WtoPv02VCdL1HOGD+U847Zv4fxVyne5KB7/IqiJJ0hVb0cI7BzxVtLZz6+n2sFKdub6ep8yeLnZxzKhH/rn7NxItmM5vEripJ0ogVYwPLWVEM99UTL9IK1gpAL0fxK8vDK41ePX1GUpGMbd/XW/PP+a9t4ddEGdu1opU//ciZNH8mBE4PXE8g29UAl+1CPX1EUJcO8/9o2ljy4lo62ru62krIippw3Ji7jryiq3KcoipLFvLpoQ4TRB+ho6+LVRRsyNCMln8mKpX4R6Q/8BRgObALONsZ8EdWnAngZKMea9yPGmJ+md6aKoijJZ9eO1kDthUoui0JlE9ni8c8C6owxo4C60PtoWoETjDGHAYcD3xCRY9I3RUVRlNTQp395oPZCJNdFobKJbDH804EFoe8XAGdEdzAWu0JvS0Ov/A9QUBQl75k0fSQlZZEfxyVlRUyaPjJDM8o+cl0UKpvIFsO/rzGmHiD0dZBTJxEpFpG3ge3A88aY19I3RUVRlNRw4MTBTDlvTLeH36d/uQb2RZHrolDZRNr2+EXkBcDpr/hHfscwxnQCh4tIFfA3ETnEGLPa5XozgZkA+++/f/AJK4qipJEDJw5WQ+9BrotCZRNpM/zGmBPdjonIJyJSbYypF5FqLI/ea6wGEXkR+AbgaPiNMfOB+WCl88U9cUVRchoNCMsPrj15tKMolEr4BidblvqfAC4IfX8BsCi6g4gMDHn6iEgv4ERgbbomqChK7qEBYfnDGeOHcsuMQxla1QvBquCnSpDxkRUCPiKyD/AwsD/wEfDvxpgdIjIEuMcYc6qIjMMK/CvGemB52Bhzs5/xVcBHUQqTY+csdlweHlrVi2WzTsjAjBQlPWS9ZK8x5nOgxqF9K3Bq6PtVwPg0T01RlBxGA8IUpSfZstSvKIqSdNwCv/I1IKyxtpb1J9SwZuxBrD+hhsba2kxPSclC1PAripK3XHvyaHqVFke0xQoIy1XjaVfl69i6FYyhY+tW6n9yY87MX0kfavgVRckrwg33wf99Pj8Z3u47ICyXjef2uXdElOIFMC0tbJ97R2YmpGQtWRHcl2o0uE9RCoNEa9GvP6HGMvpRlAwZwqjFdUmda7JZM/YgcPo8F2HsmvfSPyElo2h1PkVRCoJEvd6O+vpA7dlESXV1oHalcFHDryhK3pCo4c5l4znoqiuRioqINqmoYNBVV2ZmQkrWooZfUZS8IVHDncvGs3LaNKp/djMlQ4aACCVDhvje4lAKi6zI41cURUkGg6660nGP36/hto3k9rl30FFfT0l1NYOuujJnjGfltGk5M1clc6jhVxQlb0iG4VbjqeQ7avgVRckrog23nd6Xix68oqQCNfyKouQt0el9dl4+oMZfKVjU8CuKkrd4pfdlu+F//7VtvLpoA7t2tNKnfzmTpo/kwImDMz0tJQ9Qw68oSt6Sq3n577+2jSUPrqWjrQuAXTtaWfKgVYVcjb+SKJrOpyhK3pKrefmvLtrQbfRtOtq6eHXRhgzNSMkn1PAripK35Gpe/q4drYHaFSUIutSvKErekqt5+X36lzsa+T79yzMwGyXfUMOvKEpek4t5+ZOmj4zY4wcoKSti0vSRGZyVki+o4VcURcky7AA+jepXUoEafkVRlCyhsbY2YlvijBzYllByDzX8iqIoWYCKDSnpQqP6FUVRsgAvsSFFSSZq+BVFUbIAN1Gh9q1bWXDDMt5/bVuaZ6TkK2r4FUVRsgA3UaGW8r27lfvU+CvJQA2/oihKFuAkNtRZVMrGEd8EVLlPSR4a3KcoipIFhIsNtW/dSkv53mwc8U0+2ffo7j6q3KckAzX8iqIoWYItNrTghmWq3KekDF3qVxRFyTImTR9JSVnkx7Mq9ynJQj1+RVGULEOV+5RUooZfURQlCzlw4mA19EpK0KV+RVGUADTW1rL+hBrWjD2I9SfU0Fhbm+kpKUog1ONXFEXxicrqKvlAVnj8ItJfRJ4XkfWhr3t79C0WkZUi8mQ656goiqKyuko+kBWGH5gF1BljRgF1ofduXAGsScusFEVRwnCT1XVrV5RsJFsM/3RgQej7BcAZTp1EZBhwGnBPeqalKIqyBzdZXbd2RclGssXw72uMqQcIfR3k0u8O4IdAV6wBRWSmiKwQkRWffvpp0iaqKErh4iSrKxUVDLrqysxMSFHiIG3BfSLyAuCUm/Ijn+efDmw3xrwpIsfH6m+MmQ/MB5gwYYLxP1NFURRnwmV1O+rrKamuZtBVV2pgn5JTpM3wG2NOdDsmIp+ISLUxpl5EqoHtDt2OBb4pIqcCFUA/EfmTMea7KZqyoihKD2xZXUXJVbJlqf8J4ILQ9xcAi6I7GGOuN8YMM8YMB84BFqvRVxRFUZRgZIvhnwOcJCLrgZNC7xGRISLydEZnpiiKoih5RFYI+BhjPgdqHNq3Aqc6tL8IvJjyiSmKoihKnpEtHr+iKIqiKGlADb+iKIqiFBBq+BVFURSlgFDDryiKoigFhBp+RVEURSkg1PAriqIoSgGhhl9RFEVRCggxJv9l7EXkU+BfmZ5HEhgAfJbpSaQIvbfcRO8tN9F7y02C3Nu/GWMGOh0oCMOfL4jICmPMhEzPIxXoveUmem+5id5bbpKse9OlfkVRFEUpINTwK4qiKEoBoYY/t5if6QmkEL233ETvLTfRe8tNknJvusevKIqiKAWEevyKoiiKUkCo4c9CROQbIrJORD4QkVkOx0VE7gwdXyUiR2RinvHg497OC93TKhF5RUQOy8Q84yHWvYX1O0pEOkXkrHTOLxH83JuIHC8ib4vIuyLyUrrnGC8+/iYrRaRWRN4J3dtFmZhnUETkXhHZLiKrXY7n8udIrHvL5c8Rz3sL6xf/54gxRl9Z9AKKgQ3ACKAMeAc4KKrPqcDfAQGOAV7L9LyTeG9fBfYOfX9KPt1bWL/FwNPAWZmedxJ/b1XAe8D+ofeDMj3vJN7bDcCtoe8HAjuAskzP3ce9HQccAax2OZ6TnyM+7y0nP0f83FuoT0KfI+rxZx9HAx8YYzYaY9qAhcD0qD7TgQeMxXKgSkSq0z3ROIh5b8aYV4wxX4TeLgeGpXmO8eLn9wbwX8CjwPZ0Ti5B/Nzbd4DHjDEfARhjcuX+/NybAfqKiAB9sAx/R3qnGRxjzMtYc3UjVz9HYt5bDn+O+Pm9QYKfI2r4s4+hwMdh7zeH2oL2yUaCzvv7WB5JLhDz3kRkKPAt4O40zisZ+Pm9HQjsLSIvisibInJ+2maXGH7u7XfAWGAr8E/gCmNMV3qml1Jy9XMkKLn0ORKTZHyOlCRvOkqSEIe26NQLP32yEd/zFpEpWP+wX0vpjJKHn3u7A7jOGNNpOY85g597KwGOBGqAXsCrIrLcGPN+qieXIH7u7WTgbeAEYCTwvIgsNcZ8meK5pZpc/RzxTQ5+jvjhDhL8HFHDn31sBvYLez8My9MI2icb8TVvERkH3AOcYoz5PE1zSxQ/9zYBWBj6Zx0AnCoiHcaYx9Myw/jx+zf5mTGmCWgSkZeBw4BsN/x+7u0iYI6xNlc/EJEPgTHA6+mZYsrI1c8RX+To54gfEv4c0aX+7OMNYJSIHCAiZcA5wBNRfZ4Azg9F5R4DNBpj6tM90TiIeW8isj/wGPC9HPAWw4l5b8aYA4wxw40xw4FHgP/MAaMP/v4mFwGTRaRERPYCJgJr0jzPePBzbx9hrWQgIvsCo4GNaZ1lasjVz5GY5PDnSEyS8TmiHn+WYYzpEJHLgWexIjfvNca8KyKXho7fjRXJeSrwAbAbyyPJenze243APsDvQ0+0HSYHCm74vLecxM+9GWPWiMgzwCqgC7jHGOOZjpQN+Py9/Qy4X0T+ibU8fp0xJuurv4nIQ8DxwAAR2Qz8FCiF3P4cAV/3lpOfI+Dr3hK/Rig1QFEURVGUAkCX+hVFURSlgFDDryiKoigFhBp+RVEURSkg1PAriqIoSgGhhl9RFEVRCgg1/IqiJB0RWS0is5MwziYRuSYJU1IUJYTm8SuKknFCDwlnGWMOiTp0FNCU/hkpSv6ihl9RlKzFGPNppuegKPmGLvUrSg4TqoZ3t4jME5EvQq9fi0hR6PgMEVklIs0iskNEXgrJztrnTwtV02sRkQ9F5Bch6Vr7eI+l9tA1fxf2fpCILApd418icrHDPPcXkb+JyM7Q6zERGRY6diGWOtnBImJCrwudrh869h+h6+0WkfdFZIqIDBORZ0WkSUTeFpEjoq7/1dC97xaRLSLyBxHpl9hPX1FyEzX8ipL7nIf1vzwJ+P+AmcCVIjIYq778AqyysscB/2efJCInAw9ilZ09GLgYOAv4ZcDr3w98BTgROAM4Hxgedh0BHgf2xapwNwUYAjweOvYX4HZgHVAdev3F43o/Dt3XYcAK4CHgj8DvgfFYhWbuD7v+ocBzWNr0hwEzgMOBewPep6LkBbrUryi5Tz3wg1D1uLUiciBwNfAilsb3I8aYf4X6huvn/wj4tTHmvtD7DSJyHfAnEbnW+NDzDl3rFOBrxphlobYLiCxicyKWwR1pjNkU6vMdLI34GmPMCyKyC0tPfZuP+33AGPNQaJxfAucCzxpjFoXafgUsEZEBIU39a4G/GGNuD5v3fwArRWSQMWa7j2sqSt6gHr+i5D7Lo4z0q8BQYAPwArBaRB4NLZEPDOt3JPAjEdllv4A/A72BwT6vPRarKE93idrQQ8bWqD5bbaMf6rMx1Ocgn9cJZ1XY95+Evv7ToW1Q6OuRwHej7nNZ6NjIOK6vKDmNevyKkr8YYCpwTOjr94FbROTrxph3sB78bwL+6nCuHVTXhVWRLpzSsO+jjzkhobm4zTEo7Q7nO7UVhX29B5jrMNaWOK6vKDmNGn5FyX0mioiEef3HYHnYX4bevwq8KiI3A+8C3wbeAd4CxhhjPvAY+1OsPXcARKQCGAOsDDWtwTKsRwGvhPrsj7WHb/MeMFREhoct9Y8I9Xkv1KcNqyxuKngLODjGfSpKwaBL/YqS+wwB7hCR0SJyFtae9lwROUZEfiwiR4WM8TeB/dhjbG8GviMiN4vIISIyRkTOCu2R2ywGzhOR40XkYKyAuG6P3xizDngG+B8RmSQih2MF1jWHjfEC1oPGgyJypIhMwAoqfCs0PsAm4N9E5AgRGSAi5Un76cCtwNGh7IfxIvIVETldRP4niddQlJxBDb+i5D4PYnnLrwH/ixXhPhdoBI4FngTWY0XO/8wY8ycAY8yzwGlYUfavh16zgI/Cxr4FyzgvwoqM/weWwQ7nQuDDUL9arDiBTfbB0ErEGVirBy8CS4BtwBlhqxSPAk8DdaF+58b3o+iJMWYVVkbDcOAlrIeQW9gTC6AoBYX4CNxVFCVLEZEXgdXGmMszPRdFUXID9fgVRVEUpYBQw68oiqIoBYQu9SuKoihKAaEev6IoiqIUEGr4FUVRFKWAUMOvKIqiKAWEGn5FURRFKSDU8CuKoihKAaGGX1EURVEKiP8fvrUQvcaZrTcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "colors = plt.get_cmap(\"tab10\").colors[::-1]\n",
    "labels = df.index.unique()\n",
    "\n",
    "X = gplvm.X_loc.detach().numpy()\n",
    "for i, label in enumerate(labels):\n",
    "    X_i = X[df.index == label]\n",
    "    plt.scatter(X_i[:, 0], X_i[:, 1], c=[colors[i]], label=label)\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"pseudotime\", fontsize=14)\n",
    "plt.ylabel(\"branching\", fontsize=14)\n",
    "plt.title(\"GPLVM on Single-Cell qPCR data\", fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the first dimension of the latent $X$ for each cell (horizontal axis) corresponds well with the observed capture time (colors). On the other hand, the 32 TE cell and 64 TE cell are clustered near each other. And the fact that ICM cells differentiate into PE and EPI can also be observed from the figure!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks\n",
    "\n",
    "+ The sparse version scales well (linearly) with the number of data points. So the GPLVM can be used with large datasets. Indeed in [2] the authors have applied GPLVM to a dataset with 68k peripheral blood mononuclear cells.\n",
    "\n",
    "+ Much of the power of Gaussian Processes lies in the function prior defined by the kernel. We recommend users try out different combinations of kernels for different types of datasets! For example, if the data contains periodicities, it might make sense to use a [Periodic kernel](http://docs.pyro.ai/en/dev/contrib.gp.html#periodic). Other kernels can also be found in the [Pyro GP docs](http://docs.pyro.ai/en/dev/contrib.gp.html#module-pyro.contrib.gp.kernels)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "[1] `Resolution of Cell Fate Decisions Revealed by Single-Cell Gene Expression Analysis from Zygote to Blastocyst`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Guoji Guo, Mikael Huss, Guo Qing Tong, Chaoyang Wang, Li Li Sun, Neil D. Clarke, Paul Robson\n",
    "\n",
    "[2] `GrandPrix: Scaling up the Bayesian GPLVM for single-cell data`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Sumon Ahmed, Magnus Rattray, Alexis Boukouvalas\n",
    "\n",
    "[3] `Bayesian Gaussian Process Latent Variable Model`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Michalis K. Titsias, Neil D. Lawrence\n",
    "\n",
    "[4] `A novel approach for resolving differences in single-cell gene expression patterns from zygote to blastocyst`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Florian Buettner, Fabian J. Theis"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}