{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Le chemin des solutions de l'estimateur SLOPE" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import time\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Notions liées à l'estimateur SLOPE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**La norme $\\ell_1$ ordonnée :**\n", " $J_\\lambda(b)=\\sum_{i=1}^{p}\\lambda_i|b|_{\\downarrow i}$\n" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "def sorted_L1_norm(b,Lambda):\n", " b = - np.sort(-np.abs(b)) \n", " return np.sum(Lambda*b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**La norme $\\ell_1$ ordonée duale :** \n", "$$ J^*_\\lambda(v)=\\max\\left\\{\\frac{\\|v\\|_{(1)}} \n", "{\\lambda_1},\\frac{\\|v\\|_{(2)}}{\\sum_{i=1}^{2}\\lambda_i},\\ldots,\\frac{\\|v\\|_{(p)}}{\\sum_{i=1}^{p}\\lambda_i}\\right\\}\n", "$$" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "def dual_sorted_L1_norm(v,Lambda):\n", " v = - np.sort(-np.abs(v)) \n", " v_cum = np.cumsum(v) \n", " Lambda_cum = np.cumsum(Lambda) \n", " return np.max(v_cum/Lambda_cum)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Le schéma du SLOPE ${\\rm schm}(b)$ pour $b\\in \\mathbb{R}^p$ :**" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "def pattern(b, tol=1e-10):\n", " sign = np.sign(b)\n", " perm = np.argsort(-np.abs(b))\n", " b = np.abs(b)[perm]\n", " epsilon = min(tol*(b[0]-b[-1])/len(b), tol)\n", " jump = b - np.append(b[1:],0) > epsilon\n", " q = np.flip(np.cumsum(np.flip(jump)))\n", " m = np.empty_like(b, dtype=int)\n", " m[perm] = q * sign[perm]\n", " return m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**La matrice du schéma $U_m$ pour $m\\in \\mathcal{P}^{\\rm slope}_p$:**" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "def pattern_matrix(m):\n", " k = int(np.linalg.norm(m, ord=np.inf))\n", " U_m = np.empty((len(m),k))\n", " for j in range(k):\n", " U_m[:,j] = np.sign(m)*(np.abs(m) == k-j) \n", " return U_m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithmes 1 et 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Algorithme 1 :** " ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "def face_pattern(z, Lambda, tol=1e-10):\n", " z = z / dual_sorted_L1_norm(z,Lambda) \n", " sign = np.sign(z)\n", " perm = np.argsort(-np.abs(z))\n", " z = np.abs(z)[perm]\n", " r = np.cumsum(z) / np.cumsum(Lambda)\n", " l = np.isclose(r, 1, rtol=tol)\n", " q = np.flip(np.cumsum(np.flip(l))) \n", " m = np.empty_like(z, dtype=int)\n", " m[perm] = q * sign[perm]\n", " return m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Algorithme 2:**" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "def affine_components(X, y, Lambda, m):\n", " X_tilde = X @ pattern_matrix(m) \n", " m = - np.sort(-np.abs(m)) \n", " Lambda_tilde = pattern_matrix(m).T @ Lambda\n", " P = X_tilde.T @ X_tilde\n", " a_s = - np.linalg.solve(P, Lambda_tilde)\n", " b_s = np.linalg.solve(P, X_tilde.T @ y)\n", " a_g = - X.T @ X_tilde @ a_s\n", " b_g = X.T @ (y - X_tilde @ b_s)\n", " return a_s, b_s, a_g, b_g\n", "\n", "#######################################################################################\n", "\n", "def gamma_fuse(a_s, b_s, gamma, gamma_min = 0, tol=1e-10):\n", " a = np.append(a_s, 0)\n", " b = np.append(b_s, 0)\n", " Gamma = - (b[1:]- b[:-1]) / (a[1:] - a[:-1])\n", " epsilon = tol*(gamma+1)\n", " return np.max(Gamma[(0 < Gamma) & (Gamma < gamma - epsilon)], initial=gamma_min)\n", "\n", "#######################################################################################\n", "\n", "def gamma_split(a_g, b_g, gamma_left, gamma_right, Lambda, tol=1e-10, count_max=1e2):\n", " K = np.array([True]); count = 0\n", " L_cum = np.cumsum(Lambda)\n", " while K.any() and count < count_max:\n", " count += 1\n", " g = a_g * gamma_left + b_g \n", " sign = np.sign(g)\n", " perm = np.argsort(-np.abs(g))\n", " A_cum = np.cumsum(sign[perm]*a_g[perm])\n", " B_cum = np.cumsum(sign[perm]*b_g[perm])\n", " epsilon=tol*(gamma_left+1)\n", " K = A_cum * gamma_left + B_cum - L_cum * gamma_left > epsilon\n", " gamma_left = min(np.max(B_cum[K] / (L_cum[K] - A_cum[K]), initial=gamma_left), gamma_right)\n", " if count > count_max:\n", " print(\"attention : noeuds adjacent trop proche\")\n", " return gamma_left, count > 1\n", "\n", "#######################################################################################\n", "\n", "def pattern_kinks(X, y, Lambda, ratio=1e-2, k_max=1e3, tol=1e-10):\n", " t_in = time.time()\n", " grad = X.T @ y\n", " gamma = dual_sorted_L1_norm(grad, Lambda); Gamma = [gamma]\n", " m = face_pattern(grad, Lambda, tol); M = [m]\n", " t_out = time.time(); T = [t_out - t_in]\n", " gamma_min = ratio*gamma; k = 0; Split = [True]\n", " while gamma > gamma_min and k < k_max:\n", " k += 1\n", " a_s,b_s,a_g,b_g = affine_components(X, y, Lambda, m)\n", " gamma_f = gamma_fuse(a_s, b_s, gamma, gamma_min, tol)\n", " gamma_s, split = gamma_split(a_g, b_g, gamma_f, gamma, Lambda, tol)\n", " if gamma_s > gamma - tol:\n", " print(\"attention : noeuds adjacent trop proche\")\n", " break\n", " if split:\n", " gamma = gamma_s\n", " grad = a_g * gamma + b_g\n", " m = face_pattern(grad, Lambda, tol)\n", " t_out = time.time()\n", " else:\n", " gamma = gamma_f\n", " s = a_s * gamma + b_s\n", " z = pattern_matrix(m) @ s\n", " m = pattern(z, tol)\n", " t_out = time.time()\n", " Gamma.append(gamma)\n", " M.append(m)\n", " T.append(t_out - t_in)\n", " Split.append(split)\n", " return Gamma, M, T, Split" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Chemin des solutions de l'estimateur SLOPE :** à partir des noeuds $\\gamma_0>\\dots>\\gamma_{r+1}=0$ et des schémas\n", "$m^{(0)},\\dots,m^{(r)}$ le calcul du chemin des solutions de l'estimateur SLOPE s'effectue via le code suivant " ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "def solution_SLOPE(X, y, Lambda, Gamma, M):\n", " Sol = [np.zeros(X.shape[1])]\n", " for k, gamma in enumerate(Gamma[1:]):\n", " m = M[k] \n", " a_s,b_s, *_ = affine_components(X, y, Lambda, m)\n", " s = a_s * gamma + b_s\n", " sol = pattern_matrix(m) @ s\n", " Sol.append(sol)\n", " return Sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Chemin des solutions de l'estimateur SLOPE donné à l'exemple 3." ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD7CAYAAAB37B+tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq7ElEQVR4nO3deXhV1bn48e/KACGMYR7DVAWVUWYQPEq4Sis40AqtWEUx1mKVKu2Va0t/DrfFghWsRYWIKcSrsZZ6i614G5RUCGBAIqKASJBBZRTDEIZA1u+PEyJmIPucs8/Ze628n+fxMTvZ7zl7vefNYufdZ5+ltNYIIYQwV5zXByCEECIyMpELIYThZCIXQgjDyUQuhBCGk4lcCCEMJxO5EEIYLuKJXCnVTSlVcN5/R5RSU104NiGEEA4oN99HrpSKBz4HBmmtd7r2wEIIIaqV4PLjjQS21zSJN2/eXHfq1MnlpxZCCLutX7/+oNa6RcXvuz2RTwBermmnTp06sW7dupAffPfu3QB06NAh5Fg34v3IxjG5QfIinDCtTpRSVZ4ku9ZaUUrVAb4ALtNa76vi5+lAOkBqamq/nTtD77wEAgEAVqxYEdYxRhrvRzaOyQ2SF+GEaXWilFqvte5f8ftunpGPBt6vahIH0FrPB+YD9O/fP6x/PX71q1+Ff3QuxPuRjWNyg+RFOGFLnbh5Rv4K8JbW+sWa9u3fv78Op7UihBC1WXVn5K68j1wplQyMApa48XjVKSwspLCw0LN4P7JxTG6QvAgnbKkTV99+6FS4Z+TSI6/MxjG5QfIinDCtTmLRI4+6Rx55xNN4P7JxTG6QvAgnbKkTo87IhRCiNotqjzxWXv6/Nfxm0f+FHb9161a2bt3q4hF5z8YxuUHyIpywpU6MOiPv0msgXxadZFvBGtqnJIccb1o/zAkbx+QGyYtwwrQ6saJH/tSsJ5jy0npeXPUZv77u0pDjf/vb30bhqLxl45jcIHkRTthSJ0adkQPc/8oGlm/eT970q2mUlOjykQkhhH9Z0SPftGkTgRYnOXbqDK+8tyus+E2bNkXhyLxj45jcIHkRTthSJ0adkZ/rZ7X50Uw+O3Scf//yKhLjnf9bZFo/zAkbx+QGyYtwwrQ6saJHPmvWLACONkzljsx1/GPjl9zQt13I8TaxcUxukLwIJ2ypE6POyM8pLdWMeiqXpMR43vjZFSilXDw6IYTwJyt65AUFBRQUFBAXp5g8vAsffXGE1YWHQo63iY1jcoPkRThhS50YdUZ+fj/rZMlZrnjibXq1b8LC2weEHG8LG8fkBsmLcMK0OrGiRz5nzpzyr5MS47l1cCeeyvmET/cf5TstG4YUbwsbx+QGyYtwwpY6MeqMvKJDx04xdObb3Ni3HTPH9XLhyIQQwr+s6JHn5+eTn59fvt2sQV3G9WvPkg2fc+DoqZDjbWDjmNwgeRFO2FInRp2RV9XP2n7gGCOfzOW+kRfxwKiLQ443nY1jcoPkRThhWp1EtUeulGoCZAA9AA3cobVe7cZjn++ZZ56p9L2uLRqQdkkrstbs5J4ru1KvTnxI8aazcUxukLwIJ2ypE1fOyJVSfwbe1VpnKKXqAMla66+r29/tzyNfW3iI8fPX8PgNPZg4uKNrjyuEEH4StR65UqoRMAJ4AUBrffpCk3gk8vLyyMvLq/T9gZ2b0qt9Yxau3EFpafX/MFUXbzIbx+QGyYtwwpY6ifiMXCnVB5gPfAz0BtYD92utj1cXE401O//+wRfc9/IGFvy4P6MubRVyvKlsHJMbJC/CCdPqpLozcjcm8v7AGmCY1nqtUmoucERr/esK+6UD6QCpqan9du7cGfJznVvJo1u3bpV+duZsKVfOWkG7lHq8eveQkONNZeOY3CB5EU6YVifRnMhbA2u01p3KtocDD2mtv1ddTLTW7Mx4t5DH/7GZv987jF7tm7j++EII4aWo9ci11nuB3Uqpc/+kjSTYZnFdbm4uubm51f58/IAONKybwIJ3d4QVbyIbx+QGyYtwwpY6cetdK30Ivv2wDlAITNJaH65u/2j0yM/57T8388LKHeT+IlBpXU/T+mFO2DgmN0hehBOm1UnUWivhCHciLywsBKBLly7V7vPF1ycY8ft3uG1op0rrejqJN42NY3KD5EU4YVqdWDGROyXregohbGTFZ63k5OSQk5NT4353De9S5bqeTuNNYuOY3CB5EU7YUidGnZGH0s/64fw1ldb1NK0f5oSNY3KD5EU4YVqdWNFa2b17NwAdOnSocd+3t+zjjsx1zBnfp3xdz1DiTWHjmNwgeRFOmFYnVkzkoZB1PYUQtrGiR75s2TKWLVvmaN+q1vUMJd4UNo7JDZIX4YQtdWLUGXmo/ayK63qa1g9zwsYxuUHyIpwwrU6saK3s3bsXgNatWzuOmZuzjadyPiHngRE0KD0ecrzfhZOT2kDyIpwwrU6saK20bt065IRPHJxK3YQ4Mt7dEVa839k4JjdIXoQTttSJKysExcrSpUsBGDNmjOOYc+t6vrZ+Dz3OfkrjeokhxftdODmpDSQvwglb6sSo1kq4/axz63omvPko7VPqGdMPc8K0Hl+sSF6EE6bViRU98oMHDwLQvHnzkGMn/3kd723+jH/eP5z2bapeeMJEkeTEZpIX4YRpdRLVxZdjJZJk3zW8Mzmb95G76yS3tHHxoDxmSgHGmuRFOGFLnRh1sXPJkiUsWbIkrNiBnZvS4mABv56zkCXv77ng2p4miSQnNpO8CCdsqROjWiuR9rMGDh3O9gPHaDjucXp3aMKM6y6lX8eUsB7LL0zr8cWK5EU4YVqdWNEjLyoqAqBx48ZhPW9RURGlpZrlhcf4/bIt7D96irG92/LQ6O60bVIvrMf0WqQ5sZXkRThhWp1EdSJXSn0GHAXOAmeqeqLzxeKzVmpy/NQZnl2xnfnvFhKnIH1EV35yZReS6xh12UAIUYvEYiLvr7U+6GT/cCfy7OxsAMaPHx9ybHXxew4XM/PNLbyx8UtaN0riP0d34/re7YiLM+NDtiLNia0kL8IJ0+rEiok80n7WheLzP/uKR5d+zIefF9GnQxNmjLmUy1P93z83rccXK5IX4YRpdRLtiXwHcBjQwPNa6/lV7JMOpAOkpqb227lzZ8jPU1xcDEBycnINe4YXX1qqWbLh8/L++fV92vKf1/q7fx5pTmwleRFOmFYn0Z7I22qtv1BKtQT+BfxMa/3v6vb3Q4/8QqR/LoTwo6h+aJbW+ouy/+8H/gYMdONxK8rKyiIrKyvq8fXrJjDtmm68/eCVpF3SiqeXb+Pq2bn8bYP/3n8eaU5sJXkRTthSJxGfkSul6gNxWuujZV//C3hUa13tp7X7sUd+IX7un5vW44sVyYtwwrQ6iVprRSnVheBZOARv+f8frfV/Xygm3Im8pKQEgMTExJBjI433a/880pzYSvIinDCtTqy4IcgPpH8uhPCKFQtLZGZmkpmZ6Vk8VN8/f33D5570z90Yk40kL8IJW+rEqDNyr3rkF+J1/9y0Hl+sSF6EE6bVibRWoqi0VPPX9/fw+7e2csBH/XMhhF2saK34VVyc4gf9O7BiWoB7r/oOb27ay9VPruCpf31C8ekzXh+eEMJyRk3kCxYsYMGCBZ7F1+Rc/3z5A1cy8pJWzI1B/zzaYzKV5EU4YUudGNVaSUtLAyAnJyes5400PlSx6J/HekymkLwIJ0yrE+mRe0T650IIt0iP3CPSPxdCRJtRE/m8efOYN2+eZ/GRiFb/3Msx+ZnkRThhS50Y1VoZPXo0AG+++WZYzxtpvJvc6p/7aUx+InkRTphWJ9Ij96GK/fMb+rTll9I/F0JUQ3rkPlSxf/7Psv75nJxPOHH6rNeHJ4QwhFET+dy5c5k7d65n8dFSsX8+J2cbVz+5wlH/3K9j8prkRThhS50YNZEvX76c5cuXexYfbR2aJvOnH13OX34yhOYN6jI1u4Cbns3j/V2Hq43x+5i8InkRTthSJ9Ij9ynpnwshKop6j1wpFa+U2qCUesOtx6zNzvXP35kWYMpVXaV/LoSolputlfuBzS4+XiWzZ89m9uzZnsV7oUHdBH5xTfdq++cmjikWJC/CCVvqxJVlbZRS7YHvAf8NPODGY1Zl9erVnsZ76Vz//LYhX/HoGx8xNbuAzLzPOLo8lybJdbw+PN8x+bUWsWNLnbjSI1dKvQb8DmgITNNaX3eh/aVHHpmK/fOb+rbj8Rt7yHJzQlguaj1ypdR1wH6t9foa9ktXSq1TSq07cOBApE9bq53fP/9poCuvF3zObQvf4+jJEq8PTQjhATd65MOAsUqpz4BXgKuVUlkVd9Jaz9da99da92/RokVYTzRz5kxmzpwZ9oFGGu83DeomUFrwOsOLV7Fh19fckrGWw8dPe31YvmDbay2iw5Y6ifhvca31dGA6gFIqQLC1MjHSx61KQUGBp/F+dG5Mz//mfu556X0mzF/D4skDadkwydsD85iNr7Vwny114ur7yM+byKVH7oFVnx5k8p/X0bpxElmTB9FO3nMuhFVi8lkrWusVNU3iInqGfac5i+8cyMGjp7j5udV8dvC414ckhIgBo27Rf+yxx3jsscc8i/ejimPq36kpL6cPpvj0GW5+fjXb9h318Oi8Y+NrLdxnS50Y9X61rVu3ehrvR1WNqUe7xmTfPYRbMtZy8/OrWXznIHq0a+zB0XnHxtdauM+WOpHPWrHYZwePc0vGWo6cKCHzjgH069jU60MSQkRAPo+8FurUvD6v/mQIzRvWZWLGe6z69KDXhySEiAKjJvIZM2YwY8YMz+L9qKYxtWtSj+y7B5PaNJlJmfks37wvhkfnHRtfa+E+W+rEqB757t27PY33IydjatkwiVfSB/Pjhe9x9+L1zJnQh+t6tY3B0XnHxtdauM+WOpEeeS1y5GQJd2bms37nYWaO68XN/Tt4fUhCiBBIj1zQKCmRP98xkGHfac4vX9vIotWfeX1IQggXGDWRT58+nenTp3sW70ehjim5TgILftyftEtaMeN/P+LZFdujeHTesfG1Fu6zpU6M6pEfOnTI03g/CmdMSYnxPDvxch549QOeWLaF4tNneGDUxSilonCE3rDxtRbus6VOpEdei50t1fzXkg/JXrebO4Z15tfXXWLVZC6EbarrkRt1Ri7cFR+n+N1NPalXJ56Fq3ZQfPoM/31jT+LjZDIXwiRGTeTTpk0DCHuNvUjj/SjSMcXFKX4z5lLq143nT+9s50TJWWb/oDeJ8UZdPqnExtdauM+WOjFqIj9x4oSn8X7kxpiUUvzimu4k10lg1ltbKT59lmd+1Je6CfEuHKE3bHythftsqRPpkYtvyVy1g/+39GOGX9Sc+bf2p14dcydzIWwj7yMXjtw+rDO/H9eLVZ8elHVAhTCEG4svJyml3lNKfaCU+kgp9YgbB1aVqVOnMnXqVM/i/SgaY7p5QAfmTujL+7sOMzFjLV8Xm7cOqI2vtXCfLXXixhn5KeBqrXVvoA9wrVJqsAuPKzw0pndbnp3Yj81fHmXC/DUcOHrK60MSQlTD7TU7k4GVwD1a67XV7Sc9cnOs3HaQuxYF1wF9afIg2so6oEJ4Jqo9cqVUvFKqANgP/OtCk7gwyxUXfbMO6A+eW83OQ7IOqBB+48pErrU+q7XuA7QHBiqlelTcRymVrpRap5Rad+DAgbCeZ8qUKUyZMiXs44w03o9iMab+nZryP3cN5vjpM/zgOTPWAbXxtRbus6VOXH3Xitb6a2AFcG0VP5uvte6vte7fokWLsB6/Xr161KsX/p/2kcb7UazG1LN9Y7LTh6CB8fPXsOnzoqg/ZyRsfK2F+2ypk4h75EqpFkCJ1vprpVQ94P+AJ7TWb1QXIz1yc+04eJxbFqzh6KkzZE6SdUCFiKVo9sjbAO8opTYC+QR75NVO4sJsnZvX5y/3DKVZ/Trc+sJ75Mk6oEJ4LuKJXGu9UWvdV2vdS2vdQ2v9qBsHVpX09HTS09M9i/cjL8bUrkk9Xr17CO1T6nF7Zj5vb/HfOqA2vtbCfbbUiVF3djZr1oxmzZp5Fu9HXo2pZaMkXkkfwsWtGpC+aD3/2PhlzI/hQmx8rYX7bKkT+awVEZEjJ0u448V83t91mCfG9eIHsg6oEFEjn7UioqJRUiKL7hzI0K7N+cVrG1ks64AKEXNGTeSTJk1i0qRJnsX7kR/GlFwngYzb+pN2SUt+/b8f8Vyu9+uA+iEvwv9sqROjPo+8Q4fI/myPNN6P/DKm4Dqg/fh5dgEz39xC8akz/NzDdUD9khfhb7bUifTIhavOlmqmL9nIq+v2cOcVnfnV92QdUCHcImt2ipiIj1PMvKkXyXUSeGFlcB3Qx2+QdUCFiCajJvKJEycCkJWV5Um8H/lxTOfWAU2uE8+8Fds5cTq4DmhCDNcB9WNehP/YUidGTeTdunXzNN6P/DompRS/vLY79et+sw7oH2O4Dqhf8yL8xZY6kR65iLoXV+3gEVkHVIiIyfvIhWcmla0DulLWARUiKoyayCdMmMCECRM8i/cjU8Z0bh3Q9TFaB9SUvAhv2VInRvXI+/Tp42m8H5k0prG925KUEMe9/7OBCfPXsPjOQbRoWDcqz2VSXoR3bKkT6ZGLmHt32wHSF62nTeMkXrprEG0am//B/kLEgvTIhW8Mv6gFi+4cyAFZB1QIVxg1kY8bN45x48Z5Fu9Hpo5pQKemvHTXII6dOsPNz6/m0/3urgNqal5EbNlSJxH3yJVSHYBFQGugFJivtZ4b6eNWZciQIZ7G+5HJY+rVvgnZ6UO4JWMtNz+/hkV3DKRHu8auPLbJeRGxY0uduLFmZxugjdb6faVUQ2A9cIPW+uPqYqRHLs737XVAB9KvY4rXhySEL0WtR661/lJr/X7Z10eBzUC7SB9X1B6dm9fn1Z8MKVsHdC1522UdUCFC4WqPXCnVCegLrHXzcc8ZO3YsY8eO9Szej2wZU/uUZF69ewjtmtRj0ov5vLNlf0SPZ0teRHTZUieuvY9cKdUA+CswVWt9pIqfpwPpAKmpqWE9x8iRIyM5xIjj/cimMbVslET23UP48cK1pC9ex9wJffluzzZhPZZNeRHRY0uduPI+cqVUIvAG8JbW+g817S89cnEhR06WMOnFfDbsOsys7/dmXL/2Xh+SEL4QtR65Cq4a8AKw2ckkLkRNGiUlsvjOgQzp2owH//IBi9fs9PqQhPA1N3rkw4BbgauVUgVl/33XhcetZPTo0YwePdqzeD+ycUwQXAf0hdsGMLJ7S379+iaeD3EdUFvzItxlS51E3CPXWq8EYrL8y5gxYzyN9yMbx3ROUmI8z93aj6nZBfzuzS0cP32Wn6dd5GjpOJvzItxjS53IZ60I3ztbqnnorxv5y/o9TL6iMw/LOqCilpI1O4Wx4uMUT4zrRf26CWSs3EFxyVkev74HcbIOqBCAYRN5WloaADk5OZ7E+5GNY6rKuXVA69WJ59mydUBnfb9XteuA1pa8iMjYUidGTeTjx4/3NN6PbBxTdZRS/Oe13WlQvg7oGZ7+YdXrgNamvIjw2VIn0iMXRlq4cgePvvExV17cgucm9pN1QEWtIJ9HLqxyxxWdeWJcT/697QC3v/gex06d8fqQhPCMURN5IBAgEAh4Fu9HNo7JqfEDUpkzvg/rdh7mlgrrgNbmvAjnbKkTo3rkt99+u6fxfmTjmEJxfZ921EuMr7QOaG3Pi3DGljqRHrmwwrvbDnDXonW0bVKPlybLOqDCTlb0yEtKSigpKfEs3o9sHFM4hl/UgkV3DGL/keA6oNv3FkleRI1s+f0xaiIfNWoUo0aN8izej2wcU7gGdm7KS3cOYNDJldw0uDODenRmY+4Sjh057PWhCZ+y5ffHqB755MmTPY33IxvHFJazJbDxVXqvfIon9TaevzyRJH2YXu9M4uzbim0JXTnUrB91ugyjY9+RNGslH40r7Pn9kR65MFvJCXh/MeQ9DUW7oVVPGP4AXHo9x48dobBgBcc/eZeG+/PpemozSSr4Z/Qu1Y69TfqgOg6lba+RtO3UDRVn1B+oohaqrkdu1EReXFwMQHJycljPG2m8H9k4JkdOFkH+C7BmHhw/AB0GwfBpcNEoUKrKvJw+dZLCjSv5enMuSV++R5cTH9KI4wDspym7G/bmTPvBtLjsKjpd0p+4eLnJyHam/f5YMZGfe7/nihUrwnreSOP9yMYxXdDxg7DmWXhvAZwqgq4jYfiD0HEonPeJiE7yUnr2LDu3rGP/RytI2L2GDkcLaMlXAByhPjvq9aC4zSBSuo+gS+/h1KmbFM2RCQ+Y9vsT1U8/VEotBK4D9mute7jxmFW55557PI33IxvHVKWiPZD3DKzPhDMn4ZIxwRZK275V7u4kL3Hx8XS+bBCdLxsEgC4t5Yudn/DFxuWU7syj1dcF9C58Ggqf5uQ/EvmobneOtBxAg4uG07nvVTRolOLmCIUHbPn9cWvNzhHAMWCRk4lceuTCsUPbYeVT8MErgIZe42HYVGhxcUye/qv9n/PZhrc5XbiSZofW07lkOwmqlDM6jh0JXTjUrB91uw4jtY9cQBXRF/XWilKqE/BGNCfyoqIiABo3bhxyrBvxfmTjmADY+yG8+wf4+HWIrwOX/xiG/gyapDoKj1Zejh05zI7qLqDGtWNv477QcSjteo2kbaeL5QKqz5n2+2PFRC498sqsG9OuNcEJfNtbUKchDJwMg38KDVqG9DCxyktNF1B3NezD2faDaXFZQC6g+pBpvz+erxCklEoH0gFSU52dVVV03333RXQMkcb7kRVj0hq2Lw9O4DtXQXIzuPpXMOAuqNckrIeMVV7q1E2i+4A0GBBcoKD07Fl2bFnH/k0rSNizhtSjBbTc/DZs/q1cQPUhK35/MOyMXFimtBS2LIV3n4QvP4CGbWHYfcE2Sp36Xh+dK3RpKV/u/ITPP1iO3plH66INpJZ+DsBJnch2uYAqQmBFa+XgwYMANG/ePORYN+L9yMgxnS2BD/8SvIh58BNo2gWu+Dn0mgAJdVx5Cj/n5at9e9hZ8Dantq+k2Vfr6VKynXil5QKqB/xcJ1WJ6kSulHoZCADNgX3Ab7TWL1S3v/TI3WPUmEpOwIYsWDX3vLswfw6X3gBx7vaOTcqLkwuoquNQ2vUeSZuOcgHVTSbVCUS5R661/qEbj1OTBx980NN4PzJiTCePwLoXYPWfvrkL83t/KL8LMxqMyEuZBo1S6DniRhhxIxC8gLrlvAuo3Q+/Q6PDb0DBuTtQ+5TdgSoXUCNlUp1ciFF3dgrDHD8Ia5+DtfMveBemuLDyO1DLLqCefwdqEfXZUa8nJ9oMJOWSK+nS6wq5gGoxK27R37t3LwCtW7cO63kjjfcjX46p6HPI+6PjuzCjwZd5ccmFLqCe0HUorNudIy37ywVUB0yrEysmcumRV+arMR3aDqvmQMHLoEuDd2FeMRVadIv5ofgqLzFwaN8edhUs59T2VTQ7tJ4uZ+QCqhOm1Ynn7yN3w0MPPeRpvB/5Ykzn34UZlwj9bg/ehZnS0bND8kVeYqhZq/Y0u+Y24Dag7ALqhnc4vm0lDffn02ffEpL2Z8NquYB6PlvqxKgzcuEzu9YG3wN+7i7MAXcG78Js2MrrIxMV1PwRvsELqC17BOjYXS6g+pUVrZXdu3cD0KFDh7CeN9J4P4r5mLSG7W+X3YW5Euo1hSE/jeguzGiw8bV2k5MLqCfbDKSJ5RdQTasTKyZy6ZFXFrMxlZbCljfK7sIs8P1dmDa+1tHk5ALq0ZYDqH/xcLr0vYr6DZt4e8AuMa1OrJjIc3JyAEhLSwvreSON96Ooj+lsCXz4WtldmFvPuwtzPCTUjc5zusDG1zrWLngBNfHcGqhX0Knv1TRt2c7rww2LaXVixUQuYqj8LsynoWgXtOpRthbmDa7fhSnM8M0F1HdptD+fLqe2yB2oMWbFRF5YWAhAly5dwnreSOP9yPUxld+FOQ+O74f2A2HENLjoP4y6icfG19pvTp88QeHGVXy9JZd6X66l84lNxl1ANa1OrJjIpUdemWtjOn4I1j4L780PLmzc9eqyuzCHGTWBn2Pja+13Jl5ANa1OrJjIc3NzAbjyyivDet5I4/0o4jEVfQ6ry9bCLCkO3oV5xQPQ7nL3DtIDNr7WpjHhAqppdWLFRC5c5KO7MEXtURsuoEaTFRP51q1bAejWLbzJJtJ4Pwp5THs3wco/wEd/C96FeW4tTA/vwowGG19rG527gHqs7AJq1/MuoO6Ma8++JuddQE29yPULqKbViRUTufTIK3M8pt3vBd8D/smyWnEXpo2vdW1w6mQxOzau4vCWf5Nc4QLqPpqxu2EfzrYf5NoFVNPqxIqJPC8vD4ChQ4eG9byRxvvRBcekNRS+E7wL87N3g3dhDv5pcEHjenZ/Ip6Nr3Vt9M0F1HdI2LPW9QuoptVJtFcIuhaYC8QDGVrrmRfaX3rkUVZaClv/ETwD/2JD8C7MoT+Dfrf58i5MIZwy4QJqNEVtIldKxQOfAKOAPUA+8EOt9cfVxYQ7kW/atAmAHj1qXBY0KvF+9K0xnS2BTX8NnoEf3AopnYN3Yfae4Ou7MKPBxtdaVC2SC6im1Uk0J/IhwP/TWl9Ttj0dQGv9u+pipEfunkAgALqUFbNu/eYuzJaXfXMXZrxRn1TsGhtfa+FMKBdQf3T73YA5dVLdRO7GJeB2wO7ztveUfa9aW7duJTMzE4CSkhICgQBZWVkAFBcXEwgEyM7OBqCoqIhAIMCSJUuYNWsWDz/8MIFAgKVLlwLBFT4CgQDLli0Dgp9mFggEyj9DobCwkEAgQG5uLrPuGsW93Q8R6JZC3kP9ISONTY8MJtAthfyHB0JGGgUzBhHolkLBjEGQkUb+wwMJdEth0yODISONvIf6E+iWwtbHh0BGGrm/6EegWwqFvxsGGWnkPHA5gW4p7P79FZCRxrKpfQl0S2Hvk8MhI42lP+tNoFsKB58aARlpLPlpLwLdUih6OgAZaWTf3ZNAtxSK/3QVZKSRNfkyAt1SKHnuashII3PSpQS6pUBGGmSkMSx5J8WF78E/p0HD1szTtzJ6aWPo+X2IT2Du3LmMHTu2PPezZ89m3Lhx5dszZ85kwoQJ5duPPfYYEydOLN+eMWMGkyZNKt+ePn066enp5dvTpk1jypQp5dtTp05l6tSp5dtTpkxh2rRp5dvp6elMnz69fHvSpEnMmDGjfHvixIk89thj5dsTJkxg5sxvOnXjxo1j9uzZ5dtjx45l7ty55dujR49m3rx5zJo1i1mzZpGWlsaCBQvKfx4IBMKqPQiuuB5u7UGw7gOBQHlfdtOmTQQCAfLz8wEoKCggEAhQUFAAQH5+PoFAoPysMS8vj0AgUP5Oi9zcXAKBQPndiTk5OQQCgfJP9Fu2bBmBQKB8FZylS5cSCATKV45fsmQJgUCAoqIiALKzswkEAhQXFwOQlZVFIBCgpCQ4CWZmZpb/AwmwYMGCb31Gybx58xg9enT5tle116BRCj2vvImsDad58avBqOm72PLd17h+RVcefaeY7l+9zYAN03n4xp60PbKee7sfYOvjgxjbuzn3j2zP1scHsfXxQXyvZzMe/I8O5dv/cVlTfnltavn21d1TePh7Hcu3R1zchN+M6VS+PbRrYx67oXP59rn/tqxbjtvcOF2r6ra/Sqf5Sql0IB2gbt3w/swfMGBAeRGGrOQEA/Zk0KhpYvCzQuokQ92GUEdX2D5bYftM2Xb9su3TFbZPnbddHxJPnLedDInFZdsNoG4SJB4LbtdtAHXrQsLR87YTIaHovO0ESPi6bLshJMRBYtI320Cn1k1otOcU3P5q8C7MZ58FtoeXI4sMGDDA60MQPlE3KZnuA0fR6uK/0KxZMxr8+nEKN6/jwIoptE48xsWdWnEKOBuXwJm4OpxKCF5HqrhdqhI4e/52XHyFn8dzJq5utdvnJMW5/1eyUa2Vc2cqffr0Cf1A856hYNFDcMNz9Ll2Ys37GyKinFhM8iKcMK1OotkjTyB4sXMk8DnBi50/0lp/VF1MzHvkJSdhbi8CGYegdU9j+mFOSC+4apIX4YRpdRK1NTu11meUUvcCbxF8++HCC03ikZgzZ054gRsWw7F9zHnyD9Cun6vH5LWwc2I5yYtwwpY6MeqGoLCcOQ1/vBwatYU73jLyk/yEEAKi+66VmMnPzy+/wu/Yxmwo2g0jfkH+unWhx/tcWDmpBSQvwglb6sSoM/KQ+1lnz8CfBgTf4ZGeS+Cqq0KLN4BpPb5YkbwIJ0yrk6j1yGPpmWeeCS3go7/BV4UwPguUCj3eADaOyQ2SF+GELXVi1Bl5SEpL4dkhgIJ78kDWDxRCGM6KHnleXl75XXE12vIGHNgSXG+ybBIPKd4QNo7JDZIX4YQtdWLUGbnjfpbW8PwIOH0c7s0vX/XdtH6YEzaOyQ2SF+GEaXViRY/8+eefd7bjtn/B3o1w/Z/KJ/GQ4g1i45jcIHkRTthSJ0adkTuiNbwwCo7ug/veh/jE6DyPEELEmBU98tzc3PJPkqvWjn/DnvzgQsIVJnFH8YaxcUxukLwIJ2ypE6POyB31szKvg0Ofwn0FwU8KDDXeMDaOyQ2SF+GEaXViRY984cKFF95h15rg2pTX/K7SJO4o3kA2jskNkhfhhC11YtQZeY2yvh9co3LqRlmbUghhHSt65Dk5OeWrr1TyxQb49F8wZEq1k/gF4w1l45jcIHkRTthSJ0adkV+wn/XKLcG2ytRNkNQo9HhD2TgmN0hehBOm1YkVPfLFixdX/YN9Hwfv5LzyoWon8QvGG8zGMblB8iKcsKVOjJrIO3ToUPUP3n0yuCbmoLvDizeYjWNyg+RFOGFLnUTUI1dK/UAp9ZFSqlQpVel0323Lli0rX7G83MFP4aMlMGAyJDcNPd5wNo7JDZIX4YQtdRJRj1wpdQlQCjwPTNNaO2p8u9ojf30KbPorTP0QGrQIPd5wNo7JDZIX4YRpdRKVHrnWenPZg0fyMI698sor3/7G4Z2w8RUYcFeNk3iV8RawcUxukLwIJ2ypE6N65K1bt/72N1bNBRUHQ38WXrwFbByTGyQvwglb6qTGiVwplQNUNdqHtdb/6/SJlFLpQDpAamqq4wM839KlSwEYM2YMHPkSNiyGPrdA43ahx1vCxjG5QfIinLClTlx5H7lSagWx7pEv+y9Y+1zwEw5TOoUebwkbx+QGyYtwwrQ6seJ95K+99lrwi+MHYd1C6DXe8ST+rXiL2DgmN0hehBO21ElEE7lS6kbgj0AL4B9KqQKt9TWuHFkVmjdvHvwi5xE4cxKGPxBevEVsHJMbJC/CCVvqJNJ3rfwN+JtLx1KjJUuWwKlj3LRjAVx2IzS/KPR44KabborG4XnCxjG5QfIinLClTsz7rJWvd7HihkNwTx60uiz0eMzphzlh45jcIHkRTphWJ9X1yI2ayIv27YZnh9K4+wiY8FLo8UVFADRu3DjkWL+ycUxukLwIJ0yrEysudjbe9hpwBIY/GF68IS9WKGwckxskL8IJW+rEqM8jz15ZSPbxQdDu8vDis7PJzs52+ai8ZeOY3CB5EU7YUidGtVYi7WeZ1g9zwsYxuUHyIpwwrU6s6JEXFxcDkJycHNbzRhrvRzaOyQ2SF+GEaXViRY880mSb8mKFwsYxuUHyIpywpU6M6pFnZWWRlZXlWbwf2TgmN0hehBO21IlRrRXpkVdm45jcIHkRTphWJ1b0yEtKSgBITEwM63kjjfcjG8fkBsmLcMK0OrGiRx5psk15sUJh45jcIHkRTthSJ0b1yDMzM8nMzPQs3o9sHJMbJC/CCVvqxKjWivTIK7NxTG6QvAgnTKsTX/XIlVIHgJ1hhjcHDrp4ODaQnFRN8lKZ5KQyk3LSUWtdaYFiTybySCil1lX1L1JtJjmpmuSlMslJZTbkxKgeuRBCiMpkIhdCCMOZOJHP9/oAfEhyUjXJS2WSk8qMz4lxPXIhhBDfZuIZuRBCiPP4diJXSl2rlNqqlPpUKfVQFT9XSqmny36+USkV3moTBnGQk4BSqkgpVVD23wwvjjOWlFILlVL7lVKbqvl5bayTmnJSG+ukg1LqHaXUZqXUR0qp+6vYx9xa0Vr77j8gHtgOdAHqAB8Al1bY57vAm4ACBgNrvT5uH+QkALzh9bHGOC8jgMuBTdX8vFbVicOc1MY6aQNcXvZ1Q+ATm+YUv56RDwQ+1VoXaq1PA68A11fY53pgkQ5aAzRRSrWJ9YHGkJOc1Dpa638DX11gl9pWJ05yUutorb/UWr9f9vVRYDPQrsJuxtaKXyfydsDu87b3UDnpTvaxidPxDlFKfaCUelMpdVlsDs3XaludOFVr60Qp1QnoC6yt8CNja8Wvn36oqvhexbfXONnHJk7G+z7BW3iPKaW+C7wOXBTtA/O52lYnTtTaOlFKNQD+CkzVWh+p+OMqQoyoFb+eke8BOpy33R74Iox9bFLjeLXWR7TWx8q+/ieQqJRqHrtD9KXaVic1qq11opRKJDiJv6S1XlLFLsbWil8n8nzgIqVUZ6VUHWAC8PcK+/wd+HHZlebBQJHW+stYH2gM1ZgTpVRrpZQq+3ogwdf3UMyP1F9qW53UqDbWSdl4XwA2a63/UM1uxtaKL1srWuszSql7gbcIvltjodb6I6XUT8p+/hzwT4JXmT8FioFJXh1vLDjMyfeBe5RSZ4ATwARddjneVkqplwm+C6O5UmoP8BsgEWpnnYCjnNS6OgGGAbcCHyqlCsq+919AKphfK3JnpxBCGM6vrRUhhBAOyUQuhBCGk4lcCCEMJxO5EEIYTiZyIYQwnEzkQghhOJnIhRDCcDKRCyGE4f4/3EZmsNYG738AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "X=np.array([[1,0.5],[0.5,1]])\n", "Lambda = np.array([4,2])\n", "y = np.array([6,2])\n", "\n", "# Calcul des noeuds et schéma\n", "Gamma, M, T, Split = pattern_kinks(X, y , Lambda, ratio=0)\n", "\n", "# Calcul des solutions aux noeuds \n", "Sol = solution_SLOPE(X, y, Lambda, Gamma, M)\n", "\n", "# Représente le chemin des solutions \n", "fig, ax = plt.subplots()\n", "ax.plot([1.1*Gamma[0]] + Gamma, [Sol[0]] + Sol)\n", "for gamma in Gamma:\n", " ax.axvline(gamma, color='k', linestyle=':')\n", "ax.axhline(0, color='k', linestyle=':', xmax=0.95)\n", "plt.show() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithme d'homothopie pour le calcul du chemin des solutions de l'estimateur LASSO" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "def affine_components_LASSO(X, y, s):\n", " I=s!=0\n", " P=np.transpose(X[:,I])@X[:,I]\n", " B=np.transpose(X[:,I])@y\n", " A=-s[I]\n", " b_s=np.linalg.solve(P,B)\n", " a_s=np.linalg.solve(P,A)\n", " a_g = - X.T @ X[:,I] @ a_s\n", " b_g = X.T @ (y - X[:,I] @ b_s)\n", " return a_s, b_s, a_g, b_g\n", "\n", "#######################################################################################\n", "\n", "def gamma_zero(a_s, b_s, gamma, tol=1e-10):\n", " epsilon=tol*(gamma+1)\n", " vanish=-b_s/a_s\n", " return np.max(vanish[(0 < vanish) & (vanish < gamma - epsilon)], initial=0)\n", "\n", "#######################################################################################\n", "\n", "def gamma_gradient(a_g, b_g, gamma, tol=1e-10):\n", " K = np.array([True]); count = 0\n", " while K.any():\n", " count +=1\n", " g = a_g * gamma + b_g \n", " sign = np.sign(g)\n", " K = np.abs(g) > gamma + tol*(gamma+1) \n", " gamma=np.max(b_g[K] / ( sign[K]-a_g[K]), initial=gamma)\n", " return gamma, count > 1\n", "\n", "######################################################################################\n", "\n", "def pattern_kinks_LASSO(X, y, k_max=1e2, tol=1e-10):\n", " t_in = time.time()\n", " grad = X.T @ y\n", " gamma = np.linalg.norm(grad,ord=np.inf); Gamma = [gamma]\n", " s=np.repeat(0,X.shape[1])\n", " s[np.transpose(X)@y>gamma*(1-tol)]=1\n", " s[np.transpose(X)@y<-gamma*(1-tol)]=-1\n", " S=[s]\n", " t_out = time.time(); T = [t_out - t_in]\n", " k = 0\n", " while gamma>0 and k < k_max:\n", " k += 1\n", " a_s,b_s,a_g,b_g = affine_components_LASSO(X, y, s)\n", " gamma_z = gamma_zero(a_s, b_s, gamma)\n", " gamma_g, enter = gamma_gradient(a_g, b_g, gamma_z)\n", " if enter:\n", " gamma = gamma_g\n", " grad = a_g * gamma + b_g\n", " s = np.repeat(0,X.shape[1])\n", " s[grad>gamma*(1-tol)]=1\n", " s[grad<-gamma*(1-tol)]=-1\n", " t_out = time.time()\n", " else:\n", " gamma = gamma_z\n", " sol_support = a_s * gamma + b_s\n", " epsilon=tol*(1+np.linalg.norm(sol_support,ord=np.inf))\n", " sol_support[np.abs(sol_support)" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Le paramètre de pénalité de l'estimateur SLOPE est \n", "Lambda=np.sqrt(range(1,14))-np.sqrt(range(0,13))\n", "\n", "# Calcul des noeuds et schémas\n", "Gamma, M, T, Split = pattern_kinks(Xapp, yapp, Lambda, ratio=0)\n", "\n", "\n", "# Calcul du chemins des solutions du SLOPE et des composantes des moindres carrés en valeur absolue. \n", "Sol = solution_SLOPE(Xapp, yapp, Lambda, Gamma, M)\n", "\n", "Sol_abs = [np.abs(sol) for sol in Sol]\n", "ols = np.linalg.solve(Xapp.T@Xapp, Xapp.T@yapp)\n", "ols_abs = np.abs(ols)\n", "\n", "from cycler import cycler\n", "plt.rcParams['axes.prop_cycle'] = cycler('color', ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', \n", " '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf','b','m','k'])\n", "\n", "plt.plot([1.1*Gamma[0]] + Gamma, [Sol[0]] + Sol_abs)\n", "plt.plot(0, [ols_abs], 'x')\n", "plt.axvline(0, color='k', linestyle=':')\n", "plt.axhline(0, color='k', linestyle=':', xmax=0.95)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABZXUlEQVR4nO3dd3hT1R/H8fdJ0r0LhRYKZSPIEsteRXZlCCJDkSUyZAhYFH6AE0UFAZkiishwMAWkbChD9gZBkE1ZLZTumeT8/khARoECbW/Tntfz9CE3ubn5JNpvT8499xwhpURRFEWxfTqtAyiKoiiZQxV0RVGUXEIVdEVRlFxCFXRFUZRcQhV0RVGUXMKg1Qvnz59fFitWTKuXVxRFsUn79++/IaX0Se8xzQp6sWLF2Ldv3xM/79KlSwAUKVIksyMpiqLkeEKICw97LMNdLkIIvRDioBDiz3QeE0KIyUKI00KII0KIqk8b9nHefPNN3nzzzaw6vKIois16khb6u8AJwD2dx1oApa0/NYAZ1n8z3ahRo7LisIqiKDYvQwVdCOEPvAx8DgxNZ5c2wFxpuex0lxDCUwjhJ6W8mnlRLRo3bpzZh1QURckVMtrlMgl4HzA/5PHCwKW7tsOt991DCNFbCLFPCLEvMjLySXLecfbsWc6ePftUz1UURcnNHlvQhRAtgQgp5f5H7ZbOfQ9MEiOl/F5KGSilDPTxSfck7WP17NmTnj17PtVzFUVRcrOMdLnUAVoLIYIBR8BdCDFfStnlrn3CgbuHnfgDVzIvJly4MBM390p88sknd+6LurWTuNgjBAT0ycyXUhRFsUmPbaFLKUdIKf2llMWATsCm+4o5wAqgq3W0S00gJrP7z93cK3Hs2CAqVrKnQYMGRN3aybFjg3Bzr5SZL6MoimKznnocuhCiL4CU8jsgFAgGTgOJQI9MSXcXb69aVKgwmUOH+mI2F8HO7ioVK07F26tWZr+UoiiKTXqigi6lDAPCrLe/u+t+CfTPzGDp8faqxaZNKTRseIKCBTqrYq4oinIXm5rLJerWToKCDEREGLkesYqoWzu1jqQoipJj2ExB//iT/syd253SpQZToICB4sUGMHdudz7+JMu/GCiKotgEzeZyeVJlyzgwcGAEet0FLl+OIzY2hjGfRTBlioPW0RRFUXIEmynonTtPwNe3FU2aNMZkMuPh8RHLlq2gYcOGWkdTFEXJEWymywWgYcOGtG/fCoAyZfxVMVcURbmLzbTQATZv3szGjdvwK2Rg//6TbNq0iYYNG3Ljxg3OnDnD6dOnOXPmDGfOnOHEiROUL1+etm3b0rRpU5ydnbWOryiKkqWEZcRh9gsMDJRPMh/65s2b6dChA59+OpJduz5l7txb6HQ6HBwcSEpKumdff39/ihUrxrFjx4iOjsbJyYmmTZvSt29fmjdv/kQ5U5OTSIyOJiEmmsSYWxjsHfAvXwE7e9V3ryhK9hNC7JdSBqb3mM200Pfu3cvChQv5+OOPuHUric6d63H2bCp6vZ4OHTpQsmRJSpYsSfHixXF0dAQgLS2NrVu38scff7Bs2TJatmzJli1beLFKZRJjou8p1Ikx0STGxNyznRATjTEl5YEsBnsHilaoRImq1Sj+QjXc8z/dvDSKoiiZyWZa6LcdPHiAPXtb06TxMEqUeDfDz9uxYhltu/UgLTWVIU3r4Wxvd8/jQuhwcnfH2cMTZw9PXKz/Ont44uLpded2UmwMZw/u5ez+PcREXAfAJ6A4JapWp0TVaviWKo1Op3/i96UoipIRuaKFftsLL1QlOsYNszk5Q/ubjGls+mkmRzasoX+bFnw893fO2bkSMmDAPYXbyd09w4W4WOWqNOzWm6jL4Zw9sIezB/ayZ/kidi/7HSc3d4q/EMjzDRpT5PmKCJHeRJSKoiiZz+YK+t69ezl10kyhQkmP3Tch+hYrJ47l8j/Hqf7Ka9Tp2IVDsSn8sWEzU36ai4uLy1PnEEKQz78I+fyLUK31qyTFx3H+8AHO7t/Dmf27Ob51E16F/KncuAXPN2iEo6vrU7+WoihKRthcl0tQUBAxMQf45Zf+lCs39qH7XT97mj/GjyE5Lo5m/d7ludr1AdixYwd16tRhypQpDBgw4KnzP0paagqndm7n8PpQrv57EoO9A2Vr16Nykxb4liyjWu2Kojy1R3W52FxBP3bsGIcO9aDKC3Wp8PzEdPc5f+Qgy8eNwcndnTYhoyhYvOQ9j9eqVYuIiAhOnTqFXp+1/d3Xz53hyIbVnNgWRlpKMgWKl6RykxaUqxOEnfXkraIoSkblqoIOsGt3C5ydi1Gp4owHHrtx8Ty/jB6GZ0Ff2o/8DGcPzwf2Wbx4Ma+99hpLliyhXbt2T5XhSaUkJnJi22YOrw/lxqUL2Ds5U77+S1Ru0oL8RQKyJYOiKLYvVxX0HTt2cOLEcF58sSRVqvx0z2Nms4mf3+tPSlIib3wxATfv/Okew2QyUbp0aXx9fdmxY8dT5X9aUkqunDzB4Q2rObVzGyajkcLPPU/lJi0oXaMOBju7xx9EUZQ861EF3aYu/Qf43//+x/QZ/2AyPzg+/ORfW4m6Es5LPfo8tJgD6PV6hgwZws6dO7O9oAshKPxceYIHvEfvGT9T/40eJNyKInTKeL5/pztbf5lD9PVr2ZpJUZTc4bEtdCGEI7AVcMAyKmaxlPKj+/YJApYD56x3LZVSfvqo4z5tC/3kyZMcPzEc/8JQrdqyO/ebzSbmvNcfvcFA168mI3SP/lsVHx9P0aJFadiwIUuWLHniHJlJms1cOHqIw+tDObN/D1JKilWuSuXGLShRtRq6LO7nVxTFdjzrOPQU4CUpZbwQwg7YLoRYLaXcdd9+26SULZ817OOULVuWlFRfEhPP33P/yR3buHUlnFZDRzy2mAO4urrSt29fvvzyS86cOUPJkiUf+5ysInQ6ilWuSrHKVYm7eYOjm9ZydONalo8fg2u+/FRq1IyKDZvi6p1Ps4yKouR8GVkkWkop462bdtYfbTregS1btrB//417Liwym03sXPIb+YsEULpaxpelGzhwIAaDgYkT0x8towW3fPmp/dobvD3tJ1q/9z/yFS7CjoUL+L5/D1ZM+IILRw4hzWatYyqKkgNl6KSoEEIP7AdKAdOklB/c93gQsAQIB64AIVLKv9M5Tm+gN0DRokVfvHDhwhMHDgoKIjHxLBMm+FO3rqX/+8RfWwidPI5WQ4ZTpmbdJzpejx49WLhwIRcvXiRfvpzZAr517QpHNqzhWNgGkuNi8fIrRCXrBUtObu5ax1MUJRtl2igXIYQnsAwYKKU8dtf97oDZ2i0TDHwrpSz9qGM9bR/62bNnOXvuW+zsttGg/oE7I1t0ej1dv56Soe6Wux07doyKFSsyZswYRo4c+cR5spMxNZVTu//i8LpQrpw6gd7OjrK16lG5STB+pcuqC5YUJQ/I1GGLQoiPgAQp5fhH7HMeCJRS3njYPs8yDv3f018SHv4zDYNO8M9fW1g1eRwtBw+nbK0na53f1rx5cw4dOsSFCxdwcLCNaXEjL5zj8PrVHN+2mbTkJHwCilO5STDl6jbA3knN/a4oudUzDVsUQvhYW+YIIZyAxsA/9+3jK6zNQyFEdetxbz5j7nRt2LCBXTsvYjanYjKlsXPJb+TzL0qZGrWf+pghISFcv36dBQsWZGLSrOUTUJzGvd6h73c/07iXZaHsDT9MY2a/bmz4YTqRF8495giKouQ2GRm2WAn4GdBjKdQLpZSfCiH6AkgpvxNCDAD6AUYgCRgqpXzkAO9nmcslOfkKn39hws9+BqunTKbl4A8oW6veEx/rNiklVapUwWg0cuzYMZvsupBScvXfkxxeH8rJndswpaVRqEw5KjcNpkyNOhjs7bWOqChKJshVV4peunSJK1cWEp8wncvrGmJOs6fbuKlP3Hd+v3nz5tG1a1dCQ0Np0aLFMx1La0nxcfwdtoEjG1Zz6+oVHN3cqRDUmEqNm+PlW0jreIqiPINcVdABLl/+jX9OjuTv+aVo2mvknZkUn0VqairFixenXLlybNiw4ZmPlxNIs5mLfx/h8PpQTu/dhTSbCaj0ApWbtKDkizXUBUuKYoNy1QIXa9asITLiMIX8wdvflzI162TKce3t7Xn33Xf54IMPOHToEFWqVMmU42pJ6HQEVKxCQMUqxEfd5OjmdRzZuJYV33yBq5c3FRs1o+JLzXDL9/BpEhRFsR0200Lfs3wxviXL0HXAIGJvnmfct3qcY0PQCx+qt2mfKZmio6MpUqQIr7zyCvPmzcuUY+Y0ZpOJswf3cXh9KOcPH0AIQckXq1O5STABFas8c9eVoihZK1d0uVw8doQ/J31JrS69+GvVV/hXusLlsHK0GjKKohUqZVquwYMHM23aNM6dO4e/v3+mHTcnir5+jSMb13Bs0zqS4mLxLOhHpcbNeT6oMc7uHlrHUxQlHbmioIOlqC8f/xlSH4PZqKNJv248X7NzpuY6f/48JUuW5L333uPrr7/O1GPnVMa0NP7ds4PD60K5/M/f6A0GytSsS+UmwRQqW84mR/0oSm6Vawo6wKwPhxN78hgFq0bS/K0vyZ+/YaZn69ixI2vWrOHSpUu4u+etS+tvXLpguWBp6yZSkxLJX7QYlRu3oFy9hjg4qwuWFEVruWY+9IvHjhB54ij/xCVw47gXl/85nSWvExISQmxsLD/++GOWHD8ny18kgEY9+9L3u7k06T0QnV7PxtkzmNm3K+u/n0rE+bNaR1QU5SFspoV+uw+9Xs9+ePi7cGBLN8I3l6X10A8ztQ/9tvr163PhwgXOnDmDwWBzg4EyjZSS62f+5dD6UE7u2IYxNQW/UmWp3DSYsrXrqxWWFCWb5YoW+rUzp2g5eDgVa9ejQMFCuBVOpFbXOlw7cypLXi8kJISLFy+yePHiLDm+rRBC4FuqDM37DabPjJ9p2O1tUhITWDN9IvPeH0j4iWOPP4iiKNnCZlroty1duhSjMYF8+T+mdOlRFC3SIwvSgdlsply5cri5ubF37151YvAuUkrOHtjLpp9mEht5nUqNmlPvje44urhqHU1Rcr1c0UK/bfLkyUybNgsAsyn5MXs/PZ1Ox9ChQ9m/fz9bt27NstexRbfHrncfP40XW7bl6KZ1zBnaj1O7tqNVA0FRFBtsocfExCClZP+BQIoVe4eSJYZmQTqLpKQkihYtSq1atVixYkWWvY6tu372NOu+n0LEuTOUeLE6jXr2wz2/j9axFCVXylUtdA8PDzw9PdHpHLO0hQ7g5ORE//79WblyJf/888/jn5BHFSxRijc+n0CDLj25eOwwc957hwOrV2I2m7SOpih5is0V9N9//53ff/8dvd4Rkzkly1/vnXfewcHBIUetO5oT6fR6Alu1o/v4aRQuW47Nc2by6+hhal52RclGNlfQZ8yYwYwZM6wt9KQsf70CBQrQrVs3fv75ZyIiIrL89WydRwFf2o34hOCBIcREXGf+iMFs+/Vn0lKz/o+vouR1NlfQQ0NDCQ0NRa93wmTO+oIOMGTIEFJSUpg+fXq2vJ6tE0JQrm4QPSbMoFzdhuz5YxFzQwZw4eghraMpSq5mcwXd2dkZZ2dnS0HPhhY6wHPPPUerVq2YNm0aSUnZ85q5gZObO83fGcxroz8HAYvHjGLN9EkkxcVqHU1RcqWMrCnqKITYI4Q4LIT4WwjxSTr7CCHEZCHEaSHEESFE1ayJC/Pnz2f+/Pnodc6YTIlZ9TIPeO+997hx4wZz587NttfMLYpWqEzXcVOp0bYDJ7Zv5qchfTmxbbMa4qgomSwja4oKwEVKGS+EsAO2A+9KKXfdtU8wMBAIBmoA30opazzquM+ypijApEnFSUuLplq1ZU98jKchpaR69epcu3aN7du3ExAQkC2vm9tEXjjHuu+ncO30KYpVrkrjXu/gUcBX61iKYjOeadiitIi3btpZf+7/K9AGmGvddxfgKYTwe5bQD7N+/XrWr1+PTu+MMRtb6EIIZsyYQXx8PPXq1ePUqayZciC38wkoTufPxtGwex8unzzBnJD+7Fu5FLNJDXFUlGeVoT50IYReCHEIiADWSyl337dLYeDSXdvh1vvuP05vIcQ+IcS+yMjIpwpsZ2eHnZ0der0T5mws6ACBgYGEhYWRnJxM/fr1OXZMzWPyNHQ6PVVbtKL7N9MpWqEyW+bPZsH/hnL9bNbMnqkoeUWGCrqU0iSlrAL4A9WFEBXu2yW9iU4e6MuRUn4vpQyUUgb6+DzdlYRz5sxhzpw56PXO2TbK5W6VK1dm69at6PV6GjRowNMudK2Ae34fXhk2mlZDhpMQHcWC/w0lbO4PpCVn7QVjipJbPdEoFyllNBAGNL/voXCgyF3b/sCVZwn2MP8VdKdsPSl6t+eee45t27bh7u7OSy+9xPbt2zXJkRsIIShTsy7dJ8ygYqOm7F/1B3NC3uHcof1aR1MUm5ORUS4+QghP620noDFw/3XwK4Cu1tEuNYEYKeXVzA4LEBYWRlhYGHqdM2ZzClJq0/daokQJtm3bRqFChWjWrBkbNmzQJEdu4ejiSpO3B9Dx4y8x2NmzdOxHrJo8jsSYaK2jKYrNyEgL3Q/YLIQ4AuzF0of+pxCirxCir3WfUOAscBqYBbyTJWnvojdYlkPLrrHo6fH392fLli2UKlWKl19+mZUrV2qWJbfwL1eBN7+eQq32r3Nq11/8NKQvxzavV0McFSUDbG62xVmzLFPntmjhzMlTH1K3zi4cHLSd2S8qKooWLVpw4MAB5s+fT8eOHTXNk1vcDL/E+llTuPzPcYo8X4kmb/fHy++Bc+2KkqfkqtkW/5ucywlAs370u3l7e7N+/Xpq165N586dmT17ttaRcoV8/kXo+NGXNHl7ABHnzvDzsAHsXrYQk9GodTRFyZFsrqBv2LCBDRs2oNdbu1w0GOmSHnd3d1avXk3Tpk156623+Prrr1U3QSYQOh2VGjen+4QZlKxane2/zWX+8He5+u9JraMpSo5jcwX9tv9a6AkaJ/mPs7Mzy5cvp1OnTnzwwQf069cPo2pNZgpXL29aDR1Bm2GjSU5M4JfRIWz6aSapSdp/Q1OUnMLmlrO/PeNh59erAdqeFE2Pg4MDCxYsoHjx4owdO5YLFy6wcOFC3NzctI6WK5QKrEGR8hX56/d5HFz7J//u3Umjnv0oFfjImSYUJU+wuRb6ypUrWbly5Z0WenZfLZoROp2OL774glmzZrF+/Xrq1atHeHi41rFyDQdnZ17q0YfXPxuPo4sry8d9xooJXxB/K0rraIqiKZsr6KtXr2b16tXoddoPW3ycXr16ERoaytmzZ6lRowaHDh3SOlKu4le6LF3GTqJup66cPbCXOUP7cXj9aqTZrHU0RdGEzRT07du3c+7cf8uZ6fVOREcX5MCBS494lvaaNm3KX3/9hU6no169eqxevVrrSLmK3mCgRtsOdBs3lYIlSrLhh2n8/slwbobn7P8vFCUr2ExBL1y4MIsWLWL8+PGMHz+eS5ducuJEffLly/lvoWLFiuzevZvSpUvTqlUrvvvuO60j5TpefoVpP+pzmvUbzM3wS8x9fyA7Fi3AmJamdTRFyTY5vxpardI5U6JlW6Kioti1axdLlvyJvuI11jlkySy9ma5QoUJs3bqV5s2b069fP95//33MqmsgUwkhqBDUmB4TZlCmZh12Lv6Vee8PJPyEmhVTyRtspqBXcXfmw6hk8tVrSIUKFUgqUZbZbj0pa3dT62gZ5urqyh9//ME777zDuHHj6NSpk1rSLgs4e3jy8qBhtBvxCca0NH7/eDjrv59KckL845+sKDbMZoYt1vVy41NvR0ISPfCuVIcoV3f6JfxEZTvbaKHfZjAYmDp1KiVLliQkJITw8HCWL1/O004nrDxc8Sov0n38NP5atIADq5ZzZv9uXurRh9I16mBZiEtRchebaaGfO3eOs38uo9iVc1zx8qFk5BWMh725esX25s4WQjB06FAWL17MwYMHqVmzJidPqisfs4KdoyNBb77FG19MwMXTm5UTv2T5+DHE3ni6BVYUJSezmYJ++fJlSrRsy78+/iSHrePfAv4YKidw46btXonZrl07wsLCiIuLo1atWmzdulXrSLlWwRKleOOLCTTo0pMLRw8x5713OLB6JWazWvpOyT1spqDzfGU+jEpmTOUyOAY1pcrFU/zg8hoJAc5aJ3smNWrUYNeuXRQoUIAmTZrwyy+/aB0p19Lp9QS2akf38dMoXLYcm+fM5LfR7xN58bzW0RQlU9hMQT8Um8j3zxejpqcrAH5C8k7aXP5J89I42bMrUaIEO3bsoFatWrzxxht8/vnnamKvLORRwJd2Iz4heGAI0RHXmD/8Xbb9+jNpqSlaR1OUZ2IzBX1AQEHqerkxZ+I3JPwyG3t3d0qlnqOdYYvW0TKFt7c3a9eupUuXLowaNYpevXqRpsZQZxkhBOXqBtFjwgzK1W3Inj8WMXfYAC4eO6x1NEV5ahlZgq6IEGKzEOKEEOJvIcS76ewTJISIEUIcsv58mDVx4dTRI6SdOYnexZWUZJljps/NDA4ODsydO5cPP/yQ2bNnExwcTExMjNaxcjUnN3eavzOY10Z/DsCiz0ayZvokkuJiNU6mKE/usSsWCSH8AD8p5QEhhBuwH3hFSnn8rn2CgBApZcuMvvDTrliUYDJRcutR3tCl8OKVURQtmkzdOtue+Dg53c8//0yvXr0oW7Ysq1atIiAgQOtIuV5aagq7lvzGvpVLcXBxpWG3t3muTgM1xFHJUZ5pxSIp5VUp5QHr7TjgBKDZOmBOOktk4eiEyWTAaMx5sy1mhm7durFmzRrCw8OpWbMm+/fv1zpSrmdn70C9zt3oMnYSHgUKEjplPEvHfkRMxDWtoylKhjxRH7oQohjwArA7nYdrCSEOCyFWCyGef8jzewsh9gkh9kVGPt044M/HjCFl/ixwcMBsMmDOwbMtPqtGjRqxY8cOHBwcqF+/vlqEOpv4BBSn82fjaNi9D5dPnmBOSH/2rVyK2aSGOCo5W4YLuhDCFVgCDJZS3t/BeAAIkFJWBqYAf6R3DCnl91LKQCll4NNeGXny5EnkpQuY7ewxmQ1IUpAy9/6ilS9fnl27dlG+fHleeeUVpk6dqnWkPEGn01O1RSu6fzOdohUqs2X+bBaMHMr1s6e1jqYoD5Whgi6EsMNSzBdIKZfe/7iUMlZKGW+9HQrYCSHyZ2pSq/nz51Pmk69JFTp0whHI2XOiZwZfX1/CwsJo1aoVAwcOZMiQIZhUazFbuOf34ZVho2k1ZDgJt6JY8L+hhM37kbRk27tCWcn9MjLKRQA/AieklBMeso+vdT+EENWtx82yWbOc9ToSTGYcHD2A3F/QAVxcXFiyZAnvvvsukyZNon379iQm5s7zBzmNEIIyNevSfcIMKjZqyv4/lzEnpD/nD6nzGkrOkpEWeh3gTeClu4YlBgsh+goh+lr3aQ8cE0IcBiYDnWQWXRnz4YcfEj5rCokmM85Otwt63ihser2eSZMm8e2337J8+XKCgoK4fv261rHyDEcXV5q8PYCOH3+Jwc6OJWM/YtXkcSTGRGsdTVGADMy2KKXcDjxy3JaUciqQLZ27ly5dwngz1lLQnb0BMBoTsuOlc4xBgwZRrFgxOnfuTM2aNQkNDaVcuXJax8oz/MtV4M2vp7Dnj0XsXraQ84cP0ODNt3i+QSM1xFHRlM1cKXrbTz/9RP0x40g0m3B1ywdAbOwNjVNlv9atW7NlyxaSkpKoXbs2mzdv1jpSnmKws6P2a6/T9esp5PMvwtoZk1g8ZiS3rl3ROpqSh9lcQQdw0etINJlxd7Ocd42JjdA4kTYCAwPZvXs3hQsXplmzZsydO1frSHlOPv8idPzoS5q8PYDrZ88wN2QAu5ctxGS03VlAFdtlcwV9xIgR7P92HIkmMx4eBQCIi7OdVYsyW0BAANu3b6devXp069aNjz/+WE3slc2ETkelxs3pPmEGJapWY/tvc5k/YjBX/1Vz3CvZy+YK+s2bNzHGRFsLumUse0J83i3oAJ6enqxevZru3bvzySef0K1bN1JTU7WOlee4ennTaugI2gwbTXJ8HL+MDmHTTzNJTcobJ+0V7dnMEnS3ff/994w9e5UpF65jMLgAkJB4S+NU2rO3t2f27NmULFmS0aNHc+nSJZYuXYqXl+1PL2xrSgXWoEj5ivz1+zwOrv2Tf/fupPFb/Sj5Yg2toym5nM210AGcdTrMgFHnBEByUrSmeXIKIQSjRo1i/vz57Nixg9q1a3Pu3DmtY+VJDs7OvNSjD50/HYejswt/fP0ZKyeMJf5WlNbRlFzM5gp6SEgIq8Z+AkAKDgAkp8RpGSnHeeONN1i3bh3Xr1+nZs2a7NmzR+tIeVahMs/R5ctvqdupK2cO7GHO0H4c2bAGaTZrHU3JhWyuoCclJSFTLCvLJEsDoMdsSiRZXYp9jwYNGrBz505cXFwICgpi2bJlWkfKs/QGAzXadqDbuKkULFGS9bOm8vsnw7kZfknraEouY3MFfdq0afT/ajwAiSYzQjii0xuJjo7WNlgOVLZsWXbt2kWlSpV49dVXmThxohoBoyEvv8K0H/U5zfoN5uali8z7YCA7Fv2CUa1MpWQSmyvoYJnLBSwFXa93Rq83cuuWOjGangIFCrB582batWvH0KFDGThwIEY1RlozQggqBDWmx8TvKF2jDjsX/8K89wcS/s/fWkdTcgGbK+iDBw/mh1EjAMvqRQaDC3qdKuiP4uTkxMKFCwkJCWHatGm0bduW+Ph4rWPlac4enrw8aBjtRnyCMS2N3z/6gPWzppKcoP67KE/P5go6gME6X0aiyYydwQU7O6kK+mPodDrGjRvH9OnTCQ0NpX79+ly5oi5T11rxKi/Sffw0XmzZlqMb1zFnaD9O7dquusaUp2JzBX3SpEl8ON7ah242o9M7Ye+AKugZ1K9fP1auXMmpU6eoWbMmR48e1TpSnmfn6EjQm2/xxhcTcPH0ZuXEL1k+fgyxN55uVS8l77K5gg6WcegACSYzer0T9nZmVdCfQHBwMNu2bcNoNFK3bl3Wr1+vdSQFKFiiFG98MYEGXXpy4egh5rz3DgfXrMRsVouZKBljcwW9f//+fDJkMHDXSVGDiVu3bmFWY3sz7IUXXmD37t0EBAQQHBzMjz/+qHUkBdDp9QS2akf38dMoXLYcm36ayW+j3yfy4nmtoyk2wOYKupOTE27OzgAkWVvoOp0Rs9lMbOz9S50qj1KkSBG2b99Oo0aN6NWrFyNHjlR/FHMIjwK+tBvxCcEDQ4iOuMb84e+y7defSUtN0TqakoNlZAm6IkKIzUKIE0KIv4UQ76azjxBCTBZCnBZCHBFCVM2auDB+/HgmfjMewe0uF2eEsExEpbpdnpy7uzsrV67k7bff5osvvqBLly6kpKiikRMIIShXN4geE2ZQrm5D9vyxiLnDBnDx2GGtoyk5VEZa6EbgPSllOaAm0F8IUf6+fVoApa0/vYEZmZryPkKIO3Oi6/XOSGm5SlQV9KdjZ2fHzJkz+fLLL/n1119p0qQJN2/m7RkscxInN3eavzOY10Z/DsCiz0ay9rtvMRnVBUnKvR5b0KWUV6WUB6y344ATQOH7dmsDzJUWuwBPIYRfpqcFevfuTe/eva0LRZvQ65yRMgUhVEF/FkIIPvjgA3777Tf27NlDrVq1OH36tNaxlLsUrVCZruOmUr1Ne45tXs+mn2ZqHUnJYZ6oD10IUQx4Adh930OFgbsnpgjnwaKPEKK3EGKfEGJfZOTTDcnKly8f+fLlw0WvuzPKBcDb20UV9EzQsWNHNm7cSFRUFLVq1WLHjh1aR1LuYmfvQL3Xu1O9TXuObFjDoXWhWkdScpAMF3QhhCuwBBgspbz/7GN6K+M+cGWElPJ7KWWglDLQx8fnyZJajR07lrFjx+Jm0BNrNKHXW06Qenmpgp5Z6tSpw86dO/H09OSll15i0aJFWkdS7lOn05uUqFqNzXNmcum4upZAschQQRdC2GEp5guklEvT2SUcKHLXtj+QpZchuuv1xBn/a6F7ejqpgp6JSpcuzc6dOwkMDKRDhw58/fXX6urFHESn0xM8MATPgn6snDCWmIjrWkdScoCMjHIRwI/ACSnlhIfstgLoah3tUhOIkVJezcScd/To0YMePXrgbtATazKhs7bQ3d0dSExMVCM0MlH+/PnZsGEDnTp14oMPPqBfv35qYq8cxMHZhTbDRmM2mVg+7jPS1BTSeV5GWuh1gDeBl4QQh6w/wUKIvkKIvtZ9QoGzwGlgFvBO1sS1jJ0uUqQI7gY9cUbTnRa6m5tlsQvVSs9cjo6OLFiwgBEjRjBz5kxatWpFXJxaUCSn8C5UmJbvvs+NSxdZM11Nj5zXPXZNUSnldtLvI797Hwn0z6xQj/Lpp58C8OG/l4kxWka5ALi62gGWgu7r65sdUfIMnU7HF198QfHixenXrx/16tXjzz//xN/fX+toClCsyovUe6M7W+fPZvfS36n5aietIykasbkrRW9zM1hGuWBtoTs76wHVQs9Kb7/9NqtWreLs2bPUqFFDTeyVgwS2bEu5eg35a+F8Tu/dpXUcRSM2V9C7dOlCly5dcDdYCniStLTQ9XoTjo6OqqBnsWbNmrF9+3aEEDRs2JBDhw5pHUnBch1Bk94D8C1ZmtCp33Dj0gWtIykasLmCXrZsWcqWLYubtaAnSEvfucmUiJeXlyro2aBSpUps2bIFZ2dnGjVqxMGDB7WOpGAZo946ZCT2jo4sHzeGpHh1riOvsZmC/vXXsHkzjB49mtGjR+Nh0JN60JVp33oCYDInqYKejUqWLElYWBiurq40atSI/fv3ax1JAdy889P6vZHE3Yzkz0lfYTapqXfzEpsp6NWqQYcOlqJ+6RL8+5c90Z8Wo2QVy1n92y306OhoNWNgNilRogRbtmzB3d2dxo0bs2/fPq0jKUChMs/RuFd/Lh49xJb5s7WOo2QjmynoDRvCwoXQtGksRYsu4aMeTnh+eJ4ydcwIYcBksrTQTSaTGlaXjYoVK8aWLVvw9PSkcePG7NmzR+tIClChYROqtmjNgdDlHAvboHUcJZvYTEEHS1EPDDwKvErlQDP2L8RbL/93wmRKwMvLC4CoqChtg+YxAQEBbNmyhXz58tGkSRN2775/qh9FCw3efIuiFauwYdZUrpw6oXUcJRvYVEHfvBlOn66Dlxfs2qIj9aCrpaDrnDGZkvD29gbU0EUtFC1alLCwMHx8fGjSpAk7d+7UOlKep9PraTn4A9zy+bDimy+Ii7qhdSQli9lMQd+82dKHvnAhdOkCOh1Ef1qMA9v06PROmEyJuLu7I4RQBV0jRYoUISwsDF9fX5o2bcpff/2ldaQ8z8nVjTbDRpGanMyK8Z+rFY9yOZsp6Hv3Wor51Kmvsn//q6SkCDw6RnDqoB6D3gWzKQm9Xo+np6cq6Bry9/dn8+bNFCpUiObNm7N9+3atI+V5+YsEEDzgPa6d+Zf1309V0wPkYjZT0N9/39KHXqtWLVq2rIWLC3DZico9Y++00AE1dDEHKFy4MGFhYRQuXJjmzZuzdetWrSPleaWq1aROhy6c2LaZ/X8u0zqOkkVspqB/t+UMO87cICQkhBEjQujaFZKPm9izPcJyUtScBIC3t7c6KZoD+Pn5ERYWRpEiRWjRogVhYWFaR8rzarTrSJmaddm6YA7nD6nrBnIjmynolfw9GPDLQXacvsGN+BSavHmD/C0PEL4lP3q98z0t9KSkJJKSkjROrPj6+hIWFkaxYsUIDg5m06ZNWkfK04QQNO83mPxFA/jz26+JunJZ60hKJrOZgl67ZH6mvv4Cb8zYTLk3PuTjdQewjy3L+V+Lcu1aIUym/1rooIYu5hQFCxZk8+bNlChRgpYtW7Jx40atI+Vpdo6OvDJsNDq9nuXjPiMlMUHrSEomspmCDpai/pxjHC7l6lO/TH6qvylAL5k1q/mdFroaupjzFChQgM2bN1OqVClatmzJ+vXrtY6Up7n7FKDV0BFEX79K6JTxmM1qeoDcwmYK+oG1F1i7+TxXHSxzcG88EYF35E1adTjJsqVVuXTJUsjVxUU5k4+PD5s2baJMmTK0atWKtWvXah0pTytSviINu/fh7IG9/PXbPK3jKJnEZgp6pL3k6KIzfFW3DO6OBprl96T4lihuvBiLwWBm3rxeSGnG3t4eV1dXVdBzoPz587Np0ybKlStHmzZtWL16tdaR8rQqTYOp1Lg5e5Yv5sRfW7SOo2SCjKwpOlsIESGEOPaQx4OEEDF3LU/3YebHhHM6ExVfK8nfvxyj8o6tBJxMJKpxAc67p9L59ROsX9eWkyctayp6e3urLpccKl++fGzcuJHy5cvzyiuvsGrVKq0j5Wkv9ehD4eeeZ913k7l+9rTWcZRnlJEW+hyg+WP22SalrGL9+fTZYz2ob4OSNGtYDJPrLWqWaMQFJyhRwxdTcTe69juNwWDk8y8sb0cNXczZvL292bhxIxUrVqRt27asXLlS60h5lt5gR+uhI3Byd+eP8WNIiFYNIVv22IIupdwK5IjqGH7yFp4iAIAi8RLDecsZepcCOlq1WsAvCxz4919LP3pcXBxpaWlaxlUewcvLiw0bNlClShVeffVVli9frnWkPMvZw5NXho0mOS6OFRPGYlS/NzYrs/rQawkhDgshVgshnn/YTkKI3kKIfUKIfZGRkU/0AuEnb7F21jGCXi8LwBF7IwlLLhBwPY1EnOnUaSYODpLPPvvvxKjqdsnZPD09WbduHS+88ALt27dn2TJ1BaNWChQrQfN3BnPl5HE2zZ6hpgewUZlR0A8AAVLKysAU4I+H7Sil/F5KGSilDPTx8XmiF4k4H0uztyvQd/gbTPlzGKkCjI0KUijKSCJOeHnf5K23IlmwAG7dKgCogm4Lbhf1wMBAOnTowJIlS7SOlGeVrVWPGm07cHTTOg6tU+c2bNEzF3QpZayUMt56OxSwE0Lkf+Zk96naLAD/sl506tSR6uVewkUKYtz17CznRILZsq7ogAEXcHSEH37IB6ihi7bCw8ODtWvXUq1aNTp27MiiRYu0jpRn1enQhRIvVmfznO+5eOyI1nGUJ/TMBV0I4SuEENbb1a3HvPmsx32Yt99+m+D67fEUOm7esoxqSZD2AHh7x9K5M/zxhwFwJjo6OqtiKJnM3d2dtWvXUrNmTTp37szvv/+udaQ8Seh0BA8IwcuvMCsnfUlMxDWtIylPICPDFn8FdgJlhRDhQoi3hBB9hRB9rbu0B44JIQ4Dk4FOMos74Jw9HPDU6bl+y3K5f7y1oJtMiXTuDPHxgvDwSqrLxca4ubmxevVqateuzeuvv86vv/6qdaQ8ycHZmVeGjUKaTfwxbgypyWpeJFuRkVEunaWUflJKOymlv5TyRynld1LK76yPT5VSPi+lrCylrCml3JGVgYOCghg9vTcuZrgcaRnlkmA2AJaCHhQEvr5w5MjzqqDbIDc3N0JDQ6lbty5dunRhwYIFWkfKk7z8CtPy3Q+4eekia6ZNRKqF122CzVwpelv37t1p1aQ9hlTJ5ZuJOOkE8WY9ACZzEno9dOwIBw8W5urVJHW23ga5uroSGhpK/fr16dq1K/PmqUvTtVCsclUavNmTf/fsYNdS1QVmC2yyoL/WtjNCgoMZXHQ64s2Wt3F7gq4uXcBo1DF79utMn55CfLyWiZWn4eLiwqpVqwgKCqJbt278/PPPWkfKk6oGt6F8/ZfYsWgB/+7J0i/fSiawuYKelpaGvYsltqsUOCCINUmEMNyZQjcwED799BpGo4EBAxzx84NeveCzz+Dll2GLmrbCJjg7O7Ny5UoaNWpEjx49+Omnn7SOlOcIIWjy9gD8SpVl9dQJRF48r3Uk5RFsrqA3adKEtwZ1BMDVLNCZJHFGk2XVImsLHaBPHx39+k1n7tzTdOgAv/4KH34IGzZA48bw3XdavQPlSTg7O7NixQrLf/e33uLHH3/UOlKeY7C3p/V7/8Pe2Znl4z4jKS5W60jKQ9hcQe/Vqxc9erwFgI/OAGlmYo1m9DpnzKb/zsZ7enoiBAQEXOHHH+H6dbh2DSIioGlT6NcP3nkH1FXOOZ+TkxPLly+nWbNm9OrVi1mzZmkdKc9x9c5Hm/dGEn8rij8nfYnJaNQ6kpIOmyvoXbp04a3e3dHpBP4OdqQlG4kzmu5ZKBrA3t4eFxeXOyNdXF2hYEHw8IAVK+CDD2DGDGjSBJ5wFgJFA46Ojixbtozg4GB69+7NzJkztY6U5/iVLkuTtwdw8dgRtsxT35RyIpsr6ImJiSQnJ+Hq7UB+9CQnGYk1mSzriprvHS/r5eWV7tBFvR6+/BLmz4fdu6FaNTiiLorL8RwdHVm6dCkvv/wyffv2Zfr06VpHynOeb9CIF19+hYNrVnJ00zqt4yj3sbmCHhwcTHBwMAUC3HFLMJOQkEZsmrUP3Xjv+oheXl6PvFr0jTdg61ZLt0vt2rB0aRaHV56Zg4MDS5YsoVWrVvTv35+pU6dqHSnPqf9GDwIqvcCGH6Zz+eQJreMod7G5gt6vXz/69euHbwkPdEkmDKkmUqTEpHN7oIXu6elJTEwMJtPD10ysVg327YMKFeDVV+GTT0BdQ5GzOTg4sHjxYtq0acPAgQOZPHmy1pHyFJ1eT8t3P8Ddx4cV33xO7A3VZ5lT2FxB79ixIx07dqRgCXcAvFIs9yfrPO7pQwdLC11KSUxMzCOP6ecHYWHQrRt8/DF06IAau57D2dvbs3DhQtq2bcu7777LxIkTtY6Upzi6uvLKsNEYU1NY8c3npKWmaB1JwQYLekxMDDExMfgUcUNnEOS/XdCF+51x6Lfdnhc9I5N0OTrCTz/BhAmwbBnUqQPnz2dyeCVT2dvb8/vvv/Pqq68ydOhQvvnmG60j5Sn5/IsSPDCE6+fOsO67yeqq7BzA5gp6mzZtaNOmDXqDjgJF3SiYKgBIEm7pttAh4/OiCwFDhkBoKFy8aOmO2bo1c/MrmcvOzo5ff/2V1157jZCQEMaNG6d1pDyl5Is1qNvxTf75awt7V6i57LVm0DrAkxo0aNCd274lPCh42FLEk3B9oIXu7u6OTqd74km6mjWzjH5p3RoaNYKpU6FPn2fPrmQNOzs7fvnlF/R6Pe+//z4mk4nhw4drHSvPqP7Ka0ReOMe2X38mf9EASrxQTetIeZbttNC3T4JzW2nXrh3t2rUDoJj7CZqYFwOQiAtmcxJS/ndGU6fT4eHh8VTzopcpYynqTZpA377qIqSczmAwMG/ePF5//XVGjBjBF198oXWkPEMIQbO+7+ITUJzQyeOJuhKudaQ8y3YKeuGqsKg7MYdWcuPGDTi3lUIHBxOZWgqABOkEgNmcfM/THjYWPSM8PGDlSnj/fctFSE2bwo0bz/Y2lKxjMBiYO3cuXbp0YeTIkXz22WdaR8oz7BwdeWXYKHQGA398/RnJCWpUgRZsp6AXrw+vzUEu7MaqARVgYTdEhzlE6qsAcDPtv0Uu7ubp6flM86Lr9fDVVzBvHuzcaZn4a/fupz6cksX0ej1z5syha9eufPjhh3zyySdaR8oz3PMXoPXQEcREXCN08jjM5ocPF1ayRkZWLJothIgQQhx7yONCCDFZCHFaCHFECFE182NaFa9PhNeLdHsuCcq2gOL1KeDrCsCNlP8Wubibl5cXiYmJpKQ827CqLl1g+3bL7Xr1YNIkUCf1cya9Xs/s2bPp3r07H3/8MR999JEagZFN/MtV4KUefTl3aD/bf52rdZw8JyMt9DlA80c83gIobf3pDcx49lgPcW4rZcz/Wm7//Qec20qxMp7Yp0luJFkXuXiGoYuPExgIBw9CcLBlNEy7dqAWRcqZ9Ho9P/74Iz179uTTTz/lww8/VEU9m1Ru0oLKTYLZu2IJJ7Zt1jpOnpKRJei2AlGP2KUNMFda7AI8hRB+mRXwjnNbYVF3ohpN5JrZC4rWhEXdKed5Esc0M9FJ9y5ycZunpyeQ8aGLj+PlZRmnPmEC/PknVK0Ke/dmyqGVTKbT6Zg1axa9evVizJgxjBw5UhX1bNKwe2/8y1dg3cwpXDvzr9Zx8ozM6EMvDFy6azvcet8DhBC9hRD7hBD7Ip90isPLB+C1ObQL+ZZOixIg9gq8NgdP4z84pEoS02630J9tLHpG3B6vvm2bZZqAOnVg8mTVBZMT6XQ6Zs6cSZ8+fRg7diwjRoxQRT0b6A0GWg0ZgbOnJ8vHjyEhWn2VzQ6ZUdBFOvel+xsjpfxeShkopQz08fF5slepOxiK12f48OEM7xYMN05C4RfR1R+CnRlS5O0+9Hu7XJydnbG3t8+ULpf71axp6YJp3hzefRfat4cseBnlGel0OqZPn06/fv346quveP/991VRzwbO7h60CRlFckI8y7/5HKMa95vlMqOghwNF7tr2B65kwnHT1bx5c5q/0hGkGa4dBcBRL0jTpX9SVAjxTEMXH8fbG5Yvh/HjLfOsV61qmexLyVl0Oh3Tpk2jf//+jB8/npCQEFXUs0GBYiVo8c4Qrp76h40/TlefeRbLjIK+AuhqHe1SE4iRUl7NhOOm69KlS1zC2kV/aQ8A7g52pBisXS73zbgIzz508XGEgPfes0wTYDRaumB++CHLXk55SkIIpkyZwqBBg5gwYQJDhgxRBSYblKlZl5qvduLY5vUcXPOn1nFytcde+i+E+BUIAvILIcKBjwA7ACnld0AoEAycBhKBHlkVFuDNN98EIKxdCbiwA+oMooCHA6diLV/n7m+hg6Uf/ezZs0gpESK9HqLMUauWpQvm9dfh7bfh2DFLy91gcxMs5F5CCCZNmoROp2PSpEmYzWa+/fbbLP3/QoHa7V8n8sJ5wubOIp9/EQIqVtE6Uq702FIjpez8mMcl0D/TEj3GqFGjLDeSVsKxpWAy4u/lRHJSIhLuWVf0Ni8vL9LS0khISMDV1TVL8+XLB6tWwbBhlrHqx4/D779bRscoOYMQggkTJqDT6ZgwYQJms5kpU6aoop6FhE5H8ICh/DIqhD8nfcUbX0zEs6Cv1rFyHdu5UtSqcePGNG7cGIo3gJRYuHIQXyd7zDpBqtkRYzot9Mweuvg4BgNMnAg//miZZ71GDfjnn2x5aSWDhBCMHz+eYcOG3elbN6uVTbKUvZMzrwwbDVLyx9efkpr04O+q8mxsrqCfPXuWs2fPWqYCADgXhrud5YtGvNGT5IS4B56TFUMXM6JnT9i0yTLypWZNWLMmW19eeQwhBF999RUffPABM2bMoF+/fqqoZzFPXz9aDh5O1JVwVk+bgFSfd6ayuYLes2dPevbsCS75oWBFOLsFd+sJ0XiTJ/HprE50u4WeFUMXH6duXcuFR8WKwcsvWy5IUufhcg4hBGPHjuV///sf33//PX369FFFPYsFVKpC0JtvcXrvLnYs/lXrOLmKzZ2uu2eypRINYM8s3DACkCDdSYyPfeA59vb2uLi4ZHsL/baAAMs8MN26WUbDHD0K330HDg6axFHuI4RgzJgx6HQ6xowZg9lsZtasWeh0NtfesRkvtGhNxIVz7FryKz4BxShTo47WkXIFmyvoDRo0+G+jeAPYORWPmycAd2LN7iQnpz9LQVaORc8IV1dYtMiyCPWnn8KpU7B0KRQsqFkk5S5CCD799FP0ej2ffPIJZrOZH374Ab1er3W0XEkIQeNe/Ym6fInV0ybg5VsIn4DiWseyeTbXBDl58iQnT560bATUBp0BtyuWyVRipCtGYyLGtAen7fTy8iIq6lFT0mQ9nc5S0BcutAxvrFbN8q+SMwgh+Pjjj/n444+ZM2cOPXr0wGRSU8BmFYOdHa3fG4mjiyt/jBtDYuyjF3NXHs/mCnqfPn3oc3s9OAdXKByI+0XLvLYxuKLTpxB58cHJ9f38/IiNjSUu7sGTptnttdcsXTBSWvrYFy/WOpFyt48++ojPPvuMefPm0a1bN1XUs5Crlzdt3htJQnQUf078EpPRqHUkm2ZzBf2LL764d3mxEg1wv2K5YjRWuCL0qVw7++Bfen9/fwAuX76cLTkf5/YsjZUrWwr8xx9bJvpScoZRo0bx+eefs2DBAt58802MqtBkGd9SZWjaZxCXjh8lbO4srePYNJvrQ69du/a9dxRvgMuWrxFIUvWuCJFCxPkHT4z6+voihODy5cs899xz2ZT20Xx9YfNmywLUn3xiubL055/BxUXrZArA//73P/R6PcOHD8dsNjN//nwM6rLfLFG+XkMiL5xj38ql+AQUp1KjRy3BoDyMzf3feeyYZeGkChUqWO7wr4bOzgl3mYbR4IIwpxAT+eDVovb29hQsWDDHtNBvc3CAn36CSpUsV5eeOWOZ7KtoUa2TKQAffPABOp2O999/H7PZzIIFC7Czs9M6Vq5U7/Vu3Lh0gY0/fod34SL4P/e81pFsjs11uQwYMIABAwb8d4fBHorWwi0tFqPBBZ0+/YIOULhwYS5fvpzjxhkLAUOHWhbMOHvWsjLSX39pnUq5bdiwYXzzzTcsWrSIzp07k6amgc0SOp2elwcOw6NAAVZOGEvsjQitI9kcmyvo48aNY9y4cffeWaIB7qnRmHQOGAyppCalkpzw4C+dv78/KSkp3Lx5M5vSPpkWLWDXLvDwgIYNYfZsrRMptw0dOpSJEyeyZMkSOnbsSGpqqtaRciVHV1faDBuNMTWV5eM+Jy0lWetINsXmCnq1atWoVq3avXcWb4C7MQGjtEyuJPRpxN54sJVerFgxAA4fPpzVMZ9auXKwZw80aABvvWVpuavzcTnD4MGD+fbbb1m2bBkdOnRQRT2L5CtchJcHDSPiwlnWfjdZTXH8BGymoF/8+iK3Nt/i0KFDHDp0CIBbm29x8euL4FsJd5lCitlS0HWG9LtdvLy8KF++PHv37iUpKf1umZzAywtWr4ZBgyyTfLVsqVZCyikGDRrE1KlTWb58Oe3btyclJUXrSLlSiarVqNupKyd3bGXPcjWuN6NspqC7VXPjeIfjTOw6kQHdBhC1MYrjHY7jVs0NdDrcnVxINFuu6tMZUtJtoQPUq1ePlJQU9ubwlZ0NBvj2W5g1yzLBV40alqtLFe3179+f6dOns3LlSl599VVV1LNI9Tbtea5OA7b/Npcz+/doHccm2ExB92roRfmF5el+pjtDjwzlePvjlF9YHq+GlpkU3VzzEy+cADA5pBH7kBOjfn5+lCpVil27dtnEV+ZevWDjRoiKgurVYd06rRMpAP369eO7775j1apVtGvXjuRk1deb2YQQNO0zkALFShA6ZRw3wy89/kl5XIYKuhCiuRDipBDitBBieDqPBwkhYoQQh6w/H2Z+VEtR9+/tjzfeuL7oeqeYA7h7+ROvd0ECyfbJxNx4+C9YvXr1SExM5MCBA1kRM9PVq2e5CCkgwHLidNIkNWNjTtCnTx++//57QkNDadu2rSrqWcDOwZE2IaMw2DuwfPxnJMc/eBW48p/HFnQhhB6YBrQAygOdhRDl09l1m5SyivXn00zOyc0ffuDqtB1cnnOZaM9oYrbGcHXaDm5aF+90c8uPSehJwZF4fdJDu1wAAgICKFq0KDt27LCZKwCLFbMMZWzdGoYMsbTc1Td97b399tv8+OOPrF27ljZt2uToczO2yj2/D62H/o+YiAhWTf4as1lNxfAwGWmhVwdOSynPSilTgd+ANlkb60GpKeX5d/AtNvj9wFinz9Hrz3Pq3ZtcnbaDMy1bkjrzOwAScaZo7Enio5Ixmx/ejK1Xrx6xsbEcOXIku97CM3N1hSVLYNQoy5DGRo0gQg3V1VzPnj2ZPXs269evp3Xr1iQmqpV4Mlvh58rT6K2+nD98gG2//Kx1nBwrIwW9MHB351W49b771RJCHBZCrBZCZPolXmkJ/pSe5MXrjicY52KPV+GpeFY6gHSpikOJknh6eACWgl4wNRwp4cJn4x+6IkqpUqXw9fVl+/btOe5Co0fR6eCzz+C33+DAAcuMjdZBP4qGunfvzpw5c9i4cSOtWrVSRT0LVGrUnCrNXmbfyqUc37pJ6zg5UkYKenor597f9D0ABEgpKwNTgD/SPZAQvYUQ+4QQ+yIjI58oaNH3i+LXvzZeLVvio48iKa4+ro26UHH7APwnf0uxXj0BSMKZK8WcAbi+ajNXR4xAptOtIoSgXr16REVFcfz48SfKkhN07AjbtoHJBHXqWFruira6du3K3LlzCQsL4+WXXyYhIUHrSLlOUNe3KfJ8JdZ9P4Wrp09qHSfHyUhBDweK3LXtD1y5ewcpZayUMt56OxSwE0Lkv/9AUsrvpZSBUspAHx+fJw6bsGs3UaGhJLVqiXO+bdxavu3ORQe3l6FLxJnijucBsG/7BjHLVxA+eDDmdDqcy5UrR758+di2bZtNXrzw4ouWk6UVK0L79paFM2zwbeQqXbp0Yd68eWzdupXg4GDi1Um8TKU3GGg5+ANcPL1ZMf5z4m9pu8ZBTpORgr4XKC2EKC6EsAc6ASvu3kEI4SuEENbb1a3HzdTr6xN27ebykCF8a2fH0L17cW33EU7madz8bQsAbncKugtF7S4ikRgr1KLg6FHEb9jIpT59Md/XYtLpdNStW5fr169z+vTpzIybbfz8ICwMunaFjz6ytNxVw1Bbr7/+OgsWLOCvv/4iODg4R8zBn5s4u3vwyrBRpCQmsmL85xhtYPhxdnlsQZdSGoEBwFrgBLBQSvm3EKKvEKKvdbf2wDEhxGFgMtBJZnKTN/nYUQpPnMh7c35i5syZ+A5qQvS1XtxaZrngwMNa0JNwJc0gcLSLJvJSPN5vvEGhr74kce9eznd5k+ST935Nq1SpEh4eHmzbti0z42YrR0eYMwfGjbMsllGvHly8qHWqvK1Tp0788ssv7NixgxYtWqiinsl8AorTov9Qrp4+yYYfptvkN+yskKFx6FLKUCllGSllSSnl59b7vpNSfme9PVVK+byUsrKUsqaUckdmB83XqxcuNWtQtmxZypYti30Be1yq1+DWsSAA3AyWt5Kk8+CaUwD++mPcuGBZ6MKjTRuKTJ+GMSKCc6+2J3LqNKR1xjy9Xk/t2rW5ePEiFy5cyOzY2UYICAmxzNh45ozlZOmOTP+voDyJDh068Ntvv7F7926aNWtGbOyD8/QrT690jdrUat+Zv7ds4EDoisc/IQ+wmStFZx+bzZ6re9iyZQtbtli6Wa50uMKKYitI+DsBZ50OvYAk4U6MRwkK2p0mPiaNxFjL1zHXBg0o8edK3Fu04MbUqZx/owsp584BULVqVVxcXGy6lX5bcLBlxkY3N8uMjXPmaJ0ob2vfvj2///47e/fupVmzZsTEqHUzM1OtVztTqlottsz7kQtHDmkdR3M2U9Ar5KtAyJYQ/jfjf3z00UfsubqHse5jKXq+KJGLIxFC4K7XkyzcsHN2JElvORl149J/X3UNXl4UHvc1hSd8Q+qFC5xr9yq3fvsNg8FAzZo1OX36NFeuXHlYBJtxe8bGevWgRw947z3LaBhFG+3atWPRokXs37+fpk2bEq1mWss0QqejxYCh5PMvwp+TvuTWNdv//X0WNlPQq/tVZ3yD8ZhbmQkYFMB7W97jm4bfUM2nGpGLLUMg3Qx6koQ7ToYE9uvyARBx5sErb9yDgymxYjnOL7zAtY8/4VKfPrxQvDgODg5s3749W99XVvH2tszYOGAATJhgmbFRNQ6188orr7B48WIOHjxIkyZNuHXrltaRcg17RyfaDBsNQrB83BhS8vA1ADZT0MFS1GsVrsWBuAM0LtqY6n7V8WnvQ8KxBBL+ScDDWtDtRBQbZEXc9Ve58c/5dI9lV7AgRX6YRcGRI0ncvYcr7V+jio8Px48f50nHyOdUdnYwZQrMnAkbNlhmbPz3X61T5V2tW7dm6dKlHDlyhCZNmhAVpYbcZRbPgr60GjKcqCvhrJ72zUMvKMztbKqg77m6h20XtpEakcqa82vYc3UPPq9axrNHLo60ttBdMBtvkOJTGQ+7cCKvPHxIk9Dp8H6zC8WXLsGuUCF8p0xFLyXbw8Ky6R1lj969LQX9xg3LjI3r12udKO9q2bIly5Yt4+jRozRu3FgV9UxUtEJlgrq+zZl9u9mxaIHWcTRhMwV9z9U9hGwJgT8h/MdwOpbtSMiWEA6Lw7jXcSdyUSTuBh2J0om0tCg6VS/MTV0qscluJMc++iuYQ8mSFPvtVwr37EGJf//lyNGjXN26NZveWfZo0MByEVKRIpYZGydPVhchaSU4OJjly5dz/PhxGjVqlGOXRLRFLzRvSYWGTdm19HdO7swd3adPwmYK+rGbxxjfYDxzv5yLf29/vBy9GN9gPMduHqPAawVIOJKAX7iZBGkPwMvPO3LczhOAG/t3P/b4wt6eAu++S6MBAxDAxpkzifjmG8y56KKF4sUtMza2bAnvvmtpueeit2dTmjdvzvLly/nnn3946aWXuHHjhtaRcgUhBI3e6kehMuVYM2MiEefPah0pW9lMQe9ZoSfV/arzXInncPVx5WbyTar7VadnhZ7kb2eZZaDkhkTizQYAHPS3MDwfCMC1Iycy/Dq+tWtTuXJlzpUqxeW58zjfoSPJuWipIDc3WLoURo6EH36Axo0hl5wysDnNmjVjxYoVnDp1ipdeeinXnLvRmsHOjtbv/Q9HVzeWjx9DYmzeGQ1gMwX9trVr12L+x8zNpP++pjoWccS9pjt+axOIN+swI0hJjeTV+mUx6KKJuBj/RP0LdYOCMOt0XB00EGNkJOdfbc/Nn+bkmhMtOh2MGQO//mrphqlWDWxoFuFcpUmTJvz555+cPn2ahg0bEqHmQ84ULp5evBIyisToaFZOGIvJmKZ1pGxhcwX9yy+/5PKKy9xMvrff0ec1H1z+TqHAZRPJOJKaEsGLAV6kOqZyK9kPIjLeSs+XLx/ly5fn0PXrFFq0EJd69Yj46isu9uhJWi4Yp35bp06WGRvT0qB2bVi2TOtEeVOjRo1YtWoVZ8+epWHDhly/fl3rSLlCwRKlaNp3EOEnjrF5ziyt42QLmyvov/32Gy9/+DJRSfeODrg92qXG1jSScCYl1XKxkUvxwkSbChO+88kuDa5Xrx6pqakcOH0a/2lT8RvzGUlHj3Lu1fYPzAdjywIDYd8+eP55aNfO0nJXJ0uzX8OGDVm9ejXnz58nKCiIq1evah0pVyhXN4hqrV/l8PpQDq9frXWcLGczBT1uyyWSz0Tj6+uLfyF/bibfJPlMNHFbLGtvOAY4YqriTI2taaQa/EhNsXx1rVGrNACRhw8/0ev5+vpSunRpdu3aRVpaGp7t21N88WKEvT0Xu3Un+UTGW/w5nZ8fbNkCXbrA6NGWlnsevjZDMw0aNGD16tVcunSJhg0b5oqrlnOCup27UvyFQDb99B3hx49pHSdL2UxBt/N3I+qXE2yet5qru69S+EY+on45gZ2/25199K94UfKkCRlZipRUywmmYqUtC0mLeB2Jt56s1XP/YtIOJYoTMG8uwsmJi917kPT335n07rTn6Ahz58JXX8GiRZZpA8LDtU6V99SvX581a9Zw+fJlgoKCuHz5staRbJ5Op+flQcPwKOjHioljiY3MvecpbKagO5b0xPv1chQ8CM9vKsDwSz1xeK0ojiU97+zj2tZyub/D1vJ3WuguHg7YO8GNtBIcC1v0RK9ZtGhRAgIC7llM2r5oUQLm/oxwceZij54kHc09f/GFgPffhxUrLFeUBgZaJvpSslfdunVZu3Yt165dIygoiHD1l/WZOTi78MqwUZiNRv4YP4a05GStI2UJmynoYCnq7lX8eK1CCw64HCfG795ViDxLOXO+pA7nXQVJSf3vr7BvCW+uGksjT4Q+8Wumt5i0fZEiBMydh97NjYs9e5KUy4aItGwJO3eCi4vlgqS5c7VOlPfUrl2btWvXcv36dYKCgrh06dLjn6Q8knchf14eNIzIC+dYM2NSrpxD3aYKevKZaHTnkwCoGV+JuFP3fnVy1+s59bwBx+POpCbfRErLMMMCAe7EGgtRLvkYEbein+g1S5YsiZ+f3wOLSdv7FyZg3lz0np5c7PkWSblspebnn7fM2FinDnTrBsOGqRkbs1utWrVYv349kZGRBAUFcVGtWvLMir8QSP3Xu3Nq13b2/PFk39htgc0U9OQz0UT9coK/i95g1aVt/ON0Do9VqSSfib6zj5tBz5nn9BjidciL+UlLs8xo51PUDdCRaPTjn51Pdqb7UYtJ2xUqRMDcn9Hn8+biW71IPHDwWd9mjpIvH6xdC/37w/jx0Lq1mrExu9WoUYP169dz8+ZNGjRowPnz57WOZPMCW7WjXN0gtv8+jzMZuIrclmSooAshmgshTgohTgshhqfzuBBCTLY+fkQIUTWzg6aFx+H9ejk+n/MNP+9fiqfJnb/rR5AW/t985446wYXydpZMJ0rcOTFqKehwxVgG88knH7r03HPPkT9//nQXk7bz8yNg7lwMPj5c6tWLxH37nvYt5kh2djB1KsyYAevWQa1aYKPLr9qs6tWrs2HDBqKjowkKCuKcdWEW5ekIIWjSZyAFi5cidMp4bobnnm8+jy3oQgg9MA1oAZQHOgshyt+3WwugtPWnNzAjk3Pi1qAIjiU9Wb58Ob9+/iP+KQU553YVtwZF7s5KfHE70lyBf4qTmmK5QMPVywFHFzvOU53S0TtIMz5Z38Hdi0n/m878s3YFC1J07s8YChbkYu8+JOzZ80zvNSfq29cyS2NEhGXGxo0btU6UtwQGBrJx40ZiY2MJCgri7Nm8NUdJZrOzd6BNyEgM9g788fVnJMXnjjVfxeNODAghagEfSymbWbdHAEgpx961z0wgTEr5q3X7JBAkpXzoOEE3Nzc5ZcoUunfvTlpaGk2aNKFXr1506dKFxMREgoOD6devHx07diQmJoY2bdowaNAg2rVrx8WwE3Ts+jpda7QlsEJFbsTdYtTiCXSr25afXm9N+/cvsOjvMbxc4A0qOr3I9bQrzLw5nvae3SjvWJkraZf44eZEOnm+RRnH57mUeo6foqbwhldvSjo8x/nU08yNmk5X73coZl+KMyn/sODW9/TwHkgR++KcSv6b36J/pFe+IRSyK8Lx5MMsjv6ZPvlCKGhXiKNJ+1kWs4D++YeTz1CAQ0l7WBHzG4N8RuGp92Z/4g5WxS5msM9HuOs92JO4jTWxywgp8CnOOld2JmxmfdxKPijwBQ46R7bFr2dz/Gr+V/BrDMLAlvg1bIlfx4e+EwDYGLeKXYlhjCw4DoB1ccs5kLiT4QW/BGB17BKOJR9kWIExAPwZs5B/U48zxOdjAJbH/MqF1NMM8hkNwNLoeVwxXmJA/v8BsCj6J24aI+mb/32kUbDg6k/Em2PpYWf5svabaTKppNBVPwyABaaJALyhHwLAXNM47HGgk34QAD8Zv8RVuPOa/h0AfjCOwVsUoJ2+NwAzjR/jK4rSRt8TgOnG0RQVpWmp7wrAFOMISokKtNC/AcAk4zDKi0Ca6jsC8I1xCFVEXRrpXwXga+NAqusaE6RrA8BY4zvU1b1MPd3LGKWRcaZBNNC1prauOSkymQmmobyka0cNXWMSZTzfmt6nia4Dgbog4mQ0U03/o7muMy/o6hEtbzLDNJqXdW9SSVeLm/I635s+obWuO8/rqhMhL/Oj6XPa6nrxnK4qV+UF5pi+or2uL6V1lQiXZ5hn+oaOugGU0JXngjzFL6ZJvK4fTIAow1nzcX43T+VN/Xv4i5JsM//JT6YvsccRbwqQQhJxRONJfgzY3dn2wgc9BpJJJJ6YdLYLoEdPEgkkEIs3BdGhI4l4Eoi7s51IPInEkQ9fBIJE4kgknvz4AZBALEkk3LOdTCL58AUgnhhSSCYfBe9sp5KMt3U7jmjSSMWbAtbtWxhJw8u6HcstTBjxwse6HYUZM55Y5m+K4SYSec82gAeWEW/R3EAg7tnWocMdbwBuEYkeA+54WbcjMGCHm3U7igjssMcNT+v2dexxxBUPAG5yHYd7tq/hiDMuuANwg6s44XJn+7YyhZuzOXw6T0MIsV9KGZjeY4YMPL8wcPcp9nCgRgb2KQzcU9CFEL2xtOBxcHDIwEs/6Pfffyf6ZjSpTiZuesQR555CAqmYDGYSndN4OfIG/zaWJF0yk1ggjVjPFOKTUzHFW7fdU4hPtGwnFEi1bCekYkowk1AgjVi3FBLib2+nEuuaQkJcGqZEM3G+KUS7JBMXm4IxyUysbwrOzsnEx6Ratgul4OCYTHx0KsZkMzGFUtA7JJNw6/Z2Mtgnk3AzDWOKmdjCyZjtHEi4kYoxxUxM4RRSDQYSItMwppqJ8U/GXg+JEWkY0yzbep2BxOuW7egilqFXSdfSSDOZ/tu+mkaavOvxy0bSYu56XKSRGvffdjJppCb8t3+yNJKWdNfjZiOpKf9tp8lkkqIlV50sX/ASEgRpUnDV9b9tgKsulu3EeEGaEHdtA3dtJ8VDvE5w1dm6HSeI0/+3nRwHcQZx5/WS4yD2ru2UWEGs3b3bMfaCq47pb6fGCGIcLNsmqSM1VhDtKLjqoCP1vu1ks47UOMEt63bC7W0nHVftdcRZt6Os29EmHanx/21HWbdvOuu4aqcj8q5tVzsdkUYdqQmCGy46nAw6bty1bX/XdqSLDr1Bh11aMfwSSuOmy49B2JNgvkWSKRlPQxHshRPx5iiSTMl4GIpgLxyJM98k2ZSCp6EodsKBWPMNkk0peBmKYhD2xJgjSDGl4mUIQC8MxJiuk2JOw9tQDJ3QozNdI9W6LYQOYbpKqtmIt11xyy+k6TJGs/nOtjSFYzLLO9tm0yWk+dZd2xfBHHNn22Q6T7KMx9tg2TaaIEUm3tlOM0rSSL6znWo0YyLtznaK0YjEfNe2Zc6W29vJxlQEujvbScYU9Njd2U40JmGH453tBGMiDsIZb/3t7QQchSve+mIAxKfF46zzwFtfFIC4tDicdV546y09BbFpMTjr8uGt9wcgJi0aF11+vPWFuZuzuytZISMt9NeAZlLKXtbtN4HqUsqBd+2zChgrpdxu3d4IvC+l3P+w4wYGBsp9T9HfHBQUBEBYLluEQlEUJSOetYUeDhS5a9sfuP+a5IzskylCQ598LLmiKEpekJFRLnuB0kKI4kIIe6ATcP9MVyuArtbRLjWBmEf1nz8LZ2dnnJ2ds+LQiqIoNu2xLXQppVEIMQBYC+iB2VLKv4UQfa2PfweEAsHAaSAR6JFVgefPnw9Aly5dsuolFEVRbNJj+9CziupDVxRFeXLP2oeeo6xXS9YriqKky+YKup2dndYRFEVRciSbmcvltjlz5jBnzhytYyiKouQ4qqAriqLkEpqdFBVCRAIXnvLp+YEbmRgnu9lyfpVdGyq7dnJa/gAppU96D2hW0J+FEGLfw87y2gJbzq+ya0Nl144t5be5LhdFURQlfaqgK4qi5BK2WtC/1zrAM7Ll/Cq7NlR27dhMfpvsQ1cURVEeZKstdEVRFOU+qqAriqLkEjZX0B+3YHVOIIQ4L4Q4KoQ4JITYZ73PWwixXgjxr/Vfr7v2H2F9PyeFEM2yOetsIUSEEOLYXfc9cVYhxIvW93zaumC40Cj7x0KIy9bP/pAQIjiHZi8ihNgshDghhPhbCPGu9f4c/9k/IrutfPaOQog9QojD1vyfWO/P8Z/9Y0kpbeYHy/S9Z4ASgD1wGCivda50cp4H8t9339fAcOvt4cBX1tvlre/DAShufX/6bMxaH6gKHHuWrMAeoBYggNVAC42yfwyEpLNvTsvuB1S13nYDTlkz5vjP/hHZbeWzF4Cr9bYdsBuoaQuf/eN+bK2FXh04LaU8K6VMBX4D2micKaPaAD9bb/8MvHLX/b9JKVOklOewzClfPbtCSSm3AlH33f1EWYUQfoC7lHKntPxfPveu52R39ofJadmvSikPWG/HASewrMOb4z/7R2R/mByT3ZpZSinjrZt21h+JDXz2j2NrBf1hi1HnNBJYJ4TYLywLYwMUlNZVnKz/FrDenxPf05NmLWy9ff/9WhkghDhi7ZK5/bU5x2YXQhQDXsDSUrSpz/6+7GAjn70QQi+EOAREAOullDb32afH1gp6ev1TOXHcZR0pZVWgBdBfCFH/EfvaynuCh2fNSe9hBlASqAJcBb6x3p8jswshXIElwGApZeyjdk3nPk3zp5PdZj57KaVJSlkFy/rH1YUQFR6xe47L/zC2VtCzbTHqZyGlvGL9NwJYhqUL5br1KxrWfyOsu+fE9/SkWcOtt++/P9tJKa9bf1nNwCz+677KcdmFEHZYCuICKeVS69028dmnl92WPvvbpJTRQBjQHBv57B/F1gp6Rhas1pQQwkUI4Xb7NtAUOIYlZzfrbt2A5dbbK4BOQggHIURxoDSWEy1aeqKs1q+ncUKImtaz/F3vek62uv0LadUWy2cPOSy79bV+BE5IKSfc9VCO/+wflt2GPnsfIYSn9bYT0Bj4Bxv47B9LyzOyT/ODZTHqU1jONI/UOk86+UpgOSN+GPj7dkYgH7AR+Nf6r/ddzxlpfT8nyeaz5MCvWL4ep2Fpcbz1NFmBQCy/wGeAqVivQtYg+zzgKHAEyy+iXw7NXhfL1/MjwCHrT7AtfPaPyG4rn30l4KA15zHgQ+v9Of6zf9yPuvRfURQll7C1LhdFURTlIVRBVxRFySVUQVcURcklVEFXFEXJJVRBVxRFySVUQVcURcklVEFXFEXJJf4P+il3bvCYkG8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Gamma, S, *_ = pattern_kinks_LASSO(Xapp, yapp)\n", "\n", "# Calcul du chemins des solutions du LASSO et des composantes des moindres carrés en valeur absolue\n", "Sol = solution_LASSO(Xapp, yapp, Gamma, S)\n", "Sol_abs = [np.abs(sol) for sol in Sol]\n", "ols = np.linalg.solve(Xapp.T@Xapp, Xapp.T@yapp)\n", "ols_abs = np.abs(ols)\n", "\n", "plt.plot([1.1*Gamma[0]] + Gamma, [Sol[0]] + Sol_abs)\n", "plt.axvline(0, color='k', linestyle=':')\n", "plt.axhline(0, color='k', linestyle=':', xmax=0.95)\n", "plt.plot(0, [np.abs(ols)], 'x')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Minimisation de la somme des carrés résiduels " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Somme des carrés résiduels :** " ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [], "source": [ "def SCRV(Xval,yval,Sol):\n", " return(np.linalg.norm(yval-Xval @ Sol, ord=2)**2)" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7895.051474323932\n", "2641.6279680942725\n" ] } ], "source": [ "zero_estimator = np.zeros(X.shape[1])\n", "print(SCRV(Xval,yval,zero_estimator))\n", "ols_app = np.linalg.solve(Xapp.T@Xapp, Xapp.T@yapp)\n", "print(SCRV(Xval,yval,ols_app))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***Minimisation de la somme des carrés résiduels :*** le code suivant permet de minimiser la somme des carrés résiduels sur l'échantillon de validation de l'estimateur SLOPE " ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "def critical_points(Xval,yval,m,a_s,b_s):\n", " Xvaltilde = Xval @ pattern_matrix(m)\n", " N=np.transpose(a_s) @ np.transpose(Xvaltilde) @ (yval -Xvaltilde @ b_s)\n", " D=np.linalg.norm(Xvaltilde @ a_s, ord=2)**2\n", " c=N/D\n", " return(c)\n", "\n", "def min_SCR(Xapp,yapp,Xval,yval,Lambda,Gamma,M):\n", " SCR=list()\n", " Critical = list()\n", " Solution = list()\n", " kmax=len(Gamma)-1 \n", " for k in range(kmax):\n", " m = M[k] \n", " a_s,b_s, *_ = affine_components(Xapp, yapp, Lambda, m)\n", " c=critical_points(Xval, yval, m, a_s, b_s)\n", " gamma_0=Gamma[k]\n", " gamma_1=Gamma[k+1]\n", " Sol0=pattern_matrix(m) @ (a_s * gamma_0 + b_s )\n", " scr=SCRV(Xval,yval,Sol0)\n", " SCR.append(scr)\n", " Critical.append(gamma_0)\n", " Solution.append(Sol0)\n", " Sol1=pattern_matrix(m) @ (a_s * gamma_1 + b_s ) \n", " if cgamma_1:\n", " a= (c-gamma_1)/(gamma_0-gamma_1)\n", " Solc=a*Sol0+(1-a)*Sol1\n", " scr=SCRV(X,y,Solc)\n", " SCR.append(scr)\n", " Critical.append(c)\n", " Solution.append(Solc)\n", " m = M[-1]\n", " a_s,b_s, *_ = affine_components(Xapp, yapp, Lambda, m)\n", " gamma_0=Gamma[-1]\n", " Sol0=pattern_matrix(m) @ (a_s * gamma_0 + b_s )\n", " scr=SCRV(Xval,yval,Sol0)\n", " SCR.append(scr)\n", " Critical.append(gamma_0)\n", " Solution.append(Sol0)\n", " return Critical, SCR, Solution" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2159.045472692249\n", "357.38588754081786\n", "[-0.75067428 0.60605446 -0.2459192 0.60605446 -0.75067428 3.21594712\n", " 0. -1.70945109 0.75067428 -0.4161687 -1.54307232 0.38658049\n", " -3.21594712]\n", "[-5 4 -1 4 -5 8 0 -7 5 -3 -6 2 -8]\n" ] } ], "source": [ "Lambda=np.sqrt(range(1,X.shape[1]+1))-np.sqrt(range(0,X.shape[1]))\n", "\n", "Gamma, M, T, Split = pattern_kinks(Xapp, yapp, Lambda, ratio=0)\n", "C,S,Solution = min_SCR(Xapp, yapp, Xval,yval,Lambda,Gamma,M) \n", "\n", "i0=np.argmin(S)\n", "print(S[i0])\n", "print(C[i0]) \n", "print(Solution[i0]) \n", "print(pattern(Solution[i0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***Minimisation de la somme des carrés résiduels :*** le code suivant permet de minimiser la somme des carrés résiduels sur l'échantillon de validation de l'estimateur LASSO" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "def critical_points_LASSO(Xval,yval,s,a_s,b_s):\n", " I=s!=0\n", " N=np.transpose(a_s) @ np.transpose(Xval[:,I]) @ (yval -Xval[:,I] @ b_s)\n", " D=np.linalg.norm(Xval[:,I] @ a_s, ord=2)**2\n", " c=N/D\n", " return(c)\n", "\n", "def min_SCR_LASSO(Xapp,yapp,Xval,yval,Gamma,S):\n", " SCR=list()\n", " Critical = list()\n", " Solution = list()\n", " kmax=len(Gamma)-1 \n", " for k in range(kmax):\n", " s = S[k] \n", " a_s,b_s, *_ = affine_components_LASSO(Xapp, yapp, s)\n", " c = critical_points_LASSO(Xval, yval, s, a_s, b_s)\n", " gamma_0=Gamma[k]\n", " gamma_1=Gamma[k+1]\n", " Sol0 = np.zeros(X.shape[1],dtype=float)\n", " Sol1 = np.zeros(X.shape[1],dtype=float)\n", " I=s!=0\n", " Sol0[I] = a_s * gamma_0 + b_s \n", " scr = SCRV(Xval,yval,Sol0)\n", " SCR.append(scr)\n", " Critical.append(gamma_0)\n", " Solution.append(Sol0)\n", " Sol1[I]=a_s * gamma_1 + b_s \n", " if cgamma_1:\n", " a = (c-gamma_1)/(gamma_0-gamma_1)\n", " Solc = a*Sol0 + (1-a)*Sol1\n", " scr=SCRV(Xval,yval,Solc)\n", " SCR.append(scr)\n", " Critical.append(c)\n", " Solution.append(Solc)\n", " s = S[-1]\n", " a_s,b_s, *_ = affine_components_LASSO(Xapp, yapp, s)\n", " gamma_0=Gamma[-1]\n", " I=s!=0\n", " Sol0[I] = a_s * gamma_0 + b_s \n", " scr = SCRV(Xval,yval,Sol0)\n", " SCR.append(scr)\n", " Critical.append(gamma_0)\n", " Solution.append(Sol0)\n", " return Critical, SCR, Solution" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2498.5172935021847\n", "42.61417506197384\n", "[-1.06406262 0.72483679 0. 0.51236735 -1.02410679 3.45529462\n", " 0. -2.21375338 1.93336031 -1.15547175 -1.5784098 0.31916939\n", " -3.74940327]\n", "[-1. 1. 0. 1. -1. 1. 0. -1. 1. -1. -1. 1. -1.]\n" ] } ], "source": [ "Gamma, SGN, T = pattern_kinks_LASSO(Xapp, yapp)\n", "C,S,Solution = min_SCR_LASSO(Xapp,yapp,Xval,yval,Gamma,SGN) \n", "\n", "i0=np.argmin(S)\n", "print(S[i0])\n", "print(C[i0]) \n", "print(Solution[i0]) \n", "print(np.sign(Solution[i0]))" ] } ], "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }