{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8386fb13",
   "metadata": {},
   "source": [
    "# Lineær Regresjon (løsningsforslag til oppgave fra 22. januar)\n",
    "\n",
    "La oss fortsette med det klassiske datasettet som inneholder data om størrlsen (lende og bredde) av begerbladene (engelsk: sepal) (de ytre bladene i en blomst) og kronbladene (engelsk: petal) til tre ulike type Iris (setosa, versicolor og virginica) - på norsk hhv.: vill iris, praktiris og blått flagg iris). Vi skal bruke dette datasettet til å se nærmere på lineær regresjon. Først må vi importere noen pakker og lese inn datasettet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "40a1502f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9bc47618",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression, LinearRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.special import expit\n",
    "from scipy.stats import expon\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be254927",
   "metadata": {},
   "source": [
    "Datasettet er tilgjengelig direkte fra ScikitLearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e8da3bad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal length (cm)</th>\n",
       "      <th>sepal width (cm)</th>\n",
       "      <th>petal length (cm)</th>\n",
       "      <th>petal width (cm)</th>\n",
       "      <th>target</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)  \\\n",
       "0                5.1               3.5                1.4               0.2   \n",
       "1                4.9               3.0                1.4               0.2   \n",
       "2                4.7               3.2                1.3               0.2   \n",
       "3                4.6               3.1                1.5               0.2   \n",
       "4                5.0               3.6                1.4               0.2   \n",
       "\n",
       "   target  \n",
       "0     0.0  \n",
       "1     0.0  \n",
       "2     0.0  \n",
       "3     0.0  \n",
       "4     0.0  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()\n",
    "df = pd.DataFrame(data= np.c_[iris['data'], iris['target']],\n",
    "                     columns= iris['feature_names'] + ['target'])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67c46057",
   "metadata": {},
   "source": [
    "Nå er typen iris lagret som 0, 1 og 2. Vi vil gjerne oversette dette til faktiske navn slik at vi ikke glemmer hva de ulike tallene betyr. Dette kan vi gjøre på flere måter, men en mulighet er å bruke såkalte lambda-funksjoner som utfører en test av hvert enkelt element in en kolonne. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "89ec4edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['species'] = df['target'].apply(lambda x: \"setosa\" if x == 0.0 else (\"versicolor\" if x == 1.0 else \"virginica\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6960e646",
   "metadata": {},
   "source": [
    "Det er keitete å ha kolonnenavn som inneholder mellomrom så la oss lage nye enklere navn. Når vi setter *inplace=True* betyr det at data-framen df automatisk blir oppdatert med de nye kolonnenavnene. Hvis ikke må man bruke *df = df.rename()* for at endringene skal bli oppdatert."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c2a88bd0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal_length</th>\n",
       "      <th>sepal_width</th>\n",
       "      <th>petal_length</th>\n",
       "      <th>petal_width</th>\n",
       "      <th>target</th>\n",
       "      <th>species</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal_length  sepal_width  petal_length  petal_width  target species\n",
       "0           5.1          3.5           1.4          0.2     0.0  setosa\n",
       "1           4.9          3.0           1.4          0.2     0.0  setosa\n",
       "2           4.7          3.2           1.3          0.2     0.0  setosa\n",
       "3           4.6          3.1           1.5          0.2     0.0  setosa\n",
       "4           5.0          3.6           1.4          0.2     0.0  setosa"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.rename({'sepal length (cm)': 'sepal_length'}, axis='columns',inplace=True)\n",
    "df.rename({'sepal width (cm)': 'sepal_width'}, axis='columns',inplace=True)\n",
    "df.rename({'petal length (cm)': 'petal_length'}, axis='columns',inplace=True)\n",
    "df.rename({'petal width (cm)': 'petal_width'}, axis='columns',inplace=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91a890b3",
   "metadata": {},
   "source": [
    "Tilslutt, la oss droppe kolonnen med *target* siden vi når har laget en ny kolonne *species* som inneholder den samme informasjonen. Det gjør vi enkelt slik (merk at *axis=1\" betyr at vi vil droppe en kolonne. Dersom vi ikke spesifiserer noe er *axis* satt til 0, dvs. at en hel rad blir fjernet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "520dc0b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal_length</th>\n",
       "      <th>sepal_width</th>\n",
       "      <th>petal_length</th>\n",
       "      <th>petal_width</th>\n",
       "      <th>species</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal_length  sepal_width  petal_length  petal_width species\n",
       "0           5.1          3.5           1.4          0.2  setosa\n",
       "1           4.9          3.0           1.4          0.2  setosa\n",
       "2           4.7          3.2           1.3          0.2  setosa\n",
       "3           4.6          3.1           1.5          0.2  setosa\n",
       "4           5.0          3.6           1.4          0.2  setosa"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.drop(['target'],axis=1,inplace=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff00d0b7",
   "metadata": {},
   "source": [
    "## Forberedelser til lineær regresjon\n",
    "\n",
    "Vi skal starte med å se om det er en sammenheng mellom lengde og bredde til begerbladene til Iris Setosa. Det første vi må gjøre er å trekke ut disse to kolonnene fra data-framen og konvertere de til **numpy-vektorer**. Det gjør vi ved hjelp av funksjonen *to_numpy()*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "87aaaa80",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_setosa_sepal_length = df[df[\"species\"] == \"setosa\"].sepal_length.to_numpy()\n",
    "X_setosa_sepal_width = df[df[\"species\"] == \"setosa\"].sepal_width.to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "324ffacc",
   "metadata": {},
   "source": [
    "Naivt - uten at jeg på noen måte er biolog - vil jeg anta at det er en korrelasjon mellom lengde og bredde av begerbladene. La oss ta en titt på data-settet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ca3f849a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCCUlEQVR4nO3dfVhUdf7/8dcAKVo4hisypilpmWCaZq7UarVm3hSbW32tLbNbSzM1zV2ztoirkmqrtbaydFMzKr+baEkZZuV9mhlYImZu4k0uxBYJqIHBnN8f/ODrOKBzcGYOc+b5uK65Luczn5nz/swcZ16cm89xGIZhCAAAwCYirC4AAADAnwg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVppMuElPT5fD4dB9993XYJ9Vq1bJ4XB43b755pvgFQoAAJq0KKsLkKQvvvhCs2fPVs+ePX3qv2PHDrVq1aruftu2bQNVGgAACDGWh5uDBw/qpptu0pw5c/T444/79Jy4uDi1bt26Uctzu936z3/+o5iYGDkcjka9BgAACC7DMFReXq727dsrIuL4O54sDzfjx4/XlVdeqcsvv9zncNO7d29VVFQoMTFRf/3rX3XZZZc12LeyslKVlZV19/fv36/ExMSTrhsAAATfvn371KFDh+P2sTTcLFy4UDk5Ofriiy986u9yuTR79mxdcMEFqqys1BtvvKFBgwZp1apVGjhwYL3PSU9PV1pamlf7vn37PHZtAQCApqusrEwdO3ZUTEzMCfs6DMMwglCTl3379qlv37766KOP1KtXL0nSpZdeqvPPP18zZ870+XVSUlLkcDi0dOnSeh8/dstN7ZtTWlpKuAEAIESUlZXJ6XT69Ptt2dlSX375pYqLi3XBBRcoKipKUVFRWr16tV544QVFRUWpurrap9fp37+/du7c2eDjzZs3V6tWrTxuAADAvizbLTVo0CBt3brVo+22227Tueeeq2nTpikyMtKn18nNzZXL5QpEiQAAIARZFm5iYmLUo0cPj7ZTTz1Vbdq0qWufPn269u/frwULFkiSZs6cqc6dOyspKUlHjhxRRkaGMjMzlZmZGfT6AQBA02T52VLHU1hYqL1799bdP3LkiKZOnar9+/erRYsWSkpK0gcffKDhw4dbWCUAAGhKLDug2CpmDkgCAABNQ0gcUAwAABAIhBsAAGArhBsAAGArhBsAAGArTfpsKQAAYF6129CmghIVl1coLiZa/RJiFRkRPheLJtwAAGAj2XmFSsvKV2FpRV2byxmt1JREDe0RHpPeslsKAACbyM4r1LiMHI9gI0lFpRUal5Gj7LxCiyoLLsINAAA2UO02lJaVr/omr6ttS8vKV7Xb/tPbEW4AALCBTQUlXltsjmZIKiyt0KaCkuAVZRHCDQAANlBc3nCwaUy/UEa4AQDABuJiov3aL5QRbgAAsIF+CbFyOaPV0AnfDtWcNdUvITaYZVmCcAMAgA1ERjiUmpIoSV4Bp/Z+akpiWMx3Q7gBAMAmhvZwadaoPop3eu56indGa9aoPmEzzw2T+AEAYCNDe7g0ODGeGYoBAIB9REY4lNyljdVlWIbdUgAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFYINwAAwFairC4AAIBAq3Yb2lRQouLyCsXFRKtfQqwiIxxWl4UAIdwAAGwtO69QaVn5KiytqGtzOaOVmpKooT1cFlaGQGG3FADAtrLzCjUuI8cj2EhSUWmFxmXkKDuv0KLKEEiEGwCALVW7DaVl5cuo57HatrSsfFW76+uBUEa4AQDY0qaCEq8tNkczJBWWVmhTQUnwikJQEG4AALZUXN5wsGlMP4QOwg0AwJbiYqL92g+hg3ADALClfgmxcjmj1dAJ3w7VnDXVLyE2mGUhCAg3AABbioxwKDUlUZK8Ak7t/dSUROa7sSHCDQDAtob2cGnWqD6Kd3rueop3RmvWqD7Mc2NTTOIHALC1oT1cGpwYzwzFYYRwAwCwvcgIh5K7tLG6DAQJu6UAAICtEG4AAICtEG4AAICtEG4AAICtcEAxAFik2m1wBg8QAIQbALBAdl6h0rLyPS7s6HJGKzUlkblXgJPEbikACLLsvEKNy8jxumJ1UWmFxmXkKDuv0KLKAHsg3ABAEFW7DaVl5cuo57HatrSsfFW76+sBwBeEGwAIok0FJV5bbI5mSCosrdCmgpLgFQXYDOEGAIKouLzhYNOYfgC8EW4AIIjiYqJP3MlEPwDeCDcAEET9EmLlckaroRO+Hao5a6pfQmwwywJshXADAEEUGeFQakqiJHkFnNr7qSmJzHcDnATCDQAE2dAeLs0a1UfxTs9dT/HOaM0a1Yd5boCTxCR+AGCBoT1cGpwYzwzFQAAQbgDAIpERDiV3aWN1GYDtsFsKAADYCuEGAADYCuEGAADYCuEGAADYCgcUAwAAv6h2G03iDMAms+UmPT1dDodD991333H7rV69WhdccIGio6N11lln6ZVXXglOgQAAoEHZeYX63VOf6k9zNmrSwi3605yN+t1Tnyo7rzDotTSJcPPFF19o9uzZ6tmz53H7FRQUaPjw4RowYIByc3P14IMPauLEicrMzAxSpQAA4FjZeYUal5HjdcX7otIKjcvICXrAsTzcHDx4UDfddJPmzJmj008//bh9X3nlFZ155pmaOXOmunfvrjvvvFO33367nnnmmSBVCwAAjlbtNpSWlS+jnsdq29Ky8lXtrq9HYFgebsaPH68rr7xSl19++Qn7btiwQVdccYVH25AhQ7R582b9+uuv9T6nsrJSZWVlHjcAAOAfmwpKvLbYHM2QVFhaoU0FJUGrydJws3DhQuXk5Cg9Pd2n/kVFRWrXrp1HW7t27VRVVaUff/yx3uekp6fL6XTW3Tp27HjSdQMAgBrF5Q0Hm8b08wfLws2+ffs0adIkZWRkKDo6+sRP+P8cDs+jrg3DqLe91vTp01VaWlp327dvX+OLBgAAHuJifPsN97WfP1h2KviXX36p4uJiXXDBBXVt1dXVWrNmjV588UVVVlYqMjLS4znx8fEqKiryaCsuLlZUVJTatKn/+izNmzdX8+bN/T8AAACgfgmxcjmjVVRaUe9xNw7VXPG+X0Js0GqybMvNoEGDtHXrVm3ZsqXu1rdvX910003asmWLV7CRpOTkZK1YscKj7aOPPlLfvn11yimnBKt0AADw/0VGOJSakiipJsgcrfZ+akpiUOe7sSzcxMTEqEePHh63U089VW3atFGPHj0k1exSGj16dN1zxo4dqz179mjKlCnavn275s6dq9dee01Tp061ahgAAIS9oT1cmjWqj+Kdnrue4p3RmjWqj4b2cAW1niY9Q3FhYaH27t1bdz8hIUHLli3T5MmT9dJLL6l9+/Z64YUXdO2111pYJQAAGNrDpcGJ8U1ihmKHUXtEbpgoKyuT0+lUaWmpWrVqZXU5AADAB2Z+vy2f5wYAAMCfCDcAAMBWCDcAAMBWCDcAAMBWmvTZUgAAezlS5dYbG3ZrT8lhdYptqZuTO6tZFH9nw78INwCAoEhflq85awt09MWhn1i2XWMGJGj68ETrCoPtEG4AAAGXvixfr64p8Gp3G6prJ+DAX9gWCAAIqCNVbs1Z6x1sjjZnbYGOVLmDVBHsjnADAAioNzbs9tgVVR+3UdMP8AfCDQAgoPaUHPZrP+BECDcAgIDqFNvSr/2AEyHcAAAC6ubkzjrRtRMjHDX9AH8g3AAAAqpZVITGDEg4bp8xAxKY7wZ+w6ngAICAqz3N+9h5biIcYp4b+J3DMIwTHMNuL2YumQ4A8C9mKEZjmfn9ZssNACBomkVF6I4BZ1ldBmyOuAwAAGyFcAMAAGyFcAMAAGyFcAMAAGyFA4oBhL1fjlRrxrJ87f7psDq3aakHhyeqRbNIq8sCGq3abWhTQYmKyysUFxOtfgmxijzRTIo2QrgBENbGLPhCK/KL6+6v3Sm9sXGvBifGac7oCy2sDGic7LxCpWXlq7C0oq7N5YxWakqihvZwWVhZ8LBbCkDYOjbYHG1FfrHGLPgiyBUBJyc7r1DjMnI8go0kFZVWaFxGjrLzCi2qLLgINwDC0i9HqhsMNrVW5BfrlyPVQaoIODnVbkNpWfmqb2be2ra0rHxVu+0/d69Pu6ViY2NNvajD4VBOTo46derUqKIAINBmLMv3ud9jI84LcDXAydtUUOK1xeZohqTC0gptKihRcpc2wSvMAj6FmwMHDmjmzJlyOp0n7GsYhu655x5VV/PXDoCma/dPh/3aD7BacXnDwaYx/UKZzwcU33DDDYqLi/Op74QJExpdEAAEQ+c2LbV2p2/9gFAQFxPt136hzKdjbtxut8/BRpLKy8t11llcOwRA0/Wgj1eh9rUfYLV+CbFyOaPV0AnfDtWcNdUvwdyhJqGIA4oBhKUWzSI1OPH4f7QNToxjvhuEjMgIh1JTasL4sQGn9n5qSmJYzHfjMAzD9GHT+/fv1/r161VcXCy32+3x2MSJE/1WXCCYuWQ6APtr6HRw5rlBqLLrPDdmfr9Nh5t58+Zp7Nixatasmdq0aSOH4/8SoMPh0K5duxpXdZAQbgAcixmKYTd2nKE4oOGmY8eOGjt2rKZPn66IiNDbq0W4AQAg9Jj5/TadTg4fPqwbbrghJIMNAACwP9MJ5Y477tA777wTiFoAAABOmundUtXV1brqqqv0yy+/6LzzztMpp5zi8fhzzz3n1wL9jd1SAACEHjO/36avCj5jxgwtX75c3bp1kySvA4oBAACsZDrcPPfcc5o7d65uvfXWAJQDAMFnxzNLTsSqMYfje43gMx1umjdvrosvvjgQtQBA0Nl1TpDjsWrM4fhewxqmDyieNGmS/vGPfwSiFgAIquy8Qo3LyPG6knJRaYXGZeQoO6/QosoCx6oxh+N7DeuY3nKzadMmffrpp3r//feVlJTkdUDx4sWL/VYcAARKtdtQWla+6jujwlDNdPVpWfkanBhvm90mVo05HN9rWMt0uGndurWuueaaQNQCAEGzqaDEayvC0QxJhaUV2lRQouQubYJXWABZNeZwfK9hLdPhZt68eYGoAwCCqri84R/bxvQLBVaNORzfa1jL9DE3BQUF2rlzp1f7zp07tXv3bn/UBAABFxcT7dd+ocCqMYfjew1rmQ43t956qz777DOv9s8//5zTwwGEjH4JsXI5o9XQER4O1ZzJ0y8hNphlBZRVYw7H9xrWMh1ucnNz6z0VvH///tqyZYs/agKAgIuMcCg1JVGSvH50a++npiTa6gBXq8Ycju81rGU63DgcDpWXl3u1l5aWqrq62i9FAUAwDO3h0qxRfRTv9NwdEu+M1qxRfWw594pVYw7H9xrWMX1tqauuukotW7bU22+/rcjISEk115u6/vrrdejQIX344YcBKdRfuLYUgGOF46y5zFCMUGPm99t0uMnPz9fAgQPVunVrDRgwQJK0du1alZWV6dNPP1WPHj0aX3kQEG4AAAg9Zn6/Te+WSkxM1Ndff62RI0equLhY5eXlGj16tL755psmH2wAAID9md5yE+rYcgMAQOjx+5abr7/+Wm632+cCtm3bpqqqKp/7AwAA+ItP4aZ379766aeffH7R5ORk7d27t9FFAQAANJZPl18wDEMPP/ywWrZs6dOLHjly5KSKAhrCmRbBwZk0AEKZT+Fm4MCB2rFjh88vmpycrBYtWjS6KKA+2XmFSsvK97gAn8sZrdSURObI8COr3mc+XwD+wgHFCAnZeYUal5GjY1fW2r/pmQTMP6x6n/l8AZxIQE8FB4Kt2m0oLSvf64dPUl1bWla+qt1hldP9zqr3mc8XgL8RbtDkbSoo8dhVcSxDUmFphTYVlASvKBuy6n3m8wXgb4QbNHnF5Q3/8DWmH+pn1fvM5wvA3wg3aPLiYqJP3MlEP9TPqveZzxeAvxFu0OT1S4iVyxmthk4IdqjmrJp+CbHBLMt2rHqf+XwB+Fujws0bb7yhiy++WO3bt9eePXskSTNnztR7773n1+IASYqMcCg1JVGSvH4Aa++npiQyH8pJsup95vMF4G+mw82sWbM0ZcoUDR8+XAcOHFB1dbUkqXXr1po5c6a/6wMkSUN7uDRrVB/FOz13TcQ7ozlN2I+sep/5fAH4k+l5bhITEzVjxgyNGDFCMTEx+uqrr3TWWWcpLy9Pl156qX788cdA1eoXzHMT2pjBNjiYoRhAU2Pm99unGYqPVlBQoN69e3u1N2/eXIcOHTL7coApkREOJXdpY3UZtmfV+8znC8AfTO+WSkhI0JYtW7zaP/zwQyUmJvqjJgAAgEYzveXmz3/+s8aPH6+KigoZhqFNmzbp7bffVnp6uv75z38GokYAAACfmQ43t912m6qqqvSXv/xFhw8f1o033qgzzjhDzz//vG644YZA1AgAAOCzk7pw5o8//ii32624uDh/1hRQHFAMAEDoCegBxUf7zW9+czJPBwAPR6rcemPDbu0pOaxOsS11c3JnNYsK/FyjVi2Xs8OCh/c6vPi05aZ3795yOHxbCXJycnxe+KxZszRr1izt3r1bkpSUlKRHHnlEw4YNq7f/qlWrdNlll3m1b9++Xeeee65Py2TLDdA0pS/L15y1BTr64t8RDmnMgARNHx64kxWsWm52XqHSsvI9LhrqckYrNSWReX38jPfaHvy+5WbEiBF1/66oqNDLL7+sxMREJScnS5I2btyobdu26Z577jFVaIcOHfTkk0+qa9eukqTXX39dV199tXJzc5WUlNTg83bs2OExsLZt25paLoCmJX1Zvl5dU+DV7jZU1x6IoGHVcrPzCjUuI0fH/mVZVFqhcRk5TFzoR7zX4cn0MTd33nmnXC6XHnvsMY/21NRU7du3T3Pnzj2pgmJjY/W3v/1Nd9xxh9djtVtufv75Z7Vu3bpRr8+WG6BpOVLl1rkPf+ix5eRYEQ7pm8eG+XVXkVXLrXYb+t1Tn3psRTiaQzUzM6+b9nt2m5wk3mt7MfP7bfp/7DvvvKPRo0d7tY8aNUqZmZlmX65OdXW1Fi5cqEOHDtVtEWpI79695XK5NGjQIK1cufK4fSsrK1VWVuZxA9B0vLFh93EDhlSzJeWNDbttsdxNBSUN/thKkiGpsLRCmwpK/LrccMR7Hb5Mh5sWLVpo3bp1Xu3r1q1TdHR0Pc84vq1bt+q0005T8+bNNXbsWC1ZsqTByQBdLpdmz56tzMxMLV68WN26ddOgQYO0Zs2aBl8/PT1dTqez7taxY0fTNQIInD0lh/3ar6kvt7i84R/bxvRDw3ivw5fps6Xuu+8+jRs3Tl9++aX69+8vqeaYm7lz5+qRRx4xXUC3bt20ZcsWHThwQJmZmbrlllu0evXqegNOt27d1K1bt7r7ycnJ2rdvn5555hkNHDiw3tefPn26pkyZUne/rKyMgAM0IZ1iW/q1X1NfblyMb38E+toPDeO9Dl+mt9w88MADWrBggXJzczVx4kRNnDhRubm5mj9/vh544AHTBTRr1kxdu3ZV3759lZ6erl69eun555/3+fn9+/fXzp07G3y8efPmatWqlccNQNNxc3JnnehwhwhHTT87LLdfQqxczmg1tGiHas7k6ZcQ69flhiPe6/DVqKPkRo4cqfXr16ukpEQlJSVav369Ro4c6ZeCDMNQZWWlz/1zc3PlcnGkOxCqmkVFaMyAhOP2GTMgwe/zzli13MgIh1JTarZMH/ujW3s/NSWRA1z9gPc6fJ3UJH4n68EHH9SwYcPUsWNHlZeXa+HChVq1apWys7Ml1exS2r9/vxYsWCBJmjlzpjp37qykpCQdOXJEGRkZyszMPKkDmQFYr/Z062DPN2PVcof2cGnWqD5ec6/EM/eK3/Fehyefws3pp5/u8yR+JSW+H3X+ww8/6Oabb1ZhYaGcTqd69uyp7OxsDR48WJJUWFiovXv31vU/cuSIpk6dqv3796tFixZKSkrSBx98oOHDh/u8TABN0/Thibr/inODPlOwVcsd2sOlwYnxzJobBLzX4ceneW5ef/31un//9NNPevzxxzVkyJC6U7Y3bNig5cuX6+GHH9bkyZMDV60fMM8NAAChx8zvt+lJ/K699lpddtlluvfeez3aX3zxRX388cd69913TRccTIQbAABCT0An8Vu+fLmGDh3q1T5kyBB9/PHHZl8OAADAr0yHmzZt2mjJkiVe7e+++67atGnjl6IAAAAay/TZUmlpabrjjju0atUqjwtnZmdn65///KffCwQAADDDdLi59dZb1b17d73wwgtavHixDMNQYmKi1q9fr9/+9reBqBGoU+02wuqMhyNV7qCfxWMlq8Zr1XoVbuuzFH7rNKxh+oDiUMcBxaErO6/Qa64Kl43nqkhflh/0+VesZNV4rVqvwm19lsJvnYZ/BfSAYkn67rvv9Ne//lU33nijiouLJUnZ2dnatm1bY14OOKHsvEKNy8jxusJvUWmFxmXkKDuv0KLKAiN9Wb5eXVPgddVqtyG9uqZA6cvyrSksQKwar1XrVbitz1L4rdOwlulws3r1ap133nn6/PPPlZmZqYMHD0qSvv76a6Wmpvq9QKDabSgtK1/1bWKsbUvLylf1sd+aIepIlVtz1hYct8+ctQU6UuUOUkWBZdV4rVqvwm19lsJvnYb1GnXhzMcff1wrVqxQs2bN6tovu+wybdiwwa/FAZK0qaDE6y/coxmSCksrtKnA99mxm7I3Nuz2+uv2WG6jpp8dWDVeq9arcFufpfBbp2E90+Fm69at+uMf/+jV3rZtW/30009+KQo4WnF5wz8EjenX1O0pOezXfk2dVeO1ar0Kt/VZCr91GtYzHW5at26twkLv/cG5ubk644wz/FIUcLS4mGi/9mvqOsW29Gu/ps6q8Vq1XoXb+iyF3zoN65kONzfeeKOmTZumoqIiORwOud1urV+/XlOnTtXo0aMDUSPCXL+EWLmc0WroBFmHas4y6ZcQG8yyAubm5M460dnAEY6afnZg1XitWq/CbX2Wwm+dhvVMh5snnnhCZ555ps444wwdPHhQiYmJGjhwoC666CL99a9/DUSNCHOREQ6lptScJnrs92Pt/dSURNvMD9IsKkJjBiQct8+YAQm2mRvEqvFatV6F2/oshd86DeuZmufGMAzt3btXbdu2VVFRkXJycuR2u9W7d2+dffbZgazTb5jnJnSF27wg4TYnCPPc2Ht9lsJvnYZ/Beyq4G63W9HR0dq2bVvIhJljEW5CW7jN6Bpus7kyQ7G912cp/NZp+E/Awo0kJSUl6bXXXlP//v1PqkirEG4AAAg9AZ2h+Omnn9af//xn5eXlNbpAAACAQDG95eb000/X4cOHVVVVpWbNmqlFixYej5eUNO2Jp9hyAwBA6DHz+236quAzZ85sbF0AAAABZzrc3HLLLYGoAwAAwC9MhxtJqq6u1pIlS7R9+3Y5HA51795dV199taKiGvVyACCJs5YA+IfpNJKXl6err75aRUVF6tatmyTp22+/Vdu2bbV06VKdd955fi8SgP0x3wwAfzF9QHH//v0VFxen119/Xaeffrok6eeff9att96q4uLiJn9lcA4oBpqe7LxCjcvI0bFfRrXbTmaN6hOQoGHVcgGYF9BTwb/66iulp6fXBRup5gyqJ554Qlu2bDFdLIDwVu02lJaV7xUwJNW1pWXlq9pt6u+wJrtcAIFnOtx069ZNP/zwg1d7cXGxunbt6peiAISPTQUlHruEjmVIKiyt0KYC/04zYdVyAQSeT+GmrKys7jZjxgxNnDhRixYt0vfff6/vv/9eixYt0n333aennnoq0PUCsJni8oYDRmP6NfXlAgg8nw4obt26tRyO/ztzwDAMjRw5sq6t9rCdlJQUVVdXB6BMAHYVFxPt135NfbkAAs+ncLNy5cpA1wEgTPVLiJXLGa2i0op6j39xSIp31pyebYflAgg8n8LNJZdcEug6AISpyAiHUlMSNS4jRw7JI2jUbi9OTUn0+7wzVi0XQOBxnXkAlhvaw6VZo/oo3um5CyjeGR3Q07GtWi6AwDI9z02oY54boOlihmIADQnohTMBIFAiIxxK7tImbJYLIDDYLQUAAGyFcAMAAGzFp91SvXv39pjn5nhycnJOqiAAAICT4VO4GTFiRN2/Kyoq9PLLLysxMVHJycmSpI0bN2rbtm265557AlIkAACAr3wKN6mpqXX/vvPOOzVx4kQ99thjXn327dvn3+pwQuF2lodV4z1S5dYbG3ZrT8lhdYptqZuTO6tZFHt1/S3c1mcAgWH6VHCn06nNmzfr7LPP9mjfuXOn+vbtq9LSUr8W6G92OhU8O69QaVn5Hhf/czmjlZqSaMv5Oawab/qyfM1ZW6CjLw4d4ZDGDEjQ9OGJAVtuuAm39RmAOWZ+v03/6dmiRQutW7fOq33dunWKjuYaLMGSnVeocRk5Xlc1Liqt0LiMHGXnFVpUWWBYNd70Zfl6dY1nsJEktyG9uqZA6cvyA7LccBNu6zOAwDI9z819992ncePG6csvv1T//v0l1RxzM3fuXD3yyCN+LxDeqt2G0rLy670ejqGaqePTsvI1ODHeFpv0rRrvkSq35qwtOG6fOWsLdP8V57KL6iSE2/oMIPBMfyM/8MADWrBggXJzczVx4kRNnDhRubm5mj9/vh544IFA1IhjbCoo8foL92iGpMLSCm0qKAleUQFk1Xjf2LDba4vNsdxGTT80XritzwACr1EzFI8cOVIjR470dy3wUXF5wz8EjenX1Fk13j0lh/3aD/ULt/UZQOCxLT0ExcX4dmyTr/2aOqvG2ym2pV/7oX7htj4DCDzT4aa6ulrPPPOM+vXrp/j4eMXGxnrcEHj9EmLlckaroaMPHKo5y6Rfgj0+D6vGe3NyZ53oEI8IR00/NF64rc8AAs90uElLS9Nzzz2nkSNHqrS0VFOmTNE111yjiIgIPfroowEoEceKjHAoNaXmFORjfxBq76emJNrm4EurxtssKkJjBiQct8+YAQkcTHySwm19BhB4pr+V33zzTc2ZM0dTp05VVFSU/vSnP+mf//ynHnnkEW3cuDEQNaIeQ3u4NGtUH8U7PTfVxzujNWtUH9vNC2LVeKcPT9TdAxO8tuBEOKS7BzLPjb+E2/oMILBMT+J36qmnavv27TrzzDPlcrn0wQcfqE+fPtq1a5d69+7NJH5BFm4zujJDsb2F2/oMwHdmfr9Nny3VoUMHFRYW6swzz1TXrl310UcfqU+fPvriiy/UvHnzRheNxomMcCi5Sxurywgaq8bbLCpCdww4K+jLDTfhtj4DCAzTf3r+8Y9/1CeffCJJmjRpkh5++GGdffbZGj16tG6//Xa/FwgAAGCG6d1Sx9q4caM+++wzde3aVX/4wx/8VVfA2G23FAAA4SCgu6WO1b9//7rLMAAAAFjNp3CzdOlSn18wFLbeAAAA+/Ip3IwYMcKnF3M4HKqurj6ZeoDjCrezaTg7DADM8yncuN3uQNcBnFB2XqHSsvI9LrLockYrNSXRlvOgWDXe9GX5mrO2wOOioU8s264xA5jXB0Bo4E8xhITsvEKNy8jxunp0UWmFxmXkKDuv0KLKAsOq8aYvy9erawq8robuNqRX1xQofVl+QJYLAP7UqHDzySef6KqrrlKXLl3UtWtXXXXVVfr444/9XRsgqWbXTFpWvuo7ra+2LS0rX9XH/iKHKKvGe6TKrTlrC47bZ87aAh2pYksugKbNdLh58cUXNXToUMXExGjSpEmaOHGiWrVqpeHDh+vFF18MRI0Ic5sKSry2YBzNkFRYWqFNBSXBKyqArBrvGxt2e22xOZbbqOkHAE2Z6VPB09PT9fe//1333ntvXdvEiRN18cUX64knnvBoB/yhuLzhH/rG9GvqrBrvnpLDfu0HAFYxveWmrKxMQ4cO9Wq/4oorVFZW5peigKPFxUSfuJOJfk2dVePtFNvSr/0AwCqmw80f/vAHLVmyxKv9vffeU0pKil+KAo7WLyFWLme0GjoB2qGas4j6JcQGs6yAsWq8Nyd39rr6+bEiHDX9AKAp82m31AsvvFD37+7du+uJJ57QqlWrlJycLKnmEgzr16/X/fffH5gqEdYiIxxKTUnUuIwcOSSPA21rf4tTUxJtM9+NVeNtFhWhMQMS9Oqahg8qHjMggfluADR5Pl1bKiEhwbcXczi0a9euky4qkLi2VOhinhvr5rmJcIh5bgBYyszv90lfODPUEG5CGzMUM0MxgPBEuDkOwg0AAKEnoFcFr66u1vz58/XJJ5+ouLjY69IMn376qdmXBAAA8BvT4WbSpEmaP3++rrzySvXo0UMOh313CQAAgNBjOtwsXLhQ//rXvzR8+PBA1AMAAHBSTB8h2KxZM3Xt2jUQtQAAAJw001tu7r//fj3//PN68cUX2SV1FKvOagm3s4fCjVVnLbE+AwhlpsPNunXrtHLlSn344YdKSkrSKaec4vH44sWLfX6tWbNmadasWdq9e7ckKSkpSY888oiGDRvW4HNWr16tKVOmaNu2bWrfvr3+8pe/aOzYsWaH4VdWzUcSbvO+hJv65pt5Ytn2gM83w/oMINSZ/hOwdevW+uMf/6hLLrlEv/nNb+R0Oj1uZnTo0EFPPvmkNm/erM2bN+v3v/+9rr76am3btq3e/gUFBRo+fLgGDBig3NxcPfjgg5o4caIyMzPNDsNvsvMKNS4jx+sqzkWlFRqXkaPsvEJbLRfBkb4sX6+uKfC6SrfbkF5dU6D0ZfkBWS7rMwA7aHLz3MTGxupvf/ub7rjjDq/Hpk2bpqVLl2r79u11bWPHjtVXX32lDRs2+PT6/pznptpt6HdPfer1hVzLISneGa11037v103rVi0XwXGkyq1zH/7QK9gcLcIhffPYML/uomJ9BtCUmfn9bvQ343//+1+tW7dO69ev13//+9/Gvkyd6upqLVy4UIcOHaq7ZtWxNmzYoCuuuMKjbciQIdq8ebN+/fXXep9TWVmpsrIyj5u/bCooafALWaq5JlBhaYU2FZT4bZlWLhfB8caG3ccNNlLNFpw3Nuz263JZnwHYhelwc+jQId1+++1yuVwaOHCgBgwYoPbt2+uOO+7Q4cOHTRewdetWnXbaaWrevLnGjh2rJUuWKDGx/uMJioqK1K5dO4+2du3aqaqqSj/++GO9z0lPT/fYbdaxY0fTNTakuLzhL+TG9Gvqy0Vw7Cnx7f+Rr/18xfoMwC5Mh5spU6Zo9erVysrK0oEDB3TgwAG99957Wr16daOuCt6tWzdt2bJFGzdu1Lhx43TLLbcoP7/h4wmOPUOrdq9aQ2duTZ8+XaWlpXW3ffv2ma6xIXEx0X7t19SXi+DoFNvSr/18xfoMwC5Mh5vMzEy99tprGjZsmFq1aqVWrVpp+PDhmjNnjhYtWmS6gNp5c/r27av09HT16tVLzz//fL194+PjVVRU5NFWXFysqKgotWnTpt7nNG/evK7O2pu/9EuIlcsZrYaOAnCo5myPfgmxflumlctFcNyc3FknOrQkwlHTz59YnwHYhelwc/jwYa9dQ5IUFxfXqN1SxzIMQ5WVlfU+lpycrBUrVni0ffTRR+rbt6/XKenBEBnhUGpKzS60Y7+Ya++npiT6/SBIq5aL4GgWFaExAxKO22fMgAS/z3fD+gzALkx/OyYnJys1NVUVFf+3//uXX35RWlpagwcCN+TBBx/U2rVrtXv3bm3dulUPPfSQVq1apZtuuklSzS6l0aNH1/UfO3as9uzZoylTpmj79u2aO3euXnvtNU2dOtXsMPxmaA+XZo3qo3in5ybzeGe0Zo3qE7D5OaxaLoJj+vBE3T0wwWsLToRDuntg4Oa5YX0GYAemTwXfunWrhg0bpoqKCvXq1UsOh0NbtmxRdHS0li9frqSkJJ9f64477tAnn3yiwsJCOZ1O9ezZU9OmTdPgwYMlSbfeeqt2796tVatW1T1n9erVmjx5ct0kftOmTTM1iZ8/TwU/GjO6IhCYoZj1GUANM7/fjZrn5pdfflFGRoa++eYbGYahxMRE3XTTTWrRokWjiw6WQIUbAAAQOGZ+v01dfuHXX39Vt27d9P7772vMmDEnVSQAAEAgmNq+fcopp6iyspILZgIAgCbL9M77CRMm6KmnnlJVVVUg6gEAADgppq8K/vnnn+uTTz7RRx99pPPOO0+nnnqqx+NmrgoOAADgb6bDTevWrXXttdcGohY0AmeXAADgyXS4mTdvXiDqQCNk5xUqLSvf46KDLme0UlMSmRcEABC2TB9zk5GR0eBjf/7zn0+qGPguO69Q4zJyvK6mXFRaoXEZOcrOK7SoMgAArGU63Nx77716//33vdonT5583OAD/6l2G0rLyld9ExTVtqVl5avabXoKIwAAQp7pcLNw4UKNGjVKa9asqWubMGGC/vWvf2nlypV+LQ7121RQ4rXF5miGpMLSCm0qKAleUQAANBGmw83QoUP1yiuvaMSIEdq8ebPuueceLV68WCtXrtS5554biBpxjOLyhoNNY/oBAGAnpg8olqQbbrhBP//8s373u9+pbdu2Wr16tbp27erv2tCAuJjoE3cy0Q8AADvxKdxMmTKl3va4uDj17t1bL7/8cl3bc88955/K0KB+CbFyOaNVVFpR73E3DtVcTblfQmywSwMAwHI+hZvc3Nx627t06aKysrK6x7ksQ3BERjiUmpKocRk5ckgeAaf2E0hNSWS+GwBAWGrUVcFDmZ2uCs48NwCAcBGwq4KjaRnaw6XBifHMUAwAwFEINyEuMsKh5C5trC4DAIAmw/Sp4AAAAE0Z4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANgK4QYAANiKpeEmPT1dF154oWJiYhQXF6cRI0Zox44dx33OqlWr5HA4vG7ffPNNkKoGAABNmaXhZvXq1Ro/frw2btyoFStWqKqqSldccYUOHTp0wufu2LFDhYWFdbezzz47CBUDAICmLsrKhWdnZ3vcnzdvnuLi4vTll19q4MCBx31uXFycWrduHcDqAABAKGpSx9yUlpZKkmJjY0/Yt3fv3nK5XBo0aJBWrlzZYL/KykqVlZV53AAAgH01mXBjGIamTJmi3/3ud+rRo0eD/Vwul2bPnq3MzEwtXrxY3bp106BBg7RmzZp6+6enp8vpdNbdOnbsGKghAACAJsBhGIZhdRGSNH78eH3wwQdat26dOnToYOq5KSkpcjgcWrp0qddjlZWVqqysrLtfVlamjh07qrS0VK1atTrpugEAQOCVlZXJ6XT69PvdJLbcTJgwQUuXLtXKlStNBxtJ6t+/v3bu3FnvY82bN1erVq08bgAAwL4sPaDYMAxNmDBBS5Ys0apVq5SQkNCo18nNzZXL5fJzdQAAIBRZGm7Gjx+vt956S++9955iYmJUVFQkSXI6nWrRooUkafr06dq/f78WLFggSZo5c6Y6d+6spKQkHTlyRBkZGcrMzFRmZqZl4wAAAE2HpeFm1qxZkqRLL73Uo33evHm69dZbJUmFhYXau3dv3WNHjhzR1KlTtX//frVo0UJJSUn64IMPNHz48GCVDQAAmrAmc0BxsJg5IAkAADQNIXdAMQAAgL8QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK0QbgAAgK1EWV0AQlO129CmghIVl1coLiZa/RJiFRnhsLosAAAINzAvO69QaVn5KiytqGtzOaOVmpKooT1cFlYGAAC7pWBSdl6hxmXkeAQbSSoqrdC4jBxl5xVaVBkAADUIN/BZtdtQWla+jHoeq21Ly8pXtbu+HgAABAfhBj7bVFDitcXmaIakwtIKbSooCV5RAAAcg3ADnxWXNxxsGtMPAIBAINzAZ3Ex0X7tBwBAIBBu4LN+CbFyOaPV0AnfDtWcNdUvITaYZQEA4IFwA59FRjiUmpIoSV4Bp/Z+akoi890AACxFuIEpQ3u4NGtUH8U7PXc9xTujNWtUH+a5AQBYjkn8YNrQHi4NToxnhmIAQJNEuEGjREY4lNyljdVlAADghd1SAADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVgg3AADAVsJuhmLDMCRJZWVlFlcCAAB8Vfu7Xfs7fjxhF27Ky8slSR07drS4EgAAYFZ5ebmcTudx+zgMXyKQjbjdbv3nP/9RTEyMHA57XOixrKxMHTt21L59+9SqVSurywk4xmtvjNf+wm3MjNc/DMNQeXm52rdvr4iI4x9VE3ZbbiIiItShQwerywiIVq1ahcV/nFqM194Yr/2F25gZ78k70RabWhxQDAAAbIVwAwAAbIVwYwPNmzdXamqqmjdvbnUpQcF47Y3x2l+4jZnxBl/YHVAMAADsjS03AADAVgg3AADAVgg3AADAVgg3AADAVgg3ISY9PV0Oh0P33XefT/3Xr1+vqKgonX/++QGtK1B8HW9lZaUeeughderUSc2bN1eXLl00d+7c4BTpR76O980331SvXr3UsmVLuVwu3Xbbbfrpp5+CU+RJePTRR+VwODxu8fHxx33O6tWrdcEFFyg6OlpnnXWWXnnllSBVe/LMjnfx4sUaPHiw2rZtq1atWik5OVnLly8PYsUnrzGfca1Q/L5qzHhD+fuqMeO14vsq7GYoDmVffPGFZs+erZ49e/rUv7S0VKNHj9agQYP0ww8/BLg6/zMz3pEjR+qHH37Qa6+9pq5du6q4uFhVVVVBqNJ/fB3vunXrNHr0aP39739XSkqK9u/fr7Fjx+rOO+/UkiVLglRt4yUlJenjjz+uux8ZGdlg34KCAg0fPlxjxoxRRkaG1q9fr3vuuUdt27bVtddeG4xyT5qZ8a5Zs0aDBw/WjBkz1Lp1a82bN08pKSn6/PPP1bt372CU6xdmxlwrlL+vzI431L+vzIzXqu8rwk2IOHjwoG666SbNmTNHjz/+uE/Pufvuu3XjjTcqMjJS7777bmAL9DMz483Oztbq1au1a9cuxcbGSpI6d+4chCr9x8x4N27cqM6dO2vixImSpISEBN199916+umng1HqSYuKivL5L/lXXnlFZ555pmbOnClJ6t69uzZv3qxnnnkmZMKNmfHWjrPWjBkz9N577ykrKyukwo2ZMdcK5e8rM+O1w/eVmfFa9X3FbqkQMX78eF155ZW6/PLLfeo/b948fffdd0pNTQ1wZYFhZrxLly5V37599fTTT+uMM87QOeeco6lTp+qXX34JQqX+YWa8F110kb7//nstW7ZMhmHohx9+0KJFi3TllVcGodKTt3PnTrVv314JCQm64YYbtGvXrgb7btiwQVdccYVH25AhQ7R582b9+uuvgS7VL8yM91hut1vl5eV1P4KhwuyYQ/37ysx47fB9ZWa8Vn1fseUmBCxcuFA5OTn64osvfOq/c+dOPfDAA1q7dq2iokLvIzY73l27dmndunWKjo7WkiVL9OOPP+qee+5RSUlJSOzHNjveiy66SG+++aauv/56VVRUqKqqSn/4wx/0j3/8I8CVnrzf/va3WrBggc455xz98MMPevzxx3XRRRdp27ZtatOmjVf/oqIitWvXzqOtXbt2qqqq0o8//iiXyxWs0hvF7HiP9eyzz+rQoUMaOXJkEKr1D7NjDvXvK7PjDfXvK7Pjtez7ykCTtnfvXiMuLs7YsmVLXdsll1xiTJo0qd7+VVVVRt++fY1Zs2bVtaWmphq9evUKcKX+YXa8hmEYgwcPNqKjo40DBw7UtWVmZhoOh8M4fPhwIMs9aY0Z77Zt2wyXy2U8/fTTxldffWVkZ2cb5513nnH77bcHoWL/OnjwoNGuXTvj2Wefrffxs88+25gxY4ZH27p16wxJRmFhYTBK9KsTjfdob731ltGyZUtjxYoVQagscI435lD/vqrPiT7jUP6+qs+JxmvV9xXhpolbsmSJIcmIjIysu0kyHA6HERkZaVRVVXn0//nnn736OxyOurZPPvnEopH4xux4DcMwRo8ebXTp0sWjLT8/35BkfPvtt8EqvVEaM95Ro0YZ1113nUfb2rVrDUnGf/7zn2CV7jeXX365MXbs2HofGzBggDFx4kSPtsWLFxtRUVHGkSNHglGe3x1vvLUWLlxotGjRwnj//feDVFVgNTTmUP++asjxPuNQ/r5qyPHGa9X3VehtAwwzgwYN0tatWz3abrvtNp177rmaNm2a11HqrVq18ur/8ssv69NPP9WiRYuUkJAQ8JpPhtnxStLFF1+sd955RwcPHtRpp50mSfr2228VERGhDh06BKXuxmrMeA8fPuy1+b62nxFil4qrrKzU9u3bNWDAgHofT05OVlZWlkfbRx99pL59++qUU04JRol+daLxStLbb7+t22+/XW+//XbIHEd1PMcbc6h/X9XnRJ9xKH9f1edE47Xs+ypgsQkBc+xuiwceeMC4+eabG+wf6pt5TzTe8vJyo0OHDsZ1111nbNu2zVi9erVx9tlnG3feeacF1Z68E4133rx5RlRUlPHyyy8b3333nbFu3Tqjb9++Rr9+/Syo1pz777/fWLVqlbFr1y5j48aNxlVXXWXExMQYu3fvNgzDe6y7du0yWrZsaUyePNnIz883XnvtNeOUU04xFi1aZNUQTDE73rfeesuIiooyXnrpJaOwsLDudvQujKbO7JiPFWrfV2bHG+rfV2bHa9X3FVtubKCwsFB79+61uoygOXa8p512mlasWKEJEyaob9++atOmjUaOHOnzKfNN3bHjvfXWW1VeXq4XX3xR999/v1q3bq3f//73euqppyys0jfff/+9/vSnP+nHH39U27Zt1b9/f23cuFGdOnWS5D3WhIQELVu2TJMnT9ZLL72k9u3b64UXXgiZ08DNjvfVV19VVVWVxo8fr/Hjx9e133LLLZo/f36wy28Us2MOdWbHG+rfV2bHa9X3lcMwQmw7NgAAwHEwzw0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0AALAVwg0QQi699FLdd999liy7c+fOmjlz5nH7OBwOvfvuuye1nEcffVTnn39+g4+vWrVKDodDBw4cOKnl+IM/ann00UflcDjkcDhO+P6erEsvvbRuWVu2bAnosgArcfkFALBYUlKSPv74Y7Vq1Sqgy1m8eLG+++479evXL6DLAaxGuAFs7MiRI2rWrJnVZeAEoqKiFB8fH/DlxMbGqqysLODLAazGbikghGVnZ8vpdGrBggWSai5SN2LECKWnp6t9+/Y655xzJElbt27V73//e7Vo0UJt2rTRXXfdpYMHD9a9Tu3znnnmGblcLrVp00bjx4/Xr7/+6rG88vJy3XjjjTrttNPUvn17/eMf/zhufdOmTdM555yjli1b6qyzztLDDz/s9ZpPPvmk2rVrp5iYGN1xxx2qqKgw/T589tlnGjhwoFq0aKGOHTtq4sSJOnToUN3jnTt31owZM3T77bcrJiZGZ555pmbPnu31Gueff76io6PVt29fvfvuu167b5YtW6ZzzjlHLVq00GWXXabdu3ebrsVXBw4c0F133aV27dopOjpaPXr00Pvvvy9Jmj9/vlq3bq33339f3bp1U8uWLXXdddfp0KFDev3119W5c2edfvrpmjBhgqqrq00vGwh1hBsgRC1cuFAjR47UggULNHr06Lr2Tz75RNu3b9eKFSv0/vvv6/Dhwxo6dKhOP/10ffHFF3rnnXf08ccf69577/V4vZUrV+q7777TypUr9frrr2v+/PleV6L+29/+pp49eyonJ0fTp0/X5MmTtWLFigZrjImJ0fz585Wfn6/nn39ec+bM0d///ve6x//1r38pNTVVTzzxhDZv3iyXy6WXX37Z1PuwdetWDRkyRNdcc42+/vpr/e///q/WrVvnNb5nn31Wffv2VW5uru655x6NGzdO33zzjaSa0JaSkqLzzjtPOTk5euyxxzRt2jSP5+/bt0/XXHONhg8fri1btujOO+/UAw880KhaTsTtdmvYsGH67LPPlJGRofz8fD355JOKjIys63P48GG98MILWrhwobKzs7Vq1Spdc801WrZsmZYtW6Y33nhDs2fP1qJFi0wtG7AFA0DIuOSSS4xJkyYZL730kuF0Oo1PP/3U4/FbbrnFaNeunVFZWVnXNnv2bOP00083Dh48WNf2wQcfGBEREUZRUVHd8zp16mRUVVXV9fmf//kf4/rrr6+736lTJ2Po0KEey7v++uuNYcOG1d2XZCxZsqTB+p9++mnjggsuqLufnJxsjB071qPPb3/7W6NXr14NvsbKlSsNScbPP/9sGIZh3HzzzcZdd93l0Wft2rVGRESE8csvv9TVPmrUqLrH3W63ERcXZ8yaNcswDMOYNWuW0aZNm7r+hmEYc+bMMSQZubm5hmEYxvTp043u3bsbbre7rs+0adNM13Ks1NRUr/EuX77ciIiIMHbs2FHvc+bNm2dIMv7973/Xtd19991Gy5YtjfLy8rq2IUOGGHfffbfHcwsKCjzGBdgRx9wAISYzM1M//PCD1q1bV++Boeedd57HcTbbt29Xr169dOqpp9a1XXzxxXK73dqxY4fatWsnqeag1qO3DLhcLm3dutXjtZOTk73uH+8Mn0WLFmnmzJn697//rYMHD6qqqsrjoNnt27dr7NixXq+5cuXK47wDnr788kv9+9//1ptvvlnXZhiG3G63CgoK1L17d0lSz5496x53OByKj49XcXGxJGnHjh3q2bOnoqOj6/oc+95u375d/fv3l8Ph8Ki1MbWcyJYtW9ShQ4e63Yr1admypbp06VJ3v127durcubNOO+00j7baMQLhhHADhJjzzz9fOTk5mjdvni688EKPH1tJHiFGqvlxPbZPraPbTznlFK/H3G73Cetp6LU3btyoG264QWlpaRoyZIicTqcWLlyoZ5999oSvaYbb7dbdd9+tiRMnej125pln1v37eOOr7z0yDOO490+mlhNp0aLFCfvUN57GfoaA3RBugBDTpUsXPfvss7r00ksVGRmpF1988bj9ExMT9frrr+vQoUN1wWf9+vWKiIg47paB+mzcuNHr/rnnnltv3/Xr16tTp0566KGH6tr27Nnj0ad79+7auHGjxzFDxy7jRPr06aNt27apa9eupp53tHPPPVdvvvmmKisr1bx5c0nS5s2bPfokJiZ6zeFzbK3+qEWq2cr0/fff69tvvzX9GQHggGIgJJ1zzjlauXKlMjMzTzip30033aTo6GjdcsstysvL08qVKzVhwgTdfPPNdbukfLV+/Xo9/fTT+vbbb/XSSy/pnXfe0aRJk+rt27VrV+3du1cLFy7Ud999pxdeeEFLlizx6DNp0iTNnTtXc+fO1bfffqvU1FRt27bNVE3Tpk3Thg0bNH78eG3ZskU7d+7U0qVLNWHCBJ9f48Ybb5Tb7dZdd92l7du3a/ny5XrmmWck/d+WqbFjx+q7777TlClTtGPHDr311lteB1z7oxZJuuSSSzRw4EBde+21WrFihQoKCvThhx8qOzvb1OsA4YpwA4Sobt266dNPP9Xbb7+t+++/v8F+LVu21PLly1VSUqILL7xQ1113nQYNGnTCLT71uf/++/Xll1+qd+/eeuyxx/Tss89qyJAh9fa9+uqrNXnyZN177706//zz9dlnn+nhhx/26HP99dfrkUce0bRp03TBBRdoz549GjdunKmaevbsqdWrV2vnzp0aMGCAevfurYcfflgul8vn12jVqpWysrK0ZcsWnX/++XrooYf0yCOPSFLdcThnnnmmMjMzlZWVpV69eumVV17RjBkz/F5LrczMTF144YX605/+pMTERP3lL3/htG7ARw7Dlx3JABBm3nzzTd12220qLS316RiYxnr00Uf17rvvBu1yCLt371ZCQoJyc3OPe5kLIJSx5QYAJC1YsEDr1q1TQUGB3n33XU2bNk0jR44MaLCptXXrVp122mmm5/gxa9iwYUpKSgroMoCmgC03ACDp6aef1ssvv6yioiK5XC6NGDFCTzzxhFq2bBnQ5ZaUlKikpESS1LZtWzmdzoAta//+/frll18k1exm49IcsCvCDQAAsBV2SwEAAFsh3AAAAFsh3AAAAFsh3AAAAFsh3AAAAFsh3AAAAFsh3AAAAFsh3AAAAFv5f1Nh7Caswy6FAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(df[df[\"species\"] == \"setosa\"].sepal_length, df[df[\"species\"] == \"setosa\"].sepal_width)\n",
    "plt.xlabel(\"begerblad lengde [cm]\");\n",
    "plt.ylabel(\"begerblad bredde [cm]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70e7587f",
   "metadata": {},
   "source": [
    "Ganske riktig! Det ser ut til å være en ganske klar sammenheng!\n",
    "\n",
    "I prinsippet kan vi utføre multidimensjonal lineær regresjon. I dette tilfellet prøver vi derimot å predikere hva bredden til et kronblad vil være basert på en måling av lengden. Dette er det vi kaller enkel lineær regresjon hvor vi skal forutsi en respons (kronbladets bredde) basert på en enkelt predikator (kronbladets lengde). Vi antar altså at det er en lineær sammenheng mellom disse to verdiene som kan beskrives matematisk slik \n",
    "\n",
    "$Y \\approx \\beta_0 + \\beta_1 X + \\epsilon$\n",
    "\n",
    "Hvor Y i dette tilfellet er kronbladets bredde, X er det vi måler, altså kronbladets lengde. De to konstantene får vi ved å gjøre en tilpasning til treningsdataene våre ved å bruke *minste kvadraters metode* dvs. å minimere *residual sum of squares* (RSS). Verdiene for $\\beta_0$ og $\\beta_1$ som minimerer RSS skriver vi ofte med en hatt, $\\hat{\\beta}_0$ og $\\hat{\\beta}_1$. På samme måte skriver vi $\\hat{y}$ for å spesifisere at vi snakker om vår predksjon av Y. Vi kan da skrive\n",
    "\n",
    "$\\hat{y} = \\hat{\\beta}_0 + \\hat{\\beta}_1 x$,\n",
    "\n",
    "hvor $x = X$\n",
    "\n",
    "Siden funksjonen vi vil bruke for å tilpasse en enkel lineær regresjons til dataene våre kun inneholder én predikator så må vi legge til en ekstra dimensjon for at funksjonen ikke klager. I praksis vil det si å konvertere en vektor til en matrise med èn ekstra dimensjon. Man kan enkelt sjekke dimensjonen til et numpy-array object i pyhon ved å kalle på *shape*, den gir deg eksakt dimensjonene av arrayen. Før vi legger til en dimensjon: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "76009ba8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50,)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_setosa_sepal_length.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a081ed",
   "metadata": {},
   "source": [
    "Vi legger til en ekstra dimensjon/akse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "95fe5ea0",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_setosa_sepal_length = X_setosa_sepal_length[:, np.newaxis]\n",
    "X_setosa_sepal_width = X_setosa_sepal_width[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd9b518e",
   "metadata": {},
   "source": [
    "... og vi får en ny dimensjon i arrayen vår. Dette vil bli mer intuitivt når vi skal se på tilfeller med flere predikatorer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d980aaa8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50, 1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_setosa_sepal_length.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ced9e6a0",
   "metadata": {},
   "source": [
    "## Lineær Regresjon\n",
    "\n",
    "scikit-learn har mange nyttige funksjoner for ulike typer maskinlæringsalgoritmer. De er i tillegg relativt godt dokumenterte - så her er det bare å slå seg løs. Vi skal nå bruke [*LinearRegression()*](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html). Den tar enkelt og greit som input repsons og predikatoren i form av numpy-arrayer. Funksjonen *fit()* setter i gang prosessen med å finne den beste tilpasningen til dataene våre. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2717e360",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = LinearRegression().fit(X_setosa_sepal_length,X_setosa_sepal_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a90542a2",
   "metadata": {},
   "source": [
    "Når den er ferdig kan vi hente ut veridene av $\\hat{\\beta}_0$ (intercept) og $\\hat{\\beta}_1$ (coef). Vi kan også få en første indikasjon på hvor bra tilpasningen er ved å skrive ut $R^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1a347705",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient = 0.80\n",
      "Intercept = -0.57\n",
      "Score/R^2 = 0.55\n"
     ]
    }
   ],
   "source": [
    "print(\"Coefficient = %.2f\"%clf.coef_)\n",
    "print(\"Intercept = %.2f\"%clf.intercept_)\n",
    "print(\"Score/R^2 = %.2f\"%clf.score(X_setosa_sepal_length,X_setosa_sepal_width))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa717442",
   "metadata": {},
   "source": [
    "Eller \n",
    "\n",
    "$\\hat{y} = -0.57 + 0.8x$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c81889f",
   "metadata": {},
   "source": [
    "Husk at score (eller $R^2$) er definert som\n",
    "\n",
    "\\begin{equation}\n",
    "R^2 = 1 - \\frac{RSS}{TSS}\n",
    "\\end{equation}\n",
    "\n",
    "RSS (**r**esidual **s**um of **s**quares) er definert som\n",
    "\n",
    "$$RSS = \\sum_{i=1}^{n}\\left(y_i - \\hat{y}_{i}\\right)^{2}$$\n",
    "\n",
    "og TSS (**t**otal **s**um of **s**quares)\n",
    "\n",
    "$$TSS = \\sum_{i=1}^{n}\\left(y_i - \\bar{y}\\right)^{2},$$\n",
    "\n",
    "hvor $\\bar{y}$ er det globale gjennomsnittet i dataene våre. \n",
    "\n",
    "$R^2$ er et mål på hvor mye av variasjon i Y som kan beskrives ved X. Dersom $R^2$ er nære 1 betyr det at en stor del av variasjonen i Y kan beskrives av regresjonsanalysen. Om den derimot er nære 0 betyr det at regresjonsanalysen i liten grad evner å beskrive variasjonen i responsen (Y). Det betyr\n",
    "\n",
    "1. at den lineære modellen er feil\n",
    "2. at variansen til feilen, $\\epsilon$, $\\sigma^2 = Var(\\epsilon)$ er for høy\n",
    "\n",
    "evt. begge deleer.\n",
    "\n",
    "La oss uansett se hvordan den lineære modellen ser ut når vi plotter den sammen med datapunktene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b9ee8834",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definer en x-vektor som dekker datapunktene (i dette tilfellet mellom 4 og 7 cm)\n",
    "X_plot = np.linspace(4, 7, 300)\n",
    "# Lag så responsfunksjonen ved å bruke parameterne fra modellen vår\n",
    "fit = clf.intercept_[0] + clf.coef_[0]*X_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7213809b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff0941e2e50>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAFzCAYAAAAt54EyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZj0lEQVR4nO3deVxU9foH8M8AsumASrKYILihuKJYornljpqalSWmuF3NNcn7K6qrkhZa3dwq3EHFssU991TQzF1wSURTEDWQxASXGGTm/P6YmOvAzDAzHGb9vF8vXi/nnO8585wzB3088z3PIxEEQQARERERkRVyMHcARERERETGYjJLRERERFaLySwRERERWS0ms0RERERktZjMEhEREZHVYjJLRERERFaLySwRERERWS0ms0RERERktZzMHYCpKRQK/PHHH5BKpZBIJOYOh4iIiIjKEAQBDx48QN26deHgoPveq90ls3/88Qf8/f3NHQYRERERVeDmzZuoV6+ezjF2l8xKpVIAypPj4eFh5miIiIiIqKzCwkL4+/ur8jZd7C6ZLZ1a4OHhwWSWiIiIyILpMyWUD4ARERERkdViMktEREREVovJLBERERFZLbubM6sPQRBQUlICuVxu7lCIbJKjoyOcnJxYHo+IiCqNyWwZxcXFyMnJwePHj80dCpFNc3d3h5+fH5ydnc0dChERWTEms09RKBTIzMyEo6Mj6tatC2dnZ945IhKZIAgoLi7Gn3/+iczMTDRu3LjCgthERETaMJl9SnFxMRQKBfz9/eHu7m7ucIhslpubG6pVq4YbN26guLgYrq6u5g6JiIisFG+HaMC7RERVj79nRERWSBDMHUE5/NeEiIiIiLQTBODIEWDoUGD2bHNHU45Zk9k5c+ZAIpGo/fj6+urcJiUlBe3atYOrqysaNGiAZcuWmShaIiIiIjtSXAysXw+EhQFdugCbNwNffw1Y2EPyZr8z27x5c+Tk5Kh+Lly4oHVsZmYmIiIi0LlzZ6SmpuL999/HtGnTsGnTJhNGbN0kEgm2bt1q7jCIiIjIUuXlAXPnAvXrAyNHAmfP/m9dfj6wYYP5YtPA7Mmsk5MTfH19VT916tTROnbZsmUICAjAokWL0KxZM4wbNw5jxozB559/bsKILVNUVBQGDx5c4bicnBz069fPqPd4+k66k5MTnnnmGXTp0gWLFi2CTCYzaF/JycmQSCS4f/++UbEQERGRyM6fB8aOBQICgFmzgNxczeNWrzZtXBUwezJ79epV1K1bF0FBQXj99ddx/fp1rWOPHTuG3r17qy3r06cPTp8+jSdPnmjcRiaTobCwUO3HFOQKAceu5WNb2m0cu5YPucK8E6aLi4sBAL6+vnBxcTF6P6V30rOzs3Ho0CG8+uqriIuLQ8eOHfHgwQOxwiUiIiJTUCiAHTuAHj2A1q2BNWsAbTeoGjQAFi0C9u0zaYgVMWsy+/zzz2PdunXYu3cvVq5cidzcXHTs2BH5+fkax+fm5sLHx0dtmY+PD0pKSnD37l2N28TFxcHT01P14+/vL/pxlLXnYg5eWHAQb6w8jukb0/DGyuN4YcFB7LmYU+XvXapbt26YMmUKoqOj8cwzz6BXr14A1KcZFBcXY8qUKfDz84OrqysCAwMRFxenc7+ld9Lr1q2Lli1bYurUqUhJScHFixexYMEC1bikpCSEhYVBKpXC19cXw4cPR15eHgAgKysL3bt3BwDUqlULEokEUVFRAIA9e/bghRdeQM2aNeHl5YUBAwbg2rVrIp8dIiIiO/fgAbB0KRAcDLz0EnDwoPax3boBW7cCV64A06cDHh6milIvZk1m+/Xrh6FDh6Jly5bo2bMndu7cCQBYu3at1m3KNjEQ/ikRoa25QUxMDAoKClQ/N2/eFCl6zfZczMFbSWeRU1Cktjy3oAhvJZ01aUK7du1aODk54ejRo1i+fHm59UuWLMH27dvx/fffIyMjA0lJSQgMDDT4fZo2bYp+/fph8+bNqmXFxcWYO3cuzp07h61btyIzM1OVsPr7+6vmOWdkZCAnJweLFy8GADx69AjR0dE4deoUDhw4AAcHBwwZMgQKhcLwE0BERETqMjOBd94B6tUDpk0Dfv9d8zhnZyAqCkhNBQ4dAgYNAhwdTRqqviyqaUL16tXRsmVLXL16VeN6X19f5JaZv5GXlwcnJyd4eXlp3MbFxaVSX6sbQq4QELvjEjRNKBAASADE7riEXiG+cHSo+s5ijRo1wqeffqp1fXZ2Nho3bowXXngBEokE9evXN/q9mjZtin1Pfe0wZswY1Z8bNGiAJUuW4LnnnsPDhw9Ro0YN1K5dGwDg7e2NmjVrqsYOHTpUbb+rV6+Gt7c3Ll26hBYtWhgdHxERkd0SBOCXX5RTBLZuVU4t0MbbG5g0CZg4ESjzbbilMvuc2afJZDKkp6fDz89P4/rw8HDs379fbdm+ffsQFhaGatWqmSJEnU5m3it3R/ZpAoCcgiKczLxnknjCwsJ0ro+KikJaWhqCg4Mxbdo0tWTUUIIgqN0dT01NxaBBg1C/fn1IpVJ069YNgDKB1uXatWsYPnw4GjRoAA8PDwQFBem1HREREZWhqbSWtkS2dWsgMRHIzlbWkrWSRBYwczI7c+ZMpKSkIDMzEydOnMArr7yCwsJCjBo1CoByisDIkSNV4ydOnIgbN24gOjoa6enpWLNmDVavXo2ZM2ea6xDU5D3QnsgaM66yqlevrnN927ZtkZmZiblz5+Lvv//Ga6+9hldeecWo90pPT1clno8ePULv3r1Ro0YNJCUl4dSpU9iyZQuA/z2Ips3AgQORn5+PlStX4sSJEzhx4oRe2xEREdE/dJXWeppEAgweDCQnK6cTjBoFmOjbbDGZdZrBrVu38MYbb+Du3buoU6cOOnTogOPHj6u+7i59ar5UUFAQdu3ahRkzZuCrr75C3bp1sWTJknJfTZuLt1S//vL6jjMFDw8PDBs2DMOGDcMrr7yCvn374t69e6ppAPq4fPky9uzZg5iYGNXru3fvYv78+aoH7k6fPq22jbOzMwBALperluXn5yM9PR3Lly9H586dAQC//PJLpY6PiIjIbpw/DyxerKwDq6tkplSqLME1ZQrQsKHp4qsiZk1mN27cqHN9YmJiuWVdu3bFWW3/wzCz54Jqw8/TFbkFRRrnzUoA+Hq64rkg/RPFqrRw4UL4+fmhTZs2cHBwwA8//ABfX1+1OaxllZSUIDc3FwqFAvn5+UhOTsa8efPQpk0b/Pvf/wYABAQEwNnZGUuXLsXEiRNx8eJFzJ07V20/9evXh0QiwU8//YSIiAi4ubmhVq1a8PLywooVK+Dn54fs7Gy89957VXkKiIiIrJtCAezcqZwPq6siAaAsrTVtGjB6tMVVJKgMi5oza+0cHSSYPTAEgDJxfVrp69kDQ0zy8Jc+atSogQULFiAsLAzt27dHVlYWdu3aBQcH7ZfFb7/9Bj8/PwQEBKBbt274/vvvERMTgyNHjqBGjRoAgDp16iAxMRE//PADQkJCMH/+/HKNLZ599lnExsbivffeg4+PD6ZMmQIHBwds3LgRZ86cQYsWLTBjxgx89tlnVXoOiIiIrJINldaqLIlQWtvKThQWFsLT0xMFBQXwKPNhFhUVITMzE0FBQXB1NX4qwJ6LOYjdcUntYTA/T1fMHhiCvi00P9xGZG/E+n0jIrIrmZnAl18Cq1YBuhpBOTsDw4crk9c2bUwWnlh05WtlWVRpLlvRt4UfeoX44mTmPeQ9KIK3VDm1wFLuyBIREZEVsfHSWpXFZLaKODpIEN5Qc+1bIiIiogrJZMD33yuT2IqeF2rdGpgxA3j9dausSFAZTGaJiIiILEleHrB8OfD110CZZlFqJBJlZ66331bWkdXSDdXWMZklIiIisgSGltaaOlVZocDOMZklIiIiMheW1qo0JrNEREREpvbggbJ97JIlwO+/6x7brZtyKsGAAYCjowmCsy5MZomIiIhMxU5Ka5kSk1kiIiKiqsTSWlWKHcDshEQiwdatW80dht1JTEzU2R7YmhhzLLzuiMiuyWTA+vVAWJiy2sDmzdoT2TZtlNMOsrOB2bOZyBqAyayNiIqKwuDBg7Wuz8nJQb9+/UwXkIEkEonqp0aNGmjdujUSExPNHValDRs2DFeuXKny9+nWrRskEgnmz59fbl1ERAQkEgnmzJlT5XEQERGUpbXmzgUCA4GRI7XXiJVIgMGDgeRk5ZhRo+yuRqwYmMzaCV9fX7iY+RdEEASUlJRoXZ+QkICcnBycO3cOw4YNw+jRo7F3794qjam4uLhK9+/m5gZvb+8qfY9S/v7+SEhIUFv2xx9/4ODBg/DzYxtlIqIqd/68smRWQAAwa5b2GrFSqfKBrt9/B7ZsAbp2tdsasWJgMquLQgH8+ad5f3TNqzHA01/3ZmVlQSKRYPPmzejevTvc3d3RunVrHDt2TG2bX3/9FV26dIGbmxv8/f0xbdo0PHr0SLU+KSkJYWFhkEql8PX1xfDhw5GXl6dan5ycDIlEgr179yIsLAwuLi44cuSI1hhr1qwJX19fNGzYEO+//z5q166Nffv2qdYXFBTgX//6F7y9veHh4YEXX3wR586dU9vHvHnz4O3tDalUinHjxuG9995Dm6cmzpfewY6Li0PdunXRpEkTAMDt27cxbNgw1KpVC15eXhg0aBCysrLUjuW5555D9erVUbNmTXTq1Ak3btwAAJw7dw7du3eHVCqFh4cH2rVrh9OnTwPQ/NV8fHw8GjZsCGdnZwQHB2P9+vXlPqtVq1ZhyJAhcHd3R+PGjbF9+3at563UgAEDkJ+fj6NHj6qWJSYmonfv3uUS6r/++gsjR45ErVq14O7ujn79+uHq1atqYxITExEQEAB3d3cMGTIE+fn55d5zx44daNeuHVxdXdGgQQPExsbq/A8LEZHNUSiAHTuAHj2UXbjWrNFeI7ZBA+W82Vu3gIULWSNWJExmdcnPV07ENuePhgRCLB988AFmzpyJtLQ0NGnSBG+88YYqEblw4QL69OmDl19+GefPn8d3332HX375BVOmTFFtX1xcjLlz5+LcuXPYunUrMjMzERUVVe59/u///g9xcXFIT09Hq1atKoxLLpfj+++/x71791CtWjUAyru6/fv3R25uLnbt2oUzZ86gbdu26NGjB+7duwcA2LBhAz7++GMsWLAAZ86cQUBAAOLj48vt/8CBA0hPT8f+/fvx008/4fHjx+jevTtq1KiBw4cP45dffkGNGjXQt29fFBcXo6SkBIMHD0bXrl1x/vx5HDt2DP/6178g+ed/0ZGRkahXrx5OnTqFM2fO4L333lPFXdaWLVswffp0vPPOO7h48SImTJiA0aNH49ChQ2rjYmNj8dprr+H8+fOIiIhAZGSk6ji1cXZ2RmRkpNrd2cTERIwZM6bc2KioKJw+fRrbt2/HsWPHIAgCIiIi8OTJEwDAiRMnMGbMGEyaNAlpaWno3r075s2bp7aPvXv3YsSIEZg2bRouXbqE5cuXIzExER9//LHOOImIbMKDB8DSpUBwMPDSS7prxHbrpnzw68oVZXUC1ogVl2BnCgoKBABCQUFBuXV///23cOnSJeHvv/9WLsjLEwTlM4jm+8nL0+u4Ro0aJQwaNEjregDCli1bBEEQhMzMTAGAsGrVKtX63377TQAgpKenC4IgCG+++abwr3/9S20fR44cERwcHP53fso4efKkAEB48OCBIAiCcOjQIQGAsHXr1grjByC4uroK1atXFxwdHQUAQu3atYWrV68KgiAIBw4cEDw8PISioiK17Ro2bCgsX75cEARBeP7554XJkyerre/UqZPQunVr1etRo0YJPj4+gkwmUy1bvXq1EBwcLCgUCtUymUwmuLm5CXv37hXy8/MFAEJycrLG2KVSqZCYmKhxXUJCguDp6al63bFjR2H8+PFqY1599VUhIiJC7Vx8+OGHqtcPHz4UJBKJsHv3bo3vIQiC0LVrV2H69OnCuXPnBKlUKjx8+FBISUkRvL29heLiYqF169bC7NmzBUEQhCtXrggAhKNHj6q2v3v3ruDm5iZ8//33giAIwhtvvCH07dtX7T2GDRumdiydO3cWPvnkE7Ux69evF/z8/NSOpfS6K6vc7xsRkTW4fl0QoqMFwcND97/fzs6CEBUlCKmp5o7YKunK18rinVk79vRd0tI5laXTBM6cOYPExETUqFFD9dOnTx8oFApkZmYCAFJTUzFo0CDUr18fUqkU3bp1AwBkZ2ervU9YWJhe8SxcuBBpaWnYv38/2rRpg4ULF6JRo0aqeB4+fAgvLy+1mDIzM3Ht2jUAQEZGBp577jm1fZZ9DQAtW7aEs7Oz6vWZM2fw+++/QyqVqvZbu3ZtFBUV4dq1a6hduzaioqLQp08fDBw4EIsXL0ZOTo5q++joaIwbNw49e/bE/PnzVfFokp6ejk6dOqkt69SpE9LT09WWPf3ZVK9eHVKpVG0KhzatWrVC48aN8eOPP2LNmjV48803y90lTk9Ph5OTE55//nnVMi8vLwQHB6viSE9PR3h4uNp2ZV+fOXMGH330kdrnMX78eOTk5ODx48cVxkpEZDUEAThyBBg6FGjUCPjiC+01Yr29gTlzlFUJEhJYI9YEWGfWjj2d5JR+Za74Z46uQqHAhAkTMG3atHLbBQQE4NGjR+jduzd69+6NpKQk1KlTB9nZ2ejTp0+5h6qqV6+uVzy+vr5o1KgRGjVqhB9++AGhoaEICwtDSEgIFAoF/Pz8kJycXG67p+ekSspMoBcEodz4svEoFAq0a9cOGzZsKDe2Tp06AJQPp02bNg179uzBd999hw8//BD79+9Hhw4dMGfOHAwfPhw7d+7E7t27MXv2bGzcuBFDhgzReJyaYiy7rGwCKpFIVJ9NRcaMGYOvvvoKly5dwsmTJ8ut13ROysahbczTFAoFYmNj8fLLL5db5+rqqlesREQWTSYDvv9eOc9VW0WCUm3aKB/qev11ViQwMSazunh5KctrmDsGM2jbti1+++031Z3Rsi5cuIC7d+9i/vz58Pf3BwDVQ09iaNSoEYYOHYqYmBhs27YNbdu2RW5uLpycnBAYGKhxm+DgYJw8eRJvvvmmapk+MbVt2xbfffed6sEybUJDQxEaGoqYmBiEh4fjm2++QYcOHQAATZo0QZMmTTBjxgy88cYbSEhI0JjMNmvWDL/88gtGjhypWvbrr7+iWbNmFcapr+HDh2PmzJlo3bo1QkJCyq0PCQlBSUkJTpw4gY4dOwIA8vPzceXKFVUcISEhOH78uNp2ZV+3bdsWGRkZWq8RIiKrlZcHLF8OfP219ooEgLICwaBByiS2SxdWJDATJrO6ODgA/9yZswYFBQVIS0tTW1a7dm0EBAQYvK93330XHTp0wOTJkzF+/HhUr15d9dDU0qVLERAQAGdnZyxduhQTJ07ExYsXMXfuXJGOROmdd95B69atcfr0afTs2RPh4eEYPHgwFixYgODgYPzxxx/YtWsXBg8ejLCwMEydOhXjx49HWFgYOnbsiO+++w7nz59HgwqeFo2MjMRnn32GQYMG4aOPPkK9evWQnZ2NzZs349///jeePHmCFStW4KWXXkLdunWRkZGBK1euYOTIkfj777/x73//G6+88gqCgoJw69YtnDp1CkOHDtX4Xv/+97/x2muvqR5e27FjBzZv3oyff/5ZtPNWq1Yt5OTkaH0IrXHjxhg0aBDGjx+P5cuXQyqV4r333sOzzz6LQYMGAQCmTZuGjh074tNPP8XgwYOxb98+7NmzR20/s2bNwoABA+Dv749XX30VDg4OOH/+PC5cuFDuYTEiIqtw/jyweDGwYYP2igSAsrTW2LHA1KmsSGABOGfWhiQnJ6vuHpb+zJo1y6h9tWrVCikpKbh69So6d+6M0NBQ/Oc//1HNra1Tpw4SExPxww8/ICQkBPPnz8fnn38u5uGgZcuW6NmzJ2bNmgWJRIJdu3ahS5cuGDNmDJo0aYLXX38dWVlZ8PmnS0pkZCRiYmIwc+ZMtG3bVlVdoaKvvN3d3XH48GEEBATg5ZdfRrNmzTBmzBj8/fff8PDwgLu7Oy5fvoyhQ4eiSZMm+Ne//oUpU6ZgwoQJcHR0RH5+PkaOHIkmTZrgtddeQ79+/RAbG6vxvQYPHozFixfjs88+Q/PmzbF8+XIkJCSo5huLpWbNmjqndyQkJKBdu3YYMGAAwsPDIQgCdu3apUqAO3TogFWrVmHp0qVo06YN9u3bhw8//FBtH3369MFPP/2E/fv3o3379ujQoQO++OIL1K9fX9RjISKqUiytZfUkgj6T42xIYWEhPD09UVBQUO4r5aKiImRmZiIoKIhz/mxEr1694OvrW66WK5kff9+IyKwePFC2j12yRNm8QJdu3ZRTCQYMABwdTRAc6crXyuI0A7IZjx8/xrJly9CnTx84Ojri22+/xc8//4z9+/ebOzQiIrIUmZnAl18Cq1Zpr0gAAM7OwPDhyrqwrEhg0ZjMks0onYowb948yGQyBAcHY9OmTejZs6e5QyMiInMSBOCXX5RTBLZu1d1d09sbmDQJmDgR+GcaG1k2JrNkM9zc3ER9kIqIiKwcS2vZBSazREREZFtYWsuuMJklIiIi28DSWnaJyawGdlbggcgs+HtGRKJQKICdO5VTCQ4e1D22QQNg2jRg9GiggifkyXowmX1KaY3Nx48fw83NzczRENm2x48fAyjfupeISC+lpbUWLwauXdM9lqW1bBqT2ac4OjqiZs2ayPunha27u7uqVz0RiUMQBDx+/Bh5eXmoWbMmHPkPCxEZgqW1qAwms2X4+voCgCqhJaKqUbNmTdXvGxGRTiytRTowmS1DIpHAz88P3t7eePLkibnDIbJJ1apV4x1ZIqoYS2uRHpjMauHo6Mh/bImIiMyBpbXIAExmiYiIyDKwtBYZgcksERERmY9c/r/SWocO6R7L0lqkAZNZIiIiMj2W1iKRMJklIiIi02FpLRIZk1kiIiKqWiytRVWIySwRERFVDZkM+O475VQCltaiKsJkloiIiMSVlwcsW6YsrXXnjvZxLK1FInAwdwCl4uLiIJFI8Pbbb2sdk5ycDIlEUu7n8uXLpguUiIiINDt3DhgzBggIAGbP1p7ISqXKBPb334EtW4CuXZnIktEs4s7sqVOnsGLFCrRq1Uqv8RkZGfB4qiRHnTp1qio0IiIi0oWltcjMzJ7MPnz4EJGRkVi5ciXmzZun1zbe3t6oWbNm1QZGRERE2rG0FlkIs08zmDx5Mvr374+ePXvqvU1oaCj8/PzQo0cPHKrgf4EymQyFhYVqP0RERGSkzEzgnXeAevWUd1m1JbLOzkBUFJCaqrxjO2gQE1mqEma9M7tx40acPXsWp06d0mu8n58fVqxYgXbt2kEmk2H9+vXo0aMHkpOT0aVLF43bxMXFITY2VsywiYiI7IsgAEeOKKcSbNvG0lpkUSSCIAjmeOObN28iLCwM+/btQ+vWrQEA3bp1Q5s2bbBo0SK99zNw4EBIJBJs375d43qZTAbZU/2dCwsL4e/vj4KCArV5t0RERFRGaWmtRYuUd1h1YWktElFhYSE8PT31ytfMdmf2zJkzyMvLQ7t27VTL5HI5Dh8+jC+//BIymQyOenwd0aFDByQlJWld7+LiAhf+UhEREemPpbXIipgtme3RowcuXLigtmz06NFo2rQp3n33Xb0SWQBITU2Fn59fVYRIRERkX86dUz7Q9c03yruy2kilwNixwNSpygoFRGZktmRWKpWiRYsWasuqV68OLy8v1fKYmBjcvn0b69atAwAsWrQIgYGBaN68OYqLi5GUlIRNmzZh06ZNJo+fiIjIJrC0Flk5s5fm0iUnJwfZ2dmq18XFxZg5cyZu374NNzc3NG/eHDt37kRERIQZoyQiIrJCLK1FNsJsD4CZiyETiomIyDhyhYCTmfeQ96AI3lJXPBdUG44OnE9pETIzgaVLgdWrAV3lKp2dgeHDgenTlQ93EZmQVTwARkREtmnPxRzE7riEnIIi1TI/T1fMHhiCvi34jINZsLQW2TAms0REJJo9F3PwVtJZlP3KL7egCG8lnUX8iLZMaE2JpbXIDjCZJSIiUcgVAmJ3XCqXyAKAAEACIHbHJfQK8eWUg6rG0lpkR5jMEhGRKE5m3lObWlCWACCnoAgnM+8hvKGX6QKzJyytRXaIySwREYki74H2RNaYcaQnltYiO8dkloiIROEtdRV1HFXgwQMgIQFYsoSltciuMZklIiJRPBdUG36ersgtKNI4b1YCwNdTWaaLKoGltYjUOJg7ACIisg2ODhLMHhgCQJm4Pq309eyBIXz4yxiCABw+DLz8MtCoEbBwofZE1tsbmDMHyM5W3rllIks2jsksERGJpm8LP8SPaAtfT/WpBL6erizLZQyZDFi3DmjXDujaFdiyRXuN2DZtgLVrlUns7NmsEUt2g9MMiIhIVH1b+KFXiC87gFUGS2sR6Y3JLBERic7RQcLyW8ZgaS0igzGZJSIiMieW1iKqFCazRERE5sDSWkSiYDJLRERkSiytRSQqJrNERERVTRCAI0eUUwm2bdNekQBQViGYNAmYMIEVCYj0wGSWiIioqshkwHffKZPY1FTdY9u0AWbMAIYNA1xcTBEdkU1gMktERCQ2ltYiMhkms0RERGJhaS0ik2MyS0REVBksrUVkVkxmiYismFwhsNOWuRhSWqt7d2VVApbWIhIdk1kiIiu152IOYndcQk5BkWqZn6crZg8MQd8WfmaMzMYZUlorMlKZxLZubbr4iOwMk1kiIiu052IO3ko6C6HM8tyCIryVdBbxI9oyoRUTS2sRWSwms0REVkauEBC741K5RBYABAASALE7LqFXiC+nHFQWS2sRWTwms0REVuZk5j21qQVlCQByCopwMvMewht6mS4wW8LSWkRWg8ksEZGVyXugPZE1Zhw9pbS01oYNQHGx9nEsrUVkMZjMEhFZGW+pq6jj7B5LaxFZNSazRERW5rmg2vDzdEVuQZHGebMSAL6eyjJdpIOhpbXefhvo35+ltYgsjIO5AyAiIsM4Okgwe2AIAGXi+rTS17MHhvDhL20yM4HoaKBePWXZLG2JrLOz8g5sWhpw8CDw0ktMZIksEJNZIiIr1LeFH+JHtIWvp/pUAl9PV5bl0kQQgMOHgZdfBho1AhYu1F4j1scHiI0FsrOBNWtYI5bIwnGaARHZNFvukNW3hR96hfja7PGJgqW1iGwek1kisln20CHL0UHC8luasLQWkd1gMktENokdsuyUIaW1xo0DpkxhaS0iK8dklohsDjtk2RlDS2tNnw5ERbG0FpGNYDJLRDaHHbLsBEtrERGYzBKRDWKHLBuXmQksXQqsXq29IgGgLK0VGam8E8uKBEQ2i8ksEdkcdsiyQYIAHDminEqwbRugUGgf6+MDTJoETJig/DMR2TQms0Rkc9ghy4awtBYRVYBNE4jI5rBDlg3IywM++gioXx8YNUp7IiuRAEOGACkpwNmzwMiRTGSJ7AyTWSKySeyQZaXOnQPGjAH8/YHZs7XXiJVKlXdhf/8d2LyZNWKJ7BinGRCRzWKHLCvB0lpEVAkGJ7OPHj3C/PnzceDAAeTl5UFRZhL+9evXRQuOiKiy2CHLghlQWqugwwu4PGwMFBH98VyjOvwPCRGpGJzMjhs3DikpKXjzzTfh5+cHiUhf68TFxeH999/H9OnTsWjRIq3jUlJSEB0djd9++w1169bF//3f/2HixImixEBERCZgQGmtWxFDEOP/Io64PwvkAlhzyuZaEhNR5RiczO7evRs7d+5Ep06dRAvi1KlTWLFiBVq1aqVzXGZmJiIiIjB+/HgkJSXh6NGjmDRpEurUqYOhQ4eKFg8REYnMiNJaB7sMxtg9N9mSmIh0MvgBsFq1aqF2bfHK2Tx8+BCRkZFYuXIlatWqpXPssmXLEBAQgEWLFqFZs2YYN24cxowZg88//1y0eIiISEQyGbBuHdCuHdC1K7Bli/ZEtk0bYO1a4MYNyD/8Dz449qfWlsSAsiWxXKFpBBHZE4OT2blz52LWrFl4/PixKAFMnjwZ/fv3R8+ePSsce+zYMfTu3VttWZ8+fXD69Gk8efJE4zYymQyFhYVqP0REVMUqWVrLkJbERGTf9JpmEBoaqjY39vfff4ePjw8CAwNRrVo1tbFnz57V+803btyIs2fP4tSpU3qNz83NhU+Zbi4+Pj4oKSnB3bt34edX/uumuLg4xMbG6h0TERFVwrlzwOLFwIYNQHGx9nFSKTBuHDBlirJCQRlsSUxE+tIrmR08eLDob3zz5k1Mnz4d+/btg6ur/i0lyz5wJgiCxuWlYmJiEB0drXpdWFgIf39/IyImIiKNqqC0FlsSE5G+9EpmZ8+eLfobnzlzBnl5eWjXrp1qmVwux+HDh/Hll19CJpPB0dFRbRtfX1/k5uaqLcvLy4OTkxO8vDSX3nFxcYELu8EQEYnPgNJa6N4dePttoH9/oMzf7ZqwJTER6cvgObMNGjRAfn5+ueX3799HAw1fFWnTo0cPXLhwAWlpaaqfsLAwREZGIi0trVwiCwDh4eHYv3+/2rJ9+/YhLCys3HQHIiKqIpmZQHQ0UK+e8i6rtkTW2RkYPRpISwMOHgReekmvRBZgS2Ii0p/BpbmysrIgl8vLLZfJZLh165be+5FKpWjRooXasurVq8PLy0u1PCYmBrdv38a6desAABMnTsSXX36J6OhojB8/HseOHcPq1avx7bffGnoYRERkCCNKa2HCBOWfjVTakjh2xyW1h8F8WWeWiJ6idzK7fft21Z/37t0LT09P1Wu5XI4DBw4gKChI1OBycnKQnZ2teh0UFIRdu3ZhxowZ+Oqrr1C3bl0sWbKENWaJyOoUlyiw/lgWbtx7jPq13fFmeCCcnQz+sqzqyWTAd98pk1htFQlKtWkDzJgBDBsGiDS9iy2JiagiEqH0CaoKODgo/5KVSCQou0m1atUQGBiI//73vxgwYID4UYqosLAQnp6eKCgogAf7ehORGcTtuoSVRzLxdIlUBwkwvnMQYiJCzBfY0/LygGXLgK+/Bu7c0T5OIgEGD1bOh+3cWfmaiKiSDMnX9L4zq/jnK6WgoCCcOnUKzzzzTOWiJCKyQ3G7LmH54cxyyxUCVMvNmtCKVFqLiMhUDJ4zm5lZ/i9hIiKqWHGJAiuP6P47dOWRTLzTu6lppxxUQWktIiJTMTiZXbJkicblEokErq6uaNSoEbp06aKxGgERkT1bfywLFXVfVQjKcWM7m+BuZxWW1iIiMhWDk9mFCxfizz//xOPHj1GrVi0IgoD79+/D3d0dNWrUQF5eHho0aIBDhw6xOQER0VNu3NOvDbi+44yWmQksXQqsXg3oavHt7AxERirvxLZuXbUxEREZyeDvsT755BO0b98eV69eRX5+Pu7du4crV67g+eefx+LFi5GdnQ1fX1/MmDGjKuIlIrJa9Wu7izrOIIIAHD4MvPwy0KgRsHCh9kTWxweIjQWys4E1a5jIEpFF07uaQamGDRti06ZNaNOmjdry1NRUDB06FNevX8evv/6KoUOHIicnR8xYRcFqBkRkLsUlCjT9z26dUw0cJMDluf3EmzNr5tJaRETGqJJqBqVycnJQUlJSbnlJSYmq1WzdunXx4MEDQ3dNRGTTnJ0cML5zkMZqBqXGdw4SJ5FlaS0ishMG/43ZvXt3TJgwAalP/Q8/NTUVb731Fl588UUAwIULF0RvoEBEZAtiIkIwoUsQytb8d5AAE7qIUGf23DlgzBjA3x+YPVt7IiuVKu/C/v47sHkz0KULE1kiskoGTzPIzc3Fm2++iQMHDqBatWoAlHdle/TogfXr18PHxweHDh3CkydP0Lt37yoJujI4zYDI8v1dLMcnuy4hK/8xAr3c8X5ECNycbesJelE7gBlSWqthQ2DatCovrSVXCOzaRURGMyRfMziZLXX58mVcuXIFgiCgadOmCA4ONipYU2MyS2TZxq87hf2X8sot7xXijZUj25shIgtmoaW19lzMQeyOS8gpKFIt8/N0xeyBIejbwq9K35uIbINJkllrxWSWyHJpS2RLMaH9hwWX1tpzMQdvJZ1F2X9YSu/Jxo9oy4SWiCpUpQ+AyeVyJCYm4sCBA8jLy1O1uS118OBBQ3dJRIS/i+U6E1kA2H8pD38Xy21uyoFeBAE4ckQ5lWDbNqDM371qfHyASZOACROUfzYRuUJA7I5L5RJZABCgTGhjd1xCrxBfTjkgItEYnMxOnz4diYmJ6N+/P1q0aAEJHxggIhF8suuS3uPmDm5ZxdFYECsqrXUy857a1IKyBAA5BUU4mXkP4Q29TBcYEdk0g5PZjRs34vvvv0dERERVxENEdiorX7+uV/qOs3pWWFor74H2RNaYcURE+jA4mXV2dkajRo2qIhYismOBXu44clW/cTbt3Dlg8WJgwwaguFj7OKkUGDcOmDIFaNDAdPHp4C11FXUcEZE+DK4D884772Dx4sWws+fGiKiKva9nfVV9x1kVuRzYvh148UXlVIGEBO2JbMOGymT31i3giy8sJpEFgOeCasPP0xXa7g1LoKxq8FxQbVOGRUQ2zuA7s7/88gsOHTqE3bt3o3nz5qpas6U2b94sWnBEZD/cnB3RK8S7wmoGNvXwl4WW1jKWo4MEsweG4K2ks5AAag+ClSa4sweG8OEvIhKVwclszZo1MWTIkKqIhYjs3MqR7e2jzqwFl9aqrL4t/BA/om25OrO+rDNLRFWEdWaJyOLYZAcwKyitJSZ2ACOiyqjSOrOAsn1tcnIyrl27huHDh0MqleKPP/6Ah4cHatSoYVTQRESlnJ0cENGyrioRMrrNq8iMStBkMmDjRuU8VwsvrSUmRwcJy28RkUkYnMzeuHEDffv2RXZ2NmQyGXr16gWpVIpPP/0URUVFWLZsWVXESUR2wlJboRoc1507ytJa8fE6S2spIMG+Jh2wreurGDR1GPq2rFsV4RMR2SyDb3dMnz4dYWFh+Ouvv+Dm5qZaPmTIEBw4cEDU4IjIvpS2Qi1beD+3oAhvJZ3Fnos5lh/XuXPAmDFAQAAwZ47WRPaBsxtWhQ1C1wkrMXHIB9hTuwne2pBqtmMkIrJWRlUzOHr0KJydndWW169fH7dv3xYtMCKyL5baClWfuOZuu4Bev5+E45LFwKFDOvd3q3ZdrAodgB9b9sRDl//VzGW7VyIi4xiczCoUCsjl8nLLb926BalUKkpQRGR/LLUVqq64qsse49ULPyPqzA443q/gjmr37rg8bAwirnlA4aD5YTa2eyUiMpzByWyvXr2waNEirFixAgAgkUjw8OFDzJ49my1uicholtoKVdP71bufi6gzO/Da+f3wKNbRXrdMaa2MtNtQZKYZ9Z5ERKSZwcnswoUL0b17d4SEhKCoqAjDhw/H1atX8cwzz+Dbb7+tihiJyA5YaitU1fsJAp6/eRGjz2xHr6sn4CgYXlrLUo+RiMiaGZzM1q1bF2lpafj2229x9uxZKBQKjB07FpGRkWoPhBERGaK0FWpuQZHG+akSKAvvm7oV6nN1q2PstcN4+fCPaJ53Xffg0FBlly4tpbUs9RiJiKwZmyYQkcUorRoAaG6FGj+irenKcxlQWiuvR1/4znoP6NwZkOh+cMuijpGIyEIZkq/plcxu375d7zd/6aWX9B5rDkxmiSyb2evMnjunbHCwYQNQXKx12ANnN/zUPgJ1P5iJrv06GPQWZj9GIiILJ3oy6+CgXzlaiUSisdKBJWEyS4awh5acYh2jmOfK5OddLgd27lS2mq2gtFZRQBCuDovC38NHol2r+tZzjEREVkT0drYKXT3EiWyUPdw9E+sYxT5XJmuF+uABkJAALFkCXLume2z37sDbb8O1f3+0dNRcWssQbPdKRCQOzpkl0qB0XmPZXw5bmtco1jFa5bm6fh1YuhRYswYoLNQ+rkxpLSIiMg1D8jWD29kS2bqKOj4Byi5NcoX1/j9QrGO0qnMlCEBKCjBkCNCokXJKgbZE1scHiI0FsrOVCS8TWSIii8VklqgMQzpRWSuxjtEqzpVMBqxdC7RtC3TrBmzdqkxsNQkNVY69cQOYNUutRiwREVkmg+vMEtk6S+1EJSaxjtGiz5WepbUgkQCDByvrw+pRWouIiCwLk1miMuyhS5NYx2iR50rP0lqQSoFx44ApU4AGDUwXHxERiUqvZLZQ1wMSZfChKrJ29tClSaxjtJhzZUBpLTRsCEybBkRFAfz7iojI6uk1Z7ZmzZqoVauWzp/SMUTWztFBgtkDQwD874n8UqWvZw8MseqaoGIdo9nPVWGh8i5skybAoEG6E9nu3YFt24CMDGUyy0SWiMgm6FWaKyUlRe8ddu3atVIBVTWW5iJ9sc6s+erMVqi0tNbq1cpasdqwtBYRkVUSvQOYLWEyS4awhy5NYh3j38VyfLLrErLyHyPQyx3vR4TAzdm45gLFJQqsP5aFG/ceo35td7wZHghnRwlw+LByKsG2bdorEgDKKgSTJgETJ0L+TB2b/gzt4RolIvtjkmT28ePHyM7ORnGZByxatWql9z7i4+MRHx+PrKwsAEDz5s0xa9Ys9OvXT+P45ORkdO/evdzy9PR0NG3aVK/3ZDJLJL64XZew8kgmni4n6yABxncOQkxESKX25VzyBC9dPoyZ6bvhe/2y7o1DQ5VVCYYNA1xcbP7uuq0fHxHZL9Hb2T7tzz//xOjRo7F7926N6+Vyud77qlevHubPn49GjRoBANauXYtBgwYhNTUVzZs317pdRkaG2oHVqVNH7/ckInHF7bqE5Yczyy1XCFAt1zehfXpfzzz6C5GpuzEibRfqPLqvfSMtpbW0dSbLLSjCW0lnLbMzmQFs/fiIiPRlcNOEt99+G3/99ReOHz8ONzc37NmzB2vXrkXjxo2xfft2g/Y1cOBAREREoEmTJmjSpAk+/vhj1KhRA8ePH9e5nbe3N3x9fVU/jiL0SSciwxWXKLDySPlE9mkrj2SiuESh975C7lzHZzsX4Wj8aMw4+o32RFYqBWbMAH7/Hdi8GejSRZXIWlVnMiPY+vERERnC4DuzBw8exLZt29C+fXs4ODigfv366NWrFzw8PBAXF4f+/fsbFYhcLscPP/yAR48eITw8XOfY0NBQFBUVISQkBB9++KHGqQelZDIZZDKZ6rUhZcaISLf1x7JQUb6kEJTjxnbWUctVLkfKp6uw4ZuvEJ59Qef+Cp8NgMf/vaOztJYhncnCG3rpPgALZOvHR0RkCIOT2UePHsHb2xsAULt2bfz5559o0qQJWrZsibNnzxocwIULFxAeHo6ioiLUqFEDW7ZsQUiI5q8k/fz8sGLFCrRr1w4ymQzr169Hjx49kJycjC5dumjcJi4uDrGxsQbHRUQVu3HvceXGFRYCCQnAkiXodf26zn38GtAKa8IG4dnIoYh9WXdlAovuTCYCWz8+IiJDGJzMBgcHIyMjA4GBgWjTpg2WL1+OwMBALFu2DH5+hs/PCg4ORlpaGu7fv49NmzZh1KhRSElJ0ZjQBgcHIzg4WPU6PDwcN2/exOeff641mY2JiUF0dLTqdWFhIfz9/Q2Ok4jKq1/b3bhxepbWkjlWw9aQbkgIewmXvYMAAP+pI63w/SyyM5mIbP34iIgMYXAy+/bbbyMnJwcAMHv2bPTp0wcbNmyAs7MzEhMTDQ7A2dlZ9QBYWFgYTp06hcWLF2P58uV6bd+hQwckJSVpXe/i4gIXFxeD4yKiir0ZHoiPd6XrnGrgIFGOgyDoXVrrz+o1sT60Pza06Yf86jXL76sCFtOZrIrY+vERERnC4GQ2MjJS9efQ0FBkZWXh8uXLCAgIwDPPPFPpgARBUJvjWpHU1FSj7ggTUeU5OzlgfOcgjdUMSk3s8CycN6xXJrFpaTr3l9uwGT5t2hc/Ne2CYqdq5daP7xwEZ6eKn1st7Uz2VtJZSAC1hM8WurjZ+vERERnC4GoGH330ER4//t/8N3d3d7Rt2xbVq1fHRx99ZNC+3n//fRw5cgRZWVm4cOECPvjgAyQnJ6sS5piYGIwcOVI1ftGiRdi6dSuuXr2K3377DTExMdi0aROmTJli6GEQkUhiIkIwoUsQyuZN3o//wvpbu/F/E/oqH9bSlshKJMCQIUBKCnyv/oY6k8ahpJp6IusgASZ0Maxmbd8Wfogf0Ra+nupftft6utpE2SpbPz4iIn0Z3DTB0dEROTk5qofASuXn58Pb29ugOrNjx47FgQMHkJOTA09PT7Rq1QrvvvsuevXqBQCIiopCVlYWkpOTAQCffvopVqxYgdu3b8PNzQ3NmzdHTEwMIiIi9H5PNk0gQ1hqdyWNHbL0uGNZlR4WlWDGd6lwungeI45vRceT+yAp01RFjVQKjBsHTJkCNFCvdFDl3cTMfK7EjMtSr1Eiosqo0g5gDg4OuHPnTrlGBQcPHsSwYcPw559/Gh6xCTGZJX1ZanclMbttiWX+jgu4nvg9Rp/eVmFpLTRsCEybprW0lpjnnZ8hEZF1qpJktlatWpBIJKqdSiT/+5+/XC7Hw4cPMXHiRHz11VeVi76KMZklfWjrrlR61Zvra1xt3bZKGfpVfKUVFmL/v+PQ5PtE1L+fq3vsiy8qu3RFRABaGp2Ied75GRIRWa8qSWbXrl0LQRAwZswYLFq0CJ6enqp1zs7OCAwMrLDZgSVgMksVkSsEvLDgoNai9KVPiv/y7osm/Tq3uESBpv/ZXWHlgMtz+1X91+j/lNYSVq+GpILSWtuad8OQhE9RrW0bnbsU87zzMyQism6G5Gt6VzMYNWoUACAoKAidOnWCk5PBhRCIrIKldlcSrduWsTSU1tKWBpYtrfXgkQfGVrB7Mc87P0MiIvthcEbatWtXXLt2DQkJCbh27RoWL14Mb29v7NmzB/7+/mjevHlVxElkMpbaXanS3baMJZMBGzfqVVrrok9DrAl7qVxpLX1iEvO88zMkIrIfBn+PlZKSgpYtW+LEiRPYvHkzHj58CAA4f/48Zs+eLXqARKZmqd2VjO62Zaw7d4DYWCAgQGdpLQUk2NMkHK8Nn48BoxZhc4se5WrE6hOTmOednyERkf0wOJl97733MG/ePOzfvx/Ozs6q5d27d8exY8dEDY7IHEq7K2n7Cl0C5RPxpu6u9GZ4YLlarmXp2yFLp7Q0YPRoZRI7Zw6Ql6d5nIcH5NPfRveJKzFxyAc46d9CWTPWyJjEPO92/xkSEdkRg5PZCxcuYMiQIeWW16lTB/n5+aIERWROpd2VAJRLhszZXam025Yu+nbIKkcuV86D7d4dCA0FEhMBbTViGzYEliwBbt2C46KF6DtQ94OfhnbtAip/3u3yMyQislMG/41Zs2ZN5OTklFuempqKZ599VpSgiMzNUrsraeu2ZUyHLABAYSGweDHQpAkweDDwT4MSjV58Edi+HcjIAKZOVTY9EDkmMc+73XyGRER2zuCmCf/3f/+HY8eO4YcffkCTJk1w9uxZ3LlzByNHjsTIkSMtft4sS3ORISy1u1Klu0f9U1oLq1cDOkprwcUFiIwEpk8HWrWq2pieIuZ5t9nPkIjIhlVpB7AnT54gKioKGzduhCAIcHJyglwux/Dhw5GYmAhHLcXQLQWTWbJbGkpraeXjA0yaBEycCJRpXa2NPSSgRERkGlWazJa6du0aUlNToVAoEBoaisaNGxsVrKkxmSW7Y0BpLYSGKrt0DRumvCurJ3toQUtERKZjkmQWAEo3lWh4gtlSMZklu3HnDrBsGfD119orEgCAg4Nyvuz06UDnzhorEuhiDy1oiYjItAzJ14yaoLV69Wq0aNECrq6ucHV1RYsWLbBq1SqjgiUikRlQWgvR0cDvvwObNgFduhicyMoVAmJ3XCqXfAJQLYvdcQnyitpeibwvIiKyHwZ3APvPf/6DhQsXYurUqQgPV5bkOXbsGGbMmIGsrCzMmzdP9CCJqAJyOfDTT8qpBLoqEgDK0lrTpysbIfxTkcBY9tCCloiILJvByWx8fDxWrlyJN954Q7XspZdeQqtWrTB16lQms0SmVFgIJCQo675ev6577IsvKufDRkQAIj2oaQ8taImIyLIZnMzK5XKEhYWVW96uXTuUlJSIEhQRVaAKSmsZwx5a0BIRkWUzeM7siBEjEB8fX275ihUrEBkZKUpQRKSBIAApKcCQIUCjRsopBdoSWR8fIDYWyM5WJrxVkMgC9tGCloiILJvBd2YB5QNg+/btQ4cOHQAAx48fx82bNzFy5EhER0erxn3xxRfiRElkz0xQWstYpW1j30o6Cwmg9vCWsS1oxdgXERHZD4NLc3Xv3l2/HUskOHjwoFFBVSWW5iKrYWhprbffBl54weCKBGJgnVkiIhKTyerMWiMms5bLErs+iRmT3u1L09KAxYuBb74Biou179DDAxg3DpgyBQgKMiomMVlqO1siIrI+TGZ1YDJrmSzxbpyYMcXtuoSVRzLxdIlUBwkwvnMQYiJCzFZaSyyW+PkREZH1YjKrA5NZy2OJXZ/EjClu1yUsP5ypcV0N2WMsfHwWvX7+3iyltcRgiZ8fERFZN0PyNaMeACMSS0VdnyRQdn3qFeJrsq+ZxYypuESBlUfKJ7L+93MRdWYHXju/D9Liv7XvoIpLa1WWJX5+RERkX5jMkllZYtcnMWNafyzrf1MLBAHP37yIMae3odfVE3DQmAL+w8cHmDwZmDAB8PY2/CBMxBI/PyIisi9MZsmsLLHrk5gx3bj3GM4lTzAw/TDGnN6G5nkVTCUIDQVmzABee80kpbUqyxI/PyIisi9MZsmsLLHrk2gx3bmDwdtWYep3a1Hn8X2tw+QSB2R37oWgeR+YrbSWsSzx8yMiIvvCZJbMqrTrU25BkcYv3SUAfE3c9anSMT1VWqutjtJahc7u+K51b6xvNwA/Lx0NGFnGypws8fMjIiL7Yn3/epJNKe36BKBcG1NzdX0yKia5HNi2DejeXTlVIDFRa43YrJp+mN1zAsInJeLjF8eh38Bwo+uxmpslfn5ERGRfWJqLLIIl1inVK6bCQiAhAViypMLSWr/Wb4XVYYNwqEEYFA6O6nVmrZwlfn5ERGS9WGdWByazlssSuz5p7Wp1/TqwdCmwejXw4IH2HTxVWqs4pIVFdsgSq3OXJX5+RERknZjM6sBklvRV7m6jICDi3hXMvr4fPof2Arp+daqwtJZJO5MRERGZAZsmEFXS012tLKm0lrZuW7kFRXgr6awonckUAlTLmdASEZGlYzJLVEZpVyuvR38hMnU3RqTu0llaCw4OwODBylazVVhayxSdyZ628kgm3und1GofTiMiIvvAZJaojIs7UxD9zXy8lJ4MF3mJ1nElNaRw+td4YMoUICioyuOqss5kWigE5bixnRsYES0REZFpMJklApSltX76CVi0CK2Tk9Fax9Csmn5ICHsJz82Zgf6dgk0WotidyfSh7zgiIiJzYTJL9s2A0lpH67fCmqdKa/X1fcZEQSqJ2W2rfm13vfal7zgiIiJzYTJL9knP0loyx2rYGtINCWEv4bK3ciqBBMrqAabuaiVmt603wwPx8a50nVMNHCTKcURERJaMySzZD0EADh8GFi1SduvSUVpL5lUHXzXrjW/a9MPd6jVVy83Z1aq029ZbSWchAdQSWkPjcnZywPjOQRqrGZQa3zmID38REZHF479UZPtkMmDtWqBtW6BbN2DrVu2JbGgosG4dXG7fREj8Z6hW11dtta+nq0Hlr8TWt4Uf4ke0ha+n+lQCY+KKiQjBhC5BKJv7OkiACV1YZ5aIiKwDmyaQ7bpzB1i2DPj6ayAvT/s4HaW1LLWrlSV2ACMiIhKLIfmaWf/Fio+PR6tWreDh4QEPDw+Eh4dj9+7dOrdJSUlBu3bt4OrqigYNGmDZsmUmipY0kSsEHLuWj21pt3HsWj7kFdV7MsW+0tKA0aOBgABgzhztiayHBxAdDfz+O7BpE9C5c7kasY4OEoQ39MKgNs8ivKGXRSSygPJcXfqjAGdu/IVLfxRU6rw7OkgQUtcT7erXQkhdT6OPUcxrgYiISF9mnTNbr149zJ8/H40aNQIArF27FoMGDUJqaiqaN29ebnxmZiYiIiIwfvx4JCUl4ejRo5g0aRLq1KmDoUOHmjp8uydmW9VK7+up0lpITtY9tmFDYPp0ICoKkEoNitMSaGpB+/GudKNa0Ir1GYp5LRARERnC4qYZ1K5dG5999hnGjh1bbt27776L7du3Iz09XbVs4sSJOHfuHI4dO6bX/jnNQBza2qqW3tMzZP5mpfZlQGktvPiicipBRATg6KhXbJZGWwvaUobMdRXrMxTzWiAiIgKsaJrB0+RyOTZu3IhHjx4hPDxc45hjx46hd+/easv69OmD06dP48mTJ6YIk1BxW1VA2VZVn6+Zjd7X9evAjBlAvXrKBFVbIuviAowZA5w7Bxw4AAwcaLWJrL4taItLFBXuS6zPUMxrgYiIyBhmT2YvXLiAGjVqwMXFBRMnTsSWLVsQEqL5zlJubi58fHzUlvn4+KCkpAR3797VuI1MJkNhYaHaD1WOIW1VRd2XIAApKcCQIUCjRsopBdpqxPr4AB99BGRnK2vJtmpVYSyWzpAWtBUR6zMU81ogIiIyhtnrzAYHByMtLQ3379/Hpk2bMGrUKKSkpGhNaCVlHtApnSVRdnmpuLg4xMbGihu0nROzrao+Y5xLnsBlw3pgy1rlw126hIYq79i+9pryrqwNEbMFrVifoZjXAhERkTHMnsw6OzurHgALCwvDqVOnsHjxYixfvrzcWF9fX+Tm5qoty8vLg5OTE7y8vDTuPyYmBtHR0arXhYWF8Pf3F/EI7I+YbVV1jXnm0V+ITN2NEam7UOfxfe070VFay5aI2YJWrM9QzGuBiIjIGGZPZssSBAEymUzjuvDwcOzYsUNt2b59+xAWFoZq1app3MbFxQUuNnaHztzEbKuqaV8hd65j9OnteCk9GS7yEu0be3gA48YBU6YAQUHGHIpVEbMFrVifoZjXAhERkTHMOmf2/fffx5EjR5CVlYULFy7ggw8+QHJyMiIjIwEo76qOHDlSNX7ixIm4ceMGoqOjkZ6ejjVr1mD16tWYOXOmuQ7BLpW2VQX+98R6KUPbqpbuy0EhR++rx/HttzHYlTgNr178WXsi27ChsnrBrVvAf/9rF4ks8L8WtLro24JWrM9QzGuBiIjIGGZNZu/cuYM333wTwcHB6NGjB06cOIE9e/agV69eAICcnBxkZ2erxgcFBWHXrl1ITk5GmzZtMHfuXCxZsoQ1Zs1AtLaqhYXoe+B7XPh2GlZsnofw7Avax774IrB9O5CRAUydapU1YitLzBa0Yn2GYrbYJSIiMpTF1ZmtaqwzKy6j26pevw4sXaqsNKCtIgGgfIgrMlLZ5MAGKhKIRcwWtGK1xrXU1r9ERGR9DMnXmMyS6QgCcPiwsqTWtm3K11r8Wb0mtoQPRqMPo/Fi15ami5GIiIjMzpB8zeIeACMbJJMBGzcqk9gKSmtd9GmI1WGDsLNpZzxxqgbszka81zP8qpqIiIg0YjJLVefOHWDZMuDrr4G8PK3D5BIH7GvcAWvCXsKpes3VSmtJoOwg1SvEl19ZExERUTlMZkl8aWnA4sXAN98AxcXax3l44I9XIvFatXa4VdNX45CnO0iFN9RcS5iIiIjsF5NZEodcDvz0k3IqQXKy7rENGyof6IqKwqlrhbi1Ma3C3bODFBEREWnCZJYqp7AQSEhQ1n29fl332BdfVHbpiogAHB0BAN5SHXdun8IOUkRERKQJk1kyjkiltdhBioiIiCrDrE0TyMoIgnIKweDBQKNGyikF2hJZX1/go4+A7GxlwqulRiw7SBEREVFlMJmlislkwNq1QNu2QPfuumvEhoYC69YBWVnAf/4DeHtXuHt2kCIiIiJjcZoBaadnaS04OCjv1r79NvDCC2qltfTVt4UfeoX4soMUERERGYTJLJVnQGktjBsHTJkCBAVV+m0dHSQsv0VEREQGYTJLSkaW1oJUaoLgiIiIiDRjMmvvCguBNWuUlQmMKK1FREREZE5MZu3VtWvKBHbNmkqV1iIiIiIyJyaz9kQQgJQU5VSC7du1VyQAlKW1Jk0CJkzQqyIBERERkTkwmbUHMhmwcaMyiU1L0z02NBSYMQN47TXlXVkiIiIiC8Zk1paZsLQWERERkTkwmbVFZiqtRURERGRqTGZthVwO7NihTGJZWouIiIjsBJNZa1daWmvJEiAzU/dYltYiIiIiG8Nk1lqxtBYRERERk1mrwtJaRERERGqYzFoDltYiIiIi0ojJrCVjaS0iIiIinZjMWqK0NOVd2G+/ZWktIiIiIh2YzFqK0tJaixYp58XqwtJaRERERACYzJqfMaW1+vdXTi0gIiIisnNMZs2FpbWIiIiIKo3JrCmxtBYRERGRqJjMmgJLaxERERFVCSazVenOHSA+XvnD0lpEREREomMyW1X+/hsIDgYKCrSPYWktIiIiokrhI/FVxc0NePVVzesaNVI+/HXrFvDf/zKRJSIiIjISk9mqNH26+usXX1Q++JWRobwbyxqxRERERJXCaQZVqUULYMAAZTUCltYiIiIiEh2T2aq2fTsf6CIiIiKqIpxmUNWYyBIRERFVGSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLbMms3FxcWjfvj2kUim8vb0xePBgZGRk6NwmOTkZEomk3M/ly5dNFDURERERWQqzJrMpKSmYPHkyjh8/jv3796OkpAS9e/fGo0ePKtw2IyMDOTk5qp/GjRubIGIiIiIisiRmrTO7Z88etdcJCQnw9vbGmTNn0KVLF53bent7o2bNmlUYHRERERFZOouaM1tQUAAAqF27doVjQ0ND4efnhx49euDQoUNax8lkMhQWFqr9EBEREZFtsJhkVhAEREdH44UXXkCLFi20jvPz88OKFSuwadMmbN68GcHBwejRowcOHz6scXxcXBw8PT1VP/7+/lV1CERERERkYhJBEARzBwEAkydPxs6dO/HLL7+gXr16Bm07cOBASCQSbN++vdw6mUwGmUymel1YWAh/f38UFBTAw8Oj0nETERERkbgKCwvh6empV75mEXdmp06diu3bt+PQoUMGJ7IA0KFDB1y9elXjOhcXF3h4eKj9EBEREZFtMOsDYIIgYOrUqdiyZQuSk5MRFBRk1H5SU1Ph5+cncnREREREZOnMmsxOnjwZ33zzDbZt2wapVIrc3FwAgKenJ9zc3AAAMTExuH37NtatWwcAWLRoEQIDA9G8eXMUFxcjKSkJmzZtwqZNm8x2HERERERkHmZNZuPj4wEA3bp1U1uekJCAqKgoAEBOTg6ys7NV64qLizFz5kzcvn0bbm5uaN68OXbu3ImIiAhThU1EREREFsJiHgAzFUMmFBMRERGR6VndA2BERERERMZgMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdUyaztboqoiVwg4mXkPeQ+K4C11xXNBteHoIDF3WERERCQyJrNkc/ZczEHsjkvIKShSLfPzdMXsgSHo28LPjJERERGR2DjNgGzKnos5eCvprFoiCwC5BUV4K+ks9lzMMVNkREREVBWYzJLNkCsExO64BEHDutJlsTsuQa7QNIKIiIisEZNZshknM++VuyP7NAFATkERTmbeM11QREREVKWYzJLNyHugPZE1ZhwRERFZPiazZDO8pa6ijiMiIiLLx2SWbMZzQbXh5+kKbQW4JFBWNXguqLYpwyIiIqIqxGSWbIajgwSzB4YAQLmEtvT17IEhrDdLRERkQ5jMkk3p28IP8SPawtdTfSqBr6cr4ke0ZZ1ZIiIiG8OmCWRz+rbwQ68QX3YAIyIisgNMZskmOTpIEN7Qy9xhEBERURXjNAMiIiIislpMZomIiIjIajGZJSIiIiKrxWSWiIiIiKwWk1kiIiIislpMZomIiIjIatldaS5BEAAAhYWFZo6EiIiIiDQpzdNK8zZd7C6ZffDgAQDA39/fzJEQERERkS4PHjyAp6enzjESQZ+U14YoFAr88ccfkEqlkEhM0xGqsLAQ/v7+uHnzJjw8PEzynsTzbg485+bB824ePO/mwfNuHqY+74Ig4MGDB6hbty4cHHTPirW7O7MODg6oV6+eWd7bw8ODv3hmwPNuejzn5sHzbh487+bB824epjzvFd2RLcUHwIiIiIjIajGZJSIiIiKrxWTWBFxcXDB79my4uLiYOxS7wvNuejzn5sHzbh487+bB824elnze7e4BMCIiIiKyHbwzS0RERERWi8ksEREREVktJrNEREREZLWYzBIRERGR1WIyK6K4uDhIJBK8/fbbOselpKSgXbt2cHV1RYMGDbBs2TLTBGij9DnvycnJkEgk5X4uX75sukCt3Jw5c8qdP19fX53b8FqvPEPPO6918dy+fRsjRoyAl5cX3N3d0aZNG5w5c0bnNrzmK8/Q885rvvICAwM1nsPJkydr3caSrnW76wBWVU6dOoUVK1agVatWOsdlZmYiIiIC48ePR1JSEo4ePYpJkyahTp06GDp0qImitR36nvdSGRkZap1L6tSpU1Wh2aTmzZvj559/Vr12dHTUOpbXungMOe+leK1Xzl9//YVOnTqhe/fu2L17N7y9vXHt2jXUrFlT6za85ivPmPNeite88U6dOgW5XK56ffHiRfTq1QuvvvqqxvGWdq0zmRXBw4cPERkZiZUrV2LevHk6xy5btgwBAQFYtGgRAKBZs2Y4ffo0Pv/8c/5lZyBDznspb29vvf5SJM2cnJwqvBtbite6eAw576V4rVfOggUL4O/vj4SEBNWywMBAndvwmq88Y857KV7zxiub+M+fPx8NGzZE165dNY63tGud0wxEMHnyZPTv3x89e/ascOyxY8fQu3dvtWV9+vTB6dOn8eTJk6oK0SYZct5LhYaGws/PDz169MChQ4eqMDrbdPXqVdStWxdBQUF4/fXXcf36da1jea2Lx5DzXorXeuVs374dYWFhePXVV+Ht7Y3Q0FCsXLlS5za85ivPmPNeite8OIqLi5GUlIQxY8ZAIpFoHGNp1zqT2UrauHEjzp49i7i4OL3G5+bmwsfHR22Zj48PSkpKcPfu3aoI0SYZet79/PywYsUKbNq0CZs3b0ZwcDB69OiBw4cPV3GktuP555/HunXrsHfvXqxcuRK5ubno2LEj8vPzNY7ntS4OQ887r3VxXL9+HfHx8WjcuDH27t2LiRMnYtq0aVi3bp3WbXjNV54x553XvLi2bt2K+/fvIyoqSusYS7vWOc2gEm7evInp06dj3759cHV11Xu7sv/TKW3Cpu1/QKTOmPMeHByM4OBg1evw8HDcvHkTn3/+Obp06VJVodqUfv36qf7csmVLhIeHo2HDhli7di2io6M1bsNrvfIMPe+81sWhUCgQFhaGTz75BIDyrt9vv/2G+Ph4jBw5Uut2vOYrx5jzzmteXKtXr0a/fv1Qt25dneMs6VrnndlKOHPmDPLy8tCuXTs4OTnByckJKSkpWLJkCZycnNQmU5fy9fVFbm6u2rK8vDw4OTnBy8vLVKFbNWPOuyYdOnTA1atXqzha21W9enW0bNlS6znktV41KjrvmvBaN5yfnx9CQkLUljVr1gzZ2dlat+E1X3nGnHdNeM0b58aNG/j5558xbtw4neMs7VrnndlK6NGjBy5cuKC2bPTo0WjatCneffddjU8ch4eHY8eOHWrL9u3bh7CwMFSrVq1K47UVxpx3TVJTU+Hn51cVIdoFmUyG9PR0dO7cWeN6XutVo6LzrgmvdcN16tQJGRkZasuuXLmC+vXra92G13zlGXPeNeE1b5yEhAR4e3ujf//+OsdZ3LUukKi6du0qTJ8+XfX6vffeE958803V6+vXrwvu7u7CjBkzhEuXLgmrV68WqlWrJvz4449miNZ2VHTeFy5cKGzZskW4cuWKcPHiReG9994TAAibNm0yQ7TW6Z133hGSk5OF69evC8ePHxcGDBggSKVSISsrSxAEXutVxdDzzmtdHCdPnhScnJyEjz/+WLh69aqwYcMGwd3dXUhKSlKN4TUvPmPOO695ccjlciEgIEB49913y62z9GudyazIyiZVo0aNErp27ao2Jjk5WQgNDRWcnZ2FwMBAIT4+3rRB2qCKzvuCBQuEhg0bCq6urkKtWrWEF154Qdi5c6fpA7Viw4YNE/z8/IRq1aoJdevWFV5++WXht99+U63ntV41DD3vvNbFs2PHDqFFixaCi4uL0LRpU2HFihVq63nNVw1DzzuveXHs3btXACBkZGSUW2fp17pEEP6ZsUtEREREZGX4ABgRERERWS0ms0RERERktZjMEhEREZHVYjJLRERERFaLySwRERERWS0ms0RERERktZjMEhEREZHVYjJLRGQFoqKiMHjwYK3rExMTUbNmTZPFU5HAwEAsWrTI3GEQkR1gMktEREaztCSaiOwPk1kiIiIislpMZomIKvDjjz+iZcuWcHNzg5eXF3r27IlHjx6p1ickJKBZs2ZwdXVF06ZN8fXXX6vWZWVlQSKRYOPGjejYsSNcXV3RvHlzJCcnq8bI5XKMHTsWQUFBcHNzQ3BwMBYvXlzpuHfs2IF27drB1dUVDRo0QGxsLEpKSlTrJRIJVq1ahSFDhsDd3R2NGzfG9u3b1faxfft2NG7cGG5ubujevTvWrl0LiUSC+/fvIzk5GaNHj0ZBQQEkEgkkEgnmzJmj2vbx48cYM2YMpFIpAgICsGLFikofExFROQIREWn1xx9/CE5OTsIXX3whZGZmCufPnxe++uor4cGDB4IgCMKKFSsEPz8/YdOmTcL169eFTZs2CbVr1xYSExMFQRCEzMxMAYBQr1494ccffxQuXbokjBs3TpBKpcLdu3cFQRCE4uJiYdasWcLJkyeF69evC0lJSYK7u7vw3XffqeIYNWqUMGjQIK1xJiQkCJ6enqrXe/bsETw8PITExETh2rVrwr59+4TAwEBhzpw5qjGlcX3zzTfC1atXhWnTpgk1atQQ8vPzVbFXq1ZNmDlzpnD58mXh22+/FZ599lkBgPDXX38JMplMWLRokeDh4SHk5OQIOTk5qvNSv359oXbt2sJXX30lXL16VYiLixMcHByE9PR0UT4XIqJSTGaJiHQ4c+aMAEDIysrSuN7f31/45ptv1JbNnTtXCA8PFwThf8ns/PnzVeufPHki1KtXT1iwYIHW9500aZIwdOhQ1WtDk9nOnTsLn3zyidqY9evXC35+fqrXAIQPP/xQ9frhw4eCRCIRdu/eLQiCILz77rtCixYt1PbxwQcfqJJZTe9bqn79+sKIESNUrxUKheDt7S3Ex8drPQYiImM4mfGmMBGRxWvdujV69OiBli1bok+fPujduzdeeeUV1KpVC3/++Sdu3ryJsWPHYvz48aptSkpK4Onpqbaf8PBw1Z+dnJwQFhaG9PR01bJly5Zh1apVuHHjBv7++28UFxejTZs2Rsd95swZnDp1Ch9//LFqmVwuR1FRER4/fgx3d3cAQKtWrVTrq1evDqlUiry8PABARkYG2rdvr7bf5557Tu8Ynt63RCKBr6+vat9ERGJhMktEpIOjoyP279+PX3/9Ffv27cPSpUvxwQcf4MSJE6qEcOXKlXj++efLbVcRiUQCAPj+++8xY8YM/Pe//0V4eDikUik+++wznDhxwui4FQoFYmNj8fLLL5db5+rqqvpztWrVysWkUCgAAIIgqGIsJQiC3jHo2jcRkViYzBIRVUAikaBTp07o1KkTZs2ahfr162PLli2Ijo7Gs88+i+vXryMyMlLnPo4fP44uXboAUN65PXPmDKZMmQIAOHLkCDp27IhJkyapxl+7dq1SMbdt2xYZGRlo1KiR0fto2rQpdu3apbbs9OnTaq+dnZ0hl8uNfg8iospiMktEpMOJEydw4MAB9O7dG97e3jhx4gT+/PNPNGvWDAAwZ84cTJs2DR4eHujXrx9kMhlOnz6Nv/76C9HR0ar9fPXVV2jcuDGaNWuGhQsX4q+//sKYMWMAAI0aNcK6deuwd+9eBAUFYf369Th16hSCgoKMjnvWrFkYMGAA/P398eqrr8LBwQHnz5/HhQsXMG/ePL32MWHCBHzxxRd49913MXbsWKSlpSExMRHA/+4qBwYG4uHDhzhw4ABat24Nd3d31R1rIiJTYGkuIiIdPDw8cPjwYURERKBJkyb48MMP8d///hf9+vUDAIwbNw6rVq1CYmIiWrZsia5duyIxMbFcIjp//nwsWLAArVu3xpEjR7Bt2zY888wzAICJEyfi5ZdfxrBhw/D8888jPz9f7S6tMfr06YOffvoJ+/fvR/v27dGhQwd88cUXqF+/vt77CAoKwo8//ojNmzejVatWiI+PxwcffAAAcHFxAQB07NgREydOxLBhw1CnTh18+umnlYqbiMhQEsGQCVBERGSQrKwsBAUFITU1tVIPdFmKjz/+GMuWLcPNmzfNHQoREQBOMyAiIh2+/vprtG/fHl5eXjh69Cg+++wz1VxfIiJLwGSWiIi0unr1KubNm4d79+4hICAA77zzDmJiYswdFhGRCqcZEBEREZHV4gNgRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1/h8krms46M4PUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1, figsize=(8, 4))\n",
    "plt.clf()\n",
    "plt.scatter(df[df[\"species\"] == \"setosa\"].sepal_length, df[df[\"species\"] == \"setosa\"].sepal_width,label=\"Iris Data\")\n",
    "plt.plot(X_plot, fit, label=\"Linear Regression Model\", color=\"red\", linewidth=3)\n",
    "plt.xlabel(\"sepal length\")\n",
    "plt.ylabel(\"sepal width\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d6a6e62",
   "metadata": {},
   "source": [
    "## Alternativ framgangsmåte\n",
    "\n",
    "Vi kan også bruke et annet bibliotek i python ([statsmodels](https://www.statsmodels.org/stable/index.html)) som er et godt alternativ til scikit-learn om vi enkelt vil hente ut enda mer informasjon om modellen vår. Da må vi starte med å importere denne "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "92bc73d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32b3af73",
   "metadata": {},
   "source": [
    "Vi kan bruke OLS (**o**rdinary **l**east **s**quare) på samme måte som vi brukte scikit-learn metoden over. Den gir oss da tilgang til enda mer informasjon om modellen. Merk at intercept/$\\beta_0$ ikke er inkludert i standard-oppsettet i *statsmodels.ols()* og derfor eksplisitt må legges til ved å bruke *sm.add_constant()* på respons-vektoren."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2eb976e9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.551\n",
      "Model:                            OLS   Adj. R-squared:                  0.542\n",
      "Method:                 Least Squares   F-statistic:                     58.99\n",
      "Date:                Fri, 26 Jan 2024   Prob (F-statistic):           6.71e-10\n",
      "Time:                        09:53:01   Log-Likelihood:                -1.9002\n",
      "No. Observations:                  50   AIC:                             7.800\n",
      "Df Residuals:                      48   BIC:                             11.62\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.5694      0.522     -1.091      0.281      -1.618       0.480\n",
      "x1             0.7985      0.104      7.681      0.000       0.589       1.008\n",
      "==============================================================================\n",
      "Omnibus:                        0.680   Durbin-Watson:                   2.345\n",
      "Prob(Omnibus):                  0.712   Jarque-Bera (JB):                0.342\n",
      "Skew:                          -0.200   Prob(JB):                        0.843\n",
      "Kurtosis:                       3.060   Cond. No.                         75.0\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "3.1586748135738\n"
     ]
    }
   ],
   "source": [
    "X_setosa_sepal_length_const = sm.add_constant(X_setosa_sepal_length)\n",
    "\n",
    "model = sm.OLS(X_setosa_sepal_width,X_setosa_sepal_length_const).fit() \n",
    "  \n",
    "#display model summary \n",
    "print(model.summary()) \n",
    "  \n",
    "# residual sum of squares \n",
    "print(model.ssr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "26ed067a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definer en x-vektor som dekker datapunktene (i dette tilfellet mellom 4 og 7 cm)\n",
    "X_plot = np.linspace(4, 7, 300)\n",
    "# Lag så responsfunksjonen ved å bruke parameterne fra modellen vår\n",
    "fit = model.params[0] + model.params[1]*X_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "52ae8a70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff0940178d0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAFzCAYAAAAt54EyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZj0lEQVR4nO3deVxU9foH8M8AsumASrKYILihuKJYornljpqalSWmuF3NNcn7K6qrkhZa3dwq3EHFssU991TQzF1wSURTEDWQxASXGGTm/P6YmOvAzDAzHGb9vF8vXi/nnO8585wzB3088z3PIxEEQQARERERkRVyMHcARERERETGYjJLRERERFaLySwRERERWS0ms0RERERktZjMEhEREZHVYjJLRERERFaLySwRERERWS0ms0RERERktZzMHYCpKRQK/PHHH5BKpZBIJOYOh4iIiIjKEAQBDx48QN26deHgoPveq90ls3/88Qf8/f3NHQYRERERVeDmzZuoV6+ezjF2l8xKpVIAypPj4eFh5miIiIiIqKzCwkL4+/ur8jZd7C6ZLZ1a4OHhwWSWiIiIyILpMyWUD4ARERERkdViMktEREREVovJLBERERFZLbubM6sPQRBQUlICuVxu7lCIbJKjoyOcnJxYHo+IiCqNyWwZxcXFyMnJwePHj80dCpFNc3d3h5+fH5ydnc0dChERWTEms09RKBTIzMyEo6Mj6tatC2dnZ945IhKZIAgoLi7Gn3/+iczMTDRu3LjCgthERETaMJl9SnFxMRQKBfz9/eHu7m7ucIhslpubG6pVq4YbN26guLgYrq6u5g6JiIisFG+HaMC7RERVj79nRERWSBDMHUE5/NeEiIiIiLQTBODIEWDoUGD2bHNHU45Zk9k5c+ZAIpGo/fj6+urcJiUlBe3atYOrqysaNGiAZcuWmShaIiIiIjtSXAysXw+EhQFdugCbNwNffw1Y2EPyZr8z27x5c+Tk5Kh+Lly4oHVsZmYmIiIi0LlzZ6SmpuL999/HtGnTsGnTJhNGbN0kEgm2bt1q7jCIiIjIUuXlAXPnAvXrAyNHAmfP/m9dfj6wYYP5YtPA7Mmsk5MTfH19VT916tTROnbZsmUICAjAokWL0KxZM4wbNw5jxozB559/bsKILVNUVBQGDx5c4bicnBz069fPqPd4+k66k5MTnnnmGXTp0gWLFi2CTCYzaF/JycmQSCS4f/++UbEQERGRyM6fB8aOBQICgFmzgNxczeNWrzZtXBUwezJ79epV1K1bF0FBQXj99ddx/fp1rWOPHTuG3r17qy3r06cPTp8+jSdPnmjcRiaTobCwUO3HFOQKAceu5WNb2m0cu5YPucK8E6aLi4sBAL6+vnBxcTF6P6V30rOzs3Ho0CG8+uqriIuLQ8eOHfHgwQOxwiUiIiJTUCiAHTuAHj2A1q2BNWsAbTeoGjQAFi0C9u0zaYgVMWsy+/zzz2PdunXYu3cvVq5cidzcXHTs2BH5+fkax+fm5sLHx0dtmY+PD0pKSnD37l2N28TFxcHT01P14+/vL/pxlLXnYg5eWHAQb6w8jukb0/DGyuN4YcFB7LmYU+XvXapbt26YMmUKoqOj8cwzz6BXr14A1KcZFBcXY8qUKfDz84OrqysCAwMRFxenc7+ld9Lr1q2Lli1bYurUqUhJScHFixexYMEC1bikpCSEhYVBKpXC19cXw4cPR15eHgAgKysL3bt3BwDUqlULEokEUVFRAIA9e/bghRdeQM2aNeHl5YUBAwbg2rVrIp8dIiIiO/fgAbB0KRAcDLz0EnDwoPax3boBW7cCV64A06cDHh6milIvZk1m+/Xrh6FDh6Jly5bo2bMndu7cCQBYu3at1m3KNjEQ/ikRoa25QUxMDAoKClQ/N2/eFCl6zfZczMFbSWeRU1Cktjy3oAhvJZ01aUK7du1aODk54ejRo1i+fHm59UuWLMH27dvx/fffIyMjA0lJSQgMDDT4fZo2bYp+/fph8+bNqmXFxcWYO3cuzp07h61btyIzM1OVsPr7+6vmOWdkZCAnJweLFy8GADx69AjR0dE4deoUDhw4AAcHBwwZMgQKhcLwE0BERETqMjOBd94B6tUDpk0Dfv9d8zhnZyAqCkhNBQ4dAgYNAhwdTRqqviyqaUL16tXRsmVLXL16VeN6X19f5JaZv5GXlwcnJyd4eXlp3MbFxaVSX6sbQq4QELvjEjRNKBAASADE7riEXiG+cHSo+s5ijRo1wqeffqp1fXZ2Nho3bowXXngBEokE9evXN/q9mjZtin1Pfe0wZswY1Z8bNGiAJUuW4LnnnsPDhw9Ro0YN1K5dGwDg7e2NmjVrqsYOHTpUbb+rV6+Gt7c3Ll26hBYtWhgdHxERkd0SBOCXX5RTBLZuVU4t0MbbG5g0CZg4ESjzbbilMvuc2afJZDKkp6fDz89P4/rw8HDs379fbdm+ffsQFhaGatWqmSJEnU5m3it3R/ZpAoCcgiKczLxnknjCwsJ0ro+KikJaWhqCg4Mxbdo0tWTUUIIgqN0dT01NxaBBg1C/fn1IpVJ069YNgDKB1uXatWsYPnw4GjRoAA8PDwQFBem1HREREZWhqbSWtkS2dWsgMRHIzlbWkrWSRBYwczI7c+ZMpKSkIDMzEydOnMArr7yCwsJCjBo1CoByisDIkSNV4ydOnIgbN24gOjoa6enpWLNmDVavXo2ZM2ea6xDU5D3QnsgaM66yqlevrnN927ZtkZmZiblz5+Lvv//Ga6+9hldeecWo90pPT1clno8ePULv3r1Ro0YNJCUl4dSpU9iyZQuA/z2Ips3AgQORn5+PlStX4sSJEzhx4oRe2xEREdE/dJXWeppEAgweDCQnK6cTjBoFmOjbbDGZdZrBrVu38MYbb+Du3buoU6cOOnTogOPHj6u+7i59ar5UUFAQdu3ahRkzZuCrr75C3bp1sWTJknJfTZuLt1S//vL6jjMFDw8PDBs2DMOGDcMrr7yCvn374t69e6ppAPq4fPky9uzZg5iYGNXru3fvYv78+aoH7k6fPq22jbOzMwBALperluXn5yM9PR3Lly9H586dAQC//PJLpY6PiIjIbpw/DyxerKwDq6tkplSqLME1ZQrQsKHp4qsiZk1mN27cqHN9YmJiuWVdu3bFWW3/wzCz54Jqw8/TFbkFRRrnzUoA+Hq64rkg/RPFqrRw4UL4+fmhTZs2cHBwwA8//ABfX1+1OaxllZSUIDc3FwqFAvn5+UhOTsa8efPQpk0b/Pvf/wYABAQEwNnZGUuXLsXEiRNx8eJFzJ07V20/9evXh0QiwU8//YSIiAi4ubmhVq1a8PLywooVK+Dn54fs7Gy89957VXkKiIiIrJtCAezcqZwPq6siAaAsrTVtGjB6tMVVJKgMi5oza+0cHSSYPTAEgDJxfVrp69kDQ0zy8Jc+atSogQULFiAsLAzt27dHVlYWdu3aBQcH7ZfFb7/9Bj8/PwQEBKBbt274/vvvERMTgyNHjqBGjRoAgDp16iAxMRE//PADQkJCMH/+/HKNLZ599lnExsbivffeg4+PD6ZMmQIHBwds3LgRZ86cQYsWLTBjxgx89tlnVXoOiIiIrJINldaqLIlQWtvKThQWFsLT0xMFBQXwKPNhFhUVITMzE0FBQXB1NX4qwJ6LOYjdcUntYTA/T1fMHhiCvi00P9xGZG/E+n0jIrIrmZnAl18Cq1YBuhpBOTsDw4crk9c2bUwWnlh05WtlWVRpLlvRt4UfeoX44mTmPeQ9KIK3VDm1wFLuyBIREZEVsfHSWpXFZLaKODpIEN5Qc+1bIiIiogrJZMD33yuT2IqeF2rdGpgxA3j9dausSFAZTGaJiIiILEleHrB8OfD110CZZlFqJBJlZ66331bWkdXSDdXWMZklIiIisgSGltaaOlVZocDOMZklIiIiMheW1qo0JrNEREREpvbggbJ97JIlwO+/6x7brZtyKsGAAYCjowmCsy5MZomIiIhMxU5Ka5kSk1kiIiKiqsTSWlWKHcDshEQiwdatW80dht1JTEzU2R7YmhhzLLzuiMiuyWTA+vVAWJiy2sDmzdoT2TZtlNMOsrOB2bOZyBqAyayNiIqKwuDBg7Wuz8nJQb9+/UwXkIEkEonqp0aNGmjdujUSExPNHValDRs2DFeuXKny9+nWrRskEgnmz59fbl1ERAQkEgnmzJlT5XEQERGUpbXmzgUCA4GRI7XXiJVIgMGDgeRk5ZhRo+yuRqwYmMzaCV9fX7iY+RdEEASUlJRoXZ+QkICcnBycO3cOw4YNw+jRo7F3794qjam4uLhK9+/m5gZvb+8qfY9S/v7+SEhIUFv2xx9/4ODBg/DzYxtlIqIqd/68smRWQAAwa5b2GrFSqfKBrt9/B7ZsAbp2tdsasWJgMquLQgH8+ad5f3TNqzHA01/3ZmVlQSKRYPPmzejevTvc3d3RunVrHDt2TG2bX3/9FV26dIGbmxv8/f0xbdo0PHr0SLU+KSkJYWFhkEql8PX1xfDhw5GXl6dan5ycDIlEgr179yIsLAwuLi44cuSI1hhr1qwJX19fNGzYEO+//z5q166Nffv2qdYXFBTgX//6F7y9veHh4YEXX3wR586dU9vHvHnz4O3tDalUinHjxuG9995Dm6cmzpfewY6Li0PdunXRpEkTAMDt27cxbNgw1KpVC15eXhg0aBCysrLUjuW5555D9erVUbNmTXTq1Ak3btwAAJw7dw7du3eHVCqFh4cH2rVrh9OnTwPQ/NV8fHw8GjZsCGdnZwQHB2P9+vXlPqtVq1ZhyJAhcHd3R+PGjbF9+3at563UgAEDkJ+fj6NHj6qWJSYmonfv3uUS6r/++gsjR45ErVq14O7ujn79+uHq1atqYxITExEQEAB3d3cMGTIE+fn55d5zx44daNeuHVxdXdGgQQPExsbq/A8LEZHNUSiAHTuAHj2UXbjWrNFeI7ZBA+W82Vu3gIULWSNWJExmdcnPV07ENuePhgRCLB988AFmzpyJtLQ0NGnSBG+88YYqEblw4QL69OmDl19+GefPn8d3332HX375BVOmTFFtX1xcjLlz5+LcuXPYunUrMjMzERUVVe59/u///g9xcXFIT09Hq1atKoxLLpfj+++/x71791CtWjUAyru6/fv3R25uLnbt2oUzZ86gbdu26NGjB+7duwcA2LBhAz7++GMsWLAAZ86cQUBAAOLj48vt/8CBA0hPT8f+/fvx008/4fHjx+jevTtq1KiBw4cP45dffkGNGjXQt29fFBcXo6SkBIMHD0bXrl1x/vx5HDt2DP/6178g+ed/0ZGRkahXrx5OnTqFM2fO4L333lPFXdaWLVswffp0vPPOO7h48SImTJiA0aNH49ChQ2rjYmNj8dprr+H8+fOIiIhAZGSk6ji1cXZ2RmRkpNrd2cTERIwZM6bc2KioKJw+fRrbt2/HsWPHIAgCIiIi8OTJEwDAiRMnMGbMGEyaNAlpaWno3r075s2bp7aPvXv3YsSIEZg2bRouXbqE5cuXIzExER9//LHOOImIbMKDB8DSpUBwMPDSS7prxHbrpnzw68oVZXUC1ogVl2BnCgoKBABCQUFBuXV///23cOnSJeHvv/9WLsjLEwTlM4jm+8nL0+u4Ro0aJQwaNEjregDCli1bBEEQhMzMTAGAsGrVKtX63377TQAgpKenC4IgCG+++abwr3/9S20fR44cERwcHP53fso4efKkAEB48OCBIAiCcOjQIQGAsHXr1grjByC4uroK1atXFxwdHQUAQu3atYWrV68KgiAIBw4cEDw8PISioiK17Ro2bCgsX75cEARBeP7554XJkyerre/UqZPQunVr1etRo0YJPj4+gkwmUy1bvXq1EBwcLCgUCtUymUwmuLm5CXv37hXy8/MFAEJycrLG2KVSqZCYmKhxXUJCguDp6al63bFjR2H8+PFqY1599VUhIiJC7Vx8+OGHqtcPHz4UJBKJsHv3bo3vIQiC0LVrV2H69OnCuXPnBKlUKjx8+FBISUkRvL29heLiYqF169bC7NmzBUEQhCtXrggAhKNHj6q2v3v3ruDm5iZ8//33giAIwhtvvCH07dtX7T2GDRumdiydO3cWPvnkE7Ux69evF/z8/NSOpfS6K6vc7xsRkTW4fl0QoqMFwcND97/fzs6CEBUlCKmp5o7YKunK18rinVk79vRd0tI5laXTBM6cOYPExETUqFFD9dOnTx8oFApkZmYCAFJTUzFo0CDUr18fUqkU3bp1AwBkZ2ervU9YWJhe8SxcuBBpaWnYv38/2rRpg4ULF6JRo0aqeB4+fAgvLy+1mDIzM3Ht2jUAQEZGBp577jm1fZZ9DQAtW7aEs7Oz6vWZM2fw+++/QyqVqvZbu3ZtFBUV4dq1a6hduzaioqLQp08fDBw4EIsXL0ZOTo5q++joaIwbNw49e/bE/PnzVfFokp6ejk6dOqkt69SpE9LT09WWPf3ZVK9eHVKpVG0KhzatWrVC48aN8eOPP2LNmjV48803y90lTk9Ph5OTE55//nnVMi8vLwQHB6viSE9PR3h4uNp2ZV+fOXMGH330kdrnMX78eOTk5ODx48cVxkpEZDUEAThyBBg6FGjUCPjiC+01Yr29gTlzlFUJEhJYI9YEWGfWjj2d5JR+Za74Z46uQqHAhAkTMG3atHLbBQQE4NGjR+jduzd69+6NpKQk1KlTB9nZ2ejTp0+5h6qqV6+uVzy+vr5o1KgRGjVqhB9++AGhoaEICwtDSEgIFAoF/Pz8kJycXG67p+ekSspMoBcEodz4svEoFAq0a9cOGzZsKDe2Tp06AJQPp02bNg179uzBd999hw8//BD79+9Hhw4dMGfOHAwfPhw7d+7E7t27MXv2bGzcuBFDhgzReJyaYiy7rGwCKpFIVJ9NRcaMGYOvvvoKly5dwsmTJ8ut13ROysahbczTFAoFYmNj8fLLL5db5+rqqlesREQWTSYDvv9eOc9VW0WCUm3aKB/qev11ViQwMSazunh5KctrmDsGM2jbti1+++031Z3Rsi5cuIC7d+9i/vz58Pf3BwDVQ09iaNSoEYYOHYqYmBhs27YNbdu2RW5uLpycnBAYGKhxm+DgYJw8eRJvvvmmapk+MbVt2xbfffed6sEybUJDQxEaGoqYmBiEh4fjm2++QYcOHQAATZo0QZMmTTBjxgy88cYbSEhI0JjMNmvWDL/88gtGjhypWvbrr7+iWbNmFcapr+HDh2PmzJlo3bo1QkJCyq0PCQlBSUkJTpw4gY4dOwIA8vPzceXKFVUcISEhOH78uNp2ZV+3bdsWGRkZWq8RIiKrlZcHLF8OfP219ooEgLICwaBByiS2SxdWJDATJrO6ODgA/9yZswYFBQVIS0tTW1a7dm0EBAQYvK93330XHTp0wOTJkzF+/HhUr15d9dDU0qVLERAQAGdnZyxduhQTJ07ExYsXMXfuXJGOROmdd95B69atcfr0afTs2RPh4eEYPHgwFixYgODgYPzxxx/YtWsXBg8ejLCwMEydOhXjx49HWFgYOnbsiO+++w7nz59HgwqeFo2MjMRnn32GQYMG4aOPPkK9evWQnZ2NzZs349///jeePHmCFStW4KWXXkLdunWRkZGBK1euYOTIkfj777/x73//G6+88gqCgoJw69YtnDp1CkOHDtX4Xv/+97/x2muvqR5e27FjBzZv3oyff/5ZtPNWq1Yt5OTkaH0IrXHjxhg0aBDGjx+P5cuXQyqV4r333sOzzz6LQYMGAQCmTZuGjh074tNPP8XgwYOxb98+7NmzR20/s2bNwoABA+Dv749XX30VDg4OOH/+PC5cuFDuYTEiIqtw/jyweDGwYYP2igSAsrTW2LHA1KmsSGABOGfWhiQnJ6vuHpb+zJo1y6h9tWrVCikpKbh69So6d+6M0NBQ/Oc//1HNra1Tpw4SExPxww8/ICQkBPPnz8fnn38u5uGgZcuW6NmzJ2bNmgWJRIJdu3ahS5cuGDNmDJo0aYLXX38dWVlZ8PmnS0pkZCRiYmIwc+ZMtG3bVlVdoaKvvN3d3XH48GEEBATg5ZdfRrNmzTBmzBj8/fff8PDwgLu7Oy5fvoyhQ4eiSZMm+Ne//oUpU6ZgwoQJcHR0RH5+PkaOHIkmTZrgtddeQ79+/RAbG6vxvQYPHozFixfjs88+Q/PmzbF8+XIkJCSo5huLpWbNmjqndyQkJKBdu3YYMGAAwsPDIQgCdu3apUqAO3TogFWrVmHp0qVo06YN9u3bhw8//FBtH3369MFPP/2E/fv3o3379ujQoQO++OIL1K9fX9RjISKqUiytZfUkgj6T42xIYWEhPD09UVBQUO4r5aKiImRmZiIoKIhz/mxEr1694OvrW66WK5kff9+IyKwePFC2j12yRNm8QJdu3ZRTCQYMABwdTRAc6crXyuI0A7IZjx8/xrJly9CnTx84Ojri22+/xc8//4z9+/ebOzQiIrIUmZnAl18Cq1Zpr0gAAM7OwPDhyrqwrEhg0ZjMks0onYowb948yGQyBAcHY9OmTejZs6e5QyMiInMSBOCXX5RTBLZu1d1d09sbmDQJmDgR+GcaG1k2JrNkM9zc3ER9kIqIiKwcS2vZBSazREREZFtYWsuuMJklIiIi28DSWnaJyawGdlbggcgs+HtGRKJQKICdO5VTCQ4e1D22QQNg2jRg9GiggifkyXowmX1KaY3Nx48fw83NzczRENm2x48fAyjfupeISC+lpbUWLwauXdM9lqW1bBqT2ac4OjqiZs2ayPunha27u7uqVz0RiUMQBDx+/Bh5eXmoWbMmHPkPCxEZgqW1qAwms2X4+voCgCqhJaKqUbNmTdXvGxGRTiytRTowmS1DIpHAz88P3t7eePLkibnDIbJJ1apV4x1ZIqoYS2uRHpjMauHo6Mh/bImIiMyBpbXIAExmiYiIyDKwtBYZgcksERERmY9c/r/SWocO6R7L0lqkAZNZIiIiMj2W1iKRMJklIiIi02FpLRIZk1kiIiKqWiytRVWIySwRERFVDZkM+O475VQCltaiKsJkloiIiMSVlwcsW6YsrXXnjvZxLK1FInAwdwCl4uLiIJFI8Pbbb2sdk5ycDIlEUu7n8uXLpguUiIiINDt3DhgzBggIAGbP1p7ISqXKBPb334EtW4CuXZnIktEs4s7sqVOnsGLFCrRq1Uqv8RkZGfB4qiRHnTp1qio0IiIi0oWltcjMzJ7MPnz4EJGRkVi5ciXmzZun1zbe3t6oWbNm1QZGRERE2rG0FlkIs08zmDx5Mvr374+ePXvqvU1oaCj8/PzQo0cPHKrgf4EymQyFhYVqP0RERGSkzEzgnXeAevWUd1m1JbLOzkBUFJCaqrxjO2gQE1mqEma9M7tx40acPXsWp06d0mu8n58fVqxYgXbt2kEmk2H9+vXo0aMHkpOT0aVLF43bxMXFITY2VsywiYiI7IsgAEeOKKcSbNvG0lpkUSSCIAjmeOObN28iLCwM+/btQ+vWrQEA3bp1Q5s2bbBo0SK99zNw4EBIJBJs375d43qZTAbZU/2dCwsL4e/vj4KCArV5t0RERFRGaWmtRYuUd1h1YWktElFhYSE8PT31ytfMdmf2zJkzyMvLQ7t27VTL5HI5Dh8+jC+//BIymQyOenwd0aFDByQlJWld7+LiAhf+UhEREemPpbXIipgtme3RowcuXLigtmz06NFo2rQp3n33Xb0SWQBITU2Fn59fVYRIRERkX86dUz7Q9c03yruy2kilwNixwNSpygoFRGZktmRWKpWiRYsWasuqV68OLy8v1fKYmBjcvn0b69atAwAsWrQIgYGBaN68OYqLi5GUlIRNmzZh06ZNJo+fiIjIJrC0Flk5s5fm0iUnJwfZ2dmq18XFxZg5cyZu374NNzc3NG/eHDt37kRERIQZoyQiIrJCLK1FNsJsD4CZiyETiomIyDhyhYCTmfeQ96AI3lJXPBdUG44OnE9pETIzgaVLgdWrAV3lKp2dgeHDgenTlQ93EZmQVTwARkREtmnPxRzE7riEnIIi1TI/T1fMHhiCvi34jINZsLQW2TAms0REJJo9F3PwVtJZlP3KL7egCG8lnUX8iLZMaE2JpbXIDjCZJSIiUcgVAmJ3XCqXyAKAAEACIHbHJfQK8eWUg6rG0lpkR5jMEhGRKE5m3lObWlCWACCnoAgnM+8hvKGX6QKzJyytRXaIySwREYki74H2RNaYcaQnltYiO8dkloiIROEtdRV1HFXgwQMgIQFYsoSltciuMZklIiJRPBdUG36ersgtKNI4b1YCwNdTWaaLKoGltYjUOJg7ACIisg2ODhLMHhgCQJm4Pq309eyBIXz4yxiCABw+DLz8MtCoEbBwofZE1tsbmDMHyM5W3rllIks2jsksERGJpm8LP8SPaAtfT/WpBL6erizLZQyZDFi3DmjXDujaFdiyRXuN2DZtgLVrlUns7NmsEUt2g9MMiIhIVH1b+KFXiC87gFUGS2sR6Y3JLBERic7RQcLyW8ZgaS0igzGZJSIiMieW1iKqFCazRERE5sDSWkSiYDJLRERkSiytRSQqJrNERERVTRCAI0eUUwm2bdNekQBQViGYNAmYMIEVCYj0wGSWiIioqshkwHffKZPY1FTdY9u0AWbMAIYNA1xcTBEdkU1gMktERCQ2ltYiMhkms0RERGJhaS0ik2MyS0REVBksrUVkVkxmiYismFwhsNOWuRhSWqt7d2VVApbWIhIdk1kiIiu152IOYndcQk5BkWqZn6crZg8MQd8WfmaMzMYZUlorMlKZxLZubbr4iOwMk1kiIiu052IO3ko6C6HM8tyCIryVdBbxI9oyoRUTS2sRWSwms0REVkauEBC741K5RBYABAASALE7LqFXiC+nHFQWS2sRWTwms0REVuZk5j21qQVlCQByCopwMvMewht6mS4wW8LSWkRWg8ksEZGVyXugPZE1Zhw9pbS01oYNQHGx9nEsrUVkMZjMEhFZGW+pq6jj7B5LaxFZNSazRERW5rmg2vDzdEVuQZHGebMSAL6eyjJdpIOhpbXefhvo35+ltYgsjIO5AyAiIsM4Okgwe2AIAGXi+rTS17MHhvDhL20yM4HoaKBePWXZLG2JrLOz8g5sWhpw8CDw0ktMZIksEJNZIiIr1LeFH+JHtIWvp/pUAl9PV5bl0kQQgMOHgZdfBho1AhYu1F4j1scHiI0FsrOBNWtYI5bIwnGaARHZNFvukNW3hR96hfja7PGJgqW1iGwek1kisln20CHL0UHC8luasLQWkd1gMktENokdsuyUIaW1xo0DpkxhaS0iK8dklohsDjtk2RlDS2tNnw5ERbG0FpGNYDJLRDaHHbLsBEtrERGYzBKRDWKHLBuXmQksXQqsXq29IgGgLK0VGam8E8uKBEQ2i8ksEdkcdsiyQYIAHDminEqwbRugUGgf6+MDTJoETJig/DMR2TQms0Rkc9ghy4awtBYRVYBNE4jI5rBDlg3IywM++gioXx8YNUp7IiuRAEOGACkpwNmzwMiRTGSJ7AyTWSKySeyQZaXOnQPGjAH8/YHZs7XXiJVKlXdhf/8d2LyZNWKJ7BinGRCRzWKHLCvB0lpEVAkGJ7OPHj3C/PnzceDAAeTl5UFRZhL+9evXRQuOiKiy2CHLghlQWqugwwu4PGwMFBH98VyjOvwPCRGpGJzMjhs3DikpKXjzzTfh5+cHiUhf68TFxeH999/H9OnTsWjRIq3jUlJSEB0djd9++w1169bF//3f/2HixImixEBERCZgQGmtWxFDEOP/Io64PwvkAlhzyuZaEhNR5RiczO7evRs7d+5Ep06dRAvi1KlTWLFiBVq1aqVzXGZmJiIiIjB+/HgkJSXh6NGjmDRpEurUqYOhQ4eKFg8REYnMiNJaB7sMxtg9N9mSmIh0MvgBsFq1aqF2bfHK2Tx8+BCRkZFYuXIlatWqpXPssmXLEBAQgEWLFqFZs2YYN24cxowZg88//1y0eIiISEQyGbBuHdCuHdC1K7Bli/ZEtk0bYO1a4MYNyD/8Dz449qfWlsSAsiWxXKFpBBHZE4OT2blz52LWrFl4/PixKAFMnjwZ/fv3R8+ePSsce+zYMfTu3VttWZ8+fXD69Gk8efJE4zYymQyFhYVqP0REVMUqWVrLkJbERGTf9JpmEBoaqjY39vfff4ePjw8CAwNRrVo1tbFnz57V+803btyIs2fP4tSpU3qNz83NhU+Zbi4+Pj4oKSnB3bt34edX/uumuLg4xMbG6h0TERFVwrlzwOLFwIYNQHGx9nFSKTBuHDBlirJCQRlsSUxE+tIrmR08eLDob3zz5k1Mnz4d+/btg6ur/i0lyz5wJgiCxuWlYmJiEB0drXpdWFgIf39/IyImIiKNqqC0FlsSE5G+9EpmZ8+eLfobnzlzBnl5eWjXrp1qmVwux+HDh/Hll19CJpPB0dFRbRtfX1/k5uaqLcvLy4OTkxO8vDSX3nFxcYELu8EQEYnPgNJa6N4dePttoH9/oMzf7ZqwJTER6cvgObMNGjRAfn5+ueX3799HAw1fFWnTo0cPXLhwAWlpaaqfsLAwREZGIi0trVwiCwDh4eHYv3+/2rJ9+/YhLCys3HQHIiKqIpmZQHQ0UK+e8i6rtkTW2RkYPRpISwMOHgReekmvRBZgS2Ii0p/BpbmysrIgl8vLLZfJZLh165be+5FKpWjRooXasurVq8PLy0u1PCYmBrdv38a6desAABMnTsSXX36J6OhojB8/HseOHcPq1avx7bffGnoYRERkCCNKa2HCBOWfjVTakjh2xyW1h8F8WWeWiJ6idzK7fft21Z/37t0LT09P1Wu5XI4DBw4gKChI1OBycnKQnZ2teh0UFIRdu3ZhxowZ+Oqrr1C3bl0sWbKENWaJyOoUlyiw/lgWbtx7jPq13fFmeCCcnQz+sqzqyWTAd98pk1htFQlKtWkDzJgBDBsGiDS9iy2JiagiEqH0CaoKODgo/5KVSCQou0m1atUQGBiI//73vxgwYID4UYqosLAQnp6eKCgogAf7ehORGcTtuoSVRzLxdIlUBwkwvnMQYiJCzBfY0/LygGXLgK+/Bu7c0T5OIgEGD1bOh+3cWfmaiKiSDMnX9L4zq/jnK6WgoCCcOnUKzzzzTOWiJCKyQ3G7LmH54cxyyxUCVMvNmtCKVFqLiMhUDJ4zm5lZ/i9hIiKqWHGJAiuP6P47dOWRTLzTu6lppxxUQWktIiJTMTiZXbJkicblEokErq6uaNSoEbp06aKxGgERkT1bfywLFXVfVQjKcWM7m+BuZxWW1iIiMhWDk9mFCxfizz//xOPHj1GrVi0IgoD79+/D3d0dNWrUQF5eHho0aIBDhw6xOQER0VNu3NOvDbi+44yWmQksXQqsXg3oavHt7AxERirvxLZuXbUxEREZyeDvsT755BO0b98eV69eRX5+Pu7du4crV67g+eefx+LFi5GdnQ1fX1/MmDGjKuIlIrJa9Wu7izrOIIIAHD4MvPwy0KgRsHCh9kTWxweIjQWys4E1a5jIEpFF07uaQamGDRti06ZNaNOmjdry1NRUDB06FNevX8evv/6KoUOHIicnR8xYRcFqBkRkLsUlCjT9z26dUw0cJMDluf3EmzNr5tJaRETGqJJqBqVycnJQUlJSbnlJSYmq1WzdunXx4MEDQ3dNRGTTnJ0cML5zkMZqBqXGdw4SJ5FlaS0ishMG/43ZvXt3TJgwAalP/Q8/NTUVb731Fl588UUAwIULF0RvoEBEZAtiIkIwoUsQytb8d5AAE7qIUGf23DlgzBjA3x+YPVt7IiuVKu/C/v47sHkz0KULE1kiskoGTzPIzc3Fm2++iQMHDqBatWoAlHdle/TogfXr18PHxweHDh3CkydP0Lt37yoJujI4zYDI8v1dLMcnuy4hK/8xAr3c8X5ECNycbesJelE7gBlSWqthQ2DatCovrSVXCOzaRURGMyRfMziZLXX58mVcuXIFgiCgadOmCA4ONipYU2MyS2TZxq87hf2X8sot7xXijZUj25shIgtmoaW19lzMQeyOS8gpKFIt8/N0xeyBIejbwq9K35uIbINJkllrxWSWyHJpS2RLMaH9hwWX1tpzMQdvJZ1F2X9YSu/Jxo9oy4SWiCpUpQ+AyeVyJCYm4sCBA8jLy1O1uS118OBBQ3dJRIS/i+U6E1kA2H8pD38Xy21uyoFeBAE4ckQ5lWDbNqDM371qfHyASZOACROUfzYRuUJA7I5L5RJZABCgTGhjd1xCrxBfTjkgItEYnMxOnz4diYmJ6N+/P1q0aAEJHxggIhF8suuS3uPmDm5ZxdFYECsqrXUy857a1IKyBAA5BUU4mXkP4Q29TBcYEdk0g5PZjRs34vvvv0dERERVxENEdiorX7+uV/qOs3pWWFor74H2RNaYcURE+jA4mXV2dkajRo2qIhYismOBXu44clW/cTbt3Dlg8WJgwwaguFj7OKkUGDcOmDIFaNDAdPHp4C11FXUcEZE+DK4D884772Dx4sWws+fGiKiKva9nfVV9x1kVuRzYvh148UXlVIGEBO2JbMOGymT31i3giy8sJpEFgOeCasPP0xXa7g1LoKxq8FxQbVOGRUQ2zuA7s7/88gsOHTqE3bt3o3nz5qpas6U2b94sWnBEZD/cnB3RK8S7wmoGNvXwl4WW1jKWo4MEsweG4K2ks5AAag+ClSa4sweG8OEvIhKVwclszZo1MWTIkKqIhYjs3MqR7e2jzqwFl9aqrL4t/BA/om25OrO+rDNLRFWEdWaJyOLYZAcwKyitJSZ2ACOiyqjSOrOAsn1tcnIyrl27huHDh0MqleKPP/6Ah4cHatSoYVTQRESlnJ0cENGyrioRMrrNq8iMStBkMmDjRuU8VwsvrSUmRwcJy28RkUkYnMzeuHEDffv2RXZ2NmQyGXr16gWpVIpPP/0URUVFWLZsWVXESUR2wlJboRoc1507ytJa8fE6S2spIMG+Jh2wreurGDR1GPq2rFsV4RMR2SyDb3dMnz4dYWFh+Ouvv+Dm5qZaPmTIEBw4cEDU4IjIvpS2Qi1beD+3oAhvJZ3Fnos5lh/XuXPAmDFAQAAwZ47WRPaBsxtWhQ1C1wkrMXHIB9hTuwne2pBqtmMkIrJWRlUzOHr0KJydndWW169fH7dv3xYtMCKyL5baClWfuOZuu4Bev5+E45LFwKFDOvd3q3ZdrAodgB9b9sRDl//VzGW7VyIi4xiczCoUCsjl8nLLb926BalUKkpQRGR/LLUVqq64qsse49ULPyPqzA443q/gjmr37rg8bAwirnlA4aD5YTa2eyUiMpzByWyvXr2waNEirFixAgAgkUjw8OFDzJ49my1uicholtoKVdP71bufi6gzO/Da+f3wKNbRXrdMaa2MtNtQZKYZ9Z5ERKSZwcnswoUL0b17d4SEhKCoqAjDhw/H1atX8cwzz+Dbb7+tihiJyA5YaitU1fsJAp6/eRGjz2xHr6sn4CgYXlrLUo+RiMiaGZzM1q1bF2lpafj2229x9uxZKBQKjB07FpGRkWoPhBERGaK0FWpuQZHG+akSKAvvm7oV6nN1q2PstcN4+fCPaJ53Xffg0FBlly4tpbUs9RiJiKwZmyYQkcUorRoAaG6FGj+irenKcxlQWiuvR1/4znoP6NwZkOh+cMuijpGIyEIZkq/plcxu375d7zd/6aWX9B5rDkxmiSyb2evMnjunbHCwYQNQXKx12ANnN/zUPgJ1P5iJrv06GPQWZj9GIiILJ3oy6+CgXzlaiUSisdKBJWEyS4awh5acYh2jmOfK5OddLgd27lS2mq2gtFZRQBCuDovC38NHol2r+tZzjEREVkT0drYKXT3EiWyUPdw9E+sYxT5XJmuF+uABkJAALFkCXLume2z37sDbb8O1f3+0dNRcWssQbPdKRCQOzpkl0qB0XmPZXw5bmtco1jFa5bm6fh1YuhRYswYoLNQ+rkxpLSIiMg1D8jWD29kS2bqKOj4Byi5NcoX1/j9QrGO0qnMlCEBKCjBkCNCokXJKgbZE1scHiI0FsrOVCS8TWSIii8VklqgMQzpRWSuxjtEqzpVMBqxdC7RtC3TrBmzdqkxsNQkNVY69cQOYNUutRiwREVkmg+vMEtk6S+1EJSaxjtGiz5WepbUgkQCDByvrw+pRWouIiCwLk1miMuyhS5NYx2iR50rP0lqQSoFx44ApU4AGDUwXHxERiUqvZLZQ1wMSZfChKrJ29tClSaxjtJhzZUBpLTRsCEybBkRFAfz7iojI6uk1Z7ZmzZqoVauWzp/SMUTWztFBgtkDQwD874n8UqWvZw8MseqaoGIdo9nPVWGh8i5skybAoEG6E9nu3YFt24CMDGUyy0SWiMgm6FWaKyUlRe8ddu3atVIBVTWW5iJ9sc6s+erMVqi0tNbq1cpasdqwtBYRkVUSvQOYLWEyS4awhy5NYh3j38VyfLLrErLyHyPQyx3vR4TAzdm45gLFJQqsP5aFG/ceo35td7wZHghnRwlw+LByKsG2bdorEgDKKgSTJgETJ0L+TB2b/gzt4RolIvtjkmT28ePHyM7ORnGZByxatWql9z7i4+MRHx+PrKwsAEDz5s0xa9Ys9OvXT+P45ORkdO/evdzy9PR0NG3aVK/3ZDJLJL64XZew8kgmni4n6yABxncOQkxESKX25VzyBC9dPoyZ6bvhe/2y7o1DQ5VVCYYNA1xcbP7uuq0fHxHZL9Hb2T7tzz//xOjRo7F7926N6+Vyud77qlevHubPn49GjRoBANauXYtBgwYhNTUVzZs317pdRkaG2oHVqVNH7/ckInHF7bqE5Yczyy1XCFAt1zehfXpfzzz6C5GpuzEibRfqPLqvfSMtpbW0dSbLLSjCW0lnLbMzmQFs/fiIiPRlcNOEt99+G3/99ReOHz8ONzc37NmzB2vXrkXjxo2xfft2g/Y1cOBAREREoEmTJmjSpAk+/vhj1KhRA8ePH9e5nbe3N3x9fVU/jiL0SSciwxWXKLDySPlE9mkrj2SiuESh975C7lzHZzsX4Wj8aMw4+o32RFYqBWbMAH7/Hdi8GejSRZXIWlVnMiPY+vERERnC4DuzBw8exLZt29C+fXs4ODigfv366NWrFzw8PBAXF4f+/fsbFYhcLscPP/yAR48eITw8XOfY0NBQFBUVISQkBB9++KHGqQelZDIZZDKZ6rUhZcaISLf1x7JQUb6kEJTjxnbWUctVLkfKp6uw4ZuvEJ59Qef+Cp8NgMf/vaOztJYhncnCG3rpPgALZOvHR0RkCIOT2UePHsHb2xsAULt2bfz5559o0qQJWrZsibNnzxocwIULFxAeHo6ioiLUqFEDW7ZsQUiI5q8k/fz8sGLFCrRr1w4ymQzr169Hjx49kJycjC5dumjcJi4uDrGxsQbHRUQVu3HvceXGFRYCCQnAkiXodf26zn38GtAKa8IG4dnIoYh9WXdlAovuTCYCWz8+IiJDGJzMBgcHIyMjA4GBgWjTpg2WL1+OwMBALFu2DH5+hs/PCg4ORlpaGu7fv49NmzZh1KhRSElJ0ZjQBgcHIzg4WPU6PDwcN2/exOeff641mY2JiUF0dLTqdWFhIfz9/Q2Ok4jKq1/b3bhxepbWkjlWw9aQbkgIewmXvYMAAP+pI63w/SyyM5mIbP34iIgMYXAy+/bbbyMnJwcAMHv2bPTp0wcbNmyAs7MzEhMTDQ7A2dlZ9QBYWFgYTp06hcWLF2P58uV6bd+hQwckJSVpXe/i4gIXFxeD4yKiir0ZHoiPd6XrnGrgIFGOgyDoXVrrz+o1sT60Pza06Yf86jXL76sCFtOZrIrY+vERERnC4GQ2MjJS9efQ0FBkZWXh8uXLCAgIwDPPPFPpgARBUJvjWpHU1FSj7ggTUeU5OzlgfOcgjdUMSk3s8CycN6xXJrFpaTr3l9uwGT5t2hc/Ne2CYqdq5daP7xwEZ6eKn1st7Uz2VtJZSAC1hM8WurjZ+vERERnC4GoGH330ER4//t/8N3d3d7Rt2xbVq1fHRx99ZNC+3n//fRw5cgRZWVm4cOECPvjgAyQnJ6sS5piYGIwcOVI1ftGiRdi6dSuuXr2K3377DTExMdi0aROmTJli6GEQkUhiIkIwoUsQyuZN3o//wvpbu/F/E/oqH9bSlshKJMCQIUBKCnyv/oY6k8ahpJp6IusgASZ0Maxmbd8Wfogf0Ra+nupftft6utpE2SpbPz4iIn0Z3DTB0dEROTk5qofASuXn58Pb29ugOrNjx47FgQMHkJOTA09PT7Rq1QrvvvsuevXqBQCIiopCVlYWkpOTAQCffvopVqxYgdu3b8PNzQ3NmzdHTEwMIiIi9H5PNk0gQ1hqdyWNHbL0uGNZlR4WlWDGd6lwungeI45vRceT+yAp01RFjVQKjBsHTJkCNFCvdFDl3cTMfK7EjMtSr1Eiosqo0g5gDg4OuHPnTrlGBQcPHsSwYcPw559/Gh6xCTGZJX1ZanclMbttiWX+jgu4nvg9Rp/eVmFpLTRsCEybprW0lpjnnZ8hEZF1qpJktlatWpBIJKqdSiT/+5+/XC7Hw4cPMXHiRHz11VeVi76KMZklfWjrrlR61Zvra1xt3bZKGfpVfKUVFmL/v+PQ5PtE1L+fq3vsiy8qu3RFRABaGp2Ied75GRIRWa8qSWbXrl0LQRAwZswYLFq0CJ6enqp1zs7OCAwMrLDZgSVgMksVkSsEvLDgoNai9KVPiv/y7osm/Tq3uESBpv/ZXWHlgMtz+1X91+j/lNYSVq+GpILSWtuad8OQhE9RrW0bnbsU87zzMyQism6G5Gt6VzMYNWoUACAoKAidOnWCk5PBhRCIrIKldlcSrduWsTSU1tKWBpYtrfXgkQfGVrB7Mc87P0MiIvthcEbatWtXXLt2DQkJCbh27RoWL14Mb29v7NmzB/7+/mjevHlVxElkMpbaXanS3baMJZMBGzfqVVrrok9DrAl7qVxpLX1iEvO88zMkIrIfBn+PlZKSgpYtW+LEiRPYvHkzHj58CAA4f/48Zs+eLXqARKZmqd2VjO62Zaw7d4DYWCAgQGdpLQUk2NMkHK8Nn48BoxZhc4se5WrE6hOTmOednyERkf0wOJl97733MG/ePOzfvx/Ozs6q5d27d8exY8dEDY7IHEq7K2n7Cl0C5RPxpu6u9GZ4YLlarmXp2yFLp7Q0YPRoZRI7Zw6Ql6d5nIcH5NPfRveJKzFxyAc46d9CWTPWyJjEPO92/xkSEdkRg5PZCxcuYMiQIeWW16lTB/n5+aIERWROpd2VAJRLhszZXam025Yu+nbIKkcuV86D7d4dCA0FEhMBbTViGzYEliwBbt2C46KF6DtQ94OfhnbtAip/3u3yMyQislMG/41Zs2ZN5OTklFuempqKZ599VpSgiMzNUrsraeu2ZUyHLABAYSGweDHQpAkweDDwT4MSjV58Edi+HcjIAKZOVTY9EDkmMc+73XyGRER2zuCmCf/3f/+HY8eO4YcffkCTJk1w9uxZ3LlzByNHjsTIkSMtft4sS3ORISy1u1Klu0f9U1oLq1cDOkprwcUFiIwEpk8HWrWq2pieIuZ5t9nPkIjIhlVpB7AnT54gKioKGzduhCAIcHJyglwux/Dhw5GYmAhHLcXQLQWTWbJbGkpraeXjA0yaBEycCJRpXa2NPSSgRERkGlWazJa6du0aUlNToVAoEBoaisaNGxsVrKkxmSW7Y0BpLYSGKrt0DRumvCurJ3toQUtERKZjkmQWAEo3lWh4gtlSMZklu3HnDrBsGfD119orEgCAg4Nyvuz06UDnzhorEuhiDy1oiYjItAzJ14yaoLV69Wq0aNECrq6ucHV1RYsWLbBq1SqjgiUikRlQWgvR0cDvvwObNgFduhicyMoVAmJ3XCqXfAJQLYvdcQnyitpeibwvIiKyHwZ3APvPf/6DhQsXYurUqQgPV5bkOXbsGGbMmIGsrCzMmzdP9CCJqAJyOfDTT8qpBLoqEgDK0lrTpysbIfxTkcBY9tCCloiILJvByWx8fDxWrlyJN954Q7XspZdeQqtWrTB16lQms0SmVFgIJCQo675ev6577IsvKufDRkQAIj2oaQ8taImIyLIZnMzK5XKEhYWVW96uXTuUlJSIEhQRVaAKSmsZwx5a0BIRkWUzeM7siBEjEB8fX275ihUrEBkZKUpQRKSBIAApKcCQIUCjRsopBdoSWR8fIDYWyM5WJrxVkMgC9tGCloiILJvBd2YB5QNg+/btQ4cOHQAAx48fx82bNzFy5EhER0erxn3xxRfiRElkz0xQWstYpW1j30o6Cwmg9vCWsS1oxdgXERHZD4NLc3Xv3l2/HUskOHjwoFFBVSWW5iKrYWhprbffBl54weCKBGJgnVkiIhKTyerMWiMms5bLErs+iRmT3u1L09KAxYuBb74Biou179DDAxg3DpgyBQgKMiomMVlqO1siIrI+TGZ1YDJrmSzxbpyYMcXtuoSVRzLxdIlUBwkwvnMQYiJCzFZaSyyW+PkREZH1YjKrA5NZy2OJXZ/EjClu1yUsP5ypcV0N2WMsfHwWvX7+3iyltcRgiZ8fERFZN0PyNaMeACMSS0VdnyRQdn3qFeJrsq+ZxYypuESBlUfKJ7L+93MRdWYHXju/D9Liv7XvoIpLa1WWJX5+RERkX5jMkllZYtcnMWNafyzrf1MLBAHP37yIMae3odfVE3DQmAL+w8cHmDwZmDAB8PY2/CBMxBI/PyIisi9MZsmsLLHrk5gx3bj3GM4lTzAw/TDGnN6G5nkVTCUIDQVmzABee80kpbUqyxI/PyIisi9MZsmsLLHrk2gx3bmDwdtWYep3a1Hn8X2tw+QSB2R37oWgeR+YrbSWsSzx8yMiIvvCZJbMqrTrU25BkcYv3SUAfE3c9anSMT1VWqutjtJahc7u+K51b6xvNwA/Lx0NGFnGypws8fMjIiL7Yn3/epJNKe36BKBcG1NzdX0yKia5HNi2DejeXTlVIDFRa43YrJp+mN1zAsInJeLjF8eh38Bwo+uxmpslfn5ERGRfWJqLLIIl1inVK6bCQiAhAViypMLSWr/Wb4XVYYNwqEEYFA6O6nVmrZwlfn5ERGS9WGdWByazlssSuz5p7Wp1/TqwdCmwejXw4IH2HTxVWqs4pIVFdsgSq3OXJX5+RERknZjM6sBklvRV7m6jICDi3hXMvr4fPof2Arp+daqwtJZJO5MRERGZAZsmEFXS012tLKm0lrZuW7kFRXgr6awonckUAlTLmdASEZGlYzJLVEZpVyuvR38hMnU3RqTu0llaCw4OwODBylazVVhayxSdyZ628kgm3und1GofTiMiIvvAZJaojIs7UxD9zXy8lJ4MF3mJ1nElNaRw+td4YMoUICioyuOqss5kWigE5bixnRsYES0REZFpMJklApSltX76CVi0CK2Tk9Fax9Csmn5ICHsJz82Zgf6dgk0WotidyfSh7zgiIiJzYTJL9s2A0lpH67fCmqdKa/X1fcZEQSqJ2W2rfm13vfal7zgiIiJzYTJL9knP0loyx2rYGtINCWEv4bK3ciqBBMrqAabuaiVmt603wwPx8a50nVMNHCTKcURERJaMySzZD0EADh8GFi1SduvSUVpL5lUHXzXrjW/a9MPd6jVVy83Z1aq029ZbSWchAdQSWkPjcnZywPjOQRqrGZQa3zmID38REZHF479UZPtkMmDtWqBtW6BbN2DrVu2JbGgosG4dXG7fREj8Z6hW11dtta+nq0Hlr8TWt4Uf4ke0ha+n+lQCY+KKiQjBhC5BKJv7OkiACV1YZ5aIiKwDmyaQ7bpzB1i2DPj6ayAvT/s4HaW1LLWrlSV2ACMiIhKLIfmaWf/Fio+PR6tWreDh4QEPDw+Eh4dj9+7dOrdJSUlBu3bt4OrqigYNGmDZsmUmipY0kSsEHLuWj21pt3HsWj7kFdV7MsW+0tKA0aOBgABgzhztiayHBxAdDfz+O7BpE9C5c7kasY4OEoQ39MKgNs8ivKGXRSSygPJcXfqjAGdu/IVLfxRU6rw7OkgQUtcT7erXQkhdT6OPUcxrgYiISF9mnTNbr149zJ8/H40aNQIArF27FoMGDUJqaiqaN29ebnxmZiYiIiIwfvx4JCUl4ejRo5g0aRLq1KmDoUOHmjp8uydmW9VK7+up0lpITtY9tmFDYPp0ICoKkEoNitMSaGpB+/GudKNa0Ir1GYp5LRARERnC4qYZ1K5dG5999hnGjh1bbt27776L7du3Iz09XbVs4sSJOHfuHI4dO6bX/jnNQBza2qqW3tMzZP5mpfZlQGktvPiicipBRATg6KhXbJZGWwvaUobMdRXrMxTzWiAiIgKsaJrB0+RyOTZu3IhHjx4hPDxc45hjx46hd+/easv69OmD06dP48mTJ6YIk1BxW1VA2VZVn6+Zjd7X9evAjBlAvXrKBFVbIuviAowZA5w7Bxw4AAwcaLWJrL4taItLFBXuS6zPUMxrgYiIyBhmT2YvXLiAGjVqwMXFBRMnTsSWLVsQEqL5zlJubi58fHzUlvn4+KCkpAR3797VuI1MJkNhYaHaD1WOIW1VRd2XIAApKcCQIUCjRsopBdpqxPr4AB99BGRnK2vJtmpVYSyWzpAWtBUR6zMU81ogIiIyhtnrzAYHByMtLQ3379/Hpk2bMGrUKKSkpGhNaCVlHtApnSVRdnmpuLg4xMbGihu0nROzrao+Y5xLnsBlw3pgy1rlw126hIYq79i+9pryrqwNEbMFrVifoZjXAhERkTHMnsw6OzurHgALCwvDqVOnsHjxYixfvrzcWF9fX+Tm5qoty8vLg5OTE7y8vDTuPyYmBtHR0arXhYWF8Pf3F/EI7I+YbVV1jXnm0V+ITN2NEam7UOfxfe070VFay5aI2YJWrM9QzGuBiIjIGGZPZssSBAEymUzjuvDwcOzYsUNt2b59+xAWFoZq1app3MbFxQUuNnaHztzEbKuqaV8hd65j9OnteCk9GS7yEu0be3gA48YBU6YAQUHGHIpVEbMFrVifoZjXAhERkTHMOmf2/fffx5EjR5CVlYULFy7ggw8+QHJyMiIjIwEo76qOHDlSNX7ixIm4ceMGoqOjkZ6ejjVr1mD16tWYOXOmuQ7BLpW2VQX+98R6KUPbqpbuy0EhR++rx/HttzHYlTgNr178WXsi27ChsnrBrVvAf/9rF4ks8L8WtLro24JWrM9QzGuBiIjIGGZNZu/cuYM333wTwcHB6NGjB06cOIE9e/agV69eAICcnBxkZ2erxgcFBWHXrl1ITk5GmzZtMHfuXCxZsoQ1Zs1AtLaqhYXoe+B7XPh2GlZsnofw7Avax774IrB9O5CRAUydapU1YitLzBa0Yn2GYrbYJSIiMpTF1ZmtaqwzKy6j26pevw4sXaqsNKCtIgGgfIgrMlLZ5MAGKhKIRcwWtGK1xrXU1r9ERGR9DMnXmMyS6QgCcPiwsqTWtm3K11r8Wb0mtoQPRqMPo/Fi15ami5GIiIjMzpB8zeIeACMbJJMBGzcqk9gKSmtd9GmI1WGDsLNpZzxxqgbszka81zP8qpqIiIg0YjJLVefOHWDZMuDrr4G8PK3D5BIH7GvcAWvCXsKpes3VSmtJoOwg1SvEl19ZExERUTlMZkl8aWnA4sXAN98AxcXax3l44I9XIvFatXa4VdNX45CnO0iFN9RcS5iIiIjsF5NZEodcDvz0k3IqQXKy7rENGyof6IqKwqlrhbi1Ma3C3bODFBEREWnCZJYqp7AQSEhQ1n29fl332BdfVHbpiogAHB0BAN5SHXdun8IOUkRERKQJk1kyjkiltdhBioiIiCrDrE0TyMoIgnIKweDBQKNGyikF2hJZX1/go4+A7GxlwqulRiw7SBEREVFlMJmlislkwNq1QNu2QPfuumvEhoYC69YBWVnAf/4DeHtXuHt2kCIiIiJjcZoBaadnaS04OCjv1r79NvDCC2qltfTVt4UfeoX4soMUERERGYTJLJVnQGktjBsHTJkCBAVV+m0dHSQsv0VEREQGYTJLSkaW1oJUaoLgiIiIiDRjMmvvCguBNWuUlQmMKK1FREREZE5MZu3VtWvKBHbNmkqV1iIiIiIyJyaz9kQQgJQU5VSC7du1VyQAlKW1Jk0CJkzQqyIBERERkTkwmbUHMhmwcaMyiU1L0z02NBSYMQN47TXlXVkiIiIiC8Zk1paZsLQWERERkTkwmbVFZiqtRURERGRqTGZthVwO7NihTGJZWouIiIjsBJNZa1daWmvJEiAzU/dYltYiIiIiG8Nk1lqxtBYRERERk1mrwtJaRERERGqYzFoDltYiIiIi0ojJrCVjaS0iIiIinZjMWqK0NOVd2G+/ZWktIiIiIh2YzFqK0tJaixYp58XqwtJaRERERACYzJqfMaW1+vdXTi0gIiIisnNMZs2FpbWIiIiIKo3JrCmxtBYRERGRqJjMmgJLaxERERFVCSazVenOHSA+XvnD0lpEREREomMyW1X+/hsIDgYKCrSPYWktIiIiokrhI/FVxc0NePVVzesaNVI+/HXrFvDf/zKRJSIiIjISk9mqNH26+usXX1Q++JWRobwbyxqxRERERJXCaQZVqUULYMAAZTUCltYiIiIiEh2T2aq2fTsf6CIiIiKqIpxmUNWYyBIRERFVGSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLbMms3FxcWjfvj2kUim8vb0xePBgZGRk6NwmOTkZEomk3M/ly5dNFDURERERWQqzJrMpKSmYPHkyjh8/jv3796OkpAS9e/fGo0ePKtw2IyMDOTk5qp/GjRubIGIiIiIisiRmrTO7Z88etdcJCQnw9vbGmTNn0KVLF53bent7o2bNmlUYHRERERFZOouaM1tQUAAAqF27doVjQ0ND4efnhx49euDQoUNax8lkMhQWFqr9EBEREZFtsJhkVhAEREdH44UXXkCLFi20jvPz88OKFSuwadMmbN68GcHBwejRowcOHz6scXxcXBw8PT1VP/7+/lV1CERERERkYhJBEARzBwEAkydPxs6dO/HLL7+gXr16Bm07cOBASCQSbN++vdw6mUwGmUymel1YWAh/f38UFBTAw8Oj0nETERERkbgKCwvh6empV75mEXdmp06diu3bt+PQoUMGJ7IA0KFDB1y9elXjOhcXF3h4eKj9EBEREZFtMOsDYIIgYOrUqdiyZQuSk5MRFBRk1H5SU1Ph5+cncnREREREZOnMmsxOnjwZ33zzDbZt2wapVIrc3FwAgKenJ9zc3AAAMTExuH37NtatWwcAWLRoEQIDA9G8eXMUFxcjKSkJmzZtwqZNm8x2HERERERkHmZNZuPj4wEA3bp1U1uekJCAqKgoAEBOTg6ys7NV64qLizFz5kzcvn0bbm5uaN68OXbu3ImIiAhThU1EREREFsJiHgAzFUMmFBMRERGR6VndA2BERERERMZgMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdUyaztboqoiVwg4mXkPeQ+K4C11xXNBteHoIDF3WERERCQyJrNkc/ZczEHsjkvIKShSLfPzdMXsgSHo28LPjJERERGR2DjNgGzKnos5eCvprFoiCwC5BUV4K+ks9lzMMVNkREREVBWYzJLNkCsExO64BEHDutJlsTsuQa7QNIKIiIisEZNZshknM++VuyP7NAFATkERTmbeM11QREREVKWYzJLNyHugPZE1ZhwRERFZPiazZDO8pa6ijiMiIiLLx2SWbMZzQbXh5+kKbQW4JFBWNXguqLYpwyIiIqIqxGSWbIajgwSzB4YAQLmEtvT17IEhrDdLRERkQ5jMkk3p28IP8SPawtdTfSqBr6cr4ke0ZZ1ZIiIiG8OmCWRz+rbwQ68QX3YAIyIisgNMZskmOTpIEN7Qy9xhEBERURXjNAMiIiIislpMZomIiIjIajGZJSIiIiKrxWSWiIiIiKwWk1kiIiIislpMZomIiIjIatldaS5BEAAAhYWFZo6EiIiIiDQpzdNK8zZd7C6ZffDgAQDA39/fzJEQERERkS4PHjyAp6enzjESQZ+U14YoFAr88ccfkEqlkEhM0xGqsLAQ/v7+uHnzJjw8PEzynsTzbg485+bB824ePO/mwfNuHqY+74Ig4MGDB6hbty4cHHTPirW7O7MODg6oV6+eWd7bw8ODv3hmwPNuejzn5sHzbh487+bB824epjzvFd2RLcUHwIiIiIjIajGZJSIiIiKrxWTWBFxcXDB79my4uLiYOxS7wvNuejzn5sHzbh487+bB824elnze7e4BMCIiIiKyHbwzS0RERERWi8ksEREREVktJrNEREREZLWYzBIRERGR1WIyK6K4uDhIJBK8/fbbOselpKSgXbt2cHV1RYMGDbBs2TLTBGij9DnvycnJkEgk5X4uX75sukCt3Jw5c8qdP19fX53b8FqvPEPPO6918dy+fRsjRoyAl5cX3N3d0aZNG5w5c0bnNrzmK8/Q885rvvICAwM1nsPJkydr3caSrnW76wBWVU6dOoUVK1agVatWOsdlZmYiIiIC48ePR1JSEo4ePYpJkyahTp06GDp0qImitR36nvdSGRkZap1L6tSpU1Wh2aTmzZvj559/Vr12dHTUOpbXungMOe+leK1Xzl9//YVOnTqhe/fu2L17N7y9vXHt2jXUrFlT6za85ivPmPNeite88U6dOgW5XK56ffHiRfTq1QuvvvqqxvGWdq0zmRXBw4cPERkZiZUrV2LevHk6xy5btgwBAQFYtGgRAKBZs2Y4ffo0Pv/8c/5lZyBDznspb29vvf5SJM2cnJwqvBtbite6eAw576V4rVfOggUL4O/vj4SEBNWywMBAndvwmq88Y857KV7zxiub+M+fPx8NGzZE165dNY63tGud0wxEMHnyZPTv3x89e/ascOyxY8fQu3dvtWV9+vTB6dOn8eTJk6oK0SYZct5LhYaGws/PDz169MChQ4eqMDrbdPXqVdStWxdBQUF4/fXXcf36da1jea2Lx5DzXorXeuVs374dYWFhePXVV+Ht7Y3Q0FCsXLlS5za85ivPmPNeite8OIqLi5GUlIQxY8ZAIpFoHGNp1zqT2UrauHEjzp49i7i4OL3G5+bmwsfHR22Zj48PSkpKcPfu3aoI0SYZet79/PywYsUKbNq0CZs3b0ZwcDB69OiBw4cPV3GktuP555/HunXrsHfvXqxcuRK5ubno2LEj8vPzNY7ntS4OQ887r3VxXL9+HfHx8WjcuDH27t2LiRMnYtq0aVi3bp3WbXjNV54x553XvLi2bt2K+/fvIyoqSusYS7vWOc2gEm7evInp06dj3759cHV11Xu7sv/TKW3Cpu1/QKTOmPMeHByM4OBg1evw8HDcvHkTn3/+Obp06VJVodqUfv36qf7csmVLhIeHo2HDhli7di2io6M1bsNrvfIMPe+81sWhUCgQFhaGTz75BIDyrt9vv/2G+Ph4jBw5Uut2vOYrx5jzzmteXKtXr0a/fv1Qt25dneMs6VrnndlKOHPmDPLy8tCuXTs4OTnByckJKSkpWLJkCZycnNQmU5fy9fVFbm6u2rK8vDw4OTnBy8vLVKFbNWPOuyYdOnTA1atXqzha21W9enW0bNlS6znktV41KjrvmvBaN5yfnx9CQkLUljVr1gzZ2dlat+E1X3nGnHdNeM0b58aNG/j5558xbtw4neMs7VrnndlK6NGjBy5cuKC2bPTo0WjatCneffddjU8ch4eHY8eOHWrL9u3bh7CwMFSrVq1K47UVxpx3TVJTU+Hn51cVIdoFmUyG9PR0dO7cWeN6XutVo6LzrgmvdcN16tQJGRkZasuuXLmC+vXra92G13zlGXPeNeE1b5yEhAR4e3ujf//+OsdZ3LUukKi6du0qTJ8+XfX6vffeE958803V6+vXrwvu7u7CjBkzhEuXLgmrV68WqlWrJvz4449miNZ2VHTeFy5cKGzZskW4cuWKcPHiReG9994TAAibNm0yQ7TW6Z133hGSk5OF69evC8ePHxcGDBggSKVSISsrSxAEXutVxdDzzmtdHCdPnhScnJyEjz/+WLh69aqwYcMGwd3dXUhKSlKN4TUvPmPOO695ccjlciEgIEB49913y62z9GudyazIyiZVo0aNErp27ao2Jjk5WQgNDRWcnZ2FwMBAIT4+3rRB2qCKzvuCBQuEhg0bCq6urkKtWrWEF154Qdi5c6fpA7Viw4YNE/z8/IRq1aoJdevWFV5++WXht99+U63ntV41DD3vvNbFs2PHDqFFixaCi4uL0LRpU2HFihVq63nNVw1DzzuveXHs3btXACBkZGSUW2fp17pEEP6ZsUtEREREZGX4ABgRERERWS0ms0RERERktZjMEhEREZHVYjJLRERERFaLySwRERERWS0ms0RERERktZjMEhEREZHVYjJLRGQFoqKiMHjwYK3rExMTUbNmTZPFU5HAwEAsWrTI3GEQkR1gMktEREaztCSaiOwPk1kiIiIislpMZomIKvDjjz+iZcuWcHNzg5eXF3r27IlHjx6p1ickJKBZs2ZwdXVF06ZN8fXXX6vWZWVlQSKRYOPGjejYsSNcXV3RvHlzJCcnq8bI5XKMHTsWQUFBcHNzQ3BwMBYvXlzpuHfs2IF27drB1dUVDRo0QGxsLEpKSlTrJRIJVq1ahSFDhsDd3R2NGzfG9u3b1faxfft2NG7cGG5ubujevTvWrl0LiUSC+/fvIzk5GaNHj0ZBQQEkEgkkEgnmzJmj2vbx48cYM2YMpFIpAgICsGLFikofExFROQIREWn1xx9/CE5OTsIXX3whZGZmCufPnxe++uor4cGDB4IgCMKKFSsEPz8/YdOmTcL169eFTZs2CbVr1xYSExMFQRCEzMxMAYBQr1494ccffxQuXbokjBs3TpBKpcLdu3cFQRCE4uJiYdasWcLJkyeF69evC0lJSYK7u7vw3XffqeIYNWqUMGjQIK1xJiQkCJ6enqrXe/bsETw8PITExETh2rVrwr59+4TAwEBhzpw5qjGlcX3zzTfC1atXhWnTpgk1atQQ8vPzVbFXq1ZNmDlzpnD58mXh22+/FZ599lkBgPDXX38JMplMWLRokeDh4SHk5OQIOTk5qvNSv359oXbt2sJXX30lXL16VYiLixMcHByE9PR0UT4XIqJSTGaJiHQ4c+aMAEDIysrSuN7f31/45ptv1JbNnTtXCA8PFwThf8ns/PnzVeufPHki1KtXT1iwYIHW9500aZIwdOhQ1WtDk9nOnTsLn3zyidqY9evXC35+fqrXAIQPP/xQ9frhw4eCRCIRdu/eLQiCILz77rtCixYt1PbxwQcfqJJZTe9bqn79+sKIESNUrxUKheDt7S3Ex8drPQYiImM4mfGmMBGRxWvdujV69OiBli1bok+fPujduzdeeeUV1KpVC3/++Sdu3ryJsWPHYvz48aptSkpK4Onpqbaf8PBw1Z+dnJwQFhaG9PR01bJly5Zh1apVuHHjBv7++28UFxejTZs2Rsd95swZnDp1Ch9//LFqmVwuR1FRER4/fgx3d3cAQKtWrVTrq1evDqlUiry8PABARkYG2rdvr7bf5557Tu8Ynt63RCKBr6+vat9ERGJhMktEpIOjoyP279+PX3/9Ffv27cPSpUvxwQcf4MSJE6qEcOXKlXj++efLbVcRiUQCAPj+++8xY8YM/Pe//0V4eDikUik+++wznDhxwui4FQoFYmNj8fLLL5db5+rqqvpztWrVysWkUCgAAIIgqGIsJQiC3jHo2jcRkViYzBIRVUAikaBTp07o1KkTZs2ahfr162PLli2Ijo7Gs88+i+vXryMyMlLnPo4fP44uXboAUN65PXPmDKZMmQIAOHLkCDp27IhJkyapxl+7dq1SMbdt2xYZGRlo1KiR0fto2rQpdu3apbbs9OnTaq+dnZ0hl8uNfg8iospiMktEpMOJEydw4MAB9O7dG97e3jhx4gT+/PNPNGvWDAAwZ84cTJs2DR4eHujXrx9kMhlOnz6Nv/76C9HR0ar9fPXVV2jcuDGaNWuGhQsX4q+//sKYMWMAAI0aNcK6deuwd+9eBAUFYf369Th16hSCgoKMjnvWrFkYMGAA/P398eqrr8LBwQHnz5/HhQsXMG/ePL32MWHCBHzxxRd49913MXbsWKSlpSExMRHA/+4qBwYG4uHDhzhw4ABat24Nd3d31R1rIiJTYGkuIiIdPDw8cPjwYURERKBJkyb48MMP8d///hf9+vUDAIwbNw6rVq1CYmIiWrZsia5duyIxMbFcIjp//nwsWLAArVu3xpEjR7Bt2zY888wzAICJEyfi5ZdfxrBhw/D8888jPz9f7S6tMfr06YOffvoJ+/fvR/v27dGhQwd88cUXqF+/vt77CAoKwo8//ojNmzejVatWiI+PxwcffAAAcHFxAQB07NgREydOxLBhw1CnTh18+umnlYqbiMhQEsGQCVBERGSQrKwsBAUFITU1tVIPdFmKjz/+GMuWLcPNmzfNHQoREQBOMyAiIh2+/vprtG/fHl5eXjh69Cg+++wz1VxfIiJLwGSWiIi0unr1KubNm4d79+4hICAA77zzDmJiYswdFhGRCqcZEBEREZHV4gNgRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1mMwSERERkdViMktEREREVovJLBERERFZLSazRERERGS1/h8krms46M4PUAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1, figsize=(8, 4))\n",
    "plt.clf()\n",
    "plt.scatter(df[df[\"species\"] == \"setosa\"].sepal_length, df[df[\"species\"] == \"setosa\"].sepal_width,label=\"Iris Data\")\n",
    "plt.plot(X_plot, fit, label=\"Linear Regression Model\", color=\"red\", linewidth=3)\n",
    "plt.xlabel(\"sepal length\")\n",
    "plt.ylabel(\"sepal width\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "53736520",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmean = np.mean(X_setosa_sepal_length.squeeze())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4648c47d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Standard feil beta1 = 0.010809\n"
     ]
    }
   ],
   "source": [
    "rse = np.sqrt(model.ssr/(48))\n",
    "se_beta1 = (rse*rse)/(np.sum((X_setosa_sepal_length.squeeze()-xmean)*(X_setosa_sepal_length.squeeze()-xmean)))\n",
    "print(\"Standard feil beta1 = %.6f\"%se_beta1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61099238",
   "metadata": {},
   "source": [
    "Oppgaven fra forelesning 15. januar stopper her. Vi har ikke gått gjennom hypotesetesting ennå. Koden under er derfor å anse som \"work in progress\".\n",
    "\n",
    "# Hypotesetesting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "16de5968",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "73.8780703215139\n"
     ]
    }
   ],
   "source": [
    "tval = (model.params[1]-0)/se_beta1\n",
    "print(tval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ca3f460f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff086f8f610>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApeElEQVR4nO3df3DUdX7H8dcmJJsIZBUCCYEQoh6WmhaHoJhw1IJHuKioraOoMwIeUHOCHIJXZegdSnGiDnLUH+FHRSyVYnqWc3Amp4YbBQSuhRz2EDxFiSRAYho8k4iSkOTTPwg72c3ukg3J95PNPh8zO2M+3883+/7OR5JXPp/P97suY4wRAACAJTG2CwAAANGNMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAqn62C+iM1tZWnTp1SgMHDpTL5bJdDgAA6ARjjBoaGpSWlqaYmODzHxERRk6dOqX09HTbZQAAgC6orKzUiBEjgh6PiDAycOBASecvJikpyXI1AACgM+rr65Wenu79PR5MRISRC0szSUlJhBEAACLMxbZYsIEVAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWBV2GNm1a5emT5+utLQ0uVwuvfXWWxc9Z+fOncrOzlZCQoKuvPJKrVu3riu1AgCAPijsx8GfOXNGY8eO1YMPPqi77rrrov3Ly8t1yy23aN68eXr99de1Z88ePfzwwxoyZEinzu92DQ3SY4+F7nMpnwzck+faquti51NXeOf3xbpCne9ySYMHS9deK02aJA0YcGnvA6DPCTuM5OfnKz8/v9P9161bp5EjR2rNmjWSpDFjxujAgQNatWqVnTDSGcbYORfoq77+Wjp6VHrvPWnePOkv/9J2RQB6kR7fM7Jv3z7l5eX5tE2bNk0HDhzQuXPnAp7T2Nio+vp6nxeAPuC776SXX5a+/NJ2JQB6kR4PI9XV1UpJSfFpS0lJUXNzs2prawOeU1hYKI/H432lp6f3dJkAnNLcLP37vzOLCMDLkbtp/D862LT9EAr2kcJLly5VXV2d91VZWdnjNQJw0IkT0uHDtqsA0EuEvWckXKmpqaqurvZpq6mpUb9+/TR48OCA57jdbrnd7p4uDYBNBw5IWVm2qwDQC/R4GMnJydHbb7/t0/bee+9p/PjxiouL6+m378jlki4l6ISaWra58TUS60J0+/jj8/9/XOpdPAAiXthh5Ntvv9Xnn3/u/bq8vFwfffSRBg0apJEjR2rp0qU6efKkNm/eLEkqKCjQSy+9pMWLF2vevHnat2+fNm7cqK1bt3bfVYRjwADphRfsvDfCYytEXexc6ur8uWfPSk88EfhYQ8P5u2yCzJACiB5hh5EDBw5o8uTJ3q8XL14sSZo1a5Zee+01VVVVqaKiwns8MzNTJSUlevTRR/Xyyy8rLS1NL7zwQu+9rRe9h83ncaB7JCZKV10lffFF4OMVFYQRAHIZ0/vn0evr6+XxeFRXV6ekpCTb5QAIx3/+p/S73wU+dscd0i23OFsPAMd09vc3n00DoGelpQU/FuT2fgDRhTACoGclJwc/9n//51wdAHotwgiAnhUqjDAzAkCEEQA9bdCg4BuK6+q4/RsAYQRAD4uJCf5JvS0t0vffO1sPgF6HMAKg54W6C44PwgSiHmEEQM8jjAAIgTACoOcRRgCEQBgB0PMGDgx+rKHBuToA9EqEEQA9L9gGVokNrAAIIwAckJgY/BhhBIh6hBEAPS9UGPnuO+fqANArEUYA9LzLLgt+jJkRIOoRRgD0PGZGAIRAGAHQ85gZARACYQRAz2MDK4AQCCMAeh7LNABCIIwA6Hlud/BP7j171tlaAPQ6hBEAPc/lOh9IAjl3TmptdbYeAL0KYQSAM+Ljgx9ranKuDgC9DmEEgDOCzYxIhBEgyhFGADgjVBhpbHSuDgC9DmEEgDNCLdMQRoCoRhgB4AyWaQAEQRgB4AxmRgAEQRgB4AxmRgAEQRgB4Aw2sAIIgjACwBks0wAIgjACwBks0wAIgjACwBnMjAAIgjACwBnMjAAIgjACwBlxccGPnTvnXB0Aeh3CCABnhAojzc3O1QGg1yGMAHAGMyMAgiCMAHAGYQRAEIQRAM4gjAAIgjACwBn9+gU/RhgBohphBIAz2MAKIAjCCABnsEwDIAjCCABnEEYABEEYAeAM9owACIIwAsAZzIwACIIwAsAZbGAFEARhBIAzWKYBEARhBIAzYmMllyvwMcIIENUIIwCc4XIFnx05d04yxtl6APQahBEAzgm2b8QYqbXV2VoA9BqEEQDO4Y4aAAEQRgA4hzACIADCCADncHsvgAAIIwCcE+r23qYm5+oA0KsQRgA4J1QYaWlxrg4AvQphBIBzQoURlmmAqEUYAeCc2NjgxwgjQNQijABwDss0AAIgjABwDjMjAAIgjABwDjMjAALoUhgpKipSZmamEhISlJ2drd27d4fsv2XLFo0dO1aXXXaZhg0bpgcffFCnT5/uUsEAIhgbWAEEEHYYKS4u1qJFi7Rs2TIdPHhQkyZNUn5+vioqKgL2//DDDzVz5kzNmTNHhw8f1q9//Wvt379fc+fOveTiAUQYlmkABBB2GFm9erXmzJmjuXPnasyYMVqzZo3S09O1du3agP1///vfa9SoUVq4cKEyMzP1wx/+UA899JAOHDhwycUDiDAs0wAIIKww0tTUpLKyMuXl5fm05+Xlae/evQHPyc3N1YkTJ1RSUiJjjL766iu9+eabuvXWW4O+T2Njo+rr631eAPoAZkYABBBWGKmtrVVLS4tSUlJ82lNSUlRdXR3wnNzcXG3ZskUzZsxQfHy8UlNTdfnll+vFF18M+j6FhYXyeDzeV3p6ejhlAuitmBkBEECXNrC6XC6fr40xHdouOHLkiBYuXKhf/vKXKisr0zvvvKPy8nIVFBQE/f5Lly5VXV2d91VZWdmVMgH0NmxgBRBAiJ8MHSUnJys2NrbDLEhNTU2H2ZILCgsLNXHiRP385z+XJP31X/+1+vfvr0mTJmnlypUaNmxYh3Pcbrfcbnc4pQGIBCzTAAggrJmR+Ph4ZWdnq7S01Ke9tLRUubm5Ac/57rvvFBPj+zaxbT+QjDHhvD2ASMcyDYAAwl6mWbx4sV555RW9+uqr+uSTT/Too4+qoqLCu+yydOlSzZw509t/+vTp2rZtm9auXatjx45pz549WrhwoW644QalpaV135UA6P1CzYwQRoCoFdYyjSTNmDFDp0+f1ooVK1RVVaWsrCyVlJQoIyNDklRVVeXzzJHZs2eroaFBL730kpYsWaLLL79cU6ZM0bPPPtt9VwEgMrBnBEAALhMBayX19fXyeDyqq6tTUlKS7XIAdNWuXdKWLYGP5eVJd93lbD0AelRnf3/z2TQAnMMGVgABEEYAOIdlGgABEEYAOIcNrAACIIwAcA4zIwACIIwAcA7PGQEQAGEEgHPYwAogAMIIAOewTAMgAMIIAOewgRVAAIQRAM5hZgRAAIQRAM5hAyuAAAgjAJzDBlYAARBGADiHZRoAARBGADiHDawAAiCMAHAOMyMAAiCMAHAOMyMAAiCMAHAOYQRAAIQRAM4hjAAIgDACwDkxIX7kEEaAqEUYAeAclyt4ICGMAFGLMALAWcGWaggjQNQijABwVrAwYsz5F4CoQxgB4Cw2sQLwQxgB4CzCCAA/hBEAzuKOGgB+CCMAnMXMCAA/hBEAzgoVRlpbnasDQK9BGAHgLGZGAPghjABwFmEEgB/CCABnEUYA+CGMAHAWYQSAH8IIAGexgRWAH8IIAGcxMwLAD2EEgLN46BkAP4QRAM5iZgSAH8IIAGcRRgD4IYwAcBZhBIAfwggAZxFGAPghjABwFhtYAfghjABwFs8ZAeCHMALAWSzTAPBDGAHgLMIIAD+EEQDOIowA8EMYAeAswggAP4QRAM4ijADwQxgB4CzCCAA/hBEAzuI5IwD8EEYAOIuZEQB+CCMAnMVDzwD4IYwAcBYzIwD8EEYAOIswAsAPYQSAswgjAPwQRgA4izACwA9hBICz2MAKwA9hBICzmBkB4IcwAsBZPPQMgB/CCABnMTMCwE+XwkhRUZEyMzOVkJCg7Oxs7d69O2T/xsZGLVu2TBkZGXK73brqqqv06quvdqlgABGOMALAT79wTyguLtaiRYtUVFSkiRMnav369crPz9eRI0c0cuTIgOfcc889+uqrr7Rx40ZdffXVqqmpUXNz8yUXDyACEUYA+Ak7jKxevVpz5szR3LlzJUlr1qzRu+++q7Vr16qwsLBD/3feeUc7d+7UsWPHNGjQIEnSqFGjLq1qAJGLMALAT1jLNE1NTSorK1NeXp5Pe15envbu3RvwnO3bt2v8+PF67rnnNHz4cI0ePVqPPfaYvv/++6Dv09jYqPr6ep8XgD6CDawA/IQ1M1JbW6uWlhalpKT4tKekpKi6ujrgOceOHdOHH36ohIQE/eY3v1Ftba0efvhhff3110H3jRQWFuqpp54KpzQAkYLnjADw06UNrC6Xy+drY0yHtgtaW1vlcrm0ZcsW3XDDDbrlllu0evVqvfbaa0FnR5YuXaq6ujrvq7KysitlAuiNWKYB4CesmZHk5GTFxsZ2mAWpqanpMFtywbBhwzR8+HB5PB5v25gxY2SM0YkTJ/SDH/ygwzlut1tutzuc0gBECsIIAD9hzYzEx8crOztbpaWlPu2lpaXKzc0NeM7EiRN16tQpffvtt962zz77TDExMRoxYkQXSgYQ0QgjAPyEvUyzePFivfLKK3r11Vf1ySef6NFHH1VFRYUKCgoknV9imTlzprf//fffr8GDB+vBBx/UkSNHtGvXLv385z/XT37yEyUmJnbflQCIDIQRAH7CvrV3xowZOn36tFasWKGqqiplZWWppKREGRkZkqSqqipVVFR4+w8YMEClpaV65JFHNH78eA0ePFj33HOPVq5c2X1XASByEEYA+HEZY4ztIi6mvr5eHo9HdXV1SkpKsl0OgEvR0CA99ljgY8nJ0tNPO1sPgB7T2d/ffDYNAGfxnBEAfggjAJzFc0YA+CGMAHAWe0YA+CGMAHAWYQSAH8IIAGcFeVqzJJZpgChFGAHgLJcr+OwIMyNAVCKMAHBesDtqCCNAVCKMAHBesJkRY86/AEQVwggA57GJFUA7hBEAziOMAGiHMALAeYQRAO0QRgA4jzACoB3CCADn8fk0ANohjABwHp9PA6AdwggA57FMA6AdwggA5xFGALRDGAHgPJZpALRDGAHgPGZGALRDGAHgPO6mAdAOYQSA85gZAdAOYQSA8wgjANohjABwXqhlGjawAlGHMALAecyMAGiHMALAeYQRAO0QRgA4jzACoB3CCADnEUYAtEMYAeA8wgiAdggjAJzH3TQA2iGMAHAeMyMA2iGMAHAeYQRAO4QRAM7jU3sBtEMYAeA8ZkYAtEMYAeA8wgiAdggjAJwX6m4awggQdQgjAJzHzAiAdggjAJxHGAHQDmEEgPN46BmAdggjAJzHzAiAdggjAJxHGAHQDmEEgPMIIwDaIYwAcB5hBEA7hBEAziOMAGiHMALAedxNA6AdwggA5zEzAqAdwggA5xFGALRDGAHgvFBhhGUaIOoQRgA4j5kRAO0QRgA4jzACoB3CCADnhbqbhjACRB3CCADnMTMCoB3CCADnEUYAtEMYAeA8HnoGoB3CCADnMTMCoB3CCADnEUYAtEMYAeA8wgiAdggjAJxHGAHQTpfCSFFRkTIzM5WQkKDs7Gzt3r27U+ft2bNH/fr103XXXdeVtwXQVxBGALQTdhgpLi7WokWLtGzZMh08eFCTJk1Sfn6+KioqQp5XV1enmTNn6uabb+5ysQD6CO6mAdBO2GFk9erVmjNnjubOnasxY8ZozZo1Sk9P19q1a0Oe99BDD+n+++9XTk5Ol4sF0EcwMwKgnbDCSFNTk8rKypSXl+fTnpeXp7179wY9b9OmTfriiy+0fPnyTr1PY2Oj6uvrfV4A+pCYGMnlCnyMMAJEnbDCSG1trVpaWpSSkuLTnpKSourq6oDnHD16VE888YS2bNmifv36dep9CgsL5fF4vK/09PRwygQQCYIt1bBMA0SdLm1gdfn9RWOM6dAmSS0tLbr//vv11FNPafTo0Z3+/kuXLlVdXZ33VVlZ2ZUyAfRmwZZqmBkBok7npiraJCcnKzY2tsMsSE1NTYfZEklqaGjQgQMHdPDgQS1YsECS1NraKmOM+vXrp/fee09TpkzpcJ7b7Zbb7Q6nNACRJlgYMeb87EioTa4A+pSw/rXHx8crOztbpaWlPu2lpaXKzc3t0D8pKUmHDh3SRx995H0VFBTommuu0UcffaQJEyZcWvUAIhd31ABoE9bMiCQtXrxYDzzwgMaPH6+cnBxt2LBBFRUVKigokHR+ieXkyZPavHmzYmJilJWV5XP+0KFDlZCQ0KEdQJS52B01ndxjBiDyhf2vfcaMGTp9+rRWrFihqqoqZWVlqaSkRBkZGZKkqqqqiz5zBAC4vRfABS5jjLFdxMXU19fL4/Gorq5OSUlJtssB0B2WLZNqawMfW7VKGjjQ2XoAdLvO/v5mhxgAO5gZAdCGMALADsIIgDaEEQB2hAoj3E0DRBXCCAA7Qt3ay8wIEFUIIwDsYJkGQBvCCAA7WKYB0IYwAsAOZkYAtCGMALCDMAKgDWEEgB2EEQBtCCMA7OBuGgBtCCMA7GBmBEAbwggAO7ibBkAbwggAO1imAdCGMALADpZpALQhjACwgzACoA1hBIAdhBEAbQgjAOxgAyuANoQRAHYwMwKgDWEEgB3cTQOgDWEEgB0s0wBoQxgBYAfLNADaEEYA2EEYAdCGMALADsIIgDaEEQB2EEYAtCGMALCDu2kAtCGMALCDu2kAtCGMALCDZRoAbQgjAOxgmQZAG8IIADuYGQHQhjACwA7CCIA2hBEAdhBGALQhjACwg7tpALQhjACwg5kRAG0IIwDs4G4aAG0IIwDsYJkGQBvCCAA7WKYB0IYwAsAOwgiANoQRAHYQRgC0IYwAsIMwAqANYQSAHdxNA6ANYQSAHdxNA6ANYQSAHSzTAGhDGAFgB8s0ANoQRgDYwcwIgDaEEQB2EEYAtCGMALCDMAKgDWEEgB3cTQOgDWEEgB0u1/lXIMyMAFGFMALAnmCzI4QRIKoQRgDYEyyMsEwDRBXCCAB7goURYwgkQBQhjACwhwefARBhBIBN/foFP0YYAaIGYQSAPaHCyLlzztUBwCrCCAB7QoWR5mbn6gBgFWEEgD2EEQDqYhgpKipSZmamEhISlJ2drd27dwftu23bNk2dOlVDhgxRUlKScnJy9O6773a5YAB9CGEEgLoQRoqLi7Vo0SItW7ZMBw8e1KRJk5Sfn6+KioqA/Xft2qWpU6eqpKREZWVlmjx5sqZPn66DBw9ecvEAIhxhBIAklzHGhHPChAkTNG7cOK1du9bbNmbMGN15550qLCzs1Pe49tprNWPGDP3yl7/sVP/6+np5PB7V1dUpKSkpnHIB9Ga/+pX0pz8FPvbEE1JmprP1AOhWnf39HdbMSFNTk8rKypSXl+fTnpeXp71793bqe7S2tqqhoUGDBg0K2qexsVH19fU+LwB9EDMjABRmGKmtrVVLS4tSUlJ82lNSUlRdXd2p7/H888/rzJkzuueee4L2KSwslMfj8b7S09PDKRNApCCMAFAXN7C6/D5p0xjToS2QrVu36sknn1RxcbGGDh0atN/SpUtVV1fnfVVWVnalTAC9HWEEgKQQPwk6Sk5OVmxsbIdZkJqamg6zJf6Ki4s1Z84c/frXv9aPfvSjkH3dbrfcbnc4pQGIRIQRAApzZiQ+Pl7Z2dkqLS31aS8tLVVubm7Q87Zu3arZs2frP/7jP3Trrbd2rVIAfU9cXPBjhBEgaoQ1MyJJixcv1gMPPKDx48crJydHGzZsUEVFhQoKCiSdX2I5efKkNm/eLOl8EJk5c6b+5V/+RTfeeKN3ViUxMVEej6cbLwVAxAn2qb0SYQSIImGHkRkzZuj06dNasWKFqqqqlJWVpZKSEmVkZEiSqqqqfJ45sn79ejU3N2v+/PmaP3++t33WrFl67bXXLv0KAEQulmkAqAvPGbGB54wAfdRbb0m//W3gY/feK02e7Gg5ALpXjzxnBAC6Fcs0AEQYAWATyzQARBgBYBN30wAQYQSATaFmRs6dc64OAFYRRgDYE2rPSEuLc3UAsIowAsAelmkAiDACwCaWaQCIMALAJu6mASDCCACbQoUR9owAUYMwAsAelmkAiDACwCaWaQCIMALAJsIIABFGANhEGAEgwggAmwgjAEQYAWATYQSACCMAbAr1BFbupgGiBmEEgD3x8cGPNTU5VwcAqwgjAOwhjAAQYQSATbGxUkyQH0OEESBqEEYA2BVsdqS5WTLG2VoAWEEYAWBXqE2szI4AUYEwAsAu9o0AUY8wAsCuUGGE23uBqEAYAWAXMyNA1COMALCLmREg6hFGANjFzAgQ9QgjAOzibhog6hFGANjFzAgQ9QgjAOxiZgSIeoQRAHYxMwJEPcIIALu4mwaIeoQRAHYxMwJEPcIIALsII0DUI4wAsCtUGGlsdK4OANYQRgDYlZgY/Nj33ztXBwBrCCMA7CKMAFGPMALAroSE4MfOnnWuDgDWEEYA2MXMCBD1CCMA7AoVRpgZAaICYQSAXaGWaZgZAaICYQSAXewZAaIeYQSAXf36Bf+wvO+/l4xxth4AjiOMALAv2OyIMTyFFYgChBEA9nFHDRDVCCMA7AsVRs6cca4OAFYQRgDYN3Bg8GMNDc7VAcAKwggA+5KSgh8jjAB9HmEEgH2hZkbq652rA4AVhBEA9rFMA0Q1wggA+5gZAaIaYQSAfaH2jBBGgD6PMALAPo8n+LHaWufqAGAFYQSAfcnJwY/V1vJIeKCPI4wAsM/tDr5Uc+6cVFfnbD0AHEUYAdA7DBkS/FhVlXN1AHAcYQRA75CSEvxYeblzdQBwXJfCSFFRkTIzM5WQkKDs7Gzt3r07ZP+dO3cqOztbCQkJuvLKK7Vu3bouFQugDxs1KvixP/7x/HJNa6tj5QBwTr9wTyguLtaiRYtUVFSkiRMnav369crPz9eRI0c0cuTIDv3Ly8t1yy23aN68eXr99de1Z88ePfzwwxoyZIjuuuuubrkIAH3AlVcGP1ZeLi1Y4FwtQDS78krp8ccdfUuXMeFtU58wYYLGjRuntWvXetvGjBmjO++8U4WFhR36P/7449q+fbs++eQTb1tBQYH+93//V/v27evUe9bX18vj8aiurk5JoZ5HACBytbZK//iPPHEVsK0bw0hnf3+HtUzT1NSksrIy5eXl+bTn5eVp7969Ac/Zt29fh/7Tpk3TgQMHdO7cuYDnNDY2qr6+3ucFoI+LiZFuuMF2FQAsCCuM1NbWqqWlRSl+G81SUlJUXV0d8Jzq6uqA/Zubm1Ub5GFGhYWF8ng83ld6eno4ZQKIVFOmSP3CXj0GEOG6tIHV5XL5fG2M6dB2sf6B2i9YunSp6urqvK/KysqulAkg0iQnS3/3d7arAOCwsP4ESU5OVmxsbIdZkJqamg6zHxekpqYG7N+vXz8NHjw44Dlut1tutzuc0gD0FTfffP7Ome3buXsGiBJhhZH4+HhlZ2ertLRUf9fur5fS0lLdcccdAc/JycnR22+/7dP23nvvafz48YqLi+tCyQD6NJdLys+Xxo2TPvxQ+uIL6euvzweU5mYeDQ/0tPh4x98y7MXZxYsX64EHHtD48eOVk5OjDRs2qKKiQgUFBZLOL7GcPHlSmzdvlnT+zpmXXnpJixcv1rx587Rv3z5t3LhRW7du7d4rAdC3pKRI3P4PRIWww8iMGTN0+vRprVixQlVVVcrKylJJSYkyMjIkSVVVVaqoqPD2z8zMVElJiR599FG9/PLLSktL0wsvvMAzRgAAgKQuPGfEBp4zAgBA5OmR54wAAAB0N8IIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwKqwHwdvw4WHxNbX11uuBAAAdNaF39sXe9h7RISRhoYGSVJ6errlSgAAQLgaGhrk8XiCHo+Iz6ZpbW3VqVOnNHDgQLlcrm77vvX19UpPT1dlZWWf/cybvn6NXF/k6+vX2NevT+r718j1dZ0xRg0NDUpLS1NMTPCdIRExMxITE6MRI0b02PdPSkrqk/+DtdfXr5Hri3x9/Rr7+vVJff8aub6uCTUjcgEbWAEAgFWEEQAAYFVUhxG3263ly5fL7XbbLqXH9PVr5PoiX1+/xr5+fVLfv0aur+dFxAZWAADQd0X1zAgAALCPMAIAAKwijAAAAKsIIwAAwKqoCSNPP/20cnNzddlll+nyyy8P2KeiokLTp09X//79lZycrIULF6qpqcmnz6FDh3TTTTcpMTFRw4cP14oVKy76zH0bPvjgA7lcroCv/fv3e/sFOr5u3TqLlXfeqFGjOtT+xBNP+PTpzJj2Rl9++aXmzJmjzMxMJSYm6qqrrtLy5cs71B7J4ydJRUVFyszMVEJCgrKzs7V7927bJXVJYWGhrr/+eg0cOFBDhw7VnXfeqU8//dSnz+zZszuM1Y033mip4vA9+eSTHepPTU31HjfG6Mknn1RaWpoSExP1t3/7tzp8+LDFisMT6OeJy+XS/PnzJUXm+O3atUvTp09XWlqaXC6X3nrrLZ/jnRmzxsZGPfLII0pOTlb//v11++2368SJE91ea0Q8gbU7NDU16e6771ZOTo42btzY4XhLS4tuvfVWDRkyRB9++KFOnz6tWbNmyRijF198UdL5R+ZOnTpVkydP1v79+/XZZ59p9uzZ6t+/v5YsWeL0JYWUm5urqqoqn7Zf/OIX2rFjh8aPH+/TvmnTJv34xz/2ft2Zp+X1FitWrNC8efO8Xw8YMMD7350Z097qT3/6k1pbW7V+/XpdffXV+vjjjzVv3jydOXNGq1at8ukbqeNXXFysRYsWqaioSBMnTtT69euVn5+vI0eOaOTIkbbLC8vOnTs1f/58XX/99WpubtayZcuUl5enI0eOqH///t5+P/7xj7Vp0ybv1/Hx8TbK7bJrr71WO3bs8H4dGxvr/e/nnntOq1ev1muvvabRo0dr5cqVmjp1qj799FMNHDjQRrlh2b9/v1paWrxff/zxx5o6daruvvtub1ukjd+ZM2c0duxYPfjgg7rrrrs6HO/MmC1atEhvv/223njjDQ0ePFhLlizRbbfdprKyMp/xv2QmymzatMl4PJ4O7SUlJSYmJsacPHnS27Z161bjdrtNXV2dMcaYoqIi4/F4zNmzZ719CgsLTVpammltbe3x2i9FU1OTGTp0qFmxYoVPuyTzm9/8xk5RlygjI8P86le/Cnq8M2MaSZ577jmTmZnp0xbJ43fDDTeYgoICn7a/+Iu/ME888YSlirpPTU2NkWR27tzpbZs1a5a544477BV1iZYvX27Gjh0b8Fhra6tJTU01zzzzjLft7NmzxuPxmHXr1jlUYff62c9+Zq666irvz/ZIHz//nxWdGbNvvvnGxMXFmTfeeMPb5+TJkyYmJsa888473Vpf1CzTXMy+ffuUlZWltLQ0b9u0adPU2NiosrIyb5+bbrrJ58Ew06ZN06lTp/Tll186XXJYtm/frtraWs2ePbvDsQULFig5OVnXX3+91q1bp9bWVucL7KJnn31WgwcP1nXXXaenn37aZxmjM2MaSerq6jRo0KAO7ZE4fk1NTSorK1NeXp5Pe15envbu3Wupqu5TV1cnSR3G64MPPtDQoUM1evRozZs3TzU1NTbK67KjR48qLS1NmZmZuvfee3Xs2DFJUnl5uaqrq33G0+1266abborI8WxqatLrr7+un/zkJz4fzhrp49deZ8asrKxM586d8+mTlpamrKysbh/XqFmmuZjq6mqlpKT4tF1xxRWKj49XdXW1t8+oUaN8+lw4p7q6WpmZmY7U2hUbN27UtGnTlJ6e7tP+z//8z7r55puVmJio3/3ud1qyZIlqa2v1T//0T5Yq7byf/exnGjdunK644gr9z//8j5YuXary8nK98sorkjo3ppHiiy++0Isvvqjnn3/epz1Sx6+2tlYtLS0dxiclJSXixsafMUaLFy/WD3/4Q2VlZXnb8/PzdffddysjI0Pl5eX6xS9+oSlTpqisrCwinuw5YcIEbd68WaNHj9ZXX32llStXKjc3V4cPH/aOWaDxPH78uI1yL8lbb72lb775xuePt0gfP3+dGbPq6mrFx8friiuu6NCn2/+ddus8i8OWL19uJIV87d+/3+ecYMs08+bNM3l5eR3a4+LizNatW40xxkydOtX8wz/8g8/xEydOGElm37593XdhIXTlmisrK01MTIx58803L/r9V61aZZKSknqq/IvqyvVd8OabbxpJpra21hjTuTF1Wleu7+TJk+bqq682c+bMuej3tz1+nXXy5Ekjyezdu9enfeXKleaaa66xVFX3ePjhh01GRoaprKwM2e/UqVMmLi7O/Nd//ZdDlXWvb7/91qSkpJjnn3/e7Nmzx0gyp06d8ukzd+5cM23aNEsVdl1eXp657bbbQvaJtPGT3zJNZ8Zsy5YtJj4+vsP3+tGPfmQeeuihbq0vomdGFixYoHvvvTdkH/+ZjGBSU1P13//93z5tf/7zn3Xu3DlvckxNTe2QBi9M0/mny57SlWvetGmTBg8erNtvv/2i3//GG29UfX29vvrqK8euqb1LGdMLO9s///xzDR48uFNj6rRwr+/UqVOaPHmycnJytGHDhot+f9vj11nJycmKjY0N+O+pN9d9MY888oi2b9+uXbt2acSIESH7Dhs2TBkZGTp69KhD1XWv/v3766/+6q909OhR3XnnnZLO/yU9bNgwb59IHM/jx49rx44d2rZtW8h+kT5+F+6ECjVmqampampq0p///Gef2ZGamhrl5uZ2b0HdGm0iwMU2sLZPiW+88UaHDayXX365aWxs9PZ55plnevUG1tbWVpOZmWmWLFnSqf4vvviiSUhI8NmkGynefvttI8kcP37cGNO5Me3NTpw4YX7wgx+Ye++91zQ3N3fqnEgavxtuuMH89Kc/9WkbM2ZMRG5gbW1tNfPnzzdpaWnms88+69Q5tbW1xu12m3/7t3/r4ep6xtmzZ83w4cPNU0895d0M+eyzz3qPNzY2RuQG1uXLl5vU1FRz7ty5kP0ibfwUZANrqDG7sIG1uLjY2+fUqVM9soE1asLI8ePHzcGDB81TTz1lBgwYYA4ePGgOHjxoGhoajDHGNDc3m6ysLHPzzTebP/zhD2bHjh1mxIgRZsGCBd7v8c0335iUlBRz3333mUOHDplt27aZpKQks2rVKluXdVE7duwwksyRI0c6HNu+fbvZsGGDOXTokPn888/Nv/7rv5qkpCSzcOFCC5WGZ+/evWb16tXm4MGD5tixY6a4uNikpaWZ22+/3dunM2PaW11YmpkyZYo5ceKEqaqq8r4uiOTxM+Z8MIyLizMbN240R44cMYsWLTL9+/c3X375pe3SwvbTn/7UeDwe88EHH/iM1XfffWeMMaahocEsWbLE7N2715SXl5v333/f5OTkmOHDh5v6+nrL1XfOkiVLzAcffGCOHTtmfv/735vbbrvNDBw40DtezzzzjPF4PGbbtm3m0KFD5r777jPDhg2LmOszxpiWlhYzcuRI8/jjj/u0R+r4NTQ0eH/XSfL+zLzwB1tnxqygoMCMGDHC7Nixw/zhD38wU6ZMMWPHju30H0idFTVhZNasWQHX599//31vn+PHj5tbb73VJCYmmkGDBpkFCxZ0+Avzj3/8o5k0aZJxu90mNTXVPPnkk712VsQYY+677z6Tm5sb8Nhvf/tbc91115kBAwaYyy67zGRlZZk1a9Zc9C+C3qCsrMxMmDDBeDwek5CQYK655hqzfPlyc+bMGZ9+nRnT3mjTpk1B95RcEMnjd8HLL79sMjIyTHx8vBk3bpzPrbCRJNhYbdq0yRhjzHfffWfy8vLMkCFDTFxcnBk5cqSZNWuWqaiosFt4GGbMmGGGDRtm4uLiTFpamvn7v/97c/jwYe/x1tZW76yC2+02f/M3f2MOHTpkseLwvfvuu0aS+fTTT33aI3X83n///YD/X86aNcsY07kx+/77782CBQvMoEGDTGJiorntttt65LpdxvTCx4cCAICowXNGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVv0/KcuGCLKT9ewAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import t\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "ndf = 50-2\n",
    "\n",
    "x = np.arange(-100,100,0.001)#np.linspace(t.ppf(0.01, ndf),t.ppf(0.99, ndf), 100)\n",
    "tsurv = t.sf(x, ndf)\n",
    "\n",
    "ax.plot(x, tsurv,'r-', lw=5, alpha=0.6, label='t pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a81fa45",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsurv[cutval]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88f19831",
   "metadata": {},
   "source": [
    "# Non dependent variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee55eab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_setosa_sepal_length = df[df[\"species\"] == \"setosa\"].sepal_length.to_numpy()\n",
    "X_versicolor_sepal_width = df[df[\"species\"] == \"versicolor\"].sepal_width.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bb9c594",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(X_setosa_sepal_length, X_versicolor_sepal_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f29a04a",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_setosa_sepal_length = X_setosa_sepal_length[:, np.newaxis]\n",
    "X_versicolor_sepal_width = X_versicolor_sepal_width[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a6cf950",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = LinearRegression().fit(X_setosa_sepal_length,X_versicolor_sepal_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91386855",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Coefficient = %.2f\"%clf.coef_)\n",
    "print(\"Intercept = %.2f\"%clf.intercept_)\n",
    "print(\"Score = %.2f\"%clf.score(X_setosa_sepal_length,X_versicolor_sepal_width))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72e2fa77",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ca984a",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_plot = np.linspace(4, 7, 300)\n",
    "fit = clf.intercept_[0] + clf.coef_[0]*X_plot\n",
    "plt.figure(1, figsize=(4, 3))\n",
    "plt.clf()\n",
    "plt.scatter(X_setosa_sepal_length,X_versicolor_sepal_width)\n",
    "plt.plot(X_plot, fit, label=\"Linear Regression Model\", color=\"red\", linewidth=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecbcb79c",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_versicolor_sepal_width = sm.add_constant(X_versicolor_sepal_width)\n",
    "model = sm.OLS(X_setosa_sepal_length,X_versicolor_sepal_width).fit() \n",
    "  \n",
    "#display model summary \n",
    "print(model.summary()) \n",
    "  \n",
    "# residual sum of squares \n",
    "print(model.ssr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3d39674",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d7fed9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_plot = np.linspace(4, 7, 300)\n",
    "fit = model.params[0] + model.params[1]*X_plot\n",
    "plt.figure(1, figsize=(4, 3))\n",
    "plt.clf()\n",
    "plt.scatter(X_setosa_sepal_length,X_versicolor_sepal_width)\n",
    "plt.plot(X_plot, fit, label=\"Linear Regression Model\", color=\"red\", linewidth=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68958b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import t\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "ndf = 50-2\n",
    "\n",
    "x = np.linspace(-1,1,200)\n",
    "\n",
    "ax.plot(x, t.sf(x, ndf),'r-', lw=5, alpha=0.6, label='t pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c6e4a16",
   "metadata": {},
   "outputs": [],
   "source": [
    "tfunc = t.sf(x, ndf)\n",
    "prob = 0.0\n",
    "n = -1\n",
    "for i in x:\n",
    "    n += 1\n",
    "    if i > 0.540:\n",
    "        print(n)\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "debd043f",
   "metadata": {},
   "outputs": [],
   "source": [
    "tfunc[154]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ff20ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = abs(-0.540)\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "400715a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "tfunc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d89851f",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.arange(0,1,0.1)\n",
    "csq = c\n",
    "\n",
    "tail = t.sf(csq,ndf)\n",
    "\n",
    "fig3, ax3 = plt.subplots(1,1)\n",
    "\n",
    "ax3.bar(c, tail, color='blue')\n",
    "\n",
    "ax3.grid(True,axis='y')\n",
    "ax3.set_xlabel(r'$N$')\n",
    "ax3.set_ylabel(r'$P(\\chi_1^2>N^2$)')\n",
    "ax3.set_title(r't-stat upper tail probability (p-value)')\n",
    "ax3.set(xlim=(0, 1))\n",
    "ax3.set_yscale('log')\n",
    "\n",
    "yticks = 1/np.power(10,np.arange(1,-1,-1));\n",
    "ax3.set_yticks(yticks)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69ce8159",
   "metadata": {},
   "outputs": [],
   "source": [
    "tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b12a5ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "0.34546555*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24d02898",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}