06. Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class Overview\n",
    "<img src=\"PDF_slides/class_overview.png\"  width=\"500\">\n",
    "\n",
    "\n",
    "___\n",
    "# Agenda\n",
    " - Numerical Optimization Techniques\n",
    "  - Types of Optimization\n",
    "  - Programming the Optimization\n",
    " - **Whirlwind Lecture Alert**\n",
    "  - Entire classes cover these concepts in expanded form\n",
    "  - But we can cover them in one lecture to get a good intuition!\n",
    "  - And then you can look over this even more for better understanding.\n",
    "  - If you feel confused after this lecture, that's okay. These are not easy the first time you see them. Keep going, you got this.\n",
    "  \n",
    "<img src=\"PDF_slides/dont-be-nervous-you-got-this-meme.jpg\"  width=\"300\">\n",
    "\n",
    "\n",
    "    \n",
    "___\n",
    "\n",
    "# Last Time\n",
    "\n",
    "|Description| Equations, Derivations, Hessian Calculations, and Miscellaneous|\n",
    "|-----------|--------|\n",
    "| Sigmoid Definition | $$ p(y^{(i)}=1\\text{ | }\\mathbf{x}^{(i)},\\mathbf{w})=\\frac{1}{1+\\exp{(-\\mathbf{w}^T \\mathbf{x}^{(i)}})}$$ |\n",
    "| Log Likelihood | $$ l(\\mathbf{w}) = \\sum_i \\left( y^{(i)} \\ln [g(\\mathbf{w}^T \\mathbf{x}^{(i)})] + (1-y^{(i)}) (\\ln [1 - g(\\mathbf{w}^T \\mathbf{x}^{(i)})])  \\right)  $$ |\n",
    "| Vectorized Gradient | $$gradient =\\frac{1}{M}\\sum_{i=1}^M (y^{(i)}-g(\\mathbf{w}^T\\mathbf{x}^{(i)}))\\mathbf{x}^{(i)}$$ |\n",
    "| Regularization | $$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\left[\\underbrace{\\nabla l(\\mathbf{w})_{old}}_{\\text{old gradient}} - C \\cdot 2\\mathbf{w} \\right]$$|\n",
    "\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=\"PDF_slides/BtJXjJcCAAE7QOB.jpg\"  width=\"300\">\n",
    "\n",
    "# More Advanced  Optimization for Machine Learning\n",
    "\n",
    "## Problem One: finding the right eta\n",
    "From previous notebooks, we know that the logistic regression update equation is given by:\n",
    "\n",
    "$$ \\underbrace{w_j}_{\\text{new value}} \\leftarrow \\underbrace{w_j}_{\\text{old value}} + \\eta \\underbrace{\\left[\\left(\\sum_{i=1}^M (y^{(i)}-g(\\mathbf{w}^T\\mathbf{x}^{(i)}))x^{(i)}_j\\right) - C \\cdot 2w_j \\right]}_{\\nabla l(w)}$$\n",
    "\n",
    "Which can be made into more generic notation by denoting the objective function as $l(\\mathbf{w})$ and the gradient calculation as $\\nabla$:\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\nabla l(\\mathbf{w})$$\n",
    "\n",
    "One problem is that we still need to set the value of $\\eta$, which can drastically change the performance of the optimization algorithm. If $\\eta$ is too large, the algorithm might be unstable. If $\\eta$ is too small, it might take a long time (i.e., many iterations) to converge.\n",
    "\n",
    "\n",
    "<img src=\"PDF_slides/batch.gif\"  width=\"400\">\n",
    "\n",
    "## Optimizing Logistic Regression via Line Search\n",
    "We can solve this issue by performing a line search for the best value of $\\eta$ along the direction of the gradient.\n",
    "\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\underbrace{\\eta}_{\\text{best step?}} \\nabla l(\\mathbf{w}) $$\n",
    "\n",
    "$$ \\eta \\leftarrow \\arg\\max_\\eta l(\\mathbf{w}+\\eta\\cdot\\nabla l(\\mathbf{w})) $$\n",
    "\n",
    "<img src=\"PDF_slides/line_search copy.gif\"  width=\"400\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "import numpy as np\n",
    "from sklearn.metrics import accuracy_score\n",
    "from scipy.special import expit\n",
    "\n",
    "ds = load_iris()\n",
    "X = ds.data\n",
    "y = (ds.target>1).astype(int) # make problem binary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary Logistic Regression Object with coefficients:\n",
      "[[-0.19960106]\n",
      " [-0.45578023]\n",
      " [-0.54503119]\n",
      " [ 0.77831275]\n",
      " [ 0.51499827]]\n",
      "Accuracy of:  0.9733333333333334\n",
      "CPU times: user 4.38 ms, sys: 2.56 ms, total: 6.94 ms\n",
      "Wall time: 11 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# from last time, our logistic regression algorithm is given by (including everything we previously had):\n",
    "class BinaryLogisticRegression:\n",
    "    def __init__(self, eta, iterations=20, C=0.001):\n",
    "        self.eta = eta\n",
    "        self.iters = iterations\n",
    "        self.C = C\n",
    "        # internally we will store the weights as self.w_ to keep with sklearn conventions\n",
    "        \n",
    "    def __str__(self):\n",
    "        if(hasattr(self,'w_')):\n",
    "            return 'Binary Logistic Regression Object with coefficients:\\n'+ str(self.w_) # is we have trained the object\n",
    "        else:\n",
    "            return 'Untrained Binary Logistic Regression Object'\n",
    "        \n",
    "    # convenience, private:\n",
    "    @staticmethod\n",
    "    def _add_bias(X):\n",
    "        return np.hstack((np.ones((X.shape[0],1)),X)) # add bias term\n",
    "    \n",
    "    @staticmethod\n",
    "    def _sigmoid(theta):\n",
    "        # increase stability, redefine sigmoid operation\n",
    "        return expit(theta) #1/(1+np.exp(-theta))\n",
    "    \n",
    "    # vectorized gradient calculation with regularization using L2 Norm\n",
    "    def _get_gradient(self,X,y):\n",
    "        ydiff = y-self.predict_proba(X,add_bias=False).ravel() # get y difference\n",
    "        gradient = np.mean(X * ydiff[:,np.newaxis], axis=0) # make ydiff a column vector and multiply through\n",
    "        \n",
    "        gradient = gradient.reshape(self.w_.shape)\n",
    "        gradient[1:] += -2 * self.w_[1:] * self.C\n",
    "        \n",
    "        return gradient\n",
    "    \n",
    "    # public:\n",
    "    def predict_proba(self,X,add_bias=True):\n",
    "        # add bias term if requested\n",
    "        Xb = self._add_bias(X) if add_bias else X\n",
    "        return self._sigmoid(Xb @ self.w_) # return the probability y=1\n",
    "    \n",
    "    def predict(self,X):\n",
    "        return (self.predict_proba(X)>0.5) #return the actual prediction\n",
    "    \n",
    "    \n",
    "    def fit(self, X, y):\n",
    "        Xb = self._add_bias(X) # add bias term\n",
    "        num_samples, num_features = Xb.shape\n",
    "        \n",
    "        self.w_ = np.zeros((num_features,1)) # init weight vector to zeros\n",
    "        \n",
    "        # for as many as the max iterations\n",
    "        for _ in range(self.iters):\n",
    "            gradient = self._get_gradient(Xb,y)\n",
    "            self.w_ += gradient*self.eta # multiply by learning rate \n",
    "            # add bacause maximizing \n",
    "\n",
    "blr = BinaryLogisticRegression(eta=0.1,iterations=50,C=0.001)\n",
    "\n",
    "blr.fit(X,y)\n",
    "print(blr)\n",
    "\n",
    "yhat = blr.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y,yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary Logistic Regression Object with coefficients:\n",
      "[[-1.96717802]\n",
      " [-4.81050143]\n",
      " [-6.74146326]\n",
      " [ 8.75170244]\n",
      " [ 5.12131168]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 6.95 ms, sys: 2.52 ms, total: 9.47 ms\n",
      "Wall time: 10.9 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: divide by zero encountered in log\n",
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/scipy/optimize/optimize.py:2044: RuntimeWarning: invalid value encountered in subtract\n",
      "  r = (xf - nfc) * (fx - ffulc)\n",
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/scipy/optimize/optimize.py:2045: RuntimeWarning: invalid value encountered in subtract\n",
      "  q = (xf - fulc) * (fx - fnfc)\n",
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/scipy/optimize/optimize.py:2046: RuntimeWarning: invalid value encountered in subtract\n",
      "  p = (xf - fulc) * q - (xf - nfc) * r\n",
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/scipy/optimize/optimize.py:2047: RuntimeWarning: invalid value encountered in subtract\n",
      "  q = 2.0 * (q - r)\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# and we can update this to use a line search along the gradient like this:\n",
    "from scipy.optimize import minimize_scalar\n",
    "import copy\n",
    "from numpy import ma # (masked array) this has most numpy functions that work with NaN data.\n",
    "class LineSearchLogisticRegression(BinaryLogisticRegression):\n",
    "    \n",
    "    # define custom line search for problem\n",
    "    def __init__(self, line_iters=0.0, **kwds):        \n",
    "        self.line_iters = line_iters\n",
    "        # but keep other keywords\n",
    "        super().__init__(**kwds) # call parent initializer\n",
    "    \n",
    "    # this defines the function with the first input to be optimized\n",
    "    # therefore eta will be optimized, with all inputs constant\n",
    "    @staticmethod\n",
    "    def objective_function(eta,X,y,w,grad,C):\n",
    "        wnew = w - grad*eta\n",
    "        g = expit(X @ wnew)\n",
    "        # the line search is looking for minimization, so take the negative of l(w)\n",
    "        return -np.sum(ma.log(g[y==1]))-ma.sum(np.log(1-g[y==0])) + C*sum(wnew**2)\n",
    "    \n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        Xb = self._add_bias(X) # add bias term\n",
    "        num_samples, num_features = Xb.shape\n",
    "        \n",
    "        self.w_ = np.zeros((num_features,1)) # init weight vector to zeros\n",
    "        \n",
    "        # for as many as the max iterations\n",
    "        for _ in range(self.iters):\n",
    "            gradient = -self._get_gradient(Xb,y)\n",
    "            # minimization is in opposite direction\n",
    "            \n",
    "            # do line search in gradient direction, using scipy function\n",
    "            opts = {'maxiter':self.line_iters} # unclear exactly what this should be\n",
    "            res = minimize_scalar(self.objective_function, # objective function to optimize\n",
    "                                  bounds=(0,self.eta*10), #bounds to optimize\n",
    "                                  args=(Xb,y,self.w_,gradient,self.C), # additional argument for objective function\n",
    "                                  method='bounded', # bounded optimization for speed\n",
    "                                  options=opts) # set max iterations\n",
    "            \n",
    "            eta = res.x # get optimal learning rate\n",
    "            self.w_ -= gradient*eta # set new function values\n",
    "            # subtract to minimize\n",
    "                \n",
    "            \n",
    "\n",
    "lslr = LineSearchLogisticRegression(eta=1,\n",
    "                                    iterations=5, \n",
    "                                    line_iters=5, \n",
    "                                    C=0.001)\n",
    "\n",
    "lslr.fit(X,y)\n",
    "\n",
    "yhat = lslr.predict(X)\n",
    "print(lslr)\n",
    "print('Accuracy of: ',accuracy_score(y,yhat))         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This performs well, but was not too much faster than previously (this is because $\\eta$ was chosen well in the initial example). \n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Self Test\n",
    "How much computation (i.e., how many multiplies) are required for calculating the gradient of:\n",
    "$$ \\left( \\frac{1}{M}\\left[\\sum_{i=1}^M (y^{(i)}-g(\\mathbf{w}^T\\mathbf{x}^{(i)}))\\cdot \\mathbf{x}^{(i)}\\right] - 2C\\cdot \\mathbf{w}\\right) $$\n",
    "\n",
    "Where $M$ is the number of instance and $N$ is the number of elements in $\\mathbf{w}$.\n",
    "\n",
    "- A: $ M\\cdot N+1$\n",
    "- B: $ (M+1)\\cdot N$\n",
    "- C: $ 2N $ \n",
    "- D: $ 2N-M$ \n",
    "_____\n",
    "\n",
    "# Stochastic Gradient Descent\n",
    "Sometimes the gradient calculation is too computational:\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta\\left( \\frac{1}{M}\\left[\\sum_{i=1}^M (y^{(i)}-g(\\mathbf{w}^T\\mathbf{x}^{(i)}))\\mathbf{x}^{(i)}\\right] - 2C\\cdot \\mathbf{w}\\right) $$\n",
    "\n",
    "Instead, we can approximate the gradient using one instance, this is called stochastic gradient descent (SGD) because the steps can appear somewhat random.\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\underbrace{\\left((y^{(i)}-g(\\mathbf{w}^T\\mathbf{x}^{(i)}))\\mathbf{x}^{(i)}-2C\\cdot \\mathbf{w}\\right)}_{\\text{approx. gradient}} \\text{,   where   } i\\in M$$\n",
    "\n",
    "<img src=\"PDF_slides/SGD.gif\"  width=\"400\">\n",
    "\n",
    "Let's code up the SGD example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary Logistic Regression Object with coefficients:\n",
      "[[-0.12301927]\n",
      " [-0.27029268]\n",
      " [-0.36768569]\n",
      " [ 0.5259726 ]\n",
      " [ 0.32828145]]\n",
      "Accuracy of:  0.88\n",
      "CPU times: user 9.4 ms, sys: 2.45 ms, total: 11.8 ms\n",
      "Wall time: 12.9 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "class StochasticLogisticRegression(BinaryLogisticRegression):\n",
    "    # stochastic gradient calculation \n",
    "    def _get_gradient(self,X,y):\n",
    "        idx = int(np.random.rand()*len(y)) # grab random instance\n",
    "        ydiff = y[idx]-self.predict_proba(X[idx],add_bias=False) # get y difference (now scalar)\n",
    "        gradient = X[idx] * ydiff[:,np.newaxis] # make ydiff a column vector and multiply through\n",
    "        \n",
    "        gradient = gradient.reshape(self.w_.shape)\n",
    "        gradient[1:] += -2 * self.w_[1:] * self.C\n",
    "        \n",
    "        return gradient\n",
    "    \n",
    "    \n",
    "slr = StochasticLogisticRegression(eta=0.005, iterations=500, C=0.001) # take a lot more steps!!\n",
    "\n",
    "slr.fit(X,y)\n",
    "\n",
    "yhat = slr.predict(X)\n",
    "print(slr)\n",
    "print('Accuracy of: ',accuracy_score(y,yhat))      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"PDF_slides/hessian_second.jpg\"  width=\"300\">\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Optimizing with Second Order Derivatives\n",
    "First, let's look at the one dimensional case when we have a function $l(w)$ where w is a scalar. The optimal value of w is given by:\n",
    "\n",
    "$$ w_{optimal} \\leftarrow w_{start} - \\underbrace{[\\frac{\\partial^2}{\\partial w}l(w_{start})]^{-1}}_{\\text{inverse 2nd deriv}}\\underbrace{\\frac{\\partial}{\\partial w}l(w_{start})}_{\\text{derivative}}  $$\n",
    "\n",
    "Note that if $l(w)$ is a quadratic function, this solution converges in a single step!\n",
    "\n",
    "\n",
    "**Aside: an example with the second derivative:**\n",
    "- Say $l(w)=2w^2+4w+5$, and we want to minimize the function. We have that:\n",
    "- $\\frac{\\partial}{\\partial w}l(w)=4w+4$\n",
    "- $\\frac{\\partial^2}{\\partial w}l(w)=4$\n",
    "- Therefore, if we choose $w_{start}=0$, we have:\n",
    "- $\\frac{\\partial}{\\partial w}l(0)=4$\n",
    "- $\\frac{\\partial^2}{\\partial w}l(0)=4$ \n",
    "- So the update becomes\n",
    "- $w \\leftarrow w_{start} - \\frac{1}{4}4 = -1$\n",
    "- The solution is found in one step. This works for any initial value of $w_{start}$. Let's verify that the solution worked graphically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(-1, 2.5, '$\\\\leftarrow$found minimum')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwdklEQVR4nO3deVxVdf7H8df3ACqLIHJRA3PBpdQoU8jUDJerllnSjJm5lFlji79MbSxbRitzosyh8Tf2s2lzSqfFStRyC03cC7dc02xMMzNEEBdAhPP9/XGNaUG5wL2cey6f519zb9573l/uY94czv2e71dprTVCCCFsx7A6gBBCiMqRAhdCCJuSAhdCCJuSAhdCCJuSAhdCCJuSAhdCCJsKrO4DHjlypLoPWSEOh4Ps7GyrY1SZv4wDZCy+yF/GAfYYS0xMTJnPyxm4EELYlBS4EELYlBS4EELYlBS4EELYlBS4EELYlBS4EELYlBS4EELYlC0KXO/eirnkQ6tjCCFEhemzhZjvvYbO+tHj722TAt+GXjAXnZdrdRQhhKgQvWkdesUiOJHj8fe2RYGr63pDSQl6w0qrowghRIXotcuhYSy0auvx97ZHgTdqDC3botemIxsICSHsQv94GPbvQV3nRCnl8fe3RYHD+bPwn36Ab3ZbHUUIIdyi134GAQGozj298v72KfCErlAn2PUDEUIIH6eLz7ku+8YnoiIivXIM+xR47Tqoa65Hb16Lzj9jdRwhhLi47ZlwKg+jW2+vHcI2BQ6grusDRUXoL1dbHUUIIS7KXPMZ1IuCdh28dgxbFTjNWkJsU7mMIoTwaTonG3ZtRXXphQoI8NpxbFXgSilUtz5wcD/60LdWxxFCiDLpdemgTdR1Tq8ex1YFDqCu7Q5BtdBrllsdRQghfkebJa65323bo6IbefVY9ivw0Lqojl3QX2SgzxZaHUcIIX5t1zbIycbo1sfrhyp3T8wjR46Qmppa+jgrK4tBgwaRlJREamoqx44dIzo6mnHjxhEWFubVsD9T3fqiN65Cb1qL6urdP1GEEKIizNXLoG4EtO/k9WOVewYeExPDtGnTmDZtGi+88AK1atXimmuuIS0tjfj4eGbMmEF8fDxpaWleD1uqVVto1FguowghfIo+kQPbv3R9eRkY5PXjVegSyo4dO2jUqBHR0dFkZmaSlJQEQFJSEpmZmV4JWJbSLzO//Rr9w8FqO64QQlyMXpcOpunqp2pQ7iWUX1q3bh1du3YFIC8vj8hI191FkZGRnDx5sszXpKenk56eDkBKSgoOh6MqeUuZ/QdybP471M5cTfhV4zzyngCBgYEey2glfxkHyFh8kb+MAzw3Fm2aHN+wEuOKDtRvd6UHkpXP7QIvLi5m8+bNDBkypEIHcDqdOJ3/vU6dnZ1doddfjOrQmYLPF3O23yBUrdoeeU+Hw+HRjFbxl3GAjMUX+cs4wHNj0bu3Yv50BPOWIR7/2cTExJT5vNuXULZu3Urz5s2pV68eABEREeTmutbnzs3NJTw8vOopK0hd3xfyz6A3rav2YwshxC+ZGcsgrC7q6s7Vdky3C/yXl08AEhISyMjIACAjI4PExETPpytP6yugUSx69dLqP7YQQpynTxyHbRtRXZ2oIO9/efkztwr87NmzbN++nU6d/jstJjk5me3btzNmzBi2b99OcnKytzJekFIKlXSD68vM7w9U+/GFEAJArz3/5eX1fav1uG5dA69duzZvvvnmr56rW7cukyZN8kqoilCde6I/fge9eilq6ANWxxFC1DDaLEGvWQZtrkI1KPtatbfY7k7M31KhdVEJXdEbVqEL862OI4SoaXZscd15mXRjtR/a9gUOoJJuhLMFssysEKLamRlLICISrrqm2o/tFwVO3GXQuBk6Y6nsmSmEqDb6eBbs3Izq2hsVWKHbajzCLwq89MvMQ/+BA/usjiOEqCH0atdyHsqLu+5cjF8UOIDq1B1qB6NXLbE6ihCiBtDF51xfXsYnoBwNLcngPwUeHILq3B2duQZ9uuzb+oUQwlP0lg2uPS+797Msg98UOJz/MrP4HHrdCqujCCH8nF61GKIbQburLcvgXwXeuBm0aovOWII2TavjCCH8lD78HXyzG5V0I8qwrkb9qsABVPd+cOwo7N5qdRQhhJ/SGUsgMAjVtZelOfyvwDt0hroRmPJlphDCC3RBPnrDKlRiN1RY9S/i90v+V+CBQahufWF7pmuOphBCeJDeuArOFqB6WPfl5c/8rsABVFJfQMmUQiGER2mt0Z9/Ck1bQrNWVsfx0wKvHw1Xd0KvXY4uOmt1HCGEv/h6O/z4ParnTSilrE7jnwUOYPTsD6dPoTPXWB1FCOEnzJWfQlg4KrGb1VEAPy5wWl8BMU3QKz+R9VGEEFWmj2fBV1+iuvVBBdWyOg7gxwWulEL1uMm1Psq3X1sdRwhhcz9/p6YsWDb2Qvy2wAHUtd0hONT1pYMQQlSSLjqLXrscru6Eioq2Ok4p/y7wOsGorr3Qm9ehT+RYHUcIYVM6cw2cPoXR4yaro/yKWwvYnjlzhlmzZvH999+jlOKBBx4gJiaG1NRUjh07RnR0NOPGjSMsLMzbeStM9eiHXrEInbEUNWCI1XGEEDajtUavWAQxTeCyeKvj/IpbZ+BvvfUW7du35+WXX2batGnExsaSlpZGfHw8M2bMID4+nrS0NC9HrRzVIAau6OhaH+XcOavjCCHs5pvd8P0BVK/+PjF18JfKLfD8/Hz27NlDz549AQgMDCQ0NJTMzEySkpIASEpKIjMz07tJq8Bw3gyn8mRKoRCiwswViyC0LqpTD6uj/E65l1CysrIIDw/nlVde4eDBg8TFxTFixAjy8vKIjIwEIDIykpMny16DOz09nfT0dABSUlJwOBwejO8e3c3J8XlvoTKWUP/m2y76WzQwMNCSjJ7mL+MAGYsv8pdxwMXHUpL1I9nbNhIyYAh1Y2OrOVn5yi3wkpISDhw4wMiRI2nVqhVvvfVWhS6XOJ1OnE5n6ePs7OxKBa0qM+lG9Nz/I3vjGlSrthf8dw6Hw7KMnuQv4wAZiy/yl3HAxcdifjwHgMJOPThr4XhjYmLKfL7cSyhRUVFERUXRqpXrvv9rr72WAwcOEBERQW5uLgC5ubmEh1u7Kld5VOceEBKKuWKh1VGEEDagzxai1yyHq6/1qamDv1RugderV4+oqCiOHDkCwI4dO2jcuDEJCQlkZGQAkJGRQWJioneTVpGqXQfVrQ9s3Yg+fszqOEIIH6c3roL8Mxi9brE6ygW5NY1w5MiRzJgxg+LiYho0aMCDDz6I1prU1FRWrlyJw+Fg/Pjx3s5aZarHTejlC9Cff4oaOMLqOEIIH6VN0zV1sEkLaNnG6jgX5FaBN2vWjJSUlN89P2nSJI8H8iYV1QDVoTN69TJ0/9tRdYKtjiSE8EW7trpWHbxnnM9NHfwlv74Tsyyq9wAoOCMbHwshLshMXwD16qMSrrM6ykXVvAJvcTnEXYZesRBtllgdRwjhY/Th72D3NlSPm1CBQVbHuagaV+AARu8Bro2Pv/Ldm4+EENbQ6QuhVm1U0g1WRylXjSxwru4MUQ1cfyYJIcR5+mQu+osMVJeeqNC6VscpV40scBUQgOp5E+zbhT643+o4QggfoVctgeJzqF43Wx3FLTWywAHUdX2gTjB6uZyFCyHOr/m9aglcmYhq1NjqOG6puQUeEorq1ge9aY1rqyQhRI2mN3wOp/Iw+txqdRS31dgCB1Dn77DS6YssTiKEsJI2TfTyNGjaElq3szqO22p2gUdFoxK7odcsR+eftjqOEMIq27+ErCOovn/w6Rt3fqtGFziA6pMMZwvQq5dZHUUIYRFz2Xw4f6e2nUiBN2kBba5ybbtWLDv2CFHTFH29A/bvQfUegAoIsDpOhdT4Agcw+iTDiRz0l6utjiKEqGb5C9+FkFBUV2f5/9jHSIEDtOsAsU3Ry+ajTdPqNEKIaqKP/sDZjRmo7v1subidFDiglELd8Ac4coiizRusjiOEqCZ6+XwIDEL16m91lEqRAj9PJXSD+tGcmT/H6ihCiGqgT+SgN6wkuEc/VHik1XEqRQr8PBUYiOqTzLk9X6H377E6jhDCy/SKRVBiEpI8xOoolSYF/gvqut6ouhGYSz+yOooQwot0/hl0xhJUh84EXmKP2+bLIgX+C6p2HUL6/RG++hJ95JDVcYQQXqJXL4WCfNSNf7Q6SpVIgf9GSL+BUKs2eunHVkcRQniBPlfkWj6jzVWopi2tjlMlbu2JOXr0aOrUqYNhGAQEBJCSksLp06dJTU3l2LFjREdHM27cOMLCwryd1+uM8HquRa5WLUYPGIqKirY6khDCg/T6lZCXg3HPOKujVJlbBQ4wefJkwsPDSx+npaURHx9PcnIyaWlppKWlMWzYMK+ErG6qTzJ61RL08vmoO0ZZHUcI4SG6pAS97GNo3houv9LqOFVW6UsomZmZJCUlAZCUlERmpv9sT6bqR6M693AtcnXyhNVxhBAeojethWNHMfoNtNWiVRfi9hn41KlTAejduzdOp5O8vDwiI11zJyMjIzl58mSZr0tPTyc9PR2AlJQUHA5HVTN7VWBgIA6Hg+I77uX4+hXUWZ9O3WH3Wx2rwn4ehz+QsfgeO45DmyY5y+ejL21OVM9+KMN1/mrHsfzMrQKfMmUK9evXJy8vj+eee46YmBi3D+B0OnE6/7vGQHZ2dsVTViOHw+HKWCsY1aEL+Ys/pPD6G1Ah9rq+XzoOPyBj8T12HIfe9gXmof+g7hnH8Zyc0uftMJYLda5bl1Dq168PQEREBImJiezfv5+IiAhyc3MByM3N/dX1cX+hbhwIBfnozxdbHUUIUQVaa8zF81xLxiZeb3Ucjym3wAsLCykoKCj939u3b6dJkyYkJCSQkZEBQEZGBomJid5NagHVJA7iE9DpC9CFBVbHEUJU1p5tcGAf6oY/2m7J2Isp9xJKXl4eL730EgAlJSVcd911tG/fnhYtWpCamsrKlStxOByMHz/e62GtYNw0CDPlUXTGUlRf++yVJ4T4L/OT96FelC2XjL2Ycgu8YcOGTJs27XfP161bl0mTJnkllC9RLS53bfiwfD66Rz9UrdpWRxJCVIDeuxO+2Y0aPAoVFGR1HI+SOzHdYPS/HU6eQK9ZbnUUIUQFmZ++D+H1UN16Wx3F46TA3aBaXwGt26GXfoQ+V2R1HCGEm/S3X8Oer1B9b/XLv56lwN1k3HS7a9u1delWRxFCuMn85H0IC0cl3Wh1FK+QAndXm6ugxeXoJR+iz8nmx0L4On3gG9i52bVZce06VsfxCilwNymlMPoPhpxsOQsXwgbMRe9CaF1Uz5usjuI1UuAV0e5qiLsMvWSenIUL4cP0gX2wYxOqTzKqTojVcbxGCrwClFIYtwyRs3AhfJy58F0I8++zb5ACr7i27V3XwhfLWbgQvkj/Z6/r2nefW/367BukwCtMKYVx8x2Qm41e95nVcYQQv2EuOn/23aOf1VG8Tgq8Mn4+C/90nswLF8KH6G+/hp1basTZN0iBV4pSCmPAUDhxHL16mdVxhBDnmQvmQt0IVM/+VkepFlLglaTaXAWXxbuuhZ8ttDqOEDWe/nq7667LGwf67bzv35ICrwJjwFDXGimff2p1FCFqNK216+y7Xn1U0g1Wx6k2UuBVoFq1hSs6oJd+jC7ItzqOEDXXri2wfw/qpkF+uebJhUiBV5ExYCicOYVesdDqKELUSFprzLS5rt12rvO/FQcvRgq8ilSzVtD+WvTyNPTpsjd2FkJ40daNcHA/qv/tqED/Wu+7PFLgHmAkD4XCAvTSj6yOIkSNos0SzLQ50KgxqnNPq+NUOylwD1CxTVHXdkev/BSde9zqOELUGHrDKvjxe4zkYX6116W73C5w0zR59NFHSUlJAeD06dNMmTKFMWPGMGXKFE6fPu21kHagbhkCpon+5D2rowhRI+hz59AL/w1NW0KHzlbHsYTbBb548WJiY2NLH6elpREfH8+MGTOIj48nLS3NG/lsQzkaoq7vi177GfroD1bHEcLv6YwlkHMM4w/DUUpZHccSbhX48ePH2bJlC7169Sp9LjMzk6SkJACSkpLIzMz0TkIbUf0HQVAt9IK5VkcRwq/pwnz04nlw+ZXQpr3VcSxT7q70ALNnz2bYsGEUFBSUPpeXl0dkZCQAkZGRnDxZ9gyM9PR00tNdS6+mpKTgcDiqmtmrAgMDK5/R4eD0LYM5M282EbffTVDLNp4NVwFVGoePkbH4HqvHcfq91zlzKo/6dz9EUHR0ld7L6rFURbkFvnnzZiIiIoiLi2PXrl0VPoDT6cTpdJY+zs7OrvB7VCeHw1GljLrbDbDkY3Jefxnjkecs+9OuquPwJTIW32PlOHReLmbav6FjF/LqN4Qq5rDDZxITE1Pm8+UW+N69e9m0aRNbt26lqKiIgoICZsyYQUREBLm5uURGRpKbm0t4eLjHQ9uRCg5B9R+Mfu+fsHMLxHe0OpIQfkUveheKz2HceqfVUSxX7jXwIUOGMGvWLGbOnMnYsWO54oorGDNmDAkJCWRkZACQkZFBYmKi18PahUrqC9GNMD+ajTZLrI4jhN/QRw+j1yxHXd8X1bDss9KapNLzwJOTk9m+fTtjxoxh+/btJCcnezCWvanAINStd8IPB9EbV1kdRwi/Yc5/B4Jqo/oPtjqKT3DrS8yftWvXjnbt2gFQt25dJk2a5JVQ/kAldEUvb4VOm4tOuK5GLbAjhDfo/XtgywbULUNQ4fWsjuMT5E5ML1FKYdx2t2vrtc8WWB1HCFvTWmPOexMi6qN6D7A6js+QAvci1foK10JXSz5C5+VaHUcI29Kb1sJ/9qKSh6LqBFsdx2dIgXuZMXAEFBe5bvkVQlSYPleE/uhf0Lg5qkvNW7DqYqTAvUw1jEF174de8xn6h4NWxxHCdvSKRXA8C2PQSJRR8xasuhgp8Gqgbh4MwcGua3hCCLfpU3muW+avTHTtQyt+RQq8GqjQuq5pT7u2ondssjqOELahF8yFs4UYA++2OopPkgKvJqpHP2gYi/nBG+jic1bHEcLn6cMH0KuXo3rchLqksdVxfJIUeDVRgUEYg0bC0R/Qny+2Oo4QPk1rjfne6xAairr5Dqvj+Cwp8OoUn+DaxX7Re+hTeVanEcJ3bd0Ae3egBgxFhYZZncZnSYFXI6UUxqB7oKgQnSZrhgtRFn2uCPODNyG2KapbX6vj+DQp8GqmLrkU1eMm9Jpl6EPfWh1HCJ+jl6e5pg3efm+N3OeyIqTALaBuHgxh4Zj/fhWttdVxhPAZ+vgx9OIPoEMXmTboBilwC6iQMNQf74Jvv0Zv+NzqOEL4DPODNwBclxpFuaTALaI694TmrdEfzUbnn7E6jhCW07u3wpb1qH6DUFFV2yatppACt4gyDIyh98OpPPSi96yOI4SldPE5zHdfg+hGqD7JVsexDSlwC6mmLVHd+qJXLkIf/s7qOEJYRq9YBEcPY9wxChVUy+o4tiEFbjF16zAICcOcOwttmlbHEaLa6ePH0AvfhauuQcUnWB3HVqTALabCwlEDR8D+3ej1K6yOI0S1M9/7JwDGHaMsTmI/UuA+QHXuCS3bur7QPHXS6jhCVBv91Zew7QvUzYNRUQ2sjmM75e6JWVRUxOTJkykuLqakpIRrr72WQYMGcfr0aVJTUzl27BjR0dGMGzeOsDC55bUylGFgDHsAc8pY9EezUSPGWB1JCK/TZwsx3/0nxDRBOWWbtMoo9ww8KCiIyZMnM23aNF588UW2bdvGvn37SEtLIz4+nhkzZhAfH09aWlo1xPVfKrYpyjkAvS4dvW+X1XGE8Dr9yfuuOy6HPoAKrND+6uK8cgtcKUWdOnUAKCkpoaSkBKUUmZmZJCUlAZCUlERmZqZ3k9YA6ubBENUA852Z6HOy5KzwX/r7A+jl81FdnajW7ayOY1tu/dozTZPHHnuMo0eP0rdvX1q1akVeXh6RkZEAREZGcvJk2ddu09PTSU9PByAlJQWHw+Gh6N4RGBhoacazD07kxJTxBK9eTNjtlb8bzepxeJKMxfdUZRy6pIScFx6DuhFE3fdnjLrhHk5XMXb+TNwqcMMwmDZtGmfOnOGll17i0KFDbh/A6XTidDpLH2dnZ1c8ZTVyOBzWZmzSEnVNEmc+/BcFbTugLrm0Um9j+Tg8SMbie6oyDjN9IXr/HtSf/kzO2SI4a+3Pww6fSUxMTJnPV2gWSmhoKG3btmXbtm1ERESQm5sLQG5uLuHh1v4W9Sfq9nugdjDm2zNlbrjwK/p4FjptDsQnoBK7WR3H9sot8JMnT3LmjGutjqKiInbs2EFsbCwJCQlkZGQAkJGRQWJioneT1iAqvB7qtpGuueGrl1odRwiP0Fpjzvk/AIyh96OUsjiR/ZV7CSU3N5eZM2dimiZaazp37kzHjh1p3bo1qamprFy5EofDwfjx46sjb42huvREf7EK/eG/0PGJsriPsD294XPYuRl1+70y59tDlK7mBamPHDlSnYerMF+6Hqazf8J8+iFo2Qbj4acrdMbiS+OoKhmL76noOPSJHMzJoyGmCcaE51GG79xDaIfPxCPXwEX1Uo6GrnXDd22V2+yFbWmtMefOgqIijLse8qnytjv5Sfo4lXQjtGqLfv8NdO5xq+MIUWF60zrYthE1YAiqUWOr4/gVKXAfpwwD464xUHzOdYOPbMEmbESfzEX/+/+gWStU72Sr4/gdKXAbUA1jUH+4E3ZsQq/9zOo4QrhFa4359kwoLMS4+2HZoNgLpMBtQvXsD62vcF1Kyf7J6jhClEtv+By++hJ16zBUTBOr4/glKXCbUIaBcX6VQnP2DLnBR/g0nXMM/d4/oVVblPMWq+P4LSlwG1HRjVx3ae7dgf78U6vjCFEmbZqYs2eAaWKMeBhlyKUTb5ECtxl1XW+IT0B/OBv9g/tr0ghRXfTKRbDnK9RtI1ENLrE6jl+TArcZpRTGiIcgOATz9emy7KzwKfrwd+iP3nbtb3l9X6vj+D0pcBtS4ZGuqYWHD7gWBhLCB+hzRZivT4fgEIw7/0fWOqkGUuA2pa5KRCXdgP4sDb3nK6vjCIH++B344aBrymB4Pavj1AhS4DambhsJDWMw30yVzZCFpfSOzej0Baju/VDxCVbHqTGkwG1M1a6D8acJcPok5lsvy12awhL6RA7mWy9DbFPUbXdbHadGkQK3OdUkDjVwpOsuzRULrY4jahhtmphvpsLZAoz7HkXVqm11pBpFCtwPqJ43wVXXuNYOP7jf6jiiBtFLP3JNGRw8qtLb/4nKkwL3A66phWMgvB7mqy+i889YHUnUAHrfLvSCuajEbq77E0S1kwL3EyosHGPUn+F4Fubsv8v1cOFVJSdyMP85DRyNUMNHy5RBi0iB+xHVsq1rA4itG8lf9L7VcYSf0mYJJ/82GfJPYzzwGCo4xOpINVa5e2JmZ2czc+ZMTpw4gVIKp9NJv379OH36NKmpqRw7dozo6GjGjRtHWFhYdWQWF6F6J6O/2cPpt2diNIhFtWxjdSThZ/Si9ynasRl110Ooxs2tjlOjlXsGHhAQwPDhw0lNTWXq1KksW7aMw4cPk5aWRnx8PDNmzCA+Pp60tLRqiCvKo5TCuHsMAY6GmK++gD6Za3Uk4Uf0jk3oT9+nTs9+GHLd23LlFnhkZCRxcXEABAcHExsbS05ODpmZmSQlJQGQlJREZmamd5MKt6mQMCIeex7yT2POegFdXGx1JOEHdNYRzNemQ+NmhI/6s9VxBG5cQvmlrKwsDhw4QMuWLcnLyyMyMhJwlfzJk2XfCZienk56ejoAKSkpOByOKkb2rsDAQJ/P6I7ARo0IH/0EJ1OfpvYn7xJ+7zirI1Wav3wmYN+xmAX55Mx6ARUQQNRTLxEUGoajdh2rY3mEXT8TqECBFxYWMn36dEaMGEFIiPtfWjidTpxOZ+nj7OzsiiWsZg6Hw+czusPhcHCmbQeU8xYKPp1HYcPGGJ17WB2rUvzlMwF7jkVrjf7nNPThgxhjnybXCMJRXGy7cVyIHT6TmJiYMp93axZKcXEx06dPp1u3bnTq1AmAiIgIcnNd11dzc3MJDw/3UFThSeqPI+CyePTb/0Af2Gd1HGFDevE89Ka1qD8MR7Vtb3Uc8QvlFrjWmlmzZhEbG0v//v1Ln09ISCAjIwOAjIwMEhMTvZdSVJoKDMS47zGIiMScORWd49tnGsK36C0b0GlzUJ2SUH3/YHUc8RvlFvjevXtZvXo1O3fuZMKECUyYMIEtW7aQnJzM9u3bGTNmDNu3byc5Obka4orKUHXDMR76CxQWYr7yV/TZs1ZHEjagvz/gWuekeWvXlEG5WcfnKF3Nt+wdOXKkOg9XYXa4HuaOssahv8rEnPkcqmNX1J/+jDLscR+Xv3wmYJ+x6JO5mFP/7NrX8snpqHr1f/Xf7TIOd9hhLFW6Bi78g7oqEfXHu9Cb1spOPuKC9NmzmP/7HJw+ifE/T/6uvIXvqNA0QmF/qs+tkHUUveRDzOhGGN36WB1J+BBtlri2RTu4H+PBJ1BNW1odSVyEFHgNo5SCIfehc7LQc15B149Gtbva6ljCR+h5b8G2jajBf0K172R1HFEOuYRSA6mAAIz7HoWYppizUtCHvrU6kvABZvoCdPpCVK+bMXrdbHUc4QYp8BpK1QnBGDMJQsIw//4MOutHqyMJC5kbV6HffwM6dEENGml1HOEmKfAaTEVGYYx9BswSzJcny8JXNZTeuQU9++9wWTzGveNRRoDVkYSbpMBrOHVJY4yHJkFerutMXHbzqVH0f/ZizkqBS5q4vrQMqmV1JFEBUuACFXcZxgMT4YdDmDOeQZ8ttDqSqAb6+wOYf38awuthPDwZFRJqdSRRQVLgAgB1RUeMP/0Z/rPPdcv9uSKrIwkv0j8exkydBLWDMcZPkbneNiUFLkqpjl1QI8bAnq9cmyMXn7M6kvACfewo5t/+Akq5ytvR0OpIopKkwMWvGF16oobeD199ifnqNClxP6OPHcV86Uk4V4Qx7llUo1irI4kqkAIXv2N074e6YxRs2yhn4n7EVd5PwNlCjPHPoho3szqSqCIpcFEmo2d/1JD7YNsXrhI/JyVuZ/8t77OuyyZNWlgdSXiAFLi4IKPHTagh97tKfOZzsgytTekjhzBfmPiL8o6zOpLwEClwcVFGj36oux6C3dsw/z4ZXZBvdSRRAfrgt5jTHgc0xoS/Snn7GSlwUS7jut6oP/0Z/rMXc/pT6FNlb2AtfIvetwtz+pOuqYKPPo+KbWp1JOFhUuDCLUZiN4wHn4AjhzBTHkUfO2p1JHERest61zzviEhXeTcoe0MAYW9S4MJt6spEjPHPwumTrhI/KKsY+iLz808xZ70ATeIwHn0BVT/a6kjCS6TARYWolm0xJr4AgUGY055A79hkdSRxnjZNzA9no//9KlyZiDH+OVTdcKtjCS8qd0OHV155hS1bthAREcH06dMBOH36NKmpqRw7dozo6GjGjRtHWFiY18MK36AuuRTj8RcxZzyL+b/PoQbdjep1i2x6ayF9thDz9b+5NmNIugF1x32oAFlV0N+VewbevXt3nnjiiV89l5aWRnx8PDNmzCA+Pp60tDRv5RM+StWLwnjsBbjqGvT7b7h29ykutjpWjaRzsjFfnAhffYm6/V7U0AekvGuIcgu8bdu2vzu7zszMJCkpCYCkpCQyMzO9k074NFW7DsYDE1E3/hG9ehnm355C58ma4tVJ79uJ+dw4yPoR46G/YDjlL6GapFJ7Yubl5REZGQlAZGQkJ09eeFpZeno66enpAKSkpOBwOCpzyGoTGBjo8xndUa3jGPUIBZddwclXUuCvjxA+YSq1Lo/32Nv7y2cCnhuL1pr8Tz7g9Ox/ENAolnoTnyfw0uYeSOge+Ux8g9c3NXY6nTidztLH2dnZ3j5klTgcDp/P6I5qH0e7jhgTX8T8v+fJfWo06ra7UT37e+Rs0F8+E/DMWHRBPvqdmejMNdC+E/rusZwIDoVq/BnJZ1K9YmLKngZaqVkoERER5Oa6/lTOzc0lPFy+6RagLm2O8eTfoG179HuvYb7yV/SZU1bH8iv6u28wp4xFb1qHSh6G8cDjshFDDVapAk9ISCAjIwOAjIwMEhMTPRpK2JcKDcN46C+oQffAjs2Yzz6M3rfT6li2p00Tc3kaZspjUFKMMeGvGDcNQhkyE7gmK/cSyssvv8zu3bs5deoU999/P4MGDSI5OZnU1FRWrlyJw+Fg/Pjx1ZFV2IRSCtV7ALp1O8x/TsN86UmU8xbUrcNlz8VK0Nk/Yb71d9i3E9pfizHiIVRoXatjCR+gtNa6Og945MiR6jxchdnhepg7fGUcurAA/dFs9KolcMmlGHePRTVvVaH38JWxeEJFxqK1Rq9Zjv7gTVC4pgh2dfrELJOa+plY5ULXwL3+Jaao2VSdYNTQB9Dtr8WcPQPz+Qmonjehkoei6oRYHc9n6aOHMd95xXXWfVk8xt0Po6IaWB1L+BgpcFEtVLurMZ75B3r+O+iVn6C3bMC4YxS07+QTZ5S+QhedRS/9GL1kHtSqjbrzf1xn3XKtW5RBClxUGxUSihp6P/ra7pjvzMR85a/Q5iqM2++t8Uudaq1h8zrMD2fD8SxUYjfXJZOISKujCR8mBS6qnWpxOcZTqeiMpeiFczGfeRh1fR9U/8GoevWtjlft9LdfY378L9i3Cxo3w3jkOdTlV1odS9iA/F0mLKECAzF69ceY+iqq+w3otZ9hPjkK86N/2X7uuGmaPProo7Rr147Y2FjWr19f5r/TPxyk5B/PYaY8Cj8eRg19AOOpVI+X95133snYsWM9+p4AY8eO5c477/T6a8SFyRm4sJQKC0cNuR/tvAW98F30so/Rqxajkm5E9R5gy0sIK1as4IMPPmDevHk0bdqUevXq/eq/6wPfYC6ZB9u+gDohqORhqF43o+oEWxO4kp599lkqOomtMq8RFyYFLnyCahCDuvcR9A1/RC+eh16ehl6xCNWlJ8V/GAahEVZHdNt3331HgwYNfnWDmzZLYMcWclcvwdy+CUJCUf1uc/2Ssumc7srcgS13bXuWFLjwKapxM9SoCejkoehl89HrV3J89TJofQWqez9U+2uq7WYgrTX79u3jsssuc/s1Y8eOZd68eQDExsbSOCaGDc89TmH6Ip5fv42FR09wqtikbbt2TBrammvOl/fAgQO57LLLmDp16q/eKycnh7fffrv037Rq1Yrw8HDmzp2LYRgMHDiQp556CuP8LJWCggIef/xxPv30U0JCQrjnnnvKzTxw4EBatmxJcHAwH3zwAYZh8PDDDzN8+HCeeeYZ5s+fT1hYGI899hgDBw6sUr6yXlPRY7v7s3L3fe+//363P19fI9fAhU9SDWIwho/GePEtwu58EI5nof/5Iuaf78J85xX0/t1o0/Ta8bXWPPbYY4wcOZJz5865/bpnn32WsQ89xCVRUWwaPZxF7RqgP5zNX3cf4pO8Il6b+y7L0tNp064dQ4cO5aeffqpQrvnz5xMYGMiCBQt47rnneP3111m4cOGvjr9mzRpee+013n//fXbu3MkXX3zh1vuGhYWxaNEiRo8ezeTJk7nnnnuIi4tj8eLF3HbbbUyYMIGjRy++F2p5+bx57Mq+r6/fXHgxUuDCp6m64YTeOgzjr69iPPw0Kj4BvXEl5gsTMR8diTnnFfTOzeiisx475s/lvWLFCt5++22CgoLKf83JXMz1Kwh9538JXbmQgPxTNMj9CcfNgyicOI05u77liclP069/f1q1akVKSgrR0dHMnj27QtlatWrFhAkTaNGiBbfccgtdunRh7dq1AJw5c4b33nuPJ598ku7du3P55ZeTmppaevZ7Ma1bt+aRRx4hLi6O++67j/r16xMYGMi9995L8+bNGTduHFprNm26+BZ6F8vn7WNX9n03btxYoff1JXIJRdiCMgLgig6oKzqgC/PR275Eb9uI3rgKnbEUAoOgxeWoNleh4i6Dpi0rvUrfJ598wty5cwG4/vrrL/jvGjeIZsPE/0Hv+Qp+OOh6MtIBzVpBdgHGC6+jjAAO7t7NuXPnfnVNPCAggI4dO/LNN99UKFubNm1+9bhhw4alt4F/9913FBUV0bFjx9L/HhoayuWXX16h91VK4XA4fvW6oKAgIiIiyr3l/GL5vH3syr5vVlZWhd7Xl0iBC9tRdUJQ13aHa7u7zrz37kR//RV6z1fotDloAKWgYSzENEFd0hgaxaLqN4DIKKhX/6LX0Z1OJ926dWPfvr3MmvoskUqjc45D9lF09k9w9DCcOkmgUq5fHi3boDolodq2hyYtMF59Fb762vVLB0pnXZR1x+nPzymlfjc7o7iMLep++9eAUgrz/KWkqszuKOt9y3quvGNcLJ+njl2Vn1Vl8vkyKXBha6pWbYjviIp3nXXqM6fgu/3oA/vQ330Dh79Db9sIpsmv/i9fqzbUCYY6IfDL/SPPFVGrMJ83wgsZWVzAX8aPY0GXNgQZyvWaRo1RnbpC81ao5q2hcXNUOZdYmjdvTq1atfjyyy9Lz45LSkrYvHkzycnJAERFRf3uTHD37t00btzY7Z9F8+bNCQoKYsuWLTRt6rqzNT8/n71795Y+9gee+Fn5Cylw4VdUaF1odzWq3dWlz+lz5yD7KORko08ch9zjUHAGCgugIN81xe/n1wfVgjohBNcJ5s2b67Du0BFq9+4NUQ0g0lGpNUlCQkIYPnw4zz//PM2aNaNevXq89tprHDt2jLvuuguArl278vTTT7N8+XLi4uKYM2cOR44cqVAphYaGMnjwYKZOnUpUVBQNGzYkNTWVkpKS8l9sI574WfkLKXDh91RQEFxyKVxyKRVZNisU6OOhDE8++SQAo0aN4sSJE7Rr1465c+fSsGFDAAYPHsyePXtK19a/6667uOGGG8jJyanQcSZNmkR+fj733HMPwcHB3H333eTn53toFL7BUz8rfyDrgf+GHdYGdoe/jANkLL7IX8YB9hiLR/fEFEIIYT0pcCGEsKkqXQPftm0bb731FqZp0qtXr9Jv1IUQQnhfpc/ATdPkjTfe4IknniA1NZV169Zx+PBhT2YTQghxEZUu8P3799OoUSMaNmxIYGAgXbp0ITMz05PZhBBCXESlL6Hk5OQQFRVV+jgqKqrM24LT09NJT08HICUlBYfDUdlDVovAwECfz+gOfxkHyFh8kb+MA+w9lkoXeFmzD8u6VdjpdOJ0Oksf+/p0HTtMKXKHv4wDZCy+yF/GAfYYy4WmEVa6wKOiojh+/Hjp4+PHjxMZWf7uKRcK4kvskNEd/jIOkLH4In8ZB9h3LJW+Bt6iRQt+/PFHsrKyKC4uZv369SQkJHgymyUmTpxodQSP8JdxgIzFF/nLOMDeY6n0GXhAQAAjR45k6tSpmKZJjx49uPTSSz2ZTQghxEVUaR54hw4d6NChg6eyCCGEqAC5E/M3fvmFq535yzhAxuKL/GUcYO+xVPtiVkIIITxDzsCFEMKmpMCFEMKmZEOHMixZsoSlS5cSEBBAhw4dGDZsmNWRqmThwoXMmTOH119/nfDwcKvjVMo777zD5s2bCQwMpGHDhjz44IOEhlZu02Ir+MvCb9nZ2cycOZMTJ06glMLpdNKvXz+rY1WaaZpMnDiR+vXr23I6oRT4b+zcuZNNmzbx0ksvERQURF5entWRqiQ7O5sdO3bY9lbhn1155ZUMGTKEgIAA5syZw/z5823zi/Xnhd+eeuopoqKiePzxx0lISLDlFmABAQEMHz6cuLg4CgoKmDhxIldeeaUtxwKwePFiYmNjKSgosDpKpcgllN9Yvnw5AwYMKN29OiIiwuJEVfOvf/2LoUOHlrnMgZ1cddVVBJzffLh169a22j7LnxZ+i4yMJC4uDoDg4GBiY2Nt9Vn80vHjx9myZQu9evWyOkqlyRn4b/z44498/fXXvPfeewQFBTF8+HBatmxpdaxK2bRpE/Xr16dZs2ZWR/GolStX0qVLF6tjuM3dhd/sJisriwMHDtj2/x+zZ89m2LBhtj37hhpa4FOmTOHEiRO/e37w4MGYpsnp06eZOnUq3377LampqfzjH//w2TPYi41l/vz5PPXUU9UfqpIuNpbExEQAPv74YwICAujWrVs1p6s8dxd+s5PCwkKmT5/OiBEjCAkJsTpOhW3evJmIiAji4uLYtWuX1XEqrUYW+F/+8pcL/rfly5fTqVMnlFK0bNkSwzA4deqUz375d6GxHDp0iKysLCZMmAC4/lx87LHHeP7556lXr141JnTfxT4XgFWrVrF582YmTZpkqwKs7MJvvqq4uJjp06fTrVs3OnXqZHWcStm7dy+bNm1i69atFBUVUVBQwIwZMxgzZozV0SqkRhb4xSQmJrJz507atWvHkSNHKC4upm7dulbHqrAmTZrw+uuvlz4ePXo0zz//vM/+IirPtm3bWLBgAc888wy1a9e2Ok6F/HLht/r167N+/XrbFcXPtNbMmjWL2NhY+vfvb3WcShsyZAhDhgwBYNeuXSxatMiWn4kU+G/07NmTV155hUceeYTAwEBGjx5tq7M9f/XGG29QXFzMlClTAGjVqhWjRo2yOJV7/Gnht71797J69WqaNGlS+tfdHXfcIWsiWURupRdCCJuSaYRCCGFTUuBCCGFTUuBCCGFTUuBCCGFTUuBCCGFTUuBCCGFTUuBCCGFT/w9Q3BUBxtOWlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "w = np.linspace(-7,5,100)\n",
    "l = 2*w**2+4*w+5\n",
    "plt.plot(w,l)\n",
    "plt.text(-1,2.5,'$\\leftarrow$found minimum',fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "\n",
    "# Newton's Update Method\n",
    "<img src=\"PDF_slides/newton.png\"  width=\"600\">\n",
    "\n",
    "But how do we translate this over to objective funtions with more than one variable? We need a second derivative of a multivariate equation... enter, the hessian. Our new update is defined by Newton's method:\n",
    "\n",
    "$$ w \\leftarrow w - \\underbrace{[\\frac{\\partial^2}{\\partial w}l(w)]^{-1}}_{\\text{inverse 2nd deriv}}\\underbrace{\\frac{\\partial}{\\partial w}l(w)}_{\\text{derivative}}  $$\n",
    "\n",
    "such that, in multiple dimensions we can approximate the update as:\n",
    "\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\cdot \\underbrace{\\mathbf{H}(\\mathbf{w})^{-1}}_{\\text{inverse Hessian}}\\cdot\\underbrace{\\nabla l(\\mathbf{w})}_{\\text{gradient}}$$\n",
    "\n",
    "where the Hessian is defined as follows for any multivariate equation $l(\\mathbf{w})$:\n",
    "$$ \\nabla^2 l(\\mathbf{w}) = \\mathbf{H}(\\mathbf{w})   $$\n",
    "\n",
    "$$  \\mathbf{H}(\\mathbf{w}) =  \\begin{bmatrix}\n",
    "        \\frac{\\partial^2}{\\partial w_1}l(\\mathbf{w}) &  \\frac{\\partial}{\\partial w_1}\\frac{\\partial}{\\partial w_2}l(\\mathbf{w}) & \\ldots     & \\frac{\\partial}{\\partial w_1}\\frac{\\partial}{\\partial w_N}l(\\mathbf{w})  \\\\\n",
    "        \\frac{\\partial}{\\partial w_2}\\frac{\\partial}{\\partial w_1}l(\\mathbf{w})  & \\frac{\\partial^2}{\\partial w_2}l(\\mathbf{w}) &  \\ldots     & \\frac{\\partial}{\\partial w_2}\\frac{\\partial}{\\partial w_N}l(\\mathbf{w})  \\\\\n",
    "        &  \\vdots &\\\\\n",
    "        \\frac{\\partial}{\\partial w_N}\\frac{\\partial}{\\partial w_1}l(\\mathbf{w})  & \\frac{\\partial}{\\partial w_N}\\frac{\\partial}{\\partial w_2}l(\\mathbf{w}) &  \\ldots     & \\frac{\\partial^2}{\\partial w_N}l(\\mathbf{w}) \\\\\n",
    "     \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "____\n",
    "\n",
    "<img src=\"PDF_slides/spider_hessian.png\"  width=\"400\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "For logistic regression, we can calculate the formula for the $j^{th}$ and $k^{th}$ element of the Hessian as follows:\n",
    "\n",
    "$$ \\mathbf{H}_{j,k}(\\mathbf{w}) = \\frac{\\partial}{\\partial w_k} \\underbrace{\\frac{\\partial}{\\partial w_j}l(\\mathbf{w})}_{\\text{first derivative}} $$\n",
    "\n",
    "But we already know the result of the $j^{th}$ partial derivative from our calculation of $\\nabla l(\\mathbf{w})$: \n",
    "\n",
    "$$ \\frac{\\partial}{\\partial w_j}l(\\mathbf{w}) = \\sum_i \\left(y^{(i)}-g(\\mathbf{w}^T\\cdot\\mathbf{x}^{(i)})\\right)x_j^{(i)} $$\n",
    "\n",
    "So we can plug this back into the equation to get:\n",
    "\n",
    "$$ \n",
    "\\begin{split}\n",
    "\\mathbf{H}_{j,k}(\\mathbf{w}) & = \\frac{\\partial}{\\partial w_k}\\sum_i \\left(y^{(i)}-g(\\mathbf{w}^T\\cdot\\mathbf{x}^{(i)})\\right)x_j^{(i)} \\\\\n",
    " & = \\underbrace{\\sum_i \\frac{\\partial}{\\partial w_k} y^{(i)}x_j^{(i)}}_{\\text{no dependence on }k\\text{, zero}} -\\sum_i \\frac{\\partial}{\\partial w_k}g(\\mathbf{w}^T\\cdot\\mathbf{x}^{(i)})x_j^{(i)} \\\\\n",
    " & = -\\sum_i x_j^{(i)}\\underbrace{\\frac{\\partial}{\\partial w_k}g(\\mathbf{w}^T\\cdot\\mathbf{x}^{(i)})}_{\\text{already know this as }g(1-g)x_k} \\\\\n",
    " & =  -\\sum_{i=1}^M \\left[g(\\mathbf{w}^T\\mathbf{x}^{(i)})[1-g(\\mathbf{w}^T\\mathbf{x}^{(i)})]\\right]\\cdot{x_k}^{(i)}{x_j}^{(i)} \\\\\n",
    "\\end{split}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "Therefore the Hessian for logistic regression becomes (adding in the regularization term also):\n",
    "$$ \\mathbf{H}_{j,k}(\\mathbf{w}) =\\left( -\\sum_{i=1}^M \\underbrace{\\left[g(\\mathbf{w}^T\\mathbf{x}^{(i)})[1-g(\\mathbf{w}^T\\mathbf{x}^{(i)})]\\right]}_{\\text{scalar value for each instance}}\\cdot\\underbrace{{x_k}^{(i)}{x_j}^{(i)}}_{i^{th}\\text{ instance elements}} \\right) + \\underbrace{2\\cdot C}_{\\text{regularization}}  $$\n",
    "\n",
    "\n",
    "This equation can be calcuated in a for loop, for each $j,k$ element in the Hessian and for each instance in the dataset, but this would be **slow** in python. To vectorize this operation, we need to have each operation be linear algebra, so that it can be run efficiently with numpy. \n",
    "____\n",
    "### Calculating the Hessian for Logistic Regression using Linear Algebra\n",
    "First notice that the sum of each terms ${x_k}^{(i)}{x_j}^{(i)}$ that forms a matrix can be calculated as follows:\n",
    "\n",
    "$$    \n",
    "\\begin{bmatrix}\n",
    "        \\sum_{i=1}^M {x_1}^{(i)}{x_1}^{(i)} & \\sum_{i=1}^M {x_1}^{(i)}{x_2}^{(i)} & \\ldots     &  \\sum_{i=1}^M {x_1}^{(i)}{x_N}^{(i)} \\\\\n",
    "        \\sum_{i=1}^M {x_2}^{(i)}{x_1}^{(i)} & \\sum_{i=1}^M {x_2}^{(i)}{x_2}^{(i)} & \\ldots     &  \\sum_{i=1}^M {x_2}^{(i)}{x_N}^{(i)} \\\\\n",
    "        &  \\vdots & \\\\\n",
    "        \\sum_{i=1}^M {x_N}^{(i)}{x_1}^{(i)} & \\sum_{i=1}^M {x_N}^{(i)}{x_2}^{(i)} & \\ldots     &  \\sum_{i=1}^M {x_N}^{(i)}{x_N}^{(i)} \\\\ \\\\\n",
    "\\end{bmatrix}   \n",
    "%\n",
    "= \\mathbf{X}^T \\cdot\\mathbf{X} = \\mathbf{X}^T \\cdot\\mathbf{I} \\cdot\\mathbf{X}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{I}$ is the identity matrix of size $M\\text{x}M$. This can be seen in the following exploded view of the matrix operations: \n",
    "\n",
    "$$ \\mathbf{X}^T \\cdot\\mathbf{I} \\cdot\\mathbf{X}=\n",
    "\\begin{bmatrix}\n",
    "        \\uparrow &    \\uparrow    &    & \\uparrow  \\\\\n",
    "        \\mathbf{x}^{(1)} &  \\mathbf{x}^{(2)} & \\ldots   & \\mathbf{x}^{(M)}  \\\\\n",
    "        \\downarrow &    \\downarrow    &   & \\downarrow  \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "\\begin{bmatrix}\n",
    "        1 &  0 & \\ldots & 0 \\\\\n",
    "        0 &  1 & \\ldots & 0 \\\\\n",
    "          &  \\vdots  & &    \\\\\n",
    "        0 &  0 & \\ldots & 1 \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "\\begin{bmatrix}\n",
    "        \\leftarrow &  \\mathbf{x}^{(1)}      & \\rightarrow  \\\\\n",
    "        \\leftarrow &  \\mathbf{x}^{(2)}      & \\rightarrow  \\\\\n",
    "        &  \\vdots &\\\\\n",
    "        \\leftarrow &  \\mathbf{x}^{(M)}      & \\rightarrow  \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "$$\n",
    "\n",
    "With this equation in mind, we can revisit the calcualtion of the Hessian and use matrix operations to define the needed multiplication in an exploded view of the operations:\n",
    "\n",
    "$$ \\mathbf{H}[l(\\mathbf{w})]=\n",
    "\\begin{bmatrix}\n",
    "        \\uparrow &    \\uparrow    &    & \\uparrow  \\\\\n",
    "        \\mathbf{x}^{(1)} &  \\mathbf{x}^{(2)} & \\ldots   & \\mathbf{x}^{(M)}  \\\\\n",
    "        \\downarrow &    \\downarrow    &   & \\downarrow  \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "\\begin{bmatrix}\n",
    "        g(\\mathbf{w}^T\\mathbf{x}^{(1)})[1-g(\\mathbf{w}^T\\mathbf{x}^{(1)})] & \\ldots & 0 \\\\\n",
    "          &  \\vdots   &    \\\\\n",
    "        0  & \\ldots & g(\\mathbf{w}^T\\mathbf{x}^{(M)})[1-g(\\mathbf{w}^T\\mathbf{x}^{(M)})] \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "\\begin{bmatrix}\n",
    "        \\leftarrow &  \\mathbf{x}^{(1)}      & \\rightarrow  \\\\\n",
    "        \\leftarrow &  \\mathbf{x}^{(2)}      & \\rightarrow  \\\\\n",
    "        &  \\vdots &\\\\\n",
    "        \\leftarrow &  \\mathbf{x}^{(M)}      & \\rightarrow  \\\\\n",
    "\\end{bmatrix}\n",
    "%\n",
    "$$\n",
    "\n",
    "___\n",
    "Or, more succinctly as follows (adding in the regularization term as well):\n",
    "\n",
    "$$ \\mathbf{H}[l(\\mathbf{w})] =  \\mathbf{X}^T \\cdot \\text{diag}\\left[g(\\mathbf{X}\\cdot\\mathbf{w})\\odot(1-g(\\mathbf{X}\\cdot\\mathbf{w}))\\right]\\cdot \\mathbf{X} -2C$$\n",
    "\n",
    "___\n",
    "\n",
    "Now we can place the Hessian derivation into the Newton Update Equation, like this:\n",
    "\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\cdot \\underbrace{\\mathbf{H}[l(\\mathbf{w})]^{-1}}_{\\text{inverse Hessian}}\\cdot\\underbrace{\\nabla l(\\mathbf{w})}_{\\text{gradient}}$$\n",
    "\n",
    "Adding in the exact equations for the Hessian and gradient, we can finally get:\n",
    "\n",
    "$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\cdot \\underbrace{\\left[\\mathbf{X}^T \\cdot \\text{diag}\\left[g(\\mathbf{X}\\cdot\\mathbf{w})\\odot(1-g(\\mathbf{X}\\cdot\\mathbf{w}))\\right] \\cdot \\mathbf{X} -2C \\right]^{-1} }_{\\text{inverse Hessian}} \\cdot \\underbrace{\\mathbf{X}\\odot y_{diff}}_{\\text{gradient}}$$\n",
    "\n",
    "\n",
    "You can see the full derivation of the Hessian in my hand written notes here also:\n",
    "- https://raw.githubusercontent.com/eclarson/MachineLearningNotebooks/master/PDF_Slides/HessianCalculation.pdf\n",
    "\n",
    "\n",
    "\n",
    "So let's code this up using numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary Logistic Regression Object with coefficients:\n",
      "[[-10.33521361]\n",
      " [ -1.11842524]\n",
      " [ -0.99722044]\n",
      " [  2.33435715]\n",
      " [  5.32207095]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 18.6 ms, sys: 4.39 ms, total: 23 ms\n",
      "Wall time: 10.7 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "from numpy.linalg import pinv\n",
    "class HessianBinaryLogisticRegression(BinaryLogisticRegression):\n",
    "    # just overwrite gradient function\n",
    "    def _get_gradient(self,X,y):\n",
    "        g = self.predict_proba(X,add_bias=False).ravel() # get sigmoid value for all classes\n",
    "        hessian = X.T @ np.diag(g*(1-g)) @ X - 2 * self.C # calculate the hessian\n",
    "\n",
    "        ydiff = y-g # get y difference\n",
    "        gradient = np.sum(X * ydiff[:,np.newaxis], axis=0) # make ydiff a column vector and multiply through\n",
    "        gradient = gradient.reshape(self.w_.shape)\n",
    "        gradient[1:] += -2 * self.w_[1:] * self.C\n",
    "        \n",
    "        return pinv(hessian) @ gradient\n",
    "       \n",
    "hlr = HessianBinaryLogisticRegression(eta=1.0,\n",
    "                                      iterations=4,\n",
    "                                      C=0.001) # note that we need only a few iterations here\n",
    "\n",
    "hlr.fit(X,y)\n",
    "yhat = hlr.predict(X)\n",
    "print(hlr)\n",
    "print('Accuracy of: ',accuracy_score(y,yhat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Can we still do better? Problems With the Hessian:\n",
    "Quadratic isn’t always a great assumption:\n",
    " - highly dependent on starting point\n",
    "  - jumps can get really random!\n",
    " - near saddle points, inverse hessian is unstable\n",
    " - hessian not always invertible… or invertible with enough numerical precision\n",
    " \n",
    "The Hessian can sometimes be ill formed for these problems and can also be highly computational. Thus, we prefer to approximate the Hessian, and approximate its inverse to better control the steps we make and directions we use.  \n",
    "\n",
    "<img src=\"PDF_slides/gru_hessian.jpg\"  width=\"400\">\n",
    "\n",
    "____\n",
    "\n",
    "# Quasi-Newton Methods\n",
    "In general:\n",
    " - Approximate the Hessian with something numerically sound and efficiently invertible \n",
    " - Back off to gradient descent when the approximate hessian is not stable\n",
    " - Try to create an approximation with as many properties of the Hessian as possible, like being symmetric and positive semi-definite\n",
    " - A popular approach: Rank One Hessian Approximation\n",
    " - An even more popular appraoch: Rank Two, with Broyden-Fletcher-Goldfarb-Shanno (BFGS)\n",
    "\n",
    "\n",
    "### Rank One Hessian Approximation \n",
    "Let's work our way up to using BFGS by first looking at one quasi-newton method, the rank one Hessian approximation. **Note, I only want you to get an intuition for this process. There is no requirenemtn to understand the derivation completely.**\n",
    "Essentially, we want to update the Hessian with an approximation that is easily invertible and based on stable gradient calculations. We can define the approximate Hessian for each iteration, $\\mathbf{H}_k$. To start as simple as possible, we will assume the Hessian can be approximated with one vector. Let's start off with a few other assumptions. We wil develop some equations that characterize a family of solutions. Within this family, we will only give one popular example for the rank one family solution and one popular example for the rank two family solution. \n",
    "\n",
    "___\n",
    "#### Step One: Use the secant equation\n",
    "One property of the hessian is called the Secant equation, which relates the change in input to the change in the derivative. The Secant Equation (exact for quadratic functions) is:\n",
    "$$ \\underbrace{\\mathbf{H}_{k+1}}_{\\text{approx. Hessian}} \\cdot\\underbrace{(\\mathbf{w}_{k+1} - \\mathbf{w}_k)}_{\\text{Change in }w} = \\underbrace{\\nabla l(\\mathbf{w}_{k+1}) - \\nabla l(\\mathbf{w}_k)}_{\\text{Change in gradient}}$$\n",
    "\n",
    "or, using intermediate variables for the differences: \n",
    "$$ \\mathbf{H}_{k+1} \\mathbf{s}_k =  \\mathbf{v}_k $$\n",
    "\n",
    "where $ \\mathbf{s} = (\\mathbf{w}_{k+1} - \\mathbf{w}_k) $ and $ \\mathbf{v} = (\\nabla l(\\mathbf{w}_{k+1}) - \\nabla l(\\mathbf{w}_k)) $, the difference in the gradient. If we enforce this relationship, we can find the hessian, assuming that the current location is approximated well by a quadratic (making the Secant a good assumption also). We also want the Hessian to be symmetric and not too far away from its initial value (for stable optimization).\n",
    "___\n",
    "\n",
    "#### Step Two: Use the rank one update assumption\n",
    "For optimizing, we would like to be able to update the Hessian at $\\mathbf{w}_{k+1}$ from our previous guess, $\\mathbf{H}_k$ at $\\mathbf{w}_{k}$, and have the update be easy to calculate. Therefore, we can choose the update of the Hessian to be approximated by the rank one update (one vector). Since the Hessian is the second partial derivative, a starting approximation might be the gradient difference we already defined such that  $\\mathbf{H} \\approx \\mathbf{v}\\cdot\\mathbf{v}^T$ which would form a matrix of the differences of each partial deriviative in the gradient. In practice, we need a vector that is slightly less constrained, such that    \n",
    "$$ \\mathbf{H}_{k+1}=\\mathbf{H}_k+\\alpha_k\\mathbf{u}\\cdot\\mathbf{u}^T $$\n",
    "Where $\\mathbf{u}$ and $\\alpha_k$ can be anything we want. \n",
    "\n",
    "Substituting back into the secant formula:\n",
    "$$ (\\mathbf{H}_k +\\alpha_k\\mathbf{u}_k\\cdot\\mathbf{u}_k^T)\\mathbf{s}_{k} = \\mathbf{v}_{k} $$\n",
    "\n",
    "\n",
    "Many solutions exist for this and they are referred to as a family of rank one Hessian approximations. One solution of this equation (there are many solutions) is to use one that simplifies nicely. For example we can choose the following:\n",
    "$$ \\mathbf{u}_k=\\mathbf{v}_{k}-\\mathbf{H}_k \\mathbf{s}_{k} \\text{   and   } \\alpha_k=\\frac{1}{(\\mathbf{v}_{k}-\\mathbf{H}_k \\mathbf{s}_{k})\\mathbf{s}_{k}}=\\frac{1}{\\mathbf{u}_{k}^T\\mathbf{s}_{k}}$$\n",
    "\n",
    "___\n",
    "#### Step Three: Combine secant equation with rank one update for iterative solution\n",
    "and combining this with our initial $\\mathbf{H}_{k+1}$ formula:\n",
    "$$ \\mathbf{H}_{k+1}=\\mathbf{H}_k- \\frac{\\mathbf{u}_k\\mathbf{u}_k^T}{\\mathbf{u}_k^T\\mathbf{s}_{k}} $$\n",
    "\n",
    "This gives an update for the Hessian, which we can use in our optimization formula. However, we need to define the vectors using the secant equation assumptions, such that $ \\mathbf{v}_k $ and $ \\mathbf{s}_k $ are the difference in gradients and weights as defined, respectively.\n",
    "\n",
    "#### Step Four: Use Sherman-Morris to also iterative track the inverse of our rank one approximation\n",
    "**Now for the power of this method.** We can now assume that the inverse of the Hessian can be optimized and formulate similar equations for its update, based upon the previous inverse. Therefore, we need the inverse of $(\\mathbf{H}_k+\\mathbf{v}\\cdot\\mathbf{v}^T)^{-1}$, which luckily has a closed form solution according to the Sherman-Morrison formula:\n",
    "\n",
    "$$ (\\mathbf{A}+\\mathbf{v}\\cdot\\mathbf{v}^T)^{-1} = \\mathbf{A}^{-1} - \\frac{\\mathbf{A}^{-1} \\mathbf{v} \\mathbf{v}^T\\mathbf{A}^{-1}}{1+\\mathbf{v}^T \\mathbf{A}^{-1} \\mathbf{v}}  $$\n",
    "___\n",
    "\n",
    "\n",
    "#### Summary of Rank One Approximation\n",
    "Now the optimization can be described as a rank one approximation of the Hessian. Placing it all together, we can get the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "|Description| Equations, Derivations, Hessian Calculations, and Miscellaneous |\n",
    "|-----------|--------|\n",
    "| **Definitions with Rank 1 Approximation** |  |\n",
    "|$$ \\mathbf{w} \\leftarrow \\mathbf{w} + \\eta \\cdot \\underbrace{\\mathbf{H}[l(\\mathbf{w})]^{-1}}_{\\text{inverse Hessian}}\\cdot\\underbrace{\\nabla l(\\mathbf{w})}_{\\text{gradient}}$$ |  |\n",
    "|1. Initial Approx. Hessian for $k=0$ is identity matrix| $$\\mathbf{H}_0=\\mathbf{I}$$|\n",
    "|2. Find update direction, $\\mathbf{p}_k$ | $$ \\mathbf{p}_k = -\\mathbf{H}_k^{-1} \\nabla l(\\mathbf{w}_k) $$| \n",
    "|3. Update $\\mathbf{w}$|$$\\mathbf{w}_{k+1}\\leftarrow \\mathbf{w}_k + \\eta \\cdot \\mathbf{p}_k $$|\n",
    "|4. Save scaled direction ($\\mathbf{w}_{k+1}-\\mathbf{w}_k$)| $$\\mathbf{s}_k=\\eta \\cdot \\mathbf{p}_k$$ |\n",
    "|5a. Approximate change in derivative | $$\\mathbf{v}_k = \\nabla l(\\mathbf{w}_{k+1}) - \\nabla l(\\mathbf{w}_k) $$|\n",
    "| 5b. Define $\\mathbf{u}$ from above: | $$\\mathbf{u}_k=\\mathbf{v}_k-\\mathbf{H}_k\\mathbf{s}_k$$| \n",
    "|6. Redefine approx Hessian update| $$\\mathbf{H}_{k+1}=\\mathbf{H}_k+\\underbrace{\\frac{\\mathbf{u}_k \\mathbf{u}_k^T}{\\mathbf{u}_k^T \\mathbf{s}_k}}_{\\text{approx. Hessian}} $$ |\n",
    "|7. Approx. Inverse $\\mathbf{H}_{k+1}^{-1}$ via Sherman Morris| $$ \\mathbf{H}_{k+1}^{-1} = \\mathbf{H}_{k}^{-1} - \\frac{\\mathbf{H}_k^{-1} \\mathbf{u}_k \\mathbf{u}_k^T\\mathbf{H}_k^{-1}}{1+\\mathbf{u}_k^T \\mathbf{H}_k^{-1} \\mathbf{u}_k} $$ |\n",
    "| 8. Repeat starting at step 2| $$ k = k+1 $$| \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Rank Two Hessian Approximation: BFGS\n",
    "\n",
    "\n",
    "Although the rank one approximation is a good performer, it can be improved by adding some additional criteria to the Hessian approximation. In this case, we assume that the $\\mathbf{H}_k$ needs to also be positive semi-definite, which helps with numerical stability. One of the most popular quasi-Newton methods that does this is known as Broyden-Fletcher-Goldfarb-Shanno (BFGS). \n",
    "- https://en.wikipedia.org/wiki/Broyden–Fletcher–Goldfarb–Shanno_algorithm \n",
    "\n",
    "<img src=\"PDF_slides/bfgs_meme.png\"  width=\"300\">\n",
    "\n",
    "#### Replace Step two with a higher rank approximation of the Hessian\n",
    "In this formulation we add an additional vector matrix addition to the update equation that ensure the resulting matrix is positive semidefinite:\n",
    "$$ \\mathbf{H}_{k+1}=\\mathbf{H}_k+\\alpha_k\\mathbf{u}_k\\cdot\\mathbf{u}_k^T + \\beta_k\\mathbf{z}_k\\cdot\\mathbf{z}_k^T$$\n",
    "\n",
    "The derivation is intuitively similar to the previous rank one approximation. Again, there are many potential solutions, referred to as a family of ranke two solutions. However, it becomes easier to obtain simple solutions for $\\mathbf{u}$ and $\\mathbf{z}$. The BFGS solution takes the following form:\n",
    "$$ \\mathbf{u}_k = \\mathbf{v}_k   \\text{    and    } \\mathbf{z}_k=\\mathbf{H}_k \\mathbf{s}_k  $$\n",
    "\n",
    "After solving for the $\\alpha_k$ and $\\beta_k$ coefficients, we get the update equation as:\n",
    "$$  \\mathbf{H}_{k+1}=\\mathbf{H}_k+\\underbrace{\\frac{\\mathbf{v}_k \\mathbf{v}_k^T}{\\mathbf{v}_k^T \\mathbf{s}_k}}_{\\text{previous}} -\\underbrace{\\frac{\\mathbf{H}_k \\mathbf{s}_k \\mathbf{s}_k^T \\mathbf{H}_k}{\\mathbf{s}_k^T \\mathbf{H}_k \\mathbf{s}_k}}_{\\text{new}}  $$\n",
    "\n",
    "The complete formulation can replace steps from the previous rank update as follows:\n",
    "\n",
    "|Description| Equations, Derivations, Hessian Calculations, and Miscellaneous |\n",
    "|-----------|--------|\n",
    "| **Alternative Definitions with Rank 2 (BFGS)** |  |\n",
    "| 6. Redefine approx Hessian| $$\\mathbf{H}_{k+1}=\\mathbf{H}_k+\\frac{\\mathbf{v}_k \\mathbf{v}_k^T}{\\mathbf{v}_k^T \\mathbf{s}_k} -\\frac{\\mathbf{H}_k \\mathbf{s}_k \\mathbf{s}_k^T \\mathbf{H}_k}{\\mathbf{s}_k^T \\mathbf{H}_k \\mathbf{s}_k} $$ |\n",
    "|7. Approximate Inverse $\\mathbf{H}_{k+1}^{-1}$ via Sherman Morris| $$ \\mathbf{H}_{k+1}^{-1} = \\mathbf{H}_{k}^{-1} + \\frac{(\\mathbf{s}_k^T \\mathbf{v}_k+\\mathbf{H}_{k}^{-1})(\\mathbf{s}_k \\mathbf{s}_k^T)}{(\\mathbf{s}_k^T \\mathbf{v}_k)^2}-\\frac{\\mathbf{H}_{k}^{-1} \\mathbf{v}_k \\mathbf{s}_k^T+\\mathbf{s}_k \\mathbf{v}_k^T\\mathbf{H}_{k}^{-1}}{\\mathbf{s}_k^T \\mathbf{v}_k} $$|\n",
    "\n",
    "___\n",
    "We won't explicitly program the BFGS algorithm--instead we can take advantage of scipy's calculations to do it for us. For using this algorithm, we need to define the objective function and the gradient explicitly for another program to calculate. \n",
    "Recall that Logistic regression uses the following objective function:\n",
    "\n",
    "$$ l(w) = \\left(\\sum_i y^{(i)} \\ln g(\\mathbf{x}^{(i)}) + (1-y^{(i)})\\ln[1-g(\\mathbf{x}^{(i)})]\\right)  - C \\cdot \\sum_j w_j^2 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary Logistic Regression Object with coefficients:\n",
      "[[-0.11105493]\n",
      " [-0.26703712]\n",
      " [-0.37959666]\n",
      " [ 0.4941737 ]\n",
      " [ 0.28604208]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 29.8 ms, sys: 3.39 ms, total: 33.2 ms\n",
      "Wall time: 8.23 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# for this, we won't perform our own BFGS implementation \n",
    "# (it takes a fair amount of code and understanding, which we haven't setup yet)\n",
    "# luckily for us, scipy has its own BFGS implementation:\n",
    "from scipy.optimize import fmin_bfgs # maybe the most common bfgs algorithm in the world\n",
    "from numpy import ma\n",
    "class BFGSBinaryLogisticRegression(BinaryLogisticRegression):\n",
    "    \n",
    "    @staticmethod\n",
    "    def objective_function(w,X,y,C):\n",
    "        g = expit(X @ w)\n",
    "        # invert this because scipy minimizes, but we derived all formulas for maximzing\n",
    "        return -np.sum(ma.log(g[y==1]))-np.sum(ma.log(1-g[y==0])) + C*sum(w**2) \n",
    "        #-np.sum(y*np.log(g)+(1-y)*np.log(1-g))\n",
    "\n",
    "    @staticmethod\n",
    "    def objective_gradient(w,X,y,C):\n",
    "        g = expit(X @ w)\n",
    "        ydiff = y-g # get y difference\n",
    "        gradient = np.mean(X * ydiff[:,np.newaxis], axis=0)\n",
    "        gradient = gradient.reshape(w.shape)\n",
    "        gradient[1:] += -2 * w[1:] * C\n",
    "        return -gradient\n",
    "    \n",
    "    # just overwrite fit function\n",
    "    def fit(self, X, y):\n",
    "        Xb = self._add_bias(X) # add bias term\n",
    "        num_samples, num_features = Xb.shape\n",
    "        \n",
    "        self.w_ = fmin_bfgs(self.objective_function, # what to optimize\n",
    "                            np.zeros((num_features,1)), # starting point\n",
    "                            fprime=self.objective_gradient, # gradient function\n",
    "                            args=(Xb,y,self.C), # extra args for gradient and objective function\n",
    "                            gtol=1e-03, # stopping criteria for gradient, |v_k|\n",
    "                            maxiter=self.iters, # stopping criteria iterations\n",
    "                            disp=False)\n",
    "        \n",
    "        self.w_ = self.w_.reshape((num_features,1))\n",
    "            \n",
    "bfgslr = BFGSBinaryLogisticRegression(_,iterations=2,C=0.001) # note that we need only a few iterations here\n",
    "\n",
    "bfgslr.fit(X,y)\n",
    "yhat = bfgslr.predict(X)\n",
    "print(bfgslr)\n",
    "print('Accuracy of: ',accuracy_score(y,yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#str(bfgslr.eta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BFGS and Newton's Method for Multiclass Logistic Regression\n",
    "Now let's add BFGS and the actual Hessian to non-binary classification. As before, we will use one-versus-all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# allow for the user to specify the algorithm they want to solver the binary case\n",
    "class MultiClassLogisticRegression:\n",
    "    def __init__(self, eta, iterations=20, \n",
    "                 C=0.0001, \n",
    "                 solver=BFGSBinaryLogisticRegression):\n",
    "        self.eta = eta\n",
    "        self.iters = iterations\n",
    "        self.C = C\n",
    "        self.solver = solver\n",
    "        self.classifiers_ = []\n",
    "        # internally we will store the weights as self.w_ to keep with sklearn conventions\n",
    "    \n",
    "    def __str__(self):\n",
    "        if(hasattr(self,'w_')):\n",
    "            return 'MultiClass Logistic Regression Object with coefficients:\\n'+ str(self.w_) # is we have trained the object\n",
    "        else:\n",
    "            return 'Untrained MultiClass Logistic Regression Object'\n",
    "        \n",
    "    def fit(self,X,y):\n",
    "        num_samples, num_features = X.shape\n",
    "        self.unique_ = np.sort(np.unique(y)) # get each unique class value\n",
    "        num_unique_classes = len(self.unique_)\n",
    "        self.classifiers_ = []\n",
    "        for i,yval in enumerate(self.unique_): # for each unique value\n",
    "            y_binary = np.array(y==yval).astype(int) # create a binary problem\n",
    "            # train the binary classifier for this class\n",
    "            \n",
    "            hblr = self.solver(eta=self.eta,iterations=self.iters,C=self.C)\n",
    "            hblr.fit(X,y_binary)\n",
    "\n",
    "            # add the trained classifier to the list\n",
    "            self.classifiers_.append(hblr)\n",
    "            \n",
    "        # save all the weights into one matrix, separate column for each class\n",
    "        self.w_ = np.hstack([x.w_ for x in self.classifiers_]).T\n",
    "        \n",
    "    def predict_proba(self,X):\n",
    "        probs = []\n",
    "        for hblr in self.classifiers_:\n",
    "            probs.append(hblr.predict_proba(X).reshape((len(X),1))) # get probability for each classifier\n",
    "        \n",
    "        return np.hstack(probs) # make into single matrix\n",
    "    \n",
    "    def predict(self,X):\n",
    "        return np.argmax(self.predict_proba(X),axis=1) # take argmax along row\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "\n",
    "ds = load_iris()\n",
    "X = ds.data\n",
    "\n",
    "#X = StandardScaler().fit(X).transform(X)\n",
    "y_not_binary = ds.target # note problem is NOT binary anymore, there are three classes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MultiClass Logistic Regression Object with coefficients:\n",
      "[[ 0.25736125  0.34295145  1.1457885  -1.78999942 -0.80675923]\n",
      " [ 0.26657802  0.35690326 -1.21580818  0.38117347 -0.77120072]\n",
      " [-0.95221952 -1.24281614 -1.11187691  1.8937359   1.60791217]]\n",
      "Accuracy of:  0.9666666666666667\n",
      "CPU times: user 177 ms, sys: 4.75 ms, total: 181 ms\n",
      "Wall time: 45.8 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lr = MultiClassLogisticRegression(eta=1,\n",
    "                                  iterations=10,\n",
    "                                  C=0.01,\n",
    "                                  solver=BFGSBinaryLogisticRegression\n",
    "                                 )\n",
    "lr.fit(X,y_not_binary)\n",
    "print(lr)\n",
    "\n",
    "yhat = lr.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y_not_binary,yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MultiClass Logistic Regression Object with coefficients:\n",
      "[[  5.31819122   0.39429748   3.02935592  -4.53025718  -6.1168987 ]\n",
      " [  7.37229534  -0.24420762  -2.79505273   1.30971959  -2.76989863]\n",
      " [-35.80587165  -2.45860804  -5.49495125   8.41245733  15.17916522]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 27.2 ms, sys: 3.49 ms, total: 30.6 ms\n",
      "Wall time: 8.2 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lr = MultiClassLogisticRegression(eta=1,\n",
    "                                  iterations=10,\n",
    "                                  C=0.001,\n",
    "                                  solver=HessianBinaryLogisticRegression\n",
    "                                 )\n",
    "lr.fit(X,y_not_binary)\n",
    "print(lr)\n",
    "\n",
    "yhat = lr.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y_not_binary,yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.69103463  3.31824182 -6.53118932 -3.36892048]\n",
      " [-0.24478421 -2.79580595  1.31167976 -2.77411255]\n",
      " [-2.45880535 -6.00752588  8.8977078  16.55360299]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 325 ms, sys: 26.4 ms, total: 351 ms\n",
      "Wall time: 98.5 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# how do we compare now to sklearn?\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "lr_sk = LogisticRegression(solver='lbfgs',n_jobs=1,\n",
    "                           multi_class='ovr', C = 1/0.001, \n",
    "                           penalty='l2',\n",
    "                          max_iter=50) # all params default\n",
    "# note that sklearn is optimized for using the liblinear library with logistic regression\n",
    "# ...and its faster than our implementation here\n",
    "\n",
    "lr_sk.fit(X, y_not_binary) # no need to add bias term, sklearn does it internally!!\n",
    "print(lr_sk.coef_)\n",
    "yhat = lr_sk.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y_not_binary,yhat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.2278196   4.12445837 -6.53269249 -3.16821338]\n",
      " [-0.23894082 -2.78878647  1.30959659 -2.77331609]\n",
      " [-2.74173431 -5.4303242   7.64423933 14.01816119]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 19.1 ms, sys: 2.74 ms, total: 21.9 ms\n",
      "Wall time: 5.53 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# actually, we aren't quite as good as the lib linear implementation\n",
    "# how do we compare now to sklearn?\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "lr_sk = LogisticRegression(solver='liblinear',n_jobs=1, \n",
    "                           multi_class='ovr', C = 1/0.001, \n",
    "                           penalty='l2',max_iter=100) \n",
    "\n",
    "lr_sk.fit(X,y_not_binary) # no need to add bias term, sklearn does it internally!!\n",
    "print(lr_sk.coef_)\n",
    "yhat = lr_sk.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y_not_binary,yhat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Liblinear is a great toolkit for linear modeling (from national Taiwan University) and the paper can be found here:\n",
    "- https://www.csie.ntu.edu.tw/~cjlin/papers/liblinear.pdf\n",
    "\n",
    "Actually, this solves a slightly different problem (known as the 'dual' formulation because its a linear SVM) to make it extremely fast. So this is actually not a fair comparison to ours."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallelization across Classes [if time]\n",
    "Is the formulation of multi-class logistic regression parallelizable? Is the optimization? Or the one versus all classification?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MultiClass Logistic Regression Object with coefficients:\n",
      "[[  5.31819122   0.39429748   3.02935592  -4.53025718  -6.1168987 ]\n",
      " [  7.37229534  -0.24420762  -2.79505273   1.30971959  -2.76989863]\n",
      " [-35.80587165  -2.45860804  -5.49495125   8.41245733  15.17916522]]\n",
      "Accuracy of:  0.98\n",
      "CPU times: user 621 ms, sys: 130 ms, total: 751 ms\n",
      "Wall time: 107 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# its still faster! Can we fix it with parallelization?\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "class ParallelMultiClassLogisticRegression(MultiClassLogisticRegression):\n",
    "    @staticmethod\n",
    "    def par_logistic(yval,eta,iters,C,X,y,solver):\n",
    "        y_binary = y==yval # create a binary problem\n",
    "        # train the binary classifier for this class\n",
    "        hblr = solver(eta=eta,iterations=iters,C=C)\n",
    "        hblr.fit(X,y_binary)\n",
    "        return hblr\n",
    "    \n",
    "    def fit(self,X,y):\n",
    "        num_samples, num_features = X.shape\n",
    "        self.unique_ = np.sort(np.unique(y)) # get each unique class value\n",
    "        num_unique_classes = len(self.unique_)\n",
    "        backend = 'threading' # can also try 'multiprocessing'\n",
    "        \n",
    "        self.classifiers_ = Parallel(n_jobs=-1,backend=backend)(\n",
    "            delayed(self.par_logistic)(yval,self.eta,self.iters,self.C,X,y,self.solver) for yval in self.unique_)\n",
    "            \n",
    "        # save all the weights into one matrix, separate column for each class\n",
    "        self.w_ = np.hstack([x.w_ for x in self.classifiers_]).T\n",
    "        \n",
    "plr = ParallelMultiClassLogisticRegression(eta=1,iterations=10,C=0.001,solver=HessianBinaryLogisticRegression)\n",
    "plr.fit(X,y_not_binary)\n",
    "print(plr)\n",
    "\n",
    "yhat = plr.predict(X)\n",
    "print('Accuracy of: ',accuracy_score(y_not_binary,yhat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "Please note that the overhead of parallelization is not worth it for this problem!!\n",
    "\n",
    "**When would it make sense???**\n",
    "___\n",
    "\n",
    "\n",
    "<img src=\"PDF_slides/mark_scooter.png\"  width=\"300\">\n",
    "\n",
    "\n",
    "# Extended Logistic Regression Example\n",
    "\n",
    "In this example we will explore methods of using logistic regression in scikit-learn. A basic understanding of scikit-learn is required to complete this notebook, but we start very basic. Note also that there are more efficient methods of separating testing and training data, but we will leave that for a later lecture.\n",
    "\n",
    "First let's load a dataset and prepare it for analysis. We will use pandas to load in data, and then prepare it for classification. We will be using the titanic dataset (a very modest sized data set of about 1000 instances). The imputation methods used here are discussed in a previous notebook.\n",
    "\n",
    "Steps:\n",
    "- Load data, impute\n",
    "- One hot encode and normalize data\n",
    "- Separate into training and testing sets\n",
    "- Explore best hyper parameter, C\n",
    "\n",
    "\n",
    "## Load Titanic Data and Pre-process\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 882 entries, 0 to 890\n",
      "Data columns (total 8 columns):\n",
      " #   Column    Non-Null Count  Dtype  \n",
      "---  ------    --------------  -----  \n",
      " 0   Survived  882 non-null    int64  \n",
      " 1   Age       882 non-null    float64\n",
      " 2   Sex       882 non-null    object \n",
      " 3   Parch     882 non-null    int64  \n",
      " 4   SibSp     882 non-null    int64  \n",
      " 5   Pclass    882 non-null    int64  \n",
      " 6   Fare      882 non-null    float64\n",
      " 7   Embarked  882 non-null    object \n",
      "dtypes: float64(2), int64(4), object(2)\n",
      "memory usage: 62.0+ KB\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ericlarson/opt/anaconda3/envs/mlenv2022/lib/python3.7/site-packages/ipykernel_launcher.py:17: FutureWarning: Dropping invalid columns in DataFrameGroupBy.transform is deprecated. In a future version, a TypeError will be raised. Before calling .transform, select only columns which should be valid for the transforming function.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.read_csv('data/titanic.csv') # read in the csv file\n",
    "\n",
    "# 1. Remove attributes that just arent useful for us\n",
    "del df['PassengerId']\n",
    "del df['Name']\n",
    "del df['Cabin']\n",
    "del df['Ticket']\n",
    "\n",
    "# 2. Impute some missing values, grouped by their Pclass and SibSp numbers\n",
    "df_grouped = df.groupby(by=['Pclass','SibSp'])\n",
    "\n",
    "# # now use this grouping to fill the data set in each group, then transform back\n",
    "# fill in the numeric values\n",
    "df_imputed = df_grouped.transform(lambda grp: grp.fillna(grp.median()))\n",
    "# fill in the categorical values\n",
    "df_imputed[['Sex','Embarked']] = df_grouped[['Sex','Embarked']].apply(lambda grp: grp.fillna(grp.mode()))\n",
    "# fillin the grouped variables from original data frame\n",
    "df_imputed[['Pclass','SibSp']] = df[['Pclass','SibSp']]\n",
    "\n",
    "# 4. drop rows that still had missing values after grouped imputation\n",
    "df_imputed.dropna(inplace=True)\n",
    "\n",
    "# 5. Rearrange the columns\n",
    "df_imputed = df_imputed[['Survived','Age','Sex','Parch','SibSp','Pclass','Fare','Embarked']]\n",
    "\n",
    "df_imputed.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 882 entries, 0 to 890\n",
      "Data columns (total 12 columns):\n",
      " #   Column      Non-Null Count  Dtype  \n",
      "---  ------      --------------  -----  \n",
      " 0   Survived    882 non-null    int64  \n",
      " 1   Age         882 non-null    float64\n",
      " 2   Parch       882 non-null    int64  \n",
      " 3   SibSp       882 non-null    int64  \n",
      " 4   Pclass      882 non-null    int64  \n",
      " 5   Fare        882 non-null    float64\n",
      " 6   Embarked_C  882 non-null    uint8  \n",
      " 7   Embarked_Q  882 non-null    uint8  \n",
      " 8   Embarked_S  882 non-null    uint8  \n",
      " 9   Sex_female  882 non-null    uint8  \n",
      " 10  Sex_male    882 non-null    uint8  \n",
      " 11  FamilySize  882 non-null    int64  \n",
      "dtypes: float64(2), int64(5), uint8(5)\n",
      "memory usage: 59.4 KB\n"
     ]
    }
   ],
   "source": [
    "# perform one-hot encoding of the categorical data \"embarked\"\n",
    "tmp_df = pd.get_dummies(df_imputed.Embarked,prefix='Embarked')\n",
    "df_imputed = pd.concat((df_imputed,tmp_df),axis=1) # add back into the dataframe\n",
    "\n",
    "tmp_df = pd.get_dummies(df_imputed.Sex,prefix='Sex')\n",
    "df_imputed = pd.concat((df_imputed,tmp_df),axis=1) # add back into the dataframe\n",
    "\n",
    "# Now let's clean up the dataset\n",
    "if 'Sex' in df_imputed:\n",
    "    del df_imputed['Sex'] # if 'Sex' column still exists, delete it (as we created an ismale column)\n",
    "    \n",
    "if 'Embarked' in df_imputed:    \n",
    "    del df_imputed['Embarked'] # get reid of the original category as it is now one-hot encoded\n",
    "    \n",
    "# Finally, let's create a new variable based on the number of family members\n",
    "# traveling with the passenger\n",
    "\n",
    "# notice that this new column did not exist before this line of code--we use the pandas \n",
    "#    syntax to add it in \n",
    "df_imputed['FamilySize'] = df_imputed.Parch + df_imputed.SibSp\n",
    "df_imputed.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training and Testing Split\n",
    "For training and testing purposes, let's gather the data we have and grab 80% of the instances for training and the remaining 20% for testing. Moreover, let's repeat this process of separating the testing and training data three times. We will use the hold out cross validation method built into scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ShuffleSplit(n_splits=3, random_state=None, test_size=0.2, train_size=None)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import ShuffleSplit\n",
    "\n",
    "# we want to predict the X and y data as follows:\n",
    "if 'Survived' in df_imputed:\n",
    "    y = df_imputed['Survived'].to_numpy() # get the labels we want\n",
    "    del df_imputed['Survived'] # get rid of the class label\n",
    "    \n",
    "    norm_features = ['Age','Fare' ]\n",
    "    df_imputed[norm_features] = (df_imputed[norm_features]-df_imputed[norm_features].mean()) / df_imputed[norm_features].std()\n",
    "    X = df_imputed.to_numpy() # use everything else to predict!\n",
    "\n",
    "    ## X and y are now numpy matrices, by calling 'values' on the pandas data frames we\n",
    "    #    have converted them into simple matrices to use with scikit learn\n",
    "    \n",
    "    \n",
    "# to use the cross validation object in scikit learn, we need to grab an instance\n",
    "#    of the object and set it up. This object will be able to split our data into \n",
    "#    training and testing splits\n",
    "num_cv_iterations = 3\n",
    "num_instances = len(y)\n",
    "cv_object = ShuffleSplit(\n",
    "                         n_splits=num_cv_iterations,\n",
    "                         test_size  = 0.2)\n",
    "                         \n",
    "print(cv_object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====Iteration 0  ====\n",
      "accuracy 0.8135593220338984\n",
      "confusion matrix\n",
      " [[95 16]\n",
      " [17 49]]\n",
      "====Iteration 1  ====\n",
      "accuracy 0.847457627118644\n",
      "confusion matrix\n",
      " [[105   7]\n",
      " [ 20  45]]\n",
      "====Iteration 2  ====\n",
      "accuracy 0.7853107344632768\n",
      "confusion matrix\n",
      " [[98 17]\n",
      " [21 41]]\n"
     ]
    }
   ],
   "source": [
    "# run logistic regression and vary some parameters\n",
    "from sklearn import metrics as mt\n",
    "\n",
    "# first we create a reusable logisitic regression object\n",
    "#   here we can setup the object with different learning parameters and constants\n",
    "lr_clf = HessianBinaryLogisticRegression(eta=0.1,iterations=10) # get object\n",
    "\n",
    "# now we can use the cv_object that we setup before to iterate through the \n",
    "#    different training and testing sets. Each time we will reuse the logisitic regression \n",
    "#    object, but it gets trained on different data each time we use it.\n",
    "\n",
    "iter_num=0\n",
    "# the indices are the rows used for training and testing in each iteration\n",
    "for train_indices, test_indices in cv_object.split(X,y): \n",
    "    # I will create new variables here so that it is more obvious what \n",
    "    # the code is doing (you can compact this syntax and avoid duplicating memory,\n",
    "    # but it makes this code less readable)\n",
    "    X_train = X[train_indices]\n",
    "    y_train = y[train_indices]\n",
    "    \n",
    "    X_test = X[test_indices]\n",
    "    y_test = y[test_indices]\n",
    "    \n",
    "    # train the reusable logisitc regression model on the training data\n",
    "    lr_clf.fit(X_train,y_train)  # train object\n",
    "    y_hat = lr_clf.predict(X_test) # get test set precitions\n",
    "\n",
    "    # now let's get the accuracy and confusion matrix for this iterations of training/testing\n",
    "    acc = mt.accuracy_score(y_test,y_hat)\n",
    "    conf = mt.confusion_matrix(y_test,y_hat)\n",
    "    print(\"====Iteration\",iter_num,\" ====\")\n",
    "    print(\"accuracy\", acc )\n",
    "    print(\"confusion matrix\\n\",conf)\n",
    "    iter_num += 1\n",
    "    \n",
    "# Also note that every time you run the above code\n",
    "#   it randomly creates a new training and testing set, \n",
    "#   so accuracy will be different each time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "====Iteration 0  ====\n",
      "accuracy 0.7966101694915254\n",
      "confusion matrix\n",
      " [[92 15]\n",
      " [21 49]]\n",
      "====Iteration 1  ====\n",
      "accuracy 0.751412429378531\n",
      "confusion matrix\n",
      " [[86 19]\n",
      " [25 47]]\n",
      "====Iteration 2  ====\n",
      "accuracy 0.7627118644067796\n",
      "confusion matrix\n",
      " [[88 17]\n",
      " [25 47]]\n"
     ]
    }
   ],
   "source": [
    "# this does the exact same thing as the above block of code, but with shorter syntax\n",
    "\n",
    "for iter_num, (train_indices, test_indices) in enumerate(cv_object.split(X,y)):\n",
    "    lr_clf.fit(X[train_indices],y[train_indices])  # train object\n",
    "    y_hat = lr_clf.predict(X[test_indices]) # get test set precitions\n",
    "\n",
    "    # print the accuracy and confusion matrix \n",
    "    print(\"====Iteration\",iter_num,\" ====\")\n",
    "    print(\"accuracy\", mt.accuracy_score(y[test_indices],y_hat)) \n",
    "    print(\"confusion matrix\\n\",mt.confusion_matrix(y[test_indices],y_hat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Example, adjusting C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c71a59a63f044c13bbe3e62c3de636a7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='cost', options=(0.0001, 0.00022758459260747887, 0.0005179474679231…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.lr_explor(cost)>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ipywidgets import widgets as wd\n",
    "\n",
    "num_cv_iterations = 10\n",
    "num_instances = len(y)\n",
    "cv_object = ShuffleSplit(n_splits=num_cv_iterations,\n",
    "                         test_size  = 0.5)\n",
    "\n",
    "def lr_explor(cost):\n",
    "    print('Running')\n",
    "    lr_clf = HessianBinaryLogisticRegression(eta=0.1,iterations=10,\n",
    "                                            C=float(cost)) # get object\n",
    "    acc = []\n",
    "    for iter_num, (train_indices, test_indices) in enumerate(cv_object.split(X,y)):\n",
    "        lr_clf.fit(X[train_indices],y[train_indices])  # train object\n",
    "        y_hat = lr_clf.predict(X[test_indices]) # get test set predictions\n",
    "        acc.append(mt.accuracy_score(y[test_indices],y_hat))\n",
    "        \n",
    "    acc = np.array(acc)\n",
    "    print(acc.mean(),'+-',2.7*acc.std())\n",
    "        \n",
    "wd.interact(lr_explor,cost=list(np.logspace(-4,1,15)),__manual=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exhaustive Search for C, then Visualize with Boxplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.67 s, sys: 111 ms, total: 2.78 s\n",
      "Wall time: 2.82 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# alternatively, we can also graph out the values using boxplots\n",
    "num_cv_iterations = 20\n",
    "num_instances = len(y)\n",
    "cv_object = ShuffleSplit(n_splits=num_cv_iterations,\n",
    "                         test_size  = 0.5)\n",
    "\n",
    "def lr_explor(cost):\n",
    "    lr_clf = BFGSBinaryLogisticRegression(eta=0.1,iterations=10,\n",
    "                                            C=float(cost)) # get object\n",
    "    acc = []\n",
    "    for iter_num, (train_indices, test_indices) in enumerate(cv_object.split(X,y)):\n",
    "        lr_clf.fit(X[train_indices],y[train_indices])  # train object\n",
    "        y_hat = lr_clf.predict(X[test_indices]) # get test set predictions\n",
    "        acc.append(mt.accuracy_score(y[test_indices],y_hat))\n",
    "        \n",
    "    acc = np.array(acc)\n",
    "    return acc\n",
    "\n",
    "costs = np.logspace(-5,1,20)\n",
    "accs = []\n",
    "for c in costs:\n",
    "    accs.append(lr_explor(c))\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEpCAYAAABiNA5uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt0UlEQVR4nO3de3zcVZ3/8VcgtYJc4m4EbQuKGkC76xXjrvl5K4LYFYm3jwW8oSuyLgSVSPWndXXrpYVwi4KIXF25/D7eAusWBKmiRt0UEJUCSwQV2opYMIsgQgv5/XG+k06mSeZM5jszJzPv5+ORR2e+M+czn+8knc98zznf820bHx9HRERkJjs0OgEREUmfioWIiJSlYiEiImWpWIiISFkqFiIiUpaKhYiIlNXe6ARqSHOCRUQq1zbVxmYuFmzatGnaxzo7O9m8eXNV8auNkUIOqcRIIYdUYqSQQyoxUsghlRj1yGHBggXTPqZuKBERKUvFQkREylKxEBGRslQsRESkLBULEREpS8VCRETKUrEQkUmGhoZYsmQJO+20E0uWLGFoaKjRKUkCmvo8CxGpzNDQEKtXr2ZgYIClS5eyZs0a+vv7Aejt7W1sctJQOrKQhlq4cCELFy5k/vz5E7elcQYHBxkYGKCnp4d58+bR09PDwMAAg4ODjU5NGkxHFtJQGzduBELRKNyWxhkdHaW7u3vStu7ubkZHRxuUkaSi5Y4s1B8rMr2uri5GRkYmbRsZGaGrq6tBGUkqWqpYFPpjV65cyQMPPMDKlStZvXq1CoY0jWq79fr6+ujv72d4eJgtW7YwPDxMf38/fX19NcpY5oqW6oaarj92xYoVGryTplBtt17h/8GKFStYtmwZXV1dLF++XP8/pLWKhfpjm9NU354bMf5RbR6p7Edvby+9vb25rHIqzaOlikWhP7anp2dim/pj575UBsmrzSOV/RCZSkuNWeTVH6vpniLSalrqyCKv/thqvwFOV1zq3WWRSreHiKSvpYoFpNEfm0d3QyoxRKQ1tFyxEJmKjrJEZqZiIYKOskTKUbGQqugbuUhrULGQqugbuUhraLlioW/CkqfFixczNja23fbiv7OOjg7Wr19fx6yqo/8jMpWWKxaN/CY81QdL6X/Mch8s5WLMtQ+muW5sbGy7v6PSmXZz7TwcHS3KVFquWDRS6QfLVNN3y32wlIsx1z6YRGRuqFuxMLNDgDOAHYFz3X1VyeO7A18D9s7yGnD3C2LaisxleRxxitRaXYqFme0InAkcBGwA1pnZFe5+S9HT/hW4xd0PNbOnAP9jZhcDj0W0lQo0S1dWHvuRwnuRxxGnSK3V68iiG/i1u98JYGaXAYcBxR/448CuZtYG7ALcD2wFXhrRti5S+GDJQ7VdWakM6ubRJdcs3XrN8rcp6apXsVgI3F10fwOhCBT7InAFsAnYFXibuz9uZjFty8rjP1MKHywHnHQth11824yP11oeg7rqeslXCn+b0tzqVSzaptg2XnL/tcBNwBLgWcA1ZvajyLYAmNnRwNEA7k5nZ+fEY2NjYzzyyCMT99vb29m6devE/fnz5096/nSKn9Pe3r5dm3IxyrUvF+P6Ew8svx/HPzJV04ryqGQ/ZhOj3O8D4n4n1e5HHjGqfS9icsgjRsx7UelrljPdvtSrfTPFaHQO9SoWG4C9iu4vIhxBFDsKWOXu48Cvzew3wP6RbQFw93OAc7K746X9vsX3p+oXjllYsNoY5drnEaPe+5FHjEa9FwecdC09Z/x4xsdb5b0oF3M2ql2wM48FP5slRj1yWLBgwbSP1atYrAO6zGwfYCOwDDii5Dl3AQcCPzKzPYH9gDuBsYi2dZFCF5Dk6/oTDyzffXPk9OcalPubKDxHZK6rS7Fw961mdizwXcL01/Pdfb2ZHZM9fjawErjQzH5F6Hpa7u6bAaZqW4+8SzXLB0u1RS+P/UjlvahW6d8EVP53ITIX1O08C3dfA6wp2XZ20e1NwMGxbSuVwlFBuWID9flgqbbo5fEBmcp7ISJxWuYM7mo/IEVEWlnLFIu8zDQFsaOjo36JSBJHi6nQeyG1pmJRgdKuFy201lg6WtxG74XUmoqFtLxGHy02y2C/NLeWKhaN/lBoJuXOCJ4r72cKR4sa7Je5oGWKRQofCoXXnUnMh2yji95U71uK72e9ClazFE6RmbRMsUhBHgUrlaKXghTei5QKZ7W0GKHMpOWKRfEff+H2XPyPXa0UvpFLWrQYocyk5YpF4T9DHuuszFUpfCMvvO5MVLRE0tFyxULSkErBkm10robMRMVCBHVPgs7VkJmpWIjQ+O5JdclJ6lQsZkHfQmUqs/27yKtLLpVJC9Plof8jc5uKxSw0+luopKmRfxcpjQEVv67GoprHDo1OQERE0hdVLMzsb2udSCtZuHDhxKF64bbmsItIymK7oe42s2uA/wCucPdHa5hT01M3lojMNbHF4unA4cBy4Bwz+wbwVXef/kr3krw8Buo12J8fvZeSsqhi4e5/BAaBQTPbD3gH8B9mNg58DTjP3X9XuzSlFvI4wjnzzDMZHBxkdHSUrq4u+vr68kyxpeiIU1I2mwHup2Y/uwF3AAuBn5vZR/NMTNI3NDTE6tWrWblyJQ888AArV65k9erVDA0NNTo1EclZ1JGFmS0G3g4cCTwIXAQ8z903Zo+vBH4JrKpRnpKgwcFBBgYG6OnpYd68efT09DAwMMCKFSvo7e1tdHoikqPYMYsfApcCb3H3kdIH3f23ZnZ6nolJ+kZHR+nu7p60rbu7m9HR0brmkUdf/9DQ0HbdaXOx4GncQ2oltlg8rdwMKHf/ZA75yBzS1dXFyMgIPT09E9tGRkbo6uqqax7V9vUXutMGBgZYunQpa9asob+/H2DOFQyNe0itxI5ZDJjZy4o3mNnLdDTR2vr6+ujv72d4eJgtW7YwPDxMf3//nBvkHhwcpLe3lxUrVrDbbrtNdKMNDg42OjWRZMQeWRwO9JdsuwEYAj6YYz4yhxS+da9YsYJly5bR1dXF8uXLK/o2nkK3ye23385f/vIXTjnllIkjixNOOIENGzbUNY8UpLK+lKQntliMs/1RyI5TbJMW09vbS29v76y7PVLoNpk3bx5HHXXUpIH6o446ilWrWmu+RkrrS0l6Yj/sfwR8xsx2AMj+/VS2XWRO27JlCxdccMGk7rQLLriALVu2NDo1kWTEHlkcD3wH+L2Z/Q7YG/g9cGitEpOZpdB90yz23XdfDjnkkEndaW9605u46qqrGp2aSDJiz+DeYGYvAl4KLALuBkbc/fFaJifTS6H7pln09fVNORtq+fLljU5tTlm8eDFjY2PbbS/+YtPR0cH69evrmJXkJfp6Fllh+OlsX8jMDgHOIIx1nOvuq0oe/wjhpL9CXs8BnuLu95vZb4E/A48BW939gNnmIVIqj4F6gbGxse2Obqe8NKvMSbFncO9GGKN4JdAJtBUec/e9I9rvCJwJHARsANaZ2RXufktRnJOBk7PnHwp8yN3vLwrzanfXV2ipiWoH6kWaXewA91nAi4B/B/4GOA64Czgtsn038Gt3vzM7ue8y4LAZnn844YxxERFJQGyxOBh4s7tfDjyW/fs2wuqzMRYSxjkKNmTbtmNmOwOHAN8s2jwOXG1mN5jZ0ZGvKSIiOYkds9gB+N/s9oNm1kGYDfXsyPZtU2wbn+a5hwLDJV1QPe6+ycz2AK4xs9vc/YelDbNCcjSAu9PZ2TltQu3t7TM+HqPaGCnkkEqMFHJIJUYKORRUGqP0+VPlUUnMVN6LFGI0OofYYvELwnjFtYRzK84krD57e2T7DcBeRfcXAZumee4ySrqg3H1T9u+9ZvZtQrfWdsXC3c8Bzsnujs/U95xH33S1MVLIIZUYKeSQSowUciioNEbp86fKo5KYqbwXKcSoRw4LFiyY9rHYYvE+th0d9AGfBzqAd0a2Xwd0mdk+wEZCQTii9ElmtjuhKL29aNuTgB3c/c/Z7YMJYyciIlInZYtFNpPp3cBnYeKqef9cyYu4+1YzOxb4LmHq7Pnuvt7MjskePzt76huBq939oaLmewLfNrNCvpe4u86WEqkRnfApUylbLNz9MTP7V8LU2Vlz9zXAmpJtZ5fcvxC4sGTbncDzq3ltmVqzXMNB8qUTPmUqsbOhLgKOqWUiUl+6JKqIVCJ2zKIbOM7MTiRMgZ2YyeTur6hFYlJbuiSqiFQitlh8JfuRJpHKJVGleRxw0rUcdvFtZZ9TL1MtLaKxl9mLXUjwolonIvWVyiVRpXlcf+KBcWtDHVmfD+xCLrouRz5i14Z6z3SPufv5+aUj9VK4JGphpdXCJVG10qqkQJMv0hPbDVW6rMdTgWcBw4CKxRyklVYlVYXJF6VLxgP6+2yg2G6oV5duy442npN7RlI3WmlVUqTJF2mq5hraFwLvzSkPERFAky9SFTtmUVpUdiYsyTGWd0Ii0to0+SJNsWMWW9l+ldiNZCu8iohA+SvhdXR0lI2hyRdpii0W+5Tcf0hXrRORYlNNT53NtFVNvkhTJUcWf3H3PxU2mNmTgZ0Ky4eLiORFky/SEzvAPUS4BkWxRcC3c81GRESSFFss9nP3XxVvyO7vn39KIiKSmthica+ZTbqEanb/vvxTEhGR1MSOWZwPfNPMPg7cSTh7eyVwbq0SExGRdMQWi1XAFmCAcC3tu4DzgFNrlJeIiCQkdrmPx4GTsx8RkZrS8uLpiRqzMLOPmtlLSrZ1ZxdDEhHJ1caNGyeKQ/FtaZzYAe7jgVtKtt0CfDDXbEREJEmxxeIJhDGLYo8CT8w3HRERSVFssbgB+EDJtmOAG/NNR0REUhQ7G+pDwDVm9g7gDuDZwJ7AQbVKTERE0hF1ZOHu64F9CbOh1gEnEc7qLh3HEBGRJhR7ZIG7PwhcVsNcREQkUbEXP2onjFm8EugE2gqPufsrapOaiIikInaA+zTg/cAPgRcD3wT2ANbWKC8REUlIbLF4E/A6dz8D2Jr92wu8ulaJiYhIOmKLxc7A3dnth81sZ3e/DXhhbdISEZGUxA5w3wq8BBgBrgc+ZWYPEK7DLSIiTS62WBwPPJbd/jDwJWBX4OjYFzKzQ4AzgB2Bc919VcnjHwGOLMrrOcBT3P3+cm1FRKS2YledXVd0exR4TSUvYmY7AmcSTuLbAKwzsyuKz9Nw94lVbc3sUOBDWaEo21ZE0lG6YmzhvhYDnNuiz7OoUjfwa3e/E8DMLgMOY/vFCQsOBy6dZVsRaaDiotDZ2cnmzZsbmI3kpV7FYiHbBsghHCG8dKonmtnOwCHAsbNoezRZ15i709nZOW1C7e3tMz4eo9oYKeSQSowUckglRgo5pBSj2vZ5xEjhvWh0DvUqFm1TbBuf5rmHAsPufn+lbd39HOCcwnNm+kaTxzeeamOkkEMqMVLIIZUYKeSQUow8jkxS2I+58HexYMGCaR+LnTpbrQ2Ey7EWLAI2TfPcZWzrgqq0rYiI1ED0kYWZHQy8ANileLu7fzKi+Tqgy8z2IUy3XQYcMcVr7E5YUuTtlbYVEZHaib2s6heBrxGW+tir6GdRTHt330oYg/gu4ZwNd/f1ZnaMmR1T9NQ3Ale7+0Pl2sa8roiI5CP2yOJw4AXufnfZZ07D3dcAa0q2nV1y/0Lgwpi2IiJTWbx4MWNjY5O2FU/n7ejoYP16fd+sVGyxuA8Yq2EeIiK5GBsbm3H6bul5IBIntlicAlxsZp8H/lD8QOH8BxGRaumoIF2xxeJL2b+vL9k+TliCQ0SkajoqSFfsch/1mmIrIpKMoaEhBgcHGR0dpauri76+Pnp7exudVkNUdFKeme1NOKN6QzWD3SIiqRsaGmL16tUMDAywdOlS1qxZQ39/P0BLFozYqbNPM7PrgF8D3wLuMLMfmtn0p/uJiMxhg4ODDAwM0NPTw7x58+jp6WFgYIDBwcFGp9YQsd1LXwJ+ATzZ3Z8GPBn4OXD2jK1EROao0dFRuru7J23r7u5mdHS0QRk1Vmyx+D/ACYWT5bJ/TwReVqvEREQaqauri5GRkUnbRkZG6OrqalBGjRU7ZvEn4LmEo4uC/dC5FyKSowNOupbDLr5txsfrpa+vj/7+/okxi+HhYfr7+1m+fHndckhJbLE4CfiemZ0H/A54OnAUsKJWiYlI67n+xAPLT509sj4XUSoMYq9YsYJly5bR1dXF8uXLW3JwGyK7odz9K8DbgE7CEuKdwOHZkuAiIpMMDQ2xZMkSdtppJ5YsWcLQ0FCjU5qV3t5e1q5dy8MPP8zatWtbtlBABVNn3X0tsLaGuYhIE2j0lNOUurKaybTFwsw+7u6fzW7/+3TPi1yiXERaxHRTTlesWFGXYpFSV1YzmenIonj58b2mfZaISBFNOW1O0xYLd/+XottH1ScdEZnrClNOe3p6Jra18pTTZhF7Bvf902y/N990RGSuK0w5HR4eZsuWLRNTTvv6+hqdmlQhdoB7XukGM5uHVpwVkRKactqcZiwWZvYjwjLkTzSzH5Y8vAj4Sa0SE5G5q7e3l97e3u0Gl1tRs6xcW+7I4lygDXgJcF7R9nHCRZA0lVZEmtZU188onmlVTqOnEedpxmLh7hcBmNnP3H36icsiIk2oUBgWLlxYUZEoaPQ04jzFXvzoNjPbE+gmnL3dVvTY+TXKTURkTmumacSxs6F6gTuAfwe+DByX/fuOmmUmIjLHNdPKtbFLlH8GOMrdXwg8lP17NHBDzTITkZa0cOHCiZ/58+dPut/R0dHo9CrSTNOIY6fO7u3uXy/ZdhFwD9Cfb0oi0qpKxwVmO1aQimaaRhxbLO41sz3d/Q/Ab83sH4HN6DwLEUnQVLOYCup9dNIs04hji8VXCFfL+yZwGvB94HHglBrlJSIyK812dJKK2NlQq4tuf9XMfgA8yd1vrVViIiKSjujrWRRz97vyTkRERNI10/Us7iacqT0jd98714xERCQ5Mx1ZvL3o9kuAdwGDbLsG97HAV2NfyMwOAc4gDIqf6+6rpnjOq4DTCQsXbnb3V2bbfwv8GXgM2OruB8S+roiIVG+m61lcV7htZmcCr3X3jUXbrgSuImKQ28x2BM4EDgI2AOvM7Ap3v6XoOR3AWcAh7n6Xme1REubV7j53pxKIiMxhsSflLQAeLNn2IDD9/LTJuoFfu/ud7v4ocBlwWMlzjgC+VRgPcXddK0NEJBGxA9xXAFeY2WcIRwZ7AR/LtsdYCNxddH8D8NKS5+wLzMtmWu0KnOHuhW6uceBqMxsHvuzu50S+roiI5CC2WBwDfAo4m3CU8XvAgU9Htm+bYlvp4Hk78GLgQGAn4KfZare3Az3uvinrmrrGzG5z99Lra2BmRxOWIcHd6ezsnDah9vb2GR+PUW2MFHJIJUYKOaQSI4UcUopRbftUYsz1v4vY8yz+Cnw0+5mNwtFIwSJg0xTP2ezuDwEPZRdbej5wu7tvyvK418y+TejW2q5YZEcchaOO8ZnOlszjbMpqY6SQQyoxUsghlRgp5JBSjDzOek4hRiP/LmKvy7FgwYJpY8w0dfYVhW/vZrZkuue5e8wFkNYBXWa2D7ARWEYYoyh2OfBFM2sHnkDopjrNzJ4E7ODuf85uH0xY/VZEpCYWL17M2NjYpG3FH7gdHR2sX7++zlnNXrXX5YCZjyzOAv4uu33eNM8ZB55Z7kXcfauZHQt8lzB19nx3X29mx2SPn+3ut5rZVcAvCUuJnOvuN5vZM4Fvm1kh30vc/aqIfRMRmZWxsbFJH6ql3+pnWnuqWc00dfbvim7vU+0LufsaYE3JtrNL7p8MnFyy7U5Cd5SItIjiD+PC7VZd36naS7vmZVbLfYiI1FLhw3Cur9Sahzy6kPKg5T5ERKSs2OU+RESkhUUt9yEi0koOOOlaDrv4thkfbzXRYxZm9gLg5UAnRSfZufsn809LRKRxrj/xwPKzoY5srQH3qLWhsjOjh4ElwHLg74ETgGfXLjUREUlF7EKCJxJWg30j8HD271uALTXLTEREkhFbLPZw9x9ltx83sx3c/Urg0BrlJSIiCYktFhvM7BnZ7duBw8zs5cCjNclKRESSEjvAfRLwHOC3hHWZvkFYv6mvNmmJiEhKYledvbDo9pVm9mTgCe5eekEkEZEkaMmQfEUVCzM7HbjY3dcBZFe7UxeUiCRLS4bkK7Ybqg243MweAi4hrPz6P7VLS0REUhI1wO3uxxMuWPQBwkWMfmZmN5jZh2uZnIiIzN7ixYtZuHDhxA8w6f7ixYujY0Wfwe3ujwPXEC5rugK4gLCc+KkVZS8iInWR53U5KlnuYxegFzgceBVwHfCu6FcSEWkx030Yxwy0p3a1vtgB7q8DrwNuBC4F3uXuGjESkaY107fujo6OqBjFRaHS61GkdrW+2COL64ET3P2uWiYjIpKC0g/1Rl94aLbyXD039jyL1dERRURa3FRdSFD/bqQ8V8/VZVVFRHJW2oUEje9Gqlbs2lAiItLCdGQhIpKzcmMFhefMJSoWIiI5Kx0rgLl/tT11Q4mISFk6shARSVCe017zoGIhIpKgvKa95nFyIahYiIg0rTxPLtSYhYiIlKViISIiZdWtG8rMDgHOAHYEznX3VVM851XA6cA8YLO7vzK2rYiI1E5djizMbEfgTMLKtc8FDjez55Y8pwM4C3iDuy8G3hrbVkSkGRVfqGj+/PmT7lcyOJ2Heh1ZdAO/dvc7AczsMuAw4Jai5xwBfKuwsq2731tBWxGRpJRb+6nch31qK9/Wq1gsBO4uur8BeGnJc/YF5pnZD4BdgTPc/auRbUVEclf8gV+4HfOBPdVzGv1hX616FYu2KbaNl9xvB14MHAjsBPzUzH4W2RYAMzsaOBrA3ens7Jw2ofb29hkfj1FtjBRySCVGCjmkEiOFHFKJ0egcHnnkkYkYW7durSoPoOp9qbZ9NTHqVSw2AHsV3V8EbJriOZvd/SHgITP7IfD8yLYAuPs5wDnZ3fHiE1hKlZ7gMhvVxkghh1RipJBDKjFSyCGVGCnkkFcMoOoYtc5hwYIF0z5Wr2KxDugys32AjcAywhhFscuBL5pZO/AEQlfTacBtEW1FRKSG6jIbyt23AscC3wVuDZt8vZkdY2bHZM+5FbgK+CUwQpgie/N0beuRt4iIBHU7z8Ld1wBrSradXXL/ZODkmLYiIlI/OoNbRETK0kKCIiJNbrZTgIupWIiINLlCYahmVpeKhYhIjZSexT3bb/UpULEQEamRmS5eNNdogFtERMpSsRARkbJULEREpCwVCxGRJjc0NMSSJUvYaaedWLJkCUNDQxXH0AC3iEgTGxoaYvXq1QwMDLB06VLWrFlDf38/AL29vdFxdGQhItLEBgcHGRgYoKenh3nz5tHT08PAwACDg4MVxdGRhYhIwqo9+3p0dJTu7u5J27q7uxkdHa0oDx1ZiIgkbOPGjWzcuJFHHnlk4nYlurq6GBkZmbRtZGSErq6uiuKoWIiINLG+vj76+/sZHh5my5YtDA8P09/fT19fX0Vx1A0lItLECoPYK1asYNmyZXR1dbF8+fKKBrdBRxYiIjWVx7TVavX29rJ27Voefvhh1q5dW3GhAB1ZiIjUTF7TVlOgIwsRkRrJa9pqClQsRERqJK9pqylQsRARqZG8pq2mQMVCRKRG8pq2mgINcIuI1Ehe01ZToGIhIlJDvb299Pb26kp5IiLS/FQsRESkLBULEREpS8VCRETKUrEQEZGy2sbHxxudQ6007Y6JiNRQ21Qbm/nIom2mHzO7odxzah0jhRxSiZFCDqnESCGHVGKkkEMqMeqYw5SauViIiEhOVCxERKSsVi4W5yQQI4UcUomRQg6pxEghh1RipJBDKjEamkMzD3CLiEhOWvnIQkREIqlYiIhIWSoWIiJSloqFiIiU1TLXszCzNqAbWEg4u3sTMOLuUSP81bbPK8YMsfd399sqeP48d99Ssq3T3StecN/MPuDuZ1XYZm/gAXcfM7NnAAcAt7n7zZHtnwBsKbx3ZvZq4EXALe5+ZUU7sC3mbPYj999pvd/PWryXWZyK9sPMnufuv5zt62UxarIvFeZQ9X4UxToA2AvYCoxW8n88by0xG8rMDgbOAkaBjdnmRcCzgQ+4+9W1bJ9XjDLx73L3vSOe92rgP4D5wM+Bo939t9ljN7r7i8q0/3DJpjbgY8DnANz91IgcPgq8H3gEGAD6gWHgH4DzImP8AniVu//JzD4CvBFYA7wSuN7dP1aH/cjj76Lh72e172WO+/EY8BvgUuBSd7+lXJspYlS9L1mcaopvHvvxSuAUYAx4MeH3+WRgC/AOd787Ms5rgV4mf5m53N2vqjSnVjmyOAN4TeFDscDM9iH8IT2nxu1ziWFmg9M81AZ0ROQAcBLwWndfb2ZvAa4xs3e4+8+Y4VT/Ip/O8l1f9PwdgV0jXx/gHcBzgZ2B3wLPdPc/mtmTgP8Gyn6wADu6+5+y228DXu7uD5vZKuBGwgdVrfcjj7+LFN7Pat/LvPbjl9m+HA5cYWYPET5wLyt9j2dQ9b4UF18zKy6+nzazmC8zeezH6cDB2e9xH+BUd+8xs4OA84CDI/bjdGBf4KvAhmzzIqDPzF7n7sdH5gK0zphFO9verGIbgXl1aJ9XjKOAm4EbSn6uBx6NjPEEd18P4O7fIHzruMjM3kjc4ouLCR8CTwJOdvdPA39y909nt2M85u4PE741PQzcl+XzUGR7gAfM7O+y25uBJ2a324n7u85jP/L4nabwflb7XkI++zHu7je7+8fd/dnA+4A9gB+Z2U8iY+SxL4Xi2wOcRig47yV0N76nTvuxo7v/Mbt9F/B0AHe/hnCUEGOpuy9198vc/cfZz2XAPwFLI2NMaJUji/OBdWZ2GVA4fNsLWEao0rVun1eMdcDN7r7dH5yZfSoyxhYze6q73wOQHWEcCHwHeFa5xu5+F/AWMzuMcFRyWuTrFrvRzC4hfLBcSyhWVwFLgNhD9mOAi7Nuh3uB683sOuB5ZF0fddiPqn+nibyfVb2XOe7HpCNbdx8BRszsBOAVkTGq3hey4mtmj1JSfM0spn0e+3G9mZ1H+H0eBvwAwMx2JhTlGH81s+7s9Yu9BPhrZIwJLXFk4e6fB44k/BL/EXhZdvvI7LGats8rBvAW4KZp4u8TGeOjwJ4lbTcArwJWRcbA3S8nHAq/lKm/Xc/kn4H/JByavwM4m/Ce/A/h6Cnm9X9JGLi8lHB0dRbwXUK30CWxiWT7cRCz2I/s93YE1f1Oq86DKt/PvN7LLFY1+3HyNDHH3f26yNfPY18KxfdbbCu+R2Yf3jHFt+r9IHSD3UD4m/oe8JFs+zjw2sgY7wa+YGa3mNnV2c+twBeyxyrSEgPcInOJme3h7vc2Oo9WZWbtwFsJH8zfIBS+wwndQWdW2F3acGb2VELXVRuwodCrUKmWKBZmtjthYKsXeEq2+V7gcmCVu4/Vsn0zxShpvwfhP1Q1OcwqRpn4V7r768o856nAvwGPA58EjgPeBNwGHO/uv494nRsJ3z4vcfc7Z5nr30yx+UbghUCbu98fEWMX4ETgzYQBzEeBO4AvuftFs8mrKHbZ9zJ7XuG9uNTd75jla031O3kzcCvxv5NDCjN9sr+zUwndLjcDH3L3P8wmt0rkkYOZ7Ub4P7IIuLL4qMjMznL3D0TmktvU7lYZs3BgLWFK3T0w8Yf5buDrhMPmWravdYx31TFGLXN4d2wMM5tuim8b8IKIHC4E/ovQz/994GLg9YT+4bOzf8t5MmEW2g/M7B5C18f/c/dNEW0LNgO/K9m2kFAwxoFnRsS4GPg2oXvCCPt0GfAJM9vP3f/vTI1zeC9h23vx/SreiwvZ/nfyT1T2O/kcUJgWegrwe+BQwheBLxO+oMzIzF5CmDW4kfCBfT7hw36UMNX857XOAbgge71vAu8xszcDR7j7I4Qp0WXNNLXbwjkwFU3Xb5Vi8Qx3X128IfuQWmVmMX3k1bavdYzVZhYzSyOPGLXMoZL3Yh1wHVNP9+2IaL+nu38BJk4eK+TzBTN7b2QOf3L3fqDfzF5O6Kq4MesXvtTdY5aDPhF4DfARd/9Vls9vKhiDgvB+XpjdPtXM1rn7yuy9vAWYsVhQ/XsJ+bwXefxOih3g7i/Ibp9mZu+KbHcm4QinA/gJ4WjgoGwiyFmE8aBa5/Asd39zdnvIzD4OrDWzN1Tw2nlM7Z7QKsXid2Z2InBR4RDQzPYkfJONObml2vbNFCOFHCB0Tbzf3UdLHzCzmBjFkzu+OsNjUdz9R4SpkccRjozeRsS1A9x9IJtNdVqW979R+fXjHzKz/+PuPzazQ4H7s9iPZ90Q5VT7Xk4y2/eCfH4ne1g4QbAN2M3M2oq6XGJjzPPsbG8zW+1hijnufq2F8y7qkcN8M9vB3R/PXvuzZrYB+CGwS2SMPKZ2TwrWCt5GmAV0XfahBHAPcAXhsL3W7ZspRgo5AHyK6f/jHRfR/nIz28XdH3T3TxQ2mtmzgdsjc9juee7+GKELIvoM2Ww22luzD/prCCfXVeIY4Fwz2w/4Fdm5AGb2FMK35HI+RXXvJeTzXuTxO/kK204EvAjoBP6YdXPeFBnjr1kXzu7AuJn1uvuQhbOqH6tTDv9JmPr8vcIGd7/IzP5AmM0UI4/p+hNaYoBbZK4ws50IXRBRa2RJ/szs+YQxi8eBDwH/QhiT20gYsxhuYHoVMbPnAm+gaDYUcIXPYgmSlikWZrY/YYCseFbAFe5+az3aN1OMFHJopv2YIfZR7n5BBXksBP7b3R8s2j4xM2eGti8FbnX3B7Ji9VGyxfeAz7n7/9Y6hzzzmCLuV939nbNpO1tm1k04k3td9oF9CGFtqTX1jJGnljgpz8yWE2aHtAEjhAG9NuBSC+vA1LR9M8VIIYdm2o8yopbJMLM+wrTj44CbLZxFXRBz1vL5wF+y22cQul9WZ9tii1W1OeSVxxUlP/8JvKlwPzIPzGx/MzvQwrTk4u2HRLT9N2AQ+JKZfR74ImGc4aPZQHXM6+cRY3czW2Vmt5nZfdnPrdm2jpgYxVplzOK9wGLffknuUwkLn5U7c7na9s0UI4Uc8oiRQg6Y2XRLWbdRcqb9DN4HvNjdH7SwQuo3zOwZ7n4GcYtD7uDuW7PbB/i2lYd/bGY31SmHvPJYRDgSOZdwpNdGWDH2lMj2hcL3r4SB//PM7HgPZ6bD5Gmx03kLYcrxfMJY3KLsaOlkwsKOn41II48YeUzXn9ASRxaEvscFU2x/WvZYrds3U4wUcsgjRgo5QCgI7yTMwy/9uS8yxo6Fbp9smuSrgNdlRSvmg/pm2zZl+RcWrqGAme1LWBK7HjnklccBhGUyPg78r7v/AHjY3a/z+KU2CoWvl7AfK8yssEJrzL5sdffH3P0vwB3u/gCAh8UeY/8u8ojxDHdf7UVnbLv7Pe6+Cih7OYNSrXJk8UHgWjMbZdusgL0J1x04tg7tmylGCjnkESOFHCAs4LiLu99U+oCZ/SAyxj1m9oJCjOzb/esJ3Tp/H9H+n4EzzOwThJMEf2phyuzd2WP1yCGXPLKppqeZ2dezf/9A5Z9zkwqfmb2KcKT0dOKKxaNmtnP2Qf/iwkYLZ3PHftDnESOPKeoTWmmAewe2nfZemBWwLpveV/P2zRQjhRyaaT+qZWaLCN9Et1vzx8x6YmfvmNmuhDPG2wlrCEUvjZFXDtXmMUWs1wMv8zJnsZe0WQt8uLiAW1gv6nzCIpEzrvpqZvM9nGldur0TeJpnJ1/WIcaTCZMEDmNbl2Zhivpqj1hKpljLFAuYqKoTs1Yq/SOstn0zxUghhzxipJBDKjHyyGGauLsUz46qNzP7m0o+GPMsfLPNIWvTXhjDyQba9wfurDROXlqiWJjZCwhry+xO+ObXRhgIGyNc/vLGWrZvphgp5NBk+/FC4EtZjOL1e2abR2mMf/EyaxmVaR+VQ5n4sZf8fR7hTO+FwJXAcs+uemdmI+7eHRHjE+7+mez2c4Ehtl346G3u/t+z3hHiCp+Z9RAG2B8nnCD5GcK1YuYB5u4/jXiddxMG5e8DjiecXPkbwpXvTnT3SyPzzW1qd6uMWVxIWM5g0h+Kmf0DYUre82vcvplipJBDHjFSyIHsebXM48KIGDO1j8rBtr8Gd0Eb8ctTnEU4m/xnhDGKH5vZGzysYhu7PMWbCB/OEK4rcby7X2nhnIXTCdeHqMYtlB8cPg0wwn7/F9DrYSmWFxHOvu6JeJ0TgP0IZ4L/Anihu9+RHf1dQ1iocUYWpnYfTpjeXbgA0iLC1O7LsoHuaK0yG+pJU32j8HDd6SfVoX0zxUghhzxipJBDKjHyyOFzhJVndy352YX4z5ld3P0qdx9z9wHCJIGrsqI1my6QBZ6t8eThanE7xTQysw9P83MCcYVvnrv/KjuC+KO7/zjL4cbYHAhX69vs7r8BHswKJhV2Db4XeIm7r3L3r2U/qwhjbBUvzNgqRxZXmtl/ERYnK14j5Z3ErVtTbftmipFCDnnESCGHVGLkkcONwJC731D6gJnFzqhqM7PdPTtT292/b2Fp7m8CU133YyrPtHDyXRuwqGhGEcQfnXyOcFSydYrHYgpf8XM+VvLYEyJzuMvCyXi7AreZ2SmE64W8hrDkeYzC1O7SJfArmdo9oSXGLADM7HVs67srXiMl6tT5ats3U4wUctB+5Bsjh/b7Afe5++YpHtsz5huxmR1BGMD9Wcn2vYEV7v6+iBivLNl0g4dpvHsCb3H3sgsrmtlPgOOmKXx3u/teZdq/AfheUZEqbH8W8GZ3Pykih90IJwaOE87efi3hErm/Az7jkReCytpOObXbI5ZgKdYyxUJEJEYehS8VlufU7vHx8Zb+eetb33p0I9s3U4wUctB+6L2oZYwUcmjUfrTKAPdMYpciqFX7ZoqRQg55xEghh1RipJBDMjHM7OhG55BHDDP7TqVtWmWAe6b5xl+uR/tmipFCDtqPfGOkkENKMWYQ9UFt0yzXzvaDzdEsW2o9p/0oO/5TqiWOLKxJlrNOIUYKOWg/8o2RQg4pxSjj0Ygcql6u3XJaan06MQPkpVrlyKIplrNOJEYKOeQRI4UcUomRQg4pxZjJpyl/bY08lmvPY6n13QlTd3uBp2Sb7yUUslXuPhYbC1qnWFQ73ziP+crNEiOFHPKIkUIOqcRIIYdkYlj11xipdtVaCIXheMJS6x9x95vM7GGPX2Ydpr+exbuYxfUsWqVYfJDmWM46hRgp5JBHjBRySCVGCjmkFGNPwnkNfyrZ3gb8JKJ91cu1ez5LrT/D3VeXxL0HWG1m76kwVuucZ1HtfOM85is3S4wUctB+5BsjhRxSiWFm5wEXFJbpKHnsEnc/okz7Wqxa+09Aj1e21PrVwPeY+noWB7n7ayrJoWWKhYhIK7HJ17PYI9v8B8L1LFZ5tqJvLBULEZEWY2ZHuXu5gfpJWmLqrIiITPLpShu0ygC3iEhLyWFW1yQqFiIizanaWV2TqFiIiDSn7xAuKHVT6QNm9oNKg2mAW0REytIAt4iIlKViISIiZalYiIhIWRrgFqkTC9eY/jCwP/Bn4Cbgs1MtKyGSGh1ZiNSBmX0YOJ1wPYM9CYvbnUVYikEkeZoNJVJj2XUFNgJHufvXG52PyGzoyEKk9v4ReCLw7UYnIjJbKhYitfe3wGZ339roRERmS8VCpPbuAzrNTBNKZM5SsRCpvZ8CfyVcC1lkTtIAt0gdZLOhlgPvB64GtgCvAV7t7ic2MjeRGDqyEKkDdz+VcI7FJ4A/Eq4PfSww1MC0RKLpyEJERMrSkYWIiJSlYiEiImWpWIiISFkqFiIiUpaKhYiIlKViISIiZalYiIhIWSoWIiJSloqFiIiU9f8BN2PzfQEVNwAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now show a boxplot of the data across c\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.boxplot(accs)\n",
    "plt.xticks(range(1,len(costs)+1),['%.4f'%(c) for c in costs],rotation='vertical')\n",
    "plt.xlabel('C')\n",
    "plt.ylabel('validation accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "___\n",
    "# Next Time: Neural Networks\n",
    "___\n",
    "\n",
    "In this notebook you learned:\n",
    "- Formulation of Logistic regression with different optimization strategies\n",
    " - Line Search\n",
    " - Mini-batch\n",
    " - Stochastic Gradient\n",
    "- Newton's Approach using Hessian \n",
    "- Quasi Newton's Method\n",
    "- Use Exhaustive Searches for Finding \"C\" \n",
    "- And Training/Testing Splits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}