ML.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximating electromagnetic shielding using machine learning\n",
    "\n",
    "The goal of this notebook is to apply standard machine learning techniques to our dataset. We know that the lab already obtained a RMSE of $0.6$ using kernel ridge regression. We seek to validate this result in first place and then we eventually expect to improve it. In order to do so, we try various loss functions and regulizers.\n",
    "\n",
    "It is also important to note that we are working here on a reduced dataset. Recall that our initial dataset was pretty large and that we found different ways to reduce the feature space without influencing too much the quality of our predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from src.regressions import *\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_FOLDER = \"data/\"\n",
    "X = np.load(DATA_FOLDER + \"ML/x_train.npy\")\n",
    "y = np.load(DATA_FOLDER + \"ML/y.npy\")\n",
    "\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X: (30349, 4500)\n",
      "y: (30349,)\n"
     ]
    }
   ],
   "source": [
    "print(\"X: \" + str(X.shape))\n",
    "print(\"y: \" + str(y.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Idea 1: Good ol' least squares (MSE loss without regularizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first try a standard and naive implementation of `least_squares` on the full dataset. This will give us a first raw result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss = 26.16344670274864\n"
     ]
    }
   ],
   "source": [
    "def run_least_squares():\n",
    "    w_star = least_squares(y, X)\n",
    "    loss = rmse(y, X, w_star)\n",
    "    print(\"Loss = \" + str(loss))\n",
    "\n",
    "run_least_squares()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the loss is quite large! We hope to do better by using polynomial expansion and 4-fold cross-validation to find the best `degree`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = 4\n",
    "k_indices = build_k_indices(y, k_fold)\n",
    "\n",
    "def run_least_squares_poly():\n",
    "    rmse_tr = []\n",
    "    rmse_te = []\n",
    "\n",
    "    for degree in range(4):\n",
    "        rmse_tr_tmp = []\n",
    "        rmse_te_tmp = []\n",
    "        for k in range(k_fold):\n",
    "            try:\n",
    "                loss_tr, loss_te, _ = cross_validation(y, X, k_indices, k, degree, least_squares)\n",
    "                rmse_tr_tmp.append(loss_tr)\n",
    "                rmse_te_tmp.append(loss_te)\n",
    "            except:\n",
    "                print(\"Could not inverse the matrix\")\n",
    "        rmse_tr.append(np.mean(rmse_tr_tmp))\n",
    "        rmse_te.append(np.mean(rmse_te_tmp))\n",
    "    pd.DataFrame([rmse_tr, rmse_te]).add_prefix(\"Degree \").rename({0: \"Train error\", 1: \"Test error\"}).head()\n",
    "    return rmse_tr, rmse_te"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([3.213171631581826,\n",
       "  0.40145325242460417,\n",
       "  0.2659636431823623,\n",
       "  0.1652126444248406],\n",
       " [3.2132414414047386,\n",
       "  0.5357537373798815,\n",
       "  0.6076663576371552,\n",
       "  0.853874297553556])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_least_squares_poly()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result above is composed of two arrays. The first one contains at index $i$ the __train error__ for polynomial expansion of degree $i$. Train error does not bring much useful information, hence we rather inspect the second array, representing the __test error__. It is constructed in an identical way.\n",
    "\n",
    "We see that polynomial expansion provides better results. In particular, the best `degree` is 1.\n",
    "\n",
    "__By adding a single feature of constant values, we manage to decrease the RMSE from 26 to 0.7!__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Idea 2: Ridge regression (MSE loss with $\\mathcal{L}_2$-regularizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since least squares optimize `MSE` as loss function, we decide to tackle the question of overitting immediately and consider adding to least squares a regularizer for the complexity of the model to be part of the optimization as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_fold = 4\n",
    "k_indices = build_k_indices(y, k_fold)\n",
    "\n",
    "def run_ridge_regression():\n",
    "    for lambda_ in np.logspace(-5, 1, 7):\n",
    "        print(\"Lambda = \" + str(lambda_))\n",
    "        for degree in range(1, 3):\n",
    "            print(\"    Degree = \" + str(degree))\n",
    "            rmse_tr_tmp = []\n",
    "            rmse_te_tmp = []\n",
    "            for k in range(k_fold):\n",
    "                ridge_lambda = lambda y, X: ridge_regression(y, X, lambda_)\n",
    "                loss_tr, loss_te, _ = cross_validation(y, X, k_indices, k, degree, ridge_lambda)\n",
    "                rmse_tr_tmp.append(loss_tr)\n",
    "                rmse_te_tmp.append(loss_te)\n",
    "            print(\"        \" + str(np.mean(rmse_tr_tmp)))\n",
    "            print(\"        \" + str(np.mean(rmse_te_tmp)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambda = 1e-05\n",
      "    Degree = 1\n",
      "        0.40145325242460456\n",
      "        0.5357537372320587\n",
      "    Degree = 2\n",
      "        0.26596380103399736\n",
      "        0.607666163440107\n",
      "Lambda = 0.0001\n",
      "    Degree = 1\n",
      "        0.4014532524246221\n",
      "        0.5357537359016747\n",
      "    Degree = 2\n",
      "        0.2659683198847004\n",
      "        0.6076760153324785\n",
      "Lambda = 0.001\n",
      "    Degree = 1\n",
      "        0.4014532524263871\n",
      "        0.535753722599882\n",
      "    Degree = 2\n",
      "        0.26597880149891057\n",
      "        0.6077023419396217\n",
      "Lambda = 0.01\n",
      "    Degree = 1\n",
      "        0.40145325260288683\n",
      "        0.5357535897864518\n",
      "    Degree = 2\n",
      "        0.26598240518404337\n",
      "        0.60770831648428\n",
      "Lambda = 0.1\n",
      "    Degree = 1\n",
      "        0.40145327025267574\n",
      "        0.535752282101457\n",
      "    Degree = 2\n",
      "        0.26605512986311564\n",
      "        0.6078706202481085\n",
      "Lambda = 1.0\n",
      "    Degree = 1\n",
      "        0.40145503506146574\n",
      "        0.5357412499056584\n",
      "    Degree = 2\n",
      "        0.26837515562810677\n",
      "        0.6137603953993617\n",
      "Lambda = 10.0\n",
      "    Degree = 1\n",
      "        0.4016313108003009\n",
      "        0.535835073442734\n",
      "    Degree = 2\n",
      "        0.28025537342649093\n",
      "        0.6431991791448415\n"
     ]
    }
   ],
   "source": [
    "run_ridge_regression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that all results are quite close to each other but if we pay attention to the details, we see that the best test error is achieved with a regulizer of $\\lambda = 1$ for `degree 1` and $\\lambda = 10^{-4}$ for `degree 2`. In all cases, `degree 1` remains our best degree value and we already reach __RMSE of 0.53__!\n",
    "\n",
    "Since the results are almost equal with and without a regulizer and provided that the test error increases for larger values of $\\lambda$, we deduce that the model is not overfitted and that a $\\mathcal{L}_2$ regulizer is useless..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Idea 3: Lasso (MSE loss with $\\mathcal{L}_1$-regularizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the section above, we saw that $\\mathcal{L}_2$-regularization did not prove useful. For completeness, we try $\\mathcal{L}_1$-regularization as well. Since there exists no analytical solution to the minimum of the `lasso` loss, we decide to use `mini-batch-SGD`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_lasso():\n",
    "    for lambda_ in np.concatenate([[0], np.logspace(-5, 2, 8)]):\n",
    "        print(\"lambda = \" + str(lambda_))\n",
    "        for degree in range(1, 3):\n",
    "            print(\"    degree = \" + str(degree))\n",
    "            rmse_tr_tmp = []\n",
    "            rmse_te_tmp = []\n",
    "    \n",
    "            for k in range(k_fold):\n",
    "                lasso_lambda = lambda y, X, w: lasso(y, X, w, lambda_)\n",
    "                lasso_stoch_grad_lambda = lambda y, X, w: lasso_stoch_grad(y, X, w, lambda_)\n",
    "                loss_lambda = lambda y, X: stochastic_gradient_descent(y, X, np.zeros(X.shape[1]), 100, 1e-1, lasso_lambda, lasso_stoch_grad_lambda, batch_size = 100)\n",
    "                loss_tr, loss_te, _ = cross_validation(y, X, k_indices, k, degree, loss_lambda)\n",
    "                rmse_tr_tmp.append(loss_tr)\n",
    "                rmse_te_tmp.append(loss_te)\n",
    "            print(\"        \" + str(np.mean(rmse_tr_tmp)))\n",
    "            print(\"        \" + str(np.mean(rmse_te_tmp)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we run it, we try a `stochastic_gradient_descent` using a `lasso`, just to see what we can expect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_ = 0\n",
    "lasso_lambda = lambda y, X, w: lasso(y, X, w, lambda_)\n",
    "lasso_stoch_grad_lambda = lambda y, X, w: lasso_stoch_grad(y, X, w, lambda_)\n",
    "_, losses = stochastic_gradient_descent(y, X, np.zeros(X.shape[1]), 80, 1e-1, lasso_lambda, lasso_stoch_grad_lambda, batch_size = 100, detail = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe502b84908>]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VHXW+PHPSScdSAIJCb3XAJEiAmLFyq6sK64NG+qC667P/tb1cfvq7vqs69q7smJD1q4oYqUpJTSlCMTQWwo1BFLP74+5g2MYkoCZ3EnmvF+veeXOLTOHZLhnvl1UFWOMMaamMLcDMMYYE5wsQRhjjPHLEoQxxhi/LEEYY4zxyxKEMcYYvyxBGGOM8csShDHGGL8sQRhjjPHLEoQxxhi/ItwO4IdISUnRjh07uh2GMcY0KUuXLi1S1dS6zmvSCaJjx47k5ua6HYYxxjQpIrK5PudZFZMxxhi/LEEYY4zxyxKEMcYYvyxBGGOM8csShDHGGL8sQRhjjPHLEoQxxhi/LEGchE/W7mb97oNuh2GMMQFlCeIE5RWUcOO0XC557AuWbt5zzPHZq3dx3oPzWLfLEogxpmmzBHGCHvl0A9ER4aQmRHPVs4tZmF8MgKryyKcbmPTCUtbuPMAz8/JdjtQYY34YSxAn4NvCEt5ZuYOrh3fg1UnDaJfcgolTFzN79S5ufWU5981ez7jsDMYPyuTdr3awv7TC7ZCNMeakWYI4AY98mkd0RDg3jupMWmIM0ycNo1NKPJNeWMrMr3fy2/N68sBl2Vx/WieOVFTz2rJtbodsjDEnzRJEPX1bWMLbK7Zz1fAOpMRHA9A6PppXbhzKpYMzee6aU7h5dBdEhN4ZiQxqn8xLizajqi5HbowxJ8cSRD098mkeURFhTBrV+Xv7k2Oj+OelAxjTM+17+68c1oH8wkN8+W1xY4ZpjDENxhJEPeQ7pYerh3c8Wnqoy/n90kmOjeTFRfWaVdcYY4KOJYh6eOzzb/2WHmoTExnOpYMzmb16NwUHjgQwOmOMCQxLEHWorlY+XrubC/pl1Lv04PWzoR2orFamL9kaoOiMMSZwLEHUYUNBCftKKxjepfUJX9spJY6R3VJ4edEWFm/cQ0VVdQAiNMaYwGjSS442hsUbPY3MQzq2OqnrbxndhYlTl/DTJ78kISaC07qmMLZvWy4ekIGINGSoxhjToCxB1GHRxj20TYwhq1WLk7r+1K4p5P7+LBZsKGLO+kI+X1fIB6t28eby7fzfT/qTlhDTwBEbY0zDsARRC1VlyaY9DO3U+gd920+MieS8fumc1y8dVeWFhZu5Z+Zaxj4wj79f0o9z+7RtwKiNMaZhBLQNQkQ2icjXIrJCRHKdfa86z1c4x1f4nH+niOSJyDoROTeQsdXHlj2l7D5QxpBOJ1e95I+IcPXwjsz8xWlkJMdw0wtL+ccH3zTY6xtjTENpjEbqMaqarao5AKp6mfM8G3gdeANARHoDE4A+wFjgMREJb4T4jmvRRs9srUMbMEF4dU1L4I1bRjAuO4On5+VTeLCswd/DGGN+CNd6MYmnzuanwCvOrnHAdFUtU9WNQB4wxK34ABZv3EOruCi6psUH5PWjIsK49YyuVFUrb6/YHpD3MMaYkxXoBKHAbBFZKiKTahwbCexW1Q3O83aA74CBbc4+1yzeuIdTOrYMaG+jrmkJDMhM4vVlliCMMcEl0AlihKoOAs4DJovIKJ9jl/Nd6QHA3134mJnuRGSSiOSKSG5hYWHDRutj5/7DbNlTypBOJz7+4USNH5zJ2p0HWL1jf8Dfyxhj6iugCUJVdzg/C4A3caqMRCQCuAR41ef0bUCWz/NMYIef13xKVXNUNSc1NTVQobM4gO0PNV3UP4PIcOH1pVaKMMYEj4AlCBGJE5EE7zZwDrDKOXwW8I2q+i6Y8A4wQUSiRaQT0A1YHKj46rJ44x7ioyPolZ4Y8PdqGRfFWb3a8PaK7Tba2hgTNAJZgmgDzBeRlXhu9DNVdZZzbALfr15CVVcDM4A1wCxgsqpWBTC+Wi3euIfBHVoSHtY4o53HD8qk+FA5c9YFrtrMGGNORMAGyqlqPjDgOMcmHmf/PcA9gYqpvopLythQUMKPBjZeG/noHqm0jovi9WXbOKt3m0Z7X2OMOR6brM+PJZv2Ao3T/uAVGR7GuOx2fLK2gH2l5Y32vsYYczyWIPzI3bSH6Igw+mUmNer7jh/cjvKqat5deUzbvDHGNDpLEH5s3lNKx9ZxREc07kDuPhlJ9MlI5Nn5GymvtMZqY4y7LEH4UVxSRkpClCvv/etze7CpuJTnv9jkyvsbY4yXJQg/ig+V0zruxFaPayhjeqRxeo9UHvpkA0UlNj+TMcY9liD8KC4pp3W8OyUIgN9d0JvDFVX8a/Z612IwxhhLEDUcqaiipKzyhNefbkhd0+K5angHXl2yhTU7DrgWhzEmtFmCqMFbrZPiYgkC4JdndiepRSR/eW81qsdMSWWMMQFnCaKG4hLPGAS32iC8kmIjuf3s7izM38PsNbtdjcUYE5osQdRQfMhTgnCzDcLr8iHtSUuI5r2vdrodijEmBFmCqKHIKUG42QbhFREeRk7HlizfstftUIwxIcgSRA3eNohgKEEADMxqyba9h21JUmNMo7MEUUNxSTmxUeHERgVsHsMTkt0+GYAVW/e5HIkxJtRYgqihuKQsaEoPAH0zkogIE1ZstWomY0zjsgRRg5ujqP1pERVOr/RElm+xEoQxpnFZgqihqKTc9TEQNWVnJbNy6z6qqm08hDGm8ViCqKGopCwoejD5Gtg+mUPlVeQVlLgdijEmhFiC8FFdrew55O48TP5kZ3kaqq27qzGmMVmC8LH/cAVV1RpUbRAAnVLiSGoRae0QxphGZQnCRzCNovYlIgxsn2xdXY0xjcoShI/Cg55R1KlB1gYBnmqm9QUHOXikwu1QjDEhwhKEj+9KEMGXIAa2b4kqfL1tv9uhGGNCREAThIhsEpGvRWSFiOT67L9VRNaJyGoR+T+f/XeKSJ5z7NxAxubP0Zlcg6yKCSA702motmomY0wjaYz5JMaoapH3iYiMAcYB/VW1TETSnP29gQlAHyAD+FhEuqtqVSPECHhGUYcJtIwNvgSRFBtJ59Q468lkjGk0blQx3QL8Q1XLAFS1wNk/DpiuqmWquhHIA4Y0ZmBFh8ppFRdFeJg05tvW28CslqzYus8WEDLGNIpAJwgFZovIUhGZ5OzrDowUkUUiMkdETnH2twO2+ly7zdn3PSIySURyRSS3sLCwQYMtOlgWdF1cfWW3T6aopJxtew+7HYoxJgQEuopphKrucKqRPhKRb5z3bAkMA04BZohIZ8Df1/Zjviqr6lPAUwA5OTkN+lW6OAgHyfka6AyYW7ZlL1mtYl2OxhjT3AW0BKGqO5yfBcCbeKqMtgFvqMdioBpIcfZn+VyeCewIZHw1eWZyDd4SRM+2CSTGRLAgr6juk40x5gcKWIIQkTgRSfBuA+cAq4C3gDOc/d2BKKAIeAeYICLRItIJ6AYsDlR8/hQH4UR9viLCwxjRNYV5G4qsHcIYE3CBrGJqA7wpIt73eVlVZ4lIFPCciKwCyoFr1HO3Wy0iM4A1QCUwuTF7MB2pqOJgWWXQTdRX06juqXywahffFpbQNS3B7XCMMc1YwBKEquYDA/zsLweuPM419wD3BCqm2hQfcsZAxAVvCQLgtK4pAMxZX2QJwhgTUDaS2lFcEryjqH1ltYqlc0oc8zY0bA8uY4ypyRKEwzuKOpjbILxGdU9lYX4xZZWNVgNnjAlBliAcRU4JItjbIABGdkvhSEU1uZtsVLUxJnAsQTiKgngeppqGdW5NZLgw16qZjDEBZAnCUVxSRovIcGKjGmN6qh8mLjqCwR1aMm+9jYcwxgSOJQhH8aFyUhKCv/TgNbJbKmt2HqDwYJnboRhjmilLEI6ikuCeh6mm0d1TAZifZ9VMxpjAsAThCPZR1DX1Tk+kVVyUVTMZYwLGEoSjqZUgwsKE07qmMHdDEdXVNu2GMabhWYIAqquVPU2sDQI84yGKSsr407urySs46HY4xphmJvi77DSCA0cqqKzWJlWCALigXzrzNhTyyuItTPtyM0M6tmL84HZ0aB1Hm8QY2iRGN4leWcaY4GR3D5rWGAhfLaLCeXDCQH5/YW9eX7qNVxZv4Y7Xv/7eOV3T4nnqqsF0To13KUpjTFNlCYLv5mFqCqOo/UmJj+am0V24cWRn8osOsfvAEXYfOMLO/Ud4bv5GfvLElzx7TQ4D27d0O1RjTBNiCYKmW4KoKSxM6JoWT9e070oL5/dL5+rnFvGzpxfx6BUDOaNnGxcjNMY0JdZIDRQfatoliNp0SonjjVtG0CUtjhunLeW1pdvcDskY00RYggD2lVYAkNwi0uVIAiM1IZrpk4YzvHNrfvPaSmav3uV2SMaYJsASBFBWWUVEmBAR3nx/HfHRETx19WD6tUviF9OXs2yLzQRrjKld870jnoCyimqiI5r/ryI2KoJnJ55CWkIMNzyfy8aiQ26HZIwJYs3/rlgP5VXVREeGux1Go0iJj+b564YAMHHq4qPrYBhjTE2WIAidEoRXp5Q4nr0mh90HjnDb9OWo2lQdxphjhc5dsRZllVVEhVCCABjYviV3XdCbBXnFvLViu9vhGGOCUEDviiKySUS+FpEVIpLr7PuTiGx39q0QkfN9zr9TRPJEZJ2InBvI2HyVVYZWCcLriiHtyc5K5u731rKvtNztcIwxQaYx7opjVDVbVXN89v3b2Zetqu8DiEhvYALQBxgLPCYijdIw4EkQodEG4SssTPj7Jf3Yd7iCf3zwjdvhGGOCTDB9bR4HTFfVMlXdCOQBQxrjjcsqq0KyBAHQKz2RG07rxPQlW1m8cY/b4Rhjgkig74oKzBaRpSIyyWf/FBH5SkSeExHvBEHtgK0+52xz9gVcWUU10ZGhmSAAbjurG+2SW3DXm19TXlntdjjGmCAR6LviCFUdBJwHTBaRUcDjQBcgG9gJ/Ms5V/xcf0z3GhGZJCK5IpJbWNgwy22GahWTV2xUBH/9UR82FJTw9Lx8t8MxxgSJgCYIVd3h/CwA3gSGqOpuVa1S1Wrgab6rRtoGZPlcngns8POaT6lqjqrmpKamNkic5SHaSO3rjJ5tOKd3Gx79LI/dB464HY4xJggE7K4oInEikuDdBs4BVolIus9pPwZWOdvvABNEJFpEOgHdgMWBis9XKHZz9eeuC3pRWaXcO8sarI0xgZ3uuw3wpoh43+dlVZ0lIi+ISDae6qNNwE0AqrpaRGYAa4BKYLKqVgUwvqNCtZtrTR1ax3H9yE48/vm3XDWsg60fYUyIC1iCUNV8YICf/VfVcs09wD2Biul4Qr0NwtfkMV15fek2/vTuGt685VTCwvw1DRljQoF9bQbKKkK3m2tN8dER3DG2Jyu37uPN5TbC2phQZndFnBJECHdzrenHA9sxICuZf8z6hpKySrfDMca4JOTvilXVSmW1WhWTj7Aw4Y8X9aaopIzLnvySLcWlbodkjHFByCcI78Awq2L6vkHtW/LM1Tls3VPKhQ/P45O1u90OyRjTyGq9K4rIGT7bnWocuyRQQTWmskpPRynr5nqsM3u14b1bR5LVKpbrn8/lvg/XUVVtU4MbEyrquive57P9eo1jv2vgWFxRdrQEYVVM/rRvHcvrt5zKZTlZPPJZHj97eiG79ttAOmNCQV0JQo6z7e95k1RWYVVMdYmJDOfen/TnX5cO4Ovt+zn/oXl89k2B22EZYwKsrruiHmfb3/MmyVvFZL2Y6jZ+cCbvTDmNtIRorv3PEv763hr2l1a4HZYxJkDqGijXWUTewVNa8G7jPO90/MuaDqtiOjFd0+J5a/II7p65hmfnb2T64i1cObwD15/WibSEGLfDM8Y0oLoSxDif7ftqHKv5vEk6WoKwKqZ6i4kM5+4f9eOKoR14/PNveXpuPlMXbGL8oHb8KLsdp3RsZSOwjWkGak0QqjrH97mIRAJ9ge3ODK1NnrVBnLxe6Yk8dPlAbj+7O0/M+ZY3l2/nlcVbaZsYw4X907liWAc6pcS5HaYx5iTV1c31CRHp42wnASuBacByEbm8EeILuLIqT4Kwbq4nr2NKHP8Y35+lvzubBydk0y8ziWlfbubSJ75g7yFb69qYpqquu+JIVV3tbF8LrFfVfsBg4DcBjayRfFeCsDaIHyouOoJx2e14+uoc3po8gn2lFfz53dV1X2iMCUp1JQjfr39nA28BqOqugEXUyKwXU2D0zkhkyhldeWvFDmavbjYfF2NCSl13xX0icqGIDARGALMARCQCaBHo4BpDmU21ETA/P70rvdITueutVewrtaomY5qauu6KNwFTgKnAL31KDmcCMwMZWGOxbq6BExURxn2X9mfvoXL+8u4at8MxxpygunoxrQfG+tn/IfBhoIJqTGUVVsUUSH0ykvj56V146NM80pNj6JuRRGbLWDJbtqBlXJTb4RljalFrghCRh2o7rqq/aNhwGp9VMQXelDO6sTB/D49+9u339p/Xty33/LgfrSxRGBOU6hoodzOwCpgB7KCZzL/kyzvdd1S4JYhAiYoIY8bNw9l/uIJte0vZtvcwK7fu45l5G1myaS73ju/Hmb3auB2mMaaGuhJEOnApcBlQCbwKvK6qewMdWGMpq6wmKiIMkWaX+4JOUotIklok0ScjiXP7tOWiARn86tUVXP98LhNOyeKPF/WhRZS1BRkTLGr92qyqxar6hKqOASYCycBqEbmqMYJrDGWVth61W3qlJ/L2lBHccnoXXs3dyqQXco92OzbGuK9ed0YRGQT8ErgS+ABYGsigGlNZZbX1YHJRdEQ4d4ztyb3j+zNvQxFTXl5OhTO63Rjjrrqm2viziCwFbgfmADmqer2q1qvPoohsEpGvRWSFiOTWOPZrEVERSXGei4g8JCJ5IvKVk5QCrqyi2koQQeCnOVn8+eI+fLRmN/8zY6WtXGdMEKirDeL3QD4wwHn8zamrF0BVtX893mOMqhb57hCRLDwjs7f47D4P6OY8hgKPOz8Dqqyyyrq4BolrTu1IaXkV9876hhaR4dzz475EWOcBY1xTV4II1JoP/8Yzl9PbPvvGAdNUVYGFIpIsIumqujNAMQBWxRRsbjm9C6XllTz8aR5rdx3g3vH96ZWe6HZYxoSkuhqpN/t7ANuA0+rx+grMFpGlIjIJQEQuxjNd+Moa57YDtvo83+bs+x4RmSQiuSKSW1hYWI8QaudJEPYtNZjcfnZ3Hr58INv3Huaih+dz34frOFJhjdfGNLa62iASReROEXlERM5x2gluxVPt9NN6vP4IVR2Ep/posoiMAu4C/uDv7fzsO6YiWlWfUtUcVc1JTU2tRwi1K6+ssqm+g4yIcNGADD6+fTQXZ2fwyGd5jH1gLi8u3ExpeaXb4RkTMuq6M74A9AC+Bm4AZgM/Acap6rjaLgRQ1R3OzwLgTWA0nmqrlSKyCcgElolIWzwlhiyfyzPxDM4LKCtBBK+WcVHc/9Nspl03hISYSH731iqG//1T/v7BWnbuP+x2eMY0e3WuSe2s/4CIPAMUAe1V9WBdLywicUCYqh50ts8B/qKqaT7nbMLTM6rIWe96iohMx9M4vT/Q7Q/g6cXUOs7aIILZqO6pjOyWwtLNe3luwUaenpvPK4u2MOPm4fRsa+0TxgRKXV+dK7wbqloFbKxPcnC0AeaLyEpgMTBTVWfVcv77eKqu8oCngZ/X831+EOvF1DSICDkdW/HYFYP5+PbRxEZFcPWzi9lSXOp2aMY0W3XdGQeIyAHncRDo790WkQO1Xaiq+ao6wHn0UdV7/JzT0dsFVj0mq2oXVe2nqrnHvmrDsyqmpqdzajwvXD+E8qpqrnpuEQUHj7gdkjHNUl29mMJVNdF5JKhqhM92syjbWzfXpqlbmwSmTjyFwoNlXPPcEvYfrqj7ImPMCQn5r85lFTYXU1M1sH1LnrhyMHkFB7lt+nI8Q2iMMQ0l5O+M5VVWxdSUjeqeyh1je/L5ukI+W1fgdjjGNCshfWdUVWuDaAauHt6RTilx3D1zrU30Z0wDCuk7Y0WVogrRkdYG0ZRFRYTxv+f3Ir/wEC8t3Ox2OMY0GyGdILxrD1gJouk7q1cap3ZpzQOfbGB/qTVYG9MQQvrOaOtRNx8iwu8u6M3+wxU8+MkGt8MxplkI6TvjdwnCqpiag94ZiVyWk8W0LzeRX1jidjjGNHmhnSCcGUJtJHXzcfs53YmJDOfSJ77kPws2Ul5pjdbGnKyQvjN6SxBRtihNs5GWEMP0ScPo1iaeP727hrPun8PbK7bbCnXGnISQvjN6v11aCaJ56dsuiVduHMZ/rj2FuOgIbpu+gqF/+4Q/vL2KRfnFliyMqae6ZnNt1qwNovkSEU7vkcaobqnMXrOLd1buYEbuVqZ9uZk2idHcO74/p/dIq/uFjAlhIZ4grJtrcxcWJoztm87YvukcKqvk028KePSzPG54Ppd/jO/PTwZnuh2iMUErpO+MZRVWggglcdERXDQgg//ePJyhnVvx6/+u5NHP8mwOJ2OOI7QThLVBhKSEmEimThzCuOwM/vnhOn731ir2Hip3Oyxjgo5VMWFVTKEoKiKMf/80m7aJMTw5N5//5m7j7N5tuDQnk5HdUgkP87dEujGhJcQThNPN1RJESAoLE+48vxfjstvx36VbeWv5dmZ+vZP2rWJ5cEI2A9u3dDtEY1wV0nfGcuvFZPCMwP7jRX1Y+L9n8tgVg6hW5adPfslz8zda+4QJaSGdIKyKyfiKjgjn/H7pzLx1JKO7p/KX99Zwy4vLOHDEJv8zoSmk74zf9WIK6V+DqSEpNpKnr87hrvN78dHa3Zz3wDzmri/0e27BgSNHv2gY09yE9J2xrLKa8DAhwqbaMDWICDeO6syMm4YTHRHG1c8t5n9mrGRfaTmqyoK8Im54fglD//4J4x//goKDR9wO2ZgGF9BGahHZBBwEqoBKVc0Rkb8C44BqoACYqKo7RESAB4HzgVJn/7JAxldWaetRm9oN7tCS928bycOfbuDJOfnMWV9Aq7go1u8uoXVcFFcN68B/c7cx/vEveP7aIXROjXc7ZGMaTGPcHceoaraq5jjP/6mq/VU1G3gP+IOz/zygm/OYBDwe6MBsuVFTHzGR4fy/c3vyzpTT6JwaT4vIcO67dAALfnsGfxnXl+mThlFaVsX4x79g2Za9bodrTINp9G6uqnrA52kc4O0mMg6Ypp5uIwtFJFlE0lV1Z6BiKauoti6upt56ZyQy46bhx+wfkJXM67ecyjVTF/OzpxcyflAmwzq3ZmjnVqQlxLgQqTENI9AJQoHZIqLAk6r6FICI3ANcDewHxjjntgO2+ly7zdkXuARRWWVdXE2D6JgSx+u3nMrv31rF2yt28NKiLQD0aJPA/ZcNoE9GkssRGnPiAv31eYSqDsJTfTRZREYBqOpdqpoFvARMcc71N3T1mE7oIjJJRHJFJLew0H/Pkvoqr7IqJtNwUuKjefzKwaz4w9m8PXkEd57Xk/2HK7h26hK27zvsdnjGnLCA3h1VdYfzswB4ExhS45SXgfHO9jYgy+dYJrDDz2s+pao5qpqTmpr6g+Irq6i2eZhMg4sID2NAVjI3je7C89cN4XB5FddOXcz+wzaewjQtAbs7ikiciCR4t4FzgFUi0s3ntIuBb5ztd4CrxWMYsD+Q7Q/gbaS2KiYTOD3aJvDkVYPZWHSIW15cakugmiYlkF+f2wDzRWQlsBiYqaqzgH+IyCoR+QpP0rjNOf99IB/IA54Gfh7A2ADr5moax6ldU7h3fH+++LaY37y2kiMVNrDONA0Ba6RW1XxggJ/94/2cjtN7aXKg4vGnrLKa+OiQnq/QNJJLBmWyY99h7pu9nvl5xdw8ujM/G9qe2CjP56+6Wtm8p5TkFpG0jItyOVpjPEL67mjdXE1jmnJGN07p2IqHPt3A3TPX8vjn33JmrzTyCw+xducBDpVXkRATwX+uHcLgDjaTrHFfSN8drZuraWxDO7fmpRuG8drNw+nTLokPVu1CgfGDM7nnx309o7OfXcSX3xa7HaoxoV2CKLeR1MYlOR1bMe26mp364OxebbjimUVMnLqYJ68azOk90lyIzhiPkL47llVaN1cTXNISY5g+aRhdUuO5cVouM78KaEc+Y2oV0ndH6+ZqglHr+GhemTSM/pnJTH55GX95d411jzWuCPEEYd1cTXBKahHJyzcO5doRHXluwUYufeILtu4pdTssE2JC9u5YVa1UVKmVIEzQio4I548X9eGJKweRX3SICx6ax/wNRW6HZUJIyCaIo+tRWxuECXJj+3qWQc1IbsGN03JZblOKm0YSsndH7zKRUbaanGkC2reOZdr1Q0hNiOa6/yzh28ISt0MyISBk745lVoIwTUxaQgzTrhtCeJhw9bOL2X3Aljk1gRWyd8ejVUzWBmGakI4pcUydOIR9peVc89xi1uw4gGeWmu8Ul5Txwpeb+NM7qzlwxGaQNScvZAfKeauYrBeTaWr6ZSbx5FU5XPf8Es5/aB5pCdGM6p5Kn4xE5qwvZN6GIqqqPUlj2Za9TLtuCMmxNr+TOXEhmyCOVHhLEJYgTNNzWrcU5t8xhjnrCpmzvpCP1uzmtaXbaJfcgkmjOjMuO4Ptew9zy4vLmPDUQl68YSgp8dFuh22amJBNEN+1QVgVk2ma0hJiuDQni0tzsqiqVrbvPUxmyxaEhXkWZ+zZNpFnJ+Zw47RcLnvyS166YRhtk2yNbFN/Ifv12aqYTHMSHia0bx17NDl4jeyWyvPXDmHX/iNc8tgCPly965g2C2OOJ2Tvjt4ShE33bZq7oZ1b88qkYcRFR3DTC0uZOHUJ+dZN1tRDyN4dy6wNwoSQ/pnJvH/bSH5/YW+Wbd7LuQ/M5f6P1lNZZXM8meML2btjeZV1czWhJTI8jOtP68Qnvx7NBf3SeeiTDUx4aiHb9x0+5lyrhjIQyo3UFdYGYUJTWkIMD0wYyJieadz15irOe2Au//eTAQzr3IpPvyngozW7mbu+kLjoCHpnJNIrPZE+GYmc1asNMdapI6SEboKwkdQmxI3Lbkd2VjK3vrI+RL2oAAASKklEQVScm19cSniYUFWtpCVEc9GADMqrqlm78yDzN+RTWa30bJvAgxMG0qNtgtuhm0ZiCcKqmEwI69A6jtduPpVn5udTWlbFWb3b0L9d0vd6Q5VVVvH5ukLuevNrLnpkPnee15OJp3ZERGp5ZdMchHCCsComY8DTk+/np3c97vHoiHDO7dOWwR1a8pvXvuLP767h47W7GdMjjfatYmnfOpYOreJoEWVftpqbgCYIEdkEHASqgEpVzRGRfwIXAeXAt8C1qrrPOf9O4Hrn/F+o6oeBis3bi8lmczWmflLio3n2mhxeXLSF+2evY0Fe8dFjMZGeJDNpVGdrp2hGGqMEMUZVfVc5+Qi4U1UrReRe4E7gDhHpDUwA+gAZwMci0l1VqwIRVFllNVHhYccMLDLGHJ+IcNWwDlw1rAP7SsvZsqeULXtKmfnVTu7/aD3/XbqV31/Qm7N7t7EqqGag0b8+q+psVa10ni4EMp3tccB0VS1T1Y1AHjAkUHHYcqPG/DDJsVH0z0zmwv4ZPH7lYF66YSjREeFMemEpVz+3mK+27XM7RPMDBfoOqcBsEVkqIpP8HL8O+MDZbgds9Tm2zdn3PSIySURyRSS3sLDwpAMrr6y2HkzGNKARXVP4wBmM9/X2/Vz8yAJueiGXdbsOuh2aOUmBvkOOUNVBwHnAZBEZ5T0gIncBlcBL3l1+rj9mtI6qPqWqOaqak5qaetKBlVVWWw8mYxqYdzDevN+M4VdndeeLvGLGPjiXO177ikNllXW/gAkqAU0QqrrD+VkAvIlTZSQi1wAXAlfod0M2twFZPpdnAjsCFZsnQVgJwphASIiJ5LazujH3N2O44bROzFi6lYsemc/anQfcDs2cgIDdIUUkTkQSvNvAOcAqERkL3AFcrKqlPpe8A0wQkWgR6QR0AxYHKr6yiiqbqM+YAGsZF8VdF/TmpRuGUnKkknGPLuDFhZttKo8mIpB3yDbAfBFZiedGP1NVZwGPAAnARyKyQkSeAFDV1cAMYA0wC5gcqB5MYCUIYxrTqV1SeP+2kQzr3JrfvbWKKS8vP+5yqGWVVZZAgkTAurmqaj4wwM/+447IUdV7gHsCFZMvTy8ma4MwprGkxEfzn4mn8OTcfO6bvY5VO/bz6M8G0bddEgCbig5x3+x1vPfVTs7smcafLu5DVqtYl6MObSE8krqa+OiQ/ecb44qwMOGW07twSseWTHl5OZc89gW/GduDLXtKeXnRFiLDwxg/KJNZq3Zy1v1zmDKmK5NGd7Yvcy4J2TtkeWU10XFWxWSMG3I6tuL920Zy+4wV3D1zLeFhwuVDsvjFGd1IS4zh/53bg7/OXMO/PlrP68u2ce2ITvxoYDuSWkS6HXpIkaZc15eTk6O5ubknde1Z98+hR5sEHr1iUANHZYypr+pq5f1VO+mdnkjn1Phjjs9dX8g/P1zH19v3ExMZxvn90rlqWAcGtm/pQrTNh4gsVdWcus4L2RKEjaQ2xn1hYcKF/TOOe3xU91RGdU9l1fb9vLJ4C2+v2MEby7Zzw2md+H9je1jVU4CF7B2yrMJGUhvTVPRtl8Q9P+7Hov89k2uGd+CZ+RsZ//gXtrZ2gIXsHdI7WZ8xpumIi47gz+P68tRVg9m29zAXPjyfZ+blk1dQYl1jAyC0q5hsWmJjmqRz+rSlX2YSv3rV08h998y1tIqLIqdDS0Z0TeHcPm1pmxTjdphNXkgmCFW1gXLGNHHpSS145cZh5BcdInfTHpZs2suSTXuYvWY3f3xnNQPbJ3Ne37b0zUgiKTaSlrFRtIyNsoWNTkBIJojKakXVVpMzpqkTEbqkxtMlNZ7LTmkPQF5BCbNW7eSDVbv42/vf1Dgffnlmd247q5sb4TY5IZkgbD1qY5qvrmnxTDmjG1PO6Mb2fYfZUlzK/sPl7CutYO6GQv798XriosO5YWRnt0MNeqGZICqc9aitF5MxzVq75Ba0S25x9PmlOVnAMu6euZakFpHOc3M8oZkgjpYgLEEYE0rCw4R/X5bNwSO53PH6VyS2iGRUt1Tyi0rYWHSI0rIqLs7OsHW1HSGeIOxDYEyoiY4I54krB3Pls4u45cWlVNfoHfv0vHz+fVn20UkEQ1mIJghPFZOtB2FMaIqLjmCqM7NsXFQ4nVLi6Zwax679R/jtG1/xo0cXcNuZ3bjl9C6UV1WzfncJ63YdYPvew+wtrWBvaTn7D1fQKz2Rn5/eheTYKLf/SQERmgmiwqqYjAl1ybFR3DG25/f29UpPZPYvR/P7t1fxr4/W8+yCjew/XIF3DJ4IJLXwdJlNiIngmXn5zMjdyu1nd+dnQ9oT0cwG34ZmgrAqJmPMcSTFRvLQ5QM5p08bPl6zm04p8fRom0DPtglktYolPEyOnrt25wH++t4a/vD2al74cjN3XdCL03ukuRh9wwrJBFHuTRDWi8kYcxwX9s+odSJB8JQ4XrphKB+vLeCemWuYOHUJo7un8rsLetGtTUIjRRo4IZkgvG0QVsVkjPmhRISze7dhdPdUpn25iQc/2cDYB+dx+ZAszuzVhoykFqQnx5AY0/TWsgjRBGFVTMaYhhUVEcYNIztzyaBMHvh4PS8t2sKLC7ccPZ4YE8GQTq0Y1T2Vkd1S6dg6FhGp5RXdF5ILBuUVlPDRmt1MOCWLlnHNs/eBMcZdRSVlbC4+xI59R9i5/zAbiw6xIK+YLXtKAWjfKpafDM7kpzlZjT6xYH0XDArJBGGMMW7ZXHyIuRuKmLVqJwvyigkTOKNnGpfmZHFa1xTiogNfsRMUCUJENgEHgSqgUlVzRORS4E9AL2CIqub6nH8ncL1z/i9U9cPaXt8ShDGmKdtcfIhXl2xlRu42ikrKiAwXTunYitHdU+nQOo7Cg0fYfaCMgoNHSIyJZEBWMtlZyWS2bPGDqqeCKUHkqGqRz75eQDXwJPBrb4IQkd7AK8AQIAP4GOiuqlXHe31LEMaY5qCiqprFG/cwZ30hc9YVsm73waPHwgRS4qPZd7jiaA/MlPgobh7d5aQnHAzaNalVdS3gL/uNA6arahmwUUTy8CSLLxs3QmOMaVyR4WGM6JrCiK4p/O/5vdi5/zDFJeWkJUbTOi6a8DChoqqadbsOsnzrPlZu3UdaYuDbLQKdIBSYLSIKPKmqT9Vybjtgoc/zbc4+Y4wJKelJLUhPavG9fZHhYfRtl0TfdklcNaxDo8QR6AQxQlV3iEga8JGIfKOqc49zrr8KtWPqv0RkEjAJoH379g0XqTHGmO8J6EgxVd3h/CwA3sRTZXQ82wDfydkzgR1+XvMpVc1R1ZzU1NSGDNcYY4yPgCUIEYkTkQTvNnAOsKqWS94BJohItIh0AroBiwMVnzHGmNoFsoqpDfCm0xgdAbysqrNE5MfAw0AqMFNEVqjquaq6WkRmAGuASmBybT2YjDHGBJYNlDPGmBBT326uNludMcYYvyxBGGOM8csShDHGGL+adBuEiBQCm0/y8hSgqM6z3BGssQVrXGCxnYxgjQuCN7ZgjQtOLLYOqlrnOIEmnSB+CBHJrU8jjRuCNbZgjQsstpMRrHFB8MYWrHFBYGKzKiZjjDF+WYIwxhjjVygniNomDnRbsMYWrHGBxXYygjUuCN7YgjUuCEBsIdsGYYwxpnahXIIwxhhTi5BMECIyVkTWiUieiPzW5VieE5ECEVnls6+ViHwkIhucny1diCtLRD4TkbUislpEbguG2EQkRkQWi8hKJ64/O/s7icgiJ65XRSSqMeOqEWO4iCwXkfeCKTYR2SQiX4vIChHxruQYDJ+1ZBF5TUS+cT5vw4Mkrh7O78r7OCAivwyS2H7lfP5Xicgrzv+LBv+chVyCEJFw4FHgPKA3cLmz3Klb/gOMrbHvt8AnqtoN+MR53tgqgf9R1V7AMGCy83tyO7Yy4AxVHQBkA2NFZBhwL/BvJ669eNY2d8ttwFqf58EU2xhVzfbpDun23xPgQWCWqvYEBuD53bkel6quc35X2cBgoBTPsgWuxiYi7YBf4FnOuS8QDkwgEJ8zVQ2pBzAc+NDn+Z3AnS7H1BFY5fN8HZDubKcD64Lg9/Y2cHYwxQbEAsuAoXgGCEX4+xs3ckyZeG4aZwDv4VkIK1hi2wSk1Njn6t8TSAQ24rSHBktcfuI8B1gQDLHhWWlzK9AKz0zZ7wHnBuJzFnIlCL775XoF49KmbVR1J4DzM83NYESkIzAQWEQQxOZU4awACoCPgG+Bfapa6Zzi5t/0AeA3QLXzvDXBE5t3CeClzsqM4P7fszNQCEx1quWecdaPcTuumiYArzjbrsamqtuB+4AtwE5gP7CUAHzOQjFB1GtpU+MhIvHA68AvVfWA2/EAqGqVeor9mXhWKezl77TGjQpE5EKgQFWX+u72c6pbn7cRqjoIT/XqZBEZ5VIcviKAQcDjqjoQOIQ71VzH5dTlXwz81+1YAJw2j3FAJyADiMPzN63pB3/OQjFB1GtpU5ftFpF0AOdngRtBiEgknuTwkqq+EUyxAajqPuBzPG0kySLiXQDLrb/pCOBiEdkETMdTzfRAkMSG+l8C2O2/5zZgm6oucp6/hidhuB2Xr/OAZaq623nudmxnARtVtVBVK4A3gFMJwOcsFBPEEqCb0+Ifhafo+I7LMdX0DnCNs30Nnvr/RiUiAjwLrFXV+4MlNhFJFZFkZ7sFnv8sa4HPgJ+4FReAqt6pqpmq2hHP5+pTVb0iGGKT4y8B7OrfU1V3AVtFpIez60w8q0q6/n/Ax+V8V70E7se2BRgmIrHO/1Pv76zhP2duNvy49QDOB9bjqbu+y+VYXsFTj1iB59vU9XjqrT8BNjg/W7kQ12l4iqhfASucx/luxwb0B5Y7ca0C/uDs74xnDfM8PFUB0S7/XU8H3guW2JwYVjqP1d7Pvdt/TyeGbCDX+Zu+BbQMhric2GKBYiDJZ5/rsQF/Br5x/g+8AEQH4nNmI6mNMcb4FYpVTMYYY+rBEoQxxhi/LEEYY4zxyxKEMcYYvyxBGGOM8csShGk2ROQuZ4bLr5zZN4c6+yNE5G/OLJfemTnv8rmuytm32pkl9nYROeb/hohkiMhrzna2iJzfgLEni8jP/b2XMW6xBGGaBREZDlwIDFLV/ngG0Hnn3Lobz5QE/dQzRcdIINLn8sPqmbWzD54JCc8H/ljzPVR1h6p6ByJlO+edSIwRtRxOBo4miBrvZYwrLEGY5iIdKFLVMgBVLVLVHSISC9wI3KqqR5xjB1X1T/5eRD3TUEwCpjijVI8SkY7O/PtRwF+Ay5ySx2XOSOXnRGSJM+ncOOeaiSLyXxF5F89EefEi8omILBPP2gzjnJf/B9DFeb1/et/LeY0YEZnqnL9cRMb4vPYbIjLLKR39n7M/XET+48T6tYj8qgF/zyaE1PaNxpimZDbwBxFZD3wMvKqqc4CuwBZVPVjfF1LVfKeKKQ3Y7ed4uYj8Ac98/FMARORveKbWuM6ZCmSxiHzsXDIc6K+qe5xSxI9V9YCIpAALReQdPBPU9XVKON4ZdL0mO+/bT0R64kk03Z1j2Xhm2i0D1onIw07c7dSzVgDeqUmMOVFWgjDNgqqW4FnUZRKe6aNfFZGJNc8TkWudb+lbRSSr5nHfU08whHOA3zrTkH8OxADtnWMfqeoen9f9m4h8hSeRtQPa1PHap+GZTgFV/QbYDHgTxCequt8pHa0BOgD5QGcReVhExgJBMQuvaXosQZhmQz3TgH+uqn8EpgDj8cxL0947UZ2qTnW+pe/HsxLXMUSkM1DFic3SKcB4py0jW1Xbq6p3VblDPuddAaQCg504duNJJnW99vGU+WxX4VkwZi+eldk+x1P6eKb+/wxjvmMJwjQL4lk/uJvPrmxgs6qW4pmV9hERiXHODQf8rtcrIqnAE8AjWvtEZQeBBJ/nHwK3etstRGTgca5LwrNmRIXTltDhOK/nay6exIJTtdQez6pmfjlVV2Gq+jrwezzTZxtzwqwNwjQX8cDDTn17JZ6Sg3fVtLuAvwKrROQgcBh4nu/my2/hVA1FOte+APhOce7PZ3xXpfR35/UfAL5yksQmPL2qanoJeFdEcvHMkPsNgKoWi8gCp2H6Azzrpns9BjwhIl878U1U1bIabei+2uFZoc37BfDOOv4txvhls7kaY4zxy6qYjDHG+GUJwhhjjF+WIIwxxvhlCcIYY4xfliCMMcb4ZQnCGGOMX5YgjDHG+GUJwhhjjF//H885EwWiZtdSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel(\"SGD iterations\")\n",
    "plt.ylabel(\"RMSE\")\n",
    "plt.plot(losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "79\n",
      "510.5618316456103\n"
     ]
    }
   ],
   "source": [
    "print(np.argmin(losses))\n",
    "print(losses[np.argmin(losses)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the result, we see that convergence is really slow. Indeed, since $\\lambda = 0$, we are basically optimizing a standard `MSE` for which we know (from above) that the optimal solution is 26. We deduce it is unlikely to do better than `ridge_regression`, and we drop this idea."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Idea 4: MAE loss (with SGD) with $\\mathcal{L}_2$-regularizer "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the lectures, we also talked about `MAE` loss, which we decide to try here. From the lab, we know that `MSE` is well-suited for out task, but we decide to implement `MAE` as well and hope to be surprised!\n",
    "\n",
    "Again, the minimum of `MAE` can not be found using any kind of closed-form formula, and we naturally decide to use `mini-batch-SGD` one more time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_mae_sgd():\n",
    "    for lambda_ in [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]:\n",
    "        print(\"lambda = \" + str(lambda_))\n",
    "        for degree in range(1, 3):\n",
    "            print(\"    degree = \" + str(degree))\n",
    "            rmse_tr_tmp = []\n",
    "            rmse_te_tmp = []\n",
    "            \n",
    "            lambda_mae = lambda y, X, w: mae(y, X, w, lambda_)\n",
    "            lambda_mae_stoch_grad = lambda y, X, w: mae_stoch_grad(y, X, w, lambda_)\n",
    "            for k in range(k_fold):\n",
    "                mae_lambda = lambda y, X: stochastic_gradient_descent(y, X, np.zeros(X.shape[1]), 50, 1, lambda_mae, lambda_mae_stoch_grad, batch_size = 100)\n",
    "                loss_tr, loss_te, _ = cross_validation(y, X, k_indices, k, degree, mae_lambda)\n",
    "                rmse_tr_tmp.append(loss_tr)\n",
    "                rmse_te_tmp.append(loss_te)\n",
    "            print(\"        \" + str(np.mean(rmse_tr_tmp)))\n",
    "            print(\"        \" + str(np.mean(rmse_te_tmp)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we first try to validate the method before we run that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Degree 1\n",
      "    loss = 26.159135164285743\n",
      "Degree 2\n",
      "    loss = 23.871869495951007\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEKCAYAAAAyx7/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl4VdW9//H3NxAI8xxkDqPITEBEBbVYBsWKViniUGm14oC1rXqrv2trr0LV2l/rdSqtdQaBKhbBOoAo1qqoBBnCFEBBQYWAyCQQQtb9Y+2QgIQwnJN99jmf1/Ps55yzz/TdMX5YWXuttc05h4iIREda2AWIiMjRUXCLiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiKkcjw9t2LChy8rKisdHi4gkrZycnE3OuUblvS4uwZ2VlcW8efPi8dEiIknLzNYeyevUVSIiEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxMRlHPexWvHVdv616Av/wGz//uJ7pXYdwCjjicO+5yhfX9YTMXS0X3G44473dx/Td8T/K8r+7jC/vAyx/O9XERLxZ5iIujWvS5/W9eP6HQkV3Ks27uCht1ahy2CKSFRde2bb1Aruod2aMLTb0EM+V9ZFjQ8X8mU9VeZnlfn6sr+j7O8+ujeF+Y9VRXz30f48YvrdCdgQSMCSDivZLyoey6OrUin+PdAJFdyHU1ZXxbH9+aa/+UQkunRyUkQkYhTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMQpuEZGIUXCLiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiFFwi4hEjIJbRCRiFNwiIhGj4BYRiRgFt4hIxCi4RUQiRsEtIhIxCm4RkYhRcIuIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMQpuEZGIUXCLiESMgltEJGIU3CIiEaPgFhGJGAW3iEjEKLhFRCJGwS0iEjEKbhGRiEm84HYu7ApERBJaYgV33uvwQFfYkR92JSIiCSuxgrtmY9j6OayaFXYlIiIJK7GCu0l3qHmCb3mLiMghJVZwm0H7gbD6Tdi3N+xqREQSUmIFN0CHwbBnG3w2N+xKREQSUuIFd5uzIC0dVqq7RETkUBIvuKvWgqzTIW9m2JWIiCSkxAtugPaDYdMK2LIm7EpERBJOYgZ3h8H+Vq1uEZHvSMzgbtAW6reFvNfCrkREJOEkZnADdBgCa/4DBTvDrkREJKEkcHAPgn174JO3w65ERCShJG5wtzwNqtTSsEARkYMkbnBXrgJtz4KVs7RioIhIKYkb3OCHBW5bDxtyw65ERCRhJHhwD/K3WnRKRGS/xA7uWo2hSQ9YqfHcIiLFEju4wU/G+fxD2Lk57EpERBJC4gd3+8GAg1VvhF2JiEhCSPzgbtoTajTSsEARkUDiB3daGrQb6Fvc+wrDrkZEJHSJH9zgZ1Hu3gqffxB2JSIioYtGcLc9G9Krw8JJYVciIhK6aAR3Rm3oOhwWvwC7toRdjYhIqKIR3AAnXw2Fu2DBc2FXIiISqugEd5Nu0OIU+OhxKCoKuxoRkdBEJ7jBt7q/Xg2fzgm7EhGR0EQruDsNg+oNfKtbRCRFRSu4K1eF7B/Dildg67qwqxERCUW0ghug10/8+tw5T4VdiYhIKKIX3PVa+etR5jwNhQVhVyMiUuGiF9zgT1Lu3AjLpoddiYhIhYtmcLcdAPWydJJSRFJSNIM7LQ16XwWfvQcbloRdjYhIhYpmcAP0vBwqZ6jVLSIpJ7rBXb0+dLkIFk2B3dvCrkZEpMJEN7gBTr4KCnbAwslhVyIiUmGiHdzNekGLvvDO/4c9O8KuRkSkQkQ7uAEG3Q07voL3Hgq7EhGRChH94G7RBzpfCO89CNu+DLsaEZG4i35wA3z/d1BUCG+ODbsSEZG4S47grpcFp4yGBRPhy0VhVyMiElfJEdwA/W+BanVh5h1+ESoRkSSVPMFdrS6ceRt8+jasnBl2NSIicXPY4DazAaXutz7ouR/Gq6hj1vunUL8tzPwN7CsMuxoRkbgor8X9x1L3px703B0xruX4Va4CA++CTStg/lNhVyMiEhflBbeVcf9QjxNDx6HQ6nR46x7YvTXsakREYq684HZl3D/U48RgBoPGwreb/IxKEZEkU7mc59uY2XR867r4PsHj1mW/LWTNsqH7SHj/Eej6IzihS9gViYjEjLnDDJ0zszMP92bn3NuH2t+7d283b9684yztOH37NTzSB2o3hatnQ6X0cOsRESmHmeU453qX97rDdpU4594uvQHvAduAZWWFdsKoXh/O/SN8udBPhxcRSRLlDQccb2adg/t1gIXAM8DHZjayAuo7Pp0vgE7DYM69sHF52NWIiMREeScn+zvniq8N9hMgzznXFegF/FdcK4uVc/8IVWrCSzdA0b6wqxEROW7lBXdBqfsDgWkAzrmv4lZRrNXMhHP+AOvnwdxHw65GROS4lRfc35jZeWbWEzgdeA3AzCoD1eJdXMx0vRhOPNevHrhpVdjViIgcl/KCezQwBngS+EWplvbZwL/iWVhMmcF5f4bKVWH6GCgqCrsiEZFjVt6okjzn3BDnXA/n3FOl9r/unLs57tXFUq0TYMi98Nn78NFjYVcjInLMDjsBx8wOO47OOffz2JYTZ91HQu6LMOtOaNrTXz1HRCRiyusquRboB3wBzANyDtqixQwu+IuflDNxOGxYGnZFIiJHrbzgbgL8DRgMXAGkAy855552zj0d7+LiomYjuOKfkF4Nnr0QtqwJuyIRkaNSXh/3ZufceOfc9/DjuOsCS83sigqpLl7qtfLhXbjbh/eOjWFXJCJyxI7oCjhmlg3cBFwOvEoUu0kOlnkSXPYCbP8KJvxQS8CKSGSUN+X9LjPLAX4FvA30ds5d5ZxLjs7hFifDiGf9dPjnLoG9u8KuSESkXOW1uO/Ad490B+4B5pvZIjNbbGbJcTn1dt+HH/7VDxN8fhTs2xt2RSIih1XeetyJu+Z2LHW5CHZtgX/dDJMugR89A1VqhF2ViMghHTa4nXNrD7XfzNKAkcAhn4+kk6+GtHR4+Rfw9A/g0uehRoOwqxIR+Y7y+rhrm9ntZvawmQ0y70bgE+BHFVNiBep1JYyYABuWwBOD4ZvPwq5IROQ7yuvjfhY4EVgMXA28BVwMXOCcGxbn2sLRcagfKrhjIzw+GDYuC7siEZEDlBfcbZxzo5xzf8V3jXQCBjvnFsS/tBC1Og1+8gq4It/y/mxu2BWJiOxXXnDvH2LhnNsHrHPO7Y5vSQnihC5w1Uyo3hCeGQYrXg27IhERoPzg7m5m24JtO9Ct+L6ZbauIAkNVr5UP78yTYPJlMP+ZsCsSESl3ynsl51ztYKvlnKtc6n7tiioyVDUawpUvQ5uzYPqN8Pb94FzYVYlICjuiKe8pr2pNuHQKdLsE3hoLr9yi61eKSGjKm4AjxSql+yVhazWGd/8XdmyAH/4d0jPCrkxEUoxa3EcjLQ0G3gWD74FlM/ziVLu+CbsqEUkxCu5jcer1cNHj8PmH8PfvQ35e2BWJSApRcB+rrhfDj6f5NU4eGwDLXwm7IhFJEQru45HVD66ZAw3awuSRMOdeXUFeROJOwX286raAn77mL0Q85x6YcjnsTv4h7iISHgV3LKRX8yNOhtwHea/B38+GTSvDrkpEkpSCO1bMoO+1cOV0+PZr+NtZsOA5TdYRkZhTcMdaVj8Y/W9o0gOmXQdTr9KQQRGJKQV3PNRp5lveA34DS6bB+P6w9v2wqxKRJKHgjpe0SnDGLX6RqrQ0eOpceOse2FcYdmUiEnEK7nhr3huu/Y9f5+Tte+HJcyB/RdhViUiEKbgrQtVacOFf/GzLTXkwvp9vfRfuCbsyEYkgBXdF6noxjJkHnYb51vf4frDm3bCrEpGIUXBXtJqN4KK/w+VToXC37/ue/nM/dV5E5AgouMPS7vtw/Vw47efw8QR4uA8smKQp8yJSLgV3mKrUgEF3+/VO6raEadfCE4Ng/fywKxORBKbgTgRNusFVs/y0+S1r/WqDL42BHflhVyYiCUjBnSjS0qDHpXBjDpw2BhZOgod6wfuPwr69YVcnIglEwZ1oMmrDoLFw3ft+DPjrt8Ojp8KK17TuiYgACu7E1aiDH3kycop/PGkEPHsBfJUbbl0iEjoFdyIzgxOHwPXvwzl/gC8X+rHf02+E7RvCrk5EQqLgjoJK6XDKaPj5x9D3ej9s8KFsePsPsGd72NWJSAVTcEdJtXow5PdwwwfQ5ix4axw80A3efRAKvg27OhGpIAruKGrQFi6ZCD97E5r2hFm/gQd7woePaf0TkRSg4I6yZr3gihfhJ6/6MH/lFj+EcP4zGkIoksQU3Mmg1Wkw6l9wxT+hZqY/eflwbz+VXut/iyQdBXeyMIO2A+Dq2X4IYUYdeOkGH+ALnlOAiyQRBXeyKR5CeM3bcMkkqFrTX/vykT6wcLICXCQJKLiTlRl0PBdGvwMjJkB6NfjnaH8Sc+54KNgZdoUicowU3MnODE76gQ/wSyZB7abw2q/hz53hzXFayEokghTcqSItzbfAr3odfvo6tDwN/v0HeKALvPwr2LA07ApF5AhVDrsACUHLvn7Lz4P3HoSPn4V5j/vhhT2vgC4X+cWuRCQhmYvDinO9e/d28+bNi/nnSpzs3ASLpsD8ZyF/GaRXh04XQPYV0PJU390iInFnZjnOud7lvk7BLfs5B+tz/ASe3KlQsAMatIPsH0P3kX6MuIjEjYJbjk/BTlgyzXejfPY+pFWGDkMg+0podzakVQq7QpGkc6TBrT5uObQqNaDnZX7Lz4OPn/GrEi5/GWo19Vfr6Xk51G8ddqUiKUctbjlyhQWQ95rvSlk9G1wRtD4Dev4YTjrPjxUXkWOmFrfEXuUq0Ol8v21d76fSf/wsvHi1n2Lf9Ue+Fd6ku05oisSRWtxyfIqKYM07PsCXTod9e6BhBx/iXS+C+m3CrlAkMnRyUireri3+hObi52Htu35f85Oh63Do/EOo2Sjc+kQSnIJbwrV1HSx+wW8bFoNVgjZnQpeLoeNQqFY37ApFEo6CWxLHxmW+FZ47FbasgUpVoP0gP0OzwxCoUj3sCkUSgoJbEo9zsH4+5L4AuS/Cjq8gvQa0H+hb4e0H+utqiqQoBbcktqJ9vh8890VY8Qrs2OAn+WT1gxOH+gWx6jQPu0qRCqXglugoKvJT7Ze/DMv/BZtX+v1NevglaU/6ATQ6MdwaRSqAgluiKz+vJMTXB79HDTtAx/N8iDftqXHikpQU3JIctq73XSnLpsOad8Ht81Pu2w+EDoOh9Zn+8mwiSUDBLcnn269hxat+2v3qt6Bgux+h0uo0aD/Yj1Rp0FatcYksBbckt8IC+HwurJwJeTNh0wq/v15rH+DtB/oTnVo/RSJEwS2pZcsaWDnLb5/+Gwp3QeUMvwhWu4F+Kdr6bdQal4SmRaYktdTLgj4/89veXX6o4cpZvkW+cqZ/Td1WPsDbDvCBnlEn1JJFjpVa3JL8Nq+G1W/67dN/+yv7WCW/jkq7s/3WpKe/oLJIiNRVInIo+/bC5x8GQT4bvlgAOKjewLfE2wYt8lqNw65UUlDCBffevXtZt24du3fvjvn3CWRkZNC8eXPS09PDLiVadm7yI1RWveGDfGe+35/ZCVr29RdLbtkX6rRQ/7jEXcIF96effkqtWrVo0KABpv8BYso5x+bNm9m+fTutW+tSYsesqMivZLjqDT9mfN1HsGebf652M2hxih96mNUPGnVUkEvMJdzJyd27d5OVlaXQjgMzo0GDBuTn54ddSrSlpfmr9zTpDv1v9uupbFwKn831F0z+bC4sedG/tnpDyDodsvoryKXCVeioEoV2/OhnGwdpleCErn7r8zO/b8taWPOfYHsHlr7k91dv4FvjrU73W+MuOtkpcZOywwF/97vfUbNmTW655ZawSynX5s2bufjii/noo48YNWoUDz/8cNglpa56rfzW8zL/uDjI177rb5fN8Psz6vj+8Van+0Bv0h0q6fyDxEbKBnesFBYWUrlyfH+MGRkZ3H333eTm5pKbmxvX75KjdHCQb13n+8fXBlvea35/enVo3htangatTvVDEavUCK9uibSU+ltu3LhxdOjQgX79+rFixYr9+1evXs2QIUPo1asX/fv3Z/ny5fv39+3bl65du3LHHXdQs6ZfzGjOnDn079+f888/n06dOgEwYcIE+vTpQ48ePRg9ejT79u0DYObMmZx66qlkZ2czfPhwduzYcdR116hRg379+pGRkXG8PwKJtzrNofsIOP9BuDEHbs6D4U9Bzyv8NTnfvg+eGQb3toS/ngEv/woWPOdXRCwqCrt6iYhQWtz/M2MJS7/YFtPP7NS0Nnf+oHOZz+fk5DB58mQWLFhAYWEh2dnZ9OrVC4BrrrmG8ePH0759ez744AOuv/563nzzTW666SZuuukmRo4cyfjx4w/4vPnz55Obm0vr1q1ZtmwZU6ZM4d133yU9PZ3rr7+eiRMncu655zJ27FjeeOMNatSowX333cef/vQnfvvb3x7wWffffz8TJ078Ts1nnHEGDz74YAx+OhKaWo2h84V+A9i91Y8jX/ueX7J20T9g3uP+uap1oFk2NOvll65t2sOPZtH5CzlIynSVvPPOO1x44YVUr+6vb3j++ecDsGPHDt577z2GDx++/7V79uwB4P3332fatGkAXHrppQf0h/fp02f/0LvZs2eTk5PDySefDMCuXbvIzMxk7ty5LF26lNNPPx2AgoICTj311O/Uduutt3LrrbfG+pAlEWXU8QtgtR/oHxftg00rfYivm+dv//Nnv3wtQI1GQYgHW7NeUDMzvPolIYQS3IdrGVe0oqIi6taty4IFC47qfTVqlPRPOue48sorueeeew54zYwZMxg4cCCTJk067GepxZ3C0ipBZke/9bzc79u7CzYsgS8+LtlWvQEu6Eqp07KkZd6sl2+Zq788paRMi/uMM85g1KhR3H777RQWFjJjxgxGjx5N7dq1ad26Nc8//zzDhw/HOceiRYvo3r07ffv2ZerUqYwYMYLJkyeX+dlnn302w4YN45e//CWZmZl8/fXXbN++nb59+3LDDTewatUq2rVrx86dO1m/fj0dOnQ44P1qccsB0qv5E5nNS83DKNgJXy70F1ten+Nb5kv9X4NYmh9H3jTbh3izbD8csXLVcOqXuEuZ4M7OzmbEiBF0796dzMzM/d0aABMnTuS6665j7Nix7N27l0suuYTu3bvzwAMPcPnllzNu3DiGDBlCnTqHXk2uU6dOjB07lkGDBlFUVER6ejqPPPIIffv25amnnmLkyJH7u1/Gjh37neA+EllZWWzbto2CggKmTZvGzJkz958YlRRQpUYwTvy0kn078uGL4iCfD3mvwoIJ/rm0dGjcOZhQ1A1O6O4fV6keTv0SUxU25X3ZsmWcdNJJMf+uePr222+pVq0aZsbkyZOZNGkSL730UthllSmKP2OJIef8cMQv5vvulfXz4atFfjQL+JZ5g/Y+yJv2hGa9/X1dbCJhJNyU9yjKyclhzJgxOOeoW7cuTzzxRNgliZTNDOq28FunYX5fcZh/tQi+XORv174Hi5/3z6dV9i3xZr1LRrM0bK/JQglOwX0Y/fv3Z+HChWGXIXLsSod5x6El+7d/FXSx5PjRLIufLxmWWKkKNDrR95M37hxsXTSaJYEouEVSUa0TfJAXh3lREWxe6VvlG3L99skcWFhqRFStJnBCt5J+8ybdtdxtSBTcIuIXxGp0ot8omdPAzs0+xL9aXNLdsmpWydDEjLolrfLMTsHtSVC1ViiHkSoU3CJSthoNoM2Zfiu2dxdsWApfLvCBvnEpLJgEBdtLXlO3VUmIFwd6g3bqO48RBbeIHJ30atC8l9+KOQfffOYnDm1c4oN941LIe71kFmhaOjTsAI07+TDP7OTvq7vlqKVscEdpWddZs2Zx2223UVBQQJUqVbj//vsZMGBA2GWJlDArWSmx47kl+wv3wKY82LgsCPXgwhTFo1oAqtQKWuYdfZg36ugf12ysQC9DygZ3rFTEsq4NGzZkxowZNG3alNzcXAYPHsz69evj+p0iMVG5asnFKErbvRU2Lj+wdb7sZZj/TMlrMuoGgV7c3dLFt9AzDj0RLpWkVHCPGzeOp59+mszMTFq0aLF/dcDVq1dzww03kJ+fT/Xq1Xnsscfo2LEjq1ev5rLLLmPnzp0MGzaMBx54gB07djBnzhx+85vfUK9ePZYvX05eXh4TJkzgwQcfpKCggFNOOYVHH32USpUqMXPmTO6880727NlD27ZtefLJJ/cvD3ukevbsuf9+586d2bVrF3v27KFqVU1plojKqAMtT/FbMef8xZo3LoP85SW3uVNhd6k5FHVaHHgitNGJvgsmhSYShRPcr97mT2rE0gld4Zx7y3w6WZZ1nTp1KtnZ2QptST5mfqx4zcwDT4Y6B9u+KNV/HrTSV8+GosLiN0O9rKCbpSM0DMK8YbukbKGnTIs7GZZ1XbJkCb/+9a+ZOXPmUR+/SGSZQZ1mfuswqGR/YQF8/QnkL4P8FUErfblfSbFob8nrap7gZ4M2bB+EeQffSo/wWufhBPdhWsYVLSrLuq5bt44LL7yQZ555hrZt2x5VrSJJqXKVkiVxS9u3118LdFNesK30t7lTfd96sSo1S4V5e6jfpmRL8FZ6yrS4o7ys6zfffMPQoUO5995797feRaQMldJ9F0nDdkCpES7Ffeib8nwLvTjY17wLi6Yc+BnVG5SEeIN20KCtv63fFqoe3TmqeEiZ4I7ysq4PP/wwq1at4q677uKuu+4C/LUsMzO1doTIESvdh57V78DnCr6FLWt810vp7VChXquJD/AGbYLbtv62fusKO0GqZV0PQ8u6iogP9U9h86pgW11y++2mA19buzn0vQ5OG3NMX6VlXWNAy7qKCFWql6zHcrBd3xzYQt+82i/gFWcK7sPQsq4icljV6gbX/8yu0K9Nq9BvExGR41ahwR2P/nTx9LMVSR0VFtwZGRls3rxZARMHzjk2b95MRkZG2KWISAWosD7u5s2bs27dOvLz8yvqK1NKRkYGzZs3D7sMEakAFRbc6enp+6eIi4jIsdPJSRGRiFFwi4hEjIJbRCRi4jLl3czygbXH+PaGwKZyX5WcUvXYddypRcddtlbOuUblfVBcgvt4mNm8I5mrn4xS9dh13KlFx3381FUiIhIxCm4RkYhJxOD+W9gFhChVj13HnVp03Mcp4fq4RUTk8BKxxS0iIoeRUMFtZkPMbIWZrTKz28KuJ17M7Akz22hmuaX21TezWWa2MritF2aN8WBmLczsLTNbamZLzOymYH9SH7uZZZjZh2a2MDju/wn2tzazD4Lf9ylmViXsWuPBzCqZ2cdm9nLwOOmP28zWmNliM1tgZvOCfTH7PU+Y4DazSsAjwDlAJ2CkmXUKt6q4eQoYctC+24DZzrn2wOzgcbIpBG52znUC+gI3BP+Nk/3Y9wADnHPdgR7AEDPrC9wH/Nk51w7YAlwVYo3xdBOwrNTjVDnu7znnepQaAhiz3/OECW6gD7DKOfeJc64AmAwMC7mmuHDO/Rv4+qDdw4Cng/tPAxdUaFEVwDn3pXNufnB/O/5/5mYk+bE7b0fwMD3YHDAAeCHYn3THDWBmzYGhwN+Dx0YKHHcZYvZ7nkjB3Qz4vNTjdcG+VNHYOfdlcP8roHGYxcSbmWUBPYEPSIFjD7oLFgAbgVnAauAb51xh8JJk/X1/APgvoCh43IDUOG4HzDSzHDO7JtgXs99zXXMyATnnnJkl7XAfM6sJTAV+4Zzb5hthXrIeu3NuH9DDzOoC/wQ6hlxS3JnZecBG51yOmZ0Vdj0VrJ9zbr2ZZQKzzGx56SeP9/c8kVrc64EWpR43D/alig1m1gQguN0Ycj1xYWbp+NCe6Jx7MdidEscO4Jz7BngLOBWoa2bFjadk/H0/HTjfzNbguz4HAP9L8h83zrn1we1G/D/UfYjh73kiBfdHQPvgjHMV4BJgesg1VaTpwJXB/SuBl0KsJS6C/s3HgWXOuT+Veiqpj93MGgUtbcysGjAQ37//FnBx8LKkO27n3O3OuebOuSz8/89vOucuI8mP28xqmFmt4vvAICCXGP6eJ9QEHDM7F98nVgl4wjk3LuSS4sLMJgFn4VcL2wDcCUwD/gG0xK+s+CPn3MEnMCPNzPoB7wCLKenz/H/4fu6kPXYz64Y/GVUJ31j6h3PuLjNrg2+J1gc+Bi53zu0Jr9L4CbpKbnHOnZfsxx0c3z+Dh5WB55xz48ysATH6PU+o4BYRkfIlUleJiIgcAQW3iEjEKLhFRCJGwS0iEjEKbhGRiFFwS9yZ2X8Hq+ItClZLOyXYX9nMfh+slrYg2P671Pv2BfuWBCvr3Wxm3/mdNbOmZvZCcL9HMKw0VrXXNbPrD/VdImFRcEtcmdmpwHlAtnOuG/B9StakGQs0Bbo653oA/fELMBXbFayu1hk/aeUc/Jj3AzjnvnDOFU/o6AEcVXCXmsV3KHWB/cF90HeJhELBLfHWBNhUPMHCObfJOfeFmVUHfgbc6JzbHTy33Tn3u0N9SDB1+BpgjJVe3AS/YJWZ5QYzbu8CRgQt9RHBLLYngvWwPzazYcF7RpnZdDN7E5htZjXNbLaZzQ/WUS5emfJeoG3wefcXf1fwGRlm9mTw+o/N7HulPvtFM3st+GviD8H+Smb2VFDrYjP7ZQx/zpJCtMiUxNtM4Ldmlge8AUxxzr0NtAM+C5Z3PSLOuU+Cddsz8TNOD36+wMx+C/R2zo0BMLPf46da/zSYdv6hmb0RvCUb6Oac+zpodV8YLHrVEJhrZtPxayZ3Cf4iKF7VsNgN/mtdVzPriF8NrkPwXA/86od7gBVm9lBQdzPnXJfgs+oe6bGLlKYWt8RVsA51L3xrOR+YYmajDn6dmf0kaNV+bmYtDn7+OAwCbjO/pOocIAM/5RhgVqkpxwb83swW4f+BaUb5y272AyYAOOeW46cxFwf3bOfc1uCviaVAK+AToI2ZPWRWSNlGAAABfElEQVRmQ4BtMTg+SUEKbok759w+59wc59ydwBjgImAV0LJ4MR7n3JNBq3Yrfk2P7wjWgNjH0a2qZsBFQV95D+dcS+dc8dVYdpZ63WVAI6BXUMcGfMgfq9Jrb+wDKjvntgDd8f+AXEtwcQGRo6XglrgysxPNrH2pXT2Atc65b/ErBT5sZhnBaysBh7z+oJk1AsYDD7vDL7CzHahV6vHrwI3F/eJm1rOM99XBrx29N+irblXG55X2Dj7wCbpIWgIryios6IJJc85NBe7Ad9WIHDX1cUu81QQeCvpzC/Et7eIrgvw3cDeQa2bbgV34VfS+CJ6vFnRxpAfvfRYovRzsobxFSdfIPcHnPwAsCoYSfoof5XKwicAMM1sMzAOWAzjnNpvZu8EJyVfx10Ut9ijwl+A9hcAo59yeg86dltYMeLLUkMbbyzkWkUPS6oAiIhGjrhIRkYhRcIuIRIyCW0QkYhTcIiIRo+AWEYkYBbeISMQouEVEIkbBLSISMf8HHNTG1QPG5kIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in [1, 2]:\n",
    "    lambda_ = 0\n",
    "    print(\"Degree \" + str(i))\n",
    "    X_ = build_poly(X, i)\n",
    "    lambda_mae = lambda y, X, w: mae(y, X, w, lambda_)\n",
    "    lambda_grad_mae = lambda y, X, w: mae_stoch_grad(y, X, w, lambda_)\n",
    "    ws, losses = stochastic_gradient_descent(y, X_, np.zeros(X_.shape[1]), 50, 1, lambda_mae, lambda_grad_mae, batch_size = 100, detail = True)\n",
    "    plt.xlabel(\"SGD iterations\")\n",
    "    plt.ylabel(\"RMSE\")\n",
    "    plt.semilogy(losses, label = \"degree =\" + str(i))\n",
    "    plt.legend()\n",
    "    print(\"    loss = \" + str(losses[np.argmin(losses)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loss obtained reaches RMSE of approx. 23 for degree 2 polynomial expansion. We think that we probably won't be able to do better than ridge regression, but the effort is worth trying. Degree 3 has been tested separately and gives RMSE > 1000, hence we exclude it from the possible degree candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda = 0\n",
      "    degree = 1\n",
      "        26.355582665289\n",
      "        26.355784279496454\n",
      "    degree = 2\n",
      "        24.162021525775195\n",
      "        24.162715017406338\n",
      "lambda = 1e-06\n",
      "    degree = 1\n",
      "        26.35558775976868\n",
      "        26.355774697392782\n",
      "    degree = 2\n",
      "        24.170212141069207\n",
      "        24.170990366946498\n",
      "lambda = 1e-05\n",
      "    degree = 1\n",
      "        26.355588727236167\n",
      "        26.355778198307007\n",
      "    degree = 2\n",
      "        24.177040652403203\n",
      "        24.178066413218787\n",
      "lambda = 0.0001\n",
      "    degree = 1\n",
      "        26.355585558178895\n",
      "        26.35578251376678\n",
      "    degree = 2\n",
      "        24.157708463881253\n",
      "        24.15916725428593\n",
      "lambda = 0.001\n",
      "    degree = 1\n",
      "        26.35558779628549\n",
      "        26.355787983716368\n",
      "    degree = 2\n",
      "        24.175575073034818\n",
      "        24.176407819894198\n",
      "lambda = 0.01\n",
      "    degree = 1\n",
      "        26.355634820756176\n",
      "        26.355788383863935\n",
      "    degree = 2\n",
      "        24.231269643413476\n",
      "        24.231798904947226\n"
     ]
    }
   ],
   "source": [
    "run_mae_sgd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we see that the results are not as good as ridge regression..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "We managed to get a lower error than the lab using `Ridge Regression` with polynomial expansion, which we are really proud of. Moreover, we have tried other standard losses and regularizers that did not work as good and we conclude that we probably won't be able to do a lot better using classical machine learning approaches. In parallel, we run a few deep learning experiments that lead to promising results and therefore we decide to stop here our research on the machine learning aspects of the project."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}