{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b910d945",
   "metadata": {},
   "source": [
    "# Logistic Regression\n",
    "\n",
    "Logistic regression problem taken from: Sören Laue, Matthias Mitterreiter, and Joachim Giesen. \"GENO--GENeric Optimization for Classical Machine Learning.\" Advances in Neural Information Processing Systems 32 (2019). and https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13b5ad66",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "Given a data matrix $X$ of dimension $n\\times d$, and a label vector $y\\in\\{-1,+1\\}^n$.\n",
    "\n",
    "We have the following unconstrained optimization problem,\n",
    "\n",
    "$$\\min_{w \\in R^{d}} ||w||_1 + C \\sum_{i=1}^n \\log(\\exp(-y_i(X_i^Tw))+1), $$\n",
    "\n",
    "where $C$ is the inverse of regularization parameter\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73483897",
   "metadata": {},
   "source": [
    "## Modules Importing\n",
    "Import all necessary modules and add PyGRANSO src folder to system path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ae68ad56",
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import datasets\n",
    "from sklearn.svm import l1_min_c\n",
    "import sys\n",
    "## Adding PyGRANSO directories. Should be modified by user\n",
    "sys.path.append('/home/buyun/Documents/GitHub/PyGRANSO')\n",
    "from pygranso.pygranso import pygranso\n",
    "from pygranso.pygransoStruct import pygransoStruct\n",
    "import torch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3713c13",
   "metadata": {},
   "source": [
    "## Data Initialization \n",
    "Specify torch device, and generate data\n",
    "\n",
    "Use GPU for this problem. If no cuda device available, please set *device = torch.device('cpu')*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f80d015b",
   "metadata": {},
   "outputs": [],
   "source": [
    "device = torch.device('cuda')\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "\n",
    "X = X[y != 2]\n",
    "y = y[y != 2]\n",
    "\n",
    "X /= X.max()  # Normalize X to speed-up convergence\n",
    "\n",
    "# Demo path functions\n",
    "cs = l1_min_c(X, y, loss=\"log\") * np.logspace(0, 7, 16)\n",
    "\n",
    "X = torch.from_numpy(X).to(device=device, dtype=torch.double)\n",
    "y = torch.from_numpy(y).to(device=device, dtype=torch.double)\n",
    "[n,d] = X.shape\n",
    "y = y.unsqueeze(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "174aa2e7",
   "metadata": {},
   "source": [
    "## Function Set-Up\n",
    "\n",
    "Encode the optimization variables, and objective and constraint functions.\n",
    "\n",
    "Note: please strictly follow the format of comb_fn, which will be used in the PyGRANSO main algortihm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "76877185",
   "metadata": {},
   "outputs": [],
   "source": [
    "# variables and corresponding dimensions.\n",
    "var_in = {\"w\": [d,1]}\n",
    "\n",
    "\n",
    "def user_fn(X_struct,X,y,C):\n",
    "    w = X_struct.w\n",
    "    \n",
    "    f = torch.sum(torch.log(torch.exp(-y* (X@w)) + 1))\n",
    "    f+= torch.norm(w,p=1)/C\n",
    "\n",
    "    # inequality constraint \n",
    "    ci = None\n",
    "\n",
    "    # equality constraint\n",
    "    ce = None\n",
    "\n",
    "    return [f,ci,ce]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b21c2ec",
   "metadata": {},
   "source": [
    "## User Options\n",
    "Specify user-defined options for PyGRANSO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "54137e9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "opts = pygransoStruct()\n",
    "opts.torch_device = device\n",
    "opts.maxit = 50\n",
    "opts.opt_tol = 1e-6\n",
    "np.random.seed(1)\n",
    "opts.x0 = torch.zeros(d,1).to(device=device, dtype=torch.double)\n",
    "opts.print_level = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be9ba1d7",
   "metadata": {},
   "source": [
    "## Main Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8ce3b204",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing regularization path ...\n",
      "Problem 1 with C = 0.10007147962830593 completed\n",
      "Problem 2 with C = 0.29307379488744323 completed\n",
      "Problem 3 with C = 0.8583089764312021 completed\n",
      "Problem 4 with C = 2.5136819185942896 completed\n",
      "Problem 5 with C = 7.361680888087905 completed\n",
      "Problem 6 with C = 21.55974671940413 completed\n",
      "Problem 7 with C = 63.14083504447964 completed\n",
      "Problem 8 with C = 184.917063358914 completed\n",
      "Problem 9 with C = 541.5563525125441 completed\n",
      "Problem 10 with C = 1586.0260682241312 completed\n",
      "Problem 11 with C = 4644.906624058538 completed\n",
      "Problem 12 with C = 13603.280537740799 completed\n",
      "Problem 13 with C = 39839.17360792678 completed\n",
      "Problem 14 with C = 116674.77924601595 completed\n",
      "Problem 15 with C = 341698.95806769416 completed\n",
      "Problem 16 with C = 1000714.7962830593 completed\n",
      "This took 14.768s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABAIklEQVR4nO3dd3iUVfbA8e9JI6EltNA7EkU6EWyICiqKFbuAoiiKqwKyuKL+lHVXWRULdoqIoKKI6IqoiKAggkrovaOUQEJJgJCEJHN+f8wbN4SETMK0JOfzPHkyeWfm3jMhnLlz3/veI6qKMcaY8iMk0AEYY4zxL0v8xhhTzljiN8aYcsYSvzHGlDOW+I0xppyxxG+MMeWMJX7jdyLyroj8Xwme10hEjopIqC/iClYi8q2I3BXoOE5FRH4SkXsDHYfxjCV+c0oiskNEenizTVV9QFX/Vdy+VfVPVa2sqjnF6U9E+otIjvOmcVhEVorI1SWJPRBU9UpV/cDb7YrIJBE57vxeDorIHBE504PnjRSRD70dj/EfS/ymvFisqpWBGOBt4BMRifF2J6Xw08iLzu+lAZAETApsOMYfLPGbEhGRCiLymojscb5eE5EKee5/TEQSnfvuFREVkRbOfZNE5N/O7Zoi8rWIpDijzp9FJEREpgCNgJnOiPQxEWnitBPmPLe6iLzv9HFIRL4sKm5VdQFTgErAGXley2gR+VNE9jlTUVHFeC3viMg3IpIGXCIi9UTkcxFJFpHtIvJInrY6i0iC88ljn4i84hyPFJEPReSA87tYIiK1nfv+mkZxfjdPicgfIpIkIpNFJNq5L/f3c5fzWvaLyJOe/Huq6jHgY6C109YYEdnpxLlURLo6x3sCTwC3Ov8uK/M001hEfhGRIyLyvYjU9KRv43+W+E1JPQmcC7QH2gGdgafgr+TwKNADaAFcfIp2hgG7gFpAbdxJRVW1H/AncI0zvfNiAc+dAlQEzgZigVeLCtoZkd8NZAF/OIf/A7R0XksLoD7wdDFeyx3Ac0AVYBEwE1jptNMdGCIiVziPHQOMUdWqQHNgmnP8LiAaaAjUAB4A0gvoq7/zdQnQDKgMvJnvMRcCcU7fT4vIWYX/RtxEpDLQB1juHFqC+/dRHfcbwmciEqmq3wHPA586/y7t8v0e7sb9bxEB/L2ofk1gWOI3JdUHeFZVk1Q1Gfgn0M+57xbgfVVd64wkR56inSygLtBYVbNU9Wf1YAMpEakLXAk8oKqHnOfOP8VTzhWRFCADGA30VdUkERFgIDBUVQ+q6hHcie22YryW/6rqL86niTZALVV9VlWPq+o2YHye9rKAFiJSU1WPquqveY7XAFqoao6qLlXVwwX01Qd4RVW3qepRYARwW+6nIMc/VTVdVVfifgNqV0A7uf7u/F624H4T6Q+gqh+q6gFVzVbVl4EKuN9MTuV9Vd2kqum439DaF/F4EyCW+E1J1eN/I2ac2/Xy3Lczz315b+f3Eu6k872IbBORxz3svyFwUFUPefj4X1U1BqgGfAV0dY7Xwv2pYakzxZICfOccB89eS95jjYF6uW057T2B+9MMwADcny42ONM5uSeZpwCzcZ972CMiL4pIeAF9FfR7D8vTPsDePLeP4U7ohRmtqjGqWkdVr1XVrQAi8ncRWS8iqc5riAaKmropTr8mgCzxm5LagzvJ5WrkHANIxH2yMFfDwhpR1SOqOkxVmwHXAo+KSPfcu0/R/06genFP0Dqj5EFAPxHpAOzHPaVytpMAY1Q12jnh6elryRvnTmB7nrZiVLWKql7l9L9ZVW/HPR3yAjBdRCo5n1j+qaqtgPOBq4E7C+iroN97NrCvOL+HU3Hm8x/D/WmnmvOGmQpIAa/XlEKW+I0nwp2Tj7lfYcBU4CkRqeWcxHsayF3iNw24W0TOEpGKQKFr9kXkahFp4Uy5pAI5gMu5ex/ueeyTqGoi8C3wtohUE5FwEbnIkxejqgeBCcDTzvTMeOBVEYl1YqqfZ07e49fi+B04IiL/EJEoEQkVkdYico7Tdl8RqeX0m+I8xyUil4hIG+ccxGHcUz+uAtqfCgwVkabOvHzufHu2J6/dQ1Vwv5kkA2Ei8jRQNc/9+4AmImL5o5SyfzjjiW9wj4pzv0YC/wYSgFXAamCZcwxV/RZ4HfgR9zRO7jx2ZgFtnwH8ABwFFgNvq+qPzn2jcL+5pIhIQScK++FOkBtwL0UcUozX9BpwlYi0Bf6RG6eIHHbiiSvBa8G5xuBq3PPb23F/opiAe6oEoCewVkSO4j7Re5szJ14HmI476a8H5uOe/slvonN8gdN+BvBwMV63J2bjnu7ahHsqKYMTp7M+c74fEJFlXu7b+IFYIRbja86qkjVABS+PTP2uLL0WU37ZiN/4hIjcIO718dVwz2XPLK2Jsiy9FmPAEr/xnftxT79sxT1vPyiw4ZyWsvRajLGpHmOMKW9sxG+MMeVMWNEPKRkRmYh7dUOSqrbOc/xh4G+4PzLPUtXHimqrZs2a2qRJE1+FaowxZdLSpUv3q2qt/Md9lvhx7/L3JjA594CIXAJcB7RT1czcddNFadKkCQkJCT4J0hhjyioR+aOg4z6b6lHVBcDBfIcHAf9R1UznMUm+6t8YY0zB/D3H3xLoKiK/icj83KsZCyIiA8W9fW1CcnKyH0M0xpiyzd+JPwz3Nq/nAsOBac6l+idR1XGqGq+q8bVqnTRFZYwxpoT8nfh3ATPU7Xfce5FYsQZjjPEjfyf+L3EXkEBEWuIu1rDfzzEYY0y55svlnFNxVyuqKSK7gGdwbzA1UUTWAMeBuzwpumGMMd705fLdvDR7I3tS0qkXE8XwK+K4vkP9QIflNz5L/M6e4wXp66s+jTGmKF8u382IGatJz8oBYHdKOiNmrAYocfIvbW8kvlzHb4wxQeel2Rv/Svq50rNyGDlzLWGhQmRYKFERoUSGhxAZHkpkeChReb5XCAshJOR/a1J88Ubia5b4jTHlyp6UgmrYQ8qxLB76eHmB9+VXISzE/eYQFkry0UxyXCfOWKdn5fDS7I2W+I0xJhjUiY4kMTXjpOO1q1bgwwFdyMhykZ6VQ0ZWzl/fM7JySD+eQ0a2y/meQ8Zx9/3TEnYV2M/ulHRcLj3h00GwsMRvjCk3VJXqlcJPSvxR4aGMuPIszqhdpdht/rLlALsL+RRx+WsLuP+iZlzXvj4RYcGzJ2bwRGKMMT42/udtrN1zhN4d6lE/JgoB6sdEMap3mxJPywy/Io6o8NATjkWGh9Dv3EaEh4YwfPoqur30IxN+3sbRzOCo32MjfmNMuZCw4yAvfLeRq9rU4eVb2lPIpgHFlvuGUdCqHlVl/qZk3p2/lX/PWs/rczdz53lN6H9BE2pWruCV/kuiVBRiiY+PV9ud0xhTUgeOZtLr9YVUCA9h5sMXUjUy3O8xrNiZwrs/bWX2ur1EhIZwS3xD7uvajEY1KvqsTxFZqqrx+Y/biN8YU6a5XMrQaSs5eOw4MwadH5CkD9C+YQzv9uvE1uSjjJu/jU+W/MlHv/1Br7b1eKBbM86uF+23WCzxG2PKtLd/2sKCTck8d0NrWtf3X3ItTPNalXnhprY8enlLJi7czke//cnMlXu4qGUtHujWjPOa1UBEfHpRmE31GGPKrMVbD9Bnwq9c064er93qvXl9b0pNz+LDX//g/V92sP9oJu0aRNOxcQyf/L6T9CzXX4+LCg8t9knowqZ6LPEbY8qkpCMZ9Hp9IVUjw/jqoQupVCG4JzgysnL4fNkuxi3Yxh8HjhX4mPoxUfzy+KUet1lY4rflnMaYMifHpQyeuoIjGVm83adT0Cd9gMjwUPp0acy8YRcX+pjCrjouLkv8xpgyZ8zczSzedoB/XdeauDrFvygrkEJDhPoxUQXeV6+Q48Vlid8YU6Ys2JTMG/M2c1OnBtwc3zDQ4ZRIQReFRYWHMvyKOK+0H/yff4wxxkN7UzMY8ukKWsZW4V/XtQ50OCV2qovCvMESvzGmTMjOcfHw1GVkZOXwVp+OREWEFv2kIHZ9h/o+293TEr8xpkwY/f0mluw4xJjb2tMitnKgwwlqPpvjF5GJIpLklFnMf98wEVERsULrxpjTNm/DPt6dv5U7ujTiuvbBuQd+MPHlyd1JQM/8B0WkIXA58KcP+zbGlBO7Dh1j6KcraVW3Kk9f3SrQ4ZQKvqy5u0BEmhRw16vAY8B/fdW3McZzpa1ebF7Hs1089PFyclzK2306Ehleuuf1/cWvc/wich2wW1VXFnXptIgMBAYCNGrUyA/RGVP+lMZ6sXm98N0GVuxM4e0+HWlSs1Kgwyk1/LaOX0QqAk8AT3vyeFUdp6rxqhpfq1Yt3wZnTDlVWOHxl2ZvDFBEnvtuzV7eW7id/uc34ao2dQMdTqnizwu4mgNNgZUisgNoACwTkTp+jMEYk0dhWwB4a2sAX/nzwDGGT19JuwbRjLjqzECHU+r4LfGr6mpVjVXVJqraBNgFdFTVvf6KwRhzosK2AIiKCCXl2HE/R+OZjKwcHvx4KQK8eUdHKoTZvH5x+XI551RgMRAnIrtEZICv+jLGlEzP1id/4A4LEdKP59D95fnMWLaLYNvB97lZ61mz+zAv39KehtV9V72qLPPlqp7bi7i/ia/6NsYU7WhmNt+sTqRO1QqEhAiJKRl/reqJq1OFETNW8+i0lUxfuot/X9+aZrUCd1FU7sqj3c4U1CVxtbisVe2AxVPa2ZW7xpRTo2dvZO/hDKY/cD6dGlc76f4Zg87n49//5IXvNtBzzM88dEkL7u/WzO9TK/lXHgEs3naAL5fvLhUrj4KR7c5pTDm0YmcKHyzeQd8ujQtM+gAhIULfcxsz99FuXN6qNq/M2cSVY37m120H/BrrC99tOGnlUUaWq1SsPApWlviNKWeyclyMmLGa2CoVGN6z6G1+Y6tG8uYdHZl09zlk5bi4bdyv/P2zlRxM893J37TMbL5auYcHpiwlMTWjwMcE+8qjYGZTPcaUMxMXbmd94mHe7duRqpHhHj/v4rhYvh/SjTfmbWbcgm3MXb+PJ3u14saO9b1Sy/ZoZjZz1+/jm9WJ/LQxmcxsF7FVKlApIpS04zknPd5bRUnKI0v8xpQjOw8e49UfNnFZq9pccXbxL6GJigjlsZ5ncl37+jzxxWr+/tlKpi/dyXM3tKF5CU7+HsnIYu76JGatTmT+pmSOZ7uoXbUCt3duxFVt6hLfuBpfrdxz0hy/N4uSlEdWbN2YckJVuev9JSzdcZA5j3Y77RGzy6V8mrCTUd+sJyPLxaCLmzPo4uZF7pdzOCOLH9bt45vVe1mw2Z3s61SN5Mo2dejVpi4dG1UjJOTETxCleT+hQCqs2LolfmPKif+u2M3gT1bwzDWtuPuCpl5rN/lIJv+etY7/rthD05qVuLJ1Hf67Ys8JSfqSM2OdZJ/Iz5v3czzHRd3oSK5sXZdebevQoeHJyd6cPkv8xpRjKceO0/3l+TSoFsWMBy8g1AdJ9ufNyQz5ZDkH0rJOOJ7blUuhXnQkV7apy1Vt6tKhYYwlex8rLPHbHL8x5cCobzaQkp7FlAFdfJL0AbqeUctZ439i4ncpVK4QypQBXWjfMMYrJ4LN6bHlnMaUcb9uO8CnCTu5t2tTWtWr6tO+Clt6mZaZQ4dG1SzpBwlL/MaUYRlZOTzxxWoaVo9iSPeWPu+vsBPGtvQyuFjiN6YMe/unrWxLTuPf17chKsL3Wy0MvyKOqHyremzpZfCxOX5jyqgtSUd456ctXNe+Ht1a+qeYUe4SS1t6Gdws8RtTBrlcyogZq6kYEcb/+bkA+fUd6luiD3KW+I3xoUBdePRpwk6W7DjEize2pWblCj7vz5QulviN8ZFAFTJPOpLB89+s59xm1bk5voHP+jGlly8rcE0UkSQRWZPn2EsiskFEVonIFyIS46v+jQm0QBUyf3bmOjKzXDx3QxtbPmkK5MtVPZOAnvmOzQFaq2pbYBMwwof9GxNQgShk/uOGJL5elchDl7Yo0aZppnzwWeJX1QXAwXzHvlfVbOfHXwH7HGrKrMLm1kNChP+u2I3L5d3tUo4dz+apL9fQIrYyD3Rr7tW2TdkSyHX89wDfBrB/Y3xmxrJdHEzLJP9ES0RoCLWrVGDwJyu45s2F/Lw52Wt9vjpnE7tT0hnVuw0RYXaJjilcQP46RORJIBv46BSPGSgiCSKSkJzsvf8cxvhSjkt5/pv1PDptJZ2b1uBf17emfkwUAtSPieLFm9qy8B+X8tqt7UlNz6Lfe7/T773fWLM79bT6XbM7lfcWbuf2zo04p0l177wYU2b5dHdOEWkCfK2qrfMc6w/cD3RX1WOetGO7c5rSIDU9i0emLmf+pmTuPK8x/3d1K8JDCx9bZWbn8OGvf/LmvM0cOpbFte3q8ffL42hUo2Kx+s3OcXHD24vYeziDHx7tRnSU51W1TNkWFLtzikhP4DGgm6dJ35jSYFvyUe6dnMCfB47x/A1tuKNLoyKfUyEslAEXNuXm+AaMm7+NCQu38e2aRPp0aczDl7aghofr7yct2sHq3am8eUcHS/rGIz4b8YvIVOBioCawD3gG9yqeCsAB52G/quoDRbVlI34TzH7amMTDU5cTHhrCO3060qVZjRK1s+9wBq/9sJlpCTuJCg9l4EXNuLdrUypGFD4+23XoGJe/uoAuTaszsf85tnzTnMAKsRjjZarKhJ+3M+rb9cTVqcr4OzvRoFrxpmkKsiXpKC/N3sDstfuoVaUCg7ufwa3nNDxp2khVGfBBAou3HmDOoxd5pW9TtgTFVI8xZUXudsczlu3mytZ1ePmWdqccmRdHi9jKjO0Xz9I/DvGfb9fz1JdrmLhwO8OviKOnU9bwpdkb2e1cD3B9+3qW9E2x2IjfmGLadziD+6csZcXOFIb2aMnDl7bwWQlBVWXu+iRe+G4Dm5OO0qh6FHsPZ3I82/XXYyLDQ/hP77a2MZo5SWEjflvsa0wxrNyZwrVvLmTTviO827cTg3uc4dO6sSJCj1a1+W7IRbx4Y1t2HUo/IekDZGS5fL4NhClbLPEb46Evlu/i5rGLCQ8NYcaD59OzdR2/9R0aItxyTkMK+4Duy20gTNljc/zGFCHHpbw4ewNj52+jS9PqvNO3E9UrRQQklnoxUX/N7ec/boynbMRvzCkczsji3g+WMHb+Nvqd25gP7+0SsKQPVtrQeIeN+I3JI2/hlFpVKqAoh9Ky+Pf1rel7buNAh2elDY1XWOI3xpG/cErSkUwAHrqkeVAk/VxW2tCcLpvqMcZRUOEUgC+W7wlANMb4jiV+YxyBKJxiTCBY4jcGOJqZXege9rZixpQ1lvhNuZeYms7N7y7meLaL8NATL8ayFTOmLLKTu6ZcW7M7lXsmLeHY8Rwm3dOZQ2nHbcWMKfMs8Zty64d1+3h46nKqV4rg80FdiKtTBcASvSnzLPGbckdVef+XHfxr1jra1I9mwl3xxFaJDHRYxviNJX5TrmTnuHj263VMXvwHV5xdm9du7UBURGjRTzSmDLHEb8qNo5nZPPzxMn7cmMzAi5rxeM8zfbqzpjHBymerekRkoogkiciaPMeqi8gcEdnsfK/mq/6NySt35c6Czft57obWPHHVWZb0Tbnly+Wck4Ce+Y49DsxV1TOAuc7PxvjUmt2pXPfmL+w8eIyJ/c+hT5fg2X7BmEDwWeJX1QXAwXyHrwM+cG5/AFzvq/6NAffKnZvfde+h//mg8+nWslagQzIm4Pw9x19bVROd23uB2oU9UEQGAgMBGjVq5IfQTFliK3eMKVzArtxVd7HfQgv+quo4VY1X1fhatWyUZjyXnePima/W8uzX67i8VW0+HXieJX1j8vAo8YvIBSJSybndV0ReEZGSTJTuE5G6Tjt1gaQStGFMoY5mZnPf5AQmL/6DgRc1450+nWy5pjH5eDrV8w7QTkTaAcOACcBkoFsx+/sKuAv4j/P9v8V8vjEnyFs4pXbVCoSIsO9IJs/d0NpO4hpTCE+nerKdqZnrgDdV9S2gyqmeICJTgcVAnIjsEpEBuBP+ZSKyGejh/GxMieQWTtmdko4Cew9nsic1g3u7NrWkb8wpeDriPyIiI4C+wEUiEgKEn+oJqnp7IXd1L0Z8xhSqsMIpX69MZMSVZwUgImNKB09H/LcCmcAAVd0LNABe8llUxnjACqcYUzKejviHquo/cn9Q1T9F5GwfxWRMkTKzc4gMDy1wxG+FU4w5NU9H/JcVcOxKbwZijKf2H82kz/jfSM/KISzECqcYU1ynHPGLyCDgQaCZiKzKc1cVYJEvAzOmIOsTD3PvBwkcSMvkrTs6kpXjssIpxhRTUVM9HwPfAqM4cV+dI6qafzsGY3xqzrp9DPlkOZUjw5h2/3m0bRADWOEUY4rrlIlfVVOBVOB2EQnFvcVCGFBZRCqr6p9+iNGUc6rK2AXbeOG7DbStH824O+OpXdWuxDWmpDw6uSsiDwEjgX2AyzmsQFvfhGWMW2Z2DiNmrGbGst1c3bYuo29uR2S4XYlrzOnwdFXPECBOVQ/4MBZjTrD/aCb3T1nK0j8O8ehlLXn40haI2B76xpwuTxP/TtxTPsb4Rd6TuG/36chVbeoGOiRjygxPE/824CcRmYX7Qi4AVPUVn0RlyrXv1+5lyKcrqBoZzmf3n0+bBtGBDsmYMsXTxP+n8xXhfBnjdarKO/O38tLsjbStH834O+OJtZO4xnidR4lfVf8JICIVVfWYb0My5VFGVg5PzFjNjOW7uaZdPV66qa2dxDXGRzzdj/88EVkHbHB+bicib/s0MlNuJB/J5I7xvzJj+W6GXdaS129rb0nfGB/ydKrnNeAK3Pvpo6orReQiXwVlyo91ew5z7wdLOHQsi3f6dORKO4lrjM95XHNXVXfmW0p38u5YxhQhb+GUapUiOJKRRc3KFfjsgfNoXd9O4hqTa9a2WYxZNoa9aXupU6kOgzsOplezXl5p2+PlnCJyPqAiEg4MBtZ7JQJTbuQWTsndUfNg2nFEYFC3Zpb0jclj1rZZjFw0koycDAAS0xIZuWgkgFeSv6e7cz4A/A2oD+wG2js/l4iIDBWRtSKyRkSmiogt3SgHCiqcogpjF2wPUETGBKcxy8b8lfRzZeRkMGbZGK+07+mqnv1AH290KCL1gUeAVqqaLiLTgNuASd5o3wQvK5xiTNGyXFkkpiUWeN/etL1e6aOobZkfU9UXReQN3HvznEBVHzmNfqNEJAuoCOwpYTumlFi7J5UQgZyT/oqscIoxufan7+fv8/9e6P11KtXxSj9Fjfhz5/ETvNIboKq7RWQ07gvC0oHvVfV7b7Vvgs93axIZ+ulKKkeGkZHlIjPb9dd9VjjFGLdVyasY+tNQDmce5uaWNzNz68wTpnsiQyMZ3HGwV/oqalvmmc73D7zSGyAi1YDrgKZACvCZiPRV1Q/zPW4gMBCgUaNG3ure+JGq8taPWxj9/SbaN4xh3J2dWLTlgBVOMSaf6Zum8/xvzxNbMZYpV03hzOpn0ql2J5+t6hHVAj5753+QyBzgZlVNcX6uBnyiqlcUu0ORm4GeqjrA+flO4FxVfbCw58THx2tCgtc+dBg/yMjK4R+fr+K/K/Zwfft6/OdGuxLXmPyO5xzn+d+e5/PNn3N+vfN5oesLxETGeK19EVmqqvH5j3u6nLNWbtIHUNVDIhJbwlj+BM4VkYq4p3q648WpJBN4SYczuG/KUlbuTGH4FXE8eHFz207ZmHz2pu1l2E/DWLV/FQNaD+DhDg8TGuKfwZGniT9HRBrlVtwSkcYUcLLXE6r6m4hMB5YB2cByYFxJ2jLBZ83uVO6bnEDKsSze7duJnq29czLKmLIkYW8Cw+YPIyM7g1cufoXLGl/m1/49TfxPAgtFZD4gQFec+feSUNVngGdK+nwTnL5dncjQaSuoXjGC6YPO4+x6dlGWMXmpKh9v+JjRS0bToEoDJl4xkeYxzf0eh6fr+L8TkY7Auc6hIc7afmNQVd6Yt4VX5myiQ6MYxvbrRGwVuybPmLzSs9N5dvGzfL3tay5ueDHPX/g8VSKqBCSWotbxn6mqG5ykD/9bb9/ImfpZ5tvwTLDLyMph+PRVzFy5hxs61GdU7zZ2EteYfHYd2cXQn4ay8eBG/tb+bwxsO5AQ8XTjBO8rasT/KO4pnZcLuE+BS70ekSk19h3OYODkBFbtTuWxnnEM6mYncY3Jb9HuRTz282O4XC7e7P4mFzUI/MbGRSX+Oc73Aaq6zdfBmNJj9S73SdzDGVmM7duJy8+2k7jG5KWqvLfmPd5Y/gbNopsx5pIxNKoaHNckFfVZY4TzfbqvAzGlx6xVidw8dhGhIcL0B863pG9MPmlZaQybP4wxy8ZwWePL+Oiqj4Im6UPRI/6DIvI90ExEvsp/p6pe65uwTDBSVV6fu4VXf9hEx0YxjO0XT60qFQIdljEBl3fv/JpRNRGE/Rn7GdZpGHedfVfQTYEWlfivAjoCUyh4nt+UYXmLptSNjiS2agVW7Eyld4f6PG8ncY0BTt47Pzk9GYABrQfQv3X/AEZWuKIS/3uq2k9ExqvqfL9EZIJC/qIpe1Iz2JOawdVt6/LyLe2CbgRjTKAUtHc+wDfbv2FIpyH+D8gDRc3xdxKRekAfEakmItXzfvkjQBMYBRVNAVj+Z4olfWPyKGyPfG/tne8LRY343wXmAs2Apbiv2s2lznFTBlnRFGM8U7tibfYeOznJe2vvfF845YhfVV9X1bOAiaraTFWb5vmypF9GuVxK5QoFjwmsaIox/6OqBSZ4b+6d7wseXTqmqoNE5EIRuRtARGqKSFPfhmYCIS0zm0EfLeVIZjah+aZ0rGiKMSf6csuXrEheQY9GPahbqS6CULdSXUaeP9Jre+f7gkd79YjIM0A8EAe8D0QAHwIX+C4042+7Dh3jvslL2bj3ME/1OosalSIY/f0mK5piTAG2pWxj1O+j6FynM6O7jfbblsre4OnunDcAHXBvpYyq7hGRwOwuZHxi6R8HuX/KUjKzXLzX/xwuiXOXW7ihY4MAR2ZM8MnIzmDY/GFEhUUxquuoUpX0wfPEf1xVVUQUQEQq+TAm42efJezkyS/WUC8mkk8GxtMi1t7TjTmVF5e8yJaULbzT4x1iK5a0JlXgeJr4p4nIWCBGRO4D7gHG+y4s4w85LuU/365n/M/bOb95Dd7u05GYihGBDsuYoDZ7x2w+2/QZd599NxfWvzDQ4ZSIp/vxjxaRy4DDuOf5n1bVOUU8zQSxwxlZDJ66nB83JnPneY35v6tbER4auG1ijSkNdh7ZychFI2lbqy0Pd3w40OGUmKcjfoBVQO7GLCtPp1MRiQEmAK1xXw9wj6ouPp02jef+OJDGgA8S2L4/jX9d35p+5zYOdEjGBL2snCyGzx+OiPDiRS8SHhIe6JBKzNNVPbcALwE/4b6I6w0RGa6qJd21cwzwnareJCIRQMUStmOKadHW/Tz4kbt+zpQBnTm/ec0AR2RM6fDastdYe2Atr178KvUrl+7VbcWpuXuOqiYBiEgt4AdKsF2ziEQDFwH9AVT1OHC8uO2Y4vvw1z8Y+dVamtSsxHt3xdO4hp2jN8YTC3YtYPK6ydwadys9GvcIdDinzdPEH5Kb9B0H8PDirwI0BZKB90WkHe6tIAaralreB4nIQJyC7o0aBc8+1qVRVo6Lf329jsmL/+CSuFqMub0DVSNL78dUY/xpb9penlz4JHHV4hh+zvBAh+MVnibv70Rktoj0F5H+wCzgmxL2GYZ7q+d3VLUDkAY8nv9BqjpOVeNVNb5WrVol7MqkHDtO//d/Z/LiPxh4UTMm3HWOJX1jPJTtyubxnx8nMyeTl7q9RIXQslF/oqhi6y2A2qo6XER6A7lrlxYDH5Wwz13ALlX9zfl5OgUkfnP6tiQd5d4PlrAnJYOXbmrLzfENAx2SMaXK2FVjWbpvKc9f+DxNo8vOLjVFTfW8hlN+UVVnADMARKSNc981xe1QVfeKyE4RiVPVjUB3YF1x2zEny1s4pXqlCI5mZFElKpyP7+tCfBPbRduUXnkrXNWpVIfBHQf7fC+c3xJ/Y+zKsVzb/FquaV7sVBfUikr8tVV1df6DqrpaRJqcRr8PAx85K3q2AXefRluGkwunHEg7jgAPXtzckr4p1fJXuEpMS2TkopEAPkv+B9IP8PjPj9O4amOe7PKkT/oIpKLm+GNOcV+J9+dV1RXO/H1bVb1eVQ+VtC3jVlDhFAXeW7gjIPEY4y0FVbjKyMlgzLIxPunPpS6eXPgkhzMPM7rbaCqGl73V5kUl/gRni4YTiMi9uFfjmCBhhVNMWeXvCleT1k7ilz2/8Ng5jxFXvWxuQ17UVM8Q4AsR6cP/En087m2Zb/BhXKYYNu49QkiIkOPSk+6zwimmtFJVvtzyZaH3h4WE8e32b+nRuIfXrqJdmbySN5a9wWWNL+OWuFu80mYwOmXiV9V9wPkicgnu7RUAZqnqPJ9HZjzyw7p9DP5kORXDQzieo2Rmu/66zwqnmNLqYMZB/rnon8zbOY+mVZuyJ20PmTmZf90fFhJGlfAqPLbgMepUqsMdZ97BjS1vpGpE1RL3mZqZymPzH6N2pdqMPH9kma4t7ekmbT8CP/o4FlMMqsrYBdt44bsNtK4Xzfg74/l124G/VvVY4RRTWs3fOZ+nFz3NkeNH+Hv83+nXqh/fbv/2pFU9Vza98q8ral9Z+grvrHyH3mf0ps+ZfWhYtXhLl1WVZxY9Q9KxJD648oPTegMpDUT15OmBYBMfH68JCQmBDiNoZGbnMGLGamYs202vtnUZfVM7oiJKVyEIY/I7lnWMlxJeYvqm6bSs1pJRXUfRslpLj567/sB6pqybwrc7viXHlcOljS6lX6t+dIzt6NHIfeqGqTz/2/MM6zSM/q37n+YrCR4islRV4086bom/dEk+ksn9UxJY9mcKQ3u05JHuLcr0R1JTPqxMXskTPz/BziM76d+6Pw+1f4iI0OLXhkg6lsQnGz5h2qZppGamcnaNs+nXqh+XN7m80PMAGw5u4I5Zd9Clbhfe6v4WIVJ2tie3xF8GrNtzmPsmJ3AgLZOXb25Pr7Z1Ax2SMacly5XF2JVjGb96PHUq1uG5C58jvs5JearY0rPTmbl1JlPWTWHH4R3EVozljjPv4KaWNxFdIfqEC8JCJISKYRX5uvfXVI8sW9e8WOIv5Wav3cvQT1dQNTKc8XfG06ZBdKBDMua0bE/dzoifR7D2wFqubX4tIzqPoHJEZa/24VIXC3cvZPLayfy29zeiwqJoV7Mdy5OXn3CyOCIkgmcveNbnVwP7myX+UkpVefunrbw0eyPtGsYwvl8nYqtGBjosY0pMVflk4ye8kvAKkWGRPH3e01zW+DKf97vx4EYmr5vMV1u/KvD+upXq8v1N3/s8Dn8qLPEXpwKX8bOMrBz+8fkq/rtiD9e1r8cLN7YlMtxO4prSK+lYEk//8jS/7PmFC+tfyLPnP0utiv7ZfTeuehzPXfgcM7fORDl5wOurC8KCkSX+IJV0OIP7pixl5c4Uhl8Rx4MXN7eTuKZU+37H9zz767NkZmfyVJenuCXuloD8TdepVIfEtMQCj5cXlviD0Jrdqdw3OYGUY1m827cTPVuXnz9IU/rl30lzYNuBLNu3jJnbZtK6Rmue7xrYLY4Hdxx8wqZvAJGhkQzuODhgMfmbJf4g883qRB6dtoLqFSOYPug8zq5nJ3FN6VHQTpr/XPxPBGFQu0Hc1/a+gBcpzz2B6+9tnoOJJf4goaq8PncLr/6wiY6NYhjbL55aVcpGtR9TfhS0kyZAjcgaPNj+wQBEVLBezXqVq0SfnyX+AMpbOCUyPIT0LBe9O9ZnVO82VAizk7im9CnsBOmBjAN+jsScStm5RK2UyS2csjslHQXSs1yEhQhdW9S0pG9KrcJOkJanE6elgSX+ACmocEq2Sxn9/aYARWTM6RvccTBhcuJEQnk7cVoaBCzxi0ioiCwXka8DFUMg7bbCKaYM6t6oOxVCKxAREoEg1K1Ul5HnjyzX8+nBKJBz/IOB9UDZ3v80H1Vlws/bC73fCqeY0uzzzZ+Tlp3G+1e875U9d4xvBGTELyINgF7AhED0HyhZOS6e+GINz32znnYNqhIZfuKv3wqnmNLseM5x3l/zPh1jO1rSD3KBmup5DXgMcBX2ABEZKCIJIpKQnJzst8B8JTU9i7vfX8LU3//kwYub88WDF/Kf3m2pHxOFAPVjohjVu40VTjGl1ldbv2LfsX3c3/b+QIdiiuD3qR4RuRpIUtWlInJxYY9T1XHAOHBv0uaf6Hxj58Fj3D1pCTv2p/HiTW25Jd5dHej6DvUt0ZsyIduVzYTVE2hdozXn1Tsv0OGYIgRijv8C4FoRuQqIBKqKyIeq2jcAsfjc0j8OMXByAtkuZfKAzpzfvGagQzLG677d/i27j+7m8c6P255SpYDfp3pUdYSqNlDVJsBtwLyymvS/WrmH28f/SuXIMGY8eL4lfVMm5bhyGLdqHHHV4ujWoFugwzEesCt3fUBVeXPeFl6es4lzmlRjbL94qlcqfhk5Y0qDOX/OYcfhHYzuNtpG+6VEQBO/qv4E/BTIGLwtMzuHEZ+vZsby3dzQoT7/udG2XzBll0tdjFs1jqbRTenRqEegwzEeshG/Fx1KO879U5by+46DPHpZSx6+1Aqhm7Jt/s75bD60mecvfJ7QEBvglBaW+L1ka/JRBkxawp7UDMbc1p7r2ttqHVO2qSpjV42lQeUGXNn0ykCHY4rBEr8XLN56gAc+XEpoiDD1vi50alw90CEZ43OL9ixi7YG1jDxvJGEhlkpKE/vXOk3TEnby5BeraVS9Iu/370yjGhUDHZIxPpc72q9dsTbXNr820OGYYrLdOUvI5VJe/G4Dj01fReem1Znx4AWW9I1fzNo2i8unX07bD9py+fTLmbVtlt9jSNiXwPKk5dzT+h7CQwNbUcsUn434PZS3aErd6EhqVYlg5a7D3HZOQ/51fWvCQ+091PheQaUNRy4aCeDXHTDHrRpHjcga9D6jt9/6NN5j2coD+Yum7EnNYOWuw1zbri6jerexpG/8pqDShhk5GYxZNsZvMaxMXsmvib/S/+z+RIZF+q1f4z2WsTxQUNEUgKV/pNhyTeNXhZU2LOy4L4xbNY6YCjHcEneL3/o03mWJ3wOFFUexoinG32pXrF3gcX+VNlx/YD0Ldi2gX6t+VAy3c1qllSV+D9SoXPB2C1Y0xfhbm5ptCjx+x1l3+KX/8avHUyW8Crefebtf+ivXVk2DV1vDyBj391XTvNa0Jf4i/LQxiZRjx8k/oWNFU4y/Hcw4yKLERbSq3oq6leoiCLFRsVQMrchnGz/jUMYhn/a/5dAW5vwxh9vPup0qEVV82le5t2oazHwEUncC6v4+8xGvJX9L/Kcwa1Ui901OoGXtqjx7/dlWNMUE1PhV40nPTmdU11F8f9P3rLprFXNvmcvYy8eyN20vj8x7hMycTN/1v3o8UWFR9D2rTG6mG1zmPgtZ+aaSs9Ldx73AlnMWYtqSnTw+YxUdGlVjYv9ziI4Kp9+5TQIdlimn9hzdw6cbP+W65tfRLKbZCfe1j23PqK6jGDZ/GE8ufJIXL3qREPHumO6Pw3/w3Y7vuKvVXVSLrObVtk0BUncV73gx2Yi/ABN+3sZjn6/ighY1mTKgM9FRdoGKCay3VryFIDzY/sEC77+8yeU82ulRZu+Y7ZOlne+tfo/wkHDuPPtOr7dtClChkKm06AZead5G/HmoKq/9sJkxczdzZes6vHZbe9tS2QTcxoMbmbl1Jv3P7n/K1Tv9z+7PriO7mLhmIg2qNODmljd7pf89R/cwc+tMbom7hZpRVkzI5/ZvhsyjIKGgeZaRh0dB96e90oWN+B0ul/Ls1+sYM3czN3VqwBu3d7Ckb4LC68tfp3JEZQa0GXDKx4kII7qM4ML6F/Lcr8+xcPdCr/Q/cc1EELi79d1eac+cgirMehQiq8CVL0J0Q0Dc3695Hdp659oJvyd+EWkoIj+KyDoRWSsig/0dQ37ZOS7+8fkq3v9lB/3Pb8KLN7YlzK7GNUFg6b6lLNi1gHta30N0hegiHx8WEsbobqM5o9oZDPtpGBsPbjyt/pOOJfHF5i+4rvl1frtWoFxb8zlsX+Ae2Xe+F4augZEp7u9eSvoQmBF/NjBMVVsB5wJ/E5FWAYgDcFfMenjqcj5buotHup/BM9e0IiTErsY1gaeqvLr0VWKjYulzVh+Pn1cpvBJvXvomlSMq8+DcB9mXtq/EMUxaO4kczSny04bxgoxUmP0E1OsAnXz76SoQxdYTVXWZc/sIsB4IyLrIY8ezuW/yUr5ds5enep3Fo5e1tC0YTNCYt3MeK5NXMqj9IKLCinexYO1KtXm7+9ukZaXxt7l/Iy0rrdj9H8w4yGcbP6NXs140rNKw2M83xTTv35CWDFe/Cj6uZhbQ+QwRaQJ0AH4r4L6BIpIgIgnJycle7/twRhZ3vvc7Czcn88KNbbi3a7Oin2SMn2S7snl92es0qdqE61tcX6I24qrH8XK3l9mSsoVh84eR7cou1vOnrJtCZk6mjfb9Yc9yWDIBzrnXPeL3sYAlfhGpDHwODFHVw/nvV9VxqhqvqvG1atXyat/7j2Zy+7hfWbkrhTdu78it5zTyavvGnK6ZW2eyLXUbj3R85LSqW11Q/wKeOvcpftn9C6N+G4WqevS81MxUpm6YyuVNLqdZtA2KfMqVA18PhYo14dKn/NJlQJZzikg47qT/karO8Gffianp9JnwG3tS0hl3ZzyXxMX6s3tjipSRncFbK96iTc029GjU47Tbu6nlTew6sov31rxHgyoNPFqd8/H6j0nLSuO+Nveddv+mCEvfd4/4e0+AyKJP4HuD3xO/uCfR3wPWq+or/ux7+/40+k74jdT0LCbf04XOTa02rgk+n2z4hH3H9jGq6yivnXN6pOMj7D66m1eWvkK9yvW4oskVhT726PGjfLj+Qy5ueDFx1W0/Kp86mgQ/PAtNu0Gbm/zWbSCmei4A+gGXisgK5+sqX3e6PvEwN7+7mPSsHKbed64lfROUDh8/zPjV47mg/gWcU+ccr7UbIiH8+8J/0yG2A0/8/AQrklYU+thPN37K4eOHub/t/V7r3xTi+6cgOx16vQx+XFji9xG/qi6Ekza79Lq8pRJrVq7AkYzjxFSswIf3dqFFrO0saILTxNUTOXz8MEM6DvF62xVCKzDmkjH0/aYvj8x7hA+v+pBGVU88v5Wenc7kdZO5oN4FtK7Z2usxmDy2L4BVn8JFw6HmGX7tukxepZS/VGLy0Uwys5X7LmpqSd8EraRjSXy0/iOuanoVZ1Y/0yd9VIusxts93kZRHpz7ICkZKSfcP33TdA5mHGRg24E+6d84so/DrGEQ0xi6DvN792Uy8RdUKlGBiQt3BCQeYzzxzsp3yNZsHurwkE/7aVy1Ma9f+jqJRxMZ/OPgv7ZyzszJZNKaScTXjqdj7Y4+jaHcW/wG7N8EV41278HjZ2VykzYrlWhKm+2p2/li8xfcGnerXy6W6hDbgecufI7hC4Zzz3f3kJyeTGJaIgDXNL/G5/2Xa4d2wPyX4KxroOXlAQmhTI74CyuJaKUSTbB6Y/kbVAit4Ncplp5Ne9KzcU9W7V/1V9IH+Gj9R8zaNstvcZQrqvDtP0BCoOd/AhZGmUz8w6+IIyr8xEuerVSiCVark1cz54853HX2XdSIquHXvlfuX3nSsYycDJ/s6W+ADbNg03dw8eNe21u/JMrkVE9uScTcVT31YqIYfkWclUo0QUdVeW3Za1SPrM5dZ9/l9/73pu0t1nFzGo6nuUf7sa3g3EEBDaVMJn5wJ39L9CbYLdqziN/3/s7jnR+nUnglv/dfp1KdE6Z58h43Xjb/BTi8C26aDaGBrepXJqd6jCkNXOri1aWvUr9yfW5p6b291otjcMfBRIZGnnAsMjSSwR0DXiajbNm3Dha/BR36QqNzAx1N2R3xGxPsvt3+LRsPbWRU11GEB2gE2KtZLwDGLBvD3rS91KlUh8EdB/913HiBqnvNfoUq0OPZQEcDWOI3JiCycrJ4c/mbxFWL46qmPt+x5JR6Netlid6XVk6FPxe5SydW8u/J+8LYVI8xAfDZps/YdXQXQzoNIUTsv2GZdeygez+eBp2hQ79AR/MX+4szxs/SstIYu2os59Q5hwvqXRDocIwvzf0npKc4VbWCJ93aVI8xfjZ57WQOZhzkjY5vWKnPsmznElg6Cc57COoE14Z3wfMWZEw5cCD9AJPWTqJHox60rdU20OEYX8nJhllDoUo998VaQcZG/Mb40fjV48nMyeSRjo8EOhTjS0vGw97VcPMH7tU8QcZG/Mb4yc4jO/l046dc3+J6mkY3DXQ4xlcOJ8K856BFD2h1XYmbSZ05k82Xdmf9Wa3YfGl3UmfO9FqINuI3xk/eWvEWoRLKoHaBvVz/JKumwdxnIXWXe/+Y7k9D28BcUFYmzB4BOcfhqpdKXFUrdeZMEv/vaTQjA4DsPXtI/L+nAYi+5vR3Tw1UsfWewBggFJigql7fpu6nF+4i/PPfiTkMKVUh68bOXPyPD4KmvdIQo71m77Z3VVWoe0U9aleqXeL2AFLfepKk92eQfVQJqyzE3t2b6L89V7LGVk0jdcwwkpZHkn2sDmEVM4ldM4zowZQ4+Xs1Ph+16dP2KuYQ2zOe6OrNCn2869gxspOTyU5KIjs5mSzne3ZSMtnJyRxLSIDs7BOeoxkZJL36mlcSv6jqaTdSrA5FQoFNwGXALmAJcLuqrivsOfHx8ZqQkOBxHz+9cBcxU36nQp7fW2YYHOoTT9dHxxU75p9fGUi1jxK81p4v2gz29vwZ48G+8Vw0LE97J/yNa6HHFrzyANU/WnZye3060nXoO2jucxH3bWc05/4/5NzOc/yXVwZRs4D2UvrE023YhGK/XoDUd0eyd+wXaM7/RpISqtQZcCXRdw91jzRzssGVBTlZ4Mp2f8857hw78b7Ut55k76IwNCckT3su6pyfTfTD/4HQCAgJh9Cw/90OCYeQMPd+M6F5b0eQ+v5o9o6fdXJ8919H9H1PgrryfGm+nwv6UlI/GMPeSXNPbvPuHkT3H+re4rjQLznpWOq4f7N37Jcnt/fADUQ/8M8S/Js8w9538/2bhCjRl3UhosMlfyXz7KSkvxK96+jRk9qRiAjCatUiLDaW9OXLC+5MhLPWF5oqC3i4LFXV+JOOByDxnweMVNUrnJ9HAKjqqMKeU9zE/0vns6h++HQjNaWVAtmhoAIu50vlFD+HuG/XSIWwAv47ZAscrAohCuJ8hSiEuJyfcd/Of3+oyw/FpU3Qy5vQw2Jj/3e7Vi3CYmsR7twOiY7+a3nv5ku7k71nz0lthdWrxxnz5nredyGJPxBTPfWBnXl+3gV0yf8gERkIDARo1KhR/rtPKaaQpK/AjsvqFastgCZz9hT4H7ik7fmizWBvzxdtFtYeQGLXOogCLnV/V3WSsoKTnFFFXP+7L3TZgQLbClXQFjXJEVARdzYX0BBBQ3D/4Pzsvk9Qgfo/nrzrZe7rrXVtyUob7v9qKQW/nSg1b+wKIaEgoe7RbUiIczs0z21n5Os8bv97Hxbe3j23u0fdrhxn9J2T53buz8535+f9n/9SeHvXd3Z+V5LvO4Ucd3/fP21e4W3efLHz6U0L+E6Bx/f/d0nh7ZXg32X/V8sKba/lr7+ekNA9FTt0yAlz/AASGUns0CHFjq8ggRjx3wT0VNV7nZ/7AV1UtdBCo94a8R+sChf8vr7YMXu7PV+0Gezt+aLN8tYewOb4s8g+eZaAsMpwRkLx29x84blk7089ub2a0Zyx8NeAx+eLNoO9vVypM2eS9OprZCcmEla3LrFDhxR7fr+wEX8glnPuBvIWFW3gHPOarBs7k5nvs0xmmPt4MLTnizaDvT1ftFne2gOIvbs3EnriYE1Cldi7e5esvX88iUScuDOoRIQT+48ngyI+X7QZ7O3lir7mGs6YN5ez1q/jjHlzvXJSN1cgRvxhuE/udsed8JcAd6jq2sKeU9wRPwT/ao/SEKO95uBrD3ywIsULI0tfxueLNoO9PW8JmpO7TjBXAa/hXs45UVVP+RsqSeI3xpjyLphO7qKq3wDfBKJvY4wp72zLBmOMKWcs8RtjTDljid8YY8oZS/zGGFPOBGRVT3GJSDLwRwmfXhPY78VwfK00xVuaYoXSFW9pihVKV7ylKVY4vXgbq2qt/AdLReI/HSKSUNBypmBVmuItTbFC6Yq3NMUKpSve0hQr+CZem+oxxphyxhK/McaUM+Uh8Zds8/jAKU3xlqZYoXTFW5pihdIVb2mKFXwQb5mf4zfGGHOi8jDiN8YYk4clfmOMKWfKReIXkZtFZK2IuEQkKJdxiUhPEdkoIltE5PFAx3MqIjJRRJJEZE2gYymKiDQUkR9FZJ3zNzA40DGdiohEisjvIrLSibf4RWD9TERCRWS5iHwd6FiKIiI7RGS1iKwQkaDe8ldEYkRkuohsEJH1TtlarygXiR9YA/QGFgQ6kII4BejfAq4EWgG3i0irwEZ1SpOAnoEOwkPZwDBVbQWcC/wtyH+3mcClqtoOaA/0FJFzAxtSkQYDJS815X+XqGr7UrCWfwzwnaqeCbTDi7/jcpH4VXW9qm4MdByn0BnYoqrbVPU48AlwXYBjKpSqLgAOBjoOT6hqoqouc24fwf2fp35goyqcuuUW8gt3voJ2BYaINAB6ARMCHUtZIiLRwEXAewCqelxVU7zVfrlI/KVAQQXogzY5lVYi0gToAPwW4FBOyZk6WQEkAXNUNZjjfQ14DHAFOA5PKfC9iCwVkYGBDuYUmgLJwPvONNoEEankrcbLTOIXkR9EZE0BX0E7cjb+IyKVgc+BIapaQAn04KGqOaraHnc96s4i0jrAIRVIRK4GklR1aaBjKYYLVbUj7mnVv4nIRYEOqBBhQEfgHVXtAKQBXjv3F5AKXL6gqj0CHcNp8HkB+vJMRMJxJ/2PVHVGoOPxlKqmiMiPuM+nBOOJ9AuAa51SqpFAVRH5UFX7BjiuQqnqbud7koh8gXuaNRjP/e0CduX5tDcdLyb+MjPiL+WWAGeISFMRiQBuA74KcExlgogI7nnS9ar6SqDjKYqI1BKRGOd2FHAZsCGgQRVCVUeoagNVbYL7b3ZeMCd9EakkIlVybwOXE5xvqKjqXmCniMQ5h7oD67zVfrlI/CJyg4jsAs4DZonI7EDHlJeqZgMPAbNxn3ycpqprAxtV4URkKrAYiBORXSIyINAxncIFQD/gUmcJ3wpnhBqs6gI/isgq3AOCOaoa9MskS4nawEIRWQn8DsxS1e8CHNOpPAx85PwttAee91bDtmWDMcaUM+VixG+MMeZ/LPEbY0w5Y4nfGGPKGUv8xhhTzljiN8aYcsYSvzGAiBwt+lGnfP50EWnm3K4sImNFZKuzNcBPItJFRCJEZIGIlJkLJ03pZInfmNMkImcDoaq6zTk0AfcmdmeoaifgbqCmswHfXODWwERqjJslfmPyELeXnH2eVovIrc7xEBF529kbfY6IfCMiNzlP6wP813lcc6AL8JSqugBUdbuqznIe+6XzeGMCxj5yGnOi3rivkmwH1ASWiMgC3FcAN8FdLyEW9xXWE53nXABMdW6fDaxQ1ZxC2l8DnOOLwI3xlI34jTnRhcBUZ4fMfcB83In6QuAzVXU5+6j8mOc5dXFvoVsk5w3heO6eMcYEgiV+Y05fOu7dKQHWAu2cqmqFqQBk+DwqYwphid+YE/0M3OoUQ6mFuwrS78AvwI3OXH9t4OI8z1kPtABQ1a1AAvBPZ2dQRKSJiPRybtcA9qtqlr9ekDH5WeI35kRfAKuAlcA84DFnaudz3HukrwM+BJYBqc5zZnHiG8G9uHeC3OIUpJ+Eu5oWwCXO440JGNud0xgPiUhlVT3qjNp/By5Q1b3Ovvk/Oj8XdlI3t40ZwOOquskPIRtTIFvVY4znvnaKpEQA/3I+CaCq6SLyDO46yX8W9mSnyM6XlvRNoNmI3xhjyhmb4zfGmHLGEr8xxpQzlviNMaacscRvjDHljCV+Y4wpZ/4fF/n7cV92vw0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Computing regularization path ...\")\n",
    "start = time()\n",
    "coefs_ = []\n",
    "i = 0\n",
    "for c in cs:\n",
    "    i += 1\n",
    "    comb_fn = lambda X_struct : user_fn(X_struct,X,y,c)\n",
    "    torch.autograd.set_detect_anomaly(True)\n",
    "    soln = pygranso(var_spec = var_in,combined_fn = comb_fn,user_opts = opts)\n",
    "    print(\"Problem {} with C = {} completed\".format(i,c))\n",
    "    arr = soln.final.x.T.tolist()\n",
    "    arr = np.array(arr).ravel()\n",
    "    coefs_.append(arr)\n",
    "print(\"This took %0.3fs\" % (time() - start))\n",
    "\n",
    "coefs_ = np.array(coefs_)\n",
    "plt.plot(np.log10(cs), coefs_, marker=\"o\")\n",
    "ymin, ymax = plt.ylim()\n",
    "plt.xlabel(\"log(C)\")\n",
    "plt.ylabel(\"Coefficients\")\n",
    "plt.title(\"Logistic Regression Path\")\n",
    "plt.axis(\"tight\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}