xgrg/Interaction models.ipynb

## Interaction models.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sum_sq</th>\n",
       "      <th>df</th>\n",
       "      <th>F</th>\n",
       "      <th>PR(&gt;F)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>161.015057</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.725728</td>\n",
       "      <td>0.194616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(apoe)</th>\n",
       "      <td>919.286352</td>\n",
       "      <td>2.0</td>\n",
       "      <td>4.926366</td>\n",
       "      <td>0.010907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>910.623151</td>\n",
       "      <td>1.0</td>\n",
       "      <td>9.759881</td>\n",
       "      <td>0.002890</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(apoe):age</th>\n",
       "      <td>1702.150754</td>\n",
       "      <td>2.0</td>\n",
       "      <td>9.121660</td>\n",
       "      <td>0.000394</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <td>8.070552</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.086499</td>\n",
       "      <td>0.769825</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>4945.042466</td>\n",
       "      <td>53.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  sum_sq    df         F    PR(>F)\n",
       "Intercept     161.015057   1.0  1.725728  0.194616\n",
       "C(apoe)       919.286352   2.0  4.926366  0.010907\n",
       "age           910.623151   1.0  9.759881  0.002890\n",
       "C(apoe):age  1702.150754   2.0  9.121660  0.000394\n",
       "sex             8.070552   1.0  0.086499  0.769825\n",
       "Residual     4945.042466  53.0       NaN       NaN"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys, json, roistats, nilearn\n",
    "from nilearn import plotting as nipl\n",
    "from glob import glob\n",
    "from roistats import collect\n",
    "from nilearn import datasets \n",
    "import logging as log\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from roistats import plotting\n",
    "from roistats.contrasts import genotypes as cg\n",
    "from collections import OrderedDict\n",
    "from statsmodels.formula.api import ols\n",
    "from matplotlib import pyplot as plt\n",
    "import logging as log\n",
    "import statsmodels.formula.api as smf\n",
    "from IPython import display\n",
    "sns.set_style('darkgrid')\n",
    "log.basicConfig(level=log.INFO)\n",
    "import random\n",
    "\n",
    "N=70\n",
    "#age_0 = np.random.uniform(low=45,high=75, size=(N,))\n",
    "#age_1 = np.random.uniform(low=45,high=75, size=(N,))\n",
    "#y = beta0 * age + eps\n",
    "data = pd.DataFrame(data, columns=['y','age','sex','apoe'])\n",
    "#data['y'] = data['age'] + 1\n",
    "#data['age'] = data['age'] - data['age'].mean()\n",
    "#data['y'] = data['y'] + np.random.uniform(low=0, high=10, size=(N,))\n",
    "formula = 'y ~ age + sex + C(apoe) + C(apoe):age + 1'\n",
    "model = ols(formula, data)\n",
    "ts = model.fit()\n",
    "from statsmodels.stats.api import anova_lm\n",
    "anova_lm(ts, typ=3)\n",
    "#ts.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fe0f0b8b750>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAFjCAYAAAATuV17AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXecVNX5/9/33ultd7ZSlrIgqIAFMSAmMQZjw1gQLLGg\ngKJG0ViiorEX7AY1MbZYsAv2Bmry1djAnyaiUqWzvU3v997fH7Mzy7bZwi47LOf9evHSnblz72dm\n757PPOc5z3MkXdd1BAKBQCDoInJfCxAIBALB7okwEIFAIBB0C2EgAoFAIOgWwkAEAoFA0C2EgQgE\nAoGgWwgDEQgEAkG36HMDue666zj00EM5/vjj0495vV5mz57N0UcfzZw5c/D7/ennbr/9do466ihO\nPPFEVq9e3ReSBQKBQEAWGMjJJ5/MU0891eyxxx9/nMmTJ7N06VImTZrEY489BsCnn37K1q1bWbZs\nGbfeeis33XRTp6+zfPnyHtW9M2STFhB6OkLoyYzQk5ls09OT9LmBHHzwwbhcrmaPffLJJ0ybNg2A\nadOm8cknn6QfP+mkkwA44IAD8Pv91NbWduo6K1as6EHVO0c2aQGhpyOEnswIPZnJNj09SZ8bSFvU\n19dTUFAAQGFhIfX19QBUV1czYMCA9HHFxcVUVVX1iUaBQCDY08lKA2mPtrquSJLUB0oEAoFAYOhr\nAW2Rn59PbW0tBQUF1NTUkJeXByQjjsrKyvRxlZWVFBUVtXmO5cuXNwsd582bRzAY7V3hneT3vz8h\na7SA0NMRQk9mhJ7MOJ1OHn744fTPEydOZNKkSX2oqOfICgNpGVlMmTKF119/nblz5/LGG29wxBFH\nAHDEEUfwwgsvMHXqVP73v//hcrnSU10tmTRpUqtfUigU65030EWKigZkjRYQejpC6MmM0JOZc889\nt68l9Bp9biBXXnkly5cvx+PxcPjhhzNv3jzmzp3LZZddxpIlSxg0aBALFy4E4De/+Q2ffvopRx55\nJFarlQULFvSxeoFAINhzkfakdu41Nf6OD9oF2GymrPqGJPRkRujJjNCTmcJCZ19L6DV2qyS6QCAQ\nCLIHYSACgUAg6BbCQAQCgUDQLYSBCAQCgaBbCAMRCAQCQbcQBiIQCASCbiEMRCAQCATdQhiIQCAQ\nCLqFMBCBQCAQdAthIAKBQCDoFsJABAKBQNAthIEIBAKBoFsIAxEIBAJBtxAGIhAIBIJuIQxEIBAI\nBN1CGIhAIBAIuoUwEIFAIBB0C2EgAoFAIOgWwkAEAoFA0C2EgQgEAoGgWwgDEQgEAkG3EAYiEAgE\ngm4hDEQgEAgE3UIYiEAgEAi6hTAQgUAgEHQLYSACgUAg6BbCQAQCgUDQLYSBCAQCgaBbCAMRCAQC\nQbcQBiIQCASCbiEMRCAQCATdQhiIQCAQCLqFMBCBQCAQdAthIAKBQCDoFsJABAKBQNAthIEIBAKB\noFsIAxEIBAJBtxAGIhAIBIJusUcZiKwm+lqCQCAQ9Bv2KAORyrYjRyN9LUMgEAj6BXuUgRCPI1WU\nofi9SFJfixEIBILdmz3LQABUDWpqUGpqkHS9r9UIBALBbouhrwVkYsqUKTgcDmRZxmAwsHjxYrxe\nL5dffjllZWWUlJTw17/+FafT2bUT6zq614MSi6IVFaMZjL3zBgQCgaAfk9URiCRJLFq0iDfffJPF\nixcD8PjjjzN58mSWLl3KpEmTeOyxx7p9fj0cRiovQ4mEekqyQCAQ7DFktYHouo6mac0e++STT5g2\nbRoA06ZN4+OPP965i8TjUFGB4vWIvIhAIBB0gaw2EEmSmDNnDtOnT+e1114DoK6ujoKCAgAKCwtp\naGjY+QtpGtTVolRXI+lax8cLBAKBILtzIC+//DKFhYXU19cze/ZsSktLkToZJixfvpwVK1akf543\nbx4Wiynzi2JhqK+G4mJ0k3lnpGfEaFSw2TrQsgsRejIj9GRG6MlMy7Fo4sSJTJo0qQ8V9RxZbSCF\nhYUA5OXl8bvf/Y6VK1eSn59PbW0tBQUF1NTUkJeX1+ZrJ02a1OqXFInEOr5oJIYUCKMXFqLaHDv9\nHtrCZjMRCnVCyy5C6MmM0JMZoSczbY1F/YWsncIKh8MEg0EAQqEQn3/+OaNHj2bKlCm8/vrrALzx\nxhscccQRPX5tPZGAqkqUhnokxFJfgUAgaIusjUBqa2u55JJLkCQJVVU5/vjj+dWvfsW4ceP405/+\nxJIlSxg0aBALFy7sHQGaDvV1yaW+hUVostI71xEIBILdFEnX95xquvqvv+veC01m9KIiNLOlR3Rk\nW4gt9GRG6MmM0JOZwsIu1qntRmTtFFZWEYsiV5ShBP19rUQgEAiyBmEgnURXNaiqwlBfK1qgCAQC\nAcJAuoauozc0oFRViNbwAoFgj0cYSDfQg0HRGl4gEOzxCAPpLvF4so+WaA0vEAj2UISB7AyaaA0v\nEAj2XISB7Cyp1vAVZciJeF+rEQgEgl2GMJAeQg+Hkcu3i9bwAoFgj0EYSA+ixxNNreEztEDRdZ2o\nJhLwAoFg90YYSE+Tag1fk7k1fFWwkvpoLTqifbxAINg9EQbSG+g6us+HXF6GHG+7pYKmazSEG9jm\n30ow4RcruQQCwW7HHmUgxk8+gl25WioSQS4vQwkF2j0krsapClZRGawgoYskvEAg2H3YowzE+Zdr\ncVxxKfL2bbvsmp1pDa/rOoFYgO3+bXjjHnTRQl4gEOwG7FEGEvzTlRhXfk/OmadhefpJiO2ijp2p\n1vBVlcia2u5hqqZSG6yhPLCdiBbeNdoEAoGgm+xRBhI97Qw8Ly0m9qtfY3v8UXLOPh3D/1vR8Qt7\nCD0QQCorQ4pEkGg/6RFJRKjwl1MXrUGjfcMRCASCvmSPMhAAvaiI4B1343/gYVATuOZdhP2WG5Dq\n63aNgFgUysqwhKPIUvsfv6ZreMIetvu3EYz7QUxrCQSCLGOPM5AU8cmH4n3+VcKzzsP08TJyTp+O\n+fXXQN0F3/g1Db26CqvHhyGDicAOSfZQBXEtezbJEQgEgj3WQACwWAjPvQjv86+gjt4H+7134Zo7\nC2Xtml6/tK5pqPV1WKprMXVQCqKjE4wFKQtsoyFWL5LsAoEgK9izDaQRbdhw/A8/SuDm25ArKnDN\nPhvbg/dBsP3ltz1FIuDHWFmFJa4idVAMomoa9aE6tge2EVFFyxSBQNC3CANJIUnEjp6K9+UlRE86\nGfNrL5N7+gxMHy/r9doRNRpBrqjEFoygyB3/SmKJKBWBcmoj1ai62NhKIBD0DcJAWqC7XIT+PB/f\nE8+g5efjuGE+jsvn9XrtiKYm0KqrsdR5MWRYoZU+XtfxRryUBbaLJLtAIOgThIG0gzp2HL4nnyV4\n+VUYf1iZrB355xO9Wjui6xqqpx5LVS3mTvrBjkn2mB7tNW0CgUDQEmEgmTAYiJ76BzwvLyH268Ow\nPfGPna4d2RzazpAPfsVN6x9me6SyzWMSoQCGsgps0ThyJ5pkpZLs5f7t1MfqRINGgUCwS1Buvvnm\nm/taxK4ivL2iey+024lP+R2Jcfth/OI/WF95CXn7NhL7HwBWW5dO5TTYqVMbeKXsfV6qeI+yaBUj\nrUPJMTqbHadrKgRDmCQFzWLp1MorXdeJxMMEEyGMsoJJMXVKk9GoEI9nT8Gi0JMZoScz2abHbjf3\ntYReQxhIF9BKhhA9YRrIMua33sD8xhJ0hwN19D7QieQ3gCzJHDvoMA5zTiCuxXmn5v94qfI9tkUr\nGWkdQu6ORqLr6JEwpoSKZLWidrJjr6qpBONBYmoUi8GCLCkZj8+2PzihJzNCT2ayTU9/NhBJ1/ec\nzbzrv/6ux84lb9mM/d67MH77DYkxYwlefR3q3vt06rVms5ENdZvQNJWaWD3Plr/FkqqlxLQExxb8\nmvNKZjDcOrjZaxSzBbWwgKjJQFd+ZYqs4Lbm4TLmtNs+xWYzEQplT5Gi0JMZoScz2aansNDZ8UG7\nKcJAdgZdx7TsA2wLH0TyeojOOJXQ3IvA7sj4sh0NJEVdzMNz5W/xatWHxLQ4Rxf8kvNLTqHUWpI+\nRlYUpPwCIg4baobNqtrCYrCQbyvAIltbPZdtf3BCT2aEnsxkmx5hIP2EHjeQRiSfD+tjf0tOaeUX\nEPzTlcSn/I72dolqy0DSGuMeFpW/zSuVHxLRohyVnzSSkbYhyWtJEorLRTTPTbyLS3dlScJpzsFt\ndqNIhvTj2fYHJ/RkRujJTLbpEQbST+gtA0mh/PQj9nvuxLBuLbFDDiV05dVoJUNaHZfJQFLUx708\nX/4OL1e+T0SLcmT+ZM4vOYW9bMOS17JaUQsKiBqULrc2MSpG8ix5OIxOQMq6PzihJzNCT2ayTU9/\nNhCRRO9B9KIiosefiJabi+X9d7EsfhV0jcTY/UBpSmQbDAoNYU/GXIZVsTApd39OLj4SRZJ5v+Yz\nXqx8jw2hrZRaS3BLduRQCKPJhGoydclENF0jFA8RUSOYjWYsJnNWJR2zLQkq9GRG6MmMSKL3E3o7\nAtkRqboa20MPYP7kI9ShwwhedS2JX0wEOheBtMQT9/NCxTu8VPkeQTXMlLxDmFtyCvs4RiLn5BJ1\nu0h0oxpdlmQG5BZi0RxIHXQG3lVk2zdIoSczQk9m+nMEIgyklzF+/SW2++5GKdtO9OhjCV16OaaB\nA7psICl8iQAvVLzLixXvElBDHO6eyAVDTmVs4X4kCvOJyp1c67sDFosJLSFRYM3HZrDv0m3j2yLb\nBgChJzNCT2aEgfQT+sJAAIhEsC56BsuiZ9DNZuLzLmPtbyagdX2sT+NPBHmx4l1eqHgXvxrkN+5f\ncNHwM9l3xCFELGa0LqzSslhMRCIxJEnCaXbiNudj2CHJvqvJtgFA6MmM0JMZYSD9hPoV/03uT95H\nyFs2Y7/vboz/bwWhvUexbd5cwqNG7tQ5/YkgL1W+x/Pl7+BXgxzmPpiLxs5lr2G/INFJE0kZSAqD\nbGisHXFBJxo79jTZNgAIPZkRejIjDKSfULetGtnrQQ+FQOujflG6ju3fH6HcexcGn5/aE46lYuYf\n0Oxda4nSEn8iyCuVH/B8xdt4EwF+XTCR88dfwqjCsR2+tqWBQNI2LEYr+dYCzLJlp7R1lWwbAISe\nzAg9mREG0k+oqfEjSSBFo8g+LwT86OquNxKz2cimLT9Q/M9FFLy7lIQ7l7ILZ+M57NB2a0c6SyAR\n4pXK91nUaCSHFh/CrAMuYFzB/u2+pi0DSSFLMk6zq1XtSG+SbQOA0JMZoSczwkD6CTU1/mY/y4kE\nkt+H5Peix3fdxkxms5GN9ZtQVRXbmvWUPPwYtp834ptwINsvmUts0ICdvkZQDfNq5QcsKn+bhoSP\nSQMnM3vc+exfeECrYzMZSIqWtSO9SbYNAEJPZoSezAgD6Se0NJAUkq4h+/3g9UI82ut7M1ksJvyh\nIEE1SCgWRI3FyH/nAwY++yJSQqXq9JOpPmUausm409cKqWFeq1rKcxVvUx/z8IsBk5gz7nwOKBrf\nTE9HBgIgIWExWsiz5GE12HpttVa2DQBCT2aEnswIA+kntGcgKSRADjYaSSTSa1vZ7jhg60AkESIY\nD5KoKmfgo0/i/uxLIiWD2T5vLoED9+uRa4bVCIvr/8Uz216nPlrPhOJfMGe/uYwvOqjTBpIiuVrL\nhduc1yurtbJtABB6MiP0ZEYYSD+hIwNJIUkgh0NIXm+vJNzbG7ATuko4EUL/4lMGPvR3zBVV1E85\njPK555Jw5/bItWOyxmu+//DP9c9TH6njoKIJXDThYsa5W09tdYRBNpBvzcdpcvWo12bbACD0ZEbo\nyYwwkH5CZw0kRW8l3DvzjT8a9mF+5klyX3oF3WSifNZZ1E09sllLlO4iyTIxp42Xaz7iuVVPUxuu\n5cCig5gz7nwmFP8CqQuJfAkJu8lOvrWwx6KRbBsAhJ7MCD2ZEQaSZXz22Wfceeed6LrO9OnTmTt3\nbqde11UD2ZGeTLh3ZcpI2rIJ670LsHz7bY/VjqQwOJz4c+y8uf1t/vn9U9SGazig8EBm7zeXXxRP\n7JKRGGQDedY8nD1QO5JtA4DQkxmhJzPCQLIITdM4+uijeeaZZygqKmLGjBk88MADjBzZ8aC6MwaS\noicS7l3NOaDrmD76EOvCB5A9HupPOI7yc05HtbXe26OryGYzpkEDqdaivLX+DRateoaacDX7FRzA\nnP3OZ+KAQzptJBJgNdrIs+bvVO1Itg0AQk9mhJ7M9GcD2e268X7//fesW7eOs846C0VR8Pv9bNq0\niQkTJnT42h65qSQJ3WyBnBwkkwlUNfmvCxgMColEF14jSagjRxE7YRpSIIDzzTcp+OhT5IGDiQwd\n3OV27juiqyqGaAQTRvYZPIFpo2ZQYC3gy/LPWbL+NZZXfEWhrYgSR0mnjCSuxQnEAuiShtlgRqLr\nDRqzrZuq0JMZoScz/bkbb3a0X+0CVVVVDBw4MP1zcXEx1dXVu1yHDqh2J9rgEhg0CMluh240MuzS\nNZ1OQn++Ft+Tz6IXFDDgltsZfcPdFDVEsRltyJ3cl73VeTUVta4WS20dTtnM9NGn8trxb3L1L66j\nNlzDFf83j/OWncOXZZ93ajvdymAlf1w2lzPfPYVr/vMnqoK920ZfIBD0DbudgWTbjJuug2qxoQ4c\nhD6oBMnphG4O5J1FHTMW31PPEbzizxh/WMmAmedQ8tKbDDAWkGfNx6SYu5S/AEDXUX0+jBVVWBLJ\n6GHaqOm8evybXDvxehoiDVz56WXMWXoOn5d9lvH3sPDb+ykLbKfMX8aK8hXc/vVNRLTwTr5rgUCQ\nbfRdy9VuMmDAAMrLy9M/V1VVUVRU1Oq45cuXs2LFivTP8+bNw2Yz9a44qwncLqRIBMnrQQ8E21wC\nbDDIWCw9oOXsmUSOPgbj/fdge+IfWJZ9SGz+9RRMOoSIGiEUDxKKh1A7ahsvyUSjRsxmUHQVU10N\n1vwCok47mm7itHGnc/KY6bz78zs89b/H+fOnl7Nv/hjmjr+Q3ww5vJVZ+RO+dDSkobHRt5H6eA0u\ns5M8Sz4GOfNtZzQqvf+76gJCT2aEnsy0HIsmTpzIpEmT+lBRz7HbJdFVVeWYY47hmWeeobCwkFNO\nOWWXJtG7ghyPJZcA+33NlgB3OYneCdrad0TPy0dDI5wI4Y8GiGuxNiMHRTFSVpZA15OrhM0WCYsZ\nrAUu4vm5xKWm1yS0OB9u+oBnfnqSskAZo917M3vc+RxW0mQk1//nGsoC25GQ0NEZ7Cjhjl/fndQp\nG8iz5mdsiZJKglYEyrnhi/k0ROpwW/K57ZcLGOgY1KOfW2fItqSs0JOZbNPTn5Pou52BQHIZ7x13\n3IGu68yYMWOXLOPdGeREHMnnQ/J50VW1VwwEaLXvSPiiS4ieeHLSFSSIqTGCicb2KTtEJSkD2bHT\nvUSyDsbhtqOU5KFZFAxGLV0wmNDiLN38AU//+BRlge2Myh3N7P2SRlIbruWhb+/HE/WQY87lsglX\nUmQr2uHcEjaTjXxLAUa59TfF1ABw5pvnsHz9ZjRVQlZ0Dhk1nOdPerbnP7cOyLYBSejJTLbpEQbS\nT+grA0khqwkknw9LLEwkEOq962zdgv2eBRi//YbEvmMJXjMfde9908/r6ITjIYKJINFEFElWWhnI\njhiMBqTCAmJOGyaThsWiYzLp6DoktASL173KEysfJZQIYTc6uOTAyzhhr5OQO9giV5Flcixuck1u\npB2ikdQAMOauEwmqvvTjdsXFqmvfAtil0Um2DUhCT2ayTU9/NpDdLom+O6MpBlR3HgwdCgWFSMbe\nSUFpQ4fhf/hRAjffjlxZgWv2TGwP3gvBANAYARjtFNqKKLYPwG3NwWwwIrczpZSIJ1ArqzHXegj7\noKoKKipkfD4ZNW7kh5rvKXEMYaBtEHE1xt3f3MHZ75/Ox1uWZdwZUdU06kN1lAW2t51kD+enlyjr\n6BDOTz91wxfz2ezdhC/mZ7N3Ezd+MX8nPjGBQNAdhIH0AbqioObkog4ZBoVFYNz5rrutkCRiRx+L\n95XXiZ40HfNrr5B7+nRMHy9rahKpg1E24jK5KbQMpNBeRI7FgVFpbSWaphGrq8dWV4MVnXhcx+PR\nqa6Gar8XTZNxmnIodY2g1FWKqqvc8MV8znr/ND7avDRjIj+aiFDhL6MuWoOqN1X5j6+6F3NgJHLc\nhTk4gvFV96afa4jUpXMukiRRH6nrsY9OIBB0DmEgfYguyaiuHLQhw6C4GEzmHt9qI1078sQzaPkF\nOG6Yj/PyS5C3bWt2nKZJKJoFp5JPkWUQhfZ87EYzSosVVrFAEGNlJTlqHEWW0XVwKDmoqk4iDrE4\nFCojOWj1W4zddi81NRI3fnkdZ71/Gks3f9CukWi6jifsYbt/Gw2RenR07rupgMNrX+CgH5bxm+oX\nue+mgvTxbkt+ekGAruu4LfltnlcgEPQeIgfSB7Q3R5tuJ9/ggVi059vJqyrm11/D9o+/QyJOeOYs\nImefi8HmYPPmeKsVx7IMqhQjrAUJRYPEVTVd8y4rMga3m4DdSUWoikUb78OX8OI05BD9ZD6+8oFY\nrWC2aCh7f0jN8EfZ5NvAUNcwZo09j98NO6rd5bwWiwktIZFvycNudLT5MVQEyrnxi/nUixxInyP0\nZKY/50CEgfQBHd3gEjpyKAgeT6/sSyLV1GB76AHMHy9DHTqM4FV/YfPACe12rZckkCSdGCFCiRCR\neISEpoEEJoeTqDuPkN4k89ZbzYTDTZGL1apz401hvqn/N4u3PM7mwHpKHEOZvd8cjhx2TCsjSa1S\nkyQJW2NvLZPUd+0gsm1AEnoyk216+rOB7Ha9sHaGbLmpOu7VI6EbTeByIVktSJoGiR7cctduJz7l\nd8T32x/tky9wLnmeL9/x8Nj6I1lfm0c0JuOyq5hNTcal6xKybsKq2LGZ7JgMBnRU4uEIxkgIo9VK\nQlHQdfjuO4VAQEKSkqaSl6cz+RCNwdYRHFE8nWH2UazxruTtjUtYunEpZsnOSPdeKI3Fhzv2Cour\ncQIx/0711tpZsq23ktCTmWzT0597YYkIpA/o6jckSQIpHEb2enp8g6uIL8rPNyzmkG8eIYyV+fqd\nPMYFaCgMLoqyb2mIfUeEGTMixL4jQuQ41CZNkk6cKOFEgKgWQ8/NIWBzUlGts2iREZ9PotRQyYOh\nu3BGPMRcOfz3oquIFBah6Rrf1n/K4q1PsDm4hmJLCaeOmMPUEcfhzrGi63Fatjo2ygbc1jwcRlez\nZb+9TbZ9oxV6MpNtevpzBCIMpA/o7g0uSSBFIsg+L3og0GNGYjCYKP/yZwY/sQDnjyuoHrg/z+x/\nDx95JrJ6k42y6qZvUIMKo+w7IsS+pWHGjgyxT2mIvBwVnQRRKUTULONzOQnFkwtwJ996NfaK7aTC\nkcCgIXx9w93p8+m6znf1n7F42+NsDKymyDyYU0acz2GFx+J0mLCYdSyWpFmlMBss5Fnc7eZHepps\nG5CEnsxkmx5hIP2E3d1AdkSOx5IRScAPO7lTosFgSibRVZ3cLz5k8HMPYPA1UHvMaVScehENWg6r\nN9pYvcnK6o021myysa2quansU5qMUsaMDLLvvkFMw83UqAkOuWwOxlAgfWzc5uDT+x5rpUHXdf7b\n8DmLtz7OhsBPFJoHctKQ2fy2+ASMihGLRcJiAbM5WcQIyfyI25K3U3uPtKSiQuKGG8w0NEi43Tq3\n3RZl5EhjsrVKG88NHLjr/3yybYAUejIjDKSf0J8MJIWciCf7bfl86F3clyRF2kAafUgJ+hnw8iMU\nfLSYRG4+ZedcheeQI5NRRCO+gMKazUlDWbXRxuqNzU1lYGGMMfsmmFj9Fgf5/4/9zT9SoNQTHFjC\nVzfe064WXdf5KbCclzb8nZ8DP1JgHsBJJbM5vPgEjLIJWU6uDrNaJcxmHYtZxu1wkGtyY5B2vp7m\nvPMsbN4sp/M3paUaL7ygEQrF2nzuiSciO33NrpJtA6TQkxlhIP2E/mggKeREAsnnRfL70LuYcG9p\nICmsP//EkCfvwLZpDb4DJrN99jXEBgxt9zz+oJKOUlZttLF6k41tlTtEKuZqhu0PQ8YmGDE6zIhR\nIXLcTVrr6iSee85IMChjs6scfPKnLPM8xnr/D7hNheQY87AqDnKMbmaOuIp8cxGyDCaThNOhUOhy\nUejIRd6JRPv06VZ8viajdLl0PvhAJRSKtfnckiW7vk19tg2QQk9mhIH0E/qzgaSQ1QSS34/k83R6\n7/b2DAQATaVg6WsMfOVvSIk4VSfNpvrEc5OrxDqBP6iwriKHnyoL+GmtiTVrZLZuVdLP5xfGKN0r\nROnoMCt/ihGORzGYNHRdZ8AAjcsui7LS8zUPrb2OQMKLIhnINeazl2MMV4y5v/l7l8FqMlDkyiPf\nnoPRCGZz8/xJR7SMMoYP13jxxbYjkOHDNZ58UkQgQk9mhIH0E/YEA0khaSpyINC4d3ssYy1JRgNJ\nHVNfw+Dn7sf91TIiA4eyfc58Avt1fk8D2ShDfiERiwu/H9avl1m1RuPH1TrrVhsp325G15Pf7hVj\nHLM9gj0nzKzz/YwYFeKhstnUx6rxxGqIaGEMkpGzSv/EEcXTMCnNcyCSBBajmTyrG4vkwGyWsdmS\nyfhUE8j2qKiQuPFGM/X1bedAWj7XlzkQkZNpm2zTIwykn7AnGUgKCZBDgYxFiZ0xkBTO77+i5KkF\nmKu2U//LYymfeTmJ3IKOXwhIsoSc6yLmKiCuJaeZJAmQVTyBILfcqbHhZyMhn4VoyEI8YiLV28Xo\nqsMweCXmkh+RRn5EfMDXBKgh11jACSXncOSA6a2MRJYk7GYruaY8FN3a+F4l7PZkMt5iAYMhs6FA\n9g1IKT0iJ9M22aZHGEg/YU80kBSSBHIoiOTxoEfC7Ni7vSsGAiDFohS/9TRFbz6NbjJTfvol1B05\nHWSl4xcDis2Cml9AVLY2G7wbGiSefVHFFwmQkL0cOzWOr8HCxnU21qyV+GFVgnDVENCT5uPY7xM4\n/BYC+f/BRgHHFp3LCSOmYTFYm19PlsmxunApbnTNkL6moiTzJxZLslrebAZFaW0o2TYgpfSInEzb\nZJseYSCaZSLeAAAgAElEQVT9hD3ZQFKka0m8HvRgcsvdrhpICnP5FkqeStaOBEeOZft51xEesW+z\nY+qi1Ty38T78jX2y0slvg4yUn0fEmouqNS8KtFpNRGIhIrofT8RLLBFPD+qRsMzmDVY2rrOxab2V\njettbGcFHHYrjPgXUqiIgVvm8QvTHxi1l8SI0SHyCuJIEpgMBtzWXKxyDlqieaJdkkgn5FPLhevi\n5dy6/Fq88QZyjO4+2xExRWrKyuuVycnRCIehqkrkZFqSbXqEgfQThIE0J1VLYoxGKd8eQ1VBVfW0\nkeh6J9pw6W3Xjmg2BwAPrr6aqkjT9rbFlhIu37dxGa8Ehhwn0ZxC4npT9JLqhVVTI/HUPw0k5CCy\nzcPxJwXIyW29VDllKl9v+pGvpIfw5H4KwUL48kr45mJybGZKR4cYMSrMiNEh9h4TZ9QQJzbZiaq2\nX9H+l+9nUhbahNEoo+sapbml/PPYZ3dJ8WJbpKasFEVCVZOLDKxWsiYnky1km57+bCC9s6ORYLdA\nM5rQCoowG6DQWgt+L7qqo6kSqposdE8kkq3eVTVpJpqW/Jd8XkfTJHyHHUtgwq8Y8OIjFHz4Mrlf\nf0TZzCvxTD4Kf8KbbjsiIeFPeJsE6JDw+DHF48h5RcRkc7PB+a9/NVFWJiNJOdjtOXyshZh7aR2+\niA+looLxj96HyedNt0jZZ1wp5/Iga33f8+rmx/nxyGsxH3E3edsuofbzS1j5UjFaY7Tjyo2z195R\nxu4LY/dRGD1aY8AAfcdSFzyxekAikUj2Aqvw1LNli9RYg5JMyJvNIMsd51F6goYGKa1PkiAclnj+\n+V0/ZSUQpBAGIkA3mUjkFSC7cpC9HhS/D1lOhiFmM7TsSQU07VuiS42G4iRx/Xxqp51AzgN3MHzh\nfEL/eZu9f2Pkf5Zk5biu6ziNOcjyDnta6aAGIxji5Sj5+UTMrvQlvN6mATMYhA1rrTgowGnLZfiT\n12GqLEMDjEE/B/7j/nSLlL1dB3DD/n9jnW8lS7Y9wf/k23CMfIjpxWczOnQuFT8XsXGdjY3rbbzw\njSVtKrm5OnvvrbL33hr77KNi0vbGb1mBoiS1u4x5xOMQjyfFN017yekVXgZDMreS3COsZ13F7dbx\neqX05+Z27zGTB4IsRUxh9QHZFmK31LPT1e077Duix2O8d/gQXjs0F7s9j0vHX0WeuSg9PZacNmuM\neHQJ1ZGDljeAQCjBTTeZ2bZNTkdDgwdr3HFHFICRF/we71Y/OhqSomIdZOWLB//Rppyf/T+yeOvj\n/LfhcxyGHKYOOpNjB52GzeAkFpXYstHK5p9tbN3gYtNaOxs3KumpLaPdh3XIWtzDtnLa5Mkcsn9O\nq0hlR5JNJpMmYjYnIxWzOWksRmOyJqW7f3GpZcQeTzIH0ldTVi3J9vu5r+nPU1jCQPqAbLvB29Mj\nJ+JIPl+3qtuhrX1HriXxi4kdvAisuS4Cjny2V5t4+GETfr+E06kze3acnBydWEwndPws8hs2okvJ\n7IonfxixxQvxRpsn3XfkZ/9PLNn6ON81/Ae74mTq4DM5dtAfsBua/sAVWcYiO6nclM+61WZWr1ZY\nv15h/XopbSoul05pqUZDAxgMUFyscc01MYqL2/5TWl1WwXWfXkdErsMu5fPQ1AWMGzYQk0lHUUj/\nk+XOGczucv/0FdmmRxhIP0EYSNt0pEfW1KSRdKG6fUcMy7/Cfu9dKGXbiR51DKFLL0fPb792xGIx\nEVU1kCSsC+5ArqpAc+cTvG0B2sBBSBJc9Yd6zvn+z+QkagmY83n+oPu48oECVE3HH/VTH/YQjkbR\ndL3VYoCNgVUs2foE/6/+U2yKg6mDzmDq4DObGYlBkXFanDiUHOxmJx5PjG++UXjiCRN1dRIej9S4\nRUvSVAwGnQMPVNlnHy09BTZ4cDJSOWHRLLzKRmRJQtN1crURvHXW00BTxJL6J8sSRmPSmGQZjMam\nabHkf3Ws1t3r/tnVZJseYSD9BGEgbdNZPZKuIfv94GmAeLzp8epqbAvvR/J60HNyCV12JXpRUfMX\nRyJYFz2DZdEz6GYz4QsuJjptenJkbEFqFZb9tpuQ1QQEAuAPoJaW4n/iWaDtpoepIrrUY4FYiIaw\nl2A0RDyhE48np8piMUgkdDYF1vDalif4pu7fWBUHUwf9gamDz8RhaMrDGBQZtyMHs25nwa25rF+f\nvOaWLRKyDAUFOtFo0kzy83U2bJBJJJKm4nTqjB6tslJbDMXfIxWtRsrZjklzseystzv9+2lpMg6H\nkUQigdGYMha9MXppimQUpSma6e2K9d31ft5VCAPpJwgDaZsub3Cla8gBP3i8EI9iv+4a5LKmPT+0\nkhKCt9/d5mvlrVuw33sXxv+3gsS+YwhefR3qPs1rR1IG4rjkAqRoFIqLQdPRE3G8S94F2m450tag\nKEkQUSP4Yl6CsQCqpiFJyZVlyRVmsK5+Hc+ufpwvKv+FVbFz7KDT+X3JmTgMueg6mM1GEvEEjz6S\ni+bPx1tnY9265JTWsGEaut6Un4nFYONGmbVrZdasUVi7Vmb1Wg20xk7BJj/Ggi3M+M1o9tknGamU\nlCQNoLOkPp+W7zP139Q/g0HCZIJ77jFRVSWnV9UNHqzxwANJs+2Jv/7d/X7ubfqzgYhVWIIuo0sy\nqjMHyeFCDgWQjAYwmZJRiSQheTztvlYbOgz/Q3/H9NFSbA89gGvOTKLTTyE89yJ0R/M/ND0nF6ls\nO5SVQX4B2qjRzZ/vxOCn62CWLRRaLLjNeXhjHvxRH7qkNa6Ugv0GjuK+gffyc8N6nv7pCV7f+hQf\nVLzEySNP4+QRZ5FjLSIUAoMlSH0wQu4QG+Nz3ZRttmG36+Tk6Fx2WXLAMpkgb/h2VkSuxTuonoFH\n5zHTdgX3vPcWoZpipOr9GRg5jMWLjcTjyVHfbtcZPTppJqnpryFDumYqO65qS5FI6EQi8N13Eps3\nSyhKsgV+JCKxZYuE0ShhMJCeIpNlvTFyaVpJlmpEued8zRR0BRGB9AHZ9g1pZ/U4L/sjSkMtBEPQ\n0IBWVEzwjrYjkB2R/H6sj/0d8+uvoefnE7rsSmJHHInFaiYSiVG/ugbPtQ9iCXmI2HLJu//P5P5y\nHzSnizlzmk9hDRi1DctJV9AQqcNtyc9YNZ7QE/jjPvwRL3GtdU5ng+dnnv7xSf619WMsBgun7fsH\nTh11BmogjyefNBEMShQUwDl/MDMoLwc9ZiUW09N1Mtf/L1mAKEnJ5b8ltlJuP/C55hoSyUhlzRqZ\ntWsV1qyR+flnmVgsaSo2W5OppCKVlKm0FYFk4sQTrXg8TZ9Vbq7GW2+1XT/SMpJRlKacTGp6LGU0\nqakyp9PYJT29Tbb9ffXnCEQYSB+QbTf4zuqRK8qx33QdciSMPnAQ4TPPQXM4Ov21VVn1E/Z77sSw\ndg3xSYeQuO4GwkUDuP56M2VlcjJPrUNJicbtd8SQXDnMunowm7c1BdBlh57OoLE/pwft0pxSnjj6\n2YzX1dHwx/14Ig0k1Hirqo1N3o08/eOTfLxlGWbFzPTRp3LGvmeTZ8lreu+ShN1kx2nMRVatRKMy\nZy49Dl/Ml377doOLR37xboefQyIBmzYlTWXNGpl16xTWrWttKmPHwl57xdOm0kYaqRlz51rYtCm5\nHFpRkvmixx/f+ZYnqTqYIUOMKEr/uZ97mv5sIGIKS7DTaAMH4X/8mfTPEjpyOBmNtNcBeEfUMWPx\nPfVcunbEcMo09JmzCDZcCJI5dVI8nuRXaN3rYS+bTrxgEJUNpuRqKGsdUuPXZ0mSqI/UdahbQsZl\nzMFpdBGKB/DFfEQSYbRGvaU5I7j1l3dy0YQ/8th3/+ClNc+zZN2rTBs1g7P2nUmeNR9N1/FHAwRi\nQcyKmVxrLnsVDGNV7U+ouoam6RQ58nA4JMLhZJuY9j4OgwFGjdIYNUrj+OOTjyUSsHlzMqeyenUy\nWlmyRCYSSXYetlqTprL33k2RytChzU2lsDC59DkVgRQW9sx3xlQdj2DPRUQgfUC2fUPqLT3NOgCH\nw52KSKSaGpyP/BXDsg+psg3n4bwbWOmYDHrzQsLqaoknn7FRIw8gJNnwHnUaVbGmaaPhOaU82UEE\n0pbeiBrBE2kgFA+h6clq/NSU0RbfZp758SmWbfkQo2xk2l4zOGvMTPKtBVSHqln47f0EYn4KbYXk\nWHLxRb0YZTM3H3o7Ax2D0HWJaBQiEYlQCGIxvVsDsMFgYu3aRDpSWbtWYf16mUgkaaBWq86oUU05\nlaIijTfeMOL1Sul8TVFRz/3ZDx8uIpBM9OcIRBhIH5BtN3hv65EkkMOhRiMJNWsl3xYWi4nEp59i\nuesuTJXb+cI9lbfHXM2sq12tBz5Zhlw3FaYo139+LfWdyIF0Rm9EjeCNeAjGg5jMhmZz/Ft9W3jm\np6dYtvlDFNnASXudTJl/O7XhmnTTyKHOoSz4zX3kWHJxGBwokqHVNRKJpKGEwxLBYPNGlh19Pi1z\nDokEbN2aMhS5sQCyyVQsFp299mrKqey9t8bw4RqGHpiDEAaSme4ayMUXX0xlZSWxWIyZM2dyyimn\nMH78eE499VS++OILCgsLeeCBB3C73axevZqbb76ZSCTC0KFDufPOO3E6nWzbto1bbrmFhoYGrFYr\nt912G6WlpT323oSB9AHZdoPvKj1ttZJvi/QAGY0ma0eee7rD2hHJbkcrKELriRFxB70RNUJMClHn\nb0BtoXebfxvP/vQUH256H03XyDXnkmfJxygbsRvtPPK7xwBQZAWX2YXLlINBMrZ7rWhUIhKRCIfJ\nON3V2SS6qsLWrRJr1ijpSGXdOplwOGkqZnMyUmma/uqeqQgDyUx3DcTn8+FyuYhGo8yYMYNFixZx\nyCGHcP/993Pcccfxt7/9jYaGBv7yl79wwgkncOONN3LwwQfz0EMPEQwGmT9/Pueeey633norQ4cO\nZeXKldx///08+2zXIvNMCAPpA7LtBu8LPXIsmjSSQKCVkbQcIOWtW7DfdzfGb5a3WzsCJNedFhSg\nNraS7ylsNhPeQABPrIFA1J/OkaQoC2xn3icXUREsR0Iix5zDmLxx3P/bhc2OU2QZh8mJ0+TColgy\nzuilprvC4eR0VzyuUxEsZ+Gaa/GrDTgVN5ftcxdFlq5FWaoK27a1NpVQKGkqJlNrUyktzWwqwkAy\n010Defjhh/n4448BKC8v58knn+SMM87ghx9+QJZltm3bxqWXXsqiRYs44YQT+Ne//gXAtm3buOyy\ny3j++eeZPHkyI0aMIDXMJxIJ3n234wUdnUUk0QV9gmYyoxUWI+e6kT0e9IC/3YhEGzoM/8K/JWtH\nFmaoHYnHobISQ24uqjsPXepCIUUjckU59hvmIzfUpdunMHI4RtlEoaUYp8lFQ7iecCKc/qMc7Cjh\nH0c+xd3Lb2ddw1rqInWsqPqae79ZwMwxsyi2DwBA1TS8ES/+qB+r0YrL5MJmaNvsJCm55a7FopOX\nB/G4xC3LrqUyuglJkvFFPTy05tpWy4NTVFdLLFxoapX3UBQYPlxn+PAExxzT+PlqyUhl7VolXQD5\n4YdGXn+9yVSS019NxjJiRM9MfwnaZsWKFXz99de89tprmEwmzj77bKLRaKvjUgtH2ooDNE3D5XLx\nxhtv9JpOcQsI+hTNaEIrLELOzW0ykraQJGJHHUN88i+TtSOLX8X070/StSPpAgZdR29oQA6H0QsK\n0cyWts/XDvYb5qNs3gSShOL1Yb9xPuoLLwFQESjnhi/mE9eiDHAM5Kx9z8FudKDrOkW2Iu7/7UPp\n455d9TRvb3iTtze8ye9HnMjMsbMYaB+YfM+6RjAWJBQLYTSYyDHn4DQ6kWjb8HQ92a7EF6/DaJQa\n6y8kQlI9TmdydVdqv5YUCxem9lKBQEDioYdM3H576wEIkmmklKkcfXTj70VrilRSprJ0qYHXX09O\nwZlMOiNHJk3l8st19tmnSx+zoAP8fj8ulwuTycSGDRv4/vvvAVBVlQ8//JCpU6fyzjvvcNBBB+Fw\nOMjJyeHbb79lwoQJvPXWW0ycOBGHw0FJSQkffvghxzR+W1izZg379OAvS0xh9QHZFmJnkx45FsUW\nDRKubWg3IoHG2pF7F2BYs5r4pEMIXnkt2pAhLU6WTLBruW709vqvtyBn+u+RfE33ie5yEv9gGaFQ\njDPfPZUVlV+R0FUMksKUob/jr0c8gifiaRaRpKgIVrDop6d5Z+NbABxXejwzx85ikGNwq+saZANO\nsxOH0YlZMbc5vXXe0nPY7N2EosioqkaxqZTcj14hHJYoLNQ577w4VqtOPK5zySUW/P6m92y36zzy\nyM7VfmgalJXtOP2VrFW5/HKVSy7Jno2tsul+hu5NYcViMS6++GLKy8spLS3F7/dz8cUXc8EFF3D6\n6afz+eefk5+fz4MPPojb7WbNmjXcdNNNRCIRhgwZwoIFC3A6nZSVlXHTTTdRU1ODqqpMnTqVP/7x\njz323oSB9AHZdoNno56Ix99ujgQaGzj+9T6Un9ehlJeDLBM+ZzaRs85J7YLVhMXS6WjEed456QhE\nikWRvB4YNZp4jpuDx/wfG+zh9FJhh9HBqtkbkSQIJYLtGklVsJJFq57l7Q1voOkaU0f8nnPHzmnT\nSGRJTk9vWQ329G6OkIxsbvxiPp7GPdobXnqA7z8fiqZJyLLOr36l8uyzEWIxib/8xczmzTKRSLIU\nJz8/uX9IbyByIJnpyWW848eP57///W+PnW9nEQbSB2TbDZ4telL5B6O3gXiOm+Add4E7r00jsV+/\nQwPHeBw54EepqkIdMjS578jESS1O3rloRK4ox37jfOT6OpQNP6O5cpGtFjRV4z19FWec1rQCzG6w\ns2r2xvTPKSNpiDQQiUfQW9S2V4eqkkby8xuousoxpcdx7tjZlDhbRE4kt/81GIzkmHNaLQNO/b7G\njLETDDaPMlatCgLJZpO3325ubLuiM29eDLsdQiG9cYveDn4ZXUAYSGZ60kAOOuggvvvuux47384i\nDKQPyLYbPFv0pL79y4qMpmrN2rfL0ebLfx2XXIAUDKZfq9vtRGbOSu47sn0b0SOPJnTZFa33HelC\nNJKazpJlCU3T+TG2hcPPVtF0FVlSmDRwMi8c92obr9QJxgPURxuIJVp/668OVfPCqmd5c8MbqFqC\nY0qncs7YOQxpw0gguQzYbrInV2/JFmw2c9pAAgEpXWHucDQZSHukVneFQsnak0Sic7UnmRAGkpn+\nXEgokuiCrEFuqGvWzU+ub2pHopktaEUDkkbiaUAvKkbauCHdQl7PySUx8RC8z7+Srh0xfvVF69qR\nSASpvAzFnZeMRjLo0dz5KF4fkLxGaelEfjvE2qxYsW0k7EYnNqMDX9yLJ9xAokXTxtpwLfvk7k1D\n1MNHW5bxwab3OGrYscwaN4ehrmHNjlU1FV/Ehz/qx6iYKJLzUHQLkyerfPaZQiKRbIMyfnzHZe07\nru7Kz4dYrMlQ2krGCwSZEBFIH5Bt35CyRU+rCGR4Kf4nWxc9SRIo27Zgu+dOpM2b0W32VptYdaZ2\nRLJa0QoK0UzmlpdInqNxOsvoaZxSa9wRsask9ASeaD3+qB9N17j+P9dQFtierlovtBYyxDWM19e/\nRlyLc+Swo5k17jyGuYa3eT6LxUQ8liDkdfDs4wVUbnVgMunceuvObRTV3VYrIgLJTH+OQISB9AHZ\ndoNni56WA3bokj9he+SvzWoydhzAJQnkYADqGyAWbf3VWdfT+45IDQ1t1o5IioyeV4DmcqHTdm6k\npz6fmBbFE2lg7tJZ1Efq04+nqtbrw3W82NiwMapG+d2wo5g97nyG5zRvPbFjoaUkSVgMFlxmF3aD\no92lwF0lmVpqHZ20Nd0lDCQzwkD6CcJA2iZb9ey4Igpdb5YT2REJPWkkDR6IRWg5LyUF/Fj/0Xrf\nkabpMpCstmQ0YjS1q6cnkCS45rPLWd+wHk/YQ0yLMdhRwh2/bto/pT5Sz0urn2fJ+leJJCIcMfRI\nZo07D4fJycJv78ef8OE0uLhswpUU2ZqirtRSYLvRkbHSPVXP0pm9U3YkFpPS0Ukk0tRqRRhIZoSB\n7GIeeeQRXn31VfLz8wG4/PLLOeywwwB47LHHWLJkCYqicP311/OrX/2q0+cVBtI22aqnrZqM1Ja2\nbSEBcsAH9fXN9mxPoaxehf3uO9L7jgSvvAZtyNCm1xsM6O68VtFIT38+FYFy7l5xOzEtig6cM2YO\nhbbCVsc1RBrSLeTDiTCF1kJsRhtWow1N0yhxlHD7r1tv3CVLcjoqsRnsraKSVD1JV/ZO2ZGWjSDz\n8gxIUvbdP9lCfzaQrE2iz5o1i1mzZjV7bMOGDXzwwQe8//77VFZWMmvWLJYtW5Yu5xf0L9JJ7NRe\n6+78do9NtyDxNqANGUbk4kvRDQb0RFPyWt13TLN9R3LOOo3wzFnp2hE9kYDaapRQMBmNGNpufLiz\nDHQM4q9T/g6Ajo4/7qMhXN8q0e62uPnjgfM4c9+zeWnNCyxa9QxaWMNpcpJvLsATbXvrYE3XCMVD\nhOIhjLKBSMDFEw/nUVVmw2TSqTyqvst7p+yIroOi6NhsyU2ubDYIhbrxQQh2e3pmwrQXaCsw+uST\nT5g6dSoGg4GSkhKGDRvGypUr+0CdYFcQvG0BamkpusuJOrw02ZeqHVItSCSPF2Xl91juuysZXeTl\nI+3YtElRiJ5yOp6XlxA77HBsTz5GztmnY1ixPPm8DnowiFy2HcXvRcq4TmvnkZBwGXMocQ7BbXUj\nt9G/K8ecy4UHXMwvB/6afEs+wXiQzf5NbPVvYV3D2oznj2sJbrwjyPebygkbtxMjQKhiBIbGuhJd\n13Fb2jdmQf/G6/Vy8cUXM378eKZMmdLlRosdGsjzzz+P1+vttsDu8sILL3DiiSdy/fXX4/cnpzGq\nqqoYOHBg+pji4mKqqqp2uTbBrkEbOAj/E8/iXfIu/iefzbgCqtUS4OoqNFlBdeehDR6ClJubLCZs\nRC8sJHjbAnwL/wa6juuyP2K/8Tqkutrk84kE1FSjVFYgRTtuAVIRKOe8pecw/a3fc97Sc6gIlHfp\nvSqSgTxzASXOIThMjjaj6qsmXsuBheM5oOhAhrtKCcVDnPPBGVzz2RWsrV/T7rm9XolIVKeqPsx2\nTyXuDbMZVziOEbmljCkYm2E5sqC/c8stt2A2m/nqq6+49957ufnmm9mwYUOnX9/hFFZNTQ0zZsxg\nzJgxTJ8+nV//+tc9MmU0a9YsamtrWz1++eWXc8YZZ3DxxRcjSRIPPvggd911F3fccUebUUl7WpYv\nX86KFSvSP8+bNw+brXWCtC8wGpWs0QL9Q49cUIDk9zfVhRQU7HAOE7hsSLECpPp69ECwacXWYYcR\nnTQJ49NPYfrnk5i++pL4JfNIzDg1WTuixjFUVmB3uNDz3NBOh9+bP76Orf7NSJKE3+/nluXX89zx\nLwAglZdhuubPyWvn5RG7+170Qa3bmKS05jgcBOJ+6sN1RBNNuZyhlhLuO/JBDAaZRELDH/Xx4qoX\neP6nRXz24ZkcNuRwLhh/IWMKxjY7o9stEQwmCw4TKpj1fP406QrkxhVcTrMLi1FBljrYXL0d+sP9\n05u0HIsmTpzIpEmTMrxi1xAOh1m2bBnvv/8+FouFCRMmMGXKFN566y2uuOKKTp2jU0l0Xdf5/PPP\nef311/nxxx859thjmTFjBkOHDu3opTtNWVkZF154Ie+88w6PP/44AHPnzgVgzpw5XHrppRxwwAGd\nOpdIordNf9CzYwuSlkt+m7VoHzSY8NXz0RUDRMLNVmy1VzuSXjZrNkN+PprN3mqF0/S3fo8v1nR/\nuUxOlpyYnA7o7GqylrSXH2m5X0og5ufVtS/z8toX8cd8HDroV8zZby5j8pNGUl2d7Mbr8bS/pa1R\nNuDoxAqutugP909v0t0k+ttvw2uvJQPnP/4RetpzVq9ezR/+8Af+97//pR/75z//yTfffMOjjz7a\nqXN0KgciSRKFhYUUFBSgKAper5dLL72Ue+65p3vKO6Cmpib9/x999BGjR48GYMqUKbz//vvEYjG2\nbdvG1q1b2X///XtFg2D3ItN0Vzo/4vOjrF6N9dab0AaVQFFxchOq1Dka9x0J3HoncnUVrjkzsT1w\nDzROoRKNQkUFSnUVcqJlwjs/HSG3zCtkqrDPxI75kVxrbpv5EQCHycns/c7njRPf4YL9/8iPtT8w\nZ+lMrvj3pfxU+wNFRTq33x7lkUci3HFHtM390ONagoZwA+X+7Wzzb8MX95LQE4j1KX3DF1/AbbfB\n6tXw449wxRWwZUvPXiMYDOJwNN+Pxul0EgxmboezIx1OYT333HO8+eabuN1uZsyYwdVXX43RaETT\nNI466iiuvvrqrivvgHvvvZfVq1cjyzKDBw/m1ltvBWCvvfbi2GOP5bjjjsNgMHDTTTeJFViCDmlr\nANcB1eFCsjuQPR4krwddVZP7jhx5NPHJhyZrRxa/iv7vf6FddkWydgTQfT7kUAgpJxctJxddkrjt\nlwu48Yv5bbY56cpqsrZQJAP55kKcJhf14To0Em0eZzc6OHfcHE7Z+3SWrHuVF9cs4rxl5zJp4GTm\njJvLfoUdf9nSdJ1oIkJNIoIiy1gMVpwmZ5vLgQW9x7vvJos2JSn5z++Hf/8bzj23565ht9tbmUUg\nEMBut3f6HB1OYS1cuJAZM2YweHDrOdsNGzYwcuTITl+srxFTWG3T3/W0mkJqo0WKnEggN9S32hlR\nWb0Kx70LUFavanvfEZMZ8vPQ7I52p30yTa91FUkC1RClvKGaaCJzcj8UD/H6+td4YfVzeKIeJg6Y\nxOz95nJA4YFdvq5BNuAwO7AbHVgVa7P32t/vn52lO1NYjz4K//hHU4CcSMADD8CUKT2nKxwOM3Hi\nRN577710OuKaa66huLi4Z3Mg/QVhIG3T3/V0ZQCXo5FkhBIOpxPtFqOC/uIL2P7xd0jECc+cRfTo\nqRvzb/4AACAASURBVNgefTgZuRQUEr7metRRo7u8A2J3sNlMBEPRdutHWpI0ksW8uPo5GqINHFw8\nkTn7zeXAovFdvrYsSRgVc2ORog2DZOz398/O0h0Dicfhwgth5cpkDuS3v4UFC+jxKcUrr7wSgNtv\nv51Vq1Zx4YUX8vLLL3c6MBAG0gdk2w0u9DQn3RqlscdWKmkt1dRge+gBzB8vQ7Pa0PLz0R2O5LRU\nSQnBO+9FcrnQct1oSu/V6O74+ah6Am/MgzfiRdMz92UPJ8K8sX4xz69+joZIPROKD2b2uLkcVDyh\nWzoUWcZmsFOQk4cUNzbb/Kov6ev7pyXdTaLrOpSXJ6OQHfqE9iher5frrruOL7/8ErfbzVVXXcXU\nqVM7/XphIH1Att3gQk/bSLqG7PFgiQSJhJqmiwzLv8J5zZVI0Sia04VaXIRutaMXFiYjksFDCM//\nC4kRI9HbSXzvDG19PjE9Sn24jlAs1Gojq5ZEEmHe+Pl1Xlj1LHWROsYXTWDOfudzUNHB3copWiwm\n1JiG05Lc/MqkmPq0JXy23D8p+nMrE+Xmm2++ua9F7Cqy5aYyGhXi8U70yd5FCD3tIEnoVismtws1\nGkv319JKhqCsWY0U8CN7vcie/9/encdHVZ+LH/+cM0symUwSyL4BAi6A27VK2ktLNVBBxCoCtq9r\nBVHR9oeoaF3AqrggVr3WhV6VRVzqrQuLVkCrgIhWAbFVrqBVEdmyk5A9M5k5398fkwxkG5LJJGeY\nPO+/yCQz82TmME++53yf5znsX7Eoheb1opWVYv3ndrzn5KFbdLDZw3ruob3Xx6JZcdld2C023IY7\n6GrEqts4LeV0Lj1xKv1i+vHRwQ9Y8c1rfFb8KRnODDKdWV1KJFarBU9jI/WN9VR7qmnwNaBpYNPN\nWZVEzPHTxOlsf1xANJAEYoJIO8AlnuBssTF4YuPQ7DZ/Vbph4D39TCx792A44tAbG9GrKtFqa1Gx\nsf5zDj4fnvEToLYWrbYW3aqD3Q5h+EAN9vrYLTG4YlxoGnh8nqCrEatu5dSU07j0xKkkxybz4cHN\nrPz2NbYXbyPNmUGWM7tTicRqteD1+uNRKBp9jdR4aqj2VGPgQ9d1bHrvtd2LtOMnmhOInMIyQaQt\nsSWe4I6OR/d50cvLUdVVR6rZlcJ1zZVY//0V+HwYSUl4Tz+T2j/+95EH0TR/IWK/fkF3bHU1nmDc\nRgPl9Yeob6w/5mktALfPzVu73+SlXc9TUlfMaSlncNWp15CX+ZOgiaR1YWNruqYR0wMzSzoSaceP\nnMKKEpFyUEXaX0gST3BHx6N0HRXvRLPb/IWFTZv1G38yCv3gfrS6OvTyQ+hVlRhp6fgGDzly+srr\nhdpa9IZ6/yrFFlq3386+Ptam01o2iw2PL/hpLfCvSIYnj2DyiVNJjUvl44KPWPnt62wp/JjUuHRy\n4nPaTSRHr0DaowCv4aXWU9srq5JIO35kBRIlZAXSPoknuI7iaXc1QtPckYcfxPr1VzSOzKP293e0\nmDviv7OGFuf0FyI6HD3eOsSnvBz2VFDVUHXMRNLM4/Ow9vu3eHHXcopqCxmePIKrT72Wn2SNapFI\njrUCac/Rq5I4qxOd0PpwtSfSjp9oXoFIAjFBpB3gEk9wweLRNNDr6+DQIWg4qrDP5yNm9UoczyxC\na2xsMXekBV3zT0NM6tfpRNLV16d5AmFdYw05CQOYNvxK4qzOTp3WAmj0NbJuz1u8sHM5hbUFDOs/\nnBmnzuSn2f7GqqEkkKMdXaQYq8fS3etEkXb8SAKJEpJA2ifxBNeZeDSl0Kur0CrKWwyx0spKiXvy\nT8S893d8uQOo/f0deEe20xVP19Cc8Rj9+rc7Vrer8Rzt8rWXsbXwEwzlQ9csjBnwCx477wkONRyi\n0dd2cmNHvEYj6/as5YUvn6Og9iAn9zuFq06byS+G/AK3u/OP0xFN07C3KlIMRaQdP5JAooQkkPZJ\nPMF1JR7d60UvL0PV1LQ4rWXdtgXnIw9hObAf9y/GUXfjzajklLYPYNHRXAkYCYkdJpKuvj7DnxtM\nrfdIzyOn1cmuq75HYXDYc5jKhgp8RudOa4E/kbyz522e37mMgzUHOLn/KVw54mpG55zbYcPHrmou\nUnTFuHBYOt+bCSLv+JEEEiUkgbRP4gmuq/FoGuh1tVB2CDzuI99wu3G89DyxLy5HxcRQf90sPKN+\nRtyix/0FiIlJ1N14CyotDXQdLT4eIyERFRvbrd5Tw58bTE1jTWAGerwtnl1XfR/4vlc1UuEup8Zd\njdGFjwOv4eXvP7zNi7uWs69qL0OTTuSqU2fy89zzwpZINE3DarGRYE8g3haPTbcd8zRfpB0/kkCi\nhCSQ9kk8wXUmnhYzR5r6bamMDCwV5ajKyhYNGo+eO2K4EjD690M54o60RHngj0c9sI7mcKCSkjAc\ncSgVwimsNZexregTvMqHVbMwMuMnvDzxtTY/12DUU15/iIbG+i4N8rXadd7691s8/+Uy9lXvZUjS\nUGaceg3n5Y4JWyIB0DUdh625O3B8h0WKkXb8RHICefnll1m1ahXffPMNEydOZOHCrk2nlG28Joi0\nbYYST3CdiSd+zvX+mSMeD/rhCqyff4b7l5diOOLQ4uLQ3B7w+a+NqMQkPOMn4Bs4CPuGd9ErKvwt\n5evqUPYYPBddfOSBlfJXwNfUoNXXoVt0rI5YGr2dP+X0n1mj2FP5PUkxSZySPIKHRj+Ky972Q82q\n2Y5s+zU8nd6tZbfZOME1hEtPnMIA1wA+K9nO6m9Xsmn/RhJjEhmYMCgsiaS5SLHWU0ONpwqfZmDR\nLFhbbQeOtOMnkrfxlpSUMHLkSOLj4/F6vYwdO7ZL95cViAki7S8kiSe4zsSTOHkiWtWR40sluKhc\nuSbwtaYMLOWHUFWVYBz5L5d40Tj08vLACkUlJHD4nY0dtz7RNGIT4mhwuPwFiT3UKsTAx2HPYaoa\nDh/z+kjrXVg+w8eGfe+x/Mul/FC1h0EJJzDj1GsYM+AXWPTwbdcF/6rEP5bXFShSjLTjJ5JXIM0e\nf/xxiouLu7wCkQkxQoSB0S+5RWV666FRStPxJqdCRmaLAkJfZpZ/RG5MDOg6elUVrhtnoe/f1/4T\nKQXuRiguxlJYgN5Q3yNTA3Us9Lcnkx2fS7w9vku9sSy6hfMHjecvE17l/lEL0TWdez6+k8vXXcbf\n96zDZ4RvdWAog7rGOopritlXtZdD7tJjtrc/Xvzt33/jilVXMH31dLYe3Gp2OO2SU1gmiLQltsQT\nXGfi8Z6Th/Xzz9AMH0Zmlv8aiKvtX57KZkdzudANAzwebFs+9t/er79/C29yCpbvvyN2xatoPh/e\nEaeBteUpmkDld2MjWm01utvdVNke/spui2bBZXcRa43BYzTia+fDuaNKdF3TGZw0hEknTmZw0hC+\nKPkXq79byfq97+KyuxiUeEJYr5EYysDtc5PoSAAjcv42DuUU1j/2/YN5G+ZRWldKSW0JG/ZsYOwJ\nY0mKTeqBCGHLli3U1tZ2+RRW73U4EyKKNc9k75TiYhwP3ovu84I9BiMlFa2hPrALC10j7sk/4Vi2\nGPvf36b21jvwjvxxB0+sUDU1aHW16HFO/yySMA+1UgocFifZ8XGdHmJ1NF3TGTPgF5yXO4YP9r/P\nc18u4d5P7uK5L5cwfcRVjBt0QZvrGH3dmm/XYCgDTdPQNI1qdzXv//A+V555pdmhtRA5aVqIPsJ5\n11wsX32F9v0etEYP6oRB1Px5MbUL/ohKS0OlpFJ734NUPfFnABJunIXz7nloZaUdP2hzIik4iLW0\nGL0x/KttDY0EWyI5rlySHEldXj3oms55A8bwwgX/y8KfPUKs1cEDW+bz6zWTWbP7b3iN7hcjRosB\nCQNaJGld0xmQOCDIPcwhCUSIXqZXHGoabu6DsjL0sjJIS0NrdarKO/LHVP7lVequuQ77po0k/noy\nMSte89+vI4aBqqpCO7Af66FS/yonzCyaleSYVLJc2Thsji7fX9d0zs3N54XxL/PH0f+N0+ZkwdZ7\n+fWayfxt9xuSSIBrzrqGs7P8A74suoXxQ8dz3qDzwv48Pp8Pt9uNYRj4fD48Hg++YMdXK7ILywSR\ntktE4gku3PG4rpmO5Yc9/iSiFL5BJ1C99AV/FfuhUlRtLa2r5fT9+3A++hC2bVvxDR9Bza1z8Z0y\nLPB9raSEuCf+u01Boma1ohKTUAkJGGHeAdXMo9dSeLgk5IvXSik+Ovghz325mK/LvyLTmcX0ETOY\ncMJF2Cxda2eiaRpDUk5ANfbM7xqK0EfaKgqqC7BZbKQ5e2am7aJFi1i0aFGLTRKzZs3i+uuv79T9\nJYGYINo/ILsr2uPRCwtw3j0XvfxI0aGRmQU0VbFXVcKhMvC12j6rFPb17+J88jEoL8d96VTqr/sd\nKt6F887b0Q8eCCSl1gWJmtWKciX4E4k1tB5THYmLs1NVW8uh+jJqG2sJ9SNFKcXHBR+x7P+W8FX5\nTjKcmUwbPoOJg3/Z6UQSTQnkeCAJxATR/gHZXRIP6O4GtNJScDe0+V5soxv9ySeIWfkaqn9/6m68\nBfsbK9Hq6gI/o5xOahY92/aBLTpavAsjMemYTRs7q/n10TSo8VR3uUlja0opthR+zLL/W8zOQ1+S\nHpfOtBEzmDj4YuyW4DFLAuldkkBMIB+QwUk8fpoysJSVtZk30ly4d/TcEaNff3xJSf528UphZOdQ\nu+CPHT94kF5bXdX69Qm1SWNrSim2Fn7Csi+X8GXZDtLi0pk2/EouGnJJh4lEEkjvkgRiAvmADE7i\nOaL5lJZ2qAzVdEqrReV309yRuKefgoYGfOkZ+E46mbqbb/M3ZTwWXUOLParXVgiV7R29Po2Gh/KG\nQ906rQX+RPJp0VaWfbmYHaVfkOpI44rh0/nl0EnEWFrWWEgC6V1SSGiC47FQrjdJPC2pmFhwHOmn\n1aJwT9fxDR+B+8KL0EtLsH/xeVMB4giM7JxOPDhHem3V1qHrGtjsHbdSaUdHr09zEaLdYsPdiZG6\nHdE0jWxXDhMH/5Iz0v6D3Ye/ZfV3K1nz/ZtYdStDkk4M1JFomkb/uH7HfSHh8UJWICaQv7CDk3ja\npykDy6FDxHjqaah3t/sz1m1bcD76Ryz79/nnjtwwB5WS2rXnsVlRrkSUKwHDeuwCv868PgY+Ktzl\nXRqp2xGlFP8s2c6y/1vCv0o+Izk2mcuHT2fS0Etx2OJkBdKLJIGYIFI+kJpJPMFFUjyaBnGNDbgP\nFraYfNiC203sX17A8cJzKLud+utm4b50Cli6+KFq0f1TEl0JQcftduX1cRsNlNcfor6xvtMjdYP5\nZ/FnPPflYj4r3k7/2GQuHz6N239yO3YV1+3HDpdoTiByCssEZp8SaU3iCS7i4omPo9EWg+ZxQ3tJ\nxGrFe9aP8Iw9H+t33xC74jVsH3+E9+RTUKldWI0oBW43Wk21v5V8B6e3uvL6WDVr02ktO24j9NNa\nzTLjs5gw+CLOTh/JD5V7eOO7lWS5sjgj5axuPW44ySmsKCErkPZJPMFFajwd7dJqQSnsG94j7vFH\n0crLcU++LFA70mUaYIuBpESMeBeqqZVJqK+PgY/D7goqGyq7nUia7aveS17uSOJICMvjhUM0r0Ai\n50qTEKIFvbAA1zXTSZw8Edc109ELC1p8X2k6vrQ0SE0FSwf/lTUNz9jzqXxlJe4plxGz8jUSfz0Z\n+3t/7zjpdEThH9FbUoJl/14shyu61SpFx0L/mJRAS5RwdKUfmDAIp61rM9RF6CSBCBGhnHfN9U85\nrKrG8sMenHfPbfMzSoHPlYjKzAZ7x6dKVLyLuptvo2rZixipacTfPS/43JFjUI1eOFSGfmA/+qEy\n9I6ux3RCjB5LljObVGeadOXtRR6PhzvvvJP8/Hx+9KMfcemll7J58+YuPYYkECEiVKDpIoCmoZcf\n6vBnjZhYVHY2WryLYH/K+4YNp2rpC9TefBuWXV+S+JtfEbtsMbjb39V1LMrrRZVXoB/Yh7W8DN0b\nagW6hsuWSLYrl4TYxC4NsBKh8fl8ZGZm8vLLL/PZZ59xww03cNNNN1FQUHDsOzeRBCJEhDrWlMM2\nP69b8KWnQ3Iq6EH+a1ssuKf+ispXVuL5+XnELX2WxCt+jXXblpBjVT4fqqIC7cA+fxfgENvJWzUr\naY40MuMzibFG78XnSOBwOLj++uvJzMwE4NxzzyUnJ4edO3d2+jEkgQgRoWrvX4jvhBNQCS58g06g\n9v5jz6tWaPgSk1CZWS1G57b7s6HMHTkWn4E6fNjfTr60GN3T9ZXNkQFWuSTHpWAJlgyj2d/+Bldc\nAdOnw9aeH2lbVlbG3r17GTp0aKfvI7uwTBCpu3oihcQTXGfj0X1e9NISf3v4Y2muHXlxOcpm61Lt\nSIvWKm2C0NGcTozEpJB7bjUaHg41lFHXWHfMlihR08rkH/+Am27yZ1OlIDYW/vd/YeDA8AcIeL1e\nZs6cycCBA+lKZUcfTe1CRD/DYsWXkQn9k4Of0gKIiaHh6mup/Mtr+EacivOxh0m4ZjqWr3Z1MwgD\nVV2NVnAQS1EhekN9V7qkAGDT7WQ6s0h3pnd5Pshxa80aaG5EqWlQXQ3vv98jT6WU4tZbb8Vut3PX\nXXd16b6SQISIYgoNX7/+kJFxzFNaAEZuLtWP/5ma+xeil5aQcPU04h79I1pNN1fvhhEYuWspOIil\nrgatC5XoSoHT6iLXNYD+ccnRf1prwICWRaK67r+tB8ybN4+KigqeeuopLF3sVhDl74IQAsDncKKy\nctDiOtHio3XtyOoVodeOtKYUqq4OiorQDxzAUl2J1oUiQg2dfvb+5LgGkBCb2OW57MeNa66Bs8/2\nrz4sFhg/Hs4L/0jbu+++mz179vD0009jt3d9PoxcAzHB8XpOvbdIPMF1Jx5NKSwVh1CHD3c6GVi+\n2oXzkYVYv9pF4zl51N56B0bukb+Gg14D6UxMgWmJiZ1q3ni01r21ouYaCPjfn4IC/8qxM635u6ig\noID8/HxiYmLQm1Z0mqZx3333MXHixE49hiQQE0TTB1JPkHiCC0c8ltpqtLKyjhsyttY0d8TxzCK0\nxkbqp82g4TfTISam2wmkmWbRwZWAkZDYxWmJilpvDeUNFXh9HgZHSwI5Dpi2/nvnnXeYOHEiw4YN\na7Pv+Nlnn+X888/nggsu4KOPPgrcvnnzZsaPH8+4ceNYvHhxb4csRNTwOV0YWdn+3T2dYbHgnnJZ\ny9qR3/yqW7UjramjtwCXlXRhC7CG0+oiJz6Hfo7gtTIivExLICeddBKLFi3inHPOaXH77t27efvt\nt1m3bh1Llizh3nvvRSmFYRjcf//9LFu2jDVr1rB27Vp2795tUvRCHP8Mmx0jKxstMbHTA6Taqx2x\n33Fb92pH2gRmoCor0Q4e8NeStDMXvj0aOkn2fjhsjvDFIoIyLYEMHjyYQYMGtdnXvWHDBiZMmIDV\naiUnJ4eBAweyY8cOduzYwcCBA8nOzsZms3HhhReyYcMGk6IXIjooTceXeoyGjO3wjvwxlX95lbpr\nrsOycT2Jv55MzOuvgi+Mbe8NA1VVhVZwEGtJEbq7octbgEXPirgtDMXFxYHSeoD09HSKi4vbvb2k\npMSMEIWIKp1tyNhGU+1Iw+urw1s70lqgluSAfwtwfW2XtgCLntOjrS9nzJhBWVlZm9vnzJlDfn5+\nu/dp75q+pmkYRnjmBQgh2mfExKJnZ6OXlqJqq+nsZ7QaOJDqx/8cmDuScPU03JdOpf63/y+0uSMd\nBti0Bbi+Ht3eNJfE6ULJssQ0PZpAli9f3uX7ZGRkUFhYGPi6qKiItLQ0lFItukQWFxeTFmRr29at\nW9m2bVvg69mzZxMX1/V9zj3BZrNETCwg8RxLn4vHmYt2uALKK45UQwdhterEOmJg4kQafv5zbP+z\niJjXXiHmg414br4V3/gLOn2NpfMMqKyAuhpITICERFRTEVykvV+tP4tGjhxJXl6eiRGFT0Q03z96\n1ZGfn8/vf/97rrzySoqLi9m3bx+nn346hmGwb98+Dh48SGpqKmvXruWxxx7r8DHz8vLavEmRshUz\nGreFhpPEE1yvxBMTj6WfDiUl0Bi8RXuLbby2GBpuvAXLuAk4H36QmHm307h6VZvakbBp8EB1LZq1\nFFwujIRESHBG1PvV3mdRtDCtDmT9+vXcf//9VFRUkJCQwCmnnMLSpUsB/zbeFStWYLVaufPOO/np\nT38K+LfxLliwAKUUU6ZM4dprr+3Sc0odSPsknuD6cjy614teFrwhY4d1IEfXjng8NEybQf0VV0JM\nD7Zpt+g4kvtRZ3OgHI5uF86HQzTXgUghoQn68gdSZ0g8wfV2PBr4h1lVVoDR9uPiWIWE2qEy4p78\nEzHvvoMvJ5fa39+ON+8nPRZvbKydBk8jmiMOlZSE4YgzNZFEcwKJuF1YQojIogBf/2RIz0Czdf2s\nt0pOofbeBVQ98T+gaSTcdD3Ou+aGt3akNUP5V00FBegHD/gr73vu2fosSSBCiE7xxcVjZOV2riFj\nO7wj86h86RXqrrkO++ZNPVM70ppSUF8PxcXo+/d1uXmjCE4SiBCi0wyrFV9mNvTrD3oIf9MH5o68\n2rO1I60pBR43lJRg2b8XS0V5N+a3i2aSQIQQXRI4pZWRGdIpLQAjd0DPzB3pBNXohfJD6N2c3y4k\ngQghQuRzODGyctDiQuw91dHckXff6f7ckU5o0bxRWqWERBKIECJkhtWGkZWFltQv5GJBFe+i7ubb\nqFr2IkZaOvH33Inrxlno+/aGOdoOtB67K4mk0ySBCCG6R9PxJqdAejpaFwdCHc13yjCqljxP7S23\nY9n1JYm/+RWOpc+Cu7Nt3bspMHb3AJbCAiwNdZJIjkESiBAiLHxOF0ZmFsR0csZIe5rnjry6Cs95\nY3AsW+yfO7L1k/AFeixHbQEONG+URNIuSSBCiLAx7DEY2dloCQnd6n/Vbu3I3fN6tnakTRBNzRsL\nC/21JHU10gW4FUkgQoiwCswYSU4BvXsfMS1qRz54v6l25JWerR1prbmWpKhIihJbkQQihAg7hYYv\nMQkyQqtebyFQO/IavlNPw/nYI71TO9KaUtDQAMVF/qLEmiq0vtMJql2SQIQQPSaw1dfR/TGzRm4u\n1X9a5K8dKSvt1dqRFhT+osTiYvQD+5uq2/tmIpEEIoToUYbVhi8zGy0xqftzQZprR/66wpTakTaa\nqtv1/XuxVPW9NimSQIQQPU5pGr7UVEhLQ+vC7PUOH6+j2pH9+8IQbQgaG6G0qU1K5WE0oxev0ZhI\nEogQolcoBb74BIzM7O5t9T1Km9qRyy/D9sz/9F7tSCuq0QtlpVj27/P32/J5TYmjt0gCEUL0KiMm\ntmmrb2J4Rt021468shLPufnYnn2692tHWlHepn5bBw+AJ3p7bUkCEUL0OqXp+NLSICW121t9A4+Z\nkkrtfQ/S8PTiI7Ujf7gDrbQXa0dax+TzmXNtppdIAhFCmEIp8CUkorKywR6+MbfGj39ypHbkww9I\nMqN2pI+QBCKEMJURE4vKykKLjydsFXpH1Y54T2uqHbl6GpZdO8P0BAIkgQghIoBhseJLz4D+3a9e\nb/G4R9eOHCoj4Zrp5tSORClJIEKIiKDQ8CX1C0/1+tGaakcOv7IS99RfmV87EkUkgQghIoq/ej30\n2esdcsZTN+fWlrUjN/y/3ps7EoUkgQghIo5hteLLyELrF/qgqo4Eakd+fzuWr3b2/tyRKCIJRAgR\nkZSm4e3f/UFV7bJYcE9umjtybr45c0eigCQQIURECwyqig1P9frRVHIKtfc9aO7ckeOYJBAhRMQz\n7DEYWWGsXm8lIuaOHIckgQghjguB6vUwDKpqV6B25FVz544cRySBCCGOG0rhH1SVmQk2W488h5E7\noO3ckUceQquW2pHWJIEIIY47vtg4VFZO+Lf6Njt67sjUXxHzxkp/7cjf35bakaNIAhFCHJcCW33D\nMaiqAyre5a8dee5FjPQM4uf/wV87sveHHnm+440kECHEcSswqCo1fF192+M7eRhVS5b7a0e+3kXi\nFb/GseSZPl87IglECHFcUwp8rkRUZlaPXRcBjtSOvLISz3ljcDy3hMTLL+vTtSOSQIQQUcGIdTRd\nF3H06POo5BRq711A1ZNPg8Xirx25a26frB2RBCKEiBqG1YqR1TMtUFrznjPSXzsy87fYN28i6Vd9\nr3ZEEogQIrpoes+1QGnNbqfhqpl9du6IJBAhRFQKtECJCd+0w450OHckymtHJIEIIaKWYY9BZWWj\nxbvCN+2wI+3NHblsUg8/qbkkgQghopqhW/Clp4d92mGHjpo70vijs3v++UwkCUQIEfVaTDvs6esi\nTXynDKP2gYd65bnMIglECNFn+BxOjOwccPTsVt++QhKIEKJPMaw2jMysHmsN35dIAhFC9DmB1vAp\nPdsCJdqZ9sq98847TJw4kWHDhrFz55E90wcPHuSMM85g0qRJTJo0ifnz5we+t3PnTi666CLGjRvH\nggULTIhaCBEtlAJfQi+0QIlivXM1qR0nnXQSixYt4u67727zvQEDBrB69eo2t8+fP58FCxZw+umn\nM3PmTD788EN+9rOf9Ua4QogoZcQ60LNy0EuLUXV1ZodzXDFtBTJ48GAGDRqE6mRv/dLSUmprazn9\n9NMBuOSSS1i/fn1PhiiE6CMCreGTeq41fDSKyJN/Bw4c4NJLL+WKK65g+/btABQXF5ORkRH4mfT0\ndIqLi80KUQgRZZSm4UtJhbQ0sETkR2PE6dFTWDNmzKCsrKzN7XPmzCE/P7/d+6SlpbFp0yYSExPZ\nuXMns2bNYu3ate2uVLQgfyls3bqVbdu2Bb6ePXs2cXH2EH6L8LPZLBETC0g8xyLxBBd18ThS0FxO\nKCkGT2P3gtF1/vX553x01GfRyJEjycvL697jRogeTSDLly/v8n1sNhuJiYkAjBgxgtzcXH74ZJ3Y\nEgAADgdJREFU4QcyMjIoLCwM/FxxcTFpaWkdPk5eXl6bN6muztPleHpCXJw9YmIBiedYJJ7gojMe\nC3r/NPTSElRtbegPo+v8x5ln8h9RkjBai4h12tGri/LycgzDAGD//v3s27eP3NxcUlNTiY+PZ8eO\nHSileOONNxgzZoxZIQshopxh8V8XoV9/uS7SAdN2Ya1fv57777+fiooKfvvb33LKKaewdOlStm/f\nzpNPPonVakXXde677z4SEhIAuOeee5g7dy5ut5vRo0czevRos8IXQvQBCvD1T8YSY0crK0N5vWaH\nFFE01dltUFGgtDQyWitH55I/fCSe4CSe4HoqHt3jRispAXdDF+6k0/+MYb3SUt4MEXEKSwghIp1h\nj/FPO3S55JRWE0kgQgjRSUq34EtLh+Reag0f4eQVEEKILlBo+BKTpAUKkkCEECIkRqwDlZWD5nSa\nHYppJIEIIUSImlug9NWtvpJAhBCiG5q3+vbmtMNIIQlECCHCwBcXj5GVDbGxZofSaySBCCFEmBg2\nO0ZWdp+ZdigJRAghwigw7TA5BS3Ku/r2rRN2QgjRC5QCX2ISlhh7VK9Eojs9CiGEiXyxcWCPnFb3\n4SYJRAghREgkgQghhAiJJBAhhBAhkQQihBAiJJJAhBBChEQSiBBCiJBIAhFCCBESSSBCCCFCIglE\nCCFESCSBCCGECIkkECGEECGRBCKEECIkkkCEEEKERBKIEEKIkEgCEUIIERJJIEIIIUIiCUQIIURI\nJIEIIYQIiSQQIYQQIZEEIoQQIiSSQIQQQoREEogQQoiQSAIRQggREkkgQgghQiIJRAghREgkgQgh\nhAiJJBAhhBAhkQQihBAiJJJAhBBChEQSiBBCiJCYlkAefvhhLrjgAi6++GJmz55NTU1N4HvPPvss\n559/PhdccAEfffRR4PbNmzczfvx4xo0bx+LFi80IWwghRBPTEshPf/pT1q5dy5tvvsnAgQN59tln\nAfjuu+94++23WbduHUuWLOHee+9FKYVhGNx///0sW7aMNWvWsHbtWnbv3m1W+EII0eeZlkD+8z//\nE133P/2ZZ55JUVERABs3bmTChAlYrVZycnIYOHAgO3bsYMeOHQwcOJDs7GxsNhsXXnghGzZsMCt8\nIYTo8yLiGsiKFSv4+c9/DkBxcTGZmZmB76Wnp1NcXNzu7SUlJb0eqxBCCD9rTz74jBkzKCsra3P7\nnDlzyM/PB+Dpp5/GZrMxceJEAJRSbX5e0zQMw+jJUIUQQnRRjyaQ5cuXB/3+6tWr+eCDD3jxxRcD\nt2VkZFBYWBj4uqioiLS0NJRSFBQUBG4vLi4mLS2tw8feunUr27ZtC3w9e/ZsUlNdofwaPcLpjDE7\nhBYknuAknuAknuCeeuqpwL9HjhxJXl6eidGEkTLJBx98oCZMmKDKy8tb3P7tt9+qiy++WLndbrVv\n3z41duxYZRiG8nq9auzYserAgQPK7XarX/7yl+q7777r9PM9+eST4f4VQhZJsSgl8RyLxBOcxBNc\npMUTTj26AgnmgQceoLGxkauuugqAM844g/nz5zN06FAuuOACLrzwQqxWK/fccw+apmGxWLjrrru4\n6qqrUEoxZcoUhgwZYlb4QgjR55mWQN59990Ov3fddddx3XXXtbl99OjRjB49uifDEkII0UmW+fPn\nzzc7iN6Sk5NjdggBkRQLSDzHIvEEJ/EEF2nxhIumVDvbnoQQQohjiIg6ECGEEMcfSSBCCCFCYtpF\n9J5SVFTEbbfdRllZGRaLhalTpzJt2jQqKyuZM2cOBw8eJCcnh8cffxyXq+frQjweD5dffjmNjY34\nfD7GjRvH9ddfz4EDB7j55puprKxkxIgRPPzww1itvfd2GIbB5MmTSU9P55lnnjE1nvz8fOLj49F1\nHavVyooVK0x7vwCqq6u58847+fbbb9F1nQcffJBBgwaZEs+ePXuYM2cOmqahlGL//v3ceOONXHzx\nxabE8/zzz7NixQo0TeOkk05i4cKFlJSUmHbsvPDCC6xYsQLAlP/r8+bNY9OmTSQnJ/PWW28BBH3+\nBx54gM2bN+NwOHjooYcYNmxYj8TVa8zdRRx+JSUlateuXUoppWpqatT555+vvvvuO/Xwww+rxYsX\nK6WUevbZZ9UjjzzSazHV1dUppZTyer1q6tSp6vPPP1c33nijWrdunVJKqbvvvlv99a9/7bV4lFJq\n+fLl6pZbblHXXXedUkqZGk9+fr46fPhwi9vMfL9uv/12tWLFCqWUUo2NjaqqqsrUeJr5fD41atQo\nVVBQYEo8RUVFKj8/X7ndbqWU/5hZtWqVacfON998oyZOnKjcbrfyer1qxowZ6ocffujV1+bTTz9V\nu3btUhMnTgzc1tHzb9q0Sc2cOVMppdTnn3+upk6d2mNx9ZaoO4WVmpoayOpOp5MhQ4ZQXFzMhg0b\nmDRpEgCTJk1i/fr1vRaTw+EA/KsRr9eLpmls3bqVcePGBeJ57733ei2eoqIiPvjgA6ZOnRq4bcuW\nLabFo5q6LR/NrPerpqaG7du3M3nyZACsVisul8vU46fZxx9/zIABA8jMzDQtHsMwqK+vx+v10tDQ\nQFpammnH8u7duznzzDOx2+1YLBbOPvts3nvvPTZu3Nhrr83ZZ59NQkJCi9tavzfNTV83bNjAJZdc\nAvjr3qqrq9tt9XQ8iboEcrQDBw7w9ddfc8YZZ3Do0CFSUlIAf5KpqKjotTgMw+CSSy5h1KhRjBo1\nitzcXBISEgLdiDMyMnq1MeSDDz7IbbfdhqZpAFRUVJCYmGhaPJqmcfXVVzN58mRef/11ANPerwMH\nDtCvXz/mzp3LpEmTuOuuu6ivrzf1+Gm2bt26QM84M+JJT09nxowZnHvuuYwePRqXy8Xw4cNNO5ZP\nPPFEPv30UyorK6mvr2fz5s0UFRWZ/l6Vl5e3eP7y8nIASkpKyMjICPxcc6PY41nUJpDa2lpuuOEG\n5s2bh9PpDHxYmkHXdd544w02b97Mjh072p1j0lvxbdq0iZSUFIYNGxZoXKmUatPEsjdfr1deeYVV\nq1axZMkSXn75ZbZv327a++X1etm1axf/9V//xerVq3E4HCxevNjU4wegsbGRjRs3Mn78eKB3359m\nVVVVbNiwgffff58PP/ww8KHdWm/FNmTIEGbOnMmMGTO49tprOeWUU7BYLL3y3KFo/X8MzHkfwykq\nE4jX6+WGG27g4osvZuzYsQAkJycHloulpaX079+/1+OKj4/nnHPO4YsvvqCqqipw2qa5YWRv+Oc/\n/8nGjRsZM2YMt9xyC1u3buXBBx+kurralHjA/1caQP/+/Rk7diw7duww7f3KyMggIyOD0047DYDz\nzz+fXbt2mX78bN68mREjRgSe14x4Pv74Y3Jzc0lKSsJisTB27Fj+9a9/mXYsA0yePJlVq1bx0ksv\nkZiYyKBBg0x/rzp6/vT09MDcI+j916onRGUCmTdvHkOHDmX69OmB2/Lz81m1ahXg7wI8ZsyYXoml\nvLyc6upqABoaGvjkk08YOnQoeXl5vPPOO70ez80338ymTZvYsGEDjz32GHl5eTz66KOmxVNfX09t\nbS0AdXV1fPTRR5x00kmmvV8pKSlkZmayZ88ewH9taOjQoabF02zt2rWB01dgzvGclZXFF198gdvt\nRinFli1bOPHEE007doDA6aGCggLee+89Jk6c2OuvTeuVRUfPP2bMGN544w0APv/8cxISEgKnuo5X\nUVeJ/tlnn/Gb3/yGk046CU3T0DSNOXPmcPrpp3PTTTdRWFhIVlYWTzzxRJuLXz3h3//+N3fccQeG\nYWAYBhMmTOB3v/sd+/fv5+abb6aqqophw4bxyCOPYLPZejyeo23bto3nnnuOZ555xrR49u/fz/XX\nX4+mafh8Pi666CKuvfZaDh8+bMr7BfD1119z55134vV6yc3NZeHChfh8PtPiaWho4Nxzz2X9+vXE\nx8cDmPb6LFq0iLVr12K1Whk+fDgPPPAARUVFph3Ll19+OZWVlVitVubOnUteXl6vvjbNq/jDhw+T\nkpLC7NmzGTt2LDfeeGO7z3/ffffx4Ycf4nA4WLhwISNGjOiRuHpL1CUQIYQQvSMqT2EJIYToeZJA\nhBBChEQSiBBCiJBIAhFCCBESSSBCCCFCIglECCFESCSBCCGECIkkECGEECGJuoFSQnTVrFmzKCoq\nwuPxMG3aNKZOncrrr7/O0qVLSUxM5OSTTyYmJoY//OEPlJeXM3/+fAoLCwGYO3cuZ511lsm/gRDm\nkAQi+ryFCxeSkJCA2+1mypQpjB49mqeffpo333yTuLg4pk2bFpgxs2DBAq688krOOussCgsLufrq\nq1m3bp3Jv4EQ5pAEIvq8F154ITB0qKioiDfffJO8vLzAGNLx48ezd+9eAD755BO+//77QAO9uro6\n6urqiIuLMyd4IUwkCUT0adu2bWPLli28/vrr2O12rrjiCgYPHsz333/f7s8rpXj11Vex2+29HKkQ\nkUcuoos+rbq6moSEBOx2O7t37+aLL76gvr6eTz/9lOrqarxeL++++27g50eNGsVLL70U+Prrr782\nI2whIoJ04xV9msfjYdasWRQUFHDCCSdQVVXF7Nmz+eGHH1i6dClJSUkMHjyY9PR0brrpJioqKrjv\nvvvYvXs3hmFw9tlnM3/+fLN/DSFMIQlEiHY0X9fw+XzMmjWLKVOmBKZbCiH85BqIEO146qmn+OST\nT/B4PIwaNUqShxDtkBWIEEKIkMhFdCGEECGRBCKEECIkkkCEEEKERBKIEEKIkEgCEUIIERJJIEII\nIULy/wG5DNQ5Y7FfUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe0f0b8b990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def create_datapoint(age, sex, apoe):\n",
    "    v = 10.0\n",
    "    v = v - 0.8 * age\n",
    "    if sex == 1:\n",
    "        v = v + 1.2\n",
    "    else:\n",
    "        v = v - 1.2\n",
    "    if apoe == 2:\n",
    "        v = v - 1.5*age + 60\n",
    "    elif apoe == 1:\n",
    "        v = v - 9.5e-1*age  + 50\n",
    "    else:         \n",
    "        v = v - 0*age\n",
    "    return [v + 15 * (random.random() * 2) - 1, age, sex ,apoe]\n",
    "        \n",
    "    \n",
    "data = []\n",
    "# males\n",
    "for i in range(10):\n",
    "    data.append(create_datapoint(random.randint(45,75), 1, 2))\n",
    "    data.append(create_datapoint(random.randint(45,75), 1, 1))\n",
    "    data.append(create_datapoint(random.randint(45,75), 1, 0))\n",
    "# females\n",
    "for i in range(10):\n",
    "    data.append(create_datapoint(random.randint(45,75), 0, 2))\n",
    "    data.append(create_datapoint(random.randint(45,75), 0, 1))\n",
    "    data.append(create_datapoint(random.randint(45,75), 0, 0))   \n",
    "    \n",
    "data = pd.DataFrame(data, columns=['y','age','sex','apoe'])\n",
    "sns.lmplot('age','y', data, hue='apoe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:root:Used model for significance estimation: y ~ apoe_0_age + apoe_1_age + apoe_2_age+apoe_0 + apoe_1 + apoe_2 + sex + 1\n",
      "INFO:root:Used contrast: 1.0 * apoe_2_age - 0.5 * apoe_0_age - 0.5 * apoe_1_age - p-value: 0.005665363894268821\n",
      "INFO:root:Used contrast: 1.0 * apoe_0_age - 1.0 * apoe_2_age - p-value: 0.00020285322259629214\n",
      "INFO:root:Used contrast: 0.5 * apoe_1_age + 0.5 * apoe_2_age - 1.0 * apoe_0_age - p-value: 9.586285093932119e-05\n"
     ]
    }
   ],
   "source": [
    "from roistats.contrasts import genotypes as cg\n",
    "reload(cg)\n",
    "_ = cg.estimate(data, 'y', by='apoe', covariates=['sex'], interaction='age') #,contrasts = {'HT vs NC' : ('HT', 'NC')})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}