Coordinate Descent for Lasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coordinate Descent for Lasso\n",
    "\n",
    "I'm trying to see how good a bad coordinate descent algorithm is compared to the Cython optimised Lasso in SKLearn. If this is promising I will use this trick to do implement the most opaque part of Tibs pliable lasso."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import scipy.optimize as opt\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as graph\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.datasets import load_boston"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(506, 13) (506,)\n"
     ]
    }
   ],
   "source": [
    "x, y = load_boston(True)\n",
    "print(x.shape, y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ground Truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.25 ms ± 12.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "gt = Lasso(alpha=2, fit_intercept=False).fit(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.42991573e-02  5.31785193e-02 -0.00000000e+00  0.00000000e+00\n",
      "  0.00000000e+00  3.05467635e+00  5.11427657e-02 -0.00000000e+00\n",
      "  0.00000000e+00 -1.67773867e-03 -0.00000000e+00  2.07160117e-02\n",
      " -5.98679048e-01]\n"
     ]
    }
   ],
   "source": [
    "gt = Lasso(alpha=2, fit_intercept=False).fit(x, y)\n",
    "print(gt.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Custom Coordinate Descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "506 13\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "Ground Truth Cost 23.315123868478164\n",
      "Initial Cost 296.0734584980237\n",
      "J_1 = 94.26234922473651\n",
      "J_2 = 57.7726619655116\n",
      "J_3 = 48.41119589259517\n",
      "J_4 = 42.95568294818235\n",
      "J_5 = 38.9168862005453\n",
      "J_6 = 35.85729804244692\n",
      "J_7 = 33.61496253552244\n",
      "J_8 = 31.77520110669446\n",
      "J_9 = 30.23466944108655\n",
      "J_10 = 28.959117209399835\n",
      "J_11 = 27.893341691330264\n",
      "J_12 = 27.010537680995554\n",
      "J_13 = 26.398971254403968\n",
      "J_14 = 25.980907825440035\n",
      "J_15 = 25.628038769776367\n",
      "J_16 = 25.335610545813005\n",
      "J_17 = 25.08644816452142\n",
      "J_18 = 24.873629344605245\n",
      "J_19 = 24.690473544249592\n",
      "J_20 = 24.542195925401565\n",
      "J_21 = 24.46161058673571\n",
      "J_22 = 24.39698111652681\n",
      "J_23 = 24.34742409932715\n",
      "J_24 = 24.310168619690295\n",
      "J_25 = 24.27858310662916\n",
      "J_26 = 24.25170100275342\n",
      "J_27 = 24.227546201643726\n",
      "J_28 = 24.20517700205864\n",
      "J_29 = 24.18346778536835\n",
      "J_30 = 24.16187395389199\n",
      "J_31 = 24.140682430093335\n",
      "J_32 = 24.11652074091671\n",
      "J_33 = 24.0975789289459\n",
      "J_34 = 24.078499270659826\n",
      "J_35 = 24.059752084104815\n",
      "J_36 = 24.0407132046279\n",
      "J_37 = 24.02151159377067\n",
      "J_38 = 24.0018022201985\n",
      "J_39 = 23.980855318389807\n",
      "J_40 = 23.956103357614772\n",
      "J_41 = 23.931152725574215\n",
      "J_42 = 23.90653408419343\n",
      "J_43 = 23.881766825433388\n",
      "J_44 = 23.857168671894343\n",
      "J_45 = 23.83193695919731\n",
      "J_46 = 23.8068996826243\n",
      "J_47 = 23.78334303475591\n",
      "J_48 = 23.759168976895506\n",
      "J_49 = 23.734981479684272\n",
      "J_50 = 23.711829399544957\n",
      "J_51 = 23.689052338332512\n",
      "J_52 = 23.666437549670217\n",
      "J_53 = 23.644129599963257\n",
      "J_54 = 23.622870470885722\n",
      "J_55 = 23.60220646953369\n",
      "J_56 = 23.582588786797146\n",
      "J_57 = 23.56400053051466\n",
      "J_58 = 23.546367847435278\n",
      "J_59 = 23.52985085411929\n",
      "J_60 = 23.51440036090157\n",
      "J_61 = 23.49995041250834\n",
      "J_62 = 23.486444192572733\n",
      "J_63 = 23.47372314115652\n",
      "J_64 = 23.461742984991893\n",
      "J_65 = 23.450634261460515\n",
      "J_66 = 23.44012762722003\n",
      "J_67 = 23.430351597499435\n",
      "J_68 = 23.421157753319775\n",
      "J_69 = 23.412537282032314\n",
      "J_70 = 23.405007290646942\n",
      "J_71 = 23.397497463487824\n",
      "J_72 = 23.390361843558694\n",
      "J_73 = 23.38374351794627\n",
      "J_74 = 23.377590225188158\n",
      "J_75 = 23.371906037667433\n",
      "J_76 = 23.366633773010857\n",
      "J_77 = 23.361802970903675\n",
      "J_78 = 23.357297943269607\n",
      "J_79 = 23.353149946535854\n",
      "J_80 = 23.349410413425986\n",
      "J_81 = 23.345959387260336\n",
      "J_82 = 23.34270162352283\n",
      "J_83 = 23.339743556539844\n",
      "J_84 = 23.33700248859556\n",
      "J_85 = 23.334515584910456\n",
      "J_86 = 23.33229154344406\n",
      "J_87 = 23.330444431140158\n",
      "J_88 = 23.328498584828104\n",
      "J_89 = 23.326908665717056\n",
      "J_90 = 23.325443651143722\n",
      "J_91 = 23.32456848361574\n",
      "J_92 = 23.32359779948008\n",
      "J_93 = 23.322720370754475\n",
      "J_94 = 23.32212614536831\n",
      "J_95 = 23.32141338764808\n",
      "J_96 = 23.320756698045024\n",
      "J_97 = 23.320161112839227\n",
      "J_98 = 23.31961526875415\n",
      "J_99 = 23.31912751698447\n",
      "J_100 = 23.31867659065175\n",
      "[-2.65336120e-02  5.33883129e-02  0.00000000e+00  0.00000000e+00\n",
      "  0.00000000e+00  3.08207168e+00  5.35141490e-02  0.00000000e+00\n",
      "  6.51399000e-04 -1.69731226e-03  0.00000000e+00  2.01544610e-02\n",
      " -6.07579156e-01]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAGRVJREFUeJzt3X2wXHWd5/H35z7mERLIJYQkSwJGJTrLw0YWxXFcQQTWAq3SXVh3jG7cWLNMrc5O1SyMZTlTNU7prugsuzs4ODATXQZlFYaUg6ssQqnjErwR5CEhcHlMSCAXQp7JzX347h/965u+ne7TN/ch3efk86rq6nN+/Tunv32SfPqXX58+rYjAzMyKq63ZBZiZ2fRy0JuZFZyD3sys4Bz0ZmYF56A3Mys4B72ZWcE56M3MCs5Bb2ZWcA56M7OC62h2AQALFiyIZcuWNbsMM7Nc2bhx42sR0dOoX0sE/bJly+jt7W12GWZmuSLpxfH089SNmVnBOejNzAquYdBLmiHpYUm/kfSkpD9N7cslbZD0jKTvSepK7d1pvS89vmx6X4KZmWUZz4h+APhARJwLnAdcLuki4KvANyJiBfAGsCb1XwO8ERFvAb6R+pmZWZM0DPoo2Z9WO9MtgA8A30/t64CPpOWr0zrp8UskacoqNjOzYzKuOXpJ7ZIeBXYC9wHPArsjYih12QYsTsuLga0A6fE9wKk19rlWUq+k3v7+/sm9CjMzq2tcQR8RwxFxHrAEuBA4p1a3dF9r9H7Uz1hFxC0RsSoiVvX0NDwN1MzMJuiYzrqJiN3Ag8BFwDxJ5fPwlwDb0/I2YClAevxkYNdUFFttyyv7uPEnW3h9/8B07N7MrBDGc9ZNj6R5aXkmcCmwGXgA+Fjqthq4Jy2vT+ukx38a0/TDtM/27+e//7SP1/Yfno7dm5kVwni+GbsIWCepndIbw50R8UNJm4DvSvoz4BHg1tT/VuA7kvoojeSvmYa6AehsL71PDQ6PTNdTmJnlXsOgj4jHgPNrtD9Hab6+uv0Q8PEpqa6BzvbSxwGHHfRmZnXl+puxXeUR/ZCD3sysnlwHfWdHeepmWj4CMDMrhHwHvefozcwaynnQe47ezKyRXAd9l0f0ZmYN5Troy1M3h/1hrJlZXfkO+g6P6M3MGsl10He0leboh0Z81o2ZWT25Dvry1dOm5wILZmbFkOugLye9c97MrL5cB31b+fdMPKQ3M6sr10FfnrrxFL2ZWX35Dvo0op+mqyCbmRVCvoM+3Tvmzczqy3fQe4rezKyhfAd9GtM7583M6st10I+eXukhvZlZXbkO+vLUjZmZ1ZfvoE/3HtCbmdWX76Avn17pWXozs7ryHfTp3iN6M7P68h30vtaNmVlD+Q768umVTnozs7ryHfSjI3onvZlZPbkO+jKP6M3M6st10Ps8ejOzxvId9PjqlWZmjTQMeklLJT0gabOkJyV9LrX/iaSXJT2abldWbHODpD5JWyR9aLqK90XNzMwa6xhHnyHgDyPi15LmAhsl3Zce+0ZEfK2ys6SVwDXAO4AzgP8r6a0RMTyVhYMvU2xmNh4NR/QRsSMifp2W9wGbgcUZm1wNfDciBiLieaAPuHAqiq125IdHpmPvZmbFcExz9JKWAecDG1LT70t6TNJtkuantsXA1orNtpH9xjBhR0b0Tnozs3rGHfSS5gA/AD4fEXuBm4GzgfOAHcCN5a41Nj8qiSWtldQrqbe/v/+YCy/tI+3cOW9mVte4gl5SJ6WQvz0i7gKIiFcjYjgiRoBvcWR6ZhuwtGLzJcD26n1GxC0RsSoiVvX09EyoeP9mrJlZY+M560bArcDmiPh6Rfuiim4fBZ5Iy+uBayR1S1oOrAAenrqSq+vzh7FmZlnGc9bNxcDvAo9LejS1/TFwraTzKOXsC8BnASLiSUl3ApsonbFz3XSccVMmPHVjZpalYdBHxC+oPe9+b8Y2Xwa+PIm6xk2SP4w1M8uQ62/Ggkf0ZmaN5D/oPUdvZpYp/0GPPKI3M8uQ+6BH/sKUmVmW3Ae9wHM3ZmYZ8h/0nqM3M8uU/6BH/masmVmG/Ae9fHqlmVmW/Ac9nroxM8uS/6CXT680M8uS/6DHp1eamWXJfdDjOXozs0y5D/paV1szM7Mj8h/08umVZmZZChD0PuvGzCxL/oMez9GbmWXJf9D7h0fMzDLlPujbfNaNmVmm3Ac9iBEHvZlZXbkPevk6xWZmmfIf9HjqxswsS/6D3nP0ZmaZ8h/0+KwbM7Ms+Q96j+jNzDLlP+jxR7FmZlnyH/S+Hr2ZWabcBz34evRmZlkaBr2kpZIekLRZ0pOSPpfaT5F0n6Rn0v381C5JN0nqk/SYpAum8wXIczdmZpnGM6IfAv4wIs4BLgKuk7QSuB64PyJWAPendYArgBXptha4ecqrruCrV5qZZWsY9BGxIyJ+nZb3AZuBxcDVwLrUbR3wkbR8NfDtKHkImCdp0ZRXnghfj97MLMsxzdFLWgacD2wAFkbEDii9GQCnpW6Lga0Vm21LbdPCI3ozs2zjDnpJc4AfAJ+PiL1ZXWu0HZXFktZK6pXU29/fP94yaj6ZB/RmZvWNK+gldVIK+dsj4q7U/Gp5Sibd70zt24ClFZsvAbZX7zMibomIVRGxqqenZ6L1p+vRm5lZPeM560bArcDmiPh6xUPrgdVpeTVwT0X7J9PZNxcBe8pTPNOhNKJ31JuZ1dMxjj4XA78LPC7p0dT2x8BXgDslrQFeAj6eHrsXuBLoAw4Cn57Siqt5jt7MLFPDoI+IX1B73h3gkhr9A7huknWNmy9Hb2aWLfffjPVvxpqZZct/0OOzbszMsuQ/6H2ZYjOzTLkP+jZP3ZiZZcp90AOMOOfNzOrKfdD7evRmZtnyH/SAz680M6sv/0HvD2PNzDIVI+ibXYSZWQvLf9D7evRmZpnyH/Qe0ZuZZcp/0OM5ejOzLLkPenw9ejOzTLkPel+P3swsW/6Dvt4FlM3MDChC0OM5ejOzLPkPel/UzMwsU/6DHo/ozcyy5D/ofQkEM7NM+Q96PHVjZpYl90GPR/RmZplyH/TCl0AwM8uS/6B30puZZcp/0HuO3swsU+6Dvq3Nc/RmZllyH/RCjDjpzczqyn/Q+3r0ZmaZGga9pNsk7ZT0REXbn0h6WdKj6XZlxWM3SOqTtEXSh6ar8Eoe0JuZ1TeeEf3fApfXaP9GRJyXbvcCSFoJXAO8I23zl5Lap6rYWuTr0ZuZZWoY9BHxM2DXOPd3NfDdiBiIiOeBPuDCSdTXkMBDejOzDJOZo/99SY+lqZ35qW0xsLWiz7bUNm08R29mlm2iQX8zcDZwHrADuDG11/oZkJo5LGmtpF5Jvf39/RMsw1evNDNrZEJBHxGvRsRwRIwA3+LI9Mw2YGlF1yXA9jr7uCUiVkXEqp6enomUAfh69GZmjUwo6CUtqlj9KFA+I2c9cI2kbknLgRXAw5MrsUEteERvZpalo1EHSXcA7wcWSNoGfAl4v6TzKE3LvAB8FiAinpR0J7AJGAKui4jh6Sm9XJ+D3swsS8Ogj4hrazTfmtH/y8CXJ1PUsfHplWZmWYrxzVgP6c3M6sp/0De7ADOzFpf/oPccvZlZpvwHva9Hb2aWKf9B7xG9mVmmYgR9s4swM2th+Q965LNuzMwy5D7o8YjezCxT7oO+dJniZldhZta68h/0/uERM7NMuQ/6Nn8z1swsU+6DXsCIc97MrK78B72vR29mlin3QQ+wddebbNq+t9llmJm1pNwH/fveugCALa866M3Masl90F+4/FQADg+NNLkSM7PWlPug7+4ovYQBB72ZWU25D/quFPQe0ZuZ1Zb7oPeI3swsW+6DvqvdQW9mliX3QS+Jro42T92YmdWR+6AH6G5v48DAULPLMDNrSYUI+hUL5/APj+9gcNijejOzaoUI+ktXLmTXgcPsO+RRvZlZtUIE/YLZ3QAcPOygNzOrVoign9nVDsCbh4ebXImZWespRNDPSkF/0EFvZnaUQgR9eUS/+83BJldiZtZ6Gga9pNsk7ZT0REXbKZLuk/RMup+f2iXpJkl9kh6TdMF0Fl+24rS5dLaLnz/dfzyezswsV8Yzov9b4PKqtuuB+yNiBXB/Wge4AliRbmuBm6emzGw9c7s5/5/M51cvvnE8ns7MLFcaBn1E/AzYVdV8NbAuLa8DPlLR/u0oeQiYJ2nRVBWb5dwlJ7Np+x7/fqyZWZWJztEvjIgdAOn+tNS+GNha0W9bajuKpLWSeiX19vdPfspl3qwuBoeDw/7SlJnZGFP9YaxqtNUcYkfELRGxKiJW9fT0TPqJR8+8GfCZN2ZmlSYa9K+Wp2TS/c7Uvg1YWtFvCbB94uWN3+yuDgAODjrozcwqTTTo1wOr0/Jq4J6K9k+ms28uAvaUp3im28zREb2/HWtmVqmjUQdJdwDvBxZI2gZ8CfgKcKekNcBLwMdT93uBK4E+4CDw6Wmouaby1M0Bf2nKzGyMhkEfEdfWeeiSGn0DuG6yRU3EjM5S0A946sbMbIxCfDMW/JOCZmb1FCjo04jeQW9mNkZxgr6zPKL31I2ZWaXiBH156mbQI3ozs0oFCnpP3ZiZ1VKgoPfUjZlZLcUJ+k6fdWNmVkthgr6r3XP0Zma1FCboO9rb6GiTp27MzKoUJuihNE/vqRszs7GKFfSd7R7Rm5lVKVbQd7R5jt7MrErxgt5TN2ZmYxQq6Gd0tvOmr15pZjZGoYL+pJmd7H1zsNllmJm1lEIF/byZnew+6KA3M6tUqKA/aWYnW17dx8YXdzW7FDOzllGooL/q3DMA2LR9b5MrMTNrHYUK+lXL5gP+3Vgzs0qFCvoZHe1IcHBgqNmlmJm1jEIFfVubmNXZ7hG9mVmFQgU9wOzuDu55dDs79x1qdilmZi2hcEG/9n1n8dr+Af79ul4iotnlmJk1XUezC5hqn/ntszg0OMzXfvI0uw4c5tQ53c0uycysqQo3ogd468K5ADzbf6DJlZiZNV8hg/6fnTmf9jbx4JadzS7FzKzpChn0p87p5sxTZ/GXDz7LQ8+93uxyzMyaalJBL+kFSY9LelRSb2o7RdJ9kp5J9/OnptRj8/lL3wrAxhffaMbTm5m1jKkY0f+LiDgvIlal9euB+yNiBXB/Wj/urjr3DLo62nw1SzM74U3H1M3VwLq0vA74yDQ8x7jMm9nJ5lf2sfvg4WaVYGbWdJMN+gB+ImmjpLWpbWFE7ABI96dN8jkm7G2nz+VnT/dzxX/7OW/627JmdoKabNBfHBEXAFcA10l633g3lLRWUq+k3v7+/kmWUdvffOpd/PlHf4sdew7xZ/+wyV+gMrMT0qSCPiK2p/udwN3AhcCrkhYBpPua5zhGxC0RsSoiVvX09EymjLo62tv41+9aym8tPpnbN7zEHQ9vnZbnMTNrZRMOekmzJc0tLwOXAU8A64HVqdtq4J7JFjkZ7W3iO2suBKBv5/5mlmJm1hSTuQTCQuBuSeX9/F1E/B9JvwLulLQGeAn4+OTLnJx5s7pYvmA2t/3j8yxfMIt/e9GZpLrNzApvwkEfEc8B59Zofx24ZDJFTYfV7z6T7zz0Il+850keeWk3N/6rcx32ZnZCKOQ3Y2v51MXLue8Pfocr3nk6dz3yMt/6+XPNLsnM7Lg4YYIeSj9M8tWP/VMA/vzep/jHvtd8Jo6ZFd4JFfQAJ83oZM17lwPwib/ewL+86Re8ccBfqDKz4jrhgh7gix9eyW++dBlr3rucTTv28tn/tZHDQyPNLsvMbFqckEEPcPLMTr5w5Tmcu3QeDz+/i7d/8Uf8m289xCt7/BOEZlYsJ2zQQ2nO/vbP/HNuuvZ8rjr3DH757Otc9T9+wd2PbHPgm1lhqBU+jFy1alX09vY2uwy+/f9e4MafPM2edMXL/3jJCn7vd85mZld7cwszM6tB0saKKwfX7+egH+vQ4DAPPfc6X/nRUzz1yj4Azlowmz+6/G1ctvJ02tp87r2ZtQYH/SSNjAQ/73uN+za9wt2/fpkDh4eZP6uTs3vm8OmLl/PBlQvp6jihZ77MrMkc9FPo0OAw63+znQe37OTBLf0cPDxMZ7tYPG8mH1y5kEvPWcjZp81hwZzuZpdqZicQB/00GRwe4YGndvLLZ1+n98VdbN6xj+GR0jF8++lzefvpczln0Umcs+gk3nLaHHrmdtPZ7pG/mU298Qb9ZC5qdkLqbG/jsneczmXvOB2AvYcG+WXfazzy0m6efnUfG57fxd8/un3MNqfM7uK0ud30zO3m1NldzJvVxfxZXcyf3cm8WV2cNKODOd0dzCnfd3cwu7vDbxBmNiUc9JN00oxOLn/nIi5/56LRtjcOHGbzjr288PpBdu47xM59A+zcO8Br+wd44fUD7D4wyL6BoYb77u5oY0ZnOzM72zlpZgczOtvp7miju6Odro42OttFe5to05H70nLp8sySaB/zWKm9ra203CYhQBISiFK7NLZNotSeliFtK8ZuP7q/qn1R6qjK7Sr2V3e7ivrqtqX+jNm2omYdee1tla+zbexrY3T5SC1t6cWqettaz5deG6Ov7ejXVT52Y157dT9faM+mgYN+Gsyf3cV73rKA97ylfp/B4RH2vDnI7oOH2fPmIAcGhjkwMMS+gSH2Hxpi/0DpNjA4zMHDw+w9NMjA0AgDgyMcPDzE7jdHGBwKhiMYGQlGorwMw2l9JILhkdItAobTeukxiHRvrafemwWj7bXfLEobM7pcftNLzVXPcfSbytF9albXsM949qMJ7afxG2HN/VS1VT937T619qOGfaobG+3nmnct5TO/fVatPU0ZB32TdLa3sWBOd0t8gBtReiMYiSCACAhKbZXLlY8zZr1iu6p9lN9Qyh8FxXi3Gzm6hrFvTuWajrQFkeo6su1wxGitpTe5qm0rXtuY+kqbldrHPP+RZSreMKNiW6qPYdV65TE/6ninjWu1l9cZXR/7WHplo8vl5ym1V/+Z1/h7UNWrdp/G+6nuVXM/VW3Vz127z8T2U91Uez8xjj5Ts5/qxuORAQ56OzLVUXt8YmY550/7zMwKzkFvZlZwDnozs4Jz0JuZFZyD3sys4Bz0ZmYF56A3Mys4B72ZWcG1xNUrJfUDL05w8wXAa1NYzvHk2pvDtTdHXmtv5brPjIieRp1aIugnQ1LveC7T2Ypce3O49ubIa+15rbuSp27MzArOQW9mVnBFCPpbml3AJLj25nDtzZHX2vNa96jcz9GbmVm2IozozcwsQ66DXtLlkrZI6pN0fbPrqSZpqaQHJG2W9KSkz6X2UyTdJ+mZdD8/tUvSTen1PCbpgibX3y7pEUk/TOvLJW1IdX9PUldq707rfenxZU2ue56k70t6Kh37d+fomP9B+rvyhKQ7JM1o1eMu6TZJOyU9UdF2zMdZ0urU/xlJq5tY+39Nf2cek3S3pHkVj92Qat8i6UMV7S2dQaNKv1KTvxvQDjwLnAV0Ab8BVja7rqoaFwEXpOW5wNPASuC/ANen9uuBr6blK4EfUfr1sYuADU2u/z8Bfwf8MK3fCVyTlr8J/F5a/g/AN9PyNcD3mlz3OuAzabkLmJeHYw4sBp4HZlYc70+16nEH3gdcADxR0XZMxxk4BXgu3c9Py/ObVPtlQEda/mpF7StTvnQDy1PutOchg0ZfW7MLmMQf1LuBH1es3wDc0Oy6GtR8D/BBYAuwKLUtArak5b8Crq3oP9qvCbUuAe4HPgD8MP0Dfa3iH8Lo8Qd+DLw7LXekfmpS3SelsFRVex6O+WJgawq9jnTcP9TKxx1YVhWWx3ScgWuBv6poH9PveNZe9dhHgdvT8phsKR/3PGVQnqduyv8oyraltpaU/lt9PrABWBgROwDS/WmpWyu9pr8A/ggYSeunArsjYiitV9Y2Wnd6fE/q3wxnAf3A36Rpp7+WNJscHPOIeBn4GvASsIPScdxIPo572bEe55Y5/lX+HaX/gUD+aj9KnoO+1g+ctuQpRJLmAD8APh8Re7O61mg77q9J0oeBnRGxsbK5RtcYx2PHWwel/5LfHBHnAwcoTSHU0zK1p/nsqylND5wBzAauqNG1FY97I/VqbbnXIOkLwBBwe7mpRreWrL2ePAf9NmBpxfoSYHuTaqlLUielkL89Iu5Kza9KWpQeXwTsTO2t8pouBq6S9ALwXUrTN38BzJNU/kH5ytpG606PnwzsOp4FV9gGbIuIDWn9+5SCv9WPOcClwPMR0R8Rg8BdwHvIx3EvO9bj3ErHn/Rh8IeBT0SajyEntWfJc9D/CliRzkjoovRh1Pom1zSGJAG3Apsj4usVD60HymcXrKY0d19u/2Q6Q+EiYE/5v8HHU0TcEBFLImIZpeP604j4BPAA8LE6dZdfz8dS/6aMbCLiFWCrpLelpkuATbT4MU9eAi6SNCv93SnX3vLHvcKxHucfA5dJmp/+R3NZajvuJF0O/Gfgqog4WPHQeuCadJbTcmAF8DA5yKBRzf6QYDI3Sp/kP03pk+8vNLueGvW9l9J/5R4DHk23KynNo94PPJPuT0n9BfzP9HoeB1a1wGt4P0fOujmL0l/wPuB/A92pfUZa70uPn9Xkms8DetNx/3tKZ3Pk4pgDfwo8BTwBfIfSmR4tedyBOyh9ljBIaXS7ZiLHmdJ8eF+6fbqJtfdRmnMv/1v9ZkX/L6TatwBXVLS3dAaVb/5mrJlZweV56sbMzMbBQW9mVnAOejOzgnPQm5kVnIPezKzgHPRmZgXnoDczKzgHvZlZwf1/xY2r1gnt/10AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Helper Functions\n",
    "def model(x, m):\n",
    "    return x @ m\n",
    "\n",
    "def objective_function(x, y, m, lam):\n",
    "    return 0.5 * ((y - model(x, m))**2).mean() + lam * la.norm(m, 1)\n",
    "\n",
    "def partial_objective_function(mi, x, y, m, lam, pi):\n",
    "    m[pi] = mi\n",
    "    return objective_function(x, y, m , lam)\n",
    "\n",
    "# Class\n",
    "class CustomLasso:\n",
    "    def __init__(self, lam):\n",
    "        self.lam = lam\n",
    "        self.coef_ = None\n",
    "        self.history = []\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        self.history = []\n",
    "        \n",
    "        n, p = X.shape\n",
    "        print(n, p)\n",
    "        \n",
    "        m = np.zeros(p)\n",
    "        print(m)\n",
    "        \n",
    "        # Coordinate Descent\n",
    "        self.history.append(objective_function(x, y, m, self.lam))\n",
    "        print(f'Ground Truth Cost {objective_function(x, y, gt.coef_, gt.alpha)}')\n",
    "        print(f'Initial Cost {self.history[-1]}')\n",
    "        \n",
    "        for i in range(100):\n",
    "            for pi in range(p):\n",
    "                m[pi] = opt.fmin(partial_objective_function, m[pi], args=(x, y, m.copy(), self.lam, pi), disp=False)\n",
    "                m[abs(m) < 1e-15] = 0  # Polish\n",
    "                self.history.append(objective_function(x, y, m, self.lam))\n",
    "            print(f'J_{i+1} = {self.history[-1]}')\n",
    "        \n",
    "        print(m)\n",
    "        self.coef_ = m.copy()\n",
    "        return self\n",
    "    \n",
    "    def predict(self, X):\n",
    "        return model(X, self.coef_)\n",
    "    \n",
    "\n",
    "# Test function\n",
    "lasso = CustomLasso(lam=2).fit(x, y)\n",
    "\n",
    "graph.plot(lasso.history)\n",
    "graph.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd0AAAHVCAYAAAC0biEDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xl4VPWh//HPNyEJkS2igLJogpadXlZb5EKLIqhUrdalWOvSX7VXa1VkM4AEQcEQRVrLLeWCVdS6c1GxSNGyqFCRHcIqSyQBQxASEgJkMvn+/pDL44KSbeZ75pz363l4nmSYZ/i0I3zyOefMjLHWCgAARF6c6wAAAAQFpQsAQJRQugAARAmlCwBAlFC6AABECaULAECUULoAAEQJpQsAQJRQugAAREmdSDzo2WefbVNTUyPx0AAAeM6qVasOWGubnO5+ESnd1NRUrVy5MhIPDQCA5xhjcipzPw4vAwAQJZQuAABRQukCABAllC4AAFFC6QIAECWULgAAUULpAgAQJZQuAABRQukCABAllC4AAFFC6QIAECWULgAAUULpAgAQJZQuAABRQukCABAllC4AAFFC6QIAECWULgAgUJYtW+bsz6Z0AQCB8cYbb+j48ePO/vw6zv5kAACixFqrqVOn6t5771VCQoKzHCxdAIDvLV68WBdffLHTwpUoXQCAj1VUVGjs2LHq1q2bfvSjH7mOw+FlAIA/hcNhbdy4Uddff70aNWrkOo4kli4AwIfKy8v10EMP6ayzztIPf/hD13FOYukCAHwlFArp008/1Z133qmWLVu6jvM1LF0AgG8cP35cw4cPV926ddWmTRvXcb6FpQsA8IVjx45p69atGjZsmOcW7v9h6QIAYl4oFNKwYcN09tlne7ZwJZYuACDGlZSUaNWqVZo0aZIaNGjgOs73YukCAGKWtVYZGRnq0KGD5wtXYukCAGLUoUOHtHDhQmVlZSkuLjY2ZGykBADgG/76179qwIABMVO4EksXABBj8vPz9dprr+mhhx5yHaXKYufHAwBA4FlrNX/+fN1xxx2uo1QLSxcAEBP27Nmj//mf/9H48eNdR6k2li4AwPNKS0uVnZ2t0aNHu45SI5QuAMDTPv30Uz388MMaOHCgkpKSXMepEUoXAOBZe/bsUXFxsTIzM2WMcR2nxihdAIAnZWdn689//rN++MMfqk4df1yCROkCADxn/fr1SkpK0qRJkxQfH+86Tq2hdAEAnrJp0ya9/vrrat26dUy98UVl+Ot/DQAgpi1btkxJSUl65JFHfFe4EqULAPCIPXv2aMmSJWrdurUvLpo6FUoXAODce++9py+++ELp6em+LVyJ0gUAOFZUVKQtW7aoS5curqNEnD+uwQYAxKS33npL9evX17333us6SlSwdAEAThw7dkwlJSW65JJLXEeJGpYuACDqXnrpJdWtW1c333yz6yhRRekCAKLq4MGDatasWaAW7v+hdAEAUfPMM8+oYcOGuv76611HcYLSBQBExcaNG9W7d2+1bdvWdRRnuJAKABBx06dP1+7duwNduBJLFwAQYe+++64GDx6sRo0auY7iHEsXABAxs2bNUmJiIoV7AksXAFDrrLWaNWuW7rjjDl99NF9NUboAgGqbuyZPWQu2am/hUTVPSdbwgW31864tNGfOHHXs2JHC/QZjra31B+3Ro4dduXJlrT8uAMA75q7JU/qcDToaCp+8rW4do55FH2jWk48oMTHRYbroMsasstb2ON39OKcLAKiWrAVbv1a41loV7lyvDXGtA1W4VUHpAgCqZW/h0ZNf24qwCpc8q8RmrVWU3NxhKm+jdAEA1dI8JVmSZMPlCh3IUf3Olykuqd7J2/FtlC4AoFqGD2yrJBPWoUXPKC65oRLOaqnkhHgNHxjsN8D4Ply9DAColsvbn6XPutfXcxU36aBp9LWrl3FqlC4AoMpKS0s1cuRIjRw5UvddH7xPC6qu05auMaaupKWSkk7c/3VrbUakgwEAvKmkpETbtm3Tww8/rKZNm7qOE1Mqc073uKRLrLX/IamLpMuNMT+ObCwAgBeFQiGNHDlSLVu2pHCr4bRL13757hklJ75NOPGr9t9RAwDgaYcOHdKqVav05JNPqm7duq7jxKRKXb1sjIk3xqyVtF/SQmvtx6e4z13GmJXGmJUFBQW1nRMA4JC1Vo8//rh69OhB4dZApUrXWhu21naR1FLSRcaYTqe4zwxrbQ9rbY8mTZrUdk4AgCP79+/X3//+dz3++ONKSUlxHSemVel1utbaQkmLJV0ekTQAAM954YUXdM0118gY4zpKzKvM1ctNJIWstYXGmGRJ/SVlRjwZAMCp3Nxcvfbaa3rwwQddR/GNyrxO91xJzxlj4vXlMn7VWjsvsrEAAC6Fw2EtXbpUd999t+sovlKZq5fXS+oahSwAAA/YsWOHnn32WU2YMMF1FN/hvZcBACcVFhbqs88+09ixY11H8SVKFwAgSdq4caMmTZqkn/70p0pISHAdx5coXQCAtm/fLmOMJk6cyFXKEUTpAkDArVmzRrNnz1b79u0VHx/vOo6vUboAEGArVqxQ48aNNX78eMXFUQmRxv/DABBQGzdu1HvvvafzzjuPQ8pRQukCQAC9//77qlevntLT0yncKKJ0ASBg9uzZo/Xr1ystLY3CjTJKFwAC5O2339aBAwc0ZMgQ11ECidIFgIA4fPiw8vPz1bUrbzLoSmXeexkAEONefvllNW3aVL/97W9dRwk0li4A+FxJSYmSkpJ0ySWXuI4SeCxdAPCxv/3tbzrrrLN07bXXuo4CUboA4Ft79+5Vu3bt1KtXL9dRcAKlCwA+NG3aNJ133nm66qqrXEfBV1C6AOAzn3zyia6++mq1atXKdRR8AxdSAYCP/OlPf9Lhw4cpXI9i6QKAD1hrNXfuXN15551KTk52HQffgaULAD4wa9YsNWvWjML1OJYuAMSwiooK/eUvf9E999zD+yjHAJYuAMSwN998U71796ZwYwSlCwAxqLy8XBkZGRo0aJC6dOniOg4qidIFgBhTUVGh5cuXa/DgwUpMTHQdB1VA6QJADCkrK9PQoUPVuXNntWvXznUcVBEXUgFAjDh27Ji2bt2q+++/XykpKa7joBpYugAQA0pLSzVixAg1a9ZMqampruOgmli6AOBxR44c0Y4dO5Senq5zzjnHdRzUAEsXADzsyJEjGjFihM455xyde+65ruOghli6AOBRhYWF2r59ux577DHO4foESxcAPCgUCmnMmDG68MILKVwfYekCgMfs379fq1at0tSpU1WnDv9M+wlLFwA8xFqrP/3pT+rTpw+F60M8owDgEbm5uXr//ff16KOPuo6CCGHpAoBHzJ07VzfddJPrGIggli4AOLZjxw7NnTtXQ4cOdR0FEcbSBQCHysrKtHbtWt13332uoyAKKF0AcCQ7O1sTJkzQL37xCyUkJLiOgyigdAHAgfz8fBUVFWncuHGuoyCKKF0AiLJVq1bp6aefVq9evRQfH+86DqKI0gWAKNqwYYMaNWqkCRMmyBjjOg6ijNIFgChZsWKF3nzzTV1wwQUUbkBRugAQBUuWLFHLli01evRoCjfAKF0AiLD169dr1apVat68OYUbcJQuAETQvHnz1LBhQz344IOuo8ADKF0AiJDc3Fzl5uYqNTXVdRR4BKULABHw6quvqqCgQP/1X//lOgo8hNIFgFpWVFSkUCikrl27uo4Cj+EDDwCgFj377LM677zz9Ktf/cp1FHgQSxcAaklhYaHOOeccXXLJJa6jwKNYugBQC/77v/9baWlpuuKKK1xHgYdRugBQQzt27NDFF1+sLl26uI4Cj+PwMgDUwJNPPqldu3ZRuKgUli4AVIO1Vh988IFuvfVWNWnSxHUcxAiWLgBUw9SpU1WnTh0KF1XC0gWAKqioqNArr7yie++9VwkJCa7jIMawdAGgCp599lm1a9eOwkW1sHQBoBLC4bCmTp2qBx98kE8KQrWxdAGgEubNm6dBgwZRuKgRShcAvkdZWZkeeughXXHFFWrXrp3rOIhxlC4AfIdQKKSVK1fqd7/7nRITE13HgQ9QugBwCkePHtWDDz6oDh06KC0tzXUc+AQXUgHANxw5ckSffvqp0tPTlZKS4joOfISlCwBfUVxcrJEjR+rcc89V8+bNXceBz7B0AeCEwsJC5eTkaNy4cTr77LNdx4EPsXQBQF8u3FGjRqlVq1YULiKGpQsg8AoKCrRjxw5lZWWpXr16ruPAx1i6AAItFArp0UcfVceOHSlcRBxLF0Bg5eXlafXq1Zo6dSrvNIWoYOkCCCRrrWbNmqX+/ftTuIgali6AwNm5c6cWL16ssWPHuo6CgGHpAggUa63ef/99/frXv3YdBQHE0gUQGJs2bdK8efM0YsQI11EQUCxdAIFQWlqqnTt3aujQoa6jIMAoXQC+t3r1ak2ePFk/+9nPFB8f7zoOAozSBeBrOTk5kqSMjAzHSQBKF4CPffTRR5o9e7a6du3Ky4LgCZQuAF9asWKFzjvvPI0ZM4bChWdQugB8Z9myZfrggw/UqlUrCheewkuGAPjKu+++q27duuniiy92HQX4FpYuAN9Yv369du3apaZNm7qOApwSpQvAF1577TU1atRId999t+sowHeidAHEvNzcXJWWlur88893HQX4XpQugJg2e/ZsFRQU6LbbbnMdBTgtShdAzDp06JAaNGigrl27uo4CVApXLwOISX/5y1/Uvn17XXvtta6jAJVG6QKIOfv371fHjh3Vt29f11GAKuHwMoCY8sQTT2jDhg0ULmLSaZeuMaaVpNmSzpFUIWmGtfaPkQ4GAF9lrVV2drauvvpqtWnTxnUcoFoqs3TLJQ211raX9GNJvzfGdIhsLAD4uszMTBUVFVG4iGmnXbrW2n2S9p34utgYs1lSC0mbIpwNAGSt1T//+U/dc889atiwoes4QI1U6ZyuMSZVUldJH5/i9+4yxqw0xqwsKCionXQAAu+pp55SkyZNKFz4QqWvXjbG1Jf0hqQHrLWHv/n71toZkmZIUo8ePWytJQQQSOFwWM8995weeOABxcVxzSf8oVL/JRtjEvRl4b5orZ0T2UgAID3//PPq1asXhQtfqczVy0bSLEmbrbVTIh8JQJCFQiFNnjxZo0aN4rNw4TuV+RGyt6RfS7rEGLP2xK8rI5wLQABZa7VgwQINHjyYwoUvVebq5Q8l8V8/gIg6evSoRo0apczMTCUmJrqOA0QEJ0sAOHfs2DGtW7dOw4YNo3Dha5QuAKeKi4s1dOhQtWvXTi1atHAdB4goPvAAgDNFRUXavXu3xo0bp5SUFNdxgIhj6QJw4tChQ0pPT1erVq3UpEkT13GAqGDpAoi6AwcOKC8vT5MmTVKjRo1cxwGihqULIKoOHz6sjIwMtW7dmsJF4LB0AUTN3r17tWvXLk2ZMkVJSUmu4wBRx9IFEBWhUEhPPvmkunXrRuEisFi6ACJu165dWrdunZ544gneaQqBxtIFEFHWWr3yyiu68sorKVwEHksXQMRs3rxZH374oR566CHXUQBPYOkCiIhwOKwVK1boN7/5jesogGewdAHUujVr1mjhwoUaMWKE6yiAp7B0AdSqwsJCffHFFxo2bJjrKIDnULoAas2yZcs0bdo09e/fX3Fx/PMCfBN/KwDUii1btiglJUWjRo1yHQXwLEoXQI0tWrRI8+bNU/v27XlZEPA9uJAKQI0sWbJEHTt2VL9+/VxHATyPpQug2j744ANlZ2eradOmrqMAMYGlC6Ba/vd//1d9+/ZVnz59XEcBYgZLF0CVrV+/XkVFRTrrrLNcRwFiCqULoEqef/55NWzYULfffrvrKEDMoXQBVFpeXp6SkpKUmprqOgoQkyhdAJUyffp05efn68Ybb3QdBYhZlC6A0zpw4IBSU1PVrVs311GAmEbpAvheTz75pDZv3qzLL7/cdRQg5vGSIQCnZK3Vnj171LdvX/Xs2dN1HMAXWLoAvsVaq4kTJ+qzzz6jcIFaxNIF8DXWWq1cuVK33nqrWrVq5ToO4CssXQBf89hjjykuLo7CBSKApQtAkhQOhzVv3jwNHTpUycnJruMAvsTSBSBJ+uMf/6gLL7yQwgUiiKULBFwoFNLMmTM1ZMgQPgsXiDCWLhBwL730kgYMGEDhAlHA0gUC6tixY5o0aZLGjRtH4QJRwtIFAqiiokKLFi3SnXfeSeECUUTpAgFTUlKi+++/X5deeqlatmzpOg4QKJQuECAlJSXavHmzxo4dq8TERNdxgMChdIGAOHTokEaMGKE2bdqoSZMmruMAgcSFVEAAHDhwQHl5eZo4caIaNWrkOg4QWCxdwOcOHDigjIwMtW7dWikpKa7jAIHG0gV8bN++fcrPz9fkyZNVr14913GAwGPpAj5VXFysSZMmqW3bthQu4BEsXcCHdu/erc8++0xTpkxRnTr8NQe8gqUL+Ex5eblmzJihH//4xxQu4DH8jQR8ZMuWLcrOztbEiRNdRwFwCixdwCestZo/f76uueYa11EAfAeWLuADa9eu1ccff6whQ4a4jgLge7B0gRhXVlamrVu36s4773QdBcBpsHSBGLZ8+XJ99NFHGjZsmOsoACqBpQvEqPz8fFVUVGjo0KGuowCoJEoXiEHvv/++nn/+efXu3ZvPwwViCIeXgRizevVqtW7dWpdeeqnrKACqiKULxJAFCxZo+fLlSktLcx0FQDWwdIEY8e677+pHP/qRBg4c6DoKgGpi6QIxYOnSpdq3b5/OPPNM11EA1ABLF/C4v//977riiivUt29f11EA1BBLF/CwDRs2qE6dOixcwCcoXcCj/vrXv6pBgwa68cYbXUcBUEsoXcCDcnNz1bx5c6WmprqOAqAWUbqAxzz11FMqKCjQVVdd5ToKgFrGhVSAR1hr9fnnn6tbt27q2rWr6zgAIoClC3iAtVYTJ05UTk6OfvKTn7iOAyBCWLqAY9Zabd++Xddcc406derkOg6ACGLpAg5Za5WRkaHi4mIKFwgAli7gSDgc1rJly/SHP/xBTZo0cR0HQBSwdAFHJkyYoMaNG1O4QICwdIEoC4VCmjNnjkaPHq2EhATXcQBEEUsXiLJp06ape/fuFC4QQCxdIEqOHTum6dOn64EHHnAdBYAjLF0gSt544w1df/31rmMAcIjSBSKspKRE6enpuvnmm9WyZUvXcQA4ROkCERQKhbR8+XINGTJExhjXcQA4RukCEVJYWKgHHnhAffv2VdOmTV3HAeABXEgFRMChQ4e0c+dOPfroo0pKSnIdB4BHsHSBWrZ//36NHj1abdq00Zlnnuk6DgAPYekCtWjfvn0qKCjQ5MmTVb9+fddxAHgMSxeoJfn5+Zo4caIuvPBCChfAKbF0gVqQk5OjgwcPKisrS3Xr1nUdB4BHsXSBGiouLtbUqVPVoUMHChfA92LpAjWwdetW7d27V1OmTOF1uABOi6ULVFN5ebleeukl9enTh8IFUCksXaAa1q1bp23btmncuHGuowCIISxdoIqstVq2bJl+8YtfuI4CIMawdIEq+Pe//61169bp7rvvdh0FQAxi6QKVVFpaqv379+uuu+5yHQVAjGLpApXwr3/9S+vXr+cD6AHUyGmXrjHmGWPMfmPMxmgEArxm9+7dSklJoXAB1FhlDi8/K+nyCOcAPOmdd97RvHnz1K1bN9dRAPjAaQ8vW2uXGmNSIx8F8JaPPvpI3bt316BBg1xHAeATtXYhlTHmLmPMSmPMyoKCgtp6WMCJf/zjH9q6davOOecc11EA+EitXUhlrZ0haYYk9ejRw9bW4wLRNmfOHA0YMIBPCgJQ63jJEPAVS5cuVVlZGYULICIoXeCEmTNnqnPnzvrlL3/pOgoAn6rMS4ZekrRcUltjTK4x5v9FPhYQXRs3blTTpk115plnuo4CwMdOW7rW2sHW2nOttQnW2pbW2lnRCAZEy9SpU1WvXj1dffXVrqMA8DnekQqBZa3Vnj171KlTJ6WlpbmOAyAAOKeLQLLWatKkSTp48KD69+/vOg6AgKB0ETj/t3D79++vLl26uI4DIEAoXQRKRUWFMjIyVFBQoIsuush1HAABwzldBEZFRYWys7N122236YILLnAdB0AAsXQRCNZajRo1SsYYCheAMyxd+F4oFNLixYs1atQoNWzY0HUcAAHG0oXvPfbYY0pLS6NwATjH0oVvHTt2TK+++qrGjh2ruDh+vgTgHv8Swbdmzpypfv36UbgAPIOlC985cuSIpk2bphEjRriOAgBfwwSAr1hrNW/ePN12222uowDAt1C68I3CwkINGzZMN954o5o1a+Y6DgB8C6ULXzh69KjWrFmjMWPGyBjjOg4AnBKli5hXUFCgoUOH6uKLL+bzcAF4GhdSIabt379feXl5yszMVFJSkus4APC9WLqIWXl5eZowYYLatm2rBg0auI4DAKfF0kVMysnJUXFxsbKyslS3bl3XcQCgUli6iDmff/65nnrqKf3gBz+gcAHEFJYuYsq2bdt05MgRZWVlKSEhwXUcAKgSli5iRnFxsWbOnKnOnTtTuABiEksXMWH9+vUqKChQZmYmr8MFELNYuvC88vJyzZ8/X/369aNwAcQ0li487eOPP1ZOTo5GjhzpOgoA1BhLF54VDoe1efNm3XDDDa6jAECtYOnCkxYtWqTt27frrrvuch0FAGoNSxeeU1hYqPLycgoXgO+wdOEp77zzjnbv3q3f//73rqMAQK2jdOEZmzdvVuvWrTVo0CDXUQAgIji8DE+YM2eOli1bpvbt27uOAgARw9KFcwsXLlS/fv34LFwAvkfpwqm3335bJSUlFC6AQKB04cyLL76o6667TsnJya6jAEBUcE4XTixdulT169encAEECksXUff000/rV7/6lRo3buw6CgBEFaWLqLHWKjs7W+3ataNwAQQSh5cRFdZaZWZmql69errssstcxwEAJyhdRJy1Vjk5OerTp4/S0tJcxwEAZyhdRFRFRYUyMjJUUlKi3r17u44DAE5RuoiYiooK7dq1SzfccIM6derkOg4AOEfpIiLC4bDS09NVWlqqzp07u44DAJ7A1cuodeXl5Vq3bp3+8Ic/qGXLlq7jAIBnsHRRq6y1Sk9PV8OGDSlcAPgGli5qzfHjx/Xee+9p/PjxvNMUAJwCSxe1JjMzU507d6ZwAeA7sHRRY0eOHNHLL7+shx9+WMYY13EAwLNYuqix2bNn62c/+xmFCwCnwdJFtRUVFenpp5/WmDFjXEcBgJjA0kW1VFRU6L333tM999zjOgoAxAxKF1V24MABDR06VNdddx2fFgQAVUDpokoOHz6sLVu2aPz48ZzDBYAqonRRaXl5eRo9erR69uypBg0auI4DADGHC6lQKbm5uTp48KAmT56spKQk13EAICaxdHFaOTk5euKJJ9S2bVve+AIAaoCli++1fft2hUIhTZ48WYmJia7jAEBMY+niO+3bt08zZsxQmzZtKFwAqAUsXZzShg0bVFZWpscff1zx8fGu4wCAL7B08S0lJSV6/fXX1bVrVwoXAGoRSxdfs2LFChUVFemRRx5xHQUAfIeli5PKy8u1fPly9e/f33UUAPAlli4kSYsXL1Z+fr7uv/9+11EAwLdYulBZWZny8/N10003uY4CAL7G0g24d955R/v27dNvf/tb11EAwPco3QD7/PPPlZKSokGDBrmOAgCBQOkG1Ouvv67CwkIWLgBEEaUbQKtWrVLPnj11/vnnu44CAIFC6QbMyy+/rIqKCnXv3t11FAAIHEo3QN566y1dddVVqlevnusoABBIlG5AvPnmm7LWUrgA4BClGwAzZ87ULbfcorp167qOAgCBxptj+NzSpUvVqlUrChcAPIDS9SlrrbKystSpUycNHDjQdRwAgChdX7LWKjs7W7169VLjxo1dxwEAnEDp+oy1VuPGjVO9evX0n//5n67jAAC+ggupfKSiokK7du3S1VdfrbS0NNdxAADfwNL1iXA4rFGjRqmsrIw3vgAAj2Lp+kB5ebm2b9+u3/zmN2rTpo3rOACA78DSjXHl5eUaMWKE4uPjKVwA8DiWbgw7fvy4Vq9erfT0dDVp0sR1HADAabB0Y5S1VmPGjFHz5s0pXACIESzdGFRaWqoFCxZo4sSJSkhIcB0HAFBJLN0Y9NRTT6lXr14ULgDEGJauh81dk6dH3s7WodKQJKmBKVNvs0XTJ452nAwAUB0sXY+auyZPw19fd7JwJWnv2kX659HWmrsmz2EyAEB1UboelbVgq0JhK0kKlxap8MMX1eA/Bqoiqb6yFmx1nA4AUB2UrkftLTwqSbLhch3fu0UNe177rd8DAMQWStdj5q7JU+/H/yUrqbz4Cx1a9IySL7hIcUlnnLxP85RkdwEBANVWqQupjDGXS/qjpHhJM621j0c0VQDNXZOncW9lq/Dol+dww6VFChcfUMpPbpMx5uT9EuKMhg9s6yomAKAGTrt0jTHxkqZJukJSB0mDjTEdIh0sSOauydMDr6w9WbjlRfkqWvayEpu1VlxC0sn7pSQnKOuG/9DPu7ZwFRUAUAOVWboXSfrUWrtTkowxL0u6RtKmSAYLkgdeWXvy69ChfbLhkM7s9xuZ+C9fh2sk7Xp8kKN0AIDaUplzui0k7fnK97knbkMt+OrLf0KH9qlkzT+U0LjFycKVOIcLAH5RmdI1p7jNfutOxtxljFlpjFlZUFBQ82QBMHdNntLnbJAklRXslsLlSvnp7TJx8Sfvk5wQzzlcAPCJypRurqRWX/m+paS937yTtXaGtbaHtbYHb8BfOVkLtupoKKzy4i90ZNMS1TmrxdcKV5ImXdeZc7gA4BOVKd1PJP3AGJNmjEmU9EtJb0U2VjDsLTyq4/u2KXzkkFL63ipjvv50NEyKp3ABwEdOW7rW2nJJ90paIGmzpFettdmRDhYEzZKlY7vXKrHZBV97WZD0ZeGuf+RyR8kAAJFQqdfpWmv/IekfEc4SKEuWLNHlTYpU1HewjobCJ29PTojnkDIA+BTvSOVAeXm5tmzZooy7B2vSdZ3VIiVZRlKLlGQKFwB8jI/2i7L58+fr8OHD+t3vfidJ+nnXFpQsAAQESzeKjhz5DAFaAAAHQklEQVQ5olAopJtuusl1FACAAyzdKHnjjTdUUlKi2267zXUUAIAjlG4U7N69W2lpaerWrZvrKAAAhyjdCHvhhRckSbfccovjJAAA1yjdCPrwww912WWXqVmzZq6jAAA8gAupIuTFF19UQUEBhQsAOImlGwGvvPKKrr/+eiUlJZ3+zgCAwGDp1rI333xTKSkpFC4A4FtYurXEWqs///nPuuuuuyhcAMApsXRryYcffqguXbpQuACA70Tp1pC1VuPHj1fHjh3Vp08f13EAAB7G4eUaqKio0KZNm3TllVeqcePGruMAADyOpVtN4XBYo0ePVv369dWjRw/XcQAAMYClWw3l5eXauXOnbrnlFqWmprqOAwCIESzdKgqFQho5cqSMMerYsaPrOACAGMLSrYKysjJt2bJF9913n84//3zXcQAAMYalW0nl5eUaPny4GjVqROECAKqFpVsJpaWl+uSTTzRhwgQ1bNjQdRwAQIxi6Z6GtVYZGRlq06YNhQsAqBGW7vc4fPiw5s+fr8zMTMXF8fMJAKBmaJLvMX36dPXv35/CBQDUCpbuKXzxxRd68cUXNWLECNdRAAA+woT7Bmut5s2bp9tvv911FACAz7B0v2Lfvn2aNm2aHn30UddRAAA+xNI94ejRo1q3bp1GjRrlOgoAwKcoXUm7d+9Wenq6BgwYoDPOOMN1HACATwW+dPPy8nTw4EFNnjyZq5QBABEV6JbZunWrpkyZos6dOysxMdF1HACAzwW2dDdv3ixjjDIzM5WQkOA6DgAgAAJZutu2bdPs2bN1wQUXqE4dLuAGAERH4Ep35cqVstbqscceU3x8vOs4AIAACVTp7tu3T++8847atGnDRVMAgKgLzLHVpUuX6owzztDYsWNljHEdBwAQQIGYeyUlJVq1apW6d+9O4QIAnPH90n333XcVFxenIUOGuI4CAAg4Xy/dsrIy7d+/XwMGDHAdBQAA/y7dOXPmKBwO69Zbb3UdBQAAST4t3cLCQjVs2FD9+/d3HQUAgJN8V7ovvPCC4uLidPPNN7uOAgDA1/iqdDdt2qTu3burffv2rqMAAPAtvrmQ6plnntGWLVsoXACAZ/li6S5cuFDXXXedUlJSXEcBAOA7xfzSnT17tsLhMIULAPC8mF261lo999xzuuWWW/ikIABATIjZpfv222+rdevWFC4AIGbEXGNZa5WVlaX77rtPdevWdR0HAIBKi6mla63VRx99pH79+lG4AICYEzOlW1FRodGjR6tDhw7q2bOn6zgAAFSZpw8vz12Tp6wFW5V3sEQpZft1R+9L1bhxY9exAACoFs8u3blr8pQ+Z4NyvyjWwcV/04Hj8Zq+PqS5a/JcRwMAoFo8W7pZC7bqaCisiuNH1KDrlarTqKmOhsLKWrDVdTQAAKrFs4eX9xYelSTFn9FI8Wc0+tbtAADEGs8u3eYpyVW6HQAAr/Ns6Q4f2FbJCfFfuy05IV7DB7Z1lAgAgJrx7OHln3dtIenLc7t7C4+qeUqyhg9se/J2AABijWdLV/qyeClZAIBfePbwMgAAfkPpAgAQJZQuAABRQukCABAllC4AAFFC6QIAECWULgAAUULpAgAQJZQuAABRQukCABAllC4AAFFC6QIAECWULgAAUULpAgAQJZQuAABRYqy1tf+gxhRIyqn1Bw6GsyUdcB0CknguvITnwjt4Lk7tfGttk9PdKSKli+ozxqy01vZwnQM8F17Cc+EdPBc1w+FlAACihNIFACBKKF3vmeE6AE7iufAOngvv4LmoAc7pAgAQJSxdAACihNIFACBKKF0PMcZcbozZaoz51BjzkOs8QWWMecYYs98Ys9F1lqAzxrQyxiwyxmw2xmQbY+53nSmojDF1jTErjDHrTjwXj7jOFIs4p+sRxph4SdskXSYpV9InkgZbazc5DRZAxpi+kkokzbbWdnKdJ8iMMedKOtdau9oY00DSKkk/5+9F9BljjKR61toSY0yCpA8l3W+t/bfjaDGFpesdF0n61Fq701pbJullSdc4zhRI1tqlkg66zgHJWrvPWrv6xNfFkjZLauE2VTDZL5Wc+DbhxC9WWxVRut7RQtKer3yfK/5xAU4yxqRK6irpY7dJgssYE2+MWStpv6SF1lqeiyqidL3DnOI2fooEJBlj6kt6Q9ID1trDrvMElbU2bK3tIqmlpIuMMZx+qSJK1ztyJbX6yvctJe11lAXwjBPnD9+Q9KK1do7rPJCstYWSFku63HGUmEPpescnkn5gjEkzxiRK+qWktxxnApw6cfHOLEmbrbVTXOcJMmNME2NMyomvkyX1l7TFbarYQ+l6hLW2XNK9khboy4tFXrXWZrtNFUzGmJckLZfU1hiTa4z5f64zBVhvSb+WdIkxZu2JX1e6DhVQ50paZIxZry9HwkJr7TzHmWIOLxkCACBKWLoAAEQJpQsAQJRQugAARAmlCwBAlFC6AABECaULAECUULoAAETJ/we3r3CzMdKnkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAHVCAYAAADYaHMGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xl41eWh7fH1JhASIMzIrIAgERTkmFKc8FRF0CvDg3WorXWsPUdpOSAgGESZRAgUh+otnmpVHDg4oCAFRKCOwCXIPASTQGUMEQhJIHPe+wepJ2Igwx7e/dv7+3mePGT4uffqryGLtacYa60AAIB3RLkOAAAAaobyBgDAYyhvAAA8hvIGAMBjKG8AADyG8gYAwGMobwAAPIbyBgDAYyhvAAA8po7rAGfTokUL27FjR9cxAAAIuFOnTik3N1f79+//3lrbsqrjQ7a8O3bsqJSUFNcxAAAIqNzcXH355ZcaMGCAoqOj/1md/4abzQEAcGTNmjV67rnndNNNNykqqvqVTHkDAOBAenq66tevr6SkpBr/t5Q3AABBtnr1ar3//vvq2bOnjDE1/u9D9j5vAADC0ddff62EhAT94he/qPVlsLwBAAiSzz//XBs2bFCbNm18uhyWNwAAQbB48WJdeeWV6tevn8+XxfIGACDANm3apMzMTDVv3twvl0d5AwAQQO+8846aN2+uBx980G+XSXkDABAg+/btk7VWHTp08OvlUt4AAATAK6+8omPHjumuu+7y+2VT3gAA+NmxY8d03nnnqVevXgG5fMobAAA/+vOf/6xt27Zp0KBBAbsOnioGAICfHD58WJdddpmuvvrqgF4PyxsAAB9ZazVz5kzt3Lkz4MUtsbwBAPCJtVZbt27VsGHD1KVLl6BcJ8sbAAAfTJ8+XSdPngxacUssbwAAaqWsrEyffPKJhg8frkaNGgX1ulneAADUwpw5c9S6deugF7fE8gYAoEZKS0v1+uuva+TIkYqKcrOBWd4AANTAvHnzdOWVVzorbonlDQBAtRQVFSk5OVmPP/64jDFOs7C8AQCogrVWn3zyie666y7nxS1R3gAAnFN+fr5GjRqlG2+8UZ06dXIdRxLlDQDAWRUUFGjz5s0aPXq0YmJiXMf5AeUNAEAlcnNzNWrUKCUkJKhdu3au4/wID1gDAOAMJ06c0N69ezVp0iQ1adLEdZyfYHkDAFBBdna2xo8fr/bt26tly5au41SK5Q0AQLmjR49q3759evrpp0Nycf8LyxsAAJ2+j3vixInq3LlzSBe3xPIGAECHDx9Wenq6Zs+erdjYWNdxqsTyBgBEtJKSEs2cOVO9e/f2RHFLLG8AQAT75z//qW+++UazZ88OiVdOqy6WNwAgIllr9fbbb+vmm2/2VHFLLG8AQATavXu3/vGPf2j8+PGuo9QKyxsAEFHKysq0Zs0a3X///a6j1BrLGwAQMbZu3aolS5Zo3LhxrqP4hOUNAIgIJ06c0OHDhzV27FjXUXxGeQMAwt769ev17LPPqn///oqK8n71+eV/gTFmoDEm1RiTZoz5yW0RxphRxpgdxpgtxpiVxpgL/HG9AABUZffu3YqLi9PEiRNdR/Ebn8vbGBMt6UVJN0nqLulXxpjuZxy2UVKitbanpPckzfT1egEAqMqXX36p999/Xz169PDc08HOxR8PWOsjKc1amyFJxpj5koZI2vGvA6y1qyscv1bSb/xwvQAAnNVXX32lzp0766qrrgqr4pb8c7N5O0n7Kny8v/xzZ/OApKWVfcEY85AxJsUYk5KVleWHaACASLR27VqtX79ebdu2DbvilvyzvCs7K7bSA435jaRESddW9nVr7cuSXpakxMTESi8DAIBzWbJkiX72s5+pb9++rqMEjD+W935JHSp83F7SwTMPMsbcIClJ0mBrbaEfrhcAgB/ZuXOn9u/fr/POO891lIDyR3mvl9TVGNPJGBMj6U5JiyoeYIzpLWmuThf3ET9cJwAAP7JgwQLFxMTo97//vesoAedzeVtrSyQNl7Rc0k5JC6y1240xk40xg8sPS5bUUNK7xphNxphFZ7k4AABq7PDhwyosLNSFF17oOkpQGGtD867lxMREm5KS4joGACDEvfbaa+revbv69OnjOorPjDEbrLWJVR3n/ZeZAQBErOzsbDVt2jQsirsm+MUkAABPeumll9StWzcNGTLEdZSgo7wBAJ5z+PBh9erVS1dddZXrKE5wszkAwFNmz56trVu3RmxxSyxvAIBHWGu1adMmDRs2TJ06dXIdxymWNwDAE5KTk5WXlxfxxS2xvAEAIc5aqyVLlujhhx9Ww4YNXccJCSxvAEBIe+GFF9S6dWuKuwKWNwAgJJWWluqVV17RH/7wh7D8zWC+YHkDAELS/PnzdeWVV1LclWB5AwBCSnFxsaZPn64JEyYoKoqNWRnOCgAgZFhrtWrVKt11110U9zlwZgAAIaGgoECjRo1Sv3791KVLF9dxQhrlDQBwrqCgQNu2bdOoUaMUFxfnOk7Io7wBAE6dPHlSo0ePVseOHdWhQwfXcTyBB6wBAJzJzc3Vnj179MQTT6hFixau43gGyxsA4EROTo7GjRunNm3aqFWrVq7jeArLGwAQdMePH9fevXs1depUNW3a1HUcz2F5AwCCKj8/XxMmTFCnTp0o7lpieQMAgiYzM1M7duzQnDlzFBMT4zqOZ7G8AQBBUVZWplmzZqlPnz4Ut49Y3gCAgNu/f7++/vprzZw5k9cq9wOWNwAg4ObPn6/BgwdT3H7C8gYABEx6erqWLVum0aNHu44SVljeAICAKC0t1dq1a/XQQw+5jhJ2WN4AAL/bsWOH3n//fT3xxBOuo4QlljcAwK+ys7N1+PBhPf74466jhC3KGwDgN998841mzZqlX/ziF4qOjnYdJ2xR3gAAv0hNTVVsbKymTJnCo8oDjPIGAPhs7dq1evfdd3XxxRdT3EFAeQMAfPL111+rQ4cOSkpKoriDhPIGANTahg0btGbNGrVr147iDiLKGwBQK0uXLlXr1q316KOPuo4ScShvAECNpaenKyMjQ+3atXMdJSJR3gCAGvnggw9UWFioRx55xHWUiEV5AwCq7ejRo8rJyVH37t1dR4lovDwqAKBa5s2bp44dO+ree+91HSXisbwBAFXKyclRo0aNdM0117iOArG8AQBVePnll9WhQwcNGTLEdRSUo7wBAGd16NAhXXrppbriiitcR0EF3GwOAKjUc889p02bNlHcIYjlDQD4iQ0bNmjo0KG64IILXEdBJVjeAIAfefbZZ5Wbm0txhzCWNwBAkmSt1aJFi/TQQw+pfv36ruPgHFjeAABJpx9V3rp1a4rbA1jeABDhysrKNHfuXP3Hf/wHvxnMI1jeABDh3n//fV111VUUt4ewvAEgQpWUlGjatGlKSkpSnTrUgZewvAEgAllr9eWXX+rOO++kuD2I8gaACFNUVKRRo0bp8ssvV7du3VzHQS3wzy0AiCAFBQXatWuXRowYofj4eNdxUEssbwCIEKdOndKYMWPUpk0bdezY0XUc+IDlDQAR4OTJk0pPT9f48ePVqlUr13HgI5Y3AIS5vLw8jR07Vq1atVLbtm1dx4EfsLwBIIxlZ2crPT1dU6ZMUbNmzVzHgZ+wvAEgTBUVFSkpKUmdO3emuMMMyxsAwlBWVpY2btyoZ599VnXr1nUdB37G8gaAMGOt1Z/+9CddccUVFHeYYnkDQBg5ePCgVq1apenTp7uOggBieQNAGFmwYIFuvfVW1zEQYCxvAAgDe/fu1cKFCzVy5EjXURAELG8A8LiSkhKtW7dOjzzyiOsoCBKWNwB42K5du/T2229r8uTJrqMgiFjeAOBRR48e1ZEjRzRx4kTXURBklDcAeNDmzZs1a9YsXXPNNfw+7ghEeQOAx+zatUsxMTGaNm2ajDGu48AByhsAPCQlJUXz589XQkKCoqL4ER6p+H8eADxizZo1atWqlZ588kkWd4SjvAHAA7Zu3arPP/9c7du3p7hBeQNAqFuxYoXi4+P12GOPUdyQRHkDQEjbv3+/duzYoY4dO7qOghBCeQNAiProo4+UlZWlESNGuI6CEEN5A0AIysnJ0bFjx9S7d2/XURCCeGY/AISYd955R+edd57uu+8+11EQoljeABBC8vLyVL9+fV1//fWuoyCEsbwBIES8+uqratGihYYMGeI6CkIc5Q0AIWD//v3q0aOHfv7zn7uOAg/gZnMAcOzPf/6zNm7cSHGj2vxS3saYgcaYVGNMmjFmXCVf72eM+cYYU2KM+aU/rhMAwsG6des0dOhQDRo0yHUUeIjP5W2MiZb0oqSbJHWX9CtjTPczDvtO0r2S3vb1+gAgXLzwwgvKzc1V+/btXUeBx/jjPu8+ktKstRmSZIyZL2mIpB3/OsBau7f8a2V+uD4A8DRrrRYuXKjf/e53io2NdR0HHuSPm83bSdpX4eP95Z+rMWPMQ8aYFGNMSlZWlh+iAUDoee2119S6dWuKG7Xmj+Vd2avk29pckLX2ZUkvS1JiYmKtLgMAQlVZWZleeuklPfzww/wubvjEH989+yV1qPBxe0kH/XC5ABBWFi9erCuvvJLihs/88R20XlJXY0wnY0yMpDslLfLD5QJAWCgtLdWTTz6pgQMH6t/+7d9cx0EY8Lm8rbUlkoZLWi5pp6QF1trtxpjJxpjBkmSM+ZkxZr+k2yTNNcZs9/V6AcALysrKtHbtWt1xxx2qV6+e6zgIE355hTVr7d8l/f2Mz02s8P56nb45HQAiRnFxsR577DFNmDBBzZo1cx0HYYSXRwWAACgsLFRqaqoeeeQRiht+x6MmAMDPCgoKNGbMGDVv3lwXXnih6zgIQyxvAPCjU6dOKT09XWPHjlW7drV6yQugSixvAPCTU6dOaezYsWrRogUveYqAYnkDgB/k5OQoNTVVU6ZMUdOmTV3HQZhjeQOAj0pKSpSUlKTOnTtT3AgKljcA+ODYsWNau3at5syZozp1+JGK4GB5A0AtWWs1Z84cXX311RQ3gorvNgCohcOHD2vp0qWaMmWK6yiIQCxvAKiF9957T7fffrvrGIhQLG8AqIHvvvtO8+fP19ixY11HQQRjeQNANRUXF2v9+vUaMWKE6yiIcJQ3AFTDt99+q4kTJ+rWW2/lt4PBOcobAKqQlZWlrKwsHpyGkEF5A8A5bNu2TbNmzVLfvn15OhhCBuUNAGexc+dO1alTR9OnT1dUFD8uETr4bgSASmzevFlvv/22LrroIoobIYfvSAA4w9q1a9W4cWNNnjyZ4kZI4rsSACpITU3VqlWrdMEFF8gY4zoOUCnKGwDKrVq1SsYYjR8/nuJGSKO8AUBSZmamNm/erIsuuojiRsjjeQ8AIt6SJUvUvHlzjRw50nUUoFpY3gAi2qlTp5SZmam+ffu6jgJUG8sbQMRasGCBGjZsqPvvv991FKBGWN4AIlJeXp5iY2N18803u44C1BjLG0DEef3119WwYUPdeuutrqMAtUJ5A4go3333nXr06KHExETXUYBa42ZzABFj7ty52rBhA8UNz2N5A4gIa9as0eDBg9WmTRvXUQCfsbwBhL25c+cqJyeH4kbYYHkDCFvWWr333nu69957Va9ePddxAL9heQMIW2+//bZatWpFcSPssLwBeNqHGw8oeXmqDmTnK9oYlVqrto1j1f34V5o7PUnR0dGuIwJ+x/IG4Fkfbjyg8R9s1YHsfElSqbWSpPSNX+nTo421eMthl/GAgKG8AXhW8vJU5ReX/vCxtWXK/uIt1evQQ2rZRcnLUx2mAwKHm80BeNbB8sUtSbasVEWZ6aqfcJWi6sb+5OtAOGF5A/Cstk3iJJ0u7uzPXledRi0V07LjT74OhBvKG4BnjRnQTbFRVsVH96thrwGKbtD0h6/F1Y3WmAHdHKYDAofyBuBZN3Vvoa57P1Tblk1Vt1k7RRsjSWrXJE7Th12qob3bOU4IBAb3eQPwpPz8fKWlpemlGU+qQ4cOruMAQcXyBuA5hYWFGjt2rJo1a0ZxIyKxvAF4Sl5enrZu3aopU6aoSZMmruMATrC8AXhGWVmZkpKS1KVLF4obEY3lDcATjh8/rs8//1x/+tOfeMlTRDyWNwBPeOGFF3TttddS3IBY3gBC3JEjR/Thhx9q4sSJrqMAIYPlDSCkLVy4UHfddZfrGEBIYXkDCEkHDhzQ3/72N02YMMF1FCDksLwBhJyioiKlpKRo9OjRrqMAIYnyBhBS0tPT9fjjj2vw4MGKjY11HQcISZQ3gJCRmZmp77//XtOnT5cpf51yAD9FeQMICTt37lRycrISExNVt25d13GAkEZ5A3Bu165dioqK0owZM3geN1ANlDcAp7Zv36433nhDXbt2pbiBaqK8ATizfv16xcbGaurUqYqK4scRUF38bQHgxJ49e7R06VJ17tyZ4gZqiL8xAILus88+08mTJ/XEE0/wqHKgFihvAEF17NgxpaSkqEePHhQ3UEu8PCqAoFm2bJnq16+vRx991HUUwNNY3gCCorCwUAcOHFC/fv1cRwE8j+UNIOAWLlyoqKgoPfDAA66jAGGB8gbgdx9uPKDk5ak6mJ2vVnHS9U1yNG3EPa5jAWGD8gbgVx9uPKDxH2xVfnGp8ravVl5UtE70/IV+tvGAhvZu5zoeEBa4zxuAXyUvT1V+camKsw8rpsX5anBxP+UXlyp5earraEDYoLwB+NXB7Hzlbv5ERYfTFNPqwh99HoB/cLM5AL9qePxblXT5maIbNP3R59s2iXOUCAg/LG8AfvPqq69q6CXN1LBJix99Pq5utMYM6OYoFRB+WN4A/GL+/Pn6zW9+o5iYGF1e4dHmbZvEacyAbjxYDfAjyhuAz9599121atVKMTExkqShvdtR1kAAUd4Aas1aq+eff14PP/yw6tat6zoOEDG4zxtAra1evVo///nPKW4gyChvADVWVlamJ598UomJierbt6/rOEDEobwB1EhZWZm2bNmiYcOGqVGjRq7jABGJ8gZQbSUlJRo3bpxatGihXr16uY4DRCwesAagWoqLi5Wenq77779f7du3dx0HiGgsbwBVKioq0tixYxUTE6OEhATXcYCIx/IGcE6FhYVKTU3VyJEjdf7557uOA0AsbwDnUFxcrNGjR6tp06YUNxBC/LK8jTEDJT0nKVrSX621z5zx9XqS3pB0uaSjku6w1u71x3UD8J+eTy5TTmGpJKmsqEB1j+3RN89NVePGjR0nA1CRz8vbGBMt6UVJN0nqLulXxpjuZxz2gKTj1toukuZImuHr9QLwr4rFba1V9hfzVBTfWtf8aY3jZADO5I+bzftISrPWZlhriyTNlzTkjGOGSHq9/P33JF1vjDF+uG4AfvDhxgP/u7gLT+pU6ldqet0Dim7Q9IfPAwgd/ijvdpL2Vfh4f/nnKj3GWlsi6YSk5mdekDHmIWNMijEmJSsryw/RAFTlw40HNP6DrT98nPvNEsVe0EvG8JAYIFT5429nZQva1uIYWWtfttYmWmsTW7Zs6YdoAKqSvDxV+cWlKj11QrnffKzGV9yu6Lh417EAnIM/ynu/pA4VPm4v6eDZjjHG1JHUWNIxP1w3AB8dzM6XtVanvl2nBj2u+8nXG9WLdpAKwLn4o7zXS+pqjOlkjImRdKekRWccs0jSPeXv/1LSKmvtT5Y3gOBrEXVSJ758W/G9blRUvfo/+lqjetHaMmmgo2QAzsbnp4pZa0uMMcMlLdfpp4q9aq3dboyZLCnFWrtI0iuS5hlj0nR6cd/p6/UC8F1hYaEGtC5QzjW3q7DC5+PqRmv6sEs1tPeZD18BEAr88jxva+3fJf39jM9NrPB+gaTb/HFdAPxj7969ev755zV79mwlbjqo5OWpOpidr7ZN4jRmQDeKGwhhJlRvvU5MTLQpKSmuYwBh6dChQzpw4IB69uypmJgY13EAlDPGbLDWJlZ1HM8FASLM7t27NWvWLPXq1YviBjyKX0wCRJDU1FRJ0owZM1SnDn/9Aa9ieQMRYvfu3frb3/6mLl26UNyAx1HeQAT45ptvJElPP/20oqN53jbgdZQ3EOYOHDigjz76SF27dlVUFH/lgXDAbWdAGPvqq68UFxenp556SvwuICB88M9wIEzl5uZqzZo16t27N8UNhBmWNxCGPv30U0nS6NGjHScBEAgsbyDMFBcX67vvvtMNN9zgOgqAAGF5A2Fk0aJFKigo0P333+86CoAAYnkDYSI3N1d16tTR7bff7joKgABjeQNhYP78+SouLtbdd9/tOgqAIKC8AY9LS0tT9+7d1bNnT9dRAAQJN5sDHvThxgO66plVavl/Run/PPGaMkqbu44EIIhY3oDHfLjxgMZ/sFXH0zcq9sJEFdZvrPEfbJUkfgc3ECFY3oDHJC9PVdamT1VWmK/o+o0lSfnFpUpenuo4GYBgYXkDHrP766VqkHC1THTdH33+YHa+o0QAgo3lDXjIRx99pPNatvxJcUtS2yZxDhIBcIHyBjzAWqvnnntOAwYM0LThv1Zc3R//Ws+4utEaM6Cbo3QAgo2bzQEP+Oqrr5SYmKjY2NgfHpSWvDxVB7Pz1bZJnMYM6MaD1YAIQnkDIcxaq0mTJukPf/iDmjf/36eDDe3djrIGIhjlDYSosrIy7dixQ4MGDfpRcQMA93kDIai0tFRJSUlq0KCBLr/8ctdxAIQYljcQYkpKSpSRkaG7775bnTp1ch0HQAhieQMhpLi4WOPGjZMxRt27d3cdB0CIYnkDIaKoqEi7du3SH//4R51//vmu4wAIYSxvIASUlpZqzJgxatSoEcUNoEosb8CxU6dOad26dZo6dari4+NdxwHgASxvwLFJkyYpISGB4gZQbSxvwJHc3Fx9/PHHmj59uqKi+Hc0gOrjJwbgyF/+8hf179+f4gZQYyxvIMiOHTumN954Q2PGjHEdBYBH8U9+IIistfr444913333uY4CwMMobyBIMjMzlZSUpN/+9rdq3Lix6zgAPIzyBoIgPz9fmzdvVlJSkusoAMIA5Q0E2L59+zRu3DjdcMMNatCgges4AMIA5Q0E0MGDB5WVlaWZM2fyqHIAfsNPEyBA0tLSNHPmTF1yySWqV6+e6zgAwghPFQMCIDU1VdZaJScnq27duq7jAAgzLG/AzzIyMvTKK6+oS5cuFDeAgGB5A360efNm1a1bV9OnT1d0dLTrOADCFMsb8JOsrCy9++67SkhIoLgBBBTLG/CDtWvXylqrKVOmyBjjOg6AMMfyBnyUn5+vL774Qn379qW4AQQFyxvwwerVq3Xy5El+yQiAoGJ5A7VUUlKivXv36pZbbnEdBUCEYXkDtfD3v/9dx48f57eDAXCC5Q3U0IkTJxQdHa1f//rXrqMAiFAsb6AG3nvvPeXm5rK4AThFeQPVlJqaqosvvlg9evRwHQVAhKO8gWp46623ZIzRXXfd5ToKAFDeQFVWrVqlm266Sc2aNXMdBQAkUd7AOS1YsEB169aluAGEFMobOIs333xTt99+u2JiYlxHAYAf4aliQCWWLVumFi1aUNwAQhLLG6jAWqvnn39eDz74oBo0aOA6DgBUiuUNVJCSkqLLLruM4gYQ0ihvQPrh13l26NBB1157res4AHBOlDcinrVWu3fv1oABA9S6dWvXcQCgSpQ3IlpZWZkmTJggY4z69OnjOg4AVAsPWEPEKi0t1Z49e/SrX/1KF110kes4AFBtLG9EpJKSEj322GMqLi7WJZdc4joOANQIyxsRp7i4WNu3b9eIESPUoUMH13EAoMZY3ogo1lqNHTtW8fHxFDcAz2J5I2IUFBToiy++0NSpU3keNwBPY3kjYkyZMkUXX3wxxQ3A81jeCHt5eXn64IMPNHXqVBljXMcBAJ+xvBH2XnnlFQ0cOJDiBhA2WN4IW9nZ2frv//5vjRkzxnUUAPArljfCkrVWS5cu1YMPPug6CgD4HeWNsHPkyBGNHTtWd955p5o2beo6DgD4HeWNsHLy5Elt27ZNEydO5D5uAGGL8kbYOHDggB577DH169dP8fHxruMAQMBQ3ggLBw8eVFZWlpKTk1WnDo/DBBDeKG943t69ezVz5kxdfPHFiouLcx0HAAKOiQJPS0tLU3FxsWbOnKmYmBjXcQAgKFje8Kx9+/bpL3/5i7p27UpxA4goPpW3MaaZMWaFMebb8j8rfV6OMWaZMSbbGPOxL9cH/Mu2bdt0/PhxzZgxg/u4AUQcX5f3OEkrrbVdJa0s/7gyyZLu9vG6AEmnXzntnXfe0SWXXKLo6GjXcQAg6Hwt7yGSXi9//3VJQys7yFq7UlKuj9cFKCUlRZs3b9a0adMUFcW9PgAik68//VpZaw9JUvmf5/lyYcaYh4wxKcaYlKysLB+jIdwUFRVp9erV6tevn+soAOBUlXcWGmM+ldS6ki8l+TuMtfZlSS9LUmJiovX35cO7vvjiC2VlZfFLRgBA1Shva+0NZ/uaMSbTGNPGWnvIGNNG0hG/pgMklZSUaM+ePbr7bh42AQCS78/zXiTpHknPlP/5kc+JgAo++eQTHThwQPfdd5/rKAAQMny9z/sZSf2NMd9K6l/+sYwxicaYv/7rIGPMF5LelXS9MWa/MWaAj9eLCHD8+HFFRUVR3ABwBp+Wt7X2qKTrK/l8iqQHK3x8jS/Xg8jz4YcfKisrS7/73e9cRwGAkMOrWyDk7NixQwkJCRo6tNJnHgJAxOOJsggpCxYsUEpKihISElxHAYCQxfJGyFixYoVuvPFGNWnSxHUUAAhplDdCwsKFC1VSUkJxA0A1UN5wbt68ebrtttsUGxvrOgoAeAL3ecOp1atXq2nTphQ3ANQAyxvOPP/887rnnnvUuHFj11EAwFMobwSdtVZbtmzRJZdcQnEDQC1wszmCylqrp59+Wk2bNtV1113nOg4AeBLljaCx1mrPnj267rrrdP7557uOAwCeRXkjKMrKyjRx4kTl5+friiuucB0HADyN8kbAlZWVKSMjQ7fddpt69OjhOg4AeB7ljYAqLS3V+PHjVVBQoJ49e7qOAwBhgUebI2BKSkq0ZcsWDR8+XB06dHAdBwDCBssbAWGt1bhx4xQfH09xA4Cfsbzhd4WFhVq1apUmT56HctXBAAANYUlEQVSs+vXru44DAGGH5Q2/mz59unr06EFxA0CAsLzhN6dOndL8+fP15JNPyhjjOg4AhC2WN/zmtdde080330xxA0CAsbzhs5ycHL344osaP3686ygAEBFY3vCJtVaffPKJHnroIddRACBiUN6otaNHj+rRRx/VrbfequbNm7uOAwARg/JGreTl5WnHjh166qmnuI8bAIKM8kaNHTp0SOPGjVOfPn3UqFEj13EAIOJQ3qiRQ4cO6ciRI5o5c6bq1avnOg4ARCTKG9W2b98+PfPMM+rWrRsvwAIADvFUMVRLRkaG8vPzWdwAEAJY3qjSoUOH9OKLL+qiiy6iuAEgBLC8cU47duxQQUGBZsyYoTp1+HYBgFDA8sZZ5eXl6a233lKvXr0obgAIIfxERqU2btyo77//XtOmTXMdBQBwBpY3fqKkpEQrV67UDTfc4DoKAKASLG/8yNdff63vvvtOo0ePdh0FAHAWLG/8oKSkRBkZGbrjjjtcRwEAnAPLG5KklStXKi0tTb///e9dRwEAVIHlDR09elRRUVEUNwB4BMs7wi1evFj79u3Tww8/7DoKAKCaKO8ItnXrViUkJGjQoEGuowAAaoCbzSPUBx98oPXr16tr166uowAAaojlHYGWL1+u66+/Xo0bN3YdBQBQC5R3hFmyZIlycnIobgDwMMo7grz55psaNmwYv4sbADyO+7wjxFdffaX4+HiKGwDCAOUdAV544QV169ZNQ4YMcR0FAOAHlHeY27lzp7p166YWLVq4jgIA8BPKO0xZa/XMM8+oTp06uvHGG13HAQD4EeUdhqy12rdvn66++mqexw0AYYjyDjPWWk2aNEnHjx/X1Vdf7ToOACAAeKpYGCkrK9OePXs0bNgw9ezZ03UcAECAsLzDRFlZmZKSkpSXl0dxA0CYY3mHgZKSEm3cuFHDhw9Xu3btXMcBAAQYyzsMTJgwQfHx8RQ3AEQIlreHFRUVadmyZXrqqacUGxvrOg4AIEhY3h6WnJysSy+9lOIGgAjD8vag/Px8zZs3T48//riMMa7jAACCjOXtQW+++aZuueUWihsAIhTL20Nyc3M1Z84cTZw40XUUAIBDLG+PKCsr06pVq/Sf//mfrqMAAByjvD3g+PHjGjlypG655Ra1bNnSdRwAgGOUd4jLycnRrl27NGnSJEVHR7uOAwAIAZR3CMvMzNS4cePUq1cvNWnSxHUcAECI4AFrIerQoUPKyspScnKy6tev7zoOACCEsLxD0MGDB/X000+rS5cuatCgges4AIAQw/IOMXv37lVubq6Sk5N55TQAQKVY3iHk+++/13PPPadu3bpR3ACAs2J5h4jU1FRlZ2crOTlZderwfwsA4OxY3iGgoKBAb7zxhi6//HKKGwBQJZrCsS1btmjfvn2aNm2a6ygAAI9geTtUVlamFStW6KabbnIdBQDgISxvR9atW6fU1FQ9+uijrqMAADyG5e1AcXGxMjIydPfdd7uOAgDwIJZ3kP3jH//Qtm3bNHz4cNdRAAAexfIOoqysLEVHR1PcAACfsLyDZOnSpfr222/1xz/+0XUUAIDHUd5BsGnTJnXr1o1HlQMA/IKbzQNs8eLFWr9+vTp37uw6CgAgTLC8A2jp0qXq16+fGjdu7DoKACCMsLwDZMWKFTp69CjFDQDwO5+WtzGmmaT/kdRR0l5Jt1trj59xzGWS/q+kRpJKJU2z1v6PL9cb6t58800NGTJE8fHxrqMAAMKQr8t7nKSV1tquklaWf3ymU5J+a63tIWmgpGeNMU18vN6QlZKSori4OIobABAwvpb3EEmvl7//uqShZx5grd1trf22/P2Dko5Iaunj9Yakl156Se3atdOtt97qOgoAIIz5Wt6trLWHJKn8z/POdbAxpo+kGEnpZ/n6Q8aYFGNMSlZWlo/RgisjI0OdOnVSmzZtXEcBAIS5KsvbGPOpMWZbJW9DanJFxpg2kuZJus9aW1bZMdbal621idbaxJYtvTHOrbWaNWuW8vPzeR43ACAoqnzAmrX2hrN9zRiTaYxpY609VF7OR85yXCNJSyRNsNaurXXaEGOt1cGDB9WnTx/16NHDdRwAQITw9WbzRZLuKX//HkkfnXmAMSZG0kJJb1hr3/Xx+kKGtVZTp07VoUOH1K9fP9dxAAARxNfyfkZSf2PMt5L6l38sY0yiMeav5cfcLqmfpHuNMZvK3y7z8XqdstYqPT1dQ4cOVWJious4AIAIY6y1rjNUKjEx0aakpLiO8RPWWk2YMEHDhg3T5Zdf7joOACCMGGM2WGurXIW8PGoNlJaWav369Ro+fDiPKgcAOMPLo9bAU089pfj4eIobAOAUy7saiouLtXjxYk2YMEH16tVzHQcAEOFY3tXw7LPPqmfPnhQ3ACAksLzPoaCgQK+++qpGjx4tY4zrOAAASGJ5n9P8+fM1aNAgihsAEFJY3pU4efKkZsyYoUmTJlHcAICQw/I+Q2lpqT777DM98sgjFDcAICRR3hVkZ2drxIgR6t+/v1q1auU6DgAAlaK8y+Xk5Gj37t2aPHmy6tat6zoOAABnRXlLysrK0vjx49W9e3c1a9bMdRwAAM4p4h+wlpmZqcOHD+uZZ55Rw4YNXccBAKBKEb28MzMzNWXKFHXp0kXx8fGu4wAAUC0Ru7z37dunY8eOKTk5WXFxca7jAABQbRG5vE+cOKHZs2crISGB4gYAeE7ELe+0tDQdOnRIs2fPVnR0tOs4AADUWEQt7+LiYr322mvq27cvxQ0A8KyIWd7bt2/X7t27NXXqVNdRAADwSUQsb2utPv30Uw0ePNh1FAAAfBb2y3vDhg3auHGjRowY4ToKAAB+EdbLu6ioSBkZGXrggQdcRwEAwG/Cdnl/+eWXWr9+vUaOHOk6CgAAfhWWy/vw4cOKiorSf/3Xf7mOAgCA34Vdea9YsUJvvfWWrrzySn4fNwAgLIXVzeYbNmxQ165d1b9/f9dRAAAImLBZ3suWLdO6devUsWNH11EAAAiosFjeS5cu1RVXXKGBAwe6jgIAQMB5fnl/9tlnyszMVJMmTVxHAQAgKDy9vN966y3dfPPNuvbaa11HAQAgaDy7vLds2aI6deqoadOmrqMAABBUnizvuXPnqlGjRrrjjjtcRwEAIOg8V9779+9X27ZteVQ5ACBieaq858yZo6ysLA0aNMh1FAAAnPHEA9astcrMzFTv3r3Vu3dv13EAAHAq5Je3tVbTp09XRkaG/v3f/911HAAAnAvp5W2tVVpamgYNGqRLL73UdRwAAEJCSC/vp556SsePH6e4AQCoIGSXd15enh555BGdd955rqMAABBSjLXWdYZKGWOyJP3TdY4aaiHpe9chwhznODg4z4HHOQ48L57jC6y1Las6KGTL24uMMSnW2kTXOcIZ5zg4OM+BxzkOvHA+xyF9nzcAAPgpyhsAAI+hvP3rZdcBIgDnODg4z4HHOQ68sD3H3OcNAIDHsLwBAPAYyhsAAI+hvH1gjGlmjFlhjPm2/M+mlRxzmTFmjTFmuzFmizGGX0JeA9U5x+XHLTPGZBtjPg52Rq8yxgw0xqQaY9KMMeMq+Xo9Y8z/lH99nTGmY/BTel81znM/Y8w3xpgSY8wvXWT0umqc41HGmB3lP4NXGmMucJHTnyhv34yTtNJa21XSyvKPz3RK0m+ttT0kDZT0rDGmSRAzel11zrEkJUu6O2ipPM4YEy3pRUk3Seou6VfGmO5nHPaApOPW2i6S5kiaEdyU3lfN8/ydpHslvR3cdOGhmud4o6REa21PSe9JmhnclP5HeftmiKTXy99/XdLQMw+w1u621n5b/v5BSUckVfnqOfhBledYkqy1KyXlBitUGOgjKc1am2GtLZI0X6fPdUUVz/17kq43xpggZgwHVZ5na+1ea+0WSWUuAoaB6pzj1dbaU+UfrpXUPsgZ/Y7y9k0ra+0hSSr/85wvxG6M6SMpRlJ6ELKFixqdY1RbO0n7Kny8v/xzlR5jrS2RdEJS86CkCx/VOc/wTU3P8QOSlgY0URCE7C8mCRXGmE8lta7kS0k1vJw2kuZJusday7+wK/DXOUaNVLagz3zeaHWOwblxDgOv2ufYGPMbSYmSrg1ooiCgvKtgrb3hbF8zxmQaY9pYaw+Vl/ORsxzXSNISSROstWsDFNWz/HGOUWP7JXWo8HF7SQfPcsx+Y0wdSY0lHQtOvLBRnfMM31TrHBtjbtDpQXCttbYwSNkChpvNfbNI0j3l798j6aMzDzDGxEhaKOkNa+27QcwWLqo8x6iV9ZK6GmM6lX+P3qnT57qiiuf+l5JWWV7Vqaaqc57hmyrPsTGmt6S5kgZba8NjAFhreavlm07f/7dS0rflfzYr/3yipL+Wv/8bScWSNlV4u8x1dq+8Veccl3/8haQsSfk6/S/xAa6zh/qbpJsl7dbpx2AklX9usk7/gJOkWEnvSkqT9P8kdXad2Ytv1TjPPyv/nj0p6aik7a4ze+2tGuf4U0mZFX4GL3Kd2dc3Xh4VAACP4WZzAAA8hvIGAMBjKG8AADyG8gYAwGMobwAAPIbyBgDAYyhvAAA85v8DfZeMrPLj9RIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-2.42991573e-02  5.31785193e-02 -0.00000000e+00  0.00000000e+00\n",
      "  0.00000000e+00  3.05467635e+00  5.11427657e-02 -0.00000000e+00\n",
      "  0.00000000e+00 -1.67773867e-03 -0.00000000e+00  2.07160117e-02\n",
      " -5.98679048e-01]\n",
      "[-2.65336120e-02  5.33883129e-02  0.00000000e+00  0.00000000e+00\n",
      "  0.00000000e+00  3.08207168e+00  5.35141490e-02  0.00000000e+00\n",
      "  6.51399000e-04 -1.69731226e-03  0.00000000e+00  2.01544610e-02\n",
      " -6.07579156e-01]\n"
     ]
    }
   ],
   "source": [
    "# Coefficent Comparison\n",
    "graph.figure(figsize=(8, 8))\n",
    "graph.plot(gt.coef_, lasso.coef_, 'o')\n",
    "graph.plot([-0.6, 3.5], [-0.6, 3.5], '--', color='black', linewidth=0.5)\n",
    "graph.show()\n",
    "\n",
    "graph.figure(figsize=(8, 8))\n",
    "graph.plot(gt.coef_, lasso.coef_, 'o')\n",
    "graph.plot([-0.6, 3.5], [-0.6, 3.5], '--', color='black', linewidth=0.5)\n",
    "graph.xlim([-0.25, 0.25])\n",
    "graph.ylim([-0.25, 0.25])\n",
    "graph.show()\n",
    "\n",
    "print(gt.coef_)\n",
    "print(lasso.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CustomLasso:\n",
    "    def __init__(self, lam):\n",
    "        self.lam = lam\n",
    "        self.coef_ = None\n",
    "        self.history = []\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        self.history = []\n",
    "    \n",
    "        n, p = X.shape\n",
    "        m = np.zeros(p)\n",
    "        \n",
    "        # Coordinate Descent\n",
    "        self.history.append(objective_function(x, y, m, self.lam))\n",
    "        for i in range(100):\n",
    "            for pi in range(p):\n",
    "                m[pi] = opt.fmin(partial_objective_function, m[pi], args=(x, y, m.copy(), self.lam, pi), disp=False)\n",
    "                m[abs(m) < 1e-15] = 0  # Polish\n",
    "                self.history.append(objective_function(x, y, m, self.lam))\n",
    "        \n",
    "        self.coef_ = m.copy()\n",
    "        return self\n",
    "    \n",
    "    def predict(self, X):\n",
    "        return model(X, self.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "885 ms ± 8.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "lasso = CustomLasso(lam=2).fit(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My quick BS version is ~900x slower than sklearn's BUT the concept was proven. That is I can use scipy to force coordinate descent to work since I cannot figure out Tibs crazy fortran code. The extra problem is now I need to run iterations of the Pliable Lasso"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}