LinearRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff7e3cb4-c101-4d1b-ae26-28b7b9c0d6e6",
   "metadata": {},
   "source": [
    "## Multivariate Linear Regression from The Ground Up\n",
    "## with Minimal Math And Maximum Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3d18cbb6-87c0-41d5-896e-21e8e816ce17",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3615467b-b50d-4cbc-8ad2-a099ae682934",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "N = 10\n",
    "d = 6\n",
    "X = np.round(100*np.random.randn(N, d))/100\n",
    "kgt = np.round(100*np.random.randn(d))/2 # The answer I want to seek"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fcd23e77-bfb5-48ae-9dc8-64ab986d58bc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.76  0.4   0.98  2.24  1.87 -0.98]\n",
      " [ 0.95 -0.15 -0.1   0.41  0.14  1.45]\n",
      " [ 0.76  0.12  0.44  0.33  1.49 -0.21]\n",
      " [ 0.31 -0.85 -2.55  0.65  0.86 -0.74]\n",
      " [ 2.27 -1.45  0.05 -0.19  1.53  1.47]\n",
      " [ 0.15  0.38 -0.89 -1.98 -0.35  0.16]\n",
      " [ 1.23  1.2  -0.39 -0.3  -1.05 -1.42]\n",
      " [-1.71  1.95 -0.51 -0.44 -1.25  0.78]\n",
      " [-1.61 -0.21 -0.9   0.39 -0.51 -1.18]\n",
      " [-0.03  0.43  0.07  0.3  -0.63 -0.36]]\n"
     ]
    }
   ],
   "source": [
    "# X[row]\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c8694800-08d9-4219-ab95-c266a381ef84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-33.5 -18.  -40.5 -86.5   9.  -20. ]\n"
     ]
    }
   ],
   "source": [
    "print(kgt)\n",
    "assert(len(kgt) == len(X[0])) # This throw an error if they're not equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6003d29b-a85d-4223-9f87-10f35fd3af15",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rows = X.shape[0]\n",
    "n_cols = X.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b64c0ec-8686-4961-90f8-67e60c883a5a",
   "metadata": {},
   "source": [
    "### Step 1: Let's Manually write out a loop to do a dot product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "08db1a2f-3048-4df1-bbfb-498c2b10b3fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "j = 0\n",
    "print(j)\n",
    "j = j+1 # Overwrite j with j+1\n",
    "print(j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4ea8f7ea-1419-4cca-84d3-2370eba500a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-25.46\n",
      "-2.16\n",
      "-17.82\n",
      "-28.545\n",
      "13.41\n",
      "4.2\n",
      "-56.375\n"
     ]
    }
   ],
   "source": [
    "i = 2 # This is a row that I want to take the dot product of with kgt\n",
    "dot = -33.5*0.76 + (-18)*0.12 + (-40.5)*0.44\n",
    "dot = kgt[0]*X[i][0] + kgt[1]*X[i][1] + kgt[2]*X[i][2] # +....\n",
    "\n",
    "def get_dot_of_row(X, kgt, i):\n",
    "    \"\"\"\n",
    "    Compute and return the dot product of row i of X with kgt\n",
    "    \"\"\"\n",
    "    dot = 0 ## Accumulator variable: add up a bunch of stuff\n",
    "    j = 0 ## Column index\n",
    "    n_cols = X.shape[1] ## How many times our loop has to go through\n",
    "    while j < n_cols:\n",
    "        # \"Accumulator pattern\"\n",
    "        # I'm going to add kgt[j]*X[i][j] onto what's already in dot\n",
    "        dot = dot + kgt[j]*X[i][j] \n",
    "        print(kgt[j]*X[i][j])\n",
    "        j = j+1\n",
    "    return dot\n",
    "\n",
    "print(get_dot_of_row(X, kgt, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12747a78-d77f-4c42-8ee1-3f4ecda613a1",
   "metadata": {},
   "source": [
    "### Step 2: Use numpy to do this faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d5e453f5-7cad-4ddb-818f-f86e3f5fc637",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-25.46   -2.16  -17.82  -28.545  13.41    4.2  ]\n",
      "-56.375\n"
     ]
    }
   ],
   "source": [
    "print(X[2]*kgt)\n",
    "\n",
    "print(np.sum(X[2]*kgt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f5fdab2-494e-4416-bfa9-2975091d691e",
   "metadata": {},
   "source": [
    "### Step 3: Do all dot products at once with matrix multiplication\n",
    "\n",
    "\n",
    "Suppose I had $N$ data points, each with dimension 10\n",
    "\n",
    "Store each data point in a row of a matrix, and each column stores a dimension\n",
    "\n",
    "Now suppose I have a single observation for each data point that can be described by a linear model\n",
    "\n",
    "So my observation model is $k . x$, where $k$ is our hidden model vector in 10 dimensions\n",
    "\n",
    "$ y \\approx k . x$\n",
    "\n",
    "Written as a matrix multiplication for all $y$, stacking up each observation in $Y$, which is a $N x 1$ matrix, so the matrix multiplication:\n",
    "\n",
    "$ Y = X k $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dd445dca-5471-49a2-82ba-34381d8b349e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-263.18   -88.28   -56.375   74.505  -51.165  189.1     -2.11    54.05\n",
      "   79.44   -33.99 ]\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(X, kgt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9569468e-7c67-4ab8-97c9-d587eee09194",
   "metadata": {},
   "source": [
    "but in practice, there's some noise $\\epsilon$, which I'll put in a parallel matrix $\\epsilon$.  So a better observation model is\n",
    "\n",
    "$Y = Xk + \\epsilon$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4ed7d50e-b6ea-47ec-80c6-81ce74bc1245",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y.shape (10,)\n"
     ]
    }
   ],
   "source": [
    "eps = np.random.randn(X.shape[0])\n",
    "Y = np.dot(X, kgt) + eps ## What we actually observe: it's corrupted by noise!\n",
    "print(\"Y.shape\", Y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bede2173-f576-4322-9bde-5357a1a0983b",
   "metadata": {},
   "source": [
    "Now we want to solve for $k$ given $X$ data and noise-corrupted observation $Y$, we want to solve for $k$.  Hopefully, this will agree with the $kgt$ that we came up with.\n",
    "\n",
    "If we want to solve this as an optimization problem, we need to minimize some objective function.  This objective function should measure the difference between our observation $Y$ and what our model tells us, which is $X k$, over the variable we're trying to solve for, which is $k$\n",
    "\n",
    "One simple objective function is the L2 distance between $Y$ and $Xk$.  That's a fancy way of saying the sum of the squared differences\n",
    "\n",
    "$f(Y, Xk) = \\sum_{i=0}^{d-1} (Y_i - (Xk)_i)^2$\n",
    "\n",
    "We want to minimize this\n",
    "\n",
    "\n",
    "To minimize an objective function with gradient descent\n",
    "\n",
    "* Randomly initialize your variables\n",
    "\n",
    "* for some number of steps:\n",
    "    1) Compute the gradient of the objective function with respect to the variables\n",
    "    2) Take a step against the gradient to update your variables\n",
    "\n",
    "In the code below, I'm going to approximately compute the gradient by varying each indpendent component of $k$ by itself and seeing how the objective function changes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0863b7eb-4a6f-4136-970a-0f449f6cb0bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC60lEQVR4nO3dfXxU5Z3///dMZjKEEIaEmAyjEdCyAgarxi6CVugKgZrA+uhv600w6rc+Yq0CpsKKftv9Su1KqLpsW1lF+9utbrXG7/4Q11o2C1KKSwk3DaQSkGoVCYSEcJNMEiB3M9fvjzAHhiAgTOZMktfzsfNI5pzPnHOdE9Z597quc47DGGMEAADQDzntbgAAAIBdCEIAAKDfIggBAIB+iyAEAAD6LYIQAADotwhCAACg3yIIAQCAfstldwPiXSgU0v79+5WSkiKHw2F3cwAAwHkwxqi5uVl+v19O5xf3+xCEzmH//v3KysqyuxkAAOAC7N27V5dddtkXricInUNKSoqkrhM5ePBgm1sDAADOR1NTk7Kysqzv8S9CEDqH8HDY4MGDCUIAAPQy55rWwmRpAADQbxGEAABAv0UQAgAA/RZBCAAA9FsEIQAA0G8RhAAAQL9FEAIAAP0WQQgAAPRbBCEAANBvEYQAAEC/RRACAAD9FkEIAAD0Wzx01SZNrR1qOt6h5ESXUpMT7W4OAAD9Ej1CNvmXtX/RzT9Zq39Z+xe7mwIAQL9FELKJy+mQJHWGjM0tAQCg/yII2cTl7Dr1naGQzS0BAKD/IgjZxJ3Q1SMUpEcIAADbEIRsknCiR6gjSBACAMAuBCGb0CMEAID9CEI2STgxWbojyBwhAADsQhCyiSuh69TTIwQAgH0IQjZxWT1CBCEAAOxCELJJOAgFuXweAADbEIRs4krghooAANiNIGQT64aKDI0BAGAbgpBNTj5ig6ExAADsQhCySfiqMYbGAACwD0HIJlaPEENjAADYhiBkEyZLAwBgP4KQTRKsHiHmCAEAYBeCkE3c3FkaAADbEYRsYj1rjKvGAACwDUHIJu4T9xEKMlkaAADbEIRscrJHiCAEAIBdCEI2cScwWRoAALsRhGxiXTVGjxAAALYhCNkkfNUYN1QEAMA+BCGbhHuEuHweAAD7EIRsEr6zNJfPAwBgH4KQTVwnLp83RgrRKwQAgC0IQjYJ9whJ9AoBAGCXLx2EPvjgA82YMUN+v18Oh0PvvPOOta6jo0MLFizQuHHjlJycLL/fr3vvvVf79++P2EZbW5vmzJmj9PR0JScna+bMmdq3b19ETUNDgwoLC+X1euX1elVYWKjGxsaImurqas2YMUPJyclKT0/X3Llz1d7eHlGzfft2TZo0SUlJSbr00kv19NNPyxj7e2DCT5+XmCcEAIBdvnQQOnr0qL761a9q6dKl3dYdO3ZMW7du1T/8wz9o69atevvtt/Xxxx9r5syZEXXFxcVasWKFSktLtX79erW0tCg/P1/BYNCqKSgoUGVlpcrKylRWVqbKykoVFhZa64PBoPLy8nT06FGtX79epaWlWr58uebNm2fVNDU1aerUqfL7/dqyZYteeOEFPf/881qyZMmXPeyoCw+NSVIHV44BAGAPcxEkmRUrVpy1ZvPmzUaS2bNnjzHGmMbGRuN2u01paalVU1NTY5xOpykrKzPGGLNz504jyWzcuNGqKS8vN5LMrl27jDHGrFy50jidTlNTU2PVvPnmm8bj8ZhAIGCMMebFF180Xq/XtLa2WjUlJSXG7/ebUCh0xva2traaQCBgvfbu3WskWduMlmAwZIYveM8MX/CeOdzSFtVtAwDQ3wUCgfP6/u7xOUKBQEAOh0NDhgyRJFVUVKijo0O5ublWjd/vV3Z2tjZs2CBJKi8vl9fr1fjx462aG2+8UV6vN6ImOztbfr/fqpk2bZra2tpUUVFh1UyaNEkejyeiZv/+/fr888/P2N6SkhJrOM7r9SorKysq5+F0TqdD4dEx7i4NAIA9ejQItba26oknnlBBQYEGDx4sSaqrq1NiYqJSU1MjajMzM1VXV2fVZGRkdNteRkZGRE1mZmbE+tTUVCUmJp61Jvw+XHO6J598UoFAwHrt3bv3yx72eXOFb6rIHCEAAGzh6qkNd3R06K677lIoFNKLL754znpjjByOkxOIT/09mjXmxETpM31WkjweT0QPUk9yOR1qF3eXBgDALj3SI9TR0aE77rhDu3fv1urVq63eIEny+Xxqb29XQ0NDxGfq6+ut3hqfz6cDBw502+7Bgwcjak7v1WloaFBHR8dZa+rr6yWpW0+RHVzW88YYGgMAwA5RD0LhEPTJJ5/o/fff19ChQyPW5+TkyO12a/Xq1day2tpaVVVVaeLEiZKkCRMmKBAIaPPmzVbNpk2bFAgEImqqqqpUW1tr1axatUoej0c5OTlWzQcffBBxSf2qVavk9/s1YsSIaB/6l8bQGAAA9vrSQailpUWVlZWqrKyUJO3evVuVlZWqrq5WZ2en/u7v/k5//OMf9cYbbygYDKqurk51dXVWGPF6vXrggQc0b948rVmzRtu2bdM999yjcePGacqUKZKkMWPGaPr06SoqKtLGjRu1ceNGFRUVKT8/X1dddZUkKTc3V2PHjlVhYaG2bdumNWvWaP78+SoqKrJ6oAoKCuTxeHT//ferqqpKK1as0KJFi/TYY4994dBYLFk9QgyNAQBgjy97OdratWuNpG6v++67z+zevfuM6ySZtWvXWts4fvy4mT17tklLSzNJSUkmPz/fVFdXR+zn8OHDZtasWSYlJcWkpKSYWbNmmYaGhoiaPXv2mLy8PJOUlGTS0tLM7NmzIy6VN8aYDz/80Hz96183Ho/H+Hw+s3Dhwi+8dP5MzvfyuwsxYdH7ZviC98yf9jZEfdsAAPRn5/v97TAmDm6zHMeamprk9XoVCAQi5jpFwy3PrlX1kWN6++GJuv7y1HN/AAAAnJfz/f7mWWM2YmgMAAB7EYRsFH7wKleNAQBgD4KQjRJOPG+MHiEAAOxBELKR+0SPEE+fBwDAHgQhGyWcmCPUwbPGAACwBUHIRu4TQ2P0CAEAYA+CkI2sHiGCEAAAtiAI2chlzRFiaAwAADsQhGzksuYI0SMEAIAdCEI2sh66ShACAMAWBCEbhXuEGBoDAMAeBCEbhXuEGBoDAMAeBCEbnewRIggBAGAHgpCNwkGonRsqAgBgC4KQjdwuJksDAGAngpCN3E6ePg8AgJ0IQjZyn5gszdAYAAD2IAjZiKExAADsRRCykZunzwMAYCuCkI3c3EcIAABbEYRsdPKGivQIAQBgB4KQjdwJDI0BAGAngpCN3Dx0FQAAWxGEbMTl8wAA2IsgZKPw0FgnQQgAAFsQhGzEVWMAANiLIGQjN1eNAQBgK4KQjVxcNQYAgK0IQjZKZGgMAABbEYRsRI8QAAD2IgjZiDlCAADYiyBkI+uGiiGGxgAAsANByEbWIzY66RECAMAOBCEbWUNj9AgBAGALgpCNeOgqAAD2IgjZyOoRYmgMAABbEIRs5GJoDAAAWxGEbMTQGAAA9iII2Sh8Z2ljpCC9QgAAxNyXDkIffPCBZsyYIb/fL4fDoXfeeSdivTFGCxculN/vV1JSkiZPnqwdO3ZE1LS1tWnOnDlKT09XcnKyZs6cqX379kXUNDQ0qLCwUF6vV16vV4WFhWpsbIyoqa6u1owZM5ScnKz09HTNnTtX7e3tETXbt2/XpEmTlJSUpEsvvVRPP/20jImP0BEeGpPoFQIAwA5fOggdPXpUX/3qV7V06dIzrn/22We1ZMkSLV26VFu2bJHP59PUqVPV3Nxs1RQXF2vFihUqLS3V+vXr1dLSovz8fAWDQaumoKBAlZWVKisrU1lZmSorK1VYWGitDwaDysvL09GjR7V+/XqVlpZq+fLlmjdvnlXT1NSkqVOnyu/3a8uWLXrhhRf0/PPPa8mSJV/2sHtEeGhMIggBAGALcxEkmRUrVljvQ6GQ8fl8ZvHixday1tZW4/V6zbJly4wxxjQ2Nhq3221KS0utmpqaGuN0Ok1ZWZkxxpidO3caSWbjxo1WTXl5uZFkdu3aZYwxZuXKlcbpdJqamhqr5s033zQej8cEAgFjjDEvvvii8Xq9prW11aopKSkxfr/fhEKh8zrGQCBgJFnbjKZgMGSGL3jPDF/wnjnc0hb17QMA0F+d7/d3VOcI7d69W3V1dcrNzbWWeTweTZo0SRs2bJAkVVRUqKOjI6LG7/crOzvbqikvL5fX69X48eOtmhtvvFFerzeiJjs7W36/36qZNm2a2traVFFRYdVMmjRJHo8nomb//v36/PPPz3gMbW1tampqinj1FKfToQQnE6YBALBLVINQXV2dJCkzMzNieWZmprWurq5OiYmJSk1NPWtNRkZGt+1nZGRE1Jy+n9TUVCUmJp61Jvw+XHO6kpISa16S1+tVVlbWuQ/8IrgIQgAA2KZHrhpzOBwR740x3Zad7vSaM9VHo8acmCj9Re158sknFQgErNfevXvP2u6LlWg9gT4+JnADANCfRDUI+Xw+Sd17W+rr662eGJ/Pp/b2djU0NJy15sCBA922f/DgwYia0/fT0NCgjo6Os9bU19dL6t5rFebxeDR48OCIV09yu048gZ4eIQAAYi6qQWjkyJHy+XxavXq1tay9vV3r1q3TxIkTJUk5OTlyu90RNbW1taqqqrJqJkyYoEAgoM2bN1s1mzZtUiAQiKipqqpSbW2tVbNq1Sp5PB7l5ORYNR988EHEJfWrVq2S3+/XiBEjonnoFyw8NNZOEAIAIOa+dBBqaWlRZWWlKisrJXVNkK6srFR1dbUcDoeKi4u1aNEirVixQlVVVbr//vs1cOBAFRQUSJK8Xq8eeOABzZs3T2vWrNG2bdt0zz33aNy4cZoyZYokacyYMZo+fbqKioq0ceNGbdy4UUVFRcrPz9dVV10lScrNzdXYsWNVWFiobdu2ac2aNZo/f76KioqsXpyCggJ5PB7df//9qqqq0ooVK7Ro0SI99thj5xyqi5Xw88Y6GRoDACD2vuzlaGvXrjWSur3uu+8+Y0zXJfRPPfWU8fl8xuPxmFtuucVs3749YhvHjx83s2fPNmlpaSYpKcnk5+eb6urqiJrDhw+bWbNmmZSUFJOSkmJmzZplGhoaImr27Nlj8vLyTFJSkklLSzOzZ8+OuFTeGGM+/PBD8/Wvf914PB7j8/nMwoULz/vSeWN69vJ5Y4yZ9OzvzPAF75ktuw/3yPYBAOiPzvf722FMnNxmOU41NTXJ6/UqEAj0yHyhqUvW6ZP6Fv26aLwmXpke9e0DANAfne/3N88as5mLoTEAAGxDELJZIk+gBwDANgQhm7m5jxAAALYhCNnMRY8QAAC2IQjZ7GSPEEEIAIBYIwjZLJEgBACAbQhCNks88YiN9k6CEAAAsUYQslk4CLURhAAAiDmCkM3CQ2M8awwAgNgjCNmMoTEAAOxDELJZ+KoxghAAALFHELKZx8VVYwAA2IUgZDOGxgAAsA9ByGZMlgYAwD4EIZtx+TwAAPYhCNmMoTEAAOxDELIZV40BAGAfgpDNrB4h5ggBABBzBCGbcfk8AAD2IQjZLJGhMQAAbEMQshmTpQEAsA9ByGZcPg8AgH0IQjbjhooAANiHIGQzN0NjAADYhiBks3CPEFeNAQAQewQhm3noEQIAwDYEIZtx1RgAAPYhCNmMO0sDAGAfgpDNTs4RMgqFjM2tAQCgfyEI2Sx81ZhErxAAALFGELJZuEdIIggBABBrBCGbnRqEOpgwDQBATBGEbOZ0OuROcEiiRwgAgFgjCMUBnkAPAIA9CEJxgHsJAQBgD4JQHHAn8AR6AADsQBCKA9xUEQAAexCE4kA4CHHVGAAAsUUQigMeV4IkhsYAAIg1glAcCD+BniAEAEBsRT0IdXZ26oc//KFGjhyppKQkXXHFFXr66acVCp38kjfGaOHChfL7/UpKStLkyZO1Y8eOiO20tbVpzpw5Sk9PV3JysmbOnKl9+/ZF1DQ0NKiwsFBer1der1eFhYVqbGyMqKmurtaMGTOUnJys9PR0zZ07V+3t7dE+7ItyMggFbW4JAAD9S9SD0E9+8hMtW7ZMS5cu1UcffaRnn31Wzz33nF544QWr5tlnn9WSJUu0dOlSbdmyRT6fT1OnTlVzc7NVU1xcrBUrVqi0tFTr169XS0uL8vPzFQyeDAsFBQWqrKxUWVmZysrKVFlZqcLCQmt9MBhUXl6ejh49qvXr16u0tFTLly/XvHnzon3YF8XjPjE01kGPEAAAMWWiLC8vz3znO9+JWPatb33L3HPPPcYYY0KhkPH5fGbx4sXW+tbWVuP1es2yZcuMMcY0NjYat9ttSktLrZqamhrjdDpNWVmZMcaYnTt3Gklm48aNVk15ebmRZHbt2mWMMWblypXG6XSampoaq+bNN980Ho/HBAKB8zqeQCBgJJ13/YUoem2LGb7gPfP6xs97bB8AAPQn5/v9HfUeoZtvvllr1qzRxx9/LEn605/+pPXr1+u2226TJO3evVt1dXXKzc21PuPxeDRp0iRt2LBBklRRUaGOjo6IGr/fr+zsbKumvLxcXq9X48ePt2puvPFGeb3eiJrs7Gz5/X6rZtq0aWpra1NFRcUZ29/W1qampqaIV0+jRwgAAHu4or3BBQsWKBAIaPTo0UpISFAwGNQzzzyju+++W5JUV1cnScrMzIz4XGZmpvbs2WPVJCYmKjU1tVtN+PN1dXXKyMjotv+MjIyImtP3k5qaqsTERKvmdCUlJfrRj370ZQ/7ogxgsjQAALaIeo/QW2+9pddff12//vWvtXXrVr322mt6/vnn9dprr0XUORyOiPfGmG7LTnd6zZnqL6TmVE8++aQCgYD12rt371nbFA0ed9efobWDydIAAMRS1HuE/v7v/15PPPGE7rrrLknSuHHjtGfPHpWUlOi+++6Tz+eT1NVbM2zYMOtz9fX1Vu+Nz+dTe3u7GhoaInqF6uvrNXHiRKvmwIED3fZ/8ODBiO1s2rQpYn1DQ4M6Ojq69RSFeTweeTyeCz38C8J9hAAAsEfUe4SOHTsmpzNyswkJCdbl8yNHjpTP59Pq1aut9e3t7Vq3bp0VcnJycuR2uyNqamtrVVVVZdVMmDBBgUBAmzdvtmo2bdqkQCAQUVNVVaXa2lqrZtWqVfJ4PMrJyYnykV+4AW4unwcAwA5R7xGaMWOGnnnmGV1++eW6+uqrtW3bNi1ZskTf+c53JHUNVRUXF2vRokUaNWqURo0apUWLFmngwIEqKCiQJHm9Xj3wwAOaN2+ehg4dqrS0NM2fP1/jxo3TlClTJEljxozR9OnTVVRUpJdfflmS9OCDDyo/P19XXXWVJCk3N1djx45VYWGhnnvuOR05ckTz589XUVGRBg8eHO1Dv2DhHqFWJksDABBTUQ9CL7zwgv7hH/5BDz/8sOrr6+X3+/Xd735X/+f//B+r5vHHH9fx48f18MMPq6GhQePHj9eqVauUkpJi1fzzP/+zXC6X7rjjDh0/fly33nqrXn31VSUkJFg1b7zxhubOnWtdXTZz5kwtXbrUWp+QkKDf/va3evjhh3XTTTcpKSlJBQUFev7556N92BeFGyoCAGAPhzHG2N2IeNbU1CSv16tAINBjvUivbfhcT727Q3njhulfZl3fI/sAAKA/Od/vb541FgfoEQIAwB4EoTjgcXMfIQAA7EAQigMDrMnS9AgBABBLBKE4QI8QAAD2IAjFAeuGilw+DwBATBGE4kD4hoqtTJYGACCmCEJxgB4hAADsQRCKA1w+DwCAPQhCcWCAm0dsAABgB4JQHDi1R4gbfQMAEDsEoTgQniMUMlJniCAEAECsEITiQPg+QhI3VQQAIJYIQnEgPDQmcVNFAABiiSAUBxwOhxWG6BECACB2CEJx4mQQokcIAIBYIQjFifAl9NxLCACA2CEIxYmkRJ5ADwBArBGE4kTSiR6h4+0MjQEAECsEoTgR7hE6To8QAAAxQxCKE+EeoWPtnTa3BACA/oMgFCeS3MwRAgAg1ghCcWJAeGisnSAEAECsEITihDVZmvsIAQAQMwShODGQydIAAMQcQShOMEcIAIDYIwjFiQFu5ggBABBrBKE4Eb6P0DGCEAAAMUMQihMMjQEAEHsEoThx8qoxghAAALFCEIoTSdxHCACAmCMIxQl6hAAAiD2CUJwI9wgxRwgAgNghCMWJAfQIAQAQcwShOHHy6fMEIQAAYoUgFCesoTGCEAAAMUMQihM8awwAgNgjCMWJ8ByhzpBRR5An0AMAEAsEoTgRniMkMU8IAIBYIQjFCXeCQy6nQxI3VQQAIFYIQnHC4XAo2eOSJB1t77S5NQAA9A89EoRqamp0zz33aOjQoRo4cKCuvfZaVVRUWOuNMVq4cKH8fr+SkpI0efJk7dixI2IbbW1tmjNnjtLT05WcnKyZM2dq3759ETUNDQ0qLCyU1+uV1+tVYWGhGhsbI2qqq6s1Y8YMJScnKz09XXPnzlV7e3tPHPZFSw4/gb6NHiEAAGIh6kGooaFBN910k9xut/7rv/5LO3fu1D/90z9pyJAhVs2zzz6rJUuWaOnSpdqyZYt8Pp+mTp2q5uZmq6a4uFgrVqxQaWmp1q9fr5aWFuXn5ysYPBkSCgoKVFlZqbKyMpWVlamyslKFhYXW+mAwqLy8PB09elTr169XaWmpli9frnnz5kX7sKNiID1CAADElomyBQsWmJtvvvkL14dCIePz+czixYutZa2trcbr9Zply5YZY4xpbGw0brfblJaWWjU1NTXG6XSasrIyY4wxO3fuNJLMxo0brZry8nIjyezatcsYY8zKlSuN0+k0NTU1Vs2bb75pPB6PCQQC53U8gUDASDrv+osx84X/McMXvGfe31nX4/sCAKAvO9/v76j3CL377ru64YYb9O1vf1sZGRm67rrr9Itf/MJav3v3btXV1Sk3N9da5vF4NGnSJG3YsEGSVFFRoY6Ojogav9+v7Oxsq6a8vFxer1fjx4+3am688UZ5vd6ImuzsbPn9fqtm2rRpamtrixiqO1VbW5uampoiXrEyMDHcI8TQGAAAsRD1IPTZZ5/ppZde0qhRo/Tf//3feuihhzR37lz9+7//uySprq5OkpSZmRnxuczMTGtdXV2dEhMTlZqaetaajIyMbvvPyMiIqDl9P6mpqUpMTLRqTldSUmLNOfJ6vcrKyvqyp+CChSdLH2tjaAwAgFiIehAKhUK6/vrrtWjRIl133XX67ne/q6KiIr300ksRdQ6HI+K9MabbstOdXnOm+gupOdWTTz6pQCBgvfbu3XvWNkVTsqdrsjQ9QgAAxEbUg9CwYcM0duzYiGVjxoxRdXW1JMnn80lStx6Z+vp6q/fG5/Opvb1dDQ0NZ605cOBAt/0fPHgwoub0/TQ0NKijo6NbT1GYx+PR4MGDI16xYg2N0SMEAEBMRD0I3XTTTfrzn/8csezjjz/W8OHDJUkjR46Uz+fT6tWrrfXt7e1at26dJk6cKEnKycmR2+2OqKmtrVVVVZVVM2HCBAUCAW3evNmq2bRpkwKBQERNVVWVamtrrZpVq1bJ4/EoJycnykd+8cKXz3PVGAAAseGK9ga///3va+LEiVq0aJHuuOMObd68Wa+88opeeeUVSV1DVcXFxVq0aJFGjRqlUaNGadGiRRo4cKAKCgokSV6vVw888IDmzZunoUOHKi0tTfPnz9e4ceM0ZcoUSV29TNOnT1dRUZFefvllSdKDDz6o/Px8XXXVVZKk3NxcjR07VoWFhXruued05MgRzZ8/X0VFRTHt6TlfA605QgyNAQAQEz1xydpvfvMbk52dbTwejxk9erR55ZVXItaHQiHz1FNPGZ/PZzwej7nlllvM9u3bI2qOHz9uZs+ebdLS0kxSUpLJz8831dXVETWHDx82s2bNMikpKSYlJcXMmjXLNDQ0RNTs2bPH5OXlmaSkJJOWlmZmz55tWltbz/tYYnn5/Mvr/mKGL3jPfP+tbT2+LwAA+rLz/f52GGOM3WEsnjU1Ncnr9SoQCPR4L9LrG/foh+9UafrVPi0rjL+hOwAAeovz/f7mWWNx5ORVY8wRAgAgFghCcYSrxgAAiC2CUBwZFJ4szX2EAACICYJQHBnI5fMAAMQUQSiOJHP5PAAAMUUQiiPhHqEW5ggBABATBKE4kuJxS5LaOkNq7wzZ3BoAAPo+glAcGTTg5I2+m1s7bGwJAAD9A0EojiQ4HdaVY82tDI8BANDTCEJxJuVEr1ATPUIAAPQ4glCcCQcheoQAAOh5BKE4M3hA14Rp5ggBANDzCEJx5uTQGD1CAAD0NIJQnEk50SPUdJweIQAAehpBKM4wRwgAgNghCMWZwUnhOUIEIQAAehpBKM6c7BFiaAwAgJ5GEIoz1hwhghAAAD2OIBRnBjNHCACAmCEIxZmT9xEiCAEA0NMIQnGGR2wAABA7BKE4k0KPEAAAMUMQijOnXjVmjLG5NQAA9G0EoTgTvo9QR9CorTNkc2sAAOjbCEJxJjkxQU5H1+88ZgMAgJ5FEIozDodDgzw8eBUAgFggCMWhkxOm6RECAKAnEYTiEM8bAwAgNghCcSh8d+kAc4QAAOhRBKE4NGRgV49Q47F2m1sCAEDfRhCKQ2nJiZKkhmP0CAEA0JMIQnEodWBXEDpylB4hAAB6EkEoDoWDUANDYwAA9CiCUBxKZWgMAICYIAjFobTkrsnSDQyNAQDQowhCcWgIQ2MAAMQEQSgOpYWDED1CAAD0KIJQHArPETraHlRbZ9Dm1gAA0HcRhOLQ4AEuJZx4BH0jE6YBAOgxBKE45HA4lHri7tLcSwgAgJ7T40GopKREDodDxcXF1jJjjBYuXCi/36+kpCRNnjxZO3bsiPhcW1ub5syZo/T0dCUnJ2vmzJnat29fRE1DQ4MKCwvl9Xrl9XpVWFioxsbGiJrq6mrNmDFDycnJSk9P19y5c9XeHv/hIpV5QgAA9LgeDUJbtmzRK6+8omuuuSZi+bPPPqslS5Zo6dKl2rJli3w+n6ZOnarm5marpri4WCtWrFBpaanWr1+vlpYW5efnKxg8OWemoKBAlZWVKisrU1lZmSorK1VYWGitDwaDysvL09GjR7V+/XqVlpZq+fLlmjdvXk8edlScvKkiQ2MAAPQY00Oam5vNqFGjzOrVq82kSZPMo48+aowxJhQKGZ/PZxYvXmzVtra2Gq/Xa5YtW2aMMaaxsdG43W5TWlpq1dTU1Bin02nKysqMMcbs3LnTSDIbN260asrLy40ks2vXLmOMMStXrjROp9PU1NRYNW+++abxeDwmEAic13EEAgEj6bzro+XBf99ihi94z/x7+ecx3S8AAH3B+X5/91iP0COPPKK8vDxNmTIlYvnu3btVV1en3Nxca5nH49GkSZO0YcMGSVJFRYU6Ojoiavx+v7Kzs62a8vJyeb1ejR8/3qq58cYb5fV6I2qys7Pl9/utmmnTpqmtrU0VFRVnbHdbW5uampoiXnYI9wg1MjQGAECPcfXERktLS7V161Zt2bKl27q6ujpJUmZmZsTyzMxM7dmzx6pJTExUampqt5rw5+vq6pSRkdFt+xkZGRE1p+8nNTVViYmJVs3pSkpK9KMf/eh8DrNHhS+hP8JNFQEA6DFR7xHau3evHn30Ub3++usaMGDAF9Y5HI6I98aYbstOd3rNmeovpOZUTz75pAKBgPXau3fvWdvUU7ipIgAAPS/qQaiiokL19fXKycmRy+WSy+XSunXr9POf/1wul8vqoTm9R6a+vt5a5/P51N7eroaGhrPWHDhwoNv+Dx48GFFz+n4aGhrU0dHRracozOPxaPDgwREvO6Sd6BE6TBACAKDHRD0I3Xrrrdq+fbsqKyut1w033KBZs2apsrJSV1xxhXw+n1avXm19pr29XevWrdPEiRMlSTk5OXK73RE1tbW1qqqqsmomTJigQCCgzZs3WzWbNm1SIBCIqKmqqlJtba1Vs2rVKnk8HuXk5ET70KMqY7BHklTf1GZzSwAA6LuiPkcoJSVF2dnZEcuSk5M1dOhQa3lxcbEWLVqkUaNGadSoUVq0aJEGDhyogoICSZLX69UDDzygefPmaejQoUpLS9P8+fM1btw4a/L1mDFjNH36dBUVFenll1+WJD344IPKz8/XVVddJUnKzc3V2LFjVVhYqOeee05HjhzR/PnzVVRUZFtPz/nKSOkaVqxvbrW5JQAA9F09Mln6XB5//HEdP35cDz/8sBoaGjR+/HitWrVKKSkpVs0///M/y+Vy6Y477tDx48d166236tVXX1VCQoJV88Ybb2ju3LnW1WUzZ87U0qVLrfUJCQn67W9/q4cfflg33XSTkpKSVFBQoOeffz52B3uBLknp6hFqONah9s6QEl3cBBwAgGhzGGOM3Y2IZ01NTfJ6vQoEAjHtRTLG6K9++F/qCBr94Ym/0aVDkmK2bwAAervz/f6mmyFOORwOXTIoPE+I4TEAAHoCQSiOXTI4PE+ICdMAAPQEglAcyzgxT4ggBABAzyAIxbFwEDpIEAIAoEcQhOJY+BL6g1xCDwBAjyAIxTFuqggAQM8iCMUx5ggBANCzCEJx7BIrCDE0BgBATyAIxbHwHKFDLe0KhrjvJQAA0UYQimPpgxKV4HQoGDJcOQYAQA8gCMUxV4JTvhM3VaxpPGZzawAA6HsIQnHu0tSuZ4ztazhuc0sAAOh7CEJx7rITD1utaSQIAQAQbQShOBfuEaqhRwgAgKgjCMW5S+kRAgCgxxCE4hw9QgAA9ByCUJw7tUfIGO4lBABANBGE4pz/RBA61h5U47EOm1sDAEDfQhCKcwPcCUof1PWoDS6hBwAgughCvYA1T4ibKgIAEFUEoV4g60QQ2nOYIAQAQDQRhHqBK9KTJUmfHz5qc0sAAOhbCEK9wMhLuoLQZwcJQgAARBNBqBe4In2QJGn3IYIQAADRRBDqBUacGBqrb25TS1unza0BAKDvIAj1At4kt9IHJUqSdjM8BgBA1BCEeomRJ3qFPjvUYnNLAADoOwhCvQTzhAAAiD6CUC8RvnKMIAQAQPQQhHqJ8L2EPj7A0BgAANFCEOolxgwbLEn6tL5FHcGQza0BAKBvIAj1EpelJmmQx6X2YIgbKwIAECUEoV7C4XBotC9FkrSrrsnm1gAA0DcQhHqR0cO6gtDOWoIQAADRQBDqRUb7uuYJ7apttrklAAD0DQShXiQ8YfojeoQAAIgKglAvctWJOUL1zW061NJmc2sAAOj9CEK9yCCPS1eeuLHin/Y22tsYAAD6AIJQL3NtVqokqZIgBADARSMI9TLXXT5EEkEIAIBoiHoQKikp0de+9jWlpKQoIyNDt99+u/785z9H1BhjtHDhQvn9fiUlJWny5MnasWNHRE1bW5vmzJmj9PR0JScna+bMmdq3b19ETUNDgwoLC+X1euX1elVYWKjGxsaImurqas2YMUPJyclKT0/X3Llz1d7eHu3Djplrs4ZIkiqrGxUKGXsbAwBALxf1ILRu3To98sgj2rhxo1avXq3Ozk7l5ubq6NGTd0N+9tlntWTJEi1dulRbtmyRz+fT1KlT1dx88rLw4uJirVixQqWlpVq/fr1aWlqUn5+vYDBo1RQUFKiyslJlZWUqKytTZWWlCgsLrfXBYFB5eXk6evSo1q9fr9LSUi1fvlzz5s2L9mHHzGhfiga4nWpu69Rnh3juGAAAF8X0sPr6eiPJrFu3zhhjTCgUMj6fzyxevNiqaW1tNV6v1yxbtswYY0xjY6Nxu92mtLTUqqmpqTFOp9OUlZUZY4zZuXOnkWQ2btxo1ZSXlxtJZteuXcYYY1auXGmcTqepqamxat58803j8XhMIBA4r/YHAgEj6bzrY+HbL20wwxe8Z97aXG13UwAAiEvn+/3d43OEAoGAJCktLU2StHv3btXV1Sk3N9eq8Xg8mjRpkjZs2CBJqqioUEdHR0SN3+9Xdna2VVNeXi6v16vx48dbNTfeeKO8Xm9ETXZ2tvx+v1Uzbdo0tbW1qaKi4oztbWtrU1NTU8Qr3uSM6JowvWn3EZtbAgBA79ajQcgYo8cee0w333yzsrOzJUl1dXWSpMzMzIjazMxMa11dXZ0SExOVmpp61pqMjIxu+8zIyIioOX0/qampSkxMtGpOV1JSYs058nq9ysrK+rKH3eMmXDFUkrTxs8MyhnlCAABcqB4NQrNnz9aHH36oN998s9s6h8MR8d4Y023Z6U6vOVP9hdSc6sknn1QgELBee/fuPWub7HDDiFS5ExyqaTyu6iPH7G4OAAC9Vo8FoTlz5ujdd9/V2rVrddlll1nLfT6fJHXrkamvr7d6b3w+n9rb29XQ0HDWmgMHDnTb78GDByNqTt9PQ0ODOjo6uvUUhXk8Hg0ePDjiFW8GJrp03Yn7CW349LDNrQEAoPeKehAyxmj27Nl6++239bvf/U4jR46MWD9y5Ej5fD6tXr3aWtbe3q5169Zp4sSJkqScnBy53e6ImtraWlVVVVk1EyZMUCAQ0ObNm62aTZs2KRAIRNRUVVWptrbWqlm1apU8Ho9ycnKifegxdeOVXcNjBCEAAC6cK9obfOSRR/TrX/9a//mf/6mUlBSrR8br9SopKUkOh0PFxcVatGiRRo0apVGjRmnRokUaOHCgCgoKrNoHHnhA8+bN09ChQ5WWlqb58+dr3LhxmjJliiRpzJgxmj59uoqKivTyyy9Lkh588EHl5+frqquukiTl5uZq7NixKiws1HPPPacjR45o/vz5Kioqisueni/jpiuH6udrPtEf/nJIwZBRgvPsw4oAAOAMon25mqQzvn75y19aNaFQyDz11FPG5/MZj8djbrnlFrN9+/aI7Rw/ftzMnj3bpKWlmaSkJJOfn2+qqyMvFz98+LCZNWuWSUlJMSkpKWbWrFmmoaEhombPnj0mLy/PJCUlmbS0NDN79mzT2tp63scTj5fPG2NMe2fQjHuqzAxf8J7Zsvuw3c0BACCunO/3t8MYLjs6m6amJnm9XgUCgbjrRZr75ja9+6f9emjSlXrim6Ptbg4AAHHjfL+/edZYL3brmK7bB6z5qPukcQAAcG4EoV5s8l9lKMHp0Cf1Lfr80NFzfwAAAEQgCPVi3oFuTTxx9dhv/rTf5tYAAND7EIR6uRlf7Xp8yLt/2s9dpgEA+JIIQr3c9GyfEl1OfVLfol11zXY3BwCAXoUg1MsNHuDWN666RJL09tZ9NrcGAIDehSDUB3w7p+vBsP9fxT61dgRtbg0AAL0HQagP+MboDF06JEkNxzq0cnvtuT8AAAAkEYT6hASnQ3f/dVev0K827rG5NQAA9B4EoT7ijq9lyeV0aFt1o6pqAnY3BwCAXoEg1EdkpAzQbeOGSZJe/P1fbG4NAAC9A0GoD3nkG1+RJK3cXqePD3ApPQAA50IQ6kOu8qXom9k+SdLS39ErBADAuRCE+pjZf9PVK/SbD/fro9omm1sDAEB8Iwj1MVf7vcq7ZpiMkX70mx08dgMAgLMgCPVBT35ztDwupzZ+dkRlVXV2NwcAgLhFEOqDLksdqO9OulKS9OP3dqq5tcPmFgEAEJ8IQn3U9yZdqay0JO0PtOof3/vI7uYAABCXCEJ9VFJigp7/u6/K4ZDe+uNerfnogN1NAgAg7hCE+rDxVwzVAzeNlCTN/48/ae+RYza3CACA+EIQ6uPmT7tK11zmVcOxDj30egVPpwcA4BQEoT5ugDtBL92To7TkRO3Y36Ti0koFQ1xSDwCARBDqFy4dkqQXZ12vxASnynbU6YfvbOf+QgAAiCDUb9x4xVD97K5r5XRIb27eqx+8U0XPEACg3yMI9SPfHDdMi/+fa+RwSL/eVK1HS7epvTNkd7MAALANQaifueOGLL1w93VyJzj03oe1uudfN6m+udXuZgEAYAuCUD+Uf41f/+99X9Mgj0ubdx/RjBfWq2LPEbubBQBAzBGE+qlJf3WJ/nP2TfpKxiAdaGrTt5eVa/F/7eLyegBAv0IQ6seuvGSQ3nnkJn3ruksVMtKydZ8q7+f/o9//ud7upgEAEBMEoX5ukMelJXdeq1/ce4PSB3n06cGjuv+XW3Tfv23Wrromu5sHAECPchhuKHNWTU1N8nq9CgQCGjx4sN3N6VGB4x1a+rtP9OqGz9UR7PpnkTs2U4984yv6atYQexsHAMCXcL7f3wShc+hPQShs96Gjev6//6yVVbUK/+vIGZ6qgr++XHnXDNMAd4K9DQQA4BwIQlHSH4NQ2F/qm/XS7z/Tf1bWqPPEzRe9SW7lXTNMeeOGafzINLkSGF0FAMQfglCU9OcgFFbf3Kr/+OM+vbm5WvsajlvLhyYnKvdqnyZfdYkmXjlUKQPcNrYSAICTCEJRQhA6KRgyKv/0sH67fb/KqurUcKzDWudyOnT95an6+qh05YxI1bVZQzQw0WVjawEA/RlBKEoIQmfWEQxpw6eHteajA/qfTw5p96GjEesTnA6N9qUoZ3iqrrt8iMYO8+qKS5LlZigNABADBKEoIQidn71HjumDTw5qw6eHtW1Pg/YHuj+2IzHBqa9kDNLoYSkaO2ywrswYpJFDk3VZahJzjQAAUUUQihKC0IWpDRzX1j2NqtjToA/3NWpXXbNa2jrPWOtyOnR52kCNSE/WyPRkjUhP1mVDkjRsyAD5hyRpMHOPAABfEkEoSghC0REKGdU0HtfO2ibtqm3WrromfXbwqD4/fFRtnaGzfjbF47JCkX9IkoYNHqD0FI/SB3l0SYpH6YMSlT7Iw2X9AAALQShKCEI9KxQyqmtq1e5DR7X70FF9fuioPj98TPsbj2t/4LgaT5mQfS4pA1y6ZFBXQBo6KFFDBiZqyEC3hiS5NWSgW94kt7xJJ5YNdGtIUqIGuJ1yOBw9eIQAADuc7/d3v7is58UXX9Rzzz2n2tpaXX311frpT3+qr3/963Y3C5KcTofV03PTV9K7rT/W3qn9ja2qDRzX/sbjqmls1YFAqw61tJ14tetgc5vagyE1t3aqubVTn502cftsEl1ODR7g1iBPggYNcCk50aVBHlfX7x6XUjxdP0/9vasuQQPcXa+kxAQNcDlP/EyQ00mwAoDeos8HobfeekvFxcV68cUXddNNN+nll1/WN7/5Te3cuVOXX3653c3DOQxMdOkrGYP0lYxBX1hjjFFTa6cOtbTpYHNXQDrc0q7A8Q41HutQ4/F2BY51qPF4hxqPnVzeGTJq7wydCFTRa3Oiy6kkd4IGuMM/TwSmcGhyO+VxJSgxwSm3y6HEhAS5XQ55EpxyJziV6Dr5MzH803XqOoc8rshal9OhBKdDLqdTrgRHxPuunw4CGgCcQZ8fGhs/fryuv/56vfTSS9ayMWPG6Pbbb1dJSck5P8/QWN9kjNGx9qAajrWrubVTLW1dr6NtnWo5/X1bp1raghHrjrZ3qrUjqNaOkI53BNV+jnlO8cDhkBWQ3E6nEr4gMCWEaxIiQ5TT0XVbBKfDIYfDoQSHTv7u7Prd6ThZ27VOSrCWn1Lj0Im6M/xu1XRf53BIDnVtt+uYHHIovPzEe8fpyx2nrO96r1Prw8vPuu2T+1fE9rp//ozbPuUzCted8nexfu/2Nzt3XfeR3QvY9hdsz3Hap75oFPn05V94fOdbd9Z9fvnjuxjR+p8P0RqBP/1vclHbiqP/beQd6I76hTEMjUlqb29XRUWFnnjiiYjlubm52rBhwxk/09bWpra2Nut9UxNPYO+LHA6HNeQVDcGQUVtnUMfbg2rtDHX97Oh6HT8lMLW2B63g1B4MWT87wj+DIbV1htQRNGrvDJ74ebK24wyfae8MKRgy6gwZ6+eZGCN1BI06gkativ/gBqD/eHz6VXp48lds2XefDkKHDh1SMBhUZmZmxPLMzEzV1dWd8TMlJSX60Y9+FIvmoQ9JcDo0MNEVF3fTNsYoZKTO0CkBKWjUEX4fNKcFp8gg1RGM/Fx4ech0vYyR9d4YKWjC67omv4d/Nyb8OXUtO/X3L/jMGetCXb8HjZG6/k/GmBM/I993re9qV7h94TrJdKs/9b2s9yfqTv39xLZD1uciP6+I96d+/szbtv5Wp/3dzrQ88s35fSZyH+aMy8/0/kK3233dF+3jLG35svv/ovN1xgUXJlrDJdEaeInm8E20xoJO/3dwoVw2Dt3b/1/tGDi9i9QY84Xdpk8++aQee+wx631TU5OysrJ6tH1ANIWHrRKc3E4AAM6lTweh9PR0JSQkdOv9qa+v79ZLFObxeOTxeGLRPAAAYLM+/VyDxMRE5eTkaPXq1RHLV69erYkTJ9rUKgAAEC/6dI+QJD322GMqLCzUDTfcoAkTJuiVV15RdXW1HnroIbubBgAAbNbng9Cdd96pw4cP6+mnn1Ztba2ys7O1cuVKDR8+3O6mAQAAm/X5+whdLO4jBABA73O+3999eo4QAADA2RCEAABAv0UQAgAA/RZBCAAA9FsEIQAA0G8RhAAAQL9FEAIAAP0WQQgAAPRbBCEAANBv9flHbFys8I23m5qabG4JAAA4X+Hv7XM9QIMgdA7Nzc2SpKysLJtbAgAAvqzm5mZ5vd4vXM+zxs4hFApp//79SklJkcPhiOq2m5qalJWVpb179/Icsx7EeY4NznNscJ5jh3MdGz11no0xam5ult/vl9P5xTOB6BE6B6fTqcsuu6xH9zF48GD+nywGOM+xwXmODc5z7HCuY6MnzvPZeoLCmCwNAAD6LYIQAADotwhCNvJ4PHrqqafk8XjsbkqfxnmODc5zbHCeY4dzHRt2n2cmSwMAgH6LHiEAANBvEYQAAEC/RRACAAD9FkEIAAD0WwQhm7z44osaOXKkBgwYoJycHP3P//yP3U3qNUpKSvS1r31NKSkpysjI0O23364///nPETXGGC1cuFB+v19JSUmaPHmyduzYEVHT1tamOXPmKD09XcnJyZo5c6b27dsXy0PpVUpKSuRwOFRcXGwt4zxHT01Nje655x4NHTpUAwcO1LXXXquKigprPef64nV2duqHP/yhRo4cqaSkJF1xxRV6+umnFQqFrBrO84X54IMPNGPGDPn9fjkcDr3zzjsR66N1XhsaGlRYWCiv1yuv16vCwkI1NjZeXOMNYq60tNS43W7zi1/8wuzcudM8+uijJjk52ezZs8fupvUK06ZNM7/85S9NVVWVqaysNHl5eebyyy83LS0tVs3ixYtNSkqKWb58udm+fbu58847zbBhw0xTU5NV89BDD5lLL73UrF692mzdutV84xvfMF/96ldNZ2enHYcV1zZv3mxGjBhhrrnmGvPoo49ayznP0XHkyBEzfPhwc//995tNmzaZ3bt3m/fff9/85S9/sWo41xfvH//xH83QoUPNe++9Z3bv3m3+4z/+wwwaNMj89Kc/tWo4zxdm5cqV5gc/+IFZvny5kWRWrFgRsT5a53X69OkmOzvbbNiwwWzYsMFkZ2eb/Pz8i2o7QcgGf/3Xf20eeuihiGWjR482TzzxhE0t6t3q6+uNJLNu3TpjjDGhUMj4fD6zePFiq6a1tdV4vV6zbNkyY4wxjY2Nxu12m9LSUqumpqbGOJ1OU1ZWFtsDiHPNzc1m1KhRZvXq1WbSpElWEOI8R8+CBQvMzTff/IXrOdfRkZeXZ77zne9ELPvWt75l7rnnHmMM5zlaTg9C0TqvO3fuNJLMxo0brZry8nIjyezateuC28vQWIy1t7eroqJCubm5Ectzc3O1YcMGm1rVuwUCAUlSWlqaJGn37t2qq6uLOMcej0eTJk2yznFFRYU6Ojoiavx+v7Kzs/k7nOaRRx5RXl6epkyZErGc8xw97777rm644QZ9+9vfVkZGhq677jr94he/sNZzrqPj5ptv1po1a/Txxx9Lkv70pz9p/fr1uu222yRxnntKtM5reXm5vF6vxo8fb9XceOON8nq9F3XueehqjB06dEjBYFCZmZkRyzMzM1VXV2dTq3ovY4wee+wx3XzzzcrOzpYk6zye6Rzv2bPHqklMTFRqamq3Gv4OJ5WWlmrr1q3asmVLt3Wc5+j57LPP9NJLL+mxxx7T//7f/1ubN2/W3Llz5fF4dO+993Kuo2TBggUKBAIaPXq0EhISFAwG9cwzz+juu++WxL/pnhKt81pXV6eMjIxu28/IyLioc08QsonD4Yh4b4zptgznNnv2bH344Ydav359t3UXco75O5y0d+9ePfroo1q1apUGDBjwhXWc54sXCoV0ww03aNGiRZKk6667Tjt27NBLL72ke++916rjXF+ct956S6+//rp+/etf6+qrr1ZlZaWKi4vl9/t13333WXWc554RjfN6pvqLPfcMjcVYenq6EhISuqXX+vr6bmkZZzdnzhy9++67Wrt2rS677DJruc/nk6SznmOfz6f29nY1NDR8YU1/V1FRofr6euXk5MjlcsnlcmndunX6+c9/LpfLZZ0nzvPFGzZsmMaOHRuxbMyYMaqurpbEv+lo+fu//3s98cQTuuuuuzRu3DgVFhbq+9//vkpKSiRxnntKtM6rz+fTgQMHum3/4MGDF3XuCUIxlpiYqJycHK1evTpi+erVqzVx4kSbWtW7GGM0e/Zsvf322/rd736nkSNHRqwfOXKkfD5fxDlub2/XunXrrHOck5Mjt9sdUVNbW6uqqir+Difceuut2r59uyorK63XDTfcoFmzZqmyslJXXHEF5zlKbrrppm63gPj44481fPhwSfybjpZjx47J6Yz82ktISLAun+c894xondcJEyYoEAho8+bNVs2mTZsUCAQu7txf8DRrXLDw5fP/+q//anbu3GmKi4tNcnKy+fzzz+1uWq/wve99z3i9XvP73//e1NbWWq9jx45ZNYsXLzZer9e8/fbbZvv27ebuu+8+46Wal112mXn//ffN1q1bzd/8zd/0+0tgz+XUq8aM4TxHy+bNm43L5TLPPPOM+eSTT8wbb7xhBg4caF5//XWrhnN98e677z5z6aWXWpfPv/322yY9Pd08/vjjVg3n+cI0Nzebbdu2mW3bthlJZsmSJWbbtm3WbWGidV6nT59urrnmGlNeXm7Ky8vNuHHjuHy+t/qXf/kXM3z4cJOYmGiuv/5669JvnJukM75++ctfWjWhUMg89dRTxufzGY/HY2655Razffv2iO0cP37czJ4926SlpZmkpCSTn59vqqurY3w0vcvpQYjzHD2/+c1vTHZ2tvF4PGb06NHmlVdeiVjPub54TU1N5tFHHzWXX365GTBggLniiivMD37wA9PW1mbVcJ4vzNq1a8/43+X77rvPGBO983r48GEza9Ysk5KSYlJSUsysWbNMQ0PDRbXdYYwxF96fBAAA0HsxRwgAAPRbBCEAANBvEYQAAEC/RRACAAD9FkEIAAD0WwQhAADQbxGEAABAv0UQAgAA/RZBCADOYcSIEfrpT39qdzMA9ACCEIC4cv/99+v222+XJE2ePFnFxcUx2/err76qIUOGdFu+ZcsWPfjggzFrB4DYcdndAADoae3t7UpMTLzgz19yySVRbA2AeEKPEIC4dP/992vdunX62c9+JofDIYfDoc8//1yStHPnTt12220aNGiQMjMzVVhYqEOHDlmfnTx5smbPnq3HHntM6enpmjp1qiRpyZIlGjdunJKTk5WVlaWHH35YLS0tkqTf//73+l//638pEAhY+1u4cKGk7kNj1dXV+tu//VsNGjRIgwcP1h133KEDBw5Y6xcuXKhrr71Wv/rVrzRixAh5vV7dddddam5u7tmTBuBLIwgBiEs/+9nPNGHCBBUVFam2tla1tbXKyspSbW2tJk2apGuvvVZ//OMfVVZWpgMHDuiOO+6I+Pxrr70ml8ulP/zhD3r55ZclSU6nUz//+c9VVVWl1157Tb/73e/0+OOPS5ImTpyon/70pxo8eLC1v/nz53drlzFGt99+u44cOaJ169Zp9erV+vTTT3XnnXdG1H366ad655139N577+m9997TunXrtHjx4h46WwAuFENjAOKS1+tVYmKiBg4cKJ/PZy1/6aWXdP3112vRokXWsn/7t39TVlaWPv74Y/3VX/2VJOkrX/mKnn322YhtnjrfaOTIkfrxj3+s733ve3rxxReVmJgor9crh8MRsb/Tvf/++/rwww+1e/duZWVlSZJ+9atf6eqrr9aWLVv0ta99TZIUCoX06quvKiUlRZJUWFioNWvW6Jlnnrm4EwMgqugRAtCrVFRUaO3atRo0aJD1Gj16tKSuXpiwG264odtn165dq6lTp+rSSy9VSkqK7r33Xh0+fFhHjx497/1/9NFHysrKskKQJI0dO1ZDhgzRRx99ZC0bMWKEFYIkadiwYaqvr/9Sxwqg59EjBKBXCYVCmjFjhn7yk590Wzds2DDr9+Tk5Ih1e/bs0W233aaHHnpIP/7xj5WWlqb169frgQceUEdHx3nv3xgjh8NxzuVutztivcPhUCgUOu/9AIgNghCAuJWYmKhgMBix7Prrr9fy5cs1YsQIuVzn/5+wP/7xj+rs7NQ//dM/yens6gz/v//3/55zf6cbO3asqqurtXfvXqtXaOfOnQoEAhozZsx5twdAfGBoDEDcGjFihDZt2qTPP/9chw4dUigU0iOPPKIjR47o7rvv1ubNm/XZZ59p1apV+s53vnPWEHPllVeqs7NTL7zwgj777DP96le/0rJly7rtr6WlRWvWrNGhQ4d07NixbtuZMmWKrrnmGs2aNUtbt27V5s2bde+992rSpElnHI4DEN8IQgDi1vz585WQkKCxY8fqkksuUXV1tfx+v/7whz8oGAxq2rRpys7O1qOPPiqv12v19JzJtddeqyVLlugnP/mJsrOz9cYbb6ikpCSiZuLEiXrooYd055136pJLLuk22VrqGuJ65513lJqaqltuuUVTpkzRFVdcobfeeivqxw+g5zmMMcbuRgAAANiBHiEAANBvEYQAAEC/RRACAAD9FkEIAAD0WwQhAADQbxGEAABAv0UQAgAA/RZBCAAA9FsEIQAA0G8RhAAAQL9FEAIAAP3W/w8bHUj9O8+rzAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def obj_fn(Y, X, k):\n",
    "    \"\"\"Objective function\"\"\"\n",
    "    return np.sum((Y - (X.dot(k)))**2)\n",
    "    ## We can change our objective function!\n",
    "    #return np.sum(np.abs(Y - (X.dot(k)))) # This one is the sum of absolute differences\n",
    "    \n",
    "## Step 1: Choose a random k\n",
    "k = np.random.randn(d)\n",
    "\n",
    "lam = 0.001 # Step size; I want to take steps against the gradient\n",
    "delta = 0.001 ## This is the \"run\" in my derivative approximation\n",
    "n_steps = 1000\n",
    "f_over_iters = np.zeros(n_steps)\n",
    "for it in range(n_steps): # Keep repeating these gradient steps over and over again\n",
    "    ## Step 1: Approximately compute the gradient (NOTE: You can replace this with\n",
    "    ## an analytical expression for the gradient)\n",
    "    deriv_k = np.zeros(len(k))\n",
    "    for j in range(len(k)):\n",
    "        k_new = np.array(k) # Make a copy of k, but only change element j by delta\n",
    "        k_new[j] = k_new[j] + delta # Hold everything else fixed, but change this one element\n",
    "        rise = obj_fn(Y, X, k_new) - obj_fn(Y, X, k) ## Compute the rise\n",
    "        deriv_k[j] = rise / delta # Partial derivative is rise over run\n",
    "    ## Step 2: Go against the gradient!\n",
    "    k = k - lam*deriv_k # Go against the gradient!\n",
    "    ## Step 3: Evaluate the objective function again so we can keep\n",
    "    ## track of how it's changing\n",
    "    f_over_iters[it] = obj_fn(Y, X, k_new)\n",
    "\n",
    "plt.plot(f_over_iters)\n",
    "plt.xlabel(\"Iteration\")\n",
    "plt.ylabel(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "02e9635f-23fe-45a0-9e01-30b5fcf2089b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-33.5 -18.  -40.5 -86.5   9.  -20. ]\n"
     ]
    }
   ],
   "source": [
    "print(kgt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "422739ec-6410-4113-bc60-16e374446b59",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-32.80600634 -19.15002682 -40.57534396 -86.26318393   7.44236825\n",
      " -20.01528559]\n"
     ]
    }
   ],
   "source": [
    "print(k)  ## This is really close!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28079ebb-0820-43b6-a981-bb6bbf178b9c",
   "metadata": {},
   "source": [
    "Now that I've solved for my model, I can actually put new data points into it and get estimates from them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "067082ca-5f05-493d-8eff-b1ea7255d2c4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}