demo03.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hhhP1JoutiHm"
      },
      "source": [
        "# Demo:  Predicting Diabetes Progression using Mulitple Linear Regression"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xG0pLjGbtiHn"
      },
      "source": [
        "In this demo, you will learn how to:\n",
        "* Fit multiple linear regression models, by hand and using Python's Scikit-Learn package.  \n",
        "* Split data into training and test set.\n",
        "* Manipulate and visualize multivariable arrays.\n",
        "\n",
        "We first load the packages as usual."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "hw4RxJz6tiHo"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9VZbBnpItiHo"
      },
      "source": [
        "We will also load a few modules from [scikit-learn](https://scikit-learn.org/stable/), which is one of the most popular machine learning libraries in Python that does not focus specifically on neural nets and deep learning. Using the libary is a bit overkill at this point, but I want you to get some practice because it will be very useful later in the course, and for your projects."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "vOa4eL0BtiHo"
      },
      "outputs": [],
      "source": [
        "from sklearn import datasets, linear_model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rNT8XAsOtiHo"
      },
      "source": [
        "## Diabetes Data Example\n",
        "To illustrate the concepts, we load a well-known diabetes data set.  This dataset can be downloaded in raw form from [https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html](https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html). However, it is also included in the `sklearn.datasets` module, which prepackages a number of benchmark datasets for easy experimentation. It can be loaded as follows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "jZy_8ldOtiHp"
      },
      "outputs": [],
      "source": [
        "# Load the diabetes dataset\n",
        "diabetes = datasets.load_diabetes()\n",
        "X = diabetes.data\n",
        "y = diabetes.target"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bRB20EnEtiHp"
      },
      "source": [
        "It's always good to check what type of objects we have loaded in. You can confirm using the `type` function that our data is preloaded into `numpy` arrays, which is how we want it."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BUNB7wpatiHp",
        "outputId": "3658e9aa-4f78-46be-8e2f-5a3b0425fe5b"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "numpy.ndarray"
            ]
          },
          "metadata": {},
          "execution_count": 4
        }
      ],
      "source": [
        "type(X)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pirEMVObtiHp"
      },
      "source": [
        "The predictor variables are **age**, **sex**, **body mass index**, **average blood pressure**, and **six blood serum measurements** (think iron, cholesterol, lipid levels, etc.). Each corresponds to a column in the data matrix `X`. Note that these columns have been *mean centered* and had their *variance normalized* (i.e. scaled). This will not effect your linear regression model at all, but it's good to know."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4NVJK9nytiHp"
      },
      "source": [
        "When you look at the first 5 entries of the **age** column, we don't see age numbers we might expect:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "Nv1bb6WjtiHp",
        "outputId": "81db5244-2cf4-4a86-8a8f-56f62c5aa06f"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'[ 0.03807591 -0.00188202  0.08529891 -0.08906294  0.00538306]'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 5
        }
      ],
      "source": [
        "str(X[0:5,0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iML8iBAetiHp"
      },
      "source": [
        "And when you look at the first 5 entries of the **sex** column, we don't see binary categories like 1 for female, 2 for male. We see real numbers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "id": "9S_HUN0gtiHq",
        "outputId": "cf2ae843-94eb-4bfe-ed32-414f52db7427"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'[ 0.05068012 -0.04464164  0.05068012 -0.04464164 -0.04464164]'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 6
        }
      ],
      "source": [
        "str(X[0:5,1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S7j831MrtiHq"
      },
      "source": [
        "**NOTE**: This is one reason to be cautious about using built in datasets. Often a lot of information is lost during preprocessing which makes it harder to reason about the data than if you were working with the raw files. This is just a simple teaching demo, so we don't care for today, but keep that in mind when downloading data for your projects."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YQ3vqSWCtiHq"
      },
      "source": [
        "The target values `y` represent a \"quantitative measure of diabetes disease progression\" that we wish to predict. The exact meaning of this measure seems to have been lost in data translation, but we will work with it anyway.\n",
        "\n",
        "The number of attributes and samples are computed from the shape:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "h6cErvz2tiHq",
        "outputId": "63f2ed80-c1a3-484e-a7bd-0e1e2dca2b26"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "num samples=442 num attributes=10\n"
          ]
        }
      ],
      "source": [
        "nsamp, natt = X.shape\n",
        "print(\"num samples=\"+str(nsamp)+\" num attributes=\"+str(natt))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9zVEO1RktiHq"
      },
      "source": [
        "## Using Simple Linear Regression for Each Feature Individually\n",
        "\n",
        "As a first attempt to predict diabetes progression, we could try *one attribute at a time*.  That is, for each attribute $x_k$, we could attempt to fit a simple linear regression model:\n",
        "$$ y \\approx \\beta_{0,k} + \\beta_{1,k}x_k$$\n",
        "where $\\beta_{0,k}$ and $\\beta_{1,k}$ are the coefficients in the simple linear regression model using only the attribute $x_k$.\n",
        "\n",
        "The following code computes the $\\ell_2$ loss for each variable $k$ as well as the coefficients in the linear model, $\\beta_{0,k}$ and $\\beta_{1,k}$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "scrolled": false,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KjYkj3z7tiHq",
        "outputId": "6209f2b6-ce24-4197-8226-37afed517689"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0 loss=2.53e+06 beta0=152.13348416289594 beta1=304.1830745283059\n",
            "1 loss=2.62e+06 beta0=152.13348416289594 beta1=69.71535567841458\n",
            "2 loss=1.72e+06 beta0=152.13348416289617 beta1=949.4352603840387\n",
            "3 loss=2.11e+06 beta0=152.13348416289597 beta1=714.7382594960409\n",
            "4 loss=2.50e+06 beta0=152.13348416289594 beta1=343.25445188896646\n",
            "5 loss=2.54e+06 beta0=152.13348416289594 beta1=281.78459335245714\n",
            "6 loss=2.21e+06 beta0=152.13348416289594 beta1=-639.145279322535\n",
            "7 loss=2.14e+06 beta0=152.13348416289594 beta1=696.8830300922256\n",
            "8 loss=1.78e+06 beta0=152.13348416289585 beta1=916.1373745509136\n",
            "9 loss=2.24e+06 beta0=152.13348416289594 beta1=619.2228206843714\n"
          ]
        }
      ],
      "source": [
        "ym = np.mean(y)\n",
        "losses = np.zeros(natt)\n",
        "beta0 = np.zeros(natt)\n",
        "beta1 = np.zeros(natt)\n",
        "for k in range(natt):\n",
        "    xm = np.mean(X[:,k])\n",
        "    sxy = np.mean((X[:,k]-xm)*(y-ym))\n",
        "    sxx = np.mean((X[:,k]-xm)**2)\n",
        "    beta1[k] = sxy/sxx\n",
        "    beta0[k] = ym - beta1[k]*xm\n",
        "    errs = y - beta1[k]*X[:,k] - beta0[k]\n",
        "    losses[k] = np.sum(errs**2)\n",
        "\n",
        "    print(str(k)+\" loss=\"+\"{:.2e}\".format(losses[k])+\" beta0=\"+str(beta0[k])+\" beta1=\"+str(beta1[k]))\n",
        ""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E7pJNadktiHq"
      },
      "source": [
        "**Test yourself:** Why are all the `beta0` values the same, no matter what predictor we use?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EDvysV5stiHr"
      },
      "source": [
        "We see that the best prediction is given by BMI index, with loss 1.72e+06. If you convert this to an $R^2$ value, as described in class, we get a value of $R^2 = .344$, which on a scale from $[0,1]$ does not indicate too strong of a fit."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GpR75tY3tiHr",
        "outputId": "892f32d2-9bbd-45c8-fdc2-1673b95fbe68"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Rsqr=0.344\n"
          ]
        }
      ],
      "source": [
        "syy = np.mean((y-ym)**2)\n",
        "rsqr = 1 - np.min(losses)/(nsamp*syy)\n",
        "print(\"Rsqr=\"+\"{:.3}\".format(rsqr))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HQafzPrVtiHr"
      },
      "source": [
        "We can see this somewhat poor fit in a scatter plot as well where there is a significat variation from the regression line."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 430
        },
        "id": "cIvqtURLtiHr",
        "outputId": "c68d2fdc-9241-49c6-8419-0d6f669625c8"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACfr0lEQVR4nO2deXgT5fbHv0lJ0422lEpTtoLsZQcFKqgIZROFK9x7ERURvXhF8AfFBVBZq7K4gAriBopLQb1XLypQbFlEsCxWKpRWhFpAoS12hy5pm8zvjzAhy0wyM5lJJu35PI+PNPNm5p03ybznPe/3nKNhGIYBQRAEQRCEitD6ugMEQRAEQRCOkIFCEARBEITqIAOFIAiCIAjVQQYKQRAEQRCqgwwUgiAIgiBUBxkoBEEQBEGoDjJQCIIgCIJQHWSgEARBEAShOpr5ugNSMJvNuHTpEpo3bw6NRuPr7hAEQRAEIQCGYXDlyhW0bt0aWq1rH4lfGiiXLl1Cu3btfN0NgiAIgiAk8Mcff6Bt27Yu2/ilgdK8eXMAlhsMDw+3vl5fX4/vvvsOo0ePhk6n81X3Gh00rspBY6sMNK7KQWOrDE1lXCsrK9GuXTvrPO4KvzRQ2G2d8PBwJwMlJCQE4eHhjfoD9jY0rspBY6sMNK7KQWOrDE1tXIXIM0gkSxAEQRCE6iADhSAIgiAI1UEGCkEQBEEQqoMMFIIgCIIgVAcZKARBEARBqA5RBsrGjRvRp08fa/RMQkICdu3aZT0+fPhwaDQau/8ee+wxu3NcuHAB48ePR0hICFq1aoWnn34aDQ0N8twNQRAEQRCNAlFhxm3btsWqVavQpUsXMAyDLVu2YOLEiTh+/Dh69uwJAJg5cyZWrFhhfU9ISIj13yaTCePHj4fBYMCPP/6IgoICPPjgg9DpdHjppZdkuiWCIAiCIPwdUQbK3Xffbff3iy++iI0bN+Lw4cNWAyUkJAQGg4Hz/d999x1ycnKQnp6OmJgY9OvXD8nJyViwYAGWLVuGwMBAibdBEARBEERjQnKiNpPJhC+++AJVVVVISEiwvv7pp5/ik08+gcFgwN13343FixdbvSgZGRno3bs3YmJirO3HjBmDWbNm4dSpU+jfvz/ntYxGI4xGo/XvyspKAJbENvX19dbX2X/bvkZ4Do2rctDYKgONq3LQ2CqDXONqMjPIPF+G4qtGRIfpMTCuBQK06qlZJ+b+RBsoJ0+eREJCAmpraxEWFoavvvoK8fHxAID77rsPcXFxaN26NU6cOIEFCxbg9OnT+PLLLwEAhYWFdsYJAOvfhYWFvNdcuXIlli9f7vT6d999Z7eFxJKWlib2tggB0LgqB42tMtC4KgeNrTLIOa7FAHbnynY6WaiurhbcVrSB0q1bN2RlZaGiogL/+c9/MH36dHz//feIj4/Ho48+am3Xu3dvxMbGYuTIkcjLy0OnTp3EXsrKokWLMH/+fOvfbC7/0aNHO6W6T0tLw6hRo5pEqmBvQeOqHDS2ykDjqhw0tsrg6bim5xYh6bMsMA6vs76TtVP6IbFHjOPbvA67AyIE0QZKYGAgOnfuDAAYOHAgjh07htdffx3vvPOOU9vBgwcDAM6ePYtOnTrBYDDg6NGjdm2KiooAgFe3AgB6vR56vd7pdZ1Ox/lB8r1OeAaNq3LQ2CoDjaty0Ngqg5RxNZkZrNhxGrUm7q0cDYAVO05jdK82Pt/uEXNvHudBMZvNdvoQW7KysgAAsbGxAICEhAScPHkSly9ftrZJS0tDeHi4dZuIIAiCIAjhHM0vRUFFLe9xBkBBRS2O5pd6r1MyIMqDsmjRIowbNw7t27fHlStXkJKSgv3792P37t3Iy8tDSkoK7rzzTrRs2RInTpxAUlISbrvtNvTp0wcAMHr0aMTHx2PatGlYs2YNCgsL8fzzz2P27NmcHhKCIAiCIFxz+Qq/cSKlnVoQZaBcvnwZDz74IAoKChAREYE+ffpg9+7dGDVqFP744w+kp6dj3bp1qKqqQrt27TB58mQ8//zz1vcHBATg22+/xaxZs5CQkIDQ0FBMnz7dLm8KQRAEQXiKyczgaH4pLl+pRavmQRjUMcrn2xtK0ap5kKzt1IIoA2XTpk28x9q1a4fvv//e7Tni4uKwc+dOMZclCIIgCMGkZhdg+Tc5dtsesRFBWHp3PMb2ivVhz5RhUMcoxEYEobCi1kkkC1g0KIYIi5HmT1AtHoIgCKLRkJpdgFmf/OykySisqMWsT35GanaBj3qmHAFaDZbebdFxOvqI2L+X3h3vdx4kMlAIgiCIRoHJzGD5NzmcXgT2teXf5MBk5mrh34ztFYuNDwyAIcJ+G8cQEYSNDwzwS8+R5EyyBEEQBKEmxESzJHRq6b2OeYmxvWIxKt7QaLQ3ZKAQBEEQXkcJEWtjjWYRQ4BW02iMLzJQCIIgCK+ilIi1sUazNFVIg0IQBEF4DSVFrGw0C58fRgOLIeRv0SxNFTJQCIIgCK+gtIi1sUazNFXIQCEIglABJjODjLwSbM+6iIy8kkYZaeKNlOyNMZqlqUIaFIIgCB/TVBKLeUvE2tiiWZoqZKAQBEH4EFaT4egvYTUZjWnV700Ra2OKZmmq0BYPQRCEj2hqicVIxEqIgQwUgiAIH+ENTYaaIBErIQYyUAiCIHxEU0wsRiJWQiikQSEIgvARTTWxGIlYCSGQgUIQBOEjWE1GYUUtpw5FA4tnoTFqMkjESriDtngIgiB8BGkyCIIfMlAIgiB8CGkyCIIb2uIhCILwMaTJIAhnyEAhCIJQAaTJIAh7yEAhCIIgZMFkZsgLRMgGGSgEQRCExzSVekKE9yCRLEEQBOERbD0hx6y4bD2h1OwCH/WM8GfIQCEIgiAk09TqCRHegwwUgiD8GpOZQUZeCbZnXURGXglNhF6mqdUTIrwHaVAIgvBb+HQPS8Z382GvmhZNsZ4Q4R3Ig0IQhF/iSveQ9FmWbzrVBGmq9YQI5SEDhSAIv0OI7oFtRygLW0+IL5hYA4tXqzHWEyKUhQwUgiD8DiG6BwDIPF/mnQ41YaieEKEUZKAQBOF3CNUzFF81KtwTAqB6Qo2O2lpg1Srg2Wd92g0SyRKEiqHMnNwI1TNEh+kV7gnBQvWEGgEMA3z1FfDUU0B+PhAQADzwABAf75PukIFCECqlKWXmFGuIsbqHwopaTh0K+86BcS1kva6/4e37o3pCfswvvwDz5gH7919/zWQCkpKA1FRA4/3fBRkoBKFC2AgVx8mXzczZmNzmUgwxVvcw65OfoYG9MFbj0E7O6/oTjf3+CJn46y/g+eeB998HzGbn42lpwMmTQJ8+Xu8aaVAIQmU0pcycnqRId6V7WDuln2LX9Qca+/1JgRL6OVBXB7z2GtClC/Duu9zGya23Aj/95BPjBCAPCkGoDjGZOf3Zne7OENPAYoiNijfwekL4dA9mUwN25it3XTXT2O9PCuRNsoFhgB07gCefBH77jbtN+/bAyy8D//iHT7Z2WMhAIQiV0VQyc8pliHHpHswm5a+rJmy1JsVXjI3u/jyhKW2XuiUnx6Ip+e477uMhIcCiRRbjJTjYu33jgAwUglAZTSUzp68MscZmAHJ5B4TgL/fnCeRNukZpKbBsGfDWWxbhKxfTpgErVwJt2ni1a64gDQpBqIymkpnTV4ZYYzIA+bQmQvCH+/OUJl/IsKEBWL/eojN5801u42TwYODwYeCjj1RlnABkoBCE6mgqmTl9ZYg1FgPQlXfAFf5yf3LQ2LxlovjuO6BvX+CJJyweFEdatwY+/hj48UeLkaJCyEAhCBXSFDJz+soQY6/LN7EzCl1Xbtx5B7hoTAauEBqTt0wwv/0G3H03MGaMRXPiSFAQsHixpd0DDwBa9ZoBpEEhCJXiD5k5PU0ExhpijhoKQ1ONsBCBlFV/UxtXIQn9DI3Fm1ReDrzwAvDGG0B9PXebKVOA1auBuDivdk0qZKAQhIpRc2ZOuUI3vW2IsVsjfPiLcFLoqn/x+B6Ibq5XpYGrNEIS+vm9N8lkAjZtsiRb++sv7jYDBgDr1lnymvgRonw7GzduRJ8+fRAeHo7w8HAkJCRg165d1uO1tbWYPXs2WrZsibCwMEyePBlFRUV257hw4QLGjx+PkJAQtGrVCk8//TQaGhrkuRuCILyC3InAWENsYr82SOjUUtEJo7EIJ4VqaR4a2tEr46pWGvV26f79wMCBwL//zW2cxMRYjJejR/3OOAFEelDatm2LVatWoUuXLmAYBlu2bMHEiRNx/Phx9OzZE0lJSdixYwe++OILREREYM6cOZg0aRIOHToEADCZTBg/fjwMBgN+/PFHFBQU4MEHH4ROp8NLL72kyA0SBCEv/h662ViEk03COyAT/rBdKorffweefhr48kvu44GBlnwnzz4LhId7t28yIspAufvuu+3+fvHFF7Fx40YcPnwYbdu2xaZNm5CSkoIRI0YAAD744AP06NEDhw8fxpAhQ/Ddd98hJycH6enpiImJQb9+/ZCcnIwFCxZg2bJlCAwMlO/OCIJQBH9PdNaYhJOk4RGOmrdLBXPliiVXyauvWlLVc/G3vwGvvAJ06uTVrimBZA2KyWTCF198gaqqKiQkJCAzMxP19fVITEy0tunevTvat2+PjIwMDBkyBBkZGejduzdiYmKsbcaMGYNZs2bh1KlT6N+/v2d3QxCEFaUq2fq7B6KxCScbnXeAcMZstuQpWbQIKCzkbtO7t0Vncs1B0BgQbaCcPHkSCQkJqK2tRVhYGL766ivEx8cjKysLgYGBiIyMtGsfExODwmsDWlhYaGecsMfZY3wYjUYYjUbr35WVlQCA+vp61Nuoldl/1/MpmAlJ0Lgqh1Jjm55bhFW7fkVhpc2qOjwIC8d1R2KPGBfvdE90SDPoA9xn38grrMCh35phYFwLr0+W7sZ1yfhuSPosCwD31siS8d1gNjW4TJmvNm5qHw7A4s5Xsu/0PFAGvnHVZGRAO38+tJmZnO9jWraEeflymB9+GGjWjD+CRyWI+d6INlC6deuGrKwsVFRU4D//+Q+mT5+O77//XuxpRLFy5UosX77c6fXvvvsOISEhTq+npaUp2p+mCo2rcigxtvO7O75Shbr8TN4iemJYM0hAI+MZFOeewe5cz68nFVfjutrFPcg1To0Zeh4oAzuuQX/9hZ4ffYS2P/zA2c4cEID8O+/E6SlTUB8Wxl9fR2VUV1cLbivaQAkMDETnzp0BAAMHDsSxY8fw+uuvY8qUKairq0N5ebmdF6WoqAgGgwEAYDAYcPToUbvzsVE+bBsuFi1ahPnz51v/rqysRLt27TB69GiE2wiA6uvrkZaWhlGjRkGn04m9NYIHGlflkHtsTWYGY9YdsPOc2KIBEBMehN3zbvPIq5GeW8TpgeC7JgCsndLPY++NUISOq8nMIPN8GYqvGhEdplfE26OkN8sX0PNAGazjOnQo9K+/Du2rr0JTU8PZ1jxuHEyrV6N99+5o7+V+egq7AyIEj/OgmM1mGI1GDBw4EDqdDnv27MHkyZMBAKdPn8aFCxeQkJAAAEhISMCLL76Iy5cvo1WrVgAs1mJ4eDji4+N5r6HX66HX651e1+l0nD8QvtcJz6BxVQ65xvanvBKcLzPCOTfrdc6XGXH8zyseCQbH9WkLjTZAcJE6DYAVO05jdK82Xt3ucRxXLl3O0K7KGQmp2QV4POWXa0bc9fu+UGbE4ym/+HWYKz0PZIZh0ObAAQTPmQPNn39yt+neHVi7FtqxY/02DbyY74woA2XRokUYN24c2rdvjytXriAlJQX79+/H7t27ERERgUceeQTz589HVFQUwsPD8cQTTyAhIQFDhgwBAIwePRrx8fGYNm0a1qxZg8LCQjz//POYPXs2pwFCEGpAbrGp7fmiQ+TNlehNAautOPPQ2b+wfl8eb1s1RPbIlVhOKP4ejk14kWPHEDB3Lm7KyOA+HhkJLF8OzJoFNCGjUNTT8fLly3jwwQdRUFCAiIgI9OnTB7t378aoUaMAAGvXroVWq8XkyZNhNBoxZswYvPXWW9b3BwQE4Ntvv8WsWbOQkJCA0NBQTJ8+HStWrJD3rghCJuSe1BzPpw9gsGaQZRtgXJ+2HvfX2yG0bOim2iN72MRyjsYCm1hOCU+Gv4djE17g0iVLrpItW7g9Ilot8NhjFuMkOtrbvfM5ogyUTZs2uTweFBSEDRs2YMOGDbxt4uLisHPnTjGXJQifIPekxnc+AEj6LAsabYDHk6SvQmjVnFvEV54MtRtthA+prQVeew146SWgqoq7TWIisHYt0KuXd/umIvx1G4sgFMXdpAZYJjWTWVixe1fnYxFzPj58VSFYaNp1X+QW8VVqezUbbYSPYBjgP/8BevQAnnuO0zhhOncGtm+3ROU0YeMEIAOFIDiRe1Lz5iTpi9ojvjKMhOArT8agjlEwhPMbH7402ggfkJUF3HEH8I9/AOfOOR1mmjfHqenT0XD8ODBhAqAhXRJVMyYIDuSe1Lw9Sfoiu6ha0677ypORllOI2gbubGm+NtoIL3L5sqXS8PvvWzwojmg0wCOPoGHpUpzNzERXChixQgYKQXAg96Tmi0nSF7VH1Jh23Re6HFd6IwCIDNFh5aTefhtiTAigrg544w0gORngy/1x222W9PT9+6s+A6wvoC0eguBAbk2FmjUacsMaRhP7tUFCp5Y+9xB4e/tJiN5I30yLUfH8ySkJP4ZhgG++AXr2tFQc5jJO4uKAL74A9u+3GCcEJ2SgEAQHck9qrs7HQu5+5fCmLsed3ggACiuNsotyCRVw6hQwZoxFQ3L2rPPx0FDghReA3Fzg738nnYkbaIuHIHiQoqlwldSN73yAJQ28P7j7laqQ7A28tf1E4cVNkJISYOlS4O23ARNPlcYHHwRWrgRat/Zu3/wYMlAIwgViJjUhSd0czxcd0gzFuYf9oiaLtzOxKoE3dDkUXtyEqK8HNm4Eli0Dysq42wwZArz+OjBISIVNwhba4iEINwjRVLCiSEfPCJvULTW7gPN8/qI5EXN/TZ2mpDdq0uzeDfTtC8ydy22ctGkDfPIJ8OOPZJxIhAwUoklgMjPIyCvB9qyLyMgr8TghmuO55Uzq5iv4xqix3J9UxH531JwTprGj5O/cyunTwF13AWPHWrQkjgQFAUuWWNrdfz/pTDyAtniIRo/SWxONoeaKqzGKCA70q/tjdTKA5bMZ0rmVZGNA6ndHrTlhGjOKb0GWlwMrVgBvvgk0NHC3ufdeYPVqoH17z69HkIFCNG68USTO30WR7sZoxtAOgs6jhvtjJ6nSqzVYMwh4eMsxRIUFS5qkPP3uqCEnjD+LmsWg6O/cZLIkWXv+eaC4mLvNwIGWfCbDhkm7BsEJGShEo8VbReL8WRQpZIy2Z10SdC5f35/tJKUPuP66lElKru+OL5LlsTQGUbMQFP2d79sHzJsHnDjBfdxgsBT8mz7dUnmYkBUaUaLR4q36N94QRSq1ty5kjEqq6hAVqlO16FNunYyn3x2vaCFc0JREzYr8zn//HZg0CRgxgts4CQwEFi4EfvsNmDGDjBOFIA8K0Wjx1tYLK4qc9cnP0AB2k6QcokghK2Gprnyh935PvzbYfOicIvcnB3LrgDz57vjac+Etz6FakPV3fuWKxSPy2muWVPVcTJoEvPwycOONInpJSIEMFKLR4s2tF6VEkUL21gFInhCF3ntivAE3d4xSrehTrkmKNfTOFF0RdD7H8ROqhWCvU1hZi9KrRkSFBsIQESyLRkQuY81f9Cuy/M7NZmDLFuDZZ4HCQu42ffpYdCZ33CG+k4QkyEAhGi3eLhIntyhSyEp40ZcnUVbtXGRMqO5CzBgFaDU+F33yIcckxeX54IPruyPUc2E2M0jekct5HTk8LZ4aayYzg/V7z+KDQ/kor7n+3WL7NrJbtOS+KYHHv/ODBy06k8xM7uPR0Zb09P/6FxAQwN2GUATaOCMaLb7IRyFnoTwhK2Eu44Q9BrjWXbAr5Dt7GXgf7ID9GKmtECCLpzogPs0G37kA5++OUM/F4ynHedsVyKAR8cRYS80uwMAX0rA2/Tc74wS4bvSm5xZJ7psSSP6dX7hgCQu+9VZu46RZM2D+fODMGeDf/ybjxAeQgUI0arxZJE5uPNXGuBIHpmYXYNjqvZj63mFsOnQOAOD4/FbrGHEJUKVOUiYzg0Nni7HwvyddVh+2hW9c5AqzZuBZ4jupxlpqdgEe++RnlLsxelft+lVSv5RE1O+8qsqSSK1bN+Czz7hPOH48kJ0NvPoqEBmpXMcJl9AWD9HoUUM+CinIFbbrOHHy6SSYay88PLQDRsUbVDlG7gSorA6o9GqN9TifTkbMlg4AzLmjM4Z2juYdFznDrD1JfCdFtM1uT7mDAVBY6ft8N1y4/Z0zDJCSAixYAFy8yH2SHj2AtWstFYkJn0MGCtEk8GU+Cqm421sXiu3EKUQnsSu7EM+Nl2/rSy6xpVAB6qh4Aw6fvYzi3MPYPP1mzkyyfOdyRZeYMJffIbk+LxZPPDJiRdvutqeE4mthLe/v/OhRS82cw4e539iiBbB8OfDYY4BOp2wnCcGQgUIQKkXISjgiRIeK6nrB4kBvp+WXK+RWbOjsoI5R2JkLzgnS1blc4c5D4u7zkvt67hDjOZRje8rX4dWcXLwILFoEfPwx9/GAAGDWLEs14pb+tYBpCpAGhSBUjLu99VWTegMQprtg9RZCkGvCkitZmJzJuKR4CwzhekHRXq4+r7fuG+BSG2KLXInvhIqaxRhDhnBuca2qEsPV1AAvvgh07cpvnIwaBfzyi6W2DhknqoQ8KAShctythIW48sXqLTxdvcudLEzOZFxSjK/aBjPScgoFeQJcfV5aLTDrk59dvl8D7ye+E7o9pQGwcFx31OVfj3pRVWI4hgH+8x/g6aeB8+e523TubEnEdtddVGlY5ZCBQhB+gCsNjTsDRozeQq7cMHJvJcmZdE+K8VVRXS+qng/f58WnDWHx1ZaIq+0plhYhOqyc1Bsju0VjZ/7111VTzfv4cUs+kwMHuI+Hh1uid554wpKqnlA9ZKAQRCOAb0IUo7ew3RYCgIy8EhRW1KC0qg5RYXoYwuVPoS+0nZxJ96SIWeX0BNgalEpkkvWkX1zGU2SwDjOGdsCcEV0QoNWgvt4+DNnn1byLiiyVhjdtuh6KZotGA8ycCSQnA61aKdMHQhHIQCGIRowYvYVGA8y8tSMAYNjqvR6t8OUuMyBnvSMh3gIu5PQEeCuqTGxUjZSQfJ9V8zYagTfesBgeV3hKE9x+uyU9fb9+8l6b8ApkoBBEI0bMqtXMAO8cyAeQz9umoEL+FPpCkbPekbutFlco5gmQGaFRNVxGjBjjydslJcAwwNdfA08+CeTlcbfp0AF45RVLYT/SmfgtZKAQhIz4Og+EI7KvWq/hbqtDqQrPcibdczxX8RUjknfkun2fUmMqJ0JzxsgRGqx0NW87srOBpCQgPZ37eGgo8NxzljZB6v+cCNeQgUIQMqHGPBByJw8DhG91KFXhWc7tEdtzmcwM3tx3ljfVOwBEhujk8wTIiK1hHB2mx7KvTwkqWjg75bhbI0YISn3WVoqLgaVLgbfftlQe5mL6dOCll4DWrT27FqEayEAhCBkQumL1NlL1FkIQstXhr2UG+FBjr8WGkLMG5vPbs2UNDVbks66vB956y5JIrbycu80tt1h0JjffLP06hCohA4UgPERVeSA48ERv4QoxAld/KDNwNL/UpfcEsFSPVjxcVgRSUvazlFbx36tUQbCsn/WuXZZqwr/yFCds2xZYs8ZSkZh0Jo0SyiRLEB4iNcsp65ZnzyG1eq0QxvaKxcEFI/DpI4MRGexZrRG+arj+jlzhslzVlpVAasp+MfhEEPzrr8Cdd1r+4zJOgoMt2z2nTwNTp4o2Trz1+RCeQx4UgvAQKRMb65YvvVqDNYOAh7ccQ3BgoF2+CbkJ0GowtEs0Vk3ujcfcZDN1h7cznbqCS5gsBTnCZb2pQ5Ja4E8DICo0ECVVdW7bKiEI5hWSl5UBK1YA69cDDQ3cb546FVi9GmjXTtK11agTI/ghA4UgJGD7kC2+YhT0HvZhb+uW1wdcP15eU4+16WfwwY/nsGpSb8UemGN7xSIpsSvWpv8m+r1qe5jzTThLxncTfS5Pw2WV0iHxTehSvBusSZk8sReSd+R4LzT4GlyfV9swHd4x/oyeG18BSkq433jTTcDrr1v0Jh5cW406MYIfMlAIQiRcD1mtxpJHhAvbh70Qt3y5yLTqUpgzojO2Hj2Pwkpu44pdZT87rjvKa+pFZ5L1Bq4mnKTPsrB6kLjzeRIuq5QOydWKX4p3wzaqhq0LpHho8DW4Pq9bzmVhyZ730L2Yp26OwQCsWgVMmwZopSsS1K4TI7ghA4UgRMA3KboyToDrD/uMvBJBbnkGyj4wA7QaLJvQ01q4jmuCevGeXqpdUbqbcGzbiVHcSA2XVaIejbsV/4Zr1ZHdeUFe+XtfFFcZnaJqFA8NtsHx82pfVoDn9m3CmDOHud+g11sEsosWAc2be3x91dQLIkRBBgpBCESI98PRk+L4sBfjllf6genNCUpuhEw4AJB5vgxDu8aIOreUcFm569EIWfEn78jB4vE9MDvluEsvyNAu0bzX8VYYOPt5hRmrMSfjM8z4aTv0Jh6dyeTJwMsvAx07ynZ9n9cLIiRBBgrhVdSWaVUMQkSJZgZYPL4HopvrOe9PrFte6Qemv+YpETouxVeF6YMcERsuK3c9GqEr/haheo+NTG+EgV8ur8I/f/kOT//wEW6oKudsk9OqI0pfXI1h//qH7Nf3Wb0gwiPIQGlC+No4UFJB7417EzopRjfXY2K/NpzHWCGm0OgLbzww/SVPiS1CxyU6TK9wTywIEdhGhQaisKIGGXklsnpkJvZro24j84cfMHLWHEw8dYLzcHFIBF69dRo+6zMKn94xVJEueL1eECELZKA0EXwdXqekgt5b9ybHKowVYroL86UHpmuETDgAMDCuhVf64y5jLwOgpKoOSZ//AsD991Psd02VRub588AzzwCff44wjsN12mb4cODdeHPovbiqD1X0++7VekGEbFCitiYAaxw4rtpZ4yA1u0DR6wsRNC7/JkdSwiRv3hs7KfI9woQmMBvbKxZvPzAAkSHc8k16YLqHnXAA5/TzGod23oLV9Bgi3BsX7r6fcn3XfEJVFbB4MdC9O/D555xN0joPwuhHNuClEY/gqj4UgPLfd77PxxARRCHGKkWUgbJy5UrcfPPNaN68OVq1aoW//e1vOH36tF2b4cOHQ6PR2P332GOP2bW5cOECxo8fj5CQELRq1QpPP/00GvgS8xAeoaRxIBSpmVbd4e17EzIpCn3Iju0Vi8znR2H28M5Ox2wfmJT1kh9XE87aKf181qeDC0Zg68whWDulH6JCuY1Qd99POb9rXsNsBj75BOjaFXjhBaDW+Td/pVNXzJ2xCjMnL8G5KMs2qDcNBNvP5/V7+2HrzCE4uGAEGScqRdQWz/fff4/Zs2fj5ptvRkNDA5599lmMHj0aOTk5CA0NtbabOXMmVqxYYf07JCTE+m+TyYTx48fDYDDgxx9/REFBAR588EHodDq89NJLMtwSYYsawuuUUtD74t74Il9iwvWYOqg9jA1mQRoDwDIJzRreCTt3nsbm6TejuLrBTj+g9NaVrzRJcl6XT+RrNjVgZ77MHRcIu92SkVfiUb0bv4qyOnIEmDvX8n8uoqKA5cvR/LHH8Jo2APf6UC+jyu0wghNRBkpqaqrd3x9++CFatWqFzMxM3HbbbdbXQ0JCYDAYOM/x3XffIScnB+np6YiJiUG/fv2QnJyMBQsWYNmyZQgMDJRwGwQfagivU0pB76t7c5wUzxVXY+vRC1ibfsbaRqwhMahjFHS666ttpbNe+kqTpMR1uSYcs8mjbsqCHN9P1UdZXbwILFxo8ZxwERAAPP64pRpxlGU7KgAgA4EQhEci2YqKCgBAVJT9Puinn36KTz75BAaDAXfffTcWL15s9aJkZGSgd+/eiIm5nptgzJgxmDVrFk6dOoX+/fs7XcdoNMJovB4uWFlZCQCor69Hff31FQr7b9vXmjrRIc2gD3C/LRAd0ox33Dwd1/5tmyOuhR5FlfyCxpjwIPRv21zUNeS4N0+4qX040nNrsHHfaae09WVXazBvaybWTumHxB78eTi4xtZkZrByxykE8tybBsDKHacwvEtLSRNVem4Rkj7L4uzz3K2ZmDY4Dnd0b4WBcS1knQjdXXf28E6IaxmK6DC9x9dWw7NAzu/nTe3DAYQDAMymBp8aYPX19dAajWCSk8G8+io01dWc7cyjRsH08stAfDz7Ri/20v9Qw3fWG4i5Pw3DMJI2tc1mMyZMmIDy8nIcPHjQ+vq7776LuLg4tG7dGidOnMCCBQswaNAgfPnllwCARx99FOfPn8fu3but76murkZoaCh27tyJcePGOV1r2bJlWL58udPrKSkpdttHBEEQhIIwDFofOoSeW7Yg5K+/OJtcbd0a2TNmoOimm0RXGiYaP9XV1bjvvvtQUVGB8PBwl20le1Bmz56N7OxsO+MEsBggLL1790ZsbCxGjhyJvLw8dOrUSdK1Fi1ahPnz51v/rqysRLt27TB69Gi7G6yvr0daWhpGjRpl5y5v6rArV4A7vE7IKl+OcU3PLcKqXb+isNJmPz08CAvHdXd5fXfn9OTePOFofike3nLMbbvN02/mjbZgx1YX1w+rd5+xGxt3rJncB3f2FrctIrTPgPgxNJkZZJ4vQ/FVo5MXRMx1pVzbEU+/s3J9V335/ZSd48cRMH8+tIcOcR5mIiJgfu456B9/HANpq140TWX+YndAhCDJQJkzZw6+/fZbHDhwAG3btnXZdvDgwQCAs2fPolOnTjAYDDh69Khdm6KiIgDg1a3o9Xro9c4Jl3Q6HecHyfd6U2Vcn7bQaAM83vv3dFzH9WmL0b3ayLqfLte9SaG4ugFGk/u+F1c3uB23+V+cRK1JA+eYDX5aRYSK/jyE9plFA2DFjtMY3auNy8/JnbZE7HXFXNsVUr6zqdkFeDzll2sGxfXrni8z4vGUX0Tpf3z5/ZQL06UClMx9Gjf8NwUaLoe7VgvMnAlNcjICbrgBAc4tCBE09vlLzL2JMlAYhsETTzyBr776Cvv370dHAbUSsrKyAACxsZYfYkJCAl588UVcvnwZrVq1AgCkpaUhPDwc8exeJSE7ahHbKaGg99W9ySH+ZUNMxeyzik3iZhs1U3xFXOp3IZFQQgS9UjLi+qKAm7t6SwyAZV+fElXEUS2/PdEYjTi9cAXablyHVkZunQmGDwfWrQP69vVmzwiF8XXWcRZRBsrs2bORkpKC7du3o3nz5igsLAQAREREIDg4GHl5eUhJScGdd96Jli1b4sSJE0hKSsJtt92GPn36AABGjx6N+Ph4TJs2DWvWrEFhYSGef/55zJ49m9NLQshHYw6vE3pvcv7w5EifnXm+TNQ1xebA4PJsOBY0FAJfpInQMvbfP32Hy7GScm0lEFJvqbDSiPV7z2JuYhfB5/Wr3x7DANu3o/qJeej253nOJlUxMTj9TDIGJP2LdCaNDF9nHbdFVKK2jRs3oqKiAsOHD0dsbKz1v88++wwAEBgYiPT0dIwePRrdu3fHk08+icmTJ+Obb76xniMgIADffvstAgICkJCQgAceeAAPPvigXd4UglCC1OwCDFu9F1PfO4y527Iw9b3DGLZ6r+Rss3Ik0xJbzE5MUiu+LLtScr3xeUCE5qLJPF/GO1ZSr60EhRU1gtqtTf9N8QzMPuHkSSAxEbjnHoRwGCdXA4PxyvAHsffNN5FU3xEmyhvYqPB11nFHRG/xuKJdu3b4/vvv3Z4nLi4OO3fuFHNpgvAIpfKKeJpMKzpMj2IB15lzRycM7XyDYI+Pu60KQJgnxZ0XSGxRO66xknptJSitqhPcdvk3OaK2elTNX38BS5YA775ryQjLwRe9ErHm9gdRGdECawJNKKz07vYboSxCvaHe/M5TsUCi0aP0D88TjcHAuBbYncvvVWAn6aRR3UT1TchWhZkBFo/vgUvlNdh06BzntQHXXiCxOhyuJHfr0n8DoI4CblEiqh97Wx+jCPX1wIYNlkRq1/JaOXKsTTxWjJyJk7GWLS29zSflze03QlnUkHXcETJQiEaPN354UjUGtpOvnFVWhU4c0c31eOTWG3FzxyhJXiApOhzHsepmCFNNOndDuDLZjFXJzp3A/PmAQz01lovNb8Cq4Q/hmx638epMvLn9RiiLGrKOO0IGCtHoUeMPz5G1U/phxY7Tsk3SYjwbJjODiOBAPDO2O0qvGhEVGghDRLDgekKelrFXU5QLa3AJ2YIC/HSCzs21GCYOpUusBAfDvGAB7sdNOO9CkmMIV2b7TS0RJE0NpUqSeAIZKESjR40/PEcSe8TImiNGqGejrMqIYav3cir2xYTRCtHhuJp41BLlYmtwuZLn+EIf4zFlZcDy5ZYtHb7q8ffdB6xaBW27dlh4TbfFZ3guHNdddsNBTREkTQ05ohLlhgwUotGjxh8eF1ImacdJf2BcC2SeL8PlK7W49+b2WJf+G+8EM6FvLGanHJdFOOzOC+JPEw9rcC388iTKq53rhvhKHyOZhgaL+HXJEqCkhLvNzTcDr78OJCRYX+Kv3B0EoEr2DLhKF8gkXCOHN1RuyEAhGj1q/OHJgZAcJ5EhlqyNthOtISIIi8fHI3mHvMJhPgPLHyce1uBav/cMPjh0DuU19uOnRsOKk/R0ICkJyM7mPh4bC6xaBTzwgCUjrANchmf/ts2xO3WXrN1UYwRJU8TTqES5IQOFaBKwP7xlX59CYeX13CMx4Xosm9DTPyYbG/gmfcew4YrqejAAkhK7oEN0qNWz4S3Fvj9PPAFaDeYmdsWcEV28qomQRYNx9izw5JPA119zH9frLccXLQLCwlyeytHwVKLarhojSLhoCvoYNWnCyEAhmhh8KdX8ByE5TlhYI2DbsT9wcMEI60NGqCB417XETFIfUEpMPN6eJLypj/F4K6yyEnjhBUv6eT5D4u9/B9asAQSUKvEW/iBk96dtSk9RiyZMVCZZgvBXWI+DY7XgokrfZEj0BCE5TmyxNQJYhAqCP8o471HGXbknHrmzAasJj7J4mkzA++8DXboAL7/MbZz06wfs3w988YWqjBNA/UJ2tWVYbSqQgUI0etxtMwCWbQaTlBzwLq6ZkVeC7VkXkZFXIuu5pa4ibd/HCoeF+h2kPojlnHikThJKfhZy4dF39MABi8h15kzg8mXn461aAe+9B/z0E3D77XJ2WzbcfR81sHgrfCFk98Xzg7BAWzxEo8fb+9tSXMFH80tRXN0gaMtC6irS9n2uhMNcSNWLyBVBJVXLkp5b5JRfRo1ueUnf0XPngKefBv7zH+436XTA3LnA888DERGy91lO1Cxk9xd9TGOEPChEo8eb+9tiV/npuUUAgIe3HBO8ZSHW+8G3+mSFw4YIYQYP11aRO+QoqAiImyRsSfosyy/c8qK+o1evWoyO7t35jZMJE4BTpyzbPSo3Tlj4vo9iCmQqgT/oYxorZKAQjR5v7W+LdQWnZhcg6bMsp7buJlBXk74j7oyAsb1icXDBCGydOQQPJsS5OZsFsQ9iOSYesZMEO8b+4pYX8t3TMGbEp28HunUDXnwRMHJUwu7ZE/juO2D7dosexc+w/T6+fm8/bJ05BAcXjPCpt0vt+pjGDG3xEI0ebyVqE7PKH3St9o3U8Fu+fAWOeVCE5C+wVex/lHHezV1KexB7GroodpLIPF/msp3a3PLuvqMDLuYief8mdPnzV+4TREUBK1YA//430My/H+tqiSBh8ZdEj40R//4mE4QAvLW/LWaVzxoz+gDuNkImUK5J3zaTrKMR4C48V8yDWEqorycTj9hJovgqh3eBA7W45fm+o4bKYiz8/kP8LWc/zxsDgNmzgaVLLUYKITtq1sc0dshAIZoE3siQKGaVL9e+Ntekz5fN1Z1wV+iDOC2n0Ov5IMROEtFhehQLOK+a3PK239Gy4nI8evQrPHbkPwip5zG2xowB1q4FevTwbkebIGrLsNpUIAOFaDIonSFRzCpfqNBUjglUTKp5dw9iAD5LWy9mkhgY1wK7c/k1Omp1y4/tacDoE/tR/+bT0Bdc5G7UtSvw2mvAnXcCGlq1ews1ZVhtKpCBQjQplNzfFrPKZ42Zsqvc9ezlmkClhOfyPYgBYNjqvS6Fp89+dRI19WYYwpV5eAudJGz/9hu3fGYmMHcutIcOQc91PCLCspUzezYQGOjt3hFQnz6msUNRPAQhI0IjVmwjcRyRcwKVGp7LPogn9muDhE4tEaDVCMpgW1pVj6TPlM3wytU3W1h9DAA8PrwTYsLtp3tfh606UVgIPPywJdnaoUPOx7Va4LHHgDNnLIX/yDghmgjkQSEImRG6yh/bKxZrp/RDXX6m3ety7WubzAwOnf1LUFshmhixglIx2z5y1ddhtTalV2uwZhCwYX8eWoQGIymxKzpEh3jNLe94P5zi5TqjpWbOiy9acptwMWKERWfSp4+i/SUINUIGCtEkUargnON57+rT2uV5E3vEYGc+sHn6zYIzyQqBSxTrCiFaF7F6GKHZZ+UqwmartbGNjiqqrMW69N+w8YEBXnHPc92PXfg3w+Dei5lYsn8zQi5e4D7JjTcCr74KTJxIOhOiyUIGCtHkUKoqqSfnHdQxCjqdTvK1HfvBJWTlQozWxZ0ImAt34dJiBLyukJoKX2747oc1TrpfzseSPe/hlgsnuE8QFmbJEjtvHqDnVKIQRJOBNChEk0KpqqRqqXbqaqJ2RKzWRUwGW0e4tofkLMImVWsjJ67uJ6q6Ai/uXo8dH87lNk40GosO5cwZYMECMk4IAmSgEE0IpaqSCjnvsq9P4dDZYsUr6goRsrJIEYuKrd/DwrU9JKdR4c16KXzVkbnuR2eqxyNHv8L+dx/F/VmpCGDMTuervGkIcOwYsGkTYDB43D+CaCzQFg/RZBBblVSoTkXIeQsrjbj//SPW19itn5Hdoj25JSeETsBz7uiMpFFdJW132IqACytqkLwjF2VVdaLTgMtpVHirXoqrbTxjg43xwTAYkXcMz+3bhE6l3PlM/gy/AauGz8CoZU9gYv+2HvWLIBojZKAQTQYxE6IYPYmUVTm79fPWfX1Fv9cVQifgoZ2jPdJi2OaDCA4MkJQGXE6jQol6KY4GalmVEbNTjvPqZeYlWorzdS6+gMV738ft+T9znrdap8fGwX/Hu4MmwajT4/7wYMF9IoimBBkofohSEShqQan7Ezohniuuwrr0M4KFm1JW5axwc9WuXzG/u+i38+KLwmZS04DL2VfHJHmO5wHE5ZXhi8RxJcLdsT8baw58jEmHv0Yzjq0cAPiy5x1Yc9t0FIZHQwOL0au2bLYEoRbIQPEzlIpAUQtK3h87IbrajomNCMLWoxdERYNIiW5hz1dYKW+xOl8VNpOSBlzuvtoaSqU2GXrF5pVxF4njdB9mE+4/vhNJB1PQovYKZ5us2K5YPvJRHG9jsUZVm82WIFQEiWT9CLVEiiiF0vcXoNVgQl/Xk9TAuBYorOSvhMsl3PQkukUJ+ISsMeF6zEvsAmODWRGhrrsMr2L6KjXb69hesTi4YAQ2T78ZgCW/zMEFIwSfR0wUFAAMyz+OXZufwIr0dziNk8KwKCSNn4/JD75iNU4AFWazJQgVQh4UP0EteR6Uwhv3ZzIz+PoX10bOwTNCauA66074tjl8haNH41xxNbYevYC16WesbdTieZO7CBtb62hnLkSfR2gUVIfSi3hu3yaMOnuU8zgTFISLj8xG1n2P4p8xLbGaK5OsH/5OCcKbkIHiJ4iNQPE3vHF/Qiaf8pp6Qefi0p04TrTRoXo8+cUvKKp0obEIDwJQJeiaXLjS67AejdTsAqxL/80nFYiFopYibO4Ez82NVXji0DY8lPkNAs0N3I3+8Q9o1qxB2w4dYBubo4b7Iwh/ggwUP8GbeR58gZT7EyumFXqNyGAdKmrqJQk3HSfaZRNcaywWjuvuVItHKEL0Oo3d8yY3fIJnrdmEKSfS8OQPHyO6uoKzTUWPXoh4ewNw221KdpEgmgykQfETvJXnwVeIvb/U7AIMW70XU987jLnbhFXPFXqNGUM7AnDWk3gi3OTTWCT2iBF0HkeE6nV8mWGVL6GZmmEFz7af7uALJ/HtlnlYuXs9p3FSGhqJk8tfQcTJLDJOCEJGyIPiJ/gifNSbiLk/qfVbhF5jzojO6GYIEx02y4crjUV9vbAtJVvEeEV85Xnz12gz28iiduWFWLRvM+787UfOtuZmOhQ8/BgMq5MRFRnh5Z4SROOHDBQ/wVfho3LDty0j9P4ASN6yEDOGXEbFwGtCx+1ZF0ULHeXUWIjxivjC8yZXAUBfMTYuDOklu9H2w7ehb+AxICdOhPaVV9Cmc2fvdo4gmhBkoPgRUhNiqQV3q2oh95eRV+KRmFbMGNoaFanZBbj95X2q8AiI8Yrc1ae1Vz1vfq15MZuBjz8GFi1CpwKercJevYB164CRI73aNYJoipCB4mfIHZLpLYSuqt3dnxxbFmLHUG0eATFeEW973vw22uzHH4F58yxF+7ho2RJITgZmzgSa0WOTILwB/dL8ELWEZApF7Kra1f3JtWUhdAzV6BEQq0fypufN76LN/vgDWLAA2LqV+3izZsDs2cDSpUCLFt7tG0E0cchAIRRHzlW1t8XCnvZdibpCUrwi3vC8mcwMiq/wZ+G1xRPNiyxjWl0NvPwysHo1UFPD3WbcOOC114DuMhZLIghCMGSgEIoj56ra21sWaTmFgtpx9V3JSBYxXhHHCf2uPq1l9/Zw3SsXnhqQHo8pwwDbtlm8Jn/8wd2mWzdg7VqLgUIQhM8gA4UQjdgVrNDVcvEVI0xmxu3k6WpyXjw+HhHBgZIibRxJzS7A5kPnBLV1vEc5dCvuxlmIV8Qb4b589+qIpwakx2N67JhFZ/Ijd9gwIiOBZcuAxx8HdDrR/fM3GntVdML/EWWgrFy5El9++SV+/fVXBAcH45ZbbsHq1avRrVs3a5va2lo8+eST2LZtG4xGI8aMGYO33noLMTHXE1JduHABs2bNwr59+xAWFobp06dj5cqVaEbiM9UjZcITWu03eUcu3j+YL2jydEorH6bH0fwSPPvVSbt09VInY1Z7IoRYB4+AGN0KH0LH2ZWWxhviXjHF9TzRvHikBSooABYtArZs4T65Vgv8+9/AihVAdLTovvkj/pqnhmhaiMok+/3332P27Nk4fPgw0tLSUF9fj9GjR6Oq6notkaSkJHzzzTf44osv8P333+PSpUuYNGmS9bjJZML48eNRV1eHH3/8EVu2bMGHH36IJUuWyHdXhCKk5xZJqjYsptqvmMrF7OSsb6bF7JSf8fqes061dKRWQhZaNA4AJvSNtZsUPc3eKkdVZ3cTOmCZ0D3N7ip0nKYNaY9nxnZHRHCgoGs6ZqE9/Lvw8HIWbV0dtKtXA1278hsnI0cCWVnAW281KeOkMVdFJxoPolwWqampdn9/+OGHaNWqFTIzM3HbbbehoqICmzZtQkpKCkaMGAEA+OCDD9CjRw8cPnwYQ4YMwXfffYecnBykp6cjJiYG/fr1Q3JyMhYsWIBly5YhMDBQvrsjZGXVrl8lR7MIrfYrNjLG3faC1EgbMVEmX2T+iSdHd0dgM62o91rahdu9JlfUkBRxrxSXv9B7/fjwBXx8+AIA9yt1rtV9ZLCwLZfLV2oBhoHmyy8xYu5cBBQVcTfs1Al49VVgwgRA03S2NdQYlUYQfHi0p1JRYalLERVlcW9nZmaivr4eiYmJ1jbdu3dH+/btkZGRgSFDhiAjIwO9e/e22/IZM2YMZs2ahVOnTqF///5O1zEajTAar0cHVFZWAgDq6+vtUoWz/5aSPpzghx3Psqs10AfwP7RKr9bg8NnLvALIkd2iMbzLrUg5ch6rd592eU135wIsD9uVO04hMMD9ilzI+WyJDmkGvYDzAkBVbR1uW52OpXfHI7FHjOD3Roc0s47t4bOXUVpjQslVI0qv1kAf4Nm9XK6oEtSHyxVVqK8PR3puEVbt+hWFlTaanvAgLBzX3a5ekMnMIPN8GYqvGhEdpkdUcIDgcWIpu1qDeVszsXZKP6daROm5RUj6LAsMYDcGNXV1LseEJe78rzDPvQ/NDhzgfLgxzZvDvGgRzE88Aej1QANPReJGytH8Uo+/X/ScVYamMq5i7k/DMIwkH6/ZbMaECRNQXl6OgwcPAgBSUlIwY8YMO2MCAAYNGoQ77rgDq1evxqOPPorz589j9+7d1uPV1dUIDQ3Fzp07MY5DOb9s2TIsX77c6fWUlBSEhIRI6T5BEI2IwPJy9EhJQVxaGjQcjzRGo8GFkSORe//9MFI+E4LwGdXV1bjvvvtQUVGB8PBwl20le1Bmz56N7Oxsq3GiJIsWLcL8+fOtf1dWVqJdu3YYPXq03Q3W19cjLS0No0aNgq4JqPC9BTuui3/Swmh27fbdPP1mt16Ko/mleHgLT8ZOEefaebIAz/z3hNvziOmbLem5RZj3WZbg9hoAMeFBeOFvvfCvj35y23728E7YdOAsVtxkFjS2tri7F5OZwZh1B1BUyZ8vJiY8CDv/71bc+cYPdp4TrnbPjOmOJ7/IcjqXbai3Y9i32PsQ+r2wRWeqx7SfvsXjB7ehubGas4152DCYXn0Vrfv3R2uR/WtsyPHbo+esMjSVcWV3QIQgyUCZM2cOvv32Wxw4cABt27a1vm4wGFBXV4fy8nJERkZaXy8qKoLBYLC2OXr0qN35iq7tE7NtHNHr9dDr9U6v63Q6zg+S73XCM1qEBeNCmdFlgrQhnVu53bse0rkVosKC3SZbc3euVhGhMJqETeqxAvtmy7g+bcFAizlbj0OolvR8mRFHzlUI6tfmH/9A7TWjxGjWCHqP0LHRAVg0vidmffIzAO58MYvG98SJS1dxvswIV/Ll82VGLP4mF7U8/dMAiAjRIahZAK+hw0dxdYP1t1pc3SD48wTDYGTeUTy3dxNuLLvE2aT6hhsQuG4dmk2dCm0T0pm4Qq7fHkDPWaVo7OMq5t5ERfEwDIM5c+bgq6++wt69e9GxY0e74wMHDoROp8OePXusr50+fRoXLlxAQkICACAhIQEnT57E5cuXrW3S0tIQHh6O+Ph4Md0hvMzCcZaMmo6PLbH5LVxF9Yg5Fxu+LGTqkZp7484+rbF+6gCR7xJmzThGHLlD7DizwmRDhH2OFkNEkDXEWKjItbSqjvcYA6C8uh6v/qMvts4cgtfv7YfF43sIOq9t/hih+XK6/HUeH32+BJv+m8xtnISEwLRsGfasXw/mH/9oUiJYd8j12yMIbyDKgzJ79mykpKRg+/btaN68OQoLLVk2IyIiEBwcjIiICDzyyCOYP38+oqKiEB4ejieeeAIJCQkYMmQIAGD06NGIj4/HtGnTsGbNGhQWFuL555/H7NmzOb0khHpI7BEjW00XOerDuMoqyxIZosOqSb09yu1wZ59YvK0dgGe/OonSKvdGRcKN0fjvzxddrlIjQnQorxZnoEgdZ1fJ3DxJOe9IcZURE/u1AWDZYnr/YL6okgTu8uVE1lRi3sEUPHB8J5oxZu5OTJsGrFwJc6tWMO/c6XSYkpP5f1V0oukgykDZuHEjAGD48OF2r3/wwQd46KGHAABr166FVqvF5MmT7RK1sQQEBODbb7/FrFmzkJCQgNDQUEyfPh0rVqzw7E4IryBnTRc5zsX3sI0M0WHGLR0xZ0RnWSagsb1iMaJ7DIas3MPrTbC6xzu1dJuOf8YtHbE2/Te31502pD00Gg3iokIwLaGDNZRZCEJS3AupbdQiVCfIMLM1dqSUJOB7TzNTA+7P2oWkg58isvYq98UHDwZef93yfwDgiBSg5GTX8deq6ETTQpSBIiTgJygoCBs2bMCGDRt428TFxWEnx+qG8A/krKbseC42QZeYh6a3HrYBWg2mJ3TgNCwcJ113q9RR8QZsO3YBZVd5CtUB0GpgzR0CQHCWXUBcJlp3hsQLE3sheUeu6AKNUlbqju+57fdMLN77PrqUcNfNKQhriUsLl2DgojmWjLAuxkPprLr+hr9VRSeaHpRbnnAJuwoHLBEAYoWmYvBkhav0w9ZdMTyuSded4bT07njM25rJe01HYa7QyVTsZCzEkNBqNZIKNEoxHsf2isWoZpW4+MjjaP/jXs42tc0C8c6gSXh78N+x+d7hLo0TSk5GEP4JGSgEL+ykXHq1BmsGAQ9vOYaosGBFXOJqXuG6y1ablNgFc0Z04ZzcXBlOY3vFYu2UfqjLtzdStBpn4wQQNplKnYzdGRKe6BZEGY/l5UByMgLeeAPteZKofdP9VqwaPgOXIloJqowsJasuQRC+hwwUghPbSdk266QSBoOaV7juiuFpAGw79gfmjOgi6fyJPWKwM9+Sd6K4ugHFV4xI3pHL297dZOrJZOzOkFB0K81kAt5/H3j+eaC4mLPJyZhOWDFyJo616yUq4kRc6QGCINQCGSiEE942GNS8wvVW3wZ1jIJOp8P2rIuC2vNNpkpPxopspe3bB8ybB5zgTrpnbHkDXrn9Qbzf+XYwGstWjpiIE6GRSo7tKOKHIHwLGSiEE3JMymIe7mpe4crZN64xcUTqZCrX+73K778DTz8NfPkl9/HAQCApCfpnn8XCsOYYIdFYGNQxCpFuwrpbhOjsPg+K+CEI30MGCuGEp5Oy2If7ueIqQdfzxaQq14TPNyZLxnezayck7NeV7sLT93sDU0UlChcsgWHzRgTU8ySAu+ce4OWXLVWHAQQAinrPbMdKzXoogmhKiMokS/g/bBjv9qyLyMgrgYlDjenJpMw+3B09MOzDPTW7wKn92vQzLq+jgWUy98Wk6i5brZC+uRqTJIdaP55m+lR1plCzGSeT16KsXUe0eed1buOkd29gzx6LV+WaceIpR/NL3SbFK6+ux9H8Urfbm4Ble5Prd0MQhLyQgdKESM0uwLDVezH1vcOYuy0LU987jGGr9zoZDVInZbEPd7a9EHw1qQqd8AFwGn5CxoRtx8KXoj4qNBAzhnZARHCgywlSSIp7sQgxbF1y6BAqevdH7yXzEX2l1OlwaXA4Ti1eBfz8MzBihOj+uUKMR1DM9iZBEMpCWzxNBKFua1YnMa6XAZsPnRO1CherXXHXnmVeYlefutTZCX/Z16dQWGm0vh4TrseyCT0BAMNW7+Xc0ooIDnQ7JgCQeb4MQ7vG2F2TjZhJzynEV1kXUVJVh82HzmHzoXNu9RByRtx4pMe4cAFYsADYtg0RHIfrtQH4aMBdeGPoVIQER+OgNgABHO08QYxHUM16KIJoapCB0gQQGpVjNjNI3pFrNxE51llzFT0h9uEutH2H6BBB7ZTH2Vw7fqEM7x7I5zX8ZgztIOjMxVeNTq8FaDWoqLEYJVL0EHJE3EjWY1RVAWvWWHQkNdzZcvfeeBNeHPEI8lq2AwBUVNTicF4JtFqNrJEzYnQ5Qj0jvtBDUVQR0dQgA6UJINSz8XjKcadjtp78zdNvdplJVqx2xV/CP3kn6cpavHMgn/M9rOG3PYuj2i4H0WHOhTJ9nR9G0vUZBkhJsXhNLnKHTJ+NaovkkTPx/Y0DnY7NTvnZrsqzHJEzYuoCqVVkTFFFRFOEDJQmgKfuaPYhPjCuhd1E6GgoDIxrIerhLmUy8PaD2l2iNlcwAEqq6hAVqkNZVT3vPQKWsXXEWzlYbD/H6DA9wFgqExdfMYq7/rFjwNy5QEYGZ/sKfSjWDrsfn/S/Ew0B3I8eW+MEkC9yRmgmXFtjxhFfiYwpqohoqpCB0gTw1B3NpZPgMxQm9I116VWwfbiLrXjriwe1UJ2MK+7p18aq5+G6RwCcE5439BDuagwJoTLvHLBiPvDRR9wNtFqY//1vTIlKxOkGvShjT05PkRhdTgRH3pTIEB1WTurtVWPA1140gvAlFMXTBHAXlSMUVifhKmyWzzjhQ2jEia/CP+UQQybGG3jvce2UfrzvUzrpGt/nKBR9vRGzf/wMiffczm+cJCYCv/wC7VtvYd69CQCclTzukDNyhtXlTOzXBgmdWjpN6uyYcIUll7kJVVYCiioimjLkQWkCuPNUCJ3So8P0gsNmuZBaqA7wXTp8T7xPtltUAVoN5z2aTQ3YyWPTKamH8GTrCgyDcacP4bl9m9G28jJ3m86dgVdfBe6+26q05ttmiQzWOW3tcMEai0ppkITUXfK2t4KiioimDBkoTQRXe/CLx8cjeUeOy4kQsOgkPNny8KRQna8e1KyRIPWeHbe0HO/RbOJ/r9gtMDFI/Rzji37H0j3vYvAf2dwNwsOBxYuBJ54A9M7CXy5j1MwwuP/9I26v3ap5kKIaJDXWhPKr0gUEITNkoDQhXHkqtFq4nAgBy4QphwEg5Ry+elCzRsJjHKJJd8iRv4XPsIwJ12PqoPYwNpiRkVci2osg9jOIrirDkwc+xpQTadBymbEaDfDII8ALLwAxMc7HbXA01ExmRpCnqKyqDrNTnDVIBRW1eOyTn/HWfQNwZx/p461Gb4Vao4oIwhuQgdLE4PNUuPKwLBnfDXX5mQDkMQCknMOXD+qxvWLx1n0DMGfrzxAjcZErf4ujYXmuuBpbj16wKxEgxotgMjMovuKcd4ULnakeD/30Df7vx61oXsedzwS33QasWwf07y/onI4I8RQtHt8DyTtcb0nN2foz1qM/7uzT2u01ubaJ1OitUNKLRhBqhwwUwgqfh8VWJ+HOUHCFVCOCnUzu7GXApkPnOM8LKPugvrNPLNajP2euGD7knMhYwzI1uwDr0n+THMkkOGqHYZB49iie2/c+OpYVcLfp0MGSiG3yZECj8Ugb4i4M2F1GXsCSs+fxlON4W6vByG7RvO34tokWj49XpbdCaIg0QTQ2yEAh7HCnkxAquJVrtcc1mWg19gnkpD6oxU6od/Zpjbe1GqeU944oNZF5GnLKF6btSNe/zmHxnvdx6/ks7gahocCzzwLz5wNBQdZze6oNcbUFuT2LO+kbF8u/ycHwLrdyHnMVqj475Wc8eltHvHsgX3XeCjlLFxCEv0AGih/jq9TX7lZ0AGRZ7fFNJsy1Fx4e2gGj4g1WQyAjr0TwWEidUNmJYv3es1ib/pvTcSUnMk9EnEKidiJrKjH/h09xf9YuBDBm7kYPPgisXAm0bm09L99YSMlPw7cFKcYbVVBRi8zzZU6vCzHwvv6lABvu6+9U8kEN3go5ShcQhD9BBoqf4uvU1+5WdJ6u9oRMJruyC/Hc+Hik5RSKGgtPE74FaDWYm9gF3QxhXnW7eyLidGXcNDM14IHjOzHvUAoia69ynzQhwaIzGTTI+lJqdoFLb5KcicTERlNx1TYSauC1CNXj4IIR5K0gCB9DBoofopbU165WdJ6s9kxmBh8eyhc0mazfewbr0s8IHgs5MnOynitjgxmv/KOvNTW80hOZJyJOPuPm9t8z8fze99Gl5A/uk7VtC6xeDUydalc5Uuh2kVyhuWKjqaLD9Ch2eE2MgUfeCoLwPWSg+BmNPfW12NTrH3BU+gX4x8LTXBeuPFdKT2ieRDI5Gi03lvyJ5/e+jxG//8R9saAg4JlnLP+FhtodkpLkjc84ELNNKSSaih2DgXEtsDvX/pgao3QIguCHDBQ/Q43JpORC6KrcFlcZSLnGwpNtEl97rjwJOWWNm6qiYvzfoa2Y/vO30PFlibv3XovXpH17zsNSkrxxTfpCtikdDZgxvQxYP5U7mkroGKgtSocgCG7IQPEz1JhMSg7Erso14C7oxoXtWESHOmc35cKxnVo8V1JDTgPMJrxf8xMM776EljWVnG0q4vsg4p0NwLBhLvsg5rvFN+kLMfYAZ7E1a8C87WYM6uudvxeUU4Qg/AsyUPyMxuqmFrMqZ6ePGbd05IweccRuLITOPQ7thHquDv9eAq1Go6i4UnTI6d69wLx56HnyJOfh4uZRKHj6efR+bi6gdV8/VOx3y3HSF2LsLfryJGdxPlsDRqiQ1dELs+G+AUjeQTlFCELtkIHiZzRWN7WYVTk7mYyKN2DbsQsuk8YZwvV2Y8EV3cGFYzuh/Zv96c92205KRVYJEnHm5QFPPQX873+ch026QBQ+PAuG1SsQHREu+NpCk/Xx3bsQY4+vcrCjt8rdGKTnFmHFjtN214sK1WHF3T3RsnkQRekQhIpxv1wiVAXrpgacnQHu3NQmM4OMvBJsz7qIjLwSmMTkbVcYoavyxeN74OCCERjbK9blWLDUNpiRllMo+jqO7YS+z1ETw674U7N5srEqQWUlsHAhEB/Pa5xg0iQE/JqLNm+vQ4AI4wRw/R1kSUrsYv2cHPF0+9FWW+SOpM+ynIyh0qp6zNmWhf2nizCxXxskdGpJxglBqBAyUPwQVodgiLCfNA0RQbxCzdTsAgxbvRdT3zuMuduyMPW9wxi2eq93J04XsKtyvmlCA8uK/KGhHe0mE3YsIkJ0nO+rqK63MxCEXsfRA+XufXywJuDyb3KUNwjNZmDzZqBrV4vIta7OuU2fPpYtn//+F7jxRsmX4vsOxkYE4e0HBmBuYlfeSV+u7UdXhg471q5G/J0D+dh5Qh3ff4IgnKEtHj9FjA7B19EnQhBWMC6e835HxRuw7OscAM7bAlwCVrFCSVbDMK6XAZsPnXN6nzu8Ell18CAwdy7wM0+ekOhoS6Xhf/0LCAiQ5ZJS0697Us/JFleGDlcmWS4Wb8/GmF7+GZJPEI0dMlD8GCE6BLVEnwjBVYTKhL6xTsJGVuMQERyIwkrhoddiImG4QmE1muvp9sWgSGTV+fOWXCWff859vFkz4P/+D1i8GIiMlP3yUhKaCTESI0J0qKiul6yzEqo1Kqmq88uQfIJoCpCBolLkqrPjb3lTuFblZVVGzE45zusBmjG0g6Bz2xoIQlb/fJ4nqTs1skZWVVVZtnFefhmo5fl877oLePVVy5aPyhBSz8mTcGCuTLJ8+FtIPkE0FchAUSFy1tnxx7wptqtyk5nBsNV7XXqAtmddEnReRwPB1epfSrZUPmSNrDKbgZQUiwj2Ik+F3x49gLVrgTFjPL+egrgzEqXke2HhyiTLh7+F5BNEU4EMFJUht17E3/OmCPEAlVTVISpUh7Iq7i0BAIgM1sHMMDCZGUGeKCnZUrmQNQHYkSMWncmRI9zHW7QAli8HHnsM0HGLhtWGKyNRqsaFPa8QuATRBEGoA4riURHu9CKA+GiQgXEt4O5ZrdVY2qkRoZ6de/q1AcAf9lpeU4/73z8iOHJJLo+Sq8gqwVy8CEybBgwZwm2cBAQAc+YAZ84ATzzhN8YJF46h8ACQ0Kml5HDgh2+J4z2mAWWOJQg1Qx4UFaGEXiTzfJlbzYSZsbRTgwbFEaGencR4A27uGOW20KBQT5SnHiUNgLkjO+OJkfzhtm6pqbFoSFauBKqruduMGmXZzunZU3Jf1YKcW5ss80d3R+920Xh+ezZKq66HXSuVQI8gCPkgA0VFCF21F1bUICOvRJDb2x81KLaUVdVBq+EXptrqO9iQ48N5JZid8jNnIUGhkUuehsIyAF7fcxbdY8PFT4IMA3zxhSU65/x57jZduliMl7vusoQV+TlKhsLf2ScWY3pJ2yoiCMJ3kIGiIoSu2pN35ApeDfqzBiU1uwCzU9xXN7Z10wdoNdBqNaKrHDsiJBT2X7d2xKaD+S49VKJDuH/+GZg3D/jhB+7j4eHAkiWWrZzAQGHnVDneCIWXEg5NEIRvIQ2KihCardTWOAFcp1OXmjlVLqSm1xcSRaPVABvu649R8Qa7a7jKiWKLO6+Ru4y9I7rHuDROxKRkR1GRJYnaTTdxGycaDfDooxadyZNPNhrjBBC3tUkQRNOBPCgqwtWq3RWuVpm+LDHviaZASBSNmQHOXL6K5NV7nYrBCcHRa8SVe8ZVJMn2LJ4wXwdcGkJGIzp/+SWaTZsGXLnC3eb224F164B+/ZwOyZUvx5eofRuyMYwxQfgjog2UAwcO4OWXX0ZmZiYKCgrw1Vdf4W9/+5v1+EMPPYQtW7bYvWfMmDFITU21/l1aWoonnngC33zzDbRaLSZPnozXX38dYWFh0u+kkcCXwCoqVIfSKmnbFmIyp8qFp5oCoZPR2vQzTq+5Gifgum5lYFwLq5bnXHE1th69YOd9sTWmuLYHPNo+Yxhg+3Y0e+op9MzL435jhw7AK68AkyZx6kxSswuw7Oscuz4bwoOwbIK6xJ/uJng1b0MqIdwlCEIYog2Uqqoq9O3bFw8//DAmTZrE2Wbs2LH44IMPrH/r9Xq74/fffz8KCgqQlpaG+vp6zJgxA48++ihSUlLEdqdRwrVqL6ysRdJnWW7fyzexe5JTQixyaAqUmozYq03oG4vbX97nUcSPOyEtb4K2kyeBpCRgzx7urbfQUOC55yxtgrjHITW7AI994lx3p7CyFo998jPeVkF9JUDYBC95HBXGH2pYEURjRrQGZdy4cXjhhRdwzz338LbR6/UwGAzW/1q0uJ5jIzc3F6mpqXj//fcxePBgDBs2DG+++Sa2bduGS5eEZQRtCrCiPjb/gyHc81Wm4zmVclPLoSmQWj3YHRHBzXBXHwPeOZDvdgvJXe4ZdvsMcM6/wrl9VlwMPP64Zatmzx7ui06fDvz2G7BoEa9xYjIzWPjlSZd9X/TlSeWrJ7uBneAdx9lRMyV6HL2AEjmJCIIQhyIi2f3796NVq1bo1q0bZs2ahZKSEuuxjIwMREZG4qabbrK+lpiYCK1WiyN8GTIJn4tdxSCHpkDIpCWFytoGfHOiUHB7d8aUOyHt2F6xQH29RUPSpQuwcaMlXb0D5oQE4OhR4MMPgdatXfbp8O8lKK92vY1VVl2Pw7+XuGyjJGIneHYcY8Ltva0x4XqfeCpIuEsQvkd2kezYsWMxadIkdOzYEXl5eXj22Wcxbtw4ZGRkICAgAIWFhWjVqpV9J5o1Q1RUFAoLuScOo9EIo/F6ddLKykoAQH19Perrrz+o2X/bvtaYWDK+m3Wbh0vsumR8N5hNDTCb5L2u2HGNDmkGfYD7lWV0SDPUGuuQeb4MxVeNiA7TY2BcC+tKeWS3aLx1X1+s2vWrk85i8oC22LD/rIS7kcbliirU14dzHhvZLRrDu9zKeR8N33yDgKeegua33zjfa27TBj9PmYL4FSugCwy0GDNuOHL2sqDxPXL2MgbFRbhtpwRH80tRerUG+gD+NqVXa3D47GWrUc2YTQjUwu7eArWW18X+pj19FlyuqBI0xq6+F42Vxv6c9RVNZVzF3J+GYaQUjr/2Zo3GSSTryO+//45OnTohPT0dI0eOxEsvvYQtW7bg9OnTdu1atWqF5cuXY9asWU7nWLZsGZYvX+70ekpKCkJCQqR2nyAUI+zPP9Fr82bE/OysEwGAhsBAnL3nHpydNAkmB40WQRBEY6W6uhr33XcfKioqEB7u2rhXPMz4xhtvRHR0NM6ePYuRI0fCYDDg8uXLdm0aGhpQWloKg8HAeY5FixZh/vz51r8rKyvRrl07jB492u4G6+vrkZaWhlGjRkHnx/VI3GEyM7xeByWQMq7puUUuvT0zbonDBz+ed9oCYI+vndIPiT1iJF1DTjQAYsKDsHvebcLGuKwM2hdegHbjRmgaGjibmKdMAfPSS+jUrh3auxhbvs/5cF4J/vXxT2678vTornhgSAfFvhvpuUWc3q2F47ojPEiHh7ccc3uOzdNvxsC4Fhiz7gBv/hqxn4HJzODY73+h9LefENX1Jtx84w2ix8BkZjBm3QEUVfILd0V9LxoRTeU5622ayriyOyBCUNxA+fPPP1FSUoLYWMseckJCAsrLy5GZmYmBAwcCAPbu3Quz2YzBgwdznkOv1ztFAgGATqfj/CD5Xm8s6AAM7ep68haKmBwPYsZ1XJ+20GgDOCM4Fo+PR/KOHNSauK+jAbBix2mM7tXG5cOf7xpywV550fieCNK7SYzW0AC8+64ly2sJj/bjppuA11+H9pZbnMRfjmPrKvplVLwBwfpAtzqUF3adwfuHLmDZhJ6yazhSswvweMov1ybv65/RhTIjHk/5BRvuG4CosGC3kTlDOrfC0fxSnC8zwpW66HyZEcf/vOI2Gyw7bqVXa7BmEPCvT44jKixYdFiwDpbPfda1SCkuI1vQ96IR09ifs76isY+rmHsTbaBcvXoVZ89e3/vPz89HVlYWoqKiEBUVheXLl2Py5MkwGAzIy8vDM888g86dO2PMmDEAgB49emDs2LGYOXMm3n77bdTX12POnDm499570dqNOLAx4sskUErneOALbZazKKLjNdJyivDtCffVioUgOE/Mnj2W9PTZ2TwnMgCrVlkqEmu5delH80tRXN2AVs2DUFZlxOyU4y7DW1dN6s0ZZuxIYaURj33yM5ISu2DOiC6yfLeEhJEn78jB4vE9MDvluNsEgXIlarMNC7bVvtiOm5hQe1/kDyII4jqiDZSffvoJd9xxh/Vvdutl+vTp2LhxI06cOIEtW7agvLwcrVu3xujRo5GcnGznAfn0008xZ84cjBw50pqo7Y033pDhdvwLXyaB8laOB64aKHJnDmWvYTIzWLXrV9F9ZIlpHoj7BsehQ3SoMGPx7FngqaeA7du5j+v1wPz5wKJFMIWGcU6M6blFAICHtxyD8ZpHSavh3rKyzSFzcMEIvP3AACz7+hQKK40cre1Zm34GW4/+IUsSN6EGZotQvaAJXkyiNj6DXojRtPDLk06J7dz93uTKH0TZaAlCPKINlOHDh8OVrnb37t1uzxEVFdXkk7L5MgmUN4qzuULohFR8xQiTmRHcByHp8flISuyKOSM6C7tWZSXwwguW0GE+RfrkycDLLwMdO14zRI86GaIT+sZiy6HfsXqQ/VuF1vdhJ8/NB/Px4s5ct90urJTnuyXGwJzYr43bCV5oorayKiOGOZQ1YA2MiOBAt0aTZUvM/vMS8nvztNCgWrPRktFEqB2qxeMDfG0gyLnFIgV3ExJL8o5cvH8wX/CDXEqtFlEThclkyVPy7LOAg9DbSt++FsNl+HAA/IZoQUUt3jmQ7zIM1xXsvQZoNaioqXPT2p5lX59C8yAdiq8aJU1MYlPTu5vghdSLmtA31uW218NDOwjtvh1K/97Umo1WrUYTQdhC1Yx9gK+TQHm6xSK1QjGLqyRsjriq1OyI0Ilz8fgeeP3eftg6cwgOLhgh7IH8ww/AzTdbKg5zGSc33AC88w6QmWk1ToRUZJaK/b0Kn1QZWHQp979/BHO3ZWHqe4cxbPVeQePLokTSQFcJ7zbcNwBf/1LgMunbVwILN3Kh1O9NrdlohWb4JQhfQx4UL+DoSuULp3REqeqtnhRnk2vlxSdAdETMClfoVsFDQzsKXymfOwc88wzwxRfcx3U64P/+D1i8GIiwT4rmyZaTKyJDdHaTf0Knlli/T3rSOrGreaUqZHPpPQbGtcDHGefcfkdKq+oRFRqIsqo6yQah3L83X3squfC195YgxEAeFIVJzS7AsNV7MfW9w9YVa/K3pwS9V6mCeVJXwOm5RbKuvMb2isXBBSOweHwPl+2ErnBlrely9arF6Ojend84uftu4NQpS8VhB+MEAAoratxfRwLl1fVIy7medXnIjS0RGSI9LFHKal5Qin8J2NaLqqipw+0v70PyDvf6GgD4Wz9LFKDUaZUV4XriHbRFbjG4HPjae0sQYiAPioLw7T+XVrnOX6F09dYArQaLx/fA4ynHOa8NcE/kq3b9qsjKq8xNPg8WIQ9yj0NDzWbg00+BhQsBvuKV8fHA2rXA6NEuT1VaJU4bIhTHsQ7QagSHHfMhZTWvZIVsvt+OK0bFGzCoY5Q1DwqLIVyP2gYzKqrr3Yhw63hFuFIMLk88lUqhRqOJIPggA0UhhOoP5HSRCyU1u4B3VepqIrdsTXH3ScoEx7Vd5AqhD3LJE+fhw8DcuZaifVxERQHLlwOPPQY04/7p2G7nKWWgcI312F6x18KOcwRvIXIhdmLyNMKFC7HaHVuDPkCrwah4Aw6fvYzi3MPYPP1mDOncCmk5hQJEuPKKWYVuOXqzwKcajSaC4IMMFIUQqj8I1TfDVeP1lOhKJ4FytzJdPN6zawud4MSskKU8yEVNnH/+afGYfPopz8kCgMcfB5YtsxgpPIg1uPiIDNG5zRILOI+1rWGWnlOIr7IuuvXWOaKGiUmMdofLoA/QajCoYxR25sJqtLjyrC0e3wPJO3Jl9w4qpdXxBDUaTQTBBxkoCiF0or5qbEBSYld0iA5RPBeBu5UpmwF0TC/pAjkhE5yYFbKiD/Lqaot+ZPVqy7+5GD3asp0TH+/yVFK2JABL7RqgCpun32zNJGs2M7h/0xG37+Uaa9YwS+jUEs+Oj7d6c6JD9Xjyi19c1pZRy8QkxosjxqD3NLPx4d9LMLRztOC+sddUUzZaNRpNBMEHGSgKIWYluu3YBRxcMELxh4KnUQWG8CBcKDN6PMGJWSEr8iBnGODzzy3RORcucDb5vUVrvHX340h8agbGxrsuwSAlnJj9pBeO6466/EwM6hhlrVFhMjOyrHIdvUjLJvjHxCQmXFxURBY8y2w8+9OfsWpyb9HfRSW1OlJQm9FEEHxQFI9CsK5UIXhLNe+pQG7huO4API+QEdqPOXd0Ep6nxA1sdMb+lF2oHHQLcO+9nMZJpT4UyXc8gjGPbMB/Y/th1qfHXUYnmcwMPjyUL3pbh4124arY7C5PDAPL5CwljFeJyBu5ERplJtY44UOoQVReUy85TwhrGN3Vx2LsfnviksdRQp7ARtBtnTlEfE4ggvAS5EFRCHaSERpZ4Q3VvKcCucQeMbKsvIT2Y2jnG2SLCHkz5SCmf/su/n5yD7QcfgmTRottfUfj1VunoTTkesiwK/2BWM3JnDs6oUtMc7sVdD1Pqnx3eWKSd+RCe01bIQa1rea58PY2hNDMxixSo9XUlr1VCYEzQcgJGSgKMrZXLJISu2Bt+hm3bb0hTpRDICfHBDeoYxQM4XreIndy6iF2/3wOWfOX47OMzxBWx52XJKN9b6wYORO5rW50Osa37SVFczK08w2CJgQ2EsjYYMaUm9ph3R7n748n0SX+MDG52oa49+b2MDaYkZFXIotxZWsQuUNqcjW1prwnCDVDBorCzBnRBVuP/sEb+ulNcaLUlSnrht55sgCtIkIxqGOURxNcWk4hahvMnMdkWyEzDExffoWej87BmFJul/yFiBisG/MovuwwCNC4vpath8uTMFh3CPXKeBJd4i9F4hyN4XPFVdh69ALWpv9mbSOXB4I1iBb+9yTKa8RHULlCrdlb/eV7QDRdyEBRmACtxipOBHwvThQrkEvNLsDKHacwvzvwzH9PwGjSeDQp7DxRgMdT+FeqESE6rJokXohox4kTwLx5CNi3D205Dl8NDMaGhH9i800TYWwWKOiUth4uT8Ng+RDrlWFX8x8eyhesx1DbNoM7WG9PanYB1qWfUdQDMbZXLJrrdZIjqPhQY8p7f/seEE0TEsl6AbWJE4UK5NgJ09H7IzW1/c4TlzB7q2s3ugaWjKCS+OsvSxK1/v2BffucDpuhwee9E3HHzHewccg/rMZJZLBOVNp/MavniGAd5iV2cXtPnhQWTN6RK6jgn78WifNm0b0hnVrKXghRbdlb/fV7QDQ9yEDxEmpTzdvWPEno1JJzW0fOSSE1uwCPpxwH46Z5WXU9Nh/MFzfZ1NVZcpV06WKpKGx23j461iYeEx98Dc/cOQ9/hdlPLjOGdgQgPDpJzOq5vKYea9PPuDUgMs+XeZTgzd3kotbKukLwZv0YWes5XUNN2Vv9+XtAND3IQPEi7owCNSHnpGAyM1j2tbACiQDw4k5hHgEwDLBjB9C7NzB/PlBR4dSkMOIGPHH30/jH/atxMraL3TF2NTxnRGdOD1dEiA7zErs6eT/chcFy4c6AKL7KLRgWirvJxZ+LxHnbAyG3x1NqcU4l8OfvAdH0IAOF4ETOSWH93rO8ETt8uHU35+YC48YBd90F/Pab8/HgYGD5cpxIy8C38bdD4yCCdVwNsx6upMQuiAy2JEwrr67H2vTfnIwld3lKuHBnQESH6QWeyfU1+CYXtW0ziMEXHgg5PZ5KeGWk4s/fA6LpQQYKwYlck0JqdoFd1IVQeCf00lJLQb/evYHdu7nffP/9FqNlyRKMvrmT4NVwWk4h1qWfcYri4DKW+FbZ7u6Jz4Aor66HXPMT1+Ti6efJJrrbnnWRN8GYkDZS8JUHQk6Pp1p0aGrabiIId1AUD8GJbc4ULoSEzrL73VKxi26Ii7DoS5YssRgpXNx8M/D660BCgt3LQnK3SAkFdTzvmaIrWL8vz+19cRkQT36RBTPjegL8+4A2+M/PF92en2tykZoDx2RmsH7vWXxwKN/OcHOM+FAyKsTf6sfwhe+qIUkeFQsk/AkyUAhOHCcFW4ROCmLCcV3BpH0HrH8ROMWjY4mNBVatAh54ANByOwXdJSeTGgpqe96MvBJBBoqtAcF6GVz5GrQaYP3U/hjTKxaH8kokTS5SJvnU7AIs/PIkZ2Vl2/BeAIonIfOX+jHuDDVfJ8nzN2OPaNqQgULwwk4KK3ecAlBlfV3opODpPnaH0ot4bt9m3HKWJy+FXg88+SSwaBEQFubRteTYm5eyOs08X+b2mmYGaBGq93hyETPJp2YXuCzTwHqVLOJnjVeSkKnBA+EKf8kW6y/GHkGQgUK4ZGyvWAzv0hK7U3dhzeQ+1kyyQiYFqfvYzY1VmPPjZ5jx09cINDdwtikZOwGR69cioJMlPb2nWTHl2JuXYkAIjd5hDSMxkwvXmIjZ7nIHA7gVP8udhMzXHgg+1Jotlg+1G3sEAZCBQgiAfWjd2TsWOp1O8PuEFGFrEaJD2bUtBK3ZhH+eSMNTP3yM6GrnkGEAONXqRqwYORNH2vdG7Of5WHp3MAB4rH+Qa29e7Oo0OkyPYgH9szWMhEwu7rYaPNnukkJjjwpRY7ZYd6jV2CMIFjJQvEBTrXkhxKOwclJvnC68gh83f4kle95Dz8u/c57rr5BIvHLbNHzROxFmbQAAi+ucbxtCrFtdzr35sb1iMaJ7DD7OOIfzpdWIiwrBtIQOCGzmrI8ZGNcCu3P5w5X5DCNXk4unWw1KGBONPSqEwncJQn7IQFEYb9a8UKMh5NajEFqL3qvnY276Ds7312mb4YObJmD9LVNwRR9qd8yVsFSKW12uvXmuz/z9g/nWc9h+TtEh13+CcogW5dhqEGNMGML1ADQoqmzaUSEUvksQ8kMGioJ4UzTni+JfQg0izi2JGwIRsHoV8OqraGPk1jGkdR6MF+94GOei2kjqnxS3uqd78+4+80dv64ivfymwfk76AAZrBgEzbonD/05c9li0KMdWg5CtOeCaSHZCTwBo8lEhFL5LEPJDBopCeFM054voAbEGkXVLwmwGPv4YzKJFQAF3ltjT0e2RPGImDnbsL0tfxbrVpe7NC6lz8s6BfM73fvDjebw2ZQBahOo98oCJ2WrgMzBdbXextAjRYaVN1emmHhVC4bsEIT9koCiEt0RzvogekGwQZWRYssAeO8apuSgLao7Xbn0AKf3GwnRNZyIH3nKreyIuZWCpSnxwwQiPPieh93quuBrDVu/lNTD5trsig3WYMbQD5ozoYtdPigqh8F2CkBsyUBTCW6I5b0YPmMwMjp4pxsL/nhRnEP3xB7BwIZCSwnneBo0WHw8Yj3VD74OmZRRM1fYZS2vqTaiorne53cAFl1tdSZ2Op5+lJ58Te1+FlbWICtWhtMo5uRpgGZOIEB3Wpf/m1sAUa3RQVAgZagQhJ2SgKIS3RHPejB4Ys+4AzpeJyHsRGwy8/DKwejVQU8PZfn/HgUge8S/kRbeDBkBMMy0+/ddgFF81Wh/uaTmFvK5zhuPf7N+AvVtdaZ2OHJ6awgrucXIF131xYTtGQg1MMjrEQ2NGEPJAxQIVwlsFzrxhCKXnFgEACisFGjkMA90XnwHduwPLlnEaJ3lRbfHQ35fioX8uR150O8vbYEn8pdVo7Aq0uSq09vYDA/A2x7GIYB3mJXbBqHgDgOvbUo6TuNuqySJw95kLobSqTlR7vvviwhARhKTELpyp61lcFTRsqihVBJEgCNeQB0UhvCWa86QInBA3tMnMYNWuXzG/u7D+9C44gyV73sNNF7kzkVbqQ/H60Kn4aMB41AdwJ33j8va4c52PijfYFbUrr6nH2vQz2HbsDyweH4/kHcrrdISIS90RFabnfJ3r88K1fvNdRwMgKjQQz4/vAUNEMAZ1jMK3Jy4J6gfl67DgyutGWzkEoSxkoCiIN0RzUovACd3qYHUN7rjhaime+f4j/CM7nbuBVovCex/EnS1HoTQkwuW5pHh70nIKeXUVj6fw15QB5NXp8H3mYfoAXDWa3L7fEO5873yf1703t3OrPyqpqoMhIth6X5SvQziuxOCPffIzIkN0dt4opcP6CaKpQQaKwnhDNCe2CJyYCBx3K2l9Qx0eOfY/PH74C4TV8egnRowA1q7FDb16Q796LzQSckW4W8m6C+8VglxeA8fP/FxxFdamn3H7PsctP5OZwfq9Z7E2/TentoUVtYLOCdjfF+XrEIaQkHHHrTK1FQUkCH+HDBQv4A3RnJgicGK2OnhX0gyDMb9l4Ll9m9C+ooi7zY03Aq++CkycCGg0CAAkbXu5M6rmJXaRpXaMnF4D9jM3mRkMW73XbXsNnAW9y74+xVuMT4zhZXtflK9DGFJCxtVYFJAg/BkSyTYi2EnRVmBqi5iQZJZBHaOcth16XP4dW7c9i3f+9xK3cRIWBqxaBeTkAH/7G6Cxz5fBJ3jlWnkKWcl+cOgc7z0JQS7BMhdCJ7rHh3e23jtrkLmrFOwOvvsS+xk0RaR600hkTBDyQR4UDtRY00YOpIQkB2g1WDiuO+ryMxFVVYEn9n+CKSe+QwBjdn6jRgPMmAG8+CJgMPCeX+i2l8nM4MND+W6NqvIa/qgUITCweA3Y67kr8CcGoWMe1zIEgGuDzBVivSGUr8M13gr/JwiCHzJQHPBFTRtvIVUgmdipBU6v/R/Stn6OMGM195uGDQPWrQMGDgTg3shzt+0lNLcHS2SwDhU14pO5AUBSYhccv1CGxz/9GbYRpC/uzMXMWzti0Z3xEs5qQeiYR1+L3pGytZCU2BXbjl0QLcT2p3wd3l40CK1HxAeJjAnCc8hAscEXNW28iWiBJMMAO3agWVISep09y33S9u1hXrUaR25OxOWrRrTKK0FZVR2Sdwgz8rgmHjYxm5iJYcbQjliX/puk8N4Tf5Zjz69/Ob1uZq7XzpFqpAgZcwAYGNcCgLiVN/t5zRnRGXNGdG603hBfLBqkhoyTyJgg5IMMlGuIFZC6WtGpdYtIlEAyJwdISgK++4478VhICLBwIb4bdz+WpuWj4JcjLq/NZeRxTTyG8CDUNphETwhzRnRGN0OYKK8Ly14O48SW937Ix5Oju0va7hEy5mw7QPzK23YLx1+8IWLw5aKBtx7RtfBiEhkThLKQgXINMQLSipo63hUdAFVvEbkNSW6tB554Ati4ETDx5O144AFg5UqklgcI9nQ4Gnl8XhLB2WptYCcER11FdJgeT36ehaJKo8s+uuu/mQE+zjiHR2690emYEGPU1ZgvGd8NdfmZ1teEbi2o6TulFL4ohOkIn1YnLaeQigIShMKINlAOHDiAl19+GZmZmSgoKMBXX32Fv/3tb9bjDMNg6dKleO+991BeXo6hQ4di48aN6NKli7VNaWkpnnjiCXzzzTfQarWYPHkyXn/9dYSFhclyU1IQ6lpPzynE5kPneJM3ceHtLSJ3kybnQ7ddOALefQdYsgQoK+M8r/nmm6F94w1gyBDL5PHxXlHbKayRdzivRJIQ1BGtBlg/1X5MHXUVyyb09CizK8v5UmftjZitB76JzmxqwM786+2EbC0kJXZxqibcGBG6aFib9huGdo5WzFPJpdUhkTFBKI9oA6Wqqgp9+/bFww8/jEmTJjkdX7NmDd544w1s2bIFHTt2xOLFizFmzBjk5OQgKMjivr7//vtRUFCAtLQ01NfXY8aMGXj00UeRwlPt1hsIda1/lXVRdEIwb+ZHEDpp2j10v/sOmJBk2dbh6n/r1vj5H/9An9WrodVLF3OyZPxeLEveEjMDtAgNBMBvlPF5L8QSFxVi97eUrQeuic7M4aTi63NT8JrYInTRsH7fWazfd9br4+NPImOC8EdEGyjjxo3DuHHjOI8xDIN169bh+eefx8SJEwEAH330EWJiYvC///0P9957L3Jzc5Gamopjx47hpptuAgC8+eabuPPOO/HKK6+gdevWHtyOdISIGaNCA1Eispgbi5zp1PkQPWn+9hvw5JPAt99ynzAoCHjqKTTMn48/DxxAH+11DYZnYZTyGWiXr9S6NcpsV7s/nLmMt/b/LuoaWg0wLaGD9e+6BjOe/Spb0a0HWqGL1+PYfs9HdotWqFcEQXgLWTUo+fn5KCwsRGJiovW1iIgIDB48GBkZGbj33nuRkZGByMhIq3ECAImJidBqtThy5Ajuuecep/MajUYYjdeTVlVWVgIA6uvrUV9/PQcG+2/b18SwZHw3JH2WBYBb/HZPXwM+PnJe0rlZLldUob4+3KNzmMwMMs+X4fIVI8qqjGgREojoMD1e+vYUAgO4fTkaACt3nMLwLi0RUFkB7UsvQbthAzQ8Y2WePBmmlSuBDh2s42msq8PR/FIUXzWi5KoRep5ruUKrAXQak6T3cnH+ryt4a/9ZMAD0AddfL7tag3lbM7F2Sj8k9ogBAJRX1eC/xy6IvnaILgB7ci4hsUcM0nOLsPybU6iqrbe7niOlV2tw+Oxlt9Ec7r6zN7UPB2D5vphNDZweF3ew35fiq0ZEh+kxMK6FXxg6/ds2R1wLPYoqhYf6st/zYR0TAEh/FhD8ePqcJbhpKuMq5v40DMNInik0Go2dBuXHH3/E0KFDcenSJcTGXl+p//Of/4RGo8Fnn32Gl156CVu2bMHp06ftztWqVSssX74cs2bNcrrOsmXLsHz5cqfXU1JSEBIS4vQ6wYPJhLj0dPRISYG+ooKzSfmNNyL7kUdQ0rOnlztHEARBNHaqq6tx3333oaKiAuHhrhfrfhHFs2jRIsyfP9/6d2VlJdq1a4fRo0fb3WB9fT3S0tIwatQo6HQ60dexeiYqa1FWXYcWoXq0an59xWkyMxiz7oCoFR2LBkBMeBB2z7tN8uo1PbcISZ9lSRJ7Djp/Es+mvYcel/M5jzOtWsGUnIzQBx/E4AB710DaqUuoP5+FxT9pYTR7vvLWAIi4llgNcPZWMQBmD++Eiup6fHL0Auf7AeDx4Z2wYX+e2+u9P+0mPL89W1KEkOM1xYz97OGdMWt4J5dtPP3OuoLv+8Lei613Sc2k5xZh1a5fRX1+a+7pCVw8oci48vXJEB6EheO6+8WYeoKS39mmTFMZV3YHRAiyGiiGa+nNi4qK7DwoRUVF6Nevn7XN5cuX7d7X0NCA0tJS6/sd0ev10F8TZ9qi0+k4P0i+113hSscQpLcIMXUAFo23RIUA3BOr47/Zv3Htvey5xGIyM1ix4zRqTeIMhHblhXh232aM++1H7gY6HZCUBM1zz6EZhzVrMjNYvfsM5ncHjGYNjBzXd7zfqBAdSqtdu/Fqrzbg7wPa4uDZEruHfGSI5XN7bQ+rE9FAq4Fdhlf2czE2mDn748iR8+U4X2aEnNoXIazdk4eusRGCRJtSvrOucPd90QBYseM0Rvdqo/rtnnF92mJ0rzY4ml+KQ2f/wvp97o3S6PAQFF+Uf1wBy7Pi8ZRfrn3nr4/dhTIjHk/5xe8TOgpFibElGv+4irk3WYsFduzYEQaDAXv27LG+VllZiSNHjiAhwbInnJCQgPLycmRmXs/9sHfvXpjNZgwePFjO7giGFZc6RnmworvU7ALra64Krb39wAC8rVARNrFRM6HGajz9/Rakv/8Yv3EycaIlcmf1aoDH1XY0v9TtytVxhW4S6Gb4z88XATBISuyC1+/th6TErqiorncqY89uQj48tAO2zhyCgwtGYGyvWMEiyotlNcI6JDMMgIVfnsShM8UwmeXR3AhFSmFINcNGzCSN6obYiCBeU5MtkMhm5pUbIcUrl3+T4/XPmyAaI6I9KFevXsVZm7Tn+fn5yMrKQlRUFNq3b4958+bhhRdeQJcuXaxhxq1bt7bqVHr06IGxY8di5syZePvtt1FfX485c+bg3nvv9UkEj5RkUO4iLJSIvhAaNaNhzJicvRfPfL8Fraq485mgVy9L3ZyRI2W7ri0VIor3FVUasS79DDbc1x/bjl1w+eD/8vhFjOx+3X0uNKnZV1mXBPdHbsqr63H/piNeD4GVUhjSHxCVDdkD+MLWxRh+FIJMEJ4h2kD56aefcMcdd1j/ZrUh06dPx4cffohnnnkGVVVVePTRR1FeXo5hw4YhNTXVmgMFAD799FPMmTMHI0eOtCZqe+ONN2S4HfGIfeA4Prju6tPa6WGoRH4EId6CgX/mYMme99C38Ax3g5YtgeRkYOZMoJmwj17pomesEfj89myUVrk2bLgmeyn1UhzRwJJPxWQyo6K2QeJZXOPtZH1SC0N6E6klIdxmQ+4V61EkhKvtXmMDRxVvDvzN8CMINSLaQBk+fDhcBf5oNBqsWLECK1as4G0TFRXl06RstohZafqy0rErb0HrystYuP9DTMg9wPleplkzaGbPBpYuBVqIc30P6hgFQ3gQgCppHRcAA7g1TmxxnOw9TcRmub60/DZiruGtZH0AUFZldNLu2OLronae/paUyhPjLpfQvMSugs5D1YwJwnNk1aD4I0IfJOeKqwXrVJSAdW3bElRfi3kHP8We92bxGif7bhyIqXPeReq/Fog2TtjrLhzXHYA0iakS07DjXv/YXrE4uGAE/j6gjQJXkw9v6T5SswswO+U4r3HC4quidmI0X65gPZUT+7VBQqeWsmzruNOXbDt2AYZwvVsNDFUzJgjPafIGCuuZcPXAMYTrsfWoa32EN4RxrLcgNlyPCTn7sfe9xzDv0FYENxid2uZFtcVDf1+GGf9YjiP6Vpj1yc/YeeISMvJKsD3rIjLySgT3lw2bjAkXvyoUMyJRoYGCDRp2sv/wUD5MZgZpOYXXRLfumXNHJ6yd0g9RofIp5cVMjUq6/11NsixaDbDhPt9EmqhZZCp0u3fqoPYAnD9zqmZMEPLiF3lQlESI6G7qoPZYm86j64B3hXFja/7EmG+XQJORwXm8Qh+KdcPuw8f9x6MhoJm1fwAwZ+txznBdVxMVqxMAgBcm9oK2WTMUXzXicmUtXtz5qyz3xPZl8fh4zE4RpydJ3pGL937IR22D8BSrQzvfAEDctpIr/n1bR3z9S4HgLaZzxc6FB+VCSLSXbQ0jb6NmkalQw7FDdKhbDQxBEJ7T5A0UwL3oThXCuEuXgEWLgI8+4lytM1otPuk7Fq8Nux9lIRGcp3BclLoTbrI6gdKrNVgzCPjXxz8hKiwYS++ORysJ3hRXsA/2jVrxehIxCbxY9/u3J+SJ7ElK7IK5iV3xzNgeOPx7CWZ/+jPK3UQxrUv/Dd0MYYpMZGqP3lFz/8QIixM6tWzytZIIQmnIQLmGK9FdRl6JoHMUXzFie9ZFeR9WtbXAa68BL70EVPEIVUeOxL5/L8TiTOetHle4Em7aigVta85cFwt2EXSNqNBAlFXV8XpEtBpg/dT+1sl6bK9YmM0MHk85LupehMK63+UQMcZGBGHOCMs4BGg1GNo5Gqsm98Zj1xL5uYIdc1ukRrXYovboHaHXjQ51TsyoNEIKhtoKi6maMUEoS5PXoNjCJ7pzp1MBLBNt8o5czN2WhanvHcaw1Xs9E84yDPCf/wA9egDPPcdtnHTqBPzvf0BaGoIH9Jd2GTgLN4XoBLYevQBDuPuEWS9M7GX9m4v1Uwfgzj7X89+YzAxWfJsr8i6E8fcBbaxGgRDtUWSIDhpwaw004NYajO0ViyQ3xhvXmKfnFmHY6r2Y+t5hj75DQu7LVyJOk5mB2cwgMti99ufJL35RXHjuCLvdy2dMMyB9CUF4EzJQBGAbQcP3aOLbPpH0kM3KAu64A/jHP4Bz55yPN28OrF4N08lsZPQahu2/XILZzLg0GNxRWFlrFdB+eCjfrU6gsNIoSCx4Zx/uzLux1zLv3tnHfptj/d4zHtXMccV/fr5onfRdfabs3y/9rRfmJXZFhMOE6i4zcIfoUEH9sd3GSPosS5YIMSH35YtJNjW7AMNW78X9m4643QIDgKJK70THEQShXmiLRyB8OhW+XBOS8l5cvgw8/zzw/vvX87vbotEADz8MvPgiUv8yY/m6Q3Z9iQzRWa8rNgYi+dtTokWjHaJDBIkFheasSM0ucClGlgMhOVQMEUGY0DcWyTty7cc3WIcZQztgzoguLj9PMdssbLSKmEzG7hCSyMyb8OUWcYW388YA1z2HfHi7PwTR1CEDRQSOE23xFSOSd/BvRwiOSKirA954w5Llla/S4623WtLTDxjA+8CvuFbDJiJE51TPxh1SIlrEiAXd7de7mxzkwnHi4zKeyqrqMDuFY3xr6rEu/Qy6GZpb31NYUYPSqjpEhelhCLfcuxgtw+Gzlzla2PdXSlSLUonMxCIk7JkPb0f0qDnCiCCaImSgiMR2ot2eJSzvBm9EAsMA334LzJ8P2NQ3siMuDnj5ZeDvfwc0GkG1g4KaafHcnT3w4k5ltBxKiAXFFkP0BFcTjdnMYMW3rsd30ZcnsezrUyisdBYls6HbQuvFFF8VJmyWEtWiBhGnHJ+rtyJ61BxhRBBNETJQPMCjiIlTp4CkJCAtjftNISGWsOInnwSCg60vH84rEaQPiY8Nd1tIz1UqdD6U0jF48tCPDNYJ0jVwXZMr5borGABlLrxTBTZbSEK2WaLD9CgWcF1/TZ0ux2TurXtXewQUQTQ1yEDxALFhiQCAkhJLTZy33wZMPMnFpk0DVq4E2tinbk/NLsDC/54U1LfiKqPbQnoPJsThwx/PCzofi1Qdg7sQWqEP/b8PaIuhnVvabauYGQb3v39EVH8A4FxxFdaln5FcZNAVy7/JwcEFI9xuswyMa4Hdufzia1/XzPEUTyZzb9+7pN+zDXKEiRMEcR0yUDxAVOn3+npg40Zg2TKgrIzzfFf63YTj85dAd8stGBQbBZv0I6KFhqw+hGsVz25DRAQHCjJQFozpBpTlYPP0mzGkcyvrQ1foA1lIYTh3kwNgGdMR3Vs5Rf6YzIzb9zqexxARhJQj3OULPMVxC8nVNovteLn9DvkhQj5XQB33Lur37IAvC4kSRGOFwow9hI2YcAyjtQtFTU0F+vYF5s7lNE5qW8Vi6T8Xoc/oJXjwlNYpB4YYoaFjngu2kN7WmUPw+r39sHXmEBxcMAJje8UKzplx3+A4ALAzQNiwUXd5O4QWhuMqhugIA2B2inPoqZAwcNt7Aiyei6Ir4hLbiUXM9sbaKf1cf4f8FHdhzxpYSgWo5d4F/Z4dkKv4IUEQ9pAHRQZ4IybO/AbcNRPYsYP7jUFBODv9MUwIGYbqQPsHom04bERwoCihoeMqj08s6W7FyAAY18uAzPP2RpW7kvQbHxiAUfEGHM4rwcL/nhQcQju2Vyw23NffqWaQI1yhnkLDwNkQ4ncO5PNfQCbEbG8k9ojB6F5tGuUWgZCw52fG9lDNvYuJgBIiWqfQZIKQBhkoMmFnBJSXA089Cbz5JtDQwP2GKVNgWrkK07bloZrD+LB9uD01qqugPkQENcPqv/cRterkmzw0GkuQ0eZD5/Dp4XysGWTJdjq6Vxu3D2RXUS6O7R2jaVqE6l0aJ3wROCYzg4jgQDwzppudPmVgXAtkni+zTjQD41rg9pf3CR4fW2IjglBTb3Ibwi1VO6GGqBulcDfpq+3ehfaHQpMJQjnIQJETkwl47z1g8WKgmCc2Y8AASz6TW2/FUQEROQUVtcj6s1zQ5eeM6CLJJW47eaTnFGLToXOcRkLSZ1mYdbnabZ9dRblwYbsVIiXU09X+f2Azrd3EkOFmzG0xhOvx6j/7ofiq0TqhpuUUCtICeaKdaKxiS7UZIXJAockEoRxkoMjF3r1g5s2D5iRPlE1MjKXg30MPAVqL9EfoQ0to6vfo5vYF1sRMdAFaDQZ1jML8z7NcXuODQ+cE9UUMtlshYkM9hWw32RptYiaKZRN6YmjnaLvX+DxOLJ4KI0ls6V9QaDJBKAcZKJ7y++/AU08BX33FKdA06wKhnZ8EPPssEB5ud0zoQystx3W2URZD+PXzSZnohLirpeQb4YNrK0RMqKeQoobzP/8Ff5TWYPotHRDYTCt4zJMSu/KOk63HiSuTrFRvh1hji/A9noYmEwTBD0XxSOXKFWDhQku14a++4myS2jUBd8zYgNQH5joZJ4CwKslCsY3ckRpV4As3NJegV2ixOyFZSqvrTHhxZy66L96FlTtzBI25IVyPOSM6uzwvu11xz4C2eOTWG3FPf/sK2GIRYmwt/ybHWrunMWAyM9YClRl5JX55b2otzkgQjQEyUMRiNgObNwNdugCrV1vq6DiQe0MHTL33RTx2z3O40CKWd2IREx7LBxuqyT4E3U10DPgnOm+6oTUaYMN9/Tk9AnyhnjHhesxL7AJjgxkZeSUorKgRfD0zA7xzIB9rUnPdhr0um9DT6xNK5vkywWJLT1GDYSA0TN0fkBKaTBCEe2iLRwwHDwLz5gGZmZyHS4PD8eqtD2Bb3zEwaS1p1typ+N1pGtzhmNlViFehoKIWa9NOY2jnG+y2JIQm1ZIDhgEiQgJ5jztGfZwrrsbWoxfsqh1HhfK/n4/3fsjHr8ndvVbtV6gOSMmaPLaoQeOihq0suYXIainOSBCNCTJQhHDhAvDMM8Bnn3EertcGYMuAu/DG0KmoDArjbONqYnF8uJ0puoL1+/LcdmvOHZ2RNKqr3UNQqKB2/b48rN+XZzc5ucqLogQZeSVOIlRb2G2U1OwCrEv/zak/pVXO3it3mBng44xzeOTWGxWfUMQYA9Fhese3c+KJl0sthoGv84YoZaQ1xiglgvAltMXjiqoqYMkSoFs3XuOkbHgixjy8AS+MnMlrnADuJxb24TaxXxsM7XyDoO4N7Rzt9BAvFbgSZ3HUpfC5q5XBvQkkJouuUM6XVgOwH3NP9CNciNUBDYxr4VYb0zwoAJfKayRty6hF4yImb4gSUNZXgvAfyEDhgmGATz+1GCbJyUAtxwO1e3dg1y6E7/kONTd2dpsuXoyKX2gKeq5zit324JqcHNPjb55+MwD7KCE5SLiR33vCImTLSixxUSGyns8RKcaAED3SlVoTnvziF0l6DV8bBiy+zBuiFiONIAhhkIHiyNGjwC23AA88AFy86Hw8MhJ4/XXgxAlg7FhFVPxiz2krepSy7cE1Odl6F1hDaOG47qLPzUdkiA5DeNzhtvfzUcY52a4JWFLfT0voIOs5HZFqDIjxXhWIXPGrJaGYL/OGqMVIIwhCGKRBYbl4EVi0CPj4Y+7jWi0waxawfDnQ0n5iFVJrRCxCz8m1n+5Yf0Yo7ianxB4xSErsYidU9YS0nEKnseG6HzmZeWtHBDZT1i4XZwzYh5+P7RWLEd1jMGTlHrfGJhuRJUSvoZaEYr7MG6IWI40gCGGQgQJYcpr07s1ZaRgAkJgIrF0L9OrFewq5VfzW2jJju6P0qhFRoYEwRAQ7VRTmEj1K9VALmZzmjOiCrUf/ECzG5aOiut5JmMl3P+5oEaITlF4/sUcrLLrTdcVkORA6yfMJYzPPlwn2hAmt86KWhGLuClQCyuUNUYuRRhCEMGiLBwCaNwceftj59c6dge3bge++c2mcsMglurTNEZH0WRaSd+Rize7TqKips9vWcSceFXN1rQYoEzApBmg1WDYh3pozRCqOe/6eiGGnC9yySc+97BURpNAEfE9+noX03CKn18Wu4IW0V1NCMV/lDfFE20UQhPchA4Vl8WLghmvRM+HhwMsvA9nZwIQJlqxiXkJolIEQ8aiYyd7MALNThGka+CYYsXOb7Z6/VDHsI0M7oMFsFtxeDhGku0RntsaAK4oqjUj6LMvpdbEreKHt1ZRQzFGIvXXmEBxcMELRPqjJSCMIwj20xcMSEQGsXAkcPgy88IKluJ+XEZMjQol9cgbAsq9PCdI0OG5pFV8xInlHrqTrenIvifEGHDrLUzmaA6FbInxwaWSiQnV4YWIv3NmntfW1sb1i8ehtHfHeD/m8W262L5vMDHTX/i0mYZ7YFb+aEor5Im+IEnoxgiCUgQwUWx55xPKfjxAaZbD5YD5yLlUo0ofCSiPW7z2LuYld3La1nWC2Z3FEPAlEyp6/rWbCzDBYv++s4PdKNYj4NDKlVfV4POU4/v1nuVXjkppdgHcP5Ls1MNjjmefLMLSrxShmV/qPffKzy/faljgQQ1NPKKYmI40gCH7IQFERQifOF3dK81QIZW36b+hmCBO1mvTUyAAg2Gvg6I4fcmNLRIboUC5AKCu1r0I0Mu8cyEffti0wppdBtJ7GMdX98Qs8gu1rtAjRYeWk3rTil0hTN9IIwh8gDYqPsdUzFF8RlwWWD3YCjwzRuWznCrFaDbGVmR2NDDGFEx01EwFaDVZN6i3omlJFkEI1Ms/89wQO55WI1tPYRvTsPHEJ7xzId9le30yLUfEGUdcgCILwJ8iD4kPkzGFiC7ufDsDtNgEfYrUaYuv4cO358+oDwvWYOqg9OkSHeuyOlyqCFOrdumpswI+/C9fEsD0ZGNcCgMVgfX57ttv3FVYaPdLSEARBqB0yUHyAycxg/d6zWJv+m9MxT42TxeN74KGhHa2T8NsPDMDCL08K3v6wRaxWg8/AiI0IwuLx8WgRGuh2z59LHzAwrgUyz5fx9ofdfnGFVgOsnyo9UkXMttClshrR52fH4mh+KUqrhH1WlFDMf5C7ejJBNAXIQPEyqdkFWPb1KRRWut7OkepJiW6ut3vwsRP++r1n8cGhfJTXCDdUpGg15BYgHs0vQdJnx+3Gy7HyrJDtFzMDtBBZp8iWQR2jEBoYgKo6k9u2rSODBelpDOF6LL2rO+ryM62viTE6KKGYf6BU9WSCaOyQgeJFxGRKNTMWb8il8hpsOnRO8DWKrxixPeuinWEQoNVgbmIXzBnRGUfzS1FYWYsV35zizb7qaVZRKQJEdoWZllOI/2VdcplJlc0Jw+pQ5Eph7mqVG6DV4F+3dsTre9xHC93SKRq920a43O5KSuyKOSM6w2xqwE4buYlQo6NlaCAlFPMD+H7zjt9hgiCcIQPFS0jJlBoVpsf7B12LJW3RAHa5SBxXabaGQ7BOi1nX9CneTDfOhVCvEotjThg5UpgLWeX+38iueO+HfFS78KKEBAYAGmBUvIF3u8v2nGaHU7FiY3ceoeSJvXyyReCtrYrGsCUiJq+Rv90bQXgDMlBkwt0DVUqm1NKrRlHvEbNKU0vCqtTsAklCXtsstEISm0WG6GC+llLfcTIQusoN0Grw2j/7uuxvdZ0J979/xGqIHFwwQtBEezS/FMXVDWjVPAiLx/fA7JTjvPfy79s64s4+3l91e2urorFsiYipnkxiZ4JwhgwUGRDyQBWjLWC3WKI80EwA7ldpvk5YZTIzWPDfEx6d4/KVWkERROXV9bh/0xGnz0XsKndsr1i8/cAAtx4foS58thbPw1uOwWiyjHtsRBAeva0jvv6lwG3GWm/hra2KxrQlQtWTCcIzZM+DsmzZMmg0Grv/unfvbj1eW1uL2bNno2XLlggLC8PkyZNRVORcMM1fEFo7R4ygkQEwrpdBcEVbd+diV2lcyFXgUApv7jmDipoGj87BjitfnRlHxNY04hq/sb1icWjhSHz6yGBEBnPnmnEshshFanYBZy2egopavHsgH4vHx9vVqjn23CifGCfujDhAvhpH3riOt6DqyQThGYokauvZsycKCgqs/x08eNB6LCkpCd988w2++OILfP/997h06RImTZqkRDcUR8wDVWgiM9Y+2HzoHJJ35LpsrwF4J0hH1LZKS80uwLo9Zzw6h2PSNbYA3af/cm04MLj+uUhd5QZoNdBqNS6jolwZh+40SQyA5B05MKtgMpZixKn5Ot6CqicThGcoYqA0a9YMBoPB+l90dDQAoKKiAps2bcJrr72GESNGYODAgfjggw/w448/4vDhw0p0RVHEPFCFZkp1nI/4pif2HDOGdhDUVyVWae6q+rp6n7u8JULgEvIGaDXQalwbDoDlc/nwUL5dBldXcI2fJy58IZqkgopa3L/pCOZuy8LU9w5j2Oq9gqpNy423tioa25YIVU8mCM9QRINy5swZtG7dGkFBQUhISMDKlSvRvn17ZGZmor6+HomJida23bt3R/v27ZGRkYEhQ4Zwns9oNMJovL7fX1lZCQCor69Hff31iYj9t+1rSnK5ogr6APeT8uWKKtTXh2Nkt2i8dV9fLPv6lKh8JFwYwoOwcFx33NGtFb7MvICiSm6BqAZATHgQ+rdtLnlcuMY1PbcIq3b9isJK24yvlj4l9nBdCfpofilKr9ZAHyCpO9BqgFf+3g8ju0Vz3lNhubDPZU1qDmKa6xET1gwVNfWixy86pJmg60SHNHN6b2HZVegDGOi1lvez/3dF2dUazNuaibVT+rkdYznx5D59dR1vPwv4YH/zfL8Vvu+wmlHL2DY2msq4irk/DcMwsvqQd+3ahatXr6Jbt24oKCjA8uXLcfHiRWRnZ+Obb77BjBkz7IwNABg0aBDuuOMOrF69mvOcy5Ytw/Lly51eT0lJQUhIiJzdJwiCIAhCIaqrq3HfffehoqIC4eHhLtvKbqA4Ul5ejri4OLz22msIDg6WZKBweVDatWuH4uJiuxusr69HWloaRo0aBZ1OeqE8oZjMDMasO+DWe7F73m12qcwf3nLMo+tunn6z0751em4RVu7MRZFNwcGY5nosurOHx6tt23HVBjTDmHUH7FaDtnDdsyNSx0Coh+bbXy5h4VcnBZ9XAyAiWAd9swAUXRHnEUrPLbIKXbnyyfB5O749UYCFX56AXssg+SYzFv+khdEs3NXP9R1QEqn36avrePtZ0JSgsVWGpjKulZWViI6OFmSgKB5mHBkZia5du+Ls2bMYNWoU6urqUF5ejsjISGuboqIiGAz8lVn1ej30emetgE6n4/wg+V6XGx2AReN7ukx4tmh8TwTpr4cLF1c3WMNJpVJc3eB0fxptAOrMWrtz15m10GgDZBsLnU6Hny5U4nyZEa6UNOfLjDj+5xXe3A5DOrdCVFiwy7wlWg3wxpR+aNk8SHQItKFFmOgxrr3agE8fGQitViPqeuP6tIVGGyAob4dtrpxSh++B0awR1Weu74CSiLlPNV3HW8+CpgiNrTI09nEVc2+KGyhXr15FXl4epk2bhoEDB0Kn02HPnj2YPHkyAOD06dO4cOECEhISlO6KIohNeCaHWNXxHHy5I4oq5c8dIYeQUUjekvVTB1iTkbET+7cnLgkyHIRmY3WkuMqIif3aiHoPICyfDF/laqn4IjTVW3lzfJ2fhyAIdSC7gfLUU0/h7rvvRlxcHC5duoSlS5ciICAAU6dORUREBB555BHMnz8fUVFRCA8PxxNPPIGEhARegaw/IOaBKiTrKR9cNXI8TactNqW42NwOfOd3VfnY1rCzpMHPcRIYLpvAv5q2NYDEjLEnk76r+kN8BqSUCGJP6yR5ipQ6S2q+DkEQ6kV2A+XPP//E1KlTUVJSghtuuAHDhg3D4cOHccMNNwAA1q5dC61Wi8mTJ8NoNGLMmDF466235O6G1xH6QBXiPeCCLyzRk3TaUlKKuzOwbCdQd+d3Z9jxpcEvrKzFY5/8jLddeIb4DCAulJz0pdRgahGiQ1l1vdP3g0JTCYJoSshuoGzbts3l8aCgIGzYsAEbNmyQ+9J+gyvvweLx8Thz+Qo+OHTOLhSZb8tI6paL1JTirgws2wk0LaeQ8/wFFRbj4pGhHZAYb8CgjlGchp3JzGDhl66Frou+POmy0JqtAZSWU4jNh855fdIXWoNpwZhuaBkeYjXS0nIKfV4niSAIwpdQLR4f4dp7EIs5I7oI2nqRkk7bZGaw7Gvp20LudDej4g0YtnqvS6/BpkPnsOnQOV6PzeHfS1Be7Tpevqy6Hod/L8HQztG8bVjPVkKnlhjUMUr2Sd/dFplQA7JlmN5O/0I6DIIgmjpkoPgQV9tCQreMxGy5sKzfe4Y3TBgQVmXV1QSakVciWKDK57HJyCsR9P6MPNcGitA+S0HIFplQA5Iroy3pMAiCaMookuqe8B4BWg0Wj4/nNU4A++2L1OwCrE0XVgPH3eqfr9CgmFTk/EXghKo2xClN5SqOKLRIpJB6LAAwMK6FpH4QBEE0VshA8XNSswuQvIO7ro0hIsjOMyG2Bo7UqBax7+MqApdwozCviNB2ciKmSKSQeiwAaOuGIAjCATJQ/Bi+VTzL4vE97LZNhAo2Ac+qrAqt3OyIredlSKeWiAxxndAnMkSHIT7YAhFbdZfV7Bgi7A03Q0QQ1k7pp2BPCYIg/BfSoPgp7sJXNQCSd+RiTK9YSVsvnkS1SA2ltvW8BGg1WDWpN2eYMcuqSb194nmQEjnFp38xmxqwM1+pnhIEQfgv5EHxU8Su4gHhWy9JiV09DmXl8xpwoQG3x2Zsr1i8/cAAGMLtBaSGcL3LHChKIyVyCpBP/0IQBNEUIA+KnyJlFS8ki60hXI85IzrL0EN58pCoMdxWSuQUQRAEIQ7yoPgpUlbx7gSbGgDLJvSUdfJnvQZL7u5p8YZw6DDc1QpSm+fBdhz5oGyvBEEQnkEeFD9F6ipebHFDOVGjN0QqY3vF4tHbOuK9H/LtaupoNcDMWztStleCIAgPIQPFTxGacl5t2yaNJflYanYB3j2Q72QcMgzw7oF89G/fgowUgiAID6AtHj/GVfiqv22b+BNi8qAQBEEQ0iAPip/TmLZN/AVPKkgTBEEQwiADpRHQWLZN/AWpFaQJgiAI4dAWD0GIhKuwHxdSSwUQBEEQ5EEhCFGkZhdg2denXLahPCgEQRCeQwYKQQiErX3kSvrqLoKKIAiCEAYZKAQhAHe1j1i8kU+GIAiiKUAGCkEIQGgl6Ff+3hdDu0R7oUcEQRCNGxLJEoQAhEbkFFcZFe4JQRBE04AMFIIQgNQKxgRBEIQ0yEAhCAGwtY/4ZK8aALEUuUMQBCEbZKAQhADcVYIGKHKHIAhCTshAIQiBeFL7iCAIghAHRfEQhAio9hFBEIR3IAOFIERCtY8IgiCUh7Z4CIIgCIJQHWSgEARBEAShOshAIQiCIAhCdZCBQhAEQRCE6iADhSAIgiAI1UEGCkEQBEEQqoMMFIIgCIIgVAcZKARBEARBqA4yUAiCIAiCUB1+mUmWYRgAQGVlpd3r9fX1qK6uRmVlJXQ6nS+61iihcVUOGltloHFVDhpbZWgq48rO2+w87gq/NFCuXLkCAGjXrp2Pe0IQBEEQhFiuXLmCiIgIl200jBAzRmWYzWZcunQJzZs3h0ZzvUhbZWUl2rVrhz/++APh4eE+7GHjgsZVOWhslYHGVTlobJWhqYwrwzC4cuUKWrduDa3WtcrELz0oWq0Wbdu25T0eHh7eqD9gX0Hjqhw0tspA46ocNLbK0BTG1Z3nhIVEsgRBEARBqA4yUAiCIAiCUB2NykDR6/VYunQp9Hq9r7vSqKBxVQ4aW2WgcVUOGltloHF1xi9FsgRBEARBNG4alQeFIAiCIIjGARkoBEEQBEGoDjJQCIIgCIJQHWSgEARBEAShOvzKQCktLcX999+P8PBwREZG4pFHHsHVq1ddvufdd9/F8OHDER4eDo1Gg/LyclnO29iQMga1tbWYPXs2WrZsibCwMEyePBlFRUV2bTQajdN/27ZtU/JWfMqGDRvQoUMHBAUFYfDgwTh69KjL9l988QW6d++OoKAg9O7dGzt37rQ7zjAMlixZgtjYWAQHByMxMRFnzpxR8hZUi9xj+9BDDzl9N8eOHavkLagSMeN66tQpTJ48GR06dIBGo8G6des8PmdjRu6xXbZsmdN3tnv37gregY9h/IixY8cyffv2ZQ4fPsz88MMPTOfOnZmpU6e6fM/atWuZlStXMitXrmQAMGVlZbKct7EhZQwee+wxpl27dsyePXuYn376iRkyZAhzyy232LUBwHzwwQdMQUGB9b+amholb8VnbNu2jQkMDGQ2b97MnDp1ipk5cyYTGRnJFBUVcbY/dOgQExAQwKxZs4bJyclhnn/+eUan0zEnT560tlm1ahUTERHB/O9//2N++eUXZsKECUzHjh0b7RjyocTYTp8+nRk7dqzdd7O0tNRbt6QKxI7r0aNHmaeeeorZunUrYzAYmLVr13p8zsaKEmO7dOlSpmfPnnbf2b/++kvhO/EdfmOg5OTkMACYY8eOWV/btWsXo9FomIsXL7p9/759+zgNFE/P2xiQMgbl5eWMTqdjvvjiC+trubm5DAAmIyPD+hoA5quvvlKs72pi0KBBzOzZs61/m0wmpnXr1szKlSs52//zn/9kxo8fb/fa4MGDmX//+98MwzCM2WxmDAYD8/LLL1uPl5eXM3q9ntm6dasCd6Be5B5bhrEYKBMnTlSkv/6C2HG1JS4ujnMS9eScjQklxnbp0qVM3759ZeyluvGbLZ6MjAxERkbipptusr6WmJgIrVaLI0eOqO68/oSUMcjMzER9fT0SExOtr3Xv3h3t27dHRkaGXdvZs2cjOjoagwYNwubNmwWV2fY36urqkJmZaTceWq0WiYmJTuPBkpGRYdceAMaMGWNtn5+fj8LCQrs2ERERGDx4MO85GyNKjC3L/v370apVK3Tr1g2zZs1CSUmJ/DegUqSMqy/O6Y8oOQ5nzpxB69atceONN+L+++/HhQsXPO2uavEbA6WwsBCtWrWye61Zs2aIiopCYWGh6s7rT0gZg8LCQgQGBiIyMtLu9ZiYGLv3rFixAp9//jnS0tIwefJkPP7443jzzTdlvwdfU1xcDJPJhJiYGLvXHcfDlsLCQpft2f+LOWdjRImxBYCxY8fio48+wp49e7B69Wp8//33GDduHEwmk/w3oUKkjKsvzumPKDUOgwcPxocffojU1FRs3LgR+fn5uPXWW3HlyhVPu6xKfF7NeOHChVi9erXLNrm5uV7qTeNCDWO7ePFi67/79++PqqoqvPzyy/i///s/Ra9LEO649957rf/u3bs3+vTpg06dOmH//v0YOXKkD3tGENyMGzfO+u8+ffpg8ODBiIuLw+eff45HHnnEhz1TBp8bKE8++SQeeughl21uvPFGGAwGXL582e71hoYGlJaWwmAwSL6+UudVA0qOrcFgQF1dHcrLy+28KEVFRS7HbfDgwUhOTobRaGxUNSeio6MREBDgFMXkajwMBoPL9uz/i4qKEBsba9emX79+MvZe3SgxtlzceOONiI6OxtmzZ5uEgSJlXH1xTn/EW+MQGRmJrl274uzZs7KdU034fIvnhhtuQPfu3V3+FxgYiISEBJSXlyMzM9P63r1798JsNmPw4MGSr6/UedWAkmM7cOBA6HQ67Nmzx/ra6dOnceHCBSQkJPD2KSsrCy1atGhUxgkABAYGYuDAgXbjYTabsWfPHt7xSEhIsGsPAGlpadb2HTt2hMFgsGtTWVmJI0eOuBzjxoYSY8vFn3/+iZKSEjtjsDEjZVx9cU5/xFvjcPXqVeTl5TXe76yvVbpiGDt2LNO/f3/myJEjzMGDB5kuXbrYhcL++eefTLdu3ZgjR45YXysoKGCOHz/OvPfeewwA5sCBA8zx48eZkpISwedtCkgZ28cee4xp3749s3fvXuann35iEhISmISEBOvxr7/+mnnvvfeYkydPMmfOnGHeeustJiQkhFmyZIlX781bbNu2jdHr9cyHH37I5OTkMI8++igTGRnJFBYWMgzDMNOmTWMWLlxobX/o0CGmWbNmzCuvvMLk5uYyS5cu5QwzjoyMZLZv386cOHGCmThxYpMNM5ZzbK9cucI89dRTTEZGBpOfn8+kp6czAwYMYLp06cLU1tb65B59gdhxNRqNzPHjx5njx48zsbGxzFNPPcUcP36cOXPmjOBzNhWUGNsnn3yS2b9/P5Ofn88cOnSISUxMZKKjo5nLly97/f68gV8ZKCUlJczUqVOZsLAwJjw8nJkxYwZz5coV6/H8/HwGALNv3z7ra0uXLmUAOP33wQcfCD5vU0DK2NbU1DCPP/4406JFCyYkJIS55557mIKCAuvxXbt2Mf369WPCwsKY0NBQpm/fvszbb7/NmEwmb96aV3nzzTeZ9u3bM4GBgcygQYOYw4cPW4/dfvvtzPTp0+3af/7550zXrl2ZwMBApmfPnsyOHTvsjpvNZmbx4sVMTEwMo9frmZEjRzKnT5/2xq2oDjnHtrq6mhk9ejRzww03MDqdjomLi2NmzpzZ5CZRhhE3ruxzwPG/22+/XfA5mxJyj+2UKVOY2NhYJjAwkGnTpg0zZcoU5uzZs168I++iYZhGGPNJEARBEIRf43MNCkEQBEEQhCNkoBAEQRAEoTrIQCEIgiAIQnWQgUIQBEEQhOogA4UgCIIgCNVBBgpBEARBEKqDDBSCIAiCIFQHGSgEQRAEQagOMlAIgiAIglAdZKAQBEEQBKE6yEAhCIIgCEJ1kIFCEARBEITq+H/WF8rsVOz0ZwAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Find the index of the single variable with the best squared loss\n",
        "imax = np.argmin(losses)\n",
        "\n",
        "# Regression line over the range of x values\n",
        "xmin = np.min(X[:,imax])\n",
        "xmax = np.max(X[:,imax])\n",
        "ymin = beta0[imax] + beta1[imax]*xmin\n",
        "ymax = beta0[imax] + beta1[imax]*xmax\n",
        "plt.plot([xmin,xmax], [ymin,ymax], 'r-', linewidth=3)\n",
        "\n",
        "# Scatter plot of points\n",
        "plt.scatter(X[:,imax],y)\n",
        "plt.grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vO8AMlHqtiHr"
      },
      "source": [
        "Note that the calculations above could have been done without a for-loop using Python broadcasting:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "id": "VzuBzSLGtiHr"
      },
      "outputs": [],
      "source": [
        "# Remove the means\n",
        "ym = np.mean(y)\n",
        "y1 = y-ym  # a column vecotor each minus mean\n",
        "Xm = np.mean(X,axis=0) # averaging over column, resulting a row vector of dimension natt\n",
        "X1 = X - Xm[None,:] # minus the same mean in each column\n",
        "\n",
        "# Compute the correlations per features\n",
        "Sxx = np.mean(X1**2,axis=0) #a row vector with each element indicating the variance of one attribute\n",
        "Sxy = np.mean(X1*y1[:,None],axis=0) #a row vector with each element indicating the covarance on one attribute to the targer\n",
        "\n",
        "# Compute the coefficients and losses per feature\n",
        "beta1 = Sxy/Sxx # element wise division, resulting a row vector containing  beta1 for each attribute\n",
        "beta0 = ym - beta1*Xm # element wise multiplication, resulting a row vector containing beta0 for each attribute\n",
        "errs = ((X*beta1) + beta0) - y[:,None] # results in a matrix where every column contains the residuals for predictor k\n",
        "losses = np.sum(errs**2,axis=0) # sums up the squared size of the residuals across each column"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o4fnZ8FttiHr"
      },
      "source": [
        "Now let us see whether we get the same result:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qYVWsF_LtiHr",
        "outputId": "35bd8cc6-db7a-440a-8b3b-ff2b94b3c87c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0 loss=2.53e+06 beta0=152.13348416289594 beta1=304.1830745283064\n",
            "1 loss=2.62e+06 beta0=152.13348416289594 beta1=69.71535567841418\n",
            "2 loss=1.72e+06 beta0=152.13348416289617 beta1=949.4352603840389\n",
            "3 loss=2.11e+06 beta0=152.13348416289597 beta1=714.7382594960393\n",
            "4 loss=2.50e+06 beta0=152.13348416289594 beta1=343.2544518889671\n",
            "5 loss=2.54e+06 beta0=152.13348416289594 beta1=281.7845933524569\n",
            "6 loss=2.21e+06 beta0=152.13348416289594 beta1=-639.1452793225359\n",
            "7 loss=2.14e+06 beta0=152.13348416289594 beta1=696.8830300922268\n",
            "8 loss=1.78e+06 beta0=152.13348416289585 beta1=916.137374550914\n",
            "9 loss=2.24e+06 beta0=152.13348416289594 beta1=619.2228206843718\n"
          ]
        }
      ],
      "source": [
        "for k in range(natt):\n",
        "    print(str(k)+\" loss=\"+\"{:.2e}\".format(losses[k])+\" beta0=\"+str(beta0[k])+\" beta1=\"+str(beta1[k]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RQs1UCIQtiHs"
      },
      "source": [
        "## Improvements with a Multiple Variable Linear Model\n",
        "\n",
        "One possible way to try to improve the fit is to use multiple variables at the same time.\n",
        "\n",
        "We can manually compute the regression coefficients using the least-squares matrix formula given in class.\n",
        "\n",
        "To compute the coefficients manually, we first append an all ones columns onto `X` using the `ones` command and `hstack`.  Note that after we do this, `X` will have 11 columns -- one more column that original data matrix."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "62hNu-LntiHs",
        "outputId": "148659a1-1b47-4029-bddd-1d7c1c866dd0"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(442, 11)"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ],
      "source": [
        "ones = np.ones((nsamp,1))\n",
        "X_orig = X\n",
        "X = np.hstack((ones,X_orig))\n",
        "X.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1ZV9nIC3tiHs"
      },
      "source": [
        "Next we need to solve $\\min_{\\mathbf{\\beta}}\\|\\mathbf{y} - \\mathbf{X\\beta}\\|_2^2$. The most obvious way to do so is via the direct matrix equation derived in class: $\\mathbf{\\beta}^* = (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{y}$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SwdWIJLPtiHs",
        "outputId": "4434186f-8f54-465c-825c-80735a18def9"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 152.13348416  -10.0098663  -239.81564367  519.84592005  324.3846455\n",
            " -792.17563855  476.73902101  101.04326794  177.06323767  751.27369956\n",
            "   67.62669218]\n"
          ]
        }
      ],
      "source": [
        "Xt = np.transpose(X)\n",
        "beta=np.dot(np.linalg.inv(np.dot(Xt,X)),np.dot(Xt,y))\n",
        "print(beta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZELv7UEktiHs"
      },
      "source": [
        "Note that `beta[0]` is the value of the intercept from the regression fit."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o2Fbul9wtiHs"
      },
      "source": [
        "The notation looks a bit cleaner if we use `numpy`'s `matrix` type. Note that the `matrix` type is a subclass of `ndarray`. It really has no advantages other than the fact that the multiplication operator `*` performs actual matrix multiplication. Seems like a small thing, but can save you a lot of headache over using multiple nested `np.dot` operations."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "scrolled": false,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "K8gjSk3otiHs",
        "outputId": "ecc147ee-c043-4d18-fd1e-dd013bb04251"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 152.13348416]\n",
            " [ -10.0098663 ]\n",
            " [-239.81564367]\n",
            " [ 519.84592005]\n",
            " [ 324.3846455 ]\n",
            " [-792.17563855]\n",
            " [ 476.73902101]\n",
            " [ 101.04326794]\n",
            " [ 177.06323767]\n",
            " [ 751.27369956]\n",
            " [  67.62669218]]\n"
          ]
        }
      ],
      "source": [
        "Xmat=np.matrix(X)\n",
        "ymat=np.matrix(y)\n",
        "ymat=np.transpose(ymat)\n",
        "Xmatt=np.transpose(Xmat)\n",
        "beta=np.linalg.inv(Xmatt*Xmat)*Xmatt*ymat\n",
        "print(beta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5UaFFjNctiHs"
      },
      "source": [
        "An often slightly faster alternative is to use the `lstsq` method, which solves $\\min_{\\mathbf{\\beta}}\\|\\mathbf{y} - \\mathbf{X\\beta}\\|_2^2$ using a QR decomposition. This should find the exact same least-squares fit. You might get a warning when you run `lstqr` about some defaults changing in the future. You can safetly ignore."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_lOWXS4qtiHs",
        "outputId": "8456933c-e2aa-44f1-885c-4a0c907b1b9b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 152.13348416  -10.0098663  -239.81564367  519.84592005  324.3846455\n",
            " -792.17563855  476.73902101  101.04326794  177.06323767  751.27369956\n",
            "   67.62669218]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-16-b80b493bcb0c>:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
            "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
            "  out = np.linalg.lstsq(X,y)\n"
          ]
        }
      ],
      "source": [
        "out = np.linalg.lstsq(X,y)\n",
        "beta = out[0]\n",
        "print(beta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yy-znQlTtiHt"
      },
      "source": [
        "Finally, for very large problems (many examples and/or features), *iterative* methods can offer a much faster approach to solving $\\min_{\\mathbf{\\beta}}\\|\\mathbf{y} - \\mathbf{X\\beta}\\|_2^2$. There are several options avialable in the `scipy` library, including the populat LSQR method."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "id": "gPJTKWKstiHt"
      },
      "outputs": [],
      "source": [
        "import scipy as sp"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xpMDWfS3tiHt",
        "outputId": "5e6c79e8-b109-409e-fb4f-9fb17be236de"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 152.13348416  -10.00985784 -239.81563637  519.8459319   324.38465612\n",
            " -792.17562514  476.73903475  101.0432569   177.06325442  751.27371435\n",
            "   67.62670477]\n"
          ]
        }
      ],
      "source": [
        "out = sp.sparse.linalg.lsqr(X,y)\n",
        "beta = out[0]\n",
        "print(beta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-1NT8ayZtiHt"
      },
      "source": [
        "Hopefully you got the same values of $\\mathbf{\\beta}$ no matter which algorithm you used. The differences are negligble because this demo problem is so small, but if you ever want to compare runtimes head-to-head, you can use the `time` library to time different operations. Here's an example that times the direct approach:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0HTK-ESttiHu",
        "outputId": "4e1c5246-34bc-409b-a00f-f156a5b82dc1"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "direct matrix computation time: 0.0011584758758544922 seconds\n"
          ]
        }
      ],
      "source": [
        "import time\n",
        "t = time.time()\n",
        "beta=np.linalg.inv(Xmatt*Xmat)*Xmatt*ymat\n",
        "elapsed = time.time() - t\n",
        "print(\"direct matrix computation time: \" + str(elapsed) + \" seconds\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "df3TB8QVtiHu"
      },
      "source": [
        "Now let's compute the loss for our muliple variable model. Be careful with arrays vs. matrices! Always check what you're working with by running commands like `type(beta)`. Or multiplications may not behave as expected."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "a_943gB5tiHu",
        "outputId": "f1fd5782-78cf-46f3-c86d-c95daa3ad8ee"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "multiple variable loss=1.26e+06\n"
          ]
        }
      ],
      "source": [
        "yp = Xmat*beta\n",
        "errs = np.array(ymat - yp)\n",
        "lossm = np.sum(errs**2)\n",
        "\n",
        "print(\"multiple variable loss=\"+\"{:.2e}\".format(lossm))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0BBVANNAtiHu"
      },
      "source": [
        "We could have also done everything using matrix operations!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GFy_egaktiHu",
        "outputId": "b0508408-21bd-4531-9e78-b44fe1893249"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "multiple variable loss=1.26e+06\n"
          ]
        }
      ],
      "source": [
        "yp = Xmat*beta\n",
        "lossm = np.linalg.norm(np.array(ymat - yp))**2\n",
        "\n",
        "print(\"multiple variable loss=\"+\"{:.2e}\".format(lossm))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5DKVzCwCtiHu"
      },
      "source": [
        "Note that this loss is about $30\\%$ lower than what we acheived with the best single variable linear regression, which is significant. Let's also normalize the loss to get an $R^2$ value, which we can compare to the value of $34\\%$ we obtained earlier. Higher $R^2$ is better."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uBs0nZqytiHu",
        "outputId": "dbf06d1e-2d4f-45f1-d3ca-b20fa8a39193"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.5177484222203501"
            ]
          },
          "metadata": {},
          "execution_count": 22
        }
      ],
      "source": [
        "rsqr = 1 - lossm/(nsamp*syy)\n",
        "rsqr"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rzLalDnFtiHu"
      },
      "source": [
        "## Using built in sci-kit learn model\n",
        "\n",
        "We can also fit the multiple variable linear regression model using Sci-kit Learn's built in functionality. As mentioned, this is a bit overkill at this point, but it's good to see because you might use other built in models from `sklearn` later in the course. The disadvantage is that this approach leaves everything inside a \"black-box\". For example, how do you think this method is solving the required least-squares optimization problem? Using an iterative method like `lsqr`? Or a direct method?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5QlqYnudtiHu"
      },
      "source": [
        "To fit the linear model, we first create a regression object and then fit the data with regression object."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 75
        },
        "id": "5HJFQTVwtiHv",
        "outputId": "0cfb9bc9-f848-470a-dfa9-bbd5ba3eb4e6"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "LinearRegression()"
            ],
            "text/html": [
              "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>LinearRegression()</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LinearRegression</label><div class=\"sk-toggleable__content\"><pre>LinearRegression()</pre></div></div></div></div></div>"
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ],
      "source": [
        "X = X[:,1:11]\n",
        "regr = linear_model.LinearRegression()\n",
        "regr.fit(X,y)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QNES445OtiHv"
      },
      "source": [
        "You can see the intercept and other coefficients."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mPaxy52-tiHv",
        "outputId": "b2b234cf-8854-4586-ef68-000ed15b0a31"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "152.13348416289597"
            ]
          },
          "metadata": {},
          "execution_count": 24
        }
      ],
      "source": [
        "regr.intercept_"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "scrolled": true,
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "o7gl1BIftiHv",
        "outputId": "2497a253-b6d2-433e-b786-e91b9244833d"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([ -10.0098663 , -239.81564367,  519.84592005,  324.3846455 ,\n",
              "       -792.17563855,  476.73902101,  101.04326794,  177.06323767,\n",
              "        751.27369956,   67.62669218])"
            ]
          },
          "metadata": {},
          "execution_count": 25
        }
      ],
      "source": [
        "regr.coef_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ny5H16fGtiHv"
      },
      "source": [
        "We next compute the squared loss by using the `predict` function and you can verify that it's the same as above."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "y6XWltgCtiHv",
        "outputId": "2a21330f-fc0a-44f3-9cc2-7683fc72641c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "multiple variable loss=1.26e+06\n"
          ]
        }
      ],
      "source": [
        "y_pred = regr.predict(X)\n",
        "lossm = np.linalg.norm(y_pred - y)**2\n",
        "\n",
        "print(\"multiple variable loss=\"+\"{:.2e}\".format(lossm))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "id": "cCAn2MaOtiHv"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.4"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}