low rank AR(1) Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "#Use this line to set the number of  CPU used for execution of the program (if it's not working, use the commented line:\n",
    "#see https://stackoverflow.com/questions/17053671/python-how-do-you-stop-numpy-from-multithreading/21673595 for more information)\n",
    "\n",
    "import os\n",
    "\n",
    "os.environ['OPENBLAS_NUM_THREADS'] = '5' \n",
    "os.environ['MKL_NUM_THREADS'] = '5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true
   },
   "source": [
    "# Low rank VAR(1) Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main goal of this notebook is to implement and to test the low_rank VAR(1) algorithm. In order to do so, we will:\n",
    "- generate data according to a low-rank high dimensionnal VAR(1) process\n",
    "- implement the estimation of the low rank VAR(1) models\n",
    "- compare this model with a standard VAR(1) among differents  features\n",
    "\n",
    "\n",
    "We want to simulate an VAR(1) processus generated by a matrix A of rank r. We consider the time series $(X_t) \\in \\mathbb{R}^M$ of length generated by \n",
    "\n",
    "$$ X_{t+1} = A X_t + \\epsilon_t   $$\n",
    "    \n",
    "where $\\epsilon_t \\sim \\mathcal{N}(0,\\sigma^{2})\\$ \n",
    "\n",
    "\n",
    "## I Generation\n",
    "\n",
    "\n",
    "In order to generate A, we use the following procedure:\n",
    "\n",
    "- We uniformly generate U and V matrices of respective shape $M\\times r$ and $r \\times M$. \n",
    "- We transform U and V using Gram-Schmdit otrhonomalization\n",
    "- We generate $\\lambda_1,\\dots\\lambda_r$ r singular values according to a $\\mathcal{\\beta}(a,1)$ distribution\n",
    "- We take $A = U diag(\\lambda_1, \\dots, \\lambda_r) V^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_unit_matrix(M,P):\n",
    "    \"\"\"Generation of a MxP unit matrices using Gram Schmidt Orthonormalization\"\"\"\n",
    "    list_col = [ np.random.uniform(low = - 1, size = M) for i in range(P)]\n",
    "    \n",
    "    new_list_col = list()\n",
    "    for  col in list_col:\n",
    "        new_col = col\n",
    "        for old_col in new_list_col:\n",
    "            new_col -= old_col.dot(col) * old_col\n",
    "        new_col = new_col/np.linalg.norm(new_col)\n",
    "        new_list_col.append(new_col)\n",
    "    \n",
    "    list_col = [col.reshape(M,1) for col in new_list_col]\n",
    "    \n",
    "    return np.concatenate(list_col , axis = 1)\n",
    "\n",
    "\n",
    "def generate_matrix(M, r, a) : \n",
    "    U = random_unit_matrix(M,r)\n",
    "    V = random_unit_matrix(M,r).transpose()\n",
    "    diag = np.diag(np.random.beta(a,1, size =r))\n",
    "    A =  U.dot(diag.dot(V))\n",
    "    return A\n",
    "\n",
    "def tronc_gaussian(sigma, limit = 1):\n",
    "    \"\"\"Truncated centred gaussian distribution with deviation sigma\"\"\"\n",
    "    ret = np.random.normal(scale = sigma)\n",
    "    if abs(ret) > limit:\n",
    "        return tronc_gaussian(sigma, limit=limit)\n",
    "    return ret\n",
    "\n",
    "def vect_tronc_gaussian(n,m, sigma = 1, limit = 1):\n",
    "    A = np.array([tronc_gaussian(sigma, limit=limit) for i in range(n*m)])\n",
    "    return A.reshape((n,m))\n",
    "\n",
    "def generate_series(A,n,sigma= 1, limit_gaussian = 1 ):\n",
    "    M = A.shape[0]\n",
    "    cur =  vect_tronc_gaussian(M,1,sigma= sigma, limit=limit_gaussian)\n",
    "    list_vect = [cur]\n",
    "    for i in range(n-1):\n",
    "        cur = A.dot(cur) + vect_tronc_gaussian(M,1,sigma=sigma, limit=limit_gaussian)\n",
    "        list_vect.append(cur)\n",
    "    return np.concatenate(list_vect, axis = 1)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  II Three estimators\n",
    "\n",
    "### II.A Benchmark : Full Rank Estimator\n",
    "\n",
    "To compare the performance of our model, we use the standard VAR(1) estimator $\\widehat{A}_{std}$ defined as :\n",
    "\n",
    "$$ \\widehat{A}_{std}  =\\underset{A \\in \\mathcal{M}_{M \\times M }(\\mathbb{R})}{argmin} \\sum_{i=0}^n \\|X_{i+1} - A X_i \\|_2^2 $$ \n",
    "\n",
    "If we set :\n",
    "\n",
    "$$ Y = \\left(X_1| X_2 | \\dots | X_n \\right)$$\n",
    "$$ X = \\left(X_0|X_1 | \\dots | X_{n-1} \\right)$$\n",
    "\n",
    "We have then \n",
    "\n",
    "$$  \\widehat{A}_{std} = \\underset{A \\in \\mathcal{M}_{M \\times M }(\\mathbb{R})}{argmin} \\|Y- AX \\|_F^2 $$\n",
    "\n",
    "where $\\|\\|_F$ is the Frobenius norm. Therefore, we can have an exact expression for $\\widehat{A}_{std}$:\n",
    "$$\\widehat{A}_{std} = (X X^T)^{-1} X^T Y$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def empirical_risk(estim,series ):\n",
    "    \"\"\"Empirical Squared L2 risk for a matrix estim\"\"\"\n",
    "    n = series.shape[1]-1 \n",
    "    return( 1/n * sum([np.linalg.norm(series[:,i+1] - estim @ series[:,i])**2  for i in range(n)]) ) \n",
    "\n",
    "def std_estimator(series):\n",
    "    \"\"\" Full rank estimator\"\"\"\n",
    "    Y = series[:,1:]\n",
    "    X = series[:,:-1]\n",
    "    return  (Y @ X.transpose()) @  np.linalg.inv(X@ X.transpose() ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II.B Low Rank Estimator\n",
    "\n",
    "In this setting, we perform the minimization describe in the paper \n",
    "\n",
    "$$ \\widehat{U}_{low rank}, \\widehat{V}_{low rank}  =\\underset{U,V \\in \\mathcal{M}_{M \\times r }(\\mathbb{R})}{argmin} \\sum_{i=0}^n \\|X_{i+1} - U V^T X_i \\|_2^2 $$ \n",
    "\n",
    "And we have then $\\widehat{A}_{low rank} = \\widehat{U}_{low rank} \\widehat{V}_{low rank} $\n",
    "\n",
    "To perform this minimization, we alternate minimisation on $U$ then $V$. There are also exact expression for this form. When we suppose that we know the underlying rank of the simulation, this estimator is called an 'oracle' estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_U(V,X,Y):\n",
    "    VX =  V @ X\n",
    "    return (Y @ VX.transpose()) @  np.linalg.inv(VX@ VX.transpose() ) \n",
    "    \n",
    "def min_V(U,X,Y):\n",
    "    tUYtX = U.transpose() @ Y @ X.transpose()\n",
    "    return np.linalg.inv(U.transpose()@ U ) @ tUYtX  @  np.linalg.inv(X@ X.transpose() ) \n",
    "\n",
    "def low_rank_estimator(series, r = 10 ,  n_iter = 20, verbose = False ):\n",
    "    \"\"\"Low rank estimator with alternate minimisation\"\"\"\n",
    "    Y = series[:,1:]\n",
    "    X = series[:,:-1]\n",
    "    V_cur = random_unit_matrix(series.shape[0],r).transpose()\n",
    "    for i in range(n_iter):\n",
    "        try :\n",
    "            U_cur = min_U(V_cur, X, Y)\n",
    "            V_cur = min_V(U_cur, X, Y)\n",
    "            if verbose:\n",
    "                print(empirical_risk(U_cur @ V_cur, series))\n",
    "        except:\n",
    "            return('Fail')\n",
    "    return U_cur @ V_cur\n",
    "\n",
    "def low_rank_estimator_by_projection(series, r= 10):\n",
    "    \"\"\"Low rank estimator approached using the projection of the standard estimator on the subset of low rank matrices\"\"\"\n",
    "    tot_estim = std_estimator(series)\n",
    "    U,D,V = np.linalg.svd(tot_estim)\n",
    "    D[r:] = 0\n",
    "    return U @ np.diag(D) @ V\n",
    "\n",
    "def low_rank_estimator_thrd(series, r= 10):\n",
    "    Y = series[:,1:]\n",
    "    X = series[:,:-1]\n",
    "    U,D,V = np.linalg.svd(Y,full_matrices=False)\n",
    "    D[r:] = 0\n",
    "    Yhat = U @ np.diag(D) @ V\n",
    "    return (Yhat @ X.transpose()) @  np.linalg.inv(X@ X.transpose() ) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II.C Penalized Full Rank Estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "Another way to compare is to penalize the rank:\n",
    "\n",
    "$$  \\widehat{A}_{pen} = \\underset{A \\in \\mathcal{M}_{M \\times M }(\\mathbb{R})}{argmin} \\|Y- AX \\|_F^2  + C \\sqrt{rg(A)}$$\n",
    "\n",
    "This minimization is hard to compute in practice, beause we don't have an explicit formula for the low rank minimization. Furthermore it is impossible to use gradient-descent-like algorithm in this case because of the discontinuities.\n",
    "\n",
    "So we have to start from the problem for r:\n",
    "\n",
    "$$\\widehat{A_{lr}}(r) = \\underset{A \\in \\mathcal{M}_{M \\times M }(\\mathbb{R})\\\\ rg(A) = r }{argmin} \\|Y- AX \\|_F^2 $$ \n",
    "\n",
    "\n",
    "We are then trying to compute\n",
    "\n",
    "$$ \\widehat{r_{pen}} = \\underset{r \\in [0,M] }{argmin} \\|Y- \\widehat{A}_{lr}(r) X \\|_F^2  + C \\sqrt{r}$$\n",
    "\n",
    "and we use $\\widehat{A}_{lr}(r_{pen})$ as our objective matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log\n",
    "\n",
    "def penalized_est(series, C = 1, return_rank = False, min_rank = 0 ):\n",
    "    \"\"\"Exact version of the penalized estimator using Early Stopping (take more time)\"\"\"\n",
    "    \"\"\"\"\"\"\n",
    "    \n",
    "    M = series.shape[0]\n",
    "    n = series.shape[1]\n",
    "    err = np.zeros(M-1)\n",
    "    list_mod = list()\n",
    "    curr_err = 0\n",
    "    old_err = float('inf')\n",
    "    for r in range(M-1):\n",
    "        \n",
    "        Ar_1 = low_rank_estimator(series, r=r+1)\n",
    "        curr_err = empirical_risk(Ar_1,series) + C* sqrt((r+1)* log(9*n*(r+1)))\n",
    "\n",
    "        if curr_err >= old_err or r==M-2:\n",
    "            r_min = r +1\n",
    "            break\n",
    "        Ar = Ar_1\n",
    "        old_err = curr_err\n",
    "    if return_rank:\n",
    "        return Ar, r_min\n",
    "    return Ar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II.D Nuclear norm regularisation\n",
    "\n",
    "Another possibility is to use the nuclear norm $\\| A \\|_1 $ to penalize the rank of the matrix [2]. We have then the following problem : \n",
    "\n",
    "$$  \\widehat{A}_{pen} = \\underset{A \\in \\mathcal{M}_{M \\times M }(\\mathbb{R})}{argmin} \\|Y- AX \\|_F^2  + C \\| A \\|_1$$\n",
    "\n",
    "We use the optimisation scheme described in [3].  Let's sketch the procedure noting $f(A) = \\|Y- AX \\|_F^2$. We define :\n",
    "\n",
    "$$ A_k = A_{k-1} - \\frac{1}{t_k}\\bigtriangledown f(A_{k-1}) $$ \n",
    "\n",
    "for a decreasing sequence $t_k$ of weight.The objective can be expressed as :\n",
    "\n",
    "\n",
    "$$ \\frac{t_k}{2}\\|A -\\left(A_{k-1} - \\frac{1}{t_k}\\bigtriangledown f(A_{k-1})\\right) \\|_F^2 + C \\|A\\|_1$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.linalg import svd\n",
    "from math import sqrt\n",
    "\n",
    "def gradf(A,X,Y):\n",
    "    return ( A @ X -Y ) @ np.transpose(X)\n",
    "\n",
    "\n",
    "def step_gradient_nuclear_norm(actual, X, Y, C, t):\n",
    "    new_mat = actual - 1/t * gradf(actual,X, Y)\n",
    "    try :\n",
    "        U, D, V = svd(new_mat)\n",
    "        lmbda = C /t\n",
    "        return U @ np.diag(np.maximum(0,D-lmbda)) @V\n",
    "    except:\n",
    "        return actual\n",
    "    \n",
    "\n",
    "def nuclear_norm_estimator(series,C =10 , n_iter=10, val = None):\n",
    "    Y = series[:,1:]\n",
    "    X = series[:,:-1]\n",
    "    M = X.shape[0]\n",
    "    \n",
    "    curr = np.random.rand(M,M)\n",
    "    for i in range(n_iter):\n",
    "        curr = step_gradient_nuclear_norm(curr, X,Y,C, (i+1)*50)\n",
    "        if not val is None:\n",
    "            print('Iteration : %s' %i)\n",
    "            print(np.linalg.norm(val[1] - curr @ val[0]))\n",
    "            print(np.linalg.norm(curr, 'nuc'))\n",
    "            print('======')\n",
    "    return curr\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true
   },
   "source": [
    "### II.E Slope Heuristics\n",
    "\n",
    "##### II.E.i Slope heuristics for penalized estimator\n",
    "\n",
    "We introduce the slope heuristics proposed by Birgé and Massart[1] to compute a good value for C. The idea is to compute the complexity of the best model $\\widehat{A_C}$ for each value of C.\n",
    "\n",
    "There is normaly a value $\\widehat{C}$ of C such that the complexity (here the rank of $\\widehat{A_C}$) is large if  $C < \\widehat{C}$ and reasonnable otherwise. The slope heuristics proposed to choose $\\tilde{C} = 2 * \\widehat{C} $ as the ideal value for the parameters $C$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Cannot cast array data from dtype('float64') to dtype('<U32') according to the rule 'safe'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-58-07ea1239a411>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mr_obt\u001b[0m  \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mC\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mC_list\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m     \u001b[0mApen\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpenalized_est\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreturn_rank\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m     \u001b[0mr_obt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-57-dd5900eded01>\u001b[0m in \u001b[0;36mpenalized_est\u001b[0;34m(series, C, return_rank, min_rank)\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m         \u001b[0mAr_1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlow_rank_estimator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m         \u001b[0mcurr_err\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mempirical_risk\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mAr_1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m*\u001b[0m \u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m \u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m9\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcurr_err\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0mold_err\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-f23e9e04f68f>\u001b[0m in \u001b[0;36mempirical_risk\u001b[0;34m(estim, series)\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0;34m\"\"\"Empirical Squared L2 risk for a matrix estim\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mseries\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0;32mreturn\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mn\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mestim\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mseries\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m  \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mstd_estimator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-f23e9e04f68f>\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0;34m\"\"\"Empirical Squared L2 risk for a matrix estim\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mseries\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0;32mreturn\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mn\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mestim\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mseries\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m  \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mstd_estimator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Cannot cast array data from dtype('float64') to dtype('<U32') according to the rule 'safe'"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "\n",
    "A = generate_matrix(100,75,1)\n",
    "series = generate_series(A,1000)\n",
    "C_list = [0.02,0.05,0.1,0.2,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,2,5,10]\n",
    "r_obt  = list()\n",
    "for C in C_list:\n",
    "    Apen, r = penalized_est(series,C=C, return_rank = True)\n",
    "    r_obt.append(r)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log\n",
    "plt.style.use('bmh')\n",
    "plt.figure(figsize=(12,7))\n",
    "plt.plot(list(map(log, C_list)),r_obt)\n",
    "plt.xlabel('log(C)')\n",
    "plt.ylabel('rank')\n",
    "plt.title(\"Rank of $A_{pen}(C)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the slope change around C=1\n",
    "\n",
    "In order to find the point where the slope change, we compute the penalized estimator $\\widehat{A_{C_i}}$ for each coefficient $(C_i)_{i \\in [0,p]}$ . Then we compute the rank $r_{(C_i)}$ of each of these predictor.\n",
    "\n",
    "Then we choose:\n",
    "\n",
    "$$\\widehat{j} =  \\underset{i \\in [1,p-1]}{\\text{argmax}} (r_{C_{i+1}} - r_{C_{i}}) - (r_{C_{i}} - r_{C_{i-1}}) $$\n",
    "\n",
    "and we have $\\widehat{C} = C_{\\widehat{j}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ideal_C(C_list, series):\n",
    "    \"\"\"Find the ideal penalisation coefficient for the penalized estimator among C_list\"\"\"\n",
    "    r_obt  = list()\n",
    "    for C in C_list:\n",
    "        est = 'Fail'\n",
    "        while est == 'Fail':\n",
    "            est = penalized_est(series,C=C, return_rank = True)\n",
    "        Apen, r = est\n",
    "        if r == 1:\n",
    "            #Case where the ideal rank is 1\n",
    "            print(\"Low Rank\")\n",
    "            return 1\n",
    "        r_obt.append(r)\n",
    "    ser = pd.Series(r_obt)\n",
    "    j = np.argmax(ser.diff(1).diff(1).fillna(float('-inf')).values) - 2\n",
    "    if C_list[j] < 1:\n",
    "        return 2*C_list[j]\n",
    "    else:\n",
    "        return C_list[j]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ideal_C(C_list,series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### II.E.ii Slope heuristic for Nuclear norm estimaor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "\n",
    "A = generate_matrix(100,75,1)\n",
    "series = generate_series(A,1000)\n",
    "C_nuc_list = [1,2,10,20,100,200, 1000,2000,10000,20000]\n",
    "r_obt  = list()\n",
    "for C in C_nuc_list:\n",
    "    Apen = nuclear_norm_estimator(series,C=C, n_iter = 100)\n",
    "    r = np.linalg.norm(Apen, 'nuc')\n",
    "    r_obt.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log\n",
    "plt.style.use('bmh')\n",
    "plt.figure(figsize=(12,7))\n",
    "plt.plot(list(map(log, C_nuc_list)),r_obt)\n",
    "plt.xlabel('log(C)')\n",
    "plt.ylabel('rank')\n",
    "plt.title(\"Rank of $A_{pen}(C)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ideal_C_nuc(C_nuc_list, series):\n",
    "    \"\"\"Find the ideal penalisation coefficient for the penalized estimator among C_list\"\"\"\n",
    "    r_obt  = list()\n",
    "    for C in C_nuc_list:\n",
    "    \n",
    "        Apen = nuclear_norm_estimator(series,C=C, n_iter=100)\n",
    "        r = np.linalg.norm(Apen, 'nuc')\n",
    "        r_obt.append(r)\n",
    "    ser = pd.Series(r_obt)\n",
    "    j = np.argmax(ser.diff(1).diff(1).fillna(float('-inf')).values) - 2\n",
    "    if C_nuc_list[j] < 1:\n",
    "        return 2*C_nuc_list[j]\n",
    "    else:\n",
    "        return C_nuc_list[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III Experience"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def experience(M ,n,r,a , nb_iterations = 1000, nb_test= 100,C_list=np.geomspace(0.01,10,100)):\n",
    "    list_res = list()\n",
    "    \n",
    "    #Slope heuristics\n",
    "    A = generate_matrix(M,r,a)\n",
    "    series = generate_series(A, n + nb_test)\n",
    "    \n",
    "   \n",
    "    C=ideal_C(C_list, series)\n",
    "    if r > 40:\n",
    "        # For high value of ranks, we observe that a smaller value of C is better\n",
    "        C= C/2\n",
    "    lmd = ideal_C_nuc(C_list,series)\n",
    "    \n",
    "    for iteration in range(nb_iterations):\n",
    "        A = generate_matrix(M,r,a)\n",
    "        series = generate_series(A, n + nb_test)\n",
    "        train_set, test_set = series[:,:-nb_test], series[:,-nb_test:]\n",
    "        \n",
    "        A_std =  std_estimator(train_set)\n",
    "        A_low_rank = low_rank_estimator(train_set,r=r) \n",
    "        A_nuc = nuclear_norm_estimator(train_set,C=lmd,n_iter=100 ) \n",
    "        if A_low_rank != 'Fail':\n",
    "            A_pen = penalized_est(train_set, C=C, min_rank=0)\n",
    "            if A_pen != 'Fail':\n",
    "                min_risk =  empirical_risk(A, test_set)\n",
    "                err_1 = empirical_risk(A_std, test_set) - min_risk\n",
    "                err_2 = empirical_risk(A_low_rank, test_set) - min_risk\n",
    "                err_3 = empirical_risk(A_pen, test_set) - min_risk\n",
    "                err_4 = empirical_risk(A_nuc,test_set) -min_risk\n",
    "                list_res.append([err_1,err_2, err_3,err_4])\n",
    "    return np.array(list_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  III.A Rank Importance\n",
    "\n",
    "We fix the number n of observation, and we observe the risk excess for every estimator . We estimate the risk excess using 100 fresh values generated with A. There are some troubles with the matrix inversion, so there are some outliers. I suppose that this problem is caused by numerical problems while inversing quasi-non invertible matrics. I do not consider this matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "#parameters\n",
    "\n",
    "nb_iterations = 100 #Number of simulations for each point\n",
    "M = 100  #Numbers of series\n",
    "n_list = [200,500,1000,2000,5000] #length of series\n",
    "#n_list =  [1000]\n",
    "fontsize = 18\n",
    "reponses = []\n",
    "\n",
    "#We test the following rank\n",
    "ranks = [2,3,5,7,10,15,20,30,50,75,100]\n",
    "#ranks = [20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "hide_input": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rank: 2 \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/CDBDX/remy.garnier/.virtualenvs/python_3.6_TS/lib/python3.6/site-packages/ipykernel_launcher.py:22: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "/home/CDBDX/remy.garnier/.virtualenvs/python_3.6_TS/lib/python3.6/site-packages/ipykernel_launcher.py:24: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rank: 3 \n",
      "Rank: 5 \n",
      "Rank: 7 \n",
      "Rank: 10 \n",
      "Rank: 15 \n",
      "Rank: 20 \n",
      "Rank: 30 \n",
      "Rank: 50 \n",
      "Rank: 75 \n",
      "Rank: 100 \n",
      "Rank: 2 \n",
      "Rank: 3 \n",
      "Rank: 5 \n",
      "Rank: 7 \n",
      "Rank: 10 \n",
      "Rank: 15 \n",
      "Rank: 20 \n",
      "Rank: 30 \n",
      "Rank: 50 \n",
      "Rank: 75 \n",
      "Rank: 100 \n",
      "Rank: 2 \n",
      "Rank: 3 \n",
      "Rank: 5 \n",
      "Rank: 7 \n",
      "Rank: 10 \n",
      "Rank: 15 \n",
      "Rank: 20 \n",
      "Rank: 30 \n",
      "Rank: 50 \n",
      "Rank: 75 \n",
      "Rank: 100 \n",
      "Rank: 2 \n",
      "Rank: 3 \n",
      "Rank: 5 \n",
      "Rank: 7 \n",
      "Rank: 10 \n",
      "Rank: 15 \n",
      "Rank: 20 \n",
      "Rank: 30 \n",
      "Rank: 50 \n",
      "Rank: 75 \n",
      "Rank: 100 \n",
      "Rank: 2 \n",
      "Rank: 3 \n",
      "Rank: 5 \n",
      "Rank: 7 \n",
      "Rank: 10 \n",
      "Rank: 15 \n",
      "Rank: 20 \n",
      "Rank: 30 \n",
      "Rank: 50 \n",
      "Rank: 75 \n",
      "Rank: 100 \n"
     ]
    }
   ],
   "source": [
    "for n in n_list:\n",
    "    list_alpha = [1]\n",
    "    list_rep = list()\n",
    "    for alpha in list_alpha:\n",
    "        list_inter = list()\n",
    "        for rank in ranks:\n",
    "            print('Rank: %s ' %rank)\n",
    "            list_inter.append(experience(M, n, rank,alpha,nb_iterations=nb_iterations))\n",
    "        list_rep.append(list_inter)\n",
    "    reponses.append(list_rep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_n =0\n",
    "i_alpha = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "ab_list =list()\n",
    "type_list = ['full rank', 'oracle','penalized','nuclear']\n",
    "for rank in ranks:\n",
    "    ab_list += [rank]*(4*nb_iterations)\n",
    "\n",
    "dat = np.concatenate(reponses[i_n][i_alpha],axis = 0)\n",
    "df_l = pd.DataFrame(dat)\n",
    "df_l.columns = type_list\n",
    "df_l = df_l.stack()\n",
    "df_l = df_l.reset_index()\n",
    "df_l.drop(columns  ='level_0', inplace= True)\n",
    "df_l['rank'] = ab_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAGwCAYAAAAUiKyAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3X98FNW9//HXyQ8gJBIJifymSBARULRqG6+1tYIJarBUsVhpaV3QohEtKlYtCkijLfXKtZRCf63a3tvalgJKlOSC9UepxGpRaK8ofg1gLSRkAyYQQrLJnu8fSbZJyI8dspv9kffz8dgHzJyZ2c/Mzk72M2fOOcZai4iIiIiIiMSGuHAHICIiIiIiIsGjJE9ERERERCSGKMkTERERERGJIUryREREREREYoiSPBERERERkRiiJE9ERERERCSGhC3JM8b0M8b81Riz0xjzf8aYZU3znzbG7DXGvNP0Oj9cMYqIiIiIiESbhDC+dy1whbX2mDEmEdhmjNncVLbIWrsujLGJiIiIiIhEpbAlebZxFPZjTZOJTa+ARmZ/5ZVXbN++fUMVmoiIiIiISEQ7fvy4Z8qUKRntlcUvXbq0h8P5N2NM/LJly94GHgOetdY+vWzZshnAN5YtW3brsmXLxi1btuyVpUuXNrRcb+fOnUvvvPNOXnzxRZ577jkSEhLIzs6mtraW008/ndNOOw2fz8ewYcNoaGggMTGR4cOHU1tby8CBA0lOTsZay4gRI/B6vfTt25chQ4ZQV1fHoEGDSEpKwlrLyJEjqa2tpX///pxxxhl4vV7S09Pp06cPgL88OTmZjIwMvF4vZ5xxBgkJCRhj/OWnnXYaaWlp1NfXM2TIEIwxxMXF+csHDBhAamoqDQ0NDB06FKBVufZJ+6R90j5pn7RP2iftk/ZJ+6R90j613KeamprdY8aM+Vm7eVZjhVp4GWNOBzYAC4AKoBToA/wM+NBa+0jL5bdv327Hjx/f43GKiIiIiIhEgh07dvxtypQpF7VXFhG9a1prPwFeBqZZaw/aRrXAU8BnwhudiIiIiIhI9Ahn75oZTTV4GGOSgCuB94wxQ5vmGWAG8I9wxSgiIiIiIhJtwtm75lDgGWNMPI3J5u+ttQXGmD8ZYzIAA7wDzA9jjCIiIiIiIlElnL1r7gIuaGf+FWEIR0REREREJCZERJs8ERERERERCQ4leSIiIiIiIjFESZ6IiIiIiEgMUZInIiIiIiISQ5TkiYiIiIiIxBAleSIiIiIiIjFESZ6IiIiIhJTH42H+/PlUVFSEOxSRXkFJnoiIiIiElNvtZufOnbjd7nCHItIrKMkTERERkZDxeDwUFBRgraWgoEC1eSI9QEmeiIiIiISM2+3GWguAz+dTbZ5ID1CSJyIiIiIhU1RUhNfrBcDr9VJYWBjmiERin5I8ERERUccYEjI5OTkkJiYCkJiYyLRp08IckUjsU5InEkT6kSQi0UodY0iouFwujDEAxMXF4XK5whyRSOxTkicSRPqRJCLRSB1jSCilp6eTm5uLMYbc3FwGDRoU7pBEYp6SPJEgafkjadOmTfqRFIA9e/YwZcoUPvjgg3CHEhVUUyyhoo4xJNRcLheTJ09WLZ5ID1GSJxIkbreb+vp6AOrr6/UjKQBLliyhurqaJUuWhDuUqKCaYgkVdYwhoZaens7atWtVi9eJrKysLl8igUoIdwASWQK5gBQXF/dAJNGnsLAQn88HNN4J37x5M4sWLQpzVJFrz5497N27F4CSkhI++OADzjrrrDBHFbnaPk7ncrn0Y0m6pbPrfXV1tb9c13wRkeijmjyRIBk8eHCr6SFDhoQpkujQtvZOtXmd0+N0IiKxrbi42P9qb55uuIgTqsmTVtpeQHQnN3ClpaWtpg8ePBimSKJDcy1es5KSkjBFEh3ae5xONcXSHR1d76+//nqdWwHweDwsXryY/Px81aqLSMRRTZ5IkFx11VX+LqKNMVx99dVhjiiynXnmma2mx4wZE6ZIooPGmZKeoo4xAqM2siISyVSTJxIkLpeLTZs24fV6SUxM1A+lLixcuJA777yz1bR0zOVyUVBQAGicKQkt1Up1TW1kJdKoTwVpK+aTPJ300lPS09OZPn06GzZsYPr06fqD34VXXnnlpOmLL744PMFEgeZxpjZs2KBxpkTCrL02snrEVUQiScwneSI9yeVyUVJSolqWABQVFbWaVhuzrun8EokMaiMbGN1o7znqU0HaivkkTye99KTmcYCkfeqyvXt0fomEV0fXMF2/RHped24iBDrmYDR/n2M+yRMRiVa94Y+QiMSultcmJcEiPUtJnoj0GHXZLiLRrL2kRdcvkfDoztN6veFJPyV5IhJ2amPWvt7wRyiYVPMp4aDrl4hEIiV5It2gH5XBoZ4iRSRa6folIpFISZ6IiMQE1XyKiIg0UpIn0g36USkiIiIikSYu3AGIiIiIiIhI8CjJExERERERiSFK8kRERERERGKIkjwREREREZEYoiRPREREREQkhijJExERERERiSFK8kRERERERGKIkjwREREREZEYoiRPRERERETkFHg8HubPn09FRUW4Q2lFSZ6IiIiIiMgpcLvd7Ny5E7fbHe5QWkkIdwAiIvJvWVlZ3VqmuLg4mOGIiEiIdPd6D7rmh5vH46GgoABrLQUFBbhcLgYNGhTusAAleb1eIBeYrpbrTRcYHS8RiWb6USmhpPNLehu32421FgCfz4fb7WbRokVhjqqRkjwRkQi0detWR8tPnTo1RJFENv2oFJFo5/R6D73zmh+MG+3dWbe9vxVFRUV4vV4AvF4vhYWFSvKMMf2A14C+TXGss9YuMcacCTwLDAL+BnzdWlsX6HZV03JqdIFxRscrMPo+ikQmXcMCo5sIp0bnl/QWOTk5bNq0Ca/XS2JiItOmTQt3SH7hrMmrBa6w1h4zxiQC24wxm4G7gZXW2meNMWuBucCaMMYpIiIRTj8qRaQ36uomQyzeRO3O9X7lypWO1124cGGHZS6Xi4KCAgDi4uJwuVyOtx8qYUvybOMDrMeaJhObXha4Aripaf4zwFJOIcnTH3yRyKHvo4hEM13Dws/j8bB48WLy8/MjpmMLkfT0dHJzc9mwYQO5ubkRdW6GtU2eMSaexkcyxwKrgQ+BT6y19U2LfAwMD1N4IiIiIhIBWnZTHyltniJBsGumxDmXy0VJSUlE1eJBmJM8a20DcL4x5nRgAzA+kPUOHTrE3LlzSUhIoKGhgeuuu468vDxKS0s599xzqa2t5eDBgwwcOJCqqiqstZx++ukcPnyY/v37A3D8+HHS0tL45JNPMMYwYMAAsrOz2bdvHx6Ph+rqaoYMGUJpaSmJiYmkpqbi8XhITU2lrq6Ompoaf3mfPn047bTTqKioYODAgdTU1HDixAl/eb9+/UhKSuLIkSMMGjSIo0ePUldX5y9PSkqiT58+VFZWkp6eTmVlJV6v11+enJxMfHw8VVVVZGRkcPjwYay1ZGRkUFZWRkpKCgDHjh1j8ODBlJeXY4whLS2N8vJyBgwYQENDQ7v7lJ2dzYcffsiRI0c4ceIE6enpeDweEhMT6d+/P5WVlQwYMIATJ05QV1fnL+/Tpw+f+tSnOPvsszl27FhE7VMoP6fMzEwyMzOpqanh2LFj1NfX+49JUlIScXFxVFdXt3vuTZgwAYD9+/dH1D6F8nMaN24cffv25eDBg/7jlJCQQEpKCp988gkpKSnU19e3e+5lZ2fz/vvvU1ZWFlH7FOrPafjw4YwcObLVMevXrx8JCQkcO3aM008/vd1zb9y4cYwePZrjx49H3D6F8nPKzs7mzTffpKyszH8tP3LkCMnJyfh8Pmpqajo89y644AIyMjKora2NqH0K5ed00UUXkZaWhtfr9V/L+/XrR1VVFampqRw/fhyv19vuuTdkyBAyMzPZv39/RO1TqD+nCRMmUFZWFtDviJbnXnZ2NsXFxezfvz/i9ilUn1PzPh88eDCg3xEtz72srCwGDBjg/z623CdjTKtu6rOzsxk+fHjUn3uXXnopycnJ/u9jV78jWp57AwcO5IILLsDn82GMwVqLMab597V/HtCq3Frb6tyMpnNvwIABZGVlUV5eHtDviJbnXvM++3y+do9Jy3lty5s/p+Zzs+0+nThxgu9+97skJyezf//+Hj33OmOadybcjDEPAzXAd4Ah1tp6Y8wlwFJrbU7LZbdv327Hj28/H2x+9rg7j1ZE6zPKp0LHy5nu9NjUrDceL51fgTvVY6bjpXMsEDpezuh4OROq47VixYpWnVtce+21MVGbF4zj1Z2avGg7N3W8TrZjx46/TZky5aL2ysLZu2YG4LXWfmKMSQKuBH4AvAzMpLGHzW8Az4Urxt5EbQdEJJrpGiahpPMrvCK5m/pg0PkloRDOxzWHAs80tcuLA35vrS0wxrwLPGuM+R7wNvDLMMYo0i41wheRaKbrkUSTSO6mXnped65fvak9Yjh719wFXNDO/BLgMz0fUe+mpEVCSeeKczpmzugaJqGk88uZYO97JHdTHww6vyQUwtrxSijp5BcREWmfflRKNInkbuql50XSOHmRLGaTPBGJHPpB6dypdrwiIhJuobjmR2o39dL7BNoJX7g7aonZJE8/KkVERERiQ3p6OmvXrg13GBEnWmuZJPRiNskTEREREREJprY1dM01e+GuuWtLSZ6ISATSkwUiItKV3tTGTJxRkidyCvQD3BkdL5HIou+kiEQrDaEQGCV5IiIRpLPHPSL1kZBIoKRFQknnl4SSzq/ey+PxsHjxYvLz84Pea6ySPAF0gQlUVz+u9SO8NR0vkciimwgiEq1C+ZsiXNc/t9vNzp07cbvdLFq0KKjbVpInIiJRS0mLhJLOr1MTjBvHgXRTH+3HvjvnV6Dd+Evk8ng8FBQUYK2loKAAl8sV1No8JXm9nGpaRERERER6ltvtxloLgM/nC3ptXswmeXr8UERERKTnBLPmM1q6qRc5VUVFRXi9XgC8Xi+FhYVK8kREREREerve1FtkrMnJyWHTpk14vV4SExOZNm1aULcfc0meHj8UEREREZFI5nK5KCgoACAuLg6XyxXU7cdckiciIiIiEstUqRH90tPTyc3NZcOGDeTm5gZ9CIW4oG5NRESkDY/Hw/z586moqAh3KCIiIhHD5XIxefLkoNfigZI8EREJsZbjAImIiEij9PR01q5dG/RaPNDjmiIiEkKhHgdIREQk1AIZl7CzZcLx2KySPBERCZlQjwMkIiLdE+jA6mrfF12U5ImISMiEehwgERGRnjJz5kxHy69bty5EkXRNSZ6I9JiO7ha2na+7hbEj1OMAiYhI92jg+dikjldERCRkXC4XxhggNOMAiYiIyMlUkyciPUZ3BXufUI8DJCIiIidTkiciIiHlcrkoKSlRLZ5IBPB4PCxevJj8/HzddJGoFkgTkI5uLveG5iNK8kREJKSaxwEKtd7wRzuYdLx6p5bjVqoTJBFnwtmRilNK8kRERER6AY1bKbGkOzegesPNKyV5IiISE2Llj3ZPPU4XK8erJ7VX+xlNNZ+9ZdzKaB33LdrPr95AQyiIiEi36XG63kmP00moRPK4lYEmZp0tp2th79adG2Sx2FZVSZ5IN+hHuIgEkx6ni2zRfi3vLeNWRuu4b5EeX6Trzg2yWLy5piRPRCRC6Q9+7xPOx+li8U52KEXj8XK5XBQUFACRO27lypUrHa+zcOHCEEQi0aQ7N8hi9eaaBkMX6Ybi4uKAXtI+j8fD/PnzqaioCHcoIhGhvcfpekrLO9nStWg8Xs3jVhpjNG6lxJT2bpD1xLqRTDV5IhI20fp4RDTewZfoEK7H6WL1TnaoRPPx6slxK9WRSM/pzphxsaA77U2drBtNQyioJk9EwqLtj6Roqs2Lxjv4Eh1cLhfGGKBnH6eL1TvZoRLNx6t53MpoSUpFApGTk0NiYiKA4xtk3Vk3ksV8TV5vv7MhEqmitSvvaL6DL5Gv+XG6DRs29OjjdJHc62Ik0vEKjH5f9Zzefqy70940kHU7O76R2rGPavJEJCzC2faoO6L5Dr5EB5fLxeTJk3u0U4xYvZMdKjpeIu0LV1v77rQ3jdW2qjFfkxdpWXWkU82n9JRo7cpbd/Al1Jofp+tJ0dDrYiTR8RJpXzjb2nenvWlPtlXtKarJE5GwCFfbo+7SHXyJRbF6JztUdLxEThbutvbdaW8ai21VY74mT5xRLZ30lHC1Peou3cGXWBWLd7JDScdLpLVobWsfq4Ka5BljTrPWHg3mNkUkdkXjj6RoTU5FuhKOx0SjmY5XaGhg8+il5gyRJeDHNY0xT3ZRfhpQ1O2IRKTXiNbHI8LRMYaIiEgkU3OGyOKkJm+BMeZja+0P2xYYY/oDm4HzgxaZiIiIiPQqK1eudLyOav+c6aiTPSfLtde8R80ZIouTjlceAr5vjJndcqYxJgl4EbgIuD6IsYmIRCQNhi4iItKaOiSKLAHX5Flr840xI4BfGmPKrLVbjTH9gE3AJcBMa+3mUAUqIhIJNBi6iIjEgpkzZzpeZ926dZ2WR2Nb+1jldAiFPBofy/yjMeY/gI3A54GbrLWbgh2ciEik0WDoIiIi7YvWtvaxyFGSZ631AV8F/g68BlwBfM1a+8cQxCZySjweD/Pnz+/x8Vmkd2iv9zARERGRSNJhkmeM+Xx7L+AzwH8C1cBTQGmb8oAYY0YaY142xrxrjPk/Y8xdTfOXGmP+ZYx5p+l1dTf3UXoZtZeSUFLvYSIiIhLpOmuT9wpgOyk3wC3AvBbTFogP8L3rgXustTuahl/4mzFmS1PZSmvt4wFuR8RP7aUk1NR7mIiIiES6zpK8m0P5xtbag8DBpv8fNcbsBoaH8j0l9rXXXkoDcUowaTB0EZHQ0XAIIsHRYZJnrX2mp4IwxowGLgDeAC4F7jDGzAHeorG270jL5Q8dOsTcuXNJSEigoaGB6667jry8PEpLS0lOTiY+Pp6qqioyMjI4fPgw1loyMjIoKysjJSUFgGPHjjF48GDKy8sxxpCWlkZ5eTkDBgygoaGB6upqhgwZQmlpKYmJiaSmpuLxeEhNTaWuro6amhp/eZ8+fTjttNOoqKhg4MCB1NTUcOLECX95v379SEpK4siRIwwaNIijR49SV1fnL09KSqJPnz5UVlaSnp5OZWUlXq/XX659CnyfCgsLW7WXevHFF7nzzjujep9i8XOK9n2aMWMGu3fvZvbs2ezfvz8m9ikWPyftk/ZJ+xQ9+5SVlcWAAQMoLi4mKyuL8vJyqqqqyMzMZNeuXZSWlgIwYMAAsrKyOHjwILW1tYwePZq3336bCRMm8M9//rPdfWr28ccf63OqqyM7O5vi4mJGjBhBXV0d1dXVDBw4kE8++YSkpCT69u1LWVkZgwcPpra2lpqaGk4//XSOHDni/5xqa2sjap/C+Tk1O3HiRI/vU6f5VXOtR3cYY/paa2tPcd0U4FUg31q73hgzGPDQ+OjncmCotbbV81Dbt2+348eP727YEoNWrFjBpk2b8Hq9JCYmcu2116omT0REJMo1D8zd3iDcoVw3FjUfj+4MoaBj+W/hPL927NjxtylTplzUXlnAvWsaY64yxixtM+92Y0wVUG2M+Y0xJtFJYE3L/xH4H2vtegBrbZm1tqGpJ8+f09jRS6+m3iID53K5MMYAai8lIiIiIr1TwIOhA4uAQ80TxphzgCeBD4G9wCzgr8B/BbIx0/hL/JfAbmvtEy3mD21qrwfwZeAfDmKMSS17i1StVOfUXkpEREREQqW55q6r+eGu7XQyTt45NLaRazYLqAE+Y629Cvgd8A0H27sU+DpwRZvhElYYY/5ujNkFfBHo1S1w2/YWqdq8rrlcLiZPnqxaPBERERHplZzU5A2ksa1cs6nAn6y1VU3TrwABj2lnrd1G47ALbb3oIKaYp94inUtPT2ft2rXhDkNEREREYkzbGro9e/Zw2223sXbtWs4666wwRXUyJ0meB/gUQNO4dhcDD7YoTyTwMfIkQEVFRa16iywsLFSSJyIiIiLd0tyJinTPkiVLqK6uZsmSJfzmN78Jdzh+Th7X3A7MN8bMpLHdXQKwuUX5WJrGvZPgycnJITGxsT+bxMREpk2bFuaIRERERERkz5497N27F4CSkhI++OCDMEf0b05q8pYALwO/b5p+xlr7Lvg7UflyU7kEkcvloqCgAFBvkSIiIiISHN0ZQkEaLVmy5KTpSKnNCzjJs9a+29Sj5qVApbX2tRbFpwMraWyXJ0Gk3iJFREREOtZRb4dOlgt3T4gSnZpr8ZqVlJSEKZKTOanJw1p7GNjUzvwjNA6nICHgcrkoKSlRLZ6IiIiISIQ488wzWyV6Y8aMCWM0rTlK8iQ81FukiIiISOf0+GH4eTweFi9eTH5+fq94+mzZsmXMmTOn1XSk6DDJM8bsBXzAeGut1xgTSP2jtdZmBi06EREREREJulAkuG63m507d/aaIb/GjRvnr80bM2ZMRA2h0FnvmvuBjwDbNP1R07zOXh+FLFIREREREYlIHo+HgoICrLUUFBRQUVER7pB6xLJly0hOTo6oWjzopCbPWnt5Z9MiIiIiIhJduupkprmDGqed0bjdbqxtrBvy+Xy9qjbvpZdeCncYJwmoTZ4xJgm4AXjfWvtGaEMSERERkd6so54wW85Xj5iRpaioCK/XC4DX66WwsLBXJHmRKtDB0GuBXwAXhDAWEREREUc8Hg/z58/vNY+GiUSqnJwcEhMTAUhMTGTatGlhjqh3C6gmz1rrM8Z8BAwIcTwiIiIiAettHT30Fqqliz4ul4uCggIA4uLiNPRXmAVakwfwDPB1Y0zfUAUjIiIiEqje2tGDSCRKT08nNzcXYwy5ubm9YgiFSOYkyXsdqAfeMcYsMMZMM8Z8vu0rRHGKiIiItNJeRw8iEj4ul4vJkyerFi8COEnytgCTgbOBJ4EXgJdbvF5p+ldEREROkdqYBa69jh6ihT5nEQklJ0nezW1erjav5nkiIiJyilq2MZPORXNHD/qcJRbpvHYmlDd7Aup4BcBa+0zQ311ERET82rYxc7lcatfSiWjt6EGfc2isW7cu3CH0ajqvnQtlx1FOavJEREQkhNTGzJlo7ehBn7PEIp3XzoS646iAa/JEREQktDSYsHMul4uSkpKoqcUDfc6hMnPmTMfrqPYveHReO9NeUhzM46WaPBERkQgRzW3MwiU9PZ21a9dGTS0e6HOW2KTz2plQdxylJE9ERCRCuFwujDFAdLUxE2fC+TmrV08JFV2/nAl1UqwkT0REJEJEaxsz6VpWVpb/lZubS11dHQC1tbVcc801ZGVl9Ugc6v3QGSXFgdP1y5lQJ8VqkyciIhJBorGNmUQH9X7o3OrVq3nnnXf4yU9+wkMPPRTucEKioxsMbecXFxd3uS1dvwLXnBRv2LAhJElxwDV5xpg0Y8x5nZSfZ4wZGJywREREeqdobGMmXSsuLm71am9+qKn3Q2c8Hg9FRUUAFBYWqjYvALp+OTNjxgz69+/PjBkzgr5tJzV5K4BPN73a8xTwJjC/u0GJiIiISHCp90NnVq9ejc/nA6ChoSFma/N64gaDtG/jxo0cP36cjRs3hnWcvC8Cmzopfx6Y2r1wREQk1qhNi0hkyMnJISGh8f5+QkKCej/swpYtW1pNN9fqiQRDqMfJc5LkDQM+6qT846ZlRERE/NTRg0hkcLlc/popn8+ndlMSdLqpF7hQPz7tJMmrBj7VSfmngNruhSMiIrEk1HcqReTUNPfqJx278sorW03n5OSEKZLooZt6gYukcfLeAL5hjDmtbUHTvDnAX4MVmIiIRD919CASOdxuN3FxjT/9jDH6Pnbhq1/9aqvpG2+8MUyRRAfd1HMmksbJexwYAbxujJlpjBnb9JoJvN5U9sOgRiciIlEt1HcqRfR4WOCKioqor68HoL6+Xt/HLmzcuNFf42mMYePGjWGOKLLppp4zoR4nL+Akz1r7MnA7cBbwO+D9ptfvmubdYa3dGtToREQkqoX6TqWIHg8LnL6PzhQVFfmTFmutkuIu6KaeM6EePN5JTR7W2p8CmcC9wNqm191AprV2bVAjExGRqBfqO5XSu+nxMGf0fXRGSbEz6r3VuVCOk+coyQOw1v7LWrvSWpvX9HrSWvuvoEcmIiJRL9R3KqV30+Nhzuj76IySYmda9t5qrdXxCsBvf/tbqqurefbZZ4O+bcdJnoiIiBMul4vJkyfrD74EnR4Pc07fx8Clp6czZcoUAKZOnaqk2IHmmy/SMY/H4x97cfPmzUF/EiGhowJjjBuwwK3W2oam6a5Ya+3coEV3iurr6zl+/Hi4w5Ag6t+/v/8RABGJLunp6axdqyf6JfhycnLYtGkTXq9Xj9MFSN9HCZXm3lt9Ph9xcXG43W4WLVoU7rAi1urVq1uNW/mTn/yEhx56KGjb7+xX8zdpTPJuAxqaprtigbAmefX19Rw7dozU1FSNARMjrLVUVlaSkpKiRE9ERPxcLhcFBQWAHqeT4PN4PLz00ksAbN26ldtvv121eZ1or/dWJXkd27JlS6vpoqKioCZ5HT6uaa2Ns9bGW2vrWkx39YoPWmSn6Pjx40rwYowxhtTUVKqrq8MdioiIRBC1MZNQUptPZ9RRTWSJyTZ5SvBijzFGn6uIiJxEbcwkVNTm0xl1VOPMlVde2Wo6JycnqNsPOMkzxpQYY67tpDzXGFMSnLBEREREutbcxky1eBJsqplyRjXrzuTl5REX15iKxcfHc/vttwd1+04aOI0GUjopTwY+1a1oRERERGJAVlZWt5crLi4OVjhyCtTm0zmXy0VJSYmOVQDS09PJyclh8+bNTJs2LbyDoXdhMKAuLU+Rz+dj4cKFZGZmkpaWxrZt2wJaLy8vjxtvvLHD6Z6WlpbGc889F7b3FxEREQkG1Uw5p5p1Z/Ly8jj//PODXosHXdTkGWM+D1zeYtZ1xpix7SyaBtwIvBO80IJrxowZlJaW9tj7DRkyhI0bNwa8/JYtW/jNb37D888/z+jRoxk4cGAIoxMREZGesHLlSsfrLFy4MASRyKlQzZSEUiiHNOnqcc0vAkua/m+B65pe7fl/QMRelUpLS9m6dWtS24KOAAAgAElEQVSPvd/UqVMdLV9SUsLgwYP57Gc/G6KIOtY8vlCs8Hg8LF68mPz8fN1JCpCOmYSSzi8RiVYaV1CiVVePa/4XcCYwBjDAt5umW75GA+nW2nHW2rcCfWNjzEhjzMvGmHeNMf9njLmraX6aMWaLMeaDpn9jvkorLy+P7373u3z88cekpaUxefJkAKZPn85999130rLdeRxz27ZtpKWlsWXLFqZOncrgwYP505/+xN69e5k9ezbjx49nxIgRXH755RQVFbVad/LkyTz++OMsXLiQUaNGMXHiRH70ox91+n5PPvkkY8eO5c033zzlmJ1yu93s3LlTXR07oGMmoaTzS0LJ4/Ewf/58Kioqwh2KiEiXsrKyunwFQ6dJnrW20lq731q7j8Zavd82Tbd8fWStPXwK710P3GOtnQBkAXnGmAnA/cBL1tqzgJeapmPaY489xqJFixg2bBi7d+/2D7wZSkuXLuXBBx/kjTfe4MILL6S6upqpU6eyfv16XnvtNaZPn86cOXPYs2dPq/XWrFnDhAkTeOWVV7jrrrtYunQpf/3rX0/avrWWhx56iJ/97Gds2rSJiy++OOT7BI1/7AsKCrDWUlBQoD/6AdAxk1DS+SWhppsIIiInC7jjFWvtq9ba8mC9sbX2oLV2R9P/jwK7geHAl4BnmhZ7BpgRrPeMVAMGDCAlJYX4+HgGDx5Menp6yN/zO9/5DldccQWjR48mPT2dSZMmcfPNNzNhwgTGjBnDPffcw3nnncfzzz/far0vfvGL3HLLLYwZM4Zbb72VMWPG8Nprr7VapqGhgTvuuIPCwkIKCws555xzQr4/zTRwqXM6ZhJKOr8klHQTQUSiTXFxsf/V3rxg9arrZAgFjDEJNCZdnwUGcnKSaK21c50GYYwZDVwAvAEMttYebCoqpbHXzlYOHTrE3LlzSUhIoKGhgeuuu468vDxKS0tJSkqif//+NDQ0+MvDxefz4fV6iYuLIy4ujvr6en9M1loSExP95T6fD2stDQ0N/vLmaa/XizGGhoYGfD4fPp+Purq6Vtuqr6/3/7+5PD4+nvr6euLj4/3zfT4fABMnTqShocFffuzYMR5//HG2bt1KWVkZXq+X2tpaxo8f739fay0TJkygtrbWv09nnHEGhw4doq6uzv9D7uGHHyY+Pp7CwkLS0tJa7VPLfW7ep5afU0JCAl6vl/j4eKAxYUxMTKS+vp7jx4+TkJBAeXk5AwYMoKGhgerqaoYMGUJpaSmJiYkUFha2Grj0xRdf5M4776S0tJQ+ffpw2mmnUVFRwcCBA6mpqeHEiRP+9fv160dSUhJHjhxh0KBBHD16lLq6On95UlISffr0obKykvT0dCorK/F6vf7y5ORk4uPjqaqqIiMjg8OHD2OtJSMjg7KyMlJSGkcgOXbsGIMHD6a8vBxjDGlpaZ3uU2pqKh6Ph9TUVOrq6qipqfGXB2OfOjtm0bpPsfg5Res+bd68udX5tXnzZr7yla9E9T7F4ucUrfv061//2v93zefz8eSTT5KXlxcR+5SdnU1xcbE/PmMM1lr/vy3ntS2/6KKLSEtLo7a2NiY+p5449yZMmMCwYcOoqKggIyODY8eOAZCSkkJ5eTmDBg3CWsuRI0dIT0/n6NGjxMfH+z+n/fv3R9w+xeLnpH1qvU/NDhw4cEr71BnTfHHpijEmDXgZmERj+zzb9C8t/m+ttfEBbfDf200BXgXyrbXrjTGfWGtPb1F+xFrbql3e9u3b7fjx49vdXlVVFQMGDDhpflZWVo93vOIkE1+1ahW/+MUv2Llzp3/el770JcaNG8cPf/hD/7xbb72Vqqoqnn32WaCxjV5FRUWH021t27aNa6+9lg8++KBVBwj33HMPL730Eo888giZmZkkJSVx2223MXbsWFavXg00tsmbN28eCxYs8K83ffp0zjnnHFasWAE0DqFw0003sX79eh5//HFuuummgI9BVzr6bFtasWIFmzZt8ncmc+2117Jo0aKgxRCLVqxYwfPPP++/CfGlL31Jx6yNQJ6P13hW7dN3UkJpypQprX4oJScn90iTh0AEo11NMK4rgcYRzdew5n2cOXOm43XXrVsHRPf+S/RqPndP9fzbsWPH36ZMmXJRe2VOxsn7HjAemAdk0pjU5QDnAL8F3gQcdZtmjEkE/gj8j7V2fdPsMmPM0KbyocAhJ9uMJYMGDaKsrKzVvH/84x8hea/i4mJmzZrFtddey8SJExk2bBj79u07pW1deeWVPPXUU9x777389re/DW6gXXC5XBjTeO9BA5cGxuVy+e80W2t1zCSo9J2UUMrJyfH3Dp2YmMi0adPCHJGISGRw8rjmNcCvrLVPGWOak7kGa+37wNeMMa8AjwG3BbIx0/hX/5fAbmvtEy2Knge+AXy/6d9eO7L25z//eR588EE2b97M2LFjefrpp/nXv/7FqFGjgv5emZmZvPDCC1x99dUkJibygx/8oMtq4M7k5OTw1FNPcfPNN2OM6bEB2psHLt2wYYMGLj0Fgdbs9zZt77B1985bb6LvpISSy+WioKAAiNybCOEeJ0/XL5HeyUlN3hAaa+ugsWdMgH4tyjcC1zrY3qXA14ErjDHvNL2upjG5u9IY8wEwtWm6V5o9ezazZ89mwYIFXHXVVaSkpJCbmxuS98rPzycjI4NrrrmGG264gYsuuohLLrmkW9vMycnB7XZz9913d/j4aCi4XC4mT54ckX/sI5Hb7SYurvFSEBcXp44xJOj0nZRQab6JYIzRTQQRkRac1OQdBpKb/n8U8AIjW5R7aeyMJSDW2m38u01fW1McxBWQIUOGOB6gvLvv58SCBQtatXWDxkdPfvjDH7Zqk9dWc3u5jqbb+tznPsfhwyePeDFy5Eg2bNhwUkwttWwv2GzTpk2tpttue9q0aRw4cKDTmIJNA5c6U1RURH19432b+vp6CgsL1WZKgkrfSQkll8tFSUmJbiKIiLTgJMnbA0wAsNb6jDFvA980xjwNxANzgJKgRxgkGzduDHcIIhEpJyenVccYatMiItFENxFERE7mJMn7X+BeY8wd1tpa4AngWRpr+CyQBNwa/BBFJJSioU2LiEhHPB4PixcvJj8/PyIf1wxm+zoJvd7QG6n0Dk6SvEeBx5sSPKy1vzfG1ANfAxqAddba34Ugxl7P6/Vy4MABhg8fTkKCo6ENRbqkjjFEJJq53W527tyJ2+3Wo+YiPUxJceQKOGOwjd3u1baZtx5Y3/4aEiwVFRUcP34cj8fjuK2fSCDUpkVEopHH46GgoABrLQUFBbhcroi5UdXVj1r1chmZ1BupxIqAkrymAcsrgaXW2uWhDUla8nq9fPLJJwB88sknpKenqzZPgk5tWkQkGrndbv/QLz6fT7V5Ij1MSXHkCmgIBWvtMeATevHA5OFSUVHRatrj8YQpEhERkchSVFSE1+sFGm+KFhYWhjkiEZHI4GScvJeBL4QqEGlfZWWl/y6ltZbKysowRyQiIhIZcnJySExMBFDvwCIiLThJ8hYBnzPGLDPGDAhVQNJaamoqxjQOJ2iMITU1NcwRiYiIRAaXy+X/G6negUVE/s1JkvcS0A9YDBwxxpQaY0ravD4MTZi916BBg/w1edDYdkpERET+3TuwMUa9A4uItOCkB4+PaBwPT8KkZbIXLB999BHnn38+L730EhdccEHQty8iIhJK6h1YRORkToZQuDyEcYTcl770JcrKynrs/QYPHsxzzz3X7e1UVFRgjMFaizFGwyiISMQLZNwk9bwmwaLegUVETtZr+uIvKytj5cqVPfZ+CxcuDMp22ut4pTnJ83q9/gbn8m8ej4fFixeTn5+vR3dEREREpNdx0iZPQqy2tpYHHniAs88+myFDhvC5z32Ot99+G4AdO3Zw6aWXsm3bNi699FLOOOMM/vSnP7F3715mz57N+PHjGTFiBJdffjlFRUWttltXV8fy5cs577zzGDJkCBdccAE//elPO4zjvffeY9asWYwaNYpx48Yxb968Hq0F7S63283OnTtxu93hDkWkVyouLm71am++/FtWVlZALxERkUB1mOQZY3zGmAaHr/qeDD7WLFmyhI0bN7Jq1SqeeuopMjMzueeee1qNjbdmzRpuueUWfvOb33DhhRdSXV3N1KlTWb9+Pa+99hrTp09nzpw57Nmzx7/O7bffzrPPPsv3vvc9iouL+dGPftRhL52lpaXk5uZyzjnnsGXLFjZs2EB1dTWzZ8/G5/OF/Bh0l8fjoaCgAGstBQUFJ40zKCIiIiIS6zp7XPNXnNzRyoXAJOB9YHfTvAnAOOAfwN+CHWBvUV1dzVNPPcWTTz5JdnY2ANnZ2Vx88cWsX7+eiy66CICHH36Ya6+91r9eeno6kyZN8k/fc889FBYW8vzzz3Pvvffy4Ycfsn79en7/+98zdepUAEaPHt1hHG63m0mTJrF06VL/vDVr1jBmzBjefvttLrzwwiDudfC53W7/460+nw+3282iRYvCHJWISMfa1mw219qpxlNERE5Vh0metfabLaeNMVcCM4EZ1trn25TNAH4N3B2CGHuFffv24fV6+exnP+ufFx8fz8UXX8zevXv9SV7bHjCrq6tZsWIFRUVFlJWVUV9fz4kTJ5g4cSIAu3btIi4ujssuuyygOHbu3Mnrr7/OyJEj240x0pO8oqIivF4v0NhmsbCwUEmeiIiEXUeP3Lacr8ReRILFSccry4Gftk3wAKy1G40xPwO+B2wNVnDSqHmgV4D+/fu3Knv44Yd56aWXeOSRR8jMzCQpKYnbbrvNn+g45fP5yM7O5pFHHjmpLCMj45S22ZNycnLYtGmTv1OaadOmhTskEREREZEe5STJOw94ppPy/wfM7144vdfo0aPp06cPb7zxBmeeeSYADQ0NvPnmm3zhC1/ocL3i4mJmzZrlf4TzxIkT7Nu3j7FjxwJw7rnn4vP5+POf/+x/XLMzkydPZuPGjYwcOTIqe+50uVwUFBQAEBcXp3GTREQkIqiWTkR6kpPeNY8A2Z2UTwMquxdO75WcnMzNN9/MsmXL2LJlC++//z733HMP5eXlXHfddR2ul5mZyQsvvMDOnTt59913+da3vsWJEyf85WPHjmXGjBncddddPP/88+zfv5/t27fzu9/9rt3tzZ07l6qqKubOnctbb73Fvn37eOWVV/j2t7/N0aNHg77fwZaenk5ubi7GGHJzczWEgoiIiIj0Ok5q8n4D3GOM+SXwONDcfeM4YBGQCzwR3PB6l+bOTu644w4qKys599xz+cMf/kBqaiofffRRu+vk5+dz5513cs0115Camsr8+fOpra1ttcyaNWt49NFHeeCBB6ioqGDYsGHcdttt7W5v6NChbN68mUceeYQbbriB2tpaRowYwRe/+EX69u0b1P0NFZfLRUlJiWrxOqHBqkUkWgU6nISuYSLSmzlJ8hYDY4GbgW8Czf3pxwEG2NS0TEQaPHhw0AYoD/T9nOrbty+PPfYYjz32WKv5u3fv5tOf/jSHDx8+aZ2RI0eyYcOGVvMWLFhw0naXLVvGsmXLTlp/1KhRJ203MzOTZ57p7MncyJaens7atWvDHYaISEzTzSIRkcgVcJJnra0FvmyMyQa+BIxpKioBnrPW/m8I4gua5557LtwhiEQMddkuItFK1y8Rka45qckDoCmZi+iETkREREKrZVKlREtEJLI46XhFREREREREIpyjmjxjTAIwA/gsMJCTk0RrrZ0bpNhERERERETEoYCTPGNMGvAyMInGjlZs07+0+L8FlOSJiIiIiIiEiZPHNb8HjAfmAZk0JnU5wDnAb4E3AQ1KJiIiIiIiEkZOkrxrgF9Za58CqprmNVhr37fWfg2oAR7rcG0REREREREJOSdJ3hAaa+sA6pv+7deifCNwbTCCEhERERERkVPjpOOVw0By0/+PAl5gZItyL42dsYj0OA3KKyIiIiLSyElN3h5gAoC11ge8DXzTGNPXGNMfmEPjwOgSxfLy8rjxxhs7nA6FG2+8kby8vJC+h4iIiIhIb+GkJu9/gXuNMXdYa2uBJ4Bnaazhs0AScGvwQwyO6dOnU15e3mPvl5GRwaZNm3rs/ULlsccew1ob7jC6pEF5RUREREQaOUnyHgUeb0rwsNb+3hhTD3wNaADWWWt/F4IYg6K8vJyZM2f22PutW7eux94rlAYMGBDuEERERERExIGAH9e0jWrbzFtvrb3OWntDJCd40WD69Oncfffd3H///Zx55pmceeaZPPzww/h8PgC8Xi9Lly5l4sSJDB8+nClTpvDSSy/519+2bRtpaWm8+uqrTJ06leHDh3PFFVewc+dO/zKHDx9m3rx5TJw4kWHDhnHJJZfwP//zP53G1fJxzeb3aPuaPn26f/k33niD3Nxchg8fzsSJE7nnnnuoqqrylx8/fpy8vDxGjhzJ2WefzRNPPBGU4yciIiIiIo2ctMmTEFu3bh3WWoqKinjiiSd45plnWLNmDQD5+fn85S9/4ec//zl/+ctfuPHGG7npppv4xz/+0Woby5cvZ8mSJbz88ssMHDiQW2+91f+45YkTJzjvvPN49tlnef311/nWt77F3XffzauvvhpQfJ/5zGfYvXu3//WnP/2J1NRULr30UgDeffddZs6cybRp03jttdd45pln+Pvf/86CBQv823j44Yd55ZVXePrpp9mwYQO7du3i9ddfD8bhExERERERHDyuaYx5OIDFrLV2eTfi6dUGDx7M97//fYwxjBs3jg8//JCf/OQnjBs3jq1bt7Jz505GjBgBwC233MKrr77K008/zeOPP+7fxgMPPMBll10GwKJFi7j66qs5cOAAw4cPZ9iwYdx5553+Zb/5zW/y5z//mT/+8Y984Qtf6DK+Pn36MHjwYABqamr46le/ymWXXcZ3vvMdAFatWsWMGTO44447/Ov853/+J1/4whcoLy8nKSmJ//7v/2bVqlVMmTIFgB//+MdMmjSpm0dORERCLZBejLtaRm2lRUR6hpM2eUs7KbOAafpXSd4puuiiizDG+KcvvvhiHn30UXbt2oW1lksuuaTV8rW1tf6ErtnEiRP9/x86dCgAHo+H4cOH09DQwH/913+xYcMGDh48SF1dHXV1df6auEBZa8nLy6OhoYE1a9b4Y37nnXfYu3cvGzdubLUswN69e+nfvz91dXVcfPHF/vKUlBQmTJjg6P1FRERERKRjTpK8MztYPxNYCKQC3whGUHIyYwxbt24lMTGx1fx+/fq1mm5Z3px8Nbfr+/GPf8zq1at57LHHmDBhAsnJySxfvhyPx+MolhUrVrB9+3a2bt1KcnKyf77P5+PrX/86t91220nrDB06lA8//NDR+4iISORZuXKl43UWLlwYgkhERKQjASd51tr9HRR9aIzZArwG3Aw8GIzAeqO//e1vWGv9ydlbb73F0KFDmTRpEtZaDh06dFLNnRPFxcVMmzaNWbNmAY21bB9++CGpqakBb+O5555j1apVPPfccwwfPrxV2eTJk3nvvfcYM2ZMu+uOHj2axMRE3nrrLUaPHg1AdXU1u3fv5swz27uHICIiIiIiTgWl4xXb+EzeOhoHRJdTVFpaygMPPMAHH3zgT6bmz5/PqFGjyM7OJi8vj+eee459+/bx9ttvs2rVKkdj8WVmZvLqq69SXFzMnj17uO+++9i/v6Pc/WTvvvsueXl5LF68mBEjRlBWVkZZWRlHjhwB4K677mLHjh3cfffd7Nq1i5KSEoqKivx3cFNSUvja177G0qVLefnll9m9ezcLFiygoaHB2YESEREREZEOOXlcsyt9gEFB3F6vM3PmTHw+H1deeSXGGL72ta9x++23s2fPHr773e/ywgsvsHTpUg4cOMDAgQP59Kc/7ahm79577+Wjjz7iK1/5Cv369eOrX/0qN9xwA++//35A67/zzjscP36cBx98kAcf/HeF7aWXXsqmTZuYOHEiBQUF5Ofnk5ubi8/n41Of+hTXXHONf9lHHnmE48ePM2fOHJKSkrjllls4fvx44AdJREREREQ6FZQkzxhzEXAXsDsY2wuFjIyMHh2gPCMjw/E6CQkJrFixghUrVrRbdv/993P//fe3u+7nPvc5Dh8+3GreqFGjWs07/fTT+dWvftVpDKtXr+5w+qabbuKmm27qdP0LLrig0+OcnJzMmjVr/ENDiIiIiIhIcDkZQqGkg6I04DSgHpgXjKBCwcljjSIiItI+daIiIhL5nNTkfUTjEAktWWAHsAf4mbV2X6AbM8a4gVzgkLV2UtO8pcAtQHnTYg9aa190EKOIiIiIiEiv5qR3zcuD/N5PAz8G2j4/uNJa+/jJi8c21TSKiEg00BAKIiKBy8rK6vYyxcXFjt83KL1rngpr7WvA4S4XFBERERERkYA5aZM36lTewFr7kcNV7jDGzAHeAu6x1h5pu8ChQ4eYO3cuCQkJNDQ0cN1115GXl0dpaSlJSUn079+fhoYGfzk0dlzi9XqJj48HoKGhgcTEROrr6wGIj4+nvr6e+Ph4rLX4fD4SExPxer0YYwIuT0hIwOfztSqPi4sjLi7OX97Q0IC1tlW5MaZVzC3Lm9XW1sbMPrVX3tU+HT9+nISEBMrLyxkwYAANDQ1UV1czZMgQSktLSUxMbDXm36FDh6ipqfGX9+nTh9NOO42KigoGDhxITU0NJ06c8Jf369ePpKQkjhw5wqBBgzh69Ch1dXX+8qSkJPr06UNlZSXp6elUVlbi9Xr95cnJycTHx1NVVUVGRgaHDx/GWktGRgZlZWWkpKQAcOzYMQYPHkx5eTnGGNLS0rrcJ4/HQ2pqKnV1dSHZp2YHDhyImX0K5efU8jsZK/sUys+p2f79+2Nmn0L5OTX75z//GVH7NGTIEDIzM/H5fBhjWo3r2vz/xhGVOKk8OzubN998MyT71Pb7GOvfp+7uU8vzK1b2KSMjgwkTJjBs2DAqKirIyMjg2LFjQOPwTeXl5QwaNAhrLUeOHCE9PZ2jR48SHx9PdnY2xcXF7N+/v919ann9ipRrRCR/Tm2/j7GwT939nK644goaGhrweDyMHz+eAwcO0KdPH9LT09m1axfnnXceNTU1fPzxx5x11ln885//JCUlhYEDB7JlyxaysrI4ePBgu/vUGdN8Qe6KMcbHyW3yumStje9km6OBghZt8gYDnqb3WQ4Mtda62q63fft2O378+Ha3WVVV1epHWCzYvbux09JzzjknzJGEV6CfbXOV96lUbfdWOmbO6Hg5o+PlTKQer+a4uvO4Zij2KVKPV6SKxePVvE8zZ850vG5zj+AdHY9YPF6hpON1suZjMmzYMMfrHjhwAOj4eO7YseNvU6ZMuai9MicdrzwCTAfOB7YA7zbNnwhMAd4ButWwzFpb1vx/Y8zPgYLubE9ERERERKS3cZLk7QHOBC601r7TssAY82ngJWCPtfa3pxqMMWaotfZg0+SXgX+c6rZERERERER6IydJ3v3Aj9smeADW2h3GmNXAA0BASZ4x5rfA5UC6MeZjYAlwuTHmfBof19wHfMtBfCIiIkHh8XhYvHgx+fn5DBo0KNzhRBT1lCkiEvmc9K55FnCok/KypmUCYq39qrV2qLU20Vo7wlr7S2vt162151prz7PWXtuiVk+66aOPPiItLY2333473KGIiEQ8t9vNzp07cbvd4Q5FRETEMSc1eQeB64wxq22b3lqMMXHA9UBpMIMLpmuuuYaKiooee79Bgwbxwgsv9Nj7iYhIcHg8HgoKCrDWUlBQgMvlUm0enXekoM4WRHqPcI37Js44SfJ+DuQDRcaYJ4D3m+aPB+4GLgMWBze84KmoqDilXm1OVXNvOL1NXV1dqy75RUSijdvt9g8F4PP5cLvdLFq0KMxRiYiIBM7J45rfB35EY0+aLwD/r+lVAFwBrLbWPhb0CHuJ6dOnc++997J8+XLGjh3LuHHjeOihh/D5fABcf/31rFq16qR17rvvPv90XV0dy5cv57zzzmPIkCFccMEF/PSnP+3wPd977z1mzZrFqFGjGDduHPPmzaOszN/BKTt27OC6665j7NixjBo1iquuuoq//vWvrbaRlpbGL37xC+bMmcOIESNYvnx5MA6HiEjYFBUV+cco9Xq9FBYWhjkiEZHIM2zYMMcv6TkBJ3m20beBc2jshOVnTa/7gQnW2jtDE2Lv8Yc//IH4+HgKCwv5wQ9+wNq1a9mwYUPA699+++08++yzfO9736O4uJgf/ehHrQYHb6m0tJTc3FzOOecctmzZwoYNG6iurmb27Nn+xPLYsWPMmjWLF154ga1bt3Luuecya9YsDh8+3GpbK1asYOrUqWzbto158+ad+gEQEQmjrKwssrKyqK6ubjW/urraXyYiIhINnDyuCYC1dg/wwxDE0uudffbZPPjggwCMHTuWX/3qV7z66qtMmDChy3U//PBD1q9fz+9//3umTp0KwOjRoztc3u12M2nSJJYuXeqft2bNGsaMGcPbb7/NhRdeyOc///lW6/zgBz9g06ZNbN26la985Sv++V/+8peZM2eOgz0VEREREZFQcZzkSehMnDix1fTQoUPxeDwBrbtr1y7i4uK47LLLAlp+586dvP7664wcOfKksn379nHhhRdSXl7Oo48+yrZt2zh06BA+n4+amho+/vjjVsuff/75Ab2niEgka9kRQHOt3fXXX6/2eCIiEnWU5EWQxMTEVtPGGP+jk3FxcbTp1NTfZuRU+Hw+srOzeeSRR04qy8jIABof/ywvLyc/P5+RI0fSt29fZsyYQV1dXavl+/fvf8pxiIhEMpfLFe4QREREHFOSFyVOP/30Vp2inDhxgg8++IDzzjsPgHPPPRefz8ef//xn/+OanZk8eTIbN25k5MiRJyWXzd544w0ee+wxsrOzATh06FCrGEREYp2GThARkWjkpHdNCaNPf/rTrFu3jm3btrF7924WLFhAfX29v6qrsoUAACAASURBVHzs2LHMmDGDu+66i+eff579+/ezfft2fve737W7vblz51JVVcXcuXN566232LdvH6+88grf/va3OXr0KACZmZn84Q9/4L333mPHjh3MmzdPwyOIiIiIiEQ4JXlRYs6cOVx22WXMnj2bmTNnkpWV5a/Fa7ZmzRpmzpzJAw88wGc/+1ny8vKoqqpqd3tDhw5l8+bNGGO44YYb+I//+A/uu+8++vbtS9++fQFYtWoVx44d44orrmDevHnMnj273TZ8IiIiIiISOXrN45qDBg3q0QHKnT7is2nTppPmrV69GoDdu3eTnJzML37xi1blc+fObTXdt29fli1bxrJly07a1qhRo04a+iAzM5Nnnnmmw5gmTZrE1q1bW82bNWtWq+m22xQRERERkfAKOMkzxgwCzrDW7m4x70zgbiAN+JW1tij4IQbHCy+8EO4QREREREREQs5JTd6TwDjgMwDGmBTgz0Dz8PWzjDFXWGtfC26IIiIiEmnaGxy+7byWw1KIiEjPcZLkXQL8usX0LBoTvKuBd4AtwH2AkjwRERERiRrt3bRwUq4bGhJpnCR5g4F/tpi+CnjLWlsIYIx5msZHN0VERCTG6UetiEjgerJvEHCW5HmBpBbTXwCebjH9CaABhUREREQk6gwbNqzrhdro6R/uIoFykuTtAa43xqwGptPY2cpLLcpHAhHR1aK1FmNMuMOQILLWYq0NdxgiIiIiIo719E0EJ0neahpr7o4A/YESWid5lwF/P+VIgqR///5UVlaSmpqqRC9GWGuprKwkJSUl3KGISAzpqo1NIMvokUUREYlEASd51tpfGWMsMAOoBB611nrBP7zC6cBPQhKlAwkJCaSkpFBVVRUzSd6uXbsAGD58eJgjCQ9rLSkpKSQk9JphHUVERE5Zd29g6OaFSPRz9KvZWvtrWvew2Ty/ArgwWEF1V0JCAqmpqeEOI2hWrlwJnDwQuYiIdF/zNdaJhQsXhiCS6OXxeFi8eDH5+fkMGqTm+SIi4dbtqhFjzIU0ts/7/+3dfXRcZ33g8e/Pb4mxY2NbWjsh2aT4kLquS0mhTdpSmuImdoIKrZO2sG0JVfqSxnCKu5sDLc5uofUuS7e4L5C63e0cSLcESozBFrFFSsumnMa0i8GhgbABJyaJLdljG78kdixLz/4xI0VW9DLXGunemfl+ztHR3BeNfveZO/PM7z7PfZ5/SimdmXxIkiSpkZRKJfbu3UupVOKuu+7KOxxVZb0HyEFEpOYxo9YdI+I/RcSOEes+BvwLsAv4WkQsrXN8kiSpwMrlMl1dXaSU6Orq4siRI3mHJEktL0tL3puBLw0uRMTrq+vuozLgykYqk6H/x3oGKEnSVLLr5eSUSqWh0Y8HBgZszZOkAsiS5F3F+fPi/SxwEPjllFKKiDbgjZjkSZI0ploGxZhovyINjNHd3U1fXx8AfX197Nq1yyRPknKWJcmbB5wetvx64O/TC5OXfR34rXoFJknSdHDglclZs2YNO3bsoK+vj9mzZ7N27dq8Q5KklpclyXsG+AGAiLgSWAl8cNj2RcDz9QtNkqTmNd0T406Vzs5Ourq6AJgxYwadnZ05RyRJqnngFWAH8FsR8SHgfioJ3WeHbV8FPFm/0CRJ9VIul7njjjscFEN119bWRkdHBxFBR0eHUyhIUgFkacl7H/BK4E4qCd47U0q9ABExF/g54K/rHqEkadKmcoj7Rr/HzK6Xk9fZ2cm+fftsxZOkgqi5JS+ldCyltBp4KbAgpfSXI3b5SWBTPYOTJE2eQ9xrqrW1tbFlyxZb8SSpIDJPhp5SOjHKutPA3rpEJEmqq1KpxMDAAAD9/f1TNsR9o91jNl7r4WCrY5FGsZSkIiniPcJ6QZZ78oiISyLiP0fEFyPi8Yj40er6tur6FVMTpiTpQnV3d3Pu3DkAzp07x65du3KOSJIkTaWaW/Iioh34IvBy4FvV33MBUkrliLiNSlfO35mCOCVp2jT6PWYjve51r2Pnzp1Dy9dff31+wUiSmkKj9d7I23Qfe5bumn8ILAOuBb4DHBqx/TPA6jrFJUmS1HQXXSRpOmRJ8jqAe1JKeyJitDur9wFvq0tUklQAzXKV8qGHHjpv+Qtf+AJ33313TtFIktQ68rr/O0uS10alm+ZYBoCLJxeOpKnglfDWtmbNGrZv3865c+eYNWsWa9euzTskKbNmuegiSdMhS5LXAywfZ/s1VLpx6gLU8iV8on38Ei5pNJ2dnXR1dQEwc+ZM5zKTJKnJZUnyHgBuj4g/B84O3xAR1wJvBf6kjrFJqjOvhLemtrY2Ojo62LZtGx0dHc5lJklSk8uS5L0XeCPwFWA7kIDbIuLXgXXAAeC/1z3CFuOXcElTobOzk3379tmKJ0lSC6g5yUsp9UTEdcCHgE4ggF+hkuw9APxWSunolEQpSZqUtrY2tmzZkncYkiRpGmRpySOl9BTwpohYAHwvlUTvWyZ3kiRJklQMmZK8QSmlE8C/1jkWSZIkKRfe/qJmMqPWHSPiFyPi3nG2fzQibq1PWJIkSZKkC5GlJe/twLfH2d4PvAO4f1IRSZKUg9GmqRm5zqlqpObl4HdqJlmSvO9j/ATuK8DPTC4cqTbOKyhNrJb3ie8DSZKaT5Ykbx6V1rqxJOCSWp8sIkpAB3AopbSqum4x8AngKuBJ4BdSSscyxChJ0gUx4ZWk2tmKWWxZkrwngNdSmUJhNK8FvpPh+T5Sfa7h9/m9G/h8Sun9EfHu6vK7MjynWoxdK6SxDU9aBlv1TGQkSWp+WZK8bcC7I+LBlNJfD98QEZ3AzwN/VOuTpZQeioirRqx+E3B99fFHgS9gkidJkpSZFzU1Fca7WOgFxeLIkuS9n0oS9lcRsQH4anX9DwIrgW8C/3WS8SxNKR2sPu4Blo6206FDh7j99tuZNWsW/f39rFu3jvXr19PT08O8efOYOXMmJ06coL29naNHj5JSor29nd7eXubPnw/AqVOnWLp0KYcPHyYiWLx4MYcPH2bBggX09/fz7LPPsmzZMnp6epg9ezYLFy6kXC6zcOFCzp49y+nTp4e2z5kzh0suuYQjR46waNEiTp8+zZkzZ4a2X3zxxcydO5djx46xZMkSTp48ydmzZ4e2z507l+XLl7N8+XJ6enq4/PLLmTt3Lo888givfOUrKZfLnD17lssuu4zHHnuMl7/85cycOZPHHnuM7//+7+cb3/gGAPv37y/UMc2ZM4fjx4/T1tbG8ePH6evrG9o+2dfp6quv5qKLLuLqq68eKqfTp0/z9NNP84pXvIKnnnqK+fPns2jRoqHtp06dore3l1WrVvHNb36T3t7eQh3TVL5ON954I48++iirVq1i/vz5Q2Vy7NgxTp06xRVXXMHjjz8+6rk3f/58rrrqKp577rlCHdNUvk4rV67k8ssv5zvf+Q4rVqygv7+fffv2sWLFCg4cOMCcOXNoa2sb9dz76le/yoIFC9i/f3+hjmnQmTNn6v46XXnllbS3t3PllVcOlcnx48c5duwYV111Fd/+9rdZunTpqOfeqlWreOSRRzhw4EDDvJ8m+zrdeOONfOUrX+HVr371eZ/lvb29lYpw6VIeffTRUc+9JUuWcOmll/L8888X6pim8nW65ppraG9v54knnjjvs3z58uU8+eSTLFq0iIULF4567u3evZsrrrjivPdjEY5pql+nyXrqqacKd0wTvU4rV67ksssu48iRI7S3t3Pq1CkA5s+fz+HDh1myZAkpJY4dO0ZbWxsnT54cej/u3r2b6667jlmzZtX0PWL4ubdq1Sp2797N/v37G+L9NF3n3uBnVDMc01S+ToMOHDhwQcc0nkgp1fymj4iFwH8DfhFYVF19DLgP2JhS+m7NT1Z5vquArmH35H03pfTSYduPpZQWjfy7hx9+OK1YsSLLvyq8wSsfk+l+2EpXTSyvbCyvbJqxvKby6mozltdUsryysbyyu9Aya9TyGjzeW2/NPpPX/fe/MKag59jk2ZKXzWTLa8+ePV9evXr1a0bblmky9JTSceDOiFgPDKbr5ZQlUxxfb0RcmlI6GBGXAofq9LySJEmS1BJqngx9uFRxuPozlOBVW/omYztwW/XxbcBnJvl8kiRJktRSak7yIuLzEbFsnO0/zgv36dXyfPcBDwPfGxFPR8TtVO77uyEiHgd+urosSZIkSapRlu6aPwbsjYhfTSk9MLgyIgLYCNwNHBzrj0dKKb1ljE2rM8QkSaqqZfLzifbxPgpJkhpflu6a1wJHgR0R8cGImB0RLwP+AXgv8FngVVMQoyRJkiSpRjW35KWUHomIH6Iygfk7qbS4XQbMA96eUrpnakKUJGWxefPmzH+zYcOGKYhEkiTlIevomqcj4g7gauDHgQS8wwRPkiRJkooh0+iaEbEc+Gcq9+f9LfAd4E8i4u7qvXmSJEmSpBzV3JIXEb8E3AMMAG9OKX2yOmVCico9eT8VEb+UUqp58BVJUv3Z9VKSpNaWpSXvb4CvA9eklD4JlcnRU0q3AHdSGZhlb/1DlKQXK5fL3HHHHRw5ciTvUCRJkgolyz15fwS8J6V0buSGlNKWiPgicF/dIpOkcZRKJfbu3UupVOKuu+7KO5xCceAVSdKFKpfLbNy4kU2bNrFkyZK8w9EFyjK65rsm2P5vEfHDkw9JksZXLpfp6uoipcSOHTvo7Oy0IiqAAwcO5B2CJGmSvIjaHDINvDKeiHgJlSkVJGlKlUolzp2rdCo4d+4cpVIp54gkSWp8wy+idnV1eUtEAxu3JS8izgJvTSl9vLp8CZVRNd+TUvraiN1/DrgXmDkVgUqavGZpadm1axcDAwMADAwMsHPnTq82DpNX18vLLst+na9ZzklNPc8VTTXPscpF1JQSUKlfbc2rv+uuu27Cdbt37570/5mou+Yszm/tmwN0AH8y6f8sSRdo6dKlPPHEE0PLy5Ytm5L/Y4WvqeT5Jalouru76evrA6Cvr49du3aZ5DWoTJOhS2pszdLS0tPTc97ywYPO3ALjX/kbvEpYj6uDUh6a5fNLxTTWZ2OrfXauWbOG7du3c+7cOWbNmsXatWvzDqnpTNe5ZJInqeHcdNNNbNu2jZQSEcHNN988Jf/HL5WaSp5fkoqms7OTT3/600Clu2ZnZ2fOEelCmeQVjBV4NpZX6xitDztASomtW7eydetWoHWutkqSJI2lbqNrSpIkSWpcwwdeSSk5enUDq6Ul7+aIGBzV4CVAAn4+Il41Yr9X1zWyFmX3nWwsr9YxsoVusGXvs5/9rHPkSZJUB7t27TovyXP06sZVS5L3H6o/w/3mGPumyYUjSdmY4EnS6LyoqayWLFnCc889N7Tc3t6eYzSajImSvJ+aligkSZIk5WrkhYGnn346p0g0WeMmeSml/zNdgUiSJGnynE5FFyoixl0eaaxB0aZicm9l48ArkiRJkrjhhhvOW16zZk1OkWiynEJBkiRJEuvXr6e7u5uBgQFmzJjBnXfeOe7+ttAVl0meJDWp0brR2IVGkjSWtrY21qxZw86dO7npppsc3KyBmeRJkiRJAiqteQcPHpywFU/FZpInSU3KVjpJreT+++/PO4Sm0NbWxpYtW/IOQ5NkkiepsMYatSvLPiY6kiSp1ZjkSZIkqWFNdDHPaSPUikzyJBXe5s2bM//Nhg0bpiASSZKk4nOePEmSJElqIrbkSZKUgwMHDuQdQkOxvKTpUS6X2bhxI5s2bXIKhQZmkiep8Ox6KUnS9CiVSuzdu5dSqcRdd92Vdzi6QCZ5kqS6sKWlNg4SkY3lJU2fcrlMV1cXKSW6urro7Oy0Na9BmeRJLaRRv4Q78IokSVOvVCqRUgJgYGDA1rwGZpJXMI36JVxS67KlRZKaQ3d3N319fQD09fWxa9cuk7wGZZIn5aSWib6hPl+MG/1LuK1ykiRNvTVr1rBjxw76+vqYPXs2a9euzTskXSCTvIIY78t10b+AS5IkqfF1dnbS1dUFwIwZM+js7Mw5Il0okzwpJyOTdpP5F/PihyRJ06etrY2Ojg62bdtGR0eHg640MJM8SZIkSUClNW/fvn224jU4kzxJkiRJQKU1b8uWLXmHoUmakXcAkiRJkqT6sSVPDc0pJzSVPL8kSVIjMsmT1BDGmnJi5HoHYpEkSa3OJE8NyVEXNZUafV5BSZLU2kzyJDUEEypJkqTaFDLJi4gngZNAP3AupfSafCOS6mOsLodZ9jHZkaTWUkvdYd0gabhCJnlVP5VSKucdhCRJkiQ1kiIneVLT2rx5c+a/2bBhwxREIkkqupGtdN4XLGkiRU3yEvC5iEjAX6aU/mr4xkOHDnH77bcza9Ys+vv7WbduHevXr6enp4d58+Yxc+ZMTpw4QXt7O0ePHiWlRHt7O729vcyfPx+AU6dOsXTpUg4fPkxEsHjxYg4fPsyCBQvo7+/n2WefZdmyZfT09DB79mwWLlxIuVxm4cKFnD17ltOnTw9tnzNnDpdccglHjhxh0aJFnD59mjNnzgxtv/jii5k7dy7Hjh1jyZIlnDx5krNnzw5tnzt3LnPmzOH48eO0tbVx/Phx+vr6hrYP2r9/f9Mc01S+TsPPk6Id08qVK7n88stJKZFSAiAiSCkREZWTv/p45Pbly5ezYMEC9u/fX/djGnTgwIFpe50a+dwb9PzzzzfNMU3l6zT8M6xZjmkqX6dBTz31VNMc01S+TiPfj81wTBO9TiOPOcsxDT+/inRMU/k6DXr66aczH9Pwz68iHVMzvk4eU/ZjGk8MfpEskoh4WUrpmYj4d8CDwDtSSg8Nbn/44YfTihUr8gtwmnnFLpsil9dgbJNpyZuK4ypymRWR5ZWN5ZWN5ZVNK5bXZI7Z8pq+v5Wm2p49e768evXqUccuKWRLXkrpmervQxGxDfgR4KHx/0pqHHa9lCRJ0lQpXJIXEfOAGSmlk9XHNwLvyzksSZIkNbmxRjIdud6WPRVd4ZI8YCmwrXp/0izgYymlXfmGJNWXA69IkqaDSYvUmgqX5KWU9gE/mHcc0lQyYZMkqXhMdtUsCpfkSZIkqT5MWqTWZJInTaPxKltH8JLUqsrlMhs3bmTTpk3nTREgSbowJnmSJClXpVKJvXv3UiqVuOuuu8bd13vMJGliJnmqm7Eq3pGseCVJg8rlMjt27CClxI4dO+js7Gz51rxa69Px9rOulVqbSZ4kScpNqVTi3LlzAPT19U3YmmfyIkkTM8lT3YyseL3HTJI0kZ07d5JSAiClxAMPPDBhl81Wceutt2b+m/vvv38KIpHUaEzyJElSbpYtW8YTTzwxtHzppZfmGE2xmLBdGO/blGBG3gFIkqTW1dvbe95yT09PTpFIUvOwJU+TUsvN4d4YLkkabrx64bnnnuO6666zfsDumhfKc0cyyZNyY3eSxlHLa+XrJEmSisIkT3WxefPmTPtv2LBhiiKRlDeT4mxa8YLPWAN13XLLLQ66Ikl1YJIn5aSZvrA1O18raXp0dnbmHUKh2PVS0oUyyVNd2DInaZBJcTaW1wtafRJ0SaoXkzxJkiaplkGoTOZUq4nOFeehlTQRkzzVhffkSZIkScVgkjeOWq7MglfSJKnVjTWQiPVDRa31qVPuSFJ9mORJLcjRDyVJkpqXSd44vDJbO7tfSpIm4uTekjQ9TPLUFEZrmWrmOaYmq9XLolwus3HjRjZt2uRoftI0MmGTpOlhkqdJGS9ZsOVTRVUqldi7dy+lUsmJl8fgaJGSJDUukzw1Bb9sqlblcpmuri5SSnR1ddHZ2WlrnjJzIJELY3dNSZoeJnnDWGkXh93pNFVKpRIpJQAGBgZszRvD8M8yW+VVLyZskjQ9TPJUN7WM2Ai1fVG0O53qbbTzs6+vj61bt7J161bAJEbZ2TIlSSoik7xpkrVlqpUHEimXy+zYsYOUEjt27LA7naTCMmGrzUT1la3F43PaG0lZmeRNk1ZomapXBVMqlTh37hxQaWlp5jKbDAfGyGawLMrlMh0dHQBcdNFFfOpTn/IigiRJaiomeaPYvHlz5r8Zb564CxnooZW/nO/cuXPonqmUEg888IBJnuqmra1t6HFHR4cJnqTCa+XvBJIujEneKOo9sbcDPWSzbNkynnjiiaHlSy+9NMdoisuBMSavs7Mz7xAKo5aW4Yn28fxTrex+KElTyyRvioxVgTnQw8R6e3vPW+7p6ckpEjU7W/F0obzHTJJUZCZ5w1hp52e8FoLnnnvOsldmtkxJxTX8vVUul3njG9/IwMAAM2fOZPv27V6AkaRJMsmbIiMrMAd6kCTpxT784Q8zMDAAQH9/P/fccw933313zlFJUmMzyZsGDvQwsZEtJoMtLLfccov3L1bZMnVhnMcsG8tL0+3BBx88b7m7u9skT5ImySRvHPWc3HuQAz1kY3lpskxAsrG8LowDiUiSisQkr45qaWl5wxveMO72VvoSYHlJkm644QZ27tw5tLxmzZoco5Gk5mCSN46RCcQHPvABtm3bxrp16+xCKKmpjHfBxIGPJjbyPux169Zx9uxZ78Ouwfr16+nu7mZgYIAZM2Zw55135h2SJDW8GXkH0ChGTmh+5MiRvEOSJBXQaHOjamxtbW1DrXc33XSTCbEk1YEteTVyQnOpsdgyNTmjdaeezP3IraS7u5u+vj6gMjfqrl27rC8msH79eg4ePGgrniTViUlejay0lTeTFqkxrFmzhh07dtDX18fs2bNZu3Zt3iEVXltbG1u2bMk7DElqGiZ5NbLSrj+TFk0nW6aysSwuXGdnJ11dXQDMmDHDUYIlSdPOJK9GtVTaJi2aTiYtUjG1tbXR0dHBtm3bnBtVkpQLk7waWWlPramYk1AaznNH06mzs5N9+/bZiidJyoVJXgZW2ioSkxapuLzHTJKUJ5O8DLJU2rZMZWM5SJIkSfXhPHmSJEmS1EQK2ZIXEWuBPwVmAv8rpfT+nEPKzJYpSZIkSXkoXEteRMwEPgzcBKwE3hIRK/ONSpIkSZIaQ+GSPOBHgG+llPallM4CHwfelHNMkiRJktQQithd82XAU8OWnwauHb7DoUOHuP3225k1axb9/f2sW7eO9evX09PTw7x585g5cyYnTpygvb2do0ePklKivb2d3t5e5s+fD8CpU6dYunQphw8fJiJYvHgxhw8fZsGCBfT39/Pss8+ybNkyenp6mD17NgsXLqRcLrNw4ULOnj3L6dOnh7bPmTOHSy65hCNHjrBo0SJOnz7NmTNnhrZffPHFzJ07l2PHjrFkyRJOnjzJ2bNnh7bPnTuXOXPmcPz4cdra2jh+/Dh9fX1D2z0mj8lj8pg8Jo/JY/KYPCaPyWPymIYf03gipTRx2jWNIuJWYG1K6deqy78CXJtSevvgPg8//HBasWJFXiFKkiRJUq727Nnz5dWrV79mtG1F7K75DHDFsOXLq+skSZIkSRMoYpL3r8ArIuJ7ImIO8GZge84xSZIkSVJDKNw9eSmlcxHxdqCbyhQKpZTSozmHJUmSJEkNoXBJHkBK6QHggbzjkCRJkqRGU8TumpIkSZKkC2SSJ0mSJElNxCRPkiRJkpqISZ4kSZIkNRGTvDr6yEc+kncIDccyy8byysbyysbyysbyysbyysbyysbyysbyyqYRy8skr47uvffevENoOJZZNpZXNpZXNpZXNpZXNpZXNpZXNpZXNpZXNo1YXiZ5kiRJktREIqWUdwyZff7znz8M7M87jpGOHj3atnjx4nLecTQSyywbyysbyysbyysbyysbyysbyysbyysbyyubApfXlatXr24fbUNDJnmSJEmSpNHZXVOSJEmSmohJniRJkiQ1EZM8SZIkSWoiJnl1EBFXRMQ/RsTXI+LRiPjtvGMqsoi4OCL+JSL2VsvrvXnHVHQR8WREfC0ivhoR/zfveIosIr63Wk6DPyci4p15x1U0EVGKiEMR8W/D1i2OiAcj4vHq70V5xlgkY5TX70fEM8POtZvzjLEoxqoTPb/GNla9GBHfExFfiohvRcQnImJO3rEWxWj1oufY6MaqF/0Me0GWOjEq/qz6vnwkIn4ov8jH5sArdRARlwKXppT2RMQlwJeBn00pfT3n0AopIgKYl1I6FRGzgS8Cv51S2p1zaIUVEU8Cr0kpFXFkp8KKiJnAM8C1KaXCjcibp4h4HXAKuDeltKq67gPA0ZTS+yPi3cCilNK78oyzKMYor98HTqWU/keesRXNWHUi8DY8v0Y1Vr0I/A7wqZTSxyNiC7A3pfQXecZaFKPVi36GTWx4vQj8Kn6GAdnqxGoy/A7gZirl+KcppWvzin0stuTVQUrpYEppT/XxSeAbwMvyjaq4UsWp6uLs6o9XGzQVVgPfNsF7sZTSQ8DREavfBHy0+vijVL6YizHLS6MYp070/BrDOPXi64H7q+sts4l5jk3MenEUGevEN1FJBlO1geKl1YtbhWKSV2cRcRVwDfClfCMptoiYGRFfBQ4BD6aULK/xJeBzEfHliPiNvINpIG8G7ss7iAayNKV0sPq4B1iaZzAN4u3V7jolu4a92Ig60fNrHCPrReDbwHdTSuequzyNF5CHG61e9Byb2Mh60c+wsY11Pr0MeGrYfoV8b5rk1VFEzAe2Au9MKZ3IO54iSyn1p5ReBVwO/EhErMo7poJ7bUrph4CbgPXVbgUaR/XelTcCn8w7lkaUKn35bWEf318Ay4FXAQeBP843nGIZr070/HqxkfUisCLnkIpu3HrRc+zFRqkX/QyrUSOeTyZ5dVLtQ78V+NuU0qfyjqdRpJS+C/wjsDbvWIospfRM9fchYBuVLwAa303AnpRSb96BNJDewS4n1d+Hco6n0FJKvdUv5gPA/8T35ZAx6kTPrxoMqxd/lEo3sFnVTZdTuZdKjFkveo6N77x60c+wCY11Pj0DXDFsv0K+N03y6qB6w/RfA99IKX0w73iKLiLaI+Kl1cdzgRuAx/KNqrgiYl518AIiSt8bdwAABHxJREFUYh5wI/Bv4/+VgLdgV82stgO3VR/fBnwmx1gKb8Q9GD+H70tg3DrR82sMY9SL36CS7N1a3c0yqxqnXvQcG9959aKfYRMa63zaDry1OsrmdcDxYd06C8PRNesgIl4L/BPwNWCguvr3UkoP5BdVcUXEK6ncwDqTyoWGv0spvS/fqIorIl5O5SolwCzgYymlTTmGVHjVSv87wMtTSsfzjqeIIuI+4HqgDegF/gvwaeDvgH8P7Ad+IaXkYCOMWV7XU+nmlIAngd8sYkU/3caqE6ncl+f5NYqx6sXq5//HgcXAV4BfTik9n1+kxTBWvRgRS/AcG9Vo9WJE/A1+hgHZ6sTqhawPUemF9hzwqymlwk1vZZInSZIkSU3E7pqSJEmS1ERM8iRJkiSpiZjkSZIkSVITMcmTJEmSpCZikidJkiRJTcQkT5KkaRIRb4uIFBHX5x2LJKl5meRJkppeRFxfTa6G/5yKiD0RsSEiZuUdoyRJ9WKlJklqJfcBDwABLAPeCnwQ+D7gN3KMS5KkujHJkyS1kj0ppf89uBAR9wCPAb8WEe9JKR3OLzRJkurD7pqSpJaVUnoW2E2lZW/54PqIuDEiPhER+yLidER8NyI+FxE/OfI5IuILEfFkRFwWEfdFxLGIeC4iuiPi6lriiIj3VLuQ/nlEWDdLkibFikSS1OoGk7ujw9a9DVgM3Au8A9hMpUvn5yPiJ0Z5jnnAQ0A/8HvAh4Drgc9ExMyx/nFEzIyIvwD+EPjdlNI7UkoDkzoaSVLLs7umJKmVvCQi2njhnrw7gGuAf0kp/b9h+/16tZVvSERsAR4Ffhf4pxHP2wb8UUrpA8P2Pwx8APhpoHtkIBExF/gY8AbgtpTSvZM8NkmSAJM8SVJreW/1Z7hPAeuHrxie4EXEfOAiKq10XwKuG+V5B4A/G7HuH6q/X8GLk7zFwIPADwI/k1J6URIoSdKFMsmTJLWSvwI+CcwGfgB4F3A5cGb4ThGxHNgErAFeOuI50ijPeyCldGbEuiPV30tG2f8jwHzgdSmlL2aIX5KkCXlPniSplTyeUvr7lNLOatfKnwF+GNgyuEO15e4hYC3wp8CtVJK9G6i0zsUoz9s/zv8cbf9PUGn9u7vabVOSpLoxyZMktayU0j8DfwP8YkT8WHX1auAyYENK6fdTSltTSp9LKf09lQFW6uFvgV8GXg90RcRL6vS8kiSZ5EmSWt4fUGmJe191ebBV7rwWuIi4Ebi2Xv80pfRx4C3ATwA7qy2IkiRNmvfkSZJaWkrpWxHxceCXqtMjfBHoAf44Iq4CngZeBfwK8DUq9/LV63/fHxF9wN8B3RFxU0rpRL2eX5LUmmzJkySpMsjKAPC+lNJ3qdyD9yUqc+T9MbASuBnYU+9/nFL6DLAOeDXwuYhYWO//IUlqLZHSaIOESZIkSZIakS15kiRJktRETPIkSZIkqYmY5EmSJElSEzHJkyRJkqQmYpInSZIkSU3EJE+SJEmSmohJniRJkiQ1EZM8SZIkSWoiJnmSJEmS1ET+P7gWQB2VRZzIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "plt.figure(figsize=(15,7))\n",
    "plt.style.use('bmh')\n",
    "ax = sns.boxplot(data=df_l,hue = 'level_1', x='rank',y= 0, palette='Greys')\n",
    "#ax.set_title('Excess quadratic risk for lambda = ' +str(title[0]) +' and n = ' + str(n), fontsize = fontsize)\n",
    "ax.set_xlabel('Rank', fontsize = fontsize)\n",
    "ax.set_ylabel('Excess quadratic risk', fontsize = fontsize)\n",
    "ax.legend(fontsize =fontsize-4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_l.to_csv('Save_n='+ str(n_list[i_n]) +'_M=' + str(M) + '_' + str(list_alpha[i_alpha]) + '_v3.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\begin{tabular}{lrrrrrrrrrrr}\n",
      "\\toprule\n",
      "{} & \\multicolumn{11}{l}{0} \\\\\n",
      "rank &   2   &   3   &   5   &   7   &   10  &   15  &   20  &   30  &   50  &   75  &   100 \\\\\n",
      "level\\_1   &       &       &       &       &       &       &       &       &       &       &       \\\\\n",
      "\\midrule\n",
      "full rank & 29.35 & 29.60 & 29.65 & 29.61 & 29.68 & 29.67 & 29.70 & 29.56 & 29.54 & 29.76 & 29.59 \\\\\n",
      "nuclear   &  9.12 &  8.46 &  9.04 &  8.95 &  9.49 &  9.44 &  9.36 & 10.82 & 11.24 & 12.64 & 13.60 \\\\\n",
      "oracle    &  2.08 &  3.15 &  4.79 &  6.60 &  8.94 & 12.29 & 15.29 & 19.91 & 25.96 & 29.30 & 29.59 \\\\\n",
      "penalized &  1.09 &  1.13 &  1.27 &  1.43 &  1.65 &  2.17 &  2.63 &  3.70 & 18.85 & 23.84 & 27.18 \\\\\n",
      "\\bottomrule\n",
      "\\end{tabular}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "latex_tab = df_l.groupby(['level_1','rank']).mean().unstack().to_latex(float_format=\"%.2f\")\n",
    "print(latex_tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the perfomance of the full rank regression doesn't seem to be affected by the rank of the underlying matrices. Low rank regression perfoms faster than the other regression and better when the underlying rank of matrix are small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] Birg\u0013e, L. and Massart, P. (2006). Minimal penalties for gaussian model selection.\n",
    "Probability Theory and Related Fields, 138(1-2):33{73.\n",
    "\n",
    "[2]  Negahban, Sahand and Wainwright, Martin J. (2011) Estimation of (near) low-rank matrices with noise and high-dimensional scaling. The Annals of Statistics\n",
    "\n",
    "[3] Ji, Shuiwang & Ye, Jieping. (2009). An accelerated gradient method for trace norm minimization. Proceedings of the 26th International Conference On Machine Learning, ICML 2009. 58. 10.1145/1553374.1553434. "
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}