Logistic_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression\n",
    "\n",
    "### Joseph Melby, 2023\n",
    "\n",
    "This notebook contains an implementation and testing of a logistic classifier written (mostly) from\n",
    "scratch. This particular implementation only uses NumPy, pandas, and some visualization tools, as\n",
    "well as some data manipulation tools from scikit-learn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression Class\n",
    "\n",
    "Here we will define a logistic regression class with all the necessary functions. We will need a\n",
    "couple functions outside of this class as well.\n",
    "\n",
    "We use gradient descent as an optimization algorithm to learn the best fitting parameters for our\n",
    "model. We do this using the first derivatives of our cost function to find the direction of steepest\n",
    "descent on the parameter space. \n",
    "\n",
    "fit() is a vectorized version of the gradient descent algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standardize data function\n",
    "def standardize(X):\n",
    "    for i in range(X.shape[1]):\n",
    "        X[:, i] = (X[:, i] - np.mean(X[:, i]))/np.std(X[:, i])\n",
    "\n",
    "# Define an accuracy score\n",
    "# We will use the F1 score, which is the harmonic mean of precision and recall. \n",
    "def F1_score(y,y_hat):\n",
    "    tp,tn,fp,fn = 0,0,0,0\n",
    "    for i in range(len(y)):\n",
    "        if y[i] == 1 and y_hat[i] == 1:\n",
    "            tp += 1\n",
    "        elif y[i] == 1 and y_hat[i] == 0:\n",
    "            fn += 1\n",
    "        elif y[i] == 0 and y_hat[i] == 1:\n",
    "            fp += 1\n",
    "        elif y[i] == 0 and y_hat[i] == 0:\n",
    "            tn += 1\n",
    "    precision = tp/(tp+fp)\n",
    "    recall = tp/(tp+fn)\n",
    "    f1_score = 2*precision*recall/(precision+recall)\n",
    "    return f1_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression():\n",
    "\n",
    "    # Define sigmoid function\n",
    "    def sigmoid(self, z):\n",
    "        sig = 1/(1+np.e**(-z))\n",
    "        return sig\n",
    "    \n",
    "    # Initialization\n",
    "    def initialize(self, X):\n",
    "        # initialize weights\n",
    "        weights = np.zeros((X.shape[1] + 1, 1))\n",
    "        # Add a bias column of ones to X\n",
    "        X = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)\n",
    "\n",
    "        return weights, X\n",
    "    \n",
    "    # Define cost function for logistic regression\n",
    "    def cost(self, theta, X, y):\n",
    "        z = np.dot(X, theta)\n",
    "        cost0 = y.T.dot(np.log(self.sigmoid(z)))\n",
    "        cost1 = (1 - y).T.dot(np.log(1-self.sigmoid(z)))\n",
    "        cost = -(cost1 + cost0)/len(y)\n",
    "        return cost\n",
    "    \n",
    "    def fit(self, X, y, alpha = 0.01, num_it = 100):\n",
    "        params, X = self.initialize(X)\n",
    "        costs = np.zeros(num_it,)\n",
    "        for i in range(num_it):\n",
    "            params = params - alpha*np.dot(X.T, self.sigmoid(np.dot(X, params)) - np.reshape(y, (len(y), 1)))\n",
    "            costs[i] = self.cost(params, X, y)\n",
    "        self.weights = params\n",
    "        return costs\n",
    "    \n",
    "    # Define a prediction function\n",
    "    def predict(self, X):\n",
    "        z = np.dot(self.initialize(X)[1], self.weights)\n",
    "        pred = []\n",
    "        for i in self.sigmoid(z):\n",
    "            if i > 0.5:\n",
    "                pred.append(1)\n",
    "            else: \n",
    "                pred.append(0)\n",
    "        return pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic Synthetic Data Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate the random data\n",
    "np.random.seed(440)\n",
    "n_rows = 100\n",
    "diff = .1\n",
    "X = np.random.random((n_rows,2))\n",
    "X_prime = X[(X[:,1] - X[:,0]) <= -diff,:]\n",
    "X_2prime = X[(X[:,1] - X[:,0]) >= diff,:]\n",
    "\n",
    "del X\n",
    "X = np.append(X_prime,X_2prime,axis = 0)\n",
    "\n",
    "y = np.empty(np.shape(X)[0])\n",
    "y[(X[:,1] - X[:,0]) <= -diff] = 1\n",
    "y[((X[:,1] - X[:,0]) >= diff)] = 0\n",
    "\n",
    "X[1,:] = [.8,.1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJbCAYAAAChCjPHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNW0lEQVR4nO3de3RU9b3//9dkMtzETPCGoICRmICSaA2VmOhq7bFYbLU0bcV6VgNSz+/QoigcrSIexUuL9vZVOYC2iHBW0XKqgdPTRT1mnbbKRWmhgGnDKg1ewApFMDe0hTDZvz+2ExIy12Rm78/s/XysNStmZid8JhOTV96fz+f9CViWZQkAAABGyHN7AAAAADiBcAYAAGAQwhkAAIBBCGcAAAAGIZwBAAAYhHAGAABgEMIZAACAQfLdHoBbOjs79d577+nUU09VIBBwezgAAMDjLMtSe3u7Ro4cqby8+PUx34az9957T6NGjXJ7GAAAwGf27dunc889N+7jvg1np556qiT7C1RQUODyaAAAgNe1tbVp1KhRXRkkHt+Gs+hUZkFBAeEMAAA4JtlyKjYEAAAAGIRwBgAAYBDCGQAAgEEIZwAAAAbx7YYAAADQP5FIRB0dHW4Pw3WhUEjBYDBjn49wBgAA0mJZlg4cOKDW1lZZluX2cFwXCAQUDod19tlnZ6SxPeEMAACkpbW1VS0tLTrzzDN1yimn+PqkHcuy9OGHH+r999/X4MGDVVhY2O/PSTgDAAApsyxLBw8eVEFBgc444wy3h2OEwYMH6+jRozp48KDC4XC/wyobAgAAQMoikYgikQgN3E9SUFDQ9bXpL8IZAABI2fHjxyVJ+flMvnUX/XpEvz79QTgDAABp8/M6s1gy+fUgnAEAABiEcAYAAGAQwhkAAMDHjhw5ojvuuEMjR47UoEGDdMkll+hnP/uZo2MwIpy9+uqruu666zRy5EgFAgGtW7cu6ce88sorqqio0KBBg3T++efrqaeeyv5AAQBAVjQ0SA8+KM2da79taHBnHDU1NVq1apUeeOAB/epXv9InP/lJfe1rX9Nzzz3n2BiM2Grx4Ycf6uKLL9bNN9+sL3/5y0mvf+utt3TttdfqX/7lX/TTn/5UmzZt0re+9S2deeaZKX08AAAwQ1OTNH26tHmzFAxKeXlSZ6e0cKFUXS2tXCkVFzszlvXr16u+vl7PPfecvva1r0mSrrrqKr3zzju66667NG3atIwe0xSPEeFsypQpmjJlSsrXP/XUUxo9erQef/xxSdL48eO1detW/eAHPyCcwVUNDVJdndTSIhUWSjU1UlmZ26MCADM1NUmTJkmtrfb7kYh9i3r9dfvxLVucCWhr167V0KFD9dWvfrXH/TfffLNuuukmbdmyRVVVVVkfhxHTmul67bXXNHny5B73XXPNNdq6dWvcA1iPHj2qtra2HjcgU5qa7L/wysulhx+Wliyx35aXS1dcYT8OAOhp+nQ7mMXr2xqJ2I/PmOHMeP74xz9q/PjxvXq4lZeXdz3uhJwMZwcOHNDw4cN73Dd8+HAdP35chw4divkxixYtUjgc7rqNGjXKiaHCB6J/+W3ZYr8fiUgdHSd+2ET/8iOgwW2mrOkBJPv7b/Pm+MEsKhKRNm1y5vv18OHDOu2003rdH73v8OHD2R+EDJnW7IuTm71ZlhXz/qj58+dr3rx5Xe+3tbUR0JAR6fzlt3Gjo0MDJJm1psdPWOaQWF2d/f2YymlHwaC0dq0zX79EzWSdarybk+Hs7LPP1oEDB3rcd/DgQeXn5+v000+P+TEDBw7UwIEDnRgefCT6l18y3f/y44cznGTamh4/IAynpqXF/tqkEs7y8qTm5qwPSaeffnrM6tgHH3wgSTGratmQk9Oal19+uerr63vc9/LLL2vixIkKhUIujQp+FP3LLxXRv/wAJ5m2psfrWOaQusJCO7SmorNTGjYsq8ORJJWVlWnXrl29zsds+HhOdcKECdkfhAwJZ0eOHNGOHTu0Y8cOSXarjB07dmjv3r2S7CnJ2trarutnzZqld955R/PmzdOuXbu0YsUKPfPMM7rzzjvdGD58LPqXXyqc+ssPiDJxTY/XEYZTV1OTWtVMsq+rqcnueCTpS1/6ko4cOaIXX3yxx/2rVq3SyJEjNWnSpOwPQoZMa27dulVXXXVV1/vRtWHTp0/XypUrtX///q6gJklFRUVav3695s6dqyVLlmjkyJF68sknaaMBx5n4lx8QZeqaHq9imUN6ysqkqiq7ypjoezQYlCorJSeKVlOmTNFnP/tZffOb31RbW5uKi4v1/PPP66WXXtJPf/pTR3qcSYaEs09/+tNdC/pjWblyZa/7PvWpT+kPf/hDFkcFJFdTY68jSYVTf/kBUSau6fEywnD6Vq06sSYy1tctGJTCYXudnlPq6uq0YMEC3X///frggw80btw4Pf/887rxxhsdG4MR05pAror+5Zfsj6lg0F4I7NByBUASlV2nscwhfcXFduWsstJ+PxiUQqETP1MrK53frDJ06FA98cQT2r9/v44ePaqdO3c6GswkQypnQC4z8S8/QKKy6zTCcN8UF9tthhoa7Gpic7P9tamp8e8ftIQzoJ+if/nNmGGvI+m+dT4Ssf/yY+s83GDimh4vIwz3T1kZ07xRhDMYJVebNvKXH0xFZdc5hGFkCuEMRkjUtHHUKOmf/kk67zzzwxp/+cE0VHadRRhGJhDO4LpkHcz37bN/kOXl0WEb6Asqu84hDCMTCGdwXbKmjVHRhbYcNwP0DZVdZxCG0V+EM7gq1aaN3XGQOIBcQBhGX9HnDK5K52zK7jhuBgDgVYQzuCqdpo0n4yBxAIAXEc7gqnSaNp6MDtsAAC8inMFVNTWpnUMXCx22AQCZ1t7erm9/+9uaPHmyzjzzTAUCAS1MtbtwhrAhAK5KtWljLHTYBoAc1vyGtGe5JCuFiwPS2FukYeXZHpUOHz6sH//4x7r44os1depULV++POv/5skIZ3BdsqaNsdBhGwBy3PsbpN2LpUBQiSfyOiUrIhWMcyScjRkzRs3NzQoEAjp06JAr4YxpTbgu2rSxstJ+P9nuTTpsA4AHFNVKoUI7eFkdCW4R+7qiWkeGFQgEFAgEHPm34iGcwQjRpo1vvCHdf7/dw2zUKPuxYFAKhU6EtspKGtACQM4LnSpdeLekZEEoIF10jxQa6sSojMC0JoxyctNGOmwDgIeVzJYaH5M6WuJfEwpLF8x2bEgmIJzBaHTYBgAPi1bPdt6r2BsD/Fc1k5jWBAAAbiqZbVfHYvFh1UwinAEAADfFXXvmz6qZRDgDAABui1U982nVTGLNGQAAcFuvtWfuVs1+9atf6cMPP1R7e7skqbGxUS+88IIk6dprr9WQIUOy+u8HLMtKpTWv57S1tSkcDqu1tVUFBQVuDwcAgJzwj3/8Q2+99ZaKioo0aNCgzH3ijnZp3Wh752aoUJq6z7Vwdt555+mdd96J+dhbb72l8847r9f9qXxdUs0eTGsCAAD3da09k+trzd5++21ZlhXzFiuYZRrTmgZqaJDq6qSWFqmw0O7tRTsJAIDnlc6R8odKY2e6PRJXEc4M0tQkTZ8ubd5sd8PPy5M6O6WFC6Xqavu4IrriAwA8K3+IVHqr26NwHdOahmhqsg//3rLFfj8SkTo6ThwE/vrr9uNNTe6NEQAAZB/hzBDTp0utrSfC2MkiEfvxGTMcHRYAAHAY4cwADQ32VGa8YBYViUibNtnXAwAAbyKcGaCuzl5jlopg0D4IHAAAN/m0E1dcmfx6EM4M0NJiL/5PRV6e1Nyc1eEAABBXKBSSJH300Ucuj8Qs0a9H9OvTH+zWNEBhob0rMxWdndKwYVkdDgAAcQWDQRUWFurgwYOSpCFDhigQOPlcTP+wLEsfffSRDh48qMLCQgVTnQpLgHBmgJoau11GKiIR+3oAANxy9tlnS1JXQINUWFjY9XXpL8KZAcrKpKoqu41Gok0BwaBUWSlNmODc2AAAOFkgENCIESN01llnqaOjw+3huC4UCmWkYhZFODPEqlV2H7N47TSCQSkcthvRZgOnEgAA0hUMBjMaSmBjQ4Ahiovtylllpf1+MCiFQid2cVZW2o9n+oSApib79IHycunhh6UlS+y35eXSFVfQ9BYAAKcFLJ/uhU31ZHg3NDTY7TKam+3F/zU12ZnKjJ5KkKxal41QCACA36SaPZjWNFBZmTNTiumcSrBxY/bHAwAAmNb0LU4lAADATIQzn+JUAgAAzEQ48ylOJQAAwEyEM5/iVAIAAMxEOPOpmprk682iOJUAAADnEM58KnoqQbJ1Z8Gg3QeNUwkAAHAGrTSyJBc67rt9KgEAAOiNylmG5VLHfbdOJQAAAPFxQkAGTwjI5Y77Tp1KAACAX3FCgAtyueO+U6cSAACAxJjWzBA67gMAgEwgnGUIHfcBAEAmEM4yhI77AAAgEwhnGULHfQAAkAmEswyh4z4AAMgEwlmG0HEfAABkAuEsg1atsvuYxQtodNwHAADJEM4yiI77AACgv2hCm2HFxXaDWTruAwCAviCcZQkd95FtDQ12f72WFnu3cE0N33MA4AWEMyDHNDXZR4Vt3mxPmefl2e1ZFi60N5usXMnUOQDkMtacATmkqUmaNMleuyjZbVk6Ok60cXn9dfvxpib3xggA6B/CGZBDpk+XWlvj99SLROzHZ8xwdFgAgAxiWhPIEQ0N9lRmMpGItGmTfT1r0JBtqax9ZH0kkB7CGZAj6ursNWapnEQRDNq7hfkFiGxJZe2jxPpIoC8IZ0COaGmxf7mlEs7y8uw2LkA2RNc+trba70ciPb8vX39dmjjR/u8jR+JfE10/SUADemLNGZAjCgvtqkMqOjvt/npANqS69pH1kUDfEM6AHFFTk1rVTLKvq6nJ7njgT9G1j6l+LybSfX0kgBMIZ0COKCuTqqrin90aFQza63k4kQLZEF37mCnR9ZEATmDNGZBDVq06sdYnVuUiGJTC4ROLsdE/7DLsLZ21j6lgfSTQG5UzIIcUF9sLqCsr7feDQSkUOlHJqKxkgXUmNDXZ1cfycunhh6UlS+y35eXSFVf4u8lvOmsfU8H6SKC3gGVZltuDcENbW5vC4bBaW1tVUFDg9nB64S92JNPQYE8HNTfbv9xqapjKzITuOxETVSf9GoIbGuyQmunPyfcu/CDV7EE4MyycxesdFInQFwhwQnW1HbwSTdsFg3aVcuNG58ZlklS+Rqnw+9cR/kM4S8LEcMZf7IC70q0KvfFGdiraplfOm5rsPmbRPmexnHqq/cflkSP8PAOiUs0erDkzCOcmAu5KZydiNnYZemmtW16e9OKLrI8E+oLdmobg3ETAfW6ewpBK131TOupPn36i8388R45IDzxgT1myPhJID+HMEJybCLjPzVMY0qmcu7lGq69/SPLzCkgd05qGiP7Fngr6AgHZ4dYpDKl23U+no35Dg/Tgg9LcufbbTHXhd3vqF/ADKmeG4NxEwH3RUxhS3a2Zqam5TFbO4+34XrgwMzu+3Zz6BfyCypkhODcRMMOqVfYuwnjVoWycwpCpynl03dqWLfb7kYjU0XHiZ0t03Vp/NhbwhySQfYQzQ3BuImAGN05hyFTgcWLHN39IAtlHnzP6nAGIw6ldhun2V4vVUd/JHm006gX6JtXswZozg0T/Yp8xw170e/IJAZWVnBAAOMmpXYaZWOvm5I7vVatS+0Myk1O/gJ8QzgxTXExfIMCP+ht4nFyozx+SQHYRzgxFXyDAX/obeJxeqM8fkkD2sObMoDVnACD1LfBkYt0agOxizRkA5Ki+VM7d6tEGIPNopQEAHuFGjzYAmUc4AwCPcKNHG4DMY1oTADyEhfpA7iOcAYAHseMbyF3GTGsuXbpURUVFGjRokCoqKrRhw4aE169evVoXX3yxhgwZohEjRujmm2/W4cOHHRotAABAdhgRztasWaM77rhDCxYs0Pbt23XllVdqypQp2rt3b8zrN27cqNraWn3jG9/Qn/70J/385z/X73//e91yyy0OjxwAACCzjOhzNmnSJF166aVatmxZ133jx4/X1KlTtWjRol7X/+AHP9CyZcu0Z8+ervsWL16s733ve9q3b19K/yZ9zgAAcEHzG9Ke5ZJSiR8Baewt0rA0mvgZLGf6nB07dkzbtm3TPffc0+P+yZMna/PmzTE/pqqqSgsWLND69es1ZcoUHTx4UC+88II+//nPx/13jh49qqNHj3a939bWlpknAAAAUvf+Bmn3YikQVOIJvE7JikgF4zwTzlLl+rTmoUOHFIlENHz48B73Dx8+XAcOHIj5MVVVVVq9erWmTZumAQMG6Oyzz1ZhYaEWL14c999ZtGiRwuFw123UqFEZfR4AACAFRbVSqNAOXlZHglvEvq6o1u0RO871cBYVCAR6vG9ZVq/7ohobGzVnzhzdf//92rZtm1566SW99dZbmjVrVtzPP3/+fLW2tnbdUp3+BAAAGRQ6Vbrwbkmxf8efEJAuukcKDXViVEZxfVrzjDPOUDAY7FUlO3jwYK9qWtSiRYtUXV2tu+66S5JUXl6uU045RVdeeaUeeeQRjRgxotfHDBw4UAMHDsz8EwAAAOkpmS01PiZ1tMS/JhSWLpjt2JBM4nrlbMCAAaqoqFB9fX2P++vr61VVVRXzYz766CPl5fUcevDjFtgG7G8AAACJJK2e+bdqJhkQziRp3rx5Wr58uVasWKFdu3Zp7ty52rt3b9c05fz581Vbe2LO+brrrlNdXZ2WLVumN998U5s2bdKcOXN02WWXaeTIkW49DQAAkKqS2XZ1LBYfV80kA6Y1JWnatGk6fPiwHnroIe3fv18TJkzQ+vXrNWbMGEnS/v37e/Q8mzFjhtrb2/Uf//Ef+rd/+zcVFhbqM5/5jB577DG3ngIAAEhHtHq28171bKvh76qZZEifMzfQ5wwAAJd1tEvrRvdcexYqlKbu82Q4SzV7GDGtCQAAfKjX2jOqZpIh05oAAMBnoicFdB6V8kJS5zH77ZE3pa23nXSxt04KSIZwBgAAnNf9pIDoCqvO49KeZ0+6MAsnBRh+hBThDAAAOK+oVtp530m9zjolq7P3tZk+KcDwI6RYcwYAAJzn5kkBhh8hReUMALppaJDq6qSWFqmwUKqpkcrK3B4V4FFunRQQt43HydzZoEDlDAAkNTVJ1dVSebn08MPSkiX22/Jy6Yor7McBZJibJwUkaoIb5VIzXMIZAN9rapImTZK2bLHfj0Skjg77rSS9/rr9OAENyAK3Tgow+AgpwhkA35s+XWptPRHGThaJ2I/PmOHosAB/iBuSHAhHhh4hRTgD4GsNDdLmzfGDWVQkIm3aZF8PIMNihSQnwpGbwTABwhkAX6urk4LB1K4NBqW1a7M7HsCX3DwpwK1gmADhDICvtbRIeSn+JMzLk5qbszocwL+6hyQnw5GBR0gRzgD4WmGh1Bmj52UsnZ3SsGFZHQ7gX10hSc6HI7eCYRz0OQPQxY89vmpqpIULU7s2ErGvB5AlpXOk/KHS2JnO/rtdfc/mu141k6SAZVmpHCzlOW1tbQqHw2ptbVVBQYHbwwFc1dRk71jcvNleV5WXZ1eJIhG799fKlVJxsdujzJ7qaruNRqJNAcGgVFkpbdzo3LgAOOj4R9KeFXYwzB+SlX8i1exBOCOcweeiPb7itZIIBqVw2A4vXg1ofA0AOCHV7MGaM8Dn6PFlB64tW+zKmGSHsVDoxC7OykqCGQDnsOYM8LFoj69kuvf48uoatOJie8qyocFul9HcbC/+r6mRJkxwe3QA/IRwBvhYtMdXsgas0okeX14NZ1FlZd5/jgDMxrQm4GP0+AIA8xDOAB+jxxcAmIdwBvhYTU1qU5oSPb4AwCmEM8DHysqkqqrkZ0sGg3YvMBbGA0D2Ec4An1u1yu7hFS+gRXt8rVzp6LAAwLcIZ4DP0eMLAMxCKw0A9PgCkFua35D2LJeUyiFHAWnsLdKw8myPKmMIZwC60OMLQE54f4O0e7EUCCrxJGCnZEWkgnE5Fc6Y1gQAALmlqFYKFdrBy+pIcIvY1xXVuj3itBDOAABAbgmdKl14t6RAkgsD0kX3SKGhTowqYwhnAAAg95TMlkLhxNeEwtIFs50ZTwYRzgAAQO5JWj3LzaqZRDgDAAC5KlH1LEerZhLhDAAA5Kq41bPcrZpJtNIA0tbQINXVSS0t9sHhNTW0nwAA15TMlhofkzpaTtyXw1UziXAGpKypSZo+Xdq82e6en5cndXZKCxfa506uXEkXfQBwXLR6tvNe2U1pc7tqJjGtCaSkqUmaNMk+xkiSIhGpo8N+K0mvv24/3tTk3hgBwLe6rz3L8aqZRDgDUjJ9utTaeiKMnSwSsR+fMcPRYQEApG5rz5TzVTOJaU0gqYYGeyozmUhE2rTJvp41aADgsNI5Uv5QaexMt0fSb4Qz9JvXF8jX1dlrzOJVzboLBu2Dw730/AEgJ+QPkUpvdXsUGUE4Q5/5ZYF8S4v93FIJZ3l5UnNz1ocEAPAw1pyhT/y0QL6w0A6dqejslIYNy+pwAAAeRzhDn/hpgXxNTWpVM8m+rqYmu+MBAHgb4Qxpiy6QTxZYui+Qz2VlZVJVlT11m0gwaE/nTpjgzLgAAN5EOEPaogvkUxFdIJ/rVq2SwuH4zzsYtB9fudLRYQEAPIhwhrRFF8inwisL5IuL7fV1lZX2+8GgFAqdCGuVlfbjXtgAAQBwF7s1kTa/LpAvLpY2brSnadeutUPnsGH2GjOmMgEAmRKwLMtyexBuaGtrUzgcVmtrqwoKCtweTk5paJDKy9O7nvACAPC7VLMH05pIGwvkAQDIHsIZ+oQF8gAAZAfhDH3CAnkAALKDDQHoMxbIAwCQeYQz9FtZGQd9AwCQKUxrAgAAGIRwBgAAYBDCGQAAgEEIZwAAAAYhnAEAABiEcAYAAGAQwhkAAIBBCGcAAAAGIZwBAAAYhHAGAABgEMIZAACAQQhnAAAABiGcAQAAGIRwBgAAYBDCGQAAgEEIZwAAAAYhnAEAABgk3+0BAACAOJrfkPYsl2SlcHFAGnuLNKw826NClhHOAABIhRtB6f0N0u7FUiCoxJNdnZIVkQrGEc48gHAGAEAq3AhKRbXSzvukjhZJkcTXhgrt65HzWHMGAEAqimrtAGRFJKsjwS2SuaAUOlW68G5JgSQXBqSL7pFCQ/v/b8J1hDMAAFLhVlAqmS2FwknGFpYumJ2Zfw+uY1oTMFxDg1RXJ7W0SIWFUk2NVFbm9qgAnyqZLTU+9vE0YxyZDkrRULjzXsVe70bVzGuonAGGamqSqqul8nLp4YelJUvst+Xl0hVX2I8DcFjS6lmWglKi6hlVM88hnAEGamqSJk2Stmyx349EpI4O+60kvf66/TgBDXCBG0EpbiikauZFhDPAQNOnS62tJ8LYySIR+/EZMxwdFgDJvaAUKxRSNfMkwhlgmIYGafPm+MEsKhKRNm2yrwfgMDeCUq9QSNXMqwhngGHq6qRgMLVrg0Fp7drsjgdADG4Fpe6hkKqZZxHOAMO0tEh5Kf6fmZcnNTdndTgA4nEjKHWFQlE18zDCGWCYwkKpszO1azs7pWHDsjocAPG4FZRK50gVi6WS25z59+A4whlgmJqa5OvNoiIR+3oALnEjKOUPkUpvtd/CkwhngGHKyqSqquTrzoJBuw/ahAnOjAtADAQlZAHhDDDQqlVSOBw/oAWD9uMrVzo6LACAAwhngIGKi+0GtJWV9vvBoBQKnQhrlZX248XF7o0RAJAdnK0JGKq4WNq40e5jtnatvStz2DB7jRlTmfC95jekPcsV+6zJkwWksbdIw8qzPSogIwhngOHKyjjoHOjl/Q3S7sVSIKjEk0CdkhWRCsYRzpAzmNYEAOSeolopVGgHL6sjwS1iX1dU6/aIgZQRzgAAuSfu+ZYn44gj5B5jwtnSpUtVVFSkQYMGqaKiQhs2bEh4/dGjR7VgwQKNGTNGAwcO1NixY7VixQqHRgsAcF2s8y1PxhFHyEFGrDlbs2aN7rjjDi1dulTV1dV6+umnNWXKFDU2Nmr06NExP+aGG27Q3/72Nz3zzDMqLi7WwYMHdfz4cYdHDgBwTbR6tvNexd4YQNUMuSlgWVYqW12yatKkSbr00ku1bNmyrvvGjx+vqVOnatGiRb2uf+mll3TjjTfqzTff1GmnnZbSv3H06FEdPXq06/22tjaNGjVKra2tKigo6P+TAAA4r6NdWjda6mjp/VioUJq6j3AGY7S1tSkcDifNHq5Pax47dkzbtm3T5MmTe9w/efJkbd68OebH/OIXv9DEiRP1ve99T+ecc45KSkp055136u9//3vcf2fRokUKh8Ndt1GjRmX0eQAAXBB37RlVM+Qu18PZoUOHFIlENHz48B73Dx8+XAcOHIj5MW+++aY2btyoP/7xj1q7dq0ef/xxvfDCC5o9O/66gvnz56u1tbXrtm/fvow+DwCAS2KtPWOtGXKYEWvOJCkQ6PlXj2VZve6L6uzsVCAQ0OrVqxUO2/9D/uhHP9JXvvIVLVmyRIMHD+71MQMHDtTAgQMzP3AAgLt6rT2jaobc5nrl7IwzzlAwGOxVJTt48GCvalrUiBEjdM4553QFM8leo2ZZlt59992sjhcAYKDu1TOqZshxroezAQMGqKKiQvX19T3ur6+vV1VVVcyPqa6u1nvvvacjR4503bd7927l5eXp3HPPzep4AQAG6lp7JqpmyHmuhzNJmjdvnpYvX64VK1Zo165dmjt3rvbu3atZs2ZJsteL1dae6O5800036fTTT9fNN9+sxsZGvfrqq7rrrrs0c+bMmFOaAAAfKJ0jVSyWSm5zeyRAvxix5mzatGk6fPiwHnroIe3fv18TJkzQ+vXrNWbMGEnS/v37tXfv3q7rhw4dqvr6et12222aOHGiTj/9dN1www165JFH3HoKAAC35Q+RSm91exRAvxnR58wNqfYaAQAAyISc6XMGAACAEwhnAAAABiGcAQAAGMSIDQEAzNfQINXVSS0tUmGhVFMjlZW5PSoA8B7CGYCEmpqk6dOlzZulYFDKy5M6O6WFC6XqamnlSqm42O1RAoB3EM4AF5lejWpqkiZNklpb7fcjEfsW9frr9uNbthDQACBTaKVBKw24IF41KhIxqxpVXW0Hr+6B7GTBoFRZKW3c6Ny4ACAXpZo9qJz5mOlVG6/KlWpUQ4MdHpOJRKRNm+zr+f4BgP4jnPkQa4jcNX26HcziVaMiEfvxGTPcrUbV1dnfH4mqZlHBoLR2LeEMADKBcOYzuVK18apcqka1tNjBPZVwlpcnNTdnfUhA3zS/Ie1ZLimVVTwBaewt0rDybI8KiItw5jO5UrXxqlyqRhUW2hXVVHR2SsOGZXU4QN+9v0HavVgKBJW4vWenZEWkgnGEM7iKcOYjuVS18apcqkbV1NhT3amIROzrgYzJZLWrqFbaeZ/U0SIpyf98oUL7esBFhDMfyaWqjVflUjWqrEyqqkp9t+aECc6NDT6QyWpX6FTpwrulnfcqcdgLSBfdI4WG9n3cQAZwfJOPRKs2qXC7auNVNTWphWPJjGrUqlVSOGwHsFiCQfvxlSsdHRb8oKjWrmJZEcnqSHCLpFbtKpkthcKJrwmFpQtmZ+oZAH1GOPORXKraeFW0GhUv7EQFg/bOWberUcXFduWssvLEuEKhE+OvrGTzCLIkWu1SIMmFKVa7kn4+qmYwB01ofdSEtqFBKk9jjWtDg/vhwIu675iNVUWLVqNMCz0NDfZUd3OzHdxravj+QJZ1tEvrRn+8ViyOUKE0dV9qoSrR50vn8wB9RBNa9MIaIjNEq1EzZtgbL04+IaCy0sxec2VlrEGEw5KuFUuz2hX381E1g1monPmociblbtXGq6hGAUlkutoV6/NRNYNDqJwhplyt2ngV1SggiUxXu3p9PqpmMA+VM59VzrqjagMgJ2S62tX981E1g4OonCEpqjYAckKmq11dn28+VTMYiVYaAADzde9Tlol+ZKVzpIrFUslt/R8bkGGEMwCA+br6lCkz1a78IVLprfZbwDBMawIAckPpHCl/qDR2ptsjAbKKcAYAyA3RahfgcUxrAgAAGITKGYzU0CDV1dmHtRcW2m0+2FkKAPADwhmM0tQkTZ8ubd7cs0HuwoX2QeA0yAUAeB3TmjBG9GipLVvs9yMRqaPjxDFTr79uP97U5N4YAQDINsIZjDF9evwzPyX7/tZW++gpAAC8imlNGKGhwZ7KTCYSsc8EbWhgDRoAAzS/Ie1Zrp7nfsYTkMbeIg0rz/aokOMIZzBCXZ29xixe1ay7YNA+E5RwBsB172+Qdi+WAkElnozqlKyIVDCOcIakCGcwQkuLvfg/lXCWl2cf1u4GdpEC6KGoVtp538eHsif5ARYqtK/PJVQGXUE4gxEKC+1dmano7JSGDcvqcHphFymAmHodyh5PPw9rdwuVQVf0eUPAP/7xDzU0NOijjz7q9dimTZv6NSj4T01NalUzyb6upia74+mOXaQAEup+KHs8mTis3Q1FtXbFz4pIVkeCWyQ3K4OG6lM4e+211zRq1Ch9+tOf1plnnqlHH320x+NTpkzJyODgH2VlUlWVXZVKJBi0K1UTJjgzLoldpACS6DqUPRDnghytmkkpPLeoHH6OBupTOPu3f/s3/fCHP9Thw4e1bds21dXVaebMmer8eF7KslKZmwZ6WrVKCofjB7Rg0H585UrnxhTdRZqsqtd9FykAH0pUPcvVqlmUlyuDhupTOGtsbFRtrV26HDdunF555RUdPHhQX/nKV3Ts2LGMDhD+UVxsTx1WVtrvB4NSKHQirFVW2o87ubYruos0FdFdpAB8KG6FyQMVJS9XBg3Vp3BWUFCgv/71r13vDx48WOvWrdOgQYP0uc99rquCBqSruFjauFF64w3p/vul2bPttw0N9v1OL7qP7iJNhZu7SAEYIFaFySsVJS9XBg3Up3B29dVX69lnn+1xX35+vlavXq2xY8fq73//e0YGB/8qK7ND2f/7f/ZbJ9eYdWf6LlIABulVYfJQRcnLlUEDBaw+LBA7duyYjh8/riFDhsR8fO/evRo9enS/B5dNbW1tCofDam1tVUFBgdvDgaEaGqTyNHaFNzS4FyQBGKCjXVo32u57FiqUpu7zTnDp/tyivPYcsyzV7JFW5ezf//3fJUkDBgyIG8wkGR/MgFSZvIsUgIG6KkzyXkXJy5VBw6QVzh577DHNnTs37uN79+7t94AA05i4ixSAwUrnSBWLpZLb3B5J5nVfe8Zas6xJ64SAuro63XDDDTpy5Ih+/OMfKxCw03N7e7u+853v6Mknn4zZlBbIZdFdpDNm2O0yup8QEInYu0g5IQBAl/whUumtqV2ba8cjdZ2IMJ+qWRalFc6+8IUvaP369friF7+oDz/8UM8++6xWrFihhQsX6oMPPtDMmTOzNU7AVdFdpA0NdruM5mZ78X9NDVOZAPohF49HKp0j5Q+VxvI7P1v6tCFg69atuuqqq5SXl6cjR47o+uuv16OPPqrS0tJsjDEr2BAAAHBdrEX28bD4PudlZUOAJG3fvl333nuvPvzwQ7W3t6u6ulovvPBCTgUzAACMwPFIiCGtcHbTTTfpk5/8pBobG7VixQpt2LBBf/zjHzV16lQdPXo0W2MEAMC7OB4JJ0krnP3P//yPHnjgAf3lL3/RjBkzVF1drV//+tf6/e9/r8997nM6cuRItsYJAIA3cTwSTpLWmrMDBw7o7LPP7nX/n//8Z1199dUaMWKEfve732V0gNnCmjMAgDESrT1jrZlnZGXNWaxgJkmlpaXasGGDmjlYEACA9HE8Errp09masZx33nnasGFDpj4dAAD+4uWD05GWjIUzKX5lDQAAJMHxSPhYRsMZAADoB45HgghnAACYw8sHpyNlaR3fBGc1NEh1dVJLi1RYaB8VVFbm9qgAAFnF8Ui+16fjm7zA5FYaTU3S9OnS5s29D9muruaQbQAAclHWjm9CdjU1SZMmSVu22O9HIlJHh/1Wkl5/3X68qcm9MQIAgOxhWtMw06dLra0nwtjJIhH78RkzpI0bHR0aAMA0zW9Ie5ZLSmUSLCCNvUUaVp7tUaGfCGcGaWiwpzKTiUSkTZvs61mDBgA+9v4GafdiKRBU4smwTsmKSAXjCGc5gGlNg9TV2WvMUhEMSmvXZnc8AADDFdXaxztZEcnqSHCL2NcV1bo9YqSAyplBWlrsxf/xpjS7y8uTOC0LABxk4hRitPXGznuTjIuGtrmEcGaQwkJ7V2YqOjulYcOyOhzj0WoEgKNMnUIsmS01Phb70PSovJB05E1p620JPhFr0kxBKw2DWmk0NEjlafw/0dAgTZiQvfGYilYjAFzR0S6tG504BEWFCqWp+5yrVP3p0STVs4AUSFSP+ThQTlwilXwrCwOERCuNnFRWJlVVJV93FgzaIcSvwYxWIwBc0evsy3hcmEKMdWh6DxZr0nII4cwwq1ZJ4XD8gBYM2o+vXOnosIyRTqsRAMi4pCFI7pyJGTc4BqQRn4tx/8lYk2YSwplhiovtqlBlpf1+MCiFQifCWmWl/bgfp+2irUaSbZiIthr5+telBx+0Pw4AMiJp9czFkBMrOIbCUuUKMwMl4iKcGai42G4w+8Yb0v33S7Nn228bGuz7/RjMpPRajUjSc89JDz9sr+O74gqmOgFkSKLqmZshp1dw/DgoDh5hbqBETOzWNFhZGbsPu0un1YjUc+drdC2aX6uOADIobvsKA0JO952b3YNioh2dVM2MQ+UMOSOdViMnYy0agIyKN4Xodsjpqp6pZ1BMtCbN7UCJXghnyBk1NalXzWLpfuwVAPRLvClEE0JO6RypYrFUclJPM1MDJXohnCFnpNpqJBGOvQKQMd3DjkkhJ3+IVHqr/bY7kwMleiCcIackazWSDMdeAciYeFOIJjM1UKIHwhlySqxWI3lpfBdz7BWAjIo3hWiqXAyUPsRuTeScaKuRhgZ7ivIvf5F++tPUPjYSsdeuAUBGRKcQc0npHCl/qDR2ptsjQRyEM+Ss7q1G3nzTrqgl2jAQDNoVNz8eewUAXXIxUPoM05rwBI69AgB4BeEMnsCxVwAAr2BaE55x8lq05mZ78X9NDVOZAIDcQTiD53DsFQAglzGtCQAAYBAqZwAAoLfmN6Q9y9XzcPd4AtLYW6Rh5dkelS8QzgAAQG/vb5B2L5YCQSWeaOuUrIhUMI5wliGEMwAAnJYLVamiWmnnfVJHi6QETSQlKVRoX4+MIJwBAOC0XKhKRY962nmvEodIDlDPNDYEAADgtKJau9pkRSSrI8Et4m5VqvtB6fFwgHrGEc4AAHBa1wHkgSQXulyVSjpOqmbZQDgDAMANuVKVSjROE8bnQaw5AwDADUnXdPWxKpXpzQZxx0nVLFuMCWdLly7V97//fe3fv18XXXSRHn/8cV155ZVJP27Tpk361Kc+pQkTJmjHjh3ZHygAAJlSMlv603el4+29H8sLSUfelLbe1u3OFMJUNjYblMyWGh/7eOfmx6iaZY0R4WzNmjW64447tHTpUlVXV+vpp5/WlClT1NjYqNGjR8f9uNbWVtXW1uqf/umf9Le//c3BEQMAcp4J7SxCp0rDr5L++ovej3Uel/Y82/2O1MJUNlpg9KqeUTXLpoBlWal8V2bVpEmTdOmll2rZsmVd940fP15Tp07VokWL4n7cjTfeqAsuuEDBYFDr1q1Lq3LW1tamcDis1tZWFRQU9Gf4AIBctHuJtPXW1CtME5dIJd/K/Dg+ek9ad05q14YKpan7koeiPz2aWguMSxZ9vOA/BR3t0rrRduhLdRzoIdXs4fqGgGPHjmnbtm2aPHlyj/snT56szZs3x/24Z599Vnv27NEDDzyQ0r9z9OhRtbW19bgBAHzMlHYWQ0ZKZ1+TwoVpVKuysdmga+emqJplmevh7NChQ4pEIho+fHiP+4cPH64DBw7E/Ji//OUvuueee7R69Wrl56c2M7to0SKFw+Gu26hRo/o9dgBADjOpnUXliuTXpBOmstUCo3SOVLFYKrkt+bXoM9fDWVQg0PMbyLKsXvdJUiQS0U033aQHH3xQJSUlKX/++fPnq7W1teu2b9++fo8ZAJDjTGlnkbR61ocwlY0WGPlDpNJb7bfIGtfD2RlnnKFgMNirSnbw4MFe1TRJam9v19atW3XrrbcqPz9f+fn5euihh7Rz507l5+fr17/+dcx/Z+DAgSooKOhxAwD4nElNVi9fJeUNjv1YX8JU3OfGYn7TuR7OBgwYoIqKCtXX1/e4v76+XlVVVb2uLygoUENDg3bs2NF1mzVrlkpLS7Vjxw5NmjTJqaEDALzAlCarg4dLZfcro2Eq1nOjBYbxXA9nkjRv3jwtX75cK1as0K5duzR37lzt3btXs2bNkmRPSdbW2gsx8/LyNGHChB63s846S4MGDdKECRN0yimnuPlUAAC5xqQKU6bDVK/nRtUsFxgRzqZNm6bHH39cDz30kC655BK9+uqrWr9+vcaMGSNJ2r9/v/bu3evyKAEAnmVKhSkbYar7c6NqlhOM6HPmBvqcAQB66NEbLM0eYJmUjX5if3pU2jlfuuRRd54TJOVQnzMAAIxgSoUpG/3EaIGRU4w4vgkAANd1HVE03/11WaVzpPyh0tiZmfl80RYYyAmEMwAAojIdivqKMOVrhDMAAKIIRTAAa84AAAAMQjgDAAAwCOEMAADAIIQzAAAAgxDOAAAADEI4AwAAMAjhDAAAwCCEMwAAAIMQzgAAAAxCOAMAADAIxzcBAJCq5jekPcslWSlcHJDG3iINK8/2qOAxhDMAAFL1/gZp92IpEFTiyadOyYpIBeMIZ0gb05oAAKSqqFYKFdrBy+pIcIvY1xXVuj1i5CDCGQAAqQqdKl14t6RAkgsD0kX3SKGhTowKHkM4AwAgHSWzpVA48TWhsHTBbGfGA89hzRnS1tAg1dVJLS1SYaFUUyOVlbk9KgBwSLR6tvNexd4YQNUM/ROwLCuVLSee09bWpnA4rNbWVhUUFLg9nJzQ1CRNny5t3iwFg1JentTZKUUiUnW1tHKlVFzs9igBwAEd7dK60VJHS+/HQoXS1H2EM/SSavZgWhMpaWqSJk2Stmyx349EpI4O+60kvf66/XhTk3tjBADHxF17RtUM/Uc4Q0qmT5daW0+EsZNFIvbjM2Y4OiwAcE+stWesNUMGEM6QVEODPZUZL5hFRSLSpk329QDgeb2qZ1TNkBmEMyRVV2evMUtFMCitXZvd8QCAMbpXz6iaIUMIZ0iqpcVe/J+KvDypuTmrwwEAc3RVz0TVDBlDKw0kVVho78pMRWenNGxYVocDAGYpnSPlD5XGznR7JPAIKmdIqqYm+XqzqEjEvh4AfCN/iFR6q/0WyAAqZ0iqrEyqqrLbaCQKacGgVFkpTZjQ/3+TRrcAAL+iCS1NaFMS7XMWr51GMCiFw3aA608jWhrdAgC8iia0yKjiYjt4VVba7weDUih0YhdnZWVmghmNbgEAfse0JlJWXCxt3GhPOa5da+/KHDbMnnLMxFRmOo1uN27s/78HQFLzG9Ke5Yp9RuTJAtLYW6Rh5dkeFeBrhDOkraws8+u/oo1uk+ne6JY1aEAGvL9B2r1YCgSVeDKlU7IiUsE4whmQZUxrwgg0ugVcUlRrH9RtRSSrI8EtYl9XVOv2iAHPI5zBCDS6BVwS9wDvk3E0EeAUwhmMQKNbwEWxDvA+GUcTAY4hnMEINLoFXJS0ekbVDHAS4QxGiDa6TbbuLBi0+51lYncogG4SVc+omgGOIpzBGKtW2Y1s4wW0aKPblSsdHRbgD3GrZ1TNAKcRzmAMJxrdAkggVvWMqhngOPqcwSjZbnQLIIFo9WznvbKb0lI1A9zA2Zo+PVuTg8UBxNTRLq0bLXW02H3Npu4jnAEZkmr2oHLmM/EOFl+4kIPFAahb9Ww+VTPAJYQzH4keLN7aar8fifRsXxE9WJx1XYDPlc6R8odKY2e6PRLAlwhnPsLB4vAcDu3OjvwhUumtbo8C8C3CmU9wsDg8iUO7nUUYBhxBOPOJ6MHiqXThjx4sTjiD8YpqpZ332YvXleSbm0O7+48wDDiCPmc+wcHi8CQO7XZWUa0dcq2IZHUkuEUIw0A/EM58goPF4Vkc2u0cwjDgCMKZT3CwODyLQ7udRRgGso5w5hMcLA5P49Bu5xCGgawjnPkIB4vDszi021mEYSCrCGc+wsHi8DQO7XYOYRjIKlpp+AwHi8OzOLTbWSWzpcbHPm5j8jHCMJARhDOfKiujjxk8qHtgIChkF2EYyBqmNQF4R9d0mwgKTug+lUwYBjKGcAbAW0rnSBWLpZLb3B6J9xGGgaxgWhOAt3Bot7NK50j5Q6WxM90eCeAZhDMAQN8RhoGMY1oTAADAIFTOAMBrmt+Q9iyXvYsymYA09hZpWHm2RwUgRYQzAPCa9zdIuxdLgaAST5B0SlZEKhhHOAMMwrQmAHhNUa0UKrSDl9WR4BaxryuqdXvEALohnAGA1yQ9nDyKxrGAiQhnAOBFiQ4nj6JxLGAkwhkAeFHS6hlVM8BUhDMA8KpE1TOqZoCxCGcA4FVxq2dUzQCTEc4AwMtiVc+omgFGI5wBgJf1qp5RNQNMRzgDAK/rXj2jagYYj3AGAF7XVT0TVTMgB3B8EwD4QekcKX+oNHam2yPJHlPOFDVlHMhZhDMA8IP8IVLprW6PIrtMOVPUlHEgZzGtCQDwBlPOFDVlHMhZhDMAgDeYcqaoKeNAziKcAQC8w5QzRU0ZB3IS4QwA4B2mnClqyjiQkwhnAABvMeVMUVPGgZxDOAMAeIspZ4qaMg7kHMIZAMB7TDlT1JRxIKcQzgAA3mPKmaKmjAM5hXAGAPAmU84UNWUcyBmcEAAA8KZo1WrnfHerVaaMIxGOnDIK4QwA4F2mnClqyjji4cgpowQsy0olJntOW1ubwuGwWltbVVBQ4PZwAABwT0e7tG601NGS/NpQoTR1n5kVQMOlmj1YcwYAgN9x5JRRCGcAAIAjpwxCOAMAABw5ZRBjNgQsXbpU3//+97V//35ddNFFevzxx3XllVfGvLaurk7Lli3Tjh07dPToUV100UVauHChrrnmGodHDSBnsBsNSK5kttT4WOy1Z1TNHGNEOFuzZo3uuOMOLV26VNXV1Xr66ac1ZcoUNTY2avTo0b2uf/XVV/XZz35W3/3ud1VYWKhnn31W1113nbZs2aJPfOITLjwDAMZjNxqQXFfbj3vV8w8ZqmZOMmK35qRJk3TppZdq2bJlXfeNHz9eU6dO1aJFi1L6HBdddJGmTZum+++/P6Xr2a0J+Ay70YDUxPp/hf8nMiJndmseO3ZM27Zt0+TJk3vcP3nyZG3evDmlz9HZ2an29naddtppca85evSo2traetwA+Ai70YDUcOSU61wPZ4cOHVIkEtHw4cN73D98+HAdOHAgpc/xwx/+UB9++KFuuOGGuNcsWrRI4XC46zZq1Kh+jRtADmI3GpAajpxylevhLCoQ6PnXrGVZve6L5fnnn9fChQu1Zs0anXXWWXGvmz9/vlpbW7tu+/bt6/eYAeQYdqMBqen6f0X8P+EC1zcEnHHGGQoGg72qZAcPHuxVTTvZmjVr9I1vfEM///nPdfXVVye8duDAgRo4cGC/xwsgx7EbDUiN6UdOeZjr4WzAgAGqqKhQfX29vvSlL3XdX19fry9+8YtxP+7555/XzJkz9fzzz+vzn/+8E0MFYKp022QU1do7N9mNBsSXP0QqvdXtUfiS6+FMkubNm6evf/3rmjhxoi6//HL9+Mc/1t69ezVr1ixJ9pTkX//6V/3nf/6nJDuY1dbW6oknnlBlZWVX1W3w4MEKh5OsJwHgPem2yfjED+0qWY/daFTNAJjBiDVn06ZN0+OPP66HHnpIl1xyiV599VWtX79eY8aMkSTt379fe/fu7br+6aef1vHjxzV79myNGDGi63b77be79RQAuKmo1t7qb0UkqyPBLWJfV/z/sRsNgLGM6HPmBvqcAR7zp0djNM48WUC6ZJEdzLr3cqKHEwAH5EyfMwDIiHTbZLAbDYChCGcAvKEvbTJK50gVi6WS25wYIQCkhHAGwDsSVc9iLfiP7kbLH5L9sQFAighnALwjbvWMBf8AcgfhDIC3xKqe0SYDQA4hnAHwFg5tBpDjCGcAvIdDmwHkMCNOCADwsXSPIRp7izSsPNujyj3R6tnO+VTNAOQcwhlgknSPISoYRziLx41DmwnXfcPXDeiBcAaYpKhW2nnfx2c+RhJfGyq0r0dsbhzaTLjuG75uQA+sOQNMkrSRahSL3I2U7hmfhGsbXzegB8IZYJp0jyGCOQjXfcPXDeiBcAaYpi/HEMEchOu+4esGdCGcASZK9xgimINw3Td83YAuhDPARBxDlNsI133D1w2QRDgDzMUxRLmLcN03fN0ASYQzwFwcQ5TbCNd9w9cNIJwBRuMYotxFuO4bvm4A4QwwWtcvKvELKhcRrvuGrxt8jnAGmK50jlSxWCq5ze2RIF2E677h6waf4/gmwHRuHEOEzHHjjE8v4OsGHyOcAUA2Ea77hq8bfIxpTQAAAIMQzgAAAAxCOAMAADAI4QwAAMAghDMAAACDEM4AAAAMQjgDAAAwCH3OAADe0fyGtGe5JCuFiwPS2FukYeXZHhWQFsIZAMA73t8g7V4sBYJKPDnUKVkRqWAc4QzGIZwByByqFnBbUa208z6po0VSJPG1oUL7esAwhDMAmUPVov8IuP0TPTR9571K/DUMcKg6jEU4A5A5VC36j4DbfyWzpcbHPv4+jCMUli6Y7diQgHSwWxNA5kSrFgokuZCqRVxFtXZwtSKS1ZHgFiHgxpP0+5DvP5iNcAYgs0pm21WJRKhaxEfAzYxE34d8/8FwhDMAmUXVov8IuP0X9/uQ7z+YjzVnADIv0ZofQkVySRe1EzBSEuv70JTvPzZ+IAHCGYDMixsuCBUpI+D2X6/vQ4O+/9j4gQSY1gSQHbGm5ggVqWNaLjO6fx+a9P3Hxg8kQDgDkB29wgWhIm0E3P7r+j6UWd9/bPxAAkxrAsie7lNzhIr0mTwtl0tK50j5Q6WxM90eSU9e7cfGerp+I5wByJ6ucDGfUNFXBNz+yx8ild7q9ih68+rGD9bT9RvTmgCyq3SOVLFYKrnN7ZHkJlOn5ZAZXuzHxnq6fqNyBiC7TK1a5BJTp+XQf4mqZ+EL7apzUoZNDXK+ab8FLMtKZVLYc9ra2hQOh9Xa2qqCggK3hwMA8KuOdmnd6J5rz/IGS51/T31qcOISqeRbWR5oGmI9p5OFCqWp+3wVzlLNHkxrAgDgprg7mwtzd2qQk0L6hXAGAIDbTu7HNm5e7rfa8OJ6OocQzgAAcFusjR+5fsYqjZT7jHAGAIAJTt7Z7IWpQRop9wnhDAAAE0R3NucPOXFfrk8NclJInxDOAAAwlRemBk0939RghDMAAEyW61ODNFJOG+EMAACTeWFqkJNC0kI4AwDAdLk+NRhrPR3iIpwBAGA6pgZ9hbM1AQDIBZyx6huEMwAAckF0ahCex7QmAACAQQhnAAAABiGcAQAAGIRwBgAAYBDCGQAAgEEIZwAAAAYhnAEAABiEcAYAAGAQmtACANAfzW9Ie5ZLslK4OCCNvUUaVp7tUSGHEc4AAOiP9zdIuxdLgaAST0h1SlZEKhhHOENCTGsCANAfRbVSqNAOXlZHglvEvq6o1u0Rw3CEMwAA+iN0qnTh3ZICSS4MSBfdI4WGOjEq5DCmNQEA6K+S2VLjY1JHS/xrQmHpgtmODanPWEPnOsIZAAD9Fa2e7bxXsUNNDlXNWEPnOqY1AQDIhJLZdnUsllypmkmsoTMA4QwAgEyIu/Ysh6pmEmvoDEA4AwAgU2JVz3KpahaVqAoYlYvPK0cQzgAAyJReVaccrS4lrZ7l6PPKEYQzAAAyqXvVKZerS15ZQ5eDCGcAAGRSV9VJuV1d8soauhxEOAMAINNK50gVi6WS29weSf94ZQ1djqHPGQAAmZY/RCq91e1R9F+v/m0uV8180iCXcAbAe3zyAxxwRPfTD9yumvmkQS7hDID3+OQHOOCIrurZfPfXmhXVSjvv+/iYrEjia3O4QS5rzgB4Dx3OgcwyZQ2dTxrkEs4AeI9PfoADjomuocsf4vZIfNEgl3AGwJt88AMc8CUfNMglnAHwJh/8AAd8y+MNcglnALzL4z/AAd/yeINcwhkA7/L4D3DA1zzcIJdwBsDbPPwDHPA1rxwyHwPhDIC3efgHOOB7Xjlk/iTGhLOlS5eqqKhIgwYNUkVFhTZs2JDw+ldeeUUVFRUaNGiQzj//fD311FMOjRRAzvHoD3DA97xyyPxJjAhna9as0R133KEFCxZo+/btuvLKKzVlyhTt3bs35vVvvfWWrr32Wl155ZXavn277r33Xs2ZM0cvvviiwyMHkBM8+gMcgMxpkJtBAcuyUjl8LqsmTZqkSy+9VMuWLeu6b/z48Zo6daoWLVrU6/q7775bv/jFL7Rr166u+2bNmqWdO3fqtddei/lvHD16VEePHu16v62tTaNGjVJra6sKCgoy+GwAGOn4R9KeFdLYmWY00gTgO21tbQqHw0mzh+uVs2PHjmnbtm2aPHlyj/snT56szZs3x/yY1157rdf111xzjbZu3aqOjo6YH7No0SKFw+Gu26hRozLzBADkBpM6nANAAq6Hs0OHDikSiWj48OE97h8+fLgOHDgQ82MOHDgQ8/rjx4/r0KFDMT9m/vz5am1t7brt27cvM08AAAAgg/LdHkBUINCzD5FlWb3uS3Z9rPujBg4cqIEDB/ZzlAAAANnleuXsjDPOUDAY7FUlO3jwYK/qWNTZZ58d8/r8/HydfvrpWRsrAABAtrkezgYMGKCKigrV19f3uL++vl5VVVUxP+byyy/vdf3LL7+siRMnKhQKZW2sAAAA2eZ6OJOkefPmafny5VqxYoV27dqluXPnau/evZo1a5Yke71YbW1t1/WzZs3SO++8o3nz5mnXrl1asWKFnnnmGd15551uPQUAAICMMGLN2bRp03T48GE99NBD2r9/vyZMmKD169drzJgxkqT9+/f36HlWVFSk9evXa+7cuVqyZIlGjhypJ598Ul/+8pfdegoAAAAZYUSfMzek2msEAAAgE3KmzxkAAABOIJwBAAAYhHAGAABgEMIZAACAQQhnAAAABiGcAQAAGIRwBgAAYBDCGQAAgEGMOCHADdHeu21tbS6PBAAA+EE0cyTr/+/bcNbe3i5JGjVqlMsjAQAAftLe3q5wOBz3cd8e39TZ2an33ntPp556qgKBQEY+Z1tbm0aNGqV9+/ZxJJQBeD3MwuthHl4Ts/B6mCUbr4dlWWpvb9fIkSOVlxd/ZZlvK2d5eXk699xzs/K5CwoK+B/LILweZuH1MA+viVl4PcyS6dcjUcUsig0BAAAABiGcAQAAGIRwlkEDBw7UAw88oIEDB7o9FIjXwzS8HubhNTELr4dZ3Hw9fLshAAAAwERUzgAAAAxCOAMAADAI4QwAAMAghDMAAACDEM4AAAAMQjhL09KlS1VUVKRBgwapoqJCGzZsSHj9K6+8ooqKCg0aNEjnn3++nnrqKYdG6g/pvB51dXX67Gc/qzPPPFMFBQW6/PLL9b//+78Ojtb70v3/I2rTpk3Kz8/XJZdckt0B+ky6r8fRo0e1YMECjRkzRgMHDtTYsWO1YsUKh0brD+m+JqtXr9bFF1+sIUOGaMSIEbr55pt1+PBhh0brXa+++qquu+46jRw5UoFAQOvWrUv6MY7+PreQsp/97GdWKBSyfvKTn1iNjY3W7bffbp1yyinWO++8E/P6N9980xoyZIh1++23W42NjdZPfvITKxQKWS+88ILDI/emdF+P22+/3Xrssces3/3ud9bu3but+fPnW6FQyPrDH/7g8Mi9Kd3XI6qlpcU6//zzrcmTJ1sXX3yxM4P1gb68Htdff701adIkq76+3nrrrbesLVu2WJs2bXJw1N6W7muyYcMGKy8vz3riiSesN99809qwYYN10UUXWVOnTnV45N6zfv16a8GCBdaLL75oSbLWrl2b8Hqnf58TztJw2WWXWbNmzepx37hx46x77rkn5vXf/va3rXHjxvW471//9V+tysrKrI3RT9J9PWK58MILrQcffDDTQ/Olvr4e06ZNs+677z7rgQceIJxlULqvx69+9SsrHA5bhw8fdmJ4vpTua/L973/fOv/883vc9+STT1rnnntu1sboR6mEM6d/nzOtmaJjx45p27Ztmjx5co/7J0+erM2bN8f8mNdee63X9ddcc422bt2qjo6OrI3VD/ryepyss7NT7e3tOu2007IxRF/p6+vx7LPPas+ePXrggQeyPURf6cvr8Ytf/EITJ07U9773PZ1zzjkqKSnRnXfeqb///e9ODNnz+vKaVFVV6d1339X69etlWZb+9re/6YUXXtDnP/95J4aMbpz+fZ6f8c/oUYcOHVIkEtHw4cN73D98+HAdOHAg5sccOHAg5vXHjx/XoUOHNGLEiKyN1+v68nqc7Ic//KE+/PBD3XDDDdkYoq/05fX4y1/+onvuuUcbNmxQfj4/ijKpL6/Hm2++qY0bN2rQoEFau3atDh06pG9961v64IMPWHeWAX15TaqqqrR69WpNmzZN//jHP3T8+HFdf/31Wrx4sRNDRjdO/z6ncpamQCDQ433Lsnrdl+z6WPejb9J9PaKef/55LVy4UGvWrNFZZ52VreH5TqqvRyQS0U033aQHH3xQJSUlTg3Pd9L5/6Ozs1OBQECrV6/WZZddpmuvvVY/+tGPtHLlSqpnGZTOa9LY2Kg5c+bo/vvv17Zt2/TSSy/prbfe0qxZs5wYKk7i5O9z/lxN0RlnnKFgMNjrL5yDBw/2StNRZ599dszr8/Pzdfrpp2dtrH7Ql9cjas2aNfrGN76hn//857r66quzOUzfSPf1aG9v19atW7V9+3bdeuutkuxwYFmW8vPz9fLLL+szn/mMI2P3or78/zFixAidc845CofDXfeNHz9elmXp3Xff1QUXXJDVMXtdX16TRYsWqbq6WnfddZckqby8XKeccoquvPJKPfLII8y+OMjp3+dUzlI0YMAAVVRUqL6+vsf99fX1qqqqivkxl19+ea/rX375ZU2cOFGhUChrY/WDvrwekl0xmzFjhp577jnWbWRQuq9HQUGBGhoatGPHjq7brFmzVFpaqh07dmjSpElODd2T+vL/R3V1td577z0dOXKk677du3crLy9P5557blbH6wd9eU0++ugj5eX1/DUdDAYlnajawBmO/z7PyjYDj4pug37mmWesxsZG64477rBOOeUU6+2337Ysy7Luuece6+tf/3rX9dGtt3PnzrUaGxutZ555hlYaGZTu6/Hcc89Z+fn51pIlS6z9+/d33VpaWtx6Cp6S7utxMnZrZla6r0d7e7t17rnnWl/5ylesP/3pT9Yrr7xiXXDBBdYtt9zi1lPwnHRfk2effdbKz8+3li5dau3Zs8fauHGjNXHiROuyyy5z6yl4Rnt7u7V9+3Zr+/btliTrRz/6kbV9+/autiZu/z4nnKVpyZIl1pgxY6wBAwZYl156qfXKK690PTZ9+nTrU5/6VI/rf/vb31qf+MQnrAEDBljnnXeetWzZModH7G3pvB6f+tSnLEm9btOnT3d+4B6V7v8f3RHOMi/d12PXrl3W1VdfbQ0ePNg699xzrXnz5lkfffSRw6P2tnRfkyeffNK68MILrcGDB1sjRoyw/vmf/9l69913HR619/zmN79J+PvA7d/nAcuiNgoAAGAK1pwBAAAYhHAGAABgEMIZAACAQQhnAAAABiGcAQAAGIRwBgAAYBDCGQAAgEEIZwAAAAYhnAEAABiEcAYAMezfv19Dhw7VjTfe2OP+X/7ylwqFQlqwYIFLIwPgdYQzAIhhxIgR+va3v63/+q//0rZt2yRJv/3tb/XVr35V3/zmN/Wd73zH5REC8CrCGQDEceedd2rEiBG6++679fvf/17XX3+9vva1r+mJJ57oumbZsmW69NJLFQqFtHDhQvcGC8AzCGcAEMeQIUP0yCOP6P/+7/901VVXacqUKfrJT36iQCDQdc2IESP04IMPaurUqe4NFICn5Ls9AAAwWUlJiSQpEAho5cqVCgaDPR6PhrL//u//dnpoADyKyhkAxLFjxw594QtfUHV1tY4cOaIVK1a4PSQAPkA4A4AY/vznP+uaa67R5Zdfrt/85jf64he/qIULF6q1tdXtoQHwOMIZAJzk7bff1tVXX63S0lK9+OKLCoVCevTRR9Xc3Kzvfve7bg8PgMcRzgCgm/379+vqq6/WWWedpV/+8pcaPHiwJGncuHGaOXOmnnjiCb399tvuDhKApwUsy7LcHgQA5Krjx4/r+PHj+uY3v6lzzjlH9913n0KhUK+NAwCQKsIZAPTDwoUL9eCDD/a479lnn9WMGTPcGRCAnEc4AwAAMAhrzgAAAAxCOAMAADAI4QwAAMAghDMAAACDEM4AAAAMQjgDAAAwCOEMAADAIIQzAAAAgxDOAAAADEI4AwAAMAjhDAAAwCD/P//C2paOdujUAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Plot the data\n",
    "plt.figure(figsize = (7,7))\n",
    "\n",
    "plt.scatter(X[y == 0,0],\n",
    "            X[y == 0,1],\n",
    "            c = \"blue\",\n",
    "            s = 60,\n",
    "            label=\"0\")\n",
    "plt.scatter(X[y == 1,0],\n",
    "            X[y == 1,1],\n",
    "            c = \"orange\",\n",
    "            marker = 'v',\n",
    "            s = 60,\n",
    "            label=\"1\")\n",
    "\n",
    "plt.legend(fontsize = 12)\n",
    "plt.xlabel(\"$x_1$\",fontsize = 12)\n",
    "plt.ylabel(\"$x_2$\",fontsize = 12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's split our data into training and test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y,\n",
    "                                    shuffle=True,\n",
    "                                    random_state=435,\n",
    "                                    test_size=.3,\n",
    "                                    stratify=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 1.0\n",
      "Test F1 score is 1.0\n"
     ]
    }
   ],
   "source": [
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Messier Data Set\n",
    "\n",
    "Let's try this again on a messier data set with a less obvious decision boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate the random data\n",
    "np.random.seed(440)\n",
    "n_rows = 100\n",
    "diff = .1\n",
    "X = np.random.random((n_rows,2))\n",
    "X_prime = X[(X[:,1] - X[:,0]) <= -diff,:]\n",
    "X_2prime = X[(X[:,1] - X[:,0]) >= diff,:]\n",
    "X_3prime = [[.4,.9],[.6,.45],[.7,.9],[.3,.19],[.1,.4]]\n",
    "\n",
    "del X\n",
    "X = np.append(X_prime,np.append(X_2prime,X_3prime,axis = 0),axis=0)\n",
    "\n",
    "y = np.empty(np.shape(X)[0])\n",
    "y[(X[:,1] - X[:,0]) <= -diff] = 1\n",
    "y[((X[:,1] - X[:,0]) >= diff)] = 0\n",
    "y[-5] = 1\n",
    "y[-4] = -1\n",
    "y[-3] = 1\n",
    "y[-2] = -1\n",
    "y[-1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJbCAYAAAChCjPHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOjUlEQVR4nO3de3RU9b3//9dkMtzETPCGoICRSEBJtIZKTHRZW4vFVkvTVqxnNSD1/A4tikK1ingULy3a21flANoiwllFy6kGTk8X9Zh12ioXoYUCThtWafACViiCuaEthMn+/bGdkJCZyUwyM/uz934+1poVMrMTPnNJ8pr357Pfn4BlWZYAAABghDynBwAAAIATCGcAAAAGIZwBAAAYhHAGAABgEMIZAACAQQhnAAAABiGcAQAAGCTf6QE4pb29Xe+9955OPfVUBQIBp4cDAAA8zrIstba2avjw4crLS1wf8204e++99zRixAinhwEAAHxm3759OvfccxPe7ttwduqpp0qyH6CCggKHRwMAALyupaVFI0aM6Mggifg2nMWmMgsKCghnAAAgZ3paTsUJAQAAAAYhnAEAABiEcAYAAGAQwhkAAIBBfHtCAAAA6JtoNKq2tjanh+G4UCikYDCYse9HOAMAAGmxLEsHDhxQc3OzLMtyejiOCwQCCofDOvvsszPS2J5wBgAA0tLc3KympiadeeaZOuWUU3y9045lWfrwww/1/vvva+DAgSosLOzz9yScAQCAlFmWpYMHD6qgoEBnnHGG08MxwsCBA3X06FEdPHhQ4XC4z2GVEwIAAEDKotGootEoDdxPUlBQ0PHY9BXhDAAApOz48eOSpPx8Jt86iz0escenLwhnAAAgbX5eZxZPJh8PwhkAAIBBCGcAAAAGIZwBAAB87MiRI7rzzjs1fPhwDRgwQJdccol+/vOf53QMRoSz1157Tddff72GDx+uQCCgtWvX9vg1r776qsrLyzVgwACdf/75evrpp7M/UAAAkBWRiPTQQ9KcOfbHSMSZcVRXV2vlypV68MEH9etf/1qf/OQn9bWvfU3PP/98zsZgxKkWH374oS6++GLdcsst+vKXv9zj8W+99Zauu+46/eu//qt+9rOfaePGjfrWt76lM888M6WvBwAAZmhokKZNkzZtkoJBKS9Pam+XFiyQqqqkFSuk4uLcjGXdunWqq6vT888/r6997WuSpKuvvlrvvPOO7r77bk2dOjWj2zQlYkQ4mzx5siZPnpzy8U8//bRGjhypJ554QpI0btw4bd26VT/84Q8JZ3BUJCLV1kpNTVJhoVRdLZWWOj0qADBTQ4M0caLU3Gx/Ho3al5jNm+3bt2zJTUBbs2aNBg8erK9+9atdrr/lllt08803a8uWLaqsrMz6OIyY1kzX66+/rkmTJnW57tprr9XWrVsTbsB69OhRtbS0dLkAmdLQYL/DKyuTHnlEWrzY/lhWJl1xhX07AKCradPsYJaob2s0at8+fXpuxvOnP/1J48aN69bDraysrOP2XHBlODtw4ICGDh3a5bqhQ4fq+PHjOnToUNyvWbhwocLhcMdlxIgRuRgqfCD2zm/LFvvzaFRqazvxyyb2zo+ABqeZsqYHkOzX36ZNiYNZTDQqbdyYm9fr4cOHddppp3W7Pnbd4cOHsz8IGTKt2RsnN3uzLCvu9THz5s3T3LlzOz5vaWkhoCEj0nnnt2FDTocGSDJrTY+fsMwhudpa+/WYym5HwaC0Zk1uHr9kzWRz1XjXleHs7LPP1oEDB7pcd/DgQeXn5+v000+P+zX9+/dX//79czE8+EjsnV9POr/z45czcsm0NT1+QBhOTVOT/dikEs7y8qTGxqwPSaeffnrc6tgHH3wgSXGratngymnNyy+/XHV1dV2ue+WVVzRhwgSFQiGHRgU/ir3zS0XsnR+QS6at6fE6ljmkrrDQDq2paG+XhgzJ6nAkSaWlpdq1a1e3/TEjH8+pjh8/PvuDkCHh7MiRI9qxY4d27NghyW6VsWPHDu3du1eSPSVZU1PTcfzMmTP1zjvvaO7cudq1a5eWL1+uZ599VnfddZcTw4ePxd75pSJX7/yAGBPX9HgdYTh11dWpVc0k+7jq6uyOR5K+9KUv6ciRI3rppZe6XL9y5UoNHz5cEydOzP4gZMi05tatW3X11Vd3fB5bGzZt2jStWLFC+/fv7whqklRUVKR169Zpzpw5Wrx4sYYPH66nnnqKNhrIORPf+QExpq7p8SqWOaSntFSqrLSrjMleo8GgVFEh5aJoNXnyZH32s5/VN7/5TbW0tKi4uFgvvPCCXn75Zf3sZz/LSY8zyZBw9qlPfapjQX88K1as6HbdVVddpT/+8Y9ZHBXQs+pqex1JKnL1zg+IMXFNj5cRhtO3cuWJNZHxHrdgUAqH7XV6uVJbW6v58+frgQce0AcffKCxY8fqhRde0E033ZSzMRgxrQm4VeydX09vpoJBeyFwjpYrAJKo7OYayxzSV1xsV84qKuzPg0EpFDrxO7WiIvcnqwwePFhPPvmk9u/fr6NHj2rnzp05DWaSIZUzwM1MfOcHSFR2c40w3DvFxXaboUjEriY2NtqPTXW1f9/QEs6APoq985s+3V5H0vnU+WjUfufHqfNwgolreryMMNw3paVM88YQzmAUtzZt5J0fTEVlN3cIw8gUwhmMkKxp44gR0mc+I513nvlhjXd+MA2V3dwiDCMTCGdwXE8dzPfts3+R5eXRYRvoDSq7uUMYRiYQzuC4npo2xsQW2rLdDNA7VHZzgzCMviKcwVGpNm3sjI3EAbgBYRi9RZ8zOCqdvSk7Y7sZAIBXEc7gqHSaNp6MjcQBAF5EOIOj0mnaeDI6bAMAvIhwBkdVV6e2D108dNgGAGRaa2urvvOd72jSpEk688wzFQgEtCDV7sIZwgkBcFSqTRvjocM2ALhY4xvSnmWSrBQODkijb5WGlGV7VDp8+LB+8pOf6OKLL9aUKVO0bNmyrP+fJyOcwXE9NW2Mhw7bAOBy76+Xdi+SAkEln8hrl6yoVDA2J+Fs1KhRamxsVCAQ0KFDhxwJZ0xrwnGxpo0VFfbnPZ29SYdtAPCAohopVGgHL6stySVqH1dUk5NhBQIBBQKBnPxfiRDOYIRY08Y33pAeeMDuYTZihH1bMCiFQidCW0UFDWgBwPVCp0oX3iOppyAUkC66VwoNzsWojMC0JoxyctNGOmwDgIeNmSXVPy61NSU+JhSWLpiVsyGZgHAGo9FhGwA8LFY923mf4p8Y4L+qmcS0JgAAcNKYWXZ1LB4fVs0kwhkAAHBSwrVn/qyaSYQzAADgtHjVM59WzSTWnAEAAKd1W3vmbNXs17/+tT788EO1trZKkurr6/Xiiy9Kkq677joNGjQoq/9/wLKsVFrzek5LS4vC4bCam5tVUFDg9HAAAHCFf/7zn3rrrbdUVFSkAQMGZO4bt7VKa0faZ26GCqUp+xwLZ+edd57eeeeduLe99dZbOu+887pdn8rjkmr2YFoTAAA4r2PtmRxfa/b222/Lsqy4l3jBLNOY1jRQJCLV1kpNTVJhod3bi3YSAADPK5kt5Q+WRs9weiSOIpwZpKFBmjZN2rTJ7oaflye1t0sLFkhVVfZ2RXTFBwB4Vv4gqeQ2p0fhOKY1DdHQYG/+vWWL/Xk0KrW1ndgIfPNm+/aGBufGCAAAso9wZohp06Tm5hNh7GTRqH379Ok5HRYAAMgxwpkBIhF7KjNRMIuJRqWNG+3jAQCANxHODFBba68xS0UwaG8EDgCAk3zaiSuhTD4ehDMDNDXZi/9TkZcnNTZmdTgAACQUCoUkSR999JHDIzFL7PGIPT59wdmaBigstM/KTEV7uzRkSFaHAwBAQsFgUIWFhTp48KAkadCgQQoETt4X0z8sy9JHH32kgwcPqrCwUMFUp8KSIJwZoLrabpeRimjUPh4AAKecffbZktQR0CAVFhZ2PC59RTgzQGmpVFlpt9FIdlJAMChVVEjjx+dubAAAnCwQCGjYsGE666yz1NbW5vRwHBcKhTJSMYshnBli5Uq7j1midhrBoBQO241os4FdCQAA6QoGgxkNJbBxQoAhiovtyllFhf15MCiFQifO4qyosG/P9A4BDQ327gNlZdIjj0iLF9sfy8qkK66g6S0AALkWsHx6LmyqO8M7IRKx22U0NtqL/6urszOVGduVoKdqXTZCIQAAfpNq9mBa00ClpbmZUkxnV4ING7I/HnhI4xvSnmWSUnnvF5BG3yoNKcv2qOAFvLbgA4Qzn4rtStCTzrsSsAYNKXt/vbR7kRQIKvnqiXbJikoFY/kDitTw2oIPsObMp9iVAFlVVCOFCu0/jlZbkkvUPq6oxukRwy14bcEHCGc+xa4EyKrQqdKF90jqqTFlQLroXik0OBejghfw2oIPEM58il0JkHVjZkmhcPJjQmHpglm5GQ+8g9cWPI5w5lPV1ckb3nbGrgTolR4rHFQ20Eu8tuBxhDOfiu1K0NO6s2DQ7oPGrgTolWQVDiob6AteW/AwwlmWRCLSQw9Jc+bYHyMRp0fU3cqVdh+zRAEt27sSwAcSVjiobKCPeG3Bw2hCm+EmtA0Ndv+wTZvscJOXZ6/ZikbtCtSKFWY1dG1osPuYbdzojvHChdpapbUjpbamE9eFCqUp+/gDir7htQWXoQmtAzp33JfsgNN5XdfmzfbtJnXcLy62G8zmalcC+FCswrHzPtmNQ6lsIEN4bcGjqJxlsHJWVWUHr2QL7YNBe59MOu7DVzpXOKhsIJN4bcFFUs0erDnLkFjH/Z7OgOzccR/wjY71QaKygczitQUPYlozQ2Id91NpTxHruM92SPCVktlS/mBp9AynRwKv4bUFjyGcZUis434q4YyO+/Cl/EFSyW1OjwJexGsLHsO0ZobQcR8AAGQC4SxD6LgPAAAygXCWIXTcBwAAmUA4yyA67gMAgL4inGVQcbHd56yiwv48GJRCoRNhraLCrAa0AADAPJytmWF03AcAAH1BOMuS0lL6mCG7IhG7v15Tk322cHU1rzkA8ALCGeAyDQ3StGn2jhSdN6tfsIDN6gHAC1hzBrhIQ4M0caK9dlGy27K0tZ1o47J5s317Q4NzYwQA9A3hDHCRadOk5ubEPfWiUfv26dNzOiwAQAYxrQm4RCRiT2X2JBqVNm60j2cNGrItlbWPrI8E0kM4A1yittZeY5bKThTBoH22MH8AkS2prH2UWB8J9AbhDHCJpib7j1sq4Swvz27jAmRDbO1jc7P9eTTa9XW5ebM0YYL97yNHEh8TWz9JQAO6Ys0Z4BKFhXbVIRXt7XZ/PSAbUl37yPpIoHcIZ4BLVFenVjWT7OOqq7M7HvhTbO1jqq/FZDqvjwRwAuEMcInSUqmyMvHerTHBoL2ehx0pkA2xtY+ZElsfCeAE1pwBLrJy5Ym1PvEqF8GgFA6fWIyNvuEsw+7SWfuYCtZHAt1ROQNcpLjYXkBdUWF/HgxKodCJSkZFBQusM6Ghwa4+lpVJjzwiLV5sfywrk664wt9NftNZ+5gK1kcC3QUsy7KcHoQTWlpaFA6H1dzcrIKCAqeH0w3v2NGTSMSeDmpstP+4VVczlZkJnc9ETFad9GsIjkTskJrp78lrF36QavYgnBkWzhL1DopG6QsE5EJVlR28kk3bBYN2lXLDhtyNyySpPEap8PvjCP8hnPXAxHDGO3bAWelWhd54IzsVbdMr5w0Ndh+zWJ+zeE491X5zeeQIv8+AmFSzB2vODMK+iYCz0jkTMRtnGXpprVtenvTSS6yPBHqDszUNwb6JgPOc3IUhla77pnTUnzbtROf/RI4ckR580J6yZH0kkB7CmSHYNxFwnpO7MKRTOXdyjVZv30jy+wpIHdOahoi9Y08FfYGA7HBqF4ZUu+6n01E/EpEeekiaM8f+mKku/E5P/QJ+QOXMEOybCDgvtgtDqmdrZmpqLpOV80RnfC9YkJkzvp2c+gX8gsqZIdg3ETDDypX2WYSJqkPZ2IUhU5Xz2Lq1LVvsz6NRqa3txO+W2Lq1vpxYwBtJIPsIZ4Zg30TADE7swpCpwJOLM755IwlkH33O6HMGIIFcnWWYbn+1eB31c9mjjUa9QO+kmj1Yc2aQ2Dv26dPtRb8n7xBQUcEOAUAu5eosw0ysdcvlGd8rV6b2RjKTU7+AnxDODFNcTF8gwI/6GnhyuVCfN5JAdhHODEVfIMBf+hp4cr1QnzeSQPaw5sygNWcAIPUu8GRi3RqA7GLNGQC4VG8q5071aAOQebTSAACPcKJHG4DMI5wBgEc40aMNQOYxrQkAHsJCfcD9CGcA4EGc8Q24lzHTmkuWLFFRUZEGDBig8vJyrV+/Punxq1at0sUXX6xBgwZp2LBhuuWWW3T48OEcjRYAACA7jAhnq1ev1p133qn58+dr+/btuvLKKzV58mTt3bs37vEbNmxQTU2NvvGNb+jPf/6zfvGLX+gPf/iDbr311hyPHAAAILOM6HM2ceJEXXrppVq6dGnHdePGjdOUKVO0cOHCbsf/8Ic/1NKlS7Vnz56O6xYtWqTvf//72rdvX0r/J33OAABwQOMb0p5lklKJHwFp9K3SkDSa+BnMNX3Ojh07pm3btunee+/tcv2kSZO0adOmuF9TWVmp+fPna926dZo8ebIOHjyoF198UZ///OcT/j9Hjx7V0aNHOz5vaWnJzB0AAACpe3+9tHuRFAgq+QReu2RFpYKxnglnqXJ8WvPQoUOKRqMaOnRol+uHDh2qAwcOxP2ayspKrVq1SlOnTlW/fv109tlnq7CwUIsWLUr4/yxcuFDhcLjjMmLEiIzeDwAAkIKiGilUaAcvqy3JJWofV1Tj9IhzzvFwFhMIBLp8bllWt+ti6uvrNXv2bD3wwAPatm2bXn75Zb311luaOXNmwu8/b948NTc3d1xSnf4EAAAZFDpVuvAeSfH/xp8QkC66VwoNzsWojOL4tOYZZ5yhYDDYrUp28ODBbtW0mIULF6qqqkp33323JKmsrEynnHKKrrzySj366KMaNmxYt6/p37+/+vfvn/k7AAAA0jNmllT/uNTWlPiYUFi6YFbOhmQSxytn/fr1U3l5uerq6rpcX1dXp8rKyrhf89FHHykvr+vQgx+3wDbg/AYAAJBMj9Uz/1bNJAPCmSTNnTtXy5Yt0/Lly7Vr1y7NmTNHe/fu7ZimnDdvnmpqTsw5X3/99aqtrdXSpUv15ptvauPGjZo9e7Yuu+wyDR8+3Km7AQAAUjVmll0di8fHVTPJgGlNSZo6daoOHz6shx9+WPv379f48eO1bt06jRo1SpK0f//+Lj3Ppk+frtbWVv3Hf/yHvv3tb6uwsFCf/vSn9fjjjzt1FwAAQDpi1bOd96lrWw1/V80kQ/qcOYE+ZwAAOKytVVo7suvas1ChNGWfJ8NZqtnDiGlNAADgQ93WnlE1kwyZ1gQAAD4T2ymg/aiUF5Laj9kfj7wpbb39pIO9tVNATwhnAAAg9zrvFBBbYdV+XNrz3EkHZmGnAMO3kCKcAQCA3CuqkXbef1Kvs3bJau9+bKZ3CjB8CynWnAEAgNxzcqcAw7eQonIGAJ1EIlJtrdTUJBUWStXVUmmp06MCPMqpnQIStvE4mTMnKFA5AwBJDQ1SVZVUViY98oi0eLH9saxMuuIK+3YAGebkTgHJmuDGONQMl3AGwPcaGqSJE6UtW+zPo1Gprc3+KEmbN9u3E9CALHBqpwCDt5AinAHwvWnTpObmE2HsZNGoffv06TkdFuAPCUNSDsKRoVtIEc4A+FokIm3alDiYxUSj0saN9vEAMixeSMpFOHIyGCZBOAPga7W1UjCY2rHBoLRmTXbHA/iSkzsFOBUMkyCcAfC1piYpL8XfhHl5UmNjVocD+FfnkJTLcGTgFlKEMwC+VlgotcfpeRlPe7s0ZEhWhwP4V0dIUu7DkVPBMAH6nAHo4MceX9XV0oIFqR0bjdrHA8iSktlS/mBp9Izc/r8dfc/mOV41k6SAZVmpbCzlOS0tLQqHw2publZBQYHTwwEc1dBgn7G4aZO9riovz64SRaN2768VK6TiYqdHmT1VVXYbjWQnBQSDUkWFtGFD7sYFIIeOfyTtWW4Hw/xBWfkvUs0ehDPCGXwu1uMrUSuJYFAKh+3w4tWAxmMAIBdSzR6sOQN8jh5fduDassWujEl2GAuFTpzFWVFBMAOQO6w5A3ws1uOrJ517fHl1DVpxsT1lGYnY7TIaG+3F/9XV0vjxTo8OgJ8QzgAfi/X46qkBq3Six5dXw1lMaan37yMAszGtCfgYPb4AwDyEM8DH6PEFAOYhnAE+Vl2d2pSmRI8vAMgVwhngY6WlUmVlz3tLBoN2LzAWxgNA9hHOAJ9budLu4ZUooMV6fK1YkdNhAYBvEc4An6PHFwCYhVYaAOjxBcBdGt+Q9iyTlMomRwFp9K3SkLJsjypjCGcAOtDjC4ArvL9e2r1ICgSVfBKwXbKiUsFYV4UzpjUBAIC7FNVIoUI7eFltSS5R+7iiGqdHnBbCGQAAcJfQqdKF90gK9HBgQLroXik0OBejyhjCGQAAcJ8xs6RQOPkxobB0wazcjCeDCGcAAMB9eqyeubNqJhHOAACAWyWrnrm0aiYRzgAAgFslrJ65t2om0UoDSFskItXWSk1N9sbh1dW0nwAAx4yZJdU/LrU1nbjOxVUziXAGpKyhQZo2Tdq0ye6en5cntbdLCxbY+06uWEEXfQDIuVj1bOd9spvSurtqJjGtCaSkoUGaONHexkiSolGprc3+KEmbN9u3NzQ4N0YA8K3Oa89cXjWTCGdASqZNk5qbT4Sxk0Wj9u3Tp+d0WAAAqdPaM7m+aiYxrQn0KBKxpzJ7Eo1KGzfax7MGDQByrGS2lD9YGj3D6ZH0GeEMfeb1BfK1tfYas0RVs86CQXvjcC/dfwBwhfxBUsltTo8iIwhn6DW/LJBvarLvWyrhLC9PamzM+pAAAB7GmjP0ip8WyBcW2qEzFe3t0pAhWR0OAMDjCGfoFT8tkK+uTq1qJtnHVVdndzwAAG8jnCFtsQXyPQWWzgvk3ay0VKqstKdukwkG7enc8eNzMy4AgDcRzpC22AL5VMQWyLvdypVSOJz4fgeD9u0rVuR0WAAADyKcIW2xBfKp8MoC+eJie31dRYX9eTAohUInwlpFhX27F06AAAA4i7M1kTa/LpAvLpY2bLCnadessUPnkCH2GjOmMgEAmRKwLMtyehBOaGlpUTgcVnNzswoKCpwejqtEIlJZWXrHE14AAH6XavZgWhNpY4E8AADZQzhDr7BAHgCA7CCcoVdYIA8AQHZwQgB6jQXyAABkHuEMfVZaykbfAABkCtOaAAAABiGcAQAAGIRwBgAAYBDCGQAAgEEIZwAAAAYhnAEAABiEcAYAAGAQwhkAAIBBCGcAAAAGIZwBAAAYhHAGAABgEMIZAACAQQhnAAAABiGcAQAAGIRwBgAAYBDCGQAAgEEIZwAAAAbJd3oAAAAggcY3pD3LJFkpHByQRt8qDSnL9qiQZYQzAABS4URQen+9tHuRFAgq+WRXu2RFpYKxhDMPIJwBAJAKJ4JSUY20836prUlSNPmxoUL7eLgea84AAEhFUY0dgKyoZLUluUQzF5RCp0oX3iMp0MOBAemie6XQ4L7/n3Ac4QwAgFQ4FZTGzJJC4R7GFpYumJWZ/w+OY1oTMFwkItXWSk1NUmGhVF0tlZY6PSrAp8bMkuof/3iaMYFMB6VYKNx5n+Kvd6Nq5jVUzgBDNTRIVVVSWZn0yCPS4sX2x7Iy6Yor7NsB5FiP1bMsBaVk1TOqZp5DOAMM1NAgTZwobdlifx6NSm1t9kdJ2rzZvp2ABjjAiaCUMBRSNfMiwhlgoGnTpObmE2HsZNGoffv06TkdFgDJuaAULxRSNfMkwhlgmEhE2rQpcTCLiUaljRvt4wHkmBNBqVsopGrmVYQzwDC1tVIwmNqxwaC0Zk12xwMgDqeCUudQSNXMswhngGGamqS8FH8y8/KkxsasDgdAIk4EpY5QKKpmHkY4AwxTWCi1t6d2bHu7NGRIVocDIBGnglLJbKl8kTTm9tz8f8g5whlgmOrqntebxUSj9vEAHOJEUMofJJXcZn+EJxHOAMOUlkqVlT2vOwsG7T5o48fnZlwA4iAoIQsIZ4CBVq6UwuHEAS0YtG9fsSKnwwIA5ADhDDBQcbHdgLaiwv48GJRCoRNhraLCvr242LkxAgCyg701AUMVF0sbNth9zNassc/KHDLEXmPGVCZ8r/ENac8yxd9r8mQBafSt0pCybI8KyAjCGWC40lI2Oge6eX+9tHuRFAgq+SRQu2RFpYKxhDO4BtOaAAD3KaqRQoV28LLaklyi9nFFNU6PGEgZ4QwA4D4J97c8GVscwX2MCWdLlixRUVGRBgwYoPLycq1fvz7p8UePHtX8+fM1atQo9e/fX6NHj9by5ctzNFoAgOPi7W95MrY4ggsZseZs9erVuvPOO7VkyRJVVVXpmWee0eTJk1VfX6+RI0fG/Zobb7xRf//73/Xss8+quLhYBw8e1PHjx3M8cgCAY2LVs533Kf6JAVTN4E4By7JSOdUlqyZOnKhLL71US5cu7bhu3LhxmjJlihYuXNjt+Jdfflk33XST3nzzTZ122mkp/R9Hjx7V0aNHOz5vaWnRiBEj1NzcrIKCgr7fCQBA7rW1SmtHSm1N3W8LFUpT9hHOYIyWlhaFw+Ees4fj05rHjh3Ttm3bNGnSpC7XT5o0SZs2bYr7Nb/85S81YcIEff/739c555yjMWPG6K677tI//vGPhP/PwoULFQ6HOy4jRozI6P0AADgg4dozqmZwL8fD2aFDhxSNRjV06NAu1w8dOlQHDhyI+zVvvvmmNmzYoD/96U9as2aNnnjiCb344ouaNSvxuoJ58+apubm547Jv376M3g8AgEPirT1jrRlczIg1Z5IUCHR912NZVrfrYtrb2xUIBLRq1SqFw/YP5I9//GN95Stf0eLFizVw4MBuX9O/f3/1798/8wMHADir29ozqmZwN8crZ2eccYaCwWC3KtnBgwe7VdNihg0bpnPOOacjmEn2GjXLsvTuu+9mdbwAAAN1rp5RNYPLOR7O+vXrp/LyctXV1XW5vq6uTpWVlXG/pqqqSu+9956OHDnScd3u3buVl5enc889N6vjBQAYqGPtmaiawfUcD2eSNHfuXC1btkzLly/Xrl27NGfOHO3du1czZ86UZK8Xq6k50d355ptv1umnn65bbrlF9fX1eu2113T33XdrxowZcac0AQA+UDJbKl8kjbnd6ZEAfWLEmrOpU6fq8OHDevjhh7V//36NHz9e69at06hRoyRJ+/fv1969ezuOHzx4sOrq6nT77bdrwoQJOv3003XjjTfq0UcfdeouAACclj9IKrnN6VEAfWZEnzMnpNprBAAAIBNc0+cMAAAAJxDOAAAADEI4AwAAMIgRJwQAMF8kItXWSk1NUmGhVF0tlZY6PSoA8B7CGYCkGhqkadOkTZukYFDKy5Pa26UFC6SqKmnFCqm42OlRAoB3EM4AB5lejWpokCZOlJqb7c+jUfsSs3mzffuWLQQ0AMgUWmnQSgMOSFSNikbNqkZVVdnBq3MgO1kwKFVUSBs25G5cAOBGqWYPKmc+ZnrVxqvcUo2KROzw2JNoVNq40T6e1w8A9B3hzIdYQ+SsadPsYJaoGhWN2rdPn+5sNaq21n59JKuaxQSD0po1hDMAyATCmc+4pWrjVW6qRjU12cE9lXCWlyc1NmZ9SEDvNL4h7VkmKZVVPAFp9K3SkLJsjwpIiHDmM26p2niVm6pRhYV2RTUV7e3SkCFZHQ7Qe++vl3YvkgJBJW/v2S5ZUalgLOEMjiKc+YibqjZe5aZqVHW1PdWdimjUPh7ImExWu4pqpJ33S21Nknr44QsV2scDDiKc+YibqjZe5aZqVGmpVFmZ+tma48fnbmzwgUxWu0KnShfeI+28T8nDXkC66F4pNLj34wYygO2bfCRWtUmF01Ubr6quTi0cS2ZUo1aulMJhO4DFEwzat69YkdNhwQ+KauwqlhWVrLYkl2hq1a4xs6RQOPkxobB0waxM3QOg1whnPuKmqo1XxapRicJOTDBonznrdDWquNiunFVUnBhXKHRi/BUVnDyCLIlVuxTo4cAUq109fj+qZjAHTWh91IQ2EpHK0ljjGok4Hw68qPMZs/GqaLFqlGmhJxKxp7obG+3gXl3N6wNZ1tYqrR358VqxBEKF0pR9qYWqZN8vne8D9BJNaNENa4jMEKtGTZ9un3hx8g4BFRVm9porLWUNInKsx7ViaVa7En4/qmYwC5UzH1XOJPdWbbyKahTQg0xXu+J9P6pmyBEqZ4jLrVUbr6IaBfQg09Wubt+PqhnMQ+XMZ5WzzqjaAHCFTFe7On8/qmbIISpn6BFVGwCukOlqV8f3m0fVDEailQYAwHyd+5Rloh9ZyWypfJE05va+jw3IMMIZAMB8HX3KlJlqV/4gqeQ2+yNgGKY1AQDuUDJbyh8sjZ7h9EiArCKcAQDcIVbtAjyOaU0AAACDUDmDkSIRqbbW3qy9sNBu88GZpQAAPyCcwSgNDdK0adKmTV0b5C5YYG8EToNcAIDXMa0JY8S2ltqyxf48GpXa2k5sM7V5s317Q4NzYwQAINsIZzDGtGmJ9/yU7Oubm+2tpwAA8CqmNWGESMSeyuxJNGrvCRqJsAYNgAEa35D2LFPXfT8TCUijb5WGlGV7VHA5whmMUFtrrzFLVDXrLBi09wQlnAFw3Pvrpd2LpEBQySej2iUrKhWMJZyhR4QzGKGpyV78n0o4y8uzN2t3AmeRAuiiqEbaef/Hm7L38AssVGgf7yZUBh1BOIMRCgvtszJT0d4uDRmS1eF0w1mkAOLqtil7In3crN0pVAYd0esTAv75z38qEonoo48+6nbbxo0b+zQo+E91dWpVM8k+rro6u+PpjLNIASTVeVP2RDKxWbsTimrsip8Vlay2JJeoOyuDhupVOHv99dc1YsQIfepTn9KZZ56pxx57rMvtkydPzsjg4B+lpVJlpV2VSiYYtCtV48fnZlwSZ5EC6EHHpuyBBAe4tGompXDfYlx8Hw3Uq3D27W9/Wz/60Y90+PBhbdu2TbW1tZoxY4baP56XsqxU5qaBrlaulMLhxAEtGLRvX7Eid2OKnUXaU1Wv81mkAHwoWfXMrVWzGC9XBg3Vq3BWX1+vmhq7dDl27Fi9+uqrOnjwoL7yla/o2LFjGR0g/KO42J46rKiwPw8GpVDoRFirqLBvz+XarthZpKmInUUKwIcSVpg8UFHycmXQUL0KZwUFBfrb3/7W8fnAgQO1du1aDRgwQJ/73Oc6KmhAuoqLpQ0bpDfekB54QJo1y/4YidjX53rRfews0lQ4eRYpAAPEqzB5paLk5cqggXoVzq655ho999xzXa7Lz8/XqlWrNHr0aP3jH//IyODgX6Wldij7f//P/pjLNWadmX4WKQCDdKsweaii5OXKoIECVi8WiB07dkzHjx/XoEGD4t6+d+9ejRw5ss+Dy6aWlhaFw2E1NzeroKDA6eHAUJGIVJbGWeGRiHNBEoAB2lqltSPtvmehQmnKPu8El873LcZr9zHLUs0eaVXO/v3f/12S1K9fv4TBTJLxwQxIlclnkQIwUEeFSd6rKHm5MmiYtMLZ448/rjlz5iS8fe/evX0eEGAaE88iBWCwktlS+SJpzO1OjyTzOq89Y61Z1qS1Q0Btba1uvPFGHTlyRD/5yU8UCNjpubW1Vd/97nf11FNPxW1KC7hZ7CzS6dPtdhmddwiIRu2zSNkhAECH/EFSyW2pHeu27ZE6dkSYR9Usi9IKZ1/4whe0bt06ffGLX9SHH36o5557TsuXL9eCBQv0wQcfaMaMGdkaJ+Co2FmkkYjdLqOx0V78X13NVCaAPnDj9kgls6X8wdJo/uZnS69OCNi6dauuvvpq5eXl6ciRI7rhhhv02GOPqaSkJBtjzApOCAAAOC7eIvtEWHzvelk5IUCStm/frvvuu08ffvihWltbVVVVpRdffNFVwQwAACOwPRLiSCuc3XzzzfrkJz+p+vp6LV++XOvXr9ef/vQnTZkyRUePHs3WGAEA8C62R8JJ0gpn//M//6MHH3xQf/3rXzV9+nRVVVXpN7/5jf7whz/oc5/7nI4cOZKtcQIA4E1sj4STpLXm7MCBAzr77LO7Xf+Xv/xF11xzjYYNG6bf//73GR1gtrDmDABgjGRrz1hr5hlZWXMWL5hJUklJidavX69GNhYEACB9bI+ETnq1t2Y85513ntavX5+pbwcAgL94eeN0pCVj4UxKXFkDAAA9YHskfCyj4QwAAPQB2yNBhDMAAMzh5Y3TkbK0tm9CbkUiUm2t1NQkFRbaWwWVljo9KgBAVrE9ku/1avsmLzC5lUZDgzRtmrRpU/dNtquq2GQbAAA3ytr2TciuhgZp4kRpyxb782hUamuzP0rS5s327Q0Nzo0RAABkD9Oahpk2TWpuPhHGThaN2rdPny5t2JDToQEATNP4hrRnmaRUJsEC0uhbpSFl2R4V+ohwZpBIxJ7K7Ek0Km3caB/PGjQA8LH310u7F0mBoJJPhrVLVlQqGEs4cwGmNQ1SW2uvMUtFMCitWZPd8QAADFdUY2/vZEUlqy3JJWofV1Tj9IiRAipnBmlqshf/J5rS7CwvT2K3LADIIROnEGOtN3be18O4aGjrJoQzgxQW2mdlpqK9XRoyJKvDMR6tRgDklKlTiGNmSfWPx980PSYvJB15U9p6e5JvxJo0U9BKw6BWGpGIVJbGz0QkIo0fn73xmIpWIwAc0dYqrR2ZPATFhAqlKftyV6n682M9VM8CUiBZPebjQDlhsTTmW1kYICRaabhSaalUWdnzurNg0A4hfg1mtBoB4Ihue18m4sAUYrxN07uwWJPmIoQzw6xcKYXDiQNaMGjfvmJFTodljHRajQBAxvUYguTMnpgJg2NAGva5ONefjDVpJiGcGaa42K4KVVTYnweDUih0IqxVVNi3+3HaLtZqpKcTJmKtRr7+demhh+yvA4CM6LF65mDIiRccQ2GpYrmZgRIJEc4MVFxsN5h94w3pgQekWbPsj5GIfb0fg5mUXqsRSXr+eemRR+x1fFdcwVQngAxJVj1zMuR0C44fB8WBw8wNlIiLszUNVlrK2YedpdNqROp65mtsLZpfq44AMihh+woDQk7nMzc7B8VkZ3RSNTMOlTO4RjqtRk7GWjQAGZVoCtHpkNNRPVPXoJhsTZrTgRLdEM7gGtXVqVfN4um87RUA9EmiKUQTQk7JbKl8kTTmpJ5mpgZKdEM4g2uk2mokGba9ApAxncOOSSEnf5BUcpv9sTOTAyW6IJzBVXpqNdITtr0CkDGJphBNZmqgRBeEM7hKvFYjeWm8itn2CkBGJZpCNJUbA6UPcbYmXCfWaiQSsaco//pX6Wc/S+1ro1F77RoAZERsCtFNSmZL+YOl0TOcHgkSIJzBtTq3GnnzTbuiluyEgWDQrrj5cdsrAOjgxkDpM0xrwhPY9goA4BWEM3gC214BALyCaU14xslr0Rob7cX/1dVMZQIA3INwBs9h2ysAgJsxrQkAAGAQKmcAAKC7xjekPcvUdXP3RALS6FulIWXZHpUvEM4AAEB376+Xdi+SAkEln2hrl6yoVDCWcJYhhDMAAHLNDVWpohpp5/1SW5OkJE0kJSlUaB+PjCCcAQCQa26oSsW2etp5n5KHSDZQzzROCAAAINeKauxqkxWVrLYkl6izVanOG6UnwgbqGUc4AwAg1zo2IA/0cKDDVakex0nVLBsClmWlMuHtOS0tLQqHw2publZBQYHTw3GOG9Y9AIAXtbVKa0d+vKYrgVChNGWfs+En2ThNGJ+LpJo9WHPmd25Y9wAAXtTjmq5eVqUy/aY74TipmmWLMZWzJUuW6Ac/+IH279+viy66SE888YSuvPLKHr9u48aNuuqqqzR+/Hjt2LEj5f+PytnHUnnnFsM7JADIrLZWac050vHW7rfl9ZPOn25/7JBCmNq9WNp6W+pvuicslsZ8q+dxnvy3gr8JaXNV5Wz16tW68847tWTJElVVVemZZ57R5MmTVV9fr5EjRyb8uubmZtXU1Ogzn/mM/v73v+dwxB7C2TgA/MqEZR2hU6WhV0t/+2X329qPS3ue63xFajMY2WiB0e1vBX8TssmIytnEiRN16aWXaunSpR3XjRs3TlOmTNHChQsTft1NN92kCy64QMFgUGvXrqVy1ltuWfcAAJmUjQpTb3z0nrT2nNSOTfV38Z8fS+1N9yULP17wn4LOfyv4m9ArqWYPx8/WPHbsmLZt26ZJkyZ1uX7SpEnatGlTwq977rnntGfPHj344IMp/T9Hjx5VS0tLlws+xtk4APzIlHYWg4ZLZ1+bwoFp/C7ORguMjr8V4m9Cljkezg4dOqRoNKqhQ4d2uX7o0KE6cOBA3K/561//qnvvvVerVq1Sfn5qM7MLFy5UOBzuuIwYMaLPY/eUZD/I9LAB4EUmtbOoWN7zMen8Ls7Wm+6S2VL5ImnM7el9HdLieDiLCQS6voAsy+p2nSRFo1HdfPPNeuihhzRmzJiUv/+8efPU3Nzccdm3b1+fx+wpCX+QqZoB8DBTmqz2WD3rxe/ibLzpzh8kldxmf0TWOB7OzjjjDAWDwW5VsoMHD3arpklSa2urtm7dqttuu035+fnKz8/Xww8/rJ07dyo/P1+/+c1v4v4//fv3V0FBQZcLThLvB5mqGQAvM2lZx+UrpbyB8W/rze9i3nS7luPhrF+/fiovL1ddXV2X6+vq6lRZWdnt+IKCAkUiEe3YsaPjMnPmTJWUlGjHjh2aOHFirobuPd1+kPkBBuADpizrGDhUKn1AGQ1TvOl2JcfDmSTNnTtXy5Yt0/Lly7Vr1y7NmTNHe/fu1cyZMyXZU5I1NfZCzLy8PI0fP77L5ayzztKAAQM0fvx4nXLKKU7eFffr/IPMDzAAPzCpwpTpMMWbblcyIpxNnTpVTzzxhB5++GFdcskleu2117Ru3TqNGjVKkrR//37t3bvX4VH6BGfjAPAjUypM2QhTvOl2HSP6nDmBPmdJHP9I2rNcGj2DRZ8A/KNLb7A0e4BlUjb6if35MWnnPOmSx5y5T5Dkoj5nMBBn4wDwI1MqTNmYwaAFhqsYsX0TAACO69iiaJ7zyzpKZkv5g+0ZjEyIvemGKxDOAACIyXQo6i3ClK8RzgAAiCEUwQCsOQMAADAI4QwAAMAghDMAAACDEM4AAAAMQjgDAAAwCOEMAADAIIQzAAAAgxDOAAAADEI4AwAAMAjhDAAAwCBs3wQAQKoa35D2LJNkpXBwQBp9qzSkLNujgscQzgAASNX766Xdi6RAUMknn9olKyoVjCWcIW1MawIAkKqiGilUaAcvqy3JJWofV1Tj9IjhQoQzAABSFTpVuvAeSYEeDgxIF90rhQbnYlTwGMIZAADpGDNLCoWTHxMKSxfMys144DmsOUPaIhGptlZqapIKC6Xqaqm01OlRAUCOxKpnO+9T/BMDqJqhbwKWZaVyyonntLS0KBwOq7m5WQUFBU4PxxUaGqRp06RNm6RgUMrLk9rbpWhUqqqSVqyQioudHiUA5EBbq7R2pNTW1P22UKE0ZR/hDN2kmj2Y1kRKGhqkiROlLVvsz6NRqa3N/ihJmzfbtzc0ODdGAMiZhGvPqJqh7whnSMm0aVJz84kwdrJo1L59+vScDgsAnBNv7RlrzZABhDP0KBKxpzITBbOYaFTauNE+HgA8r1v1jKoZMoNwhh7V1tprzFIRDEpr1mR3PABgjM7VM6pmyBDCGXrU1GQv/k9FXp7U2JjV4QCAOTqqZ6JqhoyhlQZ6VFhon5WZivZ2aciQrA4HAMxSMlvKHyyNnuH0SOARVM7Qo+rqntebxUSj9vEA4Bv5g6SS2+yPQAZQOUOPSkulykq7jUaykBYMShUV0vjxff8/aXQLAPArmtDShDYlsT5nidppBINSOGwHuL40oqXRLQDAq2hCi4wqLraDV0WF/XkwKIVCJ87irKjITDCj0S0AwO+Y1kTKioulDRvsKcc1a+yzMocMsaccMzGVmU6j2w0b+v7/AZDU+Ia0Z5ni7xF5soA0+lZpSFm2RwX4GuEMaSstzfz6r1ij2550bnTLGjQgA95fL+1eJAWCSj6Z0i5ZUalgLOEMyDKmNWEEGt0CDimqsTfqtqKS1ZbkErWPK6pxesSA5xHOYAQa3QIOSbiB98nYmgjIFcIZjECjW8BB8TbwPhlbEwE5QziDEWh0Cziox+oZVTMglwhnMEKs0W1P686CQbvfWSbODgXQSbLqGVUzIKcIZzDGypV2I9tEAS3W6HbFipwOC/CHhNUzqmZArhHOYIxcNLoFkES86hlVMyDn6HMGo2S70S2AJGLVs533yW5KS9UMcAJ7a/p0b002FgcQV1urtHak1NZk9zWbso9wBmRIqtmDypnPJNpYfMECNhYHoE7Vs3lUzQCHEM58JLaxeHOz/Xk02rV9RWxjcdZ1AT5XMlvKHyyNnuH0SABfIpz5CBuLw3PYtDs78gdJJbc5PQrAtwhnPsHG4vAkNu3OLcIwkBOEM5+IbSyeShf+2MbihDMYr6hG2nm/vXhdPby42bS77wjDQE7Q58wn2FgcnsSm3blVVGOHXCsqWW1JLlHCMNAHhDOfYGNxeBabducOYRjICcKZT7CxODyLTbtzizAMZB3hzCfYWByexqbduUMYBrKOcOYjbCwOz2LT7twiDANZRTjzETYWh6exaXfuEIaBrKKVhs+wsTj6zNReV2zanVtjZkn1j3/cxuRjhGEgIwhnPlVaSh8z9JLJva46BwaCQnYRhoGsYVoTQHpM7nXVMd0mgkIudJ5KJgwDGUM4A5Ae03tdlcyWyhdJY27P7f/rR4RhICsClmWlsnDEc1paWhQOh9Xc3KyCggKnhwO4S1urtHZk1/VGJwsVSlP28Qfb645/JO1ZLo2eYW+YDiChVLMHlTMA6aPXFWLyB0kltxHMgAwinAHoHXpdAUBWcLYmgN7pdrZeDFUzx5na7gRASghnAHqPXldmMrndCYAeMa0JoPe6rT2jamYEk9udAOgR4QxA39DryjymtzsBkBThDEDf0OvKTMlO2IghTANGIpwB6Dsav5qHdieAaxHOAPQdva7MRLsTwJUIZwDgVQmrZ1TNAJMRzgDAy+JVz6iaAUYjnAGAl9HuBHAdwhkAeB3tTgBXIZwBgNfR7gRwFbZvAgA/KJkt5Q+WRs9weiTZY8qeoqaMA65FOAMAP4i1O/EyU/YUNWUccC2mNQEA3mDKnqKmjAOuRTgDAHiDKXuKmjIOuBbhDADgHabsKWrKOOBKhDMAgHeYsqeoKeOAKxHOAADeYsqeoqaMA65DOAMAeIspe4qaMg64DuEMAOA9puwpaso44CqEMwCA95iyp6gp44CrEM4AAN5kyp6ipowDrsEOAQAAb4pVrXbOc7ZaZco4kmHLKaMQzgAA3mXKnqKmjCMRtpwySsCyrFRisue0tLQoHA6rublZBQUFTg8HAADntLVKa0dKbU09HxsqlKbsM7MCaLhUswdrzgAA8Du2nDIK4QwAALDllEEIZwAAgC2nDGLMCQFLlizRD37wA+3fv18XXXSRnnjiCV155ZVxj62trdXSpUu1Y8cOHT16VBdddJEWLFiga6+9NsejBuAanI0G9GzMLKn+8fhrz6ia5YwR4Wz16tW68847tWTJElVVVemZZ57R5MmTVV9fr5EjR3Y7/rXXXtNnP/tZfe9731NhYaGee+45XX/99dqyZYs+8YlPOHAPABiPs9GAnnW0/bhPXd/IUDXLJSPO1pw4caIuvfRSLV26tOO6cePGacqUKVq4cGFK3+Oiiy7S1KlT9cADD6R0PGdrAj7D2WhAauL9rPAzkRGuOVvz2LFj2rZtmyZNmtTl+kmTJmnTpk0pfY/29na1trbqtNNOS3jM0aNH1dLS0uUCwEc4Gw1IDVtOOc7xcHbo0CFFo1ENHTq0y/VDhw7VgQMHUvoeP/rRj/Thhx/qxhtvTHjMwoULFQ6HOy4jRozo07gBuBBnowGpYcspRzkezmICga7vZi3L6nZdPC+88IIWLFig1atX66yzzkp43Lx589Tc3Nxx2bdvX5/HDMBlOBsNSE3Hz4r4mXCA4ycEnHHGGQoGg92qZAcPHuxWTTvZ6tWr9Y1vfEO/+MUvdM011yQ9tn///urfv3+fxwvA5TgbDUiN6VtOeZjj4axfv34qLy9XXV2dvvSlL3VcX1dXpy9+8YsJv+6FF17QjBkz9MILL+jzn/98LoYKwFTptskoqrHP3ORsNCCx/EFSyW1Oj8KXHA9nkjR37lx9/etf14QJE3T55ZfrJz/5ifbu3auZM2dKsqck//a3v+k///M/JdnBrKamRk8++aQqKio6qm4DBw5UONzDehIA3pNum4xP/MiuknU5G42qGQAzGLHmbOrUqXriiSf08MMP65JLLtFrr72mdevWadSoUZKk/fv3a+/evR3HP/PMMzp+/LhmzZqlYcOGdVzuuOMOp+4CACcV1din+ltRyWpLconaxxX/f5yNBsBYRvQ5cwJ9zgCP+fNjcRpnniwgXbLQDmadeznRwwlADrimzxkAZES6bTI4Gw2AoQhnALyhN20ySmZL5YukMbfnYoQAkBLCGQDvSFY9i7fgP3Y2Wv6g7I8NAFJEOAPgHQmrZyz4B+AehDMA3hKvekabDAAuQjgD4C1s2gzA5QhnALyHTZsBuJgROwQA+Fi62xCNvlUaUpbtUblPrHq2cx5VMwCuQzgDTJLuNkQFYwlniTixaTPhund43IAuCGeASYpqpJ33f7znYzT5saFC+3jE58SmzYTr3uFxA7pgzRlgkh4bqcawyN1I6e7xSbi28bgBXRDOANOkuw0RzEG47h0eN6ALwhlgmt5sQwRzEK57h8cN6EA4A0yU7jZEMAfhund43IAOhDPARGxD5G6E697hcQMkEc4Ac7ENkXsRrnuHxw2QRDgDzMU2RO5GuO4dHjeAcAYYjW2I3Itw3Ts8bgDhDDBaxx8q8QfKjQjXvcPjBp8jnAGmK5ktlS+Sxtzu9EiQLsJ17/C4wefYvgkwnRPbECFznNjj0wt43OBjhDMAyCbCde/wuMHHmNYEAAAwCOEMAADAIIQzAAAAgxDOAAAADEI4AwAAMAjhDAAAwCCEMwAAAIPQ5wwA4B2Nb0h7lkmyUjg4II2+VRpSlu1RAWkhnAEAvOP99dLuRVIgqOSTQ+2SFZUKxhLOYBzCGYDMoWoBpxXVSDvvl9qaJEWTHxsqtI8HDEM4A5A5VC36joDbN7FN03fep+SPYYBN1WEswhmAzKFq0XcE3L4bM0uqf/zj12ECobB0waycDQlIB2drAsicWNVCgR4OpGqRUFGNHVytqGS1JblECbiJ9Pg65PUHsxHOAGTWmFl2VSIZqhaJEXAzI9nrkNcfDEc4A5BZVC36joDbdwlfh7z+YD7CGYDMo2rRNwTczIj3OuT1BxcgnAHIPKoWfUfA7btur0Nef3AHwhmA7KBq0TcE3Mzo/Drk9QeXIJwByA6qFn1HwO27jteheP3BNehzBiB7OvebIlSkr1tDVQJur5TMlvIHS6NnOD0Sf6CRcp8RzgBkT0e4mEeo6C0Cbt/lD5JKbnN6FP5BI+U+Y1oTQHaVzJbKF0ljbnd6JO7EtBzchkbKfUblDEB2UbXoO6bl/MfNU4Psb9pnhDMAMB0B13/cPjXI/qZ9wrQmAACmcfvUII2U+4RwBgCAabywxyqNlHuNcAYAgIncvscqjZR7jXAGAICJvDA1SCPlXiGcAQBgKrdPDbJTSK8QzgAAMJUXpgbZ3zRthDMAAEzm9qlBGimnjXAGAIDJvDA1yE4haSGcAQBgOrdPDcYaKecPcnokrkA4AwDAdEwN+grbNwEA4AbsseobhDMAANyAPVZ9g2lNAAAAgxDOAAAADEI4AwAAMAjhDAAAwCCEMwAAAIMQzgAAAAxCOAMAADAI4QwAAMAgNKEFAKAvGt+Q9iyTZKVwcEAafas0pCzbo4KLEc4AAOiL99dLuxdJgaCST0i1S1ZUKhhLOENSTGsCANAXRTVSqNAOXlZbkkvUPq6oxukRw3CEMwAA+iJ0qnThPZICPRwYkC66VwoNzsWo4GJMawIA0FdjZkn1j0ttTYmPCYWlC2blbEi9xho6xxHOAADoq1j1bOd9ih9qXFQ1Yw2d45jWBAAgE8bMsqtj8bilaiaxhs4AhDMAADIh4dozF1XNJNbQGYBwBgBApsSrnrmpahaTrAoY48b75RKEMwAAMqVb1cml1aUeq2cuvV8uQTgDACCTOled3Fxd8soaOhcinAEAkEkdVSe5u7rklTV0LkQ4AwAg00pmS+WLpDG3Oz2SvvHKGjqXoc8ZAACZlj9IKrnN6VH0Xbf+bQ5XzXzSIJdwBsB7fPILHMiJzrsfOF0180mDXMIZAO/xyS9wICc6qmfznF9rVlQj7bz/422yosmPdXGDXNacAfAeOpwDmWXKGjqfNMglnAHwHp/8AgdyJraGLn+Q0yPxRYNcwhkAb/LBL3DAl3zQIJdwBsCbfPALHPAtjzfIJZwB8C6P/wIHfMvjDXIJZwC8y+O/wAFf83CDXMIZAG/z8C9wwNe8ssl8HIQzAN7m4V/ggO95ZZP5kxgTzpYsWaKioiINGDBA5eXlWr9+fdLjX331VZWXl2vAgAE6//zz9fTTT+dopABcx6O/wAHf88om8ycxIpytXr1ad955p+bPn6/t27fryiuv1OTJk7V37964x7/11lu67rrrdOWVV2r79u267777NHv2bL300ks5HjkAV/DoL3AAMqdBbgYFLMtKZfO5rJo4caIuvfRSLV26tOO6cePGacqUKVq4cGG34++55x798pe/1K5duzqumzlzpnbu3KnXX3897v9x9OhRHT16tOPzlpYWjRgxQs3NzSooKMjgvQFgpOMfSXuWS6NnmNFIE4DvtLS0KBwO95g9HK+cHTt2TNu2bdOkSZO6XD9p0iRt2rQp7te8/vrr3Y6/9tprtXXrVrW1tcX9moULFyocDndcRowYkZk7AMAdTOpwDgBJOB7ODh06pGg0qqFDh3a5fujQoTpw4EDcrzlw4EDc448fP65Dhw7F/Zp58+apubm547Jv377M3AEAAIAMynd6ADGBQNc+RJZldbuup+PjXR/Tv39/9e/fv4+jBAAAyC7HK2dnnHGGgsFgtyrZwYMHu1XHYs4+++y4x+fn5+v000/P2lgBAACyzfFw1q9fP5WXl6uurq7L9XV1daqsrIz7NZdffnm341955RVNmDBBoVAoa2MFAADINsfDmSTNnTtXy5Yt0/Lly7Vr1y7NmTNHe/fu1cyZMyXZ68Vqamo6jp85c6beeecdzZ07V7t27dLy5cv17LPP6q677nLqLgAAAGSEEWvOpk6dqsOHD+vhhx/W/v37NX78eK1bt06jRo2SJO3fv79Lz7OioiKtW7dOc+bM0eLFizV8+HA99dRT+vKXv+zUXQAAAMgII/qcOSHVXiMAAACZ4Jo+ZwAAADiBcAYAAGAQwhkAAIBBCGcAAAAGIZwBAAAYhHAGAABgEMIZAACAQQhnAAAABjFihwAnxHrvtrS0ODwSAADgB7HM0VP/f9+Gs9bWVknSiBEjHB4JAADwk9bWVoXD4YS3+3b7pvb2dr333ns69dRTFQgEMvI9W1paNGLECO3bt48toQzA82EWng/z8JyYhefDLNl4PizLUmtrq4YPH668vMQry3xbOcvLy9O5556ble9dUFDAD5ZBeD7MwvNhHp4Ts/B8mCXTz0eyilkMJwQAAAAYhHAGAABgEMJZBvXv318PPvig+vfv7/RQIJ4P0/B8mIfnxCw8H2Zx8vnw7QkBAAAAJqJyBgAAYBDCGQAAgEEIZwAAAAYhnAEAABiEcAYAAGAQwlmalixZoqKiIg0YMEDl5eVav3590uNfffVVlZeXa8CAATr//PP19NNP52ik/pDO81FbW6vPfvazOvPMM1VQUKDLL79c//u//5vD0Xpfuj8fMRs3blR+fr4uueSS7A7QZ9J9Po4ePar58+dr1KhR6t+/v0aPHq3ly5fnaLT+kO5zsmrVKl188cUaNGiQhg0bpltuuUWHDx/O0Wi967XXXtP111+v4cOHKxAIaO3atT1+TU7/nltI2c9//nMrFApZP/3pT636+nrrjjvusE455RTrnXfeiXv8m2++aQ0aNMi64447rPr6euunP/2pFQqFrBdffDHHI/emdJ+PO+64w3r88cet3//+99bu3butefPmWaFQyPrjH/+Y45F7U7rPR0xTU5N1/vnnW5MmTbIuvvji3AzWB3rzfNxwww3WxIkTrbq6Ouutt96ytmzZYm3cuDGHo/a2dJ+T9evXW3l5edaTTz5pvfnmm9b69eutiy66yJoyZUqOR+4969ats+bPn2+99NJLliRrzZo1SY/P9d9zwlkaLrvsMmvmzJldrhs7dqx17733xj3+O9/5jjV27Ngu1/3bv/2bVVFRkbUx+km6z0c8F154ofXQQw9lemi+1NvnY+rUqdb9999vPfjgg4SzDEr3+fj1r39thcNh6/Dhw7kYni+l+5z84Ac/sM4///wu1z311FPWueeem7Ux+lEq4SzXf8+Z1kzRsWPHtG3bNk2aNKnL9ZMmTdKmTZvifs3rr7/e7fhrr71WW7duVVtbW9bG6ge9eT5O1t7ertbWVp122mnZGKKv9Pb5eO6557Rnzx49+OCD2R6ir/Tm+fjlL3+pCRMm6Pvf/77OOeccjRkzRnfddZf+8Y9/5GLInteb56SyslLvvvuu1q1bJ8uy9Pe//10vvviiPv/5z+diyOgk13/P8zP+HT3q0KFDikajGjp0aJfrhw4dqgMHDsT9mgMHDsQ9/vjx4zp06JCGDRuWtfF6XW+ej5P96Ec/0ocffqgbb7wxG0P0ld48H3/961917733av369crP51dRJvXm+XjzzTe1YcMGDRgwQGvWrNGhQ4f0rW99Sx988AHrzjKgN89JZWWlVq1apalTp+qf//ynjh8/rhtuuEGLFi3KxZDRSa7/nlM5S1MgEOjyuWVZ3a7r6fh416N30n0+Yl544QUtWLBAq1ev1llnnZWt4flOqs9HNBrVzTffrIceekhjxozJ1fB8J52fj/b2dgUCAa1atUqXXXaZrrvuOv34xz/WihUrqJ5lUDrPSX19vWbPnq0HHnhA27Zt08svv6y33npLM2fOzMVQcZJc/j3n7WqKzjjjDAWDwW7vcA4ePNgtTcecffbZcY/Pz8/X6aefnrWx+kFvno+Y1atX6xvf+IZ+8Ytf6JprrsnmMH0j3eejtbVVW7du1fbt23XbbbdJssOBZVnKz8/XK6+8ok9/+tM5GbsX9ebnY9iwYTrnnHMUDoc7rhs3bpwsy9K7776rCy64IKtj9rrePCcLFy5UVVWV7r77bklSWVmZTjnlFF155ZV69NFHmX3JoVz/PadylqJ+/fqpvLxcdXV1Xa6vq6tTZWVl3K+5/PLLux3/yiuvaMKECQqFQlkbqx/05vmQ7IrZ9OnT9fzzz7NuI4PSfT4KCgoUiUS0Y8eOjsvMmTNVUlKiHTt2aOLEibkauif15uejqqpK7733no4cOdJx3e7du5WXl6dzzz03q+P1g948Jx999JHy8rr+mQ4Gg5JOVG2QGzn/e56V0ww8KnYa9LPPPmvV19dbd955p3XKKadYb7/9tmVZlnXvvfdaX//61zuOj516O2fOHKu+vt569tlnaaWRQek+H88//7yVn59vLV682Nq/f3/Hpampyam74CnpPh8n42zNzEr3+WhtbbXOPfdc6ytf+Yr15z//2Xr11VetCy64wLr11luduguek+5z8txzz1n5+fnWkiVLrD179lgbNmywJkyYYF122WVO3QXPaG1ttbZv325t377dkmT9+Mc/trZv397R1sTpv+eEszQtXrzYGjVqlNWvXz/r0ksvtV599dWO26ZNm2ZdddVVXY7/3e9+Z33iE5+w+vXrZ5133nnW0qVLczxib0vn+bjqqqssSd0u06ZNy/3APSrdn4/OCGeZl+7zsWvXLuuaa66xBg4caJ177rnW3LlzrY8++ijHo/a2dJ+Tp556yrrwwgutgQMHWsOGDbP+5V/+xXr33XdzPGrv+e1vf5v074HTf88DlkVtFAAAwBSsOQMAADAI4QwAAMAghDMAAACDEM4AAAAMQjgDAAAwCOEMAADAIIQzAAAAgxDOAAAADEI4AwAAMAjhDADi2L9/vwYPHqybbrqpy/W/+tWvFAqFNH/+fIdGBsDrCGcAEMewYcP0ne98R//1X/+lbdu2SZJ+97vf6atf/aq++c1v6rvf/a7DIwTgVYQzAEjgrrvu0rBhw3TPPffoD3/4g2644QZ97Wtf05NPPtlxzNKlS3XppZcqFAppwYIFzg0WgGcQzgAggUGDBunRRx/V//3f/+nqq6/W5MmT9dOf/lSBQKDjmGHDhumhhx7SlClTnBsoAE/Jd3oAAGCyMWPGSJICgYBWrFihYDDY5fZYKPvv//7vXA8NgEdROQOABHbs2KEvfOELqqqq0pEjR7R8+XKnhwTABwhnABDHX/7yF1177bW6/PLL9dvf/lZf/OIXtWDBAjU3Nzs9NAAeRzgDgJO8/fbbuuaaa1RSUqKXXnpJoVBIjz32mBobG/W9733P6eEB8DjCGQB0sn//fl1zzTU666yz9Ktf/UoDBw6UJI0dO1YzZszQk08+qbffftvZQQLwtIBlWZbTgwAAtzp+/LiOHz+ub37zmzrnnHN0//33KxQKdTtxAABSRTgDgD5YsGCBHnrooS7XPffcc5o+fbozAwLgeoQzAAAAg7DmDAAAwCCEMwAAAIMQzgAAAAxCOAMAADAI4QwAAMAghDMAAACDEM4AAAAMQjgDAAAwCOEMAADAIIQzAAAAgxDOAAAADPL/A/VmqD/F92+/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (7,7))\n",
    "\n",
    "plt.scatter(X[y == 0,0],\n",
    "            X[y == 0,1],\n",
    "            c = \"blue\",\n",
    "            s = 60,\n",
    "            label=\"0\")\n",
    "plt.scatter(X[y == 1,0],\n",
    "            X[y == 1,1],\n",
    "            c = \"orange\",\n",
    "            s = 60,\n",
    "            marker = \"v\",\n",
    "            label=\"1\")\n",
    "\n",
    "plt.legend(fontsize = 12)\n",
    "plt.xlabel(\"$x_1$\",fontsize = 12)\n",
    "plt.ylabel(\"$x_2$\",fontsize = 12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y,\n",
    "                                    shuffle=True,\n",
    "                                    random_state=435,\n",
    "                                    test_size=.3,\n",
    "                                    stratify=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 0.9841269841269841\n",
      "Test F1 score is 0.888888888888889\n"
     ]
    }
   ],
   "source": [
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty good for our small data set!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iris Data Set\n",
    "\n",
    "Now we can apply this procedure to a more interesting data set: the Iris data set, which is\n",
    "described in detail at: https://archive.ics.uci.edu/dataset/53/iris\n",
    "\n",
    "Each observation represents an iris (a type of flower) and gives it's measurements including:\n",
    "- `sepal_length`: the length of the iris's sepal in cm.\n",
    "- `sepal_width`: the width of the iris's sepal in cm.\n",
    "- `petal_length`: the length of the iris's petal in cm.\n",
    "- `petal_width`: the width of the iris's petal in cm.\n",
    "- `iris_class`: the class of the iris, can be:\n",
    "    - `0` meaning it is a setosa iris\n",
    "    - `1` meaning it is a versicolor iris\n",
    "    - `2` meaning it is a virginica iris\n",
    "\n",
    "You can take a look at some of the EDA for this data set done in the Single_Layer_Perceptron.ipynb\n",
    "notebook in this repository, so we will just load this and split up our data right away here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "## to get the iris data\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our model is a binary classifier,\n",
    "let's just see if we can separate the setosa iris (class 0) observations from the other two.\n",
    "\n",
    "To do so, we need to reassign any observation with iris_class = 2 to iris_class = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Changing all the iris_class = 2 to class 1 so that we can use our binary classifier\n",
    "iris_train.loc[iris_train.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 2, 'iris_class'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0,\n",
       "       1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0,\n",
       "       1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0,\n",
       "       0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1,\n",
       "       1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1,\n",
       "       1, 1])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1,\n",
       "       1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like our iris classes look good. Now we can train and test our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 1.0\n",
      "Test F1 score is 1.0\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like this does fairly well in separating iris class 0 from the other two. Let's also\n",
    "compare the other two and see whether we get similar discriminating power. We will need to reload\n",
    "the data set a couple times to do this of course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 0 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0\n",
      " 1]\n"
     ]
    }
   ],
   "source": [
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']\n",
    "\n",
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])\n",
    "\n",
    "# Changing all the iris_class = 1 to class 0 so that we can use our binary classifier to hopefully\n",
    "# separate class 2 from the others\n",
    "iris_train.loc[iris_train.iris_class == 1, 'iris_class'] = 0\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 1, 'iris_class'] = 0\n",
    "\n",
    "\n",
    "# We also need to re-assign the value of class 2 to 1 so that our sigmoid activation function works\n",
    "# properly \n",
    "iris_train.loc[iris_train.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "\n",
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 0.9599999999999999\n",
      "Test F1 score is 0.962962962962963\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like it did pretty well! Moving on now to see if our model separates class 1 from the other two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 1 1 0 0 1 1 0 0 0 1 0 0 1 0 0 0 1 0 1 1 1 0 1 0 0 1 0 0 0 0 1\n",
      " 0]\n"
     ]
    }
   ],
   "source": [
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']\n",
    "\n",
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])\n",
    "\n",
    "# Changing all the iris_class = 1 to class 0 so that we can use our binary classifier to hopefully\n",
    "# separate class 2 from the others\n",
    "iris_train.loc[iris_train.iris_class == 2, 'iris_class'] = 0\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 2, 'iris_class'] = 0\n",
    "\n",
    "# We can leave class 1 alone here\n",
    "\n",
    "\n",
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 0.4918032786885246\n",
      "Test F1 score is 0.5454545454545455\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0]\n",
      "[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0]\n"
     ]
    }
   ],
   "source": [
    "print(list(y_test))\n",
    "print(y_test_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly this does not perform as well in attempting to separate class 1 from the other two. As a\n",
    "final exploration, we will compare class 1 to each of the other two individually and see if there is\n",
    "anything we can learn from that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 1 0 1 1 1 0 0 0 1 0 0 1 0 1 1 0 1 0 0 1 1 1 0]\n"
     ]
    }
   ],
   "source": [
    "#### Comparing classes 1 and 0\n",
    "\n",
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']\n",
    "\n",
    "# Remove all class 2 observations\n",
    "iris_df = iris_df.loc[iris_df.iris_class < 2]\n",
    "\n",
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])\n",
    "\n",
    "\n",
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 1.0\n",
      "Test F1 score is 1.0\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like we get great separation between classes 0 and 1, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 1 0 1 1 1 0 0 0 1 0 0 1 0 1 1 0 1 0 0 1 1 1 0]\n"
     ]
    }
   ],
   "source": [
    "#### Comparing classes 1 and 2\n",
    "\n",
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']\n",
    "\n",
    "# Remove all class 2 observations\n",
    "iris_df = iris_df.loc[iris_df.iris_class > 0]\n",
    "\n",
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])\n",
    "\n",
    "# Changing all the iris_class = 1 to class 0 so that we can use our binary classifier\n",
    "iris_train.loc[iris_train.iris_class == 1, 'iris_class'] = 0\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 1, 'iris_class'] = 0\n",
    "\n",
    "\n",
    "# We also need to re-assign the value of class 2 to 1 so that our sigmoid activation function works\n",
    "# properly \n",
    "iris_train.loc[iris_train.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "\n",
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 0.9866666666666666\n",
      "Test F1 score is 0.962962962962963\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 1 0 1 1 1 0 0 0 1 0 0 1 0 1 1 0 1 0 0 1 1 1 0]\n"
     ]
    }
   ],
   "source": [
    "#### Comparing classes 0 and 2\n",
    "\n",
    "## Load the data\n",
    "iris = load_iris()\n",
    "iris_df = pd.DataFrame(iris['data'], columns = ['sepal_length','sepal_width','petal_length','petal_width'])\n",
    "iris_df['iris_class'] = iris['target']\n",
    "\n",
    "# Remove all class 2 observations\n",
    "iris_df = iris_df.loc[iris_df.iris_class != 1]\n",
    "\n",
    "## Making a statified train test split\n",
    "iris_train, iris_test = train_test_split(iris_df, \n",
    "                                            random_state=431,\n",
    "                                            shuffle=True,\n",
    "                                            test_size=.25,\n",
    "                                            stratify=iris_df['iris_class'])\n",
    "\n",
    "\n",
    "# We also need to re-assign the value of class 2 to 1 so that our sigmoid activation function works\n",
    "# properly \n",
    "iris_train.loc[iris_train.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "# also be sure to make the same changes to the test set\n",
    "iris_test.loc[iris_test.iris_class == 2, 'iris_class'] = 1\n",
    "\n",
    "\n",
    "# Define our input and output for the training data\n",
    "X_train = iris_train[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "X_test = iris_test[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']].to_numpy()\n",
    "\n",
    "y_train = iris_train.iris_class.to_numpy()\n",
    "y_test = iris_test.iris_class.to_numpy()\n",
    "print(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training F1 score is 1.0\n",
      "Test F1 score is 1.0\n"
     ]
    }
   ],
   "source": [
    "# Load and fit our model:\n",
    "\n",
    "standardize(X_train)\n",
    "standardize(X_test)\n",
    "\n",
    "# Initialize a class instance\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# Fit our model to the training data\n",
    "model = logreg.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = logreg.predict(X_train)\n",
    "\n",
    "f1_score_train = F1_score(y_train, y_train_pred)\n",
    "print(\"Training F1 score is %s\" % f1_score_train)\n",
    "\n",
    "# get prediction on test data\n",
    "y_test_pred = logreg.predict(X_test)\n",
    "\n",
    "f1_score_test = F1_score(y_test, y_test_pred)\n",
    "print(\"Test F1 score is %s\" % f1_score_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Takeaway here is that our model is very good at distinguishing class 0 from the other classes, but\n",
    "is not quite as accurate distinguishing between classes 1 and 2. This is all good though; in a\n",
    "future notebook we will build from scratch a multi-class logistic model to better handle this kind\n",
    "of data set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}