notebooks/introduction/3_model_clogit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to choice-learn's modelling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"\"\n",
    "\n",
    "import sys\n",
    "\n",
    "sys.path.append(\"../../\")\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import tensorflow as tf\n",
    "\n",
    "# Enabling eager execution sometimes decreases fitting time\n",
    "tf.compat.v1.enable_eager_execution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- [Example 1: ConditionalLogit with Swissmetro](#example-1-swissmetro)\n",
    "    - [A few words on c-MNL formulation](#the-conditional-logit-model)\n",
    "    - [Instantiation and estimation with Choice-Learn](#conditional-logit-estimation-with-choice-learn)\n",
    "- [Example 2: ConditionalLogit with ModeCanda](#example-2-the-modecanada-dataset)\n",
    "    - [Utility Formulation](#utility-formulation)\n",
    "    - [Model Specification](#model-formulation)\n",
    "    - [Lighter method for specification](#faster-specification)\n",
    "    - [Comparison with other implementations](#comparison-with-other-implementations-results)\n",
    "    - [Get utility estimation and probabilities](#estimate-utility--probabilities)\n",
    "    - [Using Gradient Descent Optimizers](#using-gradient-descent-optimizers)\n",
    "    \n",
    "\n",
    "For model customization and more explanation on ChoiceModel and the endpoints, you can go [here](./custom_model.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: SwissMetro\n",
    "\n",
    "The choice-learn package offers a high level API to conceive and estimate discrete choice models. Several models are ready to be used, you can check the list [here](../../README.md). If you want to create your own model or another one that is not in the list, the lower level API can help you. Check the notebook [here](./4_model_customization.ipynb).\n",
    "\n",
    "Let's begin this tutorial with the estimation of a Conditional Logit Model on the SwissMetro dataset[3].\n",
    "It follows the specifications described in [PyLogit](https://github.com/timothyb0912/pylogit/blob/master/examples/notebooks/Main%20PyLogit%20Example.ipynb) and [Biogeme](https://github.com/mncosta/biogeme_tutorial/blob/master/6.1-MultinomialLogitAndProbitModels/01-biogeme-basics.ipynb).\n",
    "\n",
    "\n",
    "First, we download our data as a ChoiceDataset. See the [data management tutorial](./2_data_handling.ipynb) first if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from choice_learn.datasets import load_swissmetro\n",
    "swiss_dataset = load_swissmetro(preprocessing=\"tutorial\")\n",
    "print(swiss_dataset.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Conditional Logit model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conditional Logit [2] specifies a linear utility for each features of item $i$:\n",
    "$$\n",
    "U(i) = \\sum_{features} a(i) * feature(i)\n",
    "$$\n",
    "The probability to choose $i$ among the set of available alternatives $\\mathcal{A}$ is then:\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(i) = \\frac{e^{U(i)}}{\\sum_{j \\in \\mathcal{A}} e^{U(j)}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the SwissMetro dataset we are trying to predict a customer mean of transport among train, swissmetro and car from the features:\n",
    "- TT (transit time)\n",
    "- CO (cost)\n",
    "- HE (headway)\n",
    "- Survey (where the survey took place)\n",
    "- Luggage number\n",
    "- seats configuration in the swissmetro\n",
    "- first class or not\n",
    "\n",
    "An important step is to define the right utility function for the model to fit well the dataset. Let's take the following formulation defined in PyLogit:\n",
    "\n",
    "- $\n",
    "U(train) = \\beta_{train}^{inter} + \\beta^{tt}_{train/sm} \\cdot TT(train) + \\beta^{co}_{train} \\cdot CO(train) + \\beta^{he}_{train} \\cdot HE(train) + \\beta^{survey} \\cdot SV(train)\n",
    "$\n",
    "\n",
    "- $\n",
    "U(sm) = \\beta_{sm}^{inter} + \\beta^{tt}_{train/sm} \\cdot TT(sm) + \\beta^{co}_{sm} \\cdot CO(sm) + \\beta^{he}_{sm} \\cdot HE(sm) + \\beta^{survey} \\cdot SV(sm) + \\beta^{seat} \\cdot SEAT(sm) + \\beta^{first\\_class} \\cdot FC(sm)\n",
    "$\n",
    "\n",
    "- $\n",
    "U(car) = \\beta^{tt}_{car} \\cdot TT(car) + \\beta^{co}_{car} \\cdot CO(car) + \\beta^{luggage==1} \\cdot \\mathbb{1}_{Luggage==1} + \\beta^{luggage>1} \\cdot \\mathbb{1}_{Luggage>1}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we want to estimate:\n",
    "\n",
    "- one $\\beta^{tt}_{train/sm}$ **shared** for train and sm items and one $\\beta^{tt}_{car}$ for the car item. Indeed, one can argue that customers have the same sensibility toward travel time for all public transport and a different one for private ones.\n",
    "- one $\\beta^{co}$ coefficient for **each** item.\n",
    "- one $\\beta^{inter}$ and $\\beta^{he}$ for train and sm, and zeroed for the car alternative\n",
    "- one $\\beta^{survey}$, $\\beta^{seat}$, $\\beta^{first\\_class}$, $\\beta^{luggage==1}$ and $\\beta^{luggage>1}$ shared or not by different items\n",
    "\n",
    "To build a model, we need to specify for each weight $\\beta$:\n",
    "- the name of the feature it goes with:\n",
    "    - it must match the feature name in the ChoiceDataset\n",
    "    - \"intercept\" is the standardized name used for intercept, pay attention not to override it\n",
    "- items_indexes: the items concerned, as indexed in the ChoiceDataset\n",
    "- (optionally) a unique weight name\n",
    "\n",
    "\n",
    ">  **Attention**\n",
    "> \n",
    "> - *add_coefficients* is to be used to get one coefficient by given items_indexes\\\n",
    "> - *add_shared_coefficients* is to be used to get on coefficient that is used for utility of all given items_indexes\n",
    "\n",
    "Here is how to create a model following our defined utility functions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional Logit Estimation with Choice-Learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from choice_learn.models import ConditionalLogit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialization of the model\n",
    "swiss_model = ConditionalLogit(optimizer=\"lbfgs\")\n",
    "\n",
    "# Intercept for train & sm\n",
    "swiss_model.add_coefficients(feature_name=\"intercept\", items_indexes=[0, 1])\n",
    "# beta_he for train & sm\n",
    "swiss_model.add_coefficients(feature_name=\"headway\",\n",
    "                             items_indexes=[0, 1],\n",
    "                             coefficient_name=\"beta_he\")\n",
    "# beta_co for all items\n",
    "swiss_model.add_coefficients(feature_name=\"cost\",\n",
    "                             items_indexes=[0, 1, 2])\n",
    "# beta first_class for train\n",
    "swiss_model.add_coefficients(feature_name=\"regular_class\",\n",
    "                             items_indexes=[0])\n",
    "# beta seats for train\n",
    "swiss_model.add_coefficients(feature_name=\"seats\", items_indexes=[1])\n",
    "# betas luggage for car\n",
    "swiss_model.add_coefficients(feature_name=\"single_luggage_piece\",\n",
    "                             items_indexes=[2],\n",
    "                             coefficient_name=\"beta_luggage=1\")\n",
    "swiss_model.add_coefficients(feature_name=\"multiple_luggage_piece\",\n",
    "                             items_indexes=[2],\n",
    "                             coefficient_name=\"beta_luggage>1\")\n",
    "# beta TT only for car\n",
    "swiss_model.add_coefficients(feature_name=\"travel_time\",\n",
    "                             items_indexes=[2],\n",
    "                             coefficient_name=\"beta_tt_car\")\n",
    "\n",
    "# betas TT and HE shared by train and sm\n",
    "swiss_model.add_shared_coefficient(feature_name=\"travel_time\",\n",
    "                                   items_indexes=[0, 1])\n",
    "swiss_model.add_shared_coefficient(feature_name=\"train_survey\",\n",
    "                                   items_indexes=[0, 1],\n",
    "                                   coefficient_name=\"beta_survey\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimation of the model\n",
    "history = swiss_model.fit(swiss_dataset, get_report=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the model is estimated, we can look at the weights with the .trainable_weights argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'beta_intercept:0' shape=(1, 2) dtype=float32, numpy=array([[-1.2929296, -0.5025745]], dtype=float32)>,\n",
       " <tf.Variable 'beta_he:0' shape=(1, 2) dtype=float32, numpy=array([[-0.31433573, -0.37731653]], dtype=float32)>,\n",
       " <tf.Variable 'beta_cost:0' shape=(1, 3) dtype=float32, numpy=array([[-0.5617626 , -0.2816758 , -0.51384646]], dtype=float32)>,\n",
       " <tf.Variable 'beta_regular_class:0' shape=(1, 1) dtype=float32, numpy=array([[0.5650171]], dtype=float32)>,\n",
       " <tf.Variable 'beta_seats:0' shape=(1, 1) dtype=float32, numpy=array([[-0.78244746]], dtype=float32)>,\n",
       " <tf.Variable 'beta_luggage=1:0' shape=(1, 1) dtype=float32, numpy=array([[0.42275968]], dtype=float32)>,\n",
       " <tf.Variable 'beta_luggage>1:0' shape=(1, 1) dtype=float32, numpy=array([[1.4139789]], dtype=float32)>,\n",
       " <tf.Variable 'beta_tt_car:0' shape=(1, 1) dtype=float32, numpy=array([[-0.7229831]], dtype=float32)>,\n",
       " <tf.Variable 'beta_travel_time:0' shape=(1, 1) dtype=float32, numpy=array([[-0.69901365]], dtype=float32)>,\n",
       " <tf.Variable 'beta_survey:0' shape=(1, 1) dtype=float32, numpy=array([[2.542475]], dtype=float32)>]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "swiss_model.trainable_weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily acces the negative log likelihood value for the training dataset or another one using the .evaluate() method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=5159.3047>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(swiss_dataset) * swiss_model.evaluate(swiss_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you set get_report to True in .fit, the model automatically creates a report for each of the coefficient, with its estimation, its standard deviation and more:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's add items for better readability\n",
    "items = [\"train\", \"sm\", \"train\", \"sm\", \"train\", \"sm\", \"car\", \"train\", \"sm\", \"car\", \"car\", \"car\", \"train & sm\", \"train & sm\"]\n",
    "swiss_model.report = pd.concat([swiss_model.report, pd.Series(items, name=\"item\")], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "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>Coefficient Name</th>\n",
       "      <th>Coefficient Estimation</th>\n",
       "      <th>Std. Err</th>\n",
       "      <th>z_value</th>\n",
       "      <th>P(.&gt;z)</th>\n",
       "      <th>item</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>beta_intercept_0</td>\n",
       "      <td>-1.292930</td>\n",
       "      <td>0.145000</td>\n",
       "      <td>-8.916755</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>train</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>beta_intercept_1</td>\n",
       "      <td>-0.502575</td>\n",
       "      <td>0.109712</td>\n",
       "      <td>-4.580845</td>\n",
       "      <td>4.649162e-06</td>\n",
       "      <td>sm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>beta_he_0</td>\n",
       "      <td>-0.314336</td>\n",
       "      <td>0.064734</td>\n",
       "      <td>-4.855812</td>\n",
       "      <td>1.192093e-06</td>\n",
       "      <td>train</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>beta_he_1</td>\n",
       "      <td>-0.377317</td>\n",
       "      <td>0.191370</td>\n",
       "      <td>-1.971662</td>\n",
       "      <td>4.864824e-02</td>\n",
       "      <td>sm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>beta_cost_0</td>\n",
       "      <td>-0.561763</td>\n",
       "      <td>0.094500</td>\n",
       "      <td>-5.944569</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>train</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>beta_cost_1</td>\n",
       "      <td>-0.281676</td>\n",
       "      <td>0.042252</td>\n",
       "      <td>-6.666523</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>sm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>beta_cost_2</td>\n",
       "      <td>-0.513846</td>\n",
       "      <td>0.101592</td>\n",
       "      <td>-5.057931</td>\n",
       "      <td>4.768372e-07</td>\n",
       "      <td>car</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>beta_regular_class</td>\n",
       "      <td>0.565017</td>\n",
       "      <td>0.079958</td>\n",
       "      <td>7.066461</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>train</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>beta_seats</td>\n",
       "      <td>-0.782447</td>\n",
       "      <td>0.085864</td>\n",
       "      <td>-9.112632</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>sm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>beta_luggage=1</td>\n",
       "      <td>0.422760</td>\n",
       "      <td>0.063327</td>\n",
       "      <td>6.675845</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>car</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>beta_luggage&gt;1</td>\n",
       "      <td>1.413979</td>\n",
       "      <td>0.190378</td>\n",
       "      <td>7.427222</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>car</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>beta_tt_car</td>\n",
       "      <td>-0.722983</td>\n",
       "      <td>0.044959</td>\n",
       "      <td>-16.080963</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>car</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>beta_travel_time</td>\n",
       "      <td>-0.699014</td>\n",
       "      <td>0.038635</td>\n",
       "      <td>-18.092785</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>train &amp; sm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>beta_survey</td>\n",
       "      <td>2.542475</td>\n",
       "      <td>0.098162</td>\n",
       "      <td>25.900934</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>train &amp; sm</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Coefficient Name  Coefficient Estimation  Std. Err    z_value  \\\n",
       "0     beta_intercept_0               -1.292930  0.145000  -8.916755   \n",
       "1     beta_intercept_1               -0.502575  0.109712  -4.580845   \n",
       "2            beta_he_0               -0.314336  0.064734  -4.855812   \n",
       "3            beta_he_1               -0.377317  0.191370  -1.971662   \n",
       "4          beta_cost_0               -0.561763  0.094500  -5.944569   \n",
       "5          beta_cost_1               -0.281676  0.042252  -6.666523   \n",
       "6          beta_cost_2               -0.513846  0.101592  -5.057931   \n",
       "7   beta_regular_class                0.565017  0.079958   7.066461   \n",
       "8           beta_seats               -0.782447  0.085864  -9.112632   \n",
       "9       beta_luggage=1                0.422760  0.063327   6.675845   \n",
       "10      beta_luggage>1                1.413979  0.190378   7.427222   \n",
       "11         beta_tt_car               -0.722983  0.044959 -16.080963   \n",
       "12    beta_travel_time               -0.699014  0.038635 -18.092785   \n",
       "13         beta_survey                2.542475  0.098162  25.900934   \n",
       "\n",
       "          P(.>z)        item  \n",
       "0   0.000000e+00       train  \n",
       "1   4.649162e-06          sm  \n",
       "2   1.192093e-06       train  \n",
       "3   4.864824e-02          sm  \n",
       "4   0.000000e+00       train  \n",
       "5   0.000000e+00          sm  \n",
       "6   4.768372e-07         car  \n",
       "7   0.000000e+00       train  \n",
       "8   0.000000e+00          sm  \n",
       "9   0.000000e+00         car  \n",
       "10  0.000000e+00         car  \n",
       "11  0.000000e+00         car  \n",
       "12  0.000000e+00  train & sm  \n",
       "13  0.000000e+00  train & sm  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "swiss_model.report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find the same results (estimation of parameters and negative log-likelihood) as the PyLogit package.\n",
    "One can easily interpret the coefficient. For example $\\beta_{cost}$ represent the average price elasticity of the customers. First, it's negative meaning that the most expensive it is the less likely the alternative will be chose, second we can observe that the elasticity is smaller for the SwissMetro. It is not suprising for a \"premium\" product, meaning that people choosing this alternative agreed to pay more for comfort for example and are less regarding on the price.\n",
    "\n",
    "We can also observe that $\\beta_{luggage>1} > \\beta_{luggage=1} >  \\beta_{luggage=0} = 0$ meaning that customers with luggae are more likely to use their car for transport, and it is event further the case if they have more than one piece of luggage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example #2: the ModeCanada Dataset\n",
    "\n",
    "\n",
    "### Utility formulation\n",
    "\n",
    "Let's reproduce a common example from [Torch-Choice](https://gsbdbi.github.io/torch-choice/conditional_logit_model_mode_canada/) on the ModeCanada [1] dataset:\n",
    "$$\n",
    "U(i, c) = \\beta^{inter}_i + \\beta^{price} \\cdot price(i, c) + \\beta^{freq} \\cdot freq(i, c) + \\beta^{ovt} \\cdot ovt(i, c) + \\beta^{income}_i \\cdot income(c) + \\beta^{ivt}_i \\cdot ivt(i, c) + \\epsilon(i, c)\n",
    "$$\n",
    "\n",
    "A description of the dataset and the features can be found in [1]. We want to predict the future mean of transport from [train, air, bus, car] using frequence, price, in-vehicule transport time (ivt) and out-of-vehicule transport time (ovt)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "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>case</th>\n",
       "      <th>alt</th>\n",
       "      <th>choice</th>\n",
       "      <th>dist</th>\n",
       "      <th>cost</th>\n",
       "      <th>ivt</th>\n",
       "      <th>ovt</th>\n",
       "      <th>freq</th>\n",
       "      <th>income</th>\n",
       "      <th>urban</th>\n",
       "      <th>noalt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>train</td>\n",
       "      <td>0</td>\n",
       "      <td>83</td>\n",
       "      <td>28.25</td>\n",
       "      <td>50</td>\n",
       "      <td>66</td>\n",
       "      <td>4</td>\n",
       "      <td>45.0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>car</td>\n",
       "      <td>1</td>\n",
       "      <td>83</td>\n",
       "      <td>15.77</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>train</td>\n",
       "      <td>0</td>\n",
       "      <td>83</td>\n",
       "      <td>28.25</td>\n",
       "      <td>50</td>\n",
       "      <td>66</td>\n",
       "      <td>4</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>car</td>\n",
       "      <td>1</td>\n",
       "      <td>83</td>\n",
       "      <td>15.77</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>train</td>\n",
       "      <td>0</td>\n",
       "      <td>83</td>\n",
       "      <td>28.25</td>\n",
       "      <td>50</td>\n",
       "      <td>66</td>\n",
       "      <td>4</td>\n",
       "      <td>70.0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  case    alt  choice  dist   cost  ivt  ovt  freq  income  \\\n",
       "0           1     1  train       0    83  28.25   50   66     4    45.0   \n",
       "1           2     1    car       1    83  15.77   61    0     0    45.0   \n",
       "2           3     2  train       0    83  28.25   50   66     4    25.0   \n",
       "3           4     2    car       1    83  15.77   61    0     0    25.0   \n",
       "4           5     3  train       0    83  28.25   50   66     4    70.0   \n",
       "\n",
       "   urban  noalt  \n",
       "0      0      2  \n",
       "1      0      2  \n",
       "2      0      2  \n",
       "3      0      2  \n",
       "4      0      2  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# If you want to check what's in the dataset:\n",
    "from choice_learn.datasets import load_modecanada\n",
    "\n",
    "transport_df = load_modecanada(as_frame=True)\n",
    "transport_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to estimate:\n",
    "\n",
    "- one $\\beta^{price}$, $\\beta^{freq}$ and $\\beta^{ovt}$ coefficient. They are **shared** by all items.\n",
    "- one $\\beta^{ivt}$ coefficient for **each** item.\n",
    "- one $\\beta^{inter}$ and $\\beta^{income}$ coefficient for **each** item, with **additional constraint** to be 0 for the first item (air).\n",
    "\n",
    "One notes that it makes sense to include an intercept $\\beta^{inter}$ for each item since $ivt(i, c)$ and $income(c)$ depends on each choice $c$.\n",
    "\n",
    "Additionally to previous example we manually specify the weights names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%=====================================================================%\n",
      "%%% Summary of the dataset:\n",
      "%=====================================================================%\n",
      "Number of items: 4\n",
      "Number of choices: 2779\n",
      "%=====================================================================%\n",
      " Shared Features by Choice:\n",
      " 1 shared features\n",
      " with names: (['income'],)\n",
      "\n",
      "\n",
      " Items Features by Choice:\n",
      "4 items features \n",
      " with names: (['cost', 'freq', 'ovt', 'ivt'],)\n",
      "%=====================================================================%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Loading the ChoiceDataset\n",
    "canada_dataset = load_modecanada(as_frame=False, preprocessing=\"tutorial\")\n",
    "\n",
    "print(canada_dataset.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from choice_learn.models import ConditionalLogit\n",
    "\n",
    "# Initialization of the model\n",
    "model = ConditionalLogit()\n",
    "\n",
    "# Creation of the different weights:\n",
    "\n",
    "# shared_coefficient add one coefficient that is used for all items specified in the items_indexes:\n",
    "# Here, cost, freq and ovt coefficients are shared between all items\n",
    "model.add_shared_coefficient(feature_name=\"cost\", items_indexes=[0, 1, 2, 3])\n",
    "# You can specify you own coefficient name\n",
    "model.add_shared_coefficient(feature_name=\"freq\",\n",
    "                             coefficient_name=\"beta_frequence\",\n",
    "                             items_indexes=[0, 1, 2, 3])\n",
    "model.add_shared_coefficient(feature_name=\"ovt\", items_indexes=[0, 1, 2, 3])\n",
    "\n",
    "# ivt is added for each item:\n",
    "model.add_coefficients(feature_name=\"ivt\", items_indexes=[0, 1, 2, 3])\n",
    "\n",
    "# add_coefficients adds one coefficient for each specified item_index\n",
    "# intercept, and income are added for each item except the first one that needs to be zeroed\n",
    "model.add_coefficients(feature_name=\"intercept\", items_indexes=[1, 2, 3])\n",
    "model.add_coefficients(feature_name=\"income\", items_indexes=[1, 2, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "history = model.fit(canada_dataset, get_report=True, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'beta_cost:0' shape=(1, 1) dtype=float32, numpy=array([[-0.03333882]], dtype=float32)>,\n",
       " <tf.Variable 'beta_frequence:0' shape=(1, 1) dtype=float32, numpy=array([[0.09252929]], dtype=float32)>,\n",
       " <tf.Variable 'beta_ovt:0' shape=(1, 1) dtype=float32, numpy=array([[-0.04300351]], dtype=float32)>,\n",
       " <tf.Variable 'beta_ivt:0' shape=(1, 4) dtype=float32, numpy=\n",
       " array([[ 0.05950943, -0.00678369, -0.00646029, -0.00145037]],\n",
       "       dtype=float32)>,\n",
       " <tf.Variable 'beta_intercept:0' shape=(1, 3) dtype=float32, numpy=array([[0.6983663, 1.8440988, 3.2741835]], dtype=float32)>,\n",
       " <tf.Variable 'beta_income:0' shape=(1, 3) dtype=float32, numpy=array([[-0.08908683, -0.02799302, -0.03814651]], dtype=float32)>]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.trainable_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The average neg-loglikelihood is: 0.67447394\n",
      "The total neg-loglikelihood is: 1874.3630829453468\n"
     ]
    }
   ],
   "source": [
    "print(\"The average neg-loglikelihood is:\", model.evaluate(canada_dataset).numpy())\n",
    "print(\"The total neg-loglikelihood is:\", model.evaluate(canada_dataset).numpy()*len(canada_dataset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "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>Coefficient Name</th>\n",
       "      <th>Coefficient Estimation</th>\n",
       "      <th>Std. Err</th>\n",
       "      <th>z_value</th>\n",
       "      <th>P(.&gt;z)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>beta_cost</td>\n",
       "      <td>-0.033339</td>\n",
       "      <td>0.007095</td>\n",
       "      <td>-4.699046</td>\n",
       "      <td>2.622604e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>beta_frequence</td>\n",
       "      <td>0.092529</td>\n",
       "      <td>0.005097</td>\n",
       "      <td>18.151926</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>beta_ovt</td>\n",
       "      <td>-0.043004</td>\n",
       "      <td>0.003225</td>\n",
       "      <td>-13.335685</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>beta_ivt_0</td>\n",
       "      <td>0.059509</td>\n",
       "      <td>0.010073</td>\n",
       "      <td>5.908032</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>beta_ivt_1</td>\n",
       "      <td>-0.006784</td>\n",
       "      <td>0.004433</td>\n",
       "      <td>-1.530138</td>\n",
       "      <td>1.259825e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>beta_ivt_2</td>\n",
       "      <td>-0.006460</td>\n",
       "      <td>0.001898</td>\n",
       "      <td>-3.403067</td>\n",
       "      <td>6.663799e-04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>beta_ivt_3</td>\n",
       "      <td>-0.001450</td>\n",
       "      <td>0.001187</td>\n",
       "      <td>-1.221395</td>\n",
       "      <td>2.219366e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>beta_intercept_0</td>\n",
       "      <td>0.698366</td>\n",
       "      <td>1.280191</td>\n",
       "      <td>0.545517</td>\n",
       "      <td>5.853977e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>beta_intercept_1</td>\n",
       "      <td>1.844099</td>\n",
       "      <td>0.708427</td>\n",
       "      <td>2.603088</td>\n",
       "      <td>9.238839e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>beta_intercept_2</td>\n",
       "      <td>3.274184</td>\n",
       "      <td>0.624342</td>\n",
       "      <td>5.244213</td>\n",
       "      <td>1.192093e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>beta_income_0</td>\n",
       "      <td>-0.089087</td>\n",
       "      <td>0.018347</td>\n",
       "      <td>-4.855648</td>\n",
       "      <td>1.192093e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>beta_income_1</td>\n",
       "      <td>-0.027993</td>\n",
       "      <td>0.003872</td>\n",
       "      <td>-7.228693</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>beta_income_2</td>\n",
       "      <td>-0.038147</td>\n",
       "      <td>0.004083</td>\n",
       "      <td>-9.342735</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Coefficient Name  Coefficient Estimation  Std. Err    z_value  \\\n",
       "0          beta_cost               -0.033339  0.007095  -4.699046   \n",
       "1     beta_frequence                0.092529  0.005097  18.151926   \n",
       "2           beta_ovt               -0.043004  0.003225 -13.335685   \n",
       "3         beta_ivt_0                0.059509  0.010073   5.908032   \n",
       "4         beta_ivt_1               -0.006784  0.004433  -1.530138   \n",
       "5         beta_ivt_2               -0.006460  0.001898  -3.403067   \n",
       "6         beta_ivt_3               -0.001450  0.001187  -1.221395   \n",
       "7   beta_intercept_0                0.698366  1.280191   0.545517   \n",
       "8   beta_intercept_1                1.844099  0.708427   2.603088   \n",
       "9   beta_intercept_2                3.274184  0.624342   5.244213   \n",
       "10     beta_income_0               -0.089087  0.018347  -4.855648   \n",
       "11     beta_income_1               -0.027993  0.003872  -7.228693   \n",
       "12     beta_income_2               -0.038147  0.004083  -9.342735   \n",
       "\n",
       "          P(.>z)  \n",
       "0   2.622604e-06  \n",
       "1   0.000000e+00  \n",
       "2   0.000000e+00  \n",
       "3   0.000000e+00  \n",
       "4   1.259825e-01  \n",
       "5   6.663799e-04  \n",
       "6   2.219366e-01  \n",
       "7   5.853977e-01  \n",
       "8   9.238839e-03  \n",
       "9   1.192093e-07  \n",
       "10  1.192093e-06  \n",
       "11  0.000000e+00  \n",
       "12  0.000000e+00  "
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Faster Specification\n",
    "\n",
    "A faster specification can be done using a dictionnary. It follows torch-choice method to create conditional logit models.\n",
    "The parameters dict needs to be as follows:\n",
    "- The key is the feature name\n",
    "- The value is the mode. Currently three modes are available:\n",
    "    - constant: the learned coefficient is shared by all items\n",
    "    - item: one coefficient by item is estimated, the value for the item at index 0 is set to 0\n",
    "    - item-full: one coefficient by item is estimated\n",
    "\n",
    "In order to create the same model for the ModeCanada dataset, it looks as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiation with the coefficients dictionnary\n",
    "coefficients = {\"income\": \"item\",\n",
    " \"cost\": \"constant\",\n",
    " \"freq\": \"constant\",\n",
    " \"ovt\": \"constant\",\n",
    " \"ivt\": \"item-full\",\n",
    " \"intercept\": \"item\"}\n",
    "\n",
    "# Instantiation of the model\n",
    "cmnl = ConditionalLogit(coefficients=coefficients, epochs=1000, optimizer=\"lbfgs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[<tf.Variable 'income_w_0:0' shape=(1, 3) dtype=float32, numpy=array([[-0.08908679, -0.02799301, -0.0381465 ]], dtype=float32)>, <tf.Variable 'cost_w_1:0' shape=(1, 1) dtype=float32, numpy=array([[-0.03333884]], dtype=float32)>, <tf.Variable 'freq_w_2:0' shape=(1, 1) dtype=float32, numpy=array([[0.09252931]], dtype=float32)>, <tf.Variable 'ovt_w_3:0' shape=(1, 1) dtype=float32, numpy=array([[-0.04300352]], dtype=float32)>, <tf.Variable 'ivt_w_4:0' shape=(1, 4) dtype=float32, numpy=\n",
      "array([[ 0.05950936, -0.00678372, -0.0064603 , -0.00145037]],\n",
      "      dtype=float32)>, <tf.Variable 'intercept_w_5:0' shape=(1, 3) dtype=float32, numpy=array([[0.69836473, 1.8440932 , 3.2741785 ]], dtype=float32)>]\n",
      "0.67447394\n"
     ]
    }
   ],
   "source": [
    "history = cmnl.fit(canada_dataset)\n",
    "print(cmnl.trainable_weights)\n",
    "print(cmnl.evaluate(canada_dataset).numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison with other implementations results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using L-BFGS optimizer, setting up .fit() function\n",
      "'Ground Truth' Negative LogLikelihood: tf.Tensor(1874.3633, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "\n",
    "# Here are the values obtained in the references:\n",
    "gt_weights = [\n",
    "    tf.constant([[-0.0890796, -0.0279925, -0.038146]]),\n",
    "    tf.constant([[-0.0333421]]),\n",
    "    tf.constant([[0.0925304]]),\n",
    "    tf.constant([[-0.0430032]]),\n",
    "    tf.constant([[0.0595089, -0.00678188, -0.00645982, -0.00145029]]),\n",
    "    tf.constant([[0.697311, 1.8437, 3.27381]]),\n",
    "]\n",
    "gt_model = ConditionalLogit(coefficients=coefficients)\n",
    "gt_model.instantiate(canada_dataset)\n",
    "canada_dataset\n",
    "# Here we estimate the negative log-likelihood with these coefficients (also, we obtain same value as in those papers):\n",
    "gt_model._trainable_weights = gt_weights\n",
    "print(\"'Ground Truth' Negative LogLikelihood:\", gt_model.evaluate(canada_dataset) * len(canada_dataset))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate Utility & probabilities\n",
    "In order to estimate the utilities, use the .predict_utility() method. In order to estimate the probabilities, use the .predict_probas() method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Utilities of each item for the first 5 sessions: tf.Tensor(\n",
      "[[ -4.250798   -8.238913   -3.4962516  -3.5083752]\n",
      " [ -4.250798  -10.466084   -4.196077   -4.4620376]\n",
      " [ -4.250798   -7.348045   -3.2163215  -3.1269102]\n",
      " [ -4.250798  -10.466084   -4.196077   -4.4620376]\n",
      " [ -4.250798  -10.466084   -4.196077   -4.4620376]], shape=(5, 4), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "print(\"Utilities of each item for the first 5 sessions:\", cmnl.compute_batch_utility(*canada_dataset.batch[:5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Purchase probability of each item for the first 5 sessions: tf.Tensor(\n",
      "[[0.1906133  0.00353295 0.40536723 0.40048242]\n",
      " [0.348695   0.00069692 0.3683078  0.2822966 ]\n",
      " [0.14418297 0.00651324 0.40567806 0.44362125]\n",
      " [0.348695   0.00069692 0.3683078  0.2822966 ]\n",
      " [0.348695   0.00069692 0.3683078  0.2822966 ]], shape=(5, 4), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "print(\"Purchase probability of each item for the first 5 sessions:\", cmnl.predict_probas(canada_dataset)[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Gradient Descent Optimizers\n",
    "\n",
    "For very large datasets that do not fit entirely in the memory, we have to work with data batches. In this case, the LBFGS method cannot be used, because of the large memory usage of the algorithm. In those cases, we will prefer stochastic gradient descent optimizers.\n",
    "\n",
    "In this case, it is possible to obtain the same coefficients estimation, also it is a little tricky to get it quickly. We need to adjust the learning rate over time for the optimization not to be too slow.\n",
    "L-BFGS is more efficient for small dataset - Gradient Descent for large ones !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Epoch 1999 Train Loss 0.6801: 100%|██████████| 2000/2000 [00:11<00:00, 176.83it/s]\n",
      "Epoch 3999 Train Loss 0.6776: 100%|██████████| 4000/4000 [00:20<00:00, 198.77it/s]\n",
      "Epoch 19999 Train Loss 0.6767: 100%|██████████| 20000/20000 [01:38<00:00, 202.55it/s]\n"
     ]
    }
   ],
   "source": [
    "cmnl = ConditionalLogit(coefficients=coefficients, optimizer=\"Adam\", epochs=2000, batch_size=-1)\n",
    "history = cmnl.fit(canada_dataset)\n",
    "cmnl.optimizer.lr = cmnl.optimizer.lr / 5\n",
    "cmnl.epochs = 4000\n",
    "history2 = cmnl.fit(canada_dataset)\n",
    "cmnl.optimizer.lr = cmnl.optimizer.lr  / 10\n",
    "cmnl.epochs = 20000\n",
    "history3 = cmnl.fit(canada_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be useful to look at the loss (negative loglikelihood) over time to see how the estimation goes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Second and third part of the gradient descent.')"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAAF0CAYAAACwmBP1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACTPUlEQVR4nOzdeVxU5f4H8M/MADOsI4vsqyiK4gq57yZuebVyy8otu3ktu2q2mFkuGeXteu1Wav0uaqaplUuWppILamLu5o4GyiKI7KvD9vz+gDk6MiAIOAzzeb+acM4855zvOQznzHeeTSaEECAiIiIiIiIig5IbOgAiIiIiIiIiYoJORERERERE1CAwQSciIiIiIiJqAJigExERERERETUATNCJiIiIiIiIGgAm6EREREREREQNABN0IiIiIiIiogaACToRERERERFRA8AEnYiIiIiIiKgBYILeSKxduxYymUzvY86cObhx4wZkMhnWrl1bZ/tcsWJFnW6vLt26dQsLFizA2bNn63zbZ86cQZ8+faBWqyGTybB8+XK95fLz87FgwQIcPHiwwmsLFiyATCZDampqncfXEPn6+mLSpEnS84MHD0Imk+k9N3Wprt6jMpkMCxYsqPV2GpJdu3Y1umMietAff/yBp59+Gt7e3lAqlXBxcUG3bt3wxhtvGDq0OvHgtdXQtJ9Fbty4UWW5qq4/MpkMr732Wp3tqyr18dmoLtTn9fnzzz9H8+bNYWFhAZlMhszMTL3ljh49igULFuh93dfXF0899VS9xNfQ6Pu8MmnSJPj6+tbrfuvqc2xd/J00RB999BG2b99u6DDqDRP0RmbNmjWIiorSebz++utwc3NDVFQUhg0bVmf7augJ+sKFC+slQZ8yZQqSkpKwadMmREVFYdy4cXrL5efnY+HChfWehBqjTp06ISoqCp06darX/TTk96ih7dq1CwsXLjR0GET1ZufOnejevTuys7OxdOlS7N27F5999hl69OiBzZs3Gzo8k1YX159hw4YhKioKbm5udRRVw1Ff1+ezZ8/i9ddfR79+/bB//35ERUXB1tZWb9mjR49i4cKFlSbwpmz+/PnYtm1bve6jPj/HNgaNPUE3M3QAVLeCgoIQEhKi97WuXbs+dP38/HxYWVnVdViPTUlJCYqLi+t1HxcuXMDLL7+MIUOG1Ot+DE0Igbt378LS0rLOt21nZ1et9yMR0aNaunQp/Pz8sGfPHpiZ3fu4M27cOCxdutSAkVFdaNq0KZo2bfrQcsb0uaa+Y7148SIA4OWXX0bnzp3rbT8NQUFBAVQqFWQyWZ1v29/fv863SXQ/1qCbCH3NuLTNrE+fPo1Ro0bB3t5euujExMRg3LhxcHd3l5oFDhgwQPomz9fXFxcvXkRkZKTUlP5hzX20zda++uorBAQEQKlUonXr1ti0aZNOuTt37mD69Olo3bo1bGxs4OzsjP79++Pw4cN6j2np0qX48MMP4efnB6VSiQMHDuCJJ54AAEyePFmK72HNxS5cuIARI0bA3t4eKpUKHTp0wDfffCO9rm0mVFxcjJUrV0rbrex8az84LFy4UCr7YFPE27dv47nnnoNarYaLiwumTJmCrKwsnTJCCKxYsQIdOnSApaUl7O3tMWrUKMTExFR5PFo//fQT2rVrB6VSiWbNmuGzzz6Tfvf30/5+Vq1ahcDAQCiVSun4Fy5ciC5dusDBwQF2dnbo1KkTwsPDIYTQ2UZRURHeeustuLq6wsrKCj179sTx48crxFRZE/eTJ0/ib3/7GxwcHKBSqdCxY0d8//33OmW0v4cDBw7gH//4B5ycnODo6IhnnnkGt27dkso9yns0OzsbL7/8MhwdHWFjY4PBgwcjOjpab9lr165h/PjxcHZ2hlKpRGBgIL788kudMqWlpfjwww/RsmVLWFpaokmTJmjXrh0+++wznXJXrlzBc889BxcXFyiVSnh7e2PChAnQaDRSmeTkZLzyyivw9PSEhYUF/Pz8sHDhQp0vpLR/E59++imWLVsGPz8/2NjYoFu3bjh27JhUbtKkSVKs93eHaWxN4Mi0paWlwcnJSSc515LLK3782bx5M7p16wZra2vY2Nhg0KBBOHPmTIVyf/zxB4YPHw5HR0eoVCr4+/tj5syZOmWOHDmCAQMGwNbWFlZWVujevTt27typU6a61zKg+tfWylT3Gq5turx792506tQJlpaWaNWqFVavXl1hm8eOHUOPHj2gUqng7u6OuXPnoqio6KGxVPf68+233yIwMBBWVlZo3749fvnlF53X9TXd7du3L4KCgnDo0CF0794dVlZWmDJlCoCyGskxY8bA1tYWarUaY8eORXJycnVOn7SviIgITJ48GQ4ODrC2tsbw4cMr3IsjIiIwYsQIeHp6QqVSoXnz5njllVcqdGmr7DPYo16fV69ejfbt20OlUsHBwQFPP/00Ll++rHNuXnjhBQBAly5d9H4muT+2N998EwDg5+cnxfDgPbs675Pq3Lsqo9Fo8MYbb0jv+969e+PUqVMVunZofz979+7FlClT0LRpU1hZWUGj0eD69euYPHkyWrRoASsrK3h4eGD48OE4f/58hf1duXIFgwcPhpWVFZycnDBt2jTk5ORUKKeviXt1P6tp36MnTpxAr169YGVlhWbNmuHjjz9GaWkpgLLPSI/yObYmf5PVud49LBfQ+u6779CtWzfY2NjAxsYGHTp0QHh4uE6Z3377DQMGDICdnR2srKzQo0cP7Nu3T6eM9m/i4sWLVX42lslkyMvLwzfffCOdm759+1Z5boyOoEZhzZo1AoA4duyYKCoq0nkIIURsbKwAINasWSOt88EHHwgAwsfHR7z99tsiIiJCbN++XQghRMuWLUXz5s3Ft99+KyIjI8WWLVvEG2+8IQ4cOCCEEOL06dOiWbNmomPHjiIqKkpERUWJ06dPVxkjAOHl5SVat24tNm7cKHbs2CEGDx4sAIgffvhBKnflyhXxj3/8Q2zatEkcPHhQ/PLLL+Kll14Scrlc2v/9x+Th4SH69esnfvzxR7F3715x7tw56Xy89957Unzx8fGVxnblyhVha2sr/P39xbp168TOnTvFc889JwCITz75RAghREpKioiKihIAxKhRo6Tt6nP37l2xe/duAUC89NJLUtnr16/rnPuWLVuK999/X0RERIhly5YJpVIpJk+erLOtl19+WZibm4s33nhD7N69W3z33XeiVatWwsXFRSQnJ1d5zn/99Vchl8tF3759xbZt28QPP/wgunTpInx9fcWDf/7ac9muXTvx3Xffif3794sLFy4IIYSYNGmSCA8PFxERESIiIkIsXrxYWFpaioULF+psY+LEiUImk4k333xT7N27Vyxbtkx4eHgIOzs7MXHiRKncgQMHBACd3+f+/fuFhYWF6NWrl9i8ebPYvXu3mDRpUoX3rfZ326xZMzFjxgyxZ88e8b///U/Y29uLfv36SeVq+h4tLS0V/fr1E0qlUixZskTs3btXfPDBB6JZs2YCgPjggw+kshcvXhRqtVq0bdtWrFu3Tuzdu1e88cYbQi6XiwULFkjlwsLChEKhEB988IHYt2+f2L17t1i+fLlOmbNnzwobGxvh6+srVq1aJfbt2yfWr18vxowZI7Kzs4UQQiQlJQkvLy/h4+MjvvrqK/Hbb7+JxYsXC6VSKSZNmiRtS/s34evrKwYPHiy2b98utm/fLtq2bSvs7e1FZmamEEKI69evi1GjRgkA0rmJiooSd+/erfT8EBmbqVOnCgBixowZ4tixY6KwsLDSskuWLBEymUxMmTJF/PLLL2Lr1q2iW7duwtraWly8eFEqt3v3bmFubi7atWsn1q5dK/bv3y9Wr14txo0bJ5U5ePCgMDc3F8HBwWLz5s1i+/btIjQ0VMhkMrFp0yapXHWvZUJU/9pamepew318fISnp6do3bq1WLdundizZ48YPXq0ACAiIyOlchcvXhRWVlbS/fynn34SgwYNEt7e3gKAiI2NrTSWh11/tNewzp07i++//17s2rVL9O3bV5iZmYm//vqrwvm7f199+vQRDg4OwsvLS3z++efiwIEDIjIyUuTn54vAwEChVqvF559/Lvbs2SNef/11Kd777zH6aPfl5eUlpkyZIn799Vfx9ddfC2dnZ+Hl5SUyMjKksitXrhRhYWFix44dIjIyUnzzzTeiffv2omXLljrvwco+gz3K9fmjjz4SAMRzzz0ndu7cKdatWyeaNWsm1Gq1iI6Oln5n7733nnS8938meVB8fLyYMWOGACC2bt0qxZCVlSWEqP77pLr3rso899xzQi6Xi3feeUfs3btXLF++XHh5eQm1Wq3zvtf+fjw8PMTf//538euvv4off/xRFBcXi8jISPHGG2+IH3/8UURGRopt27aJkSNHCktLS3HlyhVpG8nJycLZ2Vl4eHiINWvWiF27donnn39eeo/c/3ll4sSJwsfHRyfW6n5W69Onj3B0dBQtWrQQq1atEhEREWL69OkCgPjmm2+EEEJkZWXV+HNsTf4mq3u9e1guIIQQ8+fPFwDEM888I3744Qfp+jR//nypzLfffitkMpkYOXKk2Lp1q/j555/FU089JRQKhfjtt9+kctX9bBwVFSUsLS3F0KFDpXNzf9yNARP0RkL7h6zvUVRUVGWC/v777+tsKzU1VQAQy5cvr3Kfbdq0EX369Kl2jACEpaWlzoWquLhYtGrVSjRv3rzS9YqLi0VRUZEYMGCAePrpp6Xl2mPy9/ev8MHrxIkT1brpao0bN04olUoRFxens3zIkCHCyspKSmy0x/Hqq68+dJt37typkNhpac/90qVLdZZPnz5dqFQqUVpaKoQQ0hcC//73v3XKxcfHC0tLS/HWW29VGcMTTzwhvLy8hEajkZbl5OQIR0dHvQm6Wq0W6enpVW6zpKREFBUViUWLFglHR0cp1suXLwsAYtasWTrlN2zYIAA8NEFv1aqV6Nixo/SlktZTTz0l3NzcRElJiRDi3nt9+vTpOuWWLl0qAIikpCRpWU3eo7/++qsAID777DOd5UuWLKnwexw0aJDw9PSUPqxovfbaa0KlUknn8KmnnhIdOnSocr/9+/cXTZo0ESkpKZWWeeWVV4SNjY24efOmzvJPP/1UAJBuTNq/ibZt24ri4mKp3PHjxwUAsXHjRmnZq6++WuE9QNSYpKamip49e0r3QnNzc9G9e3cRFhYmcnJypHJxcXHCzMxMzJgxQ2f9nJwc4erqKsaMGSMt8/f3F/7+/qKgoKDS/Xbt2lU4Ozvr7KO4uFgEBQUJT09P6ZpZ3WtZTa6t1VHZNVyIssRLpVLpXGsKCgqEg4ODeOWVV6RlY8eOrfR+/rAEXYiqrz8AhIuLi/QFpRBlyZNcLhdhYWHSssoSdABi3759OttcuXKlACB++uknneUvv/xyjRL0+z+DCCHE77//LgCIDz/8UO96paWloqioSNy8ebPC/iv7DCZEza7PGRkZUrJyv7i4OKFUKsX48eMrHMeJEyceut1//etflf4uq/s+qe69S5+LFy8KAOLtt9/WWb5x48YK73vtcU2YMOGhx1VcXCwKCwtFixYtdP6m3n77bSGTycTZs2d1yg8cOPChCXpNPqtp36N//PGHTtnWrVuLQYMGSc9r+jm2un+T1b3eVScXiImJEQqFQjz//POVlsnLyxMODg5i+PDhOstLSkpE+/btRefOnaVl1f1sLIQQ1tbWNb72GRM2cW9k1q1bhxMnTug89DXvu9+zzz6r89zBwQH+/v7417/+hWXLluHMmTNSs5vaGjBgAFxcXKTnCoUCY8eOxfXr15GQkCAtX7VqFTp16gSVSgUzMzOYm5tj3759Os21tP72t7/B3Ny8VnHt378fAwYMgJeXl87ySZMmIT8/H1FRUbXafmX+9re/6Txv164d7t69i5SUFADAL7/8AplMhhdeeAHFxcXSw9XVFe3bt69yALq8vDycPHkSI0eOhIWFhbTcxsYGw4cP17tO//79YW9vX2H5/v378eSTT0KtVkOhUMDc3Bzvv/8+0tLSpFgPHDgAAHj++ed11h0zZsxD34PXr1/HlStXpHXvP9ahQ4ciKSkJV69e1VlH37kDgJs3b1a5r8pUFv/48eN1nt+9exf79u3D008/DSsrqwqx3r17V2pO3rlzZ5w7dw7Tp0/Hnj17kJ2drbOt/Px8REZGYsyYMVX2pfzll1/Qr18/uLu76+xPOw5CZGSkTvlhw4ZBoVBIz2t7boiMkaOjIw4fPowTJ07g448/xogRIxAdHY25c+eibdu2UpPjPXv2oLi4GBMmTND5+1KpVOjTp490nY2OjsZff/2Fl156CSqVSu8+8/Ly8Mcff2DUqFGwsbGRlisUCrz44otISEio8bWsNtdWrepcw7U6dOgAb29v6blKpUJAQIDO9ePAgQOV3s/rQr9+/XQGL3NxcYGzs3O1rmH29vbo37+/zrIDBw7A1ta2wrl+8Pr+MA/+Drp37w4fHx/pdwQAKSkpmDZtGry8vKTPLz4+PgCg9zPMg5/BaioqKgoFBQUVmqt7eXmhf//+FZoR15XqvE9qeu+6n/a1MWPG6CwfNWpUpe97feeyuLgYH330EVq3bg0LCwuYmZnBwsIC165d0/l9HDhwAG3atEH79u111q/Oe6Smn9VcXV0rjAHQrl27Wt2jq/s3Wd3rXXVygYiICJSUlODVV1+tNK6jR48iPT0dEydO1NlfaWkpBg8ejBMnTiAvL09nnYd9NjYFHCSukQkMDKx0kLjKPDgCqkwmw759+7Bo0SIsXboUb7zxBhwcHPD8889jyZIllY74WR2urq6VLktLS4OnpyeWLVuGN954A9OmTcPixYvh5OQEhUKB+fPn67251cUIrmlpaXq34+7uLr1eHxwdHXWeK5VKAGWDmwBlfdSFEDoX3Ps1a9as0m1nZGRUum5l29N3Do4fP47Q0FD07dsX//d//yf1I9u+fTuWLFkixao9Rw/+js3MzCoc54Nu374NAJgzZw7mzJmjt8yD/fcedu5qKi0tTW+sDx5PWloaiouL8fnnn+Pzzz+vMta5c+fC2toa69evx6pVq6BQKNC7d2988sknCAkJQUZGBkpKSuDp6VllbLdv38bPP/9c6RdR9X1uiIxZSEiIdF8sKirC22+/jf/85z9YunQpli5dKl1/tH0+H6Ttr37nzh0AqPLvVXvdrcn95GF/r7W5tgLVv4ZXFo82pvvLpaWlVXk/r63qxFAZfec+LS1N732vpvFWdsza31FpaSlCQ0Nx69YtzJ8/H23btoW1tTVKS0vRtWtXvfHX9jOMdt+VveciIiJqtf3KVOd3VNN71/20x/Xg762q972+czB79mx8+eWXePvtt9GnTx/Y29tDLpdj6tSpFd7Tfn5+Fdavznukpp/VavP+rkx1/yare72rTi5QnWuidn+jRo2qtEx6ejqsra2l5/wMwwSdAL0Dnfn4+EgDPERHR+P777/HggULUFhYiFWrVj3yvvQNyKJdpv2DXL9+Pfr27YuVK1fqlNM3UEdl8deUo6MjkpKSKizXDtTj5ORU6308CicnJ8hkMhw+fFi6QN1P3zIte3t7yGQy6eJ4v8oGxtF3Ljdt2gRzc3P88ssvOrVGD05vof39JScnw8PDQ1peXFz80C84tOd37ty5eOaZZ/SWadmyZZXbqC1HR0cp1vtvDg+eK3t7e6k2rLJvjbU3eTMzM8yePRuzZ89GZmYmfvvtN7z77rsYNGgQ4uPj4eDgAIVCodN6RB8nJye0a9cOS5Ys0fu69oM/EVXN3NwcH3zwAf7zn//gwoULAO5df3788UepplMfbSuXqv5etR/+6/J+UptrK1D9a3hNY6rqfm5I+u5jjo6OegfVq2m8lR1z8+bNAZQNNnvu3DmsXbsWEydOlMpcv369RvHWhPb9Udl7zlCfX4Da3bu0x3X79u1qv+/1ncv169djwoQJ+Oijj3SWp6amokmTJjr7e9T3dG0+q9WV6sZf3esd8PBc4P5r4oMtUB/c3+eff17p7D2VfbFhypig00MFBATgvffew5YtW3D69Glp+aN827dv3z7cvn1b+mMsKSnB5s2b4e/vL30DJ5PJKlzM/vzzT0RFRVV6AXhQTb9tGzBgALZt24Zbt27p3DDWrVsHKyurR5oSrC6+8Xvqqafw8ccfIzExsUIzr4extrZGSEgItm/fjk8//VRq5p6bm1thNNyqyGQymJmZ6TSZLigowLfffqtTTjuC5oYNGxAcHCwt//777x86WmvLli3RokULnDt3rsJNtDZq8h7t168fli5dig0bNuD111+Xln/33Xc65aysrNCvXz+cOXMG7dq10+k+UJUmTZpg1KhRSExMxMyZM3Hjxg20bt0affr0wQ8//IAlS5ZU+kHqqaeewq5du+Dv76+3C8KjuP/9WR9T6REZWlJSkt4aNW1LLO21ftCgQTAzM8Nff/1VZXPjgIAA+Pv7Y/Xq1Zg9e7beD93W1tbo0qULtm7dik8//VT62yotLcX69evh6emJgICAGh1Hba6tQPWv4TXRr18/7NixQ+/9vDoe9/WnX79++P7777Fjxw6d5rMPXt8fZsOGDTrvkaNHj+LmzZuYOnUqgHsJ4oPvja+++qpG+6nJ+enWrRssLS2xfv16jB49WlqekJCA/fv3V1lzWd0YHlVt7l29e/cGUDbaeKdOnaTlP/74Y42m09X3mXLnzp1ITEyUvlgB7n0GOHfunE4z9+q8R2rzWa0yNT3/1f2brO717kH6coHQ0FAoFAqsXLkS3bp107tejx490KRJE1y6dAmvvfZatff3MLVtcdDQMUGnCv7880+89tprGD16NFq0aAELCwvs378ff/75J9555x2pXNu2bbFp0yZs3rwZzZo1g0qlQtu2bavctpOTE/r374/58+fD2toaK1aswJUrV3SmWnvqqaewePFifPDBB+jTpw+uXr2KRYsWwc/Pr9oXZX9/f1haWmLDhg0IDAyEjY0N3N3dK/229oMPPpD6Sr3//vtwcHDAhg0bsHPnTixduhRqtbpa+72fra0tfHx88NNPP2HAgAFwcHCAk5PTQ6f6ul+PHj3w97//HZMnT8bJkyfRu3dvWFtbIykpCUeOHEHbtm3xj3/8o9L1Fy1ahGHDhmHQoEH45z//iZKSEvzrX/+CjY0N0tPTqxXDsGHDsGzZMowfPx5///vfkZaWhk8//bTCDS8wMBAvvPACli9fDnNzczz55JO4cOECPv30U9jZ2T10P1999RWGDBmCQYMGYdKkSfDw8EB6ejouX76M06dP44cffqhWvPeryXs0NDQUvXv3xltvvYW8vDyEhITg999/1/sh9rPPPkPPnj3Rq1cv/OMf/4Cvry9ycnJw/fp1/Pzzz9i/fz8AYPjw4QgKCkJISAiaNm2KmzdvYvny5fDx8UGLFi0AAMuWLUPPnj3RpUsXvPPOO2jevDlu376NHTt24KuvvoKtrS0WLVqEiIgIdO/eHa+//jpatmyJu3fv4saNG9i1axdWrVr10Gby+s4NAHzyyScYMmQIFAqF9IXDgAEDEBkZWaMPQUQNzaBBg+Dp6Ynhw4ejVatWKC0txdmzZ/Hvf/8bNjY2+Oc//wmgbGqxRYsWYd68eYiJicHgwYNhb2+P27dv4/jx47C2tsbChQsBAF9++SWGDx+Orl27YtasWfD29kZcXBz27NmDDRs2AADCwsIwcOBA9OvXD3PmzIGFhQVWrFiBCxcuYOPGjTWuMa3ttbW61/CaeO+997Bjxw70798f77//PqysrPDll19W6EtamaquP/VhwoQJ+M9//oMJEyZgyZIlaNGiBXbt2oU9e/bUaDsnT57E1KlTMXr0aMTHx2PevHnw8PDA9OnTAQCtWrWCv78/3nnnHQgh4ODggJ9//rnGzcxrcn6aNGmC+fPn491338WECRPw3HPPIS0tDQsXLoRKpcIHH3xQo30/GMNnn32GiRMnwtzcHC1btqxRN8fa3LvatGmD5557Dv/+97+hUCjQv39/XLx4Ef/+97+hVqv1TpWoz1NPPYW1a9eiVatWaNeuHU6dOoV//etfFfY7c+ZMrF69GsOGDcOHH34IFxcXbNiwAVeuXHnoPmr7WU2fmn6Ore7fZHWvd9XJBXx9ffHuu+9i8eLFKCgokKZGu3TpElJTU7Fw4ULY2Njg888/x8SJE5Geno5Ro0bB2dkZd+7cwblz53Dnzp0KLWaro23btjh48CB+/vlnuLm5wdbWFi1btsTNmzfh7++PiRMnVpjqzagYdIg6qjMPG5mzqlHc79y5o1P29u3bYtKkSaJVq1bC2tpa2NjYiHbt2on//Oc/OiND37hxQ4SGhgpbW1tpqpCqoHz08xUrVgh/f39hbm4uWrVqJTZs2KBTTqPRiDlz5ggPDw+hUqlEp06dxPbt2yuMmqk9pn/9619697dx40bRqlUrYW5uXulo6vc7f/68GD58uFCr1cLCwkK0b99e7+iZ2uOojt9++0107NhRKJVKnVFHKzv3+kalFUKI1atXiy5dughra2thaWkp/P39xYQJE8TJkycfGsO2bdtE27ZthYWFhfD29hYff/yxeP3114W9vX21j2v16tWiZcuWQqlUimbNmomwsDARHh5eIVaNRiPeeOMN4ezsLFQqlejatauIiooSPj4+Dx3FXQghzp07J8aMGSOcnZ2Fubm5cHV1Ff379xerVq2qcI4efK/r22ZN36OZmZliypQpokmTJsLKykoMHDhQXLlyRe/7JzY2VkyZMkV4eHgIc3Nz0bRpU9G9e3ed0Xz//e9/i+7duwsnJyfp/L/00kvixo0bOtu6dOmSGD16tHB0dJTKTZo0SWdanTt37ojXX39d+Pn5CXNzc+Hg4CCCg4PFvHnzRG5urhRTZX8TDx6DRqMRU6dOFU2bNhUymUznd6kdZZbImG3evFmMHz9etGjRQtjY2Ahzc3Ph7e0tXnzxRXHp0qUK5bdv3y769esn7OzshFKpFD4+PmLUqFE60wAJUTZi85AhQ4RarRZKpVL4+/tXGGH98OHDon///tI1u2vXruLnn3/WKVOTa1l1r62Vqe413MfHRwwbNqzC+n369KkwI8bvv/8uunbtKpRKpXB1dRVvvvmm+Prrr6s1intV15/K7kUPHmtlo7i3adNG7z4TEhLEs88+K2xsbIStra149tlnxdGjR2s0ivvevXvFiy++KJo0aSKNnH7t2jWdspcuXRIDBw4Utra2wt7eXowePVrExcVVuAZX9jngYeenMv/73/9Eu3bthIWFhVCr1WLEiBEVRkmvySjuQggxd+5c4e7uLuRyuc57sibvk+rcuypz9+5dMXv27Arve7VarfM3V9VxZWRkiJdeekk4OzsLKysr0bNnT3H48GG9sWp/dyqVSjg4OIiXXnpJ/PTTT9WaZk2I6n1Wq+w9qm+bNf0cW5O/yYdd76qbCwghxLp168QTTzwhVCqVsLGxER07dqzwNxUZGSmGDRsmHBwchLm5ufDw8BDDhg3TmWa5Jp+Nz549K3r06CGsrKwEAOl3qf0cZOwjvMuEEKJ+Un8iXTKZDK+++iq++OILQ4di0oqKitChQwd4eHhg7969hg6HiIioQVu7di0mT56MEydO1HggXqpbR48eRY8ePbBhw4Yaj8JPZCzYxJ2okXvppZcwcOBAuLm5ITk5GatWrcLly5fx2WefGTo0IiIiIr0iIiIQFRWF4OBgWFpa4ty5c/j444/RokWLSgeUJWoMmKATNXI5OTmYM2cO7ty5A3Nzc3Tq1Am7du3Ck08+aejQiIiIiPSys7PD3r17sXz5cuTk5MDJyQlDhgxBWFiYzowERI0Nm7gTERERERERNQDVGwKRiIiIiIiIiOoVE3QiIiIiIiKiBqBWCXpYWBhkMhlmzpxZZbnIyEgEBwdDpVKhWbNmWLVqVW12S0RERERERNToPPIgcSdOnMDXX3+Ndu3aVVkuNjYWQ4cOxcsvv4z169fj999/x/Tp09G0aVM8++yz1dpXaWkpbt26BVtbW8hkskcNmYiIqM4IIZCTkwN3d3fI5WyQVlu81xMRUUNjiHv9IyXoubm5eP755/F///d/+PDDD6ssu2rVKnh7e2P58uUAgMDAQJw8eRKffvpptRP0W7duwcvL61FCJSIiqlfx8fHw9PQ0dBhGj/d6IiJqqB7nvf6REvRXX30Vw4YNw5NPPvnQBD0qKgqhoaE6ywYNGoTw8HAUFRXB3Ny8wjoajQYajUZ6rh1oPj4+HnZ2do8SMhERUZ3Kzs6Gl5cXbG1tDR1Ko6A9j7zXExFRQ2GIe32NE/RNmzbh9OnTOHHiRLXKJycnw8XFRWeZi4sLiouLkZqaCjc3twrrhIWFYeHChRWW29nZ8aZNREQNCptj1w3teeS9noiIGprHea+vUUP6+Ph4/POf/8T69euhUqmqvd6DB6StEa/sQOfOnYusrCzpER8fX5MwiYiIiIiIiIxOjWrQT506hZSUFAQHB0vLSkpKcOjQIXzxxRfQaDRQKBQ667i6uiI5OVlnWUpKCszMzODo6Kh3P0qlEkqlsiahERERERERERm1GiXoAwYMwPnz53WWTZ48Ga1atcLbb79dITkHgG7duuHnn3/WWbZ3716EhITo7X9OREREREREZIpqlKDb2toiKChIZ5m1tTUcHR2l5XPnzkViYiLWrVsHAJg2bRq++OILzJ49Gy+//DKioqIQHh6OjRs31tEhEBERERERERm/Op/MLSkpCXFxcdJzPz8/7Nq1CwcPHkSHDh2wePFi/Pe//632FGtEREREREREpkAmtCO2NWDZ2dlQq9XIysriyK5ERNQg8N5Ut3g+iYiooTHEvanOa9CJiIiIiIiIqOaYoBMRERERERE1AEzQiYiIiIiIiBoAJuhERETU6JSWNvghdoiIiCowuQR97FdReHblUaTnFRo6FCIiIqoHy/ZeRc9P9uPa7RxDh0JERFQjJpegn7qZgVM3M1BYXGroUIiIiKgeXErKxq2su9h+NtHQoRAREdWIySXocpkMACDApm9ERESN0YgOHgCAn87eYlN3IiIyKiaXoKMsPwfv10RERI3Tk4EusFGaISGjAKfiMgwdDhERUbWZXIIuL0/QhWCGTkRE1BhZWigwqI0rAGD7GTZzJyIi42FyCbqsvAqd+TkREVHjNbKjOwBg5/kkjjtDRERGw/QSdKkG3bBxEBERPYoVK1bAz88PKpUKwcHBOHz4cJXlNRoN5s2bBx8fHyiVSvj7+2P16tU6ZZYvX46WLVvC0tISXl5emDVrFu7evVuj/QohsGDBAri7u8PS0hJ9+/bFxYsX6+agH0F3fyc0tVUiM78IkdF3DBYHERFRTZhcgq4dJK6UGToRERmZzZs3Y+bMmZg3bx7OnDmDXr16YciQIYiLi6t0nTFjxmDfvn0IDw/H1atXsXHjRrRq1Up6fcOGDXjnnXfwwQcf4PLlywgPD8fmzZsxd+7cGu136dKlWLZsGb744gucOHECrq6uGDhwIHJyDDPVmUIuw9/al9WiczR3IiIyFjJhBJ2xs7OzoVarkZWVBTs7u1ptq+0He5CjKcaBOX3h52RdRxESEZGpqct7U3V16dIFnTp1wsqVK6VlgYGBGDlyJMLCwiqU3717N8aNG4eYmBg4ODjo3eZrr72Gy5cvY9++fdKyN954A8ePH5dqyR+2XyEE3N3dMXPmTLz99tsAymruXVxc8Mknn+CVV1556LHVx/m8kJiFpz4/AqWZHCffexK2KvM62S4REZkGQ9zrTa4GXcZB4oiIyAgVFhbi1KlTCA0N1VkeGhqKo0eP6l1nx44dCAkJwdKlS+Hh4YGAgADMmTMHBQUFUpmePXvi1KlTOH78OAAgJiYGu3btwrBhw6q939jYWCQnJ+uUUSqV6NOnT6WxaTQaZGdn6zzqWht3O/g3tYamuBS7LyTX+faJiIjqmpmhA3jcZFITdwMHQkREVAOpqakoKSmBi4uLznIXFxckJ+tPPmNiYnDkyBGoVCps27YNqampmD59OtLT06V+6OPGjcOdO3fQs2dPCCFQXFyMf/zjH3jnnXeqvV/tT31lbt68qTe2sLAwLFy4sIZnoWZkMhlGdvDAvyOisf1sIkaHeNXr/oiIiGrL5GrQtdOsAczQiYjI+Gi/aNYSQlRYplVaWgqZTIYNGzagc+fOGDp0KJYtW4a1a9dKtegHDx7EkiVLsGLFCpw+fRpbt27FL7/8gsWLF9d4vzWJbe7cucjKypIe8fHxDz/4RzCigwcA4OhfabidffchpYmIiAzL5BJ01qATEZExcnJygkKhqFBbnpKSUqHmWsvNzQ0eHh5Qq9XSssDAQAghkJCQAACYP38+XnzxRUydOhVt27bF008/jY8++ghhYWEoLS2t1n5dXcvmHK9JbEqlEnZ2djqP+uDtaIWO3k0gBLDrfFK97IOIiKiumFyCLuc0a0REZIQsLCwQHByMiIgIneURERHo3r273nV69OiBW7duITc3V1oWHR0NuVwOT09PAEB+fj7kct2PAwqFAkIICCGqtV8/Pz+4urrqlCksLERkZGSlsT1OT7UrG839lz+ZoBMRUcNmcgk6wGnWiIjIOM2ePRv/+9//sHr1aly+fBmzZs1CXFwcpk2bBqCs2fiECROk8uPHj4ejoyMmT56MS5cu4dChQ3jzzTcxZcoUWFpaAgCGDx+OlStXYtOmTYiNjUVERATmz5+Pv/3tb1AoFNXar0wmw8yZM/HRRx9h27ZtuHDhAiZNmgQrKyuMHz/+MZ+lioa1dYNMBpy6mYFbmQUPX4GIiMhATG6QONagExGRsRo7dizS0tKwaNEiJCUlISgoCLt27YKPjw8AICkpSWduchsbG0RERGDGjBkICQmBo6MjxowZgw8//FAq895770Emk+G9995DYmIimjZtiuHDh2PJkiXV3i8AvPXWWygoKMD06dORkZGBLl26YO/evbC1tX0MZ6ZqrmoVnvB1wPHYdOw6n4SpvZoZOiQiIiK9TG4e9C4f/Ybb2Rr8MqMngjzUD1+BiIhID0PMjdqY1ff5XBd1A+//dBEdvZtg2/Qedb59IiJqfDgP+mMgr2Q0WSIiImq8Bge5QiYDzsRlIiEj39DhEBER6WVyCbo2PWcfdCIiItPhbKtCZ18HAMDuC/rnjSciIjI000vQy2vQmZ8TERGZlsFBZdPB7bnIBJ2IiBomE0zQy36yBp2IiMi0DGpTlqCfvJmBlJy7Bo6GiIioIpNL0LV90JmeExERmRb3JpZo79UEQgARl24bOhwiIqIKTC5Bl0nTrDFFJyIiMjWDy2vR2Q+diIgaItNL0Mt/Mj8nIiIyPdp+6FF/pSEzv9DA0RAREekyuQRd28S9lAk6ERGRyfFzskYrV1sUlwrsu5xi6HCIiIh0mFyCDjZxJyIiMmnaweJ+ZTN3IiJqYEwuQWcNOhERkWnTNnM/dO0O8jTFBo6GiIjoHpNL0KU+6BzHnYiIyCS1crWFj6MVCotLcfDqHUOHQ0REJDG5BF2aZo35ORERkUmSyWRSLfrui2zmTkREDUeNEvSVK1eiXbt2sLOzg52dHbp164Zff/210vIHDx6ETCar8Lhy5UqtA39U96ZZM1gIREREZGDa6db2X76Nu0UlBo6GiIiojFlNCnt6euLjjz9G8+bNAQDffPMNRowYgTNnzqBNmzaVrnf16lXY2dlJz5s2bfqI4daeTOqDzgydiIjIVLX3bAIXOyVuZ2twLCYNfVs6GzokIiKimtWgDx8+HEOHDkVAQAACAgKwZMkS2NjY4NixY1Wu5+zsDFdXV+mhUChqFXRt3OuDTkRERKZKLpehfysXAMD+K5xujYiIGoZH7oNeUlKCTZs2IS8vD926dauybMeOHeHm5oYBAwbgwIEDD922RqNBdna2zqOuyMuPmNOsERERmbYBrcpqzfddTuHnAiIiahBqnKCfP38eNjY2UCqVmDZtGrZt24bWrVvrLevm5oavv/4aW7ZswdatW9GyZUsMGDAAhw4dqnIfYWFhUKvV0sPLy6umYVZKBg4SR0RERECP5k5QmsmRmFmA6Nu5hg6HiIioZn3QAaBly5Y4e/YsMjMzsWXLFkycOBGRkZF6k/SWLVuiZcuW0vNu3bohPj4en376KXr37l3pPubOnYvZs2dLz7Ozs+ssSZdrB4ljI3ciIiKTZmmhQHd/Rxy4egf7rtxGS1dbQ4dEREQmrsY16BYWFmjevDlCQkIQFhaG9u3b47PPPqv2+l27dsW1a9eqLKNUKqWR4rWPOqMdJK607jZJRERExql/YHk/9Mvsh05ERIZX63nQhRDQaDTVLn/mzBm4ubnVdreP7F4NOhEREZm6/uX90E/HZSA9r9DA0RARkamrURP3d999F0OGDIGXlxdycnKwadMmHDx4ELt37wZQ1jQ9MTER69atAwAsX74cvr6+aNOmDQoLC7F+/Xps2bIFW7ZsqfsjqSbtKO6cZo2IiIg8mliilastriTnIDI6BU939DR0SEREZMJqlKDfvn0bL774IpKSkqBWq9GuXTvs3r0bAwcOBAAkJSUhLi5OKl9YWIg5c+YgMTERlpaWaNOmDXbu3ImhQ4fW7VHUgFzGQeKIiIjongGBzriSnIN9l5mgExGRYdUoQQ8PD6/y9bVr1+o8f+utt/DWW2/VOKj6JNM2cWeGTkRERAAGBLrgywN/ITL6DopKSmGuqHUPQCIiokdicncgmbYG3cBxEBERUcPQ3rMJHK0tkHO3GCdvZBg6HCIiMmGml6CX/2QfdCIiIgIAhVyGvi3LBovbf+W2gaMhIiJTZnoJutTE3bBxEBERUcMxILAsQd93hdOtERGR4Zhcgq4dJI416ERERKTVq4UTzOQyxNzJQ2xqnqHDISIiE2VyCbq2Bp2IiIhIy1Zlji7NHAAA+1mLTkREBmJyCTpr0ImIiEif/q1cALAfOhERGY7JJehazM+JiMgYrVixAn5+flCpVAgODsbhw4erLK/RaDBv3jz4+PhAqVTC398fq1evll7v27cvZDJZhcewYcOkMjk5OZg5cyZ8fHxgaWmJ7t2748SJEzr7mTRpUoVtdO3atW4Pvp4NaFXWD/2PmHTk3C0ycDRERGSKajQPemNwrwbdwIEQERHV0ObNmzFz5kysWLECPXr0wFdffYUhQ4bg0qVL8Pb21rvOmDFjcPv2bYSHh6N58+ZISUlBcXGx9PrWrVtRWFgoPU9LS0P79u0xevRoadnUqVNx4cIFfPvtt3B3d8f69evx5JNP4tKlS/Dw8JDKDR48GGvWrJGeW1hY1OXh1ztfJ2s0a2qNmDt5OHwtFUPbuhk6JCIiMjEml6DfG8WdGToRERmXZcuW4aWXXsLUqVMBAMuXL8eePXuwcuVKhIWFVSi/e/duREZGIiYmBg4OZf2rfX19dcpol2tt2rQJVlZWUoJeUFCALVu24KeffkLv3r0BAAsWLMD27duxcuVKfPjhh9K6SqUSrq6udXa8hjCglTNi7sRi3+UUJuhERPTYmVwTd20NOvNzIiIyJoWFhTh16hRCQ0N1loeGhuLo0aN619mxYwdCQkKwdOlSeHh4ICAgAHPmzEFBQUGl+wkPD8e4ceNgbW0NACguLkZJSQlUKpVOOUtLSxw5ckRn2cGDB+Hs7IyAgAC8/PLLSEkxvsHWtP3QD15NQQmb2xER0WNmejXo5T8FeNMlIiLjkZqaipKSEri4uOgsd3FxQXJyst51YmJicOTIEahUKmzbtg2pqamYPn060tPTdfqhax0/fhwXLlxAeHi4tMzW1hbdunXD4sWLERgYCBcXF2zcuBF//PEHWrRoIZUbMmQIRo8eDR8fH8TGxmL+/Pno378/Tp06BaVSWWFfGo0GGo1Gep6dnV3jc1IfQnztYasyQ1peIc4lZKKTt72hQyIiIhNicjXoMvZBJyIiIyZ7YL5QIUSFZVqlpaWQyWTYsGEDOnfujKFDh2LZsmVYu3at3lr08PBwBAUFoXPnzjrLv/32Wwgh4OHhAaVSif/+978YP348FAqFVGbs2LEYNmwYgoKCMHz4cPz666+Ijo7Gzp079cYWFhYGtVotPby8vGp6KuqFuUKOPgFNAQD7LxtfCwAiIjJuJpigl/1kE3ciIjImTk5OUCgUFWrLU1JSKtSqa7m5ucHDwwNqtVpaFhgYCCEEEhISdMrm5+dj06ZNUv/2+/n7+yMyMhK5ubmIj4/H8ePHUVRUBD8/v0rjdXNzg4+PD65du6b39blz5yIrK0t6xMfHV7qtx21AYNlo7vs4HzoRET1mJpegy8sTdM6DTkRExsTCwgLBwcGIiIjQWR4REYHu3bvrXadHjx64desWcnNzpWXR0dGQy+Xw9PTUKfv9999Do9HghRdeqDQGa2truLm5ISMjA3v27MGIESMqLZuWlob4+Hi4uekfaE2pVMLOzk7n0VD0CXCGXAZcTsrGrczK++sTERHVNZNL0GXlvdCZnhMRkbGZPXs2/ve//2H16tW4fPkyZs2ahbi4OEybNg1AWa30hAkTpPLjx4+Ho6MjJk+ejEuXLuHQoUN48803MWXKFFhaWupsOzw8HCNHjoSjo2OF/e7Zswe7d+9GbGwsIiIi0K9fP7Rs2RKTJ08GAOTm5mLOnDmIiorCjRs3cPDgQQwfPhxOTk54+umn6/GM1A8Hawup7/l+1qITEdFjZHKDxMm1X0mwBp2IiIzM2LFjkZaWhkWLFiEpKQlBQUHYtWsXfHx8AABJSUmIi4uTytvY2CAiIgIzZsxASEgIHB0dMWbMGJ2p0YCyWvUjR45g7969eveblZWFuXPnIiEhAQ4ODnj22WexZMkSmJubAwAUCgXOnz+PdevWITMzE25ubujXrx82b94MW1vbejob9WtAoAtO3szA/ispeKGrj6HDISIiEyETRjAheHZ2NtRqNbKysmrdBO7VDaex83wSFv6tDSZ2962bAImIyOTU5b2JGt75jL6dg9D/HILSTI6z74fC0kLx8JWIiKhRMcS9yfSauEuDxDX47yWIiIjIQFo428DT3hKa4lL8fj3V0OEQEZGJMMEEndOsERERUdVkMhn6tyobzf1gNPuhExHR42FyCbp2FHfm50RERFQV7Xzoh6JZg05ERI+HySXo5fk5m7gTERFRlbo2c4S5Qoa49HzcSM0zdDhERGQCTC9BL2/izvyciIiIqmKtNEOIjwMAIDL6joGjISIiU2CCCXrZz1Jm6ERERPQQvaVm7kzQiYio/plegl7eyJ3pORERET1M7wAnAEBUTBoKi0sNHA0RETV2Jpegy1mDTkRERNUU6GoHJxsl8gtLcPJmuqHDISKiRs7kEvR786AbNg4iIiJq+ORyGXq3KKtF52juRERU30wuQZdLg8QxQyciIqKHYz90IiJ6XEwuQWcNOhEREdVEz/Ia9EtJ2biTozFwNERE1JiZYIJelqGXMkEnIiKianCyUSLIww4AcPgaa9GJiKj+mF6CXv5TcBx3IiIiqqbeLdjMnYiI6p/JJehy1qATERFRDWn7oR++lopSfoggIqJ6YnIJukyqQufNlYiIiKqnk7c9rC0USMsrxKWkbEOHQ0REjZTJJeisQSciIqKasjCTo5t/2WBxkWzmTkRE9aRGCfrKlSvRrl072NnZwc7ODt26dcOvv/5a5TqRkZEIDg6GSqVCs2bNsGrVqloFXFfYB52IiIhqok+Adj50JuhERFQ/apSge3p64uOPP8bJkydx8uRJ9O/fHyNGjMDFixf1lo+NjcXQoUPRq1cvnDlzBu+++y5ef/11bNmypU6CfxSsQSciIqJHoe2HfupmBnI1xQaOhoiIGiOzmhQePny4zvMlS5Zg5cqVOHbsGNq0aVOh/KpVq+Dt7Y3ly5cDAAIDA3Hy5El8+umnePbZZx896lrgPOhERET0KHwcreHjaIWbafmI+isNA1u7GDokIiJqZB65D3pJSQk2bdqEvLw8dOvWTW+ZqKgohIaG6iwbNGgQTp48iaKiokq3rdFokJ2drfOoK3IpQWeGTkRERDXD6daIiKg+1ThBP3/+PGxsbKBUKjFt2jRs27YNrVu31ls2OTkZLi663y67uLiguLgYqample4jLCwMarVaenh5edU0zErJyqvQmZ4TERFRTfUpb+Z+6BoTdCIiqns1TtBbtmyJs2fP4tixY/jHP/6BiRMn4tKlS5WWl0nzmpXR1lw/uPx+c+fORVZWlvSIj4+vaZhVxFP2k3OYEhERUU1183eEuUKGm2n5uJGaZ+hwiIiokalRH3QAsLCwQPPmzQEAISEhOHHiBD777DN89dVXFcq6uroiOTlZZ1lKSgrMzMzg6OhY6T6USiWUSmVNQ6sWGViDTkRERI/GWmmGYB97HItJx6Frd+DrZG3okIiIqBGp9TzoQghoNBq9r3Xr1g0RERE6y/bu3YuQkBCYm5vXdtePhIPEERERUW1oR3NnP3QiIqprNUrQ3333XRw+fBg3btzA+fPnMW/ePBw8eBDPP/88gLKm6RMmTJDKT5s2DTdv3sTs2bNx+fJlrF69GuHh4ZgzZ07dHkUNaAeJK2WGTkRERI9AO1Bc1F9pKCwuNXA0RETUmNSoifvt27fx4osvIikpCWq1Gu3atcPu3bsxcOBAAEBSUhLi4uKk8n5+fti1axdmzZqFL7/8Eu7u7vjvf/9rsCnWgHtN3ImIiIgeRWs3OzjZWCA1txCnbmagm3/l3faIiIhqokYJenh4eJWvr127tsKyPn364PTp0zUKqj6xBp2IiIhqQy6XoVeLpth2JhGHrt1hgk5ERHWm1n3QjY52mjXm50RERPSIegc4AWA/dCIiqlsml6CzBp2IiIhqq1d5P/SLt7JxJ0f/YLlEREQ1ZXIJOqdZIyIiotpyslGijbsdAODIddaiExFR3TC5BF0uTbPGFJ2IiIzPihUr4OfnB5VKheDgYBw+fLjK8hqNBvPmzYOPjw+USiX8/f2xevVq6fW+fftCJpNVeAwbNkwqk5OTg5kzZ8LHxweWlpbo3r07Tpw4obMfIQQWLFgAd3d3WFpaom/fvrh48WLdHnwDc2+6tVQDR0JERI2FySXonAediIiM1ebNmzFz5kzMmzcPZ86cQa9evTBkyBCdGVQeNGbMGOzbtw/h4eG4evUqNm7ciFatWkmvb926FUlJSdLjwoULUCgUGD16tFRm6tSpiIiIwLfffovz588jNDQUTz75JBITE6UyS5cuxbJly/DFF1/gxIkTcHV1xcCBA5GTk1M/J6MB0E63dvjaHZSW8oMFERHVngkm6GUZOvugExGRsVm2bBleeuklTJ06FYGBgVi+fDm8vLywcuVKveV3796NyMhI7Nq1C08++SR8fX3RuXNndO/eXSrj4OAAV1dX6REREQErKyspQS8oKMCWLVuwdOlS9O7dG82bN8eCBQvg5+cn7VcIgeXLl2PevHl45plnEBQUhG+++Qb5+fn47rvv6v/EGEiwjz2sLRRIzS3EpaRsQ4dDRESNgAkm6GU/mZ8TEZExKSwsxKlTpxAaGqqzPDQ0FEePHtW7zo4dOxASEoKlS5fCw8MDAQEBmDNnDgoKCirdT3h4OMaNGwdra2sAQHFxMUpKSqBSqXTKWVpa4siRIwCA2NhYJCcn68SmVCrRp0+fSmPTaDTIzs7WeRgbCzO5NMXaoWvsh05ERLVncgm6XKpBN3AgRERENZCamoqSkhK4uLjoLHdxcUFycrLedWJiYnDkyBFcuHAB27Ztw/Lly/Hjjz/i1Vdf1Vv++PHjuHDhAqZOnSots7W1Rbdu3bB48WLcunULJSUlWL9+Pf744w8kJSUBgLT/msQWFhYGtVotPby8vKp3IhqYe/3QmaATEVHtmVyCXl6BDsFx3ImIyAhpu2ppCSEqLNMqLS2FTCbDhg0b0LlzZwwdOhTLli3D2rVr9daih4eHIygoCJ07d9ZZ/u2330IIAQ8PDyiVSvz3v//F+PHjoVAoHjm2uXPnIisrS3rEx8c/9NgbIm0/9FM3M5CnKTZwNEREZOxMLkHX1qCziTsRERkTJycnKBSKCjXSKSkpFWqutdzc3ODh4QG1Wi0tCwwMhBACCQkJOmXz8/OxadMmndpzLX9/f0RGRiI3Nxfx8fE4fvw4ioqK4OfnBwBwdXUFgBrFplQqYWdnp/MwRr5O1vB2sEJRiUDUX2mGDoeIiIycySXoMk6zRkRERsjCwgLBwcGIiIjQWR4REaEz6Nv9evTogVu3biE3N1daFh0dDblcDk9PT52y33//PTQaDV544YVKY7C2toabmxsyMjKwZ88ejBgxAgDg5+cnDTCnVVhYiMjIyEpja0x6BzgBYD90IiKqPRNM0NkHnYiIjNPs2bPxv//9D6tXr8bly5cxa9YsxMXFYdq0aQDKmo1PmDBBKj9+/Hg4Ojpi8uTJuHTpEg4dOoQ333wTU6ZMgaWlpc62w8PDMXLkSDg6OlbY7549e7B7927ExsYiIiIC/fr1Q8uWLTF58mQAZffWmTNn4qOPPsK2bdtw4cIFTJo0CVZWVhg/fnw9npGGQdvMnf3QiYiotswMHcDjdq8POhERkXEZO3Ys0tLSsGjRIiQlJSEoKAi7du2Cj48PACApKUlnTnQbGxtERERgxowZCAkJgaOjI8aMGYMPP/xQZ7vR0dE4cuQI9u7dq3e/WVlZmDt3LhISEuDg4IBnn30WS5Ysgbm5uVTmrbfeQkFBAaZPn46MjAx06dIFe/fuha2tbT2ciYalm78jzOQy3EjLx820PPg4Whs6JCIiMlIyYQRtvbOzs6FWq5GVlVXrPmprfo/Fwp8vYVg7N3w5vlMdRUhERKamLu9NZPznc+xXUfgjNh2LR7TBi918DR0OERHVAUPcm0yuibtc6oRu2DiIiIio8dBOtxYZnWrgSIiIyJiZXIKuzc9LG37DASIiIjISfcoT9Ki/UlFYXGrgaIiIyFiZYILOadaIiIiobrV2s4OjtQXyCktwOi7D0OEQEZGRMr0Evfwna9CJiIiorsjlMvRqUT7dGkdzJyKiR2RyCbq2DzrTcyIiIqpL2n7onA+diIgelckl6NIYcczQiYiIqA71Kp8P/UJiNlJzNQaOhoiIjJHJJehyKUFnhk5ERER1p6mtEq3dyqbhOXKNo7kTEVHNmVyCLgObuBMREVH9kJq5sx86ERE9AtNL0DnNGhEREdWT3gHlA8VdS0VpKT9rEBFRzZhggs5p1oiIiKh+hPg4wMpCgdRcDS4nZxs6HCIiMjIml6DLWYNORERE9cTCTI5uzRwBAIei2Q+diIhqxuQSdG0TdyIiIqL6wH7oRET0qEwuQdfOg84adCIiIqoP2gT95M105GmKDRwNEREZE5NL0LWYnxMREVF98HW0gpeDJYpKBI7FpBk6HCIiMiIml6CzBp2IiIjqk0wmQ+8WbOZOREQ1Z3IJurYPOvNzIiIiqi9SP/RrHCiOiIiqz+QSdDmnWSMiIqJ61t3fEWZyGWJT8xCXlm/ocIiIyEiYXIKuHcRdgBk6ERER1Q9blTk6edsDACKvsZk7ERFVj+kl6NI86IaNg4iIiBq33gFOANgPnYiIqq9GCXpYWBieeOIJ2NrawtnZGSNHjsTVq1erXOfgwYOQyWQVHleuXKlV4I9KJjVxZ4ZORERE9adPgDMAIOqvNBSVlBo4GiIiMgY1StAjIyPx6quv4tixY4iIiEBxcTFCQ0ORl5f30HWvXr2KpKQk6dGiRYtHDro2tE3cWYNORERE9amNux0crS2QqynG6ZsZhg6HiIiMgFlNCu/evVvn+Zo1a+Ds7IxTp06hd+/eVa7r7OyMJk2a1DjAuiYNEmfgOIiIiKhxk8tl6NnCCT+dvYVD1+6gSzNHQ4dEREQNXK36oGdlZQEAHBwcHlq2Y8eOcHNzw4ABA3DgwIEqy2o0GmRnZ+s86sq9adaYohMREVH9ujcfOqdbIyKih3vkBF0IgdmzZ6Nnz54ICgqqtJybmxu+/vprbNmyBVu3bkXLli0xYMAAHDp0qNJ1wsLCoFarpYeXl9ejhlkBp1kjIiKix6VX+UBxF25lIS1XY+BoiIiooatRE/f7vfbaa/jzzz9x5MiRKsu1bNkSLVu2lJ5369YN8fHx+PTTTyttFj937lzMnj1bep6dnV13Sbo0ijszdCIiIqpfzrYqBLrZ4XJSNo5cT8WIDh6GDomIiBqwR6pBnzFjBnbs2IEDBw7A09Ozxut37doV165dq/R1pVIJOzs7nUddYQ06ERERPU7a6dYiOd0aERE9RI0SdCEEXnvtNWzduhX79++Hn5/fI+30zJkzcHNze6R1a+veKO7M0ImIiKj+9Snvh374WirHwCEioirVqIn7q6++iu+++w4//fQTbG1tkZycDABQq9WwtLQEUNY8PTExEevWrQMALF++HL6+vmjTpg0KCwuxfv16bNmyBVu2bKnjQ6kebQ06ERER0eMQ7GsPS3MF7uRocDkpB63d665lIBERNS41StBXrlwJAOjbt6/O8jVr1mDSpEkAgKSkJMTFxUmvFRYWYs6cOUhMTISlpSXatGmDnTt3YujQobWL/BHJ2AediIiIHiOlmQLd/B2x/0oKDl27wwSdiIgqVeMm7voe2uQcANauXYuDBw9Kz9966y1cv34dBQUFSE9Px+HDhw2WnAP3T7NmsBCIiIge2YoVK+Dn5weVSoXg4GAcPny4yvIajQbz5s2Dj48PlEol/P39sXr1aun1vn37QiaTVXgMGzZMKlNcXIz33nsPfn5+sLS0RLNmzbBo0SKUlpZKZSZNmlRhG127dq37E2Ckerco64d+iP3QiYioCo88iruxkpX3Qmd+TkRExmbz5s2YOXMmVqxYgR49euCrr77CkCFDcOnSJXh7e+tdZ8yYMbh9+zbCw8PRvHlzpKSkoLi4WHp969atKCwslJ6npaWhffv2GD16tLTsk08+wapVq/DNN9+gTZs2OHnyJCZPngy1Wo1//vOfUrnBgwdjzZo10nMLC4u6PHyj1jugrB/6yRsZyC8shpWFyX0EIyKiajC5u4OcTdyJiMhILVu2DC+99BKmTp0KoGyclz179mDlypUICwurUH737t2IjIxETEwMHBwcAAC+vr46ZbTLtTZt2gQrKyudBD0qKgojRoyQatV9fX2xceNGnDx5UmddpVIJV1fXWh9nY+TnZA1Pe0skZBTgWEwa+rdyMXRIRETUAD3SNGvGTCa1cTdsHERERDVRWFiIU6dOITQ0VGd5aGgojh49qnedHTt2ICQkBEuXLoWHhwcCAgIwZ84cFBQUVLqf8PBwjBs3DtbW1tKynj17Yt++fYiOjgYAnDt3DkeOHKnQZe3gwYNwdnZGQEAAXn75ZaSkpDzq4TY6MplMqkU/FJ1q4GiIiKihYg06ERGREUhNTUVJSQlcXHRrXl1cXKRZVR4UExODI0eOQKVSYdu2bUhNTcX06dORnp6u0w9d6/jx47hw4QLCw8N1lr/99tvIyspCq1atoFAoUFJSgiVLluC5556TygwZMgSjR4+Gj48PYmNjMX/+fPTv3x+nTp2CUqmssC+NRgONRiM9z87OrtH5MEa9WzTFd3/EsR86ERFVyuQSdFagExGRMZM9MF2oEKLCMq3S0lLIZDJs2LABarUaQFkz+VGjRuHLL7+UpkjVCg8PR1BQEDp37qyzfPPmzVi/fj2+++47tGnTBmfPnsXMmTPh7u6OiRMnAgDGjh0rlQ8KCkJISAh8fHywc+dOPPPMMxViCwsLw8KFC2t+AoxY9+aOUMhliEnNQ3x6PrwcrAwdEhERNTAm28SdNehERGRMnJycoFAoKtSWp6SkVKhV13Jzc4OHh4eUnANAYGAghBBISEjQKZufn49NmzZJ/dvv9+abb+Kdd97BuHHj0LZtW7z44ouYNWuW3n7v9+/bx8cH165d0/v63LlzkZWVJT3i4+Mr3VZjYacyRyfvJgCASNaiExGRHqaXoJf/ZH5ORETGxMLCAsHBwYiIiNBZHhERge7du+tdp0ePHrh16xZyc3OlZdHR0ZDL5fD09NQp+/3330Oj0eCFF16osJ38/HzI5bofGRQKhc40aw9KS0tDfHw83Nzc9L6uVCphZ2en8zAFfcr7oTNBJyIifUwvQS+vQWeCTkRExmb27Nn43//+h9WrV+Py5cuYNWsW4uLiMG3aNABltdITJkyQyo8fPx6Ojo6YPHkyLl26hEOHDuHNN9/ElClT9DZvHzlyJBwdHSvsd/jw4ViyZAl27tyJGzduYNu2bVi2bBmefvppAEBubi7mzJmDqKgo3LhxAwcPHsTw4cPh5OQklaEyfVs6AwB+v54KTXGJgaMhIqKGxuT6oGsHiRPM0ImIyMiMHTsWaWlpWLRoEZKSkhAUFIRdu3bBx8cHAJCUlIS4uDipvI2NDSIiIjBjxgyEhITA0dERY8aMwYcffqiz3ejoaBw5cgR79+7Vu9/PP/8c8+fPx/Tp05GSkgJ3d3e88soreP/99wGU1aafP38e69atQ2ZmJtzc3NCvXz9s3rwZtra29XQ2jFMbdzs42yqRkqPBidgM9GzhZOiQiIioAZEJI8hUs7OzoVarkZWVVesmcOcTsjD8iyNwtVPh2LsD6ihCIiIyNXV5byLTOp9v/XgO359MwEs9/TD/qdaGDoeIiCphiHuTCTZxL/spOI47ERERGUC/8mbuB65wnngiItJlsgl6KfNzIiIiMoAeLZxgVj7d2o3UPEOHQ0REDYjJJehyDhJHREREBmSnMscTvg4AgINXWYtORET3mFyCLuMgcURERGRg/VqVTbe2/yqnWyMiontMLkGXatANHAcRERGZrv6tyvqhH4tJQ35hsYGjISKihsLkEvTyCnSUsgadiIiIDMS/qQ087S1RWFyKo9fTDB0OERE1EKaXoLMPOhERERmYTCaTatH3cTR3IiIqZ4IJetlP1qATERGRIQ0IdAEA7L9ym2PjEBERABNM0OX3JkInIiIiMpgufg6wslDgdrYGF29lGzocIiJqAEwuQWcfdCIiImoIVOYK9GzuBADYz2buREQEE0zQOYo7ERERNRQDAsv7oV++beBIiIioITC5BJ190ImIiKih6NeyLEE/l5CFlJy7Bo6GiIgMzWQTdObnREREZGjOdiq081QDAA6wmTsRkckzwQSdTdyJiIio4RjQqmw094hLTNCJiEyd6SXo5T85nQkRERE1BANblyXoR67fQUFhiYGjISIiQzK5BF0aJI75ORERETUAgW628GhiibtFpTh87Y6hwyEiIgMyuQSdg8QRERFRQyKTyRDapqwWfe8ljuZORGTKTDZBZ3pOREREDYW2mfu+y7dRXFJq4GiIiMhQTC9BB5u4ExERUcPS2dcBaktzZOQX4dTNDEOHQ0REBmJyCbpcdu/fHCiOiIiIGgIzhRwDWpXNiR7BZu5ERCbL5BJ07TRrAFDK/JyIiIgaiPv7obMSgYjINJlcgs4adCIiImqIerVoCgszOeLS83H1do6hwyEiIgMwuQRdBtagExERUcNjrTRDr+ZOAICIi2zmTkRkimqUoIeFheGJJ56Ara0tnJ2dMXLkSFy9evWh60VGRiI4OBgqlQrNmjXDqlWrHjng2pLdd8SCY7kTERFRA6IdzZ3TrRERmaYaJeiRkZF49dVXcezYMURERKC4uBihoaHIy8urdJ3Y2FgMHToUvXr1wpkzZ/Duu+/i9ddfx5YtW2od/KO4r4U7R3InIiKiBmVAoAtkMuB8YhZuZRYYOhwiInrMzGpSePfu3TrP16xZA2dnZ5w6dQq9e/fWu86qVavg7e2N5cuXAwACAwNx8uRJfPrpp3j22WcfLepakN83SBwTdCIiImpImtoq0cnbHqduZiDi0m1M7O5r6JCIiOgxqlUf9KysLACAg4NDpWWioqIQGhqqs2zQoEE4efIkioqK9K6j0WiQnZ2t86gr9+XnKGWGTkRERA3M4DauAICd55MMHAkRET1uj5ygCyEwe/Zs9OzZE0FBQZWWS05OhouLi84yFxcXFBcXIzU1Ve86YWFhUKvV0sPLy+tRw6xApwa9zrZKREREVDeGtnMDAJy4kY7b2XcNHA0RET1Oj5ygv/baa/jzzz+xcePGh5a9f+5x4N70Zg8u15o7dy6ysrKkR3x8/KOGWSXWoBMREVFD49HEEp28m0AIYBdr0YmITMojJegzZszAjh07cODAAXh6elZZ1tXVFcnJyTrLUlJSYGZmBkdHR73rKJVK2NnZ6Tzqik4NemmdbZaIiIiozgxr5w4A2PknE3QiIlNSowRdCIHXXnsNW7duxf79++Hn5/fQdbp164aIiAidZXv37kVISAjMzc1rFm0dkN9Xac9p1oiIyNisWLECfn5+UKlUCA4OxuHDh6ssr9FoMG/ePPj4+ECpVMLf3x+rV6+WXu/bty9kMlmFx7Bhw6QyxcXFeO+99+Dn5wdLS0s0a9YMixYtQmnpvW+6hRBYsGAB3N3dYWlpib59++LixYt1fwJMxNC2Zf3QT97MQFIWR3MnIjIVNUrQX331Vaxfvx7fffcdbG1tkZycjOTkZBQU3LtxzJ07FxMmTJCeT5s2DTdv3sTs2bNx+fJlrF69GuHh4ZgzZ07dHUUN3F+DXlLKBJ2IiIzH5s2bMXPmTMybNw9nzpxBr169MGTIEMTFxVW6zpgxY7Bv3z6Eh4fj6tWr2LhxI1q1aiW9vnXrViQlJUmPCxcuQKFQYPTo0VKZTz75BKtWrcIXX3yBy5cvY+nSpfjXv/6Fzz//XCqzdOlSLFu2DF988QVOnDgBV1dXDBw4EDk5OfVzMho5N7UlnvC1B8BadCIiU1KjadZWrlwJoOzb9vutWbMGkyZNAgAkJSXpfFDw8/PDrl27MGvWLHz55Zdwd3fHf//7X4NMsQYA8vuq0JmfExGRMVm2bBleeuklTJ06FQCwfPly7NmzBytXrkRYWFiF8rt370ZkZCRiYmKkGVd8fX11yjw4E8umTZtgZWWlk6BHRUVhxIgRUq26r68vNm7ciJMnTwIoqz1fvnw55s2bh2eeeQYA8M0338DFxQXfffcdXnnllbo5ASZmWFs3nLiRgZ3nkzC1VzNDh0NERI9BjZu463tok3MAWLt2LQ4ePKizXp8+fXD69GloNBrExsZi2rRpdRH7I9Pm6BwkjoiIjEVhYSFOnTpVYerS0NBQHD16VO86O3bsQEhICJYuXQoPDw8EBARgzpw5Oi3fHhQeHo5x48bB2tpaWtazZ0/s27cP0dHRAIBz587hyJEjGDp0KAAgNjYWycnJOrEplUr06dOn0tjqc0rVxmJoWzfIZMCZuEwkZOQbOhwiInoMalSD3lgo5DKUlggm6EREZDRSU1NRUlKid+rSBwdj1YqJicGRI0egUqmwbds2pKamYvr06UhPT9fph651/PhxXLhwAeHh4TrL3377bWRlZaFVq1ZQKBQoKSnBkiVL8NxzzwGAtH99sd28eVNvbGFhYVi4cGH1Dt5EOdup0NnXAX/EpmPX+ST8vbe/oUMiIqJ69sjTrBkz7fRu7INORETGRt/UpZVNW1paWgqZTIYNGzagc+fOGDp0KJYtW4a1a9fqrUUPDw9HUFAQOnfurLN88+bN0hg0p0+fxjfffINPP/0U33zzzSPH9rimVDV2T5XPic5+6EREpsEkE3RF+YcFVqATEZGxcHJygkKh0Dt16YM111pubm7w8PCAWq2WlgUGBkIIgYSEBJ2y+fn52LRpk9S//X5vvvkm3nnnHYwbNw5t27bFiy++iFmzZkn93l1dy0Ycr0ls9TmlamMyOMgNchlwLiELcWls5k5E1NiZZIKu7YPOGnQiIjIWFhYWCA4OrjB1aUREBLp37653nR49euDWrVvIzc2VlkVHR0Mul8PT01On7Pfffw+NRoMXXnihwnby8/Mhl+t+ZFAoFNI0a35+fnB1ddWJrbCwEJGRkZXGRtXT1FaJrs0cAQA7z7MWnYiosTPNBL08Q2cfdCIiMiazZ8/G//73P6xevRqXL1/GrFmzEBcXJw2++uBUp+PHj4ejoyMmT56MS5cu4dChQ3jzzTcxZcoUWFpa6mw7PDwcI0eOhKOjY4X9Dh8+HEuWLMHOnTtx48YNbNu2DcuWLcPTTz8NoKxp+8yZM/HRRx9h27ZtuHDhAiZNmgQrKyuMHz++Hs+IaRimbeZ+/paBIyEiovpmkoPEaedCZ4JORETGZOzYsUhLS8OiRYuQlJSEoKAg7Nq1Cz4+PgAqTnVqY2ODiIgIzJgxAyEhIXB0dMSYMWPw4Ycf6mw3OjoaR44cwd69e/Xu9/PPP8f8+fMxffp0pKSkwN3dHa+88gref/99qcxbb72FgoICTJ8+HRkZGejSpQv27t0LW1vbejgTpmVwG1e8/9NFXEjMxo3UPPg6WT98JSIiMkoyIRp+lpqdnQ21Wo2srKw66aPWaXEE0vMKsXdWbwS48IMDERHVXF3fm0wdz2fVXgz/A4evpeLNQS3xar/mhg6HiMgkGOLeZJpN3NkHnYiIiIyIdjT3XziaOxFRo2aiCTqbuBMREZHxCG3tCjO5DJeTsvHXndyHr0BEREbJtBP0UgMHQkRERFQN9tYW6NHcCQDnRCciasxMMkFXcBR3IiIiMjL3mrnfghEMIURERI/AJBN0mbYPOm9uREREZCRC27jCQiFH9O1cXLyVbehwiIioHphkgq6tQee3z0RERGQs1JbmGNjGBQDw46kEA0dDRET1wSQTdG0f9BL2QSciIiIjMjrYEwCw/WwiNMUlBo6GiIjqmokm6GU/2QediIiIjEmvFk3hYqdEZn4R9l9OMXQ4RERUx0w0QdeO4s4EnYiIiIyHQi7DM53KatF/YDN3IqJGxyQT9HujuBs4ECIiIqIaGlXezD0y+g5Ssu8aOBoiIqpLJpmgS33Q2cSdiIiIjIx/Uxt08m6CklKBbWcSDR0OERHVIdNM0MuPmn3QiYiIyBiNDvECUDaaO2elISJqPEwyQVewDzoREREZsWHt3KAyl+NaSi7OJWQZOhwiIqojJpmgy2Tsg05ERETGy05ljsFtXAEAP56KN3A0RERUV0wyQdcOElfCDJ2IiIiM1KjgsmbuO87ewt0izolORNQYmGSCrp0HnX22iIiIyFh193eEu1qF7LvFiLh029DhEBFRHTDRBJ1N3ImIiMi4yeUyPFs+5dr3J9nMnYioMTDpBJ3TrBEREZExG1M+mvvha6mIS8s3cDRERFRbJpmga/ugs4k7ERERGTMvByv0DmgKANjwx00DR0NERLVlkgl6eQU6B4kjIiIiozehqw8AYNOJeBQUcrA4IiJjZpIJurYGnfk5ERERGbt+rZzh5WCJrIIi/HQ20dDhEBFRLZhkgi4NEscMnYiIiIycQi7DhK6+AIC1R2+wCx8RkREz7QSdNzAiIiJqBMaEeMHSXIEryTn4Izbd0OEQEdEjMtEEvewnR3EnIiKixkBtZY6RHT0AAN8cvWHYYIiI6JGZZILOPuhERETU2EzsXjZY3N5Lt3Ers8DA0RAR0aMwyQSdfdCJiIiosWnlaoeuzRxQUiqw/hinXCMiMkY1TtAPHTqE4cOHw93dHTKZDNu3b6+y/MGDByGTySo8rly58qgx15pczj7oRERE1PhM6u4LAFgXdZNTrhERGSGzmq6Ql5eH9u3bY/LkyXj22Wervd7Vq1dhZ2cnPW/atGlNd11n5JwHnYiIiBqhJwNd4K5W4VbWXQS+vxu9WjghxMcBT/jZo6OXPSwtFIYOkYiIqlDjBH3IkCEYMmRIjXfk7OyMJk2a1Hi9+qAob+LOBJ2IiIgaEzOFHG8NboV3t51HfmEJDl9LxeFrqQAAc4UMQR5qdPZzQFc/Rzzh5wAbZY0/ChIRUT16bH3QO3bsCDc3NwwYMAAHDhyosqxGo0F2drbOoy6ZK8oOu5gJOhERGZkVK1bAz88PKpUKwcHBOHz4cJXlNRoN5s2bBx8fHyiVSvj7+2P16tXS63379tXbFW3YsGFSGV9fX71lXn31VanMpEmTKrzetWvXuj8B9FAjO3rgzw9C8cuMnlg8og2Gt3eHi50SRSUCZ+Iy8VVkDCavPYH2C/dixJe/45PdV3D42h02iSciagDq/WtTNzc3fP311wgODoZGo8G3336LAQMG4ODBg+jdu7fedcLCwrBw4cJ6i8lMUVaDXlRSWm/7ICIiqmubN2/GzJkzsWLFCvTo0QNfffUVhgwZgkuXLsHb21vvOmPGjMHt27cRHh6O5s2bIyUlBcXFxdLrW7duRWFhofQ8LS0N7du3x+jRo6VlJ06cQEnJveTtwoULGDhwoE4ZABg8eDDWrFkjPbewsKj1MdOjMVPIEeShRpCHGi9284UQAgkZBfgjNh3HY9NwLCYdcen5OBefiXPxmVh58C9YKOQI9rFHzxZO6NHcCW091NLMN0RE9HjIhHj0kdJkMhm2bduGkSNH1mi94cOHQyaTYceOHXpf12g00Gg00vPs7Gx4eXkhKytLpx/7o1qw4yLWHr2B1/o1x5xBLWu9PSIiMj3Z2dlQq9V1dm+qji5duqBTp05YuXKltCwwMBAjR45EWFhYhfK7d+/GuHHjEBMTAwcHh2rtY/ny5Xj//feRlJQEa2trvWVmzpyJX375BdeuXYOsvNvYpEmTkJmZ+dDBYytjiPNp6hIzCxD1Vxqi/krD0b9SkZR1V+d1O5UZejR3Qq8WTdGrhRO8HKwMFCkRkWEY4t5kkI5HXbt2xfr16yt9XalUQqlU1tv+zcq/DS4qZQ06EREZh8LCQpw6dQrvvPOOzvLQ0FAcPXpU7zo7duxASEgIli5dim+//RbW1tb429/+hsWLF8PS0lLvOuHh4Rg3blylyXlhYSHWr1+P2bNnS8m51sGDB6UxZ/r06YMlS5bA2dn5EY6WHgePJpYYFeyJUcGeEEIgNjUPv19PxZHrqTj6Vxqy7xbj1wvJ+PVCMgDAz8kafQKaok9AU3Rp5gArC/ZfJyKqawa5sp45cwZubm6G2DUAwNysrA96UTH7oBMRkXFITU1FSUkJXFxcdJa7uLggOTlZ7zoxMTE4cuQIVCoVtm3bhtTUVEyfPh3p6ek6/dC1jh8/jgsXLiA8PLzSOLZv347MzExMmjRJZ/mQIUMwevRo+Pj4IDY2FvPnz0f//v1x6tQpvV+662stR4Yjk8nQrKkNmjW1wYvdfFFSKvBnQmb5IHN3cDouE7GpeYhNzcPaozdgoZCjs58Degc4oW9LZ7RwtqnwhQ0REdVcjRP03NxcXL9+XXoeGxuLs2fPwsHBAd7e3pg7dy4SExOxbt06AGVN5Xx9fdGmTRvpW/ctW7Zgy5YtdXcUNWReXoNezBp0IiIyMg8mQUKIShOj0tJSyGQybNiwAWq1GgCwbNkyjBo1Cl9++WWFWvTw8HAEBQWhc+fOle4/PDwcQ4YMgbu7u87ysWPHSv8OCgpCSEgIfHx8sHPnTjzzzDMVtlPf481Q7SjkMnT0tkdHb3u8PqAFcu4W4ffraTh07Q4ir95BYmYBjpTXtn+06wo8mliid0BT9G3ZFD2aO3F0eCKiR1Tjq+fJkyfRr18/6fns2bMBABMnTsTatWuRlJSEuLg46fXCwkLMmTMHiYmJsLS0RJs2bbBz504MHTq0DsJ/NGblo7gXlbAGnYiIjIOTkxMUCkWF2vKUlJQKtepabm5u8PDwkJJzoKzPuhACCQkJaNGihbQ8Pz8fmzZtwqJFiyqN4ebNm/jtt9+wdevWh8br5uYGHx8fXLt2Te/rc+fOlT5DAPfGm6GGyVZljsFBrhgc5AohBGJS8xB59Q4io+/gWEwaEjMLsPF4HDYej4O5QoYQHwf0bdkU/Vqxdp2IqCZqnKD37dsXVY0rt3btWp3nb731Ft56660aB1aftKO4F3MUdyIiMhIWFhYIDg5GREQEnn76aWl5REQERowYoXedHj164IcffkBubi5sbGwAANHR0ZDL5fD09NQp+/3330Oj0eCFF16oNIY1a9bA2dlZZwq2yqSlpSE+Pr7SLm31Pd4M1R+ZTAb/pjbwb2qDKT39UFBYgmMxaTh4NQUHo+/gZlo+omLSEBWThrBfy2rX+7Vqir4Bzuje3JF914mIqmCSV0iL8hr0QiboRERkRGbPno0XX3wRISEh6NatG77++mvExcVh2rRpAFChm9n48eOxePFiTJ48GQsXLkRqairefPNNTJkyRW/z9pEjR8LR0VHvvktLS7FmzRpMnDgRZma6Hx9yc3OxYMECPPvss3Bzc8ONGzfw7rvvwsnJSefLBGqcLC0U6NfKGf1alQ0IGJuaV5asX71Xu77+WBzWH4uDhUKOLs0c0LelM/q2bIpmTtasXSciuo9JJujW5f2ifr2QjKXFJVCaKQwcERER0cONHTsWaWlpWLRoEZKSkhAUFIRdu3bBx8cHACp0M7OxsUFERARmzJiBkJAQODo6YsyYMfjwww91thsdHY0jR45g7969le77t99+Q1xcHKZMmVLhNYVCgfPnz2PdunXIzMyEm5sb+vXrh82bN8PW1raOjp6MhZ+TNfyc/DC5hx/yC4sR9VcaDl69g4PRKYhPLygfeC4Vi38BvB2syprCt3RG12aOsLTgZzIiMm21mgf9canr+ed+PncLMzaeAQBM6OaDRSOCar1NIiIyLZy3u27xfDZ+2r7rB66kIDL6Dv6ISddpzag0k6NLM0f0KR9sjrXrRGRoJjMPuqHZqO4d9rqom0zQiYiIiOrZ/X3Xp/ZqhjxNMY7+lSY1h0/MLMCh6Ds4FH0Hi38BvBwsy5J19l0nIhNiklc6fycbnecaNnMnIiIieqyslWYY2NoFA1u7QAiBaym5iCxvCn8iNgPx6bp915/ws0fvFk3RO6ApWrnasnadiBolk2ziDgCnbmbg2ZVHAQB7Z/VGgAv7yBERUfWxSXbd4vmk++VpyvquR0bfwYGrKUjIKNB53dlWiV4tmqJ3gBN6tWgKB2sLA0VKRI0Zm7g/RsE+9mjjboeLt7IRn57PBJ2IiIiogbBWmuHJ1i54srx2/a87eTh8raz5+7GYdKTkaLDldAK2nE6ATAYEuavRs4UTejV3QrCvPVtGEpHRMtkEHQC87K2kBJ2IiIiIGh6ZTIbmzjZo7myDyT38oCkuwckbGTh4NQWHr6XiSnIOzidm4XxiFlYe/Asqczme8HVArxZO6Nm8KQLd2ByeiIyHaSfoDmVzwMY/0GyKiIiIiBompZkCPZo7oUdzJwBASvZdHLleNnXbkeupuJOjkaZyA67A0doC3fwdy9bxd4K3o5VhD4CIqAomnqCXXaBZg05ERERknJztVHimkyee6eQpDTZ3+Foqjlwraw6flleIX/5Mwi9/JgEAPO0t0d3fEd39ndDN3xEudioDHwER0T2mnaDblyforEEnIiIiMnoymQwBLrYIcLHFSz39UFhcirPxmTj6VyqOXk/DmfgMJGQU4PuTCfj+ZAIAoFlTa3Rt5ohuzRzRtZkjmtoqDXwURGTKTDtBL2/inpCeDyEE+ycRERERNSIWZnJ09nNAZz8HzHyybHT44zfSceyvNPz+Vyou3cpGzJ08xNzJw3d/xAEAmjvboIufA7o0c0RXPwc4s4adiB4jk07QPctr0HM0xcgqKEITK07RQURERNRYWSvN0K+lM/q1dAYAZOUX4fiNdET9lYaomDRcSc7G9ZRcXE/JxYbyhN3PyRpd/BzwhK8DujRzkD4/EhHVB5NO0FXmCjS1VeJOjgbx6QVM0ImIiIhMiNrKHANbu2BgaxcAQGZ+IY7HpiMqJg3HY9NxKSkbsal5iE3Nw6YT8QAAjyaWeMLXHk/4OaCzrwOaO9uwFSYR1RmTTtABwMvesixBz8hHW0+1ocMhIiIiIgNpYmWB0DauCG3jCqCshv3kzXQcj03Hsdh0XEjMQmJmARLPFmD72Vvl65gjxMcBwT72eMLXHm091ZyHnYgeGRN0ByucjstEHEdyJyIiIqL7qK3MMSDQBQMCy2rY8zTFOBOXieM30nEiNh1n4jOQmV+E3y7fxm+XbwMo6/fezkONYB976eFow4HniKh6TD5Bd1OXDRR3O/uugSMhIiIioobMWmmGni2c0LNF2RzsRSWluJCYhVM3M3DiRjpO3cxAam4hTt7MwMmbGdJ6fk7W6ORtjxDfsoS9eVMbyOVsFk9EFZl8gu5oXdbvPD2v0MCREBEREZExMVfI0dHbHh297TG1VzMIIXAjLR8ny5P1UzczcC0lV+rHvuV02dRutiozdPS2RyfvJujkbY/2Xk2gtjQ38NEQUUNg8gm6AxN0IiIiIqoDMpkMfk7W8HOyxugQLwBlA8+dic/EqRsZOHkzHWfjM5FztxiHou/gUPSd8vWAFs426ORtj07e9ujo3QT+rGUnMklM0G2YoBMRERFR/WhiZaEztVtxSSmuJOfgdFxZDfuZ8rGQom/nIvp2rjRavF15LXvH8lr2Dt5NYKdiLTtRY2fyCTqbuBMRERHR42KmkCPIQ40gDzUmdPMFANzJ0eBMXAZOx2XidFwGzidkIftuMSKj7yDyvlp2/6Y2CPa2RyefsqSdtexEjY/JJ+j25XOfp+UVQgjBeSyJiIiI6LFqaqvUmd6tqKQUV5LKatm1iXtcej6up+TiekouNp8sq2W3VZmhg1cTdPQuG3yuA/uyExk9k0/QHcubuBcWlyK/sATWSpM/JURERERkQOYKOdp6qtHWU42J3X0BAKm5GpyNy8SpuAycvpmBPxOykHO3GIevpeLwtVRp3ebONgjxYV92ImNl8tmolYUZVOZy3C0qRXpeIRN0IiIiImpwnGyUeLK1C55sXTYnu7Yvu7aG/dTNDJ1a9gf7sncqr2Xv6N2En3eJGjD+dQJwsLLAray7yMgvhJeDlaHDISIiIiKq0v192V/sVrYsNVeDM3GZOHkzHWfiMvX2ZZfLgAAXW4T4amvZ7eHraMVunkQNBBN0APbWZQk6B4ojIiIiImPlZKPEwNYuGFhey67ty37qZrpUy56YWYAryTm4kpyD9cfiAJQNmtyxfPC5EB8HtPNUQ2WuMOShEJksJui4Nxd6Rj4TdCIiIiJqHO7vyz6pR9my29l3cepm2RRvZ+MzcT4xC2l5hfjt8m38dvk2AMBMLkMbdzt08ilrFh/sYw83taUBj4TIdDBBx72R3NPzigwcCRERERFR/XGxU2FoWzcMbesGANAUl+BCYhbOlNewn7yZgTs5GpxLyMK5hCys+f0GAMBNrSpL2L3t0cnHHq3d7GBhJjfgkRA1TkzQcV8NOpu4ExEREZEJUZopEOzjgGAfB0ztBQghkJBRgNNxGVJN+5XkHCRl3cXOP5Ow888kAICFmRztPNQI9ilL2Dt6N4GzrcrAR0Nk/Jig474adDZxJyIiIiITJpPJ4OVgBS8HK4zo4AEAyC8sxrn4LJwqH3zuVFwGMvOLcLK8xl3L28GqLGH3boJOPvZo6WILMwVr2Ylqggk6AHtrcwCsQSciIiIiepCVhRm6+Tuim78jgLJa9htp+Th5I13qy371dg7i0vMRl56PbWcSAQA2SjN08GqCjt5N0Kl8qje1lbkhD4WowWOCjvv7oDNBJyIiIiKqikwmg5+TNfycrDE6xAsAkFVQhHPxZf3YT8dl4GxcJnI0xThyPRVHrqdK67ZwtkFH77LR4jv5NEEzJxvI5ZzijUiLCTo4ijsRERERUW2oLc3RO6Apegc0BQCUlApcTc7B6bgMnInLxJm4DMSk5uFaSi6upeTi+5MJ0nodvMpq2IN97NHBuwlslExRyHTVuFPIoUOHMHz4cLi7u0Mmk2H79u0PXScyMhLBwcFQqVRo1qwZVq1a9Six1huO4k5ERMZixYoV8PPzg0qlQnBwMA4fPlxleY1Gg3nz5sHHxwdKpRL+/v5YvXq19Hrfvn0hk8kqPIYNGyaV8fX11Vvm1VdflcoIIbBgwQK4u7vD0tISffv2xcWLF+v+BBCRUVDIZWjtbocXuvrg32PaY/+cvjj13pP4vwkhmNbHH539HKAylyOroAiR0Xfwn9+i8UL4H2i3YA+GfHYYc7eex4+nEvDXnVwIIQx9OESPTY2/nsrLy0P79u0xefJkPPvssw8tHxsbi6FDh+Lll1/G+vXr8fvvv2P69Olo2rRptdZ/HO6vQRdCQCZjMxsiImp4Nm/ejJkzZ2LFihXo0aMHvvrqKwwZMgSXLl2Ct7e33nXGjBmD27dvIzw8HM2bN0dKSgqKi4ul17du3YrCwnstyNLS0tC+fXuMHj1aWnbixAmUlJRIzy9cuICBAwfqlFm6dCmWLVuGtWvXIiAgAB9++CEGDhyIq1evwtbWti5PAxEZKUcbJQa2dsHA1i4AgKKSUlxOysbpmxk4XT7NW2JmAS4nZeNyUjY2Ho8DANhbmaOTd9lI8Z187NHeswmsWctOjZRM1OIrKZlMhm3btmHkyJGVlnn77bexY8cOXL58WVo2bdo0nDt3DlFRUdXaT3Z2NtRqNbKysmBnZ/eo4VbqblEJWs3fDQA490Eo1JYcvIKIiKpW3/cmfbp06YJOnTph5cqV0rLAwECMHDkSYWFhFcrv3r0b48aNQ0xMDBwcHKq1j+XLl+P9999HUlISrK2t9ZaZOXMmfvnlF1y7dg0ymQxCCLi7u2PmzJl4++23AZTV3Lu4uOCTTz7BK6+88tD9GuJ8ElHDczv7Lk7fzMCZ+EycvpmBPxOzUFhcqlNGLgNautoh2Ode03hvBytWslGdM8S9qd6/eoqKikJoaKjOskGDBiE8PBxFRUUwN6+YDGs0Gmg0Gul5dnZ2vcaoMlfA2kKBvMISZOQVMkEnIqIGp7CwEKdOncI777yjszw0NBRHjx7Vu86OHTsQEhKCpUuX4ttvv4W1tTX+9re/YfHixbC0tNS7Tnh4OMaNG1dpcl5YWIj169dj9uzZ0ofh2NhYJCcn69zvlUol+vTpg6NHj+pN0B/3vZ6IjIOLnQpD2rphSFs3AEBhcSku3srC6fJ+7KdvZuBW1l2pln39sbJadkdrC3Qsr2UP9rFHO081rCxYy07Gp97ftcnJyXBxcdFZ5uLiguLiYqSmpsLNza3COmFhYVi4cGF9h6bD3toCeYUFSMsrhK+T/g8lREREhpKamoqSkhK999Tk5GS968TExODIkSNQqVTYtm0bUlNTMX36dKSnp+v0Q9c6fvw4Lly4gPDw8Erj2L59OzIzMzFp0iRpmXb/+mK7efOm3u0Y4l5PRMbHwkxennjbA/ADACRn3cXp8mT9dFwGLiRmIy2vEL9dvo3fLt8GUFbLHuhmh47eTdDBq6yW3deRtezU8D2Wr5Ue/EPQtqqv7A9k7ty5mD17tvQ8OzsbXl5e9RcgAGdbJRIyCpCcdbde90NERFQb+u6pld1PS0tLIZPJsGHDBqjVagDAsmXLMGrUKHz55ZcVatHDw8MRFBSEzp07V7r/8PBwDBkyBO7u7rWKzRD3eiJqHFzVKgxt64ah5bXsd4tKcPFWFs7EZZYn7plIzr6Li7eycfHWvVp2eyvzsmTfq6wveztPNWxVbDlLDUu9J+iurq4VvtlPSUmBmZkZHB0d9a6jVCqhVCrrOzQd3g5WOB2Xibj0/Me6XyIioupwcnKCQqHQe099sOZay83NDR4eHlJyDpT1WRdCICEhAS1atJCW5+fnY9OmTVi0aFGlMdy8eRO//fYbtm7dqrPc1dUVQFlN+v0t46qKzRD3eiJqnFTmCgT7OCDY595YG0lZBWUJe3l/9vOJWcjIL8L+KynYfyUFACCTAS1dbNHRuwk6etujk7c9mjlZc152Mqh6T9C7deuGn3/+WWfZ3r17ERISorf/uaF4OVgBAOIzmKATEVHDY2FhgeDgYERERODpp5+WlkdERGDEiBF61+nRowd++OEH5ObmwsbGBgAQHR0NuVwOT09PnbLff/89NBoNXnjhhUpjWLNmDZydnXWmYAMAPz8/uLq6IiIiAh07dgRQ1lc9MjISn3zyySMdLxFRbbipLeHW1lKqZS8sLsWlpGycupmBM+VzsydmFuBKcg6uJOdg4/F4APfmZe9QXsvewbMJ1FYNJ2ehxq/GCXpubi6uX78uPY+NjcXZs2fh4OAAb29vzJ07F4mJiVi3bh2AshHbv/jiC8yePRsvv/wyoqKiEB4ejo0bN9bdUdQBKUFnDToRETVQs2fPxosvvoiQkBB069YNX3/9NeLi4jBt2jQAqHAPHj9+PBYvXozJkydj4cKFSE1NxZtvvokpU6bobd4+cuTISlu3lZaWYs2aNZg4cSLMzHQ/PshkMsycORMfffQRWrRogRYtWuCjjz6ClZUVxo8fXw9ngoioZizM5FLire3LnpJ9t2zwufgMnLmZiT8TM6V52SOj70jrNne2KUvYywehC3CxhYK17FRPapygnzx5Ev369ZOea/uPTZw4EWvXrkVSUhLi4uKk1/38/LBr1y7MmjULX375Jdzd3fHf//63wcyBruVTnqDH3MkzcCRERET6jR07FmlpaVi0aBGSkpIQFBSEXbt2wcfHBwAq3INtbGwQERGBGTNmICQkBI6OjhgzZgw+/PBDne1GR0fjyJEj2Lt3b6X7/u233xAXF4cpU6boff2tt95CQUEBpk+fjoyMDHTp0gV79+7lHOhE1GA526kwOMgVg4PKuulo52U/G5+JM+Wjxt9Iy8f1lFxcT8nFj6cSAADWFgq0L0/2tSPHO9mwyw7VjVrNg/64PI7557LvFqHdgrIPJqfeexKO/CMjIqIqcN7uusXzSUQNUVquRkrYT8dl4Fx8JvIKSyqU87S3lGrYO3rbo7WbHSzM5AaImOpSo5wH3VjYqczRrKk1Yu7k4c/ELPRr6WzokIiIiIiIyIAcbZQYEOiCAYFlA16WlApcS8nB6ZuZOBtf1pf9+p1cJGQUICGjADvO3QJQ1qQ+yN1OqmHv5G0P9yaWVe2KCAATdB0dvJog5k4ejsWkMUEnIiIiIiIdCrkMrVzt0MrVDuO7eAMoa4n7Z3wWTseVDUB3Nj4TGflFOB2XidNxmdK6zrZKqVl8B68maO+lhpUF0zHSxXfEffq1dMbW04nYfzkFc4cEGjocIiIiIiJq4OxU5ujZwgk9WzgBAIQQuJGWL40WfyY+A5eTcpCSo8HeS7ex99JtAIBcBrR0tStrFu/VBB29m6CZkw2neTNxTNDv0zugKRRyGa6l5CIuLR/ejlaGDomIiIiIiIyITCaDn5M1/Jys8Uynsikt8wuLcfFWtpS0n43PRFLWXVxOysblpGx890fZAJ+2SjO0L0/WtaPOc2ws08IE/T5qS3N09nVAVEwafrt8G1N6+hk6JCIiIiIiMnJWFmZ4wtcBT/g6SMuSs+7ibHwGTsdl4mxcJs4nZiFHU4wj11Nx5HqqVM7LwRLtPe+NGt/G3Q4qc4UhDoMeAyboDxgQ6IyomDTsu8IEnYiIiIiI6oerWoXBajcMDnIDABSXlOLq7Rycjc+UBqH7604e4tMLEJ9egF/+TAIAWCjkCHSzRfvyudk7eDWBj6MVZDI2jW8MmKA/YGBrF3y48zL+iElHVkER1Jbmhg6JiIiIiIgaOTOFHG3c1WjjrsbzXXwAAFkFRTifkIUzcRk4l1A23VtaXiHOJWThXEIW1kXdBAA0sTKXatk7eDdBB88msLe2MOTh0CNigv4AH0drNHe2wfWUXBy8moIRHTwMHRIREREREZkgtWXFAeji0vPxZ4J21PhMXLqVjcz8IkRG30Fk9B1pXR9HK6kfe0dvewS62UJpxqbxDR0TdD0GtXHB9ZRc7L6QzASdiIiIiIgaBJlMBh9Ha/g4WmN4e3cAgKa4BFeScnAuoawv+5n4TMSm5uFmWj5upuXjp7Plc7Mr5Ah0t0MHTzXalyfufk7WbBrfwDBB12NIkBu+PPAXDl69g4LCElha8JsmIiIiIiJqeJRmCrT3aoL2Xk0woVvZssz88mbw8Zk6c7Ofi8/EufhMoLxpvNrSvGxO9vuSdo4ab1hM0PVo424HLwdLxKcX4Oc/b2FMiJehQyIiIiIiIqqWJlYW6BPQFH0CmgK41zT+nLY/e3wmLtzKRlZBxabxnvaWZVO9lSf9bT3UHDX+MWKCrodMJsPzXXzw8a9X8Nlv1zA4yBV2Kg4WR0RERERExuf+pvF/K28aX1hciqvJOTgTn4Fz8VnSqPEJGQVIyCjAzvJR4xVyGQJcbMv7s5fVtLdwtoVCzqbx9UEmhBCGDuJhsrOzoVarkZWVBTs7u8eyzzxNMQYtP4SEjAIMauOCFc8H801IREQSQ9ybGjOeTyIiw8sqKMKFxCycjS8bMf5cQibu5GgqlLOyUKCth7qseXz5w12tanT92Q1xb2KCXoVTNzPw3NfHUFhSipd6+mH+U60f276JiKhhY0JZt3g+iYgaHiEEbmXdxfmEssHnzsVn4nxCFvIKSyqUdbJRooOXGu08y5N2TzWaWBn3VG9M0CthyJv2jnO38PrGMwCAl3r64d2hgaxJJyIiJpR1jOeTiMg4lJQK/HUnF2fjM3E2PhN/JmTiSlIOiksrppW+jlZo59kE7TzLatuDjKw/OxP0Shj6ph1+JBaLf7kEAAh0s8Nbg1uiT4umkDNRJyIyWYa+NzU2PJ9ERMbrblEJLt7KxrnypP1cQiZupuVXKKeQy9DSxVaqYS/rz24DM4XcAFE/HBP0SjSEm/bOP5Mwd+ufyL5bDABo1tQaz3byRP9Wzmjlatvo+lsQEVHVGsK9qTHh+SQialy0U739GZ+JP8v7tevrz25prkAbdzupP3s7TzW8HawaRH7FBL0SDeWmnZFXiC8PXMfmE/HI0RRLyx2sLdDes6y/RXNnG/g5WcPH0Qq2HPmdiKjRaij3psaC55OIqHHT9mf/Mz6zLHFPyMSfCVnIvS+v0mpiZX5vEDrPJmjnpYazreqxx8wEvRIN7aadqynGjrO38Nvl2zj6VyruFpXqLWdvZQ5nWxWa2iqlh72VBWxVZvc9zGGjNION0gyWFgpYmMmhNJPDQiFvEN8aERGRfg3t3mTseD6JiExPaXl/9nMJWThX3p/9clIOCksq5leudiq01w5C59kEbT3VUFvWb4UoE/RKNOSbtqa4BJeTcnAuPhMXErMQm5qHG2l5SM0trPW2ldpk3UxR9m9zOczlcijkMulhJpdBXv7z/mX3ysjLyshkUMgBuUwGmUwGuQyQycqey8u/CJA/sFwmk5X/u/w5UL7uveX3ypSvC9l9699bjvKf2u3IZTKg7L+ybZQVgfY7Ce128ODrkJUvg7QM5WW12ypbcm97MujuS+f1+9aF3u3LKuyrwusPbEtLZ/sPvF6t2HXOQcXXtecJOs/xwD8qL3P/9u6P+f51IHtwue66usv07+9+D4u7yu3yCytqYBryvckY8XwSERFQll9dTc6Rkvaz8ZmIuZMLPWPQwdfRCh28mqBt+UB0bdztYGVhVmexMEGvhDHetHPuFiExswB3cjTSIyVHg8z8IuTcLUKuphg5d4t1/q0p1l8TT0QV6XwZIi2TPfBc+/qD32A8/IsF3S8uHvgyo6ovFh7yRYjswY3oXUd/TJXtU2edR/xypspj0rOf+7f/4HLduCr/0qW629ONszpxVb3/v/duhtA2rvoDrgFjvDc1ZDyfRERUmTxNMS4kZuHPhCycSygbhC4+vaBCObkMCHCxRWhrF8wObVnr/Rri3lR3Xy+QDluVOVq5mqNVDT4DCiFQWFIKTXEpNEWl0BSXVPh3SalAcalASWkpSkoh/SwuLXtNeoiyn8UlAqVCu46AEAJCAKUCKBUCony/pfctF+XLS0tF2XNoX7v3s+wbLIHS0nvb0b4uyl/X2X5p2Xak7eu8DqncvXNxb7/3/7tsr2X/k17HvZhx3/ags/371n9gf/q2L22v0u3fiw1VvQ5tGf3b154v7b70vq7n2Mu3KO3vvsW65/G+mBqb+49J6FuoW7qeoyFj8kwdtHAiIiKix8daaYYuzRzRpZmjtCwjr7Bs8Lm4TJxPzMT5xCzcztbgSnIOAlxsDRht7TBBb0BkMhmUZgoozRTA4x8DgUzI/V9OSMsefE16rn294jp4SBmhU0b/dh/8wkHf+pWue/96FdapPO4K233gmPRt/2Ex6a5X2ToPxHz/sjo4z1V/OVN1TA8/Jt24K/lnpeWrek88uBwP2UbF/Tx8m/eXb+PO2lkiIiJjZ29tgT4BTdEnoKm0LDnrLs4lZMLR2sKAkdUOE3QiEyQ1cdbbLJl9vYmIiIjI+LiqVXBV174bmyE1zBnhiYiIiIiIiEwME3QiIiIiIiKiBoAJOhEREREREVEDwASdiIiIiIiIqAFggk5ERERERETUADBBJyIiMiIrVqyAn58fVCoVgoODcfjw4SrLazQazJs3Dz4+PlAqlfD398fq1aul1/v27QuZTFbhMWzYMJ3tJCYm4oUXXoCjoyOsrKzQoUMHnDp1Snp90qRJFbbRtWvXuj14IiKiRo7TrBERERmJzZs3Y+bMmVixYgV69OiBr776CkOGDMGlS5fg7e2td50xY8bg9u3bCA8PR/PmzZGSkoLi4mLp9a1bt6KwsFB6npaWhvbt22P06NHSsoyMDPTo0QP9+vXDr7/+CmdnZ/z1119o0qSJzr4GDx6MNWvWSM8tLIx3HloiIiJDYIJORERkJJYtW4aXXnoJU6dOBQAsX74ce/bswcqVKxEWFlah/O7duxEZGYmYmBg4ODgAAHx9fXXKaJdrbdq0CVZWVjoJ+ieffAIvLy+d5PvB7QCAUqmEq6txzz9LRERkSGziTkREZAQKCwtx6tQphIaG6iwPDQ3F0aNH9a6zY8cOhISEYOnSpfDw8EBAQADmzJmDgoKCSvcTHh6OcePGwdrausJ2Ro8eDWdnZ3Ts2BH/93//V2HdgwcPwtnZGQEBAXj55ZeRkpLyiEdLRERkmliDTkREZARSU1NRUlICFxcXneUuLi5ITk7Wu05MTAyOHDkClUqFbdu2ITU1FdOnT0d6erpOP3St48eP48KFCwgPD6+wnZUrV2L27Nl49913cfz4cbz++utQKpWYMGECAGDIkCEYPXo0fHx8EBsbi/nz56N///44deoUlEplhX1pNBpoNBrpeXZ2do3PCRERUWNjFAm6EAIAb95ERNRwaO9J2nvU4yKTyXSeCyEqLNMqLS2FTCbDhg0boFarAZQ1kx81ahS+/PJLWFpa6pQPDw9HUFAQOnfuXGE7ISEh+OijjwAAHTt2xMWLF7Fy5UopQR87dqxUPigoCCEhIfDx8cHOnTvxzDPPVIgtLCwMCxcurLCc93oiImooDHGvN4oEPScnBwDg5eVl4EiIiIh05eTkSMlvfXJycoJCoahQW56SklKhVl3Lzc0NHh4eOvEFBgZCCIGEhAS0aNFCWp6fn49NmzZh0aJFerfTunVrnWWBgYHYsmVLpfG6ubnBx8cH165d0/v63LlzMXv2bOl5YmIiWrduzXs9ERE1OI/rXg8YSYLu7u6O+Ph42NraVlpLUF3Z2dnw8vJCfHw87Ozs6ijCx4fxGxbjNyzGb1jGHj9Qt8cghEBOTg7c3d3rKLqqWVhYIDg4GBEREXj66ael5RERERgxYoTedXr06IEffvgBubm5sLGxAQBER0dDLpfD09NTp+z3338PjUaDF154Qe92rl69qrMsOjoaPj4+lcablpaG+Ph4uLm56X1dqVTqNH23sbHhvb6B4PmrHZ6/2uM5rB2ev9q5//zZ2to+1ns9YCQJur4PErVlZ2dn1G9Yxm9YjN+wGL9hGXv8QN0dw+P6Nl1r9uzZePHFFxESEoJu3brh66+/RlxcHKZNmwagrFY6MTER69atAwCMHz8eixcvxuTJk7Fw4UKkpqbizTffxJQpU/Q2bx85ciQcHR0r7HfWrFno3r07PvroI4wZMwbHjx/H119/ja+//hoAkJubiwULFuDZZ5+Fm5sbbty4gXfffRdOTk46XyZUhff6hofnr3Z4/mqP57B2eP5qR3v+Hve93igSdCIiIirr552WloZFixYhKSkJQUFB2LVrl1STnZSUhLi4OKm8jY0NIiIiMGPGDISEhMDR0RFjxozBhx9+qLPd6OhoHDlyBHv37tW73yeeeALbtm3D3LlzsWjRIvj5+WH58uV4/vnnAQAKhQLnz5/HunXrkJmZCTc3N/Tr1w+bN2+Gra1tPZ0NIiKixocJOhERkRGZPn06pk+frve1tWvXVljWqlUrREREVLnNgICAhw6A89RTT+Gpp57S+5qlpSX27NlT5fpERET0cCY3D7pSqcQHH3ygd8oXY8D4DYvxGxbjNyxjjx9oHMdAD8ffc+3w/NUOz1/t8RzWDs9f7Rj6/MnE454fhoiIiIiIiIgqMLkadCIiIiIiIqKGiAk6ERERERERUQPABJ2IiIiIiIioAWCCTkRERERERNQAmFSCvmLFCvj5+UGlUiE4OBiHDx82dEgAgLCwMDzxxBOwtbWFs7MzRo4ciatXr+qUmTRpEmQymc6ja9f/b+/+Y9qq2jiAfzvfFrHBOkRWOgIS4rJkRSIsMhYzNjQMAtuSJTqUTIxKMiMbxplsxhhm4h8kJv7nrxhGZjTOPzYXExa1RJgaOmcoU364BV0FnZQqgYqZ0G583z/el/ty17Ltne29p+X5JCTk3NPb5zw95z4cKL0bdH3m5uawd+9eZGVlwW63Y/v27fj1118THv+hQ4eiYnM6ndpxkjh06BBcLhfS09OxefNmDA0NKRE7ANx9991R8VssFjz77LMA1Mv9l19+iW3btsHlcsFiseDEiRO64/HK99TUFHbv3g2HwwGHw4Hdu3djeno6ofFHIhEcOHAARUVFsNvtcLlcePzxx/Hbb7/pzrF58+ao16S+vt70+IH4zRez4o+1FiwWC1577TWtj5n5v5HrpeprQCSeqvXeSEbV5lRZJyrV1rGxMWzbtg12ux1ZWVnYt28fwuFwIoYdNyrVxmTLn2p1LRXzl3Tzj8vE0aNHabVa+e6773J4eJgtLS202+0cHR01OzRu3bqVHR0dHBwc5NmzZ1lbW8u8vDz+9ddfWp/GxkZWV1dzfHxc+5qcnNSdZ8+ePVy9ejU9Hg99Ph+3bNnC4uJiXr58OaHxt7a2ct26dbrYgsGgdrytrY0ZGRk8duwYBwYGuGvXLubk5PDPP/80PXaSDAaDutg9Hg8BsLu7m6R6uT958iRfeuklHjt2jAD48ccf647HK9/V1dV0u93s7e1lb28v3W436+rqEhr/9PQ0H3roIX700Uc8d+4cvV4vy8rKWFpaqjtHRUUFm5qadK/J9PS0ro8Z8ZPxmy9mxb847vHxcR4+fJgWi4U//fST1sfM/N/I9VL1NSASS+V6bySjanOqrBNVauvly5fpdru5ZcsW+nw+ejweulwuNjc3JzwH/4QqtTEZ86dSXUvV/CXb/Fs2G/T777+fe/bs0bWtXbuWBw8eNCmipQWDQQLgqVOntLbGxkbu2LFjycdMT0/TarXy6NGjWtvFixe5YsUKfvrpp4kMl62trSwuLo55bH5+nk6nk21tbVrb7OwsHQ4H3377bdNjj6WlpYWFhYWcn58nqXbury6C8cr38PAwAfD06dNaH6/XSwA8d+5cwuKP5cyZMwSg++G6oqKCLS0tSz7GzPjjMV9Uyv+OHTtYWVmpa1Ml/2T09TLZ1oCIv2Sq94lkRG1O1XViZm09efIkV6xYwYsXL2p9PvzwQ6alpTEUCiVkvPFmZm1MhfyZWddSMX9k8s2/ZfEW93A4jL6+PlRVVenaq6qq0Nvba1JUSwuFQgCAzMxMXXtPTw+ys7OxZs0aNDU1IRgMasf6+voQiUR0Y3S5XHC73YaMcWRkBC6XCwUFBaivr8eFCxcAAH6/H4FAQBdXWloaKioqtLjMjn2xcDiM999/H08++SQsFovWrnLuF4tXvr1eLxwOB8rKyrQ+GzZsgMPhMHxMoVAIFosFd9xxh679gw8+QFZWFtatW4cXXngBMzMz2jGz4/+n88Xs+BdMTEygs7MTTz31VNQxVfJ/9fUyFdeAuHHJVu8TLdG1ebmsEyPz5fV64Xa74XK5tD5bt27F3Nwc+vr6EjrORDOiNqZC/sysa6mYvwXJNP/+9X+MN2n98ccfuHLlClatWqVrX7VqFQKBgElRxUYSzz//PB544AG43W6tvaamBg8//DDy8/Ph9/vx8ssvo7KyEn19fUhLS0MgEIDNZsPKlSt15zNijGVlZXjvvfewZs0aTExM4NVXX8XGjRsxNDSkPXes3I+OjgKAqbFf7cSJE5iensYTTzyhtamc+6vFK9+BQADZ2dlR58/OzjZ0TLOzszh48CAee+wx3H777Vp7Q0MDCgoK4HQ6MTg4iBdffBHfffcdPB6P6fHHY76okv8jR44gIyMDO3fu1LWrkv9Y18tUWwPi/5NM9T7RjKjNy2WdGJmvQCAQ9TwrV66EzWZL6pwaVRuTPX9m17VUzB+QfPNvWWzQFyz+iyjwnxfx6jazNTc34/vvv8fXX3+ta9+1a5f2vdvtxvr165Gfn4/Ozs6oH54XM2KMNTU12vdFRUUoLy9HYWEhjhw5on0Aw83k3ozXp729HTU1NbrffKmc+6XEI9+x+hs5pkgkgvr6eszPz+PNN9/UHWtqatK+d7vduOeee7B+/Xr4fD6UlJQAMC/+eM0Xs/MPAIcPH0ZDQwNuvfVWXbsq+V/qehnr+ZNxDYiblwz1PtGMqs3LaZ0Yla9UzKmRtTGZ86dCXUvF/CXb/FsWb3HPysrCLbfcEvWbi2AwGPVbDjPt3bsXn3zyCbq7u5Gbm3vNvjk5OcjPz8fIyAgAwOl0IhwOY2pqStfPjDHa7XYUFRVhZGRE+8TYa+VeldhHR0fR1dWFp59++pr9VM59vPLtdDoxMTERdf7ff//dkDFFIhE88sgj8Pv98Hg8ur+ex1JSUgKr1ap7TcyMf7GbmS8qxP/VV1/h/Pnz110PgDn5X+p6mSprQNycZKn3ZkhEbV4u68TIfDmdzqjnmZqaQiQSSamcJqo2JnP+VKhrqZi/WFSff8tig26z2VBaWqq9/XKBx+PBxo0bTYrqf0iiubkZx48fxxdffIGCgoLrPmZychK//PILcnJyAAClpaWwWq26MY6Pj2NwcNDwMc7NzeGHH35ATk6O9jbYxXGFw2GcOnVKi0uV2Ds6OpCdnY3a2tpr9lM59/HKd3l5OUKhEM6cOaP1+eabbxAKhRI+poXN+cjICLq6unDnnXde9zFDQ0OIRCLaa2Jm/Fe7mfmiQvzt7e0oLS1FcXHxdfsamf/rXS9TYQ2Im6d6vTdTImrzclknRuarvLwcg4ODGB8f1/p8/vnnSEtLQ2lpaULHaaRE1cZkzJ9KdS0V8xeL8vPvhj9OLskt3Halvb2dw8PDfO6552i32/nzzz+bHRqfeeYZOhwO9vT06D7+/9KlSyTJmZkZ7t+/n729vfT7/ezu7mZ5eTlXr14ddXuF3NxcdnV10efzsbKy0pBble3fv589PT28cOECT58+zbq6OmZkZGi5bWtro8Ph4PHjxzkwMMBHH3005q0hzIh9wZUrV5iXl8cDBw7o2lXM/czMDPv7+9nf308AfP3119nf3699ynm88l1dXc17772XXq+XXq+XRUVFcbl1zrXij0Qi3L59O3Nzc3n27FndepibmyNJ/vjjj3zllVf47bff0u/3s7Ozk2vXruV9991nevzxnC9mxL8gFArxtttu41tvvRX1eLPzf73rJan+GhCJpXK9N5JRtTlV1okqtXXhNk0PPvggfT4fu7q6mJubq/Rtrkh1amMy5k+lupaK+UvG+bdsNugk+cYbbzA/P582m40lJSW6j983E4CYXx0dHSTJS5cusaqqinfddRetVivz8vLY2NjIsbEx3Xn+/vtvNjc3MzMzk+np6ayrq4vqkwgL92K0Wq10uVzcuXMnh4aGtOPz8/NsbW2l0+lkWloaN23axIGBASViX/DZZ58RAM+fP69rVzH33d3dMedLY2Mjyfjle3Jykg0NDczIyGBGRgYbGho4NTWV0Pj9fv+S62HhvvRjY2PctGkTMzMzabPZWFhYyH379kXdz9KM+OM5X8yIf8E777zD9PT0qHubk+bn/3rXS1L9NSAST9V6bySjanOqrBOVauvo6Chra2uZnp7OzMxMNjc3c3Z2NpHD/8dUqo3Jlj/V6lqq5S8Z55/lvwMTQgghhBBCCCGEiZbF/6ALIYQQQgghhBCqkw26EEIIIYQQQgihANmgCyGEEEIIIYQQCpANuhBCCCGEEEIIoQDZoAshhBBCCCGEEAqQDboQQgghhBBCCKEA2aALIYQQQgghhBAKkA26EEIIIYQQQgihANmgCyGEEEIIIYQQCpANuhBCCCGEEEIIoQDZoAshhBBCCCGEEAqQDboQQgghhBBCCKGAfwM1sspNVpOciQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(history[\"train_loss\"])\n",
    "plt.title(\"First part of the gradient descent.\")\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(history2[\"train_loss\"] + history3[\"train_loss\"])\n",
    "plt.title(\"Second and third part of the gradient descent.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'income_w_0:0' shape=(1, 3) dtype=float32, numpy=array([[-0.08402919, -0.02359896, -0.03233599]], dtype=float32)>,\n",
       " <tf.Variable 'cost_w_1:0' shape=(1, 1) dtype=float32, numpy=array([[-0.05140902]], dtype=float32)>,\n",
       " <tf.Variable 'freq_w_2:0' shape=(1, 1) dtype=float32, numpy=array([[0.09645309]], dtype=float32)>,\n",
       " <tf.Variable 'ovt_w_3:0' shape=(1, 1) dtype=float32, numpy=array([[-0.04099115]], dtype=float32)>,\n",
       " <tf.Variable 'ivt_w_4:0' shape=(1, 4) dtype=float32, numpy=\n",
       " array([[ 0.05871323, -0.00726116, -0.00368665, -0.00105638]],\n",
       "       dtype=float32)>,\n",
       " <tf.Variable 'intercept_w_5:0' shape=(1, 3) dtype=float32, numpy=array([[-1.687437  , -0.39638832,  1.1344599 ]], dtype=float32)>]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmnl.trainable_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "keep_output": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=0.67665595>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmnl.evaluate(canada_dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "[1] ModeCanada dataset in *Application and interpretation of nested logit models of intercity mode choice*, Christophier, V. F.; Koppelman, S. (1993)\\\n",
    "[2] Conditional MultinomialLogit, Train, K.; McFadden, D.; Ben-Akiva, M. (1987)\\\n",
    "[3] Siwssmetro dataset in *The acceptance of modal innovation: The case of Swissmetro*, Bierlaire, M.; Axhausen, K.; Abay, G (2001)\\"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}