13 Chapter 13/03 Oceanic tools, auto-correlation with distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from io import StringIO\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import Set2, Blues\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize(alpha=0.05):\n",
    "\n",
    "    def _summarize(ser):\n",
    "        hpd_ = pm.hpd(ser, alpha=alpha)\n",
    "        return pd.Series({\n",
    "            'mean_': ser.mean(),\n",
    "            'lower': hpd_[0],\n",
    "            'upper': hpd_[1],\n",
    "        })\n",
    "\n",
    "    return _summarize\n",
    "\n",
    "\n",
    "def agg_prediction(trace_var, source_data, group_columns, aggregate_fc):\n",
    "    df = pd.DataFrame(trace_var).stack().reset_index().drop('level_0', axis=1)\n",
    "    merged = pd.merge(\n",
    "        df, source_data[group_columns],\n",
    "        left_on='level_1', right_index=True\n",
    "    )\n",
    "\n",
    "    return (\n",
    "        merged.groupby(group_columns)[0].apply(aggregate_fc)\n",
    "        .unstack().sort_index().reset_index()\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>culture</th>\n",
       "      <th>population</th>\n",
       "      <th>contact</th>\n",
       "      <th>total_tools</th>\n",
       "      <th>mean_TU</th>\n",
       "      <th>high_contact</th>\n",
       "      <th>log_pop</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Malekula</td>\n",
       "      <td>1100</td>\n",
       "      <td>low</td>\n",
       "      <td>13</td>\n",
       "      <td>3.2</td>\n",
       "      <td>0</td>\n",
       "      <td>7.003065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Tikopia</td>\n",
       "      <td>1500</td>\n",
       "      <td>low</td>\n",
       "      <td>22</td>\n",
       "      <td>4.7</td>\n",
       "      <td>0</td>\n",
       "      <td>7.313220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Santa Cruz</td>\n",
       "      <td>3600</td>\n",
       "      <td>low</td>\n",
       "      <td>24</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0</td>\n",
       "      <td>8.188689</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Yap</td>\n",
       "      <td>4791</td>\n",
       "      <td>high</td>\n",
       "      <td>43</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1</td>\n",
       "      <td>8.474494</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Lau Fiji</td>\n",
       "      <td>7400</td>\n",
       "      <td>high</td>\n",
       "      <td>33</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1</td>\n",
       "      <td>8.909235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Trobriand</td>\n",
       "      <td>8000</td>\n",
       "      <td>high</td>\n",
       "      <td>19</td>\n",
       "      <td>4.0</td>\n",
       "      <td>1</td>\n",
       "      <td>8.987197</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Chuuk</td>\n",
       "      <td>9200</td>\n",
       "      <td>high</td>\n",
       "      <td>40</td>\n",
       "      <td>3.8</td>\n",
       "      <td>1</td>\n",
       "      <td>9.126959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Manus</td>\n",
       "      <td>13000</td>\n",
       "      <td>low</td>\n",
       "      <td>28</td>\n",
       "      <td>6.6</td>\n",
       "      <td>0</td>\n",
       "      <td>9.472705</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Tonga</td>\n",
       "      <td>17500</td>\n",
       "      <td>high</td>\n",
       "      <td>55</td>\n",
       "      <td>5.4</td>\n",
       "      <td>1</td>\n",
       "      <td>9.769956</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Hawaii</td>\n",
       "      <td>275000</td>\n",
       "      <td>low</td>\n",
       "      <td>71</td>\n",
       "      <td>6.6</td>\n",
       "      <td>0</td>\n",
       "      <td>12.524526</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      culture  population contact  total_tools  mean_TU  high_contact  \\\n",
       "0    Malekula        1100     low           13      3.2             0   \n",
       "1     Tikopia        1500     low           22      4.7             0   \n",
       "2  Santa Cruz        3600     low           24      4.0             0   \n",
       "3         Yap        4791    high           43      5.0             1   \n",
       "4    Lau Fiji        7400    high           33      5.0             1   \n",
       "5   Trobriand        8000    high           19      4.0             1   \n",
       "6       Chuuk        9200    high           40      3.8             1   \n",
       "7       Manus       13000     low           28      6.6             0   \n",
       "8       Tonga       17500    high           55      5.4             1   \n",
       "9      Hawaii      275000     low           71      6.6             0   \n",
       "\n",
       "     log_pop  \n",
       "0   7.003065  \n",
       "1   7.313220  \n",
       "2   8.188689  \n",
       "3   8.474494  \n",
       "4   8.909235  \n",
       "5   8.987197  \n",
       "6   9.126959  \n",
       "7   9.472705  \n",
       "8   9.769956  \n",
       "9  12.524526  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('../data/Kline.csv', sep=';')\n",
    "\n",
    "data = data.assign(\n",
    "    high_contact = (data.contact == 'high').astype('int'),\n",
    "    log_pop = np.log(data.population)\n",
    ")\n",
    "\n",
    "n_cultures = data.shape[0]\n",
    "\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>Malekula</th>\n",
       "      <th>Tikopia</th>\n",
       "      <th>Santa Cruz</th>\n",
       "      <th>Yap</th>\n",
       "      <th>Lau Fiji</th>\n",
       "      <th>Trobriand</th>\n",
       "      <th>Chuuk</th>\n",
       "      <th>Manus</th>\n",
       "      <th>Tonga</th>\n",
       "      <th>Hawaii</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>island</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Malekula</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.6</td>\n",
       "      <td>4.4</td>\n",
       "      <td>1.2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.2</td>\n",
       "      <td>2.8</td>\n",
       "      <td>1.9</td>\n",
       "      <td>5.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tikopia</th>\n",
       "      <td>0.5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3</td>\n",
       "      <td>4.2</td>\n",
       "      <td>1.2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.9</td>\n",
       "      <td>2.7</td>\n",
       "      <td>2.0</td>\n",
       "      <td>5.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Santa Cruz</th>\n",
       "      <td>0.6</td>\n",
       "      <td>0.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.9</td>\n",
       "      <td>1.6</td>\n",
       "      <td>1.7</td>\n",
       "      <td>2.6</td>\n",
       "      <td>2.4</td>\n",
       "      <td>2.3</td>\n",
       "      <td>5.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Yap</th>\n",
       "      <td>4.4</td>\n",
       "      <td>4.2</td>\n",
       "      <td>3.9</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.4</td>\n",
       "      <td>2.5</td>\n",
       "      <td>1.6</td>\n",
       "      <td>1.6</td>\n",
       "      <td>6.1</td>\n",
       "      <td>7.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lau Fiji</th>\n",
       "      <td>1.2</td>\n",
       "      <td>1.2</td>\n",
       "      <td>1.6</td>\n",
       "      <td>5.4</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.2</td>\n",
       "      <td>4.0</td>\n",
       "      <td>3.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>4.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Trobriand</th>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.7</td>\n",
       "      <td>2.5</td>\n",
       "      <td>3.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.8</td>\n",
       "      <td>0.8</td>\n",
       "      <td>3.9</td>\n",
       "      <td>6.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Chuuk</th>\n",
       "      <td>3.2</td>\n",
       "      <td>2.9</td>\n",
       "      <td>2.6</td>\n",
       "      <td>1.6</td>\n",
       "      <td>4.0</td>\n",
       "      <td>1.8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.2</td>\n",
       "      <td>4.8</td>\n",
       "      <td>5.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Manus</th>\n",
       "      <td>2.8</td>\n",
       "      <td>2.7</td>\n",
       "      <td>2.4</td>\n",
       "      <td>1.6</td>\n",
       "      <td>3.9</td>\n",
       "      <td>0.8</td>\n",
       "      <td>1.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.6</td>\n",
       "      <td>6.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Tonga</th>\n",
       "      <td>1.9</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.3</td>\n",
       "      <td>6.1</td>\n",
       "      <td>0.8</td>\n",
       "      <td>3.9</td>\n",
       "      <td>4.8</td>\n",
       "      <td>4.6</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Hawaii</th>\n",
       "      <td>5.7</td>\n",
       "      <td>5.3</td>\n",
       "      <td>5.4</td>\n",
       "      <td>7.2</td>\n",
       "      <td>4.9</td>\n",
       "      <td>6.7</td>\n",
       "      <td>5.8</td>\n",
       "      <td>6.7</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Malekula  Tikopia  Santa Cruz  Yap  Lau Fiji  Trobriand  Chuuk  \\\n",
       "island                                                                       \n",
       "Malekula         0.0      0.5         0.6  4.4       1.2        2.0    3.2   \n",
       "Tikopia          0.5      0.0         0.3  4.2       1.2        2.0    2.9   \n",
       "Santa Cruz       0.6      0.3         0.0  3.9       1.6        1.7    2.6   \n",
       "Yap              4.4      4.2         3.9  0.0       5.4        2.5    1.6   \n",
       "Lau Fiji         1.2      1.2         1.6  5.4       0.0        3.2    4.0   \n",
       "Trobriand        2.0      2.0         1.7  2.5       3.2        0.0    1.8   \n",
       "Chuuk            3.2      2.9         2.6  1.6       4.0        1.8    0.0   \n",
       "Manus            2.8      2.7         2.4  1.6       3.9        0.8    1.2   \n",
       "Tonga            1.9      2.0         2.3  6.1       0.8        3.9    4.8   \n",
       "Hawaii           5.7      5.3         5.4  7.2       4.9        6.7    5.8   \n",
       "\n",
       "            Manus  Tonga  Hawaii  \n",
       "island                            \n",
       "Malekula      2.8    1.9     5.7  \n",
       "Tikopia       2.7    2.0     5.3  \n",
       "Santa Cruz    2.4    2.3     5.4  \n",
       "Yap           1.6    6.1     7.2  \n",
       "Lau Fiji      3.9    0.8     4.9  \n",
       "Trobriand     0.8    3.9     6.7  \n",
       "Chuuk         1.2    4.8     5.8  \n",
       "Manus         0.0    4.6     6.7  \n",
       "Tonga         4.6    0.0     5.0  \n",
       "Hawaii        6.7    5.0     0.0  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances = (\n",
    "\"\"\"Malekula,0.0,0.5,0.6,4.4,1.2,2.0,3.2,2.8,1.9,5.7\n",
    "Tikopia,0.5,0.0,0.3,4.2,1.2,2.0,2.9,2.7,2.0,5.3\n",
    "Santa Cruz,0.6,0.3,0.0,3.9,1.6,1.7,2.6,2.4,2.3,5.4\n",
    "Yap,4.4,4.2,3.9,0.0,5.4,2.5,1.6,1.6,6.1,7.2\n",
    "Lau Fiji,1.2,1.2,1.6,5.4,0.0,3.2,4.0,3.9,0.8,4.9\n",
    "Trobriand,2.0,2.0,1.7,2.5,3.2,0.0,1.8,0.8,3.9,6.7\n",
    "Chuuk,3.2,2.9,2.6,1.6,4.0,1.8,0.0,1.2,4.8,5.8\n",
    "Manus,2.8,2.7,2.4,1.6,3.9,0.8,1.2,0.0,4.6,6.7\n",
    "Tonga,1.9,2.0,2.3,6.1,0.8,3.9,4.8,4.6,0.0,5.0\n",
    "Hawaii,5.7,5.3,5.4,7.2,4.9,6.7,5.8,6.7,5.0,0.0\"\"\")\n",
    "\n",
    "distances = pd.read_csv(StringIO(distances), header=None)\n",
    "distances.columns = ['island'] + list(distances[0].values)\n",
    "distances = distances.set_index('island')\n",
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [soc_a, soc_bar, b_pop, a]\n",
      "Sampling 2 chains: 100%|██████████| 7000/7000 [00:09<00:00, 764.39draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_varying_intercepts:\n",
    "    a = pm.Normal('a', mu=0, sd=1)\n",
    "    b_pop = pm.Normal('b_pop', mu=0, sd=1)\n",
    "    soc_bar = pm.Exponential('soc_bar', lam=1)\n",
    "\n",
    "    soc_a = pm.Normal('soc_a', mu=0, sd=soc_bar,\n",
    "                      shape=n_cultures)\n",
    "\n",
    "    tools_i = pm.Deterministic(\n",
    "        'tools_i',\n",
    "        pm.math.exp(\n",
    "            np.log(data.total_tools.mean())\n",
    "            + a\n",
    "            + soc_a[data.index.values]\n",
    "            + b_pop * (data.log_pop.values - data.log_pop.values.mean())\n",
    "        )\n",
    "    )\n",
    "\n",
    "    obs = pm.Poisson('obs', mu=tools_i,\n",
    "                     observed=data.total_tools.values)\n",
    "\n",
    "    m_varying_intercepts.name = 'varying intercepts'\n",
    "    m_varying_intercepts.trace = pm.sample(2000, tune=1500, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEfCAYAAABiR+CGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5xVdb3/8ddIaMCUGo6ahXKg/JxTdPInYmmleIGUVBJNrVQwvFXkDxWPWGqaZBpJ2vGCCohKHs1LooIXNJiDqQeiLOngB2QcGLnIcBEZNhg5c/5Ya+OecfaePbMva62Z9/Px4MHea3/X2p/Ze+39Wd/vd6/1qWhqakJERKRYdok6ABER6VyUWEREpKiUWEREpKiUWEREpKiUWEREpKiUWEREpKg+EnUAIiJdgZk1AcuBf4aLPgLMB37k7g2ttJ8OvOHuE9r5POe5+90FhlsQJZYuwsx2AW4HBgCNwPnu/nq48w4ENoRNJwLPAzOBTwAXu/sfw23MBMa4e10r2/8s8Btgr3DRK8A4d3+vnXFWAovdva+ZPQicA9wBPOLuT2W06xsuOySPbY4CBrj7uBxtTnX3R9oTax7P25c8Y5QuY7C7vwVgZrsBDwFXAZcXY+Nm1o3gMxxpYtFQWNfxTWAPd/8qMAa4OeOxK9x9cPhvFnAkwZHUCOD7AGZ2PPC3LEmlG/AocKO7DwIODR+6ukW7du1v7n6Gu29rzzodYWa7Apfk2VafGSmK8KDrGeDzOZp9wsyeMrNVZvZCeOCFmR1qZgvN7HUz+7uZHRu2nwPsHi7/lxL/CVmpx9J1fBb4E4C7/83MvhgmhNbsBawN/1WF7S4mSDStGQoscfd54fabzOw/gMbwqP23wGZgspltBK4H3gdWAqOBnsAjQA+ChAaAmdUS9LAATjCzi8PYzgE2ZrT7asttunt6uKEZM/tv4FngWIIe2QnAeOALZnY7cBFBD6k/sBtwpbvPNbPlwBPAO2b2XXc/MNzeSOCLwD3AZOAfBD3Cb7V43suBU4BuwFPu/tMsr6V0EWa2B/Btgv0xm6HAEcC7wCLgZOB+4C5gorv/1szOJtj3PgN8j2D47F9LGXtbdPTVdSwGhphZNzMbCOzNB8NWPzKzajN7yMz2Ivhy/gxgQB3BzvoQcIWZ3WNmLYd2DHg1c4G7b8sYBjsIONPdnwBuAYa7+1HAOoIv4DMJekNfA/6SJf4mdz8GuAL4SYvHWttmNo3Au2HbZwiS5cQgZP8BcAaw1t2PBobzQc+uG/CMu18L1JlZ+ihzOEFS3Bu4NNzufOC7LZ53HPAVdx8IbFTPp8uaF/YmaoBaYB7B/pfNbHffGB4o/RX4dLh8EPBgeHs+0K804XaMdu6uYzbBl/9/E3yRLwfeIzj6+Ym7HwksBH4GvETwRfmfwJ0ER0mvE3wpjwGubLHt7gRfvNksd/eNZrYnQRL6vZnNAwYTfFA+B/xP2LY6yzbmhv8vDLcBQI5t5pLuFdUBe7R4bBBwcritR4Ae4VBZ+rkBHgNONLOPEgxjvEwwR3VtuN53gd4ttvsE8JyZXQDc5+6NbcQondPgsDdxKMHn6QF33wEQDnW9bmavZ7R/N+N2Ix98zs4CFpqZEwx/VZQ+9PxpKKyLcPcmwglCM/sIcJq7vwO8kNFsFjA5/NIbHba9huCIan9ghbtvNbOPtdj834ELMxeEX7qfARoIhofS1rj74BZtbwPSV0PN9gFpynK71W22IXOYrLXn+4W7/7ZFjPDB3/F74HcEvcBnw6G/WwjmmGaHw149Mtd399FmNgA4jeALYYC7b29HzNKJuPt6M/sN8EuCXi9hj7xN4fDyZOBL7v4XM/sMsKxUsXaEeixdhJn9u5lNDe9+E3guXP47M/v3cPlXCb4s0+vsB3zW3ecC64FPm1lPgiOnTM8D/c3spHC9CuAG4DuZjdx9E8G8y4Cw3Y/C53aCX6YBHJ3lTzgi/P9LwJI8ttkejXxwkPU/hB90M9vbzH7esrG7ryZIbt8m6NVAMF/zRphQTwDSvRzMbHczu8rdF7v71QQ9xap2xiidz03A4WZ2ZDvX253ggG2JmXUnGEVI/6JyB7BLKwd/ZaXE0nW8BnzEzP6HYII6PZx1KzAlHMI5CbgmY52rMu7PI+i+zyP4WfFO7v4PYBhwvpktIBgaSoXrtzQamGpmLwJfIUgq9wFfNrMXCIa1Wu6X3QDMbBYwAbguj222xxpgVzN7iKAn8q6ZvQQ8BbyYZZ0nCH49l378FoIhsoeAXxEMh30cwN03E/wI4mUz+yNBL+dDv66TrsXdtxAcgP0qPBjL198IhraXEhwIPUIwfP0cwb78IvCmmR1e3IjzV6F6LCIiUkzqsYiISFEpsYiISFEpsYiISFEpsYiISFEpsYiISFF12hMkFy1apJ+7iSTIwIEDIzl7XN8VhWntfeu0iWXgwIE5H1+yZAn/9m//VqZo8qOY2nb33XezZs0arr766rYbl1HcXidIXEyRXpIk2/dFHF/Dcsv1GixatKjV5Z02sUjndN5557FkyZK2G4pIZDTHIiIiRaUeiyTK4MGDSaVSLFiwIOpQRCQLJRZJlFGjRrF69eqowxCRHJRYpKxWbkgx+t6F1NRvpV9VL6aOHMT+vXvmvf6oUaM0xyISc0osUhTn3LOAuV7frnWWrWvgiIlz22x3lFVxzzmHArBjxw527NjRoRhFpDwSmVjM7Au0qC/u7htzr9V1Df11NUvfbihgCzVFi6Uj5no9fcfPar7wgY5fdf7AfSp57uL2lsCQriTds15e30D/qnXt7ll3dYlMLHxQX/yVsMLhdwnK6EorCvkSzfd3/B3pseQrs8cyY8YMVq1axeWXX16S5xIBdiaVxiZYXt/A6HsXMucSHYzkK6mJZQNwo5ntBnwK+G3LBrW1tWzbti3rBrZv3x67sfpyxHThzDpWvNPeoaQ49Vj2BPbkjpY9mA44YI/uTB7ep+DtQNfdn9orW0xxOwmxpn4rjeH5+I1NwX3JX1ITS8764gB9+/bNuYE4nlFbjpiq27n9Ysc0ZFL1ziPBXSqgf1Vlu44EU6kUr7/+OgcffHDRYiqGrro/tVccY2pNv6pezfbTflW9og4pUZJ6gmTW+uISb1NHDqJ/VSXdKiroX1XJ1JGD2rX+sGHDuPDCC0sUnUggvZ+mD37au592dUnusTwGLCeoL/4bM3vQ3f8WbVjSlv179yxorPr73/8+q1atKmJEIh+W3k+T0sOKm0QmFnefAkzJWDQzqlikvE4//fTYzRuISHNJHQqTLmrz5s1s2bIl6jBEJIdE9lik6xo+fLiuFSYSc0oskigXXXQRb731VtRhiEgOSiySKCNGjNAci0jMaY5FEmX9+vVs2rQp6jBEJAf1WCRRTj31VM2xiMScEoskyqWXXkpdXccvQCkipafEIoly4oknao5FJOY0xyKJsnbtWurrS3MVZREpDvVYJFHOOOMMzbGIxJwSiyTK+PHjWblyZdRhiEgOSiySKMcdd1zs5lhWbkhxweN1rLrvTfpV9VK1QenySppYzGw68Ii7P1XK55Guo66ujjVr1kRyxdl8qmQuW9fAERPntvpYZiVMkc5MPRZJlLPOOqsocyxDf13N0rcbihRVfppXwmyfA/epLKjEtLSPat4XphyJ5QQzuxjYCzjH3f/cWiMzuw74NLBf+P8l7v6smZ0GXArsABa5+//P1rYMf4tE7Morr2TFihUFb6cjX9L59FhyUY8lOVTzvjDlSCxN7n6MmQ0DfgKckqXd+0CVu3/dzD4H3G9mfwR+ARzk7lvM7EkzO7q1tkCzxKKa98URt5g+9alP0bt37w7HdOHMOla8s6PIUeWnoz2WA/bozuThfdq9XtzeO1DN+66iHIklPeC8EPhlG22rAdz9f82sD3BgcNfTBTjmAwdladuMat4XR9xiqqmpoa6ujqFDh3Zo/eoS/C1DJlXzxroGmmBnKds4HN3G7b2DeMbUGtW8L0w5TpBsynI7n/WagIqMZRVAYxvPIZ3Y9773Pa688sqow2hm6shB9Nm9O90qKlQfvZNQzfvClKPHcgTwMPAloK1++deAiWb2eWAFsBQ40Mw+DmwBjgImAMe20la6gGuvvbYocyzFtH/vntz5zT6JOBKX/KjmfWFKnVi6AZjZLOBTwFlttG8wsznAvsBYd99qZuOBOQTzKtXu/qKZHduybcn+AomVI488MnbzBiLSXEkTi7u3lUhaesndb22xjYcJejxttpXOz9158803dRQpEmNlP4/FzB4DPtFi8WbgL+WORZLnggsuIJVKMWzYsKhDEZEsyp5Y3H1EEbZxTRFCkQS6/vrrqa2tjToMEclBZ95Lohx++OHsueeeUYchIjmoHoskyuLFi1m2bFnUYYhIDuqxSKKMGTOGVCrFSSedFHUoIpKFEoskysSJE3nzzTejDkNEclBikUQZNGgQlZWVUYchIjlojkUS5dVXX9UJkiIxpx6LJMrYsWNJpVKMGFHwr9ZFpESUWCRRbr75ZmpqaqIOQ0RyUGKRgqWr7dXUby15zfeDDjqI3XbbrSTbFpHiiHViMbO+wGsEF588IawceTTwUeACgsvl3wZ8xN0PiSzQBCis+mH+PYRcNd/bq7WKiwsXLtS1wkRiLtaJJeRhUjkKOMTdv2JmA4Db3f0IMzsDeKSUAZS3PrqGedKyV1zsxX/8uWO144tNtehFPizvxGJmQwlqobwPrAW+TVDHfhqwG/BPYLS715rZpcDJQHdgvLu3eghrZp8CZoR3uwMj3X15lhCOAmYCuPtiM9vPzEoz3tJCub444lj7IVtMmcNfu+wC/2xsoqmp9BUUFy9eTE1NjU6QlJJK79/L6xvoX7WupMO7nVF7eixjgHHu/t9m9i2gCrgWuMfd/8vMTgF+ZmY/Bb4FHA70A67gg/LELX0SuNHdnzGzUcAPgEtztP1rxv16YB9UPTKrwoa/MuXuRb3//ge3G5uC4bCO1HbP1NowGMCAAQPo1q1bQdsWaUs6qTQ2wfL6BkbfuzAW5aaToj2J5TFgspnNAB5y9zozOwQYHz4+H/gp8P+AP7l7I/AGMDrHNtcDV5rZFQSX0l+Uo+0/WtyvIEdSqa2tZdu2bVk3tn379jbPh7hwZh0r3tmRs03xaSgMcg2DpcXxdWoe0wF7dGfy8D4RxZLfPl5u2WKKW0+9pn4rjeG3S2NTcF/yl3dicffpZvYs8E3gBTP7Js1r0mfWo8/3xMtrgefc/XYzOx04PkfbNcDeGfergLcJei0f0rdv35xPnM+wU3WZd/YkDYUNmVS984iu1MNfmQYPHkwqlWLBggUlf672SNJ7F6U4xtSaflW9mu3f/ap6RR1SouR95r2ZXQWk3P0OYBbwBWAhcEzY5GjgTwQFu75sZt3MbB8zezzHZj8BvGFmFcApwK452j4NDA9jORiocffsXRIpqakjB9G/qpJuFRX0r6pk6shBZXneO++8k2uuuaYszyVdV3r/Th80lWv/7izaMxS2AnjezLYCDXwwdzLNzC4AtgPfc/dVZvZb4EWCxPWTHNu8A7gFqAV+A9xlZkPcfU7Lhu6+yMz+amZ/JvyhQDtilyLbv3fPSMaczYzGxsa2G4oUIL1/J6WHFTftGQq7D7ivxeJ3gaGttL0JuCmPbc4GZmcsyjkg7e6XA5e3Gax0WtXV1axYsUIfdpEYK8t5LNnq3Lv78PxWt6fc/YRWHhhIcIKkdBE//elPSaVSnH322VGHIiJZlCWxdLTOvbvXAh/L8fgi4MsdDEsSaNq0abzxxhtRhyEiOSThzHuRnfr168d7770XdRgikoPqsUiiPP/887z00ktRhyEiOSixSKJMmDCBO++8M+owRCQHDYVJotx///0sW7Ys6jBEJAclFkmUPn360NBQritNi0hHaChMEuWZZ55h/vz5UYchIjkosUii3HDDDUyZMiXqMEQkBw2FSaI8+OCDLF26NOowRCQHJRZJlH333ZdNmzZFHYaI5KDEIony5JNPUldXt/NaYZmVLPtV9VKlP5EYiH1iMbPBwMMEZYnvD29/z92fMrNvAD8H3nD3U6OLUjqiYxUudwEO4JetFAFbtq6BIyZmK1aav2zVK0UkP7FPLKFq4BfAJOCP6YXuPiu8jP+YqAKLo6G/rmbp28X8SW4cqzWWTtvVK7OJ/nU6cJ9KnrtYJXQLpZr3hUlKYoGgguQIYGrUgcRdMb9Yil2PomO9lPLqSI9FdTs6F9W8L0xiEou7pyAo9JSPYtS8L7dCY7pwZh0r3tlRxIjSoj8SL6ek9VgO2KM7k4d/uJRRkvbxuCVl1bwvTGISS3sVo+Z9uRUaU3UJ/p64vU6ZNe+HTKpuVpe8f1VlZEeVcXudQDEVQjXvC6MTJCVRZs6cya233gp8UJe8W0WF6pJLUanmfWE6bY9FOqfdd9+dj30sqP2WrksuUmyqeV+YxPRYzOwbZjYPOA74hZk9F3FIEoGHHnqIp59+OuowRCSHxPRY3H0W0JEZVelE7rjjDlKpFJdccknUoYhIFknpsRxpZh+68mB4guTNEcQjEZk9ezaTJ0+OOgwRySH2PRZ3nwdUZXlMvZgupmfPnvTo0SPqMEQkh9gnFpFMM2bMYNWqVZpQFYmxpAyFiQAwZcoUHn300ajDEJEc1GORRJkzZ07sziYXkebUY5FE6d69O927d486DBHJQT0WSZTp06ezevVqzbGIxJh6LJIo06dP5/HHH486DBHJQT0WSZR58+ZpjkUk5tRjERGRolJikVhauSHFkEnV9L9iNkMmVbNyQwqAu+++m4cffjji6EQkl9gPhWXUvH8KaAI+C+wGXAZUopr3sVKKCpHNa9nvB+zHtHYW4lIde5HyiX1iCVUDc4Evu/vXzOxzwH3ufkhXq3lf/Hr2+Uh+BcmOV4Vsj8JeJ9Wrjw/VvC9MUhILwIPAI+Ht9cDHI4wlMuX+4mlPPYo41rMvV09FdTs6F9W8L0xiEou7/wP4R3h3LPBArvadqeZ96WrZ5yu5PZby9FTS8nudstWoL7Yk7eNxS8qqeV+YxCSWNDP7IXAIcEKudp2p5n0patnnq5yvUz417I8//ngaGhqYP39+WWLKV5L2pyjFMabWqOZ9YRL1qzAzGw2cDAwPezDSieRTw/7pp5/mrrvuiiA66UpU874wiemxmFk/4IfA19w9+xiXJJZq2EtcqOZ9YRKTWIBzgT2AWWaWXjY0unAkCrfccgtvv/02119/fdShiEgWiUks7v5j4Mctl2ckGekCXnjhBRoayv1zaxFpj6TMsajmvQDwxBNPcNttt0UdhojkEPsei2rei4gkS1J6LCIA/OpXv2LatGlRhyEiOcS+xyKS6eWXX2bLli1RhyEiOSixSKI8+uijsTubXESa01CYiIgUlRKLJMoNN9zA3XffHXUYIpKDhsIkUV599VXefffdqMMQkRyUWCRRHnzwQc2xiMSchsJERKSolFgkUa677jruuOOOqMMQkRxiPxSWUfP+SYKC5x8FehAU+/oEqnnfpbg7mzdvjjqMdkmXua2p30q/ql4qcyudXuwTS6ga+F/geXd/wMyOBK5196FdreZ9VzdjxoyyzrG0r9xy2xUkl61r4IiJcwsLivKVXBbpiKQkFtz9Vxl3Pw28FVUsUl5Df13N0rdbXtE4ueWSiyH/ksvleZ0O3KeS5y7uPLV00r3M5fUN9K9ap15mOyUmsQCY2b4EF53sCRwVcThSJplfWFdffTXr16/n9ttvL9vzt6/XUl65ei4qUtVx6aTS2ATL6xsYfe9CFaFrh0QlFndfCwwML5c/Azg2W9va2lq2bcteaHL79u2x+9mqYsruwpl1rHhnB/Al+Dh5Hq13fm33XKLr2R2wR3cmD+/TbFm2/SluCbCmfiuNTcHtxqbgvuQvMYklnMT/m7tvdPdZZjY9V/u+ffvm3F4cj+YUU3bVGTHEJaZM+cQ0ZFL1zqPgdC31Uh4FJ/V1ioN+Vb2avVf9qnpFHVKiJOnnxicBZwKY2RfQHIskzNSRg+hfVUm3igr6V1UydeSgqEOSLNLvVfoAQO9V+ySmx0Lws+J7zewUYFfg+xHHIxG44oorWL9+fSKvF7Z/754ap0+I9HuVlB5W3CQmsbj7BuCEqOOQaG3YsCFx57GIdDVJGQpTzXsB4K677uLaa6+NOgwRySH2PRbVvBcRSZbYJxaRTOPGjWPDhg3cc889UYciIlkkZShMBIBt27bx3nvvRR2GiOSgHoskym233RaLEzZFJDv1WEREpKjUY5FEGTt2LBs3buS+++6LOhQRyUI9FhERKSr1WCRRbr75Zs2xiMSceiwiIlJU6rFIovzwhz9k06ZNPPDAA1GHIiJZKLFIovTo0YNUKtVsmWrKi8RL7BNLWIflYWCmu59rZvsArwMnA70Irnr8hrufGl2UUoh2VWj8yFGwT/ZCX4XWlFcteZHCxT6xhKrd/dzw9kTCsnhhwa+twJjIIkuw1mvJt6br1JfPv5Z8awp7nTpb3fgkU837wiQlsQBgZkcD7wKvRR1LZ5DPl1g56lHEqaZ8R3ssqtvRuajmfWESk1jMbFfgKuCbwC1ttVfN+7Z9UEe+Leqx5Kfw16m1OvEdlaR9PG5JWTXvC5OYxAKMB+5w981m1mZj1bxvW3Uez5WE16ncteTziSkOFFPHqeZ9YZJ0HsvXgUvM7BXgG8DtZvb5iGOSGFAteSk21bwvTGJ6LO7+lfRtM5sOTHf3v4e/GpMu4pxzzuGdd97h97///c5lqiUvxaaa94VJTGIRAejTpw89evSIOgwRySGRicXdR0Udg0TjZz/7WewmpEWkuaTMsRxpZlNaLjSzbwA3RxCPiIhkEfsei7vPA6qyPDYL6OhvQyWBzjzzTDZv3syTTz4ZdSgikkXsE4tIJjOjvj4eJ1OKSOuUWCRRrrrqKs2xiMRcUuZYREQkIdRjkUQ544wzePfdd5k9e3bUoYhIFkoskigHHXQQ69atizoMEclBiUUSZfz48ZpjEYk5zbGIiEhRqcciiXLKKaewZcsWnnvuuahDEZEslFgkUTXjDzvsMN5+++2owxCRHGKfWDJq3u8JrAGWhw/NAV6lC9a8L2XFxdZrxrevgFUp68aPGzdOcywiMRf7xBKqBp4CBrj7uMwHilHzPv/a7+WQ/GqNhVVhzFfz10n14kXiIymJpaTi8oVUrtoPLYe+drzfyMqNqVYrMMatHsVJJ51EQ0MDf/jDH6IORTqx9GdkeX0D/avWxXp4OI6SlliONLM5QDdgnLv/OeqAyqVUw1/L1jXvqTU2Bcua9zja7kWVcvgr0zHHHKM5Fim5dFJpbILl9Q2Mvnehism1Q5ISyyvAJnefaWaHAvcBA7I1rq2tZdu2bVk3tn379lbH6i+cWceKd3YUIdyOSuZQWHmGvwA+A3yGBzKe64A9ujN5eJ8yPHd22fanKCUppjj1igFq6rfS2BTcbmwK7kv+EpNY3P114PXw9gIz621m3bK179u3b87tZRviqY5wBy/XsNOQSdU7j8ZaDn1FFVN7KKb8KKaO61fVq9lnpF9Vr6hDSpTEnCBpZpea2Xnh7X8F6t39/YjDSqSpIwfRv6qSbhUV9K+qZOrIQVGHlLfjjz+e888/P+owpJNLf0bSB15J+ozEQWJ6LMADwAwzOwuoAM6NOJ7E2r93z8SOF5944omsXbs26jCkk0t/RpLSw4qbxCQWd18DHBN1HBKtH/zgB7GbNxCR5pIyFKaa9yIiCRH7Hotq3kumY489lq1bt/Lyyy9HHYqIZBH7xCKS6fTTT2fNmjVRhyEiOSixSKKcd955mmMRibmkzLGIiEhCqMciiTJ48GBSqRQLFiyIOhQRyUKJRRJl1KhRrF69OuowRCQHJRZJlFGjRmmORSTmNMciibJjxw527IjyIqEi0hb1WCRRhgwZojkWkZhTYpFEOffcc1m1alXUYYhIDkoskihnnnlm0eZYWlbSVJVAkeKIfWIxs8HAw8BMYAnwHaAJ+D6wN/Bz4A13PzWqGKV8UqlUswJuxaqsuWxdA0dMnFvgVmrKVklTJM5in1hC1cD1wBTgEOBgYLi7X21mW4ExUQbX1Q39dTVL325ou2Ex/S6el4grXyXN/By4zzqeuziZJRKipJr3hUlKYgE4GZjh7k3AovCfxEA5v7geeughVq1axSWXXLJzWbF6LcUQpx6Laol0nGreFyZJiaUPsK+ZPU9Q6OsSd/9rtsYdrXkfpaTGdOHMOla8U66fAFcCxm9i1CvIFIceywF7dGfy8D6J2p/ilgBV874wSUosuwEV7n6smX0NmAYMzNa4ozXvo5TUmKrLGPPmzZtxdw49tHi9giGTqpvVN+9fVdnuo9OkvnflFseYWqOa94VJ0gmSa4EXAdx9PrB/tOFIFIYPH86YMcWdUkvXN+9WUaH65gKo5n2hktRjeZagzv3DZjYAqIs4HonARRddxFtvvVXUbabrm4ukqeZ9YRKTWNz9FTM7zszmAr2AH0Qdk5TfiBEjYjdvICLNJSaxALj7NVHHINFav349mzZtijoMEckhKXMsR5rZlJYLzewbwM0RxCMROfXUUxk7dmzUYYhIDrHvsbj7PKAqy2OzgHj+7lRK4tJLL6WuTtNrInEW+8QikunEE0/UHItIzCVlKEwEgLVr11JfH4+z7EWkdeqxSKKcccYZqsciEnNKLJIo48ePZ+XKlVGHISI5KLFIohx33HGaYxGJOc2xSKLU1dWxZs2aqMMQkRzUY5FEOeusszTHIhJzSiySKFdeeSUrVqyIOgwRyUGJRYqq1HXkjz32WM2xiMRc7BNLRs37vQhKFEMwN7QvcCNwGfCUu4+LJMBOpvVqjDUd2lYhdeSzVWKsqamhrq5OV5wVibHYJ5ZQtbufmr5jZmcB+7n7VDN7HxgQXWiF+XC9+I59iXc2bVZi/ENpruRz4D6VqhEvqnlfoKQklp3MrDswBjgq6liKIfNLLMraD3GqG5/WWq+lurqaFStWcPbZZ0cUlXQFqnlfmMQlFuBk4Fl3T+VqlISa963XilePJS17r6U3V5ehrny6dnw+4rA/tZSkmOI2tKma94VJYmI5D2jzuulJqHnfskQ7TjgAAAvZSURBVFZ8HGJqqT0xFaN2fFvcneXLlzNs2LCibrdQSX/vyiWOMbVGNe8Lk6gTJM2sJ9DH3f8edSzyYeWoHX/BBRdwzTXXFH27IplU874wSeuxHAwsjToIaV05asdff/311NbWlvQ5RFTzvjBJSyyfBFZHHYRE5/DDD2fPPfeMOgwRySFRicXdHyY4p0W6qMWLF1NTU6OjSJEYS8ocS7aa96OB8RHEIxEZM2YMEyZMiDoMEckh9j2WNmreTwWmljUgidTEiRN58803ow5DRHKIfWIRyTRo0CAqKyujDkNEckjKUJgIAK+++mrsTvoTkebUY5FEGTt2LKlUihEjRkQdiohkocQiiXLzzTdTU6PL3ojEmRKLJMpBBx3EbrvtFnUYIpKD5lgkURYuXMhrr70WdRgikoN6LJIol112GalUitNOOy3qUEQkCyUWSZRbb71VcywiMafEIokyYMAAunXrFnUYIpJD7BNLRs37J4FPA93Df2OBL6Ka913KSy+9RG1tbeyuFbZmyw4umlRNTf1W+lX1Uilb6dJin1hC1cBGYIG7TzazrwA/d/evJ73mvbTPj3/8Y1KpFN/5zneKvu1ilWdetq6BIybOLWgbrZVlFkmKpCQWgPXA3uHtTwDxKtAuRTX019Usfbvhww98+TKALCWLO4/sZZnbq2PzUQfuU8lzF3fdGu8rN6R21r3vX7VOPdB2SlJiuQVYYGanA3sAR0Qcj5RQri+1UhdfKlbPpRgK6bmoSFXHpZNKYxMsr29g9L0LS17ErjNJUmK5DHjI3SeY2QnAr4CTszWura1l27ZtWTe2ffv22F1zSjEFLpxZx4p3drTRqmv8MqzwnktxXqcD9ujO5OF9Ct5Otv0pbgmwpn4rjU3B7cam4L7kL0mJ5TDgyvD2HOCOXI379u2bc2NxPJpTTIHqHM83ePBgUqkUCxYsKGNEbVuyZAkXPb1u51FuulZ6lEe52p86rl9Vr2bvZb+qXlGHlChJOvN+OXBIePsgYFmEsUhEpk2bFttCX1NHDqJ/VSXdKiroX1XJ1JGDog5JOij9XqYPEPRetk+Seiw/B6ab2bcJEuKYiOORCPTr14/33nsv6jBatX/vnhqH7yTS72VSelhxk5jE4u5rgK9HHYdE6/nnn2fFihX6sIvEWFKGwlTzXgCYMGECd955Z9RhiEgOse+xqOa9ZLr//vtZtkzTayJxFvvEIpKpT58+NDS0cuKkiMRGUobCRAB45plnmD9/ftRhiEgOSiySKDfccANTpnxouk1EYkRDYZIoDz74IEuXLo06DBHJQYlFEmXfffdl06ZNUYchIjloKEwS5cknn2Tu3MIuSS8ipaXEIoly0003MX369KjDEJEcNBQmifLII49ojkUk5tRjkUTZa6+92HPPPaMOQ0RyUI9FIpWu1JdvrfjHHnuMt956S9cKE4mx2CcWMxsMPAw8CVQC+wDbgZHAQIKrHr/h7qdGFWNnV66KivnVit8N6M+kPItfqXa8SPnFPrGEqoHFwN7uflpYQfJn7n6+mW0lRpfQz1qrPW9dozJiuRSvdnw+4vfeHbjPui5du76jVPO+MElJLACfBV4Ib/8BuDXCWLIq5EMcx9oPS5Ys4ZevbIlNDfhClaoHE9f3Lm4xJYVq3hcmSYllMXAc8Ej4f6tXPE6LouZ9frXa2xK/o97OpLQ9mHi9d30+/hHuijqIFlTzvmtIUmKZCtxqZi8Cs4AtuRpHUfM+V632fMTxCLMcMQ2ZVJ13rfg417zviu9de8Uxptao5n1hEvNzY3ff7u7nuvtXgWnAuqhjkuJoT6342bNnM3ny5DJGJ12Rat4XJjE9FjM7HjjM3a8GzgIejTgkKZL21Irv2bMnPXr0KHFE0tWp5n1hEpNYgHnAj8zsFWADcEa04UgUZsyYwapVq/RhF4mxxCQWd98GDIs6DonWlClTSKVSXH755VGHIiJZJGWO5Ugz+1B1JzP7BnBzBPFIRObMmaNCXyIxF/sei7vPI8tPi919FsEvxKSL6N69O927d486DBHJIfaJRSTT9OnTWb16teZYRGIsKUNhIkCQWB5//PGowxCRHCqampqijqEkFi1a1Dn/MJFOauDAgRVRPK++KwrT2vvWaROLiIhEQ0NhIiJSVEosIiJSVJ36V2Fm1g24HRgAVADfdfc3Mx4fCNyUscrngG8CQ4HvAqvC5fe7+9RyxBS22QH8MWPRMeH/OdcrcUynApeFd+e6+3gzGwVcBywPl89x958XIZ7rgKOBjwIXuPufMh47jOA9+yjwmLtPaGudYmgjpiOBXwBNwBvAOcARBAXq/h42e83df1TGmGqBOuD9cNF33X1VVK+TmX0K+G1G037AeGBXSrAPFYOZDQBmAr9291iW6Sg1M7seOAroDtzo7g/ns16nTizA2UCju38lLA52bbgMAHdfBAwGMLPdCapUvkKQWG4p0c6UM6bQZncfnLnAzM7JY72SxGRmHwUmAl8AtgIvhx86gIfcfVyR4sDMjgIOCWMZQJDwjshoci9Bol0VxvFfwP5trFPqmO4Cjg6/uH9HcIWIBqC6VJVN84gJ4Hh3b2jnOiWJyd1X8cFnrRtB8b4ngFMp8j5UDGbWC/hPPqgB1eWY2RHAQe5+mJl9AvgbwcFSmzr7UNhRBEccAM8S7thZjAP+090bYxRTMdYreNvuvp1gB2tw9yZgI/DxIj5/q7G4+2JgPzPrCWBm/YCN7l4Xvk9PERwEZF2n1DGFvhR+cQKsp3SvTXtiKtY6pYhpFDAzM+nF0HsEBwirow4kQi8Bp4W3NwO7mlleOaOzJ5ZPAvUA7r4D6BYeLTUTHpEPo/kVk79lZnPNbJaZ/UuZY/qomf3OzP5oZpe2528pVUzuvhnAzL4A9AHSQyhHmtkcM/uDmR1czFhC9cA+WR5bB+zbxjrFkHP77v4OgJl9EjiWIDkDfM7MnjazF81sSBHjaTOm0N3hc99oZhV5rlPqmADOB+7MuF/sfahg7v7P8PqEXVb4GqST/2hgdr4H3p1mKMzMzgXObbH4i3mufgrwh4wXbTbBMMbccG7hNjpwAcwCYhoH/BfwD6DazOaHtwtWyOtkZp8N4zrL3f8RXml6k7vPNLNDgfsI5mkK0fLvrCCYu8j1WK51iqHN7ZvZ3gQ9qIvcfYOZLQMmAA8CBwDzzOxAd3+vTDFdDcwhSL6PAt/K5+8ocUyY2VeAle7+brioFPuQFJGZDQfOA/I+OOo0icXdpwDNrk4YXrhy7/D2rsAOd3+/ldWHkXEE5e6Z5Qln0XyCv+QxufvkjPZzgc8Da/L8W0oSk5l9mmBM/Gx3/0u4rdeB18PbC8yst5l160hcGXb+naEq4O0sj+1LMFTxzxzrFEOumDCzjwPPAFe5+zOwc07hgbDJm2a2FtgPKMoPLtqKyd3vy4jvWVrsQ62tU+qYQsMI5jLTcZZiH5IiMbOvExykDE33zPPR2YfCnib4lRfACQRHcK35ErA4fcfMJplZuofy1czHSh2TmX3GzGaaWXo46rDw+fP9W4oeU2ga8AN3X5gR66Vmdl54+1+B+iJ8ITwNDA+3eTBQkx6ScPe3gO5mtn/42pwQts+6TpG0tf2bgN+EF0UlbHeGmV0T3t6LYEhoFcWTNSYz+5iZVWfMb6T34ahfJ/jwZ60U+5AUQfiDpknAMHff0J51O/WZ9+GXzzSCo7UU8B13f8vMxhMMdb0cttvg7r0z1vsccDewg+Dnmhe4+xvlisnMJgDHhc//pLtfn229csREUFjtVSCzJzcJWAjMIPgpYgVwaYveXkfjuZGg2/1PgrHdgQS/lPt9+EuVWwiGWGa4+6TW1nH31wqNI5+YCOZTNgEvZzR/IPx3P8ER/C7Ade4+uxwxha/TD8NlKeAvBEN0TVG9Tu7++/DxvwPHuPva8P4nKcE+VCj74FSEvgSfw1XACHffGGVc5WRm5wPXAEszFp/t7ivbWrdTJxYRESm/zj4UJiIiZabEIiIiRaXEIiIiRaXEIiIiRaXEIiIiRaXEIiIiRaXEIiIiRaXEIiIiRfV/X5u1+E3dRJ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_varying_intercepts.trace, varnames=['a', 'b_pop', 'soc_a'])\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "/home/petrushev/projects/jup_venv/lib/python3.7/site-packages/theano/tensor/basic.py:6592: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  result[diagonal_slice] = x\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [soc_a, rho2, eta2, b_pop, a]\n",
      "Sampling 2 chains: 100%|██████████| 8000/8000 [01:11<00:00, 111.76draws/s]\n",
      "There were 721 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 446 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_distance_ac:\n",
    "    a = pm.Normal('a', mu=0, sd=1)\n",
    "    b_pop = pm.Normal('b_pop', mu=0, sd=1)\n",
    "\n",
    "    eta2 = pm.Exponential('eta2', lam=1)\n",
    "    rho2 = pm.Exponential('rho2', lam=1)\n",
    "\n",
    "    soc_cov = pm.Deterministic(\n",
    "        'soc_cov',\n",
    "        eta2 * pm.math.exp(-rho2 * np.square(distances.values)  / 12. )\n",
    "        + np.diag(0.01 * np.ones(n_cultures))\n",
    "    )\n",
    "    soc_a = pm.MvNormal('soc_a', mu=0, cov=soc_cov, shape=n_cultures)\n",
    "\n",
    "    tools_i = pm.Deterministic(\n",
    "        'tools_i',\n",
    "        pm.math.exp(\n",
    "            np.log(data.total_tools.mean())\n",
    "            + a\n",
    "            + soc_a[data.index.values]\n",
    "            + b_pop * (data.log_pop.values - data.log_pop.values.mean())\n",
    "        )\n",
    "    )\n",
    "\n",
    "    obs = pm.Poisson('obs', tools_i,\n",
    "                     observed=data.total_tools.values)\n",
    "\n",
    "    m_distance_ac.name = 'distance autocorrelation'\n",
    "    m_distance_ac.trace = pm.sample(2000, tune=2000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/projects/jup_venv/lib/python3.7/site-packages/pymc3/plots/forestplot.py:318: MatplotlibDeprecationWarning: \n",
      "The tick1On function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use Tick.tick1line.set_visible instead.\n",
      "  ticks.tick1On = False\n",
      "/home/petrushev/projects/jup_venv/lib/python3.7/site-packages/pymc3/plots/forestplot.py:319: MatplotlibDeprecationWarning: \n",
      "The tick2On function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use Tick.tick2line.set_visible instead.\n",
      "  ticks.tick2On = False\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEfCAYAAABiR+CGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5yVZbn/8Q8RUYopCZ52pj+trizammi5aSuoRFtSKDQ1w8yzmJl52KJ51pTykLrNQ2hOSoaJFoJmYsFsTAskLW3jpTKcRFRQ5OBCzeD3x/OMLsZ5ZmbNOtzPvdb3/XrxYh3ueeYaWLOudT+H+9tj/fr1iIiIVMr7QhcgIiL1RY1FREQqSo1FREQqSo1FREQqSo1FREQqSo1FREQq6v2hCxARaQRmth6YB7ydPvR+YCbwXXdf0874JuA5d7+kxO9zrLuPL7PcsqixNAgzex9wPTAAWAcc5+5Ppy/egcAr6dDLgYeAycBHgO+7+5/SbUwGTnL3xe1s/xPAtUC/9KE/A6e7+5sl1tkHeMrdtzezicCRwA3AJHefWjRu+/Sx3bqwzW8DA9z99A7GHOTuk0qptQvfd3u6WKM0jCHu/jyAmfUG7gTOBc6sxMbNrCfJ73DQxqJdYY3jq8Bm7v6fwEnA1UXPneXuQ9I/9wGDST5JjQLGAJjZfsDfM5pKT+Bu4Efuvjvw+fSp89qMK+n15u6HuvvaUr6mO8zsA8CpXRyr3xmpiPRD1wPAZzoY9hEzm2pmS8zsD+kHL8zs82Y228yeNrN/mNnQdPw0YNP08f9X5R8hk2YsjeMTwGMA7v53M9s5bQjt6Qe8mP7pn477Pkmjac8wYK67z0i3v97M/htYl35q/yWwErjRzF4FLgX+BSwCjgY2AiYBHyJpaACY2QKSGRbA/mb2/bS2I4FXi8b9Z9ttunvr7oYNmNn/Ar8HhpLMyPYHxgKfNbPrgZNJZkg7Ar2Bc9x9upnNA+4FXjOzb7r7J9PtHQHsDNwK3Ai8RTIj/Hqb73smcCDQE5jq7udn/FtKgzCzzYBvkLweswwD9gJWAXOArwG3Az8DLnf3X5rZt0heex8HjiLZffapatbeGX36ahxPAV8ys55mNhDYgnd3W33XzJrN7E4z60fy5vxxwIDFJC/WO4GzzOxWM2u7a8eAJ4ofcPe1RbvBdgFGu/u9wDXASHffG3iZ5A14NMlsaE/g8Yz617v7vsBZwA/aPNfeNrOsA1alYx8gaZaXJyX7icChwIvuvg8wkndndj2BB9z9QmCxmbV+yhxJ0hS3AE5LtzsT+Gab73s68EV3Hwi8qplPw5qRziZagAXADJLXX5b73f3V9IPS34CPpo/vDkxMb88EdqhOud2jF3fjuJ/kzf9/Sd7I5wFvknz6+YG7DwZmAxcBj5C8Uf4PcBPJp6SnSd6UTwLOabPtXiRvvFnmufurZtaXpAn9xsxmAENIflE+DfwlHducsY3p6d+z020A0ME2O9I6K1oMbNbmud2Br6XbmgR8KN1V1vq9Ae4BDjCzD5LsxniU5BjVhenXfRPYvM127wUeNLPjgdvcfV0nNUp9GpLOJj5P8vt0h7v/EyDd1fW0mT1dNH5V0e11vPt7djgw28ycZPdXj+qX3nXaFdYg3H096QFCM3s/cLC7vwb8oWjYfcCN6Zve0enYC0g+UX0MWOjur5vZJm02/w/ghOIH0jfdjwNrSHYPtVrq7kPajP0p0LoaatYvyPqM2+1usxPFu8na+36Xufsv29QI7/4cvwF+TTIL/H266+8akmNM96e7vT5U/PXufrSZDQAOJnlDGODub5RQs9QRd19uZtcCPyaZ9ZLOyDuV7l6+EfiCuz9uZh8Hnq1Wrd2hGUuDMLN/N7Nb0rtfBR5MH/+1mf17+vh/krxZtn7NNsAn3H06sBz4qJltRPLJqdhDwI5mNiL9uh7AOOCw4kHuvoLkuMuAdNx30+/tJGemAeyT8SPslf79BWBuF7ZZinW8+yHrL6S/6Ga2hZn9sO1gd3+BpLl9g2RWA8nxmufShro/0DrLwcw2NbNz3f0pdz+PZKbYv8Qapf5cCQwys8Elft2mJB/Y5ppZL5K9CK1nVP4TeF87H/5qSo2lcTwJvN/M/kJygLp1d9Z1wM3pLpwRwAVFX3Nu0f0ZJNP3GSSnFb/D3d8ChgPHmdkskl1DhfTr2zoauMXMHga+SNJUbgP2MLM/kOzWavu67AlgZvcBlwAXd2GbpVgKfMDM7iSZiawys0eAqcDDGV9zL8nZc63PX0Oyi+xO4AqS3WEfBnD3lSQnQTxqZn8imeW85+w6aSzuvprkA9gV6Yexrvo7ya7tZ0g+CE0i2X39IMlr+WFgvpkNqmzFXddDeSwiIlJJmrGIiEhFqbGIiEhFqbGIiEhFqbGIiEhFqbGISEObM2eOzmDqpqx/u7q9QFIvFpG4DBw4MNjV43q/qKy6bSwDBw7sfFANzJ07l5122il0Gd0yfvx4li5dynnnndf54ByL+f8AGqb+oEuSZL1fxP5vXwkd/RvMmTOn3cfrtrFI+Y499ljmzp3b+UARkSI6xiIiIhWlGYtkGjJkCK+tfp3+37yClmWvs0P/jbnliN352OYbhS5NRHJMjaVBHXnrLKb7so4H7XEGAK+9nMRxP/vyGva6fHpHX/GOva0/tx75+c4HikjdUWPJiWE/aeaZl9aELqNipvsyth97X1W/xye37MOD3y91YViRzi16pcDRv5jNvGVr2LH/y5qplyjKxmJmn6VNDKy7v9rxV+VbiDfILs1ayqSZi8SotamsWw/zlq3h6F/MZtqp+hDTVVE2Ft6Ngf1zGkT1TZK0w3csWLCAtWvXhqhtA2+88UbmmVUnTF7Mwtf+WeOKaquaM5ftNuvFjSO37XRcR/8HMWiE+vN2Sm/LstdZl17Zsm59cl+6LtbG8grwIzPrDfwb8Mu2A7bffvta19Sujs4Bb87ZL1NbEyZMYMmSJTzUa493Pr29rwfs2L9PVJ/eYr8WQfXX3g79N97gNb9D/41DlxSVWE83vga4Jo2jvTlwLXVr9OjRjBgxgluO2J0d+/ehZ48e7Ni/D7ccsXvo0kSqqvU13/pBSq/50sQ6Y2kbAzszcD11qVAosHbtWnbafKOoZigi5fpY+pqPcbaVBzHPWDaIge1Gzrl0Yvjw4ZxwwgmhyxCRyEQ5Y3H3m9lwF9jkULXUszFjxrBkyZLQZYhIZKJsLFIbhxxySNRnI4lIGLHuCpMaWLlyJatXrw5dhohERjMWyTRy5EgKhQKzZs0KXYqIRESNRTKdfPLJPP/886HLEJHIqLFIplGjRukYi4iUTMdYJNPy5ctZsWJF6DJEJDKasUimgw46SMdYRKRkaiyS6bTTTmPx4sWhyxCRyFS1sZhZEzDJ3adW8/tIdRxwwAFdOsbSml2hlEkRAc1YpAMvvvgiy5Yty1wrqb08l+KUSWWxiDSmWjSW/c3s+0A/4Eh3/2t7g8zsYuCjwDbp36e6++/N7GDgNOCfwBx3/17W2Br8LA3l0EMPfecYS3cSLrOyWJT8KHmnBMny1KKxrHf3fc1sOPAD4MCMcf8C+rv7l83s08DtZvYn4DJgF3dfbWZTzGyf9sYCaiwVNnbsWBYtWgS0n3DZUQKlZisSMyVIlqcWjWV6+vds4MedjG0GcPf/M7NtgU8md711XZGZwC4ZYzcQQ4Jk3m233XZsueWWDB73YMlJl50lR3Y1/bESYv4/gMaoP29L0ytBsjw1mbFk3O7K160HehQ91oMk476j7wHEkSCZd4sXL2bp0qU0jx3W4bgvXdWc64TJmP8PQPWHoATJ8tTiAsm90r+/AHT2sWtPADP7DLAQeAb4pJl92Mx6AHsDj2WMlQo7/PDDGTt2bKfjlDAp9UYJkuWp9oylJ4CZ3UeSTX94J+PXmNk0YCvgFHd/3czGAtNIjqs0u/vDZja07diq/QQN7JxzzmHhws579seUMCl1RgmS5alqY3H3zhpJW4+4+3VttnEXcFdXxkplDR06NOp9+yISRs2vYzGze0gy64utBB6vdS3SsZaWFhYvXqxPbCJSkpo3FncfVYFtXFCBUqQTRx11FIVCgWHDOj54LyJSTFfeS6YLL7ywS8dYRESKqbFIpsGDB+sYi4iUTHksksndmT9/fugyRCQymrFIpuOPP55CocDw4cNDlyIiEVFjkUyXXnopCxYsCF2GiERGjUUyDRo0iL59+4YuQ0Qio8YimZ566ilaWlqiuo6lvdAxEamt3DQWMxsCnOTuB3Vx/EkkS8T0BH7h7v9TxfIa0kknnUShUGDEiBE1+X4dLcPfHcWhY9BSse0WUzyAyHvlorGYWUlnp5nZdsCxwECSM9vczG5399eqUV+juvzyyzs9K6w7AWD1pLN4gMqpTmOsppgD3RT0VZ6gjcXMmkiSIfsB1wB9zOyXwOeAie5+kZl9FvgpydL4q4AjgMXAnu7+drqdtcBGgBpLBe2+++706dOnwzGVeuOo9GylFmo1W9FCiLWnoK/y5GHGstzdj013hX0a+BRJXc8BF5E0nDPd/VEzO41k1ePzSJoMZjYMeNndXyjeqIK+yjd37lzeeuutDR47YfLikkO/6lXtZiuQ5xlLZ6FtCvpqPHloLLOLbv/V3Quwwe6xT7v7o+ntmcA5rYPNbA/gKuA9F1oo6Kt8Y8aMeSfzvlVzTn+WjsLGYv4/ANUfgoK+ypOHK++LPxK/3cnYdxIkzWxn4OfACHdfVKXaGtrVV1/dpaCvPFDYmFSSgr7Kk4cZS2eeMrNB7v4IsA/wmJn1BJqAA909v/sIIrfLLrvQu3fv0GV0icLGpJIU9FWeGBrLKcD1ZrYeWA4cCewLbAfcYGat4/7b3We1vwnpjtmzZzN//nz9YolISYI2Fnf/dtHtGcCMovv90r+fAvZq86UP8t6wMKmwM844g0KhwMEHHxy6FBGJSAwzFgnkuuuuo6VFexpFpDRqLJJpwIAB9OzZM3QZIhKZPJwVJjn1yCOP8Pjjj4cuQ0QioxmLZDr77LMpFAocdthhoUsRkYiosUimm266iXnz5oUuQ0Qio8YimcyMdevWhS5DRCKjYyySqbm5mdmzZ3c+UESkiGYskun888+nUCjwrW99K3QpIhKRXDYWMzvI3Sd18Py/AbcCvYF/AaPbrm4s5fv5z3/Oc889V9Y22kt0VK6FSH3LZWMBxgKZjQW4GLjZ3X9tZicCpwKn16SyBrLDDjvw5ptvdjquq1kqGyY6ZlMqo0jcQgd9bQOMBz5IsrLxMcChwM5mdg9wMMnMZDvgQ8CF7j4VOBl4I93McmCXGpfeEB566CEWLlzI9x54uaYpkdXJOdlwBYGY0w2l+pQgWZ7QM5aLgKvc/Q9m9hXgHHc/3szOdPdRZrYFMNPdDzez7YG7ganuvgYgXeX4O8AFgeqva5dccsl78liyVCMBslIzF61QK6VSgmR5QjeW3YFPmdm5QE/g5TbPrwQ+Z2Yz0/ubtz6RNpXbgRnu/p79K0qQLN95553Hm2++yeBxDwZJjazszGXDGUtnqYd5EvNrCJQg2YhCNxaAQ9x9ScZzh5GsYrwXsAXwl6LnbgXmu/v57X2hEiTLt9NOOzF37lya9yuv/o7SHWsh5v8DUP0hKEGyPKGvY/kLMBLAzPYxs2+kj7fW9RFgnruvB74OfCAd+01gnbv/oMb1NpQHHniAmTNndj6wE0p3lNgoQbI8oWcsFwBNZnYosB74dvr442b2KMmB/HvN7PMkiZELzOwHwFeAD5rZjHT8/7n7iTWsuyGMGzeOQqHAcccdV9Z2lO4osVGCZHlCB329AAxr5/F9i+7uXHR7Qvr3D6tZlyQmTpzIM888E7oMEYlM6BmL5NhWW23FihUrQpchIpEJfYxFcmzKlClMn975BY0iIsXUWCTTlVdeSVNTU+gyRCQy2hUmmSZNmqRjLCJSMs1YJFO/fv3o27dv6DJEJDJqLJLpnnvuYdq0aaHLEJHIqLFIpmuvvZYJEyZ0PlBEpIiOsUimyZMn4+6hyxCRyKix1KlKBGxtuummbLLJJlWqUETqVa4bS7pU/pNAs7vvb2YXA/uQ5LccT7IMzE+B97v7bsEKrbJyl6TvasBWsb2tP/+10XyWLFmiJS1EpCS5biwpT5vK3sBu7v5FMxsAXO/ue6XrjHWUNlkxw37S3M3Aq5bOh+TMdF/GdPoAxrUVD93qnIK4JCQFfZWny43FzIYBl5BkzL8IfAPoB/ycJHv+beBod19gZqcBXwN6AWPby0tJt/lvvLv+Vy/gCHefl1HC3sBkAHd/ysy2MbOa/k93542unEXsqhGe1ZnicK1CocDTTz/NrrvuWtMaREJT0Fd5SpmxnASc7u7/a2ZfB/oDFwK3uvuvzOxA4CIzO59kiftBwA7AWUDWfpitgR+5+wNm9m3gROC0Dsb+rej+MmBLkt1h71FK0NcJkxdXOcgqnhlLu+Fav678jKWWQVuNEJSVZwr6ajylNJZ7gBvNbAJwp7svNrPdgLHp8zOB84HPAY+5+zrgOeDoDra5HDjHzM4iyV6Z08HYt9rc70FGU4HSgr6aq/iirvWy25UM1ZowYQJLlizhzDPPrHCVtRX70ueqv/YU9FWeLl/H4u5NwL7ACuAPZrYLyRt7j3RID2Bdidu9EHjQ3QeT7GbryFKSFMlW/YGXuvh9GkYlQ7Vuvvlm7r777gpWJxIHBX2Vp5RjLOcC17r7DenB888Cs0mazR0kZ2s9BjxOMgvpSXIM5iZ3/2rGZj8CPGdmPYADSY7TZPkdSQ7LDWa2K9Di7mvNrKs/QkOoZKjWtGnTot4FI9JdCvoqTylX3i8EHkpTG7cjOZB+AfBtM2smSX88393nA78EHgbuBa7rYJs3ANcA95Nk2O9pZl9qb6C7zwH+ZmZ/BW4ETi2hdumGXr160atXr9BliEhkujxjcffbgNvaPLyK9hMgrwSu7MI27ydpKq06PJrr7mcCce/wj0hTUxMvvPCCPrGJSElqch2Lmd1Dstur2Ep3H9m1L7ep7r5/O08MJLlAUqqgqamJQqHA2WefHboUEYlITRqLu4/q5tctADLXFEl3j+3RzbKkEzNmzNAxFhEpmVY3FhGRilJjkUzjx4/nrrvuCl2GiERGjUUy3Xnnnfzud78LXYaIRCaGRSglkIceekjHWESkZJqxiIhIRamxSKbrr7+eX/3qV6HLEJHI5H5XmJkNAe4iudL/9vT2Ue4+1cy+QrLMy3PuflC4KuNRSrLklClTWLOmO/kzItLIct9YUs3AZcBVwJ9aH3T3+8zsdZIl/RtO1ZMld07+WYuX0S/OaxERaU8sjQWS1Y1HAbeELqRY11Il48lj6Uy7eS01pGRJqQUlSJYnmsbi7gWAvK1m3NmbXDVWRw2dLClS75QgWZ5oGkupSkmQLFXpiZPxz1hCzVQqkTTZCAmMeaYEycZTt42llATJUpWSOJmXPIfuJkvmpf5yxP4zqP7aU4JkeXS6cYOoZLKkSL1TgmR5opmxpKcWnwF8ChhoZie7+3uyYKR93UmWvOKKK3jppZe4/PLLq1SVSD4pQbI80TQWd78PCHc6UgN69NFHWb16degyRCQysewKG2xmN7d9MJ3FXB2gnoZw9913c80114QuQ0Qik/sZi7vPAPpnPKdZjIhIzsQyY5EAxo0bx/jx40OXISKRyf2MRcJ54oknWLVqVegyRCQyaiySaeLEiVFfmCciYWhXmIiIVJQai2S6+OKLueGGG0KXISKR0a4wyeTurFy5MnQZIhIZNZYGUEq4V7EJEyboGIuIlCz3jaUoQXIqsB74BNCbZHmXPjRYgmTVw73asbet1pL5ItJluW8sqWZgOrCHu+9pZp8GbnP33UIlSHYt4KtV3Mvmhwz3UrCXhKCgr/LE0lgAJgKT0tvLgQ8HrKXLb3aVWsQuRLgXwIfXLOLv142p+fcVCUlBX+WJprG4+1vAW+ndU4A7OhqvoK/KWNXnY8FmK5UI+YLGCMrKMwV9NZ5oGksrM/sOsBuwf0fjFPTV/XCvYvWwbHjsP4Pqrz0FfZUnqutYzOxo4GvAyHQGIx0oN9zrrLPO4qqrrqpSdSL5paCv8kQzYzGzHYDvAHu6e3X2cdWZ7oR7FXvllVd0HYs0JAV9lSeaxgIcA2wG3GdmrY8pQbKKfvazn0W9b19Ewoimsbj72cDZbR8vajIiIpIDsRxjUYJkAKeffjo//vGPQ5chIpHJ/YxFCZLhrF27ljfffDN0GSISmdw3Fgnnpz/9qY6xiEjJYtkVJiIikdCMRTKdcsopvPrqq9x2222hSxGRiGjGIiIiFaUZi2S6+uqrdYxFREqmGYuIiFRU7mcsRUFfU4BtgA8CHyJZ4fgjNFjQVy195zvfYcWKFdxxR/ZC0t1NpxSR+pX7xpJqBv4PeMjd7zCzwcCF7j4sVNBXvTvy1llM32Q4bEKXl80vJZ1yb+uvVEqROhVLY8Hdryi6+1Hg+VC1hFRacmV+1T6VMsnEUSKldIUSJMsTTWMBMLOtSK603wjYO3A5QVT7TbEWSZW1nq1ohVoplRIkyxNVY3H3F4GB6RphE4ChWWOrmSBZimqk/5WeYJkvtZ+tQNsUz0qlU9aCEiRrTwmS5YmmsaQH8f/u7q+6+31m1tTR+GomSJaiGp+WS0mwLMdZZ53F8uXLGT9+fLvPVyKhshZin7Go/tpTgmR5YjrdeAQwGsDMPkuDHmOppcsuu4xTTz018/lyEypF8koJkuWJZsZCclrxL8zsQOADwJjA9TS8chMqRfJKCZLliaaxuPsrwP6h62gkRx55JK+99hq/+c1vQpciIhGJZVeYgr4C2Hbbbdl6661DlyEikcn9jEVBX+FcdNFFUZ+NJCJhxDJjERGRSOR+xiLhjB49mpUrVzJlypTQpYhIRNRYJJOZsWxZda/CF5H6o8Yimc4991wdYxGRkukYi4iIVJRmLJLp0EMPZdWqVdx///2hSxGRiGjGIu+x6JUCX7qqmWkvfpC567Zm0SuF0CWJSERyP2MpSpCc7O7HmNmWwNPA14CNUYJkp7q7FP4me3yd9aDwLhEpSe4bS6rZ3Y9Jb19OugZ6uspxVAmS9RLU1Z5aLoevwC6pJgV9lSeWxgKAme0DrAKeDF1Ld9XyzbAWoV2gmYrUHwV9lSeaxmJmHwDOBb4KXNPZ+LwGfcUe0tWeMMFdXQ/raoSgrDxT0FfjiaaxAGOBG9x9pZl1OjivQV+1CunqjvaCu67db4vc/dKXKvalz1V/7SnoqzwxnRX2ZeBUM/sz8BXgejP7TOCa6oqCu0QSCvoqTzQzFnf/YuvtNJa4yd3/kZ41JhXQNrhrxIgRrFmzhj/+8Y8BqxKpPQV9lSeaxiK1t++++/LSSy+FLkNEIhNlY3H3b4euoRF873vfi/qgsYiEEcsxFiVIiohEIvczFiVIhrPffvuxZs0aZs6cGboUEYlI7huLhHPAAQfw4osvhi5DRCKjxiKZTjzxRB1jEZGSxXKMRUREIqEZi2QaOnQor7/+Oo8++mjoUkQkImoskumQQw5h6dKlocsQkciosUimY489VsdYRKRkuW8sRUFffYGlwLz0qWnAEzRw0FdrZkTLstfZof/GyowQkVzIfWNJNQNTgQHufnrxE7EFfWUpNzvl2ZfXVDzpcciQIRQKBWbNmtXtukSk8cTSWKqq+qmOLVXcdum6nJ+yxxkA3cpaUcKjxEwJkuWJrbEMNrNpQE/gdHf/ayU2Ws03wK6sjlqrpEcoPe1Rq7tKI1KCZHliaix/Bla4+2Qz+zxwGzAga3B3EySrk/CYnxlL99Ieu15/V1Mda6kREhjzTAmSjSeaxuLuTwNPp7dnmdnmZtYza3x3EyQrnfBYjU/87SU9VuPTVL0cY4l91qX6a08JkuWJ5sp7MzvNzI5Nb38KWObu/wpcVhC1Sno85phjOPDAA6uybZE8U4JkeaKZsQB3ABPM7HCgB3BM4HqCaZv0WC2jR4+OeheMSHcpQbI80TQWd18K7Bu6jkZSKBS6dZxKRBpbLLvCFPQVwPDhwznhhBNClyEikcn9jEVBX+GMGTOGJUuWhC5DRCKT+8Yi4RxyyCE6xiIiJYtlV5gEsHLlSlavXh26DBGJjGYskmnkyJF1cR2LiNSWGotkOvnkk3n++edDlyEikVFjkUyjRo3SMRYRKZmOsUim5cuXs2LFitBliEhkNGORTAeMHMXcpSvZ/JDLFCQmIl2W+8ZSlCA5GZgLHAasB8YAW9DACZLdUdIS/XueyWbAv9avr0qQmIjUp9w3llQzcClwM7AbsCsw0t3Pq5cESahF4FhtdG9p/spQwJhUgoK+yhNLYwH4GjDB3dcDc9I/daWab4i1ChPTbEXqgYK+yhNTY9kW2MrMHiJZ3fhUd/9b1uDuBn1VWikhTdUJGautkLOVtlpDxxohKCvPFPTVeGJqLL2BHu4+1Mz2BH4ODMwa3N2gr0orZdntSoeMdVdrkNjr85JJ4b/vMTjqT2uxL32u+mtPQV/liel04xeBhwHcfSbwsbDl1K/WkKM+O+6GfW4PhRxJw1HQV3limrH8niTc6y4zGwAsDlxP3WoNOVq8eDHPPvusDlpKw1HQV3miaSzu/mcz+y8zmw5sDJwYuqZ6d/jhh2utMBEpWTSNBcDdLwhdQyM555xzWLhwYegyRCQysRxjUYJkAEOHDmXQoEGhyxCRyOR+xqIEyXBaWlpYvHix9jGLSEly31gknKOOOopCocCwYcNClyIiEVFjkUwXXnihjrGISMnUWCTT4MGDo77iW0TCiOXgvQTg7syfPz90GSISGc1YJNPxxx9PoVBg+PDhoUsRkYiosUimSy+9lAULFoQuQ0Qik/vGUhT01Y8klwWSXXhbAT8CzgCmuvvpQQrMudZciZZlr5ecAjlo0CD69u1b5QpFpN7kvrGkmosTIs3scGAbd7/FzP4FDAhXWm2Vk6tSSgokwG7bfIhvbb9G17GISEliaSzvMLNeJImRe4euBbqa+thSk1oq7bEX1vLYCz05+ZHqX4Oq5EfJEyVIlie6xkKSJPl7dy+ELgQ6T32s9OqotUiCbE2BnD17NvPnz+fggw+u6vcTyRslSDCroz4AAAeFSURBVJYnxsZyLHBKZ4OqlSDZvZTHuGYsG6ZAbsx//7WyM5bWZMdaaYQExjxTgmTjiaqxmNlGwLbu/o/OxlYrQbLUlMeQeQ6tSZCtKXg79u9T0qeuJ554gpaWFkaNGlXFKqsv9kwN1V97SpAsT2wXSO4KPBO6iFi0puD17NGjWyl4p5xyCuPGjatSdSL5pQTJ8kQ1YwG2Bl4IXUQsWlPwuuvqq6+mpSWu3XgilaAEyfJE1Vjc/S6Sa1qkBnbZZRd69+4dugwRiUwsu8Kygr6OBsYGqKchzJ49myeffDJ0GSISmdzPWDoJ+roFuKWmBTWQM844g0KhoNONRaQkuW8sEs51112nYywiUjI1Fsk0YMAAevbsGboMEYlMLMdYJIBHHnmExx9/PHQZIhIZzVgk09lnn02hUOCwww4LXYqIRESNRTLddNNNzJs3L3QZIhIZNRbJZGasW7cudBkiEhk1FsnU3NzMwoUL273yuJwAMRGpb7lvLEUJklOAjwK90j+nADujBMmKeu+y/Jtz3tiOVzfuaoBY63L8IlLfct9YUs3Aq8Asd7/RzL4I/NDdv9wICZJdCxPLvw2X468+hYdJdynoqzyxNBaA5cAW6e2PANVNu8qRWrw5VjtATLMViYmCvsoTU2O5BphlZocAmwF7dTS4WkFfpSonpKl7oWL5VMvZStsgsUYIysozBX01npgayxnAne5+iZntD1xBElPcrmoFfZWqnGW3Sw0Vq7QhQ4ZQKBSYNWvWBo+XGyBWa7Evfa76a09BX+WJ6cr7/wB+l96eBuwWsJaGcPvtt7cb9FVugJhI3inoqzwxzVjmkTSTOcAuwLNhy6l/2267LWvWvPekgXIDxETyTkFf5YmpsfwQaDKzb5DMtE4KXE/de+CBB1i0aJF+sUSkJNE0FndfCnw5dB2NZNy4cRQKBY477rjQpYhIRGI5xqIEyQAmTpzIFVdcEboMEYlM7mcsSpAMZ6uttmLFihWhyxCRyMQyY5EApkyZwvTpnS/VIiJSTI1FMl155ZU0NTWFLkNEIpP7XWESzqRJk3jmmWdClyEikdGMRTL169ePvn37hi5DRCKjxiKZ7rnnHqZNmxa6DBGJjBqLZLr22muZMGFC6DJEJDK5P8bSJuirD7Al8AZwBDCQ5Ir859z9oFA1xqyjJMjJkyfj7oErFJHY5L6xpJqBp4At3P3gdHXji9z9ODN7HS3vApSfqZKZBHlPx8vdK2tFRIrF0lgAPgH8Ib39R+C6gLV0ybvJjy2hS6mqWmatKBVSakEJkuWJqbE8BfwXMCn9u92r8fPkwe/XdnXUaqdAtkezFalHSpAsT0yN5RbgOjN7GLgPWN3R4DwlSA4e92DdJEG2Vesc+460TY5s1QgJjHmmBMnGE01jcfc3gGMAzGxL4Bsdjc9TgmTz2GGhy+hQR4mQ9ZBHEfvPoPprTwmS5YnmdGMz28/MLkrvHg7cHbKeepKVCDlhwgTuvffewNWJ1J4SJMsTzYwFmAF818z+DLwCHBq2nPqRlQh58803UygUOPPMMwNUJRKOEiTLE01jcfe1wPDQdTSSadOmRb1vX0TCiGVXWFbQ11eAqwPU0xB69epFr169QpchIpHJ/Yylk6Cv+0jOEJMqaGpq4oUXXtCuABEpSSwzFgmgqamJ3/72t6HLEJHI9Fi/fn3oGqpizpw59fmDidSpgQMH9gjxffVeUZ72/t/qtrGIiEgY2hUmIiIVpcYiIiIVlfuzwuqBmQ0myZQ5yt2nhq6nq8zsYmAf4IPA8e7+WOCSSmZmA4DJwE/cPfcrYrdlZpcCewO9gB+5+12BS+oyM9sIaCLJUNqYJOoimqUcYn/tVEJ3X3+asVSZme0InAr8KXQtpTCzvYHd3P2LJKFqVwUuqWRmtjHwP7wbtxAVM9sL2MXd/wMYBvwkcEmlGgE85u6DgQOBKwLX02Wxv3YqoZzXnxpL9S0FRgErQxdSor1JPq3h7k8B26SfQGPyJslqDS+ELqSbHgEOTm+vBD5gZtH8zrr7RHf/cXr3o8DzIespUeyvnUro9utPu8KqzN0LAGYWupRSbQ38rej+MpJdGvPDlFM6d38beDvCf3vgnfrXpHePBu5393UBS+oWM/sLsBURLckU+2unEsp5/amxVJCZHUO6tH+R89399yHqKdNbbe73AHRuegBmNhI4FvhS6Fq6w92/YGa7AhPNbOcYm2Mj687rT42lgtz9ZuA9a5pFaimwRdH9/sBLgWppWGb2ZeA8YJi7vxa6nlKY2W7Ay+6+yN3/mu5G6Qe8HLg06aLuvv6i2V8rNfc7YCRA+mmzJV1hWmrEzDYlOWliuLu/ErqebhgEnALvhPNtAiwPWpF0WTmvP115X2XpCsxnAJ8iOU6x1N3zHSmZMrMfkUx/3waOdvcnA5dUEjMbCFwJbA/8E1gCjHL3V0PW1VVmdhxwAfBM0cPfcvdFYSoqjZn1Bm4FtgV6k5xuHMXp9rG/diqhnNefGouIiFSUdoWJiEhFqbGIiEhFqbGIiEhFqbGIiEhFqbGIiEhFqbGIiEhFqbGIiEhFqbGIiEhF/X8a7T0jfrpVowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_distance_ac.trace, varnames=['a', 'b_pop', 'rho2', 'eta2', 'soc_a'])\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAGBCAYAAADSejBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd7xcVb3+8c8kAgkkEEjoLQLyyKUTQK4UsYKgIiioCAoCUlRAlAsoYlRUmsgFJCBKFbFwFfEnUpRqoQVBinyBQOiGkgAB0pPfH2sNbMY558zMmXLmnOed135lZu+191r7nGS+s8peq7Ro0SLMzMysMcM6XQAzM7Nu5kBqZmbWDw6kZmZm/eBAamZm1g8OpGZmZv3gQGpmZtYPDqRmZmb98JZOF2AwmTx58lhgB2AqMLuzpTEzq8kIYDxw9YQJE15oRQaTJ09eAxjX4OnPT5gw4fFmlqfZHEibawfgkk4XwsysAZ8Gft7si06ePHmNWXPmPTZyicUavcRrkydPXm8gB1MH0uaaCjDxxzfy+L9fam/OI0e1N7+yYcM7ky8w6bg9O5Lvwadf15F8V33rih3J992brNKRfJ95aU5H8gW458HnOpLvOftu0fY858+dw8znn4T8+dUC40YusRjfOvdGHnumvs/FNVdehm8e8K4lSbVZB9IhYjbA4/9+iQcfb0kLSc+WWtDe/Mo6GEgXDV+iI/k+OG1WR/Kdv2xnfscvzu/M7/i52aWO5Aswdcb8juS72BIjO5Jv1tLuqMf+/TIPPjGjvpNKnfs3UA8HUjMza71Sqf7A6EBqZmZWNgxK9T4o0nN6SeOBs4FtgVeB84FjImJhL+esCjwA/CAiJuZ9NwBbA8Umn4iIjWstpQOpmZm1XokGaqTVd0sqAb8F7gNWA1YErgKmAaf2csXTgWqB9oCIuKC+wr3Bz5GamVnrlYY1tlW3ObARcGhEzIiIB4ATgAN7OkHSTsB6wO+bfWsOpGZm1m02A6ZGxPTCvjuBdSWNrkwsaSRwJnAwUG0k2SckPShphqSrJK1TT2EcSM3MrPXKg43q3aobB0yv2Fd+v3yV9McBN0XEjVWO3U9qIt4WeBswA/ijpJofC3AfqZmZtV6pVP9go54D6aJeznrTMUn/BewDbFgtcUQcUpH+QFJQ3g64tpZiukZqZmZt0EhttMdA+hwwtmLfuMIx4PVBSZOAr0fE87WUMiJeBl4AVqr1zrqqRippKnBCRJxdQ9pFwAcj4qo689geuB4YGRGeL9fMrBl6HzzU8znV3Q6sKWlsRJRnv9kSuD8iXimkW4NUs1xf0kl53yhgoaSPAO8CTgQmRsSzAJLGkoLyI7UWsy2BNAfA5YEVK24SSUcAPwD27c/wYzMzG8CaOCFDRNwl6TbgDEmHAKsDRwPfAZD0ALA/8Pd8rOhU4EngpIiYKWlL4DRJBwPDgbOAu/K5NWlnjfQV4GPAhRX79wKebWM5zMys3ZrbRwqwO3AO8DQwEzgrIiblYwJGRcQCUtB8naTXgJcj4t95127AacCjpBG9NwAf7m1ih0rtDKRXAp+hEEglrQcsB/wrvx8GfDenGwMEcHhE3FR5sdz2/U3g88BSwD9z2slV0o4BbgUujYiJlc2+kg4Cjo6I8VXOnQD8L+mZpVnAb0jPLs1r7MdgZmb9FRFPAjv3cKzHCBwR+1S8f5wUTBvWzsFGVwDvlLRaYd9ewK8r3h9AatNemjRzxWWSqs2avT+wB6mNe3ngV8BVkpYqJpL0lpzHLeUpoer0K1IVfznSQ8AfAg5q4DpmZkNXcx9/GVDaGUhnAH8kBctyjXJP3rx+5yXA2hExJVfJf0kKkmtUud7BwA8j4qGImBsRZwAvAztVpDuN1O69f4Pl3pg04mt+RDwB3EwKqGZmVqvmzmw0oLR71O5FpKbbE4BtgNdyp3H5+LLADyW9hzS0ufx1pNqDsWsDP5J0RmHfcN4cdD8P7Aqs24+m2J2Ar0laG1gsb7/u/RQzM3uzBvpIe378ZUBpdyC9EvippM1Iq7FfUnH8FFINcLuImCJpLWBKD9daCOwZEZf1kt9GwE2kwP2xXtJV/W3laaIuAb4CnBMRcyRdBCzey7XMzKzSsFLa6j2nC7S13hwRc0nNtXuQOnd/XpFkC+C8iCgHz96WsZlCxUwVeVmdosNIzcfvkbR3Yf9s3hwMK4dHl20KvBoRp+cgWiIFZzMzq8cgbtrtRCkvIs3Q/2BETK049gSwtaTF8mjZ/fL+Vatc5xzgYElbSRouaQ/g/orBTAsi4ingUOD0wrEHSfMqkvft2kNZnwBGS9pM0pKk/tb5wCo5qJqZ2RDX9kAaEbeQpnCqbNYFOIa0zM2LpNkmDgQuJ43c3bQi7U9Ji7r+hvSM6jHArnlIdGWeF5NmKzovB8AjgI9KeiRf5wyqNHPnsp6Zz30wb18m1YQvruvGzcyGsvJ6pHVtnS50bdrSR1r5fGZErFvxfvvC2w0qTi/WFl+/Tn5Y9ri8VcuzVPG++JzQn0mz/BedldPdQOHXFxGHkZqIi5atlqeZmfWkkaba7mja7aq5ds3MrEs1cYrAgcaB1MzMWq/5UwQOGA6kZmbWeq6RmpmZ9UNzl1EbULqjlGZmZgOUa6RmZtYGjUxC76ZdMzOzxIONrC4jR8FSC9qb56svtje/shGjOpMv8Nr8+R3Leyj598tzOpLvnPlt/j9UsPji1VZubL3Z89p/zwva9XP2YCMzM7N+GMSDjRxIzcys9QZx0253hHszM7MByjVSMzNrA4/aNTMza5z7SM3MzPqhvIxaved0AQdSMzNrPddIzczM+sHPkZqZmTWuRIlSnYGx1CVtuw6kZmbWdSSNB84GtgVeBc4HjomIhb2csyrwAPCDiJiY940ATgP2ABYDrgYOiojnay1LdzRAZ5LOlXRRfj1R0i0tyGNNSbMlrdvsa5uZDVWlUqmhrRpJJeC3wPPAasB2wCeAw/soxulAZaA9EdgGmAC8FRgJnFfPvQ2oGqmka0g/EEhlGwbMLSaJiMdaWYZ8/RGtzMPMbMgpUf8o3J7Tbw5sBLw3ImYAMySdAHwZOLXaCZJ2AtYDfl/Y9xZgH+CzEfFo3vcV4F+SVomIp2sp5oAKpBHxgfJrSROBHSNiq86VyMzMmqJE3X2kvQTSzYCpETG9sO9OYF1JoyNiZjGxpJHAmcC+eStbG1gamFzeEREPSHqNVEPtvkDaF0kXACMi4pNVjh0CfA3YKiKelPRR4FvAWsDjwEkRcWFO+1fg5nxsJ+BJ4EsRcW1ud38UWC//QNcCzgK2AhYA1wAHR0SHllsxM+s+vTXV9nZOD8YB0yv2ld8vD8ysOHYccFNE3CipGEjHVZxbNiNfpyZd1UfaE0nvBY4Hds5BdCPgV6Qf3ljgMOBsSTvkU+YBBwDnAMsBlwL/J2l0lcv/hPRDXRlYN2/faOHtmJkNOuVRu3VtPVdJF/WS1ZuOSfovUvPtV/tznd50fSCV9DZSINwzIu7Ouz8HXBcRv4uIuRHxJ+BK4FOFU/8eEX+OiLnAKaQO5q2rZLEz8LmImBURL5BGdG3eqvsxM7M+PUeqJBWNKxwDXh+UNAn4eg+jcMtpx1acMxZ4ttbCdFXTbhVLkzqOr4mIqwr73wo8VJH2EWDLwvvXj0fEq5KmA6uQhkYXvRM4XtJ6pEFIbwHuaE7xzcyGhiY37d4OrClpbK7gQPp8vz8iXimkW4M0gHV9SSflfaOAhZI+AmxBatadQOoCBNgQWJw6Pue7vUa6HnAr8HFJm9WQvlhVr/wNlYDZxR2SlgZ+Q+pPXTUiRgDfa7y4ZmZDVKnBrYqIuAu4DThD0hhJGwJHkwYUIekBSduQxr+sDmxS2K4gPX+6U0QsAH5Mqiy9VdIKwA+ByyJiWq231u010jsi4rOSngUukjQhIuYAU4D1K9K+Le8vW6v8QtIoUl/pUxXniFTr/V5EvJr3bdzMGzAzGwqaXCMF2J00zuVp0uCisyJiUj4mYFQOlE8WT8ojcl+OiH/nXd8ERpNG/Q4ntXIeUk85uz2QLsh/H0safftdUofyecCduep+FfA+YEdgh8K520p6P3BTPmc68Ddg1UKap0kP724v6Q/AEaTm33GS3hIR81t1Y2Zmg0pzH38hIp4kjWGpdqzHMyNin4r3c4Ev5q0h3d60C0CuhX4G+JKk7SLiXtKzQj8AXgROAPaKiBsLp10CHEQKoHsAu0fEvIrrPgUcA/yUFFRHAXuR+kpvaOU9mZkNJk0etTugDNgaaZ4HcWLFvn16Oh4Rk4ElCu8vIQXLnrwSER+rku9UCt+DIuIk4KSKZCv3VX4zM3tDC5p2B4xBUSM1MzPrlAFbIzUzs0GkuXPtDihDMpBGxPadLoOZ2VAymJt2h2QgNTOz9nIgNTMz66duCYz1ciA1M7PWcx+pmZlZ49y0a/UZNjxt7TRiVHvzK5v9St9pWuSlOfP6TtQCy4xdpiP5vvTiax3J9/mXZ/edqAVmzurM7xdg2WVHdiTfYR0IHAu7JFgNZA6kZmbWcq6RmpmZ9YMDqZmZWT+U59qt95xu4EBqZmat51G7ZmZm/dDkZdQGEk9ab2Zm1g+ukZqZWct5sJGZmVk/OJCamZn1hwcbmZmZNc6Pv5iZmfWDm3YzSYsDxwKfBFYB5gJ3A9+KiBuaUSBJRwCnR8T8Bs8fAxwH7AqsBLwI3Ax8OyLubUYZzcysTn785XXfB3YBdgOWBt4K3AhcJWnN/hZG0vLAKTRYU5Y0Gvgr8F/ATsCSwATgGeAWSRtWOce1cjMza1i9QeT9wM8LNbuXgImSHgbmAUgaC5wNvAcYDvwNODAinshBax7wMeCrwMbAQ8DewLPAE6TvIC9KOigiLpD0ZeBLwIrA48DXI+I3PZTvKFKA/2hElJeseBo4TNJsUg31Hkk3ALcAOwLTJB0IPAqsFxEP5Ps4AdgqIraXdCypJl62GHBTRLy7zp+fmdmQNJibduutkT4E7CNp0+LOiPhZRDyd354MLA+sBawKLAJOy+nKzbVfBT4DjCM1vX4nIqYBH8jHx+Qguh1wErAHMBo4Abgk11yr2Q04txBEi2U8KiKuLez6JPA54IN93XREHB8RIyJiBKkW/hJwUV/nmZlZVmpw6wL11kgPA34G3CnpCVKz7pXAbyJiTk5zMPCWiHgVQNLvgK9XXOdnEfFwPn45cGAP+f0FWCEiZuS0lwIXABsA11dJvxYQNd7LHRFxZ75ujacA8GPglog4v56TzMyGMo/azSLiSWB7SesD7wXeBfwE+I6kd0XEU6T+yZMlbQKMIjXvvlBxqamF17OBnlbRXSJf+8Okpt3i/mrm5/xq8ViN6V4n6bPAtsD69Z5rZjaUNbtpV9J4UjfitsCrwPnAMRGxsCJdCfgGsB+pFfQx4MSIuDAfvwHYGlhQOC0iYuNay9nQXLsRcV9EnB4RHwPWBhYn1VYBfgNMA9bJTaHVapsLq+yr5ijS4KZdSMF2yT7SP0ztQW5eH8ff9BuUtCqpifrw/IXBzMxqVCo1tlWTg+NvgeeB1YDtgE8Ah1dJfjjwWVLX4dKkpzrOk7RZIc0B5e67vNUcRKGOQCppNUmTJC1b3B8R/wbuAsbmvsvxwMkR8WJOUleBKmwB/F9E3BURi2q41mXA5yUtXaX8F0qq9kMGmJX/Xrywb/WKNOcCf4uIC/outpmZtdDmwEbAoRExIw8SPYHqFbe7gT0jWRARlwEzgPWaVZh6mnafJY3aXU3SV0i1v5Gk5zXfB3w8F+4V4F2S7gH2BLYElpY0KiJe6SOPckB7u6QppFG8W0pakvTc6rdJA31W7eH8H+RyXCnpi6Qf4MqkEbcfBL7Tw3nPkQY9bQv8U9IGwLvJ/a2SPgf8N6lv1szM6tVA026PVVLYDJgaEdML++4E1pU0OiJmlndGxHXl1zmW7EMaBHtd4dxPSPoaaaDsrcAXy+N4alFzjTQi5pICzTPAtaSA+SiwL/DxiLgyj8o9GDia1C+6NSnQPgX8q4Zs/kF6XOavwP6k51aHkQLdz/N1LwROy/2mlWWcBWxDmoDh16R28ztITcJb9vSDyW3qhwBHSXqIVPWfxBtfNPYmNQlMkTQ7b7UOajIzG/Ka2bRL6uucXrGv/L7qUx2SziXFhCOBj0TEM/nQ/cB9pPj2NlKF8I+SehqL8x/qHWz0DPD5PtL8jDSyt+htheOlivRnkzqMy8F664pzt6p4fxhv9MdWy/8V4Ji89ZRm+yr7LgUurdj97XzMz4uamfVDifqfC+0l9aJeTqt6LCIOyN17HydNIrR9RPwjIg4ppsvzCkwn9bteW+VS/8ELe5uZWcs1uUb6HDC2Yt+4wrGqIuLVPFr3ZtIo3mppXia1qK5U6705kJqZWcuVhpUYVudWGtZjJL0dWDPPpFe2JXB/5VgcSb+rMtB0OLBA0mhJZ0laoZB+LCkoP1LrvTmQmplZyzWzRhoRdwG3AWdIGpPnUT8aOBNA0gOStsnJ/wocKWkTScPz+Jr3AVfkQUlbksbdLCNpOeAs0pMof6/13jxhu5mZdaPdgXNI86nPBM6KiEn5mEgTAkFaCGUk6bnTFUgTMuwfEX/Ox3cjzRHwKGlSnxuAD1dO7NAbB1IzM2u5Zk8RmGfa27mHY6XC64XAt/JWLe3jpGDaMAdSMzNruT4GD/V4TjdwIDUzs5YbzMuoOZCamVnrNXdmowHFgbQFJh23J4uG1zwpRlO8Nn9+34la4KU5fc393zp77/u9juR779UndyTfL1z2z47ke+vtdS+U1PXuPmWXjuS7zJKLtT3PWa8t4pVprc9nMDft+vEXMzOzfnCN1MzMWq7JUwQOKA6kZmbWcoO5adeB1MzMWs6jds3MzPqjgRppt7TtOpCamVnLuUZqZmbWD2mwUf3ndAM//mJmZtYPrpGamVnLuWnXzMysHwbz4y+DumlX0l6SXpa0SpVjN0q6tBPlMjMbcnKNtJ6tWyLpoA6kEfEz4FbSoq2vk7Q3sCFweCfKZWY21JTjYr1bNxgKTbsHA/dIen9EXCtpGeBk4H8iYpqkcaRV1rcHhgN/AQ6MiKckrQfcT1qJ/TvAqsD1wD4RMaMD92Jm1pUG8xSBg7pGChARDwPHAz+StATwXeBB4Kc5yanAGOCtwGqkLxen5mPlpU32B7YF1gHGkwKxmZnVaDDXSAd9IM1OAuYCFwCfI9U4F+VjBwAfiYiXI+IV4Apg84rzz46I5yPiWeDHwI7tKbaZmQ10Q6Fpl4iYJ+lAUrPt8RHxr8LhDYCTJW0MLEVq3q1cne+hwuvHgZVbWV4zs8EmDSCq/5xuMFRqpETEX/PL8t9IKgGXA08Ca0XECOALVU4vVbye06pympkNRvWO2G3kudNOGTKBtAcrkfpFT4qIl/K+jaukW6vweg3gqVYXzMxsUGmkf7Q74ujQaNrtxXTgNWB7Sf8C9gYmAGMkLVVI9yVJt+bXnwd+195impl1txINNO12SSQd0jXSiJgDHAIcC7wAbAnsRuojvbeQ9BfATaS+0oeAb7W3pGZm3W0wj9odUjXSiPiPX0tEXAhcWLF7bQBJ6+T3N0bETzEzM6swpAKpmZl1xmAetetAamZmLdfsSesljQfOJk2W8ypwPnBMRCysSFcCvgHsB4wDHgNOzK2RSBpBmkZ2D2Ax4GrgoIh4vtZyOpD2Is+K1B1ficzMBrBSCYbVvYxa9f05OP4WuI/05MWKwFWk8S2nViQ/HPgs8AHgYWBX4JeS7omIO4ETgW1IA01nkrr6zgM+Ums5h/RgIzMza480126dW8+X2xzYCDg0ImZExAPACcCBVdLeDewZyYKIuAyYAawn6S3APsCxEfForoV+BfhwtVXDeuIaqZmZtVxjfaQ9HtoMmBoR0wv77gTWlTQ6ImaWd0bEdeXXkpYkBc5FwHWkgaVLA5ML6R+Q9Bqphvp0LeV0jdTMzFouNe3Wt/USSMeR5gEoKr9fvtoJks4l9aUeSZpf/Zl8neK5ZTN6uk41DqRmZtZtFtV7LCIOAEYBE4GrJG3ayHWqcSA1M7OWa/Jcu88BYyv2jSscqyoiXs2jdW8mjeItp339Wnkg01jg2VrvzX2kLXDw6dfx4LRZnS5GWywzdpmO5X3v1Z1ZFnaDHY7sSL5HfO/QjuR78DZrdiTfeQsX9p2oRcbvVjnwsz1uPf+Qtue5aN7stuTT5MdfbgfWlDQ2Il7I+7YE7s/LYb5O0u+A6yPitMLu4cAC4BFSs+4E0speABsCiwN31FpO10jNzKzlSg3+qSYi7gJuA86QNEbShsDRwJkAkh6QtE1O/lfgSEmbSBou6cPA+4ArImIBaY3p4yW9VdIKwA+ByyKicjnNHrlGamZmLVciDSCq95xe7A6cQxpZOxM4KyIm5WMi9YcCnAKMJD13ugJpQob9I+LP+fg3gdGkUb/Dgd+T5mCvmQOpmZm1XJMffyEingR27uFYqfB6IWmhkaqLjUTEXOCLeWuIm3bNzMz6wTVSMzNruWbPtTuQOJCamVnLDSuV6u4jrTd9pziQmplZ6zWyULcDqZmZWZKadutd/aXmyYU6yoHUzMxarrz6S73ndAOP2u2FpGMl3Zhf7y1pan69naTZkpboZPnMzLpF6iOtf+sGA65GmoPVCRFxdhvyuoG0oOv8KoeXiojjgeMBIuJi4OL8+iZgRKvLZ2ZmA9+AC6QdcEpEHN3pQpiZDXbdUb+sX9cFUkmfBo4F1iDNzn9SeVooSRcAIyLik/n9CGAW8O6IuKGBvCYCO0bEVpL2IdWUV5K0PXA9MDIi2jPjs5lZF+tjNZcezoE6VjPrmK7qI5U0HrgIOIo0j+J+wJmSNupkuczMrHf1Lupd3rpBt9VIHwPGRcSM/P46SdNIS+D8s3PFMjOz3jReIx34ui2QAnxZ0p7AKqQa9RJ5a9RXJR1esW9KRKzfj2uamVmFbgmM9eq2QPpZ4DDgo8CNEbFQ0uO9pK/l1+bBRmZm1rBuC6RbkFY6vx4gL8K6cuH4LGDpwvvV21g2MzPrgZt2B44ngA9KWo7UnHs28Diwaj7+ILCTpCWB2aTa64JOFNTMzN7QyOAhDzbqnzMknVax7z2kwPluUkB9lLSK+ZbARElPA+cB7wemANOAicDeDNz7NDMbEjzXbhtFxPg+kuxQ8f4m4JTC+50qji9NDyJi+z7KMpEUjMsW5f03MHifLTYza4nB+qE54ALpALY28EKnC2Fm1o0amTvXc+0OIpLOBT4BfLnTZTEz60aDefUXB9IaRMQBwAGdLoeZmQ08DqRmZtZyjT3+0h11UgdSMzNrvVIDz4V2Rxx1IDUzs9bzYCMzM7N+KDVQI+2SOOpAamZmrVeigT7SXtp287KaZwPbAq8C5wPHRMTCKmkPBo4gTSn7KPCNiLg8H7sA2AuYXzhldkSMqbWcDqQtsOpbV2T+skNjZsKXXnytY3l/4bLOrJx3xPcO7Ui+p37t9I7ku+K7duxIvtPuva8j+QJcfMreHcl3/PJLtT3P2bNKPDGj73QDiaQS8FvgPmA1YEXgKtKMdqdWpN0N+D7wQeB20mx3v5T0XxExJSc7Pk/A05CuWtjbzMy6U4kUcOrZeqm/bg5sBBwaETMi4gHgBODAKmlHAkdHxN8jYn5EnA+8DGzVlBvDNVIzM2uDxuba7fHQZsDUiJhe2HcnsK6k0RExs7wzIi4pnihpDGnq2GcLu9+Ta65rAPcCh0fEHbWW0zVSMzNrufLqL/VuPRgHTK/YV36/fE8n5Sbhc0lB98959xTSymEfJgXSW4FrJI2r9d5cIzUzs5YrNbCMWi810t6Whal6TNJiwAXA+sB7yoOSIuI7Fen+B/gUsCsp6PbJNVIzM2u58sxG9W49eA4YW7FvXOHYm0gaCfwBWBPYJiL+3dOFI2IBaanOlWq9NwdSMzNruWE00LTb8+VuB9aUVAymWwL3R8QrxYS5OfcXwBzgfcV+VUklSadKWr+wbzFgLeCReu7NzMysa0TEXcBtwBmSxkjaEDgaOBNA0gOStsnJ9wTeDuwREbMrrrOI1C96hqSVJI0CTiIF3ctrLY/7SM3MrOVaMLPR7sA5wNPATOCsiJiUjwkYlV9/jrSe9AxJxfMvzit77Q+cBvyTFBNvA94bEa/WWk4HUjMza7lSA3Pt9va4TEQ8Cezcw7FS4fV7e8sjIl4E9qmrYBUcSM3MrOXKkyzUe043GPDllDRe0iJJb6/jnGskfafvlA2V5yBJU1txbTOzwarctFvv1g2aXiOVdA2wXeH6w4C5xSQR8Viz8y2KiA+08vpmZlafZjftDiRND6TFICZpIrBjRFSd01DSWyJifrVjjZA0DFiUR2KZmdkAUaKBwUYtKUnztb2PNAfXzUmjrHYBlpQ0AjiRNJPEKOAO4AsR8VDh1C0kXUZ6oPYm4DMR8YKkfYCvAlcChwLrSPoZcEtEHJ2D63eBzwBjgCDNo3hTLs9fSA/qrp/L8yLw1Yj4ZT7+DuDHwDrA34CbW/FzMTOz7tSpPtKtgOtIEwdDCnTb5W114F/An/KDsWX7ATsCbyUtmXNK4dgqwCxgNPBURV57AQfkay9NWnrnMknD8/H5wBeAC0mB9hJgUn5QdzhwGXA1aRaN44DP9+fGzcyGoibPtTugdCqQLgR+UmjW3Q84MSKm5md3vklaY27rwjlnR8STEfE8qYZYXCRxDPCDiJhXpVn3EmDtiJiSp376JWlS4zUKaf4WEdfm478AlgVWINWcVwO+HxGzI+LvwBX9v30zs6Gl3Edaz9YtfaSdCqRPlAOepGWBZYDXm3Hzcz3TSbXPsmIz7xPASrnZFmBGRLzcQ17LAmdKekrSbOCevH+JQpqphdflmS9GkoLoyxFRXPb24T7uzczMKgzmUbudCqTzakxXrF0uqDg2t3C8t+udAmwMbBcRI4D1qqRZ2MO5S1TZN7zKPjMz64Wbdlso1/ZmAOuW90lajtQnOaWQdO3C6zWBZ2ocnbsFcF5ElK+1cR3FexpYWtIyhX3r9pTYzMx6VqrzT7cYKDMb/QQ4Ko+gfRH4PimI/o00+AjgIEl/JdVMDwR+V+O1nwC2lvQjYCNSfyzAqsADfZx7K6mJ+ShJ3yIF5Q+SBv0LhJAAACAASURBVCiZmZl1vkaafZO0LM4/Sf2Vq5CWu1kAlEfungvcQFra5nHSCNpaHENqzn2R9IjNgaRZ/S+TtGlvJ0bELOCjwEdIteaJpKbixXo5zczMKjR5GbUBpaU10oiYSAo+fe2bRXpE5YAq13iIN57L/VWV4xeQVj0v7tu+8PofwAYVp+1aeL198UBEPFDIj4i4ucr5p1WWw8zMetZIn2e39JEOlKZdMzMbzBp5nKVLhu06kJqZWcu5RmpmZtYPLVjYe8BwIDUzs5Yr0cDqL13yCEy3DIoyMzMbkFwjNTOzlmvkcRb3kZqZmWWlUv3ri7qPdAh79yar8OL89k7J+++X57Q1v7LnX57dd6IWufX2xzqS78HbrNmRfFd81459J2qBaTde1ZF8WavX+VJaap3lRnck35mzap2GvHnmzWnPRG0lSgyrM5R2Sx+pA6mZmbWca6RmZmb9MJj7SD1q18zMrB9cIzUzs5YbzM+ROpCamVnLNdRH2pKSNJ8DqZmZtVzqI60vNPbW9yhpPHA2sC3wKnA+cExELKyS9mDgCGBl4FHgGxFxeT42DPgOaa3qpUnrYB8QEY82o5xmZmZNUZ5rt96tGkkl4LfA88BqwHbAJ4DDq6TdDfg+8BlgDHAq8EtJa+ckhwH7AjsAK5HWvL4851ETB1IzM2u5Enlx7zq2XiLZ5sBGwKERMSOvI30CcGCVtCOBoyPi7xExPyLOB14GtsrH9wdOioi7I+Jl4CvAeoXjfXIgNTOzlivl9Ujr3XqwGTA1IqYX9t0JrCvpTbNpRMQlEXF2+b2kMaQm3GcljSAFzcmF9DOBh0nBuiYOpGZm1m3GAdMr9pXfL9/TSbm59lxS0P0zMJZU8a12rR6vU2nQBFJJiyS1ZA41SVMlHdSKa5uZDQWlBrceLOolq6rHJC0G/AxYH9g1D0qq+zrVdM2oXUlvByYC7wZGA9OAK4BvVVTvzcxsgBlWqn+u3V7SP0eqTRaNKxx7E0kjgd8BSwLbFGLGC8DCHq71bO3l7AKSNgFuA54gdTCPAj5M+mbx1/xDMjOzAaxJtVGA24E1JRUD4JbA/RHxSjFhbs79BTAHeF+x4hURc4B7gAmF9MsB65BiTk26pUZ6JnB1RBxZ2HevpI8CZwCr5H0rSbqe9AN9CNg7Iu6RtA9wQkSsVD5Z0i3AVRExUdIFwIiI+GQ+NgKYBbw7Im4oFkTS4qS29SkRsU/T79TMbBAqUf8k9KUeGlcj4i5JtwFnSDoEWB04mvQ8KJIeAPaPiL8AewJvBzaJiGrLVU0CJkq6DpgKnA7cERG311rOAR9IJa0AbA28q/JY/uaxb04H8HnSs0JPAX8k/VA/2uQinUP6ZnNAk69rZjZo9TEKt/o5lHrrqdyd9Hn8NDATOCsiJuVjIrVcAnwOWBuYkeNE2cURcUBEnCNpZeAaUrfh9cDH6ynngA+kwFr576gh7cUR8TCApN8CTR0gJOmrwBbA1hHR/oUDzcwMgIh4Eti5h2Olwuv31nCtiaQxOA3phkBaXnW2lpWypxZezyE9iNssHwQ+BOwQES818bpmZoNeeZKFes/pBt1QzkdJo6rWryFtzcOV6b0vu9qxdwJXAt+VVEtQNzOzskYmY+iSlb0HfCCNiBdIbdZHVh6TNFLSZElb93GZWcDiFftW6+X46lWuMRH4FOkh3a/1kZ+ZmRU0+TnSAWXAB9LsUGALSedIWlVSSdJGwFXAPNJQ6N48CCwraX0ASZ8jdSoXj0+QtGReCeAwYEHFNRYUBjd9XdKm/b8tM7OhIVUw662VdrrUtemKQBoR95MG+YwkBc1Xgf8DbgbeExFz+zj/H8APgT9Luh9YE7iBN/qIzwPuA6aQpo66FniNKn3IEXEjabj0xZKW6O+9mZkNBfVOWN9In2qndMNgIwDyaNzP9HK8VPH+bNJadeX3R5DWo6t27kvAThW7ly4cH1+R/su1ltvMzHi9j7TOk1pSlGbrloBvZmY2IHVNjdTMzLpXI4OHuqM+6kBqZmZt0MwpAgcaB1IzM2u5YTR19ZcBxYHUzMxar5H5FbojjjqQmplZ65Xyn3rP6QYOpGZm1nKNzPjXLRMyOJC2wDMvzeG52e39FzBnfuVETO0xc9bQWwRn3sKFHcl32r33dSRf1urQJF6P/KMz+QLz5te1ilbTjB65WNvznM08Xmh7roOLA6mZmbWcBxuZmZn1hwcbmZmZNc59pGZmZv2QZjaqd9Rud3AgNTOzlhsGLKozMnbLZPDdUk4zM7MByTVSMzNrg/onZOiWxl0HUjMzazkPNjIzM+sHTxFoZmbWD8MaWJB0WHfEUQdSMzNrB/eRdoSkqcDywIoR8UrFsSOAHwD7RsQFbS+cmZnVrNRAjbS3PlJJ44GzgW2BV4HzgWMi4j8mw5Y0Kqf9NLBeRDxQODYVWAUonndNRHyk1nIO6ECavQJ8DLiwYv9ewLPtL46ZmXWSpBLwW+A+YDVgReAqYBpwakXaVYDrgVt6ueQHIuKGRsvTDYH0SuAzFAKppPWA5YB/5ffDgO/mdGOAAA6PiJvy8b8AfwDWB3YBXgS+GhG/zN9qHqXwLUXSCcBWEbG9pCWBScAHgZHAPcCXImJyi+/bzGzQaKSRtpdzNgc2At4bETOAGflz+8tUBFJSq+b/AHeTYkTTdcOEDFcA75S0WmHfXsCvK94fAGwHLE36pnKZpOH5+HzgC6RgPAa4BJiUv9X05XDSt521ScH7D8C5Dd+NmdkQNKxUamjrwWbA1IiYXth3J7CupNHFhBFxd0T8ro/iHSbpMUnPS/qVpBXrurd6EnfIDOCPpGBZrtLvSQqGZZcAa0fElIhYAPyS9C1kjUKav0XEtfn4L4BlgRVqyH8FYB4wKyLmAd8DJvTvlszMhp5SnVsvxgHTK/aV3y9fZ7H+AdwGbAJsmq/9617PqNANTbsAF5Gabk8AtgFei4i7JJWPLwv8UNJ7gLG88TtYonCNqYXXs/PfI2vI+1RS2/tTkq4CLs+bmZnVo3mDcBc1eOw/RMSuhbczJH0RuE/S2yLioVqu0Q01Ukj9pCtJ2ow06uqSiuOnABsD20XECGC9Ktf4j5FcvXj91x0RjwMbkmrBz5P6S39Rx7XMzIa88uov9f3p0XOkSlPRuMKx/ng0/71SrSd0RSCNiLmk5to9gN2An1ck2QI4LyKm5Pcb13H5WfnvxQv7Vi+/yMOmh0fEnyPiK8CWwB6Slq0jDzMza57bgTUlFYPplsD9lY9K9kbSGpLOkrRYcXf++5Far9MVgTS7CDgQeDAiplYcewLYWtJikiYA++X9q9Zw3edIo3i3BZC0AfDuwvHfAD+QNDqPDn4H8ALwUqM3YmY21JTn2q13qyYi7iL1a54haYykDYGjgTMBJD0gaZsaijUN+AhwgqSRklYFTgMuj4inar23rgmkEXELKehVNusCHENqzn0ROJEUcC8njdzdtI/rLgQOAY6S9BBwHKn5ttx/vD/wVlKwfhE4AvhotYd+zcysunoHGtUw4Gh3YBngaeBPwDkRMSkfEzAKQNKxkmaTHosEuFvSbEnHRsQc0qONGwHPkILzg8A+9dzbgB5sFBHjK96vW/F++8LbDSpOL3YgF9ORnxct9oNeClxacf6387HHgQ/VXmozM/sPDcxs1Fv6iHgS2LmHY8XP9+OB43u5zj3A++ss2ZsM6EBqZmaDQ6nZkXQAcSA1M7OWa2Su3S6Jow6kZmbWHvXGxboeCO2grhlsZGZmNhC5RmpmZq3npl0zM7PG9TlXUQ9ndUPzrgOpmZm1XIneF+ruiQOpmZkZjbXsdgsH0ha458HnmDpjflvzXHzx4X0naoFll61lAZ3WuPuUXTqS7/jdKtcNbo+LT9m7I/mus9zovhO1wLz5H+9IvgDbfezrHcn3jv93QtvzXDhvdt+JmmEQ95F61K6ZmVk/uEZqZmYt1+hgo27gQGpmZq3Xy2ouPVnUHXHUgdTMzFpvEHeROpCamVkbDOJI6kBqZmYt10gfaf19qp3hQGpmZi1XaqCPtJEJHDrBj7+YmZn1g2ukZmbWcoO4i9SB1MzM2qRbImOdBnwglXQNsF1++xZSc/TcYpKIeKztBTMzs7p0y+Cheg34QBoRHyi/ljQR2DEitupciczMrF6DebDRgA+ktZB0CHAosCrwIHBMRFyTj90C/AbYBPgIMB04IiIuy8e3Bs4C1gb+lLcjI2LNfHxv4OvA6sA04ISI+HH77s7MrPsN5j7Srh+1K2lX4HhgP2As8EPgD5KUk8wHvgj8BFgG+CUwKZ87Evg9cA0wLu//OjAvH18HuBD4CjAKOBCYJGn9dtybmdmgUWpw6wJdH0hJAfTnEfHXiJgbET8D7gd2L6T5S0RcFxELSIF0nKSxwH8DY4DvRsTsiLgauLFw3hRgbET8ISIWRcS1wAvAZu24MTMzG/gGQ9PuW4FrK/Y9kveXTS28Li++NxJYDZgRES8Wjt8JbJ5fl4CvSPoUsEp+v0TezMysRoN5ZqPBUCPtyaLC64U9pClVpKtM+zlSs/DngKUiYgTwTNNKaGY2RJQHG9W7dYPBUCOdAqxbse9twG01nPtvYFlJoyLilbyv2Gy7BfCniLgRQNJKwAr9LK+Z2ZDUzLgoaTxwNrAt8CpwPmmg6X9UnCSNymk/DawXEQ8Ujo0FfgTsBCwA/g/4UkTMqrUsgyGQ/gS4UNLPSM2yewPrkPpC+/J3UlPvkZK+B2wPbMMbz6k+AbxH0rKkpuBJed+qzbwBM7NBr4nDdiWVgN8C95G66FYEriI9WXFqRdpVgOuBW3rI5TxgNKkC9hbgCuBE0pMgNen6pt2IuAL4HvBr0kCg/YH3R8QjNZz7MrAH8CngeWAv4EzeaN49i9S/+iRwNXAy6VvNkZIOaOqNmJkNYqUG//Rgc2Aj4NCImJFrmCeQnqyotDzwP8A3Kw9IWhH4MOmRyGkR8RTpyY19JS1e6711VY00IiYCE6vsP5kU5Kqds03F+3t58/ecq4Cr8oheJH0XeCqnnQ68v+KSfyF9WzEzs87YDJiaP6PL7gTWlTQ6ImaWd0bE3cDduSm40iakRyT/WXGdUYCAe2opTFcF0maTNAx4GPh5njVpLVLT8P92slxmZoNNiQZmNur50DjS5DpF5ffLAzOpzTjgpYp+1eJ1atL1Tbv9kX94nwTeC7xImtXol8AZnSyXmdlg0+T5GCqftqj1WEuuM6RrpAARcSvguXvNzFqpuXMEPkeaya5oXOFYrZ4DxkgaXu7eK1zn2VovMqRrpGZm1h5NHmx0O7BmfnSlbEvg/sKjjLW4kxQHN6q4zkukedtr4kBqZmat18hkDD3E0Yi4izRXwBmSxkjaEDia9NQFkh6QtE31s990nReAXwE/lLSipDVJT4GcExHzar01B1IzM2u5FsxZvztpIZKnSeNbzomISfmYSCNvkXSspNlA5GN3S5ot6dj8/qB8jYeBf5CeOf1GPfc25PtIzcys+0TEk8DOPRwrFV4fT1ohrKfrvATs2Z+yOJCamVnLeWFvq8s5+27BYkuMbGues+ct6DtRCwzr4L/0ZZZcrCP53nr+IR3Jd/zyS3Uk35mzau4qaqrRIzvz+wW44/+d0JF8N//Q0W3Pc901xnLecR9pQ06Dd2lvB1IzM2s510jNzMz6YfDWRx1IzcysTbqlhlkvP/5iZmbWD66RmplZy6Wm3fqqpN1SgXUgNTOz1hvEnaQOpGZm1nKDOI46kJqZWev58RczM7N+6GM1lx7P6QYetVsDSedKuii/PlbSjZ0uk5lZV2nBrPUDxYCukUqaCpwQEWdX7D8IODoixrejHBFxQOF1rxMgm5nZ0DKgA6mZmQ0eXVLBrFvXB1JJE4D/Ja1wPgv4DXAosAZphfPlIuIlSSOBF4HTIuKofO63gQ0iYjdJnwaOzec9C5xUXttO0gXAiIj4pKSJwI4RsVUbb9PMrKsN5sFGg6GP9FfA34HlgM2BDwEHRcQU4EmgHPC2AqYAxVXTtwX+JGk8cBFwFGkx2P2AMyVt1I4bMDMb7EoN/ukG3VAjPUPSaRX7hgNP5dcbA3MjYj7whKSbSQEV0krn7wSuBrYDLgEOlbQEsBB4B2l19MeAcRExI593naRpwATgn625LTOzoaNEAzXSlpSk+bohkH6pp8FG+e1OwNckrQ0slrdf52PXA3vl19sB385/bwnMB16IiJBUAr4saU9gFVJNfYm8mZmZ9airm3YlrUOqZZ5HqlGOyO/LrgPeIWkEsBlwG6kZeBtSs+6fc7rPAocBBwCj8nWeaMtNmJlZV+vqQApsCrwaEadHxJxcs3y9XzMiniANHNoH+FdEzCIF0q1JwfRPOekWwPURcX1ELJS0ArByG+/DzGxQKw82qnfrBt3QtNubJ4DRkjYDHgC+T2qyXUVSKSIWkWqlhwG/z+fcClxM6iM9sHCdD0pajtScezbwOLBqu27EzGxwa2TwUHdE0q6ukUbELcCZpL7QB/P2ZWBDUrAkH3s78Jd8zovAc6T+0WdymrOBh0gB9VrgB8AkUr/pwW25GTOzQcw10g7paeaiPPjo7Pz6MFKNs2jZQtpLgUsrzl+/4v2LwA4V17gJOCW/nlRIOxGYWNsdmJkZePUXMzOz/hnEkdSB1MzMWs6rv5iZmVlVrpGamVnLNXuu3Ty169mkOQFeBc4HjomIhVXSHgYcAaxAmq3u4Ii4Mx+7gfRI5ILCKRERG9daTgdSMzNri2Y11OY5A34L3AesBqwIXAVMA06tSLsraenLj5AefzwS+IOkdSLi1ZzsgIi4oNHyuGnXzMzao3mLem9Omnzn0IiYEREPACfwxtwARfsDP80T7rxGmip2ESmwNoUDqZmZtVyTV3/ZDJgaEdML++4E1pU0ukrayeU3eaKeu3hjcROAT0h6UNIMSVfl6Wdr5qbd5hoBMH/unLZnvGD+gr4TtcDCDj4xPeu1RR3Jd9G82R3Jd/aszvys582Z35F8ZzOvI/kCLOzQ73jdNca2Pc81Vlqm/HJEK/OZO3dW3X2kc+fO6unQOGB6xb7y++WBmTWkXT6/vh94jTSV7ALgDOCPkjaIiJo+zB1Im2s8wMznn+xwMYaGV6Z1ugTt9cSMvtMMJi90ugAdcN5xTWttbMR44G8tuO7zwGtPPT51yQbPfy1fo6i3b9GVx3pKuwggIg4p7pR0ICnQbkea6a5PDqTNdTXwaWAq0JmvtGZm9RlBCqJXt+LiEyZMeHzy5MnrkWqGjXh+woQJj1fsew6orL6PKxyrJe091TKLiJclvQCsVGsBHUibaMKECS8AP+90OczM6tSKmujrciCsDIb9cTuwpqSxEVFuvNgSuD8iXqmSdgJwEYCk4aR+05/m/tQTgYkR8Ww+PpYUaB+ptTClRYs6089kZmbWKEl/Bx4FDgFWJ9WovxMRkyQ9AOwfEX+RtCPwa9Io3duA44BPAoqI2ZLuIC14cjAwnDS3+jrAFtWeSa3Go3bNzKwb7Q4sAzxNWlv6nIgoLzAiYBRARFxFenb0AlJf6zuBnSKi3P22G6l5+1HScpwl4MO1BlFwjdTMzKxfXCM164A8M8uAMJDKYtaNHEjN2kjSmfD6Q+GdLsvW5bI4mJo1zoHU3lQjGYgfqOUySWrpA+OtJmlt3ngIvLyvI/8HJb0DuFrScdC6YDoQ/z0NRpK+KukrnS7HUOVA2gVa+WEkaVixdlR8PRA+BCWV8of8+4GfS1qmz5PoOUB1MHDtGRFTIuIT+f2+ABGxsENlegz4DrCPpG/ksjQlmJbvp/C7e5ekPft73WaUqYXXL/X2vsV5DwfmARd2qgxDnQcbDWDlD6Ja9zdw/WHlkWmSvgRsSppQ5oqIuLmZeTVYvvIH8R7AL0gznCginqrlvPz6AGBN0pRhF0fE08X7bgdJnwI+GRG75A/05YGHgHsj4p05TdvKVPi5LkEa8n8ocGZEnFo83uC19wVeAS6PiHmSFieNqPxVRJzZpFvoqwzl+1sGGF4xH2sr8iv+P1oGeC0i2jK/oaTPA3cUlgTbDnhfRByX33fs/+9Q4hrpAFXxbf57kr4oaUtoXs2h8J//ZOArpKHh6wKn5v+QHe0/y3l/Evgp8CHSpNTv6u2cYg0739dxwIak58b+Kult7awFSlo3Ii6NiF3yrv+OiGmk6ceWl3QTtL1mWv59bk16Xm4+8D1JR+eyNPQ7z+fsB3wN2FHSiIiYS3qsYPFCmpYp/L/5KOm5wuslfbqijM3Or/z/6JvApcAdkj5aa+tJP/Jem/Qox6PlspD+Xe3V7FYG650D6QBR+Y89/wf4GPA70sPG+wHflPS+wvFGPuy2qHj/DuDjwJYR8T/ALcDGwGmSdijnVa2MrSZpM+BHwN4RcSVpJq7VeytL4UNtPLAYsEEOYp8hLej7F0lrtSNwSZoIfLnw/p3AzZKOiYi7SD/31dsdTHM+HwR+Two2R5EWRf6ipGNzmrr+fRVqPtsDT5C+wOyQDz8PvK183WbdRzW53B8CLgbOAT4XEZdUHG/av+PC/43vAweQ1sK8BfghaUWRMc3Kq0reU0jraM7IX7LfAZwE/Bj4eP7352DaBg6kA8ea8Kb+pTVI/yk+GxF7A98nfds8Kn8I1h3gJB1BWuC2aElSU9wLuQn168DRpAeTT5K0raSdivm1Qx5Y9BiwfURcnnf/GXhH7hMq/5yWrnLuLqTpvf4bGAMQEfeQgtpfgAskjWlDU+olEXFwLtPKEfE34BjgUElHRcTdwEdJwfTGXM6FbfrQey+pqfv3EfFb0hqNk4CDJR2Vy1LzB3BOu1hEzCc94D4N+Fb+MjYDWErSLpJWlrRCK25I0jClKd8OAg6PiPOBKZLeIenbucZYLmvTPvskrUz6eb4vIv4E/Ib0//loYA9JSzUrr0Ke5d/LAqUp7b5B+rzYhBTMfwXs6mDaHg6kA4CkDwCPSHp/4cN9KdKSPtflJpzTgcuBV0nBdBtJb5c0vNYAFxGnRsQOOc+18+5nSR+iKwOnALvnvrLrgRVIqx98q43NjuWa6JXAKhFxj6TynNBPkubILEXEAkmHkVa6X7ziEleTaiNbkPp9yx4FziP1Ua7e4nv4BHnhgtyycIOkzSPiRNIH3VcLwXQXYFVJ/4LWf2HJH6hjSbO/kPN8hjTzy3XAccUP4Fqvm/tENwUWku7pKeBYYP38/mukeU/vVFrz8ftNvB8iYmFEzASeATaRtHm+p+NJNeVPSfpTOW1/8ysYS1rbcrqk3UmB9F2kebdPBA6UdJQk0STl30u+5xdInw/Pk/4vbwqczBvB9BvFc6z5HEgHhttJH/y/KzfdkgLmb0jB9P+Ac3PN9HxSv9avgT9S8ThFLSTtDDwkaeeIuC8ifk3qG51OapaCFLROATYAtmpjTQlgFdKH8cmSNsy1HEjBfU5EzJd0KOlD44jcD4ekFSUtGxGzc03wEuBCSe+G1z9Iyk3Eja5E0SdJ+5P6ysq/mxnAvcD/SpoQESeTag/lYPpPUh/u/bm23ezylB8fGi1pmfxzOAPYStLXy+nyIK5rSQOhdpG0Qh2tHaXcJ3gL6QN8IWlu0xmkWvedpD7Z9wFfBG4mBbl+y7WtrSSdknfdArwduAp4CTg5IrYjNTUvqtaKUatCEzaS9sn/3u4FtiF9cToG2C8P1ptEGiC3G+mLxMON5ltZhvz3upLeLWnFiLiWFDznAt/ijWB6KfB55f5vaw0H0gEgImaQvq2fD1wpaceIeDwijgbWJq2bd1ZOPgO4AvgUqSnp331dv7I2GRF/INXMLi33g5JGtW5EGtQDqf9uZEQ8nGt/Ndd861X5YR0R/w/4AelD4YeSNsyHngeGSbqa9OjGeyPi9nyNbwF/AH4v6Uf5OnsD/4/0BeXIHFB/mu/1phbdy+dJ/WP/XR5JGRHXAf9LmhP0rCrB9H8i4o6I+Fj5Z93MMuVAsxvp5/NnSXvksh0BHCHpa4Xky5B+ZttFxLN1tHYsioiXSE2c7yettLEI2JXUFzuWFGwej4jLI+K7ERHNuL/cYvFeYGdJx0bET4GdSN0Cn42Ia3LST5C+RC1oMJ/iQLa3kX5+P5c0Kjfbr0Jaequ4lOpPgL2AbfLvtl+fuXpjMNXHSa0HPwLulrR3RPwd+B5vBNNNSP8WTyN9ubEW8eMvHZb/cy4s/Ac5B9iXNKnynyT9F2ndvEMj4keSvgCsBhwbEX1+IOjNQ/M3IfWJ3pLzPA3+f3tnHiZVeaXxH4sSxGg0aoi4ozmicU3GJeMGGpKIYiIIKIKIgKgoUdBICIsgioksimhUJOA2GgWjqBMVt9HgGpYZAY+iJhoijmvMiKIi88d7rn3pdNvdtXTT8L3P4yNddau+e6vqfmd7z3sYAHRx9/8Mr/USYCmagvDdXDRYVoRB/5a735h7rFOcXytgiLsvMrOH0AbxM3d/Mo4bDZwSx+4EjEItPGfG87+L528DFgKT3f3TcA4K2lSruYYBKIr/ibv/KTb4Jh6tEPFdjkb1szPd/c9mdgEwFJjk7iVJdVZxXgcj52sqsAP6LE5HEfvpwHj0ufwNGcEfuPviOrz/9u7+Ru7vg9AmPwvoEw/PRuzpQe5+fz6yK/Ca1np91An7ISHzueGEYmYtUSS6F2KmH+EiehWMcNoM/db2Ax6PdVeiLNFbKKX7c2CNu58WrytJi1N8n/cDvd19jpndhsaEXeLuv4vnf4FKM4MzZzOhfEgRaQMjDNpBKFJo5e6nUxGZdnT3JYgBOcXM5qEo5j9qaUTz1PwJKCKZhQgYx6Io+FpgVqw1HtV6zgf2jBRqyVON+fOL/2+JWJ0zLNe4H5HzdWiKw3gz2wVtEBcAZ5nIJd9CKcS+kd5ahdK2x5nZTfE+p6I017HAvDCizUtsRLsCv0VTI/5k6tH8IuqGJ5rZvfFdTkW12imx4U1AhuxsU8tGSWFiL7dHG+poxP4egr73BHo36AAAHgNJREFUk9x9MmqheB5YTN2N6B7AbaZ6NQDu/jTQAeiCMh/Ev59BTlrR9bpwOg82EX2IOuH1yGC3N7PL4tDTgMEoGj68BEa0H3AWikaPj/9vBtyJnM8JiG9wL6oND8ydc6nIbfsBD4UR3RbVfz8BxphZn4hMZ6LRYDVmrBKKR4pI1wGYWLHjUerveHdfmYtMf+juj5tZB+RV3+/uL9fiPfO1nE5E7QYxdCcgqvyFaOOZBPRCDOE/5N6jpBFbNefZCUVkP0EGciRwqrvflDtmCkrLPYvqbUMR7b9tGOG5KPLZBaUUT0VpxZmolvwBSrFdAHQCuke6tVTX0DTe/2DEnvyjx3BhE/lkRqx5bzx2CIqaX3P3ASZW5yNoM768FCn0cFK2RM7TjsB0dx+ee/485FwM8orRU7V+b68QdNgMReE7INGFa3LH/Qy4Pf7r5+6rirys/DlsDCxAwUAHF1kqi0zPAE5CY7WuiMdbuvvHJVh3JLCLu/eJv5uhNPIlwCtAN6QytAewNBzl5qXI7JjIU5+iAdZZ1DkfXec4M5uDHNLlyGmZ4BVDrxPKiBSRNgAq1wSBB9Emvzlwj5ltEpHpDMRKPdrdH3H3K2pjRGGt1pjuqN4524Xl7t4DeAyln/Zw90GImHF2pfcotxHdGegJXOMiCI1BtZ3plmuiR2nHy5CB/ZxgL5vZji7VmjOAJSjCPtfdZ6F0+GIUYR+py/E+qDY6PVJ+JUFEGmcD81A6r3dcX2dkzLu6+72R6iVS0kPjvEERzMbAfaUworHGmthEhyHy0L9HjTh7fiJqB5pqZqfV9b1NPbFTgZYownXgJDM7I3fo/6D66EGoRlowTGSmLIPRxEUwOxJlIOZUikynIYMzJheZflLF29a0ZlX7Y0vEhN041lvt7lk/bmeUQt/cReL7IpzRUhjRzVFkfSpwY6SLewCL3H1cHLYQRfyvA7cmI1p/SIa0ARAb0ZFmdoyZbRw32lzEQt2ECmM6AN2YM8ysVRUGuFrkjm2P0mpdw1vPzuFU1Kc5Kf7ujupj9QIz2wrV6A4lej3jPMaiz2GmmY0ys7EolXaXu78dh72J0rdHxGueib83R8xjUJ3qYXffG9jPK4ge3VD6sujoJHctzdz9I2RMXwa6mNn1KBLrFPXnppEqz1o1FmWOirsvQ5mHJUWeR2ZodjCzdhGFPYocjOaoDeOI7Hh3vxJ9tvMKWG479PkPQ8ZlBHJmepgY1aBWpXno8/97QReVQ9w3hwKDzayni2jXETkheWO6AkX3wxE7uc6pZFubW7CPhaoYcvReQW1pLXIvuR+lsVcBV+YNbWFXu7bD7SJyPYlS5jvFw5sC3zWzHePvf6IsQz8XkzihnpBSu/WE3AaaRYrTgL4o1fhwVrdDxIjfIM+yV6R5W3st2bm5m//LtKxpKsQARLa5KjNIJgGGvoi483Hl9yglzGw7RBI6DBGBnoy07hjEchwZtbXs+DNQtPoZMN7dH6j0fr8DXvIcQcfMnkFtFxcgA7ExIlKtKdd15dZuFsZ6E9TT1xW4OaL9atPkpTqvXLq1C0q1fojIJp3cfb5J8vESlPab6u51Yi1bFeQgExN4JGrfGovINsNRZPYeYrAe4+qVLRqmuv7dqHXmUOBqdx8UdfIH0Xc/EPUPdwV6eg26zLVYcyJKEzdH2YzuKMKejH6bJyDHbgr6bJejumn3Yh2jWP97KMJ9JP6+FbW/HYSyLdnvfynKhOzr7q8Wu25C3ZAMaZlReaOMWmcbd7/JRIY5BtHjHwpj2gKRjbqhnr6ja+NNVzKiA1GNZgtUt5pjYuR2j/f8rbu/amZ3I1HvY6p94xLAzPZG9PsXUR3zPuAyd/8g6sOj0EYw1XMMQ1O/316of/RptDnfjuo/5yKmYvus9mZmW6PIoCmKDA53kX1KKtxdnfHLGdOWwFXAzqhGO93dV5XhPDbynDi6mR2OMhj9UKrx9+jzPtndnw1jOgb4GBgXKea6rNceaOvu03KPdUUR73z0PX6A6u/7APfWthRRi7W/jRzMG9z9URPL+06Uwjw9fisPIIZ3K+Qc/ncB6+S5BZ2RYzAQtbZcjWQmT0T91ZOR0XbEzm4XzvBCRH57tojrzWrcK5Boyo3uPiychsnAE+5+takGfxj6rQ0vldOSUDckQ1pGmFob3kGRSSYacC+q1w2JvzNjehLwmLt/bGpxeT/+rlNKLGpC3VDKtivwNXTTDTFJv2V10D+iiOWEcmzyufNpB/wJbeBXRt1oE3dfmTsmU71ZCkxx9z/nntsURbI7op7ELePvZUAbxJKdX2nN7YDlEaGVhOhRzbV9s3IdKmdMW6EopS2Koq4s5XlE+nSlu0/LrXkpsJW79zfVn+eh2uBmqCXnWZOK1s8RWavW0Vo4eCNRbXegS34ve64LYixn11lnA1bD2t9HUd6WQJ8sO2MarXcXcJu794vHdgPed/d3ilyzMxKQWOLuM8Ow7YnaWxaiuvequPa/A89F6v4MdI91qE0WqRbncSYqz+yDeA0Pon7vTYARuUzSxtkek1D/SDXS8qItat/oYhXkltdQtASsJRowEymQjEBp2KcKMKLfRanhI6L+1RHJlB1mZuNc8nRTUC1lOcGkjJuwHEa0BfLor3S1WWTIpPMy8s3dyJDuBvzS1O9KPPd/7v6ku9/i7l3dvQNK6w1FDOTfWoVgQ/aav+XSueUyov3QlI92ldZeHet+hBR83kepuFK22jRBBKW7Kj3VDHg/orcnEJtzZ+Sg3GFmNyLno1ttjKitPVd0FWq9+jXSYD4t99wsxP4+FuhlZi2tDvX8WuBtRFY6ADmHALjanX4KHG9mv4/HXi7GiJpITc2QAP2ZSORhIxd56wXELt8biXw0d/dZrnaTW03DB0ajkkzBRjTu4wyPorr7IESeOyD+683arTXJiDYgkiEtI1zM20mIhHB8PPwaImrkj+uFUpad0cbQy91fq+n97V9ZhVvGf+/F+36CGIx3AkeZGucvQ0zKHwN9TfJiZbkJY/PdHsjSTV9qosb/Pzez3czsbHd/GKn/bIZSWWvB1mZsLo4a3x4ojXeDSeO18volq4lWYRhuQTXBm6owpl+EMV2J0umDvHRDs5vEpj7M3d81tdIMNon8j3NN8BkIPOPqGwV4Ck1kaQvM92jNqeb9D4gsQH5KzPVmdjEyZqNQv+Z4M+ubc8AWoah0krt/XIxjlvuum5lGsf0VsVUfQ0azR3asSyS+J/BlT2kB6315H8Vnu5qKHtjdgX6ZIxzG9GhEtPpN7m0eRvXKA/MZlbqeh4md+6iZPWySC30RCWVciGrvU1D5Yis0YODrJXZaEgpASu2WCbY22Wci8ii7o5pVd9TO8TbwRt5omtnXXcLbdXn/rkhoYUvUjzjD3Sfljt0a1XEGufut8djYOI9r0eZXUiJO3NxbxLpnuvR8qzquF4pav+fuH5nk1qrd6HOvax6GuAWKru9x974lvITq1v2yLhnGKyM4neTuSysdm69bl4zsZBXEoq+h/uOjkK7rjPgML0fp757u/papPn4zsKK6CD2+r++hXt0L3P3ySAPPQozy7dFYuhNRhHsxitgmoLpoP5Q+frNE13Y0KlF8GxHk5kTK/irEVr3B3f8j97q1ygV1WC//HXVDtdDmyFF6C7X47Ida0WaEc4pJHOSvXoIWsdw1Z7/prYBxiKz1OkrFXwesdjH5M+LVUhfjO6GBkSLSMiBujC9vMHc/DxEVbkLkmXZIAGEO8JCZLTazG8I4fqURNbPBpraG1eGx746i3s2Q0P1jKB3VJ7f+24is82bmfbv7iDif2aU2orl130PyacfG5pC/juy390+0YWV9fh/V8r0/j41nFdps+5fmrKtHGP0pplFdWcR/ILqPbq0qMq3q30WeQ7bpdkLp1OHoO+8J9LcKoYIdgPPDiTsPkcq+Ms3t7s8jVavxZjYIta90cfefIRb0i7FmO5e4w8hYtwfq8S3KiMY5ZLrAv0eR2DJgtpmd6O5/Qw7p/wGnVPqN19mIxuvyw+0vQxHohciQZte9EKVSe+ci01e9BLrIue/zRyizchWwV2SzJqPBBy8CfwZah3ODa/xdMqLrCJIhLTFyN8YhZtbXzC42s63c/ecoFdQbTQI5HtVashTSpTV5tybW5dnAtZHyWo3qYi2BTWNjvxTR8Xub2VQz+75Ji7M18F+RrmsG6tmsTQq5EORSew+iuu1xWcowkD2/H/AXJEZfJ8JTGNNm7v5ZKTa16hB1s6bIAdoHGJFLf35CxdzTK0wj4MqG+G39GxrTNSdqsYOR8eyB6t63oEhqF1RP+1GkR6u6tmwP2NTMNkKpw4EozX4S+n1lvbq/QSILs01TeaagNoyj3X1BoddkZq0txvpFanMocKK7/wr9npuhXup+YUwHoNamzlbgJBcz+6lF/2WkUDuhKUcDkQPxJooEf4yM95I4r47596lrRJpPI4cjuMakAHUHSr9/Bxhrkvp71CWeMg3tE8egfvCyyXYmFIZkSEsEy/WJmijpf0DU+A7A42Y2wN0vQimb/VE99C1XY/6EmrxLMzsq6oKXIKr7tIhMFyPP/dNYfwWSDrsn1p+MUmEHesUUl7IqFsX5Zp/Hdci7nwgMNGmzAjQ3s7ORY/DrMIZ1rjNUivxLSejJ3xtZ+m8UIvj8OzA6F5muRCnQ9pQ5MjZJIl6IUqkrY/3VwDkoauljEs+f6u7HI6GHKo2cVQxM2AOR3Z5CTs18FIl9Fzg0d53Po6htAdKG3sfd/+HuHxZxPdsg4zzSzAxlJnYCPjDpBC9Eke8FwDUmLebdUOR8TiFrm1l/FFlnRrgN8E6kwZu6++tIYOIdxFD+HBnTmYgYWDDi8942/v25qZ1lJCp//AoJPnwPaUmfE8ddFuv/Ag1cKPv9m1A3pBppkTCzNp5jQEbt5AFUj3wgNooVqK9sRhwzAfVBdkMb8Fcqr5gEFXZw94xU0h3Vp15CJIysZ/BmtOn/M/fazYEPvcytINWcd77+NBqlytqgus87iHR1sldqX2loVDrv0xGp6R0U5X2ANvHOqK1nlKtlaTj6HhaU2KBXFvJohljLY9FnuW9mTOK5SSgDMAmxdqv8XeWM6N4oNTwZCde/i2qpfw1H5wokATgt+12ZhiwMQOSmV4q8vnaopWQlisoeRgSyh+LxZ9z9V2bWBgk/bIJqmLtHdFrX9XqjqLuDBynIzI5Etd5zPCdUYVIzehopYeXFQgp2Rk1tUXNQ32n7uJ9fRRmGNUi96E5EzDsc/eYWo3ro21W/a0JDI0WkRcBigLNJdzTDZmj49ANhVP8bpW1nmNnWZtbC1UN6KfCCiyVYkzfzKYpkAXZ095loM9gdGdDtEfljPrAg6q0zTCzdf4QRbVJfRjS3+efTyKOR0R+ANswJKFpap4worFU3uxxFCx8hx+dawFCK824Umd4fqfMuiBFbkhRzbLDZ+awxs/1MfZMHuPtjKJJ/D5hrZt+I41bHec5BwvnV/q5ykdFtwBh3H+Pu97v7M2FEm0bqdjD6rk7LpbOfBs4o1ojGey1Ftfpt0ISWHyFm+0aIPJfNLP0AGfUDACvQiPZCpKFbfW1m7cuoVt+vUmr+BWTY36p0zsU4Sp8h4lArM7srSgPDEWFuOnCTSw1rEkpfD0HO9iZFrJlQZqSItAiY2XdQqu9llJ6cF57zDBQd/hrJmF0Ux88B/u4iEtR2jdMRI7Mdal8YH2vdHBtDL7Shn4Ju/ENQbWxPNEmm3MLzWU24HfAtNLj51UrHlFWerxwI52gmUkf6u5kNRd/nLKQg85Kpaf8wpPF7RqTqir5WU1P/ADS4/d0oFVyLNvSWSKijj6nl52rkEHd06bHWZZ3DkOE9AUVDX+QjX69ghY9Gfb5jEMO7VoSwWqyfCUm0Rpq9r6Fe1KXIURmBotDrULtJd1Tv/azqd/zKtQYih+AelDqe7Rqunj1/EBr6vghFw7NRNLgTGnJeShnHlqj2OhxY7O6nhJPyHGJazzcJSwxFLS//KMRxSKg/pIi0QMQm8BK68XcGhpnZgZHmXYFk/u7NjGjgDZQ+q+0avdHm/W/u7lTUr04zs5Nco8ZmIG+2I7DM3ae5+y/d/TivEGovG2Jj6ErFNJllJpLVprljGpURDewEbBRGtCeKTM9FdeeJJrWd94Hz3b2/VxCfSnGtD6Fo7PYgxPRHEe8haHP9UUQzC1DtbBXwXKTx64IDgL3d/XPXFJO8V93EzFqYZuDORcSb84AWVb1RbRGpTaBCvAJFg1ujlO0pyGk8CGnbdkQ91pcCQws0ot1RFuEHSI7zMeBk0yi57FyeRqP4vkAG/GHkkLaP6L3O95FVCFpkTPk1pvapj909c8B7mdmdaC/+BH2/oJamFuieTkZ0HUeKSItAzqPeFXmwy4Eh7r7EzO5DTdPTkad5KCKrHBxGsab37oUiomvd/Yzc421RdPAdpJl7S2z0Z6INqbPXo8qJiUH6R1QDvtvMJqOWiAuBWyJ11ehg6lnsgDbzx4Gz3f0eU8vFeKRXuxSJwpdDFWontJm/jTbcga4e0eYoAzETzaY90zTR5QLgLK8DC9vMjkEGqjPqiaxKP/gxFAGPtiokEet4TYZkC68BnvaYzxrP7YFS0sej6Hgy+l0tRzXCd939jQLX3RPYOBwPzGwH4HTEgp3pGimXHdsKkfO2RFKeRc0TjfJOJ+AOz6kdmdkwlLYdiQz43xAp7xJkQDcDjvQi2NAJ9YcUkRYAyymJmJizy1DEsC0wyczauXsnxDg8DXnUfVGqrjZGdCBSibkd2NfMzs+ei7rUpYhoNNDUX3cL2lhfB+qNTBTYD3g2jOguKE24DKXjeplZ83wUsq7CKqnDuGQGb0SR0XtIaB9kQG9AG/GxZTKiTd39L2jeZgvUhrJHnNfnyBiNQNKPuyJj36UuRjSwGNXX++WNaKUa74eozg+KwIvBJkikoytwrpk9YmYHRi1/CZKzPNKl1ftLlMrdA6U/CzKiAC4lrAWmvusmLlbutYiBe0o+MgU+cfe33D0byl2szOQ26H7tZRJGwTQcfAgaL3c1ipa3Ra1xP0R7xf7JiDYepIi0jrCvVl5piyLTFcB57r7YzLZAkel7tfHmIw01DaXxXkAkoqP5V895VxT17YoICjfknivXKLSvoY3uvqjhfIhSb4eies9zaErFSDO7DrF0F6GU9FkeU1rWFUQd9ItI63353VY65gj0ffzG3a81jb97wkX4KorBWcX5ZL+tJoSAQkSm96Fsw888RA+iFv8Y6resdbmgijV7oHr+CM+NpIvnBiGyUYdiDFml9zwA1ZlHoTphc1RbH4tEJM5Dv7HlQfx516vpgS3BuWSRaSd0f02q4SWFrnMIijbHIZ5DX5SZWhbPt0AavpOAR70eFLoSSotkSGtANZvr8cCNKP20BSKGnOLut+aM6RtoVNgTdVyvLmmotqiG+pa7n1noNdbx/G5GqcVvIEGJ/3UJ3w9EUdEP47ghiCDyFnCFu/9PfZxfbRCGalNk5F9Es1Cfz57Lf9+mPr8RaKNrito09i2kVlfTOYUR/TFKjbdBfZ0Tga8jfdV/IEPzApLqG4lEBAquoUX9bgAi1syJNVcgRnh/RGIqKbM6nJOJiGD0NnIaB6OI9FxEsLmoPhwvM9seXX8/4Fx3v61M6xyGiIktkWPydKTpv4jItwWKRpd4mifa6JAM6VfAKnrtWgOt3P2VIHT8J2ppyfQ/X0e09kHufn3cnH9CPWGneYw6quPazQgWZQ3GtA3wZrkJPbmN/mTkRCx09/1zzw9EhJTu7r4oDOlTaIrNOvkjM03ZuAV4BQ0PfzYez651I9c8020Q8aQtcLtXEItKyoiOmuWNaDB3M2RU5qE054coMt0dpc13RFFkQQLpVaz9A5Th2Bm1mixBDlDRw6mrWa89Mt7DoiywF0rjjkIZjC6F3DcFnsvOyIjdUOrvtNI6ByInewwiIi6Pxxsdqz1hbSRDWgNiE12MIoJL0E3+GqrfZCzcySj1djliHb6CNsIVpfIua0pD1dfNGBFTO6Sk87K7d4zHj0Kb4D9RU/+xqOfxpXKfUzEIksudKDLNG9O8OH1PJIqfCRKU1IhGhPwN1NM53d1vj8d3RXW8l1EE+k3koC1Ac2RLSuSKCAnEXG1STqMS6x2BJsmc7+5/iMeaAa29DrNSS3xOZVX+isj0ZiS+P8Pd/2XSUULjQzKkNcAKV16xUm8G9ZWGqrRmFp3tj2q9/+vuC019d7NQKipL53ZB8ma7Ahd7iQc8lwuVjSka0pz1U96GCEe7lNpRqSKN/F9oCPx1tjYjfAFSERof6fym7v5yKc+l8vlUVdIoB8zscNQLOxpFafUShTYkzOxQ1LZ2I+ozT4pFjRyJtVsDvHDllZJ71EH4mI4ivyrHkpVhzaxP9Em04f3ezPoHQeenQDszmxuHv4v6Z3s0FiMKEOnLrihtOgzVfjENwt4T2M0L7CXMw9YWLG8an+2JZvYr09SWj4EDTepXq01tF8vQZ5qxdl8phxGN915T1b/LCXd/HNVHzycE8td3BG+iP/rNpZTueoAUkVZCJa+83pRXCjzXsqWhcpHolii9eCEauN0DtfRc4u7XBAvzLhSlf4LEI0rC8KxvRGR6B4oAW6OWhH2iTlqUTnGu3r4b8Jm7/yXWuxRJ1t0eNbQnEHtzlFfMvrwaWOnuQ6tdoJHDCpwn2pixIV7z+opkSAO5ja6V5yTQIopoiWo581EdayZqPWgBHAesRtT9o0rNcGxImLRdTwY+cPfB8dh2KLXcDxnTq02arb1QHXFptW/YCBDG7RE0RuuAEhvRfZCR7o4Ggl+F2Lk/zZyPIBzNRr+zRUhvdShwSGP/bBMS1lckQ5qDWcMor6yriI3/OVQfPtzdF8Xj26GotA9wjbv/usFOsgyIeveKXHq1FEZ0T0QUGu7uU+MzHIhSfLdljkq85vso1ZmN+Rrm7gsLvqCEhISyIhnSHEwi4H9GqkHL0SijYUho/g0zGwu87+4Tg2wzFRnT0fWVzi0nKqW1MwNgwLNIb3Vw1rMYxuYcxCI+DH0u69WPqdjUeSUj+jiaDtMx9/zWiIndExnTi3LPfc3dP0npv4SEdR/JkFaCrUPKK/WJXE20Par3rgbujp7QvVCk/kek2JSlIbcFPnX3dxrqvNdV5IzovsCjqK7eBLjT3SfkjvsWYmJ3QRNJxsTjWQ9rvbBnExISCkcypFXA1iHllfpEsHOnIYP5TdScf4JLp3QvlJq8D7hwfXAeyg0z2x0Nhv4l0k4eS9WiGq1Rircz8Ii7/6IBTjchIaFApPaXKuAanDwEsSe3d/crUavHC0jFaF/Ws88u1F0uRYoyPZAjsQsw28wOdkn8HYKIMuNyzfsJ1aMp6ve9OnpQr6cKoXTXVJDrUPr8YDPbqkHONiEhoSCkiPQrsC4qr5QLYUjnAgcjlvLzqBVkUyRKfxRisrYCtljXFYvWNeRSvV8l97gNkoVMqfKEhEaEZEhrwPqovBJG8yg0pPoy4FUkPNARjZd6HPiDu19sZp2RQc0i0G28iJmUCf8yiOB37j65gU8pISGhCKxX6clyYH1TXjGzvZGh/CFS8tkHaO7uL0R09G3Uu3hXvOQNpDF8KLB7MqLFwyvmYd4NnGdmZzXwKSUkJBSBVOeqBdx9rpnNa+xtCEF+mQuMdfcp8VjlFo8WSI5uP9Qf+wPg2+4+r77Pd32Gu79uZtOBT5Fmc0JCQiNFSu1uIDCzjZCSzgp3H5WTP9wOicwfh/pF70ViCxNRy8v+aOZlo9HObUwoVvAhISGh4ZEM6QYEM7sDWOzuo+PvrkA3lOZ9D2UoHnH3U03j0rZG80SXNdApJyQkJKzzSIZ0A4Bp3uXGqMXim4iR2wY4AdVC73L3e8ysHxo6fERi5SYkJCTUDolstAHA3deEeMQ4FHUeg+aG9gZ+4e73xKEPogk3jV7uMCEhIaG+kMhGGxDc/aVI534ONKmilacLsAqleRMSEhISaoGU2t2AYWYtM2NqZmcDF6G0biIWJSQkJNQSyZBuoDCzTYGLUYr3A9RT2n19mqeakJCQUB9IhnQDRoyC64baXp5KQvQJCQkJdUcypAkJCQkJCUUgsXYTEhISEhKKQDKkCQkJCQkJRSAZ0oSEhISEhCKQDGlCQkJCQkIRSIY0ISEhISGhCCRDmpCQkJCQUASSIU1ISEhISCgCyZAmJCQkJCQUgWRIExISEhISikAypAkJCQkJCUUgGdKEhISEhIQi8P9l/Hd24pICXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 600x400 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(dpi=100)\n",
    "plt.imshow(np.mean(m_distance_ac.trace['soc_cov'], axis=0),\n",
    "           cmap=Blues,)\n",
    "plt.gca().set(xticks=np.linspace(-.3, n_cultures-1.3, n_cultures),\n",
    "              yticks=np.linspace(0, n_cultures-1, n_cultures),\n",
    "              yticklabels=distances.index.values)\n",
    "plt.gca().set_xticklabels(distances.index.values, rotation=45)\n",
    "              \n",
    "plt.grid(None)\n",
    "plt.colorbar()\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/projects/jup_venv/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>var_warn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>distance autocorrelation</th>\n",
       "      <td>66.43</td>\n",
       "      <td>3.58</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.41</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>varying intercepts</th>\n",
       "      <td>69.75</td>\n",
       "      <td>4.87</td>\n",
       "      <td>3.32</td>\n",
       "      <td>0</td>\n",
       "      <td>2.48</td>\n",
       "      <td>1.06</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           WAIC pWAIC dWAIC weight    SE   dSE var_warn\n",
       "distance autocorrelation  66.43  3.58     0      1  2.41     0        1\n",
       "varying intercepts        69.75  4.87  3.32      0  2.48  1.06        1"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.compare({\n",
    "    m_varying_intercepts: m_varying_intercepts.trace,\n",
    "    m_distance_ac: m_distance_ac.trace,\n",
    "})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}