{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Regression" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Preliminaries\n", "\n", "- Goal \n", " - Introduction to Bayesian (Linear) Regression\n", "- Materials \n", " - Mandatory\n", " - These lecture notes\n", " - Optional\n", " - Bishop pp. 152-158 \n", " - In this and forthcoming lectures, we will make use of some elementary matrix calculus. The most important formulas are summarized at the bottom of this notebook in an [OPTIONAL SLIDE on matrix calculus](#matrix-calculus). For derivations, see Bishop Appendix C. \n", " " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Regression - Illustration\n", "\n", "\n", "

\n", "\n", "Given a set of (noisy) data measurements, find the 'best' relation between an input variable $x \\in \\mathbb{R}^M$ 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}^M$ and $y_n \\in \\mathbb{R}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume that we are interested in (a model for) the responses $y_n$ for **given inputs** $x_n$?, I.o.w. we are interested in building a model for the conditional distribution $p(y|x)$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, since $p(x,y)=p(y|x)\\, p(x)$, building a model $p(y|x)$ is similar to density estimation with the assumption that $x$ is drawn from a uniform distribution." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Bayesian Linear Regression\n", "\n", "- Next, we discuss (1) model specification, (2) Inference and (3) a prediction application for a Bayesian linear regression problem. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "##### 1. Model Specification\n", "\n", "\n", "- In a traditional _regression_ model, we try to 'explain the data' by a purely deterministic function $f(x_n,w)$, plus a purely random term $\\epsilon_n$ for 'unexplained noise':\n", "$$\n", " y_n = f(x_n,w) + \\epsilon_n\n", "$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In a _linear regression_ model, i.e., linear w.r.t. the parameters $w$, we assume that \n", "$$f(x_n,w)= \\sum_{j=0}^{M-1} w_j \\phi_j(x_n) = w^T \\phi(x_n)$$\n", "where $\\phi_j(x)$ are called basis functions.\n", " - For notational simplicity, from now on we will assume $f(x_n,w) = w^T x_n$, with $x_n \\in \\mathbb{R}^M$." ] }, { "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, i.e.\n", "$$\n", "\\epsilon_n \\sim \\mathcal{N}(0,\\beta^{-1}) \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Hence, given a data set $D=\\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$, the likelihood for an ordinary linear regression model is \n", "\n", "$$\\begin{align*}\n", "p(y\\,|\\,\\mathbf{X},w,\\beta) &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\beta^{-1} \\mathbf{I}) \\\\\n", " &= \\prod_n \\mathcal{N}(y_n\\,|\\,w^T x_n,\\beta^{-1}) \\tag{B-3.10}\n", "\\end{align*}$$\n", "\n", "where $w = \\left(\\begin{matrix} w_1 \\\\ w_2 \\\\ \\vdots \\\\ w_{M} \\end{matrix} \\right)$, the $(N\\times M)$-dim matrix $\\mathbf{X} = \\left(\\begin{matrix}x_1^T \\\\ x_2^T \\\\ \\vdots \\\\ x_N^T \\end{matrix} \\right) = \\left(\\begin{matrix}x_{11},x_{12},\\dots,x_{1M}\\\\ x_{21},x_{22},\\dots,x_{2M} \\\\ \\vdots \\\\ x_{N1},x_{N2},\\dots,x_{NM} \\end{matrix} \\right) $ and $y = \\left(\\begin{matrix} y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_N \\end{matrix} \\right)$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For full Bayesian learning we should also choose a prior $p(w)$:\n", "\n", "$$\\begin{equation}\n", "p(w\\,|\\,\\alpha) = \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\tag{B-3.52}\n", "\\end{equation}$$\n", "\n", "- For simplicity, we will assume that $\\alpha$ and $\\beta$ are fixed and known. \n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### 2. Inference\n", "\n", "- We'll do Bayesian inference for the parameters $w$. \n", "\n", "$$\\begin{align*}\n", "p(w|D) &\\propto p(D|w)\\cdot p(w) \\\\\n", " &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\beta^{-1} \\mathbf{I}) \\cdot \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\\\\n", " &\\propto \\exp \\big( -\\frac{\\beta}{2} \\big( {y - \\mathbf{X}w } \\big)^T \\big( {y - \\mathbf{X}w } \\big) - \\frac{\\alpha}{2}w^T w \\big) \\qquad \\text{(B-3.55)} \\\\\n", " &= \\exp\\big( -\\frac{1}{2} w^T\\big(\\underbrace{\\beta \\mathbf{X}^T\\mathbf{X} + \\alpha \\mathbf{I}}_{\\Lambda_N}\\big)w + \\big(\\underbrace{\\beta \\mathbf{X}^Ty}_{\\eta_N}\\big)^T w - \\frac{\\beta}{2}y^T y \\big) \\\\\n", " &\\propto \\mathcal{N}_c\\left(w\\,|\\,\\eta_N,\\Lambda_N \\right)\n", "\\end{align*}$$\n", "\n", "with natural parameters (see the [natural parameterization of Gaussian](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#natural-parameterization)):\n", "\n", "$$\\begin{align*}\n", "\\eta_N &= \\beta\\mathbf{X}^Ty \\\\\n", "\\Lambda_N &= \\beta \\mathbf{X}^T\\mathbf{X} + \\alpha \\mathbf{I}\n", "\\end{align*}$$\n", "\n", "- Or equivalently (in the moment parameterization of the Gaussian):\n", "$$\\begin{align*}\n", "p(w|D) &= \\mathcal{N}\\left(w\\,|\\,m_N,S_N \\right) \\qquad &&\\text{(B-3.49)} \\\\\n", "m_N &= \\beta S_N \\mathbf{X}^T y \\qquad &&\\text{(B-3.53)}\\\\\n", "S_N &= \\left(\\alpha \\mathbf{I} + \\beta \\mathbf{X}^T \\mathbf{X}\\right)^{-1} \\qquad &&\\text{(B-3.54)}\n", "\\end{align*}$$\n", "\n", "- Note that B-3.53 and B-3.54 combine to\n", "$$\n", "m_N = \\left(\\frac{\\alpha}{\\beta}\\mathbf{I} + \\mathbf{X}^T \\mathbf{X} \\right)^{-1} \\mathbf{X}^T y\\,.\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ " - (Bishop Fig.3.7) Illustration of sequential Bayesian learning for a simple linear model of the form $y(x, w) =\n", "w_0 + w_1 x$. (Bishop Fig.3.7, detailed description at Bishop, pg.154.)\n", "\n", "

" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### 3. Application: Predictive Distribution\n", "\n", "- Assume we are interested in the distribution $p(y_\\bullet \\,|\\, x_\\bullet, D)$ for a new input $x_\\bullet$. This can be worked out as (exercise B-3.10)\n", "\n", "$$\\begin{align*}\n", "p(y_\\bullet \\,|\\, x_\\bullet, D) &= \\int p(y_\\bullet \\,|\\, x_\\bullet, w) p(w\\,|\\,D)\\,\\mathrm{d}w \\\\\n", "&= \\int \\mathcal{N}(y_\\bullet \\,|\\, w^T x_\\bullet, \\beta^{-1}) \\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w \\\\\n", "&= \\int \\mathcal{N}(y_\\bullet \\,|\\, z, \\beta^{-1}) \\underbrace{\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z}_{=\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w} \\quad \\text{(sub. }z=x_\\bullet^T w \\text{)} \\\\\n", "&= \\int \\underbrace{\\underbrace{\\mathcal{N}(z \\,|\\, y_\\bullet, \\beta^{-1})}_{\\text{switch }z \\text{ and }y_\\bullet} \\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)}_{\\text{Use Gaussian product formula SRG-6}}\\,\\mathrm{d}z \\\\\n", "&= \\int \\mathcal{N}\\left(y_\\bullet\\,|\\, m_N^T x_\\bullet, \\sigma_N^2(x_\\bullet) \\right) \\underbrace{\\mathcal{N}\\left(z\\,|\\, \\cdot, \\cdot\\right)}_{\\text{integrate this out}} \\mathrm{d}z\\\\\n", "&= \\mathcal{N}\\left(y_\\bullet\\,|\\, m_N^T x_\\bullet, \\sigma_N^2(x_\\bullet) \\right)\n", "\\end{align*}$$\n", "with\n", "$$\\begin{align*}\n", "\\sigma_N^2(x_\\bullet) = \\beta^{-1} + x^T_\\bullet S_N x_\\bullet \\tag{B-3.59}\n", "\\end{align*}$$\n", "\n", "- So, the uncertainty $\\sigma_N^2(x_\\bullet)$ about the output $y_\\bullet$ contains both uncertainty about the process ($\\beta^{-1}$) and about the model parameters ($x^T_\\bullet S_N x_\\bullet$). \n", "\n", "- (See the [OPTIONAL SLIDE](#change-of-variable-derivation) below for the step in this derivation where $\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w$ gets replaced $\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z$.)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example Predictive Distribution\n", "\n", "- As an example, let's do Bayesian Linear Regression for a synthetic sinusoidal data set and a model with 9 Gaussian basis functions \n", "$$\\begin{align*}\n", "y_n &=\\sum_{m=1}^9 w_m \\phi_m(x_n) + \\epsilon_n \\\\\n", "\\phi_m(x_n) &= \\exp\\left( \\frac{x_n-\\mu_m}{\\sigma^2}\\right) \\\\\n", "\\epsilon_n &\\sim \\mathcal{N}(0,\\beta^{-1})\n", "\\end{align*}$$\n", "\n", "

\n", "\n", "- The predictive distributions for $y$ are shown in the following plots (Bishop, Fig.3.8)\n", "

\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- And some plots of draws of posteriors for the functions $w^T \\phi(x)$ (Bishop, Fig.3.9) \n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Maximum Likelihood Estimation for Linear Regression Model\n", "\n", "- Recall the posterior mean for the weight vector\n", "$$\n", "m_N = \\left(\\frac{\\alpha}{\\beta}\\mathbf{I} + \\mathbf{X}^T \\mathbf{X} \\right)^{-1} \\mathbf{X}^T y\n", "$$\n", "where $\\alpha$ is the prior precision for the weights." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The Maximum Likelihood solution for $w$ is obtained by letting $\\alpha \\rightarrow 0$, which leads to \n", "\n", "$$\\begin{equation*}\n", "\\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 if we have fewer training samples than input dimensions, i.e., if $Ngradient \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 so-called **normal equations** \n", "$\\mathbf{X}^T\\mathbf{X} \\hat w_{\\text{LS}} = \\mathbf{X}^T y$ and consequently\n", "\n", "$$\n", "\\hat w_{\\text{LS}} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T y\n", "$$\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 the (probabilistic) maximum likelihood solution ($\\hat w_{\\text{ML}}$) if the probabilistic model includes the following assumptions:\n", " 1. The observations are independently and identically distributed (**IID**) (this determines how errors are combined), and\n", " 1. The noise signal $\\epsilon_n \\sim \\mathcal{N}(0,\\,\\beta^{-1})$ is **Gaussian** distributed (determines the error metric) \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If you use the Least-Squares method, you cannot see (nor modify) these assumptions. The probabilistic method forces you to state all assumptions explicitly! " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Not Identically Distributed Data\n", "\n", "- Let's do an example regarding changing our assumptions. What if we assume that the variance of the measurement error varies with the sampling index, $\\epsilon_n \\sim \\mathcal{N}(0,\\beta_n^{-1})$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The likelihood is now (using $\\Lambda \\triangleq \\mathrm{diag}(\\beta_n)$ )\n", "$$\n", "p(y\\,|\\,\\mathbf{X},w,\\Lambda) = \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\Lambda^{-1} ) \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Combining this likelihood with the prior $p(w) = \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I})$ leads to a posterior\n", "\n", "$$\\begin{align*}\n", "p(w|D) &\\propto p(D|w)\\cdot p(w) \\\\\n", " &= \\mathcal{N}(y\\,|\\,\\mathbf{X} w,\\Lambda^{-1} \\mathbf{I}) \\cdot \\mathcal{N}(w\\,|\\,0,\\alpha^{-1}\\mathbf{I}) \\\\\n", " &\\propto \\exp \\left\\{ \\frac{1}{2} \\left( {y - \\mathbf{X}w } \\right)^T \\Lambda \\left( {y - \\mathbf{X}w } \\right) + \\frac{\\alpha}{2}w^T w \\right\\} \\\\\n", " &\\propto \\mathcal{N}\\left(w\\,|\\,m_N,S_N \\right) \n", "\\end{align*}$$\n", "with\n", "$$\\begin{align*}\n", "m_N &= S_N \\mathbf{X}^T \\Lambda y \\\\\n", "S_N &= \\left(\\alpha \\mathbf{I} + \\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \n", "\\end{align*}$$\n", "\n", "- And maximum likelihood solution \n", "$$\n", "\\hat{w}_{\\text{ML}} = \\left. m_N\\right|_{\\alpha \\rightarrow 0} = \\left(\\mathbf{X}^T \\Lambda \\mathbf{X}\\right)^{-1} \\mathbf{X}^T \\Lambda y\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This maximum likelihood solution 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: Least Squares vs Weighted Least Squares\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": null, "metadata": {}, "outputs": [], "source": [ "using Pkg; Pkg.activate(\"../.\"); Pkg.instantiate();\n", "using IJulia; try IJulia.clear_output(); catch _ end" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd3xTVf8H8HOzk2Y0adM26Uza0gFlFSp7DwFliYCKoKKCW3+ux719HtfzqKgoIO4BKiIo4gAZUvYo0paWNulOd/a8N/f+/kgJJXRRkty0+b5f/kFuM762aT6995zvORhFUQgAAAAIVwy6CwAAAADoBEEIAAAgrEEQAgAACGsQhAAAAMIaBCEAAICwBkEIAAAgrEEQAgAACGsQhAAAAMIaBCEAAICwBkEIAAAgrPXJIHzkkUdIkuzhnd1ud0CLAZciSRKW7gs+eKsHX88/iIAf+f2tjvXFDyyBQNDS0sLn83tyZ7PZLBKJAl0SaM/hcLBYLBaLRXch4QXe6sFntVoFAgGGYXQXEl78/lbvk2eEAAAAgL9AEAIAAAhrEIQAAADCGgQhAACAsAbTGQAAAISWo0ePnjp5QhIpHTdunFKpDPTLQRACAAAIFQaDYdl187mG+qtknEoCvf6U6abbVz/0+L8C+qIQhAAAAELFHTffuEBomZmT4rl5e47ivq82ZAzKmT1nTuBeFMYIAQAAhASj0VhdWjxTFe09wmJgjwxTfLzm7YC+LgQhAACAkKDT6RIlAp+DKZGCqurqgL5uuFwatdvtRUVFdFcRDEwmc8iQIbDUBQCgz4mJidGZ7T4H68wORVwcQsjtdpeWljY0NMTHx/t3ZZlwCcLvv//+oYceSklJobuQgCsqKjp69OjAgQPpLgQAAC6PTCaLVCYdrNWPjpd6D645Xb/s0RcOHTx4/x23JfOxaB6zsNmamD34vQ2fSKXSLp6t58IlCAmCmDt37saNG+kuJOCGDBkCiy8DAPqo9V9+s/ja2X/UVl8VxTW53FurzGNnzR0zdtw1E8esm6RWinieu/2qqbnpuvk7du/1y4uGSxACAAAIfQqFYt/RE7/++uuJwwclsuiPZs7MzMx8/qkn7syM8qYgQmiWOvrHPZrS0tIBAwZc+YtCEAIAAAghGIbNnj179uzZ3iMlhf9Mihb63G1gJK+kpMQvQQizRgEAAIQ0qSyq1ebyOdjiIiMjI/3y/BCEAAAAQtrCG5Z9Va5vf6TV7jreaL3qqqv88vwQhAAAAELatOnTk0ZPWfVXeX5Na7neurW0YcWu8rfWfsThcPzy/BCE9NNqtQ0NDd6bp0+frqqq6vCeJ0+e1Ol0waoLAABCxf/e//CpDz49FDv8Q724cciM3/KPTps+w19PDkFIswMHDmzZsmX+/Pmem62trQcOHEhKSurwzsOGDfviiy8IgghigQAAEBLGjB373/c//GbbL48++XRMTIwfnxmCkGbr1q2bM2fOhg0bPDffeeedZcuWdXH/MOmGBACAoKEtCF0ul15/0eCnxWLRn2c0GukqLMiKiooGDBjgWQiGIAi9Xt9+6SCr1YoQIknSZrN5jmRmZp48eZKWUgEAoF+iIQjLyspyc3OFQqFMJnM6nd7jy5YtS0xMTE1NTU1NHT16dPALCzK73f7pp586HI4ffvjBYrEghE6ePKlSqbx32Lp166lTp+65555vvvlmxYoVJSUlnuMEQXgCEgAAAkGj0Wzbtm3fvn1h8lFDQ0O9RCJ54403IiMjc3Nzfb60Zs2aW2+9NThlNDlQQQsVnNdKl6Bkoe8q2Hw+Pz09fdKkSddff73nSE1NTWxsrOffxcXFubm5iYmJ995771tvvZWYmOhdKDU6Orquri49PT04xQMAwofZbF61YlnjuaLh0XwjTt1Xb37m1dcWLV5Cd12BRUMQyuXyKVOmaDSaS7/kcDiqq6uVSiWTyQx0GfkN5HtFZKBfxWNBMuPu7A62gzh16tSQIUO8Nx0OB5vN9vw7KysLIdTS0iKXy3k83oQJE7x343A4drvvAu0AAHDlbl92wwR33bypqZ6bFhex8pnHU9SpI0aMoLewgAqtJdaee+65V1991WQyvfzyy/fdd19ndyNJcvfu3Vwu13MzOTm5F6dH85IZ85JpnitUUFCwatUq783Y2NiysjLPvz0Xjffs2TNq1CjPPb2Rqdfr4+Ligl4sAKCfa25urikpmnf1hUXLhBzWo0PjPvzfmxu++pbGwgIthILwgw8+UCqVCKH8/Pzp06fn5uaOGTOmw3sSBPH66697zxrHjRv32GOPdfa0VqsVw7D2g5Gho6ioqP1+SdnZ2Tt27PD8++GHH543b97Ro0ejo6Nra2ubm5u9dzMYDHK5vLPnpCjKZrN5Bh3p4nA4WCwWixVC765w4Hmr011FeLHZbCRJ9ptve1FRUWok3+dgZpTwv0fO0vuR4uOy3uo8Hq/bz6IQ+qjypCBCaMyYMVOnTt2zZ09nQcjhcHbu3Mnn+/7AOkRRlFAo9J4+hg6CIMRiMY93YT31uLg4b4/g888/f/LkyRdffPH06dN6vX7q1Kme4xaLRalUdvEmwDBMIBAIhb4L1AYT6zwaawhDnrc63VWEF8+vW78JwqSkpCa7b5tyvdUZGxsbUm8tv7/VQ7GPkKKo2tpaf+24GJrmzp27e/fupUuX+hyfP3/+7t27EULR0dHTp0/ncDgjRowYNGiQ9w5fffXV3XffHdRaAQDhISUlxYhxzrVeNFP04+KmxStuo6uk4KDhb3aCIDZu3Oi51vfxxx9HRESsWLHCarU+/PDDU6dO5XK5mzdvrqur886l7Jduuukmq9V68803+xyfNGnShg0bhg4dKpPJLn1UcXGxQqHobN0ZAAC4Qhu+3nTTgrnXKPkjYiKMTvxbjUl91fhF1y+mu67AoiEISZL0TBl9/PHHq6qqPFc4ORyOUqncvHkzjuODBg06efJkdHR08GsLmiVLOp2OfPvtt3d2OT4xMdEzmxQAAAIhOzs7/9Q/X37x+R+HD8pUMc8/uSgvL4/uogKOhiDkcDj/+c9/fA6y2exnn302+MX0nFar3b17d3Nzs1gsnjBhQvtJLn7X2eXvkLpMDwDol/h8/h13rkJ3rur+rv0FTGfonl6vv3vVnVt+3DokITqOz251uv/vwabxY8du+PQzuEoJAAB9HQRhN2w229SJE0TW5n3Lx8aL2mZ4ttpdLx0oHzvqqiPHTygUil48bVVV1ZEjRxISEkaNGlVXV1dYWDh9+nS/Fg4AAKBHQnHWaEh5843XydaGT2bneFMQISTjc/47LStHzHzs4Yd797R//fWXTCbbunUrQujHH3/sbANCAAAAgQZB2I2N69fdNyyBw/T9RmEIPTwy5fvz62VfrmuvvXbr1q1z5sxBCO3du3fSpElXXioAAIBegCDsislkqqzV5SoiO/xqRpSQw2KUlpb24pllMtnBgwfHjh1LUVR5eXlqauqVVQoAAKCXIAi74nA4EEI8VqffJS6L1bv1rymKEolEDAbj0KFD2dnZvS8RAADAlYEg7Ep0dHQEn6c12Dr8qtGJt1psycnJvXhmDMNWrly5c+fOjz76yLt8GgAAgOCDIOwKg8GYdfXMTWd1HX51c1FdTnZmQkJCL5553bp1KpVqxowZLS0t/XsNHQAACHHQPtGNF15+dWRu7pBo4fXZ8e2P769qeeuIdstP23r3tCRJGo3GNWvWvPrqqxEREf6oFAAAQG9AEHYjOzv7+y1bbliyeJumdVaKNFHMb7I5/6wy7NI0vvveezNmzOjd065evdrlcs2aNcu/1QIAALhccGm0e7NmzTpbem7CjSt/MvGeOqr7QocGzFp0urBw5e23X8nTcjgcf1UIAACg1+CMsEfi4uJeePHFF158ke5CAAAA+BmcEQIAAAhrcEZIj+LiYpVK1X57eoSQy+U6ceJES0uLUqnMycmBHd4BACAI4KOWHkOGDDl06NDw4cO9R8rLy6dPnx4XF5eQkFBeXs7lcvPz82msEAAAwgQEYah4+eWXp02btm7dOs/NhoYGeusBAIAwAWOEoUKv14tEIu/N2NhYGosBAIDwEcZnhFVV6PDhIL3WwIGouwVFV69evXDhwsOHD0+aNGn69OkTJkzAMCw41QEAQDgL4yA8exZ9912QXosgug3Cq6++uqSkZMuWLfv373/nnXemT5/+ww8/QBYCAECghXEQzpiBersuTIAkJiY+8MADDzzwQHl5eWZm5qFDh0aPHk13UQAA0M/BGGGocLvd3n8nJSVxuVyXy0VjPQAAECbC+IyQbh9++GFcXJzn33Pnzn3ttdeioqLy8vLYbPZXX32VkpKSl5dHb4UAABAOIAjp8dJLL7U/BWQymf/+97+3b99+7NgxhNC11167fPlyPp9PX4EAABAuIAjp8fjjj1968KGHHgp+JQAAEOZgjBAAAEBYgyAEAAAQ1iAI6XHq1CmLxeL5t81mO378uN1u99w0m82nTp1CCOl0uvLy8ksfW1paunXr1l9++aWysjJoBQMAQH8FQUiP22+//bvz7fzbt28fMWLEr7/+6rn57bffrl69GiG0YcOGRx55xOeBDzzwwOTJk3/88cdPP/101KhRO3fuDGbZAADQ/8BkGXpMmjRp7969t956K0Joz549M2bM2Lt378KFCz03J0+e3OGjzp07t3bt2urqas9KpARBeM8jAQAA9A6cEdJj4sSJu3fv9vx77969zz777J49e7w3J06c2OGjTCYTg8HwrrvGYrHar9MNAACgF8L3jLBcX75Lsys4rzUqYdTg2MHtj0yYMKGurq6iokIgEDidzrFjx5pMpubmZoPB0NDQMHbs2A6fZ9iwYVOmTMnIyJg2bdr48eMXLlyYkJAQlP8DAADot8I3CKuN1cd1x4PzWpG8SJ8glEgkQ4cO3bNnT0RExIQJExBC48eP379/f0tLy4gRIzo7z2MwGL/88suePXv+/PPPL7/88sknn9y+fXtn11EBAAD0RPgG4aSUSZNSJtFZwKRJe/fu9QbhhAkT9u7d29LS0nWwYRg2efLkyZMnv/LKK7fffvtrr70GQQgAAFcCxghpM3HixL179+7Zs8czIugJwn379nU2QIgQcrlcFEV5byYkJJAkGYxaAQCg/wrfM0LaTZgwoaqqSqlUqtVqhNCAAQMaGhqam5vbDxCeOXPG20ERGRk5efLkVatWLVq0KCUlpays7J133vnkk0/oqR4AAPoLCELaSCSS999/Pyoqynvk7bffNplMQqHQc3PmzJntBwtFIlFeXt4777xz4MCBY8eOxcbGHjhwYPDgwb7PCwAA4HJAENJp1apV7W8uXry4/c28vLxLd2KaOnXq1KlTA14ZAACEDRgjBAAAENYgCAEAAIQ1uDQKAACgb7ARlMaMJCTm3yW1IAgBAACELgohnZXSmCmNCVVakJuiFsWjRL++BAQhAACAkGMlUIWZ0pioEiNlwQP7WhCEAAAAQgJJUfU2VGpCpUZKZ0Pt1w8JqDAKwpKSknXr1tFdRcC1tLTQXQIAAFwGgwuVmyiNiSo3IYc7SOHXXrgE4ejRow8fPnz8eJBW2abRddddl5KSQncVAADQFZykqq1IY0IaM1VnpSH82guXIMzMzPzwww/prgIAAMKa3ok0ZkpjospMyEnHyV+H6AlCp9N5+vRphNDIkSPbHz99+nRBQUFGRsalK6oAAADoi3ASVVupEgNVYqAMLrqr6QgNDfUbNmwQi8UzZsx44IEH2h9/9913r7766v379y9ZsuSZZ54JfmEAAAD8Re9EhxrJz0vJ1wrIz0vJw40hmoIIISxo03K8GhoaeDzeli1b1q9fn5+f7zlos9ni4+N/++23vLw8rVY7cOBArVYbGxvb4TMIBIKWlhY+n9+TlzObzZ3tcwsCxOFwsFgsFitcLryHCHirB5/VahUIBBiG0V1IqLARlNaMNCbqnIkyBSz2FsXbB8UJ/fiENHxUdRhv+/fvF4vFniuiKpUqJyfnt99+W758edCrAwAAcBkubXinu6LLFip/s9fW1iYkJHhvxsfH19TUdHZnt9u9YcMGNpvtuZmRkeHZ5L2zO7vdbj+WCrrldrsxDIM/k4MM3urB5/meh+db3YpTlVbM0/BuJYL6HbistzqDwej2BxQqQeh2uxmMCwOWLBaLIIjO7kyS5PHjx71X3ths9ujRozu7M47jOB7gZQnAxXAcpygq+Ffdwxy81YPP8z0PnyAkKarBjp0zY+dMqN5Ow8iaB0EQPX+rs9lsJpPZ9X1CJQjj4uKampq8NxsaGmbMmNHZndls9tq1a3s4RojjOI/H80OJ4HLAGGHwwVs9+NxuN4/H6/dBaMFRuYkqNVLlJspx/kysu3AJIC6X69+3eqh8VI0aNaqysrKioiIlJcVgMBw9enT9+vV0FwUAAGEqpBreA42GIDxz5syaNWtKSko0Gs2qVauGDh161113yeXyO+64Y8GCBStXrty8efPcuXMHDBgQ/NoAACCchWbDe6DREIQSiSQ3Nzc3N9dz07se2Ntvv/3NN98UFBQsX758xYoVwS8MAADCUOg3vAcabaOdVwL6CEMc9BHSAt7qwden+wj1TlRiJEsNqMqKCLIvBUF/6CMEAABACzvRduUzoA3vfQ4EIQAA9Gf9oOE90CAIAQCgH7LiVIUFaUxUqRGZcQi/rkAQAgBAP0HXDu99HQQhAAD0bR02vIOegyAEAIC+hyBRlZUKk4b3QIMgBACAPiM8G94DDYIQAABCmqfhXWNCZw1ks4PuavojCEIAAAhFfbfhvc+BIAQAgFDhbXgvMyGjC8IvSCAIAQCANq2trZ9t/LjgVAE3Njll6g1EXDY0vAcfo/u7AAAA8DcrTn36066xucNdv32xiCgfWbrjhwfmHlz/PN11hSM4IwQAgCBp3/BebXB8df/KH69WR/E5CCGEJDPUMTfs/LJuzDXK7FyaCw0zEIQAABBYHTa8V508MF4hPJ+CCCHEwLDVWbJPd34BQRhkEIQAAOB/7RveO1ztzG5qVfN8B6diIrj2Vl2wagRtIAgBAMBvet7wLo1XnzHgPgeLW6xi9aBAFgg6AEEIAABXxNvwXmKkmuw9nfMZlzlsD8E7WKsfHS/1HGm0Ot8rap334MqAVQo6BkEIAAC9oXeiU82oykFW27BeNLxjGHbNm9tefOYG+VntUBm3xkYeb3VOfvZzkVwZiGpBFyAIAQCgp1wk0pqoUmNbwzuOM9lshFAvO/9EcuV1H+5t1hbXVZZGxsTfmD6YyeZ0/zDgbxCEAADQlUDv8B6tyopWZfn3OcFlgSAEAIAOwA7v4QOCEAAA2nga3jVmqtSIqi2w1lm4gCAEAIQ72OE9zEEQAgDCUbcN7yB8QBACAMII7PAOLgVBCADo53rX8A7CBwQhAKB/gh3eQQ9BEAIA+g+fhne6ywF9AwQhAKBvC3TDO+j3IAgBAH0SNLwDf4EgBAD0GdDwDgIBghAAEOosOFVuQtDwDgIEghAAEIpICtVYqVIjNLyDgIMgBACEEGh4B8EHQQgAoBk0vAN6QRACAOgBDe8gREAQAgCCBxreQQiCIAQABBY0vIMQB0EIAAgIaHgHfQUEIQDAb6DhHfRFEIQAgCsFDe+gT4MgBAD0BjS8g34DghAAcBnaNbxTTjj5Cw2E066v1YpjE7gRYrpr6ZMgCAEA3bjQ8G4gmxx0VwPacZgN+968r+Wfv9OjRDUmO4pRTX5yvUSRTHddfQwEIQCgY56TvxIDpTHT0PDuMBsay86wuDy5OpvNEwT51fsEiqK2/981dyuc8xZke44c0RkevX/m0i9OwHfsskAQAgAuCIWGd4ok89c+Wb17c55CYnaTuxosw1Y8MWje7bQUE8pqzxxRU/p56SneI3mKyIWJ1qLfNw2Zeyt9dfU9IRSEv/32W2VlpeffAoFg2bJl9NYDQJgItYb3Q+ufTyn+5eN5mQwMQwjZCfed371WJo1Nm3AtvYWFmqbywlEyts/BPDn/2NkjCILwcoRQEK5du9ZgMGRkZCCEIiMj6S4HgH7ORiCtmdKYqFIjZcbpruY8iqJKf/3iowVZnhRECPFZzNfHJt32+b8hCH2wBREmwveg2UWwouHz8/KEUBAihFasWHHrrfCHDACBEvoN706LMUrAYTGw9gcVQp7DWEVXSSFLNWLS9nVP3DVE2f7b9WW5OfW6BTRW1ReFVhAePHjQbDZnZWVNmzYNw7DuHwAA6AFvw7vGjOxECMbfBRx+hMnh8jnodJOI4XsNEETIYtMXPXD91jVPDovNihZWm+xvnW4iB02OH5RHd2l9TAgFYWJiotPpLC8vf/fdd1Uq1a+//spidVweQRB33XWX96ujR4++6aabOntah8PBZsOvUFA5HA4Wi9XZjw8EiM9bnaRQrQ2dM2FaC6q3Y32n4R2LSEg/XGe4Snnh+t5XxQ3JE+YRxCXXAelGEAS9f7LnLLonZvC4N755y1BUIo5LSrvjqZS8qSH4jfIvp9PpcPT044XNZjOZzK7vE4q/HhaLJTs7+/XXX1+6dGmHd+ByuW+88Yb3dz49PX3ixImdPZvZbBaJRAEpFHQCgpAWnre63om0FqQ1ozIT5SL75GUVU2Pt9gdmLkrgTokX4W7yO63pGCGZ985OFpdPd2m+XC4Xm82Gy1cIoeqC/Ko9W3BLa2TWVYPmLA/oD2tBnGVQnLCHd2YymQwGo+v7hGIQIoQWLlw4bNiwZ555psOvCgSClpYWPr9H32gIwuCDIAwyT8P7mXpbtYvXPxre3bjrn63rGwv2Mdnc+LHXZkxdFJphg+M4XHCiSPKP55cLq46vSBOLuKyD9bZNldbZb26LTskM0Csuirf3PAh7IlQ+qkiSJEnS89FpMBgOHz58ww030F0UACHNp+Hd5WJwOHTX5CdMNmfo9feg6++huxDQvcIdX6Q1n3p1qtpzc3hc5LREy73PLVvy2TF6C+u5UAlCi8WSlZU1adIkDofzxx9/jBw58rrrrqO7KABCTruGd8roO6cEABqU//r5+9kx7Y9kRAml7kZzY60oJp6uqi5LqAShSCTauXNnQUGB2+2+++67R44cSXdFAISQRjtVagyVhncA2rObWuURsT4HYwRsm6EFgvDyYBiWk5OTk5NDdyEAhIrQbHhvT1+rOfLBk81lBQwWO/GqGSNXPgu7H4ShyITU0pb6IbGS9gfL9NZMRRJdJV2ububSAACCiaSoOiv1dz25sYR8o8D9nYY83hyiKVh96sDO+6Y+LK76c45qx/T4WY17Nt2SZ9M30V0XCLZBNz763FGdxXWhZ+Or4gbhgBE8UZ9Z4CZUzggBCGd9qOHda/8bd301Xa0Q8hBCTCa2JCM2gsn4Zt2zkx9fS3dpIKgScq6y3PbynI+eHq0QydjYwQYrM3X41Kc30l3XZYAgBIAefXqHd6fFGOG2e1LQ62p19Fs79tBUEaBT5swbUifM1Z092WwxjkvPkcT1mYuiHhCEAARV/9jh3Y27OJes1sFkYCTZZ/+XwJVh8yOSho2ju4pegiAEIOD63w7v/MjoJjthw90C9oU4PKEzRKkH0lgVAL0DQQhAoNC7w3tAYRg2/JYn7vrutbcnpEh5bIRQaYvl8cO6qW+sp7s0AC4bBCEA/uRyU1oz6jcN75Un9tWf/pvJEyaNmBKTNqj9lwbNXVkijFy47pkIROBukiGNm/zvrXJ1Nl2lAlq4cZe+ViOKVnCFku7vHapCdK3RrsFaoyEuDNcabdfwTrlp+pVyuVwc/62x5rSafn5kbhrRPCOOa3NTWyqtzKxxk59cj12yfjFutzLZHAYrHJfcDOe1RnGH7cCaR2sO/JweLa63OHBJ3JSnNsqS0oPw0v12rVEA+pzQb3i/Env+vequGNvctGTPzSVZ6IXDR059996wJff73JPNjwh6dYB+O59cspBTu/y6QZ7V0P9pNN3z4NWLPz3GE0v98vwURRX/9k1d/s+Ewx49ZPzQRXd5drRguZxMkwn5NQjhjBD4Xz8+IwzlHd79eEZIuomvF6XvWXjRdU6zk1iwt2HpF6f88hL9Q9ieEbZWnTvx5LXfzkxtf/CLQt3+7KV5yx+78ufHHbat980YwzEti+FG2px11YZmnXVkWk5cQ42wud64bEXUxx9e+at49cOPKgD8zoJTlRZUaqRKjX2m4f1KOK0mKd83U0VcFmG30lIPCDWNmqKRct9TkZGxwm3FR3r9nHyTXlqrlWsK5Zpict+2e5sr4+wuiskgoiLSZMIqleRjc4vq7W3NyQMWJuFRV1a/DwhCADrWpxverxBPKGm1OX0O6h04WwBLiQKEEOLwI0y472+E2UWwe7bYrCfzpLUaaa1GXl4k1xZFVZ4jWSx9vEofr25SZ29xmB+dlVMXLabO9+eIENrx09kF8SoWk4WQn4ciIAgBuEj/aHi/QhiDqRwx5avigpuy4jxHKIRePVabtfAhegsDISIhZ9SmWqOdiOOzLjSSflVmTFyxyOeeLJdT1FQbU14k1xZJa7XSGo20VitqqtMnqBvV2fp4ddXQsccX3tGYNsjRbt7pjt+/fEwh9UlaMY/jslk8I4X+BUEIAMJJqtqK+lPD+5Wb8Mh7m59c8vufZTMUPAeJtlaao8bMGztvJd11gZDAEQhz73hx0SfPPz08JidGXG91vnemuS4qc4kqU3pkl7RGE6MpkmuKPJlnliv18aomdXZd5vDCaYv08Sq9UoUwrIvnj4hWVBptyRKB9whBUg1WB1/i32uibWCyDPC/vjJZpl3DO0WQdFdzxfzbPuFRX3Kq7p9DLJ4gafj4SKXKv0/eD/hMliGc9rL838w1ZRFKVdromRyBP2c2hhq+Sc/J/xXf/L5Mp1UhRiaDFatv9mSePkGtj1c1qbIbU7ONimSS4bsaX7e0h/88987qjVNThRwWQoggqZePVNcMWTjqzhdQANonIAiB/4VyEPazhvf2AhGEoGvtg7CmIP+vl26ZHc/PEDHLLe5tVZZxj36gGj2T3gr9wjukJ287z9NEV5S42RxP2jWlZuvj1fp4dZM6C/ffdcuzv3979KOnc6IEPBbjRL1RNWPZqDtf8LSxQhAiBEEY8pUy4s4AACAASURBVEIwCEOh4T3QIAiDzxuELpv522VDvp+hjonger6kd+ALdpxbuPGIQCqntcbLw3I5ZDUaT9pJa7Xy8sKY8kKKwfBMY/GkXZN6YEtyulMQ8M9VN+6qLSs22p24XNbqbm6xVhjtVRZHZXZU6qFlj/rxhULoowoA/+rfDe/+1Vxx9ti6Z1q1RRy+MGXydcNveCg8V4rptXMHfp2fLPamIEJIymMvT488vvuHYdetprGwLjBxl7ixRlqrjSkvlGuLvdNYzHJlkzqrUT3QM6TXpM42RysCWglBIiNOmZyowVrXaK3Q2yrN9gqrs8rlqqKICuSuYjRIeNykCG6yhJ+cJB94R9YI/xYAQQj6lVBueA9ZZXu2Frz/fy+NVAy9OtmKE5+f+nbTn5uvX/83i8Pr/sEAIYSQrb5KJfBdfE4l4uyv09BSjw+Gm5DUV3knbcZoCuWaYu80Fn2CukmV5ZnGYlCkUJesoucvDgLpXUjvolrt+martsWqMTs0DofWhesYpI7ESxgYh89TCXlqGS8uVaxWiKbFitRRwnQO86Kzz6sUdv8WBkEI+oNwa3j3I4okD7z78C+z08VcFkJIwmXfN1TBLNCd2vzeiGWP0F1dnyGIS9Ke8J1wpTU7Bdnq4BfTvjO9rVdPU+wQRTapsz3TWArm3Nykzm5OHkAy/R8BDgKZCWTGkd5JNdn1rVatya4x2bUOpwZzaxlkHUlUMjCugKsSC9RJArU0KjtGME0WoY6KGMBh0TO9CIIQ9FXh3PDuR80VZzOlfE8Kei1Mle74ezsEYc+ljbl609p/3ZwZLRe0XR01OPDPzxkX/GthQF+32870ommL9PHq5pQMl1+XhCVIZCYovRPTOykTjhvsTS02ndmhtTo0yK1lElqK0BBENYshiuCpJAJ1glgVJ8oV864R8ZTRwgw2M7TWp+04CLds2TJr1qwezkYBIJig4d2/3LiLz/K9FMZnMwmXhZZ6+ihuhHjyM58veuXWOQmCLDGr3Oz+sdI87rG1EbIYf73EpZ3pck0R22lvTVB7prF02Jl+JdwkMhGU2YVZCErvRK1Ol95WY3LUmRz1blzDJrUMUuvGNThezWKIxHyVXKDOlqmjI3KlEdeLeIpIfgqbKej+ZUJAx0F49913OxyOJUuW3HPPPYMHDw5yTQD48Da8a8xUnRXO/PwpKil9V6OJpChGuwbn/BqDPPsqGqvqixKHjVv6ZUFZ/s5zVWdFCQOWjJnJ6e28SiaBixuqPWl3hZ3p3SIpZCUoC47pXZTZhcwEanW4Wmw1ZrvG5dJxsHqmW0O5tW5c48J1HKZUyFOkRqijJGqpIFcacb1UoJLwkxhY37642HH7RHFx8WeffbZ+/frW1tbc3Nw777xz2bJlAkGoZDu0T4Q4f7VP9LOG90DrdfvEgff/FffPz89dlcBjMRFCZ5pM9/xdO//D/aKYeH/X2N/4ZfcJUVOdd9KmtFYToymS1mgu7Uy/wmksnokqZpwy40jvQmYXMuOOVpvOZC9nkVou0rEpHUloXS6NC9cJuEohN07KV0oFaqlAJY1QSwWqSEEyhi67NT4QgtpH6HA4tm/fvm7duj///DMyMnLx4sX33ntvTk6OH1++dyAIQ9yVBKHLTdXYUImBOmvobw3vgdbrIKQo6uSmd//Z9I5CyDU5cEZU/MQn1kUlZ/i9wv7ncoOws870JlVWU+rA9r16vetMbz90ZyaQBUd6J2VyOQ32Ogap4ZIVTLeGSekot87p0jhcOiFXIYtQi7gKEU8ZgoHXGXoa6k+fPv3OO+9s3LgRw7CJEyc+9NBD1157LXZl5+NXAoIwxPUiCEOn4b3qxP7Cb98y1GoliuTsxQ+m5E2hs5rLceUN9dbWBp4wksnhdn9XgBDqMgh5ZoOsRuOZxtJFZ3pz8gDXZa7E5jN05z3Pa3U4CFzHR1oe0mCEliK0pLvOhdfbXXUinlImUEkj1EKuQsRVtAUePwXDAtUpEVA07FB/5syZDRs2bNmyhcVizZkzp66ubt68ebNmzdq2bVtILR0C+hxvw/s5EzK5QmLk7+Dap/CD3780LDY1W1lhMP7n/bur82ePf/C/dNcVJBGyWLpL6JM660z3brDgGdJrTB1oiYrr4XNeOnSndyK9izLjyIHbIxgNHKqcRVYw3RoC1zpdGqdL5yb0Yu+JnVQl5GZ7zvP6buAFTadJZrfbN2/evG7duvz8fKVSef/99995553x8fEIoR07dsydO/fPP/+8+uqrg1gq6A9CueG9tepcy97NW2ane6aNZEQJN04TLv715ybN7XJ1drcPB+HApzM9uuxMbMXZDjvTeziN5dKhOwtB6Z2YwUVxGHoxQydA9QxSS+EagtC6HBrcoUGkg8FTijyBJ1FJBblCrkLEU0TyVTReqOvTOg7Cxx9/fN26dUajccqUKd999928efPan/7Pnj07NTW1trY2WEWCPs9KoAozFeIN79rDfyxKjmg/eRJDaKlK+PuBXyAIwxPfpD8/mKftsDP95Oyb9Ok5zSkZXW+w0OHQnRnHDC6KyUASlp5LaTmUjk3VE7jG7dKynRqOrdxNOgmeknn+DE8qyJUK1EJunIinDNp3IEx0HITbtm275ZZbVq9enZHR8YD5u+++O2DAgEAWBvq8Ptfw7rabRWzfK0hCNpOwmWipB3iYGmuM9dWyhFQ/9uRd6tLO9OjKUjeL7d1goShrkT5e3aTKxHkX5s+3HyMkKGTGOxq6cyI3iURsTMJu4aMKzK3B3Dq+ux5zapgOjd5WZsSYIp6CxVUIItRScVvgSQVqHjsycP+/oL2Og7CgoKDrUfeZM/vD3iIgEPROdLYVq7RSWqu7bzW8Rw0YfiD/y3kX/+13oMkpmzySporCXWt12e4XV0gcLSki3u5WG1KkT3lm45UPZLJcTllNuXeDBWmNJqa8kOVy+HSmN6TnOCPEPo91U8iGXxi6M7owI061Dd0RSMTGhGxSwDTwKC2L1PBxLYfQSXCd0a5psZbaMbZUoBLxFCKeUipUSWMWSQVqWUQql+Wf/nfQax0HIWzmAi6LT8M7QSAGAwVs5d5AUedN3bSWt6O8aXZq2745v2ub/zYyl46bTW9h4Qm3W39+aNb7o+NyYlI9R/ZWtbz44Oylnx7Fevze6qIzvf0GC56ZnO0f2DZ0Z/QdujPjFI+FRGxKysH4DD2Gl0tQhZDUOtwaB64zm3X1lhImgyMVqKQCtYinEAkVUsFYqUAti0jjsnxjFYQImPYJeq+fNbxjDMa8d3d+8uZ9b/14UCUVVhqs4gEj5q95nxGAhYlBt4r//G5RUkROzIXwmJgUNbS6qvLk/pTciR0+pNvOdE3e1Pad6e2XhzbXUz5Dd1IOErGRiI0JmXoZUxtJahVsjZ3UGiwas1N3zlbBZHAj+SmyiFSpQC0TZ4t406QCOleOBr0Gv+Hg8ngb3kuMyOAM9WG/y8WXRM186Wu3y2morxwZlxxSHXVOi/HkN283n8nniCITx83PnLm0f08RNJ47OSzKdx+oEVLmPk1hSu7EDjZY0J4lOFxvZ3pR1qID8erG1GwHm2fG207mLJ7TuwrKjJOtDuSmkIiNSbmUlItFsHApsyk2osGNa1xsjcmu1es1TTZtmb2KwxJ5zvCkArVSPDwjZo535Wi/rCwDSAJvrihxWIxydRZfLAt+AbBDPeiRy2p4JwiCwWAw+ty10RDWWHbmt8fnr86QjIuXmJzEt+WGY27Z/Hd/ax/V/WyH+oMbXpxb+dO8AXEMJ8Ey2phGG8toL9E0iznRquZ6ksn0brDQpM5uVqqrE9JbWCK9i9I7kRlHZhxdPHRHidiYlIskbBeTrMXcOgalszg0eptWb9XobVrjxYHn6cbrduVoCMIrp8n/9cD/Hhgo5UnYjBMN5qjhUyc+sqbrv0HpWVkm1EAQBoeNoLRm1IuG91AIQs2RXY1n8ll8cdLIqTFpg2isxC82rcj9MFeUKr2wec2aU3WFQ24Ysfwx75F+EISeznTPBgu8ouPcgzuGsRHT4iQkAiJK6BTx/lthinhwbUX22BphrN6F9E5kxikL7h26Q1IOJmIjERtJOC42WUvgGoLQWZw6b+CZHXU8tlTEU7QPvF6vHA1BeIXqSwvyn1rw5fRUGZ+DEKIQWldQt1uYM/35z7t4FA0ry4Cw4m1415hQhYUi+96fSchuav35/64dyDQtiOXaCPfmn9/nDJo06V8f9d0LiTZDc4TLnCq9aFGS5Vny63dvbh+EfUsXe6a3KFU1Kdnnhk47hDiN5X/fND0uRSI422J54URz/fyXRakLRCQS2ZGUg+QSJGDhbKrWjWsNVo3ZqbNYdE1WTalNa3bUiXhKT9edVKBWiIdnKxb1iYU0w8rpL/7zwog4TwoihDCEVg1R/vDT306riXvJlN3AgSAECPWRhvce2v3Sbf+X5J6pSvLcvD4LPZWf/8+P6wYvXEVvYb3mtJoieb6neiIOy2Wz0lLPpXCHrfC3b0ylJ7jyxLQpC6OSfJuML+1Mjykvsomk1SnZtbGqs7Hqn0ctOzM/qyAmo8nFvDB0N+52quL4C3+sdVdoRanD8/6ziClDeuvHepvWbKzTOHX684EnFaiEXIWIp4TA61v0lWcHTvLth8mMEuprtXEDhgStDAjC8NXnGt57gnA5zJrTM+dntT/4yHDlkm0b+m4QiuXxWoPFZ8vAwiZzVHJILGqhKz7x5zNLFicJhskFTWUHN/+ybsjwGeNHTuMf3iX552CsoSXRbCA5vLpYVVlC9pH47JLs64rHqQpiMpkREd6hOykHidloNhsTcxwuV13blUyLVi/QmGfWWZz15xxHztXskra2rRydJBvb11eOBgghnjCy1eESsC8a52pxEEl+2lu4hyAIw47eSWnMSGOiyk3IQe9GDwHgMOmjBL5TDaU8ttNipKWezlAU1aQpMtRqJXFJMWmDum6MY3K4KZOue+3YH4+NSGBiGEKo1e564lBt3gv/C1a9HWO5nNKq0paHr/0rSRDVamRV1LOazXe5qdK6L0p2bdMzcF68sDBO9KxJtM8ti39ye2RkpJSLpBxsJJsaz3RYHNq2wDNqNTaN3qYxO3QOXN9+SyC5aA6sHN2Pqabf+MnuNc/kJXqPVBptOpw1WZkSzDIgCMNC+OzwzpdENVrsPgebbE6+hIY52Z1prS7787mbEpElQ8wpt+C7HKxpL3wVkzqwi4eMu+/1/LVPTd+6aVis2Ogiz5lcYx98T5k9Img1MwmcW1vNr9bKygqjNEXyOm2cThPdqtMJo0a4bVIkcMVKbOlxbomAkPDP1La+cVj7w8I8hJCRgd+JOWN1jT//djV3+tQ6m0Zv0+htGsLtaB94yshFsHJ0GMqZd9svB35+YF/5igESCZd9SGf+qMQw4z9bglwGzBrtz+hqeKd31ujvzy1b5C5ZnNG2LiVJUQ/sq2Bf92T2rGW01OPDjbu+vmnIB6PkA+Vtb8tyvfXWPTVLvzzF6W5fOpfN0lxxlhshkiWmYZes8uyXWaOeBTNRnU5aVhRZq5XWamJ1GnVVUVKjVhepqIpKrLE0lNpaz7FYjYkZovvXkM3ahC1PPTFCaWTgVQx7PeZsZDj/cRl3G5uj4zkVTJsLkbEUN5kUFJQ7R857wtN+3m9WjoZZo35R/vcvlbu/c5pbo7JHD7n+bm5310WhfQIhCMIuhULDO71B6LJZdj61WKbXTovjWd3op0qzYvKS0Xe9Qksxlyrd97Ng81PPX5XQ/uC7J2vLp/5fzpzlV/LMPQ9CN4Vs5/e6w1v0khqtUlsYX1msrNckN2qzdGddbK42MasqOVsfrzYkqKyJqfasbJu59cfVE+8fLklT8WqY9hMW45aGVt7gDJP+HyqSYFJYDMWNo7hJbr6zmbI0EPenq5PdAjHVdtlpwa9lUz48SEu7dOBAENIC2idAx0Jnh3facQTCuf/bUV9yKv/MYbZANHXYOElcEt1FXWCsOTdK4vsnwiAp56S20O+v5bPXnbtVL6vVxNRpFfWaEbqigbri1PpzJJNVH6dqjle3qLNbRy06mpS6MyXDxY9w4Hq9TatvG7rLNxfqqsr/IpYbnscYSSQ/luLGRnJv4EX/cEAX5VQ9pxBMi4v2vKjJSczdcuSd6QNziAvT3ymEWmwunhC2UwChCIKwD+t1w3s4iMsYGpcxlO4qOsCVRNc7fH9YjXaCm9rTvct9eBfMbLYxTG7Kszy0zeoSNtYOqS8a1lCU1qxVNWkTGzTRrbompbpRnW1JVOsnjTumvnNH2iADl/QGnt6222z50nxM1+yzcjRPkSQba1p/aOfYIRJ00QnQfk3phDe2vPnKym+KK/Jk7Don2lNnUU6/+fPSXW/GSrxjfR//o0scfXXPF8sGIJggCPuYftDwHuZSR8/84ePnbhkYJ2C3DfLhbvKzc4ap987v4lHenV29C2Z6loc2uiguhWcYa7Nay8c0a7J0xanVxcp6rbSlziJX6uNVTersxrzcswnX74qTVYhJvd2TeVq99bjZqDPsr2AyuN51xeTCLHX01M5Wjj7e8qQY+V4GZDMZ/MioRR/tqy8tKCovFEbFLR08is0T5L//rznbv5+TKOQz0R86hyshZ/pDNM9xBaAzIRSEdrv9jTfe+Oeff9LT0x9//HGJBPbouqBdwztlJ4L96jZDM4vD63YqB+iJCFnMsFWvzFv31N1Z0oyoCK3B9l6RPn3pI5FKlZtEJuLCzq7eBTP1LkS4PT3mpNJpyNEVDawrTmrUKBu0iRVFstq2DRaaFcmtKar8icM18ZHVAkervUJv0+qtRXrbL8amKo7+wkKaCvHwAe1Wju5h5VHqQcd1DSMUFy5vmp1Ei4vyDPvFDRjSvgN6zD3/MV53d9npQ26nfcigvGhVVgfP2K85zIb89x6vO76LSZIkm5tzw8M5c2+Dc+LQFEKTZW688caWlpb777//s88+a2lp2bVrV2f3DJPJMiSF6m1UqQmVGmlreC/8+bPjn7ws5zGdbtLJEY17bG38oLxuHxUKa42GMgeBamprin/eaK8odMdlkGNvccpSfBbMlHKRwqFPbdaqm7QpVYXxlcXSWk10RYmbzdHHq3SpmUWpyopYYY2MWyHCW1zVequm1aoxOap7sXJ0DzVri/94ZM67YxM8810brc779lUk3vpK5vQl/viu9EmdTZYhXI5vV4x8cAB/XpocQ8iKu188Uq1Lnzrh4XeDX2T/029njVZXV6elpdXU1MjlcofDERMTs3fv3mHDhnV45/4dhKHT8H7qu/eIHe+9OS5ZyGEhhCqNtjv3VE7490+x6TldPxCCEF281137nV2NOMVhXtjrTspBUi6Sks6UpvLUqmJZnUZaq5WXF8aUF1IMRmNiyllVXHGCWBsnqoni1fAszc7qzlaOFrITosWpvVg5uueaNEV/v3mPo6GKxcDcXOFV97ymGj0zcC8X+joLwlM/fJR+YO0Dw+LbH5y1rXjmh/nCqF4OBgOvfjtr9Pjx4+np6XK5HCHE4/Hy8vKOHDnSWRD2PyHY8E5RVMHX//1j7gAeq20oK1kieH2U8oV1T89+4yd6awsdnQ7d4RQDQyIWJmJTUi4mZKMkIZJyGVIuJUNEZFONtFYbU17o3UWWo68rTZIXZCiPJUbWDGTproqs5QxuxmvNjmIRz+hdOTpGMDAj4obOtkpwuVwBTUGEkFydveCDvyiSJN0Ekx2onS7K9v98bts6c1OdLDlz8PInul5qIDQ1nty9OsH37++J8eK6ouMDxs+hpSTQhVAJwoaGhqioKO/NqKgonU7X2Z1dLtfVV1/NZLZ9QI8fP/7RRx/t7M5WqzVkF6owuDCtBZ0zYRVWLNR2eLfpm+IiON4U9BgSK2k5XOxyubp+bD87I3RTyEwgiwuzuJHBhZlxZCEwA44MTkRQmJCFIrmkkImEbCqGg1QCSsRCMi7FYSCGm4hsqJZVaWV1FYKSk9GlBbH6Opu75UyytDApclcMpzqT3ZDrbsBcRgIJuaxIQYSQEyXkKqQC1Qh+SqRAJeEndbhyNIGTCPn+FLr9ufiXOzAvt/ulW+J0p/6TI1dkRhU1aV56Ym7y4key564MxGtdORzHO7yoRlEYeclvtJtEbpIM8o+pX7LZbBZLT+/M4/FYrG6SLlSCkM/nO51O702HwxER0ekYPovFeuyxx7jctp0bk5KShMJOT5Mpiuriq8HXYcM7g4VCbR85KkJox90+B3E3iTFY3XZtM84LWHUB4em607soswuZiY73upNyUTQfqThIymnb69X78PZ7pgt055z6UqOprDSOdyhBdAZrbox1m9KoVr6bYWfFyJPiYoYIuQoRV5Hr15Wj+/p+hNqjf0XXnnxvstpzMy9e+n2cZNbXr2fOWMoTS+mtDSHUXHH2wJv3WXQajKI40thRD/w3NmtEh5dGlaNn7/j9rcGxFzopSYraU2uaO3RMKPyMSDfhtBj5kqju7xqSBAKBfz/VQyUIExMTKysrKYrynL1VVlbecMMNnd2ZwWBMmTKlh2OEIaLPNbxzhRI7O6LGZE8QX/g+by9vSrpqOo1VXbluh+68O7sqpUjIYog4VCSb4XNNgW/SSyu10lqNoK4EVZ+06AtbnNUaGfo7jq+RooYk3KRyStgxEtFoiTCt6udflkYpxwhlyW5BgolXY3TcvKli6qfb++7HUOBU7/3hjtSLpotzmYwZSZKqUwcGTLiGrqo8GssLdz167dtj43OuykQIVRpt9758U/b9a9LHdXCpc+CsGzd/vybqn/pbB8ayGFizzfXEweqUWctp/6GbGqr3vX63paJQzOM0WZ1Z8+8YcfNjDGaoBAFdQuX/f9y4cQihnTt3zpo168SJE+Xl5bNmzaK7qCvV1xveJ/xr3bJnljybGzM+QeYgyO/Lmj/VOq5f/zzddXWvh0N3Ug6Si5GIzZByKQmbwejoCjrL5RLV1UaUH0HVR/GWc3qrptlV84/QWSZnV4gJh5SUSSNl/GShbIlYkiEVqK+6eOXo6oL8lMoDjycmofNNL0kS/sqMyIM7vsi94cEgfkv6BsJqFgp8P5TETFRv7/GFsIA5vOaRt8fG58S0neQlSwSfTE29/v3HOwxCBpO1aN3+Q5+8+vWv3yPcxRZJh9/2RtqEa/1Yj93UeuKL15vOHGJz+QkT5g2ef8ely8/6cFqMW++Z8toI+eihWQgh3E2+cfzbv2rKpz79sR8L64tCJQjZbPaaNWtuvvnmYcOGnTp16q233uqjfYT9qeE9flDe3I/2b1j37Cu/H2VxeYljrrnhhX+xuKFyIt5p150TEWRb152QhYk4SMFH2ZEMERtJ2Bin7bPCJ/QwhBCTwFl1Z/Caw3hjIdF8Vm/V1GPN1Xx7WRRyslC0QBQtVkgFapH0Gl5c7qiI1Gk9WDm6taZ8qNj342mgjP9b+Wl/fR/6k8isvCNHTw+Lu+h3/3CLMyOtm4nKQWCoPpeTm9H+SLSAw8HtbpeTyeFeen8Whzd61YujV70YiGIazv3z++PzHxwoGzNE4iAcm/Z9uGnbx9et3cPmd9UVWvDD2pVpotHxbReZ2UzGk3mJs7fvsbTUh/lc1lAJQoTQokWLJk+efPbs2dTU1Li4PvZTobfhPXDEMQlTnt54uY8q/3uH7vCvhM0sHzph0DW3XMn0QpJC1vPLQ186dCdiY0I26d3ZNUmIpByGlIN4be/rTidJOXC9s/YIXp3vai21mauaidp6rFkrdDAQIxrny1CkXBovTpkojhqUgsWOihkqV/Zy4iJfJK2/ZOnzZpuLI/fdlRsghAZds/yrTW/nxRg9WUgh9E1xfYsoKWQnjpIUheiYi7fnlVs3TkpMlbbF3mO5gphC3b7P/j1m9ctdPKql8OCYeLHPwbEKcX1JQdqYPvaR618hFIQIoaioqLFjx9JdRU+FQsN7qHG7nD8/MjfdUXuPShQhYv2174Ovv3177js7u1322md56EuG7qjOh+66+hjyrBxtbfnH3nDSbiyzm6sa3boqrontJlOMDBkllnBiZbzk7JjRuTHDGelTmRFt+YQ7bPv/+0DT8fcGykX/2PAGkjvpqY8VWcMv9xuSMmLipv+Z7h9MiLhtv2skRa0vMQy7sdMh8HDGEYjmrvn96ZdXMo6VKsW80mazYuT02U+HxNpsUWmDj9bVj1RemLNTa3aQAkng2kg647QYuXZjqvSiv6WWZMR8/se2roOQyeY6CN8ZcDaCZHV0RhtWQisI+4QgN7wTLsexT/+t+esHCncJouKGr3o5efiEQL9or5346q3p7Ib7RyZ7bg6OFY9VGp9/6db57+9C7Ybu9E7KTCAL3jZ0Z3BRTMblDd355J93qwSzvsja/I/NWmlyNtQwGjgEUrdSKhNTikUms2OlgnR+7PWclEkW1TBXlyvG/f7sTQuZVTfPb1sYrMJgu+XJ6+Z99LcoJr6LR12KIxCNevDt+e8+dHeWdFC0sNpkf79YHzN9eeyAwZf1POEjUqla8MFup8VoaqwdnpgW/JjpzOj73njkvukvj6TGJ8oQQv80mv4vv2b0k58FvxKX3Sri+k5V5bOYbtzZ4f29lOPm/vjL695hToSQ003+XWdaOnCk/6vsU0JlZZnLEvyVZehqeHe7nJtWjr5RiS3PiuEwGVVG+2P51fIFDw65/p6g1XBZvrlpyE9T4jwr0XhN+vGs9YUCPSVwtw3dUZ6hOykHSbmYiI0kbMTpZpi/Tfu9gQyWMpPhrNVe00zUcwksycpJbSEHNDglrGi+KEUiHiCOHm5PyGpMHWi5zPEPU0P1/gen/jg7vf3Bn841bI2dPvaef1/WU3lYWurP/LjOWFYgVKoGXHObXJ3diyfpll825gVdMDXWHHzn4aaSExhC4oS00Q/8V5KQFvz9CCnS/fnCtL8WZLHa/ZFY2GR+qlp4zf92dP3AH+6aMoVvvHNgjITH3eJE7gAAIABJREFUPttsfupwXeLSx3Pm3xH4qv2p364sE5o8O7xrTNQ5I+Wio+H99E8fz5OTtw9qOwtJkvC/mJE27as3sq+9lc270qUjr0SHQ3d6F+W22n1SECEULeCOjjTHK4W8jieqXPrkuNXZZHHWn98bSKu3lhvN5wzOOj7FjXeJVQZsYL09s84qYUZFRCTJJFMtyTmN6QP1k1R6pQrHMBwh0xX83zVXlAyT+046GB4r/uTcyd49oTAqbtTtz15BRSAkiGMSZr6yqf0RHMeDXwbGYGZec8szBze/OCqJzWQghJpszscP1uS9sLnbB163dnfB92uX7vzCYTZIE9OGPfutMntEUKoOaRCEvkJhh3cv3aFfH0q5aC9TNpMxSinRnT2ZNDQYg6mdDd0ZXBSX1cHQ3eFYRbXJntiu9ZCkqCarfYZSfunUbpLCjfZqs0NncerOB55Gb9Ma7VU8xFe6pYl2vlqPptfbB59rkCFxhGy8Q5Guj1fpk9VNk7MbUzLqu5sv3jscgUjv8h1KMThwToQ8EC8HwOXKW/nsMTZ3yo8fZkWLHG6y0uoe/8hHPYk0jMEcuvjeoYvvDUKRfQgEYZvQbHin3AT7kiEyNoaRbn/OTL1o6A6n6vK320/9jlNYa+oUZu48EZvR86E7YsW//vXB/RumpvJZTIQQhdD/TtQmTZprdFZ5loq2OHXewDM76vgsiQyTKvHIJCt3dL1jYGXr0NJ6pU1kVqiaVNlNqdn6NLV+ovq4KhMP4hmwMmv4rgaLwYFH8i5c9fqsVJ+y9OGg1RBSbIbm/Hcfrj99gEGSGE847NanM6cvpruosIZh2MgV/8q96eHW6jI2TzBBkUx3RX1bWI8Rtmt4p0whuf7f4U9eHVfy/YqBCu8RkqJm/FQ8/9OTPFFkFw+8FEEhM36h6857ntfqRN6hOwFuMb42azzfNDdRiBDaUWM96hTOe3cnN8J3ynWH3KTT5Kg9+ef7xflfKGJZBJ+oYtrIaA7OdYl4ShEnNpoUJVp5KiNjQKMru8owpFDDczhaE9T6eM9/qiZ1dkN6jrNnLxdQmgM7jvz33kcGR+XGSRqszo+KW+vkA69+ZVPIrluLAjZG6LSaNt2S90SOZJYqGiGkd+BP5FdR45eNvO1pv79Wn9PZ7hMgoPrtNkyX5UqCsG81vDutpk23jHxusHRqShRCyIq7nz1U1Tr42rH3vtbh/TsbujPjyEEgEbtteUwpF4lYSMRpWzzz/NAd+us/q+dZji/OiPE+4feljVv4Q6c8ub79q3gCr/2Jnec8z+yoE/GUUoFKgEkjWp3xRtYgFD2kwZlTWhNdUyFqqjPLlU3qrEb1QH28Sp+g1ser9PHqQH3vrpipsebUl2+0lp6KkCvUM29O7WgBkZASoCA8vPHl8aVblrf7a8xNUVN/LF7y9Zmu27fDAQQhLSAIEepVEPbdhndLS/3f/32wufgon4XZEWvwjY/kzFtZcerwme/etdRXsGLV7FkPUUnDLedXFPMuD+0ZupNy25aHFrK7P5f5fGHqnvkZ3oU1nRhZjzkXHTw74tEX2poTHHUWZ7038KQRaiFXIccFaj1Ka8IzakxRtdoYTZG0RmOWK/Txak/aNamyG1OzDYoUqq8tw923BCgIt9037YNMsv2Sswihxw7WCFZ/lDB4lN9frm+BIKQFzBq9bCUm7FgN2WDvew3vbUN3nDj5Q9/wCGS2OW0k528c2/XS0xnF3z0/PCY5JUJj0L78/lI0+daxNz0p5VKRnE6a7rqcqEmQdrNDp7dprRnm1wTnqhj2SoatgeE0YngsxTUNpnSmE1KBSi6aE+0WqYzMzGp9TOFZaa1GWntcrj1LcLhNqqym1IH6eFVRVu6BeHVd8gCSJ+hzu0+ADmEYg6R850y7+9pvEwBd6P9BWGdD9bae/tK6bObD656tyv8VuQm2WDZi1SvqUYHdbKHToTsHclNtC2Z6dnZVirkDuRhRW3yuePP3s9I9520yPueHOZJrt38sm3tjpCil69dy4HqzU2d26DyXNL3NCYTb7jnDc8S4BW7mHHdsLMmNc7DUrRTVavniaN2UZiQv/zm27AzJZHquZzaps4umLdLHq5uTB1zamU4QBGRgv6HIm76z4Os7B184I8Td5DGd6YaMITRWBYAf9f8g7DncYdu8cszqVN7ia9QMDKu3OB59726j9p5hV7xLgJtCtm6G7toWzFTw0QAxJmJjMi7idtR1d/ToLzeoRO23BWJi2BKV6PCh34YvXOU54uk6NzvqPOd5nsBrtZa7Saf3kqZUoFJGLpIK1EJuXCQzWtxYE1Ne5Kpen/j7gTE8PctkY1qchERwksJUabl1mcMLpy1qTBtkkcEKmWFnyKJ7vvr5Mwm38foBcs/vxSMHqgYtfSh0ll8H4ApBEF5Q8P3aGxNYS89PFYkT8jZOS5u26e1BC+7sYfe6T9fdJTu7XtR1J+Uw2g3d9XQuImEzSs6vwmJk4J7LmKeSWwvdm8qP/6W3aVqtZRjGFPEUESw5y8rlOEVxMUPHDX5MKlDz2JEMNyGpr5LWaqXnNNLaEzGaL+SaYs80Fn28qlGdddhNfF12SpGW1Mrj7KkzxU9aPPruV0N5qiQINDZPsGTjwV8+embt9l8xys0RR+Xe+0Ggr5QAEEwQhBfUH/19Zrpv93qeQtJQerr9pID2O7t6F8xsv7OriI1EbEzKQXJJW9dd74buLrxiu3XFKlNOvcosf0dUrmXamBSWRPJjKW4Fg5BEpaYpZ0oFallEKpclObvzq2PrnrlOKcxyOaxNBzjOD0ZmDJc362LKi+xiaZM62zONpWDOzU3q7OaUDLJdZ7q0tVF79iRCaHbmsAhZTOd1gXDBEYgmPPQ2euhtugsBICAgCC+gSNJ3J3KESIQVtLhPVVOXDN1RnqE77153kRyM3fHIWE9Pp9oHnvfaZrOlhMngSAUqqUAt4imS1GPP/HX6dotkQbRCSLEQQrsqml89Ryy+8d+yuippoUau+VlccHDJqb0ZpBs1MAgJn5AJNVLmx5XFWa9829SDzvQIWUzamJk9rDlk4Q7b0U9erj64k3A55OlDRq5+RRrCrRoAABqFdRC2H7rTO5ExZcwflT/dlhPf7g7UwVpDfPLg6O6G7i6LT+DprRqzU2ewVTAZXE/gSQVquTBLHT1VKlBHRQzgsC6ajZJ7/9KTL6xorTmdw8IEJlsWm3+MQqwFWZ5pLPp49U7SlTsmTZypoM4v+xmJUMEuLSGIiKJ1hdKgcZj039857jY1/7pJMVwm44iu8pn7po57+tPE4RPpLg0AEHLCIgh7MnRH/v2pc9/nGyytcXzW7LRYdH4FjZwFd4xU92bzio5Wjm5bSJPDEnkDTyEePiBmjoinjBZmsJm+7clMAhc3VEtrD8eUF8q1xdIajbRW6xnSa1APr5IrbQMGn1VnH7y4M337/TOmqKKpixe/TpdwDLqqqOQMFAaOfvLyPWn86zPapvaMiZd+PUOw9PW7b/q2kN7CAAAhqJ8HIYXQrQe5OEVJOJSYjYnZKJKD4gVIzGFI2CjifI+59siu8p2vbpuf4SDINw+WvXdU4yBIB4s35v/WZE67vuuX6GLlaJ/Ay1YsEvEUkfwUNrPj0zJRU5037aS1mvOd6cr/b+/O46Iq9z+AP2cYGGaGfV9kGUDQwYUEzFARFRMTMw2vZqhpJZmvMpfqVr/KsjS1xHvLenUt2zRN700zK8VcMJQ0wQVEkU1kE5BFYJgZmHPO74+jEw67whyG83n/deaZ58z5znCYz5ztOdw4LJWKwZei42o8FZ1emS5z8iirvzrA+q6T+opUtI+jW7NadXrr6oITP5mxNJHIhse/opwS3//OhSk6nfTowwNatrjKJQ5mtKq6Akc9AcBAPw9CipCvHtKeLmcyDnxTm31WYu/iPeEfre8znrl9w4cjPaViM6nY7P3xgwkhLCET914OiJym70MzTXWa4jZHjrY0t7e2dG8ZePYyha3UW0S1+/FK62rsSwqc8y8551+2L8m3L8lvdWV63ElPvwp/pc7CsrvveuC0ZxI/XLjd3U5/yDO3RnVVxYb6Bu1++qEFniT+0QARRdVpdW/vXZdyNX3sS33iDuA9iNG1MVi5ZRfuXAoAAtTPg5AQUnD54q5l8XN95WEu8pqqzC/f35s9fGLUy5+27FN7ozAg3J+bbqKYG5S2UNRo9kDT7xdf00k0+sCztvSwkrhZW3p0PfDIncy7nXZ5Wc4FWY6FOYxYbHhlum9QUw+N3Oj9wJiKSQtjD3y5IMDWWSo+e1PzS1lTzAc/ZiX98LC9br7y9qaSjUS8KVLxyP79DfGvWDm5d/yavaS6KLf0chpFiQYMfdDWzbunXtZ54PC/Sq8/6Gmvb9HSzPU6dZRz9+4vDwBC0M/HGmVZdkhQ0NbR7t62f3deejzfesEGxZiH9SNHp/539UNe1C2L5kJRYzmldWUlPows87o2aNwiJ7vB3OXndjIfinRy9ztxk9a6ssQlL8u5IMu+pEB/SK9mgF+Fn1J/g4WKgCEaK9v7/RQ6U1Ocl5/yi6ay2C4odND4GWbmFkfemf+yLG+E212XiHyYXnpj+upBUdN7cNE6nU4kEnU8xBpD646tW6zNSol2lzGEOlhc7/hg7JjliT2yn7aqMDtp+ZTPIr0CHa0IIXVa3aqUaxZTXwiZ1Z9vw4Y71BsfxhrlBcYa7Z7Lly8PsCONDrrDospiM3WRSF1CqXOnqPKrZpODZvZyha2lt53Mx9UvXHv5/CqFzwBG6sJICCEXym/9M0c0+cVN7b3yndNYCuyL813ys5zzs/SnsXBpx43GUiSx3L/jw+q8i6KcM1RJdvhz7/sb5Ya6hBD7Af6hc17stBtfP4T+/PzN8Kr0FVMGcg+XhpA3T/1+btfmEU8sv/8Xd/QJmrD+pxfXJ7C3sqXmZrXNZMSiN5VT4u//lQGg/+nnQXjz5s0/owqXWpZ5M1JPRjqAsRzB2Fo1itedE836d6q+G6PU/ZY664vfrz7hR5daaJJLG34s0U7b/Ku+g7SuxvlO2rU+jaXGU5E/cmKFv/KWu0/LK9PrK0r2LRn3QbjLQ0MHEUIqVNrlny5T3ywdMv0ZY34Iei5hk5IOXWi5RcgScrSkbvKQkcYvJufwD/95bJD+IUXI6+EDYvd/0SNBSAhxHTg07otTdJO2Wavu7r0bAUBQ+nkQ+vv7B+8fsDtmYMvGlLIqJ7fgli0iM/HUDXsLzhzdmrLPsroiOChwbYzS9cB3zvlZ9iX5TteyaXML7nZClf5K7jSWSr/BzZ2NtXh227tvDHd46M6RKhe55Ito/4e/Xhs8bRHFx50ZgifP+eGHze5Z5fMGu9w+WeZ0kcuYx4x/gJBlaEsREd99PovM3Ixt0vTsgswsJGYWkp59TQDoZ/p5EHp6eopcvH/JuznV34lrqdPq3k+vGLt+OSFE3KRxKM7n0s6+pMA575JL3iVWJKrxzK+5drXG0+/q2KmVfsFVPgO1snu5lPDGxdRxD991doZUbOZvJ6stu8bLKCcisXncf/44tnX1tv37zFiaksiHz3tjTMxc41dCicw0NMvePSpBM83oiOjMtxvq8y9KXX0CHp7r4h/c7ksAAPSQfn6yDCHkxJWSt154RlySM8FaLK7XNpfXRw4a4UUztsUFre+ZXumnrO+5zaNd8SH7J7jJzO86xSb+SEHI2l/tPHx7ail9UFdOljm69tnHGs/9I+jv21n88/jVlNK6pSGewY7yGyrNZ1dqncbPfXDxO71fbz+Bk2WMDyfL8AIny3QTy0asXHw85zJ744bK2Vnr62s7OUw8eDAJDCQDBxJPTyuGFdOUnY7V0ETHkmaaaGiiplkdQ3QspdaxGpqodUTHEB1L1DpWQ1NqmuiYLv168AiNSrp28rGBf1/BXa/VXa9vHufu02tv2GSMeSnx65XTUv8onOIppRn2x+uqi0VVv/0j1EFqQQgZRmyiFc5zD31fOmaahzKM72IBoD/r70FIUU2rVom9vSlfXyszs9Y/IcxFlLmIWJsbnLJPtZow1MywGvp2UrZMUA1N6RjuWeK+7M0P50SYi6hH/J0pQvJqVMtTro96fkP/G8nlHljIrGZ+diz/zJEfzh2nzMTUIPvolC+5FOSIKOq5wQ7bDn4nhCDUadXp2zeWnDlM65rdQyJDn3odJ/gAGE1/D0JC6LAwYn0vR/g61qUE9XWZ9deZ1a+9+vnBY4Rh3D0HfL7rvyGhYZ0mqFpH7jzLamiqUUdoE9yD3RV+Iyf6jZxICLl0eI+7peGuVGeZRFN9g4+6jEpVXfHjkqgFvtJHRziIRdLfrx/794L/xW4+6OA9sPOZ+yRG13x+zyeFyXu1qnqXoAdCn1ndg6MlAPS4/h+E/HJwcPj351sNGu9hG7SZabW3liU65naCqnXsnZ23ppqgDgP8M2qbDRqzqlQ2/kN4qceYUj95+fVg68mK2+dzxQW6BNlZ/vODhBmfHuW3sHvTrGn8X8K4WCf6nQccrC1c/yzNevf5cWPe+Mo7NIrv0gDahiA0DdwGKCHtJWhHu3C7nqBqHatjiY6hVDq2a4dBe4zboAeSddKTxdWjBzhwLTcaNFuyqh9bzs81l8ZUeuHkpOmBLVuGutioT15mGYaXy2zu07mdibNddM8O8+AeRnk7Kh2tZq1/bt7uK/wWBtAeBGE/dz8J2vbeWqbjBCXMPW2DUhQV+9H+tW/H21/OD3G0LGqkL9Q2j1+9na9BUI3JjJDWd4S2MBMxumZTvAjy+h/7Pxjl1LLFRS4ZIBXdKiu0bXWaGMuyt8oKKZEI+06BRwhCaJs+QVuFZecJWqtiiYhiRaI2E1RDEx1DmpnbCaqhKR1LNDrWysl9xpYjVdevll/LdnL2nDtwqEgsiBPTJXbOxXXqATZ/Xw7U0KRTEzNTTEFCSLNWI7dwMGiUm5s1axoNGrN+/e6vrasVdpYMSwrrtaOWrg+a8LixygT4G4IQepi5iLI2J2IxEYvvaRt0yOBmZlA3EpQmpngtbEthCe+/+NGz/5ng5ySzIIQ0NOleTC4Y8dT/8V3XPXIOCjlTenWs199ZSLNsdlVD6N2DSGT9tr1215rfYgOsLMSEkFpN83NfvGpmLgkYG2vsikHw+v8F9fX19da9cNYodECj0YjFYrHYSD+z2t5b2+HFoNyEccrrivxTv53avMJVQszNRNfrtGHPvK18ZF53X6SPXFBffT3n0EuTt0Z5+9nLCSFamnkr9XpVaNyou8dG2DFb+dMkLxvJ3ytJZaN27sma2d+mGbvi+4AL6nmBC+oBDN3/xaBdT1Du2R7nFzHFL2JKQ9UNurkpysSPljl4Dxy/7sfn1j1rqSm1lphfv9U4dPaLD855qWUflmHMdJqWKUgIcZZJaFWtcYsFIARBCELWq8Mp3MOARFaObvf3hvoKt6CQf3z9l1ZVp1XVjXcZ0LoDJRI104Y/KFhCdCa4gwr6AQQhQLf1doLefTEoS5tmOkjkNhK5TXvP2imUf5XWhnv8PYDO0Ws33YJHGaU0gLsgCAGM594S1OBi0PYStKaBJeaUqQynEPnKZytfiE4I1E72sWcJ+3N+zTcFjTM+a/dW2AC9B0EI0Nd1djHo7Yn6etra+q4L8PvycAq27j5zvk0/+t3GHWlHRSIzjwdnPvHuCnFn9/gE6A0IQoB+y+gDEnVvOAULmfVDCe8S8u49vDWAHoQgBABDvTcgkcHFoFyCNjT34X24IAAIQgDoMfczIFG3ElQ/IBESFO4fghAA+HefCdrGDVj6+4BE0IMQhABgwto5EZfcW4LyMpwC8A5BCABC1CMDEtWpaGJuoWWozhK08+EUgEcIQgCAbmiZoLaEyGQURVG9NJyCCSUoy9DluZl15cUOXv5OvoP4Lqd7+lAQvv3226dOneKmHR0dd+3axW89AAA9qB8PSFR2Of3YmgWD5MTXSpxZo62wcJj07vcmdI/JPhSEFy5cCAkJmTx5MiFEIjHJO7EBAPS4HhmQyOBi0PYTtNsDEqnrqn9/I277RF+vO/fUPHvj1ssrYuduP0eJzLr3VnnSh4KQEKJUKqOjo/muAgDA5HVxQKLWujucwokfdy4IsPVqcWfpMDfbUJu6wnMpvqHjeuOt9bi+FYRbtmz5/vvvAwMDV61apVAo+C4HAEBwujucQl5tltLB0qAxxJZqVufO9Bvf4wMS9QajBmFmZmZZWZlBo52dXXh4OCFk1qxZdnZ2crn8hx9+GDlyZEZGhptb23el0Wq1gwcPpqjbf4+YmJiNGze2t1CVSqXvCcZh5BvzAgeruvE1NjYyDCPwj93azr7yWpNB481maoSzlY+5qiuvoNaxzaxIS7NNDGlmKbWONNGkmaWaGKKliZahtDTbRJMmhjSxlEbHNmsaGxq6Wp6lpWWn30VGvUP9xo0bjx49atCoVCo/+ugjg8aRI0cuWLBg6dKlbb6OVCpNS0uztLz9G8TV1VUul7e3UNyh3vgQhLzAqm58KpVKJpMJPAizsrISZk797uFAsej251Cv1c1JyjmRftHW1rY3ltjjq7pRv6pefvnll19+uSs9XV1d6+vr23uWoiiFQiGVYqB6AACeKZXKuGeWzN362aJAOy8b2ZVq1Tc5NWs2/auXUrA3iDrvYhRNTU2pqanc9OHDh48cOTJhwgR+SwIAgK5YtuqVr/YfLB4Ws5vx1I6d9dvJM49Of4zvorqhr+y80ul08fHxZWVlUqnU0tJyy5YtI0eO5LsoAADokqCgoNXvr+O7invUV4JQJpPl5eWpVCqapm1sbPguBwAAhKKvBCGng3NeAAAAekNfOUYIAADACwQhAAAIGoIQAAAErW8dIwQAACFrbGz87ttvMv467ejqNm1mXFhYmBEWii1CAADoEy5cuBARMrR096eTGq74Zh5+Y+GcFUuXGGG52CIEAAD+sSz7zJOzP4nw9LWTcS2T/JxX/nH05/0/TXt0eq8uGluEAADAv+zsbE8J0acg55nBTj98va23F40gBAAA/t28edNFam7Q6CKXVFRU9PaiEYQAAMA/hUKRU6M2aMyuahgYFNTbi0YQAgAA/zw9PW29FL/l39S31Gt1my6UJyxb0duLxskyAADQJ3y9a8+z8+b+92jOA47SmmZyurzhnQ2Jw4YN6+3lIggBAKBPsLOz2/Pzr7m5uZmZmY6OjomhoTKZrPPZ7huCEAAA+pCAgICAgABjLhHHCAEAQNAQhAAAIGgIQgAAEDQEIQAACBqCEAAABA1BCAAAgoYgBAAAQUMQAgCAoCEIAQBA0BCEAAAgaAhCAAAQNAQhAAAIGoIQAAAEDUEIAACChiAEAABBQxACAICgIQgBAEDQEIQAACBoCEIAABA0BCEAAAgaghAAAAQNQQgAAIKGIAQAAEFDEAIAgKAhCAEAQNAQhAAAIGgIQgAAEDQEIQAACBqCEAAABA1BCAAAgoYgBAAAQUMQAgCAoCEIAQBA0HgIwtTU1Pnz5w8bNiwuLq5le2ZmZlhYmFQqHTJkSGpqqvELAwAAAeIhCBsbG8PDwydPnlxaWtqy/cknn4yLi1OpVCtXrpw1a5ZOpzN+bQAAIDQ8BOHEiRNfeOEFpVLZsjEtLe3atWvLly8XiUQLFy4Ui8VJSUnGrw0AAISmrxwjzMnJCQwMlEgk3MPg4OCrV6920L+2trbmDmw7AgDAPRP3xoteuHDhxIkTBo0ikWjp0qXtzVJTUyOXy/UPbWxsqqur2+us1WpbblBOnz79448/bq9zQ0NDl4qGnqPRaMRisVjcK2sXtAeruvE1NjbSNE1RFN+FCEu3VnVLS0tzc/OO+/TKV9WtW7fy8/MNGs3MzDqYxdHRsa6uTv+wtrbW2dm5vc4SiaS0tFQqlXaxHmtr6y72hB5hbm6OIOQFVnUjE4lEMpkMQWh8Pbuq98pXVWRkZGRkZLdmCQoKunr1qlqtlkqlLMtmZGQsX768N2oDAABoiYdjhLdu3UpLSyssLGxoaEhLS8vLyyOEDB8+PDg4eM2aNQ0NDZs3b5ZIJBMnTjR+bQAAIDQ8BGFmZmZCQsKBAwcsLCwSEhI+++wzrn3nzp1nz5719fXdu3fvvn37Ot6VCgAA0CMolmX5rqHbZDJZVVVVF48R1tfX48CJkeFkGV5gVTc+lUqFY4TG1+Orel+5fAIAAIAXCEIAABA0BCEAAAgaghAAAAQNQQgAAIKGIAQAAEFDEAIAgKAhCAEAQNAQhAAAIGgIQgAAEDQEIQAACBqCEAAABA1BCAAAgoYgBAAAQUMQAgBAbzlz5swjUWMfCPQLUwauenHprVu3+K6oDQhCAADoFbt37Vw1f/arPuTnRwL/N0nhk5c6flR4bW0t33UZQhACAEDPYxjmvTde+88Ef397OSHEjKJmDHSJ95UnbviA79IMIQgBAKDn5ebm+ttJrSzELRsn+zr8ceR3vkpqD4IQAAB6Hk3TYhFl0GhGUTqdjpd6OoAgBACAnhcQEHD5ZkMTzbRsTC6qGTVmLF8ltQdBCAAAPc/c3HzpipXLkvOr1U1cy6nimk+vVK947Q1+C2tN3HkXAACA7lvywjI3zwGLV79Faxp1LDt46LADR//r5ubGd12G+v8W4cGDB/kuQXAuXrxYVFTEdxXC0tzcfOTIEb6rEJw///yzqqqK7yr6tBkzH//z4qXUrJxzOde+/3G/l5fXfb5gfX19SkpKj9Sm18+DUKPRzJ8/n+8qBGfbtm2//vor31UIS2Fh4cqVK/muQnA2b9588uRJvqswAWJxj+19TE9PX7t2bU+9GqefByHwhWVZvksAMAas6v0AghAAAAQNQQgAAIJGmeJ2vbm5eWRkpEjUeYozDJOcnDx+/HgjVAV6V65csbKyGjBgAN+FCIharT537lxERATfhQjLhQsXPDw8nJ2d+S5EQGpra3Nzc8PCwrrYf8aMGc8//3zpJAVHAAALiElEQVTHfUwyCHfu3Nn1Na+goEChUPRqPWCgsrJSKpVaWVnxXYiAsCxbWFjo6+vLdyHCUlpa6ujoKJFI+C5EQHQ6XVlZWdfPPlUoFP7+/h33MckgBAAA6Ck4RggAAIKGIAQAAEFDEAIAgKAhCAEAQNBMftDthoaGgwcP6nS6mJgYOzu71h1Ylk1OTi4oKIiIiAgKCtK3V1ZWJiUlyWSymJgYqVRqxJL7g0uXLp05cyYgIGDs2DbuqNLU1HTy5MmioiIvL69x48ZxF7owDHP06FF9H29v78DAQONVbPqqq6sPHTpkYWERExMjl8sNnq2trT179qz+4dChQ11dXbnpoqKiY8eOOTk5TZo0ydzc3HgV9wtnz57NyMgYMmRIeHh462ePHDnS8nxDLy+voKAgtVrdcty1gQMH+vj4GKPW/qK4uPjq1avBwcH6ddhARUXF4cOH5XJ5TEyMpaWlvj0lJSUnJycsLGzo0KHdWyRryiorK/39/WNjY2fNmuXh4VFYWNi6z/z585VK5eLFi52dnXft2sU1ZmVlOTk5xcfHR0dHDxs2rL6+3riFm7avvvrKxcUlISEhKChoyZIlrTsoFIqIiIinnnpKqVSOHj1arVazLKvRaAgh48ePj46Ojo6O/vTTT41euAnLy8tzdXWdPXv2lClTgoKCampqDDokJydbWVlF33Hs2DGu/Y8//nBwcFi0aNGoUaPGjx+v0+mMXbopW7t2rZeXV0JCgre393vvvde6Q0xMjP4zt7S0XLduHcuyeXl5YrFY3757926jF27Chg4dam1tLZFI9F/XBjIzMx0dHefNmzdhwoSQkJCGhgaufdmyZQEBAQkJCW5ubp9//nm3FmraQbhmzZrY2FhueuHChcuWLTPocPHiRVtb26qqKpZl9+3b5+fnR9M0y7Lx8fErVqxgWZam6TFjxmzZssW4hZuw5uZmDw+PQ4cOsSxbXl4ul8tzcnIM+uTm5nITarXax8dnx44d7J0gxG+Oe/Pcc88lJCRw05MnT96wYYNBh+Tk5CFDhrSecfz48Zs2bWJZVqvVBgYG7tu3r7dL7TdqamrkcnlmZibLspcvX5bJZNXV1e11LigoEIvF169fZ1k2Ly/Pzs7OeIX2L/n5+TRNDx8+vL0gfOKJJ1555RWWZWmajoiI4DLv2rVrUqm0pKSEZdnk5GRnZ2eNRtP1hZr2McIDBw7ExcVx03FxcQcOHGjdYcKECQ4ODoSQRx55pKysLCsrq+WMIpFo5syZrWeE9qSnpzc2NkZHRxNCXFxcxo4d+8svvxj00V++amlpaW9v39TUpH/q5MmTx48fr62tNVrB/cOBAwcef/xxbvrxxx9vc43VarVJSUmnT5/mfnMQQlQq1fHjx7kZLSwspk2bhlW9644dO+bl5RUcHEwIGTRokJ+f3++//95e523btj388MP6q7y5AwEpKSkNDQ1GKre/UCgUHY8apv9fEIlEM2bM4Fbp3377beTIkR4eHoSQsWPHUhR1+vTpri/UtIOwpKTE09OTm/b09OR+Dhh00A/0ZW5u7uLiUlJS0tjYWFtbazCjMcs2aSUlJe7u7vo1teNPb//+/WVlZdOmTeMeurm5ffzxx6+//rpCodi7d68xyu0XGIa5ceOGfk1u7zMXiUSffPLJokWLlEplZmYmIaS0tJQQwn07dDAjtKnltwfp8NNjGOabb75ZtGiRvsXBwWHz5s3czroTJ070eq2CUVdXV19f3/p/oeUfi6IoDw+Pbq3qpn2yDE3T+m9kMzMzmqYNOuh0upajH4nFYp1Ox3WjKEo/o06nM0q9/QFN0/qPjnT46aWlpT377LPbt293dHQkhEgkkuLiYjMzM0LIt99+u3DhwqlTp1pYWBinbJPGMAzDMB2vsaNHj75y5QohhGXZF1544cUXXzx69Cj3x8Kqfm8MVnXu26PNnklJSSqVKjY2lnvo6+ubn5/Pzbtu3bqnn346JyfHCAULQXvf3l3/Y7XJtLcI3d3dKyoquOny8nJ3d/eWnwUhxMPDQ9+BYZiKigoPDw9ra2u5XF5ZWamfUf+TGTrl7u6u/+jInY+9dbeMjIzY2NjPP/980qRJ+kYuBQkhc+bMqaurKygo6O1q+wexWOzs7NzxGqv/bCmKmjNnzvnz5wkhbm5uDMPcvHlTP2ObfyxoU8uvF9LhF8W2bdsWLFig/80tEon0X0Rz587Nzc3FDtKeYm9vb2lp2fp/oet/rDaZdhBGRUUdOnSIm05KSoqKiuKma2pqmpubuQ7Hjx/njlGlpqZKpVJuj39UVFRSUlLrGaFTISEhTU1N6enphBC1Wn3ixAnu5h5NTU36I3/Z2dlTpkzZtGnTY4891uaLnDt3TiQS6fdOQ6faW9Wrqqpa7whJT0/nDlbZ2dmFhIRwqzrLsocPH8adWLpuzJgxV65cKSsrI4SUl5dnZGRwFwtpNJq6ujp9t6qqqv379y9YsKDNF0lLS3NwcMAA9PdJq9XeunWLm27zfyEqKurUqVPcD46srKyampo2L3dp132f48On/Px8e3v7VatWvfXWW7a2thcuXODanZycfvrpJ246MjIyJiZm06ZN/v7+Gzdu5BpPnDhhY2Pz3nvvLV261M3Nrby8nJ83YJpWr14dFBSUmJg4ceLEmJgYrnHnzp0+Pj7ctKen58CBAxffwf0tduzY8eSTT37wwQevvvqqk5PTm2++yVf9pujs2bM2NjarV69esWKFk5MTd3aiVqslhJw9e5Zl2ddee23JkiUbN25MSEiwsrLSnx26Z88eZ2fnDRs2PPnkk4GBgY2NjXy+DVOzePHi0NDQzZs3h4eHP/3001zjhx9+GBERoe+TmJj44IMPtpwrMTFx0aJF69evf+mll2xsbHClULds3rx58eLFjo6OEydOXLx4cXZ2NsuyW7duVSqVXIfjx4/b2tq+//77zz//PLeDimufPn36uHHjEhMTlUrl66+/3q2FmvzdJwoKCnbs2EHT9OzZswcNGsQ1bt++PTIy0tvbmxCiVqu//vrroqKi0aNHT506VT/j+fPn9+7dK5VK58+fj12j3bVv377Tp08rFAr9HqH8/PwzZ87MmTOHEPLFF18wDKPvPGLEiLCwsBs3bvz8888FBQXW1tZjx44dM2YMb9WbpkuXLu3Zs8fCwiI+Pp5btxmG+fLLL2fMmOHk5HT58uWkpKSysjIXF5fY2NiWgxWkpKQcPHjQwcHhqaee4s6ghi5iGGbnzp3cBfVPPPEEt/85IyPj2rVr+lPADh486Ojo2HL7o7Cw8JdffikqKrK3t580adIDDzzAT/Wm6ddffy0uLtY/jI2N9fDwyM7OzsjI0F8jcO7cub1798rl8nnz5um/vZuamr799tu8vLzw8PAZM2YYHCbrmMkHIQAAwP0w7WOEAAAA9wlBCAAAgoYgBAAAQUMQAgCAoCEIAQBA0BCEAAAgaAhCAAAQNAQhAAAIGoIQAAAEDUEIAACChiAEMGF79+51cHDYvn27vuWdd95xcXHhbg8CAF2BsUYBTNvChQv37Nnz119/DR48+Pjx49HR0WvWrHnttdf4rgvAZCAIAUybWq0eNWoUwzAHDhx46KGHBg0adPjwYf19egGgUwhCAJN36dKlkSNHisViuVx+/vx5FxcXvisCMCU4Rghg8oKDgydMmFBXV7d06VKkIEB3YYsQwOTt2LEjPj5+6NChhYWF6enp/v7+fFcEYEoQhACmLTc3NzQ0dObMmR9//HFYWJiVldXJkyclEgnfdQGYDAQhgAnTarURERH19fVpaWnW1tYXL14cNWpUQkJCYmIi36UBmAwcIwQwYStXrszKytq9e7e1tTUhZNiwYevXr//Xv/61b98+vksDMBnYIgQAAEHDFiEAAAgaghAAAAQNQQgAAIKGIAQAAEFDEAIAgKAhCAEAQNAQhAAAIGj/DwTUgrnmkUBnAAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Plots, Distributions, LaTeXStrings, LinearAlgebra\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), ribbon=sqrt.(v(x_test)), label=L\"f(x)\") # 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", "scatter!(x, y, xlabel=\"x\", ylabel=\"y\", label=L\"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, color=:red, label=\"LS\") # plot LS solution\n", "\n", "# Weighted LS regression\n", "W = Diagonal(1 ./ v(x)) # weight matrix\n", "w_wls = inv(_x'*W*_x) * _x' * W * y\n", "plot!(x_test, _x_test*w_wls, color=:green, label=\"WLS\") # plot WLS solution\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Some Final Remarks\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Aside from the above example, we also recommend that you read through the [RxInfer Bayesian Linear Regression example](https://reactivebayes.github.io/RxInfer.jl/stable/examples/basic_examples/Bayesian%20Linear%20Regression%20Tutorial/). \n", "- In this lesson, we focussed on modelling the map from given inputs $x$ to uncertain outputs $y$, or more formally, on the distribution $p(y|x)$. What if you want to fit the best curve through a data set $\\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$ where both variables $x_n$ and $y_n$ are subject to errors? In other words, we must now also fit a model $p(x)$ for the inputs, leading to a generative model $p(y,x) = p(y|x) p(x)$. \n", "- While this is a very common problem that occurs throughout the sciences, a proper solution to this problem is still hardly covered in statistics textbooks. Edwin Jaynes discusses a fully Bayesian solution in his 1990 paper on [Straight Line Fitting - A Bayesian Solution](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Jaynes-1990-straight-line-fitting-a-Bayesian-solution.pdf). (Optional reading)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##
OPTIONAL SLIDES
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Some Useful Matrix Calculus\n", "\n", "When doing derivatives with matrices, e.g. for maximum likelihood estimation, it will be helpful to be familiar with some matrix calculus. We shortly recapitulate used formulas here. \n", "\n", "- We define the **gradient** of a scalar function $f(A)$ w.r.t. an $n \\times k$ matrix $A$ as\n", "$$\n", "\\nabla_A f \\triangleq\n", " \\begin{bmatrix}\n", "\\frac{\\partial{f}}{\\partial a_{11}} & \\frac{\\partial{f}}{\\partial a_{12}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{1k}}\\\\\n", "\\frac{\\partial{f}}{\\partial a_{21}} & \\frac{\\partial{f}}{\\partial a_{22}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{2k}}\\\\\n", "\\vdots & \\vdots & \\cdots & \\vdots\\\\\n", "\\frac{\\partial{f}}{\\partial a_{n1}} & \\frac{\\partial{f}}{\\partial a_{n2}} & \\cdots & \\frac{\\partial{f}}{\\partial a_{nk}}\n", " \\end{bmatrix}\n", "$$\n", " \n", "\n", " \n", "- The following formulas are useful (see Bishop App.-C)\n", "$$\\begin{align*}\n", "|A^{-1}|&=|A|^{-1} \\tag{B-C.4} \\\\\n", "\\nabla_A \\log |A| &= (A^{T})^{-1} = (A^{-1})^T \\tag{B-C.28} \\\\\n", "\\mathrm{Tr}[ABC]&= \\mathrm{Tr}[CAB] = \\mathrm{Tr}[BCA] \\tag{B-C.9} \\\\\n", "\\nabla_A \\mathrm{Tr}[AB] &=\\nabla_A \\mathrm{Tr}[BA]= B^T \\tag{B-C.25} \\\\\n", "\\nabla_A \\mathrm{Tr}[ABA^T] &= A(B+B^T) \\tag{B-C.27}\\\\\n", " \\nabla_x x^TAx &= (A+A^T)x \\tag{from B-C.27}\\\\\n", "\\nabla_X a^TXb &= \\nabla_X \\mathrm{Tr}[ba^TX] = ab^T \\notag\n", "\\end{align*}$$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Derivation Predictive Distribution\n", "\n", "- In the [derivation of the predictive distribution](#predictive-distribution), we replaced $\\mathcal{N}(w\\,|\\,m_N,S_N)\\,\\mathrm{d}w$ with $\\mathcal{N}(z\\,|\\,x_\\bullet^T m_N,x_\\bullet^T S_N x_\\bullet)\\,\\mathrm{d}z$. Here we discuss why that is allowed. \n", "\n", "- Since $z = x^T w$ (drop the bullet for notational simplicity), we have\n", "$$p(z) = \\mathcal{N}(z|m_z,\\Sigma_z)$$\n", "with\n", "$$\\begin{aligned} m_z &:= E[z] = E[x^T w] = x^T E[w] = x^T m_N \\\\ \\Sigma_z &:= E[(z-m_z)(z-m_z)^T] \\\\   &= E[(x^T w - x^T m_N)(x^T w - x^T m_N)^T] \\\\   &= x^T E[(w - m_N)(w - m_N)^T]x \\\\   &= x^T S_N x \\end{aligned}$$\n", "\n", "- Then we equate probability masses in both domains:\n", "$$ \\mathcal{N}(z|m_z,\\Sigma_z)\\mathrm{d}z = \\mathcal{N}(w|m_N,S_N)\\mathrm{d}w$$\n", "\n", "$$ \\Rightarrow \\mathcal{N}(z|x^T m_N,x^T S_N x)\\mathrm{d}z = \\mathcal{N}(w|m_N,S_N)\\mathrm{d}w $$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "open(\"../../styles/aipstyle.html\") do f\n", " display(\"text/html\", read(f, String))\n", "end\n" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "anaconda-cloud": {}, "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Julia 1.10.2", "language": "julia", "name": "julia-1.10" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.10.5" } }, "nbformat": 4, "nbformat_minor": 4 }