Comparison to regular sparse group LASSO

Compare directly to SGL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import regreg.api as rr\n",
    "%load_ext rpy2.ipython\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   user  system elapsed \n",
      "  2.371   0.031   2.424 \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "#install.packages('SGL', repos='http://cloud.r-project.org')\n",
    "library(SGL)\n",
    "load('/Users/jonathantaylor/Downloads/jon.RData')\n",
    "X = scale(x, TRUE, TRUE)\n",
    "Y = scale(y, TRUE, TRUE)\n",
    "X1 = X\n",
    "X2 = X1 / sqrt(499)\n",
    "#print(system.time(sgl1 <- SGL::SGL(list(x=X1,y=Y), g, standardize=TRUE)))   #1\n",
    "#print(system.time(sgl2 <- SGL::SGL(list(x=X2,y=Y), g, standardize=TRUE)))   #2\n",
    "print(system.time(sgl3 <- SGL::SGL(list(x=X1,y=Y), g, standardize=FALSE)))  #3\n",
    "#print(system.time(sgl4 <- SGL::SGL(list(x=X2,y=Y), g, standardize=FALSE)))  #4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [1] 0.11047717 0.09786806 0.08669806 0.07680293 0.06803717 0.06027187\n",
      " [7] 0.05339285 0.04729895 0.04190057 0.03711832 0.03288189 0.02912897\n",
      "[13] 0.02580439 0.02285925 0.02025025 0.01793903 0.01589159 0.01407783\n",
      "[19] 0.01247108 0.01104772\n"
     ]
    }
   ],
   "source": [
    "%%R -o X,Y,L,B,I\n",
    "L = sgl3$lambdas\n",
    "B = sgl3$beta\n",
    "I = sgl3$intercept\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, Y.shape\n",
    "X = X \n",
    "Y.shape = (-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13,\n",
       "       14, 15, 16, 17, 18, 19,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "groups = np.array(#[-1] + \n",
    "                  list(np.multiply.outer(np.ones(50, np.int), np.arange(20)).reshape(-1)))\n",
    "groups[:50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 0.35355339059327406,\n",
       " 1: 0.35355339059327406,\n",
       " 2: 0.35355339059327406,\n",
       " 3: 0.35355339059327406,\n",
       " 4: 0.35355339059327406,\n",
       " 5: 0.35355339059327406,\n",
       " 6: 0.35355339059327406,\n",
       " 7: 0.35355339059327406,\n",
       " 8: 0.35355339059327406,\n",
       " 9: 0.35355339059327406,\n",
       " 10: 0.35355339059327406,\n",
       " 11: 0.35355339059327406,\n",
       " 12: 0.35355339059327406,\n",
       " 13: 0.35355339059327406,\n",
       " 14: 0.35355339059327406,\n",
       " 15: 0.35355339059327406,\n",
       " 16: 0.35355339059327406,\n",
       " 17: 0.35355339059327406,\n",
       " 18: 0.35355339059327406,\n",
       " 19: 0.35355339059327406}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = 0.95\n",
    "lasso_weights = np.ones_like(groups) * alpha\n",
    "#lasso_weights[0] = 0 # unpenalized intercept\n",
    "group_weights = dict(#[(-1,0)] + \n",
    "                          [(j, (1 - alpha) * np.sqrt(50)) for j in range(20)]) # no group penalty on [-1] group -- intercept\n",
    "\n",
    "group_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = X.shape\n",
    "penalty = rr.sparse_group_lasso(groups, lasso_weights, weights=group_weights, lagrange=1.)\n",
    "alpha = 0.95\n",
    "Y = Y.reshape(-1)\n",
    "loss = rr.squared_error(X, Y, coef=1./n)\n",
    "lagrange_val = L "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "penalty = rr.sparse_group_lasso(groups,\n",
    "                                lasso_weights, \n",
    "                                weights=group_weights,\n",
    "                                lagrange=1)\n",
    "problem = rr.simple_problem(loss, penalty)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Without warm start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.34 s ± 268 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "soln = np.zeros(X.shape[1])\n",
    "for lagrange in lagrange_val:\n",
    "    penalty.lagrange = lagrange\n",
    "    problem.coefs[:] = 0\n",
    "    if hasattr(problem, 'final_step'):\n",
    "        step = problem.final_step\n",
    "    else:\n",
    "        step = 1.\n",
    "    soln = problem.solve(tol=1.e-10, start_step=step)\n",
    "    soln[:] = problem.coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With warm start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = rr.simple_problem(loss, penalty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "673 ms ± 144 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "soln = np.zeros(X.shape[1])\n",
    "for lagrange in lagrange_val:\n",
    "    penalty.lagrange = lagrange\n",
    "    problem.coefs[:] = soln\n",
    "    if hasattr(problem, 'final_step'):\n",
    "        step = problem.final_step\n",
    "    else:\n",
    "        step = 1.\n",
    "    soln = problem.solve(tol=1.e-10, start_step=step)\n",
    "    soln[:] = problem.coefs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "solns = []\n",
    "soln = np.zeros(X.shape[1])\n",
    "for lagrange in lagrange_val:\n",
    "    penalty.lagrange = lagrange\n",
    "    problem.coefs[:] = soln\n",
    "    if hasattr(problem, 'final_step'):\n",
    "        step = problem.final_step\n",
    "    else:\n",
    "        step = 1.\n",
    "    soln = problem.solve(tol=1.e-12, min_its=200)#, start_step=10*step)\n",
    "    solns.append(soln.copy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.1104770302772522, 0.11047714948654175, 0.11047716953570587)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = X.T.dot(Y) / n\n",
    "dual = penalty.conjugate\n",
    "lambda_max = dual.seminorm(G, lagrange=1)\n",
    "lambda_max, dual.terms(G).max(), L.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1366ec828>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8lfX5//HXdUZ2mAmyAglTpoyAIC4cgLWC4gCt2xa0Yv1q1Vqr1h92OCqtrQu0KNYquEXFgiioCCgBBGUKIUCYIRAgCSRnXL8/zuFwEhmBjJPkXM/H4zxy7s89uG7QNzefz33fH1FVjDHGRAdHpAswxhhTcyz0jTEmiljoG2NMFLHQN8aYKGKhb4wxUcRC3xhjooiFvjHGRBELfWOMiSIW+sYYE0VckS6gvJSUFE1PT490GcYYU6csXrx4l6qmHm+7Whf66enpZGVlRboMY4ypU0RkY0W2s+4dY4yJIhb6xhgTRSz0jTEmitS6Pv0j8Xg85ObmcvDgwUiXUmvExcXRunVr3G53pEsxxtQhFQp9ERkGPA04gZdU9bFy628EngS2BJueUdWXgutuAB4Mtv9JVaecaJG5ubkkJyeTnp6OiJzo7vWOqpKfn09ubi4ZGRmRLscYU4ccN/RFxAk8C1wI5AKLRGS6qq4st+k0VR1Xbt8mwB+BTECBxcF995xIkQcPHrTADyMiNG3alLy8vEiXYky18fuVUp8/8PH68fj8eLxKqc9HqTewLtDmpyT481BbqddPqU8PtwV/lvr8eH2KXxXVwAWUwuFlgm0KqsF2OOq2/uAXRfH7Az8D+wGh74eP4Q9OWhVqP7Svcvg4wV/rzbEDcTmrvge+Ilf6/YF1qpoNICJTgRFA+dA/kqHAp6q6O7jvp8Aw4I0TLdQCvyz7/TC1gaqyv8TL3mIPew942Hcg8DP8UxD8WVziDQawlgnyMqHs9ePxBdb7/FU7q5/gJwYvbrw4UCTsE1imTNvhdkUksP+hbcrvL4Cj3HrKHyO4DT85BjjlcPsWTWGjNqeKTz+kIqHfCtgctpwLnH6E7S4XkbOBtcBdqrr5KPu2OslajTHVQFUpLPEGArr4yMEd/tkXFub7DniOGU4uvDSkiIZSRBIHcOMlRrw0xEsMHtzBEI4RLzFhbTHixe061BbcBg9uOVKbj9jwY+ENbucps60rGKy13XPe4TzhHU3g3wdVr6oGcj8E3lDVEhEZC0wBzqvoziIyBhgD0KZNmyoqqer9+c9/5vXXX8fpdOJwOJg4cSJ9+/bl4Ycf5q233iIxMRGAK6+8kj/84Q8AJCUlUVhYGMmyTRTbe8BDzq4icvKL2LCriJxdReQXlZYJ9n0Hvce8qg4P7oYU0ViKaHto2RH42YjCwLIU0SBs20QpqcGzrR+Ew11A1aEiob8FSAtbbs3hAVsAVDU/bPEl4Imwfc8tt+/c8r+Aqk4CJgFkZmYe8VTT7/+4AqVWXs5jFx+xfcGCBXz00UcsWbKE2NhYdu3aRWlpKQ8++CDbt2/n+++/Jy4ujv379/PUU0/VSK3GABSXesnZVRwI9WC4hwf8IQ78tCCfVNlLo2BwN5AiGkkhDV1Fh4M9GNgNgj+TpP7cNVeibjw48Yd1wPh/0qkTiF1/mU4aQfUo7cc8RgWPHTy+H2GjngJENvQXAR1FJINAiI8GrgnfQERaqOq24OJwYFXw+0zgLyLSOLg8BPh9pauOgG3btpGSkkJsbCwAKSkpFBcX8+KLL5KTk0NcXBwAycnJPPLIIxGs1NRHJV4fm/KPFOzFbN8XHsrKKewhw7GdC2U7Ga5tZMh20mU7bWUnseKp0bq96mAfCezVRAqJp4QYPOqilMMfD65gm/sobU5KcePBRakG1gW2O1qbk1INrCs5dCxceHECdWcsLGLdO6rqFZFxBALcCUxW1RUiMh7IUtXpwG9EZDjgBXYDNwb33S0ijxL4iwNg/KFB3bpmyJAhjB8/nk6dOnHBBRcwatQoGjduTJs2bUhOTo50eaYe8Pr85O45wIb8IjbklQ33rQUHwvrOlSbsJ0O2caZjO+mu7WTINjJkB21le5V3qfhU2EciezWRAhLZp4nsDS4f/plUdjn4vZB4KhO0MS4HMU4HbqcQ43LgdjpCbYeW3U4h2eUkJmybn24nxDiduF2C2+HA4ZDAoKuAQyQ4UBu4QSLwXXAIoe8iwXWAw/HTttBxyhwj7Nhltg0em7B1we/htcS6nJX/wzuCCvXpq+oMYEa5tofDvv+eo1zBq+pkYHIlaqwVkpKSWLx4MV999RVz5sxh1KhRPPDAA2W2efnll3n66afJz89n/vz5pKWlHeVoJtr5/cqq7ftYsD6fbzfsZt3OQjbtLsYb1rfegCLSZTt9ZRtXOANX6xnBTwMpPqlfN08bsE2bUqBJoSAPD/ACTWJvuWAvJA49xsP7cW4HDePdoU/L+Bi6hC03jHfRMMFNcqw7EMTBYI4NC/FAKDvKBLvLIXaXWjWoE0/k1hZOp5Nzzz2Xc889lx49ejBx4kQ2bdrE/v37SU5O5qabbuKmm26ie/fu+Hy+SJdrahFVZX1eIfPX5zN/XT4LN+RTUOwBlAzZThfZyDDZTjv3oXDfRlPZf1K/VoEmkqPNydYW5PibB78HfhaScMR9ygd3i3g3ncsEt5tGCWWXGwR/VtcVqakeFvoVtGbNGhwOBx07dgTgu+++o3PnzvTu3Ztx48YxceJE4uLi8Pl8lJaWHudopr5TVTbtLmbB+nzmr89nQXY+efsD3S6tJY8LHSsY6F7JGY4VNJcTelYRgEKNI0fDAj0U7i0o4HB3o0OgdeMEMlISyUxJJL1pAq0bJ5QJ8AbxbuLcFtzRos6E/tHuqqkphYWF3HHHHRQUFOByuejQoQOTJk2iYcOGPPTQQ3Tv3p3k5GTi4+O54YYbaNmyJQDFxcW0bt06dJy7776bu+++O1KnYarR1oIDoZBfmJ3PloIDADRjD2c4VnCGKxDyaY6KPUldom42BIM9PNw3aHPyaER4X3nLhnGkpyRycUoiGSmJpDdNJCM1kbTGCcS47L2K5jDR6rov6CRlZmZq+UlUVq1aRZcuXSJUUe1lvy+Rlbe/hAXZ+SxYn8+C9bvIyQ/0szdhHwMcKxnoCIR8e8e2Yx5nnyawxN+RddqSHA2E+gZ/C7bRpExfempyLBlNE0lPSSAjJYmMlATSUxJp2ySR+Bi7Uo92IrJYVTOPt12dudI3JtIKiktZmL2bBet3sSA7n7U7Ag/dNaCI/o7VXO9awUDHSro4Nh3zOEUayyL/qSzwd2W+vxsrNB1/MNwbJbhJb5rIgJRE0oNX7RkpibRtmkBynL1R1VSehb4xx7Bm+37eXZrLvB93sXLbPlQhgYP0c6xhpGslAx0r6C4bcMrR/8Vcom4W+zsy39+NBf6uLNf2eHAR43TQK70Rv2nflAHtmtL5lGQaJ8bU4NmZaGShb0w5hSVePlq2lWlZm1m6qQA3XjIda7jLGbiS7yXrccvR787yqJPvtD0L/F1Z4O/GEn9HSojB6RB6pjVkTPumDGyXQt+2ja1bxtQ4C31jCNxts3RzAdO+3cyHy7dSXOqjk2zmIddcLnN+RRM5+vuTfCr8oBksCF7JL/J3ppg4RKBbywZc364pZ7RPITO9sXXRmIiz0DdRbXdRKe8uyWXaos38uLOQRA5wiXMho2Lm0Mex7qj7rfK3CfXJf+s/lX0EXrbX6ZQkrmqfwoB2TRnQrgmNEqy7xtQuFvom6vj9yrx1u5i2aDOzVm7H4/PTW9bxmGsOlzgXHPE1Blu0KZ/7erPA35WF/q7spgEAGSmJXNyuKWcE++VTk2Nr+nSMOSEW+iZqbC04wFtZubyZtZktBQdozD6ud85jVMwcOjm2/GT7UnUyy5/JNN9gvvZ3x4+DVo3iOa99Uwa2a8rA9k1p2Sg+AmdizMmz0D8JgenUFIfj5B568Xq9uFz2W18TSr1+Plu1g6mLNvPlj3mgfs50/MDv3XMY4sgi5ggDsmv9rZjmG8x7vjPZTQNaNYrnzsw0RvRqSdumCfY+GFOnWfJUUE5ODkOHDuX0009n8eLF3HfffbzwwguUlJTQvn17Xn75ZZKSkpgxYwZ33303iYmJDBo0iOzsbD766CMeeeQR1q9fT3Z2Nm3atOG1117j/vvvZ+7cuZSUlHD77bczduxY/H4/48aN4/PPPyctLQ23283NN9/MFVdcEenfgjpl3c79TFu0mXeXbCG/qJQW5HOH8wuucs2ltez6yfZFGsuHvoG86TuXJdoRt9PBkJ7NGZWZxqAOKTgdFvSmfqh7of9Iw2o89t5jrv7xxx+ZMmUKHTp0YOTIkcyePZvExEQef/xxJkyYwH333cfYsWP58ssvycjI4Oqrry6z/8qVK5k3bx7x8fGhVzgsWrSIkpISBg0axJAhQ1i8eDE5OTmsXLmSnTt30qVLF26++ebqO+d65KDHx4fLtjJt0WayNu7BjZfzHUsY7Z7D2Y7lOI5wL/1Sfwem+gbzkW8ARcTTsVkSD/ZLY2Sf1jSxe+ZNPVT3Qj+C2rZty4ABA/joo49YuXIlgwYNAqC0tJSBAweyevVq2rVrR0ZGBgBXX301kyZNCu0/fPhw4uMDfcCzZs1i+fLlvP322wDs3buXH3/8kXnz5nHllVficDho3rw5gwcPruGzrHt8fuW9pVuYMGsNW/cepL1s4QHXXEY6vyJF9v1k+z2axLu+s5jmO5e1mkZCjJNLerdkVP80eqc1su4bU69Z6J+AQ3PgqioXXnghb7zxRpn13333XYX2P3SMf/3rXwwdOrTMNjNmzCi/mzkKVWXOmp08/ska1uzYT09ZzxPuqZzpXHHE7b/09eBN37nM8mdSipvebRrxeL80Lu7ZkqRY+1/BRIe691/6cbpgasKAAQO4/fbbWbduHR06dKCoqIgtW7bQuXNnsrOzycnJIT09nWnTph31GEOHDuX555/nvPPOw+12s3btWlq1asWgQYOYMmUKN9xwA3l5ecydO5drrrnmqMeJVks27eGxT1bz7YbdpMs2nnG/yc+d3/xku63ahLd85/CW7xxytRmNE9xc27s1o/ql0bm5zXhmok/dC/1aIDU1lVdeeYWrr76akpLAPd1/+tOf6NSpE8899xzDhg0jMTGRfv36HfUYv/zlL8nJyaFPnz6oKqmpqbz//vtcfvnlfPbZZ3Tt2pW0tDT69OlDw4bVOI5Rx6zbWciTM1czc8UOUingUde7XO38HJf4Q9t41cGn/r5M8w3mS39PVByc2SGF+/ulcWHXU2zSDxPV7NXKVaywsJCkpCRUldtvv52OHTty1113ndQx8vPz6d+/P19//TXNmzf/yXZ16felsnbsO8g/Zq/lzaxc4v1FjHF9xC+dn5BQ7kGqj3yn85T3KjZoC1o0jOPKzDSu7NuatCZHnjHKmPrCXq0cIS+++CJTpkyhtLSU3r17M3bs2BM+xs9//nMKCgooLS3loYceOmLgR4u9BzxM/GI9k7/egN9TwvXO2YyLff8nUwku8HXlMe9olmkHTm2ezKQLO3F+l1PsVktjyrHQr2J33XXXCV/Zlzd37tyqKaYOO+jx8drCjTwzZx17i0sY7pjPPTFv/WTWqVX+NjzmvZov/D1p1SiBpy7sxKW9W1nYG3MUdSb0VdVupQtT27rlqorPr7y/dAsTPl3LloJiznYs5/6YqXR1bCyzXa6m8DfPVXzgP4OGCbE8OLgD1w5oa3O9GnMcdSL04+LiyM/Pp2nTphb8BAI/Pz+fuLi4SJdSZVSVuWvyePx/q1m9/fDtl4PK3X65W5N4xnsZr/kuwOGO5bazMxh7Tnsaxtsri42piDoR+q1btyY3N5e8vIpNKB0N4uLiyky4XpctDd5++c0xbr8s1lhe8l3Ei96fUyQJjOqfxp3nd6J5w/rzF58xNaFCoS8iw4CnASfwkqo+dpTtLgfeBvqpapaIpAOrgDXBTRaq6q0nWqTb7Q495WrqjwOlPv4yYxX/WbgxdPvlaOecMrNSedXBVN9gnvaOJI/GDOvWnHuGdqZDs6QIVm5M3XXc0BcRJ/AscCGQCywSkemqurLcdsnAnUD5J2TWq2qvKqrX1BM/bNnLnVOXkp23n186P+Fu19s/uf3yY19//uYdxQZtQf/0Jkz82an0adM4QhUbUz9U5Eq/P7BOVbMBRGQqMAJYWW67R4HHgXurtEJTr/j8ysQv1zNh1lpS/bv4r/t5znCW/U8p/PbLzqckM/mizgzu3MzGc4ypAhUJ/VbA5rDlXOD08A1EpA+Qpqofi0j50M8QkaXAPuBBVf2qMgWbuit3TzF3T1vGtzm7+bljAX+O/TcNpTi0fpU/jce81/CFvyctG8bztyGducxuvzSmSlV6IFdEHMAE4MYjrN4GtFHVfBHpC7wvIt1UdV+5Y4wBxgC0adOmsiWZWkZV+eC7rTz0/g9Qso+n3K9wuXNeaL1PhWd9I/indyRJCfH84dwOXDfQbr80pjpUJPS3AGlhy62DbYckA92BucF/fjcHpovIcFXNAkoAVHWxiKwHOgFl3rOgqpOASRB4DcPJnYqpjfYWe3jwgx/4cNlW+slq/h77XJlJTDb7U/k/z69ZrJ0Z0asl44d3p2GC3X5pTHWpSOgvAjqKSAaBsB8NhF77qKp7gZRDyyIyF7gnePdOKrBbVX0i0g7oCGRXYf2mFpu/fhf3vLmMvL2F3ON6h9uc03GGTWTytu9sHvFcj8Q14OlLuzOiV6sIVmtMdDhu6KuqV0TGATMJ3LI5WVVXiMh4IEtVpx9j97OB8SLiAfzAraq6uyoKN7VXidfHhFlrmfRVNu3Ywjsxz9HTsSG0vkATecBzCzP8Azg9owkTRvWilU0wbkyNqBNv2TR1x9od+7lz6nes2raXa52z+YPrv8RLaWj9PF837vHcSr4zhd8O6cyvzmpnA7XGVAF7y6apUarKlPk5/PWT1SR79/Bv9yTOdy4NrS9RF094RzPZN4z2zRrw0qhedG9l8wQYU9Ms9E2l7dx3kHveXs6Xa/M437GYx2NfLDM37Wp/Gv/nuZ3V2oYbBrbl9z/rYnfmGBMhFvqmUmau2M797yznYPF+/uz6L79wfVZm/Uvei3jSO4oGycm8fEVPBnduFqFKjTFgoW9OkqryzOfreOrTtfSQbP4R8yztHdtC67drY+7x3Mo8fw8u7HoKj43sQdOk2AhWbIwBC31zEjw+Pw+9/wNTF21mpONLHne/WOYlaTN8/XnAcwulMY14/JKuXJWZZq9QMKaWsNA3J6SoxMuv/7uEL9bu5A7ne/zW/XZoXaHG8Yj3Bt72nc1paY35x6heZKQkRrBaY0x5FvqmwnbuO8jNUxaxestu/up6matdc0LrVvnTGOO5my2cwp3nd2TceR1wOx0RrNYYcyQW+qZC1u3czw2TF1FQsJuX3P/kXOey0Lp5vm7c5rmLxAZNeOsXfejb1l5/bExtZaFvjuub7Hx+9WoWsQfzmBbzJN0dOaF17/jO4n7Pr2jfvDEv39SPFg3tyVpjajMLfXNM05dt5Z43l9HGv4lXYp8o87K0p72X8XfvFZzZIZXnr+1Dcpy9KM2Y2s5C3xyRqjLxy2we+2Q1AxwrmRQzgQbBd9971cEfvLcwzTeYy/u05q8jexDjsv57Y+oCC33zEz6/8sj0Ffxn4UaGO+bzpPsFYsULQJHGcrvnTub6e/Gb8zty1wUd7XZMY+oQC31TxoFSH3e8sZTZq7Zzq/ND7ndPDa3bqY24qfReVks7Hr+8O6P62YQ3xtQ1FvomZFdhCbdMyeKHzfk86prCda7ZoXU/+ltxY+l9FMQ0Z/K1fTmnU2oEKzXGnCwLfQNAdl4hN768iLzdu3nB/QwXOpeE1i30d2FM6V3EJTdl2o397O2YxtRhFvqGxRv38Mspi3AW7+KNmCfp5Tg8udl030Du8dxK22aBWzJbN06IYKXGmMqy0I9ySzbt4bp/f8MpnlymxDxGG0deaN3z3kt4wjuK/hkpTLou0+auNaYesNCPYiu37uPGyd/SzJPLmzHjSZW9APhU+KP3Rl7zXcjw01ry5JU9iXXZ+++NqQ8s9KNUdl4h10/+hqSD23kt9q+hwC/WWO7wjOMzf19uO7c99w7pjMOmMzSm3rDQj0JbCg5w7UvfQGEer8X8JfSUbbHGcn3p71jCqfzp0u5cO6BthCs1xlQ1C/0ok7e/hGtf+obCvbuYFvNX2jm2A4E5bMd67mIJp/LMNX34WY8WEa7UGFMdLPSjyN5iD9f9+xt27MrntZgn6OLYBAT68H/jGcdX/p787crTLPCNqccs9KNEUYmXG1/5lg3b8/m3+yn6ONaF1t3rGctMf38euaQrV/RtHcEqjTHVzUI/Chz0+PjVq1l8v2kXz7n/xZnOFaF1D3lu5F3/2dwzpBM3DsqIYJXGmJpgoV/PeXx+xr2+lAXr85jgnsgQ5+LQuic8V/Ef3xDGnt2O2wd3iGCVxpiaUqH34YrIMBFZIyLrROT+Y2x3uYioiGSGtf0+uN8aERlaFUWbivH7lXveWsbsVdt51PUylzm/Dq17wXsJz/lGcM3pbbj/olPtTZnGRInjXumLiBN4FrgQyAUWich0VV1Zbrtk4E7gm7C2rsBooBvQEpgtIp1U1Vd1p2CORFV56IMf+OC7rfzONZVrXZ+F1r3mPZ/HvKMZ0asVj47oboFvTBSpyJV+f2CdqmaraikwFRhxhO0eBR4HDoa1jQCmqmqJqm4A1gWPZ6qRqvLY/1bz32828WvnB9zm+jC07j3fIB7y3sQFXU7hb1eehtMevDImqlQk9FsBm8OWc4NtISLSB0hT1Y9PdN/g/mNEJEtEsvLy8sqvNifoubnrmfhFNtc7Z3Kfe1qofZavL/d6xjKwfSrPXNMHt9NmuzIm2lT6/3oRcQATgN+e7DFUdZKqZqpqZmqqvae9Ml5dkMOTM9cw0vEl491TQu3zfN24w3MH3dNSePH6TOLc9i4dY6JRRe7e2QKkhS23DrYdkgx0B+YG+4abA9NFZHgF9jVVaO6anfxx+gqGOhbxpHtiqH2JvwNjPL8lo3lTXrmpH4mxdtOWMdGqIlf6i4COIpIhIjEEBmanH1qpqntVNUVV01U1HVgIDFfVrOB2o0UkVkQygI7At1V+FoacXUX85o2ldCWHp93P4BQFYJW/DTeW3kezpk149Zb+NEqIiXClxphIOu4ln6p6RWQcMBNwApNVdYWIjAeyVHX6MfZdISJvAisBL3C73blT9YpKvIz5TxbOg7uZGDuBOPEAkO1vznWlvyepYQqv/fJ0miXHRbhSY0ykiapGuoYyMjMzNSsrK9Jl1Bmqyu2vL2Hm91t41f0Yg4JP2+7TeC4tfZS9Cem8eetA2qcmRbhSY0x1EpHFqpp5vO2sc7eOe/6L9cz4fjt/cL0RCnyAuzy/ZpO04vVr+1rgG2NC7J69Omzump08OXMNwx1f8yvXjFD7BM8VfObvy8OXdKV/RpMIVmiMqW0s9OuoQwO3XcjhcfeLofZZvr78y3cpV/ZtzXU2CYoxphwL/TqoqMTL2P8sxnlwN5NiJhAvpQCs87fkbs9t9ExrwqOX2usVjDE/ZX36dYyqcu/by1i3o4Ap7n+Fpjrcr/GM8dxNXFIjXri2jz18ZYw5Igv9OubQwO0Drqll3osfGrj9RV9aNIyPYIXGmNrMunfqkPCB2zGuw685+rvncmbbwK0xpgIs9OuIjflHHrj91NeXf/ous4FbY0yFWOjXAUUlXsa8uhjHwT1MdP89NHC73t+Cuzy30bN1Yxu4NcZUiPXp13JlB27/SZoj8OrpMgO31/W1gVtjTIVY6NdyL321gRnfb+f3Rxi43SitbeDWGHNCrHunFlu9fR9PzlzDUMe3jA0buP2Hd6QN3BpjToqFfi1V6vVz17RlNPTt5q/ul0Ltn/r68rR3pA3cGmNOinXv1FJPf7aWVdv28pL7RZpIIQBbtQm/9dxKDxu4NcacJLvSr4UWb9zN83PXc5VzLhc4l4ba7/Hcii+mAc9eY0/cGmNOjl3p1zJFJV7ufnMZLdnJw67/hNpf9g5lvr87j1/SlbQmCRGs0BhTl1no1zJ//WQVm/MLeSPmBZLkIBC4H/9x72gu6NKMqzLTjnMEY4w5OuveqUW+WJvHaws3cYtzBqc7VgPgVQd3eX5NQmIyfx3Z0/rxjTGVYlf6tURBcSn3vrWMTrKZe1xvhtqf8V3Kcm3P85d2JzU5NoIVGmPqAwv9WuKhD1awZ38Rr8Q8R6x4AVjuz+AZ76WM7N2Ki3q0iHCFxpj6wLp3aoHpy7by4bKt3Ol6h66OjQCUqJu7PL8mtWESfxzeLcIVGmPqC7vSj7Ad+w7y0Ps/0EfWcptzeqj9ce9o1msr/nvlaTSMd0ewQmNMfWJX+hEUeJnackoP7Ocp9/M4RQGY7+vKy76h3HhGOoM6pES4SmNMfWKhH0H//WYTX67N4wHX62Q4dgCwT+O5x3MrGanJ/G7YqRGu0BhT31Qo9EVkmIisEZF1InL/EdbfKiLfi8h3IjJPRLoG29NF5ECw/TsReaGqT6CuytlVxJ8/XsXZjmVc55odav9/nhvY4UhlwlW9iI+xp26NMVXruH36IuIEngUuBHKBRSIyXVVXhm32uqq+ENx+ODABGBZct15Ve1Vt2XWb3x94R36MZy9PxE4Ktf/P1493/Gfxm/M70CutUQQrNMbUVxW50u8PrFPVbFUtBaYCI8I3UNV9YYuJgFZdifXPu0u3sChnDw+6XqO57AEgTxvwgOcWerRqxB3ndYhwhcaY+qoiod8K2By2nBtsK0NEbheR9cATwG/CVmWIyFIR+UJEzqpUtfXA3gMe/jpjFX1lDVe6vgy1P+D5JYWuRvx91Gm4nTbUYoypHlWWLqr6rKq2B34HPBhs3ga0UdXewN3A6yLSoPy+IjJGRLJEJCsvL6+qSqqVJsxaQ0HRAR51vxJq+8TXj0/9mdw3tDMdmiVHrDZjTP1XkdDfAoS/5at1sO1opgKXAqhqiarmB78vBtYDncrvoKqTVDVTVTNTU1MrWnud88OWvfxn4Uaudc4OPYR1QGN41HMdXVo04MYun54FAAAUe0lEQVQz0iNboDGm3qtI6C8COopIhojEAKOB6eEbiEjHsMWLgR+D7anBgWBEpB3QEciuisLrGr9fefiDH2iie/lt2Lt1/uW9jK2k8OiIbrisW8cYU82Oe/eOqnpFZBwwE3ACk1V1hYiMB7JUdTowTkQuADzAHuCG4O5nA+NFxAP4gVtVdXd1nEht9/aSXJZsKuAp9+s0kANA4JXJL/l+xhV9W5OZbnPdGmOqX4Vew6CqM4AZ5doeDvt+51H2ewd4pzIF1gd7iz089slq+slqLnd+FWp/xHsDsXHx3H+RPYRljKkZ1p9QA/42aw17iw4w3v1yqO1jX3++8vfk3qGdSUmyVyYbY2qGhX41+z53L699s5HrnbPo4gjc+VqssfzJcx3dWjbgF6e3jXCFxphoYm/ZrEZ+v/LQBz+Qonu4y/V2qP2f3svYRlOeGdEdp8NmwjLG1BwL/Wr01uLNfLe5gAnuN8oM3v7b9zOuymxN37aNI1yhMSbaWPdONSkoLuWxT1bTX1Yx0jkv1P6w90bi4+LsDZrGmIiw0K8mT85cw/7iA4x3vxJq+8g3gK/9Pbh32Kk0tcFbY0wEWOhXg+W5Bbz+7SZucM7i1ODgbZHG8ifPL+jeqgHX9G8T4QqNMdHK+vSrmKry8AcrSNU9/J/r8CMKT3tHsp2mPG+Dt8aYCLLQr2L/+2E7320u4G/uaSQHB29/9LfiZd9FjO6XRu82NnhrjIkc696pQl6fnydnraGTbGak4/CTt3/03kBCfDz32eCtMSbC7Eq/Cr21OJfsvCJedE/DEZzk/HNfL+b7uzN+SCeaJMZEuEJjTLSzK/0qctDj4x+z19JX1nChcwkAfhWe8I4mvWkCV9vgrTGmFrAr/Sryyvwcduw7yL9ipobaPvCfwWptw7+GdLbZsIwxtYIlURXYW+zhuTnrGOz4jv6ONQCUqpMJ3ivo1rIBF/doEeEKjTEmwK70q8ALX66n8GAp94Vd5b/uO5/NegqvDjsVh92iaYypJexKv5J27DvIy19vYLhjfugtmkUayzPeyxjYrilndUyJcIXGGHOYXelX0tOf/YjfU8JvY94Ktb3k+xm7aMhLF52KiF3lG2NqD7vSr4TsvEKmLdrM1c7PSXPkAbBbk3jRezEXdW9Or7RGEa7QGGPKstCvhKdmrSXOX8wdrvdCbc96L6VYEvjtkM4RrMwYY47MQv8kLc8t4OPvt/FL5wxSZB8AuZrCa74LuCozjQ7NkiJcoTHG/JSF/kl64n9raMI+fuX6ONT2D+/l4Irjzgs6RrAyY4w5Ogv9kzDvx13MW7eLca73SZKDAKzxt+Zd31nceEY6LRrGR7hCY4w5Mgv9E6SqPDlzNa0lj184Z4fan/SOIjEuhtvObR/B6owx5tgs9E/QF2vzWJa7lzuc7xErXgCy/J2Y7e/Dree0p1GCvVTNGFN7WeifoGfnrKMVeYx0Hn518pOeUaQmx3HToPTIFWaMMRVQodAXkWEiskZE1onI/UdYf6uIfC8i34nIPBHpGrbu98H91ojI0KosvqZ9k53Popw9jHF9hFt8gTb/qXyjXRg3uAMJMfasmzGmdjtu6IuIE3gWuAjoClwdHupBr6tqD1XtBTwBTAju2xUYDXQDhgHPBY9XJz0zZx2pFDDaOfdwm/dSUpJiGdUvLXKFGWNMBVXkSr8/sE5Vs1W1FJgKjAjfQFX3hS0mAhr8PgKYqqolqroBWBc8Xp2zbHMBX/24i1tcnxArnkCbvx1f+Xvwq7MyiHPX2b/LjDFRpCL9Ea2AzWHLucDp5TcSkduBu4EY4LywfReW27fVEfYdA4wBaNOmdk428tzcdTSkkGudn4banvWOoGF8DL8Y0DaClRljTMVV2UCuqj6rqu2B3wEPnuC+k1Q1U1UzU1NTq6qkKrN2x35mrtjBjc6ZZe7L/9Tfl5sGpZMUa335xpi6oSKhvwUI77BuHWw7mqnApSe5b6303Jx1JHKAm1z/O9zmHU5CjJsbz0iPXGHGGHOCKhL6i4COIpIhIjEEBmanh28gIuHvHbgY+DH4fTowWkRiRSQD6Ah8W/mya87G/CKmL9vKL5yzaSRFgTZ/Mz7yD+TaAW3tvnxjTJ1y3H4JVfWKyDhgJuAEJqvqChEZD2Sp6nRgnIhcAHiAPcANwX1XiMibwErAC9yuqr5qOpdq8cIX63FrKb9yzQi1Pe8bjtPl5pazMiJYmTHGnLgKdUar6gxgRrm2h8O+33mMff8M/PlkC4ykbXsP8PbiXEY755IqewNt2oR3fWcxemAazZLjIlyhMcacGHsi9xgmfZmN+jyMdX10uM17MX5HDGPPsXfsGGPqHgv9o8gvLOGNbzdxmXMerWVXoE2TecN3Hpf1bkWrRvYmTWNM3WOhfxSTv95AqcfLrc4PQ23/9l5EqcTamzSNMXWWhf4R7D/o4dX5G7nI8S3tHdsA2KcJ/Mc3hJ/1aEG7VJsVyxhTN1noH8FbWbnsL/Ew1nX4Kn+Kbwj7SeDX53aIYGXGGFM5Fvrl+PzKK/NzyJQ19HRsAOCgunnZO4zzT21G15YNIlyhMcacPAv9cj5btYNNu4u5Oezp2/d8Z7KbBow5u10EKzPGmMqz0C/n5a9zaEUeQx2LDrf5htGtZQP6ZzSJYGXGGFN5FvphVm7dx4LsfK53zcIpgbdDf+XrzlpN4+ZBGYhIhCs0xpjKsdAP8/LXG0jgIFc754TaJvsuIiUplp+f1iKClRljTNWw0A/aVVjCB99t5XLnlzSQYgCy/c2Z6z+Nawe0IdZlk6QYY+o+C/2g/y7chMfn5UbnzFDbK76huJ0ufnG6TZJijKkfLPSBEq+P177ZyDmOZWUexnrbdw7De7UkNTk2whUaY0zVsNAHPl6+jbz9JdzsPHyb5lTfYIqJ46ZB6ZErzBhjqljUh76q8u95G+gouZzt/B4Anwqv+oZwekYTurVsGOEKjTGm6kR96C/K2cOKrfu4yflJqG2WP5NcTeXmM22SFGNM/RL1oT953gYasZ+RznmH27wXkdYkngu6nBLByowxpupFdehv3l3MrJXbucb5OXHiAeB7fzqLtDM3npGB02EPYxlj6peoDv3Xv92EqI/rXJ+G2iZ7LyIxxsWVma0jWJkxxlSPqA39Uq+ft7I2c55jKS1kNwB52pCP/QO4MjONBnHuCFdojDFVL2pD/7NVO9hVWMrVzs9DbW/6zqEUN9cPtIexjDH1U9SG/uvfbqIVeZzrWBZqm+obzMB2TW1mLGNMvRWVob8pv5ivftzFVa65OIJv0/zS14PNegrXnN4mwtUZY0z1icrQn7poE058jHLODbW97jufJokxDOlmt2kaY+qvCoW+iAwTkTUisk5E7j/C+rtFZKWILBeRz0Skbdg6n4h8F/xMr8riT4bH5+fNrFzOcyyluewBAgO4s/19uLJva3ubpjGmXnMdbwMRcQLPAhcCucAiEZmuqivDNlsKZKpqsYjcBjwBjAquO6Cqvaq47pM2e+UOdhWWcI37s1Dbm75z8OJidH/r2jHG1G8VudLvD6xT1WxVLQWmAiPCN1DVOapaHFxcCNTam9wPDeCe41geapvqG8wZ7ZuSkZIYwcqMMab6VST0WwGbw5Zzg21HcwvwSdhynIhkichCEbn0JGqsMscawL3arvKNMVHguN07J0JErgUygXPCmtuq6hYRaQd8LiLfq+r6cvuNAcYAtGlTfeF7tAHcpokxDO3WvNp+XWOMqS0qcqW/BUgLW24dbCtDRC4A/gAMV9WSQ+2quiX4MxuYC/Quv6+qTlLVTFXNTE1NPaETqKhjDeBe0bc1Ma6ovJHJGBNlKpJ0i4COIpIhIjHAaKDMXTgi0huYSCDwd4a1NxaR2OD3FGAQED4AXGNCA7hOG8A1xkSv43bvqKpXRMYBMwEnMFlVV4jIeCBLVacDTwJJwFsiArBJVYcDXYCJIuIn8BfMY+Xu+qkxNoBrjDEV7NNX1RnAjHJtD4d9v+Ao+80HelSmwKqwteAA89bt4v9cX/xkAPc+u8o3xkSRqOjIfm/pFlA/Ix1fhdqm+gbbE7jGmKhT70NfVXl3SS79ZA1pjjwACjSR2f6+XNa7lT2Ba4yJKvU+9Jfl7mV9XhGXOw9f5X/oG0gpbi7vU2ufITPGmGpR70P/3SW5xFHCz5zfHG7zncWpzZPp2rJBBCszxpiaV69Dv8TrY/qyrQxxZJEsBwBY72/BUu3AFX3tKt8YE33qdejPWZ1HQbGnTNfOu76zcDocDO/VMoKVGWNMZNTr0H9nSS7N2MOZju9Dbe/5zuTsjik0S46LYGXGGBMZ9Tb08wtLmLN6J5c65+EM3ps/39eVraQw0gZwjTFRqt6G/ofLtuL1+8t07bzjO5vkOBcXdrV7840x0anehv47S7bQTXLo7MgFoFhj+cTfn5/3bEGc2+7NN8ZEp3oZ+mt37Of7LXu5wvllqO0Tfz+KibN7840xUa1ehv47S3Jx4WW4c/7hNt/ZtG2aQN+2jSNYmTHGRFa9C32/X/lg6VbOcSyjqewHYKs2YaG/KyN7tyb4FlBjjIlK9S70szbuYfu+g4wIu8p/33cmfhyM7HOsWR6NMab+q3eh/+GyrcRzkAscS0JtH/jOILNtY9KaJESwMmOMibx6Ffpen58Z32/jfMdSEiQwY+NafyvWaJo9gWuMMdSz0F+QnU9+USmXOBeE2j70DcQhwkXdW0SwMmOMqR3qVeh/uGwryRRzruO7UNtH/oGc0T6F1OTYCFZmjDG1Q70J/RKvj//9sJ0hjixixQvA9/50NmgLLjnNrvKNMQbqUeh/tXYX+w56f9K143YKQ7s1j2BlxhhTe9Sb0P9w+VYas6/MGzU/9g3g7I6pNEqIiWBlxhhTe9SL0D9Q6uPTlTu4yLkIl/gByPJ3YgupXHKa3bVjjDGH1IvQ/3z1TopLfVziKNu1E+tycIG9UdMYY0LqReh/uGwrzdjD6Y5VAPhUmOE7nfO7NCMp1hXh6owxpvaoUOiLyDARWSMi60Tk/iOsv1tEVorIchH5TETahq27QUR+DH5uqMriAfYf9PD5mp1c7FyIIzhZykJ/V/JoxCU9rWvHGGPCHTf0RcQJPAtcBHQFrhaRruU2WwpkqmpP4G3gieC+TYA/AqcD/YE/ikiVvuby05U7KPX6y9614x9IYoyTwac2q8pfyhhj6ryKXOn3B9aparaqlgJTgRHhG6jqHFUtDi4uBA69tH4o8Kmq7lbVPcCnwLCqKT3gw2VbSaWAXrIeAI86+Z+vH0O6NbfJUowxppyKdHi3AjaHLecSuHI/mluAT46xb5W96vKgx0fWxj3spxFnlPyTi50LOUUKKCDZHsgyxpgjqNJRThG5FsgEzjnB/cYAYwDatGlT4f3i3E6+eeB8Plu1kw+XbeU/a1Ip9fppGO/mzA6pJ1KCMcZEhYqE/hYgLWy5dbCtDBG5APgDcI6qloTte265feeW31dVJwGTADIzM7UCNYUkxLi45LSWXHJaS/Yd9DBrxQ72H/QQ46oXNyYZY0yVqkjoLwI6ikgGgRAfDVwTvoGI9AYmAsNUdWfYqpnAX8IGb4cAv6901UfRIM7NFX1tDlxjjDma44a+qnpFZByBAHcCk1V1hYiMB7JUdTrwJJAEvBWcjnCTqg5X1d0i8iiBvzgAxqvq7mo5E2OMMcclqifUm1LtMjMzNSsrK9JlGGNMnSIii1U183jbWce3McZEEQt9Y4yJIhb6xhgTRWpdn76I5AEbI11HBaUAuyJdRA2LxnOG6DxvO+e6pa2qHvcBpVoX+nWJiGRVZOCkPonGc4boPG875/rJuneMMSaKWOgbY0wUsdCvnEmRLiACovGcITrP2865HrI+fWOMiSJ2pW+MMVHEQv8oKjBFZKyITAuu/0ZE0oPtF4rIYhH5PvjzvJqu/WSd7DmHrW8jIoUick9N1VxZlTlnEekpIgtEZEXwzzuuJms/WZX4b9stIlOC57pKRKrt5YnVoQLnfbaILBERr4hcUW5dtU77WqNU1T7lPgReLLceaAfEAMuAruW2+TXwQvD7aGBa8HtvoGXwe3dgS6TPp7rPOWz928BbwD2RPp8a+HN2AcuB04LLTQFnpM+pms/5GmBq8HsCkAOkR/qcqvC804GewKvAFWHtTYDs4M/Gwe+NI31OJ/uxK/0jO+4UkcHlKcHvbwPni4io6lJV3RpsXwHEi0hsjVRdOSd9zgAicimwgcA51xWVOechwHJVXQagqvmq6quhuiujMuesQKKIuIB4oBTYVzNlV1pFpn3NUdXlgL/cvtU+7WtNstA/sopM8xjaRlW9wF4CV3vhLgeW6OFJZWqzkz5nEUkCfgf8vxqosypV5s+5E6AiMjPYJXBfDdRbFSpzzm8DRcA2YBPwN607r0qvzNSt1Trta02r0ukSzWEi0g14nMAVYX33CPB3VS0MXvhHAxdwJtAPKAY+C77a9rPIllWt+gM+oCWBbo6vRGS2qmZHtixzIuxK/8gqMkVkaJvgP3cbAvnB5dbAe8D1qrq+2qutGpU559OBJ0QkB/g/4IHgxDu1XWXOORf4UlV3qWoxMAPoU+0VV15lzvka4H+q6tHADHlfE5gTuy6o0LSv1bBvrWOhf2ShKSJFJIbAYNb0cttMBw6N4l8BfK6qKiKNgI+B+1X16xqruPJO+pxV9SxVTVfVdOAfwF9U9ZmaKrwSTvqcCcwk10NEEoLBeA6wsobqrozKnPMm4DwAEUkEBgCra6TqyqvIeR/NTGCIiDQOTv06JNhWN0V6JLm2foCfAWsJjPj/Idg2Hhge/B5H4E6VdcC3QLtg+4ME+j2/C/s0i/T5VOc5lzvGI9SRu3cqe87AtQQGrn8Anoj0uVT3OROcEjV4ziuBeyN9LlV83v0I/AuuiMC/bFaE7Xtz8PdjHXBTpM+lMh97ItcYY6KIde8YY0wUsdA3xpgoYqFvjDFRxELfGGOiiIW+McZEEQt9Y4yJIhb6xhgTRSz0jTEmivx/we6PbgdKCd0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "SGL_vals = []\n",
    "regreg_vals = []\n",
    "rel_error = []\n",
    "for i, lagrange in enumerate(lagrange_val):\n",
    "    penalty.lagrange = lagrange\n",
    "    SGL_soln = B[:,i] #np.hstack([I, B[:,i]])\n",
    "    SGL_vals.append(problem.objective(SGL_soln))\n",
    "    regreg_vals.append(problem.objective(solns[i]))\n",
    "    rel_error.append(np.linalg.norm(solns[i] - SGL_soln) / max(np.linalg.norm(solns[i]), 1))\n",
    "SGL_vals = np.array(SGL_vals)\n",
    "regreg_vals = np.array(regreg_vals)\n",
    "plt.plot(lagrange_val, SGL_vals, label='SGL', linewidth=6)\n",
    "plt.plot(lagrange_val, regreg_vals, label='regreg', linewidth=3)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1367cb358>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8XXWd//HX596bpVnaLE1K96RtWlr2koZFEFmUIkodgR8FF0aL/BxlHMffqGV0GMVhBHV0HAEdFH6DihZkhrGjSJGiP0CgbRrWtpSGLmkLdL1Nt+z5/P64JzWkN8ltk9wleT8fjz567jnf872fr8V88j2fc77H3B0REZHehFIdgIiIpDclChER6ZMShYiI9EmJQkRE+qREISIifVKiEBGRPilRiIhIn5QoRESkT0oUIiLSp0gijcxsPvB9IAz8xN1v73E8B/gpcCawB7jG3TcHx24GFgEdwOfcfVmw/z7gA8BOdz+5W1/fBj4ItAJvAJ9w9319xTd27FivqKhIZCgiIhJYvXr1bncv66+d9beEh5mFgdeB9wLbgFXAte6+tlubzwCnuvunzWwh8Bfufo2ZzQF+CdQAE4AngJnu3mFm7wYOAj/tkSjeBzzp7u1mdgeAu3+5rxirq6u9tra2v7GKiEg3Zrba3av7a5fIpacaoN7dN7p7K7AEWNCjzQLg/mD7YeBiM7Ng/xJ3b3H3TUB90B/u/hSwt+eXufvj7t4efHwemJRAjCIiMkQSSRQTga3dPm8L9sVtE/yQbwRKEzy3L58EfhfvgJndaGa1Zla7a9euY+hSRESORdoWs83sK0A78EC84+5+j7tXu3t1WVm/l9hEROQ4JZIotgOTu32eFOyL28bMIsAYYkXtRM49ipn9JbFC90dc66CLiKRUIoliFVBlZpVmlg0sBJb2aLMUuD7YvopYMdqD/QvNLMfMKoEqYGVfXxbcYfUl4Ap3P5z4UEREZCj0myiCmsNNwDJgHfCQu68xs1vN7Iqg2b1AqZnVA18AFgfnrgEeAtYCjwGfdfcOADP7JfAcMMvMtpnZoqCvO4FC4Pdm9qKZ/WiQxioiIseh39tjM4FujxUROXaDeXvssLV83Q5++Mc3Uh2GiEhaG9GJ4ukNu7n7j/WpDkNEJK2N6ERRlJfFgeZ22js6Ux2KiEjaGtGJojgvG4B9TW0pjkREJH2N6ERRlJcFwL7DrSmOREQkfY3oRNE1o4ge1oxCRKQ3ShRA9JBmFCIivRnRieLIpSfVKEREeqVEgWoUIiJ9GdGJoiAnQiRkqlGIiPRhRCcKM6MoL1szChGRPozoRAFQnJdF9JBmFCIivVGiyMsmqhmFiEivRnyiKMrLYp9qFCIivRrxiUIzChGRvo34RFGUH5tRDIf3coiIDIURnyiK87Jp7ejkcGtHqkMREUlLIz5RFI3S09kiIn1RotB6TyIifRrxiaL4yDIemlGIiMSjRJHftdS4ZhQiIvGM+EShhQFFRPqmRDFKLy8SEenLiE8U2ZEQBTkRXXoSEenFiE8UoGU8RET6klCiMLP5ZrbezOrNbHGc4zlm9mBwfIWZVXQ7dnOwf72ZXdpt/31mttPMXu3RV4mZ/d7MNgR/Fx//8BKjZTxERHrXb6IwszBwF3AZMAe41szm9Gi2CIi6+wzge8AdwblzgIXAScB84O6gP4D/CPb1tBhY7u5VwPLg85AqystSjUJEpBeJzChqgHp33+jurcASYEGPNguA+4Pth4GLzcyC/UvcvcXdNwH1QX+4+1PA3jjf172v+4EPHcN4joteXiQi0rtEEsVEYGu3z9uCfXHbuHs70AiUJnhuT+Pc/a1g+21gXAIxDkixahQiIr1K62K2x5Z0jbusq5ndaGa1Zla7a9euAX1PUV42+5vb6OjUCrIiIj0lkii2A5O7fZ4U7IvbxswiwBhgT4Ln9rTDzMYHfY0HdsZr5O73uHu1u1eXlZUlMIzeFedl4Q6NWhhQROQoiSSKVUCVmVWaWTax4vTSHm2WAtcH21cBTwazgaXAwuCuqEqgCljZz/d17+t64NcJxDggxXlaxkNEpDf9Joqg5nATsAxYBzzk7mvM7FYzuyJodi9Qamb1wBcI7lRy9zXAQ8Ba4DHgs+7eAWBmvwSeA2aZ2TYzWxT0dTvwXjPbAFwSfB5SWsZDRKR3kUQaufujwKM99t3SbbsZuLqXc28Dbouz/9pe2u8BLk4krsFyZEZxSJeeRER6SutidrLo0pOISO+UKIi9Nxv0TgoRkXiUKIDCnAjhkGlGISIShxIFYGYUjdIyHiIi8ShRBGIryGpGISLSkxJFoDgvWzUKEZE4lCgCRVpqXEQkLiWKgBYGFBGJT4kiUJyvGYWISDxKFIGivCxa2jtpau1IdSgiImlFiSKgp7NFROJToggUBwsDKlGIiLyTEkWgKJhRqKAtIvJOShSBIs0oRETiUqII/LlGoRmFiEh3ShSBrhlFo2YUIiLvoEQRyImEycsOa0YhItKDEkU3xVrGQ0TkKEoU3RRpGQ8RkaMoUXSjGYWIyNGUKLrRjEJE5GhKFN1oRiEicjQlim6K87JobGqjo9NTHYqISNpQouimKC8bd9jfpMtPIiJdlCi60TIeIiJHSyhRmNl8M1tvZvVmtjjO8RwzezA4vsLMKroduznYv97MLu2vTzO72MzqzOxFM3vGzGYMbIiJ0zIeIiJH6zdRmFkYuAu4DJgDXGtmc3o0WwRE3X0G8D3gjuDcOcBC4CRgPnC3mYX76fOHwEfc/XTgF8BXBzbExJ0wJheAhr2HkvWVIiJpL5EZRQ1Q7+4b3b0VWAIs6NFmAXB/sP0wcLGZWbB/ibu3uPsmoD7or68+HRgdbI8B3jy+oR27meMKGZ0b4bk39iTrK0VE0l4kgTYTga3dPm8Dzuqtjbu3m1kjUBrsf77HuROD7d76vAF41MyagP3A2QnEOCjCIePsaaU8q0QhInJEOhaz/xZ4v7tPAv4v8N14jczsRjOrNbPaXbt2DdqXnzu9lG3RJrbuPTxofYqIZLJEEsV2YHK3z5OCfXHbmFmE2CWjPX2cG3e/mZUBp7n7imD/g8C58YJy93vcvdrdq8vKyhIYRmLOnTEWQJefREQCiSSKVUCVmVWaWTax4vTSHm2WAtcH21cBT7q7B/sXBndFVQJVwMo++owCY8xsZtDXe4F1xz+8Y1dVXsDYgmyefWN3Mr9WRCRt9VujCGoONwHLgDBwn7uvMbNbgVp3XwrcC/zMzOqBvcR+8BO0ewhYC7QDn3X3DoB4fQb7PwX8p5l1EkscnxzUEffDLFaneG7jHtydWE1eRGTkstgv/pmturraa2trB62/X6xo4O8feYXl/+cCppcVDFq/IiLpxMxWu3t1f+3SsZidcudMLwXQ3U8iIihRxFVRmsf4Mbk8r0QhIqJEEY+Zcc70WJ2iUyvJisgIp0TRi3Onj2XvoVbW7zgQ9/iKjXuoa4gmOSoRkeRTouhFV50i3vMUa9/cz8fvW8nXlq5JdlgiIkmnRNGLiUWjmFqad1RBu7Gpjb96YDUt7Z3U7zyoS1MiMuwpUfTh3OmlrNi058gb79ydL/7qJbZHm7hy7iQOt3bwZmNTiqMUERlaShR9OHtaKQea21nzZiMA9zy1kcfX7uDm989mYU1sBZINOw+mMkQRkSGXyOqxI1b35ykOt3Zwx2Ovcfkp4/nkuyrYF7zcqH7HQS6cVZ7KMEVEhpQSRR/KC3OpKi/gd6+8xU+e3kTF2Hxuv/IUzIzi/GzGFuSwYWf8u6JERIYLXXrqx7nTS3lpWyOHWtr50UfPpDA368ixqvIC6nXpSUSGOSWKfrwnuKz0zQ+fwsxxhe84VjWugA07DzIc1ssSEemNLj314z2zynj+5ouPvE+7u6ryAg40t7PzQAvjRh99XERkONCMoh9mFjdJAMwoj80wNuzQ5ScRGb6UKAagalxsCXIVtEVkOFOiGIDS/GyK8rL0LIWIDGtKFANgZrE7n3TpSUSGMSWKAZpRXsjrOw/ozicRGbaUKAaoqryAfYfb2HOoNdWhiIgMCSWKATpS0NblJxEZppQoBqgquEW2Xnc+icgwpUQxQONG51CYE9GdTyIybClRDJCZMWOc1nwSkeFLiWIQVJUXaEYhIsOWEsUgqCovZNeBFvYd1p1PIjL8JJQozGy+ma03s3ozWxzneI6ZPRgcX2FmFd2O3RzsX29ml/bXp8XcZmavm9k6M/vcwIY49GYEdz7p8pOIDEf9JgozCwN3AZcBc4BrzWxOj2aLgKi7zwC+B9wRnDsHWAicBMwH7jazcD99/iUwGTjR3WcDSwY0wiSoKu9a80mJQkSGn0RmFDVAvbtvdPdWYj+4F/RoswC4P9h+GLjYzCzYv8TdW9x9E1Af9NdXn38F3OrunQDuvvP4h5ccE8aMIi87rGcpRGRYSiRRTAS2dvu8LdgXt427twONQGkf5/bV53TgGjOrNbPfmVlVYkNJnVDImF5WoFVk5Yim1g7aOzpTHYbIoEjHYnYO0Ozu1cCPgfviNTKzG4NkUrtr166kBhiPXosqXdydy//tab61bH2qQxEZFIkkiu3EagZdJgX74rYxswgwBtjTx7l99bkN+K9g+xHg1HhBufs97l7t7tVlZWUJDGNozRhXwFuNzRxobkt1KJJi26JNbNx9iD+8lvZXTUUSkkiiWAVUmVmlmWUTK04v7dFmKXB9sH0V8KTHllNdCiwM7oqqBKqAlf30+d/AhcH2BcDrxze05PrzUh6aVYx0dQ1RIHZzw14tFinDQL+JIqg53AQsA9YBD7n7GjO71cyuCJrdC5SaWT3wBWBxcO4a4CFgLfAY8Fl37+itz6Cv24ErzewV4JvADYMz1KGlO5+kS92W6JHtVZv3pjASkcERSaSRuz8KPNpj3y3dtpuBq3s59zbgtkT6DPbvAy5PJK50Mrkkj+xISDMKoa5hH2dOLeaV7Y2s2rSXS086IdUhiQxIOhazM1K4686nHbrzaSQ73NrO2rf2c860Uk6fVKQZhQwLShSDqKq8gPpdmlGMZC9va6Sj05k7tYiayhJefXM/h1raUx2WyIAoUQyiqvICtkWbONyqHwwjVVch+4zJxcyrLKGj03mhYV+KoxIZGCWKQVQ1rgB32LjrUKpDkRSp2xJlWlk+xfnZzJ1SRMhg5aY9qQ5LZECUKAbRjOAWWT2hPTK5O3UN+5g7pRiAwtws5kwYzUrVKSTDKVEMoqmleWSFTWs+jVBb9hxm76FWzpxafGTfvIoSXmjYR2u7lvOQzKVEMYiywiEqx+brWYoRanXw/ETXjALgrMoSWto7eWV7Y6rCEhkwJYpBNkNrPo1YdQ1RCnMiRx6+BKiuKAH04J1kNiWKQTZr3Gi27DnE243NqQ5Fkmz1liinTykiFLIj+8YW5DCtLJ+Vm5QoJHMpUQyyD8+diJnxk6c3pjoUSaKDLe28vuPAOy47dampKKF28146Oz0FkYkMnBLFIJtckscVp03gFysbiGpBuBHjpa376HSYO/XoRDGvooT9ze2s11P7kqGUKIbAX71nOodbO/iPZzenOhRJktVbopjB6ZOLjjpWU6k6hWQ2JYohMHNcIZfMHsd/PLtZyzeMEHUNUarKCxgzKuuoY5OKRzF+TK7qFJKxlCiGyGcunE5jUxu/XNmQ6lBkiHUGy3ScGeeyE4CZMa+ihFWb9xJ7TYtIZlGiGCJzpxRzzrRSfvz0RlraO1IdjgyhjbsP0tjUxhlxCtld5lWWsGN/Cw17DycxMpHBoUQxhD5z4XR27G/hkbqeb46V4aRuS2zRv3h3PHWpCZ6n0OUnyURKFEPovBljOWXiGP79qY106NbIYWv1lihFeVlMG5vfa5uu+oUK2pKJlCiGkJnxmfdMZ9PuQ/zu1bdSHY4MkbqGKGdMfueDdj2FQl11imivbUTSlRLFELv0pBOYVpbPXX94Q4XMYaixqY0NOw/2WsjurqaymE27D7HzgJ7al8yiRDHEQiHj0xdMZ91b+/nj67tSHY4Mshcajl4IsDfzgjpFrWYVkmGUKJLgQ6dPZPyYXH74hzdSHYoMsrqGfYQMTovzoF1PJ08cw6issAraknGUKJIgOxLiU+dPY+XmvdSqmDms1G2JcuIJo8nPifTbNisc4owpRUoUknGUKJJkYc1kSvKzufuPmlUMFx2dzotb9zF3av+ziS41lSWse3s/+5vbhjAykcGlRJEkedkRPnFuBU++tlOzimFiw84DHGxpT6iQ3aWmogT3P7/kSCQTKFEk0aLzKxk/Jpd/+PUa2jv0asxMF++Ndv05Y0oxkZCxSpefJIMklCjMbL6ZrTezejNbHOd4jpk9GBxfYWYV3Y7dHOxfb2aXHkOf/2Zmw+pVcXnZEb56+RzWvbWfB1ZoDahMV7dlH6X52UwpyUv4nFHZYU6eOEZ1Csko/SYKMwsDdwGXAXOAa81sTo9mi4Cou88AvgfcEZw7B1gInATMB+42s3B/fZpZNZD4r2kZ5P2nnMC7ZpTyncfXs/tgS6rDkQGoa4gyd2oxZr0/aBdPTWUJL29rpLlNa4BJZkhkRlED1Lv7RndvBZYAC3q0WQDcH2w/DFxssf/3LACWuHuLu28C6oP+eu0zSCLfBr40sKGlJzPj61ecRFNrB3f87rVUhyPHae+hVjbtPnRMl526zKsoobWjk5e27huCyEQGXyKJYiKwtdvnbcG+uG3cvR1oBEr7OLevPm8Clrr7sF3zYkZ5IYvOr+RXq7epqJmhuh60O5ZCdpd5FbFztO6TZIq0Kmab2QTgauAHCbS90cxqzax2167Me+L5cxdVccLoXG759ataMDADrd4SJRIyTp005pjPLcrLZta4QlbqCW3JEIkkiu3A5G6fJwX74rYxswgwBtjTx7m97T8DmAHUm9lmIM/M6uMF5e73uHu1u1eXlZUlMIz0kp8T4SuXz2bNm/v5xYotqQ5HjlFdQ5Q5E0aTmxU+rvPnVRZTtyWqu98kIySSKFYBVWZWaWbZxIrTS3u0WQpcH2xfBTzpsRXwlgILg7uiKoEqYGVvfbr7b939BHevcPcK4HBQIB+WPnDqeM6dXsq3l61njwrbGaO9o5OXtjYeV32iy7yKEg62tLPurQODGJnI0Og3UQQ1h5uAZcA64CF3X2Nmt5rZFUGze4HS4Lf/LwCLg3PXAA8Ba4HHgM+6e0dvfQ7u0NJfV2H7cGsH33psfarDkQS99vYBmto6mHsc9YkuNZXBi4xUp5AM0P8CNYC7Pwo82mPfLd22m4nVFuKdextwWyJ9xmlTkEh8maxqXCGfeFcFP356EwtrJvf5Ok1JD3VHVoxNfOmOnsaPGcXkklGs2rSXRedVDlZoIkMirYrZI9XfXDKT8sIcbvn1GhW2M8DqLVHGjc5hYtGoAfUTe5HRXr2nRNKeEkUaKAgK269sb+SXK/XEdrqra4gyd8qxP2jXU01FCXsOtbJx96FBikxkaChRpIkrTpvAWZUlfHvZeqKHWlMdjvRi54Fmtu5tOq7nJ3qa11Wn0HIekuaUKNKEmfH1BSdxoLmNf33i9VSHI72o2xJ7mnowaknTxuYztiBbCwRK2lOiSCMnnjCa686aws9XNPD6Dt02mY5eaIiSHQ5x8sTRA+7LzKieWqI7nyTtKVGkmS+8dxZ52WG+8Zu1KnKmodVbopw8cTQ5keN70K6nmsoStkWbeKuxaVD6ExkKShRppiQ/m7+5uIqnN+zmD+t3pjoc6aa1vZOXtw/sQbuealSnkAygRJGGPn5OBdPG5vNPv1lHa7uWeEgXa95spLW9c1AK2V1mjx9NQU5ECwRKWlOiSEPZkRBf/cBsNu4+xE+f25zqcCRQ1xArZA/kieyewiFj7tRizSgkrSlRpKkLZ5Xz7pllfH/5Bvbqdtm0UNcQZWLRKMaNzh3Ufmsqinl9x0HdFi1pS4kiTZkZ/3D5bA63dvDd32sdqHRQtyU6qLOJLjWVpQDU6t0kkqaUKNJY1bhCPnrWFH6xooHX3t6f6nBGtDf3NfFWY/OA1nfqzamTxpAdDqlOIWlLiSLNff6SmRTmZul22RSrG8Ab7fqTmxXmtMljVKeQtKVEkeaK87P520uq+FP9Hp5Yp9tlU6Vuyz5ys0LMHj/wB+3imVdRwqvbGznc2j4k/YsMhBJFBvjI2VOZUV7Abb9dS0t7R6rDGZHqGqKcOrGIrPDQ/F9mXmUJ7Z3OC8GdVSLpRIkiA2SFQ3z18tls3nOY+5/dnOpwRpzmtg7WvNk4JIXsLmdOLSZkevBO0pMSRYZ4z6xyLpxVxg+W17Nbr01Nqle3N9LW4UNSyO4yOjeL2eNHq6AtaUmJIoN85fI5NLV18C+Pa3XZZFod3LY6lDMKiNUpXmjYR1uHnsaX9KJEkUFmlBfwsXOmsmRVA69ub0x1OCNGXUOUqaV5jC3IGdLvqaksoamtQ/+2knaUKDLM5y+ZSUleNl9buka3yyaBu1PXsG9QFwLszbwKLRAo6UmJIsOMGZXFl+bPonZLlP9+cXuqwxn2tkWb2HWgZcgvOwGUFeZQOTZfdQpJO0oUGejqMydz2qQxfPPR1zjYovvuh1LXg3ZDWcjurqaihFWbo3R2arYo6UOJIgOFQsbXF5zMzgMt/GD5hlSHM6zVbYmSlx1m1rjCpHzfvMoSGpva2LDzYFK+TyQRShQZ6vTJRfyv6knc+8wm6vVDZcisbohy+uQiIkP0oF1PNV11Cl1+kjSiRJHBvjT/REZlh/n6/6iwPRQOt7az7q0DSSlkd5lcMopxo3NYpYK2pJGEEoWZzTez9WZWb2aL4xzPMbMHg+MrzKyi27Gbg/3rzezS/vo0sweC/a+a2X1mljWwIQ5fYwty+NtLZvL0ht08vnZHqsMZdl7a2khHpw/JQoC9MTPmVZSwctNeJX9JG/0mCjMLA3cBlwFzgGvNbE6PZouAqLvPAL4H3BGcOwdYCJwEzAfuNrNwP30+AJwInAKMAm4Y0AiHuY+fM5VZ4wr5xm/W0tymdaAGU1ch+4wkFbK7nFVZwtv7m9kWbUrq94r0JpEZRQ1Q7+4b3b0VWAIs6NFmAXB/sP0wcLGZWbB/ibu3uPsmoD7or9c+3f1RDwArgUkDG+LwFgmH+NoVJ7Et2sSP/t8bqQ5nWHmhIcq0snyK8rKT+r3zKvU8haSXRBLFRGBrt8/bgn1x27h7O9AIlPZxbr99BpecPgY8lkCMI9o500u5/NTx/PCPb7B17+FUhzMsdD1od2YS6xNdZpYXMmZUlp6nkLSRzsXsu4Gn3P3peAfN7EYzqzWz2l27diU5tPTzlffPJmTGbb9dl+pQhoXNew6z91BrUh606ykUMqqnFuvOJ0kbiSSK7cDkbp8nBfvitjGzCDAG2NPHuX32aWb/CJQBX+gtKHe/x92r3b26rKwsgWEMbxOKRnHTRTN4bM3bPL1BiXOguhYCTGYhu7t5lSVs3HWIXQe0UrCkXiKJYhVQZWaVZpZNrDi9tEebpcD1wfZVwJNBjWEpsDC4K6oSqCJWd+i1TzO7AbgUuNbdtYzmMbjh/EqmlubxtaVraG3X/3QDUdcQpTA3woyygpR8f01Qp6jVrELSQL+JIqg53AQsA9YBD7n7GjO71cyuCJrdC5SaWT2xWcDi4Nw1wEPAWmK1hs+6e0dvfQZ9/QgYBzxnZi+a2S2DNNZhLycS5pYPzOGNXYf0gqMBqtsSe9AuFLKUfP/JE8aQmxXS5SdJC5FEGrn7o8CjPfbd0m27Gbi6l3NvA25LpM9gf0IxSXwXzx7HRSeW8/3lG1hw+gTKR+emOqSMc6C5jfU7DjD/5BNSFkN2JMQZk4tV0Ja0kM7FbDlOt3xgDq3tndz+2GupDiUjvbS1EXeS+kR2PPMqS1j75n4ONLelNA4RJYphqGJsPjecX8l/1W1n8X++zIYdB1IdUkZZvSWKGZye5AfteqqpKKHToa5hX0rjENFlnmHqry+q4kBzO79avZUlq7ZywcwyPnX+NN41o5TYs5DSm7qGKDPLCxmdm9rVY86YUkQ4ZKzctIcLZurOPkkdzSiGqVHZYb7xoZN5dvHF/N37ZrLmzf189N4VXPb9p/lV7VZa2rXcRzydnc4LDVHmTk3tbAIgPyfCyRPHsGpTNNWhyAinRDHMleRnc9NFVfxp8YV8+6pTAfjiwy9z3h1/4M4nNxA91JriCNPLG7sOsr+5PeX1iS41FcW8uG2fEruklBLFCJETCXN19WR+9zfn8/NFZzFn/Gi+8/jrnHP7cr7yyCu8sUvvtIBub7RL0YN2Pc2rKKG1vZOXtzWmOhQZwVSjGGHMjPOqxnJe1Vhe33GA+57ZxK9Wb+OBFQ18+IyJfONDJ5OfM3L/s1i9JUpRXhbTxuanOhQgliggtkBg17ZIsmlGMYLNHFfI7VeeyrOLL+Kv3jOd/35xOx+++1k27T6U6tBSpq5hH3OnFKdNwb84P5uq8gI9TyEppUQhjC3I4cvzT+SnnzyLnQeaueLOZ1i+buS9CKnxcBv1Ow8yN8W3xfZUU1nC6s1ROjr1IiNJDSUKOeK8qrEsvek8ppbmsej+Wr77+9fpHEE/nOq2pld9oktNZQkHWtpZ99b+VIciI5QShbzD5JI8Hv70uVw5dxL/tnwDi+5fRePhkfFk8AtbooQMTpuUXjOKrtqELj9JqihRyFFys8J85+pT+caCk3h6w26uuOsZXnt7+P82u7ohyuzxo9OumD+haBQTi0YpUUjKKFFIXGbGx86p4MH/fTZNrR38xV3PsvSlNwfU55v7mvjZc5u5/r6VnPKPy/j8khdo2JMeb+Tr6HReDArZ6aimsoSVm6LEVu8XSa70+tVJ0s6ZU0v4zV+fx2ceqONzv3yBl7buY/FlJ5IV7v93jM5O55XtjTyxbgdPrNt55Bp7RWkeF55YzmNr3ua3r7zFdTVTuOmiKsoKc4Z6OL16fccBDrV2pMUT2fHMqyjhkRe2s2n3Iaal6B0ZMnIpUUi/ykfn8otPnc1tv13Lvc9s4tXtjdx53dy4P9ibWjt4pn43y9ftYPlrO9l1oIWQxd4Ud/NlJ3Lx7HFML8vHzNixv5nvL9/Az1c08KvV21h0XiWfeve0lKyxdOSNdlPB6QF0AAAKJklEQVTS81mFrhcZrdq8V4lCkk6JQhKSHQnx9QUnc+qkIv7+kVf44A+e4YcfncsZU4p5u7GZ5a/tYPm6nfypfjct7Z0U5ES4YGYZF88u58JZ5RTnZx/V57jRufzzX5zCp86fxr88vp4fPFnPz57fwmffM4OPnTOV3Kxw0sZX1xBlbEE2k0tGJe07j8X0snxK87NZuSnKNfOmpDocGWGUKOSYXHnmJGadUMinf76aa/79eWaUF7A2uKQ0uWQU19ZM4ZLZ46ipLCE7klgJrHJsPndeN5dPX9DIt5at57ZH13Hfnzbx+UuquHLuJCIJXOYaqLotUc5IowftejIzqiv0IiNJDSUKOWYnTxzD/9x0Hl/99avsaGzmi5fO4pLZ45g5rmBAP2hPnjiGn36yhmff2M23HlvPl//zFe55aiNfvHQWl550wpD9EN9zsIXNew6zsCa9f1OfV1HCsjU72LG/mXF6c6EkkRKFHJfi/Gzuum7ukPR97vSxPPKZUh5fu4NvL1vPp39ex2mTi/jypbM4d8bYQf++F4IXA6XrHU9duuoUKzft5YOnTUhxNDKS6PZYSUtmxqUnncCyz7+bb191Krv2N3PdT1bwsXtX8Mogr6S6uiFKJGScOmnMoPY72OaMH01+dpiVm3T5SZJLMwpJa+GQcXX1ZD542gR+/vwW7vpDPR+88xkuP2U8Hzl7CjkJ1kH68mz9bk6aMDqpxfPjEQmHmDtVdQpJPiUKyQi5WWFuOH8a18ybzI+f3sRPnt7Ib195a9D6/9T5lYPW11CqqSjhu0+8TuPhNsbkpfZVrTJyKFFIRinMzeIL753J9edMZc2bg7OsSMiMM9JsxdjezKsswR1qt+zl4tnjUh2OjBBKFJKRSgtyePfMslSHkXSnTy4iK2ys3KxEIcmjYrZIBsnNCnPapCIVtCWpEppRmNl84PtAGPiJu9/e43gO8FPgTGAPcI27bw6O3QwsAjqAz7n7sr76NLNKYAlQCqwGPuburQMbpsjwMa+yhB8/tZGm1g5GZR9fAd7daetw2js7Y393dNLe6bR1dNL+jv2x7VQZlR2mMDeLgpwIBTkRwqH0fCByuOs3UZhZGLgLeC+wDVhlZkvdfW23ZouAqLvPMLOFwB3ANWY2B1gInARMAJ4ws5nBOb31eQfwPXdfYmY/Cvr+4WAMVmQ4qKko4Yd/fIPr71tJTlboyA/3ts7gB36H09YZ/MDv6Iy/P0NfSJXflThyIxTmRijMzaIwJ7ZdkBN557Hgc2FupNu+LHKzQmn7BH66SmRGUQPUu/tGADNbAiwAuieKBcDXgu2HgTst9i+xAFji7i3AJjOrD/ojXp9mtg64CLguaHN/0K8ShUjgrGklvGtGKY1NbbR2hMgOh8iOhMgLh8gKGZGwETmyHSIrbERCISJhIyscIhLszw7aRULB/rCRFbTrfn4kZJCKn6sOTW0dHGhu40BzOwea2znY0n7k88GWdhqb2tgePXzk8+HWjn67jYTsSOIoyIklktH9JJqcSJhQKHbjQzhkR/4OmxEK0W07+Lu3/aE/n2f2zuPpLJFEMRHY2u3zNuCs3tq4e7uZNRK7dDQReL7HuROD7Xh9lgL73L09TnsRAfKyIzxww9mpDiMttXd0BsnknYnlYEs7+5uD7ThJ5819ze/4nIoZ1zsST7ekErKuxMRRySYUMr754VOOvAVxqGTsXU9mdiNwI8CUKem9Ro+IJEckHKIoL5uivKNXK06Uu9PS3sn+bkmlpb2Tjk7H3elwp6PT6XSno5Nu293/jr2PpaPH/u7ndbofafPnthzV9kj7o9rGjuUdZ53qWCSSKLYDk7t9nhTsi9dmm5lFgDHEitp9nRtv/x6gyMwiwawi3ncB4O73APcAVFdXZ+YFVxFJO2ZGblaY3Kww5YWpjiY9JHJ77CqgyswqzSybWHF6aY82S4Hrg+2rgCc99s7GpcBCM8sJ7maqAlb21mdwzh+CPgj6/PXxD09ERAaq3xlFUHO4CVhG7FbW+9x9jZndCtS6+1LgXuBnQbF6L7Ef/ATtHiJW+G4HPuvuHQDx+gy+8svAEjP7J+CFoG8REUkRGw4va6+urvba2tpUhyEiklHMbLW7V/fXTk9mi4hIn5QoRESkT0oUIiLSJyUKERHpkxKFiIj0aVjc9WRmu4AtqY4jQWOB3akOIsk05pFBY848U9293xe7DItEkUnMrDaR29GGE415ZNCYhy9dehIRkT4pUYiISJ+UKJLvnlQHkAIa88igMQ9TqlGIiEifNKMQEZE+KVEMEjObb2brzazezBbHOZ5jZg8Gx1eYWUWw/71mttrMXgn+vijZsR+v4x1zt+NTzOygmf1dsmIeDAMZt5mdambPmdma4N88N5mxH68B/PedZWb3B2NdZ2Y3Jzv245XAmN9tZnVm1m5mV/U4dr2ZbQj+XN/z3Izj7vozwD/Elkp/A5gGZAMvAXN6tPkM8KNgeyHwYLB9BjAh2D4Z2J7q8Qz1mLsdfxj4FfB3qR5Pkv6tI8DLwGnB51IgnOoxDfGYrwOWBNt5wGagItVjGqQxVwCnAj8Fruq2vwTYGPxdHGwXp3pMA/mjGcXgqAHq3X2ju7cCS4AFPdosAO4Pth8GLjYzc/cX3P3NYP8aYJSZ5SQl6oE57jEDmNmHgE3ExpxJBjLu9wEvu/tLAO6+x4P3s6S5gYzZgfzgzZejgFZgf3LCHpB+x+zum939ZaCzx7mXAr93973uHgV+D8xPRtBDRYlicEwEtnb7vC3YF7eNx17z2kjsN8rurgTq3L1liOIcTMc9ZjMrIPaCqq8nIc7BNpB/65mAm9my4JLFl5IQ72AYyJgfBg4BbwENwHfcfe9QBzwIEhnzUJyblhJ5Z7YkgZmdBNxB7LfO4e5rwPfc/WAwwRgpIsB5wDzgMLA8eHHM8tSGNaRqgA5gArHLME+b2RPuvjG1Ycmx0IxicGwHJnf7PCnYF7dNMA0fA+wJPk8CHgE+7u5vDHm0g2MgYz4L+JaZbQY+D/x98GrcTDCQcW8DnnL33e5+GHgUmDvkEQ/cQMZ8HfCYu7e5+07gT0AmLHmRyJiH4ty0pEQxOFYBVWZWaWbZxIp5S3u0WQp03f1wFfCku7uZFQG/BRa7+5+SFvHAHfeY3f18d69w9wrgX4F/dvc7kxX4AB33uIm9I/4UM8sLfpheQOx98uluIGNuAC4CMLN84GzgtaREPTCJjLk3y4D3mVmxmRUTu0qwbIjiTI5UV9OHyx/g/cDrxO6U+Eqw71bgimA7l9gdPvXASmBasP+rxK7hvtjtT3mqxzOUY+7Rx9fIoLueBjpu4KPECvivAt9K9ViGesxAQbB/DbGk+MVUj2UQxzyP2CzxELHZ05pu534y+N+iHvhEqscy0D96MltERPqkS08iItInJQoREemTEoWIiPRJiUJERPqkRCEiIn1SohARkT4pUYiISJ+UKEREpE//HwZiLesJnyjNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(lagrange_val, rel_error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13a4e43c8>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4XNV97vHvbzTSjGRJ44tkNNgGETC3YsLFmIaQhABJgaTQBMIl5QRaTkmb0obTQAq5hzwl154mbUlS94RDEp6EOCRpXWLiJoSEhAMEGYzBgIkxBssYZixfZmRJo9s6f8xseXSxNZbmsvfM+3kePczsvTVaG8zrpbXXby1zziEiItUlVOkGiIhI8SncRUSqkMJdRKQKKdxFRKqQwl1EpAop3EVEqlBFw93M7jSzhJk9U4TPeruZrc/7GjCzPylGO0VEgsYqOc/dzN4K9ALfcc6dVMTPnQ9sBhY75/qK9bkiIkFR0Z67c+4hYFf+MTM72sx+ZmbrzOw3Znb8DD76MuB+BbuI1Co/jrmvBP7GOXc6cBPw9Rl8xpXA94vaKhGRAAlXugH5zKwZOAv4oZl5hyO5c+8Fbpvi27Y75/4o7zPiwDJgbWlbKyLiX74Kd7K/Sexxzp0y8YRz7sfAjwv4jMuBnzjnhordOBGRoPDVsIxzLgW8ZGbvA7CsNx7ix1yFhmREpMZVeirk94FHgOPMrNvMrgP+FLjOzJ4CNgKXHMLndQJLgF8Xv7UiIsFR0amQIiJSGr4alhERkeKo2APVtrY219nZWakfLzLJi8leQmYc1Tan0k0ROaB169btdM61T3ddxcK9s7OTrq6uSv14kUnOvP0XNDWEefCmcyrdFJEDMrOXC7nOb1MhRSpiZNSxs3eQ+rohnHPk1VmIBJLG3EWAXfsGGRl1DAyNsrdfJRISfAp3ESCRHhh7vWPvwEGuFAkGhbsIkEhnxl6/pnCXKqBwFwGSqf3hrp67VAOFuwiQ7M2Guxm8tre/wq0RmT3NlhEBEqkBWqNhGhvq1HOXqqBwFyE75r6wNcqcSJjXUgp3CT4Ny4iQC/eWCPHWqHruUhUU7iJkp0IubInQEYtqtoxUBYW71DznHIlUdlgmHovSmxkmPaBCJgk2hbvUvNTAMJnhUdqbsz130Fx3CT6Fu9S8ZK46dWFrhHisEdBcdwk+hbvUvESugKm9JUJcPXepEgp3qXne0gMLW6IsbI0A6rlL8CncpeYlvXBvjRAJ19HW3MBrKVWpSrAp3KXmJdIDROtDtESyNX0dMc11l+BTuEvNyxYwRcc26OhobdSYuwSewl1qXiKVrU71xNVzlyqgcJeal0gPjD1IheywzN7+IfoGhyvYKpHZUbhLzUukM7Q3j++5g6ZDSrAp3KWmDQyNkB4YZmFrdOyYqlSlGijcpablFzB5VKUq1UDhLjXN2xg7/4FqR64Xr3XdJcgU7lLTknnVqZ7GhjrmNtWzQ9vtSYAp3KWmJfKqU/N1tGpddwk2hbvUtER6gHDImN/UMO645rpL0CncpaYlUhnamiOEQjbueEdMVaoSbAp3qWnZjbEjk47HY1F69g0yMDRSgVaJzN604W5md5pZwsyeOcB5M7N/NrPNZrbBzE4rfjNFSmNiAZPHm+vuTZUUCZpCeu53ARcc5PyFwNLc1/XAN2bfLJHySE5YesDjValqxowE1bTh7px7CNh1kEsuAb7jsh4F5ppZvFgNFCmV4ZFRevYN0p43DdIztgSB5rpLQBVjzH0RsC3vfXfu2CRmdr2ZdZlZVzKZLMKPFpm5nb2DODe+gMnToSpVCbiyPlB1zq10zi13zi1vb28v548WmWSq6lRPcyRMSySsGTMSWMUI9+3Akrz3i3PHRHxt//Z6k4dlwNuRSWPuEkzFCPfVwAdys2b+ENjrnNtRhM8VKan9G2NP7rlDNtzVc5egCk93gZl9HzgHaDOzbuDTQD2Ac+6bwBrgImAz0Af8WakaK1JM3jTHtimmQkL2oeoLr6fL2SSRopk23J1zV01z3gF/XbQWiZRJIj3A/DkNNISn/gW2I9ZIIp1haGSU+jrV+0mw6E+s1KwDFTB54rEozu0fmxcJEoW71KwDLT3g6RgrZNK4uwSPwl1qVjI1MG4Hpom0l6oEmcJdapJzjmRvZtwmHRPFW71CJk2HlOBRuEtN2t03xNCIO+A0SIDWxjCN9XXquUsgKdylJo1Vpx5kzN3Mspt2aH0ZCSCFu9SkqfZOnYoKmSSoFO5Sk7wCpoMNy4DCXYJL4S416UAbY08Uj0V5PTXAyKgrR7NEikbhLjUpkR6gORKmqeHgRdodsUaGRx09vSpkkmBRuEtNSqQzB53j7om3qpBJgknhLjUpmSos3FWlKkGlcJealEgPTPswFfKrVFXIJMGicJealEgfvDrVM39OAw11Ic11l8BRuEvN6c0M0zc4Mu1MGcgWMmk6pASRwl1qTnKaHZgmym63p3CXYFG4S81JpLyNsacfloHsuLt67hI0CnepOYUWMHm8YZnspmMiwaBwl5oz3cbYE8VbowyOjLJr32ApmyVSVAp3qTmJ9AANdSFijfUFXd8R89Z119CMBIfCXWqOV8BkZgVdrx2ZJIgU7lJzCl16wOOFu+a6S5Ao3KXmFFqd6lnQHCEcMlWpSqAo3KXmJNKZgmfKANSFjMNaNdddgkXhLjUlMzzCnr6hgue4e1SlKkGjcJeasrM3O53xUIZlQOEuwaNwl5oyVp16CMMykJ3rvkOFTBIgCnepKYkCN8aeqCMWpX9ohFT/cCmaJVJ0BYW7mV1gZpvMbLOZ3TLF+SPM7EEze9LMNpjZRcVvqsjseeF+KFMhAeJeIVNKM2YkGKYNdzOrA+4ALgROBK4ysxMnXPYJYJVz7lTgSuDrxW6oSDEkUwOYwYI5DYf0fdqRSYKmkJ77CmCzc26Lc24QuAe4ZMI1DmjNvY4BrxaviSLFk0hnWDAnQrju0EYkVaUqQVPIn/BFwLa89925Y/k+A1xtZt3AGuBvpvogM7vezLrMrCuZTM6guSKzk92B6dCGZCA7jBMy9dwlOIr1QPUq4C7n3GLgIuC7Zjbps51zK51zy51zy9vb24v0o0UKl0gPHPJMGYD6uhDtLRFVqUpgFBLu24Elee8X547luw5YBeCcewSIAm3FaKBIMSVSM+u5Q3Z1SPXcJSgKCffHgaVmdpSZNZB9YLp6wjWvAOcBmNkJZMNd4y7iKyOjjp59g4c8DdITb1UhkwTHtOHunBsGbgDWAs+RnRWz0cxuM7OLc5d9BPgLM3sK+D5wrVO1h/jMrn2DjIy6GQ3LgKpUJVjChVzknFtD9kFp/rFP5b1+FnhzcZsmUlyJtLd36szCPR6Lks4Mkx4YoiVa2EYfIpWiClWpGTMtYPJ4c91f17ruEgAKd6kZydTMlh7wxLXdngSIwl1qhjcsM9Oee1xVqhIgCnepGYl0htZomGh93Yy+33sQq4eqEgQKd6kZiVSGha0zG5IBiITraGtuUM9dAkHhLjXjUPdOnUp2OqSqVMX/FO5SM2a6rky+jlZVqUowKNylJjjnSKZnNywD2Yeqr2kqpASAwl1qQmpgmMzwaFGGZfb0DdE/OFKklomUhsJdakJyltMgPWPruqv3Lj6ncJeakEjNrjrVs39HJj1UFX9TuEtNmOnG2BN5Vaqa6y5+p3CXmjC2aNgMV4T0dLSqSlWCQeEuNSGRyhCtD9ESKWgh1ANqbKhjblO9eu7iewp3qQnZOe5RzGzWn9XRGlXPXXxP4S41oRjVqZ7sXHc9UBV/U7hLTUikM7Meb/d0xBo1LCO+p3CXmpDMDcsUQzwWZWfvIJlhFTKJfyncpeoNDI2QHhie9Rx3jzfX3Zs7L+JHCnepesUqYPJo0w4JAoW7VL3Zbow9UVxVqhIACnepesWqTvV0qEpVAkDhLlUvkSpOdaqnORKmJRLWsIz4msJdql4inSEcMuY3NRTtM7M7Mincxb8U7lL1EukMbc0RQqHZV6d6OmJRdmjZX/ExhbtUvWIWMHni2ktVfE7hLlUvWYS9UyfqiDWSSGcYGhkt6ueKFIvCXapeMj1Ae5FmynjisSjOZf/iEPGjgsLdzC4ws01mttnMbjnANZeb2bNmttHMvlfcZorMzPDIKD37BotWwOTpUCGT+Ny0i1ubWR1wB/AOoBt43MxWO+eezbtmKXAr8Gbn3G4zW1iqBoscip29gzhXvAImz9heqgp38alCeu4rgM3OuS3OuUHgHuCSCdf8BXCHc243gHMuUdxmisxMsatTPfHWbCGTqlTFrwoJ90XAtrz33blj+Y4FjjWzh83sUTO7YKoPMrPrzazLzLqSyeTMWixyCLx1ZRa2FnfMvbUxTGN9nXru4lvFeqAaBpYC5wBXAf9uZnMnXuScW+mcW+6cW97e3l6kHy1yYPuXHihuz93MiGuuu/hYIeG+HViS935x7li+bmC1c27IOfcS8ALZsBepKG9Ypq25uOEOqlIVfysk3B8HlprZUWbWAFwJrJ5wzX+Q7bVjZm1kh2m2FLGdIjOSSGeYP6eBhnDxZ/0q3MXPpv0T75wbBm4A1gLPAauccxvN7DYzuzh32Vqgx8yeBR4EbnbO9ZSq0SKFSqSKX8DkiceivJ4aYGTUleTzRWZj2qmQAM65NcCaCcc+lffaAX+X+xLxjWRvpuhz3D0dsUaGRx09vZmiP7AVmS1VqEpVS6YGShbu8VYVMol/KdylajnnSPYWb2PsiVSlKn6mcJeqtbtviKERV9Ixd0CrQ4ovKdylao1VpxZ5uV/P/DkNNNSFNNddfEnhLlVrrDq1RMMyZqbpkOJbCnepWqWqTs3XEYtqzF18SeEuVavUwzLg7cikcBf/UbhL1UqkMjRHwjQ1FFTOMSPesEy21EPEPxTuUrWy0yBL12uH7Fz3wZFRdu0bLOnPETlUCnepWslUhrYSh3tHzFvXXUMz4i8Kd6laifRA6Xvu2pFJfErhLlUrkS5ddarHC3fNdRe/UbhLVerNDNM3OFLSmTIAC5ojhEOmKlXxHYW7VKVEqjR7p05UFzIOa9Vcd/EfhbtUpf0FTKVfildVquJHCnepSmPhXuJhGVC4iz8p3KUqlWtYBrJz3XeokEl8RuEuVSmZztAQDhFrrC/5z+qIRekfGiHVP1zynyVSKIW7VKVkOkN7cwQzK/nPinuFTCnNmBH/ULhLVUqkS7d36kTakUn8SOEuVakc1akeVamKHyncpSol0pmyzJQBaG+JEDL13MVfFO5SdTLDI+zpGyrLHHeA+roQ7S0RVamKryjcpeoky7AD00QdsUb13MVXFO5SdcpZwOSJt6qQSfxF4S5Vp9QbY09FVariNwp3qTrJdPmqUz3xWJR0Zpj0wFDZfqbIwRQU7mZ2gZltMrPNZnbLQa671MycmS0vXhNFDk0inSFk2eV4y8Wb6/661nUXn5g23M2sDrgDuBA4EbjKzE6c4roW4MPAY8VupMihSKYzzJ8ToS5U+upUT1zb7YnPFNJzXwFsds5tcc4NAvcAl0xx3eeALwL60y0Vld2BqXy9dsjbkUnhLj5RSLgvArblve/OHRtjZqcBS5xzPz3YB5nZ9WbWZWZdyWTykBsrUohEeqCsM2Vg/8wcPVQVv5j1A1UzCwH/G/jIdNc651Y655Y755a3t7fP9keLTCmRKn/PPRKuo625QT138Y1Cwn07sCTv/eLcMU8LcBLwKzPbCvwhsFoPVaUSRkYdO3tLvzH2VLLTIVWlKv5QSLg/Diw1s6PMrAG4EljtnXTO7XXOtTnnOp1zncCjwMXOua6StFjkIHr2ZRh15S1g8nS0qkpV/GPacHfODQM3AGuB54BVzrmNZnabmV1c6gaKHIr9BUzlD/d4LMprmgopPhEu5CLn3BpgzYRjnzrAtefMvlkiM+OtK9NeoWGZPX1D9A+O0NhQV/afL5JPFapSVRIVqE71jK3rrt67+IDCXarK/p57Bcbcx+a666GqVJ7CXapKIp2hNRomWl/+YRGvSlVz3cUPFO5SVRKpDAtbyz/eDtDRqipV8Q+Fu1SVcu6dOlFjQx1zm+rVcxdfULhLVanEujL5Olqj6rmLLyjcpWo453IbY1dmWAa8ue56oCqVp3CXqpHqH2ZweLSyPfdYo4ZlxBcU7lI1vDnulZgG6YnHouzsHSQzPFKxNoiAwl2qyNjG2BWoTvV4c929ZRBEKkXhLlXDLz130HRIqTyFu1QNrzq1EitCeuKqUhWfULhL1UikMkTrQ7RECloPryQ6VKUqPqFwl6qRneMexax8G2NP1BwJ0xIJa1hGKk7hLlWjktWp+bI7MincpbIU7lI1sgVM/gj3HVr2VypM4S5VI5mqzN6pE8W1l6r4gMJdqkL/4AjpzHBFp0F6OmKNJNIZhkZGK90UqWEKd6kKldyBaaJ4LIpz+6dmilSCwl2qQqKCOzBN1KFCJvEBhbtUhaQPlh7wjO2lqnCXClK4S1VI5Gan+GG2TLw1W8ikKlWppECGu3Ou0k0Qn0mkM4RDxvymhko3hdbGMI31deq5S0UFLty/99grnP3FBzUTQcZJpDO0NUcIhSpXneoxM+Ka6y4VFrhwb4mG2b6nn+d2pCrdFPERvxQweVSlKpUWuHBf3jkPgMe37q5wS8RPEil/LD3gUbhLpQUu3OOxRhbNbWTdy7sq3RTxkWQ6Q7sPZsp44rEor6cGGBnV8yGpjMCFO8AZnfPo2rpbD1YFgKGRUXr2Dfqs597I8Kijp1eFTFIZBYW7mV1gZpvMbLOZ3TLF+b8zs2fNbIOZPWBmRxa/qfud3jmfRDrDtl2aaiaws7fym3RMFG9VIZNU1rThbmZ1wB3AhcCJwFVmduKEy54EljvnTgbuBb5U7IbmW35kdty9S0Mzwv79Stub/RPuqlKVSiuk574C2Oyc2+KcGwTuAS7Jv8A596Bzri/39lFgcXGbOd6xh7XQEg3roaoA+dvr+WvMHdDqkFIxhYT7ImBb3vvu3LEDuQ64f6oTZna9mXWZWVcymSy8lRPUhYzTjpinh6oC7F9Xxk9j7vPnNNBQF9Jcd6mYoj5QNbOrgeXAl6c675xb6Zxb7pxb3t7ePqufdUbnPF54vZe9fUOz+hwJPm9FyDYfDcuYmaZDSkUVEu7bgSV57xfnjo1jZucDHwcuds6VfIrA6UfOB2DdK+q917pEOpPtKYf9NfmrIxbVmLtUTCH/NzwOLDWzo8ysAbgSWJ1/gZmdCvwb2WBPFL+Zk52yZC7hkNGlcfeal0hlfDUk44mr5y4VNG24O+eGgRuAtcBzwCrn3EYzu83MLs5d9mWgGfihma03s9UH+LiiaWyo4w8WxRTuQjI94It13CfyhmVUjyGVEC7kIufcGmDNhGOfynt9fpHbVZDlR87j7kdfJjM8QiRcV4kmiA8k0hmOWdhS6WZMEm+NMjgyyq59gyzw0fMAqQ3+GqQ8RGd0ziMzPMoz2wtbROzBTQm+/f+2lrZRUlajo46kzxYN83TEvHXdNTQj5RfocB97qFrAlMhntu/lg99dx6dXb9QUyiqyu2+Q4VHnqwImj3ZkkkoKdLi3t0ToXNA07bj7nr5B/vLudSyY08BhrRFuu+85RrWgU1UYm+Puw567F+6a6y6VEOhwh2zvfd3LB15EbHTUceMP1pNIZfjG1adz8x8dz1Pb9vCfT02azSkB5Ke9Uyda0BwhHDJVqUpFBD7cz+icR8++QV7auW/K8//8y9/zq01JPn3xiZyyZC7vPXURyxbF+OL9m+gbHC5za6XY/Fid6qkLGYe1aq67VEbww/2o7Lj7zfdu4PGt48fSH3w+wdce+D2XnraY9684AoBQyPjku0/ktdQAKx/aUvb2SnF51al+HJYBbdohlRP4cD+6vZkvvHcZr+zq433ffIRr7vwdG7r38EpPHx++50lO6GjlH95zEmb799ZccdR83rUszr/9eot2qA+4RCpDcyRMU0NBs3rLTuEulRL4cAe4csURPHTz27n1wuN5qnsPF//rw1xyx28B+ObVpxOtnzwH/pYLj2dk1PHln20qd3OliJJpf1aneuK5YRkVMkm5VUW4Q7Zi9YNvO5rffPTt/K/zj6U5GuZrV53KEQuaprx+yfwm/vzso/jxk9t5atueMrdWiiXh0+pUT0csSv/QCKl+Pd+R8qqacPe0ROv58PlL+c1Hz+Xtxy086LV//fajaWtu4HP3PaueVUAl0hlfreM+UdwrZEpp+E/Kq+rC/VC0ROv5yDuPo+vl3fz06R2Vbo4cIucciVTGlwVMHu3IJJVS0+EOcPnyJRzf0cLn1zzPo1t61IMPkN7MMP1DI76dKQP7C5leV7hLmdV8uNeFjH94z0n0Zoa5cuWjvOOfHuL/PvwSe/u1CYjfJX08x93T3hIhZOq5S/nVfLhDtsr1sY+dx5cvO5k5kTCf/a9nOfP2X/DRe59iQ7cetvpVwsfVqZ76uhDtLRFNh5Sy8+fk4AqI1tfxvuVLeN/yJTyzfS93P/oy/7n+VVZ1dXPy4hhXnLGEi994OC3R+ko3VXL8vK5Mvo5Yo9aXkbJTz30KJy2K8YVLT+axj5/HZy/+AwaGRvj4T55hxT88wEdWPcXjW3dpbN4HErnA9POwDGTnumt9GSk39dwPojVazzVndfKBNx3J+m17WNW1jdXrX+VHT3TzhvY5XL58CZeettjX86yrWTKdoSEcItbo79+mOmJRHn5xZ6WbITVG4V4AM+PUI+Zx6hHz+MS7TuSnT+9g1ePb+ML9z/OVtZs49/iFXLXiCN56bDt1IZv+A6UoEunsNMj8pSX8KB6Lkh4YpjczTHNE/8tJeehP2iGaEwlz+fIlXL58CZsTvazq2saP1nXz38++TjwW5X3Ll3DFGUtYNLex0k2teon0gO/H22H/XPfX9g5wzMLmCrdGaoXG3GfhmIXNfOyiE3jk1vP4xp+extLDWviXX/6es7/4S66583f87JkdDI2MVrqZVcvvBUwer0pVM2aknNRzL4KGcIgLl8W5cFmcbbv6+GHXNlZ1dfOXdz9BW3OEy05fzBVnLOGotjmVbmpVSaQznPmG+ZVuxrTGdmTSQ1UpI4V7kS2Z38TfvfM4/va8pfz6hST3PL6Nf//NFr756xd5+3HtXHf2G3jzMQt8P07sdwNDI+ztH/L1HHePN3SknruUk8K9RMJ1Ic474TDOO+EwXk8NcM/vtvHdR7dy9bce4/iOFv787KO45JTDiYQnL0cs09vZ6//qVE8kXEdbc4PmuktZacy9DA5rjfLh85fy278/ly9ddjIAH713A2/+wi/52i9+T08uqKRwQSlg8mjTDik39dzLKFpfx+XLl/C+0xfz8OYevvXbLfzTL17gjl9t5j2nLOK6txzFsYe1VLqZgZBI+X/pgXwdrY1s36MxdykfhXsFmBlnL23j7KVtbE70cufDL/HjJ7r5Qdc23rK0jRPirbREwjRHw7RE62mOhGmNjn/fEg1PucNUrUimg1Gd6onHoqx7edf0F4oUicK9wo5Z2Mzt71nGTe88ju899jKrurr53Uu7yAxPP4WyIbco1bJFMZYtjmX/uSjGvDkNJWlrZniEhrqQLx4GJ9IZQgYLAjAVErLDMrv7hhgYGqnpv5SlfBTuPjF/TgM3nLuUG85dCsDg8Ci9mWF6B4ZJZ4ayFY65170Dw6RyFY/du/t5ZvtefrbxtbHPWjK/MRf0czl5cYyTDo8Razp4ib5zjp59g7y6p5/tu/vZvif79eqefl7dM8D2Pf3s2jdIY30dRy5oonPBHI5sa+KoBXM4csEcOtuaOKwlSqhMFbqJVIYFzZHAVATH8wqZOjUlVsqgoHA3swuArwF1wP9xzn1hwvkI8B3gdKAHuMI5t7W4Ta0tDeEQ88MNzC+wF763f4iN2/eyYftent6+l6e797Lm6f2Bf+SCJpYtinHy4hjzmhrYsXeA7bv7eXXv/jCf+NtCU0Mdi+Y2cvjcRpYtjtHRGmVP3xAv9+zj94k0Dzz/OkMj+xdQi9aHOHJ+Nug780K/c8EcOlqLG/yJ9EAgCpg8+TsyKdylHKYNdzOrA+4A3gF0A4+b2Wrn3LN5l10H7HbOHWNmVwJfBK4oRYNlarHGes46po2zjmkbO7anbzAb9Lmwf/KVPdy3Yf92gu0tEQ6f28gJ8VbOO2HhWJAvmtfIormNxBrrDzoEMzLqeHVPP1t79rG1p4+Xd+5ja88+Xkzu48HnkwzmVedGwiGOXNDEYbn9Tp2DUedyX9nfHEbHjnnvHaOj2WP51zsH2/f086ajF5Tg32RpeFWq23b1ceoRc3EOHNm/GL0FRr2/Jp1zea8Zd9Lhprzee+/yrpvSoR3mQIufHvDzcwzDDCz7BsMIWfZ5k0HuXPYCMwhNOO79sRt/Lu8aHwwN+l0hPfcVwGbn3BYAM7sHuATID/dLgM/kXt8L/KuZmSvRurjnnHMOAAsWLOBHP/oRALfeeiuPPPLIuOsWL17M3XffDcCNN97I+vXrx50/9thjWblyJQDXX389L7zwwrjzp5xyCl/96lcBuPrqq+nu7h53/k1vehOf//znAbj00kvp6ekZd/68887jk5/8JAAXXngh/f3jZ0u8+93v5qabbhp3T/kuv/xyPvShD9HX18dFF1006fy1117Ltddey86dO7nssssmnf+rv/orrrjiCjqjA6z6zEeoBxaNjDIy6mgIh7j5ppv44z8+n02bNvHBD35w0vd/4hOf4Pzzz2f9+vXceOONk87ffvvtnHXWWWx7fj3f+tjHxp2LAx+77fM0H76U++5fyz0rv8bTwyM8MTw69j/76e+/mdZ4Jzs2/JYXfvH9se/1/sd92198lpa2Dl763X/z3C9/lAuK7DkD3nNn9r/tXXfdxV133TWpfWvWrKGpqYmvf/3rrFq1atL5X/3qVwB85Stf4b777ht3rrGxkfvvvx+Az33uczzwwAPjzh/qn70nnnyS117axQe+lz1fP38RCy74GwB6fvYvDO3aPu77Gxa+gfnnXw/Azv/6CsPp8atKRhYdz7y3XQtA8ie3M9KfGnc+euQbmfvmqwB4fdWnccPjp9s2Hr2C2JnvBeC1790y6d/NnOPfQstp72J0aIAgmPQGAAAFAElEQVTEDz8z6XzzsvNpXnY+I317Sf7H5yedbzn1Iuac8FaGU0l23vePk863rngPTcecyVBPNz1r/3XS+dhZV9LYeQqDr29h1wMrJ52f+9ZriC4+gcHtz7H7oe+MHfciv/0d19MYP5q+LevZ+dvvT/r+jnf9LZEFi0m/8Ci7Hv3xpPOHX3Iz9bF2Uht/ze51P81+dt7fJ4su+zjhphh7nvo5e5/6+aTvX3LVbYTqo+zquo/0sw9NOn/nD+/jijOOmHS8mGy6/DWzy4ALnHP/M/f+fwBnOuduyLvmmdw13bn3L+au2Tnhs64Hrs+9PQ7YVKwbKbE2oNbWbNU914ZavGcI9n0f6Zxrn+6isj5Qdc6tBCb/NexzZtblnFte6XaUk+65NtTiPUNt3HchFarbgSV57xfnjk15jZmFgRjZB6siIlIBhYT748BSMzvKzBqAK4HVE65ZDVyTe30Z8MtSjbeLiMj0ph2Wcc4Nm9kNwFqyUyHvdM5tNLPbgC7n3GrgW8B3zWwzsIvsXwDVJHBDSUWge64NtXjPUAP3Pe0DVRERCR6tCikiUoUU7iIiVaimw93MLjCzTWa22cwmVXKYWcTMfpA7/5iZdeaOv8PM1pnZ07l/nlvuts/GTO877/wRZtZrZjeVq82zNZt7NrOTzewRM9uY+28eiHWGZ/Hnu97Mvp271+fM7NZyt32mCrjnt5rZE2Y2nKvhyT93jZn9Pvd1zcTvDRznXE1+kX04/CLwBqABeAo4ccI1HwK+mXt9JfCD3OtTgcNzr08Ctlf6fspx33nn7wV+CNxU6fspw3/rMLABeGPu/QKgrtL3VOJ7fj9wT+51E7AV6Kz0PRXpnjuBk8muhXVZ3vH5wJbcP+flXs+r9D3N5quWe+5jyyo45wYBb1mFfJcA3869vhc4L7eswpPOuVdzxzcCjbnF04JgxvcNYGZ/ArxE9r6DYjb3/E5gg3PuKQDnXI9zbqRM7Z6N2dyzA+bkalYagUEghf9Ne8/Oua3OuQ3AxDW1/wj4uXNul3NuN/Bz4IJyNLpUajncFwHb8t53545NeY1zbhjYS7bnlu9S4AnnXFD2ypvxfZtZM/D3wGfL0M5ims1/62MBZ2Zrc7/Of7QM7S2G2dzzvcA+YAfwCvAV51wQdhop5J5L8b2+pPXcZ8HM/oDsCpjvrHRbyuQzwD8553praFW+MHA2cAbQBzxgZuuccw8c/NsCbQUwAhxOdojiN2b2C5dbPFCCoZZ77rNaVsHMFgM/AT7gnHux5K0tntnc95nAl8xsK3Aj8LFcgZvfzeaeu4GHnHM7nXN9wBrgtJK3ePZmc8/vB37mnBtyziWAh4EgrMNSyD2X4nt9qZbDfcbLKpjZXOCnwC3OuYfL1uLimPF9O+fe4pzrdM51Al8FbnfOTV6v1X9ms4TGWmCZmTXlAvBtjF/u2q9mc8+vAOcCmNkc4A+B58vS6tkp5J4PZC3wTjObZ2bzyP42vrZE7SyPSj/RreQXcBHwAtkn7B/PHbsNuDj3Okp2Vshm4HfAG3LHP0F2THJ93tfCSt9Pqe97wmd8hoDMlpntPQNXk32A/AzwpUrfS6nvGWjOHd9I9i+ymyt9L0W85zPI/ja2j+xvKRvzvvfPc/8uNgN/Vul7me2Xlh8QEalCtTwsIyJStRTuIiJVSOEuIlKFFO4iIlVI4S4iUoUU7iIiVUjhLiJShf4/2HWAUF+wqzYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(lagrange_val, SGL_vals - regreg_vals )\n",
    "plt.gca().set_ylim([0, (SGL_vals - regreg_vals).max()])\n",
    "plt.plot([lagrange_val[0], lagrange_val[-1]], [0,0], 'k--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the path "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from regreg.paths.sparse_group_lasso import sparse_group_lasso_path\n",
    "path = sparse_group_lasso_path.gaussian(X, \n",
    "                                        Y,\n",
    "                                        groups,\n",
    "                                        np.ones(groups.shape),\n",
    "                                        weights={})\n",
    "path.saturated_loss.coef = 1. / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sanity check to make sure problem is the same\n",
    "\n",
    "Should have same $\\lambda_{\\max}$ as before..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1104770302772522"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = path.full_gradient(path.saturated_loss, np.zeros(path.saturated_loss.shape))\n",
    "dual = path.penalty.conjugate\n",
    "dual.seminorm(G, lagrange=1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.88 s ± 249 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "path.main(lagrange_val, inner_tol=1.e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = path.main(lagrange_val, inner_tol=1.e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0,\n",
       " 1.549426481226161e-05,\n",
       " 4.084087269362177e-06,\n",
       " 5.862205934936953e-06,\n",
       " 4.015184982733448e-06,\n",
       " 9.640124049407405e-06,\n",
       " 1.0358522697299568e-05,\n",
       " 1.2429342064199581e-05,\n",
       " 7.793252203431089e-06,\n",
       " 3.7945590142103634e-06,\n",
       " 1.0275565744517125e-05,\n",
       " 4.858898451058298e-06,\n",
       " 1.2162085890149148e-05,\n",
       " 1.081814475744752e-05,\n",
       " 8.845958594638587e-06,\n",
       " 7.378774311677896e-06,\n",
       " 1.0712327198010238e-05,\n",
       " 2.702667570208341e-05,\n",
       " 1.357296010832944e-05,\n",
       " 1.24959993658787e-05]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_path = results['beta']\n",
    "solns = np.array(solns)\n",
    "[np.linalg.norm(beta_path[i] - solns[i]) / np.linalg.norm(beta_path[i]) for i in range(1, solns.shape[0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 12 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 99.8        103.09038693 105.5182928  108.54692596 112.38319025\n",
      " 116.48041016 121.36895681 126.72187208 133.17474326 140.79290915\n",
      " 149.15863939 157.58746678 165.74409741 173.11402464 179.99451821\n",
      " 187.15258132 194.34681636 201.20902089 207.51896927 213.24900553] [10.30722628  9.41984798  9.48343841 10.03790105 10.73757714 11.28839041\n",
      " 11.33699888 11.52254838 11.57455493 11.86721821 12.16465021 12.39688145\n",
      " 12.67269067 12.76330497 13.23439403 13.32916723 13.70949068 14.22810287\n",
      " 14.46204655 14.34885972]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done   5 out of   5 | elapsed:   34.0s finished\n"
     ]
    }
   ],
   "source": [
    "from regreg.paths.cross_validate import cross_validate\n",
    "cv_MSE, cv_MSE_sd = cross_validate(path,\n",
    "                                   lagrange_val, \n",
    "                                   cv=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$\\\\log(\\\\lambda)$')"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX+//HXJz0hgQBJKCH00DuRXgUVbCj2gnXFAmtB17J+3fLbZS2ICK5lUZF1bbgrCioiLB2VJtJr6IQuvbfz+2MGzeJAAubOJJn38/HIgztn7kw+XId5e+4951xzziEiInK6iFAXICIihZMCQkREAlJAiIhIQAoIEREJSAEhIiIBKSBERCQgBYSIiASkgBARkYAUECIiElBUqAv4NVJSUlzVqlVDXYaISJHy/fff73DOpea1X5EOiKpVqzJnzpxQlyEiUqSY2br87KdTTCIiEpACQkREAlJAiIhIQAoIEREJSAEhIiIBKSBERCQgBYSIiASkgBARkYDCMiBu+Md33PCP70JdhohIoRaWASEiInlTQIiISEAKCBERCUgBISIiASkgREQkIAWEiIgE5FlAmFmGmU0ysyVmttjMHvK3DzCzZWa2wMw+NbPkXK95ysyyzWy5mV3iVW0iIpI3L3sQx4FHnXP1gFZAHzOrB4wHGjjnGgErgKcA/M/dCNQHugGvmVmkh/WdN82jEJFw4FlAOOc2O+fm+rf3AUuBdOfcOOfccf9uM4BK/u0ewEfOuSPOuTVANtDCq/pEROTsgnINwsyqAk2Bmac9dRfwlX87HdiQ67mN/jYREQkBzwPCzBKBT4CHnXN7c7U/je801Pvn+H69zWyOmc3Zvn17wRYrIiI/8TQgzCwaXzi875wbmav9DuBy4BbnnPM35wAZuV5eyd/2P5xzQ51zWc65rNTUVM9qFxEJd16OYjLgbWCpc+6lXO3dgMeBK51zB3O9ZDRwo5nFmlk1IBOY5VV9IiJydlEevndboBew0Mzm+dt+DwwBYoHxvgxhhnPuPufcYjP7GFiC79RTH+fcCQ/rExGRs/AsIJxz0wEL8NSYs7ymP9Dfq5pERCT/NJNaREQCUkCIiEhACogQ0ExsESkKvLxIXWg55zhx0uW9o4hIGAvLHsSeQ8f4YcNuXhq3nD2HjoW6HBGRQiksAyI2OpLkhBiGTMym/fMTeWXCSvYfOZ73C0VEwkhYBkR8dCSZaYmMebA9LauXZeD4FbR/fiJvTFnFwaMKChERCNOAOKVexZK8eVsWo/q0pVGlZJ77ahkdXpjE29PXcPiY5uiJSHgL64A4pXFGMv+8qwWf3N+aWuWS+MsXS+g4YBL/+m4tR44rKEQkPCkgcmlepQwf3NOKD+9pReUyCTwzajEXvjiFj2at59iJk6EuT0QkqBQQAbSuUZaP723Nu3e1ICUplidHLqTLwCl88v1GDY8VkbARlvMg8sPM6FArlfaZKUxcto2Xxq/g0X/P59XJ2UQAZUrEhLpEERFPqQeRBzOjS91yfN63HW/c2oyoCCN7+wEW5uxl0vJtIalJM7FFJBgUEPkUEWF0a1CBrx7qQI3UEpx0jjvfmc1tw2axfMu+UJcnIlLgFBDnKDLCSEmMpVGlUvzfZXX5Yf0uug+eytOfLmTH/iOhLk9EpMAoIM5ThBm/aV+dKb/rzG2tq/LR7A10GjCZ1yev0hwKESkWFBC/UpkSMfzpyvp8/XAHWlYrw/Njl9H1pSl8sWATP99uW0Sk6FFAFJCaaYm8fccFvHd3SxJjo+j7wQ9c+8Z3zNuwO9SliYicFwVEAWuXmcKXD7bnuZ4NWffjQa569Rse+ugHcnYfCnVpIiLnRPMgPBAZYdzYojKXN67IG5NX8ea01YxdtIV72lfnvk41Ql2eiEi+eNaDMLMMM5tkZkvMbLGZPeRvL2Nm481spf/P0v52M7MhZpZtZgvMrJlXtQVLYmwUj11Sm4mPdaJbg/L8fVI2nV+czLZ9R3R9QkQKPS9PMR0HHnXO1QNaAX3MrB7wJDDBOZcJTPA/BugOZPp/egOve1hbUKUnxzP4xqZ8+kAbKpdJYM2OAyzatJfZa3eGujQRkTPyLCCcc5udc3P92/uApUA60AP4p3+3fwJX+bd7AO86nxlAsplV8Kq+UGhauTT/ua81NVNLcPyE47o3vuPx/8xn54GjoS5NROQXgnKR2syqAk2BmUA559xm/1NbgHL+7XRgQ66XbfS3FStmRln/RLt7O1Zn5NwcLhw4mRGz13MySAsBaqkOEckPzwPCzBKBT4CHnXN7cz/nfCfiz+lb0cx6m9kcM5uzffv2Aqw0uCIjjKe61+XLB9tTKy2JJz5ZyHX/+I6lm/fm/WIRkSDwNCDMLBpfOLzvnBvpb9566tSR/89TK97lABm5Xl7J3/Y/nHNDnXNZzrms1NTU86prxL2tGXFv6/N6bUGrXT6JEfe2YsC1jViz4wCXvzKd/l8u4YDukS0iIeblKCYD3gaWOudeyvXUaOB2//btwKhc7bf5RzO1AvbkOhVVrJkZ12VlMKFfR67PqsSb09bQ9aUpjF20WaOdRCRkvOxBtAV6ARea2Tz/z6XAc8BFZrYS6Op/DDAGWA1kA28CD3hYW6FUukQMz/ZsxCf3tyE5IYb73pvLXcNns/7Hg6EuTUTCkGcT5Zxz0wE7w9NdAuzvgD5e1VOUNK9Sms/7tmX4t2sZNH4FFw2awm8vrMk9HaoTGxUZ6vJEJExoqY1CKioygt+0r86ERzvRpW4aL45bQffB0/g2e0eoSxORMKGAKOTKl4rjtVuaM/zOCzh+wnHzWzPJ3rafo8dPhro0ESnmFBBFRKfaaYx7pAMPdslk54GjLMjZwwczgzd3QkTCjwKiCImLjqTfRbVoWKkUJWIi+f2nC7npzRms2XEg1KWJSDGkgCiC4qMjqVM+ieevaciSzXu55OWpvDY5m2MngnfaSbOxRYo/Lfd9HgrDJDsz44YLKtO5dhp/HL2YF8Yu54v5m3n+mkY0rFQq1OWJSDGgHkQRl1Yyjtdvbc4btzZnx/4j9Hh1Os+OWcqho7ovtoj8OgqIYqJbg/KM79eRGy7I4B9TV9Nt8FQNiRWRX0UBUYyUio/m2Z6N+OCelhhw81szeeI/C9hz8FioSxORIkgBUQy1qZHC2Ic7cF/HGvxn7ka6DprCVwvDYlkrESlACohiKi46kie712FUn7aUKxnL/e/P5d5/zWHr3sOhLk1EiggFRDHXIL0Unz3Qlie712Hy8u10HTiFbXsPa5VYEcmTAiIMREVGcF/HGnz9cAcapJdizY8HWbplH6u37w91aSJSiCkgQiBUNyyqmlKCD+5pSbWUEhw8eoJug6fx6qTgTrATkaJDARFmzIy0pFgapZeia900Bny9nCtemc78DbtDXZqIFDIKiDAVExXBa7c0Z2iv5uw6eJSrX/uG//d58G51qqU6RAo/BUSYu7i+b4LdzS0rM+ybNVw8aCqTl2/L+4UiUuwpIISScdH89aqG/Pu+1sRFR3DHO7N5ZMQ8dh44GurSRCSEFBDykwuqlmHMQ+15sEsmXyzYRJeBk/n0h40aEisSphQQ8j9io3z3nPjit+2pmlKCR0bM5/Z3ZrNh58FQlyYiQeZZQJjZMDPbZmaLcrU1MbMZZjbPzOaYWQt/u5nZEDPLNrMFZtbMq7okf2qXT+I/97Xhz1fW5/u1O7l40FTemraaE7qDnUjY8LIHMRzodlrbC8CfnXNNgD/4HwN0BzL9P72B1z2sS/IpMsK4vU1VxvXrSOsaZfnrl0vp+do3LNm0N9SliUgQeHbDIOfcVDOrenozUNK/XQrY5N/uAbzrfCe7Z5hZsplVcM5phbkAgj3JLj05nrdvz+KLBZv50+jFXPn36aQlxZKeHB/UOkQkuIJ9DeJhYICZbQBeBJ7yt6cDG3Ltt9Hf9gtm1tt/emrO9u3bPS1WfmZmXNG4Iv/t15Grmqazac9hFuTsYZKGxIoUW8EOiPuBR5xzGcAjwNvn+gbOuaHOuSznXFZqamqBFyhnV7pEDC9e15g65ZMwgzvfmc39733P5j2HQl2aiBSwYAfE7cBI//a/gRb+7RwgI9d+lfxtUkiVio+mYXopHru4FhOXbaPrwCm8NW01x7Wuk0ixEeyA2AR09G9fCKz0b48GbvOPZmoF7NH1h8Ivwoy+F2Yy/pGOXFCtDH/9cilX/P0b5q7f5fnv1lIdIt7zcpjrh8B3QG0z22hmdwP3AAPNbD7wN3wjlgDGAKuBbOBN4AGv6pKCV7lsAu/ccQGv39KMXQeOcs3r3/LUyIXsPqiZ2CJFmZejmG46w1PNA+zrgD5e1SLeMzO6N6xA+1qpDBq/guHfrmXc4i38/tK69GyWjpmFukQROUeaSS0FKjE2imcur8fovm3JKJPAo/+ez01vziB7275QlyYi50gBIZ6oX7EUI+9vw9+ubsiSTXvpPngaL4xdxqGjJ0JdmojkkwJCPBMRYdzcsjITH+vEFY0r8trkVVw0aAoTl20NdWkikg+eXYOQwivYM7FTEmN56fomXNc8g2dGLeKu4XMonRBN1bIlglqHiJwb9SAkaFrXKMuYB9vzu0tqs/vQMRbk7OHjORu0nLhIIaWAkKCKiYqgT+eaNEovRUJ0JI//ZwG3vzObnN2aiS1S2CggJCTioiOpWyGJP19Znzlrd3LJoKm8P3OdehMihYgCQkLGzLec+NcPd6BRpVI8/ekibnlrJut/1M2JRAoDBYSEXEaZBN7/TUv+dnVDFmzcwyUvT2X4N2s4qZsTiYSUAkIKBTPfkNhxj3SgRbUy/OnzJdww9DtWb9/vye/TWk4ieVNASKFSMTme4XdewIBrG7F8yz66D57G0KmrdKtTkRBQQEihY2Zcl5XB+H4daZ+Zyt/GLOOa179l5VYt1yESTAoIKbTKlYzjzduaM/jGJqz78QCXDZnOq5OyOaZ7TogEhWZSy3kJ1mxsM6NHk3Ta1Ejhj6MXMeDr5Xy1aDMGJMTo4yviJfUgpEhITYrltVua89otzdiy5zCLcvayafchjXQS8ZACQoqUSxtWYNwjHUlOiGbDrkPcNmwW2/YdDnVZIsWSAkKKnDIlYshMS6Ra2QRmr91J95enMXn5tlCXJVLsnDUgzOzWXNttT3uur1dFieTFzEgrGcfnv21HSmIsd7wzm/5fLuHocV3AFikoefUg+uXafuW05+4q4FpEzlmtckmM6tuWXq2q8Oa0NVzz+res2XEg1GWJFAt5BYSdYTvQ4/990myYmW0zs0Wntf/WzJaZ2WIzeyFX+1Nmlm1my83sknxVL4Jv4b+/XNWAN25tzvqdB7l8yDRGzt0Y6rJEiry8AsKdYTvQ49MNB7rlbjCzzkAPoLFzrj7wor+9HnAjUN//mtfMLDKP9xf5H90alGfMQ+2pX7EU/T6eT78R89h/5HioyxIpsvIKiDpmtsDMFubaPvW49tle6JybCuw8rfl+4Dnn3BH/PqeuLPYAPnLOHXHOrQGygRbn+pcRSU+O54N7WvJw10w+m5fD5UOmsXDjngL/PVrLScJBXjON6hbw76sFtDez/sBh4DHn3GwgHZiRa7+N/jYpprycaBcVGcHDXWvRunpZHh4xj56vf8Pjl9Th7nbViIg465lREcnlrD0I59y63D/AfqAZkOJ/fK6igDJAK+B3wMdmdk7/Ys2st5nNMbM527dvP48SJFy0rO67xWmn2mn0H7OUO4fPZsf+I6EuS6TIyGuY6xdm1sC/XQFYhG/00r/M7OHz+H0bgZHOZxZwEkgBcoCMXPtV8rf9gnNuqHMuyzmXlZqaeh4lSDgpXSKGob2a85ce9flu9Y90e3ka01bqfyxE8iOvaxDVnHOnRiHdCYx3zl0BtOT8hrl+BnQGMLNaQAywAxgN3GhmsWZWDcgEZp3H+4v8gpnRq3VVRvdtS3JCNL3ensX6nQc5qdubipxVXgFxLNd2F2AMgHNuH77/+z8jM/sQ+A6obWYbzexuYBhQ3T/09SPgdn9vYjHwMbAEGAv0cc6dOJ+/kMiZ1Clfks/7tuOmFhls3nOYpZv3sWn3oVCXJVJo5XWReoOZ/RbfqaFm+L68MbN4IPpsL3TO3XSGp24N1Oic6w/0z6MekV8lPiaSZ3s2YvaanazecYBLh0xj0PVN6FwnLdSliRQ6efUg7sY3N+EO4Abn3G5/eyvgHQ/rEvFU2cRYGqaXokKpeO4cPpvnvlqm+0yInOasPQj/PIX7ArRPAiZ5VZRIMMRFRzLi3hb8+fMlvDFlFXPW7uSVm5tSoVR8qEsTKRTOGhBmNvpszzvnrizYckSCKy46kmd7NqRV9TI8NXIhlw2ZzkvXN6ZTbZ1yEsnrGkRrYAPwITCTPNZfEimqejRJp37FUvR5fy53vDObPp1r8EjXWkRFakV8CV95ffrLA78HGgCDgYuAHc65Kc65KV4XJxJMNdMS+axPW27IyuDVSau45a2ZbN2rmxFJ+MprJvUJ59xY59zt+C5MZwOTdS8ICbUR97b2ZLmO+JhInr+2EQOva8yCjXu4dLAm1kn4yrP/7J+81hN4D+gDDAE+9bowkVC6pnklRvdtS5kSMdw2bBYvjVvOiQK8/7UW+5OiIK+lNt7FN9mtGfBn59wFzrm/OOcCLoMhUpxk+m9G1LNpJYZMzObWt2bq/tcSVvLqQdyKb9mLh4BvzWyv/2efme31vjyR0EqIiWLg9Y0ZcG0jftiwi0sHT+fb7B2hLkskKPKaB6EhHCLAdVkZNKqUzAPvf88tb8+kYql40pPjQl2WiKcUACL5VLt8EqP7tuOqJunk7D7Esi372L5Py4dL8aWAEDkHJWKjeOn6xlRLKcG+I8fpPniaTjlJsaWAEDlHZkZaUiwNKpakVHwUt7w9k0HjVxToKCeRwkABIXKeEmKiGN23HVc3TWfwhJW+UU6aWCfFiAJC5FfwnXJq8vMopyHTmL5Sp5ykeFBAiBSA67IyGN23HaUTYug1bCYDxy3nuJYPlyIur8X6RIolL5bpqOWfWPfHUYt5ZWI2s9bsZMhNTSlXUsNhpWhSD0KkACXERDHgusa8mGstpykrCn4tJy3VIcGggBDxwLX+tZzKJsZw+7BZDPh6mU45SZGjgBDxSGa5JEb1affT8uE3vzmTzXsOhboskXzzLCDMbJiZbTOzRQGee9TMnJml+B+bmQ0xs2wzW2BmzbyqSySYTi0fPuiGxiza5DvlNGn5tlCXJZIvXvYghgPdTm80swzgYmB9rubu+BYFzAR6A697WJdI0F3dtBKj+7ajXMk47nxnNut3HuSk08Q6Kdw8Cwjn3FRgZ4CnBgGPA7n/dfQA3nU+M4BkM6vgVW0ioXDqjnU3tchg857DLNm0l5Vb94W6LJEzCuo1CDPrAeQ45+af9lQ6vntfn7LR3yZSrMRFR/Jsz0ZkpiVy5PhJLntlOkOnrtIyHVIoBW0ehJkl4Lu/9cW/8n164zsNReXKlQugMpHgK1MihqS4KMomxvK3McsYv2QrL17XmCplS4S6NJGfBLMHUQOoBsw3s7VAJWCumZUHcoCMXPtW8rf9gnNuqHMuyzmXlZqa6nHJIt6JjoxgaK/mDLyuMcu27KPby9P414x1OF2bkEIiaAHhnFvonEtzzlV1zlXFdxqpmXNuCzAauM0/mqkVsMc5tzlYtYmEiplxTfNKfP1wB7KqluaZzxZx27BZbNrt7XBYTbST/PBymOuH+O5nXdvMNprZ3WfZfQywGsgG3gQe8KoukcKoYnI8797Vgr9e1YDv1+3ikpen8sn3G9WbkJDy7BqEc+6mPJ6vmmvbAX28qkWkoHmxlpOZcWurKrTPTOGxf8/n0X/P5+vFW+h/dUNSk2IL/PeJ5EUzqUUKmSplS/BR79Y8fWldJq/YziUvT+WrhTrjKsGngBAphCIjjHs6VOfL37YjPTme+9+fy0Mf/cCeg8dCXZqEEQWESCGWWS6JkQ+0od9FtfhywWYufnmKluqQoFFAiBRy0ZERPNglk8/6tCU5PoY735nN6u0HOH5Sq8OKtxQQIkVEg/RSjP5tW+7rWIPt+48wf8MePpy1XrOwxTMKCJEiJDYqkie716FBxZLER0fy1MiFXPHKdGau/jGodWgeRXhQQIgUQSVio6hbIYm/39yUPYeOccPQGfR5fy4bdh4MdWlSjCggRIooM+PyRhWZ8GhH+l1UiwnLttLlpSkMHLecg0ePh7o8KQYUECJFXFx0JA92yWTio53o3qA8r0zM5sIXp/DpDxs5qesT8isoIESKiYrJ8Qy+sSmf3N+atJKxPDJiPte88S3zNuwOdWlSRAVtuW8R+ZkXS3Wc0rxKGT57oC0jf8jh+bHLuOrVb+jZLJ0nutWhXMk4z36vFD/qQYgUQxERxrXNKzHpsU7c36kGX8zfTOcXJ/PqpGwOHzsR6vKkiFAPQqQYS4yN4oludbjpgsr0H7OEAV8v58NZ64mPjqR0QnTI6jo1RNbLnpT8eupBiISBymUT+EevLD74TUtKxESxctt+Fm/ay9hFW3QhW85IASESRtrUTOHLB9tRLaUEx0867nvvey5+eSoj527k+Akt3SH/SwEhEmaiIiNIS4qlcaVSDL6xCZFm9Pt4Pp0HTua9Get0jUJ+ooAQCVNmRo8m6Xz1UHveui2LsiVi+b/PFtHhhUm8OXU1B45osl24U0CIhLmICKNrvXJ8+kAbPvhNS2qmJdJ/zFLaPj+Rwf9dqXtQhDGNYhIRwNejaFMzhTY1U5i7fhevTcpm0H9XMHTqKm5tXYW721UjLUnzKMKJAkJEfqFZ5dK8dfsFLN28l9cnr+LNqasZ/s1abrggg94dqoe6PA2TDRLPTjGZ2TAz22Zmi3K1DTCzZWa2wMw+NbPkXM89ZWbZZrbczC7xqi4Ryb+6FUoy5KamTHi0E1c3TefDWevpNGAyq7bv14KAYcDLaxDDgW6ntY0HGjjnGgErgKcAzKwecCNQ3/+a18ws0sPaRIq0Efe2Dur/PVdLKcFz1zRiyu8606t1FXYeOMrCnL30ensmk5Zt01yKYsqzgHDOTQV2ntY2zjl36n87ZgCV/Ns9gI+cc0ecc2uAbKCFV7WJyPmpmBzPH6+oT5OMZCqVjmfl1v3cOXw2XV+awrvfrdXIp2ImlKOY7gK+8m+nAxtyPbfR3/YLZtbbzOaY2Zzt27d7XKKIBBIdGUF6cjzTnujMkJuakhQfzR9GLabVsxP425ilbNylGxcVByG5SG1mTwPHgffP9bXOuaHAUICsrCz1a0VCKDoygisbV+TKxhWZu34Xw6av4e3pa3hr2mq6NSjPnW2rkVWlNGYW6lLlPAQ9IMzsDuByoItz7tQXfA6QkWu3Sv42ESkimlUuTbObS7Np9yH+NWMdH8xcz5iFW2iYXoq72lXlsoYViYkqHFOvNAoqf4L6X8vMugGPA1c653L3QUcDN5pZrJlVAzKBWcGsTUQKRsXkeJ7oVofvnrqQ/lc34ODR4zwyYj5tn5/IKxNW8uP+I6EuUfLJsx6EmX0IdAJSzGwj8Ed8o5ZigfH+LucM59x9zrnFZvYxsATfqac+zjktCCNShCXERHFLyyrcdEFlpmXvYNj0NQwcv4JXJmVTKi6KtJJxOOd0+qkQ8ywgnHM3BWh++yz79wf6e1WPiIRGRITRsVYqHWulkr1tP8O/XcMHM9ezff9RugycQs9m6VzVNJ1KpRNCXaqcRjOpRSRoaqYl8terGrJ08152HjhGalIsL45bwYvjVtCyWhmuaVaJ7g3LkxQXupsZyc8UECJhKNQXZ6MifEuOj7i3NRt2HmTUvBxGzs3h8U8W8MyoRVxcvzw9m6XTvmYKUZGF48J2buFykVsBISIhlVEmgb4XZtKnc03mbdjNyLk5fL5gE5/P30RKYiw9mlSkZ7N06lUoqesVQaaAEJFCwcxoWrk0TSuX5pnL6zFp+TZGzt3Iu9+t5e3pa6hTPomrm/quV0hwKCBEpNCJiYrgkvrluaR+eXYdOMoXCzczcu5Gnv1qGc+PXUZibBQpibHsPXyMkrpe4RkFhIgUaqVLxNCrVRV6tarC6u37+eyHHP4xdTWrdxwg6y//pVPtVC5vXJGuddNIiCkaX2lF5RpG0TiaIiJA9dRE+l1cmxmrf+TAkRO0rF6WLxZsYtySrcRHR9KlbhpXNK5Ix1qpxEVrQehfSwEhIkWOmZEYF8UfrqjH05fVZfbanXyxYBNjFm7hiwWbSYqN4qL65biicUXa1UwhuhCOhCoKFBAiUqRFRhitqpelVfWy/OmK+ny76kc+n7+JsYu3MHJuDskJ0XRvUIErGlWgZfWyREYU/ZFQwTpFpYAQkXNWWM+dR0VG0KFWKh1qpfLXqxswbcUOPl+wiVHzcvhw1npSk2K5rGEF9h0+RmKsvv7yoiMkIsVSbFQkXeuVo2u9chw6eoJJy7fx+fxNfDhrPUeOnyQqwuj7wVzaZ6bQLjOV9OT4UJdc6CggRKTYi4+J5NKGFbi0YQX2HzlOj79PZ8+hY8xas5MvFmwGoHpqCdrXTKF9ZiqtapRVDwMFhIiEmVNzKFISY/modytWbN3PtJXbmZ69gxFzNvDP79YRFWE0rZxMu5qptK+VQqP0UoVyyQ+vKSBEJGyZGbXLJ1G7fBK/aV+dI8dP8P26XUxfuYPp2Tt4ecIKBv13BUlxUbSpUZZ2mal0yEwJddlBo4AQEfGLjYqkTY0U2tRI4XFg14GjfLNqB9NX7mDayh18vXirf78IkuKi+NeMdTTNSKZ2+aRiOZRWASEiQVdYR0GdrnSJGC5vVJHLG1XEOceaHQeYnr2DQeNXsPvgMZ75bBEAcdERNEpPpknlZJpm+P6sUKroX/RWQIiI5IOZUT01keqpiXy5YDPOOQZe34QfNuxm3vrd/LBhF8O/WcvQEycBKF8yjiYZyTStnEyTjGQaVipVZJYCOaVoVSsiUkiYGRllEsgok8CVjSsCcOT4CZZu3se89bt8wbFhN2MXbwF8E/pql0uiaeVktu87Qsn4wv/1W/grFBEpImKjImmS4esx3OFv+3H/EeZv3M0P630/o+dtYt+R4wB0e3kqneuk0aVOGk0rly50s7w9Cwhew5gUAAAILklEQVQzGwZcDmxzzjXwt5UBRgBVgbXA9c65Xea7C8hg4FLgIHCHc26uV7WJiARL2cRYLqxTjgvrlAPg5ElHj1d98zCSE6J5c+pqXp+8iuSEaDrWSuXCOml0rJVKckJMiCv3tgcxHPg78G6utieBCc6558zsSf/jJ4DuQKb/pyXwuv9PEZFiJSLCSIiJIiEmio96t2bv4WNMW7GDicu2MXn5NkbN20SEQfMqpf29i3LUKpcYkrvpeRYQzrmpZlb1tOYeQCf/9j+ByfgCogfwrnPOATPMLNnMKjjnNntVn4gUXUVlFFR+lIyL5rJGFbisUQVOnnTM37ibScu2MWHZNl4Yu5wXxi4nPTmeznV8vYs2NYI3DyPY1yDK5frS3wKU82+nAxty7bfR36aAEJGwERHx821X+11cm617D/8UFiPn5vDejPXERUcQFxVJWlKs5/WE7CK1c86ZmTvX15lZb6A3QOXKlQu8LhGRwqJcyThubFGZG1tU5sjxE8xcvZOJy7bx4az1HDp2wvPfH+yA2Hrq1JGZVQC2+dtzgIxc+1Xyt/2Cc24oMBQgKyvrnANGRKQoio2K/Gkp8yWb9hCML79gB8Ro4HbgOf+fo3K19zWzj/BdnN6j6w8i4pWifg3DzAjGJWsvh7l+iO+CdIqZbQT+iC8YPjazu4F1wPX+3cfgG+KajW+Y651e1SUiIvnj5Simm87wVJcA+zqgj1e1iIjIuSt+yw+KiEiB0FIbIiLnqKhfw8gv9SBERCQgBYSIiASkU0wiIkFWVE5RqQchIiIBKSBERCQgBYSIiASkgBARkYAUECIiEpBGMYmIFDHBGgWlHoSIiASkgBARkYAUECIiEpACQkREAlJAiIhIQAoIEREJSAEhIiIBKSBERCQgBYSIiARkzrlQ13DezGw7sO5XvEUKsKOAyinKdBx8dBx8dBx+VlyPRRXnXGpeOxXpgPi1zGyOcy4r1HWEmo6Dj46Dj47Dz8L9WOgUk4iIBKSAEBGRgMI9IIaGuoBCQsfBR8fBR8fhZ2F9LML6GoSIiJxZuPcgRETkDMIyIMzsUTNzZpZyhudPmNk8/8/oYNcXLPk4Dreb2Ur/z+3Brs9rZvYXM1vg/+88zswqnmG/Yv15OIfjUNw/DwPMbJn/WHxqZsln2G+tmS30H685wa4zmMLuFJOZZQBvAXWA5s65X4xxNrP9zrnEoBcXRHkdBzMrA8wBsgAHfO/fb1ewa/WKmZV0zu31bz8I1HPO3Rdgv2L9ecjPcQiTz8PFwETn3HEzex7AOfdEgP3WAlmBvjuKm3DsQQwCHsf3IQ9neR2HS4Dxzrmd/i+B8UC3YBUXDKe+FP1KEKafiXweh3D4PIxzzh33P5wBVAplPYVBWAWEmfUAcpxz8/PYNc7M5pjZDDO7Khi1BVM+j0M6sCHX443+tmLFzPqb2QbgFuAPZ9itWH8eIF/HISw+D7ncBXx1huccMM7Mvjez3kGsKeiiQl1AQTOz/wLlAzz1NPB74OJ8vE0V51yOmVUHJprZQufcqoKs02sFdByKvLMdB+fcKOfc08DTZvYU0Bf4Y4B9i/Xn4RyOQ5GX13Hw7/M0cBx4/wxv087/eUgDxpvZMufcVG8qDq1iFxDOua6B2s2sIVANmG9m4Os+zjWzFs65Lae9R47/z9VmNhloChSpL4QCOA45QKdcjysBkz0p1kNnOg4BvA+MIcAXY3H+PARwpuMQFp8HM7sDuBzo4s5wgTbX52GbmX0KtACKZUCEzSkm59xC51yac66qc64qvi5ys9PDwcxKm1msfzsFaAssCXrBHsnvcQC+Bi72H4/S+HocXwe5XE+ZWWauhz2AZQH2KdafB8jfcSA8Pg/d8F2Xu9I5d/AM+5Qws6RT2/iOw6LgVRlcYRMQZ2NmWWb2lv9hXWCOmc0HJgHPOeeK1RfCmeQ+Ds65ncBfgNn+n//nbytOnjOzRWa2AN8/9IcgLD8PeR6HMPk8/B1IwnfaaJ6ZvQFgZhXNbIx/n3LAdP/nYRbwpXNubGjK9V7YDXMVEZH8UQ9CREQCUkCIiEhACggREQlIASEiIgEpIEREJCAFhIiIBKSAEBGRgBQQIvlgZvt/xWvjzWyKmUX6Hzc0s3Vmdn+ufWLMbKqZFbvlb6ToUkCIeO8uYKRz7gT4ljsBbgRuO7WDc+4oMAG4ISQVigSggBA5B2bWz78sxSIzezhX+zNmttzMppvZh2b2WK6X3QKMOu2ttgH1T2v7zL+vSKGg7qxIPplZc+BOoCVgwEwzm4Lv39E1QGMgGpiL745rmFkMUN05t/a0t3sOiDWzKs65df62RcAFXv89RPJLASGSf+2AT51zBwDMbCTQHl9PfJRz7jBw2Mw+z/WaFGB37jcxs+747tz2Jb5exDoA59wJMztqZknOuX2e/21E8qBTTCLeOgTEnXpgZnHA88ADwEKgwWn7xwKHg1adyFkoIETybxpwlZkl+O8FcLW/7RvgCjOLM7NEfDecAcB//+ZIfzAA/B/wrv+U0/8EhJmVBXY4544F5W8jkgedYhLJJ+fcXDMbju8+AABvOed+ADCz0cACYCu+L/49uV46Dmjnv+fzRfhuOoR/v9/n2q8zvtNOIoWC7gchUgDMLNE5t9/MEvDdfrK3c26u/7lmwCPOuV55vMdI4Enn3ArvKxbJm3oQIgVjqJnVw3e94Z+nwgF+6nlMMrPIU3MhTucf7fSZwkEKE/UgREQkIF2kFhGRgBQQIiISkAJCREQCUkCIiEhACggREQlIASEiIgEpIEREJCAFhIiIBPT/AclGLmYlr42AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.errorbar(np.log(lagrange_val), cv_MSE, yerr=cv_MSE_sd, label='MSE')\n",
    "plt.gca().set_ylabel('MSE')\n",
    "plt.gca().set_xlabel(r'$\\log(\\lambda)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "all,-slideshow"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Compare to sglasso.ipynb

Comparison to regular sparse group LASSO.ipynb

Cox stacked block - Ruilin penalty.ipynb

Cox stacked block.ipynb

Cox stacked.ipynb

Generating an instance.ipynb