notebooks/Repeated measures ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Repeated measures ANOVA    \n",
    "```bash\n",
    "by Jeffrey Chijioke-Uche, MSIS, MSIT, CPDS (IBM Sr. Solution Architect, Hybrid Cloud & Multicloud)\n",
    "IBM Information Technology PhD Scholar at Harvard University & Walden University \n",
    "``` \n",
    "\n",
    "\n",
    "\n",
    "## Contents:    \n",
    "1. [Theory](#theory)   \n",
    "2. [Numerical example](#example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='theory'></a>\n",
    "\n",
    "Repeated measures ANOVA is the equivalent of one-way ANOVA, but for related, dependent groups, and its an extension of the [paired t-test](https://github.com/trangel/stats-with-python/blob/master/notebooks/paired t-test.ipynb) to two or more repeated samples.    \n",
    "\n",
    "## Usage\n",
    "Repeated measures ANOVA is used in two kinds of experiments    \n",
    "1. to analyze changes or consistency among measurements    \n",
    "2. to analyze differences in mean scores under 2 or more different conditions/times    \n",
    "\n",
    "### Assumptions:\n",
    "1. The dependent variable must be continuous (interval/ratio)  \n",
    "2. The observations are independent of one another.   \n",
    "3. The dependent variable should be normally or approximately normally distributed.   \n",
    "4. The dependent variable should not contain outliers.   \n",
    "\n",
    "### Hypothesis\n",
    "Null hypothesis   \n",
    "$\\mathcal{H}_0$ : $\\mu_1 = \\mu_2 = \\mu_3 = \\dots  = \\mu_k$   \n",
    "for $k$ groups   \n",
    "\n",
    "Alternative hypothesis   \n",
    "$\\mathcal{H}_1$ :  at least two means are significantly different  \n",
    "\n",
    "### Test method    \n",
    "\n",
    "Steps     \n",
    "\n",
    "$k = $ number of groups (times or conditions)    \n",
    "$n = $ number of samples (must be equal for all groups)   \n",
    "\n",
    "#### Sum of squares    \n",
    "\n",
    "#### Sum of squares for groups $SS_\\textrm{groups}$ \n",
    "A group is a time or a condition   \n",
    "\n",
    "$SS_\\textrm{groups} = \\sum_{i=1}^k n_i(\\overline{x}_i - \\overline{x})^2$\n",
    "\n",
    "$n_i = $ number of subjects under each $i^\\textrm{th}$ group  \n",
    "$\\overline{x}_i = $ mean score for each ($i^\\textrm{th}$ group  \n",
    "$\\overline{x} = $ grand mean\n",
    "\n",
    "##### Sum of squares within groups $SS_\\textrm{w}$\n",
    "This is calculated in the same way as in an independent ANOVA    \n",
    "\n",
    "$SS_\\textrm{w} = \\sum_1 ( x_{i1} - \\overline{x}_1)^2  + \\sum_2 ( x_{i2} - \\overline{x}_2)^2  + \\dots + \\sum_k ( x_{i2k} - \\overline{x}_k)^2 $\n",
    "\n",
    "\n",
    "##### Sum of squares for subjects $SS_\\textrm{subjects}$      \n",
    "We treat each subject as a level of an independent factor called subjects,   \n",
    "$SS_\\textrm{subjects} = k \\sum(\\overline{x}_i - \\overline{x})^2$   \n",
    "$\\overline{x}_i =$ mean of subject $i$    \n",
    "$\\overline{x} = $ grand mean\n",
    "\n",
    "##### Sum of squares error $SS_\\textrm{error}$    \n",
    "\n",
    "$SS_\\textrm{error} = SS_\\textrm{w} - SS_\\textrm{subjects}$\n",
    "\n",
    "\n",
    "#### Mean sum of squares    \n",
    "These are obtained by the sum of squares divided by the corresponding degrees of freedom   \n",
    "\n",
    "##### Mean sum of squares for groups $MS_\\textrm{groups}$    \n",
    "$MS_\\textrm{groups} = \\frac{\\displaystyle SS_\\textrm{groups}}{k-1}$\n",
    "\n",
    "##### Mean sum of squares for error $MS_\\textrm{error}$   \n",
    "$MS_\\textrm{error} = \\frac{\\displaystyle SS_\\textrm{error}}{(n-1)(k-1)}$\n",
    "\n",
    "#### $F$-statistic:    \n",
    "As in other tests, we calculate an $F$-statistic   \n",
    "\n",
    "$F = \\frac{\\displaystyle MS_\\textrm{groups}}{\\displaystyle MS_\\textrm{error}}$ \n",
    "\n",
    "The above will become more evindent with the following numerical example      \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='example'></a>\n",
    "## Numerical example    \n",
    "In this test, let's use a hypothethical example of blood glucose (BG) readings of persons with diabetes.    \n",
    "The test is done three times, say before, within and after a given clinical treatment and we want to know if there is a significant difference within the groups (times).     \n",
    "\n",
    "The database is not real but just to be used as a toy example.   \n",
    "This contains BG readings (in mg/dL) for 40 persons with diabetes at different times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>before</th>\n",
       "      <th>during</th>\n",
       "      <th>after</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Subject</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>89.162573</td>\n",
       "      <td>94.023517</td>\n",
       "      <td>94.594145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>90.857629</td>\n",
       "      <td>95.273755</td>\n",
       "      <td>95.040646</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>94.912999</td>\n",
       "      <td>96.612870</td>\n",
       "      <td>95.200472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>95.254064</td>\n",
       "      <td>96.818673</td>\n",
       "      <td>97.205801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>97.136291</td>\n",
       "      <td>97.760342</td>\n",
       "      <td>98.428840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>99.809999</td>\n",
       "      <td>99.169227</td>\n",
       "      <td>98.867769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>101.094087</td>\n",
       "      <td>99.579283</td>\n",
       "      <td>99.790581</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>101.531428</td>\n",
       "      <td>99.661758</td>\n",
       "      <td>100.669928</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>101.981148</td>\n",
       "      <td>100.812359</td>\n",
       "      <td>101.751155</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>101.993065</td>\n",
       "      <td>102.274035</td>\n",
       "      <td>101.751638</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             before      during       after\n",
       "Subject                                    \n",
       "0         89.162573   94.023517   94.594145\n",
       "1         90.857629   95.273755   95.040646\n",
       "2         94.912999   96.612870   95.200472\n",
       "3         95.254064   96.818673   97.205801\n",
       "4         97.136291   97.760342   98.428840\n",
       "5         99.809999   99.169227   98.867769\n",
       "6        101.094087   99.579283   99.790581\n",
       "7        101.531428   99.661758  100.669928\n",
       "8        101.981148  100.812359  101.751155\n",
       "9        101.993065  102.274035  101.751638"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "# This is to print in markdown style\n",
    "from IPython.display import Markdown, display\n",
    "def printmd(string):\n",
    "    display(Markdown(string))\n",
    "\n",
    "# Read database locally\n",
    "#df = pd.read_csv('../data/BG-db.csv',index_col=0)\n",
    "\n",
    "# Read dataset from url:\n",
    "import io\n",
    "import requests\n",
    "url=\"https://github.com/schijioke-uche/data-analysis-with-python-an-spss-model/blob/3b9249b72e078eca78a0c38c5aa215a67f21656f/data/BG-db.csv\"\n",
    "s=requests.get(url).content\n",
    "df=pd.read_csv(io.StringIO(s.decode('utf-8')),index_col=0)\n",
    "\n",
    "\n",
    "df.columns=['before','during','after']\n",
    "df.index.name='Subject'\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAFjCAYAAACgxwiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lOXV+PHvM3smM9kmK2FJSGQRZFEEapUqqHWrbxU1oGJtrb72V6xWrWitSK0Frfa1KhVtlarggrj01boD9oWiIKAIyJoQQkLIvs5MMtvz/P4YEogkmWwz2c7nuriuZJ7tJJrkzH2f+9yKpmkaQgghhBi0dL0dgBBCCCF6lyQDQgghxCAnyYAQQggxyEkyIIQQQgxykgwIIYQQg5wkA0IIIcQgJ8mAEEIIMchJMiCEEEIMcpIMCCGEEIOcJANCCCHEICfJgBBCCDHIGXo7gHBobGxk165dJCUlodfrezscIYQQIqwCgQDl5eWMHz8ei8XS6esHZDKwa9currvuut4OQwghhIioV155hSlTpnT6ugGZDCQlJQHBb0pqamovRyOEEEKEV0lJCdddd13z37/OGpDJQNPUQGpqKkOHDu3laIQQQojI6OrUeNgKCFVVZeHCheTk5DBv3jwKCgpaHF+3bh2zZ88mJyeHN954o/n15557jpycHK688kpWr14drvCEEEIIcUzYRgbWrFmD1+tl1apVbN++nUceeYRly5YB4PP5WLJkCW+++SZRUVHMnTuXmTNnkpeXx9dff81rr71GQ0MDy5cvD1d4QgghhDgmbMnAtm3bOOeccwCYNGkSu3btaj6Wl5fH8OHDiY2NBeCMM85gy5Yt7N69m1GjRvHLX/4Sp9PJPffcE67whBBCCHFM2JIBp9OJzWZr/lyv1+P3+zEYDDidTux2e/Ox6OhonE4n1dXVFBcX8+yzz1JUVMQvfvELPvroIxRFafM5Tz/9NEuXLg3XlyGEEEIMeGFLBmw2Gy6Xq/lzVVUxGAytHnO5XNjtduLi4hg5ciQmk4mRI0diNpupqqrC4XC0+ZzbbruN2267rcVrRUVFzJo1q4e/IiGEEGJgClsB4emnn8769esB2L59O6NGjWo+lpWVRUFBATU1NXi9XrZu3crkyZM544wz2LBhA5qmUVpaSkNDA3FxceEKUQghhBCEcWTgggsuYOPGjcyZMwdN01i8eDHvvfcebrebnJwc7r33Xm666SY0TWP27NmkpKSQkpLCli1buOqqq9A0jYULF0oHQSGEECLMFE3TtN4Ooqc1TROsXbtW+gwIIYQY8Lr7d082KhJCCCEGOUkGhBBCiEFOkgEhhBBikBuQexMIITqv5ONPIvas1B9eGLFnCSFCk5EBIYQQYpCTZEAIIYQY5CQZEEIIIQY5SQaEEEKIQU4KCIUQAKwv8vfo/WYMbf/Xy+bNm7njjjvIzs5G0zS8Xi+LFi3i1FNPbfOaO++8k4KCAv70pz+RlZXV7RhXrVrFlVdeidFo7PI99u3bR11dHWeeeWa34zmRx+Ph3Xff5eqrr+7R+wrRGhkZEEL0munTp7NixQpWrlzJr371K5588sl2z//888956623eiQRAHjuuedQVbVb9/jkk0/Izc3tkXhOVF5ezurVq3v8vkK0RkYGhBB9Ql1dHQkJCUDw3fbDDz8MQFxcHIsXL+bPf/4zTqeTX/ziFzz11FPcd999FBUVEQgE+OlPf8oll1zCvHnzSEhIoLa2lr/97W8sWrSIgoICVFXljjvuYNq0ac3PW716NeXl5fz617/mJz/5CY8//jhGo5FrrrmGIUOG8MQTT6DX6xk2bBgPPfQQHo+H+++/n/r6esrKyrj22muZNWsW77zzDkajkXHjxnH//fczZcoU9u3bx8iRI3E4HGzduhWTycTf/vY3Ghsbuf/++6murgbgd7/7HaNHj+bCCy/k9NNPJz8/H4fDwdNPP82zzz5Lbm4uS5cuZf78+ZH/DyIGFUkGhBC9ZtOmTcybNw+v18vevXv561//CsADDzzA4sWLyc7OZvXq1Tz//PMsWrSITz/9lGXLlrFy5UoSEhJ4/PHHcTqdXHnllUyfPh2Ayy67jAsuuIBXX32V+Ph4Fi9eTHV1Nddffz3vv/9+87Ovvvpqli1bxhNPPMH27dvxeDysXr0aTdO46KKLePXVV3E4HPzlL3/hnXfeYdy4cVx66aVceOGFlJaWMm/ePK699lquuOIKEhMTmTBhAi6Xi8suu4wHH3yQiy66iPvuu49f//rXXH/99eTm5vKvf/2L6dOnc+2113Lo0CHuu+8+XnvtNQoLC3nppZdIS0tjzpw57Ny5k1tvvZX9+/dLIiAiQpIBIUSvmT59Ok888QQABw8eZM6cOaxfv568vDx+//vfA+Dz+cjIyGhxXV5eHmeddRYANpuNrKwsCgsLAcjMzARg//79bNu2jR07dgDg9/upqqpqHn34rqbrqqqqKCsr44477gCgsbGRs846ix/84Ae89NJLfPLJJ9hsNvz+1mssxo0bB0BMTEzzdEZMTAwej4f9+/ezadMmPvzwQwBqa2sBiI+PJy0tDYC0tDQ8Hk+Hv4dC9ARJBoQQfUJiYmLzx5mZmTz66KMMGTKEbdu2UV5e3uLcrKwstm7dygUXXIDT6WT//v3NO7UpigLAyJEjSU1N5dZbb6WxsZFly5YRFxfX4j6KojTXDOh0wRKq+Ph4UlNTeeaZZ7Db7axduxar1cry5cuZNGkS1157LZs2beL//u//TrrHic9vzciRI7n88sv50Y9+RGVlZXNNQGvX6HS6btczCNFRkgwIIXpN0zSBTqfD5XJx7733YrFYWLRoEQsWLMDv96MoCn/84x9bXHfNNdfwwAMPMHfuXDweD/Pnz8fhcLQ4Z86cOfzud7/j+uuvx+l0cu211zb/wW8yZcoUbrnlFn75y182v6bT6bj//vu55ZZb0DSN6Oho/vSnP6EoCg8//DAffPABdrsdvV6P1+tl/PjxHV7dcOutt3L//ffzxhtv4HQ6250CcDgc+Hw+HnvsMX7zm9905NspRJcpmqZpvR1ET+vuvs5CDFS+ujqOvPO/lH32f5gTHdhHncLQq2ZjSoiXvQmE6Me6+3dPRgaEGCTq9+3n2wcfItDQgMFmw5V/COeBXKq+3MLYB+7v7fCEEL1I+gwIMQh4KivZs+RRAh4PGT+7kTP/8Xemv76S4ddfi6e8gp333o+3qqq3wxRC9BJJBoQY4FSfj71L/oSvuoaMG28g/b9+hM5kQmc0Muzq2Zxy+3wCbjdlaz9DCwR6O1whRC+QZECIAa70kzU4D+SSdO4PGHL5ZScdT555Hskzz8VTXkHN199EPkAhRK+TZECIASzg8VC4+k10FgsZP/1Jm8veMm/6GfpoK1Vbt+GtqYlwlEKI3ibJgBADWMkHH+GrrmHIZZdgiott8zyDLZrE758FqkrNtq8jGKEQoi+Q1QRCDFCBxkaK3noHfbSV9Cv+K+T53wzTUTUhgYC7kITaPegtlm49//vmzHaPf3fXQr/fzw033MAll1zSofvv2bOHtWvXdqtd76effsqECRNISUnp8j2Ki4vZu3cvM2fO7PI92rJy5Uquv/76Hr+vEN8lIwNCDFCVn2/CX19P2iUXY7DZOnSNdfhw0DQaDheGObqgE3ctfOGFF3j++efZs2dPh64dO3Zst/v2v/zyyzidzm7dY9OmTXz11Vfdukdbli1bFpb7CvFdMjIgxABVunYdACnnd/wdqyU5CfehAhpKSrBmjEBnMoUrvJNER0eTk5PDRx99RF1dHa+//nrzvgXf//732bhxI/feey81NTXU1NRw00038cEHH/DEE0+0uuufz+fjnnvuoaysjLS0NLZs2cJ//vOf5uf9+9//Zs+ePSxYsIDHHnuMX/3qV8TFxTFjxgxmzJhx0q6JVquVhQsXUlJSQllZGTNnzuRXv/pV826EkydP5sUXX2T06NEcOHAAq9XKlClT+M9//kNdXR3Lly/HarXy4IMPnrST4o9+9COmTp3Kvn37UBSFZ555hpUrV1JbW8uiRYtYtGhRxP47iMFJRgaEGIAajpZQt+tbYk8bjyU1teMXKgpRw9JB1Wg8Whq+ANvgcDiat/dty/Tp03n99deJiYlpfq2wsJDbb7+dVatWUVVVxc6dO1m1ahVDhw7l9ddfZ/78+VRWVra4z7nnnsvYsWN59NFHMRqNlJeX88ILL3DzzTfzwAMP8OCDD7JixQpmzJjB888/z9GjR5k0aRIvvPACb775Jq+//jp6vZ5bbrmFyy67jFmzZgEwYcIEXnrpJbxeLxaLhX/84x9kZ2ezZcsWVq9eTXx8PK+88grPPPMMDz30EAAul4tLL72UlStXkpyczPr16/nFL35BbGysJAIiImRkQIgBqGzdZwAkd2JUoIk5ORlXXj6NJSVYRwzr6dDaVVxcTGorycuJXdObdhc8UWu7/uXl5TFjxgwguLFRW7sVNhk6dCimYyMhre2aGBcXx86dO9m0aRM2mw2v19vqfU7ctTA7O7v546ZdC1vbSRHg1FNPbRG/EJEkyYAQA4ymqpSt/Qy91Yrje9M7fb3OYMCclEhjSSne6mpM8fFhiPJkTqeT1atX8+STT1JXV9e8U+GRI0eat/qF1nf4a+21UaNG8fXXX3P++edz+PDhVkccFEVpTjRO3MSotV0T3377bex2Ow899BAFBQW88cYbaJrWqd0F29tJsbWvYQBuHSP6KEkGhBhgnAdy8VZWkjzzPPRmc5fuYUlLo7GklMajJWFNBk7ctTAQCHDbbbcxcuRI/H4/drudq6++mqysrC5tvHLVVVdx7733ct111zFkyBDMrXwvJk+ezD333MMf/vCHFq+3tmtiVlYWd911F9u3b8dkMjFixAjKysoYNWoUy5Ytax4RaE9HdlI8UVZWFnfffTePP/54p79+ITpDdi0UYoA59PJKjrz1DmN+uwDHtKkdvu7EXQs1TaPw9dX4amvJuPEG9JauJRVticSuhV999RVut5uzzz6bQ4cO8fOf/5w1a9aE/blC9AbZtVAI0ULV5i3oTCbiJk3s8j0URcE+6hSqNn+JKz+fmLFjejDCyBg2bBh33nknS5cuxe/3s3Dhwt4OSYg+S5IBIQaQhuJiGoqKSJh6ZpenCJrYsrOo2vwlzryD/TIZSEpKYsWKFb0dhhD9giwtFGIAqfpyKwAJ087s9r2MsTGYk5JoKDpCoLGx2/cTQvRdkgwIMYBUbf4SFIX4KVN65H7R2SNBVXEdzO+R+wkh+iZJBoQYIPwuF3V792EfdUq7mxJ1hi07CwBnbl6P3E8I0TdJMiDEAFG7azeoKrETJ/TYPY12O+aUZBqOFMtUgRADmBQQCjFA1B7rahfXxWSg7tvdrb6ut1pB06j8fBOWlOQO3y9m3KkdPvfvf/87L730EmvXrsVsNvPNN99w9913c9FFFzFjxgzsdjtjxnSuiNHj8fDuu+9y9dVXd+q67+qJnQ1bE87dDoXoLBkZEGKAqN2xE53JhH30qB69r9nhAMDznd7+Pendd9/lkksu4f333wdgw4YN3HDDDdx111289dZblJWVdfqe5eXlrF69utux9cTOhq0J526HQnSWjAwIMQB4q6txHy4kbtJEdEZjj95bH21FZ7HgrapGUzUU3cltc7tj8+bNDB8+nDlz5vCb3/yG7Oxs3n77bYxGIz6fjw0bNvDtt9+SnZ3NN998w4svvohOp+OMM87g7rvv5umnn+brr7/G7XY3dwoEePbZZ8nNzWXp0qVomtbinM8//5x//etfKIrCJZdcwg033MD+/ft55JFHCAQCVFdXs2jRIurq6lrsbHjPPfeQlpZGUVERl156KQcOHGD37t2ce+653Hnnnezbt++k3Q53797N3//+d4xGI0VFRVxyySXccsstLXY7bNrkSIjeIsmAEANA7Y5dAD1aL9BEITg60HDkCL7aWkzxcT16/9WrV3P11VczcuRITCYTmqZxxRVXkJiYyNy5c9m/fz+XXHIJVquVp59+mrfeeouoqCh+85vfsHHjRiDY8/93v/tdi/veeuut7N+/n/nz5/P00083n5Obm8sHH3zAq6++CsBPf/pTzj77bHJzc1mwYAGjR4/mvffe4+233+bhhx9m7NixLFq0CKPRSGFhIcuXL6exsZFZs2axfv16oqKiOO+887jzzjt54IEHWLx4MdnZ2axevZrnn3+es846i+LiYt599128Xi/nnHMOv/jFL7jllls4ePCgJAKiT5BkQIgBoKapXmDCaWG5v8mRQMORI3gqK3o0GaitrWX9+vVUVVWxYsUKnE4nK1euZPjw4Sede/jwYaqqqrjllluA4La/hw8fBlrfyfC7ms7Zv38/xcXF3Hjjjc0xFBQUkJyczDPPPIPFYsHlcmGz2U66x7Bhw7Db7ZhMJhITE0/aZKi13Q4huGmSwWDAYDBgsVg68R0SIjLClgyoqsqiRYvYt28fJpOJhx9+mBEjRjQfX7duHX/9618xGAzMnj2ba665BoArrrii+Ydw6NChLFmyJFwhCjFg1O7chcFmIzozIyz3N8bGouj1eCuq0LKDowU94d1332X27NksWLAAgIaGBmbNmoXNZiMxMRE4vrPg0KFDSUtLY/ny5RiNRt5++23Gjh3LmjVrWt3s57u7CTadM3LkSLKzs3n++edRFIUXX3yR0aNH88tf/pLHH3+crKwsnnrqKY4cOdLi+U0ft6e13Q7buq4zux0KEW5hSwbWrFmD1+tl1apVbN++nUceeYRly5YBwYx5yZIlvPnmm0RFRTF37lxmzpyJ3W5H0zRpISpEJ3gqK/GUlhF/5hQUvT4sz1B0CsaEeLzlFQQaGjBERfXIfVevXs2f/vSn5s+joqK48MILWb16Nffffz8AEydO5PHHH+cvf/kLN954I/PmzSMQCJCens7FF1/c5r0dDgc+n4/HHnusxbvxMWPG8L3vfY+5c+fi9XqbVwpcfvnl3H777cTExJCamtq85XFbOxu2prXdDtsqfjxxt8NLL720Q98vIcIlbLsWLlmyhAkTJjT/T37OOeewYcMGAPbu3ctjjz3GCy+8AMDixYuZPHkyQ4YM4Z577iE9PR2/38+dd97JpEmTOv1s2bVQDCYV/9nIvsf+hxE3XM/Q2Vd0+T4n7lrYmrrdeyj/93oSz/k+saeN7/JzIDK7FgoxmPTZXQudTmeLOTe9Xo/f78dgMOB0OrHb7c3HoqOjcTqdWCwWbrrpJq6++moOHTrEzTffzEcffYTB0HaYTz/9NEuXLg3XlyFErwr1BxqgYkOwiM5XX9eh89uypzy3/ROsfkxAWe4+ilO7N++diiQDQvQlYUsGbDYbLper+XNVVZv/qH/3mMvlwm63k5mZyYgRI1AUhczMTOLi4igvLyctLa3N59x2223cdtttLV5rypCEGAwaS0pAp8OclBTeB0VHodmsKOXVoKrQyjy9EKJ/CttP8+mnn8769esB2L59O6NGHW+EkpWVRUFBATU1NXi9XrZu3crkyZN58803eeSRRwAoLS3F6XSSFO5fcEL0Y6rPh6eiEktyErp2RtB67HmpCSj+AEplbdifJYSInLD99rjgggvYuHEjc+bMQdM0Fi9ezHvvvYfb7SYnJ4d7772Xm266CU3TmD17NikpKVx11VXcd999zJ07F0VRWLx4cbtTBEIMdo2lpaBpWFJTI/I8NSUBfW4RutIqAknxEXmmECL8wvaXVqfT8dBDD7V4rakzGMDMmTNP6sltMpn485//HK6QhBhwGo+WAGBJi0wyoCXFoykKSmkVjM8KfYEQol+QST8h+rHGo6UAWFJ7dhOdNhkNaI5YlKo68Poi80whRNhJMiBEP6VpGp6yMoyxseh7aN1/R6gpCSgQHB0QQgwIkgwI0U/5ampQvV7MndhWuCdoqQkA6CQZEGLAkGRAiH7KUxpsdWtJjnAyEB+DZjQEk4Hw9CwTQkSYJANC9FONx9rcRnpkAEUJThW4G8HZENlnCyHCQpIBIfopT2lZsNlQoiPiz9ZSmqYKKiP+bCFEz5NkQIh+SPX78VRWYk50hG1zonaf35QMlEjdgBADgSQDQvRD3opKUNXITxE0iY5Cs0WhVNRI3YAQA4AkA0L0Q031ApEuHjyRmhSP4vOj1NT3WgxCiJ4hyYAQ/ZCntJeKB0+gHWtHrJRV91oMQoieIcmAEP2Qp7wCncmEMTa212JQk+MA0JVLMiBEfyfJgBD9jOr14qupwZSYiKIovRdIlCVYN1BeE9zSWAjRb0kyIEQ/4ymvAMCcnNjLkRyrG/AHUGqcvR2KEKIbJBkQop9pTgaSkno5EtCSj9UNyFSBEP2aJANC9DOe8mAb4r6QDKjHigh1UkQoRL9m6O0AhBCdc7x4MKa3Q4EoM5rdGuw3oKqg69j7i4++OBTWsE500fcyIvYsIforGRkQoh/pM8WDJzheNyD9BoToryQZEKIf8VQ01Qv0fvFgEy0puMRQKavp5UiEEF0lyYAQ/cjxlQS9Xy/gxscuXSWfpgZXElRUHOFrXTn1eHs5MiFEZ0nNgBD9iKesqXiw90YG/Kh8o6vgK10ZfkWDaJgWoyehvJHVylG+NJQyTU1hgpqIjr4xlSGEaJ+MDAjRj3gqK1GMxl7rPFiPl3cMeXypL8WInrMCaczzjSY2MRWTX+OH5fGY0fOFvoSP9AUEkE2MhOgPJBkQop9Q/X581TWYHY5eKR4sU9y8acilQmlkjBrPtf5RTFQTsWFq7jcwsjRAjv8U0tVoCnT1rNMXokpCIESfJ8mAEP2Et7IKNA1TkiPizy5V3Lynz8dDgHMCQzgvMBQT+ubj6rEiQl1FDVEYuDiQQapqJVdXy1ZdWcTjFUJ0jiQDQvQT3qaVBI7I1gtU0si/9Pn4UJkVGMZ4tZVkxGIO7lNQUQOahhEdlwQysGtGvtKVUaq4IxqzEKJzJBkQop/wVFQCYI7gyIATH+8b8vEqKjMDQzlFi2vz3O/uU2BGz3mBoWjAWn0hfmQzIyH6KkkGhOgnPBUVoNNhSkiIyPN8qHxoOIRL8TM9kMooLb7d87XEY/0GKo73G0jXbExQHdQqXnboKsIarxCi6yQZEKIf0FQVb2UVpvg4FL0+9AXdfR4a6/VHqFAaGavGM0kNPTWhJh6vGzjRFDUFi6bnK105DfjDEq8QonskGRCiH/DV1qL5/ZgSI1Mv8K2uiv26GpLVKM4JDEHpSL+AaAuaxRzcwVA7voLAjJ4pajI+RZViQiH6KEkGhOgHmusFEsNfL3BUcbFRV4xF0/PDwAj0Hf01oSioSbEoHh84G1ocOlVNIEYzsVtXhQtfGKIWQnSHJANC9APNKwnCPDLgxscn+sNowIWB4dgwdup6LfHYlsbfmSrQo2NyIAlV0aR2QIg+SJIBIfoBT3lwZMAUxpEBDY11+iLcip9pairpmq3z92jatKj85E2LRmlxRGkGduuq8KqebscrhOg5kgwI0cdpmoansgKD3Y7ebA7bc3boKijUORmm2jpUMNgaLSYazWhAV1F90jEDOiaoDryKyl7n9u6GK4ToQR1KBsrKgkU/W7du5ZVXXsHtlgYiQkRKwOVGbWgMa73AYX81m3SlRGkGZgaGdqxgsDWKgpYYh+JqhIbGkw6PUx0YNIW99V+jadKmWIi+ImQy8OCDD7Js2TJyc3O56667+Pbbb1mwYEEkYhNCcKy/AGAK006FHs3PS+4tqIrGrMBQrJ2sE/iu5iWGrUwVmNGTrcXhDNRypDG/W88RQvSckMnAzp07WbhwIR9++CFXXXUVixcvpri4OBKxCSEAb5hXErzZ8A3lqouJgUSGafZu36+5bqDi5GQAYJwabJokUwVC9B0hk4FAIICqqqxdu5YZM2bQ0NBAQ0NDqMuEED3EE8aVBFu9hWz2HmaYPo5pakqP3FOLt6PpdSgVta0eT9KicBhTKGzIxeWv75FnCiG6J2Qy8OMf/5izzz6b9PR0Jk6cyJVXXklOTk4kYhNCEOwxoLNY0EdH9+h9KwIuVrm3Y8bAT6xndryfQCg6HVpCLLpaJ3hP7imgoDDaPgkNjVzXtz3zTCFEtxhCnfDTn/6UG264Af2xFqivvPIKCRHqjS7EYBfwePDX1RE1NB1F6WJRX2v31VRecm/Bg5/rrWeQrLdR2WN3PzZVUF6NUlGDNiTppOOZ1jFsrlpDnnsXE2Km9ejXJoTovDaTgXnz5rX7A/ryyy+HJSAhxHHeyvD0F/igcQ8FgWqmGIcx1TS8R+8NwSJCPcHmQ4FWkgGzzsKwqGwONeyj0ldKoim1x2MQQnRcm8nAbbfdFsk4hBCtON6GuOfqBfb7ylnj2Y9DF83V1ok9dt8TaY4YNEVptflQk+zocRxq2Eee61tJBoToZW1OEk6dOrX5n16vJy8vj0mTJqEoClOnTg15Y1VVWbhwITk5OcybN4+CgoIWx9etW8fs2bPJycnhjTfeaHGssrKSH/zgB+Tl5XXxyxJiYOjpNsRO1cPL7q0oKNxonUKU0r1lhG0yGNDi7CjV9eAPtHpKetRIzLoo8ly7UTU1PHEIITokZMXQSy+9xF/+8hdefPFFXC4XCxcu5IUXXgh54zVr1uD1elm1ahV33XUXjzzySPMxn8/HkiVLWL58OStWrGDVqlVUHPul5/P5WLhwIRaLpRtflhADg6eiEsVgwBgX2+17aZrGq+6vqNMaudRyKiMM4a390ZLiUDQNpar1VQV6RU+mdQyNqpsST2FYYxFCtC9kMvDOO+/wwgsvEBUVRXx8PG+++SZvvfVWyBtv27aNc845B4BJkyaxa9eu5mN5eXkMHz6c2NhYTCYTZ5xxBlu2bAHg0UcfZc6cOSQnJ3f1axJiQFB9PrxV1ZgSElB03a/03+A9yC5/CaMNScwyn9IDEbavqflQe1MFGdbRABxy7wt7PEKItoX8DaPT6TCZTM2fm83m5pUF7XE6ndhsxzc60ev1+P3+5mN2+/HmJtHR0TidTt5++20SEhKakwghBjN3YRGoao80GzoSqOWfDbuIVkxcb52CLgLV+1pTJ8I2mg8BpJqHYdZFUeDeL+2JhehFIZcWTp0IZBuYAAAgAElEQVQ6lUcffZSGhgbWrFnDqlWrmD59esgb22w2XC5X8+eqqmIwGFo95nK5sNvtrFixAkVR+OKLL9izZw8LFixg2bJlJCWdXI3c5Omnn2bp0qUh4xGiv3HlB9v1mrpZL+DTArzs2ooflZusZxCri9AUnNmIGhONUlkLqgqtjG7oFB0jok5hv2sHpZ4jpFqGRiY2IUQLIUcG7rnnHkaMGMHo0aP55z//yQ9+8IMO7U1w+umns379egC2b9/OqFGjmo9lZWVRUFBATU0NXq+XrVu3MnnyZF555RVWrlzJihUrGDt2LI8++mi7iQAEVz3s27evxb+1a9eGjE+Ivs6Vfwjofhvi9xt3c1St42xTJuOMka3a1xLjUAJqsJCwDc1TBQ17IxWWEOI7Qo4MLFmyhMsvv5w5c+Z06sYXXHABGzduZM6cOWiaxuLFi3nvvfdwu93k5ORw7733ctNNN6FpGrNnzyYlpWdaoQoxULgO5oOiYHJ0vdDvgL+Czzy5JOmi+a+o8T0YXcdoiXFw8Eiw+ZCj9SLINMsITIqZAvd+psXNkgZEQvSCkMlARkYGixcvpra2lssuu4zLL7+coUNDD+XpdDoeeuihFq9lZWU1fzxz5kxmzpzZ5vUrVqwI+QwhBipN03DlH8IYG4vO2LXlfw2aj1dc2wC43joFsxLyx73HqYnBBEBXWUtbiwf1ip7h1lPIde2i3HuUZPOQyAUohAA6ME1w3XXX8dprr/H8889jNpv55S9/ydy5cyMRmxCDlqesjIDb3a0pgrcbdlClubnQPJrMMC8jbJPVgmYxBzctaqdAMCMqOI1YIKsKhOgVHVqvVF9fz+eff87GjRsJBAKcffbZ4Y5LiEHNdfBY8WBS14oHd3iL2ew9zFB9HBdZxvRkaJ2jKGiJsSgeL7ga2zxtSFQmBsVEvnufrCoQoheEHDe89dZb2b17NxdeeCG33347EyeGp32pEOI457FkwOzo/MhAndrIaw1fY0DHPOsZ6JUe2o2wi1RHLLqiMnSVNai2qFbPMSgGhkVlke/eQ5WvDIdJaoiEiKSQycA111zDjBkzmpcFCiHCr2klQVeWFb7RsB2X5uUKy2mk6WN6OLLO047VDSgVtTAirc3zMq2jyXfv4ZB7nyQDQkRYyLcMp59+OosWLeKGG26gurqa++67j9ra1tuLCiF6hutgPqaEBAzW1t9Jt+Ub7xF2+I6SpXfwA3NW6AsiQIuzo+l1wX4D7Ui3ZKJDT2GD7EkiRKSFTAYWLlzIaaedRk1NDdHR0SQnJ/Ob3/wmErEJMSj56urwVlYSnZnRqevcqpfVDTswoGOOdXJEugx2iE6HlhCDUusEr6/N04w6E2mW4VT5ynD66yIYoBAiZDJQVFRETk5Oc1viX//615SUlEQiNiEGpaYpguiRmZ267r3Gb6nTGvmhZTQpenvoCyJIc8SiAEpV+3/kh0VlA1DYkBuBqIQQTUImA3q9nvr6+uZGIIcOHULXA5umCCFa17SSIDqz48lAnr+Cjd5DpOlimGUeFfqCCFM7sE8BwPCo4NSGTBUIEVkhqwJvu+025s2bx9GjR/l//+//sX37dhYvXhyJ2IQYlI6PDGTw9uuhm28F0Ph30hEwQHZpNP9XtKNLzx2SZAt9Uhc1dR8MVTdgM8QSb0ziaGMBPtWLUWdq93whRM8ImQzMmDGD8ePHs2PHDgKBAA899BCJ3dw4RQjRNld+PvqoKCwdbNF90FaL0+gj02UnwRehTYg6y2REs1tRKuuCmxa1Y1hUNjt8X1DceIgR1r43yiHEQNRmMtDWToB79uwBYP78+eGJSIhBLODx4C46QsyY0SgdmI5r1PnZb6vBFNAxpi4+AhF2nZoYhz6/OFhI2I7hUVnsqPuCww25kgwIESEy+S9EH+I+XAiq2uGVBHvs1fh1GmPq4zFp+rDG1l0dnSpIMg3BorNS1HBQuhEKESFtjgzIO38hIs+Vf6x4sAMrCWqMHg5bncT4jIxw963VA61pKiJUKtpPBhRFYVjUSA64dlHhPUqSbFwkRNjJyIAQfUhHVxJoaOyMqQQFxtc60NFHegq0xxaFZjKGXFEAx5cYHpYlhkJEhCQDQvQhroOHUPR6rMOHtXtescVFldlDWoOVJG/nuhT2mqZNixo8KHXV7Z4q3QiFiCxJBoToI7RAAFdBAVHDhqIzGts8T0VjT0w1igbj6nppa+IuUh3BqQL9kUPtnifdCIWIrDaTAZ/Px2OPPcbq1asBOPvssxkzZgzjxo3j22+/jViAQgwWDUdLUBsbQ04RFFqduAx+RrjtRAfaThr6oqZNi/RFBSHPHSYNiISImDaTgf/5n/+htLSU888/H4DExET27t3Lk08+yd///veIBSjEYNHUbMjWTvFgAI19tmp0msKo+rgIRdZztHg7mk5BfyQ/5LnDpTWxEBHTZjKwbt06HnnkEeLjW65dPv/88zlw4EDYAxNisGleSZCZ0eY5h631NBgCZLrsRKn9cFtxvR4tPgZdaTF4Pe2e+t1uhEKI8GkzGTAajRgMx3/ZPPjgg80fm0zSIlSInnZ8JUFGq8cDqOyz16BXFU5x9r9RgSaaIxZFU9EfLQx57rCoLAIEKG4MPa0ghOi6NpMBvV5PRUVF8+eTJ08GoLS0FL2+bzc3EaI/cuUfwpycjMHW+h4B+dH1ePQBRrpiMKv992dQba4bkKkCIfqKNpOBnJwc5s+fz8GDB5tfKygo4Ne//jXXXXddRIITYrDwVlfjq6lpc1TAr6gcsNVgUBWynbERja2nNXUiDLWiACDRlIZFZ6WwIU+6EQoRRm1OOl577bXU1dVxzTXXYDQaURQFr9fLLbfcwhVXXBHJGIUY8Jx5waS7rc6Dh6OcePUqo+rj+nzb4ZAsZtQ4B/ojBaCpoLS9wlmn6E7oRlhCkjktgoEKMXi0W4F066238rOf/Yzc3OAQ3ciRI7FY+uiuaEL0Y65jyYAta+RJxzQ08my16DSFka6YSIcWFoGhmRh3bUVXUYaalNruucOisjng2kVhQ64kA0KEScimQyaTiVNPPRWPx8Nrr73Gl19+GYm4hBhUnHnBtfS2rKyTjh21uHEb/Ax32/p1rcCJAkNGAKAvDl0YmG7JQIdeWhMLEUZtjgz8+9//5re//S2JiYnceOON/PnPf2by5MmsWLGCnJwc/vu//zuScQoxoDlz8zAlJGBKaLmUV0PjgK0WNMjqp7UCmgZ1Xj2lbiM1HgO1XgN6j56LgMJtO/nKlUxMtIlYm5n0JBvxdjOKcnyvBaPOTJplGEcaD+H012EzDIzRESH6kjaTgSeeeILly5dTV1fHz372Mz7++GPS09Opq6vj6quvlmRAiB7iranBW1lF/JlTTjpWaWqkxuQhtcGKrR91G6zx6DlUZ6Gg3swRp5nGQMtBSEVTmaXoia05Sn5xy3bDcXYz2emxTDgliShz8FfUsKhsjjQeorAhj7H2yRH7OoQYLNqtGRgzZgwAw4cPJz09HYCYmBisVmv4IxNikHDmHpsiyD55iiDXFtzu95Q+PirgCSgcrjdTUG+moM5Crff4r5YYk58R9kZSrD4SLD5iTQFMo2ehrPqcpKJ8browixqfQnW9h4KjdRSU1LN1bxnf5FYweVQyE0clMjwqm03VayhsyJVkQIgwaDMZ0OmOZ/JmszkiwQjRX6zJ29Bj9zJ8uQkjkGttYP8J961uqKXU0kCC10yCr28V7qoalLiNFBx793/UZUI7to2ySaeSHdtARkwjw+0e4syBk67/xrOToclRpBVpVB7ZQP2wNBQzZCTCsFOhuBDy8wJ8ubuEHflHOW0SWPU2jjQe4tv6reiVjndfvIiMnvqyhRiw2vyJKi8vZ+nSpSd93PS5EKJn6IrKAFCHpbR4fVfZPqDv1ArUe/Xk15kpqLdwuN6M59jQv4JGqtXHiJhGMuweUqO96JQQNwOcKYkA2EorqB92fJWAXg/DMiBtKBzK0yjIV9i2WSNpQjKa6SA1vkocppQ27iqE6Io2k4E5c+a0+nFrnwshuk5XVIYWEw0x0c2veQM+cisPEeXXk9bYO9NymgalbiN5dRYO1loobzjehtxu9DMqroERMY0Mt3mwGDrfEMjVlAyUtP7mwmCA7NEQl6Dx7Q4oz03CfOpBqrzlkgwI0cPaTAbmz58fyTiEGJzqXCi1TgLjWvYXyK06hE/1k+WOQ6EDb7N7iKZBscvEniorebUWXP7gUka9ojHC3sjI2OC7/zizH6WbYflsVjz2aKJLKoIPbuOGiUkw9Xvw9dZYNJ+JSrWc7GitxYoDIUT39MNtz4QYOJqnCIYmN7+maRp7ynJRUBjhtkckDk9Ax5clNnZWRjcX/0UZApya4GpOAEz6nm8H7ExNxHGgAHOdE09s219rlBVOn6rw1ZFE1Phi9hfUMjqj/27WJERf02Yy4Ha7ZdWAEGGmKywFWtYLlLsqqWyoJiNuGJbi8Obrbr+B3fUJ5LtjCWg6jDqVsfFuTk1wM8zu6dDcf3e4UpJwHCgguqS83WQAwGKBzJRk8r3FlNRXEHs0jlRpSChEj2izA+G8efMAWLRoUaRiEWLQaW1kYHd5sNPeqUnZYXuuV9XxTW0i75dkkuuKx6Lz84P0Gm4ZX8LFGdWMiAl/IgDBkQEAW0lFiDODkqMdKCjo48vYswuc9eGMTojBo92RgbvvvpsNGzbg8XhOOr5kyZKwBibEYKArLG1RPOjxe8mrLiDGbCM9JpUyDvX4M480RLOlOpVG1UCU3sdpMaVkWOsYmtz61snh5E5KQNXpsJV2LBnQKwZijQnUWCvR9I3s/NrCmd8DQ//pxyREn9RmMrB8+XI2b97Mtm3bmDp1aiRjEmJwqHOh1LlaFA8eqMwnoAYYm5Td4wVyPlXhq5pk8t1x6FA5Laac0fZqDErvbQ2s6fW4kxKwllei+P1ohtDTIvHGJGp8lSRmlVG+ZwR7v9UYPykCwQoxgLX5k5eWlsaPf/xjxowZQ1ZWFvn5+QQCAU455RQMHfiBFUK0r7le4NgUgaZp7CnPRafoGOU4effC7nD6jWyoSKfWbybe2Mi0hKPEGb09+oyucqYmYiutILq8Cmdacsjz443J5LMXfXw5sXHDKS1RSCnVSJLVhkJ0Wci/6j6fjx/+8IfExcWhqioVFRX89a9/ZeLEiZGIT4gB67vNhsrdVVQ31jIyfjhRxp7rOFjuieI/lUPwqAZOsVUzObYsIvUAHdXUbyC6pKJDyYBFH4VVb6POX8m4cX62fm5g726ISwCjTBcI0SUhk4E//vGPPPHEE81//Ldv384f/vAH3nzzzbAHJ8RA9t3iwQOV+QCMcmT22DOONlrZUJGOhsKUuBKyj+110Jc4U5OAYPOh0g5ek2BMoSiQh9dUzsjsNPIOKBzYq3HqaeGLU4iBrM3VBE3cbneLUYBJkya1WlAohOgcXWEpWqwNYqJRVZW8qgIsBjNDY3pmvVxTIqAAMxKL+mQiAOC1R+OzWjpcRAjgMAUTqEpfGcMzwR6jcfSIQlVluKIUYmALmQzExsayZs2a5s/XrFlDXFzoZh+qqrJw4UJycnKYN28eBQUFLY6vW7eO2bNnk5OTwxtvvAFAIBDgvvvuY86cOcydO5f9+/d39usRon84VjzYNCpQWHeURr+H7IQRLTYJ66rSRiv/qQjuNHp24hHSLO5u3zNsFAVnSiImpxuj09WhS6x6O2ZdFNXeclACjBkHoHFgb7CZoRCic0L+1vnDH/7Ac889x7Rp05g2bRrPPvssDz30UMgbr1mzBq/Xy6pVq7jrrrt45JFHmo/5fD6WLFnC8uXLWbFiBatWraKiooLPPvsMgNdff5077riDJ554ohtfmhB9l67gKADq8FTg+BTBKT0wRVDrM/GfyiFowDl9PRE45vhUQcdGBxRFIcGUjEqAGl8VMbGQNgSc9QrFReGMVIiBKWTNQEZGBqtXr8btdqOqKjZbx9Yib9u2jXPOOQcITi3s2rWr+VheXh7Dhw8nNja4G9sZZ5zBli1buPjiizn33HMBKC4uJiYmprNfjxD9gq6gBAB1RCpev5eCmiLiLDEkWhO6dV9PQMeGinR8mp7pCcX9IhGAE4oISyuozh7RoWscxhSONhZQ5SslwZRE1igoLdU4eABS0oIbHQkhOqbDPy6dbU3sdDpbJA56vR6/34/BYMDpdGK3H289Gh0djdPpDAZkMLBgwQI+/fRTnnrqqZDPefrpp1tsryxEf6A7XIKmBFcSHKwuJKCpnOLI7FZvgYAGG6vScQZMnGqvJMPaf9rzuZIdaIrS4ZEBALshDqNiospbjmbVMFsURmRCfq7CoYMa2aPCGLAQA0z3JyfbYLPZcLmOz/+pqtrcn+C7x1wuV4vk4NFHH+Xjjz/mgQcewO1u/53Nbbfdxr59+1r8W7t2bQ9/NUL0IFUNFg8mJ0CUuXmKIDsho1u33VmbRJnHylBLPafFdPyPal+gmow0OOKILqtECagduqZpqsCveanzVwMwIhPMFo3CQ+BpDGPAQgwwYUsGTj/9dNavXw8ElyOOGnU8Tc/KyqKgoICamhq8Xi9bt25l8uTJ/POf/+S5554DICoqCkVReqSYSoi+RCmtQvH4UIenUu9xcdRZRpotGbs5usv3LGm0steZgM3gZVrC0W5vL9wbnCmJ6AIBoiqrO3xNgjFYgFnlDS7T1OshYySoqkJBfljCFGJACjlNsHv3bp599llqa2vRTijTffnll9u97oILLmDjxo3MmTMHTdNYvHgx7733Hm63m5ycHO69915uuukmNE1j9uzZpKSkcOGFF3Lfffdx3XXX4ff7+e1vf4vF0nPNV4ToC06sF8irCq6yOcWR0eX7eQJ6NlWloaBxVkIxRl3/LKd3piaR/O0BbCUVuJMdHbom1uhAj54qXykZ2mgURWHIUCg4qHGkMDhSIIQILWQysGDBAnJycjjllFM6NZ+p0+lOWnWQlZXV/PHMmTOZOXNmi+NWq5Unn3yyw88Qorfov9jZ9Wu37gFAqXGSX1iKToGs3Hr0uZ2/p6bBl9UpNKoGJsaUk2Dqeg+Q4nJnl6/ttFZaB7tSm4oIy4HRHbqNTtERb0qiwluCK1CPzRCDTgcZWbD322DtgBAitJDJgMVi4frrr49ELEIMCkpVHZpeR22siXJdA8NUG5aO1/K2UNhg50ijnWSzmzH2qh6ONLIa42Lwm02dKiIESDClUOEtocpbis0QXIGUlg6HDmoUF0FFTQOJcVHhCFmIASPkhPzZZ5/NihUryM/Pp7i4uPmfEKIL/H6UWidavJ2DhjoARqqxXbqVV9XxVU0yelTOjC/pl3UCLSgKrmQHltp6DA0dr/6LNyaioKPKV9b8mk4HmVnB2oG3PjsQjmiFGFBCvh353//9XwD+8Y9/NL+mKIpU7AvRBUp1PQqgJsRyUKlF0SBT61o/je21STSqBibElGM3+Ho20F7iTE0ktvAo0aWV1Gakd+gavWIgzuig2ldOQ8BFlD5YiJk6BA7manz65WHmXjiGmGhTOEMXol8LmQysW7cuEnEIMSgolcHRAHeClVJdFelqNFFdmCIo80Rx0BVHrNHT76cHTnRi86GOJgMACaZkqn3lVHnLSI8KVg3qdDA8Aw7sDfDB5/nMuaBjdQhCDEYhpwmqqqq44447mDZtGlOmTGH+/PlUVPSvNcxC9BW6quBmQQcTg2P6WVrnpwhUDbbVpAAaU+NK+tR2xN3VlAx0ZtMiOGGJ4QlTBQBDhoItysh7Gw7i8QV6JkghBqCQycDChQs57bTTWLt2LevWrWPixIncf//9kYhNiAFHqapDs5jYG+0GDTLVzk8RHHTFUuszk2mtw2EeWJ11/FEWGmPtRJdWdGrHIaPORIwhnnp/DV71+PfEYIBLvp9JncvL2i2HwxGyEANCyGSgsLCQm266CZvNRkxMDDfffLMUEArRFQ2NKA0efAl2SvQNDNGisWLs1C28qo6ddYkYFJUJseVhCrR3uVIcGDxezDWda6fsMAXXK1Z6W44OXHZ2JkaDjnf+nUugg90NhRhsQiYDiqJw9OjR5s+Li4ub2woLITpOqQrWC5Q7ggnAyC4UDu6uc+BRDYy1VxKlH5jD3s6uThUcSwaqvKUtXo+3W5h15nBKKt18ubukZ4IUYoAJ+Vf99ttvJycnh4kTJ6JpGt988w1/+MMfIhGbEAOK7ljxYO6x5nqdXVLo9BvZ74zHqvcx2t7xlr39zYlFhJVjRnb4OrPOgk0fS62/Gp/qxag7vnrgR2dn8tEXh3h3w0G+d9qQng5ZiH4vZDJw3nnnMXHiRHbs2IGqqvz+97/H4ehYq1AhxHFKVS0asCfRT4oaTXQnpwh21TlQUZgQW45BGbid9dxJCag6XadHBiA4VeBsqKXKV0aKeWjz68NTY5g0Kont+8vJL64lc0jXejsIMVC1OU2watUqAJYuXcqrr77Krl272L17N6+99ppsGSxEZ6kqSlU9nhgLHpOOjE5OEdT6TBxyxxBnbGREVP/ZmrgrNL0ed1I8URXVKP7OTYUcrxsoPenY5ecERxne23Cw+0EKMcC0mQxonajkFUK0T6l1oQQClCUGh65HdHIVwc7aREDhtJiK/t9psANcKYnoVBVreed6KFj0Vqx6O7W+Svxqy0ZMZ4xJIS0xmn9/VUSts+t7OAgxELU5TTBnzhwA0tPTueKKK1oce+WVV8IblRADjFJZA8CBJLBrRhIwd/jaSq+FokY7DlMDQyyucIXYpzhTEklhH7bSclxpSZ261mFKobChnmpfOUnm4/UBOp3CZWdn8vd/7uLjTQVcc/6odu4ixODSZjLw4osv4nQ6ef311zly5Ejz64FAgPfee4/rrrsuIgEKMRAoFcFmQ0VJBjLUGBQ6/vY+OCoAE2LKB8WoALQsIuysYDKQS6W3tEUyAHD+mcNZ+eEePvw8n9nnZaPXh1xQJcSg0OZPwogRI1p93WQy8cgjj4QtICEGIl1lLV6znhq7vlP1AhUeCyWeaJLNLlIsDWGMsG/xxNrxWcyd3sEQwKq3EaWLpsZXQUDztzxmMXLu6cOoqG1k296yNu4gxODT5sjAeeedx3nnncfFF19MVlZWJGMSYmBxN6K4GzmabsGEnjTN2uFLv60PrtwZH1MZruj6JkXBlZJIXMERDO4G/NbObUHsMKVQ1HiQat/JycTFZ2Xw4ReH+PCLQ0wdl9oz8QrRz4VcWnjzzTejtDI2KbsWCtExusrgFEFhkp5hmh196F5fAJS4jRxttJFkcpNsHjyjAk2akoHo0gpqM4d16tqEY8nAdxsQAWQOiWX0iHi27S2ltMpNSkLHkzMhBqqQycCKFSuaP/b7/Xz66ad4vd6wBiXEQKJUBIsHi5OMnKraO3zd5pLgueMG26jAMSd2IuxsMhCtt2PWRVHtLccX8GHUt+zpcPH3MthXUM0nmwuYd/HYHotZiP4q5FuU9PT05n8jRozg5z//OWvWrIlEbEIMCEplLQEdlMcbGa51LBkodxvIq43CYWogxewOc4R9kyslOEXSlSJCRVFwmFIIEGBH6d6Tjp89KZ3oKCOfbC7AL/sVCBF6ZGDLli3NH2uaxoEDB/B4ZI2uEB3iD6DU1FPiMJCki8YS6Ni+HlvKjo0K2CsHzQqC7wpYzDTExRBdWhncwbCT3wiHMYXixkNsLvyaM4ac1uKY2ahn1pnDeHf9QTbvKuH7E6VFsRjcQv5meuqpp5o/VhSF+Ph4WU0gRAcpVbUoGhxNMjGig6sI6r169lVH4bD4SBskfQXa4kpJJHHfQSzVdTQmdK6FsM0Qi0kxs6X4G/xqAINO3+L4RdMzeHf9QT78Il+SATHodapmQAjROcqx4sGjiUZOV20duuar8mg0FKYkO1EG+Qh2UzIQXVrR6WRAURQSTMmUeArZXbafCaktawOGpdgZn+XgmwMVHCl3kp7Usf8+QgxEIWsGSkpKmD9/PtOmTeP73/8+d999N1VVnWsRKsRg1dRsqNphwYEl5PmegMLOimiiDQFGxw/OWoETdXU74yZNexVsLvq61eMXfy8DgI++ONSl+wsxUIRMBn77299y1llnsXbtWj7++GPGjx/PfffdF4nYhOjfNA2lsoYamx6HpWNdB3dWRONVdUxOcmKQ5ng0JMah6nVdKiIEiDHEYzfb+PLIN6jqycMs3zstjVibibVbDuP1dW5TJCEGkpC/bqqqqrj22mux2WzYbDZuvPFGSkpKIhGbEP2aUudC5wtQnGRkaAeWFAa04BSBUacyIXFw1wo0Ce5g6MBaUY3O5w99wXcoio4z0ydS21jHvsq8k44bDXrOP3M49W4fG3cU90TIQvRLIZOBCRMm8P777zd//tlnnzF+/PiwBiXEQNBcL5BkZKgWej56f3UUTp+B8Q43FoPsGtrEmZKIomlYy7vWb2Ha0EkAbC7a3urxi45NFXz4+aEu3V+IgaDNAsIxY8agKAqapvHGG2/wu9/9DkVRcLvdxMbG8sc//jGScQrR/xxrNuRyWLFhbPdUTYNtZTYUNE5PckYiun7DlZoI34CtpALnkJROX39a8hisxii+LNrOTyZddVJH1VRHNJNHJfH1/nIKSuoYkdq57aWFGAjaTAb27j25UYcQouO0ymo8RgWbPR5CvNEvdJooazAxKs5NrFnmrk/k7MYOhgAGvYEzhpzGhoIvyasqINuRcdI5P/xeBl/vL+eTTQXc/OPTTr6JEANcm8nAqlWryMnJYenSpa0enz9/ftiCEqLfa2jE6PRwcIiJYYSuF9h2rMnQGckyKvBdXns0vihLl5MBgGlDJ7Oh4Es2F33dajIwbVwqcXYz67YW8pNLT8Vk1J98EyEGsDZrBjRN5iyF6Cpd+bH9CJJNDNGi2z23osFAfp2F9GgPadG+SITXvygKzpREzE43RlfXlltOTD0Vs97E5qKvW5B+RKgAACAASURBVP3dZtDrmDVlGM4GH59LIaEYhNocGZgzZw4AR44cYcmSJRELSIiBQK2oBqAxyY6R9t9lbisLFhfKqEDbXKmJxB8qIrq0gpqRwzt9vdlgYnLaeDYVfUVhbTHD49JPOufC6SN467NcPt5cwLlndG5jJCH6u5CrCfbv34/LJcuchOgMtbwKnx5scY52z3P5dOytthJn9pMV2xih6PofV1PzoZJuTBUMC64q2NRGA6IhiTYmZCeyK6+SorL6Lj9HiP4oZDtiRVE477zzyMzMxGw2N7/+8ssvhzUwIfotjxdLXSOHU42kK/Z2iwe3l0cT0BTOSK4ftBsSdYQr2YFG14sIAU5POw2jzsDmoq+5ZvxlrZ5z0fQMduRW8Mnmw/zsR+O6/Cwh+puQycA999wTiTiEGDCUY0sKS5LMTNCi2jzPF1D4psJGlCHAqQkNkQqvXwqYTTTGxxJdVgmqCrrOt2eMMlqYkDqWbcU7Ka4vZYj95GWK009LxW4NdiScd/EYjAYpJBSDQ8ifqI8//pipU6e2+PfWW29FIjYh+iXfseY4/qRYdO20IP62ykpjQMfERBdGnRTshuJMSUTv8xNVXdvle0wbOhn+f3t3Hh9VeTb+/3NmTTKTfScLZCHsa1hEBRVF1KqoWEEUa0ufqk+LT61Pf7Z1rdbl+fZX9am2ttpd/bZYrQulVmVR2SEgSSAkkEACZN+TSWY/5/vHQDSQQJKZ7Nf79coLmTnnPvfckplr7nPf1wXs6SYBkdHgK23c0uZi10HJtCpGj25nBh5++GFOnjzJwYMHOXr0aMfjHo+H1la5nyZEd7S6Rjw6CInsfr2AejrJkF7RmCmph3ukLSGG2MISLNV12KMj+9TGnDHT0Ss6dp/8gpsmLe3ymKvnj+W9z0r4aFcpC2eeu9BQiJGo22Dgvvvuo7y8nKeffrpTTgG9Xk9GRsaAdE6IYcftIaTJTkWMkWRd9yV3S5qDaHYZmBbdRohxlNcp7qGO5ENVddRNHt+nNqxmC1PiJpBXfZjatnpiLecGbCnxoUxJ95U2rqxrIzHm/FtDhRgJur1NkJyczPz58/nggw8YN24c8+bNQ6fTUVhYiMlkGsg+CjF81DehaFAXG0QY3f+e5Mh2wl6zR0fgNej7XM74jHmnaxXsLc/t9pilF40F4OPdZX5dS4jh4oJrBh5//HFeeeUViouLefDBBzl06BAPPfTQQPRNiGHHXlsLgDem+1mBCpuJyjYz6WF2ooJ6X4lv1NLpaI+LJrihGZ2r78mZ5oyZDsC+ivxuj7l4+hiswUY27j2BxyszN2Lku2AwkJ+fz2OPPcaHH37IrbfeyjPPPENFhWToEqIrSl0TqgLW6Nhuj5FZgb47U8HQUtO3CoYAUSERpEemUlBzhHZX17s4zEY9V8xJoanVyZ5DspBQjHwXDAa8Xi+qqrJp0yYWLVqE3W7HbpdtUEKcw+vF0tBObaSBREPXMwONTj3FzUHEh7hItroGuIPDX5ufRYvOyB4zDa+mcqCqoNtjls733Sr4SG4ViFHggsHATTfdxKWXXkpSUhIzZszglltuYcWKFRdsWFVVHnvsMVasWMHq1aspK+v8C7V582aWL1/OihUreOuttwBwu9388Ic/ZNWqVdx6661s2rSpjy9LiIHnrW9Er0JjbDDmblIQ76+xAgrZcTZJMtQHHZkI/Q4GztwqyOv2mLGJYUwcG8kXRTVUN/StJoIQw8UFkw5985vf5K677kKv9725vfnmm0RFRV2w4Y0bN+JyuVi3bh0HDhzgueee45VXXgF8H/rPPvssb7/9NsHBwdx+++0sXryYzz77jIiICH7+85/T1NTETTfdxJVXXunnSxRiYLTV1RIMqLERXT5v9+g4VB9CqNFDVoTMrvWFyxqCKyTY75mBtMgUooIj+KLyEF7Vi17XdfC29KJxFJY18smeMu68ZpJf1xRiKOs2GHj00Ud56qmnWL16NUoXX2EulI543759LFy4EICZM2dy8ODBjudKSkpITU0lPNw3lZqdnc3evXu55pprWLrUt/dX07SOAESI4UB3OvOgNTquy+fz6kLwaDqy41rQyaxA3yiKr2jRsZMYbW24rX3b9qcoCtljpvFJyVaO1B9jUmzXWxUvnTGG197P55PdJ7h9yQT0+t5nPhRiOOg2GDhzK2Dt2rV9athms2G1Wjv+rtfr8Xg8GAwGbDYboaFf1ni3WCzYbDYsFkvHuffffz/f//73L3idl156iZdffrlPfRQiYFSVsDo79eEGYk3nrhfwqPBFrRWzXmVqtEw5+8MW7wsGrFV1NGb2PQdA9pjpfFKylZzyvG6DgSCzgctnJ/OvHaXsK6xh3pSEPl9PiKGs2zB36tSpAERERNDQ0EB7ezuJiYkdKYkvxGq1dqp2qKoqBoOhy+fa2to6goPKykruuusuli1bxg033HDB66xdu5aioqJOP7LWQAw0e2MDBq9Gc2wI+i5SEB9uDKHdo2dadBsmvaQe9kegFhFOjZ+AWW867xZD8N0qAPhwZ6lf1xNiKOt2ZqC+vp7777+fo0ePMnbsWBRF4fjx48ycOZNf/OIXhIWFnbfh2bNns2XLFq677joOHDhAVlZWx3MZGRmUlZXR1NRESEgIOTk5rFmzhrq6Or71rW/x2GOPsWDBgsC9SiH6WVtdLeGA1sV6Ae106mEdGrNiZTuhv9rio9EUBWtlrV/tmPRGpiVMIqc8l8rWGhJDu769k54UTlZqBPsKq6luaCc+KsSv6woxFHU7M/DUU0+RnZ3N9u3b+fvf/85bb73F9u3bmThxIs8888wFG16yZAkmk4mVK1fy7LPP8uMf/5j169ezbt06jEYjP/rRj1izZg0rV65k+fLlxMfH85vf/IaWlhZ+/etfs3r1alavXo3DITXexdCnP71eIKyL9QLHW8w0OIxMiLQTapIENv5SjUbaYyKx1NSjeLx+tTWnB7sKAK67OA1Ngw93HPfrekIMVd3ODBQVFfHiiy92esxkMvGDH/yAZcuWXbBhnU7Hk08+2emxr9Y0WLx4MYsXL+70/COPPMIjjzzSo44LMVRomkpErZ1mq56w4HNnzHJqfLfA5sRJga9AsSXGYaltwFJbjy2x62/0PTF7zFQUFHLK87h+wlXdHrdwZhK//+AQn+w5waqlEzEZZXGzGFm6nRkwm81dPq4oCro+1BIXYqRqbm7A7NZojbWinLVeoLLNyCmbmbGhDmJDJPVwoNgSfRke/b1VEBEURmbUWArrSrC5uq8eaTLquXp+Ki1tLrblSgZWMfJ0+6ne1XbCnjwnxGjTXlfj+4/Yc8vq7q32zQrMjZdZgUBq7QgGavxuKztpOqqmcqDy0HmPu2bBOBQF/rVdbhWIkafb2wRHjx7tMuGPpmnU1voXjQsxkhhrmwEIi+k8Xd3gMFDcHERCiIsUST0cUG6rBWeoBWtVrW+Fph9fUOaMmc7f8j8gpyKfS8d2v1MqIdrCnEnx7C2opvhkE5kpXSeXEmI46jYY+OijjwayH0IMS27NS1SNHZtFT5AltNNzOadTD8+Nb5XUw/3AlhBL9NFSzM2tOCPOv7vpfFLCxxAbEsWBykN4VC+GbrIRAnztkjT2FlSzYftx/mvlrD5fU4ihptvbBElJSef9EUJAfXMtwS4NW6y10+OtLh0FDSFEmt1khMuOmP4QqHUDvmyE02l32ymsPXreY2dlxZEYY+HzL07R0iazPWLkkJWAQvjBVetLfKOLi+70+P5aK6qmMCfOJqmH+8mZXQSBWTcwDYB9FQfPe5xOp3DdxeNweVQ27jnh93WFGCokGBDCD0E1vvUCobHxHY85PAr5dRYsBi+ToiT1cH9pj47AazQQWuX/GqbJseMxG8zsrzx/NkKAK+emYjLo+HDncVRVskmKkUGCASH6qF11EVvrxGY1oA/5MitdXp0Fl6pjdpwNg/yG9R+dDlt8DMENzegdTr+aMuqNTI+fSGVrDZWt559pCA0xcdnsZKrq29lf5P+shBBDgbxVCdFH9c01mN0a7XFfLl7zqL5bBGa9yvSY7veti8DoWDdQ5V+dAoDZib56LPsvUKsA4LpL0gDYINsMxQghwYAQfeSurQfAGPvleoGCBl9BoukxbZilIFG/syUEbt3ArDGng4HK868bAMhMjmDC2Ej2FVZTUSf1JsTwJ8GAEH2goWGp8SUSCon1fSB5NV+SIb2iMVsKEg2ItoQYNPDlG/BTVHAEaREpFNQexe6+8A6QGxemo2mwfusxv68txGDrNs+AEKJ7RfUNTKx10WI10mDzgM3GsbYwml0GxlsaaW5qoXmwOzmExR48ErC23JZgrJW1xOYVQlep0pf0vK3ZY6ZxvOkk+dWFzEueed5jL54+hpjwQ2zcc4I7rpmENdjYy54LMXTIzIAQfeB0NGLyaNiifPkFVA0OtUSjQ2VSWMMg9250cYSHolNVTDb/d27MHtPzdQMGvY7rL03H4fLy8a5Sv68txGCSYECIPjA1+b73q2G+egSl7WG0eU1kWJsJ0UtBooHkDPNlfgxq9r/+Q0bUWMLMVr6oPISqXbjc9NKLxhJk0rN+6zE8XilPLYYvCQaE6CW3101Uve9bqBoW0XlWIFRmBQaaI+J0MNDU4ndbOkXHzMQpNDqaKW08ecHjrSEmrpqbSl2zgx15Us1QDF8SDAjRS0VVR0msc9ESakI1GTneHi6zAoPIG2TGHWTG3NzqK1rkp9mJvmyEPdlVAHDDonQUBd79tBgtANcXYjBIMCBELxXt24rBC22RoagaFLREoUNlcmj9YHdt1HJEhKH3eAOybmBGwiR0io79F0hNfMaYGCsLpiVSfKqZvKP+5zsQYjBIMCBELzXnnf6QCIvsmBXItDYTrPcObsdGMcfpqoWBuFVgMYUwKTaTkoYymh09a2/5FeMBeHvL+QsdCTFUSTAgRC+0Om1YymrRAEdYeMeswCSZFRhUgVw3ADArcSoaGl9UHurR8VmpkUzPjOHAkVqKTzUFpA9CDCQJBoTohfxTB0moc9MaZuaYO1pmBYaIgK8b6EU2wjOWL/bNDryzWWYHxPAjwYAQvXBs304MKrRHhFLQEoVeZgWGjECuG0gKTSDOEk1uVQEetWeB3qysWNKTwtmRVyEpisWwI8GAED2kaRpthw4DUGFKPr2DoElmBYaIQK4bUBSF2YnTsLsdFNYW9/icW68Yj6rB25tkdkAMLxIMCNFDlbYaospb0BTYqWViULxMkbwCQ0ag1w305VbBxTPGkBxnZXPOSaob/J+hEGKgSDAgRA/lluUSX+emNTqKViWESaENmGVWYMgI9LqByXFZmPUmvujhFkMAvU5hxVVZeFWNt2XtgBhGJBgQoofK9u9Cr8ERJY4gnYcJ1sbB7pI4SyDXDZj0RqbFT6S8tYpqW8+rIi6cmcSYGAsb95RR22j3ux9CDAQJBoToAbfXjfdwCQAlhiSmhNVj0Em2uaEmkOsG4KuFi3oxO6DXcdtVWXi8Gu9I3gExTEgwIEQPFNWVkFThwKNT8KakkWGRveRDUX/kGwDYX3nhKoZfdfnsZBKiQ/hoVxk1jbJ2QAx9EgwI0QN5xfuJbfJwyhLFN26ahU4Z7B6JrgR63UB0SCRjI5I5VHMUh9vR4/P0eh23Xz0Bj1dl3SdH/O6HEP1NggEheqB89x4A7IkTmDs5fpB7I84nkOsGAGYnTsWjesivKerVeZfNTiEl3srGvScor5W8A2Jok2BAiAuob2sk4qRvAdnC5UtQFJkWGMockafXDTQ2B6S97DGnqxj2Yt0A+HYW3HnNJFRV481/FwakL0L0FwkGhLiA17d+TmqlC6fZRNa8qYPdHXEB9shwAIIbAhMMZEaNI9Rs5YvKg70uUbxgWiKZyeFsPVDOsfLA9EeI/iDBgBDn0dTqpDB/D1a7inXqZBSd/MoMdarJiNMaQlBzK4rX/zwQOp2OmQmTabA3UdZ0qlfnKorC6usmA/DnDQV+90WI/mIY7A4IMZT9+V+HSLFVApC6YMEg90b0lCMyHLOtnaCmVl785D2/26tz+hYPvrT1LVKCM7o9bqJ15jmPaZpGcpyV/UU1vPpuHqkJYX7355oF4/xuQ4ivkq85QnSjqKyBTYcOMrbGlzgmcua5b/RiaLJH+W4VBGrdQIQxBlBodPU8+dAZiqJwyfQxAGzPq0QNwC4HIQJNggEhuuBVNX7zjzyMoTUk1bghLgpzbMxgd0v0kCM8FFWnC9i6AYPOSJghApu3Gbfq6vX5MRHBTBoXRUOLg8JSqWchhh4JBoTowobtxyg+1cwEcyUmj0b07NmD3SXRGzodjohQTO12TK1tAWky0hgLQKO7rk/nz5+SgEGvY/ehKlxuqWkhhhYJBoQ4S12TnTc+PIw1VCW2yjctHDt7ziD3SvSWPSoCgPCy8oC092Uw0PtbBQCWYCOzJsTS7vCQc7g6IH0SIlAkGBDiLK++l4/d6eWSSw2kVDnRdArh06YMdrdEL9mjTwcDpYEJBoL1Fsy6IJrcdaia2qc2Zk+IIzTERO7RWhpbep7RUIj+JsGAEF+xM7+SnfmVTEmPxkUp8Q0egjLTMYSEDHbXRC95goNwhQQRdqoSxeP/tLyiKEQaY/FqHlo9fatNYdDruHTmGFQNPj9Q3uu8BUL0FwkGhDitpc3Fr9/JxaDX8Z2bJ9OSm49Og/g5cwe7a6KP7FER6D1eQsurAtKev7cKANISw0iND+VUjY0SSUQkhggJBoQ47bfv5tHU6uTOaybSolSRfMK38CxybvYg90z01ZlbBREBulUQZoxCh65PWwzPUBSFhTOT0OkUtuVWyGJCMST0WzCgqiqPPfYYK1asYPXq1ZSVlXV6fvPmzSxfvpwVK1bw1ltvdXouNzeX1atX91fXhDjH9rwKPv+inAljI7np8kz2ncplXKUTJSIMS1raYHdP9JEjPBSPyehbRBiAKXm9oifCGINdbaPd2/fiQxGhZrInxNFmd7PrYKXf/RLCX/0WDGzcuBGXy8W6det48MEHee655zqec7vdPPvss/zhD3/g9ddfZ926ddTV+bbrvPbaazzyyCM4nc7+6poQnTTbnLzyTi4mg47vr/SVJy7L3UOwUyNm7lwpTDSc6XS0pI4hqMVGcH3f7vOfLdrkq1rZ4PJvR0D2xDgiQ83kl9RTWReY7Y9C9FW/BQP79u1j4cKFAMycOZODB7+s+FVSUkJqairh4eGYTCays7PZu3cvAKmpqbz00kv91S0hOtE0jV+/k0uzzcXq6yaTHBfKyeYKIo/5gtPoubJeYLhrTE8BIPLYyYC0F2mMRUGh3lXjVzt6vY4rsn1927LvJF5v33YoCBEI/RYM2Gw2rFZrx9/1ej0ej6fjudDQ0I7nLBYLNptvym3p0qUYDD0vmfDSSy8xYcKETj9XXnllgF6FGOm2HahgR14lk9OiuGFhOgA5FXmklbvQDHoiZkwb5B4KfzWPTULV6Yg4diIg7Rl0RsKN0bR5W3B42/1qKzHGwrSMaBpbneyV3ANiEPVbMGC1Wmlr+3LqS1XVjg/5s59ra2vrFBz0xtq1aykqKur0s2nTJv86L0aFxhYHr/wjF7NJz3+tnIVe57sdcKgoh9gmD2FTp6APChrkXgp/ec0mWlISsNQ1YmppDUib0cY4ABr8nB0AuGhqIqEhJvYX1lBVL7cLxODot2Bg9uzZfP755wAcOHCArKysjucyMjIoKyujqakJl8tFTk4Os2bN6q+uCHEOTdN4+e+5tLa7uftrkxkT45vFamhvgkPFAMTOmzeYXRQB1JieCkBkSYBuFZh8wUC92/9v8yajnivnpqABG/eewO2R2wVi4PVbCeMlS5awfft2Vq5ciaZpPPPMM6xfv5729nZWrFjBj370I9asWYOmaSxfvpz4+Pj+6ooQ51i/7Rh7CqqYMT6G6y7+crfArlP7SSv3FaKJnCNbCkeKprRktC2+dQPVsyb73Z5JZybMEEmLpxGX6sSkM/vVXlKslRnjY8k9WsvO/AoWzUr2u49C9Ea/BQM6nY4nn3yy02MZGV/WAV+8eDGLFy/u8tzk5ORzthsKESjFp5r44/oCwq0mfrAqG53uy90Cu4/tZXG1C3NKEkHxcYPYSxFInpBgbGPisFbUYLS147b6n1Ey2hRPi6eRelc1iUGpfrd30dQETlS3kF9ST2pCGOMSw/xuU4iekqRDYlRpd7j5+es5eLwqD9w+m6iwL9cE1Lc34igowuCFGNlFMOLUjx+HAkQdLQ1Ie1Ed6wYCs/DPoNdx9byx6HUKm/aewNbe+1LJQvSVBANiVPnNP/KoqGvjlsszyZ7Y+dbU7lNfkFbuy28RNVeqFI40jZlj0RSF6AAFA2Z9MFZ9OM2eRtxqYD64YyKCuWTGGBwuL5/sOYGqSu0CMTAkGBCjxuacE2zZd4qs1AjuvHbSOc/vLMthXIULvdVC6ISsLloQw5knOIjmlEQsNfWYm1oC0qYvAZFGfYBmBwCmpkeTkRRORV0bewsCU1NBiAvptzUDQgy0f+8s7fa5xlYHb208ismgY97kBDbt7bznvM3TQuPRI4S2qzimTOCjPYFZdS6GloasNCJOVBB9pJSKedP9bi/GlECZ/Qh1rkoSglIwHtgZgF7CEhPUGoLIKawhsfkkGZaz6hcsGBeQ6whxhswMiBHP7VH5aFcZHq/K5dnJhFvPXfld2l7E+JO+WwSeCf5/SIihqTE9BVWvJ+rI8YDUKjDrgwkzRNHiacThtQeghz5Bergh0YlB0fioykSDS1Jii/4lwYAY0TRN49P9J6lvdjA5LYrxKZFdHne87TCZJ5yoRhOetIkD3EsxUFSTkcb0FIKbWrBW9b3y4FfFmhMBqHNVBKS9L9vVuCrOhUtTWF9pxinpB0Q/kmBAjGh5xXUcOdFEfFQIi2YmdXmMzdOCWnOCCJsXb8YkMBoHuJdiINVNzgQgpqA4IO1FG+PRoaPWWYlGYBf8TQz1MivcTaNbx7+qzMh6QtFfJBgQI9apGhvb8yoICTJwzYJx6PVd/3M/3n6YzDO3CCbOGMguikHQkpyAM9RCVHEZOpfb7/YMOiORpjjsaht1+r6XNe7Owhg340K8lLXr+bTOGIi7G0KcQ4IBMSK1trv4aFcpCnDNReOwBnf9bV/TNI7a8hl/woVmMOBJl1sEI56iUDs5E73bE7CcA7Em362CI2b/axWcTafAtQlOYkwqec1GvmiWdd8i8CQYECOOx6vy4c5SHC4vC2cmkRhj6fbYenc1htoqolo8eNIngcm/tLJieKifmIEGxBYcDUh7EcYYDIqRElMNaoBvFQCYdXBjopMQvcbndUY+238q4NcQo5sEA2JE8S0YPEVto51J46KYkh593uOP2vKZWHr6FsGU2QPRRTEEuEItNI1LxlpdjyUACwl1io4YUwJ2nZtTxsYA9PBcYUaNm8Y4MOnghb/uZ1+hlDwWgSPBgBhRvjhSS1FZI3GRwSyalYSidL8ly6t5OG4rYEKZE80chCfj3EREYuSqnuG7JRSfWxiQ9mLNYwA4Yuq/D+k4s8aNiU70OoVn/7yXw8cb+u1aYnSRYECMGEdPNrEzvxJLsJFrL07D0M2CwTPK2o8QU92Ctd2Le8J0MMgugtGkNTmB9qgIIkvKMNra/W7Pqg8nwhvMcVMddsX/hYndSQ5Weeiuubg9Ko+/tpPCMgkIhP8kGBAjQsHxejbtPYHRoOP6S9K6XTD4VYW2XCaUOgDwTJFyxaOOolA9YyI6VSMuvygAzSlMciSiKhpHzP2bRnjelAT++45snG4vj7+6kyMn+ufWhBg9JBgQw15FrY2f/WEPqqZxzUVjiYkIvuA5Te566trKyDrhRg0Nx5uaPgA9FUNN/YQ03EFm4vKL0Dv9LzaU5YpHr+koMPc954DS0oShMBfTzs2YN76H+ZP3MG9ej3H/dvRlxeDxzTosnJnEg6tm43B6eOy3O+SWgfCLBANiWGu2OXnid7tobXdx2axkUhN6VgO+yHaA8SccmNxe3NPmgiK/CqORZjBQNXsKBpeb+NzDfrcXpBnJcMXQondQbmjq1blKcwPGPZ9h3rERQ+lRlJZGtKAQNIsVdDr0NZUYDx/A/OkGyt78K16nk0WzkvnBqmzsLi+PvrqDA0cCv7VRjA6yYVUMWy63l6f/uIfKuja+fuV44iJDenSeW3Vy1JbPzSUuNBTc0+f3c0/FUFYzNYuE/YeIzy2kesYkvGaTX+1NdozhiLmGg0HlJNu6Tn/dideD4chB9GXFKIA3Og7v2EzU6DjQf+Ut2t6OvqIMQ1kxp956m7ptO8j6/loum51FkEnP/7yew09/t5sf3pnNxdPH+PUaxOgjX4fEsOTxqvzPX3I4XNrAollJ3HlNz3cCHGnLx9rURkKtE++48WgRUf3YUzHUqSajb3bA6SL+gP+zA/HeMOI8oZQZG2jWnb94kdLWimnXFgxlxWghVlxzF+Geuwg1bkznQAAgOARvxiSci64l8Yav4aisJP/Hj1D9yUbmT03k8W9fhEGv8Nxf9rJ+6zG/X4cYXSQYEMOOV9X4xZv72FNQxczxsfzXilnodD2r6qZqKodacpha4sst4J55UX92VQwTNVOzcAcHkfDFIYy2Nr/bm+ZIAgXyg8q7PUZXX4Np5yZ0rc14UtJxXbLENxtwIQYD6d/+FlOefBx9SDDFL79C6V/eYHpmDM/+56VEWM28+l4+v3v/IF4pZiB6SIIBMayoqsYv133BttwKJqdF8fA352Ey6nt8/gn7URzOJqYcd6KGWPFkTunH3orhQjUZObVgFnqPl5Qd+/1uL90Vi9VrpshchaOLbYa6ypMYc7aCV8U1ba4v4ZW+5/+OASKmT2P6z58jaMwYyt95l+O//yMZyeH8//cvIiXeyvufl/CzP+zG1u7/wkgx8kkwIIYNTdP4zbt5bM45yfiUSnhu3wAAGdlJREFUCB7/9kUEmXu+7EXTNHKbdzL5mB2T04N71gIwyLIZ4VM3KYO2uGiij5RirfBvIZ4OhWnOJDyKysGzZgf0pUcx5e4GvR73nEtRk8b2+TrBiYlMe/YpQlJTqFy/geO//yOxkcH8n+8tZPaEOHIOV/PAi59xvKLZr9cjRj4JBsSwoGkaf/xnAR/uKGVcYhg//c4CQoJ6lyTolOMY9a4q5h7xoun1uGdf0k+9FcOSolC2cC4A4zbvROf2+NXcJEciQaqRfHM5TsUDmoahKA9jYS6aOQjXvMt7dlvgAkwREUz92U8JGZtK5foNlL/zLtYQE499+yJWXJVFVX07//3LrXwq9QzEeUgwIIY8TdN449+FvPtpMclxVp6652JCQ3q34lvTNL5o3k5ahQtrix3P5NloltB+6rEYrtoSY6maMZHgphaSd+zzqy0jeqY7knDpvBQYT2HM24vh+BFUSyiu+VeghUUEqNdgDA9n8uOPYI6Noez1N6netBm9TuHOayfx8DfnYdAr/OLNfbz2Xj4erxqw64qRQ+ZIxZCmaRp/WH+I9z4rITHaws/uvZiI0N5XFjzpKKHOWcHSIt+CKtfcRYHuqhiiYg8e6dXxzrBQXCHBxOcfAQ3sMT3YHnhGUEqnv05xjKFAf4KU3Xnoq12o4VG4si/xuzrmv3eWdvm4btm3CHnjZY6+9Ap5FQ68mZMBuOmyDD7cWcoHW4+x+1AVV81N7fHv0TULxvnVVzE8yMyAGLJUVePX7+Tx3mclpMRbefa7lxAdfuHsgue0o6nkNH5KSrWHmKoWPBmTfFu3hOiCptdROzkDVacQW1CMqbXvuwtMDjcrN7aQXO2iLiEc17xF/VomW42Jp/3ra0CvJ/i9v6ArLwMgMjSIWxePJys1guqGdtZtPEJ+SR2aJrsNhI8EA2JIcrm9/J83cvj3zlLSx4Tz7H9e2qdAAOCILY8mdx2XH/K98TkvXRrAnoqRyG21UDcpE0VVic8rwtDu6HUbSlM95l2bsTS3U5huYd1lZlqN/q1D6Ak1aRz2m+4Cr5fgt3+P0lgHgMmgZ8m8sSy9aCx6vcLnX5Tzz23Hsdn7r6iSGD4kGBBDjs3u5vHXdrI9t4Ip6dE8/Z+XEG7t27cpp+pgf/NWxlWrRFU34cmcjJqYcuETxajXHhtFQ+ZY9G43ifsPYW5q7dmJqor++BFMuz8Fhx33+Cm4pszCo4c9wcf7tc9neDMn41x6Czp7GyFvvQb2L2c3MpMjuH3JBFITQjlR3crfPi6i6ESjzBKMchIMiCGlqr6N/++lrRwsqeeSGWN48jsLelSBsDv7mj7D6WnjyjwvILMCondakxOoy0pD5/WSkHuYiOOnULzerg/WtNOJhDZjLMoDown33IV4MyaR6Y4n1hNKsbmWU4aBqTDonrkA50WL0TXWEfzOHzoKHAFYgo1cf0kal81OxqtqbNxzgvc/P0ZjS+9nQMTIIAsIxZCRV1zLc3/eS2u7mxsXpbPmhqk9zizYlRpnOYW2A8wtM2Ctq8U9JRs1ITmAPRajgW1MHJ5gMzGHS4goK8daWUN7bDT2yDBUgx6dV8XUasNS24CpzZd+2JM0Ds+EaR3rA3QoLGzL5N2wL9hmKebW5mwMA/BdzHXZteiaGzAePkDQP/+KY9mdHUW5FEVhano0KXFWth4op6yqlb99coSZWbHMmRSP0SDfFUcTCQbEoNM0jfVbj/GH9YdQFFh720yunt/3RCwAbtXF1vp/YXapzP+iBc1kxnnF9QHqsRhtHJHhlM+bQfiJCsLKqzp+vkpTFLwJKXjGZaJFRJ/TRqw3lCnOJA4GlZMTXMpF9gEom63ocHztdpTWFl9+A0sYzquWgfJlkB1uNfO1S9I4XtnC1gPl7C+q4cjJRhbOSCJtTM+qgIrhT4IBMahsdje/XPcFO/MribCa+dE35jIl/dw30t7a07SZZk8Dyw8Fo3fU4bjiejSrvLGJvtMMeprSU2gal0RQcyvmFhuKqqIpCi6rBWeYlbiw83/Az2sfxwljPblBp0hxR5HkCVyugW4ZDNiXf5OQ//srTPu2opnNuBZd2+kQRVFIHxNOSpyVnMM1HDhSy4c7S0mMtjA2MYzJaf7/ToqhTeaBxKApOF7PAy98ys78SqZmRPPiDy4LSCBwvL2QIlsuU6uDSD58Em9sIu45CwPQYyEAnQ5HZDjNY5NoSkuheVwy9phIVNOF17YY0bPYNhEF2GItxK4MUN2A4BDsK+5BjYzBvGMjpp2buu6fQc+CaYmsXJLFuMQwKuvbeOjlbTz1+92UVbYMTF/FoJCZATHgXG4v//ejQv7xaTEAt12VxaqrJ6DX+x+b1ruq+bx+A1annst31qPp9ThuuOPccrBCDJJ4bxhz7OPYG1LKx9YCrm+djn4Avpdp1jDaV95LyJu/wvzZv9CMpm6D5MiwIL52SRqVdW0UnWhkT0EVew9XcUV2CrddlUVSrLXf+ysGlrxDigGVV1zLK+/kcarGRmK0he/fPitgU5BtnhY21r6D5nWzYrcJfXsbjituQI1LDEj7QgTKLEcK9QYbx0x1bAspZlH7eBT6vli2p7TwSNpX3kPIm78iaON7oGm4z5ONMzHGwt3XTybncDV/+ddhNuecZMu+kyyYlsjyK8aTldqL7IxiSJNgQAyI+mY7f/pnAZ/uP4WiwPWXpHHX1yYT3Iuqg+fT5Gjhw5p1tHlauC0vFOupY7jHT8E9T9IOi6FHQeFy2wRawhwUBlVh0gxcZE8bmIAgKhb7ynsJXvdbgja9j2Jvw7Xwmk6LCjv1VVGYOzmB2RPj2ZlfwTubj7Ijr5IdeZVMy4jhlisyyZ4Yh9LN+WJ4kGBA9Ks2u5t3thzl/c+P4XJ7yUyJ4D+XT2d8SuC+UdS21fPM5y/T4mnghuIoEguK8MaN8d0eUGRZjBiajOi5rnUqH4Tlkhfsqyg4UAGBGptA+53fI2Tdq5h3bETXWIfjuhVg7L4AmF6ncOmMJC6ZPoa84jr+saWY/UU15JfUkRRr4er5Y1k8p+c1D8TQIsGA6BfNNifrtx1jw+l0p9HhQaxaOo0r56ai9yN3wNmONZTxP1tfodHexI1F4aTtL0INi8S+/Fv9mgNeiEAI1kxc3zKd9WF55AWfol3n4rK2rAHJQaBFRNN+5/cIevdPGA8fQFdfg33ZarQLlFVWFIUZ42OZMT6WY+XNvPtZMdtzK/jjPwt4/cPDzJ+SyNXzxzIzK9avPCFiYEkwIALqeEUzH+4oZVPOSVxuL2EWE9/42mRuWJiO2agP2HU0TePDo1t4I/ddcHm491g05v2FqBHRtN9+H1q43MsUw4NFM7OsZSb/Dj1IsbmGRn0bV9kmEaGG9Pu1NUso9tvvw/zJu5gO7MLyx+dxLr4B96wFPZpVS08K58FV2dxz0zS27DvFx7vL2J5Xwfa8CqLCglgwLZGLpycyJS06IAuERf+RYED4rdnmZFtuBVv2naSozJdqNSYimJsvz+Dq+WMJMgX2n9mp5kpe2/dXDtceZWybiVv261FPFuJNSMa+/FtooeEBvZ4Q/S1YM3JDy3R2hJRwOKiKd8L3M9ueynRHcv/vNNAbcF7zdbxjxxP00dsEffwPjAdzcCy5pcd1PKwhJm5YmM71l6Zx9GQTH+8uY0deBRu2H2fD9uOEWUzMn5LARdMSmZoeTUhQ31OMi/4hwYDok4o6GzkF1ew9XE1ecR2qqqEoMHtiHNcuGMfcSfEB/yZQ0VrNuwX/5vOy3ZgdXm49YSFp/ylUr5e4q66kZOYSMMibjBieDOhZ1J5FkieS7SHF7AkppdBcxSx7KuNdcf0eFHgmzaQtOQ3z5g8wHj6A5c8v4h4/hdbo1YRmje9RG4qikJUaSVZqJPfeMp2DJXXsyK9kV34ln+w5wSd7TqDTKYxPjmBaZgzTMmOYPC6KoAAtJBZ9p2gjsFTVqVOnuPLKK9m0aRPJyZKL3l92p4eyqhaOlzdTUNrAoWP11DbaO57PTIngsllJXDojiZiIvpUZ7k6bq52c8jw+L9tNftVhYpo8zD+lJ/NIMzjdmOPiSL/n20TNyebfO0sDeu3zqd/11oBdSwwfCUGBqYjpVNzkBJdRYK5EVTSCVSMTnAlkuGKJ9lrwzLw4INfpjr6sGPPnH6IvLwXAkpFO/JVXEHXRfMzRvd8K7FU1Cksb2FdYTX5xHUdPNuFVfR89ep1CSnwo6UnhX/6MCcfiR4Gy0cjfz71+CwZUVeWJJ56gqKgIk8nEz372M8aO/TLf/ObNm/nVr36FwWBg+fLl3HbbbRc8p6ckGOg9p9tLXZOdukY7tU12qhvaKatqobSihcr6tk7HhoaYmJoRTfbEOOZMiic6PDABgMfroaK1mrKmcsqayymsOUrNiRJiGlwk1bjJrAXr6SDEFB1F0s3LiL96CXqzb6GgBANisAUqGDjDpnOQF1ROkakal84DgEU1EWtNI96UTKw5kTBDFGZdUOC39mka+rKjJB//goa9OaCqvuunpxE2ZQrWzAxCkpMIThqDPrh37wF2p4eC4/XkF9dx6Fg9xytbcLo6V4OMDDWTEG0hMcbi+zM6hJiIYCLDgoiwmgkJMsh2xq/w93Ov3+ZmNm7ciMvlYt26dRw4cIDnnnuOV155BQC3282zzz7L22+/TXBwMLfffjuLFy9m//793Z4jLqy13UV+cR0utxe3R8XtVX1/dvx4sTs92OxubO1u2uxubHYXzTYXLW1dp0UNDTEyPTOGcYlhjE0MY8LYSFLiQvu0Stijellf+An19kbcXg9u1UO7q51GRzNNjhZCS+uYeqQdi91LtF3laruKXv3yfJ3ZTMRF84m7/DIi58xGZ5RvDmJks6pBXNyewdz2cZw0NnDMVEeFsYnS9iJK24s6jjMqJqyGcIL1IRgVMyadiQhjDFND5/X9A1NR8I7LYtLtV+Osb6Bh127qd+2m5XAhbceOdzrUFB2NKSoSg8WCwWpFb7Gg6HUYrFaSl9+MPiio0/HBZgPZE+PJnhgP+GYOKmptHCtv9v1UNPuyH5Y1cLi0ocvumQw6IkLNhFnNWIIMhAQZCTnzp9mA0aDDoNf5/jToMOp9fxr0Ogx6BVDQKTBxXBThVtl51G/BwL59+1i40JfqcubMmRw8eLDjuZKSElJTUwkP9y30ys7OZu/evRw4cKDbc3rDe7reeFVV1QWOHFn+9M9DbMut6PHxik7BEmTAEmRifFwQkWFmosKCiAwLIjo8iKRYC+FWc+c3E08LFRV9y1Fe197Anz//2zmPG/RGws1W0o95iCxpRdPpcIdZsCdGERwXS0hKCpa0NKzpaSgGA3bAXl19Tjv1tQP3/7uppXXAriWGD6Ozsd/aDsfALBKYSQJNEzKpd1XT4K6h3dOKzdtCnacSr+buOF5BT2xiIkal+9wBPXHq1OmPiRnTiJgxjXC3G1tpKfZT5Tgqq3BUVdFcVY3nSC2a23PO+c5xYwkZm3rB6yhARhxkxIXDLN9ng1fVqG+2U9tgp6apncZWJy02J81tvi8wLS2t1FS5cXvV8zd+HnMnJ3DPzdP6fP5Qcebz7sznX2/1WzBgs9mwWr/MX63X6/F4PBgMBmw2G6GhoR3PWSwWbDbbec/pzksvvcTLL7/c5XN33HFHAF6JGCjbBrsDQowwBUPht+pb3xzsHpzX8c3wVtcfIcNSbW1tn26v91swYLVaaWv78l6zqqodH+pnP9fW1kZoaOh5z+nO2rVrWbt2bafHHA4HM2bM4OOPP0avD9zedtHZmftTov/IGA8MGef+J2Pcv7xeL1dffTVTp07t0/n9FgzMnj2bLVu2cN1113HgwAGysrI6nsvIyKCsrIympiZCQkLIyclhzZo1KIrS7Tm9EXT6/lRfoiPRO7JAs//JGA8MGef+J2Pc/4LOWp/RU/0WDCxZsoTt27ezcuVKNE3jmWeeYf369bS3t7NixQp+9KMfsWbNGjRNY/ny5cTHx3d5jhBCCCH6V78FAzqdjieffLLTYxkZGR3/vXjxYhYvXnzBc4QQQgjRvyRZtBBCCDHK6Z944oknBrsT/WX+/PmD3YURT8a4/8kYDwwZ5/4nY9z/+jrGIzIdsRBCCCF6Tm4TCCGEEKOcBANCCCHEKCfBgBBCCDHKSTAghBBCjHISDAghhBCj3IgJBn7729+yYsUKbrnlFv7+979TVlbG7bffzqpVq3j88cdR1b5XtRK+stMPPvggK1euZNWqVZSUlMgYB1Bubi6rV68G6HZc33rrLW655RZuu+02tmzZMpjdHZa+OsaHDx9m1apVrF69mjVr1lBXVwfIGAfCV8f5jPXr17NixYqOv8s4++erY1xfX899993HHXfcwcqVKzlx4gTQhzHWRoBdu3Zp99xzj+b1ejWbzab98pe/1O655x5t165dmqZp2qOPPqp9/PHHg9zL4e2TTz7R7r//fk3TNG3btm3a9773PRnjAHn11Ve166+/Xvv617+uaZrW5bjW1NRo119/veZ0OrWWlpaO/xY9c/YY33HHHVpBQYGmaZr217/+VXvmmWdkjAPg7HHWNE07dOiQdtddd3U8JuPsn7PH+KGHHtI2bNigaZqm7dy5U9uyZUufxnhEzAxs27aNrKwsvvvd73Lvvfdy+eWXc+jQIebNmwfAokWL2LFjxyD3cnhLS0vD6/Wiqio2mw2DwSBjHCCpqam89NJLHX/valzz8vKYNWsWJpOJ0NBQUlNTKSwsHKwuDztnj/Hzzz/PpEmTAF+1N7PZLGMcAGePc2NjI88//zw/+clPOh6TcfbP2WO8f/9+qqurufvuu1m/fj3z5s3r0xiPiGCgsbGRgwcP8r//+7/89Kc/5b//+7/RNA1FUQCwWCy0trYOci+Ht5CQEMrLy7n22mt59NFHWb16tYxxgCxdurRTqe6uxtVmsxEaGtpxjMViwWazDXhfh6uzxzguLg7wvZG+8cYb3H333TLGAfDVcfZ6vTz88MP8+Mc/xmKxdBwj4+yfs/8tl5eXExYWxp/+9CcSExN57bXX+jTG/VaoaCBFRESQnp6OyWQiPT0ds9lMVVVVx/NtbW2EhYUNYg+Hvz/96U9ceumlPPjgg1RWVvKNb3wDt9vd8byMceDodF/G6GfG1Wq10tbW1unxr/6yi97717/+xSuvvMKrr75KVFSUjHGAHTp0iLKyMp544gmcTifFxcU8/fTTXHTRRTLOARQREdFR9G/x4sW88MILTJ06tddjPCJmBrKzs9m6dSuaplFdXY3dbmfBggXs3r0bgM8//5w5c+YMci+Ht7CwsI5/TOHh4Xg8HiZPnixj3A+6Gtfp06ezb98+nE4nra2tlJSUkJWVNcg9Hb7ef/993njjDV5//XVSUlIAZIwDbPr06WzYsIHXX3+d559/nszMTB5++GEZ5wDLzs7ms88+A2Dv3r1kZmb2aYxHxMzAFVdcwd69e7n11lvRNI3HHnuM5ORkHn30UZ5//nnS09NZunTpYHdzWLv77rv5yU9+wqpVq3C73TzwwANMnTpVxrgfPPTQQ+eMq16vZ/Xq1axatQpN03jggQcwm82D3dVhyev18vTTT5OYmMjatWsBmDt3Lvfff7+M8QCIjY2VcQ6ghx56iEceeYS//e1vWK1WfvGLXxAeHt7rMZZCRUIIIcQoNyJuEwghhBCi7yQYEEIIIUY5CQaEEEKIUU6CASGEEGKUk2BACCGEGOUkGBBiFNu9ezezZs1i2bJl3HjjjVx77bW88sornY7ZuXMnd955J0uXLmXJkiXcf//9nZJ6fbWtswvU+GvChAkBbU8I0bURkWdACNF3U6dO5fXXXwd8mcquu+46lixZQmZmJjk5Ofzwhz/k5ZdfZubMmQC8+eabfPe73+Wdd94ZzG4LIQJIZgaEEB0cDgd6vb4j2+Svf/1r7rvvvo5AAOCOO+7guuuuw+VyddtOWVkZ3/zmN7n55pu5/fbbKSgooLGxkUsuuaQjjfWRI0e44YYbAHjvvfe4+eabWbZsGT/5yU9wOp39+CqFEGeTYECIUe7gwYMsW7aMG264gcWLFzNv3ryOQj4HDhxg7ty555yzZs0aTCZTt20+9NBD/PCHP+Tdd9/lqaee4oEHHiAyMpLp06ezbds2ADZs2MCNN97I0aNHeeutt/jb3/7G+++/T3R0NL///e/758UKIboktwmEGOXOvk1w77338uqrr3LPPfcAdFRQdLlcfP3rXwegubmZ559/ntmzZ5/TXltbGwcPHuTHP/5xx2Pt7e00NjaybNkyNmzYwBVXXMGHH37IX/7yFzZu3EhZWRm33XYbAG63m8mTJ/fraxZCdCbBgBCig8Vi4aqrrmLHjh0ATJs2jf379zN+/HhMJhPvv/8+AKtXr+5UtfKrVFXtdCxAVVVVR3W1Z599lr1795KQkEBCQgJer5drr72WRx55BPAFE16vt59fqRDiq+Q2gRCig9frZc+ePR3fzNeuXcuvfvUrcnNzO44pLCzk5MmT6PX6LtsIDQ1l3LhxHcHA9u3bueOOOwAwmUwsXLiQZ555hhtvvBGA+fPn88knn1BfX4+maTzxxBP8+c9/7s+XKYQ4i8wMCDHKnVkzAGC325k2bRr/8R//AcCcOXN44YUXePHFF6mrq0PTNCIiInjooYfOW7L65z//OU888QS/+93vMBqNvPDCCx23G5YtW8YHH3zANddcA8DEiRP53ve+xze+8Q1UVWXSpEl85zvf6edXLYT4KqlaKIQQQoxycptACCGEGOUkGBBCCCFGOQkGhBBCiFFOggEhhBBilJNgQAghhBjlJBgQQgghRjkJBoQQQohRToIBIYQQYpT7f8HWe4P9frZzAAAAAElFTkSuQmCC",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112953400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's visualize the data\n",
    "import matplotlib as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set(style=\"ticks\")\n",
    "\n",
    "\n",
    "bins = np.arange(70,150,6)\n",
    "\n",
    "A = df['before'].values\n",
    "B = df['during'].values\n",
    "C = df['after'].values\n",
    "\n",
    "\n",
    "# Show the results of a linear regression within each dataset\n",
    "ax1 = sns.distplot(A,bins=bins,label='Before treatment')\n",
    "ax2 = sns.distplot(B,bins=bins,label='During treatment')\n",
    "ax3 = sns.distplot(C,bins=bins,label='After treatment')\n",
    "\n",
    "\n",
    "plt.pyplot.xlabel('BG level')\n",
    "plt.pyplot.ylabel('Distribution of BG levels')\n",
    "plt.pyplot.legend(bbox_to_anchor=(0.45, 0.95), loc=2, borderaxespad=0.)\n",
    "\n",
    "plt.pyplot.xlim((60,160))\n",
    "plt.pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sum of squares for groups $SS_\\textrm{groups}$ \n",
    "In this example, a group refers to a time (before, during or after treatment).      \n",
    "\n",
    "$SS_\\textrm{groups} = \\sum_{i=1}^k n_i(\\overline{x}_i - \\overline{x})^2$\n",
    "\n",
    "$n_i = $ number of subjects under each $i^\\textrm{th}$ group  \n",
    "$\\overline{x}_i = $ mean score for each ($i^\\textrm{th}$ group  \n",
    "$\\overline{x} = $ grand mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$\\overline{x}_1 = 109.727$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$\\overline{x}_2 = 107.706$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$\\overline{x}_3 = 106.735$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$\\overline{x} = 108.056$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$SS_{groups} = 186.374$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate means for each group\n",
    "mu1 = df['before'].values.mean()\n",
    "mu2 = df['during'].values.mean()\n",
    "mu3 = df['after'].values.mean()\n",
    "\n",
    "#Grand mean\n",
    "mu = df.values.mean()\n",
    "\n",
    "printmd('$\\overline{{x}}_1 = {}$'.format(round(mu1,3)))\n",
    "printmd('$\\overline{{x}}_2 = {}$'.format(round(mu2,3)))\n",
    "printmd('$\\overline{{x}}_3 = {}$'.format(round(mu3,3)))\n",
    "printmd('$\\overline{{x}} = {}$'.format(round(mu,3)))\n",
    "\n",
    "# Number of samples for each group:\n",
    "n = len(df)\n",
    "# Here the number of samples is the same for all groups:\n",
    "n1 = n\n",
    "n2 = n\n",
    "n3 = n\n",
    "\n",
    "# SS groups:\n",
    "SSgroups= n1*(mu1-mu)**2 + n2*(mu2-mu)**2 + n3*(mu3-mu)**2\n",
    "printmd(\"$SS_{{groups}} = {}$\".format(round(SSgroups,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sum of squares within groups $SS_\\textrm{w}$\n",
    "\n",
    "$SS_\\textrm{w} = \\sum_1 ( x_{i1} - \\overline{x}_1)^2  \n",
    "+ \\sum_2 ( x_{i2} - \\overline{x}_2)^2  \n",
    "+ \\dots + \\sum_k ( x_{i2k} - \\overline{x}_k)^2 $\n",
    "\n",
    "This might be more difficult to understand, so let's make this in steps.   \n",
    "Let's calculate the sum of squares for the first group   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "SS for group 1 =3903.5222013833727"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Group 1 is a column of the dataset (Before):   \n",
    "group1 = df['before'].values\n",
    "# Take mean value for this group\n",
    "mu1 = group1.mean() \n",
    "# Now calculate sum (x - mu1)^2 for all x values in this group:\n",
    "ss_group1 = np.sum((group1-mu1)**2)\n",
    "printmd('SS for group 1 ={}'.format(ss_group1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat the same operation for other groups, and the sum gives $SS_w$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "SS for group 2 =2110.5666278081503"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "SS for group 3 =2325.0030221455763"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$SS_w = 8339.092$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "group2 = df['during'].values\n",
    "# Take mean value for this group\n",
    "mu2 = group2.mean() \n",
    "# Now calculate sum (x - mu1)^2 for all x values in this group:\n",
    "ss_group2 = np.sum((group2-mu2)**2)\n",
    "printmd('SS for group 2 ={}'.format(ss_group2))\n",
    "#\n",
    "group3 = df['after'].values\n",
    "# Take mean value for this group\n",
    "mu3 = group3.mean() \n",
    "# Now calculate sum (x - mu1)^2 for all x values in this group:\n",
    "ss_group3 = np.sum((group3-mu3)**2)\n",
    "printmd('SS for group 3 ={}'.format(ss_group3))\n",
    "#\n",
    "SSw = ss_group1 + ss_group2 + ss_group3\n",
    "printmd('$SS_w = {}$'.format(round(SSw,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sum of squares for subjects $SS_\\textrm{subjects}$      \n",
    "We treat each subject as a level of an independent factor called subjects,   \n",
    "$SS_\\textrm{subjects} = k \\sum(\\overline{x}_i - \\overline{x})^2$   \n",
    "$\\overline{x}_i =$ mean of subject $i$    \n",
    "$\\overline{x} = $ grand mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>before</th>\n",
       "      <th>during</th>\n",
       "      <th>after</th>\n",
       "      <th>Subject means</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Subject</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>89.162573</td>\n",
       "      <td>94.023517</td>\n",
       "      <td>94.594145</td>\n",
       "      <td>92.593412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>90.857629</td>\n",
       "      <td>95.273755</td>\n",
       "      <td>95.040646</td>\n",
       "      <td>93.724010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>94.912999</td>\n",
       "      <td>96.612870</td>\n",
       "      <td>95.200472</td>\n",
       "      <td>95.575447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>95.254064</td>\n",
       "      <td>96.818673</td>\n",
       "      <td>97.205801</td>\n",
       "      <td>96.426179</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>97.136291</td>\n",
       "      <td>97.760342</td>\n",
       "      <td>98.428840</td>\n",
       "      <td>97.775158</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            before     during      after  Subject means\n",
       "Subject                                                \n",
       "0        89.162573  94.023517  94.594145      92.593412\n",
       "1        90.857629  95.273755  95.040646      93.724010\n",
       "2        94.912999  96.612870  95.200472      95.575447\n",
       "3        95.254064  96.818673  97.205801      96.426179\n",
       "4        97.136291  97.760342  98.428840      97.775158"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's calculate the subject means, i.e., the means for each row in our dataframe  \n",
    "\n",
    "def subject_mean(a,b,c):\n",
    "    '''Returns mean value of 3 numbers'''\n",
    "    return float((a+b+c)/3.0)\n",
    "\n",
    "df['Subject means'] = df.apply(lambda row: subject_mean(row['before'],row['during'],row['after']), axis = 1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$SS_{subjects} = 7996.418$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k=3.0\n",
    "subject_means = df['Subject means'].values\n",
    "SSsubjects = k * np.sum( (subject_means - mu)**2 )\n",
    "\n",
    "printmd('$SS_{{subjects}} = {}$'.format(round(SSsubjects,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sum of squares error $SS_\\textrm{error}$    \n",
    "\n",
    "$SS_\\textrm{error} = SS_\\textrm{w} - SS_\\textrm{subjects}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$SS_{error} = 342.674$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "SSerror = SSw - SSsubjects\n",
    "printmd('$SS_{{error}} = {}$'.format(round(SSerror,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Determine MS and the F-statistic\n",
    "\n",
    "##### Mean sum of squares for groups $MS_\\textrm{groups}$    \n",
    "$MS_\\textrm{groups} = \\frac{\\displaystyle SS_\\textrm{groups}}{k-1}$\n",
    "\n",
    "##### Mean sum of squares for error $MS_\\textrm{error}$   \n",
    "$MS_\\textrm{error} = \\frac{\\displaystyle SS_\\textrm{error}}{(n-1)(k-1)}$\n",
    "\n",
    "#### $F$-statistic:    \n",
    "As in other tests, we calculate an $F$-statistic   \n",
    "\n",
    "$F = \\frac{\\displaystyle MS_\\textrm{groups}}{\\displaystyle MS_\\textrm{error}}$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$MS_{groups} = 93.187$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$MS_{error} = 4.393$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$F$-statistic $= 21.211$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "DF1 = 2, DF2 = 78"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Degrees of freedom for groups:\n",
    "df1 = k-1\n",
    "# Degrees of freedom for error\n",
    "df2 = (n-1)*(k-1)\n",
    "\n",
    "MSgroups = SSgroups/df1\n",
    "\n",
    "MSerror = SSerror/df2\n",
    "\n",
    "Fstatistic = MSgroups/MSerror\n",
    "\n",
    "printmd('$MS_{{groups}} = {}$'.format(round(MSgroups,3)))\n",
    "printmd('$MS_{{error}} = {}$'.format(round(MSerror,3)))\n",
    "printmd('$F$-statistic $= {}$'.format(round(Fstatistic,3)))\n",
    "printmd('DF1 = {}, DF2 = {}'.format(int(df1),int(df2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, look at the table for the probability associated to the F-statistic with degrees of freedom DF1 and DF2, and significance of 0.05  \n",
    "\n",
    "\n",
    "\n",
    "In this case,\n",
    "$F(2,78) = 21.211$   \n",
    "From the [F-distribution calculator](http://stattrek.com/online-calculator/f-distribution.aspx), the    \n",
    "p-value for $F > 21.211 = 4\\times10^{-8}$, or virtually $0$\n",
    "Therefore, we can safely reject the null hypothesis.\n",
    "\n",
    "To my knowledge there is no scipy implementation of repeated measures ANOVA.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}