{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HW3 - Pierre Fernandez"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 1\n",
"\n",
"Derive the dual problem of LASSO $$\\textrm{minimize} \\quad \\frac{1}{2} \\| Xw-y \\|_2^2 + \\lambda \\| w\\|_1 \\qquad \\textrm{(LASSO)}$$\n",
"\n",
"### Answer\n",
"\n",
"First note that \n",
"$$ \n",
"\\begin{aligned} \n",
"\\textrm{min}_w \\quad \\frac{1}{2} \\| Xw-y \\|_2^2 + \\lambda \\| w\\|_1 \\quad \\Leftrightarrow \\quad \\textrm{min}_{w,z} \\quad &\\frac{1}{2} \\| z \\|_2^2 + \\lambda \\| w\\|_1 \\\\ &Xw-y = z \n",
"\\end{aligned} \n",
"$$\n",
"\n",
"The associated Lagragian is \n",
"$$\n",
"\\begin{aligned} \n",
"L(w,z,\\nu) &= \\frac{1}{2} \\| z \\|^2_2 + \\lambda \\| w\\|_1 + \\nu^T(Xw-y-z)\\\\\n",
" &= \\frac{1}{2} \\| z \\|^2_2 - \\nu^Tz + \\lambda \\| w\\|_1 + \\nu^TXw - \\nu^Ty\n",
"\\end{aligned}\n",
"$$\n",
"Therefore $$ g(\\nu) = \\inf_{w,z} L(w,z,\\nu) = \n",
"-\\sup_z \\left(\\nu^Tz - \\frac{1}{2}\\| z \\|^2_2 \\right)\n",
"-\\lambda \\sup_w \\left( -\\frac{1}{\\lambda}\\nu^TXw - \\| w \\|_1 \\right) - \\nu^Ty$$\n",
"\n",
"Now, recall (from example 3.26 from the book and the fact that the dual norm of the $L_1$-norm is the $L_\\infty$-norm) that the conjugate function of $f=\\|.\\|$ is: $$f^\\star(y)=\n",
"\\begin{cases}\n",
"0 \\qquad \\textrm{if } \\|y\\|_\\infty \\leq1\n",
"\\newline \\infty \\qquad \\textrm{othewise}\n",
"\\end{cases} $$\n",
"and that (from example 3.27 from the book and the fact that $L_2$-norm is self-dual) that the conjugate of $f=\\frac{1}{2}\\|.\\|_2^2$ is $f^\\star=\\frac{1}{2}\\|.\\|_2^2$.\n",
"\n",
"- if $ \\| -\\frac{1}{\\lambda}\\nu^TX \\|_\\infty \\leq 1 \\Leftrightarrow \\| \\nu^TX \\|_\\infty \\leq \\lambda $, $g(\\nu) = -\\infty$\n",
"- otherwise, $g(\\nu) = -\\frac{1}{2} \\| \\nu \\|_2^2 - \\nu^Ty $\n",
"\n",
"Therefore the dual problem is:\n",
"$$ \n",
"\\begin{aligned} \n",
"\\textrm{min}_{\\nu} \\quad &\\frac{1}{2} \\| \\nu \\|_2^2 + \\nu^Ty \\\\ &\\| \\nu^TX \\|_\\infty \\leq \\lambda \n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Now let us rewrite this as a general Quadratic Problem.\n",
"\n",
"Let $ A = \\begin{pmatrix}X^T \\\\ -X^T \\end{pmatrix} $ and $ b = \\lambda \\mathbb{1}$ (where $\\mathbb{1}$ is the vector of $\\mathbb{R}^{2d}$ full of ones). We have that $\\| \\nu^TX \\|_\\infty \\leq \\lambda \\Leftrightarrow \\| X^T\\nu \\|_\\infty \\leq \\lambda \\Leftrightarrow A\\nu \\leq b$.\n",
"\n",
"Let $Q = \\frac{1}{2} I_n$ and $p=y$. We have that $\\frac{1}{2} \\| \\nu \\|_2^2 = \\frac{1}{2} \\nu^TI_n\\nu$, therefore $\\frac{1}{2} \\| \\nu \\|_2^2 + \\nu^Ty = \\nu^TQ\\nu + p^T\\nu$.\n",
"\n",
"With these notations, the dual problem becomes a QP: \n",
"\n",
"$$ \n",
"\\begin{aligned} \n",
"&\\textrm{min}_{v} \\quad &&\\nu^TQ\\nu + p^T\\nu \\\\ &\\textrm{subject to} && A\\nu \\leq b\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 2\n",
"Implement the barrier method to solve QP."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Write a function v_seq = centering_step(Q,p,A,b,t,v0,eps) which implements the Newton method to solve the centering step given the inputs (Q; p; A; b), the barrier method parameter t (see lectures), initial variable v0 and a target precision eps. The function outputs the sequence of variables iterates (vi). Use a backtracking line search with appropriate parameters."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $\\Phi_t(v) = t(v^TQv+p^Tv) + \\phi(v)$, where $\\phi(v) = -\\sum_{i=1}^{2d} \\log(b_i-A_i v)$ is the logarithmic barrier function (and $A_1, ..., A_{2d}$ are the rows of $A$)\n",
"\n",
"In order to perform the backtracking line search we first have to derive the gradient and Hessian of $\\Phi_t$.\n",
"\n",
"We have, with the notation $d(v)_i = \\frac{1}{b_i-A_iv}$, that $\\nabla\\phi(v) = A^Td(v)$ and $\\nabla^2\\phi(v) = A^Tdiag(d(v))^2A$. Therefore $$ \\nabla\\Phi_t(v) = t(2Qv+p)+A^Td(v) $$ $$ \\nabla^2\\Phi_t(v) = 2tQ + A^Tdiag(d(v))^2A $$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import sklearn"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# backtracking line search parameters, chosen as they are in Boyd's book\n",
"alpha = 0.01\n",
"beta = 0.5\n",
"assert 0 < alpha and alpha < 0.5 and 0 < beta and beta < 1.0\n",
"\n",
"def f(Q,p,v):\n",
" return v.T@Q@v + p.T@v\n",
"\n",
"def phi(A,b,v):\n",
" assert (b-A.dot(v)>0).all()\n",
" return -np.sum(np.log(b-A.dot(v)))\n",
"\n",
"# Phi_t at v\n",
"def phi_t(Q, p, A, b, t, v):\n",
" return t * f(Q,p,v) + phi(A,b,v)\n",
"\n",
"def centering_step(Q,p,A,b,t,v0,eps): \n",
" v_seq = [v0]\n",
" while True:\n",
" v = v_seq[-1]\n",
" # compute Newton step and decrement\n",
" d = 1/(b-A.dot(v))\n",
" D2 = np.diag(d.flatten()**2)\n",
" nabla = t*(2*Q@v+p) + A.T.dot(d) # gradient of Phi_t at v\n",
" nabla2_inv = np.linalg.inv(2*t*Q + A.T @ D2 @ A) # hessian of Phi_t at v\n",
" delta_v_nt = - nabla2_inv @ nabla # newton step\n",
" lambda2 = nabla.T @ nabla2_inv @ nabla # newton decrement\n",
" if lambda2/2 < eps: # stopping criterion\n",
" break \n",
" # backtracking line search to find t\n",
" gamma = 1\n",
" while not((b-A.dot(v+gamma*delta_v_nt)>0).all()) or phi_t(Q,p,A,b,t, v + gamma*delta_v_nt) > phi_t(Q,p,A,b,t,v) + alpha*gamma*nabla.T@delta_v_nt:\n",
" gamma *= beta\n",
" # update v\n",
" v_seq.append(v + gamma*delta_v_nt) \n",
" return v_seq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2. Write a function v seq = barr method(Q,p,A,b,v0,eps) which implements the barrier method to solve QP using precedent function given the data inputs (Q; p; A; b), a feasible point v0, a precision criterion eps. The function outputs the sequence of variables iterates (vi)."
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [],
"source": [
"def barr_method(Q,p,A,b,v0,eps):\n",
" # barrier method parameter\n",
" assert mu>1\n",
" v_seq = [v0]\n",
" niter_seq = [0]\n",
" m = A.shape[0]\n",
" eps_cs = 1e-6\n",
" t = 1\n",
" while True:\n",
" # centering step\n",
" v_seq_cs = centering_step(Q,p,A,b,t,v_seq[-1],eps_cs)\n",
" # update\n",
" v_seq.append(v_seq_cs[-1]) \n",
" niter_seq.append(len(v_seq_cs))\n",
" if m/t < eps: # stopping criterion\n",
" break\n",
" # increase\n",
" t *= mu\n",
" return v_seq, niter_seq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 3\n",
"Test your function on randomly generated matrices $X$ and observations $y$ with $\\lambda = 10$. Repeat for different values of the barrier method parameter $\\mu$ and check the impact on $w$. What would be an appropriate choice for $\\mu$?"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [],
"source": [
"# Generate data\n",
"from sklearn.datasets import make_regression\n",
"n = 100\n",
"d = 10\n",
"\n",
"np.random.seed(42)\n",
"X, y, coef = make_regression(n_samples=n, n_features=d, coef=True, noise = 0.1)\n",
"y = y.reshape((-1,1))"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [],
"source": [
"# Generate the dual QP parameters\n",
"lbd = 10\n",
"A = np.vstack((X.T,-X.T))\n",
"b = lbd * np.ones((2*d,1))\n",
"Q = 0.5*np.eye(n)\n",
"p = y\n",
"v0 = np.zeros((n,1)) # stricly feasible\n",
"eps = 1e-6"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, note that strong duality holds because the primal problem is convex and that $0$ is stricly feasible so Slater's condition holds. The KKT stationarity condition implies that $Xw^\\star-y-v^\\star=0$, therefore, if $X^TX$ is invertible (which is often the case if $n>d$), then $w^\\star=(X^TX)^{-1}X^T(y+v^\\star)$"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [],
"source": [
"K = np.linalg.inv(X.T @ X) @ X.T\n",
"\n",
"def get_w(v):\n",
" return (K @ (y+v)).flatten()\n",
"\n",
"def get_lasso(w):\n",
" return 0.5* np.linalg.norm(X@w - y)**2 + lbd*np.linalg.norm(w, ord=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can plot the gap $f(v_t)-f^*$ for different values of $\\mu$ and check the impact on $w$."
]
},
{
"cell_type": "code",
"execution_count": 88,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Barrier method for mu = 2\n",
"Barrier method for mu = 15\n",
"Barrier method for mu = 50\n",
"Barrier method for mu = 100\n",
"Barrier method for mu = 200\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAE+CAYAAAANqS0iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dfZhV5X3v//cX0IM1YxMVjcxgkII48liD2CTiwRSVoEfjI6AnicApnt/PpGmS/k5se8yDtY05xVNN9ZzWJiLJFfF4TAxJimiuqJVjE1FkGBUZoCqHwQciNBVMqAr374+9xgzzAHtm1szea/b7dV1ce/a9177Xd98w24/3WutekVJCkiRJxTGk0gVIkiSpZwxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGACdJklQwBjhJkqSCGVbpAgbKsccem0aPHl3pMiRJkg5p7dq1r6eURnT3emEDXEQcCTwGfDml9ONDbT969Gieeuqp/i9MkiSpjyJi68Fer5pDqBFxZ0TsiIhnO7TPjoiWiNgSEde1e+mLwL0DW6UkSVLlVU2AA+4CZrdviIihwO3Ax4BTgfkRcWpEzAI2AK8NdJGSJEmVVjWHUFNKj0XE6A7N04EtKaUXACLiHuAi4D3AkZRC3a8jYmVKaf8AlitJklQxVRPgulEPbGv3vBU4I6X0aYCIuBp4vbvwFhGLgcUAJ554Yv9WKkmSePvtt2ltbWXv3r2VLqUQhg8fTkNDA4cddliP3lftAS66aEvv/pDSXQd7c0rpDuAOgGnTpqWDbStJkvqutbWVuro6Ro8eTURX/xlXm5QSO3fupLW1lZNOOqlH762mc+C60gqMave8AXi5QrVIkqRD2Lt3L8ccc4zhrQwRwTHHHNOr2cpqD3BPAuMi4qSIOByYB/ywwjVJkqSDMLyVr7djVTUBLiKWAz8DxkdEa0QsSim9A3waeBB4Hrg3pfRcJeuUJEmqtKo5By6lNL+b9pXAygEuR5IkqWpVzQzcYHDnghu4c8ENlS5DkiQNcga4HKX9x5P2H1/pMiRJUs62bdvG2WefTWNjIxMmTODWW2+taD1VcwhVkiSpWg0bNoybb76Z0047jd27d/PBD36Qc845h1NPPbUi9TgDJ0mSBp2ZM2fS0tICwM6dO5k4cWKf+jvhhBM47bTTAKirq6OxsZHt27f3uc7ecgZOkiT1i6/+6Dk2vPxGrn2eOvIovvwfJhxyuy1btjBu3DgAmpubmTRpUqdtZsyYwe7duzu1L1myhFmzZnXb90svvcS6des444wzelB5vgxwkiRpUNm6dSv19fUMGVI60Njc3MzkyZM7bbd69eoe971nzx4uvfRSbrnlFo466qg+19pbBjhJktQvypkp6w9NTU0HBLa1a9cyd+7cTtv1dAbu7bff5tJLL+Wqq67ikksuybfoHjLASZKkQWX9+vXv3p5q8+bNrFixghtvvLHTdj2ZgUspsWjRIhobG/n85z+fW6295UUMkiRpUGlqamL//v1MmTKFG264gcbGRpYtW9anPh9//HG+853v8PDDDzN16lSmTp3KypWVu8+AM3CSJGlQaW5uZt26ddTV1eXW55lnnklKKbf++soZOEmSNGjs3r2bIUOG5BreqpEBTpIkDRp1dXVs2rSp0mX0OwOcJElSwRjgJEmSCsYAJ0mSVDAGOEmSpIIxwEmSJBWMAU6SJKlgDHCSJEkFY4CTJEkqGAOcJElSGRYuXMhxxx3HxIkTD2gfPXo0kyZNYurUqUybNm1AajHASZIkleHqq69m1apVXb72yCOP0NTUxFNPPTUgtRjgJEnSoDNz5kxaWloA2LlzZ6dZs94466yzOProo/vcTx6GVboASZI0SD1wHbz6TL59vn8SfOymQ262ZcsWxo0bB0BzczOTJk3qtM2MGTPYvXt3p/YlS5Ywa9asskuKCM4991wigmuuuYbFixeX/d7eMsBJkqRBZevWrdTX1zNkSOlAY3NzM5MnT+603erVq3PZ3+OPP87IkSPZsWMH55xzDqeccgpnnXVWLn13xwAnSZL6RxkzZf2hqanpgMC2du1a5s6d22m7vGbgRo4cCcBxxx3HxRdfzJo1awxwkiRJPbF+/Xr27t0LwObNm1mxYgU33nhjp+3ymIF788032b9/P3V1dbz55ps89NBDfOlLX+pzv4fiRQySJGlQaWpqYv/+/UyZMoUbbriBxsZGli1b1ud+58+fz4c+9CFaWlpoaGjgW9/6Fq+99hpnnnkmU6ZMYfr06Zx//vnMnj07h09xcM7ASZKkQaW5uZl169ZRV1eXa7/Lly/vsn39+vW57qcczsBJkqRBY/fu3QwZMiT38FZtDHCSJGnQqKurY9OmTZUuo98Z4CRJkgrGACdJklQwBjhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBWOAkyRJKpjCBriI+HhE/H1ErIiIcytdjyRJ0kCpqgAXEXdGxI6IeLZD++yIaImILRFxHUBK6QcppT8ArgbmVqBcSZJUQ0aPHs2kSZOYOnUq06ZNe7d91apVjB8/nrFjx3LTTTcNSC3Vdi/Uu4DbgG+3NUTEUOB24BygFXgyIn6YUtqQbfJfs9clSZL61SOPPMKxxx777vN9+/Zx7bXX8pOf/ISGhgZOP/10LrzwQk499dR+raOqZuBSSo8Buzo0Twe2pJReSCm9BdwDXBQlXwceSCk9PdC1SpKk6jVz5kxaWloA2LlzJxMnTuyX/axZs4axY8cyZswYDj/8cObNm8eKFSv6ZV/tVdsMXFfqgW3tnrcCZwCfAWYBvx0RY1NKf9vxjRGxGFgMcOKJJw5AqZIkqc3X13ydjbs25trnKUefwhenf/GQ223ZsoVx48YB0NzczKRJkzptM2PGDHbv3t2pfcmSJcyaNatTe0Rw7rnnEhFcc801LF68mO3btzNq1Kh3t2loaOCJJ57oyUfqlSIEuOiiLaWUvgF842BvTCndAdwBMG3atNQPtUmSpCqzdetW6uvrGTKkdKCxubmZyZMnd9pu9erVPer38ccfZ+TIkezYsYNzzjmHU045hZQ6x4uIrqJLvooQ4FqBUe2eNwAvV6gWSZJUpnJmyvpDU1PTAYFt7dq1zJ3b+XrHns7AjRw5EoDjjjuOiy++mDVr1vCRj3yEbdt+c6CwtbX13e36UxEC3JPAuIg4CdgOzAOurGxJkiSpWq1fv569e/cCsHnzZlasWMGNN97YabuezMC9+eab7N+/n7q6Ot58800eeughvvSlL3H66aezefNmXnzxRerr67nnnnu4++67c/ss3amqABcRy4GZwLER0Qp8OaX0rYj4NPAgMBS4M6X0XAXLlCRJVaypqYkjjjiCKVOmMHnyZBobG1m2bBnXX399r/t87bXXuPjiiwF45513uPLKK5k9ezYAt912G+eddx779u1j4cKFTJgwIZfPcTBVFeBSSvO7aV8JrBzgciRJUgE1Nzezbt066urqcutzzJgxrF+/vsvX5syZw5w5c3LbVzmqahkRSZKkvti9ezdDhgzJNbxVIwOcJEkaNOrq6ti0aVOly+h3BjhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGACdJklQwBjhJkqQyLFy4kOOOO46JEyce0L5q1SrGjx/P2LFjuemmmw7ZngcDnCRJUhmuvvpqVq1adUDbvn37uPbaa3nggQfYsGEDy5cvZ8OGDd2258UAJ0mSBp2ZM2fS0tICwM6dOzvNmvXGWWedxdFHH31A25o1axg7dixjxozh8MMPZ968eaxYsaLb9rwMy60nSZKkdl79y7/k357fmGuf/67xFN7/p396yO22bNnCuHHjAGhubmbSpEmdtpkxYwa7d+/u1L5kyRJmzZpVVj3bt29n1KhR7z5vaGjgiSee6LY9LwY4SZI0qGzdupX6+nqGDCkdaGxubmby5Mmdtlu9enWf95VS6tQWEd2258UAJ0mS+kU5M2X9oamp6YDAtnbtWubOndtpuzxm4BoaGti2bdu7z1tbWxk5cmS37XkxwEmSpEFl/fr17N27F4DNmzezYsUKbrzxxk7b5TEDd/rpp7N582ZefPFF6uvrueeee7j77rsZP358l+15McDl7O3D6vnWp/6u0mXoIIa+dytX3/qXlS5DktRPmpqaOOKII5gyZQqTJ0+msbGRZcuWcf311/ep3/nz5/Poo4/y+uuv09DQwFe/+lUWLVrEbbfdxnnnnce+fftYuHAhEyZMAOi2PQ8GuBwNfe9W+GWlq9DBvH1YvX9HkjTINTc3s27dOurq6nLtd/ny5V22z5kzhzlz5pTdngcDXI6c1al+zo5K0uC2e/duhgwZknt4qzauAydJkgaNuro6Nm3aVOky+p0BTpIkqWAMcJIkSQVjgJMkSSoYL2JQ7Ukw9+9+llt3F02t58ozTsytP0mSDsUZONWe/O5kwoZX3mBF0/b8OpQkqQzOwKkm/a9rPpRLP3nO5EmSVC5n4CRJkgrGACdJklQwBjhJkqRD2LZtG2effTaNjY1MmDCBW2+99d3XVq1axfjx4xk7diw33XTTIdvzYICTJEk6hGHDhnHzzTfz/PPP8/Of/5zbb7+dDRs2sG/fPq699loeeOABNmzYwPLlyw/anhcDnCRJGnRmzpxJS0sLADt37mTixIl96u+EE07gtNNOA0q362psbGT79u2sWbOGsWPHMmbMGA4//HDmzZvHihUrum3Pi1ehSn204ZU3XFdOkrqw+t5NvL5tT659HjvqPcy44uRDbrdlyxbGjRsHQHNzM5MmTeq0zYwZM9i9e3en9iVLljBr1qxu+37ppZdYt24dZ5xxBg899BCjRo1697WGhgaeeOIJtm/f3mV7XgxwUh9cNLU+1/42vPIGgAFOkvpg69at1NfXM2RI6UBjc3MzkydP7rTd6tWre9z3nj17uPTSS7nllls46qijSCl12iYium3PiwFO6oMrzzgx17DlunKSBpNyZsr6Q1NT0wGBbe3atcydO7fTdj2dgXv77be59NJLueqqq7jkkkuA0szatm3b3t2mtbWVkSNHdtueFwOcJEkaVNavX8/evXsB2Lx5MytWrODGG2/stF1PZuBSSixatIjGxkY+//nPv9t++umns3nzZl588UXq6+u55557uPvuuxk/fnyX7XnxIgZJkjSoNDU1sX//fqZMmcINN9xAY2Mjy5Yt61Ofjz/+ON/5znd4+OGHmTp1KlOnTmXlypUMGzaM2267jfPOO4/GxkauuOIKJkyY0G17XpyBkyRJg0pzczPr1q2jrq4utz7PPPPMLs9rA5gzZw5z5swpuz0PzsBJkqRBY/fu3QwZMiTX8FaNCjkDFxFHAv8DeAt4NKX03QqXJOUmz2VJXJJEUq2pq6tj06ZNlS6j31XNDFxE3BkROyLi2Q7tsyOiJSK2RMR1WfMlwH0ppT8ALhzwYqV+ctHUek494ahc+trwyhusaNqeS1+SpOpSTTNwdwG3Ad9ua4iIocDtwDlAK/BkRPwQaACeyTbbN7BlSv0nz2VJXJJEkgavqpmBSyk9Buzq0Dwd2JJSeiGl9BZwD3ARpTDXkG3T7WeIiMUR8VREPPWLX/yiP8qWJEkddHeyvzrr7VhVTYDrRj2wrd3z1qzt+8ClEfE/gR919+aU0h0ppWkppWkjRozo30olSRLDhw9n586dhrgypJTYuXMnw4cP7/F7q+kQale6uudESim9CSwY6GIkSdLBNTQ00Nraike+yjN8+HAaGhoOvWEH1R7gWoFR7Z43AC9XqBZJknQIhx12GCeddFKlyxj0qj3APQmMi4iTgO3APODKypak4kuw9PzKljDpMpjW/5PIeS5JAi5LIknVomrOgYuI5cDPgPER0RoRi1JK7wCfBh4EngfuTSk9V8k6VXARdH1kfgC9+gw8c1+/7ybPJUnAZUkkqZpUzQxcSml+N+0rgZUDXI4GrSGl/LbgHypXwgDN/uW5JAm4LIkkVZOqmYGTJElSeQxwkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgqmaq5ClVT9erOunGvHSVL+DHCSynLR1Poev2fDK28AGOAkKWcGOEll6c26cq4dJ0n945DnwEXEkdnje/q/HEmSJB1KOTNw74uIBcAWYFU/1yP1u/1pHwtW9f99SLsVr5Uec6xhzpg5XH7y5bn1J0mqbuVchfr7wNXAmIg4rn/LkfrXYUOHMSSGVrqMXLXsamHlC95tTpJqSTkzcGuAhcColNKOfq5H6leHDTmcw4YcztLZSytXRNu9UHOqoaKziZKkijhkgEspPZ/92NzPtUgahFx6RJLyV9ZCvhFxa/Z4RP+WI2kwuWhqPaeecFSP3rPhlTdY0bS9nyqSpMGh3GVEfj97/D/AB/upFkmDjEuPSFL/KPdWWqsi4mfA+yNiYUR8MCKG92dhkiRJ6lpZM3AppT+OiDHAo8BJwIXAhIh4C3g2pTS3/0qUJElSe2XfiSGl9EJEzEopbWpryxb3ndgvlUn9ZP+vfsXWT3yycgW8+nLp8eF8api3ayPPf/BYmJ1Ld5KkAujRrbTah7fs+R7g57lWJPWjYcccwzsAeytdSX6O2/4r4PVKlyFJGkDeC1U1ZdiIEQwbMYIPfOHKyhXRtg7cgm/n0t3G86fn0k816c3SI+DyI5JqR7kXMRwgIv5D3oVIEvRu6RFw+RFJtaW3M3B/Afwoz0IkCXq39Ai4/Iik2tKrGTggcq1CkiRJZettgEu5ViFJkqSyeRGDas7rrXu4/+anK1fAq1eUHnOqYdeIxQCV/UzAydOPZ8KM+orWIEm1wgCnmnLy9OMrXULJW2/Cq8/k01faV3rsaX9HjoC69+dSwuutewAMcJI0QHob4F7LtQppgEyYUV/5kPHUenjmH3Lr7sGmlwA4b+bo8t/06jPw/kmwIJ86Kj37J0m1plcBLqV0Tt6FSDVj2oLSn7zcl60D15Mw1rYW3SDTm/XjXDtOUhF5CFXSoHDR1J7PrG545Q0AA5ykwjHASRoUerN+nGvHSSqq3t6J4ciIGJp3MZIkSTq0sgJcRAyJiCsj4h8iYgewEXglIp6LiL+KiHH9W6YkSZLalDsD9wjwO8CfAO9PKY1KKR0HzAB+DtwUEf+xn2qUJElSO+WeAzcrpfR2x8aU0i7ge8D3IuKwXCuTJElSl8oKcF2FN4CI+HxK6b9nT8cALXkVJkkDwaVHJBVRr65CjYj3An8NjI+IvUAzsAjIcXErSepfLj0iqah6u5DvL4EFEXEe8DowGfh+noVJUn9z6RFJRdWjABcR/wT8WUrpEYCU0oPZS2vzLkySJEld6+k6cIuBT0fETyPiQ/1RULki4uMR8fcRsSIizq1kLZIkSQOpRwEupfRsSulS4P8D/iwifhwRU3u604i4MyJ2RMSzHdpnR0RLRGyJiOsOUcsPUkp/AFwNzO1pDZIkSUXVqzsxAFuAPwdagad68f67gNntG7I7O9wOfAw4FZgfEadGxKQsKLb/c1y7t/7X7H2SJEk1oafnwD0MjAP2AhuyP1f3dKcppcciYnSH5unAlpTSC9m+7gEuSil9Dbigi1oCuAl4IKX0dE9rkKTecukRSZXW06tQ/xh4PqX0636opR7Y1u55K3DGQbb/DDAL+O2IGJtS+tuOG0TEYkrn7XHiiX5xSuo7lx6RVA3KCnAREamk25mutm36UEt00dZtfymlbwDfOFiHKaU7gDsApk2b1pfaJAlw6RFJ1aHcc+AejojPRMQB31oRcXhEfDQilgGf6mMtrcCods8bgJf72KckSdKgU+4h1M3APuD+iDgB+CUwHBgKPAT8dUqpqY+1PAmMi4iTgO3APODKPvYpSZI06JQb4D6cUlocEf8JOBEYAfw6uyNDj0XEcmAmcGxEtAJfTil9KyI+DTxIKRjemVJ6rjf9S5IkDWblBrgHI+JnwPHAJ4H1QK/DVUppfjftK4GVve1XkiSpFpQV4FJKX4iIMcCjwEnAhcCEiHgLeDal5EK6knQQLj0iKU9lLyOSUnohImallDa1tUXEe4CJ/VKZJA0SLj0iKW89WgeufXjLnu8Bfp5rRZJ67Ffv/JoFqxaU/4Z4rfTYk/ccxPhd5wCwYNXfdHptzpg5XH7y5bnsp6hcekRS3np7Ky1JVeKYI47ht4YdUekyutSyq4WVL3haqyTlrad3YpBUZUYcMYIRR4xg6eyl5b9p6fmlx5685yDuf660xvd1sw9c+adHs4KSpLI5AydJklQwzsBJterVZ34zE9fnvq4oPS69/sD2tnPtutvPpMtgmrN0ktRTBjipFk26rNIVlAIkGOAkqRcMcFItmrYg3+B0c+kcOBZ0uCVy2zlwXZ1rl9fs3yDWm7XjwPXjpFpggJOkKtSbtePA9eOkWmGAk6Qq1Ju148D146Ra4VWokiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGixgkaZDpzfIjLj0iFYsBTpIGkd4sP+LSI1LxGOCkQWDvxo1s/cQnK7f/4bMB2PqJWw5on7drY6n9u13U9urLpceHO7921AUX8L65V+RbZI3ozfIjLj0iFY8BTiq4oy64oNIl5GrvxlLoM8BJUvcMcFLBvW/uFRUPO09nt9L6wBeuPKD9K9mttJYe7FZaC759QHMlZxIlqSi8ClWSJKlgDHCSJEkF4yFUSZJLj0gFY4CTpBrn0iNS8RjgJKnGufSIVDyeAydJklQwBjhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBeNVqJKkXnHtOKlyDHCSpB5z7TipsgxwknLxeuse7s9uat9m/K5zALj/uac7v+HVK0qPHd6zd/hsAJ6+uYv3DKCTpx/PhBk9Dym1wrXjpMoywEnqs5OnH1/pEnL1euseAAOcpKplgJPUZxNm1HcZdhas+hsArpt9Zec3Lb0eXn0Gjp50QPPWl14G4AMfHpl/oV2ZdBlMW3BAU8eZREmqNgY4SZUx6bJKV1AKkNApwElStTPASaqMaQu6Dk4Pf7L0uODb/V/D0vP7fx+S1A8McJKkAePSI1I+DHCSpAHh0iNSfgod4CLiSOAx4MsppR9Xuh5JUvdcekTKT0VupRURd0bEjoh4tkP77IhoiYgtEXFdGV19Ebi3f6qUJEmqTpWagbsLuA149yzliBgK3A6cA7QCT0bED4GhwNc6vH8hMBnYAAwfgHolSZKqRkUCXErpsYgY3aF5OrAlpfQCQETcA1yUUvoacEHHPiLibOBI4FTg1xGxMqW0v8M2i4HFACee6PkTkiRpcKimc+DqgW3tnrcCZ3S3cUrpzwAi4mrg9Y7hLdvmDuAOgGnTpqU8i5UkSaqUagpw0UXbIUNXSumu/EuRlJeWXS0sWFX+Qrnzdm0E4Cs9eE+vxWvMSUdyef/vSX3g0iNSZ9UU4FqBUe2eNwAvV6gWSTmYM2ZOpUs4qBbegsAAV8VcekTqWjUFuCeBcRFxErAdmAd0cQNFSUVx+cmXc/nJPYtHW79buhPD0tlL+6OkAyy4a1q/70N949IjUtcqtYzIcuBnwPiIaI2IRSmld4BPAw8CzwP3ppSeq0R9kiRJ1axSV6HO76Z9JbBygMuRJEkqlIrMwEmSJKn3DHCSJEkFU00XMUiSlIveLD0CLj+i4jDASZIGld4sPQIuP6JiMcBJkgaV3iw9Ai4/omLxHDhJkqSCMcBJkiQVjAFOkiSpYAxwkiRJBWOAkyRJKhivQpUkKdOb9eNcO06VYICTJInerR/n2nGqFAOcJEn0bv04145TpXgOnCRJUsE4Ayeptr31Jiw9/8C2V68oPS69fmBqmHQZTFswMPuSNCg4Ayepdh05Ag4/srI1vPoMPHNfZWuQVDjOwEmqXXXvL/2ZvfTA9pufLj0u+FT/19Bx9k+SymCAk1R19m7cyNZPfLLf9zNv10YAtn73wH3tHT671P6JW3Lb11EXXMD75l6RW3+qHi49okowwEmqKkddcEGlS8jd3o2loGiAG3xcekSVYoCTVFXeN/eKAQs6X1lVunBgaYdDqE9nh1A/8IUrc9nPQMwmqjJcekSV4kUMkiRJBWOAkyRJKhgDnCRJUsEY4CRJkgrGACdJklQwXoUqSdIAc+049ZUBTpKkAeTaccqDAU6SpAHk2nHKg+fASZIkFYwBTpIkqWAMcJIkSQVjgJMkSSoYA5wkSVLBeBWqJEkF4Npxas8AJ0lSlXPtOHVkgJMkqcq5dpw6MsBJqmktu1pYsGrBAW3jd50DwIJVf5PLPubt2gjAVzrsB4B4rfTYxWtzxszh8pMvz6UGSYOLAU5SzZozZk6lS+hWy64WAAOcpC4VNsBFxBDgz4GjgKdSSssqXJKkgrn85Mu7DEj3P/c0r7fu4cLnPpPLfnaNLM3AXfjcKZ1ffPWZ0uPOSQc0t2Szdvc/93QuNQCcPP14Jszo+blUkqpPRZYRiYg7I2JHRDzboX12RLRExJaIuO4Q3VwE1ANvA639Vauk2nPy9OM5tuE9lS4jV6+37mHTmtcqXYaknFRqBu4u4Dbg220NETEUuB04h1IgezIifggMBb7W4f0LgfHAz1JKfxcR9wE/HYC6JdWACTPqc52p2vqJWwD4wBeu7Pzi0utLjws+dUBz2/l3183u4j29cP/N+c3kSaq8igS4lNJjETG6Q/N0YEtK6QWAiLgHuCil9DXggo59REQr8Fb2dF//VStJklRdqukcuHpgW7vnrcAZB9n++8DfRMQM4LGuNoiIxcBigBNPdB0cSVJtcfHfwauaAlx00Za62zil9Ctg0cE6TCndAdwBMG3atG77kqSKevUZWHr+gW1ty4t0bO/1Pq7I+ru+69cnXQbTuljmRIXl4r+DWzUFuFZgVLvnDcDLFapFkgbGpMsqXcFvroQ1wA0qLv47uFVTgHsSGBcRJwHbgXlAPmfvSlK1mrag6+DUtrDv7KX57KftIoYOF0sA+c3ySRowlVpGZDnwM2B8RLRGxKKU0jvAp4EHgeeBe1NKz1WiPkmSpGpWqatQ53fTvhJYOcDlSJIkFUpFZuAkSZLUewY4SZKkgjHASZIkFYwBTpIkqWAMcJIkSQVTTevASZKkCuvN7bfAW3ANNOb5+UwAAAsLSURBVAOcJEkCenf7LfAWXJVggJMkSUDvbr8F3oKrEjwHTpIkqWAMcJIkSQVjgJMkSSoYA5wkSVLBGOAkSZIKxgAnSZJUMAY4SZKkgjHASZIkFYwBTpIkqWAMcJIkSQVjgJMkSSqYSClVuoYBERG/ALYOwK6OBV4fgP1UO8fBMWjjOJQ4DiWOQ4njUOI4lHQ1Dh9IKY3o7g01E+AGSkQ8lVKaVuk6Ks1xcAzaOA4ljkOJ41DiOJQ4DiW9GQcPoUqSJBWMAU6SJKlgDHD5u6PSBVQJx8ExaOM4lDgOJY5DieNQ4jiU9HgcPAdOkiSpYJyBkyRJKhgDXE4iYnZEtETEloi4rtL1DJSIuDMidkTEs+3ajo6In0TE5uzxfZWscSBExKiIeCQino+I5yLis1l7TY1FRAyPiDURsT4bh69m7SdFxBPZOPyviDi80rX2t4gYGhHrIuLH2fOaGwOAiHgpIp6JiKaIeCprq7Xfi/dGxH0RsTH7jvhQDY7B+OzfQNufNyLij2ptHAAi4nPZ9+OzEbE8+97s8feDAS4HETEUuB34GHAqMD8iTq1sVQPmLmB2h7brgJ+mlMYBP82eD3bvAF9IKTUCvwdcm/0bqLWx+DfgoymlKcBUYHZE/B7wdeCvs3H4F2BRBWscKJ8Fnm/3vBbHoM3ZKaWp7ZZJqLXfi1uBVSmlU4AplP5d1NQYpJRasn8DU4EPAr8C7qfGxiEi6oE/BKallCYCQ4F59OL7wQCXj+nAlpTSCymlt4B7gIsqXNOASCk9Buzq0HwRsCz7eRnw8QEtqgJSSq+klJ7Oft5N6Qu6nhobi1SyJ3t6WPYnAR8F7svaB/04REQDcD7wzex5UGNjcAg183sREUcBZwHfAkgpvZVS+iU1NAZd+H3gn1NKW6nNcRgGHBERw4DfAl6hF98PBrh81APb2j1vzdpq1fEppVegFGyA4ypcz4CKiNHA7wJPUINjkR06bAJ2AD8B/hn4ZUrpnWyTWvj9uAX4L8D+7Pkx1N4YtEnAQxGxNiIWZ2219HsxBvgFsDQ7pP7NiDiS2hqDjuYBy7Ofa2ocUkrbgSXA/6UU3P4VWEsvvh8McPmILtq8vLcGRcR7gO8Bf5RSeqPS9VRCSmlfdpikgdLsdGNXmw1sVQMnIi4AdqSU1rZv7mLTQTsGHXwkpXQapVNMro2Isypd0AAbBpwG/M+U0u8CbzLIDxMeTHZu14XA/650LZWQneN3EXASMBI4ktLvRkeH/H4wwOWjFRjV7nkD8HKFaqkGr0XECQDZ444K1zMgIuIwSuHtuyml72fNNTkWANlhokcpnRP43uxwAQz+34+PABdGxEuUTqf4KKUZuVoag3ellF7OHndQOudpOrX1e9EKtKaUnsie30cp0NXSGLT3MeDplNJr2fNaG4dZwIsppV+klN4Gvg98mF58Pxjg8vEkMC67iuRwStPDP6xwTZX0Q+BT2c+fAlZUsJYBkZ3j9C3g+ZTSf2/3Uk2NRUSMiIj3Zj8fQenL6nngEeCybLNBPQ4ppT9JKTWklEZT+i54OKV0FTU0Bm0i4siIqGv7GTgXeJYa+r1IKb0KbIuI8VnT7wMbqKEx6GA+vzl8CrU3Dv8X+L2I+K3svxtt/x56/P3gQr45iYg5lP4veyhwZ0rpLypc0oCIiOXATOBY4DXgy8APgHuBEyn9Y708pdTxQodBJSLOBFYDz/Cb857+lNJ5cDUzFhExmdIJuEMp/Q/ivSmlGyJiDKXZqKOBdcB/TCn9W+UqHRgRMRP445TSBbU4Btlnvj97Ogy4O6X0FxFxDLX1ezGV0gUthwMvAAvIfj+okTEAiIjfonS++JiU0r9mbTX1bwEgW15pLqXVC9YB/4nSOW89+n4wwEmSJBWMh1AlSZIKxgAnSZJUMAY4SZKkgjHASZIkFYwBTpIkqWAMcJIOKSJSRNzc7vkfR8RXcur7roi47NBb9nk/l0fE8xHxSIf20dnn+0y7ttsi4uo+7m9mRHy4L3106O+fssfREXFlXv1mff5pV/uSVL0McJLK8W/AJRFxbKULaS8ihvZg80XA/5tSOruL13YAn80W4s7LTEorrOcipdTW12igRwGujHE6IMC125ekKmWAk1SOd4A7gM91fKHjDFpE7MkeZ0bEP0bEvRGxKSJuioirImJNRDwTEb/TrptZEbE62+6C7P1DI+KvIuLJiGiOiGva9ftIRNxNaeHkjvXMz/p/NiK+nrV9CTgT+NuI+KsuPt8vgJ/ymxXh2/f3OxGxKrsZ++qIOCWr7YUoeW9E7G+7x2e2zVjgPwOfi4imiJgRER+IiJ9mn+WnEXFiu/H7RkT8U9Znl7ORbeMK3ATMyPr9XE/GKSJ+kH2O5yK7sXxE3AQckfX33Q5/h5H1/Ww2pnPb9f1oRNwXERsj4rsREW39RcSGrJYlXX0WSX037NCbSBIAtwPNEfHfevCeKZRuZr+L0gr030wpTY+IzwKfAf4o22408O+B3wEeyQLQJ4F/TSmdHhH/Dng8Ih7Ktp8OTEwpvdh+ZxExEvg68EHgX4CHIuLj2d0gPkrprghPdVPrTcADEXFnh/Y7gP+cUtocEWcA/yOl9NGI2AScSumm1GsphaongIaU0paI+FtgT0ppSVbbj4Bvp5SWRcRC4BvAx7N9nEApYJ5C6dZC9x1kTK/LPkdb0F3cg3FamFLaFaXbnD0ZEd9LKV0XEZ9OKU3tYl+XAFMp/T0em73nsey13wUmULpn4+PARyJiA3AxcEpKKUV2WzVJ+XMGTlJZUkpvAN8G/rAHb3sypfRKdkuYfwbagsUzlEJbm3tTSvtTSpspBb1TKN0385MR0UTplmTHAOOy7dd0DG+Z04FHsxtFvwN8FzirzM/3IrCGdocnI+I9lA6D/u+sjr+jFLagdOu0s7I/X6MUwE6ndG/krnwIuDv7+TvZ9m1+kH3+DcDx5dTbTk/G6Q8jYj3wc2BUu+26cyawPKW0L7v5+D9S+oxtfbemlPYDTZT+Pt8A9gLfjIhLgF/18LNIKpMBTlJP3ELpXLIj27W9Q/Zdkh1Ga38eWft7+e1v93w/Bx4B6HhPvwQE8JmU0tTsz0kppbYA+GY39UW5H6Qbfwl8kd98Nw4BftmuhqkppcbstdXADEqzXCuB91I67+0xytP+M7cfp55+hrLGKUr3ZZ0FfCilNIXS/RaHl9F3d9rXvA8YloXm6cD3KM0ururRJ5FUNgOcpLJlN5m+l1KIa/MSpUOWABcBh/Wi68sjYkh2XtwYoAV4EPh/IuIwgIg4OSKOPFgnlGag/n1EHBulE/fnU5o1KktKaSOwAbgge/4G8GJEXJ7VEBExpd2+PgzsTyntpTQLdQ2lYAewG6hr1/0/AfOyn68C/k+5dXXQsd9yx+m3gX9JKf0qIk4Bfq/da2+3vb+Dx4C52Xl2IyjNNq7prrBsxvK3U0orKR0e7+qwrKQcGOAk9dTNlM6HavP3lELTGuAMup8dO5gWSkHrAUrnm+0FvkkpTD0dEc9SOnx50PN2U0qvAH8CPAKsB55OKa3oYS1/ATS0e34VsCg79PgcpZBKdlh4G6XDkVAKbnX85sKKHwEXt13EQOnQ84KIaAY+AXy2h3W1aQbeiYj1EfE5yh+nVcCwbP9/3q5uKJ3n19x2EUM792f7Ww88DPyXlNKrB6mtDvhxto9/pIuLXiTlI1LqeORCkiRJ1cwZOEmSpIIxwEmSJBWMAU6SJKlgDHCSJEkFY4CTJEkqGAOcJElSwRjgJEmSCsYAJ0mSVDD/P7NrLrI2RXcSAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10,5)) \n",
"plt.xlabel('Number of Newton iterations')\n",
"plt.ylabel('$f(v_t)-f^*$')\n",
"\n",
"w_seq_mu = {}\n",
"\n",
"for mu in [2,15,50,100,200]:\n",
" print('Barrier method for mu =',mu)\n",
" v_seq, niter_seq = barr_method(Q,p,A,b,v0,eps)\n",
" w_seq_mu[mu] = [get_w(v) for v in v_seq]\n",
" y_axis = [(f(Q,p,v) - f(Q,p,v_seq[-1])).flatten() for v in v_seq]\n",
" x_axis = np.array(niter_seq).cumsum()\n",
" plt.step(x_axis, y_axis, label='$\\mu$ = '+str(mu))\n",
"\n",
"plt.semilogy()\n",
"plt.legend(loc = 'upper right')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAocAAAE9CAYAAAB5rchsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5xV1X338c+XmxodTFRMgMGCgSAgF+OAWi8hiXihjQRvYKzRqMG02Gq0faJpY6yPMbFV45NoWk3RUKsiNVpINaiJNjHGgiDDyEVgKqKDRAkaHU1RLr/nj704OQ5z2eA5c2aG7/v1mtecvfZaa//OOYz5Za299lJEYGZmZmYG0K3SAZiZmZlZx+Hk0MzMzMwKnByamZmZWYGTQzMzMzMrcHJoZmZmZgVODs3MzMysoEelA+gqDjjggBg4cGClwzAzMzNr06JFi34bEX2aO+fksEQGDhzIwoULKx2GmZmZWZskrW3pnKeVzczMzKzAyaGZmZmZFZQ1OZR0kqSVkuolXdHM+T0k3ZfOz5c0sOjclal8paQTU9kASU9IWiFpmaRLiuqfkcq2SaopKu8l6U5Jz0laIml80bn/Sv3Xpp8D24rLzMzMrCsr2z2HkroDtwITgAbgGUlzI2J5UbULgDciYrCkqcD1wBRJw4GpwAigH/AzSZ8AtgCXR8SzkqqARZIeS30uBU4FbmsSypcBImJkSv5+KmlsRGxL58+OiKY3CzYbVwk+FjMzM9tFmzdvpqGhgU2bNlU6lE5jzz33pLq6mp49e+ZuU84FKeOA+oh4AUDSLGASUJwcTgKuTq/vB26RpFQ+KyLeBdZIqgfGRcTTwHqAiGiUtALoDyyPiBXpOk3jGA78PLV5TdLvgBpgQSuxNxtXRMROfQJmZmZWMg0NDVRVVTFw4MDm/vfemogINm7cSENDA4MGDcrdrpzTyv2Bl4uOG1JZs3UiYgvwJrB/nrZpqvcwYH4bcSwBJknqIWkQcDgwoOj8nWlK+Rv6w7+0luJ6H0nTJC2UtHDDhg1thGFmZmYfxKZNm9h///2dGOYkif3333+nR1rLmRw29801HXlrqU6rbSXtA/wYuDQi3mojjjvIksuFwM3Ar8mmpyGbUh4JHJt+ztmJ2ImI2yOiJiJq+vRp9lFBZmZmVkJODHfOrnxe5UwOG3j/CF018EpLdST1APYFXm+traSeZInh3RHxQFtBRMSWiPhqRIyJiEnAh4HV6dy69LsRuIdsKry1uMzMzMy6tHImh88AQyQNktSLbIHJ3CZ15gLnptenA4+n+/rmAlPTquFBwBBgQZr2nQGsiIib8gQh6UOS9k6vJwBbImJ5mmY+IJX3BP6UbFFLa3GZmZmZdWllW5ASEVskXQw8AnQH7oiIZZKuARZGxFyyRO+utODkdbIEklRvNtnilS3A9IjYKukYsqnf5yTVpkt9PSIeljQZ+D7QB3hIUm1EnAgcCDwiaRuwjj9MHe+Rynum+H4G/DCdazaurub+v/8+b67Zq4Q9boNOkENv676OXx/3cMn6O4RefC0+UrL+Su5jI+Hk71Q6inZzz/yXmFO7rtJhmFkZTD9sL3pteLvSYZTcK+sa+JuLp/Hb115F3bpxznkXcNUVf12xeMq6fV5EPAw83KTsqqLXm4AzWmj7LeBbTcp+RfP3AxIRDwIPNlP+IjC0mfJ3yBanNNdXi3F1JW+u2YPNPT5Gz80l+h/S7Ymh7wfpMF5t3MSiDeuZ+dLTlQ6l3cxfk90BcsSg/SociZlZPj169ODKv7+OQ0eN4e23G5k84ThOP2Uiw4cPr0w8Fbmq7bRlT65j1YJXS9rnez370eu9Bi4Y/5PSdTrydKj5Uun6K5Mvc02lQ2hWqUe95r+SEqV9S9Zlh3fEoP2YNKY/XzjioEqHYmYltmLFCj7eZ59Kh8H48eO57bbbGDp0KBs3buRTn/oUS5cubbthCz7eZzAcOjg76LMPI0cMZ926dU4OrXXLH1rK628EvbeVbl3MPm+/xR7vLoIvPVSyPu2DmVO7juXr32J4394l6c+JkplZ6dXX1zNkyBAA6urqGDly5A51jj32WBobG3cov+GGGzj++ONb7PvFF19k8eLFHHHEEaULeCc5OewktmzcyD6//z1/3O2XJevz+d8uYsUI/xP4IEo90rc9MbzvoqNK1qeZWVf09z9ZxvJX2nqa3c4Z3q833/zciFbrrF27lv79+9OtW7amt66ujlGjRu1Q78knn9zp67/99tucdtpp3HzzzfTuXZpBgl3hzKCT2LztPTb3hKundi9Znyt/052h5N9Ox3ZU6pG+4X17M2lM02fFm5lZR1FbW/u+ZHDRokVMmbLjDrs7O3K4efNmTjvtNM4++2xOPfXU0ga9k5wcdhKbt25hW2wtaZ9D6cXE2Lukfe6OPNJnZtb+2hrhK5clS5YUdhxZvXo1c+bM4dprr92h3s6MHEYEF1xwAcOGDeOyyy4rWay7yslhJ9JN3bnzpDtL1+Gdf1K6vszMzHYDtbW17LXXXowePZpRo0YxbNgwZs6cyTe+8Y1d7vOpp57irrvuYuTIkYwZMwaA6667jokTJ5Yq7J3i5NB2K+W6R9DMzHYPdXV1LF68mKqqqpL1ecwxx9CR9too5w4pZh3O9nsES8X3CJqZ7T4aGxvp1q1bSRPDjsgjh7bb8T2CZma2K6qqqli1alWlwyg7jxyamZmZWYGTQzMzMzMrcHJoZmZmZgVODs3MzMyswMmhmZmZmRV4tbJ1WKV+JiH4uYRmZmZtcXLYaWyDiNLuavKb5+BjI0vXX4mVet9i8HMJzczM2uLksLOIAEr89PSPjYSRp5e2zxLzMwnNzMzal5PDTkXwpYcqHYSZmZmV2Pnnn89//ud/cuCBB7J06dJC+cCBA6mqqqJ79+706NGDhQsXlj0WL0gxMzMzq7DzzjuPefPmNXvuiSeeoLa2tl0SQ3ByaGZmZrZTxo8fz8qVKwHYuHEjhx566Afu87jjjmO//fb7wP2UgqeVzczMzHZCfX09Q4YMAaCuro6RI3dc3HnsscfS2Ni4Q/kNN9zA8ccfn/takjjhhBOQxEUXXcS0adN2PfCcnByamZlZ5/PTK7KnbpTSx0bCyd9ptcratWvp378/3bplk691dXWMGjVqh3pPPvlkSUJ66qmn6NevH6+99hoTJkzgkEMO4bjjjitJ3y1xcmhmZmaWU21t7fuSwUWLFjFlypQd6pVq5LBfv34AHHjggUyePJkFCxY4OTQzMzPbQRsjfOWyZMkSNm3aBMDq1auZM2cO11577Q71SjFy+M4777Bt2zaqqqp45513ePTRR7nqqqs+cL9t8YIUMzMzs5xqa2vZtm0bo0eP5pprrmHYsGHMnDnzA/d71llncdRRR7Fy5Uqqq6uZMWMGr776KscccwyjR49m3Lhx/Mmf/AknnXRSCd5F6zxyaGZmZpZTXV0dixcvpqqqqqT93nvvvc2WL1mypKTXycMjh2ZmZmY5NDY20q1bt5Inhh2Nk0MzMzOzHKqqqli1alWlwyi7siaHkk6StFJSvaQrmjm/h6T70vn5kgYWnbsyla+UdGIqGyDpCUkrJC2TdElR/TNS2TZJNUXlvSTdKek5SUskjU/lH5L0kKTnU7vvFLU5T9IGSbXp58KyfEBmZmZmHUzZkkNJ3YFbgZOB4cBZkoY3qXYB8EZEDAa+C1yf2g4HpgIjgJOAH6T+tgCXR8Qw4EhgelGfS4FTgV82ucaXASJiJDABuFHS9vd9Q0QcAhwGHC3p5KJ290XEmPTzLx/kszAzMzPrLMq5IGUcUB8RLwBImgVMApYX1ZkEXJ1e3w/cIkmpfFZEvAuskVQPjIuIp4H1ABHRKGkF0B9YHhEr0nWaxjEc+Hlq85qk3wE1EbEAeCKVvyfpWaC6hO9/t3PP/JeYU7uuZP0tX/8Ww/v2Lll/ZmZm1rZyTiv3B14uOm5IZc3WiYgtwJvA/nnapinow4D5bcSxBJgkqYekQcDhwIAmfX0Y+BwpiUxOk1Qn6X5J76tvzZtTu47l698qWX/D+/Zm0pim/2TMzMysnMo5crjDEB4QOeu02lbSPsCPgUsjoq1s5A5gGLAQWAv8mmx6entfPYB7ge9tH+UEfgLcGxHvSvoKMBP4TNOOJU0DpgEcdNBBbYSxexjetzf3XXRUpcMwMzOzXVTOkcMG3j9CVw280lKdlKTtC7zeWltJPckSw7sj4oG2goiILRHx1XTv4CTgw8Dqoiq3A6sj4uaiNhvTlDbAD8lGG5vr+/aIqImImj59+rQVipmZmVmHV87k8BlgiKRBknqRLTCZ26TOXODc9Pp04PGIiFQ+Na1mHgQMARak+xFnACsi4qY8QaRVyXun1xOALRGxPB1fS5aQXtqkTd+iw1OAFXnftJmZmVlnVrZp5YjYIuli4BGgO3BHRCyTdA2wMCLmkiV6d6UFJ6+TJZCkerPJFq9sAaZHxFZJxwDnAM9Jqk2X+npEPCxpMvB9oA/wkKTaiDgROBB4RNI2YF1qj6Rq4G+B54Fn00KWW9LK5L+SdEq69uvAeeX6nMzMzMw6krJunxcRDwMPNym7quj1JuCMFtp+C/hWk7Jf0fz9iETEg8CDzZS/CAxtpryhlb6uBK5s7pyZmZlZqQ0cOJCqqiq6d+9Ojx49WLhwIQDz5s3jkksuYevWrVx44YVcccUOj40uOe+tbGZmZtYBPPHEExxwwAGF461btzJ9+nQee+wxqqurGTt2LKeccgrDhzd9bHRpefs8MzMzs50wfvx4Vq5cCcDGjRs59NBDy3KdBQsWMHjwYA4++GB69erF1KlTmTNnTlmuVcwjh2ZmZmY7ob6+niFDhgBQV1fHyJEjd6hz7LHH0tjYuEP5DTfcwPHHH79DuSROOOEEJHHRRRcxbdo01q1bx4ABf3h4S3V1NfPnt/V45w/OyaGZmZl1OtcvuJ7nX3++pH0est8hfG3c11qts3btWvr370+3btnka11dHaNGjdqh3pNPPrlT137qqafo168fr732GhMmTOCQQw4he4DL+zWzE1zJOTk0MzMzy6m2tvZ9yeCiRYuYMmXKDvV2duSwX79+ABx44IFMnjyZBQsWcPTRR/Pyy3/YMK6hoaFQr5ycHJqZmVmn09YIX7ksWbKETZs2AbB69WrmzJnDtddeu0O9nRk5fOedd9i2bRtVVVW88847PProo1x11VWMHTuW1atXs2bNGvr378+sWbO45557SvZeWuLk0MzMzCyn2tpa9tprL0aPHs2oUaMYNmwYM2fO5Bvf+MYu9/nqq68yefJkALZs2cIXvvAFTjrpJABuueUWTjzxRLZu3cr555/PiBEjSvI+WuPk0MzMzCynuro6Fi9eTFVVVcn6PPjgg1myZEmz5yZOnMjEiRNLdq08/CgbMzMzsxwaGxvp1q1bSRPDjsjJoZmZmVkOVVVVrFq1qtJhlJ2TQzMzMzMrcHJoZmZmZgVODs3MzMyswMmhmZmZmRU4OTQzMzOzAieHZmZmZlbQZnKozJ9JuiodHyRpXPlDMzMzM7P2lmfk8AfAUcBZ6bgRuLVsEZmZmZlZxeRJDo+IiOnAJoCIeAPoVdaozMzMzHYj559/PgceeCCHHnro+8rnzZvH0KFDGTx4MN/5znfaLC+FPMnhZkndgQCQ1AfYVtIozMzMzHZj5513HvPmzXtf2datW5k+fTo//elPWb58Offeey/Lly9vsbxU8iSH3wMeBA6U9C3gV8C3SxaBmZmZWScyfvx4Vq5cCcDGjRt3GO3bFccddxz77bff+8oWLFjA4MGDOfjgg+nVqxdTp05lzpw5LZaXSo+2KkTE3ZIWAZ8FBHw+IlaULAIzMzOzTqS+vp4hQ4YAUFdXx8iRI3eoc+yxx9LY2LhD+Q033MDxxx+f6zrr1q1jwIABhePq6mrmz5/fYnmptJkcSrorIs4Bnm+mzMzMzKzd/ea663h3xfNtV9wJeww7hI99/eut1lm7di39+/enW7ds8rWuro5Ro0btUO/JJ5/8wPFExA5lklosL5U2k0NgRJOLdwcOL1kEZmZmZp1EbW3t+5LBRYsWMWXKlB3qlWLksLq6mpdffrlw3NDQQL9+/VosL5UWk0NJVwJfB/aS9BbZlDLAe8DtJYvAKuae+S8xp3Zdyfpbvv4thvftXbL+zMzMWtLWCF+5LFmyhE2bNgGwevVq5syZw7XXXrtDvVKMHI4dO5bVq1ezZs0a+vfvz6xZs7jnnnsYOnRos+Wl0uKClIj4dkRUAf8YEb0joir97B8RV5YsAquYObXrWL7+rZL1N7xvbyaN6V+y/szMzDqa2tpatm3bxujRo7nmmmsYNmwYM2fO/MD9nnXWWRx11FGsXLmS6upqZsyYQY8ePbjllls48cQTGTZsGGeeeSYjRoxosbxU8ixIuVLSR4AhwJ5F5b8sWRRWMcP79ua+i46qdBhmZmadQl1dHYsXL6aqqqqk/d57773Nlk+cOJGJEyfmLi+FPAtSLgQuAaqBWuBI4GngM2WJyMzMzKwDamxspFu3biVPDDuaPM85vAQYC6yNiE8DhwEbyhqVmZmZWQdTVVXFqlWrKh1G2eVJDjdFxCYASXtExPPA0DydSzpJ0kpJ9ZKuaOb8HpLuS+fnSxpYdO7KVL5S0ompbICkJyStkLRM0iVF9c9IZdsk1RSV95J0p6TnJC2RNL7o3OGpvF7S95TWgUvaT9Jjklan3x/J837NzMzMOrs8yWGDpA8D/wE8JmkO8EpbjdIjb24FTgaGA2dJGt6k2gXAGxExGPgucH1qOxyYSvYYnZOAH6T+tgCXR8Qwsunt6UV9LgVOBZreC/llgIgYCUwAbpS0/X3/EzCN7H7KIelaAFcAP4+IIcDP07GZmZlZl9dmchgRkyPidxFxNfANYAbw+Rx9jwPqI+KFiHgPmAVMalJnErB9ic/9wGfT6N0kYFZEvBsRa4B6YFxErI+IZ1NcjcAKoH86XhERK5uJYzhZgkdEvAb8DqiR1BfoHRFPR/Y0yX8tel/Fcc3M+X7NzMzMOr08I4dI+oikUUAj0ADk2USwP/By0XFDKmu2TkRsAd4E9s/TNk1BHwa0tV/MEmCSpB6SBpE9wHtA6q+hhWt8NCLWp7jWAwe2cQ0zMzOzLiHPauX/C5wHvABsS8VB26uVm9vHpel+Ly3VabWtpH2AHwOXRkRbD+q7AxgGLATWAr8mm57OE1+rJE0jm5bmoIMO2pmmZmZmZh1Snu3zzgQ+nqaGd0YD2QjddtXseK/i9joNknoA+wKvt9ZWUk+yxPDuiHigrSDSiORXtx9L+jWwGngj9dtcfK9K6hsR69P082st9H07abeYmpqanUoszczMzDqiPNPKS4EP70LfzwBDJA2S1ItsgcncJnXmAuem16cDj6f7/+YCU9Nq5kFki0UWpPsRZwArIuKmPEFI+pCkvdPrCcCWiFieposbJR2Z+v0iMKeZuM4tKjczMzPr0vKMHH4bWCxpKfDu9sKIOKW1RhGxRdLFwCNAd+COiFgm6RpgYUTMJUv07pJUTzZiODW1XSZpNrCcbAp4ekRslXQMcA7wnKTadKmvR8TDkiYD3wf6AA9Jqo2IE8nuF3xE0jZgXWq/3Z8DPwL2An6afgC+A8yWdAHwEnBGjs/JzMzMrNPLkxzOJHvEzHP84Z7DXCLiYeDhJmVXFb3eRAuJV0R8C/hWk7Jf0fy9gkTEg8CDzZS/SAvPZYyIhTSzuCYiNgKfba6NmZmZWSm9/PLLfPGLX+Q3v/kN3bp1Y9q0aVxySfYo53nz5nHJJZewdetWLrzwQq644opWy0shT3L424j4XsmuaGZmZmYFPXr04MYbb+STn/wkjY2NHH744UyYMIGhQ4cyffp0HnvsMaqrqxk7diynnHJKi+XDhzd9nPSuyXPP4SJJ35Z0lKRPbv8pydXNzMzMOpnx48ezcmX2aOWNGzdy6KF5nvDXsr59+/LJT2apVVVVFcOGDWPdunUsWLCAwYMHc/DBB9OrVy+mTp3KnDlzWiwvlTwjh4el30cWleV5lI2ZmZlZl1NfX8+QIUMAqKurY+TIkTvUOfbYY2lsbNyh/IYbbuD4449vse8XX3yRxYsXc8QRR/Doo48yYMAfHt5SXV3N/PnzWbduXbPlpdJmchgRny7Z1czMzMxK4MnZq/jty2+XtM8DBuzDsWd+otU6a9eupX///nTrlk2+1tXVMWrUqB3je/LJnb7+22+/zWmnncbNN99M7969yR7g8n6SWiwvlRaTQ0l/FhH/Jumy5s7nfZSMmZmZWVdRW1v7vmRw0aJFTJkyZYd6OztyuHnzZk477TTOPvtsTj31VCAbEXz55T9sGNfQ0EC/fv1aLC+V1kYO906/q5o55wc+m5mZWcW0NcJXLkuWLGHTpk0ArF69mjlz5nDttdfuUG9nRg4jggsuuIBhw4Zx2WV/GJMbO3Ysq1evZs2aNfTv359Zs2Zxzz33MHTo0GbLS6XF5DAibksvfxYRTxWfk3R0ySIwMzMz6yRqa2vZa6+9GD16NKNGjWLYsGHMnDmTb3zjG7vc51NPPcVdd93FyJEjGTNmDADXXXcdEydO5JZbbuHEE09k69atnH/++YwYMQKgxfJSyLMg5ftA09XJzZWZmZmZdWl1dXUsXryYqqrmJlZ3zTHHHNPsfYQAEydOZOLEibnLS6G1ew6PAv4Y6NPkvsPeZDuemJmZme02Ghsb6datW0kTw46otZHDXsA+qU7xp/AW2T7IZmZmZruNqqoqVq1aVekwyq61ew5/AfxC0o8iYi2ApG7APhHxVnsFaGZmZmbtJ88OKd+W1FvS3sByYKWkvylzXGZmZmZWAXmSw+FppPDzwMPAQcA5ZY3KzMzMzCoiT3LYU1JPsuRwTkRsxs85NDMzswpoaVWvNW9XPq88yeFtwItkD8X+paQ/IluUYmZmZtZu9txzTzZu3OgEMaeIYOPGjey555471S7P3srfA75XVLRWkvdbNjMzs3ZVXV1NQ0MDGzZsqHQoncaee+5JdXX1TrVpMzmU9FHgOqBfRJwsaThwFDBjl6I0MzMz2wU9e/Zk0KBBlQ6jy8szrfwj4BFg+47Oq4BLyxWQmZmZmVVOnuTwgIiYDWwDiIgtwNayRmVmZmZmFZEnOXxH0v6kFcqSjgTeLGtUZmZmZlYRbd5zCFwGzAU+LukpoA/ePs/MzMysS8qzWvlZSZ8ChgICVqZnHZqZmZlZF5Nn5HD7fYbLyhyLmZmZmVVYnnsOzczMzGw34eTQzMzMzAraTA4l3SXpy5IOaY+AzMzMzKxy8owc3gn0Bb4v6X8k/VjSJWWOy8zMzMwqIM9q5ccl/QIYC3wa+AowAvh/ZY7NzMzMzNpZnr2Vfw7sDTwNPAmMjYjXyh2YmZmZmbW/PNPKdcB7wKHAKOBQSXvl6VzSSZJWSqqXdEUz5/eQdF86P1/SwKJzV6bylZJOTGUDJD0haYWkZcXT25LOSGXbJNUUlfeUNFPSc6ndlal8qKTaop+3JF2azl0taV3RuYl53q+ZmZlZZ5dnWvmrAJL2Ab5Edg/ix4A9WmsnqTtwKzABaACekTQ3IpYXVbsAeCMiBkuaClwPTJE0HJhKNn3dD/iZpE8AW4DL04O5q4BFkh5LfS4FTgVuaxLKGcAeETFS0oeA5ZLujYiVwJiiWNcBDxa1+25E3NDW52NmZmbWleRZrXyxpPuAWuDzwB3AyTn6HgfUR8QLEfEeMAuY1KTOJGBmen0/8FlJSuWzIuLdiFgD1APjImJ9RDwLEBGNwAqgfzpekRK+pgLYW1IPYC+yUdC3mtT5LPA/EbE2x/syMzMz67Ly7JCyF3ATsCjtlJJXf+DlouMG4IiW6kTEFklvAvun8v9u0rZ/ccM0BX0YML+NOO4nSzbXAx8CvhoRrzepMxW4t0nZxZK+CCwkG618o43rmJmZmXV6bY4cRsQ/RsT8nUwMIduHeYfuctZptW2a4v4xcGlENB0FbGocsJVsenoQcLmkg4v66gWcAvx7UZt/Aj5ONu28HrixuY4lTZO0UNLCDRs2tBGGmZmZWcdXzh1SGoABRcfVwCst1UnTvvsCr7fWVlJPssTw7oh4IEccXwDmRcTmtMr6KaCm6PzJwLMR8er2goh4NSK2RsQ24IdkCeYOIuL2iKiJiJo+ffrkCMXMzMysYytncvgMMETSoDQ6NxWY26TOXODc9Pp04PGIiFQ+Na1mHgQMARak+xFnACsi4qaccbwEfEaZvYEjgeeLzp9FkyllSX2LDieTLXYxMzMz6/LKlhymaeiLgUfIFo7Mjohlkq6RdEqqNgPYX1I9cBlwRWq7DJgNLAfmAdMjYitwNHAOWbL3vsfMSJosqQE4CnhI0iPpGrcC+5AleM8Ad0ZEXWrzIbLV1E1HIP8hPfqmjuzB318t6YdjZmZm1kHleQj2qWSPmDmQ7F5AARERvdtqGxEPAw83Kbuq6PUmskfNNNf2W8C3mpT9iubvRyQiHuT9j6LZXv52K9f4PdkCmKbl5zRX38zMzKyry7Na+R+Az0XEinIHY2ZmZmaVlWda+VUnhmZmZma7hzwjhwvTQ7D/A3h3e2HOlcJmZmZm1onkSQ57A78HTigqC3ZcxGFmZmZmnVyevZW/1B6BmJmZmVnl5dlbuVrSg5Jek/SqpB9Lqm6P4MzMzMysfeVZkHIn2UOp+5Htb/yTVGZmZmZmXUye5LBPRNwZEVvSz48A7xVnZmZm1gXlSQ5/K+nPJHVPP38GbCx3YGZmZmbW/vIkh+cDZwK/AdaT7YF8fjmDMjMzM7PKyLNa+SXglLbqWfuYctvTJetr+fq3GN63zV0QzczMbDeSZ+TQOooobXfD+/Zm0pj+pe3UzMzMOrU8D8G2jkJw30VHVToKMzMz68JaHTmU1E3Sme0VjJmZmZlVVqvJYURsAy5up1jMzMzMrMLy3HP4mKS/ljRA0n7bf8oemZmZmZm1uzz3HG5/bM30orIADi59OGZmZmZWSXkeZTOoPQIxMzMzs8prMzmU1BP4c+C4VPRfwG0RsbmMcZmZmZlZBeSZVv4noCfwg3R8Tiq7sFxBmZmZmVll5EkOx0bE6KLjxyUtKZVDiBUAABinSURBVFdAZmZmZlY5eVYrb5X08e0Hkg4GtpYvJDMzMzOrlDwjh38DPCHpBUDAHwFfKmtUZmZmZlYRrSaHkroB/wsMAYaSJYfPR8S77RCbmZmZmbWzVpPDiNgm6caIOAqoa6eYzMzMzKxC8txz+Kik0ySp7NGYmZmZWUXluefwMmBvYIukTWRTyxERvcsamZmZmZm1u7buORQwIiJeaqd4zMzMzKyCWp1WjogAHmynWMzMzMyswvLcc/jfksbuSueSTpK0UlK9pCuaOb+HpPvS+fmSBhaduzKVr5R0YiobIOkJSSskLZN0SVH9M1LZNkk1ReU9Jc2U9Fxqd2XRuRdTea2khUXl+0l6TNLq9Psju/L+zczMzDqbPMnhp8kSxP+RVJeSqTZXLkvqDtwKnAwMB86SNLxJtQuANyJiMPBd4PrUdjgwFRgBnAT8IPW3Bbg8IoYBRwLTi/pcCpwK/LLJNc4A9oiIkcDhwEXFSSjw6YgYExE1RWVXAD+PiCHAz9OxmZmZWZeXZ0HKybvY9zigPiJeAJA0C5gELC+qMwm4Or2+H7gl3ec4CZiVnqe4RlI9MC4ingbWA0REo6QVQH9geUSsSNdpGkcAe0vqAewFvAe81Ubsk4Dx6fVM4L+Ar+V942ZmZmadVZsjhxGxFhgAfCa9/n2edmRJ28tFxw2prNk6EbEFeBPYP0/bNPp3GDC/jTjuB94hSypfAm6IiNfTuSB7VM8iSdOK2nw0IrYnoeuBA5vrWNI0SQslLdywYUMbYZiZmZl1fG0meZK+STZqtv1evZ7Av+Xou7nnIkbOOq22lbQP8GPg0ohoaxRwHNle0P2AQcDlaX9ogKMj4pNko6PTJR3XRl/vDyji9oioiYiaPn367ExTMzMzsw4pzwjgZOAUstE3IuIVoCpHuwayEcftqoFXWqqTpn33BV5vra2knmSJ4d0R8UCOOL4AzIuIzRHxGvAUUFP0XkjlD5IlkgCvSuqbrtcXeC3HdczMzMw6vTzJ4XvpkTYBIGnvnH0/AwyRNEhSL7IFJnOb1JkLnJtenw48nq41F5iaVjMPItvbeUG6H3EGsCIibsoZx0vAZ5TZm2why/OS9pZUVfSeTiBb1NI0rnOBOTmvZWZmZtap5UkOZ0u6DfiwpC8DPwN+2FajdA/hxcAjwApgdkQsk3SNpFNStRnA/mnByWWkVcERsQyYTbZ4ZR4wPSK2AkcD55Ale7XpZyKApMmSGoCjgIckPZKucSuwD1ni9wxwZ0TUAR8FfiVpCbAAeCgi5qU23wEmSFoNTEjHZmZmZl2esoG6NipJE8hG1gQ8EhGPlTuwzqampiYWLlzYdsVdNOPc2wC4YOZFZbuGmZmZ7R4kLWryGL+CPI+yISWDTgjNzMzMurg808pmZmZmtptwcmhmZmZmBbmSQ0l7SRpa7mDMzMzMrLLyPAT7c0At2aphJI2R1PSRNGZmZmbWBeQZObya7OHQvwOIiFpgYPlCMjMzM7NKyZMcbomIN8seiZmZmZlVXJ5H2SyV9AWgu6QhwF8Bvy5vWGZmZmZWCXlGDv8SGAG8C9wDvAlcWs6gzMzMzKwy8owcDo2IvwX+ttzBmJmZmVll5Rk5vEnS85L+r6QRZY/IzMzMzCqmzeQwIj4NjAc2ALdLek7S35U7MDMzMzNrf7kegh0Rv4mI7wFfIXvm4VVljcrMzMzMKiLPQ7CHSbpa0lLgFrKVytVlj8zMzMzM2l2eBSl3AvcCJ0TEK2WOx8zMzMwqqM3kMCKObI9AzMzMzKzyWkwOJc2OiDMlPQdE8SkgImJU2aMzMzMzs3bV2sjhJen3n7ZHIGZmZmZWeS0uSImI9enlX0TE2uIf4C/aJzwzMzMza095HmUzoZmyk0sdiJmZmZlVXmv3HP452QjhwZLqik5VAU+VOzAzMzMza3+t3XN4D/BT4NvAFUXljRHxelmjMjMzM7OKaDE5jIg3gTeBswAkHQjsCewjaZ+IeKl9QjQzMzOz9pJnh5TPSVoNrAF+AbxINqJoZmZmZl1MngUp1wJHAqsiYhDwWXzPoZmZmVmXlCc53BwRG4FukrpFxBPAmDLHZWZmZmYVkGdv5d9J2gf4JXC3pNeALeUNy8zMzMwqIc/I4STgf4GvAvOA/wE+V86gzMzMzKwy2kwOI+KdiNgaEVsiYmZEfC9NM7dJ0kmSVkqql3RFM+f3kHRfOj9f0sCic1em8pWSTkxlAyQ9IWmFpGWSLimqf0Yq2yappqi8p6SZkp5L7a7M0dfVktZJqk0/E/O8XzMzM7POrrWHYDcCUVyUjgVERPRurWNJ3YFbyXZYaQCekTQ3IpYXVbsAeCMiBkuaClwPTJE0HJgKjAD6AT+T9Amy6ezLI+JZSVXAIkmPpT6XAqcCtzUJ5Qxgj4gYKelDwHJJ9wLvttIXwHcj4obW3qOZmZlZV9Pa3spVEdG76Keq+HeOvscB9RHxQkS8B8wim6IuNgmYmV7fD3xWklL5rIh4NyLWAPXAuIhYHxHPpvgagRVA/3S8IiJWNvdWgL0l9QD2At4D3mqtLzMzM7PdVZ57DpF0jKQvpdcHSBqUo1l/4OWi4wZ2TL4KdSJiC9lDt/fP0zZNQR8GzG8jjvuBd4D1wEvADU13eGmhr4sl1Um6Q9JH2riGmZmZWZeQ5yHY3wS+BlyZinoB/5ajbzVTFjnrtNo2rZ7+MXBpRLzVRhzjgK1k09ODgMslHdxGX/8EfJzskT3rgRub61jSNEkLJS3csGFDG2GYmZmZdXx5Rg4nA6eQjb4REa8AVTnaNQADio6rgVdaqpOmffcFXm+traSeZMnc3RHxQI44vgDMi4jNEfEa2QO8a1rrKyJeTYtwtgE/JEswdxARt0dETUTU9OnTJ0coZmZmZh1bnuTwvYgI0sidpL1z9v0MMETSIEm9yBaYzG1SZy5wbnp9OvB4utZcYGpazTwIGAIsSPcjzgBWRMRNOeN4CfiMMnuT7fbyfGt9SepbdDiZbLGLmZmZWZeXJzmcLek24MOSvgz8DPiXthqlewgvBh4hW+wxOyKWSbpG0imp2gxgf0n1wGXAFantMmA2sJzs2YrTI2IrcDRwDlmy977HzEiaLKkBOAp4SNIj6Rq3AvuQJXjPAHdGRF1rfQH/kB59Uwd8muwZj2ZmZmZdnrKBujYqSROAE8juBXwkIh4rd2CdTU1NTSxcuLBs/c84N3tCzwUzLyrbNczMzGz3IGlRRNQ0dy7P9nmkZPCx1Fl3SWdHxN0ljNHMzMzMOoAWp5Ul9U67lNwi6YR0z97FwAvAme0XopmZmZm1l9ZGDu8C3gCeBi4E/obsMTaTIqK2HWIzMzMzs3bWWnJ4cESMBJD0L8BvgYPSbiJmZmZm1gW1tlp58/YXaaXwGieGZmZmZl1bayOHoyVt3zFEwF7pWEDk3F/ZzMzMzDqRFpPDiOjenoGYmZmZWeXleQi2mZmZme0mnByamZmZWYGTQzMzMzMrcHJoZmZmZgVODs3MzMyswMmhmZmZmRU4OTQzMzOzAieHZmZmZlbg5NDMzMzMCpwcmpmZmVmBk0MzMzMzK3ByaGZmZmYFTg7NzMzMrMDJoZmZmZkVODk0MzMzswInh2ZmZmZW4OTQzMzMzAqcHJqZmZlZgZNDMzMzMytwcmhmZmZmBU4OzczMzKzAyaGZmZmZFZQ1OZR0kqSVkuolXdHM+T0k3ZfOz5c0sOjclal8paQTU9kASU9IWiFpmaRLiuqfkcq2SaopKu8paaak51K7K9uKT9KgFM/qFF+v0n86ZmZmZh1P2ZJDSd2BW4GTgeHAWZKGN6l2AfBGRAwGvgtcn9oOB6YCI4CTgB+k/rYAl0fEMOBIYHpRn0uBU4FfNrnGGcAeETESOBy4SNLANuK7HvhuRAwB3khxmpmZmXV55Rw5HAfUR8QLEfEeMAuY1KTOJGBmen0/8FlJSuWzIuLdiFgD1APjImJ9RDwLEBGNwAqgfzpeERErm4kjgL0l9QD2At4D3mopvnT9z6R4SPF9/oN+GGZmZmadQTmTw/7Ay0XHDams2ToRsQV4E9g/T9s0BX0YML+NOO4H3gHWAy8BN0TE661cY3/gdymeluI2MzMz65J6lLFvNVMWOeu02lbSPsCPgUsj4q024hgHbAX6AR8BnpT0s129djFJ04BpAAcddFAbYZiZmZl1fOUcOWwABhQdVwOvtFQnTfvuC7zeWltJPckSw7sj4oEccXwBmBcRmyPiNeApoKaVa/wW+HCKp6W4AYiI2yOiJiJq+vTpkyMUMzMzs46tnMnhM8CQtPK3F9kCk7lN6swFzk2vTwcej4hI5VPTauZBwBBgQbofcAawIiJuyhnHS8BnlNmbbCHL8y3Fl67/RIqHFN+cnX73ZmZmZp1Q2ZLDdM/excAjZAtHZkfEMknXSDolVZsB7C+pHrgMuCK1XQbMBpYD84DpEbEVOBo4hyzZq00/EwEkTZbUABwFPCTpkXSNW4F9yFYzPwPcGRF1LcWX2nwNuCzFtX+K08zMzKzLUzZQZh9UTU1NLFy4sGz9zzj3NgAumHlR2a5hZmZmuwdJiyKiprlz3iHFzMzMzAqcHJqZmZlZgZNDMzMzMytwcmhmZmZmBU4OzczMzKzAyaGZmZmZFTg5NDMzM7MCJ4dmZmZmVuDk0MzMzMwKnByamZmZWYGTQzMzMzMrcHJoZmZmZgVODs3MzMyswMmhmZmZmRU4OTQzMzOzAieHZmZmZlbg5NDMzMzMCpwcmpmZmVmBk0MzMzMzK3ByaGZmZmYFPSodgOWjbq9WOgQzMzPbDTg57CTOv/OqSodgZmZmuwFPK5uZmZlZgZNDMzMzMytwcmhmZmZmBU4OzczMzKzAyaGZmZmZFTg5NDMzM7MCJ4dmZmZmVlDW5FDSSZJWSqqXdEUz5/eQdF86P1/SwKJzV6bylZJOTGUDJD0haYWkZZIuKap/RirbJqmmqPxsSbVFP9skjZFU1aT8t5JuTm3Ok7Sh6NyF5fyczMzMzDqKsj0EW1J34FZgAtAAPCNpbkQsL6p2AfBGRAyWNBW4HpgiaTgwFRgB9AN+JukTwBbg8oh4VlIVsEjSY6nPpcCpwG3FcUTE3cDdKaaRwJyIqE2nxxTFuwh4oKjpfRFxcUk+DDMzM7NOopwjh+OA+oh4ISLeA2YBk5rUmQTMTK/vBz4rSal8VkS8GxFrgHpgXESsj4hnASKiEVgB9E/HKyJiZRsxnQXc27RQ0hDgQODJXXifZmZmZl1GOZPD/sDLRccNqazZOhGxBXgT2D9P2zQFfRgwfydimkIzySFZ0nhfRERR2WmS6iTdL2nATlzDzMzMrNMq597KaqYsctZpta2kfYAfA5dGxFu5gpGOAH4fEUubOT0VOKfo+CfAvRHxrqSvkI1ufqaZPqcB09Lh25KajlweAPw2T3zWrvy9dDz+Tjomfy8dj7+Tjqezfid/1NKJciaHDUDxiFs18EoLdRok9QD2BV5vra2knmSJ4d0R8QD5TaX5KeXRQI+IWLS9LCI2FlX5Idm9kDuIiNuB21u6oKSFEVHT0nmrDH8vHY+/k47J30vH4++k4+mK30k5p5WfAYZIGiSpF1lyNrdJnbnAuen16cDjaWp3LjA1rWYeBAwBFqT7EWcAKyLipryBSOoGnEF232NTO9yHKKlv0eEpZPc2mpmZmXV5ZRs5jIgtki4GHgG6A3dExDJJ1wALI2IuWaJ3l6R6shHDqantMkmzgeVkK5SnR8RWSceQTf8+J2n7iuOvR8TDkiYD3wf6AA9Jqo2IE1Od44CGiHihmVDPBCY2KfsrSaeka78OnPfBPxEzMzOzjk/vX4NhpSRpWpp6tg7E30vH4++kY/L30vH4O+l4uuJ34uTQzMzMzAq8fZ6ZmZmZFTg5LJO2tg609ifpRUnPpS0RF1Y6nt2VpDskvSZpaVHZfpIek7Q6/f5IJWPc3bTwnVwtaV3RNqJN7822Mmppu1j/rVRWK99Ll/p78bRyGaStA1dRtHUgcFaTrQOtnUl6EaiJiM74PKouQ9JxwNvAv0bEoansH4DXI+I76f9MfSQivlbJOHcnLXwnVwNvR8QNlYxtd5WemtG3eLtY4PNkCyT9t1IhrXwvZ9KF/l48clgeebYONNstRcQvyZ4CUKx4K82ZZP+xtXbSwndiFdTKdrH+W6mg1rbx7UqcHJZHnq0Drf0F8KikRWl3G+s4PhoR6yH7jy/ZXudWeRenbUTv8PRl5TTZLtZ/Kx1EM9v4dpm/FyeH5ZFn60Brf0dHxCeBk4HpaSrNzJr3T8DHgTHAeuDGyoaze9qV7WKt/Jr5XrrU34uTw/LIs3WgtbOIeCX9fg14kGz63zqGV7fvTJR+v1bheHZ7EfFqRGyNiG1k24j676WdtbBdrP9WKqy576Wr/b04OSyPPFsHWjuStHe6eRhJewMnAEtbb2XtqHgrzXOBORWMxdhhG9HJ+O+lXbWyXaz/Viqope+lq/29eLVymaRl7Dfzh60Dv1XhkHZrkg4mGy2EbNvIe/ydVIake4HxwAHAq8A3gf8AZgMHAS8BZ0SEF0i0kxa+k/FkU2QBvAhctP1eNyu/tF3sk8BzwLZU/HWy+9v8t1IhrXwvZ9GF/l6cHJqZmZlZgaeVzczMzKzAyaGZmZmZFTg5NDMzM7MCJ4dmZmZmVuDk0MzMzMwKnByaWYclKSTdWHT815KuLlHfP5J0ein6auM6Z0haIemJEvR1qaQPfcA+viLpi+n1eZL6fdC4ivoeL+mPm7uWmXUeTg7NrCN7FzhV0gGVDqSYpO47Uf0C4C8i4tMluPSlwE4lh01jjYh/joh/TYfnATuVHErq0crp8UAhOWxyLTPrJJwcmllHtgW4Hfhq0xNNR/4kvZ1+j5f0C0mzJa2S9B1JZ0taIOk5SR8v6uZ4SU+men+a2neX9I+SnpFUJ+mion6fkHQP2QNwm8ZzVup/qaTrU9lVwDHAP0v6xyb1la6zNLWbUnSd/yyqd0sa4fsrskTuie2jkJJOkPS0pGcl/Xva7xVJL0q6StKvgDOaXPfqNAJ7OlAD3C2pVtJekg5Pn90iSY8UbdP2X5Kuk/QL4BJJn5M0X9JiST+T9FFJA4GvAF9N/R27/VqpjzGS/jt9pg9K+khR39en72eVpGNT+YhUVpvaDGnpH4mZlZaTQzPr6G4Fzpa07060GQ1cAowEzgE+ERHjgH8B/rKo3kDgU8CfkCVwe5KN9L0ZEWOBscCXJQ1K9ccBfxsRw4svlqZmrwc+Q7ZLwlhJn4+Ia4CFwNkR8TdNYjw11R0NHA/8Y5MtuN4nIr5Htkf7pyPi02k09e+A4yPik+k6lxU12RQRx0TErBb6u78otjFkifj3gdMj4nDgDqB4F6EPR8SnIuJG4FfAkRFxGDAL+D8R8SLwz8B3I2JMRDzZ5JL/CnwtIkaRJdffLDrXI30/lxaVfwX4fym2GrI9682sHbQ2PWBmVnER8ZakfwX+CvjfnM2e2b51laT/AR5N5c8BxdO7syNiG7Ba0gvAIWT7bo8qGpXcFxgCvAcsiIg1zVxvLPBfEbEhXfNu4DiybQFbcgxwb0RsBV5No3JjgbdyvscjgeHAU5IAegFPF52/L2c/2w0FDgUeS/11B4q3/yrurxq4LyWzvYDmPpOClNh/OCJ+kYpmAv9eVOWB9HsRWcIO2Xv5W0nVwAMRsXon34+Z7SInh2bWGdwMPAvcWVS2hTT7oSyb6VV07t2i19uKjrfx/v/uNd0/NAABfxkRjxSfkDQeeKeF+NTmO8jfpvC+kj1baf9YRJzVwvmWYm0tnmURcVSO/r4P3BQRc9PncvVOXqup7d/PVtL3ExH3SJpPNqr7iKQLI+LxD3gdM8vB08pm1uFFxOvAbLIp3+1eBA5PrycBPXeh6zMkdUv3IR4MrAQeAf5cUk8ASZ+QtHcb/cwHPiXpgLQA5CzgF220+SUwJd3j2IdspHEBsBYYLmmPNOL22aI2jUBVev3fwNGSBqc4PyTpEznfd3P9rQT6SDoq9ddT0ogW2u0LrEuvz22hv4KIeBN4Y/v9hGRT/a1+PpIOBl5I0+lzgVFtvx0zKwWPHJpZZ3EjcHHR8Q+BOZIWAD9n50fKIEuIfgF8FPhKRGyS9C9kU5vPphHJDcDnW+skItZLuhJ4gmwE7uGImNPGtR8EjgKWkI1Y/p+I+A2ApNlAHbAaWFzU5nbgp5LWp/sOzwPulbRHOv93wKpc7zzzI7J7Lf83xXI68L2UlPYgG7Fd1ky7q4F/l7SOLEndfk/mT4D7JU3i/fd2QpZE/rOyR/G8AHypjdimAH8maTPwG+CanXhfZvYBKKLprIqZmZmZ7a48rWxmZmZmBU4OzczMzKzAyaGZmZmZFTg5NDMzM7MCJ4dmZmZmVuDk0MzMzMwKnByamZmZWYGTQzMzMzMr+P9wozh7OZ72XwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10,5)) \n",
"plt.xlabel('Number of outer iterations')\n",
"plt.ylabel('Relative error on w estimate')\n",
"\n",
"norm_coef = np.linalg.norm(coef)\n",
"\n",
"for mu in [2,15,50,100,200]:\n",
" w_seq = w_seq_mu[mu]\n",
" y_axis = [np.linalg.norm(w-coef)/norm_coef for w in w_seq[1:]]\n",
" x_axis = np.arange(1,len(w_seq))\n",
" plt.step(x_axis, y_axis, label='$\\mu$ = '+str(mu))\n",
"\n",
"plt.legend(loc = 'upper right')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can note that $w$ converges towards a *biased* estimation of the initial parameter.\n",
"\n",
"And we can check that we minimize the initial lasso function:"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAApgAAAE9CAYAAACvCBIJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de7iWdZ3o//cHkMj2UkeMRlm6kYFwCcjynKMQGZZS6XhIIJtJLcvftrLdrymb2ZPl5Z6xxmZrY5eNiUVN4qhNP8oIdbaHzJ2ixGKFIIcUNssTRqWogwf4/P54buxhsQ4PeD/r+H5d13Ot5/7e38PnXjfL6+P3ex8iM5EkSZLKMqS3A5AkSdLAYoIpSZKkUplgSpIkqVQmmJIkSSqVCaYkSZJKZYIpSZKkUg3r7QD0R/vtt1+OGTOmt8OQJEnq1pIlS36bmW/taJ8JZh8yZswYHn744d4OQ5IkqVsRsb6zfS6RS5IkqVQmmJIkSSqVCaYkSZJK5TWYkiRp0Hj11Vdpa2tjy5YtvR1KvzFixAgaGxvZY489am5jgilJkgaNtrY2GhoaGDNmDBHR2+H0eZnJpk2baGtr4+CDD665nUvkkiRp0NiyZQsjR440uaxRRDBy5MhdnvE1wZQkSYOKyeWu2Z3flwmmJEmSSmWCKUmSpFKZYA4iN5x3GTecd1lvhyFJkkq2YcMG3vWud9HU1MTEiRO5+uqrezUe7yIfRHLb23o7BEmSVAfDhg3j61//OkcccQSbN2/myCOP5KSTTuLQQw/tlXicwZQkSeph06dPZ9WqVQBs2rSJSZMmvaH+9t9/f4444ggAGhoaaGpq4oknnnjDce4uZzAlSZJ62Nq1axk/fjwAra2tTJ48eac6U6dOZfPmzTuVX3nllcyYMaPTvtetW8fSpUs59thjywt4F5lgSpKkQekrP3mEFU8+X2qfhx6wF5d+YGKXddavX8/o0aMZMqSykNza2sphhx22U7377rtvl8d/4YUXOPPMM7nqqqvYa6+9drl9WUwwJUmSelBLS8sOCeWSJUuYNWvWTvV2dQbz1Vdf5cwzz+Scc87hjDPOKDfoXWSCKUmSBqXuZhrrZdmyZa+/GWfNmjUsWLCAyy+/fKd6uzKDmZl89KMfpampic9+9rOlxbq7vMlHkiSpB7W0tLBt2zamTJnCZZddRlNTE/PmzXtDfd5///18//vf56677qK5uZnm5mYWLlxYUsS7zhlMSZKkHtTa2srSpUtpaGgorc8TTjiBzCytvzfKGUxJkqQesnnzZoYMGVJqctkXmWBKkiT1kIaGBlavXt3bYdSdCaYkSZJKZYIpSZKkUplgSpIkqVQmmJIkSSqVjymqs4g4CLgG+C2wOjOv6OWQJEmS6qpuM5gRMSIiFkfEsoh4JCK+0km9fSLi1oh4NCJWRsRxVfuGRsTSiLitljYRcXFELC/G+0xRNiEiWqo+z2/ft5vHdUNEbIyI5e3KT46IVRGxNiIuqdr1duCnmXk+cOjujitJktRf1HOJ/GXgxMycAjQDJ0fEOzqodzWwKDMPAaYAK6v2Xdxuu9M2ETEJuAA4pih7f0SMz8xVmdmcmc3AkcBLwI/adxgRoyKioV3ZuA7G/i5wcrt6Q4FvAqdQSSLnRMT2ZHIpMDsi7gLu7qA/SZKkAaVuCWZWvFBs7lF8dnjEfETsBUwD5hZtXsnMPxT7GoH3AdfX2KYJeCAzX8rM14B7gdPbhfVu4DeZub6DkN8JLIiIEcU4FwDf6OC4fg78rl3xMcDazHwsM18BbgJOK/adB1yamScWxyNJkjSg1fUmn2KJuwXYCNyZmQ+2qzIWeBb4TrEUfn1EvKXYdxXweWBbjW2WA9MiYmRE7AnMBA5s13Y2ML+jWDPzFmARcFNEnAOcD5xd46GOBjZUbbcVZRR9fjoivgWs66hxRHwgIq577rnnahxOkiRpR+effz6jRo1i0qRJO5SPGTOGyZMn09zczFFHHdUjsdQ1wczMrcXSdCNwTLGMXW0YcARwbWYeDrwIXBIR7wc2ZuaSDrrtsE1mrgS+CtxJJalbBry2vVFEDAdOBW7pIt6vAVuAa4FTq2ZguxMddVf0uTwzz8rMCzPzc52M+5PM/Pjee+9d43CSJEk7Ovfcc1m0aFGH++6++25aWlp4+OGHeySWHnlMUbGEfQ/trl2kMtPXVjWzeSuV5PF44NSIWEdlufnEiPjXbtqQmXMz84jMnEZlGXtN1VinAL/KzGc6izMipgKTqFyjeekuHGIbO86WNgJP7kJ7SZI0iEyfPp1Vq1YBsGnTpp1mHXfHtGnT2Hfffd9wP2Wo513kb42IfYrvbwZmAI9W18nMp4ENETGhKHo3sCIzv5iZjZk5hsqy9l2Z+eGu2hTjjCp+HgScwY7L4XPoZHm8aHM48G0q106eB+wbEZfXeLgPAeMj4uBipnQ28OMa20qSpEFm7dq1jB8/HoDW1lYmT568U52pU6fS3Ny80+c//uM/dmmsiOA973kPRx55JNddd10p8Xenns/B3B+YV9xhPQS4OTNvA4iIhcDHMvNJ4FPAD4rE7DEqyV13Omvzw4gYCbwKXJSZvy/G2xM4CfhEF33uCXwwM39TtPkIcG77ShExH5gO7BcRbVRu4JkbEZ8EbgeGAjdk5iM1HIckSeotP7sEnv51uX3+6WQ4petHXq9fv57Ro0czZEhlnq+1tZXDDjtsp3r33XdfKSHdf//9HHDAAWzcuJGTTjqJQw45hGnTppXSd2fqlmBmZitweCf7ZlZ9bwE6veI0M++hsrxeXdZhm8yc2kkfLwEju4n3/nbbr1KZ0Wxfb04n7RcCC7saQ5IkqaWlZYeEcsmSJcyaNWunelOnTmXz5s07lV955ZXMmDGj5vEOOOAAAEaNGsXpp5/O4sWL+2+CKUmS1Kd1M9NYL8uWLWPLli0ArFmzhgULFnD55TtflVfGDOaLL77Itm3baGho4MUXX+SOO+7gS1/60hvutzu+i1ySJKkHtbS0sG3bNqZMmcJll11GU1MT8+bNe8P9zpkzh+OOO45Vq1bR2NjI3LlzeeaZZzjhhBOYMmUKxxxzDO973/s4+eT291yXzxlMSZKkHtTa2srSpUtpaGjovvIumD+/43uZly1bVuo4tXAGU5IkqYds3ryZIUOGlJ5c9jUmmJIkST2koaGB1atX93YYdWeCKUmSpFKZYEqSJKlUJpiSJEkqlQmmJEmSSmWCKUmSpFKZYEqSJKlUJpiSJEkqlQmmJEmSSmWCKUmSNACMGTOGyZMn09zczFFHHfV6+aJFi5gwYQLjxo3jiiuu6JFYfBe5JEnSAHH33Xez3377vb69detWLrroIu68804aGxs5+uijOfXUUzn00EPrGoczmJIkST1s+vTprFq1CoBNmzYxadKkuoyzePFixo0bx9ixYxk+fDizZ89mwYIFdRmrmjOYkiRJPWzt2rWMHz8egNbWViZPnrxTnalTp7J58+adyq+88kpmzJixU3lE8J73vIeI4BOf+AQf//jHeeKJJzjwwANfr9PY2MiDDz5Y4pF0zARTkiQNSl9d/FUe/d2jpfZ5yL6H8IVjvtBlnfXr1zN69GiGDKksJLe2tnLYYYftVO++++7bpbHvv/9+DjjgADZu3MhJJ53EIYccQmbuVC8idqnf3WGCKUmS1INaWlp2SCiXLFnCrFmzdqq3qzOYBxxwAACjRo3i9NNPZ/HixRx//PFs2LDh9TptbW2v16snE0xJkjQodTfTWC/Lli1jy5YtAKxZs4YFCxZw+eWX71RvV2YwX3zxRbZt20ZDQwMvvvgid9xxB1/60pc4+uijWbNmDY8//jijR4/mpptu4sYbbyztWDpjgilJktSDWlpaePOb38yUKVM47LDDaGpqYt68efzd3/3dbvf5zDPPcPrppwPw2muv8aEPfYiTTz4ZgGuuuYb3vve9bN26lfPPP5+JEyeWchxdMcGUJEnqQa2trSxdupSGhobS+hw7dizLli3rcN/MmTOZOXNmaWPVwscUSZIk9ZDNmzczZMiQUpPLvsgEU5IkqYc0NDSwevXq3g6j7kwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSBoDzzz+fUaNGMWnSpB3KFy1axIQJExg3bhxXXHFFt+VlMMGUJEkaAM4991wWLVq0Q9nWrVu56KKL+NnPfsaKFSuYP38+K1as6LS8LCaYkiRJPWz69OmsWrUKgE2bNu0067g7pk2bxr777rtD2eLFixk3bhxjx45l+PDhzJ49mwULFnRaXpZhpfUkSZKkmqxdu5bx48cD0NrayuTJk3eqM3XqVDZv3rxT+ZVXXsmMGTNqGueJJ57gwAMPfH27sbGRBx98sNPysphgSpKkQenpv/97Xl75aKl9vqnpEP70b/6myzrr169n9OjRDBlSWUhubW3lsMMO26nefffd94bjycydyiKi0/KymGBKkiT1oJaWlh0SyiVLljBr1qyd6pUxg9nY2MiGDRte325ra+OAAw7otLwsJpiSJGlQ6m6msV6WLVvGli1bAFizZg0LFizg8ssv36leGTOYRx99NGvWrOHxxx9n9OjR3HTTTdx4441MmDChw/KyeJOPJElSD2ppaWHbtm1MmTKFyy67jKamJubNm/eG+50zZw7HHXccq1atorGxkblz5zJs2DCuueYa3vve99LU1MTZZ5/NxIkTOy0vizOYdRYRBwHXAL8FVmdmuQ+akiRJ/UpraytLly6loaGh1H7nz5/fYfnMmTOZOXNmzeVlqNsMZkSMiIjFEbEsIh6JiK90Um+fiLg1Ih6NiJURcVzVvqERsTQibqulTURcHBHLi/E+U5RNiIiWqs/z2/ft5nHdEBEbI2J5u/KTI2JVRKyNiEuqdr0d+Glmng8curvjSpKk/m/z5s0MGTKk9OSyr6nnEvnLwImZOQVoBk6OiHd0UO9qYFFmHgJMAVZW7bu43XanbSJiEnABcExR9v6IGJ+ZqzKzOTObgSOBl4Afte8wIkZFREO7snEdjP1d4OR29YYC3wROoZJEzomI7cnkUmB2RNwF3N1Bf5IkaZBoaGhg9erVvR1G3dUtwcyKF4rNPYrPDvfER8RewDRgbtHmlcz8Q7GvEXgfcH2NbZqABzLzpcx8DbgXOL1dWO8GfpOZ6zsI+Z3AgogYUYxzAfCNDo7r58Dv2hUfA6zNzMcy8xXgJuC0Yt95wKWZeWJxPJIkSQNaXW/yKZa4W4CNwJ2Z2f4JnmOBZ4HvFEvh10fEW4p9VwGfB7bV2GY5MC0iRkbEnsBM4MB2bWcDHV6gkJm3AIuAmyLiHOB84OwaD3U0sKFqu60oo+jz0xHxLWBdR40j4gMRcd1zzz1X43CSJEl9V10TzMzcWixNNwLHFMvY1YYBRwDXZubhwIvAJRHxfmBjZi7poNsO22TmSuCrwJ1UkrplwGvbG0XEcOBU4JYu4v0asAW4Fji1aga2Ox09mTSLPpdn5lmZeWFmfq6TcX+SmR/fe++9axxOkiSp7+qRxxQVS9j30O7aRSozfW1VM5u3UkkejwdOjYh1VJabT4yIf+2mDZk5NzOPyMxpVJax11SNdQrwq8x8prM4I2IqMInKNZqX7sIhtrHjbGkj8OQutJckSRow6nkX+VsjYp/i+5uBGcAO72PKzKeBDRExoSh6N7AiM7+YmY2ZOYbKsvZdmfnhrtoU44wqfh4EnMGOy+Fz6GR5vGhzOPBtKtdOngfsGxE7P/W0Yw8B4yPi4GKmdDbw4xrbSpIkDSj1fA7m/sC84g7rIcDNmXkbQEQsBD6WmU8CnwJ+UCRmj1FJ7rrTWZsfRsRI4FXgosz8fTHensBJwCe66HNP4IOZ+ZuizUeAc9tXioj5wHRgv4hoo3IDz9yI+CRwOzAUuCEzH6nhOCRJkgacuiWYmdkKHN7JvplV31uAo7ro5x4qy+vVZR22ycypnfTxEjCym3jvb7f9KpUZzfb15nTSfiGwsKsxJEmSBgNfFSlJkqRSmWBKkiT1cxs2bOBd73oXTU1NTJw4kauvvvr1fYsWLWLChAmMGzeOK664otvyMphgSpIk9XPDhg3j61//OitXruSBBx7gm9/8JitWrGDr1q1cdNFF/OxnP2PFihXMnz+/y/KymGBKkiT1sOnTp7Nq1SoANm3axKRJ7R8Vvmv2339/jjjiCKDyOsqmpiaeeOIJFi9ezLhx4xg7dizDhw9n9uzZLFiwoNPystTzLnJJkiR1YO3atYwfPx6A1tZWJk+evFOdqVOnsnnz5p3Kr7zySmbMmNFp3+vWrWPp0qUce+yx3HHHHRx44B8f1d3Y2MiDDz7IE0880WF5WUwwJUnSoHTfzav57YZaX9pXm/0O/C9MPfvtXdZZv349o0ePZsiQykJya2srhx122M7x3XffLo//wgsvcOaZZ3LVVVex1157kZk71YmITsvLYoIpSZLUg1paWnZIKJcsWcKsWbN2qrerM5ivvvoqZ555Jueccw5nnHEGUJmZ3LBhw+t12traOOCAAzotL4sJpiRJGpS6m2msl2XLlrFlyxYA1qxZw4IFC7j88p1fHrgrM5iZyUc/+lGampr47Gc/+3r50UcfzZo1a3j88ccZPXo0N910EzfeeCMTJkzosLws3uQjSZLUg1paWti2bRtTpkzhsssuo6mpiXnz5r2hPu+//36+//3vc9ddd9Hc3ExzczMLFy5k2LBhXHPNNbz3ve+lqamJs88+m4kTJ3ZaXhZnMCVJknpQa2srS5cupaGhobQ+TzjhhA6vqwSYOXMmM2fOrLm8DM5gSpIk9ZDNmzczZMiQUpPLvsgEU5IkqYc0NDSwevXq3g6j7kwwJUmSVCoTTEmSJJXKBFOSJEmlquku8oj4c2BMdf3M/F6dYpIkSaqbzCz1rTUDXWd3p3el2wQzIr4P/BnQAmzdPhZggilJkvqVESNGsGnTJkaOHGmSWYPMZNOmTYwYMWKX2tUyg3kUcGjuTvoqSZLUhzQ2NtLW1sazzz7b26H0GyNGjKCxsXGX2tSSYC4H/hR4aneCkiRJ6iv22GMPDj744N4OY8CrJcHcD1gREYuBl7cXZuapdYtKkiRJ/VYtCeaX6x2EJEmSBo5uE8zMvDci3gYcXRQtzsyN9Q1LkiRJ/VW3z8GMiLOBxcAHgbOBByPirHoHJkmSpP6pliXyvwWO3j5rGRFvBf4DuLWegUmSJKl/quVNPkPaLYlvqrGdJEmSBqFaZjAXRcTtwPxiexawsH4hSZIkqT+r5Safv46IM4HjgQCuy8wf1T0ySZIk9Us1vYs8M38I/LDOsUiSJGkA6DTBjIhfZOYJEbGZyrvHX98FZGbuVffoJEmS1O90mmBm5gnFz4aeC0eSJEn9XS3Pwfx+LWWSJEkS1Pa4oYnVGxExDDiyPuFIkiSpv+s0wYyILxbXXx4WEc8Xn83AM8CCHotQkiRJ/UqnCWZm/kNx/eU/ZuZexachM0dm5hd7MEZJkiT1I7UskS+OiL23b0TEPhHxF3WMSZIkSf1YLQnmpZn53PaNzPwDcGn9QpIkSVJ/VtO7yDsoq+kB7ZIkSRp8akkwH46If4qIP4uIsRHxv4Al9Q5MkiRJ/VMtCeangFeAfwNuAbYAF9UzKEmSJPVf3S51Z+aLwCU9EIskSZIGgG4TzIh4O/A5YEx1/cw8sX5hSZIkqb+q5WadW4BvAdcDW+sbzsATEQcB1wC/BVZn5hW9HJIkSVJd1XIN5muZeW1mLs7MJds/3TWKiBERsTgilkXEIxHxlU7q7RMRt0bEoxGxMiKOq9o3NCKWRsRttbSJiIsjYnkx3meKsgkR0VL1eX77vt0RETdExMaIWN6u/OSIWBURayOi+pKCtwM/zczzgUN3d1xJkqT+opYE8ycR8d8iYv+I2Hf7p4Z2LwMnZuYUoBk4OSLe0UG9q4FFmXkIMAVYWbXv4nbbnbaJiEnABcAxRdn7I2J8Zq7KzObMbKbyDvWXgB+17zAiRkVEQ7uycR2M/V3g5Hb1hgLfBE6hkkTOiYjtyeRSYHZE3AXc3UF/kiRJA0otCeZHgL8G/g+VxxMtAR7urlFWvFBs7lF8srpOROwFTAPmFm1eKR7kTkQ0Au+jsjRfS5sm4IHMfCkzXwPuBU5vF9a7gd9k5voOQn4nsCAiRhTjXAB8o4Pj+jnwu3bFxwBrM/OxzHwFuAk4rdh3HpWH1Z9YHI8kSdKA1m2CmZkHd/AZW0vnxRJ3C7ARuDMzH2xXZSzwLPCdYin8+oh4S7HvKuDzwLYa2ywHpkXEyIjYE5gJHNiu7WxgfifHeQuwCLgpIs4BzgfOruU4gdHAhqrttqKMos9PR8S3gHUdNY6ID0TEdc8991xHuyVJkvqVbhPMiPirjj61dJ6ZW4ul6UbgmGIZu9ow4Ajg2sw8HHgRuCQi3g9s7ORazw7bZOZK4KvAnVSSumXAa1XHMRw4lcpNS53F+zUqz/m8Fji1aga2O9FRd0WfyzPzrMy8MDM/18m4P8nMj++9994d7ZYkSepXalkiP7rqMxX4MpVErWbFEvY9tLt2kcpMX1vVzOatVJLH44FTI2IdleXmEyPiX7tpQ2bOzcwjMnMalWXsNVVjnQL8KjOf6SzOiJgKTKJyjeauvG+9jR1nSxuBJ3ehvSRJ0oBRyxL5p6o+FwCHA8O7axcRb42IfYrvbwZmAI+26/tpYENETCiK3g2syMwvZmZjZo6hsqx9V2Z+uKs2xTijip8HAWew43L4HDpZHi/aHA58m8q1k+cB+0bE5d0dZ+EhYHxEHFzMlM4GflxjW0mSpAGlludgtvcSML6GevsD84o7rIcAN2fmbQARsRD4WGY+SeVVlD8oErPHqCR33emszQ8jYiTwKnBRZv6+GG9P4CTgE130uSfwwcz8TdHmI8C57StFxHxgOrBfRLRRuYFnbkR8ErgdGArckJmP1HAcPS9h1r/8stQuT2sezYeOPajUPiVJUv9Vy5t8fsIf7/4eQuUxPDd31y4zW6nMdna0b2bV9xbgqC76uYfK8np1WYdtMnNqJ328BIzsJt77222/SmVGs329OZ20Xwgs7GqMPqGjq0XfgBVPPQ9ggilJkl7XaYIZEW/KzJeBK6uKXwPWZ2Zb3SNT3fzbJ47rvlKNyp4NlSRJ/V9XM5i/pHLzzMcy8y97KB5JkiT1c10lmMOL6xD/PCLOaL8zM/+9fmFJkiSpv+oqwbwQOAfYB/hAu30JmGBKkiRpJ50mmJn5C+AXEfFwZs7twZgkSZLUj9XyHEyTS0mSJNWsljf5SJIkSTUzwZQkSVKpanqTT0ScCkwrNu/NzJ/ULyRJkiT1Z93OYEbEPwAXU3nf9wrg00WZJEmStJNaZjDfBzRn5jaAiJgHLAW+WM/AJEmS1D/Veg3mPlXf965HIJIkSRoYapnB/AdgaUTcDQSVazGdvZQkSVKHuk0wM3N+RNwDHE0lwfxCZj5d78AkSZLUP9Vyk8/xwPOZ+WOgAfh8RPzXukcmSZKkfqmWazCvBV6KiCnAXwPrge/VNSpJkiT1W7UkmK9lZgKnAd/IzKupzGRKkiRJO6nlJp/NEfFF4MPAtIgYCuxR37AkSZLUX9UygzkLeBn4aHFzz2jgH+salSRJkvqtmmYwgaszc2tEvB04BJhf37AkSZLUX9Uyg/lz4E0RMRr438B5wHfrGZQkSZL6r1oSzMjMl4AzgH/OzNOBifUNS5IkSf1VTQlmRBwHnAP8tCgbWr+QJEmS1J/VkmBeTOXVkD/KzEciYixwd33DkiRJUn9Vy6sif07lOszt248Bn65nUJIkSeq/uk0wI+KtwOepXHc5Ynt5Zp5Yx7gkSZLUT9WyRP4D4FHgYOArwDrgoTrGJEmSpH6slgRzZGbOBV7NzHsz83zgHXWOS5IkSf1ULQ9af7X4+VREvA94EmisX0iSJEnqz2pJMC+PiL2B/xf4Z2Av4DN1jUqSJEn9Vi13kd9WfH0OeBdARJhgSpIkqUO1XIPZkc+WGoUkSZIGjN1NMKPUKCRJkjRg7G6CmaVGIUmSpAGj02swI2IzHSeSAby5bhFJkiSpX+s0wczMhp4MRJIkSQPD7i6RS5IkSR0ywZQkSVKpTDAlSZJUKhNMSZIklaqWV0XqDYiIg4BrgN8CqzPzil4OSZIkqa7qNoMZESMiYnFELIuIRyLiK53U2ycibo2IRyNiZUQcV7VvaEQsjYjbamkTERdHxPJivM8UZRMioqXq8/wbedVlRNwQERsjYnm78pMjYlVErI2IS6p2vR34aWaeDxy6u+NKkiT1F/VcIn8ZODEzpwDNwMkR8Y4O6l0NLMrMQ4ApwMqqfRe32+60TURMAi4AjinK3h8R4zNzVWY2Z2YzcCTwEvCj9h1GxKiIaGhXNq6Dsb8LnNyu3lDgm8ApVJLIORGxPZlcCsyOiLuAuzvoT5IkaUCpW4KZFS8Um3sUnx0e3B4RewHTgLlFm1cy8w/FvkbgfcD1NbZpAh7IzJcy8zXgXuD0dmG9G/hNZq7vIOR3AgsiYkQxzgXANzo4rp8Dv2tXfAywNjMfy8xXgJuA04p95wGXZuaJxfFIkiQNaHW9yadY4m4BNgJ3ZuaD7aqMBZ4FvlMshV8fEW8p9l0FfB7YVmOb5cC0iBgZEXsCM4ED27WdDczvKNbMvAVYBNwUEecA5wNn13ioo4ENVdttRRlFn5+OiG8B6zpqHBEfiIjrnnvuuRqHkyRJ6rvqmmBm5tZiaboROKZYxq42DDgCuDYzDwdeBC6JiPcDGzNzSQfddtgmM1cCXwXupJLULQNe294oIoYDpwK3dBHv14AtwLXAqVUzsN2Jjror+lyemWdl5oWZ+blOxv1JZn587733rnE4SZKkvqtHHlNULGHfQ7trF6nM9LVVzWzeSiV5PB44NSLWUVluPjEi/rWbNmTm3Mw8IjOnUVnGXlM11inArzLzmc7ijIipwCQq12heuguH2MaOs6WNwJO70F6SJGnAqOdd5G+NiH2K728GZgCPVtfJzKeBDRExoSh6N7AiM7+YmY2ZOYbKsvZdmfnhrtoU44wqfh4EnMGOy+Fz6AGSgPsAABDkSURBVGR5vGhzOPBtKtdOngfsGxGX13i4DwHjI+LgYqZ0NvDjGttKkiQNKPV8Dub+wLziDushwM2ZeRtARCwEPpaZTwKfAn5QJGaPUUnuutNZmx9GxEjgVeCizPx9Md6ewEnAJ7roc0/gg5n5m6LNR4Bz21eKiPnAdGC/iGijcgPP3Ij4JHA7MBS4ITMfqeE4JEmSBpy6JZiZ2Qoc3sm+mVXfW4CjuujnHirL69VlHbbJzKmd9PESMLKbeO9vt/0qlRnN9vXmdNJ+IbCwqzEkSZIGA18VKUmSpFKZYEqSJKlUJpiSJEkqlQmmJEmSSmWCKUmSpFKZYEqSJKlUJpiSJEkqlQmmJEmSSmWCKUmSpFLV81WRGiRWPPU8s/7ll6X1d1rzaD507EGl9SdJknqWCabekNOaR5fa34qnngcwwZQkqR8zwdQb8qFjDyo1GSxzJlSSJPUOr8GUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVCoTTEmSJJXKBFOSJEmlGtbbAQx0EXEQcA3wW2B1Zl7RyyH1eSueep5Z//LL0vo7rXk0Hzr2oNL6kyRJXavbDGZEjIiIxRGxLCIeiYivdFJvn4i4NSIejYiVEXFc1b6hEbE0Im6rpU1EXBwRy4vxPlPLGLtxXDdExMaIWN6u/OSIWBURayPikqpdbwd+mpnnA4fu7riDxWnNozl0/71K62/FU8+zoOWJ0vqTJEndq+cM5svAiZn5QkTsAfwiIn6WmQ+0q3c1sCgzz4qI4cCeVfsuBlYC7TOOndpExCTgAuAY4BVgUUT8NDPXdDMGABExCvjPzNxcVTYuM9e2q/pdKjOS36uqNxT4JnAS0AY8FBE/zswVwFLgbyNiFvD9rn9l+tCxB5U621jmTKgkSapN3WYws+KFYnOP4pPVdSJiL2AaMLdo80pm/qHY1wi8D7i+xjZNwAOZ+VJmvgbcC5ze1RjtvBNYEBEjinEuAL7RwXH9HPhdu+JjgLWZ+VhmvgLcBJxW7DsPuDQzTyyOR5IkaUCr600+xRJ3C7ARuDMzH2xXZSzwLPCdYin8+oh4S7HvKuDzwLYa2ywHpkXEyIjYE5gJHNjNGK/LzFuARcBNEXEOcD5wdo2HOhrYULXdVpRR9PnpiPgWsK6jxhHxgYi47rnnnqtxOEmSpL6rrglmZm7NzGagETimWMauNgw4Arg2Mw8HXgQuiYj3Axszc0kH3XbYJjNXAl8F7qSS1C0DXuusfifxfg3YAlwLnFo1A9ud6Ki7os/lmXlWZl6YmZ/rZNyfZObH99577xqHkyRJ6rt65DFFxZL0PcDJ7Xa1AW1VM5u3UkkGjwdOjYh1VJabT4yIf+2mDZk5NzOPyMxpVJax13RVv72ImApMAn4EXLoLh9hGZbZ0u0bgyV1oL0mSNGDU8y7yt0bEPsX3NwMzgEer62Tm08CGiJhQFL0bWJGZX8zMxswcA8wG7srMD3fVphhnVPHzIOAMYH5X9dvFezjwbSrXTp4H7BsRl9d4uA8B4yPi4OImotnAj2tsK0mSNKDU8y7y/YF5xR3WQ4CbM/M2gIhYCHwsM58EPgX8oEjMHqOS3HWnszY/jIiRwKvARZn5+27qV9sT+GBm/qaI8SPAue0rRcR8YDqwX0S0UbmBZ25EfBK4HRgK3JCZj9RwHJIkSQNO3RLMzGwFDu9k38yq7y3AUV30cw+V5fXqsg7bZObUTvrocoyizv3ttl+lMqPZvt6cTtovBBZ2NYYkSdJg4KsiJUmSVCoTTEmSJJXKBFOSJEmlMsGUJElSqUwwJUmSVKp6PqZIfVLCd0p+Jfrks+CoWp4uJUmSBgMTzMEkoniBZYme/nXlZx9OMFc89Tyz/uWXpfZ5WvNoPnTsQaX2KUnSQGGCOagMqbw1/byfltdl2bOhJTuteXTpfa546nkAE0xJkjphgqkB7UPHHlR6Ilj2bKgkSQONN/lIkiSpVCaYkiRJKpUJpiRJkkplgilJkqRSmWBKkiSpVN5FLu2Gsp+t6XM1JUkDiQmmtIvKframz9WUJA00JpiDzLbcynmLSnzrTjzDzHwLHyyvxz6v7Gdr+lxNSdJA4zWYg8geQ4cxJIaW2ucqXmFhvFhqn5IkqX9zBnMQ2WPIcPYYMpzvnPyd0vo877tHldbXYOY1nZKkgcQEc5DZ9tJLrP/Lvyqtv9lP/ycrJ/rP6I3wmk5J0kATmdnbMahw1FFH5cMPP1y3/m+55HZ+9/tkr22/K63P1154nje9/DBnf+Cp0vpk8llwVInXiQ4ys/7ll6x46nkO3X+v3g6lS86ySlL/FhFLMrPDpUwTzD6k3gnmI/c9werFz5Ta55Nr/gDAiP9cU06H2/89RpTTX51sG/oE/2fawtL6O4ThfCH/pJS+ntm8hd++8HIpfW23ftifMW/vC0vr78HHK/+Tc+zB+5bWpyTpjw49YC8u/cDEuo7RVYLp2uYgMnHqaCZOLXc59tav/DPPPf5mKOvmodj2xyRTu+VtDSN4W8OIUvuc+Kf7M/OU40rr78YH/y8LWp4orT9JUt/iDGYfUu8ZTEmSpLJ0NYPpY4okSZJUKhNMSZIklcoEU5IkSaUywZQkSVKpTDAlSZJUKhNMSZIklcoEU5IkSaUywZQkSVKpTDAlSZJUKhNMSZIklcoEU5IkSaXyXeR9SEQ8C6xvV7wf8NteCEdd87z0PZ6Tvsdz0jd5Xvqe/npO/mtmvrWjHSaYfVxEPNzZi+TVezwvfY/npO/xnPRNnpe+ZyCeE5fIJUmSVCoTTEmSJJXKBLPvu663A1CHPC99j+ek7/Gc9E2el75nwJ0Tr8GUJElSqZzBlCRJUqlMMPuwiDg5IlZFxNqIuKS34xFExLqI+HVEtETEw70dz2AVETdExMaIWF5Vtm9E3BkRa4qff9KbMQ42nZyTL0fEE8XfS0tEzOzNGAebiDgwIu6OiJUR8UhEXFyU+7fSi7o4LwPq78Ul8j4qIoYCq4GTgDbgIWBOZq7o1cAGuYhYBxyVmf3xeWUDRkRMA14AvpeZk4qyrwG/y8wriv8h+5PM/EJvxjmYdHJOvgy8kJlX9mZsg1VE7A/sn5m/iogGYAnwF8C5+LfSa7o4L2czgP5enMHsu44B1mbmY5n5CnATcFovxyT1CZn5c+B37YpPA+YV3+dR+Q+2ekgn50S9KDOfysxfFd83AyuB0fi30qu6OC8Diglm3zUa2FC13cYA/AfYDyVwR0QsiYiP93Yw2sHbMvMpqPwHHBjVy/Go4pMR0VosobsU20siYgxwOPAg/q30Ge3OCwygvxcTzL4rOijzeobed3xmHgGcAlxULAtK6ti1wJ8BzcBTwNd7N5zBKSL+C/BD4DOZ+Xxvx6OKDs7LgPp7McHsu9qAA6u2G4EneykWFTLzyeLnRuBHVC5lUN/wTHFt0/ZrnDb2cjyDXmY+k5lbM3Mb8G38e+lxEbEHlSTmB5n570Wxfyu9rKPzMtD+Xkww+66HgPERcXBEDAdmAz/u5ZgGtYh4S3FBNhHxFuA9wPKuW6kH/Rj4SPH9I8CCXoxFvJ68bHc6/r30qIgIYC6wMjP/qWqXfyu9qLPzMtD+XryLvA8rHlFwFTAUuCEz/2cvhzSoRcRYKrOWAMOAGz0nvSMi5gPTgf2AZ4BLgf8PuBk4CPi/wAcz05tOekgn52Q6leW+BNYBn9h+7Z/qLyJOAO4Dfg1sK4r/hsr1fv6t9JIuzsscBtDfiwmmJEmSSuUSuSRJkkplgilJkqRSmWBKkiSpVCaYkiRJKpUJpiRJkkplgilpQIuIjIivV21/LiK+XFLf342Is8roq5txPhgRKyPi7hL6+kxE7PkG+7gwIv6q+H5uRBzwRuOq6nt6RPx5R2NJ6j9MMCUNdC8DZ0TEfr0dSLWIGLoL1T8K/LfMfFcJQ38G2KUEs32smfmtzPxesXkusEsJZkQM62L3dOD1BLPdWJL6CRNMSQPda8B1wH9vv6P9DGREvFD8nB4R90bEzRGxOiKuiIhzImJxRPw6Iv6sqpsZEXFfUe/9RfuhEfGPEfFQRLRGxCeq+r07Im6k8pDl9vHMKfpfHhFfLcq+BJwAfCsi/rFd/SjGWV60m1U1zm1V9a4pZho/TSUZvHv7bGhEvCcifhkRv4qIW4r3IxMR6yLiSxHxC+CD7cb9cjETfBZwFPCDiGiJiDdHxJHF725JRNxe9UrCeyLi7yPiXuDiiPhARDwYEUsj4j8i4m0RMQa4EPjvRX9Tt49V9NEcEQ8Uv9MfRcSfVPX91eL8rI6IqUX5xKKspWgzvrN/JJLKZYIpaTD4JnBOROy9C22mABcDk4G/BN6emccA1wOfqqo3Bngn8D4qSeAIKjOOz2Xm0cDRwAURcXBR/xjgbzPz0OrBimXmrwInUnmbx9ER8ReZeRnwMHBOZv51uxjPKOpOAWYA/9judXM7yMxvAE8C78rMdxWzuv8DmJGZRxTjfLaqyZbMPCEzb+qkv1urYmumksz/M3BWZh4J3ABUv+1qn8x8Z2Z+HfgF8I7MPBy4Cfh8Zq4DvgX8r8xszsz72g35PeALmXkYlQT90qp9w4rz85mq8guBq4vYjgLaOvvdSCpXV8sUkjQgZObzEfE94NPAf9bY7KHtr2mLiN8AdxTlvwaql6pvzsxtwJqIeAw4hMp76g+rmh3dGxgPvAIszszHOxjvaOCezHy2GPMHwDQqr8DszAnA/MzcCjxTzA4eDTxf4zG+AzgUuD8iAIYDv6za/2819rPdBGAScGfR31Cg+lV31f01Av9WJMTDgY5+J68r/udgn8y8tyiaB9xSVeXfi59LqCT9UDmWv42IRuDfM3PNLh6PpN1kgilpsLgK+BXwnaqy1yhWcqKSEQ2v2vdy1fdtVdvb2PG/ne3ft5tAAJ/KzNurd0TEdODFTuKLbo+g9javH1dhRBft78zMOZ3s7yzWruJ5JDOPq6G/fwb+KTN/XPxevryLY7W3/fxspTg/mXljRDxIZXb59oj4WGbe9QbHkVQDl8glDQqZ+TvgZirL19utA44svp8G7LEbXX8wIoYU12WOBVYBtwP/T0TsARARb4+It3TTz4PAOyNiv+KmmjnAvd20+Tkwq7jm861UZjwXA+uBQyPiTcXM37ur2mwGGorvDwDHR8S4Is49I+LtNR53R/2tAt4aEccV/e0RERM7abc38ETx/SOd9Pe6zHwO+P326yupXLbQ5e8nIsYCjxWXBvwYOKz7w5FUBmcwJQ0mXwc+WbX9bWBBRCwG/je7PmMHlaTqXuBtwIWZuSUirqeyTPurYmb0WeAvuuokM5+KiC8Cd1OZCVyYmQu6GftHwHHAMiozp5/PzKcBIuJmoBVYAyytanMd8LOIeKq4DvNcYH5EvKnY/z+A1TUdecV3qVx7+p9FLGcB3ygS22FUZo4f6aDdl4FbIuIJKonu9mtUfwLcGhGnseO1rlBJRL8VlccsPQac101ss4APR8SrwNPAZbtwXJLegMhsv7ojSZIk7T6XyCVJklQqE0xJkiSVygRTkiRJpTLBlCRJUqlMMCVJklQqE0xJkiSVygRTkiRJpTLBlCRJUqn+fycogck/HdktAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(10,5)) \n",
"plt.xlabel('Number of outer iterations')\n",
"plt.ylabel('Lasso function')\n",
"\n",
"norm_coef = np.linalg.norm(coef)\n",
"\n",
"for mu in [2,15,50,100,200]:\n",
" w_seq = w_seq_mu[mu]\n",
" y_axis = [get_lasso(w) for w in w_seq[1:]]\n",
" x_axis = np.arange(1,len(w_seq))\n",
" plt.step(x_axis, y_axis, label='$\\mu$ = '+str(mu))\n",
"\n",
"plt.semilogy()\n",
"plt.legend(loc = 'upper right')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first plots show the trade-off in the choice of $\\mu$. For $\\mu=2$, the treads are short; the number of Newton steps required to re-center is around $2$ or $3$. But the risers are also short, since the duality gap reduction per outer iteration is only a factor of $2$. At the other extreme, when $\\mu=200$, the treads are longer, typically around $7$ Newton steps, but the risers are also much larger, since the duality gap is reduced by $200$ in each outer iteration. \n",
"\n",
"The plots also point out that for a too small value of $\\mu$ the convergence of both the dual and primal solution takes more global Newton iterations (around $70$). Therefore an appropriate choice for $\\mu$ would be somewhere between $20$ and $200$: $50$ for instance. "
]
}
],
"metadata": {
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}