{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Smooth Bilevel Programming for Sparse Regularization\n", "\n", "This tour is a tutorial reviewing the method developped in the paper:\n", "\n", "> _[Smooth Bilevel Programming for Sparse Regularization](https://arxiv.org/abs/2106.01429)_\n", "Clarice Poon, Gabriel PeyrĂ©, 2021\n", "\n", "We present a surprisingly simple reformulation of the Lasso as a bilevel program, which can be solved using efficient method such as L-BFGS. This gives an algorithm which is simple to implement and is in the same ballpark as state of the art approach when the regularization parameter $\\lambda$ is small (it can in particular cope in a seamless way with the constraint formulation when $\\lambda=0$). For large $\\lambda$, and when the solution is very sparse, the method is not as good as for instance coordinate methods with safe screening/support pruning such as the [Celer](https://github.com/mathurinm/celer) solver which we warmly recommend." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time\n", "import warnings\n", "import progressbar" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Lasso problem\n", "\n", "We aim at solving the [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) problem\n", "$$\n", " \\min_x \\frac{1}{2\\lambda}\\| A x -y \\|^2 + \\|x\\|_1. \n", "$$\n", "We generate here a synthetic example using a random matrix $A \\in \\mathbb{R}^{n \\times p}$ and $y \\in \\mathbb{R}^p$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "n = 10;\n", "p = n*2;\n", "A = np.random.randn(n,p)\n", "y = np.random.randn(n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simplest (but arguably not the fastest) algorithm to solve this problem is the iterative soft thresholding\n", "$$\n", " x_{k+1} = \\text{Thresh}_{\\tau \\lambda}( y - \\tau A^\\top (Ax_k-y) )\n", "$$\n", "where the step size should satisfies $\\tau < 2/\\|A\\|$, and Thresh is the soft thresholding\n", "$$\n", " \\text{Thresh}_\\sigma( x ) = \\text{sign}(x) \\max(|x|-\\sigma,0)\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhhklEQVR4nO3deXxU5dn/8c9lTAg7CGENCCiiiBJCRNGK4oK7KI9FEapP6yMStS6tVqxVf1qX2qe1VkWorba1REBFBAvuS90VEsK+GlkiSEKAsIaQ5Pr9weiDmkjITOZkZr7v1ysvZuacnPs6Qb+5uefMdczdERGR+HdQ0AWIiEh0KPBFRBKEAl9EJEEo8EVEEoQCX0QkQRwcdAE/pG3btt6tW7egyxARiRm5ubkb3T2tum0NOvC7devGnDlzgi5DRCRmmNnqmrZpSUdEJEEo8EVEEoQCX0QkQTToNfzq7Nmzh8LCQsrKyoIuJe6lpqaSnp5OcnJy0KWISATEXOAXFhbSvHlzunXrhpkFXU7ccndKSkooLCyke/fuQZcjIhEQkSUdM3vazIrMbGEN2081s1Izyw993VXXscrKymjTpo3Cvp6ZGW3atNG/pETiSKRm+P8AHgee+YF93nf38yMxmMI+OvRzFokvEZnhu/t7wKZIHEtEJJF99sUm/vZ+AfXRuj6aV+kMNLN5ZvaKmR1d005mNtrM5pjZnOLi4iiWJyISrKJtZVz3bB45n65h157KiB8/WoGfBxzq7n2Bx4CXatrR3Z909yx3z0pLq/bTwSIicaeisoqfPzuXbWV7GD8qkyYpkb+mJiqB7+5b3X176PEsINnM2kZj7EgrKSkhIyODjIwMOnToQOfOnb95vnz5cvr06RPW8VetWhX2Mb62a9cuTjnlFCora54plJeXM2jQICoqKiIypojUzR9eX86nX2zigYuP4cgOLepljKgEvpl1sNA7gGY2IDRuSTTGjrQ2bdqQn59Pfn4+Y8aM4eabb/7meUpKStDlfcvTTz/NsGHDSEpKqnGflJQUTj/9dKZMmRLFykRkX28s3sCE/3zO5cd3ZVhmer2NE6nLMicBHwO9zKzQzK4yszFmNia0yyXAQjObBzwKXOZxejPdyspKrr76ao4++miGDBnCrl27AJg4cSIDBgwgIyODa665hsrKSnbs2MF5551H37596dOnzzehW9MxDlROTg5Dhw795vngwYN54403APjNb37DDTfcAMBFF11ETk5OOKctInW0umQHv3gun2M6t+Su83vX61gRWSRy9xH72f44ey/bjKh7Xl7E4nVbI3rM3p1acPcFNb6nvF8rVqxg0qRJ/PWvf2X48OFMnTqV/v37M2XKFD788EOSk5O59tprycnJoWnTpnTq1ImZM2cCUFpayubNm6s9xqhRow6ojvLycgoKCti3vfQ999zDXXfdRVFREXPnzmXGjBkA9OnTh9mzZ9f5nEWkbsr2VJI9MY+DzHhiZCapyTX/azwS1Esnwrp3705GRgYA/fv3Z9WqVbz11lvk5uZy3HHHkZGRwVtvvUVBQQHHHHMMb775Jrfddhvvv/8+LVu2rPEYAAUFBVx11VVccskl3xv33Xff5eSTT2bMmDG8++67bNy4kVatWn1rn0GDBuHuPPzww0yePPmbpZ6kpCRSUlLYtm1bvfxMRKR6d09fxOL1W/nTpX3pckiTeh8v5lor7CucmXh9adSo0TePk5KS2LVrF+7OlVdeyYMPPvi9/XNzc5k1axa33347Q4YM4Yorrqj2GAA9evTgqaeeqjbwzYxmzZpRVlZGeno6jRs3/t6nZBcsWMD69etp27YtzZs3/9a23bt3k5qaGta5i0jtPTd7LVPmrOX6wYdz2pHtozKmZvhRcPrpp/PCCy9QVFQEwKZNm1i9ejXr1q2jSZMmjBo1iltuuYW8vLw6j3HyySfzyiuv8NBDD3H33XfTunVrKisrvwn99evXM3LkSKZPn07Tpk157bXXvvnekpIS0tLS1CRNJEoWrSvlzukLOenwNtx85hFRGzemZ/ixonfv3tx3330MGTKEqqoqkpOTGTduHKWlpdx6660cdNBBJCcnM378+DqPcdBBe393t27dmt27dwMwZMgQPvjgA0488USGDRvGH//4R4466ijuvPNObrvtNs466ywA3nnnHc4999zwT1RE9qt01x6yJ+bRukkKf76sH0kHRbGFibs32K/+/fv7dy1evPh7ryWKjRs3+jXXXOM9evTwBx54wD/66CMfN26cu7tPnTrVR48e7cOHD/d33nnH3d3z8vJ81KhR+z3uxRdf7EuXLq12WyL/vEUiraqqyv/nn7P9sNtn+pxVJfUyBjDHa8hUzfBjSJs2bZgwYcK3Xhs4cCAAw4YNY9iwYd/a1q9fPwYPHkxlZWWN1+KXl5dz0UUX0atXr/opWkS+8Zf3Cnhj8QbuOr83/Q89JOrjK/Dj3M9+9rMf3J6SksIVV1wRpWpEEtcnBSX8/tWlnHdsR356UrdAatCbtiIi9axoaxnXPzuXbm2b8tB/HRtY63HN8EVE6lFFZRXXT5rLjt0VPHv18TRrFFzsxuQM3+OzK0ODo5+zSPj+97VlfPbFJh4cdgxHtG++/2+oRzEX+KmpqZSUlCiM6pmH7mmrD2OJ1N1ri77iL+8VMOqErlzUr3PQ5cTekk56ejqFhYXo5ij1LzU1lfT0+uvcJxLPVm3cwS3PzaNvekvurOemaLUVc4GfnJxM9+7dgy5DRKRGZXsqyc7JIynJGDcyk0YH129TtNqKucAXEWnI3J3fvLSQpV9t5en/Po701vXfFK22Ym4NX0SkIZsyey0v5Bby88GHM7hXu6DL+RYFvohIhCz8spS7Zizi5J5tufGM6DVFqy0FvohIBJTu3EN2Ti5tmqbwyKUZ0W2KVktawxcRCVNVlfPL5/P5qrSMKdcMpE2zRvv/pgBohi8iEqYJ733Om0uKuOPco8js2jrocmqkwBcRCcNHn2/kD68t44K+nbjyxG5Bl/ODFPgiInW0YWsZN0yaS/e2TfndsGMCa4pWW1rDFxGpgz2VVVz/bB47yyuZdPUJNA2wKVptNfwKRUQaoN+/upTZqzbz58sy6BlwU7Ta0pKOiMgBenXhev76/hdcMfBQhmYE3xSttiIS+Gb2tJkVmdnCGrabmT1qZivNbL6ZZUZiXBGRaCso3s4tz8+nb5dW3HHeUUGXc0AiNcP/B3D2D2w/B+gZ+hoNjI/QuCIiUbOrvJJrc/JITjKeaEBN0WorIoHv7u8Bm35gl6HAM6Gbqn8CtDKzjpEYW0QkGtydO15awLIN23jksn50btU46JIOWLTW8DsDa/d5Xhh67XvMbLSZzTGzOep5LyINxaTP1vJi3pfccFpPTjkiLehy6iRagV/dxanV3rLK3Z909yx3z0pLi80fqojElwWFpfy/UFO0G07vGXQ5dRatwC8EuuzzPB1YF6WxRUTqbMvOcrJzcmnbLIU/X9avQTZFq61oBf4M4IrQ1TonAKXuvj5KY4uI1ElVlfOL5+axYWsZT4zqzyFNU4IuKSwR+eCVmU0CTgXamlkhcDeQDODuE4BZwLnASmAn8NNIjCsiUp/G/+dz3l5axL1DjyajS6ugywlbRALf3UfsZ7sD10ViLBGRaPhw5Ub++PoyLuzbiZ+ccGjQ5USEPmkrIvIdX5XubYrWI60ZD8ZAU7TaUi8dEZF9fN0UbdeeSqaMyoyJpmi1FT9nIiISAb97ZSlzVm/msRH9OLxdbDRFqy0t6YiIhMxasJ6nPviC/z6xGxf07RR0ORGnwBcRAT4v3s6tz8+jX9dW/Prc2GqKVlsKfBFJeDvLK8iemEuj5CTGXZ5JysHxGY1awxeRhObu3DFtISuKtvPMzwbQKQabotVWfP4aExGppZxP1zBt7pfcdPoRnNwzvvt3KfBFJGHNL9zCvS8v5pQj0vj5aYcHXU69U+CLSELavKOc7Il5pDVvxCOXZnBQDDdFqy2t4YtIwqmqcm5+Lp/ibbt5fsxAWsd4U7Ta0gxfRBLOuHdW8u6yYu68oDd946ApWm0p8EUkoXywYiMPv7mcizI6Mer4rkGXE1UKfBFJGOtLd3HD5Ln0bNeMB+KoKVptKfBFJCGUV1RxXU4eu/dUMn5Uf5qkJN5bmIl3xiKSkB58ZQl5a7Yw7vJMDktrFnQ5gdAMX0Ti3svz1vH3D1fx05O6cd6xHYMuJzAKfBGJayuLtjN26nwyu7bi9nPisylabSnwRSRu7di9T1O0kfHbFK22tIYvInHJ3fn1tAWsLN7Ov352PB1bxm9TtNpK7F93IhK3Jn6ymun56/jFGUfwo55tgy6nQVDgi0jcyV+7hXv/vZjBvdK4bnD8N0WrLQW+iMSVzTvKuS4nj3bNU/lTgjRFqy2t4YtI3Kiqcm6asrcp2gvZA2nVJDGaotVWRGb4Zna2mS0zs5VmNraa7aeaWamZ5Ye+7orEuCIi+3rs7ZX8Z3kxd1/Ym2PTWwVdToMT9gzfzJKAccCZQCEw28xmuPvi7+z6vrufH+54IiLVeW95MY+8tZxh/Tpz+YDEaopWW5GY4Q8AVrp7gbuXA5OBoRE4rohIrazbsosbJ8/liHbNuf/ixGuKVluRCPzOwNp9nheGXvuugWY2z8xeMbOjazqYmY02szlmNqe4uDgC5YlIPCuvqOLanDz2VDrjR2XSOCUp6JIarEgEfnW/Sv07z/OAQ929L/AY8FJNB3P3J909y92z0tLi+4bCIhK++2cuJn/tFn5/ybH0SNCmaLUVicAvBLrs8zwdWLfvDu6+1d23hx7PApLNTJ+EEJGwzJi3jn9+vJqrftSdc49J3KZotRWJwJ8N9DSz7maWAlwGzNh3BzPrYKFFNTMbEBq3JAJji0iCWrFhG2Onzifr0NaMPefIoMuJCWFfpePuFWZ2PfAakAQ87e6LzGxMaPsE4BIg28wqgF3AZe7+3WUfEZFa2bG7guycPJqkJPH45ZkkJ+kzpLURkQ9ehZZpZn3ntQn7PH4ceDwSY4lIYnN3xr64gILi7Uy86ng6tEwNuqSYoV+LIhJTnvl4NS/PW8cvh/TixMP1VuCBUOCLSMzIW7OZ+2Yu5vQj25F9ymFBlxNzFPgiEhM27Sjn+pw82rdI5eHhaopWF2qeJiINXmWVc+PkuWzcUc6L2SfSskly0CXFJM3wRaTBe/StFby/YiP3XHg0fTq3DLqcmKXAF5EG7d1lRTz69gr+KzOdy47rsv9vkBop8EWkwSrcvJObpuTTq31z7ruoj5qihUmBLyIN0u6KSq7LyaOy0hk/qr+aokWA3rQVkQbpvn8vYV5hKRNGZdK9bdOgy4kLmuGLSIMzPf9L/vXJaq4+uTtn91FTtEhR4ItIg7J8wzbGTl3Acd1a86uz1RQtkhT4ItJgbN9dwZiJuTRtdLCaotUD/TRFpEFwd26bOp9VG3fw2Ih+tG+hpmiRpsAXkQbhHx+tYub89dxyVi8GHtYm6HLikgJfRAKXu3oz989cwhlHtWfMIDVFqy8KfBEJVMn23Vz/bB6dWjXmj8P7qilaPdJ1+CISmL1N0fIp+bopWmM1RatPmuGLSGD+/OZyPli5kd8OVVO0aFDgi0gg3llaxKNvr+TH/dO59LiuQZeTEBT4IhJ1azftbYp2VMcW/PaiPkGXkzAU+CISVbsrKrnu2TyqqpzxIzNJTVZTtGjRm7YiElX3vryY+YWl/OUn/emmpmhRpRm+iETNtLmF5Hy6hmsG9eCsozsEXU7CUeCLSFQs+2obt7+4gAHdD+HWs3oFXU5Cikjgm9nZZrbMzFaa2dhqtpuZPRraPt/MMiMxrojEhm1le8iemEuzRsk8PqIfB6spWiDC/qmbWRIwDjgH6A2MMLPe39ntHKBn6Gs0MD7ccUUkNrg7Y6cuYPWmnTx+eT/aqSlaYCLxpu0AYKW7FwCY2WRgKLB4n32GAs+4uwOfmFkrM+vo7usjML6EadJna3hveXHQZUic2rJzDx8XlDD2nCM5oYeaogUpEoHfGVi7z/NC4Pha7NMZ+F7gm9lo9v4rgK5d9WGM+vbqwvXc/uICOrdqTNNGujxO6seYUw7jmkE9gi4j4UUi8KvrdOR12Gfvi+5PAk8CZGVlVbuPRMYXG3dw6/Pz6dulFc9dcwKNDlbgi8SzSLxzUgh02ed5OrCuDvtIFO0qryR7Yi4HJxlPjMxU2IskgEgE/mygp5l1N7MU4DJgxnf2mQFcEbpa5wSgVOv3wXF3fvPSQpZt2MYjl/Wjc6vGQZckIlEQ9pKOu1eY2fXAa0AS8LS7LzKzMaHtE4BZwLnASmAn8NNwx5W6mzx7LVPzCrnpjJ6cckRa0OWISJREpLWCu89ib6jv+9qEfR47cF0kxpLwLCgs5e7pixh0RBo3nNYz6HJEJIr06YcEsmVnOdk5ubRtlsIjl2bozkIiCUbN0xJEVZXzi+fmsWFrGc+POZFDmqYEXZKIRJlm+AniiXdX8vbSIu46vzcZXVoFXY6IBECBnwA+XLmRh99YztCMTow64dCgyxGRgCjw49xXpWXcMGkuh6U148Fhx2CmdXuRRKXAj2PlFVVcm5NL2Z5Kxo/qT5MUvWUjksiUAHHswVeWkLdmC49f3o/D2zULuhwRCZhm+HHq3/PX8fcPV/HTk7px/rGdgi5HRBoABX4cWlm0ndtemE9m11bcfs5RQZcjIg2EAj/O7NhdQfbEXFKTkxg3MpOUg/VXLCJ7aQ0/jrg7v562gM+Lt/Ovq46nY0s1RROR/6PpXxyZ+Mlqpuev4xdnHsFJh7cNuhwRaWAU+HEif+0W7v33Yk47sh3Xnnp40OWISAOkwI8Dm3aUc+3EXNq3SOXh4X3VFE1EqqU1/BhXWeXcNCWfjdvLmZp9Iq2aqCmaiFRPgR/jHnt7Be8tL+aBi4/hmPSWQZcjIg2YlnRi2LvLivjzWysYltmZEQO67P8bRCShKfBj1JdbdnHTlHx6tW/O/RepKZqI7J8CPwbtrqjk2pw8Kiud8aP60zglKeiSRCQGaA0/Bt0/cwnz1m5hwqhMurdtGnQ5IhIjNMOPMdPzv+SZj1czelAPzu7TMehyRCSGKPBjyPIN2xg7dQEDuh3Cr87qFXQ5IhJjFPgxYvvuCsZMzKVpo4N5/PJ+HJykvzoROTBKjRjg7tz2wnxWbdzBYyP60a5FatAliUgMCutNWzM7BJgCdANWAcPdfXM1+60CtgGVQIW7Z4UzbqL5+4ermLlgPWPPOZKBh7UJuhwRiVHhzvDHAm+5e0/grdDzmgx29wyF/YHJXb2JB2Yt4cze7blmUI+gyxGRGBZu4A8F/hl6/E/gojCPJ/vYuH031+XMpXPrxvzhx3314SoRCUu4gd/e3dcDhP5sV8N+DrxuZrlmNvqHDmhmo81sjpnNKS4uDrO82FVZ5dw4eS6bd5bzxMhMWjZODrokEYlx+13DN7M3gQ7VbLrjAMY5yd3XmVk74A0zW+ru71W3o7s/CTwJkJWV5QcwRlz50xvL+XBlCb+/5FiO7qSmaCISvv0GvrufUdM2M9tgZh3dfb2ZdQSKajjGutCfRWY2DRgAVBv4Am8v3cDj76zk0qwuDM9SUzQRiYxwl3RmAFeGHl8JTP/uDmbW1Myaf/0YGAIsDHPcuLV2005umpxP744tuGfo0UGXIyJxJNzA/x1wppmtAM4MPcfMOpnZrNA+7YEPzGwe8Bkw091fDXPcuFS2p5LsnFwAJozqT2qymqKJSOSEdR2+u5cAp1fz+jrg3NDjAqBvOOMkinteXszCL7fytyuy6NqmSdDliEic0SdtG4gXcguZ9Nkask89jDN6tw+6HBGJQwr8BmDJ+q3cMW0BA3u04ZdnHhF0OSISpxT4Adtatofsibm0bJzMoyPUFE1E6o9ugBIgd+fW5+exdvMuJo8+gbTmjYIuSUTimKaTAfrr+wW8tmgDt59zJMd1OyTockQkzinwA/JpQQkPvbqMc/p04KofdQ+6HBFJAAr8ABRtK+P6SXM59JAm/P6SY9UUTUSiQmv4UVZRWcXPn53LtrI9/OuqATRPVVM0EYkOBX6U/e/ry/j0i008PLwvR3ZoEXQ5IpJAtKQTRa8v+oq//KeAkcd3ZVhmetDliEiCUeBHyaqNO/jl8/M4Nr0ld13QO+hyRCQBKfCjYG9TtDwOMmPc5Zk0OlhN0UQk+rSGHwV3vrSQJeu38vf/Po4uh6gpmogEQzP8ejZl9hqezy3khtMOZ/CRNd0BUkSk/inw69HCL0u5c/oiTu7ZlhvPUFM0EQmWAr+elO7cQ3ZOLm2apvDIpRkkHaQPV4lIsLSGXw+qqpxfPp/PV6VlTLlmIG2aqSmaiARPM/x6MOG9z3lzSRF3nHsUmV1bB12OiAigwI+4jz7fyB9eW8YFfTtx5Yndgi5HROQbCvwI+qq0jBsmzaV726b8btgxaoomIg2K1vAjZE9lFdc/m8fO8komXX0CTRvpRysiDYtSKUIeemUpc1Zv5tER/ejZvnnQ5YiIfI+WdCJg1oL1/O2DL7hy4KFc2LdT0OWIiFRLgR+mguLt/OqF+WR0acUd56kpmog0XGEFvpn92MwWmVmVmWX9wH5nm9kyM1tpZmPDGbMh2VleQfbEPJKTjCdGZpJysH5/ikjDFW5CLQSGAe/VtIOZJQHjgHOA3sAIM4v5qbC7c8e0hSwv2sajI/rRqVXjoEsSEflBYb1p6+5LgP1dfjgAWOnuBaF9JwNDgcXhjB20nE/XMG3ul9x8xhGc3DMt6HJERPYrGmsQnYG1+zwvDL1WLTMbbWZzzGxOcXFxvRdXF/MLt3Dvy4s5tVcaPz/t8KDLERGplf3O8M3sTaBDNZvucPfptRijuum/17Szuz8JPAmQlZVV435B2byjnOyJeaQ1b8SfhmdwkJqiiUiM2G/gu/sZYY5RCHTZ53k6sC7MYwaiqsq5+bl8irft5vkxA2ndNCXokkREai0aSzqzgZ5m1t3MUoDLgBlRGDfiHn9nJe8uK+bOC3rTt0uroMsRETkg4V6WebGZFQIDgZlm9lro9U5mNgvA3SuA64HXgCXAc+6+KLyyo+/9FcX86c3lXNyvM6OO7xp0OSIiB8zcG9wy+TeysrJ8zpw5QZfBui27OP+xD2jbLIWXrjuJJinqSCEiDZOZ5bp7tZ+L0ieF9qO8ooprc/Ior6hi/Kj+CnsRiVlKr/14YNYS8tdu4YmRmRyW1izockRE6kwz/B8wY946/vHRKq76UXfOPaZj0OWIiIRFgV+DFRu2MXbqfLIObc3Yc44MuhwRkbAp8KuxY3cF2Tl5NElJ4vHLM0lO0o9JRGKf1vC/w90Z++ICCoq3M/Gq4+nQMjXokkREIkJT1+945uPVvDxvHb8c0osTD28bdDkiIhGjwN9H3prN3DdzMacf2Y7sUw4LuhwRkYhS4IeUbN/NdTl5dGiZysNqiiYicUhr+EBllXPTlHxKdpTzYvaJtGySHHRJIiIRpxk+8Oc3l/P+io3ce+HR9OncMuhyRETqRcIH/jvLinj07ZVc0j+dS4/rsv9vEBGJUQkd+IWbd3LzlHyO6tiC3w7ts79bNYqIxLSEDfzdFZVcm5NHZaUzfmQmjVOSgi5JRKReJeybtve+vJj5haX85Sf96da2adDliIjUu4Sc4U+bW0jOp2u4ZlAPzjq6utv1iojEn4QL/KVfbeX2FxcwoPsh3HpWr6DLERGJmoQK/G1le8iemEfz1GQev7wfB6spmogkkIRZw3d3fvXCfNZs2smz/3M87ZqrKZqIJJaEmeI+9cEXvLLwK351Vi+O79Em6HJERKIuIQJ/9qpN/O6VpQzp3Z7Rg3oEXY6ISCDiPvCLt+1tipbeujF/GN5XH64SkYQV12v4FZVV3DBpLqW79vCPnw6gRaqaoolI4orrwH/4jeV8XFDCH37cl96dWgRdjohIoMJa0jGzH5vZIjOrMrOsH9hvlZktMLN8M5sTzpi19cbiDTzx7ueMGNCFS/qnR2NIEZEGLdwZ/kJgGPCXWuw72N03hjlerawp2ckvnsunT+cW3H3B0dEYUkSkwQsr8N19CdCg3ggt21NJdk4uBowf2Z/UZDVFExGB6F2l48DrZpZrZqN/aEczG21mc8xsTnFx8YEP5NCrfXP+dGkGXQ5pUtd6RUTizn5n+Gb2JlBdh7E73H16Lcc5yd3XmVk74A0zW+ru71W3o7s/CTwJkJWV5bU8/jcapyTx8KUZB/ptIiJxb7+B7+5nhDuIu68L/VlkZtOAAUC1gS8iIvWj3pd0zKypmTX/+jEwhL1v9oqISBSFe1nmxWZWCAwEZprZa6HXO5nZrNBu7YEPzGwe8Bkw091fDWdcERE5cOFepTMNmFbN6+uAc0OPC4C+4YwjIiLhi/teOiIispcCX0QkQSjwRUQShAJfRCRBmPsBf7YpasysGFhdx29vC0Sld08UxMu5xMt5gM6lIYqX84DwzuVQd0+rbkODDvxwmNkcd6+xg2csiZdziZfzAJ1LQxQv5wH1dy5a0hERSRAKfBGRBBHPgf9k0AVEULycS7ycB+hcGqJ4OQ+op3OJ2zV8ERH5tnie4YuIyD4U+CIiCSKuA9/Mfmtm80M3T3/dzDoFXVNdmNn/mtnS0LlMM7NWQddUV7W98X1DZWZnm9kyM1tpZmODriccZva0mRWZWUy3KzezLmb2jpktCf23dWPQNdWVmaWa2WdmNi90LvdE9PjxvIZvZi3cfWvo8Q1Ab3cfE3BZB8zMhgBvu3uFmT0E4O63BVxWnZjZUUAVe298f4u7zwm4pFozsyRgOXAmUAjMBka4++JAC6sjMxsEbAeecfc+QddTV2bWEejo7nmhe2/kAhfF4t+L7b1BeFN3325mycAHwI3u/kkkjh/XM/yvwz6kKXvvrRtz3P11d68IPf0ESA+ynnC4+xJ3XxZ0HXU0AFjp7gXuXg5MBoYGXFOdhW4zuinoOsLl7uvdPS/0eBuwBOgcbFV143ttDz1NDn1FLLfiOvABzOx+M1sLjATuCrqeCPgZ8ErQRSSozsDafZ4XEqPBEq/MrBvQD/g04FLqzMySzCwfKALecPeInUvMB76ZvWlmC6v5Ggrg7ne4excgB7g+2Gprtr/zCO1zB1DB3nNpsGpzLjHKqnktJv/VGI/MrBkwFbjpO/+6jynuXunuGez9l/wAM4vYcltYd7xqCA7gJuvPAjOBu+uxnDrb33mY2ZXA+cDp3sDfeInEje8bqEKgyz7P04F1AdUi+witd08Fctz9xaDriQR332Jm7wJnE6H7gMf8DP+HmFnPfZ5eCCwNqpZwmNnZwG3Ahe6+M+h6EthsoKeZdTezFOAyYEbANSW80BudTwFL3P3hoOsJh5mlfX0Vnpk1Bs4ggrkV71fpTAV6sfeqkNXAGHf/MtiqDpyZrQQaASWhlz6JxauNYO+N74HHgDRgC5Dv7mcFWtQBMLNzgUeAJOBpd78/2IrqzswmAaeytxXvBuBud38q0KLqwMx+BLwPLGDv/+sAv3b3WcFVVTdmdizwT/b+93UQ8Jy73xux48dz4IuIyP+J6yUdERH5Pwp8EZEEocAXEUkQCnwRkQShwBcRSRAKfBGRBKHAFxFJEP8fs4LjUOhKSyEAAAAASUVORK5CYII=\n", "text/plain": [ "