Tutorial_02_jump_distance_mixture_modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "410dff76",
   "metadata": {},
   "source": [
    "# Ermine Tutorial 02\n",
    "## Jump-distance mixture modeling of single-particle tracking data with ermine\n",
    "\n",
    "by Sebastian Malkusch (c) 2021  \n",
    "malkusch@med.uni-frankfurt.de  \n",
    "Data Science| Klinische Pharmakologie  \n",
    "Institut für Klinische Pharmakologie  \n",
    "Goethe - Universität  \n",
    "Theodor-Stern-Kai 7  \n",
    "60590 Frankfurt am Main"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aae860f",
   "metadata": {},
   "source": [
    "## Abstract\n",
    "In this tutorial, we will use unsupervised machine learning to build a multimodal model that describes the probability density distribution of jump widths from Tutorial 01:  \n",
    "$p(r) = \\sum_{i=1}^{n_{max}} \\omega_{i} \\frac{2 r}{\\mu_{i}} e^{-\\frac{r^2}{\\mu_{i}}}$   \n",
    "with   \n",
    "$\\mu_{i} = 4 D_{i} \\tau$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c24d0e11",
   "metadata": {},
   "source": [
    "## Prerequisite\n",
    "* Tutorial 01 must be completed\n",
    "* The file Fab_CS6_cell20_jumps.csv needs to be saved in folder results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df7a2814",
   "metadata": {},
   "source": [
    "## Import libraries\n",
    "You need to install Ermine prior to usage.  \n",
    "Run  \n",
    "python -m pip install pyErmine  \n",
    "from the command line  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d94d5818",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import ermine as em\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e71984a1",
   "metadata": {},
   "source": [
    "## Import Jump_df\n",
    "Store the path to the jump data set in the variable \"file_name\"\n",
    "The routine imports the jump data to a data frame (data_df)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ee9ed296",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>id_departure</th>\n",
       "      <th>frame_departure</th>\n",
       "      <th>x [nm]_departure</th>\n",
       "      <th>y [nm]_departure</th>\n",
       "      <th>sigma [nm]_departure</th>\n",
       "      <th>intensity [photon]_departure</th>\n",
       "      <th>offset [photon]_departure</th>\n",
       "      <th>bkgstd [photon]_departure</th>\n",
       "      <th>uncertainty_xy [nm]_departure</th>\n",
       "      <th>...</th>\n",
       "      <th>seg.id_destination</th>\n",
       "      <th>seg.loc_count_destination</th>\n",
       "      <th>seg.track_pos_destination</th>\n",
       "      <th>track.noise_destination</th>\n",
       "      <th>track.t_start_destination</th>\n",
       "      <th>track.lifetime_destination</th>\n",
       "      <th>track.id_destination</th>\n",
       "      <th>track.loc_count_destination</th>\n",
       "      <th>track.seg_count_destination</th>\n",
       "      <th>jump_distance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>149</td>\n",
       "      <td>1</td>\n",
       "      <td>31496.6</td>\n",
       "      <td>32117.6</td>\n",
       "      <td>167.4</td>\n",
       "      <td>665.8</td>\n",
       "      <td>34.8</td>\n",
       "      <td>9.3</td>\n",
       "      <td>22.5</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>13.285330</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>357</td>\n",
       "      <td>2</td>\n",
       "      <td>31489.3</td>\n",
       "      <td>32106.5</td>\n",
       "      <td>168.6</td>\n",
       "      <td>916.5</td>\n",
       "      <td>33.0</td>\n",
       "      <td>10.8</td>\n",
       "      <td>19.2</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>7.964923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>570</td>\n",
       "      <td>3</td>\n",
       "      <td>31496.3</td>\n",
       "      <td>32110.3</td>\n",
       "      <td>151.3</td>\n",
       "      <td>860.2</td>\n",
       "      <td>38.8</td>\n",
       "      <td>8.0</td>\n",
       "      <td>14.1</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>50.420631</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>783</td>\n",
       "      <td>4</td>\n",
       "      <td>31470.3</td>\n",
       "      <td>32067.1</td>\n",
       "      <td>182.5</td>\n",
       "      <td>594.9</td>\n",
       "      <td>32.0</td>\n",
       "      <td>7.3</td>\n",
       "      <td>24.0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>26.654268</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>996</td>\n",
       "      <td>5</td>\n",
       "      <td>31472.0</td>\n",
       "      <td>32093.7</td>\n",
       "      <td>182.9</td>\n",
       "      <td>630.6</td>\n",
       "      <td>31.6</td>\n",
       "      <td>9.5</td>\n",
       "      <td>27.6</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>57.015437</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 74 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  id_departure  frame_departure  x [nm]_departure  \\\n",
       "0           0           149                1           31496.6   \n",
       "1           1           357                2           31489.3   \n",
       "2           2           570                3           31496.3   \n",
       "3           3           783                4           31470.3   \n",
       "4           4           996                5           31472.0   \n",
       "\n",
       "   y [nm]_departure  sigma [nm]_departure  intensity [photon]_departure  \\\n",
       "0           32117.6                 167.4                         665.8   \n",
       "1           32106.5                 168.6                         916.5   \n",
       "2           32110.3                 151.3                         860.2   \n",
       "3           32067.1                 182.5                         594.9   \n",
       "4           32093.7                 182.9                         630.6   \n",
       "\n",
       "   offset [photon]_departure  bkgstd [photon]_departure  \\\n",
       "0                       34.8                        9.3   \n",
       "1                       33.0                       10.8   \n",
       "2                       38.8                        8.0   \n",
       "3                       32.0                        7.3   \n",
       "4                       31.6                        9.5   \n",
       "\n",
       "   uncertainty_xy [nm]_departure  ...  seg.id_destination  \\\n",
       "0                           22.5  ...                 1.0   \n",
       "1                           19.2  ...                 1.0   \n",
       "2                           14.1  ...                 1.0   \n",
       "3                           24.0  ...                 1.0   \n",
       "4                           27.6  ...                 1.0   \n",
       "\n",
       "   seg.loc_count_destination  seg.track_pos_destination  \\\n",
       "0                        9.0                        0.0   \n",
       "1                        9.0                        0.0   \n",
       "2                        9.0                        0.0   \n",
       "3                        9.0                        0.0   \n",
       "4                        9.0                        0.0   \n",
       "\n",
       "   track.noise_destination  track.t_start_destination  \\\n",
       "0                      0.0                        1.0   \n",
       "1                      0.0                        1.0   \n",
       "2                      0.0                        1.0   \n",
       "3                      0.0                        1.0   \n",
       "4                      0.0                        1.0   \n",
       "\n",
       "   track.lifetime_destination  track.id_destination  \\\n",
       "0                         8.0                   1.0   \n",
       "1                         8.0                   1.0   \n",
       "2                         8.0                   1.0   \n",
       "3                         8.0                   1.0   \n",
       "4                         8.0                   1.0   \n",
       "\n",
       "   track.loc_count_destination  track.seg_count_destination  jump_distance  \n",
       "0                          9.0                          1.0      13.285330  \n",
       "1                          9.0                          1.0       7.964923  \n",
       "2                          9.0                          1.0      50.420631  \n",
       "3                          9.0                          1.0      26.654268  \n",
       "4                          9.0                          1.0      57.015437  \n",
       "\n",
       "[5 rows x 74 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "file_name = \"results/Fab_CS6_cell20_jumps.csv\"\n",
    "data_df = pd.read_csv(filepath_or_buffer = file_name)\n",
    "data_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dcbfa3d",
   "metadata": {},
   "source": [
    "## Screen for optimal model parameters\n",
    "The routine trains multiple models and evaluates the model quality by means of log Likelihood, BIC, AIC and AICc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b2312b20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>classes</th>\n",
       "      <th>dof</th>\n",
       "      <th>instances</th>\n",
       "      <th>log_likelihood</th>\n",
       "      <th>BIC</th>\n",
       "      <th>AIC</th>\n",
       "      <th>AICc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>130621</td>\n",
       "      <td>-695942.642841</td>\n",
       "      <td>1.391897e+06</td>\n",
       "      <td>1.391887e+06</td>\n",
       "      <td>1.391887e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>130621</td>\n",
       "      <td>-687172.149107</td>\n",
       "      <td>1.374380e+06</td>\n",
       "      <td>1.374350e+06</td>\n",
       "      <td>1.374350e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>130621</td>\n",
       "      <td>-686949.356781</td>\n",
       "      <td>1.373958e+06</td>\n",
       "      <td>1.373909e+06</td>\n",
       "      <td>1.373909e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>7</td>\n",
       "      <td>130621</td>\n",
       "      <td>-686949.305839</td>\n",
       "      <td>1.373981e+06</td>\n",
       "      <td>1.373913e+06</td>\n",
       "      <td>1.373913e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5</td>\n",
       "      <td>9</td>\n",
       "      <td>130621</td>\n",
       "      <td>-686952.643902</td>\n",
       "      <td>1.374011e+06</td>\n",
       "      <td>1.373923e+06</td>\n",
       "      <td>1.373923e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   classes  dof  instances  log_likelihood           BIC           AIC  \\\n",
       "0        1    1     130621  -695942.642841  1.391897e+06  1.391887e+06   \n",
       "0        2    3     130621  -687172.149107  1.374380e+06  1.374350e+06   \n",
       "0        3    5     130621  -686949.356781  1.373958e+06  1.373909e+06   \n",
       "0        4    7     130621  -686949.305839  1.373981e+06  1.373913e+06   \n",
       "0        5    9     130621  -686952.643902  1.374011e+06  1.373923e+06   \n",
       "\n",
       "           AICc  \n",
       "0  1.391887e+06  \n",
       "0  1.374350e+06  \n",
       "0  1.373909e+06  \n",
       "0  1.373913e+06  \n",
       "0  1.373923e+06  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed  = 42\n",
    "np.random.seed(seed)\n",
    "x = data_df[\"jump_distance\"].values\n",
    "\n",
    "jmm = em.JumpDistanceMixtureModel(n_components=1,\n",
    "                                  degrees_of_freedom = 4,\n",
    "                                  tau = 0.02,\n",
    "                                  init_params = \"wm\",\n",
    "                                  params = \"wm\")\n",
    "jmm.fit(x)\n",
    "model_df = pd.DataFrame(jmm.evaluate(x))\n",
    "for i in np.arange(start=2, stop=6, step=1):\n",
    "    jmm = em.JumpDistanceMixtureModel(n_components=i,\n",
    "                                      degrees_of_freedom = 4,\n",
    "                                      tau = 0.02,\n",
    "                                      init_params = \"wm\",\n",
    "                                      params = \"wm\")\n",
    "    jmm.fit(x)\n",
    "    model_df = model_df.append(pd.DataFrame(jmm.evaluate(x)))\n",
    "        \n",
    "model_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eb45e29",
   "metadata": {},
   "source": [
    "## Identify optimal model by BIC\n",
    "Plots the BIC and idententifies the optimal model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ee286dcb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwkklEQVR4nO3de3xV9Z3v/9c7N0LCVRIQSCKoKApixCSVVqlTtbWOFfFGeLSjrW2tZ9rptD09p3WmM+3U6fn1Nmc608s41nq000rQItaq9TLaqm1VCIgYBCsKQiASCHK/5fL5/bG+GzbJzo3sZOXyeT4e+5G1v991+axvdvYn6/tdF5kZzjnnXLKMuANwzjnX/3hycM4514YnB+ecc214cnDOOdeGJwfnnHNteHJwzjnXhieHIU7SNyT9Iu44UpG0T9KpvbDej0v6Q7rX28Vtb5R0aRfmmyLJJGX1ZD2DhaTfS/pU3HEMJZ4cBrnwBZt4tUg6mPT+o3HHl5Dqj9/MRpjZW3HFNNRJqpT0uqTdkuol3StpVNxxub7hyWGQC1+wI8xsBLAJ+EhS2S/jjs/1a38E3mdmo4FTgSzgn+MNyfUVTw4OIEfSzyXtlbRGUlmiQtIkSUskbZe0QdLn21uJpNFhPdslvS3pa5IyQt3HJf1R0g/Df6LrJF0S6r4FXAT8KBzR/CiUm6TTw/Q9kn4i6bdhnj9KOlnSDyS9G9Z3XlIsX5X0Ztin1yTN70pDJHXnfELS5rDuWyWVS1otaVcivjB/RtjPt8N/1z+XNDqp/q9CXYOkv2+1rYykOBsk3S/ppK7EGZSGmHZLWiwpN6x3rKRHwu/h3TBdlLTdj0t6K7TNhvaOIM1ss5ntSCpqBk7voO1M0l9LeiOs+3ZJp0l6QdKesH85SfN/WtJ6STslPSxpUlLdZeF3uju0t1pt62ZJa8P+PSHplFAuSf8afhe7Q/vM7EabugQz89cQeQEbgUtblX0DOARcAWQC/x/wYqjLAFYA/wjkEP33+BbwoXbW/3Pg18BIYArwZ+CToe7jQBPwRSAbWADsBk4K9b8HPtVqfQacHqbvAXYA5wO5wDPABuDGEPc/A79LWvZ6YFLYhwXAfmBiUix/aGcfpoTt3hG288HQPg8B44HJQD3w/jD/zcD60DYjgAeB/wp1ZwP7gLnAMOD/hja4NNR/AXgRKAr1/wksahVHVge/y2VhH08C1gK3hrpxwLVAXvhdPAA8FOrygT3AmeH9RGBGB5+ZC8PvyUIbfrCDeQ14GBgFzAAOA0+HthkNvAbcFOb9QPh9zg77/kPguVBXEGK8LnxWvhja7VOh/urQ5mcRHc18DfhTqPsQ0Wd2DFFCOSvxe/dXN78v4g4gbTsCd4c/2pouzn9D+LCuAe6LO/4+aqONpE4O/530/mzgYJh+D7Cp1fy3Af8vxbozw5fB2UllnwF+H6Y/DmwFlFS/DPirMP17Ok8OP02q+xtgbdL7c4BdHez7KmBeUiydJYfJSWUNwIKk90uAL4Tpp4G/Tqo7E2gMX1r/CFQl1eUDRziWHNYClyTVT0xaNhFHR8nhY0nvvwvc0c68pcC7STHsIkoew7vx2ZkcPitndDCPEXVDJd6vAL6S9P5fgB+E6Z8B302qGxH2fQpRwn8xqU5ALceSw28J/3SE9xnAAeAUoqTzZ+ACICPuv7mB/BpM3Ur3AJd3ZUZJ04i+5N5nZjOI/oMbyt5Jmj4A5Co6S+YUYFLoStklaRfwd8CEFOsoIDq6eDup7G2iL5WELRb+mpPqJ9F125KmD6Z4PyLxRtKNklYlxT0zxJjubU2i7T5nEbXRJGBzosLM9hMlmoRTgKVJMa4l6rpJ1b6ptP69jQCQlCfpP0N31h7gOWCMpMwQwwLgVqBO0qOSpne2ITPbAjwOVHUy6wm1m5ntI2qbybRtN0t+T9Ru/5bUbjuJEshkM3sG+BHwY2CbpDvlg+gnZNAkBzN7juhDclTo73xc0gpJzyf9EXwa+LGZvRuWre/jcAeKzcAGMxuT9BppZlekmHcH0X9+pySVlQBbkt5PlqRW9VvDdNpuDxz6n38KfA4YZ2ZjgBpa9VunyVba7nMT0ZdiHVCcFFceUZdPwmbgw63aNzd8EffE/yQ6gnmPmY0i6taCsP9m9oSZXUZ0pLKOqK26Igs4rYexJRzXbpLyidpmC23bTcnvidrtM63abbiZ/QnAzP7dzM4n6to6A/hfaYp5SBk0yaEddwJ/Ez4oXwZ+EsrPAM4Ig5ovSurSEccQtAzYI+krkoZLypQ0U1J56xnNrBm4H/iWpJHhC/pLQPI1FOOBz0vKlnQ9UX/wY6FuG1HfdDrkEyWb7QCSPkF05NAbFgFflDRV0gjg/wCLzawJ+BVwpaQLw0DsNzn+b+4OovZKDKYWSpqXhphGEv2XvisMcH89USFpgqSrwpfxYaIxkeZUK5H0UUklYZD3FOBbRN1o6XAf8AlJpZKGEbXbS2a2EXgUmCHpmnAE+3ng5KRl7wBukzQjxDk6fJ5QdOLAeyRlE42RHGpv/1zHBm1yCH+o7wUekLSKaLBvYqjOAqYBFwMLgbskjen7KPu38IX/EaI+6w1ERwd3EQ0upvI3RH+QbwF/IPoCuDup/iWidt9B9EVznZkluln+DbgunH3y7z2M+zWi/u0XiJLOOUSnZfaGu4H/Iuq62UD0ZfQ3IY41wGeJ2qEOeJeo7zzh34gGcJ+UtJdocPo9aYjpB8BwonZ+kag7KCGD6MhiK9GR9vuBv25nPWcDfyJKIH8EXic66u4xM3sa+Aei8Zs6oiOSylC3g+iEgm8TdTVNI+n3Z2ZLge8AVaHbrAb4cKgeRXQk9C5Rt1UD8P10xDzU6Pgu4IFN0hTgETObGfoZXzeziSnmu4NowOue8P5p4Ktmtrwv4x1KJH2caEDxwrhjcc51btAeOZjZHmBD0uGmJJ0bqh8C/iKUFxB1M/mVuM45Fwya5CBpEVE3wpmSaiV9Evgo8ElJrxCdsproz30CaJD0GvA74H8ldW8459yQN6i6lZxzzqXHoDlycM45lz4pbwecTNLdwJVAvZm1OR0wnHp3O9BCdH73F8zsD6Hub4nObhDR1a0/COXfIzoL5gjwJvAJM9sVBpTXEp0VAdGg8a2dxVhQUGBTpkzpbDbnnHNJVqxYscPMClPVddqtJGku0alsP28nOYwA9puZSZoF3G9m08PNrqqACqIk8DjwP8zsDUkfBJ4xsyZJ3wEws68kn23UnR0sKyuz6urq7izinHNDnqQVZlaWqq7TbqVUVx63qt+XdEuExMVHEF3g9KKZHQgXBD0LzA/LPBnK4NiNx5xzzvUTaRlzkDRf0jqiKxtvDsU1wFxJ48JtA67g+EvgE24mupFWwlRJL0t6VtJF6YjPOedc93Q65tAV4YrFpaEL6naiu06uDV1GTxF1S71CNCZxlKL72zcBiYfO1AElZtYg6XzgIUkzwjULtFr2FuAWgJKSknTshnPOuSAtySHBzJ4LN7srMLMdZvYzolvzIun/kHTrAEk3EQ10X5LoljKzw0T3e8HMVkh6k+gCtTYDCmZ2J9G9kygrK/PzcZ1zadXY2EhtbS2HDh2KO5Qey83NpaioiOzs7C4v0+PkoOhJXW+GAenZRLdtbgh1482sXlIJcA0wJ5RfDnyF6IEpB5LWVQjsNLNmRQ+Wn4Zfueyci0FtbS0jR45kypQpHH8z4YHFzGhoaKC2tpapU6d2ebmunMq6iOgGdQWSaonu8JgdNnoH0UNDbpTUSHQnyAVJA9RLJI0jupXzZxO3yCa63/ow4KnQ6IlTVucC35TURHQnxVvNrN3BcOec6y2HDh0a8IkBQBLjxo1j+/bt3Vqu0+RgZgs7qf8O0R0SU9WlHFA2s5TPoTWzJUR3aXTOudgN9MSQcCL7MaSvkK7bfZBvPfoaDfsOxx2Kc871K0M6Oew71MRPn9/Agyt7+uAt55xLv8zMTEpLSzn33HOZPXs2f/rTnwDYuHEjM2ceu1Z42bJlzJ07lzPPPJPp06fzqU99igMHDrS32i5J69lKA820CSM5/5SxLFq+iU9dNHXQHEI65waH4cOHs2rVKgCeeOIJbrvtNp599tnj5tm2bRvXX389VVVVzJkzBzNjyZIl7N27l7y8vBPe9pA+cgBYUF7MW9v3s3zju53P7JxzMdmzZw9jx45tU/7jH/+Ym266iTlz5gDR+MJ1113HhAkTerS9IX3kAHDlrInc/pvXqFq2iYqpJ8UdjnOuH/qn36zhta1trsXtkbMnjeLrH5nR4TwHDx6ktLSUQ4cOUVdXxzPPPNNmnpqaGm666aa0xgZ+5EBeThZXlU7i0Vfr2H2gMe5wnHPuqES30rp163j88ce58cYb6atn8Az5IweAhRUl/PKlTTy0ags3vXdK3OE45/qZzv7D7wtz5sxhx44dba5XmDFjBitWrGDevHntLHlihvyRA8DMyaOZOXkUi5Zt6rOs7Jxz3bFu3Tqam5sZN27cceWf+9znuPfee3nppZeOlv3iF7/gnXfe6dH2/MghqCwv4WsP1fBK7W5Ki8fEHY5zzh0dc4DoNhj33nsvmZmZx80zYcIEqqqq+PKXv0x9fT0ZGRnMnTuXa665pkfb9uQQzCudxLceXcvi5Zs8OTjn+oXm5uaU5VOmTKGmpubo+zlz5vD888+nddverRSMzM3mylkTeXjVVvYfbup8AeecG8Q8OSSprChh/5FmfvPK1rhDcc65WHlySDK7ZAxnTBjBouWb4w7FOdcPDJYTVE5kPzw5JJFEZXkJr2zexdq69F7w4pwbWHJzc2loaBjwCSLxPIfc3NxuLecD0q1cM3sy3358HVXLNvFP82Z2voBzblAqKiqitra2289B6I8ST4LrDk8OrYzJy+HyGSez9OUt3HbFWeRmZ3a+kHNu0MnOzu7Wk9MGG+9WSqGyopg9h5p47NW6uENxzrlYeHJIYc6p45gyLo+qZT4w7Zwbmjw5pCCJBeUlLNu4k/X1++IOxznn+lynyUHS3ZLqJdW0Uz9P0mpJqyRVS7owqe5vJdVIWiPpC0nlJ0l6StIb4efYpLrbJK2X9LqkD/Vw/07YdecXkZUhFi/fFFcIzjkXm64cOdwDXN5B/dPAuWZWCtwM3AUgaSbwaaACOBe4UtK0sMxXgafNbFpY/qthmbOBSmBG2OZPJMUyIlw4chiXnjWBJSu3cLgp9SXszjk3WHWaHMzsOWBnB/X77NiJwPlAYvos4EUzO2BmTcCzwPxQNw+4N0zfC1ydVF5lZofNbAOwnii5xKKyopid+4/w36/VxxWCc87FIi1jDpLmS1oHPEp09ABQA8yVNE5SHnAFUBzqJphZHUD4OT6UTwaSR4FrQ1mqbd4SurGqe+s85IumFTJ5zHCqvGvJOTfEpCU5mNlSM5tOdARweyhbC3wHeAp4HHgF6OyOdkq1+na2eaeZlZlZWWFh4YmG3qHMDHFDWTHPv7GDzTsP9Mo2nHOuP0rr2UqhC+o0SQXh/c/MbLaZzSXqmnojzLpN0kSA8DPRb1PLsaMLgCIg1rvg3VBeRIZgsd9vyTk3hPQ4OUg6XZLC9GwgB2gI78eHnyXANcCisNjDQOKJ2DcBv04qr5Q0TNJUYBqwrKcx9sTE0cO5+MzxPLBiM03NLXGG4pxzfabT22dIWgRcDBRIqgW+DmQDmNkdwLXAjZIagYPAgqQB6iWSxgGNwGfN7N1Q/m3gfkmfBDYB14f1rZF0P/AaURfUZ80s9lOFKsuLueW/6vnd69u57OwJcYfjnHO9TgP9joMAZWVlVl1d3Wvrb2xu4b3ffoZZk0fzs4+X99p2nHOuL0laYWZlqer8CukuyM7M4Przi/jd6/XU7T4YdzjOOdfrPDl00YLyYloMHqiujTsU55zrdZ4cuuiUcfm87/RxLF6+mZaWgd8V55xzHfHk0A2V5SVs2XWQ59fviDsU55zrVZ4cuuGDMyYwNi+bqmV+xbRzbnDz5NANw7IyuXZ2EU+9to0d+w7HHY5zzvUaTw7dVFlRTFOLsWSFD0w75wYvTw7ddPr4kZRPGcvi5ZsZDNeIOOdcKp4cTkBleQlv7djPSxvavZO5c84NaJ4cTsAV50xkZG6WD0w75wYtTw4nYHhOJvPPm8xjNe+w68CRuMNxzrm08+RwghaUF3OkqYWlL2+JOxTnnEs7Tw4naMak0cwqGk3VMh+Yds4NPp4ceqCyvITXt+3l5c274g7FOefSypNDD1xVOom8nEwfmHbODTqeHHpgxLAsPjJrEr95pY69hxrjDsc559LGk0MPVVYUc7CxmYdfifVR1845l1aeHHqotHgM008eyeLlm+MOxTnn0saTQw9JorK8mNW1u1mzdXfc4TjnXFp0mhwk3S2pXlJNO/XzJK2WtEpStaQLk+q+KGmNpBpJiyTlhvLFYf5VkjZKWhXKp0g6mFR3R5r2s1fNP6+IYVkZVC3zowfn3ODQlSOHe4DLO6h/GjjXzEqBm4G7ACRNBj4PlJnZTCATqAQwswVmVhqWWQI8mLS+NxN1ZnZr93YnHqPzsrninIk8tGoLB480xx2Oc871WKfJwcyeA9q9w5yZ7bNjV4HlA8lXhGUBwyVlAXnAcaO2kgTcACzqZtz9TmV5MXsPNfHoq3Vxh+Kccz2WljEHSfMlrQMeJTp6wMy2AN8HNgF1wG4ze7LVohcB28zsjaSyqZJelvSspIs62OYtoRurevv27enYjR6pmHoSpxbm+zUPzrlBIS3JwcyWmtl04GrgdgBJY4F5wFRgEpAv6WOtFl3I8UcNdUCJmZ0HfAm4T9KodrZ5p5mVmVlZYWFhOnajRxID09Vvv8sb2/bGHY5zzvVIWs9WCl1Qp0kqAC4FNpjZdjNrJBpXeG9i3tDVdA2wOGn5w2bWEKZXAG8CZ6Qzxt50zewisjNFlZ/W6pwb4HqcHCSdHsYOkDQbyAEaiLqTLpCUF+ovAdYmLXopsM7MapPWVSgpM0yfCkwD3uppjH2lYMQwLjt7Ag+urOVwkw9MO+cGrq6cyroIeAE4U1KtpE9KulVS4kyia4GacDrqj4EFFnkJ+BWwEng1bOvOpFVX0nYgei6wWtIrYdlbzWxAPW6tsryEdw808sSabXGH4pxzJ0yD4XbTZWVlVl1dHXcYALS0GHO/9ztKTsrjvk9fEHc4zjnXLkkrzKwsVZ1fIZ1mGRliQVkxf3qzgbcb9scdjnPOnRBPDr3g+rJiMoTfb8k5N2B5cugFJ4/O5QPTx/PAiloam1viDsc557rNk0MvqSwvYfvewzyzrj7uUJxzrts8OfSSi88sZMKoYX7FtHNuQPLk0EuyMjO4oayYZ/+8na27DsYdjnPOdYsnh150Q1kxBtxf7QPTzrmBxZNDLyo+KY8LTy/g/uWbaW4Z+NeTOOeGDk8OvayyvIStuw/x3Bvx3znWOee6ypNDL7vs7AmMy8/xgWnn3IDiyaGX5WRlcO35RTy9tp76vYfiDsc557rEk0MfWFBeTFOL8asVtZ3P7Jxz/YAnhz5wWuEIKqaexOLlm2nxgWnn3ADgyaGPLKwo5u2GA7z4VkPcoTjnXKc8OfSRD8+cyKjcLH9KnHNuQPDk0EdyszO5ZnYRj9e8w7v7j8QdjnPOdciTQx+qrCjmSHMLD768Je5QnHOuQ54c+tD0k0dRWjyGqmWbGAxP4HPODV5deYb03ZLqJdW0Uz9P0mpJqyRVS7owqe6LktZIqpG0SFJuKP+GpC1hmVWSrkha5jZJ6yW9LulD6djJ/mRhRTFv1O9j5aZ34w7FOefa1ZUjh3uAyzuofxo418xKgZuBuwAkTQY+D5SZ2UwgE6hMWu5fzaw0vB4Ly5wd5pkRtvkTSZnd2qN+7spZk8jPyWTRMh+Yds71X50mBzN7DtjZQf0+O9ZHkg8k95dkAcMlZQF5wNZONjcPqDKzw2a2AVgPVHQW40CSPyyLq0on8cjqrew51Bh3OM45l1JaxhwkzZe0DniU6OgBM9sCfB/YBNQBu83syaTFPhe6o+6WNDaUTQaS/6WuDWWptnlL6Maq3r59YN3UrrK8hEONLfx6VWe50jnn4pGW5GBmS81sOnA1cDtA+MKfB0wFJgH5kj4WFvkP4DSglChx/EsoV6rVt7PNO82szMzKCgsL07EbfWZW0WjOmjjKb8bnnOu30nq2UuiCOk1SAXApsMHMtptZI/Ag8N4w3zYzazazFuCnHOs6qgWKk1ZZROddUQOOJBZWFLNm6x5erd0ddzjOOddGj5ODpNMlKUzPBnKABqLupAsk5YX6S4C1Yb6JSauYDyTOhHoYqJQ0TNJUYBqwrKcx9kfzSieTm53BouV+9OCc63+yOptB0iLgYqBAUi3wdSAbwMzuAK4FbpTUCBwEFoQB6pck/QpYCTQBLwN3htV+V1IpUZfRRuAzYX1rJN0PvBaW+ayZNadlT/uZ0cOzueKciTy8ait/f8VZ5A/r9FfhnHN9RoPhYqyysjKrrq6OO4xuW75xJ9ff8QLfvW4WN5QVd76Ac86lkaQVZlaWqs6vkI5R2SljOX38CB+Yds71O54cYiSJyvJiVm7axZ+37Y07HOecO8qTQ8yumV1ETmYGi/zowTnXj3hyiNlJ+Tl8cMYElr68hUONg3Ls3Tk3AHly6AcWVpSw60AjT6x5J+5QnHMO8OTQL8w5dRzFJw33riXnXL/hyaEfyMgQleUlvPjWTjbs2B93OM4558mhv7j+/CIyM0SVXzHtnOsHPDn0E+NH5fKB6eNZsqKWI00tcYfjnBviPDn0Iwsritmx7whPr90WdyjOuSHOk0M/8v4zxjNxdC6LlvtT4pxz8fLk0I9kZojry4p5/o3t1L57IO5wnHNDmCeHfuaGsiIA7q+ujTkS59xQ5smhnykam8fcaYU8UL2Z5paBf8dc59zA5MmhH1pYUUzd7kM8++f6uENxzg1Rnhz6oUvOmkDBiGEsWuYD0865eHhy6IeyMzO47vwinllXT/2eQ3GH45wbgjw59FMLyotpbjEeWOED0865vufJoZ+aWpDPBaeeRNXyTbT4wLRzro91mhwk3S2pXlJNO/XzJK2WtEpStaQLk+q+KGmNpBpJiyTlhvLvSVoXllsqaUwonyLpYFjXKkl3pGk/B6SFFSVs3nmQP73ZEHcozrkhpitHDvcAl3dQ/zRwrpmVAjcDdwFImgx8Higzs5lAJlAZlnkKmGlms4A/A7clre9NMysNr1u7sS+DzodmnMyYvGwW+c34nHN9rNPkYGbPATs7qN9nZol+j3wguQ8kCxguKQvIA7aGZZ40s6Ywz4tA0QnEPujlZmcy/7zJPLnmHRr2HY47HOfcEJKWMQdJ8yWtAx4lOnrAzLYA3wc2AXXAbjN7MsXiNwO/TXo/VdLLkp6VdFEH27wldGNVb9++PR270S8trCihsdl4cOWWuENxzg0haUkOZrbUzKYDVwO3A0gaC8wDpgKTgHxJH0teTtLfA03AL0NRHVBiZucBXwLukzSqnW3eaWZlZlZWWFiYjt3ol86YMJLZJWOoWr6JYwdozjnXu9J6tlLogjpNUgFwKbDBzLabWSPwIPDexLySbgKuBD6a6JYys8Nm1hCmVwBvAmekM8aBqLKihDe376f67XfjDsU5N0T0ODlIOl2SwvRsIAdoIOpOukBSXqi/BFgb5rsc+ApwlZkdSFpXoaTMMH0qMA14q6cxDnRXzprIyGFZ/oxp51yf6cqprIuAF4AzJdVK+qSkWyUlziS6FqiRtAr4MbDAIi8BvwJWAq+Gbd0ZlvkRMBJ4qtUpq3OB1ZJeCcveambtDoYPFXk5WVxVOonHXq1j98HGuMNxzg0BGgz92GVlZVZdXR13GL2qZsturvzhH/jmvBncOGdK3OE45wYBSSvMrCxVnV8hPUDMnDyamZNHsWjZZh+Yds71Ok8OA8iC8hLW1u1hde3uuENxzg1ynhwGkHmlkxienUmVXzHtnOtlnhwGkFG52fzlrIk8vGor+w83db6Ac86dIE8OA8zCimL2H2nmN69sjTsU59wg5slhgJldMpZp40ewaLk/Jc4513s8OQwwkqisKOGVzbtYW7cn7nCcc4OUJ4cB6JrzJpOTmcFiP3pwzvUSTw4D0Nj8HC6feTIPrqzlUGNz3OE45wYhTw4DVGVFMXsONfHbmrq4Q3HODUKeHAaoOaeOY8q4PBYt864l51z6eXIYoCSxoLyEZRt28ub2fXGH45wbZDw5DGDXnj+ZrAz5wLRzLu08OQxg40fmcslZ41myopYjTS1xh+OcG0Q8OQxwlRUlNOw/wlOvbYs7FOfcIOLJYYCbO62QyWOG+834nHNp5clhgMvMENeXFfH8GzvYvPNA5ws451wXeHIYBG4oKyZD+MC0cy5tuvIM6bsl1Uuqaad+nqTV4VnQ1ZIuTKr7oqQ1kmokLZKUG8pPkvSUpDfCz7FJy9wmab2k1yV9KB07OdhNGjOc959RyAMrNtPU7APTzrme68qRwz3A5R3UPw2ca2alwM3AXQCSJgOfB8rMbCaQCVSGZb4KPG1m08LyXw3LnB3mmRG2+RNJmd3bpaGpsqKEbXsO8/vXt8cdinNuEOg0OZjZc8DODur32bGHGucDyQ84zgKGS8oC8oDEQwjmAfeG6XuBq5PKq8zssJltANYDFV3blaHtA9PHUzhymA9MO+fSIi1jDpLmS1oHPEp09ICZbQG+D2wC6oDdZvZkWGSCmdWF+eqA8aF8MpDccV4bylJt85bQjVW9fbv/t5ydmcH15xfxzLp63tl9KO5wnHMDXFqSg5ktNbPpREcAtwOEcYR5wFRgEpAv6WOdrEqpVt/ONu80szIzKyssLDzh2AeTBeXFtBg8UO0D0865nknr2UqhC+o0SQXApcAGM9tuZo3Ag8B7w6zbJE0ECD/rQ3ktUJy0yiKOdUW5TpwyLp/3nT6OxdWbaWlJmVOdc65LepwcJJ0uSWF6NpADNBB1J10gKS/UXwKsDYs9DNwUpm8Cfp1UXilpmKSpwDRgWU9jHEoqy0uoffcgf1i/I+5QnHMDWFZnM0haBFwMFEiqBb4OZAOY2R3AtcCNkhqBg8CCMED9kqRfASuBJuBl4M6w2m8D90v6JFESuT6sb42k+4HXwjKfNTN/mk03fHDGBMbmZVO1fBNzz/DuNufcidGxE40GrrKyMquuro47jH7j9kde4+cvbOSF2y6hYMSwuMNxzvVTklaYWVmqOr9CehBaWFFMY7OxZEVt3KE45wYoTw6D0OnjR1J2ylgWL9/MYDgydM71PU8Og1RlRQlv7djPSxvavX7ROefa5clhkPrLcyYyMjeLqmV+xbRzrvs8OQxSw3Myubp0Mo/VvMPuA41xh+OcG2A8OQxilRXFHGlqYenLPjDtnOseTw6D2IxJo5lVNJoqH5h2znWTJ4dBrrK8hHXv7GXV5l1xh+KcG0A8OQxyV5VOIi8nk6plfjM+51zXeXIY5EYMy+Ijsybxm9Vb2Xe4Ke5wnHMDhCeHIaCyopgDR5p5eJXf4NY51zWeHIaA0uIxnDlhpD8lzjnXZZ4chgBJVFYUs7p2N2u27o47HOfcAODJYYiYf95kcrIyfGDaOdclnhyGiDF5OVwx82QeWrWFg0f8ERnOuY55chhCKitK2HuoiUdfrYs7FOdcP+fJYQh5z9STOLUgn8U+MO2c64QnhyFEEgvKi1m+8V3W1++NOxznXD/myWGIufb8IrIz5QPTzrkOdZocJN0tqV5STTv18yStlrRKUrWkC0P5maEs8doj6QuhbnFS+UZJq0L5FEkHk+ruSN+uOoCCEcO47OwJLFlZy+EmH5h2zqWW1YV57gF+BPy8nfqngYfNzCTNAu4HppvZ60ApgKRMYAuwFMDMFiQWlvQvQPLJ92+aWWm39sJ1S2V5CY+9+g5PrtnGR86dFHc4zrl+qNMjBzN7Dmj3WZNmts+O3Q86H0h1b+hLiL70304ulCTgBmBRlyN2PXbh6QUUjR3uV0w759qVljEHSfMlrQMeBW5OMUslqRPARcA2M3sjqWyqpJclPSvpog62eUvoxqrevn17j+IfajIyxIKyYv64voG3G/bHHY5zrh9KS3Iws6VmNh24Grg9uU5SDnAV8ECKRRdyfNKoA0rM7DzgS8B9kka1s807zazMzMoKCwvTsBdDy3VlRWQIFi/3gWnnXFtpPVspdEGdJqkgqfjDwEoz25Y8r6Qs4BpgcdLyh82sIUyvAN4EzkhnjC4ycfRw/uLM8TywopbG5pa4w3HO9TM9Tg6STg9jB0iaDeQADUmztD46SLgUWGdmRx9wLKkwDF4j6VRgGvBWT2N0qVVWlLB972GeWVcfdyjOuX6m07OVJC0CLgYKJNUCXweyAczsDuBa4EZJjcBBYEFigFpSHnAZ8JkUq041DjEX+KakJqAZuNXM2h0Mdz3zF2cWMmHUMKqWbeJDM06OOxznXD/SaXIws4Wd1H8H+E47dQeAce3UfTxF2RJgSWcxufTIyszg+vOL+cnv17N110EmjRked0jOuX7Cr5Ae4haUF9Ni8EB1beczO+eGDE8OQ1zxSXlcNK2A+6s309yS6hIV59xQ5MnBUVlewpZdB3n+Db9exDkX8eTguOzsCYzLz/Gb8TnnjvLk4MjJyuDa84v477Xb2L73cNzhOOf6AU8ODogGpptajF+t8IFp55wnBxecVjiCiqknsXj5Jo7dR9E5N1R5cnBHVZYXs7HhAC+81dD5zM65Qc2TgzvqinMmMio3ywemnXOeHNwxudmZzD9vMo/XvMO7+4/EHY5zLkaeHNxxKitKONLcwoMvb4k7FOdcjDw5uOOcNXEU5xaPoWqZD0w7N5R5cnBtLCwv5o36fazctCvuUJxzMfHk4Nr4yLmTyM/JpGqZP2PauaHKk4NrI39YFleVTuKR1XXsPdQYdzjOuRh4cnApVZaXcLCxmV+v2hp3KM65GHhycCnNKhrNWRNHUbXcu5acG4o8ObiUJLGwopiaLXuo2bI77nCcc32s0+Qg6W5J9ZJq2qmfJ2m1pFWSqiVdGMrPDGWJ1x5JXwh135C0JanuiqT13SZpvaTXJX0oTfvpTsC80skMy8pgkQ9MOzfkdOXI4R7g8g7qnwbONbNS4GbgLgAze93MSkP5+cABYGnScv+aqDezxwAknQ1UAjPCNn8iKbNbe+TSZvTwbP7ynIn8etVWDhxpijsc51wf6jQ5mNlzwM4O6vfZsaul8oFUV05dArxpZm93srl5QJWZHTazDcB6oKKzGF3vqawoYd/hJh5ZXRd3KM65PpSWMQdJ8yWtAx4lOnporRJY1Krsc6E76m5JY0PZZCD5rm+1oSzVNm8J3VjV27f74y17S/mUsZxWmO/XPDg3xKQlOZjZUjObDlwN3J5cJykHuAp4IKn4P4DTgFKgDviXxOypVt/ONu80szIzKyssLOxR/K59kqgsL2Hlpl38edveuMNxzvWRtJ6tFLqgTpNUkFT8YWClmW1Lmm+bmTWbWQvwU451HdUCxUnLFgF+on3Mrpk9mexM+cC0c0NIj5ODpNMlKUzPBnKA5KfFLKRVl5KkiUlv5wOJM6EeBiolDZM0FZgGLOtpjK5nxo0YxgdnnMzSl7dwqLE57nCcc30gq7MZJC0CLgYKJNUCXweyAczsDuBa4EZJjcBBYEFigFpSHnAZ8JlWq/2upFKiLqONiXozWyPpfuA1oAn4rJn5t1E/sLC8hEdX1/HEmneYV5pyGMg5N4hoMNyWuayszKqrq+MOY1BraTHe//3fUTQmj0W3XBB3OM65NJC0wszKUtV1euTgHEBGRjQw/b0nXuc/fv8mBSNyGDU8m1G52YzMzWJ0mB6Rm0VmRqrzCpxzA4knB9dl159fxJ3PvcV3Hl/X4XwjhmUxKjeLkbnZjBqedTSBjBoefuZmdzCdxbAsv+7Rubh5cnBdNn5ULiv/4TL2Hmpkz8Em9hxqjF4Hm6KyQ03H1SWm39lziDfqw/wHG2nppCczJyvjaKIYmZvNqJBYRqVIJiOP1h2bzs/JJJwj4Zw7QZ4cXLdkZogxeTmMycs5oeXNjANHmlslleMTzJ6D4eehRvaG91t2HTw6fbippcNtZIi2Ry2dHK2Myj2WYEbmZpGV6fekdEObJwfXpySRPyyL/GFZTBx9Yus43NR8NFG0d7SSnFj2Hmpi084DR6f3Hu78PlH5OZkpu8Xanz6+Gy03e/B0jZkZZtGphYkTWKJpMKI6kt4fmw7Lhvd0Um/RDG3WbckxpCpPrD+UtyTPkyLGlqPLtF3WiE6+6CymxDpos61ovS2t1tnu/qRYNrHeFrOUyxLKE3FOLcjn4jPHp/eXjicHNwANy8pk2IhMCkYMO6Hlm1uMfYdSd4slEkiiCywxvWPfEd7asf9oWVMnfWM5mRlHu8US4/PHvkzafikmf3nRUf3Redr7gjz+CyZl/XFfNK3ep5jf9W9XzproycG5dMjMEKPzshmdl31Cy5sZBxub23SL7UmRYPYeaqKlxUDH7g0jCQEKZYn3iRmEkuqS3h8dRmmnPqzr6Fxt6o7VJ7aXalvHlo0mUm2H1utrrz6Utd7vRAySOt3W8W2QNE+K9WekWDa5vTKk45ZFSXEk16dsn2g6o9U6lSru9qaT2q29ZSE6OzBlu4f3GUnryMnqnS5QTw7OdZMk8nKyyMvJ4uTRuXGH41yv8FE355xzbXhycM4514YnB+ecc214cnDOOdeGJwfnnHNteHJwzjnXhicH55xzbXhycM4518ageNiPpO3A2z1YRQGwI03hpJPH1T0eV/d4XN0zGOM6xcwKU1UMiuTQU5Kq23saUpw8ru7xuLrH4+qeoRaXdys555xrw5ODc865Njw5RO6MO4B2eFzd43F1j8fVPUMqLh9zcM4514YfOTjnnGvDk4Nzzrk2hkxykHS3pHpJNe3US9K/S1ovabWk2f0krosl7Za0Krz+sQ9iKpb0O0lrJa2R9Lcp5unz9upiXHG0V66kZZJeCXH9U4p54vp8dSW2Pm+zsN1MSS9LeiRFXSzt1YW4YmmrsO2Nkl4N261OUZ/eNosedD34X8BcYDZQ0079FcBviZ7GdwHwUj+J62LgkT5uq4nA7DA9EvgzcHbc7dXFuOJoLwEjwnQ28BJwQdzt1Y3Y+rzNwna/BNyXattxtVcX4oqlrcK2NwIFHdSntc2GzJGDmT0H7OxglnnAzy3yIjBG0sR+EFefM7M6M1sZpvcCa4HJrWbr8/bqYlx9LrTBvvA2O7xan+kR1+erK7H1OUlFwF8Cd7UzSyzt1YW4+rO0ttmQSQ5dMBnYnPS+ln7wxRPMCd0Cv5U0oy83LGkKcB7Rf5zJYm2vDuKCGNordEWsAuqBp8ys37RXF2KDvm+zHwD/G2hppz6u9voBHccF8f09GvCkpBWSbklRn9Y28+RwjFKUxf4fFrCS6P4n5wI/BB7qqw1LGgEsAb5gZntaV6dYpE/aq5O4YmkvM2s2s1KgCKiQNLPVLLG1Vxdi69M2k3QlUG9mKzqaLUVZr7ZXF+OK7e8ReJ+ZzQY+DHxW0txW9WltM08Ox9QCxUnvi4CtMcVylJntSXQLmNljQLakgt7erqRsoi/gX5rZgylmiaW9OosrrvZK2v4u4PfA5a2qYv98tRdbDG32PuAqSRuBKuADkn7Rap442qvTuOL8fJnZ1vCzHlgKVLSaJa1t5snhmIeBG8OI/wXAbjOrizsoSSdLUpiuIPqdNfTyNgX8DFhrZv+3ndn6vL26EldM7VUoaUyYHg5cCqxrNVssn6+uxNbXbWZmt5lZkZlNASqBZ8zsY61m6/P26kpccXy+wrbyJY1MTAMfBFqf4ZjWNss64WgHGEmLiM40KJBUC3ydaHAOM7sDeIxotH89cAD4RD+J6zrgf0hqAg4ClRZOTehF7wP+Cng19FUD/B1QkhRXHO3VlbjiaK+JwL2SMom+LO43s0ck3ZoUVyyfry7GFkebtdFP2quzuOJqqwnA0pCXsoD7zOzx3mwzv32Gc865NrxbyTnnXBueHJxzzrXhycE551wbnhycc8614cnBOedcG54cnOsGSd+Q9OW443Cut3lycM4514YnB+c6IOnGcG/8VyT9V6u6T0taHuqWSMoL5ddLqgnlz4WyGYqeq7AqrG9aKP9YUvl/KrpJXqake8I6XpX0xb7fczfU+UVwzrUj3HHzQaIbnu2QdBLweWCfmX1f0jgzawjz/jOwzcx+KOlV4HIz2yJpjJntkvRD4EUz+6WkHCATmAJ8F7jGzBol/QR4EVgDfNvMLgvrHhPui+Rcn/EjB+fa9wHgV2a2A8DMWj93Y6ak50My+CiQuH3zH4F7JH2aKAkAvAD8naSvEN3V8yBwCXA+sDzcDuQS4FTgLeBUST+UdDnQ+s6zzvU6Tw7OtU90fMvje4DPmdk5wD8BuQBmdivwNaI7ZK4KRxj3AVcR3Y/nCUkfCOu/18xKw+tMM/uGmb0LnEt0B9XPMjAfPOMGOE8OzrXvaeAGSeMAQrdSspFAnaLbiH80USjpNDN7ycz+EdgBFEs6FXjLzP6d6O6Zs8L6r5M0PrF+SacougV0hpktAf6B6DGyzvWpIXNXVue6y8zWSPoW8KykZuBlouf4JvwD0ZPo3gZeJUoWAN8LA84iSgCvAF8FPiapEXgH+KaZ7ZT0NaKne2UAjURHCgeB/xfKAG7rxd10LiUfkHbOOdeGdys555xrw5ODc865Njw5OOeca8OTg3POuTY8OTjnnGvDk4Nzzrk2PDk455xr4/8HmhkdPaSLatUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "optModel_idx = np.argmin(model_df[\"BIC\"].values)\n",
    "optModel_components = model_df[\"classes\"].values[optModel_idx]\n",
    "model_df.plot.line(x=\"classes\", y=\"BIC\")\n",
    "plt.title(str(\"The optimal model has %i modes\" %(optModel_components)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "754c43d2",
   "metadata": {},
   "source": [
    "## Train optimal Model\n",
    "Parameterize the optimal model by training with the jump distance sequence.\n",
    "Store the optimized model parameters to optModel_diffusion_coefficients and optModel_weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "170eaa5d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done.\n"
     ]
    }
   ],
   "source": [
    "x = data_df[\"jump_distance\"].values\n",
    "jmm = em.JumpDistanceMixtureModel(n_components = optModel_components,\n",
    "                                  degrees_of_freedom = 4,\n",
    "                                  tau = 0.02,\n",
    "                                  init_params = \"wm\",\n",
    "                                  params = \"wm\")\n",
    "jmm.fit(x)\n",
    "optModel_diffusion_coefficients = jmm.diffusion_coefficients()\n",
    "optModel_weights = jmm._weights\n",
    "print(\"done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11cc5037",
   "metadata": {},
   "source": [
    "## Plot optimal model\n",
    "plot the optimal model and superimpose it to the pdf of the original jumpdistance distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "74c8f0ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEHCAYAAAC5u6FsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABq8klEQVR4nO3dd3iN5xvA8e+TvfeQRexN7L1qb0pr1q5dW1Gl1KyiZqk9i1ZbtBSlVGuv2CtmEitChkTmeX5/ZPyQHTk5STyf6zqX5Lzrfs+Rc593PPctpJQoiqIoSnrp6ToARVEUJXdRiUNRFEXJEJU4FEVRlAxRiUNRFEXJEJU4FEVRlAwx0HUA2cHBwUF6enrqOgxFUZRc5ezZs8+klI5vP/9eJA5PT0/OnDmj6zAURVFyFSHE/eSeV6eqFEVRlAxRiUNRFEXJEJU4FEVRlAx5L65xKHlHdHQ0fn5+RERE6DoURckzTExMcHd3x9DQMF3zq8Sh5Cp+fn5YWlri6emJEELX4ShKrielJDAwED8/PwoWLJiuZdSpKiVXiYiIwN7eXiUNRckiQgjs7e0zdBSvEoeS66ikoShZK6N/UypxKIqiKBmiEkcu9zwsStchKEqiFi1aEBQUlOo8FhYWyT7fq1cvtm/froWolKymEkcu5u0bRJsl/+k6DEVBSolGo2HPnj3Y2NjoOhxFy1TiyMXOP3iB34tXBIdH6zoUJY8YN24c33//feLvU6ZMYerUqTRs2JCKFStStmxZdu7cCcC9e/coWbIkgwcPpmLFivj6+uLp6cmzZ88AaNeuHZUqVaJ06dKsWLHije2MHj2aihUr0rBhQwICApLEcfbsWerVq0elSpVo2rQpjx490uJeKxmlEkcudu7+CwCuPQ7RcSRKXtG5c2e2bduW+PtPP/1E7969+e233zh37hyHDh1i9OjRJLScvnHjBj169OD8+fMUKFDgjXWtWbOGs2fPcubMGRYtWkRgYCAAYWFhVKxYkXPnzlGvXj2mTp36xnLR0dF89tlnbN++nbNnz9KnTx8mTpyo5T1XMkKN48iFpJQIIbjgF0wpF0uuPgyheiF7AK4/DsHa1BAXa1MdR6nkRhUqVODp06c8fPiQgIAAbG1tcXFxYeTIkRw5cgQ9PT38/f158uQJAAUKFKB69erJrmvRokX89ttvAPj6+nLr1i3s7e3R09OjU6dOAHTv3p0PP/zwjeVu3LjB5cuXady4MQCxsbG4uLhoa5eVTFCJI5c5eSeQmXuusbZ3VQJCI2lU0oPL/sGJ0xceuIWVqSHfdCinwyiV3Kxjx45s376dx48f07lzZzZv3kxAQABnz57F0NAQT0/PxHv+zc3Nk13H4cOHOXDgAMePH8fMzIz69eunOE7g7VtBpZSULl2a48ePZ+2OKVlGnarKJZ4+fcqsWbP4sFUz/pzShVLlKxB2YAlB109w2e954nx3n4Wx++IjIqJjdRitkpt17tyZrVu3sn37djp27EhwcDBOTk4YGhpy6NAh7t9PttL2G4KDg7G1tcXMzIzr169z4sSJxGkajSbx7qkff/yR2rVrv7Fs8eLFCQgISEwc0dHRXLlyJQv3UHlX6ogjh5NSsnz5csaOHUtYWBgm+QpTqlRpbvg+RV48wncndmPkkJ8/S6+gWdMmPHgejou1CQevPaVlOXV4r2Rc6dKlCQ0Nxc3NDRcXF7p160br1q2pXLkyXl5elChRIs11NGvWjOXLl1OuXDmKFy/+xuksc3Nzrly5QqVKlbC2tn7jmgqAkZER27dvZ9iwYQQHBxMTE8OIESMoXbp0lu+rkjki4SJXXla5cmWZGxs5SSn57LPPWLp0KWWr1aHixyMxdvCgZ01PZu65RrNSjgRdPcqMrycTEfiQ9h0/5kbRLvSoU4Jrj0JY16eqrnchy127do2SJUvqOgxFyXOS+9sSQpyVUlZ+e151xJGDffHFFyxdupSmnfsR7tWZGCtz2pV1QQjBFy1Kxp0bLtiOq0bFiD2/g+2rFmDqeJJ8tX9kzZ1gomM1GOqrs5GKomQt9amSQ/3666/Mnj2bAQMGYFWvNx0r52dIg6I4W5kAb15QLO5mh1vDT5i4ZDOx4SEM69YK4xd3ufssDID34ahSUZTsoxJHDvT06VM+/fRTKleuzLRv5nH5YQhl3axTnL9kPitO3nmOuacX3WduxMLSimtrxvLLH/t5GRlDma/2UWX6AX4755eNe6EoSl6lEkcO9PnnnxMaGsqGDRs4djeI0q5WmBjqpzi/u60pLyNjOHb7GcWLFmHZlj+wsndmytAebNr1FwXszelUxYNtZ3yzcS8URcmrVOLIYc6dO8f69esZNWoUJUuWZPfFR5T3sEl1GSEEpVytuOAXTD5rExzzuTBs7gaMrewY1aczthEP8fKw4YJfMJEx/79N96JfkDqNpShKhqnEkcNMnjwZW1tbxo8fz4IDN7nyMIRqnvZpLlfMyRIAF+u4ayBlinpSpNc3YGDM7m+HE/4iAHcbUy76xQ0WPHD1CW2WHMXvxSvt7YyiKHmSShw5yPnz59m9ezdjxoxh360Qfj3nz5ctS2JhkvbNb6VdrbA1M8TMKG7efFYmhBna4vzxFKJevWRs/64UtjHg+O1AnoRE8PkvF7ExNeRpqOrd/T7YtWsXs2fPBmDHjh1cvXo1cdrkyZM5cOCArkJTciE1jiMH6dGjB7/++iuXb9yh/WpvPm9WAk/75Es6JOdlRMwbSebLHZfQSGhp84hxA7pTpkZD3DpOJEoDlQvYci8wjB41PGlRNvcMFHyfxnFIKZFSoqeXtd/vevXqRatWrejYsWOWrlfJ3TIyjkMdceQQT548YevWrfTp04dlJx5TvZB9hpIGkOTIxMPWjGLOFtRq0ITPvpjGxaN/cWzbErzcbWhT3jXuiCNEHXFkVFhYGC1btqR8+fKUKVOGbdu2vVFO/MyZM9SvXx+IK0v+ySef8MEHH1C0aFFWrlyZuJ5vv/2WKlWqUK5cOb766isg+VLlFhYWyZYh9/b2pnr16pQrV4727dvz4kVcteRFixZRqlQpypUrR+fOnQFYt24dQ4cO5dixY+zatYuxY8fi5eXF7du332igdPDgQSpUqEDZsmXp06cPkZGRAHh6evLVV18llla/fv269l9oJcdSAwBziC1bthAdHU2XHn3pt+shCzpXeOd1tiznkjgA8OOe/Xlw5xa//bgOu4CWCOGOlakRj3Nx4pg3bx43btzI0nUWL16c0aNHpzrP3r17cXV1Zffu3UBcXaZx48alOP/Fixc5ceIEYWFhVKhQgZYtW3L58mVu3brFqVOnkFLSpk0bjhw5Qv78+blx4wZr165N7IuRUIZ83rx5fP3110ydOpUlS5bQo0cPFi9eTL169Zg8eTJTp05lwYIFzJ49m7t372JsbJykG1/NmjVp06ZNskccERER9OrVi4MHD1KsWDF69OjBsmXLGDFiBAAODg6cO3eO77//nrlz57Jq1aoMvrpKXqGOOHKIjRs3UqlSJXxibCjvYYOF8bvndHdbszcGDA6fOIPS5SsxY9xnPLjrg62ZIY+Cc2/i0JWyZcty4MABxo0bx7///ou1dcpjbADatm2LqakpDg4ONGjQgFOnTrF//372799PhQoVqFixItevX+fWrVtA0lLlb5ch/++//wgODiYoKIh69eoB0LNnT44cOQJAuXLl6NatG5s2bcLAIP3/j27cuEHBggUpVqxYknUCieXPK1WqxL1799K9XiXvUUccOcCVK1c4d+4cCxcuZPsZX5qUyqeV7RgZGzN98Rp6tW3AF0N6M2zBNp6GRGplW9khrSMDbSlWrBhnz55lz549TJgwgSZNmmBgYIBGowFIUj787bLhQgiklEyYMIEBAwa8Me3evXsplipPaX1v2717N0eOHGHXrl1MmzYt3ZVl07reaWxsDIC+vj4xMTHpWqeSN6kjjhxg48aN6OvrU7Nxa+4/D6ecR+rfYN+Fs6sbU79bwd1b19m+6CuehKjbcTPq4cOHmJmZ0b17d8aMGcO5c+fw9PTk7NmzAPzyyy9vzL9z504iIiIIDAzk8OHDVKlShaZNm7JmzRpevnwJgL+/P0+fPk12e8mVIbe2tsbW1pZ///0XiPs/VK9ePTQaDb6+vjRo0IA5c+YQFBSUuI0ElpaWhIaGJtlOiRIluHfvHj4+Pm+sU1Hepo44dCw2NpZNmzbRvHlzdvuEU7uIAwZZfBfN26rWrk+/EeNZ+d0sXPVcYUwDrW4vr7l06RJjx45FT08PQ0NDli1bxqtXr+jbty8zZ86kWrVqb8xftWpVWrZsyYMHD5g0aRKurq64urpy7do1atSoAYCFhQWbNm1CXz9phYCUypCvX7+egQMHEh4eTqFChVi7di2xsbF0796d4OBgpJSMHDkSGxubN9bXuXNnPv30UxYtWpSYkABMTExYu3YtH330ETExMVSpUoWBAwdm8aun5AXqdlwdO3r0KLVr12bF2g0suufE7A/LYWdupPXtajQaxn7aleP/HebE8RNUq1JJ69vMCrntdtwpU6ZgYWHBmDFjMr0OCwuLJEcNipLV1O24ucju3bvR19fHz7wYNQvZZ0vSgLgLrpO+XYqhmTXdunYlLCwMKSUxsZps2b6iKLmXShw6tnv3bmrXrs1vV4Nons0D8Wzs7Cn/yZfcuX2LESNG8MfFRwzadC5bY8jrpkyZ8k5HG4A62lByHHWNQ4d8fX25ePEin0+exj9GBom3zmanIuWrU6jvZ6xatYgnlsW4aV6GyJhYjA1SrsarKMr7TR1x6NCePXsAcC5dg6JOFjqJwcrUgEbdh1CtWjX2LJuKDH3K2fsvdBKLoii5g0ocOrRnzx48PT15LBwo5JCx8iJZxcbUiGdhsSxbvR4pJaF/fseha490EouiKLmDVhOHEKKZEOKGEMJHCDE+melCCLEofvpFIUTFtJYVQngJIU4IIbyFEGeEEFW1uQ/aEhERwYEDB2jRogXnfYMo6mypkzjcbEw5evsZT7GmQpexPLp5gTWL5xIeFcOTXFyORFEU7dFa4hBC6ANLgeZAKaCLEKLUW7M1B4rGP/oDy9Kx7BxgqpTSC5gc/3uuc/z4ccLDw6nfqAm+L8IpYGemkzgqedoSGa1h9p/XadCiPc0/7MydvzZRftBi+qw7rZOYFEXJ2bR5xFEV8JFS3pFSRgFbgbZvzdMW2CDjnABshBAuaSwrAav4n62Bh1rcB605cuQIQggem3pS0MEcA33dnDXUE4IeNTy5+yyMYs6WjJo0C3tnV0L3L+TOw2cEhUfpJC5FUXIubX5auQGvN7n2i38uPfOktuwI4FshhC8wF5iQ3MaFEP3jT2WdSShDnZMcOvwPRk4F2XczhBqF0u7wp01FnCyY2KIkJV2sMLe05Ov53/PU/wGRxzZw/HagTmN73+3du5fixYtTpEiRxEZMb+vTpw9OTk6UKVPmjecTSo+ULFmS0qVLs3DhwjS3p6+vj5eXF6VLl6Z8+fLMnz8/sQbXu0opzgQRERFUrVqV8uXLU7p06cRS85ndl/RI7fWNjY2lQoUKtGrVKtll7927h6mpKV5eXlkSS4KYmBhatmyJg4MDly9fTnw+rddgwIABHD16NNl1vnr1Ci8vL4yMjBLL/78LbSaO5CqxvT1MPaV5Ult2EDBSSukBjARWJ7dxKeUKKWVlKWVlR0fHdIacPaKiojhx4jjuJSsysWUpGmupqGFGlHGzRl8v7mWvULUmXfoO5u6/O1i77TcdR/b+io2NZciQIfz5559cvXqVLVu2vNG5L0GvXr3Yu3dvkucNDAyYN28e165d48SJEyxdujTZ5V9namqKt7c3V65c4a+//mLPnj1MnTo1S/YnpTgTGBsb8/fff3PhwgW8vb3Zu3cvJ06cyPS+JDh8+DC9evVK8nxar+/ChQvTrFJQuHBhvL290xVHeg0aNIjixYuzc+dOOnXqhJ+fH5D2a3Dy5Mk3qiq/LuF9dXV1zZIYtZk4/ACP1353J+lppZTmSW3ZnsCv8T//TNxprVzl3LlzREZEUL5KDV2HkqJPR07Ao3Bxdi6aTGCgOupIza5du5L0tli2bBnDhg17p/WeOnWKIkWKUKhQIYyMjOjcuTM7d+5MMl/dunWxs7NL8ryLiwsVK8bdb2JpaUnJkiXx9/dP9/adnJxYsWIFS5YsSbNybnqkFGcCIQQWFnG3pUdHRxMdHZ1YCTi1fdm0aRNVq1bFy8uLAQMGEBsbm654Unt9/fz82L17N/369Uv3/tWvXz+xP0xgYGDikdW9e/coUaIE/fr1o0yZMnTr1o0DBw5Qq1YtihYtyqlTpxLXMXXqVKytrZk/fz61atVi1apVdOnSheDg4FRfg2vXrlGsWDH09fWTbTSW1bQ5APA0UFQIURDwBzoDXd+aZxcwVAixFagGBEspHwkhAlJZ9iFQDzgMfADc0uI+aEVCj4NGDerrNI7UGBub8PX8ZfRp34SPe/TjwB+/plnOWxf69++f7PMrVqwAUm72NHr0aIoXL87vv//O77//nmS5jJg4cSJbtmx547nChQsnqZILUKdOnWQr086dO5dGjRq98Zy/vz8eHv///uTu7s7JkyczHB/EfXidP38+SQHGtBQqVAiNRsPTp09xdnZ+Y1pG9iW9YmNj4/rS+PgwZMiQZON9fV+uXbvGtm3bOHr0KIaGhgwePJjNmzfTo0ePNLeV2us7YsQI5syZk+z+pcTHx4eiRYsCcc27ypYt+8a0n3/+mRUrVlClShV+/PFH/vvvP3bt2sXMmTPZsWMHwBun5wBq1KiRWAE5pdcA4M8//6RZs2ZA8o3GsprWEoeUMkYIMRTYB+gDa6SUV4QQA+OnLwf2AC0AHyAc6J3asvGr/hRYKIQwACKIuxsrV/nr78MY27tTrmh+XYeSquKly9G0x1D2rltAzU+/5tjKyTkyeejShQsX0Gg0lClThvv377Nnzx4GDRr0xrfl1yX3IZCS5L7lZ+b1f/nyJR06dGDBggVYWVmlvUA64oCM7Ut66evr4+3tTVBQEO3bt+fy5ctvXBN5e182bNjA2bNnqVKlChB3Lt/JyQmAatWqERkZycuXL3n+/HnitYhvvvmGpk2bpvj6/vHHHzg5OVGpUiUOHz6crrjv37+Pm5tbYn/4ixcvUq5cucTpBQsWTEwkpUuXpmHDhgghKFu2bIabYiX3fu7bt4+1a9cCcY3GxowZw7hx42jVqhV16tTJ0PrTQ6slR6SUe4hLDq8/t/y1nyUwJL3Lxj//H5A7SrkmIzY2lmNHj+JWoQF6ueBD+Ivx4/H1/o8zP87ll57t6VinXNoLZaO0jhDSavbUunVrWrdunente3t7U6lS3H/Hv/76K7GL39WrVylfvnyS+TPyLd3d3R1f3//fI+Ln55fhc9TR0dF06NCBbt26JXbwy4g7d+6gr6+f+GH8Om0ccSSwsbGhfv367N27NzFxJLcvUkp69uzJrFmzkqwj4ejh8OHDrFu3jnXr1r0xPaXX9+jRo+zatYs9e/YQERFBSEgI3bt3Z9OmTSnG6+3t/UaiOHv2bGLXRvh/EyyIKzCa8Luenl6GmmIl9xqEh4cTFBSU+H8juUZjkydPTvc20kPVqspmV65cIfxlCOUrZeyUga4YGBgwee5Suresz+D+/fjwyonEb1VKXHn6ly9fEhsby6+//oqbmxuvXr1i3bp1bNy4Mcn8GfmWXqVKFW7dusXdu3dxc3Nj69at/Pjjj+leXkpJ3759KVmyJKNGjXpjWsOGDdmwYQNubm/f6Ph/AQEBDBw4kKFDh77z0VN6BAQEYGhoiI2NDa9evUpsz5vavjRs2JC2bdsycuRInJyceP78OaGhoRQoUCDN7aX0+pYuXToxER0+fJi5c+emmjQg7sgzofPjrVu32LlzJ9OnT8/sS5GslF6DQ4cO0aDB/3vqPHz4EDs7O7p3746FhUWShJkV1CdANku4EFa/Ti0dR5J+Hp6FGTZhKgHXTzPpm0W6DidHadGiBXfu3MHLy4uBAwdy5coVKleuTP/+/RMvZGaWgYEBS5YsoWnTppQsWZKPP/6Y0qVLJ2734cO4+0W6dOlCjRo1uHHjBu7u7qxeHXej4dGjR9m4cSN///03Xl5eeHl5sWfPHjQaDT4+PsleqE64bbN06dI0atSIJk2aJDnvnlkpxZmwL48ePaJBgwaUK1eOKlWq0Lhx48RbYVPal1KlSjF9+nSaNGlCuXLlaNy4MY8epa9kTmqvb0Z5e3uj0WgoX748X3/9NSVLlmT9+vWZWldKUnoNXr++AXGNxhJuFpgxYwZffvlllsYBqpFTtuveux9btm7lyMW7ueqbu0aj4eO2LXh+/zo3r19946JidsptjZxyosuXL7NmzRrmz5+v61BytXv37tGqVSsuX75MkSJFOH/+PJaW2V86qGLFipw8eRJDQ8M05/X09OTMmTM4ODgkmaYaOeVgx0+ewtGzZK5KGhB3LrbjiGnExsYwYMCALLk9U9GNMmXKqKSRBfT19QkODqZw4cLo6enpJGlA3O39aSWNhCPJ6OjoLPnsyV2fXrlcZGQk925eo2jppBdNc4MyJYpRsvUA/vzzz2TP3yvK+8TDwwNfX19u377NzZs3dR1OqhIGAPr7+6c6lia9VOLIRhcvXkQTG0O1qrluzCIA7ramyFLNqFWrFsOHD0/3uWRFUfIWlTiy0b/H4m4PrFmtio4jyRxrU0M0COYtXk5ERASDBg1Sp6wU5T2kEkc2OnLsBIbm1uRz082F5XclhMDDzgysXfj666/ZuXOnVsoZKIqSs6nEkY28z5/FoWCpXD362tXahFtPX3LOqiYmrsXo9ekgAgICiIrR8MM/t3UdnqIo2UAljmwSHh7Og9s3yV8s+ZLSuYWLtSkr/73D49BoFny/gshXLxkyZCh3n4Ux68/rhEWmfxSsoii5k0oc2cTb2xup0VC0VM4q2ZFRbrYm3A8M55PqBShXpgyu9bvz888/seXnuIJ+N5+kvyicoii5k0oc2eTixYsAlM/ipi/ZrbizFf3rFKKkS1xxtXIte1C4eCkWTpuAJjKcG49V4lCUvE4ljmxy6dIl9EzMKVm0kK5DeSemRvo0KPH/gnf5bMzpPW4WQc+eEHF8M1cfhegwOkVRsoNKHNnkvPcFjB0KYGtmpOtQspSDhTFGrsUo3qADAad28d+x47oOKU9Jq3Vsau1WAYKCgujYsSMlSpSgZMmSHD+e+vujzdaxnp6elC1bFi8vLypXTlLFIsXWqGntY2al9dqm1jo2p7WNhZRbx2Z121hQiSNbSCm5fPkStu5FcvUdVclxtDTm3rNw7Ov3xMbBmePrZhEVFaXrsPKE9LSOTa3dKsDw4cNp1qwZ169f58KFC2nW+dJm61iIq+Tq7e1NcrXjUmqNmtY+piW51rHpeW3Tah2bk9rGQsqtY7O6bSyoxJEt/Pz8CA0JwaVQcV2HkuWcrUy4FxhGQKQ+o7+aTcTTu0ybNUfXYWUrXbaOTa3dakhICEeOHKFv374AGBkZYWNjk+7tZ3Xr2LSk1Bo1tX2EzLWOTeu1zarWsdnRNhb+3zo2IiJC621jQSWObJFwYdyzSN5LHI6Wxlx/HIqliQENm7UiX7m6zJk1g9u3s2dMR//+/RNbv2bVzxk1ceJEpkyZ8sZzhQsXTvKNEOKaHyWUxH79ceDAgSTzJtfaNLme4bGxsXh5eeHk5ETjxo0T24neuXMHR0dHevfuTYUKFejXrx9hYWEZ2rfXW8e+y75AXJJr0qQJlSpVSrMB19utUVPax9dbx3p7e6Ovr8/mzZvT3K+0XtuE1rEZKQiYUutYHx8fhg8fzsWLF7l+/Xpi29i5c+cyc+bMxOW/+uor5s6dm/h7QttYa2vrVF8b+H/r2IS2sRcuXODy5ctvlFvPSqqRUza4dOkSAFUqeek2EC2wMTVEX0+Qz8oEgCb9JrBtfEcGDhzInj/3IhEYGeTd7yc5oXVsSu1WY2JiOHfuHIsXL6ZatWoMHz6c2bNnM23atAztY1a1jj169Ciurq48ffqUxo0bU6JECerWrZtkvuRao6a0jwcPHsxU69jUXtusbh2r7bax8P/WsS9fvtR621hQiSNbnDnvjaG1E9WL5+we45khhMDFygQny7hWmC2rl+bvBr05sHsJTYbMoHGbjnzRspTWtv/6N9es+jkjclLr2Lfbrbq7u+Pu7p74zbRjx47JXgROTVa2jk2I3cnJifbt23Pq1KkkiSOtVrdv72NmW8ceP348xdc2q1vHarNtLCRtHavttrGgTlVli1NnvXHMXyTPfvN2sjLGwSLuj8HDzowqzT/GuUhZ/t00j71nb+k4Ou16u3VsaGhoYuvYrl27Jpn/33//xdvbO8kjuQ/a11ubRkVFsXXrVtq0afPGPAEBAQQFBQEktlstUaIEAPny5cPDwyPxvPvBgwcpVSouiTds2DDZ015vrzut1rHp3ZewsLDEJBMWFsb+/fsTe4knSKk1amr72LBhQ7Zv3554Ku358+fcv38/1f2C1F/bWbNm4efnx71799i6dSsffPBBplrHJhxlZIXU2gC/3jr24cOHmJmZ0b17d8aMGcO5c+eyLIbX5c1PshwkKioKv3s+lCiVuZaUuUHNwg54edgk/t6ynBt6dQcgI8O59OsSnoRE6C44LdN169jU2q0CLF68mG7dulGuXDm8vb354osvdNI69smTJ9SuXZvy5ctTtWpVWrZsmXj+PWFfUmqNmto+ZrZ1bFa2jQXtt45N6bUB3mgdmx1tY0G1jtW6M+e8qVKpAhPmfE/rDz/WSQzZTUqJz9OXHNy4kA3LFzD5+y1MHdQ5S9atWse+O9U69t293jYWyBWtY1NrGwuqdWyOcvLcBQBK5uEjjrcJISjqbEnvoaOxc/Fg0dfjEg/jFd1TrWPfXULbWC8vL0JDQ3N069isbhsLKnFo3eWr10AIPAoW1nUo2c7YxJTPvpxN0OMHzJiZ9OKlouRWCW1jvb29sbS0zNGtY7O6bSyoxKF1N25cx8zeBWNjE12HohNNmzYlf5XGzJw1i2vXrus6HEVRsoBKHFp21+cWti6eug5Dp76bPx89A2Padu2FlJKfz/jyyeqThEep3h2KkhupxKFFGo0G/3u3sXcvqOtQdMrFxYUewydwy/sk69dvYN2xezwJieCT1aeIismaAnqKomQflTi06MGDB0RHRZLPI3eXUs8Kffp+ilWBUnw2chR+j54yqVUpAkIjuabKsCtKrqMShxZdvx53Tt+9YFEdR6J7enp6dB01jZfBQUSf2IiBnh757cy48+ylrkNTFCWDVOLQooTEUbCIShwAHzashW319lw7vBPv08fJZ22CzxOVOBQlt1GJQ4uuX7+OobkV+ZyT1vl5H1mYGLBh8RzyuXkw58vROJrpcStAJQ5FyW1U4tCi69evY+rggbVp6iM63yf2NtaM/mo2927f5Nzujdx5mrEy34qi6J5KHFp0/fp1DOzcsTRRieN1tT5oSv2mrdmxZjE+d24Tq8n7ZW8yK632ppByS9bMtFzVZuvYPn364OTklKS44eu+++47SpcuTZkyZejSpcsbFQcy2gY3Lam9tmm1uc3u1rGQdvvYlFrHQta3j1WJQ0tevHjBkydPwMYNC2NVvf5tIyfNxMDAgBd/LcPvebiuw8mR0tPeNEFyLVkz03JVm61je/Xqxd69e1Oc7u/vz6JFizhz5gyXL18mNjaWrVu3Jk7PaBvcBJltHZtam1vI3taxkHb72JRax0LWt49ViUNLEkpZWzrnR18vb/UZzwqO+VzoP+oLQm6dYdWG1EtW53S6bB2bmrRarqYlq1vH1q1bN82SFzExMbx69YqYmBjCw8MTP+hSa4OrjdaxmZFc61gg3e1jU2sdCym31oU81jpWCNFMCHFDCOEjhBifzHQhhFgUP/2iEKJiepYVQnwWP+2KECJHNrhOqF1j7+ap20BysA7d++JYsCSLZnyZ2G8ho/r37//GI6H967p16+jfv39i857ff/89ybwAz549S/w9s4fwum4dm1pL1pRarqZXVraOTYubmxtjxowhf/78uLi4YG1tTZMmTYCU2+Bqq3VsRtrcJkipdWzCtLTax6a3dSwkbR+bZ1rHCiH0gaVAY8APOC2E2CWlfP2vqTlQNP5RDVgGVEttWSFEA6AtUE5KGSmEyJG3LPn4+KCnp4eji0faM7+n9PX1aT9sKitHdeaLL77g+++/13VIGZYTWsem1pI1pZarGZFVrWPT8uLFC3bu3Mndu3exsbHho48+YtOmTXTv3j3FNrjOzs5Z3joW0t/mNkFqrWNB++1j81Lr2KqAj5TyDoAQYitxH/ivJ462wAYZ9y6eEELYCCFcAM9Ulh0EzJZSRgJIKZN+FcoBfHx8cHJxx9rCVNeh5GjtGtZmV412LF++nJ49e2b4G3FK3wZ79er1xnnt1q1b07p16yTzOTg4ZLptLOSM1rHpacn6dsvV9MrK1rFpOXDgAAULFsTR0RGADz/8kGPHjtG9e/cU2+A2a9Ysy1vHQvpe09el1joWtNs+Nq+1jnUDfF/73S/+ufTMk9qyxYA6QoiTQoh/hBBVktu4EKK/EOKMEOJMQEDAO+xG5vj4+ODgVgBLE3VhPDU2ZkZ8N2cWBpb29Ojdl+joaF2HlCG6bh2bWkvWtFquZmfr2PTInz8/J06cIDw8HCklBw8eTLwAnlIbXG20jk1Pm9u3abt1LKTcPjavtY5N7irc28eHKc2T2rIGgC1QHRgL/CSS+V8tpVwhpawspayc8A0mO/n4+GCbzwNzdUdVmgrks6d5/y+4ee0KCxcuJPhVdJZcjM0Oum4dm1pL1pRaruqidSxAly5dqFGjBjdu3MDd3Z3Vq1e/sS/VqlWjY8eOVKxYkbJly6LRaBKvRUHybXC10To2tdc0JdpuHQspt4/NU61jhRA1gClSyqbxv08AkFLOem2eH4DDUsot8b/fAOoTd6oq2WWFEHuJO1V1OH7abaC6lDLFw4rsbh37/Plz7O3tadRrDNXa9KRlOZds23Zu9fBFOH27f0y070UKDvyBzzvUpkdNzyTzqdax7061jn13ubF1LKTePjantI49DRQVQhQUQhgBnYFdb82zC+gRf3dVdSBYSvkojWV3AB8ACCGKAUbAu49oyUI+Pj4A3Im0wCu/jW6DySVcbc0o3XEEkTEaAvcvZ/5fNwgKj9J1WHmSah377nJT61jI+vaxWjuPIqWMEUIMBfYB+sAaKeUVIcTA+OnLgT1AC8AHCAd6p7Zs/KrXAGuEEJeBKKCnzGHnNRISR8HChXGzURfH0+uDyqW5W/8THu1fSZkGF5m7z5Xp7bP2PLGiZIWE1rEJcnLrWPj/AMCsotUT8FLKPcQlh9efW/7azxIYkt5l45+PArpnbaRZ69atWyAEbeu82znu903jUs6Um/YFE32P4/3TdwTaFKd2UUealcmn69AURXmNGjmuBRev3sDQyoFKhXLkEJMcy1BfD3d7Sz6fNo/Ap49xvLGD8b9c5EGgKkmiKDlJuhKHEOIXIURLIYRKNOlw49YtrJw8MlTeQfm/0l6VaNelFwe2r6e86Qtm7rn2xvQcdmZSUXK9jP5NpTcRLAO6AreEELOFECUyGtj75MG9uzi45td1GLnawDFfYmvvwIkNszh1JwBv3yAATExMCAwMVMlDUbKIlJLAwEBMTEzSvUy6rnFIKQ8AB4QQ1kAX4C8hhC+wEtgkpcxdo7a0KDg4mNAXgeTz8NR1KLmapZU1w7+cweThn9Ky1hHm7bdlY99quLu74+fnhy4GdSpKXmViYoK7u3u650/3xXEhhD1xF6U/Ac4Dm4HaQE/ixl4owO3btwHwLFhIx5Hkfg1btGP39i0c2rwYd7PSSFkVQ0NDChYsqOvQFOW9lt5rHL8C/wJmQGspZRsp5TYp5WeAhTYDzG0SbsUtWqSIjiPJ/YQQjJk6h5iYGJ7sX87jkIi0F1IURevSe41jlZSylJRyVvwAPYQQxgDJjSp8n925cweAUiWK6jiSvMEtvyd9ho4h+Op/rN/2m67DURSF9CeO6ck89259G/Oomz530De1wtnORteh5Bld+g7G1q0Qc7/6nLAw1aNcUXQt1cQhhMgnhKgEmAohKgghKsY/6hN32kp5y41btzGzd1G34mYhQyMjPvxsCoGP/bOsjamiKJmX1sXxpkAvwB14vbhNKPCFlmLK1e7du4eNs6euw8hzatWuzf4arZg/fz7du3d/o/eBoijZK9UjDinleillA6CXlLLBa482UspfsynGXEOj0fD0kS8OqutflnO3NcO0Vg9sbO1o16UnGo1G1yEpynsrrVNVCTWhPIUQo95+ZEN8ucqjR4+IiY7G1V0ljqxmYWyAhZUNTo37cfeqNzPmLdZ1SIry3krr4rh5/L8WgGUyD+U1d+/eBcAjv6duA8mjPOzMKFm7Je6lqzBj6iQeP36s65AU5b2U6jUOKeUP8f+qK5LpkJA4Cnh66jaQPKpv7YJYGBtw3W0eQzo05LNhw/n5p226DktR3jvpHQA4RwhhJYQwFEIcFEI8e+00lhIvIXEUKuip20DyKAcLY0wM9fEqU4qSLXqz/eefWLZ2s67DUpT3TnrHcTSRUoYArQA/oBhx/b6V19y9excDS3scrNVZPG2bN20SzgVLMuyzoTx68lTX4SjKeyW9iSOhL2ELYIuU8rmW4snV7ty9i4G1MyaGqvq8tllbmDJ36QpiI17So+8AXYejKO+V9H7C/S6EuA5UBg4KIRwBVTjoLXfu3MVcDf7LNoWLlaRGxwEc2L2D7du36zocRXlvpCtxSCnHAzWAyvEl1MOAttoMLLeJjo7mob8f1k5uug7lvdJ74HCsPIozePBgVWpdUbJJRs6plAQ6CSF6AB2BJtoJKXfy9fVFo9Fgm08ljuxU3M0Gh5YjCQ4OZsiQZNvXK4qSxdJ7V9VGYC5x/TeqxD9UVdzXJNxR5ahGjWcrAz09KnuV46P+o/j5558Z9+1ytp/x1XVYipKnpbeRU2WglFT9OlN07949AJzddJs4wsNeYmJqRsDjRwQGPKFIidIYGRvrNCZtq17Inn9jmlO16t/Mn/I5hR5Z4WjVmHrFHHUdmqLkSek9VXUZyKfNQHK7Bw8egBC4umb/qapX4WFIKfnzt200rVSE58+ecvzIAfp1aELQ80BOHPmbWV+MICToRbbHlh0qe9py/3kkbYfPRBMdheHR5QzbfJZXUbG6Dk1R8qT0Jg4H4KoQYp8QYlfCQ5uB5TYPHjzA3MYRO8vsrTZ/7dJ52tf14tyJ/yhWuhxd+w7BxNSMuo2aM+eHTdg7OeP/4B6nj/6DmbkFh/f9wfF/DmRrjNpmoKdH/eKOrL4USYt+n3Pl1L/EXP6TE3cDdR2aouRJ6T1VNUWbQeQFDx48wNjWCWtTw7Rnzort3fXh6aOHlK1UlVoNmmBjZ0/hYiUZNHZS3AyWVtRu2AyADt370L5rL/T09Ni+cRWRr15Rrc4HaGJjMTDMnni17YMSzpy485zhnwzmxfWTnPh9OT83aECD4u10HZqi5DkivZcthBAFgKJSygNCCDNAX0oZqtXoskjlypXlmTNntLqNokWLEmLuweKV6/Gw0+5Rh5SSgZ1aEvT8GT/uO46+vn66l42OiiLoRSCxsbEM7daOr+Yto2zFKlqMNvtIKRFC8DwwgK7N6hBrYkXg3asY5/FrPIqiLUKIs8m1B0/vXVWfAtuBH+KfcgN2ZFl0uZxGo8HX1xdp7oCVlo84/ju4l4DHj5i+aDWLNvyWoaQBcd30HJ1diImKwsHJGXtHJyIjXmkp2uyVMPDSzt6RibMXEvboDsPHjNNxVIqS96T3GscQoBYQAiClvAU4aSuo3CYgIIDIyEg05vZYGqf37F/GhYYEM23sEFYsmIVjPhec3+FCvLtnIZZt/QMbO3sGdWnNqgWzszBS3avTsClF67VnxdJF/P3332w/48vyf27rOixFyRPS+ykXKaWMSvhGJ4QwANStufEePHgAgJVDPvT0tFNuxO/+XfK5urNk805c3QtkyTqFEBibmFK6fCVKlq9IwmnLvFIy5fOvZjD049O0+agLRQZ8j4GZFSVdrNRtuoryjtJ7xPGPEOILwFQI0Rj4Gfhde2HlLgmJw0FLo8bDQkMZ3KUV304eQ9GSZTC3zLrqu/r6+oye8g016jVi9hcjWDTjS/LKcJ3S+Z2YvXg1r0JeoH9kKQPqFGTUNm9eRsboOjRFydXSmzjGAwHAJWAAsAf4UltB5TYJiaNggfxaWb+5pSWDxn7FRz37a2X9EHeUYWpugZmFhda2oQvVqlRi5JczOHf0EN57NuJiY8K5+3lzPIuiZJd0naqSUmqEEDuAHVJKVUnuLb6+vhiZmJLfNesv+2xeuZiXoaF8OmI8enraK9cuhGD4xOkIIfhj+4/o6enR4sPOWttedvqwW2/OnfiPH+bNoP2Xnpx/4ExddbpKUTIt1U8iEWeKEOIZcB24IYQIEEJMzp7wcocHDx5gauuMm23W3oYrpcT33h387t/JlusOQgg0Gg1//f4LB/fsyDOnrIQQTJi1AGdXDw58/wVHr9zTdUiKkqulOo5DCDGSuOZN/aWUd+OfKwQsA/ZKKb/LlijfkbbHcVStWhWfYPhh8y+4Z1HyiI2N5cblC5QqX5HY2NgM33b7LsJCQzE2MeHOretYWlnj4q6dU3DZ7fplb/p/1AJTTy9eXDuWra+pouRGmR3H0QPokpA0AKSUd4Du8dMU4o44Ys3syWdlkmXr/P3nTfTr0ISrF89l+wecuaUlGo2Gz/t3Y/bEUdm6bW0qUcaLYRO+JvTWKSZN/0bX4ShKrpVW4jCUUj57+8n46xxpjnQTQjQTQtwQQvgIIcYnM10IIRbFT78ohKiYgWXHCCGkEMIhrTi0KSIigidPnmDj5IKBftZdg2japiPjZ3xHybIVsmydGWFkbMz0xWv4at73REVG6iQGbejwSV88Kjbgm2mT+eeff3QdjqLkSml90kVlchpCCH1gKdAcKAV0EUKUemu25kDR+Ed/4k6BpbmsEMIDaAw8SCN+rfPz8wMgXxZWxV25YBbH/zlIm06f6HRMRZkKlYmNjaVnm/r8+ds2ncWRlYQQfDxyBtbOHnz00Uf4+qreHYqSUWkljvJCiJBkHqFA2TSWrQr4SCnvSCmjgK0kbTfbFtgg45wAbIQQLulY9jvgc3LAIMSEW3E9PLLmOkBUZCRH/97P5fOns2R978rOwYkiJUrjmM9F16FkmcZeBXH7aBIvw1/RoUMHIiIidB2SouQqqSYOKaW+lNIqmYellDKtU1VuwOtf5/zin0vPPCkuK4RoA/hLKS+ktnEhRH8hxBkhxBlt9qJOSByFPd99NLeUksjICFb/doABo7545/VlBX19faYtXEWpchX5Yf7MPFHXysrUkPGdG+LcajSnT59myJAheeYOMkXJDtobGADJnWN5+68zpXmSfT6+Ku9EIM3bgaWUK6SUlaWUlR0dtXfPfmLiKPjuRxxnT/zLh3W9uHn1IsYmpu+8vqx09eI5Nv2wkDPH/9V1KFmikKMFbdu1pWbH/qxZs4Yuo6bhH5T7k6KiZAdtJg4/4PU+qu7Aw3TOk9LzhYGCwAUhxL34588JIXTWnfDBgwcYWtjhbPPuZUAcnFz4oHkbChUrmQWRZa3KNeqy7cApajVowovAJPdL5EpNS+cjqERbSlerx7ZF0xg878fEaY+DI1h22EeH0SlKzqXNxHEaKCqEKCiEMAI6A293DdwF9Ii/u6o6ECylfJTSslLKS1JKJymlp5TSk7gEU1FK+ViL+5Gq+/cfoGfpgK250Tutx+f6FZCSCbMWYmycdbf1ZiVXjwIsnTOVPu0b8TI0RNfhvDMbMyNqF3PiVe0hOLu68+eCz9n130UAtpx6wPJ/7qhTWIqSDK3VAJdSxgghhgL7AH1gjZTyihBiYPz05cTVvGoB+ADhQO/UltVWrO/i3v37mNg6Y/iOt+J+P+dr7t66zs+HzmJgoL3S7O+qXuOWmJlbYGKavS1yteXDiu5UL2iPUcNN9O3YlO6dOnD/6lm2n/XjZWQMAaGROGXh+BxFyQvS3QEwN9PWyHEpJWZm5jhWbcW2NcveaV1BzwPxvXebshWrZlF02iOl5PefNuGUz5Xq9RqmezmNRoOenh6P/H159uQx5haWFCpWgheBzwgMeIKpmTm29g6Ymeum0OLxfw4wul9XCnjVwuPjSdiYmzK2aXFV10p5b71TB0AleYGBgUREvMLWyfWd1rPiu1lcu3Q+VyQNgJjoaH5av4K9O39KcR4pJRuWLWDKqIH0ad+IJhUK8d3XEwD4ddNqBnzcnCWzvwLg0N5d9GhVj48+qEyj8p409irIs6eP0Wg0/DB/Jnt3/ERIkPYr2tao14hPP/+a++f/JeLoBtxtTbnxOFd0R1aUbJVzz4nkAgl3VDm5ZD5xvAoP46/ff0WIuA+u3MDQyIhFG37Bxs6BoOeB2NjZExsby/dzpnLiyN+MnDyTyjXqcmD3b4SGBJO/YBGatv2IyjXrAtCm0ydUqlEXWzt7AKrXbciMJWt5FR7G82dPefroITa29jx78phNKxYRGxPDlv3HsbS24auR/SlSogyVa9ShRFmvLK8Y3PvTgdy57cPB7etwyV8QY8P2Wbp+RckLVOJ4BwmJI5+bRxpzpszYxJQt+48TE53qQPwcx8bOga1rlvH9t18zduq3tO3cg4tnTuLglA8Dg7ghPmt+O4CBYdLhPh6ehfHwLJz4u6tHAVw9ko6DcXJx5e+LD/C9d4f8BYvw4vkzbl69xIE/fgPA3tGJsV/PpW7jFkgps2yU/ZQZ3/Aq0J/dP8ykQqwVdNJN2RdFyalU4ngH/x81nrnEEfQ8kH4dmjJ6yuxcc7QRFRnJ5lVL2LV1A08e+WNgYEhUVFwtqx9+/vONI4DkkkZGGRoZUahYCQDs7B3Zuv8ELwKfceq/Q/z39z6cXd2RUvJpx6Z4Fi5G9/7D8CxS7J22qa+vz9TvVjCwUyvOr53M2f4NqVRRJQ9FSaCucbwDPz8/9AwM8XDJ3DCS0JAgChQqgpOLdlrOZqUHd30IePIIA0ND9u38mfyFijJ1wQr2n79N2YpV+HXzGq02mnqdrb0DTdt+xLSFqyheuhxRkREULl6Kf/bvJiIiHICf1q/g8UO/TG/D3MKSeau2YGBqQfPmzbl7927aCynKe0LdVfUOunXrxq97D7F+70k87DJ2e6pGo+FFYAD2js5ZHldWWzxrMtvWLqdL38EMGTeF8LCXb9z5NHP8cE4fPcyP+45hamauszgjI15hZGzC44d+fNSgEgB1GjWn4yf9qFi9dqZOZU3buI8j3w3GxcmRo0eP4uSU9V0eFSWnUndVaYGfnx96FvbYZWLw35ljR2hXpzznTx3TQmTvLiQ4iItnTwJx13A+6tmfLn0HAyS5XXbYxGls3PMvenp6aDSabI81gbGJKUIIXNw8+PnQWbp9+hnep4/z2Sft8T59HCDDA/o6NKxO4e7T8Pf3p0WLFoSGqrusFEUljnfg6+uHvqU9ZkYZb7TkUbAw3fsPo3T5SlqILPNiY2PZsWUdnRpV5YshvYmKjOSjHp8yfOJ07ByS/7ZtYWnFq/AwujWvza5tG7I54uS5uHkwaOwkdvx7ga8XrsSrSg38H9yjc+Nq7Niyjuio9N2MUNLFCtdiXvT4cgHnzntTp3ELIvNQfxJFyQyVODJJSon/Q38s7ZwzfArk+bOnPPK9T/+REzAyNtZShBnne+82Az5uzpxJYyhUtAQL1v2c7vgcnPJRuWZdChQqquUoM8bYxJRGLdsjhOBlaAiW1jas/z6u43FkxCtioqNTXV4IQdvyrhwMc6dCt/FcOPkfPXv2JDY2NjvCV5QcSSWOTAoMDCQqMhKbTFyj+P2nzQzt3o6Hvve1EFnGxcTEIKXE0MiYoOeBTJm/nCWbd1KkROl0r0MIwfgZ31GoWEn++HmzFqPNvOKly7Fy+z5W/rIPQyMjNv6wiM5NarBv58+pnmIr72HDD90r8+WIQbg3/ZRt27YxcOBAnZ6WUxRdUokjk/z9/QFwcM54g6NOvQcwb/VW3PJ7ZnFUGed3/y4DO7Xgr99/IZ+rO1v/OkmTNh0zPSbi1x/X8s2k0TkmKb5NCIGDU9xdcGUrVsXcwoKpowfRs3V9jh7an+Jypkb6OFsZk6/OxwweOZZVq1YxfPhwVQRReS+pxJFJCS1jnZwzNmr80rlT7P5lKxWr1dJGWBm24ruZ+N69nVi08F0LLHbtO5h1uw7h4p4/x38jr1anAWt3/s3UBSuIjHiVeAE9MjL5joBCCEq7WlH9o8GMHj2aJUuWUKl1D16EqWseyvtFJY5MSkgcrm4ZG4Pxz/7drFo4Gz29jF9QzypRkZHs2LoeKSUjJ81i3e+Hqdu4RZas29jEFI8Chfi8fzfWL/suS9apTXp6ejRu9SE/7j1Gv2GfExUZSc/W9fl28hiCngcmmb9EPiuO+gQy+5s5FKnfgfO7N1Gr4wCiYnJ2klSUrKQSRyb5+/sjhB7ubhk7VTV0/FTW7TyUOKr66aOHvAoP00aIyQoMeMJnn7RnzpejOX/yKLb2Dri8Q8mU5BgZG2Nr74ilpXWWrlebDAwNMTYxJSYmmmp1PmD3L1sTG1a9fuRU1s2aoz7PqD3nECXaDaP1x59wbe962vYbpavQFSXbqcSRSX5+fhhb2eFgmf6Bf3/98SvDe3agX4fGXD5/GoDtG1fRukYp1iz+ligt3+Z59eI5+rRrxK1rl5m+aDUVq9fW2ra+mL2Q9t16c+SvPbnqOoCZuQUjJ83ktyPeFCxanFP/HaZXmwZcOhf3ftmZG/F9t4oM+6Aoo5qWYNz0edRt0YG96xfy9fTp/HX1MQeu6qyvmKJkC5U4MsnPzw+DDHT+k1Ly6+Y1nD76D04ubhgZxd3m2rBFW6rVaciqhd8wpGsbAgOeaC3mLau/R9/AgB9+2sMHLdpqbTsJDu7+jfGDenDm2BGtbyur2do7JP4cGhzEgI+bM3viSEKCXmCgr0dBB3NMDPXR09Njxnff416lCV9NmkS3gaP54Z/bOoxcUbRPlRzJpNKlS3M/xoY/dv2GsUHa1yvWLp3Hyu9m0b5rb0Z9NRt9/TeXObT3d6aNHYKzqxtLN+9McbBdZvx3cC81GzThVXgYUZGRb3woalNsbCxH/95HnUbNs6xyrS6Eh71kzeJv2bZ2OdXqfsDclVuSzHPydgCTxg7n5cX92NXoyN2Dm7Eyfbd2woqia6rkSBbz9fPD2MYhXUnj1H+HWfndLGztHRn2xbQkSQOgQbPWzFu9lScP/dmy+vssiVFKyZLZX/H5gO78+ds2zC0ssy1pQFyV2bqNW7Bu6TwWz5qcbdvNambmFgwdP5W1O/9m8Ni45lM/rlrCg7s+ifNULujAoC+/oV3X3jw/vp1P+g3OVafoFCUjVFn1TAgNDSU0JISC9umriiuEwLNIMUqU8cLYJOX+1RWq1mTZ1j8yNPAuLdHR0XT8pB/N2n2cZevMqBeBzwgNCU5sHZtbJbwvLwKfsX7ZAlbMn8UPP+2heJny6OsJ2lf0QFaYw6OXsez6cTWfmunx6YSZPA2NolmZfO/cl15Rcgp1qioTrl+/TsmSJanYczJLJg1Lc/6w0FCMTU0zNEbiivdZnj97Sp1GzTMcX2RkBKsXzqH30NGJhf90eaooNjYWfX19bly+QOESpd95rEhOEBjwhF83r6HvsHGEBL3g2dPHiYnF/0U4g0d+TsB/W8lXoRGFO46lU/WCfPZBzirHoihpUaeqslDCGA7HfKnfivvI7wFTRg1k6TdT+LCeV7oL6wEs+eYrZn0xkuAXzzMUW2RkBOMH9mDzysWcP3UMPT09nV9fSEgafdo3yjFFEN+VvaMzn46YgJ6eHisXzKJ3u7gbHKKjonC1MaXPiAl0GjyOx+cPEL57FisOXlP9y5U8QyWOTEgoN+KcRq/xDcsXcOjPXXhVrUmnXgMwNEr/xdIxU+YQGhKUoWsDkRGvGDfgE079d4gJMxdQs37jdC+rbcVKl2Ps19/StO1Hug4ly306YgINW7RjzeJv6d2uIdcunaetlxvDR41lwswFXDhxhKDtk/hy63Fdh6ooWUIljkxIOOJwS6Vz39NHD9n9yxYatmyLV5UadPv0swxto3DxUnT79DP2/LqV00f/Sdcyu3/Zwumjh/li9iJafdQtQ9vTNiEE7br04sKZE8yfOl7X4WQpGzt7psxfzrcrfyQ0OIgr3mcTp7X+uDszlq7jyb0b/D6jH8cv3NBhpIqSNVTiyAQ/Pz+MLaxxsrVMcZ4dW9cTGxODoaERHzeskqnR4b2HjsbVw5OFM75MtYx3VGQkkZERtO/amx9++pOWHbpkeFvZxef6Vc6e+I/QkGBdh5LlajVowo97j9Ghe19eBD5jcJfWXPE+S73GLViw9idk2HOaN6rPlStXdB2qorwTlTgywd/fH0MrBxwsk+9VERMdzR8/b6J6vYZ8MmA4E+cszlRLVWNjEwZ/PpnAgCc8uOOT7DwxMTFMGTWAMf26oNFoKFMhyXWsHKVr38Gs33WI6OiobC21kl3MLS3R09Pj6SN/Hvk9YHS/zoSHvaRCtVrMWvMrYRFR1KhZiwMHDug6VEXJNJU4MsHPzw9paoebjWmy029dv0xIUBC1P2iKv+99GjRrk+ltNWjWmu1/n6Vg0eLJTj97/F8O7/uD2g2bJTs+JKcxMDQkPOwl3ZrVYs3ib3UdjtYUL1OejXv+ZfayDZiZW3D0733YmxvR79vNxJjY0qRpM+YsXKbrMBUlU1TiyIQHvr6Y2jphYpj8B3XJshXYefQSvvfu8Hn/bkS8Cs/0toQQmFtaEvziOWeO/790h5SS8LCXiaXBO/UakOltZDcrG1t6DR5Fyw5ddR2KVllYWuFVpQaxsbEs/WYqfT9sTOzNI/zyxz4KlqvKuBGD6T9sTI4vP68ob1OJI4MiIyMJfPYMB+fkB/+FhYYSFRmJta0dA0ZPZOmPu7CwtHrn7c6bOo4vBvdKvDaw8YdF9GnXkBeBzyheutw7rz+7deo9EANDQxbPmpznPzj19fVZtvUPGjRrzcrvZjGqbyemz55HlaYdWbl4Hh916sykX87h+zzzXzAUJTupxJFBDx8+BCCfa/K34v60/gc+rFeBy+dP8+umNRQskvwppozq3n8YL0ND2LZ2OXt3/MTyudMoUdYLa1u7LFm/Llw6d5qdWzdw//YtXYeidda2dny9YCXTF63mZUgwZuYWzF+0lJJtBvDr9p+ZN7wbPx25qOswFSVdVOLIoIQxHPndk+9hcWjv77gXKMjFsydZNncakqwZmV+sVFnqNWnJllVLmTl+OBWr12bi7MW5uoRHs3Yfse3ASTwKFiYk6IWuw8kWH7Roy5b9J3BycWXfzp8RvudpM+hLYgMfMLl3a06ePMnIbd4c83mm61AVJUW591NHRxLGcBQtmD/JtAd3ffC5foUGzdvQtd9Qdvx7IUtOUyXo89lYXr0Kx8LKmplL12VoQGFOJITAzsGJET078NXIAe9NUcCEmxiMjI154HOdCP/rzNu4ixgMqFu3Lps3rufQjQAdR6koKcv9RYOyma+vLwClihZMMu3vP3cBUKxkGWaOH06vIVnXFe5F4DPsnZyp37QVp4/+g2F8B8HcTghB8w87Y2pqpvPSKNmtUcv2lKtYDUMjI6xt7bAtVRPjF3d5sHsBa8L8+LzJj3nmfVbyFnXEkUG37/siDI3J55D02sKhP3+nTIUqhIYE88/+PzAwyJo/+sjICD4f0J3Puren/8gJLNm0I1PjQnKqlh26UOuDJqz4bhbPnr5f3fOcXFyxtXfgns9NAs/+ie+dGxT2qsmdw9tp3KQJz56pU1ZKzqMSRwb53PPF1MYxybfjyIhX2Dk40qR1B+o2bsHuUzdwSqOWVXrdvHKJOzev0W/4ODyLFKd4mfJIKYmJjs6S9ecETx8/ZMvq7/n3wJ+6DkUnChUrwfClOzBxKsRt72NYuBbm2NFjVKpUiZMnT+o6PEV5g0ocGfTw4UMsbB2TPG9sYsp3a3+ifrNWfNqxGVcvnE1m6YwLD3tJ2YpV+PnQGRo0aw3Aq/Aw+n7YmB+zqOFTTuDhWZhtf52kfdfeGa4InFc0rlaO8Uu3MHzidFwLlWTU4q0A1K5dm0WLFr0314CUnE+riUMI0UwIcUMI4SOESFLZTsRZFD/9ohCiYlrLCiG+FUJcj5//NyGEjTb34W0BTx5hbZ+0revtm9eIiYkhMOApGk0s5hYp17FKr9NH/6Fjg8pcv+yNnf3/k5WpmTk2tvZsXbPsnQYX5jROLq78uGoJXZrW5EXg+3eKxsrUkMalXOjUeyADJs3jAc64u3tgb2/P8OHD+eijjwgOzns1vpTcR2uJQwihDywFmgOlgC5CiFJvzdYcKBr/6A8sS8eyfwFlpJTlgJvABG3tw9uklAQ9e4q905uD/8JehtK77QesXjSHoiXLsPrXvyhc/O1dzRi/+3f5clhfbO0d8PAsnGR6z8EjCXr+jJ3bNr7TdnKaanUb0rJDZ8zM8841nMwo627N09BIAhwrEBQcTJ8+fdixYweVKlXi/Pnzug5Pec9p84ijKuAjpbwjpYwCtgJt35qnLbBBxjkB2AghXFJbVkq5X0oZE7/8CcBdi/vwhqCgIGKiInHO9+a1izPH/yUmOpqyFarQtlYZ9u746Z22Ex72knEDP0EIwTfLNyZ79FK+cnUqVK3J5pWLiYqMfKft5SSFi5VkyLgpHDt8gIO7d+g6HJ2xMDbg63ZlGT58OFXHbyGwYh9WrlmDv78/1atXp8Owqaw8clvXYSrvKW0mDjfA97Xf/eKfS8886VkWoA+Q7NVUIUR/IcQZIcSZgICsuSc+YdS461ujxk/8cwAzcwsKFStB9XqNcMuf9FbdjLh94xrPnjzi6wUrcS+Q8rp6DRnNsyeP2fPrlnfaXk6j0Wj4af0P7Ny24b0+ry+EoGJ+W2Z3q42xgR4HbjxHT0+PWI2GXxdPYfSn3bl02zftFSlKFtNm4kjupvy3PwVSmifNZYUQE4EYYHNyG5dSrpBSVpZSVnZ0THoxOzMSEoeH+/9zmJSSE/8cpEqtephZWDJx9iLKVqyS6W1EvAqnbMUq/PLPearWrp/qvJVr1qVBszaYmVtkens5kZ6eHjOXrGP+6m08e/I41V4k7wMhBL1qenLVuBRT1u7BoVBZAF7eOkm1ShVYu20HMbF5u96XkrNoM3H4Aa/X5XAHHqZznlSXFUL0BFoB3WQ2fiX184srN1Iw//8Tx91bN3jyyJ/Kterx0QeV+WH+zEyv/+rFc3xYryLep4+na8S5EIIZS9bQpE3HTG8zp7K1d+Dp44d0bVaTrWvyzt1jmWVjZsTwhsW4GGqC18D5fPblTIZNmoWRmRV9OreneKPOPHmuepor2UObI8dPA0WFEAUBf6Az8HYd7V3AUCHEVqAaECylfCSECEhpWSFEM2AcUE9Kma23FPncizst4JzPJfG5iIhwqtZuQOUadYgcNJKS5Spkat1SShZMm4ixiUmGCyM+ffSQP7ZvpufgUbmiJ0d6ubjnp2u/Ie/UzyQvKZ7PkuL5Eq53xR11VKleix6t6nHn8M8UKnmSI3t3UqmCl85iVN4PWjviiL+APRTYB1wDfpJSXhFCDBRCDIyfbQ9wB/ABVgKDU1s2fpklgCXwlxDCWwixXFv78LY7D3wxMLXAxNQs8blS5SqyYN3PCCHo8ElfKlStmal1CyGY9f065q78McMVby+eO8mqhd/wz/7dmdp2TiWEoPfQuH4Vk4b1Jeyl+kb9Ns8ixRn79Vzc8nsSHR5CtapVmTZtGisO3yIi+v0+xadoj3gfLj5WrlxZnjlz5p3XU7txSy5ducref+PWFRMTw80rFylWuhzdW9TGo0Ahvl35Y4bXu3nlYkxMzPiwe59M1WuKjY2lW7OaGJuYsm7XoTxX8+n8yaNMHNqHuau2UKp8xbQXeA9pNBoOnPdh3vAuhD6+j4GdO98uW8uIjxvpOjQlFxNCnJVSJulHrUaOZ8DjR4+wtPv/4L+bVy/Sr0MTDu3dxYiJM+jSd3CG13n+1DGWfTuNi+cyX1ZCX1+fHoNGcuvaZY4d/ivT68mpKlSrxS//nKNIidKcP3VM1+HkSHp6enzgVQSHhv0wMLMi5rkfIzs3ZdKkScTExKS9AkXJAJU4MiAw4PEbo8bPnfgPABtbO/K556di9doZWl9I0AumjhqIq0cBPp82752OFJq07oCLe37WLZ2fJ29hNTUzZ9ncaYzs9REBTx7pOpwcyUBfjzatWtFi+i+07dwTpIbp06dTvXp1Ll26pOvwlDxElVVPJ41GQ0hgAOVeaxl77sRRPAsXY+ua5fjdv8O2A6cytM5lc6fxPDCAH37a884lSgwMDene/zPmTRnH3Vs3KFSsxDutLyfqOWgklWvUwdHZhZiYGAwM1H/ft31UyR2NdMOgTQXMyjXh5o3r+Oz6nvLly9O/f3/mfDsXK8u8dfu2kv3UNY50evr0Kc7Oznz82WRGDB9GTHQ0TSsVoXn7TnTr/xkBjx9SrlK1DK0z+MVzLpw5Qd3GLd4ptgRRkZE8fuhL/oJFsmR9OdXC6RN55O/LzKXrcnUHRG2LitHwn88zdu/dz/Hvx4AmBn0rRybMmMe0oZ/oOjwlF1DXON5R4qhxt7hR49cuefMqPAxXjwLcuHyBkmXTfxvu3Vs32LpmGZbWNlmWNCCuo1z+gkWIjY0l6Hlglq03p8nnnh8X9/x58pRcVjIy0OODEk7MG9GdjQfO0rbXZzjbWjH9sx5UqFw1sQ2yomSUShzplJA43N3iBv9JqaFanQ/wu3eHr8cMTveHWGRkBJNHfMqG5Qu10mdbSsnQbm2Z9vmQLF93TtGp1wCGT5zOP/v+4PC+P3QdTq5QOL8b4778ih/3HMHZsxjeZ0/jkT8/w0d/ri6eKxmmEkc6JXw784hPHOUqVeO7tT8x6qvZrPxlH0bGxulaj9RoKFWuApO+XYKNnX2WxymEoFrdDzh++AA3Ll/I8vXnFBqNhq1rl/Pbj+vUkUcGGBkbs33fv3QdPBaEHovmf0v+/Pk5cOCArkNTchGVONLpnq8fAPny5UOj0XD9sjf37/iw8YeF2KYzAdy5eR0DA0MmzFpIjXrau7++4yf9sLC0Yu3SeVrbhq7p6ekx54dNfLtiM/d8bqo7rTJAX1+foaPG8ePf57EoUhmNhMaNG1O1alV195WSLipxpNP9B/4YmNtgaGTE/du36NOuET+uWsKqhd8QFRWV5vLPnj5mSLc2zJk0RuuxWlha0bnPII78tYdL505rfXu6knDENqrPx8yaMEK3weRCBdxc+GTKcip/vhHX2h05ffo05cqVo3///oSEhOg6PCUHU4kjnfz8/TGxjvugung2brBe9/6fsefUDfK5pt4SRErJrAkjiHj1iq6fDtV6rACd+wzCzsGJ7+dMzdOncoyMjflyzhLGfv0tYS9D81Qf9uzQprwbZfM7MuDzqZTu8TVFihRh5cqVFCpUiOad+xAa/oqXkTEcuZk1rQmUvEEljnR69OgR5jZx5dkvnTuFta0dqxZ8wyP/B2kuu3PbBo7/c4Ahn3+FZ+Gi2g4VADNzC0ZMmkG7Lj2zZXu6VKlGHewcnRjStQ3ffjVW1+HkKtamhrTxcqNRSScKVKrPrK2HOH36NLF6huzdthYbaxuKNe9L/w2nOXv//ewFrySlEkc6BTx5hKVdXOK4ePYkhYqW4OyJ/3gVFpbmsn737lClVn0+7N5H22G+oVHL9jRt+1Geq12VHGNjEz5o0Y4GTVuj0Wjy9FGWNggh+LiyBzP3XCPUMj9u/VfQqF0X0MTy6NAGHizsysCZq4hVfT8UVOJIl5iYGIKeP8PawZnnz57id/8uNRs0Ztexy5SvUiPF5WJjY4mOimLo+KnMXfmjTgarRbwK5+sxg/n9p03Zvu3s1mPgcKrUrs/MCcPz/Ck6bSjmbEnvmp4M3HiW2sVd+HruYvaevcUHLdpiZGjAmRXjKFC0BOvWqTvZ3ncqcaTD06dPkRoNtg7OhIYEU7lmPf7Zv4d9O7enmgy2rF7Kpx2bEhoSjKGRUTZG/H/GJqb4P7jHD/Nn8jI071/w1NPTw8TEFFMzc12HkitV9rRj9ofl6FQlPxB3o8X0RavZdfQSfcZOw//ebXr37k3x4sXZv39/YgLRaCTf7rtO8Ct1jel9oBJHOiSM4XB0zkeBQkWZuXQtRsbGxESnfDfVi8BnrF70La4enunq5qctQghGTp7Ji8AAVi38RmdxZBchBKOnfEOfz8ayc9sGln07TX07zqB81iaYGL7ZEMzQyIh+AwYxdu1fWBauyP3HgTRt2pR8+fLxzTff8Nt5P5b/c4dZe67pKGolO6nEkQ4Jo8bzubhw7/YtTEzNWLJpB60/7p7iMrb2Dixcv52xX3+r82sMJcp40a5LL7ZvWMmta5d1Gkt2SHi9b129jM+NK8SqkdFZpl2tcize8AuLdxzFqX4Pnj17xvjx4/m4ZjEK393J/kv+nLiTd8vdKHFU4kiHhMTh5OBAz1b1+KRFHeZ+9XmK8/+8YSVXL5yjXKVq2No7ZFeYqRow6gssrW2Y+9XnaDTvxwXOMVPnMGvpem5evcT8qeNVaY0sIISgmLMlZQs48unwMTSZtRuPKo2RMdH8tXUFtxf3pOPgLzh/466uQ1W0SCWOdHj48CEIPSKCnhIdHUWBwkVxeK28+usunTvFwukT2bltQzZHmTorG1uGfTGdEmW93puxDkIIjIyNOXviP44d/ouXIcG6DilPaVo6H0XcnRg2bTF7zt6h99AxFCpcBP99K6hYohDFy1bA+7I6dZUXqbLq6dC3b182/ryT7n0HsnbBDP44cQ07B8ck84WHvaRn6wZoNLFs+P0fzC3frceGknXCQkOJjo5izqQxjJw0E8d8LroOKc/a//svzJ8xmZBnTwBBkyaNadasGcOGDUNfXz/N5ZWcQ5VVfwf+/g/RM7fD58p5LKysGdy1dbKnPZbPnc5D33t8OWdJjk4ae37dyucDuhMbG6vrULKNuaUlD+764H36OE8e+asL5lrUpHUH9p64wraDpyncqAt/H/qHUaNGYWxiSvOWrTh00pu7z9Ie/6TkXCpxpIOfvz+GVvZcu3iOQsVKUL9p62S7zzVu3YEh476iQtWaOogy/aSU/HdwLxt/WKjrULJVuUrV+OXwWUqVr8jEob3ZsXW9rkPK0zwKFOTbb+ZQcNh67Cs1R+gbsnfPbj6oXoGKDVqxefvO9+rLS16iEkc6PHr0EDNLWzw8C9OwRTsGjp74xvRX4WFcvXiOshWr0LVf9tSiehctPuxM49YfsnrhN1w4c0LX4WQrUzNzIiNeERkZQWx0DBqNhvCwl7oOK8/KZ23CzM61+GXzeo5cecC0xWup26g5r+6epftH7TAxMaFp06Zcu6auheQm6hpHGiIjIzExMaFwsz50bViZmROGsWnPvxQo9P+aUwumT+SXjavYduAUrh4FsipsrQoLDaVnmwZER0WyZscB7B2ddR1SttJoNAgh2LJ6Kb/9uJ4fft6DnX3S61aKdoSFhdF76Ggen91PTHgoIPH09KRixYpMnTqVMmXK6DpEBXWNI9MeP34MgJm5OUVKlKLnoJF4eBZOnB4VGYn3qWO07dwz1yQNiDvnP3PpWkJDgpk2NucfJWU1PT09hBCUKl+JGvUbYWvnwLkT/xGdjhL5yrszNzfnx9Xfs+XIFdrN/Z0KHw3nxYsgfv31V8qWLYudnR19+vRJvBVeyVnUEUcajh8/Ts2aNbFxK4yjtRlrdxxMLDMipUQIQXRUFLGxMZiYmmVl2Nni2OG/cHZ1p3CxkroORaf8H9yjc+NqdB8wnAGjvtB1OO8VjUay6O9bVPa0o+ira4ydPJ271y4SGxONEALXAgUpX6Y0rT8ZSOMGdSnsaKHrkN8b6ogjkxK+8YQ+9eXW1Uvs37U9cdr3c6Yyb8o49A0McmXSAKhZvzGFi5UkKjKS3b9seW/vNnLL78mcFZvp3HsgZ4//y6wJw9+L2l45gZ6eoFu1/Gw6cR/nsnUw7zCTurP2UXrIcvI37sXTJ0/Z88fvDOrUkjJlylK8VjMWLFzEq1evdB36e0sljjQk1KmKjY6ibeeeVK3TAIDL58+wZfX3xMRE66TqbVbb8+sWZoz7jJULZr+3yaNGvUZY29rhc/0K3qePY2RkzI0rF4l4Fa7r0PI8R0sTqha0o+eaU3So5M6UNmUY2akxm5fM4Z9L91i44Rd6DBxByYLu3Dy2j5EjhmNmZoZDPlf69R+Ir6+vrnfhvaJOVaVh/PjxfDPnW5Aalm/bTblK1YiMjKBXmwZEhIezac9/OXrMRnpJKZn9xQh+/3kznXoP5LMJX+eJhJhZUZGRAHRsUJHylaszbdFqNBrNe/2aaFvgy0h2XXhIzxqe6OmlXN/tzt37jJkyixfXTxD53B/iP8PMrO0QxhZ4eNXh5+9n8DjShOBX0XSolHqHTiVlKZ2qSjoYQXnDw4cPEQZGyOgIjvy1h3KVqvHb5rXcv32L+Wt+yhNJA+LKc4yb8R3GpmZsW7ucF4HPmDBrAcbGJroOTSeMjI0B+HrhKoyMjQl4/IgBnVowfsZ3VK1dX7fB5VH2Fsb0rlUwzfkKFSzAogXzCXkVTYl8lvx95CjHTpzgzN6feOZ3l+v7N1O2yGYMrRwQhqZ81LIhXTt9RJMmTfAPjsTZKmn1XyVj1BFHGho2asShY2eoUL4sE2YuwC2/J9FRUZw4cpA6jZpncaS6J6Vk/bLvWDF/Jn2GjqHfiPG6DilHuHf7FgumfcG46fN48sif/bu203/kF9jY2es6NOU1ERERXDx7grNnz3Fo51b87t9JnCaEwLRgBVwLlaBTreL06t6VIkWKJE4PCI0kNCKaQurie6KUjjhU4khDsRIleRiux7y531KijBebVy7m414DMLfIG0caKTlx5G/KV66GqZk5j/x9yefqrvPy8DnFbz+uZf3337HtwElO/HuI4OeBtP64u3p9cqAHASFMWLQBfb+zPLp7C1sjDX53biZOd3Nzp3iJ4lx+8JxYx6JYFqvCgmGdqVHYgQPXntKynAvWpoY63APdUokjk4nDwtKKyFioUrUaJcqWZ+2SecxfvY3q9RpmcZQ5U8CTR3RtWpMqterx2RfTcHHz0HVIOUJMdDQGhoZ83r8bzwMDWPXLfratXY6TixsNmrXWdXjKa/ZceoSpkT7VCtphZmRAwJNHbN3yI0dPniY8SkP4g8uEv3j6xjImzoVwLlyaoKAXtKpXnc8HfkKxEqWQCEyN4k5zvYqKJTQyGifLuNO5UTEaDPREqtdnchuVODKROMLCwrCwsEDoG7BuxwHOnfiPe7dvMm76fC1EmTPFxMTw46qlrFn8LUjJR73607XvkBzTZ0TXpJSEBL3AysaWT1rWpVS5inwxeyFfjxlM9bof0KRNx8TxPkrOc+7+C84+eE6bYpac+ucvzp44wrXLlzCztOTx/bsEBz1PnFfP0AQpBMZmFuhbOqJn44p5/lJ0bteCBhVK8OXv14mI1lCvuCPfd6uIoX7uv5FCJY5MJA4fHx+KFi2KvpEJu/67gO17fD77yUN/fpg/g307f8bQyJixX39Lyw5ddB1WjqLRaHgVFoa+gT5DurWlcasPad+1F+3reDFg9ETadPqEs8f/pWTZCnnmpoq87v6dW/z11wGO+QTgLIK5cHg3oSFByXaV1NPXxzGfKyb5y2FmpI+ZvsTMwZUezWtSrGQpbkbacOLuCwrYm9OibD7KuFljbPD/i/RSSm4+eUkRJwv0c8hRi0ocmUgcR44coV69egg9fewdHVm/6/B7/0373u1bbFm1hA+79aF4mfIc+WsPd2/doEb9RhQpUVrdrvoaKSXBL56zcsFsGrZoi52jM12b1mDiN4upVqcBsyeOpOegURQrVYb7d3woUKgIxiamug5bSYfIiAiuXfbm+dMnvHwZyr8H93HnxhWklMTExPDsyaOkC+npYWXjwKuwUISRKcYFvDCzdUQT9JAShTx4gSXB+taYOxeib4saDG1aOsW7v569jORFWBRFnCy0ejSrk8QhhGgGLAT0gVVSytlvTRfx01sA4UAvKeW51JYVQtgB2wBP4B7wsZTyRWpxZDZxLF++nEGDBiGEHo1atWPqdysyvI68bsH0ify07gcAbGztqVi9NiXLefFBi3bqeshbIl6Fc+HMSQoVK0Hwi+dMHT2IMVPnADCocyvmrd6Kg1M+vvlyFKOnfIOdvROH9/3OB83bYmZhQWDAE1zc8mNoZKTjPVHSEvHqFbdvXuXOzetcvHyV8KAAPDzyc+/2TbxPHyciPBwzS0uCXzxHptTKWeghBFha22Fo5wpR4cRGhqNvZEq4NMTIzArnMjWoUMAO+fIZATEmuBQtT1lPJxoVs8PJyREnJydMTDJ/S322Jw4hhD5wE2gM+AGngS5SyquvzdMC+Iy4xFENWCilrJbaskKIOcBzKeVsIcR4wFZKOS61WDKbOFq1asXu3bsxNbfgl0Nn1a2XKQh48ogzx/7lzLF/OHfyKE8e+rF0804qVKvFltXf8+dv28jn5oGLe35s7R2wsrahfOXqFC5eirCXoQQ8foSRsXHcwyjuX0Mj4/emW9yLwGd4nz5Ohao1eeT/gOVzpzPqq9k8uOPDuIGfsOrX/bwIfMbYT7vyw09/8io8jGljBzNv9TaiIiNZu2QuIybNIDYmlj9/20bHHv2IiozkxJG/+aB5G2JjY7h28TwVqtVCo9HwyPc+nkWKERMTw/NnT3F2cUciCQsNwcrGDoCY6CiMjE3UEaSWxcTE8NDvPv737/LY35fHD/0wM7dACD2OHD6I3wNfDM0ssLOx5tF9H8JDXqDRaBIHPabH7t27adGiRabi08UAwKqAj5TyTnwAW4G2wNXX5mkLbJBx2euEEMJGCOFC3NFESsu2BerHL78eOAykmjgy6+rVuFA/HDAOjbElz8NU5dTk6FvYU61JO6o1aQdA8ItAzMwteB4WhbGVHXbOrvg9uM/ZE//xKr73xadjp2LrXoQThw8xfUSfJOusWq8xkxet52VIML2bVkbo6SGEXlxVWz2BuaU1K38/CsCkgV3wu3srbh4ExB+6T1myEY9CRdmyfD4Hdv2UZBud+4+gcbvOnD16mO9nJB2vUqFGPYZO+oaw0BCGd26aZLqZuQWLfvorbltDP8Hvrk+Seb5avCEuhh++4+/ff04yvdOnw2nUthOx+sb07fD/bYzq2wUpJQ3bdsIhf3GihTHW9o5MHjWA2JgYwsNf8fnAHhgYGGBuaUXwqxjmjBuEz9WL/LlzOzEx0QQ9e8r65Qtp32MAK+dMpsXHPTn2958EPXuKQz63xHm6D/mc/IWKMXN0vzeeT/g5ODCAxT8f4O6Nq8yb+BkO+dyIjYkmKDAAe2dXDAwMCAl6zrfrd7Fg8gh8rl3C3tmF2JgYggMDsHN2oVPfYWz54TtqNmzB0YO7CXn+DDunfMTExBDyPJDuQ8ZSsHhpFkwagaGxceLzdk7OxMTEEBEexvzNu7l/6wbfjh+c+HzCPPoGhsjYWKZ8v4n5Xw7nzrXLSebp1G84P69ZQpW6jTh6YPcb60+IoVCJMiyYNAIDI6O4518EYucYN09keBjzf9zD/Vs3mDN+UOLzCfMYGBgipWTK0o3MnziM29cvJ5mnU7/h/LLueyrXbsh/B/54Y/0hLwLpPngsroVLsfCrUegbGKAvYwjyu4VepDOGBvqYmJoxf/Nu7t26zpxxg7C1dyQ6KprQkBeYmVtgZGyCRhOLuXsJHt/wxtHOmuLFiyf5f/eutHnE0RFoJqXsF//7J0A1KeXQ1+b5A5gtpfwv/veDxCUBz5SWFUIESSltXlvHCymlbTLb7w/0j/+1OHAjM/uhZ2ZdUt/CLk/cyK15FaKnZ2qVwnFx7qL2JWdS+5JzCD19Q82r0ICY4CcPAAfgWSZWU0BKmaRRjTaPOJK7YvN2lkppnvQsmyop5QrgnS9KCCHOxIYFJTlUy42EEGdiQp6pfclh1L7kTHltX5I75ZRZ2jyB6Qe8fnXUHXi7K0tK86S27JP401nE//vmyB1FURRFq7SZOE4DRYUQBYUQRkBnYNdb8+wCeog41YFgKeWjNJbdBfSM/7knsFOL+6AoiqK8RWunqqSUMUKIocA+4m6pXSOlvCKEGBg/fTmwh7g7qnyIux23d2rLxq96NvCTEKIv8AD4SFv7EC8v3YOr9iVnUvuSM6l9ScF7MQBQURRFyTrqJm1FURQlQ1TiUBRFUTJEJY5UCCGaCSFuCCF84kep5ypCiHtCiEtCCG8hxJn45+yEEH8JIW7F/5tkDExOIIRYI4R4KoS4/NpzKcYuhJgQ/z7dEEIkHa2nIynsxxQhhH/8++IdX0EhYVqO3A8AIYSHEOKQEOKaEOKKEGJ4/PO58X1JaV9y3XsjhDARQpwSQlyI35ep8c9r732RUqpHMg/iLsrfBgoBRsAFoJSu48rgPtwDHN56bg4wPv7n8cA3uo4zhdjrAhWBy2nFDpSKf3+MgYLx75u+rvchlf2YAoxJZt4cux/x8bkAFeN/tiSuLFCpXPq+pLQvue69IW7cm0X8z4bASaC6Nt8XdcSRssSSKVLKKCCh7Elu15a4Ui3E/9tOd6GkTEp5BHj+1tMpxd4W2CqljJRS3iXuLr2q2RFnWlLYj5Tk2P0AkFI+kvFFSKWUocA1wI3c+b6ktC8pycn7IqWUL+N/NYx/SLT4vqjEkTI3wPe13/1I/T9WTiSB/UKIs/ElWACcZdxYGeL/ddJZdBmXUuy58b0aKoS4GH8qK+EUQq7ZDyGEJ1CBuG+3ufp9eWtfIBe+N0IIfSGEN3EDov+SUmr1fVGJI2XvXPYkB6glpawINAeGCCHq6jogLclt79UyoDDgBTwC5sU/nyv2QwhhAfwCjJBShqQ2azLP5aj9SWZfcuV7I6WMlVJ6EVdlo6oQokwqs7/zvqjEkbL0lEzJ0aSUD+P/fQr8RtzhaG4u2ZJS7LnqvZJSPon/Q9cAK/n/aYIcvx9CCEPiPmg3Syl/jX86V74vye1Lbn5vAKSUQcRVDG+GFt8XlThSlp6SKTmWEMJcCGGZ8DPQBLhM7i7ZklLsu4DOQghjIURBoChwSgfxpUvCH3O89sS9L5DD90MIIYDVwDUp5fzXJuW69yWlfcmN740QwlEIYRP/synQCLiONt8XXd8RkJMfxJVDuUncXQcTdR1PBmMvRNydExeAKwnxA/bAQeBW/L92uo41hfi3EHeqIJq4b0h9U4sdmBj/Pt0Amus6/jT2YyNwCbgY/0fsktP3Iz622sSd0rgIeMc/WuTS9yWlfcl17w1QDjgfH/NlYHL881p7X1TJEUVRFCVD1KkqRVEUJUNU4lAURVEyRCUORVEUJUNU4lAURVEyRCUORVEUJUNU4lAURVEyRCUO5b0khDim4+33EkIsif95oBCiRyrz1hdC1My+6BQldVrrOa4oOZmUMsd8EEspl6cxS33gJaDTZKcoCdQRh/JeEkK8jP8m/8drzy0RQvSK//meEGKmEOK4EOKMEKKiEGKfEOK2EGJg/Dz1hRBHhBC/CSGuCiGWCyFS/JsSQvQWQtwUQvwD1Hrt+SlCiDHxPw+LX9dFIcTW+MqtA4GR8Y2F6gghWgshTgohzgshDgghnF9bzxohxGEhxB0hxLDXttEjfp0XhBAb459zFEL8IoQ4Hf+ohaKkgzriUJSU+UopawghvgPWEfdhb0JcCZeEo4SqxDXGuQ/sBT4Etr+9ovgaSFOBSkAwcIi4MhFvGw8UlFJGCiFspJRBQojlwEsp5dz4ddkC1aWUUgjRD/gcGB2/fAmgAXHNiW4IIZYBxYgrMVFLSvlMCGEXP+9C4Dsp5X9CiPzAPqBkZl4o5f2iEoeipCyhqOUl4jqshQKhQoiIhKJywCkp5R0AIcQW4mogJUkcQDXgsJQyIH7ebcR9oL/tIrBZCLED2JFCXO7AtvhkZATcfW3abillJBAphHgKOAMfANullM8ApJQJjaUaAaXi6v0BYCWEsIzfT0VJkTpVpbzPYnjzb8DkremR8f9qXvs54feEL11vF3tLrfhbegrDtQSWEndkclYIkdyXu8XAEillWWDAW3G/HmdsfJwihW3rATWklF7xDzeVNJT0UIlDeZ/dJ+4bt7EQwhpomIl1VI0vva8HdAL+S2G+k0B9IYR9fB+Ij96eIX4dHlLKQ8SdfrIBLIBQ4k49JbAG/ON/7knaDgIfCyHs47eTcKpqPzD0te17pWNdiqISh/LeklJKX+An4k8Pkfw1h7QcB2YTV876LnENs5Lb2CNgSvz8B4BzycymD2wSQlyKj+U7GdeY53egfcLF8fj1/CyE+Bd4llaAUsorwAzgHyHEBSCh/8QwoHL8RfOrxF2EV5Q0qbLqynsn/pv3OSllgXdcT31gjJSyVVbEpSi5hTriUN4rQghX4r71z9V1LIqSW6kjDkXJYkKIk4DxW09/IqW8pIt4FCWrqcShKIqiZIg6VaUoiqJkiEociqIoSoaoxKEoiqJkiEociqIoSob8DwA9996xeSgoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_df = pd.DataFrame({\"r\": np.arange(0,300,1),\n",
    "                        \"superposition\": np.zeros(300)})\n",
    "for i in np.arange(0, optModel_components, 1):\n",
    "    unimodal_judi_model = em.JumpDistanceModel(diffusion_coefficient = optModel_diffusion_coefficients[i],\n",
    "                                               degrees_of_freedom = 4,\n",
    "                                               tau=0.02)\n",
    "    component_string = str(\" $\\\\omega$ = %.2f, D = %.2e [$\\\\mu$m^2/s]\" %(optModel_weights[i], optModel_diffusion_coefficients[i]))\n",
    "    model_df[component_string] = optModel_weights[i] * unimodal_judi_model.pdf(distance = model_df[\"r\"])\n",
    "    model_df[\"superposition\"] = model_df[\"superposition\"] + model_df[component_string]\n",
    "\n",
    "sns.kdeplot(data=data_df, x=\"jump_distance\", fill = True, bw_adjust = 0.1, clip = [0, 300])\n",
    "sns.lineplot(data=model_df.melt(id_vars=['r']), x=\"r\", y=\"value\", color=\"black\", style=\"variable\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}