L2_Optmization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "from scipy import stats\n",
    "from sklearn.datasets import make_regression\n",
    "from numpy import linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Basic optimization\n",
    "\n",
    "- Prior to 1950 many people used maximization rather than minimization\n",
    "\n",
    "`Model fitting:`\n",
    "$$F(d_i, w) =  y_i \\\\ y_i= \\text{label data} \\\\ \\\\ F= \\text{model} \\\\ y_i = \\text{label/outputs}$$\n",
    "\n",
    "- Optimization problems usually use x for the unknown rather than $d_i$\n",
    "\n",
    "\n",
    "`Linear model example:`\n",
    "\n",
    "$$d^T_iw = b_i$$\n",
    "\n",
    "`Loss function:`\n",
    "- Sample data vectors d and labels y then measure the loss function\n",
    "\n",
    "$$min \\sum_i \\ell(d_i, w, y_i$$\n",
    "\n",
    "`Least-squares loss function:`\n",
    "$$min ||D w- b||^2 \\\\ \\ell (d_i, w, b_i) = (d^T_i w-b_i)^2$$\n",
    "\n",
    "Where $d_i^T$ is rows of matrix d, then d*w is matrix vector product (inner product) minus $b_i$\n",
    "\n",
    "`Penalized regressions: ridge penality`\n",
    "$$min ||w||_2^2 + ||D w- b||^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Eigenvalue decomposition in optimization\n",
    "- Eigvalues values show the correct basis\n",
    "- Action of a matrix is described by eigenvalues\n",
    "\n",
    "\n",
    "$$A=UDU^{-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python example of decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8, 18,  4],\n",
       "       [ 4,  9, 13],\n",
       "       [17, 19,  5]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randint(20, size=(3,3))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 8. 18.  4.]\n",
      " [ 4.  9. 13.]\n",
      " [17. 19.  5.]]\n"
     ]
    }
   ],
   "source": [
    "values, vectors = np.linalg.eig(A)\n",
    "U = vectors\n",
    "U_inv = np.linalg.inv(vectors)\n",
    "D = np.diag(values)\n",
    "eigen = U@D@U_inv\n",
    "print(np.real(np.round(eigen,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example of eigenvalue decomposition in optimization:**\n",
    "$$x = \\beta_1 e_1 + \\beta_2 e_2 \\\\\n",
    "Ax =  A(\\beta_1 e_1 + \\beta_2 e_2) \\\\\n",
    "= A\\beta_1e_1 + A\\beta_2 e_2 $$\n",
    "\n",
    "Mulitply A by eigenvector = eigenvalue eigenvalue\n",
    "$$= \\beta_1 \\lambda_1 e_1 + \\beta_2 \\lambda_2 e_2 \\tag{2.1}$$\n",
    "\n",
    "**Inverse of 2.1**\n",
    "$\\beta_1 \\lambda_1^{-1} e_1 + \\beta_2 \\lambda_2^{-1} e_2 \\tag{2.2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1. Eigenvalue review\n",
    "\n",
    "Using slightly different notation:\n",
    "\n",
    "$$AV= c IV$$\n",
    "\n",
    "Where:\n",
    "- c = eigenvalue\n",
    "- IV = eigenvector\n",
    "\n",
    "Rearranged:\n",
    "\n",
    "$$AV- c IV = 0 \\\\\n",
    "(A-c I)V=0$$\n",
    "Where A-cI is called the characteristic matrix of A\n",
    "\n",
    "- With |A-c I|=0 there is infinite solutions.\n",
    "- Thus to find a uniqaue solution we must normalize by requiring the elements of $v_i$ of V so that $\\sum v_i^2=1$\n",
    "\n",
    "Matrix A is:\n",
    "\n",
    "1. Positive definte $\\Rightarrow$ eigenvalues>0 \n",
    "- Negative definite $\\Rightarrow$ eigenvalues<0\n",
    "- Postive semi-definite $\\Rightarrow$ all eigenvalues are nonnegative and at least one=0 \n",
    "- Negative semi-definite $\\Rightarrow$  all eigenvalues are nonpositive and at least one=0 \n",
    "- Indefinite $\\Rightarrow$ some eigenvalues>0 and some <0\n",
    "\n",
    "### 2.1.1 Finding eigenvalues:\n",
    "\n",
    "$$Det|A-cI| =0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.1.1 Python example: Eigenvalue with numpy\n",
    "\n",
    "- Using np.linalg.eig(MATRIX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3., -9.])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[-6,3],[3,-6]])\n",
    "np.linalg.eig(A)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.1.2 Python example: Eigenvalue with user defined function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eigen(matrix):\n",
    "    A = np.diag(matrix) #Diagonal\n",
    "    I = np.diag(np.fliplr(matrix)) #Opposite Diagonal\n",
    "    #Quadratic equation\n",
    "    c = np.prod(A)-np.prod(I) #C-term \n",
    "    b = -np.sum(A) #B-term\n",
    "    root_term = math.sqrt(b**2-4*(c))#Quadrtic formula root term\n",
    "    solu1 = (-(b)+root_term)/2 #Quadratic formula solution1\n",
    "    solu2 = ((-b)-root_term)/2 #Quadratic formula solution2\n",
    "    vec1 = matrix-np.diag(np.repeat(solu1,2))\n",
    "    vec2 = matrix-np.diag(np.repeat(solu2,2))\n",
    "    print(\"Original matrix: \\n\",matrix)\n",
    "    print(\"Eigen-values: {}, {}\\n\".format(round(solu1,3),round(solu2,3)))\n",
    "    #Classification of matrix\n",
    "    if solu1>0 and solu2>0:\n",
    "        print('Pos definite')\n",
    "    if solu1<0 and solu2<0:\n",
    "        print('Neg definite')\n",
    "    if (solu1==0 or solu2==0) and (solu1>0 and solu2>0):\n",
    "        print('Pos semi-def')\n",
    "    if (solu1==0 or solu2==0) and (solu1<0 and solu2<0):\n",
    "        print('Pos semi-def')\n",
    "    if (solu1<0 and solu2>0) or (solu1>0 and solu2<0):\n",
    "        print('Indefinite')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original matrix: \n",
      " [[0 0]\n",
      " [0 4]]\n",
      "Eigen-values: 4.0, 0.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = np.random.randint(5, size=(2,2))\n",
    "eigen(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Estimation problem\n",
    "\n",
    "Suppose:\n",
    "\n",
    "$$\\lambda_1 =1, \\ \\lambda_2=0.1, \\ \\lambda_3 = 0.01 \\\\\n",
    "Ax = b= \\beta_1 e_1 +0.1\\beta_2 e_2 +0.01 \\beta_3 e_3 \\\\\n",
    "\\hat{b} = b + \\eta$$\n",
    "\n",
    "$\\eta$, noise, would show up at every frequency uniformly\n",
    "\n",
    "Example of adding random Guassian 2-d noise with no directionality at 0.02\n",
    "\n",
    "$$AX = \\beta_1 1.02 +\\beta_2 0.12+0.03 \\beta_3 $$\n",
    "Inversion would blow up the matrix\n",
    "\n",
    "**This issue crops up all the time**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Condition number\n",
    "\n",
    "`Condition number:` Ratio of smallest to largest singular value\n",
    "\n",
    "$$\\kappa = \\frac{\\sigma_{\\text{max}}}{\\sigma_{\\text{min}}} \\\\\n",
    "\\kappa = \\frac{\\lambda_{\\text{max}}}{\\lambda_{\\text{min}}}\n",
    "\\\\ \\kappa = ||A||||A||^{-1}\n",
    "$$\n",
    "- A measure of how good you matrix is\n",
    "- Singular values are always positive\n",
    "- $\\lambda$ is for the symmetric matrix\n",
    "- The condition number tells you how much you matrix will react to additional noise\n",
    "- This a higher condition number will negatively impact convolutions, linear operators, etc\n",
    "- The condition number for a convolution is usually very large\n",
    "    - This is a problem for images\n",
    "    - Radar and sonar systems are designed to have low condition numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Python examples: Condition number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.1 Python example: Condition number blows up the noise\n",
    "Code written by: CMSC 764 hmwk2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. Recovery error with no noise = 1.63e-06\n",
      "2. Recovery error with noise of 1.79e-07 = 178.930\n"
     ]
    }
   ],
   "source": [
    "def condition_blowup(n):\n",
    "    A = sp.linalg.hilbert(n)# construct an n by n matrix\n",
    "    x = np.random.randn(n,1)# construct n by 1 signal \n",
    "    b = A@x \n",
    "    x_recovered = la.inv(A)@b #Can also use la.solve(A,b) \n",
    "    # 1. Recovering X: inverting A dot b \n",
    "    print('1. Recovery error with no noise = {:0.3}'.format(la.norm(x_recovered-x)))\n",
    "    \n",
    "    # 2.1 Now add small amount of noise to vector b\n",
    "    b_noise = b+np.random.randn(n,1)*0.0000001\n",
    "    # 2.2 The error blows up when we try to recover X\n",
    "    x_recovered_withnoise = la.inv(A)@b_noise\n",
    "    print('2. Recovery error with noise of {:.03} = {:.03f}'.format(la.norm(b_noise-b), la.norm(x_recovered_withnoise-x)))\n",
    "    \n",
    "np.random.seed(0)\n",
    "condition_blowup(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.2 Python example: Function that returns Mxn matrix with specific condition number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The 2x2 matrix A with condition number=2:\n",
      "[[8 8]\n",
      " [3 7]]\n",
      "Constructed matrix which will yield a condtion number=[[1.46829282 0.97228825]\n",
      " [0.02360495 1.37775707]]\n",
      "inf \n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def buildmat(m,n,condition_number):\n",
    "    A = np.random.randint(10, size=(m,n))\n",
    "    print(\"The {}x{} matrix A with condition number={}:\\n{}\".format(m,n, condition_number,A))\n",
    "    #print(\"Np.linalg.cond\",np.linalg.cond(A))\n",
    "    U, S, V = np.linalg.svd(A)\n",
    "    S = np.array([[S[j] if i==j else 0 for j in range(n)] for i in range(m)]) #Creates diagnoal matrix with S entries\n",
    "    S[S!=0]= np.linspace(condition_number,1,min(m,n)) #Creates matrix with condition number in 1,1\n",
    "    A = U@S@V\n",
    "    print(\"Constructed matrix which will yield a condtion number={}\".format(A, condition_number))\n",
    "m,n,condNumber = 2,2,2\n",
    "\n",
    "np.random.seed(100)\n",
    "buildmat(m,n,condNumber)\n",
    "print(np.linalg.cond(A),\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Deblurring\n",
    "\n",
    "- Blur with linear operator then invert to get back\n",
    "- Now blur with 2-d convolution\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.2 Python example: Deblurring\n",
    "\n",
    "- Blur with linear operator then take inverse to recover the image\n",
    "- If there is noise then you will be unable to recover the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 4. Under-determined systems\n",
    "\n",
    "What if a matrix the matrix is not full-rank? \n",
    "\n",
    "$$A \\in R^{MxN}, \\ \\ M<N \\\\\n",
    "b= Ax + \\eta$$\n",
    "\n",
    "\n",
    "\n",
    "If the error is bounded $||\\eta||<\\varepsilon$ solve:\n",
    "\n",
    "$$\\text{Min} ||x|| s.j.t. \\ \\ ||Ax-b||\\leq \\varepsilon$$\n",
    "\n",
    "- If noise is bounded then finding the simplist solution that satisfies the following constraint\n",
    "- Finds the shortest vector min 2-norm subject to constraint\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 Geometric interpretation\n",
    "\n",
    "- If two norm is too small there is no solution\n",
    "- 2-norm ball is tangent to the plane defined by the linear system\n",
    "- Picking different balls for constrain changes where the ball intersects line\n",
    "\n",
    "<img src='Images/l2_4.1_geometric.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Ridge regression\n",
    "\n",
    "- Ridge regression solves linear system with 2-norm penalty\n",
    "\n",
    "If the error is bounded $(||\\eta||\\leq \\varepsilon)$ \n",
    "\n",
    "$$b= Ax+\\eta$$\n",
    "\n",
    "$$Min ||x|| s.j.t ||Ax-b||\\leq \\varepsilon$$\n",
    "\n",
    "This is equivalent to:\n",
    "\n",
    "$$Min \\ \\ \\lambda||x||^2 +||Ax-b||^2$$\n",
    "\n",
    "- As $\\lambda \\uparrow error \\uparrow$\n",
    "- If $\\lambda=0$ then it is linear\n",
    "\n",
    "One way to analyze ridge penalty is to look at what at it does to the condition number of the system\n",
    "- The solution is a lot less noise sensative\n",
    "\n",
    "The new condition number:\n",
    "\n",
    "$$\\frac{\\sigma^2_{\\text{max}}+\\lambda}{\\sigma^2_{\\text{min}}+\\lambda}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1 Thihonov regularization\n",
    "\n",
    "- Anything that involves a smooth l-2 penalty\n",
    "- Easy to solve these linear systems\n",
    "    - Convex\n",
    "    - Smooth\n",
    "\n",
    "- Improved condition number, less noise sensitivity\n",
    "- Parameter $\\lambda$ can be set in two ways:\n",
    "    1. Empirically\n",
    "    2. Use noise bounds and theory (BIC, etc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Bayesian model for ridge regression\n",
    "\n",
    "### 5.2.1 Bayes rule review\n",
    "\n",
    "$$P(H|E) = P(\\text{Hypothesis}) \\cdot \\frac{P(\\text{Evidence}|\\text{Hypothesis})}{\\text{Evidence}} $$\n",
    "\n",
    "\n",
    "#### 5.2.1.1 Example: Baye's theorem\n",
    "\n",
    "#### The Cookie problem:\n",
    "\n",
    "- Bowl 1: 30 vanilla, 10 chocolate\n",
    "- Bowl 2: 20 vanilla, 20 chocolate\n",
    "\n",
    "**Suppose we randomly draw a vanilla cookie. What is the probability it came from bowl 1?**\n",
    "\n",
    "- This answer would be straightforward IF the question was P(Vanilla|Bowl 1).\n",
    "   - 3/4 because 30/40 cookies are vanilla in bowl one\n",
    "   \n",
    "Instead the answer can be found using Bayes\n",
    "\n",
    "$$P(\\text{Bowl 1}|\\text{Vanilla})= \\frac{p(B_1)p(V|B_1)}{p(V)}\\\\\n",
    "= \\frac{1/2\\cdot3/4}{5/8}\\\\\n",
    "=.6$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2.2 Bayes in ridge\n",
    "\n",
    "1. Choose distribution of data given parameters\n",
    "2. Observe data from random process\n",
    "3. Use Bayes rules\n",
    "4. Find most likely parameters given data\n",
    "5. Uncertainty qunatification and confidence\n",
    "\n",
    "\n",
    "$$Min \\ \\ \\lambda||x||^2 +||Ax-b||^2$$\n",
    "\n",
    "Assumptions:\n",
    "- Prior: Expected signal power: $E\\{x^2_i\\}=E^2$\n",
    "- Linear measurement model: $b=Ax+\\eta$\n",
    "- Noise is i.i.d Gaussian: $eta = . N(0,\\sigma)$\n",
    "- MAP (maximum a-posteriori) estimate: $P(x|b)$\n",
    "\n",
    "\n",
    "**Bayes' rule**\n",
    "$$P(X|Y) \\propto P(X|Y)\\times P(X)$$\n",
    "\n",
    "$$Max x P(x|b) \\propto p(b|x)p(x)$$\n",
    "- Unknown parameter is x\n",
    "- Model: $b \\sim N(Ax, \\sigma)$\n",
    "- Prior: $x \\sim N(0, E)$\n",
    "\n",
    "\n",
    "Probability of data:\n",
    "$$P(b|X) = \\prod_i \\frac{1}{\\sigma \\sqrt{2\\pi}}e ^{-\\frac{1}{2\\sigma^2}(b_i-(Ax)_i)^2} \\\\ \n",
    "= (2\\pi\\sigma^2)^{{-m/2}} e ^{-\\frac{1}{2\\sigma^2}||b_i-Ax||^2}\n",
    "$$\n",
    "\n",
    "Prior:\n",
    "$$P(x) = (2\\pi E^2)^{-n/2} e^{-1\\frac{1}{2E^2}}||x||^2$$\n",
    "\n",
    "\n",
    "Max: \n",
    "\n",
    "$$Exp(\\frac{-||b-Ax||^2}{2\\sigma^2})exp(\\frac{-||x||^2}{2E^2})$$\n",
    "\n",
    "- Ignore all the constants, we are trying to find the x that maximizes\n",
    "\n",
    "Negative Log-Likelihood\n",
    "\n",
    "$$-\\frac{||b-Ax||^2}{2\\sigma^2}-\\frac{||x||^2}{2E^2}$$\n",
    "\n",
    "Multiply by -1 to make into minimization problem\n",
    "\n",
    "$$\\frac{\\sigma^2}{E^2}||x||^2+||b-Ax||^2$$\n",
    "\n",
    "- $\\lambda = \\frac{\\sigma^2}{E^2}$\n",
    "    - $\\sigma$ = noise\n",
    "    - E = expected value of signal\n",
    "    - If we have more noise then the larger the $\\lambda$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2.3 Priors\n",
    "\n",
    "- Priors can add information\n",
    "- Ridge/Tikhonov priors require many assumptions\n",
    "- One general prior is sparsity\n",
    "- Usually we want weaker priors so our model doesn't break if we are wrong\n",
    "- Spare priors are often used\n",
    "\n",
    "### 5.2.3.1 Sparse priors\n",
    "\n",
    "- Signal has very few non-zeros: small $\\ell_0$ norm\n",
    "- Low density signals have a rapid decay of coefficients\n",
    "- Fast decay: small and weak $\\ell_p$ norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Overfitting\n",
    "\n",
    "- Sparsity used a lot in machine learning\n",
    "- Try to do signal estimation on high-dimensional vectors without over-fitting\n",
    "\n",
    "\n",
    "$$Ax = b$$\n",
    "- Features/data=A\n",
    "- Model parameters= x\n",
    "- Labels = b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1 Python example: Overfitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEKCAYAAAD3tSVSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfXwU5bnw8d9FSCBIBRWs8lZ4LFIVJEiwtuBDD8pBjxgp1pf2qcVqq9bTo9Yjiq0V9GiLpa3W1lax2loPHokvYNBjrcKRKlolAQygpYrYQwKWIAZfiBCS6/ljJmGzmX2f2Z3dvb6fTz6ws7Mz92Th2nuv+57rFlXFGGNMYeqR6wYYY4wJjgV5Y4wpYBbkjTGmgFmQN8aYAmZB3hhjCpgFeWOMKWAW5I0xpoBZkDfGmAJmQb6Iicg7InJqAMcdJSJrReRDEbnC7+Nn0hYR2SgiX0rh9Sntb0zYWJAPiBtAW0Tko4ifX+W6XVlyLfC8qn5KVe8MU1tU9ThVfR68P+Sit0Xub0w+siAfrDNVtW/Ez3dz3aAs+QywMZcNEJGeYWlLoRCReSIyr9DOVegsyGeZiBwlIrtE5AT38SAR2dmREhCROSKy2U0vvC4iX456/TsiMltE6kXkYxG5T0Q+LSJPu695TkQOidj3evc474vI70Skd4x2DRKRx0SkSUS2xEuziMgxIvK8iDS76YyqiOdWAP8E/Mr99nJ0sq93r/3RqH1/ISJ3JtNG93qvE5F64GOvtnT01EXkQWAYsMx97toY27r07N3H17i//90isrjjdyoiJ0Skhh5xn7sl1u8xqu19RaRNRI6M2DZaRLaLyKeSOUacYx8uIjUi8g8R+UBElonIwZkcMwwK9bp8p6r2E8AP8A5waoznvg28AfQBngF+GvHcOcAgnA/g84CPgSOjjvsX4NPAYGAHsAYYB/QCVgBzI/bdAAwFDgVWAbdEt9E9Vx1wI1AG/B/gbWCaR9tLgbeA77v7TgE+BEZF7PM88K0Y1x7z9Ti97j3Awe6+JcB24KRk2uhezzr3esu92hL5vni9R9HbYjx+1X2PDnXfx8vcNv0duNK9xpnAvqjf96+BX8f5N7MROCPi8ZPAv0Xt8yTQHOPnyRjH/Sww1f33cSjwMjA7jX/T84B5Wfr/k/Bcfl1Xof9YTz5YS93easfPtwFU9V7gTeAV4EjgBx0vUNVHVHWbqrar6mJ3vxOjjvtLVf2HqjYCLwCvqOpaVd0LLMEJ+B1+papbVXUXcCvwVY92TgAGqurNqrpPVd8G7gXO99j3JKAvMN/ddwVO4PE6rpeYr1fVv+N8YM1w950C7FHVv6TQxjvd621Jsj3puNN9j3YBy4AK97p6us+1qurjOB8GnVT1clW9PM5xVwMd3/D+L3AscE/UMaarav8YP9O9Dqqqb6nqs6q6123zs8AhItJPRF51v7WMTu9X0Z2IVInIjGQfpyvWdbnn+IKIvCwiK0Xkv0Sk1N1+q4i8ICKPikifTNuQD3om3sVkYIaqPhfjuXuBGuASNzgDICLfAK4Ghrub+gIDol77j4i/t3g87hvxeGvE3/+O0wON9hlgkIg0R2wrwfkAiTYI2Kqq7VHHHeyxr5dEr38I5wPjD8DX3MeptHErwXs34u97cK5pENCobhczzbasxvlmBfAT4Iequi/tVrpE5BzgKmAkzjeOPsAlOG0/A1gQ57VPApPchx1pqavcxy/G+GCpwoktS5N8nNa54lwXOP+mpqhqi4jcCpwlIn8FjlLVk0XkMuAioOAnQ1iQzwER6QvcAdwHzBORx1R1l4h8Bif4nwK8rKptIrIOkAxONzTi78OAbR77bAW2qOrIJI63DRgqIj0iAvUw4G9JtifR6x8BfiYiQ4AvA19IsY2pLJDgtW+6CyxsBwaLiEQE+qHA5hSOsRq4VkTOBsqB/4reQUSeBk6O8foXVPX0qP2nALfhpP7WupvfAdapaivQJBL7n1dkYBV3IFRV58W7CFX9ViqP0zlXvOtyXxf573w/0I7ze3va3fa0+/qCD/KWrsmNXwB17j/2p4C73e0H4QSZJgAR+SaQ6dfofxWRISJyKE4efLHHPq8CH7iDluUiUuIO+k3w2PcVnHGCa0WkVJwB4zOBh5NsT9zXq2oTTh79dzhB/Y002pisf+Dk9hNtS8bLQBvwXRHpKSJn0T3NlshrwBHAz4A5Ud92AFDV07XrjK3In9O7HRHG4nxA/hUnlXE/cDjweoptC5ukrktERgCn46QEDwF2u0/txsnjFzwL8sHqmKXR8bPE/c9/Gs5gHTipmRNE5P+p6us4/8Ffxgk2Y3AGSzPxEPAnnEHKt4Fusz1UtQ0n0FYAW4CdwG+Bfh777sP5un26u9+vgW+o6l+TaUySr38IJ23xUMTrkm5jCn4M3OCOl1wTZ1tC7nXNBC7GGQT9Ok5giUzF3S0id3sfAdy03XrgHVV9OtZ+KVqEMxD8rtueN4HX/UgD5VjC63Jn2jwAXOBuf58D/176Abuy2uIcka4pRFNIROQdnJklscYFTIBE5BXgblX9XZL7l+HMPDrXHWzOChH5Pc4Mrw3ZOmfQxLlP4gngZ+7gPiIyBrheVb8mIpcAvVT1l7lsZzZYT94Yn4jIZBE5wk3XzAKOB/6YwiHmAquyHOD/G/hn4F4RuTBb582CrwKfB24U556M81R1PfB3EXkBmIaT4il4NvBqjH9GAdU4s5s2A19R1e2JXiTOjXH/A9TjDDZnjar+SzbPly2q+iDwoMf263PQnJzyLV0jIiVALc40sunugMfDOIMbaziQFzPGGJMlfqZrrsS5+6/DbcDt7pS393EGpIwxxmSRLz15d07zAzh3VF6NMwuiCThCVfeLyBdwblGeFu84AwYM0OHDh2fcHmOMKSZ1dXU7VXWg13N+5eTvwCnp2lFI6TCgWVX3u48biHFHpDvKfQnAsGHDqK2t9alJxhhTHETk77GeyzhdIyLTgR2qWhe52WNXz68MqrpQVStVtXLgQM8PImOMMWnyoyc/EagSkX/BqTdxME7Pvr+I9HR780Pwvp3eGGNMgDLuyavq9ao6RFWH41QEXKGq/w9nSthX3N1m4dyYYIwxJouCnCd/HfCwOIsmrMUpxpWy1tZWGhoa+OSTT3xtnElf7969GTJkCKWlpbluijEmAV+DvDprYT7v/v1tUi/Q1E1DQwOf+tSnGD58OPGq5ZnsUFXee+89GhoaGDFiRK6bY4xJIPR3vH7yyScW4ENERDjssMNoamrKdVNMBpaubWTBM5vY1tzCoP7lzJ42ihnjkl0SwOST0Ad5wAJ8yNj7kd+Wrm3k+sfX09LaBkBjcwvXP74ewAJ9AbICZcYUmQXPbOoM8B1aWttY8MymHLXIBMmCfJIaGho466yzGDlyJEcddRRXXnkl+/bt4/e//z3f/e53c908li5dyuuvH1gv4cYbb+S556zCsOluW7P38rextpv8ZkE+CarKzJkzmTFjBm+++SZ/+9vf+Oijj/jBD36Q+MVp2L9/f+KdokQH+ZtvvplTTz01zitMsRrUvzyl7Sa/FVyQX7q2kYnzVzBizlNMnL+CpWsbMz7mihUr6N27N9/85jcBKCkp4fbbb+f+++9nz549bN26ldNOO41Ro0Zx0003AfDxxx9zxhlnMHbsWEaPHs3ixc6qe3V1dUyePJnx48czbdo0tm93KtF+6Utf4vvf/z6TJ0/m1ltvZfjw4bS3O6u/7dmzh6FDh9La2sq9997LhAkTGDt2LGeffTZ79uzhpZdeoqamhtmzZ1NRUcHmzZu58MILefTRRwFYvnw548aNY8yYMVx00UXs3essVjR8+HDmzp3LCSecwJgxY/jrX53FmVauXElFRQUVFRWMGzeODz/8MOPfoQmP2dNGUV5a0mVbeWkJs6eNylGLTJAKKsh3DCg1NregHBhQyjTQb9y4kfHjx3fZdvDBBzNs2DD279/Pq6++yqJFi1i3bh2PPPIItbW1/PGPf2TQoEG89tprbNiwgdNOO43W1lb+7d/+jUcffZS6ujouuuiiLt8GmpubWblyJXPnzmXs2LGsXLkSgGXLljFt2jRKS0uZOXMmq1ev5rXXXuOYY47hvvvu44tf/CJVVVUsWLCAdevWcdRRR3Ue85NPPuHCCy9k8eLFrF+/nv379/Ob3/ym8/kBAwawZs0avvOd7/DTn/4UgJ/+9KfcddddrFu3jhdeeIHycuvhFZIZ4wbz45ljGNy/HAEG9y/nxzPH2KBrgSqoIB/UgJKqes4o6dg+depUDjvsMMrLy5k5cyYvvvgiY8aM4bnnnuO6667jhRdeoF+/fmzatIkNGzYwdepUKioquOWWW2hoaOg83nnnndfl7x29/4cffrjzuQ0bNnDyySczZswYFi1axMaNG+O2fdOmTYwYMYKjjz4agFmzZvHnP/+58/mZM2cCMH78eN555x0AJk6cyNVXX82dd95Jc3MzPXvmxSQsk4IZ4wazas4Utsw/g1VzpliAL2AFFeSDGlA67rjjulXH/OCDD9i6dSslJSXdPgBEhKOPPpq6ujrGjBnD9ddfz80334yqctxxx7Fu3TrWrVvH+vXr+dOf/tT5uoMOOqjz71VVVTz99NPs2rWLuro6pkyZAsCFF17Ir371K9avX8/cuXMT3gmcqJR0r169ACcF1TEWMGfOHH7729/S0tLCSSed1JnGMcbkn4IK8kENKJ1yyins2bOHP/zhDwC0tbXx7//+71x44YX06dOHZ599ll27dtHS0sLSpUuZOHEi27Zto0+fPnz961/nmmuuYc2aNYwaNYqmpiZefvllwCnZEKsn3rdvX0488USuvPJKpk+fTkmJk0P98MMPOfLII2ltbWXRokWd+3/qU5/yzJ1/7nOf45133uGtt94C4MEHH2Ty5Mlxr3fz5s2MGTOG6667jsrKSgvyxuSxggryQQ0oiQhLlizhkUceYeTIkRx99NH07t2bH/3oRwBMmjSJCy64gIqKCs4++2wqKytZv349J554IhUVFdx6663ccMMNlJWV8eijj3LdddcxduxYKioqeOmll2Ke97zzzuM///M/u6Rx/uM//oPPf/7zTJ06lc997nOd288//3wWLFjAuHHj2Lx5c+f23r1787vf/Y5zzjmHMWPG0KNHDy677LK413vHHXcwevRoxo4dS3l5Oaeffnq6vzpjTI75tsarHyorKzU6LfLGG29wzDHHJH0Mu107O1J9X0x22f+D4iIidapa6fVcwY2ozRg32P4xm6JmZQtMpIJK1xhjrGyB6cqCvDEFxsoWmEgW5I0pMFa2wESyIG9MgbGyBSZSwQ28GlPsOgZXbXaNAQvySSkpKWHMmDGdj5cuXcrOnTv5wx/+wJ133snzzz9PWVkZX/ziFzufP/roozn22GNTOk/fvn356KOPfG27KU42y8x0sCCfhPLyctatW9dl2/Dhw6msdKalPv/88/Tt27dLkJ8+fXrKQd4YU9hW19zD0DULOFyb2CED2XrCbCZUXRroOQsvJ19fDbePhnn9nT/rqwM5zfPPP8/06dN55513uPvuu7n99tupqKhg5cqV3cr+bt68mdNOO43x48dz8sknd5YJ2LJlC1/4wheYMGECP/zhDwNppzEmHFbX3MPouhs4giZ6CBxBE6PrbmB1zT2BnrewevL11bDsCmh1p4rt3uo8Bjj+3LQP29LSQkVFBQAjRoxgyZIlnc8NHz6cyy67jL59+3LNNdcATnGx6dOn85WvfAVwat/cfffdjBw5kldeeYXLL7+cFStWcOWVV/Kd73yHb3zjG9x1111pt88YE35D1yygXPZ12VYu+xi6ZgEE2JsvrCC//OYDAb5Da4uzPYMg75WuSdZHH33ESy+9xDnnnNO5rWPRjlWrVvHYY48BcMEFF3Ddddel3UZjTLgdrk3QvWI5h+vOQM9bWEF+d0Nq27Ogvb2d/v37x/yQ8KpTb0xRqq92OmS7G6DfEDjlxow6Z2GzQwZyBE0e2wdwRIDnLaycfL8hqW33SXSZ38jHBx98MCNGjOCRRx4BnPrur732GuAszvHwww8DdCkbbEzR6Ui17t4K6IFUa0Bjarmw9YTZtGhZl20tWsbWE2YHet7CCvKn3AilUXf1lZY72wN05plnsmTJEioqKnjhhRe6lf1dtGgR9913H2PHjuW4447jiSeeAOAXv/gFd911FxMmTGD37t2BttGYUIuXai0QE6ouZcP4W3iXgbSr8C4D2TD+lsBn12RcalhEegN/BnrhpH8eVdW5IjICeBg4FFgDXKCq+2IfyZ9Sw4X+lS8srNSw8dW8/oBXLBKY15zt1uSdoEsN7wWmqOpHIlIKvCgiTwNXA7er6sMicjdwMfCbeAfyxfHnWlA3Jt/0G+Kmajy2m4xknK5RR8dtmqXujwJTgEfd7Q8AMzI9lzGmQOUo1VoMfMnJi0iJiKwDdgDPApuBZlXd7+7SAKR9j3WYVq8y9n6YABx/Lpx5J/QbCojz55l32rdyH/gyhVJV24AKEekPLAG8krWekUFELgEuARg2bFi353v37s17773HYYcdZtMNQ0BVee+99+jdu3eum2IKjaVaA+HrPHlVbRaR54GTgP4i0tPtzQ8BtsV4zUJgITgDr9HPDxkyhIaGBpqaus8vNbnRu3dvhgyxXKkx+SDjIC8iA4FWN8CXA6cCtwH/A3wFZ4bNLOCJdI5fWlrKiBEjMm2mMcYUJT968kcCD4hICU6Ov1pVnxSR14GHReQWYC1wnw/nMsYYk4KMg7yq1gPjPLa/DZyY6fGNMcVh6dpGW+gkAIVVu8YYk5eWrm3k+sfX09LaBkBjcwvXP74ewAJ9hizIG1NEwtpbXvDMps4A36GltY0Fz2wKRfvymQV5Y4pEmHvL25pbUtpukldYBcqMMTHF6y3n2qD+5SltN8mzIG9MkQhzb3n2tFGUl5Z02VZeWsLsaaNy1KLCYekaY0LK7/z5oP7lNHoE9DD0ljuuK4zjBfnOgrwxIRRE/nz2tFFdjgnh6i3PGDfYgnoALF1jTAgFkT+fMW4wP545hsH9yxFgcP9yfjxzjAXWAmc9eWNCKKj8ufWWi4/15I0JIZttYvxiQd6YELLZJsYvlq4xJoRstonxiwV5Y0LK8ufGDxbkTU6EtYZKWGT991NfDctvht0NzuLZp9xoqzQVCAvyJuvCXEMlDLL++6mvhmVXQKs7c2f3VucxWKAvADbwarIuzDVUwiCd38/StY1MnL+CEXOeYuL8FSxd25j4RPXVcPtoePzbBwJ8h9YWp2dv8p715E3WhbmGSrr8TK+k+vtJq+f/5NVQez/QbVnlA3Y3JN1mE17WkzdZV2hzwDuCbGNzC8qBIJtUb9pDqr+flHv+9dWJAzywp/yIhG014WdB3mRdNuaAp5W+SJPf6adUfz8pfzNafjMJA7yW8ZPW8xK2tRisrrmHd+d9lva5/Xh33mdZXXNPrpuUEgvyJuuCrqHid886Eb/TT6n+flLq+ddXOwOrMahCQ/sA5rR+iwc+siWaV9fcw+i6GziCJnoIHEETo+tuyKtAbzl5kxNBzgHP9lJyQZTwTeX3k3R1yY5ZNDG0K1zVejk17ZMA58Ol2A1ds4By2ddlW7nsY+iaBVB1aY5alRrryZuCk+2B3VyXIEi657/85u6zaFztCg+2ndoZ4K2EguNwbYqxfWeWW5I+68mbgpPtxTHCUIIgqZ5/nNkydeN/wsLXRyJ2c1oXO2QgR9A90O+QAeTLsLQFeVNwcrE4Rl6UIOg3xDsf328oE6ouZVVV9psUdltPmE2/uhu6pGxatIyt42fnTZC3dI0pOLY4Rgyn3AilUd9mSsud7cbThKpL2TD+Ft5lIO0qvMtANoy/hQl5ko8HENX4U6myqbKyUmtra3PdDGMKl9WoKUgiUqeqlV7PWbrGmGJy/LkW1ItMxukaERkqIv8jIm+IyEYRudLdfqiIPCsib7p/HpJ5c40xxqTCj5z8fuDfVfUY4CTgX0XkWGAOsFxVRwLL3cfGmCzI5h2/JtwyTteo6nZgu/v3D0XkDWAwcBbwJXe3B4DngesyPZ8xQSik+vYdd/xObVvJ4rJqBrXsZPvSAazeem1eDRgaf/iakxeR4cA44BXg0+4HAKq6XUQOj/GaS4BLAIYNG+Znc4xJSqHVt1/wzCbm6L1cUPocPcTZNpidHLrmhzD8EMvJFxnfplCKSF/gMeAqVf0g2dep6kJVrVTVyoEDB/rVHFOggkhDFFp9+8oPnuWCkgMBvkM5e32pEW+poPziS09eREpxAvwiVX3c3fwPETnS7cUfCezw41ymeAXV4y60+vY3lT0Yu/eWYY34QvvWUwz8mF0jwH3AG6r684inaoBZ7t9nAU9kei6TvkLofQXV4y60+vb9+DDOk0MyOnahfespBn6kayYCFwBTRGSd+/MvwHxgqoi8CUx1H5scyHbp3aAE1ePOdYExv0m8JzO5u7W+msV7vs3bvb7Gi2VXUNXjxc6n8vVbTzHwY3bNi8T+d3VKpsc3mct26d2gBFV4LAwFxnxVfii07Oq+veyg9Add3TLFQ3o4v/8hspP5pb+FVqhpn5T2e1BIs5rCyu54LQKFknMOsvBYXhQYS9bpt8ET/wptEXXQS8pg+h1ddkspwHqUKe4j+7i2ZzXP6uS03gPL72eHFSgrAoWSc7bCY0k6/lw46y7oNxQQ58+z7urSi085hRdjwHZQj/foXdqD7y1el/JYTz7m9/NxKUDryReBXJTeDUqYe9yJesaBpiYeqIItKw88HjEZvrch5u4pp/BilCnerofx/p5WIPWeeL59w+xYCrBc9oG7FGC/uhtYDaG+ycx68kXAesDBS9Qz9nr+e4vXccPS9ZmfPDrAg/P4gdgF4lMOsB5lilvoxW2tXXP8qfTE8+0bZtylAEPMevJFIsw94EKQqGfs9bwCi/7yv1R+5tDM3pvoAJ9oO2kMYnekeiLKFM9pOrNzucBIyfbE8+0b5uHa5DnFJOxLAVpP3hgfJOoZx3peIeUcdOQ9D4/ddB7prAiR1rTR4891UkDzmuF7G6g9eKrnbsn2xPPtG+YO8b4jf4cMyHJLUmM9eWN8kKhnHOt5SC0HHT0j5az2PyFxJ8Z782PaqB898Xz6hpmvSwFakDfGB4kC3uxpo/je4nWeve5UctCRaZ+bet5PCe2xdx4xOe6x0g6wT14Ndb9nhrZR1bMHS3r+M9e0fKPg57lPqLqU1Ti5+cN1JztkAFvHzw71oCtYkDfGF4l6xjPGDab277tY9Jf/7RLoU+35bmtuoarHi/yo9D4OYm/8XvysmjSuJIEnr4ba+zof9tB2zuaPnD1pMEz/eZwXFoYJVZeCG9SPcH/CzoK8MT5J1DO+ZcYYKj9zaEYpkll9X+X7rQspk/3xd6y8OOljpqTu97G3F0GQz0cW5I3Jokxz0NfL75ML8EEFXG1LbbvJOQvyxmSBXzdClbU2x99BSoLtUUuJd0CXku7bTCjYFEpjApbVKqDjL/T/mMkcP+jzmrRZkDcmYH7WaNnV3tdzuyps/sz5wefFp//czfdHjPiWHQTDTgr2vCZtFuSNCVjaNVrqq+H20TCvv/NnfTW/LPsWe7VramSvlnBl6+VMf/vL2VkjYNhJUNr7wON9H8OyK5z2psrjGo2/LMgbE7C0arQ8eTU8folbFEydP5ddwRnHD+IH+h0a2gfQrkJD+wBmt15KTfuk7FVw9Cg7TGtL6uvHujXqo6/RAr2/bODVhF6+LyzhdaNUaYnw8d79jJjzVPdrqq+G2vsh+tap1hYmbP4ljV9+hkmLv+h5rsbmFpaubWTGuMG+/946jvdCy9Zui4QDqa8fG+/DIt3FTUw3FuRNqIV9YYlkAmn0jVL9+5Ty0Sf7aW6JUaJ3+c10C/Addjd0FjyLVSbh+sfXU/v3XTxW1+jb7y3yfdhWNoAh4lGUK9X1Y2N9KGS42LjpytI1JtT8GLQMahHzVGbNzBg3mFVzprBl/hn0KetJa3vXIN7lmuIFOTeQehUYizzWf72y1dcFOSLfh5/sP5c9WtZ1h9Ly1NePjfWhkOFi46YrC/Im1DJdWCLI6YuxPoBuWrYx7usSXlPMICedgXTGuMGcPT52j7xNvb8JpLsgR+TratonMaf1W53jAvQbCmfemXqKxaNGfVofFiYuC/Im1DJdWCLdbwLJ9P5jBcz397TG/RDpV17qub3zmryCHwKVF3UG0qVrG3msLvY5YtW0SXdBjujX1bRPYtK+Ozm5/HGn/HA6OfTjz3U+HCKXKUznw8LEZUHehFpadc8jpPNNINnef7yAGetDZOnaRj7e55QlqOrxIi+WXcHbvb7Gql5XcMexbzo7eQW/mQu7zIH3+vCKJDiDu5EyWZAj0/chpqga9Rbg/WcDrybUMq17nvIKSCS//unsaaO4avE6z2PE+hBZ8MwmTtcX+FGv+zmITzp73IPZyeD1c2H4IU6g6/iJIVHapV3h4LKeHNSrpy+za/yoP29yw4K8Cb1Minqls7BFsr3/GeMGM69mY+csmUixPkQqP3iWBaUxqkimMH0w3iIkHXa3tLJu7j8nPFayOt6HjhlF31u8jnVPLeTa0sX0aXnXGUs45UbrjYeMpWtMQUtniblUxgHmVR2XUhrjprIH41eRTHL6YLzZNR2CWBA7MpV1Zo8Xubb11/Rp2Y7dzBRe1pM3BS/VbwKp9P5TSmM8eTX9+DD+yZOcPhh5Xq8efWkPCWRB7MhU1rU9q+kTsRQeYDczhZAFeWOipJp/TupDpL4aau8j/nKsktL0wcj0yexHX6O1LWLaZBrrviYjMmU1yOuGKLCbmULGlyAvIvcD04Edqjra3XYosBgYDrwDnKuq7/txPmOC5vsC009fl3ifiCmSqVjwzKauAR5obdNuA8V+iBwL2KY+3flqAuVXTv73wGlR2+YAy1V1JLDcfWxM8amvhpZd8feZeW/aZYIzvWEsFZFjAb7d+WoC5UuQV9U/A9H/is8CHnD//gAww49zGZN3EvXiyw/NKIed6Q1jqYgcyF7WPomflF7OnvIjsZuZwivInPynVXU7gKpuF5HDvXYSkUuASwCGDRsWYHOMyYEnryTEqWUAABYtSURBVE7Qi+8Bp9+W0SnSmSaaia6prDOAmwI5j/FHzgdeVXUhsBCgsrIyRuk9Y8LNsxplySq3ZHAsAjPvybjnazcqmXiCDPL/EJEj3V78kcCOAM9lTM7EKof8z31vpE+sksHglCrwKbXh+0CxKRhBBvkaYBYw3/3ziQDPZUzWRPfaP967v1sZhKltKylv2R77IBnm4Y1Jll9TKP8L+BIwQEQagLk4wb1aRC4G/hc4x49zmcKUL6s/efXao1X1eJH5pb+NOVW9XaHumDlMCLCdQcuX98v4FORV9asxnjrFj+Obwhb21Z8iJar+CDHuBHW1KzzYdioLXx/JqqogWhi8fHq/jNWuMSHgx+pPmUhl5ahk5p4PjnEnqCpc1Xo5c/dfFMgc9mzJ9ftlUpPz2TXGZPNmnmip9kpjVX88pE8pfcp6UvnBs6h4VxVo1AHUtE/qPE6+yuX7ZVJnPXmTc9m8mSdaqr3SWItnzD3zOFbNmcIvBi7z/E/Vrs4doh37z542KrC1Z4OWy/fLpM6CvMm5TFYdyjRQptorjVu6uL7aKbfrRWBZ+6TO/YHA1p4NWmCrRJlAWLrG5FwyN/N4zeYAMh4ATGflKM856fXVTi31GHr0G8qWeWd0Pp44f0VSq0+Fkd18lV9EY6zqnguVlZVaW1ub62aYkInOm4PTc+xd2oP393RflWlw/3JWzZmS0bETLSzSzW0jYpcvKC3vVtNlxJynPG+TEmDL/DM8njEmNhGpU9VKr+esJ29CL1bePNZUxlQGAH3plSaqMulRtCudbxDGpMOCvI/sBpFgpDprwytQRr43/cpLEYHmPa2d71OyPX9Py66K/Vy/oZ53tma7qJgpXhbkfWI3iAQnVq+3f3kpe/e3JwyU0e9N5MLbGb9P9dXQ+nHs52PUVre8tskWy8n7ZOL8FZ6BKJX8cNDC8E0jnTYsXdvIi0t+zVU8zCDZyTYdwB2cz6QvXw4kDpSx3ptIab9Pt4+OPaMGYN7u1I9pTIosJ58FYb9BJAzfNNJtw4zGn3FWjwProw6RnSzgV8ifHoLTb2PGnPiFvpJ5D9J+n+IF+PJD0zumMT6yefI+CfsNImG4FT2tNtRXQ+393e4gFXAGO5dd4ewTRzLvQVrvU301cVfMznAxEGP8YEHeJ2G/QSQM3zTSasPymyFeTfbWFnj82zCvH9x0qLMSUxSv9yZS2u9TvLZVXmylhE0oWLrGJ2EfSEt3yp6fefyk21Bf7QTQ3Q3EDfDRtA1q74P33oJZNV2e6tWzR+e3iD6lPehVWtJldk1a17S7IfZzaS7KbYzfLMj7KMyr86QzZc/vPH5SbXjyanfJvAwmBGxZ6Rxn+s89b3ZShLlnHpf5e9VviHdOvt/QzI5rjI8sXVMk4tZcicHvPH7CNjx5tdMTzyTAu7T2PqivTvoa0qqBc8qNzt2skUrLY06bNCYXrCdfRFL9phErV97Y3MLStY1p9YRjtqEzwMcjUH6I89eW94n3YSDAnqdvZFvzzzyfj7y2lL+xRKaTyg+BnuVOe/oNcQK85eJNiFiQNzHFyqED/k6/dGfQxNVvKHxvQ9dtCT4Yere8m9Q4QLzefsxCZK3uMVt2Ob13HxflNsZPlq4xMcWblZJM2ibpFEiiGTSIdwpk+s9hxGRi3c8nqrz4yZdZ0+tSqnq82Lk9ehwgpVk/y28+EOA7tLa412BM+FiQNzF15NBjiTf1sSMFklS99HizVAAqL4rdS55Vw+M9TqPdI9CLu0LTofIhPy1byFk9XvQci0jpHodYbU10DcbkiAV5E9eMcYMZnMaNXikN2vYbErsBlRcnnI5YUvVzrtXv0tA+IGavvoz9/GLgMlbNmdItBZPSPQ6x2hrvGozJIQvyWZKvS71Bejd6pZQC8ZqlgiQV4MH5IJr05cs5r8+98efl7N7qeXdsUjOP6qsj6tRE3eVqM2pMiFmBsizwbWGKHEr1pqikCrZFz1KBjGepvDvvsxxBU8znFZARk7vdLBVX9GAr4AR6dQaEbUaNyTErUJZjKc3eCKlUp18mvPEp+qYnn2ap/HjfOfyk9B56ifeCIgLolpWIe7NUUrwGWzsCfPSMH8JR7dOYDpauyYIw1I3JtoQLXnvd1erDLJXag6cyu/VS3mvvG3vWDUDd75I/aAqDrSkNOBuTBdaTz4JiXeotZu8/3pTJDGepON8g9lGzbxJben0t9o7a7nzYJPOtIWb5gu6DrYXwrc0UlsB78iJymohsEpG3RGRO0OcLo1BUqKyvhv8Y6FRr7Ph5oCp7548UJ5C/y4CMBqcjv0G8T9/4O3dUr/T6PXQMtM7rD/s+pk269of2l/T2HGwtxm9tJtwCDfIiUgLcBZwOHAt8VUSODfKcYZRO3Rhf1VfD45dA276u27eszH6gr68G8f5n167wo33nZJzmmDFuMKvmTGHz+Btp0zj13jtE/x46Blp3bwUUWnbR1g7vtfelXYWG9gHMaf0WS9smdjtU2NcVMMUn0Nk1IvIFYJ6qTnMfXw+gqj/22r9QZ9fkXFiWqPOcpeJoBx7cfypz91/UZXumyyeurrmH8Wuu9aU309A+gEn77ozbNq+ZVO48HAbbIKwJSC5n1wwGIqNLA/D5gM9pooXlbkzPWSqAlPC9vZfyRPukbk9lmuaYUHUpDD8k5odLKgbJe10ee7Utcl2BxuaWzgAPtri7yY2gc/Je35W7fHUQkUtEpFZEapuaYs9vzguRedzbRydcli5rbYmRHsm6WB822k7twVM9n/IlzXH8uXDmnRnXed+mh3V5HKttHemiwf3Luw0vZ3vJRWOC/t/fAET+zxoCbIvcQVUXqmqlqlYOHDgw4OYEKDqPu3trUuuPZqUt6j1nHIARk7PWrHglAQIfnD7+XGdOe5rX26Jl/GT/gZk4ybTNBmFNGAQd5FcDI0VkhIiUAecDKdxqmEfCVJ0wVlokWqp3fmYqziIbWRucnlWTYqAX6DeUDeNvoe7gqSm1zQZhTRgEmpNX1f0i8l3gGaAEuF9VNwZ5zpwJU3XCmOcUmNec1aZ00TEnvaOUQVT5gmwtn7j0+N+wYPsmKj94lgVlCyljv/eOEbVzJgCrUpyIlM6Si8b4LfCboVT1v4H/Dvo8OZfCDTNF1ZZox5+b0zovkbNfGpmE7oN5pQ9yiHx4YABJesD4b2a8GHfYF3c3xcHuePXLKTd2n8GRq+qEYWpLyETfkVrTPomavZMynqoZS5gXdzfFwYK8XxKkIoqyLZFVJnPUhuhiYbGWM8znwVAriGbisSDvpxynIrrIdVuib3zqmG3U0bYs8FqgO3LeeqR8HQxNeRFyU3RCMoG6AOViznyY5umHYLaRV7EwpfvNG/k8GJrSClymKFlPPgi56MUGdc50Uy4ZzDbyK/0QKwXTUWKgENIbNhffJGJBPgjxerFBBfkgzpnJB0eaM3z8TD/EysEHNciaC8Vaxtokz9I1QcjFnPkgzplOyiWFtVC91r31M/0Q6y7af/rcwJytt+v3Wr+hKGNtQs2CfBDi3L6f7XM2tB+WfjCJ8QGhzVt57Kbzuh/zyaudksadPfiIDHi/oU79GPcbQKwVlPyc/eJ1F+3Z4wfzWF1jTlZuCmLVqJyXsTahZwt5B8GrpG5peZcgl41z7tEy5rR+i5r2SektHB6nRLEqtCK8Nv42p9JjR816r7krHmuhxlrou0SENo9/k36lWJJaYDwguTy3KWzxSg1bTz4IXaoeSrdebNDnbOfAwhY1bvnetFIeXrVmXCJQJkrlmmud1ZWWXEoqS/rF6pm3qQaafsjlQKUNkppcsIHXoORinrp7zqPmPOUZbhubW5g4f0Xys0o62v/4t2Pu0pl11/bYx/FIJcUbFJ09bVRgN/fkcqDSBklNLlhPPmg5mLseK2gIpJ4PPv5ckJL4+8QlnuUU4g0YdtRj3zL/DFbNmeJrfjmXA5U2SGpywYJ8kHJUY94rmHjd6Zl0Cmf8hbESMQkIVF7k+Y0mVwOGuRyotEFSkws28BqkWAOXHgORfku2ZosAW+afkfB4b80bzVG6FUm0LraUOKmbXNbuMabI5HKN1+KWwxrz0dUPY83sSDYfPPWT25jX836+XrKcHm6/vlvA71EKM35tgd2YELF0TZByMV8+htnTRlHao2tULu0hSeeDB/UvZ+7+izhq7yJG7H2IK1svZ5f2RdVJA+0t7WcB3pgQsiAfpDjL3eVEdM87Ueolwj99bmCX3WvaJ3HC3oWM2PsQIz55iIpPFrK0baIfrTTG+MiCfJByMV8+hgXPbKK1rev4S2ubJjXwunRtI4/VNcYdfLXKh8aEk+Xkg5bruu6uTG7E8aonk+6xjDHZZUG+SGRyI06ywXtQ/3JbpciYkLF0TZHI5EacZD4IOqo7+l2AyxiTGQvyRSKTG3G8PiAiHdKnlB/PHMNT9dttlSJjQsbSNdkSZ4WlbKU4oufOp/I6cHLzjc0tnZUiB0e0denaRt7f0+r5esvVG5M7FuSzIc4KS0vbJubFQsyJPiDi9datAJcxuWPpmmyIs8JSoSzEHK+3bgW4jMkdC/LZEKe8QaHUGI/VW+9fXhqqbyTGFJuCDfKra+7h3XmfpX1uP96d91lW19yTu8bEKW8QKzjmW4oj1uydeVXH5ahFxhjIMMiLyDkislFE2kWkMuq560XkLRHZJCLTMmtmalbX3MPouhs4giZ6CBxBE6PrbshdoI9T3qBQaoznaxldvxfWNiZsMh143QDMBLpETxE5FjgfOA4YBDwnIkerauLbJn0wdM0CymVfl23lso+haxZA1aXZaEJXHXe8esyumeHuUgg3EKU7eydXOhbWDvugtzGZyCjIq+obANK9yPhZwMOquhfYIiJvAScCL2dyvmQdrk2exbcO153ZOL23OOUN8i04Fop4g972fphCEVROfjAQuVpGg7utGxG5RERqRaS2qanJl5PvkIExtg/w5fimMBTKoLcx8SQM8iLynIhs8Pg5K97LPLZ5FjFU1YWqWqmqlQMHegfnVG09YTYtWtZlW4uWsfWE2b4c3xSGQhn0NiaehOkaVT01jeM2AEMjHg8BtqVxnLRMqLqU1Ti5+cN1JztkAFvHz2ZCLvLxWVSsxcHSve7Z00Z1yclDfg56GxNPUHe81gAPicjPcQZeRwKvBnQuTxOqLu0cZD3C/SlkxTqImMl1R5ZrKLYPRlM8MgryIvJl4JfAQOApEVmnqtNUdaOIVAOvA/uBf83WzJpiVayDiJletw16m0KX6eyaJcCSGM/dCtyayfFN8op1ELFYr9uYZBXsHa/FplgHEYv1uo1JlgX5gGXrjspCuXM2VcV63cYkq6hKDa+uucedcdPEDhnI1hOCm3GzdG0jNy3b2KXGepCDocU6iFis121MskTVc/p6TlRWVmptbW0gx+6oZxNZ7qBFy9gw/hbfA330jI9og/uXs2rOFF/PaYwpXiJSp6qVXs8VTbombj0bn3nN+Ihkg4LGmGwpmnRNNuvZJAriNihY+Ir1xjQTPkXTk89mPZt4QdxzULC+Gm4fDfP6O3/WV/veJpM9Hem6xuYWlANjMVbG2ORC0QT5bNaz8ZrxAc4qSd1qrHes/7p7K6DOn0svh9tGWNDPU4WypKMpDEWTrslmPZuUZnx4rf/a3gotu5y/Ryz6HatUsQkXu0HLhEnRBHnIbj2bpG+Xj7X+ayR30W8L8vlhUP9yGj0Cuo3FmFwomnRNaMVa/zVaMh8GJhTsBi0TJhbkc81r/VcvyX4YmJzL1/VuTWEqqnRNKEWv/1p+COz7CNoi5vS7i36b/GHVLU1YWJAPg+j1X+urPRf9NsaYVFmQD6M4i36bruymI2PisyBv8laxroZlTCosyIeQ9U6TU6yrYRmTCgvyIZNPvdNcfxjZTUfGJGZTKEMmX26JD0N9FlsVypjELMiHTL70TsPwYWQ3HRmTmKVrQiZfbokP6sMolRSQrQplTGIW5ENm9rRR3VaVCmPvNIgPo3TGI+ymI2Pis3RNyOTLLfFBpErCkAIyptBYTz6E8qF3GkSqJF/GI4zJJxbkTdr8/jDKl/EIY/KJpWtMaNhsGWP8Zz15Exo2W8YY/2UU5EVkAXAmsA/YDHxTVZvd564HLgbagCtU9ZkM22qKQD6MRxiTTzJN1zwLjFbV44G/AdcDiMixwPnAccBpwK9FpPvK1sYYYwKVUZBX1T+p6n734V+AjuWLzgIeVtW9qroFeAs4MZNzGWOMSZ2fA68XAU+7fx8MbI14rsHd1o2IXCIitSJS29TU5GNzjDHGJMzJi8hzwBEeT/1AVZ9w9/kBsB9Y1PEyj/3V6/iquhBYCFBZWem5jzHGmPQkDPKqemq850VkFjAdOEVVO4J0AzA0YrchwLZ0G2mMMSY9ciAup/FikdOAnwOTVbUpYvtxwEM4efhBwHJgpKq2eR7owOuagL+n3aDcGwDszHUjssCus7DYdea/z6jqQK8nMg3ybwG9gPfcTX9R1cvc536Ak6ffD1ylqk97H6VwiEitqlbmuh1Bs+ssLHadhS2jefKq+tk4z90K3JrJ8Y0xxmTGyhoYY0wBsyDvr4W5bkCW2HUWFrvOApZRTt4YY0y4WU/eGGMKmAV5Y4wpYBbkMyAih4rIsyLypvvnIR77VIjIyyKyUUTqReS8XLQ1E8lcp7vfH0WkWUSezHYbMyEip4nIJhF5S0TmeDzfS0QWu8+/IiLDs9/KzCVxnf9XRNaIyH4R+Uou2uiHJK7zahF53f3/uFxEPpOLdmaLBfnMzAGWq+pInBu+uv2DAvYA31DVjoqcd4hI/yy20Q/JXCfAAuCCrLXKB2511LuA04Fjga+6VVQjXQy8704Zvh24LbutzFyS1/m/wIU4NzLmpSSvcy1Q6VbPfRT4SXZbmV0W5DNzFvCA+/cHgBnRO6jq31T1Tffv24AdgOedaSGW8DoBVHU58GG2GuWTE4G3VPVtVd0HPIxzvZEir/9R4BQR8arPFGYJr1NV31HVeqA9Fw30STLX+T+qusd9GFk9tyBZkM/Mp1V1O4D75+HxdhaRE4EynAVW8klK15lnkqmY2rmPW1p7N3BYVlrnn6Qrw+a5VK/zYg5Uzy1ItvxfAvGqcKZ4nCOBB4FZqhq6npJf15mHkqmYmnRV1RArhGtIRtLXKSJfByqByYG2KMcsyCcQrwqniPxDRI5U1e1uEN8RY7+DgaeAG1T1LwE1NSN+XGeeSqZiasc+DSLSE+gH7MpO83xTLJVhk7pOETkVpwMzWVX3ZqltOWHpmszUALPcv88CnojeQUTKgCXAH1T1kSy2zU8JrzOPrQZGisgI9706H+d6I0Ve/1eAFRFltfNFMtdZCBJep4iMA+4BqlS1kDos3lTVftL8wcnLLgfedP881N1eCfzW/fvXgVZgXcRPRa7b7vd1uo9fAJqAFpwe1bRctz3J6/sXnDWKN+MshgNwM04QAOgNPIKzjOWrwP/JdZsDus4J7vv2MU5l2Y25bnNA1/kc8I+I/481uW5zkD9W1sAYYwqYpWuMMaaAWZA3xpgCZkHeGGMKmAV5Y4wpYBbkjTGmgFmQN8aYAmZB3hhjCtj/B10BIaZv6fX3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_regression\n",
    "np.random.seed(1)\n",
    "x, y = make_regression(n_samples=100, n_features=1, noise=10,effective_rank=1, bias=5)\n",
    "x = pd.DataFrame(x)\n",
    "# Create polynomial model\n",
    "for i in range(0,20):\n",
    "    x['x'+str(i)] = x.iloc[:,[0]]**i\n",
    "x=x.drop(0, axis=1)\n",
    "results = sm.OLS(y, x.iloc[:, 1:]).fit()\n",
    "plt.scatter(x['x1'],y)\n",
    "plt.scatter(x.x1, results.predict())\n",
    "plt.title(r'Example of overfitting: $y=a_1+\\cdots+a_{20}$')\n",
    "plt.legend(['Observations','Fitted'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2 Overfitting & condition number\n",
    "\n",
    "**A large condition number $\\rightarrow$ risk of overfitting**\n",
    "\n",
    "- When the condition number is big you will likely overfit if there is any type of noise or rounding error.\n",
    "\n",
    "\n",
    "- The rounding error blows up when inverting matrixs\n",
    "\n",
    "\n",
    "- If the condition number is very big and you start solving the linear system by row elimination you a little rounding error which become part of the solution.\n",
    "\n",
    "\n",
    "- The condition number to a linear system is usually large so the errors explode\n",
    "\n",
    "\n",
    "- If you start solving the system you are accumulating the rounding errors, and when you invert the system this blows up the error\n",
    "\n",
    "\n",
    "- If you have a well conditioned problem then it is difficult to overfit because the X has weak dependence on the noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2.1 Python example: Overfitting & condition number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition number: 9.285999564131141e+17\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxcVf3/8dcnW5Nu6Zau6QqlG6UUQkVQKHsBAQF/UgQFRHHDBVf8yhcR9esu+vWLyiKrKGJBLFApIBQUgSbdNwpd0jZdQku2pku2+fz+uDd0GCbJJM1ksryfj8c8cu8959z5zHQ6n7n33HuOuTsiIiKx0lIdgIiIdE5KECIiEpcShIiIxKUEISIicSlBiIhIXEoQIiISlxKEiIjEpQQhXZaZjTMzN7OMBOpebWb/bqJsjJlVm1l6K5//FjP7Y2vatAcL3Gtm5Wa22Mw+aGbrOzoO6f6UIKRDmFmxmdWa2ZCY7cvDL/lxqYkM3H2ru/d194ZUxdAowaT3AeAsIN/dZ7n7v9x9UtQ+is3szASeq7eZfdPMlphZmZltM7OHzGxmC+1GmNl8M9sR79/OzO4L/62rox6tSr7SOShBSEfaDFzeuGJm04Gc1IXTZY0Fit19X1t3YGbDgVeAI4BPAsOBKcBjwB/N7JpmmkeAp4FLm6nz0zDp9u0syVdaTwlCOtKDwCei1q8CHoiuYGa5ZvaAme02sy1mdpOZpYVl6Wb2czPbY2abgPPjtP2Dme00s+1m9oNEfrnG/moPT0dtMrO9ZrbZzK5opnm2mf0lrLvUzGZE7XekmT0avpbNZvalqLJZZlZkZlVmVmpmvwyLXgr/VoS/vN8fE+u1wN3A+8Py75nZbDMrCcsfBMYAT4Tl32wi7j8Bv3b3z7j7Cnevdfdqd3+U4Ajla2Z2RLyG7l7q7r8FCpt5X6Q7cHc99Ej6AygGzgTWE/xSTQe2EfwadmBcWO8B4O9AP2Ac8AZwbVj2WeB1YDQwCHghbJsRlj8O3AH0AYYCi4HPhGVXA/9uIrZxjfsJ21YBk8KyEcC0JtrdAtQBHwEyga8THCVlEvz4WgLcDGQBE4BNwDlh21eAj4fLfYETY2Np5r1812sBZgMlse91M+1PBRaGy6OB54EdwO3A4nD7FcDPW/g3zYj+t4vafh9QFj6WAJem+vOnR9seOoKQjtZ4FHEWwZf99saC8Nf+ZcC33X2vuxcDvwA+Hlb5KPArd9/m7mXAj6LaDgPOBb7i7vvc/S3gNmBuG2KMAEebWY6773T3Nc3UXeLu89y9DvglkA2cCJwA5Ln7rR78Ot8E3BUVTx1wpJkN8eCX+6ttiLOtzgIeDpd/DvyH4KjjSaAg3L4cmNzG/f8vMJEgSf83cJ+ZndzmaCVllCCkoz0IfIzgV/ADMWVDCH5tb4natgUYFS6PJDjqiC5rNJbgl/tOM6swswqCo4mhrQnOg/P6lxEcrew0s6fMrLkvynficfcIUBLGORYY2RhLGM9/AcPC6tcCRwGvm1mhmX2oNXEepqEcSszTgT+5e727/wPYE24fHVWnVdx9qbu/He5zAfAQcMnhBi0dTwlCOpS7byE4DXMeQYdotD0Ev6zHRm0bw6Evqp0EX1zRZY22ATXAEHcfED76u/u0NsS40N3PIji99DrBL/+mvBNP2FeST3C6ZhuwOSqWAe7ez93PC5/jTXe/nODL+ifAPDPrQ3DK5nC1tI89BK8NYBXwMTPLMLM5wBAzOxL4IUFfR3twwNppX9KBlCAkFa4FTveYq3A8uNLlEeCHZtbPzMYCXwUa7zV4BPiSmeWb2UDgxqi2O4FngF+YWX8zSzOzI8zs1NYEZmbDzOzC8Mu6BqgGmrsC53gzuyTs4P5K2OZVgv6PKjP7lpnlhB3sR5vZCeHzXGlmeeFRR0W4rwZgN8EprgmtiTtGaQvtnyfoN4Gg3+QkgoR2GfAiQWL4prs32QltZtlAr3C1V7jeWPYRM+sb/hucDVwJzG/ri5HUUYKQDufuG929qIniLwL7CDp0/01wtc09YdldwEJgBbCU9x6BfILgFNVaoByYx6FfyolKA75GcBRQRtCh+/lm6v+d4Iu1nKCv5BJ3rwuT3QXAsQRHTHsIvnhzw3ZzgDVmVg38Gpjr7gfdfT/Br/eXw1NTJ7Yyfgj6Zm4K2389ttDdnwMGmtkVYX/O6e4+wt2vcffTCDq4/9nCcxwgSJ4QHGUdiCr7MsFRXwXwM+DT7r6oDa9DUszcNaOcSE9jZqMIjrieJEi8mwiS6ReBGe5+bgrDk05CRxAiPZC7bwfeDxwkuDy4jOD0UgaHrhqTHk5HECIiEpeOIEREJK4WR8HsKoYMGeLjxo1LdRgikmzrw4FrJ01qvp4kZMmSJXvcPS9eWbdJEOPGjaOoqKkLY0Sk25g9O/i7aFEqo+g2zGxLU2U6xSQiInEpQYiISFxKECIiEpcShIiIxKUEISIicSlBiIhIXEoQIiISV7e5D0JEegYH9h6s47W1pTREnIg79REnEnEaGh9NbGuIWu+Vmcb0Ubkckz+A3JzMVL+sTkkJQkS6lB0VB9hWtp9PP9B+N8ZOGNKHY0cPYEb4mDKiH70y0ttt/12VEoSIdBnl+2rZUXGA3N5ZzL/+ZNLTLHiYHVqOsy0tzchIM9LsUNnemnpWlVSyoqSC5dsq+NeGPTy2LJi8MCs9jSkj+3Nsfi7HjhnAjPwBjBvch7S0njUxnhKEiHQZv120gTMizthBvemdP+Cw9pWbk8kHJg7hAxOHAODu7Kw8yIptQcJYvq2Cvy4p4f5XgpEo+mdnMGP0gOBIIz840sjr16u5p+jylCBEpEvYXnGA+1/ZwmX9etE7q/1P/5gZIwfkMHJADudODyYibIg4b761N0walSzfVsFvF22kIRJMk5Cbk0njQUXjxAnRMyhET6fwrokVPO4isdMvvLuMmLJDG47JH8Ajn3l/yy+ylZQgRKRL+NWzb4DD6IG9O+w509OMycP7M3l4fy47Idi2v7aeNTuqWLGtgq1l+99Vv/EElFn8U1HRmw1rYnszbWL227g2ckBOC6+kbZKaIMxsDsF8u+nA3e7+45jysQTzDecRzGh1pbuXhGUNwKqw6lZ3vzCZsYpI5/Vm6V4eXVrCNSePp9fLqb06v3dWBieMG8QJ4walNI6OkLR32szSgduBc4GpwOVmNjWm2s+BB9z9GOBWgsnWGx1w92PDh5KDSA/2s4Xr6ZOVwRdOOzLVofQoyUzFs4AN7r7J3WuBh4GLYupMBf4ZLr8Qp1xEerglW8p5Zm0p150ygUF9slIdTo+SzAQxCtgWtV4Sbou2Arg0XL4Y6Gdmg8P1bDMrMrNXzezD8Z7AzK4L6xTt3r27PWMXkU7A3fnJ068zpG8vrv3g+FSH0+MkM0HE66WJ6Yfn68CpZrYMOBXYDtSHZWPcvQD4GPArMzviPTtzv9PdC9y9IC8v7ox5ItKFLVq/m8Wby/jyGUfSO0vX1HS0ZL7jJcDoqPV8YEd0BXffAVwCYGZ9gUvdvTKqDHffZGaLgJnAxiTGKyKdSCQSHD2MHdybubPGpDqcHimZRxCFwEQzG29mWcBcYH50BTMbYmaNMXyb4IomzGygmfVqrAOcDKxNYqwi0sn8fcV2Xt+1l6+dPYnMdI0rmgpJe9fdvR64HlgIrAMecfc1ZnarmTVelTQbWG9mbwDDgB+G26cARWa2gqDz+sfurgQh0kPU1Dfwi2feYNrI/nwovGlNOl5ST+q5+wJgQcy2m6OW5wHz4rT7DzA9mbGJSOf1p9e2UlJ+gP+5eHqPG/+oM9Fxm4h0KtU19fzf8xs46YjBfDAcJ0lSQwlCRDqVu17axNv7avnmnMlNDlkhHUMJQkQ6jT3VNdz9r02ce/Rwjh19eKO1yuFTghCRTuP/nt/AwfoIXz9nUqpDEZQgRKST2Pr2fh56bQsfLcjniLy+qQ5HUIIQkU7itufeIM2ML59xVKpDkZAShIik3LqdVTy+fDvXnDye4bnZqQ5HQkoQIpJyP336dfr1yuBzp75nyDVJISUIEUmp1za9zQvrd/P5044kt3dmqsORKEoQIpIy7s6Pn36d4f2zufqkcakOR2IoQYhIyjyztpRlWyv4ypkTyc5MT3U4EkMJQkRSor4hws8WrueIvD585Pj8VIcjcShBiEhKPLZ0OxvequYb50wiQ8N5d0r6VxGRDnewroHbnnuDGaMHcM604akOR5qgBCEiHe6BV4rZWXmQb82ZpAH5OjElCBHpUJUH6rj9hY2cclQeJx2h4bw7MyUIEelQd7y4kcoDdXxTA/J1ekoQItJhSqsOcs/Lm7lwxkiOHpWb6nCkBUoQItJhfvKP14lE4Gtna0C+rkAJQkQ6xJItZTy2bDufPmU8Ywf3SXU4kgAlCBFJuoaIc/Pf1zC8fzafn31kqsORBClBiEjS/aVwG2t2VPFf50+hT6+MVIcjCVKCEJGkqthfy88Wvs6s8YO44JgRqQ5HWkEJQkSS6pfPvkHlgTpuuWCaborrYpQgRCRp1u2s4o+vbuHKE8cydWT/VIcjraQEISJJ4e58d/4acnMy+epZuqy1K1KCEJGkeGLlThZvLuPr50xiQO+sVIcjbaAEISLtbn9tPf/z1DqOHtWfuSeMSXU40ka63kxE2t3tL2xgV9VBbr9iJulp6pjuqpJ6BGFmc8xsvZltMLMb45SPNbN/mtlKM1tkZvlRZVeZ2Zvh46pkxiki7ad4zz7uemkzl8wcxfFjB6U6HDkMSUsQZpYO3A6cC0wFLjezqTHVfg484O7HALcCPwrbDgK+C7wPmAV818wGJitWEWk/P3hqLZnpxo3nTk51KHKYknkEMQvY4O6b3L0WeBi4KKbOVOCf4fILUeXnAM+6e5m7lwPPAnOSGKuItIMX1r/Fc+ve4ktnTGRo/+xUhyOHKZkJYhSwLWq9JNwWbQVwabh8MdDPzAYn2BYzu87MisysaPfu3e0WuIi0Xk19A7c+sZYJQ/pwzcnjUx2OtINkJoh4PVMes/514FQzWwacCmwH6hNsi7vf6e4F7l6Ql5d3uPGKyGG49+ViNu/Zx80XTCUrQxdIdgfJvIqpBBgdtZ4P7Iiu4O47gEsAzKwvcKm7V5pZCTA7pu2iJMYqIoehtOogv/nnm5w5ZRizJw1NdTjSTpKZ5guBiWY23syygLnA/OgKZjbEzBpj+DZwT7i8EDjbzAaGndNnh9tEpBP60YJ11EWcmz8Uex2KdGVJSxDuXg9cT/DFvg54xN3XmNmtZnZhWG02sN7M3gCGAT8M25YB3ydIMoXAreE2EelkiorLeHz5Dq774ATGDO6d6nCkHSX1Rjl3XwAsiNl2c9TyPGBeE23v4dARhYh0Qo0TAY3Mzebzpx2R6nCknaknSUTa7M+Lt7J2ZzARUO8sDczQ3ShBiEibVOyv5efPrOfECYM4f7omAuqOlCBEpE1+8cwb7D1Yzy0XaiKg7koJQkRabe2OKh56bQsfP3Esk4drIqDuSglCRFrF3bll/hoG9M7ihjM1EVB3pgQhIq0yf8UOFheX8Y1zJpHbOzPV4UgSKUGISML21dTzPwvWMX1ULh8tGN1yA+nSdF2aiCTs9hc2UFpVw++uPF4TAfUAShAi3Vx9Q4RXN5Xx3LpSqg7UURdx6hsi1DU49ZEI9Q1OXUOE+njb3ykPtlUeqOPS4/I5boymZ+kJlCBEuqH6hgiLN5fx5KqdPL16F2X7asnJTGdIvywy0tLISDMy0tPITLd3lrMz08jolRFuSyMj3chMf3fd3JxMPvWBCal+edJBlCBEuomGiLN4cxlPrdrB06t3sae6lt5Z6ZwxZRjnTx/B7El5ZGempzpM6UKUIES6sIaIU1RcxlOrdrJg1S72VNeQk5nO6VOG8qHpI5g9aSg5WUoK0jZKECJdTCTiLNlazlMrd7Jg1U7e2ltDdmYap08eyvnTR3La5DyNiyTtQp8ikS4gEnGWbSvnyTAplFbV0CsjjdMmDeX8Y0Zw+uSh9Oml/87SvvSJEmlndQ0R3ijdy6qSSlZur2T19kqqDtThgDtE3PFwAt3GZcfDMuCdZX+nTV1DhP21DWRlpDH7qDzOP2YEZ0wZRl8lBUkifbpEDkNDxNm0u5oVJZWsKqlg5fZK1u6ooqY+AkC/7AyOHpnLuMF9MIM0s2DCdQPDSDOwxuW0oCCoF2x7p43BjPwBnDFlKP2ydfeydAwlCJEERSLOlrL9rCypYGVJJatKKlm9o5L9tQ0A9M5K5+hRuXz8xLFMz8/lmPwBjB3UmzTdUCZdlBKESDOWbS3n6TW7WFVSyartlew9WA9Ar4w0po3sz0cLRjN9VC7H5OcyIa+v7i6WbkUJQqQJT63cyZcfXoYZTBnRnwtnjOSY/FymjxrAxGF9yUzXUGbSvSlBiMTx6JISvjFvBceNGcgfrj6B3Byd95eeRwlCJMYfX93CTY+v5uQjB3PXJwp0T4H0WPrki0S5+1+b+MFT6zh98lB+e8VxGppCejQlCBGCWdJ+8/wGfvnsG5w/fQS3XXYsWRnqY5CeTQlCejx35ydPr+f3L27kkuNG8dNLjyFDHdAiShDSs0UizveeWMP9r2zhiveN4fsXHa37FkRCzSYIM7ukuXJ3f6x9wxHpOA0R59uPreSRohI+9YHxfOf8KZgpOYg0aukI4oLw71DgJOD5cP00YBGgBCFdUl1DhK8+soInVuzgS2dM5IYzJyo5iMRoNkG4+zUAZvYkMNXdd4brI4Dbkx+eSPurqW/g+j8t49m1pdx47mQ+e+oRqQ5JpFNKtA9iXGNyCJUCRyUhHpGkOlDbwHUPFvGvN/fwvQuncdVJ41IdkkinleilGovMbKGZXW1mVwFPAS+01MjM5pjZejPbYGY3xikfY2YvmNkyM1tpZueF28eZ2QEzWx4+ft+qVyUSR3VNPVfdu5h/b9jDTy89RslBpAUJHUG4+/VmdjFwSrjpTnf/W3NtzCyd4DTUWUAJUGhm8919bVS1m4BH3P13ZjYVWACMC8s2uvuxib8UkaZV7q/jE/cuZvX2Sn49dyYXzhiZ6pBEOr3WXOa6FNjr7s+ZWW8z6+fue5upPwvY4O6bAMzsYeAiIDpBONA/XM4FdrQiHpGE7Kmu4eN/WMzGt6r53RXHcfa04akOSaRLSOgUk5l9GpgH3BFuGgU83kKzUcC2qPWScFu0W4ArzayE4Ojhi1Fl48NTTy+a2QebiOs6Mysys6Ldu3cn8lKkhymtOshld7zC5j3V3HVVgZKDSCsk2gfxBeBkoArA3d8kuPS1OfGuGfSY9cuB+9w9HzgPeNDM0oCdwBh3nwl8FfiTmfWPaYu73+nuBe5ekJeXl+BLkZ6ipHw/H73jFXZVHuT+a2Zx6lH6jIi0RqIJosbdaxtXzCyD937ZxyoBRket5/PeU0jXAo8AuPsrQDYwxN1r3P3tcPsSYCO6akpaYU91DR+76zXK99Xyx0+9j/dNGJzqkES6nEQTxItm9l9AjpmdBfwVeKKFNoXARDMbb2ZZwFxgfkydrcAZAGY2hSBB7DazvLCTGzObAEwENiUYq/RwB2ob+NT9RZRWHeS+T85i5piBqQ5JpEtKNEHcCOwGVgGfIegvuKm5Bu5eD1wPLATWEVyttMbMbjWzC8NqXwM+bWYrgD8DV7u7E1wttTLcPg/4rLuXte6lSU/UEHG+8pdlrCip4NdzZ3KckoNIm7V4FVP4S/5+d78SuKs1O3f3BQTJJHrbzVHLawn6NmLbPQo82prnEgH44VPrWLimlP/+0FTmHK0OaZHD0eIRhLs3AHnhaSKRTuvelzdzz8ubuebkcVz7gfGpDkeky0v0Pohi4GUzmw/sa9zo7r9MRlAirbVwzS5ufXIt50wbxk3nT011OCLdQqIJYkf4SAP6JS8ckdZbvq2CLz+8jBn5A/jVZTNJ13wOIu0i0aE2vpfsQETaYuvb+7n2vkKG9svm7qsKyMnSHNIi7SWhBGFmecA3gWkEl6IC4O6nJykukRZV7K/l6vsW0+DOvdecwJC+vVIdkki3kuhlrg8BrwPjge8R9EkUJikmkRYdrGvgugeWUFJ2gDs/XsAReX1THZJIt5Noghjs7n8A6tz9RXf/JHBiEuMSaVIk4nxj3koWF5fxi4/OYNb4QakOSaRbSrSTui78u9PMzifosM5PTkgizfvpwvU8sWIH35ozmQs0bLdI0iSaIH5gZrkEdz7/hmCI7huSFpVIEx56bQu/f3EjV7xvDJ89dUKqwxHp1hK9iunJcLESOC154Yg07YXX3+K/H1/N6ZOH8r0Lp2Gmy1lFkinRq5juJc7orWFfhEjSrd5eyRf+tJSpI/vzm8tnkpGeaPeZiLRVoqeYnoxazgYuRrO/SQfZXnGAa+4rZGDvLO656gT69GrNRIgi0laJnmJ618B5ZvZn4LmkRCQSpfJAHdfcu5iDdQ089Kn3MbR/dsuNRKRdtPU4fSIwpj0DEYlVWx/hc39cwuY9+7jjyuM5aphGeRHpSIn2Qewl6IOw8O8u4FtJjEt6OHfnxsdW8p+Nb/OL/zeDk44ckuqQRHqcRE8x6aebdIiGiPOfjXt46NWtPL1mFzeceRSXHq9bbkRSIdEjiOOaK3f3pe0TjvRUG3dX8+iSEv62bDs7Kw/SPzuDL51+JF8648hUhybSYyV6OchvgeOAlQSnmY4BXiO4w9oBDdonrVa5v44nVu5g3pISlm+rID3NOGXiEG46fypnTBlKdqZGZhVJpdZMGPRpd18FYGZHA19396uTFJd0U/UNEf715h7mLSnh2XWl1NZHmDSsH985bwoXzRzJ0H66Skmks0g0QUxuTA4A7r7azI5NUkzSDa3ftZdHlwankHbvrWFg70w+NmsMHzk+n2kj++uuaJFOKNEEsc7M7gb+SHBK6UpgXdKikm6hbF8t85dvZ97SElZvryIjzTh98lAuPT6f0yYNJStDd0OLdGaJJohrgM8BXw7XXwJ+l5SIpEurrqln0fq3eGLFDp5//S3qGpyjR/XnuxdM5cIZIxmsSX1EuoxEL3M9CNwG3GZmg4D8cJsIe6preG5tKQvX7OLlDW9T2xBhSN9eXH3SOC49Pp/Jw/unOkQRaYNEL3NdBFwY1l8O7DazF939q0mMTTqxbWX7WbhmF8+sKaVoSxkRh/yBOXz8/WM5Z9pwjh87kPQ09SuIdGWJnmLKdfcqM/sUcK+7f9fMViYzMOlc3J11O/fyzNpdLFxTyrqdVQBMHt6PL54+kXOmDWfKiH7qbBbpRhJNEBlmNgL4KPCdJMYjnUhDxFm6tZyFq3fxzNpStpbtxwwKxg7kpvOncPbU4YwZ3DvVYYpIkiSaIG4FFgL/dvdCM5sAvJm8sKSjNUScvQfrqDxQx8bd1TyzppTn1pWyp7qWrPQ0Tj5yMJ+bfQRnThlGXj91NIv0BIl2Uv8V+GvU+ibg0mQFJW23r6ae8v21VB6oo3J/8IVfeaCOigOHlhu3Vxw4VG9vTT0eNSVU314ZzJ6UxznThjN7Uh79sjNT96JEJCVaPfOKmS1192bHZpKOtWl3Nf9YvYuFa3axsqSyyXoZaUZuTia5vTPJzckkr28vjszrG27LCv7mZDK8fzYnjB9IrwwNdSHSk7Vlaq6EeyHNbA7wayAduNvdfxxTPga4HxgQ1rnR3ReEZd8GrgUagC+5+8I2xNotuTtrd1axcPUunl6zizdKqwGYMXoAN5x5FMNze4Vf9lnvJIQBOZn0zkpXJ7KIJKwtCeKpRCqZWTpwO3AWUAIUmtl8d18bVe0m4BF3/52ZTQUWAOPC5bnANGAk8JyZHeXuDW2It1uIRJxl2ypYuGYXT6/exday/aQZzBo/iFsumMrZ04YzckBOqsMUkW6k1QnC3W8Kv/znuvtDzVSdBWwI+ysws4eBi4DoBOFA411UuRya5/oi4GF3rwE2m9mGcH+vtDberqy+IcLizWU8vSY4fVRaVUNmunHykUP4/OwjOHPqMIbozmQRSZJmE4SZ9Qe+AIwC5gPPhuvfILhhrrkEMQrYFrVeArwvps4twDNm9kWgD3BmVNtXY9qOihPfdcB1AGPGdI8ZUGvqG3h5wx7+sWoXz60rpXx/HdmZacw+aihzjh7O6VOG0l8dxiLSAVo6gngQKCf45f4pgsSQBVzk7stbaBvvZLfHrF8O3OfuvzCz9wMPhkOJJ9IWd78TuBOgoKDgPeWdUW19hLJ9teypruHtfbWU7avh7epa9lTXsq1sPy++sZvqmnr6ZWdw5pRhnDNtOKcelUdOljqMRaRjtZQgJrj7dIBwNNc9wBh335vAvkuA0VHr+Rw6hdToWmAOgLu/YmbZwJAE23YqW9/ez7by/eyprqFsXy1vV9fydvjlHySCICnsPVgft31mujG0XzYXzBjBOdOGc9IRQzTaqYikVEsJoq5xwd0bzGxzgskBoBCYaGbjge0Enc4fi6mzFTgDuM/MpgDZwG6C01l/MrNfEnRSTwQWJ/i8HW7DW9Wc+csX37UtzWBQn14M7pPF4L5ZHD0qN1juk8WgvlkM7tOLIX2zGNQni8F9e9E/O0NXGIlIp9JSgphhZlXhsgE54boB7u5NDtPp7vVmdj3BHdjpwD3uvsbMbgWK3H0+8DXgLjO7geAU0tXu7sAaM3uEoEO7HvhCZ76C6ZVNbwPw+yuP58ihfRjUpxcDcjJJ02B1ItKFNZsg3P2wTnyH9zQsiNl2c9TyWuDkJtr+EPjh4Tx/RykqLmNY/16cM22YjgJEpNvQSe52UFRcTsG4QUoOItKtKEEcpu0VB9hecYCCsQNTHYqISLtSgjhMRcVlAJwwblCKIxERaV9KEIepqLicPlnpTB7eL9WhiIi0KyWIw1RYXMZxYweSka63UkS6F32rHYaqg3WsL91LwVidXhKR7kcJ4jAs3VKOO5wwTh3UItL9KEEchqLictLTjGPHDEh1KCIi7U4J4jAUFpcxbWR/eme1ZVoNEZHOTQmijWrrIyzfVqH+BxHptpQg2mj1jkpq6iPqfxCRbksJoo2WFJcDcLwShIh0U0oQbVRYXMa4wb0Z2i871aGIiCSFEkQbuDtFW4IB+kREuisliDbYtGcfZftqNUCfiJFqFLIAAA8PSURBVHRrShBt0DhAn44gRKQ7U4Jog8Licgb2zuSIvD6pDkVEJGmUINpgyRZNECQi3Z8SRCvt3lvD5j37dP+DiHR7ShCttGSL+h9EpGdQgmilwuJyemWkcfTI3FSHIiKSVEoQrVRUXMaM0QPIytBbJyLdm77lWmF/bT2rd1Sp/0FEegQliFZYvq2Choir/0FEegQliFYoKi7HDI4boyMIEen+lCBaobC4jEnD+pGbk5nqUEREkk4JIkH1DRGWbinnBJ1eEpEeQgkiQa/v2su+2gYK1EEtIj2EEkSCNECfiPQ0SU0QZjbHzNab2QYzuzFO+W1mtjx8vGFmFVFlDVFl85MZZyIKt5QzMjebUQNyUh2KiEiHyEjWjs0sHbgdOAsoAQrNbL67r22s4+43RNX/IjAzahcH3P3YZMXXGu5OUXEZ7xs/ONWhiIh0mGQeQcwCNrj7JnevBR4GLmqm/uXAn5MYT5uVlB+gtKpGN8iJSI+SzAQxCtgWtV4SbnsPMxsLjAeej9qcbWZFZvaqmX24iXbXhXWKdu/e3V5xv0eRBugTkR4omQki3mQJ3kTducA8d2+I2jbG3QuAjwG/MrMj3rMz9zvdvcDdC/Ly8g4/4iYUFpfTr1cGRw3rl7TnEBHpbJKZIEqA0VHr+cCOJurOJeb0krvvCP9uAhbx7v6JDlVUXMZxYweSnqYJgkSk50hmgigEJprZeDPLIkgC77kaycwmAQOBV6K2DTSzXuHyEOBkYG1s245Qsb+WN0qr1f8gIj1O0q5icvd6M7seWAikA/e4+xozuxUocvfGZHE58LC7R59+mgLcYWYRgiT24+irnzrS0q3lgPofRKTnSVqCAHD3BcCCmG03x6zfEqfdf4DpyYwtUYXF5WSmGzPyB6Q6FBGRDqU7qVtQVFzG0aNyyclKT3UoIiIdSgmiGQfrGlixrZKCsep/EJGeRwmiGau3V1LbEFH/g4j0SEoQzSgsDjuodQQhIj2QEkQzlmwpY0JeHwb37ZXqUEREOpwSRBMiEadoSzknjNXpJRHpmZQgmrBxdzUV++s0QZCI9FhKEE14p/9BHdQi0kMpQTShqLiMIX2zGDe4d6pDERFJCSWIJhRuKaNg7CDMNECfiPRMShBxlFYdZFvZAfU/iEiPpgQRR1HY/3CC+h9EpAdTgoijsLiMnMx0po7sn+pQRERSRgkijqItZcwcM4DMdL09ItJz6RswRnVNPWt3VGl4DRHp8ZQgYizbWk7Edf+DiIgSRIzC4nLSDGaO0QRBItKzKUHEWLKljCkj+tMvOzPVoYiIpJQSRJS6hgjLtlbo8lYREZQg3mXdzir21zboBjkREZQg3uXQBEE6ghARUYKIUlRcRv7AHIbnZqc6FBGRlFOCCLk7hcXl6n8QEQkpQYS2lu1nT3WN+h9EREJKEKFCDdAnIvIuShChouIycnMyOTKvb6pDERHpFJQgQoXFZRw/diBpaZogSEQElCAAeLu6ho2796n/QUQkihIEsGSL+h9ERGIlNUGY2RwzW29mG8zsxjjlt5nZ8vDxhplVRJVdZWZvho+rkhnnki3lZKWnMX1UbjKfRkSkS8lI1o7NLB24HTgLKAEKzWy+u69trOPuN0TV/yIwM1weBHwXKAAcWBK2LU9GrIXFZRyTn0t2Znoydi8i0iUl8whiFrDB3Te5ey3wMHBRM/UvB/4cLp8DPOvuZWFSeBaYk4wgD9Y1sGp7peZ/EBGJkcwEMQrYFrVeEm57DzMbC4wHnm9NWzO7zsyKzKxo9+7dbQqy6mAd500fwSkTh7SpvYhId5W0U0xAvOtFvYm6c4F57t7QmrbufidwJ0BBQUFT+27W0H7Z/HruzLY0FRHp1pJ5BFECjI5azwd2NFF3LodOL7W2rYiIJEEyE0QhMNHMxptZFkESmB9bycwmAQOBV6I2LwTONrOBZjYQODvcJiIiHSRpp5jcvd7Mrif4Yk8H7nH3NWZ2K1Dk7o3J4nLgYXf3qLZlZvZ9giQDcKu7lyUrVhERea9k9kHg7guABTHbbo5Zv6WJtvcA9yQtOBERaZbupBYRkbiUIEREJC4lCBERiUsJQkRE4rKoi4e6NDPbDWxJdRwtGALsSXUQCegqcULXiVVxtq+uEid0/ljHuntevIJukyC6AjMrcveCVMfRkq4SJ3SdWBVn++oqcULXijWWTjGJiEhcShAiIhKXEkTHujPVASSoq8QJXSdWxdm+ukqc0LVifRf1QYiISFw6ghARkbiUIEREJC4liHZmZqPN7AUzW2dma8zsy3HqzDazSjNbHj5ujrevDoi12MxWhTEUxSk3M/tfM9tgZivN7LgUxDgp6n1abmZVZvaVmDopez/N7B4ze8vMVkdtG2Rmz5rZm+HfgU20vSqs86aZXZWCOH9mZq+H/7Z/M7MBTbRt9nPSAXHeYmbbo/59z2ui7RwzWx9+Xm9MZpzNxPqXqDiLzWx5E2077D09LO6uRzs+gBHAceFyP+ANYGpMndnAk50g1mJgSDPl5wH/IJjh70TgtRTHmw7sIrixp1O8n8ApwHHA6qhtPwVuDJdvBH4Sp90gYFP4d2C4PLCD4zwbyAiXfxIvzkQ+Jx0Q5y3A1xP4bGwEJgBZwIrY/3cdEWtM+S+Am1P9nh7OQ0cQ7czdd7r70nB5L7COJubi7gIuAh7wwKvAADMbkcJ4zgA2ununuWPe3V8CYucquQi4P1y+H/hwnKbnAM+6e5m7lwPPAnM6Mk53f8bd68PVVwlmbkypJt7PRMwCNrj7JnevBR4m+HdImuZiNTMDPsq7Z8rscpQgksjMxgEzgdfiFL/fzFaY2T/MbFqHBnaIA8+Y2RIzuy5O+ShgW9R6CalNdrFT00brDO9no2HuvhOCHwzA0Dh1Ott7+0mCo8V4WvqcdITrw1Nh9zRxyq6zvZ8fBErd/c0myjvDe9oiJYgkMbO+wKPAV9y9KqZ4KcFpkhnAb4DHOzq+0MnufhxwLvAFMzslptzitEnJddHhtLUXAn+NU9xZ3s/W6Ezv7XeAeuChJqq09DlJtt8BRwDHAjsJTt3E6jTvZ+hymj96SPV7mhAliCQws0yC5PCQuz8WW+7uVe5eHS4vADLNbEgHh4m77wj/vgX8jeAwPVoJMDpqPR/Y0THRvce5wFJ3L40t6CzvZ5TSxlNx4d+34tTpFO9t2Dn+IeAKD0+Ox0rgc5JU7l7q7g3uHgHuauL5O8X7CWBmGcAlwF+aqpPq9zRRShDtLDz3+Adgnbv/sok6w8N6mNksgn+HtzsuSjCzPmbWr3GZoMNydUy1+cAnwquZTgQqG0+dpECTv8g6w/sZYz7QeFXSVcDf49RZCJxtZgPDUyZnh9s6jJnNAb4FXOju+5uok8jnJKli+r0ubuL5C4GJZjY+PNqcS/DvkApnAq+7e0m8ws7wniYs1b3k3e0BfIDg0HYlsDx8nAd8FvhsWOd6YA3BlRavAielIM4J4fOvCGP5Trg9Ok4Dbie4OmQVUJCi97Q3wRd+btS2TvF+EiStnUAdwa/Ya4HBwD+BN8O/g8K6BcDdUW0/CWwIH9ekIM4NBOftGz+nvw/rjgQWNPc56eA4Hww/fysJvvRHxMYZrp9HcNXgxmTH2VSs4fb7Gj+bUXVT9p4ezkNDbYiISFw6xSQiInEpQYiISFxKECIiEpcShIiIxKUEISIicSlBSKdkZg3hSJerzeyvZta7hfrVHRVbzPMWmNn/tlBntpk9mcC+FplZp5jc3sz+nwUjEr8Q/RrD13JSquOTjpGR6gBEmnDA3Y8FMLOHCO57iHvjYSq5exHQqYZrNrMMPzQIX1tdC3ze3V8I1xtf42ygGvjPYe5fugAdQUhX8C/gSAAz+2p4VLHaYuaFCMsfNLOLotYfMrMLzexqM3vMzJ62YP6Fn0bVuTwcm3+1mf0kanu1mf0kHFDtOTObFf7K32RmF4Z13jk6CMv/Y2bLwr+TmntRZpZjZg+Hg9D9BciJKjvbzF4xs6XhEVTfcPt5Fszh8G8L5upofO5bzOxOM3sGeMDM0i2Y76Ew3P9novb9jajt34sT180EN3z+PtzHbDN70oLBJz8L3BAe3X2wudcn3UCq79TTQ494D6A6/JtBMFTF54DjCe6o7QP0JbgLdWZM/VOBx8PlXGBzuI+rCeZcyAWygS0EY/eMBLYCeWG954EPh+0dODdc/hvwDJAJzACWh9tnE85FAfTn0PwKZwKPxtaJeY1fBe4Jl48hGDCvABgCvAT0Ccu+Bdwcxr0NGB9u/3PUc98CLAFywvXrgJvC5V4ERwDjCYZ1uJPgLvk04EnglDixLSK8cz7mNd5CC3Mz6NF9HjrFJJ1Vjh2ajetfBONbfQ74m7vvAzCzxwiGVV7W2MjdXzSz281sKMGAaY+6e304VNM/3b0ybLsWGEswLMYid98dbn+IYCKYx4Fa4Olw16uAGnevM7NVwLg4MecC95vZRILkktnCazwF+N8w7pVmtjLcfiIwFXg5jDsLeAWYDGxy981hvT8TJIJG8939QLh8NnCMmX0kKraJ4fazo96zvuH2l1qIVXogJQjprN7pg2jUOCBfAh4EriAYsO2TUdtropYbCD7/ze2zzt0bx6KJNLZ394gFI3bG+j7wgrtfHJ6OWZRArPHGujGCyYQuf9dGs5kt7GtfzD6+6O7vGgDQzM4BfuTudyQQm/Rw6oOQruQl4MNm1jscBfNigqOLWPcBXwFw9zUt7PM14FQzG2Jm6QSjxr7Yxvhyge3h8tUJ1H+JIJFhZkcTnGaCYMDBk82ssd+lt5kdBbwOTAiTD8Blzex7IfA5C4aex8yOCt+zhcAno/o0RoVHW4naSzCVrvQAShDSZXgwlet9wGKCL/a73X1ZnHqlBFO93pvAPncC3wZeIBhdc6m7xxueOxE/BX5kZi8TzJHckt8BfcNTS98keF2Ep7uuBv4clr0KTA5PH30eeNrM/g2UApVN7PtuYC2w1MxWA3cQ9I88A/wJeCU8VTaP1n3hPwFcrE7qnkGjuUq3E94zsQo4rrHPobsws77uXh2ebrsdeNPdb0t1XNI96QhCuhUzO5PgVMxvultyCH067LxfQ3BKS30JkjQ6ghARkbh0BCEiInEpQYiISFxKECIiEpcShIiIxKUEISIicf1/NpRj03ZHQnwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "x, y = make_regression(n_samples=20, n_features=1, noise=10,effective_rank=1, bias=5)\n",
    "x = pd.DataFrame(x)\n",
    "# Create polynomial model\n",
    "\n",
    "rsq = np.zeros(20)\n",
    "degree = np.arange(1,20)\n",
    "for i in degree:\n",
    "    x['x'+str(i)] = x.iloc[:,[0]]**i\n",
    "    results = sm.OLS(y, x.iloc[:, 1:]).fit()\n",
    "    rsq[i] = results.rsquared\n",
    "\n",
    "    \n",
    "print('Condition number: {}'.format(la.cond(x,2)))\n",
    "ax = plt.plot(degree, rsq[1:])\n",
    "plt.ylabel('R-squared')\n",
    "plt.xlabel('Polynomial degree fit')\n",
    "#Find the polynomial when r-squared decreases\n",
    "for idx,i in enumerate(np.diff(rsq[1:], axis=0)):\n",
    "    if i<0:\n",
    "        plt.axvline(idx, color='red')\n",
    "        plt.title('Model is best fit @ '+str(idx))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.3 Spare recovery problems\n",
    "\n",
    "**Used to control over-fitting**\n",
    "\n",
    "\n",
    "Usually used when matrix is fat (M>N)\n",
    "\n",
    "$$Min ||x||_0 \\ \\ s.j.t \\ \\ AX=b$$\n",
    "- $||x||_0$ = sparse solution\n",
    "- AX = Fat matrix (underdetermined)\n",
    "- b = dense measurements\n",
    "\n",
    "\n",
    "The problem is that L-0 is not convex but we can find a solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3.1 Convex relaxation\n",
    "\n",
    "Keep trying larger norms until it touches the plane\n",
    "\n",
    "<img src='Images/l2_6.4_Convex_relaxation.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.3.1.1 Python example: Convex relaxation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.4 Sparse optimization problems\n",
    "\n",
    "`Basic Pursuit:`\n",
    "- Min 1-L solution\n",
    "$$Min |x| \\ \\ s.j.t \\ \\ Ax=b$$\n",
    "\n",
    "`Basic Pursuit Denoising:`\n",
    "- 2-L with L-1 penalty\n",
    "- Sometimes also called Lasso\n",
    "\n",
    "$$Min \\lambda|x|+ \\frac{1}{2} ||Ax=b||^2$$\n",
    "\n",
    "`Lasso:`\n",
    "- 2-L with \n",
    "$$Min \\frac{1}{2} ||Ax=b||^2\\ \\ s.j.t \\ \\ |x|\\leq \\lambda$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4.1 Sparse optimization with Bayes\n",
    "\n",
    "Basic Pursuit Denoising:\n",
    "\n",
    "$$\\lambda|x|+ \\frac{1}{2} ||Ax=b||^2$$\n",
    "\n",
    "\n",
    "\n",
    "- The smaller the $\\lambda$ the denser the solution\n",
    "- The larger the $\\lambda$ the sparser the solution\n",
    "\n",
    "\n",
    "\n",
    "Prior is often Laplace distribution which is more robust to outliers\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "`Laplace distribution:`\n",
    "\n",
    "$$P(x) = e^{-|x|}$$\n",
    "\n",
    "`Laplace PDF:`\n",
    "$$f(x; \\mu, \\lambda)= \\frac{1}{2\\lambda}\\exp(-\\frac{|x - \\mu|}{\\lambda})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6.4.1. Python example: Laplace vs. Gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXzcdb3v8ddnsk32ttmarUm6723atIWyF8QqUAQExAMWjooLBZWjiJ4DB7j3ejz1KnqvXBWRo6LIYacgUChQFqEb3fc1bbZmbfZ95nv/+M0kkzRtJ8kkv8zk83zQx2R+85uZT2j6nm8+3+/v9xNjDEoppYKfw+4ClFJKBYYGulJKhQgNdKWUChEa6EopFSI00JVSKkSE2/XGycnJJjc31663V0qpoPTpp59WGWNS+nrMtkDPzc1ly5Ytdr29UkoFJRE5fqbHtOWilFIhQgNdKaVChAa6UkqFCNt66Eqp0NPR0UFxcTGtra12lxL0nE4nWVlZRERE+P0cDXSlVMAUFxcTHx9Pbm4uImJ3OUHLGEN1dTXFxcXk5eX5/Ty/Wi4islxEDojIYRG5v4/HbxeRShHZ7vnztX7UrpQKEa2trSQlJWmYD5KIkJSU1O/fdM45QheRMOAx4DNAMbBZRNYYY/b22vW/jTGr+vXuSqmQo2EeGAP5/+jPCH0xcNgYc9QY0w48A1zb73dSKggYY3h2cxGtHS67S1Gq3/wJ9EygyOd+sWdbbzeIyE4ReV5Esvt6IRG5U0S2iMiWysrKAZSr1NB6Z18F972wk0fXHbS7FDVAcXFxg36NwsJCZs+ePaDnPvjgg6xbt27QNQyEP4He17i/91UxXgVyjTFzgXXAn/p6IWPM48aYAmNMQUpKn0euKmWrupYOAMpqdZWGGphHHnmEK664wpb39ifQiwHfEXcWUOq7gzGm2hjT5rn7e2BhYMpTani1eFoteh2v0PLqq6+yZMkS8vPzueKKKygvLwfgoYce4rbbbmPZsmVMmTKF3//+96c9t7CwkIsuuogFCxawYMECPv74467HVq9ezZw5c5g3bx7332+tF7n99tt5/vnnASvcFy1axOzZs7nzzjvxXiHu0ksv5Yc//CGLFy9m6tSpfPjhhwH5Pv1ZtrgZmCIieUAJ8CXgy747iEi6MabMc3cFsC8g1Sk1zCoarHGJy+22uZLg9/Cre9hbWh/Q15yZkcC/XzOr38+78MIL2bBhAyLCE088werVq/n5z38OwM6dO9mwYQNNTU3k5+dz1VVX9Xhuamoqb7/9Nk6nk0OHDnHLLbewZcsW3njjDV5++WU2btxITEwMNTU1p73vqlWrePDBBwG47bbbeO2117jmmmsA6OzsZNOmTbz++us8/PDDAWnTnDPQjTGdIrIKWAuEAU8aY/aIyCPAFmPMGuAeEVkBdAI1wO2DrkwpG1Q2WK2WqsZ2mytRgVRcXMzNN99MWVkZ7e3tPdZ2X3vttURHRxMdHc1ll13Gpk2bmD9/ftfjHR0drFq1iu3btxMWFsbBg9b8yrp167jjjjuIiYkBYNy4cae973vvvcfq1atpbm6mpqaGWbNmdQX69ddfD8DChQspLCwMyPfp14FFxpjXgdd7bXvQ5+sfAT8KSEVK2aiivs1zqz30wRrISHqo3H333dx7772sWLGC9evX89BDD3U91nt5YO/7jz76KGlpaezYsQO3243T6QSsFVFnW1rY2trKt7/9bbZs2UJ2djYPPfRQj3XlUVFRAISFhdHZ2TnYbxHQc7ko1UO5Z4ReXt/W1e9Uwa+uro7MTGtx3p/+1HPNxiuvvEJrayvV1dWsX7+eRYsWnfbc9PR0HA4HTz31FC6XNc9y5ZVX8uSTT9Lc3AxwWsvFG97Jyck0NjZ29dWHkh76r5SPk3XWCL2lw0VDWycJTv/Po6FGhubmZrKysrru33vvvTz00EPceOONZGZmct5553Hs2LGuxxcvXsxVV13FiRMneOCBB8jIyOjRAvn2t7/NDTfcwHPPPcdll11GbGwsAMuXL2f79u0UFBQQGRnJ5z//eX7yk590PW/MmDF8/etfZ86cOeTm5p72QTEUxK5RSEFBgdELXKiRpMPlZuq/vUFeUixHq5pYd+/FTE6Nt7usoLJv3z5mzJhhdxl+e+ihh4iLi+P73/++3aX0qa//nyLyqTGmoK/9teWilEdFQxvGwNysRMBquygVTLTlopRHuWcidE7WGF7eXsrJOp0YDXW+k6OhQEfoSnmUewJ8nneE3qCBroKLBrpSHic9I/S85Fjio8K7ljAqFSw00JXyOFnfSmSYg3GxkaQlOrtaMEoFCw10pTzK61pJTYhCREhLiNJAV0FHA10pj5P1rYxPsI4CTIt36iqXIFZeXs6Xv/xlJk6cyMKFCzn//PN56aWXhvQ9t2zZwj333DOk73EuuspFKY+K+jZmpCcAkJrgpKKh9ZyHd6uRxxjDF77wBVauXMnTTz8NwPHjx1mzZs2Qvm9BQQEFBX0uDx82OkJXCisETta3kuYdoSdE0eEynGrusLky1V/vvvsukZGRfPOb3+zalpOTw913333GU+GuX7+eq6++umv/VatW8cc//hGA+++/n5kzZzJ37tyuA5Cee+45Zs+ezbx587j44otPe41NmzaxdOlS8vPzWbp0KQcOHADgj3/8I9dffz3Lly9nypQp3HfffQH93nWErhTQ0NZJc7uL8YnWCZO8wV5e38q42Eg7Swteb9wPJ3cF9jXHz4HP/fSsu+zZs4cFCxb0+diZToV7JjU1Nbz00kvs378fEaG2thawznO+du1aMjMzu7b5mj59Oh988AHh4eGsW7eOH//4x7zwwgsAbN++nW3bthEVFcW0adO4++67yc7u8yJv/aaBrhTda9B9R+hgBbq3DaOC01133cVHH31EZGQk69at6/NUuGeSkJCA0+nka1/7GldddVXXCPyCCy7g9ttv56abbuo6Da6vuro6Vq5cyaFDhxAROjq6f9O7/PLLSUy0jnWYOXMmx48f10BXKpC8a9C9k6Kp8d0jdDVA5xhJD5VZs2Z1jYYBHnvsMaqqqigoKDjjqXDDw8Nx+1zUxHumxPDwcDZt2sQ777zDM888w69//Wveffddfvvb37Jx40b+/ve/M3/+fLZv396jhgceeIDLLruMl156icLCQi699NKux7ynzYXAnjoXtIeuFEDXYf7eEXqqZ4Repof/B51ly5bR2trKb37zm65t3lPcnulUuDk5Oezdu5e2tjbq6up45513AGhsbKSuro7Pf/7z/PKXv+wK7iNHjrBkyRIeeeQRkpOTKSoq6lGD7+l6vb344aAjdKXovvTc+EQr0KPCw0iOi9KLRQchEeHll1/me9/7HqtXryYlJYXY2Fj+8z//kwULFvR5Ktzs7Gxuuukm5s6dy5QpU8jPzwegoaGBa6+9ltZWa8XTo48+CsAPfvADDh06hDGGyy+/nHnz5vH+++931XDfffexcuVKfvGLX7Bs2bLh+9719LlKwQMv72bNjlJ2/PuVXdtW/PojEqMjeOqrS2ysLLgE2+lzRzo9fa5SA+B7UJFXeqJTWy4qqGigK4U1+ZmW2DvQoymrbdFL0amgoYGuFNbk5/iEqB7bMsY4aWp3Ud8auFUIo4F+AAbGQP4/aqCrUa+t00VlQxuZY2J6bE9PjAagrK7FjrKCktPppLq6WkN9kIwxVFdXdy2r9JeuclGjnnfJYsaYnv94vPfLaluZPl4PLvJHVlYWxcXFVFZW2l1K0HM6nT0udu0PDXQ16pXUWiPwzDHRPbZ7R+ilOkL3W0REBHl5eXaXMWppy0WNeqW13hF6z0BPjY/CIehadBU0NNDVqFfqGaGP77XKJTzMQVqCU0foKmhooKtRr7S2heS4KJwRYac9lp7o1BG6Choa6GrUK6ltIXNM36sJMsZE6whdBQ0NdDXqlda2nNY/98oYE01ZXasuw1NBQQNdjWrGGEprW88Y6OmJTto73VQ3tQ9zZUr1nwa6GtVqmzto6XCdJdA9BxdpH10FAb8CXUSWi8gBETksIvefZb8viogREXuvlKqUn7rXoJ+ph25t1z66CgbnDHQRCQMeAz4HzARuEZGZfewXD9wDbAx0kUoNlbK6vtege3WP0DXQ1cjnzwh9MXDYGHPUGNMOPANc28d+/wNYDejvpipoeNegZ4yJBlcn7HkJ1j0M2/4C7c0kxUYSGe6gVE+jq4KAP4f+ZwK+11cqBnqc8V9E8oFsY8xrIvL9M72QiNwJ3AkwYcKE/lerVICV1rYQGe4giXr4w01QuhXEAcYNH/xvHF9+lqwx0RSfara7VKXOyZ8RuvSxrWsNl4g4gEeBfznXCxljHjfGFBhjClJSUvyvUqkhUlLbwsREB/LUdVCxD274AzxQBV95Bdqb4E/XMCOhleJT2nJRI58/gV4MZPvczwJKfe7HA7OB9SJSCJwHrNGJURUMSmtb+BeegvJdcPNTMOeL4AiDiZfCbS9Bay3fqf85RdVNdpeq1Dn5E+ibgSkikicikcCXgDXeB40xdcaYZGNMrjEmF9gArDDG6AVD1YgXX7OHzzS9Cud9G6Z8pueD42fDlf+TqY2bOK/tHzS26YUu1Mh2zkA3xnQCq4C1wD7gWWPMHhF5RERWDHWBSg2Vtk4XX2v/My3hiXDpGVbjFvwzdQnT+FH40xRX1w9vgUr1k1/r0I0xrxtjphpjJhlj/pdn24PGmDV97Hupjs5VMKja9yEXOXZxaOqd4EzseydHGNWL/oUJjko6dr40vAUq1U96pKgatSK2/J56E0NH/sqz7pcwfwWH3Jlk7v4t6Dld1Aimga5Gp/oykk+8ybOuS8hMPfuKq6Q4J09xFeMaD0Lx5mEqUKn+00BXo9O2p3CYTv7GZ0mNjzrrriLC9sRltInTOuBIqRFKA12NPsbAruc46JyHGZuHw9HXoRY9JY1L4oOIC2H3i9b6dKVGIA10Nfqc3AVVB3kr7EKyx8b49ZTscTE83X4htDfA/r8PcYFKDYwGuhp9dj0HjnCebVrIhHH+BXrW2GjWt07GHZcOe18Z4gKVGhgNdDW6GAO7X6Qj7zJOtDrJHtf3WRZ7yx4bg8HBqZzPwuF10NY4xIUq1X8a6Gp0KdsB9cWUZy0H8LvlkuXZ70jyMuhstUJdqRFGA12NLofeAoSD8dYJQ7P9bLl4R/K7wmZCTDLse3WoKlRqwDTQ1ehycC1kLuBwsxXk/gZ6YnQEcVHhFNW2w5Qr4cg74HYNZaVK9ZsGuho9Giuh5FOY8lmKalpIjI4gMTrCr6eKCBPGxXCiphkmXw4tp6Bk6xAXrFT/aKCr0ePw24CBqVdyoqbZ7wlRr9zkGAqrmmDSMusiGNpHVyOMBroaPQ69BXFpMH4eRaea/Z4Q9cpJiqXoVDOdUWMgc6EGuhpxNNDV6OB2w7EPYNIy3AjFp1r8XoPulZsUQ4fLWBeWnnyF1b5pqh6igpXqPw10NTpU7oPmasi9iJP1rbR3uv2eEPXKTYoFoLDa03bBQOGHQ1CsUgOjga5Gh2MfWLd5F1l9cGBicmy/XiLXs39hVRNk5ENELBR+FNAylRoMDXQ1Ohz7EMbmwpgJHPUEem4/Az01PgpnhIPC6mYIi4AJ52mgqxFFA12FPrcLjn8EuRcB1gjbGeFgfIKzXy8jIuQmxXLce8Ho3AusVk5TVaArVmpANNBV6Du5E1rrIO8SwOqB5ybF+nXa3N5ykmKsETp0fUBw/B+BqlSpQdFAV6HvmGfiMs8K4KNVTeT1s93ilZscy4nqZlxu4+mjx2jbRY0YGugq9B3/GJImQ/x4Ol1uimqa+90/98pNiqXd5aasrkX76GrE0UBXoc0YKNoI2ecBUFLbQofLkJc0sEDPSbKWOh7vartcCBV7dT26GhE00FVoqz4CLTWQvRiAY54VLnkpAx+hg2ctOsCE863b4k2Dq1OpANBAV6GtaIN1O8EaoXvXoOcOcIQ+PsFJZLij63XIyAdHuPVbgFI200BXoa1oIzjHQNIUwBqhx0WFkxwXOaCXcziEnHExHKvytFwiomH8XCjaHKiKlRowDXQV2oo2We0Wh/Wjfqy6mbzkWET6v2TRa1JKHEerfC5Bl70YSreCq2Ow1So1KBroKnQ110Dlfshe0rWpsKppwCtcvCalxnK8upn2Tre1IWsRdDRD+e5Bva5Sg6WBrkJX8Rbr1hPo7Z1uik81k5fUv5Ny9TY5NQ6X23CipqnH62vbRdlNA12FrqKNIGGQuQCAEzVNuM3AV7h4TUqJA+BwhSfQE7MgPl0nRpXtNNBV6CraCOlzIdIK8MMVVt97Smr8oF7WG+hHKj19dBGr7aJLF5XN/Ap0EVkuIgdE5LCI3N/H498UkV0isl1EPhKRmYEvVal+cLusa35mLeradKjcCuCJgxyhx0aFk57o5EiF78ToEqg9AQ0nB/XaSg3GOQNdRMKAx4DPATOBW/oI7KeNMXOMMfOB1cAvAl6pUv1RdQg6miBjQdemQxWNZI2NJiYyfNAvPzk1rnuEDl0HLlGko3RlH39G6IuBw8aYo8aYduAZ4FrfHYwx9T53YwETuBKVGoDSbdZtRn7XpkMVjUxOjQvIy09KieNIZRPGeH7U0+eBI8K6LJ1SNvEn0DOBIp/7xZ5tPYjIXSJyBGuEfk9fLyQid4rIFhHZUllZOZB6lfJP6TaIjINk64Ail9twtLKRKYEK9NQ4Gts6Ka9vszaER0HaLGs9ulI28SfQ+zoC47QRuDHmMWPMJOCHwL/19ULGmMeNMQXGmIKUlJT+VapUf5Ru84yawwAoPtVMW6d70BOiXpNSek60AtZvA6XbrQtSK2UDfwK9GMj2uZ8FlJ5l/2eALwymKKUGxdVpXdTCt93imRCdFKAR+uTeK13AWh7ZVg81RwLyHkr1lz+BvhmYIiJ5IhIJfAlY47uDiEzxuXsVcChwJSrVT5X7obO1R6Af9gRvoHroKfFRxDvDewa6dwK2RNsuyh7nDHRjTCewClgL7AOeNcbsEZFHRGSFZ7dVIrJHRLYD9wIrh6xipc6lrwnR8kZS46NIjI4IyFuICJNS4nq2XFKmQ3h09/srNcz8Wr9ljHkdeL3Xtgd9vv5OgOtSauBKt0JUIozN69p0uKKBKWmBGZ17TU2L4939Fd0bwsKtvr1OjCqb6JGiKvSUboOMeV1nWDTGcLiiMWATol7TxidQ1dhOVWNb98bMBVC20+rjKzXMNNBVaOlsg5O7e7RbyupaaWp3BWxC1Gv6eOsD4sDJhu6NGfnQ2QKV+wL6Xkr5QwNdhZaKveDu6BHoB8qtwA3UGnSvaZ5A398j0HViVNlHA12Flj4mRPeXWYE7Y3xCQN8qOS6K5LhI9pf5HCg9bqLVv9c+urKBBroKLaXbIHosjMnp2rSvrJ6MRCeJMYFZ4eJr2vj4rt8AAKtvnzFfR+jKFhroKrSUbrNG5z6XmNt/sp4Z6YEdnXtNS0vgYHkDLrfPwdOZC6zWT0frkLynUmeiga5CR0cLVOzr0W5p7XBxpLJpyAJ9+vh4WjvcnKhp7t6YsQDcnXpJOjXsNNBV6CjfYwWp7xGiFY243Ibp6YFdsug1rWuli08f3fv+2nZRw0wDXYWOPiZE93kmLIdqhD41LR6RXitdErMgNkWPGFXDTgNdhY7SbVaQJnSf3Xn/yQacEQ5ykwZ3laIziY4MIzcptmslDWD17zPyNdDVsNNAV6GjjwnRfWX1TEuLJ8zR11mgA2NaWjz7fVsuYPXRqw5AW2PfT1JqCGigq9DQ3mSdZdGn3WKMYV9ZPdMDvP68t5kZCRRWN9PQ2tG9MSMfjBtO7hrS91bKlwa6Cg1lO60A9Qn0yoY2TjV3MGOIJkS95mQmArC31HdidL51q20XNYw00FVo8AZn+vyuTXs9E6LTh2hC1GtWpvX6u30DPX48xGfoEaNqWGmgq9BQug3i0yEhvWvT7pI6wGqJDKXUeCep8VHs8bxfF50YVcNMA12FBu+EqI+dxXVMTI4lwRn4Q/57m5OZyK6+Ar36MLTW9f0kpQJMA10Fv9Z6qD50WqDvKqljTlbisJQwKzORI5WNNLf7nAfdW0/ZjmGpQSkNdBX8vIHpPXUt1oRoWV1r14TlUJuTmYjbwL6yXudGB227qGGjga6CX9cRot0Tot7++dysMcNSwmzvxKhv2yU2CcZM0EBXw0YDXQW/0m2QOAFik7s27SyuQwRmDfGEqNf4BCdJsZE9Ax10YlQNKw10FfxKt/UYnQPsKqllckocsVF+XQd90ESE2WeaGD1VCM01w1KHGt000FVwazkFp471ucJluPrnXnMyEzlU0UhLu6t7o/bR1TDSQFfBrXS7desT6OX1rVQ0tA3bChev+dljcLlNz1F6+jzrVgNdDQMNdBXc+pgQ3VnsnRAd5kCfYE3Abjtxqntj9FjrOqMa6GoYaKCr4Fa6DcbmWcHpsb3oFOEOYWb68AZ6clwUOUkxbDtR2/OBjAXdv0koNYQ00FVw62NC9NPjp5iVkUB0ZNiwl5OfPYatJ05hjM81RjPyob4YGiuGvR41umigq+DVWAl1RT0OKOpwudlRVMeCnLFneeLQyZ8wlgrPQU1duiZGdZSuhpYGugpe3r50Zneg7y9roKXDxULbAt3qo2/17aOnzwVE++hqyGmgq+BVuhWQ7pUkwKfHrfXeCybYE+gz0hOICnf07KNHxUPyVA10NeQ00FXwKt1mBWVU9wUsPj1RS3qik4wx0baUFBHmYG5WYs+VLuA5YnQr+PbWlQowDXQVnIyBkq092i0AW4+fsq1/7pU/YSy7S+pp6/Q5wChzATSWQ0OZfYWpkOdXoIvIchE5ICKHReT+Ph6/V0T2ishOEXlHRHICX6pSPupLoamixwFFZXUtlNS2sNCmdovXwpyxtLvcXevhAT1iVA2Lcwa6iIQBjwGfA2YCt4jIzF67bQMKjDFzgeeB1YEuVKkevJd281nhsvW41be2a0LUa1HuOAA2HfM5f0vabJAwDXQ1pPwZoS8GDhtjjhpj2oFngGt9dzDGvGeMafbc3QBkBbZMpXop3QaOcBg/u2vT5sIanBEOZgzxNUTPZVxsJNPS4tlwtLp7Y2QMpM7QQFdDyp9AzwSKfO4Xe7adyVeBN/p6QETuFJEtIrKlsrLS/yqV6q1kqxWQEd2TnxuOVlOQM47IcPunhpZMHMenx0/R4XJ3b8yYbwW6ToyqIeLPT770sa3Pn0gRuRUoAH7W1+PGmMeNMQXGmIKUlBT/q1TKlzGeI0S72y01Te3sP9nA+ZOSbCys25K8JJrbXT3Pj56RD83V1sFQSg0BfwK9GMj2uZ8FlPbeSUSuAP4VWGGMaQtMeUr14dQxaK3tscJlo6e9cd7EkRHoi/KsPv5G3z66ToyqIeZPoG8GpohInohEAl8C1vjuICL5wO+wwlxPWKGGVol3QrR7hcsnR6uJiQwb9jMsnklqvJOJKbGnT4w6IjTQ1ZA5Z6AbYzqBVcBaYB/wrDFmj4g8IiIrPLv9DIgDnhOR7SKy5gwvp9TglW6DcCekdi+2+uRINQW544gIs79/7rUkL4nNx2pwuT0dyvAoSJvV/YGkVID5dX0uY8zrwOu9tj3o8/UVAa5LqTMr3Qbj50BYBABVjW0cqmjkugVnm6sffudNHMffNp1gT2ld98WqM/Jh94vWPID0NT2l1MCNnOGMUv5wu6yzFvq0W7zLA88fIf1zr6WTrItWf3ioqntjRj601UHNUZuqUqFMA10Fl6pD0NHUY4XLx0eqiY0MG/ZriJ5LSnwUM9MT+PCQzxJdnRhVQ0gDXQUX7xGinhUuxhg+OFjJ0snJhI+g/rnXRVOT+fT4KZraOq0NqTMgLEoDXQ2JkfcvQKmzKdoEUYmQNAWAo1VNFJ9q4ZKpI/O4hounpNDhMt1HjYZFWP1/vdiFGgIa6Cq4FG+GrIXgsH503z9gtTNGaqAvzBmLM8LRs4+euQDKtoPbfeYnKjUAGugqeLQ1QMVeyFrcten9g5VMTIkle1yMjYWdmTMijPMmJvFB7z56eyNUH7KvMBWSNNBV8CjZCsYN2YsAaO1wseFo9YgdnXtdNCWFo5VNFJ/ynL8uc6F1W7zFvqJUSNJAV8GjeJN16wnEjcdqaOt0j/hAv2SqtXzxPU97iKQp4BwDRRttrEqFIg10FTyKNkPyNIi2zpOy/kAFUeGOEXP+ljOZlBJHblIM6/aWWxscDshebE3wKhVAGugqOBhjTYh62i3GGNbtK+f8SUk4I8JsLu7sRIQrZqTxyZFqGr3LF7MXQ+U+aKk9+5OV6gcNdBUcao5CSw1kWYG+/2QDRTUtfHbWeJsL888VM9Nod7n58KCn7ZK9xLrVProKIA10FRy87QnPCpe1e04iAlfMSLOxKP8V5IxlTEwEb+/ztF0yFliXpNM+ugogDXQVHIo3Q1QCpEwHYO2echZOGEtKfJTNhfknPMzBsmmpvLe/gk6XG6LirMvnaaCrANJAV8GheJN1QI7DQVFNM/vK6oOm3eJ1xcw0TjV38OnxU9aG7CVQ8im4Ou0tTIUMDXQ18rU1QPmeHu0WgCtnBUe7xeviqSlEhTt4fVeZtSF7iXWAUcVeewtTIUMDXY18RZusA4pyzgesQJ8+Pp6cpFibC+ufuKhwlk1P5e+7TloXvcj2HPGqbRcVIBroauQ7/rE1gZi1mJLaFjYXnuKqOel2VzUgV8/NoKqxzboGamI2xKfrenQVMBroauQ7/jFkzIeoOF7bYV2ffMX8DJuLGphl01OJiQzj1Z1l1hWLshfDiQ12l6VChAa6Gtk6WqFkC0yw2i2vbC9lfvaYoGu3eEVHhnHFjDTe3F1Gh8sNORdA3QmoPWF3aSoEaKCrka3kU3C1Q84FHCpvYG9ZPdcG6ejc6+q56Zxq7uAfh6sg90JrY+E/7C1KhQQNdDWyHf/Yup1wHmt2lOIQuGpucPbPvS6ZlkJidAQvbi2BlBkQPQ4KP7K7LBUCNNDVyHb8H5A6C7dzLC9vL2HppGRS4512VzUoUeFhXDs/g7V7TlLX5oLcC6DwQ7vLUiFAA12NXK4OawVIzlI2HK2mqKaFLy7MsruqgLhxYTZtnW5e3SrxYuwAABXcSURBVFEKuRdB7XHto6tB00BXI1fZTuhogpyl/G1zEYnRESyfHVxHh57J7MwEpo+P57ktRdpHVwGjga5GrmPvA1Cbsoi1u09yXX7miD9Vrr9EhBsLstlRXMdBk6V9dBUQGuhq5Dr6HqTO4oVDnbS73Ny8KNvuigLqC/MziAgTntlcon10FRAa6Gpkam+GExswEy/lmU0nmJc9hhnpCXZXFVBJcVEsn53Oc58W0Za1VPvoatA00NXIdOITcLWzL3oBhyoa+afFE+yuaEjcvjSHhtZO3mqeZm048p69BamgpoGuRqaj68ERwa+PpjIuNjJoD/U/lwUTxjI7M4H/uysME58BR96xuyQVxDTQ1ch0dD0t4wt442ADty6ZEDKTob2JCCvPz+VgRRMVqRfAkfV6fnQ1YH4FuogsF5EDInJYRO7v4/GLRWSriHSKyBcDX6YaVZqq4OROPjZzCHcIt56XY3dFQ+qaeRmMjYngxYbp0FZnne5AqQE4Z6CLSBjwGPA5YCZwi4jM7LXbCeB24OlAF6hGoaPrAfh9SQ7XzMsgNSG4jww9F2dEGLcvzeM3J7Ix4tC2ixowf0boi4HDxpijxph24BngWt8djDGFxpidgHsIalSjzaG3aQlPZFN7Dl+/aKLd1QyLlUtzcEUmUhg1HQ6vs7scFaT8CfRMoMjnfrFnW7+JyJ0iskVEtlRWVg7kJVSoc7twH3qLdZ3zuGJmesgtVTyTMTGR3HpeDmsap2NKtkJzjd0lqSDkT6BLH9vMQN7MGPO4MabAGFOQkpIykJdQoa54M46WGt5sn8/dy6bYXc2w+uqFeXxEPoKBw9p2Uf3nT6AXA76H6GUBpUNTjhrt2ve9TgdhmEnLmJOVaHc5wyo1wcmMhZdQaRJp2rnG7nJUEPIn0DcDU0QkT0QigS8B+tOmhkT9jlfZ5JrOnVfm212KLe66fBrvmYWEHX0bOtvsLkcFmXMGujGmE1gFrAX2Ac8aY/aIyCMisgJARBaJSDFwI/A7EdkzlEWr0FRVfIDk5qOUpV3M/Owxdpdji7QEJ46Z1+B0t1C45Q27y1FBxq916MaY140xU40xk4wx/8uz7UFjzBrP15uNMVnGmFhjTJIxZtZQFq1C08bXnwJgyfJbba7EXldefRNNODn60bN2l6KCjB4pqkaEwxUNpBe/SWn0VLInz7a7HFslxMVRnnoxcxo+Yv2+MrvLUUFEA13ZzhjDL59/jwWOQyQW3Gh3OSNC9tIbSZE6nn/lJVo7XHaXo4KEBrqy3SvbS0kreROA2Pk32FzNyBAx/bO4HZEsbFzP4x8ctbscFSQ00JWt6lo6+J9/38tN0Vsw4+dC0iS7SxoZnIk4pi3nhqhN/Pa9A5yobra7IhUENNCVrX76xj6cTSVM6zyAzLrO7nJGlrk3keA6xYWO3fzwhZ243QM6nk+NIhroyjbv7i/nb5uK+MkkzyrX2dfbW9BIM+VKcCby46xdfHK0mj99Umh3RWqE00BXtqhpaue+53cxIy2WixrfgtyLYGyu3WWNLOFRMPNacirfY/nUeH76xn6OVDbaXZUawTTQ1bAzxvCjF3dS19LOby5qQ2oLIX90rz0/o7k3Ix1NrJ5ZSHRkGN99ZjttnbrqRfVNA10Nuyc+PMbaPeXc99np5Ba/ApHxMGOF3WWNTDkXQNJkEnY/xeob5rKrpI7/8dpeu6tSI5QGuhpWG49W89M397N81ni+VjAW9rwEs6+DyBi7SxuZRKDgn6F4E1cmVfKNiyfylw0neHlbid2VqRFIA10Nm5LaFu56ehs5STH87Ma5yPa/QEczLPq63aWNbPNugXAnbHmSH3x2GovzxnH/izvZWVxrd2VqhNFAV8OirqWDO/5rE22dLn5360LiIx2w6XGrpZA+1+7yRraYcTDretj5LOEdjTz25QUkx0Xxz3/cQlGNrk9X3TTQ1ZBr63Txjae2cKyqid/dtpApafFw8E2oPQFLvmF3ecFh8degvRG2/pmU+Cj+eMciOlxubv+vTdQ2t9tdnRohNNDVkGrvdLPq6W1sOFrDz744j6WTksEY+Mf/gYQsmHaV3SUGh8yF1tLOT34NnW1MTo3n8dsWUlTTwsr/2kx9a4fdFaoRQANdDZkOl5tVT2/l7b3lPLxiFl/I91yKtvBDKNoAF34XwsLtLTKYXPg9aCiDnf8NwJKJSfy/f1rA3tI6vvKHTRrqSgNdDY3WDhd3/XUrb+0t59+vmcnKpbndD76/GuLGQ/5tttUXlCYtg/Fz4aNHwWWF9xUz03jsywvYXWKF+qkmbb+MZhroKuBONbVz6xMbeWtvOQ9dM5M7LsjrfrDwI2uEfuF3IcJpX5HBSAQuvR9qjsK2p7o2XzlrPI/90wL2ltVzw28+1onSUUwDXQVUYVUTN/z2Y3YW1/HrL+dzu2+Yu92w9seQkAkLVtpXZDCb9nnIPg/W/xTam7o2f3bWeP7y1SVUN7Vz3f/7mG0nTtlYpLKLBroKmLV7TnLN//2I6sZ2nvrqYq6em9Fzhx1PQ9kOuOJhPZBooETgM49AY7k1sexjcd44XvjW+TgjHNz0u0946pNCjNEzNI4mGuhq0No6XfzH6/v4xlOfkpscy2t3X8iSiUk9d2qphXcegaxFMOeL9hQaKiYssdalf/QoVB3u8dDk1HheXXUhF0xO5oFX9vC9/96uk6WjiAa6GpTdJXVc++t/8LsPjvLlJRN47pvnkz2uj9H32n+Fpir43GprlKkGZ/lPraNHX/uutQzUx9jYSJ5cuYh7PzOVNTtKWf7oB3xwsNKmQtVw0kBXA9LU1snqN/fzhcf+QU1TO39YWcBPrpuDMyLs9J0PvgXb/wIXfAcyFwx/saEoPg0+87A1wbz5idMedjiEey6fwvPfWkp0ZBhfeXITP3x+J9WNbTYUq4aL2NVjKygoMFu2bLHlvdXAGWN4ZXsp//HGPsrr27h+QSYPXj2TMTGRfT+hvhR+dwnEJME33rfO8a0Cw+2Gp2+CY+/D19454ykUWjtcPPr2QZ746BgxkWF894qpfOX8HCLCdDwXjETkU2NMQZ+PaaArfxhjeHtvOb965xB7SuuZm5XIv18zk4U54878pM42+ONVUL4Xvv4upE4fvoJHi6Yq+O2FEBFthXrMmf8+DpU38Mhre/nwUBW5STHcddlkvpCfqcEeZDTQ1YC1d7p5c89Jfrv+CHvL6slJiuGeZVO4Lj8Th+MsvXC3G17+pnVU441/gllfGL6iR5sTG+BP10DGAvjKy1a4n4ExhvcOVPDztw6yp7Se7HHRfPOSSVyXn0lMpB61Gww00FW/naxr5W+bTvD0phNUNrSRlxzLqssmc+38DMLPNaIzBl7/vtXbXfZvcPEPhqfo0Wz3i/D8HTD1c3DTn87Z2jLG8O7+Cv7PO4fYUVxHvDOcGxdmc+t5E5iYEjdMRauB0EBXfqlr6WDt7pO8vL2ET45WYwxcOi2FlefncsnUlLOPyL1cHbDmHmvN+dJ7rDXTuqpleGz6vfVBOvFSuPmvEHXuYDbGsOX4Kf78yXHe2FVGp9swP3sMK+ZlcPXcdFIT9GjekUYDXZ1RYVUT7x2oYP2BSj45Uk27y01OUgzXzs/k+vxMcpNj/X+xxgp44WvWJN2lP4JLfqhhPty2Pw2v3AUp0+Hmv0DSJL+fWtHQyotbS1izvZS9ZfU4BApyx3HZtFQunZbC9PHxiP592k4DXQHWaOxETTObC0/x6fEaNhyt4ViVdfj4xORYLpueyjXzMpiXldi/f7jGWOc3X3MPtNXDVT/Xiz7b6fA71gerq8Na2rjwDnD0b+LzcEUDa7aXsm5fBXvL6gFIT3Ry/sQkCnLHUZA7lskpcf791qYCSgN9FHK7DcdrmtlfVs++snr2ljWwvaiWKs865HhnOAU5Y7nUM/rKSerHSNzXyd2w7iE4/DakzoQb/gBpMwP3jaiBqS2Cl79lrVPPWgyXP2CdT30AI+zy+lbeP1DJ+oMVbDpWQ1WjdUbHxOgI5mWPYUZ6PDPTE5iRnsDE5Nhzz7GoQdFAD1GdLjcVDW0cr27mRE0Tx6ubOV7TzInqZg5XNNLS4QLAIZCXHMuczMSu0dXU1PiBj64626xR4OYn4Mg7EJVgnQVw8Z0QFhHA71ANijGw42+w7mFoPAnZS6zR+swVEDmwD3BjDMerm9lcWMOWwlPsKqnjUEUDHS4rRyLDHeQlxZKTFENuciy5SbHkJsWQPS6G1IQoosL7OPBM9cugA11ElgO/AsKAJ4wxP+31eBTwZ2AhUA3cbIwpPNtraqCfzuU2NLR2UNvcQV1L95/alg4qG9qoqG+loqGNcs9tdWMbbp+/vnCHkDU2mglJsUxMjmVmegLT0+OZmhbf9xGc/jIGqo/A8X9Yp789uBba6iAuzbqE3MI7zrr+WdmsowW2/hk2/tY69W5ELORdDJMvhwnnQ8q0QX0Qd7jcHKlsZG9pPftPNnC0sonj1U0cr2mmvdPdY9+xMRGkJTgZn+gkLd5JWkIUY2IiGRsbwZjoSMbERDA2xrpNcEZoS6cPgwp0EQkDDgKfAYqBzcAtxpi9Pvt8G5hrjPmmiHwJuM4Yc/PZXne4A90YgzHgNga357b7vrXN+DzmfbzTbeh0uelwuelwGc9t99edLkO759b3sU63m9YOFy3tbpo7Omltd9Hc7qK5w3Xa103tndS1dNDQ2nnG+kUgKTaKtIQoUuOjSPX8Y0hLdDJhXAy5SbGkJzr9/3XX7bL+oXc0W3/am6G1DpoqrMnNxgrrmp9VB6HqELQ3WM+LTYHJn4HZN8DES3REHkyMgROfwO4X4NDbUHvc2h4WBakzYNxEGJMNidnWh3X0GHCOAWciOBOsc8eERYLDv8GB220oq2+lsKqJktoWyutaOVnfSnm997aNqsa23qei6eIQiI0KJzYynJioMOs2MszaFhVObGQYMZHhOCMcRIZ7/oQ5iAr3vR/W47HIcAcRYYJDhDCHEO4QHA4hzHO/xx8RwsJ6PSZi+4fMYAP9fOAhY8xnPfd/BGCM+Q+ffdZ69vlERMKBk0CKOcuLDzTQP37uUcbvsc5dIRi8byDG+lo8W7q+NvTYbj2PrufT59fdRE7/Fvx6nme7eN5dBBwYRATB+mG1Xt86oY4IiEjXPg4xSNfzPI93fWceXf97fbf5VtrXvoCrHVznOqeHQHw6pEyF5GlWXzznAkiarCtXQoEx1mi9ZCuc3AEnd8Gp41BfYv18nI2EWcEeHukJ+HBAPD8XvreAOM7wmHVrsAZQLrf1p+fXPQdcPQZh7u5BmOG085MND+n5b96zqc870mtD5YLvsmTF1wf2tmcJdH8ODcsEinzuFwNLzrSPMaZTROqAJKCqVyF3AncCTJgwwa/ie0tIGk9j4hTrZ6UrWMSTMZ5b8fnfJ4J4/0d6wtS7r/WfdG3vvu1+LcH6RBYRwgREHF33HSI4BBziQBx47jtweL92COEOBw7fAOwRhufYPlT7OsKsHmpENETEeP5EW6OwuDSITbXOvaLX+wxdItaSxqRJMPfG7u1ut+e3tHLrNzbfP51t1soZV5tnUNBhbXN30jVywnjS1ffW3cc274DH6uMOtrPuxgp3l7vXh4DbnHbf0P0be9fXWCX1+BrPPr5fe/b3fKfd37ZH78d8bqzhp2f/hHEpg/yO++bPv9i+hmO9Pw/92QdjzOPA42CN0P1479PMXnYLLLtlIE9VSp2LwwHx460/QcTbaBztU67+NFyLgWyf+1lA6Zn28bRcEoGaQBSolFLKP/4E+mZgiojkiUgk8CVgTa991gDei0R+EXj3bP1zpZRSgXfOlounJ74KWIv1G82Txpg9IvIIsMUYswb4A/CUiBzGGpl/aSiLVkopdTq/Zr2MMa8Dr/fa9qDP163Ajb2fp5RSavjoMbpKKRUiNNCVUipEaKArpVSI0EBXSqkQYdvZFkWkEjg+wKcn0+so1BFC6+q/kVqb1tU/Wlf/DKauHGNMn4ea2hbogyEiW850LgM7aV39N1Jr07r6R+vqn6GqS1suSikVIjTQlVIqRARroD9udwFnoHX130itTevqH62rf4akrqDsoSullDpdsI7QlVJK9aKBrpRSISJoA11E5ovIBhHZLiJbRGSx3TV5icjdInJARPaIyGq76/ElIt8XESMiyXbXAiAiPxOR/SKyU0ReEpExNtez3PN3d1hE7rezFi8RyRaR90Rkn+dn6jt21+RLRMJEZJuIvGZ3LV4iMkZEnvf8bO3zXErTdiLyPc/f4W4R+ZuIOAP5+kEb6MBq4GFjzHzgQc9924nIZcC1WBfNngX8b5tL6iIi2VgX+z5hdy0+3gZmG2PmYl2M/Ed2FeK5IPpjwOeAmcAtIjLTrnp8dAL/YoyZAZwH3DVC6vL6DrDP7iJ6+RXwpjFmOjCPEVCfiGQC9wAFxpjZWKcjD+ipxoM50A2Q4Pk6kdOvomSXbwE/Nca0ARhjKmyux9ejwH30cXlAuxhj3jLGdHrubsC6IpZdFgOHjTFHjTHtwDNYH862MsaUGWO2er5uwAqnTHursohIFnAV8ITdtXiJSAJwMdZ1GjDGtBtjau2tqks4EO25slsMAc6tYA707wI/E5EirFGwbSO7XqYCF4nIRhF5X0QW2V0QgIisAEqMMTvsruUs/hl4w8b37+uC6CMiOL1EJBfIBzbaW0mXX2INEtx2F+JjIlAJ/JenFfSEiMTaXZQxpgQrq04AZUCdMeatQL7HiL6su4isA/q6Wu2/ApcD3zPGvCAiN2F9Gl8xAuoKB8Zi/Wq8CHhWRCYOxyX5zlHXj4Erh7qGvpytLmPMK559/hWrtfDX4aytF78udm4XEYkDXgC+a4ypHwH1XA1UGGM+FZFL7a7HRziwALjbGLNRRH4F3A88YGdRIjIW6ze+PKAWeE5EbjXG/CVQ7zGiA90Yc8aAFpE/Y/XuAJ5jGH/lO0dd3wJe9AT4JhFxY52Ip9KuukRkDtYP0Q4RAautsVVEFhtjTtpVl099K4GrgcttvhatPxdEt4WIRGCF+V+NMS/aXY/HBcAKEfk84AQSROQvxphbba6rGCg2xnh/i3keK9DtdgVwzBhTCSAiLwJLgYAFejC3XEqBSzxfLwMO2ViLr5ex6kFEpgKR2Hy2N2PMLmNMqjEm1xiTi/UDv2A4wvxcRGQ58ENghTGm2eZy/Lkg+rAT61P4D8A+Y8wv7K7HyxjzI2NMludn6ktYF4e3O8zx/FwXicg0z6bLgb02luR1AjhPRGI8f6eXE+DJ2hE9Qj+HrwO/8kwutAJ32lyP15PAkyKyG2gHVto86hzpfg1EAW97fnvYYIz5ph2FnOmC6HbU0ssFwG3ALhHZ7tn2Y8+1flXf7gb+6vlgPgrcYXM9eNo/zwNbsdqL2wjwKQD00H+llAoRwdxyUUop5UMDXSmlQoQGulJKhQgNdKWUChEa6EopFSI00JVSKkRooCulVIj4/7kInExvmPgzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "loc, scale = 0, 1\n",
    "laplace = np.random.laplace(loc, scale, 1000)\n",
    "x = np.arange(-8., 8., .01)\n",
    "laplace_pdf = np.exp(-abs(x-loc)/scale)/(2.*scale)\n",
    "plt.plot(x, laplace_pdf, label='Laplacian')\n",
    "gaussian_pdf = (1/(scale * np.sqrt(2 * np.pi)) *\n",
    "     np.exp(-(x - loc)**2 / (2 * scale**2)))\n",
    "plt.plot(x,gaussian_pdf, label='Gaussian')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4.2 Sparse optimization with Cross validation\n",
    "\n",
    "- K-fold CV\n",
    "- Leave-one-out CV\n",
    "- Random sampling CV\n",
    "\n",
    "Most straightforward CV is train vs. test\n",
    "- Fitting with different $\\lambda$\n",
    "\n",
    "<img src= 'Images/l2_6.4.2_Sparse optimization.png'>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5 After model selection\n",
    "\n",
    "## 6.5.1 Co-sparsity\n",
    "\n",
    "$$\\lambda|\\phi x|+ \\frac{1}{2} ||Ax=b||^2$$\n",
    "Where $\\phi$ is some transformation\n",
    "\n",
    "- Sometimes signal is sparse under a transform\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}