basic_demo_revised.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstartes the steps in Stability-based model selection for learning PDEs from data (PDE-STRIDE)\n",
    "$$\n",
    "u_t + uu_x = 0.1u_{xx}\n",
    "$$\n",
    "In this demo notebook, we show stability selection in combination with three different sparsity promoters LASSO, STRidge and IHT with debiasing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (12, 8)\n",
    "import numpy as np\n",
    "import math\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "#%matplotlib notebook\n",
    "from Iterative_solvers_l0_l1 import *\n",
    "import scipy.io as sio\n",
    "import itertools\n",
    "from operator import itemgetter, attrgetter\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import lars_path, enet_path, lasso_path\n",
    "from itertools import cycle\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import datasets, linear_model\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "import seaborn as sns\n",
    "c = sns.color_palette(\"Blues\",20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " shape of u  (256, 101)\n",
      " dt is  0.1  dx is  0.0625\n"
     ]
    }
   ],
   "source": [
    "X_APR    = loadtxt('Burgers/X_burgers.dat', dtype='double')\n",
    "T_APR    = loadtxt('Burgers/T_burgers.dat', dtype='double') \n",
    "data = sio.loadmat('Burgers/burgers.mat')\n",
    "u = real(data['usol'])\n",
    "dx = X_APR[1] - X_APR[0]\n",
    "dt = T_APR[1] - T_APR[0]\n",
    "print(\" shape of u \", u.shape)\n",
    "print(\" dt is \", dt, \" dx is \", dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create the design matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " The value of D is changed to 4  4\n",
      " length of description  20\n",
      "i : 1 u\n",
      "i : 2 u^2\n",
      "i : 3 u^3\n",
      "i : 4 u_{x}\n",
      "i : 5 uu_{x}\n",
      "i : 6 u^2u_{x}\n",
      "i : 7 u^3u_{x}\n",
      "i : 8 u_{xx}\n",
      "i : 9 uu_{xx}\n",
      "i : 10 u^2u_{xx}\n",
      "i : 11 u^3u_{xx}\n",
      "i : 12 u_{xxx}\n",
      "i : 13 uu_{xxx}\n",
      "i : 14 u^2u_{xxx}\n",
      "i : 15 u^3u_{xxx}\n",
      "i : 16 u_{xxxx}\n",
      "i : 17 uu_{xxxx}\n",
      "i : 18 u^2u_{xxxx}\n",
      "i : 19 u^3u_{xxxx}\n",
      " shape of the truncated description  19\n",
      "shape 19116 20\n",
      "shape of no_intercept (19116, 19)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/lustre/projects/jupyterhub/maddukon/Data_driven/Stability_selection_results /PDE-STRIDE/Iterative_solvers_l0_l1.py:35: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  X_no_intercept[:,i-1] = X[:,i]\n"
     ]
    }
   ],
   "source": [
    "numpy.random.seed(0)\n",
    "noise_level = 0.01\n",
    "un = u + noise_level*std(u)*np.random.randn(u.shape[0],u.shape[1])\n",
    "\n",
    "#y, X, description = build_linear_system(u, dt, dx, D=4, P=3, time_diff = 'FD', space_diff = 'FD')\n",
    "# In this demo, we don't use SVD based denoising but simply use gaussian blur + polynomial differentiation\n",
    "y, X, description = build_linear_system(un, dt, dx, D=4, P=3, time_diff = 'poly', space_diff = 'poly',deg_x = 5, deg_t = 5, \n",
    "                                      width_x = 10, width_t = 10)\n",
    "print(\" length of description \", len(description))\n",
    "description_trunc, X_no_intercept = remove_intercept(X, description) # This step removes the constant column in the design\n",
    "u_t = y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute regularization path for a random sub-sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " lambda max is  [0.06333404]  lambda min is  [0.0071062]\n",
      " lambda max is  [1.]  lambda min is  [0.11220185]\n",
      " norm is  1805.231355509598\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/sw/apps/python3/3.6.5/lib/python3.6/site-packages/ipykernel_launcher.py:13: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  del sys.path[0]\n",
      "/lustre/projects/jupyterhub/maddukon/Data_driven/Stability_selection_results /PDE-STRIDE/Iterative_solvers_l0_l1.py:419: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if lam != 0: coeff[0,biginds] = np.linalg.lstsq(X[:, biginds].T.dot(X[:, biginds]) + lam*np.eye(len(biginds)),X[:, biginds].T.dot(y))[0].reshape(len(biginds),)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " done with alpha  0\n",
      " converged at iteration  4  frac_iter  0.0  flag  1\n",
      " done with alpha  1\n",
      " converged at iteration  4  frac_iter  0.0  flag  1\n",
      " done with alpha  2\n",
      " converged at iteration  4  frac_iter  0.0  flag  1\n",
      " done with alpha  3\n",
      " converged at iteration  5  frac_iter  0.0  flag  1\n",
      " done with alpha  4\n",
      " converged at iteration  5  frac_iter  0.0  flag  1\n",
      " done with alpha  5\n",
      " converged at iteration  12  frac_iter  0.07692307692307693  flag  1\n",
      " done with alpha  6\n",
      " converged at iteration  10  frac_iter  0.09090909090909091  flag  1\n",
      " done with alpha  7\n",
      " converged at iteration  11  frac_iter  0.08333333333333333  flag  1\n",
      " done with alpha  8\n",
      " converged at iteration  11  frac_iter  0.08333333333333333  flag  1\n",
      " done with alpha  9\n",
      " converged at iteration  11  frac_iter  0.08333333333333333  flag  1\n",
      " done with alpha  10\n",
      " converged at iteration  14  frac_iter  0.06666666666666667  flag  1\n",
      " done with alpha  11\n",
      " converged at iteration  17  frac_iter  0.05555555555555555  flag  1\n",
      " done with alpha  12\n",
      " converged at iteration  16  frac_iter  0.058823529411764705  flag  1\n",
      " done with alpha  13\n",
      " converged at iteration  23  frac_iter  0.041666666666666664  flag  1\n",
      " done with alpha  14\n",
      " converged at iteration  24  frac_iter  0.04  flag  1\n",
      " done with alpha  15\n",
      " converged at iteration  27  frac_iter  0.03571428571428571  flag  1\n",
      " done with alpha  16\n",
      " converged at iteration  29  frac_iter  0.03333333333333333  flag  1\n",
      " done with alpha  17\n",
      " converged at iteration  31  frac_iter  0.03125  flag  1\n",
      " done with alpha  18\n",
      " converged at iteration  32  frac_iter  0.030303030303030304  flag  1\n",
      " done with alpha  19\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(50)\n",
    "reduced_size = 200\n",
    "ordered = np.arange(0, len(u_t))\n",
    "np.random.shuffle(ordered)\n",
    "\n",
    "X_reduced = np.zeros((reduced_size, len(description_trunc))) # random sub-sample from the design matrix\n",
    "y_reduced = np.zeros((reduced_size,1))                       # random sub-sample from the measurement vector\n",
    "X_stand   = np.zeros((reduced_size, len(description_trunc))) # standardized sub-design matrix\n",
    "y_stand   = np.zeros((reduced_size,1))                       # standardized sub-measurement vector\n",
    "\n",
    "for i in range(0, reduced_size):\n",
    "    X_reduced[i,:] = X_no_intercept[int(ordered[i]),:]\n",
    "    y_reduced[i,:] = u_t[int(ordered[i])]\n",
    "\n",
    "X_stand = preprocessing.scale(X_reduced)\n",
    "y_stand = preprocessing.scale(y_reduced, with_std=False)\n",
    "y_stand = y_stand.reshape(reduced_size,1)\n",
    "\n",
    "mul_vector = np.matmul(X_stand.transpose(), y_stand)     \n",
    "lambda_max = (1.0/reduced_size)*np.linalg.norm(mul_vector, np.inf)  # lambda_{max} in paper \n",
    "M             = 20 # M in the paper\n",
    "factor        = np.exp(np.log(10.)/M)\n",
    "lambda_regl   = np.zeros((M,1))\n",
    "coeff_alpha   = np.zeros((len(description_trunc),M)) # stores the output\n",
    "\n",
    "alphas_regl = np.zeros((M,1))\n",
    "for alpha in range(0, M):\n",
    "    alphas_regl[alpha] = lambda_max/factor**alpha\n",
    "    \n",
    "normalized_alphas_sample = np.zeros((M,))\n",
    "normalized_alphas_sample = alphas_regl/alphas_regl[0]\n",
    "neg_log_stable_lambda = -np.log10(normalized_alphas_sample)\n",
    "\n",
    "print(\" lambda max is \", alphas_regl[0], \" lambda min is \", alphas_regl[M-1])\n",
    "print(\" lambda max is \", normalized_alphas_sample[0], \" lambda min is \", normalized_alphas_sample[M-1])\n",
    "\n",
    "p = X_stand.shape[1]\n",
    "coefs_lasso   = np.zeros((p, M))\n",
    "coefs_ihtd    = np.zeros((p, M))\n",
    "coefs_STR     = np.zeros((p, M))\n",
    "\n",
    "print(\" norm is \", np.linalg.norm(X_stand.T.dot(X_stand),2))\n",
    "\n",
    "for alpha in range(0, M):\n",
    "    clf_l = linear_model.Lasso(alphas_regl[alpha], fit_intercept=False, max_iter=10000, tol=1e-8) #L1\n",
    "    clf_l.fit(X_stand,y_stand)\n",
    "    coefs_lasso[:,alpha] = clf_l.coef_\n",
    "        \n",
    "    clf_ihtd = Iterative_hard_thresholding_debias(X_stand, y_stand, alphas_regl[alpha], max_iter=10000, sub_iter=25, tol=1e-8, htp_flag = 1, print_flag=True)\n",
    "    coefs_ihtd[:, alpha] = clf_ihtd[:]\n",
    "    \n",
    "    clf_STR = STRidge(X_stand, y_stand, alphas_regl[alpha], lam = 10**-5, max_iter=10000)\n",
    "    coefs_STR[:,alpha] = clf_STR[:]\n",
    "\n",
    "    print(\" done with alpha \", alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,' $ -log(\\\\lambda^{*}) $')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8IAAAEQCAYAAAB7gs2gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmczfX+wPHXZ/Z9xizGDGaxZ5dJsl0UIaFSUaFF4ib100K3kq0u6gqpKFc3uaFSoc0lpGhhJCJjm8HYx4yZYfaZz++Pc2Y645wz6zlzZnk/H4/zmHO+7+/yPoyv8z6fTWmtEUIIIYQQQggh6gonRycghBBCCCGEEEJUJSmEhRBCCCGEEELUKVIICyGEEEIIIYSoU6QQFkIIIYQQQghRp0ghLIQQQgghhBCiTpFCWAghhBBCCCFEnSKFsBBCCCGEEEKIOkUKYSGEEEIIIYQQdYoUwkIIIYQQQggh6hQphIUQQgghhBBC1Ckujk6gKgUHB+uoqChHpyGEqEZiY2OTtNYhjs7DluReJ4S4li3vdUqpAcBCwBlYprWec028F7AAaA+M0Fp/ahIbA7xofDlba/2BcXtn4D+AJ/A18KTWWpeUh9zrhBCWlPV+V6cK4aioKHbv3u3oNIQQ1YhS6oSjc7A1udcJIa5lq3udUsoZeAvoByQCu5RS67XWB012Owk8CDxzzbGBwMtADKCBWOOxKcA7wKPALxgK4QHANyXlIvc6IYQlZb3fSddoIYQQQghRVl2Ao1rr41rrHGA1MNR0B611gtZ6H1BwzbG3Apu01snG4ncTMEApFQb4aa1/NrYCrwCG2f2dCCHqNCmEhRBCCCFEWTUETpm8TjRuq8yxDY3PK3JOIYSoECmEhRBCCCFEjaCUGqeU2q2U2n3x4kVHpyOEqMGkEBZCCCGEEGV1Gmhs8rqRcVtljj1tfF7qObXW72qtY7TWMSEhtWqeQyFEFatTk2UJIWqnR1btIzUrz2y7v4cL/x7Z3gEZidqkIr9f5T2mtlyjuuZVl69hB7uA5kqpaAzF6gjgvjIeuxF4VSlVz/i6P/C81jpZKZWmlOqKYbKs0cCbNs5bCCGKkRZhIUSNZ+mDYUnbhSiPivx+lfeY2nKN6ppXXb6GrWmt84CJGIraP4GPtdYHlFIzlVJDAJRSNyilEoG7gaVKqQPGY5OBWRiK6V3ATOM2gL8Dy4CjwDFKmTFaCCEqS1qEhRBCiAradfKy3Y+pLdeoyDFyjepJa/01hiWOTLdNM3m+i+JdnU33Ww4st7B9N9DWtpn+5cqVK2zZssVse/v27YmKiuLy5cts377dLH799dfTqFEjkpKS2Llzp1m8S5cuNGjQgHPnzvHrr7+axbt160ZwcDCJiYns2bPHLN6rVy8CAgJISEhg3759ZvG+ffvi4+PD0aNHOXjwoFm8f//+eHh4cOjQIQ4fPmwWHzRoEC4uLhw4cIBjx46ZxYcMGQLA77//zokTxVeccXZ25rbbbgMgNjaW06eL91Z3d3fn1ltvBWDXrl20bNkSPz8/s2sIUV2pUtYqr1ViYmK0rDcnRO0z/H3zDxeFPn3o+hKPVUrFaq1jbJ2TI8m9zrZK+v0SorqQe13J4uLiaNWqldn2pUuXMm7cOHbv3s0NN9xgFv/oo48YOXIkW7dupW/fvmbxDRs2MHjwYDZs2FBUVJratm0bf/vb3/joo4+4//77zeK7d++mc+fOLF26lPHjx1vMu0WLFrz++us8++yzZvEzZ84QFhbGyy+/zMyZM83i6enp+Pj48PTTTzN//nyzeGEdMG7cON57771iMR8fH9LT0wEYOXIkq1evLhYPCwvjzJkzAHTo0IGBAwcyZ84cs2sIUdXKer+TFmEhRI2WkJzh6BREHTbvdvMP1gDPbThUrmPKu391vUZ1zasuXEOULDIyktjYWLPtERERAFx33XUW41FRUQDExMRYjDdt2hSAnj17Wow3b94cgAEDBliMFxbnd955p8VCvDC/UaNGWSzEg4ODAXjssccYOnSoWdzT0xOA//u//7NYiBd64YUXzApxJ6e/RlC++uqrZoW4q6tr0fNz585x5MgRq+cXojqSQlgIUSOlZ+exZs9Z/hcny2cIx2kS7GX3Y2rLNSpyjFxD2IqHhwfXX2+91dzb27vEuK+vb4nxgICAEuOBgYEEBgZajYeEhFDSLNihoaGEhoZajYeHhxMeHm413qhRIxo1sthbHTB8URAZGWk1Hh0dTXR0tNX49ddfb9a1WojqTibLEkLUKPkFmv8dusiktQf4X9xFbm0Vgp+H5e/0/K1sF6I8rP0elfT7Vd5jass1qmtedfkaQlSFiIgITp486eg0hCgXGSMshKgx/jx/hX//fIqE5ExaN/Dh4RsbERVYuVaTuj5uTpTNt39eZNnPp3jn7raE+Lg5Oh0hyk3udcKeXn31VV544QWuXr2Kl5f0ZhCOJWOEhRA1lrW1MgGCvFyZ3Duam6ICUEpVcWairopPzsDH3Zlgb9fSdxZCiDqmcPbqgoICR6ciRJlJISyEqHZKWhNz4Z2t8XB1rsJshID4S5lEB3rJly9CCGFBTEwMMTG1qsOBqANkjLAQokaRIlhUtbwCzanLmUQHeTo6FSGEqJYKCgo4duwYZ8+edXQqQpSZFMJCiGolMzff0SkIUczpy1nk5muiKzkeXQghaquCggJatmzJ22+/7ehUhCgz6RothKg29p9J550dsvyCqF7ijWtVS4uwEEJY5uLiQsOGDWUJJVGjSCEshHC4jJx8Ptx9mk1xSYT5uTs6HSGKib+UiZuzIszPw9GpCCFEtRUZGSlLKIkaRbpGCyEcau/pNCZ/cZDvDidxe5v6vD70OlkrU1Qr8ckZRAZ64uwkE2UJIYQ1ERER0iIsahT5VCmEcIgr2Xl8sOs0W49copG/B7MHNaFFfW8A/j2yvYOzE8JAa01CciY9ous5OhUhhKjWIiMjWbNmDfn5+Tg7y8SWovqTQlgIYVfW1gRWgFJwR/tQ7u4QhpuLdFAR1c+FKzlk5OQTHSQTZQkhREnuueceOnbsSEFBgRTCokaQQlgIYVfW1gTWwJzBrWgaLAWGqL7iLxknygqUibKEEKIkHTp0oEOHDo5OQ4gykyYYIYTDSBEsqrv45EycFETUk0JYCCFKkpOTw7Zt24iPj3d0KkKUiRTCQgghhBXxlzJo5O8hXfeFEKIUubm59OnThzVr1jg6FSHKRP5nF0LYzbm0bEenIESlJCRnEiXjg4UQolTe3t4EBQXJzNGixnBoIayUGqCUilNKHVVKTbUQd1dKrTHGf1FKRRm3BymltiqlriilFld13kKI0v14PJln1//p6DSEqLDUzFySM3JlfLAQQpRRRESErCUsagyHFcJKKWfgLWAg0BoYqZRqfc1ujwApWutmwBvAXOP2LOAl4JkqSlcIUUbZeQUs2XGCBd8n0DjAEz93yzNHyprAorqLT84EkBmjhbhGJRoy7ldK7TV5FCilOhpj24znLIzVr9p3JWwhMjJSWoRFjeHIT6JdgKNa6+MASqnVwFDgoMk+Q4HpxuefAouVUkprfRX4USnVrArzFUKU4tTlTOZvjefU5SzuaBfKvdeH4+KkHJ2WEBVSOGN0lLQIC1HEpCGjH5AI7FJKrddam35+K2rIUEqNwNCQca/W+r/Af43naQd8obXea3Lc/Vrr3VXyRoRdREZGsnnzZrTWKCX//4vqzZFdoxsCp0xeJxq3WdxHa50HpAJB5bmIUmqcUmq3Umr3xYsXK5GuEMIarTXfHU5iyvpDpGbl8WL/Ztwf01CKYFGjxSdnUt/HDR936b0ghImihgytdQ5Q2JBhaijwgfH5p8DNyrwqGmk8VtQiEyZM4H//+5+j0xCiTGr9/+5a63eBdwFiYmK0g9MRosZ7ZNU+q2sDtw3z5cleUdTzcq3irISwvYRLGdIaLIQ5Sw0ZN1rbR2udp5QqbMhIMtnnXswL6PeVUvnAWmC21lo+t9UwLVu2dHQKQpSZI1uETwONTV43Mm6zuI9SygXwBy5VSXZCCIusFcEAL/VvJkWwBZUYT9dPKRWrlNpv/Nm3qnOvqzJz8zmbli3jg4WwA6XUjUCG1voPk833a63bAT2Nj1FWjpWeftVYeno6H374IXFxcY5ORYhSObIQ3gU0V0pFK6XcgBHA+mv2WQ+MMT4fDmyRbweFqL6cpSu0mUpODJgE3G78cDgG+LBqshYnkjPRIDNGC2HOFg0ZI4BVpgdorU8bf6YDH2Hogm1Ga/2u1jpGax0TEhJSibch7CEzM5PRo0dL92hRIzisEDaO+Z0IbAT+BD7WWh9QSs1USg0x7vZvIEgpdRSYDBS1pCilEoD5wINKqUQLHyyFEKI6qPB4Oq31b1rrM8btBwBPpZR7lWRdxxVOlCUtwkKYqVRDhlLKCbgHk/HBSikXpVSw8bkrMBj4A1HjhISE4OHhITNHixrBoWOEtdZfA19fs22ayfMs4G4rx0bZNTkhhJnUzFxHp1AT2Wo83V3AHq11th1zFUYJyZn4ubsQKF39hSjGeI8qbMhwBpYXNmQAu7XW6zE0ZHxobMhIxlAsF+oFnCpcNcTIHdhoLIKdgc3Ae1XwdoSNKaWIiIiQQljUCLV+siwhhG0kJGcwZ/Px0ncUNqeUaoOhu3T/EvYZB4wDiIiIqKLMaq/45Ayigjxl+Q8hLKhkQ8Y2oOs1264CnW2eqHCIyMhITp486eg0hCiVI8cICyFqiF9PXObFrw5ToDU+7s4W9/H3kO/VrKjUeDqlVCPgc2C01vqYtYvIuDnbySvQnEzJIjpQukULIUR5SYuwqCnkk6sQwiqtNZ/vP8+q2DM0DfbiuZubEOjl5ui0apqi8XQYCt4RwH3X7FM4nu4nTMbTKaUCgK+AqVrrHVWYc52WeDmTvAJNdJBMlCWEEOU1ffp0ZsyY4eg0hCiVFMJCCIty8gp4Z8cJfjieQo8m9ZjQPRJ3F+lEUl6VHE83EWgGTFNKFXY77K+1vlC176Juib+UCSAtwkIIUQGNGjVydApClIkUwkIIMykZucz97hhHkzIYeX04d7YPlbGSlVDR8XRa69nAbLsnKIqJv5SBu4sTDfxkgm4hhCiv8+fPs2zZMu68806uu+46R6cjhFVSCAshijmelMHc745xNSefZ/s24cbIAEenJESVik/OJLKep6yLLYQQFXD16lVefPFFwsLCpBAW1ZoUwkLUYY+s2kdqVp7ZdgW8NrQVUdI1VNQxBVqTkJxBr6aBjk5FCCFqpEaNGqGUkpmjRbUnA/6EqMMsFcEAGqQIFnXShfQcMnMLiA6S338hhKgINzc3wsPDZeZoUe1JISyEEEIYxSdnADJRlhBCVIYsoSRqAimEhahCWmtyc3MdnYYQwor4Sxk4K2gc4OHoVIQQosaKjIwkMTHR0WkIUSIZIyxEJZ0+fZqtW7dy+fLlokdqairPPfcczZs354svvuDZZ58tiuXl5REQEMD3339P+/bt2b59O19++SVhYWGEh4cTFhZGWFgY0dHRuLjY75/o1Zx8u51biJoq/lImDQM8cJOlwoQQosKWLl2Kl5f0rBHVmxTCQpRTRkYGH3/8MV26dKF169bs3buXUaNGFcW9vLwICAjgwQcfpHnz5gQHB3PDDTcQEBBAQEAAHh4eXLhwgQYNGgCwb98+Fi1aRHZ2drHrnD59mvDwcJYsWcK6deto164d7du3p127drRq1Qp394ov7ZKamcus/x2t8PFC1FbxyRl0CPdzdBpCCFGj+fnJfVRUf1IIC1FGf/zxB0uXLuXDDz8kNTWVF198kVmzZtGrVy8OHz5MQEAA/v7+uLm5FTuuR48e9OjRw+p5J06cyOOPP05KSgpnz57lzJkznD17lvr16wPg4uLC2bNn2bJlCzk5OQB4enqSmpqKq6srmzZtIisri/bt2xMREVHqer8Xr+Qwa+MRkq7m4O3mbLFl2N9Dbg2i7knJyOVyZp5MlCWEEJV07NgxXn/9dSZNmiRLKIlqSz7tClEKrTWDBg3i22+/xc3NjeHDh/PYY4/Rs2dPAHx9ffH19a3UNZRSBAYGEhgYSJs2bYrFxo4dy9ixY8nNzeXIkSPs37+fs2fP4urqCsCcOXPYsmULYPgGtl27dvTu3ZvZs2cDkJmZiaenJwCnL2cxc+MRMvMKmHZrc1qF+lQqbyFqk4TCibKCPB2ciRBC1GyZmZksWbKE3r17SyEsqi0phIWw4MCBA3z++ee88MILKKXo2bMnt9xyC2PGjCE4ONghObm6utK6dWtat25dbPtnn33GH3/8wf79+9m3bx/79+/n8OHDRfFOnTqRnp5O05atSfFqSL3GzZg2epAUwUJcIz45E5Clw4QQorIiIiIAZC1hUa1JISyEUVZWFp988glLly5lx44duLm5ce+999K8eXP+8Y9/ODo9q/z9/enevTvdu3e3GH/00UfZ/vNuftj1G2lntpOfm8N/s45xy03vo7Xm/vvvp3nz5rRt25aWLVvSvHnzohZkIeqS+EsZhPq64e3m7OhUhBCiRvPz88Pf31+WUBLVmhTCQgC///47/fv358KFC7Ro0YLXX3/doa2/ttT7nrHEhvRl9Eh3nr85isvnTuHkZJgRNzk5mdjYWNasWUNBQQFg6Kb92muv8fTTT5OWlsaqVato1aoVLVu2JDQ0tNQxyELUVPGXMqU1WAghbCQyMlIKYVGtSSEs6qzDhw+TmJhI3759adWqFbfccgsPP/wwffv2rTXF3g/Hkln8QwJRQV680K8Zfh4uhPq3LIoHBQURFxdHZmYmcXFxxMXFcejQIbp27QoYuoiPHz++aH9/f39atWrFnDlz6N27N8nJyZw4cYKmTZvKDJGiRsvIyedceja9mwU6OhUhhKgVmjRpwtWrVx2dhhBWSSEs6hStNT/++CP/+te/WL9+PS1atODPP//E3d2d//73v45Or9IeWbWP1Kw8s+1JV3LwK2EmaE9PTzp27EjHjh2Lbb/xxhs5ceIEhw4dKiqSDx06VLQ24Hfffcc999wDQEhICE2bNqVZs2bMmDGDJk2akJKSQm5uLiEhIbXmywVROyUYxwfLjNFCCGEbn332mfzfL6o1KYRFnbF161amTJnCrl27CAoK4sUXX+Txxx+vVTdpS0VwSdtL4+TkREREBBEREfTv398s3qNHDz755BOOHTtW9Ni+fXtR/P333+fpp5/G19eXpk2b0rRpUyIjI5k2bRr+/v5cunQJFxcX/P39K5SfELYiM0YLUXZKqQHAQsAZWKa1nnNN3B1YAXQGLgH3aq0TlFJRwJ9AnHHXn7XW443HdAb+A3gCXwNPaq213d+MsJva9PlK1E5SCItaLT09nfz8fAICAkhPT+fy5cu88847jB49uqhVU1RcWFgYw4cPtxrv168fCxcu5NixYxw9epQDBw7wzTff8MorrwAwe/ZsFixYQEBAAJGRkURFRREVFcX8+fNxcnLizJkzuLq6EhwcLP+hCruKT87E38OFep6ujk5FiGpNKeUMvAX0AxKBXUqp9Vrrgya7PQKkaK2bKaVGAHOBe42xY1rr4t2PDN4BHgV+wVAIDwC+sdPbEFVgz549TJs2jfnz59OiRQtHpyOEGSmERa104cIF5s+fz5IlS5g4cSKzZ89m8ODBDB48uGiiqNomMzff0SmYadeuHe3atSu2TWtdVNTefffdNGzYkISEBE6cOMGxY8fYv38/CxYsAGDSpEmsXbsWd3d3GjZsSKNGjWjbti1vvfUWAL/++itOTk40atSI+vXr19q/W2F/8ZcyiA7yki9chChdF+Co1vo4gFJqNTAUMC2EhwLTjc8/BRarEv5xKaXCAD+t9c/G1yuAYUghXKPl5OTw1VdfMX78eCmERbUkhbCoVZKSknjttddYvHgxWVlZDB8+nDvuuAOgVhdJyRk5/HPTMUenUSamn4W6detGt27drO47ceJEevXqRWJiYtHDdE3CJ554gl9//RUwrLMcHh7OggULGDZsmP3egKh1cvMLSLycRaeGMuGbEGXQEDhl8joRuNHaPlrrPKVUKhBkjEUrpX4D0oAXtdY/GPdPvOacDe2Qu6hCkZGRgKwlLKovKYRFrfLEE0+wZs0a7rvvPl566SVatmxZ+kE13MmUTF7ddJQr2dWvRbiyevfuTe/eva3Gly1bxvHjx4sVyqGhoVWXoKgVTl3OIq9AExUo44OFsLOzQITW+pJxTPAXSqk25TmBUmocMA4gIiLCDikKWwkNDcXNzU2WUBLVlhTCokZLTk5m/vz5jB49mhYtWjBz5kymTZvGdddd5+jUqsS+M2m8vuU47i7OzBzYglc2HbU4MZZ/CTNG12SWul4LUV7xlwonypJ5A4Qog9NAY5PXjYzbLO2TqJRyAfyBS8bJr7IBtNaxSqljQAvj/o1KOSfG494F3gWIiYmRybSqMScnJxo3biwtwqLaqp2fjkWtl5KSwhtvvMHChQtJS0ujQYMGtGjRgubNmzs6tSqz7cgl3tlxgob+HjzfrxkhPm78e2R7R6clRI2TkJyJh4sTDfzcHZ2KEDXBLqC5UioaQ7E6Arjvmn3WA2OAn4DhwBattVZKhQDJWut8pVQToDlwXGudrJRKU0p1xTBZ1mjgzSp6P8KOunTpgp+fDDsR1ZMUwqLGmTt3Lv/85z9JTU3lrrvu4uWXX65TrYJaaz7ee5ZP9p6jfbgvT/dpgrebs6PTEqLGik/OICrQEyeZKEuIUhnH/E4ENmJYPmm51vqAUmomsFtrvR74N/ChUuookIyhWAboBcxUSuUCBcB4rXWyMfZ3/lo+6Rtkoqxa4aOPPnJ0CgDk5uaabXNycsLZ2RmtNXl55r3p7B13dnbGycmJgoIC8vPNh7fZO+7i4oJSivz8fAoKCqo87upqWKUhLy8PSyulFcYLCgrsNs9P7Z09qJIOxl+g8+jFHIy/YJf9q+s1qqukpKSi52fOnKFv377s3buXTz/9tE4Vwbn5BSz+4QSf7D1Hn2aBPH9LUymChaiEAq1JuJQp3aKFKAet9dda6xZa66Za61eM26YZi2C01lla67u11s201l0KZ5jWWq/VWrfRWnfUWl+vtd5gcs7dWuu2xnNOlDWEhS15eHjg5uZW7PHMM88AkJmZaRZzc3Nj+vTpgGElEkvxf/3rXwAcO3bMYvzdd98FYO/evRbjhV8S7Nixw2J8wwbDP49vv/3WYnzr1q0AfPLJJxbju3fvBmD58uUW44cOHQJg0aJFFuOnTxtGJ7z66qsW42lpaQA8//zzFuOFxe/EiRPNYv7+/kV/Nw8++KBZvFGjv0ZKfP/99zb+bfiLQ1uEK7oguzH2PIZ16vKBSVrrjbbI6ZFV+4rGWEb36ca0bYmwLRF/DxeL3U7Lu39NuIapsl7DlscUFBRw6NAhduzYwbyVX3L20B6uXkik13NLCWkVg24/inpd3ejQwTZ5VcV7t+U1AO7tFMbwDg1kqRchKul8WjZZeQUyUZYQQtjBpk2beOaZZ1i3bh1RUVFVeu2vv/6a5ORk7r//fmbNmmXW6tilSxfA0PI4e/Zss+N79OgBgLe3t8V4z549AQgMDLQYLzx/WFiYxXiHDh0Aw+zaluKF8920aNHCYrxJkyaAYb4US/HCYrJz584W4yEhIYBhBQ9L8cIu7b1797YYd3c3DCe69dZbixW2hQo/ow4ZMqRYYQuG1uJCd999t9ncPt7e3kXPo6Ojzc5tKw4rhCuzILtSqjWGbjZtgHBgs1Kqhda60tPmWis8bLW9Ll/DWiw/J4ujcQc51MmNVq1asWvXLrp27QqAu289gpp1oEnvu/AIMPyDVU5ONe692/IaAHd3DLMaE0KUXXxyJgDRgdIiLIQQtubk5MS+fftISEio8kJ43rx57N69m7feeguA4ODgolbWp59+mmnTpjFt2rSi/SMjI1m9ejUAEyZMYOrUqcXOd91117F8+XIAxowZw5dfflks3rlzZxYvXgwYirtr4z179mTevHkA3HbbbWbx/v37M2PGDAD69OljFh82bBhTpkyhoKCA7t27m8Xvu+8+nnjiCa5cuUK/fv3M4o8++igPP/wwFy9eZMiQIWbxJ598khEjRpCQkMDIkSPN4lOnTmXo0KEcPHiQRx55xCw+c+ZM+vXrx+7du3niiSfM4q+//jrdu3dn+/btTJkyhS+++KJY/O2336ZTp05s3LiR6dOns2bNGgDuvfdennrqKezBkS3ClVmQfSiwWmudDcQbx6B0wTApg01kXDrLN88OLrZNPWz4+ebX+/lk2Rts/2y5xWPvWh7LL+++hHr4a6vxLa8+TMrR381i6mHFm1/vY8rIv5GRmmwWd3L14I6lO9jwZD9y0pOLcirk7ObOgi9288WEnuRnZ5gd71EvlF/7HuLG6CDQ5v31fRs2p/+s1awdewMYuzSYXqNZ2w7898tt3BhVz+J7i+57L9c/8JzVlsrO414hqusA1j7c2Sx23Vz4JSGFvMC/vvnJTk/hzG/bOPPbNvZ/vJC7lscCsPbhzmbvHQzHF8Yt+bVvCtnZ2fRq2cBi3PT8gNk1fklIIf1KutXz37U8ll9PXLb653Pt+a+9RmH+JR2f8PO3KCvX11pz++23m918TOPt2rXjjz/+sBoPDQ3lwgXzrvRKKQoKCvD19eXKlStmcVdXV3JycnBzc7M4FsfHx4f09HScnJwsjgWpX78+58+ft/q707ZtW/bv3281Pnjw4KL/4IQoq/hLGTgraFzPw9GpCCFErVO4lnBVL6GUm5vLTz/9RE5ODt7e3jg7O+Pr61sU9/LyMpvEy7QVsrS4t7e3WdzL668vVH18fEqM+/r6mo0b9vT0LBZ3c3OzGrc0AZmHh+H/MaWUxXjh+UqLOzs7VyheOKbXxcXFYrywFbi0uKura7F4YcuzPTiyEK7MguwNgZ+vOdbiwusVXW/uYtweq7Ft53I4fq7kMbdXLliZKt74IT7rcpLluJMz287lkJNXcuN2XpZ5kQugXNzYdi6H/Nxsq8fO23IcSht6Y2FQO0Ca8jYcbyeVPXdpx9v7/JW9hj3/bIUQ5uKTM2lczxNXZ5kyQwghbK2wS2xVF8KxsbHk5OSuR6WXAAAgAElEQVRQv359Nm/ebBafNWtWiccXjv+15u233y4x/v7775cYL2x5tmb9+vVWY05OTmzcaH1EqLe3d4nx4ODgEuONGzcuMd6yZcsS4x07diwx3q1btxLjffv2pW/fvlbjtqQcNReBUmo4MEBrPdb4ehRwo9Z6osk+fxj3STS+PoahWJ4O/Ky1Xmnc/m/gG631pyVdMyYmRhcOHLdm+PuGAjjnSioZyeeKxbr6ZdCx8434+vpx+tRJjh+L4+e04t3pfMOicXZ1Iyv1Eh21eTHc5aaerPjzCpmXL5KdZt7i+3+33YiLswtHDx/i7BnD9wSm1wiIaAkYWqxzrqbR1e+vgtjFyZWbevUGYMHmfeRmpBc7t5OzC34Nm3Lo+185m3gCnWcslhVc3/sGnFzc8As3tMamnz9Jfrahy+AdkYbi3c/fnx69/gbA9u+38uH+tGLnd3H3wifUsLTgIK9TZGZmFouH1K/Px+cDAEg7E493SEOcXf/6tmvyDf6EhRm6/B4+fBitNcvirP9+jm1ZvGUwMDCQkJAQnl1/yOoxrw1pRW5uLseP/1VwlvUaDRo0wN/fn8zMTKZtsn5D/3s7d7y9valfv37RttJysqSkYz596HqrMVE+SqlYrXWMo/OwpbLc64TB2NX76NTQj8d7Rjk6FSHsSu51wlEaNGjA4MGDWbZsWZVdc+7cuUydOpXffvuNjh07Vtl1RfVQ1vudI1uEK7wgexmPrRQ3H3/cfIoP/J5tWnx0igB6FBXO1/LwD2L2Q/0sxlb8uQfPgBA8jWNeTd0WYxj4fmunv1qvLV3DKygMr6Cw4jmZeG+vlRZn4I//jOdKRjZHE5M5fDKJo6eS2EPx9+ob+tf1H3/Q/BrRd97GhlTrreYP3zvU4vaPje+lsOA21a1t07/Of5PhprWshJb5Ad06WY1ZUzgzbIsGfx1b/mt4AdYL4ddeeoZNmzYxePBgxo8fT//+/cuUkxCiaqVk5HI5M48o+TcohBB2M3jwYLtOeGRJXFwcrVu3liJYlMiRfcGKFmRXSrlhmPzq2n4AhQuyg8mC7MbtI5RS7sYF3ZsDv9oiKX8Py98N2Gp7dbmGj5c7HVuEcc8t7fjHQ31wV5ZbRfNzcog/Y95yba+8asM13nzzTZ5++ml27tzJwIEDadq0KW7afMxsZa4hhKi8+EuGHjUyUZYQQtjPsmXLeOGFF6r0mj179qRr164W188VoojW2mEPYBBwGDgGvGDcNhMYYnzuAXwCHMVQ6DYxOfYF43FxwMCyXK9z5866rA4cP6+vGz5PK68Q/eqrr5Zp/+tHvakPHD9frmuU5xh7X+P+l1brgL4ztEePl3RA3xm6y4Nv6eB+s3RA3xn6leVbdGZWTpmvK7TOysrSq1ev1n369NFz587VWmudmZmpt2zZogsKChycnSgE7Nb2v9cNMN6rjgJTLcTdgTXG+C9AlEnseeP2OODWslyvvPe68t5XaotP957Rdy2P1Vez8xydihB2VxX3uqp+lOdeJ+qWW2+9Vbdp08bRaQgHKev9zmFjhB2hImNJZs6cycCBA7nhhhvslFX1cTUzh06j3iTxQiqNQwPYs2Iil69kMWXxt6zd8gdNGwUy/6nb6H9jc0enWuNorVFKsWrVKu677z5uueUW5s+fT7t27RydWp1n73FzxqXiDmOyVBwwUpssFaeU+jvQXms93rhU3B1a68Kl4lZhmBU/HNgMlLpUXFnudbV9vW1bHyNETSdjhIWjfPTRR0ycOJE///yT0NBQu1/v7bffZtKkSTz22GNFSyeJuqWs9zuZJrMU06ZNqxNFMIC3pxtfvDaK66Lq8/m8B/D2dKNhiB8rZ9zDV2+MwUkphj7zISNeXM2p86lFxx2Mv0Dn0Ys5GF/yTNp1WeGSP3fccQcLFy4kNjaWjh07MmHCBC5evOjg7ISdFS0Vp7XOAQqXijM1FPjA+PxT4OZrl4rTWsdjaBnuYoukavt627Y+RgghRMV4e3uTkpLCyZNWVlSxsZUrV5Kfn0+/fpbn6hGiSFmajWvLo6JdaBISEvSkSZP0lStXKnR8bZGVnavnfrBN17t5pg7qN0u/vnK7Tk69qpvf9br27DlNtxj+L30lI9vRadYIly5d0pMmTdLOzs76pptucnQ6NUJBQYG+evWqPnv2rI6Li9PHjx8viq1cuVIvXLhQz5o1S//xxx/lOi927i6IYX6DZSavRwGLr9nnD6CRyetjQDCwGHjAZPu/geGlXbMs97q7lscWPQCLj3X7z+n7nnzJavyu5bE6osewEuNBLTuXeH6fesEWY8rZVd+1PFa7+QZajDu5uOp1+89pZ3cvi3E330C9bv85jXKyGPeuH2l471bijRs3Lvz9sPjo1q1bifEXXnihxLjJ75/ErcRfeOGFEuPdunUrMd64ceMS4z4+PiXGXV1dLf9uKqW11lopZTHu6upa4nvz8fEpMV7W373yQrpGCwf57bffNKA/+eQTu18rLy9Pu7u7a0AnJyfb/Xqieirr/U5m3SmDxMREFi1aRGRkJJMnT3Z0Og7j7ubCc6P/xr392vPMom94cckm5q7YTnZOHlprLiRfYfycL/hwxj2OTrXaCwwMZOHChYwfP570dMMyV2lpaWzbto3bb7+9qAW5NtFak5aWxuXLl7l8+TJa66LZHFeuXMnRo0eLYpcvX6Zp06ZF6/h17tyZffv2FVt8fsiQIaxbtw6AZ555hnPnDMudRUVF0aZNmyp+d45X0TXTS7Ji12mOp9up45ByYsWu02hXb8D6LPdWD3d2ZcWukhcLWLHrNMrJGZ1veV10IYQQ9hcZGQnAyZMn7T6cZe/evWRnZ3P99ddTr149G2QvajMZI1xGN998MwcOHCA+Ph5PT08bZ1YzPbPwa95e+zOmv0JeHq7Mf2oQY27r7LjEaqh58+YxZcoU+vbtyxtvvEH79tVzrGJBQQGpqalcvny5aDmErVu3sn//fpKTk4seHh4eRWsG3nHHHWzYsKHY7I1t27Zl//79AHTv3p2dO3fi5+dHQEAAAQEBdO3alaVLlwIwZ84c0tPT8fPzw8/PD39/f6Kjo7npppsAOHXqFF5eXvj6+uLm5kZ5VMEY4ZuA6VrrW42vnwfQWv/TZJ+Nxn1+Mi4Vdw4IAaaa7mu6X0nXLM+a6ZZ8+EAHs22jVv5erv0rckxVXKO0Y2SNblFbyRhh4Shaa/z8/Hj44Yc53XGM1f2s3X9L+v/q2mPmz5/P008/TWJiIg0bNqxYwqLGqwnrCNco06ZNo3fv3rz33ntMmjTJ0elUCx9/t59rv0fJyMrlhXc2SSFcAZMnT8bb25tp06bRqVMnHnnkEWbNmlUlE0skJSVx4sQJkpKSSEpK4uLFiyQlJTFr1iyUUsybN48PPviApKQkLl26RH5+Pl5eXly9ehWA5cuXs3LlSgACAgIIDAwstmbgoEGDaNOmDfXq1aNevXoEBATQoEGDovjGjRvx9PTE2dnZYn5Tp04tMf/GjRuXGHewoqXiMKx3PgK475p9CpeK+wmTpeKUUuuBj5RS8zFMlmWzpeJK4ulq+e/BVvtX12sIIYSwPaUUjz76KB06dOB0CR10nvrsoPVgGfn4+DBs2DApgkWZSCFcRn/729/o1asXc+fOZdy4cXh4eDg6JYeb9dgtTF7wNRlZxdfITc/I5tX3tzLp3m74eLk7KLuax8XFhccff5z77ruPmTNnsnjxYi5dusTatWvJyckp6kJd2G1aKYWPjw+urq7k5OSQmZlZFCsoKCApKYnGjRvj7u5ObGwsX375JefPn+fChQtFP3fs2EFwcDALFy5k9uzZxfJxcnLi2Wefxd/fH39/f1q1akVwcDDBwcGEhIQQFBSE1obZsN944w0WLFhAQECAxWL20UcfLfG9+/j42OKPsFrSWucppSYCGwFnYLnW+oBSaiaGMSzrMYz9/VApdRRIxlAsY9zvY+AgkAc8rkuZMbqs/D1crHY1s8X+1fUaFT1GCCFExc2fPx+ADSW07jauZ/mzdWJqVpmvc/z4cfbu3UteXh4uLnJPFyUrtWu0UqopkKi1zlZK9QbaAyu01perID+bqmwXmq1bt7J69WrmzJkj4w6MHpi2hq92xJGVk4eHmwt/uz4aT3dXvvj+IKGBPvzjwd48dHtnXF2kZaa84uLicHZ2plmzZmzYsIEhQ4aY7bN161Z69+5dtCzTtfbs2UOnTp1YunQp48ePJzAwkNDQUOrXr09oaCgLFy6kQYMGHDhwgKNHjxISElJU7AYEBODkVPsnlq/r3QXHzPiET7/7g10fPE7r6Pp2zkwI4Sh1/V4nHO/q1auM+TjOarwiXaNnDGxOmwa+AGRkZNCrVy88PT354YcfKpesqNFs2TV6LRCjlGoGvAusAz4CBlUuxZqnT58+9OnTx9FpVCtLn7+jaO3h+oE+/HfmvXh7uvHLgVO88M7/eHL+lyz6eCczxt3Cnb3bFLVYHoy/wKiXP+bDGffIh28rWrZsWfS8bdu2LFq0qHDm4KKfzZo1A6BTp05FE0sVttIGBwcXdRl+8MEHeeihh6yOn23Tpk2dnGBKQNumoXy8eT+RDQIcnYoQQoha6o033mDy5MkMfedHXNxtN9fOy98coXt0PUbf0JD33lxEbGwszzzzjM3OL2q3shTCBcaufXcAb2qt31RK/WbvxKqzn376idTUVAYMGODoVByucO3hwqLW29NQaN3YpjGb3nyYb386zEtLN/HAtI/pfF1DXp3Qn86tGjLs2Q9JvJDGHc+tZM+KiUXHCcuio6N54oknrMZbtWpFq1atrMbd3aWLurAsPNgPgDMX02geEezgbIQQQtRGhfOCOKVfAPdIs3hFh7P0bxXMF/vPE3sqlV2rPgVg8ODBNspa1HZlKYRzlVIjMUzkcrtxm6v9UqretNZMnjyZs2fPcuTIEVxd6+wfRZHW0fWJXTHRbLtSioHdWtL/xuZ8tPF3Zvz7O26d9D6hgT6kpGfKkktCVAPhwYYuZWcvpUshLIQQwi4Kl/Ubc50ba5IVfZsHM/amsk10aW1ZpUK9mwWx/OeTrPpzP8rJmffivXjzuHl3alst0VTRY+Qa9j+mvMoyAPAh4CbgFa11vHHm0w9tcvUaSCnFtGnTOHHiBCtWrHB0OjWCs7MTowZ1Yv9HT3JnnzacT75CTq5hvp+snDy+2hHHB1/FOjhLIeqm8BBDi/Dpi2kOzkQIUVMopQYopeKUUkeVUmbT+iul3JVSa4zxX5RSUcbt/ZRSsUqp/caffU2O2WY8517jQ8ZN1SKFawn/uDeOnHxN9ya2m2sn1NedwQ2yKMjNoXn3AWRpy/PSWCqqSovZ8hi5hv2PKa+ytAj301oXrRdkLIbLPn1bLTRgwABiYmJ45ZVXGD16tLQKl5Gnuys/7E0w256ZnctT87+iacMguneILBpHbErGFAthH4WF8BkphIUQZaCUcgbeAvoBicAupdR6rbXp2jePACla62ZKqRHAXOBeIAm4XWt9RinVFsNs+qbr3NyvtZbZr2qhsLAwXFxc2H3gCC2a96NlfW+bnv/7778H4Mv3F/H8duvz+W45cqnc566KY+QajlGWQngMsPCabQ9a2FZnFLYKDxkyhP/+9788+OCDjk6pxrC05JKzs6JAa/o9sZzo8Ho8MKAj9w/oSGSY4dvCq5k5MqZYCDvx9XLH18udM0npjk5FCFEzdAGOaq2PAyilVgNDMSzzVmgoMN34/FNgsVJKaa1N55g5AHgqpdy11tn2T1s4krOzM1P/8SLbrwbTLboeThYaPSojIiKCSZMm0TQ6Crbvtbrf2z+eKPe5q+IYuYZjWC2EjeOC7wOilVLrTUK+GNa5rNMGDx5M3759yc6We3d5jLmtM5t+OVpsyaXBPVqxZOow1n1/kA+/+Y3Z729j1vKt9OoUxQMDOvHVjkNcTLkqY4qFsJPwEF9pERZClFVD4JTJ60TgRmv7GCdcTQWCMLQIF7oL2HNNEfy+Uiofw4ols3Vpa3yKGqXPAxP588cTdI+2/RKkX331FWvWrClaQcOad+62vELGhE8O2P0YuYbtjrGVklqEdwJngWDA9LcqHdhnz6RqAqUUmzdvttiNV5Ts2iWXlkwdhrenG/cN6Mh9Azpy4txlVm3cy8pv9zLun58XOzYrJ4+vdxrGFI+5rXOJ15Hu1EKUTXiwH2eSpBAWQlQNpVQbDN2l+5tsvl9rfVop5YuhEB4FmE3GopQaB4yDvyZgEjXD94fO4nHlHE2DO9n0vGfOnGHDhg306NEDF5eSO7uG+JR/FY2qOEau4RhWJ8vSWp/QWm/TWt+ktf7e5LFHa227Uco1mFKKgoICNmzYQH5+vqPTqTEKl1y6Lqo+n897wKybc2SDAKaO6c3+j54kwNfD7PiMrFwm/etLZiz7jnXbD3LibArXfmlc2J36z4SL3PHcSq5m5tj1PQlRk4UF+3LmonSNFkKUyWnAdLrfRsZtFvdRSrkA/sAl4+tGwOfAaK31scIDtNanjT/TgY8wdME2o7V+V2sdo7WOCQkJsckbEvaXmpXLJ+8t4KOnbrf5Z+bly5dz9uxZOnUyFNjWlmIqbYkmex8j17D/MeWlSut1opS6E8O3dvUBZXxorbWfzbKoIjExMXr3btvOwfDll19y++23s2rVKkaMGGHTcwv44KtY8zHFTorQQB/OJV+hoMDw+1vP15MOzRvQoUUYHZqHsWrj7/ywN6FY9+vSulNLC3LdpJSK1VrHODoPWyrvve6lJZtYsHoHqVum4eRUlsUEhBA1ja3udcbC9jBwM4aCdxdwn9b6gMk+jwPttNbjjZNl3am1vkcpFQB8D8zQWn92zTkDtNZJSilXYBWwWWu9pKRc7PG5TtjHxkMX+cfcRez5z2wSEhKKZpG2hS5durBr1y5+//132re3zbI6omYr6/2uLJ945gFDtNb+Wms/rbVvTSyC7WXQoEG0bt2aWbNmUVBQ4Oh0ap0xt3Vm4E0t8HAzfPvj4ebCHb3bcOzzZ0n634tsXzqORU/fzp192pCemcOSz37l4Vlr2fTrUbJyDB0XsnLyWLf9IFPe/JYzSWlmrccgLciibgsP8SMvv4CLlzMcnYoQopoz9gqciGHG5z+Bj7XWB5RSM5VSQ4y7/RsIUkodBSYDhUssTQSaAdOuWSbJHdiolNoH7MVQYL9Xde9K2NuO+JSiruwnT5602Xm11uzfvx93d3fatWtns/OKuqEsbcvntdZ/2j2TGsrJyYmXXnqJkSNH8tlnnzF8+HBHp1TrWBpTDIblmG5o3YgbWjcq2jc3L5/IofNIScssdo7cvAIWfbyTRR/vpJ6vJ62b1KdNdH1aNwmlTXR9Fq3ZWaEJuaQVWdQG4SG+gGEJpdBAHwdnI4So7rTWXwNfX7NtmsnzLOBuC8fNBmZbOW3JE3+IGis5I4c/z12hd6dWfAGcOHGCnj172uTccXFxZGVlMW3aNJm3R5RbWVqEdxsXRR+plLqz8GH3zGqQu+++m1atWjFz5kwZK2wHpY0pNuXq4sw//94fL4/iazt7ursy+b4ezH/qNu7q24aCAs2azft5av6X9HtiORt+PFSsBXn99j+ZuvhbTp6/bLWlX1qRRW0RHmxcS1gmzBJCCGFjP8VfRgO33WiYHdiWLcI//PADAPfff7/NzinqjrK0CPsBGRSf2U8Dn1neve5xdnZm+vTp/OMf/yAxMdGm4x6EQevo+sSumFimfS0t0XRb95a8MqF/sf201py+mEbMmLdIvZJVLJaTl8/CNTtZuGYnXh6uNG8cRPPGwbSICKZFZDAtGgczd8X3sqyTqBXCQ4yFsEyYJYQQwsZ2xKcQHehJs7BA3n77bbp27WqzcwcFBTF27Fj57C0qpNRCWGv9UFUkUtPdc889DBs2DHd3x00BLv5irTu1KaUUjer7M3firWYTcnm6uzDhrq5Eh9cj7kQSR04lEXvoNJ9tO1A0QZeprJw8vtpRtmWdhKhuQut54+SkpEVYCCGETZ1Pz+bwxas8EBMOwIQJE2x6/lWrVvHLL7/g5ma9t6AQ1pTaNVop1UIp9Z1S6g/j6/ZKqRftn1rNopTC3d2djIwMduzY4eh06rzydKe2NCHXbd1b8cqE/owdegOvTRrIF6+N4uCa/+PS/15k9weP4+dt/oVHZrZhWafVm/aRdjXLLA6GMcWdRy/mYPwF27xRIWzAxcWZ+vV8OJskLcJCCCFsZ2d8CgDdouoBcPr0aZt9Tj5+/DgbNmwgJiZGxgeLCinLGOH3gOeBXACt9T5A1gmyYtKkSQwcOJALF6TQcbTC7tRlmcRq6fN3EFLPG6Ww2oIM4OHuSpsmocx7YoDZOGQXZyc83Fx4aOanNL59LndNWcmHX/9GcpphJl4ZUyyqs/AQX85clBZhIYQQtrMjPoUWId7U9zU0ILz++uv079/f4goe5bVixQqys7Pp0sXiktNClKoshbCX1vrXa7bl2SOZ2uDZZ58lIyODGTNmODoVUQ7laUEGy63Iw/7WmrNfP8+Wt8cy/s4b2X/sPOP++TmRQ+YxePIH9H9iOReuGVMsRHURHuwnXaOFEELYzOnULBKSM+keXa9oW2RkJBkZGSQnJ1f6/F9++SUgE2WJiitLIZyklGqKYYIslFLDgbN2zaoGa9myJY899hhLly7l0KFDjk5HlEN5WpDBciuyk5MTN7WLYO7EAcR9Mpkf33uMp0Z0Z9+Rc+yJO0O2yczUhWOKhagOGob4yWRZQgghbGZnfAoKuCk6oGhb4aRWJ06cqNS5tdYcOXKEevXq0bhx40qdS9RdZSmEHweWAq2UUqeBpwDbjnSvZV5++WW8vLyYMmWKo1MRdlRaK7JSis6tGjJrfD+wMHQlMzuXp+Z/xe4/T9uki5AQlREe7EtKeiaZ2bml7yyEEEKUQGvNj8eTad3Ah0Cvvz4f2aoQTkhIIC0tjZkzZ1bqPKJuK7UQ1lof11rfAoQArbTWPbTWCXbPrAarX78+U6dOJSUlhYyMDEenI+yorK3Isx67xWxMsbOTIj+/gJ7jlnLjQ2/z1qc/F40nNiUTbImq8NcSStI9WgghROWcSMnkdGp2sW7RABEREUDl1xI+efIkjRo1onfv3pU6j6jbrBbCSqkHjD8nK6UmA48Bj5q8rjClVKBSapNS6ojxZz0r+40x7nNEKTXGZPsrSqlTSqkrlcnDnp577jm+//57vLy8HJ2KqAYsjSm+o3cbTn05lTefuR1XV2eeWfg10cNeY/T0T9iy+xgFBQUywZaoMlIICyGEsJUd8Sk4KegaVfwjflBQEGvXrmXYMMuTkpbVzp07CQoKolmzZpU6j6jbSmoR9jb+9LXyqIypwHda6+bAd8bXxSilAoGXgRuBLsDLJgXzBuO2asvFxQWlFGfOnOHnn392dDqiGrA0ptjfx4OxQ29gx3vj+fX9vzN26A1s/vUot/3fB7QesYBej73LhZQrMsGWsLuwIMNt/cwlGScshBCi4rTW7IxPoX24H34eLsViSinuvPPOoi7SFbV582YAPDw8KnUeUbdZLYS11kuNP2dYelTyukOBD4zPPwAsfS10K7BJa52stU4BNgEDjDn9rLWuERN23XvvvYwYMYKsLMvryoq6o7Qxxe2aNeBfTw7i+OfP8MHLd+Ph5sLB+Atk5+QDhgm2vt4pE2wJ+5AWYSGEELZwLCmD8+k5Zt2iC/3222+sW7euwuePi4tjy5YtNGnSpMLnEALKMEZYKfWBUirA5HU9pdTySl431KSQPQeEWtinIXDK5HWicVu5KKXGKaV2K6V2X7x4sfyZVtL06dM5ceIEixcvrvJri+qnLGOKPdxdueeWdiSnZZrFMrJyeWnpZnumKOooP293vD3dZOZoIYQQlbIjPgUXJ0WXCH+L8SVLlvDoo49W+Pz/+c9/AOjTp0+FzyEElG3W6PZa68uFL4yts51KO0gptVkp9YeFx1DT/bRhuly7TZmrtX5Xax2jtY4JCQmx12Wsuvnmmxk4cCCvvPIKly5dqvLri5rL0gRbnu6uzB7fz0EZidpMKUV4sK+sJSyEEKLCCrRmR3wKnRr54e3uYnGfiIgILl68WOEJZb/55hsARo0aVeE8hYCyFcJOppNZGcfuWv7NNqG1vkVr3dbCYx1wXikVZjxfGGBpOtzTgOnCYI2M22qcefPmkZaWxuzZsx2diqhBrp1gC6BBkA+jB13vwKxEbRYe4iddo4UQQlTYofNXSc7ItdotGv5aQunUqVNW9ynJ+fPnadGiBQEBAaXvLEQJylII/wv4SSk1Syk1G9gJzKvkddcDhbNAjwEsDRTYCPQ3dsWuB/Q3bqtx2rZty9ixY7ly5YqsFyvKxXSCLV8vd+LPpLBl9zFHpyVqqfBgP84kSddoIYQQFbMzPgU3Z0Xnxpa7RcNfSyhVZC3hU6dOce7cOSZMmFDhHIUoVJZ1hFcAdwLnMYznvVNr/WElrzsH6KeUOgLcYnyNUipGKbXMeN1kYBawy/iYadyGUmqeUioR8FJKJSqlplcyH7t75513eO+991BKOToVUYOYTrD17aIHad44iAlz15Geke3o1EQtFB7iy9mkdAoKChydihBCiBomv0DzU0IKMY398XR1trpfYYtwRQrhxMRERo0aRf/+/SucpxCFSlpH2M/4MxBDAfyR8XHOuK3CtNaXtNY3a62bG7tQJxu379ZajzXZb7nWupnx8b7J9ue01o201k7Gn9Mrk09VcHIy/FH/9ttv/P777w7ORtQkhRNsXd+yIUumDuPU+VReWrLJ0WmJWigs2I/cvHwupZpP1CaEEEKU5MC5dFKz8ujepOQyoWHDhvzyyy/cc8895b7G2rVr+fjjj2XGaGETJbUIf2T8GQvsNnkUvhbllJuby+DBg5kwYYJ0kRYV0q19JI8P78rSz3/l+z3HHZ2OqGXCg41rCcuEWUIIIYNIOs0AACAASURBVMppx/EUPF2d6NTQr8T9XFxc6NKlC/7+1rtPW/PVV19x0003yfrBwiZKKoTnGH9ep7VuYvKI1lrL1zAV4OrqyowZM/jpp5/47LPPHJ2OqKFmjLuZJg0DGT9nHVczcxydjiiBUipQKbVJKXXE+NPi7CFKqTHGfY4opcYYt3kppb5SSh1SSh1QSs2xdKwtyVrCQoiyUEoNUErFKaWOKqWmWoi7K6XWGOO/KKWiTGLPG7fHKaVuLes5RfWWm1/Azycuc0NEAG4upU9B9M033xQtg1RWBw4c4NChQ3h6elYwSyGKK2n254VAZwyTY8k0tTby0EMPsWDBAqZMmcLtt9+Om5ubo1MSNYyXhxtLpw6j3xPLeWnpJuY/dZujUxLWTQW+01rPMX6wmwpMMd3BONTkZSAGw1JysUqp9UA28LrWeqtSyg34Tik1UGv9jb2SDQ82FsIyYZYQwgqllDPwFtAPSAR2KaXWa60Pmuz2CJCitW6mlBoBzAXuVUq1BkYAbYBwYLNSqoXxmNLOWSGPrNpHalae2XZ/Dxf+PbJ9pfevqmNqyjW2H0tm+7HkEo8BWLlyJTt37uTBBx+0us+1li1bBsBdd91V5mOEKElJhXCuUupdoJFSatG1Qa31JPulVXs5Ozvz2muvMWjQIJYsWcKkSfLHKMqvR8coJtx1I++s/YU7e7ehR8coR6ckLBsK9DY+/wDYxjWFMHArsMlkMsBNwACt9SpgK4DWOkcptQfDMnJ20yDIB6WUtAgLIUrSBTiqtT4OoJRajeFeZ1q0DgWmG59/CixWhtlChwKrtdbZQLxS6qjxfJThnBVSWKTl52aTdvqvIUUpwLrv0glr2Iig4BAyrl7l+NHDJBxKKHa8d3A4bj7+JKWksu67HWbnj4iKJjUrj5yMdK5eSCwWSwF+TwjFx8eXyynJJJ78a3Kowuv4hkXh4u5JdloKGcnnivIq1Lzldbh7eHD+/DkyUy6aXb+gUTOOJWVw/txZLp4/VyyWcCiBgIgWKCfn/2/vvsOjqNYHjn/fkISE0EkoAaQrXRBEQeyAIAgiysUCqNeCol7FcgHFAihg/92LIFyFC9hRQaToBcQuCghK6FKkQyAoECD1/f2xG0zChmzCbmbL+3meebIzc86cd7Jwsu/OzDkcS9lL2uFDec4doPm5rQHYtWM729avzlNfRKhYpzF/nshkyfI1HP7zjzz7IyMj+fNEFABH9+0g4/jRPOeelBRB8+bNAdiwYQNHj/61Pyoqiu3btzNjxgxq1apF9erVSUlJYenSpSfnF46JiaFGjRrUqVOHSpUqMWfOHABuuummU34PxhSLqnpcgHhc39r9jmuKozxLQfUCeWnTpo0GguzsbO3Vq5e+9NJLTodigtjRY2na+IaXtenfXtHU42lOhxO0gOXqpz4H+CPXa8m9nmv7I8ATudZHAI/kK1MR2ALU96bdM+nr6vYcp4PGzCp2fWNMYPJVXwdcD7yRa70/MD5fmSSgVq71ze7PleOBW3Jtf9N9vEKP6Wnxpq/rM2WF9pmyQrs8+5Hiuusm3yLaZ8oKbdHvYc/7xbW/docenvdHRGqfKSu0yjltPe6PiCqtfaas0LLV63rcHxlbVvtMWaGlKyZ43F+6Qrz2mbJCI8uU87i/bPV62mfKCo2IKu1xf+WzW2ufKSuUiEiP+2t36O7aL+Jxf+Neg1z7Pf7u0PYP/t9p9+f691fgUq6c53PLWapWrVrs4wfD/pdeeum0+3v0KODfnnt/8+bNT7u/oN+fiKiqatmyZT3uj4qKUlXVqKgoz//2ypZVVVUp4N9O1apVT3vuzZs3P+3+Hj16nNxfVHjZ353uJv5HVfU94DVVnZZ/OU09UwgRYdasWQwZMsTpUEwQi4uN5vWhvdiyK4Wn/7PY6XDClogsEpEkD0uv3OVydfZFPX4k8C7wL3VfLSmg3F0islxElicnn3rVwFs14svZYFnGmIBV3L4utlLVAvYoQ6+szyX1PQ7hAOL6qFy6fBXP+7OzXPvLVfR89KwMhl5Zn2jNKODwruNHxcR5rp92jKFX1qdUVGnP9UtFMPTK+lDgtHeuaTsjIj3fBFomNdlVv8h/nfI3U7zpQZ944gnGjx9/2v3/93//V9yojDktcX0287BDZDXQElihqiHxjHDbtm11+fLAGfBaVZk9ezb16tWjVatWTodjgtQ/XvqU/3yynEXjb6dDyzpOhxN0RGSFqrb107E3AJep6h4RqQF8qarn5Ctzo7vM3e71Se5y77rXpwBHtQiPo5xJX3f90LfZvvcPfvrv4GLVN8YEJl/1dSLSHnhaVa9yrw8DUNUxucp87i7zg/vLvL1AAq5xEk6WzSnnrnbaY3riTV93/dSfC9z34W2nfrwtavmSqhMqbeTo1KkTR44c4ccffyywTH73338/06ZN488//0SKmXib8OBtf3e6K8Kf4brFv6WIHBaRI7l/+izSMJaamsqgQYO46667yMrKcjocE6RG39OF2tUqMGjsbI6nef7G2ThmDq7HSXD//MRDmc+BLiJSyT2qdBf3NkRkNFABeLAEYgVcA2bZYFnGmNNYBjQSkXrugfz64errcsvd910PfOG+K2YO0M89qnQ9oBHwk5fHNCFkx44dJ58f9laDBg04cuQIBw4c8FNUJtwUmAir6qOqWhGYp6rlVbVc7p8lGGPIKlu2LK+88grLli1j0qRJTodjglS5MqWZ+M9ebNpxkJFvfOF0OCavsUBnEdkEdHKvIyJtReQNAHUNkjUK1wfBZcBIVU0RkVrA40BT4GcRWSUid/g74MSEchz88xgn7EsVY4wHqpoJ3IfrC7t1wAequkZERopIT3exN4Eq7sGwhvDXleA1wAe4BsH6DBisqlkFHdMX8VaI8XxLsK+2l1SdUGkjx/r16/n3v/992jL5NWjQAIDNmzcXqZ4xBSnw1ug8hUTqAI1UdZGIxAKRqhp0lwwC7dZocN0e3aVLF3766SfWr19PjRo1nA7JBKn7XpjD1Lkr+GLCHZQrU5r+T33AjGf60rReQc9FGfDvrdFOOZO+bvq8n7l77GzWvv8g9RIr+zgyY4xTrK8zwW7t2rVccsklzJgxg27dujkdjglgvrg1OudAd+Ia+j7nkmUtYPaZhWdyiAgTJkwgLS2Nhx56yOlwTBB77t4uJMaX585nP6bXIzNYty2Z3o+9RerxdKdDM0EkMcF1w88um0LJGGOMH8ycOZOBAwdy/PjxItVr0qQJBw4csCTY+EyhiTAwGLgIOAygqpsAu8TkQ40aNeL555/nmmuucToUE8TKx8Uw4bGebNpxkD0HD6Oq7E85yqCx9r2V8V5OIrw7Oehu+jHGGBMEFi9ezKeffkpMTEyR6tkAWcbXvEmE01T15CUl9+h/ZzrIusnngQce4Oabb3Y6DBPkdh84TGSpCLKyXP9FT6RnMv/7DUybt8LhyEywqBFfDsCmUDLGGOMXa9asoXnz5sVKbEeOHMndd9/th6hMOPImEf5KRIYDsSLSGZgJfOrfsMKTqvLKK6/w7LPPOh2KCVIjJi0iMyvvXILHTmQwYtIihyIywaZi2RhiS0exx0aONsYY42OqSlJSEs2aNStW/a1btzJ37lwfR2XClTeJ8FAgGVgN3A3MB57wZ1DhSkRYtWoVzzzzDGvXrnU6HBOERt3diTIxUXm2lYmJYvSgzg5FZIKNiJCYUI7d9oywMcYYH9u9ezd//PFHkadOytGgQQN2795d5OeLjfGk0ERYVbOBacAzuCY9/696M9S0KZYXX3yRsmXLcs8992C/ZlNUA7u3oVv7s4mJ/mvagotb1WXA1QVPam9Mfq65hC0RNsYY41sHDhygadOmtGzZslj169evD7iuDBtzprwZNfoyYBPwGjAB2Cgil/g5rrCVkJDA888/z9dff81///tfp8MxQWjSsN4kVIpDBETgz6MnyM7OLryiMW6JCeVtsCxjjDE+d+6557JmzRouvvjiYtW3uYSNL3lza/RLQBdVvVRVLwGuAl7xb1jh7fbbb+eiiy7iscceIzU11elwTJCJi41m9gv9aVK3KiNuv4KlSTt4c44NlmW8lxhfjj0Hj9hdKcYYYwJKw4YNOffcc20EaeMTkYUXIUpVN+SsqOpGEYk6XQVzZiIiIvjPf/7DgQMHiIuLczocE4Sa1qvKiun3oap8s2obj0/8H1dfdA413VPjGHM6iQnlSUvP5OCfx4ivaH2QMcYY3+jRowdNmzbl+eefL1b9KlWqsGrVKh9HZcKVN1eEl4vIGyJymXv5D7Dc34GFuyZNmpy8bSQtLc3haEywEhHGP9qTzKxsHnx5rl3hM17JmULJRo42xhjjK9nZ2SxZsoT09PTCCxtTArxJhO8B1gIPuJe17m2mBIwZM4a2bdtap2GKrX7Nyoz4++XM/XY9Hy1Z43Q4JggkxrvuHLABs4wxxvjKtm3bOHbsWLFHjM7x+OOPF/sZY2Ny8yYRjgT+T1WvU9XrgH8BpfwblsnRqlUrkpKSeOGFF5wOxQSx+29oz3nnJPLwq/NIOXzM6XBMgEt030JvA2YZY4zxlTVrXF/Gn2kinJmZyU8//URWVpYvwjJhzJtEeDEQm2s9Fljkn3BMft26deOGG25g9OjRNkKeKbbIyFJMHHotKYePM3T8506HYwJcjSplAbsibIwxxneSkpIAaNq06Rkdp0GDBqSnp7Nr1y5fhGXCmDeJcIyqHs1Zcb8u47+QTH6vvvoqUVFR3HvvvfaMpym2lg2rM+SmjsxYsJLFy35zOhwTwKKjIqlaKY7dyZYIG2OM8Y3atWtz4403Ur78mQ3cmTOF0pYtW3wRlglj3iTCqSJyXs6KiLQBjvsvJJNfYmIizz77LF999RXr1693OhwTxIYNvJRGtatw3wufknrcnjs3BUuML2+JsDHGGJ+55ZZbeOedd874OPXr1wdsLmFz5rxJhB8EZorINyLyLfA+cJ9/wzL53Xvvvaxbt44mTZo4HYoJYjGlo5jwWC+27TnEyDe/cDocE8ASE8qx20aNNsYY4wPZ2dk+mwWldu3a9OzZk2rVqvnkeCZ8FZoIq+oyoDGukaIHAU1UdYW/AzN5lSpVinr16gHwzTff2C3Sptg6tqrLnb3OZ/zMH1i+zp6vMZ7ViC9v0ycZY4zxiY0bNxIXF8esWbPO+FiRkZF88skn9OjRwweRmXDmzRVhVDVDVZPcS4a/gzIFmzVrFpdccolPbi0x4WvUoM5Uq1yWe8bOJiPTRl00p0qML0fyH6mkpWc6HYoxxpggl5SURFZWFnXq1PHZMTMz7e+TOTNeJcK+JiKVRWShiGxy/6xUQLmB7jKbRGSge1sZEZknIutFZI2IjC3Z6J3Vs2dPOnTowODBg9m+fbvT4ZggVaFsDP96+BqStuzj5Xe+Ze3W/bQZMJ61W/c7HZoJEDlTKO09aFeFjTEu/vr8JiK3ikiyiKxyL3eU1DmZkpGUlISI+OwRv2HDhlG9enWfHMuEL0cSYWAosFhVG+Ganmlo/gIiUhl4CrgAaAc8lavDfVFVGwOtgYtEpFvJhO28UqVKMWPGDLKysrj11lvJzs52OiQTpHp0bEyfK5rz7NQl9HhoGuu2JdP7sbdsEC0D/JUI77IBs4wxf/Hn57f3VbWVe3nDr2dhSlxSUhINGzYkNja28MJeqFKlCgcPHuTQoUM+OZ4JT04lwr2Aae7X04BrPZS5ClioqimqeghYCHRV1WOqugRAVdOBn4FaJRBzwKhfvz6vvvoqS5Ys4dVXX3U6HBPEXvrH1QDsTTmCqrI/5SiDxs52OCoTCBLjXYmwDZhljMnFPr+ZYlmzZg3NmjXz2fFyplCykaPNmXAqEa6mqnvcr/cCnoZ9qwnsyLW+073tJBGpCFyD61tJj0TkLhFZLiLLk5OTzyzqAHL77bdz2223UbduXadDMUHssx82ICLkjL12Ij2T+d9vYNo8Gw8v3CUmlAOwKZSMMbn58/NbHxH5VUQ+FJHaPozZBIC7776b/v37++x4lggbX4j014FFZBHg6eb9x3OvqKqKSJGHQBaRSOBd4F+qWuCM2qo6GZgM0LZt25AZallEmDJlitNhmCA3YtIi0jPyDpZ17EQGIyYtYmD3Ng5FZQJBpXKxxERHsvuAJcLGhBOHPr99Cryrqmkicjeuq81XFFD/LuAugLPOOquozRuHPPjggz49Xs5cwlu2FJgCGFMov10RVtVOqtrcw/IJsE9EagC4f3oaoWcXkPsbwVrubTkmA5tUNazvDVZVXn75ZZ588kmnQzFBaNTdnSgTE5VnW5mYKEYP6uxQRCZQiAg14svZFErGhBknPr+p6kFVzZlk9g2gwG9iVXWyqrZV1bYJCQnFO0lTonbv3s327dt9OvVn2bJleeSRR2jTxr60N8Xn1K3Rc4CB7tcDgU88lPkc6CIildyDLHRxb0NERgMVAN9+vRSERIT169czevRovv76a6fDMUFmYPc2dGt/NjHRf90cUqtqBQZcfZ6DUZlAkZhQ3m6NNsbk5pfPbznJtVtPYJ2P4zYOGj9+PA0bNiQjw7czsL7wwgt06dLFp8c04cWpRHgs0FlENgGd3OuISFsReQNAVVOAUcAy9zJSVVNEpBau23OaAj/bMPvw8ssvU79+fQYMGMDhw/ah1RTNpGG9SagUhwiUjY1m4/YDLF72m9NhFcnhw4ftOSE/SIwvb4NlGWNy89fntwfcUyr9AjwA3FqSJ2X8KykpiXPOOYfo6GifHjczM5Pff//dp8c04cWRRNh9C8yVqtrIfQtOinv7clW9I1e5Kara0L1MdW/bqaqiqk1smH2XsmXLMmPGDHbs2MEDDzzgdDgmyMTFRjP7hf40qVuVz/91G03qJnDH6I9JPpTqdGhemT17Nk2bNqVfv34+ve3KuAbM2p182H6vxhjAf5/fVHWYqjZT1XNV9XJVXe/MGRp/SEpK8umI0TlGjx5NvXr1SEtLK7ywMR44dUXY+Fj79u0ZPnw406dPZ82aNU6HY4JM03pVWTH9Ps5rXJNpT93AoaMnuHvMrIBOgHbt2sV1111H7969iY+P57XXXkNEnA4rpCTGl+dEeiaHjhx3OhRjjDFB6OjRo2zdupXmzZv7/NgNGjRAVdm2bZvPj23CgyXCIeTJJ59k6dKlfvnWzYSPFg2r89w9XVjww0YmfvSj0+F4tHz5cpo0acJnn33GuHHjWLZsGe3atXM6rJCTmOCeSzjZbo82xhhTdGvXrgXwSyKcM3K0PRplissS4RASFRV1MhlYuXJlQF/NM4Htnj4XcHWHsxk+8X+s/m2v0+GcdOLECQBatmzJzTffzOrVq3nssceIiooqpKYpjhrx7rmEbQolY4wxxdCgQQPeffddOnTo4JdjgyXCpvgsEQ5B3377LW3atGHSpElOh2KClIjw+tDeVCobw4CnZ3LsRLqj8Rw/fpzhw4fTtGlTDh8+THR0NBMnTjz5R9D4R8141xVhm0LJGGNMcVSpUoV+/fpRtWpVnx+7WrVqxMXFWSJsis0S4RDUoUMHOnfuzMMPP8zGjRudDscEqYRKcbzxxHWs/z2Zf47/zLE4Fi9eTIsWLRgzZgyXXnop2dnZjsVSHCJSWUQWisgm989KBZQb6C6zSUQGetg/R0SS/B/xX05eEbYplIwxxhTDggUL+PXXX/1ybBHhlVde4frrr/fL8U3os0Q4BEVERDBlyhRKly5N//79yczMdDokE6SuPL8hD914EW98spxPvl5bom0fP36c2267jU6dOiEiLF68mKlTp1KxYsUSjcMHhgKLVbURsNi9noeIVAaeAi4A2gFP5U6YReQ64GjJhPuX0tGRxFcoY7dGG2OMKZY777yTl156ya/H79ixo9+Ob0KbJcIhqmbNmrz++uv89NNPPPXUU06HY4LY03deyXnnJHLvuE/Yuf/PEmu3dOnS7N+/n+HDh/Prr79yxRVXlFjbPtYLmOZ+PQ241kOZq4CFqpqiqoeAhUBXABEpCwwBRpdArKdITChvV4SNMcYU2aFDh9i1a5dfBsrKkZKSwldffWXj4phisUQ4hPXt25e7776bWrVqOR2KCWLRUZFMe+oG0jKyuH3UR2Rl+e/W5IyMDEaNGsWOHTuIiIjg008/5dlnnyU2NtZvbZaAaqq6x/16L1DNQ5mawI5c6zvd2wBGAS8Bx/wW4WkkxpezUaONMcYUWc50nv6czeS9997jsssuY8+ePYUXNiYfS4RD3MSJE7nnnnuAv0bcNaaoGtauwisPduebVdt48e1v/NLGb7/9RseOHXnyySeZOXMm4LrNPxiIyCIRSfKw9MpdTl1fWXv9tbWItAIaqOosL8reJSLLRWR5cnJy0U+iAIkJ5e3WaGOMMUWWlOQa1sKfV4Rt5GhzJoLjU6YpNhEBXAMONWjQgNWrVzsckQlWt3RrxQ1XtmDUlCUsTdrO2q37aTNgPGu37j+j46oqU6ZMoVWrVmzatIkPPviAIUOG+CjqkqGqnVS1uYflE2CfiNQAcP/09AvbBdTOtV7Lva090FZEtgHfAmeLyJcFxDBZVduqatuEhASfnVuN+HLsP5RKeoaNNWCMMcZ7SUlJlCtXjtq1axdeuJhyEuEtW7b4rQ0TuiwRDhNnn302AN27d2f37t0OR2OCkYjw70euoVbV8gx4eia9HpnOum3J9H7sLVKPF396pddee42///3vtGvXjl9//ZUbbrjBh1EHhDlAzijQA4FPPJT5HOgiIpXcg2R1AT5X1YmqmqiqdYGOwEZVvawEYj4p0T2F0t6UEh+ryxhjTBAbNWoUX3311cmLMv5w1llnERERYVeETbFYIhwmateuzbx58zh06BDdu3fnyBF75s8UXYWyMUwdcT079v3JngNHUFX2pxxl0NjZRT5Wzq36AwYMYPz48SxatChUn2cfC3QWkU1AJ/c6ItJWRN4AUNUUXM8CL3MvI93bHJeY4EqEbcAsY4wxRVGpUiVat27t1zaio6M566yzLBE2xWKJcBhp1aoVM2fOZPXq1fTt25eMjAynQzJBaOP2ZKIiI8jKdj3qeiI9k/nfb2DavBVe1U9LS+Phhx+mffv2pKWlUb58eQYPHhw0zwMXlaoeVNUrVbWR+xbqFPf25ap6R65yU1S1oXuZ6uE421TVfw9aFeCvRNi+PDPGGOOdgwcP8swzz7Bp0ya/tzV58mSGDj1lZkJjChWanzxNgbp27crEiROpU6eOX29VMaFrxKRFZGTmHTn62IkMRkxaVGjdNWvW0K5dO15++WU6dOhAdrb/RqA2vpEYXw7ABswyxhjjtVWrVvH000+zfft2v7fVuXNnWrRo4fd2TOixRDgM3Xnnnbz++utERkZy/Phxp8MxQWbU3Z0oExOVZ1ts6UhGD+p82nqTJ0+mbdu27Nmzh08//ZTXXnst2KdFCgtVKpShdHSk3RptjDHGazlTJ/lzxOgcO3fu5O233+bYMUdmGTRBzBLhMPb777/TrFkz3nnnHadDMUFkYPc2dGt/NjHRkSe3xZaOou+VBX8bm56ezoQJE7jkkktYvXo1PXr0KIlQjQ+ICDWqlGP3Abs12hhjjHeSkpKoUqUKVatW9XtbP/74I7fccgsbN270e1smtFgiHMaqV6/OWWedxW233cbXX3/tdDgmiEwa1puESnGIQHzFMqQcPs4DL83FNU3uX3bt2sWRI0eIjo5m4cKFzJ8/n2rVqjkUtSmuGvHl7IqwMcYYryUlJdG8efMSeQyvfv36gM0lbIrOEuEwVrp0aWbNmkX9+vW59tprWb9+vdMhmSARFxvN7Bf606RuVT7/1+0Mv/UyZixYyYSPfjxZ5ptvvqFNmzYMHjwYgISEBEqVKuVUyOYMJCaUZ49dETbGGOMFVWXbtm0lcls0/DWXsCXCpqgiCy9iQlmlSpWYP38+F154Id26dWPp0qV2xc54pWm9qqyYfh8AjW+7jF9/28s/x39G07oJrFm6gIceeoh69erZSI4hIDG+HAu+34Cq2iB7xhhjTktE2LlzZ4mNQ1O+fHni4+PZsmVLibRnQoddETbUq1ePuXPn0rhxY6KiogqvYEw+ERERvPnEdTSsVZlrH57CA4+MoFu3bixbtoymTZs6HZ45Q4kJ5Tl2IoM/j55wOhRjjINEpLKILBSRTe6flQooN9BdZpOIDMy1/UsR2SAiq9xLVff20iLyvoj8JiI/ikjdkjkj4y8RERHExcWVWHv169e3K8KmyCwRNgCcf/75LFiwgMqVK5OWlsaJE/aB1xRN+bgYxv/jCtIzMql1+f28/e4HVKhQwemwjA8kxrvnErbbo40Jd0OBxaraCFjsXs9DRCoDTwEXAO2Ap/IlzDeraiv3st+97e/AIVVtCLwCjPPnSRj/evfdd7nrrrvIzMwssTanTZvGtGnTSqw9ExosETZ5pKenc9VVV9G8eXMWLFjgdDgmSKxcuZLs7Gw6nt+cd0bfyMFjwqBxn5wyeJYJTokJ7rmEbcAsY8JdLyAn25gGXOuhzFXAQlVNUdVDwEKgaxGO+yFwpdhzGEHrs88+Y968eURGltwTmI0bNyYxMbHE2jOhwRJhk0d0dDRPPPEEkZGRXH311Vx33XUlMhm6CU6qyssvv8z555/Pa6+9BkDvK1oxelBnPl6yhhdm2GjkoaBmgvuKsCXCxoS7aqq6x/16L+BpUJGawI5c6zvd23JMdd8WPSJXsnuyjqpmAn8CVXwauSkxOSNGl6StW7cyZswY9u7dW6LtmuBmibA5RadOnfjll1947rnn+Oyzz2jSpAmLFi1yOiwTYFJTU7n55pt5+OGH6dWrF7feeuvJfQ/2u4h+nVvy9BtfMP+7Dc4FaXyiehXXFeE9B+3WaGNCnYgsEpEkD0uv3OXUdctPUW/7uVlVWwAXu5f+xYjvLhFZLiLLk5OTi1rd+FlWVhZr164t8UR4x44dNTceVgAAIABJREFUDB8+nF9//bVE2zXBzRJh41Hp0qUZNmwY69ato1+/fpx//vkAHD161OHITCDYuXMnF110Ee+99x7PPfccH374IeXKlTu5X0SY8M9etGpUnVtHfsiG35NZu3U/bQaMZ+3W/ac5sglEsaWjqFw+1q4IGxMGVLWTqjb3sHwC7BORGgDun5469F1A7VzrtdzbUNWcn0eAd3A9Q5ynjohEAhWAgwXEN1lV26pq24SEhDM9XeNjW7du5cSJEzRr1qxE27W5hE1xWCJsTqtOnTq8+eabVKhQgYyMDC688EJuuukm9uzZU3hlE7I2b97Mrl27mDdvHsOGDfM4pU5s6Sjef+4mYktH0uefb9Pz4ems25ZM78feIvV4ugNRmzORGF/eBssyxswBckaBHgh84qHM50AXEankHiSrC/C5iESKSDyAiEQBPYAkD8e9HvhCbZCJoHTw4EEaNmxIixYtSrTdxMRESpcubYmwKRJLhI3XsrOzuf766/n4448555xzePXVV0t0REDjvN27dwNw6aWXsnXrVrp163ba8rWrVeDtUX9jy64U9hw4gqqyP+Uog8bOLolwjQ8lJpSzK8LGmLFAZxHZBHRyryMibUXkDQBVTQFGAcvcy0j3ttK4EuJfgVW4rgL/x33cN4EqIvIbMAQPo1Gb4HDBBRewadOmk3cSlpSIiAibQskUmSXCxmulS5fm6aefJikpiYsuuoiHHnqINm3asHPnTqdDMyXgww8/pEGDBsye7Upiy5Yt61W9zTsPEhVVimz3l/sn0jOZ//0Gps1b4bdYje8lxpe3RNiYMKeqB1X1SlVt5L6FOsW9fbmq3pGr3BRVbeheprq3papqG1VtqarNVPUfqprl3ndCVW9wl2+nqlucOUMTzBo0aMC2bducDsMEEUcSYR9MyP6ZiPwiImtE5HURKVVy0ZuGDRsyf/58PvroI2rXrk316tUBSEtLczgy4y/jx4+nb9++nHfeeVx88cVFqjti0iLSM7LybDt2IoMRk2wAtmCSmFCefYdSycjMKrywMcaYsNS1a1dGjhzpSNvTp09n2bJljrRtgpNTV4TPdEL2vqp6LtAcSABuKJGozUkiwnXXXcfcuXOJjIzkzz//pGHDhjz66KOkpKQ4HZ7xEVVl2LBh3H///fTs2ZNFixZRpUrRZrQYdXcnysRE5dlWOqoUowd19mWoxs8SE8qjquw9aAPmGWOMOVV6ejqLFy/m+PHjjrRfqVKlEp272AQ/pxLhM5qQXVVz7s+LBKIp+vD9xscyMjLo3LkzL7/8MvXr12fcuHGOdYTGdxYuXMjYsWMZNGgQH330EbGxsUU+xsDubejW/mxiol1/nEQgPTOLujU83ghiAlSNnCmUDtjt0cYYY061adMmMjMzS3zqpBxbt25l8ODBrF271pH2TfBxKhE+4wnZReRzXMP2HwE+LKghm2+uZMTHxzNlyhR++eUXLrnkEoYOHUqjRo2w33lwyhmss0uXLixevJgJEyZQqlTxn0CYNKw3CZXiEIGaCeU556wEbhj2Dr9sstHHg0ViQnkAGznaGGOMR0lJrkHAnUqE09LSmDBhAj///LMj7Zvg47dE2M8TsqOqVwE1cI1CeMVpytl8cyWoefPmzJkzh6+//pqbb76ZnN/5mjVrTiZXJrDt3buXyy67jBUrXINZXXHFFR6nRyqKuNhoZr/QnyZ1q/LJiwP49OUBlI+LodcjM9i6226lDwaJCa4rwjZgljHGGE+SkpIoVaoU55xzjiPt161bFxGxkaON1/yWCPtzQvZcbZzANYddnuTaOO/iiy9m3LhxAGzbto3WrVtz8cUX89133zkcmTmdjRs30r59e5YvX86BAwd8euym9aqyYvp9NK1XlVpVK/DpSwPIyMyix5Dp7Eux504DXXyFMkRFlmK33RptjDHGgzp16nDjjTcSExPjSPsxMTHUrFnTEmHjNadujT6TCdnL5kqiI4HuwPoSiNkUU61atRg/fjxbtmyhY8eOXHvttfb8RgBaunQpHTp0IDU1lS+//JKrrrrKr+01rpvArOdvYe/BI1z7yAwOp57wa3vmzERERFAjvhy7k+3WaGOMMae64447mDFjhqMxNGjQwBJh4zWnEuEzmZA9DpiTa0L2/cDrJX8KxluRkZHcddddbNq0iWeffZYlS5Zw3nnn2ejSAWTZsmVcccUVVKxYke+//57zzz+/RNpt16w274z6G0lb9vG34e+Slp5ZIu2a4kmML2dXhI0xxpwiKyuLzEzn/4Y3atTIpvM0XnMkET7DCdn3qer57gnZm6vq/arq/P88U6i4uDiGDx/O5s2b+e9//0vlypUBGDZsGAsWLLBniB3UqlUrHnzwQb7//nsaNmxYom1fdeHZTBp6LV/+vJXbR39EVlZ2ibZvvFcjvrw9I2yMMeYUv/zyC3FxcXz++eeOxjF58mSWL18OwNqt+2kzYDxrt3p6AtOzkqhjbfi/jrecuiJswlh8fDz9+vUD4I8//uCtt97i6quvpmXLlkybNo309HSHIwwfixYtYv/+/URFRfHcc89RtWpVR+K4qWsrxgy+io+XrOHhf823L0UCVGJCOfbYqNHGGGPySUpKIj09nTp16nhV3l8JUc7gnqnH07n20Rms25ZM78feIvV44Z8tS6KOteH/OkVhibBxVMWKFdm8eTPTpk1DRLj11lupV68eP/74o9Ohhby5c+fSvXt3HnnkEadDAeDBfhfx0I0XMenjnxg77SvAv98CmqJLjC/P0ePp9jy3McaYPJKSkoiOjvbqrjJ/JkQ7d+6kZ8+e9Lp7NMmHUlFV9qccZdDY2YW2cfeYWX6vY234v05RSDhdeWnbtq3m3C5hAo+q8r///Y/x48czY8YMKlasyIoVK6hatSq1a9cu/ADGa3PnzqVPnz60bNmShQsXUrFiRadDAiA7O5u7xszm7c9W8eID3fi/979n5/7D1K5WgZ+n30dcbLTP2xSRFara1ucHdpC/+rr3Fv7KbSM/ZOWM+2lc16ajMyaYWF9n/Onqq69m165d/PLLL4WWveXJ95n33QZOpGcSEx1Jj46NmfFM3yLV6d7xHCYP682RY+mkHk/n6PF0jh5LY8/e/Vx/VTtKN+iC1LjgZP2oyAi6d2xMm8Y1PR5/xfpdzPt2AxmZWX6rY22ceZ0yMVG8/ODVDOzexmOdHN72d5YIm4B2/vnns2rVKm688UYeeughWrdu7XRIQS9Qk+AcGZlZ/G34uyz4YSNRkaXIyMzy+g9lcdiHQ+99s3IrXR6YyrxXBnJF2wY+P74xxn/Cva9b/N1Kel9/I80bVMvzpeqDDz5I9+7d2bhxI4MHDz65PfV4Okmb9/Hqi89y+819WLVqFY8++ugpxx01ahQXXnghP/zwAw8/OpSkzfvytPHiiy9y7rnn8sUXXzBmzJhT6j80dCQjpq3g75dUZNb7007ZP3XqVGrVqsXMmTOZPHnyybhy2njvvfeoUqUK06ZN46233jql/pw5c9i69whdb3yAWpF7TvlCeeHChQC88sorzJ8/P8+5t2lWhyWLXM/8PvfccyxZsiRP3cqVK/P+++8DUKFCRUpVbsS3X8yjab2CH7N69d3veOo/i0jP+Cu5EaBy+VjKxXmedulI6glSjhzH25TlxA8vUKpqc6IadPOuggkqCZXi2D7nn6ct421/F+mzqIzxg5kzZ/LKK6/w5ptvMmPGDFq3bs2IESPo3bu306EFpaysLJ544omATYIBoiJL0a3D2Xy+dNPJbwFPpGcy//sNTJu3otBvAY3/JCaUB7ABs4wxQSX1eDp3jP6Q1GOprFy/nVaNahAR4XqWNGek46ysLI4dOwZAdrayatMe0jMyeXLS//jbddfk2Z9bVpbr79TRYydYuX476RmZedrIzs4+2U7++tnZyr3jZrP/RAzP/raCszT1ZFx/lXHVz8jI4OjR1JNx5bSRc0ErIyPDY3xHj6Vx7aMz2Jv8J8l/5j333NLT0zl27Fiec/9p9VZSj6cTFxtNWlraKccvU6bMyd9vqeqtSY1OpPdjb+W5g+t4WgbfrtrGwp9+Y9FPv7FuW/IpbStw9EQGV7U/+5R9AB8tWeMxCY6LjebZe7pQNjaauNhoysZGU7ZMNH16vM/+o3/kjTUmiufv68pNXVt5bOPtz1bxz/GfcexEht/qWBu+qTN6UGeP5YvDrgiboJCSksI777zD1KlTGTx4MLfffjsHDx7k+++/p2vXrkRFRTkdYtDYu3cvMTExAZkE5zir5ziSD6West2bbwGLKtyvkhTFsRPpVOk8mmfuvJLHBlzq8+MbY/wnnPu6ot6K64tbd/1RJ1jauLhVXa44vwGLfvqNb3/5nbT0TEpHR9Lx3DpUKBvDfHfZHGVionjloe4MuPo8j8efNm8FQ16df0pCVFCdvn378tnibyjV6h57D4O4jeLWAbs12iNLhENDdnY2ERERTJgwgcGDB1OtWjVuueUWBg4cSIsWLZwOLyDNnz+f9957jylTphAZGfg3gnj6o1c6qhT/euSaAv9QFlc4fzgsjhrdnuNvnVvy6pAefjm+McY/wrWv8/T3JDqyFH07t6BDi1NHOP5+9e98sGh1nlt3T1e+pOoEUxs5GtdJoFO7hnS+oCEdz61DmRjXVWJ/J0Tjxo1j3vwF7Invwa7kw9SuVtGrcUZSj6fTuv+/2bn/T7/VsTb8XwcsEfbIEuHQkpGRwWeffcbUqVP59NNPyczMpG3btnz99dfExsY6HV7AmD9/Pr1796ZFixYsXryYChUqOB2SV3L/0XPPhsBHY2+mW4dzfNqOPz8cikhl4H2gLrAN6KuqhzyUGwg84V4drarT3NujgfHAZUA28LiqflRYu/7s65r1e5XkQ0f58vW7TvscWG5rt+6n/1MfMOOZvl7VKWr5kqxjTLAK10S4oDuMjH9VqVCGnXOHetxXUglRoP4tsTb8X8cSYQ8sEQ5dycnJvPvuu6xbt46JEycC8Pzzz9OoUSO6du0atolx7iR44cKFVKpUyemQvJb7j17NhPJUqVCGtVuTeXtkX665uInP2vFzIvw8kKKqY0VkKFBJVf+Zr0xlYDnQFtejUiuANqp6SESeAUqp6hMiEgFUVtUDhbXrr74u9Xg6tXqM5UR6JmdVL+qHF+9G/y5q+ZKsY0wwC9dE2NMV4djSkTx1x5X0uaL5KeU/XJzEyDe/4Hiad+VLqk4wtVHYrc5gX14a/7JE2ANLhMPHiRMnaNCgAbt376Zs2bJcc8019O3bl65duxIT43lUwlCzYMECrr322qBMgnPk/qOXGF+Oax6ezqqNe5j+9A30vqyZT9rwcyK8AbhMVfeISA3gS1U9J1+ZG91l7navT3KXe1dEdgCNVbVIlzP81dfd8uT7zP5qLVnZGrbPHRkTrMI1EYbQ6YdCpQ1/O3jwIJ06dWLIkCH079/fsTiMcywR9sAS4fCSmZnJl19+yQcffMDHH3/MwYMHGTNmDEOHDiU9PZ3s7OyQToq/++47nnjiCT7++OOgTII9+fPoCXo9OoPl63YxdUQfbrjyzJ8J93Mi/IeqVnS/FuBQznquMo8AMao62r0+AjgOvAGsBmbiujV6M3Cfqu4rrF1/9HWerqqUihCa1a/GWdU9D7y2fe8frNmyj6xs9apOUcv7so63cxMaE6x80dedyeMeIlIO+CZXsVrAW6r6oIjcCrwA7HLvG6+qbxQWj7d9Xag8/xgqbfhbVlYWsbGxDBkyhLFjxzoWh3GOt/1dREkEY4wTIiMj6dSpE5MnT2bPnj3873//45ZbbgFc8+pVrVqV/v37M2fOHNLS0hyO1nfWr18PwEUXXcQXX3wRMkkwQIWyMXz60gAubF6bW0d+yLv/+8XpkBCRRSKS5GHplbucur51LMo3j5G4Pih+r6rnAT8AL54mjrtEZLmILE9OPnV6ijM1YtKiPEkwQFa2smbrfrbv/cPjsmbr/jzJZmF1ilrel3WOnchgxKRFPvyNGROShgKLVbURsNi9noc7WX4KuABoBzwlIpVU9YiqtspZgN+Bj3NVfT/X/kKT4KKIi41m9gv9aVK3KrOev6XQJK2o5UuqTqi04W+lSpWiXr16bN682dE4TOCzK8ImLK1cuZLXXnuNjz/+mEOHDlG+fHl69erFxIkTiYuLczq8Yps4cSL33XcfH3zwAX369HE6HL9JPZ5On6Fv8/XKbUwedi23dGtd7GMF6q3RwHvAUaCcqmaLSG3gM1Ut9J7wkroi7OspL0qijYLqxERH8n8P9/D5qOTGBAofXRE+o8c9cpU5G1cifZaqqvuKcFtVva8o8djnOlOQbt26sW/fPn7++WenQzEOsCvCxpxG69ateeONN9i3bx8LFizg+uuvZ8OGDScnh3/11VeZNGkSO3fudDhS76gqTz31FPfeey9XX3013bp1czokv4qLjebjcTdzRdv63DVmNv+duwJwPVPcZsB41m7d73CEJ80BBrpfDwQ+8VDmc6CLiFQSkUpAF+Bz9xXkT3HdFg1wJbDWv+EWbGD3NnRrfzYx0a7pt2KiI7m6wzmnTRyLWqck2vBUR4CMrCzqJ1Yu+BdgjAGopqp73K/3AtU8lKkJ7Mi1vtO9Lbd+uK4A574a00dEfhWRD91f/BlTbA0aNGDz5s2E0wU/U3SWCJuwFhUVRdeuXXnzzTdZunQp4p6n55133mHQoEHUrl2b8847jyeffJJffnH+NlxPMjMzGTRoECNHjuT2229n1qxZJxP6UFYmJpoPx9xE53YNuGfcJ4yf+QPXPjqDdduS6f3YW6QeT3c6RICxQGcR2QR0cq8jIm1F5A0AVU0BRgHL3MtI9zaAfwJPi8ivQH/g4RKOP49Jw3qTUCkOEahauSyvD73W53VKoo38dRITylM/sTI9H5nBop9+K7SuMaHMj4975NYPeDfX+qdAXVVtCSwEpp0mPr8+BmJCQ8eOHenatSvp6QHxWcAEKlUNm6VNmzZqjDeys7N1zZo1Om7cOO3YsaNGREToP/7xD1VVzcjI0NmzZ+uRI0ccjtJl4cKFCujw4cM1Ozvb6XBK3Im0DL3usRka03GElr30KY3pOEIrXvGM3vLk+17VB5ZrAPRPvlz82det2bJPz+v/b12zZZ/f6pREG/nr7Es5ou1ufU3LX/60zv12ndfHMCZY+KKvAzYANdyvawAbPJS5EZiUa30ScGOu9XOBjadpoxTwpzfx2Oc6Y4wn3vZ39oywMV44ePAg6enp1KhRg2+//ZaLL76Y6OhoLr/8cq666iouv/xyWrZsSUREyd1koaonr2CvWLGCNm3Cd7TbKXOW88BLnxZrFOBwnlLE5JVy+Bg9H5nBLxv3MPXJ67m+gHkzjQlGPnpG+AXgoP41N3plVX0sX5nKuOZDz3k+4Wdcc6OnuPePBdJU9alcdWqo+5ZrEekN/FNVLywsHuvrTGGys7NL9LOZCQz2jLAxPlSlShVq1KgBwAUXXMAXX3zBfffdx9atWxkyZAitW7fmm29cs0Js27aNpKQk/Pkl086dO2nXrt3JNsM5CQZ4+o3FNgqwOWOVy5dh/isDadesFgOfmclbC1Y6HZIxgeZMH/cA6Eve26IBHhCRNSLyC/AAcKtfz8KEvPT0dBISEnjuueecDsUEMEuEjSmiqKgoLr/8cl566SU2bNjAjh07mD59OhdccAEAr7/+Oi1atKBatWr07duXiRMnsn79ep8lxuvWraNDhw5s2LCBjIyMwiuEgVF3d6JMTFSebWViohg9qLNDEZlgVT4uhjkvDuCy8+px53Oz+M/sZU6HZEzAUNWDqnqlqjZS1U45Ca6qLlfVO3KVm6KqDd3L1HzHqK+q6/NtG6aqzVT1XFW9PP9+Y4oqOjqaqKgotmzZ4nQoJoBZImzMGapVqxb9+/cnJiYGgPvuu4+pU6fStWtXfvjhB+69917at29PdnY2AN9++y2rV68mKyuryG398MMPdOzYkfT0dL766iuuuOIKn55LsCrOyMHGFCQuNpqPxt7M1R3O5oGXPuVf738PFG9U8qLWCZU2AjWucG7DmHBTv359m0vYnJY9I2yMH6kqmzdvZvPmzVx11VUANG3alHXr1lGuXDnatWtH+/bt6dSpE5deeulpj7V69WouuOACatasyeeff079+vVL4hSCRurxdFr3/zc79/9J7WoV+Xn6fcTFRhdaz54RNgVJz8jktlEf8fGSNQy/9TJmLFjJzv2HqV2tglf/vv76N+ldnaKWD9Q2AjWucG4DrK8z4WfgwIEsXrw4aKbCNL5jzwgbEwBEhIYNG55MggE++eQTpk+fTv/+/UlJSWHMmDG8+eabgCtxvueee5gwYQI///xznlufmzVrxqOPPsp3331nSbAHcbHRzH6hP03qVmXW87d49cHQmNOJjopk2pPXc9NV5/Lcf79k94HDqCr7U44yaOzsQuvfPWYWyYdSva5T1PKB2kagxhXObRgTjho0aMCuXbs4fvy406GYAGVXhI1xWGpqKkeOHKF69eokJyfTokUL9u3bB0BsbCyNGzdm1qxZ1KlTx+FIQ5NdJTGF+e/c5dz/4lwys7JPbisVIbRtUpMGtap4rLN550GWr99Nlpd1ilo+UNsI1LhCuQ0bId/6OuPZ119/zZw5c3jiiSeYPn06+f+txMfH8/LLLwPw/PPPk5SUlGd/rVq1Tg62NXLkSH77Le888w0bNuTJJ58EYPjw4adceW7evDmPPeYaVH3IkCEcOHAgz/62bdvywAMPAHDvvfdy9OjRPPs7duzIXXfdBcBtt912yiN1nTp1YsCAAWRmZnL77befcv7XXHMNN9xwA0eOHGHw4MGn7L/++uvp2bMnycnJPPzww6fs79+/P507d2bHjh08/vjjp+y/8847ufjii9m0aROjRo06Zf/999/P+eefz+rVq3nhhRdO2f/oo4/SokULli1bxr///e9T9o8YMYJGjRqdst0b3vZ3kcU6ujHGZ+Li4oiLiwMgISGBPXv2sH37dpYuXcrSpUv58ccfmTlzJo888ojDkRoTnp6cvDhPEgyQla0sW7eLfSmpHuts3/cH2flGMj9dnaKWD9Q2AjWuUG4jZ4T8whJhY8LNJZdcwiWXXAK4Bhr99ttv8+yvWbPmyddJSUmn7D/nnHNOvv7111/5+eef8+w/duzYydcrV65kw4YNBcaybNkydu3alWdbzmc/gKVLl/LHH3/k2V+tWrWTr7///vtTBkjNuTswOzv7lNgBWrRoAUBmZqbH/Rde6JqhLC0tzeP+K6+8EnBdsPG0v2fPngAcPnzY4/6bbroJgD/++MPj/jvvvBNwTVHqaf/hw4dP2eZrdkXYGBPW7CqJKcy0eSsY8up8jp3460NImZgoXnmoe4EDshW1Tqi0EahxhXMbOayvM8aEC3tG2BhjjPGB4oxKXtQ6odJGoMYVzm0YY4zxzBJhY4wxphCThvUmoVIcIlC1clleH3qtz+uEShuBGlc4t2GMMeZUlggbY4wxhSjOqORFrRMqbQRqXOHchjHGmFM58oywiFQG3gfqAtuAvqp6yEO5gcAT7tXRqjot3/45QH1Vbe5Nu/YsiTEmP3tuzhgTDqyvM8aEi0B/RngosFhVGwGL3et5uJPlp4ALgHbAUyJSKdf+64Cj+esZY4wxxhhjjDGn41Qi3AvIubo7DfD0gMtVwEJVTXFfLV4IdAUQkbLAEGB0CcRqjDHGGGOMMSaEOJUIV1PVPe7Xe4FqHsrUBHbkWt/p3gYwCngJOJa/Un4icpeILBeR5cnJyWcQsjHGGGOMMcaYUBDprwOLyCKguoddj+deUVUVEa8fVBaRVkADVX1IROoWVl5VJwOTwfUsibftGGOMMcYYY4wJTU4NlrUBuExV94hIDeBLVT0nX5kb3WXudq9PAr4EKgIjgHRciXxV4HtVvcyLdpOB34sQajxwoAjlA5WdR2AJlfOA0DiXOqqa4HQQvmR9XdALlfOA0DmXUDgP6+tC430EO49AFCrnEirn4VV/51Qi/AJwUFXHishQoLKqPpavTGVgBZAzS/zPQBtVTclVpi4w19tRo4sR5/JQGGHRziOwhMp5QGidSzgLlffRziPwhMq5hMp5hLtQeR/tPAJPqJxLqJyHt5x6Rngs0FlENgGd3OuISFsReQPAnfCOApa5l5G5k2BjjDHGGGOMMaY4/PaM8Omo6kHgSg/blwN35FqfAkw5zXG2AX65GmyMMcYYY4wxJjQ5dUU4WEx2OgAfsfMILKFyHhBa5xLOQuV9tPMIPKFyLqFyHuEuVN5HO4/AEyrnEirn4RVHnhE2xhhjjDHGGGOcYleEjTHGGGOMMcaElbBPhEWkq4hsEJHf3CNY599fWkTed+//0Zu5i53ixbkMEZG1IvKriCwWkTpOxFmYws4jV7k+IqIiEpCj23lzHiLS1/2erBGRd0o6Rm948e/qLBFZIiIr3f+2rnYiTnN61tcFHuvrAo/1d6EhVPo76+sCi/V1IUhVw3YBSgGbgfpANPAL0DRfmXuB192v+wHvOx33GZzL5UAZ9+t7AvFcvDkPd7lywNfAUqCt03EX8/1oBKwEKrnXqzoddzHPYzJwj/t1U2Cb03HbUqz30fq6ADsPdznr6wLrXKy/C/AlVPo76+ucj70Y74f1dUG2hPsV4XbAb6q6RVXTgfeAXvnK9AKmuV9/CFwpIlKCMXqr0HNR1SWqesy9uhSoVcIxesOb9wRcU2uNA06UZHBF4M153Am8pqqHAFR1fwnH6A1vzkOB8u7XFYDdJRif8Y71dYHH+rrAY/1daAiV/s76usBifV0ICvdEuCawI9f6Tvc2j2VUNRP4E6hSItEVjTfnktvfgQV+jah4Cj0PETkPqK2q80oysCLy5v04GzhbRL4TkaUi0rXEovOeN+fxNHCLiOwE5gP3l0xopgisrws81tcFHuvvQkOo9HcsZ4igAAAFZklEQVTW1wUW6+tCkCPzCBtnicgtQFvgUqdjKSoRiQBeBm51OBRfiMR1G81luL7F/VpEWqjqH45GVXQ3Av9V1ZdEpD0wQ0Saq2q204GZ8GZ9XcAIlb4OrL8zAcj6uoBhfV2QCfcrwruA2rnWa7m3eSwjIpG4bg84WCLRFY0354KIdAIeB3qqaloJxVYUhZ1HOaA58KWIbAMuBOYE4MAK3rwfO4E5qpqhqluBjbg60EDizXn8HfgAQFV/AGKA+BKJznjL+rrAY31d4LH+LjSESn9nfV1gsb4uFPs6px9SdnLB9c3NFqAefz0s3ixfmcHkHVDhA6fjPoNzaY3r4fhGTsd7JueRr/yXBOagCt68H12Bae7X8bhuU6nidOzFOI8FwK3u101wPUciTsduS5HfR+vrAuw88pW3vi4wzsX6uwBfQqW/s77O+diL8X5YXxdki+MBOL0AV+P6xmYz8Lh720hc36yB6xuQmcBvwE9AfadjPoNzWQTsA1a5lzlOx1yc88hXNiA7TC/fD8F1O9BaYDXQz+mYi3keTYHv3B3pKqCL0zHbUqz30fq6ADuPfGWtrwuMc7H+LgiWUOnvrK8LrMX6Oufj9vUi7pM1xhhjjDHGGGPCQrg/I2yMMcYYY4wxJsxYImyMMcYYY4wxJqxYImyMMcYYY4wxJqxYImyMMcYYY4wxJqxYImyMMcYYY4wxJqxYImyMMcYYY4wxJqxYImyMMcYYY4wxJqxYImxMEYlIUxG5VURqi0g5p+Mxxhh/sL7OGBMOrK8LX5YIm5AiIleKyFvFrBsrIl+JSCn3egsR+V1E7slXNAq4H+gNHBWRaBH5WkQizyx6Y4zxjvV1xphwYH2d8SdLhE2oORdYVcy6twMfq2oWgKquBvoBA/KVqw1MBX4DyqlqOrAY+Fsx2zXGmKKyvs4YEw6srzN+Y4mwCTXnAr8AiEhjEflCRFaJyCIRiXdvb+L+pu9XEXlURH5z170Z+CTf8fYDzXJvUNW5wIeqOl9VD7s3z3bXN8aYkmB9nTEmHFhfZ/zGEmETas4FVolIaeAjYIiqtgIWAg+5b3N5G/iHqrYE6gNJIhIN1FfVbfmONxYoLSJ1cm9U1b35yiUB5/v8bIwxxjPr64wx4cD6OuM3du+7CSoisgio7mHX48B8oIKqJovI34BvVTXndpq1QE/gOuAXVV2Za/t+IB74I19b3YA4YB6ubw9/LyguVc0SkXQRKaeqR4p9gsYYg/V1xpjwYH2dcZIlwiaoqGqngvaJSEtgnXu1KbA61+4WuDrHluR91qQ58BlwHIjJdawYYByuTvY2d7n5hYRXGjjhzXkYY8zpWF9njAkH1tcZJ9mt0SaUnHyOBNiFq9NEROoD/YHpwEHgbPf2VsAtuL5JPASUcneUAE8A09231KzG1WEWSESqAAdUNcOXJ2SMMR5YX2eMCQfW1xm/sivCJpScCyxzv54BXC0iq3F9K3i7qh4UkRnAPPf2L4FtqrrFXed/QEcR2QF0Bi5yb18NDC+k7ctx3WpjjDH+Zn2dMSYcWF9n/EpU1ekYjCkxIlJWVY+6Xz+K69mTJ9zr5wEPqWr/Yhz3Y2Coqm70acDGGFMM1tcZY8KB9XXmTNit0SbcPCQia0RkFVAXGJWzQ1V/BpbkTLzuLffIhLOtszTGBBDr64wx4cD6OlNsdkXYGGOMMcYYY0xYsSvCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCiiXCxhhjjDHGGGPCyv8DC79bI31tsBkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.rcParams['figure.figsize'] = (16, 4.)\n",
    "covar_1 = 4 # advection term\n",
    "covar_2 = 7 # diffusion term\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_lasso[i,:],'-d',color=c[19])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_lasso[i,:],'-s',color=c[11])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, coefs_lasso[i,:],'--',color='black')\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $\")\n",
    "plt.ylabel(\"  coefficients \")\n",
    "        \n",
    "plt.subplot(1,3,2)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_ihtd[i,:],'-d',color=c[19])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_ihtd[i,:],'-s',color=c[11])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, coefs_ihtd[i,:],'--',color='black')\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $\")\n",
    "        \n",
    "plt.subplot(1,3,3)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_STR[i,:],'-d',color=c[19])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, coefs_STR[i,:],'-s',color=c[11])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, coefs_STR[i,:],'--',color='black')\n",
    "        \n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stability selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/lustre/projects/jupyterhub/maddukon/Data_driven/Stability_selection_results /PDE-STRIDE/Iterative_solvers_l0_l1.py:453: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  y_reduced[i]   = y[int(ordered_reduced[i])]  # changed\n",
      "/lustre/projects/jupyterhub/maddukon/Data_driven/Stability_selection_results /PDE-STRIDE/Iterative_solvers_l0_l1.py:419: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if lam != 0: coeff[0,biginds] = np.linalg.lstsq(X[:, biginds].T.dot(X[:, biginds]) + lam*np.eye(len(biginds)),X[:, biginds].T.dot(y))[0].reshape(len(biginds),)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " done with sample  0  size  1500\n",
      " done with sample  5  size  1500\n",
      " done with sample  10  size  1500\n",
      " done with sample  15  size  1500\n",
      " done with sample  20  size  1500\n",
      " done with sample  25  size  1500\n",
      " done with sample  30  size  1500\n",
      " done with sample  35  size  1500\n",
      " done with sample  40  size  1500\n",
      " done with sample  45  size  1500\n",
      " done with sample  50  size  1500\n",
      " done with sample  55  size  1500\n",
      " done with sample  60  size  1500\n",
      " done with sample  65  size  1500\n"
     ]
    }
   ],
   "source": [
    "B = 100; reduced_size = 1500; rescale = 1;\n",
    "subsample_info_lasso, subsample_info_iht_d, \\\n",
    "            subsample_info_STR, alphas_sample = stability_selection(X_no_intercept, u_t, reduced_size, M, B, rescale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covar_1 =  4#PP + 1\n",
    "covar_2 =  7#2*PP + 1\n",
    "\n",
    "normalized_alphas_sample = np.zeros((B, M))\n",
    "for sample in range(0, B):\n",
    "    normalized_alphas_sample[sample,:] = alphas_sample[sample,:]/alphas_sample[sample,0]\n",
    "neg_log_stable_lambda = -np.log10(normalized_alphas_sample[0,:])\n",
    "\n",
    "importance_plot = np.zeros((p, M))\n",
    "for alpha in range(0, M):\n",
    "    for i in range(0, p):\n",
    "        for sample in range(0, B):\n",
    "            if(subsample_info_lasso[sample, i, alpha] != 0):\n",
    "                importance_plot[i,alpha] = importance_plot[i,alpha] + 1\n",
    "\n",
    "        importance_plot[i,alpha] = importance_plot[i,alpha]/B\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1 ):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-d',color=c[19],label='$cov = %s$' %description_trunc[i])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-s', color=c[11],label='$cov = %s$' %description_trunc[i])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '--', color='black')\n",
    "\n",
    "\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.title(\" lasso \")\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $ \")\n",
    "plt.ylabel(\" $ \\Pi $\")\n",
    "plt.ylim(0, 1.05)  \n",
    "plt.legend(loc='best', prop={'size':13})\n",
    "plt.text(0.3, 0.85,'$\\pi_{th}=0.8$',color='red')\n",
    "plt.axhline(y=0.8, color='r', linestyle='-')\n",
    "pylab.rcParams['figure.figsize'] = (15, 7.5)\n",
    "\n",
    "importance_plot = np.zeros((p, M))\n",
    "for alpha in range(0, M):\n",
    "    for i in range(0, p):\n",
    "        for sample in range(0, B):\n",
    "            if(subsample_info_iht_d[sample, i, alpha] != 0):\n",
    "                importance_plot[i,alpha] = importance_plot[i,alpha] + 1\n",
    "\n",
    "        importance_plot[i,alpha] = importance_plot[i,alpha]/B\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1 ):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-d',color=c[19],label='$cov = %s$' %description_trunc[i])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-s', color=c[11],label='$cov = %s$' %description_trunc[i])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '--', color='black')\n",
    "        \n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.title(\" Iterative hardthresholding with debias (IHT-d) \")\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $ \")\n",
    "plt.ylabel(\" $ \\Pi $\")\n",
    "plt.ylim(0, 1.05)  \n",
    "plt.legend(loc='best', prop={'size':13})\n",
    "plt.text(0.8, 0.85,'$\\pi_{th}=0.8$',color='red')\n",
    "plt.axhline(y=0.8, color='r', linestyle='-')\n",
    "pylab.rcParams['figure.figsize'] = (15, 7.5)\n",
    "\n",
    "importance_plot = np.zeros((p, M))\n",
    "for alpha in range(0, M):\n",
    "    for i in range(0, p):\n",
    "        for sample in range(0, B):\n",
    "            if(subsample_info_STR[sample, i, alpha] != 0):\n",
    "                importance_plot[i,alpha] = importance_plot[i,alpha] + 1\n",
    "\n",
    "        importance_plot[i,alpha] = importance_plot[i,alpha]/B\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "for i in range(0, p):\n",
    "    if(i == covar_1 ):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-d',color=c[19],label='$%s$' %description_trunc[i])\n",
    "    elif(i == covar_2):\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '-s', color=c[11],label='$%s$' %description_trunc[i])\n",
    "    else:\n",
    "        plt.plot(neg_log_stable_lambda, importance_plot[i,:], '--', color='black')\n",
    "        \n",
    "\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "plt.title(\" STRidge \")\n",
    "plt.xlabel(\" $ -log(\\lambda^{*}) $ \")\n",
    "plt.ylabel(\" $ \\pi $\")\n",
    "plt.ylim(0, 1.05)  \n",
    "plt.text(0.2, 0.85,'$ \\sigma =0.8$', color='red')\n",
    "plt.legend(loc='best', prop={'size':13})\n",
    "plt.axhline(y=0.8, color='r', linestyle='-')\n",
    "pylab.rcParams['figure.figsize'] = (17, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}