IHME_data_ex.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from src.IterativeFitting import IterativeFitting as IF\n",
    "from src.CorrFuncs import covariance_matrix, trend_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "study_df = pd.read_excel(\"SBPvPAD_data.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>114.0</td>\n",
       "      <td>21.000000</td>\n",
       "      <td>5428.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>1.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>128.0</td>\n",
       "      <td>39.000000</td>\n",
       "      <td>5410.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>0.410</td>\n",
       "      <td>0.2200</td>\n",
       "      <td>0.590</td>\n",
       "      <td>0.093434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>139.0</td>\n",
       "      <td>51.000000</td>\n",
       "      <td>5398.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>0.440</td>\n",
       "      <td>0.2700</td>\n",
       "      <td>0.620</td>\n",
       "      <td>0.088384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>159.0</td>\n",
       "      <td>116.000000</td>\n",
       "      <td>5333.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>1.080</td>\n",
       "      <td>0.9900</td>\n",
       "      <td>1.170</td>\n",
       "      <td>0.045455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>123.9</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>786.333333</td>\n",
       "      <td>1036.0</td>\n",
       "      <td>0.230</td>\n",
       "      <td>0.1600</td>\n",
       "      <td>0.310</td>\n",
       "      <td>0.037879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>129.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>1.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>134.6</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>139.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>143.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>149.5</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Author   dose       cases     controls  subjects  logOR     low  \\\n",
       "0  Stoekenbroek  114.0   21.000000  5428.500000    5449.5  1.000  1.0000   \n",
       "1  Stoekenbroek  128.0   39.000000  5410.500000    5449.5  0.410  0.2200   \n",
       "2  Stoekenbroek  139.0   51.000000  5398.500000    5449.5  0.440  0.2700   \n",
       "3  Stoekenbroek  159.0  116.000000  5333.500000    5449.5  1.080  0.9900   \n",
       "4         Itoga  123.9  249.666667   786.333333    1036.0  0.230  0.1600   \n",
       "5         Itoga  129.8  249.666667  2738.333333    2988.0  1.000  1.0000   \n",
       "6         Itoga  134.6  249.666667  6048.333333    6298.0 -0.062 -0.1500   \n",
       "7         Itoga  139.0  249.666667  9341.333333    9591.0  0.077  0.0004   \n",
       "8         Itoga  143.7  249.666667  8543.333333    8793.0  0.110  0.0220   \n",
       "9         Itoga  149.5  249.666667  4401.333333    4651.0  0.190  0.1100   \n",
       "\n",
       "    high  std_error  \n",
       "0  1.000   0.000000  \n",
       "1  0.590   0.093434  \n",
       "2  0.620   0.088384  \n",
       "3  1.170   0.045455  \n",
       "4  0.310   0.037879  \n",
       "5  1.000   0.000000  \n",
       "6 -0.023   0.032071  \n",
       "7  0.150   0.037778  \n",
       "8  0.210   0.047475  \n",
       "9  0.280   0.042929  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1.000\n",
       "1    0.410\n",
       "2    0.440\n",
       "3    1.080\n",
       "4    0.230\n",
       "5    1.000\n",
       "6   -0.062\n",
       "7    0.077\n",
       "8    0.110\n",
       "9    0.190\n",
       "Name: logOR, dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df[\"logOR\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>114.0</td>\n",
       "      <td>21.000000</td>\n",
       "      <td>5428.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>128.0</td>\n",
       "      <td>39.000000</td>\n",
       "      <td>5410.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>0.410</td>\n",
       "      <td>0.2200</td>\n",
       "      <td>0.590</td>\n",
       "      <td>0.093434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>139.0</td>\n",
       "      <td>51.000000</td>\n",
       "      <td>5398.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>0.440</td>\n",
       "      <td>0.2700</td>\n",
       "      <td>0.620</td>\n",
       "      <td>0.088384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Stoekenbroek</td>\n",
       "      <td>159.0</td>\n",
       "      <td>116.000000</td>\n",
       "      <td>5333.500000</td>\n",
       "      <td>5449.5</td>\n",
       "      <td>1.080</td>\n",
       "      <td>0.9900</td>\n",
       "      <td>1.170</td>\n",
       "      <td>0.045455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>123.9</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>786.333333</td>\n",
       "      <td>1036.0</td>\n",
       "      <td>0.230</td>\n",
       "      <td>0.1600</td>\n",
       "      <td>0.310</td>\n",
       "      <td>0.037879</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>129.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>134.6</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>139.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>143.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>149.5</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Author   dose       cases     controls  subjects  logOR     low  \\\n",
       "0  Stoekenbroek  114.0   21.000000  5428.500000    5449.5  0.000  0.0000   \n",
       "1  Stoekenbroek  128.0   39.000000  5410.500000    5449.5  0.410  0.2200   \n",
       "2  Stoekenbroek  139.0   51.000000  5398.500000    5449.5  0.440  0.2700   \n",
       "3  Stoekenbroek  159.0  116.000000  5333.500000    5449.5  1.080  0.9900   \n",
       "4         Itoga  123.9  249.666667   786.333333    1036.0  0.230  0.1600   \n",
       "5         Itoga  129.8  249.666667  2738.333333    2988.0  0.000  0.0000   \n",
       "6         Itoga  134.6  249.666667  6048.333333    6298.0 -0.062 -0.1500   \n",
       "7         Itoga  139.0  249.666667  9341.333333    9591.0  0.077  0.0004   \n",
       "8         Itoga  143.7  249.666667  8543.333333    8793.0  0.110  0.0220   \n",
       "9         Itoga  149.5  249.666667  4401.333333    4651.0  0.190  0.1100   \n",
       "\n",
       "    high  std_error  \n",
       "0  0.000   0.000000  \n",
       "1  0.590   0.093434  \n",
       "2  0.620   0.088384  \n",
       "3  1.170   0.045455  \n",
       "4  0.310   0.037879  \n",
       "5  0.000   0.000000  \n",
       "6 -0.023   0.032071  \n",
       "7  0.150   0.037778  \n",
       "8  0.210   0.047475  \n",
       "9  0.280   0.042929  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df.iloc[0,5:8] = np.log(study_df.iloc[0,5:8].to_numpy().astype(np.float64))\n",
    "study_df.iloc[5,5:8] = np.log(study_df.iloc[5,5:8].to_numpy().astype(np.float64))\n",
    "study_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "study_df_s = study_df.loc[study_df[\"Author\"] == \"Stoekenbroek\"]\n",
    "study_df_i = study_df.loc[study_df[\"Author\"] == \"Itoga\"]\n",
    "study_df_i = study_df_i.iloc[1:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roll_func(study_df_j,col):\n",
    "    col_s_run = []\n",
    "    col = study_df_j[col].to_numpy().astype(np.float64)\n",
    "    for i in range(col.shape[0]-1):\n",
    "        col_s_run.append(col[i+1]-col[0])\n",
    "    return np.array(col_s_run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_s_run = roll_func(study_df_j=study_df_s, col=\"dose\")\n",
    "OR_s_rise = roll_func(study_df_j=study_df_s, col=\"logOR\")\n",
    "low_s_rise = roll_func(study_df_s,\"low\")\n",
    "high_s_rise = roll_func(study_df_s,\"high\")\n",
    "\n",
    "std_s_rise = roll_func(study_df_s, \"std_error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_slopes_s = OR_s_rise / dose_s_run\n",
    "low_slopes_s = low_s_rise / dose_s_run\n",
    "high_slopes_s = high_s_rise / dose_s_run\n",
    "\n",
    "std_s_rise = std_s_rise / dose_s_run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma_s = np.diag(std_s_rise)**2\n",
    "inv_Sigma_s = np.linalg.inv(Sigma_s)\n",
    "beta_s = np.linalg.inv(np.ones((3,1)).T@inv_Sigma_s@np.ones((3,1)))@np.ones((3,1)).T@inv_Sigma_s@log_slopes_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_beta_s = np.linalg.inv(np.ones((3,1)).T@inv_Sigma_s@np.ones((3,1)))@np.ones((3,1)).T@inv_Sigma_s@low_slopes_s\n",
    "high_beta_s = np.linalg.inv(np.ones((3,1)).T@inv_Sigma_s@np.ones((3,1)))@np.ones((3,1)).T@inv_Sigma_s@high_slopes_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_covar_matrix(study_df,std_rise,Lx,dose_index=1):\n",
    "    Nx_s = study_df[\"subjects\"].to_numpy().astype(np.float64)\n",
    "    Lx_s = study_df[Lx].to_numpy().astype(np.float64)[dose_index:]\n",
    "    M1x_s = np.sum(study_df[\"cases\"].to_numpy().astype(np.float64))\n",
    "    A0_s = M1x_s*Nx_s[1:]/(Nx_s.sum())\n",
    "\n",
    "    it_fit_ex = IF(Lx_s,A0_s,Nx_s,M1x_s)\n",
    "    A_s, B_s, a0_s, b0_s = it_fit_ex.convexProgram()\n",
    "\n",
    "    C = covariance_matrix(A_s,B_s,a0_s,b0_s,std_rise**2)\n",
    "    inv_C = np.linalg.inv(C)\n",
    "\n",
    "    return C, inv_C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/reductions/solvers/solving_chain.py:354: FutureWarning: \n",
      "    You specified your problem should be solved by ECOS. Starting in\n",
      "    CXVPY 1.6.0, ECOS will no longer be installed by default with CVXPY.\n",
      "    Please either add an explicit dependency on ECOS or switch to our new\n",
      "    default solver, Clarabel, by either not specifying a solver argument\n",
      "    or specifying ``solver=cp.CLARABEL``.\n",
      "    \n",
      "  warnings.warn(ECOS_DEP_DEPRECATION_MSG, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "C_s, inv_C_s = create_covar_matrix(study_df_s,std_s_rise,\"logOR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_C_s, low_inv_C_s = create_covar_matrix(study_df_s,std_s_rise,\"low\")\n",
    "\n",
    "high_C_s, high_inv_C_s = create_covar_matrix(study_df_s,std_s_rise,\"high\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_s_GL = np.linalg.inv(np.ones((3,1)).T@inv_C_s@np.ones((3,1)))@np.ones((3,1)).T@inv_C_s@log_slopes_s\n",
    "low_beta_s_GL = np.linalg.inv(np.ones((3,1)).T@low_inv_C_s@np.ones((3,1)))@np.ones((3,1)).T@low_inv_C_s@log_slopes_s\n",
    "high_beta_s_GL = np.linalg.inv(np.ones((3,1)).T@high_inv_C_s@np.ones((3,1)))@np.ones((3,1)).T@high_inv_C_s@log_slopes_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02414612])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_s_GL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_i_run = roll_func(study_df_j=study_df_i, col=\"dose\")\n",
    "OR_i_rise = roll_func(study_df_j=study_df_i, col=\"logOR\")\n",
    "low_i_rise = roll_func(study_df_i,\"low\")\n",
    "high_i_rise = roll_func(study_df_i,\"high\")\n",
    "\n",
    "std_i_rise = roll_func(study_df_i, \"std_error\")\n",
    "\n",
    "log_slopes_i = OR_i_rise / dose_i_run\n",
    "low_slopes_i = low_i_rise / dose_i_run\n",
    "high_slopes_i = high_i_rise / dose_i_run\n",
    "\n",
    "std_i_rise = std_i_rise / dose_i_run\n",
    "\n",
    "Sigma_i = np.diag(std_i_rise)**2\n",
    "inv_Sigma_i = np.linalg.inv(Sigma_i)\n",
    "beta_i = np.linalg.inv(np.ones((4,1)).T@inv_Sigma_i@np.ones((4,1)))@np.ones((4,1)).T@inv_Sigma_i@log_slopes_i\n",
    "\n",
    "C_i, inv_C_i = create_covar_matrix(study_df_i,std_i_rise,\"logOR\")\n",
    "\n",
    "beta_i_GL = np.linalg.inv(np.ones((4,1)).T@inv_C_i@np.ones((4,1)))@np.ones((4,1)).T@inv_C_i@log_slopes_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.01291667,  0.00836957,  0.00791367,  0.00964467])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_slopes_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4.46410700e-05, 1.97051483e-05, 1.62839694e-05, 9.57345970e-06],\n",
       "       [1.97051483e-05, 1.68615370e-05, 1.08005842e-05, 6.34973910e-06],\n",
       "       [1.62839694e-05, 9.57345970e-06, 1.16652950e-05, 5.24730674e-06],\n",
       "       [1.08005842e-05, 6.34973910e-06, 5.24730674e-06, 4.74870311e-06]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4.46410700e-05, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 1.68615370e-05, 0.00000000e+00, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 1.16652950e-05, 0.00000000e+00],\n",
       "       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 4.74870311e-06]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sigma_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_fn(beta,x):\n",
    "    return beta*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(115,165)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x17e885550>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACRSklEQVR4nOzdeXxU5fU/8M+sd/YtM1kJIRCIbGGzfFkqEAUBq+BWl9IKLrQqtaJ1Kf7EBZdAXbFYcavYVuoutihWVAIIqMhSUHbMwpI9s+8z9/7+mJnLTPZJJplkct6vV14md+6982QMmZPnOec5Ao7jOBBCCCGEpAhhsgdACCGEEJJIFNwQQgghJKVQcEMIIYSQlELBDSGEEEJSCgU3hBBCCEkpFNwQQgghJKVQcEMIIYSQlELBDSGEEEJSCgU3hBBCCEkpFNwQQgghJKX0muBm5cqVEAgEWLp0aavnrFu3DgKBIOZDJpP13CAJIYQQ0uuJkz0AANi9ezdefvllFBUVtXuuRqPB0aNH+a8FAkF3Do0QQgghfUzSZ24cDgcWLFiAV199FXq9vt3zBQIBMjMz+Y+MjIweGCUhhBBC+oqkz9wsWbIEv/jFLzBz5kw8/vjj7Z7vcDiQl5cHlmUxfvx4PPnkkxg5cmSr53u9Xni9Xv5rlmXR2NiItLQ0mvUhhBBC+giO42C325GdnQ2hsO25maQGN2+//Tb27t2L3bt3d+j8wsJC/O1vf0NRURGsViuefvppTJkyBT/++CMGDBjQ4jUlJSV49NFHEzlsQgghhCTJqVOnWn3PjxBwHMf10HhinDp1Cueffz42b97M59rMmDEDY8eOxfPPP9+he/j9fgwfPhzXX389HnvssRbPaTpzY7VaMXDgQJw6dQoajabL3wchhBBCup/NZkNubi4sFgu0Wm2b5yZt5mbPnj2ora3F+PHj+WPBYBDbtm3DmjVr4PV6IRKJ2ryHRCLBuHHjcOLEiVbPYRgGDMM0O67RaCi4IYQQQvqYjqSUJC24ueiii3Dw4MGYYzfeeCPOO+883H///e0GNkAoGDp48CAuueSS7homIYQQQvqYpAU3arUao0aNijmmVCqRlpbGH7/hhhuQk5ODkpISAMCKFSswadIkFBQUwGKx4KmnnkJFRQVuueWWHh8/IYQQQnqnpFdLtaWysjImI9psNmPx4sWorq6GXq/HhAkTsHPnTowYMSKJoySEEEJIb5K0hOJksdls0Gq1sFqtbebcBINB+P3+HhwZSQUSiaRDS6qEEELi09H3b6CXz9wkA8dxqK6uhsViSfZQSB+l0+mQmZlJ+ygRQkiSUHDTRCSwSU9Ph0KhoDco0mEcx8HlcqG2thYAkJWVleQREUJI/0TBTZRgMMgHNmlpackeDumD5HI5AKC2thbp6em0REUIIUmQ9N5SvUkkx0ahUCR5JKQvi/z8UM4WIYQkBwU3LaClKNIV9PNDCCHJRcENIYQQQlIKBTekxy1atAiXX355XNcMGjSowz3HCCGE9G8U3KSIRYsWQSAQYOXKlTHHN2zY0OuWSVavXo1169Yl9J7l5eUQCATYv39/Qu8LAHV1dbjtttswcOBAMAyDzMxMzJ49Gzt27Ej4cxFCCOk6qpZKITKZDKtWrcLvfvc76PX6ZA+nVe11c+1trrrqKvh8Prz55psYPHgwampq8OWXX6KhoSHZQyOEENICmrlJITNnzkRmZibfi6s1H3zwAUaOHAmGYTBo0CA888wzrZ5rtVohEonw/fffAwBYloXBYMCkSZP4c/75z38iNzeX//rUqVO45pproNPpYDAYMH/+fJSXl/OPN12WstvtWLBgAZRKJbKysvDcc89hxowZWLp0acxYXC4XbrrpJqjVagwcOBCvvPIK/1h+fj4AYNy4cRAIBJgxYwYAoLS0FBMnToRSqYROp8PUqVNRUVHR5usTzWKxYPv27Vi1ahWKi4uRl5eHiRMnYtmyZZg3b16H70MIIf1FIOCB3+9K6hgouGkHx3Hw+ZxJ+Yi3M4ZIJMKTTz6Jv/zlLzh9+nSL5+zZswfXXHMNrrvuOhw8eBCPPPIIli9f3uoykVarxdixY1FaWgoAOHjwIAQCAfbt2weHwwEA2Lp1K6ZPnw4gVP48e/ZsqNVqbN++HTt27IBKpcKcOXPg8/lafI67774bO3bswL///W9s3rwZ27dvx969e5ud98wzz+D888/Hvn37cPvtt+O2227D0aNHAQDfffcdAOCLL75AVVUVPvzwQwQCAVx++eWYPn06Dhw4gF27duG3v/1tXMt0KpUKKpUKGzZsgNfr7fB1hBDSX3AcC5erAXV1h3Hy5GYcPLgeZWVfJXVMtCzVDr/fhZISVVKee9kyB6RSZVzXXHHFFRg7diwefvhhvP76680ef/bZZ3HRRRdh+fLlAIBhw4bh0KFDeOqpp7Bo0aIW7zljxgyUlpbinnvuQWlpKWbNmoUjR47g66+/xpw5c1BaWor77rsPAPDOO++AZVm89tprfBDxxhtvQKfTobS0FBdffHHMve12O958802sX78eF110EX9+dnZ2s3FccskluP322wEA999/P5577jls2bIFhYWFMJlMAIC0tDRkZmYCABobG2G1WnHppZdiyJAhAIDhw4fH9XqKxWKsW7cOixcvxtq1azF+/HhMnz4d1113HYqKiuK6FyGEpIpAwAOnsxYORw0sljK4XA3w+50QCsVgWRZSqTqp46OZmxS0atUqvPnmmzh8+HCzxw4fPoypU6fGHJs6dSqOHz+OYDDY4v2mT5+Or7/+GsFgEFu3bsWMGTP4gOfs2bM4ceIEvwz0v//9DydOnIBareZnPQwGAzweD06ePNns3j/99BP8fj8mTpzIH9NqtSgsLGx2bnQwIRAIkJmZybc6aInBYMCiRYswe/ZsXHbZZVi9ejWqqqoAhDrOR8anUqnw5JNPYvv27THH3nrrLQChnJuzZ8/i3//+Nx/MjR8/PuFJ0YQQ0luFZmfqUVd3CCdPfo6DB9fj0KEPUFGxDS5XI2QyHQyGYdDrh8T9R3l3oJmbdkgkCixb5kjac3fGtGnTMHv2bCxbtqzV2Zh472e327F3715s27YNTz75JDIzM7Fy5UqMGTMG2dnZGDp0KADA4XBgwoQJfGAQLTK70lkSiSTma4FAAJZl27zmjTfewB/+8Ad89tlneOedd/Dggw9i8+bNOP/882MqqwwGA+RyecyxjIwM/nOZTIZZs2Zh1qxZWL58OW655RY8/PDDCXl9CSGkN/L73TGzM253A3w+F0QiCWQyLXS6QRAKe2cY0TtH1YsIBIJeEYXGa+XKlRg7dmyzGZDhw4c3K2HesWMHhg0b1mofJJ1Oh6KiIqxZswYSiQTnnXce0tPTce2112Ljxo18vg0AjB8/Hu+88w7S09PbbUkPAIMHD4ZEIsHu3bsxcOBAAKEk5mPHjmHatGkd/n6lUikAtDj7NG7cOIwbNw7Lli3D5MmTsX79ekyaNAkFBQXNzm3pWEtGjBiBDRs2dHh8hBDS20VyZ5zOWthsp2G3n4HHYwHHhZaZZDID1OoBvW57kZbQslSKGj16NBYsWIAXXngh5vgf//hHfPnll3jsscdw7NgxvPnmm1izZg3uueeeNu83Y8YMvPXWW3wgYzAYMHz4cLzzzjsxwc2CBQtgNBoxf/58bN++HWVlZSgtLcUf/vCHFpOc1Wo1Fi5ciHvvvRdbtmzBjz/+iJtvvhlCoTCuf0Dp6emQy+X47LPPUFNTA6vVirKyMixbtgy7du1CRUUFPv/8cxw/fjyuvJuGhgZceOGF+Oc//4kDBw6grKwM7733Hv785z9j/vz5Hb4PIYT0Rn6/CxZLOU6f/gY//vgODh16FydObEJ9/VEIBCLodPlISyuEWp0NqVTZJwIbgIKblLZixYpmyzbjx4/Hu+++i7fffhujRo3CQw89hBUrVrS7vDJ9+nQEg0E+twYIBTxNjykUCmzbtg0DBw7ElVdeieHDh+Pmm2+Gx+NpdSbn2WefxeTJk3HppZdi5syZmDp1KoYPHw6ZTNbh71UsFuOFF17Ayy+/jOzsbMyfPx8KhQJHjhzBVVddhWHDhuG3v/0tlixZgt/97ncdvq9KpcL//d//4bnnnsO0adMwatQoLF++HIsXL8aaNWs6fB9CCOkNWDYIp7MOtbU/4vjxTTh48F84fPhDVFbugMdjg1xuRFpaIfT6fMjlhl677NQeARdvvXEfZ7PZoNVqYbVam73ZejwelJWVIT8/P643VpJYTqcTOTk5eOaZZ3DzzTcnezhxo58jQkhv4vM5w7kz1bBYyuF2N8Lvd0EkYiCT6cAwGgiFLacldIbdfhZyuQEjR/4yYfcE2n7/bqpvhmQkpezbtw9HjhzBxIkTYbVasWLFCgCgZR9CCOkElg3C7Q7lzlitlbDbz8LrtYHjOEilaigUxk4XrPQVFNyQXuHpp5/G0aNHIZVKMWHCBGzfvh1GozHZwyKEkD7B53PA6ayF3V4Fi6UcHo8ZgYAHIhEDhtFCpxuc0NmZ3o6CG5J048aNw549e5I9DEII6TNYNgCXqx5OZx2s1grY7VXweu0AODCMBgpFOiQSebKHmTQU3BBCCCF9gNdrD+fOVIVzZ8wIBLwQi2WQybRQKEz9anamLRTcEEIIIb1QZHbG4aiB1VoJh6MaXq8NgAAMo4ZKlQmxmIoWWkLBDSGEENJLeL22cO7MWVgsFeHcGR8kEjkYhmZnOoqCG0IIISRJgkF/OHemFhZLOZzOGni99vDu+BqoVNkQi5lkD7PPoeCGEEII6SEcx/GzMw5HaHbG7baAZX0Qi0OzM0plBgQC2mO3Kyi4IYQQQrpRaHamLjw7UwansxZerwMCgRAMo4FaTbMziUbBDWlRaWkpiouLYTabodPpOnTNI488gg0bNsR01iaEkP4mNDtjDTegPAubLTQ7Ewz6IZHIIZPpoFRm0uxMN6JXNgWsXbsWarUagUCAP+ZwOCCRSGL6PgGhoEUgEODkyZNt3nPKlCmoqqqCVqtN6FhnzJiBpUuXJvSeEa+++irGjBkDlUoFnU6HcePGoaSkpFueixBCogWDPthsZ1BVtQ+HD3+EH354B0eP/hvV1fsQDAah0QyA0VgIrXYgGEZDgU03o5mbFFBcXAyHw4Hvv/8ekyZNAgBs374dmZmZ+Pbbb+HxePgeR1u2bMHAgQMxZMiQNu8plUqRmZnZ7WNPlL/97W9YunQpXnjhBUyfPh1erxcHDhzADz/8kOyhEUJSEMdx8HgsMZVNXq8FLBuAWKyATKaFSkWzM8lCr3oKKCwsRFZWFkpLS/ljpaWlmD9/PvLz8/HNN9/EHC8uLgbLsigpKUF+fj7kcjnGjBmD999/P+Y8gUAAi8XCH3v11VeRm5sLhUKBK664As8++2yLS1b/+Mc/MGjQIGi1Wlx33XWw2+0AgEWLFmHr1q1YvXo1BAIBBAIBysvLYTabsWDBAphMJsjlcgwdOhRvvPFGXK/Bv//9b1xzzTW4+eabUVBQgJEjR+L666/HE088Edd9CCGkNYGAFzbbaZw9uweHD3+IH398B8eObURNzf/AcSw0mlykpRVCq82l2Zkko5mbdnAcB1+g/fO6g1QMCASCDp1bXFyMLVu24E9/+hOA0AzNfffdh2AwiC1btmDGjBlwu9349ttvcdNNN6GkpAT//Oc/sXbtWgwdOhTbtm3Dr3/9a5hMJkyfPr3Z/Xfs2IFbb70Vq1atwrx58/DFF19g+fLlzc47efIkNmzYgI0bN8JsNuOaa67BypUr8cQTT2D16tU4duwYRo0axTfHNJlMuPPOO3Ho0CFs2rQJRqMRJ06cgNvtjuu1yszMxNatW1FRUYG8vLy4riWEkJaEZmfM4dyZM7BaT8HrtSAYDEAqVYJhtFCpsjv8e5r0HApu2uELAL9/1ZyU516zWA9G0rFzi4uLsXTpUgQCAbjdbuzbtw/Tp0+H3+/H2rVrAQC7du2C1+vFjBkzMGLECHzxxReYPHkyAGDw4MH4+uuv8fLLL7cY3PzlL3/B3Llzcc899wAAhg0bhp07d2Ljxo0x57Esi3Xr1kGtVgMAfvOb3+DLL7/EE088Aa1WC6lUCoVCEbPkVVlZiXHjxuH8888HAAwaNCiu1wkAHn74YVx55ZUYNGgQhg0bhsmTJ+OSSy7B1VdfDaGQ/noihHRMIOCB01kHh6MGFksZXK56+HxOiERiMIwWGs1AiEQd/MVMkoZ+66eIGTNmwOl0Yvfu3di+fTuGDRvGz8JE8m5KS0sxePBgOBwOuFwuzJo1CyqViv/4+9//3mqi8dGjRzFx4sSYY02/BkKBSSSwAYCsrCzU1ta2OfbbbrsNb7/9NsaOHYv77rsPO3fu5B+bO3cuP76RI0cCAEaOHMkfmzt3Lv88u3btwsGDB3HnnXciEAhg4cKFmDNnDliW7diLSAjpdziOhcvVgLq6wzh5cjMOHvwXDh/+ABUVpXC56iGT6ZCWNgx6/RAoFEYKbPoImrlph1QcmkFJ1nN3VEFBAQYMGIAtW7bAbDbzsy/Z2dnIzc3Fzp07sWXLFlx44YVwOBwAgE8++QQ5OTkx92GYru21IJHE/sMXCATtBhdz585FRUUFPv30U2zevBkXXXQRlixZgqeffhqvvfYav0QVufenn34Kv98PAJDLY7vejho1CqNGjcLtt9+OW2+9FRdccAG2bt2K4uLiLn1fhJDUEZqdqY2anWmA3++EUCgGw+hodiYFUHDTDoFA0OGloWQrLi5GaWkpzGYz7r33Xv74tGnTsGnTJnz33Xe47bbbMGLECDAMg8rKyhaXoFpSWFiI3bt3xxxr+nVHSKVSBIPBZsdNJhMWLlyIhQsX4oILLsC9996Lp59+ulnwBaDDOTUjRowAADidzrjHSQhJHRzHwu1u5HNnbLZT8Hhs4LgAJBIVZDId1Oocyp1JIRTcpJDi4mIsWbIEfr8/JmiZPn06fv/738Pn86G4uBhqtRr33HMP7rrrLrAsi5///OewWq3YsWMHNBoNFi5c2Ozed9xxB6ZNm4Znn30Wl112Gb766its2rQp7l8GgwYNwrfffovy8nKoVCoYDAY88sgjmDBhAkaOHAmv14uNGzdi+PDhcd33tttuQ3Z2Ni688EIMGDAAVVVVePzxx2Eymfi8IkJI/+H3u2NmZ9zuBvh8LohEEshkWuh0eRAK6S0wVdH/2RRSXFwMt9uN8847DxkZGfzx6dOnw2638yXjAPDYY4/BZDKhpKQEP/30E3Q6HcaPH48HHnigxXtPnToVa9euxaOPPooHH3wQs2fPxl133YU1a9bENcZ77rkHCxcuxIgRI+B2u1FWVgapVIply5ahvLwccrkcF1xwAd5+++247jtz5kz87W9/w0svvYSGhgYYjUZMnjwZX375JdLS0uK6FyGk74nkzoT2nTkDm+00PB4rOC4IqVQNmcwAtXoAzc70EwKO47hkDwIAVq5ciWXLluHOO+/E888/3+p57733HpYvX47y8nIMHToUq1atwiWXXNLh57HZbNBqtbBardBoNDGPeTwelJWVIT8/n9/0jrRu8eLFOHLkCLZv357sofQq9HNESM/w+13h2Znq8OxMI/x+N4RCKWQyLRhGQ7MzSWC3n4VcbsDIkb9M6H3bev9uqlf8X9+9ezdefvllFBUVtXnezp07cf3116OkpASXXnop1q9fj8svvxx79+7FqFGjemi0/dfTTz+NWbNmQalUYtOmTXjzzTfx17/+NdnDIoT0E9GzM1ZrJez2s/B6reA4DlKpCnK5ERqNItnDJL1A0oMbh8OBBQsW4NVXX8Xjjz/e5rmrV6/GnDlz+GTZxx57DJs3b8aaNWv4vVxI9/nuu+/w5z//GXa7HYMHD8YLL7yAW265JdnDIoSkMJ/PGTU7Ux6enXFBJGIgk+mg0w2GUChK9jBJL5P04GbJkiX4xS9+gZkzZ7Yb3OzatQt33313zLHZs2djw4YNrV7j9Xrh9Xr5r202W5fG25+9++67yR4CISTFsWwQbnfT2RlbeHZGDYXCCImEZmdI25Ia3Lz99tvYu3dvh0uKq6urYxJlASAjIwPV1dWtXlNSUoJHH320S+MkhBDSfXw+Bz87YzaXweMxw+93QyxmwDA0O9OX+AJCVFoyIXOrMDKJ40hacHPq1Cnceeed2Lx5c7cmXS5btixmtsdmsyE3N7fbno8QQkjbWDYQlTtTAbu9Gl6vDQIBwrMz6ZBI5O3fiCRdkBXgrFWNSrMWFY06VNtUYDkhMtRWXHJB8saVtOBmz549qK2txfjx4/ljwWAQ27Ztw5o1a+D1eiESxUbqmZmZqKmpiTlWU1MT06eoKYZhurzrLiGEkK7xeu3h2ZmqcO6MGcGgL5w7o4FCQbMzfQHHAbV2JSrMOlQ2anHaokGAjf3/JhN7IJf4wHFc0krvkxbcXHTRRTh48GDMsRtvvBHnnXce7r///maBDQB+35KlS5fyxzZv3kybtBFCSC8Tmp2ph9NZC4ulAg5HaHYGEIBh1FCpMiEW01YJvR3HAWaXHBVmLSobtThl1sITiN22XyoKIE3pQobagTyDFTrxMejUKggEhUkadRKDG7Va3ax8W6lUIi0tjT9+ww03ICcnByUlJQCAO++8E9OnT8czzzyDX/ziF3j77bfx/fff45VXXunx8RNCCInl9drCm+iFZmc8HjMCAR8kEjkYRgOFwkSzM32A3SPll5kqG7Vw+GJXP0TCIAwKN9LVDuTqbBhksEDB+CEMT9LY7T4ke6/EpFdLtaWyshJC4bnG5VOmTMH69evx4IMP4oEHHsDQoUOxYcMG2uOGEEKSIBj0R83OlMPprIHXa4dAIIBUqoFKlQ2xmNICeju3X4xT4ZmZCrMWZldsNZpAwEIv98CkciJHZ8NgYyPUjB8iYa/YA7hFvSq4KS0tbfNrAPjlL3+JX/4ysbsekp61aNEiWCyWNkv4mxo0aBCWLl0asyRJCOl5Ho81nDtzFhZLJdxuM1jWB7FYDobRQqnMgEAgbP9GJGl8QSHOWDThYEaHWrsSQPRUCwetzAujyolMjQNDjQ3QKryQiNhkDTluvSq4IZ23aNEivPnmmygpKcGf/vQn/viGDRtwxRVXoJd02QAQ2owx0eMpLy9Hfn4+9u3bh7Fjxyb03p0JxghJFaHZmbrw7EwZnM5aeL0OCARCMIwGajXNzvR2QVaAapsqtMxk1uKsVQ2Wiw1AlVIvjKpQ3swQYyNMKhek4r4TzDRFwU0KkclkWLVqFX73u99Br9cnezit0mq1yR4CIaQVHMfB67XC6ayD3X4WVmsF3G4LgsFQ7oxMpoNSmUmzM70YxwF1DiUqGrWoNGtx2qKFP9ikoknih1HpQobajkFpFmRpHGDEwYTkyni9NgQCnq7fqAvopzOFzJw5E5mZmXwCdms++OADjBw5EgzDYNCgQXjmmWdaPddqtUIkEuH7778HALAsC4PBgEmTJvHn/POf/4zZO+jUqVO45pproNPpYDAYMH/+fJSXl/OPL1q0CJdffjn/td1ux4IFC6BUKpGVlYXnnnsOM2bMaLYE5XK5cNNNN0GtVmPgwIExieT5+fkAgHHjxkEgEGDGjBkAQkubEydOhFKphE6nw9SpU1FRUdHm60NIfxMM+mCznUFV1T4cObIBP/zwDo4e/RhVVXsRDAagVmfDaCyEVjsQDKOhwKaXCVU0yfC/0xn498FC/HX7RPz9u7HYeiIfZQ0G+IMiSERBZGpsKMquwmWjDuPG/9uLa8f/gOJhFchPs0Im6Xxg4/Xacfr0LuzZ8zI2bfo9vvjiPhw48I/EfpNxopmbdnAcB7c/mJTnlktEce0RIBKJ8OSTT+JXv/oV/vCHP2DAgAHNztmzZw+uueYaPPLII7j22muxc+dO3H777UhLS8OiRYuana/VajF27FiUlpbi/PPPx8GDByEQCLBv3z44HA6oVCps3boV06dPBwD4/X7Mnj0bkydPxvbt2yEWi/H4449jzpw5OHDgAKRSabPnuPvuu7Fjxw78+9//RkZGBh566CHs3bu32fLSM888g8ceewwPPPAA3n//fdx2222YPn06CgsL8d1332HixIn44osvMHLkSEilUgQCAVx++eVYvHgx/vWvf8Hn8+G7775L2r4LhPQWHMfB47HA5aqDzXYWFks5vF4rWNYPsVgBmUwLlYpmZ3ozh1fKJwBXNupg9zataGJhULiQHk4Czk8zQyn1Q5iA/6V+vwt1dT+ipuYAamsPwmwuA8A1OyeZKLhph9sfxIiH/puU5z60YjYU0vj+F11xxRUYO3YsHn74Ybz++uvNHn/22Wdx0UUXYfny5QCAYcOG4dChQ3jqqadaDG4AYMaMGSgtLcU999yD0tJSzJo1C0eOHMHXX3+NOXPmoLS0FPfddx8A4J133gHLsnjttdf4IOKNN96ATqdDaWkpLr744ph72+12vPnmm1i/fj0uuugi/vzs7Oxm47jkkktw++23AwDuv/9+PPfcc9iyZQsKCwthMpkAAGlpafymjo2NjbBarbj00ksxZMgQAMDw4cPjej0JSRWBgJfPnTGby+By1cHnc0IoFEEq1UCjGQCRqPkfH6R38PhFOGUJVzQ16tDYSkWTUeVEtsaOAlMj1DJfQiqaAgEv6usP8cFMY+MJcFxsPo5Gk4v09NHIyCiCQmGEVjuwy8/bFRTcpKBVq1bhwgsvxD333NPsscOHD2P+/Pkxx6ZOnYrnn38ewWCwxc0Tp0+fjtdffx3BYBBbt27FxRdfjMzMTJSWlqKoqAgnTpzgl4H+97//4cSJE1Cr1TH38Hg8OHnyZLN7//TTT/D7/Zg4cSJ/TKvVorCw+eZPRUVF/OcCgQCZmZmora1t9XUwGAxYtGgRZs+ejVmzZmHmzJm45pprkJWVhcrKSowYMYI/94EHHsAFF1yAuXPn8sdefvllLFiwoNX7E9KbhWZnzOHcmTOwWCrh9VoQDAYglSrBMFqoVNk0k9lL+YNCnLGoUWnWoaJRi1q7ClyTiiaNzAuT0okMjQNDTQ3QJaiiKRj0o6HhaDiYOYCGhmNg2UDMOSpVJtLTi5CRMRrp6aMhlxv4x+z2s10eQ1dRcNMOuUSEQytmJ+25O2PatGmYPXs2li1b1upsTLz3s9vt2Lt3L7Zt24Ynn3wSmZmZWLlyJcaMGYPs7GwMHToUAOBwODBhwgS89dZbze4TmV3pLIkkdldMgUAAlm37H/Ibb7yBP/zhD/jss8/wzjvv4MEHH8TmzZtx/vnnY//+/fx5BoMBcrk85ljTJq2E9HaBgDdcpl0Di6UMLlc9fD4nRCJxeHZmIEQiSfs3Ij2OZYFquzqUBNyow5kerGhi2QAaG0/wMzP19YcRDPpizlEojFHBTBGUyq79Pu9uFNy0QyAQxL001BusXLkSY8eObTYDMnz4cOzYsSPm2I4dOzBs2LAWZ20AQKfToaioCGvWrIFEIsF5552H9PR0XHvttdi4cSOfbwMA48ePxzvvvIP09HRoNJp2xzl48GBIJBLs3r0bAweGpjGtViuOHTuGadOmdfj7jeTyBIPN86PGjRuHcePGYdmyZZg8eTLWr1+PSZMmoaCgoNm5LR0jpLfiOBYejyXcgPIUbLbT4dyZICQSBWQyHdTqHJqd6YU4Dqh3KEI7AZt1OGXWwB+Mfa+Rif1ICwcz+QYLsrT2hFQ0sWwQFksZH8zU1R1CIOCOfW6ZLiaYCeVg9Z2fo773rk06ZPTo0ViwYAFeeOGFmON//OMf8bOf/QyPPfYYrr32WuzatQtr1qzBX//61zbvN2PGDPzlL3/B1VdfDSA00zF8+HC88847ePHFF/nzFixYgKeeegrz58/HihUrMGDAAFRUVODDDz/Efffd1yzJWa1WY+HChbj33nthMBiQnp6Ohx9+GEKhMK5/SOnp6ZDL5fjss88wYMAAyGQyNDY24pVXXsG8efOQnZ2No0eP4vjx47jhhhs6fF9CeptAwBM1O1MOl6sefr8TAoEYMpmOZmd6MYubQWWjLlyirYPbH/v/SSIKIk3hgincoylXZ4VCGuhyMMNxLKzWSj6Yqa39AX6/M+YcqVSN9PRRyMgoQnp6ETSaAX0qmGmKgpsUtmLFCrzzzjsxx8aPH493330XDz30EB577DFkZWVhxYoV7S5fTZ8+Hc8//zyfWwOEAp7//e9/MccUCgW2bduG+++/H1deeSXsdjtycnJw0UUXtTqT8+yzz+LWW2/FpZdeCo1Gg/vuuw+nTp2CTNbxpnpisRgvvPACVqxYgYceeggXXHAB3nnnHRw5cgRvvvkmGhoakJWVhSVLluB3v/tdh+9LSLJxHAu32wynsxY222nYbKfg8djAcQFIJCqanenFnF4JKs1aPm/G5on9nSYSsNAr3EhXOZCjt2FwmgVKqa/LFU0cx8FuP8PnzNTW/hBuWnqORKKAyTSKn5nR6fJSqjpOwPWmrWt7gM1mg1arhdVqbfZm6/F4UFZWhvz8/LjeWEliOZ1O5OTk4JlnnsHNN9+c7OHEjX6OSFf5/e6Y3Bm3uwE+nwsikRgMowPDaGh2phfyBkShHk1mLSoatWhwKmMeF4CDTu6GUeVCltaOAmMjNHIvxF2saOI4Dk5nDT8zU1NzEB5PY8w5IhEDk2kEPzOj1w/utiamdvtZyOUGjByZ2FZJbb1/N0UzNyTp9u3bhyNHjmDixImwWq1YsWIFADSr6iIkVYVmZxqjZmdOw+OxguNYSKUqyGQGqNV9e5kgFQWCApy1avi9ZqptTSuaALXME6poUjtRkN4Ag8KTkIoml6s+Kpg5AJerLuZxoVACo/E8vjzbYBjarwJiCm5Ir/D000/j6NGjkEqlmDBhArZv3w6j0ZjsYRHSbfx+V3h2pjo8O9MIv98DoVAMmUwLnW4QhEL6Fd2bsCxQY1ehwqxDZaMWZ6waBNnYpRyF1Aej0okMtQODjWZkqF2Qiru+EazHY0ZNzUE+mHE4qmIeFwhESEsbFp6ZGQ2j8bx+vW8R/cshSTdu3Djs2bMn2cMgpFtxHAuXq4GvbLLbz8DrtYLjOEilKsjlRmg0ivZvRHoMxwENTjmfM3PKrIWvSUUTIw4gLRzMDDJYkaOzJaSiyeu1o67uB352xmqtjHlcIBBCrx8SFcwMh0Qi79qTphAKbgghpJv4fM6o2Zny8OyMCyIRE56dyafZmV7G6mbCOTOh2RmXP3b2QywMIk0ZamuQq7dioMEKhdQPYReDmY60NNDp8vlgxmQaCalU2fLNCAU3hBCSKCwbhNsdmZ2phN1+Fl6vLTw7o4ZCYYREQrMzvYnTJ8GpRi0f0FibVDQJBSwMCjdMUT2aVIy/y20NOtrSIFLNlJ4+CgzT/t5hJISCG0II6YLQ7EwNHI5qmM3l8HjM8PtdEIsZMIwOOl33VaWQ+PkCIpyyaPj9ZupbrGgK9WjKDPdo0iagoqmrLQ1IfCi4IYSQOLBsEC5XfXh2pgJ2ezW8XhsEAoRnZ0yU+9CLBFgBzlrVfDBTbVM3r2hiPHxbgwJTI9KU7i5XNIVaGhznk4Bbbmlg4quZ0tNH9/qWBn0JBTeEENIOr9cOl6sOdnsVLJbQ7Ewg4OFzZxQKmp3pLVgOqLWrwj2aQhVNATb2/41C4oNR5US62onBaWZkaJyQirqWBBxqafATH8zU1f2IQMATc45Mpo8JZvpaS4O+hIIbQghpgmUD/OyMxVIBh6Oa3+GVYTRQKNJpdqaX4Dig0SUPVzPpUGnWwhuIfWuTigKhHk0qB/LSLBigs0Mm7lpbg1BLg4pwMHMAtbU/ttjSIJIzk5FRRDtJ9yAKbkiLSktLUVxcDLPZDJ1O16FrHnnkEWzYsCGmszYhfYXXa2tS2WRGIOCFWCwLz86YaHaml7B5pKhs1PE7ATt9TMzjYmEQBoUL6epQRVOe3goF07WKpo61NFAiPX0kPzuj1aZWS4O+hIKbFLB27Vrce++9MJvNEItD/0sdDgf0ej2mTp2K0tJS/txI0HLixAkMGTKk1XtOmTIFVVVV0Gq1CR3rjBkzMHbsWDz//PMJvW9ngjHSvwWD/qjZmXI4nTXweu0QCASQStVQqTIhFlP7jN7A5RNHtTXQweKOnTUTCljo5W6YVC5k62wYYmzsckVTR1oaiMUymEwj+CRgSh7vPSi4SQHFxcVwOBz4/vvvMWnSJADA9u3bkZmZiW+//RYej4fvcbRlyxYMHDiwzcAGAKRSKTIzM7t97IT0JK/XBoejBg7HWVgslfB4zAgGfRCL5WAYLZTKDPpLuxfwBYQ4bdHwm+fVOVRNzghXNCmdyNQ4UWCqh07uhVjUtYomp7Mu3DW77ZYGGRlFfEsD2qeod6L/KymgsLAQWVlZKC0t5YOb0tJSzJ8/H1999RW++eYbvnN3ZIaDZVmsWrUKr7zyCqqrqzFs2DAsX74cV199dcx50TMhr776KlasWIGGhgbMnj0bF1xwAVasWAGLxRIznn/84x9Yvnw5zGYz5s6di1dffRVqtRqLFi3C1q1bsXXrVqxevRoAUFZWBq1Wi9///vf4/PPP4XA4MGDAADzwwAO48cYbe+T1I6krNDtTF56dKYPTWRuenRGBYTRQqbIhFjPt34h0q2Ckoik8M1NtU4HlYoNMFePl2xoUmBqRpnJD2sWKpkhLg8jsTMstDQr5vBmjsbBftzToSyi4aQ/HAX5Xcp5bokBHM96Ki4uxZcsW/OlPfwIQmqG57777EAwGsWXLFsyYMQNutxvffvstbrrpJpSUlOCf//wn1q5di6FDh2Lbtm349a9/DZPJhOnTpze7/44dO3Drrbdi1apVmDdvHr744gssX7682XknT57Ehg0bsHHjRpjNZlxzzTVYuXIlnnjiCaxevRrHjh3DqFGj+OaYJpMJd955Jw4dOoRNmzbBaDTixIkTcLvdXXjhSH/FcRyfO2O3n4XVWgG324Jg0AeJRA6ZTAelMpNmZ5IsVNGkRGW4R9NpS/OKJrnEB6MyVJ6db7QgU+PockWT12tDbe0PqK09gJqag7DZTsU8HmppUMAHMybTcFqa7KMouGmP3wU8mZ2c537gLNDB7bWLi4uxdOlSBAIBuN1u7Nu3D9OnT4ff78fatWsBALt27YLX68WMGTMwYsQIfPHFF5g8eTIAYPDgwfj666/x8ssvtxjc/OUvf8HcuXNxzz33AACGDRuGnTt3YuPGjTHnsSyLdevWQa1WAwB+85vf4Msvv8QTTzwBrVYLqVQKhUIRs+RVWVmJcePG4fzzzwcADBo0KL7XifRrwaAPTmddVO5MLXw+OwQCIRhGB7WaZmeSjeMAs0se7p4d6tHkCcR2qJaKAqG2BmoH8gxW5OptXa5o8vmcqKs7FA5mDsBiKUfrLQ2KYDKNoJYGKYKCmxQxY8YMOJ1O7N69G2azGcOGDeNnYW688UZ4PB6UlpZi8ODBcDgccLlcmDVrVsw9fD4fxo0b1+L9jx49iiuuuCLm2MSJE5sFN4MGDeIDGwDIyspCbW1tm2O/7bbbcNVVV2Hv3r24+OKLcfnll2PKlCkAgLlz52L79u0AgLy8PPz4448YOXIkKioqAAAXXHABNm3a1IFXiKSK0OyMFU5nLWy2s7Bay+HxWMGyfojFCshk2vD+ITQ7k0x2j/RcjyazFg5vSxVNbpjUDgzQ2ZBvsHS5oikQ8KCu7jA/M2M2U0uD/oqCm/ZIFKEZlGQ9dwcVFBRgwIAB2LJlC8xmMz/7kp2djdzcXOzcuRNbtmzBhRdeCIfDAQD45JNPkJOTE3MfhunaX7gSSexfYwKBACzb9rr43LlzUVFRgU8//RSbN2/GRRddhCVLluDpp5/Ga6+9xi9RRe796aefwu/3AwDkctprpD8Izc7UwumshdlcBperDj6fMzw7o4VGM4ByIZLM7Q9VNFWE+zSZXbG/vyIVTUaVC9laG4aYGqHuYkVTMOhDff1RPphpbGyppUEWv2leqKWBvtPPR/oOCm7aIxB0eGko2YqLi1FaWgqz2Yx7772XPz5t2jRs2rQJ3333HW677TaMGDECDMOgsrKyxSWolhQWFmL37t0xx5p+3RFSqRTBYLDZcZPJhIULF2LhwoW44IILcO+99+Lpp59uFnwBoRkckto4joPHYwnnzpyBxVIJr9cClg1CIlGAYbRQqbJpQ7Qk8gWFOGPRoLJRiwqzDrV2JRDT1oCDVhZqa5CpcaDA2AC9wtOliiaWDaCh4Rhfmt3QcKTFlgbRwQy1NOifKLhJIcXFxViyZAn8fn9M0DJ9+nT8/ve/h8/nQ3FxMdRqNe655x7cddddYFkWP//5z2G1WrFjxw5oNBosXLiw2b3vuOMOTJs2Dc8++ywuu+wyfPXVV9i0aVPcby6DBg3Ct99+i/LycqhUKhgMBjzyyCOYMGECRo4cCa/Xi40bN2L48OFdfj1I3xIIeMOb6NXAYimD290An88JoVAEqVQDjWYgRCJJ+zci3SLIClBlU/E9mqps6lYqmkJ5M0OMjTCpXV2qaGLZIMzmn/icmfr6wy22NIgEMxkZReFyfgp6+zsKblJIcXEx3G43zjvvPGRkZPDHp0+fDrvdzpeMA8Bjjz0Gk8mEkpIS/PTTT9DpdBg/fjweeOCBFu89depUrF27Fo8++igefPBBzJ49G3fddRfWrFkT1xjvueceLFy4ECNGjIDb7UZZWRmkUimWLVuG8vJyyOVyXHDBBXj77bc7/0KQPiE0O2MO586chtV6Cl6vlZ+dkcl0NDuTRBwH1DmU/DLTabMW/iYVTTKJP1zRZMegNAuytV2raOI4FhZLBb/MVFdHLQ1I5wg4juvarkd9jM1mg1arhdVqhUYTm0jm8XhQVlaG/Px8ftM70rrFixfjyJEjfMIvCaGfo9YFAp6o2ZlyuFz18PudEAjEkMm0YBgtzc4kCccBFreMTwCubGyjoknlwECDFQN0Niikna9o4jgONttpftO82tqD8PnsMedQS4O+x24/C7ncgJEjf5nQ+7b1/t0UzdyQDnv66acxa9YsKJVKbNq0CW+++Sb++te/JntYpBfjOBZu97nZGZvtFDweKzguCIlEBZlMR395J5HDKwktM5m1qGzUwd6kokkkZEM9mlRO5OhsyE8zQyn1Q9jJ2ILjODgc1THBjMdjjjmHWhqQRKDghnTYd999hz//+c+w2+0YPHgwXnjhBdxyyy3JHhbpZfx+N1yuOtjt1bBay+FyhXJnRCIxGEYHrTaPZmeSxOMX4ZQlNCtT0ahDY5OKJgE46BVumFROZGnsGGI0QyP3dqmiKdTS4ADfPdvlqo95XCSSwmgcHp6ZGU0tDUhC0E8Q6bB333032UMgvVBodqYxanbmdHh2hoVUqoJMpqfZmSTxB4U4a1WHlpoataixq8A1qWjSyLwwKZ3I0Dgw1NQAncILSReSgN1uc0ww43BUxzwuFIqRljaMz5lJSxtGZfwk4Si4IYTEze93wemsg8NRHa5saoTf74ZQGMqd0ekG0V/fScCyQLVdHUoCbtThjLV5RZNSGqpoytA4MNhoRrrKCam488FMpKVBZJmppZYGBkMBnzNjNFJLA9L96LdPC/pZjjVJsFT8+eE4Fi5XA5zOWlitp2C3n4HXawXHcZBKVZDL06DRdHzTSZIYHAfUOxShnYDNOpwya+APxv5al4n9oSRgjQODDBbkaO1gxJ2vaAq1NPiRD2YslrImZwig1+fzwYzJNBKSODYkJSQRKLiJEtkB1+Vy0c63pNNcrlCj1aa7Nfc1odmZ2hZmZ6Th2Zl8mp1JAoub4feaqTTr4PbH/pxJREGkKV0wqRwYqLdioN7apYqmUEuDQ3wSsNl8sllLA612YFQwMwoMo27lboT0jKT+ZnrppZfw0ksvoby8HAAwcuRIPPTQQ5g7d26L569btw433nhjzDGGYeDxeFo8P14ikQg6nY7vhaRQKChPgHQYx3FwuVyora2FTqeDSNS3KjxYNgi3OzI7Uwm7/Sy8XlvU7IyRZmeSwOmVhEqzzaGAxuaJXdIRCVjoFW6kqx3I0dqRn2aGivF1uqIp1NLgCB/MNDYeb9bSQK3O5oOZ9PRRkMmopQHpXZIa3AwYMAArV67E0KFDwXEc3nzzTcyfPx/79u3DyJEjW7xGo9Hg6NGj/NeJDj4i3arba/ZISGt0Ol1M1/PezOdz8rMzZnMZPB4z/H43xGIpGEZHZbhJ4A2IcCq8z0yFWYsGZ2z7FwE46MI9mrK0dhQYG6GReyHuZEVTMOhHY+NxPpiprz8ClvXHnBPd0iAjowgKhbHT3x8hPSGpwc1ll10W8/UTTzyBl156Cd98802rwY1AIOjWNw6BQICsrCykp6fzzRkJ6SiJRNKrZ2xYNgiXqz5qdqYKXq8NAAepVA2Fwkj5ET0sEBTgjFUT3jhPh2pb04omQCPzwKh0IkPtRIGpEQalu9MVTR1raWAI7wJMLQ1I39RrFsyDwSDee+89OJ1OTJ48udXzHA4H8vLywLIsxo8fjyeffLLVQAgAvF4vvF4v/7XNZuvQeEQiUa9+kyKko3w+R7gBZRUslnJ4PGYEAh6IRAxkMi0UCpqd6UksC9TYVagwh8qzz1g0CDapaFJIfeFgxoHBaWZkaFyQips3nO2IjrQ0YBgN0tNHUUsDkjKSHtwcPHgQkydPhsfjgUqlwkcffYQRI0a0eG5hYSH+9re/oaioCFarFU8//TSmTJmCH3/8EQMGDGjxmpKSEjz66KPd+S0Q0quwbICfnbFYKuBwVIdnZ0JvYkplBpXi9iCOAxqccj5n5pRZC1+TiiZGHEBaOJjJC7c16GxF07mWBpG9ZlpraTAqqqXBQGppQFJK0ntL+Xw+VFZWwmq14v3338drr72GrVu3thrgRPP7/Rg+fDiuv/56PPbYYy2e09LMTW5ubod6UxDSV3i99nDuTGh2xu02IxDwQiyWQSbTQipV0+xMD7K6GX6ZqcKshcsXu0mdWBgM92hyIldvxUCDFUqpv9PBjMNRFc6ZOdhKSwM5TKYRfN5MqNKNfh5I96DeUgCkUikKCgoAABMmTMDu3buxevVqvPzyy+1eK5FIMG7cOJw4caLVcxiGAcMwrT5OSF/EsgE4nXXh3JnQ7IzP5wDHAQyjhkqVSbMzPcjlE6MyvMxU0aiDtUlFk1DAwhBua5Cts2Fwmhkqxt/ptgbnWhqE9pppu6VBEQyGAirbJz0iGPQhEPC2f2I363U/7SzLxsy0tCUYDOLgwYO45JJLunlUhCSf12sL586chcVSAY/HjGDQB7FYDobRhpM+aWmhJ/gCIpyyaFAZ3mumztG8okkr98CkciJTY0eBqRHaLlQ0xd/SoJD6d5Fux3Ec/H4XfD47fD4nOC4IoVAChlFBrc5K6tiSGtwsW7YMc+fOxcCBA2G327F+/XqUlpbiv//9LwDghhtuQE5ODkpKSgAAK1aswKRJk1BQUACLxYKnnnoKFRUV1LyRpKRg0B+VO1MOp7MGXq8dAoEQUqkaKlU2xGKalewJAVaAs1Y1v3letU3drKJJzXhgVLmQoXZgiLERRlXnK5oiLQ0iAU1bLQ3S04tgMlFLA9L9WDYAn88Bn88Bvz+yWakCDKOBwTAUKlU65HIDZDJ90n83JTW4qa2txQ033ICqqipotVoUFRXhv//9L2bNmgUAqKyshDBqJyqz2YzFixejuroaer0eEyZMwM6dOzuUn0NIb8dxHD8743CEZmfcbguCQR8kEjlkMh3NzvQQlgNq7arwLsChiqYAG5ujopD4kBYOZvLTzMjUOMF0sqIp0tIgstdM6y0NipCRMZpaGpAeEQh44PXa4fM5wLJ+CIUiSCQqqFSZUKtzoFCkQS7XQybT9brfS0lPKO5p8SQkEdLdQrMzdeHZmTI4nbXweh0QCIRgGA0YRpv0v4D6A44DGl3ycDWTDpVmLbyB2L/9pKIAH8zk6S0YoLdDJu5cW4NIS4NIzkzLLQ3ywjkzo6mlAel2LBvkl5j8fic4joNYLINUqoRanQ2VKpOflZFKle3fsBv0qYRiQvqT0OyMFU5nLWy2s7DZQrMzLBuAWCyHTKaFUpnZ6/4KSkU2jxSVjaFApqJRC6cvNogUC4MwRCqadDbkGSxQMH4IOxHMRFoaRIKZhoZj4LjYWR61OqdJSwNdF747QtoWDPr4JaZAwA2BQAiJRAmZTIv09JFQKEzhYEbXJ/O3KLghpJsFgz6+sslsLoPLVQefzw6BQASG0UKjGQCRSNr+jUiXuP3icAJwqKLJ4o5tjisM92gyKZ3I1tkx2NgIdScrmiItDSLBTEstDZTK9KhgpggKRVqXvj9CWsNxHAIBN7/EFJ34q9PlQa3O5mdlGEaTEhs4UnBDSIJxHAePxxJT2eT1hmZnJBIlGEYDlYpmZ7qbLyDEaYuG3zwvVNEU/Uubg07ugVEVamsw1NgAndLTqYqmUEuDk3zOTF3dIQSDsVWfcrkhKpgZDZWqb/QfI31PKPHXGV5iiiT+ysOJv0OgVIYSf+VyfcomolNwQ0gCBAJePnfm3OyME0KhCFKpBhpNLs3OdLNgpKIp3HSyyqYG26StgYrx8m0NCkyNSFO5Ie1ERVOopUF5VDDzI/8mEhFqaXAumKGWBqS7BAKe8BKTHcGgP1xRqYJSmQ6NZgAUijTIZKHE3/6yeSMFN4R0Qmh2xhzOnTkDq/UUvF4LgsEApFIlGEYLlSqb3sy6EccBtXYl36PpdAsVTXKJD0alC+lqB/LTLMjSOiAVxd/WINTS4BQfzNTW/kAtDUhScBwbnpVxwO93hBN/GUgkKqSlFfKJv3K5IWmJv70BBTeEdFAg4IHTWQeHowYWSxlcrnr4fE6IROJw7szAPpl411dwHGB2yVERnpk5ZdbCE4h9vaWiQKitQbhHU67e1qmKptiWBgfCLQ0sMeeEWhqM5LtnU0sD0h2CQX94k7xQ4i8ghFSqhEymQXr6CCgURj5fhn7/nEPBDSGt4DgWbndkduY0bLbT8HqtfO6MTKajpYZuZvdI+QTgSrMWDm8LFU0KN0xqBwbobMjvQkWT01nL7wBcU3MQbje1NCA9K5L4G1liYtlAOPFXDa02FxrNAMhkesjlhpRJ/O0u9C+TkCih2ZnaqNmZBvj9TgiFYjCMjmZnupnbL8Yps4ZvOGl2xW5UJxSw0IXbGmRr7Rhiauh0RZPb3cg3mqypOQCns6WWBoV8MEMtDUiisWwwasdfJwAOEokCUqkaGRlFUCozIJeHgplUTfztLhTckH4tNDvTGDM74/HYwHEBSCQqmp3pZv5gpKIp1EG7xt68okkr98CodCFT40CBsQF6hQdiUfzBTKilwbllJpvtdMzjoZYGQ/lgxmgcThsokoQKBLxRib9eCAQiSKUqKBRp0GrHQi5P4/eWoSXOrqHghvQ7fr87ZnbG7W6Az+eCSCSBTKaFTpdHyw3dJMgKUG1T8ctMZ62tVTSF8maGGBthUrs6VdHk8zlQV/cjH8xYLOVNzhBArx8ctQswtTQgicNxbHjHX0e4fUEQYjEDqVQNg6EgvLdMaFZGIlHSH1AJRr/BScrjOBYuVwM/O2O3n4HHYwXHBSGVqiGTGaBWD6BfLt2A44A6h4Lfa+a0RQt/MPYvUpnEHy7PdmKQwYxsXecqmvx+N+rro1sa/NRGS4MimEwjqaUBSZhQ4m9kx18XAAFfOWk0FkKpTA/ny+hpW4geQMENSUl+vys8O1Mdnp1phN/vhlAoDc/ODKLZmW7AcYDFLUNloxYVZh1OmbVw+2PzVCThiqYMlQO5BhtydVYopPFXNAUCXjQ0HOGTgBsajrfR0iBU0UQtDUgixCb+OsKJvyJIpWpoNNlQqwdAoTBEJf7SlgA9jX67k5TAskE+d8ZqrYTdfhZerxUcx0EqVUEuN0KjoSWH7uDwSlAZ3mumolEHe5OKJpGQhUHhgknlRI7OhsFpZig7UdEUDPrR0HCMz5tpaDgClg3EnBNqaVAUDmaopQFJjFBTSWc4mIkk/soglWqQnj4SKlUmX8UkkcjbvR/pfhTckD7L53NGzc6Uh2dnXBCJGMhkOuh0gykprxt4/CKcsmj5YKaxSUWTQMBBL3fDpHIiU2PHEKMZWrk37oqmSEuDmpoDqKk5gPr6w620NCjiZ2aopQFJhHOJv45w4q8g/EeSARkZReG9ZSI7/tLbaG9E/1dInxGanWloMjtjC8/OqKFQGCkhtBv4g0KctapDScCNWtTYVeCaVDRpZKG2BpkaB4aYGpGmcMdd0RRqaVDGLzPV1R1qoaWBNmqZqQhqNe0CTbomOvHX73eCZYMQiSSQStXQ6/Oh0eTwm+RJpSr6eesjKLghvZrP5+BnZ8zmMng8Zvj9bojFDBiGZme6A8sC1XY1KhpDszNnrRoEm1Q0KaXhiiaNA0OMZqSrnJCK46toirQ0CCUAH0Bt7Y+ttDQ4F8yEWhrQmwvpvGDQD7/fAa83tOOvQCCAWCwHw2iRljYs3FQytMREib99FwU3pFdh2QBf2WS1VsBur4bXa4NAgPDsTDqtaScYxwH1TgWfBHzarIEvGPurQSb2h5KANaG2BjlaGxhxfBVNkZYGkWqmmpoD8HqtMedQSwOSSKHEX0+THX9De8toNFnhxN9IU0ktJf6mEApuSNJ5vfbw7ExVOHfGjGDQG86d0UKhoNmZRLO4mdAuwI1aVJp1LVQ0BUM9mlQO5OptGKi3dKqiKdTS4Fww43Y3xDweaWkQ6ZxNLQ1IV7SU+CsWy8AwaphMI6KaSuppCTvF0W8R0uNCszP1cDprYbFUwOEIzc4AAjCMGipVJm01nmBOryS0C3B4vxmbJ/b1FQlY6JVupKscyNHakZ9mhorxQRjnH7Jud0OTlgY1MY9TSwOSSIGAN2qJyQOhMNRUUi7XIyNjNBQKI7+3DAXN/Qv93yY9wuu1wemshd1+FhZLBTweMwIBHyQSORhGA4XCRLMzCeQNiHAq3D27wqxFg1MZ87gAHHQKN0xKFzK1dgwxNkLXiYomj8eK2tpzwYzdfib2ecItDSIzM9TSgHRW0x1/OS4IkUgKqVQFvT4/vONvaG8ZSvwlFNyQbhEM+qNmZ8rhdNbA67WHSyo1UKmy6U0ugQJBAc5Yz/VoqrY1rWgCNDIPvxNwgakRBqUbkjjbGvh8DtTW/sAHM1ZrRZMzQi0NIsGMyTSCpv9Jp7BsIKqpZKhqTiJRgGE04cRfE1/FRL9LSFMU3JCE8Xis4dyZs7BYKuF2m8GyPr4SQanMoIS9BGFZoMauQkV487wzluYVTQqpLxzMOJCfZkamJv6KJr/fhbq6w6itPYCamoOwWFpuaRAJZtLTR0EqVXX5+yP9j98fveOvD0KhGBKJCipVJtTqHCgUkaaSlPhL2kfBDem00OxMXXh2pgxOZ214dkYEhtFArabZmUThOKDBKedzZk6Ztc0qmhhxAEalE+nqUEXTAF38FU2BgBf19Uf4YKax8VizYEatzuFLs0MtDbSJ+BZJPxKd+Ov3O8FxocRfqVQJk2k4VKqMqL1llO3fkJAmKLghHcZxHLxeK5zOOtjtZ2G1lsPttiIYDOXOyGQ6KJWZ9FdVgljdDL/MVGHWwuWL3XNDLIxUNDkxQG9FnsEKpdQfVzBzrqVBKJhpuaVBBp8AnJ4+mloakLgFg76oppKe8PK0EjKZDunpo8I7/hoo8ZckDP0UkTYFgz44nXVRuTO18PkiszNamp1JIJdPHNWjSQurJ3Y/H6GAhUERamuQHe7RpI6zoinU0uBEuKXBQdTXH0Iw6Is5Ry5Pi1pmGg2VKiMR3x7pJziOi0n8ZdkARCIpGEYFnS4PanUOv0meVKqmxF/SLSi4ITE4joPHYwlXNoX2nfF6rWBZP8RiBWQyLVQqmp1JBF9AhFMWDSrDe83UOVqoaJJ7YFSF2hoUmBrirmhq2tKgtvZHBALumHPOtTQoQkZGEVSqLHrDIR3WPPGXg0SiBMNoYDAU8Im/crmetnggPYaCG4JAwMvnzpjNZXC56uDzOcM7eWqg0QygbcgTIMAKcNaq5jfPq7apm1U0qRkPjCoXMtXhHk1xVjQ1b2nwA3w+R8w5UqkK6emjwg0ni6DR5FIwQzosEPDA67XD73cgEPCFf0+ooVSmQ6PJhUIRKsdmGC1t70CSJu7gxu/3QyJpedOt+vp6GI3GLg+KdK/Q7Iw5at+ZSni9FgSDAUilSjCMFioVNSTsKpYDau2q8C7AoYqmABv7y14h8cGociE9XNGUpXFCKg52+Dk4joPdfjZmr5mWWhqkp4/ku2frdPk080Y6hONY+HxO+Hz2cOIvG078VSEtrTBqx18DJf6SXiXu4Oa6667D+++/3+yNr6amBhdddBF++OGHhA2OJE4g4IHTWQeHowYWSxlcrnr4fE6IRGIwjBYazUDaKbaLOA5odMnD1Uw6VJq18AZi/4lJ+YomJ/L0FgzQ2yETx9fWwOGo4QOZ2tqDrbQ0GBFeZhoNvb6A/oImHRKb+OsGIAwn/mqRnj6ST/yVyfT0+4L0anEHN5WVlbjlllvw+uuv88eqq6tRXFyMkSNHJnRwpPM4juVzZ6zWU7DZTodzZ4KQSBSQyXRQq3NodqaLbB4pKhtDgUxFoxZOX2xytVgYhCFS0aSzYZDBAiUTX0VTpKVBJJhpuaXBeXwScFraMHrjIe0KNZV0w+u18zv+CoViMIwaWu1AaDQ54dYFBjCMhn5XkD4l7uDm008/xbRp03D33Xfj2WefxdmzZ1FcXIwxY8bg7bff7o4xkg4Kzc7URs3ONMDvd0IgEEMm09HsTAK4fOJQWwOzFhWNOljczSua9Ao3TEonsrR2DDE1QhNnRZPHYwnvAhyqaGq5pcEwvnM2tTQgHRFK/A0tMQUCbnAcC4lEAalUDYNhMJTKDL6KiRJ/SV8Xd3BjMpnw+eef4+c//zkAYOPGjRg/fjzeeustCOPtske6hONYuN2NcDrrYLOdhs12Ch6PDRwXgESiotmZBPAFhDht0fCb54UqmqJfz3MVTRlqJ4YaG6BTeiCOo6KpYy0NhvAb55lMw6mlAWlXIOAJLzHZEQz6w3vLRBJ/c/imkjKZjpYtScrpVLVUbm4uNm/ejAsuuACzZs3CP/7xD3oD7SF+vztmdsbtboDP5wrnzuig1dLsTFcEWQGqrGpUhJtOVtnUYJu0NVAxXhiVTmSqHRhsbIRJHV9FU6ilwSE+mDGbfwIQGwxRSwMSj0hTyUgVE8dx4aaSahgMQ6OaSurpZ4n0Cx0KbvR6fYvBi8vlwn/+8x+kpZ3bsbSxsTFxoyNRszO14dmZ0/B4rOA4FlKpCjKZAWr1AAouO4njgFq7ku/RdLqFiia5xAejMlLRZEGW1gEmjoqmUEuDw3ww09h4vJWWBuc2zqOWBqQtwaCfn5UJBDwAwCf+mkznQalM5/Nl6I8d0h91KLh5/vnnu3kYJJrf7wrPzlSHZ2ca4fe7IRRKIZNpoNMNoi3KO4njALNbxu81c8qshScQ+8tfKgqE2hqEezTl6m1xVTSFWhoc5YOZhoajLbQ0yORzZjIyRkMup5YGpGWRxN9IMMOyocRfqVQFjWYANJoB/KxMKPGX0gMI6dA75MKFC7t7HP0ax7FwuRr4yia7/Qy8Xis4joNUqoJcboRGQzkWnWX3SEM9msJJwA5vCxVNCjdMagcGaO3ITzNDwfgh7GAww7JBNDYe54OZ+vrDbbQ0CJVnK5Xpifr2SIph2WDUjr9OhHb8DSX+RtphRMqxJRJ5u/cjpD+K+8//vXv3QiKRYPTo0QCAjz/+GG+88QZGjBiBRx55BFIp7WTbET6fM2p2pjw8O+OCSMRAJtNBpxtMSX6d5PaLccqs4RtOml2xgaFQwEIn9yBd5UCW1oEhxkaoZb4OtzVg2SAslnK+mqmuruWWBpFlJmppQNoSCHijEn+9EAhEkEqVUCjSoNGMCe8to4dMpqffCYR0UNzBze9+9zv86U9/wujRo/HTTz/h2muvxZVXXon33nsPLpcrriWsl156CS+99BLKy8sBACNHjsRDDz2EuXPntnrNe++9h+XLl6O8vBxDhw7FqlWrcMkll8T7bfQ4lg3C7Y7MzlTCbj8Lr9cWnp1RQ6EwUgVMJ/mDQpyxaFARnpmptTevaNLKPTAqQ20NCkwN0MdR0cRxHKzWSj6Yqa09GP6L+hxqaUA6IpL4G5mZ4Tg2nPirgsFQALU6i9/xVyJR0s8QIZ0Ud3Bz7NgxjB07FkAo0Jg+fTrWr1+PHTt24LrrrosruBkwYABWrlyJoUOHguM4vPnmm5g/fz727dvX4oaAO3fuxPXXX4+SkhJceumlWL9+PS6//HLs3bsXo0aNivdb6XY+n4OfnTGby+HxmOH3uyAWM2AYmp3prCArQLVNhYrw5nlnra1VNIXyZgYbzUhXOyHtYEXTuZYG54IZamlAOuNc4q8DgYALgAASiQIMo4XROAwKRTqfL0P92whJHAHHcR3fkAOARqPBnj17MHToUMyaNQuXXnop7rzzTlRWVqKwsBBut7v9m7TBYDDgqaeews0339zssWuvvRZOpxMbN27kj02aNAljx47F2rVrO3R/m80GrVYLq9UKjUbTpbE2xXEcnM7a8OxMBez2ani9NggEgFSqBsNoaY28EzgOqHMo+V2AT1u08Adjg0KZxA+jMrTXTJ7BjBydA1JRsMNJwKGWBgf47tlud2zVH7U0IO05l/jr5PeWEYlCib8KhRFqdQ4UirSoHX8pGCYkHvG8f8c9c3P++efj8ccfx8yZM7F161a89NJLAICysjJkZGR0bsQAgsEg3nvvPTidTkyePLnFc3bt2oW777475tjs2bOxYcOGTj9vIjkcVTh+/FN4PJZw7owWCgXNzsSL4wCLW4bKRi0qzDqcMmvh9sdWNElEgdDMjMqBXL0NuXorFNKOVzS5XA0xwYzTWRvzuFAohtE4nM+ZMRiGUkkticGyQfj9zvDMTCTxVwapVA2TaURMU0n6o4aQnhV3cPP8889jwYIF2LBhA/7f//t/KCgoAAC8//77mDJlStwDOHjwICZPngyPxwOVSoWPPvoII0aMaPHc6urqZgFURkYGqqurW72/1+uF1+vlv7bZbHGPsaNYNgiv1waDYSj9VRYnh1eCyvBeMxWNWti9sdu/i4QsDAoXTConcrR2DDaaoWR8Ha5oar+lgQhpaUP5Zaa0tPOopQGJcS7x1xFO/BWEqxkNyMgo4mdlQjv+0lYNhCRT3P8Ci4qKcPDgwWbHn3rqKYhE8c9QFBYWYv/+/bBarXj//fexcOFCbN26tdUAJ14lJSV49NFHE3IvkjgevwinLNpwMKNDY5OKJoGAg17uhlHlRJYmVNGklXs7XNEUaWkQaTbZtKWBQCCEXj+Yn5kxGkfQX9eEF5346/c7wbIBfsdfvT4/pqmkVKqixF9CepmE/Xkhk3Wu0ZpUKuVnfyZMmIDdu3dj9erVePnll5udm5mZiZqa2I7INTU1yMzMbPX+y5Yti1nKstlsyM3N7dRYSef5g0KctapDScCNWtTYVeCaVDRpZOG2BhoHhpgakaZwQyzqWDATaWkQCWZaammg0w3iZ2ZMppG0DT3hBYN++P2h5SW/3wUAfOJvWtowKJXpfFNJSvwlpPfrUHBjMBhw7NgxGI3GVlsxRHS1/QLLsjHLSNEmT56ML7/8EkuXLuWPbd68udUcHQBgGAYMQ8sLPY1lgWq7GhWNWlQ26nCmhYompTRc0aRxYHCaGRmajlc0RVoaRIKZlloaaDQD+GAmPX00GCaxCeSk7/L73fwSE8v6+B1/VapMaDQDoFCkhZtKammJmZA+qEPBzXPPPQe1Wg0gsa0Yli1bhrlz52LgwIGw2+1Yv349SktL8d///hcAcMMNNyAnJwclJSUAgDvvvBPTp0/HM888g1/84hd4++238f333+OVV15J2JhI53AcUO9QhCqazDqcNmvgC8b+eMnEfj6YyTNYkaO1gRF3rKIp0tIgEsy01NJApcpEenpRVEsDQyK/RdJHRSf+hpaYWEgkcjCMGibT8KjEXz3tNUVIioir/UIgEIBAIMDs2bO7VBkVUVtbixtuuAFVVVXQarUoKirCf//7X8yaNQsAUFlZCaHw3F9NU6ZMwfr16/Hggw/igQcewNChQ7Fhw4ZeucdNf2BxM3yPpkqzroWKpiDSFKG9ZnL1NgzUWzpc0cSyATQ2nuCDmZZbGhiRkTGab2ugVJoS+e2RPioY9MHns8PncyIQ8IQTf5WQy/XIyBgNhcIYzpfRU+IvISkq7n1uFAoFDh8+jLy8vO4aU7fqzn1urNZTOHToPej1Q1JyKtvplYR7NIXyZqyeJhVNAhZ6hQsmdaiiKT/NDLWsYxVNoZYGZXwwU1d3qFlLA5lMxycAh3rsUEuD/o7juHDibyiY4bgghEIJGEYFpTIdanUOPysjlarp54WQPqxb97mZOHEi9u3b12eDG9Jx3oAIp80aVISDmXqnMuZxATjoFG6YlC5kauwYYmqEroMVTRzHwmqt5IOZ2tofWmhpoEZ6+ig+mKGWBoRlA1FNJaMTfzUwGIZCpUrnm0pSKT8h/Vfcwc3tt9+OP/7xjzh9+jQmTJgApTL2Da+oqChhgyM9KxAU4Kw11KOpslGLapu6SUUToJF5+J2Ah5gaYVS6OlTRFGppcCYqmDkIrzd2zyGJRAGTaSSfM0MtDUgg4IHXaw8n/vohFIogkaigVGbwib+hppI6+lkhhPDiDm6uu+46AMAf/vAH/phAIADHcRAIBAgGg4kbHelWLAvU2FXhtgY6nLFqEGRj3yAUUl84mHEgP82MTI0TUnHHKpocjmo+mKmpOQiPp2lLAwYm0wh+ZkavH0K7OfdjocTf0BKT3+8Ex3EQixlIpSqYTOfxib8ymR5SqbL9GxJC+q24g5uysrLuGAfpARwHNDjlqDSHkoBPmbXNKpqYSEWT2oE8vRUD9B2vaHK56vl2BjU1B+FyNW1pIIHReB61NCAAIom/kaaSbggEQkgkynBu1UgoFCZ+x1/6OSGExCPu4IZybfoWq5sJJQE36lBh1sLli92ATCwMIk3pQrrKiQE6G/LSLFBK/R0KZkItDQ7ynbOppQFpzbnE31AwE534q9PlQa3O5vswUeIvIaSrqA4yxbh84pgeTVZPbEsBoYCFQeGGSeVEttaOwcZGqBkfhB1IV/B67air+4GfnbFaK2MeP9fSIBTMUEuD/iuU+OsMLzFFJ/6qYTAMCe/4G6piEos7t7s5IYS0hoKbPs4XEOGURYPKRh0qzVrUOZpXNGnlbphULmRoHCgwNUDfwYqmUEuDH/lgxmwuA7U0IC0JBDzhWRk7AgEfhEIRpNJQObZGkwuFIjQrwzBayqsihHQ7Cm76mAArwFmrmt88r6WKJjXjgVHlQoY61HDSqHJD0oG2BqGWBoeiWhqcaKGlQW7UXjOjqKVBP8RxbHhWxgG/3wGOYyEWyyCRqJCWVhi146+BEn8JIUlBwU0vx3JArV0V3gVYizMWDQJs7F++oYqmUBJwfpoFWRoHpOL2q9ZiWxocQEPDsVZbGkT6M1FLg/4nGPRHzcq4AQghlSohk2mQnj4CCoWRr2KixF9CSG/QqeDGYrHg/fffx8mTJ3HvvffCYDBg7969yMjIQE5OTqLH2K9wHNDokqMy3NKg0qyFN9C0oimANKUT6WonBuotyNXbIRO339agIy0NFApjOJgJlWdTS4P+heM4BAJuPphh2UA48VcNrTY3asdfAxhGQ4m/hJBeKe7g5sCBA5g5cya0Wi3Ky8uxePFiGAwGfPjhh6isrMTf//737hhnSrN5pHzOTEWjFk5fbEVRpKLJFK5oGmSwQMm0X9HUsZYGen7TvPT0IqhUmfSG1Y+wbDBqx18nAA4SiQJSqRoZGUV84q9MpqfkcEJInxF3cHP33Xdj0aJF+POf/8x3CgeASy65BL/61a8SOrhU5faLwzMzoc3zLO7mFU16hRsmpRNZ2lBbA00HKpo63tJgdFRLgwEUzPQjgYA33IfJgWDQB4FACKlUBYUiDRrNmKglJh0l/hJC+qy4g5vdu3fj5ZdfbnY8JycH1dXVCRlUqvEFhDht0fCb54UqmqIDCg46uQdGlROZageGmBqhV3ggbqeiqeMtDUbxMzM6XR5tU99PcBwLv98Frze04y/LBsM7/qphMBSE95bRQy43QCJRUpBLCEkZcQc3DMPAZrM1O37s2DGYTJSfAQBBVoBq27keTVU2NVguNqBQMV4YlaFgZrCxESZ1+xVNHMfB6ayJo6VBEfT6wfQXeD/RPPFXAKlUCYbRwmQ6D0plOmQyPeRyPUQiabv3I4SQviru4GbevHlYsWIF3n33XQChvlKVlZW4//77cdVVVyV8gH1Fgz2Irw+Jse/ETNQ6s+BvUtEklzSpaNI6wHSgoinU0iASzByAy1UX83ikpUFkmYlaGvQPsYm/jnDirwhSqRoaTQ7U6gFRe8toaLaOENKvCDiOa383tyhWqxVXX301vv/+e9jtdmRnZ6O6uhqTJ0/Gp59+2qxLeG9js9mg1WphtVqh0SRuj5ZDp/x47j92/mupKACj0gWT2oE8gw25emuHKpo8HjPfzqCm5gAcjqqYx0MtDYbxwYzReB79Fd4PhJpKOsNLTC6EEn9lkEo1UKuzoFJlhmdlDJT4SwhJSfG8f8c9c6PVarF582Z8/fXXOHDgABwOB8aPH4+ZM2d2esCpYEimGOdlB8B5jyPPJEJ+mgWqDlQ0eb021Nb+wDebtNlOxTweamkwhA9mTKYRtF19PxBK/HWEE3+9EAgELST+6sOJv7RdFSGERIt75ubUqVPIzc3trvF0u+6auQEAq/UUDh16D3r9kFaXAfx+F2prfwwHMwdgsZSjeUuDfD5nxmQaQbu8prhI4m+kHJtlgxCJJJBK1VCpMvimkjKZHlKpihJ/CSH9UrfO3AwaNAg///nP8etf/xpXX3019Hp9pwfaHwQCHtTVHeZnZszmllsaRKqZqKVB6gsG/fD7HfB6HQgE3BAIBBCL5WAYLYzGYVAo0vkqJlpyJISQ+MUd3Hz//fdYv349VqxYgTvuuANz5szBr3/9a1x22WVgGKb9G6S4YNDP7zFTU3MQjY0ttTTIitlrRi6nADFVhRJ/PU12/A01ldRosqFW50ChSINMpodMpqXEX0IISYC4l6UiOI5DaWkp1q9fjw8++AAsy+LKK6/E3/72t0SPMaESvSwVefnOnPkWn3/+R5w5812zYEahMPGBDLU0SG2RxN/IEhPHcRCLZWAYNVSqrKimknpIJIpkD5cQQvqMeN6/Ox3cRNu7dy9uvvlmHDhwAMFg++XNyZTI4Cby0gkEAlRV7cUrr0wAEGppEAlmQlvYZ1CeRIoKBLz8EtO5xF8lZDI9NJoBTXb8pcRfQgjprG7NuYk4ffo01q9fj/Xr1+OHH37A5MmT8eKLL3b2dn1SdMCSkTEGF174BLxeK3JyJtPGeSkoOvHX53OA49hw4q8Ken0+n/grl+shlaopoCWEkCSJO7h5+eWXsX79euzYsQPnnXceFixYgI8//hh5eXndMb4+QygUoajoNzh06D16U0sRLBsI92FyhveWCbWzYBgN0tKGQak08VVMYjHlmxFCSG8Rd3Dz+OOP4/rrr8cLL7yAMWPGdMeYCEkKvz96x18fhEIxJBIVVKpMPvE3FMxQ4i8hhPRmcQc3lZWVNDNB+rzWEn+lUiVMpuFQqTLCS0wGSvwlhJA+pkPBzYEDBzBq1CgIhUIcPHiwzXOLiooSMjBCEikY9PFLTIGABwKBABKJEjKZDunpo2ISf6k3FyGE9G0dCm7Gjh2L6upqpKenY+zYsRAIBIgusop8LRAIen21FEl9HMeFE39DwQzLBiASScEwKuh0eVCrc/hN8ijxlxBCUk+HgpuysjKYTCb+c0J6k1DiryO8xBRpKqkEw2hgMAzlE3/lcj315SKEkH6gQ8FNdCVURUUFpkyZArE49tJAIICdO3f2+6op0v0CAU+4O7YDwaAfAoEQUqkaSmVGeG8ZA7/ERIm/hBDS/8SdUFxcXIyqqiqkp6fHHLdarSguLqZlKZJQHMfC53PC57OHE3/ZcOKvCmlphVE7/hqowSghhBAAnQhuIrk1TTU0NECppDcX0jWhxN/QElMg4AYgDO/4q0V6+kgoFCbI5XrIZHpK/CWEENKiDgc3V155JYBQ8vCiRYtimmQGg0EcOHAAU6ZMSfwIScoKNZV0w+u1h3f8DUAolIBh1NBqc6HRDIBMFkr8ZRgNJf4SQgjpkA4HN1qtFkDoDUmtVkMul/OPSaVSTJo0CYsXL078CEnKCCX+hpaYAgE3OI6DRCKHVKqGwTAYSmUGX8VEib+EEEI6q8PBzRtvvAEAGDRoEO655x5agiLtCgQ84SUmezjxVxBO/DXxTSVlMn24qST14iKEEJIYcefcPPzww90xDtLHRZpKRqqYQjv+MpBIVDAYhkGtzopqKqlK9nAJIYSksE51BX///ffx7rvvorKyEj6fL+axvXv3JmRgpHcLBv38rEwg4AEggFSqgEymhck0PKapJCX+EkII6UlxBzcvvPAC/t//+39YtGgRPv74Y9x44404efIkdu/ejSVLlnTHGEmSRRJ/I8EMywYhFIohlaqg1eaGd/w1hBN/1bS3DCGEkKSKO7j561//ildeeQXXX3891q1bh/vuuw+DBw/GQw89hMbGxrjuVVJSgg8//BBHjhyBXC7HlClTsGrVKhQWFrZ6zbp163DjjTfGHGMYBh6PJ95vhbSCZYNRO/46AQASiQxSqQYZGUVQKtP5WRmJRN7O3QghhJCe1amu4JGSb7lcDrvdDgD4zW9+g0mTJmHNmjUdvtfWrVuxZMkS/OxnP0MgEMADDzyAiy++GIcOHWozYVmj0eDo0aP811Qi3DWBgDcq8dcLgUAEqVQJhSINWu1YyOVp/I6/lPhLCCGkt4s7uMnMzERjYyPy8vIwcOBAfPPNNxgzZgzKyspimml2xGeffRbz9bp165Ceno49e/Zg2rRprV4nEAiQmZkZ79AJziX+RmZmOI6FSCSFVKqCwVAQlfhrgESipMCREEJInxN3cHPhhRfi3//+N8aNG4cbb7wRd911F95//318//33/EZ/nWW1WgEABoOhzfMcDgfy8vLAsizGjx+PJ598EiNHjuzSc6eqc4m/DgQCLoQSf5WQSjUwGodBoUjnq5hEImmyh0sIIYR0mYCLc7qFZVmwLMs3znz77bexc+dODB06FL/73e8glXbuDZJlWcybNw8WiwVff/11q+ft2rULx48fR1FREaxWK55++mls27YNP/74IwYMGNDsfK/XC6/Xy39ts9mQm5sLq9UKjUbTqbG2xmo9hUOH3oNePyQpSbXnEn+d/N4yIlEo8VehMEKtzoFCkRa14y8l/hJCCOkbbDYbtFpth96/4w5uusttt92GTZs24euvv24xSGmN3+/H8OHDcf311+Oxxx5r9vgjjzyCRx99tNnxVAhuWDYIv98ZnplxAuDCib9qqNXZ4R1/Q7MyEomi28dDCCGEdJd4gpsOLUsdOHCgw09eVFTU4XMjfv/732Pjxo3Ytm1bXIENAEgkEowbNw4nTpxo8fFly5bh7rvv5r+OzNz0RecSfx3hxF8BpFIV5HIDMjKK+FmZUOJvp7YwIoQQQvq8Dr0Djh07FgKBoN2EYYFAgGAw2OEn5zgOd9xxBz766COUlpYiPz+/w9dGBINBHDx4EJdcckmLjzMME9Pks6+ITvz1+51g2UA48VcNvT4fGk0O31RSKlVR4i8hhBAS1qHgpqysrFuefMmSJVi/fj0+/vhjqNVqVFdXAwg16Yw05rzhhhuQk5ODkpISAMCKFSswadIkFBQUwGKx4KmnnkJFRQVuueWWbhljTwkG/fD7Q8tLfr8LACCRKMAwWqSlDQvvLaOHTKaHWNz3gjVCCCGkp3QouMnLy+uWJ3/ppZcAADNmzIg5/sYbb2DRokUAQvvqCIXn8lfMZjMWL16M6upq6PV6TJgwATt37sSIESO6ZYzdxe9380tMLOuHUCiCVKqCWp0JtXoAX44tk2kp8ZcQQgiJQ6cSiv/xj39g7dq1KCsrw65du5CXl4fnn38e+fn5mD9/fneMM2HiSUiKV2sJxdGJv6ElJhYSiRwMo4ZKlQWVKpMSfwkhhKSGYADw2gBF29u6xCvhCcXRXnrpJTz00ENYunQpnnjiCT7HRqfT4fnnn+/1wU1PiF5iCgQ84cRfJeRyPdLTR0GpNIXzZfSU+EsIIaRv8jmB+mNA3bHQf+uPhj5v/AkYciGw4N2kDS3ud9a//OUvePXVV3H55Zdj5cqV/PHzzz8f99xzT0IH1xeJxTLY7WfBMCrodIOgVmfzszJSqZoSfwkhhPQdHAc460OBCx/IhIMY2+nWr7NU9twYWxB3cFNWVoZx48Y1O84wDJxOZ0IG1VcplekYPHgmGEZDib+EEEL6DpYFLBVA/fFw8BIOZuqPAW5z69cpjICpEDAOBYyFgGkYYBwGaOLb1iXR4g5u8vPzsX///mZJxp999hmGDx+esIH1RWIxg7S0YckeBiGEENIyvwdoOHEucKk7GgpoGo4DAU8rFwkA3cBwEDPs3IepMOF5NYkSd3Bz9913Y8mSJfB4POA4Dt999x3+9a9/oaSkBK+99lp3jJEQQggh8XCbm+fC1B8Lzc5wbMvXiBggreDc7EskgDEMAaR9q9gl7uDmlltugVwux4MPPgiXy4Vf/epXyM7OxurVq3Hdddd1xxgJIYQQ0hTHAbYzsbkw9cdDszHO2tavk2ljl5Ain+vyAKGo58bfjeIKbgKBANavX4/Zs2djwYIFcLlccDgcSE9P767xEUIIIf1b0B+qQOKXkSIzMscBn6P16zQ553JhjEPDy0qFgCodSPHilriCG7FYjFtvvRWHDx8GACgUCigUfWuqihBCCOmVvPZzQUt0ENP4E8AGWr5GKAYMg2OXkYzDQsEMo+7Z8fcicS9LTZw4Efv27eu2XYsJIYSQlMVxgKO2eS5M/bHQElNrpKqWZ2EM+YBI0nPj7yPiDm5uv/12/PGPf8Tp06cxYcIEKJXKmMc70xWcEEIISSlsMJS8W9c0iDkKeKytX6dMD8/ANEnq1eSk/FJSIsXdfiG6zxN/k3DH8Hi7gidDd7ZfIIQQ0s/43aHS6uhlpLpjoWNBbysXCQB9XlRSb9RSUi8tre4NurX9Qnd1CCeEEEJ6LVdj84TeuqPhnXhbmSMQy4C0oVHLSOFlpbQCQCLr0eH3N3EHN5RrQwghJCWxbKilQEul1a761q+T6c4l8kZyYYxDQxvfpUhpdV9DXRsJIYT0LwEf0HiySRATrlLyu1q/TpvbvM2AsRBQGikfppeh4IYQQkhq8ljDvZKaLCc1lgFcK/mhQgmQNiQqiAnPwqQNBRhVz46fdBoFN4QQQvoujgPs1U16JYWXk+xVrV8nVUflwkRVJekHUWl1CqDghhBCSO8XDADm8nPl1PxGd8cBbxul1arM2H1hIp+rs2gpKYVRcEMIIaT38DnDpdVRuTB1x0I5MkFfy9cIhIA+P2p/mKjSarmuR4dPeoe4gxu9Xg9BC9GuQCCATCZDQUEBFi1ahBtvvDEhAySEEJKCnA3hze2a7A9jrWz9GrEcMBbE5sIYC0M5MmKm58ZOer24g5uHHnoITzzxBObOnYuJEycCAL777jt89tlnWLJkCcrKynDbbbchEAhg8eLFCR8wIYSQPoJlQ8FK015JdUcBd2Pr18kNscFLJC9Gmwu0sJEsIU3FHdx8/fXXePzxx3HrrbfGHH/55Zfx+eef44MPPkBRURFeeOEFCm4IIaQ/CHijduk9HlVafQIIuFu5SADocs8tIcWUVqf16PBJ6om7/YJKpcL+/ftRUFAQc/zEiRMYO3YsHA4HTp48iaKiIjidzoQONhGo/QIhhHSS23IueImeiTGXAxzb8jUiaWhH3qZdq9MKAKmiJ0dP+rhubb9gMBjwn//8B3fddVfM8f/85z8wGEI9MZxOJ9Tq/ttqnRBC+iyOA2xnm5RWhz931LR+HaNt3uzROAzQ5QEiql0hPSvun7jly5fjtttuw5YtW/icm927d+PTTz/F2rVrAQCbN2/G9OnTEztSQgghiRP0hzazi5RW812rjwM+e+vXqbNbDmJUGVRaTXqNuJelAGDHjh1Ys2YNjh49CgAoLCzEHXfcgSlTpiR8gIlGy1KEkH7F6zgXtPDLScdDpdVsoOVrBCLAMDi2tNo0LLRLr4x+b5LkiOf9u1PBTV9GwQ0hJOVwHOCsi92dNxLE2E63fp1Eca4iKTqQMQwGxNKeGz8hHdCtOTcAEAwGsWHDBhw+fBgAMHLkSMybNw8iEXU/JYSQbsMGAUtF1BJS1HKSx9L6dQpjk67V4YBGk0Ol1SQlxR3cnDhxApdccgnOnDmDwsJCAEBJSQlyc3PxySefYMiQIQkfJCGE9Ct+d6i0Orprdd2x0LGgt5WLBIA+r0kuTDiQURh6dPiEJFvcy1KXXHIJOI7DW2+9xVdHNTQ04Ne//jWEQiE++eSTbhlootCyFCGk13A1tlJaXQGglV/NIiY88zIUMb2S0goAibxHh09IT+rWZamtW7fim2++4QMbAEhLS8PKlSsxderU+EdLCCGpjOMA25nmbQbqj4byZFoj055L5I1uN6DLA4SUAkBIW+IObhiGgd3evEzQ4XBAKqUENEJIPxXwAY0/tVxa7W9jQ1PNgNjS6siSktJEpdWEdFLcwc2ll16K3/72t3j99df5fW6+/fZb3HrrrZg3b17CB0gIIb2KxxbbYiASxDT+BHDBlq8RigHDkHNLSNGl1YyqZ8dPSD8Qd3DzwgsvYOHChZg8eTIkEgkAIBAIYN68eVi9enXCB0gIIT2O40K78TZt9lh/DLBXtX6dVBXV7HHYuRJrQz4gkvTc+Anp5+IObnQ6HT7++GMcP34cR44cAQAMHz68Wa8pQgjp9YKBcGl10yDmOOC1tn6dKqP5Dr3GYYAmm5aSCOkFOt3wY+jQoRg6dGgix0IIId3D5wIajjffH6bxJBD0tXyNQAjoB7UQxAwF5PoeHT4hJD4dCm7uvvvuDt/w2Wef7fRgCCGkS5wNTXJhwkGMtbL1a8RywBjpWl14LrnXMASQyHpu7ISQhOlQcLNv374O3UxA07GEkO7GsqGWApHgJTqQcTW0fp3cELs7b+Rz7UDapZeQFNOh4GbLli3dPQ5CCIkV8AINJ5sn9DacAPyu1q/TDoyqSopaTlIae27shJCk6nTOTSKUlJTgww8/xJEjRyCXyzFlyhSsWrWKb+vQmvfeew/Lly9HeXk5hg4dilWrVuGSSy7poVETQhLKY215FsZcDnBsy9cIJUDakOYJvcahgFTZo8MnhPQ+SQ1utm7diiVLluBnP/sZAoEAHnjgAVx88cU4dOgQlMqWf0Ht3LkT119/PUpKSnDppZdi/fr1uPzyy7F3716MGjWqh78DQkiHcFyohDpSiRTdbsBR0/p1jKZ5s0dTYWiXXlFSf30RQnqxuHtLdae6ujqkp6dj69atmDZtWovnXHvttXA6ndi4cSN/bNKkSRg7dizWrl3b7nNQbylCulEwAJjLWi6t9jXf2Zynzmq5tFqdSaXVhBAA3dxbqjtZraF9JaL7VjW1a9euZtVbs2fPxoYNG1o83+v1wus910XXZrN1faCE9Hc+57nWAnVHo0qrfwJYf8vXCISAPr/lpF6ZtmfHTwhJab0muGFZFkuXLsXUqVPbXF6qrq5GRkZGzLGMjAxUV1e3eH5JSQkeffTRhI6VkH6B4wBnffNcmPrjgPVU69dJFKEO1ZE2A5HkXsNgQMz03PgJIf1WrwlulixZgh9++AFff/11Qu+7bNmymJkem82G3NzchD4HIX0aGwQslc1zYeqPAW5z69cpjOFlpKg2A6ZhoUaQVFpNCEmiXhHc/P73v8fGjRuxbds2DBgwoM1zMzMzUVMTm4BYU1ODzMzMFs9nGAYMQ38tEgK/J1RGHZl9iS6tDnhauUgA6AY2z4UxFQKK1pePCSEkmZIa3HAchzvuuAMfffQRSktLkZ+f3+41kydPxpdffomlS5fyxzZv3ozJkyd340gJ6UPc5lZKqysAtFI/IJKGlpKaBjFpBYBU0aPDJ4SQrkpqcLNkyRKsX78eH3/8MdRqNZ83o9VqIZfLAQA33HADcnJyUFJSAgC48847MX36dDzzzDP4xS9+gbfffhvff/89XnnllaR9H4T0OI4DbGealFaHgxhnXevXybSxS0iRnBj9IEAo6rHhE0JId0pqcPPSSy8BAGbMmBFz/I033sCiRYsAAJWVlRBGrd9PmTIF69evx4MPPogHHngAQ4cOxYYNG2iPG5Kagv5QBVL0Dr2RgMbvbP06TU7UEtKwc5VJShOVVhNCUl6v2uemJ9A+N6RX8trDMzDHYhN6G38C2EDL1wjFUaXVTbpWM+qeHT8hhHSzPrvPDSEpjeMAR204cIksI4U/bGdav06ibLlXkj4fEEt7bvyEENJHUHBDSKKxwVBfpJhcmHBA47G2fp3SFF4+GhZblaTJoaUkQgiJAwU3hHSW331uKSk6F6bhBBD0tnKRANDnNQliwkm9VFpNCCEJQcENIe1xNTbPhak7Gtr4rtXSaibcYqBJUm/aEEAi79HhE0JIf0PBDSEAwLKhvJemuTB1RwFXfevXyXQtJ/Tq8qi0mhBCkoSCG9K/BHzh0uomuTD1J9ourdbmRjV7jFpOUhopH4YQQnoZCm5IavLYziX0Ru/S21gGcMGWrxFKQstGTTtWpw0FGFXPjp8QQkinUXBD+i6OA+zVTRJ6w5/bq1q/TqpuXpFkHBbapVck6bHhE0II6R4U3JDeLxgIl1Y33R/mOOBto7RalRm7hBT5XJ1FS0mEEJLCKLghvYfPGV5KiuwPEy6tbjwJBH0tXyMQhjazi+mVFE7qlet6dPiEEEJ6BwpuSM9z1jcvq64/DlgrW79GLAeMBbG5MJHSajHTc2MnhBDS61FwQ7oHy4aClfrj4eDl6LnP3Y2tX6dIO7epnanw3OfaXCCqgSohhBDSGgpuSNcEvKEdeaMrkuqPhUqrA+7Wr9MNbJILE15OUqb13NgJIYSkJApuSMe4Lc2XkeqPhhJ9Obbla0RSIK2g5dJqqaInR08IIaQfoeCGnMNxgO1sy6XVjprWr2O0saXVkfJqXR4goh8xQgghPYveefqjoD8049I0F6b+OOCzt36dOjsqFyYqiFFlUGk1IYSQXoOCm1Tmc0blwkTtEdP4E8D6W75GIAIM+S13rZZpenb8hBBCSCdQcNPXcVyotDp6X5hIEGM73fp1EsW5XJjoPWIMgwGxtOfGTwghhCQYBTd9BRsELJXNc2HqjgIeS+vXKYxNlpHCQYwmh0qrCSGEpCQKbnobv7tJaXX4o+EEEPC0cpEA0Oc175VkHAYoDD06fEIIISTZKLhJFrf53L4wMaXVFQC4lq8RMeGlpKGI2egurQCQyHt0+IQQQkhvRcFNd+I4wHamyTJSOKBx1rV+nUzbfHM707BQabVQ1HPjJ4QQQvogCm4SxdkAVOxoXlrtd7Z+jWZAy6XVShOVVhNCCCGdRMFNotQcBN79TfPjQnGoAonPhYk0fRwGMKqeHychhBCS4ii4SRRjIZA9rvn+MIZ8QCRJ9ugIIYSQfoOCm0TRZAG/LU32KAghhJB+jzY6IYQQQkhKoeCGEEIIISmFghtCCCGEpBQKbgghhBCSUii4IYQQQkhKoeCGEEIIISmFghtCCCGEpBQKbgghhBCSUii4IYQQQkhKoeCGEEIIISmFghtCCCGEpBQKbgghhBCSUpIa3Gzbtg2XXXYZsrOzIRAIsGHDhjbPLy0thUAgaPZRXV3dMwMmhBBCSK+X1ODG6XRizJgxePHFF+O67ujRo6iqquI/0tPTu2mEhBBCCOlrxMl88rlz52Lu3LlxX5eeng6dTpf4ARFCCCGkz+uTOTdjx45FVlYWZs2ahR07drR5rtfrhc1mi/kghBBCSOrqU8FNVlYW1q5diw8++AAffPABcnNzMWPGDOzdu7fVa0pKSqDVavmP3NzcHhwxIYQQQnqagOM4LtmDAACBQICPPvoIl19+eVzXTZ8+HQMHDsQ//vGPFh/3er3wer381zabDbm5ubBardBoNF0ZMiGEEEJ6iM1mg1ar7dD7d1JzbhJh4sSJ+Prrr1t9nGEYMAzTgyMihBBCSDL1qWWpluzfvx9ZWVnJHgYhhBBCeomkztw4HA6cOHGC/7qsrAz79++HwWDAwIEDsWzZMpw5cwZ///vfAQDPP/888vPzMXLkSHg8Hrz22mv46quv8PnnnyfrWyCEEEJIL5PU4Ob7779HcXEx//Xdd98NAFi4cCHWrVuHqqoqVFZW8o/7fD788Y9/xJkzZ6BQKFBUVIQvvvgi5h6EEEII6d96TUJxT4knIYkQQgghvUM87999PueGEEIIISQaBTeEEEIISSkU3BBCCCEkpVBwQwghhJCUQsENIYQQQlIKBTeEEEIISSkU3BBCCCEkpVBwQwghhJCUQsENIYQQQlIKBTeEEEIISSkU3BBCCCEkpVBwQwghhJCUQsENIYQQQlIKBTeEEEIISSniZA+AEEIIIX1DkOVgcfnQ6PShwelDg8OHRqcX9Y7QsUanD/UOLwrSVXjiitFJGycFN4QQQkg/xbIcLG4/GhxeNISDk9jPfWhwesNBjA9mlw8s1/593f5g9w++DRTcEEIIISmio8FK5POOBitNaeUSpCmlSFNJYVBKkaZiQl8rpTCoGOToZIn/5uJAwQ0hhBDSS0WCleiln4ZwwNL088hHZ4IVjUzMBygxwUokeFEySFOFghe9UgqJqHen7FJwQwghhPQQluVgdftjlnrqnT40RnJXwp9HZle6EqwYVQwMfLASClCafm5U9Y1gJV4U3BBCCCGddC5YaTqbEgpWzn0eCljMLj+CnYhWms+s9K9gJV4U3BBCCCFhLMvB5vFHVf80rwQ6l7sSylnpTLCibjKzYmyy/HPuOAO9QgqpuH8HK/Gi4IYQQkjKigQrrc+mnJtxqXckLlg5l68Sm7tCwUrPoOCGEEJIn8FxHGzuAOr5ip/mwUpjOJ+lwemD2elDoJPBSlo4sbatYCVNyUCvlIARi7rhuyWdRcENIYSQpOE4DjZPIGb2pOlyUHTybWNngxVGHLXcw5xbBmohj4WClb6PghtCCCEJ0zRYaW85yOzywR/sXLBiiM5TUUphUJ2bYTmXYEvBSn9EwQ0hhJBWRYKV5pvBNd0YrmvBiipqZqWlYKXpchAFK6QtFNwQQkg/Eh2sxCz9tBCsNIbzWjobrDSrBIpaAorMqqSppNArpJBJKFghiUPBDSGE9GEcx8HuDZxb+uEDk+Zb7Td0IVhRSkXhmZS2d7GNzLZQsEKSiYIbQgjpRSLBSqPD16GN4RqdPviCbNzPEwlWDEoGRmXrybUUrJC+iIIbQgjpRhzHweENtBmgRG8M19lgRSEVtdDAMPzf8MZwRgpWSD9BwQ0hhMSB4zg4fcFzOSrh5Z5znzffa8UXiD9YkUtEfKNCfq+VJsEKP8uiZCCXUrBCSAQFN4SQfi0SrPBBSks710aSax2hzzsTrPAzKy3mqzTPXaFghZDOo+CGEJJSOI6DyxdsIZm2pYTb0IyLt4szK83zVZrnrlCwQkjPoeCGENKrRYKVmNyUNpJr6x3eTgUrMomQb1oYWfoxqmLLlg1RAYtCSr8+Cemt6F8nIaRHRQcrDVGzJy3ttZKoYKVpt+Wmu9hSsEJIaqF/zYSQLok3WGlweuHxdz1YORekMFEbxZ3LXaFghZD+i/71E0KacfkCreapRJaGojeG60ywwoiF/FJP7IwKE/M5BSuEkHjRbwtC+gG3LxhTCRQdoNQ7zm2zH9mK3+0Pxv0cUrEwtBlc1C62fCVQVOJtJKBRSEUQCATd8N0SQvq7pAY327Ztw1NPPYU9e/agqqoKH330ES6//PI2ryktLcXdd9+NH3/8Ebm5uXjwwQexaNGiHhkvIb1FW8FKS7vadjVY4Zd7wl8blUyzXWyVFKwQQnqJpAY3TqcTY8aMwU033YQrr7yy3fPLysrwi1/8ArfeeiveeustfPnll7jllluQlZWF2bNn98CICekeHn8wJl+lpUqg6Mc6G6xEz6ZEttznOy+HgxQjBSuEkD4uqcHN3LlzMXfu3A6fv3btWuTn5+OZZ54BAAwfPhxff/01nnvuOQpuSK/SNFhpbxdbl68TwYpIGJOvYmxjF9s0FUPBCiGk3+hTOTe7du3CzJkzY47Nnj0bS5cubfUar9cLr9fLf22z2bpreCSFRYKVprvY1od3re2OYKWlSqBzn0uhYsQUrBBCSAv6VHBTXV2NjIyMmGMZGRmw2Wxwu92Qy+XNrikpKcGjjz7aU0MkfYTHH4yp9mltF9uGcPDi7ESwIhEJYvZRaWk2JXrPFQpWCCEkMfpUcNMZy5Ytw913381/bbPZkJubm8QRke7QXrASXQnU6PTB4Q3E/RwSkSBmM7joYKXZxnAqKdQUrBBCSFL0qeAmMzMTNTU1Mcdqamqg0WhanLUBAIZhwDBMTwyPJFAkWOlI2XKigpXWdrGNBDAUrBBCSN/Qp4KbyZMn49NPP405tnnzZkyePDlJIyId5Q0EW136iU62jeSydCZYEQsF5yqBWtjFloIVQgjpH5Ia3DgcDpw4cYL/uqysDPv374fBYMDAgQOxbNkynDlzBn//+98BALfeeivWrFmD++67DzfddBO++uorvPvuu/jkk0+S9S30W9HBSmRL/dhg5dzyUKPDB3sng5WWmha2lFybpmKgkVGwQgghJMnBzffff4/i4mL+60huzMKFC7Fu3TpUVVWhsrKSfzw/Px+ffPIJ7rrrLqxevRoDBgzAa6+9RmXgCeALsC0sAbWQaNtTwYqSgUZOwQohhJD4CTiO45I9iJ5ks9mg1WphtVqh0WiSPZxuEwlWOrLlfkMngxWRUNA8NyVqF9voxFsKVgghhHRFPO/ffSrnpj/zBViYXbHVQK0l19Y7vLB7ui9YMShDu9hSsEIIIaQ3ouAmSfxBtsWclWb5K+Fdbm2dDFb0iti8lEjzwpb2WtHIJBAKKVghhBDSt1FwkyCBIBszg9JSsBLd1LArwcq5mZW2c1coWCGEENIfUXDTCZE0JYFAwH/+bVkjFrz2bVz3EQrAJ89GGhga29jFViunYIUQQghpDwU3nRCdZxL5PE0l5YMVPmBp1m05ekmIoWCFEEII6QYU3CTIsHQ1TjxxCQUrhBBCSJJRcJMgFNQQQgghvYMw2QMghBBCCEkkCm4IIYQQklIouCGEEEJISqHghhBCCCEphYIbQgghhKQUCm4IIYQQklIouCGEEEJISqHghhBCCCEphYIbQgghhKQUCm4IIYQQklIouCGEEEJISqHghhBCCCEphYIbQgghhKSUftcVnOM4AIDNZkvySAghhBDSUZH37cj7eFv6XXBjt9sBALm5uUkeCSGEEELiZbfbodVq2zxHwHUkBEohLMvi7NmzUKvVEAgEyR5O0tlsNuTm5uLUqVPQaDTJHk7Kote5Z9Dr3DPode459Fqfw3Ec7HY7srOzIRS2nVXT72ZuhEIhBgwYkOxh9Doajabf/8PpCfQ69wx6nXsGvc49h17rkPZmbCIooZgQQgghKYWCG0IIIYSkFApu+jmGYfDwww+DYZhkDyWl0evcM+h17hn0Ovcceq07p98lFBNCCCEktdHMDSGEEEJSCgU3hBBCCEkpFNwQQgghJKVQcEMIIYSQlELBTQratm0bLrvsMmRnZ0MgEGDDhg38Y36/H/fffz9Gjx4NpVKJ7Oxs3HDDDTh79mzMPRobG7FgwQJoNBrodDrcfPPNcDgcPfyd9G5tvc4A8Mgjj+C8886DUqmEXq/HzJkz8e2338acQ69zx7T3Wke79dZbIRAI8Pzzz8ccp9e6fe29zosWLYJAIIj5mDNnTsw59Dq3ryM/z4cPH8a8efOg1WqhVCrxs5/9DJWVlfzjHo8HS5YsQVpaGlQqFa666irU1NT04HfRu1Fwk4KcTifGjBmDF198sdljLpcLe/fuxfLly7F37158+OGHOHr0KObNmxdz3oIFC/Djjz9i8+bN2LhxI7Zt24bf/va3PfUt9Altvc4AMGzYMKxZswYHDx7E119/jUGDBuHiiy9GXV0dfw69zh3T3msd8dFHH+Gbb75BdnZ2s8fotW5fR17nOXPmoKqqiv/417/+FfM4vc7ta+91PnnyJH7+85/jvPPOQ2lpKQ4cOIDly5dDJpPx59x11134z3/+g/feew9bt27F2bNnceWVV/bUt9D7cSSlAeA++uijNs/57rvvOABcRUUFx3Ecd+jQIQ4At3v3bv6cTZs2cQKBgDtz5kx3DrfP6sjrbLVaOQDcF198wXEcvc6d1dprffr0aS4nJ4f74YcfuLy8PO65557jH6PXOn4tvc4LFy7k5s+f3+o19DrHr6XX+dprr+V+/etft3qNxWLhJBIJ99577/HHDh8+zAHgdu3a1V1D7VNo5obAarVCIBBAp9MBAHbt2gWdTofzzz+fP2fmzJkQCoXNllVIx/h8PrzyyivQarUYM2YMAHqdE4llWfzmN7/Bvffei5EjRzZ7nF7rxCktLUV6ejoKCwtx2223oaGhgX+MXueuY1kWn3zyCYYNG4bZs2cjPT0d//d//xezdLVnzx74/X7MnDmTP3beeedh4MCB2LVrVxJG3ftQcNPPeTwe3H///bj++uv5pmzV1dVIT0+POU8sFsNgMKC6ujoZw+yzNm7cCJVKBZlMhueeew6bN2+G0WgEQK9zIq1atQpisRh/+MMfWnycXuvEmDNnDv7+97/jyy+/xKpVq7B161bMnTsXwWAQAL3OiVBbWwuHw4GVK1dizpw5+Pzzz3HFFVfgyiuvxNatWwGEXmepVMr/QRqRkZFBr3NYv+sKTs7x+/245pprwHEcXnrppWQPJyUVFxdj//79qK+vx6uvvoprrrkG3377bbM3ANJ5e/bswerVq7F3714IBIJkDyelXXfddfzno0ePRlFREYYMGYLS0lJcdNFFSRxZ6mBZFgAwf/583HXXXQCAsWPHYufOnVi7di2mT5+ezOH1GTRz009FApuKigps3ryZn7UBgMzMTNTW1sacHwgE0NjYiMzMzJ4eap+mVCpRUFCASZMm4fXXX4dYLMbrr78OgF7nRNm+fTtqa2sxcOBAiMViiMViVFRU4I9//CMGDRoEgF7r7jJ48GAYjUacOHECAL3OiWA0GiEWizFixIiY48OHD+erpTIzM+Hz+WCxWGLOqampodc5jIKbfigS2Bw/fhxffPEF0tLSYh6fPHkyLBYL9uzZwx/76quvwLIs/u///q+nh5tSWJaF1+sFQK9zovzmN7/BgQMHsH//fv4jOzsb9957L/773/8CoNe6u5w+fRoNDQ3IysoCQK9zIkilUvzsZz/D0aNHY44fO3YMeXl5AIAJEyZAIpHgyy+/5B8/evQoKisrMXny5B4db29Fy1IpyOFw8H9JAUBZWRn2798Pg8GArKwsXH311di7dy82btyIYDDIr9EaDAZIpVIMHz4cc+bMweLFi7F27Vr4/X78/ve/x3XXXddiiW1/1dbrnJaWhieeeALz5s1DVlYW6uvr8eKLL+LMmTP45S9/CQD0Osehrdd64MCBzQJ0iUSCzMxMFBYWAqDXuqPaep0NBgMeffRRXHXVVcjMzMTJkydx3333oaCgALNnzwZAr3NHtffzfO+99+Laa6/FtGnTUFxcjM8++wz/+c9/UFpaCgDQarW4+eabcffdd8NgMECj0eCOO+7A5MmTMWnSpCR9V71Mssu1SOJt2bKFA9DsY+HChVxZWVmLjwHgtmzZwt+joaGBu/766zmVSsVpNBruxhtv5Ox2e/K+qV6ordfZ7XZzV1xxBZednc1JpVIuKyuLmzdvHvfdd9/F3INe545p67VuSdNScI6j17oj2nqdXS4Xd/HFF3Mmk4mTSCRcXl4et3jxYq66ujrmHvQ6t68jP8+vv/46V1BQwMlkMm7MmDHchg0bYu7hdru522+/ndPr9ZxCoeCuuOIKrqqqqoe/k95LwHEc1wMxFCGEEEJIj6CcG0IIIYSkFApuCCGEEJJSKLghhBBCSEqh4IYQQgghKYWCG0IIIYSkFApuCCGEEJJSKLghhBBCSEqh4IYQ0mGLFi3C5Zdfzn89Y8YMLF26tNP3GzRoEJ5//vk2zxEIBNiwYUOnn6MjysvLIRAIsH///m59HkJIz6DghpBeqK6uDrfddhsGDhwIhmGQmZmJ2bNnY8eOHQm5f9MgpbM+/PBDPPbYY10fECGEJBD1liKkF7rqqqvg8/nw5ptvYvDgwaipqcGXX36JhoaGZA8thsFgSPYQejW/3w+JRJLsYcTojWMiJNFo5oaQXsZisWD79u1YtWoViouLkZeXh4kTJ2LZsmWYN28eAOCmm27CpZdeGnOd3+9Heno6Xn/9dQDA+++/j9GjR0MulyMtLQ0zZ86E0+nEI488gjfffBMff/wxBAIBBAIB35Dv4MGDuPDCC/lrfvvb38LhcLQ61qbLUl6vF/fffz9yc3PBMAwKCgr48bTGbrfj+uuvh1KpRE5ODl588cU2z29vjCzLYsWKFRgwYAAYhsHYsWPx2Wefxdzju+++w7hx4yCTyXD++edj3759bT4nEFpCe+yxx9ocq0AgwEsvvYR58+ZBqVTiiSeeAAB8/PHHGD9+PGQyGQYPHoxHH30UgUAAAMBxHB555BF+li47Oxt/+MMf+Hv+9a9/xdChQyGTyZCRkYGrr746ZkxNl/XGjh2LRx55pEtjIqTPS3JvK0JIE36/n1OpVNzSpUs5j8fT4jk7duzgRCIRd/bsWf7Yhx9+yCmVSs5ut3Nnz57lxGIx9+yzz3JlZWXcgQMHuBdffJGz2+2c3W7nrrnmGm7OnDlcVVUVV1VVxXm9Xs7hcHBZWVnclVdeyR08eJD78ssvufz8/JhmfgsXLuTmz5/Pfz19+nTuzjvv5L++5ppruNzcXO7DDz/kTp48yX3xxRfc22+/3er3mpeXx6nVaq6kpIQ7evQo98ILL3AikYj7/PPP+XMAcB999BHHcVyHxvjss89yGo2G+9e//sUdOXKEu++++ziJRMIdO3aM4ziOs9vtnMlk4n71q19xP/zwA/ef//yHGzx4MAeA27dvX5fHmp6ezv3tb3/jTp48yVVUVHDbtm3jNBoNt27dOu7kyZPc559/zg0aNIh75JFHOI7juPfee4/TaDTcp59+ylVUVHDffvst98orr3Acx3G7d+/mRCIRt379eq68vJzbu3cvt3r16pgxNW0QOmbMGO7hhx/u0pgI6esouCGkF3r//fc5vV7PyWQybsqUKdyyZcu4//3vfzHnjBgxglu1ahX/9WWXXcYtWrSI4ziO27NnDweAKy8vb/H+TYMUjuO4V155hdPr9ZzD4eCPffLJJ5xQKOQ7P7cV3Bw9epQDwG3evLnD32deXh43Z86cmGPXXnstN3fuXP7r6OCmI2PMzs7mnnjiiZh7/uxnP+Nuv/12juM47uWXX+bS0tI4t9vNP/7SSy91KLjpyFiXLl0ac85FF13EPfnkkzHH/vGPf3BZWVkcx3HcM888ww0bNozz+XzNnvODDz7gNBoNZ7PZWh1TR4KbeMdESF9Hy1KE9EJXXXUVzp49i3//+9+YM2cOSktLMX78eKxbt44/55ZbbsEbb7wBAKipqcGmTZtw0003AQDGjBmDiy66CKNHj8Yvf/lLvPrqqzCbzW0+5+HDhzFmzBgolUr+2NSpU8GyLI4ePdrumPfv3w+RSITp06fH9b1Onjy52deHDx/u1BhtNhvOnj2LqVOnxlw3depU/p6HDx9GUVERZDJZq2PoyljPP//8mK//97//YcWKFVCpVPzH4sWLUVVVBZfLhV/+8pdwu90YPHgwFi9ejI8++ohfHpo1axby8vIwePBg/OY3v8Fbb70Fl8vVobF2ZUyE9HUU3BDSS8lkMsyaNQvLly/Hzp07sWjRIjz88MP84zfccAN++ukn7Nq1C//85z+Rn5+PCy64AAAgEomwefNmbNq0CSNGjMBf/vIXFBYWoqysrNvGK5fLu+3efUl04AUADocDjz76KPbv389/HDx4EMePH4dMJkNubi6OHj2Kv/71r5DL5bj99tsxbdo0+P1+qNVq7N27F//617+QlZWFhx56CGPGjIHFYgEACIVCcBwX83x+v7/LYyKkr6PghpA+YsSIEXA6nfzXaWlpuPzyy/HGG29g3bp1uPHGG2POFwgEmDp1Kh599FHs27cPUqkUH330EQBAKpUiGAzGnD98+HD873//i3mOHTt2QCgUorCwsN3xjR49GizLYuvWrXF9X998802zr4cPH97iue2NUaPRIDs7u1nJ/I4dOzBixAj+HgcOHIDH42l1DIkYa8T48eNx9OhRFBQUNPsQCkO/guVyOS677DK88MILKC0txa5du3Dw4EEAgFgsxsyZM/HnP/8ZBw4cQHl5Ob766isAgMlkQlVVFf9cNputQwFsR8ZESF9GpeCE9DINDQ345S9/iZtuuglFRUVQq/9/e/cTClsUxwH8+3iGGd78If8azIjYmcnCyIIolEYkf4o0JbNULBTRzCiriQ17cxcsJNn4G7KQNGWBpGyE5YyNCJuZ31uZmtd7Db16r7l9P8tzz7n9zuJ2v3XOufcHzs7O4Pf70dHREdd3eHgYTqcTkUgELpcr1h4MBnF4eIiWlhbk5eUhGAwiHA7HXsRWqxV7e3u4ublBTk4ODAYDBgYG4PV64XK54PP5EA6HMTIygsHBQeTn5yes22q1wuVyYWhoCAsLC7DZbLi/v0coFEJvb+8fx52cnMDv96OzsxP7+/tYW1vD1tbWb/t+psbx8XF4vV6UlZXBbrcjEAjg/PwcKysrAID+/n5MTU3B7XZjcnISd3d3mJubSzi/r9b6wePxwOl0oqSkBN3d3UhJScHFxQWurq4wOzsLRVEQiUTgcDig0+mwvLwMrVYLi8WCzc1N3N7eor6+HiaTCdvb24hGo7Gw2dTUBEVR0N7eDqPRCI/Hg9TU1ITzSFQTUdL735t+iCje+/u7TExMSHV1tRgMBtHpdFJZWSnT09Py+voa1zcajYrFYpG2tra49uvra2ltbZXc3FxJT0+XiooKWVxcjF0PhULS3NwsWVlZAkCOjo5EROTy8lIaGxslIyNDsrOzxe12y/Pzc2xcotNSb29vMjY2JoWFhaLRaKS8vFyWlpb+OFeLxSIzMzPS09MjOp1OCgoK4k4DicRvKP5MjZFIRHw+n5jNZklLSxObzSY7Oztx9zw9PRWbzSYajUbsdrusr69/akPxV2v9sLu7K3V1daLVakWv10tNTU3sRNTGxoY4HA7R6/WSmZkptbW1cnBwICIix8fH0tDQICaTSbRarVRVVcnq6mrsvk9PT9LX1yd6vV6Ki4tFUZTfbij+ak1Eye6byC8LtkSUNF5eXmA2mxEIBNDV1fW/y1E1q9WK0dHRv/rdBBH9G1yWIkpC0WgUj4+PmJ+fh9FojH3cj4iIGG6IktLDwwNKS0tRVFQERVHw/TsfZSKiD1yWIiIiIlXhmT8iIiJSFYYbIiIiUhWGGyIiIlIVhhsiIiJSFYYbIiIiUhWGGyIiIlIVhhsiIiJSFYYbIiIiUhWGGyIiIlKVn3brndE98jN8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,slope_fn(beta_s,x),label=\"No weights--S\",c=\"olive\")\n",
    "plt.fill_between(x,slope_fn(low_beta_s,x),slope_fn(high_beta_s,x),alpha=0.3,color=\"olive\")\n",
    "plt.plot(x,slope_fn(beta_s_GL,x),label=\"Weights--S\",c=\"cornflowerblue\")\n",
    "plt.fill_between(x,slope_fn(low_beta_s_GL,x),slope_fn(high_beta_s_GL,x),alpha=0.3,color=\"cornflowerblue\")\n",
    "plt.plot(x,slope_fn(beta_i,x),label=\"No weights--I\")\n",
    "plt.plot(x,slope_fn(beta_i_GL,x),label=\"Weights--I\")\n",
    "plt.ylabel(\"log relative risk\")\n",
    "plt.xlabel(\"Systolic blood pressure\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Itoga"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "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>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>129.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>134.6</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>139.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>143.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>149.5</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Author   dose       cases     controls  subjects  logOR     low   high  \\\n",
       "5  Itoga  129.8  249.666667  2738.333333    2988.0  0.000  0.0000  0.000   \n",
       "6  Itoga  134.6  249.666667  6048.333333    6298.0 -0.062 -0.1500 -0.023   \n",
       "7  Itoga  139.0  249.666667  9341.333333    9591.0  0.077  0.0004  0.150   \n",
       "8  Itoga  143.7  249.666667  8543.333333    8793.0  0.110  0.0220  0.210   \n",
       "9  Itoga  149.5  249.666667  4401.333333    4651.0  0.190  0.1100  0.280   \n",
       "\n",
       "   std_error  \n",
       "5   0.000000  \n",
       "6   0.032071  \n",
       "7   0.037778  \n",
       "8   0.047475  \n",
       "9   0.042929  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/reductions/solvers/solving_chain.py:354: FutureWarning: \n",
      "    You specified your problem should be solved by ECOS. Starting in\n",
      "    CXVPY 1.6.0, ECOS will no longer be installed by default with CVXPY.\n",
      "    Please either add an explicit dependency on ECOS or switch to our new\n",
      "    default solver, Clarabel, by either not specifying a solver argument\n",
      "    or specifying ``solver=cp.CLARABEL``.\n",
      "    \n",
      "  warnings.warn(ECOS_DEP_DEPRECATION_MSG, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "Nx_i = study_df_i[\"subjects\"].to_numpy().astype(np.float64)\n",
    "M1x_i = np.sum(study_df_i[\"cases\"].to_numpy().astype(np.float64))\n",
    "A0_i = M1x_i*Nx_i[1:]/(Nx_i.sum())\n",
    "L_i = study_df_i[\"logOR\"][1:].to_numpy()\n",
    "L_i_p = L_i.copy()\n",
    "L_i_p[0] = -0.062\n",
    "x_i = study_df_i[\"dose\"][1:].to_numpy()\n",
    "x_i = x_i - 129.8\n",
    "v_i = study_df_i[\"std_error\"][1:].to_numpy()\n",
    "v_ip = v_i.copy()\n",
    "# v_ip[0] = \n",
    "\n",
    "it_fit_ex = IF(L_i_p,A0_i,Nx_i,M1x_i)\n",
    "A_i, B_i, a0_i, b0_i = it_fit_ex.convexProgram()\n",
    "\n",
    "C = covariance_matrix(A_i,B_i,a0_i,b0_i,v_ip**2)\n",
    "inv_C_i = np.linalg.inv(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([213.99462598, 372.59450971, 352.59471619, 201.36317414])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "107.7863073121539"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a0_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00102853, 0.00087018, 0.00108647, 0.00090527],\n",
       "       [0.00087018, 0.00142716, 0.00138118, 0.00115083],\n",
       "       [0.00108647, 0.00090527, 0.00225385, 0.00143687],\n",
       "       [0.00138118, 0.00115083, 0.00143687, 0.00184292]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Itoga corrected estimates\n",
    "vb_star = 1/(np.dot(x_i,np.dot(inv_C_i,x_i)))\n",
    "b_star = vb_star*(np.dot(x_i,np.dot(inv_C_i,L_i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.018745169004179014"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Itoga unadjusted estimates\n",
    "vb = 1/(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_ip)**2),x_i)))\n",
    "b = vb*(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_ip)**2),L_i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0077151064467757356"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a jump of size of order of magnitude. Why is the jump so big maybe? How do we know the difference induced by our method returns the results that are \"most correct?\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat Itoga w/out negative term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx_i = study_df_i[\"subjects\"].to_numpy().astype(np.float64)[1:]\n",
    "M1x_i = np.sum(study_df_i[\"cases\"].to_numpy().astype(np.float64)[1:])\n",
    "A0_i = M1x_i*Nx_i[1:]/(Nx_i.sum())\n",
    "L_i = study_df_i[\"logOR\"][2:].to_numpy()\n",
    "x_i = study_df_i[\"dose\"][2:].to_numpy()\n",
    "v_i = study_df_i[\"std_error\"][2:].to_numpy()\n",
    "\n",
    "it_fit_ex = IF(L_i,A0_i,Nx_i,M1x_i)\n",
    "A_i, B_i, a0_i, b0_i = it_fit_ex.convexProgram()\n",
    "\n",
    "C = covariance_matrix(A_i,B_i,a0_i,b0_i,v_i**2)\n",
    "inv_C_i = np.linalg.inv(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([322.48570552, 305.22983423, 174.39182048])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00142716, 0.00110126, 0.00087457],\n",
       "       [0.00110126, 0.00225385, 0.00108721],\n",
       "       [0.00087457, 0.00108721, 0.00184292]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_nn = C\n",
    "inv_C_nn = np.linalg.inv(C_nn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0008633607306968095"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vb_star = 1/(np.dot(x_i,np.dot(inv_C_nn,x_i)))\n",
    "b_star = vb_star*(np.dot(x_i,np.dot(inv_C_nn,L_i)))\n",
    "b_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.000859255738818118"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vb = 1/(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_i**2)),x_i)))\n",
    "b = vb*(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_i**2)),L_i)))\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so now the estimates are almost exactly the same. I really think the reason behind the jump here is the fact that we have that negative in the first non-reference exposure level. But why does that make such a big difference? I'm not entirely sure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Itoga simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.146112\n",
      "         Iterations 8\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:               outcomes   No. Observations:                40000\n",
      "Model:                          Logit   Df Residuals:                    39995\n",
      "Method:                           MLE   Df Model:                            4\n",
      "Date:                Thu, 23 May 2024   Pseudo R-squ.:                0.002329\n",
      "Time:                        13:39:07   Log-Likelihood:                -5844.5\n",
      "converged:                       True   LL-Null:                       -5858.1\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.740e-05\n",
      "==================================================================================\n",
      "                     coef    std err          z      P>|z|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------\n",
      "categories_3.0     0.0853      0.088      0.970      0.332      -0.087       0.258\n",
      "categories_4.0     0.1749      0.084      2.070      0.038       0.009       0.341\n",
      "categories_5.0     0.3418      0.079      4.343      0.000       0.188       0.496\n",
      "categories_6.0     0.6223      0.192      3.243      0.001       0.246       0.998\n",
      "ones              -3.5435      0.062    -56.760      0.000      -3.666      -3.421\n",
      "==================================================================================\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Performing the Logit regression to recover the log odds ratios\n",
    "cats_df = pd.DataFrame(cats_out_df,columns=[\"outcomes\",\"categories\"])\n",
    "cats_df = pd.get_dummies(cats_df, columns=[\"categories\"], dtype=np.float32)\n",
    "cats_df[\"ones\"] = np.ones(len(cats_df))\n",
    "to_reg = cats_df.loc[:, [x for x in cats_df.columns if x not in (\"outcomes\", \"categories_2.0\")]]\n",
    "\n",
    "results = sm.Logit(endog=cats_df[\"outcomes\"], exog=to_reg).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sanity check on logit, if our estimate is actually closer than unadjusted estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.146110\n",
      "         Iterations 8\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                40000\n",
      "Model:                          Logit   Df Residuals:                    39998\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Thu, 23 May 2024   Pseudo R-squ.:                0.002343\n",
      "Time:                        15:48:56   Log-Likelihood:                -5844.4\n",
      "converged:                       True   LL-Null:                       -5858.1\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.611e-07\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -3.6303      0.060    -60.823      0.000      -3.747      -3.513\n",
      "x1             0.0255      0.005      5.222      0.000       0.016       0.035\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "\n",
    "x = sm.add_constant(df[:,1])\n",
    "results = sm.Logit(endog=df[:,0], exog=x).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "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>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>129.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>134.6</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>139.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>143.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>149.5</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Author   dose       cases     controls  subjects  logOR     low   high  \\\n",
       "5  Itoga  129.8  249.666667  2738.333333    2988.0  0.000  0.0000  0.000   \n",
       "6  Itoga  134.6  249.666667  6048.333333    6298.0 -0.062 -0.1500 -0.023   \n",
       "7  Itoga  139.0  249.666667  9341.333333    9591.0  0.077  0.0004  0.150   \n",
       "8  Itoga  143.7  249.666667  8543.333333    8793.0  0.110  0.0220  0.210   \n",
       "9  Itoga  149.5  249.666667  4401.333333    4651.0  0.190  0.1100  0.280   \n",
       "\n",
       "   std_error  \n",
       "5   0.000000  \n",
       "6   0.032071  \n",
       "7   0.037778  \n",
       "8   0.047475  \n",
       "9   0.042929  "
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "study_df_i[\"rounded_cases\"] = np.round(study_df_i[\"cases\"].to_numpy())\n",
    "study_df_i[\"rounded_controls\"] = np.round(study_df_i[\"controls\"].to_numpy())"
   ]
  },
  {
   "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></th>\n",
       "      <th>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "      <th>rounded_cases</th>\n",
       "      <th>rounded_controls</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>0.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>250.0</td>\n",
       "      <td>2738.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>4.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "      <td>250.0</td>\n",
       "      <td>6048.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>9.2</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "      <td>250.0</td>\n",
       "      <td>9341.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>13.9</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "      <td>250.0</td>\n",
       "      <td>8543.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>19.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "      <td>250.0</td>\n",
       "      <td>4401.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Author  dose       cases     controls  subjects  logOR     low   high  \\\n",
       "5  Itoga   0.0  249.666667  2738.333333    2988.0  0.000  0.0000  0.000   \n",
       "6  Itoga   4.8  249.666667  6048.333333    6298.0 -0.062 -0.1500 -0.023   \n",
       "7  Itoga   9.2  249.666667  9341.333333    9591.0  0.077  0.0004  0.150   \n",
       "8  Itoga  13.9  249.666667  8543.333333    8793.0  0.110  0.0220  0.210   \n",
       "9  Itoga  19.7  249.666667  4401.333333    4651.0  0.190  0.1100  0.280   \n",
       "\n",
       "   std_error  rounded_cases  rounded_controls  \n",
       "5   0.000000          250.0            2738.0  \n",
       "6   0.032071          250.0            6048.0  \n",
       "7   0.037778          250.0            9341.0  \n",
       "8   0.047475          250.0            8543.0  \n",
       "9   0.042929          250.0            4401.0  "
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df_i[\"dose\"] -= study_df_i[\"dose\"].to_numpy()[0]\n",
    "study_df_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ref_list = []\n",
    "for i in range(1,len(study_df_i)):\n",
    "    df_ref_cases = np.array([[study_df_i[\"dose\"].to_numpy()[i] for _ in range(int(study_df_i[\"rounded_cases\"].to_numpy()[i]))],\n",
    "                    [1 for _ in range(int(study_df_i[\"rounded_cases\"].to_numpy()[i]))]]).T\n",
    "    df_ref_controls = np.array([[study_df_i[\"dose\"].to_numpy()[i] for _ in range(int(study_df_i[\"rounded_controls\"].to_numpy()[i]))],\n",
    "                    [0 for _ in range(int(study_df_i[\"rounded_controls\"].to_numpy()[i]))]]).T\n",
    "    df_ref_list.append(np.vstack((df_ref_cases,df_ref_controls)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ref = np.vstack(df_ref_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "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>dose</th>\n",
       "      <th>outcome</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4.8</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.8</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.8</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.8</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.8</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   dose  outcome\n",
       "0   4.8      1.0\n",
       "1   4.8      1.0\n",
       "2   4.8      1.0\n",
       "3   4.8      1.0\n",
       "4   4.8      1.0"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_ref = pd.DataFrame(df_ref, columns=[\"dose\",\"outcome\"])\n",
    "df_ref.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.148448\n",
      "         Iterations 7\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                outcome   No. Observations:                29333\n",
      "Model:                          Logit   Df Residuals:                    29331\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Thu, 23 May 2024   Pseudo R-squ.:                0.001614\n",
      "Time:                        13:33:25   Log-Likelihood:                -4354.4\n",
      "converged:                       True   LL-Null:                       -4361.5\n",
      "Covariance Type:            nonrobust   LLR p-value:                 0.0001752\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -3.6289      0.084    -43.278      0.000      -3.793      -3.465\n",
      "dose           0.0246      0.007      3.767      0.000       0.012       0.037\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "# Step 4: Build the logistic regression model without an intercept\n",
    "X = df_ref['dose']\n",
    "y = df_ref['outcome']\n",
    "# Add a constant column to X if needed (here we will not add it to exclude the intercept)\n",
    "X = sm.add_constant(X)\n",
    "logit_model = sm.Logit(y, X).fit()\n",
    "# Step 5: Evaluate the model\n",
    "print(logit_model.summary())\n",
    "# # Optional: Make predictions\n",
    "# predictions = logit_model.predict(X)\n",
    "# predicted_classes = (predictions > 0.5).astype(int)\n",
    "# # Optional: Print model coefficients\n",
    "# print(f'Coefficient (slope) for dose: {logit_model.params[\"dose\"]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/p0/hl6rzzgd2v5802nvrkllqms40000gp/T/ipykernel_8866/3760040021.py:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
      "  np.exp(logit_model.params[1])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0248559494152192"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(logit_model.params[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/p0/hl6rzzgd2v5802nvrkllqms40000gp/T/ipykernel_8866/2289857680.py:1: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n",
      "  logit_model.params[1]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.024552065558216256"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logit_model.params[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.2 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}