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",
    "\n",
    "## COMPAS data\n",
    "\n",
    "*Following data filtering follows similar procedures as in the COMPAS analysis (link TBA)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7214, 53)"
      ]
     },
     "execution_count": 79,
     "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",
    "%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": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6172, 13)"
      ]
     },
     "execution_count": 80,
     "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": 81,
   "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": 81,
     "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": 82,
   "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": 83,
   "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": 84,
   "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": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Low       3421\n",
       "Medium    1607\n",
       "High      1144\n",
       "Name: score_text, dtype: int64"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compas.score_text.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "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": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tab = compas.groupby(['sex', 'race']).size()\n",
    "tab.unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "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": 88,
   "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": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x20c43c2b780>"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "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))"
   ]
  },
  {
   "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": 89,
   "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": 90,
   "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": 90,
     "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": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(25000, 8)\n",
      "(25000, 8)\n",
      "(12059, 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": 91,
     "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",
    "\n",
    "print(train_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": 111,
   "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).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGtCAYAAABJFTPQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4VNXWx/HvTgdCC50ECL330JGOIgKKqKCgoNKbNBEREbHRpHdQ5FJEUKooSJdOEnpCCzWhhJBASG+z3z8GeCMgBJnJmWTW53nuQ87MKb8MxKx79j57Ka01QgghhBDCNjkYHUAIIYQQQvw7KdaEEEIIIWyYFGtCCCGEEDZMijUhhBBCCBsmxZoQQgghhA2TYk0IIYQQwoZJsSaEEEIIYcOkWBNCCCGEsGFSrAkhhBBC2DAnowNYSt68ebW3t7fRMYQQQgghnsrf3/+W1jpfWvbNNMWat7c3fn5+RscQQgghhHgqpdTltO4rw6BCCCGEEDZMijUhhBBCCBsmxZoQQgghhA3LNHPWHicpKYmQkBDi4+ONjmI33Nzc8PLywtnZ2egoQgghRKaQqYu1kJAQsmfPjre3N0opo+NkelprwsPDCQkJoXjx4kbHEUIIITKFTD0MGh8fT548eaRQSydKKfLkySN3MoUQQggLsmqxppRqpZQ6o5QKUkqNeMJ+byiltFLKJ9Vrn9477oxS6qXnyPBfDxX/gXzeQgghhGVZbRhUKeUIzAJaAiGAr1JqvdY68KH9sgMDgYOpXqsAdAIqAoWBrUqpMlrrFGvlFUIIIYSwRda8s1YbCNJaX9BaJwIrgFcfs99XwAQg9djZq8AKrXWC1voiEHTvfBnSjRs36NSpEyVLlqRChQq0bt2as2fPPvN5pk6dSmxsrEUyXbp0ieXLlz/Y9vPzY+DAgRY5txBCCCEsx5rFmicQnGo75N5rDyilqgNFtNa/P+ux947vqZTyU0r5hYWFWSa1hWmtad++PU2aNOH8+fMEBgby7bffEhoa+sznelKxlpLybDcdHy7WfHx8mD59+jNnEkIIIYR1WbNYe9zkJf3gTaUcgCnA0Gc99sELWs/XWvtorX3y5UtTe610t2PHDpydnendu/eD16pVq0bDhg35+OOPqVSpEpUrV+aXX34BYOfOnTRp0oQ33niDcuXK0blzZ7TWTJ8+nWvXrtG0aVOaNm0KgLu7O6NHj6ZOnTrs37+fsWPHUqtWLSpVqkTPnj3R2vyRBQUF0aJFC6pWrUqNGjU4f/48I0aMYPfu3VSrVo0pU6awc+dO2rRpA0BERASvvfYaVapUoW7duhw/fhyAMWPG8MEHH9CkSRNKlCghxZ0QQgiRDqy5dEcIUCTVthdwLdV2dqASsPPepPSCwHqlVLs0HPvMvtwQQOC1u89zikdUKJyDL9pWfOI+J0+epGbNmo+8vnr1ao4ePcqxY8e4desWtWrVolGjRgAcOXKEgIAAChcuTIMGDdi7dy8DBw5k8uTJ7Nixg7x58wIQExNDpUqVGDt2rDlPhQqMHj0agHfffZfff/+dtm3b0rlzZ0aMGEH79u2Jj4/HZDIxbtw4Jk2axO+/m29q7ty580G2L774gurVq7N27Vq2b9/Oe++9x9GjRwE4ffo0O3bsICoqirJly9KnTx9ZU00IIYSwImveWfMFSiuliiulXDA/MLD+/pta60itdV6ttbfW2hs4ALTTWvvd26+TUspVKVUcKA0csmLWdLdnzx7efvttHB0dKVCgAI0bN8bX1xeA2rVr4+XlhYODA9WqVePSpUuPPYejoyMdOnR4sL1jxw7q1KlD5cqV2b59OwEBAURFRXH16lXat28PmBetzZo161OzvfvuuwA0a9aM8PBwIiMjAXjllVdwdXUlb9685M+f/z8N5wohhBAi7ax2Z01rnayU6g9sBhyBH7XWAUqpsYCf1nr9E44NUEqtBAKBZKDf8z4J+rQ7YNZSsWJFfv3110devz9E+Tiurq4PvnZ0dCQ5Ofmx+7m5ueHo6AiY15Tr27cvfn5+FClShDFjxhAfH//E6/ybxx1zf0mOtGYTQgghhGVYdZ01rfUfWusyWuuSWutv7r02+nGFmta6yb27ave3v7l3XFmt9Z/WzGlNzZo1IyEhgQULFjx4zdfXl9y5c/PLL7+QkpJCWFgYf//9N7VrP/mB1+zZsxMVFfXY9+4vRJs3b16io6MfFIg5cuTAy8uLtWvXApCQkEBsbOwTz9WoUSOWLVsGmIdH8+bNS44cOZ7tGxdCCCEyooQouHHC6BT/kKnbTdkCpRRr1qxh0KBBjBs3Djc3N7y9vZk6dSrR0dFUrVoVpRQTJkygYMGCnD59+l/P1bNnT15++WUKFSrEjh07/vFerly56NGjB5UrV8bb25tatWo9eG/JkiX06tWL0aNH4+zszKpVq6hSpQpOTk5UrVqVbt26Ub169Qf7jxkzhvfff58qVaqQNWtWFi9ebPkPRgghhLAlMbfg4Fw4NB+y5IYBh8HB0ehUAKj/Mkxmi3x8fLSfn98/Xjt16hTly5c3KJH9ks9dCCFEhnHnCuybAYeXQHI8lG8DDQaD16MPB1qSUspfa+3z9D3lzpoQQggh7FFoIOydBidWgVJQpRM0GAj5yhqd7BFSrAkhhBDCflw5CHumwNk/wTkb1OkN9fpCTi+jk/0rKdaEEEIIkblpDee2mIu0K/sgiwc0GQm1e0BWD6PTPZUUa0IIIYTInFKSIXCtuUgLPQk5vKDVOKjxHrhkMzpdmkmxJoQQQojMJSkOji6DvdPhzmXIWxZemwOV3gAnF6PTPTMp1oQQQgiROcTdAb8f4MAciAkDTx9o9R2UeRkcrLq0rFVJsWZl7u7uREdHP9j+6aef8PPzY+bMmQamEkIIITKRqFA4MBv8foSEu1CyOTQcDN4NzU96ZnBSrAkhhBAiYwo/b14j7ehyMCVBhdeg4SAoVNXoZBaVce8JZgLdunX7R99Qd3d3wNziqXHjxrz11luUKVOGESNGsGzZMmrXrk3lypU5f/48ABs2bKBOnTpUr16dFi1aPGiqPmbMGD744AOaNGlCiRIlmD59evp/c0IIIYS1XD8Gq96HmT7muWnV3oH+fvDmokxXqIE93Vn7c4Tle30VrAwvj3viLnFxcVSrVu3BdkREBO3atXvqqY8dO8apU6fw8PCgRIkSdO/enUOHDjFt2jRmzJjB1KlTadiwIQcOHEApxcKFC5kwYQLff/89AKdPn2bHjh1ERUVRtmxZ+vTpg7Oz8/N9v0IIIYRRtIbLe81PdgZtBZfsUH8A1O0L2Qsanc6q7KdYM0iWLFk4evTog+37c9aeplatWhQqVAiAkiVL8uKLLwJQuXLlB31BQ0JC6NixI9evXycxMZHixYs/OP6VV17B1dUVV1dX8ufPT2hoKF5etrvgnxBCCPFYJpN5Ads9UyDEF7Llg+ajwedDyJLL6HTpwn6KtafcATOCk5MTJpMJAK01iYmJD95zdXV98LWDg8ODbQcHB5KTkwEYMGAAQ4YMoV27duzcuZMxY8Y89nhHR8cHxwghhBAZQkqSuRXUnqlw6wzkKgavfA/VOoNzFqPTpSv7KdZskLe3N/7+/rz11lusW7eOpKSkZzo+MjIST09PABYvXmyNiEIIIUT6SowxN1XfNwPuhkD+ivD6QqjYHhzts2yxz+/aRvTo0YNXX32V2rVr07x5c7Jle7bVlMeMGcObb76Jp6cndevW5eLFi1ZKKoQQQlhZbAQcWgAH50JcBBStD22mQOmWmWL5jeehtNZGZ7AIHx8f/fBcsFOnTlG+fHmDEtkv+dyFEEKkWeRV2D8L/H+CpBgo08q8RlrRukYnsyqllL/W2ict+8qdNSGEEEKkv7CzsHcaHP8FtAkqvwENPoICFY1OhkmbcFC2s7qZFGtCCCGESD9X/c1Pdp76HZxcwed9qNcfchczOhlaa+Ycm8OlyEuMazTOZgo2KdaEEEIIYV1aw4WdsGcyXPwb3HLCC0OhTm9wz2d0OgBSTCl8ffBrfj37K6+Ves2m7q5JsSaEEEII6zClwKkN5jtp14+Ce0Fo+RXU7AZuOYxO90BCSgKf/P0J265so3vl7gysPhBlQw81SLEmhBBCCMtKToBjK8xz0iLOg0dJaDsdqnYyD33akLuJdxmwbQBHbh5hRO0RdC7f2ehIj5BiTQghhBCWkRAFfovMT3dG3zD36XxzMZRvCw6ORqd7RGhMKH229eFi5EUmNJpAq+KtjI70WFKsWVF4eDjNmzcH4MaNGzg6OpIvn3ls/tChQ7i4uBgZTwghhLCMmFvm9dEOzYf4SCjeCNrPgRJNbXaNtAuRF+i9pTeRCZHMbj6beoXrGR3pX0mxZkV58uR50Bd0zJgxuLu7M2zYsH/so7VGa42Dg21MYhRCCCHS7PZl2D/T3HEgOR7Kt4EGg8GrptHJnuh42HH6beuHg3JgUatFVMhTwehITyQVggGCgoKoVKkSvXv3pkaNGgQHB5Mr1/83o12xYgXdu3cHIDQ0lNdffx0fHx9q167NgQMHjIothBBCmIUGwuqeML26edizUgfodwg6LrX5Qm13yG66/9Wd7C7ZWfryUpsv1MCO7qyNPzSe0xGnLXrOch7l+KT2J//p2MDAQBYtWsTcuXOf2GR94MCBDB8+nLp163Lp0iXatGnDyZMn/2tkIYQQ4r+7ctD8ZOfZP8E5m3npjXp9IaeX0cnSZP359YzeO5oyucswu8Vs8mbJa3SkNLGbYs3WlCxZklq1aj11v61bt3LmzJkH27dv3yYuLo4sWbJYM54QQghhpjWc22Iu0q7sgywe0GQk1O4BWT2MTpcmWmsWBSxiiv8U6hSqw7Sm08jm/Gz9uI1kN8Xaf70DZi2pm7Y7ODiQukdrfHz8g6+11vIwghBCiPSXkgyBa81FWuhJyOEFrcZDjXfBJeMUOiZtYqLvRJaeWsrL3i/zdcOvcXHMWL9TZc6aDXBwcCB37tycO3cOk8nEmjVrHrzXokULZs2a9WD7/gMLQgghhFUkxYHvQphRA377EFKS4LU5MPAI1O2doQq1pJQkRuwewdJTS+lSvgvjGo3LcIUa2NGdNVs3fvx4WrVqRdGiRalQoQIJCQkAzJo1iz59+rBo0SKSk5Np2rTpP4o3IYQQwiLi7oDfD3BgDsSEgacPtPoOyrwMGXDFgpikGAbvGMz+6/sZVGMQH1T6wKa6EjwLlXr4LSPz8fHRfn5+/3jt1KlTlC9f3qBE9ks+dyGEyECibpgLNL8fIeEulGwODQeDd0ObXSPtacLjwum7rS9nIs4wpv4YXiv1mtGRHqGU8tda+6RlX7mzJoQQQtij8POwbwYcXQ6mJKjwGjQcZO46kIEFRwXTe0tvbsbeZHqz6TTyamR0pOcmxZoQQghhT64fNz80ELgWHJygWmeoPwDylDQ62XM7FX6KPlv7kKyTWfDiAqrlr2Z0JIvI9MWa1jrDjlFnRJllWF0IITKd5ATY/hXsmwku7uYCrW5fyF7Q6GQWcfD6QT7a8RHZXbLzY4sfKZGrhNGRLCZTF2tubm6Eh4eTJ08eKdjSgdaa8PBw3NzcjI4ihBAitZun4LceEHoCar4PLcZAllxPOyrD2HRpEyN3j6RYjmLMaTGHgtkyRwF6n1WLNaVUK2Aa4Ags1FqPe+j93kA/IAWIBnpqrQOVUt7AKeD+arAHtNa9n/X6Xl5ehISEEBYW9t+/CfFM3Nzc8PLKGCtZCyFEpqe1ubn6ltHmu2lvr4CyLxudyqKWn1rOuEPjqJ6/OtObTSena06jI1mc1Yo1pZQjMAtoCYQAvkqp9VrrwFS7Lddaz723fztgMtDq3nvntdbPNdjs7OxM8eLFn+cUQgghRMYUFQrr+kLQVij9Irw6C9zzG53KYrTWzDgygwUnFtC0SFMmNJqAm1PmHNmx5p212kCQ1voCgFJqBfAq8KBY01rfTbV/NkAmPAkhhBDP6/RGWD8AEmOg9SSo1T3DLsPxOMmmZMbuH8uaoDV0KN2BUXVH4eSQeWd2WfM78wSCU22HAHUe3kkp1Q8YArgAzVK9VVwpdQS4C4zSWu9+zLE9gZ4ARYsWtVxyIYQQIiNKjIFNn8LhxVCwCnRYCPnKGp3KouKS4xi+azg7Q3bSu2pv+lbtm+nnpVuzWHvcJ/fInTOt9SxgllLqHWAU0BW4DhTVWocrpWoCa5VSFR+6E4fWej4wH8yL4lr6GxBCCCEyjKv+5ocIIi5Ag0HQ9DNwynitlZ4kMiGS/tv6cyzsGKPqjKJjuY5GR0oX1izWQoAiqba9gGtP2H8FMAdAa50AJNz72l8pdR4oA/j9++FCCCGEHTKlwJ7JsHMcuBeErhug+AtGp7K4GzE36L2lN1eirjCp8SRe9H7R6EjpxprFmi9QWilVHLgKdALeSb2DUqq01vrcvc1XgHP3Xs8HRGitU5RSJYDSwAUrZhVCCCEyntuXYHUvCD4AlTrAK99DltxGp7K4oNtB9N7am5ikGOa1nEetgrWsdq3bMYlcvRNHJU/bearUasWa1jpZKdUf2Ix56Y4ftdYBSqmxgJ/Wej3QXynVAkgCbmMeAgVoBIxVSiVjXtajt9Y6wlpZhRBCiAxFazj+C2wcZn5w4PUFUPnNTPUQwX1Hbx6l37Z+uDi68FOrnyjrYb05eCevRtJriT8AO4Y1wcXJNhrYZ+pG7kIIIUSmE3cbfh8MAWugaD1oPw9yFzM6lVXsDN7JsF3DKJitIHNbzMUru/XW8Vx9OIRPV5/AI5sLc7vUpGoR6y4aLI3chRBCiMzo4t+wpjdEh0Kzz6HhYHBwNDqVVaw+t5qx+8dS3qM8s1rMwsPNwyrXSUw28c3GQBbvv0zdEh7MfKcGed1drXKt/0qKNSGEEMLWJSfA9q9h3wxzw/UPt4BnDaNTWYXWmgUnFjDjyAwaFG7A5CaTyeqc1SrXuhkVT79lh/G9dJseLxTnk1blcHK0jaHP1KRYE0IIIWzZzdOwujvcuNfX86VvwCWb0amsIsWUwnjf8fx8+mfalGjD2PpjcXZ0tsq1/C/fps9Sf6Lik5n+dnXaVS1sletYghRrQgghhC3SGg4tgC2fm4uzTj9DudZGp7KaxJREPt39KX9d/ouuFboyxGcIDsryd7m01iw9eIWxGwIonCsLiz+oTflCOSx+HUuSYk0IIYSwNVGhsK4fBG2BUi3NfT2zFzA6ldVEJ0YzaMcgDt44yDCfYXSt2PXpB/0H8UkpjFp7kl/9Q2haNh9TO1YnZ1br3LmzJCnWhBBCCFty+g9Y3z/T9vV82K24W/TZ2oeg20F82/Bb2pZsa5XrXL0TR+8l/py4GsnA5qUZ1Lw0Dg4Z43OVYk0IIYSwBYkxsHkk+P8EBSvD6wshfzmjU1nV5buX6bWlFxHxEcxoPoOGng2tcp29QbcY8PMRkpJNLHzPhxYVMtZdSinWhBBCCKNdPQyre0D4eWjw0b2+nra1fISlBdwKoO+2vmit+eHFH6icr7LFr6G1Zv7fFxi/6TQl87kz792alMjnbvHrWJsUa0IIIYRR/tHXswB0XQ/FGxmdyur2XdvHoB2D8HDzYG6LuXjn9Lb4NWISkhn+23E2Hr9O68oFmfhGVbK5ZsyyJ2OmFkIIITK625dhTS+4sh8qvg5tJmfKvp4P23hhI6P2jqJEzhLMbTGXfFnzWfwaF2/F0GuJH0E3o/n05XL0bFQClYHn/UmxJoQQQqQnreH4SvhjmHm7/Xyo8lamfojgvv8F/I+JfhPxKeDD9GbTye6S3eLX2BoYyuBfjuLkqPjfB3VoWDqvxa+R3qRYE0IIIdJL3G34fQgErM70fT1T01oz5fAUFp1cRMtiLfnuhe9wdbTsnDyTSTN12zmmbztHJc8czO1SE6/c1ul8kN6kWBNCCCHSw8Xd9/p63sj0fT1TSzIlMWbfGNafX0/Hsh35tPanOFr4+46MS2LwL0fZfvomHWp48U37Srg5Z57PVoo1IYQQwpqSE2HH17B3OniUgA//As+aRqdKF7FJsQzdNZQ9V/fQr1o/elXpZfG5Y6dv3KXXEn+u3o7jq1cr0qVusQw9P+1xpFgTQgghrCXsDPz24b2+nt3gpW8zbV/Ph92Ov03/bf05GX6SL+p9wRtl3rD4NTYcu8bwX4/j7ubEL73qUrOYh8WvYQukWBNCCCEsTWvwXQh/jbKLvp4PuxZ9jV5benE95jqTm0ymedHmFj1/coqJ8ZtOs2D3RXyK5WZ25xrkz+Fm0WvYEinWhBBCCEuKvmnu63nuLyjVAl6dnan7ej7s7O2z9NnSh7iUOOa1nEfNApYd8g2PTqD/8iPsvxBO13rF+OyVCrg4Wb7huy2RYk0IIYSwlDN/wrr+kBgNL0+E2j3sYkmO+/xu+DFw+0CyOGdhcavFlM5d2qLnPxZ8hz5L/QmPSWTSm1V5o6aXRc9vq6RYE0IIIZ5XYgxs/gz8F0GBytAh8/f1fNi2y9sY/vdwPLN7Mq/FPAq5F7Lo+Vf6BjNq3UnyubvyW5/6VPLMadHz2zIp1oQQQojnkbqvZ/2B0GxUpu/r+bCVZ1byzcFvqJS3ErOazSKXWy6LnTshOYUvNwSy/OAVGpbKy/S3q+ORzcVi588IpFgTQggh/gtTCuyZAju/s6u+nqlprZl7bC6zj82mkVcjJjaaSFZnyy1EeyMynj7L/Dly5Q69G5fk45fK4uhgP8PK90mxJoQQQjyr25fNC9xe2QcV20ObKXbR1zO1FFMK3xz8hlVnV9GuZDvG1B+Ds4Ozxc5/6GIEfZcdJjYxmdmda9C6smWHVTMSKdaEEEKIZ3F8JWwcal6eo/08qNLRrh4iAEhISWDE3yPYemUrH1b6kI9qfGSxhWi11vy07xLfbDxFUY+s/NyjDqULWL6HaEYixZoQQgiRFnF3zEXayV+hSF14fR7k9jY6Vbq7m3iXgdsH4h/qzye1PqFLhS4WO3dcYgoj15xgzZGrtChfgMkdq5LDzXJ36zIqKdaEEEKIp/lHX89R0GAwONrfr9CbsTfpvbU3FyMvMqHRBF4u/rLFzh0cEUuvJf6cunGXoS3L0K9pKRzscH7a49jfvzQhhBAirey4r+fDLkZepPeW3txJuMPs5rOpV7iexc6962wYA38+gtaaH7vWomm5/BY7d2YgxZoQQgjxOGFn4LfucOM41Ohq7uvp6m50KkMcDztOv239cFAO/NjqRyrmqWiR82qtmb3zPJP+OkPZAtmZ925NiuWxj96pz0KKNSGEECK11H09nbNCp+VQ7hWjUxlmd8huhu4aSh63PMxvOZ8iOYpY5LxR8UkMW3WMzQGhtK1amPEdKpPVRcqSx5FPRQghhLjvkb6esyB7QaNTGWb9+fV8sfcLSucuzewWs8mbJa9Fzht0M5peS/y4FB7LqFfK82HD4hZ7mjQzkmJNCCGEADizyVyo2Wlfz9S01vwU8BOT/SdTp1AdpjaZiruLZYaAN528wbBVx3B1cmDph3WoVzKPRc6bmUmxJoQQwr4lxsJfn4Hfj/f6ei6A/OWNTmUYkzYxyW8SSwKX8JL3S3zb8FtcHJ+/vVOKSTN5yxlm7ThPVa+czOlSk8K5slggceYnxZoQQgj7de0I/NYDwoOg/gBo9rnd9fVMLSkliVF7R/HHxT/oXL4zw2sNx0E5PPd578QmMnDFUf4+G0anWkUY064ibs6OFkhsH6RYE0IIYX9MKbB3Kuz4FrLlh/fWQYnGRqcyVExSDEN2DmHftX18VOMjPqz0oUXmkQVci6T3Un9CIxP47vXKvF27qAXS2hcp1oQQQtiXO1dgda//7+v5ymTI6mF0KkOFx4XTb1s/TkecZmz9sbQv3d4i51175CojVh8nVxYXVvSqS42i9tU/1VKkWBNCCGE/Uvf1fG0uVO1ktw8R3BccFUzvLb25GXuTaU2n0bjI899hTEox8e0fp1i09xK1i3sw650a5Mtuv8PLz0uKNSGEEJmf9PV8rNMRp+m9pTdJpiQWvLiAavmrPfc5b0bF03/ZEQ5diuCDBsX5tHU5nB2ff96bPbPqp6eUaqWUOqOUClJKjXjM+72VUieUUkeVUnuUUhVSvffpvePOKKVesmZOIYQQmdilPTC3IQSsgaajoNtGKdSAQ9cP0W1TN5wdnVny8hKLFGqHr9ym7Yw9HL96h6kdqzG6bQUp1CzAanfWlFKOwCygJRAC+Cql1mutA1PttlxrPffe/u2AyUCre0VbJ6AiUBjYqpQqo7VOsVZeIYQQmUxyIuz8FvZMBY/i8OEW8LLPvp4P23xpM5/u/pRiOYoxp8UcCmZ7voV/tdYsP3SFMesDKJjTjdV9GlChcA4LpRXWHAatDQRprS8AKKVWAK8CD4o1rfXdVPtnA/S9r18FVmitE4CLSqmge+fbb8W8QgghMouws7C6O1w/Zvd9PVPTWvPz6Z8Zd2gc1fJXY0azGeR0zflc54xPSuGLdQH84hdMozL5mN6pGrmyPv+6bOL/WbNY8wSCU22HAHUe3kkp1Q8YArgAzVIde+ChYz0fc2xPoCdA0aLyKLAQQti9B309PwfnLNBxGZRvY3Qqm3A64jQTfSdy6MYhmhRpwsRGE3Fzcnuuc167E0efpf4cC4mkf9NSDG5ZBkcH+35gwxqsWaw97m9LP/KC1rOAWUqpd4BRQNdnOHY+MB/Ax8fnkfeFEELYkdAA2DQCLv4NJZvDa7Ptuq/nfWGxYcw4MoO1QWvJ6ZqTkXVG8maZN3FyeL4SYN/5WwxYfoSEZBPz3q3JSxXls7YWaxZrIUCRVNtewLUn7L8CmPMfjxVCCGGvYsJhxzfgvwjccprXTav5PjjY98T2+OR4Fgcs5oeTP5BkSqJrxa70qNKDHC7PN5dMa80Pey7y3Z+n8c6TlXnv+lAqvwwxW5M1izVfoLRSqjhwFfMDA++k3kEpVVprfe7e5ivA/a/XA8uVUpMxP2BQGjhkxaxCCCEympQkOLQAdo2DhGio3RMaf2L3C9yatIk/Lv7BtMPTuBFzgxZFWzCk5hCK5Cjy9IOfIjYxmU9+O8GGY9doVbEgk96qirurrAJmbVb7hLXWyUqp/sBmwBH4UWsdoJQaC/hprdcD/ZUBGSqqAAAgAElEQVRSLYAk4DbmIVDu7bcS88MIyUA/eRJUCCHEA2f/gs0jIfycecjzpW8hfzmjUxnu6M2jTPCdwIlbJyjvUZ7vGn6HT0Efi5z70q0Yei3x59zNKIa3KkufxiUt0o5KPJ3SOnNM9fLx8dF+fn5GxxBCCGFNYWfMRVrQVshTylyklX7R7rsQXI2+yhT/KWy+tJn8WfIzsMZA2pZsa5Em7ADbT4fy0YqjODoopneqTqMy+SxyXnumlPLXWqepkpZ7l0IIIWxf3G3YOR58F4BzNnORVqsHONn3EhHRidEsPLGQJYFLcFAO9Knah24Vu5HVOatFzm8yaWZsD2LqtrOUL5iDee/WpIiHZc4t0k6KNSGEELYrJdn84MCObyH+jnnNtGajIFteo5MZKtmUzJqgNcw8MpOI+AjalWzHwOoDKZCtgMWucTc+iSG/HGXrqZu8Xt2Tb9pXJouLo8XOL9JOijUhhBC26fwO85DnzUDwfgFajYOClYxOZbh91/Yx0XciQXeCqJG/BrObz6Zi3ooWvcbZ0Ch6LfEnOCKWMW0r0LW+t8xPM5AUa0IIIWxL+HnzorZn7vXw7LgUyrWx+3lpF+5cYJLfJHZf3Y2XuxeTm0ymRdEWFi+iNh6/zse/HiOrixPLe9SldnH7frrWFkixJoQQwjbE34W/J8KBOeDkCs2/gLp9wfn5VtnP6G7H32bOsTmsPLOSLE5ZGFpzKO+UfwcXR8vO10tOMTFx8xnm/X2BGkVzMadLTQrksO/P3lZIsSaEEMJYphQ4shS2fwUxt6BaZ2j+ud13H0hKSWL56eXMOz6PmKQY3izzJn2r9cXDzfJ3uiJiEhnw82H2BoXTuU5RvmhbERcn+15U2JZIsSaEEMI4l/aaW0TdOA5F6sI7K8GzhtGpDKW1ZvuV7Uz2n8yVqCs08GzAsJrDKJW7lFWudyIkkt5L/QmLTmBChyq8Vev5F88VliXFmhBCiPR3+zJsGQ2BayGHF3T4ASp1sPt5aafCTzHRbyK+N3wpmbMkc1rMoaFnQ6tdb5VfMJ+tPUnebC782rseVbxyWe1a4r+TYk0IIUT6SYiGPVNg3wxQDtBkJNQfAC72vXbXzdibTD88nfXn15PLNRej6oyiQ5kOz91s/d8kJpv46vdAlhy4TP2SeZjxdnXyuLta5Vri+UmxJoQQwvpMJjixEraOgajrUPktaDEGcnoaHMxYcclx/BTwE4tOLiLZlEy3it3oUaUH2V2yW+2aoXfj6bPUn8NX7tCzUQmGv1QWJ0eZn2bLpFgTQghhXcG+sOkTuOoPhWvAW/+DIrWNTmUokzax8cJGph2eRmhsKC2LtWRwzcEUyW7d+WK+lyLou+wwMQnJzHynOm2qFLbq9YRlSLEmhBDCOiKvmu+knVgJ7gXhtblQpSM42PddnMOhh5noO5GT4SepmKci4xuNp2aBmla9ptaaJQcuM3ZDIF65s7D0wzqULWi9u3fCsqRYE0IIYVmJseY5aXunmpfleGEYNBwMru5GJzNUSFQIU/yn8Nflv8ifNT/fNvyWV0q8YrFm6/8mPimFkWtOsPrwVZqXy8/kjtXImcXZqtcUliXFmhBCCMvQGgJWw5YvIDIYKrwKLceauxDYsajEKBacWMDSwKU4OTjRt1pfulboarFm608SHBFL76X+BFy7y6AWpRnYrDQODvb9xG1GJMWaEEKI53ftCGz6FK7sh4KVof1c8LbekhMZQbIpmdXnVjPr6CyrNVt/kt3nwhj48xGSTZofuvrQvHz6XFdYnhRrQggh/ruoUNg2Fo4ug2x5oe10qN4FHByNTmaovVf3MslvEkF3gqhZoCazW8ymYh7LNlv/N1pr5u66wMTNpymV35157/pQPG+2dLm2sA4p1oQQQjy75AQ4MBv+nmT+un5/aPQxuOU0Opmhzt85zyS/Sey5ugcvdy+mNJlC86LNLd5s/d9EJyTz8apj/HnyBq9UKcSEDlXI5iq/6jM6+RsUQgiRdlrD6d/hr1Fw+xKUbQ0vfg15ShqdzFC3428z6+gsfj37K1mdsjLMZxhvl3vb4s3Wn+R8WDS9lvhzISyaz1qXp/sLxdOtSBTWJcWaEEKItLlx0tzH89JuyFce3l0DJZsZncpQiSmJLD+1nPnH5xObHPug2Xput9zpmuOvgBsMXXkMZycHln5Yh/ql8qbr9YV1SbEmhBDiyWJuwfav4fBi8zBn60lQ831wtN9fIVprtl3ZxmT/yQRHBfOC5wsM9RlKyVzpe4cxxaSZuvUsM7YHUcUrJ3O61MQzV5Z0zSCsz35/0oQQQjxZciL4LoCd4yExGmr3hMafQFYPo5MZKiA8gIm+E/EP9adUrlLMazGP+p710z1HZGwSH/1yhJ1nwnizphdfvVYJN2f7frAjs5JiTQghxD9pDef+gs0jITwISjaHVt9BvrJGJzNUaEwo049MZ8P5DeR2y83ndT/n9dKvW63Z+pOcun6XXkv8uR4Zx9evVaJznaIyPy0Tk2JNCCHE/ws7Y14v7fw2yFMK3lkJpV8EOy4EYpNiWRywmEUB95qtV+pGj8rWbbb+JOuOXuWT346TM4szK3rWo2ax9J0fJ9KfFGtCCCEgNgJ2jYdDC8DFHV76Fmr1AKf0e5rR1pi0id8v/M60w9O4GXuTF4u9yOCag/HK7mVInqQUE+P+PM0Pey5S29uDmZ2rkz+7myFZRPqSYk0IIexZSjL4L4Id30B8JNTsBk0/My9wa8f8Q/2Z4DuBwPBAKuapyMRGE6lRoIZhecKiEui//DAHL0bQrb43n71SHmdH6/YUFbZDijUhhLBX57fDppEQdgq8X4BW46BgJaNTGSo4Kpgp/lPYcnkLBbIWSLdm609y5Mpt+iw9zO3YRCa/VZXXaxhzZ08YR4o1IYSwN+HnzYvanvnD3GS941Io18au56VFJUYx//h8lp1ahpODE/2q9aNrxa5kcTJ2GYyfD13hi3UB5M/hym996lPJ0747RNgrKdaEEMJexEfC3xPhwFxwcoUWY6BuX/PXdirZlMxvZ39j1tFZ3Em4w6ulXmVA9QHkz5rf0FwJySmMWR/Az4eCeaF0XqZ3qk7ubPY7f9DeSbEmhBCZnSkFjiwxL2wbcwuqd4ZmoyF7AaOTGWrP1T1M8p3E+cjz+BTw4eNaH1MhTwWjY3E9Mo7eSw9zLPgOfZuUZOiLZXF0sN+7nkKKNSGEyNwu7TG3iLpxAorUhc6roHB1o1MZKuh2EJP8JrH32l6KZi/K1KZTaVakmU2sU3bgQjj9lh0mPimFuV1q0KpSIaMjCRsgxZoQQmRGty/Dls8hcB3kLAJv/AgVX7freWkR8RHMOjKLX8/9SjbnbHzs8zFvl3sbZ0dno6OhtebHvZf49o9TFMuTlV/erUup/Mas4yZsjxRrQgiRmSREw57JsG8mODial+GoPwCc7bdfZGJKIstOLWP+8fnEJcfRsWxH+lTtk+7N1v9NbGIyn64+wbqj12hZoQCT36pKdjfjC0hhO6RYE0KIzMBkguMrYOuXEH0DKr9lfoAgp6fRyQyjtWbL5S1M9p/M1eirNPJqxFCfoZTIWcLoaA9cDo+h1xJ/zoRG8fFLZenTuCQOMj9NPESKNSGEyOiCD8Gfn8C1w+BZ07wUR5FaRqcyVMCtACb4TuDwzcPmZust51G/cPo3W3+SHWdu8tHPR1BKsahbLZqUNfYJVGG7pFgTQoiMKjIEto6BE6sgeyFoP898R83Bfle2D4sNY4r/FDZc2ICHmwej642mfan2hjRb/zcmk2bWjiAmbz1LuYI5mNelJkXzZDU6lrBhtvOvVwghRNokxsK+6bBnKmgTvDAMGg4GV3ejkxnK94Yvw3YNIzoxmg8rfUj3yt1xd7Gtz+RufBJDVx5jS2Aor1YrzLjXq5DFxdHoWMLGWbVYU0q1AqYBjsBCrfW4h94fAnQHkoEw4AOt9eV776UAJ+7tekVr3c6aWYUQwuZpDQFr4K/P4W4IVHgNWo6F3MWMTmYorTXLTi1jkt8kimQvwo8v/UjJXCWNjvWIc6FR9Friz+WIWEa3qcD7DbxtYrkQYfusVqwppRyBWUBLIATwVUqt11oHptrtCOCjtY5VSvUBJgAd770Xp7WuZq18QgiRoSTGwMahcOxnKFgFXp8P3g2MTmW4uOQ4xuwbwx8X/6BpkaZ80/AbsrvY3pIXf564zrBVx8ji4siy7nWoWyKP0ZFEBmLNO2u1gSCt9QUApdQK4FXgQbGmtd6Rav8DQBcr5hFCiIzp5mlY1RXCzkCTT6HRx+ZlOexccFQwg3cM5uztswyoPoDulbsb2nD9cVJMmombzzB313mqFcnFnC41KJTTfpdREf+NNYs1TyA41XYIUOcJ+38I/Jlq200p5Yd5iHSc1nqt5SMKIYSNO/ozbBwCLtngvbVQoonRiWzCnqt7+OTvT9BoZjWfxQteLxgd6RG3YxIZuOIIu8/d4u3aRRnTrgKuTlJki2dnzWLtcQPx+rE7KtUF8AEap3q5qNb6mlKqBLBdKXVCa33+oeN6Aj0BihYtapnUQghhCxJj4c+P4chS8H4BOiyE7AWNTmU4kzax8MRCZh6ZSencpZnaZCpFchQxOtYjTl6NpNcSf8KiEhj3emU61ZbfUeK/s2axFgKk/gnyAq49vJNSqgXwGdBYa51w/3Wt9bV7f15QSu0EqgP/KNa01vOB+QA+Pj6PLQSFECLDCTtrHva8eco85Nl4BDjKw/vRidF8tucztgdvp3Xx1nxR7wuyOtvekhe/+Ycwcs0JPLK5sLJ3PaoVyWV0JJHBWfOn3xcorZQqDlwFOgHvpN5BKVUdmAe00lrfTPV6biBWa52glMoLNMD88IEQQmRux1fBho/A2Q26/AalmhudyCZcuHOBj3Z8RHBUMMNrDadL+S429yRlYrKJbzYGsnj/ZeqW8GDmOzXI6+5qdCyRCVitWNNaJyul+gObMS/d8aPWOkApNRbw01qvByYC7sCqez9095foKA/MU0qZAAfMc9YCH3shIYTIDJLiYNMI8P8JitaHN36AHIWNTmUTtl7eymd7PsPNyY0FLy6gVkHb685w8248/ZYfxvfSbbo3LM6Il8vh5GhbDzuIjEtpnTlGD318fLSfn5/RMYQQ4tmFn4eVXSH0hHlx26ajZNgTSDGlMPPoTBaeWEiVvFX4vsn3FMxme/P2/C9H0GfpYaLikxn/RhXaVZUiWzydUspfa+2Tln3lvwZCCGGkk6th/UBzcfbOKijzotGJbMKd+DsM/3s4+6/vp0PpDoysMxIXRxejY/2D1pqlB68wdkMAhXNlYfEHtSlfKIfRsUQmJMWaEEIYISke/voMfBeCV214cxHk9DI6lU04FX6KwTsHczP2Jl/U+4I3yrxhdKRHxCelMGrtSX71D6Fp2XxM7VidnFmdjY4lMikp1oQQIr1FXIBV3eD6Mag/AJp/AY7yix5gw/kNfLn/S3K55mJxq8VUzlfZ6EiPCLkdS5+lhzlxNZKBzUszqHlpHBxs62EHkblIsSaEEOkpcB2s6w/KAd5eAWVfNjqRTUgyJTHJdxLLTy/Hp4APkxpPIk8W22vJtDfoFgN+PkJSsomF7/nQokIBoyMJOyDFmhBCpIfkRNjyORycC5414c2fIJcslApwK+4WQ3cO5fDNw7xb4V0G1xyMs4Nt3WnUWjP/7wuM33SakvncmfduTUrkczc6lrATUqwJIYS13b5sHva8dhjq9oUWX4KTbU2WN8rRm0cZsnMIUYlRjH9hPK1LtDY60iNiEpIZ/utxNp64TuvKBZnwRlXcXeXXp0g/8q9NCCGs6fRGWNvH3Gyv41Io39boRDZBa82qs6v47tB3FMxakDmt51DWo6zRsR5x8VYMvZb4EXQzmk9fLkfPRiVsbjFekflJsSaEENaQkgRbx8D+mVComnnY06O40alsQkJKAl8f+Jq1QWtp6NmQcS+MI6drTqNjPWJrYCiDfzmKk6Pifx/UoWHpvEZHEnZKijUhhLC0O8Hw6/sQ4gu1e8KLX4OTtB0CuB59nUE7BxEYHkivKr3oU7UPjg6ORsf6B5NJM3XbOaZvO0clzxzM6VyTIh6214NU2A8p1oQQwpLOboY1vSAl2Xw3rWJ7oxPZjIPXD/Lxro9JNCUyrek0mhVtZnSkR0TGJTH4l6NsP32TDjW8+KZ9JdycbauYFPZHijUhhLCElCTY/hXsnQYFK8ObiyFPSaNT2QStNYsDFjPl8BS8c3gztelUiue0vSHh0zfu0muJP1dvx/HVqxXpUreYzE8TNiHNxZpSKpvWOsaaYYQQIkOKvAq/fgDBB8DnA3jpO3B2MzqVTYhNimX0vtFsvrSZlsVa8lWDr8jmnM3oWI/YcOwaw389jrubEyt61sXH28PoSEI88NRiTSlVH1gIuANFlVJVgV5a677WDieEEDbv3FZY0xOSE6DDD1DZ9lojGeXy3csM2jGIC5EXGFRjEB9U+sDm7lQlp5gYv+k0C3ZfxKdYbmZ3rkH+HFJoC9uSljtrU4CXgPUAWutjSqlGVk0lhBC2LiUZdn4Lu7+H/BXhrcWQt7TRqWzGruBdfLr7UxwcHJjTYg71C9c3OtIjwqMT6L/8CPsvhPNevWKMeqUCLk4ORscS4hFpGgbVWgc/9P+GUqwTRwghMoC71+G37nB5D9R4D16eAM5ZjE5lE0zaxNxjc5lzbA7lPcozpekUPN09jY71iGPBd+iz1J/wmEQmvVmVN2p6GR1JiH+VlmIt+N5QqFZKuQADgVPWjSWEEDbq/A5zoZYUC+3nQdVORieyGXcT7zJy90h2heyiXcl2fF73c9ycbG9IcaVvMKPWnSSfuyu/9alPJU/bW+NNiNTSUqz1BqYBnkAI8Bcg89WEEPbFlAK7xsOuCZCvnHnYM5/trbhvlHO3zzFoxyCuRV9jZJ2RdCrbyebmpyUkp/DlhkCWH7xCw1J5mf52dTyySdsvYfvSUqyV1Vp3Tv2CUqoBsNc6kYQQwsZEhcLq7nDxb6jWGVpPBBfbe6LRKJsubWL03tFkc87Gj61+pHr+6kZHesSNyHj6LPPnyJU79G5ckmEvlsHJUeaniYwhLcXaDKBGGl4TQojM5+Lf5mHP+Lvw6myo3vnpx9iJZFMy0w5P46eAn6iWrxrfN/me/FnzGx3rEYcuRtB32WFiE5OZ3bkGrSsXMjqSEM/kX4s1pVQ9oD6QTyk1JNVbOQBZzlkIkbmZTLB7Euz8DvKUgnfXQoEKRqeyGRHxEQzfNZyDNw7SsWxHPqn1Cc6OzkbH+getNT/tu8Q3G09RxCMrP/eoQ+kC2Y2OJcQze9KdNRfMa6s5Aan/dd8FZCEhIUTmFR0Gq3vAhR1Q+S1oMwVc3Y1OZTMCbgUwaOcgIuIi+KrBV7xW6jWjIz0iLjGFkWtOsObIVVqUL8DkjlXJ4WZbxaQQafWvxZrWehewSyn1k9b6cjpmEkII41zaC799CLER0Ha6eWkOG5sob6Q159bw9YGvyZMlD/9r/T8q5qlodKRHBEfE0muJP6du3GVoyzL0a1oKBwf5OxQZV1rmrMUqpSYCFYEHz2BrrW2vA68QQvxXJhPsnQLbv4bcxaHHKnOPTwFAUkoS4w6NY+XZldQpVIcJjSbg4WZ7LZn+CrjBx78eR2vNj11r0bSc7c2hE+JZpaVYWwb8ArTBvIxHVyDMmqGEECJdxYTDml4QtAUqdYC208BV5jbddzP2JoN3DuZ42HHer/Q+A6sPxMkhza2l08X1yDjGrA9gc0Ao5QvlYE7nGnjnlSd2ReaQlp+2PFrrH5RSH6UaGt1l7WBCCJEurhwwN2GPCYNXvgefD2XYMxX/UH+G7hxKbHIskxpP4iXvl4yO9A/JKSYW77/M5L/OkKI1n7QqR/cXiuMsy3KITCQtxVrSvT+vK6VeAa4B0pdDCJGxmUywfwZs/RJyFYEPt0Dhakanshlaa5afXs4k30l4Zvdk4YsLKZW7lNGx/uFY8B1GrjlBwLW7NC2bj7GvVqKIR1ajYwlhcWkp1r5WSuUEhmJeXy0HMNiqqYQQwppiI2BtHzi7Ccq3g1dngpu0HLovLjmOr/Z/xYYLG2ji1YRvXviGHC45jI71wN34JL7ffIb/HbhMPndXZneuwcuVCtpcxwQhLOWJxZpSyhEorbX+HYgEmqZLKiGEsJZgX/j1fYi6YW7AXrunDHumEhIVwuCdgzkTcYa+1frSq0ovHJRtDClqrfnjxA2+3BBAWHQC79UtxtCXysqSHCLTe2KxprVOUUq1A6akUx4hhLAOreHAbNgyGnIUhg83g2dNo1PZlH1X9zF893BMJhMzm8+kkVcjoyM9EBwRy+frTrLzTBgVC+dgwXs+VC2Sy+hYQqSLtAyD7lNKzcT8RGjM/Re11oetlkoIISwp7jas6w+nf4dybczDnllyG53KZmit+eHkD0w/PJ2SuUoyrek0iuYoanQsAJJSTCzYfYHp287hqBSft6lA13rFpK+nsCtpKdbq3/tzbKrXNCDrrAkhbN/Vw7CqK9y9Bi99B3X7yLBnKjFJMYzaM4qtV7bSyrsVX9b/kqzOtjFJ3+9SBCPXnOBsaDQvVSzAmHYVKZQzi9GxhEh3Ty3WtNYyT00IkfFoDYfmw+bPIHtBeH8TFKlldCqbciHyAoN2DOLK3SsM8xnGexXes4lJ+ndiExm/6TQ/HwrGM1cWFr7nQ4sKBYyOJYRhbGtVQyGEsIT4SFg/AALXQZlW8NocyGp7q+0baduVbXy25zNcHV2Z33I+tQvVNjoSWmvWHr3K17+f4k5cEj0bleCj5qXJ5iq/qoR9k58AIUTmcu0orOoGd65Ay7FQbwA4yPym+1JMKcw6OosFJxZQMU9FpjadSsFsBY2OxYWwaD5fd5K9QeFUK5KLJe0rU6Gw7SwXIoSRpFgTQmQOWoPfD7DpU8iaF97/A4rWNTqVTYlMiOST3Z+w9+pe2pdqz2d1zXfWjJSQnMLcnReYtTMIVycHvn6tEu/ULiqN14VI5anFmlIqK+YFcYtqrXsopUoDZe+tvSaEEMZLiIL1AyFgNZRqAe3nQ7Y8RqeyKWcizvDRjo8IjQ1ldL3RvFH6DcPnp+07f4tRa09yISyGtlUL83mb8uTP7mZoJiFsUVrGBhYBCUC9e9shwNdpOblSqpVS6oxSKkgpNeIx7w9RSgUqpY4rpbYppYqleq+rUurcvf91Tcv1hBB26MYJmNcYAtdC89Hwziop1B7y+4Xf6fJHF5JSkvip1U+8WeZNQwu18OgEhqw8yjsLDpKcoln8QW1mvF1dCjUh/kVahkFLaq07KqXeBtBax6k0/JTf634wC2iJucDzVUqt11oHptrtCOCjtY5VSvUBJgAdlVIewBeAD+ZlQvzvHXv7mb47IUTmpTUcXgx/fgJuuaDr7+DdwOhUNiXJlMRkv8ksPbWUGvlr8H2T78mbJa9heUwmzSr/YL778zQxCcn0b1qK/s1K4ebsaFgmITKCtBRriUqpLJiLJpRSJTHfaXua2kCQ1vrCveNWAK8CD4o1rfWOVPsfALrc+/olYIvWOuLesVuAVsDPabiuECKzS4iGjUPg+C9Qoim8vgDc8xmdyqbcirvFsF3D8A/1p0v5LgzxGYKzg3Ftmc6GRvHZmhP4XrpNbW8PvmlfidIFshuWR4iMJC3F2hfAJqCIUmoZ0ADolobjPIHgVNshQJ0n7P8h8OcTjvVMwzWFEJldaKB5kdvwIGj6GbwwFBzkzkxqx8KOMWTHEO4m3uXbht/StmRbw7LEJaYwY/s55v99AXc3JyZ0qMIbNb3kAQIhnsHTGrkr4DTwOlAXUMBHWutbaTj3434S9b9cpwvmIc/Gz3KsUqon0BOgaFHbaI0ihLCiI8tg41BwzQ7vrYPittO70lasOruK7w5+R/6s+VnSegnlPMoZlmXnmZt8vu4kwRFxdKjhxcjW5cjjbuzTp0JkRE9r5K6VUmu11jWBjc947hCgSKptL+DawzsppVoAnwGNtdYJqY5t8tCxOx+Tbz4wH8DHx+exhaAQIhNIjIU/hsHRZeD9AnT4AbLLivapJaQk8N3B7/jt3G80KNyA8Y3Gk9M1pyFZbt6N58vfA9l4/Dol8mXj5x51qVdSHvoQ4r9KyzDoAaVULa217zOe2xcorZQqDlwFOgHvpN5BKVUdmAe00lrfTPXWZuBbpdT9TssvAp8+4/WFEJlB2BlY2ZX/Y+++w6Ou8v7/P096JaGEkJAOoQVCCx0UpEVEUEFFFLHciq6ia9l1FWworu12Leu66uqKIIKC0quI0ktCCYQaQnohBNLrzJzfH5O9f6xfhEGZ+UyS9+O6vMwkn+G8rsNk5s35nELRMbj2Get/ctvzvxRUFvDE5ic4XHyYB3o8wCO9HsHVgD4yWzRf7c7krXXHqTVbeHJ0J2ZcG4Onm/x9CfF72FKsjQBmKKUygUqstyi11jr+Uk/SWpuUUo9iLbxcgc+11qlKqTlAktZ6BfAW4Ad827DANEtrPUFrfU4p9QrWgg9gzn8WGwghmpGDi2DVE+DuA9O+gw7XGZ3I6ewt2MvTPz9NrbmWd4e/y8jIkYbkSM0r5bnvD3Mwu4ShHdvwyk3diW7ja0gWIZoapfWl7x5euPfZhbTWmXZJ9BslJCTopKQko2MIIa6G+mpY8yfYPx8ih1hve7YIMTqVU9Fa8+WRL/lb8t+IaBHBuyPeJSYgxuE5KmtN/G3jCf69I4OWPu48P74bE3qGGr7hrhDOTimVrLVOsOVaW0bWZC6YEMJxzp603vY8k2pd6Tn8OXCVk/EuVFVfxUs7XmJtxlpGRozk1SGv4ufh5/AcG1ILeGlFKnmlNUwdEMEzY7sQ4GPc9iBCNFW2vAOuxlqwKcALiAaOA3F2zCWEaG7MJusmtxtfAFcPuHMpxI4yOpXTyS7L5vGfHiftfBqP93mc+7rfhzUoC7IAACAASURBVIty7EH1eSXVvLgilY1HCunSzp8Ppvahb2TLyz9RCPGbXLZY01r3uPCxUqoPMMNuiYQQzU/aJtgwG84csa72vPljCJCtFS+ktWZV+ir+uvuvKKX4aNRHDGnv2BMbTGYLX+zI4J2NJ9Aanr2+C/cNjcbd1bHFohDNzRXfW9Ba71NK9bNHGCFEM3PmmLVIS9sILaPgtvnQ9UaQ+U7/5Wz1WebsnMPm7M30CurFX4f9lTD/MIdmOJBdwnPfHeJIfhnXdWnLyxPiCG/l49AMQjRXly3WlFJPXvDQBegDFNktkRCi6as8Cz/9FZL+DR5+MOZV6P8guMmGqb+0IWMDr+x6har6Kp7q+xTTuk1z6LYcZTX1vL3+OPN3ZdLW35OP7uxDYvd2soBACAeyZWTtwsPbTFjnsC21TxwhRJNmqoXdH8OWt6GuAhLug+HPgq9smPpLJTUlvLb7NdZmrCWudRxzh86lQ2AHh7WvtWb1oXzmrDzC2Ypapg+K4qkxnfD3kgUEQjiaLXPWXnZEECFEE6Y1HFluXTxQkgmxY6yjaUGdjU7mlH7O/pmXdr5ESU0Jj/R6hPt73O/QQ9iziqt4fvlhfj5RRPf2LfjX9ATiwwId1r4Q4r/9arGmlFrJJbbt0FpPsEsiIUTTkpsM62dB1k5o2w3u+g46GrNxq7Mrryvnzb1vsixtGbEtY/nHyH/QtXVXh7VfZ7Lw6dZ03t90EjcXxQvju3H3oEjcZAGBEIa61Mja2w5LIYRoekpzYdPLkLIYfINg/LvQe5rsmfYrdubt5IUdL3Cm6gwP9HiAh3o+hIerh8Pa35txjlnfH+JEYQWJce14cUI3QgK8Hda+EOLX/eq7ptb6Z0cGEUI0EbUVsP092PEBaAsMfRKGPgFeLYxO5pSq6qt4J/kdFh9fTFSLKOZfP5/4oEue5ndVlVTV8dc1x1iclE37QG8+m57AyK7BDmtfCHF5l7oN+o3W+jal1CEucjv0cmeDCiGaGYsFDi6ETa9ARQF0nwQjX4SWFz2xTgDJhcnM3jab3IpcpnWbxmO9H8PLzcshbWut+W5fLnPXHKW0up4Z18Tw+KhYfDxk5FMIZ3Op38rHG/4/3hFBhBCN2OmtsP45KEiB9glw+3wI7290KqdVY6rhg/0fMP/IfEL9Qvl87OcktLPpiMCr4lRRBc8vO8yOU8X0jgjktZt70DVERj6FcFaXug2a3/B/pzqwXQjhRIpPWVd4HlsFAeHWA9e7T5JNbS/hUNEhZm2fxenS09ze+Xae7PskPu6O2Vy2pt7MRz+d4qOfTuHl7sLcm7tzR78IXFzk70sIZ2bLprgDgQ+AroAH4ApUaq3ln2FCNFfV5+HnN2HPJ+DmBSNfgIF/AHeZkP5r6s31fHTwIz4//DltvNvw8aiPGdx+sMPa35F2llnLDnP6bCUTe4Uy+4ZuBPnLJsRCNAa2TE74OzAF+BZIAO4GOtozlBDCSZnrYe9n8PPrUFNqXd153Wzwa2t0Mqd2/Nxxntv2HCfOn2Bih4n8uf+faeHhmH/vnq2oZe7qo3y/P5eo1j7Mv78/w2KDHNK2EOLqsGkmqdY6TSnlqrU2A/9WSu2wcy4hhDPRGk6ss57jWZwG0dfC2NegXXejkzk1k8XE54c/56ODHxHgEcAH133A8PDhDmnbYtEsTsrm9bXHqKoz8dh1HfnDiI54uTvuqCohxNVhS7FWpZTyAA4opd4E8gFf+8YSQjiN/BTYMAtOb4HWsTD1G+sJBDIv7ZLSS9KZtW0Wh4sPkxiVyKwBswj0cswpAMcLypn1/SGSMs8zILoVc2/uTse2/pd/ohDCKdlSrE3DeoD7o8ATQDgwyZ6hhBBOoLwAfnwV9i8A70C4/i1IuBdc5WzISzFbzCw4uoD3972Pj7sPb137FolRiQ5pu7rOzHubTvKvren4e7nx1uR4JvcNk0PXhWjkLrXPWoTWOuuC1aA1gJwTKkRTV18NO/4O2/4G5joY9Ahc8zR4tzQ6mdPLLstm9vbZ7Duzj+Hhw3lx0Iu08W7jkLY3HzvD88sPk3O+mlv7hvHsuK608nXcCQhCCPu51MjaMqAPgFJqqdZaRtOEaMosFji8BH54GcpyoOuNMOplaN3B6GROz6ItfHP8G95Jfgc35cbcoXO5MeZGh4xoFZbVMGflEVYfyqdjWz8WPziQATGt7d6uEMJxLlWsXfguE2PvIEIIA2Xtsm5qm5sMIT3hlo8haqjRqRqF/Ip8nt/xPLvzdzM4dDAvD36Zdr7t7N6u2aJZsCuTt9cfp9Zs4ekxnXjwmg54uMmh60I0NZcq1vSvfC2EaCrOZ8APL0Hq9+AfAjf9E+JvBxf5wL8crTXL0pbx5t43MWszzw98nls73eqQ0bTDuaXM+v4QB3NKGRbbhlcmdieqjaz7EqKpulSx1lMpVYZ1hM274WsaHmvZFFeIRqymFLb+L+z6CFzcYPizMHgmeMgHvi2Kqop4aedLbMnZQkJwAnOGzCHcP9zu7daazLy17jifbz9NK19P3pvSiwk9Q2UBgRBN3KWOm5LNeIRoaswm2DcPNr8GVWeh51QY+Ty0CDU6WaOgtWbt6bXM3T2XWnMtz/R7hqldp+Ki7D8SmV9azUML9nEwu4SpAyJ4JrELAd6yMleI5sCmTXGFEE1A2g+wfjYUHYXIITB2CYT2NjpVo3Gu5hyv7nqVjZkbiQ+K59UhrxIdEO2QtnelF/PIV/uoqTfzz7v6kNg9xCHtCiGcgxRrQjR1Z45ZN7VN+wFaRsPtC6DLeNnU9gpsytrEnJ1zKK8r5/E+j3NP3D24udj/7VNrzefbM3htzVEiW/uweNpA2dxWiGZIijUhmqrKs9bbnclfgIcfjJkL/R8ANzm821altaW8vud1VqWvomurrnw65lM6tezkkLar6kw8+90hlh/IY0y3YP73tp74e8ltTyGaIynWhGhqTLWw+5+w5W2oq4R+98O1fwFf2XvrSmzL3caL21+kuKaYh3s+zAPxD+Du4phiKbO4khnzkzleWM6fxnbm4Ws74OIiI6FCNFdSrAnRVGgNR5bBxhehJBNix8KYVyHIMSNBTUVlfSVv7X2LpSeX0jGwI++PfJ+41nEOa3/z8TM8/vV+lFJ8cW9/ru0U5LC2hRDOSYo1IZqCnGTrprbZu6BtHEz7HjpcZ3SqRmdP/h5e2PECeRV53Nv9Xh7p9Qiero65bWyxaD7cnMY7P5ygS7sWfHxXXyJa+zikbSGEc5NiTYjGrDTHejzUoW/ANwhufA96TwMX2XnnSlSbqnlv33t8dfQrIvwj+PL6L+nVtpfD2i+rqeepbw6y8UghE3uF8vot8Xh7yN+hEMJKijUhGqPaCtj+Huz4ALQFhj0FQ58AT1kpeKUOnDnA7O2zySzL5M6ud/J4n8fxdvN2WPsnC8uZMT+ZzHNVvHhjN+4ZHCWb3Aoh/osUa0I0JhYzHFgIP74CFYXQfTKMehECI4xO1ujUmmv58MCHzEudRzufdnw25jP6h/R3aIa1h/J5+tuDeHu4svB/BsgB7EKIi5JiTYjG4vQW67y0gkMQ1h9u/wrC+xmdqlFKLU5l9rbZpJWkMSl2Ek8nPI2fh5/D2jdbNG+tP84/fz5Fr/BAPrqrDyEBjhvNE0I0LlKsCeHszqbBxufh+BoIiIDJn0PcLbKp7W9Qb6nn05RP+STlE1p7teYfI//BsLBhDs1wvrKOxxbtZ+vJs0wdEMGLN3bD003mpwkhfp0Ua0I4q6pz8PObsPdTcPOGkS/CwIfBXUZgfouT508ya9ssjp47yviY8fyl/18I8AxwaIbDuaXMmJ9MUXktb0zqwe395Pa1EOLy7FqsKaUSgfcAV+BfWuvXf/Hza4B3gXhgitZ6yQU/MwOHGh5maa0n2DOrEE7DVAdJn8FPr0NtGfS5G0bMAr+2RidrlMwWM1+kfsGHBz7E38Ofd4e/y8jIkQ7PsTQ5h+e+P0QrXw++eWgQvcIDHZ5BCNE42a1YU0q5Ah8Co4EcYK9SaoXW+sgFl2UB9wBPX+SPqNZaO27tvBBG0xqOr4UNs+HcKYgZAWPnQrDjNmRtajJKM5i1fRYpRSmMjhzN7IGzaeXVyqEZ6kwW5q4+wrydmQyMacXfp/ahjZ8c+SWEsJ09R9b6A2la63QApdQiYCLwf8Wa1jqj4WcWO+YQwvnlp1gXD2RshTadYOq3EDta5qX9RhZtYeHRhby37z08XD14Y9gbXB99vcO3xDhTVsMjC/exN+M8DwyL5pnELri5ujg0gxCi8bNnsdYeyL7gcQ4w4Aqe76WUSgJMwOta62VXM5wQTqG8wLoNx/6vwLsljHsb+t4DrnJg92+VU57DCzteYG/BXq4Ju4YXB71IWx/H30JOzjzHwwv2UV5j4v07ejOhZ6jDMwghmgZ7FmsX+yesvoLnR2it85RSMcCPSqlDWutT/9WAUg8CDwJERMhEXdGI1FXBzg9h29/AXAeDHoFr/gTeMo/pt9Jas+TkEt7e+zZKKeYMnsNNHW9y+Gia1poFu7OYszKV0EBv5t3Xn64hLRyaQQjRtNizWMsBwi94HAbk2fpkrXVew//TlVI/Ab2BU7+45hPgE4CEhIQrKQSFMIbFAoeXwA8vQVkudJ0Ao1+GVjFGJ2vUCioLeGnHS2zP286AkAHMGTyHUD/Hj2TV1JuZvewwS5JzGNE5iHdv702Aj4ySCiF+H3sWa3uBWKVUNJALTAGm2vJEpVRLoEprXauUagMMAd60W1IhHKEwFVY8BrlJENILbvkUooYYnapR01qzMn0lr+9+HZM2MWvALG7rfBsuyvHzwnLOV/Hwgn0cyi3lsZGx/HFkLC4uMudQCPH72a1Y01qblFKPAuuxbt3xudY6VSk1B0jSWq9QSvUDvgdaAjcqpV7WWscBXYGPGxYeuGCds3bkV5oSwrmZ62HrO7DlLfAKgJv+CfG3g4tMNP89zlafZc7OOWzO3kyftn14ZcgrRLQwZjrE9rSzzPx6P/UmC/+6O4FR3YINySGEaJqU1k3j7mFCQoJOSkoyOoYQ/y3vACx/FAoPWc/xvP5N8JXzH3+v9RnreXXXq1TVV/FYn8e4q+tduLo4/hQArTWfbEnnjXXH6BDkx8fT+hIT5Lhjq4QQjZdSKllrnWDLtXKCgRD2YKq1nj6w7W/gGwRTFkKXG4xO1eiV1JTw2u7XWJuxlu6tuzN36FxiAo2Z71dZa+LPS1JYfSifcT3a8dbknvh6yluqEOLqk3cWIa62nCRY/ggUHYNed1o3tvVuaXSqRu+n7J94eefLlNSWMLP3TO7rfh9uLsa8hZ0+W8mM+Umknang2eu78OA1MQ5fdSqEaD6kWBPiaqmvhs1zrVty+IfAnUshdpTRqRq98rpy3tjzBstPLadTy058NOojurTqYlieH44U8sTiA7i5Kr68bwBDY9sYlkUI0TxIsSbE1ZC50zqadu4U9L0XRs8BL9lb6/fakbeDF3e8SFFVEQ/0eICHez6Mu0EbBlssmnc3neT9TSfp3r4F/7yrL2EtfQzJIoRoXqRYE+L3qK2ATXNgzycQGA53L4eY4UanavSq6qt4J/kdFh9fTHRANPOvn0+PoB6G5SmtrueJxQf48dgZJvUJY+7N3fFyd/yCBiFE8yTFmhC/VfrPsGImlGRC/xkw8gXwlJWAv1dyYTKzt80mtyKXu7vdzczeM/Fy8zIsz7GCMmbMTyb3fDWvTIzjroGRMj9NCOFQUqwJcaVqymDj85D8BbTqAPeuhcjBRqdq1Krqq9iSs4V1Gev4MetH2vu159+J/6ZvcF9Dc608mMefl6Tg7+XG4hkD6RvZytA8QojmSYo1Ia7EyR9g5WNQng+DZ8Lw58BD5i39FjWmGrbmbmXd6XVsydlCjbmGIO8g7u1+LzPiZ+Djbly/mswW3lh3jE+3niYhsiX/uLMPbVsYN7onhGjepFgTwhbV52H9LDjwFbTpDPdvhDCb9jIUF6gz17EtdxvrMtbxU/ZPVJuqaeXViokdJzI2aix92vYxZHPbCxVX1PLowv3sTC9m+qBIZt3QDQ83OW1CCGEcKdaEuJxja2DVE1BZBMOegmufATdPo1M1GvXmenbm72R9xnp+zPqRivoKAj0DGRc9jsToRBKCEwzbL+2XDmaX8PCCZIor63j71p5M7htmdCQhhJBiTYhfVVkMa/8Mh5dAcHeYuhhCexmdqlEwWUzsyd/Duox1bMraRFldGf4e/oyKHEViVCL9Q/rj7mLMFhy/5pu92cxefpggP0+WPjyY7u0DjI4khBCAFGtCXFzq97D6aagptc5LG/oEuHkYncqpmS1mkgqTrAVa5ibO157H192X68KvY2zUWAaHDjZsj7RLqTWZeXnlERbuzmJoxza8f0dvWvnK37UQwnlIsSbEhSrOwOqn4OgKCO0NE1dAcJzRqZyWRVvYV7iP9Rnr2Zi5keKaYrzdvBkeNpyx0WMZ2n4onq7Oe8u4oLSGh79KZn9WCQ9d24E/je2Mq4tsyyGEcC5SrAkBoDWkfAPrnoG6Khj1EgyaCa7yK/JLWmsOFh1kfcZ6NmRs4Ez1GbxcvRgWNozEqESGhQ3D283b6JiXtef0Of7w1T6q6kx8dGcfru8RYnQkIYS4KPkkEqIsz7qA4MQ6COsPEz+EoE5Gp3IqWmtSi1NZd3odGzI3kF+Zj4eLB0PbD2Vs1FiGhw83dKuNK6G15osdGcxdfZSIVj58/cAAYoP9jY4lhBC/Soo10XxpDfsXWLfkMNfB2L/CgBlg8NYRzkJrzbFzx1iXsY71GevJrcjFzcWNwaGDmdl7JiPCR+Dn0bhObKiuM/Pc94f4fn8uo7oG887tPWnh5Xzz6IQQ4kJSrInmqSQLVjwG6ZshcihMeB9adzA6lVM4ef7k/xVomWWZuCpXBoYMZEb8DK6LuI4Az8a5SjL7XBUz5idztKCMp0Z34pERHXGR+WlCiEZAijXRvFgskPw5bHzROrI27m1IuB9cmvemp+ml6aw/vZ71Ges5VXoKF+VCv3b9mB43nVERo2jp1dLoiL/LzyeKeOzr/Wit+Xx6P0Z0aWt0JCGEsJkUa6L5OJduHU3L2AoxI+DG96BlpNGpDJNVlsX6jPWsy1jHifMnUCj6BPdhVpdZjIocRRvvNkZH/N201vzjp1O8veE4nYP9+XhaXyJb+xodSwghrogUa6Lps5hh98ewaQ64usOED6D3NFDN7xZYbkWutUA7vY6j544C0CuoF3/p/xdGR46mrU/TGXEqr6nn6W8Psj61kAk9Q3l9Ug98POQtTwjR+Mg7l2jaik7AikchezfEjoHx70JAe6NTOVRBZQHrM6y3OA+dPQRAjzY9eDrhacZEjiHEr+ltWZF2poIZ85PIKK5i9g1duX9oNKoZFudCiKZBijXRNJlNsPMD2PxXcPeGmz+G+NubzWhaUVURGzI3sD5jPfvP7Aega6uu/LHPHxkbNZYw/6Z75uX61AKe+uYgnm4uLLh/AIM6tDY6khBC/C5SrImmp/AILP8D5O2HLuPhhnfAP9joVHZXXF3MD5k/sC5jHcmFyWg0sS1jmdl7JmOjxhLZomnPzzNbNO9sPM6Hm0/RMzyQf97Vh5AA59+cVwghLkeKNdF0mOth29/g5zfBqwVM/jfE3dykR9NKakr4IctaoO0t2ItFW4gOiOahng+RGJVITGCM0REdoqSqjscWHWDLiSKm9Avn5YlxeLrJfnlCiKZBijXRNOQfhGWPQOEh6D4Jrn8TfBv/asaLKasrY1PmJtZnrmd33m5M2kSEfwT3d7+fxOhEYgNjm9X8rCN5ZcxYkERhaS1/vaUHd/SPMDqSEEJcVVKsicbNVGsdSdv2N2txdvtX0HW80amuuoq6CjZnb2Z9xnq2523HZDHR3q89d8fdTWJUIl1adWlWBdp/LNufy1++SyHQ24PFMwbSO6Jx7wcnhBAXI8WaaLxykq1z04qOQc+pMHYu+LQyOtVVU1Vfxc85P7Pu9Dq25W6jzlJHsE8wU7tMJTEqke5tujfLAg2g3mzhtTVH+ff2DPpHt+LDqX0I8vc0OpYQQtiFFGui8amvhs2vwc6/g38I3LkEYkcbneqqqDZVszVnK+sz1rMlZws15hqCvIO4tfOtJEYlEh8Uj4tq3qctFJXX8sjCfew5fY77hkTz7LguuLs27z4RQjRtUqyJxiVzJyx/BM6dgr73wOg54NU4z6r8j1pzLdtzt7MuYx0/Zf9EtamaVl6tmNhxIolRifRu2xtXOVwegH1Z53l4QTKl1fW8N6UXE3s1rz3zhBDNkxRronGoq4QfXoY9n0BgONy9HGKGG53qN6s317MzfyfrTq9jc/ZmKuorCPQMZFz0OBKjE0kITsDNRX49/0NrzcI9Wby0IpV2AV589/AQuoW2MDqWEEI4hHwaCOeX/jOsmAklmdB/Box8ATz9jE51xSzawq78Xaw9vZZNWZsoryvH38OfUZGjSIxKpH9If9xd3I2O6XRq6s28uDyVxUnZXNspiPem9CLQx8PoWEII4TBSrAnnVVMGG1+A5H9Dqxi4dy1EDjY61RUrrCxkWdoyvjv5HXmVefi6+3Jd+HWMjRrL4NDBuLtKgfZr8kqqeXhBMgdzSpl5XUf+OKoTri7Nc1GFEKL5kmJNOKeTP8DKx6E8DwY9CiNmgYeP0alsZrKY2Ja7jaUnlrIldwsWbWFAyACe6PsEIyJG4OkqKxcvZ8eps8xcuJ9ak4WPp/VlbFw7oyMJIYQhpFgTzqX6PKyfBQe+gjad4b4NEN7P6FQ2y6vI47uT3/F92vecqTpDa6/W3Bt3L5NiJxHeItzoeI2C1prPtp3mr2uPEd3Gl4+n9aVDUOO77S2EEFeLFGvCeRxbA6uegMoiGPYUXPNncPcyOtVl1Vvq+Tn7Z5acWMKOvB0ADGk/hOf6P8c14dfIPLQrUF5Tz3PfH2blwTyu796Ot27tiZ+nvE0JIZo3eRcUxqsshnXPwKFvIbg7TF0Eob2NTnVZWWVZLD25lOVpyymuKSbYJ5gZPWdwc8ebCfULNTpeo6G1JiWnlEV7s1hxII/qejPPJHbhoWtjmu2mv0IIcSG7FmtKqUTgPcAV+JfW+vVf/Pwa4F0gHpiitV5ywc+mA7MbHr6qtZ5nz6zCIKnLYM3T1tufw5+FoU+Cm/Ou9Ksz17EpaxNLTixhT8EeXJUrw8KGcWunWxkSOkT2Q7sCZTX1LN+fy9d7sjmSX4a3uyvj40OYNiiS+LBAo+MJIYTTsFuxppRyBT4ERgM5wF6l1Aqt9ZELLssC7gGe/sVzWwEvAgmABpIbnnveXnmFg1WcgdVPwdEVENILpi2Ddt2NTvWr0kvSWXJyCStPraSktoT2fu2Z2XsmN3W8ibY+bY2O12horUnOPM/Xe7JZfSiPmnoLcaEteOWm7kzsFUoLL7llLIQQv2TPkbX+QJrWOh1AKbUImAj8X7Gmtc5o+JnlF88dC2zUWp9r+PlGIBH42o55hSNobb3dufbP1o1uR74Igx8DV+e7I19tqmZj5kaWnljKvjP7cFNujIgYweROkxkYMrDZH/t0Jc5X1vHd/lwW7cni5JkKfD1cuaVPGHf0i6BHWOM+gUIIIezNnp+Q7YHsCx7nAAN+x3PlXJnGriwPVj0JJ9ZCWD+Y+CEEdTY61f/j+LnjLDmxhNXpqymvLyeyRSRP9n2SCR0m0Nq7tdHxGg2tNbvSz7FobxZrDxdQZ7LQKzyQNyb1YHx8KL6ycEAIIWxiz3fLi80M1lfzuUqpB4EHASIiImxPJhxLa+tWHOueA3MtjJkLAx8GJ5rfVVVfxdrTa1l6cimHzh7Cw8WDUZGjmNxpMgnBCTLR/QqcrahlSXIOi/dmc/psJS283LijXzhT+kfQNUSOiBJCiCtlz2ItB7hwY6kwIO8Knjv8F8/96ZcXaa0/AT4BSEhIsLUQFI5UkmXd3PbUjxA5BCZ8AK07GJ0KsI78HCk+wrcnvmXt6bVUmaroENCBZ/o9w/iY8QR6ySR3W1ksmm1pZ1m0N4sNqYWYLJr+Ua2YeV1HxvUIwcvdeQpzIYRobOxZrO0FYpVS0UAuMAWYauNz1wOvKaVaNjweAzx79SMKu7FYIPlz2PiidWRt3NuQcD+4GD/Pq7yunNXpq1l6cinHzh3Dy9WLsVFjmdxpMj2Desoo2hUoKK3h26RsFidlk3O+mpY+7twzOIop/cPp2Nbf6HhCCNEk2K1Y01qblFKPYi28XIHPtdapSqk5QJLWeoVSqh/wPdASuFEp9bLWOk5rfU4p9QrWgg9gzn8WG4hG4Fw6rHgMMrZCzHC48X1oGWloJK01B4sOsuTEEtZnrKfGXEOXVl2YNWAW42LG0cJDbs/ZymS28POJIr7ek82PxwqxaBjcoTXPJHZhTFwwnm4yiiaEEFeT0rpp3D1MSEjQSUlJRsdo3ixm2PMJbJoDLm4w5lXoczcYOFJVUlPCyvSVLD2xlFOlp/Bx82FczDgmx06mW+tuMop2BXLOV/FNUg7f7M2moKyGNn6e3JoQxu0J4US18TU6nhBCNCpKqWStdYIt18pyLPH71VfDgYWw80M4dwpix8D4dyHAmAW8WmuSCpNYcmIJP2T+QJ2ljh5tevDSoJe4Pvp6fNwbz4HwRqs3W9h0tJCv92Sz5WQRANfEBvHShDhGdm2Lu6vxt7WFEKKpk2JN/HaVxbD3U+toWlWx9Yio276ErhMMGU0rri5mxakVLD25lMyyTPzd/bkl9hYmd5pM51bOt0WIM8ssrmTR3my+TcrhbEUt7Vp4MfO6WG5LCCOspRS7QgjhSFKsiStXfAp2/QP2fwWmauiUCINnWld7OrhIs2gLu/J2seTkEjZnb8ZkMdGnbR8ejH+Q0ZGj8XbzdmiexqzWZGZ9aiGL9mSx41Qxri6KEZ3bckf/cK7tFISbjKIJIYQhpFgTvrZ1GQAAH8JJREFUtstJgu3vwdGV4OoO8bfBoJnQtovDo5ypOsOytGV8d/I7cityCfQM5I4udzApdhIdAp1ja5DGIu1MBYv2ZLF0Xw7nq+oJa+nN02M6cWtCOMEtvIyOJ4QQzZ4Ua+LSLBY4sQ52vA9ZO8ErAIY+AQNmgH87h0YxW8xsy93GkpNL2JqzFbM2079dfx7v8zgjI0bi4eq8B8A7m5p6M2sO5fP1niz2ZpzHzUUxJi6YKf0iGNqxDS4usvBCCCGchRRr4uLqayBlEez4OxSfhIAISHwdek8DTz+HRsmvyOe7tO/4/uT3FFYV0sqrFdPjpjMpdhIRLeTkiitxNL+MRXuy+H5/LmU1JqLb+PKX67swqU8YQf6eRscTQghxEVKsif9WdQ72fgZ7PobKIgjpCZM+g243OfSw9XpLPT9n/8ySk0vYkbsDgMGhg3mm/zMMDxuOu6u7w7I0dpW1Jlal5LFwTzYHs0vwcHXh+h7tmNIvgoExrWT7EiGEcHJSrAmr8xnWrTf2L4D6Kug4GoY8BlHDHLpoILssm6Unl7IsbRnFNcW09WnLg/EPcnPszbT3M2YrkMbqUE4pC/dkseJALpV1ZmLb+vH8+G7c0rs9LX3llrEQQjQWUqw1d7nJsOMDOLIclGvDooFHIbibwyLUmev4MetHlpxcwu783bgoF65pfw2TOk1iaPuhuLnIy9RWZTX1LD+Qx6I9WaTmleHl7sINPUKZOiCcPhEtZRRNCCEaIfkUbI4sFji5wVqkZW4DzwAY/Jh10UCLUIfFSC9NZ+mJpaw8tZLztecJ9Q3l0V6PclPHmwj2DXZYjsZOa82+rBK+3pPF6pR8quvNdAtpwSsT45jQqz0B3nLLWAghGjMp1poTUy2kLLYuGjh7HFqEwZi51iOhvBxzNmaNqYaNmRtZcmIJ+87sw025MSJiBJNiJzEwZCCuLnKupK1Kqur4bl8ui/ZmcaKwAl8PV27qHcod/SPo0T5ARtGEEKKJkGKtOag+37Bo4BOoKITgHnDLpxB3s3W/NAc4cf6EdRQtfSXldeVE+Efwxz5/ZGLHibTxbuOQDE2B1prdp8/x9Z4s1h4uoM5koWdYAK/f0oPxPUPx85RfaSGEaGrknb0pO58Juz6CfV9CfSV0uA5u/hhihjtk0UBVfRXrMtax9MRSUs6m4O7izqjIUUyOnUxCuwRclOyIb6uzFbUsTc5h8d5s0s9W4u/lxpR+4UzpF0G3UMeMigohhDCGFGtNUd4B6ya2qcusRVn3ydbjoNp1d0jzqcWpLD2xlDWn11BZX0lMQAx/SvgTN3a4kZZeLR2SoSmwWDTbT51l0Z5sNhwpoN6sSYhsySMjOjKuRwjeHnLLWAghmgMp1poKrSHtB+txUBlbwcMfBv0BBjwEAWF2b77OXMfq9NV8fexrjp47iqerJ2OjxjK502R6BfWS+VNXoLCshm+TslmclE32uWoCfdy5e1AUU/qFExvsb3Q8IYQQDibFWmNnqoND31pXdhYdBf9QGP0K9J1uPRrKzkprS/nm+DcsPLaQs9Vn6RjYkecGPMcNMTfQwkNuz9nKbNH8fOIMX+/J5sdjZzBbNINiWvP0mM6MjWuHl7uMogkhRHMlxVpjVV0Cyf+G3R9DeT60jYOb/gndJ4Gb/Tc8zS7PZv6R+SxLW0a1qZohoUOYO3Qug0IGySjaFcgtqWbx3my+Tcomv7SGNn4ePDAshtv7hRPdxtfoeEIIIZyAFGuNTUl2w6KBeVBXYV0sMPHv0GGkQxYNpBSl8EXqF2zK2oSLcmFc9Dju7nY3nVt1tnvbTUW92cKmo2dYtDeLn08UATAsNogXxndjZNdgPNxk4YUQQoj/nxRrjUV+ivVW5+Gl1sfdb7EuGgjpafemzRYzP+X8xLzUeew/sx9/d3/ujbuXqV2n0tanrd3bbyoKy2r4cmcG3yTlUFReS3ALTx4d0ZHbEsIJb+VjdDwhhBBOSoo1Z6Y1nPrRurIz/Sfw8LMuGBj4MASG2735alM1K9JWMP/ofDLLMmnv155n+j3DzbE34+sut+hsdbygnE+3prP8QC5mi2ZE57bc0T+C4Z2DcHOVUTQhhBCXJsWaMzLVQep31pG0wsPg1w5GvQR97wXvQLs3X1xdzKLji1h0bBEltSV0b92dt659i1ERo+ScThtprdmeVswnW9PZcqIIb3dXpvaP4L6h0US2lkJXCCGE7eST15nUlEHyF9Y5aeV5ENQVJv4DekwGN0+7N59ems6XqV+y8tRK6ix1DA8fzvRu0+kb3FcWDdiozmRhVUoen2xJ51hBOW38PPnT2M7cOSCCQB/7L/wQQgjR9Eix5gxKc2H3R5A8D2rLIGoY3PgexI62+6IBrTVJhUl8mfolP+X8hIeLBxM7TmRat2lEB0Tbte2mpLS6nq/3ZPHF9gwKymqIbevHm5Pimdg7FE832XZDCCHEbyfFmpEKDjcsGlhinZ8WdxMMehTa97F70yaLiR8yf+CL1C9ILU6lpWdLHu75MLd3vp3W3q3t3n5TkXO+in9vz2DRniwq68wM7tCav07qwbWxQbi4yGikEEKI30+KNUfT2rpYYMf71sUD7r7Q7wHrooGWkXZvvrK+ku9OfseCIwvIq8wjskUkzw98ngkdJuDl5mX39puKQzmlfLI1nTWH8gG4MT6E/xkWQ/f29t+IWAghRPMixZqjmOsh9XtrkVZwCPyCYeQL1kUDPq3s3nxhZSELjy3k2+PfUl5fTp+2fXim/zMMDx8uB6rbyGLRbD5+hk+3prMr/Rx+nm7cNySKe4ZE0z7Q2+h4Qgghmigp1uytttw6F23XR1CWA206wYQPIP52hywaOH7uOF8e+ZI1p9dg0RZGRYxietx04oPi7d52U1FTb2bZ/lw+3ZrOqaJKQgK8mDWuK7f3D6eFl7vR8YQQQjRxUqzZS1m+ddFA0hdQWwqRQ+GG/4XYMeBi35EsrTU783byReoX7MzfibebN7d3vp27ut5FmL/9D3VvKs5V1rFgVyZf7szgbEUd3UJa8O7tvbghPgR32R9NCCGEg0ixdrUVHoGdf4eUb0CboesEGPwYhPW1e9P15nrWnF7DvCPzOHn+JEHeQTze53Fu7XQrAZ4yl8pWGWcr+Wzbab5Nzqam3sLwzkE8OCyGQR1ayxYmQgghHE6KtatBa8jYCtvfh7SN4O4DCffCwD9AK/tvf1FWV8a3x79l4dGFnKk+Q8fAjrwy5BXGRY/Dw1X29rJVcuY5PtmSzoYjhbi7uHBT71D+Z1gMnYL9jY4mhBCiGZNi7fcwm+DIMuuigfyD4BsEI2ZDv/sdsmggtyKXBUcW8N3J76gyVTEwZCBzhsxhcOhgGQGykdmi2ZBawCdb09mfVUKAtzuPDO/I3YMjaesvq2OFEEIYT4q136K2AvbPh53/gNIsaN3Ruolt/BRwt/8H/OGzh5mXOo8NmRtwwYXro6/n7ri76dKqi93bbiqq6kwsSc7hX1tPk3WuiohWPrw8IY5bE8Lw8ZBfCyGEEM5DPpWuRHkB7P4Ykj6DmlKIGATXvw6drrf7ogGLtrAlZwtfpH5BcmEyfu5+TO82naldp9LOt51d225KzpTX8OWOTBbszqSkqp7eEYE8e30XxsS1w1U2sRVCCOGEpFizVW0FfJAAdRXQ9UbrooHwfnZvtsZUw8r0lXyZ+iUZZRmE+Ibwp4Q/cUvsLfh5+Nm9/abiZGE5n25NZ9n+POotFsZ0C+bBa2LoG2n/29VCCCHE7yHFmq08/axbb4QlQOsOdm/ufM15Fh1fxKJjizhXc46urbryxrA3GB01GncX2dvLFlprdp4q5tOt6Ww+XoSnmwu39Qvj/qExRLfxNTqeEEIIYRMp1q5Ez9vt3kRGaQbzj8xn+anl1JpruSbsGu6Ju4eE4ARZNGCjerOFNYfy+WRLOql5ZbT29eDJ0Z24a2AkrXxldawQQojGxa7FmlIqEXgPcAX+pbV+/Rc/9wS+BPoCxcDtWusMpVQUcBQ43nDpLq31Q/bMaiStNfvP7Gde6jw2Z2/GzcWNCR0mcHe3u4kJjDE6XqNRXlPPoj3Z/Hv7afJKa+gQ5Mvrt/Tgpt7t8XJ3NTqeEEII8ZvYrVhTSrkCHwKjgRxgr1Jqhdb6yAWX3Q+c11p3VEpNAd4A/jN8dUpr3cte+ZyB2WJmU9Ym5qXOI+VsCgGeATwQ/wB3dLmDNt5tjI7XaOSVVPPFjgy+3p1Fea2JAdGteOWm7ozo3BYXWTQghBCikbPnyFp/IE1rnQ6glFoETAQuLNYmAi81fL0E+LtqBvf6quqr+D7te+YfmU9uRS7h/uHMGjCLCR0m4OPuY3S8RuNwbin/2prOqpR8NDCuRwgPDIsmPizQ6GhCCCHEVWPPYq09kH3B4xxgwK9do7U2KaVKgdYNP4tWSu0HyoDZWuutdszqEEVVRSw8tpBvjn9DWV0ZvYJ68aeEPzE8fDiuLnKbzhZaa346UcSnW9LZcaoYXw9X7h4Uxb1DoghvJYWuEEKIpseexdrFRsi0jdfkAxFa62KlVF9gmVIqTmtd9l9PVupB4EGAiIiIqxDZPtLOpzHvyDxWp6/GZDExMmIk0+Om06ttk77Le1XVmswsP5DHv7amc6KwguAWnvzl+i7c0T+CAG9ZHSuEEKLpsmexlgOEX/A4DMj7lWtylFJuQABwTmutgVoArXWyUuoU0AlIuvDJWutPgE8AEhISflkIGkprze6C3XyR+gXbc7fj5erFpNhJTOs2jYgWzltYOpuSqjq+2p3FFzsyKCqvpUs7f965rSfj40PxcLPvRsRCCCGEM7BnsbYXiFVKRQO5wBRg6i+uWQFMB3YCk4EftdZaKRWEtWgzK6VigFgg3Y5Zr5p6Sz3rTq/jyyNfcuzcMVp7tWZm75nc1uk2Ar1kLpWtsoqr+GxbOt8k5VBdb+aaTkG8c1s0Qzu2kS1MhBBCNCt2K9Ya5qA9CqzHunXH51rrVKXUHCBJa70C+AyYr5RKA85hLegArgHmKKVMgBl4SGt9zl5Zr4byunKWnljKgqMLKKwqJCYghpcHv8wNMTfg6eppdLxGY3/WeT7dms66wwW4uigm9GzP/wyLpmtIC6OjCSGEEIZQ1juOjV9CQoJOSkq6/IVXWX5FPguOLmDpyaVU1lfSv11/psdNZ2j7obgouU1nC7NF88PRQj7dkk5S5nlaeLlx58BI7hkcRXALL6PjCSGEEFedUipZa51gy7VygsFvdKT4CPNS57E+Yz0AY6PGMj1uOt1adzM4WeNRXWdmyb4cPt92mtNnKwlr6c2LN3bjtoRwfD3lpSmEEEKAFGtXxKItbMvdxrzUeewp2IOvuy93dr2Tu7reRYhfiNHxGo2zFbV8uTOT+TszOF9VT8+wAP4+tTeJce1wc5XRSCGEEOJCUqzZqKq+ijvX3ElaSRptfdryVN+nmNRpEv4e/kZHazTSzlTw2bZ0lu7Lpc5kYVTXYB4YFk3/6FayaEAIIYT4FVKs2cjH3YeBIQO5v8f9jI0ai7uL7O1lC601u0+f419b0/nh6Bk83FyY3DeM+4dG0yHIz+h4QgghhNOTYu0KPNP/GaMjNBoms4W1hwv4dGs6KTmltPL14PGRsUwbFEkbP1kdK4QQQthKijVxVVXUmli8N5vPt50mt6Sa6Da+zL25O5P6hOHlLkdqCSGEEFdKijXxu5ktmj2nz7EqJY8VB/MorzHRL6olL97YjVFdg3FxkfloQgghxG8lxZr4TSwWTXLWeVYdzGPN4QKKymvxdndldLdg7h0SRe+IlkZHFEIIIZoEKdaEzbTW7M8uYdXBfNYcyqegrAZPNxeu69KW8fGhjOgShI+HvKSEEEKIq0k+WcUlaa05nFvGqpQ8VqXkk1tSjYerC9d2DuLZ+C6M7BqMn2xgK4QQQtiNfMqK/4fWmqP55axKyWP1oXwyi6twc1EMi23Dk6M7MToumBZesnWJEEII4QhSrIn/c6KwnFUH81h1KJ/0okpcXRSDO7TmkeEdGRMXTKCPh9ERhRBCiGZHirVm7lRRBatT8lmVkseJwgpcFAyMac39Q6NJjGtHa9kTTQghhDCUFGvNUFZxFSsb5qAdzS9DKegX2Yo5E+NI7N6Otv5eRkcUQgghRAMp1pqJnPNVrE7JZ/WhfFJySgHoExHIC+O7Ma5HCO0CpEATQgghnJEUa01Yfmk1aw4VsColj/1ZJQD0DAvguXFduCE+lPaB3gYnFEIIIcTlSLHWxJwpr2FtQ4G2N+M8AN1CWvDnxM6M7xFKRGsfgxMKIYQQ4kpIsdYEFFfUsvZwAatT8tl9uhiLhs7B/jw1uhM3xIcQE+RndEQhhBBC/EZSrDVSJVV1rE8tYFVKPjtOFWO2aDoE+TLzuljGx4cQG+xvdEQhhBBCXAVSrDUiZTX1bEgtZFVKHttOnsVk0US29uGha2MYHx9Kl3b+KCWHpgshhBBNiRRrTq6i1sQPR6wF2pYTZ6kzW2gf6M39w6K5MT6UuNAWUqAJIYQQTZgUa06oqs7Ej8fOsOpgPpuPn6HWZCEkwIu7B0VyQ3wIvcIDpUATQgghmgkp1pxETb2Zn46fYWVKPj8ePUN1vZkgf0/u6B/B+PgQ+kS0xMVFCjQhhBCiuZFizUC1JjNbTpxldUoeG48UUllnprWvB5P6tmd8fCj9olrhKgWaEEII0axJseZgdSYL20+dZdXBfDYcKaC8xkSgjzs39gxlfHwoA2Na4ebqYnRMIYQQQjgJKdYcwGS2sDO9mFUH81mXWkBpdT3+Xm6MjWvH+PgQhnRsg7sUaEIIIYS4CCnW7MRs0ew+XczqlHzWHS6guLIOP083RncLZnx8CENj2+Dp5mp0TCGEEEI4OSnWriKLRZOcdZ5VB/NYc7iAovJavN1dGdUtmBt6hDC8cxBe7lKgCSGEEMJ2Uqz9Tlpr9meXsOpgPmsO5VNQVoOnmwvXdWnL+PhQruvSFm8PKdCEEEII8dtIsfYbaK05lFvKqpR8Vqfkk1tSjYerC9d2DuLZ+C6M7BqMn6d0rRBCCCF+P6kobKS15mh+OatS8lh9KJ/M4ircXRXDYoN4akwnRnULpoWXu9ExhRBCCNHESLFmo1qThVv/uYMak4XBHVrzyPCOjI1rR4CPFGhCCCGEsB8p1mzk5e7KJ3cn0KWdP639PI2OI4QQQohmQoq1KzCkYxujIwghhBCimZGdWIUQQgghnJgUa0IIIYQQTsyuxZpSKlEpdVwplaaU+stFfu6plFrc8PPdSqmoC372bMP3jyulxtozpxBCCCGEs7JbsaaUcgU+BK4HugF3KKW6/eKy+4HzWuuOwN+ANxqe2w2YAsQBicA/Gv48IYQQQohmxZ4ja/2BNK11uta6DlgETPzFNROBeQ1fLwFGKqVUw/cXaa1rtdangbSGP08IIYQQolmxZ7HWHsi+4HFOw/cueo3W2gSUAq1tfK4QQgghRJNnz2JNXeR72sZrbHkuSqkHlVJJSqmkoqKi3xBRCCGEEMK52bNYywHCL3gcBuT92jVKKTcgADhn43PRWn+itU7QWicEBQVdxehCCCGEEM7BnsXaXiBWKRWtlPLAumBgxS+uWQFMb/h6MvCj1lo3fH9Kw2rRaCAW2GPHrEIIIYQQTsluJxhorU1KqUeB9YAr8LnWOlUpNQdI0lqvAD4D5iul0rCOqE1peG6qUuob4AhgAh7RWpvtlVUIIYQQwlkp60BW45eQkKCTkpKMjiGEEEIIcVlKqWStdYIt18oJBkIIIYQQTkyKNSGEEEIIJybFmhBCCCGEE5NiTQghhBDCiTWZBQZKqSIg0wFNtQHOOqCdxkr65/Kkjy5N+ufypI8uTfrn8qSPLs0R/ROptbZpk9gmU6w5ilIqydbVG82R9M/lSR9dmvTP5UkfXZr0z+VJH12as/WP3AYVQgghhHBiUqwJIYQQQjgxKdau3CdGB3By0j+XJ310adI/lyd9dGnSP5cnfXRpTtU/MmdNCCGEEMKJyciaEEIIIYQTk2LtVyilEpVSx5VSaUqpv1zk59copfYppUxKqclGZDSSDf3zpFLqiFIqRSm1SSkVaUROo9jQPw8ppQ4ppQ4opbYppboZkdNIl+ujC66brJTSSimnWZnlCDa8hu5RShU1vIYOKKX+x4icRrLlNaSUuq3hvShVKbXQ0RmNZMNr6G8XvH5OKKVKjMhpJBv6KEIptVkptb/h82ycETnRWst/v/gPcAVOATGAB3AQ6PaLa6KAeOBLYLLRmZ2wf0YAPg1fPwwsNjq3k/VPiwu+ngCsMzq3s/VRw3X+wBZgF5BgdG5n6h/gHuDvRmd18j6KBfYDLRsetzU6tzP1zy+unwl8bnRuZ+sjrHPXHm74uhuQYURWGVm7uP5AmtY6XWtdBywCJl54gdY6Q2udAliMCGgwW/pns9a6quHhLiDMwRmNZEv/lF3w0BdobpNHL9tHDV4B3gRqHBnOCdjaP82ZLX30APCh1vo8gNb6jIMzGulKX0N3AF87JJnzsKWPNNCi4esAIM+B+f6PFGsX1x7IvuBxTsP3hNWV9s/9wFq7JnIuNvWPUuoRpdQprMXIYw7K5iwu20dKqd5AuNZ6lSODOQlbf8cmNdyaWaKUCndMNKdhSx91AjoppbYrpXYppRIdls54Nr9PN0xTiQZ+dEAuZ2JLH70E3KWUygHWYB2BdDgp1i5OXeR7zW3k41Js7h+l1F1AAvCWXRM5F5v6R2v9oda6A/AMMNvuqZzLJftIKeUC/A14ymGJnIstr6GVQJTWOh74AZhn91TOxZY+csN6K3Q41pGjfymlAu2cy1lcyefYFGCJ1tpsxzzOyJY+ugP4QmsdBowD5je8PzmUFGsXlwNc+K/UMAwa+nRSNvWPUmoUMAuYoLWudVA2Z3Clr59FwE12TeR8LtdH/kB34CelVAYwEFjRjBYZXPY1pLUuvuD36lOgr4OyOQtbfs9ygOVa63qt9WngONbirTm4kvehKTS/W6BgWx/dD3wDoLXeCXhhPTfUoaRYu7i9QKxSKlop5YH1hbzC4EzO5LL903AL62OshVpzmicCtvXPhR8YNwAnHZjPGVyyj7TWpVrrNlrrKK11FNZ5jxO01knGxHU4W15DIRc8nAAcdWA+Z2DL+/QyrIudUEq1wXpbNN2hKY1j0+eYUqoz0BLY6eB8zsCWPsoCRgIopbpiLdaKHJoSKdYuSmttgv+vvbsLsaIO4zj+/WlYyGpplJVlG0kvYhGSldmLxSIIFdobbRBKWYSImaZdVKAS5YVU0IXZhVaUL5kJFqRWbJGmoamradFN0ctVaphpdRFPF/McPKy7crLYGd3f5+bMPjPzn//8d8/h2f+cmYepwDqKD8C3I2K3pHmS7gCQNDKvYd8DLJK0u7wed69GxofismcTsDJvC+8xyW6D4zM1HyWwA5gBTCypu6VocIx6rAbHZ1r+DbVTfOdxUjm9LUeDY7QO2CdpD9AGzIqIfeX0uHv9i/dYK7A88nbHnqTBMZoJPJzvs2XApDLGyhUMzMzMzCrMM2tmZmZmFeZkzczMzKzCnKyZmZmZVZiTNTMzM7MKc7JmZmZmVmFO1sysEiRNkBSSLuvGY06X1Le7jnc8JI2RdH3Z/TCz8jhZM7OqaAU2UDyYsrtMB0pP1iSdcozVYwAna2Y9mJM1MyudpCZgNEVpl/s6rJstaZekdknzMzZU0kcZ2ybp4ozPkrQli5vPzVizpG8kvV5X9LyvpGnAeUCbpLbcdqGkrfmw2bl1ffhe0tw81q7a7J+kJklLMrZT0l0ZHytpU26/Ms+v4zl/Iuk5SZ8Cj0m6XdIXkrbnuQ2S1Aw8CjyeD5e+UdJZklbleW6RNPr//W2YWdUc6785M7PuMh5YGxHfStovaUREbJM0LtddGxGHJQ3M7d8C5kfEakmnAb0kjaWo+3gNRYHmNZJuoigXcynwUERslLQYmBIRCyTNAG6JiL3Z7lMRsV9Sb+BjSVdGxM5ctzciRkiaAjwBTAaeAQ5ExBUAkgZkWaOngZaIOCTpSYoqFfM6Oe8zIuLm2r7AdRERkiYDsyNipqRXgN8jYkFutxR4MSI2SBpC8fT1y//j+JtZhTlZM7MqaAVeyuXl+fM2oAVYEhGHATKR6gcMjojVGfsTitksYCywPdtpokjefgB+jIiNGX+TojzTgk76ca+kRyg+G88FhgG1ZO3dfP0SuDOXW6ibCYyIXyXdlvttlATQh67rLq6oWz4fWKGi5mcf4Lsu9mkBhmXbAP0l9YuIg11sb2YnOCdrZlYqSWcCtwLDJQXQGwhJsylmyDrWxBOdE/B8RCzq0H5zJ20cVWdP0kUUM2YjM+l6jaJoc81f+fo3Rz47u+rfhxHR2kU/6x2qW34ZeCEi1kgaA8zpYp9ewKiI+KOB9s3sJODvrJlZ2e4G3oiICyOiOSIuoJhVugFYDzxYu2NT0sCI+A34SdL4jJ2a69fltk0ZHyzp7DzGEEmjcrl2IwPAQaBfLvenSJ4OSBoEjGug7+spCkGTxxwAbAZGSxqasb6SLmmgrdOBn3N5Yl28vo+dHfOqBto2sxOYkzUzK1srsLpDbBVwf0SsBdYAWyXtoJj5AngAmCZpJ/A5cE5ErAeWApsk7QLe4UiS8zUwMbcfCCzM+KvAB5LaIqKd4hLqbmAxULtseizPAgMkfSWpneL7b78Ak4BlebzNQCOPI5kDrJT0GbC3Lv4eMKF2gwHFJdyr84aGPRQ3IJjZSUwRR10NMDM7aeRl0PcjYnjJXTEzOy6eWTMzMzOrMM+smZmZmVWYZ9bMzMzMKszJmpmZmVmFOVkzMzMzqzAna2ZmZmYV5mTNzMzMrMKcrJmZmZlV2D+Tm6vMmW9ZIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00909091 0.05670447 0.04090909 0.        ]\n",
      " [0.03636364 0.09335408 0.07727273 0.        ]\n",
      " [0.07272727 0.14629259 0.10909091 0.        ]\n",
      " [0.10454545 0.19167265 0.13636364 0.        ]\n",
      " [0.16818182 0.25233309 0.19090909 0.        ]\n",
      " [0.21818182 0.29859214 0.25454545 0.        ]\n",
      " [0.28181818 0.35540211 0.3        0.        ]\n",
      " [0.34545455 0.39316675 0.37272727 0.        ]]\n"
     ]
    }
   ],
   "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]\n",
    "\n",
    "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, False)\n",
    "\n",
    "    ## Human error rate - Jotain väärin viel'\n",
    "    # Get judges with correct leniency as list\n",
    "    correct_leniency_list = test.loc[test['acceptanceRate_R'].round(1) == r / 10, 'judgeID_J']\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",
    "    # Get their failure rate\n",
    "    failure_rates[r-1, 1] = np.sum(released.result_Y == 0) / np.sum(test.judgeID_J.isin(correct_leniency_list))\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, False)\n",
    "    ## Dec tree\n",
    "    failure_rates[r-1, 2] = contraction(test[test.decision_T==1], 'judgeID_J', 'decision_T',\n",
    "                                   'result_Y', 'B_prob_0_tree', 'acceptanceRate_R',  r / 10, False)\n",
    "    \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')\n",
    "#plt.plot(np.arange(0.1,0.9,.1), failure_rates[:,3], label = 'Tree')\n",
    "\n",
    "\n",
    "\n",
    "plt.title('')\n",
    "plt.xlabel('Acceptance rate')\n",
    "plt.ylabel('Failure rate')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "print(failure_rates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.05670446964643095\n",
      "0.09335407868415203\n",
      "0.1462925851703407\n",
      "0.19167264895908112\n",
      "0.25233309404163673\n",
      "0.29859214120613575\n",
      "0.3554021121039805\n",
      "0.3931667516573177\n"
     ]
    }
   ],
   "source": [
    "for r in np.arange(1, 9) / 10:\n",
    "    # Get judges with correct leniency as list\n",
    "    correct_leniency_list = test.loc[test['acceptanceRate_R'].round(1) == r, 'judgeID_J']\n",
    "    \n",
    "    # Released are the peopöe they judged and released, T = 1\n",
    "    released = test[test.judgeID_J.isin(correct_leniency_list) & (test.decision_T == 1)]\n",
    "\n",
    "    # Get their failure rate\n",
    "    print(np.sum(released['result_Y'] == 0) / np.sum(test.judgeID_J.isin(correct_leniency_list)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  26.,  284., 1436., 4093., 6920., 6846., 3888., 1230.,  252.,\n",
       "          25.]),\n",
       " array([-3.77161989, -3.01127888, -2.25093788, -1.49059687, -0.73025587,\n",
       "         0.03008514,  0.79042614,  1.55076715,  2.31110815,  3.07144916,\n",
       "         3.83179016]),\n",
       " <a list of 10 Patch objects>)"