GL_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as np\n",
    "from jax import random, jit\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from src.IterativeFitting import IterativeFitting as IF\n",
    "from src.CorrFuncs import covariance_matrix, trend_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Platform 'METAL' is experimental and not all JAX functionality may be correctly supported!\n",
      "2024-06-13 16:15:45.720583: W pjrt_plugin/src/mps_client.cc:563] WARNING: JAX Apple GPU support is experimental and not all JAX functionality is correctly supported!\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Metal device set to: Apple M1 Pro\n",
      "\n",
      "systemMemory: 16.00 GB\n",
      "maxCacheSize: 5.33 GB\n",
      "\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "uniform() got an unexpected keyword argument 'low'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[6], line 17\u001b[0m\n\u001b[1;32m     14\u001b[0m n \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m40000\u001b[39m\n\u001b[1;32m     16\u001b[0m \u001b[38;5;66;03m# Constructing covariate vectors for additional covariates to adjust\u001b[39;00m\n\u001b[0;32m---> 17\u001b[0m age_vector \u001b[38;5;241m=\u001b[39m \u001b[43mrandom\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43muniform\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlow\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mhigh\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43msize\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mn\u001b[49m\u001b[43m,\u001b[49m\u001b[43mkey\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     18\u001b[0m income_vector \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39muniform(low\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m,high\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m,size\u001b[38;5;241m=\u001b[39mn)\n\u001b[1;32m     19\u001b[0m sex_vector \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mrandom\u001b[38;5;241m.\u001b[39mbinomial(n\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m,p\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.6\u001b[39m,size\u001b[38;5;241m=\u001b[39mn)\n",
      "\u001b[0;31mTypeError\u001b[0m: uniform() got an unexpected keyword argument 'low'"
     ]
    }
   ],
   "source": [
    "# Setting seed\n",
    "golden_seed = 6\n",
    "# np.random.seed(golden_seed)\n",
    "key = random.key(golden_seed)\n",
    "\n",
    "# Initializing simulation slope and intercept parameters\n",
    "beta0 = -0.1\n",
    "beta1 = 0.02\n",
    "beta_age = 0.05\n",
    "beta_income = -0.05\n",
    "beta_sex = 1.0\n",
    "\n",
    "# Number of samples\n",
    "n = 40000\n",
    "\n",
    "# Constructing covariate vectors for additional covariates to adjust\n",
    "age_vector = random.uniform(minval=0,maxval=1,size=n,key=key)\n",
    "income_vector = np.random.uniform(low=0,high=1,size=n)\n",
    "sex_vector = np.random.binomial(n=1,p=0.6,size=n)\n",
    "\n",
    "# Generating plausible x values to throw into probability generator\n",
    "xs = np.sort(np.random.uniform(low=0,high=20,size=n))\n",
    "\n",
    "# Getting quantiles of sorted xs\n",
    "quantile_indices = np.array([n//4, n//2, 3*n//4])\n",
    "\n",
    "# Getting exposure levels (randomly chosen from quantile ranges)\n",
    "exposure_levels = [np.random.choice(xs[:quantile_indices[0]])]\n",
    "for i in range(1,len(quantile_indices)):\n",
    "    exposure_levels.append(np.random.choice(xs[quantile_indices[i-1]:quantile_indices[i]]))\n",
    "exposure_levels.append(np.random.choice(xs[quantile_indices[-1]:]))\n",
    "exposure_indices = np.array([np.where(xs == elem)[0][0] for elem in exposure_levels])\n",
    "\n",
    "# Function to generate probabilities of being a case v non-case\n",
    "other_covs = beta0 + (beta_age*age_vector) + (income_vector*beta_income) + (sex_vector*beta_sex) \n",
    "p = lambda x: (np.exp(beta1*x + other_covs))/(1 + np.exp(beta1*x + other_covs))\n",
    "\n",
    "# Actually calculating probabilities on exposures as defined\n",
    "px = p(xs)\n",
    "\n",
    "# Actually assigning to case or not\n",
    "outcomes = np.array([np.random.binomial(n=1,p=p,size=1)[0] for p in px])\n",
    "\n",
    "# Constructing and sorting dataframe of outcomes and exposure\n",
    "df = np.stack([outcomes,xs],axis=1)\n",
    "df = df[np.argsort(df[:, 1])]\n",
    "\n",
    "# Getting xs corresponding to each exposure interval\n",
    "categories = [df[:exposure_indices[0]]]\n",
    "for i in range(1,len(exposure_indices)):\n",
    "    categories.append(df[exposure_indices[i-1]:exposure_indices[i]])\n",
    "categories.append(df[exposure_indices[-1]:])\n",
    "\n",
    "# Getting case numbers\n",
    "cases = np.array([np.sum(df[:,0]) for df in categories])\n",
    "N = np.array([df.shape[0] for df in categories])\n",
    "non_cases = N - cases\n",
    "M = np.sum(cases)\n",
    "a_0, A = cases[0], cases[1:]\n",
    "b_0, B = non_cases[0], non_cases[1:]\n",
    "\n",
    "# # Getting crude estimates to compare\n",
    "crude = lambda x,y: (x/y) / (a_0/b_0)\n",
    "crudeors = crude(A,B)\n",
    "\n",
    "# Constructing dummy variables with respect to exposure level\n",
    "# dummy_list = []\n",
    "# for i in range(1,len(exposure_indices)):\n",
    "#     exp_cat = np.zeros(n)\n",
    "#     exp_cat[exposure_indices[i-1]:exposure_indices[i]] = 1\n",
    "#     dummy_list.append(exp_cat)\n",
    "# last_cat = np.zeros(n)\n",
    "# last_cat[exposure_indices[-1]:] = 1\n",
    "# dummy_list.append(last_cat)\n",
    "\n",
    "# # Getting regression matrix\n",
    "# to_adjust_matrix = np.column_stack([np.ones(n), age_vector, income_vector, sex_vector] + dummy_list)\n",
    "\n",
    "# # Performing logit\n",
    "# cols = [\"intercept\", \"age\", \"income\", \"sex\"]\n",
    "# cols += [\"exposure_\" + str(i+1) for i in range(len(exposure_levels))]\n",
    "# to_adjust_df = pd.DataFrame(to_adjust_matrix, columns=cols)\n",
    "to_adjust_df = pd.DataFrame({\n",
    "    'age':age_vector,\n",
    "    'income':income_vector,\n",
    "    'sex':sex_vector,\n",
    "    'exposure':np.array(exposure_levels)\n",
    "})\n",
    "to_adjust_df['intercept']=1\n",
    "\n",
    "# to_adjust_df = pd.concat([to_adjust_df,pd.get_dummies(to_adjust_df['exposure'],prefix='exposure').astype(int).iloc[:,1:]],axis=1).drop('exposure',axis=1)\n",
    "\n",
    "# results = sm.Logit(endog=df[:,0], exog=to_adjust_df).fit()\n",
    "# true_C = results.cov_params().to_numpy()[4:,4:]\n",
    "# true_C_inv = np.linalg.inv(true_C)\n",
    "\n",
    "# # Get log-ORs and corresponding variances\n",
    "# L = results.params.to_numpy()[4:]\n",
    "# v = (results.bse.to_numpy()[4:])**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 727,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3838, 16437, 26021, 33225])"
      ]
     },
     "execution_count": 727,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exposure_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 724,
   "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>age</th>\n",
       "      <th>income</th>\n",
       "      <th>sex</th>\n",
       "      <th>intercept</th>\n",
       "      <th>exposure_1</th>\n",
       "      <th>exposure_2</th>\n",
       "      <th>exposure_3</th>\n",
       "      <th>exposure_4</th>\n",
       "      <th>exposure_5</th>\n",
       "      <th>exposure_6</th>\n",
       "      <th>...</th>\n",
       "      <th>exposure_10</th>\n",
       "      <th>exposure_11</th>\n",
       "      <th>exposure_12</th>\n",
       "      <th>exposure_13</th>\n",
       "      <th>exposure_14</th>\n",
       "      <th>exposure_15</th>\n",
       "      <th>exposure_16</th>\n",
       "      <th>exposure_17</th>\n",
       "      <th>exposure_18</th>\n",
       "      <th>exposure_19</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.892860</td>\n",
       "      <td>0.870276</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.331980</td>\n",
       "      <td>0.569688</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.821229</td>\n",
       "      <td>0.746632</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.041697</td>\n",
       "      <td>0.231306</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.107657</td>\n",
       "      <td>0.319958</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39995</th>\n",
       "      <td>0.404027</td>\n",
       "      <td>0.217178</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39996</th>\n",
       "      <td>0.945486</td>\n",
       "      <td>0.739423</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39997</th>\n",
       "      <td>0.644698</td>\n",
       "      <td>0.015777</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39998</th>\n",
       "      <td>0.680413</td>\n",
       "      <td>0.701499</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39999</th>\n",
       "      <td>0.550457</td>\n",
       "      <td>0.309621</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>40000 rows × 23 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            age    income  sex  intercept  exposure_1  exposure_2  exposure_3  \\\n",
       "0      0.892860  0.870276    1          1           0           0           0   \n",
       "1      0.331980  0.569688    1          1           0           0           0   \n",
       "2      0.821229  0.746632    1          1           0           0           0   \n",
       "3      0.041697  0.231306    0          1           0           0           0   \n",
       "4      0.107657  0.319958    1          1           0           0           0   \n",
       "...         ...       ...  ...        ...         ...         ...         ...   \n",
       "39995  0.404027  0.217178    1          1           0           0           0   \n",
       "39996  0.945486  0.739423    1          1           0           0           0   \n",
       "39997  0.644698  0.015777    1          1           0           0           0   \n",
       "39998  0.680413  0.701499    0          1           0           0           0   \n",
       "39999  0.550457  0.309621    0          1           0           0           0   \n",
       "\n",
       "       exposure_4  exposure_5  exposure_6  ...  exposure_10  exposure_11  \\\n",
       "0               0           0           0  ...            0            0   \n",
       "1               0           0           0  ...            0            0   \n",
       "2               0           0           0  ...            0            0   \n",
       "3               0           0           0  ...            0            0   \n",
       "4               0           0           0  ...            0            0   \n",
       "...           ...         ...         ...  ...          ...          ...   \n",
       "39995           0           0           0  ...            0            0   \n",
       "39996           0           0           0  ...            0            0   \n",
       "39997           0           0           0  ...            0            0   \n",
       "39998           0           0           0  ...            0            0   \n",
       "39999           0           0           0  ...            0            0   \n",
       "\n",
       "       exposure_12  exposure_13  exposure_14  exposure_15  exposure_16  \\\n",
       "0                0            0            0            0            0   \n",
       "1                0            0            0            0            0   \n",
       "2                0            0            0            0            0   \n",
       "3                0            0            0            0            0   \n",
       "4                0            0            0            0            0   \n",
       "...            ...          ...          ...          ...          ...   \n",
       "39995            0            0            0            0            0   \n",
       "39996            0            0            0            0            0   \n",
       "39997            0            0            0            0            0   \n",
       "39998            0            0            0            0            0   \n",
       "39999            0            0            0            0            0   \n",
       "\n",
       "       exposure_17  exposure_18  exposure_19  \n",
       "0                0            0            0  \n",
       "1                0            0            0  \n",
       "2                0            0            0  \n",
       "3                0            0            0  \n",
       "4                0            0            0  \n",
       "...            ...          ...          ...  \n",
       "39995            0            0            1  \n",
       "39996            0            0            1  \n",
       "39997            0            0            1  \n",
       "39998            0            0            1  \n",
       "39999            0            0            1  \n",
       "\n",
       "[40000 rows x 23 columns]"
      ]
     },
     "execution_count": 724,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "to_adjust_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 720,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Cannot broadcast dimensions  (4,) (19,)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[720], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Running GL\u001b[39;00m\n\u001b[1;32m      2\u001b[0m it_fit_ex \u001b[38;5;241m=\u001b[39m IF(L,A,N,M)\n\u001b[0;32m----> 3\u001b[0m Ap, Bp, a0p, b0p \u001b[38;5;241m=\u001b[39m \u001b[43mit_fit_ex\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mconvexProgram\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m      5\u001b[0m C \u001b[38;5;241m=\u001b[39m covariance_matrix(Ap,Bp,a0p,b0p,v)\n\u001b[1;32m      6\u001b[0m inv_C \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlinalg\u001b[38;5;241m.\u001b[39minv(C)\n",
      "File \u001b[0;32m~/Documents/CorrelationCorrection/src/IterativeFitting.py:52\u001b[0m, in \u001b[0;36mIterativeFitting.convexProgram\u001b[0;34m(self, constraints, A_const, N_const, M1_const, OR)\u001b[0m\n\u001b[1;32m     47\u001b[0m constraints_eval \u001b[38;5;241m=\u001b[39m [c(A,N,M1) \u001b[38;5;28;01mfor\u001b[39;00m c \u001b[38;5;129;01min\u001b[39;00m constraints]\n\u001b[1;32m     48\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m OR:\n\u001b[1;32m     49\u001b[0m     obj \u001b[38;5;241m=\u001b[39m cp\u001b[38;5;241m.\u001b[39mMinimize(\n\u001b[1;32m     50\u001b[0m         cp\u001b[38;5;241m.\u001b[39mscalar_product(\u001b[38;5;241m-\u001b[39mL, A) \u001b[38;5;241m-\u001b[39m\n\u001b[1;32m     51\u001b[0m         cp\u001b[38;5;241m.\u001b[39mentr(M1 \u001b[38;5;241m-\u001b[39m cp\u001b[38;5;241m.\u001b[39msum(A)) \u001b[38;5;241m-\u001b[39m (M1 \u001b[38;5;241m-\u001b[39m cp\u001b[38;5;241m.\u001b[39msum(A)) \u001b[38;5;241m+\u001b[39m\n\u001b[0;32m---> 52\u001b[0m         cp\u001b[38;5;241m.\u001b[39msum(\u001b[38;5;241m-\u001b[39mcp\u001b[38;5;241m.\u001b[39mentr(\u001b[43mN\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m:\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mA\u001b[49m) \u001b[38;5;241m-\u001b[39m (N[\u001b[38;5;241m1\u001b[39m:] \u001b[38;5;241m-\u001b[39m A)) \u001b[38;5;241m+\u001b[39m\n\u001b[1;32m     53\u001b[0m         cp\u001b[38;5;241m.\u001b[39msum(\u001b[38;5;241m-\u001b[39mcp\u001b[38;5;241m.\u001b[39mentr(A) \u001b[38;5;241m-\u001b[39m A) \u001b[38;5;241m-\u001b[39m\n\u001b[1;32m     54\u001b[0m         cp\u001b[38;5;241m.\u001b[39mentr(N[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m-\u001b[39m M1 \u001b[38;5;241m+\u001b[39m cp\u001b[38;5;241m.\u001b[39msum(A)) \u001b[38;5;241m-\u001b[39m (N[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m-\u001b[39m M1 \u001b[38;5;241m+\u001b[39m cp\u001b[38;5;241m.\u001b[39msum(A)) \n\u001b[1;32m     55\u001b[0m     )\n\u001b[1;32m     56\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m     57\u001b[0m     log_N \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlog(N[\u001b[38;5;241m1\u001b[39m:])\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:777\u001b[0m, in \u001b[0;36mExpression.__array_ufunc__\u001b[0;34m(self, ufunc, method, *args, **kwargs)\u001b[0m\n\u001b[1;32m    773\u001b[0m ufunc_handler \u001b[38;5;241m=\u001b[39m __BINARY_EXPRESSION_UFUNCS__[ufunc]\n\u001b[1;32m    774\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m kwargs \u001b[38;5;241m==\u001b[39m {} \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    775\u001b[0m         \u001b[38;5;28mlen\u001b[39m(args) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    776\u001b[0m         args[\u001b[38;5;241m1\u001b[39m] \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28mself\u001b[39m:\n\u001b[0;32m--> 777\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mufunc_handler\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    778\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m kwargs\u001b[38;5;241m.\u001b[39mkeys() \u001b[38;5;241m==\u001b[39m {\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mout\u001b[39m\u001b[38;5;124m'\u001b[39m} \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    779\u001b[0m         \u001b[38;5;28mlen\u001b[39m(args) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    780\u001b[0m         args[\u001b[38;5;241m1\u001b[39m] \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28mself\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    781\u001b[0m         \u001b[38;5;28misinstance\u001b[39m(kwargs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mout\u001b[39m\u001b[38;5;124m'\u001b[39m], \u001b[38;5;28mtuple\u001b[39m) \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    782\u001b[0m         \u001b[38;5;28mlen\u001b[39m(kwargs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mout\u001b[39m\u001b[38;5;124m'\u001b[39m]) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m1\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \\\n\u001b[1;32m    783\u001b[0m         args[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;129;01mis\u001b[39;00m kwargs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mout\u001b[39m\u001b[38;5;124m'\u001b[39m][\u001b[38;5;241m0\u001b[39m]:\n\u001b[1;32m    784\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(__INPLACE_MUTATION_ERROR__)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:86\u001b[0m, in \u001b[0;36m<lambda>\u001b[0;34m(self, a)\u001b[0m\n\u001b[1;32m     71\u001b[0m __INPLACE_MUTATION_ERROR__ \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\"\"\u001b[39m\n\u001b[1;32m     72\u001b[0m \u001b[38;5;124mYou\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mre trying to mutate a CVXPY expression inplace. This is prone to errors or\u001b[39m\n\u001b[1;32m     73\u001b[0m \u001b[38;5;124mcode that doesn\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mt behave as expected. Consider alternatives. For example, replace\u001b[39m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     76\u001b[0m \u001b[38;5;124m    x = x + 1\u001b[39m\n\u001b[1;32m     77\u001b[0m \u001b[38;5;124m\"\"\"\u001b[39m\n\u001b[1;32m     79\u001b[0m __ABS_ERROR__ \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\"\"\u001b[39m\n\u001b[1;32m     80\u001b[0m \u001b[38;5;124mYou\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mre calling the built-in abs function on a CVXPY expression. This is not\u001b[39m\n\u001b[1;32m     81\u001b[0m \u001b[38;5;124msupported. Consider using the abs function provided by CVXPY.\u001b[39m\n\u001b[1;32m     82\u001b[0m \u001b[38;5;124m\"\"\"\u001b[39m\n\u001b[1;32m     84\u001b[0m __BINARY_EXPRESSION_UFUNCS__ \u001b[38;5;241m=\u001b[39m {\n\u001b[1;32m     85\u001b[0m         np\u001b[38;5;241m.\u001b[39madd: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__radd__\u001b[39m(a),\n\u001b[0;32m---> 86\u001b[0m         np\u001b[38;5;241m.\u001b[39msubtract: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__rsub__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m)\u001b[49m,\n\u001b[1;32m     87\u001b[0m         np\u001b[38;5;241m.\u001b[39mmultiply: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__rmul__\u001b[39m(a),\n\u001b[1;32m     88\u001b[0m         np\u001b[38;5;241m.\u001b[39mdivide: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m__rdiv__(a),\n\u001b[1;32m     89\u001b[0m         np\u001b[38;5;241m.\u001b[39mmatmul: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__rmatmul__\u001b[39m(a),\n\u001b[1;32m     90\u001b[0m         np\u001b[38;5;241m.\u001b[39mpower: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__rpow__\u001b[39m(a),\n\u001b[1;32m     91\u001b[0m         np\u001b[38;5;241m.\u001b[39mleft_shift: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__rlshift__\u001b[39m(a),\n\u001b[1;32m     92\u001b[0m         np\u001b[38;5;241m.\u001b[39mright_shift: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__rrshift__\u001b[39m(a),\n\u001b[1;32m     93\u001b[0m         np\u001b[38;5;241m.\u001b[39mequal: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__eq__\u001b[39m(a),\n\u001b[1;32m     94\u001b[0m         \u001b[38;5;66;03m# <= and >= are backwards because this is only called for code of the\u001b[39;00m\n\u001b[1;32m     95\u001b[0m         \u001b[38;5;66;03m# form ndarray <= Expression\u001b[39;00m\n\u001b[1;32m     96\u001b[0m         np\u001b[38;5;241m.\u001b[39mless_equal: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__ge__\u001b[39m(a),\n\u001b[1;32m     97\u001b[0m         np\u001b[38;5;241m.\u001b[39mgreater_equal: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__le__\u001b[39m(a),\n\u001b[1;32m     98\u001b[0m         np\u001b[38;5;241m.\u001b[39mless: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__gt__\u001b[39m(a),\n\u001b[1;32m     99\u001b[0m         np\u001b[38;5;241m.\u001b[39mgreater: \u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;28mself\u001b[39m, a: \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__lt__\u001b[39m(a),\n\u001b[1;32m    100\u001b[0m }\n\u001b[1;32m    103\u001b[0m ExpressionLike \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mExpression | np.typing.ArrayLike\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    106\u001b[0m \u001b[38;5;28;01mclass\u001b[39;00m \u001b[38;5;21;01mExpression\u001b[39;00m(u\u001b[38;5;241m.\u001b[39mCanonical):\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:50\u001b[0m, in \u001b[0;36m_cast_other.<locals>.cast_op\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"A wrapped binary operator that can handle non-Expression arguments.\u001b[39;00m\n\u001b[1;32m     48\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m     49\u001b[0m other \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcast_to_const(other)\n\u001b[0;32m---> 50\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mbinary_op\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mother\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:619\u001b[0m, in \u001b[0;36mExpression.__rsub__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    615\u001b[0m \u001b[38;5;129m@_cast_other\u001b[39m\n\u001b[1;32m    616\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__rsub__\u001b[39m(\u001b[38;5;28mself\u001b[39m, other: ExpressionLike) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mExpression\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m    617\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Expression : The difference of two expressions.\u001b[39;00m\n\u001b[1;32m    618\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 619\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mother\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:50\u001b[0m, in \u001b[0;36m_cast_other.<locals>.cast_op\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"A wrapped binary operator that can handle non-Expression arguments.\u001b[39;00m\n\u001b[1;32m     48\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m     49\u001b[0m other \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcast_to_const(other)\n\u001b[0;32m---> 50\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mbinary_op\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mother\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:613\u001b[0m, in \u001b[0;36mExpression.__sub__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    609\u001b[0m \u001b[38;5;129m@_cast_other\u001b[39m\n\u001b[1;32m    610\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__sub__\u001b[39m(\u001b[38;5;28mself\u001b[39m, other: ExpressionLike) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mExpression\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m    611\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Expression : The difference of two expressions.\u001b[39;00m\n\u001b[1;32m    612\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 613\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43mother\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:50\u001b[0m, in \u001b[0;36m_cast_other.<locals>.cast_op\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"A wrapped binary operator that can handle non-Expression arguments.\u001b[39;00m\n\u001b[1;32m     48\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m     49\u001b[0m other \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcast_to_const(other)\n\u001b[0;32m---> 50\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mbinary_op\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mother\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/expressions/expression.py:599\u001b[0m, in \u001b[0;36mExpression.__add__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    597\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\n\u001b[1;32m    598\u001b[0m \u001b[38;5;28mself\u001b[39m, other \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbroadcast(\u001b[38;5;28mself\u001b[39m, other)\n\u001b[0;32m--> 599\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mcvxtypes\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_expr\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mother\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/atoms/affine/add_expr.py:33\u001b[0m, in \u001b[0;36mAddExpression.__init__\u001b[0;34m(self, arg_groups)\u001b[0m\n\u001b[1;32m     32\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, arg_groups: Iterable[Expression]) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m---> 33\u001b[0m     \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mAddExpression\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43marg_groups\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     34\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39margs \u001b[38;5;241m=\u001b[39m []\n\u001b[1;32m     35\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m group \u001b[38;5;129;01min\u001b[39;00m arg_groups:\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/atoms/atom.py:51\u001b[0m, in \u001b[0;36mAtom.__init__\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39margs \u001b[38;5;241m=\u001b[39m [Atom\u001b[38;5;241m.\u001b[39mcast_to_const(arg) \u001b[38;5;28;01mfor\u001b[39;00m arg \u001b[38;5;129;01min\u001b[39;00m args]\n\u001b[1;32m     50\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvalidate_arguments()\n\u001b[0;32m---> 51\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_shape \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape_from_args\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     52\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_shape) \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m     53\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mAtoms must be at most 2D.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/atoms/affine/add_expr.py:41\u001b[0m, in \u001b[0;36mAddExpression.shape_from_args\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     38\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mshape_from_args\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m Tuple[\u001b[38;5;28mint\u001b[39m, \u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m]:\n\u001b[1;32m     39\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Returns the (row, col) shape of the expression.\u001b[39;00m\n\u001b[1;32m     40\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m---> 41\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mu\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msum_shapes\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[43marg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43marg\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/utilities/shape.py:50\u001b[0m, in \u001b[0;36msum_shapes\u001b[0;34m(shapes)\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m t \u001b[38;5;129;01min\u001b[39;00m shapes[\u001b[38;5;241m1\u001b[39m:]:\n\u001b[1;32m     48\u001b[0m     \u001b[38;5;66;03m# Only allow broadcasting for 0D arrays or summation of scalars.\u001b[39;00m\n\u001b[1;32m     49\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m shape \u001b[38;5;241m!=\u001b[39m t \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(squeezed(shape)) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(squeezed(t)) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m---> 50\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m     51\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCannot broadcast dimensions \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m\n\u001b[1;32m     52\u001b[0m             \u001b[38;5;28mlen\u001b[39m(shapes)\u001b[38;5;241m*\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m \u001b[38;5;28mtuple\u001b[39m(shapes))\n\u001b[1;32m     54\u001b[0m     longer \u001b[38;5;241m=\u001b[39m shape \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(shape) \u001b[38;5;241m>\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlen\u001b[39m(t) \u001b[38;5;28;01melse\u001b[39;00m t\n\u001b[1;32m     55\u001b[0m     shorter \u001b[38;5;241m=\u001b[39m shape \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(shape) \u001b[38;5;241m<\u001b[39m \u001b[38;5;28mlen\u001b[39m(t) \u001b[38;5;28;01melse\u001b[39;00m t\n",
      "\u001b[0;31mValueError\u001b[0m: Cannot broadcast dimensions  (4,) (19,)"
     ]
    }
   ],
   "source": [
    "# Running GL\n",
    "it_fit_ex = IF(L,A,N,M)\n",
    "Ap, Bp, a0p, b0p = it_fit_ex.convexProgram()\n",
    "\n",
    "C = covariance_matrix(Ap,Bp,a0p,b0p,v)\n",
    "inv_C = np.linalg.inv(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 711,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recovering slope estimate for corrected correlation\n",
    "x = np.array(exposure_levels)\n",
    "vb_star = 1/(np.dot(x,np.dot(inv_C,x)))\n",
    "b_star = vb_star*(np.dot(x,np.dot(inv_C,L)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 712,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01948094498080707"
      ]
     },
     "execution_count": 712,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 713,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026950187506593885"
      ]
     },
     "execution_count": 713,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Recovering standard slope estimate\n",
    "vb = 1/(np.dot(x,np.dot(np.linalg.inv(np.diag(v)),x)))\n",
    "b = vb*(np.dot(x,np.dot(np.linalg.inv(np.diag(v)),L)))\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 714,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.016708103411178617"
      ]
     },
     "execution_count": 714,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rootWols = np.diag(1/np.sqrt(v))\n",
    "X_intercept_ols = np.vstack([np.ones(len(x)),x]).T\n",
    "intercept_ols_slope = np.linalg.lstsq(rootWols@X_intercept_ols,rootWols@L,rcond = 1e-15)[0][1]\n",
    "intercept_ols_slope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 715,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01946475967133422"
      ]
     },
     "execution_count": 715,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Getting Oracle estimates\n",
    "vb_oracle = 1/(np.dot(x,np.dot(true_C_inv,x)))\n",
    "b_oracle = vb_oracle*(np.dot(x,np.dot(true_C_inv,L)))\n",
    "b_oracle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 716,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01947493896280948"
      ]
     },
     "execution_count": 716,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Getting estimates from covariance matrix generated from raw counts\n",
    "C_raw = covariance_matrix(A,B,a_0,b_0,v)\n",
    "inv_C_raw = np.linalg.inv(C_raw)\n",
    "\n",
    "vb_r = 1/(x.T@inv_C_raw@x)\n",
    "b_r = vb_r*(x.T@inv_C_raw@L)\n",
    "b_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 717,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x30c0cbd10>"
      ]
     },
     "execution_count": 717,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAACbqElEQVR4nOzdd3xV9f3H8de5OzfJvdlkQhISBET2UFDA9cPWUayzWlcdba1aS4OIA0Flx9ZaV7XWVQeu4saBshGQpewkBLLJvjfJ3fee3x83hIQhQ5KbhM/z8eDxIDfn3vu5JOR+8j3f9/koqqqqCCGEEEKEiCbUBQghhBDi1CbNiBBCCCFCSpoRIYQQQoSUNCNCCCGECClpRoQQQggRUtKMCCGEECKkpBkRQgghREhJMyKEEEKIkNKFuoBjEQgEKCsrIzIyEkVRQl2OEEIIIY6Bqqo0NDSQnJyMRnPk9Y8u0YyUlZWRlpYW6jKEEEIIcQKKi4tJTU094ue7RDMSGRkJBF+MxWIJcTVCCCGEOBZ2u520tLSW9/Ej6RLNyP5TMxaLRZoRIYQQoos52hYL2cAqhBBCiJCSZkQIIYQQISXNiBBCCCFCqkvsGTkWfr8fr9cb6jJEB9Pr9Wi12lCXIYQQ4mfoFs1IY2MjJSUlqKoa6lJEB1MUhdTUVCIiIkJdihBCiBPU5ZsRv99PSUkJZrOZ+Ph4uSjaKURVVaqqqigpKSE7O1tWSIQQoovq8s2I1+tFVVXi4+MJCwsLdTmig8XHx7Nnzx68Xq80I0II0UV1mw2ssiJyapKvuxBCdH3dphkRQgghRNckzYgQQgghQkqaESGEEEKElDQjIVZRUcGf//xnsrKyMJlM9OjRgzFjxvDcc8/hcDgASE9P58knnwxtoUIIIUQ76fJpmq5s9+7djBkzhqioKGbNmsUZZ5yB0Wjkxx9/5IUXXiAlJYXLLrss1GUKIYToxr7e+zWf7v6U3HG5aDWhSSV222bE4fEd8XMaRcGk157UY82G4/+nvPPOO9HpdHz//feEh4e33J6ZmcmvfvUruYibEEKIduP2u3ni+yd4a8dbACzMX8gVfa4ISS3dthnpP+2LI37u3NPiefmWkS0fD3vsa5xe/2GPHZURw4Lfn9Xy8dlzv6W2yXPIcXvmXHxc9dXU1PDll18ya9asNo1IaxJbFUII0R6K7EXkLM1he+12AG4ZcAuXZYVuJV72jIRIfn4+qqpy2mmntbk9Li6OiIgIIiIimDJlSoiqE0II0V0tKlzE1Z9czfba7UQbo3n2/GeZNGwSeo0+ZDV125WRbY9OOOLnNAetOKx/+IJjPnbFlHN/XmFHsXbtWgKBANdffz1ut7tdn0sIIcSpw+VzMXfdXN7b9R4AQxOGMm/sPHqE9whxZd24GTmePRztdexPycrKQlEUdu7c2eb2zMxMALm0vRBCiJOm0FZIztIcdtXtQkHh9oG388dBf0Sn6RxtgJymCZHY2FguvPBCnn76aZqamkJdjhBCiG7q44KPueaTa9hVt4sYUwzPX/g8dw+5u9M0ItCNV0a6gmeffZYxY8YwfPhwpk+fzsCBA9FoNKxbt44dO3YwbNiwlmNLS0vZtGlTm/v36tWL6OjoDq5aCCFEV+D0OZm1ZhYL8xcCMDJxJHPOmUO8OT60hR2GonaB/KjdbsdqtWKz2bBYLG0+53K5KCwsJCMjA5PJFKIKT1x5eTmzZs3i008/paSkBKPRSP/+/bnqqqu48847MZvNpKens3fv3kPu+/rrr/Pb3/42BFV3Hl396y+EEO0hvy6fnKU5FNgKUFD446A/csfAOzr8OiI/9f7dmqyMhFhSUhL//Oc/+ec//3nEY/bs2dNxBQkhhOiyVFVlYf5CZq2ZhcvvIi4sjrnnzGVk0sij3zmEpBkRQgghugGH18Hj3z3Ox7s/BuCspLOYfc5sYsNiQ1zZ0UkzIoQQQnRxO2t3krM0hz32PWgUDXcNvotbz7gVjdI1cirSjAghhBBdlKqqvJf3HnPXzsXtd5NgTmDe2HkM6zHs6HfuRKQZEUIIIbqgRk8jj65+lM/3fA7AOSnnMPPsmUSbul7KUpoRIYQQoovZXrOdnKU5FDUUoVN03DP0Hm46/aYuc1rmYNKMCCGEEF2Eqqq8vfNt5q+bjzfgJSk8iXlj5zE4YXCoS/tZpBkRQgghugC7x870VdP5au9XAIxPG8/jYx7HarSGuLKfT5oRIYQQopPbUr2FnKU5lDaWotPomDRsEr/t91uUg4a5dlXSjAghhBCdlKqq/Hf7f/nb+r/hC/hIiUghd1wuA+IGhLq0k6pr7nQRJ0V6ejpPPvlkqMsQQghxGDa3jXu+vYd56+bhC/i4sNeFvHPpO92uEQFpRkKuuLiY3/3udyQnJ2MwGOjVqxd//vOfqampCXVpx83j8TBv3jwGDRqE2WwmLi6OMWPG8PLLL+P1ekNdnhBCdBmbKjdx1cdXsaR4CXqNngdGPcAT457AYjjyfJeurHudplmWC2og1FWAooGxOUc9bPfu3Zx11ln06dOHt956i4yMDLZu3crkyZP5/PPP+e6774iJiTnkfh6PB4PB0B6VnzCPx8OECRPYvHkzjz32GGPGjMFisfDdd9+Rm5vLkCFDGDx4cKjLFEKITi2gBnh166s8teEpfKqPtMg0csfl0j+2f6hLa1fda2VEDUDAH/o/x9gQ/elPf8JgMPDll18ybtw4evbsyS9+8Qu+/vprSktLefDBB4Hg6ZTHHnuMG2+8EYvFwh133AHAlClT6NOnD2azmczMTB5++OFDViA+/vhjRowYgclkIi4ujssvv/yI9dTX13PbbbcRHx+PxWLhvPPOY/Pmzcf0Wp588kmWLVvG4sWL+dOf/sTgwYPJzMzkuuuuY82aNWRnZwPgdru55557SEhIwGQycfbZZ7Nu3ToAAoEAqampPPfcc20ee+PGjWg0msNOLhZCiO6izlXHXYvvCu4PUX1clH4R71zyTrdvRKC7NSNdSG1tLV988QV33nknYWFhbT6XmJjI9ddfz4IFC1BVFYDc3FwGDRrExo0befjhhwGIjIzklVdeYdu2bfzjH//gxRdf5O9//3vL43z66adcfvnl/PKXv2Tjxo0sXryYkSOPPLnxqquuorKyks8//5z169czdOhQzj//fGpra4/6et544w0uuOAChgwZcsjn9Ho94eHhANx33328//77vPrqq2zYsIGsrCwmTJhAbW0tGo2G3/zmN7z55puHPPaYMWPo1avXUesQQnQNqqricflwNLjxuHwtP+tOVev3refKj69keelyjFoj086axryx84gwRIS6tA7RvU7TdCF5eXmoqkq/fv0O+/l+/fpRV1dHVVUVAOeddx5//etf2xzz0EMPtfw9PT2dnJwc3n77be677z4AZs6cybXXXsuMGTNajhs0aNBhn2/FihWsXbuWyspKjEYjEGyAFi5cyHvvvdeyGvNTr2f8+PE/eUxTUxPPPfccr7zyCr/4xS8AePHFF/nqq6946aWXmDx5Mtdffz1PPPEERUVF9OzZk0AgwNtvv93mtQohui6v20fxrhryN1ZQWWIj4FfRaBUSUq1kDUkkrU8seuOp89YUUAO89ONLPLPpGfyqn3RLOrnjcjkt5rRQl9ahTp2veCd1rL8NDB8+/JDbFixYwFNPPUVBQQGNjY34fD4slgObmzZt2sTtt99+TI+/efNmGhsbiY1tO2ra6XRSUFBw1Psfy+soKCjA6/UyZsyYltv0ej0jR45k+/btAAwePJh+/frx5ptvcv/997N06VIqKyu56qqrjul1CCE6r8piGys/3EFlsR1FoxBhMaIzaQn4AuzZXkXh1koS0iyM+VVfEtK6/oW8jqbGWcPU5VNZXb4agEszL+WhMx/CrDeHuLKOJ81IiGRlZaEoCtu3bz/sPo7t27cTHR1NfHw8QMtpjv1Wr17N9ddfz4wZM5gwYQJWq5W3336bJ554ouWYg0///JTGxkaSkpJYsmTJIZ+Lioo66v379OnDjh07jvn5fsr111/f0oy8+eabXHTRRYc0SUKIrqWy2Mbit7Zgr3GQkGZBp9e2+XxkTBg+r5/KYjvfvL2F864d0K0bkrXla5myfArVzmpMWhMPjHqAiVkTu81FzI5X99ozomhAow39n2MYVBQbG8uFF17Is88+i9PpbPO5iooK3njjDa655pojfmOuWrWKXr168eCDDzJ8+HCys7MP2eA5cOBAFi9efEz/dEOHDqWiogKdTkdWVlabP3FxcUe9/3XXXcfXX3/Nxo0bD/mc1+ulqamJ3r17YzAYWLlyZZvPrVu3jv79+7d5rC1btrB+/Xree+89rr/++mN6DUKIzsnr9rHywx3YaxwkZUQd0ojsp9NrScqIwlbtYOWHO/C6fR1cafvzB/w8t+k5bv/qdqqd1fS29uati9/i8uzLT9lGBLrbysgxxGk7k6effprRo0czYcIEHn/88TbR3pSUFGbOnHnE+2ZnZ1NUVMTbb7/NiBEj+PTTT/nf//7X5phHHnmE888/n969e3Pttdfi8/n47LPPmDJlyiGPd8EFF3DWWWcxceJE5s2bR58+fSgrK2vZBHu400St3XvvvXz66aecf/75PPbYY5x99tlERkby/fffM3fuXF566SUGDx7MH//4RyZPnkxMTAw9e/Zk3rx5OBwObr311pbHSk9PZ/To0dx66634/X4uu+yy4/yXFUJ0JsW7aqgstpOQZjnqG66iKCSkWqgstlOSV0vGgIQOqrL9VTmqmLp8Kmsq1gBwedblTB01lTDdsa9id1fda2Wki8nOzub7778nMzOTq6++mt69e3PHHXdw7rnnsnr16sNeY2S/yy67jL/85S/cddddDB48mFWrVrWkbPYbP3487777Lh999BGDBw/mvPPOY+3atYd9PEVR+Oyzzxg7diy33HILffr04dprr2Xv3r306NHjqK/FaDTy1Vdfcd999/Gvf/2LM888kxEjRvDUU09xzz33MGBA8IqBc+bM4YorruCGG25g6NCh5Ofn88UXXxAdHd3m8a6//no2b97M5Zdfflynm4QQnYuqquRvrEDRKEdcETmYzqBFURTyNpR3m5TNqrJVXPnxlaypWEOYLoxZZ8/i0TGPSiPSTFG7wFfabrdjtVqx2WxtNmgCuFwuCgsLycjIwGQyhahCESry9Reic/O4fLzzt1VoFIXImGN/422odRJQVa6eNBqDqesu4vsCPp7d9Cz//vHfqKhkR2eTOy6XTGtmqEvrED/1/t1a1/0KCyGE6PR8Xj8Bv4rOdGyrIvtpdBp8Lh8+r7/LNiMVTRVMWTaFDZUbALiqz1XcN+I+TDr5xelgcppGHJPTTz+diIiIw/554403Ql2eEKKT0um1aLQKAd/xjeoI+AJotMd+aqezWV6ynKs+vooNlRsI14czb+w8pp01TRqRI+ia7abocJ999tkRh90dy54SIcSpSW/UkpBqZc/2quM6TdNoc5PePx69sWs1I96Al39u/Ccvb3kZgH4x/Zg/bj69LHIF6Z9yQisjzzzzDOnp6ZhMJkaNGnXETZEHe/vtt1EUhYkTJ57I04oQ6tWr1yGR3/1/IiMjQ12eEKKTUhSFrCGJqAEVn9d/TPfxefyoqkr20KQuFXctbyznlkW3tDQiv+n7G17/5evSiByD425GFixYwKRJk3jkkUfYsGEDgwYNYsKECVRWVv7k/fbs2UNOTg7nnHPOCRcrhBCi60nrE0tCWjCue7TMhKqqVJYEY8Cp2UdOFHY23xZ9y5UfX8nmqs1E6iP52/i/8cCoBzBqjaEurUs47mbkb3/7G7fffju33HIL/fv35/nnn8dsNvOf//zniPfx+/0tVwvNzDw1dhALIYQI0ht1jPlVX6xxZsoL6/F5Dr9C4vP4KS+sxxpn5uyJ/brEjBqv38u8dfO459t7sHvsDIgdwIJLF3BhrwtDXVqXclxfaY/Hw/r165k6dWrLbRqNhgsuuIDVq1cf8X6PPvooCQkJ3HrrrSxfvvyoz+N2u3G73S0f2+324ylTCCFEJ5OQZuW8awccmE2jKERYjWh0GgK+AI02N6qqkpBm4eyJ/YhPPXIMtLMoaSjhvmX38WP1jwD8tt9vmTRsEnqtPsSVdT3H1YxUV1fj9/sP2bDYo0ePI84lWbFiBS+99BKbNm065ueZPXt2m0mzQgghur6ENCuX3D6Mkrxa8jaUU1liw+fyodEqpPePJ3toEqnZMV1iReTrvV8zbeU0GrwNWAwWHh/zOOf2PDfUZXVZ7foVb2ho4IYbbuDFF188pvkm+02dOpVJkya1fGy320lLS2uPEoUQQnQgvVFHxoAE0k+Px+v24/P60em16I3aLrFZ1eP3kPt9Lm/teAuAgfEDmT92PskRySGurGs7rmYkLi4OrVbLvn372ty+b98+EhMTDzm+oKCAPXv2cOmll7bcFggEs+Y6nY6dO3fSu3fvQ+5nNBoxGmXTz8mwZMkSzj33XOrq6o5p+q4QQnQERVEwmHRd6oJmRfYicpbmsL12OwC3DLiFu4fcjV4jp2V+ruPawGowGBg2bFibSbCBQIDFixdz1llnHXJ83759+fHHH9m0aVPLn8suu4xzzz2XTZs2nbKrHYqi/OSf6dOnn9Djjh8/nnvvvfek1iqEEAIWFS7i6k+uZnvtdqKMUTxz/jPB/SHSiJwUx92STpo0iZtuuonhw4czcuRInnzySZqamrjlllsAuPHGG0lJSWH27NmYTKaWAWn77f/t/ODbTyXl5eUtf1+wYAHTpk1j586dLbdFRES0/F1VVfx+Pzpd1/ntQQghuguXz8W8dfN4d9e7AAxNGMrcsXNJDD/0bIA4cccd7b3mmmvIzc1l2rRpDB48mE2bNrFo0aKWTa1FRUVt3mxDxe1Vj/uPP3Ag/+4PBG/z+NRjetzjkZiY2PLHarWiKErLxzt27CAyMpLPP/+cYcOGYTQaWbFiBTfffPMhF4u79957GT9+PAA333wzS5cu5R//+EfLCsuePXtajl2/fj3Dhw/HbDYzevToNs2PEEKIQxXaCrn+s+t5d9e7KCjcfsbtvDThJWlE2sEJ/bp91113cddddx32c0uWLPnJ+77yyisn8pTH7a4X6477Pr//vwiGZxkA2Ljby7++bKRPso7JEw9EzO5/vZ5G16HNx4t3ntyL89x///3k5uaSmZlJdHT0UY//xz/+wa5duxgwYACPPvooAPHx8S0NyYMPPsgTTzxBfHw8f/jDH/jd737HypUrT2rNQgjRXXxc8DGPffcYTp+TGFMMs8+ezeiU0aEu66RTVbVTbCSWtf9O6tFHH+XCC4/9ojlWqxWDwYDZbD7sZuKZM2cybtw4INjoXHzxxbhcLkwmGdokhBD7OX1OZq+Zzf/y/wfAyMSRzDlnDvHm+BBXdnJ53T6Kd9WQv7GCyhIbAb+KRquQkGola0giaX1iOzRi3W2bkadvP/pqwsF0reYxDcnU8/Tt0RzcIM65IernFXaMhg8fflIfb+DAgS1/T0pKAqCyspKePXue1OcRQoiuqqC+gJylOeTX56Og8IdBf+D3A3+PVtO1hvUdTWWx7cDF5zQKERYjOpOWgC/Anu1VFG6tJCHNwphf9SUhzdohNXXbZsSo/3nLTFqNgvYwO2p+7uMeq/Dw8DYfazSaQ2Y6HGmK7uHo9Qd2fO9fgtsfsxZCiFOZqqoszF/IrDWzcPldxIXFMeecOYxKGhXq0k66ymIbi9/agr3GQUKaBZ2+baMVGROGz+unstjON29v4bxrB3RIQ3JCU3tFx4uPjz9kY/DBV7U1GAz4/cc2FVMIIQQ4vA4eXPEg01ZNw+V3cVbSWbx76bvdshHxun2s/HAH9hoHSRlRhzQi++n0WpIyorBVO1j54Q68bl+71ybNSBdx3nnn8f333/Paa6+Rl5fHI488wpYtW9ock56ezpo1a9izZw/V1dWy8iGEED9hZ+1Orv30Wj7e/TEaRcPdQ+7m+QufJy7s2K8Y3pUU76qhsjg4Eflom1QVRSEhNThpuSSvtt1rk2aki5gwYQIPP/ww9913HyNGjKChoYEbb7yxzTE5OTlotVr69+9PfHw8RUVFIapWCCE6L1VVeXfXu1z/2fUU2gpJCEvgpf97iTsG3oFG6Z5vi6qqkr+xAkWjHHFF5GA6QzBZk7eh/JBtAieborb3M5wEdrsdq9WKzWbDYmk7ydHlclFYWEhGRoYkQ05B8vUXQhyPRk8jj65+lM/3fA7A2SlnM/PsmcSYTu7lGTobj8vHO39bhUZRiIwJO+b7NdQ6CagqV08afUKX7v+p9+/Wuu0GViGEEKK17TXbyVmaQ1FDEVpFyz1D7+Hm02/utqshrfm8fgJ+FZ3p+JJBGp0Gn8uHz+tv1zlC0owIIYTo1lRVZcHOBcxbNw9vwEtieCLzx85ncMLgUJfWYXR6LRqtQsB3fHsJA74AGu2xn9o5UdKMCCGE6LYaPA08suoRvtr7FQDjU8fz+NmPYzV2zPUzOgu9UUtCqpU926uO6zRNo81Nev949EZpRoQQQojjtqV6CzlLcyhtLEWn0fGXoX/hhv43hORy56GmKApZQxIp3FrZcun3o/F5/KiqSvbQpHb/N5NmRAghRLeiqir/3f5f/rb+b/gCPlIiUpg/dj5nxJ8R6tJCKq1PLAlpwbhuUkbUTzYYqqpSWRKMAadmt//m3u6/a0cIIcQpw+a28edv/8y8dfPwBXxc0PMC3rn0nVO+EQHQG3WM+VVfrHFmygvr8XkOf5FMn8dPeWE91jgzZ0/s1yEzamRlRAghRLewuWozk5dOprypHL1GT87wHH7T9zen5GmZI0lIs3LetQMOzKZRFCKsRjQ6DQFfgEabG1VVSUizcPbEfsSnHjmOezJJMyKEEKJLC6gBXt36Kk9teAqf6iMtMo3ccbn0j+0f6tI6pYQ0K5fcPoySvFryNpRTWWLD5/Kh0Sqk948ne2gSqdkxMrVXtI/p06ezcOHCQ2baCCFEV1XnquOhlQ+xrGQZABelX8QjZz1ChCEixJV1bnqjjowBCaSfHo/X7W/Z1Ko3akOykiR7RkKsoqKCu+++m8zMTIxGI2lpaVx66aUsXrw41KUJIUSntmHfBq78+EqWlSzDoDHw8JkPM2/sPGlEjoOiKBhMOsyRRgwmXchOacnKSAjt2bOHMWPGEBUVxfz58znjjDPwer188cUX/OlPf2LHjh2H3Mfr9aLX60NQrRBCdA4BNcB/tvyHpzc+jV/1k25JJ3dcLqfFnBbq0sQJ6lbNyPJZswj4D787uCNptFrOeeCBox535513oigKa9euJTw8vOX2008/nd/97ndAsGt99tln+fzzz1m8eDGTJ08mPT2de++9l/r6+pb7LFy4kMsvv7zNMKM5c+bw97//HYfDwdVXX018fPwhNfz73//miSeeoLCwkPT0dO655x7uvPPOn/HqhRCi/dQ4a3hgxQOsKlsFwCWZl/DwmQ9j1ptDXJn4ObpVMxLw+1E7QTNyLBfbra2tZdGiRcycObNNI7JfVFRUy9+nT5/OnDlzePLJJ9HpdHzzzTdHffx33nmH6dOn88wzz3D22Wfz+uuv89RTT5GZmdlyzBtvvMG0adN4+umnGTJkCBs3buT2228nPDycm2666VheqhBCdJh1FeuYsmwKVc4qTFoTD4x6gIlZEyUt0w10q2akK8nPz0dVVfr27XvUY6+77jpuueWW43r8J598kltvvZVbb70VgMcff5yvv/4al8vVcswjjzzCE088wa9//WsAMjIy2LZtG//617+kGRFCdBr+gJ8XfniB5394noAaINOayRPjniArOivUpYmTRJqREGl9OuVohg8fftyPv337dv7whz+0ue2ss87i22+/BaCpqYmCggJuvfVWbr/99pZjfD4fVuupNbNBCNF5VTuruX/Z/aypWAPAxKyJTB05VU7LdDPSjIRIdnY2iqIcdpPqwQ4+jaPRaA5pZrxe73E9f2NjIwAvvvgio0aNavM5rbZ9ByIJIcSxWF22mvuX30+tq5YwXRgPn/kwl/a+NNRliXYg0d4QiYmJYcKECTzzzDM0NTUd8vnWm1MPFh8fT0NDQ5v7HXztkH79+rFmzZo2t3333Xctf+/RowfJycns3r2brKysNn8yMjJO7EUJIcRJ4Av4eGrDU/z+q99T66olOzqbty95WxqRbqxbrYxotNpj2jza3jTHuLLwzDPPMGbMGEaOHMmjjz7KwIED8fl8fPXVVzz33HNs3779sPcbNWoUZrOZBx54gHvuuYc1a9bwyiuvtDnmz3/+MzfffDPDhw9nzJgxvPHGG2zdurXNBtYZM2Zwzz33YLVaueiii3C73Xz//ffU1dUxadKkE379QghxovY17WPK8ims37cegCv7XMmUEVMw6Uwhrky0p27VjBxLnLYzyczMZMOGDcycOZO//vWvlJeXEx8fz7Bhw3juueeOeL+YmBj++9//MnnyZF588UXOP/98pk+fzh133NFyzDXXXENBQQH33XcfLpeLK664gj/+8Y988cUXLcfcdtttmM1m5s+fz+TJkwkPD+eMM87g3nvvbc+XLYQQh7WidAUPLH+AOncdZp2Z6aOn84uMX4S6LNEBFPV4dlKGiN1ux2q1YrPZsFjaDu1xuVwUFhaSkZGBySSd86lGvv5CdH3egJenNz7Nf7b8B4C+MX3JHZdLL0uvEFcmfq6fev9urVutjAghhOhayhvLuW/ZfWyq2gTAtaddS86IHIxaY2gLEx1KmhEhhBAhsaR4CQ+tfAib20aEPoIZo2fwf+n/F+qyRAhIMyKEEKJDef1entzwJK9tew2A02NPZ/64+aRFpoW4MhEq0owIIYToMCUNJdy37D5+rP4RgN/2+y1/GfYXDFpDiCsToSTNiBBCiA7x9d6vmbZyGg3eBiINkTw+5nHO63leqMsSnYA0I0IIIdqVx+8h9/tc3trxFgAD4wcyf+x8kiOSQ1yZ6CykGRFCCNFuiuxF5CzNYXtt8CKOt5x+C3cPvRu9Rh/iykRnIs2IEEKIdrFozyKmr5pOk7eJKGMUM8+eydjUsaEuS3RC0owIIYQ4qVw+F/PXzeedXe8AMDRhKHPHziUxPDHElYnOSpoRIYQQJ02hrZCcpTnsqtsFwG1n3MafBv8JnUbebsSRydTeELn55puZOHHiMR+vKAoLFy5st3pOtunTpzN48OBQlyGE6ECf7P6Eaz65hl11u4gxxfD8Bc/z56F/lkZEHJV8h5xivF4ver1sHBNCnDxOn5M5a+fwQd4HAIxIHMGcc+aQYE4IcWWiq5CVkU5g/Pjx3HPPPdx3333ExMSQmJjI9OnTWz6fnp4OwOWXX46iKC0fA3z44YcMHToUk8lEZmYmM2bMwOfztXxeURSee+45LrvsMsLDw5k5cyYAH3/8MSNGjMBkMhEXF8fll1/ech+3201OTg4pKSmEh4czatQolixZ0vL5V155haioKBYuXEh2djYmk4kJEyZQXFzc8vkZM2awefNmFEVBURReeeWVk/7vJoQIvYL6Aq779Do+yPsABYU/DPoDL174ojQi4rh0u5URVVXxB0IziFirCb7xnohXX32VSZMmsWbNGlavXs3NN9/MmDFjuPDCC1m3bh0JCQm8/PLLXHTRRWi1WgCWL1/OjTfeyFNPPcU555xDQUEBd9xxBwCPPPJIy2NPnz6dOXPm8OSTT6LT6fj000+5/PLLefDBB3nttdfweDx89tlnLcffddddbNu2jbfffpvk5GT+97//cdFFF/Hjjz+SnZ0NgMPhYObMmbz22msYDAbuvPNOrr32WlauXMk111zDli1bWLRoEV9//TUAVqv1hP5dhBCd18L8hcxaMwunz0msKZa5Y+cyKmlUqMsSXVC3a0b8AZVPv9sbkue++Mxe6LQn1owMHDiwpYHIzs7m6aefZvHixVx44YXEx8cDEBUVRWLigd3oM2bM4P777+emm24CIDMzk8cee4z77ruvTTNy3XXXccstt7R8fO2113LttdcyY8aMltsGDRoEQFFRES+//DJFRUUkJwcvSJSTk8OiRYt4+eWXmTVrFhA83fP0008zalTwB8+rr75Kv379WLt2LSNHjiQiIgKdTtemXiFE9+DwOpi5ZiYfFXwEwJlJZzL7nNnEhcWFuDLRVXW7ZqSrGjhwYJuPk5KSqKys/Mn7bN68mZUrV7acegHw+/24XC4cDgdmsxmA4cOHt7nfpk2buP322w/7mD/++CN+v58+ffq0ud3tdhMbG9vysU6nY8SIES0f9+3bl6ioKLZv387IkSN/sm4hRNe1q24XOUtzKLQVolE03DnoTm474za0Gm2oSxMnoKLOz5KtLq4ebUajObFfpk+GbteMaDUKF5/ZK2TPfaIO3lSqKAqBQOAn79PY2MiMGTP49a9/fcjnTCZTy9/Dw8PbfC4sLOwnH1Or1bJ+/fqW00H7RURE/GQ9QojuS1VV3s97nzlr5+D2u0kIS2DO2DmMSBxx9DuLTsfuCPDxOifLtrkJqJAWq2NMP2PI6ul2zYiiKCd8qqQz0+v1+P3+NrcNHTqUnTt3kpWVdVyPNXDgQBYvXtzm1M1+Q4YMwe/3U1lZyTnnnHPEx/D5fHz//fctqyA7d+6kvr6efv36AWAwGA6pVwjRNTV5m5ixegafF34OwJiUMcw6exYxppgQVyaOl9ur8vVmF4s2OnF5g7cNSteTmahDVdUT3vf4c3W7ZqS7Sk9PZ/HixYwZMwaj0Uh0dDTTpk3jkksuoWfPnlx55ZVoNBo2b97Mli1bePzxx4/4WI888gjnn38+vXv35tprr8Xn8/HZZ58xZcoU+vTpw/XXX8+NN97IE088wZAhQ6iqqmLx4sUMHDiQiy++GAg2R3fffTdPPfUUOp2Ou+66izPPPLOlOUlPT6ewsJBNmzaRmppKZGQkRmPoum4hxInZXrOdnKU5FDUUoVW03D3kbm4ZcAsaRcKYXc2GAg9vrWiivikY8ugVr+Wq0WayE6D0q3fZuG45gx/8Bxpdx7cG8t3URTzxxBN89dVXpKWlMWTIEAAmTJjAJ598wpdffsmIESM488wz+fvf/06vXj99mmr8+PG8++67fPTRRwwePJjzzjuPtWvXtnz+5Zdf5sYbb+Svf/0rp512GhMnTmTdunX07Nmz5Riz2cyUKVO47rrrGDNmDBERESxYsKDl81dccQUXXXQR5557LvHx8bz11lsn+V9ECNGeVFXl7R1vc/1n11PUUERieCKvXPQKt55xqzQiXZTXr1LfpBIbqeG2C8OZ+usIrLsXs2bKTRQseBH77h1ULP4gJLUpqqqGJgd7HOx2O1arFZvNhsViafM5l8tFYWEhGRkZbfZJiPbzyiuvcO+991JfXx/qUuTrL0Q7aPA08MiqR/hq71cAjE8dz2NjHiPKFBXawsRxKa3xUdcUYEBPAwABVWXVDg8js/Q0bFnD7vdeoqmksM19DJYoRs3/L1rTkfcWHo+fev9uTU7TCCGEaLG1eit/XfpXShtL0Sk67h12Lzf2vzFkewnEidlW7OXJTxqwhCk8fr0ek15BoygM1OWzdd4L2PK2QvNahNsXwKgLrnZ57PWUfPJfel15+MRle5FmRAghBKqq8sb2N3hi/RP4Aj6Sw5OZP24+A+MHHv3OotPJTtYRF6mhZ7wWj1clUFXM7ndfoHrjGlCDSU2n109RrYN6h5eBqVbC9MEEZdHXH5J04ZUYrNEdVq80I+K43Xzzzdx8882hLkMIcZLY3DamrZzGN8XfAHB+z/OZMXoGVqNcObkr8AdUVm53sybPw18ujUSnVdBrFR6+2oqmqZY9bz1HxarFqM2jQlSgsLqJygZXy2PYXb6WZsTvdFK08D9k3fTXDnsN0owIIcQpbHPVZu5beh9lTWXoNXr+OvyvXNf3Ojkt0wWoqsqWIi/vrnJSXhe8lMJ3uzyc3c+Iz+Gg/OPXKPnmEwIuZ5v7KQT3jwBEmw30jDG3NCIA4fExxGS3vfBle5NmRAghTkEBNcBrW1/jHxv+gU/1kRqRSu74XE6PPT3UpYljUFTl491VDnaUBlc7wo0KlwwPY2SGQvFnb7P303fwNdqAYONRYXcTbda3NB1p0WZ6WExEGg+0AcbIcDJ++Wt6XHAtiqFjAwHSjAghxCmm3lXPgysfZFnJMgAmpE/gkbMeIdIQGeLKxNHUNvj53xona3Z5UAGdBs4faOIXQww0fP8V6x94HXfNPiB4Oqa60U1xrQOPP0Cj20CfhODX2KjTtGxa1ZmM9Dp/AimX3IzGHJpTc9KMCCHEKWTDvg3ct+w+9jn2YdAYmDJyClf1uUpOy3RyDneARRtcfP2DC2/zxa1HZhuYONKEZvc6tj3+b5pK9gDBJqTe4aGo1onTG1w5MWg1RIUZ2jymRqcldczZ9Lz8dnTRSR34ag4lzYgQQpwCAmqA/2z5D09vfBq/6ifdkk7uuFxOizkt1KWJn+Dzqyzf5uajdU4aXcF9Hn2SdVw12kyMbScFT7eN6Ta4fRTVOmhovta7VqMhJcpEosWEprnhVDQKSUOH0uuKOzAmH984kfYizYgQQnRzNc4aHljxAKvKVgFwcebFPHzmw4Trw49yTxFKAVVl9vt2iqqDSyGJURquOMtMtqGcwrfmsadVTHc/u9NLg8uLoigkWUwkR4WhazXENb5vHzKuvA1z9rAOfS1HI82IEEJ0Y+sq1jFl2RSqnFWYtCamjprK5VmXy2mZLkCjKAxM11PXFOCyEWGMSGqk5H9/5/vV37TEdD3+AF6/SrghuDE10WrC61dJijJh1B64bH9Uz1Qyr7gZy6Dx0Am/9jJgIMSKi4v53e9+R3JyMgaDgV69evHnP/+ZmpqalmPGjx/Pvffee8THWLp0Keeddx4xMTGYzWays7O56aab8Hg8HfAKhBCdkT/g57nNz3Hbl7dR5awi05rJmxe/ya+zfy2NSCdVZffzwpeN7Crzttx20ZAwZlyup+fWV1g/9WbKl3+J6vPhC6gU1TrYWFxPQWUj++e6aBWF9FhzSyMSHh/LGbfdxaDpL2EZfG6nbERAVkZCavfu3Zx11ln06dOHt956i4yMDLZu3crkyZP5/PPP+e6774iJ+ekR3du2beOiiy5qmaAbFhZGXl4e77//Pn6/v4NeiRCiM6l2VnP/svtZU7EGgF/1/hUPjHoAs94c4srET/lyk4t1+R6qG/xM/bUF1eej8sv32PvZu4fEdEvrnfgDwVM0Wo2Czx9A32olJJQx3RMhzUgI/elPf8JgMPDll18SFhYcStSzZ0+GDBlC7969efDBB3nuued+8jG+/PJLEhMTmTdvXsttvXv35qKLLmrX2oUQndPqstVMXT6VGlcNYbowHjrzIS7rfVmoyxKH4fWrON0qFnOwibh0eBh1jQEuHW5k34pFFP7vNdw1lcChMV2AML2WtBgz0WYD+9c7OkNM90R0u2ZEVVWcPufRD2wHYbqwY17+rK2t5YsvvmDmzJktjch+iYmJXH/99SxYsIBnn332Jx8nMTGR8vJyli1bxtixY0+4diFE1+YL+Hhu83O8+MOLqKhkRWXxxLgnyIzKDHVp4iCqqvJ9gYcPvnOSFKXlnkuC1/6IDFP4TdIWdj/1Ek3Fbafp1ju9FFQ1AsGYbmq0mfhIY0sTEozpnkPPy28LeUz3RHS7ZsTpczLqzVEhee4116055mXQvLw8VFWlX79+h/18v379qKuro6qq6icf56qrruKLL75g3LhxJCYmcuaZZ3L++edz4403/uS4ZiFE97GvaR9Tlk9h/b71AFyRfQX3j7wfk67zL8+favLKvby7ykHhvuBpdK9Pxe4IQOl2Ct75F7a8bS0xXW9ARd+chIkK0xMVZsASpuv0Md0T0e2aka5GVdWjH/QTtFotL7/8Mo8//jjffPMNa9asYdasWcydO5e1a9eSlNT1OmQhxLFbUbqCB5Y/QJ27DrPOzLSzpnFx5sWhLkscpKLezwerHWwsDG5ONepgwpAwzk6opOiFJ6je1HaabnGtA7vLx+C0KHQaBQU4LTGS1mvv8X37kHHV7Zizhnb8CzrJTqgZeeaZZ5g/fz4VFRUMGjSIf/7zn4wcOfKwx37wwQfMmjWL/Px8vF4v2dnZ/PWvf+WGG274WYUfSZgujDXXrWmXxz6W5z5WWVlZKIrC9u3bufzyyw/5/Pbt24mOjiY+Pv6YHi8lJYUbbriBG264gccee4w+ffrw/PPPM2PGjGOuSQjRdXgDXp7e+DT/2fIfAPrG9GX+2PmkW9NDW5hoo8EZ4OPvnSzb6sYfCIZZzuln5KJsB7WfP8UPB8V0S+qcVDa4oTkfY3d6iQkPXjl1fyPS2WO6J+K4m5EFCxYwadIknn/+eUaNGsWTTz7JhAkT2LlzJwkJCYccHxMTw4MPPkjfvn0xGAx88skn3HLLLSQkJDBhwoST8iJaUxSlS+wYj42N5cILL+TZZ5/lL3/5S5t9IxUVFbzxxhvceOONJxTBi46OJikpiaamppNZshCik6hoqmDy0slsqtoEwDWnXcPkEZMxao2hLUy08PhUFv/g4vMNLpyeYGNxRk89vxqs4l/5Otvf/ISAO7i/0RdQKat3UmF3tZmmmxZjxtxmmm4smb/6DTGjL0XRdq8TG4p6nOcJRo0axYgRI3j66acBCAQCpKWlcffdd3P//fcf02MMHTqUiy++mMcee+yYjrfb7VitVmw22yH7IFwuF4WFhWRkZGAyda3zo3l5eYwePZp+/frx+OOPt4n2ut3ulmjv+PHjSUlJYfLkyW3un5SUxMKFC9m0aROXX345vXv3xuVy8dprr5Gbm8s333zDuHHjQvTqOkZX/voLcSKWFC/hoZUPYXPbiNBHMH30dCakn/xf7MSJCagqa3Z5WLjGSW1j8LRLWpyWK0bosWz7sE1MF4KNyKbienzNMd1Io56eMWYiTa2m6VoiyPjF5V0mptvaT71/t3ZcrZXH42H9+vVMnTq15TaNRsMFF1zA6tWrj3p/VVX55ptv2LlzJ3Pnzj3icW63G7fb3fKx3W4/njK7jOzsbL7//nseeeQRrr76ampra0lMTGTixIk88sgjba4x8uabb/Lmm2+2uf9jjz3GxRdfzIoVK/jDH/5AWVkZERERnH766SxcuLDbNyJCnEq8fi9PbniS17a9BkD/2P7kjs0lzZIW4spEayXVfv6zOLgqHR2uYeJII+lV37L3n69T1Sqm2xLF1ShEmfU0uX1HiOleRMolN3WpmO6JOK5mpLq6Gr/fT48ePdrc3qNHD3bs2HHE+9lsNlJSUnC73Wi1Wp599lkuvPDCIx4/e/bsU2avQ69evXjllVd+8pglS5b85Odff/31k1eQEKLTKW0sZfLSyfxY/SMAv+33W/4y7C8YtIaj3FN0hCZXgHBT8FohPeN1nNPPSJxFYai6meL//ptdB03TLa5zkp0QQVjzKZj02HC0zZtUATR6Hamjz+6yMd0T0SEnnSIjI9m0aRONjY0sXryYSZMmkZmZyfjx4w97/NSpU5k0aVLLx3a7nbQ06f6FEKeexXsX8/Cqh2nwNBBpiOSxMY9xfs/zQ12WIHjRsreXO1izy82M31iJjQw2FxNT91Dwzr/Y0Sqme/A03dJ6J1nxEQAtg+wOxHR/jzG5dwheUegcVzMSFxeHVqtl3759bW7ft28fiYmJR7yfRqMhKyuYfx48eDDbt29n9uzZR2xGjEYjRqNsxBJCnLo8fg9PfP8Eb+4Inp4dGDeQeePmkRKREuLKxH46Deyz+XH7YONuL6PjS9j9zouHjenWOoKzwhRFIdFiIiWqbfoyvt9pwWm63SCmeyKOqxkxGAwMGzaMxYsXM3HiRCC4gXXx4sXcddddx/w4gUCgzZ4QIYQQBxTbi8lZlsO2mm0A3NT/Jv489M/otfoQV3ZqCwRUVu/0MChDT4RJg6IoXHu2GVtlLbqV/2Bdq5guwN5aB+U2F/tjuvERRlKjzRh1B0/TvQXLoHHdJqZ7Io77NM2kSZO46aabGD58OCNHjuTJJ5+kqamJW265BYAbb7yRlJQUZs+eDQT3fwwfPpzevXvjdrv57LPPeP311486c0UIIU5FX+z5gumrptPobcRqtDJzzEzGpclm9FDbUuThvVVOSmv9XFBr5Jox4fgcTXi+fo2axZ+2xHRbC55+UX8ipnsdMaMv6XYx3RNx3P8C11xzDVVVVUybNo2KigoGDx7MokWLWja1FhUVodEc6Pqampq48847KSkpISwsjL59+/Lf//6Xa6655uS9CiGE6OLcfjfz181nwc4FAAxJGMK8sfNIDD/yKXDR/kqqfby72sG24uCKh9moEGsKUPzpW4dM091nd2M2aLGGBVewEq0mIk16LN0kptuejvs6I6FwLNcZSU9PP2TgnOj+nE4ne/bskeuMiC5tj20POUtz2Fm3E4DbzriNPw3+EzqN/MYcKnWNAT5c62DVDg8qoNXA+NP1jPIsp+LT1w8zTdeJx+/HbNBxRoqVg0+4HIjp3ozGfOrMDWuX64x0RlptcNnL4/FIM3IK8niCm8L2fx8I0dV8svsTHl39KE6fkxhTDLPOnsWYlDGhLuuU5fKoLNro5KvNLjzN2z+GZeo5L2ILts//zd42MV0vxbUOHN7ggQathkRL21+KTsWY7ono8s2ITqfDbDZTVVWFXq9vc4pIdG+BQICqqirMZjM6XZf/VhanGKfPyZy1c/gg7wMAhvcYztyxc0kwHzpWQ7Q/f0Bl+TY3H61z0uAMnjDIStTxi+S9+Be/SNmurS3HNrp97G0V09VqFFKiwg4zTXdY8zTdUyumeyK6/E9wRVFISkqisLCQvXv3hroc0cE0Gg09e/Y8oRk+QoTK7vrd/HXpX8mvz0dB4feDfs/vB/5eTsuEyLZiL2+vaKK8LhjHTbBquKR3HZHf/Zvqj9a2xHT3c/sCNLi8bWK6+68VAhLTPRHd4jvfYDCQnZ3dsmQvTh0Gg0FWw0SX8mH+h8xcMxOnz0msKZY5Y+dwZtKZoS7rlGZzBCivCxBhUrjoNBdp216j6sVvqG41Tdfp8bdsTI0JN5AaZSY+0igx3ZOkWzQjEPwNWTYwCiE6K4fXwcw1M/mo4CMARiWNYs45c4gLiwtxZaeemgY/1fYAp6UEm4tRfQzY65roWbCAqlc/obL1NF2bkwqbC42iMDgtCl3zZdtTow/sUZSY7s8n/2pCCNHOdtXtImdpDoW2QjSKhj8O+iO3n3E7Wo1svO5ou8q8/P3jBiJMCo9fF4UeH6Vfvkf4Z++wrzE4lHV/TLe03tkyTdds0OILBNC1+ppJTPfkkWZECCHaiaqqfJD3AbPXzsbtd5MQlsCcsXMYkTgi1KWdsjJ66IgK1xBjhj3fLqJ20WtHjOkCmPRaeh5umu4FF5Fy8akV021P0owIIUQ7aPI28ejqR/ms8DMAxiSPYdY5s4gxxYS4slOHqqps2O1l1Q43d/4iAq1GQaeB32dto+qjlygv3dPmeKfXT0FVIxCM6aZGB/eFtJ2me05zTLd7XYzO6wug14Vu/500I0IIcZLtqN1BztIc9tr3olW03DXkLn434HdoFNls3VEKKny8u8pBQUVwE+rKHW4G6/PZ/c4L1LeK6bp9gZZNqGa9lkSLCb1WQ6LVhPbgmO6Vv8eYlNnxL6Yd1TW4yS+tp7LexYXDUzHoQnPqUJoRIYQ4SVRV5Z2d7zBv3Tw8AQ89zD2YP24+QxKGhLq0U0alzc8H3zlZXxBMVxp0cF5aPZGL/saGHw6dplvn8DIo1YqpeW5Memx4m8cLxnRvx5zVfb6GqqpSWe8kv8RGtd3VcntlnZPU+IiQ1CTNiBBCnAQNngamr5rOl3u/BGBc6jgeH/M4Uaao0BZ2imhyBfhkvZNvf3TjD4ACnJnmZFDJmzS++zU2/4GYbkmdk8oGN/un6dpc3pZmZL/uGNMNBFRKq5vIL63H7ghesE1RIDUugqwUK5ZwQ8hqk2ZECCF+pq3VW8lZmkNJYwk6Rce9w+7lxv43ysX4OoDXr/Ltjy4+Xe/C4Q42F/17+Bhj+xDXxwtpcAd/828d0w00j2SLNhtIizZjNrSappsQS+avrifmrIu7TUzX5w+wd18DBWU2nO7gxlytRiE9MZLeyVbCjKF/naGvQAghuihVVXlzx5vkfp+LL+AjOTyZ+ePmMzB+YKhL6/ZUVWVdvocPvnNS0xA89ZIapXKu/xs037yJozmmC8H1jx9Lbbh9wTfiCKOOnjHhh07T/eWv6XH+Nd0mpuv2+Nldbqewwo7XF/w3Muo1ZCZZSU+KDNn+kMORZkQIIU6AzW1j2sppfFP8DQDnpZ3Ho2MexWq0hriy7k9VVZ74sIGdZcFTL1FhKucb12NZ9zKe2koCBBuQ/etSChAfaaS60X1KxHSbnF7yy2wUVTYSCARXgcJNOrJSrKQlRKDthFetlmZECCGO0w9VPzB56WTKmsrQaXTkDM/hur7XyWmZDqIoCr2TdOyp9DLekkfK5n/jLivEQ/M0XaeX4hoHPWPNRDVfwj3ZGpwh051juvWNbvJKbJTVNLXcFhVhIDsliqRYc6f+/pRmRAghjpGqqry27TWeXP8kPtVHakQqueNyOT3u9FCX1q01OAN8vM7JyGwDWUnB5mKMpYDEiv/gXvkD7ubjGt0+imod2Jun6ZbVO1uakbbTdIfT68o7ukVMV1VVquqd5JXaqLYdSMYkRIeRnWIl1mLq1E3IftKMCCHEMah31fPQyodYWrIUgP/r9X9MHz2dSENkiCvr/j5e5+TbLW72VPn485l1FL77ItWb1raN6dY5qG0Kxnn3T9NNjgpr8zjx/fo2T9Pt+jHdgKpSVt1EXqkNe8vrhpTmZIw1hMmYEyHNiBBCHMXGyo1MXjqZfY59GDQG7htxH1efdnWX+I2zKwqoKk63SrgpuLfhl8PCKCmpZfje9/j+00WozTFdgNJ6J8V1TvbHdOMjjKRGmw+dpnvl77AMHNvlY7oHkjF2nO7gv4NWo9CrRyS9U6yYO0Ey5kR0zaqFEKIDBNQA/9nyH57e+DR+1U8vSy9yx+XSN6ZvqEvrtraXeHl3lYPocA13XxyJz9FE7cevceaSTwi4Xc0txwHB64OoR4npXoKi7TzJkRPh9vopLLezu/xAMsbQnIzJSIzEoO/ar0+aESGEOIxaVy0PrHiAlaUrAfhlxi+ZdtY0wvXhR7mnOBGlNT7eW+1kS1Fwv0d1nYdtH3xC7eJ38B00TVerUUiINAIQE25gQLKVCONhYroXXIuiN3b8izmJmlxe8ksPTcb0TrbSMyECrbbzJWNOhDQjQghxkHUV67h/2f1UOisxao08MOoBLs+6XE7LtIP6pgAfrXWyYocbVQUtAc7UbiJ928tU1lcAh07T1Wk0xEYY0CoKCrQ0It0pplvf6Ca/1EZp9YFkjDXcQHZqFMmdPBlzIqQZEUKIZv6Anxd/fJHnNj9HQA2QYc3giXFPkB2dHerSuh2XV+XLjU6+2OTC4wNUlYH6fE7PfwXNvl1A25iuwxvcH6HXakiLNrekY2B/THcsPS+/tUvHdFVVpcrmIr+knqrWyZioMLJSrMRZu0Yy5kRIMyKEEEC1s5r7l9/PmvI1AFzW+zIeHPUgZr05xJV1L/6Aysodbj5a68TmCJ52SFdKGV76OsaS9S3HOTx+9tQ0tcR0tRqFZGvYYabpdv2Y7v5kTH6pDdv+ZAyQEh/enIzp2qeajoU0I0KIU9535d9x/7L7qXHVEKYL48FRD/KrrF+FuqxuRVVVthR5eW+1k7La4GXZ46nhzMp3iSz8BuWgrakqKnaXt01MV685sCoQ368vGVfdjrn34I58GSeVzx+gqLKRglIbjoOTMckWzCZ9iCvsONKMCCFOWf6An+c2P8cLP7yAikpWVBa543LpHdU71KV1O6W1fp76tBEAs9rAmLpPiCv4GE0guPLh8QdocPmIbb4+RrhBR0ZcBFFh+oNiumlkXnlLl47p7k/GFJbb8exPxug0ZCZbyEi0dPlkzImQZkQIcUqqdFQyZdkUvt/3PQBXZF/BlJFTCNOFHeWe4li5vSpGfbBhSI3VMaqXH+2Pi+iZ/y5ab3Bjpi+gUm5zUm5zoaoQnmbFqNOiAD0iD5ye6A4x3SaXl4IyO0X7GvA3J2PMJh1ZycGZMbpukow5EdKMCCFOOStLVzJ1+VTq3HWYdWYePuthLsm8JNRldRs+v8rH65ws2epm2tUWok0BSr98l76fvouvqW1Mt7TeiS8QXB2IMOoIBA4MuIP9Md0r6HHBNV02pnvkZIyVpNjwNptxT1XSjAghThm+gI+nNz7NS1teAuC06NPIHZdLujU9tIV1M1oN5Ff4cDj9fPPeN6T++BLumirgQEy3pM6J2xfcO2LSa48wTfcXpFx8U5eM6bYkY0rrqao/kIyJjwrOjOnOyZgTIc2IEOKUUNFUwX3L7mNj5UYArjntGiaPmIxR2zV/2+5MVFVl8x4v2Um6lku4X2TdQtaSBUTVbWkZZAfB0zKF1U0EVLUlphsfaTxomu5Yev76NnRRPTr8tfxcAVWlvHlmTOtkTHJcONkpVqwR8v12ONKMCCG6vWUly3hgxQPY3DYi9BE8MvoRLkq/KNRldQuF+3y8u8pBXrmP/xtkYkJCAbvfeZH6XVuJbj7G6fUT1rwpU69RSIs2E1DVbhXT9fkDFFc2kl9mw+E6kIzp2ZyMCT+FkjEnQpoRIUS35Q14+cf6f/DqtlcB6B/bn9yxuaRZ0kJcWddXZffzv++crMsP/vZvcNbg+GQhG8q/ADW4OTM4TddJbZObfokWrGHBN+Qkq6nNY3XlmK6n1cyY1smYjCQLGUkWjKdgMuZESDMihOiWShtLuW/pffxQ/QMA1/e7nknDJmHQdq3R6p1NkyvAZ+tdfPOjC18ANK4GBtZ8Qe+SjzHhAIIx3dI6J/sa3Oyfptvo9rU0I/t15Ziuw+Ul/+BkjFFH75TgzJhTORlzIqQZEUJ0O4uLFvPwyodp8DQQaYjksdGPcX6v80NdVpfm86ss2eLm4++dONwqeF1kVy+jX8l7RAbqgse0iukGmldHosIM9IzpPtN0ba2SMfsv02YNN5CdYiUpTpIxJ0qaESFEt+Hxe/jb+r/xxvY3ADgj7gzmj5tPSkRKiCvrulRVZX2Blw++c1BlD4DfR3L1Os4oeZtYX1mbGO6OCjuNzVcSjTDq6BljxtJqr0RXjemqqkq1zUVeqY2qemfL7fFWE1mpUcRLMuZnk2ZECNEtFNuLyVmWw7aabQDc1P8m/jz0z+i1snHwRO3e52PBiiZ27/NDIEB0zRYGlSwgyZOHBhWV4EmY/W/DidYwSuoc9Iw2Ex1+mJjuJTejCYsMzYs5AaqqUlbTRF7JgWQMQEpccGZMlCRjThppRoQQXd4Xe75g+qrpNHobsRqtPD7mccanjQ91WV1eeZ2f3RU+zLUFDCx7n15NG9ApgZZpukW1DhIijSRaghtSY8MNxLZqQjR6HaljxtLz8q4V0/U3z4w5JBmTEEHvFKskY9qBNCNCiC7L7Xczf918FuxcAMDg+MHMHzefxPCuO0Y+lBqcASptAXonBt8aTld2MTr/A5KqV2PSeEEJbkQtqnW0TNOtsKn0sJhQOLBComgUkoaNoNcVd2BMygjNizkBHq+fwormZIw3mIzR6zRkSjKm3UkzIoTocKqq4nX78Xn96PRa9EbtcZ9z32vfS87SHHbU7gDg1gG38qchf0Kvkd9aT0ThPh9//7gBox4eGG+j7H8vUrVpDRmqCpq2MV0ABYVEa3CabuuvXDCmewfm3oNC80JOgKN5ZszeVsmYMKOOrGQLPXtESjKmA0gzIoToMF63j+JdNeRvrKCyxEbAr6LRKiSkWskakkhan1j0xqP/WPps92fMWD0Dh89BtDGaWefM4uyUszvgFXRfqbFazP4GwrZ9w+pvPsDMgTkqFXYXe2sczbtEIC7CSFq0uc00XWvPNHpf+TssA8/pMjFdW5OH/NJ6SqvaJmOyUqwkSzKmQ0kzIoToEJXFNlZ+uIPKYjuKRiHCYkRn0hLwBdizvYrCrZUkpFkY86u+JKRZD/sYLp+LOWvn8H7e+wAM6zGMuefMpUd419mP0FnsKPWyfJub350fjupsovjj1xi/9kuMvgYOfguONOlQUY8Q040j81fXdZmY7v5kTH6pjcpWyZg4q4lsScaEjDQjQoh2V1lsY/FbW7DXOEhIs6A76Nx7ZEwYPq+fymI737y9hfOuHXBIQ7K7fjc5y3LIq8tDQeGOgXfwh0F/QKeRH2PHo7zWz/vfOdi8xwt+H4l7FhG96Q18jXZMBGerVDS48flV0qLDAAg36BiUGtVySXdojulefCU9zr+6S8R0VVWlvMZBXmk99Y2SjOls5H+xEKJded0+Vn64A3uNg6SMqCP+1qnTa0nKiKK8sJ6VH+7gktuHtZyy+ajgIx7/7nGcPiexplhmnzObs5LP6siX0eXZHQE+Wudk+TY3AX8AQ/lW+pcuxODfhk/jQwVqGt0UN0/TVVCIjzRiaj4Vs78R6Wox3f3JmIIyG03NyRhNczImK9lKeJjsMeoMpBkRQrSr4l01VBbbSUizHHX5W1EUElItVBbbKcmrpcdpEcxaM4sPCz4EYFTiKOaMnUNcWFxHlN4tuL0qX212sWijE7dHRVe1m+ziT+jr/p4IrQtVcyCm6/AE36z1Wg2p0WYMrfaEdLWYbjAZ00BhuQ23JGM6PWlGhBDtRlVV8jdWoGiUQ07NHInOEEzWLFu/lg/yX2S3bTcaRcMfB/2R28+4Ha1G3kSORSCgsnqnh4VrHdQ3qWjrSum59wv6N64gVt8AWnD5/OyuamqJ6WoVheSosEOn6XahmK7D7aOg1HZQMkZL72QrvSQZ02lJMyKEaDdet5/KEhsRlmM/H6+iUpD0PV+q7+CzeYkPi2fu2LmMSBzRjpV2L1uLvLy32kFJjR9NYw0Jhd/Qv+5rko21KK3OSmgVhUa3DwWFHhYTKdFh6DUHVq+6UkzX3uQhr9RGaXXj/qHBWJpnxiTHhqPRyKbUzkyaESFEu/F5/QT8KjrTsa1muHHyof41Npm+A2BUwpnMHT+H2LDY9iyz21BVlWcXNbKp0IviasBasJx+NV+QbihHa1Tx+APUNnlarpiq12rISogg3KA7NKZ71a1Yzji7U8d0VVWlxh6cGVNZd1AyJsVKfFSYJGO6CGlGhBDtRqfXotEqBHyBox5brhTxlv45qjUVKKqGsfZLyb16GqYwQwdU2j0oikIPk5uwncs5bd/nZOv3YjD68KsqxXUuym1OAqpKmF6LtXnjZoz5wL9veEIcmROvJ+bMizt1TPdAMsZGfaO75fbk2GAyJjpSkjFdjTQjQoh2ozdqSUi1smd7FZExYYc9RkVlrXYJn+rexKf4sKrRXLD3t4zOGIVRZoD8JJdH5YtNTgb01JMRC6VfLCDu0w+Y0OAi3OgmoAYvWFZS58QXCDaE4UYd2oNOWXSVmK4/EKC4spH80kOTMb2TrURIMqbLkmZECNFuFEUha0gihVsrWy793poLB//Tv8KP2nUAnOYfxETHzTQ5VLKHJskS+1EsXOtg8SYnW5f8yLnlz+OurQbArIXqVjFdAJNOS1qMmZguOE3X4/Ozp7yB3eV23N7g69HrNGQkWshMsmA0dN5VHHFspBkRQrSrtD6xJKQF47qtrzNSquzhLf1z1Goq0ahaJviuZLTvQvY1x4BTs2NCXHnno6oqHh8Y9QqqqjJKs5Gda5eQ7tmEy1jX0mSoKi2NyP6Ybnykkf27QoIx3XHNMd2EUL2co3K6fRSU2dhT0SoZY9DSO0WSMd2NNCNCiHalN+oY86u+fPP2FsoL64lPjWRd2BI+1y3Ar/iIUmO51vNHkt3p7CuxY40zc/bEfsc0o+ZUsrfKx7urHJj0Cjdm72H3Oy9Qv2sbYxVQjNDk8RFm0KEBNAr0ijHj9PoPiekmDh1O+lV/wJiYHtLX81PsTR7yS22UtE7GmPVkpUSREifJmO5I/rcLIdpdQpqV864dwNcfrudV9Sl2638EoI97EL+svQF/vZZKNbgicvbEfsSnWkJccedR0+Bn4Ron3+3yoGmswZy3lBWuDzBrghs33V4/Rc3TdDNiw+nRnJSJCW+78Te+Xz8yrrq908Z09ydj8ktt7DsoGZOVYiVBkjHdmjQjQogOsc9UxKuJcylrKkOLlvH1lzPEPg6tVkNCfyvZQ5NIzY6RFZFmDneAzze4+PoHF/6mBkx5K8moWkLfsGLMGg8ef4DSOieVDe6Wabou76Gppc4e01VVlfJaB/klNuraJGPMZKVESTLmFCH/64UQ7UpVVV7b9hpPrn8Sn+ojNSKV+WPn0yeyb8umVr1RK7/1NvP5VZZudfPJ904aG5wYC9aQVrqE/qbdRIU7DonpAkSFGUiLMRN+8DTdThzTPZCMsdPUfAVYjdKcjEmRZMypRpoRIUS7qXfV89DKh1hashSAC3tdyIzRM4g0BJMbBpP8CNpPVVU2Fnp5f7WDyloPhr0bSdrzDf31+fSIqG/ZnFpQ1URtU3AFIdygo2esGWurCLTREtkc072qU8Z0vb7gzJjdZa2SMVoNGUmRZCRZMUky5pQkPwmEEO1iU+UmJi+bTEVTBQaNgftG3MfVp10tKyCHsXufj3dXOsgv86Av20Zs/rf0VbbTM7waBRVVVVv+3ZKtJhxu3xFiur8k5ZKbOmVM93DJGJNBS1aylZ49ItHrJBlzKpNmRAhxUgXUAC9veZl/bvwnftVPL0svcsfl0jemb6hL65Qq6vzMfs+Grmo31p3fkuX9kSzzPnSKv2WarsWkJz3WDECEUcegtKiWJqSzx3TtjuZkTNWBZEykWU92ipWUuAhJxghAmhEhxElU66rlgRUPsLJ0JQC/yPgFj5z1COH68BBX1rn4/Co6bfBN2Fy9jf5bP0GpyOe08DLCDB4a3T6Kax3YmvdSeHwB0mLCDkR06dwxXVVVqbW7ySutb5OMibUEZ8YkREsyRrQlzYgQ4qT4vuJ7piybQqWzEqPWyP0j7+eK7CvkTaeVQEDl6x9cfLXZxb2j66n/9EWqNq1lsApKpIrL6yev2klN856Q1tN0tUrrabr9yLj6DsyZA0P1Ug5LVVUqaoMzY+oaDiRjkmLNZKVYiYk0hbA60ZlJMyKE+Fn8AT///vHfPLv5WQJqgAxrBrnjcukT3SfUpXU6igIbtlTj+m45H3/2JX1NxcHbgZomD/mVjS0x3bgII6nRZkxdYJquP6BSUtlIfpmNRueBZExaQgRZkowRx0CaESHECat2VjN1+VS+K/8OgMt6X8aDox7ErDeHuLLOI6/cS0qMFoPPQdFHr9L/6yUkOzWkmWraHBdp0qEoYDUdLqYb3xzT/WWniul6fQH2VNgpaJWM0TUnYzKTLJgM8hYjjo18pwghTsia8jVMWTaFGlcNYbowHhj1ABOzJoa6rE6jot7P+6sdbMp3co66gswdr+FrasACRBhhn91Fo9tHVnwEAAathkGpURhbrYQciOlejaI3HOGZOp7T7WN3uZ09FXZ8/gPJmN7JwZkxkowRx0uaESHEcfEH/Dz/w/P8a/O/UFHJisoid1wuvaN6h7q0TqHBGeDjdU6WbXGiKdlGxK7l2NU8fBENqEDNQdN0EyJNWJqvt7K/EdGFGel1/i+bp+lGhOqlHMLu8FBQaqO4dTImTE9WqpVUScaIn+GEmpFnnnmG+fPnU1FRwaBBg/jnP//JyJEjD3vsiy++yGuvvcaWLVsAGDZsGLNmzTri8UKIzqvSUcn9y+9nXcU6AH6d/WvuH3k/YbqwEFcWeh6fytebXXy+3oGvbDemHUtJcufRP7yESK2rJabr8PgA0Gs0pEaHEdHqwm8avY7Us8fTc+KtnSqmW2N3kV9io6LO0XJbrCU4M6aHJGPESXDczciCBQuYNGkSzz//PKNGjeLJJ59kwoQJ7Ny5k4SEQ//zLFmyhN/85jeMHj0ak8nE3Llz+b//+z+2bt1KSkrKSXkRQoj2t6p0FVNXTKXWVUuYLoxpZ03jksxLQl1WyAVUlTW7PCxc48RWXIxxx1ISbTvoH1FCvKUBjz/A9vLGlpiuVlFIjgo7dJrusBGkX/n7ThPT3Z+MyS+1Uds6GRNjJitVkjHi5FJUdf9i27EZNWoUI0aM4OmnnwYgEAiQlpbG3Xffzf3333/U+/v9fqKjo3n66ae58cYbj+k57XY7VqsVm82GxSLTPIXoSL6Aj2c2PcO/f/w3AKdFn0buuFzSremhLawT2F7i5b1VDkr2VGLcuRxL5Q/0DS8jxVjTclGygKqyucSGxxcIxnSjTOi1B/ZUxPfvH5ym20liuv6ASklVI/mlrZMxkJYQSe9kC5HmzrN3RXR+x/r+fVwrIx6Ph/Xr1zN16tSW2zQaDRdccAGrV68+psdwOBx4vV5iYmKO56mFECFQ0VTBlGVT2FC5AYCr+1zN5BGTMelO7d+Ky2r9vLfawZadtRjzVhFVspEsUymZ0ZUEAj5K6tykRIehIRhxzYqPQK/TdOqY7hGTMYmRZCZLMka0r+P67qqursbv99OjR482t/fo0YMdO3Yc02NMmTKF5ORkLrjggiMe43a7cbsPLAva7fbjKVMIcRIsK1nGgysepN5dT7g+nOmjp3NR+kWhLqtTyNtjY9cXS7HsWUcvfRl9osrRKV7K6g9M0zXoFHo0n8qIbLUvpLPFdCUZIzqDDm1158yZw9tvv82SJUswmY78m9Xs2bOZMWNGB1YmhNjPG/Dy1IaneGXrKwD0i+lH7rhcelp6hrawEHJ7VSptflIsAUq/fAfNZ+/Tt9JKuqUSs8ZNZYOLkjonvkAACE7TDdO3/fFqtESSccmV9Divc8R0G1rNjGmeW0dEWHBmTGq8JGNExzquZiQuLg6tVsu+ffva3L5v3z4SExN/8r65ubnMmTOHr7/+moEDf/rc6NSpU5k0aVLLx3a7nbS0tOMpVQhxAsoay5i8bDI/VP0AwHV9r+Ovw/+KQRv6N89QKan28eRHNnTFP3JJ/b/w1VYCcHq4nZomD3m1jpaYrkmnPXSabpiRXhdcTMrFN3WKmG6t3UVeqY2K2gPJmBiLkeyUKEnGiJA5rmbEYDAwbNgwFi9ezMSJE4HgBtbFixdz1113HfF+8+bNY+bMmXzxxRcMHz78qM9jNBoxGo3HU5oQ4mf6pugbHlr5EA2eBiL1kTw65lEu6HXk06mnAlVV0e1ejfGbFWhsFdRZ7US2OrNSaXfh9vlbYrrxFhP7T2q0xHQvvw2dNT4k9e+nqir76pzkldS3ScYkxpjJTrESYzm19wCJ0Dvu0zSTJk3ipptuYvjw4YwcOZInn3ySpqYmbrnlFgBuvPFGUlJSmD17NgBz585l2rRpvPnmm6Snp1NRUQFAREQEERGh/y1BiFOd1+/lb+v/xn+3/xeAM+LOYN7YeaRGpoa4stAorvaxdKubS3rsZs97/6I+bztnYsIc7cbl9eJVtOg1CgrQM8ZMvdNLUieN6QYCKsWHScakxgdnxkgyRnQWx92MXHPNNVRVVTFt2jQqKioYPHgwixYtatnUWlRUhEZzYMPTc889h8fj4corr2zzOI888gjTp0//edULIX6W4oZiJi+dzNaarQDc2P9G7h16L3rtqTfYrLYxwIdrHaz5vhzDzuX4G7+kl6kKAL3axO6q4DTdREsY6bHB2TsRRh0RxgM/RuP69SPz6t9jzjwjJK9hP68vwN59DRSU2XB59idjFNITLWQmWQgzSjJGdC7HfZ2RUJDrjAhx8n2550seWfUIjd5GLAYLM8+eyfi08aEuq8O5PCqfb3Sy+LsqlB2r0Bf/QIqhmr7mUgw4Kal3Uml3t0zTjY8wkhkfQeudFdZePYMx3QFjQhrTdXl8FJS1TcYY9Vp6J1tIT7RIMkZ0uHa5zogQoutz+93MXzefBTsXADA4fjDzxs4jKSIpxJV1LH9AZfk2Nx+tqsOz5TsMe74nVqmjv7UEi7aR8noXZc0xXQBrmIGeMWGEGw6K6V7+W2JG/SKkMd0Gh4eCMhvFlW2TMVnNyRitJGNEJyfNiBCnkL32veQszWFHbfC6QL8b8DvuGnIXes2pc1pGVVU27/Hy/go7tZvXY8xfRUygnv7hJfQw1KMAe2udlNucQDCm2zPGjDXswL9RZ4npHjYZE2kkK8VKYoxZkjGiy5BmRIhTxGe7P2PG6hk4fA6ijdHMPHsm56SeE+qyOtSeSh/vrGhkz/c/YNy1Aourmj7mcnqaqlDVAErzyZckqwmb00tKVFini+nuT8bkl9ZTY2+bjMlKsRIryRjRBUkzIkQ35/K5mLN2Du/nvQ/A0IShzBs7jx7hPY5yz+5DVVVe+aaRtct3YNyxlPCGfWSG7SMrpgKHy8W2Mgd6rYa+iZEAGLQazkixtjQhnSGmGwiolFQHkzENjmAyRlEgTZIxohuQZkSIbmy3bTc5S3PIq8tDQeH2gbfzx0F/RKc5tf7rN+RtwffRB5jzm0g11tA3poyAz0H+Pgc254Fpuh5/AEPzEDuFzhHTlWSMOBXId7EQ3dRHBR/x+HeP4/Q5iTHFMPuc2YxOHh3qsjqEz6+yZIubFMoIfPsiVZvWkR7QEB9lxEQDxdXBmC6AgkIPi5GUqLA203Tj+vUn8+o7QhbTdXmCM2MKyxvw+YOXmZdkjOiupBkRoptxeB3MWjOLDws+BGBk4kjmnDOHeHNorwLakd5fXMaK/y0jsWIloy3bUQC94kfx2dhc3tAS040LN5IaY247TTfEMd1Gp5f8UhvFlQ2SjBGnDGlGuhFVVfG6/fi8fnR6LXqjVnbTn2Ly6/LJWZpDga0AjaLhD4P+wB1n3IFWE/rpsO3N51fB1UTRh68QvfgroirTSTVXtTkm0qjDoNNg0ms7XUy3tsFFfqmN8poDyZjoSCPZkowRpwBpRroBr9tH8a4a8jdWUFliI+BX0WgVElKtZA1JJK1PLHo5r9ytqarKwvyFzFozC5ffRXxYPHPHzmVE4ohQl9buKm1+3l9hx/PjcoYV/QdfUwMGYHz0FqoaXGyr89A/yYICaBSFASlW9K1WF4Ix3avocd5VHR7TPZCMsVFjd7XcnhhtJitVkjHi1CHvUF1cZbGNlR/uoLLYjqJRiLAY0Zm0BHwB9myvonBrJQlpFsb8qi8JadZQlyvagcPr4NHvHuXT3Z8CMDp5NLPOnkVsWGyIK2tfja4AH69pYuVXG9DvWInGWU9WtAezFmqaPBS3mqZb3egmPiI4fHN/IxKM6V5CysU3dnhM90jJmP0zYyySjBGnGGlGurDKYhuL39qCvcZBQpoFnb7t0nJkTBg+r5/KYjvfvL2F864dIA1JN7Ozdic5S3PYY9+DVtFy15C7+N2A36FRuu/mRq9PZfEPThZ9tgW2LMPYUEW83kb/qBJ8nka21Dpo8vgA0Gs0pESHERtxYAp4KGO6Xl+Aon0N5B+cjOlhITNZkjHi1CXf+V2U1+1j5Yc7sNc4SMqIOuL5ZJ1eS1JGFOWF9az8cAeX3D5MTtl0A6qq8u6ud5m7di6egIcEcwLzx85naI+hoS6t3QRUlXV5Hv736S6cG5aiqy3GonXQ31JCjN7Grn0NLTFdjaKQbA0jKeow03Sv+gPGHr06tHaXx8/ucht7yhvwtkrGZCZbyEiMRK/r/nt6hPgp8q7URRXvqqGy2E5CmuWoG9sURSEh1UJlsZ2SvFoyBiR0UJWiPTR6Gpm+ejpf7PkCgHNSzmHm2TOJNkWHuLL2s7PUyzuf76Fq9XL0FTuJ0HjoG1FKqqkWBRUVBVXtfDHdRqeXglIbRZWNLTNuwk16slOspCaEo9V03xUsIY6HNCNdkKqq5G+sQNEoh5yaORKdIZisydtQTvrp8bIzv4vaVrONnKU5FDcUo1N0/Hnon7nx9Bu77WmZ8jo/731dzq5vlqMv/gETXrLMFfQ0lLPP7sBnCEOvCV7EPSMuHEVROkVMt67BTV5pfdtkTISRrFQrSZKMEeIQ0ox0QV63n8oSGxEW49EPbiXCaqSyxIbX7cdgki99V6KqKm/teIvc73PxBrwkhScxb+w8BicMDnVp7aa83M6sOYvRF36P0e+lp6mKrLAyau0N/FB5YJpurxgzAGGtGvNgTPcGYs78BUoHrT6oqkplvZP8EhvVrZIxPaLDyE6JIsZilCZEiCOQd6QuyOf1E/Cr6EzHd55Zo9Pgc/nweaUZ6UrsHjuPrHyEr4u+BmB82ngeH/M4VmP324wcUFXwein9YgF7P3+ftPIE0Cr0jSzB4ahne4kTbyC45yLcoCMqrO204VDEdAMBldLmZIz94GRMshVLuCRjhDgaeUfqgnR6LRqtQsAXOK77BXwBNNpjP7UjQu/Hqh+ZvGwypY2l6DQ6/jrsr1zf7/pu9xt2QFVZtc3JooVrOG/fv9HYKgAYFtlIvcPN7nIHruaYrlGnpWeM+TDTdDs2puvzN8+MKbXhbE7GaDUK6YmR9E62SjJGiOMg/1u6IL1RS0KqlT3bq4iMCTvm+zXa3KT3j0dvlGaks1NVlde3vc7fN/wdX8BHSkQKueNyGRA3INSlnXSqqlK7YQVLn/oGZ7WD7WE6Tg8Pfk6rqNicXlw+f0tMNyHSxP5rlmn0OlLPOZeeE2/tsJiu2+MPzoypsOP1tU3GpCdGYpBkjBDHTZqRLkhRFLKGJFK4tbLl0u9H4/P4UVWV7KFJ3e636u7G5rbx0IqHWFKyBIALe13I9NHTsRgsoS3sJCut8aEr30b5wn9Rn7eD/r5wasIjSNCU4/SqLXtAUqPDMOg0JFkPiukOHxmcpttBMd1Gp5eCsuZkTGB/MkZHVoqVtIQIScYI8TNIM9JFpfWJJSEtGNf9qeuMQPPGupJgDDg1O6YDqxTHa1PlJiYvm0xFUwV6jZ77RtzHNadd060ayPqmAP9btJtNny3lNNsy+oeXABBGA2rDPrY1uYk2GzitRyQABq2G1KgDK4Bx/fuTefXvMWd0zCpRXYOb/NJ6ylolY6IijGRLMkaIk0aakS5Kb9Qx5ld9+ebtLZQX1pOQakFnOHSFxOfxU1lixxpn5uyJ/eSCZ51UQA3wytZXeGrDU/hVPz0je5I7Lpd+sf1CXdpJ4/KqfL6sjGUfLEHZuwWdGsBl1OPxByird7LP7m6ZpqtRFAIqtB5Q25Ex3ZZkTKmNalvbZExWSnBmjDQhQpw88s7UhSWkWTnv2gEHZtMoChFWIxqdhoAvQKPNjaqqJKRZOHtiP+JTu9cyf3dR66rlwRUPsqJ0BQC/SP8F086aRoShY+eltBd/QGXZxloWvfUt3l0b0fh9ROsa6WsuxuWoZnOxE39zTNcapqdnjPkw03Q7JqYbTMY0NSdjPEBzMiaueWaMJGOEaBfSjHRxCWlWLrl9GCV5teRtKKeyxIbP5UOjVUjvH0/20CRSs2NkRaST+r7ie6Ysm0KlsxKj1sj9I+/niuwrusVv3aqq8kO+g/dfX0LDj+vQeJyEa1z0iywlyVjHPruLkrrgqY9wg46eMWasraK6RkskGZdeTY9zr2z3mG5LMqbMhtPdNhmTmWzFLP9/hGhX8j+sG9AbdWQMSCD99Hi8bn/Lpla9Udst3tS6o4Aa4N8//ptnNj1DQA2Qbkknd1wup8WcFurSToq9FW4W/Hcl5d99h8Zpw6j46BNeRrJ+HyZd8HsyIdJIbZOHBIuJ2DYxXRO9LryYlItvQmMKb9c6D5+M0ZCRZCUjSZIxQnQUaUa6EUVRMJh0ckGzTq7aWc0Dyx9gdflqAC7NvJSHznwIs94c4sp+vhq7j3cXrGPn4hVoGqrRESAjrJIeSjFltXZ2qCqDUqNQCO4L6Z904NThgZjubeisce1aZ5PTS/6RkjHxEWi1kowRoiPJu5YQHWhN+RruX34/1c5qTFoTD4x6gIlZE7vFCpZ91w8sfnEBedt8aIBkYw299EVU1deT5wzuv9AoCg6Pr82ekI6M6dY3uskrsVFW09RyW1SEgeyUKJJiJRkjRKhIMyJEB/AH/Pzrh3/x/ObnUVHpbe3NE+OfoHdU71CX9rP4/CrFOwtxfPECVZu/J16FNGNPkvSVNNiryK92Az8xTbcDYrqqqlJV7yTvoGRMQnQY2ZKMEaJTkGZEiHZW5ahiyvIprKtYB8DlWZczddRUwnTHfvXczmhvYQX//scXBAp/4DzrD2gVFQU4zbSbH0psLTHd2HAjaTFhmFrtvwjGdG/DMmB0u8V0A6pKWXUTeaU27E3NyRggJT6crJQorJKMEaLTkGZEiHa0qnQVU1dMpdZVS5gujIfPfJhLe18a6rJ+Fl9TA0UfvsLebxYR2NeHgKrQ4DcSpXMCYNJriTDp0CgcGtPtkRCM6Y66qN1iuj5/gKJ9DeSX2XG6fUAwGdOrRyS9ky2YTfqjPIIQoqNJMyJEO/AFfDy76Vn+/eO/UVHpE92H3HG5ZFgzQl3aCdtX7eSLNz4nc+dr+JoaUIDhlnwcDht7KpoYkGxFp1FQgL6JkS2XbodWMd3zrkLRtU8z4Pb6KSy3s7v8QDLGoNeQmWQhI9GCQQZECtFpSTMixElW0VTBlGVT2FC5AYCr+lzFfSPuw6QzhbiyE9PY5OWD/y5h89drwNEAEQbSTFDb5KG4tr5lmm6F3dVy2fb9jUhHxHSbXF4KSoPJGH9zMsZs0pGVbKVngiRjhOgKpBkR4iRaVrKMB1c8SL27nnB9ONPPms5FGReFuqwT4vEF+PyD1axcuAK/rRaAOL0dnd/OllIbTZ7gKZDW03T364iYbn2jm/xSG6XVB5Ix1nAD2alRJEsyRoguRZoRIU4Cb8DLPzf8k5e3vgxAv5h+5I7LpaelZ4grO36qqrJ88Y98+eZinPsqAIjUOugXXkJdbSlFrWK6ydYwkqIOnqY7ivQr72iXmK6qqlTZXOSX1FPVOhkTFZwZE2eVZIwQXZE0I0L8TGWNZUxeNpkfqn4A4Dd9f0PO8BwM2q6X1vhxQwEfvryIut1FABg1Hk4zl9HTVIOCSpNOE5KY7v5kTH6pDdshyRgr1nDjSX9OIUTHkWZEiJ/hm6JveHjlw9g9diL1kcwYM4MLe10Y6rKOW1FhBe+/8AnFPxagqAG0+MkwlWN07yFeo0chuPkzNTqM5CjTQTHdXs3TdE9+TNfnD1BU2UhBqQ2HJGOE6LakGRHiBHj9Xv62/m/8d/t/ARgQO4D54+aTGpka4sqOj6+pgdef+h9bV+9A9ftRUEkzVmH17aa62o5NVfH6DJzWIxIAQ6uVkPaM6e5PxhSW2/HsT8boNGQmSzJGiO5ImhEhjlNJQwmTl05mS80WAG7ofwN/GfoX9Nqu81t6wOOhdNFbFH3xPxz7rKj+ROL19fRgN7W1NewLBBsAs0FHD0vbFJDRaiHj4qvocf7Jj+k2ubwUlNkp2tdwIBljbJ4ZkxCBTpIxQnRL0owIcRy+2vsVj6x8hAZvAxaDhcfHPM65Pc8NdVnHzO/zsfidr/CtXIC5sRiAbLODsEAN9vpy9jXHdI06LWkx5sNM072ElItvPOkx3SMmY1KsJMWFo5FNqUJ0a9KMCHEM3H43uetyeXvn2wAMih/EvLHzSI5IDnFlx0ZVVWo3LGfhi5+wrVghRm9mtDW4CdSg+DEF6qj0+dFpNKQ2x3Q1ze//wZjuefSceOtJjem2JGNK66mqP5CMiY8KzoyRZIwQpw5pRoQ4iiJ7ETlLc9heux2AWwbcwt1D7kav6RqnZep3bGbPey9Qn7eD5ICBAqUfUZp67C4fVlPwR0CixQio9LC0f0w3oKqUN8+MaZ2MSY4LJmOiIiQZI8SpRpoRIX7C54WfM2P1DJq8TUQZo5h59kzGpo4NdVnHpCKvgP+98CGOgu0Mi9wNgBJwkR5YTW21kz16LQNTo1A4cM2Q/eL6n07m1Xec1Jiuzx+guLKR/DIbDteBZEzP5mRMuCRjhDhlSTMixGG4fC7mrpvLe7veA2BowlDmjp1LYnhiiCs7Ovu+fXz0r/f4Yd1e/AEFiCHTVIrNXsc+u6tlmq7ZoMMfUNFpDpwKaYnpnjHmpNXjaTUzpnUyJiPJQkaSBaMkY4Q45UkzIsRBdtt2k7M0h7y6PBQUbjvjNu4cfCc6Tef+7+JpsLPopXf5bslO3F4FULDqGklQC8gv3YdfDTYhVpOetBgzEcZW03QT4sn89U3EjJpw0mK6juZkzN6DkjG9U4IzYyQZI4TYr3P/dBWig31c8DGPffcYTp+TGFMMs8+Zzejk0aEu6ycFPB6Wv/k+3366iQYHgEKYxk2/8FLM/n3s3GcHgishPWPMRIUdOB1itFrIuORqepx35UmL6dpaJWPU5tus4QayUqwkSzJGCHEY0owIATi8Dmavnc3C/IUAjEwcyZxz5hBvjg9tYT9B9fv54dPP+fydFVTVBd/29YqPdGMp2eHVaBUVFT1xEUaiwvTERhjbLaarqirVNhd5pTaq6p0tt8dbTWSlRhEvyRghxE+QZkSc8vLr8slZmkOBrQAFhT8O+iN3DLwDraZz7mVQVZU9K5fzyauL2FPmAxQUAiTrK9A25tPQ6EY1R4GioABZ8REt9z3ZMV1VVSmraSKv5EAyBiBFkjFCiOMgzYg4ZamqysL8hcxaMwuX30VcWBxzz5nLyKSRoS7tiOw7N/PDG//ho/V6/GgBhQRdDWGuXTgagqdjNIpCo9vX5nTMyY7p+ptnxhySjEmIoHeKVZIxQojjIs2IOCU5vA4e++4xPtn9CQBnJZ3FrHNmERd28i7qdTI5SnZT+M6/qPphPagqcYYs3H4tUb48muqqcBBcH0lonqbbeoZM3OkDgjHd9NN/dh0er5/CiuZkjDeYjNHrNGRKMkYI8TNIMyJOOTtrd5KzNIc99j1oFA13Db6LW8+4FY3S+dIdruoKlr78JmtX7WFExC5MmuDekDPCC/ihuIam5oRMbLiRtOgwTPqDp+nehuWMn78B93DJmDCjjqxkCz17REoyRgjxs0gzIk4ZqqryXt57zFkzB0/AQ4I5gXlj5zGsx7BQl3YIX1MDRR++TMm3i1hbmUG9L4I8RyJnRATnyYRpVeIjjLi8fnoeEtNNIPPXN56UmK6tyUN+aT2lVQeSMZbmmTGSjBFCnCzSjIhTQqOnkUdXP8rnez4H4OyUs5l19iyiTdEhrqytgMfDtg8XUP3NB+AI7gHpH15MYWME7rp8HEYz5ubVj/RYM0rzJlU4eTHd/cmY/FIbla2SMXFWE9kpVuKjwiQZI4Q4qaQZEd3e9prt5CzNoaihCK2i5c9D/8xNp9/UqU7LqH4/e7/5lC/e/Ib8Sj2ZYZH0DbdT2+ShuLYeX/M03Qqbhsy4YBR3/6qELsxEr/+7lJRf3ojGZD7xGlSV8hoHeaX11DceSMYkx4aTnSrJGCFE+5FmRHRbqqry9s63mb9uPt6Al6TwJOaNncfghMGhLq2FqqpUrVvGl69+yNa94FGDb/g1HhNbbDaa3MGkSutpuvtp9PpgTPfyW9FZYk+4hv3JmIIyG03NyRhNczImK9lKeJgkY4QQ7UuaEdEt2T12pq+azld7vwJgfNp4Hh/zOFajNcSVHWDbsYlv//M2G7Y7aQoEm4wIrZMobx6uhnKaCK5+JFnDSI46zDTdq36PMaHnCT9/MBnTQGG5DXerZExGooXMJAtGgyRjhBAdQ5oR0e1sqd5CztIcShtL0Wl0TBo2id/2+22n2efgKC7gu1dfY9W6aup8EYAJo+LltPAy0kzVVNgcFDe1X0zX4fZRUGo7KBmjpXeylV6SjBFChIA0I6LbUFWV/27/L39b/zd8AR8pESnkjstlQNyAUJcGgLtmHz+8+TLLlhRS5ooCItDiJ5YSepsqiDMFVyISLSZizIaTHtO1N3nIK7VRWt1IcyIYizk4MyYlLhyNpnM0a0KIU480I6JbsLltPLTyIZYULwHggp4XMGPMDCwGS0jrAvA12tn1/qt8+/kPFDZEoRIFqMQq+9A07kQNuCh1aYkNj0IheGpmfyMSnpBA5hU3EjPyxGK6qqpSYw8mY/bVSTJGCNE5STMiurxNlZu4b9l9lDeVo9fomTxiMteedm3I32QDHg+li96k6Iv/kV9tYndjOgAWpQ6jYwf42k7Tbc1otZBx6dX0OPfEYroHkjE26hvdLbcnxwZnxkRHSjJGCNF5SDMiuqyAGuDVra/y1Ian8Kk+0iLTyB2XS//Y/iGtS/X7KV/2Cds/eAetvQKANFMTZc4IFEcxeCoBMOq0pEWHHTpN92fEdP2BAMWVjeSXtkrGKAo9e0TQO9lKhCRjhBCdkDQjokuqc9Xx4IoHWV66HICL0i/ikbMeIcIQcZR7th9VValdv5Qf3v4vy/NMeALxnBtdiVYJoEHlNOMudtU3oNNoSIkOo0ekseVaIT83puvx+dlT3sDucjtub/CaJPuTMRlJFkySjBFCdGLSjIguZ/2+9dy37D4qHZUYNAbuH3U/V2ZfGdLTMvadm9j9zgvU5+/Er2pwBwbgDWgpbtSTHhk8TRJtNpARF05suBGd5uTEdJ1uHwVlNvZUtErGGLT0TpFkjBCi65BmRHQZATXASz++xDObnsGv+km3pJM7LpfTYk4LWU2O4gK2vvUSG9cVkW0uRwH8fh89/FtobLJR1eQjNTwKnSZ42fYerS5a9nNiuvYmD/mlNkraJGP0ZKVESTJGCNHlSDMiuoQaZw1Tl09ldflqAC7JvISHz3wYs/7EL3/+c7hr9pH/7r9Z9e1O8psS8KkpGBQ3iqOECrsLFRUNYDUbCagqcKA5ONGY7hGTMRYTWalWEiQZI4TooqQZEZ3e2vK1TFk+hWpnNSatiQdGPcDErIkheeP1NdrZu/BlVi9ax3ZbPK5AEgBmGiitrkbrCzYJFpP+8NN0TyCmq6oq5bUO8kts1LVKxiTFmslOiZJkjBCiyzuhZuSZZ55h/vz5VFRUMGjQIP75z38ycuTIwx67detWpk2bxvr169m7dy9///vfuffee39OzeIU4Q/4eeGHF3j+h+cJqAF6W3uTOy6XrOisDq8l4HFTuugt1n34FT9WxWD3pwAQpnGTaSymsnI3oLbEdK1h+rbTdE8gpnsgGWOnyeUFmpMxCRH0TpFkjBCi+zjuZmTBggVMmjSJ559/nlGjRvHkk08yYcIEdu7cSUJCwiHHOxwOMjMzueqqq/jLX/5yUooW3V+Vo4r7l9/P2oq1AEzMmsjUkVM7/LSM6vezb9nHbH7/AzaWmKnypgKgw0d2eAUZpmBaRh9lIkyvPSkxXa8vODNmd1mrZIxWQ3pSJJlJVknGCCG6HUVV929/OzajRo1ixIgRPP300wAEAgHS0tK4++67uf/++3/yvunp6dx7773HvTJit9uxWq3YbDYsltBfUVO0r1Vlq5i6fCq1rlrCdGE8fObDXNr70g6tYX9Md+uC19i4G4pccYCCQoBIfyl6Zz5DUsyE6Q9tDE40pnu4ZIzJoCUr2UrPHpHodZKMEUJ0Lcf6/n1cKyMej4f169czderUlts0Gg0XXHABq1evPvFqhQB8AR/PbnqWf//4b1RUsqOzyR2XS6Y1s0PrsO/YyO53X2TtD3byHYn4CTYcEYFK9I6d6FQHGkWhyeNv04woGk1zTPeO44rp2h3NyZiqA8mYSLOerBQrqXERkowRQnR7x9WMVFdX4/f76dGjR5vbe/TowY4dO05aUW63G7f7wEY9u91+0h5bdE4VTRVMWTaFDZUbALiyz5VMGTEFk850lHuePI7iAgrf/RdVP2wAVcWrpuJHS5hqw+jcgd5fDyj0iDSREn24abq/x5x+bFd/VVWV2gY3eSX1bZIxsZbgzJiEaEnGCCFOHZ0yTTN79mxmzJgR6jJEB1lespwHVjxAvbses87M9NHT+UXGLzrs+d01+9jz3otsXrYJIy6suuDyRFZYOXX1lWg9weuHxIQbSYsOa7MaYu3Vi95X345lwFnH9FyqqlJRG5wZU9fQNhmTlWIlJrLjmi8hhOgsjqsZiYuLQ6vVsm/fvja379u3j8TExJNW1NSpU5k0aVLLx3a7nbS0tJP2+KJz8Aa8/HPjP3l5y8sA9Ivpx/xx8+ll6dUhz+9rtFP04cuULFnEzvoYdjh6E6OzMzpqFwpg1PjJjGik3vHzY7r+gEpJZSP5ZTYanfuTMZCWEEmWJGOEEKe442pGDAYDw4YNY/HixUycOBEIbmBdvHgxd91110krymg0YjTKtRO6s/LGciYvm8zmqs0AXHvateSMyMGobf+v+/6Y7t7P/4fP0QhAirGWPEciHmctdQYfMebgf40kaxhJ1rCDYrrX0OPcK44ppuv1BdhTYaegVTJGp9WQkRRJZpIFk6FTLk4KIUSHOu6fhJMmTeKmm25i+PDhjBw5kieffJKmpiZuueUWAG688UZSUlKYPXs2ENz0um3btpa/l5aWsmnTJiIiIsjK6vjrRYjQ+7boWx5a+RB2j51IfSQzxszgwl4XtvvzBmO6H5G3cAFbSg24AnEMjmykrslDUV09Ud4qNPioVAzEmCMBDorpXkbKL284ppiu0+1jd7mdPRV2fP4DyZjeycGZMZKMEUKIA467Gbnmmmuoqqpi2rRpVFRUMHjwYBYtWtSyqbWoqAhNq2XrsrIyhgwZ0vJxbm4uubm5jBs3jiVLlvz8VyC6DK/fy983/J3Xt70OwIDYAcwbN4+0yPY9BReM6S4h/91X2FboYZcjGY8aXNUINBbgdTUAYNBoSIkKp4flwOrM8cZ0G5qTMcWtkzFherJSJRkjhBBHctzXGQkFuc5I11fSUMLkpZPZUrMFgN/2+y2Thk1Cr23fvRL2HRspeOcFdmyvZHtTKk3+4AZRIw6Mzp0YfJVoFYUkq4kka1irabrHF9PdPzOmotbRclusxUhWShQ9JBkjhDhFtct1RoQ4EV/v/ZppK6fR4G0g0hDJ42Me57ye57XrczqK8yl89wV2rd/BtsZUan3BU4IGxUsfcxmWQBl5DQ0/K6a7PxmTX2qjtnUyJqY5GWORZIwQQhwLaUZEu/H4PeR+n8tbO94CYGD8QOaPnU9yRHK7Pae7poK97/+b/JXfsb0hiTJ3PwAU/CRqyhgUXY1e8aOiZ3BaFKZWezes6enBabpHien6AyolVY3klx6cjImgd7KVSLOh3V6fEEJ0R9KMiHZRZC8iZ2kO22u3A3DL6bdw99C70Wva57TM/phu4bdfsaM+lj3OfgTQACph3nLC3HmgeNBERxG8rDstjUgwpnsTMSP/7ydjuodPxiikJ1ronSzJGCGEOFHy01OcdIsKFzF99XSavE1EGaOYefZMxqaObZfnCsZ036Ro0UIaGtwsqTsdrxr8tjb6awhz7UQfaMCs19EzNgJNq70bxxrTlWSMEEK0L2lGxEnj8rmYt24e7+56F4ChCUOZO3YuieEn74J4++2P6e758E1cdbUAmDQQqWmgwasnzLUTg78ak05LamwEcYdM0z16TLeh1cyY5rl1RIQFZ8akxUsyRgghThZpRsRJUWgrJGdpDrvqdqGgcNsZt3Hn4DvRaU7ut9j+mO7u916muMjGTkcyQyIbMGmCezdON+9me1kNeo1CSkwwprt/NUSj15M69jx6TrwNnSXmiM9Ra3eRd1AyJibSSHaqJGOEEKI9SDMifraPCz7mse8ew+lzEmOKYfbZsxmdMvqkP499x0Z2v/MC9QW7UIGtTX2p90WwuT6eUTFlAEQZVU5LiMASpm8b0x1xJulX3oEx4fDXNFFVlX11TvJK66m1H0jGJDYnY2IlGSOEEO1GmhFxwpw+J7PXzOZ/+f8DYETiCOacM4cEc8JJfZ5gTPdflGz8Aa0SQKeAxxcgypuHyxuH31OAM/LAALuY8ANplrjTz2iO6fY77GMHWiVjGpqTMYoCafERZKVIMkYIITqCNCPihBTUF5CzNIf8+nwUFP4w6A/8fuDv0Wq0R7/zMQrGdF+kZPUKdjfFkecYQC9jJeGefCrsLlRVJZJKYsINh5w6CcZ0b8cy4MzDPrbXF2DvvgYKymy4PG2TMZlJFsKM8l9DCCE6ivzEFcdFVVUW5i9k1ppZuPwu4sLimHPOHEYljTppz+FrtFP80csUfbuI4oYIdjT1wxkIXqK9qMmM1eFCQcViOtI03SPHdF0eHwVlbZMxRr2W3skW0hMtkowRQogQkGZEHDOH18Hj3z3Ox7s/BuDMpDOZfc5s4sLiTsrjt47pltcrbG3qjd0XDoBJ48Hg3IXBU0a4XkvPWDPWMH2rabrW5pjurw8b021weCgos1FceWgyJjU+Aq0kY4QQImSkGRHHZGftTiYvm0yhrRCNouFPg//EbWfchkb5+SsJrWO6VdUOtjelss8TBYBO8ZMVVk5mWCU1eieKEn6YmO6vSPnlbw8b061tcJFfYqO8VTImOtJIdoqVxBizJGOEEKITkGZE/CRVVXkv7z3mrp2L2+8mISyBuWPnMjxx+El57Nrvv2X3+69QW1rJTkcSRa4MVBRQA4R5SxgYtY8Ec7BhSIg8aJruEWK6+5Mx+aX11LROxkSbyUqVZIwQQnQ20oyII2r0NPLo6kf5fM/nAIxJGcOss2cRYzryNTqOlX3HBna/8yLV+fnsdvYg3zkAvxrc/GrwVhLu3oUBB6rfDBxoHn4qphsIqJRUNydjHG2TMb1TrFgkGSOEEJ2SNCPisLbXbCdnaQ5FDUVoFS33DL2Hm0+/+Wefltkf063avAGAUnccOx0pAOj8NsLdOzH460mINJIaHXXQNN3Dx3QlGSOEEF2b/JQWbaiqyoKdC5i3bh7egJfE8ETmj53P4ITBP+tx98d0y79bjsurYGzuMVJNNeTZI9C7SzH49hEbrict2tpyzRAAa3pG8zTdtjFdlyc4M6awvAGfPwC0TsZEotedvJixEEKI9iPNiGjR4GngkVWP8NXerwAYnzqex8Y8RpQp6oQf09doo/jDlyle8gVNLpXNDZk0+Y2Mi9qKTqOiRWVo+C4q/C7SEixEHkNMt9HpJb/URnFlgyRjhBCiG5BmRACwpXoLOUtzKG0sRafR8Zehf+GG/jeccNok4HFT+vmbFH3xP7xNTQDoFA313jC8qo4Cu4HTooKbS6PDDUSHGw4T070CRXfgW7S2wUV+qY3yGknGCCFEdyLNyClOVVXe2P4GT6x/Al/AR0pECvPHzueM+DNO7PH8fvYt/ZA9H71FY62NYlcsvcKaqHd4KK51YA78gDbgxKP3oEZZUeDQmO7FN6AxhrXUV1nnJK/URo3d1fI8PaLDyE6JIsZilCZECCG6OGlGTmE2t42HVz7Mt8XfAnBBzwuYMWYGFoPluB/rQEz3ZRrKyyl2xbGzaQBuVc8+WwM4awAwa+pIiQ6jh8XS0oQcLqZ7pGRMavPMGEnGCCFE9yHNyClqc9VmJi+dTHlTOXqNnpzhOfym729OaJXBvnMjuxe8QF3BLio9VrY19afRH1zZ0AaacHtchCkKiVYTydawg6bpnkX6lbe3xHS9vgBFzckYZ3MyRqtRSE+MpHeyVZIxQgjRDclP9lNMQA3w6tZXeWrDU/hUH2mRacwfN5/TY08/7sdyFOVR+O4LVP2wgXqfmW1NfajxBldVDIqXXsZS6qvySIg0HCamO5DMq+9oiem6PH4Ky20UljfgbZWMyUy2kCHJGCGE6NakGTmF1LnqeGjlQywrWQbAhPQJTD9rOhGGiON6HHdNBXvfe4Hy75bT5NOxoymDUncsAAoBMsMqyTaXo1f8eMKsbZoQa3oGva++HcvpwcF6jU4vBaU2iiobCajBaEy4SUdWipW0hAi0hxl2J4QQonuRZuQUsWHfBiYvm0yloxKDxsCUkVO4qs9Vx3VapnVM1+32k+dIotDZgwDBhsHoLSPck0/vKBP65ouj7W9EwhN6kHnFjS0x3boGN3ml9W2TMRFGslKtJEkyRgghTinSjHRzATXAf7b8h6c3Po1f9ZNuSSd3XC6nxZx27I9xUEx3jyuenU3JeNTgdFy9r5Zw906sWgdp8WYMugOrGa1jumi1VNY7yS+xUX1QMiYrJTgzRpoQIYQ49Ugz0o3VOGt4YMUDrCpbBcDFmRfz8JkPE64PP6b7t47puupqW26v9xjwqHq0/kbC3buIpJa02DDiIqytYrph9Pq/y0i5+AbQmyipbiK/tB5762RMXHMyJlySMUIIcSqTZqSbWlexjinLplDlrMKkNfHAqAeYmDXxmFYegjHdb9j9/is0lZdR5wtHq4Rh0ToB6GMux9ZQh9lXRmqUiR4WK5rmxw3GdM////buPDqqMk38+LeqUkv2IvueEAhgh0BYI7hgS5pFW8mMLWC3rSjaHg7aOjQj6IgRPQot2uNp9Ui3o9D+bNx+I9hjKw5EggoIsihbiCHEBMjCnlQqqdRy3/mjoKQgCQRJKgnP55ycQ2497837vvXeqodb9dxLWsFMCLVSUWej/PBRqYwRQgjRJnk36GU8moe/7vorS79biqY0MiMzeWHcC2T1ybqo9t676f6VU+VlAFQ0x7HbnkYYp7ghZj86ICRIY1h0PSGmSF+ZLnodiaPGkvGr+8GaRFlNAxUlB3G5z1TG6MlMjCQjMRyTVMYIIYQ4iyQjvcix5mPM/2I+m2s3A1DQv4DHRj9GiDHkgm3PLtM9o8Wt4WisBpJxu+2cbHITFeJdMhEWoy/uTJmuFt+ffdX1VJUfRNOkMkYIIcTFkWSkl9hUvYn5X87nhOMEwUHBLLh6Abf0u+WC7VqO1fjupuv2KH5wxGNzm+nj2kdtgwOlFNG6L4gJhuBzkpozZbqetFz2HK6nevsh32PWMBNZyVYSo6UyRgghRPskGenh3Jqb1757jdd3vo5CkdUnixfGvUBmZGb77c4q0/U4nVS3RLHPnkyTZgbAZi/HqJoJtxhJiwo+52668fT917vxDLyWXdUNHNtZ7Xssrk8wWVIZI4QQogMkGenB6ux1zPtyHtvqtgFwW9ZtzB89H0uQpc023jLdv1P12SpcdjvHXWHstWdyyu298JlBOQhx7CfC0ERabATWEKPf3XTTb5mOJ2ciu2oaaSg5AngrY5JPV8ZESmWMEEKIDpJkpIf66vBXPP7l45xsOUlIUAiFYwq5KfOmNuO9Zbqr+OEf7+I4eYJGj4W99n7UOfsAYNB56B9cSxSHUMHu88p0UyYU4B5ZwO4jzTSXe8t8DXod6fHh9EuOJEQqY4QQQlwieQfpYVyai1d2vMKbu98EYFDUIF4Y9wLpEemtxitN48S2db4y3RYtiO+b0qh0xKDQg9KINdQxzHoEs94FGE7/eMt0E26YjCvvdvaecOE6aAPAdLoypm9COCajVMYIIYT4aSQZ6UFqGmt49ItH+fbotwBMGziNfx/175gN5lbjzy7T9Sg9B5oTKGtKxHM62TC5jxDeUkaaVcOsD/6xoV5H9DUTcF89lVIbaEdaAG9lTL+kSNLiwjAYpDJGCCHE5SHJSA9RfLCYJzY8QX1LPWHGMBaOXciEjAmtxp5bptuiBfHFyatwKG/SEuSpJ7Tle1JCmkiJDfa7kV143nhcV0/lQLMB6r3luVIZI4QQojNJMtLNuTwuXtr+Em/tfQuA7OhsloxbQmp46nmxZ5fpKk3zbTfr3Ri0RvRKI7SljETTCdISgwk+fVl4BZhH3Ig771cccpvAe6FV4qzee8bEREpljBBCiM4jyUg3dsh2iEe/eJRdx3YBcOdVd/JvI/4Nk8G/YuXsMl3N6aTBE0ypPZGfhVQSGuS9DHtuRCXVJxpIj7MQbvZWziidHn3ueNyjpnAUC7hBByTHhp6ujGn94x8hhBDicpJkpJsqqixiwYYF2Fw2wk3hPHPNM4xPG+8Xo7U4OLx6BVWrV+JqagK8Zzm2nUqlUUXQ4nRwbYz3GiBRZo2oxNNJSJAJLfcXuHJvokXvTTh8lTFJEYScdXVVIYQQorNJMtLNOD1OXtz6Iiv2rQBgSOwQlly/hKSwJF/MuWW6bqUH9DQ0Ozh4ogmDtg+TMR2zpwJNmXw3sdPMYbhzJ+DMycej955dMQXpyUyKoG9ChFTGCCGECAhJRrqRqoYq5q6fS8mJEgBmZM/g98N/j1HvPVOhNI0TW9dx4L+XYa+tQaGjyhHDPnsSIe7DGJu+ByBEX09WSBnxEWb0Oh1aeAzO3Ek4B16L0nuf8hBLEP2TvPeMCZLKGCGEEAEkyUg3sfqH1Ty18SnsLjtWs5Vnr32W61Ou9z3eULKVAx+8zqny/SjgqDOSvfZkbB7v/WI8xBCtKyMxwkKSNZggvQ5PdCpNuZNxZY4EnTfhiAw1kZUcSWJMqO+MiRBCCBFIkowEmMPtYMk3S3j/+/cBGBY3jOevf56E0ATgTJnuXzi6cwcA9e4Q9tpTOOaKAMCoc2NxlJNiriMt1YrRoMeTfBX2oZNwp2T7/k6s1XvPGKmMEUII0d1IMhJAFfUVzF0/l+9Pej9euS/nPmbnziZIH3RemW6zZmJvYxLVzmhAhx6NjOAjZIXUoFduDIYwXJkjaBw6CS3GezVWHZAcc7oyJkwqY4QQQnRPkowEyMcHPubpTU/T7G4myhLFc9c+xzXJ1+BurKfio2UcLF6N5nThUgbKmpI40ByPOn3lVLOrhlxrLXHBHpTBhHPg9TQNmYCKiAW8lTFp8eH0l8oYIYQQPYAkI12s2d3M4i2L+bDsQwBGJYxi8XWLiTFEcPCjN31luho6Kptj2WdPxI238sXoPkGUVk6W1UWEtQ+O7J/jzB6PCg4HTlfGJEaQkRiBWSpjhBBC9BCSjHSh8lPlzF0/l/2n9qNDxwNDH+CB7Ps59sX/sOV/3sFx8qQv9gd7H/Y0ez9uMWh2Il1lZEXYiU5MwjlkAo0DrwWj96OXEJOBfilW0qQyRgghRA8kyUgXWbV/Fc9tfo5mdzPRlmgWX7uIrMNN7HjiXuy1NQC4lZ4gnfcy7ukhJymzn8LiqqF/2AniBvXDlTuZxsyRoPee9YiwGBiQFiWVMUIIIXo0SUY6WZOriWc3P8s/yv8BwNWJV/MfCdM59cZydpXvB7w3svvOlsJJVwjjo0sI0ikMOsU11lKM6dm4c++nKXWwb5/RIToG9I0nVipjhBBC9AKSjHSi709+z9z1c6mor0Cv03N/xjTGfXOcH/7+nC+mxaNRdcLGURWOpjOxv97CwD4tuPoOh6GTcMRmeAOVRkKIYuCAVKxSGSOEEKIXkWSkEyil+LDsQxZtWUSLp4VYczQPNeQQtWwdxzUNDR0Hm/tA8yFqGxwopQgz7KZPaBApeWNoHDYJLSIOAJ3HRVKwm6tysgiVyhghhBC9kCQjl5ndZWfhpoV8WvEpAMNJZup6D6FNe9GAmpY+7LYl0UIw4c2NWFQNkdY+JI+bDEPz8QR7L2amb7GTbHHws7yhWMyShAghhOi9JBm5jEqOlzB3/VyqbFXo0TGlysoN+53o0XHCFcZeewon3d475+q1FsKio+k/bgq6wTegTlfG6BuPk2pqJPuasRjNpkAORwghhOgSkoxcBkop3it9jyXfLMGpOYlyBnH3znAyG4w0eszsbEjiuCcaAAMeMpM1YkeNImhQnq8yxnDiIKmGU2TfmE+QJSSQwxFCCCG6lCQjP5HNaaNwQyFrqtYAMPioiTtLIjA6TWy3xXPYGX/6JnWKAVmhJOaNQKX9eM+YoMMlpOiP87MJN2OMiArQKIQQQojAkWTkJ9hzbA9z1j5EdctR9BpMKQ/j+qpQyuyxHGhORNMZQa8joX8Mfa8dgTExAwWgaRgrtpGs1TLo5tswxyYHeihCCCFEwEgycgmUUizb+BJ/LluGR6eIatZzz+5IdEdiWWtPwYUFvdFAcnY6aWNyMUZ5K2NwOzGVfkWCo4qBU35NSPrUwA5ECCGE6AYkGemgozWlzP/092wxVIMOhhwx85t94YS49Wy0B6MskaQPH0jyqByCQr33jNE5GjHt+ZwY2wEG/MudRPxsZoBHIYQQQnQfkoxcJHdjPWs+XMzi5k85YfEQpMGkUiujKkMJMbnQQqMYnDcO91XXYTB5K2N0tuOYd/4v1pPf07/g10SNegidXu4dI4QQQpxNkpEL0FocHPzkbZbteouVqSfRLBDTZGDc10nYjmZTEmNk2A3ZuPqPBn0QBkB/rArzd6sJO15G5i3TiL/hUXRBMtVCCCFEa+Qdsg3K46GueBW7P3mb/0qqYk+6E4DB1UaGr4eghHjSp91MVP90XKfbGA6XYP72UywnK8iYUEDyTU+iNwcHbhBCCCFEDyDJyDmUpnFi6zoO/PcydjVXsiy7gVMWDYNHx9DtiYzyZBAz418JTh3gbXC6Msb03WpM9dWkjPsFaQXPEBRuDeg4hBBCiJ7iik5GHC4Hjc5GwkxhWIwWGkq2cuD91zlxoIy16U38c5AdTQ9hthBuOH4bQyb8gpBoq7ex24mpdAOmXWsIajxGwuixZPzqGSnTFUIIITroikxGth7ayhtb32Bt+Vo0pRFv9zC9tg99TxkhNJT/N7SBkmjvxzL97IP5Vb8/ED7MewVVHHbMez7HtOdz9A4bMTlDybz9KULSBwZwREIIIUTPdUmlHa+++ioZGRlYLBby8vLYsmVLu/EffPABgwYNwmKxkJOTwyeffHJJnb0c/v7t35n+7nSKyosIb3IyeVc9t355nObv97OOAzw78hgl0U6ClJFb+8xixpBnCY+IBttxLBvfIWLFo1i2fUSfxDiGz1vM4Ll/kkRECCGE+Ak6fGbkvffeY86cOSxdupS8vDxeeuklJk6cSGlpKXFxcefFb9y4kTvuuINFixbxy1/+khUrVlBQUMD27dsZPHjwZRnExdp6aCuFawsxO92MPWAns9KBzqPQdLA7J4rduZEoHcQaU5ieMI94czocO0Twzk8xlm9FpzyExieQeds9RI26Ucp0hRBCiMtAp5RSHWmQl5fHqFGjeOWVVwDQNI3U1FQeeugh5s+ff178tGnTsNvtfPzxx75tV199Nbm5uSxduvSi/mZDQwORkZHU19cTERHRke76mf3/76dx3ccMKG/E4PQOu9kSxMYJ2dT2aQBgePh4fhn7ALpDlUTu/oSgQ3vQAWarlb63TCf+hn+RMl0hhBDiIlzs+3eH3lWdTifbtm3jscce823T6/Xk5+ezadOmVtts2rSJOXPm+G2bOHEiq1atavPvtLS00NLS4vu9oaGhI91slcPloOTbz5hQYvNu0Ok4Nfo6igfVY6cBo87MrbGz6FsVwuGV9zBEF4xepycoJJj0CQUk33QnerPlJ/dDCCGEEP46lIwcO3YMj8dDfHy83/b4+Hj27dvXapva2tpW42tra9v8O4sWLWLhwoUd6doFNTobqbQaOZpgIrbWidugp/iqU9iVjXhjOhOO/Rz9P/6TmhM/AKDS80i78SbSCmZKma4QQgjRibrl5w2PPfaY39mUhoYGUlNTf9I+w0xh6HV6NgwIo+DISYLcHn5eHsUP8THkfFaKZitGA9DrqE62MPK5/yIyse9P+ptCCCGEuLAOJSMxMTEYDAbq6ur8ttfV1ZGQkNBqm4SEhA7FA5jNZsxmc0e6dkEWo4X8fvkUlRdxKNlMykEHkV8WM/SsmGPxJr4eGMmQYZMlERFCCCG6SIfKQUwmEyNGjKCoqMi3TdM0ioqKGDNmTKttxowZ4xcPsGbNmjbjO9PMkTPRlMbG/qEog863vSHKyP9ebWXliD7Uhhm4d+S9Xd43IYQQ4krV4drUOXPm8Prrr/O3v/2NkpISZs2ahd1u55577gHgrrvu8vuC68MPP8zq1at58cUX2bdvH0899RRbt27lwQcfvHyjuEgjU0bydP7TNAQb+SEjhOZwA1+OjOS90VYORYegQ8fT+U8zMnlkl/dNCCGEuFJ1+Dsj06ZN4+jRozz55JPU1taSm5vL6tWrfV9SraqqQn/W9TfGjh3LihUreOKJJ3j88cfJyspi1apVXX6NkTN+nftrBsQOYPnXf+Wd8iI8etDr9OT3z+fekfdKIiKEEEJ0sQ5fZyQQLtd1Rs517r1phBBCCHH5dMp1Rnobi9EiSYgQQggRYHI9cyGEEEIElCQjQgghhAgoSUaEEEIIEVCSjAghhBAioCQZEUIIIURASTIihBBCiICSZEQIIYQQASXJiBBCCCECqkdc9OzMRWIbGhoC3BMhhBBCXKwz79sXuth7j0hGbDYbAKmpqQHuiRBCCCE6ymazERkZ2ebjPeLeNJqmUV1dTXh4ODqd7rLtt6GhgdTUVA4ePHhZ73nTU8l8+JP58Cfz8SOZC38yH/5kPn6klMJms5GUlOR3E91z9YgzI3q9npSUlE7bf0RExBW/YM4m8+FP5sOfzMePZC78yXz4k/nwau+MyBnyBVYhhBBCBJQkI0IIIYQIqCs6GTGbzRQWFmI2mwPdlW5B5sOfzIc/mY8fyVz4k/nwJ/PRcT3iC6xCCCGE6L2u6DMjQgghhAg8SUaEEEIIEVCSjAghhBAioCQZEUIIIURA9fpk5NVXXyUjIwOLxUJeXh5btmxpN/6DDz5g0KBBWCwWcnJy+OSTT7qop51r0aJFjBo1ivDwcOLi4igoKKC0tLTdNsuXL0en0/n9WCyWLupx53rqqafOG9ugQYPabdNb1wZARkbGefOh0+mYPXt2q/G9aW188cUX3HLLLSQlJaHT6Vi1apXf40opnnzySRITEwkODiY/P5+ysrIL7rejrz3dRXvz4XK5mDdvHjk5OYSGhpKUlMRdd91FdXV1u/u8lOOtu7jQ+pgxY8Z5Y5s0adIF99tT10dn6dXJyHvvvcecOXMoLCxk+/btDB06lIkTJ3LkyJFW4zdu3Mgdd9zBzJkz2bFjBwUFBRQUFLB79+4u7vnlt379embPns3XX3/NmjVrcLlcTJgwAbvd3m67iIgIampqfD+VlZVd1OPOl52d7Te2r776qs3Y3rw2AL755hu/uVizZg0At99+e5ttesvasNvtDB06lFdffbXVx59//nn+/Oc/s3TpUjZv3kxoaCgTJ07E4XC0uc+OvvZ0J+3NR1NTE9u3b2fBggVs376dDz/8kNLSUm699dYL7rcjx1t3cqH1ATBp0iS/sb3zzjvt7rMnr49Oo3qx0aNHq9mzZ/t+93g8KikpSS1atKjV+KlTp6qbb77Zb1teXp564IEHOrWfgXDkyBEFqPXr17cZs2zZMhUZGdl1nepChYWFaujQoRcdfyWtDaWUevjhh1W/fv2UpmmtPt5b1wagVq5c6ftd0zSVkJCglixZ4tt26tQpZTab1TvvvNPmfjr62tNdnTsfrdmyZYsCVGVlZZsxHT3euqvW5uPuu+9WU6ZM6dB+esv6uJx67ZkRp9PJtm3byM/P923T6/Xk5+ezadOmVtts2rTJLx5g4sSJbcb3ZPX19QBERUW1G9fY2Eh6ejqpqalMmTKFPXv2dEX3ukRZWRlJSUlkZmbym9/8hqqqqjZjr6S14XQ6efvtt7n33nvbvTFlb14bZ1RUVFBbW+v33EdGRpKXl9fmc38prz09WX19PTqdDqvV2m5cR463nqa4uJi4uDgGDhzIrFmzOH78eJuxV9r6uFi9Nhk5duwYHo+H+Ph4v+3x8fHU1ta22qa2trZD8T2Vpmk88sgjXHPNNQwePLjNuIEDB/Lmm2/y0Ucf8fbbb6NpGmPHjuXQoUNd2NvOkZeXx/Lly1m9ejWvvfYaFRUVXHfdddhstlbjr5S1AbBq1SpOnTrFjBkz2ozpzWvjbGee344895fy2tNTORwO5s2bxx133NHuDeE6erz1JJMmTeKtt96iqKiIP/7xj6xfv57Jkyfj8Xhajb+S1kdH9Ii79orLa/bs2ezevfuCn9mOGTOGMWPG+H4fO3YsV111FX/5y1945plnOrubnWry5Mm+fw8ZMoS8vDzS09N5//33mTlzZgB7FnhvvPEGkydPJikpqc2Y3rw2xMVxuVxMnToVpRSvvfZau7G9+XibPn267985OTkMGTKEfv36UVxczPjx4wPYs56l154ZiYmJwWAwUFdX57e9rq6OhISEVtskJCR0KL4nevDBB/n4449Zt24dKSkpHWprNBoZNmwY+/fv76TeBY7VamXAgAFtju1KWBsAlZWVrF27lvvuu69D7Xrr2jjz/Hbkub+U156e5kwiUllZyZo1a9o9K9KaCx1vPVlmZiYxMTFtju1KWB+XotcmIyaTiREjRlBUVOTbpmkaRUVFfv+jO9uYMWP84gHWrFnTZnxPopTiwQcfZOXKlXz++ef07du3w/vweDzs2rWLxMTETuhhYDU2NlJeXt7m2Hrz2jjbsmXLiIuL4+abb+5Qu966Nvr27UtCQoLfc9/Q0MDmzZvbfO4v5bWnJzmTiJSVlbF27Vqio6M7vI8LHW892aFDhzh+/HibY+vt6+OSBfobtJ3p3XffVWazWS1fvlzt3btX/e53v1NWq1XV1tYqpZT67W9/q+bPn++L37BhgwoKClIvvPCCKikpUYWFhcpoNKpdu3YFagiXzaxZs1RkZKQqLi5WNTU1vp+mpiZfzLnzsXDhQvXZZ5+p8vJytW3bNjV9+nRlsVjUnj17AjGEy+oPf/iDKi4uVhUVFWrDhg0qPz9fxcTEqCNHjiilrqy1cYbH41FpaWlq3rx55z3Wm9eGzWZTO3bsUDt27FCA+tOf/qR27Njhqw5ZvHixslqt6qOPPlI7d+5UU6ZMUX379lXNzc2+fdx4443q5Zdf9v1+odee7qy9+XA6nerWW29VKSkp6ttvv/V7LWlpafHt49z5uNDx1p21Nx82m03NnTtXbdq0SVVUVKi1a9eq4cOHq6ysLOVwOHz76E3ro7P06mREKaVefvlllZaWpkwmkxo9erT6+uuvfY+NGzdO3X333X7x77//vhowYIAymUwqOztb/fOf/+ziHncOoNWfZcuW+WLOnY9HHnnEN3fx8fHqpptuUtu3b+/6zneCadOmqcTERGUymVRycrKaNm2a2r9/v+/xK2ltnPHZZ58pQJWWlp73WG9eG+vWrWv12DgzXk3T1IIFC1R8fLwym81q/Pjx581Renq6Kiws9NvW3mtPd9befFRUVLT5WrJu3TrfPs6djwsdb91Ze/PR1NSkJkyYoGJjY5XRaFTp6enq/vvvPy+p6E3ro7PolFKqC07ACCGEEEK0qtd+Z0QIIYQQPYMkI0IIIYQIKElGhBBCCBFQkowIIYQQIqAkGRFCCCFEQEkyIoQQQoiAkmRECCGEEAElyYgQQgghAkqSESGEEEIElCQjQgghhAgoSUaEEEIIEVCSjAghhBAioP4PEwMVOxfU260AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_valsp = np.hstack([0,x])\n",
    "plt.scatter(x,L,alpha=0.5,c=\"rebeccapurple\",s=100)\n",
    "plt.scatter(0,0,s=50,c=\"forestgreen\")\n",
    "plt.plot(x_valsp,b_star*x_valsp,label = 'GL',linestyle = 'dashed')\n",
    "plt.plot(x_valsp,b_oracle*x_valsp,label = 'Oracle_Cov', alpha = 0.5,lw = 5)\n",
    "plt.plot(x_valsp,beta1*x_valsp,label = 'Truth',c=\"cornflowerblue\",linestyle=\"dashdot\")\n",
    "plt.plot(x_valsp,b_r*x_valsp,label=\"Crude\",c=\"maroon\",lw = 5,alpha = 0.5)\n",
    "plt.plot(x_valsp,intercept_ols_slope*x_valsp,label=\"Intercept\",c=\"lightsteelblue\")\n",
    "plt.plot(x_valsp,b*x_valsp,label = \"OLS\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "correlationCorrection",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}