week05b-mult-regression-example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filename = 'data/week03_wcoa_cruise/WCOA2013_hy1.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>DATE_TIME</th>\n",
       "      <th>EXPOCODE</th>\n",
       "      <th>SECT_ID</th>\n",
       "      <th>LEG</th>\n",
       "      <th>LINE</th>\n",
       "      <th>STNNBR</th>\n",
       "      <th>CASTNO</th>\n",
       "      <th>BTLNBR</th>\n",
       "      <th>BTLNBR_FLAG_W</th>\n",
       "      <th>LATITUDE</th>\n",
       "      <th>...</th>\n",
       "      <th>TCARBN</th>\n",
       "      <th>TCARBN_FLAG_W</th>\n",
       "      <th>ALKALI</th>\n",
       "      <th>ALKALI_FLAG_W</th>\n",
       "      <th>PH_TOT</th>\n",
       "      <th>PH_TOT_FLAG_W</th>\n",
       "      <th>PH_TMP</th>\n",
       "      <th>CO32</th>\n",
       "      <th>CO32__FLAG_W</th>\n",
       "      <th>CHLORA</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2013-08-05 02:12:20</td>\n",
       "      <td>317W20130803</td>\n",
       "      <td>WCOA2013</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>48.2</td>\n",
       "      <td>...</td>\n",
       "      <td>2370.2</td>\n",
       "      <td>2</td>\n",
       "      <td>2369.0</td>\n",
       "      <td>2</td>\n",
       "      <td>7.294</td>\n",
       "      <td>2</td>\n",
       "      <td>25.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>9</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2013-08-05 02:12:53</td>\n",
       "      <td>317W20130803</td>\n",
       "      <td>WCOA2013</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>48.2</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>9</td>\n",
       "      <td>NaN</td>\n",
       "      <td>9</td>\n",
       "      <td>7.295</td>\n",
       "      <td>2</td>\n",
       "      <td>25.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>9</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2013-08-05 02:19:58</td>\n",
       "      <td>317W20130803</td>\n",
       "      <td>WCOA2013</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>48.2</td>\n",
       "      <td>...</td>\n",
       "      <td>2349.6</td>\n",
       "      <td>2</td>\n",
       "      <td>2343.7</td>\n",
       "      <td>2</td>\n",
       "      <td>7.282</td>\n",
       "      <td>2</td>\n",
       "      <td>25.0</td>\n",
       "      <td>43.521</td>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2013-08-05 02:27:01</td>\n",
       "      <td>317W20130803</td>\n",
       "      <td>WCOA2013</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>48.2</td>\n",
       "      <td>...</td>\n",
       "      <td>2318.7</td>\n",
       "      <td>2</td>\n",
       "      <td>2311.9</td>\n",
       "      <td>2</td>\n",
       "      <td>7.287</td>\n",
       "      <td>2</td>\n",
       "      <td>25.0</td>\n",
       "      <td>45.641</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2013-08-05 02:30:53</td>\n",
       "      <td>317W20130803</td>\n",
       "      <td>WCOA2013</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>48.2</td>\n",
       "      <td>...</td>\n",
       "      <td>2300.0</td>\n",
       "      <td>2</td>\n",
       "      <td>2299.7</td>\n",
       "      <td>2</td>\n",
       "      <td>7.308</td>\n",
       "      <td>2</td>\n",
       "      <td>25.0</td>\n",
       "      <td>47.741</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 42 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            DATE_TIME      EXPOCODE   SECT_ID  LEG LINE  STNNBR  CASTNO  \\\n",
       "0 2013-08-05 02:12:20  317W20130803  WCOA2013    1    2      11       1   \n",
       "1 2013-08-05 02:12:53  317W20130803  WCOA2013    1    2      11       1   \n",
       "2 2013-08-05 02:19:58  317W20130803  WCOA2013    1    2      11       1   \n",
       "3 2013-08-05 02:27:01  317W20130803  WCOA2013    1    2      11       1   \n",
       "4 2013-08-05 02:30:53  317W20130803  WCOA2013    1    2      11       1   \n",
       "\n",
       "   BTLNBR  BTLNBR_FLAG_W  LATITUDE   ...    TCARBN  TCARBN_FLAG_W  ALKALI  \\\n",
       "0       1              2      48.2   ...    2370.2              2  2369.0   \n",
       "1       2              2      48.2   ...       NaN              9     NaN   \n",
       "2       3              2      48.2   ...    2349.6              2  2343.7   \n",
       "3       4              2      48.2   ...    2318.7              2  2311.9   \n",
       "4       5              2      48.2   ...    2300.0              2  2299.7   \n",
       "\n",
       "   ALKALI_FLAG_W  PH_TOT  PH_TOT_FLAG_W  PH_TMP    CO32  CO32__FLAG_W  CHLORA  \n",
       "0              2   7.294              2    25.0     NaN             9     NaN  \n",
       "1              9   7.295              2    25.0     NaN             9     NaN  \n",
       "2              2   7.282              2    25.0  43.521             3     NaN  \n",
       "3              2   7.287              2    25.0  45.641             2     NaN  \n",
       "4              2   7.308              2    25.0  47.741             2     NaN  \n",
       "\n",
       "[5 rows x 42 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(filename,header=31,na_values=-999,parse_dates=[[8,9]])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple regression exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal is to predict pH based on the T,S,p,O2 using a standard multiple regression\n",
    "\n",
    "First create a subset of data:\n",
    "* Finite data (not NaN)\n",
    "* 30 > pressure > 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 464\n"
     ]
    }
   ],
   "source": [
    "# CTDOXY, CTDSAL, CTDPRS, CTDTMP\n",
    "subdf = df[ (np.isfinite(df['CTDOXY']) &\n",
    "             np.isfinite(df['CTDTMP']) & \n",
    "             np.isfinite(df['CTDPRS']) & \n",
    "             np.isfinite(df['CTDSAL']) & \n",
    "             np.isfinite(df['PH_TOT']) &\n",
    "             (df['CTDPRS'] >= 30) & \n",
    "             (df['CTDPRS'] <= 300 ) )]\n",
    "\n",
    "index, = np.where((df['CTDPRS'] > 30) & (df['CTDPRS'] < 300) & np.isfinite(df['CTDOXY'] + df['CTDTMP'] + df['CTDPRS'] + df['CTDSAL']))\n",
    "\n",
    "print('N =',len(index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[   1.      197.8       6.7236   33.904   149.1   ]\n",
      " [   1.      153.7       7.205    33.817   124.4   ]\n",
      " [   1.      128.5       7.3828   33.67    155.6   ]\n",
      " ..., \n",
      " [   1.       60.5      10.7555   33.647   169.6   ]\n",
      " [   1.       40.4      11.437    33.606   221.9   ]\n",
      " [   1.       30.2      12.4522   33.589   247.6   ]]\n"
     ]
    }
   ],
   "source": [
    "# vector of pH\n",
    "y = subdf['PH_TOT']\n",
    "\n",
    "# matrix of predictor variables\n",
    "pres = subdf['CTDPRS']\n",
    "temp = subdf['CTDTMP']\n",
    "sal = subdf['CTDSAL']\n",
    "oxy = subdf['CTDOXY']\n",
    "A = np.array([np.ones(len(pres)),pres,temp,sal,oxy])\n",
    "A = np.transpose(A)\n",
    "\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the unknown coefficients that provide the best least squares fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coeff = np.linalg.lstsq(A,y)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the modeled vs. observed values and the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-Squared: 0.966327583834\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAFkCAYAAACuFXjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl4lNX9///nYQsIIcomkU0WrbiAhn0RrcWtEpaASkRF\nEMSlVtHu1a+trfbX6kdt1YpFBFQaRQ2yFRVFZSeRCLiACgiyJhgEwpb1/P64M2RmMjOZmSwzmbwe\n15UL577Pfc+5L9T7nXPe532MtRYRERERf+pFugMiIiIS3RQsiIiISEAKFkRERCQgBQsiIiISkIIF\nERERCUjBgoiIiASkYEFEREQCUrAgIiIiASlYEBERkYAULIiIiEhAIQULxpjvjDElPn6eDXDNPcaY\nr4wxx40xm40xt1S+2yIiIlJTGoTYvjdQ3+3zRcD7wFxfjY0xdwGPAZOAT4F+wHRjzEFr7eLQuysi\nIiI1zVRmIyljzDPAz6215/o5vwpYaa39rduxJ4G+1tohYX+xiIiI1JiwcxaMMQ2BccCMAM3igJNe\nx04CfY0x9X20FxERkSgT6jSEu1FAAjA7QJv3gEnGmPnW2ixjTG/gdqAh0ArI9nWRMaYlcDWwg/LB\nhoiIiPjXGDgbeM9am1sVN6xMsDARWGKt3R+gzV+AM4E1xph6wH5gFvAboCTAdVcDcyrRNxERkbpu\nHPDfqrhRWMGCMaYjMBQYGaidtfYkzsjCFJygYR8wBciz1h4IcOkOgNdee43u3buH08WoMnXqVJ5+\n+ulId6PK6HmiVyw9C+h5olksPQvE1vNs3ryZm2++GUrfpVUh3JGFiThTCP8LprG1thjYC2CMGQss\nrOCSkwDdu3cnKSkpzC5Gj4SEhJh4Dhc9T/SKpWcBPU80i6Vngdh7nlJVNo0fcrBgjDHAbcAsa22J\n17nHgXbW2vGln88B+gLrgBbAA8AFwK2V67aIiIjUlHBGFoYCHYCZPs4llp5zqQ88CJwLFAIfAQOt\ntd+H8b0iIiISASEHC9bapXgWZnI/N8Hr8xYg5sZ1REQkduw+spujBUcj3Y2opr0hakBqamqku1Cl\n9DzRK5aeBfQ80SxWnmXd7nX0md6H/Z0DLeyTSlVwrC7GmCRg/fr162Mx4URERKLAa5teY9KCSfQ6\nqxfpN6RzZrMzI92lKpGVlUWvXr0Aellrs6rinhpZEBGROqXElvD7D37PLfNuYeyFY1l267KYCRSq\nS2WKMomIiNQqefl53DzvZhZ+vZAnrnyCBwc8iLPITwJRsCAiInXCjkM7GJ42nB2HdrAwdSHXnXtd\npLtUayhYEBGRmLdi5wpS5qbQPK45a25fwwVtLoh0l2oV5SyIiEhMm5E1g5+98jMuaH0B6yatU6AQ\nBgULIiISk4pKinjgvQeYtHASEy6ewPu3vE+r01pFulu1kqYhREQk5hw+eZixb49l6balPHvts9zT\n5x4lMlaCggUREYkp3+Z+y/DXh7P/6H6WjFvClV2vjHSXaj1NQ4iISMz4cPuH9HupHyW2hHWT1ilQ\nqCIKFkREJCb8O/PfXP3a1fRp14d1k9ZxbstzI92lmKFgQUREarXC4kLuXnw39/zvHn7R9xcsvmkx\npzc+PdLdiinKWRARkVrr4ImDXP/m9SzfuZz/DPsPk3tN9tkuOxtGj4Z9+yAxEdLToU2bGu5sLaZg\nQUREaqXNBzaTnJbMoZOH+OCWD7js7Mv8th09Glatcv55+3ZISYGVK2uoozFA0xAiIlLrLPl2Cf1n\n9Kdxg8ZkTs4MGCiAM6IQ6LMEpmBBRERqDWstT615imFpwxjSaQirb19N5zM6V3hdYmLgzxKYpiFE\nRKRWyC/K567FdzFzw0x+O+i3PHbFY9SvVz+oa9PTnakH95wFCZ6CBRERiXo5x3JIeSOFzL2ZvDLy\nFW7peUtI17dpoxyFylCwICIiUW1T9iaS05LJL8rn4/EfM6DDAJ/ttOKh+ihnQUREotY7W95h4IyB\ntGjSgozJGX4DBShb8bB9u/NnSkoNdjTGKVgQEZGoY63l8RWPM+qNUVzT7RpWTlhJx4SOAa/Riofq\no2BBRESiyonCE9w872b+uOyP/L8h/4+518+laaOmFV6nFQ/VRzkLIiISNfbl7WPkGyP5PPtz3hjz\nBjdccEPQ12rFQ/VRsCAiIlFh/d71jHh9BAArJqyg11m9QrpeKx6qj6YhREQk4uZ+OZdLZ15Ku+bt\nyJycGXKgINVLwYKIiFSr7GwYPBi6dnX+zMkpO1diS3jko0e48a0bSemewsfjPyYxXskG0UbTECIi\nUq38beJ0rOAY498Zz9ub3+bxKx7nd4N/hzEmsp0VnxQsiIhItfK1pHHX4V0Mf3043+Z+yzs3vsOI\n80ZEpnMSFAULIiJSrRITnREFl6bnraHP9FE0btCY1bevpseZPSLXOQmKchZERKRapafDoEHQpQuc\nM+YVvu5/Oee0PIeMyRkKFGoJBQsiIlKt2rSBT5YXM2bab/n2wvGM6zGOD275gDZNq2/jhkBJlRI6\nTUOIiEi1OpJ/hHHp4/jft//jqaue4v7+91d7IqO/pEoJT0gjC8aY74wxJT5+ng1wzThjzAZjzDFj\nzF5jzAxjTIvKd11ERKLd9h+3M3DGQJbvXM6i1EVMHTC1RlY8aJ+IqhXqNERvoK3bz5WABeb6amyM\nGQTMBqYD5wNjgL7Af8Lsr4iI1BKf7PiEvtP7crLoJGtvX8u151xbY9+tfSKqVkjBgrU211qb4/oB\nkoFt1toVfi7pD3xnrX3eWrvTWrsaeBEnYBARkVognPn/6eunM/TVofRs25N1k9bRvXX36u+oG/ek\nykGDtE9EZYWd4GiMaQiMA2YEaLYG6GCMubb0mjOB64HF4X6viIjULNf8//btzp8pKf7bFpUUcd+S\n+7hj0R3ckXQH7457l5antQRqNunQtU/Etm3On22qL5eyTqhMguMoIAFnmsEna+1qY8zNwBvGmMal\n37cA+EUlvldERGpQsPP/P574kRvfupFl3y3j+Z8/z9197vY475102K0bbN2qF3ltUJmlkxOBJdba\n/f4aGGPOB/4J/AlIAq4GOuNMRYiISBWrjt/eg5n//yb3G/rP6M+nez/l/VveLxcoQPkgIy8v8CiF\nRI+wRhaMMR2BocDICpr+DlhlrX2q9PMXxpi7gRXGmD9aa7MDXTx16lQSEhI8jqWmppKamhpOt0VE\nYl51LBlMT3fus2+fEyh4z/8v3baUG966gbbN2pIxOYNuLbr5vI93JUfQKoXKSktLIy0tzePY4cOH\nq/x7jLU29IuM+RMwGehgrS0J0O4toMBae5PbsQHASqCdv1EJY0wSsH79+vUkJSWF3D8Rkbqqa1fP\nF3KXLs68fVXLzoaU0Zavmz9Pbp/7uaLTlaTf9DoJjRP8XpOT40w95OWVHRs0SPUPqlpWVha9evUC\n6GWtzaqKe4Y8DWGcBbK3AbO8AwVjzOPGGPcchoXAaGPMncaYzqVLKf8JrAs0fSEiIuGpqSWDo8YU\nsvqMu8jtdy+svY+TLy8KGCiAk5uwdatWKdRG4UxDDAU6ADN9nEssPQeAtXa2MaYZcA/wJHAI+BBn\nekJERKpYRVMGVSH3eC6fXTgGWq+C+TPgs4ns7+KMNowe7fnd1pY/ppGE2iesaYjqpmkIEZHo4R4E\nJHT7koPXDGfPD0coem0efD8YcEYJoCxfItAxBQvVqzqmIbQ3hIiIBHQqafLcRdDrJk7bczZr7vyQ\n+z85m30NykYMBgzwvM5X8qISGmsnBQsiIhLQ3n0WBj4JV/4Wvk6mTdZr9H4yvtwIgfdqB1e+hK9j\nUrsoWBAREb/yi/I5NnQKnDUbVvwelv2VdgN958b7y5eo7hwKqX4KFkRExKfso9mMemMUh9pn0Xjh\naxR9No4mzeBFP2X1XCWWvSlHofarTAVHERGJURv2b6DP9D5sO/gdDV79hJOZ4ygqcmokTJkS6d5J\nTVOwICIiHtI3pzPo5UG0btqaju9mcvzbfh7nlaRY9yhYEBERAKy1/HX5Xxk9dzTDzh3GigkrOLij\nfbl2SlKsexQsiIgIJwpPcFP6TTz80cM8evmjvD76dU5reFq5wCA+XkmKdZGCBRGRCKuOnSJDsefI\nHobMGsL8LfN58/o3efiyh3Eq+zuBgXt5Zm0pXTdpNYSISIRV1U6RvsotB3qxZ2fDVRMy+fKiEZQU\n16ftslU8s+QShrhd52+Fg9QtChZERCLMO2Ew3ARC76CjW7fAIwFD7knjm6SJsP9ieGMe+462ZR9V\ns621xBZNQ4iIRFhV7RTpHWTk5Tkvfm8ltoSHlz3MNxfdBF9dD7M/gqNt/d5HRCMLIiIRVlU7RXqX\nW4byL/6jBUe5dd6tvLPlHTp983d2zvs1YMrdR8SdggURkQirqryA9HRn6iEvr+yY+4t/56GdDH99\nONt/3M78sfPpd0YyKTth9244eBBatID27bXaQcpTsCAiEiPatHFyFHyNUqz6fhWj3hhF00ZNWXP7\nGi5scyGg3AQJjoIFEZEY4muUYtaGWdyx8A4GdBjAvy9/iykprYNeMSECSnAUEYlZxSXF/Or9XzFh\n/gTG9xzP0luWMuXm1qxa5eQ2rFpVPgEy0jUfJDppZEFEJAYdPnmYm9Jv4t2t7/LM1c/wy36/xBhT\n4TLNqqr5ILFFwYKISIzZdnAbyWnJ7M3by/9u+h9Xd7v61DnvFRPeKx+qquaDxBZNQ4iIxJC3sz7i\nvP/ry9btRXRZto5Lml/tcd67fLP3yoeqqvkgsUUjCyIitZh7iWfT+0W2n/cL7K7L4c25fHbyDJKT\noWFDz9URgaYVqqrmg8QWBQsiIlEo2H0eRo+GVWuK4OqpcP5z1P/0Xor/9xSUOP9737gR8vOdtsHk\nIGgvCPFFwYKISBQKNtFwd+6PMO4GOPtjWDiNBl9MobjE/32VgyDhULAgIhKFfCUaeo82/H3GFrKT\nh4PJhVeXwo7L6dEHGjUqq8roGlVwUQ6ChEMJjiIiNSyYWga+Eg1dow3bt8Oq/e9x2av96dS+IUlZ\nGXSpdzmDBsGiRc4IRPv2TtnnggLn+rg43wmNIsHQyIKISA0LNMXgGj3YvRvi4z33axgwAMBCv3/B\n1Q/QaO81ZDycRvNfNi/3Hd4jE+3aKRdBwqeRBRGRSgin4mGgWgauQGLnTmdkoH175yVvLeTkFkDy\nHXDt/bDmAS7ZvIDmceUDBdASSKlaGlkQEamEcCoeBiqM5B1IrFvnBCHHzQGOjhoD7dfCOzOJ33Yb\n87b6/w73JZAtW0JhoRPQaD8ICYdGFkREguRrFCGcioeBCiN5jwAUFcGqbz/ns959odUWmL0MNtxG\n69aBX/iuJZDbtjkJjxkZ/veDEKmIRhZERILkaxShovLJvgSqZeAaEVi3zgkU+MkCSBmHOdQVO/Nj\nONwp6O9xUQlnqSwFCyIiQfL10l2zpvIVD30VYBqVYlnNP2Do72HLSC7Z+QpNLmwW1veEE9CIuFOw\nICISJF8v3aqoeOg9YtGh80nib5oM7V/j9E0PcX72n5k3v17YeQYq4SyVpWBBRCRIob50gy3Z7DFi\n0Ww/BWNHknvmRs7dlMbXb4+tdL9VwlkqS8GCiEiQgn3puoKETz/13JehWzdo3bp84HBqxCIxC8aO\ngHrFMHM5RY37VNuziIQipNUQxpjvjDElPn6e9dN+Zun5Yq/2n1dN90VEoo9rWsG71HJenu8VCenp\n0KDHWzBxMBxtC//JhL19lFsgUSPUpZO9gbZuP1cCFpjrp/0vS9sllv7ZHjgYoL2ISK0XzGoDV5v9\n+y29f/UoRSnXw5YRMHM55LUjPl65BRI9QgoWrLW51toc1w+QDGyz1q7w0z7Pq31f4HRgVmU7LiIS\nrbxHBOLinNLN3m2OFx7nwj+PZVfXR+DDv8Lb/yWufhMGDYKtW1U4SaJH2DkLxpiGwDjgyRAumwh8\nYK3dFe73iohEO3+JkO7HnntlN5fOHMHBllvg9XTYMsrjegUKEk0qk+A4CkgAZgfT2BiTCFwLVD61\nV0QkivlLhHQdW7d7Hde+MZKG9RrS49NVbNxy8ak2+fnBlYwWqUmVCRYmAkustfuDbH8b8CMwP9gv\nmDp1KgkJCR7HUlNTSU1NDfYWIiJRZc6mOdy+4HZ6ndWL9BvSMePOpGNHz2RIVViUYKWlpZGWluZx\n7PDhw1X+PcZaG/pFxnQEtgMjrbWLgrzmG2CBtfZXQbRNAtavX7+epKSkkPsnIhJtSmwJDy17iL+t\n/Bvje47nxWEvEtcgDnD2mXAVZQJnvwiNLEi4srKy6NWrF0Ava21WVdwz3JGFiUA28L9gGhtjLge6\nAjPC/D4RkVorLz+PW+bdwoKvF/DElU/w4IAHMcacOq8KixLtQg4WjPNv+G3ALGttide5x4F21trx\nXpfdDqyz1m4Ot6MiIlUl2MqKVXGvHYd2MDxtODsO7WBh6kKuO/e6cvdQhUWJduGMLAwFOgAzfZxL\nLD13ijGmOU4y5C/D+C4RkSrna/dI75d1sAFFuX0dOjhLJVu0gPgLV7Dv0hQSmjRnze1ruKDNBdX7\nYCLVJORgwVq7FKjv59wEH8eOAM1C75qISPUIZsvmYAIKX9cWFDg/ed1mQNJdNN87iHV/e5NWp7Wq\nms6LRECoFRxFRGo976JJBw5A165OomFOjnMsmIDC172oVwRXPwAjJsFnEzk5/X1KjipQkNpNwYKI\n1Dnp6c6Kgy5dnMqKvvZs8A4C/O3T4LpXXBwQdxhSk6Hfv2Dxc7DoBQpONPTYBwKcKY7Bg8sHKCLR\nSsGCiMQ875czOFMK27Y5u0C6c40guAcUgwb5X6HgSk78eNO3NLm3P/U6roXX3oXMewDjcU8X1xSH\nr02lRKKRggURiXmBXs7+RhDatIG333Y+79vnXONvBODD7R/y83n96NjJsvnBdQxKHOrzni4VTXFo\n5EGijYIFEYl5gV7OgUYQghkB+Hfmv7n6tavp064Payet5dyW51Y4KlHRFIdGHiTaVKbcs4hIrZCY\n6Lx43T+7BKpxECjIKCwu5L537+OFT1/gvn738eRVT9KgXoOgllxWVIQp2ORKkZqiYEFEYt60aTBw\nIJw4AU2awIsvBnedvyDj4ImDXP/m9SzfuZz/DPsPk3tNPtUmmCWXFRVhChTciESCggURiXl33ums\neADnzylTAr+sXaMDu3c7qyVatHBe8IWF0CFpMweGJlPc8BCtP/qA2e9dxgi30YOqGBVQ+WeJNgoW\nRCTmhfoCdx8dAOjRw/lzVfYSGDMWjnSA/2ay71Bn9uE5elAVowIq/yzRRgmOIhKzXKsK9uzxPF7R\nC9w7mNi7z7L59KfhpmGwcwjMWA2HOp86v3t3Wdtgl1yK1CYaWRCRmOU9QhAXB717OzkMgweXDfNP\nm+ZMVbg+t2rlNjpQP59jP7uLg+1mwsrfwoePgfWseH/wYNk/a1RAYpGCBRGJWe6/8QO0beu8yPv1\ng4wM59j27TBgABw9Wva5aVPo2xeyj+Zw8MoUDrXM5LnLXyFt6S182gjy8z3v26JF9T+LSCQpWBCR\nmOX+G7/7540bPY+7AgWXY8fgiwObaHFPMqfZfN4f+wn92/fnnpVOgaRu3coSJgHat6/6votEEwUL\nIhKzWrTwfKkHPQJw3jscT7kZ9pzLlkfm0yGhw6lTbdrA1q1arSB1i4IFEYlZ7dvDzp1ln88808lV\nKCnxd4WFS/8GP/sjfDWa4/Nmk5rRtFxhJeUlSF2j1RAiErO8VyZY6yQ8Fhb6aNzgBKTc7AQKHz8C\nb86FwqYqtyyCRhZEpBarqLSy9whA165+btRsH4wdCWd+Dm++AV/e4HFa5ZalrtPIgojUWsOHe264\nlJwcuL3P+gqJ6+GOPtB8Dz9ZtYK+TW8gLi6I60TqEI0siEit5b2qwfuzN1cZ5d27nZURcUlzOTjk\nNnomXsTim98hMd6JCnJylMAo4k7BgojUGa5piRJbwp8//jOPLn+UcReNY3rydJo0bFKunYg4NA0h\nIlHPVba5a1fnz5wc57hrzwYX78++7tH5J8doc88NPLr8UR6/4nFeHfWqR6AgIuVpZEFEop6/bZ8X\nLQp+umD0aFj1+S5IHQ5nfMt5G97h94+MqJkHEKnlFCyISNTzt2tkRdMF7qsldrEG7hgFRY1hxmoK\nmgYYhhARD5qGEJGo570aIdjVCa4Rie3NXqXwpssh9xyYngHZPbTCQSQEChZEJOqFu+3z3v3FMPR3\nkHIrbLoZXv2ATq3a0LcvFBSUz4EQEd8ULIhI1HIlJQ4Y4LzcW7VyphRSUip+wW/ddYTvB46EgU/A\nu0/DgpegKI727aFhQ8jMLKvPoAqNIoEpZ0FEIs5fJUbvxEYX9yRHX/c4o8t2vuo5nOJ2u+C/i2Hr\nNafa+arGqAqNIoEpWBCRiPO32iHQS9z73Kl7dPoEkkZjTp4Bs9fCD9092rlyFdyDD+UviASmYEFE\nIs7XaofsbDhwwP813i/4ffuApOlw3d2wcwgN579JweGyPanj4qB377J8B1VoFAmeggURibjExPK/\n6Y8eDXl5ZceaNoULLoAffij/gi8qKeLE5Q9Cx39Bxt3w7jNc3LshDRv632RKFRpFgqdgQUQizrVn\ng/uLfcAAzzZnngnr1pW/9scTP3LjWzeS02kZXb56HrbcTeKA8sGBiIRPwYKIRJyv4kq+Rhu8fZP7\nDclpyRw4doD3b3mfKzpfUb0dFamjtHRSRKJSRbUVlm5bSr+X+lHP1CNjcgZXdL7C7x4SIlI5IQUL\nxpjvjDElPn6eDXBNI2PMY8aYHcaYk8aY7caY2yrdcxGJepV5ebtGG7Ztc/50TSlYa3ku4zmunXMt\n/dv3Z0HyWm4b3o2uXeGcc0orNqp+gkiVCnUaojdQ3+3zRcD7wNwA17wJtAYmANuARDSiIVJr+KuB\nEAx/SyLDVVhcyL1L7uXF9S/yQP8H+MeV/+CyIfVPfYc31U8QqRohBQvW2lz3z8aYZGCbtXaFr/bG\nmGuAS4Eu1tpDpYe/D6ejIhIZlXnh+9sAKhy5x3MZ8+YYVn2/ihnDZzDxkokV3lP1E0SqRti/4Rtj\nGgLjgBkBmiUDnwK/NcbsNsZ8bYx5whjTONzvFZGaVZkXfrgbQHn7MudL+r7Uly9yvmDZ+GWnAgVf\n94yPD30PCREJrDKrIUYBCcDsAG264IwsnARGAq2AF4AWwO2V+G4RqSHBrErwx9eSSG8VTXMs/mYx\nqW+ncvbpZ/PhrR9y9ulnV/gdWjIpUrWMtTa8C415F8i31o4I0OY9YDBwprX2aOmxUTh5DE2ttfl+\nrksC1g8ZMoSEhASPc6mpqaSmpobVZxEJXU5OxS/jYPIa/LXp1w8yMsra9e3r1FPYv98y4MH/Y8c5\nv6HBtuEkrnmNjm2bKRgQcZOWlkZaWprHscOHD7N8+XKAXtbarKr4nrCCBWNMR2A7MNJauyhAu1nA\nQGvtuW7HzgO+BM611m7zc10SsH79+vUkJSWF3D8RCY/3C33aNLjzzop/ax88GI8kw0GDyuc1eLeJ\nj4fWrT1HLQAaNYIjx/LpeM8Ucs6aDcv/AB/9BWw9v/cWkTJZWVn06tULqjBYCHcaYiKQDfyvgnar\ngDHGmNOstcdLj/0EKAF2h/ndIlJNvJMZBw4sK7kcKLnRO4/h00+d5ZLuAYZ3m7w8z3LOLgUNs2n9\n61HktcmCt+fA5zcF/C4RqX4hJzgaYwxwGzDLWlvide5xY4x7DsN/gVxgpjGmuzFmCPAPYIa/KQgR\niRzvF/GJE4HPu3jnMeTnl691EFSuQ9sNcEcf8up/R5PXl5cLFIK+j4hUqXBWQwwFOgAzfZxLLD0H\ngLX2GHAlcDqQCbwKzAfuC+N7RaSaeb+ImzQJfN7FvdpiXJznOVeA4d4mPt7HTbqnw8RBcKwNTM+k\n+Pu+DBoEnTo57Tt10goHkUgJeRrCWrsUz8JM7ucm+Dj2DXB16F0TkZrmvbLgxRdhypSKt3J239vB\nOzfBFWC4t8nJgWHDYNMmKLEWLn2Mwksfhi9ugPkzofA0TJxyE0SihTaSEpFTfG3oFOoLO5jlktbC\nli2QX3wCRkyEi16HZY/C8ocAA0CPHuE9g4hUPQULIhI2f8shKwowkpMhjz0wYSS0+RLmvglfjTl1\nPi4OFvldZyUiNU3BgoiELdxS0BsOZMIdI6CkPsxYBfsv8Tjfu7dqKYhEE23oJCJhC6cUdNrnaRSO\nGwKHOsH0TI9AIS4uupIYteW1iEPBgoiELdi9H7KzYdDgEs4Y/TA3pd/EGfuvh9kfwdG25dq6b0cd\naa6RE215LXWdggURCZv7cshBg5yKj75+Ex95w1FWtx/DoYseg6V/p9vnsxnUr/x+cgUFNfwAFajK\nXTNFajPlLIjEqGD2a6gs92TG7Gw45xzPio/dusHpZ+9k9+ARcPo2SJsP3yST2wW2bXOmHdwDhIYN\nq7Z/lVWZTbREYomCBZEYFW7yYTB8BSKjR5cv35x3+mryrh4FhU1hxhrIuRCAVq2cDaMKCz3b9+xZ\nNf2rKsEsAxWpCxQsiMSo6hxC9xWIlLv/xbNg2BTY3R/mvk1ccSvadXFeugUFkJnp2Tw+PvqWSwaz\nDFSkLlDOgkiMCjb5MBy+ApFT9zfFcNWvYOQE2HgrvLoUjreid29n6mHlSsjNLX/P1q2jJ7FRRDwp\nWBCJUd7Jh1U5hO4rEElPh36XHcaMGw79n4Yl/4SF/6FR/Ublvt9X4KJ8AJHopWkIkRhVnUPovuby\n8xps48j1ydjde2HOEth2FeDkJXgnV6anl+0NAU6ugvIBRKKXggUR8SnQagrvQOSj7z5izJtjaNmk\nJQ1nraNw/09OnbO2fHJlmzaQkVEzzyEiladpCJEYF24VwmAKEmVnQ7cbX+SKWVdRvDuJRSPXEXf0\nJ+XaqT6BSO2mkQWRGBfuEsqKVlMUlRRxyR+msu/852DdvRx+7ykmZjWgZUs4etSzrfIRRGo3jSyI\nxDjvl/yuXcGNNHi/4A8cKLvmm10/cu2ca9nXfhosnAZL/gUlDdi3D9q397wuPl75CCK1nUYWRGKc\ndxXCH3+wKZfbAAAgAElEQVSE7793/jnQSIN7EuOBA07Bpbw82H5kCz2fHc5prXK5MGspX6y/3OO7\nfCU/akmkSO2mYEEkxnm/vHfv9qy0uHo1nH22MyLgL4mxa9fSa7q+B9ffSPHJdmRMyiD+tq4+AwMV\nMhKJLZqGEIlRrsTGAQOcz2vWOC9x72kCa2HnzsC7KrZNtNDvnzDu5/D9IJKy1tC1RVfatIG333YC\nhX37nOtd0xra3lkkdihYEIlR/lYzpKc7Gzj54mvVQkFxAZ1/eQdcez8Jmx9g4M4FLHizeYXfo+2d\nRWKHpiFEYpSv1Qyu2gn+7N4NzZtDixbOCMT0OT9w50ejWbt7LbNGzGL8I+OD+p5Ax0Wk9lGwIBKj\nfG2v7L6MEpwRhkaNID/f2dzJ9ZOXBztPfMElLyST0Oo4H43/iIEdBgb9PYGOi0jto2kIkRjla28I\nX7/dt24NxngdPHch3D4AeyKBzMmZDOww0G8Ogr89KKpzbwoRqVkaWRCJUb5WJXj/tp+f7/kZLAz6\nBwz9PWwZSdK+V+iY0AzwX9zJ3+qH6lwV4a8UdaAS1SISPo0siNQh7r/tl0tybHASRt0KV/6O+qse\nYuDut5j/ZrNTp6MpB0FJlSI1SyMLInWI+2/7gwe75S802w9jR8KZG+GtNDodH8uqbZ7XRlMOgpIq\nRWqWggWRGBHqEHx6OnTrBnnNsmDsCKhXDDOXw94+ZDeFxo2dduefDw0aODkK8fFlKyUimYOgpEqR\nmqVgQaQW8hUYhLJhlOv6Y53eghG3Qs4F8Po7kNeO+vXh2LGytp995nltjx6Rr9Doq6R0oOMiUjkK\nFkRqIV+BQShD8MnDLZlN/gJjHoHPx8L8l6GoCeBjZYSXaBjaj0RSpUhdpgRHkVrIV2DgPeTubwj+\neOFx1p89Fn76CHz4V3j7v6cCBSibfvBHQ/sidY9GFkRqIV9z84GG4F3TDrsO7yb3qhGUdNsCr6fD\nllGn2sTFQe/ezhTEhg3+v/v4cSd/QUsSReoOBQsiUco9L6FVK2fDp9xcJxB48UWYMiX43R6TkyFz\n7zpnxUNJQ3h5Fey/2KNNu3bO9V27el5rjPPdLp99FjgfQkRij4IFkSjlnZfgsn079OrljAKsWRPc\nb/ifFc+BCbfD3l7wRjocO7NcG38rClzloN1FQ96CiNSckHIWjDHfGWNKfPw866f9ZT7aFhtjNIAp\nUoFAL+T8fCeQGDYs8DbQJbaEP3z4B4qG3wyfp8LsZT4Dhbg4/2Wae/Ys//3KWxCpW0IdWegN1Hf7\nfBHwPjA3wDUWOBfIO3XAWu1sL1KBVq28SzGXt2lT2W/93ssl8/LzuGXeLSz4egEdNj/JrvkPAL6X\nOvTuXTZC4T2dkZPjTGNs3Oh87tFDSxJF6pqQggVrba77Z2NMMrDNWruigksPWGuPhNo5kdqmKvcm\ncM8TAGjaFIqKyk8JuFu3zhlheGbWDm5dMpxvcnZw5ieLaHvy57S6BL76ymnnKrTkyoEI9PJv08a5\nr4jUXWHnLBhjGgLjgCcragpsMMY0Br4A/mStXR3u94pEs1AKI/njCji8iyGdeaaTo+C+4qGwEDIy\nytoUFcGqXSsZMHMU9QubUzxrLfsPnM9+oG9fOHmyUo8nInVUZRIcRwEJwOwAbfYBU4BPgThgMvCx\nMaavtTbA4iyR2qkq9iZwDzjcJSb6niJISXF+8y8qAi6ZAcPuon7OIPJffQuOtzzVNjPTyW3Qbowi\nEqrKBAsTgSXW2v3+GlhrvwG+cTu01hjTFZgKjK/Ed4tEparYm8A7wGjQAPr18z1V4AoeBl1axOqm\nv4EBT8OnU2i44lnyjzf0aGut07ft2509IbZuVcAgIsEJK1gwxnQEhgIjw7g8AxgUTMOpU6eSkJDg\ncSw1NZXU1NQwvlak+vkqjBRqHoN3wNGvX+CpjMMnD9N4wljYuZSW657jJ0fuZs/phqOH/V+Tl6da\nCSKxIC0tjbS0NI9jhw8H+I8/TMZ6Z1EFc5Exf8KZUuhgrS0J8dr3gSPW2jEB2iQB69evX09SUlLI\n/ROpTqG+/D22gsZZjrhypf/7uKYW/N3f/brTu35L3nXDOXByP29e/yZDuwz1+Z3160NxsWe/unSB\nbV7bUItI7ZeVlUWvXr0Aellrs6riniGPLBhjDHAbMMs7UDDGPA60s9aOL/18H/Ad8CXQGCfA+Clw\nZeW6LRI5oSYx+stjSE528ghc9+na1QkKKgpATn1/5w8h6Xqa7GnDht+u49yW555qM20aDBwIJ05A\nkyawaJFTkyEvr+w+qpUgIsEKZyOpoUAHYKaPc4ml51waAf8HbAI+xqnL8DNr7cdhfK9IVAg1idHf\nBk+bNnkeP3rUCRpWrXICEH927wb6/BtuuRr29KHVO2s9AgWA2293AoOiIufPX//ayVFwL7akWgki\nEqyQRxastUvxLMzkfm6C1+cngCfC65pIzQtmiiHUJMZAGzz5s2uXM5Xg3Y/C4kL2XnIfXPwCrL0P\n3n+SvaZBuVUOrgJKLhs3avtmEQmf9oYQcRPMFEOoL39/L+mePT1rJLj7/nvnx9WPbt0g84uD3P3x\n9RRetAIWTIesSQAUU7bKITlZBZREpOopWBBxE8wUQ0W/oQebALlwYVnQ0bIlbN7sTEX4khe3mYuf\nT6ZZq8Ocl/EBW7KG+GznXpLZlQ/h+iwiEq5wchZEYpa//IJQuEYnKso/cAUdq1c7Ozv6ra7YbQlM\n6k9xfhMyJmXQ/EffgYK7RYs88xMWLQr9OSQ6ZWcH3jxMpDpoZEHETTj5Bd4qGp3wHnk4dgw2+Kxn\naqH/M3DVr+Db62i4ZA49X473WNFQ7grrvEDS05WfEKuqoqS4SKgULIi4qYokQO/dIlu18jzv/T97\n42sjyPr5cN3dkPQyrPgdLPsrJ6hfbnMpbwUFZaMZeoHEpqooKS4SKgULIlXM+4Xu/dn7f+7lAoCm\nOXDDaGiXAemvwKZbnHY+vis+Hlq0gP37PXej1AskdlVFSXGRUClnQSRM/uaOc3M923l/9v6fez33\n/wrP3AST+0LLb2HWJ6cCBV+aNYMjR2DHDujdO/B3SOxIT1e9DKl5GlkQCZO/ueOKfvNLT3eWQrpy\nD0pKnBGC05Lmc+DScZB7LiUvz4cjHQikWzfPe1Y210JqB9XLkEhQsCASJO/ExN27Pc+7hv4renG3\naQOtW7uXXrY0+On/R07SH2mRncJp78/mx5KmHDM+pijcfPFFWTKjXiAiUp00DSESJO8lkQcPep53\njSC0aQNvv+183rfPCRy8l7edGm1ocAJSbubHpD9Qb8X/I3faXHZta8rRo56BQrNm0LSp5z2Kiiou\nDS0iUhU0siASJO+kwRYtnGJHvkYQKlrelp4Ow8buY9P5I8k//XN48w2Kv7zB73cXFsL69TBlilOh\nsajIf79ERKqaggURP7ynHbyXRLZv73/ov6LlbRty1rOx3wiKi6Hhayso3NkrYF/y851AYeXK8ttP\nHzhAub0hRESqkqYhRPzwnnawtnwWumtFRIcO0KABNGwIzZs7P+5cxxo2hCa953LNG5dS8EM7il/I\nrDBQcHHPiXD1Iz7eyX0IZrdKEZFwaWRBxA/v0YDcXNi2zfOY92/54Ly8XdtBu0YlNm6Eo8dK4LI/\nU3T5o7BpnLMZVFGToPvjnhPhGtHo2hWPio6akhCR6qBgQcSPYIrf+Hs5Hz1aFihMmwYXJR2D68dD\n93T44G+w8rdAWelG47XyIS7O2ZXSWidI8bccUgV6RKQmKFgQ8SOY2gXeL2t3rm2j+1+9CyYOdwot\nvT4Pvh5Rrq0rUDDGWfXQsqUzZVFRDoLqK4hITVCwIOJHMFtRFxQ4owDWOisU6tWD4mK3UYL2azg2\ndhQUNYYZqyE78F7R1jqjEkePws6dkJzsrH4It48iIlVBCY4ipbzLN7uKHvnbCnj0aMjMdFYqFBTA\ngAHOEseBA0sb9HgVbrscDp4D/8k8FSg0axZ8nzIztQWxiESeRhZESnnXRhgwwPkN3/XZ+7d873yF\nVaugb1946eVirv77H9nf7e+QNREWvwDFjYCyjZ9c962ItdpBUkQiTyMLIqW8X/7eL/SMDM8Rhz17\nyt8jc2Mel08bxf6uT8C7T8OCl04FCgAJCeXLRIfaLxGRmqaRBZFS3smK3isUwBk9GDjQc7niKad/\nB6nDOZTwPW0/XMz+tdeUa7Jvn5PTEGq/REQiSSMLIqXS051phLg458cfn4FCp+VwRx9ocJKLP11L\nV1s+UGja1AlAghEfry2IRSR6aGRB6hzvMs7p6c4IwujRTvGk/PwQb5g0Ha67G3YO4bTFb/LuVy0A\n6NwZjh8va2ZtcKMK8fFOUSeVbRaRaKFgQeocX5s8QflKjBWp37AIrnqQ4j7/goy74d1naNW+IcOG\nwaZN5YMO98DBW1wctGun/R1EJDopWJBaz3ukYNo0uPNOz5ED95dvRZs8uYuPh9atnfscOwYbNpSe\naPwjxdffCJ2XweLnIfNuwKm2+P33FffZOx+id2+teBCR6KVgQWo975EC9wREX9tD+yuR7H4sLs55\ngbsHGv36lZ5s+Q2kJkPTA/Dq+/DdFaeuc986OpDevaFRI1VeFJHaQcGC1HreIwMnTgQ+76tEck6O\nE2ScOAFNmsCaNXDBBU5718hFVhbQZSlcfwMcbQvTM+BgN497FxRU3N/4eFi0SFMNIlJ7aDWE1Hre\nSwubNAl83ns5JDjTFnl5zshAXh5MmVJ2zhm5sBQlPQc3Xwu7B8BLa8sFCv7u7a11awUKIlK7KFiQ\nWi893Vli6FpquGaN52fvIX7XtMX27c6fKSnlCyW5f967vxCG3QU/vxfW3g//XQj5CX77U9HyyD17\nfJePFhGJVpqGkFrP12ZKgZIFfSU4HjzoeWz3bmdPiFYdc/nh52Pg9FXwzsuwYUKF/bHW2VCqpKTs\nWHy8M0WRn+/8uIIUJTWKSG2gkQWpc7ynJfbsKZ9rUFwM2/O+JKNHX/Iaf0mTN5cRt7niQMGlpMSz\nsNLWrc7SSHcq4ywitYWCBalzXNMWriqNrt/2PZyzGCYNgMKmMD2DE1sGh/w9rVvDtm3O6EGbNuWD\nFJVxFpHaQsGC1DmuaYu2bT2Px8VBozgLA5+Em5KdJZEzVsOhs33ex5jA+QnewYB3boWWS4pIbRFS\nzoIx5jugk49Tz1tr763g2kHAx8Dn1tqkUL5XpDp45yk0bJKPvW4KBefMhuV/gI/+ArYsnu7RwwkO\nNm50pi3cVz7ExUGDBmV1Fnr2LB8M+MqtEBGpDUIdWegNtHX7uRKwwNxAFxljEoDZwAdh9FEkbNnZ\nTjGlxo2dn759nVUI2dleUw9NsylIvYJjnV+Ht+fAssdOBQqNGkGfPk5thIYNneu8l0i2a+dsaX3y\npPOzbp2WR4pI7AhpZMFam+v+2RiTDGyz1q6o4NJpwBygBBgRUg9FKiE5GTIzyz5nZkLHjk4AcCqp\n8cyNkDqcgvoFMHM57OnrcY+CAqd9mzb+kxKVfyAisSzsnAVjTENgHDCjgnYTgM7An8P9LpFwbdpU\n/lh+vts20+fNg9sHwvHWMD2zXKDg4goSvIOCuDjlH4hI7KtMnYVRgGt6wSdjzDnA48Bga22Jqaha\njUgleW8q5Z+FIY/BFQ/DFzfQaMlMCo6d5re1616+SkVrukFEYl1lgoWJwBJr7X5fJ40x9XCmHh6x\n1m5zHa7E94lUyHtTqWbNfCyLbHACRkyEi16n3sePctr6hzi7s+GLL3zfs2nTspEDJSmKSF0UVrBg\njOkIDAVGBmgWj5MQebEx5vnSY/Wcy00BcJW19uNA3zN16lQSEjzL6qamppKamhpOt6UO8M4pKChw\npgpOrV6I3wNjR0Lrr+CNtyjZPJqjwI4dntfVqwdnnQWHD0OLFs5ogkYRRCTapKWlkZaW5nHs8OHD\nVf49xgaz8433Rcb8CZgMdLDWlvhpY4DuXofvAX4KjAZ2WGtPlLvQuTYJWL9+/XqSkrTKUoI3eHDZ\nyEI5Z2VC6ggoaQBp82H/JX7v06WLM83gfq9BgzSqICLRLysri169egH0stZmVcU9Qx5ZKA0CbgNm\neQcKxpjHgXbW2vHWiUK+8jqfA5y01m4Ov8si/rnnFOzZ4zYFcWGaM/Ww/2J4Y56zxXQAiYm+95AQ\nEamLwlkNMRToAMz0cS6x9JxIRLhyClavdpY7Ykrgpw/DmJvgyxtg9kcVBgrgBB0qzywi4gh5ZMFa\nuxSo7+dcwJ12rLV/RksopQp88QUMHAgnTjibNtWv7+QZdO/uBAkbN0K+PQo33ArnvQPv/wNW/4pg\ncmzj452gw9fKBxGRukhbVEutNHCgW60EyraD3rCh9EDCTic/4YxtkLYAvhnm8z7GeFZjNAbWrHH+\nWSsfREQc2khKolZ2tpOw2LWr82dOTtm5Ez5TY0t1WA2T+0LcEZixxm+gAKVTFW4GDoQLLqhcv0VE\nYo1GFiRivAsoeS9N9K6ZkJJS9pt+kyaeIwunXDwLhk2B3f1h7ttwvFXAPvTo4QQMmmoQEfFPwYJE\nTKBgAPyvRsjOdkYbNm50m0IwxTD0dzDoSciaBIufh2KvYYNS9etD+/bOj2oniIhUTMGCRExFSxMT\nE50gwv0zOEHGqdwEcKYbRqdCt3dhyT9h3b0ESmQsLnYCBeUjiIgER8GCRIy/YMDF32oEj6DijG2Q\nOhya74E5S2DbVadOGeNsKV1UVJYA6aKaCSIiwVOCo0RMerpTFbFLF987N7ZpA2+/XVYgKSXFSXI8\nFVSc/ZGTyFi/EKav8wgUwJmiKCgoHygAHDjgmTApIiL+aWRBIsbX0kT3+glNmkC3bvDZZ845V15D\nejoMvPdFtp33C5rnXk671XPZcvAMQilcnpcHycnOyIN2kBQRCUzBgkQV9/oJeXleuQnA3v1F/H7F\nVLad/xzNN9/LhXueYuM3DQhjixOncFNpOWhfCZYiIuJQsCAR4W/ZpHf9BI8goPGPHPz5Dczc9DEs\nfpEj6+9gdRX2SXkMIiK+KViQiPC3bNK7fkLTpnDxxbDz2BZ+uHI49Vrk0nbBUvatv/xUG+8qjO7i\n4z3vFx8PrVs7AUpBAWRmlp3T3g8iIr4pwVEiYudO35/XrHFe6A0aOH+uWwcPz3qPvBv70+Xshnw6\nJZMu9S73uPbii53aCe7i4pykyTVrPJMot26FbducwGTRosAJliIi4tDIgkSE95D/7t3QuDH07Ons\nGHnnnbB3nyX5r/9iZ/cHuLbbtfx39H9pHtfc75JK72OuZEV/eQja+0FEJDgKFiRq5OdDRoYzUlBM\nAfz8HjjvJdp992vmP/Q36tdzhg/8veT14hcRqR4KFiTqFMf9ADeMhvZrYd4sGhwaT/3SCbOK9pMQ\nEZGqp2BBIiIx0Zl6KKfNF5CaDA2Pw+yPYNdADsaXna5oPwkREal6SnCUiOjUycfBcxfC7QMgPwGm\nZ8KugQC0aFHWpKL9JEREpOopWJCImDbNWe3gsDDoH5A6ArZfCS+vhMMdT7U9eNDZZXLwYGjlteO0\nljuKiFQ/BQtS47Kz3So1NjgJo8bDlb+F5Q9R7623aNaoGZdc4ixpdNVJ2L7dmX6wFvr2dZZGxsU5\ntRK0x4OISPVSzoLUuOHDSwOFZvvhxlHQdgNtlqeRvWxsubZdu3oWVdq0yfnTVaY5M1N5CyIi1U3B\ngtS4jRuBtp85W0vXK4aZyzmncx+PNq5VD3v2eF7rChLcKW9BRKR6aRpCalzxT96C2wfB0UT4z6ew\ntw+FhZ7TCa5VD67gwBhn2sEX5S2IiFQvBQtSY6y1PPrJoxSlXA9bRsLMTyDvLMApxtStW1nA4D1a\nYC00alT+nvHxKtMsIlLdNA0hNeJ44XEmzJ/A3C/n8vu+j/HJ0t+TgaHIrU1eHgwb5gQF3tMP4Cyh\nLCjwnIpo3VpFmUREqpuCBal2u4/sZsTrI9hy4GvO25DOG6+MIjERkpKcEQV3mzb5zksAaN/e+XEV\nZQJNQYiI1AQFC1Kt1u1ex8g3RtKofiPOXbGKDe/1BJylkH36lN9CuqDA8/q4OGjXLvCGUSIiUr2U\nsyB+ZWc7hZBcBZG++MIZDahXz0k4jI+HL7/0f/2cTXO4bNZldD2jK5mTMznybU+P87m5zpbR8W7l\nnK31vEfv3mVbSlvrf2dJERGpPhpZEL+892E4VUip1NGjMGAAHDnieV2JLWHqwof412d/o9nW2yh5\nbxoMiyMx0bmPS2Ki87Jv3drzvr5GE3z1R/UVRERqhoIF8ct7RcKJE+XbeB/Ly8/jlnm3MH/LAnj/\nSY6ueYA1GFJSnBe/rykE7yCid2/fQYD2hRARiQwFC+KX90u8SRPPEQDXMZcdh3YwPG04Ow7toO2y\nRexf8/NT5/btc0YRfAUB/oKIivqj5EYRkZqhnAXxKz0dBg1y9mgYNAjWrIFLLnHyFQCaNXOOAaz8\nfiV9pvfhWOEx1k5aS9eSn3vcK9CL3RVEuHIT/OUhePdHyY0iIjVDIwvil6+RgKys8u1e/uxl7lx0\nJ4M6DuKt69+i5Wktgx4tqGx/RESk+ilYkLBkZ0PK6GK+POs3HL7gKW49fwovpTxLw/oNAc8Xe3Z2\nxasYXHtBaKWDiEj00TSEVMh7CWVODoy48TCrzx7G4e7/hMXPsfWfL5wKFLy5VjG4tplOSQmvjYiI\nREZIwYIx5jtjTImPn2f9tB9kjFlpjPnBGHPcGLPZGHN/1XRdaor3i/zam78l65L+0H4tvPYuZN7D\n/n3G7/XBrGLQSgcRkegV6shCb6Ct28+VgAXm+ml/DHgWuBQ4D/gL8FdjzKSweisR4fHi7vwhG3r3\no0FDC9MzYPtQIHACo/c5X22DaSMiIpERUs6CtTbX/bMxJhnYZq1d4af9BmCD26H/GmNG4wQPL4XY\nV4mQU0sW+/wbrv0lzQ8OJeM3rzNh9ek+Exi98w9efBGmTAmc7FgdCZEiIlI1wk5wNMY0BMYBT4Zw\nzSXAAOCP4X6v1Lw33iyk98P3sb/DCyTuvJ9PH3+Cs9o28LsywbvS4pQpFa9i0EoHEZHoVZkEx1FA\nAjC7oobGmF3GmJNABvC8tXZmJb5XatDBEwe5dek15HZ6ia5fTafJJ09zw5gG5OT4v0b5ByIisaUy\nSycnAkustfuDaDsYaAb0B/5ujNlqrX2jooumTp1KQkKCx7HU1FRSU1PD6a+EaPOBzSSnJXM4/zA/\nWfcBXyweAlS8L4MqLYqI1Iy0tDTS0tI8jh0+fLjKv8dY723+grnImI7AdmCktXZRiNf+EbjZWts9\nQJskYP369etJSkoKuX9SsYrqGiz5dglj3x5Lx4SOLBi7gKG9O3sEAF26OBUXfcnJ0e6QIiKRkpWV\nRa9evQB6WWt9lNILXbgjCxOBbOB/YVxbH4gL83ulivjbwdFayzNrn+FXS3/Fdedcx5yUOcTHxYc0\nWqD8AxGR2BJyzoIxxgC3AbOstSVe5x43xsx2+3y3MWaYMaZb6c/twIPAq5Xst1SSr7yC/KJ8Ji2Y\nxAPvP8BvBv6GeTfOIz4uHijbl6FTJ4iPh927ywo0iYhIbAtnZGEo0AHwlaSYWHrOpR7wN+BsoAjY\nBvzaWvufML5XqpD3SEHLjjkMfXU0GXsyeGXkK9zS8xaP9q7RgsGDYedOZ/fJnTsD5y6IiEhsCDlY\nsNYuxZlK8HVugtfn54DnwuuaVCf3ugbNz9nE/quGU5B7kk9u+4T+7fv7vc7XiIT2dRARiW3aG6KO\nco0UPLV4Pt8OGUirpi3InJwZMFCA8rkKrVrBOedoXwcRkVimYKGOstbytxV/Y9Qbo7j2nGtZMWEF\nHRI6VHjdtGlOzkKDBs6fBQXOlIQ71VUQEYkt2qK6DjpZdJJJCyYx5/M5/OmyP/HwZQ9TzwQXN95+\ne1lwkJcHmzeXb6O6CiIisUXBQh2zL28fI98YyefZnzN3zFyuv+D6oK5z5SVkZHgeLyjw/Bwfr30d\nRERijYKFOmT93vWMeH0EACsnriQpMfiCV+51Gdw1agS9eyu5UUQklilYqCPmfjmX2965jR5n9mDe\njfNIjC8/VxBoVYO/PIQePbR0UkQk1inBMcaV2BIe+egRbnzrRlK6p/DxbR/7DRS8VzV07FhWeMk7\nDyEuzinStCikYt8iIlIbaWQhhm3ffYy+j48nt006nbb+jf+787c0bmB8th09uvyqhvz8sqWQ7nUZ\nNN0gIlK3KFiIUbsO76LH08M5dsa38Po8dn49gtE7/E8ZBFruuG+f9nsQEanLNA1Ry2VnO1MFXbuW\nTRms2bWGPtP7cNL8CDNWw9dOUuPu3f7vE2i5o5ZCiojUbRpZqOW8d48cfM+r7Ow5ib7t+nL0iXSO\n5bQ+1fbgwbLrvJMZX3wRpkxxPrdsCcbADz+UTTmIiEjdpWChljs1fWCK4Wd/5NsL/87EiybywrAX\nOPfPjTjm1rZFi7J/9g4ypkzRNIOIiPimaYhayH3q4cABoFEejB0FA5+g89dP89Lwl2hUvxHt23te\n5/7Z14ZQIiIivmhkoZZwnzY4cMBt5cLp31Fv8nBswvd037CYNx67hksvddolJEDTps6qhiZNnKkG\nF+8tqpWXICIi/ihYqCV8VlDstBxuTKFe0Rlsun8t3Vt3Z/Bg35UW8/I8pxp8LYUUERHxRcFCLVFu\nmiDpJbjuLth5GUk759K9dQvf7fzcQ0shRUQkWMpZqCVOTRPUK4Jr7ofhk2n+7RQGbF/CwrktyrcL\ndA8REZEQaGShlkhPh+E3HGLTT27kRNsP+ful/+Y3j9zls51reqFVK7AWcnM11SAiIuFTsFBLHKr/\nDYfGDKfxsRwW3fA+V3S+wmc7TS+IiEhV0zRELbB021L6vdSPeqYeGZMz/AYKIiIi1UHBQhSz1vJc\nxnNcO+daBrQfwJrb19CtRbdId0tEROoYBQtRqrC4kLsW38W9S+7l/v73szB1IQmNEyLdLRERqYOU\nsxSzuxIAAApzSURBVBCFco/nMubNMaz6fhUvD3+ZCZdMiHSXRESkDlOwEGW+zPmS4a8PJy8/j2Xj\nlzG44+BId0lEROo4TUNEkcXfLGbAjAE0bdiUjMkZChRERCQqKFiIAtZanlz9JMlpyVzR+QpW376a\ns08/O9LdEhERATQNEXH5RflMWTSF2Rtn84fBf+AvV/yFekYxnIiIRA8FCxGUfTSblLkprN+7njkp\nc7jpopsi3SUREZFyFCxEyMb9Gxn++nAKiwtZPmE5fdv1jXSXREREfNJ4dwTM2zyPgS8PpPVprcmc\nnKlAQUREopqChRpkreWx5Y+RMjeF5HOTWT5hOe2at4t0t0RERALSNEQNOVF4gokLJvL6F6/z6OWP\n8tCQhzDGRLpbIiIiFdLIQg14/uXnGTJrCAu+XsBb17/Fw5c9XKsDhbS0tEh3oUrF0vPE0rOAniea\nxdKzQOw9T1ULKVgwxnxnjCnx8fOsn/ajjDHvG2NyjDGHjTGrjTFXVU3Xa4dP937Kg//3INlHs1k5\nYSWjzx8d6S5VWqz9RxVLzxNLzwJ6nmgWS88Csfc8VS3UkYXeQFu3nysBC8z1034I8D5wLZAEfAQs\nNMb0DKu3tZC1lvhG8WRMzuCSxEsi3R0REZGQhZSzYK3Ndf9sjEkGtllrV/hpP9Xr0B+NMSOAZGBj\nKN9dW/Vp14cBHQbQtlnbSHdFREQkLGHnLBhjGgLjgBkhXGOAeOBguN8rIiIiNasyqyFGAQnA7BCu\n+TXQFP/TFi6NATZv3hxez6LM4cOHycrKinQ3qoyeJ3rF0rOAnieaxdKzQGw9j9u7s3FV3dNYa8O7\n0Jh3gXxr7Ygg298EvAgMt9Z+FETbOWF1TERERADGWWv/WxU3CmtkwRjTERgKjAyy/VjgP8CYigKF\nUu/hTHHsAE6G00cREZE6qjFwNs67tEqENbJgjPkTMBnoYK0tqaBtKvAScKO1dlE4nRQREZHICTnB\nsTRJ8TZglnegYIx53Bgz2+3zTTg5DQ8CmcaYM0t/mleu2yIiIlJTwlkNMRToAMz0cS6x9JzLZKA+\n8Dyw1+3nmTC+V0RERCIg7ARHERERqRu0N4SIiIgEpGBBREREAqrxYCGMzajaGmPmGGO+NsYUG2Oe\nquk++xNrG2uF8TyDjDErjTE/GGOOG2M2G2Pur+l++xPq83hdO8gYU2iMiYoqLWH83Vzmo22xMaZN\nTffdl3D+bowxjYwxjxljdhhjThpjthtjbqvBbvsVxt/PTLe/E/f2n9d03330LZy/m3HGmA3GmGPG\nmL3GmBnGmBY12W9/wnyee4wxX7n9f+2WmuxzIMaYesaYv5T++3/cGLPVGPNQENddboxZX/rfzjfG\nmPEhfbG1tkZ/gJZAG7efnwHFwKV+2ncCngZuBtYDT9V0n6vwWZ4GfgX0AroCjwH5QM9IP0uYz3Mx\ncCPQHegI3AQcBSZF+lnCeR636xKArcASICvSzxHm381lpee7ul8X6eeozN8NMB9YDfy09N+3fsCA\nSD9LmH8/8V7tzwJ+AB6uhc8yCCgC7in9//VA4HPgrUg/S5jPcxdwCBiDU6vgRuAIcF2kn6W0f38A\ncoBrSv87SCnt3y8CXHN26f+b/wH8pPTvqhC4MujvjYIHfwb4Jsi2HxFFwUJlnsXtmi+AhyLd9yp8\nnreB2ZHue2WeB0gD/gw8Ei3BQqjP4hYsNI90X6voea7B2VPm9Ej3tSqex0f7kaUv3A6R7nsYfzcP\nAt96HfsF8H2k+x7m86wC/u517ElgeaT7XtqXhcB0r2NvAa8EuObvwCavY2nA/4L93ojmLJgwNqOK\nVuE8izHRu7FWmM9zCTAA+LiauhW2YJ/HGDMB6IwTLESlEP5uDLChdFj4fWPMwOrvXeiCfJ5k4FPg\nt8aY3caZlnzCGFNlte+rSpj/X5sIfGCt3VU9vQrP/9/euYXWVURh+FttJZBKHkSJFSxiE+sN6kMJ\nQmsoeMFKEUXFIoqSFq2tIAr6EsSKGLQUoU/65K0vjYKg4qWiog9NRKX4VDX1ghdqi0LVxltKzviw\nJnA8rfvsOSc7s4X/g/2Q2Rn2+s+aM3udmTUzJbVMAmeb2fpYpx+4CXi9egvTKKmnhxN3Dv4LGDKz\nxVXZlsAEcLmZDQKY2Sp8dOeNgjqXAu+0lO3F++tSdHOQ1HzQyWFUdaXKg7VyUFqPmX0PnIHvqbE9\nhHCyPThy01ZP/PKNAWtDCA2P5WpJGd/8CNyFv2B78D1P3jezoRDCp9WbmEQZPecCl+Gd9nXA6cBT\nwGnApqoNTCSpLzCzZcB6YGOVRnVIWy0hhAkzuxUYj8HbEuBVfHShbpTxzV5gs5m9EkLYb2ar8TZ2\nCt7ujlRvZiGPA33A52Y2i+cejoYQ9hTUOZMT7T4C9JlZTwjh73YPzR0sjABvhhAOZ7ZjPkjSYr67\n5UP4wVo/V2pZZ6ToWQucikevT5jZlyGE8UqtS6dQj5ktwg8veziE8NVc8UIZl0hb34QQpoCppqIP\nzWwFcB+QlthUPWXa2iKgAdwSQpgGMLP7gZfMbGuZzm4BSe3X7gCO4jkZdaOtFjO7ENgFbAfexjfn\n24kfHLh5AWxMoYxvHgX6gcnYLxwGngMexNtgbm7G88M2Agfw3LFdZnYohLC7sqdmnHdZjs/RbUio\nU8uchVQt0cnTwNW5bZ8v3zTVHQU+y60hVQ/+a6MBzOCJP8fxOf+5snW5dcyDb3YA+3Jr6EQP3llP\ntZSdH320IreObvyDB3U7c9vehW9eAF5sKVsTvzv9uXV06ht8pPQs/EfDFuCX3BqiXd8Bd7eUjQIH\nCup80PruJAapZZ+bc2RhBB8GKZpn+b9QWov9+2Ctt6o2rEO68c1ifNi7TpTR8xtwcUvZNjzz/gb8\nBNQ60I1vLsGnJ+pEWT37gBvNrDeE8EcsW4m/kH6o0L5UkvxjZuvwFSt1zNsqq6UXD6ibaQCBeo3O\nJfkmhDCLH08wd3Lya9WZlkQvHiQ306B4K4RJfKqrmatieTkyRUaGd76PneTeGC3Z9MAqvKP7GNgd\n/74gd4SXqgUfOprBo9T+pqs2GeuJerYCG4CBeG0CfgUeya2j07bWcr9WqyESfXMvcC3+IroIzwA/\nTk1GSDrQsxT4FhjHl+oOA18AT+fW0U1bi/3ZRG7bu/TN7fgS8C14cvAa4KM66UrUM4gnQQ4AQ8Ae\n4CdgeW4d0b5n8dGFa/ClqtfjSynHCjSdAxzDV0WsjH33DHBF6edmEnslHhkN/McH8V5LWSP+f/P1\ndW6npWrBp1FadcwCz+TW0aGee/D11MfwOddPgDtza+imrbXcr1uwkOKbB4CDwO+xo3sXGM6toRvf\nAOfhyWfTeOCwA+jJraMLPX1Ry0hu2+dBy7bYF0zjIz3PA8ty6+hEDz69tT9qOQq8DAzm1tBk31Lg\nSeCb+P0+iK/eWtLGR8P4XkV/xjq3pTxXB0kJIYQQohCdDSGEEEKIQhQsCCGEEKIQBQtCCCGEKETB\nghBCCCEKUbAghBBCiEIULAghhBCiEAULQgghhChEwYIQQgghClGwIIQQQohCFCwIIYQQohAFC0II\nIYQo5B/2BkRbQQOSLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1198b1470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAFkCAYAAAAdXVDGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X+cHXV97/HXh02ykh9EIGwSARcjCBQMsEuCIchFUdEr\n4tW1rYvUFvvoBX/VGx9eebS2t157LdVaqb+4QalIH8J6FbdF0DYtWLBKCLILaDDBavghkiyBhCUY\nDCH53j++M2bO7MycOWfPnJlz9v18PM5j98z5zo8zO3vmc77fz/f7NeccIiIiImU4qOwDEBERkZlL\ngYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaB\niIiIiJSm8EDEzN5rZg+a2bNmdqeZrcgou8TMrjOzB8xsn5l9OqXcb5vZpmCb95nZG4p7ByIiIlKU\nQgMRM/td4G+BvwBOA+4D1pnZopRVeoHHgb8E7k3Z5pnA9cCXgFOBG4F/MrPfau3Ri4iISNGsyEnv\nzOxOYINz7gPBcwN+AXzWOffJOuv+O3CPc+6DseVfA+Y65y6ILFsflH1Pq9+DiIiIFKewGhEzmw0M\nAreGy5yPem4BVk1j06uCbUStm+Y2RUREpASzCtz2IqAHmIgtnwCOn8Z2l6Rsc0naCmZ2OHAe8BDw\n62nsW0REZKZ5AXAMsM4592SrN15kIFIl5wHXlX0QIiIiHewd+BzNlioyEHkC2Acsji1fDGybxna3\nNbHNhwC++tWvcuKJJ05j1zPPmjVruOKKK8o+jI6ic9YcnbfG6Zw1R+etMZs2beKiiy6C4F7aaoUF\nIs65vWY2BpwLfAt+k6x6LvDZaWx6fcI2XhssT/NrgBNPPJGBgYFp7HrmWbhwoc5Zg3TOmqPz1jid\ns+bovDWtkNSGoptmPg18JQhI7gLWAHOBrwCY2eXAi5xzvx+uYGanAAbMB44Inj/nnNsUFPkMcJuZ\nfRD4NjCMT4r9o4Lfi4iIiLRYoYGIc+7rwZghH8M3n9wLnOec2x4UWQIcHVvtHiDsUzwAXAg8DCwL\ntrnezC4EPh48/hN4s3PuJ0W+FxEREWm9wpNVnXNXAlemvHZxwrK6XYqdc98Evjn9oxMREZEyaa4Z\nyTQ8PFz2IXQcnbPm6Lw1TuesOTpv1VLoyKpVYWYDwNjY2JgSlERERBowPj7O4OAgwKBzbrzV21eN\niIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGI\niIiIlEaBiIiIiJRGgYiIiIiURoFIl5qYgLPOgpe+1P98/PGyj0hERGQqBSJdamgIfvAD2LLF/3zr\nW8s+IhERkakUiHSprVuzn4uIiFSBApEutXRp9nMREZEqmFX2AUgxRkd9c8zWrT4IGR0t+4hERESm\nUiDSpfr64PvfL/soREREsqlpRkREREqjQESky6krt4hUmQKRLqQbj0SpK7eIVJkCkS6kG085qhoA\nqiu3iFSZApEuVPaNpwo35DKOoaoBoLpyi0iVKRDpQmXfeKpwQy7jGMoOANOMjsLq1bBsmf+prtwi\nUiXqvtuFyh5DpAo35DKOYelSH/hEn1eBunKLSJUpEOlCZd94qnBDLuMYyg4ARUQ6kQIRabkq3JDL\nOIayA0ARkU5UeI6Imb3XzB40s2fN7E4zW1Gn/DlmNmZmvzazn5rZ78de/30z229m+4Kf+81sd7Hv\nopqqkBSaJLwh//zn/mdf38w8BhERqa/QQMTMfhf4W+AvgNOA+4B1ZrYopfwxwM3ArcApwGeAq83s\ntbGik8CSyKO/gMOvvCokhYqIiExH0TUia4CrnHP/4JzbDFwK7AbelVL+3cAW59yHnXMPOOe+ANwQ\nbCfKOee2O+ceDx7bC3sHFRZPwLz77urUioiIiORRWCBiZrOBQXztBuCjB+AWYFXKaq8IXo9al1B+\nvpk9ZGaPmNk/mdlvteiwO0o8AXPPHtWKiIhIZymyRmQR0ANMxJZP4JtTkixJKX+ImfUGzx/A16hc\nALwD/x7uMLMXteKgO8noKPT21i6rytgVIiIieXRcrxnn3J3AneFzM1sPbAIuweeipFqzZg0LFy6s\nWTY8PMzw8HABRzp9ExM+DyTa8yOadOkczJnja0JC0VqSeuuLiIhEjYyMMDIyUrNscnKy0H2aby0p\nYMO+aWY3MOSc+1Zk+VeAhc65tySsczsw5pz7YGTZHwBXOOcOzdjX14G9zrl3pLw+AIyNjY0xMDDQ\n5Dtqv7PO8kmoodWra7uHxl9fsAB+9rMDwUa99UVEROoZHx9ncHAQYNA5N97q7RfWNOOc2wuMAeeG\ny8zMgud3pKy2Plo+8LpgeSIzOwh4OdB1jRL1RgeNPz/iiNoaj+mMLlrVrsGdRudRRCRb0b1mPg38\nkZm908xOANYCc4GvAJjZ5WZ2baT8WmCZmX3CzI43s/cAbwu2Q7DOn5vZa83sJWZ2GnAd8GLg6oLf\nS9vVmzNmus+zXHBBbdfgN70p/7pygLpYi4hkKzRHxDn39WDMkI8Bi4F7gfMi3W2XAEdHyj9kZm8E\nrgD+GHgU+EPnXLQnzaHAF4N1d+JrXVYF3YO7Sr3RQaf7epb77st+LvlUYd4dEZEqKzxZ1Tl3JXBl\nymsXJyz7Hr7bb9r2Pgh8MO31blJvyPCs1zslUbVTjrNZVZh3R0Skygof4l2mp9kcg3iTwIteBPff\nn3+/y5dnP2+Vbm+6GB31ScLLlvmfmghPRKSWApGKa/ZGHW8C2LcPVqUNI5fg5ptrb6A335x/3ah6\ngVSzTRedkgSqOW9ERLJ13DgiM8nEhB+2PSrvjTreJADw7LP5992qmWTDQAr88bz1rbXbbbbpot52\nRUSkM6hGpMKGhmoHK4P8N+q1a6cuO/jg6R9To+rVeDTbdKEk0Po6pdZIRGY2BSIVFr+59vbmv1Ff\nemnt854eWJ86GkvrhTfBX/6ydnk8kGq26SJP1+SZfiPu9vwbEekOCkQqLH5zPf30/DfqeBDT3w8n\nndSa48ojvAmGNTq9vY0na2YFElk1KeF6/f3tvxFXKfhRrZGIdALliFTYdMYBKbvbaPymd+SRjedw\nZOWBZOWwRNfLOqYiVCl3pexrQEQkD9WIVNh0elyMjsLKlb4morcXnnuusW/nGzfCIYfA7Nn+ZyNd\nf2F6o7qGmv1Gn1auHTfiRx/Nft5O6josVaqhE0mjQKTCphMM9PX59fbs8Y8f/rCxpokzz4Rdu+D5\n5/3PRrr+wvQDIWg+mImXa6ZZqFk7dmQ/byd1He4uzQQVyhOSTqBApMKmGwxMJ0cg3tW3ka6/AM7B\npk3NB0LQ/Df6cL3+fj8j8ZIlje13Og47LPu5ePqm3rhmggrlCUknUCBSYc0GA3l7rGSJd/VttOvv\n0JAPnqIa/RBs9ht9uN5RR/ljePjh9n0bPOqo7OfiZd1UFaQkywoq0s5ZK5pIRYqmQKTCmg0G4j1W\nAE49tbGmiZtu8l1+wf9sZGTVjRvhjjumLi87YTbPB/d0KS8jn6y/jZoTkmUFFWnnTNejdAL1mqmw\nm26Cc8/1w7M3Egwk1Tz8/OeN5Qh85CN+v+B//umf5u/9ceaZvmkmav781n4I5pksL6vXyAUXwF13\n+d+3bIE3vQk2bJj+cbVqRNpul/W3UXNCsqxedGnnTNejdALViFRYPBg4//x839yTah4azfGIf7Dd\nfXf+2oPdu6cuiwcm0zExAccdl/2teWIC9u49kCy7YkXtB/d999WWjz/vZJ3QtJH1TT1+/W7fXs33\n0G5ZTZVqgpFOpkCkwuLBwK5d+aqpR0cPNKuEnGvswzz+QbZnz/Sqyn/1q9ZVsefJPxka8jUeYbLs\nnDkzp9dIvJr+2GOrdyPPuqmOjvok41De634mUxOMdDIFIhWW9K0mTzV1X5//hn9Q5K+7bx+cd17+\nfa9d628Gs2aBWWPHkPZtrFVV7EnbCfcZ1gbEm1ni6yxfnv28WeH+jznGd7k+5pj210o0G8BWRV8f\nHHFE7TI1z2RTV23pZApEKiz+zRDyV7medNLU5pC8zQ8TE7Vdh+PbqXcM/f3Jy1tVXRzfzoIFB74B\nXnCBrwV4/vnsdW6+ufYbZCPJuFHxZpBw/w8/3P4eO6FmA9gqUVNDd+uE5kNpI+dc1z+AAcCNjY25\nTvPjHzu3YIFzs2b5nxs35ltv2zbnfAhR+8hj9eqp65k519/vX5uYyF5/YsK5lSud6+31682Z49yK\nFfXXC4979Wrnli1L39fERHqZ3t6px57nmJsVP1dJ+wd/rO0yMeGvlfg5qLro337lSn/NZF0H0rni\n/zedcH3OROH/5JFHjjnAAQOugHu0es1U3LvedSAfYtcuuPjiA709sgwNJS8/66zkHiahtK63zvkx\nMfJk4EdHdQU/qmreHI08c7U00hOgtze9bJ6eN/XkrWlo5zf6vj742c+an6eoLPG//erVvqlBakWv\n28MP902nTzzR/DVcBvWM6gxp83a1mgKRivvRj7Kfp0n7xw6bCdJuzkldb+ttM0/ZZueJafQDavly\nP4praPZsX+2b9OHcignq4t1QZ83y52/vXv98zhx/TO0OBDqx26ZuTvnEr9tQ2ZMsNkITMnaGdv0P\nKkekS2X9Y2ddXFndfBv5sGjVPDGN7HNiwn87jCbXPvNMeq+RVtz4or0VFizwvYOee+5ApfNM67Ez\nHcoLySfrOu2U4E29fKotbXTuoigQqbCNGw98sw6dcsrUckmJX6OjfuyM3t6pvV6yPuDTRm81y/dh\nMTEBZ5zhayXM/E145crG54lp5gMq7LIbr9FJ6zWSduNrJJEu2lsh3tMjdPfdSsbLI9pTa8ECuOqq\nso+omrL+fxct6owkUPXy8aqWtBtOtLpkydTRuQtVROJJ1R50aLLq/Pm1CV0HHZSctJeW+JWUSFkv\n8W/jxqmJjuG6eSQlurYrEW3ZsuRE0bRk0bSk12YT6ZLeu5Lx8oufvwULlKSaJCkZOfzfXrFC113V\nRZOyq5ZUnvTZ7x9KVp2xnnmm9vn+/cnfHJKaGG67bWqS0ZFH1m8/PukkePppX4sRzbXIO85GUtVw\nO6qLJyayv00kfYtMy6NotslmdBRe/OLkbxHhyLSdlFDYbknjnxx7rE+81fk6IJ6MHCasbt06tSq9\nU5pqZpKsBNCy/16NjsDdKmqaqbB4k0r8eSipieE1r5labtGi/Pv+8pdrq8mvuSbfekk3/Ha09Q8N\nTQ3cwDcNNdrE02yuQl8fnH568mvTHZk2TdWqdqcj6Tx32mBs7RJt2pgzxzdJbtkyNQhWnk31ZAUb\nZf+9Gp1lvVUUiFRY3tl3k/IqwjlqolzO+V5uuw1e/vIDA5qF3YbzGB31tSlz5hxYdt99cP/9+dZv\nVto/t1njbdB581TCfJgXvMA/Vq70eQ2rV/tB3RYs8D97e/MdazO6aabapAH8oPxviVUXPz+9vUoC\nzaOsID5pQMaq/L3Wrz/wBXT+fDjtNH9sSbmJLVVEe0/VHnRojkh/f207XX9//nUPOihfnkTedfPm\niITtn2ZT2/uLlJWf0a7BzNLyGoocvCmeFxP9G+cZHK5qOnUwtqLU+xtu26bz1ayyBlXLGpCxCuKD\naN52m3OnnFJsjkg7goD3Ag8CzwJ3AivqlD8HGAN+DfwU+P2EMr8NbAq2eR/whjrb7MhApJlk09DJ\nJ0+9Sa5YkW+/STfzPIHIj3/sXE9P8vqzZuXbd7PC0VzbmSialhy7YkXtB83GjcV98GR9mHZq8mfV\nP6jbKf437Ok5MMJxeMOIv5539OWZLiuIn8mSrqmik1WLDkJ+Nwgo3gmcAFwF7AAWpZQ/BngG+CRw\nfBDE7AVeGylzZrDsg0GZjwF7gN/KOI6ODETCD+T4sOH1bqzbtiUPNb5yZb79Jt1cTzut/nrpGdfF\n14iE4r0GivyAyaqFadc3raybdlKgVMVvy51Yc9MuWT3B0v7fqvg3riINM59s1qyk66rYQKToHJE1\nwFXOuX9wzm0GLgV2A+9KKf9uYItz7sPOuQecc18Abgi2E/pj4J+dc58OyvwvYBx4X3FvoxxhQtqR\nR9Yur9dmPjSU3HNjfDxfW+js2VOXzcrRvyot43r+fN/2WLSJianjrkBxCWBhPkxaEnGo0RyHZscx\niefCNDP5XTvazeP7eP3ra/Ncli3r7KTbVsq6dtP+35RTky567e3d68daqkp+RlXEcxF7eorfZ2GB\niJnNBgaBW8NlzjkH3AKsSlntFcHrUeti5VflKNNVGu3FkfZB9Pzz+RIaTzpp6rI8N4a0ZNrdu+uv\nO10TE76r57331i7v6SluYKy+PtiwAbZt8wFJmkYDoVYloDYze3M7kl/j+4j/zX71q85Oum2ltARe\nSP9/K7vnRZVFr7277vJJ9TN9ULW4aMLqggXw3e8Wn6xaZI3IIqAHmIgtnwCWpKyzJKX8IWbWW6dM\n2jY7Vtgr4+67D4xSumJF/ci92UAlFO3xEtqxI3sdOHABx+3fD694RbHftNO67+7bBxdd1Np9xTkH\nmzYlvzZ/vh/yvZH33ao5V5zzwVlYYzNvXv2grB3zveTZpkaj9cIxQ+I9sVavTv5/W7BA3+yzaD6j\n+sKxpPbu9T/PPtsP51Akdd+tsHDI8j17/E0l7yy2o6NTu4xGJQUq0SrL+DdUgIUL6wcS4QWc5Jln\niv2mnfWBct99rd1X3NDQgRmS4/bu9QPDNfK+WzXnytAQ3HOPv3bA1zRcckl79t3IPpLs2aNaEfD/\nl+ef7wOzbdsO/C3BTykQn1bgiCP0zT6L5jOqpiJHVn0C2Acsji1fDGxLWWdbSvmnnXN76pRJ2+Zv\nrFmzhoULF9YsGx4eZnh4uN6qpWh2lNK+vqmz0IbSvjHVm+5550549FH/e9Ysnxs3Jq9vVvsh2qpv\nIhs3+hmD0wIBqN1vEaYzK3GS0dEDo2aGI7G26rjq7b9V+663j2OPrf2bzZ8/tUbr4Ydbv+9OMzRU\n+3+8Z48/Tw8/7AOU7dtry2/fnj7btLTn+u50IyMjjIyM1CybnJwsdqdFZMCGD3x33c9EnhvwC+B/\nppT/a+C+2LLrge9Enn8NuDFW5gfAlRnH0ZG9ZqYzb0taV9b4WCRhj4XkTOn0R1pPlKyeM0V0Jc2z\nv7Q5elolrffMggXlzv1R5rw/9ST19ol3/T7oIPWmyeo1k9QzDvwcVe0+V+r51N3Gxjq718yngT8y\ns3ea2QnAWmAu8BUAM7vczK6NlF8LLDOzT5jZ8Wb2HuBtwXZCnwFeb2YfDMp8FJ8U+/mC30vbhb0y\nenv9I09+SOixx5KXx3M9wpqQ559v7NjiVZph005WzUS0PbtVQ3fnmRth//7iqvnDnjpJTWGHHup/\nhn+/RmYhbuY4kmZgbvb6KUp4nKuC1PL1633NmnNTRwPev797Ro1tVlbTwXPPJS9/5pn2n6tuGuFX\n2q/QSe+cc183s0X4sT4WA/cC5znnwgrFJcDRkfIPmdkbgSvw3XQfBf7QOXdLpMx6M7sQ+Hjw+E/g\nzc65nxT5XsrQ1wff+pb/J9+6NTmJNG5iwpcPm1HiDjus9nm8qn7WLH/T+tWvapfPn+8zp9OqNOs1\n7YT7jgYqrWiecTmbXcLkx1ZXWYd5PEl27oRHHjnwfPbs4qrMo+c/2nS2YUMx+2tW/DjDSe2Ghuqv\nOxMTC0dHfRPMj37kn8+e7Se527EjO+hv97lSEqhMR+Gz7zrnrgSuTHltygwmzrnv4Ws4srb5TeCb\nLTnACgu7pIZt51u2wJvelH1zueCC9BsjwFFH1T5futRvN3TGGf7D7/Wv90mrzvkg5M47k7v1hvJ8\n8MRrY1qRKLZ0aXrQFRUmP9abfbhRSfN8HHnkgeNqdeCV9ziS9hUGqWEwuXYtXHppbXBZdG5B0gy7\nYZt93Jw5td/6Z2JiYV9f8v/zS1+aHYi0+1zFP0dm4t9Kmld4ICLNe/3rpybw/fCH2dPJ1+shMjpa\ne0M6/HBfff/EE7XbHB9v7FjjH0RJDjvMJ9G2MlGsvz9fIAKtDQTCcxifdn35cn8D3bq1mMArTdqN\nIFpD9uijB5o/tmypTfLNSkAu8jjB11YlOfVUXwOgxMKpsv7fDjqo/edKSaAyLUUknlTtQYcmqzYz\ndHhaAluY/Odc9tDGYdJZf79PtgzntaiXfBZNPmzkeKcrOkHTvHnOzZ3b+P6bSbRLmwconiS8YEF7\nEvjShnrPGoY+nhzajrk2suYjij7mzvWJvkp+9OLXaDh/0ezZU89d0fM6ycxTdLJq6UFCOx7dGoj0\n9k79oE6aayV8zJnjy2RN9pR242okiEj6cOztLeZmkhRUJQUjWTezZuacSDuHVZtIKyswjAcE7ehR\nk3d+nnYEsZ0k7RpN+n9v17xOMnN0eq8ZKdCePVOz1G++2Y+6mNSLIxxhM2tQn7RmjnrNLlHx3g+Q\nbyC2ZsSbW37xi6ll5s7Nzn/Im2gX7ZkSH79h0SL/Wvz8PfQQ3H9/6uG3VFLPmazmoKVL/bXSzrk2\n8jajxc305Me0a/TJJ2uXm7VnXieRVlIg0mHCG0c80Ag/mMJJ0JIsX+5/jo6m34DShnJv5EZwUMJV\nFe+t0yrxG+3OnVPnttm92yf55t3G4YcnjyIb7aK4a5fvjhyeQ+f8a/Eulfv3H+iqWrSkLpRZc5X0\n96dPmFeU+PWVdK3A1IkEZ3ryY9KXh4mJqQHxmWdmJ5WL147JHSU/JatW2OzZtbPJzp59IMg466za\n7rLRD6qJiak3RDNfWwLZwUq8i20zDj546jbivXVaJZ4kF++pEspK4o1vY+/e5K6w8WDsiCP8TRzg\n6KNJlWesk1aIH9+GDf7Y43/TWbMO9I5qt/ixvOhFPiC6++7aGaNPPdXXZCn50UtKBj3//Npzaeav\nXY2sWl9ad3cphwKRCpszpzYQiY4jkpWlPjTkv6FHzZ2bb59HHTX9obVvugle/WpfGwBw8snF3Uji\nQdVZZzV+/PFtvPSlta+HN/isLopZNUbOZfd0apX48YWzLcdrHc44o7wP3fj1FdbKPP64r7UKA8ZZ\ns9rTnbhTJH15CMcWCTnnu/rqplqfxj2pFjXNVFg0CIk/Dz+YkqrVk/Ik8k6tnhYwpE05Hjcx4W8o\nYRACPkfi/POLr/6cmJjaLBM68cT820nLoclq0oo3JUTt23egueT88+vvv9lq4/D44qJ/CzNfW1ZW\nVXT8HK5de2Ck1Z/8xNeK7Nnju6lnNadJOt1U69PkdxVTRAZs1R50aK+ZeFdcs3w9T9LmX4n34Ejr\ntpq0fnyOmjRZvSJWrGjo7Tcsa98rV+bfTlpX2Cx559jp7W38fTTaYySrC3e7eqHk7RKd9TfLe73P\nVGnzSc30HkZ5NPM/PpOp18wMdsoptc+d8yOt1vs2m5YYGo/60+aHWL8eenpqy+bN8cj6NhavSk7T\nbI1A1r6feCLfNiC7tinNTTdNPWfNmm61cZ7an6K/Neede6Rek5bmLEl3002+Vqm/3yck9/e3r/dT\np2vmf1yKo0Ckwm66aWrvmDyTxcWDBrPkCc/SbngnneSHd1+wwLfVL1gAV12V75hbUcXZ7ARa9bqq\nFukjH0nuthwXDy6TTLfaOD4nkVn7e6HkDabqHYeaGdKFN9MNG3yPuFYFwiLtpkCkwvr64PTTpy6v\n9+Ec77LpXPI4Hlk3vEsv9UHP88/7n6tW5auhCHMAkm4weW7C0HyNQDjbbNprRap3jL29/rzcdFP9\nbWXlouQRr/0JK+2j3Y2LPh95g6noe50/v/52ZCrNfCudToFIxSWNA1Hvw7mvz3ctjUq6UWbd8OIJ\nr7t25fugc0FvnW3bapf39OS7CUPzNQJ9ff7b4bx5tcvnzSu+6vXww5OXhwHII4/kqwKemJjaG6rR\nY087X2F343ZURecNpsJv9Xfc4ZuU5szxtTe9vT6oVDNDfeoBIp1O3XcrLJywLMz5OOww3+yS58M5\nz2yYWeOJ7NyZvu2sD7po//yofft8rUqem2uzE2iF5+vQQ/3YHfv3+5vay15W/NgK8aaP3l5fm9Vo\nINGK8Q3C8xcfm6OdtQtZ11aSoSHfUyZ0+unqgpqXZr6VTqdApMLiN/Xly/N/OK9d60dZfPZZ3/U2\nb45HKGtgs6wPuqwgZcuWfDfXRm9ioaQgyDm4557ix1ZISobduvXA6KZ5g5FWfLsNz9/jj3fOjKj6\nVt88zXwrnU6BSIXFm0camacjzPEA//Pkk+H22+Hss/OtHx94asECX7Vf74MuaXpys9oB1oq6yWRt\nt+gbW/x9h/MANVqr0cpvt80GdGXQt/rmddLfWSSJckQqLN48kjYPTJKkG++rX51//Xgb/89+li+/\nIL7exISvmYkq6iZTZq+Z6PtOmweo0e3MpK6Ya9c210tLamkOFelEqhGpsIULa5tHFi7Mv25SzUSe\n7qWhZr9lJa3Xrqrj6H4OP9zXxDzxRHuqq6PvO2seoEa2M5PEa/AuuWRmnofp0hwq0okUiFTY5GTt\n81/+0t/k8uQcrF0LL3957bKyxhlo1821Kjdxtdk3TjkiraHzKJ1ITTMVFq8BcS7/OAGXXjp12Xe/\n25rjkmwatbFxmvujNXQepRMpEKmwp55KXp7nW068zLJl+RNVQW3N0l4zNTem1XQepROpaabCnnsu\neXmebznxHJHt2xubil5tzdJOVWlWq7JwnJyswe50HqUTqUakwpKmlo9OnZ5VWxHthdDTk39k1JDa\nmkWqRUO5S7dSIFJhy5fXPl+xwn/bufTS+h9I0bli4r1l8gQVamsWqRZ9OZBupUCkwm6+uba99+ab\n/fI8H0hZH1J5goqi2pqVeyLSHH05kG6lHJEKS2vvzTMKZbxM3pFR6+17upR7ItIcdQuXbqVApAPl\n+UBKKlOFbqSqXhZpjhJRpVspEOlAeT6QqvqhpTlFREQkqrAcETM71MyuM7NJM9tpZleb2bwc633M\nzB4zs91m9m9mdmzs9dvMbH/ksc/MrizqfUhraZwDERGJKrJG5HpgMXAuMAf4CnAVcFHaCmZ2GfA+\n4J3AQ8BfQobXAAAXpklEQVT/AdaZ2YnOuXBUDQd8EfhzIOzgurv1hy9FqGpNjYiIlKOQQMTMTgDO\nAwadc/cEy94PfNvMPuSc25ay6geAv3TO3Rys805gAvhvwNcj5XY757YXcewiIiLSPkU1zawCdoZB\nSOAWfG3GGUkrmNlLgCXAreEy59zTwIZge1HvMLPtZvZjM/srMzu4pUdfEerqKiIi3a6oppklQM1t\n0zm3z8x2BK+lrePwNSBRE7F1rgMeBh4DlgOfBF4GvG36h10t6uoqIiLdrqFAxMwuBy7LKOKAE6d1\nRHU4566OPL3fzLYCt5rZS5xzDxa573ZTV1cREel2jdaIfAq4pk6ZLcA2oGbUCjPrAQ4LXkuyDZ98\nupjaWpHFwD2Ja3h3BesdC2QGImvWrGHhwoU1y4aHhxkeHs5arTTq6ioiIu00MjLCyMhIzbLJyclC\n92nOudZv1Cer3g+cHklWfR3wHeCotGRVM3sM+Bvn3BXB80PwQck7nXPfSFlnNfA94BTn3MaUMgPA\n2NjYGAMDA9N7c230+OPVHJRMRERmjvHxcQYHB8F3QBlv9fYLyRFxzm02s3XAl8zs3fjuu58DRqJB\niJltBi5zzt0YLPo74M/M7Gf47rt/CTwK3BiUXwZciA9ongROAT4N3J4WhHQydXUVEZFuV+Q4IhcC\nn8f3ltkP3IDvnht1HPCbthLn3CfNbC5+vJEXAv8BvCEyhshzwGuC7cwDfgF8A/h4cW9DREREilJY\nIOKce4qMwcuCMj0Jyz4KfDSl/KPAOdM/OhEREamCwoZ4FxEREalHgYiIiIiURoGIiIiIlEaBiIiI\niJRGgYiIiIiURoGIiIiIlEaBSMVpBl4REelmCkQqLpyBd8sW//Otby37iERERFpHgUjFaQZeERHp\nZgpEKi4+465m4BURkW5S5Fwz0gKjo1Nn4BUREekWCkQqTjPwiohIN1PTjIiIiJRGgYiIiIiURoGI\niIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgUiFacI7ERHpdgpEKkwT3omISLdTIFJh\nmvBORES6nQKRCtOEdyIi0u0010yFacI7ERHpdgpEKqyTJ7ybmPA5LtEgqq+v7KMSEZGqUdOMFEKJ\ntiIikocCESmEEm1FRCQPBSJSCCXaiohIHoUFImZ2qJldZ2aTZrbTzK42s3l11nmLma0zsyfMbL+Z\nLU8o02tmXwjK7DKzG8xM2QcVMzoKq1fDsmX+pxJtRUQkSZHJqtcDi4FzgTnAV4CrgIsy1pkH/Afw\n/4AvpZT5O+ANwBDwNPAF4JvAK1tx0NIanZxoKyIi7VNIIGJmJwDnAYPOuXuCZe8Hvm1mH3LObUta\nzzn31aBsP2AJ2z0EeBfwdufc7cGyi4FNZrbSOXdXEe+nLOp5IiIi3a6opplVwM4wCAncAjjgjGls\ndxAfPN0aLnDOPQA8Euyzq6jniYiIdLuiApElQM0Ubc65fcCO4LXpbPc559zTseUT09xuJanniYiI\ndLuGmmbM7HLgsowiDjhxWkdUoDVr1rBw4cKaZcPDwwwPD5d0RNmWLvW1IdHnIiIiRRkZGWFkZKRm\n2eTkZKH7bDRH5FPANXXKbAG2ATXZDGbWAxwWvNasbcAcMzskViuyOM92r7jiCgYGBqax+/bSEO8i\nItJOSV/Ox8fHGRwcLGyfDQUizrkngSfrlTOz9cALzey0SJ7IufgE1A15d5ewbAx4PtjWPwb7Oh54\nMbA+53Y7hnqeiIhItyuk14xzbrOZrQO+ZGbvxnff/RwwEu0xY2abgcucczcGzw/FBxVH4oOWE8zM\ngG3OuQnn3NNm9vfAp81sJ7AL+Czwg27rMSMiIjITFDmy6oXAZnxvmZuB7wGXxMocB0STNi4A7gFu\nwteIjADjsfXWBNu7AbgNeAw/poiIiIh0mMIGNHPOPUX24GU453piz68Frq2zzh7g/cFDREREOpjm\nmhEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQ\nERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBAR\nERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBER\nEZHSKBARERGR0igQERERkdIoEBEREZHSFBaImNmhZnadmU2a2U4zu9rM5tVZ5y1mts7MnjCz/Wa2\nPKHMbcFr4WOfmV1Z1PsQERGR4hRZI3I9cCJwLvBG4GzgqjrrzAP+A/gw4FLKOOCLwGJgCbA0KC8i\nIiIdZlYRGzWzE4DzgEHn3D3BsvcD3zazDznntiWt55z7alC2H7CMXex2zm1v8WGLiIhImxVVI7IK\n2BkGIYFb8LUZZ7Rg++8ws+1m9mMz+yszO7gF2xQREZE2K6RGBN9k8nh0gXNun5ntCF6bjuuAh4HH\ngOXAJ4GXAW+b5nZFRESkzRoKRMzscuCyjCIOnxdSGOfc1ZGn95vZVuBWM3uJc+7BrHXXrFnDwoUL\na5YNDw8zPDxcwJGKiIh0lpGREUZGRmqWTU5OFrpPcy4tJzShsNnhwOF1im0Bfg/4lHPuN2XNrAf4\nNfA259yNdfbTDzwInOqc+1GdsnOBZ4DznHP/llJmABgbGxtjYGCgzuGLiIhIaHx8nMHBQfB5n+Ot\n3n5DNSLOuSeBJ+uVM7P1wAvN7LRInsi5+ATUDXl3l7PcaUHZrTnLi4iISEUUkqzqnNsMrAO+ZGYr\nzGw18DlgJNpjxsw2m9mbI88PNbNTgJPwQcsJZnaKmS0OXl9mZn9mZgNm1m9mFwDXArc75zYW8V5E\nRESkOEWOI3IhsBnfW+Zm4HvAJbEyxwHRpI0LgHuAm/C1HCPAeGS954DX4IOcTcDfAN8I1hMREZEO\nU1SvGZxzTwEX1SnTE3t+Lb6GI638o8A5rTg+ERERKZ/mmhEREZHSKBARERGR0igQERERkdIoEBER\nEZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERER\nkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR\n0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERERkdIUFoiY2aFm\ndp2ZTZrZTjO72szmZZSfZWafMLMfmdkzZvZLM7vWzJbGyvWa2RfM7Akz22VmN5hZX1HvQ0RERIpT\nZI3I9cCJwLnAG4Gzgasyys8FTgX+N3Aa8BbgeODGWLm/C7Y3FGzzRcA3W3ngIiIi0h6zitiomZ0A\nnAcMOufuCZa9H/i2mX3IObctvo5z7ulgneh23gdsMLOjnHOPmtkhwLuAtzvnbg/KXAxsMrOVzrm7\ning/IiIiUoyiakRWATvDICRwC+CAMxrYzguDdZ4Kng/ig6dbwwLOuQeAR4J9ioiISAcpKhBZAjwe\nXeCc2wfsCF6ry8x6gb8GrnfOPRPZ7nNB7UnURN7tioiISHU0FIiY2eVmtj/jsc/MXjbdgzKzWcA3\n8LUh75nu9kRERKSaGs0R+RRwTZ0yW4BtQE1PFjPrAQ4LXksVCUKOBl4dqQ0hWHeOmR0SqxVZXG+7\nAGvWrGHhwoU1y4aHhxkeHq63qoiISNcbGRlhZGSkZtnk5GSh+zTnXOs36pNV7wdOjySrvg74DnBU\nUrJqUCYMQpYBr3LO7Yi9fgiwHZ+s+o/BsuOBTcAr0pJVzWwAGBsbG2NgYKAVb1FERGRGGB8fZ3Bw\nEHwHlPFWb7+QHBHn3GZgHfAlM1thZquBzwEj0SDEzDab2ZuD32fhu+EOABcBs81scfCYHWz3aeDv\ngU+b2TlmNgh8GfiBesyIiIh0nkK67wYuBD6P7y2zH7gB+ECszHFA2FZyJHB+8Pu9wU/D54m8Cvhe\nsGwNsC/YXi/wL8B7W3/4IiIiUrTCAhHn3FP4mo2sMj2R3x8GejKKh+X2AO8PHiIiItLBNNeMiIiI\nlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiU\nRoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRG\ngYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaBiIiIiJRGgYiIiIiURoGIiIiIlEaB\niIiIiJRGgYiIiIiURoGIiIiIlKawQMTMDjWz68xs0sx2mtnVZjYvo/wsM/uEmf3IzJ4xs1+a2bVm\ntjRW7jYz2x957DOzK4t6HyIiIlKcImtErgdOBM4F3gicDVyVUX4ucCrwv4HTgLcAxwM3xso54IvA\nYmAJsBT4cCsPXERERNpjVhEbNbMTgPOAQefcPcGy9wPfNrMPOee2xddxzj0drBPdzvuADWZ2lHPu\n0chLu51z24s4dhEREWmfompEVgE7wyAkcAu+NuOMBrbzwmCdp2LL32Fm283sx2b2V2Z28PQOV0RE\nRMpQSI0Ivsnk8egC59w+M9sRvFaXmfUCfw1c75x7JvLSdcDDwGPAcuCTwMuAt7XguEVERKSNGgpE\nzOxy4LKMIg6fFzItZjYL+EawvffU7MC5qyNP7zezrcCtZvYS59yD0923iIiItE+jNSKfAq6pU2YL\nsA3oiy40sx7gsOC1VJEg5Gjg1bHakCR3AQYcC2QGImvWrGHhwoU1y4aHhxkeHq6zCxERke43MjLC\nyMhIzbLJyclC92nOudZv1Cer3g+cHklWfR3wHeCopGTVoEwYhCwDXuWc25FjX6uB7wGnOOc2ppQZ\nAMbGxsYYGBho5i2JiIjMSOPj4wwODoLvgDLe6u0XkqzqnNsMrAO+ZGYrgmDhc8BINAgxs81m9ubg\n91nAN4EB4CJgtpktDh6zgzLLzOzPzGzAzPrN7ALgWuD2tCBEREREqquoZFWAC4HP43vL7AduAD4Q\nK3McELaVHAmcH/x+b/DT8Hkir8LXejwHvCbYzjzgF/galI8X8g5ERESkUIUFIs65p/A1G1lleiK/\nPwz0ZBQnGEvknFYcn8hMMTEBQ0OwdSssXQqjo9DXV389EZF20FwzIl1uaAh+8APYssX/fOtbyz4i\nEZEDFIiIdLmtW7Ofi4iUSYGISJdbujT7uYhImYpMVhWRChgd9c0x0RwREZGqUCAi0uX6+uD73y/7\nKEREkqlpRkREREqjQERERERKo0BERERESqNAREREREqjQERERERKo0BERERESqNAREREREqjQERE\nRERKo0BERERESqNAREREREqjQERERERKo0BERERESqNAREREREqjQERERERKo0BERERESqNARERE\nREqjQERERERKo0BERERESqNAREREREqjQERERERKo0BERERESqNAREREREqjQEQyjYyMlH0IHUfn\nrDk6b43TOWuOzlu1FBaImNmhZnadmU2a2U4zu9rM5tVZ5y/MbJOZPWNmO8zs38xsZaxMr5l9wcye\nMLNdZnaDmfUV9T5mOv3DNk7nrDk6b43TOWuOzlu1FFkjcj1wInAu8EbgbOCqOus8ALwXOBlYDTwE\n/KuZHR4p83fB9oaCbb4I+GYrD1xERETaY1YRGzWzE4DzgEHn3D3BsvcD3zazDznntiWt55z7Wmw7\nHwT+EFgO/LuZHQK8C3i7c+72oMzFwCYzW+mcu6uI9yMiIiLFKKpGZBWwMwxCArcADjgjzwbMbDZw\nCfAUcF+weBAfPN0alnPOPQA8EuxTREREOkghNSLAEuDx6ALn3D4z2xG8lsrM3gh8DZgLPAa81jm3\nI7Ld55xzT8dWm6iz3RcAbNq0KfcbEG9ycpLx8fGyD6Oj6Jw1R+etcTpnzdF5a0zk3vmCQnbgnMv9\nAC4H9mc89gEvA/4E2JSw/gRwSZ19HAwsA1YCXwK2AIuC14aBZxPW2QBcnrHNC/G1MXrooYceeuih\nR3OPCxuJGfI+Gq0R+RRwTZ0yW4BtQE1PFjPrAQ4LXkvlnHs22MYW4C4z+yk+T+QTwbpzzOyQWK3I\n4jrbXQe8A5/8+us6xy8iIiIHvAA4Bn8vbbmGAhHn3JPAk/XKmdl64IVmdlokT+RcwPC1F404COgN\nfh8Dng+29Y/Bvo4HXgysr3Pc1ze4XxEREfHuKGrDFjRdtH7DZt/B14q8G5gDfBm4yzn3e5Eym4HL\nnHM3mtlc4CPAt4CtwCLgfcDb8b1vNgXrXAm8AbgY2AV8FtjvnHtlIW9EREREClNUsir4vIzP43vL\n7AduAD4QK3McsDD4fR9wAvBOfBDyJPBD4KwwCAmsCcregK8p+Rf82CMiIiLSYQqrERERERGpR3PN\niIiISGkUiIiIiEhpZkQgYmbvNbMHzexZM7vTzFaUfUxVEUw0uD/2+EmszMfM7DEz2x1MRHhsWcdb\nBjN7pZl9y8x+GZyfCxLKZJ6jmThZY73zZmbXJFx734mVmVHnzcz+xMzuMrOnzWzCzP7RzF6WUE7X\nWyDPOdO1NpWZXWpm9wUT006a2R1m9vpYmbZcZ10fiJjZ7wJ/C/wFcBp+uPh1Zrao1AOrlo34sViW\nBI+zwhfM7DJ876X/jh9k7lf48zenhOMsyzzgXuA9+EF9auQ8RzNxssbM8xb4Z2qvveHY6zPtvL0S\n+Bx+KozXALPxE38eHBbQ9TZF3XMW0LVW6xfAZcAAfvqU7wI3mtmJ0ObrrIhR0qr0AO4EPhN5bsCj\nwIfLPrYqPPAB2njG648BayLPDwGeBX6n7GMv6XztBy5o5BwFz/cAb4mUOT7Y1sqy31OJ5+0aYDRj\nHZ0334NwP773oK635s+ZrrV85+5J4OLg97ZdZ11dIxJMnDdI7SR5Dt+lWJPkHXBcUH3+czP7qpkd\nDWBmL8F/c4iev6fxg9Lp/JH7HJ2OJmtMc05Qnb7ZzK40s8Mir2mSS3ghvjZpB+h6y6nmnEXoWkth\nZgeZ2dvxc7zd0e7rrMhxRKpgEdCDn+MmagIfuYmvMfoD4AFgKfBR4HtmdjL+QnQkn7/MyQtnkDzn\naDHNTdbY7f4ZX437IPBS/FxW3zGzVcEXhmYnuewKZmb4qu/vO+fCvC1dbxlSzhnoWksUfM6vxw/h\nvgtfu/GAma2ijddZtwciUodzLjp3wEYzuwt4GPgdYHM5RyUzgXPu65Gn95vZj4GfA+cA/17KQVXL\nlcBvAavLPpAOknjOdK2l2gycgh9Y9G3AP5jZ2e0+iK5umgGewI/Cuji2vN4keTOWc24S+ClwLP4c\nGTp/WfKco99M1phRZsZzzj2I/58NM/Nn7Hkzs88D/xU4xzm3NfKSrrcUGedsCl1rnnPueefcFufc\nPc65j+A7c3yANl9nXR2IOOf24ifKOzdcFlTdnUuBE/h0MjObj//nfCz4Z91G7fk7BJ+drvPHbz7Q\n6p2j6GSNYZm6kzXONGZ2FHA4fq4pmKHnLbihvhl4lXPukehrut6SZZ2zlPK61pIdBPS2/TorO0u3\nDVnAvwPsxs9hcwJwFT4z+Iiyj60KD+Bv8N2u+oEzgX/Dt/EdHrz+4eB8vQl4OfBPwH8Cc8o+9jae\no3n46stT8Rnh/yN4fnTec4SvMn4QXxU8CPwA+I+y31tZ5y147ZPBB1t/8GF2N7AJmD1Tz1vwfnfi\nu6QujjxeECmj662Bc6ZrLfW8/VVwzvqBk/F5M88Dr273dVb6yWjTCX8P8BC+69F64PSyj6kqD2AE\n3535WXy28/XAS2JlPorvyrUbWAccW/Zxt/kc/ZfgRrov9vhy3nOEn6Dxc/jq4F3AN4C+st9bWecN\nnxz3L/hvXb8GtgD/l9gXhJl23lLO1z7gnbFyut5ynjNda6nn7ergXDwbnJt/JQhC2n2dadI7ERER\nKU1X54iIiIhItSkQERERkdIoEBEREZHSKBARERGR0igQERERkdIoEBEREZHSKBARERGR0igQERER\nkdIoEBEREZHSKBARERGR0igQERERkdL8fzdWaam7ga2hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119d1ab00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_1 = coeff[0] + coeff[1]* pres + coeff[2]*temp + coeff[3]*sal + coeff[4]*oxy\n",
    "#OR\n",
    "y_1 = np.dot(A,coeff)\n",
    "\n",
    "residual = y_1 - y\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(y,y_1)\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(y,y_1,'.')\n",
    "x_line = np.linspace(7.2,7.9,299)\n",
    "y_line = slope*x_line + intercept\n",
    "ax.plot(x_line,y_line)\n",
    "\n",
    "# Plot the residual\n",
    "fig,ax = plt.subplots()\n",
    "ax.plot(subdf['CTDPRS'],residual,'.')\n",
    "\n",
    "print('R-Squared:',r_value**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `statsmodels` package can provide statistical analysis on the significance of the least squares fit, including the F statistic and adjusted R-squared.\n",
    "\n",
    "See http://statsmodels.sourceforge.net/devel/examples/notebooks/generated/ols.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 PH_TOT   R-squared:                       0.966\n",
      "Model:                            OLS   Adj. R-squared:                  0.966\n",
      "Method:                 Least Squares   F-statistic:                     2992.\n",
      "Date:                Sat, 11 Mar 2017   Prob (F-statistic):          1.75e-305\n",
      "Time:                        09:33:41   Log-Likelihood:                 962.59\n",
      "No. Observations:                 422   AIC:                            -1915.\n",
      "Df Residuals:                     417   BIC:                            -1895.\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
      "------------------------------------------------------------------------------\n",
      "const          7.0659      0.142     49.878      0.000         6.787     7.344\n",
      "x1             0.0004   2.78e-05     13.148      0.000         0.000     0.000\n",
      "x2             0.0252      0.001     20.759      0.000         0.023     0.028\n",
      "x3            -0.0030      0.004     -0.697      0.486        -0.011     0.005\n",
      "x4             0.0021   3.14e-05     65.418      0.000         0.002     0.002\n",
      "==============================================================================\n",
      "Omnibus:                      314.923   Durbin-Watson:                   1.454\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):            14169.009\n",
      "Skew:                           2.642   Prob(JB):                         0.00\n",
      "Kurtosis:                      30.891   Cond. No.                     2.08e+04\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 2.08e+04. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n"
     ]
    }
   ],
   "source": [
    "res = sm.OLS(y, A).fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes\n",
    "* The R-squared and adjusted R-squared are high, but it would also be useful to compute the variance inflation factor to check for multicollinearity.\n",
    "* The warning about the high condition number occurs because we used variables with very different magnitudes. This can lead to problems like round-off errors when subtracting small number from big number. Standardizing the input variables and transforming them to z-scores would help solve this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}