lessons/notebooks/06_Linear-Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Preliminaries\n",
    "\n",
    "- Goal \n",
    "  - Maximum likelihood estimates for various linear regression variants\n",
    "- Materials        \n",
    "  - Mandatory\n",
    "    - These lecture notes\n",
    "  - Optional\n",
    "    - Bishop pp. 140-144 \n",
    "    - [G. Deng et al., _A model-based approach for the development of LMS algorithms_', ISCAS-05 symposium, 2005](./files/Deng-2005-A-model-based-approach-for-the-development-of-LMS-algorithms.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Regression - Illustration\n",
    "\n",
    "\n",
    "<img src=\"./figures/fig-bishop12.png\" width=\"300px\">\n",
    "\n",
    "Given a set of (noisy) data measurements, find the 'best' relation between an input variable $x \\in \\mathbb{R}^D$ and input-dependent outcomes $y \\in \\mathbb{R}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Regression vs Density Estimation\n",
    "\n",
    "\n",
    "- Observe $N$ IID data **pairs** $D=\\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$ with $x_n \\in \\mathbb{R}^D$ and $y_n \\in \\mathbb{R}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- <span style=\"color:blue\">[Q.]</span> We could try to build a model for the data by density estimation, $p(x,y)$, but what if we are interested only in (a model for) the responses $y_n$ for **given inputs** $x_n$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- <span style=\"color:blue\">[A.]</span> We will build a model only for the conditional distribution $p(y|x)$. \n",
    "  - Note that, since $p(x,y)=p(y|x)\\, p(x)$, this is a building block for the joint data density.\n",
    "  - In a sense, this is density modeling with the assumption that $x$ is drawn from a uniform distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Next, we discuss model (1) specification, (2) ML estimation and (3) prediction for the linear regression model. \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  1. Model Specification for Linear Regression\n",
    "\n",
    "\n",
    "- In a _regression_ model, we try to 'explain the data' by a purely deterministic term $f(x,w)$, plus a purely random term $\\epsilon_n$ for 'unexplained noise',\n",
    "\n",
    "    $$\n",
    "    y_n  = f(x_n,w) + \\epsilon_n\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  In _linear regression_, we assume that \n",
    "\n",
    "$$f(x,w)=w^T x \\,.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  In _ordinary linear regression_, the noise process $\\epsilon_n$ is zero-mean Gaussian with constant variance $\\sigma^2$, i.e.\n",
    "\n",
    "$$\n",
    "y_n  = w^T x_n  + \\mathcal{N}(0,\\sigma^2) \\,,\n",
    "$$\n",
    "or equivalently, the likelihood model is \n",
    "$$\n",
    "p(y_n|\\,x_n,w) = \\mathcal{N}(y_n|\\,w^T x_n,\\sigma^2) \\,.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- For full Bayesian learning we should also choose a prior $p(w)$; In ML estimation, the prior $p(w)$ is uniformly distributed (so it can be ignored).\n",
    "\n",
    "<!--- \n",
    "%$y_n$ is Gaussian around $\\theta^T x_n$ with `unexplained' noise $\\sigma^2$.\n",
    "% -   Note that 1st term RHS is purely deterministic; 2nd term is purely random and independent of $x_n$\n",
    "%-  Choose a prior over the parameters (possibly flat: $\\sigma\\rightarrow \\infty$)\n",
    "-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  2. ML Estimation for Linear Regression Model\n",
    "\n",
    "- Let's work out the log-likelihood for multiple observations\n",
    "$$\\begin{align*}\n",
    "\\log p(D|w) &\\stackrel{\\text{IID}}{=} \\sum_n \\log \\mathcal{N}(y_n|\\,w^T x_n,\\sigma^2) \\propto -\\frac{1}{2\\sigma^2} \\sum_{n} {(y_n - w^T x_n)^2}\\\\\n",
    "    &= -\\frac{1}{2\\sigma^2}\\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X} w } \\right)\n",
    "\\end{align*}$$\n",
    "where  we defined $N\\times 1$ vector $y  = \\left(y_1 ,y_2 , \\ldots ,y_N \\right)^T$ and $(N\\times D)$-dim matrix $\\mathbf{X}  = \\left( x_1 ,x_2 , \\ldots ,x_n \\right)^T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Set the derivative $\\nabla_{w} \\log p(D|w) = \\frac{1}{\\sigma^2} \\mathbf{X}^T(y-\\mathbf{X} w)$ to zero for\n",
    "the maximum likelihood estimate\n",
    "$$\\begin{equation*}\n",
    "\\boxed{\\hat w_{\\text{ML}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y}\n",
    "\\end{equation*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  The matrix $\\mathbf{X}^\\dagger \\equiv  (\\mathbf{X}^T \\mathbf{X})^{-1}\\mathbf{X}^T$ is also known as the **Moore-Penrose pseudo-inverse** (which is sort-of-an-inverse for non-square matrices)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Note that size ($N\\times D$) of the data matrix $\\mathbf{X}$ grows with number of observations, but the size ($D\\times D$) of $\\mathbf{X}^T\\mathbf{X}$ is independent of training data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 3. Prediction of New Data Points \n",
    "\n",
    "- Now, we want to apply the trained model. New data points can be predicted by\n",
    "$$\\begin{equation*}\n",
    "p(y_\\bullet \\,|\\, x_\\bullet,\\hat w_{\\text{ML}}) = \\mathcal{N}(y_\\bullet \\,|\\, \\hat w_{\\text{ML}}^T x_\\bullet, \\sigma^2 ) \n",
    "\\end{equation*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Note that the expected value of a predicted new data point\n",
    "\n",
    "$$\n",
    "\\mathrm{E}[y_\\bullet] = \\hat w_{\\text{ML}}^T x_\\bullet = x_\\bullet^T \\hat{w}_{\\text{ML}} = \\left( x_\\bullet^T \\mathbf{X}^\\dagger \\right) y\n",
    "$$\n",
    "\n",
    "can also be expressed as a linear combination of the observed data points \n",
    "\n",
    "$$y  = \\left( {y_1 ,y_1 , \\ldots ,y_N } \\right)^T \\,.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Deterministic Least-Squares Regression\n",
    "\n",
    "-  (You may say that) we don't need to work with probabilistic models. E.g., there's also the deterministic **least-squares** solution: minimize sum of squared errors,\n",
    "$$\\begin{align*} \\hat w_{\\text{LS}} &= \\arg\\min_{w} \\sum_n {\\left( {y_n  - w ^T x_n } \\right)} ^2 \n",
    "  = \\arg\\min_{w} \\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X} w } \\right)\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Setting the gradient \n",
    "$ \\frac{\\partial \\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X}w } \\right)}{\\partial w} = -2 \\mathbf{X}^T \\left(y - \\mathbf{X} w  \\right)\n",
    "$ to zero yields the **normal equations** \n",
    "$\\mathbf{X}^T\\mathbf{X} \\hat w_{\\text{LS}} = \\mathbf{X}^T y$  and consequently\n",
    "$$\n",
    "\\boxed{\\hat w_{\\text{LS}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y} \n",
    "$$\n",
    "which is the same answer as we got for the maximum likelihood weights $\\hat w_{\\text{ML}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- $\\Rightarrow$ Least-squares regression ($\\hat w_{\\text{LS}}$) corresponds to (probabilistic) maximum likelihood ($\\hat w_{\\text{ML}}$) if \n",
    "  1. **IID samples** (determines how errors are combined), and\n",
    "  1.  Noise $\\epsilon_n \\sim \\mathcal{N}(0,\\,\\sigma^2)$ is **Gaussian** (determines error metric)      \n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Probabilistic vs. Deterministic Approach\n",
    "\n",
    "-  The (deterministic) least-squares approach assumed IID Gaussian distributed data, but these assumptions are not obvious from looking at the least-squares (LS) criterion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  If the data were better modeled by non-Gaussian assumptions (or not IID), then LS might not be appropriate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  The probabilistic approach makes all these issues completely transparent by focusing on the **model specification** rather than the error criterion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Next, we will show this by two examples: (1) samples not identically distributed, and (2) few data points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Not Identically Distributed Data\n",
    "\n",
    "-  What if we assume that the variance of the measurement error varies with the sampling index,  $\\epsilon_n \\sim \\mathcal{N}(0,\\sigma_n^2)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Let's make the log-likelihood again (use $\\Lambda \\triangleq \\mathrm{diag}[1/\\sigma_n^2]$): \n",
    "$$\\begin{align*}\n",
    "\\mathrm{L(w)} &\\triangleq \\log p(D|w) \n",
    "  \\propto -\\frac{1}{2} \\sum_n \\frac{(y_n-w^T x_n)^2}{\\sigma_n^2} = -\\frac{1}{2} (y- \\mathbf{X}w)^T \\Lambda (y- \\mathbf{X} w)\\,.\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Set derivative $\\partial \\mathrm{L(w)} / \\partial w = -\\mathbf{X}^T\\Lambda (y-\\mathbf{X} w)$\n",
    "to zero to get the **normal equations** \n",
    "$\\mathbf{X}^T \\Lambda \\mathbf{X} \\hat{w}_{\\text{WLS}} = \\mathbf{X}^T \\Lambda y$ \n",
    "and consequently \n",
    "$$ \\boxed{\\hat{w}_{\\text{WLS}} = \\left(\\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \\mathbf{X}^T \\Lambda y}$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  This is also called the **Weighted Least Squares** (WLS) solution. (Note that we just stumbled upon it, the crucial aspect is appropriate model specification!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Note also that the dimension of $\\Lambda$ grows with the number of data points. In general, models for which the number of parameters grow as the number of observations increase are called **non-parametric models**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### CODE EXAMPLE\n",
    "\n",
    "We'll compare the Least Squares and Weighted Least Squares solutions for a simple linear regression model with input-dependent noise:\n",
    "\n",
    "$$\\begin{align*}\n",
    "x &\\sim \\text{Unif}[0,1]\\\\\n",
    "y|x &\\sim \\mathcal{N}(f(x), v(x))\\\\\n",
    "f(x) &= 5x - 2\\\\\n",
    "v(x) &= 10e^{2x^2}-9.5\\\\\n",
    "\\mathcal{D} &= \\{(x_1,y_1),\\ldots,(x_N,y_N)\\}\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAG2CAYAAACd5Zf9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3XtcVHX+P/DXMDDcBLwioKiId03FUNNUwLzkJWF2t7LV0rbbllnaxUwXRMbiW+uudvdnF93WWmu3Bs3UTAUlL5v3TF3ygnkD8TrIRQZmPr8/pjlwYIAB5s7r+XjMY5s5Z868QZbz4nPen89RCCEEiIiIiDyMl7MLICIiIrIHhhwiIiLySAw5RERE5JEYcoiIiMgjMeQQERGRR2LIISIiIo/EkENEREQeiSGHiIiIPBJDDhEREXkkhhwiIiLySG4TcioqKvCXv/wFUVFR8Pf3R9euXZGWlgaj0ejs0oiIiMgFeTu7AGu98cYbWLFiBf7xj3+gb9++2L9/Px599FGEhITg+eefd3Z5RERE5GLcJuTs2bMHiYmJmDRpEgCgS5cu+Ne//oX9+/c7uTIiIiJyRW4TckaMGIEVK1bgl19+QY8ePXDkyBH88MMPWL58ea3vKSsrQ1lZmfTcaDTi+vXraNOmDRQKhSPKJiIioiYSQuDWrVuIiIiAl1cDOm2EmzAajWL+/PlCoVAIb29voVAoxOuvv17nexYtWiQA8MEHH3zwwQcfHvA4f/58g7KDQggh4AbWrl2Ll19+GX/961/Rt29fHD58GHPmzMHf//53zJgxw+J7qo/k6HQ6dOrUCefPn0dwcLCjSiciIqImKCwsRGRkJG7evImQkBCr3+c2IScyMhLz58/HrFmzpNeWLFmCNWvW4H//+59VxygsLERISAh0Oh1DDhERkZto7PnbbaaQl5SU1LgOp1QqOYWciIiILHKbxuP77rsPr732Gjp16oS+ffvi0KFD+Pvf/44//elPzi6NiIiIXJDbXK66desWkpOTodVqUVBQgIiICDz00ENISUmBSqWy6hi8XEVEROR+Gnv+dpuQYwvWfJOMRiP0er2DK6PmzsfHB0ql0tllEBG5pMaGHLe5XOUIer0eubm57PMhp2jZsiXCwsK4hhMRkY0w5PxGCIG8vDwolUpERkY2bLEhoiYQQqCkpAQFBQUAgPDwcCdXRETkGRhyflNRUYGSkhJEREQgICDA2eVQM+Pv7w8AKCgoQGhoKC9dERHZAIcrfmMwGADA6iZmIlszh+vy8nInV0JE5BkYcqphPwQ5C3/2iIhsiyGHiIiIPBJDjpsTQuDJJ59E69atoVAocPjwYVy7dg2hoaE4e/asVccoKytDp06dcODAAfsWS0RE5EAMOW5u8+bNWL16NTZs2IC8vDz069cP6enpuO+++9ClSxerjuHr64uXXnoJr7zyin2LJSIiciCGHDd3+vRphIeHY/jw4QgLC0N5eTk+/vhjPP744w06zrRp05CdnY0TJ07YqVIiIiLHYshxYzNnzsTs2bNx7tw5KBQKdOnSBZs2bYK3tzeGDRsm7ZeWloaIiAhcu3ZNem3KlCkYNWqUtPBhmzZtMHz4cPzrX/9y+NdBRERkD1wnpx7FxcW1blMqlfDz87NqXy8vL2ktlLr2DQwMtLq2t956C9HR0Vi5ciX27dsHpVKJ1157DbGxsbL9Fi5ciM2bN+Pxxx+HVqvFihUrsHPnThw5ckS26OGQIUOQnZ1t9ecTERG5MoacerRo0aLWbRMnTsS3334rPQ8NDUVJSYnFfePi4pCVlSU979KlC65evVpjv4bcSiwkJARBQUFQKpUICwsDAJw9exYRERGy/ZRKJdasWYOBAwdi/vz5eOedd7By5Up07txZtl+HDh2sblYmIiJydQw5Hqa0tFQ2umTWtWtXLF26FE899RQefPBBTJs2rcY+/v7+tYY0IiIid8OQU4+ioqJat1Vfet987yFLqt8Ly14jJm3btsWNGzcsbtu5cyeUSiXOnj2LiooKeHvL//mvX7+Odu3a2aUuIiIiR2PjcT0CAwNrfVQfMalr36r9OHXt21QxMTE4fvx4jde/+OILfP3118jKysL58+eh0Whq7PPzzz8jJiamyTUQERG5AoYcDzN+/HgcO3ZMNppz4cIFPP3003jjjTcwYsQIrF69Gunp6di7d6/svdnZ2Rg3bpyjSyYiIrILhhwPc8cddyA2NhZffvklAFMj88yZMzFkyBA8++yzAICxY8fi2WefxfTp06XLcXv27IFOp8Mf/vAHp9VORERkSwrRkOk8bq6wsBAhISHQ6XQIDg6Wbbt9+zZyc3MRFRVlsXHXnWzcuBEvvfQSfv755xq9QLW5//77ERMTgwULFti5OqqNJ/0MEhHZUl3n77qw8dgDTZw4ESdPnsTFixcRGRlZ7/5lZWUYMGAA5s6d64DqiIjIHlJTU6FUKpGcnFxjm0ajgcFgQGpqquMLcyJervJQzz//vFUBBzDdu+ovf/lLjeZoIiJyH0qlEikpKTUmlmg0GqSkpNSYEdwccCSHiIjIA5hHcFJSUqTn5oCTlpZmcYTH0zHkEBEReYiqQWfJkiXQ6/XNNuAAvFxFRETkUZKTk6FSqaDX66FSqZptwAEYcoiIiDyKRqORAo5er7e4+GtzwZBDRETkIar24JSVlSEtLc1iM3JzwZ4cIiIiD2CpydhSM3JzwpBDRETkAQwGg8UmY/Nzg8HgjLKciiGHiIjIA9S10F9zG8ExY0+Om5s5cyYUCgUUCgV8fHzQvn17jB07Fp988gmMRqOzyyMiInIahhwbSU1NrbWxS6PR2HUp7XvvvRd5eXk4e/YsNm3ahISEBDz//POYPHkyKioq7Pa5RERErowhx0acuZy2r68vwsLC0KFDBwwaNAgLFizAunXrsGnTJqxevdpun0tEROTK2JNjI662nPbo0aMxYMAAfP3113j88ccd+tlERESugCHHhlxtOe1evXrhp59+cspnExERORsvV9mYKy2nLYSAQqFw2ucTERE5E0OOjbnSctonTpxAVFSU0z6fiIjImRhybMiVltPevn07jh49it///vcO/2wiIiJXwJ4cG3HmctplZWXIz8+HwWDA5cuXsXnzZqSnp2Py5Ml45JFH7PKZREREro4hx0acuZz25s2bER4eDm9vb7Rq1QoDBgzA22+/jRkzZsDLi4N1RETUPCmEEMLZRThKYWEhQkJCoNPpEBwcLNt2+/Zt5ObmIioqCn5+fk6qkJoz/gwSEVlW1/m7Lvwzn4iIiDwSQw4RERF5JIYcIiIi8kgMOUREROSRGHKIiIjIIzHkEBERkUdiyCEiIiKP5FYh5+LFi5g+fTratGmDgIAADBw4EAcOHHB2WUREROSC3GbF4xs3buDuu+9GQkICNm3ahNDQUJw+fRotW7Z0dmlERETkgtxmJOeNN95AZGQkVq1ahSFDhqBLly645557EB0d7ezSPE6XLl2wfPly6blCoUBGRoYTKyIiImo4twk569evR2xsLO6//36EhoYiJiYGH374YZ3vKSsrQ2FhoezhaWbOnImkpKRatx86dAiTJ09GaGgo/Pz80KVLFzz44IO4evWq1Z+Rl5eHCRMm2KJcIiIih3GbkHPmzBl88MEH6N69O7777jv8+c9/xnPPPYdPP/201vekp6cjJCREekRGRjqwYucrKCjAmDFj0LZtW3z33Xc4ceIEPvnkE4SHh6OkpMTq44SFhcHX19eOldZPCIGKigqr9i0vL7dbHfY8NhER2ZbbhByj0YhBgwbh9ddfR0xMDJ566ik88cQT+OCDD2p9z6uvvgqdTic9zp8/78CKnW/37t0oLCzERx99hJiYGERFRWH06NFYvnw5OnXqZPVxql6uOnv2LBQKBb7++mskJCQgICAAAwYMwJ49e2p89qhRo+Dv74/IyEg899xzKC4ulravWbMGsbGxCAoKQlhYGP74xz+ioKBA2p6VlQWFQoHvvvsOsbGx8PX1RXZ2do3azPV8+eWXiI+Ph5+fH9asWWNVDXl5eZg0aRL8/f0RFRWFzz//3OKluhUrViAxMRGBgYFYsmSJ1d83IiJyLrcJOeHh4ejTp4/std69e+PcuXO1vsfX1xfBwcGyh7WEAIqLnfOw1X3hw8LCUFFRAa1WC1vfbH7hwoV46aWXcPjwYfTo0QMPPfSQNNJy9OhRjB8/Hr/73e/w008/4YsvvsAPP/yAZ599Vnq/Xq+HRqPBkSNHkJGRgdzcXMycObPG58ybNw/p6ek4ceIE+vfvX2s9r7zyCp577jmcOHEC48ePt6qGRx55BJcuXUJWVha++uorrFy5Uha0zBYtWoTExEQcPXoUf/rTn5rwXSMiIocSbuKhhx4SI0aMkL02Z84cMWzYMKuPodPpBACh0+lqbCstLRXHjx8XpaWlQgghioqEMMUNxz+Kiqz/vsyYMUMkJibWun3BggXC29tbtG7dWtx7773izTffFPn5+XUes3PnzmLZsmXScwBCq9UKIYTIzc0VAMRHH30kbT927JgAIE6cOCGEEOLhhx8WTz75pOyY2dnZwsvLS/r+Vvfjjz8KAOLWrVtCCCEyMzMFAJGRkVFnreZ6li9fLnu9vhpOnDghAIh9+/ZJ20+ePCkA1Pja58yZU2cNtlL9Z5CIiEzqOn/XxW1GcubOnYu9e/fi9ddfx6lTp/D5559j5cqVmDVrlrNLc2mvvfYa8vPzsWLFCvTp0wcrVqxAr169cPTo0SYdt+qoSnh4OABIoyAHDhzA6tWr0aJFC+kxfvx4GI1G5ObmAjA1RCcmJqJz584ICgpCfHw8ANQYmYuNjbWqnur71VdDTk4OvL29MWjQIOk93bp1Q6tWreo9NhERuQe3WSdn8ODB0Gq1ePXVV5GWloaoqCgsX74c06ZNs8vnBQQARUV2ObRVn21Lbdq0wf3334/7778f6enpiImJwdKlS/GPf/yj0cf08fGR/luhUAAw9U2Z//epp57Cc889V+N9nTp1QnFxMcaNG4dx48ZhzZo1aNeuHc6dO4fx48dDr9fL9g8MDLSqnur71VdDTk6OxeMIC5f1rK2BiIhci9uEHACYPHkyJk+e7JDPUigATzy3qVQqREdHyxpwbW3QoEE4duwYunXrZnH70aNHcfXqVfzf//2fNONt//79Dq2hV69eqKiowKFDh3DnnXcCAE6dOoWbN2/atA4iInIet7lcRbXT6XQ4fPiw7HHu3Dls2LAB06dPx4YNG/DLL78gJycHS5cuxcaNG5GYmGi3el555RXs2bMHs2bNwuHDh3Hy5EmsX78es2fPBmAaSVGpVHjnnXdw5swZrF+/HhqNxqE19OrVC2PGjMGTTz6JH3/8EYcOHcKTTz4Jf39/aWSqNgaDAb169cI333xj05qJiMi23GokhyzLyspCTEyM7LUZM2YgJSUFAQEBePHFF3H+/Hn4+vqie/fu+Oijj/Dwww/brZ7+/ftjx44dWLhwIUaOHAkhBKKjo/Hggw8CANq1a4fVq1djwYIFePvttzFo0CAsXboUU6ZMcVgNAPDpp5/isccew6hRoxAWFob09HQcO3YMfn5+dR5bCIGcnBzodDqb1UtERLanEJaaEDxUYWEhQkJCoNPpakwnv337NnJzcxEVFVXvSY4804ULFxAZGYmtW7finnvucfjn82eQiMiyus7fdeFIDjVb27dvR1FREe644w7k5eVh3rx56NKlC0aNGuXs0oiIyAYYcqjZKi8vx4IFC3DmzBkEBQVh+PDh+Oyzz2Qzx4jcVWpqKpRKJZKTk2ts02g0MBgMSE1NdXxhRA7ExmNqtsaPH4+ff/4ZJSUluHz5MrRaLTp37uzssohsQqlUIiUlpUZTv0ajQUpKCpRKpZMqI3IcjuQQEXkg8whOSkqK9NwccNLS0iyO8BB5GoYcIiIPVTXoLFmyBHq9ngGHmhVeriIi8mDJyclQqVTQ6/VQqVQMONSsMOQQEXkwjUYjBRy9Xm/zhTeJXBlDDhGRh6rag1NWVoa0tDSLzchEnoo9OUREHshSk7GlZmQiT8aRHLKJs2fPQqFQ4PDhwwBMt5pQKBS84SWRkxgMBotNxsnJyUhLS4PBYHBSZUSOw5DjxlasWIGgoCBUVFRIrxUVFcHHxwcjR46U7ZudnQ2FQoFffvkFANClSxcsX7681mN/9dVXGDp0KEJCQhAUFIS+ffvixRdftLq24cOHIy8vDyEhIQ38qojIFlJTU2sdqUlOTuZCgNQsMOS4sYSEBBQVFWH//v3Sa9nZ2QgLC8O+fftQUlIivZ6VlYWIiAj06NGj3uNu3boVU6dOxR/+8Af8+OOPOHDgAF577TXo9Xqra1OpVAgLC6v3jt72Vl5ebtV+BoMBRqPRLjUIIWRBlIiIHIMhx4317NkTERERyMrKkl7LyspCYmIioqOjsXv3btnrCQkJVh13w4YNGDFiBF5++WX07NkTPXr0QFJSEt555x2ra6t+uWr16tVo2bIlvvvuO/Tu3RstWrTAvffei7y8PNn7Vq1ahd69e8PPzw+9evXC+++/L9v+yiuvoEePHggICEDXrl2RnJwsCzKpqakYOHAgPvnkE3Tt2hW+vr6wdA9acz0bNmxAnz594Ovri19//dWqGnbv3o2BAwfCz88PsbGxyMjIsHip7rvvvkNsbCx8fX2RnZ1t9feOiIhsg43HtRECqDIS4lABAYCVIyDx8fHIzMzE/PnzAQCZmZmYN28ejEYjMjMzMWbMGOj1euzZs8fqkBIWFobPP/8cP//8M/r169foL6O6kpISLF26FP/85z/h5eWF6dOn46WXXsJnn30GAPjwww+xaNEivPvuu4iJicGhQ4fwxBNPIDAwEDNmzAAABAUFYfXq1YiIiMDRo0fxxBNPICgoCPPmzZM+59SpU/jyyy/x1Vdf1bl0fUlJCdLT0/HRRx+hTZs2CA0NrbeGW7du4b777sPEiRPx+eef49dff8WcOXMsHn/evHlYunQpunbtipYtW9rs+0hERFYSzYhOpxMAhE6nq7GttLRUHD9+XJSWlppeKCoSwhR1HP8oKrL6a1q5cqUIDAwU5eXlorCwUHh7e4vLly+LtWvXiuHDhwshhNixY4cAIE6fPi29r3PnzmLZsmUWj1lUVCQmTpwoAIjOnTuLBx98UHz88cfi9u3btdaRm5srAIhDhw4JIYTIzMwUAMSNGzeEEEKsWrVKABCnTp2S3vPee++J9u3bS88jIyPF559/LjuuRqMRw4YNq/Vz33zzTXHnnXdKzxctWiR8fHxEQUFBre+pWs/hw4dlr9dXwwcffCDatGlT+XMihPjwww8tfu0ZGRl11lBdjZ9BIiISQtR9/q4LR3LcXEJCAoqLi7Fv3z7cuHEDPXr0QGhoKOLi4vDwww+juLgYWVlZ6NSpE7p27WrVMQMDA/Htt9/i9OnTyMzMxN69e/Hiiy/irbfewp49exAQENCoWgMCAhAdHS09Dw8PR0FBAQDgypUrOH/+PB577DE88cQT0j4VFRWy5uX//Oc/WL58OU6dOoWioiJUVFQgODhY9jmdO3dGu3bt6q1HpVKhf//+0nNrasjJyUH//v3h5+cnbR8yZIjF48fGxtZbAxER2Q9DTm0CAoCiIud9tpW6deuGjh07IjMzEzdu3EBcXBwA0yWnqKgo7Nq1C5mZmRg9enSDy4iOjkZ0dDQef/xxLFy4ED169MAXX3yBRx99tMHHAgAfHx/Zc4VCIfXLmJt+P/zwQwwdOlS2n/mS0969ezF16lQsXrwY48ePR0hICNauXYu//e1vsv0DAwOtqsff31/WGG1NDUKIGs3U5q+hOmvrICIi+2DIqY1CAbjJSSohIQFZWVm4ceMGXn75Zen1uLg4fPfdd9i7d2+jg4lZly5dEBAQgOLi4qaWa1H79u3RoUMHnDlzBtOmTbO4z65du9C5c2csXLhQes3cLOyoGnr16oXPPvsMZWVl8PX1BQDZ7DYiInIdDDkeICEhAbNmzUJ5ebk0kgOYQs7TTz+N27dvW5xZdfHiRWlGkFmnTp3w9ttvo6SkBBMnTkTnzp1x8+ZNvP322ygvL8fYsWPt9nWkpqbiueeeQ3BwMCZMmICysjLs378fN27cwAsvvIBu3brh3LlzWLt2LQYPHoxvv/0WWq3WoTX88Y9/xMKFC/Hkk09i/vz5OHfuHJYuXQoA9U6X37NnDx599FHs2LED7du3t2ndRERUE6eQe4CEhASUlpaiW7duspNnXFwcbt26hejoaERGRtZ439KlSxETEyN7rF+/HnFxcThz5gweeeQR9OrVCxMmTEB+fj62bNmCnj172u3rePzxx/HRRx9h9erVuOOOOxAXF4fVq1cjKioKAJCYmIi5c+fi2WefxcCBA7F7926bL0tfXw3BwcH45ptvcPjwYQwcOBALFy6Ulsiv2qdjSXFxMXJycqxeu4eIiJpGIWprKPBAhYWFCAkJgU6nq9Gsevv2beTm5iIqKqrekxVRVZ999hkeffRR6HQ6+Pv7N/o4/BkkIrKsrvN3XXi5iqiBPv30U3Tt2hUdOnTAkSNH8Morr+CBBx5oUsAhIiLbY8ghaqD8/HykpKQgPz8f4eHhuP/++/Haa685uyxqhNTUVCiVSouXPTUaDQwGA+/xROTG2JND1EDz5s3D2bNnpctLy5Yta/TaQeRcSqUSKSkp0Gg0stc1Gg1SUlLqXDGbiFwfR3KIqNkyj+CYm8eTk5OlgJOWlmbzxnYiciyGnGqaUR82uRj+7DlH1aCzZMkS6PV6BhwiD8HLVb8xD0vr9XonV0LNVclvN4StvjI02V9ycjJUKhX0ej1UKhUDDpGH4EjOb7y9vREQEIArV67Ax8cHXl7Mf+QYQgiUlJSgoKAALVu2ZB+IE2g0Ging6PV6aDQaBh0iD8CQ8xuFQoHw8HDk5uba9FYBRNZq2bIlwsLCnF1Gs1O9B8f8HACDDpGbY8ipQqVSoXv37rxkRQ7n4+PDERwnsNRkbKkZmYjcE0NONV5eXlxtlqiZMBgMFpuMzc8NBoMzyiIiG+FtHYiIiMilNfb8ze5aIiIi8kgMOUREROSRGHKIiIjIIzHkEBERkUdiyCEiIiKPxJBDREREHokhh4iIiDwSQw4RERF5JIYcIiIi8kgMOUREROSRGHKIiIjII7ltyElPT4dCocCcOXOcXQoRERG5ILcMOfv27cPKlSvRv39/Z5dCRERELsrtQk5RURGmTZuGDz/8EK1atXJ2OUREROSi3C7kzJo1C5MmTcKYMWPq3besrAyFhYWyBxERETUP3s4uoCHWrl2LgwcPYt++fVbtn56ejsWLF9u5KiIiInJFbjOSc/78eTz//PNYs2YN/Pz8rHrPq6++Cp1OJz3Onz9v5yqJiIjIVSiEEMLZRVgjIyMDarUaSqVSes1gMEChUMDLywtlZWWybZYUFhYiJCQEOp0OwcHB9i6ZiMilpKamQqlUIjk5ucY2jUYDg8GA1NRUxxdGVI/Gnr/dZiTnnnvuwdGjR3H48GHpERsbi2nTpuHw4cP1BhwiouZOqVQiJSUFGo1G9rpGo0FKSgp/j5LHcZuenKCgIPTr10/2WmBgINq0aVPjdSIiqsk8gpOSkiI9NwectLQ0iyM8RO7MbUIOERE1XdWgs2TJEuj1egYc8lhu05NjC+zJISIy8fX1hV6vh0qlQllZmbPLIaqTx/fkEBGRbWg0Ging6PX6Gj06RJ6CIYeIqBmp2oNTVlaGtLQ0i83IRJ6APTlERM2EpSZjS83IRJ6CIYeIqJkwGAwWm4zNzw0GgzPKIrIbNh4TERGRS2PjMREREVEVDDlERETkkRhyiIiIyCMx5BAREZFHYsghIiIij8SQQ0RERB6JIYeIiIg8EkMOEREReSSGHCIiIvJIDDlERETkkRhyiIiIyCMx5BAREZFHYsghIiIij8SQQ0RERB6JIYeIiIg8EkMOERERNVpJSQlOnTrl7DIsYsghIiJycampqdBoNBa3aTQapKamOrSe69ev49NPP4VarUbbtm0xbdo0h36+tbydXQARERHVTalUIiUlBQCQnJwsva7RaJCSkoK0tDS713DhwgVkZGQgIyMDWVlZMBgM0rYrV66guLgYgYGBdq+jIRhyiIiIXJw52FQNOlUDTtXgYy9z587Ff/7zH+n5HXfcAbVaDbVajQEDBkChUNi9hoZiyCEiInIDVYPOkiVLoNfrbR5wjEYj9u3bB61WC61Wi3Xr1qFXr14AgN/97nfIy8tDUlIS1Go1oqOjbfa59qIQQghnF+EohYWFCAkJgU6nQ3BwsLPLISIiajBfX1/o9XqoVCqUlZU1+Xjl5eXIysqSQs2lS5ekbenp6Zg/f36TP6OpGnv+5kgOERGRm9BoNFLA0ev10Gg0TRrJOXToEEaPHo2bN29KrwUFBWHixIlISkrCxIkTbVG20zDkEBERuYHqPTjm5wCsCjpXr17FN998A5VKJc2G6t27N8rLyxEaGorExESo1WqMHj0avr6+dv1aHIUhh4iIyMVZajK21Ixc3a+//oqMjAxotVpkZ2fDaDSiV69eUsjx8/PD/v370b17dyiVSgd9NY7DkENEROTiDAaDxSZj8/Oq07kB4J133sHq1atx8OBB2esDBw5EUlISKioq4O1tigDmxmJPxJBDRETk4upa7G/hwoX48ccfIYSQpnEfOnQIBw8ehJeXF0aMGIGkpCQkJSUhKirKQRW7BoYcIiIiN6PX67F9+3ZkZGRg3bp1yM/Px48//ojBgwcDAJ566incfffdmDJlCtq1a+fkap2HIYeIiMgNFBUVYdOmTdBqtfj2229RWFgobQsODsbp06elkDN06FAMHTrUWaW6DIYcIiIiF1X1EtT+/fvxwAMPSNvCwsKkGVEJCQlQqVTOKtNlMeQQERG5kNzcXGlG1KBBg7B8+XIAwIgRI3DXXXdh5MiRUKvVGDp0KLy8eJ/tunDFYyIiIicSQuCnn36Sgs2RI0ekbR07dsS5c+dc8r5QjsQVj4mIiNzQPffcg8zMTOm5l5cXRo0aBbVajaSkpGYfcJqCIYeIiMgBysrKsG3bNmzatAnLli2T1qnp378/9uwGHoiEAAAgAElEQVTZg3HjxkGtVmPy5Mlo27atk6v1DLxcRUREZCeFhYXYuHEjtFotNm7ciKKiIgDAtm3bMHr0aADAlStXEBAQgMDAQGeW6tJ4uYqIiMhF7N+/H8nJydi2bRvKy8ul1yMiIpCUlISwsDDptea8jo29MeQQERE10enTp2E0GtG9e3cApr6azZs3AwB69uwp9dcMHjyYM6IciCGHiIiogYQQOHz4MLRaLTIyMnD06FHMmDEDq1evBgDExMRg2bJlGD9+PHr37u3cYpsxhhwiIiIrCCGwc+dOKdj8+uuv0jalUomSkhLpuUKhwJw5c5xRJlXBkENERFQLg8EApVIJwBRcZs2ahWPHjgEA/P39MX78eGlGVOvWrZ1ZKlnAkENERFTFzZs3pRlRO3fuRG5uLgICAgAAM2fOxNGjR6FWqzFu3DjpdXJNDDlERNTsXbp0CevXr4dWq0VmZqZsRtT27dsxefJkAMBLL73krBKpEdymxTs9PR2DBw9GUFAQQkNDkZSUhJycHGeXRUREbu7TTz9Fhw4d8PTTT2PLli0oLy9H7969sWDBAuzbtw+TJk1ydomur7AQ2LgRSE0FXGj5PbcZydmxYwdmzZqFwYMHo6KiAgsXLsS4ceNw/PhxLqBERET1EkLg4MGD0Gq1GDZsmBRehg0bBgAYOnSoNNW7Z8+ezizV9d26BfzwA5CVBWRmAgcOAEajadu0acBvU+mdzW1Cjnm9AbNVq1YhNDQUBw4cwKhRo5xUFRERubKKigpkZ2dLM6LOnz8PAFCr1VLI6d69O/Lz89G+fXtnlurabt0Cdu0yhZqsLGD/fsBgkO8THQ3ExzuhuNq5TcipTqfTAUCd3exlZWUoKyuTnhcWFtq9LiKqX2pqKpRKJZKTk2ts02g0MBgMSE1Ndbljk/swGAx44oknsG7dOly/fl16PSAgABMmTMDUqVNl+zPgVFNUJA81+/bVDDVdu5pCTUICEBcHREY6odC6uWXIEULghRdewIgRI9CvX79a90tPT8fixYsdWBkRWUOpVCIlJQUAZGFEo9EgJSUFaWlpLnlscl03btzAvn37MG7cOACmn4MTJ07g+vXraNOmDaZMmQK1Wo0xY8bA39/fydW6oOLimqGmokK+T1SUPNR06uSEQhtIuKFnnnlGdO7cWZw/f77O/W7fvi10Op30OH/+vAAgdDqdgyolotqkpaUJACItLc3ic1c9NrmOCxcuiHfffVeMGTNGeHt7C29vb3H9+nVp+9atW0VWVpYoLy93YpUuqrhYiO+/F2LBAiGGDxfC21sIU8tw5aNLFyFmzhRi9Wohzp51ark6na5R52+3CznPPvus6Nixozhz5kyD39vYbxIR2Yc5fKhUKpuHEHsem5znzJkz4vXXXxdDhgwRAGSPfv36icOHDzu7RNdUXCzE1q1CLFwoxN13C+HjUzPUdOokxIwZQqxaJURurpMLlmvs+VshhAvN9aqDEAKzZ8+GVqtFVlaWdBO0hmjsrdqJyH58fX2h1+uhUqlkPXSufmxyDKPRiIqKCqhUKgDABx98gGeeeQaAaQXiYcOGSTOiunXr5sxSXUtpKbBnj2nmU1YW8N//AlXW/gFg6qFJSDBdgoqPN12OclGNPX+7TU/OrFmz8Pnnn2PdunUICgpCfn4+ACAkJITXV4nclEajkUKIXq+HRqOx2DDsascm+yovL8eOHTuQkZGBjIwMLFiwQAo2SUlJWLduHdRqNRITExEWFubkal1EaSmwd6881Oj18n06dqwZahQKJxTrQPYYVrIHVBuWND9WrVpl9TF4uYrIdbAnh6oqKioSX331lXj44YdFq1atZL/nJ02a5OzyXE9pqRCZmUIsWiTEqFFCqFQ1Lz916CDE9OlCfPSRWD57tkhbvNjiodLS0sSiRYscWX2DNZuenKZgyCFyDbWFDluEEXsem+yjtLRUBAcHy4JNu3btxGOPPSY2bNggSktLnV2i85WWCpGVJURqqhBxcUL4+tYMNRERQkybJsSHHwpx8qQQRqP0dnf//0Vjz99uc7mKiDyHwWBAWlpajctH5ueG6utxuMixqenOnTuHjIwM5OTk4L333gMA+Pn5YejQoTh58iTUajXUajWGDx8u3f27WSorM11yMk/p3rMHuH1bvk94eOWU7vh4oFu3Wi8/mX/+qy6vUHVZBU+9lOs2jce2wMZjIiLHEkLg+PHjyMjIgFarxYEDB6RtFy5cQIcOHQCY7vwdEhIChaf3iNSmrAz48cfKULN7d81QExYmDzXduze4p8YcbMy9au4ScBp7/mbIISIiu/j888+RmpqKkydPSq8pFAqMGDECSUlJmDlzZp2r1ns0vb5mqCktle/Tvn1lk3BCAtCjh00ahd1x1qHHz64iIiLXpdfrkZmZib59+6Jjx44ATKM4J0+ehEqlwtixY6FWq3HfffchNDTUydU6gV5vut+TefbTrl01Q01oqDzU9Oxp89lPzW3WYYNDzsyZM/GnP/2JN8UkImrmioqKsHnzZmi1Wnz77bfQ6XRIT0/H/PnzAQCTJ0/GF198gQkTJiAoKMjJ1TpYeXnNUFNSIt+nXbvKUBMfD/Tubdcp3dV7cMzPAXhs0GlwyLl16xbGjRuHyMhIPProo5gxY4Z0TZWIiDxbaWkp/vWvf0Gr1eL777+XXe5o3769rFk4JCQEDzzwgDPKdLzycuDAAXmoKS6W79O2rTzU9OnjsHVqLDUZW2pG9jQNDjlfffUVrl27hjVr1mD16tVYtGgRxowZg8ceewyJiYnw8fGxR51EROQkRUVFaNGihfT8ueeeQ/FvJ/Do6GhpRtRdd90FLy8vZ5XpWBUVplCTlWUKNj/8UDPUtGlTM9Q46fvTXGcdNrnx+NChQ/jkk0/w0UcfoUWLFpg+fTqeeeaZRt12wd7YeExEVD8hBH7++WdotVpotVro9XocO3ZM2v7qq68iICAAarUaffv2bR4zoioqgIMH5aGmqEi+T+vW8lDTt6/TQo2ncUrjcV5eHrZs2YItW7ZAqVRi4sSJOHbsGPr06YM333wTc+fObcrhiYjIQYxGI/bs2QOtVouMjAycPn1a2ubl5YXz588jMjISAJCenu6sMh2nogI4dEgeam7dku/TqhUQF1c5pbtfP4YaF9PgkFNeXo7169dj1apV2LJlC/r374+5c+di2rRpUmPZ2rVr8fTTTzPkEBG5ieeffx7vvvuu9NzX1xfjxo2TZkS1bdvWabWlpqZCqVRa7BnRaDQwGAxITU1t2ocYDJWhJisLyM4GCgvl+7RsKQ81d9zBUOPiGhxywsPDYTQa8dBDD+HHH3/EwIEDa+wzfvx4tGzZ0iYFEhGR7RQWFmLTpk3QarV45ZVXEBMTAwAYO3Ys/vnPf2Ly5MlQq9UYP368rA/HmZRKpcXm2KrNtA1mMACHD1eGmp07LYeaUaPkoaY5r8LshhoccpYtW4b7778ffn5+te7TqlUr5ObmNqkwIiKyjcuXL2P9+vXQarXYtm0b9L/dnTo6OloKORMmTEBBQQFUKpUzS7XIJrckMBiAI0fkoUank+8TEmIKNeZ1avr3Z6hxc1zxmIjIQ126dAkPPPAAdu/ejaq/6nv06AG1Wo2pU6daHI13VQ26JYHRCPz0U+WU7p07gZs35fsEB8tDzYABDDUuird1sAJDDhF5KiEEjhw5gkuXLmHixIkAgIqKCoSFheHatWuIjY2FWq1GUlISevfu7bYzomq9JYHRCBw9Kg81N27I3xwUVBlq4uOBmBiGGjfB2zqQW3JIQyGRhzIYDNi1a5c0I+rs2bOIjIzEr7/+CoVCAW9vb6xduxY9e/aUZka5s6q3JCjX6/H/Zs3CUz17moLNjh2WQ83IkfJQ483TXnPCf21yKrs0FBJ5uMzMTKxZswbr16/H1atXpdf9/PwwaNAg6HQ6afLHmDFjnFWmTWkWL8aXqanYNHEi7vX1RcnmzQh4/335Ti1ayEPNoEEMNc0c//XJqWzSUEjk4XQ6HQIDA+H92wk7IyMDn3zyCQDTRI/77rsPSUlJGDduHAIDA51Zqu0YjcDx40BWFk588AGePn4cyQCwcSMAIACAXqXCNr0eqrFjcY9GYwo1XHWfqmDIIaerGnSWLFlSf0MhUTOQl5cnzYjavn07Nm/ejNGjRwMApk6dCoPBALVajVGjRnnG7XSEkEKN9PhtlKq3eZ+AAGDECGlKt+rOO7H///4PBoMB9wwd6py6yaWx8ZhcRq0NhUTNxMmTJ6X+mr1798pmRKWmpmLRokVOrM7GhABOnJCHmitX5PsEBAB33125Tk1sLEdqmik2HpNbq9pQqNfrodFoOJJDzcrRo0fRv39/2WtDhgyRZkT16tXLSZXZiBDA//4nDzUFBfJ9/P1NocY8pTs2FnDBdXvIfTDkkNNV78ExPwfAoEMep6KiAj/88AO0Wi0CAgKk+0D169cP0dHRiIqKglqtxpQpU9CxY0cnV9sEQgA5OfJQc/myfB8/P3moGTyYoYZsiiGHnMpSk7GlZmQid1ZaWorvv/8eWq0W33zzDa5duwYAaN26NTQaDby9vaFQKHDixAn37a8RAjh5snKdmqwsID9fvo+fHzB8uDzU+Po6oVhqLhhyyKkMBoPFJmPzc4PB4IyyiGzm5Zdfxvvvv4+SkhLptdatW2PKlClISkqS7etWAUcI4NQpeajJy5Pv4+tbGWri44GhQxlqyKEYcsip6lrojyM45G4uXryI9evXY+bMmfD39wdgaqgvKSlBZGSk1F8zcuRIaTq42xACOH3aFGbMwebSJfk+vr7AsGHyUFPHfQ6J7M3N/l9GRORacnJypBlR//3vfwEAERERSExMBAA8+eSTUKvVGDRokHvdSkEI4MwZeai5eFG+j0olDzV33cVQQy6FIYeIqIEuXLiA999/H1qtFv/73/9k24YNGya7k3enTp3QqVMnR5fYcEIAubnyUHPhgnwflcoUZKqGmt9GrIhcEUMOEVE9ysvLcePGDYSGhgIASkpKpFlRPj4+GD16tDQjKjw83JmlNow51JiDzfnz8u0+PqZLTuZ1au66y7R2DZGbYMghIrKgpKQEW7ZskWZEJSQk4KuvvgIA9OjRA88//zyGDBmCiRMnSveJcnlnz8qndP/6q3y7jw8wZEhlqBk2jKGG3BpDDhE1K6mpqVAqlRYb2xcsWIATJ04AAL777juUlpZK2w4ePAiDwQClUgkAWL58uWMKbopff5WHmrNn5du9vU2hxjyle9gwwFPufUUEhhwip6nrZKvRaGAwGOqcfUaNo1QqLa7BpNFopEtQZl26dJFmRN19991SwHFZ58/Lp3Tn5sq3e3ub1qYxh5rhwxlqyKMx5BA5SV0nW/MCiWR7ycnJuHLlClJSUrBy5Ur89NNPePfdd5GSkoKxY8ciPz8farUaarUaAwYMcO0ZURcuyEPNmTPy7UplZaiJjzetLtyihePrJHIShhwiJ7G0srOlFaCp6YxGI/bt2ydN9c7JyQFgmiUVGhqKiooKpKWlYeHChfDy8nJytXW4eFE+++n0afl2pdJ0v6eqoSYoyPF1ErkI3oWcyMnMwcZ8c1IGHNvKzs7G1KlTcanKwnUqlQpjxozBli1bUFFRAZXKRe98f+mSPNScOiXf7uUlDzUjRjDUkEdq7PmbIYfIBfj6+kp3YXfJk62bKCoqwubNm9GyZUuMGTMGAHDp0iV06NABQUFBmDhxItRqNSZMmIC33nrL9cLlpUvAjh2VwebkSfl2Ly/gzjvloYa/y6gZaOz5m5eriJxMo9FIAUev10Oj0Tj/ZOtGrl69im+++QZarRbff/89bt++jbFjx0ohJyIiAtnZ2Rg8eDB8f7tvUvXLgubngINvJ5KXJw81v/wi3+7lBcTEVE7pHjECCAlxXH1E7k40IzqdTgAQOp3O2aUQCSGESEtLEwBEWlqaxedUu/fff1/ExcUJLy8vAUB6dO3aVcyfP18YjUaL76vte+yQ731enhBr1wrx5z8L0auXEKZ1hisfCoUQgwYJ8eKLQnzzjRA3btivFiI30tjzN0dyiJzEUpOxpWZkAoQQOH36NLp16ya9ptVqsWPHDgBATEwMkpKSoFar0a9fvzpnRBkMDrzz/eXLlSM1WVnAb2vwSBQKYODAypGakSMBd1lYkOg3paWmQcm8PCA/H7jvPtMdQFwBe3KInITr5NTNaDRi79690Gq10Gq1OH36NC5cuIAOHToAANavX48zZ84gKSkJXbp0cW6xZgUF8lBz/Lh8u0IBDBhQuU7NyJFAq1ZOKJSobkIAOl1leKn6yM+XP9fp5O89cwaIirJtPWw8tgJDDpFr0+v12L59O7RaLdatW4fLly9L23x9ffGf//wHkydPdmKF1Vy5Ig81x47V3Kd6qGnd2sFFElUyGEw/ttWDiqUgc/u29cf19wfCw02Pjz8Geva0bd1sPCYit/fvf/8b06dPl56HhIRg0qRJUKvVuPfee9HC2QvZXb0qDzU//1xzn/79K2c/jRoFtGnj2BqpWSorqxlcLAWZggJT0LFWy5aV4aXqIyxM/jw42DRQ6WoYcojI4QoKCrB+/XpotVqMHTsWc+bMAQBMnDgRnTp1kqZ6x8fHQ+XMi/vXrgE7d1auU3P0aM197rhDHmratnVwkeSphACKiuoebTH/9/Xr1h/XywsIDa0ZVKqHmLAw0wiNO2PIISKHyM3NlVYc3rVrF4xGIwDg+vXrUshp1aoVzp4967xbKVy/Lg81P/1Uc59+/SpDTVwcQw01mNFoys91jbiYHyUl1h9Xpap7tMX8aNfOdBuz5qCZfJlE5CxCCIwcORK7du2SvX7nnXdKN7+syqEB58YNU6gxr1Pz00+mP5+r6ttXPlITGuq4+sitlJebJtTV16ibnw9UVFh/3KCgukdczP/dqpVrXjJyJoYcIrIZg8GA3bt344cffsCrr74KwBRa2rdvD6VSiVGjRkGtViMxMRGdOnVyfIE3bgDZ2ZWh5siRmqGmd+/KKd1xcQw1hOLi+kdc8vNNLVsNmcrTrl3dIy7mbbxRfONxdhU1e5zK3TRlZWXYtm0btFot1q9fj4KCAgBATk4OevToAQA4c+YMQkJC0MbRTbg3b1aGmqws4NChmmehXr3koaZ9e8fWSE4hhCnz1teom5cH3Lpl/XG9vWuGFkshpn17wMfHfl+fp2k2s6vef/99/PWvf0VeXh769u2L5cuXY+TIkc4ui9yYUqm0uPhe1cX6qKa9e/di2bJl2LhxI4qKiqTXW7ZsicmTJ0s9NwDQtWtXxxSl09UMNVXqAGCa21o11ISFOaY2coiKCtMU6fqadfPzTTOSrBUQUPtoS9Ug06aNqbGXXINbhZwvvvgCc+bMwfvvv4+7774b/+///T9MmDABx48fd87QN3kES6sMW1qNuLnLz8+XLj0BphlSX375JQDT/aHMKw7HxcXBx1F/ohYWykPNwYM1Q02PHpXr1MTFmc5E5HZu37auUffKlZo/AnVp3br+Rt3wcKBFC/a7uCO3ulw1dOhQDBo0CB988IH0Wu/evZGUlIT09PR638/LVVQXc7BxqbtSO9mpU6ekGVF79uzBggULsGTJEgBAaWkplixZgsTERMTGxsLLEX++FhYCP/xQGWoOHKh5RuveXR5qIiLsXxc1ihCmf1JrVtW9edP643p5mS4H1be2S1gY8Ns9WxuNl7sdw+MvV+n1ehw4cADz58+XvT5u3Djs3r3bSVWRJ0lOTsaSJUukO4I3x4AjhMChQ4eQkZEBrVaLn6stdpebmyv9t7+/P1577TX7FnTrVs1QU30ls27dKmc/xccDv932gZzHaKx/VV3zttJS64/r61v3JSNzcGnXDlAq7ff1VcXL3a7NbULO1atXYTAYpKFys/bt2yM/P9/ie8rKylBW5aJrYWGhXWsk96bRaKSAo9frodFoml3QMRgMGDduHK5duwbA9As8Pj5emhHVsWNH+xZQVATs2lW5Ts3+/TVDTXS0PNTYuyaS6PXWrap7+XLDVtUNCbFufZeQENe7ZMTL3a7NbUKOWfU1NIQQta6rkZ6ejsWLFzuiLHJz1X8pmZ8D9d8J3B2Hq0tLS7F161ZotVocPnwYBw4cgEKhgLe3N6ZOnYpLly5BrVZj0qRJaG3Pey0VFQG7d1eGmn37ap4du3aVh5rISPvV00zVtapu1SDzW/a1ikJROUW6rkbdsDBTU687qxp0zKPBDDiuwW1CTtu2baFUKmuM2hQUFNQY3TF79dVX8cILL0jPCwsLEclfkFSNpb+6LP11Vht3Ga6+efMmvv32W2i1WmzevBnFxcXStgMHDiA2NhYA8O6779qviOJiU6gxr1Ozb1/NVdGiouQrCnfubL96PJgQ1q+qW+VHoV4+PtY16oaGNp9VdQFe7nZVbvMjqFKpcOedd+L777+HWq2WXv/++++RmJho8T2+vr7wbWpXGdmVK4yCGAwGi391mZ8b6hl3d4fh6vfeew9z5sxBRZVA0bFjR2nF4YEDB9rng0tK5KHmxx9rhprOneVTurt0sU8tHqKiwvpVdcvLrT9uixb1N+qGh5tmI7naJSNXwMvdLkq4kbVr1wofHx/x8ccfi+PHj4s5c+aIwMBAcfbsWaver9PpBACh0+nsXClZKy0tTQAQaWlpVr3uysw1q1Qqm9S+aNGiWo+RlpYmFi1aZHFbTk6OeOONN8S+ffuk17Zv3y4AiD59+oiFCxeKffv2CaPR2KT6LCouFmLrViH+8hchRowQwsdHCNOgQuUjMlKIRx4RYtUqIXJzbV+DmyouFuLUKSF++EGIf/9biLffFuLVV4WYOVOI8eOF6N9fiNBQIRSKmt/Suh5t2gjRr58QY8eavu2vvCLE8uVCfPGFEDt2CPHLL0LcuuXsr969Vf995Y6/v1xdY8/fbhVyhBDivffeE507dxYqlUoMGjRI7Nixw+r3MuS4Jk/6BWEOOCqVqsnHsjYAGo1GsW/fPrFgwQLRp08fAUAAELNnz5beU15eLnJycppcUw0lJUJs2yZEcrIQI0cKoVLVPMt27CjEww8L8cknQpw5I4Q9wpWLMhqFuH5diGPHTNnvn/8U4q9/FeKFF4R46CEh4uOF6NlTiODghgUXpVKIDh2EiI0VYvJkIZ54QoiUFCE++ECIjAwh/vtfIc6dE6KszNnfAc/nSX+oubLGnr/d5nKV2TPPPINnnnnG2WWQDbl60561l9RsPVxd32WwF198EbNnz0ZGRgYuXLggvc/b2xsJCQkYNmyY7DXzLRaa5PZtYM+eyinde/eaptxU1aFD5eWnhARTj42HXd8wGKxbVTcvr2Gr6vr7138TxvBw043Puaqua2jq5W6yL7daDLCpuBiga/P19ZVCQllDzgx2Vlt/TfXGYkuzs2wR1mpbpFAIgc6dO+P8+fMIDAzEhAkTkJSUhEmTJqFly5ZN+kzJ7dumIFM11FT/t4mIkIearl3dNtTcvi1f9r+2EFNQ0LBVdVu1qrtJ1xxkgoPd9ltHZFcevxggeTZnNO1ZO0JT34iKeVtjZ2fV5fr16+jSpQu8vLyg/23EpOrdvd944w20aNECY8aMgb+/f6M+Q6asDPjvfyundO/ZUzPUhIfLQ010tEufmYUwrSlozaq6N25Yf1wvL9MMovoadcPCAD8/+319RFQ7hhxyuqasUdMUDZn6XdcltdTUVJsOV1+8eFFacTgrK6vG+19++WUsW7YMAPDQQw816Ng1lJWZZjxVDTW3b8v3CQurDDXx8abbJrhAqDEagatX6x5xMQeZkhLrj6tS1b+2i3mKtKNW1SWiRrJHg5CrYuOx63F2015Dm55t2VhcVdWZTgsXLpSah82P+Ph4cfDgQbF48eKmfV/KyoTIzhYiLU2I0aOF8POr2dXavr0QU6cKsWKFEP/7n8MbhfV6U9Psf/9raqL94AMhFi0S4sknhbjvPlOzbYcOQnh7N6xZNzjY1OQbF2f68ubOFeLNN03NwFu3mpqDr11rVn3RRG6j2TQek2dxdtNeQ5qebXlJzWg0Yv/+/dBqtdBqtXjzzTcxZcoUAIBarUZmZiZatGiBLVu2yOqJiYmBQqGwfqRLrzctuGdep2b37po3C2rfXr6icM+edhmpKS62rlH36tWGHdfSqrqWbsQYGGjzL4mIXBwbj4lQf9NzbZfUGtJYXF5ejh07dkCr1WLdunW4ePGitO3RRx/FJ598Itu/UQsl6vWm+z2ZQ82uXTVDTWioPNT06tXoUCMEcP26davqFhVZf1xvb+saddu3N63AS0SejY3HRI1U3whNU2/7AABXrlxBjx49cPPmTem1Fi1aYNKkSUhKSsLEiRNrvKeulZ6lzysvrww1WVmmO3ZXb0Bp27aySTg+Hujdu95QU1FhmkFkzaq61WeQ1yUw0LobMbZuzSnSRNR0DDkewhVuj+COrGl6bugltWvXruGbb75BQUEB5s2bBwBo164dIiIi4OPjg8TERKjVatxzzz0Nv+1IeTlw4IA81FS/8VCbNvJQ06ePFGpKS4H8s/U36hYUmEZprNW6tXXruwQFNezLJSJqCoYcD+EuN4l0JdaO0FgzonLu3DlkZGQgIyMDO3fuhMFgQEBAAGbPni1N7d68eTMiIiKgbMiUnIqKmqGm2nUf0aYNyofH4fod8fi1awJOqfogv8ALeblA3m55gNHprP9opdJ0Oaiuy0ZhYaYHbxFHRK6IIcdDuMNNIl2NLZqeP/30U7z11ls4ePCg7PUBAwZArVZDr9dLIScyMrL+oioqgIMHgawsiMwsiOxseBXLQ02xb2scbROHPb7x2KJPwI6rfVH6jRfwTf2HB0xrtlhzyahtW06Rdkcc1SWqxJDjQVz99giuxqqelyqMRiP27t2Lfv36SY1veXl5OHjwIBQKBUaMGCHd1TsqKqrWY5eVVfa2XL5YgfJ9h6uiOD0AAB3VSURBVBG0PxPhv2ShW142AipuAQAUvz2uoxV2IA5ZiEcmEvBzWT+ISzUbVkJCrFvfJSTEJZa5ITvhqC5RJc6u8kCuensEd6TX65GZmSnNiMrPz8c///lPTJ8+HQBw9uxZbN26FVOmTIG/f2i906MvXzKg043DSEAm4pGFkchGCApln3kDLaVQswPxyA/tj7AIr3ovG9liwWPyDLaYDUjkShp7/mbI8TC13eeIrFdaWooNGzZAq9Viw4aNuHXLB0AYgHD4+XXFPfdMR3T0iBohpnr/LwB4wYABOIJ4ZCEBmRiFnTVCTbF3CE51jEN+z3gUD46HKrY/wjsqpVV1vTneSo3A3wXkSRhyrODpIYd/vVmvvBy4fLlyhOXiRSMKCryQlwfk5pZiy5ajAMIBtAegsvq4wYEGxLf+CWN9snDX7Uz0vboT/np5t68xOASKuFFQmNepGTCAzS9kFxzVJU/BdXKaOVus5eIJSkqsX1VXHu+r9rj4AxgiO27btrU06rY3olvJT4g8nYVWhzPhvXsncP6m7L0IDgZGjpSmdHsNHMhQQ3bnjJveErkahhwP4ezbI9iTEKa7Q9e3om5enulu09YrB3AZQB4UisuYOXMcOnVS1Qgy7dubbtoIwHRXyKNHK6d079hR89bVQUGyUIOBA3nNiRzKWTe9JXI1vFxFTmMwWL+qbkNG2gMCLE+N3rDhQ+ze/R8AeQDyoVBcR1zcSCQlJSEpKQmdO3eueTCjEfj5Z3mouX5dvk+LFvJQExPDUENOU9tlal6+JnfGy1XkMm7frn/ExbyqrtFo/XFbtbJufReVqgzbt2+DVqvF4sWLERERAQAICSnHwYM7MW7cOKjVL2Dy5Mlo27at/EOMRuDYMXmouXZNvk9goCnUmFcVHjSIoYZcRn2jutu2bZM9r4rr6JCn4W9msooQQGFh3SMu5sfNm/Ufz8zLyzSDqL71XcLCTIvY1aawsBCbNm3CokVabNy4Ebd+u241aNAgPP300wCARx55BDNmzEBg1dtRC1Ez1FS/DXZgIDBihDzU8K6QXHTORVmz/hPX0aHmgiGnmTMaTed0a5p1q9/Mui6+vvWPuISHA+3aNa0H99SpU3juueewbds26KvcKTI8PBxJSUkYPHiw9FqLFi1Moeb4cdMdus2h5soV+UEDAipDTXw8EBvLUGMBF51zT1wdnZoThhwPpddXBpS6Lh1dvmzqjbFWcLA8pOTm7kZwcAkefniM7PWWLYElS2z/1/yZM2dw9epVDBlimv3Upk0bfP/996ioqECPHj2gVquhVqsxePBgeHl5mULNiROVoSYry3KouftueahRWT9tvLniydJ9cXV0ai7YeOxmiorqH3HJy6vZRlIXhcI0olLfirrh4aY8UJW9mxyFEDhy5Ai0Wi20Wi2OHj2Ku+66C3v27JH2+eyzzzBo0CD07t3bFGpyckxhxhxsCgrkB/X3l4eawYMZapqAi865L66jQ+6CjcduTAhTKKlrxMW8rdoNqOvk4yMPLrWFmNDQxl+Nsddf87t378a///1vZGRk4OzZs9LrSqUSLVq0kH4xQwhMi401hZnFi03/e/my/GB+fjVDjQ1um82eFJPk5GRpNEClUjHguAmuo0PNAUOOHVVUyFfVrS3E5OebVuC1VmBg/TdhDA8HWrc2Nfbamy2GvsvKyuBbJXj87W9/w9dffw0A8Pf3x/jx46FWqzF50iS0vnYNWL268vJTXp78YL6+wPDhlVO6hwyxSaipjj0pJjxZuh+uo0PNhmhGdDqdACB0Op1Nj7typRCvvirEzJlC3HuvEAMGCBEaKoRCIYRpnMa6R5s2QvTtK8SYMUI8/LAQ8+YJsWyZEGvXCrFjhxC//CJEYaFNS7cplUolAAiVSmXV/jdv3hSfffaZ+MMf/iACAwNFTk6OtO2rr74SjzzyiNB+/bUoOXLE9E3+4x+FiIio+Y3z9RUiPl6IxYtN36jSUnt9iTWkpaUJACItLc3ic0/X3L9+d1TbvxH/7ciVNfb8zZEcG3jzTeDUKcvblErTirn1re/Svr1dBhus1tRLL9b+NZ+Xl4d169ZBq9UiMzMT5VWGsDZt2oQe3bsDp0/jd9ev43cGAzB7NnDxovwgKhUwbFjllO6hQ+ueX25HzbmBk7cScU+evDo6UXUMOTYwdappbRhLIaZtW/e4TVFTLr1YO/S9Y8cOxMfHy97bu1cvPJaQgN+1bo0u+/cDnToBFy7IP0ClAu66Sx5q/P2b9gXbUHPtSeHJ0rYc1eNlzTo6RB7DTiNLLslel6s8RWMuPVjax2g0ij//+c8CgJg8ebL0eklJiQgMCBBJAwaIb3//e3EzMVGIyMial598fIQYOVKI5GQhtm8XoqTEfl+0DZi/B+bLdRzup8bgZSSi2jX2/M2QQzINPWEvWrRIpKWlifLycrF9+3Yxe/ZsERkZKQAIAKJdu3ZC5OYKsWqVEI88IgwdO1oONSNGCPGXvwixdasQxcUO+Vpr+1osSUtLE4sWLbL4ekODIVFt+PNEZBlDjhUYcqzT0AbiuXPnitatW0vBphMgnlCpxPbOncWttm1rhhpvbyHuvluIhQuF+P57p4Wa6hr6lzT/8iZ74MggUU0MOVZgyKlffb9gr1+/LtauXSuMRqP02tzf/148DIjPfH1FQYsWlkPN8OFCLFggxJYtQhQVOfrLslpD/pJuzMgPkTUa+ocGkadjyLECQ07dajvB///27j44qvre4/h3s8kmykMGylOAIBCtAakVE6EUGCDDFacOCrZDmWIGFNtSQ4tWtFibzZrVQH3kXgdxjLRDtaOOoaDSINFiVKAP4k0oc4MPkEaiMQIWk7SWXXbzvX8wu8OSDewu2XP2nH2/ZnZ0T85uvvmyu+ezv/M756xevVo3bNigc+fO1czMTB0tos333696662q48ZFDzXTpp0+rn7nTtWurrjqMDs88E0aZuL1B/REyIkBIad3Zwec1tZWXbt2rY4aNUpHiegSEa0W0UPRTvDjdKpOnaq6Zo3qa6/FHWrOV8v5licD36RjZ3YotRPm5ADREXJiQMjpndvt1l/96len73zyif7fL3+pT4voh72FmilTVH/xC9UdO5JyhkIzP+z5Jh2fVAildkAfgd4RcmKQjiHnXN+yKyoqtLS0VH+5bJmWDRqke664QvWyy3qEmu6MDNVrrlG9+27V2lpVA/pXUVGhJSUlUcNGMkcH+CadGPp24VJ5RCyVa0N6IOTEIB1Dztkbm3/961/6x2ee0V+MGaNPiej70UZqMjJUi4tVV69W/eMfDQk1vdXtdDojdhslc+PJN+kLwwiYffHegNkIOTFIx5Cjqvro3Xfr90V0X3GxHna5eoSaoMOhJy69VE+tWqX66quqX35pdsmqquGRnFDQCd1P1gcq31YvHHOZ7IvROpiJkBODtAk57e16bMMG/fv06XokyiHdQRE92K+fti5apIGtW1VPnEi5DXzoAzRa0EFqYiTH/vg3hlkIOTGwbcj5/HPtfvFF/eL739ejQ4f2nFPjcKhOnqzrMzJ0vogOy8rq8RSpNBx99u8MfaCGgg4frKmHb/npg9E6mIGQEwPbhJyjR1Vfekn19ttVJ06MOlLzvyL64siRuu2WW/SzpqaYvoGlyobqzFGls+suKSlht1GKSaWAjORiJAdmIeTEwLIh5+hR1Zoa1ZUrNRgl1KiINojo/2RkaOXVV+vmxx/X9vb28MPjCS+p9CGWKqEL55ZquzqRHLwfYSZCTgwsE3KOHVPdskX1pz9VnTQpaqg5MmiQ6s9+pvqHP2j3sWO6ZcuWqH9XIt+yU2E4mtEBwHi9BdbQ+2727NlRl/N+RLIluv3OFJjviy9E3n5bpL7+9O3vf++xygEReVNE6kXk/aFD5YZbb5V169aJiIhDRG666aaoTx0MBqWyslLKy8sjlofuB4PBiOVer1f8fr+4XC7x+/3i9Xp7PNYI8dYN4MI5nU5xu90iIhHvvV27domISElJScT6vB+R8pIUulJSyozkfPGF6tatqqtWqX7zm6oOR8/Rmiuu0OcGDdKbRHSIiBYUFOjq1at19+7dGggEklIWw9EA+BxAKmIkJ5WdOBE5UrN//+koc4bDOTnyp0BAlm7eLNlz54oMGyZdTz0l3zx6VDwLF8qkSZPE4XAkrUSv1ytutzti9CT032jf7ADY05nv+wceeED8fn/UUVXACgg5yfDll5GhprGxR6j5fPBgeSMQkJc7O+UtETl68qRkZGTIhPx8mTlsmIiIrFixwrCSg8GgzJ49u8fy0Afbrl27JBgMisfjMawmAOYoLy8PBxyXy0XAgWVlmF1ALFpaWmT58uUybtw4ueiii6SgoEAqKirE7/ebXdppHR0i27eL3HWXSFGRyODBIjfeKPL44yINDacDTmGhyIoVsv3mm2WEiIz45z/l5s5OeSU7W6bOny+/+c1vpL29XWbOnGnKn+DxeKSkpETcbrd4vd4eP6+vrxen02lCZQCMFm1uHmBJSdp91qd27Nihy5Yt0507d+rhw4f15Zdf1mHDhuldd90V1/MkZU5OIKCam9tjTk3gssv0o7lz9b+nTdPt1dXh1d9//33Nzc3VJUuWaE1NjXZ1dfVdLX2A/fFAeuMzAKko7Q4hf+ihh3TcuHFxPSZpE4+vu07161/Xr0pL9U+33aY3l5SED8MWEf3ud78bsbrf7+/b39/HUulcOQCMw6kbkKrSbuJxR0eHDB48+Jzr+Hw+8fl84fudnZ1JqcX/4ovyX/PnyzvPPSd6xtybyy67TBYuXCjf+973ItbPyspKSh19hf3xQHri1A2wG0uGnMOHD8sTTzwhjz766DnXW7t2rdx///1Jr8c1cKD4fD5RVSkuLpYFCxbIwoULZcKECUk9IipZUuVcOQCMda4DC/gMgCUlZVwpRhUVFeFdOr3d3n333YjHfPrpp3rppZfq8uXLz/v8J0+e1I6OjvCttbU1aefJ2bdvnx45cqTPn9do7I8HAKQaS+6uWrlypSxevPic64wdOzb8/21tbTJnzhyZNm2aPP300+d9/uzsbMnOzr7QMmNSVFRkyO9JJs6VAwCwE1NDzpAhQ2TIkCExrfvpp5/KnDlzpKioSH77299KRoYljn63FPbHAwDsxKF61lnqUlBbW5vMmjVLxowZI7/73e8iztcyYsSImJ+ns7NTcnNzpaOjQwYOHNhn9Xk8HnE6nVFHObxeLyfRAwDgAiS6/bbEcEhdXZ0cOnRIdu3aJaNHj5a8vLzwLRWELmp39gmzQrt/OIkeAADGs8TRVcuWLZNly5aZXUavos1biTa/BQAAGMcSu6v6SrJ2V4WEgk3o0GsCTnzY7QcAiMbWu6usory8PBxwOIle/Njt17c8Hk+v1xzyer0ERgC2R8jpQ1zU7sKUl5dLZWVlRNBht1/iCI0A0l4SztmTspJ27SrlJHp9iWtn9R1elwDsIO0u0JmIZIUcLmrX90IBx+VymV2K5REaAVhdottvdlf1gXOdRK+yspKT6MWJ3X59i7liANKVJQ4hT3Vc1K7vnD0HJ3RfhF4miguuAkhXhBwTcch0JK6d1fcIjQDSGSHHRKGjX0QiNzhnbpjSCdfO6luERgDpjpBjIs6UHIndfn2L0Agg3XHG4xTAmZIBAOhdottvQk6KyM7ODk8O9fl8ZpcDAEDK4LIOFsYh0wAA9D1CjsnOnIPj8/l6XNYAAAAkhonHJuLoFwAAkoeQYyKOfgEAIHmYeAwAAFIaE48BAADOQMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBTOTxeMTr9Ub9mdfrFY/HY2xBAGAjhBzARE6nU9xud4+g4/V6xe12i9PpNKkyALC+TLMLANJZeXm5iIi43e7w/VDAqaysDP8cABA/h6qq2UUYpbOzU3Jzc6Wjo0MGDhxodjlAWCjYuFwu8fv9BBwAOEOi229CDpAisrOzxe/3i8vlEp/PZ3Y5AJAyEt1+MycHSAFerzcccPx+f6+TkQEAsSPkACY7cw6Oz+eTysrKqJORAQDxYeIxYKJok4yjTUYGAMSPkAOYKBgMRp1kHLofDAbNKAsAbIGJxwAAIKWlzcRjn88nV111lTgcDmlsbDS7HAAAkKIsF3LuueceGTlypNllAACAFGepkLNjxw6pq6uTRx55xOxSAABAirPMxOPPP/9cfvjDH8q2bdvk4osvjukxPp8v4qRqnZ2dySoPAACkGEuM5KiqLFu2TFasWCHFxcUxP27t2rWSm5sbvuXn5yexSgAAkEpMDTkej0ccDsc5b/v27ZMnnnhCOjs75d57743r+e+9917p6OgI31pbW5P0lwAAgFRj6iHkx48fl+PHj59znbFjx8rixYvl1VdfFYfDEV4eDAbF6XTKkiVLZPPmzTH9Pg4hBwDAemx9gc4jR45EzKdpa2uTefPmSU1NjUydOlVGjx4d0/MQcgAAsJ5Et9+WmHg8ZsyYiPv9+/cXEZGCgoKYAw4AAEgvlph4DAAAEC9LjOScbezYsWKBvWwAAMBEjOQAAABbIuQAAABbIuQAAABbIuQAAABbIuQAAABbIuQAAABbIuQAAABbIuQAAABbIuQAQAI8Ho94vd6oP/N6veLxeIwtCEAPhBwASIDT6RS3290j6Hi9XnG73eJ0Ok2qDECIJS/rAABmKy8vFxERt9sdvh8KOJWVleGfAzCPQ9PoIlCJXqodAHoTCjYul0v8fj8BB0iCRLffhBwAuEDZ2dni9/vF5XKJz+czuxzAdhLdfjMnBwAugNfrDQccv9/f62RkAMYj5ABAgs6cg+Pz+aSysjLqZGQA5mDiMQAkINok42iTkQGYh5ADAAkIBoNRJxmH7geDQTPKAnAGJh4DAICUxsRjAACAMxByAACALRFyAACALaXVxOPQ9KPOzk6TKwEAALEKbbfjnUacViGnq6tLRETy8/NNrgQAAMSrq6tLcnNzY14/rY6u6u7ulra2NhkwYIA4HI4+e97Ozk7Jz8+X1tZWjtpKIvpsHHptDPpsDPpsjGT2WVWlq6tLRo4cKRkZsc+0SauRnIyMDBk9enTSnn/gwIG8gQxAn41Dr41Bn41Bn42RrD7HM4ITwsRjAABgS4QcAABgS06Px+Mxuwg7cDqdMnv2bMnMTKs9gIajz8ah18agz8agz8ZItT6n1cRjAACQPthdBQAAbImQAwAAbImQAwAAbImQAwAAbImQE6Mnn3xSxo0bJzk5OVJUVCTvvPPOOdffsmWLTJw4UbKzs2XixImydetWgyq1tnj6XF1dLTNnzpRBgwbJoEGDZO7cufK3v/3NwGqtK97Xc8gLL7wgDodDFixYkOQK7SPeXn/55ZdSVlYmeXl5kpOTIxMmTJDa2lqDqrWuePu8fv16ufzyy+Wiiy6S/Px8ufPOO+XkyZMGVWtNb7/9tsyfP19GjhwpDodDtm3bdt7HvPXWW1JUVCQ5OTkyfvx4eeqppwyo9AyK83rhhRc0KytLq6urtampSVetWqX9+vXTjz/+OOr6e/fuVafTqVVVVXrw4EGtqqrSzMxM/ctf/mJw5dYSb59/8IMf6IYNG7ShoUEPHjyot9xyi+bm5uonn3xicOXWEm+fQ1paWnTUqFE6c+ZMvfHGGw2q1tri7bXP59Pi4mL9zne+o7t379aWlhZ95513tLGx0eDKrSXePj/33HOanZ2tv//97/Uf//iH7ty5U/Py8vSOO+4wuHJrqa2t1fvuu0+3bNmiIqJbt2495/rNzc168cUX66pVq7SpqUmrq6s1KytLa2pqDKpYlZATgylTpuiKFSsilhUWFuqaNWuirr9o0SK97rrrIpbNmzdPFy9enLQa7SDePp8tEAjogAEDdPPmzckozzYS6XMgENDp06frM888o0uXLiXkxCjeXm/cuFHHjx+vfr/fiPJsI94+l5WVaUlJScSyn//85zpjxoyk1Wg3sYSce+65RwsLCyOW/fjHP9ZvfetbySwtArurzsPv98t7770n1157bcTya6+9Vvbu3Rv1MX/+8597rD9v3rxe10difT7bV199JadOnZLBgwcno0RbSLTPlZWVMnToUFm+fHmyS7SNRHr9yiuvyLRp06SsrEyGDx8ukyZNkqqqKgkGg0aUbEmJ9HnGjBny3nvvhXdvNzc3S21trVx//fVJrzed9LYt3Ldvn5w6dcqQGlLjlIQp7Pjx4xIMBmX48OERy4cPHy7t7e1RH9Pe3h7X+kisz2dbs2aNjBo1SubOnZuMEm0hkT7v2bNHNm3aJI2NjUaUaBuJ9Lq5uVl27dolS5YskdraWvnoo4+krKxMAoGAuN1uI8q2nET6vHjxYjl27JjMmDFDVFUCgYD85Cc/kTVr1hhRctrobVsYCATk+PHjkpeXl/QaCDkxcjgcEfdVtceyC1kfpyXat4ceekief/55qa+vl5ycnGSVZxux9rmrq0tuvvlmqa6uliFDhhhVnq3E85ru7u6WYcOGydNPPy1Op1OKioqkra1NHn74YULOecTT5/r6ennwwQflySeflKlTp8qhQ4dk1apVkpeXJ+Xl5UaUmzai/btEW54shJzzGDJkiDidzh7fCI4ePdojoYaMGDEirvWRWJ9DHnnkEamqqpI33nhDrrzyymSWaXnx9vnw4cPS0tIi8+fPDy/r7u4WEZHMzEz54IMPpKCgILlFW1Qir+m8vDzJysoSp9MZXjZhwgRpb28Xv98vLpcrqTVbUSJ9Li8vl9LSUrnttttEROQb3/iG/Pvf/5Yf/ehHct9990lGBjM5+kJv28LMzEz52te+ZkgN/Eueh8vlkqKiInn99dcjlr/++uvy7W9/O+pjpk2b1mP9urq6XtdHYn0WEXn44YfF6/XKa6+9JsXFxcku0/Li7XNhYaEcOHBAGhsbw7cbbrhB5syZI42NjZKfn29U6ZaTyGt6+vTpcujQoXCQFBH58MMPJS8vj4DTi0T6/NVXX/UIMk6nU/T0wThJqzXd9LYtLC4ulqysLGOKMGyKs4WFDk/ctGmTNjU16R133KH9+vXTlpYWVVUtLS2NmMW/Z88edTqdum7dOj148KCuW7eOQ8hjEG+ff/3rX6vL5dKamhr97LPPwreuri6z/gRLiLfPZ+PoqtjF2+sjR45o//79deXKlfrBBx/o9u3bddiwYfrAAw+Y9SdYQrx9rqio0AEDBujzzz+vzc3NWldXpwUFBbpo0SKz/gRL6Orq0oaGBm1oaFAR0ccee0wbGhrCh+qvWbNGS0tLw+uHDiG/8847tampSTdt2sQh5Klqw4YNeskll6jL5dKrr75a33rrrfDPZs2apUuXLo1Y/6WXXtLLL79cs7KytLCwULds2WJwxdYUT58vueQSFZEet4qKCuMLt5h4X89nIuTEJ95e7927V6dOnarZ2dk6fvx4ffDBBzUQCBhctfXE0+dTp06px+PRgoICzcnJ0fz8fL399tv1xIkTJlRuHW+++WbUz9xQb5cuXaqzZs2KeEx9fb1OnjxZXS6Xjh07Vjdu3GhozQ5VxuYAAID9MCcHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHAADYEiEHgGUdO3ZMRowYIVVVVeFlf/3rX8XlckldXZ2JlQFIBVy7CoCl1dbWyoIFC2Tv3r1SWFgokydPluuvv17Wr19vdmkATEbIAWB5ZWVl8sYbb8g111wj+/fvl3fffVdycnLMLguAyQg5ACzvP//5j0yaNElaW1tl3759cuWVV5pdEoAUwJwcAJbX3NwsbW1t0t3dLR9//LHZ5QBIEYzkALA0v98vU6ZMkauuukoKCwvlsccekwMHDsjw4cPNLg2AyQg5ACzt7rvvlpqaGtm/f7/0799f5syZIwMGDJDt27ebXRoAk7G7CoBl1dfXy/r16+XZZ5+VgQMHSkZGhjz77LOye/du2bhxo9nlATAZIzkAAMCWGMkBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC2RMgBAAC29P/eIgd3JtLvCgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x12a416f10>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using PyPlot\n",
    "\n",
    "# Model specification: y|x ~ 𝒩(f(x), v(x))\n",
    "f(x) = 5*x - 2 \n",
    "v(x) = 10*exp.(2*x.^2) - 9.5 # input dependent noise variance\n",
    "x_test = [0.0, 1.0]\n",
    "plot(x_test, f(x_test), \"k--\") # plot f(x)\n",
    "\n",
    "# Generate N samples (x,y), where x ~ Unif[0,1]\n",
    "N = 50\n",
    "x = rand(N)\n",
    "y = f(x) + sqrt.(v(x)) .* randn(N)\n",
    "plot(x, y, \"kx\"); xlabel(\"x\"); ylabel(\"y\") # Plot samples\n",
    "\n",
    "# Add constant to input so we can estimate both the offset and the slope\n",
    "_x = [x ones(N)]\n",
    "_x_test = hcat(x_test, ones(2))\n",
    "\n",
    "# LS regression\n",
    "w_ls = pinv(_x) * y\n",
    "plot(x_test, _x_test*w_ls, \"b-\") # plot LS solution\n",
    "\n",
    "# Weighted LS regression\n",
    "W = diagm(1./v(x)) # weight matrix\n",
    "w_wls = inv(_x'*W*_x) * _x' * W * y\n",
    "plot(x_test, _x_test*w_wls, \"r-\") # plot WLS solution\n",
    "ylim([-5,8]); legend([\"f(x)\", \"D\", \"LS linear regr.\", \"WLS linear regr.\"],loc=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Too Few Training Samples\n",
    "\n",
    "-  If we have fewer training samples than input dimensions, $\\mathbf{X}^T\\mathbf{X}$ will not be invertible. (Why?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  As a general recipe, in case of (expected) problems, **go back to full Bayesian!** Do proper model specification, Bayesian inference etc. Let's do this next. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  **Model specification**. Let's try a Gaussian prior for $w$ (why is this reasonable?)\n",
    "\n",
    "$$\n",
    "p(w) = \\mathcal{N}(w|0,\\Sigma) = \\mathcal{N}(w|0,\\varepsilon I)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  **Learning**. Let's do Bayesian inference,\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\log p(w|D) &\\propto \\log p(D|w)p(w) \\\\\n",
    "   &\\stackrel{IID}{=} \\log \\sum_n p(y_n|x_n,w) + \\log p(w)\\\\\n",
    "   &= \\log \\sum_n \\mathcal{N}(y_n|\\,w^Tx_n,\\sigma^2) + \\log \\mathcal{N}(w|0,\\varepsilon I)\\\\\n",
    "   &\\propto \\frac{1}{2\\sigma^2}\\left( {y - \\mathbf{X}w } \\right)^T \\left( {y - \\mathbf{X}w } \\right)  + \\frac{1}{2 \\epsilon}w^T w\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  **Done!** The posterior $p(w|D)$ specifies all we know about $w$ after seeing the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Too Few Training Samples, cont'd: the MAP estimate\n",
    "\n",
    "- As discussed, for practical purposes, you often want a point estimate for $w$, rather than a posterior distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  For instance, let's take a **Maximum A Posteriori (MAP) estimate**. Set derivative \n",
    "$$\\nabla_{w} \\log p(w|D) = -\\frac{1}{\\sigma^2}\\mathbf{X}^T(y-\\mathbf{X}w) + \\frac{1}{\\varepsilon} w\n",
    "$$ \n",
    "to zero, yielding\n",
    "$$\n",
    "\\boxed{ \\hat{w}_{\\text{MAP}} = \\left( \\mathbf{X}^T\\mathbf{X} + \\frac{\\sigma^2}{\\varepsilon} I \\right)^{-1}\\mathbf{X}^T y }\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Note that, in contrast to $\\mathbf{X}^T\\mathbf{X}$, the matrix $\\left( \\mathbf{X}^T\\mathbf{X} + (\\sigma^2 / \\varepsilon) I \\right)$ is always invertible! (Why?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Note also that $\\hat{w}_{\\text{LS}}$ is retrieved by letting $\\varepsilon  \\rightarrow \\infty$. Does that make sense?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  Adaptive Linear Regression\n",
    "\n",
    "-  What if the data arrives one point at a time?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Two standard _adaptive_ linear regression approaches: RLS and LMS. Here we shortly recap the LMS approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  **Least Mean Squares** (LMS) is gradient-descent on a 'local-in-time' approximation of the square-error cost function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  Define the cost-of-current-sample as  \n",
    "$$\\begin{equation*}\n",
    "E_n(w) = \\frac{1}{2}(y_n - w^Tx_n)^2\n",
    "\\end{equation*}$$ \n",
    "and track the optimum by gradient descent (at each sample index $n$):\n",
    "$$\\begin{equation*}\n",
    "w_{n+1} = w_n - \\eta \\, \\left. \\frac{\\partial E_n}{\\partial w} \\right|_{w_n}\n",
    "\\end{equation*}$$\n",
    "which leads to the LMS update:\n",
    "$$\n",
    "\\boxed{ w_{n+1} = w_n + \\eta \\, (y_n - w_n^T x_n) x_n }\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "-  <span style=\"color:red\">(OPTIONAL)</span> This is not a probabilistic modelling derivation. Is there also a Bayesian treatment of LMS? Sure, e.g., have a look at [G. Deng et al., _A model-based approach for the development of LMS algorithms_', ISCAS-05 symposium, 2005](./files/Deng-2005-A-model-based-approach-for-the-development-of-LMS-algorithms.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "---\n",
    "The cell below loads the style file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<!--\n",
       "This HTML file contains custom styles and some javascript.\n",
       "Include it a Jupyter notebook for improved rendering.\n",
       "-->\n",
       "\n",
       "<!-- Fonts -->\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<!-- Custom style -->\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.container {\n",
       "    min-width: 960px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "    margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running {\n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256);\n",
       "    border-radius: 0px;\n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Alegreya Sans' sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 125%;\n",
       "    font-weight: 400;\n",
       "    width:800px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: 400;\n",
       "    font-size: 45pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    margin-top:16px;\n",
       "    font-size: 22pt;\n",
       "    font-weight: 600;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: rgb(102,102,0);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-size: 14pt;\n",
       "    text-align: center;\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Nixie One', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 16pt;\n",
       "    color: rgb(163,0,0);\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.8em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"PT Mono\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       ".boxed { /* draw a border around a piece of text */\n",
       "  border: 1px solid blue ;\n",
       "}\n",
       "\n",
       "h4#CODE-EXAMPLE,\n",
       "h4#END-OF-CODE-EXAMPLE {\n",
       "    margin: 10px 0;\n",
       "    padding: 10px;\n",
       "    background-color: #d0f9ca !important;\n",
       "    border-top: #849f81 1px solid;\n",
       "    border-bottom: #849f81 1px solid;\n",
       "}\n",
       "\n",
       ".emphasis {\n",
       "    color: red;\n",
       "}\n",
       "\n",
       ".exercise {\n",
       "    color: green;\n",
       "}\n",
       "\n",
       ".proof {\n",
       "    color: blue;\n",
       "}\n",
       "\n",
       "code {\n",
       "  padding: 2px 4px !important;\n",
       "  font-size: 90% !important;\n",
       "  color: #222 !important;\n",
       "  background-color: #efefef !important;\n",
       "  border-radius: 2px !important;\n",
       "}\n",
       "\n",
       "/* This removes the actual style cells from the notebooks, but no in print mode\n",
       "   as they will be removed through some other method */\n",
       "@media not print {\n",
       "  .cell:nth-last-child(-n+2) {\n",
       "    display: none;\n",
       "  }\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n",
       "<!-- MathJax styling -->\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"],\n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "open(\"../../styles/aipstyle.html\") do f\n",
    "    display(\"text/html\", readstring(f))\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "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": 1
}