Bachelors_thesis_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bachelors thesis' analyses\n",
    "\n",
    "*This Jupyter notebook is for the analyses and model building for Riku Laine's bachelors thesis*\n",
    "\n",
    "**Contents**\n",
    "\n",
    "1. [Compas]()\n",
    "* [Creation of synthetic data]()\n",
    "* [Implementation of competing algorithm]()\n",
    "* [Implementation of our model]()\n",
    "* [Comparisons and other validation analysis??]()\n",
    "\n",
    "*etc etc...*\n",
    "\n",
    "## COMPAS data\n",
    "\n",
    "*The following data filtering procedure follows the one described in the [ProPublica methodology](https://www.propublica.org/article/how-we-analyzed-the-compas-recidivism-algorithm).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7214, 53)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scs\n",
    "%matplotlib inline\n",
    "\n",
    "# Read file\n",
    "compas_raw = pd.read_csv(\"../data/compas-scores-two-years.csv\")\n",
    "\n",
    "# Check dimensions, number of rows should be 7214\n",
    "compas_raw.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6172, 13)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Select columns\n",
    "compas = compas_raw[['age', 'c_charge_degree', 'race', 'age_cat', 'score_text', 'sex', 'priors_count'\n",
    "                     , 'days_b_screening_arrest', 'decile_score', 'is_recid', 'two_year_recid', 'c_jail_in', 'c_jail_out']]\n",
    "\n",
    "# Subset values\n",
    "compas = compas.query('days_b_screening_arrest <= 30 and \\\n",
    "                      days_b_screening_arrest >= -30 and \\\n",
    "                      is_recid != -1 and \\\n",
    "                      c_charge_degree != \"O\"')\n",
    "\n",
    "# Drop row if score_text is na\n",
    "compas = compas[compas.score_text.notnull()]\n",
    "\n",
    "compas.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                length_of_stay  decile_score\n",
      "length_of_stay        1.000000      0.207478\n",
      "decile_score          0.207478      1.000000\n",
      "[[1.         0.20747808]\n",
      " [0.20747808 1.        ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.2074780847803181, 5.439585463018677e-61)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out = pd.to_datetime(compas.c_jail_out, format=\"%Y-%m-%d %H:%M:%S\")\n",
    "in_ = pd.to_datetime(compas.c_jail_in,  format=\"%Y-%m-%d %H:%M:%S\")\n",
    "\n",
    "compas['length_of_stay'] = (out - in_).astype('timedelta64[D]')\n",
    "\n",
    "# Correlation should be 0.2073297, but R uses n-1 \n",
    "# as denominator in variance. Reference:\n",
    "# https://stackoverflow.com/questions/53404367/why-pearson-correlation-is-different-between-tensorflow-and-scipy\n",
    "print(compas[['length_of_stay', 'decile_score']].corr())\n",
    "\n",
    "print(np.corrcoef(compas.length_of_stay, compas.decile_score))\n",
    "\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "pearsonr(compas.length_of_stay, compas.decile_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25 - 45            3532\n",
      "Less than 25       1347\n",
      "Greater than 45    1293\n",
      "Name: age_cat, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "print(compas.age_cat.value_counts())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "African-American    3175\n",
      "Caucasian           2103\n",
      "Hispanic             509\n",
      "Other                343\n",
      "Asian                 31\n",
      "Native American       11\n",
      "Name: race, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "print(compas.race.value_counts())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Black defendants: 51.44%\n",
      "White defendants: 34.07%\n",
      "Hispanic defendants: 8.25%\n",
      "Asian defendants: 0.50%\n",
      "Native American defendants: 0.18%\n"
     ]
    }
   ],
   "source": [
    "print(\"Black defendants: %.2f%%\" %            (3175 / 6172 * 100))\n",
    "print(\"White defendants: %.2f%%\" %            (2103 / 6172 * 100))\n",
    "print(\"Hispanic defendants: %.2f%%\" %         (509  / 6172 * 100))\n",
    "print(\"Asian defendants: %.2f%%\" %            (31   / 6172 * 100))\n",
    "print(\"Native American defendants: %.2f%%\" %  (11   / 6172 * 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Low       3421\n",
       "Medium    1607\n",
       "High      1144\n",
       "Name: score_text, dtype: int64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compas.score_text.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>race</th>\n",
       "      <th>African-American</th>\n",
       "      <th>Asian</th>\n",
       "      <th>Caucasian</th>\n",
       "      <th>Hispanic</th>\n",
       "      <th>Native American</th>\n",
       "      <th>Other</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</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>Female</th>\n",
       "      <td>549</td>\n",
       "      <td>2</td>\n",
       "      <td>482</td>\n",
       "      <td>82</td>\n",
       "      <td>2</td>\n",
       "      <td>58</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Male</th>\n",
       "      <td>2626</td>\n",
       "      <td>29</td>\n",
       "      <td>1621</td>\n",
       "      <td>427</td>\n",
       "      <td>9</td>\n",
       "      <td>285</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "race    African-American  Asian  Caucasian  Hispanic  Native American  Other\n",
       "sex                                                                         \n",
       "Female               549      2        482        82                2     58\n",
       "Male                2626     29       1621       427                9    285"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab = compas.groupby(['sex', 'race']).size()\n",
    "tab.unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGsCAYAAACCSeB1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHR1JREFUeJzt3X20XWddJ/Dvj4byPhToi9A2RGx4nREopeCgDlBUSh3bUaq8DI2dSsexKiwUjTDjIAud4loKMiBOpUCZQaEg2ko7DggUB5S3QqlC0QamJrHQFmlLS6lS/M0fZ0eP6U1yb3Jv7pObz2ets87ez37Ovr+d5D75nmfvc3Z1dwAAGNNdVrsAAAB2TVgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJa2tEVd2jqv6wqm6uqnfsos9zq+o9+7u21VZVv1VV/2W16wBYblV1a1U9ZLXrYGWV71k78FTVZUkeneRbuvvvprbnJfmpJP+6u+9YxfIWVFWV5PNJbu/uR652PQC7U1XPSfKiJA9PckuSK5L8cnd/aFUL46BkZu0AU1UbknxXkk7yA3ObHpzkr3YV1Kpq3YoXt3vfneTIJA+pqsfvrx9aVYfsr58FrA1V9aIkr07yK0mOSrI+yW8mOXU16+LgJawdeM5I8pEkb06yKUmq6peS/GKSH5mmxM+qqh+tqg9X1auq6itJXja1/eO7wqp6VFW9t6q+UlXXVdVLpvYTq+rPquqmqvpiVb22qg6de11X1Y9X1dVVdWNVvW6aOdudTUkuSnLpjrrn9ndZVb2iqv50qv8Pq+oBVfXWqvpqVX18Cqk7+j98ru6/rKofntv25qp6fVVdWlVfS/KUqe0Vc31Oraorpn1/vqqePrWfWVVXVdUtVfWFqvqPc695clVtr6qfqarrpz+XMxfzFwYcOKrqvklenuSc7n5Xd3+tu7/R3X/Y3S/e3fhYVRum8XHd3P4uq6ofm1t//tw489mqOn5q3zyNRzva/93ca46rqg9Ol7l8uarePretq+q4afmUqvrUNLZtq6qXzfXbUdumqto67eelK/hHyXLqbo8D6JFkS5KfSPK4JN9IctTU/rIk/2uu348muSOzU6PrktxjavvQtP0+Sb6Y5GeS3H1af8K07XFJnji9bkOSq5K8cG7fneTdSQ7L7B3nDUmevpua75nkq0mekeSHknw5yaFz2y+bjuvbktw3yWeT/FWSp001vCXJm6a+90qyLcmZ07bjp/09atr+5iQ3J3lSZm9G7j61vWLafuK0/Xum7Ucnefi07ZSphkryb5LcluT4aduTpz/Plye563QstyW532r/m/Dw8Fi+R5KnT7/r63axfZfj47Te86+dxrcfm5ZPT/I3SR4/jTPHJXnw3LYHTePSjyT5WpIHTtt+N8lL58a075zbfyc5blp+cpJ/NfX79iTXJTltp9p+O7P/Dx6d5O+SPGK1/8w99vwws3YAqarvzOx054XdfXlm14A9Zzcvuba7/3t339HdX99p2/cn+VJ3/1p3397dt3T3R5Okuy/v7o9Mr7smyf/ILLzMO7e7b+rurUk+kOQxu6njBzMbFN6TWchbl1kwmvem7v58d9+c5H8n+Xx3/3HPTuu+I8lj5+q+prvfNNX3ySS/l+SZc/u6qLs/3N3/0N237/Rzzkryxu5+77T9b7r7c9NxXzLV0N39wane75p77TeSvLxn77IvTXJrkoft5riBA88Dkny5d3FJySLHx135sSS/2t0fn8aZLd3919N+39Hd107j0tuTXJ3Zm8tkNvY8OMmDpvF6wevmuvuy7v7zaR9XZhbydq7tl7r769396SSfziy0MThh7cCyKcl7uvvL0/rvZKdTijvZtpttx2YW9u6kqh5aVe+uqi9V1Vczu27j8J26fWlu+bYk955e+5npVOatVbUj6GzKLGDe0bMPRLxrgbqvm1v++gLr956WH5zkCdMpiJuq6qYkz03yLXP99/a4T66qj0ynV2/KbPZs/rj/dqcB/B+PG1gz/jbJ4bu6zneR4+Ou7G78OWO6PGPHuPYv5/b7c5nNxH1sGmP/wy728YSq+kBV3VBVNyf58QVqW3DsZmyrfdE5i1RV90jyw0kOqaodv2x3S3JYVe3qndHuPuq7Lcmzd7Ht9Uk+leTZ3X1LVb0w/3zmape6+1E71X1MkqcmObGqfmhqvmeSu1fV4XPBc7G2Jflgd3/P7srYw+u/befGqrpbZjN0Z2Q2M/eNqvqDzAZI4ODxZ0luT3JakncusH134+PXpucdl34kd34judD48+DMTk+elOTPuvubVXVFpvGnu7+U5PlT3+9M8sdV9SfdvWWnXf1OktcmObm7b6+qV2fxQZKBmVk7cJyW5JtJHpnZKcfHJHlEkv+bWcBYqncn+ZaqemFV3a2q7lNVT5i23SezgebWqnp4kv+0D3U/L7Przx42V/dDk2zPrsPinup+aFU9r6ruOj0eX1WPWOTrz09yZlWdVFV3qaqjp2M8NLPwe0OSO6rq5CTfuxf1AQew6VKMX0zyuqo6raruOY0zJ1fVr2Y342N335DZNWn/vqoOmWbA5sPZG5L8bFU9rmaOm4LavTJ7k3lDMvuwU2Yza5nWT5/e+CbJjVPfby5Q/n2SfGUKaidm95fJcAAR1g4cmzK7rmtrd39pxyOzd1HPzRJnSbv7lswusv+3mU2LX53kKdPmn83sl/yWzN7tvX2hfSyh7t+cr3mq+7ey+1O4u6v7e5M8K8m1U+2vzCxoLeb1H8vswwmvyuyDBh/M7ALfW5L8dJILMxsMn5Pk4qXWBxz4uvvXM/uOtf+cWYDaluQnk/xB9jw+Pj/JizM7nfqoJH86t993JPnlzGbAbpn2d//u/mySX8tsVu+6zD4k8OG5fT4+yUer6tbMxqUXdPf/W6D0n0jy8qq6JbPAeeHe/QkwGl+KCwAwMDNrAAADE9YAAAYmrAEADExYAwAYmLAGADCwIb4U9/DDD+8NGzasdhnACrv88su/3N1HrHYdBwLjIqx9ix0ThwhrGzZsyCc+8YnVLgNYYVX116tdw4HCuAhr32LHRKdBAQAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBLSqsVdVhVfXOqvpcVV1VVd9RVfevqvdW1dXT8/2mvlVVr6mqLVV1ZVUdv7KHAACwdq1bZL/fSPJH3f3Mqjo0yT2TvCTJ+7r73KranGRzkp9PcnKSjdPjCUlePz0DsAI2bL5kxfZ9zbmnrNi+gcXZ48xaVf2LJN+d5Pwk6e6/7+6bkpya5IKp2wVJTpuWT03ylp75SJLDquqBy145AMBBYDGnQR+S5IYkb6qqT1XVG6rqXkmO6u4vJsn0fOTU/+gk2+Zev31qAwBgiRYT1tYlOT7J67v7sUm+ltkpz12pBdr6Tp2qzp6ud7th69atiyoWYC0zLgILWUxY255ke3d/dFp/Z2bh7bodpzen5+vn+h879/pjkly78067+7zu3tjdR6xfv35v6wdYM4yLwEL2GNa6+0tJtlXVw6amk5J8NsnFSTZNbZuSXDQtX5zkjOlToU9McvOO06UAACzNYj8N+lNJ3jp9EvQLSc7MLOhdWFVnJdma5PSp76VJnpFkS5Lbpr4AAOyFRYW17r4iyQkLbDppgb6d5Jx9rAsAgLiDAQDA0IQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAANbVFirqmuq6s+r6oqq+sTUdv+qem9VXT09329qr6p6TVVtqaorq+r4lTwAAIC1bCkza0/p7sd09wnT+uYk7+vujUneN60nyclJNk6Ps5O8frmKBQA42OzLadBTk1wwLV+Q5LS59rf0zEeSHFZVD9yHnwMAcNBabFjrJO+pqsur6uyp7aju/mKSTM9HTu1HJ9k299rtUxsAAEu0bpH9ntTd11bVkUneW1Wf203fWqCt79RpFvpenOSwI444YpFlAKxdxkVgIYuaWevua6fn65P8fpITk1y34/Tm9Hz91H17kmPnXn5MkmsX2Od53b2xu49Yv3793h8BwBphXAQWssewVlX3qqr77FhO8r1J/iLJxUk2Td02JbloWr44yRnTp0KfmOTmHadLAQBYmsWcBj0qye9X1Y7+v9Pdf1RVH09yYVWdlWRrktOn/pcmeUaSLUluS3LmslcNAHCQ2GNY6+4vJHn0Au1/m+SkBdo7yTnLUh0AwEHOHQwAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADGzRYa2qDqmqT1XVu6f1b62qj1bV1VX19qo6dGq/27S+Zdq+YWVKBwBY+5Yys/aCJFfNrb8yyau6e2OSG5OcNbWfleTG7j4uyaumfgAA7IVFhbWqOibJKUneMK1XkqcmeefU5YIkp03Lp07rmbafNPUHAGCJFjuz9uokP5fkH6b1ByS5qbvvmNa3Jzl6Wj46ybYkmbbfPPUHAGCJ9hjWqur7k1zf3ZfPNy/QtRexbX6/Z0/Xu92wdevWRRULsJYZF4GFLGZm7UlJfqCqrknytsxOf746yWFVtW7qc0ySa6fl7UmOTZJp+32TfGXnnXb3ed29sbuPWL9+/T4dBMBaYFwEFrLHsNbdv9Ddx3T3hiTPSvL+7n5ukg8keebUbVOSi6bli6f1TNvf3913mlkDAGDP9uV71n4+yYuqaktm16SdP7Wfn+QBU/uLkmzetxIBAA5e6/bc5Z9092VJLpuWv5DkxAX63J7k9GWoDQDgoOcOBgAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBrVvtApZqw+ZLVmzf15x7yortGwBgb5hZAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMLAD7nZTK8mtrACA0ZhZAwAYmJk11hwzpACsJWbWAAAGJqwBAAxMWAMAGJhr1tYA12gBwNplZg0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiY71ljt1bqO9x8fxsALI6wBsAu+dJtWH1OgwIADMzMGqtiJd+tA8BaIqzBIFwfCMBCnAYFABiYsAYAMDBhDQBgYMIaAMDAhDUAgIHtMaxV1d2r6mNV9emq+kxV/dLU/q1V9dGqurqq3l5Vh07td5vWt0zbN6zsIQAArF2LmVn7uyRP7e5HJ3lMkqdX1ROTvDLJq7p7Y5Ibk5w19T8ryY3dfVySV039AADYC3sMaz1z67R61+nRSZ6a5J1T+wVJTpuWT53WM20/qapq2SoGADiILOpLcavqkCSXJzkuyeuSfD7JTd19x9Rle5Kjp+Wjk2xLku6+o6puTvKAJF9exrphVbjzAgD726I+YNDd3+zuxyQ5JsmJSR6xULfpeaFZtN65oarOnq53u2Hr1q2LrRdgzTIuAgtZ0u2muvumqrosyROTHFZV66bZtWOSXDt1257k2CTbq2pdkvsm+coC+zovyXlJcsIJJ9wpzK01ZmSAPTnYxkVgcRbzadAjquqwafkeSZ6W5KokH0jyzKnbpiQXTcsXT+uZtr+/uw06AAB7YTEzaw9McsF03dpdklzY3e+uqs8meVtVvSLJp5KcP/U/P8n/rKotmc2oPWsF6gYAOCjsMax195VJHrtA+xcyu35t5/bbk5y+LNUBABzk3MEAAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAa2brULAIDltGHzJSu272vOPWXF9g27YmYNAGBgwhoAwMCcBgVgVazk6UpYS8ysAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAfBoU1jhfEApwYDOzBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwPYY1qrq2Kr6QFVdVVWfqaoXTO33r6r3VtXV0/P9pvaqqtdU1ZaqurKqjl/pgwAAWKsWM7N2R5Kf6e5HJHliknOq6pFJNid5X3dvTPK+aT1JTk6ycXqcneT1y141AMBBYo9hrbu/2N2fnJZvSXJVkqOTnJrkgqnbBUlOm5ZPTfKWnvlIksOq6oHLXjkAwEFgSdesVdWGJI9N8tEkR3X3F5NZoEty5NTt6CTb5l62fWoDAGCJ1i22Y1XdO8nvJXlhd3+1qnbZdYG2XmB/Zyd5cZLDjjjiiMWWAbBmGRfHt2HzJSu272vOPWXF9s2BbVEza1V118yC2lu7+11T83U7Tm9Oz9dP7duTHDv38mOSXLvzPrv7vO7e2N1HrF+/fm/rB1gzjIvAQhbzadBKcn6Sq7r71+c2XZxk07S8KclFc+1nTJ8KfWKSm3ecLgUAYGkWcxr0SUmel+TPq+qKqe0lSc5NcmFVnZVka5LTp22XJnlGki1Jbkty5rJWDABwENljWOvuD2Xh69CS5KQF+neSc/axLgAA4g4GAABDE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxs3WoXAAAkGzZfsmL7vubcU1Zs36w8M2sAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAAMT1gAABiasAQAMbN1qFwAArKwNmy9Zkf1ec+4pK7Jf/jkzawAAAxPWAAAGJqwBAAxsj2Gtqt5YVddX1V/Mtd2/qt5bVVdPz/eb2quqXlNVW6rqyqo6fiWLBwBY6xYzs/bmJE/fqW1zkvd198Yk75vWk+TkJBunx9lJXr88ZQIAHJz2GNa6+0+SfGWn5lOTXDAtX5DktLn2t/TMR5IcVlUPXK5iAQAONnt7zdpR3f3FJJmej5zaj06yba7f9qkNAIC9sNwfMKgF2nrBjlVnT9e83bB169ZlLgPgwGNcBBayt2Htuh2nN6fn66f27UmOnet3TJJrF9pBd5/X3Ru7+4j169fvZRkAa4dxEVjI3t7B4OIkm5KcOz1fNNf+k1X1tiRPSHLzjtOlAMDaslJ3RkjcHWHeHsNaVf1ukicnObyqtif5r5mFtAur6qwkW5OcPnW/NMkzkmxJcluSM1egZgCAg8Yew1p3P3sXm05aoG8nOWdfiwIAYMYdDAAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAA9vb200BAKwYt7L6J2bWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBrVvtAgAA9qcNmy9ZsX1fc+4py75PM2sAAAMT1gAABiasAQAMTFgDABiYsAYAMDBhDQBgYMIaAMDAhDUAgIEJawAAAxPWAAAGJqwBAAxMWAMAGJiwBgAwMGENAGBgwhoAwMCENQCAgQlrAAADE9YAAAYmrAEADExYAwAYmLAGADAwYQ0AYGDCGgDAwIQ1AICBCWsAAANbkbBWVU+vqr+sqi1VtXklfgYAwMFg2cNaVR2S5HVJTk7yyCTPrqpHLvfPAQA4GKzEzNqJSbZ09xe6+++TvC3JqSvwcwAA1ryVCGtHJ9k2t759agMAYInWrcA+a4G2vlOnqrOTvDjJYUlur6rPrEAt++LwJF9e7SKWwVo5jsSxDKdeueTjePBK1bIWDD4urol/sxPHMp61chxLHRcXNSZW951y1D6pqu9I8rLu/r5p/ReSpLv/27L+oBVWVZ/o7hNWu459tVaOI3EsI1orx8GeraW/a8cynrVyHMnKHMtKnAb9eJKNVfWtVXVokmcluXgFfg4AwJq37KdBu/uOqvrJJP8nySFJ3tjdI03lAwAcMFbimrV096VJLl2Jfe9H5612ActkrRxH4lhGtFaOgz1bS3/XjmU8a+U4khU4lmW/Zg0AgOXjdlMAAAMT1gAABiasrXFV9ZbVruFgVlWHVtUZVfW0af05VfXaqjqnqu662vXBwcaYuPqMi0vnmrVJVT08szstfLS7b51rf3p3/9HqVbZ4VbXzV6RUkqckeX+SdPcP7Pei9lJVPSHJVd391aq6R5LNSY5P8tkkv9LdN69qgYtUVW/N7IM890xyU5J7J3lXkpMy+/3btIrlwS4ZE8eyVsbExLi4N4S1JFX100nOSXJVksckeUF3XzRt+2R3H7+a9S1WVX0ys1/cN2R214hK8ruZfddduvuDq1fd0kzf3P7o6atgzktyW5J3ZvbL/Oju/sFVLXCRqurK7v72qlqX5G+SPKi7v1lVleTT3f3tq1wi3IkxcTxrZUxMjIt7Y0W+uuMA9Pwkj+vuW6tqQ5J3VtWG7v6NLHz7rFGdkOQFSV6a5MXdfUVVff1AGpDm3KW775iWT5j7z+FDVXXFahW1F+4yfTn0vTJ7F3nfJF9JcrckB9R0f1XdN8kvJDktyRFT8/VJLkpybnfftFq1seyMieNZK2NiYlxcMteszRyyY5q/u69J8uQkJ1fVr+cAGpi6+x+6+1VJzkzy0qp6bQ7cQP4XVXXmtPzpqjohSarqoUm+sXplLdn5ST6X5IrM/sN4R1X9dmZ3+njbaha2Fy5McmOSJ3f3A7r7AZmdUroxyTtWtTKWmzFxPGtlTEyMi0vmNGiSqnp/khd19xVzbeuSvDHJc7v7kFUrbh9U1SlJntTdL1ntWpZqerfyG0m+K7Mb4h6fZNv0+Onu/vQqlrckVfWgJOnua6vqsCRPS7K1uz+2upUtTVX9ZXc/bKnbOPAYE8ezlsbExLi45J8jrCVVdUySO7r7Swtse1J3f3gVyiJJVd0nyUMyeze8vbuvW+WSDlpV9Z4kf5zkgh1/D1V1VJIfTfI93f20VSyPZWRMHJcxcSz7a1wU1oBFqar7ZfYJtFOTHDk1X5fk4syuzbhxtWoDWA37a1wU1oB9VlVndvebVrsOgFEs57gorAH7rKq2dvf61a4DYBTLOS4eqJ+KAfazqrpyV5uSHLU/awEYwf4aF4U1YLGOSvJ9mX0kfV4l+dP9Xw7Aqtsv46KwBizWu5Pce/7rHHaoqsv2fzkAq26/jIuuWQMAGJg7GAAADExYAwAYmLAGADAwYQ0AYGDCGgDAwP4/nYj7cn+WBzwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = compas.query(\"race in ['Caucasian', 'African-American']\").hist(\"decile_score\", by = \"race\",\n",
    "                                                                     figsize=(10,7), sharey=True,\n",
    "                                                                    grid = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>69</td>\n",
       "      <td>34</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c_charge_degree</th>\n",
       "      <td>F</td>\n",
       "      <td>F</td>\n",
       "      <td>F</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>race</th>\n",
       "      <td>Other</td>\n",
       "      <td>African-American</td>\n",
       "      <td>African-American</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age_cat</th>\n",
       "      <td>Greater than 45</td>\n",
       "      <td>25 - 45</td>\n",
       "      <td>Less than 25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>score_text</th>\n",
       "      <td>Low</td>\n",
       "      <td>Low</td>\n",
       "      <td>Low</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <td>Male</td>\n",
       "      <td>Male</td>\n",
       "      <td>Male</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>priors_count</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_b_screening_arrest</th>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>decile_score</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>is_recid</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>two_year_recid</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c_jail_in</th>\n",
       "      <td>2013-08-13 06:03:42</td>\n",
       "      <td>2013-01-26 03:45:27</td>\n",
       "      <td>2013-04-13 04:58:34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c_jail_out</th>\n",
       "      <td>2013-08-14 05:41:20</td>\n",
       "      <td>2013-02-05 05:36:53</td>\n",
       "      <td>2013-04-14 07:02:04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>length_of_stay</th>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                           0                    1  \\\n",
       "age                                       69                   34   \n",
       "c_charge_degree                            F                    F   \n",
       "race                                   Other     African-American   \n",
       "age_cat                      Greater than 45              25 - 45   \n",
       "score_text                               Low                  Low   \n",
       "sex                                     Male                 Male   \n",
       "priors_count                               0                    0   \n",
       "days_b_screening_arrest                   -1                   -1   \n",
       "decile_score                               1                    3   \n",
       "is_recid                                   0                    1   \n",
       "two_year_recid                             0                    1   \n",
       "c_jail_in                2013-08-13 06:03:42  2013-01-26 03:45:27   \n",
       "c_jail_out               2013-08-14 05:41:20  2013-02-05 05:36:53   \n",
       "length_of_stay                             0                   10   \n",
       "\n",
       "                                           2  \n",
       "age                                       24  \n",
       "c_charge_degree                            F  \n",
       "race                        African-American  \n",
       "age_cat                         Less than 25  \n",
       "score_text                               Low  \n",
       "sex                                     Male  \n",
       "priors_count                               4  \n",
       "days_b_screening_arrest                   -1  \n",
       "decile_score                               4  \n",
       "is_recid                                   1  \n",
       "two_year_recid                             1  \n",
       "c_jail_in                2013-04-13 04:58:34  \n",
       "c_jail_out               2013-04-14 07:02:04  \n",
       "length_of_stay                             1  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\stats\\stats.py:1713: 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",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl81NW9//HXZ7InQCAhCYFAEssaFhcQUdzFiktFW6zaqtTa2vur/ry2ve3V9tfeq7e9rVdv7ea19cpt1Vp3a7nWlqpo6wqyCMimrCGEkATIRvbk/P6YicY0gUkyk+8s7+fjkQcz3zmT+TB+ec/xzPmeY845REQkPvi8LkBERIaOQl9EJI4o9EVE4ohCX0Qkjij0RUTiiEJfRCSOKPRFROKIQl9EJI4o9EVE4kii1wX0NHr0aFdUVOR1GSIiUWXNmjXVzrmcY7WLuNAvKipi9erVXpchIhJVzGxPMO00vCMiEkcU+iIicUShLyISRxT6IiJxRKEvIhJHFPoiInFEoS8iEkcU+h7TdpUiMpQi7uKseOCc4zdv7ubVbVW8u7eG/MxUfnnNbIpGZ3hdmojEOPX0PfDs2n3c8b+bKa9p4oLpeRyoa+ay/3qDt3ce9Lo0EYlxCv0hVtPYyg9e2MJJE0ay/NYz+Y/Fx/PcTfPJzkjm2qUrWVd62OsSRSSGKfSH2F1/3kptUxvfv2wmPp8BUJidwbP/Zz6jh6Vw+7Mbaevo9LhKEYlVCv0htLb0MI+t2ssX5xdRMnbExx7LTE/ijkuns7Wingdf2+VRhSIS6xT6Q+jxVaUMT03k1gWTe338k9PHcMH0PH768vuUHmwc4upEJB4o9IdIZ6djxdYqzpqcQ0ZK35Om7rh0Bok+H//+wpYhrE5E4oVCf4hs3FdLdUML503LPWq7MZmpfHF+EX/eVMH7B+qHqDoRiRcK/SHy8tZKfAZnTT566ANcP7+Y9OQE/uuV7UNQmYjEE4X+EFmx9QAnTRhFVkbyMduOykjmmnmFLFtfzu7qI0NQnYjEC4X+EKiobea9fXWce4yhne6+dHoxiQk+fvnXHWGsTETijUJ/CLyyrRKA86bmBf2c3BGpXDlnPM+sLaOyrjlcpYlInFHoD4GXt1QybmQak/OG9et5Xzy9mLYOx2Or9oapMhGJNwr9MHPOsWrXQc6cPBoz69dzi0dncObkHB5duUdX6YpISCj0w2xfTRN1ze1MH5s5oOcvObWQyvoWlm+qCHFlIhKPFPphtmW/f659z2UXgnX2lFzGZ6Xx8Ft7QlmWiMQphX6YbS6vwwymjhk+oOcn+Ixr5xWyatchtuyvC3F1IhJvFPphtmV/HcXZGaQnD3y/ms/OGU9yoo/HVpWGsDIRiUcK/TDbvL+OafkDG9rpMjI9mQtnjOG5dftobusIUWUiEo+CCn0zW2hm28xsu5nd1svjKWb2RODxlWZW1OPxCWbWYGb/FJqyo0N9cxulhxoHPJ7f3ZVzxlPX3K4vdEVkUI4Z+maWANwHXAiUAFebWUmPZjcAh51zE4F7gbt6PH4v8KfBlxtdtlb4v8Sdlj+w8fzu5h2XzfisNJ54R3P2RWTggunpzwW2O+d2OudagceBRT3aLAIeCtx+GjjPApPSzewyYCewKTQlR4/N5f4vXkvyBzZdszufz/js7PG8ueOg1toXkQELJvTHAd27l2WBY722cc61A7VAtpllAP8M3DH4UqPP5vI6sjKSyRuREpLft3hOAT6Dp9aoty8iAxNM6Pd2GakLss0dwL3OuYajvoDZjWa22sxWV1VVBVFSdNhSUce0/OH9vhK3L/mZaZw5OYdn1pTR2dnzP4GIyLEFE/plwPhu9wuA8r7amFkikAkcAk4B/sPMdgO3At82s5t7voBz7gHn3Bzn3JycnJx+/yUiUXtHJ1sr6ikZ5Mydnj59UgHltc2s2n0opL9XROJDMKH/DjDJzIrNLBm4CljWo80yYEng9mJghfM7wzlX5JwrAn4C/Ltz7hchqj2i7aw+Qmt756Cna/Z0/rQ8MpITeG7dvpD+XhGJD8cM/cAY/c3AcmAL8KRzbpOZ3WlmlwaaLcU/hr8d+Drwd9M6482OSv+I1uS8wc/c6S4tOYGFM/L548b9mrMvIv0W1GWizrkXgBd6HPtet9vNwBXH+B3/OoD6olbpIf8MmwnZ6SH/3ZefOI5n1paxYmslF83MD/nvF5HYpStyw2TPoUZGpicxIjUp5L/71E9kkzs8hd9riEdE+kmhHyZ7DzVSmBX6Xj74F2FbdMJYXt1WyeEjrWF5DRGJTQr9MCk91Mj4MIU+wKITxtHW4bQsg4j0i0I/DNo7Otl3uIkJYQz96WNHUJSdzh837g/ba4hI7FHoh8H+2mbaOx2FYfgSt4uZcdHMfN7ccZBDGuIRkSAp9MOga+ZOOId3AC6elU9Hp+MvGuIRkSAp9MPgw+maYQ79knwN8YhI/yj0w6D0UCNJCUZ+ZlpYX8fMuHiWhnhEJHgK/TAoPdhIwah0EnyhWWjtaC6eOZaOTs3iEZHgKPTDINzTNbublj+c4tEZ/Ok9hb6IHJtCPwxKDzUyISu8QztdzIzzS/J4a0c19c1tQ/KaIhK9FPohVtvYRm1TG4VZGUP2mgum5dHW4fjb+9VD9poiEp0U+iE2VNM1u5tdOIqsjGRe3KwhHhE5OoV+iA3VdM3uEnzGuVNzWbG1kraOziF7XRGJPgr9EAvnkspHs2BaHnXN7byjHbVE5CgU+iFWeqiR7IxkhqUEtVVByJw5eTTJiT5e3HxgSF9XRKKLQj/EymuaGDdqaGbudJeenMjpE0fz0pYDOKdN00Wkdwr9EKuobSZvRKonr31+SR57DzWx7UC9J68vIpFPoR9iFXXN5Gd6E/rnTcsF4MVNGuIRkd4p9EOoqbWD2qY2z3r6ucNTOWH8SF7aotAXkd4p9EOooq4ZwLOePviHeNaX1XIgUIuISHcK/RCqqPUH7RiPevoAnyzJA9AsHhHplUI/hCrqmgDI87CnPzF3GIXZ6RriEZFeKfRDqKK2BfC2p29mnD8tjze3H6Shpd2zOkQkMin0Q6iitonhqYlkDPGFWT2dX5JHa0cnf3u/ytM6RCTyKPRDqKKu2dNefpfZhaMYmZ7ESxrXF5EeFPohVFHXwhgPx/O7JCb4OHdKLiu2VdKuBdhEpBuFfghV1DZFRE8f/EM8NY1trNlz2OtSRCSCKPRDpL2jk6r6Fk/n6Hd3xuQckhN8msUjIh+j0A+RqoYWOp230zW7G5aSyKmfyObFzVqATUQ+otAPkUi4MKunBSV57D7YyI6qBq9LEZEIodAPka5lDyLhi9wuC7oWYNtc6XElIhIpFPohsj8Ce/r5mWnMHJepcX0R+ZBCP0Qq6ppJTvCRlZHsdSkfs2BaHmtLD1Pd0OJ1KSISART6IVJR20xeZgpm5nUpH7OgJBfnYMUWDfGIiEI/ZCpqI+Nq3J5K8kcwNjOVFzXEIyIo9EPmQF0zYzKHfm/cYzEzFpTk8doHVTS3dXhdjoh4LKjQN7OFZrbNzLab2W29PJ5iZk8EHl9pZkWB43PN7N3Az3ozuzy05UcG5xz7a5sZMyLF61J6dX5JHs1tnbyxvdrrUkTEY8cMfTNLAO4DLgRKgKvNrKRHsxuAw865icC9wF2B4+8Bc5xzJwALgV+ZmbdLUIZBXVM7Le2dnm2TeCynFGczLCVRs3hEJKie/lxgu3Nup3OuFXgcWNSjzSLgocDtp4HzzMycc43Oua5F3VOBmLw0tKrBP10zZ3hk9vSTE32cNSWHl7ZU0tkZk/8JRCRIwYT+OGBvt/tlgWO9tgmEfC2QDWBmp5jZJmAj8A/dPgRiRmW9fzpkpIY+wPnT8qiqb2Hd3hqvSxERDwUT+r3NQezZXeyzjXNupXNuOnAycLuZ/d0YiJndaGarzWx1VVX0bfxR3dAKQM6wyA39c6flkpzg44WN+70uRUQ8FEzolwHju90vAMr7ahMYs88EDnVv4JzbAhwBZvR8AefcA865Oc65OTk5OcFXHyGqoqCnPyI1iTMnj+ZPG/driEckjgUT+u8Ak8ys2MySgauAZT3aLAOWBG4vBlY451zgOYkAZlYITAF2h6TyCFJV30JSgpGZluR1KUd10cx8ymubebdMQzwi8eqYoR8Yg78ZWA5sAZ50zm0yszvN7NJAs6VAtpltB74OdE3rPB1Yb2bvAr8Hvuqci7l5g9UNLeQMi7yrcXtaUJJHcoKPP27QEI9IvApq+qRz7gXghR7HvtftdjNwRS/PewR4ZJA1Rryq+hZGR/DQTpcRqUmcMck/xPOdi6bh80X2h5SIhJ6uyA2BqvqWiP4StzsN8YjEN4V+CFQ3tET0l7jddQ3xPL9eQzwi8UihP0gdnY6DR1oZHSU9/cy0JM6eksOy9eW0d3R6XY6IDDGF/iAdbmylo9NFTU8f4NMnjaO6oYXXtRaPSNxR6A9S1+Yk0RT650zNJTMtiefW7fO6FBEZYgr9Qeq6MCtahncAUhITuHhWPss3HeBIS8ytiiEiR6HQH6RouBq3N58+cRxNbR38+b0Kr0sRkSGk0B+kaBzeAZhdOIoJWen8XkM8InFFoT9IVfUtpCb5yEhO8LqUfjEzLj9xHG/sqKbscKPX5YjIEFHoD1JVvX+OfqQvwdCbK+YUAPDk6jKPKxGRoaLQH6Sqhui5GrenglHpnDkph6dW79WcfZE4odAfpOr66LkwqzdXzx3P/tpm/vp+9O1jICL9p9AfpKooWoKhN+dNy2P0sBQeW7X32I1FJOop9AehraOTQ0daozr0kxJ8LJ5dwCvbKqmobfa6HBEJM4X+IBw64t8mMZqHdwCuOnk8HZ2Ox98p9boUEQkzhf4gROuFWT0Vjc7gzMk5/G5lKW36Qlckpin0ByFWQh9gyamFVNa38JdNB7wuRUTCSKE/CFVdV+NG+fAOwNlTchmflcZDb+32uhQRCSOF/iBE42JrfUnwGdfOK2TVrkNs2V/ndTkiEiYK/UGoqm9heEoiaVG2BENfPjtnPCmJPh5+a4/XpYhImCj0B6G6ITo2RA/WyPRkFp0wlufW7aO2qc3rckQkDBT6gxBNG6IH67pTi2hq6+DpNVqPRyQWKfQHIdqvxu3NjHGZzC4cxSNv7aaz03ldjoiEmEJ/EKrrYy/0Aa47tZDdBxv52wdaj0ck1ij0B6i5rYO65nZGD0v2upSQu3BGPqOHpfCIvtAViTkK/QGK1h2zgpGc6ONzc8ezYlslpQe1wYpILFHoD1B1g3/dnVgMfYDPnVKIz4zfrlRvXySWKPQHKJYuzOrNmMxULpiexxPv7KWptcPrckQkRBT6AxRL6+705bpTi6htauN/15d7XYqIhIhCf4C6xvSzM2I39E8pzmJK3nB+8+ZunNP0TZFYoNAfoKr6FkamJ5GcGLtvoZlx3WmFbN5fx9rSw16XIyIhELuJFWaxeDVuby47YRzDUxK1Ho9IjFDoD1B1DF6N25uMlEQWzynghY37qazXdooi0U6hP0BVDS0xO3Onp2vnFdLW4Xhcm6eLRD2F/gBVxegSDL05LmcYZ07O4dGVe7SdokiUU+gPwJGWdhpbO+Im9AGum1fIgboWXtys7RRFoplCfwC6pmvGy/AOwDlTcykYlcZDb+72uhQRGQSF/gDEw4VZPXVtp7hy1yG2Vmg7RZFoFVTom9lCM9tmZtvN7LZeHk8xsycCj680s6LA8fPNbI2ZbQz8eW5oy/dGdQxtiN4fXdspavVNkeh1zNA3swTgPuBCoAS42sxKejS7ATjsnJsI3AvcFTheDXzKOTcTWAI8EqrCvfThujvDY29Z5aMZlZHMpceP5dm12k5RJFoF09OfC2x3zu10zrUCjwOLerRZBDwUuP00cJ6ZmXNunXOua+GWTUCqmUV997iqvgWfxfYSDH1Zcpp/O8VntJ2iSFQKJvTHAd0naJcFjvXaxjnXDtQC2T3afAZY55xr6fkCZnajma02s9VVVZG/W1NVQytZGSkk+MzrUobcjHGZnDRhJI+8vUfbKYpEoWBCv7dk6/mv/ahtzGw6/iGfr/T2As65B5xzc5xzc3JycoIoyVtV9S0xuWNWsJacVsSu6iPaTlEkCgUT+mXA+G73C4Cea+1+2MbMEoFM4FDgfgHwe+A659yOwRYcCWJxQ/T+uHBGPrnDU1j6+i6vSxGRfgom9N8BJplZsZklA1cBy3q0WYb/i1qAxcAK55wzs5HAH4HbnXNvhKpor8XqhujBSk708YX5Rbz2QTVb9mv6pkg0OWboB8bobwaWA1uAJ51zm8zsTjO7NNBsKZBtZtuBrwNd0zpvBiYC3zWzdwM/uSH/Wwwh5xxV9S3kDk/1uhRPfX5uIenJCfz3azu9LkVE+iExmEbOuReAF3oc+163283AFb087/vA9wdZY0SpaWyjtaMzrnv6AJnpSXx2znh++/YevnnBFPIz07wuSUSCoCty+6kqcGFWbpyHPsANpxfT6Ry/eWO316WISJAU+v1UWafQ7zI+K52LZ43lkbf3fHiVsohENoV+P3VtJJI7Ir7H9LvcumASzW0d3P9qTEzMEol5Cv1+qozDxdaO5hM5w1g8u4BH3t5DeU2T1+WIyDEo9Pupqr6F9OQEhqUE9R14XLjlvEng4OcrPvC6FBE5BoV+P1XWt2g8v4eCUel87pQJPLm6TMsui0Q4hX4/VdY1a2inF7ecN4mRaUl86+kNtGtLRZGIpdDvJ12Y1busjGTuWDSdDWW1/PdrWp5BJFIp9PspnjZE76+LZ+azcPoY7n3pfbZXNnhdjoj0QqHfD02tHdS3tJM7QqHfGzPjzsumk5GcwI0Pr+bQkVavSxKRHhT6/dA1Rz/etknsj9zhqTxw3Rz21TRx/W/eobG13euSRKQbhX4/dM3R14VZR3dyURY/u/pENpbV8JVH1lDXrK0VRSKFQr8fuvbG1ZTNY7tg+hh+9JlZvLXjIJf87HXe21frdUkigkK/XyrrAkswKPSD8tk543n8xnm0tnfy6fvf5J7l27ShuojHFPr9UFnfQqLPGJUev1sl9tecoiz+eMvpXDB9DL94ZTtn3LWCn7z0/ocfoCIytBT6/VBZ38LoYSn44nBD9MHIHpbCz68+kT/ecjpzi7P4yUsfcNqPVnDT79aycudBnNMG6yJDRQvI9ENVfYumaw7C9LGZPLjkZHZVH+G3b+/hqdV7+eOG/UzJG84Npxdz2YnjSE5UP0QknPQvrB+07k5oFI/O4LuXlLDy2wu46zMz8fmMbz2zgbPvfoVH3tpNm5ZxEAkbhX4/VNVr3Z1QSktO4MqTJ/DCLafz6+tPZuzINL77h01c8JO/8dLmAxr2EQkDhX6Q2js6OXiklRytuxNyZsY5U3J56h9O5cHr5gDwpYdX8/kHV7KpXFM9RUJJoR+kg0dacU7TNcPJzFhQksfyW8/kjkuns2V/HZf8/HW++dR6Dmi2j0hIKPSDpL1xh05Sgo8lpxXx6jfP4ctnHMcf3i3n7Ltf5d4X36emUev5iAyGQj9IFYGeZp6WYBgymWlJfPuiabz09bM4d2ouP335A+b98GW+8/uNrCs9TGenxvxF+ktTNoO0v9a//2v+SIX+UJuQnc59nz+J/1tRx69f381Ta8p4dGUpucNTOHNyDscXZDJ9XCYl+SNITUrwulyRiKbQD1J5TTNJCcboDA3veGXqmBHctXgWt180lVe2VfLi5gO8vOUAT68pAyDBZ0zKHcbswlGcMSmH0yZmMyI1yeOqRSKLQj9I+2ubGJOZqqtxI8DI9GQuP7GAy08swDlHeW0zG8tqeW9fLRv31fLcun08urKU5AQflxyfz/WnFTOzINPrskUigkI/SPtrmsnPTPO6DOnBzBg3Mo1xI9NYOGMMAG0dnawrreH5DeU8vaaMZ9fuY8G0XP7lU9MZn5XuccUi3tIXuUEqr21ibKbG86NBUoKPucVZ3LloBm9/+zy+tXAKb+44yIIf/5X7XtlOh74Aljim0A9CZ6fjQF0z+SPV0482I1KT+OrZE3n5G2exYFoedy/fxjUPrtS8f4lbCv0gVDe00NbhyFdPP2rlZ6Zx3+dP4p4rjufdvTVc9NPXeHvnQa/LEhlyCv0glNf6e4Ua049+i2cXsOzm+WSmJ3Ht0pU88U6p1yWJDCmFfhD21wTm6KunHxMm5Q3n91+dz7zjsvnnZzbywz9t0YVeEjcU+kHo6umP1Zh+zMhMS+LXXziZa+cV8qu/7uSfnlqvJZ0lLmjKZhD21zSRkuhjVLou9IkliQk+7lw0nbwRKdzzl/c51NjKL6+Zrat6Jaappx+E/bXNjB2ZhpkuzIo1ZsbN507ih5+eyV/fr+JLD62mqbXD67JEwkahH4Ty2iaN58e4q+dO4O7Fx/PGjmq+/LCCX2KXQj8Iuho3PiyeXcA9geC/6XdradcYv8SgoELfzBaa2TYz225mt/XyeIqZPRF4fKWZFQWOZ5vZK2bWYGa/CG3pQ6O9o5PK+mbGanXNuPCZ2QXcuWgGK7ZW8t0/bNKWjRJzjvlFrpklAPcB5wNlwDtmtsw5t7lbsxuAw865iWZ2FXAXcCXQDHwXmBH4iTqV9S10Os3RjyfXziukvKaJ+1/dQcGoNG46Z6LXJYmETDA9/bnAdufcTudcK/A4sKhHm0XAQ4HbTwPnmZk55444517HH/5RSevox6dvfnIKi04Yy93Lt/HKtkqvyxEJmWBCfxywt9v9ssCxXts459qBWiA7FAV6rbwmMEdfPf244vMZP/r0LKaOGc7XnniXssONXpckEhLBhH5v8xR7DnQG06bvFzC70cxWm9nqqqqqYJ82JNTTj19pyQncf81sOjocNz26lpZ2zeiR6BdM6JcB47vdLwDK+2pjZolAJnAo2CKccw845+Y45+bk5OQE+7QhUV7TzLCURO3AFKeKR2dw9xWzWF9Wy70vfuB1OSKDFkzovwNMMrNiM0sGrgKW9WizDFgSuL0YWOFiZNpD2eEmzdyJcwtn5HP13PH86m87WL076L6MSEQ6ZugHxuhvBpYDW4AnnXObzOxOM7s00GwpkG1m24GvAx9O6zSz3cCPgS+YWZmZlYT47xBWu6obKB6d4XUZ4rHvXFxCwag0vvHUeo60tHtdjsiABTVP3zn3gnNusnPuE865HwSOfc85tyxwu9k5d4VzbqJzbq5zbme35xY557Kcc8OccwU9pnpGtI5Ox95DTRQp9OPesJRE/vOKEyg91MgP/7TF63JEBkxX5B5FeU0TrR2dHKfQF2BucRY3zC/mt2+XagMWiVoK/aPYWX0EgKJshb74feOTU5iQlc5tz2zQ+jwSlRT6R7E7EPoa05cuackJ/OgzM9l9sJGfvPS+1+WI9JtC/yh2VR8hIzmBnOEpXpciEeS0T4zm6rkT+O/XdrJ+b43X5Yj0i0L/KHZVH6FodIbW0Ze/c/tFU8kdnso/P7OB1natxinRQ6F/FLsPHtHMHenViNQkfnD5DLZW1HP/qzu8LkckaAr9PrS2d1J2uEkzd6RP503LY9EJY/nFKx+wraLe63JEgqLQ78Pew410dDrN3JGj+t4lJQxPTeJbz2ygozMmLkKXGKfQ70PXzB0N78jRZA9L4V8vnc76vTX8+o1dXpcjckwK/T7sCoS+hnfkWD41K58F03K55y/b2HPwiNfliByVQr8Pu6qPkJmWxKiMZK9LkQhnZnz/spkk+Xx862kN80hkU+j3QTN3pD/GZKbyvU+VsHLXIX75V83mkcil0O/D7upGDe1IvyyeXcAls/L58Yvvs7b0sNfliPRKod+LxtZ2ymubNHNH+sXM+MHlM8nPTOWWx9ZR29TmdUkif0eh34tN5XU4B9PHjvC6FIkymWlJ/PSqE6mobeZrT7xLp8b3JcIo9HvRtZ7KrPGZHlci0Wh24Sj+5VMlrNhayb1alE0iTKLXBUSiDWW15Gemkjtc2yTKwFwzr5D39tXx8xXbmZY/gotm5ntdkgignn6vNpTVMKtAvXwZODPjzsumc9KEkXztiXdZtUt760pkUOj3UNvYxu6DjcwqGOl1KRLlUhITWLrkZApGpfGlh95ha0Wd1yWJKPR72rDPP55/vEJfQmBURjIP33AK6cmJXLd0FTuqGrwuSeKcQr+HDWW1AMzU8I6EyLiRaTx8w1w6nePKX72lHr94SqHfw7t7aygenUFmWpLXpUgMmZw3nMdvPJUEn3HVA2+zThdviUcU+j3oS1wJl4m5w3jyK6cyPDWRKx94m9+vK/O6JIlDCv1uDtQ1c6CuReP5EjaF2Rn84abTA7N61vNvz2+mua3D67Ikjij0u+m6KOt4XZQlYZSVkcwjN5zCklMLWfr6Li79xeu8t6/W67IkTij0u3ljezWpST6mj1XoS3glJfi4Y9EMfn39ydQ0trHovjf43h/e42BDi9elSYxT6Ad0djr+9F4FZ0/OJTUpwetyJE6cMyWXF792FledPJ5HV5Zy1t2v8uO/bKNa4S9hotAPWLf3MJX1LVw4c4zXpUicyUxP4geXz2T5rWcyf2I2P1uxndN+tILbn93ApnIN+0hoae2dgBc2VpCc4OPcqblelyJxamLuMH517Ry2Vzaw9PWdPLt2H4+t2ssJ40fy+VMmcMmssaQl6/9CZXDMucha+nXOnDlu9erVQ/qazjlOv+sVpo4ZztIvnDykry3Sl9rGNp5ZW8ajK/ewo+oII1IT+czsAj5/ygQm5g73ujyJMGa2xjk351jt1NMHNu6rZV9NE7cumOR1KSIfykxP4ounF3P9/CJW7TrEoytL+e3be/j1G7s5pTiLz50ygYUzxpCSqN6/BE+hj39oJ9FnnF+S53UpIn/HzDjluGxOOS6b6oYSnl5TxmOrSvnHx98lKyOZq+eOZ8lpRVoKXIIS98M7Ta0dnHPPq0weM5yHvzh3yF5XZDA6Ox1v7Kjmkbf28OKWAyT5fFx+4ji+fGaxhn7ilIZ3gvTgazupqGvmZ1ef6HUpIkHz+YwzJuVwxqQcdlcf4cHXd/LU6jKeWL2XBdNy+fIZxzG3OAsz87pUiTBx3dM/UNfMOfe8ylmTc7j/mtlD8poi4XKwoYVH3t7Dw2/t4dBNB8+WAAAHU0lEQVSRVo4vyOTGMz/BwhljSPAp/GNdsD39uJ6n/59/2UZbRye3XTjV61JEBi17WAq3LpjMm7edy/cvm0FtUxs3/W4tZ9/zCr95YxeHj7R6XaJEgLjt6f9lUwVf+e0avnR6Md+5uCTsrycy1Do6HS9uPsADf9vB2tIaEn3G/ImjOb8kj/kTR1OUna7hnxgS0jF9M1sI/BRIAB50zv2ox+MpwMPAbOAgcKVzbnfgsduBG4AO4Bbn3PJ+/D3C4pk1ZXzrmQ3MKhjJLedpmqbEpgSfsXDGGBbOGMN7+2p5fsN+nt9Qzl+fqwIgb0QKJfkjmJo/goJRaYwZkUpe4Cc7IxmfhoRi0jF7+maWALwPnA+UAe8AVzvnNndr81VglnPuH8zsKuBy59yVZlYCPAbMBcYCLwGTnXN9riUbzp5+eU0TS1/fxdLXdzF/YjYPXDuHjJS4/y5b4ohzjl3VR3hjx0HW7jnMlv11bK9soL3z4zmQ4DMy05J6/RmZ7v9zRC/HMtOSSEtK0P9BeCCUPf25wHbn3M7AL34cWARs7tZmEfCvgdtPA78w/3/1RcDjzrkWYJeZbQ/8vreC/Yv0V2eno7Wjk/rmdg43tlJR28ym8jrWlR7m5a2VOOf47JwC/u2yGbqoReKOmXFczjCOyxnGtfMKAWjr6KS6oYWK2uYP95Q4UNdMbVPbhz+HG1vZffAItU1t1DW10XmUvmJSgpGZlkxmWmKPD4vkDz8ohqcmkpzgIzHBSPT5SEowEhN8JPrM/5MQOObramMk9Wif0HXM578dqg8a5xxdfWHXdR9wDhzdHgvc77rtf3/BMP+fgds+87/v1vW4xx+IwYT+OGBvt/tlwCl9tXHOtZtZLZAdOP52j+eOG3C1R7F+bw2fvv9NOvo4G8eNTOOL84tYcloRBaPSw1GCSFRKSvCRn5lGfmZaUO07Ox31Le3UdftQqG1qo6bx4/e7Hq9qaGF7VQO1jW3Ut7QTrq8RkxK6Bf+Hof3xUO566a4Rjo/uh6emvvg/EMBn9rEPiotn5vPjK08I62sHE/q9fSz1fIv6ahPMczGzG4EbA3cbzGxbEHX1yx7gTeD/hfoXw2igOvS/NirpvfiI3ouP6L34yFHfi/eBe68a8O8uDKZRMKFfBozvdr8AKO+jTZmZJQKZwKEgn4tz7gHggWAKjjRmtjqYcbR4oPfiI3ovPqL34iOR8F4EM0//HWCSmRWbWTJwFbCsR5tlwJLA7cXACuf//6dlwFVmlmJmxcAkYFVoShcRkf46Zk8/MEZ/M7Ac/5TN/3HObTKzO4HVzrllwFLgkcAXtYfwfzAQaPck/i9924GbjjZzR0REwiviLs6KNmZ2Y2B4Ku7pvfiI3ouP6L34SCS8Fwp9EZE4Etdr74iIxBuFfpDMbLyZvWJmW8xsk5n9Y+B4lpm9aGYfBP4c5XWtQ8XMEsxsnZk9H7hfbGYrA+/FE4Ev/mOemY00s6fNbGvg/Dg1Xs8LM/ta4N/He2b2mJmlxtN5YWb/Y2aVZvZet2O9ngvm9zMz225mG8zspKGoUaEfvHbgG865acA84KbAMhO3AS875yYBLwfux4t/BLZ0u38XcG/gvTiMf82lePBT4M/OuanA8fjfk7g7L8xsHHALMMc5NwP/xI+riK/z4jfAwh7H+joXLsQ/o3ES/uuU7h+SCv2XHOunvz/AH/CvR7QNyA8cywe2eV3bEP39CwIn8LnA8/gvxKsGEgOPnwos97rOIXgfRgC7CHw/1u143J0XfHRlfhb+mYHPAxfE23kBFAHvHetcAH6Ffx2zv2sXzh/19AfAzIqAE4GVQJ5zbj9A4M9c7yobUj8BvgV0Bu5nAzXOufbA/bAtuRFhjgOqgF8HhroeNLMM4vC8cM7tA+4BSoH9QC2whvg8L7rr61zobYmbsL83Cv1+MrNhwDPArc65Oq/r8YKZXQJUOufWdD/cS9N4mBqWCJwE3O+cOxE4QhwM5fQmMFa9CCjGv6puBv4hjJ7i4bwIhif/ZhT6/WBmSfgD/1Hn3LOBwwfMLD/weD5Q6VV9Q2g+cKmZ7QYexz/E8xNgZGAZDuhjyY0YVAaUOedWBu4/jf9DIB7PiwXALudclXOuDXgWOI34PC+66+tcCGqZmlBT6AcpsFT0UmCLc+7H3R7qvgTFEvxj/THNOXe7c67AOVeE/4u6Fc65zwOv4F+GA+LnvagA9prZlMCh8/BfgR535wX+YZ15ZpYe+PfS9V7E3XnRQ1/nwjLgusAsnnlAbdcwUDjp4qwgmdnpwGvARj4ax/42/nH9J4EJ+E/6K5xzhzwp0gNmdjbwT865S8zsOPw9/yxgHXCN8++lENPM7ATgQSAZ2Alcj79DFXfnhZndAVyJf7bbOuBL+Mep4+K8MLPHgLPxr6Z5APgX4Dl6ORcCH4y/wD/bpxG43jkX9r1iFfoiInFEwzsiInFEoS8iEkcU+iIicUShLyISRxT6IiJxRKEvIhJHFPoiInFEoS8iEkf+P7zJNsjtUue4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "display(compas.head(3).T)\n",
    "sns.kdeplot(np.array(compas_raw.age))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate synthetic data set\n",
    "\n",
    "In the chunk below, we generate the synthetic data as described by Lakkaraju et al.\n",
    "\n",
    "**Variables**\n",
    "\n",
    "* M = number of judges\n",
    "* subj = number of subjects assigned to each judge\n",
    "* betas $\\beta$ are coefficients\n",
    "* R = acceptance rates\n",
    "* X = invidual's features observable to all (models and judges)\n",
    "* Z = information observable for judges only\n",
    "* W = unobservable / inaccessible information\n",
    "* T = decisions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.random as npr\n",
    "\n",
    "npr.seed(0)\n",
    "\n",
    "nJudges_M = 100\n",
    "nSubjects_N = 500\n",
    "\n",
    "beta_X = 1.0\n",
    "beta_Z = 1.0\n",
    "beta_W = 0.2\n",
    "\n",
    "judgeID_J = np.repeat(np.arange(0, nJudges_M, dtype = np.int32), nSubjects_N)\n",
    "\n",
    "acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 10)\n",
    "\n",
    "acceptanceRate_R = np.repeat(acceptance_rates, nSubjects_N)\n",
    "\n",
    "X = npr.normal(size = nJudges_M * nSubjects_N)\n",
    "Z = npr.normal(size = nJudges_M * nSubjects_N)\n",
    "W = npr.normal(size = nJudges_M * nSubjects_N)\n",
    "\n",
    "probabilities_Y = 1 / (1 + np.exp(-(beta_X * X + beta_Z * Z + beta_W * W)))\n",
    "\n",
    "# 0 if P(Y = 0| X = x;Z = z;W = w) >= 0.5 , 1 otherwise\n",
    "result_Y = 1 - probabilities_Y.round()\n",
    "\n",
    "probabilities_T = 1 / (1 + np.exp(-(beta_X * X + beta_Z * Z)))\n",
    "probabilities_T += npr.normal(0, .1, nJudges_M * nSubjects_N)\n",
    "\n",
    "decision_T = np.zeros(nJudges_M * nSubjects_N) - 1\n",
    "\n",
    "tmp = pd.DataFrame(np.column_stack((judgeID_J, acceptanceRate_R, X,\n",
    "                                    Z, W, result_Y, probabilities_T, decision_T)),\n",
    "                   columns = [\"judgeID_J\", \"acceptanceRate_R\", \"X\",\n",
    "                              \"Z\", \"W\", \"result_Y\", \"probabilities_T\", \"decision_T\"])\n",
    "\n",
    "# Sort by judges then probabilities\n",
    "df = tmp.sort_values(by = [\"judgeID_J\", \"probabilities_T\"], ascending = False)\n",
    "\n",
    "# Iterate over the data. Subject is in the top (1-r)*100% if\n",
    "# his within-judge-index is over acceptance threshold times\n",
    "# the number of subjects assigned to each judge. If subject\n",
    "# is over the limit they are assigned a zero, else one.\n",
    "for i in range(nJudges_M * nSubjects_N):\n",
    "    index = i % nSubjects_N\n",
    "    if index < (1 - df.acceptanceRate_R[i]) * nSubjects_N:\n",
    "        df.decision_T[i] = 0\n",
    "    else:\n",
    "        df.decision_T[i] = 1  # TARKISTA!!!!!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic stats of the created data set.  We see that sensitivity is XX% and specifity YY%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0    26137\n",
      "1.0    23863\n",
      "Name: decision_T, dtype: int64\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>decision_T</th>\n",
       "      <th>0.0</th>\n",
       "      <th>1.0</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>result_Y</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.0</th>\n",
       "      <td>13119</td>\n",
       "      <td>12056</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.0</th>\n",
       "      <td>13018</td>\n",
       "      <td>11807</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "decision_T    0.0    1.0\n",
       "result_Y                \n",
       "0.0         13119  12056\n",
       "1.0         13018  11807"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(df.decision_T.value_counts())\n",
    "\n",
    "tab = df.groupby(['result_Y', 'decision_T']).size()\n",
    "tab.unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(25000, 8)\n",
      "(25000, 8)\n",
      "(12059, 8)\n",
      "(11804, 8)\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>decision_T</th>\n",
       "      <th>1.0</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>result_Y</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.0</th>\n",
       "      <td>6105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.0</th>\n",
       "      <td>5954</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "decision_T   1.0\n",
       "result_Y        \n",
       "0.0         6105\n",
       "1.0         5954"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Shuffle and split data set to test and train\n",
    "train, test = np.split(df.sample(frac = 1, random_state = 0), 2)\n",
    "\n",
    "print(train.shape)\n",
    "print(test.shape)\n",
    "\n",
    "train_labeled = train[train.decision_T == 1]\n",
    "test_labeled = test[test.decision_T == 1]\n",
    "\n",
    "print(train_labeled.shape)\n",
    "print(test_labeled.shape)\n",
    "\n",
    "tab = train_labeled.groupby(['result_Y', 'decision_T']).size()\n",
    "tab.unstack()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In his report Lakkaraju says that they used logistic regression.\n",
    "\n",
    "### Machine evaluation\n",
    "\n",
    "Next we train a logistic regression to predict the ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the class\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# instantiate the model (using the default parameters)\n",
    "logreg_machine = LogisticRegression(solver='lbfgs')\n",
    "\n",
    "# fit, reshape X to be of shape (n_samples, n_features)\n",
    "logreg_machine.fit(train_labeled.X.values.reshape(-1,1), train_labeled.result_Y)\n",
    "\n",
    "# predict probabilities and attach to data \n",
    "label_probabilities_machine = logreg_machine.predict_proba(test.X.values.reshape(-1,1))\n",
    "\n",
    "test['B_prob_0_machine'] = label_probabilities_machine[:, 0]\n",
    "test['B_prob_1_machine'] = label_probabilities_machine[:, 1]\n",
    "\n",
    "from sklearn import tree\n",
    "\n",
    "clf = tree.DecisionTreeClassifier()\n",
    "clf = clf.fit(train_labeled.X.values.reshape(-1,1), train_labeled.result_Y)\n",
    "\n",
    "preds = clf.predict_proba(test.X.values.reshape(-1,1))\n",
    "\n",
    "test['B_prob_0_tree'] = preds[:, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of contraction algorithm\n",
    "\n",
    "Below is an implementation of Lakkaraju's team's algorithm presented in [their paper](https://helka.finna.fi/PrimoRecord/pci.acm3098066). Relevant parameters to be passed to the fucntion are presented in the description."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contraction(df, judgeIDJ_col, decisionT_col, resultY_col, modelProbS_col, accRateR_col, r, binning = False):\n",
    "    '''\n",
    "    This is an implementation of the algorithm presented by Lakkaraju\n",
    "    et al. in their paper \"The Selective Labels Problem: Evaluating \n",
    "    Algorithmic Predictions in the Presence of Unobservables\" (2017).\n",
    "    \n",
    "    Parameters:\n",
    "    df = The (Pandas) data frame containing the data, judge decisions,\n",
    "    judge IDs, results and probability scores.\n",
    "    judgeIDJ_col = String, the name of the column containing the judges' IDs\n",
    "    in df.\n",
    "    decisionT_col = String, the name of the column containing the judges' decisions\n",
    "    resultY_col = String, the name of the column containing the realization\n",
    "    modelProbS_col = String, the name of the column containing the probability\n",
    "    scores from the black-box model B.\n",
    "    accRateR_col = String, the name of the column containing the judges' \n",
    "    acceptance rates\n",
    "    r = Float between 0 and 1, the given acceptance rate.\n",
    "    binning = Boolean, should judges with same acceptance rate be binned\n",
    "    \n",
    "    Returns:\n",
    "    u = The estimated failure rate at acceptance rate r.\n",
    "    '''\n",
    "    # Sort first by acceptance rate and judge ID.\n",
    "    sorted_df = df.sort_values(by = [accRateR_col, judgeIDJ_col], ascending = False)\n",
    "\n",
    "    if binning:\n",
    "        # Get maximum leniency\n",
    "        max_leniency = sorted_df[accRateR_col].values[0].round(1)\n",
    "\n",
    "        # Get list of judges that are the most lenient\n",
    "        most_lenient_list = sorted_df.loc[sorted_df[accRateR_col].round(1) == max_leniency, judgeIDJ_col]\n",
    "\n",
    "        # Subset to obtain D_q\n",
    "        D_q = sorted_df[sorted_df[judgeIDJ_col].isin(most_lenient_list.unique())]\n",
    "    else:\n",
    "        # Get most lenient judge\n",
    "        most_lenient_ID = sorted_df[judgeIDJ_col].values[0]\n",
    "        \n",
    "        # Subset\n",
    "        D_q = sorted_df[sorted_df[judgeIDJ_col] == most_lenient_ID]\n",
    "    \n",
    "    R_q = D_q[D_q[decisionT_col] == 1]\n",
    "\n",
    "    R_sort_q = R_q.sort_values(by = modelProbS_col, ascending = False)\n",
    "\n",
    "    number_to_remove = int(np.round((1 - r) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])))\n",
    "\n",
    "    R_B = R_sort_q[number_to_remove:R_sort_q.shape[0]]\n",
    "\n",
    "    return np.sum(R_B[resultY_col] == 0) / D_q.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validation analysis\n",
    "\n",
    "Below we try to replicate the results obtained by Lakkaraju et alii in their paper linked above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAG5CAYAAADRUnNdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XlY1NUawPHvYREQccUdFTUtFRQVxKXct0pNy8r0umaIu7mkmSlZpqmV+35R00zTTM1WTdxKE9xz31ABwQVFkJ05948ZubiDMAzg+3kensv8fueceecHV97OqrTWCCGEEEKInMnK0gEIIYQQQohHk2RNCCGEECIHk2RNCCGEECIHk2RNCCGEECIHk2RNCCGEECIHk2RNCCGEECIHk2RNiDxMKbVEKTXW9H1LpVSwhUMSQgiRQZKsCZELKKWClVJxSqmYNF9lnlRPa91Xa/15dsSYWUqpz5RSyywdx8OYEl2tlBpu4TieU0rl+s0xlVI2pufpaulYhMgNJFkTIvdor7UukOYrLLveWCllY8n6OUBPINL0v+IJ8sDPW4gcRZI1IXIxpZSVUmqdUipcKXVLKbVdKVUtzf2VSim/h9R7oGcjbdm7Q6ZKqbFKqXBgsel6B6XUYdN77VZKuT0irrvtD1BKnQVOmq7PUUqFKKVuK6UClVINTdfbAR8A3Uy9hvtN1wsrpZYqpa6Y6k1USj3w75ZSqpxSKlYpVSjNNS+l1FVTLFWVUjuVUlFKqetKqVUZeMYFgNeB/kB1pZTHffcbK6X2mtq+rJTqbrqeXyn1tVLqkuneTqWUneleI1OdW0qpQ0qpxmna262UmqSUCjLV+1EpVcR0e6epzN3eVS+lVBWlVIBS6obps6247zmEKKWGK6WOmtr77m4cpvuvm2K4rZQ6q5RqnZFnbyr7mVJqjantaOA/SqkGaT7jFaXULKWUbdrPARwzfY43TO2k6/dLiGeNJGtC5H6bgSpAKeBfYEUWtesCFADKAwOUUl4Yk7a+QDHAH9iolMr3mDY6AF6Au+n1P0BNoCiwDlirlLLTWm8GpgLfmnoN65rKrwTigMqAJ/Aq0Pv+N9FaXwaCMCZVd3UFvtdaJwOTgJ+BIqbPNTcDz+FN4KYp3q1Aj7s3lFIVTe1+hfGZ1AaOmm5/bfqs3qbPOxYwKKXKAZuACabrY4D1Sqliad6zh+mrDKBMbQE0Nn3eu72rgab7nwGlgepAJeDj+z7DW0Ar0726wN2EsiHGn+MIoDDQDLhoqpOuZ59GJ2AVUAhYAyQDQwFnoBHQFuiX9nMANUyf44en/P0S4pkgyZoQuccGU4/DLaXUBgCttUFrvUxrHa21jgf8gLpKKccseL9kwE9rnai1jgN8gHla60CtdYrW2t9UzusxbXyutb5pqo/WeoXWOtKUQE0FCgLPPayiUqos0AJ4X2sdq7UOB2YAXR7xXquAd0x1rYC3TdcAkgBXoLTWOl5r/Vd6HoBJT2C11tpgaq+b+v8w33+A37TW32utk7XW17XWh5RS1kAvYIjW+orpee3WWidhTMI2aa1/N/38fgMOY0xm7lqutT6utb4DjAe6KKXUw4LTWp/WWv9p+jldxZjYNbmv2AytdbjW+gbG5P5u7+C7wGJTfYPW+rLW+tRTPHuA3Vrrn0ztxJl+T/4xPZfzwKKHxJXW0/x+CfFMkGRNiNyjo9a6sOmrI4BSylopNVUpdV4pdRs4ayrrnAXvF6G1TkzzugIwOk3CeAtjb07Zx7RxOe0LpdQHSqmTSqkojL1Vjo+JtQJgB0Skeb+5QMlHlF8LvKSUKomxhyhea/236d4IwBYIMg0HpmvumTIOEzcGvjVd+hFjb+PdxKoccO4hVUsC+R5xrwLwzn3PsT7GXrS70j63ixifQ9FHxFhKKfW9UirU9DuwjAefaXia72NNn+Fx8Wf02d8fM0qpF5RSPyvjEP1tYOJD4rr/PTP6+yXEM0EmgQqRu/UAXgGaY/yjXgy4hnFo7JG01slKqQQgf5rLpfh/sgdw/6rDy8AnWusvMhBfahtKqWbAcIw9NsdNl6PSxPqw94sFipp6tR7/RlrfUEptwzhsWRv4Ls29KxiH1zDND9uilNqptb7whGZ7mOL7NU3HVj7T9c2mGGs+pF4EkIhxCPHYQz7XUq11/8e8b7k035cHEjAucCj8kLJfmO67a60jlVKdgemPafv+WCo/4nq6n73J/T+/hcBe4G2tdYxSaiTQ7hFl775nRn+/hHgmSM+aELmbE8Y/1DcwJl6TMlD3MMYhPWul1KvAi08ovwgYaJrUrpRSBZRS7TMw5OqEcWj1OsZeLj+MPWt3RQCud4f7TPPQdgDTlVIFlXExxXNpJ+M/xCqMw5av8/8hUJRSb5mG9gBuYUwWUtIRcw+Mw5Aeab7eBjqYJv2vBNoqpd5QxoUMzkqpWlrrFIw9XDNMPV/WyriowBbjnMJOSqlWpuv2Sqlm6t6tWHqYeqYcgU8wzr3TwFVAK6UqpSnrBNwBokzz4Uam43Pd9V+gr+n9rZRSLkqp55/y2d/PCWMyfkcZF73cna+G6fncwDiH7q7M/n4JkWdJsiZE7rYUCDN9HQP+fnzxewzBOCn8FsbeqE2PK6y1/gfjisj5GIcwT2Ocs5Vev2CcoH8GCAZuA1fS3F+DsdcqUim1z3TtPxgTuuOm91yLsQfwUTZgnGR/SWudtkfLGwhUSt0B1gMDtdaXAJRSp5RSb9/fkFLqRYxDk3NN873CTXO3fjTF/7apZ649MBpjz9cB/r+Y4n3gBLDfdO9zQGmtgzE+948x9oJewjhMm/bf4xUYE8ErgDUwDEBrHQ1MBv4xDRV6YlyoUA9jYrQJ+OExz+cepmHi94BZpvoB/L9XL6PP/n4jMCbO0Rh72dbcd38CsMr0OV7Pgt8vIfIsZfyPNSGEEDmBUmo3sERrvczSsQghcgbpWRNCCCGEyMEkWRNCCCGEyMFkGFQIIYQQIgeTnjUhhBBCiBwsz+yz5uzsrF1dXS0dhhBCCCHEE+3fv/+61rp4esrmmWTN1dWVoKAgS4chhBBCCPFESqmLTy5lJMOgQgghhBA5mCRrQgghhBA5mCRrQgghhBA5WJ6Zs/YwSUlJhISEEB8fb+lQhJnY29vj4uKCra2tpUMRQgghzCJPJ2shISE4OTnh6uqK6WxokYdorblx4wYhISFUrFjR0uEIIYQQZpGnh0Hj4+MpVqyYJGp5lFKKYsWKSc+pEEKIPM2syZpSqq1S6pRS6qxSasxjynVWSmmllGeaax+a6p1SSrXJRAxPW1XkAvLzFUIIkdeZbRhUKWUNzAVaASFAoFJqk9b6+H3lnIAhwD9prlUHugA1gDLAVqVUVa11irniFUIIIYTIiczZs1YPOKu1Pq+1TgRWA689pNynwFQg7VjWa8BqrXWC1voCcNbUXq4UHh5Oly5dqFy5MtWrV+eVV17h9OnTGW5nxowZxMbGZklMwcHBrFq1KvV1UFAQQ4YMyZK2hRBCCJF1zJmslQUup3kdYrqWSilVGyintd6c0bq5hdaaTp060bRpU86dO8fx48f5/PPPiYiIyHBbj0vWUlIy1ul4f7Lm6enJrFmzMhyTEEIIIczLnMnawyYT6dSbSlkBXwMjMlo3TRs+SqkgpVTQtWvXnjpQcwoICMDW1hZfX9/Uax4eHrz44ouMGjUKNzc33N3dWbNmDQDbt2+nadOmdO7cmRdeeIFu3bqhtWbWrFmEhYXRrFkzmjVrBkCBAgUYP3483t7e7Nmzh4kTJ+Ll5YWbmxs+Pj5obXxkZ8+epWXLltSqVYs6depw7tw5xowZw65du/Dw8ODrr79m+/bttGvXDoDIyEg6duxIzZo1qV+/PkeOHAHAz8+PPn360LRpUypVqiTJnRBCCJENzLl1RwhQLs1rFyAszWsnwA3YbpokXgrYpJTqkI66AGitFwGLADw9PR9I5tL65KdjHA+7nfFP8RjVyxRkQvsajy3z77//Urdu3Qeur1+/nkOHDnH48GGuX7+Ol5cXjRs3BuDgwYMcO3aMMmXK0KhRI/766y+GDBnCV199RUBAAM7OzgDcuXMHNzc3Jk6caIynenXGjx8PQPfu3dm8eTPt27enW7dujBkzhk6dOhEfH4/BYGDKlClMnz6dzZuNnZrbt29PjW3ChAnUrl2bDRs2sG3bNnr06MGhQ4cAOHnyJAEBAURHR/P888/Tv39/2eNMCCGEMCNz9qwFAlWUUhWVUvkwLhjYdPem1jpKa+2stXbVWrsCe4EOWusgU7kuSik7pVRFoAqwz4yxZrvdu3fzzjvvYG1tTcmSJWnSpAmBgYEA1KtXDxcXF6ysrPDw8CA4OPihbVhbW/PGG2+kvg4ICMDb2xt3d3e2bdvGsWPHiI6OJjQ0lE6dOgHGTWTz58//xNi6d+8OQPPmzblx4wZRUVEAvPrqq9jZ2eHs7EyJEiWeajhXCCGEEOlntp41rXWyUmoQ8DtgDfhrrY8ppSYCQVrrTY+pe0wp9T1wHEgGBmZ2JeiTesDMpUaNGqxbt+6B63eHKB/Gzs4u9Xtra2uSk5MfWs7e3h5ra2vAuKfcgAEDCAoKoly5cvj5+REfH//Y93mUh9W5u0VGemMTQgghRNYw6z5rWutftNZVtdaVtdaTTNfGPyxR01o3NfWq3X09yVTvea31r+aM05yaN29OQkICixcvTr0WGBhIkSJFWLNmDSkpKVy7do2dO3dSr97jF7w6OTkRHR390Ht3N4Z1dnYmJiYmNUEsWLAgLi4ubNiwAYCEhARiY2Mf21bjxo359ttvAePwqLOzMwULFszYBxdCCCFyo4RoCD9q6SjukaePm8oJlFL8+OOPDBs2jClTpmBvb4+rqyszZswgJiaGWrVqoZRi6tSplCpVipMnTz6yLR8fH15++WVKly5NQEDAPfcKFy7Me++9h7u7O66urnh5eaXeW7FiBf369WP8+PHY2tqydu1aatasiY2NDbVq1aJXr17Url07tbyfnx+9e/emZs2a5M+fn+XLl2f9gxFCCCFykugI2LcQApeAfWEYchCsrC0dFQDqaYbJciJPT08dFBR0z7UTJ05QrVo1C0Uksov8nIUQQjy162dhz2w49B2kJEK19tBoKLh4PrluJiil9mut0/Um0rMmhBBCiGdPSBDs/hpO/gzW+cCjKzQcDMUqWzqyB0iyJoQQQohng8EAZ/6Av2fBxb/AvhC8NAK8+0GBEpaO7pEkWRNCCCFE3pacCEfXGpO0ayehoAu0mQx1eoBdAUtH90SSrAkhhBAib4q/DfuXwd75EB0GJWpAp0Xg9jpY554N3SVZE0IIIUTeEh1uTNCC/CHhNri+BK/NhsotQD3sRMucTZI1IYQQQuQN104bhzqPrAFDMlR/DRoOgbJ1LB1ZpkiyZmYFChQgJiYm9fWyZcsICgpizpw5FoxKCCGEyEMu/QN/zYRTP4ONvXEuWoOBULSSpSPLEpKsCSGEECL3MRjg9G/GJO3yXnAoAk1Gg9d7UKC4paPLUmY9bko8Xq9eve45N7RAAeOKlO3bt9OkSRPeeustqlatypgxY/j222+pV68e7u7unDt3DoCffvoJb29vateuTcuWLVMPVffz86NPnz40bdqUSpUqMWvWrOz/cEIIIYQ5JCfAgRUwzxtWvwO3w+DlqfD+MWg2NtOJWty1cCKP7M+iYLPGs9Oz9uuYrD/rq5Q7vDzlsUXi4uLw8PBIfR0ZGUmHDh2e2PThw4c5ceIERYsWpVKlSvTt25d9+/Yxc+ZMZs+ezYwZM3jxxRfZu3cvSimWLFnC1KlT+fLLLwE4efIkAQEBREdH8/zzz9O/f39sbXPPyhchhBDiHvFRELTUuHAgJtz4N/iN/0L1jmCdNelM7JnTHO3dBUNKMnW37SWfQ/4saTeznp1kzUIcHBw4dOhQ6uu7c9aexMvLi9KlSwNQuXJlWrduDYC7u3vquaAhISG8/fbbXLlyhcTERCpWrJha/9VXX8XOzg47OztKlChBREQELi4uWfnRhBBCCPO7HWZa2bkUEqOhUlPoNB8qNcvSlZ0xe/dyboAPBpII/bgHDXJIogbPUrL2hB4wS7CxscFgMACgtSYxMTH1np2dXer3VlZWqa+trKxITk4GYPDgwQwfPpwOHTqwfft2/Pz8Hlrf2to6tY4QQgiRK1w9CX/PNq7s1ClQo5NxZWcZjyfXzaBbP24gdNxHXC1s4MLH3fBp+2GWv0dmPDvJWg7k6urK/v37eeutt9i4cSNJSUkZqh8VFUXZsmUBWL58uTlCFEIIIbKP1nBpr3HRwOlfwcYBPHsbV3YWcTXD22muz57N9XnzOVZBcX5UZ0a3+CjL3yezJFmzoPfee4/XXnuNevXq0aJFCxwdHTNU38/PjzfffJOyZctSv359Lly4YKZIhRBCCDMyGIzbbvw1E0ICwaEoNP3QuLLTsZh53jIxkStjP+L25s0E1FRc7PcynzfzQ+XATXOV1trSMWQJT09Pff9csBMnTlCtWjULRSSyi/ychRAil0qKhyOrjcOdN85C4QrQcDB4dIN85pszlnzzJiGDBhO3fz+rm1gT0bkRs5vPwTYbj6BSSu3XWnump6z0rAkhhBAie8XdNB4FtXcB3LkKpWtBZ3+o9lqWrex8lMTgYC7160diWBhzOtoS3bgWC5t+na2JWkZJsiaEEEKI7BEVYlzZuX8ZJMZA5ebQaChUbJItZ3bGBgURMnAQyRiY1M2WuOquLG0xh/y2OWfl58NIsiaEEEII84o4bjyz8+ha4yICt9eNKztL18y2EKJ+2syVsWPRpUvw0Wt3iC9VjG9aLqSQXaFsi+FpSbImhBBCiKynNVz8y7ho4MwfYJsfvPpC/QFQpEI2hqG5Pn8+12fNxrpOTUa2vkpMfgeWt15E8fy541gqSdaEEEIIkXUMKXByszFJC90P+Z2h2TjwehfyF83WUHRiIlfGTyBqwwbsXm3NMO8zXE9JYFmrZZRzKpetsWSGJGtCCCGEyLykODj8nXFlZ+R5KFIRXv0KPLqCrUO2h5MSFUXI4CHE7tuHk29f3q+0h/Do6yxqtYiqRapmezyZIQe5m9GNGzfw8PDAw8ODUqVKUbZs2dTXaU8ryKytW7dSqFCh1LY9PDxSj6TKiOTkZAoXLvzUcfj7+xMeHp76unfv3pw6deqp2xNCCJELxEbCzmkwwx02vw/2heDN5TB4v7E3zQKJWuLlywS/05XYgwdx/nwiHz1/hLNR55jRdAYeJbL+BARzk541MypWrFjquaB+fn4UKFCAkSNH3lNGa43WGiurzOXNzZo1Y8OGDZlqI7P8/f2pU6cOpUqVAmDp0qUWjUcIIYQZ3boMe+fB/uWQdAeea2Vc2en6Yras7HyU2IMHCRk4CJ2SQtklC/kwZiUHQw4ytclUGpZtaLG4MkN61izg7NmzuLm54evrS506dbh8+fI9PVqrV6+mb9++AERERPD666/j6elJvXr12Lt3b7rfZ8SIESxatCj19bhx45g5cya3b9+mefPm1KlTh5o1a7J58+YH6m7dupWOHTumvvb19WXlypUATJgwAS8vr9TPoLVmzZo1HDp0iLfffju15/DFF19MTVZXrlyJu7s7bm5ujB07Fvh/T96YMWOoVasWDRo04OrVqxl4kkIIIbJd+L+w3gdm1oJ9i6Bae/D9C/6zDiq+ZNFE7fZvv3GpZy+sHB2p8N23TErYwM6QnXzc4GPaura1WFyZ9cz0rH2x7wtORp7M0jZfKPoCo+uNfqq6x48fZ+nSpSxYsOCxh6wPGTKEDz74gPr16xMcHEy7du34999/HygXEBCAh8f/u3Y3bNhAly5dGDNmDD4+PgCsXbuWgIAAHBwc2LhxI05OTly9epVGjRrRrl27dMc+dOhQPvnkE7TWdO3ald9++423336b2bNnM2fOnHviAAgJCWHcuHEEBQVRqFAhWrZsyebNm2nbti1RUVE0adKEKVOmMHz4cPz9/RkzZky6YxFCCJENtIbgXcZFA2e3gq0jePtC/f5Q2PIT9bXW3FiyhGtffoVD7dqUnTObaWcX8suFXxhaZyhvVn3T0iFmyjOTrOU0lStXxsvL64nltm7des+8r5s3bxIXF4eDw71zAB42DOrq6srly5eJiIggJCSEUqVKUaZMGRITExk9ejS7d+/GysqKy5cvc/369XTPV/vzzz+ZNm0a8fHxXL9+nbp16/Lyyy8/svw///xD8+bNcXZ2BqBr167s3LmTtm3b4uDgkFq3bt267Nq1K10xCCGEyAaGFDixyZikhR0Ex+LQ/GPjXDSHIpaODgCdlET4xIncWruOgq+8TOnJk5l/YgnfnfyOXjV68a7bu5YOMdOemWTtaXvAzCXtoe1WVlakPaM1Pj4+9XutNfv27SNfvnxP9T5vvPEGP/zwA8HBwXTp0gWAb775hqioKA4cOICNjQ0uLi73vCeAjY0NBoPhgZhiY2MZNGgQBw4coGzZsowbN+6Buvd73PmzaT+XtbX1Y3sZhRBCZJPEWDj0LeyZAzeDoWhlaDcDar0DtvaWji5VSnQ0oUOHcefvvynWrx/Fhw5h1anvWHB4AZ2e68TwusNz5MHsGSVz1nIAKysrihQpwpkzZzAYDPz444+p91q2bMncuXNTX9+dA5ZeXbp0YfXq1axfv5433ngDgKioKEqUKIGNjQ1btmwhNDT0gXoVKlTg2LFjJCYmcvPmTbZt2wZAXFwcVlZWODs7Ex0dzQ8//JBax8nJiejo6Afaql+/PgEBAdy4cYPk5GRWr15NkyZNMvQ5hBBCZIPYSNj+Bcxwg19GQv5i8NYKGBQInr1zVKKWFBrKxa5dubNvH6UnfUaJ94ex+cLPTNk3hRblWzC+wfg8kajBM9SzltN98cUXtG3blvLly1O9enUSEhIAmDt3Lv3792fp0qUkJyfTrFmze5K3u+6fszZhwgQ6depErVq1uHbtGhUrVqREiRIAdO/enfbt2+Pp6UmdOnWoUqXKA+1VrFiRjh074u7uTtWqValTpw5gXOHas2dP3NzcqFChAt7e3ql1evfuTd++fXFwcGDfvn2p111cXJg4cSJNmzZFa0379u159dVXpRdNCCFyipsXYc9cOLgCkmKhShvjys4KDS26YOBR4o4e5XL/AeiEBMovXoRjgwYEXArg478+xru0N180/gIbq7yT4qjHDVHlJp6enjooKOieaydOnKBatWoWikhkF/k5CyHEU7pyGP6aBcd+NCZl7m9Bw8FQsrqlI3uk21u2EDbqA2yKFaPcwgXYPfccgeGB+G7xpWqRqixpswRHW8cnN2RhSqn9WmvP9JTNO2mnEEIIIZ5Mazi/3bho4HwA5HMyruqsPwAKlbV0dI+ktSZy2XKuTp2Kvbs75ebNxcbZmeM3jjN422BcnFyY33J+rkjUMkqSNSGEEOJZceUw/DTUuLKzQElo6Qd1e4PD059ekx10cjLhkyZx67vVOLVqRZmpX2Dl4MCFqAv039qfgvkKsrDVQgrb5+zP8bQkWRNCCCHyupRk+GsGbJ9sPFi9/Syo+XaOWjDwKCkxdwgd/j53du6i6Lt9KDFiBMrKivA74fhsMe4juqjVIko5lrJwpOYjyZoQQgiRl0Weh/X9IGQf1OhkPFw9f1FLR5UuSeHhXO7nS8LZs5Ty86NIl7cBuBl/E58tPsQkxuDfxh/XQq6WDdTMzLp1h1KqrVLqlFLqrFLqgW3plVK+SqmjSqlDSqndSqnqpuuuSqk40/VDSqkF5oxTCCGEyHO0hv3LYP6LcP0UvL4EOi/NNYla/PHjBL/1NkkhIZRbsCA1UYtJjKH/1v6ExYQxp8UcqhXL+wvMzNazppSyBuYCrYAQIFAptUlrfTxNsVVa6wWm8h2Ar4C7h3ed01rfe26REEIIIZ4sOgJ+GgKnf4OKTaDj/By9eOB+0QEBhI4YiXWhQlRYtQr756sCkJCSwNCAoZyMPMnMZjOpW7KuhSPNHubsWasHnNVan9daJwKrgdfSFtBa307z0hHIG/uIpBEeHk6XLl2oXLky1atX55VXXuH06dMEBwfj5uaWZe8zfvx4tm7dCsCuXbuoUaMGHh4ehIaG0rlz56dqc9myZYSFhaW+7tu3L8ePH39MDSGEEBZ3fBPMq29c8dn2C+i+IVclapErVhIycBB2rq64rl6dmqglG5IZtWMU+8L38dmLn9Gk3LOzubo556yVBS6neR0CeN9fSCk1EBgO5AOap7lVUSl1ELgNjNNaP3BopFLKB/ABKF++fNZFnkW01nTq1ImePXuyevVqwHgCQUREBOXKZe3BtxMnTkz9/ttvv2XkyJH07t0bgHXr1j1Vm8uWLcPNzY0yZcoAsGTJkswHep/k5GRsbGTqpBBCZFp8FPw6Bg6vgtK14PXFUPx5S0eVbjolhYgvvuDmNyso0Lw5ZadPwyp/fgAM2oDf334EXA7gw3of0q5SOwtHm73M2bP2sC2PH+g501rP1VpXBkYD40yXrwDltda1MSZyq5RSBR9Sd5HW2lNr7Vm8ePEsDD1rBAQEYGtri6+vb+o1Dw8PXnrppXvKBQcH89JLL1GnTh3q1KnD33//DcCVK1do3LgxHh4euLm5sWvXLlJSUujVqxdubm64u7vz9ddfA9CrVy/WrVvHkiVL+P7775k4cSLdunW7pwcvJSWFkSNH4u7uTs2aNZk9ezZgTPS8vLxwc3PDx8cHrTXr1q0jKCiIbt264eHhQVxcHE2bNuXuxsPfffcd7u7uuLm5MXr0/89dLVCgAB999BG1atWifv36REREPPBc/Pz88PHxoXXr1vTo0YOUlBRGjRqFl5cXNWvWZOHChallp02blnp9woQJWfFjEUKIvCd4t3Fu2pHV0HgUvLs1VyVqhthYQgYP4eY3Kyjaswcus2elJmpaa6YHTWfjuY0MqDWArtW6Wjja7GfOLo0QIG33kQsQ9oiyYBwmnQ+gtU4AEkzf71dKnQOqAkGPrv544Z9/TsKJk09b/aHsqr1AqbFjH3n/33//pW7dJ4+nlyhRgi1btmBvb8+ZM2d45513CAoKYtWqVbRp04aPPvqIlJQUYmNjOXToEKGhofz7778A3LppGTRDAAAgAElEQVR16562+vbty+7du2nXrh2dO3cmODg49d6iRYu4cOECBw8exMbGhsjISAAGDRrE+PHjAeNRVJs3b6Zz587MmTOH6dOn4+l57wbLYWFhjB49mv3791OkSBFat27Nhg0b6NixI3fu3KF+/fpMmjSJDz74gMWLFzNu3Djut3//fnbv3o2DgwOLFi2iUKFCBAYGkpCQQKNGjWjdujVnzpzhzJkz7Nu3D601HTp0YOfOnTRu3PiJz1QIIZ4JSfEQ8Bn8PQeKVoQ+f0A5L0tHlSFJEVcJ6d+f+JMnKTluHEX/0+2e+4uPLmbF8RV0q9YN31q+j2glbzNnshYIVFFKVQRCgS7APemwUqqK1vqM6eWrwBnT9eJApNY6RSlVCagCnDdjrBaVlJTEoEGDOHToENbW1pw+fRoALy8v+vTpQ1JSEh07dsTDw4NKlSpx/vx5Bg8ezKuvvkrr1q3T/T5bt27F19c3ddixaFHjiqCAgACmTp1KbGwskZGR1KhRg/bt2z+yncDAQJo2bcrd3sxu3bqxc+dOOnbsSL58+WjXztg9XbduXbZs2fLQNjp06ICDgwMAf/zxB0eOHEkdro2KiuLMmTP88ccf/PHHH9SuXRuAmJgYzpw5I8maEEIAhB+F9T5w9Th49oHWn0G+3LV7f/ypU1z27U9KVBQuc+fg1KzZPffXnFzD7IOzaVepHR94fZBnDmbPKLMla1rrZKXUIOB3wBrw11ofU0pNBIK01puAQUqplkAScBPoaareGJiolEoGUgBfrXVkZuJ5XA+YudSoUSNd88W+/vprSpYsyeHDhzEYDNjbGzcpbNy4MTt37uTnn3+me/fujBo1ih49enD48GF+//135s6dy/fff4+/v3+64tFaP/CLHh8fz4ABAwgKCqJcuXL4+fkRHx//xHYexdbWNvU9rK2tH3lYu6Pj//9B0Voze/Zs2rRpc0+Z33//nQ8//JB+/fo9Nh4hhHimGFLg71mwbZJxG46ua6Fq+v/DPaeI2bWb0GHDsHJ0xHXlCuyr33se6a8XfmXSP5No6tKUiY0mYqXMuttYjmbWT661/kVrXVVrXVlrPcl0bbwpUUNrPVRrXUNr7aG1bqa1Pma6/oPpei2tdR2t9U/mjNNcmjdvTkJCAosXL069FhgYyI4dO+4pFxUVRenSpbGysmLFihWkpKQAcPHiRUqUKMF7773Hu+++y4EDB7h+/ToGg4E33niDTz/9lAMHDqQ7ntatW7NgwYLUBCoyMjI1MXN2diYmJuae5NLJyYno6OgH2vH29mbHjh1cv36dlJQUvvvuO5o0efpVOW3atGH+/PkkJSUBcPr0ae7cuUObNm3w9/cnJiYGgNDQUK5evfrU7yOEELnezWBY9ips9YPnX4b+e3JlonZz9Rou+/piW64crt+veSBR2x26m7G7xlKnZB2mNZmGrZWthSLNGWQZnhkppfjxxx8ZNmwYU6ZMwd7eHldXV2bMmHFPuQEDBvDGG2+wdu1amjVrltrrtH37dqZNm4atrS0FChTgm2++ITQ0lN69e2MwGACYPHlyuuPp27cvp0+fpmbNmtja2vLee+8xaNAg3nvvPdzd3XF1dcXL6/9zHXr16oWvry8ODg7s2bMn9Xrp0qWZPHkyzZo1Q2vNK6+8wmuvvfawt0x3XMHBwdSpUwetNcWLF2fDhg20bt2aEydO0KBBA8C4eGHlypWUKFHiqd9LCCFyJa3h4Ar47UNQVtBpofG4qFw2LKgNBq5O/5JIf38cmzSm7JdfYV3g3qHbg1cP8n7A+1QpUoXZzWdjb5Pzj8QyN/W4Ia3cxNPTU99dqXjXiRMnqFYt7+9s/KyTn7MQIk+LuWbc4PbUL+D6knGD28JZu/1TdjDExRH2wWiit2yhSNd3KDl2LOq+rZtORZ6i9++9KWZfjGVtl1HMoZiFojU/pdR+rbXnk0tKz5oQQgiRc538GTYNgYRoaPM5ePcHq9w3dyv5+nUuDxhI/NGjlBgzmqI9ez4wh/rS7Uv029KP/Db5WdhqYZ5O1DJKkjUhhBAip0mIht/GwMGVUModXt8MJXLnCELC2bNc9ulHcmQkLrNn4dSy5QNlrsZexWeLDyk6Bf9W/pQpUMYCkeZceT5Ze9gKSJF35JVhfCGESHVxD/zYD6Iuw4vDoemHYJPP0lE9lTt79hAyZCjKzo4KK1bg4P7gMYtRCVH029KPm/E3+W+b/1KpcCULRJqz5b6+1Aywt7fnxo0b8gc9j9Jac+PGjdStToQQIldLToAt42Hpy8aFA71/hZYTcm2iduuHH7j0ng+2pUpScc3qhyZqsUmxDPhzABdvX2RW81m4OWfdmdl5SZ7uWXNxcSEkJIRr165ZOhRhJvb29ri4uFg6DCGEyJyIY8YNbiP+hTo9jfPT7ApYOqqnog0Grs2YyY1Fi3Bs2JCyM2dg7eT0QLnElESGBQzj3+v/8lXTr/Au/cDx4cIkTydrtra2VKxY0dJhCCGEEA9nSIE9c2Hbp2BfCN5ZA8+3tXRUT82QkMCVDz/k9i+/UvjNNyk1/mOU7YN7pKUYUvhw14fsubKHiQ0n0qJ8CwtEm3vk6WRNCCGEyLFuXYIf+8PF3fBCO2g/ExydLR3VU0uOjCRk4CDiDh6kxMgRFH333YfOGdda8+neT/nj4h+M9BxJpyqdLBBt7iLJmhBCCJGdtIZDq+DX0cbXr80Dj665boPbtBLOX+Cyry/JERGUnTGDgm3bPLLszAMz+eHMD7zn/h49a/R8ZDlL2Xn6Ghdv3KF7A1dLh5JKkjUhhBAiu9y5Dj8NhZOboUIj4wa3RSpYOqpMubNvHyGDh6CsramwfBkOHh6PLLv036X899//8mbVNxlce3A2Rpk+Gw+FMuL7wzxfyom3vcqTzyZnrMOUZE0IIYTIDqd+g02DIf4WtPoUGgwEK2tLR5UpURs3EjbuY/KVK0e5hQvIV+7RJyusP7Oer/Z/RVvXtnzk/VGO21bLf/cFJm4+jnfFoizu6ZljEjWQZE0IIYQwr4QY+H0sHFgOJd2g+49QKndvUaG15vqcuVyfO5f83t64zJqJdaFCjyy/9eJWPtnzCY3KNuLzFz/HOgclqVprpv1+innbz9G2RilmdPHA3jbnxAeSrAkhhBDmc+kf+NEHbl6ERsOg2ViwsbN0VJliSEzkyrhx3N70E4U6daL0J36ofI/eC25P2B4+2PkBNZ1r8lWTr7C1fnB1qKUkpxgY++NRvg8Koat3eT59zQ1rq5zV4weSrAkhhBBZLzkRdkyB3V9DIRfo/QtUaGjpqDIt5dYtQgYNJjYoiOLDhlKsX7/HDmceuXaEoQFDcS3kypwWc8hvmz8bo328uMQUBn93gK0nrjKkRRXeb1klxw3N3iXJmhBCCJGVrp4wbnAbfgRq/wfaTAb7gpaOKtMSL17kcj9fkkJDKTN9OoXavfrY8mdvnmXAnwMoZl+MhS0XUsju0cOk2S0qNol3lwey/9JNPn2tRo5a+fkwkqwJIYQQWcFggH/mw9ZPwM4JuqyCFx6f0OQWsQcOEDJgIADlly0lf926jy0fGhNKvy39yGeVj0WtF1E8f/HsCDNdwqPi6eH/D8HXY5nbtQ6vuJe2dEhPJMmaEEIIkVm3LsOG/hC8C55/BdrPggI5J0HJjKiff+bKh2OxLV3auOLT1fWx5a/HXcfnDx/iUuJY1nYZ5ZwevUI0u529GkNP/31ExSWxrLcXDZ/LHZsQS7ImhBBCPC2t4cj38MtI0AboMBtqd8/VG9zepbXmxsKFXJsxEwfPurjMno1NkSKPrXM78Ta+W3y5FneNRa0WUbVI1WyK9skOXb5F76X7sLayYrVPfdzK5pxh2SeRZE0IIYR4GrGRsHkYHN8I5epDpwVQNG+cR60TE7ni9wlR69dTsH17Sk/6DKvHrPgEiEuOY/CfgzkXdY45zefgUeLRm+Nmt+2nrtJ/5QGKO9mx4t16VCjmaOmQMkSSNSGEECKjzmyBjQONCVtLP2g4JNdvcHtXyu3bhAwZSuzevTgPHIjzoIFPXCWZZEhi5I6RHLx6kKlNptKobKNsivbJNhwMZeTaw1Qt6cSyPl6UcLK3dEgZJsmaEEIIkV6Jd+CPcRDkDyWqQ7d1ULqmpaPKMokhIVzu50vipUuUnjKZwh07PrGOQRsYt3scO0N2Mr7BeNq6ts2GSNNnya7zfPbzCRpUKsaiHnVxss85e7xlhCRrQgghRHqEBBm35Ig8Dw0GQfOPwTb39dI8Stzhw1weMBCdlET5JUtw9K73xDpaayb/M5lfLvzC0DpDebPqm9kQ6ZNprZny20kW7jjPK+6l+OqtnHcqQUZIsiaEEEI8TkoS7JgKu76EgmWg12ZwfdHSUWWp27//QdgHH2BTogTlFi7ArlKldNWbd3geq0+tpleNXrzr9q6Zo0yf5BQDY9YfZd3+EP5TvzyfdMiZpxJkhCRrQgghxKNcO2XsTbtyCGp1hZengH3uWUX4JFprIv39uTptOg4eHrjMm4tN0aLpqrvy+EoWHF5Ap+c6Mbzu8Byx+39cYgqDVh3gz5NXGdayCkNb5NxTCTJCkjUhhBDifgYD7FsEWyeAbX54awVU72DpqLKUTkoi/NPPuPX99zi93JYykydjZZ++Yd1N5zbxReAXtCjfgvENxueIhOhWbCLvLg/iwKWbfNbRjf/Ur2DpkLKMJGtCCCFEWlGhsHEAnN8OVVpDhzngVNLSUWWplOhoQoe9z52//qKYjw/Fhw1FWVmlq27ApQDG/zUe71LefNH4C2ysLJ9KXImKo8d/93HxRizzutbh5VxwKkFGWP4JCyGEEDnF0XXw83BISYZ2M6BurzyxwW1aSWFhXO7nS8KFC5T+7FMKd+6c7rqB4YGM3DGSakWrMbP5TOys7cwYafqcvRpNj//uIzo+meV96tGgcjFLh5TlJFkTQgghYiPh5xFwbD241DNucFussqWjynKxBw8SMmQIOi6e8osW4tiwYbrrHr9xnMHbBuPi5ML8lvNxtLX8xrIHLt2kz7JAbKysWN2vPjXK5J35hGlJsiaEEOLZdnYrbBwEd64Zt+NoNAys896fx5ur1xA+aRK2pUpRzt8fuypV0l33QtQFfLf4UjBfQRa2Wkhh+8JmjDR9Ak5dpf/K/ZQsaM+KPt6UL5bf0iGZTd77bRRCCCHSIzEWtoyHwMVQ/AV4ZzWUyTlHJGUVQ2IiEZ9+yq2163B86SXKTpuKdeH0J1vhd8Lx2eKDUopFrRZRyrGUGaNNn/UHQvhg3RGeL+XEst71KO5k+eFYc5JkTQghxLMndL9xS44bZ6H+AGgxHmwdLB1VlksKDydk6FDiDx+hWL9+FB8yGGWd/s1hb8bfxGeLDzGJMfi38ce1kKv5gk2nxTvPM+mXEzSsXIyF3XPvqQQZIcmaEEKIZ0dKknFz2x1TwakU9NgIlZpaOiqziA0KImToMHRcHGVnzaRg69YZqh+TGEP/rf0JiwljYauFVCtWzUyRpo/BYDyVYNHO87zqXpqv3q6FnU3uPZUgIyRZE0II8Wy4fsbYmxZ2AGq+DS9PBQfLz73Kalprbn67iogpU8hXtiwuy5dh99xzGWojISWBoQFDORl5kpnNZlK3ZF0zRZs+SSkGRv9whPUHQunRoAIT2tfI9acSZIQka0IIIfI2rSFwCfxhOsvzzWVQo5OlozILQ3w84RP8iNq4kQLNmlFm6hdYOzllqI1kQzKjdoxiX/g+Jr80mSblmpgp2vSJTUxm4LcHCDh1jRGtqjKo+XM5YhPe7JS+HfCeklKqrVLqlFLqrFJqzEPu+yqljiqlDimldiulqqe596Gp3imlVBtzximEECKPun0FVr4Bv4yECg2h/548m6glhYVxsWs3ojZuxHnQIFzmzslwombQBvz+9iPgcgBj6o2hXaV2Zoo2fW7eSaTbkn/Ycfoan3dyZ3AeOT4qo8zWs6aUsgbmAq2AECBQKbVJa308TbFVWusFpvIdgK+AtqakrQtQAygDbFVKVdVap5grXiGEEHnMvz/A5uGQkgivfgme7+a5DW7vurP3H0Lffx+dlITLvLk4NW+e4Ta01kwPms7GcxsZUGsA3ap1M0Ok6Rd2K44e/vu4FBnLvG51aOuWt04lyAhzDoPWA85qrc8DKKVWA68Bqcma1vp2mvKOgDZ9/xqwWmudAFxQSp01tbfHjPEKIYTIC+Juwi+j4OhaKFsXOi0C54zN2cottNZELlvO1enTyefqisvs2dhVqvhUbS0+upgVx1fQ9YWu+NbyzeJIM+ZMRDQ9/PcRE5/MN33qUb9S3juVICPMmayVBS6neR0CeN9fSCk1EBgO5APu/qdAWWDvfXXLPqSuD+ADUL58+SwJWgghRC6lNZz5Aza/D9Hh0OwjeHF4ntzgFsAQF8eVcR9z++efcWrVitKTJ2Nd4OlOFVhzcg2zD86mXaV2jK432qJDjfsvGk8lyGdjxZp+DahepqDFYskpzPkb/LCftH7ggtZzgblKqa7AOKBnBuouAhYBeHp6PnBfCCHEMyLiGPz+EZwPAOeq0HeLsVctj0oMCSFk0GASTp2i+LBhFOvn89QJ1q8XfmXSP5No4tKEiY0mYqXMOp39sbadjGDAtwcoVdCeFe96U65o3j2VICPMmayFAOXSvHYBwh5TfjUw/ynrCiGEeBbFXIWASXDgG7ArCG2nGOem2eSzdGRmE7P7L8JGjEBrTbmFCyjQuPFTt7UrZBdjd42lTsk6TG8yHVsry20wu25/CKN/OEL10gVZ2tsL5wJ5+1SCjDBnshYIVFFKVQRCMS4Y6Jq2gFKqitb6jOnlq8Dd7zcBq5RSX2FcYFAF2GfGWIUQQuQmSfGwdx7s+gqS46BeP2jyAeQvaunIzEZrzY0lS7j29QzsnnsOlzmzyZeJKUAHrx5k+PbhVClShdnNZ2NvY5+F0WbMwh3nmPzrSV58zpkF3etSwC5vDl0/LbM9Da11slJqEPA7YA34a62PKaUmAkFa603AIKVUSyAJuIlxCBRTue8xLkZIBgbKSlAhhBBoDcfWwxY/iLoEz78CrT7NswsI7jLcuUPY2I+I/v13nF5uS5lJk7DK//RDhKciTzHwz4GUdCzJ/JbzccqXsS0+sorBoJn86wkW77pAu5ql+fKtZ+dUgoxQWueNqV6enp46KCjI0mEIIYQwl5Ag+O1DCNkHJd2hzSSoZNkNW7ND4sWLhAwaRMK585QYMZyiffpkagHApduX6PFrD6ytrFnx8grKFCiThdGmX1KKgQ/WHeHHg6H0aujK+HbVsXqGTiVQSu3XWnump6z0MwohhMjZbl2GPz8xbsVRoCR0mAMeXcEq7/fARG/fTtioD1BWVpRfshjHhg0z1V5IdAg+W3xI0Sn4t/K3WKIWm5hM/5UH2HH6GqPaPM+AppWfyc1u00uSNSGEEDlTQjTsngF75hhfvzQSXhwGdpYZsstO2mDgxsKFXJs1G7sXXsBl9mzyuTywg1WGnIo8he9WXxJTElnYaiGVClfKomgzJvJOIr2XBXI05BZTXnenSz3ZeutJJFkTQgiRsxhS4NC3sO0ziIkA97egxXgoXO7JdfOAlJgYwsaMIWbrnxRs357SEz/BysEhU20GhgcyZNsQHG0d+eblb6hcuHIWRZsxITdj6eG/j5Cbccz/T13a1ChlkThyG0nWhBBC5Bzndxj3S4s4CuW8ocsqcEnXtJ48IeH8eUIGDSbx4kVKjv2QIt27Z3p4cMvFLYzeOZryTuVZ0GoBpRwtkyCdjoimx3/3cScxmRV96uH9jJ9KkBGSrAkhhLC862fhj3Fw+lcoXB46LzUeuP4MzWOK/vNPwj4YjbKzo7y/P47e9TLd5pqTa5j0zyRqFa/FnBZzKGRXKAsizbig4Ej6LAvE3taa7/s1oFppOZUgIyRZE0IIYTmxkbBjKgQuBhsHaOkH3v3B1nJ7fmU3bTBwfc4crs+bj72bGy6zZ2FbOnOHlmutmXd4HgsOL6CJSxOmNZmGg03mhlKf1tbjEQxcdYCyhR1Y3qeenErwFCRZE0IIkf2SEyFwCez4AhJuQ52e0GwsFChh6ciyVcrt24SN+oCYHTso9PrrlJowHiu7zO3cn2xIZtI/k1h3eh0dn+vIhAYTsLGyzJ/774Mu8+H6o9QoU5ClvbwoJqcSPBVJ1oQQQmQfreHUr8Yhz8hzUKmZcb+0kjUsHVm2SzhzhsuDBpEUGkbJ8R9T5J13Mj0/LT45ntE7R7Pt8jb6uvdlSO0hFtkSQ2vN/B3nmPrbKV6q4syC/9TFUU4leGry5IQQQmSPK0fg97EQvMt42HrXtVCl1TM1L+2u27/9TtjYsVjlz0+F5cvIXzfzh87fTrzN4D8Hc/DqQcbUG0O3at2yINKMMxg0n/18Av+/LtChVhmmv1mLfDaWOxw+L5BkTQghhHlFh8O2T+Hgt+BQBF6ZDnV7gbXlDg23FJ2SwrUZM7mxeDEOtWpRdtZMbEuWzHS7V2Ov0m9LP4JvBzO18VTaVmybBdFmXGKygVHrDrPxUNgzeSqBuUiyJoQQwjyS4uDvObD7a0hJhAYDofEocChs6cgsIuXWLUJHjOTOX39R+K23KDnuI6zy5ct0uxeiLuC7xZdbCbeY12IeDco0yIJoM+5OQjL9vz3AztPX+KDt8/RvIqcSZBVJ1oQQQmQtgwH+XQdbP4HbIVCtPbSaCEUts2N+ThB/8iQhgwaTHBFBqYmfUOStt7Kk3SPXjjDwz4FYKSuWtl1K9WLVs6TdjLoRk0CfZYEcDY1i6hs1ecvr2djAOLtIsiaEECLrXNprnJcWuh9K14LXF4FrI0tHZVFRm3/myrhxWBcsSIUV3+Dg4ZEl7e4K2cWIHSMoZl+MRa0WUa6gZRKkkJux9PjvPkJvxbGwuyetqmd+WFfcS5I1IYQQmXczGLb6wbEfwak0dFwANd8Gq2d3YrlOTubq9C+JXLYMh7p1cZnxNTbFi2dJ25vObWLCXxOoUqQK81rOw9nBOUvazaiT4bfp6b+PuMQUVvb1xsu1qEXiyOskWRNCCPH04m/Dri9h7zxQ1tBkDDQaAvkcLR2ZRSVHRhI6fASxe/dSpGtXSo4ZjcqC+WkAy/5dxpf7v8S7lDczms2gQL4CWdJuRgUGR/LuskAc8lmz1rchz5dyskgczwJJ1oQQQmRcSjIc/Aa2TYLY61CrK7T4GAqWsXRkFhd37BghgweTcv0GpT//nMKvd8qSdg3awJdBX/LN8W9o49qGz1/8nHzWWZMAZtSW4xEMWnWAskUc+KZPPVyKyKkE5iTJmhBCiIw5+6fxsPVrJ6B8Q2i7DsrUtnRUOcKtDRsIn+CHddGiVPj2Wxzc3bKk3aSUJD7++2N+Pv8z77zwDmPqjcFKWWaIeU3gJT5cfxR3l8Is7eVFUUfLJIzPEknWhBBCpM/Vk8aTB85ugSKu8NYK40pP2Z4BnZRExBdTublyJfnr1aPs119hU6xYlrQdmxTL+9vf5++wvxlSewh93fta7FSCedvPMe33UzSuWpz53erIqQTZRJ6yEEKIx7tzA7Z/DkFLIV8BaP0Z1PMBGznnESD5+nVCh71PbFAQRXv2pMSokSibrPnzGhkfycCtAzkeeZxPGn7C61Vez5J2M8pg0EzcfJxlfwfT0aMMUzvLqQTZSZI1IYQQD5ecAP8shJ3TITEGPPtA0zHgaJmVhzlR3JEjhAweQkpUFGWmTaVQ+/ZZ1nZoTCj9tvQj/E44M5rOoFn5ZlnWdkYkJhsYsfYwPx0O490XK/LRK9XkVIJsJsmaEEKIe2kNJzbBlvHGLTmqtIZWn0KJFywdWY5ya906wj+ZiE2JErh+twr7atWyrO1Tkafov7U/CSkJLG69mNolLDMnMCYhmf4r97PrzHVGt30B3yaV5FQCC5BkTQghxP+FHjAuHrj0NxSvBv9ZD8+1sHRUOYpOTCT888+5tXoNjg0bUObLL7EpUiTL2g8MD2TItiE42jqyvO1ynivyXJa1nRE3YhLovSyQY2G3mdq5Jm95yqkEliLJmhBCCIgKNR62fvg7yO8M7b6G2j3AWv5MpJUUcZXQYcOIO3iQYn3fpfiwYVk2Pw1g68WtjN45GhcnFxa2Wkgpx1JZ1nZGXI6MpYf/Pq5ExbGoe11aVJNTCSxJ/l8ohBDPssQ78Ncs+Gsm6BRoNAxeGgH2BS0dWY4Te+AgIUOHYIi5Q9mvv6Lgyy9nafvfn/qez/Z+Rs3iNZnbYi6F7AplafvpdeKK8VSChGQD3/b1pm4FOZXA0iRZE0KIZ5HBAEdWw58TIfoK1OgELf2MW3KIe2itubVmDeGTPse2dGnKL/kv9s9XzdL25x+ez/zD82ns0pjpTabjYOOQZe1nxD/nb9D3myAc89mw1rcBVUvKqQQ5gSRrQgjxrAnebTxs/cphKFsX3lwO5b0tHVWOZEhIIHziRKJ+WI9j45coO20a1oX+x959h0dVbQ8f/570DgmhJ6H3DqFLr6KAigUQLBeVIh2kCITeO0pHFALSFQEFKdIEERJ67yWNAOk9mdnvHwP3zeWHGGXOTBLW53l8ksnMOWvdSzJZ2efstcy34mUwGpj05yQ2X91Mx1IdGdtgLPY29mY7/z/x64UI+q07ha+nM6t71KVoXusUjOL/kmJNCCFeFo9umHZ4Xt4BHj7w1gqo3OmlHrb+POnh4YT0H0DKuXPk692L/H37otnamu38qYZUhh8azr67+/ikyif0r9Hfajst1x+/y5c/nqOqT15WylSCbEeKNSGEyO2SY+DQTFPPNFsHaDYa6n8ODjLP8a8knThByMBBqORkin61AI9Wrcx6/ri0OPr/1p/g+8EMrz2cbhW7mfX8WaWUYuH+68zafZWm5fKz6P2auDhIaZDdyL+IEELkVoZ009SBA1MhORpqdIPmo8HdOjsMcwKlFNGBa7g/YwYOvoEdh7wAACAASURBVL74rF6FY6lSZo0RmRRJr729uBV7ixmNZ/BqCfNuVMgqo1ExfvsFVv1xhzdrFGXG21Wxt5VV1uxIijUhhMhtlIJru01zPB9eheKNoM0UKFzV2plla8aUFCLGjiX2p224NW9OkenTsHU37w32t2Jv0WtPL2JSY1jUYhH1i9Q36/mzKjXDwOCNZ/j5bDifNirByFdlKkF2JsWaEELkJvcvmJra3twPXqWg8zoo96oMW/8b6aGh3OvXj9SLl/Du1xfv3r3RzHwv37kH5+izrw82mg0r266kUr5KZj1/ViWkZtAzMIgj1x8x8tXy9Gxi3pVDYX5SrAkhRG6QEAn7J8PJ1eDoAW2ngX8PsJMbxf9O4rFjhA4chMrIwGfxItybmX8G5++hvzP4wGDyOeVjaaul+Hn4mT1GVtyPS6HHqhNcCo9n1jvVeLuWj1XyEP+MFGtCCJGTpafAsUVweA5kJEOdntBkGLhII9O/o5Qi6tvviJw1C4eSJfD56iscS5Qwe5ztN7YTcCSA0p6lWdxyMd7O3maPkRWn78Xw2eogElIzWP5BLZqXl6kEOYUUa0IIkRMpBRd+gD3jIPYulGsHrSaAdxlrZ5YjGJOSCB89hrhffsG9VSsKT52KrZur2eN8d/47ZgfPpm6husxrNg83Bzezx8iKH0+FMHzLOQq4O/JDjwaULyQTKnISKdaEECKnCQmCXSMh5DgUrAIdf4KSTa2dVY6Rdu8eIX37kXr1KvkHDSLfZ5+avb+ZURmZEzSHVRdX0bpYa6Y2moqDreUvSRuMihm7LrP00E3qlvBicbda0kMtB9K1WNM0rS0wH7AFViilpj31/GDgEyADeAD8Ryl15/FzBuDc45feVUp10DNXIYTI9pKjYedwOLsBXAtAh6+g+vtgY75GrbldwuHfCR06FJTCd9lS3Bo1MnuMdGM6AUcC2HFzB53LdWZEnRHYWuHfKC4lnQHrTrH/ygO61fNjbPtK0pojh9KtWNM0zRZYCLQCQoATmqZtU0pdzPSyU4C/UipJ07TewAzgvcfPJSulquuVnxBC5Cj3jsPm/0B8BDQaCq8MBEeZ25hVSikeLV/Bg7lzcSxTBp+vv8LBz/w3+SelJzH44GCOhB6hX41+fFrF/Kt2WXHrYSKfrDrBnUdJTHyjMt3rFbN4DsJ89FxZqwNcV0rdBNA0bT3QEfhvsaaU2p/p9ccA67RwFkKI7MpohKPzYd9EyOMD//kVfGpZO6scxZCQSPioUcT/+ise7dpReNJEbFzMP70hOiWaz/d9zoVHFxhXfxydynYye4ysOHztAZ+vPYmtjUZgj7rUL5XPKnkI89GzWCsK3Mv0OAR43qTgHsDOTI+dNE0LwnSJdJpSauvTB2ia9hnwGYCfDn8hCSGEVSU8gB97wo19UPEN6LAAnMw3RPxlkHb7Nvf69iXt5i0KDBuG18cf6bLSFZoQSq89vQhPDGde03k08zN/+4+/o5Ri5ZHbTP75ImULurP8A398vWSkWG6gZ7H2rJ8G9cwXalo3wB9okunLfkqpME3TSgK/aZp2Til1439OptQyYBmAv7//M88thBA50s2D8MOnprmer8+FWh9LY9t/KP7AAcK+GIZma4vfiuW4NmigS5wrUVfovbc3KYYUlrdeTo0CNXSJ8zypGQZG/3ieTcEhtKlUkDnvVsfVUfYQ5hZ6/kuGAL6ZHvsAYU+/SNO0lsAooIlSKvXJ15VSYY8/3tQ07QBQA7jx9PFCCJGrGDLg4HTT4HXvMtDtByhU2dpZ5SjKaOTh4sU8/HohjhXK47PgKxx8iuoS60TECQb8NgBne2dWt11Nac/SusR5nsj4FHoFBnPybgz9W5RhYIsyMjoql9GzWDsBlNE0rQQQCnQGumZ+gaZpNYClQFulVGSmr3sCSUqpVE3TvIGGmDYfCCFE7hUXBls+gTtHoFpXaDcTHK3TlyunMiQkEDZsOAm//YZHh/YUnjABGycnXWLtvbOX4YeG4+Puw5KWSyjsVliXOM9zPjSWT1cHEZ2UxsKuNXmtquVzEPrTrVhTSmVomtYX+BVT646VSqkLmqZNAIKUUtuAmYAbsOnxPQRPWnRUAJZqmmYEbDDds3bxmYGEECI3uLobtvYyTSR4YwlU72LtjHKctLt3ude7D2m3b1Pwyy/x7N5Nt52YG69sZPKfk6nsXZmFzReS1ymvLnGeZ/uZML7YfAYvFwc292pA5aJyP2NupSmVO2718vf3V0FBQdZOQwgh/pmMNPhtAhz9CgpWhre/hfxlrZ1VjpN4/Dih/QeglMJn/nxc6z1vP9u/p5RiyZklLDqziMY+jZnVZBbOds66xPorRqNi9p4rLNx/A/9inizuVov87o4WzUG8OE3TgpVS/ll5rdx9KIQQ1hJ929Q7LTTYNHS9zWSwt+wv/twgZvNmwsdPwMHHB9/Fi3AoXlyXOAajgSl/TmHj1Y10KNWBcQ3GYW9jr0usv5KQmsHA9afZe+k+7/n7MvGNyjjYSaPb3E6KNSGEsIaLP8FP/QAF76yCSm9YO6McRxkMRM6aTdS33+LaoAFF583F1kOfmZephlRGHBrB3rt76VG5BwNqDrB4s9u7j5L4ZPUJbjxIZFz7inzYoLhVGu4Ky5NiTQghLCk9BX79EoK+gaK14O2V4Fnc2lnlOIaEBMKGDCXh4EE833+fgiNHoNnp8ystLi2OAb8NIOh+EMNqD6N7xe66xHmeozce0mftSZSCVR/X4ZUy3hbPQViPFGtCCGEpD6/Bpo/h/jmo3xdajAU7Gar9T6WFhBLSuzepN29SMGAMXl27/v1B/1JkUiS99/bmZuxNpjeaTruS7XSL9SxKKQKP3WH89ouU9HZl+Qf+FPd2tWgOwvqkWBNCCEs4sx52DAY7R+i6Ecq2sXZGOVJScDAh/fqjMjLwW75Mt0a3ALdjb9NzT09iUmNY2GIhDYroF+tZ0jKMjN12gXXH79KifAHmda6Ou5Nl75ET2YMUa0IIoae0RPjlCzi9FvwaQKcVkEefBq25XcyPW4kICMC+SBF8Fi/GsWQJ3WKde3COz/d9jqZprGy7kkr5KukW61keJaTSe81Jjt+Ook/TUgxpXQ5baXT70pJiTQgh9BJxHjZ/bLr82XgYNBkOtvK2+08po5EHc+fyaPkKXOrVw2feXGzz6tfX7EjoEQYdGISXkxdLWy2lmEcx3WI9y8WwOD5dHcTDhFTmd65Ox+pS3L/s5F1DCCHMTSkI/hZ2jTQNXv9gK5Rsau2sciRjYiKhw4aTsG8fed97j0KjR6HZ63cpcPuN7QQcCaC0Z2kWt1yMt7Nlb+TfeS6cwRvPkMfZnk296lPVx/LNdkX2I8WaEEKYU0osbB8AF36EUs3hzaXgVsDaWeVI6WFh3OvzOalXr1Jw1Cg8u72va6uKVRdWMStoFnUK1WF+s/m4OVhu1JfRqJi/7xrz912jhl9elnarRQEPfcZkiZxHijUhhDCX0GBTk9uYe6adng0Hgo00LP03kk6dMm0kSEnBd+kS3Bo10i2WURmZGzyX7y58R6tirZjWaBoOtpbbpZuYmsGQjWfYdSGCTjV9mPxmZZzsbS0WX2R/UqwJIcSLUgqOLYI9Y8G9EHy8E/z0GXf0Mojdvp3wUaOxK1gQ31Xf4ViqlG6x0o3pjD0ylu03t9O5XGdG1BmBrY3lCqV7UUl8ujqIq/fjGf1aBXq8UkIa3Yr/Q4o1IYR4EUlRsLUPXN0J5V6Djl+Di5e1s8qRlNHIgwULeLRkKS61a1N0wXzsPD11i5eUnsSQg0P4PfR3+lbvy2dVP7NooXT8VhS91gSTbjCy8qPaNC0nl8vFs0mxJoQQ/9adP2BLD0h8AG2nQ92eIKsi/4oxKYmw4SOI37OHPG93onBAAJqDfpcio1Oi+Xzf51x4dIFx9cfRqWwn3WI9y/d/3iXgp/P4ebmw/EN/SuW33P1xIufJcrGmaZqrUipRz2SEECJHMBrg9zmwfyrk9YMeu6FIDWtnlWOlR0QQ0udzUi5dosCI4Xh9+KGuK1xhCWH03NOT8MRw5jadS3O/5rrFelq6wcjEHRdZ/ccdmpTNz4IuNcjjLI1uxfP9bbGmaVoDYAXgBvhpmlYN6KmU6qN3ckIIke3E34cfP4ObB6ByJ3h9HjjpMzz8ZZB89iwhn/fFmJiIz+JFuDdtqmu8q9FX6b2nN8mGZJa1WkbNgjV1jZdZdGIafdae5I+bj/iscUmGty0vjW5FlmRlZW0u0AbYBqCUOqNpWmNdsxJCiOzoxm/ww2eQmgDtF0DND+Sy5wuI++UXwkZ+iZ23N8W+WYFT2bK6xguKCKL/b/1xtndmVdtVlPEso2u8zK5ExPPJ6hPcj01l9jvV6FTLx2KxRc6XpcugSql7Ty1JG/RJRwghsiFDBhyYAofnQP5y8ME2KFjR2lnlWEopHn69kIcLF+JcqxY+Xy3AzkvfTRn77uxj2KFhFHUvytKWSynsVljXeJntvhDBoA2ncXG0Y33PetT002/ThMidslKs3Xt8KVRpmuYA9Acu6ZuWEEJkE7EhsLkH3DsGNbrDqzPAwcXaWeVYxpQUwkaOJH7nLvK88QaFJozHRseNBACbrm5i0rFJVPauzMLmC8nrZJmpAEopFu6/zqzdV6nqk4dl3f0plEca3Yp/LivFWi9gPlAUCAF2A3K/mhAi97uyE7b2BkM6vLUCqr5j7YxytPT7kYT07UvK+fMUGDoErx49dN1IoJRiydklLDq9iMY+jZnZeCYu9pYptJPTDHyx+Qw7zobTsXoRpneqKo1uxb+WlWKtnFLq/cxf0DStIXBEn5SEEMLKMtJg71hTo9tCVeGd7yCffo1ZXwbJFy4Q0udzDPHx+Cz8Gvfm+u7ANBgNTD0+lQ1XNtChVAfGNRiHvY1ldl2GxSTzWWAQF8LiGN62PL2alJRGt+KFZKVY+wp4ervMs74mhBA5X9RN08iosFNQpye0ngh2jtbOKkeL2/UrYSNGYOvlSfF13+NUrpyu8VINqYw8PJI9d/bwn8r/YWDNgRYrloLvRNEz8CQp6QZWfOBPiwoFLRJX5G5/WaxpmlYfaADk1zRtcKanPABZyxVC5D7nt8C2AaZ5nu+tgQrtrZ1RjqaU4tGSJTyYvwDn6tXx+for7Ly9dY0ZnxZP/9/6E3Q/iGG1h9G9Yndd42W28cQ9Rm89T+G8Tqz7tC5lCrpbLLbI3Z63suaAqbeaHZD5Oy4OeFvPpIQQwqLSk2HXCAj+Dnxqw9srTc1uxb9mTE0lfNRo4nbswKN9ewpPmoiNo74rlDdjbzJw/0Duxd9jWqNpvFbyNV3jPZFhMDLll8usPHKLV0p783XXGuR1sdwgeJH7/WWxppQ6CBzUNO07pdQdC+YkhBCW8+AKbPoIIi9Cw4HQfDTYSkf5F5Hx4AH3+vYl5cxZ8g8cSL6e+s/c3HNnD6N/H42TnRPLWi2jdqHausZ7IjYpnb7rTnL42kM+blicUe0qYGdrY5HY4uWRlXvWkjRNmwlUAv6751gpZbn5HEIIYW5Kwenv4ZehYO8C72+BMi2tnVWOl3LpEvf6fI4hJoaiC+bj0bq1rvEyjBl8deorVp5fSVXvqsxuOptCroV0jfnE9ch4PlkVRGhMMtM7VeG92rIaK/SRlWJtLbABeB1TG48PgQd6JiWEELpKjYefh8DZDVC8Eby1HDws1yQ1t4rfu5fQL4ZhmycPxdeuwamivo2Do1KiGHZoGH+G/8m7Zd9leJ3hONha5vLj/suR9F93Ckd7G9Z9Wg//4vo29RUvt6wUa/mUUt9omjYg06XRg3onJoQQugg/C5s/Nu36bPolNB4KNrJn6kUopXi0YgUP5szFqUoVfL7+CvsCBXSNeeHhBQYeGEhUchQTGkzgzTJv6hrvCaUUSw/dZPquy1Qo5MHyD/0pmtfZIrHFyysrxVr644/hmqa9BoQBMtRMCJGzKAUnVsCvo8DFCz7cDsVfsXZWOZ4xLY2IMQHE/vQTHu1epfCUKdg46dul/4drPzD52GS8nb1Z3W41lfJV0jXeEynpBkZsOcvW02G8VqUwM9+piotDlqY2CvFCsvJdNknTtDzAEEz91TyAQbpmJYQQ5pQcA9v6wqXtULoVvLkEXPVtIfEyyHj0iJB+/Uk+eRLvfn3x7tNH140EaYY0ph6fyuarm6lXuB4zGs/A08kyczYjYlPoGRjEmZBYhrYuy+fNSkujW2Exzy3WNE2zBcoopXYAsUAzi2QlhBDmcu+EqcltfBi0mgj1+5r6qIkXknLlKiG9e5MRFUXReXPxaNtW13gRiREMPjCYcw/P8UmVT+hbvS+2Frp8ffpeDJ+tDiIhNYOl3WvRppJlNjAI8cRzizWllEHTtA7AXAvlI4QQ5mE0wh9fw77x4FEE/vMr+PhbO6tcIX7/fsKGDMXGzY1igYE4V6msa7zj4cf54tAXpBpSmdd0Hi2KtdA1XmY/nAxhxA/nKODuyA89GlC+kIfFYgvxRFYugx7VNO1rTDtCE598USl1UreshBDiRSQ+hB97wfU9pikEHb4G57zWzirHU0oR9e13RM6ciVOFCvgsXoR9Qf3GKSmlWHVhFXNPzqW4R3HmNptLyTwldYuXmcGomLHrMksP3aRuCS8Wd6uFl6s0uhXWkZVircHjjxMyfU0B0mdNCJH93P4dtnwCSY+g3Syo/QnIvUUvTKWlET5+PLFbfsC9TRuKTJuKjbN+uyAT0xMJOBLA7ju7aVWsFRMbTsTV3lW3eJnFpaTTf90pDlx5QPd6xQhoXxF7aXQrrOhvizWllNynJoTI/owGODQTDk4Hr5LQdSMUrmrtrHKFjOhoQvv1JykoCO8+vfHu2xdNx/v+bsXeYtD+QdyKu8XgWoP5qNJHFruZ/+aDBD5ZHcTdR0lMeqMy3eoVs0hcIZ5H9hwLIXK+uHD44VO4fRiqvgevzQZHGaJtDqnXr3Ovdx8y7t+nyKxZ5Hld33mb++7uY9Tvo3CwcWBZq2XULVxX13iZHbr6gL7fn8TWRmPNJ3WpVzKfxWIL8Ty6rutqmtZW07QrmqZd1zRtxDOeH6xp2kVN085qmrZP07RimZ77UNO0a4//+1DPPIUQOdi1vbDkFQgNho6L4M2lUqiZScKhQ9zu3AVjcjLFAlfrWqgZjAYWnFzAwP0DKeFRgg2vb7BYoaaUYsXhm3z07XGK5HVmW99XpFAT2YpuK2uP234sBFoBIcAJTdO2KaUuZnrZKcBfKZWkaVpvYAbwnqZpXsBYwB/T/XHBj4+N1itfIUQOY0iH3ybBkXlQoCK88x3kL2ftrHIFpRTRgYHcnzYdx7Jl8V28CPvC+o3jikmJYfjh4RwNO0qnMp0YWXckjraOusXLLDXDwKgfz7M5OIQ2lQoy593quDrKRSeRvfztd6SmaS6YGuL6KaU+1TStDFDuce+156kDXFdK3Xx8nvVAR+C/xZpSan+m1x8Duj3+vA2wRykV9fjYPUBbYF2W/lcJIXK3mLuwuQeEHIdaH0PbqWAvI3/MQaWnEzFxEjEbN+LWsgVFp0/HxlW/G/svPrrIoP2DeJD8gHH1x9GpbCfdYj0tMj6FXoHBnLwbw4AWZRjQogw2NrIZRWQ/Wfnz4VsgGKj/+HEIsAn4u2KtKHAv0+MQ4Hlr2j2Anc85tujTB2ia9hnwGYCfn9/fpCOEyBUu7YCf+pj6qL29Eipb7pd7bmeIiSFkwECS/vyTfJ99Rv6BA3TdSLD1+lYmHZuEp5Mnq9quokr+KrrFetq5kFg+CwwiJimdRe/XpF0V/VYOhXhRWSnWSiml3tM0rQuAUipZy9q2nGe9Rj3zhZrWDdMlzyb/5Fil1DJgGYC/v/8zzy2EyCXSU2DPGDi+DApXh3e+Ne36FGaRevMW93r3IiMsnCLTp5GnY0fdYqUb0pl+YjobrmygbqG6zGgyAy8nL93iPW3bmTC+2HQGbzdHNveuT6UieSwWW4h/IyvFWpqmac48LpY0TSsFpGbhuBDAN9NjH0xD4P+HpmktgVFAE6VUaqZjmz517IEsxBRC5EaPbsCmjyDiLNTrAy3HgZ1l7ml6GSQcOULowEFo9vb4rVqFS80ausW6n3ifwQcHc/bBWT6u/DH9a/THzsYy94gZjYrZe66wcP8N/It5sqR7Lbzd5PtIZH9Z+QkZC+wCfDVNWws0BD7KwnEngDKappUAQoHOQNfML9A0rQawFGirlIrM9NSvwBRN055M6G0NjMxCTCFEbnN2E+wYCLb20GU9lHvV2hnlKlFr13J/ylQcS5UybSQo+n/uODGbExEnGHpwKCkZKcxuMpvWxVvrFutp8SnpDNpwmr2XIulc25cJHSvjYCeNbkXO8HeD3DXgMvAWUA/T5ckBSqmHf3dipVSGpml9MRVetsBKpdQFTdMmAEFKqW3ATMAN2PT4yupdpVQHpVSUpmkTMRV8ABOebDYQQrwk0hJh5zA4tQZ868Hb30AeH2tnlWuojAzuT5lC9PfrcGvWjCIzZ2Lrps9GAqUUgRcDmRM8B193X75t8y0l81ruEvadR4l8ujqIGw8SGd+hEh/UL2axJrtCmIOm1PNv9dI0LVgpVctC+fxr/v7+KigoyNppCCHMIfKS6bLngyvQaAg0HQm20k7BXAyxsYQOGkTi0T/w6vEfCgwejGZrq0uspPQkxh4dy67bu2jh14JJDSfh5uCmS6xnOXL9IZ9/fxKlYNH7NWlY2ttisYV4nsf1lX9WXpuVd79jmqbVVkqd+PuXCiHEC1AKTq6GncPB0Q26/wClZAyxOaXdvs293n1ICwmh8OTJ5O30lm6x7sTdYeD+gdyMvcmAmgPoUbmHxVa0lFKs/uMOE3ZcpKS3Kys+9KdYPsvMFhXC3LJSrDUDemqadgdIxHQpVCmlZOieEMJ8Ii/B7jFwfQ+UaAJvLQf3gtbOKldJPHaMkAED0TSNYt+uxMU/S3/U/ysH7h1g5OGR2NnYsbjlYhoUaaBbrKelZRgZu+08647fo0X5AszrXB13J3uLxRfC3LJSrMndvEII/cTfh/2T4VQgOLhD68lQrzfY6HNZ7mUVvX4DEZMm4VC8GL6LF+Pg6/v3B/0LBqOBxWcWs/TsUip4VWBes3kUcSuiS6xneZiQSu81wZy4HU2fpqUY0rocttLoVuRwWSnWpH+ZEML80hLh6NdwZD4YUqFOT2gyDFws12/rZaAyMrg/YwbRqwNxbdyIonPmYOumzz1jsamxDD88nCOhR3ij9BuMqjsKJzsnXWI9y56L9xn5w1niUzKY37k6Havrt7NVCEvKSrH2M6aCTQOcgBLAFaCSjnkJIXIrowFOr4XfJkNCBFTsCC3GQr5S1s4s1zHExxM6eAiJhw/j9eGHFBj2hW4bCS5HXWbg/oHcT7rPmHpjeKfsOxa7Py0uJZ0J2y+yOTiECoU9COxRjQqFPSwSWwhL+NtiTSn1P/M/NE2rCfTULSMhRO51fS/sDoDIC+BTG95dDX7Pm0In/q20u3dNGwnu3KHQhPF4vvuubrG239jO+D/Gk8cxD9+1/Y5q+avpFutpR64/5ItNZ4iIS6Fvs9L0b1FG+qeJXOcf74VXSp3UNK22HskIIXKpiPOmUVE3fgPP4vDOd1DxDZBeV7pIPH6c0P4DQCn8vvkG17p1dImTbkhnZtBM1l1eh39Bf2Y2mYm3s2VaYySnGZi28xKr/rhDSW9XtvRuQA0/z78/UIgc6G+LNU3TBmd6aAPUBB7olpEQIveICzNd7jy9FpzyQJupULuHjIrSUczmzYSPn4CDry++ixfhUKyYLnEikyIZcmAIpx+c5oOKHzCo1iCLjY06eTeaIRvPcOthIh81KM7wtuVxdpANKSL3yspPlnumzzMw3cO2RZ90hBC5Qmo8HFkAR78CZYD6n0PjoeAsKx96UQYDkbNmE/Xtt7g2bEjRuXOw9dDnvq3g+8EMPTiUxPREZjaeSdsSbXWJ87TUDAPz915jycEbFM7jzPef1KWBNLkVL4Gs3LM23hKJCCFyAUMGnFoN+6dCYiRU7gQtAkyXPoVuDAkJhA0ZSsLBg3h260bBEcPR7My/yqWU4vvL3zPrxCyKuhdleavllPYsbfY4z3IpPI5BG05zOSKed2r5MKZ9RTykd5p4SfzlT7Omadt5TtsOpVQHXTISQuQ8SsG13aamtg+vgF996LIOfPRruipM0kJCCendm9SbNyk0NgDPLl10iZOckcy4o+P45dYvNPVtypRXpuDu4P73B76gDIORpYduMm/vVfI4O7DiA39aVpRmyeLl8rw/vWZZLAshRM4Vdhp2j4bbh8GrFLy3Fsq/JpsHLCApOJiQfv1RGRn4LV+GawN9pgTci7vHwAMDuRZ9jX41+vFJlU+w0fTfcXnzQQJDNp3h1N0Y2lUpxKQ3quDl6qB7XCGym78s1pRSBy2ZiBAih4kNgX0T4ex6cPaCV2eC/8dgK5emLCHmx61EBARgX6QIPksW41iihC5xDoUcYsThEWhoLGq5iFeKvqJLnMyMRkXgsTtM3XkJRztb5neuTodqRSzWt02I7OZ5l0E3KqXe1TTtHM+4HCqzQYV4SaXEwu9z4dhi0+XPhgOh0WDTbk+hO2U08mDuXB4tX4FL/Xr4zJuHbR7z/39vVEaWnlnK4jOLKedVjrlN5+Lj7mP2OE8LjUlm2OYzHLn+iCZl8zPj7aoU9LDcFAQhsqPnXQYd8Pjj65ZIRAiRzRnSIfg7ODAVkh5B1feg+WjI62ftzF4aGVFRhI8eQ8Jvv5G383sUGjUKzd78K5mxqbF8+fuXHAo5RIdSHRhdbzTOds5mj5OZUootJ0MZv+0CBqWY8mYVutTxldU0IXj+ZdDwxx/vWC4dIUS2oxRc+QX2BMCj61C8EbSeCEVqWDuzl4ZSirjt27k/ZSqGxEQKjhqFZ7f3dSlkrkRdYdCBQYQnhDOq7ijeAFwsiwAAIABJREFUK/ee7gXTg/hUvvzxHHsu3qdOcS9mvVMNv3wuusYUIifJSlPcesBXQAXAAbAFEpVSMnhNiNwuNNi0w/POEfAuC13WQ9m2snnAgtJDQwkfN57Ew4dxrlaNwpMm4limjC6xfr75M+OOjsPdwZ1v235L9QLVdYmT2c5z4Yzaep6E1AxGtavAf14pga2NfH8JkVlWGvF8DXQGNgH+wAeAZRrrCCGsI/oO7JsA5zeDa354bQ7U/BBsLdOhXpia3Eav/Z7IefMATKtpXbvoMog93ZjOnKA5rLm0hpoFajK76Wzdx0bFJqUzdtt5tp4Oo3JRD+a+W50yBfVvBSJETpSld16l1HVN02yVUgbgW03TjuqclxDCGpKj4fBs+HMpaLbQ+AtoOAAc5ZeoJaVeu0b46DEknzmDa+NGFB43DvsiRXSJ9TD5IUMODOFk5Em6VejGYP/B2Nvou6P34NUHDN98lgcJqQxoUYa+zUtjbyvD14X4K1kp1pI0TXMATmuaNgMIB1z1TUsIYVEZaRD0DRycDskxUL0rNBsFeYpaO7OXijEtjUdLl/Fw2TJsXV0pMnMGHq+/rts9Y6cjTzP4wGDi0+KZ3mg67Uq20yXOE4mpGUz55RJr/7xLmQJuLP/Anyo+sotYiL+TlWKtO6YB7n2BQYAv0EnPpIQQFqIUXPwJ9o6D6FtQsim0mgiFpTOPpSWdOkX4mDGkXb+BR/v2FBw5AjsvL11iKaVYf2U9M07MoLBrYRa3NLXn0NOJ21EM2XiGe9FJfNqoBENal8PJXoavC5EVz+uz5qeUuptpN2gKIHNChcgt7h2HX0dByHHIXwHe3wKlW8jmAQszJCTyYN48oteuxa5QIXyXLcWtcWPd4qVkpDDx2ES23dhGY5/GTHllCnkc9VvdSkk3MHfPVZYdvomPpzPrP61H3ZL5dIsnRG70vJW1rUBNAE3TtiilZDVNiNwg6ibsHQ8Xt4JbQWi/AKq/L5sHrCDh0CHCx40jIzwCz/ffJ//Agdi66XeXSUh8CIMODOJK1BX6VO9Dz6o9dR0bdT40lsEbT3P1fgJd6vgx6rUKuDnK95kQ/9Tzfmoy/3ldUu9EhBA6S4qCQzPh+HLTSKimI6F+X3B0s3ZmL52MqCjuT51G3PbtOJQqRbHv1+JSQ9++db+H/s7wQ8NRKL5u8TWNffRbvcswGFl04AYL9l3Dy9WBbz+uTbNyBXSLJ0Ru97xiTf3F50KInCQjFY4vMxVqqfFQo5tp84B7IWtn9tJRShG3Y4epuW1CAt6ff06+np9h46DfcHKjMrL87HIWnl5IGc8yzGs6D18PX93iXY+MZ8jGM5wJiaVDtSJM6FiJvC4yfF2IF/G8Yq2apmlxmFbYnB9/zuPHSpriCpHNKQXnt8C+8RBzF0q3glYToGBFa2f2UkoPDSV8/HgSD+nf3PaJ+LR4vjz8JQdCDvBaydcYW3+sbmOjjEbFt0dvM2PXZVwcbFnYtSavVS2sSywhXjbPGzcl23SEyKnuHIXdo00TCApWhu4/Qqnm1s7qpfR/mtt++SWe73fVpbltZteirzHowCBC40MZUWcEXct31a0FyL2oJIZuOsOft6JoWaEAU96qQgF3Gb4uhLnInZ5C5CYPr8PesXB5B7gXgY6LoFpnsJG/vazhf5rbNmpE4XFjsS+qf++6Xbd2EXA0AFd7V75p8w01C9bUJY5Sig0n7jFxx0U0TWPG21V5p5aPDF8XwsykWBMiN0h8aGpoG7QS7Jyg+Wio9zk4yDBsa7B0c9sn0o3pzA2eS+DFQGoUqMHsJrPJ75Jfl1iRcSmM+OEcv12OpH7JfMx8pyo+nvL9JoQepFgTIidLT4Y/l8DhOZCWCLU+NO3ydJOdd9Ziyea2mT1MfsgXB78g6H4QXcp34Qv/L7C31Wds1PYzYYz56TzJaQbGtq/Ih/WLYyPD14XQjRRrQuRERiOc22Qath4XAmVfhVbjIb++XejFXzMmJhI5N1Nz26VLcGvSxCKxzzw4w+ADg4lNjWXKK1NoX6q9LnGiE9MY89N5dpwNp5pvXua8W41S+aX1ixB6k2JNiJzm1iHT5oHwM1C4Gry5GEro1zNL/L2Ew4cJHzvWYs1tn1BKsenqJqYen0pBl4KsabeG8l7ldYm1/3Ikw7acJSYpjaGty9KrSSnsZPi6EBYhxZoQOcWDK7BnLFzdCR4+8OYyqPIO2MgvTGvJiI7m/pSp/7+57dq1uNTUt7ntEykZKUz+czJbr2+lYdGGTG80XZexUQmpGUzacZH1J+5RvpA7331cm0pFZPi6EJYkxZoQ2V1CJByYCsGrwMEVWo6Dur3AXp9+WeLvmZrb/sz9KVMs1tw2s9CEUAbtH8SlqEv0rNqT3tV6Y6vDjt8/bjzii81nCItJpnfTUgxsWQZHO9lZLISlSbEmRHaVlgTHFsLv8yAjBWp/Ak2Ggau3tTN7qaWHhZma2x48hFO1qvhNnIhT2bIWi3807CjDDg3DaDTyVfOvaOrb1OwxUtINzNh1hZVHblE8nwubetWnVjH9N0kIIZ5NijUhshujAc6sh98mQXwYlH8dWo4H79LWzuylpoxGor9fx4M5c1BYrrntE0npSSw4tYDvL31PqbylmNdsHsU8ipk9zpl7MQzeeJobDxL5oH4xRrxaHhcH+VUhhDXp+hOoaVpbYD5gC6xQSk176vnGwDygKtBZKbU503MG4Nzjh3eVUh30zFWIbOHGftg9Bu6fgyI14e1voFgDa2f10ku9ft3U3Pb0aVxfeYVC48bh4KN/c9snTkScIOBIACEJIXQu15lBtQbhYm/enmbpBiNf7bvGwgM3KODuyJoedXmljKziCpEd6FasaZpmCywEWgEhwAlN07YppS5metld4CNg6DNOkayUqq5XfkJkK/cvwp4xcH0v5PWDTt9Apbdk84CVqbQ0Hi5fzqMlS7FxdaXIjOl4tG9vsQ79SelJzAmew4YrG/B192Vlm5XULlTb7HGuRMQzeONpLoTF8VbNooxtX4k8zvr0aBNC/HN6rqzVAa4rpW4CaJq2HugI/LdYU0rdfvycUcc8hMi+4iNg/2Q4tQYc3aH1JKjzGdg5Wjuzl17y6dOEjxlD6rXreLz+OgW/HGmR5rZPHAs/xtgjYwlPDKdbhW70q9HP7KtpBqNixeGbzN59FXcnO5Z2r0WbSoXMGkMI8eL0LNaKAvcyPQ4B6v6D4500TQsCMoBpSqmtT79A07TPgM8A/Pz8XiBVISwsNQH++BqOzAdDuml3Z+MvwEVu4rY2Y2IikfPmE71mDXaFCuGzZDHuTZtaLH5CWgKzg2ez+epminkUY9Wrq6hRwPztQO48SmTopjOcuB1Nm0oFmfxmFbzd5I8EIbIjPYu1Z10nUP/geD+lVJimaSWB3zRNO6eUuvE/J1NqGbAMwN/f/5+cWwjrMBpMq2j7J0PCfaj4BrQcC14lrZ2Z4Knmtl27kn/QIIs0t33iSOgRxv0xjsikSD6q9BGfV/8cJzsns8ZQSrH2z7tM+eUStjYac96txps1isrwdSGyMT2LtRDAN9NjHyAsqwcrpcIef7ypadoBoAZw47kHCZGdhQTBtv4QeQF868J7a8C3jrWzEjxubjt1KnHbLN/cFiAuLY5ZJ2bx4/UfKZmnJIGvBlI1f1Wzx4mITWHYlrMcuvqARmW8mfF2VQrnkX59QmR3ehZrJ4AymqaVAEKBzkDXrByoaZonkKSUStU0zRtoCMzQLVMh9JSWaGrDcWwxeBSBd1ZBxY4gKxlW93+a2/bpQ75ePS3W3BbgUMghxv8xnofJD+lRuQe9q/fG0da8lyOVUvx0OoyAn86TblBMfKMy3er6yWqaEDmEbsWaUipD07S+wK+YWnesVEpd0DRtAhCklNqmaVpt4EfAE2ivadp4pVQloAKw9PHGAxtM96xd/ItQQmRfN/bD9v4Qcxf8e5imDzh5WDsrgfWb28amxjL9+HS239xO6bylWdBsAZW8K5k9zqOEVEZvPc/O8xHUKubJ7HeqUdzbcpd2hRAvTlMqd9zq5e/vr4KCgqydhhAmSVGmYeun10K+0tDhK+mXlk38T3NbpSgwaCCe779vsea2AL/d/Y2JxyYSnRJNjyo96Fm1Jw625l/N23PxPiN/OEtccgaDW5fl00YlsbWR1TQhsgNN04KVUv5Zea20pRbCnJSCiz/BL19A0iN4ZTA0GQ725r1JXPw7qTdumJrbnjqFa8OGFBo/3qLNbaNTopl6fCo7b+2knGc5FrVYRIV8FcweJy4lnQnbL7I5OISKhT1Y80k1yheSFV0hciop1oQwl7hw+GUoXN4BhatBty1Q2Pw3iYt/7n+a27q4UGT6NDw6dLDoPVt77uxh0rFJxKXG0adaHz6p8gn2tuZvPHvk+kO+2HSG+/Gp9Gtemn7Ny+BgJ82VhcjJpFgT4kUpBSdXm8ZEGVKh1QSo9znYyo9XdpB85gzho8eQeu0aHq+9Zmpumy+fxeI/Sn7ElD+nsPvObip4VWBZq2WU8ypn9jjJaQam7bzEqj/uUDK/K1t6N6C6b16zxxFCWJ78NhHiRTy6AdsHwO3DUOwV6LAA8pWydlaCx81t588nOtA6zW2VUvx6+1em/DmFhPQE+tXox8eVP8bexvyraSfvRjNk4xluPUzkPw1LMKxtOZzsLXcPnhBCX1KsCfFvGDLg2CLYPwVs7eH1eVDzQ5nlmU0kHP6diLFjSQ8Px7NLF/IPHmzR5rYPkx8y+dhk9t7dS+V8lZnYcCKlPUubPU5qhoH5e6+x5OANCudxZt2n9ahfynKrhkIIy5BiTYh/KuI8bOsLYaegXDt4bbapf5qwuozoaCKnTSP2p204lCxJsbVrcKlZ02LxlVL8fOtnph2fRnJ6MoNqDeKDih9gZ2P+t9pL4XEM2nCayxHxvOfvy+jXK+DuJMPXhciNpFgTIqvSU+DQTDgyD5w94e1vodKb0tw2G1BKEffzL6bmtnFxePfpTb5evSza3DYyKZKJf0zkQMgBquavysSGEymZx/xjxDIMRpYeusm8vVfJ6+LANx/606JCQbPHEUJkH1KsCZEVd4/Btn7w8CpU6wJtpsjQ9WwiPTyciHHjSTh4EKeqj5vblrNcc1ulFNtubGP6iemkGdIY6j+UbhW6YWtj/nvGLoXH8eWP5zh1N4bXqhZmUsfKeLpariAVQliHFGtCPE9qPOybAMeXQx5fUzuO0i2tnZXgcXPbdet4MNvU3LbgyBF4dutm0ea2EYkRjP9jPL+H/k7NAjUZ32A8xfMUN3uchNQM5u25yrdHb5PH2Z4FXWrQoZpcehfiZSHFmhB/5doe2D4Q4kKhbk9oPgYc3aydlcD6zW2VUvx4/UdmnphJhjGDEXVG0KV8F2w0824wUUqx83wEE7Zf5H58Cl3q+DGsTTnyushqmhAvEynWhHha4iPYNQLObQTvctBjN/jWsXZWgsfNbVes4NHiJdi4uFB42lTydOxo0ea24QnhjPtjHEfDjuJf0J8JDSbg6+Fr9ji3HyYSsO0Ch64+oFIRDxZ3q0kNP0+zxxFCZH9SrAnxhFJwfgvsHAYpcaYxUY2GgJ2jtTMTPNXctl07Co760qLNbZVSbLq6idlBs1EoRtUdxbvl3jX7alpKuoElB2+w6MANHG1tGNe+It3qFcPOVtrCCPGykmJNCIDYUPh5MFzdBUVrQYevoWBFa2clAGNSEg/mzydqdSB2BQvis3gR7s2aWTSHkPgQxh0dx58Rf1K3cF3G1R+Hj7uP2eMcuvqAgJ/Oc/tREh2qFWH0axUo4CFzZYV42UmxJl5uRiMEr4Q940AZTLs86/YCHXbyiX8u4fcjpua2oaF4du1K/sGDsHWz3H2DRmVkw5UNzA2ei41mQ0D9AN4u87bZL7tGxKYwccdFfj4XTklvV9Z+UpeGpb3NGkMIkXNJsSZeXg+vm9px3D0KJZuaphB4lbB2VoInzW2nE/vTT6bmtt+vtWhzW4B7cfcIOBpA0P0gGhZpyNj6YynsVtisMTIMRr47epu5e66SYVQMbV2WTxuXxNFO/lgQQvx/UqyJl48hHY5+BQemgb0TdFwI1d+X5rbZgFKKuF9+4f7kTM1te/bExtFy9w0ajAa+v/w9C04uwN7GngkNJvBG6TfMvpoWfCeKUT+e53JEPM3K5Wd8h8r45XMxawwhRO4gxZp4uYSdNo2KijgHFTpAu1ngLt3frU2lpxO3ezfRqwNJPnPGKs1tAW7F3iLgSACnH5ymsU9jAuoFUNDVvN8f0YlpTNt5mQ1B9yiSx4ml3WvRumJBi+5oFULkLFKsiZdDejIcmApHvwZXb3g3ECp2sHZWL72M6GhiNmwket06Mu7fx76YH4XGBpD33Xct2tzWYDQQeDGQr09/jYOtA5NfmUz7ku3NWkAZjYpNwfeYtvMy8SkZ9GxSkv7Ny+DqKG/DQojnk3cJkfvd/h229YeoG1CjO7SeaJrtKawm5epVogMDid22HZWaimuD+hQaNxa3Jk3QbCzbouJmzE3GHBnD2YdnaerblIB6AeR3yW/WGJfC4xi99TzBd6KpU9yLiW9Uplwhd7PGEELkXlKsidwrJRb2jIXgb8GzOHzwk2kjgbAKZTCQcPAgUasDSTp2DM3JiTwdO+LVvRuOZcpYPJ8MYwbfXfiOxacX42zvzLRG02hXop1ZV9MSUjOYu+cq3x29TV5ne2a/U423ahaVS55CiH9EijWRO13+BX4eAgkRUL8vNBsFDnLztjUYEhKI3bKFqDVrSb93D7tChcg/ZDB5334bO0/rrHBei77GmCNjuPDoAq2KteLLul/i7Wy+VhlKKX45F8GEHReIjE+lax0/vpAxUUKIf0mKNZG7JDwwTSC48AMUqASd15ia3AqLS7tzh6g1a4ndsgVjUhLONWpQYMhg3Fu2RLOzzltPujGdledWsuTsEtzt3ZnVZBZtircxa4xbDxMJ+Ok8h689pFIRD5Z0qyVjooQQL0SKNZE7KAVnN5hmeqYlQrPR0HAA2MlKhiUppUj64w+iVgeScPAg2Nnh8WpbvLp/gHOVylbN7UrUFcYcGcOlqEu0Ld6WkXVH4uXkZbbzp6QbWHzgBosPmsZEje9QiW71imFrI5c8hRAvRoo1kfPF3IXtA+HGPvCtCx2+gvzlrJ3VS8WYnEzstu1Erwkk9dp1bPPlw7t3b/J2fg/7AgWsmlu6IZ3l55az/OxyPBw9mNt0Li2LtTRrjANXIhm77QJ3HiXRsXoRRrWTMVFCCPORYk3kXEYDnFgBe8ebGtq+OhNqfwIW3k34MksPDyf6+3XEbNyIITYWx4oVKDx1Kh7tXrVoI9u/cvHRRcYcGcPV6Ku8VvI1RtQeQV6nvGY7f3hsMhN3XOSXcxGUzC9jooQQ+pBiTeRMkZdNo6JCjkPplvD6XMjrZ+2sXgpKKZJPnSYqcDXxu/eAUri3aIHXB91x9vfPFjsd0wxpLDmzhJXnV+Ll5MWCZgto5me+4e/pBiOrZEyUEMJCpFgTOUtGGhyZB4dmgoMrvLkMqr4ro6IsQKWlEbdrF1GrA0k5fx4bDw+8PvwQz65dcfApau30/uv8w/OMOTKG6zHX6VCqA8NqDyOPYx6znT/odhSjt5rGRDUvX4DxHSrh6yU7jYUQ+pFiTeQcIcGm1bTIC1C5E7SdDm7mbV4q/q+MR4+IXr+e6PXrMTx4iEPJkhQaG0CeDh2wcXW1dnr/lWpIZdHpRXx34Tu8nb1Z2GIhjX0am+38UYlpTNt5iY1BITImSghhUVKsiewvLRH2T4Fji8CtEHRZD+VetXZWuV7KpUtErQ4kbscOVHo6ro0a4TWlO64NG1p8ysDfOR15moCjAdyKvcVbZd5iqP9Q3B3MMyHAaFRsDLrHtF2XSXg8JmpAizK4OMjbpxDCMuTdRmRvNw+YRkXF3AH//0DLceBkvkta4n8pg4H4ffuIXh1IUlAQmrMzed95G89u3XAsWdLa6f0fKRkpfH3qa1ZfXE0h10IsbbmUBkUbmO38F8PiGL31HCfvxlCnhBeT3qhM2YIyJkoIYVlSrInsKTkado+GU2vAqxR89DMUf8XaWeVahthYYjZvIXrtWtLDwrAvUoQCw4aR9+1O2Hp4WDu9Zzp5/yQBRwO4E3eHd8u+y6Bag3BzcDPLueNT0pm75xrfHb2Fp4uDjIkSQliVFGsi+7m4DX4ZCokP4ZVB0GQ42DtbO6tcKfXmTaICA4nd+hMqORkXf38KjBiOe/PmVpsy8HeS0pP46tRXrL20liJuRVjRegV1C9c1y7mVUvx8LpyJOy7+d0zUsDblyeNib5bzCyHEv5E9343Fyyk+wlSkXdoOharC+5ugcDVrZ5XrKKORxCNHiFodSOLhw2j29ni8/jpeH3THqUIFa6f3XCciThBwJICQhBC6lO/CwJoDcbE3z07MzGOiKhf1YGl3f6r7mq8nmxBC/FtSrAnrU8p0uXP3KEhPMd2XVr8v2MpqhjkZExOJ+eknogPXkHbrFrb5vfHu3w/P997DLl8+a6f3XEnpScwJnsOGKxvwdfdlZZuV1C5U2yznTkk3sOjADZYcuIGjnQ0TOlbi/boyJkoIkX1IsSasK+oWbB8Atw5CsYbQfgF4l7Z2VrlKWkgo0WvXErN5M8b4eJwqV6bIjOl4tG2L5pD9Z6ceCz/G2CNjCU8Mp1uFbvSr0c9sq2mZx0S9Ub0IX75WgQLuMiZKCJG96FqsaZrWFpgP2AIrlFLTnnq+MTAPqAp0VkptzvTch8Doxw8nKaVW6ZmrsDCjAY4tht8mgY2daQJBzY9kVJSZKKVIOnGC6MBA4vf9BpqGe+tWpoHqNarniBvlE9ISmB08m81XN1PMoxirXl1FjQI1zHLusBjTmKid5yMold+V7z+tS4NSMiZKCJE96VasaZpmCywEWgEhwAlN07YppS5metld4CNg6FPHegFjAX9AAcGPj43WK19hQfcvmJrbhgZD2bbw2hzIk3064OdkxtRU4n7+hajAQFIvXcI2Tx7y9eiBZ9cu2BcubO30suxI6BHG/TGOyKRIPqr0EZ9X/xwnuxdf8Uo3GPn2yC3m7b2GUSm+aFOOTxuVxMFO/kgQQmRfeq6s1QGuK6VuAmiath7oCPy3WFNK3X78nPGpY9sAe5RSUY+f3wO0BdbpmK/QW0YqHJoFv88Bp7zQ6RvTJIIcsMqT3aVHRhKzfj3R6zdgiIrCsUxpCk0YT5727bFxzjk7aePS4ph1YhY/Xv+REnlKsPrV1VTLb55NJiduRzH6x/NcuR9Pi/IFGCdjooQQOYSexVpR4F6mxyFAVvfXP+vY/7P0omnaZ8BnAH5+MsQ7W7v7p2k17eEVqNoZ2k4FFy9rZ5XjJZ87Z5oysGsXZGTg1qQJXh90x6V+/RxxqfMJpRQ7bu5gVtAsYlJj6FG5B72r98bR1vGFz/0oIZVpOy+zKTiEonmdWda9Fq0rFTJD1kIIYRl6FmvP+k2hzHmsUmoZsAzA398/q+cWlpSaAPsmwPFlkMcH3t8CZVpaO6scTaWnE793L1GrVpN8+jQ2Li54du6MV7f3cShWzNrp/WPXoq8x+c/JBN8Ppop3FRa3XEzFfBVf+LxGo2L9iXtM33WZxNQMejctRb/mpWVMlBAix9HzXSsE8M302AcI+wfHNn3q2ANmyUpYzvW9sH0QxN6DOp9BizHgKKN6/q2M6GhiNm0m+vvvyYiIwN7Xl4JfjiTPW29h62aezv2WlJieyOLTi1lzaQ1uDm6MrT+Wt8q8hY324vePnQ+NZfTW85y+F0Pdx2OiysiYKCFEDqVnsXYCKKNpWgkgFOgMdM3isb8CUzRN83z8uDUw0vwpCl0kRcGukXB2PXiXhf/8Cn7m6TD/Mkq9do2o1YHEbt+OSknBpV49CgWMwa1JEzRbW2un948ppfj1zq/MPD6TyORIOv2/9u48PsryWuD472Syb5MVEkjYQsIethAQARMFVFZtXVu3utWqtWq1m+2ttdXa1mptsVVbbdXe1ta21iSoLBoEEQJhC7KGNQkQyEJmsk2SmXnuHzNyqWUZlmSG5Hw/Hz/OvPPOOydPJjOHZzuZX+Qb475BfHj86Z98Go2ODp5dspPXPtlHQlQoz143mqvHapkopdSFrdOSNWOMU0Tux5N4WYBXjTFbROQJoNQYUyAiE4C3gXhgroj8yBgzwhhTLyI/xpPwATzx2WIDFcCMgS3/gne/BY4GmPYtmPYIBJ/7vKOexrjdNH30EfWvv07LqtVIWBjWeXOJv+lmwodk+Tu8s7bXtpeflvyUVYdWMSxhGM/mP3teFhAYYygq85SJqmlq46aJ/Xlk5hAtE6WU6hbEmO4x1SsnJ8eUlpb6O4yeyeWE7UWwagFUrYU+Y2HeAkgZ6e/ILjiupiZs/3qb+v/9Mx37Kwju3Zv4L32JuOuuJTj+3Hue/KXV2crvy37PH7f8kQhLBPePvZ/rh1yPJejcewb31DTxP+9s4eNdtYzqa+UnV41ktJaJUkoFOBFZZ4zJ8eVcnWmrzl5bo6dM1OrfQkMFxA+E2b/0bG5r0bfWmWivqKD+z3/G9s9/4W5uJmL0aJIfeIDYmTORkAu7d6i4opin1zzNweaDzB00l4dzHiYp4tw3oHV0uPht8S5e/GgPYSFaJkop1X3pN6o6c7YqKHkJ1r0GbTbodxFc/lMYciWch56SnsIYQ0tJCfWvvU7TsmVgsRB7xRUk3HIzEdnZ/g7vnFU2VvKzNT/jo6qPGBw3mD9e/kdyUnz6R+RpFe84wg/f2UJFfQtXj+3Ld2cN1TJRSqluS5M15buDG2DVC7Dlbc/8tOHzPQXX08b7O7ILitvhwFZYyNHX36CtvBxLfDyJ93yV+BtuJKR3L3+Hd87aXG28+umrvLL5FSxi4ZGcR/jSsC8REnTuPYQHG1p5onAr72/RMlFKqZ5DkzXtT4aLAAAgAElEQVR1am437Hzfk6Tt/xhCY2DiPTDxqxCnGxGfibbycmwFBTS89Q9cDQ2EDR1K6pNPEjtnNkFh3WMRxscHPuanJT+lorGCywdcziM5j5ASde4b0GqZKKVUT6bJmjqx9hbY9BdY9Vuo3w3WdJj5JIy7BcJj/R3dBaPj0CHsCxdiK1pI2/btEBRE9KX5JNx8C5G5E7rNlhLVzdX8bM3PWFqxlAGxA3hpxktM7jP5vFz7+DJR04f14odztUyUUqpn0WRN/afGw55qA6WvQOtR6DMOrnkVhs3XRQM+cjU0YF+0GHthIS3eFcrho7Pp/dhjxF55BcFJ3WfYrsPVwRvb3uDFTS9ijOGBsQ9w64hbCbWEnvO165ra+Ol72/mHt0zU72/JYcbw3uchaqWUurDot6/yOLzFM9S5+S1wdcDQ2Z75aP0maaF1H7hbW2kqLsZWtJCmFSugo4PQgQNJeuDrWGfPviDLQJ3OmkNreLLkSfbY9pCfns+3c79N3+j/KuF7xrRMlFJK/Sf99OvJjIHdH8AnC2BPMYREwvjbPHPSEjP8HV3AM04nzatWYy8qonHJEtwtLQQnJ5Pw5S8TO3cO4cOHd5thzuPVtNTwTOkzvLv3XfpG92XBpQu4JP2S83LtDRVH+VHhVjZWNjBpUAI/nq9lopRSSpO1nqjD4elBW/UC1GyDmFS47IeeRC0ywd/RBTRjDI6yMmyFRdjfew9XXR1BMTHEzLoS65w5RE6YcEGWgPKF0+3kze1v8sLGF2hztXHP6Hu4Y+QdhAef+5YZnx6w8dySnXyw/QhJ0aE8d/1orhqjZaKUUgo0WetZmmuh9FXPnLTmGug9Cq5+CUZ8AYLPfY5Rd9a2Zw/2oiJsRQvpqKhAQkOJzssjdu4coqdN6zarOU9m45GN/GT1T9hxdAcX97mY7078Lv1jz31ot/xwI88t3cm7m6uJDQ/m0cuHcOvkAUSH6UeTUkp9Rj8Re4KanbD6Bdj0JjgdkDnTMx9t4DSdj3YKHYcPY3/3PeyFhTi2bgURIidNJOmrXyVm5gwsMd1/eK7eUc9z657j37v+Te/I3jyb9yzT+00/5x6vvbXNPL90J+9sOkhkiIUHLsvkjikDsUZc2NUalFKqM2iy1l0ZA/tWeIY6d74PweEw+gaYdC8kD/F3dAHLZbfTuHgxtqKFtJSUgDGEjxxJr+98m9hZswjpdeFvWusLl9vFP8v/yfPrn6elo4WvjPwK92TfQ2TIuW2ZUXW0hd98sIt/rK8ixCLcPW0QX52WQUKU9uwqpdTJaLLW3TjbPRUGVi2A6jKITIK878GEOyCq+2wZcT6529poWvYR9qJCmpZ9hOnoIKR/P5LuvZfY2bMJGzTQ3yF2qS21W/jx6h+zpW4LE1Im8NjEx8iIO7cFJ4ftDhZ8uIs311YgCLdc1J+v5WVoiSillPKBJmvdRetRWPcnKHkZGg9C0hCY+2vIvh5C9Avx84zLRcuaNdgKi2hcvBh3UxOWpCTibrwB65w5hI8a1eMmt9vabPx6/a95a+dbJEYk8vTUp5k1cNY5tUNtUxsvLtvNG6v343IbrpuQzv35g+kTF3EeI1dKqe5Nk7ULXf0eWP0ibPgzdDTDoDyY92vIuAyCtBTP8YwxOD7dgr2oCPu77+KsqSEoKoqYGTOInTuHqIkTkeCe9yfhNm7e2fUOz617Dlu7jS8P+zL3jrmXmNCzn5PX0NLOy8v38KdP9uHocPGFcWk8cGkm/RK18oBSSp2pnvfN1B0YA5UlnqHObUUQFAyjroWL7oOUkf6OLuC079uHrWgh9qIi2vftQ0JCiLpkGtY5c4jOyyMovOf2PO6o38GTJU+y4cgGxiSP4fuTvs+QhLOf09jo6OCVj/fyyoq9NLU7mZvdh29MzyQjOfo8Rq2UUj2LJmsXEpcTthV4Fg0cKIXwOJj6MOTeDTHnXiy7O3HW1GB/7z1shUU4Nm/2rOScMIGEO24nduZMLFarv0P0q6b2Jl7Y+AJ/3f5XYkNjeWLyE8wfPJ8gObve2JZ2J699sp+Xlu+moaWDy0f05qEZWQxN0TqySil1rjRZuxA47LDhDc9wp60CEgbBrGdgzJcgNMrf0QUMV1MTjUuWYi8spHn1anC7CRs+jF6PPkrs7FmEpGhCa4zhvb3v8YvSX1DXWse1WdfywLgHsIadXfLq6HDxvyUV/G7ZLmqb2skfkszDM4YwKq1nJ8NKKXU+abIWyBoqoeRFWP86tNmh/8Vw5dOQdQUEdc9d8s+Uu72d5uXLPTU5i4sxbW2EpKWRePddWOfOJSxDy2Z9Zk/DHp4seZI11WsYkTiC31z6G0Ymnd2webvTzd9LK1nw4S6q7Q4mZyTy0s1ZjO+vFTCUUup802QtEB1Y7xnq3PK25/6Iqz3z0fqO829cAcK43bSsLcVeVIh90WLcdjuWhATirrmG2DmziRgzpset5DyVlo4WXix7kTe2vEFkSCQ/mPQDvpj5RSxnkfA7XW7+teEAv/6gnKqjrYzvH8+z149mcoZuC6OUUp1Fk7VA4XZ5Nq/9ZAFUfAJhsXDRvZD7VYhL93d0fmeMoW37dk9NzoULcR4+jERGEjP9Mqxz5xI1aRISorvfH88Yw9KKpfx87c+pbq7mqsFX8dD4h0gIP/PeL7fbUFh2kF8tLWdvbTOj+lr5yVUjuSQrWRNjpZTqZJqs+Vt7M2z8C6z+rWcbDms/uPynMPYmCNfJ2e2VldgXLsRWWET77t0QHEz01KnEfutRYi69lKAI3a/rRCrsFTy15ilWHlhJVnwWP5/2c8b2GnvG1zHGsGhLNc8u2cnOw00MTYnh5ZvHM2N4b03SlFKqi2iy5i+N1Z6C6qWveja07ZsD1/4PDJ0Llp79a3HW1WF/733sRUW0btwIQETOeFIef5yYy2cSHB/v5wgDl8Pp4JVPX+HVza8SYgnh2xO+zQ1DbyA46MzeU8YYincc4dklO/n0gJ1ByVH85saxzB6VSlCQJmlKKdWVenZW4A/Vn3rmo21+C9xOGDYHLvo6pOf26KLq7uZmGj/4AFthEc2ffAIuF2FDhpD8zYexzp5NSJ8+/g4x4C2vWs5TJU9xoOkAswbO4pGcR0iOTD6jaxhj+GR3Hc8s3sGGigbSEyL45bWjmT+mD8EW3WRZKaX8QZO1rmAM7Frq2cR2zzIIiYKc22HSPZ5tOHoo095O08qV2AuLaPzwQ4zDQUifPiTecQexc2YTnpXl7xAvCAebDvL0mqcprixmoHUgf5j5ByamTjzj66zdV88vF+9g9Z56Uq3hPHX1KK7NSSNEkzSllPIrTdY6U4cDyv7mmY9Wsx1iUmH64zD+NojomUN5xu2mdf16bEVFNL73Pi6bDUtcHNarr8I6Zw4RY8ciWibLJ+2udl7b8hovl72MiPDguAe5ZfgthFjObKHFpsoGfrlkJ8t31pAUHcbjc4dzQ24/wkN0exillAoEmqx1huZaWPsHWPN7aKmFlFFw9cueLTiCQ/0dnV84duzEXlSEbWERzoOHkIgIYi69lNi5c4iePBkJ7ZntcrZWHVzFUyVPsc++jxn9Z/BozqOkRqee0TW2HrTz7JKdLN12mPjIEL43ayg3TxpARKgmaUopFUg0WTufanZ45qNtehNcbZ7Nay+6DwZM7ZHz0ToOHMC28F3sRUW07dwJFgtRF0+m10MPeVZyRmn1hTN1uPkwvyj9BYv2LSI9Jp3fTf8dU/pOOaNr7DrSyHNLy1lYdoiY8GC+OSOLr0wZSHSYfhwopVQg0k/nc2UM7F3umY9WvhiCwz1loCbdC8k9b86V8+hRGhctwlZYROu6dQBEjB1L7x98n9grryQ4QXe4Pxsd7g7+su0v/Hbjb3EZF/eOuZfbR95OmCXM52vsr2vm+aXl/HvjASJCLHz90sHcOWUQ1kjdn04ppQKZJmtny9kOW/7lSdKqN0NUMuQ/5lk4ENWzdnN3Oxw0LVuG7Z0CmlasAKeT0MEZJD/4ILFzZhOalubvEC9opdWlPFnyJLsadjEtbRrfyf0O6TG+b5R8oKGVBR+W8/fSKoKDhDunDuKr0waRGO17oqeUUsp/NFk7Uy31sO5Pnj3SGg9B8jCYtwBGXQsh4f6Orst8VvLJVlhA4/uLcDc1EdyrFwm33IJ13lzChgzRTVPPUW1rLc+WPkvhnkL6RPXh+fznyU/P97ldj9gdvFC8i7+uqQTg5kn9uTcvg16xPed9qpRS3YEma75qb4GlP4QNf4aOFhiUD/MXQMZlPWo+WtuuXdjeKcBWVITz0CGCIiOJmTkT6/x5RObmIhadnH6uXG4Xf9vxNxZsWECrq5W7Rt3FXdl3ERHsW7WGuqY2Xlq+h9c+2YfTbbguJ437L82kb5xWe1BKqQuRJmu+ComAitUw4guemp29R/g7oi7jrKnBtnAh9oJCHFu3/v9CgW9+k5jLtOTT+bSpZhNPrn6SbfXbmJQ6ie9N/B4DrQN9eq6tpYPfr9jDqyv34uhwcdXYvnzjskz6J+pCDqWUupBpsuYrEbiruMeUgnK3tHgqCrxT4Kko4HYTPnIkvb/3XWJnzSI4qWfNy+tsDY4GfrX+V/yz/J/0iujFLy75BZf3v9ynIc9GRwd/XLmP36/YQ6PDyezsVB6ansngXjFdELlSSqnO1jMyj/OlmydqxuWiefVq7AUF2JcsxbS0eCoK3H0X1rlzCcvI8HeI3Y7buHm7/G2eW/8cTe1N3Dr8Vr425mtEhZy+N6y13cXrq/bx4ke7OdrSwYzhvXloehbD+8R2fuBKKaW6TKdmHyJyBfA8YAH+YIx5+nOPhwGvA+OBOuB6Y8w+ERkAbAN2eE9dbYy5pzNj7amMMbTt2IHtnQLsRUU4a2oIionBOns21nlziRg/XisKdJJtddv4SclPKKspY1yvcXx/0vfJjM887fMcHS7+uqaCF4p3U9vUxiVZyTw8I4vR6XFdELVSSqmu1mnJmohYgBeAGUAVsFZECowxW4877Q7gqDFmsIjcAPwMuN772G5jzJjOiq+n66iu9lQUeKeAtvJyCAkheto0rPPmEZ13CUFhuq1DZ7G321mwYQF/2/E34sLieGrKU8wZNOe0Q54dLjdvlVbxmw/LOWRzMHFgAr+7aRwTBujedUop1Z11Zs9aLrDLGLMHQETeBOYDxydr84HHvbf/ASwQ3e+h07iammhcvARbQQEtJSVgDBFjxpDyw/8h5oorCI7vmfVKu4q93c6SfUv49YZf09DWwPVDruf+sfcTG3rqYUuny82/Nx7k+Q92Ulnfyth+cTxz7WgmZyTq9ihKKdUDdGay1heoPO5+FTDxZOcYY5wiYgMSvY8NFJENgB34vjFmxedfQETuBu4G6Nev3/mNvpswHR00rVyJvaCQxg8/xDgchPTrR9J992GdN5dQbbdOdbDpIMWVxSyrXEZpdSlO4yQ7KZvfTf8dwxOHn/K5brehaPMhfrV0J3tqmhnRJ5Y/3jaSvCHJmqQppVQP0pnJ2om+TYyP5xwC+hlj6kRkPPBvERlhjLH/x4nGvAy8DJCTk/P5a/dYxhgcn36KraAQ+8KFuOrrsVitxH3haqzz5hE+erR+2XcSYwxb67YeS9B2HPVMuxxkHcQtI24hPz2f7ORsguTk8wCNMSzeepjnluxke3UjWb2jefGm8Vw+orf+3pRSqgfqzGStCji+Jk4acPAk51SJSDBgBeqNMQZoAzDGrBOR3UAWUNqJ8V7w2qsOYC8swFZQSPvevUhoKNH5+VjnzyN6yhQkNNTfIXZL7a521lSvobiimGVVyzjScoQgCWJsr7E8kvMIeel59I/tf9rrGGNYtrOG55bspKzKxsCkKJ6/YQxzsvtgCdIkTSmleqrOTNbWApkiMhA4ANwAfOlz5xQAtwKrgGuAD40xRkSS8SRtLhEZBGQCezox1guWy2bD/v4ibIUFtJZ6CqdH5uSQcPtXiL38ciyxuo1DZ2hwNLDiwAqKK4tZeWAlLc4WIoIjmNJ3CnnpeUztO5X4cN/nAH6yu5ZfLt7Juv1HSYuP4BfXZHP12L4EW3QlrlJK9XSdlqx556DdDyzCs3XHq8aYLSLyBFBqjCkAXgHeEJFdQD2ehA5gGvCEiDgBF3CPMaa+s2K90Jj2dpqWL/cUTl+2DNPRQeigQSQ/+CDWuXMI6dvX3yF2SxX2imPDmxuObMBlXCRHJDN70Gzy0/PJTc0lzHJmq2jX7a/nmUU7WbWnjpTYcH5y1Uiuy0knNFiTNKWUUh7iGXG88OXk5JjS0u47SmqMoXXDRmwF79D43vu4bDYsiYnEzp6Fdd58wkcM1/lM55nbuNlcu9kzvFm5jN223QBkxWeRl55Hfno+wxOHn3L+2clsrrLxyyU7WLajhqToUO7NG8yXJvYjPERrqyqlVE8gIuuMMTm+nNu9t+TvBtr37cNWUIitsJCOykokPJyY6dOxzptL1OTJSLD+Cs8nh9PB6kOrWVa5jGWVy6hz1GERCzm9c7h2yLVcknYJaTFpZ3VtYwxbDtr59QflLN56mLjIEL59xVBundyfyFD9PSqllDox/YYIQM6jR7G/+y62ggIcm8pAhKiLJpF0773EzJiBJVoLc59Pda11LK9aTnFlMasOrsLhchAdEs2UvlPIT8/n4r4XYw2zntW165vbWbmrlhXlNawor+WQzUFMWDAPTc/i9ikDiAkPOc8/jVJKqe5Gk7UA4XY4aFq2zDMPbcUKcDoJGzKEXo8+Suyc2YT07u3vELsNYwx77XuPDW9uqtmEwZAalcrVmVeTl57HhN4TCLGceSLV7nSzvuLoseRs8wEbxkBseDBTMpP4+uBkZo1KIS5SV+YqpZTyjSZrfmTcblrWlmIrLKDx/UW4m5oI7tWLhFtv8eyHNmSIv0PsNlxuFxtrNrKschnFlcXst+8HYFjCML425mvkp+czJH7IGc/7M8awu6aZFeU1fFxey6o9dbS0uwgOEsb1i+fh6VlMzUpmVF+rbr+hlFLqrGiy5gdtu3Z55qEVFeI8eIigyEhiZs7EOn8ekbm5iEUnmZ8PLR0tfHLwE4ori1letZyGtgaCg4KZmDKRm4fdzCXpl5ASlXLG1z3a3M7K3bWs2OkZ3jxocwAwMCmKa8anMTUzmUmDEnSIUyml1HmhyVoXcdbUYFu4EHtBIY6tW8FiIeriyfR6+JvEXHYpQRER/g6xWzjScoSPqj6iuKKYkkMltLvbiQ2NZVraNPLS87i4z8VEh0af0TXbnW42VBxlRbknOSvzDm3GhAczZXAS91+azNTMJNITIjvpp1JKKdWTabLWidwtLTR+8AG2gkKaV64Et5vwESPo/b3vEjtrFsFJSf4O8YJnjKG8odwzvFlRzKd1nwKQFp3G9UOvJz89n7G9xhIc5Ptb3RjD3trmY8nZqt11NLe7sAQJY9PjePCyLKZmJZHd16qb1iqllOp0mqydZ8blonn1auwFBdiXLMW0tBDcJ5XEu+7COm8uYRkZ/g7xgtfh7mD94fXH5p8daDoAQHZSNt8Y9w3y0vLIiMs4o/lntpYOz9BmeQ3Ld9ZyoKEVgP6JkVw9ri9TM5O5KCORWB3aVEop1cU0WTtPHNu3Y3unAHtREc6aGoJiYrDOnoV13jwixo9HgrQH5lw0tTfx8cGPKa4oZsWBFTS2NxJmCWNS6iTuHHUnl6RdQnJkss/X63C52VjZwIqdNSwvr6WsqgG3gZiwYCYPTuRreRlMzUyif6Juk6KUUsq/NFk7Bx3V1diLirAVFNK2cycEBxM9bRrWefOIzs8jKOzMSg+p/3So6RDLqjzDm2sPr8XpdhIfFs9l/S4jLz2Pi1IvIjLEt3lixhj217V4es7Ka1m1u46mNidBAmPS4/j6pZlMy0pidFqcDm0qpZQKKJqsnSFXUxONi5dgKyigpaQEjCFizBh6/88PiL3ySoLjfS/erf6TMYbt9duP1d/cVr8NgAGxA7h5+M3kp+eTnZSNJci31bK21g5W7a5luXfuWWW9Z2gzLT6CeWP6MC0ziYsykrBG6NCmUkqpwKXJmo/cDgeHvvcYjR9+iHE4COnXj6T77sM6dw6h/fv7O7wLVoerg7XVa/mw8kOWVS7jcMthgiSIMcljeHj8w+Sl5zHQOtCnazldbjZVNbDcu6XGxkrP0GZ0WDCTMxK5e+ogpmYm0z8xUuuoKqWUumBosuajoPBwnHV1xH3hamLnziVizBj9wj9LtjYbKw6sYFnlMj4+8DHNHc1EBEcwuc9k7k+/n2lp00gIT/DpWhV1LSwvr2FFeQ2f7Kqj0Tu0OTo9jvsvzWRaZhKj0+MI0aFNpZRSFyhN1s5A/9f+5O8QLliVjZXHiqOvO7wOl3GRFJHElQOvJD89n9yUXMKDw097Hbujg0921XkqBuyqZX9dCwB94yKYM9oztDk5IwlrpA5tKqWU6h40WVOdwm3cbKndQnFlMcWVxexq2AXA4LjB3D7ydvLT8xmRNIIgOXWPl2do03as1ubGygZcbkNUqIWLMpK4Y8pApmYmM0CHNpVSSnVTmqyp86bV2cra6rUUVxbzUeVH1LTWYBEL43uP51sTvkVeeh7pMemnvU5lvXdoc2ctK3fX0uhwIgLZaXHcm5fB1MxkxvbToU2llFI9gyZr6qy1u9opqyljTfUaSg6VUFZbhtPtJCokiov7XEx+v3ym9p2KNcx6yus0OjpYtbvuWMWAfd6hzT7WcGaPSmVqZjKTMxKJjwrtih9LKaWUCiiarCmfOd1OttVto6S6hDWH1rDhyAYcLgdBEsTwhOHcMvwWJqZOJKd3DqGWkydWLrehrKrhWHK2vsIztBkZauGiQYncNnkAU7OSGZQUpUObSimlejxN1tRJuY2b8qPlrKlew5pDayg9XEpTRxMAmfGZXJN1DbkpuYxPGU9saOwpr1V1tOVYcrZyVx221g5EYFRfK/dc4tlSY1y/eEKDdWhTKaWUOp4ma+oYYwz77PtYc2gNJdUllFaXcrTtKAD9Y/tz5cAryU3NZULvCSRGJJ7yWk1tTlbvrju2MGBPbTMAqdZwLh/Rm6mZyVw8OIkEHdpUSimlTkmTtR7uYNNBSg6VHOs9O9J6BICUqBSmpU0jNzWX3JRcUqJSTnmdpjYnnx6wUbqvnuXltazffxSn2xARYmHSoARumtSfaVlJZCRH69CmUkopdQY0WethaltrWXNozbFFAVVNVQAkhCcwMWUiuam5TEyZSFpM2kmTqnanm+3VdjZV2dhU2UBZVQPlR5owxvP4qL5W7p42iCmZSYzvH09YsG/loZRSSin13zRZ6+ZsbTbWVq891nu2x7YHgJjQGCb0nsBNw29iYspEMuIyTpicud2GPbXNx5KyjVU2th200+5yA5AYFcro9Dhmj+pDdrqVMWlxumpTKaWUOo80WetmmjuaWXd43bHes+312zEYIoIjGNd7HFcNvorc1FyGxg/9r4LoxhgO2RyepKzSRllVA5urbDS2OQGICrUwKs3KVy4ewOj0OLLTrPSNi9BhTaWUUqoTabJ2gXM4HWyq2XSs5+zT2k9xGRehQaGM6TWG+8bcR25qLiMTRxJi+c8STA0t7WyqslFW2cCmqgY2VdmoaWwDIMQiDEuN5aqxfclOszImPY5BydFYgjQxU0oppbqSJmsXmA53B1tqtxxLzjYe2Ui7ux2LWBiZNJLbR97OxNSJjE4e/R+1NlvbXWyqrGdjZQNlVTY2VTUcq6spAoOSopiamcSY9Diy0+IYlhqjc82UUkqpAKDJWoBzuV1sP7r92HYa6w+vp9XZiiAMTRjKjUNvJDc1l/G9xxMVEgVAh8vNzsONlFUdYVOlp8ds5+FGXG7PCoA+1nCy0+K4YUI/RqdZGZlmJTZcC58rpZRSgUiTtQBjjGF3w+5jVQLWHl5LY3sjAIOsg5ifMf9YlYC48DiMMeyva2HplgY2Ve5nU1UDWw7acHR4FgBYI0IYnR7H9GG9GJ0WR3a6lV4x4acKQSmllFIBRJM1PzPGUNlYeWyfs5LqEuod9QCkRacxs/9MclNymZAygeTIZI7YHWyqsvGHjw6zqWoHZVU2bK0dAISHBDGyj5UvT+x/bJ5Zv4RIXQCglFJKXcA0WfOD6ubqY/ucraleQ3VzNQC9Inoxuc9kclNyyU3NJSa4F5urbGyqaOCdlfspqyrjkM0BgCVIGNI7hlmjUjw9ZmlxZPWOJtii5ZqUUkqp7kSTtS5Q11rH2sNrj22nsd++H4C4sDgmpEzgzpF3Mjo5h+ameDZX2Vi+3sZvqsrZU7Px2DUGJEaSOzCB7LQ4xqRbGZ5qJSJUFwAopZRS3Z0ma53A3m6ntLrUsxltdQnlR8sBiA6JZnzv8VyTeR2poSNpaEik7ICdNzY38NihHTi9CwCSY8IYnRbHF8b2JTvNs59ZXKRuNKuUUkr1RJqsnQctHS1sOLLh2KKAbfXbcBs34ZZwxvYay8VDZxDpHkJtXTKbdzXx849stLQfAg4RExZMdrqVu6YNYnRaHKPTraTEhus8M6WUUkoBmqydlXZXO5tqNh1bFFBWW4bT7SQ4KJjhCaOY2efLBDkyOVTTi9KSFhY3twOthAZXMTw1luty0slOszI6PY6BiVEE6UazSimllDoJTdZ81OHq4LWtr1FyqIQNRzbQ5mojSILoF5VFdvQ82psGsf9gL1Zu8WyZIQJZvVxcNrQXo9PjGJ0Wx5CUGEKDdQGAUkoppXzXqcmaiFwBPA9YgD8YY57+3ONhwOvAeKAOuN4Ys8/72HeBOwAX8IAxZlFnxno6QVh4Y8tfsBBFiuRhb+xP5aFUNrsiAEiLj2Bsehy3XWRldFocI/taiQrTXFgppZRS56bTsgkRsQAvADOAKmCtiBQYY7Yed9odwFFjzGARuQH4GXC9iAwHbgBGAH2ApSKSZYxxdVa8pyVC7bYHaW6zkBgVSnaalS/kx3m3zW+tYAMAAAqUSURBVLCSGB3mt9CUUkop1X11ZtdPLrDLGLMHQETeBOYDxydr84HHvbf/ASwQz8z6+cCbxpg2YK+I7PJeb1UnxntKliDh97dMJj0+krT4CF0AoJRSSqku0ZkTqPoClcfdr/IeO+E5xhgnYAMSfXwuInK3iJSKSGlNTc15DP3EJmckka4VAZRSSinVhTozWTtRRmN8PMeX52KMedkYk2OMyUlOTj6LEJVSSimlAltnJmtVQPpx99OAgyc7R0SCAStQ7+NzlVJKKaW6vc5M1tYCmSIyUERC8SwYKPjcOQXArd7b1wAfGmOM9/gNIhImIgOBTGBNJ8aqlFJKKRWQOm2BgTHGKSL3A4vwbN3xqjFmi4g8AZQaYwqAV4A3vAsI6vEkdHjP+zuexQhO4D6/rgRVSimllPIT8XRkXfhycnJMaWmpv8NQSimllDotEVlnjMnx5VzdTl8ppZRSKoBpsqaUUkopFcA0WVNKKaWUCmCarCmllFJKBTBN1pRSSimlApgma0oppZRSAUyTNaWUUkqpAKbJmlJKKaVUANNkTSmllFIqgGmyppRSSikVwLpNuSkRqQH2d8FLJQG1XfA6Fyptn9PTNjo1bZ/T0zY6NW2f09M2OrWuaJ/+xphkX07sNslaVxGRUl9refVE2j6np210ato+p6dtdGraPqenbXRqgdY+OgyqlFJKKRXANFlTSimllApgmqyduZf9HUCA0/Y5PW2jU9P2OT1to1PT9jk9baNTC6j20TlrSimllFIBTHvWlFJKKaUCmCZrSimllFIBTJO1kxCRK0Rkh4jsEpHvnODxaSKyXkScInKNP2L0Jx/a52ER2SoiZSLygYj090ec/uJD+9wjIptFZKOIfCwiw/0Rpz+dro2OO+8aETEiEjDL6LuCD++h20Skxvse2igid/ojTn/y5T0kItd5P4u2iMhfujpGf/LhPfTcce+fnSLS4I84/cmHNuonIsUissH7fTbLH3FijNH/PvcfYAF2A4OAUGATMPxz5wwAsoHXgWv8HXMAtk8+EOm9/TXgb/6OO8DaJ/a42/OA9/0dd6C1kfe8GGA5sBrI8XfcgdQ+wG3AAn/HGuBtlAlsAOK993v5O+5Aap/Pnf914FV/xx1obYRnocHXvLeHA/v8Eav2rJ1YLrDLGLPHGNMOvAnMP/4EY8w+Y0wZ4PZHgH7mS/sUG2NavHdXA2ldHKM/+dI+9uPuRgE9baXPadvI68fAzwFHVwYXAHxtn57Mlza6C3jBGHMUwBhzpItj9KczfQ/dCPy1SyILHL60kQFivbetwMEujO8YTdZOrC9Qedz9Ku8x5XGm7XMH8F6nRhRYfGofEblPRHbjSUYe6KLYAsVp20hExgLpxpiirgwsQPj6N/ZF79DMP0QkvWtCCxi+tFEWkCUiK0VktYhc0WXR+Z/Pn9PeaSoDgQ+7IK5A4ksbPQ7cJCJVwLt4eiC7nCZrJyYnONbTej5Oxef2EZGbgBzgF50aUWDxqX2MMS8YYzKAbwPf7/SoAssp20hEgoDngG92WUSBxZf3UCEwwBiTDSwFXuv0qAKLL20UjGcoNA9Pz9EfRCSuk+MKFGfyPXYD8A9jjKsT4wlEvrTRjcCfjDFpwCzgDe/nU5fSZO3EqoDj/5Wahp+6PgOUT+0jItOBx4B5xpi2LootEJzp++dN4KpOjSjwnK6NYoCRwDIR2QdMAgp60CKD076HjDF1x/1d/R4Y30WxBQpf/s6qgHeMMR3GmL3ADjzJW09wJp9DN9DzhkDBtza6A/g7gDFmFRCOp8h7l9Jk7cTWApkiMlBEQvG8kQv8HFMgOW37eIewXsKTqPWkeSLgW/sc/4UxGyjvwvgCwSnbyBhjM8YkGWMGGGMG4Jn3OM8YU+qfcLucL++h1OPuzgO2dWF8gcCXz+l/41nshIgk4RkW3dOlUfqPT99jIjIEiAdWdXF8gcCXNqoALgMQkWF4krWaLo0STdZOyBjjBO4HFuH5APy7MWaLiDwhIvMARGSCdwz7WuAlEdniv4i7li/tg2fYMxp4y7ssvMckuz62z/3erQQ2Ag8Dt/opXL/wsY16LB/b5wHve2gTnjmPt/knWv/wsY0WAXUishUoBh41xtT5J+KudQZ/YzcCbxrvcseexMc2+iZwl/fv7K/Abf5oKy03pZRSSikVwLRnTSmllFIqgGmyppRSSikVwDRZU0oppZQKYJqsKaWUUkoFME3WlFJKKaUCmCZrSqmAICJXi4gRkaFd+JoPikhkV73e2RCRPBGZ7O84lFL+o8maUipQ3Ah8jGdjyq7yIOD3ZE1Egk/xcB6gyZpSPZgma0opvxORaOBiPKVdbvjcY98Skc0isklEnvYeGywiS73H1otIhvf4oyKy1lvc/EfeYwNEZLuIvHZc0fNIEXkA6AMUi0ix99zfiUipd7PZHx0Xwz4R+ZH3tTZ/1vsnItEi8kfvsTIR+aL3+EwRWeU9/y3vz/f5n3mZiDwlIh8B3xCRuSJSIiIbvD9bbxEZANwDPOTdXHqqiCSLyD+9P+daEbn4/P42lFKB5lT/mlNKqa5yFfC+MWaniNSLyDhjzHoRudL72ERjTIuIJHjP/1/gaWPM2yISDgSJyEw8dR9z8RRoLhCRaXjKxQwB7jDGrBSRV4F7jTHPiMjDQL4xptZ73ceMMfUiYgE+EJFsY0yZ97FaY8w4EbkXeAS4E/gBYDPGjAIQkXhvWaPvA9ONMc0i8m08VSqeOMHPHWeMueSz5wKTjDFGRO4EvmWM+aaIvAg0GWOe8Z73F+A5Y8zHItIPz+7rw86x/ZVSAUyTNaVUILgR+JX39pve++uB6cAfjTEtAN5EKgboa4x523vMAZ7eLGAmsMF7nWg8yVsFUGmMWek9/mc85ZmeOUEc14nI3Xg+G1OB4cBnydq/vP9fB3zBe3s6x/UEGmOOisgc7/NWighAKCevu/i3426nAX8TT83PUGDvSZ4zHRjuvTZArIjEGGMaT3K+UuoCp8maUsqvRCQRuBQYKSIGsABGRL6Fp4fs8zXxhBMT4KfGmJc+d/0BJ7jGf9XZE5GBeHrMJniTrj/hKdr8mTbv/138/2fnyeJbYoy58SRxHq/5uNu/AZ41xhSISB7w+EmeEwRcZIxp9eH6SqluQOesKaX87RrgdWNMf2PMAGNMOp5epSnAYuD2z1ZsikiCMcYOVInIVd5jYd7HF3nPjfYe7ysivbyv0U9ELvLe/mwhA0AjEOO9HYsnebKJSG/gSh9iX4ynEDTe14wHVgMXi8hg77FIEcny4VpW4ID39q3HHT8+xhO95hgfrq2UuoBpsqaU8rcbgbc/d+yfwJeMMe8DBUCpiGzE0/MFcDPwgIiUAZ8AKcaYxcBfgFUishn4B/+f5GwDbvWenwD8znv8ZeA9ESk2xmzCM4S6BXgV+GzY9FR+AsSLyKcisgnP/Lca4Dbgr97XWw34sh3J48BbIrICqD3ueCFw9WcLDPAM4eZ4FzRsxbMAQSnVjYkx/zUaoJRS3YZ3GLTIGDPSz6EopdRZ0Z41pZRSSqkApj1rSimllFIBTHvWlFJKKaUCmCZrSimllFIBTJM1pZRSSqkApsmaUkoppVQA02RNKaWUUiqA/R/qCK/tzPfzdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00197628 0.05670447 0.0034227  0.02964427]\n",
      " [0.01910408 0.09335408 0.02452938 0.05928854]\n",
      " [0.04150198 0.14629259 0.04962921 0.08827404]\n",
      " [0.08366271 0.19167265 0.09925841 0.12384717]\n",
      " [0.13372859 0.25233309 0.15345123 0.17654809]\n",
      " [0.18972332 0.29859214 0.21791215 0.24505929]\n",
      " [0.256917   0.35540211 0.29207074 0.31027668]\n",
      " [0.33135705 0.39316675 0.38049059 0.38076416]]\n"
     ]
    }
   ],
   "source": [
    "failure_rates = np.zeros((8, 4))\n",
    "\n",
    "for r in np.arange(1, 9):\n",
    "    failure_rates[r-1, 0] = contraction(test[test.decision_T==1], 'judgeID_J', 'decision_T',\n",
    "                                   'result_Y', 'B_prob_0_machine', 'acceptanceRate_R',  r / 10, True)\n",
    "\n",
    "    ## Human error rate - Jotain väärin viel'\n",
    "    # Get judges with correct leniency as list\n",
    "    correct_leniency_list = test.judgeID_J[test['acceptanceRate_R'].round(1) == r / 10]\n",
    "    \n",
    "    # Released are the people they judged and released, T = 1\n",
    "    released = test[test.judgeID_J.isin(correct_leniency_list) & (test.decision_T == 1)]\n",
    "    \n",
    "    #print(released.shape, correct_leniency_list.shape[0], released.result_Y.value_counts())\n",
    "    # Get their failure rate, aka ratio of reoffenders to number of people judged in total\n",
    "    failure_rates[r-1, 1] = np.sum(released.result_Y == 0) / correct_leniency_list.shape[0]\n",
    "\n",
    "    ## True evaluation\n",
    "    failure_rates[r-1, 2] = contraction(test, 'judgeID_J', 'decision_T',\n",
    "                                   'result_Y', 'B_prob_0_machine', 'acceptanceRate_R',  r / 10, True)\n",
    "    ## Dec tree\n",
    "    failure_rates[r-1, 3] = contraction(test[test.decision_T==1], 'judgeID_J', 'decision_T',\n",
    "                                   'result_Y', 'B_prob_0_tree', 'acceptanceRate_R',  r / 10, True)\n",
    "    \n",
    "\n",
    "# klassifikaatioille scity.stats semin kautta error bati xerr ja yerr argumentttien kautta    \n",
    "\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.plot(np.arange(0.1,0.9,.1), failure_rates[:,0], label = 'Contraction')\n",
    "plt.plot(np.arange(0.1,0.9,.1), failure_rates[:,1], label = 'Human')\n",
    "plt.plot(np.arange(0.1,0.9,.1), failure_rates[:,2], label = 'True Evaluation')\n",
    "plt.plot(np.arange(0.1,0.9,.1), failure_rates[:,3], label = 'Classification ree')\n",
    "\n",
    "plt.title('Failure rate vs. Acceptance rate')\n",
    "plt.xlabel('Acceptance rate')\n",
    "plt.ylabel('Failure rate')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "print(failure_rates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of our model\n",
    "\n",
    "As defined in somewhere, our model is defined by the probabilistic expression $$ P(Y=0 | \\text{do}(R=r)) = \\sum_x P(Y=0|R=r, X=x, T=1)P(T=1|R=r, X=x)P(X=x)$$\n",
    "\n",
    "we approximate $P(X=x)$ with <!--frequencies and **Gaussian kernel**.-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "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",