{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "from definitions import RESULTS_FOLDER, FIGURE_FOLDER\n",
    "\n",
    "figure_location = os.path.join(FIGURE_FOLDER,'twostep')\n",
    "if not os.path.exists(figure_location):\n",
    "    os.makedirs(figure_location)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_lesion_hpc = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_lesion_hpc.csv'))\n",
    "data_lesion_dls = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_lesion_dls.csv'))\n",
    "data_control = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep', 'results_control.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>Action1</th>\n",
       "      <th>Action2</th>\n",
       "      <th>Agent</th>\n",
       "      <th>DLS reliability</th>\n",
       "      <th>HPC reliability</th>\n",
       "      <th>P(SR)</th>\n",
       "      <th>Reward</th>\n",
       "      <th>StartState</th>\n",
       "      <th>State2</th>\n",
       "      <th>Terminus</th>\n",
       "      <th>Trial</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.057327</td>\n",
       "      <td>0.087327</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.115198</td>\n",
       "      <td>0.167028</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.175993</td>\n",
       "      <td>0.239769</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.224014</td>\n",
       "      <td>0.306158</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.252047</td>\n",
       "      <td>0.366749</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  Action1  Action2  Agent  DLS reliability  HPC reliability  \\\n",
       "0           0      0.0      0.0    0.0         0.057327         0.087327   \n",
       "1           1      0.0      0.0    0.0         0.115198         0.167028   \n",
       "2           2      0.0      0.0    0.0         0.175993         0.239769   \n",
       "3           3      0.0      0.0    0.0         0.224014         0.306158   \n",
       "4           4      0.0      0.0    0.0         0.252047         0.366749   \n",
       "\n",
       "   P(SR)  Reward  StartState  State2  Terminus  Trial  \n",
       "0    1.0     1.0         0.0     3.0       5.0    0.0  \n",
       "1    1.0     1.0         0.0     3.0       5.0    1.0  \n",
       "2    1.0     1.0         0.0     3.0       5.0    2.0  \n",
       "3    1.0     0.0         0.0     3.0       5.0    3.0  \n",
       "4    1.0     1.0         0.0     4.0       7.0    4.0  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_lesion_dls.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_common_or_rare(action, out):\n",
    "    left_outcomes = (5, 6)\n",
    "    right_outcomes = (7, 8)\n",
    "    if action == 0 and out in left_outcomes:\n",
    "        return 'common'\n",
    "    elif action == 0 and out in right_outcomes:\n",
    "        return 'rare'\n",
    "    elif action == 1 and out in left_outcomes:\n",
    "        return 'rare'\n",
    "    elif action == 1 and out in right_outcomes:\n",
    "        return 'common'\n",
    "    else:\n",
    "        raise ValueError('The combination of action and outcome does not make sense')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_relevant_columns(dataframe):\n",
    "    dataframe['PreviousAction'] = dataframe.groupby(['Agent'])['Action1'].shift(1)\n",
    "    dataframe['PreviousStart'] = dataframe.groupby(['Agent'])['StartState'].shift(1)\n",
    "    dataframe['PreviousReward'] = dataframe.groupby(['Agent'])['Reward'].shift(1)\n",
    "    dataframe['Stay'] = (dataframe.PreviousAction == dataframe.Action1)\n",
    "    dataframe['Transition'] = np.vectorize(is_common_or_rare)(dataframe['Action1'], dataframe['Terminus'])\n",
    "    dataframe['PreviousTransition'] = dataframe.groupby(['Agent'])['Transition'].shift(1)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "add_relevant_columns(data_control)\n",
    "add_relevant_columns(data_lesion_dls)\n",
    "add_relevant_columns(data_lesion_hpc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>Action1</th>\n",
       "      <th>Action2</th>\n",
       "      <th>Agent</th>\n",
       "      <th>DLS reliability</th>\n",
       "      <th>HPC reliability</th>\n",
       "      <th>P(SR)</th>\n",
       "      <th>Reward</th>\n",
       "      <th>StartState</th>\n",
       "      <th>State2</th>\n",
       "      <th>Terminus</th>\n",
       "      <th>Trial</th>\n",
       "      <th>PreviousAction</th>\n",
       "      <th>PreviousStart</th>\n",
       "      <th>PreviousReward</th>\n",
       "      <th>Stay</th>\n",
       "      <th>Transition</th>\n",
       "      <th>PreviousTransition</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.057327</td>\n",
       "      <td>0.087327</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>common</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.115198</td>\n",
       "      <td>0.167028</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>True</td>\n",
       "      <td>common</td>\n",
       "      <td>common</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.175993</td>\n",
       "      <td>0.239769</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>True</td>\n",
       "      <td>common</td>\n",
       "      <td>common</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.224014</td>\n",
       "      <td>0.306158</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>True</td>\n",
       "      <td>common</td>\n",
       "      <td>common</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.252047</td>\n",
       "      <td>0.366749</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>True</td>\n",
       "      <td>rare</td>\n",
       "      <td>common</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  Action1  Action2  Agent  DLS reliability  HPC reliability  \\\n",
       "0           0      0.0      0.0    0.0         0.057327         0.087327   \n",
       "1           1      0.0      0.0    0.0         0.115198         0.167028   \n",
       "2           2      0.0      0.0    0.0         0.175993         0.239769   \n",
       "3           3      0.0      0.0    0.0         0.224014         0.306158   \n",
       "4           4      0.0      0.0    0.0         0.252047         0.366749   \n",
       "\n",
       "   P(SR)  Reward  StartState  State2  Terminus  Trial  PreviousAction  \\\n",
       "0    1.0     1.0         0.0     3.0       5.0    0.0             NaN   \n",
       "1    1.0     1.0         0.0     3.0       5.0    1.0             0.0   \n",
       "2    1.0     1.0         0.0     3.0       5.0    2.0             0.0   \n",
       "3    1.0     0.0         0.0     3.0       5.0    3.0             0.0   \n",
       "4    1.0     1.0         0.0     4.0       7.0    4.0             0.0   \n",
       "\n",
       "   PreviousStart  PreviousReward   Stay Transition PreviousTransition  \n",
       "0            NaN             NaN  False     common                NaN  \n",
       "1            0.0             1.0   True     common             common  \n",
       "2            0.0             1.0   True     common             common  \n",
       "3            0.0             1.0   True     common             common  \n",
       "4            0.0             0.0   True       rare             common  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_lesion_dls.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_mean_stay_prob(data):\n",
    "    means = data[data['Trial']>0].groupby(['PreviousTransition', 'PreviousReward'])['Stay'].mean()\n",
    "    sems = data.groupby(['PreviousTransition', 'PreviousReward'])['Stay'].sem()\n",
    "    return means, sems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_lesion_hpc, sem_lesion_hpc = compute_mean_stay_prob(data_lesion_hpc)\n",
    "mean_lesion_dls, sem_lesion_dls = compute_mean_stay_prob(data_lesion_dls)\n",
    "mean_full, sem_full = compute_mean_stay_prob(data_control)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PreviousTransition  PreviousReward\n",
       "common              0.0               0.558473\n",
       "                    1.0               0.836690\n",
       "rare                0.0               0.563003\n",
       "                    1.0               0.823360\n",
       "Name: Stay, dtype: float64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_lesion_hpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_daw_style(ax, data, yerr, title=''):\n",
    "    lightgray = '#d1d1d1'\n",
    "    darkgray = '#929292'\n",
    "\n",
    "    bar_width= 0.2\n",
    "\n",
    "    bars1 = data[:2][::-1]\n",
    "    bars2 = data[2:][::-1]\n",
    "    errs1 = yerr[:2][::-1]\n",
    "    errs2 = yerr[2:][::-1]\n",
    "\n",
    "    # The x position of bars\n",
    "    r1 = np.array([0.125, 0.625]) \n",
    "    r2 = [x + bar_width + .05 for x in r1]\n",
    "    \n",
    "    plt.sca(ax)\n",
    "    \n",
    "    plt.bar(r1, bars1, width=bar_width, color='blue',  capsize=4)\n",
    "    plt.bar(r2, bars2, width=bar_width, color='red',  capsize=4)\n",
    "    plt.xticks([r+ bar_width/2 +.025 for r in r1], ['Rewarded', 'Unrewarded'], fontsize=12)\n",
    "    plt.yticks(fontsize=12)\n",
    "    plt.title(title, fontsize=16)\n",
    "    plt.ylim([0.4, 1])\n",
    "    plt.xlim([0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAC8CAYAAADIKoJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+cVVW9//HXmx8GDiCigCEo3ZTqkr+Kr5b5A70aivk7vCBg3jJT8iZp5c8UEfNqZfdWpPkrEtMyU9SQ1CwQSlJM8MpNqRRU8AeIIIMigp/vH3vPuDmcM7MH5pw5M/N+Ph77MXuvvfZZ6+yz15x11t5rLUUEZmZmZlZ9OrR0BszMzMysOFfUzMzMzKqUK2pmZmZmVcoVNTMzM7Mq5YqamZmZWZVyRc3MzMysSrmiVkUkHSfpEUmvSXpb0hJJ0yQdkYkzVNIESbk/O0mnSgpJA5uYn55pWp9oynHWfmWutd2K7OuU7ptQEHdghbNpVrUy5aLYctgWvF59mUu3J0iq6nG5JM2UNHMLjmuT/1NcUasSkr4G3A38HfgScBQwKd19aCbqUOBSmvbZTQc+DbzcxGz1TNNyRc3KYUuvS7P2YARJ+cguj7VojqxFdGrpDFi9bwDTIuJLmbA/ADc0pfUsS1JnYENELAeWN0MezZqNr0uzBs2PiH+0dCas5blFrXr0Al4ptiMi3oOkyZqkhQvg3brm8HTfwHR7nKSrJS0D3gF6FmsOljRS0h8kLZdUK+lJSV/I7B8IPJ9u3pBpej813b9Y0pTCvJZqZpf0UUkPSFor6QVJ/5HuHyvpmTQPf5T04SafOWuVSlyXiyXdKunLkv4haZ2kv0o6pODYKZJekrS/pMfTeIsl/WeRdPaV9Pv0Glsr6WFJ+xaJd7CkhyStTuMtkPSlzP4Gy0wmXkiaJOnc9PGFtZKmS+qTLnekabwo6bwS5+Sg9LGHWkmvS5osqWsm3tA03tAc5/TkNK+1abr/K+krDX44VtVK3eJrztuambI4VtKzSh7HmS1pd0k1kn6aXpuvSvq+pE4Fx39E0t2SVqXHzlXmMZ5MvJHpd8A7khZKOr5EfnaUdK2kpWncZySd3hzvtdq5Ra16PAZ8QdJzwD0RsahInBuB/iS3Rg8ANhaJcxHwOHA60BFYVyK9fwHuBP4LeA84CLhRUteIuI7kdtQJwF3AlcC96XH/bPpbA+DXwA3A94BxwM2Sdie5lXs+0Bn4H+A2YL8tTMOqR8fCf9wk12MeBwOfJLmW3wHOA2ZI2isins3E6wH8CrgK+AcwEvihpDURMQVA0p7ALOD/gFOBILneZkn6VEQsSOMdC/wG+BPwFWAFMBjYNZNeY2UmayzwNMm13hf4b+AWoDswA7ie5NbWf0n634i4v+D4W4E7gJ8A+wKXADXpe8hN0gHpa/0Q+CbJj/OPkjzWYNWtsAxFRBT7n19OBwEfJimD25Bcx78BnuP9MncQcDHJd8NPACT1A+YAa4CzgNXAV4Hpkj4XETPSeIeR/M+fDpwL9Cb5HugM1Jd1ST1IymZXYAJJI8Iw4FpJH4iIH5XrBFSFiPBSBQswCHiK5IskSL4obgc+WxBvQrq/U0H4wDT8r4AK9p2a7htYIu0OJJX2G4AFRV7ztCLHLAamFAkPYEKR/J6SCdse2AC8DvTIhH8tjbtrS38eXrZsyVxrDS0TCuIOzBy/GFgP7JIJ6w6sBKZmwqakx44sSP8hYEldGSCpWK0Cembi9Ehf7650W2m684AOOd9n0TKT7gtgUbaMAtek4RdnwjoBrwE/K3L+rit4zYtIfpgNSreHpvGGljj/A9PtbwArW/q68JJ/aaAMzSn1OWfCJ5BU6AqvxwkNxSmRj8VpOdkuE1b3P/rGgrh/Bf6Y2f4eyf/43TJhHUkqX3/NhP2J5EdUh0zYfmkaMzNh3yZpdNi9IN0bSL4rOzV0Xlr74lufVSKSFrR9SFoTrgDmA8cDD0i6uAkvNS3SK7YhafP17ZKWAu+my2nAR5qc+Xxm1K1ExBskX1BzI+LNTJxn0r8DypQHq5zjgf9XsHwq57FzI+KFuo2IWMP7HQ+yNpL8us/6JbALsHO6fRDw24hYlXm9N0laiA9Ogz5C0nJ2Y6SPGRTTxDLzUERsyGzXXdsPZPKxgaRVotj1fkeR99WBpHWtKR4Htk9vYX1OklvSWo/CMvSlhqOXxaMRsTqzvdl1nAnPXscHkZTj+mfsImkNvB3YW1IPSR1J3ted2XIXEX8hqSRmHQH8BXheSe/xTmlr4wPADsC/bukbbA1867OKpBfyI+lS13z8O+BSSZPTCk5jGu1BJ6kbScvDWyS3gf5J0opxJvDFLct9owrzvr5EGECXMuXBKufpKHgQusit0FJeLRG2c0HYGxHxboljdwZeInn2s1iZeIWkZReSf/Sk8YvagjJT6touFl7sei88B9n3lVtEzJI0AvhPkl7lSJoFnBMRTzXltaziNitDLWBLr+NewJNFXu8Vkhbs7UluY3amdHnP6gPsRvLjqJgdSoS3Ca6oVbGIWCbpRpJ79ruTr2t2ngdJP03SgnBgRMypC2zCFykkzdDbZAMk9WrC8Wal9C0RtrQgbHtJnQsqa3XH1sVdCexU5PV2SvdBcusEGq4ENUeZaYq+wMKCbXj/fdU9e7pJGaTIF1ZE3AncmVY2h5I80/c7Sf0bakG0qpb7828hDZW7SPe/RVLxKlXel2S2Xye5C3N2ifSeLRHeJvjWZ5WQVOp230fTv3U9Qt9J/3YtEjevbdO/9V9wkrYHji2I11BaS4CPF4R9bivyZFbnU9nyIKk7ybiCjxbE6wicWBA2EniB9ys0s4Cj0tfIvt7R6T5InidbDJwmSSXylLfMNJeTCrZHknRgqPuxVvclVlgGh5d6wYiojYjfAj8FPkj1fKlb0232+ac/Gj7bMtnZzCyScjywLiC91fnvwJMRsSa9g/Q48HllhqCStB/J89FZvyP5LnwhIuYVWdaU9+20LLeoVY+nJf2R5PbE8yQPPA8HzgDuyDyz83/p33MlzQA2RsS8Jqb1Z+BNYLKkS0l6k11M0rKwXSbeqyS/ZEZKegpYCzwfEa+TPDNzs6QfAL8F9qKJPdLMSngVeFDJMC91vT5rgMsL4q0Brpa0I8lA0aOAw4BTM89pXk7yA+JhSVeR/Jo/j6TiNRGSp6oljSfp4fwHSdeRjO/2MaBPRFxK/jLTXIZL+i7wIMlzaZcCt6TPshIRL6e3MC+QtIKktWEMSQ+9epImkrRO/BFYRtJr/GskY3R5DLvW63GS2+/fTSs575D0MP5Ai+bqfT8g+T54KC0vb5LkbxDJj646l5Jc49Mk/ZSk1+dlbD5U1Q9IKnmz0++cZ0nK4EdJWrnL9YOpKrhFrXqcR/J5TCS5cH9FcrvlfJKu/nV+S9IFehxJC8PjTU0o/Qd9PEmLxJ0kw2/cSNKNPxvvPZKHpbcHfp+mdXS6++ckhewE4D6SrtJFx78xa6JZwPeB75CUgy7AkbH5kDVvkrQ0fQG4BzgEODsifl4XIX0Oa2ga9+fAVKAWODjSoTnSePcAh6ebN5F0Njid9KHmvGWmGY0h+VK7m2TYghtIynxhnLkkQ29MIWlJnFQQ5y8krRM/IHnG7irSVsbyZNsqIe2IcizwIslnP5nk853Scrl6X0QsIxlCaiFwLUmZ6QUcFRG/y8T7PTCapEPOXSRDyIyn4FZm2qFhf+B+ku/KB4CbSc7BH8v8dlqcovEOgmZmFSFpMckwBGMaiTcFOCwi+lciX5WiZEDpn5EMQ9DSD5KbWRVwi5qZmZlZlapYRU3SWZLmpVM/TGkk7tclvaJkupObJVXLfXczMzOziqnYrU9JJ5D0WhoGdI2IU0vEG0Yy1cqhJA+/3k0ycN75FcmomZmZWZWoWItaRNwVEdNIehE25AvATRGxMB3g9XLcm9DMzMzaoWp8Rm0wsCCzvQDoK8lj/piZmVm7Uo3jqHUDsnOL1a13p6A1TtLpJF3oqamp+eRHP/pRzKrBE088sSIielc6XZcJq1YuE2abylsmqrGiVksy2GuduvXNRh6OiOuB6wGGDBkS8+Y1ddxXs/KQtKTxWM3PZcKqlcuE2abylolqvPW5kGSU+zp7Aa+mo+GbmZmZtRuVHJ6jk6QuJCN7d5TUpcSExrcAX5L0r+lcehdTJaMtm5mZmVVSJVvULgbeJpkSaUy6frGkXSTVStoFIJ1e4mqSaSGWpMulFcynmZmZWVWo2DNqETEBmFBid7eCuNcA15Q5S2ZmZmZVrRqfUTMzMzMzXFEzMzMzq1quqJmZmZlVKVfUzMzMzKqUK2pmZmZmVcoVNTMzM7Mq5YqabUJqeGmq2267jSFDhtCtWzc++MEPcuSRRzJnzpzmz7iZmVkb5Iqalc0111zD+PHjufDCC3n11Vd54YUXGDduHPfcc09LZ83MzKxVcEXNymL16tVccsklTJ48mRNOOIGamho6d+7M0UcfzXe/+13eeecdxo8fT79+/ejXrx/jx4/nnXfeAWDmzJn079+fq6++mj59+vDBD36QadOmcf/99zNo0CB69erFd77znfq0JkyYwIgRIxgzZgzdu3dnjz32YNGiRVx55ZX06dOHAQMG8OCDD9bHX7ZsGccccwy9evVit91244YbbtjktU466SROOeUUunfvzuDBg/EkzmZm1lJcUbOyePTRR1m3bh3HH3980f1XXHEFc+fOZf78+SxYsIDHHnuMSZMm1e9/5ZVXWLduHUuXLmXixIl8+ctf5tZbb+WJJ55g9uzZTJw4keeee64+/n333cfYsWN544032GeffRg2bBjvvfceS5cu5ZJLLuErX/lKfdxRo0bRv39/li1bxp133smFF17Iww8/XL//3nvvZeTIkaxatYpjjjmGs846qwxnyMzMrHGuqFlZvP766+y444506lR8lrJf/OIXXHLJJfTp04fevXtz6aWXMnXq1Pr9nTt35qKLLqJz586MHDmSFStWcPbZZ9e3cg0ePJinnnqqPv6BBx7IsGHD6NSpEyNGjGD58uWcf/759ccvXryYVatW8eKLLzJnzhyuuuoqunTpwt57781pp522SdoHHHAAw4cPp2PHjowdO5YFCxaU70SZmZk1wBU1K4sddtiBFStWsGHDhqL7ly1bxq677lq/veuuu7Js2bJNju/YsSMAXbt2BaBv3771+7t27UptbW39duG+HXfccbPja2trWbZsGb169aJ79+6bpL106dL67Z122ql+fdttt2XdunUl34eZmVk5uaJmZfHpT3+aLl26MG3atKL7+/Xrx5IlS+q3X3jhBfr161f2fPXr14+VK1eyZs2aTdLeeeedy562mZlZU1Wsoiapl6S7Ja2VtETSySXi9ZT0c0mvpcuESuXRms92223HxIkT+epXv8q0adN46623ePfdd5kxYwbf+ta3GDVqFJMmTWL58uWsWLGCiRMnMmbMmLLna8CAAey///5ccMEFrFu3jqeeeoqbbrqJ0aNHlz1tMzOzpir+AFF5TAbWA32BvYHpkhZExMKCeD8AtgUGAn2AhyUtiYifVTCv1gzOOecc+vbty6RJkxg9ejTdu3fnk5/8JBdddBGf+MQnePPNN9lzzz0BGDFiBBdffHFF8nX77bdzxhln0K9fP7bffnsuu+wyDj/88IqkbWZm1hSKiPInItUAbwAfj4hFadhUYGlEnF8QdwVwZEQ8nm5fmG4f2FAaQ4YMCQ+jYNVC0hMRMaQl8+AyYdXEZcJsU3nLRK5bn5KelDReUt/GYxc1CNhYV0lLLQAGl0qyYP3jW5iumZmZWauV9xm1ScBBwHOSZkg6WVLXJqTTDVhdELYa6F4k7u+A8yV1l7Qb8EWSW6GbkXS6pHmS5i1fvrwJ2TFrm1wmzDblMmGtXa6KWkT8JiJOAAYA9wDjgJcl3Szp0BwvUQv0KAjrAawpEvdrwNvA39O0bgdeKpGv6yNiSEQM6d27d563YtamuUyYbcplwlq7JvX6jIiVwC3AdcALwInA9ZIWSTqsgUMXAZ0k7Z4J2wso7EhARKyMiNERsVNEDE7z+FhT8mlmZmbWFuR9Rq2DpGGSbgWWAaOB/wJ2iojdgAuAW0sdHxFrgbuAiZJqJH0GOBaYWhhX0ocl7SCpo6QjgdNJbr2amZmZtSt5W9SWAd8HngL+NSKOjIjbIuJtSG6NAn9r5DXGAV2B10huZ54ZEQslHSipNhPvk8D/ktwWvRIYXWQIDzMzM7M2L+84ap+LiAb7NEfEIY3sXwkcVyR8Nklng7rtO4A7cubLzMzMrM3K26L2YLFASa81Y17MzMzMLCNvRa1zYYCkzkDH5s2OmZmZmdVp8NanpNlAAF0kPVKwuz/w53JlzFqI1PD+JsxkMXDgQF599VU6duxIt27dOOKII/jxj39Mt27dGj/YzMzMGm1RuxG4GdgA3JRZbgTOBE4oa+6s1bvvvvuora1l/vz5PPnkk1x55ZVNfo0NGzaUIWdmZmbVr8GKWkT8PCKmAPuk63XLLRHxQES8W5lsWmu30047MWzYMObPnw/A9OnT2WeffejRowcDBgxgwoQJ9XEXL16MJG666SZ22WUXDj00GVN57ty57L///vTs2ZO99tqLmTNntsA7MTMzq5yStz4ljY2IunHO9pe0f7F4EXFzWXJmbcpLL73EjBkz6itdNTU13HLLLQwePJinn36aww8/nL333pvjjnu/Y/CsWbP429/+RocOHVi6dClHHXUUU6dO5YgjjuDhhx/mxBNP5JlnnsGjjZuZWVvVUIvaqMz62BLLmPJlzdqC4447ju7duzNgwAD69OnDZZddBsDQoUPZY4896NChA3vuuSejRo1i1qxZmxw7YcIEampq6Nq1K7feeivDhw9n+PDhdOjQgcMPP5whQ4Zw//33t8TbMjMzq4iSFbWIGJ5ZP6TEkmeeT2vHpk2bxpo1a5g5cybPPPMMK1asAOAvf/kLhxxyCL1792a77bbjuuuuq99XZ8CAAfXrS5Ys4de//jU9e/asX+bMmcPLL79c0fdjZmZWSSUraum0UY0ulcxsU0gtv9j7Dj74YE499VS+8Y1vAHDyySdzzDHH8OKLL7J69WrOOOMMoqBHqTInccCAAYwdO5ZVq1bVL2vXruX888+v6PswMzOrpIYqWhuAdxtY6vab5TJ+/Hgeeugh5s+fz5o1a+jVqxddunThscce47bbbmvw2DFjxnDffffxwAMPsHHjRtatW8fMmTN56aWXKpR7MzOzymtoHLUPVSwX1i707t2bU045hcsvv5yf/OQnnHvuuZx11lkcfPDBnHTSSaxatarksQMGDOCee+7hW9/6FqNGjaJjx47su+++XHvttRV8B2bWllTDnY8mDE1p7VTJilpELKlkRtqlavwv0Yz/NRYvXrxZWLZi9fnPf77ocQMHDtzsNijAfvvtt1mHAzMzs7asoeE5ro+I09P1qSQzFGwmIk7Jk5CkXiSD5X4WWAFcEBGb3e+S9AHgf4DjSaau+hNwRkQszZOOmZmZWVvR0K3P5zPr/2iGtCYD64G+wN7AdEkLImJhQbyzgU8DewKrgRuAH+FZEMzMzKydaejW55WZ9cu2JhFJNcCJwMcjohaYI+lekrHYCrvtfQh4ICJeTY/9JXDN1qRvZmZm1hrlHl5D0qGSbpA0Pf37b01IZxCwMSIWZcIWAIOLxL0J+IykfpK2BUYDM5qQlpmZmVmbkKuiJukc4JfASmA68Dpwm6Rzc6bTjeQ2ZtZqoHuRuIuAF4ClwJvAx4CJJfJ1uqR5kuYtX748Z1bM2q6GykRLjytYDX1nrP3x94S1dnlb1M4FDo2I8yLiJxFxPnBoGp5HLdCjIKwHsKZI3GuBLsAOQA1wFyVa1CLi+ogYEhFDPN+jmcuEWSGXCWvtmjKzQGGHguco0RO0iEVAJ0m7Z8L2Ago7EtSFT4mIlRHxDklHgn0l7diEvJqZmZm1ermmkAImADdJ2l1SV0mDgOuBS/MkEhFrSVrGJkqqkfQZ4FhgapHojwOnSNpOUmdgHLAsIlYUiWtmZmbWZjU0PMcG3m8xq3u6ZFRB2MnAjTnTGgfcDLxG8ozbmRGxUNKBwIyI6JbG+wbwQ+DvwDbA0yRjqpmZmZm1KxWbQioiVgLHFQmfTdLZoG77dZKenmZmZmbtmqeQMjMzq1bV0F3aE5K2qIZa1DYh6RjgYGBH3r8VmnsKKTMzMzNrmrzjqF0K/DSNP4LkGbNhwKryZc3MzMysfcs7PMcXgcMj4uvA+vTv0cDAcmXMzMzMrL3LW1HrGRFPp+vrJXWOiMdIboWamZmZWRnkfUbtn5IGR8RCkuEyzpT0BvBG+bJmZmZm1r7lrahdTDKlE8AFwC9IhtQYV45MmZmZmVnOilpE3J9Z/wuwW9lyZGZmZmZA04bn2B04CegHLAPuiIi/lytjZmZmZu1d3uE5TgaeBPYE1gJ7AH9Nw83MzMysDPK2qE0ChkfEI3UB6RydU4HbypExMzMzs/Yu7/Ac3YFHC8LmAjXNmx0zMzMzq5O3onYN8B1JXQAkdQWuSMPNzKyaSC2/mFmzKHnrU9KLQN1MrAJ2As5Ox0/bPg17GbgyT0KSegE3AZ8FVgAXRMRmt00lzQAOzARtAzwbEXvkScfMzMysrWjoGbUxzZzWZGA90BfYG5guaUE6iG69iDgyuy1pJvCHZs6LmZmZWdUrWVGLiFnNlYikGuBE4OMRUQvMkXQvMBY4v4HjBpK0rv1Hc+XFzMzMrLXIOzxHZ0mXSXpO0rr072WStsmZziBgY0QsyoQtAAY3ctwpwOyIeD5nOmZmZmZtRt7OBFcDhwFnAHulfw8Frsp5fDdgdUHYapLepA05BZhSaqek0yXNkzRv+fLlObNi1na5TJhtymXCWru8FbURwDER8WBEPBsRDwLHk8xUkEct0KMgrAewptQBkg4g6cBwZ6k4EXF9RAyJiCG9e/fOmRWztstlwmxTLhPW2uWtqJXqa523D/YioFM6DVWdvYCFJeIDfAG4K32mzczMzKzdyVtR+zVwn6Rhkj4m6QhgGnBHnoMjYi1wFzBRUo2kzwDHksxssJl0nLYRNHDb08zMzKyty1tR+xbwe5IhNp4AfgT8ETivCWmNA7oCrwG3A2dGxEJJB0oqbDU7juQZtj824fXNzMzM2pRG5/qU1JFkTLXvRMQlW5pQRKwkqYAVhs8m6WyQDbudpDJnZmZm1m412qIWERuBayJiXQXyY2ZmZmapvLc+75N0dFlzYmZmZmabaPTWZ6oLcKekR4HsHKBExCnlyJiZmZlZe5e3ovZ0upiZmZlZheSqqEXEZeXOiJmZmZltKm+LGpIOBUYB/YBlwC8j4uFyZczMzMysvcs7Kfs5wC+BlcB04HXgNknnljFvZmZmZu1a3ha1c4FDI6L+OTVJU4GHgO+XI2NmZmZm7V3e4TkA/lGw/RyZ3p9mZmZm1rzyVtQmADdJ2l1SV0mDgOuBSyV1qFvKlkszMzOzdijvrc+fpn9HkbSiKd0ene5TGt6xWXNnZmZm1o7lrah9qKy5MDMzM7PN5B1HbUm5M2JmZmZmm6rYc2WSekm6W9JaSUskndxA3E9IekRSraRXJZ1dqXyamZmZVYvcA942g8nAeqAvsDcwXdKCiFiYjSRpR+B3wNeBO4FtgP4VzKeZmZlZVahIi5qkGuBE4NsRURsRc4B7gbFFop8DPBARv4iIdyJiTUT8rRL5NDMzM6smeWcmOEbS1rS+DQI2RsSiTNgCYHCRuJ8CVkr6s6TXJN0naZcS+Tpd0jxJ85YvX74V2TNrG1wmzDblMmEASC2/bKG8LWqXAy9L+rGk/bYgnW7A6oKw1UD3InH7A18AzgZ2AZ4Hbi/2ohFxfUQMiYghvXv33oJsmbUtLhNmm3KZsNYuV0UtIvYCDgPeBn4j6VlJF0samDOdWqBHQVgPYE2RuG8Dd0fE4xGxDrgM2F/SdjnTMjMzM2sTcj+jFhELIuKbwADgq8AI4J9p78zRjcxMsAjoJGn3TNhewMIicZ9i06mp6ta3vN3QzMzMrBVqUmcCSR8GLgGuBbqk6zcAZ5H00CwqItYCdwETJdVI+gxwLDC1SPSfAcdL2ltSZ+DbwJyIWNWUvJqZmZm1drk6CEj6KkkPzd2AO4CxETE3s/83wGuNvMw44OY03uvAmRGxUNKBwIyI6AYQEX+QdCEwHdgWmAOUHHPNzMzMrK3K25PzSOD7wD0Rsb5wZ0S8JemEhl4gIlYCxxUJn03S2SAbdi1Jq52ZmZlZu5V3CqnP5Yjz4NZnx8zMzMzq5B4bTdIxwMHAjmQe7I+IU8qQLzMzM7N2L++At5cCP03jjyB5xmwY4Af8reW19CCGWzGQoZmZWUPy9vr8InB4RHwdWJ/+PRoYWK6MWXVo6fqP60BmZtae5a2o9YyIp9P19ZI6R8RjJLdCzczMzKwM8j6j9k9JgyNiIfA0cKakN4A3ypc1MzMzs/Ytb0XtYmCHdP184DaSITW+Wo5MmZmZmVn+4Tnuz6w/RjLwrZmZmZmVUd5enytLhDc2G4GZmZmZbaG8nQk6Fwak83B2bN7smJmZmVmdBm99SpoNBNBF0iMFu/sDfy5XxszMzMzau8aeUbuRZBaC/wfclAkP4FXgD2XKl5mZmVm712BFLSJ+DiBpbkQ8szUJSepFUtn7LLACuCAibisSbwJwEfBOJnjPiHhua9I3MzMza20afEZN0iclfbyukiapt6RfSFog6TpJ3ZqQ1mRgPdAXGA1cK2lwibi/iohumcWVNDMzM2t3GutM8N/ATpntG4FBwPXAx4Gr8yQiqQY4Efh2RNRGxBzgXmBsk3NsZmZm1k40VlH7GDAbQFJP4EhgdERMBkaRzPeZxyBgY0QsyoQtAEq1qB0taaWkhZLOzJmGmZmZVVhLzwfd1ueEbqyi1onkdiXAp4BX6ipbEfEi0DNnOt2A1QVhq4HuReLeQVJB7A18GbhE0qhiLyrpdEnzJM1bvnx5zqyYtV0uE2abcpmw1q6xitpCYES6PhL4fd0OSTuzeeWrlFqgR0FYD2BNYcSI+L+IWBYRGyPiz8D/AJ8v9qIRcX1EDImIIb17985yQu60AAAHVElEQVSZFbO2y2XCbFMuE9baNTY8x3nAfZKuAzYCB2T2/Tvwp5zpLAI6Sdo9Iv6ehu1FUhFsTJAMEWJmZmbWrjTYopY+9L8LcDjwLxHxbGb3dODreRKJiLXAXcBESTWSPgMcC0wtjCvpWEnbK7Ev8DXgnlzvxszMzKwNaXQKqYhYExFPRMSagvBnI2JZE9IaB3QFXgNuB86MiIWSDpRUm4k3EvgHyW3RW4Cr6sZzMzMzM2tPGrv12WwiYiVwXJHw2SSdDeq2i3YcMDMzM2tv8k7KbmZmZmYV5oqamZmZWZVyRc3MzMysSrmiZmZmZlalXFEzMzMzq1KuqJmZmZlVqYoNz2Fm1qhqmF05oqVzYGZWzy1qZmZmZlXKFTUzMzOzKuWKmpmZmVmVckXNzMzMrEq5omZmZmZWpVxRMzMzM6tSFauoSeol6W5JayUtkXRyI/G3kfSMpJcqlUczs8ZILb+YWftRyXHUJgPrgb7A3sB0SQsiYmGJ+N8EXgO6VSh/ZmZmZlWlIi1qkmqAE4FvR0RtRMwB7gXGloj/IWAMcGUl8mdmZmZWjSp163MQsDEiFmXCFgCDS8T/EXAh8Ha5M2ZmZmZWrSp167MbsLogbDXQvTCipOOBThFxt6ShDb2opNOB09PNWknPNkNem02RZ0l2BFZUPicNaAUPvLTS87hry2TDZWKruUw0D5eJXFrpZ1l1Wul5zFUmFBWY107SPsCfImLbTNi5wNCIODoTVgPMB4ZHxN/TitqtEdG/7JmsAEnzImJIS+ejtfN5bDv8WTYPn8e2w59l82hL57FSLWqLgE6Sdo+Iv6dhewGFHQl2BwYCs5XUPLcBtpP0CvCpiFhcmeyamZmZtbyKVNQiYq2ku4CJkk4j6fV5LLB/QdSngQGZ7f2BHwOfAJZXIq9mZmZm1aKSA96OA7qSDLlxO3BmRCyUdKCkWoCI2BARr9QtwErgvXR7YwXzWi7Xt3QG2gifx7bDn2Xz8HlsO/xZNo82cx4r8oyamZmZmTWdp5AyMzMzq1KuqFUxSVMkTar0se2Jz3Hr4s+r/HyOWxd/XuXX0ue43VbUJC2W9LakWkmvpCfT01UVISkk7VYQNkHSrS2VJ2t+LhP5uUy0Dy4T+blMlE+7railjo6IbiS9UPcBLmiJTEiq5JyrLU5SxxZIs12d463gMtECXCaqmstEC3CZeF97r6gBkPYwfYCkICLpA5K+J+kFSa9Kuk5S13TfLEknpusHpL8ihqfbh0man65/WNIfJL0uaYWkX0jqWZdm+kvtPElPAWsldZK0j6S/Sloj6VdAl2w+JX1O0nxJqyT9WdKemX0NHltOkoZKeknSuZJek/SypP/I7J8i6VpJ90taCxzic1zdXCa2jstE2+MysXVcJracK2qApP7AkcA/0qCrSOYn3RvYDdgZuCTdNwsYmq4fBDwHHJzZnlX3siSTyvcDPkYyPtyEgqRHAUcBPUk+i2nAVKAX8GuSiezr8vgJ4GbgK8AOwE+Be9MLeZuGjq2QnYDtSM7Vl4DJkrbP7D8ZuIJk2rA5+BxXNZeJZuEy0Ya4TDQLl4ktERHtcgEWA7XAGiCAh9MPScBa4MOZuJ8Gnk/X/w14Kl3/HXAaMDfdngWcUCK944AnC9L/Ymb7IGAZ6ZApadifgUnp+rXA5QWv+SzJhdngsc1wrgLYrSBsAsn0XpAUlrdJ5mit2/8ayWwSAFOAWzL7fI6rcHGZcJmopnNcDYvLhMtENZzjqrwfW0HHRcTvJR0M3EYyies2wLbAE3p/AlUBdffLHwUGSepLUss/BrhM0o7AvsAjAJL6AD8EDiT5ddABeKMg/Rcz6/2ApZF+sqklmfVdgS9I+s9M2DbpcdHIsVtrI9C5IKwz8G5m+/WI2JDZfgvIPnSbfa+98TmuVi4T+bhMuEz489qUy0SZzrFvfQIRMYukNv89YAVJrX9wRPRMl+0ieZiUiHgLeAI4G3g6ItaT1JjPAf4ZESvSl72S5EPbMyJ6AGNILrJNks6svwzsrMwVCeySWX8RuCKTp54RsW1E3J7j2K31AskcrFkfomkXYPa9+hxXOZeJRrlMuEz489qUy0SZzrErau/7b+BwYE/gBuAHaU0cSTtLGpaJOws4i/fvgc8s2Iak5l4LrJK0M/DNRtJ/FNgAfC19mPEEkl8FdW4AzpC0nxI1ko6S1D3HsVvrV8DFkvpL6iDpMOBo4M4tebGIeA+f49bAZaI0lwmXCX9em3KZKNM5dkUtFRHLgVuAbwPnkTwwOlfSm8DvgY9kos8i+fAfKbENcBnJZPKrgenAXY2kvx44ATiVpFn237PHRMQ84Mskk9S/kebv1DzHNoOJJL9U5qSvfzUwOiKe3orX9Dmuci4TDXKZcJnw57Upl4kynWPP9WlmZmZWpdyiZmZmZlalXFEzMzMzq1KuqJmZmZlVKVfUzMzMzKqUK2pmZmZmVcoVNTMzM7Mq5YqamZmZWZVyRc3MzMysSrmiZmZmZlal/j+04rKFQHEKtgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x180 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize= (10,2.5), sharey=True)\n",
    "\n",
    "plot_daw_style(axes[0], list(mean_lesion_hpc), list(sem_lesion_hpc), title='Striatum')\n",
    "plot_daw_style(axes[1], list(mean_lesion_dls), list(sem_lesion_dls), title='Hippocampus')\n",
    "plot_daw_style(axes[2], list(mean_full), list(sem_full), title='Full model')\n",
    "\n",
    "\n",
    "\n",
    "leg = axes[0].legend(['Common', 'Rare'], fontsize=12, frameon=False, handlelength=0.7, title='Previous transition')\n",
    "plt.sca(axes[0])\n",
    "plt.ylabel('Stay probability', fontsize=12)\n",
    "\n",
    "plt.savefig(os.path.join(figure_location, 'StayProbability.pdf'))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Doll et al analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "figure_location = os.path.join(FIGURE_FOLDER,'twostep_deterministic')\n",
    "if not os.path.exists(figure_location):\n",
    "    os.makedirs(figure_location)\n",
    "\n",
    "# load data\n",
    "data_lesion_hpc = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_lesion_hpc.csv'))\n",
    "data_lesion_dls = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_lesion_dls.csv'))\n",
    "data_control = pd.read_csv(os.path.join(RESULTS_FOLDER, 'twostep_deterministic', 'results_control.csv'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>Action1</th>\n",
       "      <th>Action2</th>\n",
       "      <th>Agent</th>\n",
       "      <th>DLS reliability</th>\n",
       "      <th>HPC reliability</th>\n",
       "      <th>P(SR)</th>\n",
       "      <th>Reward</th>\n",
       "      <th>StartState</th>\n",
       "      <th>State2</th>\n",
       "      <th>Terminus</th>\n",
       "      <th>Trial</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.030000</td>\n",
       "      <td>0.030000</td>\n",
       "      <td>0.681250</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.059100</td>\n",
       "      <td>0.059100</td>\n",
       "      <td>0.657567</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.087327</td>\n",
       "      <td>0.087327</td>\n",
       "      <td>0.644885</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.084707</td>\n",
       "      <td>0.114707</td>\n",
       "      <td>0.632511</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.109166</td>\n",
       "      <td>0.141266</td>\n",
       "      <td>0.636400</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Unnamed: 0  Action1  Action2  Agent  DLS reliability  HPC reliability  \\\n",
       "0           0      0.0      0.0    0.0         0.030000         0.030000   \n",
       "1           1      1.0      1.0    0.0         0.059100         0.059100   \n",
       "2           2      1.0      1.0    0.0         0.087327         0.087327   \n",
       "3           3      1.0      1.0    0.0         0.084707         0.114707   \n",
       "4           4      1.0      1.0    0.0         0.109166         0.141266   \n",
       "\n",
       "      P(SR)  Reward  StartState  State2  Terminus  Trial  \n",
       "0  0.681250     0.0         0.0     2.0       4.0    0.0  \n",
       "1  0.657567     0.0         1.0     3.0       7.0    1.0  \n",
       "2  0.644885     0.0         1.0     3.0       7.0    2.0  \n",
       "3  0.632511     1.0         1.0     3.0       7.0    3.0  \n",
       "4  0.636400     0.0         1.0     3.0       7.0    4.0  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_control.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_relevant_columns(dataframe):\n",
    "    dataframe['PreviousAction'] = dataframe.groupby(['Agent'])['Action1'].shift(1)\n",
    "    dataframe['PreviousStart'] = dataframe.groupby(['Agent'])['StartState'].shift(1)\n",
    "    dataframe['PreviousReward'] = dataframe.groupby(['Agent'])['Reward'].shift(1)\n",
    "    dataframe['Stay'] = (dataframe.PreviousAction == dataframe.Action1)\n",
    "    dataframe['SameStart'] = (dataframe.StartState == dataframe.PreviousStart)\n",
    "\n",
    "\n",
    "add_relevant_columns(data_control)\n",
    "add_relevant_columns(data_lesion_dls)\n",
    "add_relevant_columns(data_lesion_hpc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_mean_stay_prob(data):\n",
    "    means = data[data['Trial']>0].groupby(['PreviousReward', 'SameStart'])['Stay'].mean()\n",
    "    sems = data.groupby(['PreviousReward', 'SameStart'])['Stay'].sem()\n",
    "    return means[::-1], sems[::-1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_lesion_hpc, sem_lesion_hpc = compute_mean_stay_prob(data_lesion_hpc)\n",
    "mean_lesion_dls, sem_lesion_dls = compute_mean_stay_prob(data_lesion_dls)\n",
    "mean_full, sem_full = compute_mean_stay_prob(data_control)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PreviousReward  SameStart\n",
       "1.0             True         0.953491\n",
       "                False        0.956948\n",
       "0.0             True         0.475589\n",
       "                False        0.475560\n",
       "Name: Stay, dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_lesion_dls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_doll_style(ax, data, yerr, title=''):\n",
    "    lightgray = '#d1d1d1'\n",
    "    darkgray = '#929292'\n",
    "\n",
    "    bar_width = 0.2\n",
    "\n",
    "    bars1 = data[:2]\n",
    "    bars2 = data[2:]\n",
    "    errs1 = yerr[:2]\n",
    "    errs2 = yerr[2:]\n",
    "\n",
    "    # The x position of bars\n",
    "    r1 = np.arange(len(bars1)) * .8 + 1.5 * bar_width\n",
    "    r2 = [x + bar_width for x in r1]\n",
    "\n",
    "    plt.sca(ax)\n",
    "\n",
    "    plt.bar(r1, bars1, width=bar_width, color=lightgray, yerr=errs1, capsize=4)\n",
    "    plt.bar(r2, bars2, width=bar_width, color=darkgray, yerr=errs2, capsize=4)\n",
    "    plt.ylabel('Stay probability', fontsize=15)\n",
    "    plt.xticks([r + bar_width / 2 for r in r1], ['same', 'different'], fontsize=15)\n",
    "    plt.yticks(fontsize=15)\n",
    "    plt.title(title, fontsize=18)\n",
    "    plt.ylim([0.40, .96])\n",
    "    plt.xlim([0, 1.6])\n",
    "\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAACsCAYAAABmUVoTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYFNXZ9/HvDQIiguyLskXRGDCAMklM0ICEqCA+KiGoGJWIMeADRty36GCMRkVNwAhRVAQxmkdFRRGVGIjEVwQUiYIgMSyC7Isgynq/f5yasWm6Z3pmepnl97muurr71Kmq09V9nz5ddeqUuTsiIiIiIhJUy3UBRERERETKEzWQRURERERiqIEsIiIiIhJDDWQRERERkRhqIIuIiIiIxFADWUREREQkhhrIUmJmNtDM3My657osItlmZjPMbFmuyyEimWFmT5rZnri0O6LfvZa5KldpmNksM1tahuUP2BdVhRrIlZCZHWlmD5vZx2a2w8w2m9lCM3vCzE6JyZdvZmfnoHxto213zva2RWKZWffoR++aIvK4mb2czXKJSPFi4jfZdGKuyygV10G5LoCkl5nlATOB3cAE4COgNnAMcCawDfhHlP024AnghRJuZiLwNLCrlMVsG217GTC/lOsQyZVTAct1IUSk0F+BqQnSS33kVEQN5MrnNuAQ4Hh336/xaWZDgealXbGZ1XX3be6+F9hbtmKKVEzuXto/hiKSGe+5+5O5LoRULupiUfkcDWyMbxwDuPs+d18ddXEouMf4xbGnpAryRq/Hm9lPoj5M24Ep0bwD+iCbWd2oj9ZsM9tgZjvNbKmZ/cHMDonJN5BvjmA/HrPtGcnWHbPsAX0/zWxZlN7JzKab2XYzW2dmI83sIDM7OHq+ysy+NrN/mtl3SrNjRSDp93BG9F080sxeNLOtZvaFmU02syPj8hacFh5oZsPMbEn03VxiZsOSbPPHZvZGtN6vzOw9MxuUJG87M3vczD4zs11mtjoqU5eYPKea2TNm9mm0vi1m9rqZdUv2fqN6Y3KUd3NUPxxqZtXM7CYz+2/0Pt4zs65lec8FcZ0gvXA9MWkHW+iytdhCl7ItZvZvM7s30f6RqsfMekbfm18kmJf2PrbRb+ZSM/tWTH2wycweNbNDzKy6md0SEzPzLEF3kCi+/hDF6U4zWxPFXasEeRtG698Y/Q6+aWbHF1HG70dl2xite7GZ3Whm1dO5LyoyHUGufP4DfNvM+rr780nyrAcuJHSVeAt4OEm+POBnwCOErhhFOQK4FHgOeArYA3QDrgOOB06L8v0TuBO4KdruW1H62mLWX5SWwBvAM8CzhFPgVxOOcncgdDH5A9AYuAZ4wcy+4+77yrBNqVwOMbPGZVxHHcKfv3eBGwl/Vi8HTjSz4919TVz+YYQzOn8hdH06HxhlZg3dfURBJjM7E5gMrAHui/KeB4wzsyPd/eaYvHnA34EawKPAh0BDQiz+CJgXZR0YpU8APuOb+P27mZ3i7gVxGfve3iTE7w3A94BLgIOBjcAPgNHRdq8BpphZG3ffVpr3XEJ/jsoyAXgAqE7Y9z1KuT6peBLF784E379sqkuoD94ErifEyCVALeBL4ARgVPT6WuBlM2vr7tsBzKwG4XftROBvwEhCV8khwKlmlufuq6O8NaO8JxB+q9+Nnv8d2ALs91tnZv9D+K1cDNwLbAa6Ar8HOhLiUtxdUyWagB8S+gY7sAR4jBBQ30mQ14HxSdbj0dQzwbyB0bzuMWk1gRoJ8v4uyvv9mLTuUdrAVNYdM28GsCwubVmU/+dx6fMIlcKLgMWkXxHlPy3Xn5Wm3E8x38Xippdjlkn0PZwR5ftjXPo5UfrYBNvcBrSMSa9J+GHbXZBOaOwtJ/zIHR6X91+EP4FHR2lGaBB/DXRM8F6rxTyvk2B+M2ADMDXJe7s2Lv35KMbmxsY+8D9R/l+X5j1H6cuAGUV8XgNj0jbFl1lT1ZiKid+nY/L1jNJ+kWAdTwJ7Uki7I1pHyxTKNSvKOzwu/aUoZmYDB8Wk943yD4pJGxKl3Rm3jrOi9Mdj0i6P0n4bl/eaKH1pTNohhINk/wCqx+W/Nsp/UlH7oqpM6mJRybj7/wO6EP5FHgb8EngIWGhmb8Wf7i3GB+4+PcXt7nL33QAWujY0iP7RFyz/gxJst6RWufv/xaXNIjQYRnsU5ZGCI2NHZ7A8UvE8DPw0yVQSf4h94e6TCUdpEo0WM8ndP4vJu4twBPQgwgW1EGK5NfCYR0eLYvLeS+gmd1aU3JlwxuRxd18QvzGPOWPi7l8WPI9O4zYiNLZnkzhW9xKOEMd6ixBjYwtiPyYdEsdYKu+5pLYCHczsuFIuLxVfovi9I6clCmdR/xyXVhAzY9x9T1w67B8z50TruDt2Be7+IuGP8NlmVnCx8NmEP5kPxG3vQcLR6linEc6mPgY0MLPGBRPfXOh4avFvr/JTF4tKyN3/TTgSi5m1IZxevRQ4GXjRzLp4ahcaLSnJds3scmAw4Uc6/s9Xg5Ksq4T+myBtc5J5BemNMlccqYA+SfZn8JvfoGJt8QO7UQAsIvyY1YltmEbp8RZGjwV/ZL8VPX6UIO+HcXkLflzfL66gZnYU4XTqaUD9uNl+4BJ87u5fx6UljDF33xzts0Qxlsp7LqkrCd3F/m1mnxKOjE0Bpri6UVUVSeM3hz5L8Dtbkt+lb0Xr2Jpg3R8BxxF+VzcRYmeVR90zCrj712b2X0I3wwIF1+BMKKLszYqYV2WogVzJuftyYIKZFfQ37gp8n3CEtTg7Ut2OmV1F6B/5OqFf1WpCV48jgPGkfkFooh/nAsm+r0WNqJFsnobpknRL9t1N9l1LlD8+b0m+pwV5i4ohzOxQQl/iOsAfgX8Tuj7sI/SdTtR3N10xlsp7TpYPEtQB7v6imbUFehMOBvQEBgFvmVnPFA8GSOVWmt+VsiprzJQ09lOtfwpeX0WI/URWlWDblZYayFWEu7uZzSY0kI/IwCYuJPQb7BV71MbMTk9UnCLWsyl6bJhg3rcIp5FEyqMGZtY8wVHkY4F1cUePAdonWEfB0Z1Po8f/RI8dEuRtH5d3cfSY9Mr1yE+Aw4FL3P3x2BlmlunT0qm8Zwj1QKI6IOFRZnffROgr+WR02vkPhAuEzwLiu19J1VPU70ppz1xk2n+AHmZWz92/iJvXnnBdwuaYvN3N7NDYo8hmdjDhvgOxddIn0eP2cnjUvVxRH+RKxsx+amYH/PExs9p806+o4JTmdhJXGKWxl9DwLfy3GpXjhgR5CwI40bYLunX0jE00s/MJP+oi5dl+33czOwf4NolvxnOBxdy2NroSfTghlgru3PcesAL4pZk1j8lbg28uqHkxSv6AcOr1EjM7oEEd01+x4OiVxc0/lcxeKwCpvWcI9cCxZnZETN5awP/GriwaLmu/LiLRNQcF3UzSVb9JxfYp4TsW/7tyMmG0pvLoBcJBzOtiE6NRbb4LvBBzfc2LhBFkhsetYyjhTFGsqYSRZ240swO6PppZbTOrW/biV3w6glz5PAA0MrOXCKdPdgCtgAGEIWImRH2UAd4BeprZ9YQfYXf3p0u53WeBu4BXzex5oF60zURHfBcSTulebmY7CP+E17n7m+6+2MymA7+OftDnEy4+OodwV6QapSyfSKZtAPqa2eGEkR8KhnlbC+QnyL8EmG1mYwnxMIAwfNrv3H0lgLvvtXCDn8nAHDN7OMp7LmH4pzvd/ZMor5vZLwlDO71rZgXDvNUndD2YRrjQbhbRkHFR14TPCDF2IaHO+G7a9kgp3nPkQcJQdtOjvDWj8sV3+6oLfB7Vd+8D6whnmoYQjq5NyeB7kQrC3bdG3QwHmtmThO6GxxCu1fk3ic/Q5NqjwEXAzdHF9QVlvhz4HLg5Ju844FfA7dH1BbMJF/j25cBrBLab2UWEUWgWm9njhN/WBoSzXX2BPqTWDbNSUwO58rmKcFrxJMIYxvUJV3kvIFwNOz4m7+WEq2xvJvzQQLiFdGncSzgiNQj4E+EH+Bngcb45Yg2Au39lZucRrjL+I2EcyJmE8SIh/BCOBi6Inr8FnAKMIZwuEimPviT0332AcIrfCI3Sq9398wT5RxP+SA4jjFSxArjS3f8Um8ndp5jZT4BbCEeNaxIudvuVu4+LyzvHzL4H/BboT7hodgNhKLV/RXm2mNlpwD3Rtg8iDIvYmxC/mWwgp/qe/2XhZiA3EeqWVYT4n0v4A1BgB6EO+Qnh6OChhMbDS8BdsSN/SJX3G8IZl7MJjcC5wBmEo6zlroHs7rvM7Kd8E8v9CH/6niYM5xY7qs1OM+tJiJWzgJ8TGsk/IcRc87h1T43qiRsIv7GNo3UvjdaR6KLgKsf2HwFLRERKysJd39q6e9sU8nYnjLTwS3cfn9GClRNV8T2LSMWmPsgiIiIiIjHUQBYRERERiaEGsoiIiIhIDPVBFhERERGJUSWOIJ9++ulOuHpVk6aqOuWc4lCTptxTHGrSlJqUhnkzs7nAY8Bf3X1zcfnLmw0bNuS6CCIVWjrqgFzHYX5+PiNGjEg6/7bbbiM/Pz97BaqEFixYUOT8MWPGMHbs2KTzBw8ezJAhQ5LO79ixY6nLVhlUhjgUqShSHQf5I8IYuvdFA7I/Crzh6p8hUlWU+zqguMbZ2rVri51f1DqqeuNMyoVyH4cilUXKfZDNrA7h7k0XAycTBm6fAIwvuJNTeZWXl+dz587NdTFEcsmKz1LMCspYB2Q6DotrIJeVGsjF02dQrEofhyIVQEpxmHIfZHf/0t0fc/duhFuoPk6409nHZvZPMxtoZgeXrqwiUt6pDhDJPcWhSHaU9iK9fXzT0XkvoTX+ELAsujWiiFRuqgNEck9xKJIhKTeQzewQM7vYzP4BfEI4xfMQ0MrdTwZaAm8Cf8lISUUkp1QHiOSe4lAkO1JqIJvZo8Aa4M/AcuAUdz/W3e9x97UA7r4J+BPQNkNlFZEcUR0gknuKQ5HsSXUUi+8C1xCGltlWRL6PgFPKXCoRKW9UB4jknuJQJEtSbSD3Az53993xM8zsIOBwd1/h7tuBmeksoIiUC6oDRHJPcSiSJan2Qf4vcHySeZ2i+SJSeakOEMk9xaFIlqTaQC5qzLiDgZ1pKIuIlF+qA0RyT3EokiVJu1iYWUegc0xSbzM7Ni7bwUB/YEkGyiYiOaQ6oGp46623uPTSS1m8eHGuiyIJKA5FcqOoPsjnALdFzx24NUm+/wK/TmehRKRcUB2QQ23btmXt2rVUr16dOnXq0Lt3b0aPHs2hhx6a1u2cfPLJFapx3LZtW8aNG0fPnj1zXZRsURyK5EBRXSzuBOoC9QindXpEr2OnWu5+lLtPz3RBRSTrVAfk2JQpU9i+fTvvvfcec+bM4Y477jggj7uzb9++HJROskRxKJIDSRvI7r47uqXldnev5u4zotex0wFX0opI5aA6oPw44ogj6NWrFx9++CEA3bt35+abb6Zr164ccsghfPrpp2zdupVBgwbRokULjjjiCG655Rb27t3Lzp07qV+/fuGyAOvXr6d27dqsW7eOGTNm0LJly8J5ixYtonv37tSvX58OHTrw0ksvFc7r3r0748aNK3w9fvx4TjrpJCA01O+99166d+9O165d6devH5988knC97Nu3TquuOIKTj75ZPr06cNzzz1XOO+3v/0tDz74YOHrOXPm8NOfhpvCXXjhhaxYsYIzzzyTQw89lHvuuQeAWbNm8aMf/Yj69evTqlUrxo8fD8DWrVu56KKLaNKkCW3atOGOO+4o/DMxfvx4unbtyvDhw6lfvz5HHnkkb7/9NuPHj6dVq1Y0bdqUJ554orAcO3fu5JprrqF169Y0a9aMwYMH89VXX5XgUywdxaFIbiRtIJtZezOrFfO8yCl7RRaRbFAdUH6sXLmSqVOncvzx3wxgMHHiRB5++GG2bdtGmzZtuPjiiznooINYunQp77//Pq+//jrjxo2jVq1a9O3bl7/+9a+Fy/7tb3+jW7duNG3adL/t7N69mzPPPJNTTz2VdevWMXr0aC644IKUumC8/vrrzJs3j5deeolZs2Zxzz33UL9+/YR5b7jhBpo1a8b06dMZOXIko0ePZvbs2cVuY+LEibRu3brwyPp1113HihUr6NWrF8OGDWP9+vXMnz+fzp1Dl91hw4axdetWPv30U2bOnMmECRN4/PHHC9c3e/ZsOnbsyMaNGxkwYADnnXcec+bMYenSpTz55JMMHTqU7du3A3D99dezZMkS5s+fz9KlS1m1ahW33357sWUuK8WhSG4U1cXiQ8KwMQXP/51kKpgnIpWL6oAcO/vss6lfvz4nnXQS3bp146abbiqcN3DgQDp06MBBBx3Epk2bePXVV/njH/9InTp1aNq0KcOHD+fpp58GYMCAAfs1kJ966ikGDBhwwPbeeecdtm/fzg033EDNmjXp0aMHffr02W/ZZGrUqMGXX37JsmXLcHeOPPJImjRpckC+NWvW8P7773PllVdSq1Ytjj32WM455xxefvnl0uwiJk2aRM+ePTn//POpUaMGjRo1onPnzuzdu5dnnnmGu+66i7p169K2bVuuvvpqJk6cWLjst771LX75y19SvXp1zj33XFauXMmtt95KrVq1OPXUU6lZsyZLly7F3XnkkUd44IEHaNiwIXXr1uWmm24q3L8ZpjgUyYGiLtI7BVgY81xEqhbVATn2wgsvJL0YrVWrVoXPly9fzu7du2nRokVh2r59+wrz9OjRg6+++orZs2fTvHlz5s+fzznnnHPAOlevXk2rVq2oVu2bYydt2rRh1apVxZa1R48enHfeedx5552sWbOGHj16cNVVVx1wUeG6des47LDDqFOnTmFaixYtWLhwYfwqU7Jy5UqOOuqoA9I3bNjArl27aNOmTdL30qxZs8LntWvXTpi2fft21q9fz44dO+jSpUvhPHdn7969pSpzCSkORXIgaQPZ3Wcmei4iVYPqgPLN7JshcVu1akWtWrXYsGEDBx10YLVerVo1+vfvz1//+leaNWtGnz59qFu37gH5Dj/8cFauXMm+ffsKG8krVqzgmGOOAaBOnTrs2LGjMP+aNWv2W/6CCy7gggsuYOPGjVx33XWMHz+eoUOH7penadOmbN26lS+//LKwkbxmzZrC7h61a9fer2/vhg0bkr7vgvf+7rvvHvBeGjduTI0aNVi+fDnt27cvfC9HHHHEAXmL07hxY2rXrs1HH31UquXLorLFYX5+PiNGjEg6/7bbbiM/Pz97BRJJItUbhYiISDnVokULTj31VK6++mq++OIL9u3bx3/+8x9mzvymPTVgwACeeeYZJk2alLB7BcAPfvAD6tSpwz333MPu3buZMWMGU6ZM4bzzzgOgc+fOPP/88+zYsYOlS5fy6KOPFi47Z84cFixYwO7du6lduzY1a9akevXqB2yjefPmdOrUiT/96U/s3LmTJUuWMHnyZHr37g3At7/9bWbNmsXWrVvZsGEDkyZN2m/5Zs2a8emnnxa+vuCCC5g+fTp/+9vf2LNnDxs3bmT+/PlUr16d/v37c/PNN7Nt2zaWL1/O/fffzy9+8YsS799q1arxq1/9iuHDh7Nu3ToAVq1axWuvvVbidVV2CxYsKHLq27cvH3zwAR988AF5eXnk5eUVvv7ggw/o27dvkcuLZEtRF+mtN7N1qU7ZLLSIZJ7qgIplwoQJ7Nq1i/bt29OgQQP69evH559/Xji/oPG7evVqevXqlXAdNWvW5KWXXuLVV1+lcePGXH755UyYMIFjjw33pRg+fDg1a9akWbNmXHzxxVxwwQWFy37xxRfcfvvtnHzyyfTq1Yv69etz8cUXJ9zO3XffzerVq+nZsyfDhw9nyJAh/PCHPwSgT58+HHPMMfTq1YvBgwdz2mmn7bfsjTfeyB133EH9+vUZOXIkrVu3ZurUqdx33300bNiQzp0788EHHwAwevRo6tSpw5FHHslJJ53EgAEDuOSSS0q1f++++27atWvHiSeeSL169ejZs2dWxo+ubHE4ZswYOnXqRKdOnZg7dy5z584tfN2pUyfGjBmT6yKKAGDunniGWT5hUPKUuHvycyY5lpeX53Pnzs11MURyqahb1CZeIM11QKbjMNNHlzp27JjR9VcG+gyKpTgso0rwHZDcSykOi+qDnJ+2oohIhaM6QCT3FIeSbuoHnpqiRrEQERERkQqkuKP4ffv2pW/fvgAMGjQIYL/rCYpbR1U5ip+0gWxmfwNudPf/RM+L5O7901oyEckp1QEiuac4lHQbM2YMY8eO3S+tU6dOhc8HDx7MkCFDsl2scqeoI8hNgBrR86aUoA+UiFQKqgNEck9xKGk1ZMgQNYBTUFQf5FNinnfPSmlEpNxQHSCSe4pDkdzI+jjI0f3i/25mO8xstZndbmYHDpZ54HJ5Zva6mW00s01mNt3MfpCNMouIiIhI1ZFyA9nMvmtmT5nZUjP7Mnp8ysxS7q1tZg2A6YRTRGcBtwNXA0UOS2NmraLlDgIuAi6Mnr9uZm2KWlZE0iMddYCIlI3iUCQ7UhrFwszOBv4G/Ad4FlhH6At1FjDXzPq7+wsprGowUBvo6+5fAG+YWT0g38zuidISOQOoGy23JSrT28AGoDegkcVFMiiNdYCIlJLiUCR7Uh3m7W7gRaC/x9xZxMxuJATpPUAqQdkLeC2uIfx0tP5uwJQky9UA9gDbY9K2R2klHnhdREosXXWASFoMHDiQli1bcscdd+S6KNmkOBTJklQbyK2AK2IDEsDd3cweBianuJ5jgTfj1rHCzHZE85I1kJ8jdMe4z8x+H6XdCmwG/i/FbYtI6aWrDpAMSeUOZqmOX9q2bVvWrl1L9erVOfTQQzn99NN58MEHOfTQQ8taTCkbxaFIlqTaB3ku0CHJvOOA91JcTwNgS4L0zdG8hNx9NXAK8DNgbTT1BU5z9/UpbltESi9ddYBUEFOmTGH79u3Mnz+f999/n7vuuisn5dizZ09OtltOKQ5FsiRpA9nMDimYgKuAy83sejP7tpk1iB5vAIYAV5Zgm4nGcLQk6QVlaUE4fTSP0E2jV/T8FTNrnWSZy8xsrpnNXb9ebWiRkkpHHaA4rPiaN2/Oaaedxvz58wHYuXMn11xzDa1bt6ZZs2YMHjyYr776CoBLLrmE6dOnA/Dee+/RqVMn3nrrLQDeeecd+vcP97BYuXIll156KT/+8Y/p1q0bN954I1988U3Pu169evHYY4/Rr18/TjzxRPbs2cOiRYs44YQTqFu3Lueeey5ff/11NndDzigORXKjqCPI24Ft0TQbOBK4C1hIuDhuIXBnlD47xe1tBuonSD+MxEeWC1xL6A7Sz92nufs0wtHkvcA1iRZw94fdPc/d85o0aZJi8UQkRpnrAMVhxffZZ5/x6quv0q5dOwCuv/56lixZwvz581m6dCmrVq3i9ttvB6BLly7MmTMHCA3kli1bMnfuXADmzZtHly5dAHB3Bg0axPTp05k8eTJr1qw54M5e06ZN48EHH2TWrFm4O8OHD+fCCy9k06ZN/PznP+e5557L1i7INcWhSA4U1Qf5EtJ/x56PCX2NC0VDuNWJ5iVzLPCRu+8uSHD3XWb2EXBUmssoIkEm6gCpIM4++2zMjO3bt9OjRw9GjBiBu/PII4+wYMECGjZsCMBNN93EgAEDuOuuu8jLy+Pee+8FQoN40KBBPP/884WvBwwYAEDr1q1p3Tqc/GvYsCEXXnjhAQ3k888/n+bNmxcuu2fPHq688krMjH79+nH//fdnZT+UA4pDkRwo6k564zOwvVeBa82srrtvi9LOBb4CZhax3HKgt5nVdPddAGZWi9DnKtmFfSJSBhmqA6SCeOGFF+jZsyczZ85kwIABbNiwgV27drFjx47CI8EQjgbv3bsXCBcBLl++nI0bN7J48WJGjRrFQw89xObNm/nwww8Ll9u4cSN333037733Hjt27GDfvn3Uq1dvv+0XNI4B1q1bR9OmTTH7ZtCiNm2qxhD4ikOR3Mj2nfTGAjuB582sp5ldBuQD98cO/RYNfP5ozHLjgMOByWZ2hpn1IQxl0wJ4OGulFxGpYrp168bAgQO55ppraNy4MbVr1+ajjz5iy5YtbNmyha1bt7J9exiBs3bt2rRv355JkybRrl07atSoQefOnZk4cSItW7akQYNwLfaoUaMwM5599lnefvtt7rzzTuIGZtivMdykSRPWrVu3X54VK1Zk4d2LSFVVkjvpnRvd3nmFma2Ln1JZh7tvBn4CVCcc+R0BPADcFpf1oChPwXLzgNMJNwuZCEwADgF+6u4fpPoeRKT00lEHSMV05ZVX8sYbb7BgwQJ+9atfMXz4cNatCx/5qlWreO211wrzdunShaeffrrwaHFeXh5PP/00eXl5hXl27NjBIYccQt26dVm7di1PPPFEkdvv1KkT1atXZ9SoUezZs4fnn3+ed999NwPvtPxTHIpkR0oNZDMbADwBLAVaAi8BL0fLfwE8mOoG3X2hu/dw99ru3sLdf+vue+PytHX3gXFpf3f3H7t7w2jq5u4zUt2uiJReOusAyYyOHTsWO5VWkyZNuOiii/jd737H3XffTbt27TjxxBOpV68ePXv2ZPHixYV5u3TpwpdfflnYQC54fcIJJxTm+fWvf82iRYvo2rUrw4YNo0ePHkVuv0aNGtx///2MHz+eBg0a8Mwzz9C3b99Sv5+KSnEokj0Wf1orYSaz9wnDrP0B2A3kuft7ZlYXeAN41t1HZrSkZZCXl+cFV1KLVFFluuNkOuqATMdhKjfKKIuyNDCrCn0GxVIcllEl+A5knD6DYqUUh6l2sTga+Fd0pHcvUA8gutDubmBoaUooIhWG6gCR3FMcimRJqg3krUCt6Pkq4Dsx8wxolM5CiUi5ozpAJPcUhyJZUtQ4yLHmAh2B1wh9nm41sz3ALuBWUr9RiIhUTKoDRHJPcSiSJakeQb4LKBhT51bgXeAh4HHCnXwuS3/RKo/8/HzMLOmUn5+f6yKKFEd1gEjuKQ5FsiSlI8ju/g7wTvR8C3BWdKOOWrHjF1dVxXWI79u3b+EV14MGDQLg0Ucf3S9PUeuoBB3ipYJTHSCSe4pDkexJtYtFIQujtzcGNrgkEronAAAYNUlEQVT7zvQXqfIZM2bMAbdR7dSpU+HzwYMHM2TIkGwXS6RUVAeI5J7iUCSzUm4gm1lv4BagS7TcHjObB/ze3V/JUPkqhSFDhqgBLBWe6gCR3FMcimRHqjcK+TXhznfbgd8AP48etwMvRfNFpJJSHSCSe4pDkexJ9QjyTcDD7h5/GHSsmY0Fbgb+ktaSiUh5ojpAKqW2bdsybtw4evbsmeuipEJxKJIlqTaQGwHPJ5n3HPCL9BRHRMop1QHl3OjRo4vNM2zYsJTW1bZtW7766is+/fRT6tSpA8C4ceN48sknmTFjRlmKKWWjOBTJklSHefsH0C3JvG7AP9NTHBEpp1QHVDF79uzhT3/6U0bWmwu52m6aKQ5FsiRpA9nM2hdMwCjgQjMbY2anmdnx0eNY4ELggWwVWESyQ3VA1XbttdcycuRItmzZknD+22+/zfe+9z0OO+wwvve97/H2228nXVevXr147LHH6NevHyeeeCJ79uxh3bp1XHXVVXTv3p1evXoxadIkAHbu3Mn3v/99Nm/eDMDDDz/MCSecwPbt2wG45ZZbuPLKKwF45ZVXOP7446lXrx6tWrXab0z5ZcuWYWY8+uijtG7dmh49egAwceJE2rRpQ6NGjfj9739f5v2UaYpDkdwo6gjyh8C/o2ka0Ar4NfAq4W4+rxIGJW8VzReRykV1QBWWl5dH9+7dGTly5AHzNm3axBlnnMEVV1zBxo0bueqqqzjjjDOSNqYBpk2bxoMPPsisWbOoVq0aV1xxBccccwxvvPEGjzzyCJMmTeJf//oXtWrVokOHDsybNw+AefPm0aJFC+bPnw/AP//5T7p1CwdR69Spw4QJE9iyZQuvvPIKY8aM4YUXXthvuzNnzmTRokW89tprLFy4kCFDhjBx4kRWr17Nxo0b+eyzz9K1yzJFcSiSA0X1QT4la6UQkfJIdUAVd/vtt9O1a1d+85vf7Jf+yiuvcPTRR3PhhRcCcP755zNq1ChmzpzJWWedlXBd559/Ps2bNwfCjZE2b97M4MGDAWjZsiU/+9nPmDZtGl27dqVLly7MnTuX7t2788knnzBo0CDmzp3L119/zZw5czj55JMB6N69e+H6O3bsyPnnn8/MmTM5++yzC9Pz8/ML+1E/++yz9OnThx//+McA/O53v+PBBx9Mw57KKMWhSA4kbSC7+8xsFkREyhfVAXLcccfRp08f/vCHP/Cd73ynMH316tW0adNmv7xt2rRh3bp1SddV0DgG+Pzzz1m/fj0nnXRSYdrevXs54YQTgHD0euTIkSxatIijjz6aE088kfz8fN555x3atWtH48aNAZg9ezY33HADH374Ibt27WLnzp38/Oc/32+7rVq12q/csa/r1KlDo0aNSrJLsk5xKJIbJbqTnpn9ADgJaAhsAma5++xMFExEyh/VAVXPiBEjOOGEE7j66qsL0w4//HCWL1++X74VK1bQoUOHpOsJN34LmjdvzhFHHMGUKVMS5u3UqRPLli3jzTffpEuXLhx11FF8/vnnvPLKK4XdKwAGDBjA0KFDefXVVzn44IO58sor2bBhQ9LttmjRgkWLFhW+3rFjBxs3bixmD5Q/ikORzEv1RiF1zGwq8P+Au4BLose3zewVMzskg2UUkRxTHVB1tWvXjnPPPZdRo0YVpvXu3ZslS5bw1FNPsWfPHp555hkWLlxY2HWhOMcddxx16tThscce4+uvv2bv3r188sknfPjhhwDUrl2b9u3b88wzz5CXlwdA586d+ctf/rJfA3nbtm00bNiQgw8+mHfffZennnqqyO3269ePl19+mVmzZrFr1y5uvfVW9u3bV9JdkjOKQ5HsSfUI8j3AD4FzgefcfZ+ZVQN+RhiU/G4gtQE2RaQiUh1QzqU6xnFp3HrrrUycOLHwdaNGjXj55Zf5zW9+w5AhQ2jXrh0vv/wy9erVS2l91atXZ9SoUdx333307t2bXbt20bZtW4YOHVqYp0uXLnz88cccd9xxha/feOON/RrhDz30EFdffTVDhw6lW7du9O/fv8gLBTt06MCf//xnBgwYwJdffslVV11Fy5YtS7o7cklxKJIl5u7FZzJbA9zq7g8nmHcZcLu7Nz9wyfIhLy/P586dm7H1L1iwIGPrhnDxiUgZWfFZilg4DXWA4rDy02dQLMVhGVWC70DG6TMoVkpxmOqNQg4DViaZtxJI7bCBiFRUqgNEck9xKJIlqTaQPwCGWOzVDkD0ekg0X0QqL9UBIrmnOBTJklT7IN9EGIz8YzObDKwFmgLnAG2BXhkpnYiUF6oDRHJPcSiSJSk1kN39TTM7HrgV+DnQAvgcmA30dfeFmSuiiOSa6gCR3FMcimRPsQ3k6ArZFsAKdz8v80USkfJEdYBI7ikORbIrlT7I1YBlhEHJRaTqUR0gknuKQ5EsKraB7O57gOWABiAvp/Lz8zGzpFN+fn6uiygVmOoAkdxTHIpkV6oX6d0N3Gxmb7n7+kwWSA40evToIuc3atSo8C5XBY9XXHFFyuvI5A0GpNJQHSCSe4pDkSxJtYF8KqHv0zIzm0e4cjb2DiPu7uemu3AiUm6oDhDJPcWhSJak2kBuDCyOe10qZtYeGE24XeYWYBwwwt33prBsX+BG4DhgBzAH+Jm7f1na8lQGU6dOZdq0afulxR5BPv300+ndu3e2iyWVS9rqABEpNcWhSJakOszbKenYmJk1AKYDC4GzgKOA+wh9oW8pZtlLgQcJ96K/FmgA9CD1Rn6l1bt3bzWAJaPSVQeISOkpDkWyJ9uNy8FAbcJ4jV8Ab5hZPSDfzO6J0g5gZo2BB4Bh7v5IzKzJGS+xiIiIiFQpqd5qGjP7rpk9ZWZLzezL6PEpM+tYgu31Al6Lawg/TWg0dytiuf7R4xMl2JaIpFGa6gARKQPFoUh2pNRANrOzgXnA8cCzwG+jx+OBudH8VBwLfByb4O4rCP2Jjy1iuR8Q+l0NMrPPzGy3mc02sx+luF0RKYM01gEiUkqKQ5HsKckwby8C/d298IpZM7uREJz3AC+ksJ4GhAvz4m2O5iXTHPg2oZ/ydcDG6HGamR3t7mvjFzCzy4DLAFq3bp1C0USkCKWqAxSHImmlOBTJklS7WLQCxsUGJITxZICHgZYl2KYnSLMk6QWqAYcCg9x9krtPA84G9gJDE27E/WF3z3P3vCZNmpSgeCKSQKnqAMWhSFopDkWyJNUG8lygQ5J5xwHvpbiezUD9BOmHkfjIcoFN0eOMgoSoH/M8oH2K2xaR0ktXHSAipac4FMmSVLtYXAU8bWY1CKdv1gFNgXOAS4HzzKzw9pfuviPJej4mrq+xmbUC6hDXNznOIsIRZotLN2Bfiu9BREovXXWAiJSe4lAkS1JtIL8bPd4F3BmTXtBgnR2Xv3qS9bwKXGtmdd19W5R2LvAVMLOI7b8M3AacAkwFMLPDgC7AyFTegIiUSbrqABEpPcWhSJak2kC+hKL7CKdqLHAF8LyZ3Q0cCeQD98cO/WZmS4GZ7j4IwN3nmtmLwKNmdgOwgXCR3m7gz2kol4gULV11gIiUnuJQJEtSvZPe+HRszN03m9lPCHfEm0Lod/wAoZEcX674f76/AO4F7gcOAf4F9HD3zekom4gkl646QERKT3Eokj1Zv02zuy8k3CK6qDxtE6RtB4ZEk4iIiIhIRqR8Jz0RERERkapADWSpEvLz8zGzpFN+fn6uiygiIiLlRNa7WIik24IFC4rN07dvX/r27QvAoEGDAHj00UdTXk/Hjh3LUEIRERGpSFJqIJtZdXffm+nCiGTKmDFjGDt27H5pnTp1Knw+ePBghgxR9/ZkVAeI5J7iUCR7Uj2CvMrMJgCPu/uiTBZIJBOGDBmiBnDZqA4QyT3FoUiWpNoH+S9AP+BDM5ttZpeZWb0MlktEyhfVASK5pzgUyZKUGsjufpu7Hwn8FFhMGIv4czObZGY9M1lAkYqgsl8EqDpAJPcUhyLZU6JRLNz9TXe/CGgODAO+DbxmZsvMLN/MDs9EIUXKu/z8fNwdd6dbt25069at8LW7V/gGcgHVASK5pzgUybzSjmKRB/wYOBbYDLwFXApcZ2aXufuTaSqfSLkwevTolPOuWrWqxMsMGzasxGXKMdUBIrmnOBTJkJQbyGbWBhgIXAS0BaYT7gv/grvvMrPqwEjC7aAVlFKlTJ06lWnTpu2XdsUVVxQ+P/300+ndu3e2i5VWqgNEck9xKJIdqQ7z9ibhX+pnwHjCFbTLY/O4+14zewr4TboLKVLe9e7du8I3gIuiOkAk9xSHUhHk5+czYsSIpPNvu+22CtHtMNUjyBuA3sAb7u5F5JsPfKvMpRKR8kZ1gEjuKQ4l54rrPtioUSNGjRoFUPgYe0a1uHWUly6HKTWQ3b1/ivl2A8uLzSgiFYrqAJHcUxxKRVBZuhyW6CI9M2sJHAMcHD/P3aemq1AiUj6pDhDJPcVh0SrLKf6KqrJ0OUy1D3Jd4G/AqQVJ0WPsKZ7qaSyXiJQjqgNEck9xmNroQJXlFL/kVqrjIN8FtAZOJgTkOUB34FHgv8CJmSiciJQbqgNEck9xKJIlqXax6A3cAsyOXq929znAP83sPuBaIKW+USJSIakOEMk9xWEKKksfWMmtVBvIzYCV0fAxXwINY+ZNBZ5Le8lEpDxRHVAM9XvMvSrwGSgOU1BZ+sBKbqXaQF4JNI6efwL0AV6LXv8A+DrN5RKR8qXK1wFVZWij8kyfgeJQJFtSbSC/AfQEJgMPAE+YWRdgJ2HQ8vsyUzwRKSdUBxRDp3Vzrwp8BopDkSxJtYF8PXAIgLtPNLPtQD+gNjAU+Etmiici5YTqgGLotG7uVYHPQHEokiWp3ihkB7Aj5vVkwj9YEakCVAeI5J7iUCR7Uhrmzcz2mtn3k8zrYmZ701ssESlPVAeI5J7iUCR7Uh0H2YqYVwPYk4ayiEj5pTpAJPcUhyJZkrSLhZm1BtrGJB1vZvG3tTwYuJgwQLmIVCKqA0RyT3EokhtF9UH+JXAb4RaWDoxJku8r4NI0l0tEck91gEjuKQ5FcqCoBvJDwLOEUzoLgAuix1i7gBXuvjMzxRORHFIdIJJ7ikORHEjaQHb39cB6ADP7FvC5u+/KVsFEJLdUB4jknuJQJDdSHeZtecFzMzsEGAQcC6wBJsTOF5HKR3WASO4pDkWyJ+koFmZ2n5ktiUurC7wH/BE4F7gV+MDMjkl1g2bW3sz+bmY7zGy1md1uZtVLsHw1M5tnZm5mfVJdTkRKJlN1gIikTnEokhtFDfN2CvBkXNo1wDHAr9y9MXA4sAz4bSobM7MGwHTChQZnAbcDVwMjSlDmS4EjSpBfREon7XWAiJSY4lAkB4pqILcF5sWl/QxY6O6PQWHfqPuArilubzDhlph93f0Ndx9LaBxfZWb1ils4amD/Hrg5xe2JSOm1Jf11gIiUTFsUhyJZV1QD+SDg64IXZtYQ+A7wZly+ZUDzFLfXC3jN3b+ISXua0GjulsLyvwP+Bfw9xe2JSOllog4QkZJRHIrkQFEN5CVA95jXBf19X4vL1xTYlOL2jgU+jk1w9xWEe8sfW9SCZtaRMB7kNSluS0TKJhN1gIiUjOJQJAfM3RPPMBsIPEIYlHwtcAWwDfiOu++OyfcXoI27n17sxsx2A9e6+x/j0j8jXIF7UxHLzgRmu/t1ZtaWcMegM9395ST5LwMui15+G1hcXPnKscbAhlwXooqr6J/BhlRiNFY66gDFoaRZRf8MFIdlV9G/A5VBRf8MUorDosZBHm9mLYD/BeoTrpj937iAbEK42K4kF9klapFbkvSC7ZxHCOozU96I+8PAwyUoV7llZnPdPS/X5ajKquJnkI46QHEo6VQVPwPF4f6q4negvKkqn0GR4yC7+13AXUXMX0/J+jxtJgR4vMOALYkWMLMawL3A3UA1M6sPFFzQV8fM6rr7thKUQURSlIE6QERKSHEokn1F9UHOhI+J62tsZq2AOsT1TY5RB2gJ3E9oYG8GPojmPQ28n5GSioiIiEiVlNKd9NLoVeDauKO+5wJfATOTLLOdMA5krObAX4GbOPBK3sqoUpwaq+D0GYi+A7mnz0D0Hci9KvEZJL1ILyMbC+MYLwQ+JHSZOJJwZPiP7n5LTL6lwEx3H5RkPW0p5iI9EREREZHSyOoRZHffbGY/AR4EphD6HT8A5CcoV8q3nxYRERERSZesHkEWERERESnvsn2RnkhamdlxZuZm1j167WY2NGZ+NTP7s5mtjeblR+lnmdkiM9tlZstyUvgkzKx/NPapSIWgOBTJPcVhemX7Ij2RTPshoX96gb7A5cAgQv/3z8ysOjCBcNHor4Avs13IYvQnDMQ+PsflECktxaFI7ikOy0ANZKlU3P2duKRjgc3u/lhBgpm1JIyl/ZS7zyrL9qJxuve5+96yrEekMlEciuSe4rCM3F1TBiagAzAN2ET4R7aIcPcjgDOAN4B1wBfAO8CpccvnE27l+ANgLmEovFnAt4CmwAuEIfAWAT0SbP9S4CNgJ7AcuC7X+yRN+/VyYGW0T6cAPyXchbF7NN+BodHzGdHr2GlggrT8KH814AZgabTflgAXx21/BvAs4bat/wH2Aq2ieccBrxBuA7sN+D+gecyy3QvKGs3bDnwKXB6TZ3yy8mkq9XdGsZj+fao41FTS74ziMP37VHGYyf2b6w+4sk7Rl+UVoDfwk+iLfEM0byhwBXBa9IW+P/pidY1ZPh/YQbgpygXA2cCKqEL4O3ANcCowHdgIHBKz7LXAbuD30fpviL7gQ3O9X8q4T8+KAmRMtO/ujCqHZBVCe2AcYbSUE6OpGXBOlO/qKK1llP/PUZBeB/QkDEW4F+gTU4YZwOeEG9T0iz7fekA7YGv02ZwN/IxwCmsO31wMW1AhfALcEn02j0Vp34/yHEUY2/u9mDK3zPW+r8iTYlFxqDjM/aQ4VBxWtDjM+YdcGSdCfxkHvptC3mqEri6vAY/FpOdH6+gWk3Z5lHZrTFr7KK1X9Lpe9KW+LW47twNrgOq53j9l2K/vAq/GpT2SrEKI2Y8b4pZpG+WLDfR2wD4O/Ic8AZgT83oG4chF87h8E4HFQM2YtKOjCuWM6HVBhXB7TJ4awHrgDzFpzwIzcr2/K8OkWMzIPlUcairpd0ZxmP59qjjM8KRRLDJjE+Gf3FgzO9fMmsbONLOWZvaEma0C9hD+2Z4KHBO3nl3AWzGvl0aPbyZIOyJ6/CHh9tz/Z2YHFUzRMs0It+2ucKILCY4HXoyb9XyaNvETQoUwOW6//R3oHG2/wDx3XxO3fE9gMrAvZtn/AsuAvLi8rxc8cffdhH/QFfJzqQAUi2mkOJRSUhymkeIwO9RAzgB330cI7jWEUwZrzOwtMzvezKoBLwE/Am4l3Eb7e4QrSA+OW9W2aF0FdkWPW2K2VZBWsGzj6PEjQiVTMP0jSm9VtneXM00IRxXWxaXHvy6txoSb02xl//02Ptpui5i8a5Msf33csrsJd4uM3+db4l7v4sDPXtJAsZh2ikMpMcVh2ikOs0CjWGSIu38M/Cy6qvNkQv+dVwinFY4nnP6ZVpDfzGqnadObosc+JP7iLk7TdrJtPeHIQtO49PjXpbUpWn9Xwj/neLEVjydZfjKhj1e8DWUunZSaYjGtFIdSKorDtFIcZoEayBkWnTJ408zuB57im39eOwvymFkbwhdxQRo2+f8IfYIOd/dX0rC+csHd95rZfMKFCWNjZvVN0ybeJPxjPszd3yjF8n8nXLU7z6OOU2VQbv5BVyaKxbJTHEpZKQ7LTnGYHWogZ4CZdQRGAs8Qhi1pQDjd8AFh+JrPgPvM7LdAXWAEsCod23b3LdHdcf4UVTL/JHSlOQY4xd3PScd2cuRO4HkzG0P4d9oNOD0dK3b3xWY2FnjazO4hDCN0MGFoomPc/dJiVpFPuGjiFTN7jPAv+QjClbnj3X1GCYrzMXCWmZ1N+K6sdvfVJXk/EigWM0JxKCWiOMwIxWGGqQ9yZqwhnMq5mdCP6iHC2Iz/4+47Cf/y9hCuzvwdcBcwM10bd/d7COMS9iJ04v8rYVict4parrxz98nAMOBMwpiXxxPuCJQu/0v4PC4CphL6W51BqFCLK9sSwhA0O4CHCZ/7CMJRkaVFLJrIQ4QLFx4jDItzWQmXl28oFtNMcSiloDhMM8Vh5lnZj36LiIiIiFQeOoIsIiIiIhJDDWQRERERkRhqIIuIiIiIxFADWUREREQkhhrIIiIiIiIx1EAWEREREYmhBrKIiIiISAw1kEVEREREYvx/roPGSyZG1coAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x180 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize= (10,2.5), sharey=True)\n",
    "\n",
    "plot_doll_style(axes[0], list(mean_lesion_hpc), list(sem_lesion_hpc), title='Striatum')\n",
    "plot_doll_style(axes[1], list(mean_lesion_dls), list(sem_lesion_dls), title='Hippocampus')\n",
    "plot_doll_style(axes[2], list(mean_full), list(sem_full), title='Full model')\n",
    "\n",
    "\n",
    "\n",
    "leg = axes[1].legend(['Reward', 'No reward'], fontsize=12, frameon=False, handlelength=.7)\n",
    "leg.set_title('Previous outcome', prop={'size': 12})\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig(os.path.join(figure_location, 'StayProbability_DeterministicTask.pdf'))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PreviousReward  SameStart\n",
       "1.0             True         0.742813\n",
       "                False        0.551995\n",
       "0.0             True         0.632391\n",
       "                False        0.525771\n",
       "Name: Stay, dtype: float64"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_lesion_hpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}