M368K Workbook - least square approximations- two homework problems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least square approximations\n",
    "\n",
    "by [Richard Tsai](https://www.researchgate.net/profile/Richard_Tsai3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QR factorization\n",
    "\n",
    "QR factorization finds orthonormal basis of the subspace spanned **successively** by columns of $A$.\n",
    "\n",
    "$A=QR$ where the columns of $Q$ are orthonormal and $R$ is an upper triangular matrix.\n",
    "\n",
    "Let $a_j$ be the jth column of $A$ and $q_j$ be the jth column of $Q$.\n",
    "\n",
    "span$\\{a_1,\\cdots,a_j\\}=$span$\\{q_1,\\cdots,q_j\\}, j=1,2,\\cdots,n$\n",
    "\n",
    "### Use QR factorization to solve the least square problem (=to solve the normal equations):\n",
    "\n",
    "$A=QR$ \n",
    "\n",
    "1. $y:= Q^\\top b$\n",
    "2. Solve $Rx=y$ by back-substitutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pivoted QR factorization and truncation and applications\n",
    "\n",
    "QR-factorization can be used to solve both **underdetermined**, \n",
    "**square** and **overdetermined**\n",
    "linear systems of equations. \n",
    "\n",
    "It can also be used to compute a numerically \n",
    "stable solution of an **ill-conditioned** system.\n",
    "\n",
    "\n",
    "With pivoted QR-factorization a matrix can be factorized into\n",
    "$$AP=QR$$\n",
    "where $P$ is a permutation matrix, making the diagonal elements $r_{i,i}$ in \n",
    "the uppertriangular matrix $R$ be ordered in decending order, i.e.\n",
    "$\\vert r_{i,i} \\vert \\le \\vert r_{i-1,i-1} \\vert $. \n",
    "\n",
    "$Q$ is orthonormal, i.e. $Q^\\top Q = I$. Assume the matrix $A$ is $n\\times n$. \n",
    "\n",
    "In Julia, $Q,R,$ and $P$ can be computed as follows:\n",
    "\n",
    "```\n",
    "F=qrfact(A, Val{true}) #Val{true} tells Julia to compute teh pivoted QR factorizaton\n",
    "Q=F[:Q]\n",
    "R=F[:R]\n",
    "P=F[:P]\n",
    "```\n",
    "\n",
    "The solution of a linear system of equations $A{\\bf x}={\\bf b}$ is obtained from \n",
    "$$R{\\bf y}=Q^T{\\bf b}, \\quad and \\quad {\\bf x}=P{\\bf y}$$\n",
    "If the matrix $A$ is ill-conditioned $R$ can be written on block form\n",
    "$$R=\\pmatrix{R_{11}&R_{12}\\cr\n",
    "               0   &E\\cr}$$\n",
    "where \n",
    "* $R_{11}$ is an uppertriangular matrices of dimension $r\\times r$, \n",
    "* $E$ is an uppertriangular matrices of dimension $n-r\\times n-r$, and **the elements of \n",
    "$E$ are very small**. \n",
    "    - The smaller the elements of $E$ are the more \n",
    "ill-conditioned is the matrix $A$. \n",
    "    - If the elements of $E$ are exactly zero the matrix $A$ is rank-deficient with $rank(A)=r$.\n",
    "        \n",
    "By setting small but nonzero elements in $E$ to zero we regard the matrix $A$ \n",
    "as rank-deficient even if it isn't that mathematically. Then we cannot \n",
    "find a unique solution. But we can, instead, **find a solution with the smallest 2-norm** to\n",
    "$\\hat A \\hat{\\bf x} = {\\bf b}$, where\n",
    "$$\\hat A = Q \\pmatrix{R_{11}&R_{12}\\cr\n",
    "                         0  & 0 \\cr} P^T$$\n",
    "according to the following algorithm:\n",
    "\n",
    "1. Compute the vectors $\\hat {\\bf d}$ and ${\\bf d}_E$ from\n",
    "$$Q^T {\\bf b} =\\pmatrix{\\hat {\\bf d}\\cr\n",
    "                         {\\bf d}_E \\cr}$$\n",
    "2. Solve the system \n",
    "$$R_{11}\\hat {\\bf y}=\\hat {\\bf d}$$\n",
    "3. Finally compute\n",
    "$$\\hat {\\bf x}=P\\pmatrix{\\hat {\\bf y}\\cr\n",
    "                            0\\cr}$$\n",
    "\n",
    "\n",
    "The goal of this exercise is to illustrate the following relations:\n",
    "\n",
    "1. the sensitivity of the solution $\\hat {\\bf x}$  from the experimentally computed condition number of $R_{11}$ (see below) as a\n",
    "function of $r$\n",
    "2. the condition number of $R_{11}$ computed with Julia's function 'cond( )' as a function of $r$\n",
    "3. the norm of the matrix $E$ as a function of $r$\n",
    "4. the norm of the residual ${\\bf res}=A \\hat {\\bf x} - {\\bf b}$ as a \n",
    "function of $r$\n",
    "\n",
    "What is important in practice is to find the rank $r$ for\n",
    "which both the residual ${\\bf res}$ and the sensitivity of disturbancies in\n",
    "**b** are \"small\". Often this is a compromize; the norm of the residual\n",
    "increases while the sensitivity (condition number) decreases as the rank $r$\n",
    "decreases.  \n",
    "\n",
    "\n",
    "In the following box, a function named 'gen_illposed_system()' is defined. \n",
    "It creates an ill-conditioned system $A{\\bf x}={\\bf b}$.\n",
    "\n",
    "1. Compute a system of size $n=18$, i.e. $A$ is $18\\times 18$ and ${\\bf b}$ is\n",
    "$18\\times 1$. \n",
    "2. Estimate the condition number $\\kappa (A)$ experimentally\n",
    "from the relation\n",
    "$${\\Vert \\Delta {\\bf x} \\Vert_2 \\over \\Vert {\\bf x} \\Vert_2} = \n",
    "\\kappa (A) {\\Vert \\Delta b \\Vert_2 \\over \\Vert {\\bf  b} \\Vert_2}$$\n",
    "by disturbing ${\\bf b}$ with some small perturbances $\\Delta {\\bf b}$ and see\n",
    "how ${\\bf x}$ is changed, $\\Delta {\\bf x}$. Compare your result with \n",
    "Julia's function `cond()`.\n",
    "3. Write a program showing a table of the following quantities \n",
    "as functions of the rank $r=8,7,6,5$: $\\tt norm(E)$, $\\tt cond(R_{11})$,\n",
    "the experimentally computed condition number of $R_{11}$, $\\tt norm(x)$,\n",
    "and $\\tt norm(res)$. Express your results in a way that describes the \n",
    "behaviour of the relations mentioned above as functions of the rank $r$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gen_illposed_system (generic function with 1 method)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function gen_illposed_system(n)\n",
    "# ILLPOSED Creates an ill-posed system Ax=b.\n",
    "#    ILLPOSED(N) creates an ill-posed system of size N.\n",
    "#    N >= 2\n",
    "    \n",
    "    r = max(convert(Int,round(0.25*n)),2)\n",
    "    U,S,V = svd(rand(n,n))\n",
    "    S[n-r+1:n]=linspace(1.5e-9,2.45e-12,r)\n",
    "    A = U*diagm(S)*V'\n",
    "    b = A*rand(n,1)\n",
    "    \n",
    "    return A,b\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Solving least square problems with QR factorization\n",
    "\n",
    "The least squares method is used to fit functions to measured data. In the\n",
    "following box, measurements of the $CO_2$-concentration in the athmosphere is\n",
    "tabulated. \n",
    "The concentration is measured in parts per million and is\n",
    "available for each month from 1980 to 2003. \n",
    "The data are tabulated following teh following format:\n",
    "\n",
    "1980 338.01 338.36 340.08  ........................ 338.21\n",
    "\n",
    "1981 339.23 340.47 341.38  ........................ 339.61\n",
    "\n",
    "Hence each row starts with the year, then follows the 12 measurements for\n",
    "each month. \n",
    "\n",
    "Pick out the monthly values (columns 2 to 13) and store them columnwise in Q\n",
    "\n",
    "```\n",
    "Q=maunaloa[:,2:13]'```\n",
    "\n",
    "Put the measurements  in time increasing order in one column\n",
    "\n",
    "```\n",
    "b=Q[:] \n",
    "using PyPlot\n",
    "plot(b,\"x\")```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "maunaloa=\n",
    "[1980 338.01 338.36 340.08 340.77 341.46 341.17 339.56 337.60 335.88 336.01 337.10 338.21;  \n",
    "1981 339.23 340.47 341.38 342.51 342.91 342.25 340.49 338.43 336.69 336.85 338.36 339.61;\n",
    "1982 340.75 341.61 342.70 343.56 344.13 343.35 342.06 339.82 337.97 337.86 339.26 340.49; \n",
    "1983 341.37 342.52 343.10 344.94 345.75 345.32 343.99 342.39 339.86 339.99 341.16 342.99;\n",
    "1984 343.70 344.51 345.28 347.08 347.43 346.79 345.40 343.28 341.07 341.35 342.98 344.22;\n",
    "1985 344.97 346.00 347.43 348.35 348.93 348.25 346.56 344.69 343.09 342.80 344.24 345.56;\n",
    "1986 346.29 346.96 347.86 349.55 350.21 349.54 347.94 345.91 344.86 344.17 345.66 346.90;\n",
    "1987 348.02 348.47 349.42 350.99 351.84 351.25 349.52 348.10 346.44 346.36 347.81 348.96; \n",
    "1988 350.43 351.72 352.22 353.59 354.22 353.79 352.39 350.44 348.72 348.88 350.07 351.34;\n",
    "1989 352.76 353.07 353.68 355.42 355.67 355.13 353.90 351.67 349.80 349.99 351.30 352.53;\n",
    "1990 353.66 354.70 355.39 356.20 357.16 356.22 354.82 352.91 350.96 351.18 352.83 354.21;\n",
    "1991 354.72 355.75 357.16 358.60 359.34 358.24 356.17 354.03 352.16 352.21 353.75 354.99;\n",
    "1992 355.98 356.72 357.81 359.15 359.66 359.25 357.03 355.00 353.01 353.31 354.16 355.40;\n",
    "1993 356.70 357.16 358.38 359.46 360.28 359.60 357.57 355.52 353.70 353.98 355.33 356.80;\n",
    "1994 358.36 358.91 359.97 361.26 361.68 360.95 359.55 357.49 355.84 355.99 357.58 359.04; \n",
    "1995 359.96 361.00 361.64 363.45 363.79 363.26 361.90 359.46 358.06 357.75 359.56 360.70;\n",
    "1996 362.05 363.25 364.03 364.72 365.41 364.97 363.65 361.49 359.46 359.60 360.76 362.33;\n",
    "1997 363.18 364.00 364.57 366.35 366.79 365.62 364.47 362.51 360.19 360.77 362.43 364.28;\n",
    "1998 365.32 366.15 367.31 368.61 369.29 368.87 367.64 365.77 363.90 364.23 365.46 366.97;\n",
    "1999 368.15 368.87 369.59 371.14 371.00 370.35 369.27 366.94 364.63 365.12 366.67 368.01;\n",
    "2000 369.14 369.46 370.52 371.66 371.82 371.70 370.12 368.12 366.62 366.73 368.29 369.53;\n",
    "2001 370.28 371.50 372.12 372.87 374.02 373.30 371.62 369.55 367.96 368.09 369.68 371.24;\n",
    "2002 372.43 373.09 373.52 374.86 375.55 375.41 374.02 371.49 370.71 370.25 372.08 373.78;\n",
    "2003 374.68 375.63 376.11 377.65 378.35 378.13 376.62 374.50 372.99 373.00 374.35 375.70];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGgCAYAAABfSOayAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3X901NWdP/7nZGxYFRNMUogsqZiQYFGwtrBkstuuhxAFWZduPaeHhiLbdZPqkZMWun4Uba3aQ3FtD1moXVqpXf3DlNZa7NYvSAMu2DYZDFYq1dbARIquBDYJIYgCMvP+/hHumzt37vvnvCczkzwf53BKZub9nve8N+w8vfd1XzdkGIYBIiIiojxVkO0LICIiIkoHwwwRERHlNYYZIiIiymsMM0RERJTXGGaIiIgorzHMEBERUV5jmCEiIqK8xjBDREREeY1hhoiIiPIawwwRERHlNYYZIiIiymsXZfsC/EgkEnj33Xdx2WWXIRQKZftyiIiIyAXDMHDy5ElMnjwZBQUBjqcYHvznf/6nMXPmTOOyyy4zLrvsMqO2ttbYunWr+fyRI0eML37xi8akSZOMSy65xLj++uuNn//850nn6O/vNxobG43LLrvMKC4uNv7lX/7FOHnypJfLMN5++20DAP/wD//wD//wD//k4Z+3337b0/e+E08jM1OmTMEjjzyCadOmAQCeeuopLF68GK+++iquueYa3HbbbRgcHMR///d/o6ysDG1tbfj85z+PvXv34vrrrwcALF26FEeOHEF7ezs+/PBDfOlLX0JzczPa2tpcX8dll10GAHj77bdRVFTk5SMQERFRlgwNDaGiosL8Hg9KyDAMI50TlJSU4Dvf+Q5uv/12jB8/Hhs3bsSyZcvM50tLS/Hv//7v+Nd//Vf86U9/wowZM9DV1YXZs2cDAF544QXcfPPNeOeddzB58mTte5w5cwZnzpwxfxY348SJEwwzREREeWJoaAjFxcWBf3/7nrCKx+PYvHkzTp06hUgkAgCoq6vDT3/6UwwMDCCRSGDz5s04ffo0brjhBgBAZ2cnJkyYYAYZAJg/fz4KCgqwZ88ey/dau3YtiouLzT8VFRV+L5uIiIhGGc9hZv/+/Rg/fjzGjRuHO+64A1u2bMGMGTMAAD/72c/w4YcforS0FOPGjcOXv/xlbNmyxZyW6u3txcSJE5POd9FFF6GkpAS9vb2W77l69WqcOHHC/PP22297vWwiIiIapTyvZpo+fTr27duHwcFBPPvss1i+fDl2796NGTNm4Bvf+AYGBwexY8cOlJWV4bnnnsPnP/95/OY3v8HMmTMtz2kYhu2qpHHjxmHcuHFeL5WIiIjGAM9hprCw0BxpmT17Nrq6urB+/Xr8v//3//DYY4/hj3/8I6655hoAwHXXXYff/OY3+P73v48f/OAHKC8vx7Fjx5LOd+7cORw/fhyTJk0K4OMQERHRWJP2Iu9EIoEzZ87g/fffHz6hsm48HA4jkUgAACKRCAYHB/HKK6+Yz7/44otIJBKYO3duupdCREREY5CnkZn77rsPCxcuREVFBU6ePIm2tjbs2rUL27dvx9VXX41p06bhy1/+Mr773e+itLQUzz33HNrb2/H8888DAD7+8Y9jwYIFaGpqwg9+8AN8+OGHWLFiBZYsWWK5komIiIjIjqcwc/ToUSxbtgxHjhxBcXExZs2ahe3bt6OhoQEAsHXrVtx777245ZZb8N5772HatGl46qmncPPNN5vnePrpp7FixQrU19ejoKAAt956KzZs2BDspyIiIqIxI+0+M9mQqXXqRERElDk512eGiIiIKBcwzBAREY1Sre3d2LDzgPa5DTsPoLW9e4SvKDMYZoiIiEapcEEI6zSBZsPOA1jX3o1wgXWPt3zCMENERDRKtdRXo7ayJCnQiCCzqqEGAEbF6AzDDBER0ShWV1UGAFjX3o2a+7clBZnRMjrDMENERDSKtdRXm+HlbDyBwvDwV78INS311dm8vEAwzBAREY0hZ+OJURVkAIYZIiKiUafVokZGjMqMNp43miQiIqLcJlYxRXv60RHrT5pmEtadL/wdDaMzozOiERERjWGiTqYj1o+6qlIAF8LLqoYaM9x0xPqydo1B4sgMERHRKCRGXNa1d6Mj1g8AKXUyYsl2vo/OMMwQERGNUi311fiPHd1IGEBhuCAptIi/xxN5t0VjCoYZIiKiUWrDzgNmkDkbT6SMwuT7iIzAMENERDQKyauYWuqrzZ+B0RNiBIYZIiKiUUYNMkByDY3882jA1UxERESjTDxhpBT7tkqrmdQ6mXzfQZthhoiIaJRZqenuK3rPiOeF0bCDNqeZiIiIxgDdNJNuOiofMcwQERGNEXKgeezFgzgbT+R9kAEYZoiIiEaN1vPTRbpwsmHnAcQTBlY21JhBRu09k69YM0NERDRKiLoYscmkINfFbNh5wAwyovdMvuPIDBER0SjhVBcjnhttvWcYZoiIiEYRq7oY8ZhV75mOWB/qqsocp6hyEaeZiIiIRoFWaXqppb7anEYqDA9/1f/uYJ+22FfssA3AcYoqV3FkhoiIKIe4LeJVyX1kACTVxTgtv26pr9ZOO+XL0m2GGSIiohwiQkm0px9tTbXm4yJY1FWVorW9OyXQqPUy8miLW/m6dJvTTERERDmkpb4adVWl6Ij1o3FTFEBykOmI9bua8nnsxYPmqMqqhhrtFJLV+8tTVLkeZACOzBAREeWctqZaNG6KoiPWj6rVWxE3DDPI2I2UiD2ZrPrIqHsy6eiWbud6oOHIDBERUQ5qa6pFOBRC3BgOIE5BBriw55Kuj0xLfbXjaiS5RqZ7zUJPIzrZxJEZIiKiHKAW/m7YecAMMgAQgnMvGLVg10sfGV2xr65vTS5imCEiIsoB6mokuUYGAAwAjZuiSUXBsnTDiJii0i3dFs/nqpBhGLl7dRaGhoZQXFyMEydOoKioKNuXQ0REFAh5JEUOMqsaahDt6UdHrB91VaVJgUaM6MQTRsqSbrGUWzyf7aZ3mfr+Zs0MERFRjmipr0ZtZQkAJAWZlvpqtDXVmgFHrmERIzq6ICM/vrKhJqmxnmrDzgNo9bCMO5dwmomIiCiHbG6OoHL1/4eEgZTVSG1NteZoi+C0H5N8vDyVJR5vbe9G16EBs8BYluvbGAgMM0RERDlkw84DZpDRLY3W1b24bXanCz4iyNRVlWpHdtSAk4s4zURERJQj0lka7bbZndiLaV17N2ru32YGGXn6Kl+2MRAYZoiIiNIURC2KLkDEE8PN8qw2gJTPq2t2Z0UNPm1NtUkBJ5+CDMAwQ0RElDZRi5LOjtO6pdHhgpA5ciLXyajn9TqiY9XlN9+2MRBYM0NERJQmL0W4VnRFtvJ5aytLASQHl3jCMLc9UPvLRHv6tf1lrBrrRXv6824bA4FhhoiIKACZ2nHa7rwbdh6wLN61G9HRBR+5h42XzsG5gGGGiIgoIC311ZabPNpRtzKQieXRuikgOeiIkRS7ESHdVJYcfOZMLUk5r/xzrmKYISIiCoDo16KbqnHq16Lr/wJcGEmpqyq1nALyMiKke/983sZAYJghIiIKgNyvRZ6qEVM4dv1a7GpuxLJpu80j/Y4IAfqAo15XrmOYISIiSpM8VSP6tai1KLpgIE8v6UZYplx+sba4V7xO/Gy1OmmsYJghIiJKkzxVI0ZORCCRa1FU6vSSPMICAFMuvxifn11hOwVktTpJfp2qtb0be97qR11VWcprhoNZH+ZeVZrz2xgIDDNEREQ+iZEV+UtfDiThUChph2uVOsoCwAwyALRhQz7WanWSfE7d8eGCEKI9A4j2DCS9JnnX7jKHT587GGaIiIh80hXuiikfAIgbhuOUjy7QiPoapxEWv8W78nvK7yv+nk/dfwGGGSIiIt90QUQOBPLPbsOBWsBrd3w6xbtWgSbfggzAMENERJQWq5EVu0Ci9pXpiPUBuLBT9pLHO7G5OZLx5dFqjU6+bWMgMMwQERGlqaW+Gv+xoxsJIzUQ6AKJPD0FANGegaSRnGjPgDk9lclwIU+JAcjblVAMM0RERGnasPOAGWR0gcCqpkU3JeV3esrPNeveP1+6/soYZoiIiHwQU0UAtEujO2J92NwcsTy+pb4aHbE+RHsGLDv3Wk0vudn+wK6eRg0y6pRYvgUahhkiIiIf5KkiXdGsPFVkZXNzBDX3b9N27rU7zmn7A7tuw8BwSKqtLElZ+i3X8OTDNgYCwwwREY1J6Y5uyCMr8nFymHAKBH4799ptf+BmNZLT58qXERkhZBhG/kSv84aGhlBcXIwTJ06gqKgo25dDRER5RJ0ekkcoRPdbuSDX7ZSNCCRulzZbde71sjTa73tnS6a+vwsCOxMREVEeUKeHoj0DWNfejcZNUXMlkVwQK4KPlZb6ajNMFIYLzC0GdDbsPIDW9m7Lzr2rGmqw7vzzbqjvnctBJpMYZoiIaEyRQwNwYSVPR6w/6WerUZJWJWyoU0XP/v4dbSARASZcELLt3LuqoSZpekp9P1njpmjKNNVYxJoZIiIac+Sak8Jw8n/XW60sEtQeMepy6neOf4C6qlLf9Szq81bFvo2bouaO3G1Nta42mBytGGaIiGhM0nW/PRtPOE7ZuO0RIwKNGo68Fh7rin3VIGP1urGCYYaIiPJCuquPdMfI3W9nT73cnGpyWlnktkfM3kPHU8KRn2XVclAR7ycHGfV1+bSsOgismSEiorwgQkDjpmjS4yIEdB0aQKs0/WNHbRpXV1VqjnTIoyt2NSibmyOWxbct9dUIF4S09Sy6Ql8301Bqsa8aZOTXeQl1owHDDBER5YWW+mozdIhAI0KAeNxp5ZF8TG1lSVLxrzgHMBxwaitLbAONrkeM+h6rGmrQvWZhSniRA03N/dtc1dPYvd9Yx2kmIiLKG21NtWa9SNXqrYgbhhlC3PZYUVcStSq9XsR0lfyzyqpHjKBbdi0eFz/LNTtOy6rl0KYr9vUzzTaaMMwQEVFeaWuqNYMMAE9BRld3IwKAXfGtzKpHDICkER+rzSVFOHLb/VcdfZJfN9wXp9+8B2MVwwwREeWVDTsPmEEGAMIhfVGwTrp7GgGpIzuCHFacVkLZjeyox8rvp75ODjJjafWSimGGiIjyhvgyn3L5xXjn+AcIh0KIGwYaN0XNgli7KZd09zQKYkWV08iO/DMAy2XaTv1wxhKGGSIiygvqdAsAfGX+hdGJxk1R1FaWOo6wpBMIRmpkx+n4cCikrbUZq7UzDDNERJQX4onkYl8guVFdR8z9lItafCv2U3IacUl3ZAdw3rHaSdehAcSNC4FGXLeXQDXaeFqavXHjRsyaNQtFRUUoKipCJBLBtm3bAACHDh1CKBTS/nnmmWfMcxw+fBiLFi3CJZdcgokTJ+Luu+/GuXPngv1UREQ06qxsqMGcqReKa+XlzXsPHQcA1FaWuN6xWi6+7To04LifkuBnWXVQhnf1Hl5GLlZyyZtkjtUpJ09hZsqUKXjkkUewd+9e7N27F/PmzcPixYvx+uuvo6KiAkeOHEn689BDD+HSSy/FwoULAQDxeByLFi3C2bNn0dHRgaeeegpPPvkkHnjggYx8OCIiGl1WKl/WaiO5zc0Rx3PoesCIgOC2kV02dquWr6etqda87nAoZF7/WAwyABAyDCOtnsclJSX4zne+g9tvvz3lueuvvx6f/OQn8cQTTwAAtm3bhn/4h3/Au+++i0mTJgEAfvCDH+Cee+7B//3f/6GwsNDVew4NDaG4uBgnTpxAUVFROpdPRER5THzBi2DhpvGcLqCo9ThO5/P6vkD6xcO642vu34az8QTCoRBWzJuW87Uymfr+9l0zE4/H8cwzz+DUqVOIRFKT8CuvvIJ9+/bh+9//vvlYZ2cnZs6caQYZALjppptw55134vXXX8f111+vfa8zZ87gzJkz5s9DQ0N+L5uIiPKMVQiQO/lubo642jXaTfGtbj8l+RrUQCSmeNT3Va9b3W1bhBe3tS5qUFGnytx0Px6tPIeZ/fv3IxKJ4PTp0xg/fjy2bNmCGTNmpLzuiSeewMc//nHU1dWZj/X29iYFGQDmz729vZbvuXbtWjz00ENeL5WIiEYB3QoiObjUVZUlPWcXLHQjF/KoiFUjO3ENal8XuYZFfV/1unW7bXstHpav2W2fmrHAc5iZPn069u3bh8HBQTz77LNYvnw5du/enRRoPvjgA7S1teEb3/iG6/OGQtaJcvXq1Vi1apX589DQECoqKrxeOhERZUG60yu6kNIR6wOAlDoR8feOWF/Sed0sqXYTEMRUlG6ERt36wC5cAfDdJ8Zrn5qxwHOYKSwsxLRp0wAAs2fPRldXF9avX48f/vCH5mt+/vOf4/3338dtt92WdGx5eTlefvnlpMeOHj0KACkjNrJx48Zh3LhxXi+ViIhyQBC9WXS9YXTt/YVoz4A5YiPClFiBJM4nT1OJc7sJCGIVkxpEdAHCqqeN2z2ZdNLtUzMapb1rdiKRSKpnAYanmP7xH/8RH/3oR5Mej0Qi2L9/P44dO2Y+1t7ejqKiIu1UFRER5T95KbOblUJ255FXEIkVPU7nlcOUuqQaGJ6msgsIqxpqzIDgZxWTegyAtHa/Vld0qe+V60XAmeBpZOa+++7DwoULUVFRgZMnT6KtrQ27du3C9u3bzdccPHgQL730ErZu3Zpy/I033ogZM2Zg2bJlePTRR9Hb24uvf/3ruOuuuzjyQkQ0iompIauuu7rpJnV6Sq1nWfJ4p7kU266brzw6sqqhxjwegOsGe4LbzSFl6jGsdQmepzBz9OhRLFu2DEeOHEFxcTFmzZqF7du3o6GhwXzNj3/8Y/z1X/81brzxxpTjw+Ewnn/+edx5552IRCK49NJLsXz5cjz88MPpfxIiIsppdVVliPYMpIxqWE03qat/5NcMF+MOmGHCadpGnTLyw0/RrfrZ1Pcf67UuQfEUZkS/GDvf/va38e1vf9vy+SuvvFI7akNERGOHPL1i15hOPA9Au4WBWGHkdbSkMFyAFfOmuQ4Rfopu1WNaNaFGXuU0FmtdgsK9mYiIKGPEVBGQGkLkkGJXA9IR60O0Z0A7jaRbKm1XbAwgZZrJTaDxU3SrHqPWsuhWPpE/DDNERJQxavGt+NJev+MA4i4b0G9ujpirh9TpKXVjSavRErGUWw09qxpqkgp8rfjZHDLdDSXJvbRXMxEREVlpqa82lz4LjZui5q7PwIWgYUVXdAvYj5bIAWXDzgOI9gykhB55VMYqeLRqNp+Ur6s1jRocCg5HZoiIKKPkrQbEiExdVSnammrNx9U6F930lDyi0hHrs91UUj6XlymioLcgoJHBMENERBnXUl9tBplwKIS2plrzcSB1WshqekqQVzJZcbuFgSxTWxCk2wWZ7DHMEBFRxm3YeSHIxA0jKYjoRkjkwl/5HPJoiFOdi5/Ow0FvQSBCjNX+UuIzcoQnPQwzRESUUW76s+iCgTw95WcfI10wcTOqEuQWBCLEiEJjtW8O4K5xH9kLGYbLcvIcMjQ0hOLiYpw4cQJFRUXZvhwiIrKg67Wiq4URr9VNt8grmbrXLPR9DaKA2G14kN9X9KTxeg75/XWN88ZakMnU9zdXMxERUcaoxbdqLYy84midFHQEq5VMXvjZT8lqC4LuNQtT9oNy8/7imMdePGg+7meTSdJjmCEiooxRN0XULYm2mvqRH/cTIuTzeAlE8vuumDct5XndxplO5EAFIK1wRqlYM0NERIGzW70DALWVJba1MH62D9Dxup9SprYgEIFKECGJezIFg2GGiIgC52Yl0e//Mmg59eNn+wCVn0CUiS0I5ACl7i/ldjsFsscwQ0REgXNaSQTAdnPIILYC8BOIgt6CQHzm2soS1FWVJZ3Dy3YKZI+rmYiIKGN0K4kAfVdf8aXutblcLjeky+Vry4ZMfX9zZIaIiFIE9SXcUl+d1J8FSF2SLY/i1FWVoiPWn/S4eE+rRnd+muO5EcQ94GaTI4NhhoiIUgQVENSVRB2xPsepn9rKUk+N7nRTWkse70zZXFK+Jrsgwq69+YdhhoiIUvjtniuzWklUV1Vm+56Cl86/uq69TtdkhV178w/DDBERJVGnV+SAUFdVmrGVRDJ1espNcFCPEV17xXNuw5h8nbpAwyCTexhmiIgoiTq9IgJCOBRCR6wftZWljudId2m1rtGdU4BQjwGQ1HnXyzYEcqARtT4Au/bmKnYAJiKiJHKH28ZNUTPIxA0DdVWljl/m6RbO+un8a3UMAM9bGQjs2ps/ODJDREQpWuqrEe3pR0esPynIdMT6HUdJ0ike9jM95eYYLyM88nnZtTc/MMwQEVGK4VU7F4JMYbgAbU21jtsByI/7qVfxMz1ldYxQW1mCzc0RV9cusGtvfmGYISKiFPHEhZEY3aiGU82LrnhYDhxWU1Fi48nW9u6UqSir4KCbskqnAJlde/MPwwwREaUIF4hi3xJzKbWu54pd/YvdiqRMNboT0ilAdnPsWOramw8YZoiIKIkaKOx6rjjVv1itSAq60Z0qnc677NqbfxhmiIhGkSBa8OtGJuRAI3fAdZqqURvmAdAGmnQb3dHYxjBDRDSKBDF9Y1Wrom4YqQYZEaTEa9XX1FaWaANNEI3uaGxjmCEiGkWC2IbA6rxOHXlFkKqtLEl6L/n966rKkupVgmx0B3CX6rGKYYaIaJQRK5F0gcDvF7qbjrzJu1+XmcdZBSmrqahVDTW+G91lurCYchPDDBHRKCNWIoVDoaRA4PcL3U39i+C0JNvqnOqxgL9Gd5kamaLcxjBDRDQKyNMravfes/EE/u7fX8Q7xz/w/IXup1+LmympTDS6k99fXJ+fqSrKPwwzRESjgDy9AgAdsX6z6R0AvHP8A8d9lXT1JiJ0iOfF9JRdvxY3U1JBN7pTr13dIJMN7kY3hhkiolFA/dIXAUSEGbHjtd2Uja7eRHTkHa6DSd4tW9c8z8uUlCqdRnfqtYtAJbZj6Do0YPvelN8YZoiIRomW+mqzB4zct0WMkNRVldoGArt6E6tNJuXwks7ICpBeszr5fcQUm7hmtxtkUv5imCEiGkU2N0dQc/+2pGXOciiprSy1Pd6u3kQdZZH3MAL0Iyti1MZuL6OgllOrtUIdsX5fI0SUfxhmiIhGETG9ovIyQmJVwGsVdMRjuiDjZhVRkMup50wtwZ6eAXOnb/XaWTszOjHMEBHlMXlUQ93tOdrTnxQS7L7Q1fPIBbyNm6KYM7UEK8+HEquVSn6XQwe5nDpcEDKDjNV+UDT6MMwQEeUxMaohplfkaRVRL6ILNG7P07gpahYRA9Yrlfwuh1anmOTjnWp8VOkUH1N+Y5ghIsqydGpGkrvulmpHNcQ57Mj1JvJ55AJaEWyswoKb/jIqdYpJXk7dEet3rPGR71M6xceU3xhmiIiyLN2aEflLWxT/ioJbqxU8upA0Z+pwIW9HrD/pPPIIjV1YAODYX8bu2qM9/UnLqZ364sjSWdZN+Y9hhogoS3RTLAAsw4Md3aiI1TSLVUgSwUYEGXl0Zc7UEtRWpoYL8bNYEu5nikddhSSCjJfl1Oks66b8xzBDRJQluhGZde3dWL/jgOeRCbuuu14Ka63O47Q0Wg4y6ueRf7a6djnIFIYL0NZUy5oXco1hhogoS+Qv/I5YH+qqyswv9HAohLamWgDWdTNiZEecQx0V6Yj1YXNzxHxeV5irWw0lTy25CRPpTvGIXb47Yv3aMMYpInLCMENElEXJNSPDLfdFoNmw84D5nK5uRh7Z0YWJaM+AGQqsCnO9robSSXeKR+zyzVVI5BfDDBFRlrXUV5tTSwDwlfmp+yzpvtTl7QsEtR5GFAFbTUEFtRrKL65CoiAwzBARZVnjpqg5tRQ3DKxr70ZhuMDVsZubI2Yg0E0juem9YrUaaiSazXEVEgWBYYaIKIvk6Ry56FWMpKyYN81VjxjdNJKXUQ8/PWKCwFVIFAR30Z+IiAInhw1R7CsTeyw5rSTSTSMB9qMe6saPVucgygccmSEiyhI1bKj1Lh2xPse6EbtpJLejHm63AdB1KpZXVKkrrrzseE2UDoYZIqIsUb/4dVNCdit70i2ebW3vRtehAW1nX3WTSkDfF0ddUaW7NqJMY5ghIsoBfgph0y2eFUui1eZ8ch2PfA63QcnvjtdEfoUMw8i7UvGhoSEUFxfjxIkTKCoqyvblEBGNqHQ2ptS9XjfFZBdExGtEfY0YfVEfY5AhVaa+vzkyQ0SUZ9LdmFImj7bolnZbHaNb+ZSN1VBEAMMMEVFO8TLq4mXPJTtel2VbrXzyumM2UVAYZoiIfAhyqkfmdtTFz4iKFbsOwbrXWq184nYElC3sM0NE5IMIHWo/FvFFLpYrW2nVHAsMf/mL/ZDE81ajLi311WYA8Tu1I5+7e81CrGqo0X4uu+vQfQa78xAFjSMzREQ+6Fb2eJnqsRuBkTd4tBt1cTOiYjeC1Lgpql2WrX4uQbd6Sjwm/q7eH25HQCOBYYaIyKd0pnrchCGxT5Ju1MVtozs3ocnt0m7dtBm3I6BcwDBDROSBOtIhF8+GQyFPIxF2Ychu1MXrnkvq425GkBhEKJ8wzBAReaCOdIjQIXa87jo04Ol8upVETqMuXpvlBVksTJSLGGaIiDyQg8HP9r6Nd45/gLqqUnPKpiPWn7TRozoNo47sqCMwf/fvL+Kd4x/Yjrr4mdrJ1q7YRCOBYYaIyCOxd1FHrB8Akopo1aXKKnlkB0DScut17d1mOPK7RYEVL8uvifINl2YTEfkwZ2oJwiHr5de1lSXasCAvW1aDzKqGGqxqqEka3VGPBWC53HnDzgNolYKS/Li8/Lq2ssR2+bXuHCqrpeVezkEUFI7MEBH5EC4YrpERIx1etwPoiPUh2jNgeYzVCIzXrQx0xb51VWWI9gyknMfLdghBbqlAlC6GGSIij6wKdL3Uo2xujlguvXYKQoD71Um6YmH5HB2xPl/bIaTbZ4coSAwzRDRmBLEFgdMXttt6lHRqWLysTrL6PPI5RKjyGkK4SopyBWtmiGjMSHcLAiB1pEMON6saamzrUdT3c7OFgJUgtjLwew65XkY9h3ieaCRxZIaIxowgpkbkkQ6r5nV2Gy16aXhnJ4jVSX6HIYnGAAAgAElEQVTPoa7Iks/BehnKBk8jMxs3bsSsWbNQVFSEoqIiRCIRbNu2Lek1nZ2dmDdvHi699FIUFRXhM5/5DD744APz+YGBASxduhRFRUWYMGECbr/9drz33nvBfBoiIgfyaqKa+7d5CjLqCh55lEZewSPeQy7iFcfqaljkx90svQ5iZCedc+hWZK2YN831exMFzdPIzJQpU/DII49g2rThX9qnnnoKixcvxquvvoprrrkGnZ2dWLBgAVavXo3vfe97uOiii/CHP/wBBQUXMtPSpUtx5MgRtLe348MPP8SXvvQlNDc3o62tLdhPRkRkwW8DOXUFjxil0a3gUc8pjtUFGS+BKoiRnaBGhwS5XsbvOYjS4SnM3HLLLUk/r1mzBhs3bkQ0GsU111yDlStXoqWlBffee6/5munTp5t//9Of/oQXXngBXV1dmD17NgDge9/7Hm6++WZ897vfxeTJk9P5LEQ0CgVRtKs7zs/0SjrTVEGt/vG6lUGmz2EVCrlbNo0k3zUz8XgczzzzDE6dOoVIJIJjx45hz549WLp0Kerq6hCLxXD11VdjzZo1+Lu/+zsAw1NQEyZMMIMMAMyfPx8FBQXYs2cP/umf/kn7XmfOnMGZM2fMn4eGhvxeNhHlmaD7mYjj6qpK0dZUm1Lf4hSQdCt4aitLbN9PnC+d1T8i1Omuy22oC+IcwsqGGstQyBEZGmmeVzPt378f48ePx7hx43DHHXdgy5YtmDFjBnp6egAADz74IJqamvDCCy/gk5/8JOrr63HgwPAcbG9vLyZOnJh0vosuugglJSXo7e21fM+1a9eiuLjY/FNRUeH1sokoT8n1GaKew+uIhqhXkYOM3GVXrEBq3BR1tapJXcFTV1XmepWUmxVEuu66ItQ1boomrRZyuxKrtb0bXYcGtNfp9nPrPls6dTtEQfE8MjN9+nTs27cPg4ODePbZZ7F8+XLs3r0biUQCAPDlL38ZX/rSlwAA119/PXbu3Ikf//jHWLt2reU5DcNAyKYt+OrVq7Fq1Srz56GhIQYaohwW9NRQuv1MRBCorSyx3EOpIBRK2mPJjjoiIc7hZgrJzRSXbjRK3Q9KnMttqAsXhMzNMOVzN26Kmo+7vZ9B19wQpctzmCksLDQLgGfPno2uri6sX7/erJOZMWNG0us//vGP4/DhwwCA8vJyHDt2LOn5c+fO4fjx45g0aZLle44bNw7jxo3zeqlElCVBTQ3JoSieMBAOhVL6mYjpDrdTQ3VVZSnPewlIVt1/RZ8Zu8Bldax6n6zqa+Sdub02uku+B8OBZv2OA4gbhjnl5lYQNTdEQUq7z0wikcCZM2cwdepUTJ48GW+++WbS893d3Vi4cCEAIBKJYHBwEK+88go+9alPAQBefPFFJBIJzJ07N91LIaIcEVSxqxyKug4NIG5cCDTiXG4DktXojpdVTU4jEqsaarRTSGKKRx35EaMtToFGDUdW2yA4Uf/vIu6nlyADWHcVVj8D0UjxVDNz33334Te/+Q0OHTqE/fv3Y/Xq1di1axeWLl2KUCiEu+++Gxs2bMDPf/5zHDx4EN/4xjfw5z//GbfffjuA4VGaBQsWoKmpCS+//DJ+97vfYcWKFViyZAlXMhGNMun0c9GdQ4xKxI0L/9UvgoCX0Qm1W6065aPjtkdMR6xPez55ikc91vxcmtEMXX2NbprKi5b66qTdvuOGwToXynueRmaOHj2KZcuW4ciRIyguLsasWbOwfft2NDQ0AAC++tWv4vTp01i5ciUGBgZw3XXXob29HVVVVeY5nn76aaxYsQL19fUoKCjArbfeig0bNgT7qYgoJ/jt52Jl76Hj5t/DIX1AsKMGATdTPoC7HjEAEO0ZsD2fKJB1O1KlXq+ob3FzzVYaN0XNERkxxcQ6F8p3nsLME0884fiae++9N6nPjKqkpIQN8ojGiCBa7uv6mQDDIyrhUAhzplovi1avRQ4ecjt+wL6A1WnaTDznpiDWbRGzGnZ0hbpei27lc8jL0hloKN9xbyYiCpQo2gWgHfnoiPVhc3PE9fl0/UwAmH93s5xYDQatmlAj90exmvIRr1fDSKvFCIt6PrcjVbpRGxHaxJJyNdA4Fd3KU1qiRkYtCmbhLuUrhhkiCpRctKv7go/2DHgaobEaURF7ATmNKLS2d2PPW8mFt/I2BLWVJUlf4k5dfHVhxG1BrNuRKl1tjnzNauhwcy/drEDy2kmZKFcwzBBRoFrqq9ER60O0Z8B8TA0kbkcArKZy5L+r/V1U4YIQoj0DKUuy3a6ukpeH62pY5kwtcRUC3C7LBjKzWogrkGg0Y5ghosBtbo6YX9ail4nVKiC7L1l5NEE3lSMfbxWQ0l0mLkaaRMM6tYbFDTaZI8qskGEYeTdJOjQ0hOLiYpw4cQJFRUXZvhwislC1equ5cia29mbzcT89Z9Il3lOMqnh5b7vCWTddg/12RM7EJptE2ZSp72+OzBBRRmzYeSBpCXDjpmhSEBjJIAOkt0xcLr5VO+/qalhUfqd4gt5kk2i04sgMEQXOalmxCDaZGslwc01+RmYEufNu95qFno71y2qXbzVMcYSG8kGmvr8975pNRGRHN/LS1lRrBplwSB9SZGJEws0u1F6vye8Oz+l23vWrpb7anM6qWr1VW0Ts9X4QjTYMM0QUKMuW/9KUk1MQkLcxEK/1Oz1lVXxbW1liGWg27DyAVqmpnt8w1GrzGvU97OjCYLam64hyEWtmiChQ6nSHlyXJMq8dc61Y9VepqypDtGcAHbE+23oULyuR1OkxueZFXItoAuil5kUNg6KwmkGGaBjDDBFlTLpLkoPY28mqlsTtfklums0JasGu+lnVXb7dfB719fIKMQYZomEMM0QEIP2iW93xIgiI58XxXlrw+93bycvnsRv98bISySmoeR1h0hb7alaIEY11rJkhIgDpF93qjpfDgnp8S321bVCwqlNp3BS1fL1cg+L287TUV5thKYidveV6n5r7t5mfwc97yKNC8v2Irb3ZLAoeqUJkolzGMENEANIvuh2Jol3xBa4GGl3gcns9mVilpAYkAL7eY6UmyMgrxPysyiIajTjNRESmdItuM12029ZUa/assatzcXs9dsXJ8YThe9pNDUh+CqDd3A+303VEox3DDBElSbfoNp3jRZ2LLiSIACE3jrMKTHK9jHo98vN2xcliFEh+XFyH3Uokq12+5fsjP+7m/nCTSCJ7DDNElCSdott0j3fbvt8pMKlLonWjJG5GO2orSz1tUKk+36oJNfIqJ46oEAWDYYaITH57wgR1vG7Uwm2di3x+3ZJo+WfA22iH22kzNSCp7yGHF46oEAWHYYaIAKTfEybd44V06lyszi+fx8/1uJ0243QQUXZwNRMRobW9Gx2xPnP6RV0dU1tZYj5u1YLfatomnjBQV1WqnVKxOp/VcmmrwKRb1SOuRz2PeL3bKZ5s7clERO5xZIaIEC4IIdozgLqqspR6E7WY1arw1WpUIlwQQkesH7WVpUmPq3UwctGuGiAaN0UxZ2oJALhe1SO2DdBNR7kdJUln2iwTO38TkR7DDBElTe2saqgxRzoApAQZP6ubxPHiZ90IiwhR0Z5+dMT6zefEUmwAtt1u1esKuv7H6rNYcVvMTETpY5ghIgDJX9RiCTPgrgW/m1EIEZCsztdSX20GmbqqUjOAiJ/l3jJOgqjf8dvbRb0X8vuJYMYNIomCxTBDlOeCnM7Q9WRxU/jqZhTCTSGtmErqiPWj5v5tSaFHfBY3gmgy57eYV3cv1rV3Y/2O4X2VRFAjouAwzBDluSCnM9QaEwCu+sW4mUpy039GBAgRZOTQ4yUAeAkiQde26O6FCDLhUIgbQxJlAMMMUZ4KYjpDLbrVda5dMW9ayvl17JZUe6lfGamdsoVM1LboRmTETtdemxASkTOGGaI8FcR0hq7oVpwHAOqqSrVFwXaBRp1K8lK/km7Rrp9gol5LPGGg69CANgx6GalRR2Ria2/2tS8TETljmCHKU/KXcEesb3hZ9fn/+penM+y+gOVziPCz5PFOAEipVRHH29Wb6EZV3NavBFG063bllN1x4h6qYdDrSE3jpqjliAwDDVGwQoZh5N3mIENDQyguLsaJEydQVFSU7cshyir5v/YBmF+eXpZUi3OIEOJntY3VqIrbqS7dLtUiDInnRaBymk7qiPUh2jPg+fOIWh35Hnr5LOq9qKsqTdoY0yogEo0Vmfr+5sgMUZ6TpzMA4CvzU/clcjOikc5O2emMquh2r7Y6p3qMel75mN//ZdDT51FHlcQUm5ul6UHdCyLyh2GGKM+p0xlqnxg3rKaHnIppAZijJlb1JU5bB/iZGnI6BoCnImKrUaVwKOQ54AWxLJyIvGGYIcpjclM5eTpDfAGvmDfNsbkbAO0X+ZTLL8Y7xz8AYF9M63VURcduJZTXY6w+j/o57K5VbuAnAo3bVUjcbJJo5DHMEOUpp8Ag+sTY7ZlkNxX1zvEPzKkWwH7ExGvBrY6fqS5dkz+vUzxWo0oiJM6ZWmI5rUVEuYFhhihPqV/C6ohJR6zP9gu4pb7aLJQV1HPEEwZqK+1rR/yMquj46S+jHiN2/vYyxaOGPbswxkBDlJsYZoiyIIius/LzVlMlTlMsm5sj5mvsgog6YqJevzxCEg6FPNeF+OkvY3VMXVWZ9vWseSEavRhmiLIg6K6z6XwBO03v6EZM1OsXrxFFyF2HBnRvpeVn9Y/umHjCSJkWk1/vdik0a16I8g/DDFEW+G3uZiWdL2C76R273jGiI7C807XdDtdWo1EihKiByy6I6cJbuCBkvr98jJeAGPQ+TUQ0MhhmiLIkqFoTmdcvY7vpHXFtViMmqxpqzOASDoWS2v/rpomsRqNECKmtLLW8RyqnbsbiXF4DYib2aSKizGOYIcoSETzEiIg8xeN3FMDLl7HT9E5tZYnj1NWcqSXY0zOAuGFod7iWR0h0o1FLHu9EtGdA+z7p7ljtJSAGsWknEWUPwwzRCJK/NOXgIQLNksc7UVdVZjsKYDf6AgC1lSW201fy9gFOje6s3iMpdJ0PMuoUlVOzOxE2dNLdsdrrEu8gNu0kouxhmCEaQfKXpmz21MvREetHtGfAcqRCdw6r0RcRiHSjE362D9Dxu8O1GjZWzJsWWO2QuC6vS7x1IzLyjtdi004iyk0MM0QjSP3SXNVQk9RpNu5i31e3xcNWoxNBFB+ns/+QGjbEfQiidshvwFKvXw4y6o7XRJR7GGaIRpjcrC6pN8v56Rq7LQjkcwDWtSFOoxPpFh/7XQputzpKVzvkhtW2DII67WZHHZGJrb3ZUyAiouxgmCHKgs3NEdTcvy0lyDhtQSDXy6jTNfLzbkYn/NSWCH6WgrsZzfEyNSSIz6sWLKvTbm6a3ambdsrXwUBDlLsYZoiyQG0yp24UCei/NNWaG3n0Zbj7bWnKyhurL2M/tSXpsBrNEWorS5I6EsvXqqNbgSS6//pZgWS1aad8Hez+S5SbGGaIRtiFtvvWTeasvsx1NTfyzwnD3fSP29oS3copeVpH17fGajm17rF0am+CXIGUznUQUfYxzBCNIPlLU2ziaBUm5FEAq+XYYpqotrLEXMHktDeRly9uXWBQd9vWfTa30t2GQb1mvyuQuB8TUX5jmCEaQW6/NNXn5QAhziH3aamrKnP9xevli9vLCiU/y6nT3QcpqBVI3I+JKL+FDMPFWtAcMzQ0hOLiYpw4cQJFRUXZvhyiEaGOfMi1M25DhN+9h8R7i/oa+Rrkx7L1xV+1eqt2BRK79hLllkx9f3NkhijL3AYMdZQESK2Zcfri9rv3kNXKJ7+roYLEFUhExDBDNEKsQosIGB2xPmxujpiPOwUMNUC47b6rvtbNKIZu5RMA36uhgtqdmiuQiAhgmCEaMVajIkK0Z8AMBFYBoyPWByC1H4uXL26vDfPsVj756bTb2t6NrkMD6Ij1p7xeXlLthCuQiEhgmCEaIU6jIuI5q4CxYeeBpH2b1ADh5YvbbcM8t7UnXkJEuCBkjqbodqd2u6SaK5CISGCYIXIpiKkRp1ERq4AR9CiE24Z5usAgHhN/Vz+bl60YRKCRe8O4XVLNFUhEJDDMELnkt3hWDUHyqEg4FEpqZGcVMIIchfCyGaPb5nfq9ThRgxh3pyaidDDMELnkt3hWDUHqVgZdhwYcA4aXAGE3gqRr85+tOhO5yR0A7k5NRL4xzBB54LV4Vj0m2tOftIWB+N8gA4bdCJJVTUo26kzUJdVqDQ0RkVsMM0QeuS2eVY8RQSYcCiWFlyWPdyLaM6A9BgB+d7Av6WeZrlZHF4TEe+iCl5el0EGRi33lJdUMNETkR0G2L4Aon7S2d6NxU1Tbc2XDzgNolRraqeZMLTFHIeQQtLk5Yu7VpGqpr8bfThvec0m8jyACgNj0UT1uVUMN1rV3o+b+bdqw5HQO3WdXr0E+j91nV18rBxn5esXjXIlERF5wZIbIA9EfRR1REKMudv1RwgUXgoxa4JuJRnfqCNKKedM8n0O9fj8F0Co3xcwjOUpERPmPYYbIQas0ciHXuYhRiimXX+zYH8XLCiIdP7U66uooAOZojdtzqPdAHC+uSXyO2soS19NCXFJNREFjmCFyIEYkaitLtGEEAOqqSjFnaon2+KB6xHip1bEKT6saasxw42U/JXEPVjXUpASi4c9f5uo8RESZwDBD5CC5yVvql3ZhuMC2P0pQPWLcNrpzE5687qckHy8HIsD9jt1ERJnCMEPkgm6aB3AXCoKYVvEyTWUVnoTayhJsbo6kNdVFRJRLGGaIXJKneQDnTRa9bH8QZKM7XXjKxKaMuqJiIqJsYJghcklM86isQoGX1T+ZbnQXxFSXHNp0RcXqtRMRjRSGGRrVgtgcUrxWFAGLuhndjtW6jRfdLIdOZ/m1G0FMdXXEhpv36YqKrfrkEBGNBIYZGtWC6I1iFyp0gUbmZUm17rVWy71Humvvhp0HUjoIq0XB7A1DRNniqQPwxo0bMWvWLBQVFaGoqAiRSATbtm0zn7/hhhsQCoWS/txxxx1J5zh8+DAWLVqESy65BBMnTsTdd9+Nc+fOBfNpiBRyJ1y5U6+XZnF2UzRuRiRa6qtdL4eWXyu2PfDS+VcVVNfedO8BEVEmeRqZmTJlCh555BFMmzYNAPDUU09h8eLFePXVV3HNNdcAAJqamvDwww+bx1xyySXm3+PxOBYtWoTy8nJ0dHTgyJEjuO222/CRj3wE3/72t4P4PEQp/DSck6U7ReN2SbXutepeRV6CWGt7t9mxWL1WuajYDTa6I6Jc5inM3HLLLUk/r1mzBhs3bkQ0GjXDzCWXXILy8nLt8b/+9a/xxhtvYMeOHZg0aRI+8YlP4Fvf+hbuuecePPjggygsLPT5MYjs+dkcMghellRbvVYEGq9BLFwQMouH5feUN3lkECGi0cD3RpPxeBybN2/GqVOnEIlEzMeffvpplJWV4dprr8Xq1avx/vvvm891dnZi5syZmDRpkvnYTTfdhKGhIbz++uuW73XmzBkMDQ0l/SGyo06vLHm8M2nEY8njneZzXqZbvLBaDq1Oezm9Vuy07TWIqZs3rmvvRtXqrSmbPBIR5TvPBcD79+9HJBLB6dOnMX78eGzZsgUzZswAADQ2NuLKK6/E5MmT8dprr+Gee+7Bm2++iV/84hcAgN7e3qQgA8D8ube31/I9165di4ceesjrpdIYJhf+AjB3jZ499XJ0xPoR7Rkww4RdIbDVaigxhTNnaknKFIwozgX03XF1K5/salLEJpZeu/bK7yXuRdwwEA6FGGSIaFTxHGamT5+Offv2YXBwEM8++yyWL1+O3bt3Y8aMGWhubjZfN3PmTFxxxRWor69HLBZDVVWV7XlDIetixtWrV2PVqlXmz0NDQ6ioqPB66TSGqF/iqxpqzFBQV1WK2srSpOeswoHVaii5FkXmtqZFfc6qJkX0mfG7QaV43fodBxA3hsNT3DA8BSIiolznOcwUFhaaBcCzZ89GV1cX1q9fjx/+8Icpr507dy4A4ODBg6iqqkJ5eTlefvnlpNccPXoUAFJGbGTjxo3DuHHjvF4qjXEt9dXoiPUh2jOQtNS5I9aPvYeOA4Djbs9W/V/U3bO9Fue6EVTX3sZNUXNEJm4YKTU0RET5Lu0+M4lEAmfOnNE+t2/fPgDAFVdcAQCIRCJYs2YNjh07hokTJwIA2tvbUVRUZE5VEQVpc3MENfdvM+tN2ppqk37e3BxxPIfdaigROPysknISRNdeudi3rak2pahYPh8RUb7yFGbuu+8+LFy4EBUVFTh58iTa2tqwa9cubN++HbFYDG1tbbj55ptRWlqK1157DStXrsRnPvMZzJo1CwBw4403YsaMGVi2bBkeffRR9Pb24utf/zruuusujryQKaiuveL1cuFv46ao62XSMqvVUFaPB/EZglgSrhb7ysGsrqqU/WGIaFTwtJrp6NGjWLZsGaZPn476+np0dXVh+/btaGhoQGFhIXbs2IEbb7wRV199Nb72ta/h1ltvxa9+9Svz+HA4jOeffx7hcBiRSARf/OIXcdtttyX1pSESdSrpNIuTX7+qoQbdaxaa00J1VaXoXrNQu6rI7lxqCLJ7PN3PEESzOzGyoxb7ilVOuuJlIqJ85Glk5oknnrB8rqKiArt373Y8x5VXXomtW7d6eVsaY/zuUyTT9WzR1bmo7+P2XOvau82CYrviXK+fQYzo6AqPhz9Dn7mtgBM2uiOisYJ7M1FOSrdrr1pvIv8sL512qj+xKsKVV0Y5Fed6+QwixIjNG+Xl5W5WXxERjUUhwzDybtJ8aGgIxcXFOHHiBIqKirJ9OZRBcrFu95qFI/a+YoQknjBSal827DyA3x3sQ7ggZNtnRjzu9TOom2DKgYZBhojyWaa+vzkyQzlFLpzVFe+OVJ2HPEJitdu2VahQX++14Fge0SkMXyhrG8ltGIiI8gnDDOUUESLkaZy2plpzibHgdVWTV5mq23G7HFpeJQXAV/dfIqKxgmGGcopcj1Jx+cXoiPUn9UqRf3a747NX6rJqueZFXc5stQRbBBe5KZ+XhndiREdYMW+a62OJiMYahhnKOXOmlgCAucGiboTGT+2I294v6koiEWTEtdRWlprHWW130BHrAwDUVZUlvY+bhnfyCI5cN6MWBTPQEBENY5ihnKMWzooQIX52CjJWoUUEj45YX1LnX7XgVh5Bifb0m9cgtgKQz5uykWPCQNehAXP5tFo/4zQ1Jo/o1FWVJR0vBxo2uyMiuoBhhnKSWjgbDoVSuuxasRotEcSO2XZ1MPJ0lxxk1B418nusa++2DD1qYLLiZgsDNrojIkrGpdmUc9SAIaaWRFBwM8VkVXwrT9uIoKQ7n9mp9/x7imXVdkXA8kiSfJ1Bb0BJRJSvuDSbclq6exGJ4wGkhJCOWD+mXH4x3jn+gesNEp0a1un2U5LFExdGYnQridRpHnUkSVxnJjagJCKiZAwzFAirqR230yvi+NrKEu1S5s/PrgAwHCJqK90HGl1ocdP7JVwQctyqQPcZ5dd6mRojIiL/GGYoEOn2ZUnezXl4BZBYEWR1vFMRrNUmkE69X6y2MFA/n91r5Xob9ochIsoshhlKm5e+LHa87GXkt2YGgGNIcVOEK+heK29qOWdqiWNBMhERpYdhhtLmpS+LE6upITes6m6Ev57wV/jfwdPa9wSA3x0cHgnS1fZY1f3o9mWyGo1ioCEiygyGGUqbPLrxs71vW/ZlcdtnxeteRoKu7kacU67b0Y0U2XXndVv3I87tdlSHiIiCwTBDgZDrRACk9GUB4BgI0tnLSH6NXHfjt25H/Ox1WbVdUOOIDBFRZjDMUCBEnYgYkQGA2srSpJVHdoHAqpC2I9Znu4pIHenxUnejk+7xREQ08gqyfQE0Ooi+LKLBHHAhEABI2nDR6nhdaFBXNgnm8ufzNTKylvpqc5rKz7Jor8e3tnebo0+qDTsPoPV8GCMiosxgmKFAyH1ZutcsNKeTRCCQ90LSWWmzamlVQ425BQGgH8WRA4Vad9O4KeopUFgt6bb77Ovau7Hk8c6k18qBi6GGiChzOM1EaXOqK0m3z4qbqR8RKETdjroVgt/P4qZuJ3ljygHzcXX7BDcFxERE5B3DDKVNnSJSV//Y1b243QbBacm2XIAsVlDJ/V50G0Tq3s9tszyV/Dq5Rkg8xrobIqLMYZihtMkFuFaBwGqEw+02CG6WbM+ZWgIA6Ij1m5s+yiMsYlm0VYAy636U5dNul1XLgQsAC4iJiEYIw8wYlu7mkDpe+6y4WQ7tdupHXKsIMvIIjpsAJep+dE3+3IQRNXBxXyYiopHBMDOGpbs5pI6fPit2NTFep37cjOCox8cTBroODSTV2sjncxPq1Hsmzs19mYiIMo9hZgwLoklckNeiq4lxGukRWxDorr1xU9Sy1kX+7LpuxYD7UKcLMvLfuY0BEVFmMcyMcZloEudn+spqRMXNSI9uFZNc/GsXaNR9pMT7egl1InCJa1GPsSuAJiKi9DHMEOIJA+FQSFvjYTfNYhVaxPRVR6wvqb+M1UhHOtsYJG9hUKoNInLxr/q+coASwcdrqBP3plUTZFrqq22vgYiI0scwQ+g6NIC4cSHQ6EYndKxqbgTR6M5upCOd5dCC/Hp1FZPV8VYByirUucF9mYiIsoNhZoxTe7GI0Ql12kbHqeZGPGc30hHULtNOfWjUz6wLUOIzq6GOiIhyG8PMGNXa3p2ygkd8yYeApOZzdpxqbuwCRpBLw92sYhLiCQO1lSUpx4vPnDi/USbrXIiI8gP3ZhqjRE8VObC01FcjHArBABDChSZ0Tlrq9RszOu1xJKap1MftNpHUkUdaxL5QuvMKKxtqUFdVZr5GPr62shTRngHUVZU5noeIiHIDR2bGKHlERYxiNG6KmrUzccPwFCZ0ocWpqNfL0nCrURzxenlXbnkFGMoAABtxSURBVDc1N/JraitLbFcjsXCXiCi3McyMYfIX+vodB8xeK21NtY4rilqlkROr0OKmqNft0nCrYuOO2HCfmbqqMu1nswsi8nv//i+D2vfmFBMRUe5jmBnjWuqrzSATDoXQ1lRrPg5Yj27I4cLtEma7LQ2cinetRnGiPQOW7+/2mtwWDhMRUW5izcwYt2HngaSpJbk+pKW+GqsaarSjGy311doiWjFKY3eceK18nDxNteTxTu21iusRS7Dl3jK6z9V6Pvg4fX67uh4iIsp9HJnJQ0GtAnLTrM5upGJzc8Q8xkujOXlUB0ht/y/3p1HZde3VfS43n99qas3vZptERDSyGGbyUBAbROoKbeMJw7L9v9UXu59pGnXKyGpPI/UaxHXouvaK17rdhkAOMmogUvvsEBFRbmOYyUNeVgFZ0TWrk5dry1NEdiHJS38X9TN0xPoQ7RmwHNVRp6msRpL8bEMgf351RMZNw0AiIsodIcMw8m7d6dDQEIqLi3HixAkUFRVl+3KyRnwJiyARxJevVWDQndvLa62I7QcKwwXoXrPQ9XWpj4uaH6fzOJ0/yHtJRETJMvX9zQLgPGbVrC7dc6pFtvIXe6um0Zz8vrWVJa4bzXkpvrXb9qCuqtQMMn6LeDNxL4mIaGQwzOSRViUkuF0F5JXdF7u8I7YcLi5M+ZRZrmSS+enaa1XwLKaE3JzH7nq4qomIKD+xZiaPyIW/ot5EtwoIQFqrcOzqYOR6HdGozuv0UhA7ZWfyPE4NA4mIKLcwzOQR9YsaQMqqG3WFkFdelms7Fd1aLSEXq6Z0zfPE824EseN2UIGIiIiyh2Emz8irgESfFdE4TozKyPsUqex61DRuiqas4rHbhsBpSbbVEnKxaqq2slT7+ayo1y6PPKlLx90GkCACERERZRfDTB7a3BwxVwGJQCN+dprqsetRo+6iLYgA1RHrS6qRUet1NjdHUo4DkoPQksc7LbcgcGpSF0R/HZXdVBxHZIiI8gPDzAgJqmuveL0cJMKhkOtVOH571NRVlSUV1rrt2qubkrL6TE6BJIj+OkRENPowzIyQoEYV1C9vMTUkAo2bpnVua16sjgG8de1Vp6RWzJvmO5CIUSKra+cWBEREYw/DzAhJZ1RBjOqI4+Xi3I5YP6ZcfjHeOf6B5VYEVtfjZxsCr117gdSRJADmEmovXXuFuqoyRHsGUq49nekmIiLKXwwzGSZPL+lGROyKdQUxqlNbWaJdZfT52RUAhoNEbaV1oJGvRQ0YjZuimDO1xHFEQ67XUUOQ1RSaVZfgIJrUyT1hON1ERDQ2McxkmDq9JI+IADB7tdjR9XbpiPUBgOWXt26ERFyLuveQmKoCnGt7OmJ9rvdicrPs2eueTuroy7r27qSpLwYZIqKxh2Emw+Qv745YHwrO17aonGo9vNS5WH2ht9Rf2ERRXs4tfhaBRvyvrrYHgHZ0yCpQ2QWM2soSbG6OuG5SpwtHdkXFREQ0NjDMjABds7vklUDJje/szuO1zkU1Z2oJAGiXc4tApU5V6YKM7nOp16MLZuk0qVPDkTpVJi8dJyKisYNhJot0IyV2dNsMxBOGpyXf4u+6uhf1HGptT11VWdrN5dJpUqc2ydPV4ridriIiotGDYWaEiBoXEUTkoFBXVWqOmFix+vKWp4fcLvm223tJUEeB1IZ46mvdCqJJHbcgICIiGcPMCNiw80BS11vxZSyCQltTrePxdl/e6pJsuyXfbjdVdBN4soVbEBARkYxhJsOcesk4BYXW9m7seavf8vjayhLMmVpi1rnYFQe7HdHI9V2kuQUBERHJGGZsBLEFga5oVZ7+Ed1sAf0XcbgghGjPQMoSbruVPVbFwW5GNDiFQ0RE+YZhxka6WxCoYUg9Lp4wHJcmu+0cbDctJK7DanWRHMpaLUaR0p3CCXJvKiIiIhnDjI10tiAAUsOQGBkR5xR/dwoKTj1mnKaFvISyTE3hZGLHayIiIoBhxlE8YZgFtmqQ8NLoDhgOClZhyCkoWPWY8TIt5DeUpTuqIo4X+zGp1+BmSwciIiIrDDMOwgWhpF2pRZBwO6LgpXOvHatpJLcre9K5jnRHVcTxqxpqUjaYBNxt6UBERGSFYcaB3NhOBBqxl5HbMJBu5167aSQv00J+ryPdkR35eHmDSYD7KRERUfoYZjTU3aXlvYtCgOuOvYKfni3iGgD9btC1lSWeVxel0zsm3REmNRAREREFpSDbF5CLxLRI46aoGSTammoRDoVgAGag2bDzgOO55BGM7jULzWkWp2PFNXTE+rTFvnVVZVjVUON6dZHf65C11FebQcjv3lBCYbjA1zUQERGpODKjoXbXbamvRuOmKOKGgXAohLgxXBRsFSScRlX+esJfWY6qqAW1IriI59xO7cjX0nVoIGVaTEyfeRndSWdkR54aU6eZ2L+GiIjSwZEZCy311VjVUIOOWD+qVm81p5Zia282HxeBReU0qnJl6aUALuzXJIjnxXnFNaxr70bN/ds8BxlxLbppMXn6zM3oTrojO+KzqseLx7gFARER+cUwY6OlvtociQmHQuYeSnLI0H2Zi+ejPQPmY3IYaGuqNZ8Xx9st2U5nakcOZbr3amuqdWxWZ7X8222gUfemUo8H7AuZiYiI7HCaSaEW/8pTS42bopgztQQrpS9lv43u3BbUBrHhY7rFu+lu7MiNIYmIKJM8jcxs3LgRs2bNQlFREYqKihCJRLBt27aU1xmGgYULFyIUCuG5555Leu7w4cNYtGgRLrnkEkycOBF33303zp07l96nCJCu+De29mZzNVPXoQujLS311Y5Lo+1GVZyeD6Jo1+176bSef6+VFiGrtb3b8R4A0B4vXxdHZYiIKB2ewsyUKVPwyCOPYO/evdi7dy/mzZuHxYsX4/XXX0963X/8x38gFEqtJ4nH41i0aBHOnj2Ljo4OPPXUU3jyySfxwAMPpPcpAtRSX20GF1Fnoi7PdhsmdKMqbp/3MrXTahNwROhwuhaVKBzWhScR9KxqhoiIiEaSp2mmW265JennNWvWYOPGjYhGo7jmmmsAAH/4wx+wbt06dHV14Yorrkh6/a9//Wu88cYb2LFjByZNmoRPfOIT+Na3voV77rkHDz74IAoLC9P8OMGYM7UEwPDy65r7t2m3MHDitF+S0/NepmacOvSKEGb1Xjpy4bD8WtEw0EufHSIiokzyXTMTj8fxzDPP4NSpU4hEIgCA999/H1/4whfw2GOPoby8POWYzs5OzJw5E5MmTTIfu+mmm3DnnXfi9ddfx/XXX699rzNnzuDMmTPmz0NDQ34v2xUx7SGCTGG4APGEYVmvoi6ndtovSXQUtnq+I9aHuqoy13sh2XXoVYOM1etV6vL0de3dWL/jgLksXRRDExERZZvnMLN//35EIhGcPn0a48ePx5YtWzBjxgwAwMqVK1FXV4fFixdrj+3t7U0KMgDMn3t7ey3fc+3atXjooYe8Xmpa1GkZ0asFcN6fyGlU5XcH+2yf74j1ed4LyarIN54wUFuZOoripvhWDT3qqi4iIqJcEDIMw9NSkrNnz+Lw4cMYHBzEs88+ix/96EfYvXs3Dh48iK997Wt49dVXMX78+OGTh0LYsmULPvvZzwIAmpub8Ze//AXbt283z/f+++/j0ksvxbZt27BgwQLte+pGZioqKnDixAkUFRV5/tBOrKaArKZrMrG/kNU1OL2XPJrUvWZhINdStXor4tKvCfdTIiIiP4aGhlBcXBz497fnkZnCwkJMmzYNADB79mx0dXVh/fr1uPjiixGLxTBhwoSk199666349Kc/jV27dqG8vBwvv/xy0vNHjx4FgJQRG9m4ceMwbtw4r5fqi9MUkZhy0S1xlpd1686rTg/Z8bOcWh1NatwU1Y6ieLkWXedjduwlIqJcknbTvEQigTNnzuDee+/Fa6+9hn379pl/AKC1tRX/9V//BQCIRCLYv38/jh07Zh7f3t6OoqIic6oq2+ymiFY11GDO1BLLJc6iEFe3asnP6h8vy6nVZdxiFKlxU9T3tcjFvnLnYxFouKcSERHlAk8jM/fddx8WLlyIiooKnDx5Em1tbdi1axe2b9+O8vJybdHvxz72MVx11VUAgBtvvBEzZszAsmXL8Oijj6K3txdf//rXcdddd43YyIsTp74xdk3s7Apx/UzNuG2Yp3uPtqZaM4yIERov1yIvR5c7H4vP53YbBCIiokzzFGaOHj2KZcuW4ciRIyguLsasWbOwfft2NDQ0uDo+HA7j+eefx5133olIJIJLL70Uy5cvx8MPP+zr4kea03Jq+X/9dtv18l6C1WiSHGjUJeZO3CwNZ7M7IiLKBZ4LgHNBpgqI7FiNalg9nk4hrtf3cpKJomAiIiKvcqYAeKzy0sQu3f2UgtzLKIi9nYiIiHIZw4xLTrU0gtvpIbuVT+GCkOMGlm54maoiIiLKVwwz5wWxrNppWbf8s24LArEfkuhl4+ca/FwLERFRPmOYOc9pfyNd111Za3s39rzVb1nPUltZkjTaogsWIsio+x65vQZZkFNVREREuYxh5rx0l1WHC0KI9gygrqos6XG7c+hWPsk7c6eztNvttBgREVG+Y5hB6hSTGi7cjGL4DUMt9dXmexWGC5L6wXhd2h1kB2IiIqJ8kXYH4NFA7twrd90Nh0LoiPW77twrugSva+9Gzf3bXI2oWK02ctv5V/c5rLr+dh0aQOv5sEVERDRaMMwgOYQ0boqaQUbsReRlWiadLQjUa5ADjtv3VrcxUDfJ9LqlAhERUa5jmDlPDgJykBH1K0KrzZ5EG3YewJLHO10FEavVRuI966pKkwKO20DT1lRrnqNq9Vbtbt9ERESjCWtmzhN7EYkgo9avAMNhw2nVEwDfWxDI+yHNmVqSdJyX5dRtTbWoWr0V8fPNnRlkiIhoNGOYOS+euDASo+uWK4qA7Qp9Abju66IrxA1qOfWGnQfMIAMA4ZC+KJiIiGg0YJg5L1wQShrBsBtV0a16qq0sQV1VWVpBJIjl1Ck1MudHmsTO2URERKMNwwz8dcuVl1SHQyFsbo5ozzuSy6HVICM+j9g5m4GGiIhGIxYAw356Z1VDjXZURSyplkc+1OfXne/7MlLkqTL588hFwW4LiYmIiPJFyDCMvOtrn6ktxN1SR3LEyEddVWlS0bCXotugGt6xcR4REeWqTH1/M8x4ZBVURKARIzVeu/YCSDnv8OqmPkR7BrgaiYiI8l6mvr9ZM+OR1ZSUvBzaT9feVQ01Zk8ZQbdCioiIiJIxzHhkNUUjlkPrlnXbkQuNdYGGQYaIiMgew0wA1Kknu2XdOnKgKQxfqMn2MsJDREQ0VnE1U5qslnV73YZA3tMJgOd9mYiIiMYqjsykobW9G3vest4qoLayxFPXXhFkAGDFvGkAvG1jQERENBYxzPggViCFC0KI9gygrqrMfM7PCiR1OwRAX0PDQENERJSKYcaHIFcgiSCj2w5Bfg+3IzxERERjDcOMD0GuQHKzuSSb3BEREVlj07w0iFEVtXC3e83CrF0TERFRrsrU9zdXM6WBK5CIiIiyj9NMaeAKJCIiouzjyIxP6gokuW7Ga48ZIiIi8o8jMz5wBRIREVHuYJjxgSuQiIiIcgdXMxEREdGI4GomIiIiIg2GGSIiIsprDDNERESU1xhmiIiIKK8xzBAREVFeY5ghIiKivMYwQ0RERHmNYYaIiIjyGsMMERER5TWGGSIiIsprebk3k9iBYWhoKMtXQkRERG6J7+2gd1LKyzBz8uRJAEBFRUWWr4SIiIi8OnnyJIqLiwM7X15uNJlIJPDuu+/isssuQygUCvTcQ0NDqKiowNtvv81NLD3gffOP984f3jf/eO/84X3zT9y7w4cPIxQKYfLkySgoCK7SJS9HZgoKCjBlypSMvkdRURF/WX3gffOP984f3jf/eO/84X3zr7i4OCP3jgXARERElNcYZoiIiCivhR988MEHs30RuSYcDuOGG27ARRfl5Sxc1vC++cd75w/vm3+8d/7wvvmXyXuXlwXARERERAKnmYiIiCivMcwQERFRXmOYISIiorzGMENERER5jWGGiIiI8hrDjOT73/8+pk6dir/6q7/C3Llz8fLLL2f7knLOgw8+iFAolPTn6quvNp8/ffo07rrrLpSWlmL8+PG49dZbcfTo0SxecXa89NJLuOWWWzB58mSEQiE899xzSc8bhoEHHngAV1xxBS6++GLMnz8fBw4cSHrNwMAAli5diqKiIkyYMAG333473nvvvZH8GFnhdO/++Z//OeV3cMGCBUmvGYv3bu3atZgzZw4uu+wyTJw4EZ/97Gfx5ptvJr3Gzb/Pw4cPY9GiRbjkkkswceJE3H333Th37txIfpQR5ea+3XDDDSm/c3fccUfSa8bafQOAjRs3YtasWWZH5Egkgm3btpnPj+TvG8PMeT/96U+xatUqfPOb38Tvf/97XHfddbjppptw7NixbF9azrnmmmtw5MgR889vf/tb87mVK1fiV7/6FZ555hns3r0b7777Lj73uc9l8Wqz49SpU7juuuvw2GOPaZ9/9NFHsWHDBmzcuBF79uzBpZdeiptuugmnT582X7N06VK8/vrraG9vx/PPP4+XXnoJzc3NI/URssbp3gHAggULkn4Hf/KTnyQ9Pxbv3e7du3HXXXchGo2ivb0dH374IW688UacOnXKfI3Tv894PI5Fixbh7Nmz6OjowFNPPYUnn3wSDzzwQDY+0ohwc98AoKmpKel37tFHHzWfG4v3DQCmTJmCRx55BHv37sXevXsxb948LF68GK+//jqAEf59M8gwDMP4m7/5G+Ouu+4yf47H48bkyZONtWvXZvGqcs83v/lN47rrrtM+Nzg4aHzkIx8xnnnmGfOxP/3pTwYAo7Ozc6QuMecAMLZs2WL+nEgkjPLycuM73/mO+djg4KAxbtw44yc/+YlhGIbxxhtvGACMrq4u8zXbtm0zQqGQ8b//+78jd/FZpt47wzCM5cuXG4sXL7Y8hvdu2LFjxwwAxu7duw3DcPfvc+vWrUZBQYHR29trvmbjxo1GUVGRcebMmZH9AFmi3jfDMIy///u/N77yla9YHsP7dsHll19u/OhHPxrx3zeOzAA4e/YsXnnlFcyfP998rKCgAPPnz0dnZ2cWryw3HThwAJMnT0ZlZSWWLl2Kw4cPAwBeeeUVfPjhh0n38eqrr8bHPvYx3kfJW2+9hd7e3qT7VFxcjLlz55r3qbOzExMmTMDs2bPN18yfPx8FBQXYs2fPiF9zrtm1axcmTpyI6dOn484770R/f7/5HO/dsBMnTgAASkpKALj799nZ2YmZM2di0qRJ5mtuuukmDA0Nmf+1Pdqp9014+umnUVZWhmuvvRarV6/G+++/bz7H+zY8yrJ582acOnUKkUhkxH/f2I8ZQF9fH+LxeNINBYBJkybhz3/+c5auKjfNnTsXTz75JKZPn44jR47goYcewqc//Wn88Y9/RG9vLwoLCzFhwoSkYyZNmoTe3t4sXXHuEfdC9/smnuvt7cXEiROTnr/oootQUlIy5u/lggUL8LnPfQ5XXXUVYrEY7rvvPixcuBCdnZ0Ih8O8dwASiQS++tWv4m//9m9x7bXXAoCrf5+9vb3a30vx3Ginu28A0NjYiCuvvBKTJ0/Ga6+9hnvuuQdvvvkmfvGLXwAY2/dt//79iEQiOH36NMaPH48tW7ZgxowZ2Ldv34j+vjHM2DAMA6FQKNuXkVMWLlxo/n3WrFmYO3currzySvzsZz/DxRdfrD2G99EdwzBQUGA/WMp7CSxZssT8+8yZMzFr1ixUVVVh165dqK+vtzxuLN27u+66C3/84x+T6tmsuL0vY+HeWd03ud5q5syZuOKKK1BfX49YLIaqqirbc472+zZ9+nTs27cPg4ODePbZZ7F8+XLs3r3b8vWZ+n3jNBOAsrIyhMPhlCrrY8eOpaRGSjZhwgTU1NTg4MGDKC8vx9mzZzE4OJj0Gt7HZOXl5QBg+/tWXl6eUnx+7tw5HD9+nPdSUVlZibKyMhw8eBAA792KFSvw/PPP43/+538wZcoU83E3/z7Ly8tTfi/Fz6P93lndN525c+cCQNLv3Fi9b4WFhZg2bRpmz56NtWvX4rrrrsP69etH/PeNYQbD/8f41Kc+hZ07d5qPJRIJ7Ny5E5FIJItXlvvee+89xGIxXHHFFfjUpz6Fj3zkI0n3sbu7G4cPH+Z9lFx11VUoLy9Puk9DQ0PYs2ePeZ8ikQgGBwfxyiuvmK958cUXkUgkzP9HSsPeeecd9Pf344orrgAwdu+dYRhYsWIFtmzZghdffBFXXXVV0vNu/n1GIhHs378/KQy2t7ejqKgIM2bMGJkPMsKc7pvOvn37ACDpd26s3TcriUQCZ86cGfnft7TKlkeRzZs3G4WFhcaTTz5pvPHGG0Zzc7MxYcKEpCprMoyvfe1rxq5du4y33nrL+N3vfmfMnz/fKCsrM44dO2YYhmHccccdxsc+9jHjxRdfNPbu3WtEIhEjEolk+apH3smTJ41XX33VePXVVw0Axrp164xXX33V+Mtf/mIYhmE88sgjxoQJE4xf/vKXxmuvvWYsXrzYuOqqq4wPPvjAPMeCBQuM66+/3tizZ4/x29/+1qiurja+8IUvZOsjjRi7e3fy5Enj3/7t34zOzk7jrbfeMnbs2GF88pOfNKqrq43Tp0+b5xiL9+7OO+80iouLjV27dhlHjhwx/7z//vvma5z+fZ47d8649tprjRtvvNHYt2+f8cILLxgf/ehHjdWrV2fjI40Ip/t28OBB4+GHHzb27t1rvPXWW8Yvf/lLo7Ky0vjMZz5jnmMs3jfDMIzVq1cbL730kvHWW28Zr732mnHvvfcaoVDI+PWvf20Yxsj+vjHMSP7/9u0eRWEggOJ4tnBCJAhKgkjANJ7B0iYgWFoFK1GwsQ2WnsIDeBRbryBYiI2VhQhaWLwtFgRx2Y8mcfD/g1QJw8xjAi8hWSwWajabMsao3W5rvV4XPaWXk6apGo2GjDGKokhpmmq73d7PX69XTadTVatVlctl9ft9HQ6HAmdcjNVqJcdxno7hcCjp6/fs+Xyuer0u13WVJIk2m83DGMfjUYPBQL7vq1KpaDQa6Xw+F7CafP2U3eVyUbfbVRiGKpVKiuNYk8nk6aHjHbP7LjPHcbRcLu/X/OX+3O126vV68jxPQRAoyzLdbrecV5Of33Lb7/fqdDqq1WpyXVetVkuz2Uyn0+lhnHfLTZLG47HiOJYxRmEYKkmSe5GR8t1vH5L0v3c5AAAAr4NvZgAAgNUoMwAAwGqUGQAAYDXKDAAAsBplBgAAWI0yAwAArEaZAQAAVqPMAAAAq1FmAACA1SgzAADAapQZAABgtU8qiLSwj5jh0AAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x137a24350>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x137d6fe10>"
      ]
     },
     "execution_count": 277,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q=maunaloa[:,2:13]'\n",
    "\n",
    "b=Q[:] \n",
    "\n",
    "using PyPlot\n",
    "\n",
    "plot(b,\"x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the data seems quite oscillatory.\n",
    "A suitable function to fit to the given data could be\n",
    "$$y(t)=c_1+c_2e^{\\alpha t} + \\sum_{k=1}^{np} \\left(a_k cos({2\\pi kt\\over 12}) +\n",
    "b_k sin({2\\pi kt\\over 12})\\right)$$\n",
    "where $t$ is in months, $\\alpha = 0.00037$. \n",
    "**The other parameters to be computed with the least squares method.** \n",
    "\n",
    "1. Solve the least square problem for $np=1,2,...,12$ by solving directly the normal equation with the backlash operator. This means `x=(A'A)\\(A'b)`. Report also $cond(A'A)$.\n",
    "\n",
    "2. Use the function `qrfact()` to solve the least squares problem for $np=1,2,\\cdots,12$, truncating $A$ to rank $r$, for $r=1,2,\\cdots,14$. This means, following the notation in teh previous discussion, $R_{11}$ is is an $r\\times r$ matrix . Report also $cond(R_{11})$.\n",
    "\n",
    "3. Present the computed values of the parameters in each case. Plot out the solutions $y(t)$ against the data points.\n",
    "Compare your numerical solutions computed in 1. and in 2. Try to explain the discrepancies that you observe.\n",
    "\n",
    "4. Compute also the residual and write out its 2-norm in each case. Did the residual always decrease as $np$ increases?\n",
    "\n",
    "5. Perturb the data by a small amount: ` b.+=1e-5*rand(size(b))`, and compare your least square solutions with the ones that you computed without the perturbation. \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}