{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "078900c6-d71d-44f0-9e76-cc5ce66f0ca9",
   "metadata": {},
   "source": [
    "# Logistic regression: using general fit\n",
    "\n",
    "Exploring different fitting techniques and using MLE asymptotic approx of parameter distribution to estimate model CI intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e18cec5e",
   "metadata": {},
   "source": [
    "## Common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3b351af7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import os\n",
    "\n",
    "# our libs\n",
    "import data_utils\n",
    "import logit_utils_gen\n",
    "\n",
    "np.set_printoptions(precision=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0de7f02-942b-42ff-874c-3c0e2ac1e0dc",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3373f87f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_path = \"./results\"\n",
    "if not os.path.exists(results_path): os.makedirs(results_path)\n",
    "\n",
    "data_path = \"../../data/\"\n",
    "suv_filename = os.path.join(data_path, \"suv_percentilesSLOthenUWM.mat\")\n",
    "flags_filename = os.path.join(data_path,\"flags_combined.mat\")\n",
    "normal_range_filename = os.path.join(data_path, \"normal_range.mat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "811d990d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Loading all data\n",
    "\"\"\"\n",
    "suv_dict = scipy.io.loadmat(suv_filename)\n",
    "flags_dict = scipy.io.loadmat(flags_filename)\n",
    "normal_range_dict = scipy.io.loadmat(normal_range_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4a292c31",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get data\n",
    "perc = 95\n",
    "organ = \"lung\"\n",
    "\n",
    "x0, y = data_utils.get_data(organ, perc, suv_dict, flags_dict)\n",
    "\n",
    "scales = [\"plain\", \"log\"]\n",
    "xs = [x0, np.log(x0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e9bb7b-5a50-4cae-9af4-261c11453dd4",
   "metadata": {},
   "source": [
    "## Fit\n",
    "\n",
    "Ref: logistic regression documentation\n",
    "* https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression\n",
    "* https://stats.stackexchange.com/questions/186830/what-is-scikit-learns-logisticregression-minimizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "239ae789",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "scale",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "pars",
         "rawType": "object",
         "type": "unknown"
        },
        {
         "name": "cost",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "success",
         "rawType": "bool",
         "type": "boolean"
        }
       ],
       "ref": "c29dae99-1d2d-4ec3-ad24-d349a78fccd9",
       "rows": [
        [
         "plain",
         "[-24.952835540044568    0.7197040288775033   1.3756115410561167\n  -4.742671756961299 ]",
         "6.958343476437286",
         "True"
        ],
        [
         "log",
         "[-2.2050225858748224e+01  1.6444831830139661e-06  9.0151795404966251e+00\n -8.5482042264396672e+00]",
         "5.608180093799321",
         "True"
        ]
       ],
       "shape": {
        "columns": 3,
        "rows": 2
       }
      },
      "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>pars</th>\n",
       "      <th>cost</th>\n",
       "      <th>success</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>scale</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>plain</th>\n",
       "      <td>[-24.952835540044568, 0.7197040288775033, 1.37...</td>\n",
       "      <td>6.958343</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log</th>\n",
       "      <td>[-22.050225858748224, 1.644483183013966e-06, 9...</td>\n",
       "      <td>5.608180</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                    pars      cost  success\n",
       "scale                                                                      \n",
       "plain  [-24.952835540044568, 0.7197040288775033, 1.37...  6.958343     True\n",
       "log    [-22.050225858748224, 1.644483183013966e-06, 9...  5.608180     True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fits\n",
    "lg = logit_utils_gen.LogisticPolyRegression(3, mono = True, lam = (0, 1e-3))\n",
    "ress = [lg.fit(x, y, method=\"diff_evol\") for x in xs]\n",
    "\n",
    "pd.DataFrame(ress, index = pd.Index(scales, name = 'scale'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d9044f0b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pars:\n"
     ]
    },
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "scale",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "0",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "1",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "2",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "3",
         "rawType": "float64",
         "type": "float"
        }
       ],
       "ref": "9dff318b-17a2-4baf-ab23-0932b8097b53",
       "rows": [
        [
         "plain",
         "-24.952835540044568",
         "0.7197040288775033",
         "1.3756115410561167",
         "-4.742671756961299"
        ],
        [
         "log",
         "-22.050225858748224",
         "1.644483183013966e-06",
         "9.015179540496625",
         "-8.548204226439667"
        ]
       ],
       "shape": {
        "columns": 4,
        "rows": 2
       }
      },
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>scale</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>plain</th>\n",
       "      <td>-24.952836</td>\n",
       "      <td>0.719704</td>\n",
       "      <td>1.375612</td>\n",
       "      <td>-4.742672</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log</th>\n",
       "      <td>-22.050226</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>9.015180</td>\n",
       "      <td>-8.548204</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               0         1         2         3\n",
       "scale                                         \n",
       "plain -24.952836  0.719704  1.375612 -4.742672\n",
       "log   -22.050226  0.000002  9.015180 -8.548204"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parss = np.array([r[\"pars\"] for r in ress])\n",
    "\n",
    "print(\"pars:\")\n",
    "pd.DataFrame(parss, index = pd.Index(scales, name = 'scale'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9a99322a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------\n",
      "pars=[np.float64(-24.952835540044568), np.float64(0.7197040288775033), np.float64(1.3756115410561167), np.float64(-4.742671756961299)]\n",
      "beta=[np.float64(-24.952835540044568), np.float64(23.010909283460883), np.float64(-6.524074004316852), np.float64(0.6307690372955947)]\n",
      "nllf = 6.310796258553427, jac = [np.float64(0.04990566392214985), np.float64(-0.001439408016133653), np.float64(-0.0027512214980024985), np.float64(0.009485341724485175)]\n",
      "------------------------\n",
      "pars=[np.float64(-22.050225858748224), np.float64(1.644483183013966e-06), np.float64(9.015179540496625), np.float64(-8.548204226439667)]\n",
      "beta=[np.float64(-22.050225858748224), np.float64(73.07179549692368), np.float64(-77.06359585018566), np.float64(27.091154049129646)]\n",
      "nllf = 4.9676223757332, jac = [np.float64(0.044100490125984805), np.float64(1.5062468017637075e-07), np.float64(-0.01803014829444692), np.float64(0.017096548916171614)]\n"
     ]
    }
   ],
   "source": [
    "for x, pars in zip(xs, parss): \n",
    "    r = lg.nllf(x, y, pars, jac = True)\n",
    "    beta = lg.get_beta(pars)\n",
    "    print(\"------------------------\")\n",
    "    print(f\"pars={list(pars)}\")\n",
    "    print(f\"beta={list(beta)}\")\n",
    "    print(f\"nllf = {r[0]}, jac = {list(r[1])}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9cf6fec1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cov_pars = array([[88.66981718045956  , -1.2461162929855938, -7.2032252157241405,\n",
      "        12.643615608184488 ],\n",
      "       [-1.2461162929796519,  5.875091109892306 ,  1.086039589044909 ,\n",
      "         0.2688880961870854],\n",
      "       [-7.203225215732767 ,  1.08603958904545  ,  0.9368560858745143,\n",
      "        -1.105349253065405 ],\n",
      "       [12.643615608183422 ,  0.2688880961869379, -1.1053492530639295,\n",
      "         1.9662718284625142]])\n",
      "cov_pars = array([[ 4.0189309259208073e+01,  1.5283870335397518e-05,\n",
      "        -1.0539487960613590e+01,  8.4965481264772382e+00],\n",
      "       [ 1.5283870335398561e-05,  1.0456138304267554e+01,\n",
      "        -1.5255470929999245e-06,  5.7878951772164524e-06],\n",
      "       [-1.0539487960613663e+01, -1.5255470929995728e-06,\n",
      "         5.8583150404942561e+00, -3.0423947835076719e+00],\n",
      "       [ 8.4965481264773377e+00,  5.7878951772162592e-06,\n",
      "        -3.0423947835076688e+00,  2.0221168062093726e+00]])\n"
     ]
    }
   ],
   "source": [
    "# asymptotic covariance matrix of parameters\n",
    "cov_parss = np.array([lg.cov(x, y, pars) for x, pars in zip(xs, parss)])\n",
    "\n",
    "for cov_pars in cov_parss: print(f\"{cov_pars = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cbc4c71e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[9.1084353421439957e+01 6.0574470013571613e+00 3.0447099119140469e-01\n",
      " 1.7647907003263518e-03]\n",
      "[4.5051054919554510e+01 3.0080446689362188e+00 1.0641517426858928e-02\n",
      " 1.0456138304261696e+01]\n"
     ]
    }
   ],
   "source": [
    "# check eigenvalues\n",
    "for cov_pars in cov_parss: print(np.linalg.eigvals(cov_pars))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2359d3c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# defining two-sided confidence intervals\n",
    "alpha = 0.05   # significance level\n",
    "probs = [alpha/2, 1 - alpha/2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4ba4d5d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pars_CI = array([[-43.40876822677736  ,  -4.030971772150962 ,  -0.5214635957804892,\n",
      "         -7.491008028113301 ],\n",
      "       [ -6.496902853311777 ,   5.4703798299059665,   3.272686677892722 ,\n",
      "         -1.9943354858092968]])\n",
      "pars_CI = array([[-34.47542511475169 ,  -6.337728605324418,   4.271291263280167,\n",
      "        -11.33529562100979 ],\n",
      "       [ -9.625026602744763,   6.337731894290783,  13.759067817713081,\n",
      "         -5.761112831869545]])\n"
     ]
    }
   ],
   "source": [
    "# CI of params (assuming asymptotic distr of parameters) : Wald approximation\n",
    "pars_CIs = np.array([lg.get_pars_quantiles_normal(probs, pars, cov_pars) for pars, cov_pars in zip(parss, cov_parss)])\n",
    "\n",
    "for pars_CI in pars_CIs: print(f\"{pars_CI = }\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "96acc3ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "index",
         "rawType": "int64",
         "type": "integer"
        },
        {
         "name": "scale",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "coef",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "value",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "LCL",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "UCL",
         "rawType": "float64",
         "type": "float"
        }
       ],
       "ref": "bcfba9b0-99bf-4501-a183-fca1e53e25d6",
       "rows": [
        [
         "0",
         "plain",
         "p0",
         "-24.952835540044568",
         "-43.40876822677736",
         "-6.496902853311777"
        ],
        [
         "1",
         "plain",
         "p1",
         "0.7197040288775033",
         "-4.030971772150962",
         "5.4703798299059665"
        ],
        [
         "2",
         "plain",
         "p2",
         "1.3756115410561167",
         "-0.5214635957804892",
         "3.272686677892722"
        ],
        [
         "3",
         "plain",
         "p3",
         "-4.742671756961299",
         "-7.491008028113301",
         "-1.9943354858092968"
        ],
        [
         "4",
         "log",
         "p0",
         "-22.050225858748224",
         "-34.47542511475169",
         "-9.625026602744763"
        ],
        [
         "5",
         "log",
         "p1",
         "1.644483183013966e-06",
         "-6.337728605324418",
         "6.337731894290783"
        ],
        [
         "6",
         "log",
         "p2",
         "9.015179540496625",
         "4.271291263280167",
         "13.759067817713081"
        ],
        [
         "7",
         "log",
         "p3",
         "-8.548204226439667",
         "-11.33529562100979",
         "-5.761112831869545"
        ]
       ],
       "shape": {
        "columns": 5,
        "rows": 8
       }
      },
      "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>scale</th>\n",
       "      <th>coef</th>\n",
       "      <th>value</th>\n",
       "      <th>LCL</th>\n",
       "      <th>UCL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>plain</td>\n",
       "      <td>p0</td>\n",
       "      <td>-24.952836</td>\n",
       "      <td>-43.408768</td>\n",
       "      <td>-6.496903</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>plain</td>\n",
       "      <td>p1</td>\n",
       "      <td>0.719704</td>\n",
       "      <td>-4.030972</td>\n",
       "      <td>5.470380</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>plain</td>\n",
       "      <td>p2</td>\n",
       "      <td>1.375612</td>\n",
       "      <td>-0.521464</td>\n",
       "      <td>3.272687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>plain</td>\n",
       "      <td>p3</td>\n",
       "      <td>-4.742672</td>\n",
       "      <td>-7.491008</td>\n",
       "      <td>-1.994335</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>log</td>\n",
       "      <td>p0</td>\n",
       "      <td>-22.050226</td>\n",
       "      <td>-34.475425</td>\n",
       "      <td>-9.625027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>log</td>\n",
       "      <td>p1</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>-6.337729</td>\n",
       "      <td>6.337732</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>log</td>\n",
       "      <td>p2</td>\n",
       "      <td>9.015180</td>\n",
       "      <td>4.271291</td>\n",
       "      <td>13.759068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>log</td>\n",
       "      <td>p3</td>\n",
       "      <td>-8.548204</td>\n",
       "      <td>-11.335296</td>\n",
       "      <td>-5.761113</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   scale coef      value        LCL        UCL\n",
       "0  plain   p0 -24.952836 -43.408768  -6.496903\n",
       "1  plain   p1   0.719704  -4.030972   5.470380\n",
       "2  plain   p2   1.375612  -0.521464   3.272687\n",
       "3  plain   p3  -4.742672  -7.491008  -1.994335\n",
       "4    log   p0 -22.050226 -34.475425  -9.625027\n",
       "5    log   p1   0.000002  -6.337729   6.337732\n",
       "6    log   p2   9.015180   4.271291  13.759068\n",
       "7    log   p3  -8.548204 -11.335296  -5.761113"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# nr of parameters\n",
    "d = parss.shape[-1]\n",
    "\n",
    "# present table of results\n",
    "df_pars_CI_wald = pd.DataFrame({\n",
    "    \"scale\" :[lab for lab in scales for _ in range(d)],\n",
    "    \"coef\" : [ f\"p{i}\" for _ in scales for i in range(d)],\n",
    "    \"value\" : parss.flatten(), \n",
    "    \"LCL\" : pars_CIs[:,0].flatten(), \n",
    "    \"UCL\": pars_CIs[:,1].flatten()\n",
    "})\n",
    "\n",
    "df_pars_CI_wald"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6ef6f46b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.microsoft.datawrangler.viewer.v0+json": {
       "columns": [
        {
         "name": "scale",
         "rawType": "object",
         "type": "string"
        },
        {
         "name": "LLF",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "AIC",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "BIC",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "A",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "chi2",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "p-value(chi2)",
         "rawType": "float64",
         "type": "float"
        },
        {
         "name": "n",
         "rawType": "int64",
         "type": "integer"
        },
        {
         "name": "k",
         "rawType": "int64",
         "type": "integer"
        },
        {
         "name": "dof",
         "rawType": "int64",
         "type": "integer"
        }
       ],
       "ref": "b01a6f6b-bc09-4d1d-a601-170cf17e85d7",
       "rows": [
        [
         "plain",
         "-6.310796258553427",
         "20.621592517106855",
         "28.86336455929253",
         "0.9655172413793104",
         "13.812457417951638",
         "0.9999999944452796",
         "58",
         "4",
         "54"
        ],
        [
         "log",
         "-4.9676223757332",
         "17.9352447514664",
         "26.177016793652076",
         "0.9655172413793104",
         "8.60216708864235",
         "0.9999999999998123",
         "58",
         "4",
         "54"
        ]
       ],
       "shape": {
        "columns": 9,
        "rows": 2
       }
      },
      "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>LLF</th>\n",
       "      <th>AIC</th>\n",
       "      <th>BIC</th>\n",
       "      <th>A</th>\n",
       "      <th>chi2</th>\n",
       "      <th>p-value(chi2)</th>\n",
       "      <th>n</th>\n",
       "      <th>k</th>\n",
       "      <th>dof</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>scale</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>plain</th>\n",
       "      <td>-6.310796</td>\n",
       "      <td>20.621593</td>\n",
       "      <td>28.863365</td>\n",
       "      <td>0.965517</td>\n",
       "      <td>13.812457</td>\n",
       "      <td>1.0</td>\n",
       "      <td>58</td>\n",
       "      <td>4</td>\n",
       "      <td>54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log</th>\n",
       "      <td>-4.967622</td>\n",
       "      <td>17.935245</td>\n",
       "      <td>26.177017</td>\n",
       "      <td>0.965517</td>\n",
       "      <td>8.602167</td>\n",
       "      <td>1.0</td>\n",
       "      <td>58</td>\n",
       "      <td>4</td>\n",
       "      <td>54</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            LLF        AIC        BIC         A       chi2  p-value(chi2)   n  \\\n",
       "scale                                                                           \n",
       "plain -6.310796  20.621593  28.863365  0.965517  13.812457            1.0  58   \n",
       "log   -4.967622  17.935245  26.177017  0.965517   8.602167            1.0  58   \n",
       "\n",
       "       k  dof  \n",
       "scale          \n",
       "plain  4   54  \n",
       "log    4   54  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# goodness of fit measures\n",
    "pd.DataFrame(\n",
    "    [lg.goodness_of_fit(x, y, pars) for x, pars in zip(xs, parss)], \n",
    "    index = pd.Index(scales, name = 'scale')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "54f463d0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHIAAAHaCAYAAABsPG7IAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAA+H1JREFUeJzs3XecXHW5+PHPmbazvffd7G42bdMLJCQBEiCQEHqvUpRmR36+BK5KVVC8cvWKilxRLIioV1FEQQG5FpDQFUJJr9t3dno95/v74+xusskm2TKz0543r2GzM2fOeXYzmfOd53y/z6MppRRCCCGEEEIIIYQQIuVZkh2AEEIIIYQQQgghhBgdSeQIIYQQQgghhBBCpAlJ5AghhBBCCCGEEEKkCUnkCCGEEEIIIYQQQqQJSeQIIYQQQgghhBBCpAlJ5AghhBBCCCGEEEKkCUnkCCGEEEIIIYQQQqQJSeQIIYQQQgghhBBCpAlJ5AghhBBCCCGEEEKkCUnkiKywevVqVq9enewwjuiqq66ioKAg7vt9+umnWbhwIU6nE03T6O/vj/sxstlVV11Fc3NzUo6taRp33HHHmJ/3wgsvoGkaL7zwQtxjEkKIVCRjARkLJFK8xwJ33HEHmqbFbX9CiMwiiRwhMlxvby8XXnghubm5fPvb3+YnP/kJ+fn5k3Ls73znO2iaxrJlyw65jaZph7zdcMMNkxKnEEIIkclkLCCEEJnFluwAhBCJ9corr+D1ern77rtZs2bNpB770Ucfpbm5mQ0bNrB582amTZs24nYnn3wyV1xxxUH3z5gxI9EhZq3jjz+eYDCIw+FIdihCCCESTMYCQgiRWSSRI0QChUIhHA4HFkvyJr91dXUBUFJSErd9+v3+I17J27ZtGy+++CK//vWvuf7663n00Ue5/fbbR9x2xowZXH755XGLTxyZxWLB6XQmOwwhhMh4MhaQsYAQQsSbLK0SCbNnzx4+8pGPUFdXR05ODi0tLXz0ox8lEokMbbN161YuuOACysrKyMvL45hjjuGpp54atp/BWh6/+MUv+PKXv0xDQwNOp5OTTjqJzZs3H3Tchx56iNbWVnJzc1m6dCl/+9vfJvRzrF69mrlz5/Laa6+xYsUKcnNzaWlp4cEHHxwxzp///Od84QtfoL6+nry8PDweDwC//OUvWbJkCbm5uVRUVHD55ZezZ8+eEY+5detW1q5dS35+PnV1ddx1110opcYV+5VXXgnA0UcfjaZpXHXVVUOPjyamwbX6W7ZsYf369RQWFnLZZZcd8diPPvoopaWlnHbaaZx//vk8+uijY45/f4FAgPfee4+enp4jbrtp0ybOO+88ampqcDqdNDQ0cPHFF+N2u4dt99Of/pSlS5eSl5dHaWkpxx9/PH/605+GHv/tb3/LaaedNvQabm1t5e6770bX9SPGYBgG3/jGN5gzZw5Op5Pq6mquv/56XC7XEZ87+Dsfz+tgx44dfOxjH2PmzJnk5uZSXl7OBRdcwPbt24dtN1KNnMHX+saNGznhhBPIy8ujvr6e++6774gxCyHESGQsIGMBGQuMbywwklgsxt13301rays5OTk0NzfzH//xH4TD4YOOe8cdd1BXV0deXh4nnHACGzdupLm5edjfvRAifcmMHJEQe/fuZenSpfT393Pdddcxa9Ys9uzZw69+9SsCgQAOh4POzk5WrFhBIBDgU5/6FOXl5fzoRz/izDPP5Fe/+hXnnHPOsH1+5StfwWKx8NnPfha32819993HZZddxssvvzy0zcMPP8z111/PihUruPHGG9m6dStnnnkmZWVlNDY2jvvncblcrF+/ngsvvJBLLrmEX/ziF3z0ox/F4XDw4Q9/eNi2d999Nw6Hg89+9rOEw2EcDgePPPIIV199NUcffTT33nsvnZ2dfPOb3+Qf//gHb7zxxrArZLqus27dOo455hjuu+8+nn76aW6//XZisRh33XXXmOL+/Oc/z8yZM3nooYe46667aGlpobW1FWBMMcViMdauXcuxxx7Lf/7nf5KXl3fEYz/66KOce+65OBwOLrnkEr773e/yyiuvcPTRRx+0bSgUGnFQVlRUNLT0Z8OGDZxwwgncfvvthy3uG4lEWLt2LeFwmE9+8pPU1NSwZ88efv/739Pf309xcTEAd955J3fccQcrVqzgrrvuwuFw8PLLL/P8889zyimnDP2OCgoKuOmmmygoKOD555/ntttuw+Px8LWvfe2wP//1118/9Dv+1Kc+xbZt23jggQd44403+Mc//oHdbj/s88f7OnjllVd48cUXufjii2loaGD79u1897vfZfXq1WzcuPGIf3cul4t169Zx7rnncuGFF/KrX/2Km2++mXnz5nHqqace9rlCCLE/GQvIWEDGAhMbCxzommuu4Uc/+hHnn38+/+///T9efvll7r33Xt59911+85vfDG136623ct9993HGGWewdu1a3nrrLdauXUsoFBrT8YQQKUwJkQBXXHGFslgs6pVXXjnoMcMwlFJK3XjjjQpQf/vb34Ye83q9qqWlRTU3Nytd15VSSv3lL39RgGpra1PhcHho229+85sKUP/+97+VUkpFIhFVVVWlFi5cOGy7hx56SAFq1apV4/pZVq1apQD19a9/fei+cDisFi5cqKqqqlQkEhkW59SpU1UgEBjadjCuuXPnqmAwOHT/73//ewWo2267bei+K6+8UgHqk5/85LDf12mnnaYcDofq7u4ec/w//OEPFTDs72I8Md1yyy2jPuarr76qAPXnP/956GdoaGhQn/70pw/aFjjk7bHHHhvabvD3e/vttx/22G+88YYC1C9/+ctDbrNp0yZlsVjUOeecM/Q6GzT4+lRKDft7HHT99dervLw8FQqFhu678sorVVNT09D3f/vb3xSgHn300WHPffrpp0e8/0BjeR0c+DsZKeaXXnpJAerHP/7x0H2Dv8+//OUvQ/cNvtb33y4cDquamhp13nnnHTZmIYQ4kIwFZCwgY4HxjwVuv/12tf9HtTfffFMB6pprrhm23Wc/+1kFqOeff14ppVRHR4ey2Wzq7LPPHrbdHXfcoQB15ZVXHva4Qoj0IEurRNwZhsETTzzBGWecwVFHHXXQ44OtFP/whz+wdOlSjj322KHHCgoKuO6669i+fTsbN24c9ryrr756WGHW4447DjCnHgO8+uqrdHV1ccMNNwzb7qqrrhq68jJeNpuN66+/fuh7h8PB9ddfT1dXF6+99tqwba+88kpyc3OHvh+M62Mf+9iwmiSnnXYas2bNOmj6OMAnPvGJoT9rmsYnPvEJIpEIzz777IR+jonE9NGPfnTU+3/00Ueprq7mhBNOAMyf4aKLLuLnP//5iFORzzrrLP785z8fdBt8PphTw5VSR2y1Pfh3/cwzzxAIBEbc5oknnsAwDG677baDahbs3+pz/79Hr9dLT08Pxx133NDU7kP55S9/SXFxMSeffDI9PT1DtyVLllBQUMBf/vKXw/4Mg8bzOtg/5mg0Sm9vL9OmTaOkpITXX3/9iMcsKCgYVqPA4XCwdOnSoX9nQggxGjIWkLGAjAXiMxYY9Ic//AGAm266adj9/+///T+Aob+v5557jlgsxsc+9rFh233yk58c0/GEEKlNEjki7rq7u/F4PMydO/ew2+3YsYOZM2cedH9bW9vQ4/ubMmXKsO9LS0sBhtYZD24/ffr0YdvZ7XamTp06hp/gYHV1dQcV9BvsonBg7ZGWlpZh3w/GNdLPOmvWrIN+TovFclC8hzrWeI01JpvNRkNDw6j2res6P//5zznhhBPYtm0bmzdvZvPmzSxbtozOzk6ee+65g57T0NDAmjVrDrpVV1eP+WdraWnhpptu4vvf/z4VFRWsXbuWb3/728PWxG/ZsgWLxcLs2bMPu6933nmHc845h+LiYoqKiqisrBxKchy4xn5/mzZtwu12U1VVRWVl5bCbz+cbKjp5OON9HQSDQW677TYaGxvJycmhoqKCyspK+vv7DxvzoIaGhmEDWDD/rY13Pb8QIjvJWEDGAjIWmPhYYH87duzAYrEc1PWrpqaGkpKSob+vwa8HbldWVjb070UIkf6kRo5IG1ardcT71TgK/yXS/lduMkVOTs6ou208//zztLe38/Of/5yf//znBz3+6KOPDq07T5Svf/3rXHXVVfz2t7/lT3/6E5/61Ke49957+ec//znqQWh/fz+rVq2iqKiIu+66i9bWVpxOJ6+//jo333wzhmEc8rmGYVBVVXXIoo6VlZXj+rlG45Of/CQ//OEPufHGG1m+fDnFxcVomsbFF1982JgHpcu/MyFEdkqX9ygZC8hYIFFjgQMvtgghspMkckTcVVZWUlRUxNtvv33Y7Zqamnj//fcPun9wmmpTU9OYjju4/aZNmzjxxBOH7o9Go2zbto0FCxaMaX/727t370FtNj/44AMAmpubRxXX+++/PyyuwfsO/DkNw2Dr1q1DV97GcqzRGmtMY/Hoo49SVVXFt7/97YMe+/Wvf81vfvMbHnzwwYQPcufNm8e8efP4whe+wIsvvsjKlSt58MEH+dKXvkRrayuGYbBx40YWLlw44vNfeOEFent7+fWvf83xxx8/dP+2bduOeOzW1laeffZZVq5cOe6fc7yvg1/96ldceeWVfP3rXx+6LxQK0d/fP644hBBiPGQsMHJcMhaQscB4NTU1YRgGmzZtGpqxBtDZ2Ul/f//Q39fg182bNw+bGdbb2yuza4XIILK0SsSdxWLh7LPP5sknn+TVV1896PHBq2br169nw4YNvPTSS0OP+f1+HnroIZqbm4841fVARx11FJWVlTz44IPD2po+8sgjI36Ife+999i5c+eo9h2Lxfje97439H0kEuF73/selZWVLFmy5IhxVVVV8eCDDw5rD/nHP/6Rd999l9NOO+2g5zzwwANDf1ZK8cADD2C32znppJOG7t+yZQtbtmwZVfzxiGk0gsEgv/71rzn99NM5//zzD7p94hOfwOv18rvf/W7M+x5ty1GPx0MsFht237x587BYLEM/69lnn43FYuGuu+466Gra4Otz8Krv/ld5I5EI3/nOd44Y64UXXoiu69x9990HPRaLxUadVBnN6+BAVqv1oCvT3/rWt0bVJlUIIeJFxgIHxyVjARkLDBrLWGDQ+vXrAfjGN74x7P77778fYOjv66STTsJms/Hd73532Hb7v56EEOlPZuSIhLjnnnv405/+xKpVq7juuutoa2ujvb2dX/7yl/z973+npKSEW265hccee4xTTz2VT33qU5SVlfGjH/2Ibdu28b//+7+jnr47yG6386UvfYnrr7+eE088kYsuuoht27bxwx/+cMR18W1tbaxatYoXXnjhiPuuq6vjq1/9Ktu3b2fGjBk8/vjjvPnmmzz00ENHbB1pt9v56le/ytVXX82qVau45JJLhtp7Njc385nPfGbY9k6nk6effporr7ySZcuW8cc//pGnnnqK//iP/xg2DXdwIDeetfJjjWm0fve73+H1ejnzzDNHfPyYY46hsrKSRx99lIsuumjo/g8++ICf/vSnB21fXV3NySefDIy+5ejzzz/PJz7xCS644AJmzJhBLBbjJz/5CVarlfPOOw8w141//vOf5+677+a4447j3HPPJScnh1deeYW6ujruvfdeVqxYQWlpKVdeeSWf+tSn0DSNn/zkJ6Oavr9q1Squv/567r33Xt58801OOeUU7HY7mzZt4pe//CXf/OY3Of/88w+7j9G+Dg50+umn85Of/ITi4mJmz57NSy+9xLPPPkt5efkR4xZCiHiSscDwuGQsYJKxwOjHAvtbsGABV155JQ899NDQkq8NGzbwox/9iLPPPnuoKHR1dTWf/vSn+frXv86ZZ57JunXreOutt/jjH/9IRUWFLM0SIlNMep8skTV27NihrrjiClVZWalycnLU1KlT1cc//vFh7UC3bNmizj//fFVSUqKcTqdaunSp+v3vfz9sP4OtJg9sIblt2zYFqB/+8IfD7v/Od76jWlpaVE5OjjrqqKPUX//6V7Vq1aqDWo4yyjakq1atUnPmzFGvvvqqWr58uXI6naqpqUk98MADo4pz0OOPP64WLVqkcnJyVFlZmbrsssvU7t27h21z5ZVXqvz8fLVlyxZ1yimnqLy8PFVdXa1uv/32g1pjNjU1DWtzeSgjtRwdT0yjccYZZyin06n8fv8ht7nqqquU3W5XPT09SqnDtxzd/+9ntC1Ht27dqj784Q+r1tZW5XQ6VVlZmTrhhBPUs88+e9C2P/jBD4Z+/tLSUrVq1aqhNqlKKfWPf/xDHXPMMSo3N1fV1dWpz33uc+qZZ545qG33gS1HBz300ENqyZIlKjc3VxUWFqp58+apz33uc2rv3r2H/RnG8jo48HficrnU1VdfrSoqKlRBQYFau3ateu+991RTU9OwlqOHaj8+Z86cEeMZzWtNCCEOJGOB4WQsYJKxwJHHAge2H1dKqWg0qu68807V0tKi7Ha7amxsVLfeeuuwNuhKKRWLxdQXv/hFVVNTo3Jzc9WJJ56o3n33XVVeXq5uuOGGwx5XCJEeNKVSrDqcEClm9erV9PT0HHGdvxDxctVVV/GrX/0Kn8+X7FCEEEIgYwGR/vr7+yktLeVLX/oSn//855MdjhBigqRGjhBCCCGEEEJkiGAweNB9g7V1Vq9ePbnBCCESQmrkCCGEEEIIIUSGePzxx3nkkUdYv349BQUF/P3vf+exxx7jlFNOYeXKlckOTwgRB5LIEUIIIYQQQogMMX/+fGw2G/fddx8ej2eoAPKXvvSlZIcmhIgTqZEjhBBCCCGEEEIIkSakRo4QQgghhBBCCCFEmpBEjhBCCCGEEEIIIUSayMoaOYZhsHfvXgoLC9E0LdnhCCGEECil8Hq91NXVYbHIdZaJknO9EEKIVCPnehEvWZnI2bt3L42NjckOQwghhDjIrl27aGhoSHYYaU/O9UIIIVKVnOvFRGVlIqewsBAw/wEVFRUlORohhBACPB4PjY2NQ+coMTFyrhdCCJFq5Fwv4iUrEzmDU6yLiopkcCeEECKlyDKg+JBzvRBCiFQl53oxUbIwTwghhBBCCCGEECJNSCJHCCGEEEIIIYQQIk1IIkcIIYQQQgghhBAiTUgiRwghhBBCCCGEECJNSCJHCCGEEEIIIYQQIk1IIkcIIYQQQgghhBAiTUgiRwghhBBCCCGEECJNSCJHCCGEEEIIIYQQIk1IIkcIIYQQQgghhBAiTdiSHYDIUoYOO14EXycUVEPTCrBYkx/T9r/Dtr+BBjQdCy3HHTquVPwZhNifvEaFENks2e+BoxlXJDvGbCC/YyFEBkp6Iuevf/0rX/va13jttddob2/nN7/5DWefffZhn/PCCy9w00038c4779DY2MgXvvAFrrrqqkmJV8TBxt/B0zeDZ++++4rqYN1XYfaZyYvpyU9DsG+/O78GuaVwxn8fHFcq/gxC7E9eoyKFyLleTLpkvweOZlyR7BizgfyOhRAZKulLq/x+PwsWLODb3/72qLbftm0bp512GieccAJvvvkmN954I9dccw3PPPNMgiMVcbHxd/CLK4afUAE87eb9G3+XpJg+dMBga0DQZT62f1yp+DMIsT95jYoUI+d6MamS/R44mnHFn74o79OJluzXgRBCJJCmlFLJDmKQpmlHvEp3880389RTT/H2228P3XfxxRfT39/P008/ParjeDweiouLcbvdFBUVTTRsMVqGDt+Ye/AJdYhmXiW58d+TN+X1iDENKKyDzwy85lLtZxBif6n470yMSracm+RcLxIq2e+Box1XaBZQxqEelPfpiUr260CIQ5Bzk4iXpC+tGquXXnqJNWvWDLtv7dq13HjjjYd8TjgcJhwOD33v8XgSFZ44nB0vHmFgo8Czx9yu5bgUiWmAd6+5LaTezyDE/lLx31kaMgxFMKrjD8cIRHSCUZ1ARCcU1QnHdEJRY+jrsdMqaCzLS3bIGUXO9WLckv0eONpxxSGTOCDv03EQp9eBUooeX4QdvX56fBH6/BF6fWH6g9Ghc0IwohOK6cR0RcwwiOkKXSkMZT7fUAqlMG8D+xTpZ+2cGj5z8oxkhyHEkLRL5HR0dFBdXT3svurqajweD8FgkNzc3IOec++993LnnXcmJiCPB77+dXC7E7P/TNKzCTaHjrzd5q9Axf8mPh44KCaPTfHzGoONBSOcZL95hfm133fk/X7zCsgtiU+MQoxFsF9eowOUghg5hMknQt7A13wi5BIljwhOouQRxTlwyyGGkxg5xLAz2tXHtyyJcsMFZyf0Z8k2KXeuF6MTDsNLL0FfH7hc0NNjjo9iMTAM8zb4iTZRkjDWeM/Sxy8cm+jTQqN/Dx6NLHifTphxngvDKo8uptHDVLxU4aWCKJKoF6D735BEjkgpaZfIGY9bb72Vm266aeh7j8dDY2NjfHbu88H27eagJCcnPvvMVFENIoe7ArXfdn5/4uMZPNZATDtyFT9oNPhXEZRHFNpBG+vmF20UPwM6xAJxDFSI0dKz4jVqYCFMMSFKCWklA38uIawVE6aQCEWEKcTQHBM7kDKwEsaqIlhVBAsRrESxqCgWouiWCNZwVXx+KDEhCT3XiyMLheCRR+C558wxkaaZ4yK7HSwDSVFNM2+JNIljjRgGf3bu4VeObXhUlBLdzujfg0cjvd+nk2v058JQzM4e7Rh6tDY8NJjL3vanDJy4yMGLAx925cNBACvhoXOCQRi/JUyECLloFBgWcpSGAwsONBxKw4IGmONLm7JgVaDpOpaojkXXQTfQDAPNYgGLdeDrwM1qMeOyWMCimX/WtP3+DKCZjw/++VAG/x0Objf4b1Jjv3+f2r5tGX7XYfd94HMyhc3K6pXrkx2FEMOkXSKnpqaGzs7OYfd1dnZSVFQ04hU6gJycHHISnWSpqoKCgsQeI93V1UL/PyHsw5xceiANcgph1uKDT6IJjEn1v8QrDg8/rlF0OKAtAA51wPGtOVDSYv7Z1wOxCIf8GWyOgW0z8EQm0oDKiNeoUhAmF48qxatKzJtRjE8V41XFBClAjXLGjE1FyFE+nMqPU/nJIUAOIXII4tBCOLQIDi2MwxrDbo1ht+nYbQZ2u4HVpqHZrGC1DgymrfsNrK28E9jOzCU3J/i3kX1S9lwvRhYIwMMPw1/+As3NkMy6D5M01ugjyKP8i7+yk2LymU0h2sAHdfM9OHyEPQy+/6bv+3RqO/K5sD9nKm/nX81WfQ4G++rklGpd1Fu2U27poFjro0hzYdNiIx4lSJTt9KOhMZM8ZlPJTCpopoRK8ijAgaYwk4Z+v3kBOBjETDRp4MyF3FzIz4eKCqishPJy8zNFfj7k5Zm33FxwOs3kqNMJDse+BKkQIiulXSJn+fLl/OEPfxh235///GeWL1+epIjEqGkWmHYyvPMbGBzs7HvQ/DJtzeQlcQZierFtEd8LP48BzPYzcMXkAOVT98VYNhW63uOQP0PZftsKMem0tHqNGsqCR5XgNspwG6XmV1WGh3Ii2sgf2AdZVIx83OTjIU/zkWcLkG8LkuuIkJcTJdcRwZkTxe7QwGYzkzD7f7XZwFayLzkznquIXT5zkC3iSs71acTng//5H/jb32Dq1ORf1JqEsYZC8W028DrtTKWUQvZPIA6+B797+J0U14N7z6FjTKH36fR06HOh31LBP4uuY5fzmKHJ1tXaLqZZ36beuo087cgztULE2EE/BjCHKtYxjQVUk4PNXELY3w/uvWbyRikzKZOfDzNmmP9OqqvNhE1pKZSVQWFhZs5kEUIkTNITOT6fj82bNw99v23bNt58803KysqYMmUKt956K3v27OHHP/4xADfccAMPPPAAn/vc5/jwhz/M888/zy9+8QueeuqpZP0IYiwqZ8Kcc2DznyHs3Xd/TqE5sKqcOekh/bvESTBcz5zObiA6/EGLHSqmQV7FvvvyKqBqFvRtHX7FzeYwBw37bytEMqTga1QphU8vwBUrx6VX4KKSfirxWCoxtANORfuNZfM0L0V2L4UOP4U5QQqcQQrzIuTnR8nNA81h3+/KpAXIH7iJVCLn+gxlGPC978Hf/w7TppkfVFNBgscaPQTYTj9TKD4giTMgrwKq2qBnMxiHGVfkFKbU+3TGGeFcuD1nJS8Wf5yIpQBQNFk2Mce2gSpL+6h36ybEDtzMpYpTmcZR1GEPRaG300zgKAXFxdDQALNnQ1OT+eeaGnM2jRBCxEHSEzmvvvoqJ5xwwtD3g+vbr7zySh555BHa29vZuXPn0OMtLS089dRTfOYzn+Gb3/wmDQ0NfP/732ft2rWTHrsYp8qZUDEd+ndBxA+OfChpnNyZOAMUiq24KMypgClTIeSGoNv8IJlTDLnFjHhFLK8C8srN7fUoWO3gPMS2QiRDMl6jhgGxGHrMoD9WRq+qpk+roc9Sh8taS3Rwhs0BnV5tWpRip4/ivADFBWGKCyMUF0UpzI9gsx5Y48ACOAduIl3IuT5DuVywcSM0NqZOEmdQAscau/HgIUIthYfeKK8CppQfflwhY4nEG/gdR0IhXjbWs8VyFAAVWjvH2v9IiaV3TLvrJUAHPtYzncvDM3F2u8C10awFVVkJp5xiJm+mTze/l1k2QogESXoiZ/Xq1Ydtw/fII4+M+Jw33ngjgVGJhNMsUNqU7CjoJ4SLIAU4AA2cJeZtVLQxbCtEMiTiNapANyAagUiUWFThoopeSwO91gZ6bQ24tGqUw3rQMy2aQXFegNKiIKWFIUoKgpQUBCnIjchYN8PJuT5DuVxmkePKymRHMrIEjTV240FhYDtina7RjCtkLJFoPlXMM9q1eC2laBjMs/6ThbaXsIyxKHUnPnpVgLO9dVy0y47NstOcZXPCCTBvnrlsSmbcCCEmSdITOUIkUyd+fESokNaSQoxAQTRmthSOhDFiCpelhh5bEz32KfQ6GnA5KlEHTrEBHLYY5UUBygZupYUBivNDWC0JbDsshJhcLpf5/pBlRaa34MI2wvueSD0+VcjTkYvwqRLycXO84ymqLXvGvJ92vR9v0M0lPVWcZW/GsuZoWLoUZs0yCxELIcQkk0SOyGod+Iig45ABmRBg6OaHsnAYFYngtVTSbW+lx95ET8EU+lQVujr4tOF0RCkv8lNeFKC82Pya75RZNkJkPJfL/JpF/9gj6Gymj6KRauOIlOJXBTwTuRifKqFQc7HO8XPyNd/YdqLrBDy99Fk9XOFYyunnfxjt6KPNYsVCCJFEksgRWa0DH9rAf0JkFwUx3VwWEQ4RNHLpsTXR7WihJ7eJntw6Imq/DyoDM9DtthgVxX4qigJUFPspL/ZL0kaIbNXXl+wIJt1evLgJUU2Su3OJw/KrAp6OXIxXlVCo9bPO8fjYkjiGDm4PKhphS7nBMTNOZ925X0PLP0xdJCGEmESSyBFZbSsunPLPQGSFgWVSoRDRsEGvVjewRKqJ7oIp+FXRsE0BLBaD8oGEjXnzUZQXlqSNEMLU3m52jcsiu/EQJEY+9mSHIg4hquz8KXIBXlVKgdbPWsfPyde8R34iAAp8frNteGkJexvrKa+p4cITbsUuSRwhRAqRT7Aia4WJsRvPQKFjITJQLIYeiuCKlNCr1ZtLpGxT6LdXovYv0qnM/5UUhKgo9g0lbsoKg1ikpo0QYiRKmYmcLCvuugsPgMzkTVFKwUvRU3CrCvLwstb+OAWjTeJEItDvAmcuzJ1LqLkBl38XV80+l6aS5DfoEEKI/UkiR2StLvz4iVBJirVMFWKcDAP6Q4X0Rivo1WrpsU6hz1qLkXfwW31eToSKEnOWTeVAXRuHXU9C1EKItBQMgseTVYVeFYr36SFPZuOkrM36XLYas9EwWOV4kkKLZxTPUtDvhlgUGhqhrQ1VUsLm7ndYVLOIk1tPTnjcQggxVpLIEVmrEz9+ojTJgEykIV1Z6Vfl9MYq6Y2W06tqcFnr0C12DqzB6bDHqCjy77dEyk+eM5qcwIUQmcHlMpM5FRXJjmTSeAizF68UOk5R/UY5/4ytAWCR7e+j606l69DbC/l5sGABNDWBxUKnr4OinCIunHMhOTb5+xZCpB5J5Iis1YEPUFhkerRIcSGVi8uooE9V0WdU0WdU0q8qUNpAt7X9mq7ZrbGB7lFmB6mKYj+FuVLXRggRZy6XWSw9i2bk7MaDhzBTKU12KOIAMWXjhegZ6Nips2xnnvXlIz8pHIL+fqiqhoULodT8e1VK0enr5MK5FzK9fHpC4xZCiPGSRI7IWrtwY92/TogQSRZTNtyqDJeqpN+owKUqcBmVBBihwKIGDkuY8nwv5WVhykpClBf5pRixEGJyuFzmbAZb9gwl9+Alho5j/+y5SAkbYifQrypx4uc4+1NHPg/6vBAMQes0mDcPcvbNuukL9lHiLOHYxmMTG7QQQkxA9px9hdiPQrEVF/lS6FgkQVTZcasy3KqcfqOcflVOv6rAp4qHFyHeT6HRS6mtl7LiIGUVOmU1Gvl5MUnaJFGHrwOH1UG+XepsiSzkciU7gkm3DRcWuQCUcrqMOj7QFwJwvP335GqBw2ytwO02i8otWADTpoFl+N9ph6+D45uOp6GoIXFBCyHEBEkiR2QlFyFchKRjlUgYQ2n4VDEeVYJHleFRZbhVKW6jjABFh3xeDgFKjA5Ko3sp1XooLfBRWu/AXl0KxcVI5iY17PHswR/1c9Gci1hUuyjZ4Qgx+bq7s+r9SMfgfXoplHFDSjGUxj+jZl2c6dZ/UWfdeZitFfT1gdUGixab9XAOeA0HogGsFivHNR2HlkWvbyFE+pFEjshKnfjwE6FKOlaJCYgpKz5VjFeVDLt5VCk+VYxxmOn3TvwUa32UWHopoZuSyG6KAzvJtQTQioqgpRFqaqCkLqs+LKWDne6dRPUoVy64klNaT5HBvshOe/dmVX2cTvz0EaSU7Gq3nuo+0OfTp6qxE2Kx7W+H3lAps6hxrtNM4tTVjbjZXu9eppVNY27V3ARFLIQQ8SGJHJGVOvETwcAuU6TFYRhKI0ABPlVsJmyM4n1/ViUEKIDDFMu2EKNIc1GsuSjS+sw/W/oo1vrI0YIQCoPbZ07xzs+HaXXm4LKiAqxSgyEVbXNtQ9M0PrL4I6xqWiVJHJGddB26usCZPUmNPXjwE6HxMDMqxeQKKSdvxI4DYLHt74deUqUU9HRDYSEsOQoqK0fcLGbECMaCnNB8AjaLfEQSQqQ2eZcSWWkvXjRAk45VWc1QFvyqEJ8qGrgV46doIFlThF8Voo5Q1NJGhEKtnyKtn8KBW5HmosjiIg/vwZNplAE+PwQC4HCYs24aG6G6Oqs+FKWjXe5dAFy7+FpWTlmZ5GiESCK322w9np89s1p348EAaZKQQl6PHU+YXEq1LmZa3xx5I2VAd4+5NPnoo6Gs7JD76/R1UltQy9H1RycmYCGEiCNJ5IistI1+cuXln/F0ZR1I1BQPJWf2fS0a6AZ1+GSeBZ18zUOB5qZQc1MwcDOTNm5yCI5u5VM0Al4vxGKQXwBtbdDQYLY7lVkdKa/D10EwFuQjiz4iSRwhXC4zkVNenuxIJs0m+nBKt6qU0WNU84E+H4Bj7M9i0dTBGxmGOROntBSOXgolJYfcn1KK3mAvF865kKIcmXUlhEh98klWZJ0wMfbgkY5VGUApCJK/rz6NUYJX7Vv+dKSlT2AufyrQPAM391DSZvD7XPwjDxBHG2AgAH6/uVSqrMwsrlhXJ7Nv0khvoBdXyMWlcy/lxJYTkx2OEMnnckE4nDXvYzEMduGWcUOKUAo2RE8ENKZaNlJt2XPwRoNJnLJycyZOcfFh9+kKuSjOKWZ5w/LEBC2EEHEmiRyRdTrx45NCx2lDKQiRh2eg45NHleJRZQPJm2JiRxhY24gMJWUKNM9BiRongfhPiNFj5uybcBjy8qC11Zx9U1l5UJtTkdrcITcdvg7OmXUOZ8w8Q2riCAH7Wo9nyb+HHgL4iFAihY5TQrvRRJdqwEKMJfb/O3gDpaCnx7x4smyZWRvnCLr8XRxddzRTiqckIGIhhIg/SeSIrNOJjwBR8rAnOxSxn8HZNf1GBS5VQb8qx21U0K/KiHDozigaBvmaZ2ipU8HQV/fYlj5N/Ccwixd7vea3JcUwew7U12dVHYlMEoqF2OHewanTTuXCuRdi0SQJJwRgdv/JIp348BGhXgodJ51S8HrsWABmWd8kX/MduIXZYrygwCxsPIokjlKKqB5lYc1CSdYLIdKGJHJE1unCD4BFCh0njaE0PKqUPlVNr1FNn6qiz6gkTN4hnqEo0NwDHaDM7k+FmosirZ8CzYNFMyY1/mEMA/w+CAQhx2EmbpqazOLFNnmLTVe6ofNB7wccXXc0l82/TDqYCLG/9nazWHuW6MKPQmGTQsdJt8dooUfVYSXKXNuGgzdwu8Fug8WLzdo4o+AOuynKKWJWxaw4RyuEEIkjI1ORddrxSRJnkgVUPt1GHd1GLT2qlh6jZsQlURoGhVo/pVo3JVovxZZeSrReijQXNi2WhMgPw9DB44FIBAoKYfZsc/lUSUnWLDfIVEopPuj9gKmlU/nwog+TZz9UglGILKSUmcjJkvo4AB34OFK9NZF4SsEbA7Nx2qxvkKf5h2/g94Guw+IlZkfIUeoN9DK1dCp1hXXxDFcIIRJKEjki62ynX5ZVJZjXKKLDaKRLNdBhNOJVB18VsxGhTOumzNJJudZJqcVM3qRcwuZAug4eN0RjZtJm3nwzgZOTk+zIRJzs8uyiMKeQDy/6MNUF1ckOR4jUEgyaSezcQy95zTTbZNyQEnYZrfSqGmxEDp6NEwmDPwBz55qzYkdJKYU/6ueouqNkWZUQIq1IIkdklQBRegjIgCzOIspBu9HEXqOZvUbTCIkbRanWQ6VlLxVaO5WWdoq13vF3g0oGwzATOJGIWUBxais0NoJdXkuZpC/YRzAa5Nol1zK7cnaywxEi9Qy2Hq+oSHYkkyJEjHa85Mu4IamGz8Z5HacW3PegYZivy5apMGvWmGbFBmNBcm25zKyYGe+QhRAioSSRI7JKN34CRKmWjlUT5jWK2WW0ssuYRofRgMI69JiGToXWQY1lF9WW3VRa9pKjhZMY7QQoZRYwDgbNAsbzF0gCJ0OFY2F2uXdxTts5nNB8QrLDESI1uVwQCmXNjJwu/PiJSqfLJNthzMClqrATZo7tlf0eUWbx7YoKmDdvzJ0hewI91BTW0FLSEt+AhRAiwSSRI7JKDwGCRMmVK2vj4jWK2WbMZJvehktVDXusSOul3rKdOssOaiw7sWvRJEUZJ8oAr89M4BQUwIIF0NIiS6gylFKKTX2bWFizkPPazpMp9kIcistlLjHNkmLuXQMXgGQmb/IoBf+KHQNAm/U1nFpo34Mer3leXrBwXMlFd9jNmqlrsFvl71cIkV6y4ywsxIBuAiikY9VYhFUOW/XZbNFn06P2FQLUMKi27KLRsoVGyxaKLP3JCzKeBosYhyNm29L582HKFGkhnuF2eXZRnlvOhxZ8iHyH/F0LcUguV7IjmFTS6TL52o0m+lQ1NiLMtr2274Fw2LwtWjSupX4RPYJNs9FW0RbHaIUQYnJIIkdkFbPzhDgSpaDDmMImfR7bjRkYA28VGgY1ll20WN5linXT8Kti6S4WMxM4sZi5hGruPLOIcRZ1ZslWnrAHX8THdUuuY2rp1GSHI0Rq6+7Oqs58e/FKEifJ/q0vA2C69d/7xh3KgH6XWa9u6vjet3sDvVTkVTCjfEa8QhVCiEkjiRyRVaRj1eHFlI0t+hw26otxq31Xt0q1bqZb/0Wz9f2D232mu2jULGJsqIEixlPNBI7UwMkKMSPGNtc2Tm49WeriCDEae/dmTX0chWIbLhk3JFGPUU270YSGwWzrq/secPVDaanZpWqMdXGGdhF0cdLUk2QWphAiLUkiR2SNMDE68UnniREEVR7vxhbzvr6QMOYA3U6YFuu7zLD+m3KtI/MuwIbD5gwciwblFdDaCrW1WVP3QZi29G1hetl0Lp57MVaL9chPECLbud1Zk+j2EaGHAAU4kh1K1no7thSAFst7FFo85p3hsDkjp232uJOKuqFjYDC3am68QhVCiEkln1hE1ugmQIAoFeQlO5SUEVR5vB07mvf1hcQGBqoFWj9t1teZbv03Di2S5AgTIBQ0iyPabFBXZ87Aqa4e9xU9kb76gn1YLVYumnsRZbllyQ5HiPQQjYI1O5Kegx2r6pEltsngMUrYYZjLnubaNph3KgX9/dDcbM6eHaf+UD8lzhJpOy6ESFuSyBFZw2w9HpMp0kBYOflXbNmwBE6F1s4828s0WjZj0VSSI0yAcAjcHvNKcnOz2YGqoiKraj2IfWJGjF3uXZw+43QW1y5OdjhCpAfDMOuIZUniuws/QaI4ZbicFO/oR6GwUG/ZSpml27zT44bCApg9e0Lnb1fIxczymZTnlscpWiGEmFxyZhJZw+xYpbCSHQPQkRjKwvv6At6IrSQysISqQmtnoe0f1Fu2ZWZOIxI2lwJYbWb3qWnToLxcEjhZbptrG9PKpnFu27nSalyI0YrFzGROFiVyADQpdjzpgiqPzbq57GmudWA2TjQKkYjZjKCwcGL7jwWZUzVH3v+FEGlLEjkia3Thy+qh2G69hVdiJ+BW5tWnEq2bJba/0mDZmpk5jUjYnIFjsUBdPUyfDpWVksARuIIuFIoL5lxAaW5pssMRIn0MJnKypJbYLjzYyY5lZKnm3dgidOxUaHupsewCFLhc5pLo5uYJ7Vs3dCyahSnFU+ISqxBCJEN2nImFAHbgzsrp0UGVx8vRE9lutAGQQ4DFtr8z3fqvDF1CNVjE2GIO+KZNg6oqSeAIwBzA7/TsZP209Rxdd3SywxEivWTRjByFYjv90iAhCWLKxvv6QgDm2l4xT9+BAOQ4oK1twjWa3GE3RY4imoqbJh6sEEIkSfZ9qhVZKYJOO76sqo+jFGw12tgQPZEweQOtO19jge3FDCxirCA0kMAZLGIsCRwxgu3922kuaea82efJlHohxioWM08uWZDIcRGinxD50rFq0m3R5xAmjwKtnymWTeZrzuuFWW3m0ugJcofcNBY3UpFXEYdohRAiOSSRI7JCLwF8RChjfG0q001I5fL36Dp2G9MAKNW6WGl/mgpLZ5IjS4BwaKAdrsOsgTN1qiyhEiPyR/yE9TDnzDpHulQJMR66njUzcsyOVREqpdPlpFLKLHIMMNv6mjlz2O2BwiJziXQc+KN+ZlfOlmS+ECKtSSJHZIXB1uMNFCU7lITrMBr4a+R0AhRiIcYC20vMs27AohnJDi2+YlGzBanFAk3N5gycsjJJ4IgRKaXY6trKMQ3HsKJxRbLDESI9ZdHSqi78RDBwSI2cSbXbmIpHlWEnxHTrv0GPmRds5syBvIkn1QxljoVaSlomvC8hhEgmSeSIrNCNHwMDWwZ3rDKUxr/0Y3grtgKFhWKth9X2Jym19CQ7tPjSdbP9aCwGVdUwcyZUV0sCRxxWp7+TEmcJ580+D5tFTn1CjEuWJXI0lHSsmmSDs3FmWv+FXRu4YFNRCU3xqWfjDXspdBTSVCL1cYQQ6U1GsyIrmC1EM3cwFlY5vBA9k3ajGYBp1n+zzPacOQjKFIMJnGgUSkth+gxobJxw0UOR+aJ6lC5/FxfPvZippVOTHY4Q6SuLEjl7pGPVpOs1qugwmtDQabO9bnafRDMv2NjjU+PQHXZTkVdBTUFNXPYnhBDJIokckRUyuWOVxyjh2ei5eFQ5NiIst/+ZVuvGZIcVP7GoWcQ4pkNZKbROg4aGuA3qRObb3r+d6WXTOXXaqckORYj0liWJHIViD15ys6hBQip4J2bOxmmxvE++5oF+t3nBpq4ubsfwhD0c33Q8Fi2zX8NCiMyXmZ9shdhPDIM9eDKyY1WnUc/zkbMJk0c+Hk5y/JoyS3eyw4qPcBi8HvPPZeVmEeOGBrMrlRCj5A17MZTBOW3nUJhTmOxwhEhvWZLICRDFTYhcGSZPGr8qYJsxC4DZtlchGASHA2bMiNvSaaUUCiUzM4UQGUHOUCLj9RHET5QicpIdSlxt1Wfx9+ipGNgo19o5yfEb8jR/ssOauFBoXxvxmlpoaYGaGllCJcZMKcWO/h2snLKSpfVLkx2OEOkvSxI5fQQJEqOK/GSHkjXejS1GYaVa20WFpQM83n1NDOLEH/WTb8+nqVjq4wgh0p8kckTG6x5oIVpLQbJDiZsPYvN4MbYW0GiyfMBx9qewabFkhzUxgwkcu90sajh1KlRUSBFjMW7dgW6KnEWcOfNMmUYvRDxkVSInKjNyJklEOXhfXwjAXNsr4PdDbq6ZyInjGKA/1E9Zbhn1RfVx26cQQiSLnKFExushQAwjY4oWvhdbyD9jJwMwy/o6y2zPpXGuQ0EgCD7fvgROayuUl0sCR0yIoQzave2cN/s8Wstakx2OEJkhNnDBIMPfn/sIogBrBne6TCWb9PlEyaFY66VB2ww+P8yeDcXFcT2OO+Tm6LqjpXOhECIjyDuZyHi9BMmUjlUbY4vZEDsJgNnWVzna9pf0HE8rZV5x8/vB6TRn3zQ3SwJHxM0u9y4aixulwLEQ8aTryY5gUvQRTHYIWcNQFjbGlgAwx/oKms8LBQXmuCCOlFIYymBa2bS47lcIIZJFEjki4+3Fiz0Drqq9E1vCK7ETAZhrfZkltr+mX87D0MHrhVAY8vOgrc2chRPnq24iu4VjYbwRLxfPvZjyvPJkhyNE5oil+RLeUerAhy0Dxg3pYJsxEz9FOPEzVXvbLHK8YAHkx7c+UTAWxGlzMqV4Slz3K4QQySKJHJHRFIrdGdCxaos+eyiJs8D6Igtt/0ivJE4sBh63+bW4GGbOMluKxnmgJgTAtv5ttFW0sap5VbJDESKzxGJZMWtyL16pjzMJlIK3Y2Yh+jbb69i8LnOM0Nwc92N5wh6KncU0FjfGfd9CCJEMcpYSGS1AFBfBtB6Q7dWb+Ht0HWAup1pk/0eSIxoDPQZutzkdv7zC7EBVX2+2FBUiAbxhLxbNwpkzz8RpcyY7HCEySxbMyAkSpY8guWl+ASgdtBtNuFQVNiLM0l6HSBjmzjWXXMeZN+xlXvU8OS8IITJG+n66FWIUeggQJEZ1mrYQ7TWqeD56NgorzZZ3Odr2l2SHNDp6DNwec9BfXm52nmhokBbiIuF2uHdwTMMxLKlbkuxQhMg8WZDIcREiSIxycpMdSsZ7Wz8agOnWf5Pj7TJn4zQmZsZMKBZiRvmMhOxbCCGSQRI5IqP1EiRANC2vrHmNYp6NnEcMBzWWHRxn/2Pqz2iPhM0W4gooLYFp080Ejk3eakTi9QX7yLfnc/qM06XduBCJEIuZ62Ey2L7W40XJDiWj9RmV7DVa0DCYbdmwbzZOTk7cj2UoA03TqCusi/u+hRAiWeTTlchovQQAhSXNulZFlZ3no+cQpIBSrYsT7U9g1VK0W4hSEAyYLcRtNqiuMZdQ1dRIAkdMGqUUuz27Wdu6lpnlM5MdjhCZKQtm5PQRREdJseME+1fsGACaLO9T6N0FJSUwJTGFiP0RP/n2fOoL6xOyfyGESAb5lCUyWg8B0q31uFLwUvRkXKoSJ37WOP4XhxZJdlgH0wc6UIXDkJsLU1vNDlQVFVlRDFOklk5/J2W5Zayfvh5NXn9CJEaWJHJEYvUbZWw3zIT7AsuLEInA3HkJq5/njXgpchZRU1CTkP0LIUQySCJHZLTdeMghveqyvK8vZKsxBw2D1Y7fka/5kh3ScLGBAsaGbq5nb2szCxgXFCQ7MpGlDGXQ6e/kwtkXSkcSIRJJT9GZoXHUiQ9rml0ASjf/ii0HNKZYPqDUtyWhs3HALHTcVtGG3Zp+y+yFEOJQJJEjMpaBYg/etGo93mXUsmGgzfgS21+psexOckT7GUrgGFBeBq3ToK4O7Onz+xWZaY9nD/WF9ZzSekqyQxEis0VScHZonO3Bm5Z19dKFxyhhmzELgAWWf5ivqXnzE9rNMqJHaC1rTdj+hRAiGSSRIzKWmxA+ImmTyAmpXF6InImBlSbL+8yxvpLskEx6DPoHEjgV5fsSOFL/RqSAmBGjP9TPObPOoTyvPNnhCJHZIhGwZG7tmAg6PQRwyvA4Yf6lH4PCQoNlC+W+TeZsnAR1qgLQDR2LZpFCx0KIjCNnKpGxzNbjUSrIS3YoR6QU/CO6jgBFFGm9rLQ/nfwyM7oO7n7za3mF2UK8vl5aiIuUstO9k+bSZla3rE52KEJkvgxP5LgGOlaV4Ex2KBnJaxSxRZ8NwPxJmo3ji/gocBRIoWMhRMaRRI7IWL0ECaOnRY2cLcYcdhnTsKCz2v5kcosb67q5hCoWM5dQTZtuJnBkBo5IMRE9QiAaYP209RQ4pEaTEAkXDmd0IsdsPR6jNk1m8qabf+vLUFiptWynyvdewmfjgFnouCS3hKr8qoQeRwghJpt8MhMZa7D1uJbiRQv9qpCXoycBsND2d8os3ckJ5MAETus0aGiQBI5IWTv6dzC9bDorGlckOxQhskOGz8jpI0gMA7u0Ho87nypksz4X2K82TgI7VQ3yhr0srl2M1ZL6F/WEEGIs5BOayFhd+FM+iWMuqVpLlBwqtL3MTUZdHGO/BE6ZJHBEeghGg0SNKKfNOI1ce26ywxEiO2RBIgdI+bFDOnojehwGNmosO6jxbTS7XiawU9Ug3dBpLmlO+HGEEGKyySc1kbF24Un5zhMf6AvYa7RgJcpx9j9g0dTkHdwwwOOBSBhKy8waOA0N0oVKpIUd7h3MrpzNsvplyQ5FiOwRjWZ0IqcLPxZJ4sRdr1HFFsOsjXOU5S/muGPu3ITPxokZMSwWi9THEUJkpJQ4G3/729+mubkZp9PJsmXL2LBhw2G3/8Y3vsHMmTPJzc2lsbGRz3zmM4RCoUmKVqSDKDqd+FK6Y5XXKOaV2GoAFtv+RrHFNUlHVuD3QXc35DphyVGwahW0tEgSR6QFf8ScbXfa9NOwW+U1my7kXJ8BMjyRswevdKyKM6Xg1dhqQGOqZSMVvk3mbJwE18YBc1lVoaNQOlYJITJS0s/Gjz/+ODfddBO33347r7/+OgsWLGDt2rV0dXWNuP3PfvYzbrnlFm6//XbeffddHn74YR5//HH+4z/+Y5IjF6msjyABouSm8IBsQ+xEYjio1nYx2/ra5Bw0GoGubojpMGcOrFoNra0JvyomRDzt6N/BgpoFLK5dnOxQxCjJuT5DZPDSqhgGXfhTetyQjvYYLbQbTViIscj6f+ZsnKmtkJOT8GN7I17Kc8upyKtI+LGEEGKyJf1sfP/993Pttddy9dVXM3v2bB588EHy8vL4wQ9+MOL2L774IitXruTSSy+lubmZU045hUsuueSIV/ZEdukd6DyRqjNydust7DKmoaGz3P6nSWg1rsxW4v1uqKuD444zpzU7pcWqSC+esAe71c766euleGUakXN9hsjgGTmuoQtAqTluSEeG0ng1tgqANusbFPp2Q1HRpNTGAXNGzoyKGWiJH2QJIcSkS+rZOBKJ8Nprr7FmzZqh+ywWC2vWrOGll14a8TkrVqzgtddeGxrMbd26lT/84Q+sX7/+kMcJh8N4PJ5hN5HZeggMdJ5IvQ96urKyIXYiALOtr1Ni6UvsAWMx6OoyixcfdRSsWGEWNRYiDe1072Rx7WLmVs1NdihilORcnyEMwzyfZGgip48gIWIyIyeOtuhz6FeVOAgy3/oPCE3ebBwAQxlMKZ6cpJEQQky2pJ6tenp60HWd6urqYfdXV1fz3nvvjficSy+9lJ6eHo499liUUsRiMW644YbDTre+9957ufPOO+Mau0htvQRQTGLh4DHYqC/Bo8rIxccC24uJPZjfB34/1NbB/PnmunQh0lR/qJ9cWy6nTj8Vi5aZHyYzkZzrM0QsZiZzMjSR4yJEBB1HCl4ASkcR5eCN2LEAzLf9kxxfDxQWTtpsnIgewW61S6FjIUTGSruz8QsvvMA999zDd77zHV5//XV+/etf89RTT3H33Xcf8jm33norbrd76LZr165JjFgkQzs+bCk4GPOrAt6KLQdgif2vOLRIYg6k69DTDdEYzJ0Hy5dLEkekNaUUu9y7WFq/lLaKtmSHIxJMzvUpKMMTOdJ6PL5eix1PgEIKNRdtltchGISpUydtSbcv4qPAUSCFjoUQGSupM3IqKiqwWq10dnYOu7+zs5OampoRn/PFL36RD33oQ1xzzTUAzJs3D7/fz3XXXcfnP/95LCMMMHJycsiZpGmcIjXswk1eCk6Pfi26ihgOKrU9tFreScxBAn7w+qCy0qyDU1WVmOMIMYlcIRcFjgLWTVsn9Q7SjJzrM8RgIidD//31Ix3R4qXDaOB9fREAK2zPYPX3Q0EhNDVNWgzesJfqgmpKnCWTdkwhhJhMSb2s4nA4WLJkCc8999zQfYZh8Nxzz7F8+fIRnxMIBA4awFmt5swLpVJzKY2YXAGiuAilXKHjLqOOrcZsQHGM/dn4j4WNgVk4kQjMng3HHitJHJERlFLs9uxmeeNyppVNS3Y4YozkXJ8hYjGzl3SGzsjpxC/LquIgpmy8GF0LwAzrW9RadkIgAC0tkJs7aXH4Ij6ml02XxL8QImMlfcrCTTfdxJVXXslRRx3F0qVL+cY3voHf7+fqq68G4IorrqC+vp57770XgDPOOIP777+fRYsWsWzZMjZv3swXv/hFzjjjjKFBnshuvQQIEKWK/GSHMszrA2vFp1v/Tbll5Ja74xYMgMdjzsKZMzALRwYvIkP0BHoocZbIbJw0Juf6DJDhS6u68JEjiZwJezO2Ao8qIw8vR9leMGv15edP6mwcMBO+jcWNk3pMIYSYTElP5Fx00UV0d3dz22230dHRwcKFC3n66aeHiiLu3Llz2FW5L3zhC2iaxhe+8AX27NlDZWUlZ5xxBl/+8peT9SOIFNNHkCDRlOo80a5PocNowkIsvgWODQNcfebAelYbzJo1ad0ghJgMSinafe2cPuN0mkuakx2OGCc512cAXc/YRE4EHTdhnCk0bkhHPUY17+hHA7Dc/mccWthsuDBnrpnMmSRRPYrVYqW2oHbSjimEEJNNU1k4R9nj8VBcXIzb7aaoqGhiO9u7F2691ZwJUVAQnwDFhPyFbXybV5hLaiwrUgr+GLmULlVPm/U1ltmfj8+OoxHo64OycpgzB2prZRaOyDidvk6iRpQ7Vt9BQ1FDssNJqLiem4T8PuNt1y74/OehunpSP5RPhm783MKzFJJDCZNTjDfTRJWdpyKX0a8qabG8yyrH783ZOIYBq08wO1ZNElfQhS/i46snf5WKvIpJO64QoyHnJhEvmXdZRWQ9s/NE6uQn9xgtdKl6rESZZ3s5PjsNBsHlguYWsxZOXZ0kcUTGMZRBh6+D1c2rMz6JI0TKy+ClVR7ChIjJ0qpxUgr+EV1Hv6rEiZ+l9ucBBT4/NE6Z1CQOgDfipSy3jPLc8kk9rhBCTKbMOxuLrLc3hVqPKwVvDNTGmWV9gzzNP/Gdej3mVOWZs+Cooya1eKAQk6nT10llfiVrpq5JdihCiAxO5LgJE0YnR5ZWjcu/9aVsN2ahoXOC47fkagHzgpMzB5qbJz0eX8TH9HIpdCyEyGyZdzYWWW8vnpSpj7PTmE6vqsFGhHm2DRPcmzKXUuk6LFgA8+eDFP0UGcpQBl2BLk6aehI1BSO3qBZCTKKMTuSYrcctyAf/sdqtt/B67HgAltmeo9qyx3zA64O6eigunvSYDGUwpXjKpB9XCCEmU2p82hUiToIDrcdzU6D1uDkbZyUAs62v4tSCE9kb9PZCjhMWLoQGWWYiMtte715qC2o5qeWkZIcihIAMT+SEkx1CWnIbpfxf9HRAY4b1LWbZ3jIfCIfAbjNbjk/yrJiYEcOiWaTQsRAi42Xe2VhktVTqWLXTmEa/qsROiDm2V8e/I6Wguwdy82DpUkniiIynGzp9wT5OnnqyFKoUIlVkcCKnj6DMxRkjr1HMn6PnE8VJlbaHZbZn9z3o8UJ1DZRPfo0aX8RHgaOA2kJJ5AghMlvyP+0KEUdmIieWEjNy3omZLThnWd8kRxvn1T7DgJ5uKC4x6+EkYVAkxGTb491DXWEdq5tXJzsUIcSgWMz8moF1R7rwS32cMeg3yvlT5AICFFKouVjt+C1WzTAfjEbAoiVlNg6AN+ylNLdULgIIITJe5l1WEVmtjyAGCluSX9rdRi1dqgELOm2218e3EzWQxCkrg2XLJIkjskLMiNEf6mdt61pKc0uTHY4QYpCuJzuChFAouvDjlETOqPQaVfwxcjEBCinRujnV8djwRg5uD1RUmm3qk8AX8dFa2opFk484QojMJmctkVHM1uPJ907sKACmWjeOr1PV4HKq0jJYugyKiuIcoRCpaZd7F00lTTIbR4hUMzgjJ8MEiOInIq3HR6Fdn8Lz0bOJkkO51s7Jjl/h1EL7NojFzPHL1KlJW4KnK52m4qakHFsIISaTJHJERunAhzXJs3G8RjE7jBkAzLGOpzaOgt4eKCo0l1NJEkdkiYgewR/1c/m0yynMKUx2OEKI/cViGbmsyjPQejwfR7JDSVlRZef12HG8qy8BoFrbxUmOX+PQIsM39LihvAzq6pIQpVlfTUOT+jhCiKwgiRyRUfbgTXqh4436EhQW6izbKLX0jPHZCnr7zMLGS46CUllaIrLHTvdOWktbWTllZbJDEUIcKENn5LgJEyYmS6sOodOo5+/RU/Eqczwyw/oWS23PY9MOeD3ounlrnQbW5MxuGix0XFeYnESSEEJMJjlriYwRIkYvgaQmcsIqh036PADmWl8Z+w76+8Fuh8WLobIyvsEJkcJCsRBhPcz66evJs+clOxwhxIEyNpETIoKOXcpGDlEKulUdG2NL2G7MBDTy8LDS/gz11u0jP8njNi8+1ddPZqjDeCNeip3FVObJ+EkIkfkkkSMyhmugY1UZuUmL4X19ATEclGpd1Fp2jO3Jfh8YChYtgFqZFiyyy87+ncwqn8UxDcckOxQhxEgyNpETRhv4L9uFVQ67jam8G1tMj9o3q2Wa9d8cbfvLoTtwGjpEY+ZsHFvyPlr4Ij4WVC/AapF6R0KIzCeJHJExzNbjUXJJTm0NQ1l4L7YYgDm2V8dWSiAcBn8A5s6FKVMSE6AQKSoQDaCjs376enJsOckORwgxksFCthnGwyGSExlMKQiRh18V4lPFdBn1dBiN9KkqGEhoWYjRat1Im/V1yizdh9+hxwMlxdDQkPjgDyNmxGguaU5qDEIIMVkkkSMyRh9BdBT2JHWe2GW0EqAQJ35aLO+O/ol6DPpdMLUVZs3KyGKSQhzOjv4dzK2ay9H1Ryc7FCHEoWTojJxu/ElvkjAWUWXHpSrxqSL8qgi/KiSk8tCxmTdlxcDK/ik3hWXgcSu6shEiF+MQHwGKtD5arRuZYX2LXC1w5IAMAyIRmDvPXBqeJFLoWAiRbSSRIzJGsluPv68vAGC69d9YNWN0T1IG9PRCTS3Mn5+0dp1CJIsn7MFqsXL6jNOxWeSUJETKytBETid+nCncetxQGt2qjr16M+3GFLpVLSou8Spy8ZOveSi3dFFt2UWNZRd5mn9su/F6oKgYGhvjENP4SaFjIUS2kVGzyBid+LEkaY271yhir9EMmImcUevrg5ISs7hxjiwpEdlFKcVO905WNq5kQc2CZIcjhDicaDTZEcSdgaKXADkpOByOKStb9Lm8rR891DFqUB5eCrV+8jUPBZqHXM2PlRhWLYYVHQsGDMzJ0QANY+BxHQs6Ti1IHj4so73odCjKgFAY2mYnfQwzWOi4Or86qXEIIcRkSb0zlxDjtAdP0jpWbdbnARq1lh0UWfpH96SAHyxWcyZOYXLq+giRTK6Qi3x7PqfNOA2LJrPRhEhpkUjGzRr1EiZEjFyStyToQDFlY6O+hI2xJYTIB8BBiHrLNmotO6m17KDQ4k5ylAO8XnP8kgK1/bxhL4tqFkmhYyFE1pBEjsgIEXS6CSRlMGYobajl+AzrW6N7kh4Drw9mz4aamgRGJ0RqUkqx27ObddPWMb1serLDEUIcSQYmctyECaNTgjPZoQDgMip4IXoGblUBQD4e5theYbr139i1FJsRpQwIhmDBTHAm//enG7oUOhZCZBVJ5IiM4BroWJWMwdgeYyoBCskhwBTL5lE8Q0Fvn5nAkeLGIkt1+jspyy1j/fT1aPJvQIjUFw5nYCInRIhY0pdWKWXW2XsldgI6dnLxcZT9/2ixvDfx5U+J4vNBQQE0NSU7EmJGDE3TqC+qT3YoQggxaSSRIzKC2Xo8Rm0SZuR8oM8HYJr1bayafuQneDyQm2suqUpihwchksVQBp2+Ti6ccyENRcltVyuEGKUMnZFjoLAlsWtVVNn5e/RUdhgzAai3bOU4+x9waslt4HBYSkEgAPPmm+OZJPOGvRTlFEmhYyFEVpFEjsgIfQSJYWCf5MGYXxWw25gKwAzrv478hEjEvKq5eAmUlh55eyEy0G7PbuqK6ji59eRkhyKEGK0MTOR4CCf1+Lqy8nz0HNqNJizoLLH9ldnWV1N/om4KzcYBs9BxeW45FXkVyQ5FCCEmTWadkUXWGmw9rk1y16rN+jwUFqq1XRRbXEfYWkG/CxoaoaVlUuITItVE9AjukJvTpp8mg24h0kk0mnGJnH5CSTu2oTT+Gj2ddqMJGxHWOX7OHFsaJHGUMps1NDVDfn6yowHMGTkzymdI0XwhRFaRdzyREbrwT3oSRynYpM8FYIZtFEWOfT5zCnJbW8YNhoUYre3925lRPoPVzauTHYoQYiwycEZOJ34cTH6XI6XgxdhadhgzsBDjJPtvqLLsnfQ4xsXvh7x8aG5OdiRDlFI0laTG7CAhhJgsmXVGFllrD95Jbz3eperxqRLshGmybDr8xoYO/gC0ToOSkkmJT4hU44/40Q2dM2eeSZ49L9nhCCHGIgNn5HThI2eSEzlKwaux1WzW56FhsMr+JLXWnZMaw7gpZSZympvNpVUpIKJHsFltUh9HCJF1MuuMLLJSFJ1OfJOeyNmmtwEwxbIJmxY7/Maufigvg2nTEh+YEClIKcW2/m0srl3M0vqlyQ5HCDFWGZbIiaDjJoxzsscORhvv6EcDsML2NE3W0XS7TBE+r7mcKoWWh3vDXgodhdQXSscqIUR2yZwzssha/QPtQ3MnsWOVoSxs12cAMNX67uE3DocBBbPawOFIfHBCpKC+YB/59nzOmnUWVsvkL2UQQkyAYUAsllGJHA9hwpPcejyg8nk5ehIAC6z/YLrtnUk79oQpAwJBaG1Nmdk4YBY6rsqvosRZkuxQhBBiUmXOGVlkrcHW45M5I2evMYUQ+TjxU2vZcZgtFfT3Q30D1MvVIpGdDGWw27Ob45qOY2b5zGSHI4QYq1jMTOZkUCLHPXARaLKWVikF/4yeTJhcyrROFtj+OSnHjRu3B4qLUmo2DphLdmeUz0BL+SrRQggRX5lzRhZZy0WIMPqkFiwcXFbVbH0fi6YOvaHPB3kDBY5lkCGy1F7vXmoKajht+mky2BYiHWVkIidMGH3SZuRsM2ax05iOhs6x9j9i0YxJOW5c6DpEwjB9BjidyY5miFIKQxk0FjcmOxQhhJh0mXNGFlnLRRCNyWs9HlM2dhijWFY1WBSwZSoUF09KbEKkmnAsjCvk4rTpp1FbWJvscIQQ4zGYyMmgRKyHMACWSRg7BFUeL0fXALDA9hJllu6EHzOu3P1QXg5TpiQ7kmHCepgcW47UxxFCZCVJ5Ii010NwUo+325hKDAcFWj+V2mHahXo9UFgEU6dOXnBCpJht/duYXTGbk6aelOxQhBDjlYEzcrwDiZxEO3BJ1Xzry5Ny3LiJRUE3zNk49smrRTgag4WOpWOVECIbZc4ZWWStvXgntX3o1oFlVS2W9w59cdIwIBgyiwLmSZtlkZ36Q/1YNSvntp1Lrj032eEIIcYrFjMzEhmUyPFMUiKnw2hkhzEDDZ2V9qfTa0kVmHX+qqpSss6fN+KlrrCOwpzCZIcihBCTLnPOyCIrKRQdeCetY1VY5bDbMGfYtFjfO/SGHjeUlkBz86TEJUSqMZTBTvdOjms6joU1C5MdjhBiInQ942bk9BLEPgnD4Lf0FQDMtL5FuaUr4ceLq2AQrFaYNcv8mmL8ET/Ty6cnOwwhhEiKzDkji6zkIYyf6KR1rNqpT8fARonWTal2iDXuegwiUZg2HXJyJiUuIVLNHs8eagpqOHvW2VLgWIh0l4FLq3oJJLxJQofRQIcxBQs682wbEnqsuFPKvCjV1GzOyEkxSinQoLFICh0LIbJT5pyRRVZyESJIdNJm5GwzZgHmbJxDfjbtd0NF6hUFFGKyhGNh+kP9nDHjDGoKapIdjhBiojIskWOgcBFKeMeqt2LLAZhm/Tf5mjehx4o7txuKis3ZOCmYjPdFfOTb85lSLGMtIUR2yowzsshaLoKEiOGchBk5YZVDu2EOGFosh1hWFYuag91p08E2ObOEhEglSik2921mQc0CTmg5IdnhCCHiIcMSOQGihIgldEZOl1FHu9GMhs68dCtwHI1CNAIzZ6ZsnT932E15bjn1RalXu0cIISZDZpyRRdbqI4hCTUr70N3GVBRWSrRuiiz9I2/k8ZizcVKwKKAQk6HL30VhTiEXzbkIp82Z7HCEEPGQYYkcL2Ei6AltlLBvNs47FFo8CTtO/ClwuaCmFpqakh3MIblDbtoq27BZ5KKZECI7ZcYZWWQtFyGYhCQOmPVxAKZYNo+8gaFDTIeWqSlZFFCIRIvoETr9nZw6/VRmVsxMdjhCiHjJuEROhAh6wmbk9Bg17DGmomGk32wcf8Cs7zd7dsqOZZRSGMpgWtm0ZIcihBBJkxlnZJG12vFim4SXcUxZ2WO0ADDFumnkjbxeKC6CurqExyNEKtrSt4U5lXNYP219skMRQsRTxiVywglN5AzOxplq2XjoGbypyNDB54PWVigvT3Y0hxSMBcm150p9HCFEVsuMM7LIWnvxTkrHqg5jCjEc5OGlXOs8eAOlIBQ2Z+M4HAmPR4hU0+3vxmlzctHci8h35Cc7HCFEPMVi5tcULHo7Hj4iAGgJmNHrNYrZZUwDFPNt/4z7/hPK5TKXh09P7Zbe/aF+Spwl0rFKCJHVJJEj0laIGP2EJqVj1U5jYFmVdfPI41i/DwryoaEh4bEIkWoieoR2XzuntJ7CnMo5yQ5HCBFvg4mcDOEdSOQkwmbdfA+ss2yn2OJK2HHiLhgwZ1zNmWsurUph7pCbmeUzybGldpxCCJFIUiFMpK0+ggSJUUZuQo+jFOzUzXXYjSPWx1HmmvLZs1O2u4MQiTLYpWpO5RzOnHkmWoZcsRdC7EfXkx1BXHkIJ2S/SsEWw0zkTLO+k5BjJIRhmM0aZsyE6upkR3NEMSPGjPIZyQ5DiEmj6zrRaDTZYYhJYLfbsY6yPpkkckTachEkSJRcChN6nG5VR4h87ISpsew8eINgEJxOmCJrtUX22evdS1FOER9a8CEKcxL7b1EIkSSxWMYsqwLoJYg9AZPSO4xGfKoEO2GmWA5RTy8V9bugrAza2lL+7zkUC+GwOmgqSd2OWkLEi1KKjo4O+vv7kx2KmEQlJSXU1NQc8eKoJHJE2nIRQkclvNjx4GycBstWrJpx8AZeH0ydCsXFCY1DiFTji/joD/Vz1cKr5OqoEJksw5ZW9RJISKHjwdk4Ldb3sGlp8jsLhUABbbPNi1Ipzh1yU+IskULHIisMJnGqqqrIy8uTWc8ZTilFIBCgq6sLgNra2sNuL4kckbb6CKJITLHC/e00zETOiN2qImGwWaFJrgyJ7KIbOltdWzluynGc0npKssMRQiRSBiVyFAoXwbgncqLKznZ9JgCt6bKsKhYFt9ssblxfn+xoRqU/1M/iusXk2WUpu8hsuq4PJXHKU7iLnIiv3FyzZEhXVxdVVVWHXWYlxY5F2urGjyXBSZx+owyPKseCTr1l28EbeL1QUQkVFQmNQ4hUs9W1laaSJi6bfxl2a+ILjgshkiiDEjl+ooSIkRPna5k7jBnEcFCouajS9sR13wlhGNDbayZw5s1L+SVVgyJ6hLaKtmSHIUTCDdbEyZP6m1ln8O/8SHWRJJEj0tZefAlvPb5rYDZOjWUnDu2ALheGDrphzsZJkwGQEPHQ4evAarFy6dxLqcqvSnY4QohEi8XMSr4ZwEeEMHrcZ+QMdquaZn07DYYECnp6zItQixeDPT2S8RE9gs1io6lYZkGL7CHLqbLPaP/OJZEj0pKOQRc+nIlO5OitAEwZqVuVzwdFhXCE9YtCZBJP2ENvoJdz287lqLqjkh2OEGIyZNCMHC9hIujkxDGR4zWK6DCaAJUGy6oU9PVBQQEsXpJW3TbdITfFzmKpjyOEEEgiR6SpfkIEiCZ0Rk5Y5dCt6gBosG494FEFwRA0NYPDkbAYhEglET3CNtc2Tmw5kdNnnC5XiYTIFhmVyIkQQccex0TO1oEix7WWnRRo3rjtN/4U9PeD1QaLFkFpabIDGhN32M2U4ikUO6W5hBBCSCJHpCUXIULEyCVx04HbjSYUFoq1Xgo0z/AHB1uOp0lxQCEmylAG7/e8z6LaRVw+/3JsFqmVL0TWOMI6/XTiJQwQtxp7SsEWfTZgLqtKXQMzcSwWM4mThrOJg9EgcyrnJDsMIcQRXHXVVWiaxle+8pVh9z/xxBPDLgIqpXjooYdYtmwZBQUFlJSUcNRRR/GNb3yDQCAw2WGnHUnkiLTUR5BwnKdGH2iP0QxA3UhFjn0+cxBUVJSw4wuRKpRSbOrdRENRAx9e9GEKcwqTHZIQYjJFImYCIAP4iBx5ozHwqFI8qgwLMaZYRuhumRKUWdjYkQNHHQ1T0m9pUigWwm61M718erJDESLt6IbipS29/PbNPby0pRfdSHzNM6fTyVe/+lVcLtcht/nQhz7EjTfeyFlnncVf/vIX3nzzTb74xS/y29/+lj/96U8JjzHdySVVkZZcBIHEtR5XCvbqLQAHd6uKRkGzpOVASIjx2Na/jQJHAVcvupqGooZkhyOEmGwZlMjxEAbi9yFmt2HW0qux7MKupeDMJTVQ2Dg/H446CqrSs0B9b6CXyvxKppVNS3YoQqSVp99u584nN9LuDg3dV1vs5PYzZrNubuJm5q1Zs4bNmzdz7733ct999x30+C9+8QseffRRnnjiCc4666yh+5ubmznzzDPxeDwHPUcMlxlnZZF1egkSz4HYgdyqHD9FWIlSY9k9/EGvF8rK0nYwJMRY7HLvwqJZuHrR1SysWZjscIQQyRAOZ0wip5cgtjjO5t1tTAWgwXJgLb0UoOvQ3QXFRbBsWVqPW1whF4tqFuG0OZMdihBp4+m32/noT18flsQB6HCH+OhPX+fpt9sTdmyr1co999zDt771LXbv3n3Q448++igzZ84clsQZpGkaxcVSC+tIMuOsLLLOXrzkJHBC2eCyqmrLbmzafkUelYJYFJqbM2ZQK8ShtHvbCcVCXD7/clY0rkh2OEKIZMmgGTm9BOO2LDuq7HQa5izFlEvkRCPQ0w1V1XDMcigvT3ZE4xYzzHHY7MrZSY5EiPShG4o7n9w44mXvwfvufHJjQpdZnXPOOSxcuJDbb7/9oMc2bdrEzJkzE3bsbJAZZ2WRVRSKjgS3Ht9jHGJZVTAAublQU5OwYwuRCrr93bjDbi6acxEntZyU7HCEEMmUIYkchaKPII44JXL2Gs0YWCnS+iiy9Mdln3ERDIDLBc0tcMwxaV/PzxV0UZZbxsxy+dAnxGht2NZ30Eyc/Smg3R1iw7a+hMbx1a9+lR/96Ee8++67w4+vEl+nJ9Ol/1lZZB0fEbyEE9axKqZsQ1fY6i3bhz/oD0B1DeTlJeTYQqSCDl8HPcEezm07l9NnSptxIbJehiRyAkQJEY1bImdwWVV9qszGMQyzM1UgCLPazJo4zvRfitQb7GVmxUxKc9OrXboQydTlPXQSZzzbjdfxxx/P2rVrufXWW4fdP2PGDN57772EHjvTpf9ZWWQdFyGCxMhN0IycDqMBHTv5eCjWevc9YBiASsuWnUKM1m7PbjxhD5fOvZTzZ5+PRZPThBBZLxrNiESOlwgR9LgszVYKdutmIqcxFRI5wYBZD6cgH5YuhXnzwJq4zp6TRSlFVI+ysHphskMRIq1UFY4uiTva7SbiK1/5Ck8++SQvvfTS0H2XXnopH3zwAb/97W8P2l4phdvtTnhc6S79z8oi67gIEiSasKVVeweWVdVZtzFsIkLAD3n5aV0sUIhDUUqxzbWNqB7l6oVXc+bMMyWJI4QwZUgix0eEMHpcZuT0qSqCFGAjQvWBTREmk6GbtXBCYZg5C447HhobIUNmUnrCHgpzCplZIcuqhBiLpS1l1BY7D9nfV8PsXrW0pSzhscybN4/LLruM//7v/x6678ILL+Siiy7ikksu4Z577uHVV19lx44d/P73v2fNmjX85S9/SXhc6S79z8oi67gIoQBrgl6+++rjbB/+QCBozsbJyUnIcYVIFt3Qeb/3fexWO9ctuY6TW0+W5VRCiH0yJJHjJUwkTomcwWVVtZYdWDV9wvsbl0AAuruhrBxWrIAFC8w6fhmkJ9BDY1Ej9YX1yQ5FiLRitWjcfoZZIPzAEd3g97efMRurZXLGe3fddReGYeyLQdP42c9+xv33388TTzzBqlWrmD9/PnfccQdnnXUWa9eunZS40lniqsUKkSB9BBO2b58qwq3K0TCotezY94AeA4smy6pExglGg2zq20RTSRNXLbiKBTULkh2SECKVGAbEYhmSyIkAYDnkNerR26Mnse34YC0cqwXaZsOsWeBwTH4ckyAQDXBU3VFycUGIcVg3t5bvXr6YO5/cOKzwcU2xk9vPmM26uYn5XPPII48cdF9zczPhcHjYfRaLhRtuuIEbbrghIXFkupQ4K3/729+mubkZp9PJsmXL2LBhw2G37+/v5+Mf/zi1tbXk5OQwY8YM/vCHP0xStCLZOvElbjaO3gxApbaXHG2/Nxt/AAoKoKIiIccVIhl6A71s6tvE0vql3LzyZkniiISSc32aisXMxEEGJHJ8A4mciQopJ93K/ADUYN12hK3jSCnwes1aOCXFsHyFWQsng5M4ObYcWVYlxASsm1vL328+kceuPYZvXryQx649hr/ffGLCkjhi8kxoRk44HCZngstMHn/8cW666SYefPBBli1bxje+8Q3Wrl3L+++/T9UItUgikQgnn3wyVVVV/OpXv6K+vp4dO3ZQUlIyoThE+mjHl7BCx3uNZgDqrNuHPxAKQksL2BPTKUuIyaQbOtv6t2Eog3PbzuX82efjtKV/ZxOROBM938u5Po0NJnIy4PznIYzZdHdi9hotKCyUal3ka96JB3ZEyryg5PNCfgHMmw9Tp2ZER6rD6Q30Ul1QzdTSqckORYi0ZrVoLG8tT3YYIs7G9Gn4j3/8Iz//+c/529/+xq5duzAMg/z8fBYtWsQpp5zC1VdfTV1d3ZgCuP/++7n22mu5+uqrAXjwwQd56qmn+MEPfsAtt9xy0PY/+MEP6Ovr48UXX8Q+MKhobm4e0zFF+oqg00MgIYWOlYIOoxGAWsvOfQ9Eo2bnh+rquB9TiMnWG+hlt3c3U0uncn7b+RzTcIxMWRcHiff5Xs71aSyDZuT0EcQWj/o4k7msKhyCfjfkOs2W4q2t5gzhLNAf6md182oc1syccSSEEBMxqrPyb37zG2bMmMGHP/xhbDYbN998M7/+9a955pln+P73v8+qVat49tlnmTp1KjfccAPd3d2jOngkEuG1115jzZo1+wKyWFizZs2w9mT7+93vfsfy5cv5+Mc/TnV1NXPnzuWee+5B1w9daC4cDuPxeIbdRHoa7FiViBk5/aqCMHnYiFChte97wO+HomIol0y2SF+hWIh3u9/FHXZz5owz+cJxX2B543JJ4ohhEnG+l3N9mhtM5GTAe0UPgQkXOlYK2o0pANQnclmVrkNPD/j80NwMx68yixlnSRLHH/HjtDlZUrck2aEIIURKGtWn4fvuu4//+q//4tRTT8UywhWZCy+8EIA9e/bwrW99i5/+9Kd85jOfOeJ+e3p60HWd6gNmOlRXV/Pee++N+JytW7fy/PPPc9lll/GHP/yBzZs387GPfYxoNMrtt98+4nPuvfde7rzzziPGI1KfixAhYuQS/ynegwOzKsserNpgVXUF4TDMnGnOyhEizYRjYXa4dxDVo7RVtnH+7POZVzVPEjhiRIk438u5Ps3FYmb2Is1n5CgUfQTJmWAix6eKCVKABZ0KrSNO0e1PgdcHwYBZl29Wm9loIcves9t97bSWtTKzXOrjCCHESEaVyDnUFbMD1dfX85WvfGVCAR2JYRhUVVXx0EMPYbVaWbJkCXv27OFrX/vaIQd3t956KzfddNPQ9x6Ph8bGxoTGKRLDRZAIOvYEFDvuGEjk1By4rMpuhxFqOAiRygLRAHs8e4gYEWaUzeDU6aeytH6pTFEXh5Uq53s516eQDFlaFSRGiNiEZ+R0KbMNdrnWiU2LxSO0fQwdevsgxwELFpp1cDKgNtFYGcogEA2wsnElVotcRBNCiJGMeX1KKBTCeYjiau3t7dSOoT1zRUUFVquVzs7OYfd3dnZSU1Mz4nNqa2ux2+1Y95sd0dbWRkdHB5FIBMcIlftzcnImXJRZpAYXIbSB/+LJrI/TAECtZde+B/wBKCyE0tK4Hk+IRDCUQbe/m+5AN3aLnamlU1k7bS3L6peRY5P3QDE28Trfy7k+zel6RiRyvIQJE6OEiRUI7hwYK1RZdscjrH1CIXD3Q2UlzF+Q1V0yewO9lOeWs7h2cbJDEUKIlDXms/LixYt58803D7r/f//3f5k/f/6Y9uVwOFiyZAnPPffc0H2GYfDcc8+xfPnyEZ+zcuVKNm/ejGEYQ/d98MEH1NbWjjiwE5mlB39C9tunqoiQi50w5ftPlY6Eob4+7QewInNF9Sidvk42dm/kna530JXOumnruPW4W7nzhDs5vul4SeKIcYnX+V7O9WkuQ2bkeIkQQSdngjX2ugyzyHeVZU88wjJ53ODzwfQZsPLYrE7iAHT6O1lcu5jK/MpkhyKEEClrzGfl1atXc8wxx/DVr34VAL/fz1VXXcWHPvQh/uM//mPMAdx00038z//8Dz/60Y949913+ehHP4rf7x/qbHHFFVdw6623Dm3/0Y9+lL6+Pj796U/zwQcf8NRTT3HPPffw8Y9/fMzHFulnL76EdKwaXFZVbdmNRRtoTRqLgtVmXh0TIgUopQjFQnT5u9jSt4V/d/2bD/o+QKE4bspxfGLpJ7j3pHu5ZvE1zK+ej80S/38rInvE83wv5/o0liGJHN9AImciS6vCKod+ZY4JquOVyHH3m7/fRYvMW5bPKgvFQtgsNpY1LEt2KEIIkdLGPMr/zne+w2mnncY111zD73//e9rb2ykoKGDDhg3MnTt3zAFcdNFFdHd3c9ttt9HR0cHChQt5+umnh4oi7ty5c1jBxcbGRp555hk+85nPMH/+fOrr6/n0pz/NzTffPOZji/RioOjAl5COVe0DbcdrDlpWVQBlZXE/nhCHYiiDiB4hokcIx8IEY0GC0SARPYKmadgtdopyiphZMZOZ5TNpLWtletl0CnMKkx26yDDxPN/LuT6NZVAiRwGWCSzN7jLM+jhFWh9OLTjxoPr7za8LF5mdqQQdvg4aihqYUzkn2aEIIURKG9cn4lNPPZVzzz2X7373u9hsNp588slxJXEGfeITn+ATn/jEiI+98MILB923fPly/vnPf477eCI9eQnjJxL3GTmG0ugcSuTsV+g4HDILDUq3qqymlMJQxiFvCvPxwe0Uatifh7bbbz+60s2boaNQWLCgUEPHdFgd5FhzyLHlUJ5bTmNNI/WF9VTkVVBTUEN9UT0FjuxoQSuSK57neznXp6kMSeT4iUx4H4OJnIkvq1JmEkfTzCROU9OEY8sESik8YQ9nzzpblgQLkQFeeukljj32WNatW8dTTz01dP/27dtpaWk55HOOOeaYyQoxrY35E/GWLVu49NJL6ejo4JlnnuH//u//OPPMM/n0pz/Nl7/8ZexZWF1fTA4XIYLEKCU3rvvtVdVEycFBiDKty7xTj4HFKt2qMoBSipgRI6JHiBpRonp06GvMiA3dDGWgaWYh7cGkikKhoWHRLCPeBh/TNO2g7+0WOxbNgs1iw26x47A6hm45thzy7Hk4bU6cNie5tlxy7bk4bU7y7HkU5RRR6CikMKeQXFuutAoXSSHnewFkTCLHF8dETrU2wULHbo+ZxFm0GKZMmXBcmaI/1E9xTjFLapckOxQhMouhw44XwdcJBdXQtML8nJNgDz/8MJ/85Cd5+OGH2bt3L3V1dcMef/bZZ5kzZ/jsu/Ly8oTHlSnGnMhZuHAhp512Gs888wwlJSWcfPLJrF+/niuuuII///nPvPHGG4mIUwj6CBIiFvcZOfvq4+zaVx/HH4CCfJA3k5QX1aMEY0FCsRDhWNj8qoeHHlco7Bb7UDLFbrXjtDmpyKugwFFAYU4hhY5C8ux5Q7NgHFbHsO3tFruZkLGaX20WG1bNitViPejrYPJmMNkjRLqS870AzEQOmImHNNZPaELLqnRloUeZXdYmNCMnGDBr8C1eIkmcA7T72llWv4yGooZkhyJE5tj4O3j6ZvDs3XdfUR2s+yrMPjNhh/X5fDz++OO8+uqrdHR08MgjjxxUX6+8vPyQ3SvFkY2rRs6HPvShYfetWLGCN954gxtvvDFecQlxEBfmevSJDMRG0j6QyKndf1lVKGROdbZJsdhUoJQiGAsSiAbwR/z4o34MZXazsWrWodksZbllVOVXUZ1fTbGzeChJU+AoIN+RT549j3x7Pk6bU2a5CHEEcr4XwL5ETprrI4h97D0+hvSqanTs5BCgSHONbyfRKHi9MKtNauIcIBgNoqGxqnmVnJ+FiJeNv4NfXAH7Ld8HwNNu3n/hjxOWzPnFL37BrFmzmDlzJpdffjk33ngjt956q/z7jqMxf0o9cFA3qLCwkIcffnjCAQlxKC5Ccd+noSxDU6WHCh3rOlg0WVaVJEopAtEAnrAHX8RHSA+Bglx7Lnn2PKoLqmkuaaa2oJay3LKhW4mzRBI0QsSRnO8FYJ4TM4CL0IQ6VnUZ5iyRKsue8U1OMgzo64XGKTB7dtrPcIq3Pd49TCubxoLqBckORYjMYOjmTJwDkzgwcJ8GT98Cs05LyDKrhx9+mMsvvxyAdevW4Xa7+b//+z9Wr149tM2KFSuGNToAcyaPGB2ZbiDSRju+CV1NG0mPqiGGgxwClGrd5p2BAOTlybKqSaKUwhvx0h/qxxv2YmCQZzPrxCyoWUBraSu1hbXUFNRQmVdJgaNAkjVCCDFZYrG0TzroGHgITzCRM1AfZ1zLqhT09kBFBSxcKLN9DxAzYgRjQU5oOQG7VWpvCREXO14cvpzqIAo8e8ztWo6L66Hff/99NmzYwG9+8xsAbDYbF110EQ8//PCwRM7jjz9OW1tbXI+dTeRMItLGXjwJqI+zr+340Dg1GIRp08DhiOuxxD7hWJjeYC/9oX50Q6cwp5CqvCqOm3IcU0un0lhsdmmSrhVCCJFkGbC0yk+UCDo540zkKAVdhlmkc1z1cbxeyM01ixvn5Y0rhkzW4eugtqCWpfVLkx2KEJnD1xnf7cbg4YcfJhaLDSturJQiJyeHBx54YOi+xsZGpk2bFvfjZwtJ5Ii0ECJGPyFyie+Vms6BqdLVloEOFAN1V2RZVfwFo0G6A924Q25sFhtV+VWc2HIisypm0VraSn1RvRQHFkKIVJMRiZwIEXQKGd8FGo8qJUQ+FmKUa2P80GMYEAjCggVQVjau42cypRS9wV4unXcpRTlFyQ5HiMxRUB3f7UYpFovx4x//mK9//euccsopwx47++yzeeyxx1i3bl1cj5mtJJEj0oKLIEFilMWx9bihtKErbEOJnGAInE4ZbMVJVI/S5e+iL9hHji2H2oJaTm49mTmVc5hRPgOnzZnsEIUQQhxOBiRyfAOJnPEurRpcVlWhdWDVxlgzyN0PpSXQ0jKuY2e6nkAPZc4yVjauTHYoQmSWphVmdypPOyPXydHMx5tWxPWwv//973G5XHzkIx+huLh42GPnnXceDz/88FAip7e3l46OjmHblJSU4HTK54PRkESOSAsuQgSJkkth3PbZryqI4sROeF99nGAQqqvNKdBiXJRSuMNu2r3tGMqguqCas2edzeLaxUwvn47DKkvWhBAibcRi5tqiNDaYyLGPN5GjBuvj7B7bE2Mx8zZjJuTIUuGRdPo6Obn1ZOoK6468sRBi9CxWs8X4L64ANIYncwbqSaz7StwLHT/88MOsWbPmoCQOmImc++67D4/HA8CaNWsO2uaxxx7j4osvjmtMmWpc6xhOO+002tvbD/qzEInSRxAdhS2OxY4Hr7BVWvZi0RSgzAFXdXXaF3ZMBkMZtHvb+XfXv3EFXSytX8pNy2/iK2u+wocWfIg5VXMkiSNEmpHzvciERI6fKAAWxndu7zFqAKiwjPH13++CyipobBzXcTOdN+zFaXdyXFN8C60KIQbMPtNsMV5UO/z+orqEtR5/8skneeqpp0Z8bOnSpSilmD9/PkqpEW+SxBm9cc3I+etf/0owGDzoz0IkigvzNaaNcxA2kn31cQYKF0ajYLdLt6ox0g2dPd499If6qcqv4pxZ57CicQVTS6dKdykh0pyc70UmdK3yERn3c3VlpV+Z44JyS9fonxgOg2aBmTPBGv/Wvplgr3cvc6vmMqtiVrJDESJzzT7TbDG+40WzsHFBtbmcKgEtx8XkkqVVIi10E4hjCse8uDiYyKnSBqZKB4JQkA8lJXE8UubSDZ293r24Qi7qC+s5e+bZrJyykvI8SYQJIUTGiIw/CZIqJpLI6VflKKw4CJKPZ5TPUmZtnClNUFMz7mNnsqgeJWbEOL7peGl0IESiWaxxbzEukk8SOSIt7MUb145VfooIUIiGToVloMhWOGwWI7TIgOJwlFJ0+jvp9HVSX1TPWTPP4vjm4ylxliQ7NCGEEPEWjab9ebGf0LiXVfUZZhfLMkvX6CcmBUPgyIEZM9J+NlOidPo7qS2sZVHtomSHIoQQaUkSOSLl6Rh04sMZx5dr50B9nHKtE7sWNduDakBFRdyOkYk8YQ87+ndQllvGpfMvZU3LGkpzS5MdlhBCiESJRNI+kdNHcNwdq3qV2Zq3XBvDsiqfz6yLUyrnx5EopegL9nFK6ykUOAqSHY4QQqQlSeSIlNdPiABRCohfodyuA+vjBIPgzJW244cQ1aNsdW1F0zRWN6/mrFlnMaV4SrLDEkIIkWjhcNonclwTSOTsPyNnVGIx88LQlCkyG+cQ+kP9FOUUsbR+abJDEUKItCWJHJHyXIQIEaOS/Ljtc3BGTtX+iZy6OnA643aMTNET6GGPZw+zKmZx/uzzWVS7SNazCyFEtkjzpVU6Bh7C2MfR9VIp6FMDiRytc3RP8nmhuMTsgClG1O5rZ1n9MpqKm5IdihBCpC1J5IiU10eQMDo547yadqCwyqFfVQKDiRwFui6DrgNE9ShbXFtwWB2c23YuZ886m8KcwmSHJYQQYjKl+dIqP1Ei6ONanu1VJcRwYCVKsdZ35CcoBeEIzGmWTlWHEI6FATh2yrHS2VIIISZgXImcpqYm7Hb7QX8WIhH64tx6vGtgNk6R1kuuFjAHqQ6HtB3fjzvkZnv/dmZVzOLiuRczv3q+DLiEyEJyvhfpPiPHR4QIOkXkjPm5g/VxSrUeLJo68hMCAcjLNWf4ihG1+9ppLG5kQc2CZIcihBBpbVyJnLfffnvEPwuRCL0DiZx46TywPk4gCIWFUFQU1+OkI6UUuz278UV8rJu2jkvmXUJRjvxehMhWcr4XRCJpXevFR4Qoxrhq5Iy5Po7fD9OmQX78loJnEqUU7pCbc2adg9MmS9mFEGIi0vcSi8gae/HEbVkV7JuRU6XtNu+IhKGmJq2vOMZDzIjxbve7WDQL1y65lmuXXCtJHCGEyHZpPiPHPzAjxz6OccTgjJxR1ceJRszlVA0NYz5OtugL9lGWW8bR9UcnOxQhRAKtXr2aG2+8Ma77vOOOO1i4cOGE9qFpGk888URc4kkF6XtmFllBoWjHRy7xmc4fU1Z6VA0wMCPHMMwrjVnerSoYDfJO9zu0lrXy2RWfZc3UNVLQWAghRNoncnxEALCMY3n24Iyc8tHMyPF4zSXaFRVjPk626PJ3MadqDnWFsvRMCDE2n/3sZ3nuuedGte2hkj7t7e2ceuqpcY4seaTYsUhpXiL4iJAbp5dqr6rBwIYTP4VaPwRDZqeq0tK47D8duUNudrh3sKJhBR9Z/BHKcrM7qSWEEGKAUmY77TRP5IxnYVhA5RMiHw2DUq378BsrA/QYNDWl9e8qkXRDR1c6R9UdlexQhBBpqKCggIKCggnto6amJk7RmKLRaFJrB8rZRqS0PoIEiY6r28RIBpdVVVv2mEv+Q0GzTWiWth3v9nezy7OLta1r+fjSj0sSRwghxD6xmDlzNY2TEz4iGIyiUPEBBmfjFGl92LTY4TcOhiA3V7pfHkZvsJey3DLmVs1NdihCiEnmcrm44oorKC0tJS8vj1NPPZVNmzYN2+Z//ud/aGxsJC8vj3POOYf777+fkpKSoccPnGXzwgsvsHTpUvLz8ykpKWHlypXs2LGDRx55hDvvvJO33noLTdPQNI1HHnkEOHhp1e7du7nkkksoKysjPz+fo446ipdffnnEn2H79u1omsbjjz/OqlWrcDqdPProowB8//vfp62tDafTyaxZs/jOd74z7LkvvvgiCxcuxOl0ctRRR/HEE0+gaRpvvvnmuH+nIDNyRIpzESRELG5Lq7oNczpvpWWveUcsZg680riQ43jt9uzGH/Fz4ZwLObftXGwWeTsQQgixH103Ezlp3Eq7nxDWcVy3HKyPU66NYllVIGB2qsrLG/NxskVPoIfjm46XC0ZCTIBSimBUT8qxc+3WcXewveqqq9i0aRO/+93vKCoq4uabb2b9+vVs3LgRu93OP/7xD2644Qa++tWvcuaZZ/Lss8/yxS9+8ZD7i8VinH322Vx77bU89thjRCIRNmzYgKZpXHTRRbz99ts8/fTTPPvsswAUFxcftA+fz8eqVauor6/nd7/7HTU1Nbz++usYhnHYn+WWW27h61//OosWLRpK5tx222088MADLFq0iDfeeINrr72W/Px8rrzySjweD2eccQbr16/nZz/7GTt27Ihb/aBRf3J79tlnWbNmzSEfNwyDe+65hy984QtxCUwIMGfkKMa3tv1ASkHX/okcPQYWa1Yuq9rl3kVEj3D1oqs5eerJ0lpcCDFEzvdiSCxmnjzT+BzRRyjBHasUGLrZNEGMKGbEMJTBktolyQ5FiLQWjOrMvu2ZpBx7411ryXOM/aLvYALnH//4BytWrADg0UcfpbGxkSeeeIILLriAb33rW5x66ql89rOfBWDGjBm8+OKL/P73vx9xnx6PB7fbzemnn05raysAbW1tQ48XFBRgs9kOu5TqZz/7Gd3d3bzyyiuUDdRKnTZt2hF/nhtvvJFzzz136Pvbb7+dr3/960P3tbS0sHHjRr73ve9x5ZVX8rOf/QxN0/if//kfnE4ns2fPZs+ePVx77bVHPNaRjPoSxfr16/nEJz5BIBA46LG3336bo48+mu9+97sTDkiI/bkIxW1fPlU0sN5dp1zrHJgK7YT9pu1lg53unUSNqCRxhBAjkvO9GJIBS6tcBMeXyFEDiZwjdawKh8GRYxY6FiPqCfRQkVfBnKo5yQ5FCDHJ3n33XWw2G8uWLRu6r7y8nJkzZ/Luu+8C8P7777N06dJhzzvw+/2VlZVx1VVXsXbtWs444wy++c1v0t7ePqa43nzzTRYtWjSUxBmto47aV+fL7/ezZcsWPvKRjwzV8CkoKOBLX/oSW7ZsAcyfbf78+Tj3K+NxuJ9tLEadVvvb3/7GVVddxYIFC3jkkUdYuXLl0FW5u+++m/POO29o+pIQ8bIXL/Y4lXLqVuZsnHKty1zvHgpCQyM4HHHZfzrY3r8dpRQfXvRhVjWtkiSOEOIgcr4XQwaXVqVpIieGgZfwmBM5EeXAq8zZukeckRMIQkkxFBWNN8yM1xvo5eTWkynKkd+REBORa7ey8a61STt2KvnhD3/Ipz71KZ5++mkef/xxvvCFL/DnP/+ZY445ZlTPz83NHddx8/Pzh/7s8/kAs77P/okqAOskLEkedSJn2bJlvPHGG9xyyy2ccMIJXHfddfzzn/9k165dPPbYY8OmGAkRL3vxJqg+jgLdgMrKuOw7Hex07wTgI4s/wvFNxyc5GiFEqpLzvRiS5jNy/ESIoI95HDE4GycfD07tCDODIxGorUvr5WeJFNWjaJrGoppFyQ5FiLSnadq4ljclU1tbG7FYjJdffnloaVVvby/vv/8+s2fPBmDmzJm88sorw5534PcjWbRoEYsWLeLWW29l+fLl/OxnP+OYY47B4XCg64evJTR//ny+//3v09fXN+ZZOYOqq6upq6tj69atXHbZZSNuM3PmTH76058SDofJyckBRvezjcaYXglOp5P/+q//oquri+985zvk5+fz6quvMnPmzLgEI8T+gkTpJxS3jlXDEjnRKNjtWVMfp93bTkSP8JFFksQRQhyZnO8FkPYzcnwDiZwicsb0vH31cY6wrCoaBZsNKirGG2LG6w50U5VfxezK2ckORQiRBNOnT+ess87i2muv5Xvf+x6FhYXccsst1NfXc9ZZZwHwyU9+kuOPP57777+fM844g+eff54//vGPh1w5sG3bNh566CHOPPNM6urqeP/999m0aRNXXHEFAM3NzWzbto0333yThoYGCgsLh5Iogy655BLuuecezj77bO69915qa2t54403qKurY/ny5aP++e68804+9alPUVxczLp16wiHw7z66qu4XC5uuukmLr30Uj7/+c9z3XXXccstt7Bz507+8z//E2DCKyPGdGbesmULxx9/PM8//zwPPvggc+fOZfXq1fz2t7+dUBBCjGSw9XhuHBI5MWWjd+AKW6W216yPk5cLI1QxzzS9gV7cYTcXzbmI1c2rkx2OECINyPleABkwIydKBH3MS6tcypytW6Z1H37DYAAK8rPmotB49AX7WFK7hHxH/pE3FkJkpB/+8IcsWbKE008/neXLl6OU4g9/+AN2uzlbcuXKlTz44IPcf//9LFiwgKeffprPfOYzw+rK7C8vL4/33nuP8847jxkzZnDdddfx8Y9/nOuvvx6A8847j3Xr1nHCCSdQWVnJY489dtA+HA4Hf/rTn6iqqmL9+vXMmzePr3zlK2NeEnXNNdfw/e9/nx/+8IfMmzePVatW8cgjj9DS0gJAUVERTz75JG+++SYLFy7k85//PLfddhvAIX++0Rr1J+QHHniAW265hbVr1/LrX/+ayspKrrnmGr72ta9x8cUXc/755/Otb31rWL93ISbCTOTEp/V4r6pGYcWJnwLNA+EQNDSkdUvV0XCH3HT4Oji37VxOm3Ga1MQRQhyRnO/FkDTvWuUjQhRjzImcfsMsXFxi6Tn8hqEwNLdk/FhivMKxMFbNysKahckORQgxiV544YVh35eWlvLjH//4sM+59tprh3Vyuvbaa4d1kbrjjju44447AHNJ029+85tD7isnJ4df/epXB92vlBr2fVNT04jbjaS5ufmg5w+69NJLufTSSw/53BUrVvDWW28Nff/oo49it9uZMmXKqI59KKO+xHLbbbfxve99j//93/+lcqCuiMVi4eabb+bVV1/l3XffZc4cqUYv4sdFCAOFLQ7FjgeXVVVZ9qKhQJHxHSYC0QA73Ds4pfUULphzARYtPa+oCiEml5zvxZAMWFqlARqjT0QpBf1qIJGj9R56Q0MHDVlWdRi9wV4q8yuZVTEr2aEIIVLcf/7nf/LWW2+xefNmvvWtb/GjH/2IK6+8MtlhxcWPf/xj/v73v7Nt2zaeeOIJbr75Zi688MJxF1weNOoZOe+88w61tbUjPjZnzhxefvll7rnnngkFI8T++gjGbV/D6uNEwmanqgyeCh0zYmzq3cTKKSv50IIPYbOkV2E0IUTyyPleDEnzpVU+Iox8/fTQguQTxYmGQZHmOvSGgSDk5mX8RaGJcAVdnNJ6Crn2iX1YEUJkvg0bNnDffffh9XqZOnUq//3f/80111yT7LDioqOjg9tuu42Ojg5qa2u54IIL+PKXvzzh/Y76092hBnWDrFYrX/ziFycckBCDuvBjGcNVtENRCrr2T+T4QlBQAIWFE953KlJK8UHvB8yqmMXVC6/GaZvY+kshRHaR870YkuaJHD+RMT9ncFlVodaPVTtM15NgEJqaIGdshZSzhW7oKBRzqmT2nhDiyH7xi18kO4SE+dznPsfnPve5uO83Pc/MIivswROXQsd+CglSgIZBudZptgqtqkrbNf9Hsr1/O2W5ZVy96GrK8+RKoRBCiHFK86VVLkJYx3hByD2aZVUoUIY5lhAjcoVclDpLmVkune6EECIRRnVmtlgsWK3WMd/uuuuuRMcvMlQUnW4CcSl0PLisqkzrwq5FAAUZWqSz299N1Ihy2bzLmFE+I9nhCCHSjJzvxTBpXuy4jyD2sRY6VmbNm+LDJXIikYxfoj1RvYFeppdPlwtKQgiRIKOa7rBt27Zx7Vw6WojxchEiRIxiJj5leXh9nAjYHRnZdjwQDdDh6+CC2Rdw7JRjkx2OECINyfleDBOLpW0SB8BFcAIdqw6TyAmFITc3Y5doT5RSirAelm5VQgiRQKNK5DQ1NSU6DiGGMVuPR6mhYML72lcfpx38YcjLvMGXoQw2921mecNyzp19rrQZF0KMi5zvxTD6YWrEpLgYBj4iY07kuFUZACXaYVqPh8NQX5+2S84SzRvxUuAokJnBQgiRQKM+A3V1dR328VgsxoYNGyYckBBgJnKiGNgnWMZJV1b6lLmGvVLbaw6+yisybvC1zbWNhqIGLp9/OQ6rI9nhCCHSmJzvxZBYLNkRjJuPCBH0MSVyQiqXEPmAoljrO8RWA/VxZFnVIfUGemkoamBK8ZRkhyKEEBlr1J9ma2trhw3u5s2bx65du4a+7+3tZfny5fGNTmQt10DrcW2CXat6VRUGNnIIUKi5zMFXWVk8QkwZfcE+dKVz8dyLqS08fLcZIYQ4EjnfiyFpvLTKP45ETv9AoeMCzY1NO0QSKxoDmz1ja+3FgzfiZXHtYixaZl00E0KIVDLqlkBKqWHfb9++nWg0eththBivbgJxaDwOPYaZ2Ki07EXTY2C1ZVR9nIgeYbdnN2fMOIPlDfLBSggxcXK+F0PSeGnVeGbkDNXHOVyh43AInE4oKppoiBkpFAvhsDqkW5UQieR2QyAwecfLy8uoz0+JdMcdd/DEE0/w5ptvJvxYE+/tvB+pyyHiZQ8enHF4ee5L5HSYxQmdORkz+FJKsalvE3Or5nLBnAvk358QYtLI+02WSPOlVWNdoj3YevywHatCYWhoAFtch9AZozfQS2V+JdPKpiU7FCEyk9sNd98NPYep4xVvFRXwxS+OOZnT0dHBl7/8ZZ566in27NlDVVUVCxcu5MYbb+Skk06iubmZG2+8kRtvvDExcWc4OQuJlKNj0IEvPq3HlZnIqdDazatotXVgn/h+U0Gnv5MiRxGXz7+cAsfEi0ILIYQQw6R5IgfGtkR7cGnVYWfkGDqUS0vtQ+kP9bO2YS259txkhyJEZgoEzCRObq45U2ayjhcIjCmRs337dlauXElJSQlf+9rXmDdvHtFolGeeeYaPf/zjvPfeewkM+vCi0Sj2DPg8OOrLFJqm4fV68Xg8uN1uNE3D5/Ph8XiGbkLEg4sQAaLkTjDPGFJOvMosRlhh6TCniFdUxCPEpIvoEbr93ayfsV66Qggh4krO92KIrkOaLqPzEz3yRgdwDyytKj5U63E9BharLDE4BP3/t3ffcXJX9f7HX9+Z2dnZ3lt67z2QAtKDQRBBhR9yRYp6rwqIiCiiaASFgIoXRASUdlG5IFexoRHpvZgQpASE9La9z+zU7/n9Mdklk22zs2Vmdt9PH3tz91vPfNnsOfl8z/l87OhSvPll85PcEpExIDs7WoV3uL8SDBZdeOGFWJbFK6+8wic/+UlmzZrF/Pnzueyyy3jppZfiusa9995LYWEhf//735k7dy65ubmcdNJJ7N+/v+sY27a55pprmDBhApmZmSxZsoQNGzZ07d+xYweWZfHggw9yzDHH4PF4+M1vfsP555/P6aefznXXXUdFRQWFhYVcc801hMNhvv71r1NcXMyECRO45557Ytp0xRVXMGvWLLKzs5k2bRrf+c53ui0/HykDypEza9asmO+XLl0a872mWstQiJYeD1NOzqCu07msKt9qJNN4wXKMmsHX1satzC+fz0dmfCTZTRGRUUb9vXQZBTNy4hU0bnzkAX3MyOlcoj1KxhJDrcnfRKGnkJklM5PdFBFJosbGRjZs2MC1115LTk73f88V9pIs/vzzz2fHjh089dRTXdt8Ph8//vGP+dWvfoXD4eCcc87h8ssv5ze/+Q0AN998MzfeeCN33HEHS5cu5e677+ZjH/sYb731FjNnfvC76Jvf/CY33ngjS5cuxePx8NRTT/HEE08wYcIEnnnmGZ5//nk+97nP8cILL3D00Ufz8ssv8+CDD/KFL3yBE088kQkTJgCQl5fHvffey7hx43jjjTf4z//8T/Ly8vjGN74xdA8wTnEHcp588snhbIdIlwZ8BAgPKEFhT+oPXlY1igZfdd46PC4PZ80/i+yMEZhSKSJjivp76RIOp+2MnCY6cA4gP07nsqps2nBbvQSBAn4oK4fMzKFo4qjT2NHIoopFlGRp6ZnIWPb+++9jjGHOnDkDOq+qqgrbtmO2hUIhbr/9dqZPnw7AxRdfzDXXXNO1/8c//jFXXHEFn/rUpwC44YYbePLJJ7npppu49dZbu4679NJL+cQnPhFz7eLiYn7605/icDiYPXs2P/zhD/H5fHzrW98C4Morr+T666/nueee67r+VVdd1XX+lClTuPzyy3nggQdSO5BzzDHH9HtMY2PjoBojAtEZOdaB/w1GvV0JQJnjQH6c4uK0H3yFIiH2te3jjHlnMK9sXrKbIyKjkPp76RIKgSM9S0g30IF7IImOOytW9basCqKBrbKywTZtVDLG4A/7WVyxWDP2RMa4RCtbrl+/vtu27OzsriAORIM9tbW1ALS2trJv3z6OPPLImHOOPPJIXn/99Zhthx12WLdrz58/H8dBfVxFRQULFizo+t7pdFJSUtJ1P4AHH3yQn/70p2zdupX29nbC4TD5SSqkMyS986OPPsr/+3//j/Hjxw/F5WSMq8E76NLjxkDdgaVVpY790cFoaRmk+eBiW9M2ZpXM4tTZp2qgJCIjTv39GBMMpm2/2UTHwEqPm2gOvV4rVtl29FmMgpm9w8EX8pGdka1lVSLCzJkzsSxrSBIaH5qU2LKshAJFPS3x6unaPW3rnCX04osv8ulPf5qTTz6Zv/zlL7z22mt8+9vfJhgc2FLeoZJwIGfnzp2sW7eOKVOmcOaZZ+JwOLjvvvuGsm0yRu2hddAVq9pNAQGycRCmmFrAgl7WY6aLFn806eiZ888kP3N0lFAXkdSn/n4MS9MZOUEieAkNMJDTT8WqgB8yPQrk9KKxo5GK3AomF0xOdlNEJMmKi4tZu3Ytt956K16vt9v+5ubmIblPfn4+48aN4/nnn4/Z/vzzzzNv3tCvXHjhhReYPHky3/72tznssMOYOXMmO3fuHPL7xGtAZYGCwSC///3vufPOO3n++edZs2YNe/bs4bXXXmPhwoXD1UYZQ8LY1OIddMWqzrLjxVYtzpAPMt1pPfgyxrCzZSfHTTmOZVXLkt0cERnl1N8LAIFAWgZyvAQJEiFnAC+FWkwxAIWO+p4P8AegsCBa8le6aQm0cNyU48hwpn9JXxEZvFtvvZUjjzySFStWcM0117Bo0SLC4TD/+Mc/uO2229iyZUu3c6688kr27t07oJdFX//611m3bh3Tp09nyZIl3HPPPWzevLkrGfJQmjlzJrt27eKBBx7g8MMP55FHHuHhhx8e8vvEK+5/LX/5y1/mf//3f5k5cybnnHMODz74ICUlJWRkZOB0Di4prUinRjrwESKfweWyqe9aVlUNvkC0dF4PU+rSRXV7NSVZJZw25zQcVvoNqkUkfai/ly5pOiOn/UAgpwhPXMeHTAbtphDoY2lVKBRNdJymS82GU9gOY2Exu3R2spsiMnb4fCl9n2nTprFp0yauvfZavva1r7F//37KyspYvnw5t912W4/n7N+/n127dg3oPpdccgktLS187Wtfo7a2lnnz5vGnP/0ppmLVUPnYxz7GV7/6VS6++GICgQCnnHIK3/nOd/je97435PeKh2XiXGTmcrm44oor+OY3v0leXl7X9oyMDF5//fVhmb40XFpbWykoKKClpWXwyYn27YMrr4wmv8vNHZoGjmHvUM/VPMVUigZVteqvgf+g1oznqIxHmN74NEyfDsvScyZL2A7zdt3b/MfC/+CMeWckuzkiMkyGtG8ahNHS36fK80xr110H77wD06YluyUDsoU6ruZppsU5lqi3K/hL8Fw8ePmU5+c9HGGgtg5WrYJJk4a+wWmuzltHyA5xw5obKMoqSnZzRFJavH2T3+9n+/btTJ06FY/noKB0Swt8//tQ38vsweFQWgrf+U5ar25IJ73+tz9E3DNyfvWrX3H33XdTVVXFKaecwmc+8xk+8pGPDEljRTo14COETcYg8nDbxkGDKQeg1NoHxk7r/Dg7mncwvWg6a6evTXZTRGQMUH8vXdI02bGXEOEBjCVaDuTH6X02ThhcLjgosCkfaOxoZPm45QriiIyEgoJoUGWkZuRAdGWDgjgpJ+5Aztlnn83ZZ5/N9u3buffee7nooovw+XzYts3bb7+dNm/oJLU10gEwqNLjTaaUCBm48ZMfqQenC9L0baw36CUUCXHanNPIy9QAUkSGn/p76ZLGS6sg/rFEix3Nj1PgaOz5gGAgmmtPM6+7McYQiARYVLEo2U0RGTsKChRYkYFPe5g6dSpXX301O3bs4Ne//jWf/OQnOeecc5gwYQKXXHLJcLRRxpBo6fHBvf2rs8cB0bLjVjAAmZlpG8jZ0byDZVXLWDVhVbKbIiJjjPp7SedAzkBGEq0mOpOkwGrq+YBgEPLyIUOJfA/lDXnJychhZrHKjouIjKSESwNZlsXatWtZu3YtjY2N3Hfffdxzzz1D2TYZg/bQSvYgK1bVm0oASq390YobpaXgdg9F80ZUU0cTWRlZfGz2x3A5BvdMREQSpf5+DAsG0zKQ00YAQ1wpIIEPAjn5Vi8zckIhKC4eiqaNOg2+BipzK5lUoNxBIiIjaUh65+LiYi699FJef/31obicjFEhItThJWsA5UJ70lmxqsyxH8LhaCAnzRhj2NO6h9UTVjOndE6ymyMiAqi/H1OMifahaRjIacaPK86CCcYcHMjpaUbOgYCQllX1qDXQyrJxy3A6VNFORGQkxdU7X3/99XR0dMR1wZdffplHHnlkUI2Ssamz9HjWIGbkBI2b5gNJC0ut/dGNabisqt5XT4GngJNmnISVhokmRSQ9qb+XLpEI2HZaBnIa6Ii78qWPXMK4sbDJtVq6HxCORHPtKdFxN6FICIfDwewSlR0XERlpcfXOb7/9NpMmTeLCCy/kb3/7G3V1dV37wuEw//rXv/j5z3/OEUccwVlnnRVTrlQkXg100EF4UDNyGuxKwCKHFrJCzdElVWkWyDHGsK9tH0dNOoqpRVOT3RwRGUPU30uXcDg6XSXNXiYYDM344w7ktNrR2Ti5VgtOy+5+gBId96qho4GSrBJmlcxKdlNERMacuAI59913H4899hihUIj/+I//oLKyErfbTV5eHpmZmSxdupS7776bc889l3feeYejjz56QI249dZbmTJlCh6Ph5UrV/LKK6/Edd4DDzyAZVmcfvrpA7qfpKZGOgZULrQnnflxyhwH8uNkedJu8FXdXk1ZThlrZ6jcuIiMrOHs79XXp5lwOC1n5ASJ4CUY91ii72VVRPME5eSmZa694dbY0ciiikXkutNrnCUiMhrEvYZl8eLF/PKXv+SOO+7g9ddfZ9euXXR0dFBaWsqSJUsoTTAPyYMPPshll13G7bffzsqVK7nppptYu3Yt7777LuXl5b2et2PHDi6//HKOOuqohO4rqacBHzC40uOd+XFKHdXRt2hVVWk1CLWNTa23lk8t/BTj8sYluzkiMgYNR3+vvj4NpenSKi8hgkTIJb7AS9yJjtNsZtJws42NbWwWli9MdlNERMakAScjcTgcLF26lKVLlw5JA37yk5/wn//5n1xwwQUA3H777TzyyCPcfffdfPOb3+zxnEgkwqc//Wmuvvpqnn32WZqbm4ekLZJcNXhxDLL0eL19oGKVY390AFpUNBRNGzF7W/cyLn8ca6auSXZTRGSMG8r+Xn19GkrTGTntBAkSiX9pVVcgp7nnAwzKj9ODpo4mCjILmF2q/DgiI63F34Iv5Bux+2VnZFPgKRix+/Xl3nvv5dJLLx3QmMCyLB5++OFRN7M37kBOJBLhxz/+MX/6058IBoOccMIJrFu3jqysrIRvHgwG2bhxI1deeWXXNofDwZo1a3jxxRd7Pe+aa66hvLycz33uczz77LP93icQCBAIBLq+b21tTbjNMnyipccTz4/jMzl4yQcMJfb+6OAzjfLjROwITf4mzpt9HiXZJclujoiMUUPd36uvT1NpHMgJYQ8gkBMtK17Q04ycSCT6+RXI6aaho4GF5Qspyy5LdlNExpQWfwvff+b71PvqR+yepdmlfOfo7ww4mFNdXc21117LI488wt69eykvL2fJkiVceumlnHDCCUyZMoVLL72USy+9dHga3oMdO3YwdepUXnvtNZYsWTJi9x0OcQdyrrvuOr73ve+xZs0asrKyuPnmm6mtreXuu+9O+Ob19fVEIhEqKipitldUVPDOO+/0eM5zzz3HXXfdxebNm+O+z/r167n66qsTbqcMvwBh6vENqmJVw4HZOIVWPRmhNsjMTKtAzv72/VTlVXH05IHlmBIRGUpD3d+rr09TnUur0mxJUTtBwti44siRYxuLNlMIQL6jhxw5SnTcI2MM/rCfpVVLVVlTZIT5Qj7qffVkubLIzsgesfv5Qr4BBXJ27NjBkUceSWFhIT/60Y9YuHAhoVCIv//971x00UW99v8Sv7hfs9x33338/Oc/5+9//zt/+MMf+POf/8xvfvMbbLuHDP/DpK2tjc985jP88pe/HNAa/SuvvJKWlpaur927dw9jKyURH5QeT3xGTl3XsqrqaKLj3LxoMCcN2MamoaOBE6aeQFFWei0HE5HRJdn9vfr6FNFZtSoNZ+RAfPn2vCYfGycOwuTQwwyuYBCyssHjGepmprX2YDu57lyVHRdJouyMbPIy84b9K9Fg0YUXXohlWbzyyit88pOfZNasWcyfP5/LLruMl156Ke7r3HvvvUyaNIns7Gw+/vGP09DQ0O2YP/7xjyxbtgyPx8O0adO4+uqrCYfDPV5v6tRoReClS6OB6GOPPRaAV199lRNPPJHS0lIKCgo45phj2LRp08A/+AiKe/rDrl27OPnkk7u+X7NmDZZlsW/fPiZMmJDQzUtLS3E6ndTU1MRsr6mpobKystvxW7duZceOHZx66qld2zoHli6Xi3fffZfp06d3Oy8zM5PMNPkH/VjVSAcdhMgi8enL9eZAomOrGoIhKC1NmzeJ1e3VVOZUajaOiCTdUPf36uvTVJourfIeCOTE4+CKVT0OF4JBmKhEx4eq99UzIX8CkwomJbspIpKCGhsb2bBhA9deey05OTnd9hcWFvZ43vnnn8+OHTt46qmnAHj55Zf53Oc+x/r16zn99NPZsGED69atiznn2Wef5dxzz+WnP/0pRx11FFu3buW//uu/ALodC/DKK6+wYsUKHnvsMebPn4/7QEXCtrY2zjvvPG655RaMMdx4442cfPLJvPfee+Sl6PLauAM54XAYzyFvJDIyMgiFQgnf3O12s3z5ch5//PGu5EO2bfP4449z8cUXdzt+zpw5vPHGGzHbrrrqKtra2rj55puZOHFiwm2R5GqggwiGjDjXtB/KmIMTHe+LbixIjaRc/bGNTZ23jk8t/BSl2YlVfxMRGSpD3d+rr09TaVq1qi2hQE5zzwcY0mqJ9khpD7azvGo5TkdiYzYRGd3ef/99jDHMmTNnQOdVVVXFzP69+eabOemkk/jGN74BwKxZs3jhhRfYsGFD1zFXX3013/zmNznvvPMAmDZtGt///vf5xje+0WMgp6wsmterpKQk5mXS8ccfH3PcL37xCwoLC3n66af56Ec/OqDPMVLiDuQYYzj//PNj3nb5/X6++MUvxkTafv/73w+oAZdddhnnnXcehx12GCtWrOCmm27C6/V2VbY499xzGT9+POvXr8fj8bBgwYKY8zsjeodul/TSWXo8UW2mkCBZOAhTFN4HGa60GXzVemspyynj2MnHJrspIiLD0t+rr09DaTojp5EOMuLMHNByINFxj6XH7Uh0Jk6KvolNFn/Yj9vpZk7pwP6BJiJjhzEmofPWr18f8/2WLVv4+Mc/HrNt9erVMYGc119/neeff55rr722a1skEsHv9+Pz+cjOjm9pWE1NDVdddRVPPfUUtbW1RCIRfD4fu3btSuizjIS4AzmdUa6DnXPOOYNuwFlnnUVdXR3f/e53qa6uZsmSJWzYsKErKeKuXbtwpNkgQgZuP+0440/Z1E29iUZUi61anMEOyPSkxeDLGEOtt5Yz5p5BRW5F/yeIiAyz4ejv1denoTSdkRMN5Ay09HhPiY6DSnTcg3pfPeU55cwonpHspohIipo5cyaWZY1IQuP29nauvvpqPvGJT3Tbd+js4r6cd955NDQ0cPPNNzN58mQyMzNZvXo1wWD8szxHWtyBnHvuuWfYGnHxxRf3OL0a6Foj15t777136BskI8pg2EULOYNIdFxvR/PjlDn2RxMdjxsPztSf8lvnq6M4q5hjpx6b7KaIiADD19+rr08znTNy0ig/jMHQjH8ApccPBHJ6rFgVjCY5jvNt7ljR7G/myIlHkulSPioR6VlxcTFr167l1ltv5ZJLLumWJ6e5ubnXPDkHmzt3Li+//HLMtkMTJS9btox3332XGTPiCy535sSJRCIx259//nl+/vOfd+UI3L17N/X1I1fiPRHp9ZpFRiUvIZroIHtQgZyDKlZFbCguHqrmDRtjDNXt1Rw58UjG5Y1LdnNEREQ+0Fm1Ko0COQEi+AjFFciJGCdeE12C3euMnCIlOj5YMBLEYTlYWLEw2U0RkRR36623EolEWLFiBb/73e9477332LJlCz/96U9ZvXp1j+dceeWVnHvuuV3fX3LJJWzYsIEf//jHvPfee/zsZz+LWVYF8N3vfpf77ruPq6++mrfeeostW7bwwAMPcNVVV/V4j/LycrKystiwYQM1NTW0tLQA0VlEv/rVr9iyZQsvv/wyn/70p8nKyhqipzE8FMiRpKvHh49QwoEc2zhoMNHp+aXsA4u0yI/TEmghz52nSlUiIpJ6IgdyxKRRIKOdICEicQVy2kwBBgcugmTh7X6AbdKmaMJIqfXWUpFbwbyyecluisiY5wv5aAu0DfuXL5RYHtNp06axadMmjjvuOL72ta+xYMECTjzxRB5//HFuu+22Hs/Zv39/TE6aVatW8ctf/pKbb76ZxYsX8+ijj3YL0Kxdu5a//OUvPProoxx++OGsWrWK//7v/2by5Mk93sPlcvHTn/6UO+64g3HjxnHaaacBcNddd9HU1MSyZcv4zGc+wyWXXEJ5eXlCn32kxL20SmS41OOjgzCeBH8cm0wpETLIIEB+qBrcmWmRH2dv216OmHAE04qmJbspIiIiscLhZLdgwNoJEowzkNN6INFxgdXYQ6zqQKJOLauK0dTRxHFTjiM7Q89FJFmyM7IpzS6l3ldPR7hjRO5Zml2a0N/7qqoqfvazn/Gzn/2sx/07duyI+b6nZdSf/exn+exnPxuz7Wtf+1rM92vXrmXt2rW9tuPQ5Muf//zn+fznPx+zbenSpbz66qsx284444xer5kKFMiRpKs/ULHKQWJv/RoOWlZldQQgOwsOWYuZanwhH07LybFTjsVKo7edIiIyRoTDaTUbB8A7oEBOH4mOQ2FwuVJ+LDGSfCEfbqebxZWLk90UkTGtwFPAd47+TsIzZRKRnZFNgUczFFONAjmSdDW0JxjCiao7ULGq1NoPwQBMnJjyg889rXuYUzpH68xFRCQ1HZIIMh20EySMwRnHqKLvQE4Q3BkK5Byk1lvLxPyJzC6ZneymiIx5BZ4CBVZEOXIk+XbSQtYQVKwqdVRHZ0On+Jr2UCREMBLk+KnH43IolioiIimoM9lxGmkniAOw4gjktNh9VawKRZdVZSQ+NhlNjDG0BlpZNXEVGU49ExGRVKBAjiRVkAjVtCec6DhsXDSbUgBKzR5wOlI+0fG+tn1MLJjIYeMOS3ZTREREepaGS6vaCRJv6Kmtzxk5ISgoTLvPP1zag+3kunNZXKFlVSIiqUKBHEmqhkFWrGowFRgcZNNGTrAeMjNTOpBjG5uWQAsnTD1ByQJFRCR1RSJpNyOnlWBcx4VMBj6iRRF6DOQYkxZFE0ZKjbeGqYVTmVo0NdlNERlzDk3UK6NfvP/NFciRpBps6fG6g5dVBYLRgVdm5lA2cUjVemspyy5j1YRVyW6KiIhI79KwalUTHWTEMbTtzI+TiY9Myx+7s3MArfw4QPQfFL6Qj9UTV+Ow9M8GkZGScWBpp883ckmNJTV0/jfP6Gd5rxJ0SFLV48PG4EowpvhBfpz90anQxSVD2bwhV+ut5bTZp1GaXZrspoiIiPQuDQM5jXQMvmJVOAQZqljVqcnfRKGnkIXlKs4gMpKcTieFhYXU1tYCkJ2drUq3o5wxBp/PR21tLYWFhTidffdnCuRIUtUxuChz/YGKVWXWvuiGFE503BpoJScjh9UTVye7KSIiIn1Ls6VVBkMz/rgCOW2mEOglkBMMQYZbgZwDar21LKtaxri8ccluisiYU1kZ/XdOZzBHxobCwsKu//Z9USBHkmoPrXENunrSYbJpN4WAoSS8J/oGLYXz4+xr28eiikXMKJ6R7KaIiIj0LRhMq2S/HYTxEx7YjBxHc/edoRAU5KtiFRCMBInYEY6ceKRmAogkgWVZVFVVUV5eTigUSnZzZARkZGT0OxOnkwI5kjQ2hj20Jpwfp96ORioLrAbcoVbI9KRscsJQJETYDnP05KO1xlxERFJfMAiO9Omv2gkSJEIB/efJaz1Qejyvx4pVwWjFKmFv616mFE5h+bjlyW6KyJjmdDrj/se9jB3p00PLqNOCn1YCgwjkRPPjlDn2QyAARUWQor/kqturGZ83nqWVS5PdFBERkf6FQmkVyGkjQIAwmXG8o+xzaZUhZV8KjaSIHaEt2MYJ007A4/IkuzkiInKI9OmhZdQZdMUqcyDRsVUdXctfXDyUzRsyxhia/E0cNfkoctxacy8iImkgzWbktBEkQKTfpVUhk0EHuQDkWc2xO1Wxqkt1ezWVuZWqsikikqLSp4eWUaceH0EiZCaQI8eYD5ZWlVp7o+v4UzQ/TpO/iYLMAlaMX5HspoiIiMQnDWfkADjoO5dL52ycaOnxQOzOUCiaGyc7eziamDaMMTR0NHDslGMp9BQmuzkiItKD9OmhZdSpw4cBrH4GXT1pM4UEycJBmOLQXnBnpuxU6Oq2ahZXLmZi/sRkN0VERCQ+aZbsuI1gXMd9UHq8ufvOzkDOGJ+RU++rp8hTxJETj0x2U0REpBcK5EjS7KcNZwJBHPhgWVWJVYMj0AHZWSk58PKH/ViWxYcmfUgVH0REJH2k2dKqVgJEE9z0c9yBGTm9JjrOyRnzFauqvdWsnLCS8fnjk90UERHpRfr00DLq7KKFnEEmOi517I8ONktKU/LNYXV7NZMKJrGoYlGymyIiIhK/cDitAjn1+HDFsVS7rb/S44WFQ9uwNNPibyHblc0xk49JdlNERKQP6dNDy6jiI0QDHYknOj64YpWxU3LgZYyhNdDKUZOPwu10J7s5IiIi8UuzGTn1+OLKuddn6fExXrHKGMPu1t0srlzMrJJZyW6OiIj0IX16aBlVBlOxKmIcNJpyAErtPdGS4ymY6LjJ30Shp5BlVcuS3RQREZGBSaNkxwZDEx39VqyCPkqPd1asGsOJjmu8NRRkFnD6nNO1HFxEJMWlRw8to05nICcrgUBOkynHxkUmPvJCNZCZmZKBnJr2GuaXzWd8ntaYi4hIGjEmGshJk3/MewnhJ0wmrj6PC5kMfERn3HRLdjzGEx0HI0FqvbWcNPMkZhTPSHZzRESkHwrkSFLU4QX6LxPa47ld+XGqsQIByMsHd2otXQpFQtjG5oiJR+itloiIpJdIBGw7bWbktBEgQKTfpVXtpgAANx1kWv7YnaEguMduIGdb0zbmls7l5BknJ7spIiISh/TooWXU2U97gvWqoN6uBKDU2h99g1ZaOnQNGyI13hqq8qqU5FhERNJPOBydlZM2gZwgAcL9Lq36oPR4D/lxgiHIyQVX37N6RqMGXwMZjgzOnH8mOe6xGcgSEUk3Y6+3kpSwnabEEx0fKD1eZu0HrJRbVmWMobGjkTXT1mhAJCIi6adzRk6alOFuI0AQewCBnObuO8PJq1gVsSP4Qj4MBmMMBoPb6SbLlTXss3ojdoQ9rXs4dfapLK5YPKz3EhGRoaNAjow4HyFq8JLDwJdDBYyHVlMCQGlkZ3SQmWKBnLZgGzkZOSyvWp7spoiIiAxcOJxmS6uCWIDVz1zfzkTHeT0FcgwjvqyqI9TB3ra9dIQ6yMvMw8LCsiwsLPxhPx3hDtxON8VZxRRnFeNyDO2w3TY279a/y/Ti6Xx8zse1FFxEJI0okCMjrhYvXoJUkjvgczvz4+RbjXiCzeDxQO7ArzOcqturmVs6l+nF05PdFBERkYFLu0BOIK7jWjsrVjkOXVo1shWrfCEf25u347JcTCuaxnFTj2Ne2TxcDhcWFg7LQWNHI1ubtvJW7Vv8u+HfbKnbgtvpZnz+eHLdgx/3ROwIW+q3MDF/Iv+1/L8oyioagk8mIiIjRYEcGXE1tNNBOKGKVXX2OADKrH0QCEJ5RbT8eIqI2BECkQBHTjoSh5UeA2AREZEYnUur0mSGRnO8gRw7GqzIOzRHTigczY0zAoGc9mA725q2sWrCKk6cdiILyheQ4ew+HirJLmFmyUxOmnESLf4WXqt+jad3PM2/G/5NIBygPLec0uzShMYaYTvMlrotTCuaxpcO/xLTiqYNxUcTEZERpECOjLhavBhMYhWrzIFAjmM/2BEoSq03SPW+ekqzSllauTTZTREREUlMmiU7rsfXb36csHHhI7oUu9fS48McyGkLtLG9eTvHTTmOzy37HNkZ8d2vwFPAsVOO5ahJR/FO/Ts8v/t5Xt37Km/WvklBZgFVeVW4nfEtV+8IdfBe43vMLpnNhYdfyMSCiYP5SCIikiQK5MiI200rrgQKphnzwdKqMnaD5Ui5/Dj1vnrWTFujKcoiIpK+0mxpVUMcgZy2A6XHM/CTSUfszlAIsrPAPfDcffFq8bewq2UXJ047kQuWXoDH5RnwNZwOJ/PL5zO/fD4fm/0xXt7zMs/sfIb3Gt7DYCj0FFKSVUJWRlbMeZ1FGKrbq3E6nCyuWMx/Lf8vqvKqhurjiYjICFMgR0aUwbCdpoQSHTebEkJk4iJIUWgvZGamVCAnEA7gsBwcNu6wZDdFREQkcZ1Lq9IgkBPBphk/mQMoPd5txVg4BPmVw7aUrD3Yzq6WXZw04yTOXXwuma7MQV+zMreS0+acxoenf5jXa17n7bq3eb3mdXa37u4aj5gDuX9sY1PoKeToyUfzoUkf6nU5l4iIpA8FcmRENeOnGT85g8iPU+qoxuHzR5McZ2X1c9bIqfXWMi5vHPPK5iW7KSIiIolLoxk57QQJxJF3r62v0uMRe9heDBlj2Na0jaMnH835S84f8gBKVkYWqyasYtWEVQTCAbY3b2dn805CdqjrGJfDxYLyBYzPG6/KVCIio4QCOTKiavHSTpBSBr4OvSs/jrUPgkEoLU2pRIxN/iZOnH5itynNIiIiaSWNAjltBAkQoYC+lyq1dpUeP7Ri1QHD9GJoX9s+SrNL+eTcTw77LJhMVyZzSucwp3TOsN5HRESSL/V7aBlVavASJNLvWvaedFWscuwDDBQUDHHrEtcebCc7I1tJjkVEJP2lUdWqNgIEifS7tKprRo6jOXaHHQGHNSyJjgPhAI0djXx05keVVFhERIaUAjkyompoxwKsAVasCphMmk0pAGX2rmiZ0BTKj1PTXsOUwinMKJ6R7KaIiIgMThpVrWojSIhIv0UUWu1CIJojJ8YwVqza1rSN+eXzOXH6iUN+bRERGdtSv4eWUWUHLbgTWNFXf6BaVZ7VTFagKaUSHRtj8Ia8HDHxCJyOgc80EhERSSlpNiPHOvC/3oSNE++B0uN5PZYedw/50qrGjkbcTjdnzDtDS65FRGTIKZAjIyaMzW5ayE2gYlVMfpxAAAqLom/QUkCTv4lCTyGLKxYnuykiIiKDFw5H/0yLQE6w32PaTQFgkUEAD77YnaFQtHiCc+hexNjGZnfLbo6dciwLyxcO2XVFREQ6KZAjI6YeH20EEqpYVXtwfpxIGEpKhrp5Catpr2Fu6VzG5Y1LdlNEREQGLxxOiyAOQBMd/S7WbjXFQHQ2TvfS4+Ehz7lX015DZW4lH5v9MVWJEhGRYaFAjoyYWrx4CZEzwBk5xnywtKrM2hMdXKbIsqqIHSFsh1k5YaUGayIiMjpEIsluQdzq8fVbQKH1QKLjAqux+04D5OQMWXuMMdT56jhi4hGU5ZQN2XVFREQOpkCOjJga2rEx/SYkPFSLKSaIBychikN7wZ06+XHqffWUZpdq6rSIiIwenUur0kA9PjL7yb3X0lmx6tBEx8ZE/xzCRMetgVZy3bmsnrh6yK4pIiJyKAVyZMRUH6hYNVCdZcdLrWocgY7ogGsI354NRr2vnqVVSynKKkp2U0RERIZGmiytChGhjWC/pcdb7ejSqnzHITNywiHIcA1pIGd/+37ml81netH0IbumiIjIoRTIkRGzjSayEsiP05Xo2LEPggEoLU2JAWYoEgJgaeXSJLdERERkCEUiH8xWSWFtBAkQHsDSqkNLj4eHtPR4IBwgYiIcM+UYLbcWEZFhpUCOjIgOQlTTPshEx3ujA8vCwiFuXWLqffWU55Qzv3x+spsiIiIydNJkaVUbAYJE+lxaFTRuOsgFIK9bICcEnqwhq4K5v30/kwsms6RyyZBcT0REpDcK5MiISDTRsd94aDbRZIEV9i5wuVImP05DRwPLxy0n152b7KaIiIgMnbQJ5AQJEOlzRk7nbBwPXjKtQOzOUChasWoIZs/YxqbF38KxU47F4/IM+noiIiJ9USBHRkQtXnyEyB7gjJxaezwABVYDnmATZKZGouNAOIDTcuqtm4iIjD5ps7Qq0G8RhdbeEh0D2DbkDs3LmM7iByvGrxiS64mIiPRFgRwZEftow2BwDDDdcY09AYAKxx4IBKCgcMimQA9Gna+OytxK5pbOTXZTREREhlaaJDtuI9jvMa0mmui44NBExwxtxaqa9hpWTlhJeU75kFxPRESkLwrkyIjYRnO/yQh7EhPICYejiY5TQFNHEyvGryArIyvZTRERERlawf4DJKmgjUC/x7TavczIidjgcAxJIMcX8pHpymT1BJUcFxGRkaFAjgy7EBG20UQemQM7z2TQYCoAKLd2R98OpsCyqo5QB26nm8WVi5PdFBERkaEXCkWDHCmugY5+Z/q29La0KhQasopV9b56qnKrmFUya9DXEhERiUfq99KS9mrw0oqfvAEmOq63qzA4yaaN3GAduFMjP06tt5bx+eM1YBMRkdEpGEyLQE49PjL7mO1rzAdLq/KtQ5ZWhULgdkPW4GfWtgRaWF61nAxn8pd+i4jI2JD6vbSkvX20JVSxqsZEEx1XOPZgBQKQnQU5OcPRxAFpCbSwasIq3M6BfR4REZG0kAaBHIOhkY4+S4/7ySZEJmDIs5pjd4ZCkJc36M8ZjARxWS7mlc8b1HVEREQGIrV7aRkVhiTRcTAAJaVJT77oDXrJcmWxsHxhUtshIiIybEKhpPe3/ekgTDvBPmfkdFasyrVacFmR2J3h8JDM8u2sVqVZuiIiMpIUyJFht5XGPt+Y9cQ2FnX2OAAqHLuj86MLC4ehdQNT56tjQv4EphdPT3ZTREREhkcazMhpxk8HITx9jC9a7M5lVT2UHochKT3e1NHE0qqlZGcMTfUrERGReKR2Ly1pL0iE7TQPOD9Oo6kgjBs3HRSGq8GVkRKBnNZAKysnrMTlGFhgSkREJG2kQbLjJjrwE+4zkNPaW6JjY4PFoBMdh+0wAAvKFwzqOiIiIgOV2r20pL1q2mklMOCKVZ3Lqsode7ECfvB4kp7o2Bv0kp2RrQGbiIiMbmkQyGnGTwSDq4+hbGei44KeEh1nuAcdyGnsaKQ4q5g5pXMGdR0REZGBSu1eWtLePtrwESKbgVVy+CA/zl4I+KG0FFzJnQVT661lQv4EphVNS2o7REREhlWaLK0CsPrIv9frjJzg0JQeb/A1ML98PoWewkFdR0REZKBSu5eWtJdIomNjoMburFi1G2wbiouHq4lxtsnQFmxjxfgVWlYlIiKjWzic8smOG/D1ud82Fq2mEOghkBMKRfPjDOIFkW1swibM4orFCV9DREQkUQrkyLB6n8Y+16/3pMUUEyAbJyFKInvB6Up6fhxfyKdlVSIiMvoZkxYzcvbj7XN84TX52LhwECbHao3dGQ5BQcGg7t/sb6Yws1DLqkREJClSu5eWtBYgzE6aB5wfp/bAsqoyx36cAR9keQY94BqszmVV04tUrUpEREaxSCQazEnhQI7BUEN7nImOm3FY5tALDLpiVb2vnhnFMyjPKR/UdURERBKRur20pL39XYmOB1axqroz0bG1B/x+KCqOrmVPks5lVSvHr8TpcCatHSIiIsMuEokuaU7hQI6XEO0EE6xYdSCoM4j8OMYY/GE/y6qWYaX4EjQRERmdUqKXvvXWW5kyZQoej4eVK1fyyiuv9HrsL3/5S4466iiKioooKipizZo1fR4vybM/gUTHxsB+ezIAVY7dYEeiiY6TyBtStSoRkcFSX58mwuGUD+TEU3q85UDFqvweK1YNLtFxR7iD7IxsZhTPSPgaIiIig5H0XvrBBx/ksssuY926dWzatInFixezdu1aamtrezz+qaee4uyzz+bJJ5/kxRdfZOLEiXz4wx9m7969I9xy6c9e2ojWk4j/bVWLKaGDXJyEKDM7weFM+rKqOm8dE/MnqlqViEiC1NenkTQI5DTj7zeQ0+uMnM5ATk5O4vf3N1OUVcTEgokJX0NERGQwkt5L/+QnP+E///M/ueCCC5g3bx6333472dnZ3H333T0e/5vf/IYLL7yQJUuWMGfOHO68805s2+bxxx8f4ZZLf6KJjge2FGm/PQmIlh13BXzg8SQ10XHXsqoJWlYlIpIo9fVppHNpVQovGWrGj43B1ccwtjOQU+DoYUZOTvaglmw3+5uZWzoXt3NgS8dFRESGSlIDOcFgkI0bN7JmzZqubQ6HgzVr1vDiiy/GdQ2fz0coFKK4j/LUgUCA1tbWmC8ZXv4EEx3v61pWtfNAfpwicCdvoNS5rGp+2fyktUFEJJ2pr08z4XDKJztuwt/n/rBx0m6is3l7nJFTUJjwvY0xROwIs0pmJXwNERGRwUpqL11fX08kEqGioiJme0VFBdXV1XFd44orrmDcuHExA8RDrV+/noKCgq6viRM1FXa47aeNdoIDSnRsG4vqAzNyqhw7IRJOen6cOm8dE/InaFmViEiC1NenmTRYWlWPr89F262mGLBw48eDL3anMYOqWNUR7iArI4sphVMSvoaIiMhgpW4vHYfrr7+eBx54gIcffhiPx9PrcVdeeSUtLS1dX7t37x7BVo5Nu2gZcKLjelNJiEzcdFBs748OIpO4rAqgLdjGivErtKxKRCRJ1NePsDSoWlXdT+nxZhN9CVRo1R+yQmzwFau68uPkK1AoIiLJ03svOAJKS0txOp3U1NTEbK+pqaGysrLPc3/84x9z/fXX89hjj7Fo0aI+j83MzCQzc2BLfGRwttGEY4CJjj+oVrULR7ADMj1JTXTsDXrJcmVpWZWIyCCor08zKT4jx8ZQi7fPQE6THQ3kFDnqY3eEwuByDTrR8bFTjiXTpZ81ERFJnqT20m63m+XLl8ckL+xMZrh69epez/vhD3/I97//fTZs2MBhhx02Ek2VAYhg8xZ1A86Psz8mP04gGsTp4+3rcKvz1TE+bzzTi6cnrQ0iIulOfX2a6QzkpGiy4zYCeAnGPSMnxiBLjxtjsI2t/DgiIpJ0SZ2RA3DZZZdx3nnncdhhh7FixQpuuukmvF4vF1xwAQDnnnsu48ePZ/369QDccMMNfPe73+X+++9nypQpXevrc3NzyR3EmmcZOtW0U4+PYrLiPidkMqi1xwEwzrETwiEoKxuuJsalNdDKR2d9FJcj6X9NRETSmvr6NJLiS6s6S48X9THGaOorkJOTnXARhY5wB5nOTCYXTE7ofBERkaGS9H+hnnXWWdTV1fHd736X6upqlixZwoYNG7qSIu7atQvHQYOJ2267jWAwyBlnnBFznXXr1vG9731vJJsuvdhJC+0EmUT8y6Jq7fHYuMihhTzTAJYjWrEqSXwhHx6XR8uqRESGgPr6NJLiVauaDgRyMuk5d13IZNBuCoGellYFoWBcwrONmv3NFGcVM6lgUkLni4iIDJWkB3IALr74Yi6++OIe9z311FMx3+/YsWP4GySDsp0mDAZHAvlxxjl3YgX8kOVJaqLjOm8d4/LGMaN4RtLaICIymqivTxMpniOnGT82Bmcv2QFaTAkAHrx4rI7YnbaBvLyE793ib+HoyUcrP46IiCRdavbSkrYMJqH8OPsOSnRMhx8Ki5KaH6cl0MLh4w8nwxl/1S0REZG017m0KkVz5DTj77OQQq+JjgdZscoYQ9gOKz+OiIikBAVyZEjV4aOadgoGEMjxmywaTXR6fZVjJ0TCUF4+XE3sV0cougZ+QdmCpLVBREQkKcLhaBAnRQM5tXj7nPHba6LjcASczoQrVnWEO/C4PEwpnJLQ+SIiIkNJgRwZUjtppo0A+QMI5FTbEwEosurIslujA60kLquq99VTmVfJzJKZSWuDiIhIUkQiyW5Bn/bT3nfp8QOBnKIhrljV7G+mKKtI+XFERCQlKJAjQ2oHzRjode16T/baU4EDs3E6/NElVUlMdNzsb+bwqsNxOxOraiEiIpK2wuFkt6BXEWzq8ZLVV+nxA0urCntKdOzxJLxsu8XfwtzSucqPIyIiKUGBHBkynflx+hpgdTvHwO7IdAAmOLaC3w8lpdG3ZkngD/vJcGawoELLqkREZAzqXFqVgloI0EG41xk5AZOJj2gy425Lq4LB6EuiBD9bxEQ0U1dERFKGAjkyZJrws5dWCon/bVedqcJPDhkEqHDsBjsCpaXD2Mq+1fvqqcitUDJDEREZmyKR6FuWFNR8oPR4b4Gczvw4ObTgtoKxO20DBQUJ3TcQDuByuJiQPyGh80VERIaaAjkyZHbRQgsBCgYQyNlzYDbOeMc2nOEguFzJXVbV0cxhVYfhcSWvYpaIiEjSpPDSqs5ATmYvgZxeK1Z1BqZycxO6b1uwjfzMfMbnjU/ofBERkaGmQI4MmZ00E8HGNYAfq932DAAmOg8sq8rKSviN2WAFI0GcDicLKxYm5f4iIiJJl8JLq5roAOi1atUHFasaYneEguDOSLhiVVugjcrcSvIz8xM6X0REZKgpkCND5m3q+6wkcag2O58mU4aFzQTH9mggp6w8OisnCeq8dZTnlDO7ZHZS7i8iIpJ0Kb60CnpvW1NviY6DQXBnJjwjxxvyMrt0NlaKBrhERGTsUSBHhkQbAXbQNKD8OHvs6LKqcmsPmVYHGBtKSoarif1q8jexrGoZWRlZSWuDiIhIUqXw0qpavLhw9rjPmA9m5BRZdbE7g0HIzwdnz+f2xRiDMYaJ+RMHfK6IiMhwUSBHhsRWmmjGP6D8OLsOLKua5NwKwRBkuJOWHycUCeGwHCyqWJSU+4uIiKSEcDhlZ+Tspa3Xmb9+sgmQDRgKrMbYneFIwuOLjnAHWRlZyo8jIiIpRYEcGRLvUE8YG3cvb8oOFTRuauzo262JjvfB3wHZ2dE3ZklQ76unLLuMOaVzknJ/ERGRlBAIpGSOnHaC1OElh4we93cuq8q3mnBZB88qMtGvBJdVtQZayXPnMS5vXELni4iIDAcFcmTQIti8xn7yyIz7nL32VGyc5FsN5DuawR+A8nJwJOdHsqGjgSWVS8hxJ5YIUUREZFRoaAC3O9mt6KYWL15C5NBz23pNdBwOR3PvJVqxKtDGlMIpWnYtIiIpRYEcGbRdtLCfdkqIf5Cz+0DZ8UmOrWDbYAFlZcPUwr6F7TAWFosrFyfl/iIiIimjsTElAzk1tOMnTFZvpce7AjmH5scJRT9PgoEcf8TPrJJZCZ0rIiIyXBTIkUH7Nw20EyS3l7dkh7KNxV57GgATne9DR0e07Hhp6XA2s1f1vnpKsku0rEpERMY224aWlpQM5NTixWCweis9bkdfBhX1VLEqOyehz2QbGwuL8fnKjyMiIqlFgRwZtH9RQybOXgdXh6o14wmQRSY+yqx94PNFy4574k+UPJQafNFlVfmZycnPIyIikhK83miOnBQM5OykhYw+K1ZFq152W1oVCkJxcUJ5f9qD7eS6c5XoWEREUo4COTIojXTwbxooITvuc3ZEZgMwwbENB3a07HhFxXA1sU9hO5oQcXGFllWJiMgY194encGSYoEcG8MOmnud+dtiSgiRiYsgBYcGcgyQl5fQfVsDrRR6CqnITc4YRUREpDcK5Mig/JsGmvBTFGfZcds42B6JLmGa5twCfn90Jk6SllU1+BoozipmXtm8pNxfREQkZaRoIKeJDlrw91qxqs6uAqDEqsZhHVQ63bajM3EGkeh4ZvFMXI6e8/KIiIgkiwI5MihbiCYVdMb5o7TXnkKAbDx4qXLsjC6rKiqGnORUi2rwNbC4cjEFnoKk3F9ERCRlpGggp6afilV1JloavMyxP3ZH52dJMJATtsNMLZqa0LkiIiLDSYEcSViAMJupoTDO2TgA2yLRmS9TnVtwWHa0LGhVVUJr1wcrbIexsVlSuWTE7y0iIpJy2tujfzpSa3hYQzshbNy95MipszsDOftidwSD0Vm/2fEv/+4UioRwOpzKjyMiIikptXpqSSvbaKIOb9xlx0Mmg132DACmO9+GQBAy3UlbVtXY0UhJVomWVYmIiMAHgZwUU4O313IKIZPRVXq8xxk5hYUJvSxqC7aR585TxSoREUlJCuRIwt6jET9hsnpZs36onfYsImSQbzVSYtVEl1XlF0BBcpY11fvqWVSxiEJPYVLuLyIiklLa26MloFLMNprw0HOemnq7CrDIoYVsyxu7045EAzkJaAu0UZZTRklWSULni4iIDCcFciQhBsNG9vWaeLAn2yJzgehsHMsi+qZs3LikLKuK2BFso2VVIiIiXZqawNnz8qVkCRJhH229Jzo20UTH3ZZVGQMknui4PdjOrJJZWEkYo4iIiPRHgRxJyG5a2UlL3GXHfSaH/fZkAKY6tkAoBBmu5FWr6lC1KhERkRiNjSmX6LgOL20Eei09/kF+nEOWVYVCkJGRcCDHGMOkgkkJnSsiIjLcFMiRhGymmlYCFJAZ1/HbI3MwOCiz9pLvaI4uq8rNhaKi4W1oLxp8DSwoX0BRVnLuLyIiknJSMJBTgxcfIbJ7mJFjzAelx7vNyAkEwJOZUCAnGAnicrqoyq1KqM0iIiLDTYEcGbAQEZ5nF3m4sXpNPxirs1rVNOcWwIDfD+PGJ2UKd8SOELbDLKtaNuL3FhERSUmhUDRHTsoFctoxgLOHIWu7KcBPDg7ClFi1sTsDfigtS2ic0RZoI9edy7i8cQm2WkREZHgpkCMD9i4N7KKFCuJ7y9VsF9NgKrGIMNX5TjQ3jjsjWnY8CRo7GinJVrUqERGRLu3tB/rn1Ark7KOt11dGdSYaaCm2anFakYP2GLBtKC5O6J7twXZKs0tVDEFERFKWAjkyYJvYTxC7x2nOPXknshSACY5teKyO6GCxqDjhAdZg1fvqmV82n+Ks5NxfREQk5aRgIMdg2EYTOb3mx+llWVU4DE5XwhWr2oPtzCyeqUTHIiKSshTIkQFpI8DL7KGErLiOD5hM3o8sAGCuc1N0QXsoDBMnJq1aVcREWD5u+YjfW0REJGWlYCCnnSD1+HqvWNVbomN/ALI8kJ+f0H1tYzOxYGJC54qIiIwEBXJkQN6gllq8lMVZrer9yELCuCm06qhy7IIOH2RnQWXlMLe0Z/W+ekqySlhQviAp9xcREUlJ7e3RmSwZ8c22HQm1eGkn2GPFqrBx0WDKASizDpmR4/dDcUlCnyUUCeF0OJXoWEREUpoCORI3g+EV9uLAIoP+kwfaxmLLgWVV85ybohNwvF6oqIScnGFubc8afA0srlysde8iIiIHa2+PzpRNoeVENXjxE8aDq9u+RlOOwYkHL7lWa+xOOwIlJQndsz3YTp47j6o8BXJERCR1KZAjcavByxvUUE58QZjd9gzaTSGZ+JjmfBsiEcCC8eOHt6G9CNthDEbVqkRERA7V3p7sFnRTQztg9Vgh84NlVftiY0+RMDicCefHaQu2UZRVpDx6IiKS0hTIkbhtppom/BTFmR/n7XA0YDLL+S9cVhi87ZCXC+Xlw9nMXtX76inLKVO1KhERkUO1tUXz2KWQrTTh7mWo2pnouPzQRMd+P3g8UFCQ0D07Ex07LA2RRUQkdamXkriEsXmB3WSTgaPXQqAfaLDLqTGTsIgwx/UaYKDDDxMmJm39fYOvgaWVS8nPTCz5oYiIyKjV0JBS+XFaCfAeDT2+PDIGqu1oMuJuFav8ASgqSjhps21sJhVMSuhcERGRkaJAjsRlM9X8mwbGkRfX8Vsi0dk4Uxz/Jsc6UAkjIwOqkrPmPBgJYlkWSyqXJOX+IiIiKa2xMaUqVm2jiWb8FOHptq/OjMNPDhkEuic6joShtDShe4btMBaW8uOIiEjKUyBH+mVjeIxtGAzZvZQAPZjP5LAtMheAea6N0Y1tbdHEg8XJWXNe76unPKdcy6pEREQOZQw0NaVUIOd9Gglj91hcYVdkBgATHFtxWvYHOyIRcDgSzo/THmwnLzNPFatERCTlKZAj/XqLWv5FDROIb0nS5vAR2Lgot/ZQ5tgfLWdqG5g2LWnVMBo7GlletZwcd3KqZYmIiKQsvz/6lSKBHINhM9Xk9FB2HGCXHQ3kTHK+H7sj4IfMweXHKfQUUpZTltD5IiIiI0WBHOmTOTAbJ0SEPDL7Pb7ZLua9yCIAlmc8E93Y2hKdjTNu3HA2tVeBcACn5dSyKhERkZ60H1gCnSKBnGra2UcbxT3kx2mxi2g1JTiIMN6xPXanPxAN4ni6L8eKR1uwjelF05XoWEREUp56KunTv2lgE/sZH+dsnE3hozE4mOh4jwrH3uha9UgEpk8HZ/fp0SOhzldHZV4lc0rnJOX+IiIiKS3FAjlbaaKVAPk9vEDaZc8EoNKxC7cVjN0ZDkFZ4rNpIpEIkwsnJ3y+iIjISFEgR3plMDzBdnyEKOwh2eChauzx7LJnYmGz3NU5G6c1Wj1i/Phhbm3vmjqaWDVhFVkZ8ZVNFxERGVNSLJDzLvVY0GOVzM78OJMchyyrsu3o8u0E8+NE7AhYKD+OiIikBQVypFc7aOYV9sZVqcoY+GfoGABmOt+g0NEIdgRCYZg2HVyu4W5uj9qD7XhcHpZWLk3K/UVERFJee3s0EJKkmbMHCxHhDWop6OEFUofJps5El2lPPDQ/js8H2dnRl0cJ8Ia85LpzGZeXnGXgIiIiA6FAjvSoczZOC4Ee16gfapc9kzozHhdBlriej25sbYXCApgwYZhb27ua9hqmFk1lZvHMpLVBREQkpbW3J7sFXXbSQj2+HsuO745MByxKrf3kWIe02eeDqnGJ58cJtFHgKVCiYxERSQsK5EiP3qaOp9nJePKwepjafLCwcbIxfDQA85z/JNvyRt/sBYLR2ThJmqptjMEb8nLkxCNxOpL/llFERCQlpVAgZyuN+AiRTUa3fZ35cbrNxolEosuqKisTvm9bsI2phVNxOZIzg1hERGQgFMiRbvyEeYi38ROmhOx+j98UPppWU4wHLwtcr0Y3trREZ+NMnDjMre1dk7+JQk8hiysXJ60NIiIiKa+lJRoISQFvUocbZ7eXSCGTwT47moi4W34crxfycgeV6DgUCTGtaFrC54uIiIwkBXKkm8fYyhvUMp3+15nvi0zi7chhAByZsSFaQSIYjFarmjMXMvsvWT5catprmF82X4kLRURE+tLQkBKJjtsI8B4NPS7p3mdPwcZFntVEoVUfu9PfAePGQ0b3WTzxiNgRHJZD+XFERCRtKJAjMXbTwp/5N0V4yKTv6cUBk8lzoZMBmO3czETnNsBAc1N0QJXE2ThhO4xtbFZOWImVIm8ZRUREUlJTU0oEcrbRRBMdPebH2RmZBURn48R06+FQNEnzIJZVeUNectw5CuSIiEjaUCBHukSw+R1bqMfH+H4qVRkDL4ZOxEce+VYjh7meiu5oa4esLJg3DxzJ+/Gq89ZRllPGwvKFSWuDiIhIyotEokurUiCQ8zb1hLHJIDavXbvJZ7s9G4CpzndiT2r3Qn4+lJQkfN+2QBuFnkLKc8oTvoaIiMhIUiBHurzIHl5iD1Mp6jfB8TZ7LjvsuVjYHJXxCBlWKDoY9Plg5iwoLByZRveioaOBw8YdRoGnIKntEBERSWk+HwQCSQ/kNNHBs+zscVnVm+EVGJxUOXZQ6qg+aI+Jtn3CxEGVTm8LtjG9aLoSHYuISNpQIEeAaJWI/+VN3DjJpe/BXJ1dyYuhDwOw2PUCZZ2DqqYmKC+D6dOHu7l98of9uBwullctT2o7REREUl57ezS3XZIDOS+yh2raqSQ3ZrvP5PDvSHR27WLXi7EnBYPgzoCKikHdOxwJM7Vo6qCuISIiMpIUyBHq8XEnm6ilnakU9nlss13MY8EzCOOmyrGDRc6Xoju8XnA6YO68hJMNDpX9bfuZXDCZuWVzk9oOERGRlJcCgRwfIR5nG/lk4jxkaPpmeAU2LsqtPVRYe2JPbPdCYREU9V+coTcRO4JlWcqPIyIiaUWBnDHOR4i72MS7NDCb0j6XVLWbPB4NnkmALEqtfRyf8QccloGAH3xemD1n0G/FBss2Nm3BNo6efDRuZ/LX+4uIiKS0zkBOEl/CvMJedtLSLT+f32TxbmQxEJ2NE5Pk2BgIhWDChEGVTm8PtpPrzlUgR0RE0ooCOWNYBJvf8C9eYg8zKcbVx4+D32TxaPBMfORTYNWzxv37aF6ccAiaW2DadJgzZ1CDqaFQ76unJKuEw8cfntR2iIiIpIX29uifSSpQECLCY2wjE2e3JMdvhQ8jQgYl1n7GOXbEntjaCnm5MG5wAZi2YBuFWYWUZZcN6joiIiIjSYGcMSqMzUO8zaNsZQqFZNH7m7hWu5ANwU/RakrIoYUPux/CY3WAHYGGhujbsEWLklqlqlOtt5bDxx+uyhMiIiLxaGpK6u03U817NDCB/JjtAePhnchSABa7Xop9TxQJR5Mcz5oNOTmDun9boI2ZxTNxOhJPliwiIjLSlJ5/DPIR4le8zj/YRiW5FODp9di9kSk8HTqVIB6yaOfD7ofIsdrB2FDfAGVlsHRp0vPiAHiDXtxON6snrE52U0RERFKf3w/PPjvoYEiibAyPsx2g2wul18OrCZFJkVXLRMf7sSc2NUfHH1OmDLoNERNhSuHgryMiIjKSFMgZYxrp4C428SJ7mEIh+WT2eJwx8GZkBZvCR2FwUGbt5Tj3H8m2vNE3YfUNUFwMy5ZDdvYIf4qe7Wvbx4ziGUpyLCIiEo+NG2HHDpg5Mym3f4ta3qCm22ycrZF5vB05DIClrudiZ+ME/OCwosu5XYMbxobtMJZlUZVbNajriIiIjDQFcsaQbTRxD6/xJnXMorjX5VQtdhGvhI9nrz0NgJnO11nlehynFYlOZW5uhqoqWLYMcnN7vMZIi9gR/GE/R08+GpdDP9YiIiJ9ikTg8cejwZAkVKyqw8v9vEEIm7yDXirV2uN4PrQWgEXOF5nk3HrQWSaal2/qVKisHHQb2oPt5LnzlOhYRETSjv7FOwZ4CfI33mcD79OMn3mUdksoCBA0bl4Pr2ZLZDk2ThxEWOF6gtnOzdG3Yd528Ppg2jRYvDippUoPVeutpTynnMPGHZbspoiIiKS+N9+ELVuiee5GWLRi5mu8SwPz+CDJcLvJ54ng6di4mOT4N0tdz8We2NYenQU8RMUV2gJtFGUVUZajRMciIpJeFMgZxQyGTezn92xhC/WUkc18yrqVGPcbD+9HFvJm+HD8RNfJT3C8zwrXk+Q7mqPlPZubonlwFiyIDqBSILHxwep99Zw25zQKPYXJboqIiEhqMwaefBLC4RHPjxPG5le8zivsZRYlXRUzg8bN48GP4yeHYquGozL+Ghur6eiIfi1cCPn5PV98gNqCbRw+/nAcVmqNaURERPqjQM4gRGzDazsaKW/04XT5qMrJwXHIGyLbGPY2deANhslxuxhflNV1TOe+Vn+IrbXthCKGDKdFUY4bh2XhyXCS7XaSmxl73sHntvlDVLf6MQa8gTAG6LBD7M9rZmP2XrZmNmAwTLYL6LAs3gq1ggGn06LFGk+dZwVN7kXYVvRHIdc0MNd+lKLAezSFwtQ1NOGwIzTll9BYPo3a1hyc/9yNy+kk2+3AsixyM114XE7q2gO0+kO4HBaV+VlMLMlmwiHtHqi+nl+nOm8due5cVk1YlfB9UknENryyvZHaNj/leR5WTC3G6UhuWXcRkbFqqH4nD/Q6hx6/cHwBN2zYwo4GH5OLszlxTgWNHUEavUGKczOpzPewfHIRG3c29X+PbduwN71GdXYRrdWtZLmc1HuDtHaEyPO4sIDmjhDeQBgAt8tBaW4mHcEwu5s6AKjI9zCtNId9zX5qWv24nBZVhVlgDNvrvYDFtLIclk4qwumwsI1hT5OPh613+HPmu4yzc9nZ7qMjFKHVPZs9eacTdBTgstuY3HYf20wrBrANFNkBPOEg3inT2dyRS2DjbjKcDmaU55LnyehxbBAP29hMyJ/Ei1sb0rLP1XhBRGTsSolAzq233sqPfvQjqqurWbx4MbfccgsrVqzo9fiHHnqI73znO+zYsYOZM2dyww03cPLJJ49gi2HDm/u5+s9vE9mzl29sb6SxJoxjRzvHzi5jRnkeAO/XtvHUu3W0HxgIAeRmujh2dnQK76H7utR7u23qPG9GeV6P1zUYWjM72J/bwnvFNTRmtQMWBd4sPJEMWohgW0687pl4PbNpz5xD2FXYdX5mcA/FHS+T792IPxKgNdiB09i0ZuawtWwC+/PKMCELQoG4ns/upg5e3dmIJ8PJmrnlXc9kIPp6fp3X6wh1UN1ezZnzzmRmcXKSNQ6lzp+r/S3+rm1VBR7WnTqPkxYoGaOIpK907usH+zt5oNfp6fiDPQv8+uVd3bY7rGjgo797/Ov+P9OwcRuvF04Aqz3OT9EW811tW4A39rbEbNt2yPhlX0sHz71fz7SyHHZ423g1bxdvle0ly+umNhgm7MilpuBU2rIWA5ARbmB80/10hBqIhosMuQEfLbbNq6WT2BYqgsaObvc7dGwQj7AdpqY1wLce2kNr2wcPLV36XI0XRETGNssYY/o/bPg8+OCDnHvuudx+++2sXLmSm266iYceeoh3332X8vLybse/8MILHH300axfv56PfvSj3H///dxwww1s2rSJBQsWxHXP1tZWCgoKaGlpIT+B6bkb3tzPl369CQOUtzXwjaf/h8bsfHzuLAA+uijagf7lX/sHfO3+dL5tAwg6wjR7fNRlt7E7r5HGLC9+V4iMiIvcUDnGUUbAVYU/Yzx+93gCrgqwPsiNY5kQeR1vUOR9gXz/NjzhEJnhECGni6asPPbll1GbW0zQOfjS4h9dVDWgAdb7tW19Pr+PLqpiWlkOb9a+yYrxK7h01aV4XL2XUU8HB/9cHazz3dpt5yzT4ExkFBts35TK0r2vP9hAfycP9Dq9HZ+Inu7x5GMbqf7atwg7nDTkFA7BXfoWsWy2F9bzZukeGrO95AbysKwFtGfNp80zH9uRBSZCsfc5Stsew2FCALjsMHkBL2GHiy3lU9lTUNHvvQYy1ti8Zw+PvbuLotD5OCno2p4Ofa7GCyLpazT39TKykh7IWblyJYcffjg/+9nPALBtm4kTJ/LlL3+Zb37zm92OP+uss/B6vfzlL3/p2rZq1SqWLFnC7bffHtc9B/MXKGIbPnTDE11vQHoK5ORmRic69TjbJk4GsC2LiOUk4IS2TIPXbdPuhqZMQ1OWTVumRdCVRdiRg0UBDgqIOIoIOUswjp4TEWeEG8jr2EKh7y2Kfe/gjnRgYRF0uujI8FCTW0RtbjEtnlxg6Kbn5ma6+OyHpsY19dk2hruf297n88vNdHHMvAgl2SV880PfTPuKE4f+XB3KAioLPDx3xfGaNi0ySo3mwV269/WHivd38kCv09/xiei6x7lz4bnn+MNPfk12Ux3/LpuMGYbcMAYIOhzU5USozQ6zs9CmLjcL2yrHyTh87ukxY5TM4F6qWn6HJ7QPMHhCQXJCHUQsB81Z+WwtHk9dbnFc9453rNHsb+G+V14lEqyiMHwuFrHPIZX7XI0XRNLbaO7rZWQldWlVMBhk48aNXHnllV3bHA4Ha9as4cUXX+zxnBdffJHLLrssZtvatWv5wx/+0Ot9AoEAgcAHS4JaW1sTbvMr2xtjOs+G7AJ+dPS5hJ0OjOXoejsS/dPCWB98bywruv3An8aysC0Lw4E/LceBbU6M1b2q1IAYG3ekmazgfvL8u8gN7CQ3sBN3uImww0XY6aTd7aElq4QWTy4tnhy87iyGMnhzsPZAmL1NHUwszu732L1NHf0GwVoCjdS3e/jyyv9I+yAOdP+5OpQB9rf4eWV7I6unl4xcw0REBmk09PWH6vydfOyPnyLH3Xt/7Q1GBnSd/o7vfoEP3sVZB65omeifTmPjMDYO2+CtbmPtjdUQidC+/JOEnb232fQwDjAf3OCD8QtgLLAtBxGHA/vA+CXicPU7hskIN1Lge5MC3xsU+v9NRiSE09gABJwZ7CqoZF9BGY1ZBZgB5L7pb6zhDXrZ2bKTZp+NCcwl117VLYjT+RlTtc/VeEFERCDJgZz6+noikQgVFbHTZSsqKnjnnXd6PKe6urrH46urq3u9z/r167n66qsH32Cgti2284w4XdTklw7JtftkglgHvhymA8v4cBg/lvHisFujX6YNyzThiNQDjUQcEHY4qcnKYG+Ok7DDSdhRRNDpJuhyEekaaPkOfA2vt+ubaQ5l9XvcvpYOglZLn8dYOJlfchIrx68cquYl1aE/V4M9TkQkVYyGvr43uxuHpu9M+DoHBTkOfYFk80EwpQMPXWGt/rvhoWFsLNOB027BFanGGanFZdeSEXofZ2QXWBZtbovmTCcdGZl0ZGQRdGbQkZGJ3+UGWg58DUxfYw23082SyiVk24fz7hZ/j0Gcg6Vin6vxgoiIQIokOx5uV155ZcybvdbWViZOnJjQtcrzDsnDYiIUeTcRdDYe2GAdeBvWWeQ7+r2FhWUsOv+HceAwDiwsHMZxYAR2YBhmbA5cBLAx2Adm8URn7dgHZvJEvxxErEJCGaVELAdhh5OQ04Xdx3RpN+COJPTxB+XzixeycEJBv8e9saeFq7a/0c9RTk6b9QmsQVTESiXdfq4GeZyIyFgzrH09kNuxmSLfO1gHzYKZW5VPflbveeRaO0Js2d//zKDO6/R1fMxMGeugwE3nfoeFfWCcYLAIulyEHK6uGS2rpkVfOr20rR56zMBzaH9qffB/Dx3XHBi7WMaBw1g4jQun7Yi20kSiX9gHLmMRsSxsKwtjTSXimoY/M/pCKejMIOjKiPlsmUBm4ivT+xxrZLmymFc2j1e2N2PxUr/XSsU+V+MFERGBJAdySktLcTqd1NTUxGyvqamhsrKyx3MqKysHdDxAZmYmmZmZg28wsGJqMVUFHqpb/AfGNU6acg/v2t+5NtkYQ01rYEiSFR7s0IoUPcmAnsdoSVRV4OH8w+Nbr728ynDnE44PnvEhOp/xh2aMnkR+3X6uDtH5mVdMjS9PgIhIqhgVfT3QnrWE9qwlwAe/k38eZ46c/n63//yQHDm9HZ+Iznv88svHAwz59QfLPUQN6fyc8Yw10rnPTee2i4jI0Bn6LHcD4Ha7Wb58OY8//njXNtu2efzxx1m9enWP56xevTrmeIB//OMfvR4/1JwOi3WnzgN6e3cF606dx/c+Nr/PYxJhAf951NRhymIzvNadOi/upHvxPuPRlMRvLH5mERkbRnNf39/v5IFep6/jE3HoPYb6+qkm3n4ynfvcdG67iIgMnaQGcgAuu+wyfvnLX/I///M/bNmyhS996Ut4vV4uuOACAM4999yYBIlf+cpX2LBhAzfeeCPvvPMO3/ve9/jnP//JxRdfPGJtPmlBFbeds4zKgthpq5UFnq6Sj30dc/s5y7j9nGVUFcQ/7bXqwLWvPHketw3w3Hhl95GwMVFF2RncnkAZzHie8WgzFj+ziIwNo7WvH47r9HZ8PA79t3tP9xjM9RN14rzyYRm3dKpKoJ9M5z43ndsuIiJDI+nlxwF+9rOf8aMf/Yjq6mqWLFnCT3/6U1aujCaxPfbYY5kyZQr33ntv1/EPPfQQV111FTt27GDmzJn88Ic/5OSTT477fkNV9i1iG17Z3khtm5/yvOg01kPfgPR1TOe+PU0+Hn2rGl8wQpbbwfSyXFwOB4XZbkpz3VQWZHW7due5+5p8bN7TTMQ21LYFsI0hEIpQnO2m3hukvi2AZcGsijxyM128V9uGMZCf5aLVH8FhWZw4r4ILjpyK02HxyvZGqlv91LcFqGnt4LEttYTCEbLdLsYXZlHvDeJ2gCczg9KcDBwOB1UFHgo8GWzZ38aeZi+eDCeLJxRx5MxSVk0rGdRboXie8WgzFj+ziIz+kqSjua8fjuscevzC8QXcsGELOxp8TC7O5sQ5FTR2BGn0BinOzaQy38PyyUVs3NkU1z0Ovn5xlpt3atrY3eRjfGEWloGdTT5qWjvAQK4ngzmVeTR6A7ywNZoTcNGEAk6cU8GrO5v4155mstxOlk0uAtvw+Lu1wAfjC7fL0XW/6lY/je0Bcj0u/vFWNVvrvLgcMKEoG8uyqG0LUJ7nxmFZGGPwhw0LxxdQkptJQVYG/3i7Gl8wTLbbxdr5lYwvyh5UP5nOfW46t11krBrtfb2MnJQI5Iw0/QUSEZFUo75paOl5iohIqlHfJEMl6UurREREREREREQkPgrkiIiIiIiIiIikCQVyRERERERERETShAI5IiIiIiIiIiJpQoEcEREREREREZE0oUCOiIiIiIiIiEiaUCBHRERERERERCRNKJAjIiIiIiIiIpImFMgREREREREREUkTCuSIiIiIiIiIiKQJV7IbkAzGGABaW1uT3BIREZGozj6ps4+SwVFfLyIiqUZ9vQyVMRnIaWtrA2DixIlJbomIiEistrY2CgoKkt2MtKe+XkREUpX6ehksy4zBcKBt2+zbt4+8vDwsy0p2c2K0trYyceJEdu/eTX5+frKbM+ro+Q4fPdvhpec7vFLh+RpjaGtrY9y4cTgcWvk8WKnc1w+FVPiZTWV6Pn3T8+mdnk3f9Hz61t/zUV8vQ2VMzshxOBxMmDAh2c3oU35+vn45DiM93+GjZzu89HyHV7Kfr97ODZ106OuHQrJ/ZlOdnk/f9Hx6p2fTNz2fvvX1fNTXy1BQGFBEREREREREJE0okCMiIiIiIiIikiYUyEkxmZmZrFu3jszMzGQ3ZVTS8x0+erbDS893eOn5SrrRz2zf9Hz6pufTOz2bvun59E3PR0bKmEx2LCIiIiIiIiKSjjQjR0REREREREQkTSiQIyIiIiIiIiKSJhTIERERERERERFJEwrkiIiIiIiIiIikCQVyUsQzzzzDqaeeyrhx47Asiz/84Q/JbtKosX79eg4//HDy8vIoLy/n9NNP59133012s0aN2267jUWLFpGfn09+fj6rV6/mb3/7W7KbNWpdf/31WJbFpZdemuymjArf+973sCwr5mvOnDnJbpYIALfeeitTpkzB4/GwcuVKXnnllT6Pf+ihh5gzZw4ej4eFCxfy17/+dYRamhwDeT733ntvt7/rHo9nBFs7chIZUz711FMsW7aMzMxMZsyYwb333jvs7UyWgT6fp556qtvPjmVZVFdXj0yDR1CiY+ax8rsnkeczln73yMhSICdFeL1eFi9ezK233prspow6Tz/9NBdddBEvvfQS//jHPwiFQnz4wx/G6/Umu2mjwoQJE7j++uvZuHEj//znPzn++OM57bTTeOutt5LdtFHn1Vdf5Y477mDRokXJbsqoMn/+fPbv39/19dxzzyW7SSI8+OCDXHbZZaxbt45NmzaxePFi1q5dS21tbY/Hv/DCC5x99tl87nOf47XXXuP000/n9NNP58033xzhlo+MgT4fgPz8/Ji/6zt37hzBFo+cgY4pt2/fzimnnMJxxx3H5s2bufTSS/n85z/P3//+92FuaXIkOuZ+9913Y35+ysvLh6mFyZPImHks/e5J9N8UY+V3j4wwIykHMA8//HCymzFq1dbWGsA8/fTTyW7KqFVUVGTuvPPOZDdjVGlrazMzZ840//jHP8wxxxxjvvKVryS7SaPCunXrzOLFi5PdDJFuVqxYYS666KKu7yORiBk3bpxZv359j8f/v//3/8wpp5wSs23lypXmC1/4wrC2M1kG+nzuueceU1BQMEKtSx3xjCm/8Y1vmPnz58dsO+uss8zatWuHsWWpIZ7n8+STTxrANDU1jUibUkk8Y+ax9rvnYPE8n7H6u0eGn2bkyJjT0tICQHFxcZJbMvpEIhEeeOABvF4vq1evTnZzRpWLLrqIU045hTVr1iS7KaPOe++9x7hx45g2bRqf/vSn2bVrV7KbJGNcMBhk48aNMX/fHQ4Ha9as4cUXX+zxnBdffLHb74e1a9f2enw6S+T5ALS3tzN58mQmTpyomaMHGUs/O4OxZMkSqqqqOPHEE3n++eeT3ZwREc+YeSz//MT7bwr97pHhoECOjCm2bXPppZdy5JFHsmDBgmQ3Z9R44403yM3NJTMzky9+8Ys8/PDDzJs3L9nNGjUeeOABNm3axPr165PdlFFn5cqV3HvvvWzYsIHbbruN7du3c9RRR9HW1pbspskYVl9fTyQSoaKiImZ7RUVFr3k5qqurB3R8Okvk+cyePZu7776bP/7xj/z617/Gtm2OOOII9uzZMxJNTmm9/ey0trbS0dGRpFaljqqqKm6//XZ+97vf8bvf/Y6JEydy7LHHsmnTpmQ3bVjFO2YeS797Dhbv89HvHhkurmQ3QGQkXXTRRbz55pvKgTHEZs+ezebNm2lpaeH//u//OO+883j66acVzBkCu3fv5itf+Qr/+Mc/lBxvGHzkIx/p+v8XLVrEypUrmTx5Mr/97W/53Oc+l8SWichQWr16dcxM0SOOOIK5c+dyxx138P3vfz+JLZNUN3v2bGbPnt31/RFHHMHWrVv57//+b371q18lsWXDS2PmvsX7fPS7R4aLAjkyZlx88cX85S9/4ZlnnmHChAnJbs6o4na7mTFjBgDLly/n1Vdf5eabb+aOO+5IcsvS38aNG6mtrWXZsmVd2yKRCM888ww/+9nPCAQCOJ3OJLZwdCksLGTWrFm8//77yW6KjGGlpaU4nU5qampittfU1FBZWdnjOZWVlQM6Pp0l8nwOlZGRwdKlS/V3nd5/dvLz88nKykpSq1LbihUrRnWAYyBj5rH0u6fTYP5Nod89MlS0tEpGPWMMF198MQ8//DBPPPEEU6dOTXaTRj3btgkEAsluxqhwwgkn8MYbb7B58+aur8MOO4xPf/rTbN68WUGcIdbe3s7WrVupqqpKdlNkDHO73SxfvpzHH3+8a5tt2zz++OO95h9bvXp1zPEA//jHP0ZlvrJEns+hIpEIb7zxhv6uM7Z+dobK5s2bR+XPTiJj5rH08zMU/6bQ7x4ZKpqRkyLa29tjIrPbt29n8+bNFBcXM2nSpCS2LP1ddNFF3H///fzxj38kLy+va81uQUGB3jQNgSuvvJKPfOQjTJo0iba2Nu6//36eeuqpUVu2dKTl5eV1W3udk5NDSUmJ8jwNgcsvv5xTTz2VyZMns2/fPtatW4fT6eTss89OdtNkjLvssss477zzOOyww1ixYgU33XQTXq+XCy64AIBzzz2X8ePHd+XO+spXvsIxxxzDjTfeyCmnnMIDDzzAP//5T37xi18k82MMm4E+n2uuuYZVq1YxY8YMmpub+dGPfsTOnTv5/Oc/n8yPMSz6G1NeeeWV7N27l/vuuw+AL37xi/zsZz/jG9/4Bp/97Gd54okn+O1vf8sjjzySrI8wrAb6fG666SamTp3K/Pnz8fv93HnnnTzxxBM8+uijyfoIwyaeMfNY/t2TyPMZS797ZIQluWqWHNBZ2vDQr/POOy/ZTUt7PT1XwNxzzz3Jbtqo8NnPftZMnjzZuN1uU1ZWZk444QTz6KOPJrtZo5rKjw+ds846y1RVVRm3223Gjx9vzjrrLPP+++8nu1kixhhjbrnlFjNp0iTjdrvNihUrzEsvvdS175hjjuk2Rvjtb39rZs2aZdxut5k/f7555JFHRrjFI2sgz+fSSy/tOraiosKcfPLJZtOmTUlo9fDrb0x53nnnmWOOOabbOUuWLDFut9tMmzZtVI+RBvp8brjhBjN9+nTj8XhMcXGxOfbYY80TTzyRnMYPs3jGzGP5d08iz2cs/e6RkWUZY8ywR4tERERERERERGTQlCNHRERERERERCRNKJAjIiIiIiIiIpImFMgREREREREREUkTCuSIiIiIiIiIiKQJBXJERERERERERNKEAjkiIiIiIiIiImlCgRwRERERERERkTShQI6IiIiIiIiISJpQIEdEREREREREJE0okCMiIiIiIiIikiYUyBGRtNHQ0EB5eTk7duwY0Hlf+MIX+PSnPz3kxybiU5/6FDfeeOOwXV9ERAbv8ssv5/TTTx/Saybah0nvEnmmqTQmgJ7HBRoriEh/FMgRkYTt2bOHL37xi8yYMQOPx0NFRQUf/vCHeeONNwAwxlBYWMgtt9zS7dwLL7yQFStWADB//nzWrVvX4z3Wr19PSUkJDQ0NXHvttZx22mlMmTJlQO1cv349v/jFLxI69qtf/Sqf+MQnBnS/vlx11VVce+21tLS0DNk1RURkaG3evJklS5YM6TUT7cPSxUiPCSCxZ5pKYwLoeVygsYKI9EeBHBFJyI4dO1i6dCkNDQ386le/4p133uH//u//mDdvHpmZmQBs3bqVlpYWDjvssG7nb9y4keXLlwOwcOFC3nzzzW7H7N+/n+uuu45rrrmGrKws7rrrLj73uc8NuK3FxcXk5OQkdOwrr7zSY/sTtWDBAqZPn86vf/3rIbumiIgMrddff31IAzk+ny/hPiwdjPSYoKSkJOFnmkpjAuh5XKCxgoj0R4EckQP+93//l6ysLPbv39+17YILLmDRokVD9kZkx44dWJbF7373O44++miysrI4/PDD2bVrF88++yyrVq0iOzubE044gebm5phz161bx8KFC8nJyaGiooIvfelLhEKhEWv7oW655RZycnJ48MEHWb16NVOmTOGoo47ipptuYtasWUB0YOZyuboNhkOhEP/617+6Bm2LFi3qcdD2rW99i6lTp/LFL36Rv/71r2RmZrJq1aqYY8rLy7nzzjtjtr366qt4PB62b9/e9cw7p13bts11113HzJkzu94Ynn/++QAxxwaDQTIyMnjhhRf49re/jWVZ3e59qOrqaizL4uabb2bp0qV4PB7mz5/Pc889F3PcqaeeygMPPNDntUREpLuR6O/27NlDfX09ixcv7tr25ptvcvLJJ5Ofn09lZSVf+9rXCAaDXftffvllPvShD5GVlcWSJUt45plnsCyrq2/rrQ9LdFzQ15hgpJ7TwUZ6TAA9P9NUGhPA4MYFGiuISJ+MiBhjjLFt2yxatMhcfPHFxhhjvvvd75oJEyaYPXv2dDv22muvNTk5OX1+7dy5s9t5f/jDHwxgTjjhBPPss8+aTZs2mYkTJ5qjjjrKnHzyyebVV181L730kikpKTE/+clPYtr2ne98xzz//PNmx44d5q9//aspKyszP//5zwfc9qFywQUXmIqKCrN9+/Zej/n6179uFi1a1G37a6+9ZgCzadMmY4wxf/rTn4zT6TQdHR1dx/zzn/80DofDPPnkk8YYYy655BJz0kkndbvW8ccfb7761a/GbDvuuOPMJZdcYoyJPvPCwsKufT/4wQ/MwoULzRNPPGF27Nhhnn/+eXPXXXd1OzYSiZiXX37ZAGbz5s1m//79pqmpqc9n8re//c0AZtGiReapp54yW7ZsMSeddJKZNGmSiUQiMce53W7j9/v7vJ6IiMQaib76z3/+sykoKOj6ftOmTSYvL898+9vfNu+995558sknTVVVlbnmmmuMMca88cYbJicnx3z72982W7ZsMf/3f/9nysvLTWZmpgmFQsaY3vuwRMYF/Y0JBvqchsJIjwmM6fmZptKYwJjBjQs0VhCRviiQI3KQP//5zyYzM9P84Ac/MEVFRebNN9/s8biGhgbz3nvv9fnVOXg72Pe+9z1TXFxs6uvru7adc845ZsqUKcbr9XZtO+mkk8w3vvGNPtt69tlnm6985SsDbvvBHn/8cfPjH/84rv1bt241f/zjH7v2bdy40UyaNMlYlmUOO+wwc8UVV5i33nor5vzjjz/efPazn+123TvvvNNkZmaaYDBojDFmx44dBjCvvfZa1zEf+tCHzJlnntn1/WmnndbjtS655BKzdu3aru83bNhg8vLyTF1dnTEm+syPPvrorv1HHXWU+da3vtXj5z302IcfftiUlJT0eGxPrr/+epORkREzkP3nP/9pALNr166uba+//roBzI4dO+K+toiIRA13X/39738/pi9Yvny5ufDCC2OO+da3vmVWrFhhjIkGCs4+++yY/aeeeqpZunRp1/e99WFDNS44dExgzNCPC1JpTGBMz880lcYExgxuXKCxgoj0xZWESUAiKeujH/0o8+bN45prruHRRx9l/vz5PR5XXFxMcXHxgK//+uuv8/GPf5ySkpKubbt27eKss84iOzs7Zttpp53W9f3OnTv54Q9/yNNPP83evXsJhUL4/X6uv/76Abf9YMcffzzHH398XPv/9re/0dbWxsc+9jEAli1bxrZt23juued49NFHeeihh7jxxhv5/e9/z6mnngrApk2bOOOMM7pdd+PGjSxcuJCMjAwAJk+eTEFBAW+++SZLlizhwQcfZOPGjbzzzjtd53R0dODxeLpda+HChTz88MNANJHilVdeyde//nVKS0uB7nkOPvaxj3HFFVfwz3/+kzPPPJNPfvKTFBUV9Xjsa6+9FjO1vj+bN2/mE5/4REzSxfz8/G7HZWVlAdGcCSIiMjDD3Vdv3ry563f/O++8w8aNG7vlKnG73QQCAXbu3MmTTz7ZbSlQZmZmTP/RWx+WyLggnjEBDP24IJXGBNDzM02lMQEMblygsYKI9EU5ckQOsmHDBt555x0ikQgVFRW9HnfdddeRm5vb59euXbu6nbd582ZWrlwZs+3111+PWWft9/t59913uwYLdXV1HH744TQ0NPCTn/yE5557jhdeeAGHwxEzoOit7S+//HLXIKrzuM985jNAdADTWU3il7/8JcuWLWPBggWcddZZMfuffvppvvOd73DXXXexdOlSvF4vAE6nk2OOOYZrr72Wt956i/Lycu6//34Adu/eTXNzc48Dx8cee4wjjjgiZtuCBQt488038fv9XHHFFVxxxRVMmjSpa39paSlNTU3drrVgwQL27NlDe3s7DzzwAPv37+eyyy6LeeYHP6fLL7+cLVu2cMIJJ/Df//3fzJgxg+3bt/d47KHf96enKicvvvgipaWljB8/vmtbY2MjAGVlZXFfW0REokair+78Xf7WW2+RkZHRleel09tvv83ChQvZvHkzbre7W1+3ZcuWmP6jtz5soOOCeMcEfT2neMYFqT4m6O2ZptKYoPOcRMcFGiuISJ+SPSVIJFVs3LjR5OXlmfvvv998+MMfNmeccUavxyYyXbulpcVYlmVeffXVrm3btm3rNm32lVdeMQ6Hw7S1tRljjLnrrrtMcXGxsW2765hbbrnFAKa2trbftre0tJiZM2d2fX/EEUeYf//738YYY2bOnGmCwaBpbGw0ixcvNuFw2BhjutZ9d+43xphjjjmmz7Xvfr/flJSUmC9/+cvGmOi0a8A88sgjMcc9+uijBjDPP/98zPYvfelL5pRTTjHf//73zeTJk43P54vZ/6Mf/cgsXry4233b2tqMZVnm+eefN9OnT4/JEdD5zDdu3Nhjmzs6OkxGRob5y1/+0uOxU6dONf/zP//T62c+mM/nM06n0/zgBz/o2haJRMzSpUvN1772tZhj77zzTjNhwoS4risiIh8Y7r66tbU1pi/4+9//bhwOR0yekm3btpmMjAzzt7/9zfz5z382DocjJp/LY489ZgDzxBNPdG3rqQ9LZFwQz5igv+fU37igtrY25ccExvT8TFNlTGDM4McFGiuISF8UyBExxmzfvt1UVlaa9evXG2OMeemll/rs7BPxzDPPGJfLFTPY+/3vf2+Ki4tjjvvFL34RM8D6wx/+YFwul/nDH/5g/v3vf5sbb7zRlJaWmvHjx8fd9smTJ5tgMGgeeeQRc9555xljooPVhQsXdv3/kyZNMpdddlnXGvqD9xsTHcB0Ouecc8x1111nXnrpJbN9+3bz+OOPmxNOOMGUlJSYrVu3GmOiiRbnzJljFi1aZB577DGzefNmc8cdd5jS0lJz/vnnd3s+P//5z01ZWZnJyckxDz30ULf9//rXv4zL5TKNjY3d9k2ZMsWsWrXKzJw5M2ZQ3vnMOwfgN9xwg/mf//kf8/bbb5t33nnHfPWrXzWVlZWmsbGx27Gdz+1b3/qW2bt3r2lubu5234O9/PLLxuVymTlz5pgXXnjBvP322+aMM84wU6dO7ZYQ8bzzzusxT4CIiPRuJPrqZ599NqYvaG5uNsXFxebSSy81W7duNY8//riZO3eu+cxnPmOMMWbfvn0mMzPTXHzxxWbr1q3mT3/6k5k2bZoBTENDQ9d1e+rDEhkX9DcmiPc59TUuSIcxQW/P1JjUGBMYM/hxgcYKItIXBXJkzGtoaDCzZ882X/jCF2K2n3zyyTEJ8wbrlltuMfPnz4/Ztm7dOnPCCSfEbLvoooti3pxFIhHzhS98weTl5Zny8nJz2WWXmQsvvNCccsopcbd9zZo15u233zarV68277//vjHGmBdffNF86lOf6jqmtbXV3HfffWbevHnm4Ycfjtm/e/duc+SRR3Yd+5Of/MQcccQRprS01Hg8HjNz5kxz0UUXmd27d8e0Y+vWrebjH/+4KSkpMfn5+WbZsmXmF7/4RddbvoM999xzBjDHHXdcr89wxYoV5vbbb++2/dRTTzWA+e1vfxuz/ZZbbjELFizo+v7qq682s2bNMh6Px5SWlprTTjvNvP322z0ea4wxv/rVr8y4ceMMYC6//HJjjDH33HOP6Wky4x133GEWLFhg7rvvPlNVVWWys7PNxz/+8ZhkhsZE3/gVFBSYF198sdfPKSIisUayrz60L3jmmWfMsmXLjMfjMdOmTTPr16+P6cd+85vfmIkTJ5qcnBzz8Y9/3Fx99dVmxowZ3a59aB+WyLigrzGBMfE/p/7GBekwJujpmRozsmMCY4ZnXKCxgoj0xzLGmJFeziUiI+uSSy7B5/NhjOGuu+4Cojlx6urq+Na3vsV7773HzJkzAbjwwgs55phjaG1t7dr//PPPc9NNN/HQQw8l82PwyCOP8PWvf50333wThyM5Kb7WrVvH008/zVNPPRWz/aKLLqKpqakrH0BvbrvtNh5++GEeffTRYWyliIgkg23bHHvssXzoQx/iuuuui9mXCn1Yp77GBWeeeWZajAkgNZ7pcIwLNFYQkf4o2bHIGDB37lzuu+8+rrrqqq5tb7zxBgsWLADgBz/4AbNnz2bp0qVYlsWZZ54Zs3/BggVs27aNhQsX8vbbbyflMwCccsop/Nd//Rd79+5NWhv+9re/8cMf/rDb9s2bN7No0aJ+z8/IyOCWW24ZjqaJiMgIe+aZZ/jd737Htm3beOWVVzjrrLPYuXMnl19+ebdjU6EP69TXuCBdxgSQGs90OMYFGiuISH80I0dEZJCMMRQUFPDAAw9w8sknJ7s5IiIyQh566CG++c1vsnfvXioqKlizZg3XXXddn9W0ZPTTuEBEhpsCOSIiIiIiIiIiaUJLq0RERERERERE0oQCOSIiIiIiIiIiaUKBHBERERERERGRNKFAjoiIiIiIiIhImlAgR0REREREREQkTSiQIyIiIiIiIiKSJhTIERERERERERFJEwrkiIiIiIiIiIikCQVyRERERERERETShAI5IiIiIiIiIiJp4v8DFCx3BTboW4kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting results\n",
    "fig, axs = plt.subplots(ncols = len(scales), figsize = (6*len(scales), 5))\n",
    "\n",
    "for ax, scale, x, pars, cov_pars in zip(axs, scales, xs, parss, cov_parss):\n",
    "\n",
    "    ax.set_title(f\"cond. prob. for AE: scale {scale}\")\n",
    "\n",
    "    if scale == \"plain\":\n",
    "        xlab = r\"$x = max_{visit} SUV(visit, p)$\"\n",
    "    else:\n",
    "        xlab = r\"$x = log(max_{visit} SUV(visit, p))$\"\n",
    "    \n",
    "    ax.set_xlabel(xlab)\n",
    "    ax.set_ylabel(\"P(AE|X = x)\")\n",
    "\n",
    "    # plotting data points\n",
    "    for i, lab in enumerate([\"NC\", \"AE\"]): ax.scatter(x[y == i], y[y == i], label = lab)  \n",
    "\n",
    "    # plotting fitted model \n",
    "    #xp = np.linspace(min(x), max(x), 100)\n",
    "    xp = np.linspace(min(x), max(x)+1, 100)\n",
    "    yp = lg.model(xp, pars)\n",
    "\n",
    "    ax.plot(xp, yp, label = \"logistic reg\")\n",
    "\n",
    "    for lab, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
    "        \n",
    "        # define quantile model\n",
    "        fname = f\"lg.get_model_quantiles_{lab}\"\n",
    "        \n",
    "        # get result of model quantiles at probs\n",
    "        # pars and cov_pars are obtained via MLE method\n",
    "        res = eval(fname)(xp, probs, pars, cov_pars)\n",
    "\n",
    "        # make plot of quantiles\n",
    "        ax.fill_between(xp, *res, color = c, alpha = 0.5, label = f\"CI:{lab}\")\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "fig.legend(handles, labels, bbox_to_anchor=(1.02, 0.5), loc='center right')\n",
    "\n",
    "plt.savefig(os.path.join(results_path, \"logit_fit.pdf\"))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67fa2c6c",
   "metadata": {},
   "source": [
    "##  Model CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92080804",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model CI:scale:plain,method:normal\n",
      "model CI:scale:plain,method:delta\n",
      "model CI:scale:plain,boots-method:normal\n",
      "model CI:scale:plain,boots-method:nonparam_boots\n"
     ]
    }
   ],
   "source": [
    "# plotting results\n",
    "fig, axs = plt.subplots(ncols = len(scales), figsize = (6*len(scales), 5))\n",
    "\n",
    "boots_pars_results = dict()\n",
    "\n",
    "for ax, scale, x, pars, cov_pars in zip(axs, scales, xs, parss, cov_parss):\n",
    "\n",
    "    ax.set_title(f\"cond. prob. for AE: scale {scale}\")\n",
    "    \n",
    "    if scale == \"plain\":\n",
    "        xlab = r\"$x = max_{visit} SUV(visit, p)$\"\n",
    "    else:\n",
    "        xlab = r\"$x = log(max_{visit} SUV(visit, p))$\"\n",
    "    \n",
    "    ax.set_xlabel(xlab)\n",
    "    ax.set_ylabel(\"P(AE|X = x)\")\n",
    "\n",
    "    # plotting data points\n",
    "    for i, lab in enumerate([\"NC\", \"AE\"]): \n",
    "        ax.scatter(x[y == i], y[y == i], label = lab)  \n",
    "\n",
    "    # plotting fitted model\n",
    "    dx = max(x) - min(x)\n",
    "    xp = np.linspace(min(x) - 0.1*dx, max(x) + 0.1*dx, 100)\n",
    "    yp = lg.model(xp, pars)\n",
    "\n",
    "    ax.plot(xp, yp, label = \"logistic reg\")\n",
    "\n",
    "    for method, c in zip([\"normal\", \"delta\"], [\"red\", \"green\"]):\n",
    "        \n",
    "        print(f\"model CI:scale:{scale},method:{method}\")\n",
    "\n",
    "        # define quantile model\n",
    "        fname = f\"lg.get_model_quantiles_{method}\"\n",
    "        \n",
    "        # get result of model quantiles at probs\n",
    "        # pars and cov_pars are obtained via MLE method\n",
    "        res = eval(fname)(xp, probs, pars, cov_pars)\n",
    "\n",
    "        # make plot of quantiles\n",
    "        ax.fill_between(xp, *res, color = c, alpha = 0.5, label = f\"CI:{method}\")\n",
    "\n",
    "    # plot quantiles of models at bootstapped parameters\n",
    "    for method, c in zip(\n",
    "        [\"normal\", \"nonparam_boots\", \"nonparam_stratified_boots\", \"parametric_boots\"],\n",
    "        [\"orange\", \"pink\", \"purple\", \"cyan\"]):\n",
    "    \n",
    "        print(f\"model CI:scale:{scale},boots-method:{method}\")\n",
    "\n",
    "        # generate boostrapped parameters\n",
    "        m = 10000\n",
    "        fname = f\"lg.get_{method}_pars\"\n",
    "        bpars = eval(fname)(x, y, m)\n",
    "    \n",
    "        boots_pars_results[(scale, method)] = bpars\n",
    "\n",
    "        # quantiles\n",
    "        quant = np.quantile([lg.model(xp, p) for p in bpars], probs, axis = 0)\n",
    "\n",
    "        # make plot of quantiles\n",
    "        ax.fill_between(xp, *quant, color = c, alpha = 0.5, label = f\"CI:boots, {method}\")\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "fig.legend(handles, labels, bbox_to_anchor=(1.17, 0.5), loc='center right')\n",
    "\n",
    "plt.savefig(os.path.join(results_path, \"logit_fit.pdf\"))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce16edf",
   "metadata": {},
   "source": [
    "## Parameter CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18dce684",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pars CI:Wald\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'parss' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[1], line 4\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpars CI:Wald\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m      3\u001b[0m \u001b[38;5;66;03m# nr of parameters\u001b[39;00m\n\u001b[0;32m----> 4\u001b[0m d \u001b[38;5;241m=\u001b[39m \u001b[43mparss\u001b[49m\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m]\n\u001b[1;32m      6\u001b[0m \u001b[38;5;66;03m# CI of params:\u001b[39;00m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;66;03m#   - using Wald approximation \u001b[39;00m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m#   - assuming asymptotic MLE distr of parameters \u001b[39;00m\n\u001b[1;32m      9\u001b[0m \u001b[38;5;66;03m#   - cov = hessian(nllf)^{-1}\u001b[39;00m\n\u001b[1;32m     10\u001b[0m \u001b[38;5;66;03m#   - this could be considered as delta method\u001b[39;00m\n\u001b[1;32m     12\u001b[0m df_pars_CI_wald \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame({\n\u001b[1;32m     13\u001b[0m     \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmethod\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mWald [MLE asymp. SE]\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m     14\u001b[0m     \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mscale\u001b[39m\u001b[38;5;124m\"\u001b[39m : [ scale \u001b[38;5;28;01mfor\u001b[39;00m scale \u001b[38;5;129;01min\u001b[39;00m scales \u001b[38;5;28;01mfor\u001b[39;00m _ \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(d)],\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     18\u001b[0m     \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mUCL\u001b[39m\u001b[38;5;124m\"\u001b[39m: pars_CIs[:,\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mflatten()\n\u001b[1;32m     19\u001b[0m })\n",
      "\u001b[0;31mNameError\u001b[0m: name 'parss' is not defined"
     ]
    }
   ],
   "source": [
    "# nr of parameters\n",
    "d = parss.shape[-1]\n",
    "\n",
    "# CI of params:\n",
    "#   - using Wald approximation \n",
    "#   - assuming asymptotic MLE distr of parameters \n",
    "#   - cov = hessian(nllf)^{-1}\n",
    "#   - this could be considered as delta method\n",
    "\n",
    "print(f\"pars CI:scale:{scales},method:Wald/delta\")\n",
    "lst_pars_CI = [df_pars_CI_wald]\n",
    "\n",
    "for scale, x, pars, cov_pars in zip(scales, xs, parss, cov_parss):\n",
    "\n",
    "    # plot quantiles of models at bootstapped parameters\n",
    "    for method in [\"normal\", \"nonparam_boots\", \"nonparam_stratified_boots\", \"parametric_boots\"]:\n",
    "    \n",
    "        print(f\"model CI:scale:{scale},boots-method:{method}\")\n",
    "\n",
    "        # get boostrapped parameters\n",
    "        bpars = boots_pars_results[(scale, method)]\n",
    "    \n",
    "        # quantiles\n",
    "        quant = np.quantile(bpars, probs, axis = 0)\n",
    "\n",
    "        df_tmp = pd.DataFrame({\n",
    "            \"method\": f\"{method}\",\n",
    "            \"scale\" : scale,\n",
    "            \"coef\" : [f\"p{i}\" for i in range(d)],\n",
    "            \"value\" : pars, \n",
    "            \"LCL\" : quant[0], \n",
    "            \"UCL\": quant[1]\n",
    "        })\n",
    "        \n",
    "        lst_pars_CI.append(df_tmp)\n",
    "\n",
    "df_pars_CI = pd.concat(lst_pars_CI, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa12e13",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base (3.12.3)",
   "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": 5
}