path: root/jax_optimization_tutorial/linear-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3387e989",
   "metadata": {},
   "source": [
    "# Linear Regression with JAX\n",
    "Source: https://danielrothenberg.com/blog/2020/Sep/jax-first-steps-pt1/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "36893279",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import jax.numpy as jnp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "565bff63",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import random\n",
    "\n",
    "def make_key():\n",
    "    \"\"\" Helper function to generate a key for jax's parallel PRNG \n",
    "    using standard numpy random functions. \n",
    "\n",
    "    \"\"\"\n",
    "    seed = np.random.randint(2**16 - 1)\n",
    "    return random.PRNGKey(seed)\n",
    "\n",
    "n = 100\n",
    "rands = random.uniform(make_key(), shape=(n, ), minval=-1, maxval=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "11fa0b53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]]\n",
      "[[ 0.56719136  1.        ]\n",
      " [ 0.90619564  1.        ]\n",
      " [ 0.8181107   1.        ]\n",
      " [-0.07588315  1.        ]\n",
      " [ 0.65934706  1.        ]\n",
      " [ 0.70064425  1.        ]\n",
      " [ 0.28244472  1.        ]\n",
      " [-0.07778668  1.        ]\n",
      " [ 0.9808245   1.        ]\n",
      " [ 0.05257344  1.        ]\n",
      " [-0.38835692  1.        ]\n",
      " [ 0.74025965  1.        ]\n",
      " [-0.8808477   1.        ]\n",
      " [-0.80188036  1.        ]\n",
      " [ 0.7314439   1.        ]\n",
      " [-0.69432974  1.        ]\n",
      " [ 0.6456201   1.        ]\n",
      " [-0.13038063  1.        ]\n",
      " [ 0.7454684   1.        ]\n",
      " [-0.17364717  1.        ]\n",
      " [ 0.8420925   1.        ]\n",
      " [ 0.32458067  1.        ]\n",
      " [-0.23971272  1.        ]\n",
      " [-0.841902    1.        ]\n",
      " [-0.9334514   1.        ]\n",
      " [ 0.8027196   1.        ]\n",
      " [-0.76277256  1.        ]\n",
      " [-0.3472433   1.        ]\n",
      " [-0.18131876  1.        ]\n",
      " [-0.81412005  1.        ]\n",
      " [ 0.742816    1.        ]\n",
      " [ 0.29102564  1.        ]\n",
      " [-0.12259459  1.        ]\n",
      " [ 0.7445142   1.        ]\n",
      " [ 0.4648075   1.        ]\n",
      " [ 0.97792006  1.        ]\n",
      " [ 0.98037124  1.        ]\n",
      " [-0.74772596  1.        ]\n",
      " [-0.9600477   1.        ]\n",
      " [ 0.29590702  1.        ]\n",
      " [ 0.8539684   1.        ]\n",
      " [ 0.7292776   1.        ]\n",
      " [-0.35773468  1.        ]\n",
      " [ 0.06983089  1.        ]\n",
      " [-0.1719327   1.        ]\n",
      " [ 0.3244648   1.        ]\n",
      " [-0.66061735  1.        ]\n",
      " [-0.9326389   1.        ]\n",
      " [-0.48183155  1.        ]\n",
      " [ 0.60397077  1.        ]\n",
      " [-0.24433494  1.        ]\n",
      " [ 0.9065745   1.        ]\n",
      " [-0.53917074  1.        ]\n",
      " [-0.26952457  1.        ]\n",
      " [-0.03242016  1.        ]\n",
      " [-0.6970415   1.        ]\n",
      " [-0.8621006   1.        ]\n",
      " [ 0.64250565  1.        ]\n",
      " [ 0.3600428   1.        ]\n",
      " [ 0.1736505   1.        ]\n",
      " [ 0.8175428   1.        ]\n",
      " [-0.4976089   1.        ]\n",
      " [ 0.99707437  1.        ]\n",
      " [ 0.61124754  1.        ]\n",
      " [ 0.4626403   1.        ]\n",
      " [ 0.7739856   1.        ]\n",
      " [ 0.18812919  1.        ]\n",
      " [ 0.9244573   1.        ]\n",
      " [ 0.55314875  1.        ]\n",
      " [ 0.99845624  1.        ]\n",
      " [ 0.9773083   1.        ]\n",
      " [ 0.87139964  1.        ]\n",
      " [ 0.32497478  1.        ]\n",
      " [ 0.76745677  1.        ]\n",
      " [-0.75871634  1.        ]\n",
      " [-0.10237694  1.        ]\n",
      " [-0.10543108  1.        ]\n",
      " [-0.19929266  1.        ]\n",
      " [-0.67447877  1.        ]\n",
      " [-0.54395914  1.        ]\n",
      " [ 0.4002776   1.        ]\n",
      " [-0.7200577   1.        ]\n",
      " [ 0.17349315  1.        ]\n",
      " [ 0.5890646   1.        ]\n",
      " [ 0.32220316  1.        ]\n",
      " [ 0.870471    1.        ]\n",
      " [-0.4356978   1.        ]\n",
      " [-0.54046607  1.        ]\n",
      " [ 0.60255     1.        ]\n",
      " [ 0.08106732  1.        ]\n",
      " [ 0.7314334   1.        ]\n",
      " [-0.59192777  1.        ]\n",
      " [-0.43610024  1.        ]\n",
      " [ 0.36425257  1.        ]\n",
      " [-0.91518784  1.        ]\n",
      " [ 0.3074093   1.        ]\n",
      " [-0.43152308  1.        ]\n",
      " [-0.32482266  1.        ]\n",
      " [ 0.7852454   1.        ]\n",
      " [ 0.57705736  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "# We have to use functions to update an array unlike numpy as JAX does not mutate arrays in-place\n",
    "\n",
    "x = jnp.ones((n, 2))\n",
    "print(x)\n",
    "x = x.at[:, 0].set(rands)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9a1e2fb",
   "metadata": {},
   "source": [
    "# Simple first model\n",
    "\n",
    "$$\n",
    "    y = X\\beta + \\epsilon\n",
    "$$\n",
    "\n",
    "where $\\beta = [\\beta_0, \\beta_1]$ where $\\beta_0$ is an offset bias and $\\beta_1$ a relating x and y, and $\\epsilon$ is an uncorrelated noise modeled as a normal distribution, $\\mathcal{N}(0, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7bf67529",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix our true slope and bias\n",
    "slope, bias = 3.0, 2.0\n",
    "beta_true = jnp.array((slope, bias))\n",
    "eps = random.normal(make_key(), shape=(n,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5a6d3ce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = (x @ beta_true) + eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "462c84ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fbd01976f40>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa+ElEQVR4nO3df4wcZ3kH8O/D5RzOUHIxOVGyydlOG5ymdcnBiqY9CYhBOPxQck0CGEEbflQulCKg1PTcUAlVRb42KpQKJBqFH22JklA7XF0F6iaco6oRdrnrGUziODFJQ7IJxBAuNOSSnJ2nf+ysb25vZnZ25p133nfm+5Es383uzr47u/fMu8/7vO+IqoKIiPz1vLIbQERE+TCQExF5joGciMhzDORERJ5jICci8txpZTzpWWedpRs2bCjjqYmIvDU3N/cTVR3p3l5KIN+wYQNmZ2fLeGoiIm+JyINR25laISLyHAM5EZHnGMiJiDzHQE5E5DkGciIiz5VStUJEZMP0fAvX7juKRxYWcfbwEHZs3YSJsUbZzTKOgZyIKml6voWdtxzG4tJJAEBrYRE7bzkMAKuCue8Bn4GciCrp2n1HTwXxjsWlk7h239EVQTou4M8++Dj233Pci+DOQE5ElfTIwmKq7XEB/4YDP0Tnag1JvXkXGBnsFJFhEdktIveIyBER+W0T+yUiyurs4aFU2+MCfvcldzq9eReZqlr5LIB/V9ULALwcwBFD+yUiymTH1k0YGhxYsW1ocAA7tm5asS0u4EeJC/plyx3IReQMAK8G8EUAUNVnVXUh736JiPKYGGtg1xWb0RgeggBoDA9h1xWbV6VGogK+xOyzn6Bvk4kc+UYAxwF8WUReDmAOwIdV9RfhO4nIdgDbAWB0dNTA0xIRJZsYa/TMaXduD1etXHLBCPbMtVbkzqN6866QvBdfFpEmgAMAxlX1oIh8FsDPVfUv4h7TbDaVqx8SkctcLEkUkTlVbXZvN9EjfxjAw6p6MPh9N4BJA/slIipNmt68K3LnyFX1RwAeEpHOd47XAbg7736JiCgdU3XkHwJwg4isAXA/gPcY2i8REfVgJJCr6iEAq/I2RERUPK5+SETkOU7RJyKyoMgqGAZyIqKC9bMSYxZMrRARFSxpJUYT2CMnokpwcQJPR9qVGLNiICci76VNXYSD/fDaQagCTywuFR74zx4eQisiaJtau4WpFSLyXprURSfYtxYWoQB+9tQSFhaXoFgO/NPzrULal3YlxqwYyInIe2lSF1HBPqzI9cbTrsSYFVMrROS9NKmLNPnoItcbL3LtFvbIich7aVIXafLRrq433gsDORF5L03qIirYhxW53vj0fAvjUzPYOHkrxqdmjOfimVohokrolbrovoBEr6qVqHLG8OPPGBqECLDwVHLVS9GTgQADF5bIgheWICKXdQdfABgcEECBpeeiY+bQ4MCqbwHT8y187GvfxcmIONsYHsKdk1v6aleRF5YgIiqdyQlBURUuSyeTO72dqpfOc35i+jBuOPBDxD3K5MAqAzkRec90+iJrkO08bnq+lRjEAbMDqxzsJCLvmV7LJGuQ7Tzu2n1HE4O46YFVBnIi8p7ptUx6VbhECQfnpOcdEDE6GQhgICeiCojrQWftWUeVMw4PDcbev7vcMe55BcDfvu3lxicGMUdORN7bsXXTqiqTvOmL7nLGqEqWqEqVuPYIgHdePFrI7E4GciLyXneNeBGrGfbzHDbaE8Y6ciKqHZfXLk/COnIiItiZaWkbAzkRVVJcrzupVJGBnIioROHAfcbQIH7x7IlTszHDve6iL7tWBpYfEpH3uq/+s7C4tGpKfafXbbpU0QUM5ETkvV5X/+l4ZGGx8MuudRS9dG0YUytE5LXp+Vbk1YGinD08ZKU00PaAKgM5EXmrEzDTCPe6i7zsGpC89ksRz8vUChF5KymlMvg8wZlrBwu52HEvtgdU2SMnIm8lBcZr32p+TZO00lwM2iT2yInICVkGB+MCYyOUCy+DrQHVDgZyIipdd/lgZ3CwVzC3HTDTSnMxaJOYWiGi0mUdHLS9OFU/ih5QDWMgJ6LS5RkctBkwXcVATlRhvqzyZ3twsGqYIyeqqKx55zK4muv2BQM5UUWZviBxkWwPDlaNsdSKiAwAmAXQUtW3mNovEWXj2yp/zHVnZzJH/mEARwC8yOA+iSgj5p3758uYQjcjqRUROQfAmwFcb2J/5A6bK7iRWcw798enMYVupnrkfwfg4wB+ydD+yAFVvCRWVSX1JH3sYca9niJ7zD5fOSh3IBeRtwB4TFXnROS1CffbDmA7AIyOjuZ9WrLA5w92nfQ64fr2XsW9ntkHH8eeuVZhHQvfxhTCTKRWxgFcJiL/C+AmAFtE5Kvdd1LV61S1qarNkZERA09LRfP5g10nPlWnpBH3em48+FChr7OIKwfZSk3mDuSqulNVz1HVDQC2AZhR1XflbhmVroqXxKqiqp1w49p9UjVye2th0UiwND2mYDPnzjpyisXBMj+4fMI1uaLhgEjkdgESg2XaNpiuZbf5TcnoFH1VvQPAHSb3SeXxebCsTnZs3bQipwyUe8LtDEi2FhYhADr96NbCIj568yF85OZDaMR8lqbnW3jq2ROr9jk0OIArX9lYkSMHsGL/HeFxnH4H7E2OKdj8psS1ViiRj4NldTMx1sDsg4/jxoMP4aQqBkRw5SvLed+6A2d3kA0H9e6A2v3YjuGhQXzysl/HxFgDzfXrVnQs4q7V2QmWZQ7Y26zjZ2qFyHPT8y3smWudyiGfVMWeuVYp9c9pr2YPrE4zxD32BaefdiroTow1cOfkFjww9WbcObkFjR5ppTLHD2ymJhnIiTznUtVKvwEyfP8sQbdXsCxz/MDm+jFMrRB5zqWqlaR0R9z9ez02Kej2Gscpe/zAVmqSgZzIcy6tqRIVODsDkt0Dk90BNWvQTQqWdRmwZyAn8lzZvc6wpMDZa3p9UUG3DgP2ojFF9kVqNps6Oztr/XmJqsrXVfuoPyIyp6rN7u3skRNVgG+9Tp54zGIgJyKruKqmeSw/JCKrXCqXrAoGciKyyqVyyapgICciq1xe5MtXDOREZBVX1TSPg51EZFVdJunYxEBORNb5Vi7pOgZyckod6ovr8BqT1P31F4GBnJxRh/riOrzGJHV//UXhYCc5ow71xXV4jUnq/vqLwh45OaMO9cV1eI1Jinj9TNWwR04OqUN9cR1eYxLTrz/qSvUfvfkQPjF9OEcr/cNATs6oQ31xHV5jkksuGOlrey9RqRoFcMOBH5ZyqbuyMLVCzqhDfbErr7GsdMT+e473tb2XuJSMAlYusOwKBnJySh3qi8t+jWVWjpjOkSddWq4u4w4AUytEtVNm5YjpHPmOrZsgfT5XFTGQE9VMmZUzpscIJsYaeOfFo6uCeZ3GHQCmVoi812++u8yLNRcxRvBXE5vRXL+u9HGHMvGanUQe6853A+3e6K4rNscGsiyPITfEXbOTqRUij2XJd0+MNbDris1oDA9BADSGhxjEPcfUClEfXJtFmDXfnbZyxrXXS9EYyIlScnHBp7h89/Dawdz7dvH1UjSmVohScnHBpx1bN2FwYHUB3pNPn8g9s9HF10vRGMiJUnJxwauJsQZesGb1F+ul5zR3wHXx9VI0BnKilFxd8OqJxaXI7XkDrquvl1ZjICdKydUFr4oKuK6+XlqNg51UCyaqL1xZ8Krbjq2bIuvC8wZcV18vrcZATpVnovqi+0Twmbdf5ExAKzLglr3AF6XDQF4xLtT9utCGsKTqi7S11KZPBKaPCQNuvTGQV4gLdb8225A2OOatvnDhRECUJPdgp4icKyL7ReRuEblLRD5somHUPxfqfm21IeoSXztvORxZO513MLDIE4HLpudbGJ+awcbJWzE+NVOrK+74xkTVygkAH1PVCwFcDOCDInKhgf1Sn1yo+7XVhn6CY97qi7JPBGXo50RJ5csdyFX1UVX9n+Dn/wNwBAC/L5bAhbpfW23oJzjmXSQq74kgbrq8iWn0RfH1W0RdGc2Ri8gGAGMADkbcth3AdgAYHR01+bQUKKoMzcU29Lumdp7BwLxVIXErRRe9gnSeAVYfv0XUmbFALiIvBLAHwEdU9efdt6vqdQCuA9rrkZt6XlrmQt2vrTbYPmnlORHEzbyM225C3gHWMi8+Qf0zEshFZBDtIH6Dqt5iYp+UjQtlaDba4MJJK60ygmLeShuTJ0rXylGrKHcgFxEB8EUAR1T10/mbRJSOCyetNMpIeeVNjZg6UbL00g4TPfJxAL8H4LCIHAq2/bmqfsPAvom8V8a3BxPfAkycKPN+M6B0cgdyVf0vYNVFrIlSqcvXbtvfHlwY+AY4aGoLZ3YifTCpS9CxJe5r9+yDj2P/Pcd5nHNwZQyBg6Z2iBZdAxWh2Wzq7Oys9eeNkvaK4rzyuHnjUzORf+QCIPyp5HFeyacOBf9uzBKROVVtdm+v/XrkaSc+cIKEeXFfr7u7FjzOy3ybcZl3MhalU/vUStocHnN95sV97Y7C49zm4+ChL9VFPqt9jzztlHIXpr9XTdTU97hRcx7ndm887sTHE1291T6Qp11Hg5e9Mi/qa/c7Lx7lcY7QSanE4Ymu3mqfWkk7uu9KFUDVRH3tbq5fx+PcJSqlEnbJBSMWW0OuqX3VCpEPNk7eumoQOKwxPIQ7J7dYaw+VI65qpfY98rrxqXSNlvUaGGaOvN5qnyOvE99K12hZ1BhNGHPk9cYeuQWu9IJ9LF3rlyvH2rTOa/jk3ruw0LX8LQeDiYG8YC6t/hb31TxtLbfrXDrWRegMDFf1ZEXZMZAXzKVe8IAITkYMbg9I8ppnLgWOpLa4dKyLxAk21I2BvGAmZ4TmDahRQTxpe+c5Xenl9moLZ99SXXGws2CmZoSaGKhsxDxn3HbArTVmerWFs2+prhjIC2ZqRqiJgJqlLS71cnu1hbNvqa6YWimYqRmhJgJqlraYWk/aRJ69V1uKmH3r0vgAURwGcgtMDE6ZCqj9tsXElWZM5dnTtMXkQKBL4wNESZha8URZaQMT60mbyrPbXtvapfEBoiTskXuizEW78vZyTebZbZbeuTQ+QJSEgZwK5+t1G11vN/P31MHUiidsrJMyPd/C+NQMNk7eivGpGWP79rWaxOV2c90cCmOP3BNFz1rsZ2Cv356gqbSQ7R6oy2vQ12UWK6XDQO6JovO1aQND1kqOvLntsipIXJ0Oz/w9hTG14omiZy2mDQxlVXKwgmQlzmKlMAZyTxSdr00bGMrqCbIHupLL+Xuyj4HcE0XXUKcNDEX2BJMGW9kDXcl2TT25jdfspFPSDCZ256qBdsDPG0R67beo5yXyCa/Z6YGy64LTDOwVVcnRa7DV5QoSorIxkDsibVVG2cG+0x7Tz5kmB+5qBQlR2SobyLMEvLjH2Aieacr/qryIk+uzKIlcVsnBziyz3uIe84npw1Zm0KXpkVa5BI9VGETZVTKQZwl4cY+58eBDVoJnmqqMKpfgsQqDKLtKplayBLy42+KuZ2k6eKZZa7vq6QfmwImyqWSPPEvNcdxtcVeYNx080/RImX4goiiV7JFnuapN3GOufGUDe+Zaua6Qk1avHilL8IgoSiUDeZaAl/SY5vp1zgRPph+IqJuRmZ0icimAzwIYAHC9qk4l3Z8zO4mI+hc3szN3jlxEBgB8HsAbAVwI4B0icmHe/RIRUTomUiuvAnBMVe8HABG5CcDlAO42sG+iWC7MciVygYlA3gDwUOj3hwH8VvedRGQ7gO0AMDo6auBpl/EPun6qPMuVqF/Wyg9V9TpVbapqc2RkxNh+ee3CeqryLFeifpkI5C0A54Z+PyfYZgX/oOupyrNcifplIpB/B8D5IrJRRNYA2AZgr4H9psI/6HrihSaIluUO5Kp6AsAfA9gH4AiAr6nqXXn3mxb/oNuSrq5TRZzlSrTMSI5cVb+hqi9T1V9R1U+Z2Gda/IOu5zgBF9kiWub9zE5OW0+3lnkVdc9y7XwrqevngOrL+0AOcNo6xwlYjkj1VolA7oqy6tldX97WlSssEVUVA3kO4QA1vHYQTz59AkvPtdeusdkjzLLaoy22esr8VkJ1Vsn1yG3oHmD82VNLp4J4h616dpcH/mzV+bN6ieqMPfKMogJUFFs9QlfHCWz1lF3+VkJUNPbIM0obiOreI7TVU3b5WwlR0dgjzyhugDGMPUK7PWVXv5UQFY098oyiJiINDgiGhwbZIwxhT5moeOyRZ8SJSOmxp0xULAbyFOLqoH0MUFy7nah6GMh7qNKMwSq9FiJaxhx5D1Va77xKr4WIlrFH3kOeOmjX0hic/UhUTeyR95C1DtrFpWU5+5GomhjIe8i63rmLaQyu3U5UTd6mVmylLbKWGbqYxmDJJFE1eRnIbVdfZCkzdHVpWR9LJokomZepFRfTFt2YxiAiW7zskbuYtujGNAYR2eJlIHc1bdGNaQwissHL1ArTFkREy7zskTNtQUS0zMtADjBtQUTU4WVqhYiIljGQExF5joGciMhzDORERJ5jICci8hwDORGR5xjIiYg8x0BOROQ5byYEuXbZNCIiV3gRyHn1dyKieF6kVnxYf5yIqCxeBHIf1h8nIiqLF6kVX9YfZx6fiMqQq0cuIteKyD0i8j0R+bqIDBtq1wo+rD/eyeO3FhahWM7jT8+3ym4aEVVc3tTKbQB+Q1V/E8C9AHbmb9JqE2MN7LpiMxrDQxAAjeEh7Lpis1O9XebxiagsuVIrqvofoV8PALgqX3Piub7+OPP4RFQWk4Od7wXwzbgbRWS7iMyKyOzx48cNPq0b4vL1ruXxiah6egZyEbldRL4f8e/y0H2uAXACwA1x+1HV61S1qarNkZERM61PaXq+hfGpGWycvBXjUzOF5K19yOMTUTX1TK2o6uuTbheRdwN4C4DXqaoaapcxtiYT8TqiRFSWXDlyEbkUwMcBvEZVnzLTJLOSBiFNB1nX8/hEVE1568g/B+B0ALeJCAAcUNX3525VCmlrtjkISURVl7dq5VdNNaQf/aRLfJlMRESUlRdT9Lv1U7MdNQg5OCD4xTMnCh38JCKyxYsp+t36SZd0D0IOrx3Ek0+fwMLiEgCupEhE/vOyR95vzfbEWAN3Tm7BA1Nvxto1p2HpuZXFNZyBSUQ+8zKQ56nZ5uAnEVWNl4E8z9ornIFJRFXjZY4cyF6zvWPrphUVLwBnYBKR37wN5FlxBiYRVU3tAjnAGZhEVC1e5siJiGgZAzkRkecYyImIPMdATkTkOQZyIiLPSRnXghCR4wAe7HG3swD8xEJzsmDbsnG1ba62C2Dbsqpq29ar6qpLrJUSyNMQkVlVbZbdjihsWzauts3VdgFsW1Z1axtTK0REnmMgJyLynMuB/LqyG5CAbcvG1ba52i6AbcuqVm1zNkdORETpuNwjJyKiFBjIiYg8V2ogF5G3ishdIvKciMSW44jIpSJyVESOichkaPtGETkYbL9ZRNYYbNs6EblNRO4L/j8z4j6XiMih0L+nRWQiuO0rIvJA6LaLbLYtuN/J0PPvDW0v5LilPGYXici3g/f9eyLy9tBtxo9Z3GcndPvpwTE4FhyTDaHbdgbbj4rI1rxtydC2PxGRu4Pj9C0RWR+6LfK9tdi2d4vI8VAb/iB029XBZ+A+Ebnacrs+E2rTvSKyELqt6GP2JRF5TES+H3O7iMjfB23/noi8InRbvmOmqqX9A/BrADYBuANAM+Y+AwB+AOA8AGsAfBfAhcFtXwOwLfj5CwA+YLBtfwNgMvh5EsBf97j/OgCPA1gb/P4VAFcVdNxStQ3AkzHbCzluadoF4GUAzg9+PhvAowCGizhmSZ+d0H3+CMAXgp+3Abg5+PnC4P6nA9gY7GfActsuCX2ePtBpW9J7a7Ft7wbwuYjHrgNwf/D/mcHPZ9pqV9f9PwTgSzaOWbD/VwN4BYDvx9z+JgDfBCAALgZw0NQxK7VHrqpHVLXXVY9fBeCYqt6vqs8CuAnA5SIiALYA2B3c7x8BTBhs3uXBPtPu+yoA31TVpwy2IU6/bTul4OPWs12qeq+q3hf8/AiAxwCsmqlmSORnJ6HNuwG8LjhGlwO4SVWfUdUHABwL9metbaq6P/R5OgDgHIPPn6ttCbYCuE1VH1fVnwG4DcClJbXrHQBuNPTcPanqf6LdmYtzOYB/0rYDAIZF5KUwcMx8yJE3ADwU+v3hYNuLASyo6omu7aa8RFUfDX7+EYCX9Lj/Nqz+0Hwq+Ar1GRE5vYS2PV9EZkXkQCflg2KPW1/HTERehXbP6gehzSaPWdxnJ/I+wTF5Au1jlOaxRbct7H1o9+Y6ot5b2227MnivdovIuX0+tsh2IUhDbQQwE9pc5DFLI679uY9Z4VcIEpHbAfxyxE3XqOq/Fv38SZLaFv5FVVVEYus0g7PqZgD7Qpt3oh3M1qBdN/pnAP7SctvWq2pLRM4DMCMih9EOVJkZPmb/DOBqVX0u2JzrmFWViLwLQBPAa0KbV723qvqD6D0U4t8A3Kiqz4jIH6L9rWaLxefvZRuA3ap6MrSt7GNWmMIDuaq+PucuWgDODf1+TrDtp2h/NTkt6El1thtpm4j8WEReqqqPBkHnsYRdvQ3A11V1KbTvTs/0GRH5MoA/td02VW0F/98vIncAGAOwBzmOm4l2iciLANyK9sn8QGjfuY5ZhLjPTtR9HhaR0wCcgfZnK81ji24bROT1aJ8kX6Oqz3S2x7y3poJSz7ap6k9Dv16P9vhI57Gv7XrsHbbaFbINwAfDGwo+ZmnEtT/3MfMhtfIdAOdLu9JiDdpv0F5tjxLsRzs3DQBXAzDZw98b7DPNvlfl4oJA1slJTwCIHMkuqm0icmYnNSEiZwEYB3B3wcctTbvWAPg62rnC3V23mT5mkZ+dhDZfBWAmOEZ7AWyTdlXLRgDnA/jvnO3pq20iMgbgHwBcpqqPhbZHvreW2/bS0K+XATgS/LwPwBuCNp4J4A1Y+U210HYFbbsA7UHDb4e2FX3M0tgL4PeD6pWLATwRdF7yH7MiR3F7/QPwu2jng54B8GMA+4LtZwP4Ruh+bwJwL9pnz2tC289D+4/rGIB/AXC6wba9GMC3ANwH4HYA64LtTQDXh+63Ae0z6vO6Hj8D4DDaweirAF5os20Afid4/u8G/7+v6OOWsl3vArAE4FDo30VFHbOozw7a6ZrLgp+fHxyDY8ExOS/02GuCxx0F8MYCPv+92nZ78HfROU57e723Ftu2C8BdQRv2A7gg9Nj3BsfzGID32GxX8PsnAUx1Pc7GMbsR7SqsJbTj2vsAvB/A+4PbBcDng7YfRqhSL+8x4xR9IiLP+ZBaISKiBAzkRESeYyAnIvIcAzkRkecYyImIPMdATkTkOQZyIiLP/T9fM9C1mtrtXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x[:,0], y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "547a51fd",
   "metadata": {},
   "source": [
    "# Analytical Model Fitting\n",
    "\n",
    "\n",
    "If the problem is simple enough, we can find a closed form analytical solution to the problem. We seek to find the $\\beta$ that predicts the relationship between x and y. If we use $\\mathcal{L}(\\beta)=||\\mathbf{X}\\beta-\\mathbf{Y}||$, the analytical solution is:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial \\mathcal{L}(\\beta)}{\\partial \\beta} = -2 \\mathbf{Y}^T\\mathbf{X}+2\\beta^T\\mathbf{X}^T\\mathbf{X}\n",
    "$$\n",
    "\n",
    "And solving for 0 gives:\n",
    "\n",
    "$$\n",
    "\\beta = (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{Y}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "61ad3658",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.0753026 1.8672338]\n"
     ]
    }
   ],
   "source": [
    "beta_ols = jnp.linalg.inv(x.T @ x) @ (x.T @ y)\n",
    "print(beta_ols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2430ba06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3oUlEQVR4nO3dd3hU1fbw8e9OI6GGXkIJHUKREoUrigUVBISAomIDC6ACKiFRvHrvT99bQDME4WJBQZGrol6ECCqighXp0oTQEUikQ4CQQNp+/ziTMIRJMsmc6evzPDwkZ86cs+dMsrJn7bX3UVprhBBC+K4gTzdACCGEcySQCyGEj5NALoQQPk4CuRBC+DgJ5EII4eNCPHHSOnXq6OjoaE+cWgghfNaGDRtOaK3rFt/ukUAeHR3N+vXrPXFqIYTwWUqpA/a2S2pFCCF8nARyIYTwcRLIhRDCx3kkR25Pbm4uaWlpXLhwwdNN8Qrh4eE0btyY0NBQTzdFCOHlvCaQp6WlUa1aNaKjo1FKebo5HqW15uTJk6SlpdG8eXNPN0cI4eW8JpBfuHBBgriVUoratWtz/PhxTzdFCJ+WsjGdpGU7+TMjm0aREST2bUtc1yhPN8t0XhPIAQniNuRaCOGclI3pPL9wK9m5+QCkZ2Tz/MKtAHaDuS8Hfa8K5EIIYZakZTuLgnih7Nx8kpbtvCJA2wv6if/bzMtLtpGRlev1gV0CudXJkyfp06cPAEeOHCE4OJi6dY0JVGvXriUsLMyTzRNClNOfGdkOb7cX9HMLNKezcoGye/OeZkogV0pFArOBjoAGHtFarzLj2O5Su3ZtNm3aBMBLL71E1apVSUhIKHo8Ly+PkBD5uyeEr2gUGUG6naDdKDLiim0lBX1bJfXmvYFZkWk68LXW+i6lVBhQ2aTjetTIkSMJDw9n48aN9OrVi+rVq18W4Dt27MgXX3xBdHQ0H3zwATNmzCAnJ4cePXrwxhtvEBwc7OFXIETgSuzb9rJ0CUBEaDCJfdtesW9JQb84RwK+Jzg9IUgpVQPoDcwB0FrnaK0znDymS/5VRFpaGr/++ivJyckl7pOamsonn3zCypUr2bRpE8HBwXz44YcVfflCCBPEdY1i8tBOREVGoICoyAgmD+1kt0ed2LctEaFld7zs9ea9gRk98ubAceA9pdRVwAbgaa31edudlFKjgdEATZs2NeG07jFs2LAye9bLly9nw4YNXH311QBkZ2dTr149dzRPCFGKuK5RDqVCCvcprFqpERHK+Zw8cvMv3dO4pN68NzAjkIcA3YDxWus1SqnpwCTgb7Y7aa3fBt4GiI2NLfWOz950Q+gqVaoUfR0SEkJBQUHR94WzULXWjBgxgsmTJ7u9fUIIcxQP+r5UjmhGIE8D0rTWa6zfL8AI5H4nOjqaL774AoDffvuN/fv3A9CnTx8GDx7MhAkTqFevHqdOneLcuXM0a9bMk80VQjjB0d68N3A6R661PgIcUkoVfuboA2x39rje6M477+TUqVN06NCBmTNn0qZNGwBiYmL45z//yW233Ubnzp259dZbOXz4sIdbK4QIFMqMNIZSqgtG+WEYsA94WGt9uqT9Y2NjdfEbS6SmptK+fXun2+JP5JoIIWwppTZorWOLbzel/FBrvQm44uBCCCFcT9YjF0IIHydTFYUQwsVcXQEjgVwIIVyovKswVoSkVoQQwoVKW4XRLNIjF0L4DW+cxFOeVRgrSnrkNtLS0hg8eDCtW7emZcuWPP300+Tk5ADwww8/MHDgwCue88UXX9C1a1euuuoqYmJimDVr1hX7XLx4kVtuuYUuXbrwySef8Nhjj7F9u1Fq/+9//9u1L0qIAFGYwkjPyEZzKYWRsjH9iv16TVlB9KQvafn8V0RP+pJeU1ZcsZ9ZSlqfxcx1WySQW2mtGTp0KHFxcezevZtdu3aRmZnJCy+8UOJzcnNzGT16NEuWLGHz5s1s3LiRG2+88Yr9Nm7cCMCmTZu45557mD17NjExMYAEciHM4kgKwzbYA+Rb59GUFPTNYG9BLrPXbZFAbrVixQrCw8N5+OGHAQgODmbatGm8++67ZGVl2X3OuXPnyMvLo3bt2gBUqlSJtm0vf3OOHTvGAw88wLp16+jSpQt79+7lxhtvZP369UyaNIns7Gy6dOnC/fff79oXKISfcySFYS/YFzI7b12oPKswVpTkyK22bdtG9+7dL9tWvXp1mjZtyp49e+w+p1atWgwaNIhmzZrRp08fBg4cyPDhwwkKuvT3sV69esyePRuLxVK0TkuhKVOmMHPmzKIbWgghKs6RG0mUlZd21Xrjrl63xXsD+UcuuPnwfeavqjh79my2bt3Kd999h8Vi4dtvv2Xu3Lmmn0cIUTpHbiRR1g0kvHW98bJ4byB3QdAtTUxMDAsWLLhs29mzZzl48CCtWrVi7dq1JT63U6dOdOrUiQcffJDmzZtLIBfCA4qvKW6vasVesC/kyvXGZUKQm/Tp04dJkyYxb948HnroIfLz85k4cSIjR46kcmX7d67LzMxk/fr1RQOcmzZtKvfStaGhoeTm5hIaGursSxAi4JWVwrAN9ukZ2QQrRb7WRNkJrvaCb+FzC28+oRRkZOWWGpzdMSFIArmVUopFixbx5JNP8o9//IOCggL69+9/WVXJ8uXLady4cdH38+fP59VXX2XMmDFERERQpUqVcvfGR48eTefOnenWrZvcHk4IN3AkX20v+CYu2AwacguMbEFGdm7R/iUF55SN6Uz8dHNRdUwhs2/kbMoytuUly9g6Rq6JEI4zM33Ra8oKh27GXFxUZAQrJ90MwIspW/lw9UFKirAK2D9lQLmO79JlbIUQwpPMTl9UtHql8HkpG9NLDeIgE4KEEOIyZq9nUtEgW/i8pGU7Sw3ifj0hyJtuuuxpci2EcJzZ65nYm41ZFtvgbO+8OccPcOKr18g7us9/JwSFh4dz8uRJateujVIuqCH3IVprTp48SXh4uKebIoRPcGQyUHnYK2U8fzHvsgFOW8WrXgrbo7Xm4sGtnF27kOx9xrhg98ZVies6vkLtKonXBPLGjRuTlpbG8ePHPd0UrxAeHn5ZhYwQomSOTAYqr+LVLcXz8IXnsNe7ju/Tkqcmv8WJVQvIOWLMDO/bJZRn76lP9J3/rHCbSuI1gTw0NJTmzZt7uhlCCB/kyGQgd5wjMzOTOXPmMG3aNP48cIDgIBh+QwR/vTOCqPq1qHn136FZG9PaVMhryg+FEMKdzCxXPHz4MP/5z3948803ycjIIDwUnhtWj/jbC6hStw3BHf8KjfqDk2ljKT8UQggrs8oVt2/fztSpU/nggw/IyckhsjK8Nb4pI3qepVKja1AdJkHdXi55DbYkkAsh/FJpPe7SyhXLCuRaa3788UcsFgtffvklAFG14PXxrenf7iihzW6E9s9CZAeXvC57JJALIfxGYfBOz8hGQVEtd/Eed0XKFfPy8li4cCFJSUkUpoY7R4fxxvgW9Iw6THCrAdBuAlRpauZLcogEciGEXyieLik++mfb4y5PuWJmZibvvfce06ZNY//+/QDc1r0Gr41qSNvaJwhqOxzajIVKtU1/TY6SQC6E8Aul3f2nUGGP25FyxSNHjjBz5kzeeOMNTp8+DcDD/Rryz+FVaVDtIkExT0LLRyCkit1zufNG0BLIhRA+L2VjukOLXBX2uEsrJUxNTSU5OZl58+aRk5NDcBD8/aHWTOiXS43qVVExz0GzeyCo5KWn3bF0rS0J5EIIn1YYNMtSvMdtO+FHa83PP//MoEFPsGTJEgAqV4L/JFzFQ1cfJ7xmfYh5Hhrd7lAJoTODqRUhgVwI4dNKS6kUDnjau3EEGAOYixYtIikpiXXr1gHQoFYlZk3sTP9W+wmp3wxi3oC615arTWav/VIWCeRCCJ9WWnCcdk8Xuz3g8+fP895775GcnFw0gNmxZSSznmlHz3o7CGraAdrPhRoxFWqT2Wu/lEUCuRDCa1RkgLCkoBkVGXHFc48ePVo0gHnq1CkAbuvRhOmPR9G2yk5Ui79Au0+hShOnXocr1n4pjQRyIYRXqOgAoSNBc8eOHUUDmBcvXgTg0cEd+H/3VqZh6AFUm/7Q5kuoVMuU1+KOtV9sSSAXQniFig4QlhQ0B3dpxM8//0xSUlLRAKZS8PKYnjx96wVqhGRAu8etJYT2b7DuDEfuDWoWCeRCCK/gzAChbdDMz89n0aJF9Hw8ibVr1wJQOSKMGROv4/5u6YRXOg8xk6Dp3RDkHyHQP16FEKJU7pycUlHODhAWDmBOmzaNffv2AdC4QS3eevYa+jXbRnD1PIiZBg37Ob0KobeRQC6En3P35JSKqugAob0BzG4do3n9mY70iFyDqlsJ2n8Cdf/i0vZ7kgRyIfycuyenVFR5Bwh37tzJ1KlTLxvAvKNPF6aOakyr4JWoxnWg/Y9Qo73bXoOnmBbIlVLBwHogXWs90KzjCiGc4+7JKc4oa4BQa83KlStJSkpi8eLFRdufvP8m/j4snHq5q1EtboZ2b0LlwLlVopk98qeBVKC6iccUQjjJ3ZNTXCE/P5+UlBSSkpJYs2YNAJUqVeKlcX0Ze1Mm1S7+Di2fgjYfQljNCp/HF8YS7DElkCulGgMDgH8B8WYcU3gnX/1BD2TunpxipqysLObOnUtycjJ79+4FoFatmkxL7Mu9nf8gLHcLtE6AFl9AiHN/mHxlLMEes3rkrwHPAtVMOp7wQr78gx4oyrorjq/8AT527BjjX5zMwg/fJS/rLACNmjThvRdv4pqaazh+fi0JW+5kE7cS3zSGOCeDOPjOWII9TgdypdRA4JjWeoNS6sZS9hsNjAZo2tT9d9AQzvPlH/RAUNYfWl94j3bt2sXUqVN5b+775OYYA5g1mrbi8WHNGd9xMxfDU0nc+wDfZnTBWBIrx7TOhC+NJRRnRo+8FzBIKdUfCAeqK6U+0Fo/YLuT1vpt4G2A2NjY4jfvED7Al3/QA4Ev/6G1HcDU2ggPUZ268dSQWjzcYj2/ZdVh3MFn2ZLdnnx9efgw6zW6YizBXalIpwO51vp54HkAa488oXgQF/7BHwbN/Jmv/aHNz8/n888/x2KxsGrVKsAYwHz6sSG0bnKUuxquY9nZv3DvvsnsvVi4iJX9PmB6Rja9pqxwKmCaPZbgzlSk1JELh/nyoFkg8NY/tMV7pU/1bsqpzd+SnJzMnj17AKhZsyYvT7ibx3qdIuLUN3x0og99d83kaF6dy44VrNQVPXIwkiyFr91ewHSkZ2z2WII7PyGZGsi11j8AP5h5TOE9fHHQLJB42x/alI3pvLxkG6ezcgHIzzrDtl8+5L7/9yUF2cYAZr1GTZj54p0MabuLkIwUqDeeL6v/H0nb0zmdl3vZ8SJCg7mzexSfbUi/7DUW3jzClm3ALE/P2MyxBHd+QpIeuSgXXxk0C0RxXaNYf+AU89ccIl9rgpXizu6eeb9sg2fuqXTOrkvh/O/L0Xk5AIQ1aMWdg7sQ33077SMXEtJkEvReQMrWUzz/+dYrerKREaG8NKgDcV2jiG1W67LOREn36iwMmJ4aO3DnJyQJ5EL4iZSN6Xy2Ib0o9ZCvNZ9tSCe2WS23B/OkZTvJ+ON3zqxdSPau1RT2mau17s5DcS0Z13YNBWzgzWN3sfnkrfw07Nai59m7bVuVSiFFr6F4Z6LXlBWlBkxPjR248xOSBHIh/IQ3VK3k5+ezePFiNsx8gYvpqcbG4BDqdenN6IG1GRP9A4dycph8+BF+zOwGKBQ5Rc+vSNAtK2B6auzAnalICeRC+AlPVq1kZ2fz/vvvk5yczO7duwEICq9K07/0YXzfIB5quIJ15zsw9sAkNmdf3iO1DagVCbplBUxPjh24KxUpgVwIP+GJnueJEyd4/fXXmTlzJidOnAAgOjqaYfcNoWO9XQyu+R1Lz/bi7r2vsO/ilYtYFQ+oFQ26pQXMQBikl0AuhJ9wZ89zz549JCcnM3fuXLKzjT8e3bt3558J93Jbk00EHZ7L7qp389DWOWw5WZVGkRG8Fme0o7SA6qqg6++D9Erbqcl0tdjYWL1+/Xq3n1cIf+fqmYSrVq3CYrGwaNGiohmYAwYM4B/j+9Kl0jeoU+ug7VPQ+gmnViEU9imlNmitY4tvlx65EH7EFT3PgoICFi9ejMViYeXKlQCEhYXx4AP38/fHutM0cz5kvwYtEuG6Tx1ehVBW0jSPBHIhhF3Z2dnMmzePqVOnFg1gRkZGMu7JMUy8qxGRf74NGZsg5jlocme5bmQsK2maSwK5EOIyJ06c4I033mDmzJkcP34cgGbNmpE44Ukeuwkq7ZsJp1pDt2RocGuFbmTsDaWS/kQCuRACgL179zJt2jTefffdogHMbt268ULCGOLapxG01wJnesP1n0Htq506l68t8OXtJJALEeDWrFlDUlISCxcuLBrAvP3223lxwoP8peavqD8mwcU74dZfoHobU87prQt8+aogTzdACOF+hQOY119/PT179uSzzz4jJCSEhx9+mF1rP+OrF2tz7ZlxqODK0P936PGOaUEcjFLJiNDgy7bJSpoVJz1yIQLIhQsXigYwd+3aBUCNGjV44onHmfhgD+ocmwOHxkLbpyH2PxAW6ZJ2BMIkHXeSQC5EADh58mTRAOaxY8cA45aLEyY8zegBUVTePwP2LICYZ+H6BRAc7vI2+fskHXeSQC68mj/XGrvjtdkbwOzatSvPJUzgrmtyCd45FfaE25QQBpdxRHP48/vqCRLIhdfy51pjV7+2NWvWYLFYWLhwIQUFBYAxgPlc/Fh6R+1C7XwBDrSB7tOhfp8KlRBWlD+/r54ig53Ca5VWa+zrXPHaCgoKWLJkCb1796Znz54sWLCA4OBgRowYwfaNP/LVlFhuODMSdXIVXL8Q+nwHDW5xaxAH/35fPUV65MJr+XOtsZmv7cKFC3zwwQdMnTqVHTt2AMYA5uOPP86EUUOof/oD2D0Ymt0Nt/4K1Vs71XZnmf2+SppGeuTCi5VUU+wPtcZmvLaTJ0/yr3/9i+joaEaNGsWOHTto0qQJycnJpG9bypRBadTf2B9CqsDA7XDNLI8HcTD3fS1M06RnZKMx0jTPfLKJLi9/Q8rGdCdb6jskkAuv5c+1xs68tn379jF+/HiaNm3Kiy++yNGjR+nSpQsffPBf9q2ey4Su31Fl7Z0Q2QkG7YMuUyCioateSrnd1K5uubaXpqRbw2Vk5/L8wq0BE8wltSK8lj/XGlfkta1bt46kpCQ+++yzogHMfv36kTAxnpvbnkelvgrrj0NMojGNvowSQk+lJL7fcbxc20tTWjomkNZukUAuvJo/1xo78toKCgr46quvSEpK4qeffgIgNDSUBx98kIQJT9Gx2hbY/jRsq2yUEDYe6lAJoScrR8zMkZc01d+ZY/oiSa0I4YUuXLjAnDlz6NixI3fccQc//fQT1atX59lnn+WPPb8zd9JVdNw/GP740JiB2XcdNB3mcB24JytHzMyR20tROXtMXyQ9ciG8yKlTp3jrrbeYMWMGR48eBaBJkyY888wzjHowjmp/vgfrekH9m6B3CtTqXqHzeLIiyMxb0hV+enh5yTZOZ+Ve9pi/jKc4QgK5EF5g//79TJs2jTlz5pCVlQXAVVddRWJiIncPuJrQ3TPgh+6llhCWJ+ftydUHzR77KExRBXIZotyzUwgPWrduHRaLhQULFhQNYPbt25eJEydyS/e6qNQkOPw1tBptLGQV0cDucYrnvMHokU4e2sluMCvv/sI7yD07hfASBQUFLF26lKSkJH788UcAQkJCeOCBB5gYH0/nBhmwfQr8uBnaPgNXvwFhNUo9ZnnvuOPPFUGBSAK5EBVU3o/yFy9e5MMPP8RisZCamgpA9erVGTNmDE+NH0djNsD2MXDwFLRPNHLgwZUcaktFct6OVgQFcsrCV0ggF6ICylO+d/r06aIBzCNHjgDQuHFjYwDz0RFUP7kENveFkKoQMwkax5V7FcKSct6RlUMr8OoukQWufIOUHwpRAY6U7/3xxx8888wzNGnShL/+9a8cOXKEzp0789///pd9OzczsT9U/6ErHPgYYl+HvmuhacWWkk3s25bQ4CsXv8q8kOfU7EZZ4Mo3SCAXogJKS2Vs2LCB4cOH06pVK6ZPn8758+e55ZZbWLZsGZtWf8MDnXYQurQtnFwHvT+Hm5dBg5udWoUwrmsUVcKu/ICdW6CdCrr+vHCZP5HUihAVUDyVobXmwr71XPjtc2Jf2QTYDGBOnEiXVtUh1QJf3gtN74HbVkO1lqa26Ux2rt3tzgRduUmyb5AeuRAVUDijUOflkrnlWw7PGcuxBS9zdt8mqlWrxsSJE9m3bx//nTGRLudfgWXXQGgkDEiFa940PYiDa1aL9OeFy/yJ9MhFwDGjCuOG6MrEnvmRT+a+Q+65kwDUrteA5xLiGT1qFDUuboLtoyBjK7SbYCwhG1rdBa/mEjNnTBaSMkXfIIFcBBRnqzAOHDjA+L/+ky8XfEBBzgUAmrVqxz/+/lfuuXsYYceWwpq+kJNh3Mi49+cOlxA6y1VB158XLvMXEsgDgKfrgD19flvlnThT6LfffiMpKYlP//c/CvKN54c360L1HkOp1uYqOjTeTdi3XSCkOnSYBFGDS6w+ceX1kKAbmCSQ+zlP1wG78/yOBMjyVGForfn666+xWCysWLECABUUTJUON1H96iHUbNiA+2p9zSN1pnNgXwsY+CbUu7HU6hNPvx/CPzkdyJVSTYB5QH1AA29rrac7e1xhjor2QH3t/I4GSEeqMHJycvjoo4+YOnUqv//+OwBVq1Zl9OjRfJrdiXo1QxlZZwn3117KL5ldeOyPv7P9Qkv217+pzHZ6+v2oCG/6RCXsM6NHngdM1Fr/ppSqBmxQSn2rtd5uwrGFkzxdB+yu8zsaIEsbEMzIyGDWrFnMmDGDP//8E4BGjRrx9NNPM3r0aCJDTtHrwwn0iVjOkjO9idszlYM5xi3UohysDPH0+1Fe8gnCNzgdyLXWh4HD1q/PKaVSgShAArkX8HQdsLvO72iAtDcgOKJTFX7671QefOcdMjMzAejYsSMJCQkMHz6csPPbYdsTcORbOkTfx4BV75B24VIFSnkqQyIrh16xbnbhdm/ki58gApGpdeRKqWigK7DGzmOjlVLrlVLrjx8v/735RMV4ug7YXecvTw11XNcoVk66mYX3NCJ66xyeHNyLadOmkZmZSZ8+fVi6dClbNm9mRL+mhP1yB/ww0LiBw6B9tO07g4TBvYmKjEBh9MTLs/RrSatGu3I16ZSN6fSasoLmk76k15QV5Zqy72ufIAKVaYOdSqmqwGfAM1rrs8Uf11q/DbwNxnrkZp1XlM7TdcDuOr+jNdRaa5YtW4bFYmH58uUABAcHc99995GQkEDXqzpD+ufwTU/IOwvtn4Xo+y8rIXSmMqSk2ZclbXeWs6kRT3+iE44xJZArpUIxgviHWuuFZhxTmMfTJWnuOH9ZfzBycnKYP38+FovlsgHMUaNG8cwzz9A0qj7s/68xhT4sEjo8D40HgzJ38rO7A6OzqRGzJhnJgKlrmVG1ooA5QKrWOtn5JglRMfb+YGRkZPD2228zffr0ogHMhg0b8vTTTzNmzBgiqwTB7lmw+DWI7GTMwKx3g1MLWJXGFbMvS+NsasSMT1QyYOp6ZvTIewEPAluVUpus2/6qtf7KhGMLUSEHDx5k+vTpvPPOO5w7dw6ADh06kJCQwH333UdY/inY+QrsfQca3Ao3fgk1u7i8Xe5OdZnxCcDZT1QyYOp6ZlSt/AK4pvsiAo6zH8E3bdqExWLh448/Jt86A/Pmm28mMTGRvn37ojL3wqan4OCn0Gy4sQZ41Rauejl2uTPV5e5PAPbIgKnryczOEpQnoEj+zxz2PoInLtjMS4u3cSY7t8Rrq7Xmm2++ISkp6bIBzOHDh5OQkEC3bt3g1G+w8l44uhxaPQEDd0B4Pbe/Rnfz9GA3yICpO0ggt6M8OT3J/5nH3kfw3HxNhrWio/i1zcnJ4eOPP8ZisbB1q7G9SpUqRQOYzZo2haMrYMVtcGY7tIuHHrMhtJp7X5iLONqB8PRgtzd8KvB3EsjtKE9OT/J/5nHko3Z2bj6TP9/AnuXzee2110hPN2qiGzZsyFNPPcWYMWOoWaM6pC2CZXdBXqZNCWGYq1+C2/hSB8IbPhX4OwnkdpQnpyf5P/OU9BG8UN7Z45zbsISDm5ayNsfYLyYmpmgAs1IIsH8e/JIEYbWgwwvQeJDpJYTewNc6EJ7+VODvJJDbUZ6cnuT/zGPvIzhAzrF9nF27iPOpP0GB8dhNN91EQkIC/fr1Iyg/E3ZPh52vQWQXuOYdqNfbZSWEnpayMb3EP3jSgQhM/tdVMUF5ppV7egq8P4nrGsXkoZ2MBai0JuTwFo598jcOv/cU57d9D1pTLeYGkj74khUrVtD/pm4EbXkBFreA05vhxqVw01dQ33V14J5WmFIpiXQgApP0yO0oT05P8n/mGtCxHud+X4FlvoUtW7YAEBQWTpVOt9G6zz387d4biGuVDWsft5YQ3gd910HV5h5uuXvYS6nYuqldXTe2RngLCeQlKE9OT/J/zjt79mzRDMy0tDQAGjRowPjx43niiSeoWbMmnNoA2yfAN99Dq8cDpoTQVlmpk+93yIJ0gUgCufBoHXxaWhrTp0/n7bff5uxZY6219u3bk5CQwP3330+lsDBrCeEUOLvDWkI4x29KCMurrAFhyZEHJgnkAc5TZWxbtmzBYrEwf/588vLyALjhhhtITEzk9ttvJwgNaQth+yuQl2XcyLjZfX5VQlgRJQ0IF5IceWCSQO4h3jIb9OUl29xWxqa1Zvny5SQlJfHNN98AEBQUxN13301CQgJXX3015F+AfbNhexKE14WOf4OoOypUQugt19hMhe1/afG2oolShWSQPXBJIPcAb5nMkbIx3e7dasDcj+i5ubl8+umnWCwWNm3aBEDlypV59NFHmTBhAs2bN4ecM0bve+d0qNkVer4Lda+rcPWJt1xjVygck/HHP1SiYiSQe4C3TOZIWrazxMfK+ojuSBA5e/Yss2fP5rXXXuPQoUMA1K9fv2gAs1atWpB9GDZNgr2zoWE/uPFrqNnZ4ddQUju85Rq7kgyyi0ISyD3AzNmgzvTKSjtfaR/Ry+rtpqenM336dGbNmlU0gNmuXbuiAczw8HA4uxvWTIJDC4zp833XQ9VoB1912e2QGbcikEgg9wCzZoO66jZekRGhpT6/pN7uS+8vZdFrP/PRRx9dNoCZkJBA//79CQoKgpPrYd0rcOwHaP0kDNxp5MIroLRet8y4FYFEZnZ6gFmzQUsLZM6046VBHUp9nm2vVmtN9h+bOPrp39k8fRTz5s2joKCAYcOGsWbNGn744QcGDhhA0NHlsPwW+Hko1O0Fg/ZD55crHMSLt6P4dplxKwKJ9Mg9wKzZoJ66jVejyAjSTp4ja+cvnF27iJyjewEICg1n7OPGErItWrQw1kU58KkxiFlwwViFsNlwCA6zpoTWOvX6S+t1mz3jVgYWhTeTQO4hZgxUeeI2XufOnaPt8Z9YO/tN8s4eAyCociS1rhnE1L8n8tBNHY0Swt2zIDUJwutDp/+DqIFFJYRmVZSUtc61WYOB/lwBI/yDBHIf5s4F+9PT05kxYwazZs3izJkzxrnqNqFy9zhaX9uf5wZ2Iq5DVdg2xSghrNUdes6FetddcSyzKkrctc5NIFTACN8mgdyHuSOQ/f7771gsFj766CNyc42a8969e5OQkMCAAQOMAczsw7DjNVg8Gxr1h5u/Me5IXwIzK0rcUYInFTDC20kgF1fQWvP999+TlJTE119/DRgzMIcNG8bEiRPp0aOHsePZ3Ub65NACiH4A+m1wqITQ1ypKvL29kr8XUrXiwwpzt+kZ2Wgu5W5TNqZX6Hh5eXnMnz+f2NhY+vTpw9dff01QaDjVug2gS/z73Pf8NCOIn1wHP98F314LEQ1h4C6IneFwHbivVZR4c3vN/hkQvkl65D7MrNztuXPnmDNnDtOmTePgwYMA1KhVh9BO/Qm/qh/BEdU5iWbJsrn0OrCUuhyEdhONHHhoVWuPcIXDPUKzUkLu6ol685rzkr8XIIHcpzmbu/3zzz/5z3/+w1tvvUVGRgYAbdq0ISEhgXePNOHw+XyCyef2Gj/xeL3PCFW5vHVoOH8b+w8ICgUqXtHhbG7b3ZUk3jodXvL3AiS14tNKytGWlbvdtm0bjzzyCNHR0UyZMoWMjAyuu+46UlJSSE1NZdSoUZzOyuL+Wl+xou0YRtT5guQj99Nv10zeTb++KIiD85OSKspT5/U2Ff0ZEP5FArkPK0/utnAAs3///nTs2JH33nuPvLw87rzzTlatWsXPP//M4MGDCco7C9v+zS8xo7ip+jomHprAsL2vsuLcNWiCrggQnuoRSk/U4M35e+E+klrxYY7kbvPy8liwYAEWi4UNGzYAEBERwSOPPMKECRNo2bKlsWPWn8Zd6PfOgUYD2NzqU8Z/lVtmjborKzpKy4F7eyWJu3hz/l64j9Jau/2ksbGxev369W4/byDJzMwsGsA8cOAAAHXr1mXcuHE8+eST1KlTx9jx7E5rCeFCiH4Q2sdDlWaAY4OJxXPVYAT8yUM7mZoDL35cV51XCG+mlNqgtY4tvl165F6uvJUZhw8fZubMmbz55pucPn0agNatWzNx4kQeeughIiKsPdYTayH1FTj2M7QZC3fshkq1LzuWIwN8ruoRllWNIT1RIS6RQO7FHKnMKAz0f+zZSd7mxZzavIK83BwArr32WhITExk0aJAxA1Nr+HOZEcDP7YX2E+Ev8yCkilPtdEVFhyM5cG+tJBHC3QIukJe3h1vS/u6oYS6rV7rotzSeee0jjv+6gOy966x7KHrefDtT//Ei1157rbGpIA/++NQI4AV51hsZ33tZ9Ym3kRy4EI4LqEBe3trjkvZff+AUn21Id3kNc0m90vRTmXzyySc8OvH/OJ9ulNupkDCqdLqF6rGDCWrRygjiedmwf65xI+PKUdD5n8ZaKBW8D6Y7uXNBMCF8XUAF8vLOgitp//lrDpFfbJDYFbPpivdKC3KyydzyLVm/LebepCMABEVUp1q3gVTrNoDgyjUAyDx3HLb9G3bOgNo94NoPoO61prXLHSQHLoTjAiqQl7f2uKTtxYN4WftXVGGvNPP0Cc7+toTMjV9RcCETgFatWlHQcSC5La4nKLQSAPVDTvBo3c+5p9Z3cHYI3LwcIku/2483kxy4EI4JqAlB5Z0FV9L24BJSE2bnb9uGnyVq61zS33qEs6s+peBCJu2uimXhwoXs2LGDqX9PoErlyrSsdIhXGk9nWZtxVArSrG77HfzlfZ8O4kIIxwVUIC/vLLiS9h/eo4nLZtNprfnxxx+54447iImJYfnnn0BBHkOGDGHlypWkblrHkCFDCA4OJq7JIZZf/RoLWj1Pek49hh+ZR43rZtL3mmucbocQwncEVGqlvHnX0vaPbVbL1PxtXl4eCxcuxGKxsG6dUYESHh7OyJEjiY+Pp3Xr1saOWsPhr437YGbup1H7idDyM+JDqhBf4bMLIXyZKTM7lVL9gOlAMDBbaz2ltP1lZucl58+f591332XatGns378fgNq1azNu3DjGjh1L3brWu8wX5MHBT2H7q6DzIeY5aHaPV5cQCiHM5bKZnUqpYOB14FYgDVinlFqstd7u7LH92ZEjR5g5cyZvvPFG0QzMVq1aER8fz4gRI6hcubKxY14W7HsPUqdC5cZw1b+h0e0+UUIohHAPM1Ir1wB7tNb7AJRSHwODAQnkdqSmppKcnMy8efPIyTFmYPbs2ZPExEQGDx5McLA1955zGna9Drtm+mwJoTvJ7c5EIDMjkEcBh2y+TwN6FN9JKTUaGA3QtGlTE05bNm/55dZa8/PPP2OxWFiyZAkASini4uJITEy8NAMTICsNdkwzeuGNB0OfFVAjxu1t9iXuvsmEEN7GbYOdWuu3gbfByJG7+nze8Mudn59fNIC5du1awBjAHDFiBPHx8bRp0+bSzmdSjVUI01Kg+Ui4fTNUaeKWdvo6ud2ZCHRmBPJ0wDbiNLZu8yhP/nKfP3+euXPnkpyczL59+wBjAHPs2LGMHTuWevXqXdr5xGrYPgVOrILW4+COPVCplkvb52/kJhMi0JkRyNcBrZVSzTEC+L3AfSYc1yme+OU+evQor7/+Oq+//jqnTp0CoGXLlsTHxzNy5MhLA5haw59LjUWszltvZHztRxBS2WVt82eywJYIdE4Hcq11nlJqHLAMo/zwXa31Nqdb5iR3/nLv3LmT5ORk3n//fS5evAhAjx49SExMJC4u7tIAZkEeHPgEUl81vo95DpreDUGuyXB5yxiBq8kCWyLQmRJBtNZfAV+ZcSyzuPqXW2vNypUrSUpKYvHixYAxgDl48GASEhLo1asXqrBEMC8L9r4LO6Yad9/pMgUa9nNpCaE3jBG4iyywJQKd387sdNUvd35+PikpKVgsFlavXg1ApUqVigYw27a1+UNx8ZRRQrh7JtS5FnrNhzo9nTq/owJtANDeAluB8olECL8N5GDu6nlZWVlFA5h79+4FoFatWkUDmPXr17+08/lDRgnh/rnQOA76/AA12pvSDkcF+gBgIH0iEcKvA7kZjh07VjSAefLkSQBatGhRNIBZpYrNbdLObDem0KcvhhYPQ/8tpOxUJL25kz8z9rm1V+jNA4DecHclIfyJBPIS7Nq1q2gA88KFCwBcffXVJCYmMnTo0EsDmADHVxklhCdXQ5un+DL6Z/79zVHSP92MAgqL5t3ZK/TWAUB39ZQD/ROJCCwBtYytI1auXMmQIUNo164ds2bN4sKFCwwaNIiffvqJNWvWMGzYMCOIaw3pX8K3veHX+6BhXxj0Bym5I0n4/FBRb7j4zKfCXqGrxXWNYvLQTkRFRqCAqMgIJg/t5PHeaGk9ZTOVd415IXyZ9MgxBjA///xzLBYLq1atAowBzIceeoj4+HjatWt3aeeCXJsSwiBrCeGwohJCe4GqOHf1Cr3xDjvu6il76ycSIVwhoAN5VlYW77//PsnJyezZswcwBjCfeOIJxo8ff/kAZl4W7J1jLSFsDl1eNXrhxUoIHQlIgdwrdFfuXkoSRSAJyEB+/PjxogHMEydOANC8eXPi4+N5+OGHLx/AvHjSWkL4urWE8BOoc8WaYEVKClSFAr1X6M6esjd+IhHCFQIqkO/evZvk5GTmzp17xQDmkCFDCAmxuRznD8GOZNj/PjQeAn1+hBrtSjjyJfYCVeGAZ5T0CqWnLIQLBEQg//XXX7FYLKSkpFB4R6Q77riDhIQErr/++kszMKFYCeEj0H8rVHY8yEigKpv0lIUwl98G8vz8fBYvXkxSUlLRAGZYWFjRAGb79sUm6Bz/1VpCuBbajIdBeyGsZqnnKKke2pcClcx+FML3+V0gz87OLhrA3L17NwA1a9bkySefZNy4cTRo0ODSzroA/vzKuJFxVjq0TzBy4CFlD7z5w8xBf3gNQgg/CuT2BjCjo6OLBjCrVq16aeeCXDjwsZFCUcEQMwma3lWuVQj9YeagP7wGIYQfBHJ7A5ixsbFFMzAvG8DMO2+UEKZOhWotoasFGt5WoVUIK1oP7U2pDJn9KIR/8NlAvmrVKpKSki4bwBw4cCAJCQn07t378gHMiyeNmxjveh3qXQ/X/Q/qXOPU+StSD+1tqQxvXo9FCOE4n5qiX1BQQEpKCr169eLaa69l0aJFhIaG8uijj7Jt2zaWLFnCDTfccCmInz8IG56BJa0h6xDc+jNc/5nTQRyMMsOI0ODLtpVVD+2u6emOqshrEEJ4H5/pkWdnZ9OtWzd27NgBQGRkZNEAZsOGDYv2S9mYzoIVXzM0/EP6VF/Psbr30bqcJYSOqEiZobelMqRUUgj/4DOBPCIigvbt25OdnU18fDyPPPLI5QOYwE8rF1Hj9ykk19nN3BN38FL6GHJ31WByfYjran6byltm6I2pDF8qlRRC2OczgRxg1qxZ1KxZ8/IBTF1grEKY+got/9zHG2eG8Pgfz3FRVzIeL/CeKgxZyEkI4Qo+Fcjr1q176ZuCXPjjI2MVwqBKEPMcNyyPII/gK57nLVUYksoQQriCTwVywCgh3DPbWIWwWmvo9ho0uAWUon7kCq9LXRQnqQwhhNl8qmqFoz/A583h+C9G9Umf5dDw1qI6cKnCEEIEIt/qkUd2hlt/gept7D4sqQshRCDyrUBeqZbxrxSSuhBCBBrfSq0IIYS4ggRyIYTwcRLIhRDCx0kgF0IIHyeBXAghfJwEciGE8HESyIUQwsdJIBdCCB/nWxOC7PCmW6cJIYQn+HQg97ZbpwkhhCf4dGrF226dJoQQnuDTgdzbbp0mhBCe4NOpFW+8dZo9kscXQriSUz1ypVSSUmqHUmqLUmqRUirSpHY5xBfWHy/M46dnZKO5lMdP2Zju6aYJIfyEs6mVb4GOWuvOwC7geeeb5Li4rlFMHtqJqMgIFBAVGcHkoZ28qrcreXwhhKs5lVrRWn9j8+1q4C7nmlN+3r7+uOTxhRCuZuZg5yPA0pIeVEqNVkqtV0qtP378uImn9W4l5eu9LY8vhPBdZQZypdR3Sqnf7fwbbLPPC0Ae8GFJx9Fav621jtVax9atW9ec1jspZWM6vaasoPmkL+k1ZYVL8ta+kMcXQvi2MlMrWutbSntcKTUSGAj00Vprk9rlcu6aTCT3ERVCuJpTOXKlVD/gWeAGrXWWOU1yj9IGIc0Ost6exxdC+DZn68hnApWAb5VSAKu11o873SonOFqzLYOQQgh/4WzVSiuzGmKG8qRLfGUykRBClMWnp+gXV56abXuDkAoj+Ltq4FMIIVzBp6foF1eedIntIGR6RjYKKByplVUUhRC+xK965OWt2Y7rGsXKSTcTFRlB8XIbmX0phPAVfhXIK1qzLQOfQghf5leBvKJrr8jsSyGEL/OrHDlUrGY7sW/by6pdQGZfCiF8h98F8oqQ2ZdCCF8mgdxKZl8KIXyVX+XIhRAiEEkgF0IIHyeBXAghfJwEciGE8HESyIUQwscpT9wLQil1HDhQwafXAU6Y2ByzSLvKR9pVPtKu8vHWdoFzbWumtb7iFmseCeTOUEqt11rHerodxUm7ykfaVT7SrvLx1naBa9omqRUhhPBxEsiFEMLH+WIgf9vTDSiBtKt8pF3lI+0qH29tF7igbT6XIxdCCHE5X+yRCyGEsCGBXAghfJxXBnKl1DCl1DalVIFSqsQyHaVUP6XUTqXUHqXUJJvtzZVSa6zbP1FKhZnUrlpKqW+VUrut/9e0s89NSqlNNv8uKKXirI/NVUrtt3msi7vaZd0v3+bci222e/J6dVFKrbK+31uUUvfYPGbq9Srp58Xm8UrW17/Hej2ibR573rp9p1KqrzPtqEC74pVS263XZ7lSqpnNY3bfUze1a6RS6rjN+R+zeWyE9X3frZQa4eZ2TbNp0y6lVIbNY668Xu8qpY4ppX4v4XGllJphbfcWpVQ3m8ecu15aa6/7B7QH2gI/ALEl7BMM7AVaAGHAZiDG+tinwL3Wr98CnjCpXa8Ck6xfTwJeKWP/WsApoLL1+7nAXS64Xg61C8gsYbvHrhfQBmht/boRcBiINPt6lfbzYrPPk8Bb1q/vBT6xfh1j3b8S0Nx6nGA3tusmm5+hJwrbVdp76qZ2jQRm2nluLWCf9f+a1q9ruqtdxfYfD7zr6utlPXZvoBvwewmP9weWAgroCawx63p5ZY9ca52qtS7rzsfXAHu01vu01jnAx8BgpZQCbgYWWPd7H4gzqWmDrcdz9Lh3AUu11lkmnb8k5W1XEU9fL631Lq31buvXfwLHgCtmrpnA7s9LKe1dAPSxXp/BwMda64ta6/3AHuvx3NIurfX3Nj9Dq4HGJp3bqXaVoi/wrdb6lNb6NPAt0M9D7RoOzDfp3KXSWv+E0XEryWBgnjasBiKVUg0x4Xp5ZSB3UBRwyOb7NOu22kCG1jqv2HYz1NdaH7Z+fQSoX8b+93LlD9G/rB+rpimlKrm5XeFKqfVKqdWF6R686Hoppa7B6GXttdls1vUq6efF7j7W63EG4/o48lxXtsvWoxi9ukL23lN3tutO6/uzQCnVpJzPdWW7sKagmgMrbDa76no5oqS2O329PHaHIKXUd0ADOw+9oLX+3N3tKVRau2y/0VprpVSJtZvWv7SdgGU2m5/HCGhhGLWkzwH/z43taqa1TldKtQBWKKW2YgSrCjP5ev0XGKG1LrBurvD18kdKqQeAWOAGm81XvKda6732j2C6JcB8rfVFpdQYjE8zN7vp3I64F1igtc632ebJ6+UyHgvkWutbnDxEOtDE5vvG1m0nMT6yhFh7VYXbnW6XUuqoUqqh1vqwNfAcK+VQdwOLtNa5Nscu7J1eVEq9ByS4s11a63Tr//uUUj8AXYHP8PD1UkpVB77E+CO+2ubYFb5edpT082JvnzSlVAhQA+PnyZHnurJdKKVuwfjjeIPW+mLh9hLeUzMCU5nt0lqftPl2NsaYSOFzbyz23B9MaJND7bJxLzDWdoMLr5cjSmq709fLl1Mr64DWyqi4CMN40xZrY/Tge4z8NMAIwKwe/mLr8Rw57hW5OWswK8xLxwF2R7dd0S6lVM3C1IRSqg7QC9ju6etlfe8WYeQOFxR7zMzrZffnpZT23gWssF6fxcC9yqhqaQ60BtY60ZZytUsp1RWYBQzSWh+z2W73PXVjuxrafDsISLV+vQy4zdq+msBtXP7J1KXtsratHcbA4Sqbba68Xo5YDDxkrV7pCZyxdlacv16uGsF15h8wBCNPdBE4Ciyzbm8EfGWzX39gF8Zf1BdstrfA+EXbA/wPqGRSu2oDy4HdwHdALev2WGC2zX7RGH9lg4o9fwWwFSMgfQBUdVe7gGut595s/f9Rb7hewANALrDJ5l8XV1wvez8vGKmaQdavw62vf4/1erSwee4L1uftBG43+ee9rHZ9Z/09KLw+i8t6T93UrsnANuv5vwfa2Tz3Eet13AM87M52Wb9/CZhS7Hmuvl7zMaqucjHi16PA48Dj1scV8Lq13Vuxqchz9nrJFH0hhPBxvpxaEUIIgQRyIYTweRLIhRDCx0kgF0IIHyeBXAghfJwEciGE8HESyIUQwsf9f/BkZ0pxs5YsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x[:, 0], y)\n",
    "\n",
    "xs = np.linspace(-1, 1)\n",
    "xs = np.stack([xs, np.ones_like(xs)], axis=1)\n",
    "ys_true = xs@beta_true\n",
    "ys_fit = xs@beta_ols\n",
    "plt.plot(xs[:, 0], ys_true, color='k', lw=2, label=\"True\")\n",
    "plt.plot(xs[:, 0], ys_fit, color='orange', lw=1, label='OLS fit')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380ba33a",
   "metadata": {},
   "source": [
    "# Gradient Descent\n",
    "\n",
    "Minimize the error metric, in this case MSE:\n",
    "$$\n",
    "    MSE = \\frac{1}{n}\\Sigma^n_{i=1}(\\hat{y_i}-y_i)^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cbd1bd0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse(y_hat, y): return jnp.power(y_hat - y, 2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "97540a87",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.145107\n"
     ]
    }
   ],
   "source": [
    "test = jnp.array([-1., 1.])\n",
    "print(mse(x@test, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "699cb172",
   "metadata": {},
   "source": [
    "Now we can use the power of JAX for automatic differentiation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5c497e02",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "48ec93e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(beta, x): return x @ beta\n",
    "def loss_fn(beta, x, y=y): return mse(predict(beta, x), y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0ab92e2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_loss_fn = grad(loss_fn, argnums=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4dad4a6b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.  1.] [-3.3981533 -2.7396417]\n"
     ]
    }
   ],
   "source": [
    "print(test, grad_loss_fn(test, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ca4a041",
   "metadata": {},
   "source": [
    "We will find the optimal value of $\\beta$ by iteratively updating it with the old value subtracted by a learning rate times the derivative of the loss function w.r.t. $\\beta$, i.e.:\n",
    "\n",
    "$$\n",
    "    \\beta_{n+1} = \\beta_{n} - \\gamma \\frac{\\partial\\mathcal{L}(\\beta)}{\\partial\\beta}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "af4d03ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_descent(beta, lr=1e-1, x=x, y=y): return beta - lr * grad_loss_fn(beta, x, y=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "70d53e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_gd = jnp.array([-1., 1.])\n",
    "lr = 5e-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b2ec519b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 [0.17539409 1.7678423 ]\n",
      "10 [1.766862  2.0548797]\n",
      "20 [2.4689775 1.9785113]\n",
      "30 [2.7926652 1.9215833]\n",
      "40 [2.9433815 1.8928522]\n",
      "50 [3.0137112 1.8792199]\n",
      "60 [3.0465453 1.8728325]\n",
      "converged -> relative error: 0.0009619885240681469)\n"
     ]
    }
   ],
   "source": [
    "for epoch in range(200):\n",
    "    if (epoch % 10) == 0: print(epoch, beta_gd)\n",
    "    \n",
    "    beta_new = gradient_descent(beta_gd)\n",
    "    \n",
    "    rel_err = np.sqrt(np.mean((beta_gd - beta_new)**2))\n",
    "    \n",
    "    if rel_err < 1e-3:\n",
    "        print(f\"converged -> relative error: {rel_err})\")\n",
    "        break\n",
    "        \n",
    "    beta_gd = beta_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9cb25f99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([3.0570934, 1.8707792], dtype=float32)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_gd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "079e59f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7pUlEQVR4nO3deVxUZdvA8d8NgqCmuO+KmvuGimWSmlppWoiaT2WLVk+WmeUCZfW8vdVTbxYjpmmpTy4tVpYp6pNli5Vm5k6m4m6W5IIiKorKcr9/nGFCHJZhzqxc38/Hj3DmzDn3nBmuuc91b0prjRBCCN8V4OkCCCGEcI4EciGE8HESyIUQwsdJIBdCCB8ngVwIIXxcOU+ctEaNGjo8PNwTpxZCCJ+1ZcuWk1rrmgW3eySQh4eHs3nzZk+cWgghfJZS6rC97ZJaEUIIHyeBXAghfJwEciGE8HEeyZHbk5WVxZEjR7h48aKniyK8UEhICA0aNCAoKMjTRRHC63hNID9y5AjXXHMN4eHhKKU8XRzhRbTWnDp1iiNHjtCkSRNPF0cIr+M1gfzixYsSxIVdSimqV69Oamqqp4sifEzithTiV+3hr/RM6oWFEtevJTGd6nu6WKbzmkAOSBAXhZLPhnBU4rYUnl3yG5lZOQCkpGfy7JLfAOwGc18O+l4VyIUQwizxq/bYgniezKwc4lftuSpA2wv6cZ/9yksrdpJ+IcvrA7sEcqtTp07Rt29fAI4dO0ZgYCA1axoDqDZu3EhwcLAniyeEcNBf6Zkl3m4v6Gflak5fyAKKr817mimBXCkVBrwLtAM08JDWer0Zx3aX6tWrk5SUBMCLL75IpUqViI2NtT2enZ1NuXLyvSeEr6gXFkqKnaBdLyz0qm2FBf38CqvNewOzItM04Cut9Z1KqWCggknH9aiRI0cSEhLCtm3biIqKonLlylcE+Hbt2vHf//6X8PBwPvzwQ6ZPn87ly5e5/vrrefvttwkMDPTwKxCi7Irr1/KKdAlAaFAgcf1aXrVvYUG/oJIEfE9wekCQUqoK0BOYC6C1vqy1TnfymC75VxpHjhzh559/JiEhodB9kpOTWbRoEevWrSMpKYnAwEAWLlxY2pcvhDBBTKf6vDakPfXDQlFA/bBQXhvS3m6NOq5fS0KDiq942avNewMzauRNgFRgvlKqI7AFeEprfT7/TkqpUcAogEaNGplwWvcYNmxYsTXr7777ji1bttC1a1cAMjMzqVWrljuKJ4QoQkyn+iVKheTtk9drpUpoEOcvZ5OV8/eaxoXV5r2BGYG8HNAZGKu13qCUmgZMAv4n/05a6znAHIDIyMgiV3z2pgWhK1asaPu5XLly5Obm2n7PG4WqtWbEiBG89tprbi+fEMIcBYO+L3VHNCOQHwGOaK03WH9fjBHI/U54eDj//e9/Adi6dSuHDh0CoG/fvgwaNIjx48dTq1Yt0tLSOHfuHI0bN/ZkcYUQTihpbd4bOJ0j11ofA/5USuXdc/QFdjl7XG80dOhQ0tLSaNu2LTNmzKBFixYAtGnThldeeYVbb72VDh06cMstt3D06FEPl1YIUVYoM9IYSqkIjO6HwcBB4EGt9enC9o+MjNQFF5ZITk6mdevWTpdF+C/5jIiyTim1RWsdWXC7Kd0PtdZJwFUHF0II4XoyH7kQQvg4GaoohBAu5uoeMBLIhRDChRydhbE0JLUihBAuVNQsjGaRGrkQwm944yAeR2ZhLC2pkedz/Phxhg8fTtOmTenSpQs33HADS5cudeqYL774IhaLBYAXXniBb7/9tlTHSUpKYuXKlXYf++GHH6hSpQqdOnWiZcuW9OzZ0zZwyVN+//13Pvroo0Ifj4uLo23btsTFxTFr1izef/99ABYsWMBff/3lrmIKP5KXwkhJz0TzdwojcVvKVftFTV5N+KQvaPbsSsInfUHU5NVX7WeWwuZnMXPeFqmRW2mtiYmJYcSIEbYAdPjwYZYvX37VvqWd0vbll18udfmSkpLYvHkzAwYMsPt4jx49bME7KSmJmJgYQkNDbXOsu1teIB8+fLjdx+fMmUNaWtpV89gsWLCAdu3aUa9ePXcUU/iRkiwkUTBfnWMdR+PK+cYdmYWxtKRGbrV69WqCg4N57LHHbNsaN27M2LFjASPAREdH06dPH/r27UtGRgZ9+/alc+fOtG/fnmXLltme9+qrr9KiRQtuvPFG9uz5Ow82cuRIFi9eDMCWLVvo1asXXbp0oV+/fraRoDfddBPPPPMM1113HS1atGDt2rVcvnyZF154gUWLFhEREcGiRYuKfC0RERG88MILzJgxA4DU1FSGDh1K165d6dq1K+vWrQPgxx9/JCIigoiICDp16sS5c+cAeP3112nfvj0dO3Zk0iRjtoUDBw7Qv39/unTpQo8ePdi9e7ftNT355JN0796dpk2b2l7fpEmTWLt2LREREUydOvWK8kVHR5ORkUGXLl1YtGiR7a5l8eLFbN68mXvvvZeIiAgyM71zylDhnUqSwrAX7POYnbfO48gsjKUlNXKrnTt30rlz5yL32bp1K9u3b6datWpkZ2ezdOlSKleuzMmTJ+nWrRvR0dFs3bqVTz75hKSkJLKzs+ncuTNdunS54jhZWVmMHTuWZcuWUbNmTRYtWsTzzz/PvHnzAKPGv3HjRlauXMlLL73Et99+y8svv8zmzZttwbk4nTt3Jj4+HoCnnnqK8ePHc+ONN/LHH3/Qr18/kpOTsVgszJw5k6ioKDIyMggJCeHLL79k2bJlbNiwgQoVKpCWlgbAqFGjmDVrFs2bN2fDhg08/vjjrF69GoCjR4/y008/sXv3bqKjo7nzzjuZPHkyFovFbopn+fLlVKpU6YqFPADuvPNOZsyYgcViITJSxpcJx5RkIYni8tKumm/c1fO2eG8g/8gFi+0OL/l0BGPGjOGnn34iODiYTZs2AXDLLbdQrVo1wEjFPPfcc6xZs4aAgABSUlI4fvw4a9euZfDgwVSoYKytER0dfdWx9+zZw44dO7jlllsAyMnJoW7durbHhwwZAkCXLl34/fffS/VS80+98O2337Jr19/T35w9e5aMjAyioqKYMGEC9957L0OGDKFBgwZ8++23PPjgg7byV6tWjYyMDH7++WeGDRtmO8alS5dsP8fExBAQEECbNm04fvx4qcorhLNKksIobgEJb51vvDjeG8gdCLpmaNu2LZ9//rnt95kzZ3Ly5Mkraob5p7RduHAhqampbNmyhaCgIMLDw23T2hZHa03btm1Zv97+anjly5cHIDAwkOzs7NK8HLZt22ablyQ3N5dffvmFkJCQK/aZNGkSAwcOZOXKlURFRbFq1Sq7x8rNzSUsLMxWgy6svOBdUxCLsqXgnOL2eq3YC/Z5XDnfuKt700iO3KpPnz5cvHiRd955x7btwoULhe5/5swZatWqRVBQEN9//z2HDx8GoGfPniQmJpKZmcm5c+dYsWLFVc9t2bIlqamptkCelZXFzp07iyzfNddcY8thF2f79u38+9//ZsyYMQDceuutvPXWW7bH8wLygQMHaN++Pc888wxdu3Zl9+7d3HLLLcyfP9/22tPS0qhcuTJNmjThs88+A4xg/euvv5pWXjOeJwQYwXzdpD4cmjyQdZP6XBUs8+erAQKtK4fZy1vn9W5pkq9XS/5tES99TaeXv77icXtK2pvGGd5bI3czpRSJiYmMHz+eN954g5o1a1KxYkVef/11u/vfe++93HHHHbRv357IyEhatWoFGLnpu+66i44dO1KrVi3bqkH5BQcHs3jxYp588knOnDlDdnY248aNo23btoWWr3fv3kyePJmIiAieffZZ7rrrriseX7t2LZ06deLChQvUqlWL6dOn23qsTJ8+nTFjxtChQweys7Pp2bMns2bN4s033+T7778nICCAtm3bctttt1G+fHmSkpKIjIwkODiYAQMG8H//938sXLiQ0aNH88orr5CVlcXdd99Nx44dCy1vhw4dCAwMpGPHjowcOZLx48cX+x6A0Xj62GOPERoayvr16wkN9c1bXeG9SpKvtjcaM27xr6AhK9e460zPzLLtX1ivl8RtKUz89Fdb75g8Zi/kbMo0to6SaWxFachnRBTFzPRF1OTVJVqMuaD6YaGsm9QHgH8l/sbCX/6gsAirgEOTBzp0fJdOYyuEEJ5k9nwmpe29kve8xG0pRQZxMLdhVXLkQgifZ/Z8JqUNsnnPi1+1p8gg7tcDgqTHgyiMfDZEUcyezySuX0tCgwKL3zGf/MHZ3nkvpx7m5Mo3yT5+0H8HBIWEhHDq1CmqV6+OUi7oQy58ltaaU6dOXdV9Uog8JRkM5Ah7XRnPX8q+ooEzv/oFcvJ55dFac+mP3zi7cQmZB412wS4NKhHTaWypylUYrwnkDRo04MiRI6Smpnq6KMILhYSE0KBBA08XQ3gpV8xnUrB3S8E8fN457NWuJ/RtxpOvzeLk+sVcPrbf2Dc0mCcfiOSR2PhSl6kwXhPIg4KCaNKkiaeLIYTwQSUZDOSOc2RkZDB37lymTp3KX9axJdVqVObpB5oxtuufVKhTDRrXNq1Mebym+6EQQriTmd0Vjx49yltvvcU777xDeno6AD26NGHGE81oX3ELqkE0tI6FsHZOlVm6HwohhJVZ3RV37drFlClT+PDDD7l8+TIAwwdG8H/3V6NR4DZUk/bQah5UbGj+i8hHArkQwi8VVeMuydzlhdFa8+OPP2KxWPjiiy8AY2T4s4/0InZADtVy90GLodD8MyhfzTUvrgAJ5EIIv5EXvFPSM1Fg68tdsMZdmu6K2dnZLFmyhPj4ePJSw6Gh5Zkyvjcjrz9BKH8a6ZMmI6Gce6eWkEAuhPALBdMlBVv/8te4HemumJGRwfz585k6dSqHDh0CoG7t6rz9dBS3N9tNueAT0OYZaDgUAhzre24WCeRCCL9Q1Oo/efJq3CXprnjs2DFmzJjB22+/zenTpwGIaNuUmRM60S1sPQFhmdBmJtTuC3bGvrhzIWgJ5EIIn5e4LaVEk1zl1biL6kqYnJxMQkIC77//vq0Bc2DfLkwZ1ZgWAT+i6gZB6xVQrfAVxcye+6U4EsiFED4tL2gWp2CNO/+AH601a9euJTp6tG0NAaUUj99/C/9zZwXqXPoRwq+H1huhUtNiz+VMY2ppSCAXQvi0olIqeQ2eBYfQ58lbezc+Pt62pGNISAj/+8QAxvQ5zzUXtkD4aGgxB0JqlbhMZs/9UhwJ5EIIn1ZUcJx6V4TdGvD58+eZP38+CQkJtgbM6tWrkRDXn7s7HCb40iZoMgGaLYagSg6Xyey5X4ojgVwI4TVK00BYWNCsHxZ61XOPHz9ua8BMS0sDoEXzpsx8uhd96mwmQG2HVk9D47shIKjUr8MVc78URQK5EMIrlLaBsCRBc/fu3bYGzEuXLgFw041dmTq2Ax2Dv0VVOgCtX4N6A+z2QHGUO+Z+yU8CuRDCK5S2gbCwoDkooh5r164lPj7+igbM+4b159UR9Wh0cQXUbABtPoEa3Ux/PSVZG9QsEsiFEF7BmQbC/EEzJyeHpUuX0u2xeDZu3AhA+fLlmfjoUGJvD6DqmS+gxhBovRYquybV4W4SyIUoA9w5OKW0nG0gzGvAnDp1KgcPHgSgevXqvDxhKA91O0lI2ldQ458QtQMq1DO17J4mgVwIP+fuwSmlVdoGQnsNmM2aNSU+7naimyUTeG4F1BkHN86D4CqufAkeI4FcCD/n7sEppeVoA+GePXuYMmXKFQ2YN3S7jinje9Kt8g+orK8gPA6arIDA8m57HZ5gWiBXSgUCm4EUrfXtZh1XCOEcdw9OcUZxDYRaa9atW0d8fDzLly+3bR8aM5DXHm3DtVmJqOC10OY5aDAIlFetL+8yZtbInwKSgcomHlMI4SR3D05xhZycHBITE4mPj2fDhg2A0YD56EN38fw/alAr7SMI0dD5XajZo9RdCH2hLcEeUwK5UqoBMBB4FZhgxjGFd/LVD3pZ5u7BKWa6cOECCxYsICEhgQMHDgBQrVo1nnnyPsb0zaLisUUQPAD6fA1h7Z06l6+0JdhjVo38TeBp4BqTjie8kC9/0MuK4lbF8ZUv4BMnTjD2X6+xZOE8si+cBaB2/UZMee5+bmiwmxpn5rNoS2/+e3kGDzTrSUyY86/FV9oS7HE6kCulbgdOaK23KKVuKmK/UcAogEaNGjl7WuEBvvxBLwuK+6L1hfdo7969TJkyhfkL3iPrstGAGVy3BX1vu46JXfdxfeV3mHdoIPNOzCY9x8jiJptUmfCltoSCzKiRRwHRSqkBQAhQWSn1odb6vvw7aa3nAHMAIiMjCy7eIXyAL3/QywJf/qLN34CptREeKjS/jiEDWzO29SbqBH3Df04O5sk/nuZ87pU9UMx6ja5oS3BXKtLpQK61fhZ4FsBaI48tGMSFf/CHRjN/5mtftDk5OSxbtgyLxcL69esBowFz5APDCamdw+gmP5CtTzLrxFBWnrmRHApfRi0lPZOoyaudCphmtyW4MxUp/chFiflyo1lZ4K1ftAVrpU/2bETar9+QkJDA/v37AahatSrjx/6TcbdX5JqUd9mcVpNX/vonazM6YcwqbghUihx99Q29AttrtxcwS1IzNrstwZ13SKYGcq31D8APZh5TeA9fbDQrS7ztizZxWwovrdjJ6QtZAORcOMPOnxYy/OUvyM00GjBr1WvIq88+xgPd0gk+PA8u9eGHWrMZv72c7Xl5QoMCGdqlPp9vSbniNeYtHpFf/oDpSM3YzLYEd94hSY1cOMRXGs3KophO9dl8OI2PN/xJjtYEKsXQLp55v/IHz6y0FM5uSuT8ju/Q2cYamMF1mtOhd28mdP+ToTXeIDj3Hrj1FxL3h1qfd2UQDwsN4sXotsR0qk9k42pXVCYKW6szL2B6qu3AnXdIEsiF8BOJ21L4fEuKLfWQozWfb0khsnE1twfz+FV7SP99B2c2LiFz7y/k1ZlDm3Ul6pZuPNl+G1GVPuGjU/0Z8se7/PfuO63PW2132baK5cvZXkPBykTU5NVFBkxPtR248w5JArkQfsIbeq3k5OSwfPlytsx4nkspycbGwHJUatubgbe2Zsy1a2hW/mPmnhzEpCNPcj63AvnHYJYm6BYXMD3VduDOVKQEciH8hCd7rWRmZvLee++RkJDAvn37AAgIqUSVzv255+ZajG60imC1h9mpQ1me3pMs/fcyavkDammCbnEB05NtB+5KRUogF8JPeKLmefLkSWbOnMmMGTM4efIkAOHh4fT/x33UqX6ckTVXcDy7GgnH7uX7c5ForpzEqmBALW3QLSpgloVGegnkQvgJd9Y89+/fT0JCAgsWLCAz0/jy6NKlC8/Hjia69RECD7zN0XKdeHH/c3x9oin1wkKZepdRjqICqquCrr830ittp0+mq0VGRurNmze7/bxC+DtXjyRcv349FouFpUuX2kZgDhw4kOfH3U+3sJ9Rv38ADWKgdRxUaW3aeYVBKbVFax1ZcLvUyIXwI66oeebm5rJ8+XIsFgvr1q0DIDg4mPvuu4/nHo+m2aXP4a/HofpDMOA3qFCy88tMmuaRQC6EsCszM5P333+fKVOm2Boww8LCePzx0Uy4ryvVj78Lhx+Dlk9B5HQIDivxsWUmTXNJIBdCXOHkyZO8/fbbzJgxg9TUVAAaN27M+PFPMWpAbUIPTof9i430SY/PITDE4XN4Q1dJfyKBXAgBwIEDB5g6dSrz5s2zNWB27tyZZ2LHMbTLBQL3JMCBKtDmGSMPHlD4JFbF8bUJvrydBHIhyrgNGzYQHx/PkiVLbA2Yt912G5MmjqZHnZ2ovc9ASke4bhbUuqnUy6jl560TfPmqsrEyqRDiCnkNmD169KBbt258/vnnlCtXjgcffJDkrd+x8tX29Dw9EpX+G9z0JfT+Emr3NiWIg9FVMjToyhq9zKRZelIjF6IMuXjxoq0Bc+/evQBUqVKF0aNHM+7h26h96n3YPxSa3A/9t0ClcJeUoywM0nEnCeRClAGnTp2yNWCeOHECMJZcHD9+PI8MaU/F32fCr3Oh+eNwxz4IqeHyMvn7IB13kkAuvJo/9zV2x2uz14DZqVMn4mJj+UePawjcY4Ftb0KridD9AyhX0dTzF8af31dPkEAuvJY/9zV29WvbsGEDFouFJUuWkJubCxgNmHETx3FTk+Oo3ZPhN2X0QGk0DAKCijmiefz5ffUUaewUXquovsa+zhWvLTc3lxUrVtCzZ0+6devG4sWLCQwMZMSIEexM2sDKhH70PvcI6tB8iHgDbkuC8OFuDeLg3++rp0iNXHgtf+5rbOZru3jxIh9++CFTpkxh9+7dgNGA+dhjjzHusXupc/Yz2DcQavWCGz+DGtc5VXZnmf2+SppGauTCixXWp9gf+hqb8dpOnTrFq6++Snh4OI888gi7d++mYcOGJCQkcGT3GiYPzaDO5l5w8Tjc8jP0WOzxIA7mvq95aZqU9Ew0Rppm3KIkIl76msRtKU6W1HdIIBdey5/7Gjvz2g4ePMjYsWNp1KgR//rXvzh+/DgRERF8+OGHHNi8mPHXbaTSuj4QVBkG7oLrZkPl5q56KQ7r3aqmQ9uLYi9NA5CemcWzS34rM8FcUivCa/lzX+PSvLZNmzYRHx/P559/bmvA7N+/P7ETJ9KnLajkN2DdLmg13gjeQZWLLIOnUhLf7051aHtRikrHlKW5WySQC6/mz32NS/LacnNzWblyJfHx8axZswaAoKAg7r//fiZOGEf7KnsheRJsOQ+tn4bweyEwuNhze7LniJk58sKG+jtzTF8kqRUhvNDFixeZO3cu7dq144477mDNmjVUrlyZp59+mkP7d7Hg+W60//1O2DMN2r0AA3dCswdLFMTBsz1HzMyR20tROXtMXyQ1ciG8SFpaGrNmzWL69OkcP34cgIYNGzJu3Dj++cBQKh/7EDbfCNWvgxveg5pRpTqPJ3sEmbkkXd7dw0srdnL6QtYVj/lLe0pJSCAXwgscOnSIqVOnMnfuXC5cuABAx44diYuL4x+330DQ/hnwQyeoHw19voOwtlcdw5GctydnHzS77SMvRVWWuyHKmp1CeNCmTZuwWCwsXrzY1oDZr18/Jk6cyM1d66KSLZCyHJo+aDRiVmhg9zgFc95g1EhfG9LebjBzdH/hHWTNTiG8RG5uLl9++SXx8fH8+OOPAJQrV4777ruPiRMm0KHuOdj1OqzeBC3GQvQBCK5a5DEdXXHHn3sElUUSyIUoJUdv5S9dusTChQuxWCwkJycDULlyZR599FGeHPsEDdQ22DUaDp+A1hPhxk+hXMlSHaXJeZe0R1BZTln4CgnkQpSCI933Tp8+bWvAPHbsGAANGjRg3LhxPPLwCCqfWg6/9jNmHmzzDDQY4vAyaoXlvMMqODePikxw5Ruk+6EQpVCS7nu///4748aNo2HDhjz33HMcO3aMDh068MEHH3BwTxITB2gqfx8BhxdB5Azot8k6E6Hja2HG9WtJUODVq/dkXMx2anSjTHDlGySQC1EKRaUytmzZwj333MO1117LtGnTOH/+PDfffDOrVq0iaf1X3NduJ0FftoS0LdBrBfRZBXX6OrWMWkyn+lQMvvoGOytXOxV0/XniMn8iqRUhSqFgKkNrzcWDm7m4dRmRrycB+RowJ04komlF2G2BlXdD4+FG7btSE1PLdCYzy+52Z4KuLJLsG6RGLkQp5I0o1NlZZGz/hqNzx3Bi8UucPZjENddcw8SJEzl48CAfvPkkEedegW+6Q0gduH0PdJ1hehAH18wW6c8Tl/kTqZGLMseMXhi9wisQeeZHFi34D1nnTgFQvVYdnomdwKhHHqHKhV9g1wOQcRBaTYBuCyCokgtezd/MHDGZR7op+gYJ5KJMcbYXxuHDhxn73Ct8sfhDci9fBKDxta349wvPcdewoQQfXQo/9wK0MYlV47vctgKPq4KuP09c5i8kkJcBnu4H7Onz5+fowJk8W7duJT4+nk8/+4zcHOP5IY0jqHz9EKq2aEfHejsIXtUGKoZDxGSo27/QxktXXg8JumWTBHI/5+l+wO48f0kCpCO9MLTWfPXVV1gsFlavXg2ACgikYtveVO46mNr1qjOixn+5r/oUdhzqADGLoMb1xZZR+mULszkdyJVSDYH3gdqABuZorac5e1xhjtLWQH3t/CUNkCXphXH58mU++ugjpkyZwo4dOwCoVKkSo0aN4tPM9jSunss/ay4lJuwHvjzTnX8ceJ1DlxpwqJggDp5/P0rDm+6ohH1m1MizgYla661KqWuALUqpb7TWu0w4tnCSp/sBu+v8JQ2QRTUIpqenM3v2bKZPn85ff/0FQL169XjqqacYNWoUYfow/RZPpGv5jXyS1o9b9r5NanY1AOqXsGeIp98PR8kdhG9wOpBrrY8CR60/n1NKJQP1AQnkXsDT/YDddf6SBkh7DYIj2ldkzQdTuP8//yEjIwOAdu3aERsbyz13303w6XWw7W5I3069Jg9y69rRpF4KsR3TkZ4hYRWCrpo3O2+7N/LFO4iyyNR+5EqpcKATsMHOY6OUUpuVUptTUx1fm0+Ujqf7Abvr/I70oY7pVJ91k/qw5K56hP82l8cHRTF16lQyMjLo27cvX375JduTtjGid0WCv78RNj9uDJ2PPkS7m1/l+Zhu1A8LRWHUxB2Z+rWwWaNdOZt04rYUoiavpsmkL4iavNqhIfu+dgdRVpnW2KmUqgR8DozTWp8t+LjWeg4wB4z5yM06ryiap/sBu+v8Je1DrbVm1apVWCwWvvvuOwACAwMZPnw4sbGxdOrQGg6+B1+0hvLVoe1z0GAQqL/rPM70DCls9GVh253lbGrE03d0omRMCeRKqSCMIL5Qa73EjGMK83i6S5o7zl/cF8bly5f5+OOPsVgsVzRgPvLII4wbN45GdSrDvndg2QCo1hm6zYWaPZya/8QedwdGZ1MjZg0ykgZT1zKj14oC5gLJWusE54skROnY+8JIT09nzpw5TJs2zdaAWbduXZ566ikeffRRwoLPw543YeM8qDcQ+nwNYe1dVkZXjL4sirOpETPuqKTB1PXMqJFHAfcDvymlkqzbntNarzTh2EKUyh9//MG0adP4z3/+w7lz5wBo27YtsbGxDB8+nODMA5A8AY4kQpMH4LZtULGRy8vl7lSXGXcAzt5RSYOp65nRa+UnwNz7T1FmOXsLnpSUhMVi4ZNPPiHHOgKzT58+xMXF0a9fP9TJ9bB+GJzaAC2egDv2Q/lqrno5drkz1eXuOwB7pMHU9WRkZyEcCSiS/zOHvVvwuMW/8uLynZzJzCr02mqt+frrr4mPj7+iAfOee+4hNjaWzp0i4K+V8G1PuJACrWMh6pMSL6Pmyzzd2A3SYOoOEsjtcCSnJ/k/89i7Bc/K0aRbe3QUvLaXL1/mk08+wWKx8NtvxvaKFSvaGjAbN6gLhz+GlfdDQHljGbWGQyHAPz72Ja1AeLqx2xvuCvydf3yiTeZITk/yf+Ypya12ZlYOry3bwv7vPubNN98kJcXoE123bl2efPJJHn30UapWKgf7/wMrpkLlVtBlGtR2bgUeb+NLFQhvuCvwdxLI7XAkpyf5P/MUdgueJ/tsKue2rOCPpC/ZeNnYr02bNrYGzPK56bDXAvtnQ+2boecyoyuhH/K1CoSn7wr8nQRyOxzJ6Un+zzz2bsEBLp84yNmNSzmfvAZyjcd69+5NbGws/fv3J+D8Qfj1KWMR4/B74NYNcE0zT7wEt0jcllLoF55UIMomWerNDkeGlXt6CLw/ielUn9eGtDcmoNKacke3c2LR/3B0/pOc3/k9aM01bXoR/+EXrF69mgHX1yJg3V3w9Q1QvibcsQe6vu33QTwvhWKPVCDKJqmR2+FITk/yf+Ya2K4W53asxvKxhe3btwMQEBxCxfa30rzvXfzPXT2JqbsLvusL5/ZZl1Gb7/Jl1LyFvZRKfr1b1XRjaYS3kEBeCEdyepL/c97Zs2dtIzCPHDkCQJ06dRg7diyjR4+mapVr4I/PIPkOOJplLKMWfo/bllHzFsWlTr7fLRPSlUUSyIVH+8EfOXKEadOmMWfOHM6eNeZaa926NbGxsdx7772UD8yBA/NgzRSo0BA6vAL1brtiEquypLgGYcmRl00SyMs4T3Vj2759OxaLhY8//pjs7GwAevXqRVxcHLfddhsBWadh7+uwbybUuAG6fwQ1b3BZeXxFYQ3CeSRHXjZJIPcQbxkN+tKKnW7rxqa15rvvviM+Pp6vv/4agICAAP7xj38QGxtL165d4fxh2Doefv8AGgyGvj9AldalOp+3XGMz5ZX/xeU7bQOl8kgje9klgdwDvGUwR+K2FLur1YC5t+hZWVl8+umnWCwWkpKSAKhQoQIPP/ww48ePp0mTJpD+G/x8vzGUvtlDMOA3qFD6a+Et19gV8tpk/PGLSpSOBHIP8JbBHPGr9hT6WHG36CUJImfPnuXdd9/lzTff5M8//wSgdu3atgbMalWrwok18MMTkLYVWj4FkW9BcFiJX0Nh5fCWa+xK0sgu8kgg9wAzR4M6Uysr6nxF3aIXV9tNSUlh2rRpzJ4929aA2apVK1sDZkj5YGP62E1vwOU0aB0HPT6HwJDCTulwOWTErShLJJB7gFmjQV21jFdYaFCRzy+stvvie1+y9M21fPTRR1c0YMbGxjJgwAACdBYc+gCS441ad5tnoP4gCAi0c5biFVXrlhG3oiwpm324PMys0aBFBTJnyvFidNsin5e/Vqu1JvP3JI5/+gK/TnuE999/n9zcXIYNG8aGDRv44YcfuP3WHgTsjoflTeDPJXDdHLj1F2g4pNRBvGA5Cm6XEbeiLJEauQeYNRrUU8t41QsL5cipc1zY8xNnNy7l8vEDAAQEhTDmMWMK2aZNm8KFv2Db03BgrtH3+6avoGoHwJweJUXVus0ecSsNi8KbSSD3EDMaqjyxjNe5c+dombqGje++Q/bZEwAEVAij2nXRTHkhjgd6t4Oze2DDP43ad/j9cNtWqNjYdgyzepQUN8+1WY2B/twDRvgHCeQ+zJ0T9qekpDB9+nRmz57NmTNnjHPVbEiFLjE07z6AZ25vT0zDP2HNYEhdZ11GbR+Ur37VsczqUeKueW7KQg8Y4dskkPswdwSyHTt2YLFY+Oijj8jKMvqc9+zZk9jYWAYOHEiAUkbf7133wB9/GsuodV8I5SoUekwze5S4owue9IAR3k4CubiK1prvv/+e+Ph4vvrqK8AYgTls2DAmTpzI9ddfD7lZcHgh7HrDWDqt9dPQaFiJllHztR4l3l5eyd8LCeQ+zOzcbXZ2Np999hkWi4WtW7cCRgNmxfZ9ad7nboYP78317arA7jdhdwJc0xw6T4E6tzi0jJqvreHozeWV/L0ACeQ+zazc7blz55g7dy5Tp07ljz/+AKBKtRoEtR9ASMf+BIZWhsB0jq6N49LeVZSv3xd6LIHqkYDjNUKzUkLuqol685zzkr8XIIHcpzmbu/3rr7946623mDVrFunp6QC0aNGC2NhY5h1ryNHzOTQMPsaoGh8SHfYj/z3Tg/v/fJNP77rfdozS1gidzW27uybqrcPhJX8vQAYE+bTCcrTF5W537tzJQw89RHh4OJMnTyY9PZ0bb7yRxMREkpOTeeSRR6ieu4e3Gr3OsmsncDa3IjfvncXzKU+w6WS1K47l7KCk0vLUeb1NaT8Dwr9IIPdhjoxezGvAHDBgAO3atWP+/PlkZ2czdOhQ1q9fz9q1axkUHU3AidWw+hbmNXmFXy80p8fud4k/NoLU7KrA1QHCUzVCqYkaZASrAEmt+LSS5G6zs7NZvHgxFouFLVu2ABAaGspDDz3E+PHjadasGeRmw+FPYdfrkJMJbZ5mfdgcFi7dQ2Zu0Q18ruzRUVQO3Nt7kriLN+fvhftIIPdxheVuMzIybA2Yhw8fBqBmzZo88cQTPP7449SoUQOyM2HfO5BsgdC60P5/of7toAIYBGgVXGyAcFWPjuJy4N7ck8TdvDV/L9xHArmXc7RnxtGjR5kxYwbvvPMOp0+fBqB58+ZMnDiRBx54gNDQULiUBjtegb0zoPr1cMP7UDPqqmOVJEC4qkZYXG8MqYkK8TcJ5F6sJD0z8gL97/v3kP3rctJ+XU121mUAunfvTlxcHNHR0QQEBMD5P2HLVDi0ABoMgr6roUobp8vpihphSXLgUhMVwlDmArmjNdzC9ndHH+biaqVLtx5h3JsfkfrzYjIPbLLuoejW5zam/PtfdO/e3diUvsOYAzxlBTR9CAZshwoNTC2r2SQHLkTJlalA7mjf48L233w4jc+3pLi8D3NhtdKUtAwWLVrEwxP/l/MpRnc7VS6Yiu1vpnLkIAKaXkv3G26AE2uNBsy0zdDySejyJgRXNa18riQ5cCFKrkwFckdHwRW2/8cb/iRH6xIfp7QK1kpzL2eSsf0bLmxdzt3xxwAICK3MNZ1v55rOAwmsUAVFLm1zV8M3/4KLqcYkVj0WO7yMmqdJDlyIkitTgdzRvseFbS8YxIvbv7TyaqUZp09ydusKMratJPdiBgDXXnstue1uJ6tpDwKCyhOssogJ+5pHay7hsqoIrf4NDZxbgcfTJAcuRMmUqQFBjo6CK2x7YCETRJmdv20Zcpb6vy0gZdZDnF3/KbkXM2jVMZIlS5awe/duprwQS61KmlE1P2dNq4cZGPYTLx9/gj0dv7XOROi7QVwIUXJlKpA7OgqusP3vub6hy0bTaa358ccfueOOO2jTpg3fLVsEudkMHjyYdevWkZy0icGDBxN4+QQx6i3WtXmEyMp/8NChF3nudDyD+99HTGfvbsgUQpirTKVWHM27FrV/ZONqpuZvs7OzWbJkCRaLhU2bjB4oISEhjBw5kgkTJtC8eXNjx7N7jQE8fy6G8HspP3ALt1Zqwq2lPrMQwtcpXUi+16GDKNUfmAYEAu9qrScXtX9kZKTevHmz0+f1B+fPn2fevHlMnTqVQ4cOAVC9enWeeOIJxowZQ82aNY0dT26A5DeMnijNH4cWYyCkpgdLLoRwN6XUFq11ZMHtTtfIlVKBwEzgFuAIsEkptVxrvcvZY/uzY8eOMWPGDN5++23bCMxrr72WCRMmMGLECCpUqABaw19fGl0IMw5B64nGKMxyFT1ceiGENzEjtXIdsF9rfRBAKfUJMAiQQG5HcnIyCQkJvP/++1y+bIzA7NatG3FxcQwaNIjAwEBjGbVDHxo1cIDWz0Djf0BAkAdL7t1kuTNRlpkRyOsDf+b7/QhwfcGdlFKjgFEAjRo1MuG0xfOWP26tNWvXrsVisbBixQoAlFLExMQQFxf39wjM7POwZ66xjFrFJhDxOtTt79AyamWRLHcmyjq3NXZqrecAc8DIkbv6fN7wx52Tk2NrwNy4cSNgNGCOGDGCCRMm0KJFC2PHiyeNCaz2vQ21ekDUIqhx1XehKIQsdybKOjMCeQrQMN/vDazbPMqTf9znz59nwYIFJCQkcPDgQcBowBwzZgxjxoyhVq1axo4Zh4za9+8LoeGdcMtPULmFS8vmj2SRCVHWmRHINwHNlVJNMAL43cBwE47rFE/8cR8/fpyZM2cyc+ZM0tLSAGjWrBkTJkxg5MiRRgMmwOkk2PUGHPsamj0CA3ca84GLUpEJtkRZ53Qg11pnK6WeAFZhdD+cp7Xe6XTJnOTOP+49e/aQkJDAe++9x6VLlwC4/vrriYuLIyYmxmjA1BqOrTZ6oJzZAS3HwXWzIKiy6eXJ4y1tBK4mE2yJss6UHLnWeiWw0oxjmcXVf9xaa9atW0d8fDzLly8HjAbMQYMGERsbS1RUFEopyM2BPz4zauDZGdA6DsKXQ2B5U8pRGG9oI3AXmWBLlHV+O7LTVX/cOTk5JCYmYrFY+OWXXwAoX768rQGzZUvrF0XORTj4njEKM6QmtPsX1L8DlHtmRShrDYD2JtgqK3ckQvhtIAdzZ8+7cOGCrQHzwIEDAFSrVs3WgFm7dm1jx8unjXUw97wF1bpAt3lQ80a3dyEs6w2AZemORAi/DuRmOHHihK0B89SpUwA0bdrU1oBZsaJ1lOWFI7B7KhycD/Wjoc83ENbOqBXO+97ttUJvbgD0htWVhPAnEsgLsXfvXlsD5sWLFwHo2rUrcXFxDBkyxGjABDizy1hG7cgyaDISbvuVxL0BxM/aQ0r6Fyggr9O8O2uF3toA6K6aclm/IxFlS5maxrYk1q1bx+DBg2nVqhWzZ8/m4sWLREdHs2bNGjZs2MCwYcOMIH7iJ/gxGr7rA5WawR37oUsCiXsDeHbJb7bacMGRT3m1QleL6VSf14a0p35YKAqoHxbKa0Pae7w2WlRN2UyOzjEvhC+TGjlGA+ayZcuwWCysX78eMBowH3jgASZMmECrVq2MHXUuHFlhdCG8eMxYRi1qEZT7OzjYC1QFuatW6I0r7LirpuytdyRCuEKZDuQXLlzgvffeIyEhgf379wNGA+bo0aMZO3bs3w2YOZeN0ZfJ8cbal22egYZDIeDqy1eSgFSWa4Xuyt1Ll0RRlpTJQJ6ammprwDx58iQATZo0YcKECTz44IN/N2BmnYP9c4xGzCqtocs0qHNzkT1QCgtUecp6rdCdNWVvvCMRwhXKVCDft28fCQkJLFiw4KoGzMGDB1OunPVyZB6HvdNh/2yofTP0Wg7VOpfoHPYCVV6DZ32pFUpNWQgXKBOB/Oeff8ZisZCYmEjeikh33HEHsbGx9OjRwxiBCXBuvzGA549PofE90G8jVGrq0LkkUBVPaspCmMtvA3lOTg7Lly8nPj7e1oAZHBxsa8Bs3br13zuf2mw0YJ74AZqPhtt3Q0itYs9RWH9oXwpUMvpRCN/nd4E8MzPT1oC5b98+AKpWrcrjjz/OE088QZ06dYwdtYajXxur8JzbB60mQLf5EFSpROfxh5GD/vAahBB+FMjtNWCGh4fbGjArVbIG6NxsYxKr5DeMn9s8DY3vdngZNX8YOegPr0EI4QeB3F4DZmRkpG0Epq0BM/sCHJgHu6dAxUbQ4VWod1up50ApbX9ob0plyOhHIfyDzwby9evXEx8ff0UD5u23305sbCw9e/b8uwHz0qm/l1GrEQVRH0ONbk6fvzT9ob0tleHN87EIIUrOp4bo5+bmkpiYSFRUFN27d2fp0qUEBQXx8MMPs3PnTlasWEGvXr2MIH7+MGx+ClY0hwt/ws1roOcSU4I4GN0MQ4MCr9hWXH9odw1PL6nSvAYhhPfxmRp5ZmYmnTt3Zvfu3QCEhYXZGjDr1v17mbTV678ha8frXFd+I19cGECNzt/Q/7ouppenNN0MvS2VIV0lhfAPPhPIQ0NDad26NZmZmUyYMIGHHnro7wZMreHEjxz/5SXapW9n7tloYk+N4lxuRUJTUrkYlOKS4ORoN0NvTGX4UldJIYR9PhPIAWbPnk3VqlX/bsDMzYEjiUYf8KwzzP/jdub99RSX9d89ULypF4ZM5CSEcAWfCuQ1a9Y0fsi5CIfeN0ZhBleDts9Cg0HMfvbLq6aNBe/phSGpDCGEK/hUIOdyOuybBXumGXOfXP8u1Oxh60LojamLgiSVIYQwm0/1WiFtK5zZCX2+hpu+gFo9r+gHLr0whBBlkW/VyOv0Mf4VQlIXQoiyyLcCeQlI6kIIUdb4VmpFCCHEVSSQCyGEj5NALoQQPk4CuRBC+DgJ5EII4eMkkAshhI+TQC6EED5OArkQQvg4nx8Q5E1LpwkhhCf4dCD3tqXThBDCE3w6teJtS6cJIYQn+HQg97al04QQwhN8OrXiC/OPg+TxhRCu5VSNXCkVr5TarZTarpRaqpQKM6lcJeIL84/n5fFT0jPR/J3HT9yW4umiCSH8hLOplW+AdlrrDsBe4Fnni1RyMZ3q89qQ9tQPC0UB9cNCeW1Ie6+q7UoeXwjhak6lVrTWX+f79RfgTueK4zhvn39c8vhCCFczs7HzIeDLwh5USo1SSm1WSm1OTU018bTerbB8vbfl8YUQvqvYQK6U+lYptcPOv0H59nkeyAYWFnYcrfUcrXWk1jqyZs2a5pTeSYnbUoiavJomk74gavJql+StfSGPL4TwbcWmVrTWNxf1uFJqJHA70FdrrU0ql8u5azCRrCMqhHA1p3LkSqn+wNNAL631BXOK5B5FNUKaHWS9PY8vhPBtzvYjnwGUB75RSgH8orV+zOlSOaGkfbalEVII4S+c7bVyrVkFMYMj6RJfGUwkhBDF8ekh+gU50mfbXiOkwgj+rmr4FEIIV/DpIfoFOZIuyd8ImZKeiQLyWmplFkUhhC/xqxq5o322YzrVZ92kPtQPC6VgdxsZfSmE8BV+FchL22dbGj6FEL7MrwJ5aedekdGXQghf5lc5cihdn+24fi2v6O0CMvpSCOE7/C6Ql4aMvhRC+DIJ5FYy+lII4av8KkcuhBBlkQRyIYTwcRLIhRDCx0kgF0IIHyeBXAghfJzyxFoQSqlU4HApn14DOGliccwi5XKMlMsxUi7HeGu5wLmyNdZaX7XEmkcCuTOUUpu11pGeLkdBUi7HSLkcI+VyjLeWC1xTNkmtCCGEj5NALoQQPs4XA/kcTxegEFIux0i5HCPlcoy3lgtcUDafy5ELIYS4ki/WyIUQQuQjgVwIIXycVwZypdQwpdROpVSuUqrQbjpKqf5KqT1Kqf1KqUn5tjdRSm2wbl+klAo2qVzVlFLfKKX2Wf+vamef3kqppHz/LiqlYqyPLVBKHcr3WIS7ymXdLyffuZfn2+7J6xWhlFpvfb+3K6XuyveYqdersM9LvsfLW1//fuv1CM/32LPW7XuUUv2cKUcpyjVBKbXLen2+U0o1zveY3ffUTeUaqZRKzXf+f+Z7bIT1fd+nlBrh5nJNzVemvUqp9HyPufJ6zVNKnVBK7SjkcaWUmm4t93alVOd8jzl3vbTWXvcPaA20BH4AIgvZJxA4ADQFgoFfgTbWxz4F7rb+PAsYbVK53gAmWX+eBLxezP7VgDSggvX3BcCdLrheJSoXkFHIdo9dL6AF0Nz6cz3gKBBm9vUq6vOSb5/HgVnWn+8GFll/bmPdvzzQxHqcQDeWq3e+z9DovHIV9Z66qVwjgRl2nlsNOGj9v6r156ruKleB/ccC81x9vazH7gl0BnYU8vgA4EtAAd2ADWZdL6+skWutk7XWxa18fB2wX2t9UGt9GfgEGKSUUkAfYLF1v/eAGJOKNsh6vJIe907gS631BZPOXxhHy2Xj6eultd6rtd5n/fkv4ARw1cg1E9j9vBRR3sVAX+v1GQR8orW+pLU+BOy3Hs8t5dJaf5/vM/QL0MCkcztVriL0A77RWqdprU8D3wD9PVSue4CPTTp3kbTWazAqboUZBLyvDb8AYUqpuphwvbwykJdQfeDPfL8fsW6rDqRrrbMLbDdDba31UevPx4Daxex/N1d/iF613lZNVUqVd3O5QpRSm5VSv+Sle/Ci66WUug6jlnUg32azrldhnxe7+1ivxxmM61OS57qyXPk9jFGry2PvPXVnuYZa35/FSqmGDj7XleXCmoJqAqzOt9lV16skCiu709fLYysEKaW+BerYeeh5rfUyd5cnT1Hlyv+L1lorpQrtu2n9pm0PrMq3+VmMgBaM0Zf0GeBlN5arsdY6RSnVFFitlPoNI1iVmsnX6wNghNY617q51NfLHyml7gMigV75Nl/1nmqtD9g/gulWAB9rrS8ppR7FuJvp46Zzl8TdwGKtdU6+bZ68Xi7jsUCutb7ZyUOkAA3z/d7Auu0Uxi1LOWutKm+70+VSSh1XStXVWh+1Bp4TRRzqH8BSrXVWvmPn1U4vKaXmA7HuLJfWOsX6/0Gl1A9AJ+BzPHy9lFKVgS8wvsR/yXfsUl8vOwr7vNjb54hSqhxQBePzVJLnurJcKKVuxvhy7KW1vpS3vZD31IzAVGy5tNan8v36LkabSN5zbyrw3B9MKFOJypXP3cCY/BtceL1KorCyO329fDm1sgloroweF8EYb9pybbQefI+RnwYYAZhVw19uPV5JjntVbs4azPLy0jGA3dZtV5RLKVU1LzWhlKoBRAG7PH29rO/dUozc4eICj5l5vex+Xooo753Aauv1WQ7crYxeLU2A5sBGJ8riULmUUp2A2UC01vpEvu1231M3lqtuvl+jgWTrz6uAW63lqwrcypV3pi4tl7VsrTAaDtfn2+bK61USy4EHrL1XugFnrJUV56+Xq1pwnfkHDMbIE10CjgOrrNvrASvz7TcA2Ivxjfp8vu1NMf7Q9gOfAeVNKld14DtgH/AtUM26PRJ4N99+4RjfsgEFnr8a+A0jIH0IVHJXuYDu1nP/av3/YW+4XsB9QBaQlO9fhCuul73PC0aqJtr6c4j19e+3Xo+m+Z77vPV5e4DbTP68F1eub61/B3nXZ3lx76mbyvUasNN6/u+BVvme+5D1Ou4HHnRnuay/vwhMLvA8V1+vjzF6XWVhxK+HgceAx6yPK2Cmtdy/ka9HnrPXS4boCyGEj/Pl1IoQQggkkAshhM+TQC6EED5OArkQQvg4CeRCCOHjJJALIYSPk0AuhBA+7v8BPOkKLUGxzdEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x[:, 0], y)\n",
    "\n",
    "xs = np.linspace(-1, 1)\n",
    "xs = np.stack([xs, np.ones_like(xs)], axis=1)\n",
    "ys_true = xs@beta_true\n",
    "ys_fit = xs@beta_gd\n",
    "plt.plot(xs[:, 0], ys_true, color='k', lw=2, label=\"True\")\n",
    "plt.plot(xs[:, 0], ys_fit, color='orange', lw=1, label='Gradient Descent fit')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2afcb6e",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}