{
 "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 13.968s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABAdElEQVR4nO3dd3xUZdbA8d9JDy2hJBBqQJDeQkCsq2IXG669N3bdIvq6trUsii6WVWF3lXVtWNdeUFTsomJCQu8dSUIahAQICSSZ8/4xNxpChkzCtCTnu5/szNyZ+zxnJnJy57nPPY+oKsYYY1qOsGAHYIwxJrAs8RtjTAtjid8YY1oYS/zGGNPCWOI3xpgWxhK/Mca0MJb4TcCJyH9E5N5G7NdTRHaLSLg/4gpVIvKpiFwV7DgORkS+FZHrgx2H8Y4lfnNQIrJZRE7yZZuq+ntVndLQvlV1i6q2UdWqhvQnIleLSJXzR2OniCwRkfGNiT0YVPV0VX3J1+2KyEwR2ed8LkUi8oWIDPBiv8ki8qqv4zGBY4nftBQ/qWobIB54GnhDROJ93UkT/DbyqPO5dAcKgJnBDccEgiV+0ygiEi0i00Rkq/MzTUSiazx/u4jkOs9dLyIqIn2d52aKyIPO/U4i8rGIFDtHnd+LSJiIvAL0BD5yjkhvF5Fkp50IZ98OIvKi08cOEfmgvrhV1QW8ArQG+tV4L/8QkS0iku8MRcU24L3MEJFPRKQUOEFEuorIuyJSKCKbROSmGm2NEZFM55tHvog84WyPEZFXRWS781lkiEhn57lfhlGcz+YeEflZRApE5GURiXOeq/58rnLeyzYRudub36eq7gFeB4Y4bU0XkSwnzgUicqyz/TTgr8BFzu9lSY1meonIjyKyS0Q+F5FO3vRtAs8Sv2msu4GxwAhgODAGuAd+SQ7/B5wE9AWOP0g7twLZQALQGXdSUVW9AtgCnOUM7zxax76vAK2AwUAi8GR9QTtH5NcAFcDPzuaHgcOd99IX6Abc14D3cinwENAWmAd8BCxx2hkH3CwipzqvnQ5MV9V2wGHAW872q4A4oAfQEfg9UFZHX1c7PycAfYA2wL9rveYYoL/T930iMtDzJ+ImIm2Ay4BFzqYM3J9HB9x/EN4WkRhV/Qz4O/Cm83sZXutzuAb37yIK+Et9/ZrgsMRvGusy4AFVLVDVQuB+4ArnuQuBF1V1hXMkOfkg7VQASUAvVa1Q1e/ViwJSIpIEnA78XlV3OPt+d5BdxopIMVAO/AO4XFULRESAicAtqlqkqrtwJ7aLG/BePlTVH51vE0OBBFV9QFX3qepG4Nka7VUAfUWkk6ruVtW0Gts7An1VtUpVF6jqzjr6ugx4QlU3qupu4C7g4upvQY77VbVMVZfg/gM0vI52qv3F+VzW4/4jcjWAqr6qqttVtVJVHweicf8xOZgXVXWtqpbh/oM2op7XmyCxxG8aqyu/HjHj3O9a47msGs/VvF/bY7iTzucislFE7vSy/x5Akaru8PL1aaoaD7QHZgHHOtsTcH9rWOAMsRQDnznbwbv3UnNbL6BrdVtOe3/F/W0G4Drc3y5WO8M51SeZXwHm4D73sFVEHhWRyDr6qutzj6jRPkBejft7cCd0T/6hqvGq2kVVz1bVDQAi8hcRWSUiJc57iAPqG7ppSL8miCzxm8baijvJVevpbAPIxX2ysFoPT42o6i5VvVVV+wBnA/8nIuOqnz5I/1lAh4aeoHWOkm8ErhCRkcA23EMqg50EGK+qcc4JT2/fS804s4BNNdqKV9W2qnqG0/86Vb0E93DII8A7ItLa+cZyv6oOAo4CxgNX1tFXXZ97JZDfkM/hYJzx/Ntxf9tp7/zBLAGkjvdrmiBL/MYbkc7Jx+qfCOB/wD0ikuCcxLsPqJ7i9xZwjYgMFJFWgMc5+yIyXkT6OkMuJUAV4HKezsc9jn0AVc0FPgWeFpH2IhIpIsd582ZUtQh4DrjPGZ55FnhSRBKdmLrVGJP3+r045gO7ROQOEYkVkXARGSIio522LxeRBKffYmcfl4icICJDnXMQO3EP/bjqaP9/wC0i0tsZl68eb6/05r17qS3uPyaFQISI3Ae0q/F8PpAsIpY/mij7xRlvfIL7qLj6ZzLwIJAJLAWWAQudbajqp8A/gW9wD+NUj2PvraPtfsCXwG7gJ+BpVf3GeW4q7j8uxSJS14nCK3AnyNW4pyLe3ID3NA04Q0SGAXdUxykiO514+jfiveBcYzAe9/j2JtzfKJ7DPVQCcBqwQkR24z7Re7EzJt4FeAd30l8FfId7+Ke2F5ztc532y4E/N+B9e2MO7uGutbiHksrZfzjrbed2u4gs9HHfJgDEFmIx/ubMKlkORPv4yDTgmtN7MS2XHfEbvxCR88Q9P7497rHsj5pqomxO78UYsMRv/Od3uIdfNuAet78xuOEckub0XoyxoR5jjGlp7IjfGGNamIj6XxJ8nTp10uTk5GCHYYwxTcqCBQu2qWpC7e1NIvEnJyeTmZkZ7DCMMaZJEZGf69puQz3GGNPCWOI3xpgWxhK/Mca0MJb4jTGmhbHEb4wxLUyTmNVjjDEtzQeLcnhszhq2FpfRNT6W207tz7kju/mkbUv8xhgTYj5YlMNd7y2jrKIKgJziMu56bxmAT5K/DfUYY0yIeWzOml+SfrWyiioem7PGJ+1b4jfGmBCztbisQdsbyhK/McaEmK7xsQ3a3lCW+I0xJsRce0zyAdtiI8O57dT+PmnfEr8xxoSYXeXudX46t4tGgG7xsUydMNRm9RhjTHPkcinvLszmmL6dePX6I/zShx3xG2NMCEnfVERWURkXpHb3Wx+W+I0xJoS8syCbttERnDKoi9/6sMRvjDEhonRvJZ8uz2X88CRio8L91o/fEr+I9BCRb0RkpYisEJFJzvYOIvKFiKxzbtv7KwZjjGlKZi/LZc++Kn47yn/DPODfI/5K4FZVHQSMBf4oIoOAO4GvVLUf8JXz2BhjWrx3FmTTp1NrUnr693jYb4lfVXNVdaFzfxewCugGnAO85LzsJeBcf8VgjDFNxc/bS5m/qYjzR3VHRPzaV0DG+EUkGRgJpAOdVTXXeSoP6Oxhn4kikikimYWFhYEI0xhjgubdhTmECUxI8c1c/YPxe+IXkTbAu8DNqrqz5nOqqoDWtZ+q/ldVU1U1NSHhgEXijTGm2XC5lHcXZHNMvwSS4nxTluFg/Jr4RSQSd9J/TVXfczbni0iS83wSUODPGIwxJtSlbdxOTnGZ30/qVvPnrB4BngdWqeoTNZ6aBVzl3L8K+NBfMRhjTFPwzoJs2sZEcMqgOke+fc6fJRuOBq4AlonIYmfbX4GHgbdE5DrgZ+BCP8ZgjDEhbVd5BZ8sz2VCSndiIv03d78mvyV+Vf0B8HRqepy/+jXGmKbkk2W5lFe4uCBAwzxgV+4aY0xQvbMgm8MSWjOiR3zA+rTEb4wxQbJ5WykZm3fw21E9/D53vyZL/MYYEyTvLswO2Nz9mizxG2NMEFQ5c/ePOzyBzu1iAtq3JX5jjAmCnzZsZ2tJecDm7tdkid8YY4LgnQVZtIuJ4KSBgZm7X5MlfmOMCbCd5RV8ujyPs0d0Ddjc/Zos8RtjTIDNXprL3koXF4zqEZT+LfEbY0yAvbMgm36JbRjWPS4o/VviN8aYANpYuJsFP+/gtwGou++JJX5jjAmgdxZkEx4mnDcysHP3a7LEb4wxAVLlUt5bmMNvDk8gMcBz92uyxG+MMQHy4/pt5O0Mztz9mizxG2NMgLyzIJv4VpGMG5gY1Dgs8RtjTACUlFUwZ0Ue5wzvSnRE4Ofu12SJ3xhjAuDjpVvZW+nit0Gau1+TJX5jjAmAdxZk079zW4Z0axfsUCzxG2OMv60v2MWiLcVBnbtfkyV+Y4zxs3cW5BAeJpwbxLn7NVniN8YYP6pyKe8vyuaE/gkktI0OdjiAJX5jjPGr79cVkr9zb9Dn7tdkid8YY/zo7QXZtG8VyYkDAl933xNL/MYY4ycleyr4YkU+54zoRlRE6KTb0InEGGOamVlLt7KvyhVSwzxgid8YY/zmnQXZDOjSlsFdgz93vyZL/MYY4wfr8nexJKuYC1J7hMTc/Zos8RtjjB+8syCbiDDhnBFdgx3KASzxG2OMj1VWuXhvUQ4nDEikU5vQmLtfkyV+Y4zxse/XbaNwV2jN3a/JEr8xxvjY2wuy6Ng6ihMHBLfuvieW+I0xxod2lO7jy5UFnDOiG5HhoZliQzMqY4xpoj4K0bn7NVniN8YYH3o7M5tBSe0YFGJz92uyxG+MMT6yOm8ny3JKuCA1dI/2ASKCHYAxxgTaB4tyeGzOGrYWl9E1PpbbTu1/SLXyq9vLKS4DCNmx/WqW+I0xLcoHi3K4671llFVUAZBTXMad7y1FXcp5jRiXr90ewEOzV9EmOiJkFl6pTVQ12DHUKzU1VTMzM4MdhjGmGTj64a9/OTKvLSJMiIoIc/+Eh+13Pzpi/8eRzvNfrSrYL+lX6xYfy493ntjoOGdvnM30hdPJK82jS+suTEqZxJl9zmxQGyKyQFVTD3ifjY7KGGOaoK0ekj7AxOP6sK/Sxb4ql/u20sXeGvf3Vboor3Cxs6zyl9fVlfTr66c+szfOZvK8yZRXlQOQW5rL5HmTARqc/Ovit8QvIi8A44ECVR3ibJsM3AAUOi/7q6p+4q8YjDGmpvUFuxGBugY6usXHcvtpAxrcpqdvEF3jYxsTIgDTF07/JelXK68qZ/rC6T5J/P48AzETOK2O7U+q6gjnx5K+MSYgVuXu5KJnfqJVVDjRtRZFiY0M57ZT+zeq3dtO7U9sZLjP2gPIK81r0PaG8lviV9W5QJG/2jfGGG8tyy7hkmfTiAwPY9afjuGR84fRLT4WwX2kP3XC0EafiD13ZDemThjqs/YAurTu0qDtDRWMMf4/iciVQCZwq6ruqOtFIjIRmAjQs2fPAIZnjGlOFvy8g6tfmE9cq0hev34sPTu2ok9CG5/OuDl3ZDeftnfTyJu464e79tsWEx7DpJRJPmk/0JNNZwCHASOAXOBxTy9U1f+qaqqqpiYkJAQoPGNMc/LThu1c8Xw6HdtE8dbvjqRnx1bBDskrAzq4zzXERcchCEmtk5h81GSfjO9DgI/4VTW/+r6IPAt8HMj+jTEtx9y1hdzwciY9O7TiteuPILFdTLBD8lp6XjoAb45/k25tfH8tQECP+EUkqcbD84DlgezfGNMyfLkyn+tfyqRPQhvemDi2SSV9gLStafRo28MvSR/8O53zf8DxQCcRyQb+BhwvIiMABTYDv/NX/8aYlmn20lwmvbGIwV3b8dK1Y4hvFRXskBqk0lVJRn4Gp/c+3W99+C3xq+oldWx+3l/9GWPM+4uyufWtJaT0bM+L14ymbUxksENqsOXbllNaUcrYpLF+68Ou3DXGNAtvzN/CXe8v48g+HXn2ylRaRzfN9Jae6x7fH9NljN/6aJqfjDHG1PDSvM38bdYKju+fwH8uH0VMrQuqmpK03DQGdhhI+5j2fuvDEr8xLZyvSxQH2jPfbWDqp6s5eVBn/n3pSKIjmm7S31Oxh8WFi7l84OV+7ccSvzEtWF0liu96bxlAyCd/VeVfX6/niS/WMn5YEk9eNCLk6+DXZ1HBIipdlX4d3wdbgcuYFu2xOWsOqC5ZVlHFY3PWBCki76gqj81ZwxNfrOX8lO5Mv3hkk0/64B7miQiLYGTiSL/2Y0f8xrRgnkoH5xSXsa/SRVRE6CVTVWXKx6t44cdNXHpETx48ZwhhYRLssHwiPTedEQkjaBXp3yuMLfEb04K1bx1FUem+Op879tGvuebo3lx6RE/aBXlaZM3zELFR4ezZV8U1Rydz3/hBiDSPpL+jfAerilbxpxF/8ntfoffn3BgTEEuyitlZto/aeTM2MoyJx/Wmb2IbHv50NUdN/ZoHP155SAuLHIrq8xA5xWUosGdfFRFhwrBucc0m6QPMz5sPwBFJR/i9L6+O+EXkaGCxqpaKyOVACjBdVX/2a3TGGL/IKtrDdS9lkBQfyw3H9uGZ7zbWOatneU4Jz36/kRfnbWbmvM2MH5bEDcf1YXDXuIDF+shnqw84D1HpUv7x+VrOS2n4GrmhKi03jdaRrRnSaYjf+/J2qGcGMFxEhgO3As8BLwO/8Vdgxhj/KCmr4JqZGeyrdPHGxLH0TWzLlUcm1/naId3imH7xSG47tT8v/riZN+Zv4YPFWzmmbycmHteHY/t18vlRd8meCjJ/LmL+piLSNxWRW1Je5+uC9Q3EX9Jz0xndeTQRYf4fgfe2h0pVVRE5B/i3qj4vItf5MzBjjO/tq3Rx46sL+Hl7KS9fewR9E9t6tV/39q24d/wgbhrXj9fTt/Dij5u48oX5DOjSlonH9WH8sK6NPhFcsKucjE07mL9pO+mbiliTvwtViAoPY3iPONpER7B7b+UB+x3K0oahJmd3Dlm7srhs4GUB6c/bxL9LRO4CLgeOE5EwoOkVwTCmBVNV7npvGfM2bOfxC4Zz5GEdG9xGXGwkNx5/GNcek8ysxVt59vuN/N9bS3j0szVce0wyF4/pyderCg56QVj2jj3M31T0y8/GbaUAtIoKZ1Sv9pwxNIkxvTswokc8MZHhB1xrAIe+tGGoqS7TcEQX/4/vA4jWtepw7ReJdAEuBTJU9XsR6Qkcr6ov+ztAgNTUVM3MzAxEV8Y0W//8ah1PfLGWm0/qx80nHe6TNlWVb9cW8t/vNvLTxu1EhwtV6h6DrxYdEcbZI7pSUeli/qYitjpDN+1iIhjTu4Pz05HBXdt5nIvf1K8urs/t391ORn4GX1/wtU+HzkRkgaqm1t7u7RH/Lap6R/UDVd0iIoN9Fp0xxq8+WJTDE1+sZUJKNyaN6+ezdkWEE/onckL/RJZll3DhMz+xt9aJ2L2VLt7OzKZTm2iO6N2B3znJvn/ntl7Pv/f10oahxKUu0vPSObLrkQGbpeRt4j8ZuKPWttPr2GaMCTFpG7dz+ztLGdunAw9PGOa35DK0exzltZJ+NQEy7h7XrKZf+sq6HesoKi/ye5mGmg6a+EXkRuAPQB8RWVrjqbbAPH8GZow5dBsKd/O7VxbQo0Msz1ye6vcrcbvGx5JTx2ybrvGxlvQ9SMtNAwho4q/vv4LXgbOAWc5t9c8oVQ3M6WdjTKNs372Xa17MIDJcmHnNGOJa+X8+xm2n9ie2Vknk5nYi1tfSc9NJbpdMl9ZdAtbnQY/4VbUEKAEuEZFwoLOzTxsRaaOqWwIQozFNVrBOSpZXVHH9y5nk7yznjYlj6dHBv7VfqlW/t+Z8ItaXKqoqyMzP5OzDzg5ov95eufsnYDKQD7iczQoM809YxjR9wSp57HIpt7y5mMVZxcy4LIWRPf23oEddmvOJWF9btm0ZZZVlAR3mAe9P7t4M9FfV7X6MxZhm5WAlj/2ZGB/5bDWfLs/j7jMGctqQJL/1Yw5dWm4agjC6y+iA9uvtmZ4s3EM+xhgveSopsLW4DG+un2mM19J/5pm5G7libC+uP7a3X/owvpOem86gjoOIiw5c7SPw/oh/I/CtiMwG9lZvVNUn/BKVMU3cxsLdREWEsbfSdcBzCpz4+HdcNLoH56d0J6FttE/6/GZNAfd9uIIT+ifwt7OaT7ni5mpPxR6WFi7lqsFXBbxvbxP/FucnyvkxxtRhz75K/v31ep77fhOCEhkuVFT9enQfExnG+SndWJdfysOfruYfc9Zw0sDOXDymB8f2SyC8kQuKrNy6kz+9tpD+ndvy70tTiGgGq1E1d5n5mVRqZUDKMNfmVeJX1fsBRKSVqu7xb0jGND2qyqfL89x160vKmZDSjTtPH8C89ds9znBZX7CbNzO28O7CHD5bkUe3+FguTO3BhaO7kxTnfQGy3JIyrp2ZQduYSF64ejSto219paYgLTeNqLAovy+zWBdva/UcCTwPtFHVnk555t+p6h/8HSBYrR4T2tYX7GbyrBX8sH4bA5Pa8cA5gxmd3MHr/fdWVvHFynzemJ/FD+u3ESZwfP9ELh7dgxMHJB706H333kou+M9PZBXt4e3fH8nApHa+eEsmAM6fdT7to9vz3KnP+a2PQ63VMw04FfeFXKjqEhE5znfhGdP0lO6t5J9fr+OFHzYRExnO/WcP5rIjejZ4mCU6Ipzxw7oyflhXtmzfw5uZW3g7M5uJqwtIbBvNBanduSi1Jz07uufi17w2ICoijH2VLmZeO8aSfhOyvWw7a3esZVLKpKD07/V3QlXNqnWyqO6iHMY0c6rKx0tzeWj2KvJ2lnPBqO7ccfoAOrU59JO0PTu24rZTB3DLSYfzzZpC3pi/hRnfbuCpbzZwTN9OHJbQmjczsyivcJ803lvpIjJc2OFh3VwTmn5ZZjFAZZhr8zbxZ4nIUYCKSCQwCVjlv7CMCU3r8nfxt1krmLdhO4O7tuOpy1IY1cv3F0hFhIdx8qDOnDyoM7klZbydmc2bGe6hoNoqqtTv1wYY30rLTaNtZFsGdRwUlP69Tfy/B6YD3YAc4HPgj/4KyphQs3tvJdO/XMuLP26mdXQEU84dwqVjejZ6Fk5DJMXFctO4fvzxhL4c9tdP6nxNc1uGsDlTVdK2pjEmaQzhYeH17+AH3s7q2QZYUTbT7NWurfOXUw4nLEx4aPYqCnbt5eLRPbjt1P509MGwTkOFhwndDlL90jQN2buy2Vq6lauHXB20GOory3y7qj4qIv/Cfd3JflT1Jr9FZkyA1VVb59a3l+BSGNotjmeuGBXwuje13XZq/2a/DGFz91PuT0BgyzDXVt8Rf/U4vs2lNM1eXbV1XArxsZF88MejAzKsUx+rftn0peemk9gqkeR2yUGLob6yzB85ty8FJhxjgsfTOHlJWUVIJP1qVv2y6XKpi/l58zmu+3FBLanh1YRjEflCROJrPG4vInP8FpUxAZa9Y4/H1als/Nz4ypqiNRTvLQ7qMA94P6snQVWLqx+o6g4RSfRPSMYETmWVixd/3MwTX6zF5Tqwto6Nnxtfql5mMRj1eWryNvFXiUjP6hW3RKQXdZzsNaYpWZxVzF/fW8bK3J2MG5DI/ecMJnPzDhs/N36TnptOn7g+JLYK7nGzt4n/buAHEfkOEOBYYKLfojLGj3aVV/CPOWt4Oe1nEttGM+OyFE4b0gURoXv7VpbojV/sq9rHgvwFTOg3IdiheD2P/zMRSQGqB6Zudub2eyQiLwDjgQJVHeJs6wC8CSQDm4ELVXVH40I3pmFUlc+W5zH5oxUU7NrLlWN78ZdT+9M2xv+LkBuzpHAJ5VXlQR/mgXpO7orIAOc2BegJbHV+ejrbDmYmcFqtbXcCX6lqP+Ar57ExfpdTXMYNL2dy42sL6dA6mvf/cDT3nzPEkr4JmLTcNMIkLODLLNalviP+/8M9pPN4Hc8pcKKnHVV1rogk19p8DnC8c/8l4FvgDi/iNKZRKqtczJznPnmrCnefMZBrjk62hUpMwKXnpjOk0xDaRrX1boelb8FXD0BJNsR1h3H3wbALfRJLfYn/C+f2OlXd6IP+OqtqrnM/D+jsgzaNqdPS7GLuem8ZK7bu5MQBiTxwzmC6t28V7LBMC7Rr3y6Wb1vOtUOu9W6HpW/BRzdBhXNtSUmW+zH4JPnXd9hzl3P7ziH3VIu6V4DxODNIRCaKSKaIZBYWFvq6e9OM7SqvYPKsFZz71I8U7trLjMtSeP6qVEv6Jmgy8zKp0iqO7Hqkdzt89cCvSb9aRZl7uw/Ud8RfJCKfA31EZFbtJ1X17Ab2ly8iSaqaKyJJQIGnF6rqf4H/gnsFrgb2Y1qI/YuqxXDK4C58uiyP/F3lXOGcvG1n4/gmyNLz0okJj2F4wnDvdijJbtj2Bqov8Z8BpACvUPc4f0PNAq4CHnZuP/RBm6aFOrCoWjkv/riZpHbRvHfjUUEvqGZMtbStaaR0TiEqPMq7HeK6u4d36truA/UN9TyvqmnAs6r6Xe2fg+0oIv8DfgL6i0i2iFyHO+GfLCLrgJOcx8Y0Sl1F1QAkTCzpm5BRuKeQDSUbGjaNc9x9uC+ZqiEy1tl+6Oo74h8lIl2By0Tk2dqRqGqRpx1V9RIPT41rWIjG1M1TUbXc4vIAR2KMZ9VlGhpUn6f3cYBCTDyUlwR8Vs9/cM+37wMsYP/Er852YwJq++69PPzpao8zA6yomgkl6bnpxEXHMaDDAO93ys5w3172NvQY4/OY6ivL/E/gnyIyQ1Vv9HnvxjSAy6X8L2MLj362htK9lYwbkMiPG7b9svA4WFE1E1pUlbTcNMZ0GUOYNODakewMCI+CJC9PBjeQtyUbbhSRY4B+qvqiiHQC2qrqJr9EZUwty7JLuOfD5SzJKmZsnw5MOWcI/Tq3PWCpRCuqZkLJzzt/Jn9PfsPLMGdlQJdhEOGfJT69Svwi8jcgFegPvAhEAa8CR/slKmMcJWUVPP75Gl5J+5mOraOZdtEIzhnR9ZdFLGxREhPKGjW+X1UBWxfBqKv9ExTeV+c8DxgJLARQ1a0i4uV1x8Y0nKry/qIc/v7JKopK93HVkcnccvLhxMXanHzTdKTnptO1dVd6tO3h/U75K6CyDHr4r6aPt4l/n6qqiCiAiLT2W0SmxVubv4t7PljO/E1FjOgRz8xrxjCkW1ywwzKmQapcVaTnpXNSz5Matsxi9Ynd7sFP/G+JyDNAvIjcAFwLPOu3qEyLVLq3kulfreOFHzbRJiaCqROGclFqD8JCaL1bY7y1umg1u/btavj4fnYGtOkMcQ34ltBA3p7c/YeInAzsxD3Of5+qflHPbsZ4pbpO/gMfryS3pJyLUntwx+kD6NDay6scjQlBP+X+BMCYpAZOx8zOcB/t+3Exdm+P+AGWAtWnmJf4IRbTAm3eVsp9s1Ywd20hA5Pa8e9LUxjVy666NU1fWm4a/dr3o1NsJ+93Kt0GRRsh5Sr/BYb3s3ouBB7DXT9fgH+JyG2q6vOqnaZ5qzn9sk1MBHv2VhIbFcHfzhrEFWN7WZ180yyUV5azKH8RF/Zv4JW22ZnuWz+O70PD1twdraoFACKSAHyJH8o1m+ardlG1XeWVhItw+6mHc+VRvYMcnTG+s7hwMftc+7wvw1wtOwMkHLqO9E9gDm8Pr8Kqk75jewP2NQaAqZ+uOqCoWpUqz8y16wBN85Kem06ERDCq86iG7Zg9H7oMgSj/rh3h7RH/ZyIyB/if8/gi4BP/hGSam32VLp7/YRP5O/fW+bynYmvGNFVpW9MYmjCU1pENmPnuqoKchTDcU31L3zlo4heRvriXS7xNRCYAxzhP/QS85u/gTNM3b8M27vtwBesLdhMTEUZ5peuA11hRNdOclOwtYWXRSiYOm9iwHQtWwb7dfh/fh/qP+KfhLL+oqu8B7wGIyFDnubP8GJtpwgp2lvPQJ6v4cPFWenSI5YWrU9lZVrnfGD9YUTXT/GTmZeJSV+Pm7wN0T/V9ULXUl/g7q+qy2htVdZmIJPsnJNOUVVa5eCXtZ574fC17K13cNK4ffzj+MGIiw395jRVVM81ZWm4asRGxDOs0rGE7ZmdCq47Qwf/V7utL/PEHec6+n5v9LPh5B/d+sJyVuTs57vAE7j97ML077T/GaUXVTHOXlpvGqM6jiAxvYF2p7Pl+v3CrWn2JP1NEblDV/coziMj1uBdmMYai0n088ulq3szMoku7GJ6+LIXTh3RpWH0SY5qBvNI8Nu/czG8P/23DdizbAdvW+myFrfrUl/hvBt4Xkcv4NdGn4i7LfJ4f4zJNgMulvJmZxSOfrWZ3eSUTj+vDTeP60Sa6IReEG9N8pOemAw0swwyQ46TX7r5fbasu9a3AlQ8cJSInAEOczbNV9Wu/R2ZC2vKcEu75YDmLs4oZk9yBKecOoX8Xq9RtWrb03HQ6xHSgX/t+DdsxKwMQ6Jbil7hq87ZI2zfAN36OxTQBJWUVPOEsjNKhdRRPXDic80Z2s2Ed06LN3jib6Qunk1uaS0x4DJ9u+pQz+5zpfQPZGZA4CKIDc/Bk38mNR/svaxjD8QMSmLO8gO2le7libC9uPaW/LYxiWrzZG2czed5kyqvKASivKmfyvMkA3iV/lwtyMmFw4EbPLfGbOtWuq5NTXM5raVn06BDLrKuPYWh3WxjFGIDpC6f/kvSrlVeVM33hdO8S//Z1UF4SkAu3qlm9HVOnx+asOaCuDkCVSy3pG1NDXmleg7YfIAArbtVmid8cQFXJ8VA/J7e4vM7txrRUbaPqHpfv0rqLdw1kZ0BMHHRs4AnhQ2CJ3+xnY+Furnxhvsfnra6OMb96ZeUr7Ny3kzDZP5XGhMcwKWWSd41kZUC3VAgLXDq2xG8AKNtXxT/mrOG0ad+zeEsxE0Z2JSZy//88rK6OMb96acVLPJrxKCf1PIkHjnqApNZJCEJS6yQmHzXZu/H9vbugYGVAh3nATu62eKrKFyvzuf+jleQUl3HeyG7cdcYAEtvGcNzhOVZXx5g6PL/seaYtnMYpvU7h4eMeJjIsknP6ntPwhnIWAgo9LPGbANmyfQ9/m7Wcb9YUcnjnNrw5cSxH9On4y/NWV8eYAz279Fn+ueifnJZ8GlOPnUpE2CGk0WxnWLVbAxdsOUSW+Fug8ooq/vPdBp7+dgORYcLdZwzk6qOTibT1bo05qP8s+Q9PLX6KM3qfwUPHPHRoSR/cFTk79YfY9r4J0EuW+FuYb1YX8LdZK9hStIfxw5K458xBdImLCXZYxoS8GYtn8PSSpzmrz1lMOXoK4WHh9e90MKruGT2Hn+6bABvAEn8Lkb1jDw98tJLPV+bTJ6E1r153BMf06xTssIwJearKU4uf4pmlz3DOYedw/1H3H3rSByjaCHu2B2Thldos8TdzeyureO77Tfzr63UIwu2n9ef6Y/oQFWHDOsbUR1X516J/8eyyZzmv73lMPmryAVM3Gy07033bIzAVOWuyxN+M7F9bJ5azh3dlzoo8Nm4r5bTBXbj3rEF0s3n4pgmqLoKWV5pHl9ZdmJQyqWFF0BpBVXly4ZO8uPxFzu93PvcdeZ/vkj64T+xGtYGEAb5r00uW+JuJA2vrlDHjuw10bB3JzGtGc3z/xCBHaEzj1C6Cllua27AiaI2gqjye+TgvrXyJi/pfxF+P+Ktvkz64x/e7pYAvho0ayL7vNxOeautER4Rb0jdN2sGKoPmDqvJoxqO8tPIlLu5/MXcfcbfvk/6+PZC3PGALr9RmR/zNxFZPtXVKrLaOabpc6iK3NLfO57wugtYAqsojGY/w2qrXuGzgZdwx+g7/rDWxdRFoVcCv2K1mR/xNXMGucm55czHq4XmrrWOaIpe6+Hzz55w/63yPr1GUCbMmMGPxDNbtWIeqp38F3lFV/p7+d15b9RpXDLrCf0kfalTkDPyMHgjSEb+IbAZ2AVVApaoG5903YVUu5dW0n/nHnDWUV1Zx8sBEvl+/jfIK1y+vsdo6pqlRVb7O+pqnFz/N2h1r6R3Xm0v6X8L769/fb7gnOjyak3uezNbSrcxY4p5fn9wumXE9x3Fyr5MZ1HFQg5K2S108lPYQb619i2sGX8Mto27x76py2RnQoQ+0Ds6U6mAO9ZygqtuC2H+TtXDLDu79YDkrtu7k2H6duP/swfRJaHPArB6rrWOaClXl+5zv+feif7OqaBW92vVi6rFTOT35dMLDwhmeONzjrJ5tZdv4esvXfPHzF8xcMZPnlz9P19ZdGdfL/UdgeMLwg47Ru9TFAz89wLvr3uW6IdcxKWWSf5N+9YVbfY73Xx/1kEP9etSoTt1H/KneJv7U1FTNzMz0b1BNwI7SfTzy2WreyMiic7to7hs/mDOGdrH1bk2TparM2zqPpxY/xbJty+jepju/H/57zuxzZqPKIRSXF/Nt9rd8+fOXzNs6jwpXBQmxCZzY80RO6nUSqZ1TiQiL2G96aExEDGWVZdww9Ab+PPLP/v/3VLwFpg2FM/4BY27wa1cisqCuEZVgHfEr8LmIKPCMqv639gtEZCIwEaBnz54BDi+0uFzKW5lZPPzZanaVV3LDsb2ZdNLhtIm2c/OmaVJV0vPSeWrRUywuXEzX1l25/6j7Oeuws4gMa/w6zvEx8Zzb91zO7Xsuu/ftZm72XL7c8iWzNszizTVvEh8dT9+4vizZtoQKVwUAZZVlREgEfeL6BOYgyssVt0o++oiCJ6dRmZtLRFISibfcTNxZZ/kkhGAd8XdT1RwRSQS+AP6sqnM9vb4lH/Evzynh3g+Xs2hLMWOSOzDl3CH071L3ij/GNAUZeRk8tfgpFuQvoHOrzkwcNpHz+p5HZHjjE359yirLmJczjy+2fMEnGz9B65gOkdQ6ic9/+7nfYvjFp3fCgplwVxZ4eM8lH31E7r33oeW/nteQmBiSpjzQoOQfUkf8qprj3BaIyPvAGMBj4m+JdpZX8MTna3n5p810aB3FExcO57yR3WxYxzRZiwoW8dSip0jPSychNoG7xtzF+YefT3R4tN/7jo2IZVyvcYzrNY5PNn5S52v8MT20TtkZ0HWkx6QPUPDktP2SPoCWl1Pw5DSfHPUHPPGLSGsgTFV3OfdPAR4IdByhSlX5YHEOD81eTVHpXi4f24tbT+lPXKz/joaM8aXa5RXO63seSwqX8OPWH+kQ04HbR9/OBYdfQExEcKrCdmndpc5rA7xeI/dQVO6FvKUw9saDvyy37msXPG1vqGAc8XcG3neOXCOA11X1syDEEXS1Z+FcPrYn364pJH1TEcN7xDPzmtEM6RYX7DCN8Vpd5RWeXvI0seGx/N+o/+Oi/hfRKrJVUGOclDJpvxihgWvkHorcJVC1r97x/YiEBCoLCg7cnpTkkzACnvhVdSMwPND9hpq6aus88tkaYiPDmDphKBel9iAszIZ1TNNSV3kFgLjoOK4Zck0QIjpQ9TTQQBd9A7w+sRuZnHxA4peYGBJvudknYdi0kCDxVFsnvlUUl4xp2bOYTNPlaZw8f09+gCM5uDP7nBmYRF9b1nyI6wltPQ8rVRYWUr5oEbFjx1KxZYtfZvVY4g8ST7V18qy2jmnCgjp+3hRkZ9Zbf7/olVfRykq6Tv4bUcnJfgnDavUEWEWVixnfbrDaOqZZOrLrkQdsC9j4eajbuRV2Zh808VftLmXHG2/Q9uST/Zb0wY74A2rBz0X89b3lrMnfxdBu7ViXv5vySqutY5qHvNI8Pt/8Ob3b9aa8qjzw4+ehzovx/eJ33sa1cycdr7vWr6FY4g+Akj0VPPzZav43fwvd4mN57spUThrU2WrrmGZDVZn802SqtIqnTnqKHm17BDuk0JOdAeFR0GVonU9rRQVFL71Mq9GjiR3u3/kvlvj9SFX5cPFWHpy9kh17Krjh2N7cfNLhtHZKLZw7spsletMsfLD+A37M+ZG7xtxlSd+T7ExIGgERdV+wtvPTT6nMzaXL3+7zeyiW+P1k07ZS7v1gOT+s38aIHvG8fO1QBnVtF+ywjPG5vNI8Hs14lFGdR3HxgIuDHU5oqtznXnwl9bo6n1ZVtj/3PNH9+tLmuOP8Ho4lfh/bW1nFM99t5N/frCc6Iowp5w7h0jE9Cbc5+aYZUlXu/+l+Kl2VTDlqiu+XKGwu8pdDZbnHhVdKf/iBvWvXkjR1KhLm/8/QEr8P/bRhO3d/sIyNhaWMH5bEfeMHkdguOJelGxMIH274kB9yfuDOMXfSo50N8XiU7RSZ9DCjZ/vzLxDRuTNxZ54RkHAs8ftAUek+Hpq9incXZtOzQyteunYMvzk8IdhhGeNX+aX5PDr/UVISU7hkwCXBDie0Zc+HtknQ7sBzemXLlrMnLY3E225DoqICEo4l/kOgqry9IJu/f7KK3eWV/PGEw/jzif2IiQwPdmimGatdBC0Y0yVVlQfSHqDCVcGUo22Ip17ZGe5hnjqq625/4XnC2rQh/qILAxaOJX4v1Z56eeWRPflqdSHzNxUxOrk9D503lMM7W5184191FUGbPG8yQECT/0cbP2Ju9lzuGH0HPdtZiZGD2l0IOzbXeWJ3X1YWu+Z8TsfrriW8TZuAhWSJ3wt1FVSb+qm7oNoj5w/lglFWUM0ERl1F0Mqrypm+cHrAEn/BngIenv8wKYkpXDrw0oD02aQd5MKtohdnQng47S+/IqAhWeL3gqeCanGtorhotB3tmMDxVAQtUIuIVM/i2Ve1jweOfsCGeLyRnQFhEdB1xH6bK4uKKH7vPeLOPovIzokBDcl+a17wVFAt3wqqmQDrGNuxzu2BKoL28caPmZs9l5tG3kSvdr0C0meTl53hvlo3cv86XDteex0tL6fjtf4tz1AXS/z1WJZd4nEOvhVUM4G0p2KPx+faRbWjtKLUr/0X7ilk6vypjEwcyWUDL/NrX81GVSXkLDxgmMdVVsaO116jzQknEH3YYQEPyxK/BxVVLqZ9uZbznv6RVlHhRIXv/1FZQTUTaA/Pf5jtZdu5bsh1JLVOQhCSWidx7mHnsr54PVd8egVbd2/1S9+qygM/PeAe4jnqAcLDbOaaVwpXQUUpdN9//n7xe+9RVVxMx+vrvpLX32yMvw7rC3bxf28tYWl2CeeN7MbkswbzzZoCK6hmgmb2xtm8v/59bhh6Azel3MTNo27e7/nT+5zOX779C5fMvoR/nvhPhif4tsjX7E2z+Tb7W/6S+heS45J92nazljXffVvjil2trKToxZnEjhhBbEpKUMKyxF+Dy6W8OG8zj362mlZR4cy4LIXTh7rXuLSCaiZYtuzcwpS0KYxMHMkfRvyhztcc1fUoXj3jVf741R+59rNrefCYBzm99+k+6b9wTyFT06cyPGE4lw+83CdtthjZmdCqE7RP/mXTri++oCI7m8533oHUMa8/EGyox5FVtIdLn0tjyscrObZfJ+bcctwvSd+YYKmoquC2ubcRJmE8cuwjRIR5PlbrE9+H1898nSGdhnD73Nt5evHTqHpa8sc71RdqlVeWM+XoKTbE01DZGe4yDU6Cry7GFpWcTJsTTwxaWC0+8asqb2Vkcfr071mes5NHfzuMZ69MJbGt1dgxwTdt4TRWbl/JlKOnkNSm/gOR9jHtefaUZzn7sLOZsWQGd8y9g/LKxs8++2TTJ3yb9S1/Hvlnesf1bnQ7LdKeIti+br9hnj3p6ZSvWEGHa68JSDE2T1r0UE/BrnLuencZX60uYGyfDjz22+H06NAq2GEZA8Dc7Lm8vPJlLhlwCeN6jvN6v6jwKB48+kH6xPVh2sJp5OzOYfqJ0+kU26lB/W8r28bU+VMZljCMKwYF9gKjZiFngfu2xoye7c89T3inTsSdc06QgnJrsUf8nyzL5dQn5/LD+m3cO34Qr18/1pK+CRn5pfnc88M99G/fn1tTb23w/iLCdUOvY9rx01hXvI5LZ1/KmqI1Xu+vqkz5aQplFWU2xNNY2RkgYdDVfQK3fM0aSn/4gQ5XXEFYdN2LsQRKi0v8JXsquPmNRfzhtYX06NCK2Tcdw3XH9LaSCyZkVLmquOuHuyivKuex3zxGdHjjk8S4XuOYedpMqlxVXPnplXyX9Z1X+322+TO+zvqaP438E33i+jS6/xYtaz4kDoZodw2e7c8/j7RqRfuLLwpyYC0s8c9dW8ip0+by8dJcbjnpcN698Sj6JlphNRNanl32LBl5Gdx9xN0+GVcf1HEQr5/5Or3a9eLPX/+Zl1e8fNCTvtvKtvH39L8zrNMwrhx05SH33yK5XO6hHmd8v2LrVnbO/oT2F1xAeFxckINrxmP8NatpdomLoU+nVvy4oYh+iW149spUhnYP/odvTG0L8hcwY8kMxvcZz9mHne2zdju37szM02Zy9w9381jmY2zauYm/HvFXIsMi93udqvJg2oPsqdhjQzyHYtta2Lvzl4VXil56GUTocFVo/CFtlom/djXN3JJyckvKOaF/J2Zcnmr18k1IKi4v5o65d9C9TXfuGXuPz+d4t4psxePHP86/Fv2L55Y9R9bOLB4//nHion89CJqzeQ5fbfmKm1Nupk+8DfE0Wnb1hVujqSopYcfbb9PujNOJ7No1uHE5mmXi91RNc21+qSV9E5JUlXvn3cv28u28dsZrtI5s7Zd+wiSMSSmTSG6XzOSfJnP5J5dz/uHn8/qq13+p8Nm9TXeuGnyVX/pvMbIzICYeOhzGjv8+i+7ZQ8frglOeoS7NcozfUzVNT9uNCbbXV7/Ot1nfcuuoWxnUcZDf+zun7zk8d8pz5Jfm83jm4+SW5qLO/wrLCpmzeY7fY2jWsjOh+2hcFRUUvfoqrY89lpj+oVPbq1kmfk9VM62apglFq7av4vHMx/lN998EtOrlqM6jaBt14OSGvVV7mb5wesDiaHbKS6BgFXQfTcmHH1K1bVtIHe1DM038t53an9haQzpWTdOEotKKUm6bexvtY9oz5egpAa/dUlhWWOf2QC3s0izlLAQUTUqh6IUXiRk8mFZHjKl3t0Bqlon/3JHdmDphKN3iYxGgW3wsUycMtSJrJuT8Pf3vZO3K4uFjH6Z9TPuA9+9pAZdALezSLGVnAsKuDXvYt3kzHa+/LmjF2Dxplid3wappmtA3a8MsZm2YxR+G/4HRXQ5cjzUQJqVM2m/xdoCY8BgmpUwKSjzNQvZ8tNPhFL30PyJ79KDtyScHO6IDNMsjfmNC3eaSzTyY9iCpnVOZOGxi0OI4s8+ZTD5q8n4Lu0w+anLAFm5vdlQhO4OyysMpW7KEDldfhUSE3vF16EVkTDO3r2oft829jejwaKYeOzXoF0md2edMS/S+sn0DlO1g+5ISwuPjiZ8wIdgR1cmO+I0JsCcWPMHqotU8ePSDNpbe3GRnsLckgt0L19P+8ssJiw3NmYSW+I0JoG+2fMNrq17j8oGX85sevwl2OMbXsjPYvq49EhND+8suDXY0HlniNyZA8krzuHfevQzsMJBbRt0S7HB+tfQteHIITI533y59K9gRNVkVq9Mo2RhJ/IQJRLQP/CwtbwVljF9ETgOmA+HAc6r6sK/7+PaRq4h8dz7xO6G4HVScP4bj73gpZNprCjHae/Zte9vbwYjjw7ntnteICo9qdJslT91NwYvvUblbiWgjJF4zgbg/PtS4xpa+Rcn0WylYFEPlni5EtNpL4vJbiZsEDLsw+PH5qU1/tgcQGbG90W0FghzqmpwN7lAkHFgLnAxkAxnAJaq60tM+qampmpmZ6XUf3z5yFfGvzCe68tdteyOg+IrG/SP2dXtNIUZ7z6HXHrgTTO7T76JVv84Ll3Al6cbziJt4L2gVuCrBVfu2xv0arymZcim5PwhaFVajPRdJxypx974BYRE1fsJr3da+H0HJjMnkzqgjvj+cv39iVQV1/XqL58cl/32Q3P9+fGCbvzvL/Z4lzFnTVn69L2G1Hv+6r8fPsHaMh/o7aWR7viQiC1Q19YDtQUj8RwKTVfVU5/FdAKo61dM+DU38P44ZSIedB25XYF8jvuNEVUJdl180tj1/tBnq7fmjTU/tuYCyaFBx/7icF/1yv8Y2l/MagI47IcJ1YHuVYe6jdVHnhxr3FcKcf0JhzmPUPYYas7fusVQFwsIb9+9Oq8DTpyiNaNOdrBranueLkQ4eX4PD86LNhrxn8aK9hkbnub2INtAvc1XDG/QhT4k/GEM93YCsGo+zgSNqv0hEJgITAXr27NmgDuLrSPrVcsZ2aFBbAL1/KPJpe/5oM9Tb80ebntoToDClA6IKNZKxUOux4r7j3A9fUlJne+EuqEqOQ0WcrO9O3gioCEL1H5kazwv0nOf5/XY4YUCD3y/A9i9Xe3yu/UkjnCPcsFpHvDWOgkX2e03Re196bu/scfsfhde4L6J1POdi+5ylnt/ziYOcmGr+sXE+s+rHvxy5//r89k8XeI7xZPeyhqjzS/3lPlT/bt3/9+v9oi89Di7Q/viG/16Kvqr7d1I97BOKgnHE/1vgNFW93nl8BXCEqv7J0z6+OuIvagdHz2/4X2Bft+ePNkO9PX+02dLaA1iXOpDK3Qdub+zR5bpjxlK57cA/eBGd4uj3Q1rQ4/NHm6Heni95OuIPxqyeHKBHjcfdnW0+U3H+GPbW+i6zN8K9PRTa80ebod6eP9psae0BJF4z4YDhDQlXEq9p3IVCiXfcjUTtvwqXREWSeMfdIRGfP9oM9fYCIRhH/BG4T+6Ow53wM4BLVXWFp30aesQPoT/boynEaO859NoDP8xI+egjCp6cRmVuLhFJSSTecjNxZ50VMvH5o81Qb89XQubkrhPMGcA03NM5X1DVg35CjUn8xhjT0oXSyV1U9RPgk2D0bYwxLZ1duWuMMS2MJX5jjGlhLPEbY0wLY4nfGGNamKDM6mkoESkEfm7k7p2AbT4Mx9+aUrxNKVZoWvE2pVihacXblGKFQ4u3l6om1N7YJBL/oRCRzLqmM4WqphRvU4oVmla8TSlWaFrxNqVYwT/x2lCPMca0MJb4jTGmhWkJif+/wQ6ggZpSvE0pVmha8TalWKFpxduUYgU/xNvsx/iNMcbsryUc8RtjjKnBEr8xxrQwLSLxi8gFIrJCRFwiEpLTuETkNBFZIyLrReTOYMdzMCLygogUiMjyYMdSHxHpISLfiMhK57+BScGO6WBEJEZE5ovIEife+4MdU31EJFxEFonIx8GOpT4isllElonIYhEJ6ZK/IhIvIu+IyGoRWeUsW+sTLSLxA8uBCcDcYAdSF2cB+qeA04FBwCUiMii4UR3UTOC0YAfhpUrgVlUdBIwF/hjin+1e4ERVHQ6MAE4TkbHBDalek4DgLjXVMCeo6ogmMJd/OvCZqg4AhuPDz7hFJH5VXaWqa4Idx0GMAdar6kZV3Qe8AZwT5Jg8UtW5gOcFZUOIquaq6kLn/i7c/3i6BTcqz9SteiG/SOcnZGdgiEh34EzguWDH0pyISBxwHPA8gKruU9ViX7XfIhJ/E1DXAvQhm5yaKhFJBkYC6UEO5aCcoZPFQAHwhaqGcrzTgNsBV5Dj8JYCn4vIAhGZGOxgDqI3UAi86AyjPScirX3VeLNJ/CLypYgsr+MnZI+cTeCISBvgXeBmVa1jCfTQoapVqjoC93rUY0RkSJBDqpOIjAcKVHVBsGNpgGNUNQX3sOofReS4YAfkQQSQAsxQ1ZFAKeCzc39BWYHLH1T1pGDHcAj8vgB9SyYikbiT/muq+l6w4/GWqhaLyDe4z6eE4on0o4GznaVUY4B2IvKqql4e5Lg8UtUc57ZARN7HPcwaiuf+soHsGt/23sGHib/ZHPE3cRlAPxHpLSJRwMXArCDH1CyIiOAeJ12lqk8EO576iEiCiMQ792OBk4HVQQ3KA1W9S1W7q2oy7v9mvw7lpC8irUWkbfV94BRC8w8qqpoHZIlIf2fTOGClr9pvEYlfRM4TkWzgSGC2iMwJdkw1qWol8CdgDu6Tj2+p6orgRuWZiPwP+AnoLyLZInJdsGM6iKOBK4ATnSl8i50j1FCVBHwjIktxHxB8oaohP02yiegM/CAiS4D5wGxV/SzIMR3Mn4HXnP8WRgB/91XDVrLBGGNamBZxxG+MMeZXlviNMaaFscRvjDEtjCV+Y4xpYSzxG2NMC2OJ3xhARHbX/6qD7v+OiPRx7rcRkWdEZINTGuBbETlCRKJEZK6INJsLJ03TZInfmEMkIoOBcFXd6Gx6DncRu36qOgq4BujkFOD7CrgoOJEa42aJ35gaxO0xp87TMhG5yNkeJiJPO7XRvxCRT0Tkt85ulwEfOq87DDgCuEdVXQCquklVZzuv/cB5vTFBY185jdnfBNxXSQ4HOgEZIjIX9xXAybjXS0jEfYX1C84+RwP/c+4PBharapWH9pcDo/0RuDHesiN+Y/Z3DPA/p0JmPvAd7kR9DPC2qrqcOirf1NgnCXcJ3Xo5fxD2VdeMMSYYLPEbc+jKcFenBFgBDHdWVfMkGij3e1TGeGCJ35j9fQ9c5CyGkoB7FaT5wI/A+c5Yf2fg+Br7rAL6AqjqBiATuN+pDIqIJIvImc79jsA2Va0I1BsypjZL/Mbs731gKbAE+Bq43RnaeRd3jfSVwKvAQqDE2Wc2+/8huB53Jcj1zoL0M3GvpgVwgvN6Y4LGqnMa4yURaaOqu52j9vnA0aqa59TN/8Z57OmkbnUb7wF3quraAIRsTJ1sVo8x3vvYWSQlCpjifBNAVctE5G+410ne4mlnZ5GdDyzpm2CzI35jjGlhbIzfGGNaGEv8xhjTwljiN8aYFsYSvzHGtDCW+I0xpoX5f74ivhAs3hnOAAAAAElFTkSuQmCC\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
}