{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Factor Graphs" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "### Preliminaries\n", "\n", "- Goal \n", " - Introduction to Forney-style factor graphs and message passing-based inference \n", "- Materials \n", " - Mandatory\n", " - These lecture notes \n", " - Loeliger (2007), [The factor graph approach to model based signal processing](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Loeliger-2007-The-factor-graph-approach-to-model-based-signal-processing.pdf), pp. 1295-1302 (until section V)\n", " - Optional\n", " - Frederico Wadehn (2015), [Probabilistic graphical models: Factor graphs and more](https://www.youtube.com/watch?v=Fv2YbVg9Frc&t=31) video lecture (**highly recommended**)\n", " - References\n", " - Forney (2001), [Codes on graphs: normal realizations](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Forney-2001-Codes-on-graphs-normal-realizations.pdf)\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Why Factor Graphs?\n", "\n", "- A probabilistic inference task gets its computational load mainly through the need for marginalization (i.e., computing integrals). E.g., for a model $p(x_1,x_2,x_3,x_4,x_5)$, the inference task $p(x_2|x_3)$ is given by \n", "$$\\begin{align*}\n", "p(x_2|x_3) = \\frac{p(x_2,x_3)}{p(x_3)} = \\frac{\\idotsint p(x_1,x_2,x_3,x_4,x_5) \\, \\mathrm{d}x_1 \\mathrm{d}x_4 \\mathrm{d}x_5}{\\idotsint p(x_1,x_2,x_3,x_4,x_5) \\, \\mathrm{d}x_1 \\mathrm{d}x_2 \\mathrm{d}x_4 \\mathrm{d}x_5}\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Since these computations (integrals or sums) suffer from the \"curse of dimensionality\", we often need to solve a simpler problem in order to get an answer. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Factor graphs provide a computationally efficient approach to solving inference problems **if the probabilistic model can be factorized**. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Factorization helps. For instance, if $p(x_1,x_2,x_3,x_4,x_5) = p(x_1)p(x_2,x_3)p(x_4)p(x_5|x_4)$, then\n", "\n", "$$\\begin{align*}\n", "p(x_2|x_3) &= \\frac{\\idotsint p(x_1)p(x_2,x_3)p(x_4)p(x_5|x_4) \\, \\mathrm{d}x_1 \\mathrm{d}x_4 \\mathrm{d}x_5}{\\idotsint p(x_1)p(x_2,x_3)p(x_4)p(x_5|x_4) \\, \\mathrm{d}x_1 \\mathrm{d}x_2 \\mathrm{d}x_4 \\mathrm{d}x_5} \n", " = \\frac{p(x_2,x_3)}{\\int p(x_2,x_3) \\mathrm{d}x_2}\n", "\\end{align*}$$\n", "\n", "which is computationally much cheaper than the general case above." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In this lesson, we discuss how computationally efficient inference in *factorized* probability distributions can be automated by message passing-based inference in factor graphs." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Factor Graph Construction Rules\n", "\n", "- Consider a function \n", "$$\n", "f(x_1,x_2,x_3,x_4,x_5) = f_a(x_1,x_2,x_3) \\cdot f_b(x_3,x_4,x_5) \\cdot f_c(x_4)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The factorization of this function can be graphically represented by a **Forney-style Factor Graph** (FFG):\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- An FFG is an **undirected** graph subject to the following construction rules ([Forney, 2001](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Forney-2001-Codes-on-graphs-normal-realizations.pdf))\n", "\n", " 1. A **node** for every factor;\n", " 1. An **edge** (or **half-edge**) for every variable;\n", " 1. Node $f_\\bullet$ is connected to edge $x$ **iff** variable $x$ appears in factor $f_\\bullet$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- A **configuration** is an assigment of values to all variables.\n", "\n", "- A configuration $\\omega=(x_1,x_2,x_3,x_4,x_5)$ is said to be **valid** iff $f(\\omega) \\neq 0$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Equality Nodes for Branching Points\n", "\n", "\n", "- Note that a variable can appear in maximally two factors in an FFG (since an edge has only two end points)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Consider the factorization (where $x_2$ appears in three factors) \n", "\n", "$$\n", " f(x_1,x_2,x_3,x_4) = f_a(x_1,x_2)\\cdot f_b(x_2,x_3) \\cdot f_c(x_2,x_4)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For the factor graph representation, we will instead consider the function $g$, defined as\n", "$$\\begin{align*}\n", " g(x_1,x_2&,x_2^\\prime,x_2^{\\prime\\prime},x_3,x_4) \n", " = f_a(x_1,x_2)\\cdot f_b(x_2^\\prime,x_3) \\cdot f_c(x_2^{\\prime\\prime},x_4) \\cdot f_=(x_2,x_2^\\prime,x_2^{\\prime\\prime})\n", "\\end{align*}$$\n", " where \n", "$$\n", "f_=(x_2,x_2^\\prime,x_2^{\\prime\\prime}) \\triangleq \\delta(x_2-x_2^\\prime)\\, \\delta(x_2-x_2^{\\prime\\prime})\n", "$$\n", "\n", "

\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Equality Nodes for Branching Points, cont'd\n", "\n", "- Note that through introduction of auxiliary variables $X_2^\\prime$ and $X_2^{\\prime\\prime}$ and a factor $f_=(x_2,x_2^\\prime,x_2^{\\prime\\prime})$ each variable in $g$ appears in maximally two factors. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The constraint $f_=(x,x^\\prime,x^{\\prime\\prime})$ enforces that $X=X^\\prime=X^{\\prime\\prime}$ **for every valid configuration**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Since $f$ is a marginal of $g$, i.e., \n", "$$\n", "f(x_1,x_2,x_3,x_4) = \\iint g(x_1,x_2,x_2^\\prime,x_2^{\\prime\\prime},x_3,x_4)\\, \\mathrm{d}x_2^\\prime \\mathrm{d}x_2^{\\prime\\prime}\n", "$$\n", "it follows that any inference problem on $f$ can be executed by a corresponding inference problem on $g$, e.g.,\n", "$$\\begin{align*}\n", "f(x_1 \\mid x_2) &\\triangleq \\frac{\\iint f(x_1,x_2,x_3,x_4) \\,\\mathrm{d}x_3 \\mathrm{d}x_4 }{ \\idotsint f(x_1,x_2,x_3,x_4) \\,\\mathrm{d}x_1 \\mathrm{d}x_3 \\mathrm{d}x_4} \\\\\n", " &= \\frac{\\idotsint g(x_1,x_2,x_2^\\prime,x_2^{\\prime\\prime},x_3,x_4) \\,\\mathrm{d}x_2^\\prime \\mathrm{d}x_2^{\\prime\\prime} \\mathrm{d}x_3 \\mathrm{d}x_4 }{ \\idotsint g(x_1,x_2,x_2^\\prime,x_2^{\\prime\\prime},x_3,x_4) \\,\\mathrm{d}x_1 \\mathrm{d}x_2^\\prime \\mathrm{d}x_2^{\\prime\\prime} \\mathrm{d}x_3 \\mathrm{d}x_4} \\\\\n", " &= g(x_1 \\mid x_2)\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- $\\Rightarrow$ **Any factorization of a global function $f$ can be represented by a Forney-style Factor Graph**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Probabilistic Models as Factor Graphs\n", "\n", "- FFGs can be used to express conditional independence (factorization) in probabilistic models. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For example, the (previously shown) graph for \n", "$f_a(x_1,x_2,x_3) \\cdot f_b(x_3,x_4,x_5) \\cdot f_c(x_4)$ \n", "could represent the probabilistic model\n", "$$\n", "p(x_1,x_2,x_3,x_4,x_5) = p(x_1,x_2|x_3) \\cdot p(x_3,x_5|x_4) \\cdot p(x_4)\n", "$$\n", "where we identify \n", "$$\\begin{align*}\n", "f_a(x_1,x_2,x_3) &= p(x_1,x_2|x_3) \\\\\n", "f_b(x_3,x_4,x_5) &= p(x_3,x_5|x_4) \\\\\n", "f_c(x_4) &= p(x_4)\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This is the graph\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Inference by Closing Boxes\n", "\n", "- Factorizations provide opportunities to cut on the amount of needed computations when doing inference. In what follows, we will use FFGs to process these opportunities in an automatic way by message passing between the nodes of the graph. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume we wish to compute the marginal\n", "$$\n", "\\bar{f}(x_3) \\triangleq \\sum_{x_1,x_2,x_4,x_5,x_6,x_7}f(x_1,x_2,\\ldots,x_7) \n", "$$\n", "for a model $f$ with given factorization \n", "$$\n", "f(x_1,x_2,\\ldots,x_7) = f_a(x_1) f_b(x_2) f_c(x_1,x_2,x_3) f_d(x_4) f_e(x_3,x_4,x_5) f_f(x_5,x_6,x_7) f_g(x_7)\n", "$$ \n", "- The FFG for $f$ is (we will discuss the usage of directed edges below):\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Due to the factorization and the [Generalized Distributive Law](https://en.wikipedia.org/wiki/Generalized_distributive_law), we can decompose this sum-of-products to the following product-of-sums:\n", "$$\\begin{align*}\\bar{f}&(x_3) = \\\\\n", " &\\underbrace{ \\Bigg( \\sum_{x_1,x_2} \\underbrace{f_a(x_1)}_{\\overrightarrow{\\mu}_{X_1}(x_1)}\\, \\underbrace{f_b(x_2)}_{\\overrightarrow{\\mu}_{X_2}(x_2)}\\,f_c(x_1,x_2,x_3)\\Bigg) }_{\\overrightarrow{\\mu}_{X_3}(x_3)} \n", " \\underbrace{ \\cdot\\Bigg( \\sum_{x_4,x_5} \\underbrace{f_d(x_4)}_{\\overrightarrow{\\mu}_{X_4}(x_4)}\\,f_e(x_3,x_4,x_5) \\cdot \\underbrace{ \\big( \\sum_{x_6,x_7} f_f(x_5,x_6,x_7)\\,\\underbrace{f_g(x_7)}_{\\overleftarrow{\\mu}_{X_7}(x_7)}\\big) }_{\\overleftarrow{\\mu}_{X_5}(x_5)} \\Bigg) }_{\\overleftarrow{\\mu}_{X_3}(x_3)}\n", "\\end{align*}$$\n", "which is computationally (much) lighter than executing the full sum $\\sum_{x_1,\\ldots,x_7}f(x_1,x_2,\\ldots,x_7)$\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that the new factor $\\overrightarrow{\\mu}_{X_3}(x_3)$ is obtained by multiplying all enclosed factors ($f_a$, $f_b, f_c$) by the red dashed box, followed by marginalization (summing) over all enclosed variables ($x_1$, $x_2$).\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This is the **Closing the Box**-rule, which is a general recipe for marginalization of latent variables (inside the box) and leads to a new factor that has the variables (edges) that cross the box as arguments. For instance, the argument of the remaining factor $\\overrightarrow{\\mu}_{X_3}(x_3)$ is the variable on the edge that crosses the red box ($x_3$).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Hence, $\\overrightarrow{\\mu}_{X_3}(x_3)$ can be interpreted as a **message from the red box toward variable** $x_3$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We drew _directed edges_ in the FFG in order to distinguish forward messages $\\overrightarrow{\\mu}_\\bullet(\\cdot)$ (in the same direction as the arrow of the edge) from backward messages $\\overleftarrow{\\mu}_\\bullet(\\cdot)$ (in opposite direction). This is just a notational convenience since an FFG is computationally an undirected graph. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sum-Product Algorithm\n", "\n", "- Closing-the-box can also be interpreted as a **message update rule** for an outgoing message from a node. For a node $f(y,x_1,\\ldots,x_n)$ with incoming messages $\\overrightarrow{\\mu}_{X_1}(x_1), \\overrightarrow{\\mu}_{X_1}(x_1), \\ldots,\\overrightarrow{\\mu}_{X_n}(x_n)$, the outgoing message is given by ([Loeliger (2007), pg.1299](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Loeliger-2007-The-factor-graph-approach-to-model-based-signal-processing.pdf)): \n", "\n", "$$ \\boxed{\n", "\\underbrace{\\overrightarrow{\\mu}_{Y}(y)}_{\\substack{ \\text{outgoing}\\\\ \\text{message}}} = \\sum_{x_1,\\ldots,x_n} \\underbrace{\\overrightarrow{\\mu}_{X_1}(x_1)\\cdots \\overrightarrow{\\mu}_{X_n}(x_n)}_{\\substack{\\text{incoming} \\\\ \\text{messages}}} \\cdot \\underbrace{f(y,x_1,\\ldots,x_n)}_{\\substack{\\text{node}\\\\ \\text{function}}} }\n", "$$\n", "\n", "

\n", "\n", "- This is called the **Sum-Product Message** (SPM) update rule. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that all SPM update rules can be computed from information that is **locally available** at each node." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If the factor graph for a function $f$ has **no cycles** (i.e., the graph is a tree), then the marginal $\\bar{f}(x_3) = \\sum_{x_1,x_2,x_4,x_5,x_6,x_7}f(x_1,x_2,\\ldots,x_7)$ is given by multiplying the forward and backward messages on that edge:\n", "\n", "$$ \\boxed{\n", "\\bar{f}(x_3) = \\overrightarrow{\\mu}_{X_3}(x_3)\\cdot \\overleftarrow{\\mu}_{X_3}(x_3)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- It follows that the marginal $\\bar{f}(x_3) = \\sum_{x_1,x_2,x_4,x_5,x_6,x_7}f(x_1,x_2,\\ldots,x_7)$ can be efficiently computed through sum-product messages. Executing inference through SP message passing is called the **Sum-Product Algorithm**." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Message Passing Schedules\n", "\n", "- In a tree graph, start with messages from the terminals and keep passing messages through the internal nodes towards the \"target\" variable ($x_3$ in above problem) until you have both the forward and backward message for the target variable. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In a tree graph, if you continue to pass messages throughout the graph, the Sum-Product Algorithm computes **exact** marginals for all hidden variables." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- If the graph contains cycles, we have in principle an infinite tree by \"unrolling\" the graph. In this case, the SP Algorithm is not guaranteed to find exact marginals. In practice, if we apply the SP algorithm for just a few iterations (\"unrolls\"), then we often find satisfying approximate marginals. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Terminal Nodes and Processing Observations \n", "\n", "- We can use terminal nodes represent observations, e.g., add a factor $f(𝑦)=𝛿(π‘¦βˆ’3)$ to terminate the half-edge for variable $π‘Œ$ if $𝑦=3$ is observed.\n", "\n", "

\n", "\n", "- Terminal nodes that carry observations are denoted by small black boxes." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", " \n", "- The message out of a **terminal node** (attached to only 1 edge) is the factor itself. For instance, closing a box around terminal node $f_a(x_1)$ would lead to $$\\overrightarrow{\\mu}_{X_1}(x_1) \\triangleq \\sum_{ \\stackrel{ \\textrm{enclosed} }{ \\textrm{variables} } } \\;\\prod_{\\stackrel{ \\textrm{enclosed} }{ \\textrm{factors} }} f_a(x_1) = f_a(x_1)\\,$$\n", "since there are no enclosed variables. \n", "\n", "

\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The message from a half-edge is $1$ (one). You can verify this by imagining that a half-edge $x$ can be terminated by a node function $f(x)=1$ without affecting any inference issue." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Automating Bayesian Inference by Message Passing\n", "\n", "- The foregoing message update rules can be worked out in closed-form and put into tables (e.g., see Tables 1 through 6 in [Loeliger (2007)](./files/Loeliger-2007-The-factor-graph-approach-to-model-based-signal-processing.pdf) for many standard factors such as essential probability distributions and operations such as additions, fixed-gain multiplications and branching (equality nodes).\n", "\n", "- In the optional slides below, we have worked out a few of these update rules, eg, for the [equality node](#sp-for-equality-node), the [addition node](#sp-for-addition-node) and the [multiplication node](#sp-for-multiplication-node).\n", "\n", "- If the update rules for all node types in a graph have been tabulated, then inference by message passing comes down to executing a set of table-lookup operations, thus creating a completely **automatable Bayesian inference framework**. \n", "\n", "- In our research lab [BIASlab](http://biaslab.org) (FLUX 7.060), we are developing [ForneyLab](https://github.com/biaslab/ForneyLab.jl), which is a (Julia) toolbox for automating Bayesian inference by message passing in a factor graph.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example: Bayesian Linear Regression by Message Passing\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume we want to estimate some function $f: \\mathbb{R}^D \\rightarrow \\mathbb{R}$ from a given data set $D = \\{(x_1,y_1), \\ldots, (x_N,y_N)\\}$, with model assumption $y_i = f(x_i) + \\epsilon_i$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "#### model specification \n", "\n", "- We will assume a linear model with white Gaussian noise and a Gaussian prior on the coefficients $w$:\n", "$$\\begin{align*}\n", " y_i &= w^T x_i + \\epsilon_i \\\\\n", " \\epsilon_i &\\sim \\mathcal{N}(0, \\sigma^2) \\\\ \n", " w &\\sim \\mathcal{N}(0,\\Sigma)\n", "\\end{align*}$$\n", "or equivalently\n", "$$\\begin{align*}\n", "p(w,\\epsilon,D) &= \\overbrace{p(w)}^{\\text{weight prior}} \\prod_{i=1}^N \\overbrace{p(y_i\\,|\\,x_i,w,\\epsilon_i)}^{\\text{regression model}} \\overbrace{p(\\epsilon_i)}^{\\text{noise model}} \\\\\n", " &= \\mathcal{N}(w\\,|\\,0,\\Sigma) \\prod_{i=1}^N \\delta(y_i - w^T x_i - \\epsilon_i) \\mathcal{N}(\\epsilon_i\\,|\\,0,\\sigma^2) \n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "#### Inference (parameter estimation)\n", "\n", "- We are interested in inferring the posterior $p(w|D)$. We will execute inference by message passing on the FFG for the model." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- The left figure shows the factor graph for this model. \n", "- The right figure shows the message passing scheme. \n", " \n", "

\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### CODE EXAMPLE\n", "\n", "Let's solve this problem by message passing-based inference with Julia's FFG toolbox [ForneyLab](http://forneylab.org)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "using Pkg;Pkg.activate(\"probprog/workspace/\");Pkg.instantiate();\n", "IJulia.clear_output();" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGwCAYAAABLvHTgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df3RU9Z3/8dckkAToZGqAZJJDxCBuNUaoJGJBRCsEQzXU9azWapSeth6hiKZUtyLdjXGVtNqK3aXGxfVH1yzFPVWE9GtzTLcKRbDBxBRDeqzSVCLMbL4SnAQkCUzu9w9O5suYhMzvO3Pn+Thnztm5c/OZN3POel/9/LQZhmEIAADAYlLMLgAAACAaCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSxpldgFkGBwd1+PBh2e122Ww2s8sBAAABMAxDvb29ysvLU0rK2ftqkjbkHD58WPn5+WaXAQAAQtDZ2alp06ad9Z6kDTl2u13S6R8pMzPT5GoAAEAgenp6lJ+f73uOn03ShpyhIarMzExCDgAACSaQqSZMPAYAAJZEyAEAAJZEyAEAAJZEyAEAAJZEyAEAAJZEyAEAAJZEyAEAAJZEyAEAAJZEyAEAAJaUtDseAwCAyPMOGmrq6FZXb5+y7RmaW5Cl1BRzDsIm5AAAgIhoaHOpur5dLk+f71quI0NV5YUqK8qNeT0MVwEAgLA1tLm0sq7FL+BIktvTp5V1LWpoc8W8JkIOAAAIi3fQUHV9u4wRPhu6Vl3fLu/gSHdEDyEHAACEpamje1gPzpkMSS5Pn5o6umNXlAg5AAAgTF29owecUO6LFEIOAAAIS7Y9I6L3RQohBwAAhGVuQZZyHRkabaG4TadXWc0tyIplWYQcAAAQntQUm6rKCyVpWNAZel9VXhjz/XIIOQAAIGxlRbmqrZgjp8N/SMrpyFBtxRxT9slhM0AAABARZUW5Ki10suMxAACwntQUm+adP9nsMiQxXAUAACyKkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACwp7kNOTU2NbDabKisrfdf6+/u1evVqTZkyRZMmTdKyZcv08ccfm1glAACIN3Edcvbu3atNmzZp1qxZftcrKyu1detWbdmyRbt27dKxY8d0/fXXy+v1mlQpAACIN3Ebco4dO6bbbrtNzzzzjM455xzfdY/Ho2effVY/+9nPtHjxYl166aWqq6vTe++9p9/97nejttff36+enh6/FwAAycg7aGjPgSPa1npIew4ckXfQMLukqIjbkLNq1Spdd911Wrx4sd/15uZmnTx5UkuWLPFdy8vLU1FRkXbv3j1qezU1NXI4HL5Xfn5+1GoHACBeNbS5tOAnv9c3n3lb925p1TefeVsLfvJ7NbS5zC4t4uIy5GzZskUtLS2qqakZ9pnb7VZaWppf744k5eTkyO12j9rm2rVr5fF4fK/Ozs6I1w0AQDxraHNpZV2LXJ4+v+tuT59W1rVYLuiMM7uAz+vs7NS9996r119/XRkZGQH/nWEYstlso36enp6u9PT0SJQIAEDC8Q4aqq5v10gDU4Ykm6Tq+naVFjqVmjL68zSRxF1PTnNzs7q6ulRcXKxx48Zp3Lhx2rFjh/71X/9V48aNU05OjgYGBnT06FG/v+vq6lJOTo5JVQMAEN+aOrqH9eCcyZDk8vSpqaM7dkVFWdyFnEWLFum9995Ta2ur71VSUqLbbrvN93+PHz9ejY2Nvr9xuVxqa2vT/PnzTawcAID41dU7esAJ5b5EEHfDVXa7XUVFRX7XJk2apMmTJ/uuf+c739EPfvADTZ48WVlZWbrvvvt0ySWXDJukDAAATsu2BzYFJND7EkHchZxAbNiwQePGjdPNN9+sEydOaNGiRXrhhReUmppqdmkAAMSluQVZynVkyO3pG3Fejk2S05GhuQVZsS4tamyGYVhzcfwYenp65HA45PF4lJmZaXY5AABE3dDqKkl+QWdomnFtxRyVFeXGvK5gBPP8jrs5OQAAIDrKinJVWzFHTof/kJTTkZEQASdYCTlcBQAAQlNWlKvSQqeaOrrV1dunbPvpISqrLBs/EyEHAIAkk5pi07zzJ5tdRtQxXAUAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACyJkAMAACxpnNkFAACQSLyDhpo6utXV26dse4bmFmQpNcVmdlkYASEHAIAANbS5VF3fLpenz3ct15GhqvJClRXlmlgZRsJwFQAAAWhoc2llXYtfwJEkt6dPK+ta1NDmMqkyjIaQAwDAGLyDhqrr22WM8NnQter6dnkHR7oDZiHkAAAwhqaO7mE9OGcyJLk8fWrq6I5dURgTIQcAgDF09Y4ecEK5D7FByAEAYAzZ9oyI3ofYIOQAADCGuQVZynVkaLSF4jadXmU1tyArlmVhDIQcAADGkJpiU1V5oSQNCzpD76vKC9kvJ84QcgAACEBZUa5qK+bI6fAfknI6MlRbMYd9cuIQmwECABCgsqJclRY62fE4QRByAAAIQmqKTfPOn2x2GQgAw1UAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkAAMCS2PEYAGAJ3kGD4xbgh5ADAEh4DW0uVde3y+Xp813LdWSoqryQgzOTGMNVAICE5B00tOfAEf1L/X6tqGvxCziS5Pb0aWVdixraXCZVCLPRkwMASDgj9dx8niHJJqm6vl2lhU6GrpIQPTkAgITS0ObSyhF6bkZiSHJ5+tTU0R39whB36MkBACQM76Ch6vp2GUH+3Vsf/l8mJCchQg4AIGE0dXQH1IPzeRvfOOD7v5mQnDwYrgIAJIyu3uADzucxITl5EHIAAAkj254RdhtDQ13V9e3yDgY78IVEQsgBACSMuQVZynVkKNwZNUxITg6EHABAwkhNsamqvFCShgWdofffvuI83f3VmQG1F4nhL8QvQg4AIKGUFeWqtmKOnA7/oSunI0NPV8zRP5dfrCtmTgmorUgMfyF+sboKAJBwyopyVVroHPWsqqFhLbenb8Tl5jadDkVzC7JiWjdii5ADAEhIqSk2zTt/8qifVZUXamVdi2ySX9AZGtaqKi9kvxyLY7gKAGBJZxvWqq2Ywz45SYCeHACAZY01rAVrI+QAACztbMNasDaGqwAAgCURcgAAgCXFZcipra3VrFmzlJmZqczMTM2bN0+//e1vfZ/39/dr9erVmjJliiZNmqRly5bp448/NrFiAAAQb+Iy5EybNk0//vGP9c477+idd97RNddco69//evav3+/JKmyslJbt27Vli1btGvXLh07dkzXX3+9vF6vyZUDAIB4YTMMIyFOJ8vKytLjjz+uf/iHf9DUqVP14osv6hvf+IYk6fDhw8rPz9drr72ma6+9dsS/7+/vV39/v+99T0+P8vPz5fF4lJmZGZN/AwAACE9PT48cDkdAz++47Mk5k9fr1ZYtW3T8+HHNmzdPzc3NOnnypJYsWeK7Jy8vT0VFRdq9e/eo7dTU1MjhcPhe+fn5sSgfAACYJG5DznvvvacvfOELSk9P14oVK7R161YVFhbK7XYrLS1N55xzjt/9OTk5crvdo7a3du1aeTwe36uzszPa/wQAAGCiuN0n50tf+pJaW1v16aef6uWXX9by5cu1Y8eOUe83DEM22+ibO6Wnpys9PT0apQIAgDgUdk/O4sWL/VY+DQl3EnBaWppmzpypkpIS1dTUaPbs2fr5z38up9OpgYEBHT161O/+rq4u5eTkhPWdAADAOsIOOe+8847OO+88SVJHR4fv+rPPPqvbb7893OZ9DMNQf3+/iouLNX78eDU2Nvo+c7lcamtr0/z58yP2fQAAILGFPVw1MDAgu90uSZo9e7ZaW1s1Y8YMzZ8/Xw899FBIbT744INaunSp8vPz1dvbqy1btujNN99UQ0ODHA6HvvOd7+gHP/iBJk+erKysLN1333265JJLtHjx4nD/OQAAwCLCDjkzZ87UH//4R9ntdh0/flyffvqpJMlut6u7uzukNv/3f/9Xt99+u1wulxwOh2bNmqWGhgaVlpZKkjZs2KBx48bp5ptv1okTJ7Ro0SK98MILSk1NDfefAwAALCLsfXKefvpprV27VtOnT1dKSormzp2rp59+WnV1dVq3bp0++uijSNUaUcGsswcAAPEhmOd32D05K1as0NSpU/XBBx/ozjvv1C233KIZM2bI5XLp7rvvDrd5AACAkER8x+NTp05p69atGhgY0C233BK3Q0j05AAAkHhi2pMzrMFx43TTTTdFulkAAICgxO2OxwAAAOEg5AAAAEuKaMh56623/E76BgAAMEtEQ87SpUt16NChSDYJAAAQkoiGnAgv1AIAAAgZc3IAAIAlhbWE/D//8z/93p86dUqvvPKKsrOzfdfuuOOOcL4CAAAgJGGFnOeff97v/cmTJ/XrX/9aEyZMkCTZbDZCDgAAMEVYIeeNN97we2+327V582bNmDEjrKIAAADCxZwcAABgSYQcAABgSRENOQ8++KCysrIi2SQAAEBIIn4KeaLgFHIAABJPMM9vhqsAAIAlEXIAAIAlEXIAAIAlhbVPzpCTJ0/K7Xbrs88+09SpU5l8DAAATBdyT86xY8f07//+77r66qvlcDh03nnnqbCwUFOnTtX06dN15513au/evZGsFQAAIGAhhZwNGzbovPPO0zPPPKNrrrlGr7zyilpbW/X+++9rz549qqqq0qlTp1RaWqqysjJ98MEHka4bAADgrEJaQn7TTTfpn//5n3XJJZec9b7+/n49++yzSktL03e/+92Qi4wGlpADAJB4gnl+s08OIQcAgITBPjkAACDphbS6avv27UH/TWlpqSZMmBDK1wEAAAQtpJBzww03BHW/zWbTBx98oBkzZoTydQAAAEELebjK7XZrcHAwoNfEiRMjWTMAAMCYQgo5y5cvD2roqaKigsm9AAAgplhdxeoqAAASRlRXV504cUKHDh0adn3//v3BNgUAABA1QYWcX//61/q7v/s7fe1rX9OsWbP0xz/+0ffZ7bffHvHiAADxyTtoaM+BI9rWekh7DhyRdzApBwUQ54JaXfXII4+opaVFU6dO1TvvvKPly5dr3bp1uvXWW5Wko14AkHQa2lyqrm+Xy9Pnu5bryFBVeaHKinJNrAzwF1TIOXnypKZOnSpJKikp0c6dO3XjjTfqww8/lM1mi0qBAID40dDm0sq6Fn3+f9a6PX1aWdei2oo5BB3EjaCGq7Kzs7Vv3z7f+8mTJ6uxsVF//vOf/a4DAKzHO2iour59WMCR5LtWXd/O0BXiRlAh58UXX1R2drbftbS0NP3qV7/Sjh07IloYACC+NHV0+w1RfZ4hyeXpU1NHd+yKAs4iqJAzbdo0OZ1Ov2tPPvmkXC6XrrjiiogWBgCIL129owecUO4Doi3sAzrXrFmjK664Qh9//LHf9YGBAe3duzfc5gEAcSLbnhHR+4Boi8gp5EuXLtXChQv9gs7Ro0f1la98JRLNAwDiwNyCLOU6MjTaMhObTq+ymluQFcuygFGFHXJsNpuqqqq0fPnyYUGHZeUAYB2pKTZVlRdK0rCgM/S+qrxQqSmstkV8iEhPjqQRgw7LygHAWsqKclVbMUdOh/+QlNORwfJxxJ2g9skZyZm9NVVVVZKkhQsX6qWXXgq3aQBAHCorylVpoVNNHd3q6u1Ttv30EBU9OIg3YYecRx99VJMmTfK9Hwo61113XbhNAwDiVGqKTfPOn2x2GcBZhR1y1q5dO+xaVVWVUlNT9dOf/jTc5gEAAEJiM5J0dnAwR7UDAID4EMzzO2ITjwEAAOIJIQcAAFgSIQcAAFhSVENOSkqKrrnmGjU3N0fzawAgIXgHDe05cETbWg9pz4EjvtO6R7sOIDxhr646m+eee04fffSR7rnnHr311lvR/CoAiGsNbS5V17f7neKd68jQstm52v4n17DrVeWFAW+s5x002LMGGAGrq1hdBSDKGtpcWlnXokD/YzsUTwLZQXi08BRMSAISSUxXVz355JM6fPhwuM0AgCV5Bw1V17cHHHAk+e6trm8/69DVUHg6M+BIktvTp5V1LWpocwVfMGAhYYecNWvW6Morr/Q7mFOSBgYGtHfv3nCbB4CE1tTRPSyEBMKQ5PL0qamje8TPzxaeAg1JgNVFZOJxWVnZsBPIjx49qq985SuRaB4AElZXb/ABJ5C/Hys8jRWSgGQQ9sRjm82mqqoqZWdna+HChdq5c6emTZsmyf/wTgBIRtn2jLFvCuHvAw1P4YYsIJFFbHXVmSeQ79y5U+PHj5fNxux+AMltbkGWch0Zcnv6gpqXY5PkdJxeKTWSQMNTuCELSGRhD1ed2VtTVVWl5cuXa+HChTp48GC4TQNAwktNsamqvFDS/181NZah+6rKC0ddCj4UnkZr06bTq6xGC0lAMgg75Dz66KOaNGmS7/1Q0LnuuuvCbRoALKGsKFe1FXPkdPj3quQ6MnTXwgLlfu6605Ex5vLxs4WnQEISkAyitk/OI488op/+9Kf69NNPo9F82NgnB0CsjbZpXzib+bFPDpJNMM9vNgMk5ABIcOx4jGQSzPM7qsc6AACiLzXFpnnnTza7DCDucAo5AACwJHpyACAEDBEB8S+qISclJUVXX321Hn/8cRUXF0fzqwAgZpjsCySGqA5XPffcc7rqqqt0zz33RPNrACBmOBQTSBxRDTnf+ta3VFVVpbfeeiuov6upqdFll10mu92u7Oxs3XDDDXr//ff97unv79fq1as1ZcoUTZo0ScuWLRt2SCgARBKHYgKJJS4nHu/YsUOrVq3S22+/rcbGRp06dUpLlizR8ePHffdUVlZq69at2rJli3bt2qVjx47p+uuvl9frNbFyAFbGoZhAYglrTs7HH3+s2tpa7d69W263WzabTTk5OZo/f75WrFih/Pz8kNptaGjwe//8888rOztbzc3NWrhwoTwej5599lm9+OKLWrx4sSSprq5O+fn5+t3vfqdrr702nH8WAIyIQzGBxBJyT86uXbt00UUXaevWrZo9e7buuOMOVVRUaPbs2Xr11Vd18cUXBz1MNRqPxyNJyso6fQZLc3OzTp48qSVLlvjuycvLU1FRkXbv3j1iG/39/erp6fF7AUAwOBQTSCwh9+R8//vf13e/+11t2LBh1M8rKyu1d+/ekIuTTh8AumbNGi1YsEBFRUWSJLfbrbS0NJ1zzjl+9+bk5Mjtdo/YTk1Njaqrq8OqBUByG+tE8bFODgcQWyH35LS1tWnFihWjfn7XXXepra0t1OZ97r77bu3bt0+/+tWvxrzXMAzZbCPvU7F27Vp5PB7fq7OzM+zaACQXDsUEEkvIISc3N3fUoSFJ2rNnj3Jzw9svYvXq1dq+fbveeOMNTZs2zXfd6XRqYGBAR48e9bu/q6tLOTk5I7aVnp6uzMxMvxcABGu0E8UDOTkcQGyFPFx13333acWKFWpublZpaalycnJks9nkdrvV2Nio//iP/9CTTz4ZUtuGYWj16tXaunWr3nzzTRUUFPh9XlxcrPHjx6uxsVE333yzJMnlcqmtrU2PPfZYqP8kAAhIWVGuSgud7HgMxLmwTiF/6aWXtGHDBjU3N/uWbqempqq4uFhr1qzxBZBgfe9739PmzZu1bds2felLX/JddzgcmjBhgiRp5cqV+s1vfqMXXnhBWVlZuu+++3TkyBE1NzcrNTV1zO/gFHIAABJPMM/vsELOkJMnT+qTTz6RJE2ZMkXjx48Pq73R5tU8//zz+ta3viVJ6uvr0/3336/NmzfrxIkTWrRokZ566qmAl60TcgAASDwxDzmJiJADAEDiCeb5HbUdjxcvXqwZM2ZEq3kASAjeQUN7DhzRttZD2nPgCEc+ADEUtVPI//7v/943hAUAyYjTygFzMVzFcBWAKBg6rfzz/4EdmnHIcnMgNME8v6PWkwMAseYdNKK+rDuQ7xjrtHKbTp9WXlroZNk5EEVRCzmdnZ2qqqrSc889F62vAACfWAwNBfodwZxWPu/8yRGpDcBwUZt43N3drV/+8pfRah4AfIaGhj4fLNyePq2sa1FDmytq3+Ea4Ts4rRyIDyH35Gzfvv2sn//1r38NtWkACFgshobO9h1D33Pmd3BaORAfQg45N9xwg2w2m842b3m0Tf0AIFJiMTQ01nfoc9/BaeVAfAjrgM6XX35Zg4ODI75aWloiWScAjCgWQ0Nuz4mg7uO0ciA+hBxyiouLzxpkxurlAYBIiMXQUPfxgaDv47RywHwhD1fdf//9On78+Kifz5w5U2+88UaozQNAQGIxNJT1hfSQ7uO0csBcIYecK6+88qyfT5o0SVdddVWozQNAQIaGhlbWtcgm+QWdSA0NHTwy+v+gO5Mzc3hvUWqKjWXigElCGq7at2+fBgcHA75///79OnXqVChfBQBjOtvQ0C9unSPHhLSQz47yDhr6VdPBMe/LZSIxEHdC6sm59NJL5Xa7NXXq1IDunzdvnlpbWzmwE0DUjDQ0dPR4v/7l/4S3QWBTR7fcPf1j3nfLZecyDAXEmZBCjmEY+qd/+idNnDgxoPsHBgKbtAcA4ThzaKihzaVVm98dNk9naIPAQCf/Broq67wpgf33EEDshBRyFi5cqPfffz/g++fNm6cJEyaE8lUAELRIbhDIxn5A4gop5Lz55psRLgMAIieSGwSysR+QuKJ2dhUAmCWSGwSysR+QuEIKOQcPjr3S4EyHDh0K5WsAICSRHmJiYz8gMYU0XHXZZZdp2bJluvPOOzV37twR7/F4PPrv//5v/fznP9ddd92l1atXh1UoAAQqGkNMbOwHJJ6QQs7Xv/512e12lZWVafz48SopKVFeXp4yMjJ09OhRtbe3a//+/SopKdHjjz+upUuXRrpuABhVtDYIZGM/ILHYjBAOmEpLS1NnZ6cyMzOVk5Ojm2++WUeOHNGJEyc0ZcoUXXrppbr22mtVVFQUjZojoqenRw6HQx6PR5mZmWaXAyAKGtpcqq4Pb58cAPElmOd3SCGnoKBAtbW1KisrU0pKitxut7Kzs0Mu2AyEHCA5eAcNhpgACwnm+R3ScNV9992nZcuWqaSkRDabTf/1X/+lBQsWqKioiP1wAMQVhpiA5BVST450+jyqbdu26Uc/+pFmzJihv/3tb7LZbJo5c6Zmz56tL3/5y5o9e3bczsehJwcAgMQT9eGqM82cOVNvv/22Jk2apH379qm1tdX3amtrU29vbzjNRw0hBwCAxBPTkHM2hmHIZovPsW9CDpB4mF8DIOpzcgIVrwEHQOJhpRSAYHGsA4C419Dm0sq6lmHnUQ2dKN7Q5jKpMgDxjJADIK6NdaK4dPpEce9g1EbeASQoQg6AoHgHDe05cETbWg9pz4EjUQ8XwZwoDgBniuqcHADWYsa8mEieKA4gudCTAyAgZs2LifSJ4gCSByEHwJjMnBczdKL4aGs1bTrdmxTMieIAkgMhB8CYzJwXM3SiuKRhQSecE8UBWB8hB8CYzJ4XU1aUq9qKOXI6/IeknI4M1VbMYZ8cACNi4jGAMcXDvJiyolyVFjrZ8RhAwAg5AMY0NC/G7ekbcV6OTad7VaI9L4YTxQEEg+EqAGNiXgyARETIARAQ5sUASDQMVwEIGPNiACQSQg6AoDAvBkCiYLgKAABYEiEHAABYEiEHAABYEnNygATmHTSiOgk42u0DQDQRcoAE1dDmUnV9u9+ZUrmODFWVF0ZkOXe02weAaGO4CkhADW0uraxrGXZoptvTp5V1LWpoc8V1+wAQC4QcIMF4Bw1V17ePeLzC0LXq+nZ5B0e6I3btewcN7TlwRNtaD2nPgSMh1wMAoWK4CkgwTR3dw3pYzmRIcnn61NTRHdJ+NpFon6EuAPGAnhwgwXT1jh5AQrkv0u0z1AUgXhBygASTbc8Y+6Yg7otk+9EeSgOAYBBygAQztyBLuY6MYaeBD7Hp9NDQ3IKsmLcfzFAXAEQbIQdIMKkpNlWVF0rSsCAy9L6qvDDk/WzCaT/aQ2kAEAxCDpCAyopyVVsxR06H/5CR05Gh2oo5YU/uDbX9aA+lAUAwWF0FJKiyolyVFjqjtiNxKO0PDXW5PX0jzsux6XRQCnUoDQCCQcgBElhqii2kZeLRan9oqGtlXYtskl/QicRQGgAEg+EqABEV7aE0AAgUPTkAIi7aQ2kAEAhCDoCoiPZQGgCMheEqAABgSYQcAABgSYQcAABgSYQcAABgSYQcAABgSXEZcnbu3Kny8nLl5eXJZrPp1Vdf9fvcMAw99NBDysvL04QJE3T11Vdr//79JlULAADiUVyGnOPHj2v27NnauHHjiJ8/9thjeuKJJ7Rx40bt3btXTqdTpaWl6u3tjXGlAAAgXsXlPjlLly7V0qVLR/zMMAw9+eSTWrdunW688UZJ0i9/+Uvl5ORo8+bNuuuuu2JZKgAAiFNx2ZNzNh0dHXK73VqyZInvWnp6uq666irt3r171L/r7+9XT0+P3wuINO+goT0Hjmhb6yHtOXBE3sGRjqkEAMRCXPbknI3b7ZYk5eTk+F3PycnRRx99NOrf1dTUqLq6Oqq1Ibk1tLlUXd8ul6fPdy3XkaGq8kLOawIAEyRcT84Qm83/DBzDMIZdO9PatWvl8Xh8r87OzmiXiCTS0ObSyroWv4AjSW5Pn1bWtaihzWVSZQCQvBIu5DidTkn/v0dnSFdX17DenTOlp6crMzPT7wVEgnfQUHV9u0YamBq6Vl3fztAVAMRYwoWcgoICOZ1ONTY2+q4NDAxox44dmj9/vomVIVk1dQsiXJUAAA7XSURBVHQP68E5kyHJ5elTU0d37IoCAMTnnJxjx47pww8/9L3v6OhQa2ursrKydO6556qyslLr16/XBRdcoAsuuEDr16/XxIkTdeutt5pYNZJVV+/oASeU+wAAkRGXIeedd97RV7/6Vd/7NWvWSJKWL1+uF154Qf/4j/+oEydO6Hvf+56OHj2qyy+/XK+//rrsdrtZJSOJZdszInofACAybIZhJOVEgZ6eHjkcDnk8HubnICzeQUMLfvJ7uT19I87LsUlyOjK064fXKDVl9MnxAICxBfP8Trg5OUC8SU2xqaq8UNLpQHOmofdV5YUEHACIMUIOEAFlRbmqrZgjp8N/SMrpyFBtxRz2yQEAE8TlnBwg3nkHDTV1dKurt0/Z9gzNLchSWVGuSgudw67TgwMA5iDkwJJGCiGRChtj7Ww87/zJEfkeAEB4CDmwnGgerzC0s/HnJxgP7WzM0BQAxA/m5MBSonm8AjsbA0BiIeTAMqIdQtjZGAASCyEHlhHtEMLOxgCQWAg5sIxohxB2NgaAxELIgWVEO4TMLchSriNj2IZ/Q2w6PcF5bkFWSO0DACKLkAPLiHYIYWdjAEgshBzEjHfQ0J4DR7St9ZD2HDgS8VVIsQgh7GwMAImDAzo5oDMmorl3jRnfFc3NBgEAowvm+U3IIeRE3Wgb6A1Fgmj0gBBCAMCagnl+s+MxomqsvWtsOr13TWmhM6IhJDXFxvEKAJDkmJODqGIDPQCAWQg5iCo20AMAmIWQg6hiAz0AgFkIOYgqNtADAJiFkIOoCmTvmlsuO1e/2Xc4KnvnAACSF0vIWUIeEyPtXXPOxPEyJH362UnftWjtnQMAsAb2yQkAISf2zty75m+fHNeG330w7J5o7p0DAEh8wTy/Ga5CzAztXXP9rDxt2ds54j1Dibu6vp2hKwBAWAg5iDn2zgEAxAIhBzHH3jkAgFgg5CDm2DsHABALhBzEHHvnAABigZCDmAtk75yq8kJODQcAhIWQA1OUFeWqtmKOnA7/ISmnI4Pl4wCAiBhndgFIXmVFuSotdPr2zsm2nx6iogcHABAJhByYamjvHAAAIo3hKgAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmEHAAAYEmcQh5h3kFDTR3d6urtU7Y9Q3MLspSaYjO7LAAAkg4hJ4Ia2lyqrm+Xy9Pnu5bryFBVeaHKinJNrAwAgOTDcFWENLS5tLKuxS/gSJLb06eVdS1qaHOZVBkAAMmJkBMB3kFD1fXtMkb4bOhadX27vIMj3QEAAKKBkBMBTR3dw3pwzmRIcnn61NTRHbuiAABIcoScCOjqHT3ghHIfAAAIHxOPIyDbnhHWfazIAgAg8gg5ETC3IEu5jgy5PX0jzsuxSXI6ToeXz2NFFgAA0cFwVQSkpthUVV4o6XSgOdPQ+6rywmG9M6zIAgAgegg5EVJWlKvaijlyOvyHpJyODNVWzBnWK8OKLAAAoovhqggqK8pVaaEzoPk1wazImnf+5ChWDQCANRFyIiw1xRZQKGFFFgAA0cVwlUnCXZEFAADOjpBjkqEVWaMtFLfp9CqrkVZkAQCAsRFyTHK2FVlDRlqRBQAAAkPIMdHQiizHxPHDPhvpGgAACBwhJw58+tnJYdc8n51krxwAAMJAyDHR0F45I2GvHAAAwkPIMRGnlwMAED2EHBOxVw4AANFDyDERe+UAABA9CR1ynnrqKRUUFCgjI0PFxcX6wx/+YHZJQWGvHAAAoidhQ85LL72kyspKrVu3Tu+++66uvPJKLV26VAcPHjS7tICFeno5AAAYm80wjIRcunP55Zdrzpw5qq2t9V276KKLdMMNN6impmbMv+/p6ZHD4ZDH41FmZmY0Sx1TQ5tL1fXtfpOQcx0ZqiovHHZ6OQAAySyY53dCHtA5MDCg5uZmPfDAA37XlyxZot27d4/4N/39/erv7/e97+npiWqNwQjm9HIAABCYhAw5n3zyibxer3Jycvyu5+TkyO12j/g3NTU1qq6ujkV5IQn09HIAABCYhJ2TI0k2m39Ph2EYw64NWbt2rTwej+/V2dkZixIBAIBJErInZ8qUKUpNTR3Wa9PV1TWsd2dIenq60tPTY1EeAACIAwnZk5OWlqbi4mI1Njb6XW9sbNT8+fNNqgoAAMSThOzJkaQ1a9bo9ttvV0lJiebNm6dNmzbp4MGDWrFihdmlAQCAOJCwIecb3/iGjhw5oocfflgul0tFRUV67bXXNH36dLNLAwAAcSBh98kJVzztkwMAAAITzPM7IefkAAAAjIWQAwAALImQAwAALClhJx6Ha2gqUjwd7wAAAM5u6LkdyJTipA05vb29kqT8/HyTKwEAAMHq7e2Vw+E46z1Ju7pqcHBQhw8flt1uH/UoiLPp6elRfn6+Ojs7WZ0VQ/zu5uG3Nwe/uzn43c0RyO9uGIZ6e3uVl5enlJSzz7pJ2p6clJQUTZs2Lex2MjMz+X8AE/C7m4ff3hz87ubgdzfHWL/7WD04Q5h4DAAALImQAwAALCn1oYceesjsIhJVamqqrr76ao0bl7SjfqbgdzcPv705+N3Nwe9ujkj+7kk78RgAAFgbw1UAAMCSCDkAAMCSCDkAAMCSCDkAAMCSCDkheuqpp1RQUKCMjAwVFxfrD3/4g9klWVpNTY0uu+wy2e12ZWdn64YbbtD7779vdllJp6amRjabTZWVlWaXYnmHDh1SRUWFJk+erIkTJ+rLX/6ympubzS7L8k6dOqUf/ehHKigo0IQJEzRjxgw9/PDDGhwcNLs0S9m5c6fKy8uVl5cnm82mV1991e9zwzD00EMPKS8vTxMmTNDVV1+t/fv3B/09hJwQvPTSS6qsrNS6dev07rvv6sorr9TSpUt18OBBs0uzrB07dmjVqlV6++231djYqFOnTmnJkiU6fvy42aUljb1792rTpk2aNWuW2aVY3tGjR3XFFVdo/Pjx+u1vf6v29nb97Gc/0xe/+EWzS7O8n/zkJ3r66ae1ceNG/fnPf9Zjjz2mxx9/XP/2b/9mdmmWcvz4cc2ePVsbN24c8fPHHntMTzzxhDZu3Ki9e/fK6XSqtLTUd+5kwAwEbe7cucaKFSv8rl144YXGAw88YFJFyaerq8uQZOzYscPsUpJCb2+vccEFFxiNjY3GVVddZdx7771ml2RpP/zhD40FCxaYXUZSuu6664xvf/vbftduvPFGo6KiwqSKrE+SsXXrVt/7wcFBw+l0Gj/+8Y991/r6+gyHw2E8/fTTQbVNT06QBgYG1NzcrCVLlvhdX7JkiXbv3m1SVcnH4/FIkrKyskyuJDmsWrVK1113nRYvXmx2KUlh+/btKikp0U033aTs7GxdeumleuaZZ8wuKyksWLBA//M//6O//OUvkqQ//elP2rVrl772ta+ZXFny6OjokNvt9nvOpqen66qrrgr6Ocs2jkH65JNP5PV6lZOT43c9JydHbrfbpKqSi2EYWrNmjRYsWKCioiKzy7G8LVu2qKWlRXv37jW7lKTx17/+VbW1tVqzZo0efPBBNTU16Z577lF6erruuOMOs8uztB/+8IfyeDy68MILlZqaKq/Xq0cffVTf/OY3zS4taQw9S0d6zn700UdBtUXICZHNZvN7bxjGsGuIjrvvvlv79u3Trl27zC7F8jo7O3Xvvffq9ddfV0ZGhtnlJI3BwUGVlJRo/fr1kqRLL71U+/fvV21tLSEnyl566SXV1dVp8+bNuvjii9Xa2qrKykrl5eVp+fLlZpeXVCLxnCXkBGnKlClKTU0d1mvT1dU1LHUi8lavXq3t27dr586dmjZtmtnlWF5zc7O6urpUXFzsu+b1erVz505t3LhR/f39Sk1NNbFCa8rNzVVhYaHftYsuukgvv/yySRUlj/vvv18PPPCAbrnlFknSJZdcoo8++kg1NTWEnBhxOp2STvfo5Obm+q6H8pxlTk6Q0tLSVFxcrMbGRr/rjY2Nmj9/vklVWZ9hGLr77rv1yiuv6Pe//70KCgrMLikpLFq0SO+9955aW1t9r5KSEt12221qbW0l4ETJFVdcMWyLhL/85S+aPn26SRUlj88++0wpKf6PxtTUVJaQx1BBQYGcTqffc3ZgYEA7duwI+jlLT04I1qxZo9tvv10lJSWaN2+eNm3apIMHD2rFihVml2ZZq1at0ubNm7Vt2zbZ7XZfT5rD4dCECRNMrs667Hb7sHlPkyZN0uTJk5kPFUXf//73NX/+fK1fv14333yzmpqatGnTJm3atMns0iyvvLxcjz76qM4991xdfPHFevfdd/XEE0/o29/+ttmlWcqxY8f04Ycf+t53dHSotbVVWVlZOvfcc1VZWan169frggsu0AUXXKD169dr4sSJuvXWW4P7oois/0pCv/jFL4zp06cbaWlpxpw5c1jKHGWSRnw9//zzZpeWdFhCHhv19fVGUVGRkZ6eblx44YXGpk2bzC4pKfT09Bj33nuvce655xoZGRnGjBkzjHXr1hn9/f1ml2Ypb7zxxoj/TV++fLlhGKeXkVdVVRlOp9NIT083Fi5caLz33ntBf4/NMAwjEqkMAAAgnjAnBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhBwAAWBIhB4BlrF+/XjabbdjriSeeMLs0ACbggE4AltHb26vjx4/73j/88MN67bXXtGvXLk2bNs3EygCYYZzZBQBApNjtdtntdklSdXW1XnvtNe3YsYOAAyQphqsAWE51dbWef/557dixQ9OnTze7HAAmIeQAsBQCDoAhhBwAlkHAAXAm5uQAsIRHHnlEGzdu1G9+8xulp6fL7XZLks455xylp6ebXB0AM7C6CkDCMwxDX/ziF9XT0zPss7fffluXX365CVUBMBshBwAAWBJzcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCURcgAAgCX9P5SOAQN4IiXNAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using PyPlot, ForneyLab, LinearAlgebra\n", "\n", "# Parameters\n", "Ξ£ = 1e5 * Diagonal(I,3) # Covariance matrix of prior on w\n", "Οƒ2 = 2.0 # Noise variance\n", "\n", "# Generate data set\n", "w = [1.0; 2.0; 0.25]\n", "N = 30\n", "z = 10.0*rand(N)\n", "x_train = [[1.0; z; z^2] for z in z] # Feature vector x = [1.0; z; z^2]\n", "f(x) = (w'*x)[1]\n", "y_train = map(f, x_train) + sqrt(Οƒ2)*randn(N) # y[i] = w' * x[i] + Ο΅\n", "scatter(z, y_train); xlabel(L\"z\"); ylabel(L\"f([1.0, z, z^2]) + \\epsilon\");" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "Now build the factor graph in ForneyLab, perform sum-product message passing and plot results (mean of posterior)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeWzU973/++dsnrHHHu87NraxDQZs43VsAwHCkgBJmqbZmiaNdNLt156e9qhVpbTSuT1HPYlOpbbRSU+rm+U0CWRpSoGEJCQhIQGCt8E2XjEY430F72OPZ/3eP/j5e+sAbRIIY8z7IY1UxvZ8Px8U1S8+y/utURRFQQghhBBiEdH6ewBCCCGEENeaBBwhhBBCLDoScIQQQgix6EjAEUIIIcSiIwFHCCGEEIuOBBwhhBBCLDoScIQQQgix6Oj9PQB/8fl89Pf3ExISgkaj8fdwhBBCCPEZKIrC1NQUCQkJaLVXXqe5aQNOf38/SUlJ/h6GEEIIIb6Anp4elixZcsWv37QBJyQkBLj4F2SxWPw8GiGEEEJ8FpOTkyQlJam/x6/kpg04c9tSFotFAo4QQghxg/lHx0vkkLEQQgghFh0JOEIIIYRYdCTgCCGEEGLRkYAjhBBCiEVHAo4QQgghFh0JOEIIIYRYdCTgCCGEEGLRkYAjhBBCiEVHAo4QQgghFh0JOEIIIYRYdCTgCCGEEGLRkYAjhBBCiEVHAo4QQgghrimn08nY2JhfxyABRwghhBDXjKIo1NTUUF5eTl9fn9/GIQFHCCGEENdMS0sLIyMjaLVaQkND/TYOCThCCCGEuCb6+vo4d+4cAGvWrCE4ONhvY5GAI4QQQoirNjk5SX19PQAZGRnEx8f7dTwScIQQQghxVVwuFzabDa/XS3R0NMuXL/f3kCTgCCGEEOKLUxSF2tpaZmZmCAoKoqCgAI1G4+9hScARQgghxBfX2trK+fPn0el0FBUVYTAYcLlcck1cCCGEEDemgYEBzp49C0Bubi4WiwVFUThx4gTHjh2jv7/fb2OTgCOEEEKIz21qaoq6ujoAli1bRmJiIgCnTp2ivr6e5uZmv25VScARQgghxOfidrvVQ8VRUVFkZWUBF6+JV1dX09/fT0JCAj6fz29jlIAjhBBCiM9s7lDx9PQ0gYGB6qHiyclJysvL6ejoICEhgaKiInVVxx8k4AghhBDiMzt9+jTDw8NotVqKiooICAjA7XZz/PhxWltbCQ4OJi8vT13V8RcJOEIIIYT4TPr7+2lrawMuVioODQ1FURRsNhsNDQ1oNBpyc3MXxFVxCThCCCGE+IcmJyc5efIkMP9QcWtrKydOnGBmZoasrCzKysrQ6XQ4HA5/DlcCjhBCCCH+vk9XKp7bfurv76e8vJzz58+TmprK+vXrCQ4Opqqqiv379zM6Ouq3Mev99mQhhBBCLHhXqlQ8NTXF0aNH6erqIi4ujnXr1hETE0NLSwsff/wxDoeD5cuXExER4ZdxywqOEEIIIa7o1KlTl1QqdrvdfPLJJ5w+fZqQkBBKSkpIT0+np6eHd999l5mZGTIyMvx60FgCjhBCCCEuq6+vj/b2dgDy8vLUSsU2m42TJ0+i0+nIz89nzZo1jIyM8PbbbzM2NkZCQgJbtmwhMDDQb2OXgCOEEEKIS0xMTKiHijMyMoiPjwcuruhUVlYyOzvLqlWrKCsrw+Fw8M4779DX10dERATp6ek0NjYyMTHht/FLwBFCCCHEPE6nE5vNhs/nIzY2luXLlwMXV3SOHDnC2NgYy5YtY8OGDWi1Wg4dOkRbWxvBwcEsX76c2dlZnE6nX29SScARQgghhMrn83HixAkcDodatE+j0TAxMcHhw4fp7+8nPj6eTZs2ERoaytGjR2loaCAgIIAVK1bg8XgAyM7OJi4uzm/zkIAjhBBCCFVjYyOjo6MYDAaKi4sxGAy4XC4+/vhjzp49S2hoKLfccgtLliyhuroam82GoihkZGQAoNFoSE9PJyUlxa/zkIAjhBBCCAA6Ojro7u5Go9GQn5+P2WzG5/NRXl5OfX09BoOB0tJSVq1aRXNzM8ePH2d2dpZly5ZhMpnQarUkJiayYsUKf09FAo4QQggh4MKFCzQ3NwOQlZVFTEwMAE1NTVRVVeHz+cjPz6e4uJjOzk4+/vhjJiYmSE5OJiQkBK1WS2RkJGvWrPF7mwaQgCOEEELc9Kanpzlx4gSKorBkyRKWLVsGQHd3N4cPH8Zut5OZmcmGDRsYGRnhww8/ZGhoiPj4eMLDw9HpdISEhFBUVIRWezFazJ3F8RcJOEIIIcRNzOPxYLPZcLvdhIWFkZubC8DY2Bjvv/8+IyMjLFmyhK1bt+JyuTh8+DDd3d1ERkYSGRlJQEAAJpMJq9WKwWAA4PDhw/zmN79hamrKb/OSgCOEEELcpBRFoa6ujqmpKUwmk7oCMzs7y/vvv09XVxfh4eFs3bqVoKAgPv74Y86cOYPFYiE6OpqgoCD0ej1Wq1Ut6lddXc2uXbuoq6vjww8/9NvcJOAIIYQQN6nTp08zODiIVqulqKgIk8mE1+vlyJEjnDp1isDAQDZu3MiSJUs4duwYTU1NGI1GYmJiCA4OVn/OYrEAF8/rPPvss2qrho0bN/ptbgsy4PT19fHwww8TGRlJUFAQa9asoaamRv26oij88pe/JCEhQf3LnzsYJYQQQoh/rL+/n7a2NgByc3MJCwsDwGazYbPZ0Gg0lJSUkJ2dTUVFBfX19Wg0GuLi4tRwk5ubS1RUFADt7e384Q9/YHJykuTkZH74wx+qn+kPCy7gjI2NsXbtWgwGAwcPHqSlpYXf/OY38/6Sfv3rX/Pb3/6W3//+99hsNuLi4ti6datf9/qEEEKIG8X4+Dh1dXUALFu2jCVLlgBw5swZjhw5gsvlIicnh7Vr11JXV0dtbS0Oh4O4uDjMZjMGg4GVK1eqP9fb28vTTz/NyMgIsbGx/PCHP1RvYfmL3q9Pv4z/+q//IikpiT/96U/qe39bLEhRFJ566il+8YtfcM899wDw4osvEhsbyyuvvMJ3v/vd6z1kIYQQ4oYxOzs7rw3DXMfvoaEh3nnnHex2O+np6WzdupW2tjZqamoYGxsjMTERs9lMYGAgaWlp6k2r8+fP8/TTTzMwMEBERAQ/+MEPSE5O9ucUgQW4gvPmm29SWFjIfffdR0xMDHl5eTz77LPq1zs6OhgcHGTbtm3qe0ajkQ0bNlBeXn7Fz3U6nUxOTs57CSGEEDcTr9eLzWZjdnaWkJAQ8vPz0Wg0TE1N8dZbb6krMDt27GB4eBibzaa2ZggMDCQ4OJiEhARWrlwJXGzI+fTTT9PV1YXFYuE73/mO2rfK3xZcwDl37hx//OMfycjI4L333uN73/se//Iv/8JLL70EwODgIACxsbHzfi42Nlb92uU8+eSThIaGqq+kpKQvbxJCCCHEAlRfX8/4+DgBAQEUFxej1+txu9289957dHZ2EhISwo4dO3C73VRWVtLR0UFsbCwBAQGEh4cTGRmp9qaanp7mf/7nfzhz5gyBgYE8+uij5OXl+XuKqgUXcOYqJT7xxBPk5eXx3e9+l29/+9v88Y9/nPd9n66SqCjK362c+PjjjzMxMaG+enp6vpTxCyGEEAtRW1sbfX19aDQaCgsLCQoKQlEUPv74YxobGzEajWzdupXQ0FAqKio4c+YMERER6PV6oqOjsVgs6jVyl8vFs88+q/7cgw8+yNq1a/09xXkWXMCJj49Xl77mZGVl0d3dDaB2Jv30as3w8PAlqzp/y2g0YrFY5r2EEEKIm8Hg4CCtra0A5OTkEBkZCUBtbS0VFRUAlJaWkpaWRnl5OadOnSI4OBi9Xk9cXBxBQUFqIT+3283//u//YrPZ0Ol03H333WzZsmVBtGf4Wwsu4Kxdu5bTp0/Pe+/MmTMsXboUgNTUVOLi4jh06JD6dZfLxZEjRygrK7uuYxVCCCEWusnJSWpra4GLv0PnDgB3dHRw6NAh3G43OTk5FBUVUVlZSUtLC3q9HqPRSEJCAiaTiZKSEkwmEz6fj5dffplPPvkERVG4/fbbueOOO9DpdP6c4mUtuFtU//qv/0pZWRlPPPEE999/P9XV1TzzzDM888wzwMWtqR//+Mc88cQTZGRkkJGRwRNPPEFQUBAPPfSQn0cvhBBCLBxOp5Pq6mq8Xi/R0dGsWrUKuNhY84033mB6epqUlBQ2b95MTU0NLS0teDweQkJCiI2NxWg0UlxcTHBwMIqisGfPHg4fPozH42HLli3ce++9anuGhWbBBZyioiL27dvH448/zn/8x3+QmprKU089xTe+8Q31e372s5/hcDj4/ve/z9jYGFarlffff5+QkBA/jlwIIYRYOHw+HzabDYfDgdlspqCgQD0cvG/fPkZGRoiKiuKuu+6ipaWFU6dOMTU1RVhYmFpot6CggIiICBRF4a233uLtt9/G5XKxbt06vv71r2MymS77bI/Hw/j4uFoE0B80iqIofnu6H01OThIaGsrExIScxxFCCLHo1NbW0tfXh8FgYP369ZjNZtxuN3v37qWpqQmz2cwDDzzA6Ogo9fX1dHd3ExERQXBwMNHR0eTm5qrbWYcOHeLVV19lZmaGwsJCvvOd71zxd6fH46GqqoqxsTEKCwvVs7PXymf9/b3gzuAIIYQQ4up8+saU2WxGURQ++ugjmpubMRgMbN++nZmZGU6dOkVXVxeRkZEYDAaio6NZsWKFGm6OHz/O66+/zvT0NNnZ2fzTP/3TFYOF1+ulurqa0dFR9Hq92oDTHyTgCCGEEItIf3+/emMqOztb3SY6ceIE5eXlaDQaNmzYgNlspqmpibNnzxIZGYmiKCxZsoTU1FQyMjKAi6tAL7/8MpOTk6xYsYLHHnuMiIiIyz53bktsZGQEvV5PSUkJoaGh12fSlyEBRwghhFgk/rbHVFpamnoD+ezZs7z33nt4vV7y8vJISUmhvr6eM2fOEBYWhtvtZunSpSQkJKgHkZubm3nhhRcYGRkhNTWVxx577IrbTXPh5vz58+h0OqxWq18bbYIEHCGEEGJRcDgcVFdXqz2m5mrKnT9/nn379uF0Olm2bBkFBQWcPHmSM2fOYDabcblcpKWlERsbq1Ypbmtr4/nnn2doaIikpCQee+yxK/aX8vl81NTUMDw8rIabK63yXE8ScIQQQogbnMfjobq6GqfTicViUXtMTU9P8/rrrzM5OUlMTAybN2+moaGBtrY29Ho9TqeT1NRUIiMjKSwsRKvV0tnZyXPPPUdfXx+xsbE8+uij6pbVpymKQm1tLYODg2i1WoqLi9Uigv624K6JCyGEEOKzUxSFuro6Jicn1bo1cz2m9u7dy+DgIMHBwezYsYOWlhba29vxeDwoikJKSgqRkZFYrVb0ej19fX08++yzdHV1ERUVxSOPPEJ2dvbffe7AwABarZaioiK/Xgv/NFnBEUIIIW5gra2t6gpKUVERgYGB+Hw+Dh48yJkzZ9QbU52dnXR2djI1NQVAYmIikZGRlJSUYDQaGRoa4tlnn6W9vZ2wsDAeeughCgsLL/tMRVE4efLkvJtaMTEx13Pa/5AEHCGEEOIG1d3dzdmzZwFYs2YN4eHhwMWr3SdOnECr1bJp0ybGxsbo7u5maGgIvV5PVFQUMTExlJSUEBQUxMjICM899xynT5/GYrFw3333UVpaetn+Uoqi0NDQQG9vLxqNhoKCgr/bC9JfJOAIIYQQN6ALFy7Q0NAAQGZmJomJiQA0NTXx4YcfoigKVqsVrVZLb28vPT09mEwmgoODSUhIoLi4GIvFwvj4OM899xzNzc0EBQVx9913s3HjRrTaSyPCXLjp7u5Go9GQn59PfHz8dZ33ZyUBRwghhLjB2O12Tpw4gaIoJCYmsnz5cgB6enrYv38/Xq+XlStXEhUVRW9vL+3t7QQGBmIwGEhOTqawsJCIiAjsdjt/+tOfaGhowGg0ctddd7F169YrNs9sbGxUw01eXh4JCQnXc9qfiwQcIYQQ4gbicrmoqqrC7XYTERHBmjVrgIs1cP785z/jdDpJSkoiIyODnp4e2traCAoKQlEU0tLSyM/PJyYmBofDwYsvvkhNTQ16vZ7t27ezffv2KzbPbGxspKurC7i4HTa3YrRQScARQgghbhBzBfVmZmYICgqiqKgIrVbL7Owsr776KhMTE0RERJCXl0dvby9nzpzBZDLhdDrJyMggOzubxMREZmdn2bVrF1VVVWg0GrZu3cqdd95JQEDAZZ/b1NREZ2cncDHcLFmy5DrO+ouRgCOEEELcIE6ePMno6CgGgwGr1UpAQABer5c9e/bQ39+PyWSirKyM/v5+2traMBgMzMzMsHz5crKyskhLS8PpdPLaa69x/PhxfD4fmzZt4qtf/eoV+0Y1NzfT0dEBQG5uLklJSddzyl+YBBwhhBDiBnD69Ol517KDg4NRFIV3332X06dPo9PpWL9+PUNDQ2ogmZ6eJjMzk8zMTFasWIHL5WLPnj0cOXIEj8fD+vXruffeezGbzZd9ZktLC+fOnQMgJyfnitWMFyIp9CeEEEIscHPbTXBxFWWuoF5FRYW6zVRSUsLU1BS9vb04nU4cDgfp6enq1pTb7Wb//v0cPnwYl8tFaWkp999//xU7g586dYr29nbgYtPOub5WNwpZwRFCCCEWsJGREerr6wFIT09Xt4iam5t5//33URSF7OxsfD4f/f39TExM4HA4WLp0KRkZGeTl5eH1ejlw4ADvv/8+DoeDoqIiHnzwQbVuzqe1traq9XVWr15NSkrKdZnrtSQBRwghhFig7HY7NpsNn89HQkICK1asAC5eB9+7dy9er5fU1FSCg4M5f/48w8PDOJ1O4uPjyczMpLCwEEVROHjwIO+++y7T09Pk5+fz9a9//YptFVpbW2lrawMuhpvU1NTrNt9rSQKOEEIIsQA5nU71Onh4eDhr1qxBo9EwOjrKq6++isvlIi4ujoSEBMbGxujt7cXr9RIZGcmKFSsoLi5Go9Hw/vvv88477zA1NUVOTg4PPvjgFSsPL5ZwAxJwhBBCiAXH6/Vech1cp9PhcDh4+eWXmZqaIjQ0lNTUVCYmJujo6MDr9RIcHMzKlSspKSlBp9Px4Ycf8vbbbzM2NsaqVat46KGHrnjFezGFG5CAI4QQQiwoc126x8bG1OvgRqMRr9fLK6+8wvDwMCaTiczMTBwOB+3t7fh8PoxGIytXrqS0tBSDwcCRI0c4cOAAFy5cICsriwcffPCKB4UXW7gBCThCCCHEgtLS0sLAwIDaHXzuOvi+ffvo7OxEp9ORlZWF1+ulra0Nn8+HVqtl1apVrFu3jsDAQD755BPefPNNhoeHyczM5P777yc9Pf2yz1uM4QYk4AghhBALRkdHh1p3Zs2aNURGRgLw4YcfUl9fj0ajITMzE61WS1tbG16vF5/Px8qVK1m3bh1ms5ny8nLeeOMNBgYGSE9P59577yUrK+uyz1us4QYk4AghhBALwtDQEM3NzQCsWLFC7fVUU1PD0aNHURSFpUuXEhwczLlz53C5XLjdbrKysrjlllsIDQ2lqqqK/fv309fXR1paGvfccw85OTmXfd6pU6cWbbgBCThCCCGE342Pj1NTU4OiKCQnJ5ORkQFAW1sbBw4cQFEUYmNjiYmJoauri5mZGVwuF5mZmdxyyy1ERkZis9nYt28fPT09pKSk8NWvfpW8vLzLPq+lpWVenZvFFm5AKhkLIYQQfuVwOKiursbr9RIdHU12djYA/f39/PnPf8br9WKxWEhOTqa3t5eJiQm1eebGjRuJjY2lpqaG/fv309XVxdKlS/nKV75CYWEhGo3mkuc1Nzer22DZ2dk3ZBG/z0ICjhBCCOEnbrebyspKnE4nFouFwsJCtFot4+Pj7N69G6fTSUBAAMuWLePChQuMjIzgdDpJSUlh48aNJCQkcPLkSfbu3cu5c+dYsmQJd955J1ar9R+Gm5ycnBuu/cLnIQFHCCGE8AOfz4fNZsNut2MymbBarej1emZnZ3nppZeYmppCo9GwfPly7HY7/f39uFwuEhIS2LRpE8nJyTQ0NKjhJjExkTvuuIOysjK02ktPoDQ1NalNOBd7uAE5gyOEEEJcd3O1bkZGRtDr9VitVkwmE16vl927d3P+/Hm8Xi/Lly/H7XbT3d2Ny+UiOjqazZs3k5aWRmNjI3v37qWtrY24uDh27NjB+vXrLwk3iqLQ2Niohpvc3NxFH25AAo4QQghx3Z06dYr+/n611s1cR+89e/bQ1dWFy+UiIyMDvV5PZ2cns7OzhIeHs3nzZjIzM2lubmbfvn2cOXOG2NhYtm/fzoYNG9DpdPOeMxduOjs7gYtXz5OTk6/3dP1CAo4QQghxHXV0dNDe3g5cDBxzTS8PHjxIU1MTTqeTpUuXEhISQnt7O7Ozs4SEhLBp0yZWrVpFS0sLe/fu5fTp00RHR3PbbbexadMm9Pr5p04URaGhoYGuri4A8vLy1E7kNwMJOEIIIcR1Mjg4SFNTEzC/1k1FRQXl5eU4nU7i4uKIi4vjzJkzOBwOgoKC2LRpE2vWrKG1tZV9+/bR2tpKZGQk27ZtY/PmzRgMhnnPURSF+vp6uru70Wg05OXlXbEH1WIlh4yFEEKI62B0dJSamhoAli5dqta6aW5u5uDBg7hcLiwWC6mpqWq4MRqNbNiwgaKiItra2ti7dy8tLS1ERESwbds2tmzZQkBAwLznzJ3v6evrU8PNXJC6mUjAEUIIIb5k09PT2Gw2fD4fsbGxaq2bzs5O/vKXv+ByudDr9axevZozZ85gt9sxGAysW7eOkpISzp49y1//+ldaWlrUszhbtmzBaDTOe47P56O2tpaBgQE0Gg0FBQXEx8f7Y8p+J1tUQgghxJfI6XRSWVmJy+UiLCyMgoICNBoNQ0NDvPzyyzidTrxeL4WFhZw7d46pqSl0Oh1lZWWsX7+e9vZ29u3bR0tLC2FhYdx6663cfvvtmEymec/x+XzU1NSojToLCwtv2nADEnCEEEKIL43H46GqqoqZmRnMZjPFxcXodDomJiZ48cUXsdvtTE9PU1JSQldXF2NjY2g0GoqKiti4cSPnzp1j3759NDY2qgeNrxRuTpw4weDgoHozKy4uzk+zXhgk4AghhBBfgrnQMTExQUBAAFarFaPRiMPh4MUXX2R8fJzx8XHWrl3L4OAgIyMj6pmZbdu20dnZyf79+2lqapoXboKCguY9x+v1Ul1dzdDQEDqdjuLiYmJiYvw064VDAo4QQghxjc3dYjp//jw6nQ6r1YrZbMbj8bB7926GhoYYGRmhrKyM8fFxBgYGgIuNL3fu3ElXVxf79++nsbGR4OBgNmzYwPbt2zGbzfOeMxdu/vY50dHR/pjygiOHjIUQQohrrLW1ld7eXjQaDYWFhYSFheHz+fjzn/9MZ2cnw8PDlJSUMDs7S09Pj9qS4Stf+Qrd3d3s37+fhoYGNdzs2LGD4ODgec+Y2/4aHR1VqyFHRET4acYLjwQcIYQQ4hrq6Ojg7NmzwMW2CHPbRQcOHKClpYWhoSHy8vLQaDR0dXWh0WhIS0vja1/7Gj09Pbzxxhvqys369evZsWMHISEh854x16RzfHwcg8GA1WolPDz8us91IZOAI4QQQlwj/f398wr5zVUO/uijj7DZbAwNDZGVlUVwcDBtbW1oNBqSkpJ44IEH6Ovr48CBAzQ0NGA2m1m/fj07d+68JNy4XC4qKiqYnJwkICCAkpISQkNDr/tcFzoJOEIIIcQ1cOHCBerq6gBISUlRC/lVV1dz+PBhhoaGSElJITY2ltOnT6PRaIiPj+ehhx6iv7+fAwcOcPLkSYKDg1m7di07d+5Ue1TNmZ2dpaKiArvdjtFopLS09JIAJC5acIeMf/nLX6LRaOa9/vaqm6Io/PKXvyQhIYHAwEA2btxIc3OzH0cshBDiZjc5OakW8ouPj2f16tUANDY28vbbbzM0NERMTIxapRggOjqahx9+mIGBgXnhpqys7LLhxuFwUF5ejt1ux2QysXbtWgk3f8eCCzgAq1atYmBgQH01NjaqX/v1r3/Nb3/7W37/+99js9mIi4tj69atTE1N+XHEQgghblYzMzNUVVXh8XiIiIggPz8fjUajFugbGhrCYrGwatUqzpw5g6IohIeH881vfpPBwUHeeust6uvrCQ4OprS0lJ07d16y5TQ9Pc3x48eZnp4mKCiItWvXXnKjSsy3ILeo9Hr9ZQsUKYrCU089xS9+8QvuueceAF588UViY2N55ZVX+O53v3u9hyqEEOImNleleHZ2FovFQnFxMVqtlr6+Pl577TUGBgYICAigoKCA06dP4/P5sFgsfPOb32RoaIh33nlHPXNTUlLCzp07CQsLm/cMu91ORUUFs7OzmM1mysrKLin0Jy61IFdw2traSEhIIDU1lQcffJBz584BF0+mDw4Osm3bNvV75xqRlZeX/93PdDqdTE5OznsJIYQQX9TcNe3p6WkCAwOxWq0YDAYuXLjAyy+/TF9fHwAlJSW0tbXh9Xoxm808+uijXLhwQQ03QUFBlJSUcMcdd1xyE2pycpLy8nJmZ2cJCQlh7dq1Em4+owUXcKxWKy+99BLvvfcezz77LIODg5SVlTEyMsLg4CAAsbGx834mNjZW/dqVPPnkk4SGhqqvuZPtQgghxOf16SrFpaWlmEwmpqam2LVrF11dXbhcLsrKyjh79ixutxuTycQjjzzCyMgIBw8enLdyc7lwMz4+Tnl5OU6nk9DQUMrKyi5primubMEFnO3bt/O1r32N7OxstmzZwttvvw1c3Iqao9Fo5v2MoiiXvPdpjz/+OBMTE+qrp6fn2g9eCCHEoqcoCnV1dZdUKXY4HOzatYv29namp6dZt24dnZ2duN1ujEYjjzzyCBMTE7z77rvqmRur1XrZcDMyMkJFRQVut5uIiAhKS0sJCAjw04xvTAsu4Hya2WwmOzubtrY29VzOp1drhoeHL1nV+TSj0YjFYpn3EkIIIT6v5uZm+vv71aaWYWFhuN1uXn31VU6dOsXExATr1q2jt7cXl8uFwWDg4YcfZnJykvfee4/6+npCQkIoLi6+bLgZHh6msrISj8dDVK8bKC8AACAASURBVFSUuvUlPp8FH3CcTienTp0iPj6e1NRU4uLiOHTokPp1l8vFkSNHKCsr8+MohRBC3Aza2tro6OgAYM2aNURHR+PxePjLX/7CyZMnGR0dpbS0lOHhYZxOJ3q9nm984xvY7XYOHTr0D8PNwMCAet08NjYWq9WKXr8g7wMteAvub+2nP/0pd955J8nJyQwPD/OrX/2KyclJHn30UTQaDT/+8Y954oknyMjIICMjgyeeeIKgoCAeeughfw9dCCHEItbd3U1raytwsSlmYmIiPp+PN998U214abVaGR8fx+FwoNPpeOCBB7Db7Rw+fJj6+nosFgtFRUWXDTe9vb2cPHkSRVFISEggLy8PrXbBr0MsWAsu4PT29vL1r3+dCxcuEB0dTUlJCZWVlSxduhSAn/3sZzgcDr7//e8zNjaG1Wrl/fffl2JHQgghvjQDAwM0NDQAkJGRQWpqKoqi8O6773L06FGGhobIz8/H6XSq4ebee+9ldnaWw4cP09DQQGhoKEVFRezcufOScNPV1aV+fnJyMjk5Of/wbKn4+zSKoij+HoQ/TE5OEhoaysTEhJzHEUIIcUXnz5+nuroan89HcnIyubm5KIrCRx99xIEDBxgYGCA3Nxe9Xs/ExAQ6nY6vfe1ruN1uPv74YxobGwkLC6OwsPCy4aa9vZ2WlhYAUlNTWbVqlYSbv+Oz/v5ecCs4QgghxEIxPj4+rwVDTk4OABUVFbz99tsMDAywYsUKDAYD4+PjaLVa7rzzTpxOJ0ePHlXDTVFRETt27Lgk3LS2ttLW1gZcXBlasWLFdZ/jYiUBRwghhLiMqakpKisr8Xq9REdHqy0Y6urq2Lt3L319fSxbtozg4GA13OzYsQOv18vRo0dpbm4mIiKCwsJCduzYMa9CsaIoNDU10dnZCUBWVhbp6el+muniJAFHCCGE+BSHw0FlZSVut1vdXtJqtbS0tPDaa6/R29tLUlIS4eHharjZtm0biqKo4SYqKkoNN3/bW0pRFE6ePElvby8AOTk56jlTce1IwBFCCCH+htPpVHs/zRXj0+v1tLe3s3v3brq6uoiLiyM2NlYNN1u2bEGr1XLs2DFaWlqIjo6mqKiI7du3zzsn4vP5qKmpYXBwEI1GQ15eHomJiX6c7eIlAUcIIYT4v9xu97z+UnMVhHt6enjhhRc4d+4cUVFRJCQkMDY2hk6nY9OmTWg0Go4ePUpraytRUVGXDTcejwebzcaFCxfQarUUFhb+wyK14ouTgCOEEEIAXq8Xm812SX+pgYEBnn/+edra2ggLC2Pp0qXqban169ej0Wg4fvw4ra2txMTEUFRUxO233z6vfMlccBobG0Ov11NUVERUVJQfZ7v4ScARQghx05vbOhoZGUGv11NSUoLZbGZ4eJjnnnuO1tZWLBYLqampTExMoNVqKS0tRaPRUF5ezpkzZ4iLi6OwsJDbb7+d4OBg9bNnZ2epqqpicnISg8FASUnJvAPH4sshAUcIIcRNbe7Q79DQEFqtluLiYkJDQxkdHeX555+nubkZs9lMWloadrtd/R6dTkdlZSVnz54lPj6eoqIitm3bNi/cTE9PU1lZyczMDEajkdLSUilMe51IwBFCCHFTa2xspK+vD41GQ1FREZGRkUxOTvL8889TX1+PyWQiLS2N6elpNBoNBQUF6PV6KisrOXfuHImJiRQVFbF161bMZrP6uZOTk1RWVuJ0OgkKCqK0tJSgoCA/zvTmIgFHCCHETevUqVN0dXUBkJ+fT0xMDHa7neeff56amhoCAgJIS0vD6XSqt570ej1VVVV0dHSwZMkSiouL2bp1K4GBgernjo6OUl1djdvtxmKxYLVaMZlM/prmTUkCjhBCiJtSW1sbZ8+eBSA3N5eEhARmZmZ44YUXqKqqwmAwkJqaisfjUb9Hp9NRXV1Nd3c3SUlJWK1WtmzZMi+8nD9/HpvNhtfrJTw8HKvVisFg8Mscb2YScIQQQtx0Ojs71c7gq1atIjk5mdnZWXbt2sUnn3yCXq8nOTmZuXaNOTk56HQ6bDYbfX19JCcnU1JSwq233orRaFQ/t7+/n7q6Onw+HzExMRQWFqLT6fwyx5udBBwhhBA3ld7eXhobGwHIzMxUt6B2797NRx99hFarJTExEb1ej6IorFy5EoATJ07Q399PcnIypaWlbNq0iYCAAPVz/7YjeEJCAnl5eWi12us/QQFIwBFCCHETGRwc5OTJk8DFzt3Lly/H5XLx6quv8sEHH6DRaEhMTFS3nJYvX45Wq6W2tpahoSFSUlIoKytjw4YN87ad2tra1BWhpUuXkp2dLR3B/UwCjhBCiJvC+fPnqampQVEUkpKSWJG1kmOtg/x1z59prjpChNlIQny8etNp2bJl6HQ6amtruXDhAqmpqZSVlbF+/Xo13CiKQnNzMx0dHYB0BF9IJOAIIYRY9EZHR7HZbPh8PuLj4xnQxfC9Jz/gbNX7zLTbQFGwhEexORwsFg0pKSnqys3Y2BhpaWmsW7eOsrIy9PqLvzp9Ph8nT56kr68PuHiWJy0tzZ/TFH9DNgeFEEIsauPj41RVVeH1eomJiWE4IJ7/s+sEZ6s/wNF+AhQfusAwHBojbzcOMqG1oNVqqaurY3x8nGXLlrFhwwbWrl2rhhuv10t1dbVaPyc/P1/CzQIjKzhCCCEWrbliex6Ph6ioKPLyC7jl14exnynHcbYKRfGhNVnQBprRaDTogiN4s2EQ32gPzlkH6enp3HLLLRQXF6sHhl0uF9XV1WqzzcLCQmJiYvw8U/FpEnCEEEIsSna7ncrKStxuN+Hh4RQVFVHVMUa77WMcZyrxKT50AcHoAi1oNFq05jDQ6hjubKXbGMnavFVs3LiRgoIC9cCww+GgsrISu92OwWDAarUSHh7u55mKy5GAI4QQYtGZmZmhoqICp9NJaGgoVqsVrVbLWwf2M9NWgc/rQWcKRmcOBa0WbZAFtFpcQ+2gKITHJbNlyxZyc3PVcGO326moqGB2dhaTyURJSYn0lVrAJOAIIYRYVGZnZ9UgEhISQklJCTqdjjfeeIO6o+/j87jRGYPQ/d8VG50pBDRa3EPnQKNFHx7PrZu3smbNGvUzx8bGqKqqwu12ExwcTElJybzWDGLhkYAjhBBi0XA6nVRUVDAzM4PZbKakpASDwcD+/fvZu3cvQTqF4OBgXMa5cGMGrQbXcCcanZ6A8DiWrC7h67evVT9zaGiImpoavF4vYWFhWK3WeQX+xMIkt6iEEEIsCi6Xi4qKCux2O4GBgZSWlmI0Gtm/fz/79u3D4XAQFmphS8EKNDod2oBAfID7fCdavYGAiAQC04r4r8duR6e9uC3V3d2t9pWKiYmhrKxMws0NQlZwhBBC3PDcbjcVFRVMTU1hMpkoLS3FZDKxf/9+9u/fz/T0NCEhISQkJGAwGAgMDubD1guM9HeiDQhCHx7P0ty1PPHNTdy+Oh6YX504KSmJnJwcab1wA5GAI4QQ4obmdruprKxkcnISo9FIaWkpQUFB7Nu3jzfffJOpqSksFosabkwmE1lBQWinhnAkZhCesJRt225nZ+kqdFrNJdWJ09PTycrK8vMsxeclAUcIIcQNy+PxUF1dzfj4OAEBAZSWlmI2m9m7dy9vvfUWExMTargxmUwEBASgKAptbW1ERUaSlZXF7bffTmxsLHCxOnFdXR39/f0ArF69mtTUVH9OUXxBV73WtmXLFg4ePHjJ+16v92o/WgghhLiiuWrCo6OjGAwGSkpKCA4O5q9//Stvv/024+PjhISEEB8fj9FoRK/X4/F4OHfuHJGRkWRnZ7Nz50413Ljdbqqqqujv70er1VJQUCDh5gZ21QHnxIkTpKSkAKjLeQDPP/88jzzyyNV+vBBCCHEJn8+HzWZjZGQEvV6v1qTZs2cPBw8eZGRkRD1zYzKZ1HDT3d1NdHQ0ubm57Nixg6ioKODi1fLy8nIuXLiAXq/HarWSkJDg51mKq3HVAcflcqmFjnJzczl37hwAZWVlfPjhh1f78UIIIcQ8c+Hm/Pnz6HQ6rFYrFouFv/zlL7z33nucP38ei8VCfHw8AQEBaLVanE4nfX19xMXFkZ+fz44dO9QKxFNTU3zyySfqGZ6ysjI1+Igb11WfwUlPT6eqqoqQkBCmp6cZHx8HICQkhNHR0aseoBBCCDHH5/NRU1PD8PCwGm5CQ0N5/fXX+eCDDxgeHiYsLIzY2FiMRiM6nQ6Hw8HY2Bjx8fHk5+ezefNmzGYzcLHLeHV1tVrAz2q1EhQU5OdZimvhqldwvv/97/Otb32LDRs2kJubyzPPPAPAsWPH1H1NIYQQ4mr5fD5qa2sZHBxEq9VSVFREaGgor776Kh988AGDg4OEhYURExODyWRCq9UyMzPD+Pg4iYmJlJSUsG3bNjXcDAwMUFFRgdvtJiIigrVr10q4WUSuegXne9/7HtHR0bS1tfHtb3+bBx98kLS0NAYGBvjnf/7nazFGIYQQN7m5cDMwMKCGm7CwMF555RWOHTvGwMAAUVFRREVFYTKZ0Ol02O12HA4HSUlJWK1WbrnlFgwGA3DxzGhTUxOAum2l0+n8OUVxjWkURVGu5Qd6PB727duHy+XiwQcfXLD/wUxOThIaGqpeIRRCCLEwKYpCbW2tertpLty8/PLLlJeX09fXR3R0NBEREQQGBqLX65mcnMTr9bJkyRLWrl1LSUkJer0eRVE4deoU7e3tAKSkpLB69Wq1oaZY+D7r7+9rXgdHr9dz3333XeuPFUIIcRNSFEWtS6PVaiksLCQ0NJTdu3dTUVFBT08PsbGxhIeHq2duxsbG0Ol0JCcns2HDBgoLC9Fqtfh8Pk6ePElfXx8AWVlZpKen+3mG4ssihf6EEEIsSHPhpq+vD41GQ0FBARaLhd27d1NZWUlPTw9xcXFYLBb1zM3Y2BgBAQEkJSWxceNG1qxZg0ajwe12q9fKNRoNa9asYcmSJf6eovgSScARQgix4CiKoq62aDQaCgsLCQkJYdeuXdhsNrq7u4mLiyMkJITAwEAURWF8fJygoCCSk5PZsmULK1euBMDhcFBZWYndbkev11NYWEh0dLSfZyi+bNe0a9jx48dxOp3X8iOFEELcZBRFob6+nt7eXnXlJigoSA03cys3ZrOZwMBAvF4vk5OThISEsGzZMu644w413ExMTHDs2DHsdjsmk4m1a9dKuLlJXNOAs337dnVvUwghhPi85sJNT08PGo2G/Px8TCYTL730EidOnFBXboKCgjCbzXg8Hux2O2FhYWRkZHDnnXeSlpYGwPDwsPoPb4vFwvr16+VSyU3kmm5RXeMLWUIIIW4inw43eXl5GAwGdu3apb6fkJCA0WjEbDbjcrlwOp1ERkayfPlybrvtNrX+Wnd3Nw0NDSiKQnR0NAUFBeoVcXFzkDM4Qggh/G7uzM3ctlR+fj4ajYbdu3fT0NBAf38/8fHxGAwGgoODmZ2dxev1Eh0dzapVq9i2bRsREREAtLa20tbWBkBSUhI5OTlotdd0w0LcAK4q4Lz00kvz/uzxeNi7dy8xMTHqe9/85jev5hE8+eST/PznP+dHP/oRTz31FABOp5Of/vSnvPrqqzgcDjZv3swf/vAHOREvhBA3oMuFG6/Xy+uvv05DQwNDQ0PExcWh1+sJDQ1lenoarVZLbGwsubm5bNmyhZCQkEuugWdmZrJ8+XI/z074y1UFnD/96U/z/ux2u9mzZw+BgYEAaDSaqwo4NpuNZ555hpycnHnv//jHP+bAgQO89tprREZG8pOf/IQ77riDmpqaBVtYUAghxKUuF26cTid79uyhsbGR4eFhoqOj0el0hIWFYbfbMRgMxMTEkJ+fz6233kpgYOAl18BzcnJITk729/SEH11VwPnoo4/m/TkkJIRXXnlFPeB1Nex2O9/4xjd49tln+dWvfqW+PzExwfPPP8+uXbvYsmULALt37yYpKYkPPviA22677aqfLYQQ4st3uTo3ExMTvPnmmzQ0NDAyMkJkZKS6cjMxMUFwcDBRUdFYkrNwxKzkZP8Mq2K8nLBVMz09LdfAhWrBnsH5wQ9+wM6dO9myZcu8gFNTU4Pb7Wbbtm3qewkJCaxevZry8vIrBhyn0znvCvvk5OSXN3ghhBB/16fDTX5+PsPDw7z77rvU19czPj5OREQEer0ei8WilucfJ4hDPcHYJ/VoWlrwztoJnujggYJ41q1IxGq1EhIS4u/piQVgQQac1157jdraWmw22yVfGxwcJCAggPDw8Hnvx8bGMjg4eMXPfPLJJ/n3f//3az5WIYQQn8/lwk1vby+HDx+mtraWmZkZLBaLeqB4amqKiIgIJjRmDg6HExCXhkajxTM1gmv4HC7Fx/9bMUh+6S0SboTqmh4r//nPf66eYv+ienp6+NGPfsTu3bsxmUyf+ecURfm7zdIef/xxJiYm1FdPT89VjVMIIcTn5/P5qKmpoa+vD61WS15eHh0dHRw6dAibzcbs7CxBQUEEBgYSFBTE9PQ00dHRLElKptqdjDE+A41Gi3usH9fQWVB8aIPCMCZm8eT75/D6pFyJuOiaBpzHH3+csLCwq/qMmpoahoeHKSgoQK/Xo9frOXLkCP/93/+NXq8nNjYWl8vF2NjYvJ8bHh5W6x9cjtFoxGKxzHsJIYS4fubCzcDAAFqtltzcXM6cOcOHH37IiRMn8Hq9GAwGQkJC0Ov1OBwO4uLiSE1NJSl3HROmOBTFh2u4A/fIxX+k6kNjCYjPBK2OgYlZqjtG/TxLsVAsuMIAmzdvprGxkZMnT6qvwsJCvvGNb6j/22AwcOjQIfVnBgYGaGpqoqyszI8jF0IIcSVerxebzcbg4KAablpaWvj444+pqalBo9GgKArh4eEoioKiKCQkJJCZmcldd92FPiwWxevB1X8Gz+QwAIaopQREp8xbvR+emvXXFMUCs+DO4ISEhLB69ep575nNZiIjI9X3H3vsMX7yk58QGRlJREQEP/3pT8nOzlZvVQkhhFg45sLN+fPn0el0rF69moaGBiorK6mrqyMgIACn00l8fDxOpxOTyURsbCwrVqzgtttuIyoqCstoD86+FnwuBxqNloDYdHTB4Zc8Kybksx9tEIvbggs4n8Xvfvc79Ho9999/v1ro74UXXpAaOEIIscB4PB6qq6sZGRlBp9OxcuVK6urqqK6upr6+HpPJxPT0NElJSdjtdkJDQ4mLi1P/0WqxWBgbG8PR00yY3su4N4CAhEy0RvO852iAuFATxalXdw5ULB4a5Ro0kHK73QwODjIzM0N0dPRVHzS+HuauHE5MTMh5HCGE+BJ4PB6qqqoYHR1Fr9eTkZFBfX09NpuNpqYmjEYj09PTJCcnMzk5SXR0NLGxsRQUFLBx40YCAwPp7++nrq4On89H64iH3zeBVh/A3/7imtug+uPD+dy+Ot4fUxXX0Wf9/f2FV3Dsdjsvv/wyr776KtXV1fNqzCxZsoRt27bxne98h6Kioi/6CCGEEDcot9tNVVUVY2NjGAwGUlJSqKmpwWazcerUKUwmE3a7neTkZCYmJoiPjycuLo6SkhLWrl2LwWCgra2N1tZWAOLi4ti+PY/s1vP8+4EWBib+/7M2caEm/p87V0q4EfN8oRWc3/3ud/znf/4nKSkp3HXXXRQXF5OYmEhgYCCjo6M0NTVx7Ngx9u3bR0lJCU8//TQZGRlfxvi/MFnBEUKIL4fL5aKyspKJiQkMBgOJiYk0NDRQXV1Ne3s7AQEB6rbU5OQkSUlJxMfHs379evUfxfX19fT29gKQlpbGypUr1cPEXp9Cdccow1OzxIRc3JbSaa9cJkQsLp/19/cXCjj33Xcf//Zv/0Z2dvbf/T6n08nzzz9PQEAA3/rWtz7vY75UEnCEEOLam52dpbKykqmpKYxGI5GRkTQ3N1NZWUlPTw9arVY9UDw9PU1KSgpLlizh1ltvZdWqVWpPqdHRUTQaDdnZ2SxdutTf0xILyJcacBYDCThCCHFtORwOKioqmJ6eVmuPnTp1ioqKCgYHB/F6vXi9XqKiovB4PKSmprJ06VK2bdtGWloaU1NTVFdXMzMzIz2lxBV96WdwhBBCiDnT09NUVFTgcDgwGo0EBgbS1NTE8ePHGRsbw+VyodPp1Do3y5YtIz09ndtvv534+HiGh4epqanB4/EQFBSE1WolODjY39MSN7AvFHDefPPNz/0zW7duJTAw8Is8TgghxAJmt9upqKhgdnYWo9GIXq+nqamJ8vJypqammJqawmw2YzQaCQgIIDU1lRUrVrBt2zYiIiLo6OigubkZRVGIjIyksLCQgIAAf09L3OC+UMC5++67P9f3azQa2traSEtL+yKPE0IIsUBNTk5SUVGBy+UiICAARVFoamqiqqqK6elpxsfHCQ0NRafTYbFYWLp0KXl5edx6660EBgbS2NhIZ2cnAElJSeTk5KDVLrgi++IG9IW3qAYHB4mJiflM3yvdXYUQYvEZGxujqqoKt9uNXq/H4/HQ0tKidgQfHx9X+xNGR0eTkpJCcXEx69evR6PRUFVVxYULFwBYuXIly5Yt8+d0xCLzhQLOo48++rm2mx5++GE5yCuEEIvIyMgI1dXVeDwetFotLpeLxsZGmpubsdvtTE1NERYWhqIoJCcns3TpUvUa+OzsrLrCo9PpyM/PJy4uzt9TEouM3KKSW1RCCPG5DA0NceLECXw+H4qi4PV6qa2tpa2tTT1zExoaikajYdmyZaSmpqrXwC9cuEBNTQ1ut5vAwECKi4vl/4PF5/Kl3aJyOByMjo6SmJg47/3m5mZWrVr1+UcqhBDihtHX10ddXR2KouDxePD5fFRVVdHd3c3o6CgulwuLxYJeryc9PZ309HRuu+02UlJS5h0mDg8Pp6ioCKPR6O8piUXqc53k2rNnD5mZmezYsYOcnByqqqrUrz3yyCPXfHBCCCEWjq6uLmpra/H5fMzOzuJ2uzl27Bg9PT0MDQ3h8XgwmUyYTCaysrLIzs7m7rvvJjk5mfr6epqamlAUhaSkJMrKyiTciC/V51rB+dWvfkVtbS3R0dGcOHGCRx99lF/84hc89NBD3KQ7XUIIcVNob2+npaUFr9fLzMwMWq2WI0eOMD4+Tm9vL0ajEa1WS2hoKJmZmaxZs4YtW7YQGBhIZWUlIyMjgBwmFtfP5wo4brdbrSpZWFjI0aNHueeeezh79qzaI0QIIcSN7dO9nizOYc61n8XtdjMzM4OiKBw+fJjp6Wm6u7sJCgoCICYmhszMTIqLi9m4cSMul4tjx46plYnz8/OJjY318+zEzeJzBZyYmBgaGhrIyckBIDIykkOHDvHoo4/S0NDwpQxQCCHE9fNu04DarVtRFNwXurC4x7gnN4asmMB5oaW3t5egoCB122n58uWsX78eq9XKhQsXqKurw+PxYDabKSoqkpIh4rr6XLeoent70ev1l73Od/z4cdauXXtNB/dlkltUQggx37tNA/yf3bUogKL4cA934Jm6gOJyoPi83J8bhb2nlampKQYHBwkMDESr1bJs2TKysrK49dZbyc7O5uzZs7S2tgIQFRVFQUGBVCYW18yXcotqyZIll7z31FNP8cADD9xQ4UYIIcR8Xp/Cvx9ouRhufF5cg2fxzozjm50GxYd3ZoLX32nEmmBkbGwUk8mEwWAgIyOD1atXc9ttt5GcnPz/sXff0XHe953v38/0ATAz6DPoBAkQlQRBEmxiUahCVVuyEtnx8drJOtlYLrmOLe+9dnajo13bOo58Y2dTtE68LrEiS3djyZISWY4qaZEUQQJsIFhAECTaABhM7+157h9cjkWxiAIJAiC/r3P4BzCD5/k9PBLxmV/5funt7WV8fByARYsW0dbWJpWJxZy46v/qvvKVr3DLLbcwOjp63vdTqRR79+692ssLIYS4DrqHfGeXpbIZkuMnyMYCZGNBQCMdnCI1dZqIf5ph9xQGgwGr1UpbWxtr167lwQcfxOVy8c477zA+Po5Op6Ojo4Nly5ZJuBFz5pr8l3f33XezefPm80KO3+9n3bp11+LyQgghZtlUOIGWSZMcO0o2HiQbmkZT9KSmh8n4x8mEPKiZJBlVo7CwkOXLl7Nx40Y++tGPYjAY2LFjB6FQCLPZzPr166mtrZ3rRxI3uRn3ojpHURQee+wxysvL2bx5Mzt27MgtZcnRcSGEWBjsBpXEWD9qMkY25EGxFJCeGEBNRMgGpwANBXA6XXR2drJx40Y2b96M2+3O1bdxOBx0dXV9qFY+QsyWqw445zz22GMAuZBjNBrl6LgQQiwA4XCY5NhRHLok04EJFIuN5OhRtHScjN+NpuhQFI2C4gru2LKB22/bysqVKzl69ChnzpwBoKqqio6ODvR6/Rw/jRBnXXXAee8szXtDznPPPXe1lxZCCDHLznUEDwaDbK7QeN5rJTl6BDWdIB2cRFF0KAoYCiv43Xtu5Xcf+hiLFi2iu7sbn88HQEtLCw0NDXP8JEKc76oDzre+9S3y8/NzX58LOffee+/VXloIIcQsmpqaYt++fXi9XiYnJ1mxqJwzA8fYl4kTDU6hKAqKTo/DWcMn77+dr/3RJ8jLy2PHjh0kEgmMRiMrV66kvLx8rh9FiAvMWjfxb37zm3z3u98lEAjMxuWvmtTBEULczEZHRzlw4AATExMEAgFUVaWvrw+/38/k5CSxlIZiNNLU1MJD993JXdvuJBwOc/jwYVRVxWaz0dXVdd4HXCGuhyv9/T1rAWe+k4AjhLhZDQ4OcuTIEUZGRkgmk4RCIQYHB/F4PPh8PhRFwW63s2LFCrZu3crmzZs5deoUp0+fBsDlOrvR2GC4Zts4hbhis1LoTwghxMKlaRr9/f0MDAxw8uRJDAYDk5OTjI6OMjk5SSgUQlEUSkpKWLNmDdu2baOjo4P9+/fn9ts0NzfT0NAgh0jEvCcBRwghbgKqqnLgwAFOnTrFDBYrAgAAIABJREFUwMAAVquVwcFBpqencbvdRKNRdDodlZWVbNiwgfvuuw+Xy8XOnTtJJBLSLFMsOBJwhBDiBpfJZNi3bx9DQ0MMDAxQUFDA0aNHCQQCjI+PE4/HMRgM1NXVcdttt3HfffeRyWTYtWsXqqpSUFBAV1cXBQUFc/0oQlyxWa2hrdPp2Lp1Kz09PbN5GyGEuGllVY3dg15ePDDG7kEvWfX8bZWpVIrdu3dz7Ngxjh8/jtls5vDhw3i9Xs6cOUM8HsdkMtHU1MTHPvYxfu/3fg+Px8PBgwdRVRWXy8WmTZsk3IgFZ1ZncH70ox9x5swZ/vRP/5SdO3fO5q2EEOKmklU1/vbNAX688zSBeDr3/QqHhcfub+Wu9gqi0Sh79uzhxIkTuN1udDodhw8fJhKJMD4+jqqqWK1W2tvbeeCBB1izZg19fX2506+y30YsZHKKSk5RCSEWmFf73Pw/zx8mEEtf8Nq5KPLkR5aQHxzi+PHj+P1+0uk0Q0NDhMNhJiYmACgoKGDt2rU89NBD1NfX09vbSyqVwmg0smrVKsrKyq7jUwlxZa7bKarvf//7PPzww1RWVl7tpYQQQnyAV/vcPPJ0L5f6ZKoB2WiAr/3tz9lYksasZCnQpxgfG8uFG6PRiMPhYOvWrTz88MMoisKePXvQNA273U5XVxd5eXnX87GEuOauegZHp9NRX1/P9u3bc0024ey678GDB+nq6rrqQc4GmcERQiw0WVVj43fexB1MXPI9mdAUibFjZHxjaGoWNRbCmI5QW6CixUNYrVZKSkr4yEc+wkc+8hHcbjdjY2MAVFdXs3z5cuknJea161oH56677rqgk7jf72fdunVks9lrcQshhLjpdQ/5Lhtu0t5REuPHyAQmzoabaAA1ESUVD9I3GaWhsoTG6mo++clPcsstt3D06FHC4TCKotDW1kZ9ff11fBohZtdVBxxFUXjssccoLy+/IOTcpNt7hBBiVkyFLx5uNE0jNXWK5PgxMiEvmppBjfjIppNkIz60VBK9xYpXX8IXvvglGpYspre3l0wmg8ViYdWqVRQXF1/npxFidl2zU1Tv7SS+Y8cOjEaj7LwXQohrqNxmueB7mpolOX6CpPs42WgALZ0iG/WjqmnU4BQAekseJtcSDMvvYjxuIHLgAAAlJSWsWrUKs9l8XZ9DiOvhqgPOe2dp3htynnvuuau9tBBCiPdYU19MhcPCRDCBBmiZFImRIyQnBtASEdRUHDUeQs2kyIQ8KDojepMZy6IO8tu2kg376D8xwLrFJSxZsoSWlhb5ICpuWFcdcL71rW+d1032XMi59957r/bSQggh3kOvU3js/tazp6iSMWKnD5CePo2aiKEmo6jJCGoqQSY8jd5kRTHnkde8CeuizrObjrMpSmw1rF69moqKirl+HCFm1azVwfnmN7/Jd7/73VzBqPlGTlEJIRaq537Txzf++md4xs+gJaOoiSgGLUU6GiYV8aM356HPL6RgxV0YCytI+8dR0CgpcrD7//1POOy2uX4EIWbsSn9/S6E/CThCiAXkzJkzvPzyy4yOjnH8jBuv348OlURgmgmvD3dcj8HhpGDVR1AANeZHAfQFJfzjlz/GvStq5voRhLgq1/WYuBBCiNl3+PBhfvWrXzE9Pc309DQmNU5pnp7Tp0fQNI0aVzmdi1oYKtuA2z2Omk6AoqNyUQNPfOY27mqXZSlx85iXAeepp57iqaee4vTp0wC0tbXxF3/xF9x9990AJJNJHn30UX7+858Tj8e57bbb+Pu///vzCg0KIcSNQlVVdu7cyfbt2/F6vXi9XuLxOMlkkqGhIcxmM/n5+dx7773cf//9DJwc5LjbSFzVs3ZNF7evWIxeJ5uJxc1lXi5Rvfzyy+j1ehoaGgD46U9/ypNPPsn+/ftpa2vjkUce4eWXX+YnP/kJJSUlfPWrX8Xn89HT03PFFThliUoIsRCk02l+9atfsW/fPjweD36/n0wmQzgcZmRkBLvdjt1u59Of/jTNzc243W4AysvL6ezsxGQyzfETCHFtzYs9ODqdjltvvZUnn3ySVatWXdW1iouLefLJJ/nd3/1dysrK+NnPfsbHP/5xAMbHx6mpqeGVV15h27ZtF/35ZDJJMpnMfR0KhaipqZGAI4SYtyKRCP/yL//CwMAA4+PjRKNRVFVlamqK6elpCgsLcTqdfOlLX0Kn0+WqEjc1NUkXcHHDutKAo5vNQfzoRz9iy5Yt/Omf/umMr5HNZnn22WeJRqOsX7+enp4e0uk0d955Z+49lZWVtLe3s2vXrkte54knnsDhcOT+1NTIRjshxPw1OTnJD3/4Q44ePcqpU6eIRCJks1mGh4fx+/0UFxfT2trKN77xDZLJJOFwGLPZzLp162hsbJRwI25683KJCs5uplu/fj2JRIKCggKeeeYZ7rnnHp555hn+8A//8LzZGIA777yT+vp6fvCDH1z0ejKDI4RYKI4dO8bzzz+Px+NheHgYRVFIp9OcPn0anU6H3W5n8+bN3HfffUxNna1WXFxczKpVq7BYLqx2LMSNZMGfompqauLAgQMEAgF+8Ytf8JnPfIbt27df8v2apl32E4vZbJZy5EKIeW/Hjh28/vrrTE5OMjExgclkIhqNMjw8TF5eHna7nYceeojm5uZcuGloaKC5uVlmbYR4j6sKOKOjozz11FPs2rWLiYkJFEXB6XSyYcMGPve5z13VMpDJZMptMl69ejV79+7lr//6r/n4xz9OKpXC7/dTVFSUe//U1BQbNmy4mscRQog5k81m+eUvf0lPTw8jIyOEQiGMRiMejwePx4PdbqesrIzPfOYz5OXlEYlEMJlMdHZ2Ul5ePtfDF2LemfEenHfeeYeWlhZeeOEFOjo6+PSnP82nPvUpOjo6+OUvf0lbWxs7d+68ZgPVNI1kMsmqVaswGo289tprudfcbjd9fX0ScIQQC1IkEuGHP/whe/bs4dixY4TDYfR6PWNjY7nNxI2NjTzyyCOYTCYymQzFxcVs3rxZwo0QlzDjGZw/+7M/44/+6I/43ve+d8nXv/zlL7N3794Pfe1vfOMb3H333dTU1BAOh3n22Wd5++23efXVV3E4HHz2s5/lq1/9KiUlJRQXF/Poo4+ybNkybr/99pk+jhBCzImRkRF+/vOfMzo6yqlTp7BYLKTTaUZGRshkMpSUlLBmzRq2bt1KJpMBzi5JNTU1odPN6jkRIRa0GQecvr4+nn766Uu+/id/8if8z//5P2d07cnJSf7Df/gPuN1uHA4Hy5cv59VXX+WOO+4A4Hvf+x4Gg4GHH344V+jvJz/5yRXXwBFCiPmgt7eXF198kbGxMUZHR7HZbASDQdxuN2azGafTye23305rayuZTAaj0UhnZydOp3Ouhy7EvDfjgFNRUcGuXbtoamq66Ou7d++ecbfa//W//tdlX7dYLPzN3/wNf/M3fzOj6wshxFxSVZVXXnmFnTt3curUKUKhEAUFBUxNTeHz+cjPz8flcnHffffhdDrRNI2ioiJWrVqF1Wqd6+ELsSDMOOA8+uijfO5zn6Onp4c77rgDp9OJoihMTEzw2muv8cMf/pDvf//713KsQgix4MViMZ577jkOHz7MiRMnUBQFk8nE+Pg44XAYh8PBkiVLuPPOO3NHYBsbG1m6dKksSQnxIcw44Hz+85+npKSE733ve/zgBz8gm80CoNfrWbVqFf/0T//Eww8/fM0GKoQQC93Y2BjPPfccJ06cYHBwkIKCAhKJBB6Ph2w2S2lpKStWrGDdunXk5eVhNpvp7OykrKxsrocuxIJzTQr9pdNppqenASgtLcVoNF71wGab9KISQlxP+/fv56WXXuLEiRNMTk5it9sJBAIEAgEURaG8vJyuri5WrFiB0WiktLSUzs5OKdwnxPtc10J/RqNxxvtthBDiRpZOp3njjTd46623OHbsWK46u8fjIRwOY7FYqK2tZe3atTQ0NKDX66WXlBDXwKxVMr799ts5deoUp06dmq1bCCHEvBYIBHjxxRfp7u5mYGAAo9GITqfD7XaTSCRwOBw0NDSwbt06KioqsFgsrFy5kpKSkrkeuhAL3qwFnAcffDC3bCWEEDebU6dO8ctf/pKDBw8yMjKCzWYjHA4TDofJZDKUl5fT2trKmjVrsNvtuFwuOjo6MJlMcz10IW4IsxZwvvCFL8zWpYUQYt7KZrN0d3fz61//moMHDxIOh8nPzycQCBAOh1EUhZqaGpYvX86KFSuwWq20tbWxaNGiuR66EDeUedtsUwghFppoNMobb7zB22+/TV9fHzqdDr1ej8/nIxaLYbFYqKuro7Ozk6amJhwOB6tWrZKDDkLMglkrqjAyMsJ//I//cbYuL4QQ88r4+DjPPfccL7zwAvv378dgMJBIJAiFQsRiMQoLC2ltbeXWW2+ltbWV+vp6Nm/eLOFGiFkyazM4Pp+Pn/70p/zoRz+arVsIIcScU1WVw4cP89prr7Fr1y6CwSBms5lQKEQqlSKTyeByuWhubmblypWUlZWxfPlyKisr53roQtzQZhxwXnrppcu+LqenhBA3umg0yu7du3n99dfp7e3NfT8cDpNMJtHpdCxevJjW1laWL19ORUUFnZ2d5OXlzeGohbg5zDjgPPDAAyiKwuXqBEoNByHEjcrtdrN9+3Zef/11BgcHMZvNxONx0uk0iUSCvLw8lixZQkdHB0uWLKGlpYXGxkb5d1GI62TGe3AqKir4xS9+gaqqF/3z3k8zQghxo1BVlb6+Pp5//nmeffZZBgcHMRgMhEIhkskk8XickpISOjo62LhxIx0dHWzatImlS5dKuBHiOprxDM6qVavo7e3lgQceuOjrHzS7I4QQC00sFmPv3r288cYb7Nq1C51ORzabJZlMks1myWQy1NbW0traSnt7Oy0tLSxbtmxBtK8R4kYz44Dzta99jWg0esnXGxoaeOutt2Z6eSGEmFfGx8fZs2cPv/rVrzhx4gQmk4lIJIKiKCSTScxmM42NjbS0tNDW1kZXVxdVVVVzPWwhblrXpNnmQiTNNoUQVyKbzXLkyBG6u7v513/9V0KhEJqmEYvFMBgMxONxioqKaGpqorm5mRUrVrB69WrZSCzELJnVZpuHDh2ivb0dne7KtvAcOXKEpqYmDAapKyiEWDjC4TA9PT3s2LGDt99+O1fbRtM0dDod8XicmpoampubaWlp4ZZbbpEmmULMEzNKHJ2dnUxMTFBWVnZF71+/fj0HDhxg8eLFM7mdEEJcd8PDw/T09PDKK69w4sQJzGYzfr8fo9FIJpNBr9fT1tZGQ0MDK1euZNOmTdjsDt495WMqnKDcZmFNfTF6nYQdIebCjAKOpmn81//6X694CjaVSs3kNkIIcd1lMhkOHTpEb28vL730ErFYDDhbvNRsNueqEre3t1NfX8+WLVtYsWIFrx2d4vGXe3AHE7lrVTgsPHZ/K3e1V8zV4whx05pRwNm8eTPHjx+/4vevX78eq9U6k1sJIcR1EwgE6O3tZceOHWzfvh2TyUQwGESv12M0GnNLUi0tLTQ3N3PPPffgdDp5tc/NI0/38v4NjRPBBI883ctTn1opIUeI62xGAeftt9++xsMQQoi5o2kag4ODHDx4kFdeeYWTJ09iNBqZnJwkPz+fVCqFTqejo6ODuro6brnlFjZv3ozZbCarajz+cv8F4QZAAxTg8Zf7uaPVJctVQlxHsutXCHFTSyQSHDhwgAMHDvDyyy+TTCZJp9OEQiHy8/OJxWIUFxfT3t5ObW0t999/P0uXLs39fPeQ77xlqffTAHcwQfeQj/VLSq7DEwkhYIYBZ3h4mNra2it+/9jYmNSDEELMO5OTk/T29rJ9+3beffddzGYzXq8Xo9GY22/T0NCQa7lw//33U1BQcN41psKXDjczeZ8Q4tqYUauGrq4u/viP/5ju7u5LvicYDPKP//iPtLe38/zzz894gEIIca2da7fw+uuv8/TTT7Nnzx4URWF8fJy8vDyy2Sw6nY6uri7a2tp48MEH+cQnPnFBuAEot1mu6J5X+j4hxLUxoxmcj370o9hsNu666y6MRiOrV6+msrISi8WC3++nv7+fI0eOsHr1ap588knuvvvuaz1uIYSYkXA4TG9vL/v27eO1115DVVUikQjpdBqbzUYsFsPpdNLa2kp9fT0f//jHKS8vv+T11tQXU+GwMBFMXHQfjgK4HGePjAshrp8ZVTI2mUyMjIxgt9txOp08/PDDeL1e4vE4paWldHZ2sm3bNtrb22djzNeEVDIW4uYzNDTEgQMHeOutt+jr60On0zE5OZk75amqKs3NzdTU1LBhwwbuuece9Hr9B1733Ckq4LyQc25LsZyiEuLaudLf3zMKOPX19Tz11FPcdddd6HQ6JiYmLvsJZz6SgCPEzSORSHDw4EEOHTrEr3/9awKBANFolHA4jM1mI5FIYLfbaW9vp6amhgceeIDm5uYPdY9X+9w8/nK/1MERYpbNaquGRx99lI985COsXr0aRVH453/+ZzZu3Eh7e7vUuxFCzCsTExP09PSwe/duuru7UVUVj8eDTqcjLy+PZDJJfX09DQ0NLFu2jAcffPCie20+yF3tFdzR6qJ7SCoZCzEfzLjZ5pEjR3jxxRf5L//lv7B48WJOnz6Noig0NDTQ0dHBihUr6OjomLf7b2QGR4iFL6tqlwwU2WyWvr4+Dh48yJtvvsno6CipVAqv14vVakVVVcxmM21tbdTU1LBt2zbWrl0rfaSEmOdmdYnqvRoaGnj33XfJz8/n0KFDuXoSBw4coK+vj3A4fDWXnzUScIRYWN4fZvzRJP/9345edEloXbWVnp4e9u7dy86dO0kkEvh8PlKpFBaLhUAkjq24jKVLm9nUtZyHPvYgLpdrDp9OCHGlrlvAuRxN0+btpyEJOEIsHBfb33JRmkrG7+bTS1WCQ4cZHBwkHo8zPT2NyWQilEgzGc6g2V0YilyYK1uoX76Wxx/skH0yQiwQ8yLgzGcScIRYGC7V5+n91FSc1OQgqcmTKJMn6HQaCQaDxGKxsz2l4inG40YMpdUYi6qw1q/EVFItJ52EWGBmdZOxEEJcD5fr83SOpmlkgpOkJgZJjh0l7RtDTccZSmXJNxswGM7+CRhsmKsqMJUvwVq/Ap05/+zPI/2ihLgRScARQsxbH9TnSU0nSU8NkZw4SXL8GJl4CDUehnSCjDkPRVEoKiqisKKOgUkjeTXtmJyLUXTn17aRflFC3Hgk4Agh5q3L9W/KhKdJTgySGjtK2jdCNhlHjYXQtCyKosNqtlBfv4iWlhZS9ir2jTkwFFy+mrD0ixLixiEBRwgxb12sf5OWTZOaOk1y8uySVDYWJBsLQTqJqmXR6fRYCuys61rJ8mXtrF69GrVkMb/42cEZ3U8IsTBJwBFCzFvv7/OUjfhJTAyQcp8gPT1MNhklGwtAVgU0dEYLhoIitm3q4rat61i3bh0NDQ2oGlQ4jku/KCFuIjPqJi6EENeDXqfw2P2taNkMqYmTxAa7iZ94l+TUEOmwl2zEj5bOgKKgz3NQWLGITz14D//p9z/KfffdR2NjI4qi5K4Dv+0Pdc65rx+7v1U2GAtxA5EZHCHEnLpcNWKAzjIdn10c5gc9+4iPn0ZNxciE/RgVlaJ8PXZbIQVFpTQ1t3DnprV0ruigra3tgiaZd7VX8NSnVl5QT8cl/aKEuCFJwBFCzJnLNai8ramUI0eO0NPTw3Dvu6wpjHHCnySUiUKBQp7ZgtPppLq6mq6uLpqamli5ciVOp/OS95N+UULcPOZlob8nnniC559/nmPHjmG1WtmwYQPf+c53aGpqyr0nmUzy6KOP8vOf/5x4PM5tt93G3//931NdXX1F95BCf0LMrUsV8FOAbCzIf2rRSE+c4NixY4TDYcbGxshmsyQSCRwOB06nk5aWFlatWkVjYyPLly/HbDbPxaMIIa6jK/39PS/34Gzfvp0vfOELvPvuu7z22mtkMhnuvPNOotFo7j1f/vKXeeGFF3j22Wd55513iEQi3HfffWSz2TkcuRDiSlyqgJ+WzZCcPEVsYA//40fPcORIP6OjowwPDxOPx0mn01RXV9PU1MS2bdu49dZb2bhxI11dXRJuhBDnmZczOO/n8XgoLy9n+/btbN68mWAwSFlZGT/72c/4+Mc/DsD4+Dg1NTW88sorbNu27QOvKTM4Qsyd3YNefv8f3z3ve9mon4T7JCn38bPViFMJakxRjIqWm7WpqqqiubmZrq4uqqqq6OzsJD8/f46eQggxF26oVg3BYBCA4uKzRzh7enpIp9PceeedufdUVlbS3t7Orl27LhpwkskkyWQy93UoFJrlUQshLuW9BfW0bJr09DBJ9wBJ93HUZIxsxEc2ESViVSiwGKmpqaG6upoNGzZQV1dHc3MzS5YsmbfNfIUQc2/eBxxN0/jKV77Cxo0baW9vB2BiYgKTyURRUdF573U6nUxMTFz0Ok888QSPP/74rI9XCPHBzhXUy0R8JCdOkho/Tto/jppKkA5No6hp1HQCW7mLpYvraG8/W7CvrKyMzs5OmXUVQnygeR9wvvjFL3Lo0CHeeeedD3yvpmmX/ET39a9/na985Su5r0OhEDU1NddsnEKIK9dRmY8tMsz4wJGzszapGJmIDzUeRcumURTIL6mkq7OdzZs3U1VVRWNjI42Njeh083LroBBinpnXAedLX/oSL730Ejt27DjvdJTL5SKVSuH3+8+bxZmammLDhg0XvZbZbJZNiELMAyMjI/T29tKa7Gdw+BBaKkE65IFsGi2TRG+1Y7CVcN/W9Tx01xaKioro7OyksLBwrocuhFhA5uVHIU3T+OIXv8jzzz/Pm2++SX19/Xmvr1q1CqPRyGuvvZb7ntvtpq+v75IBRwgxt2KxGO+++y4vvvgiL7zwAmpggkWWBFp4EjUVB03F4HBRWN3I5z/z+3zy/ttpamrilo2bOOrN8uKBMXYPesmq8/5chBBiHpiXMzhf+MIXeOaZZ3jxxRex2Wy5fTUOhwOr1YrD4eCzn/0sX/3qVykpKaG4uJhHH32UZcuWcfvtt8/x6IUQ76VpGkNDQ/T09LBnzx5GRkZydW0y8TjlFrCWOSksc7K8czV3beqiqNDBihUr6B5P8kff3X7RQoBSeVgIcTnz8pj4pfbR/PjHP+YP/uAPAEgkEnzta1/jmWeeOa/Q35Xuq5Fj4kLMvlAoxIEDB9i3bx/79+8nHo8zOjpKKBQiFothNJqwlTopq6xj/S2b6GyopLFhCU1NTbx2dOqShQABnvrUSgk5QtyErvT397wMONeDBBwhPpwP6hl13nuzWQYGBnKzNhMTE4RCIcbHx4nFYqTTacwFhXiwQXkThuJK9KY8ahpb+O8Pr+WOVhcbv/PmeTM373Wu+/c7//dWabMgxE3mhqqDI4SYW5frGfX+WRSPx8OBAwfo7u6mv7+fSCTCxMQEwWCQWCxGXl4eBaUVnMoWY6xsRW80YSyswFBShSel45Gne/ny7Y2XDDcAGuAOJuge8rF+SclsPbYQYgGTgCOEuKxL9YyaCCZ45One3FJRKpXiyJEj9Pb20t3djdfrxe/3MzU1RSQSAaCiooLFixezO+7EYrCjM1kxlS9GZynIXVcDfrzz9BWN7b0FA4UQ4r0k4AghLulSPaPgbBBRgMdf7qfZlubQwbOzNgMDA0SjUSYnJwkGg8Tjcex2OzU1NXR1dVFQtZQ33hjEWFSJoagCRbnwMGcgnr6i8Z0rGCiEEO8nAUcIcUndQ77LLhVlUwmGjhzlr3+0l8DwcYLBIF6vF6/XSzgcxmAwUF1dTUtLC7feeitFRUX8sj+ApboNnTnvsvd2WI2E4umLhqtze3DW1Bdf3QMKIW5YEnCEEJd0qSUgTc2S8btJTp4i5T7Oq6diOPRpkpEAwWCQTCZDUVERtbW1bNy4kdbWVgwGA1O6El73KejMH7wx+I6Wcn7RO4YC54Wccz/52P2tssFYCHFJEnCEEJd0sSWgbDRAauoUqYlBkp7TaMkYJ+MhtEQEJZvEWWSnvq6OlStXsmnTJsxmM+Xl5bS1L+OO/7H7ihtk3tJYxu2tzgs2N7ukDo4Q4gpIwBFCXNKa+mIqHBYmggnUTIqU5wypyUGSEyfJxgKoyRhqLEQ2EUXR6VAsdvzWSj72O/dze1cLJpOJ9vZ2qqqq2D3ovexy1/u57BbWLynhjlbXFR9PF0KIcyTgCCEuSa9T+Iv7Wvjjv/sVyclTJN0DZAJjqMkEmUQILR5By6bRm/LRFxRirm7D7Grg30+n+IMHa1jW3obRaAQ+3Imnivfsr9HrFDkKLoT40CTgCCEuyefzYZk+zt32EV7Y0006FCSTiKIlQqipBDqDEX1BCaayOsyLVqA3WdCZrESL6knaa3LhBj7ciSfZXyOEuFoScIQQF0gmk/T393P48GG6u7uZGhlhWZHKaCLGVNhPJJNFb7Wht5dhrVuOsagCFN15R7/fP2Pz3uWuS5VP1ynwt78vLRiEEFdPAo4QIudcY8zDhw/T09PD0aNHCYfD+P1+fD4fyWSSYruNeNKIpbIZc1Uzik6H3urAWLYInem3szTvn7HR6xQeu7+VR57uveBk1Dl/+/ud3LNcwo0Q4upJwBFCAOD1ejl8+DAHDx6kt7eX6elpQqFQrqaN1WrF5XLR0tLCO7FKIpoR9EaMpbUYbKW561yuRs1d7RU89amVV9z2QQghZkoCjhA3uUQiQX9/P319fezdu5czZ84QDofx+XwEg0EASktLqaysZOvWrSxdupTFw37+oTeMsbga9L/9Z+RKatTc1V4hJ6OEELNOAo4QNylVVTl16hRHjhyhp6eHI0eO5IKN3+8nnU5js9koLS1l1apVbNq0CYPBQFFREX+2ZQtrt0RnXKNGTkYJIWabBBwhbkKTk5McPnyYvr4+ent7mZqaIhQK4fP5iMX4+tXLAAAgAElEQVRiaHojVnsJNQ0NPPzReykrLcFkMtHS0kJNTQ2KonBXu11mYoQQ85YEHCFuIpFIhL6+Po4ePZpbjopEIkxPTxOJREhkNMKaCdVYhKW0nYChhtNvjfJ/fbSGz27bhMlkOu96MhMjhJivJOAIcRNIp9OcOHGC/v5+ent76e/vz52OCgaDqKqKYrIS0EyYnQ2Yq5tR9AZ0lgJixXX85d4E9U1e2QQshFgwJOAIsYBkVe1DLQlpmsbw8DBHjhzh8OHD7N+/P7fHxu/3k8lksFgsFBUVMZC0U7C0Fb3VhqI3YiypRW8ryfWOevzlfu5odckSlBBiQZCAI8QC8Wqf+0Mdr/Z4PPT393Ps2DG6u7sZGxsjGAzi9/uJx+OYTCbKysqoqamhcmkHJwey6BQFQ6ELY3EVik6fu5YGuIMJuod8siQlhFgQJOAIsQC82ufmkad7LyiONxFM8MjTvTz1qd9W/41EIvT39zMwMEB3dzcDAwOEw2G8Xi+xWAydTkdRUREVFRV0dnaycuVK9o9FMIx5MZbVoTNZLzmOD9NPSggh5pIEHCHmuayq8fjL/Ret/KtxtvbM4y/3c2tjCScHTnD8+HEOHTrE/v37cyejIpEImqZRUFBAeXk5zc3NrFu3jqKiIvLz81lf3co/u4c+cCwfpp+UEELMJQk4Qsxz3UO+85al3k/VVIbPnObvnp6GwCh79uzB6/Xi9/sJBAIAmEwmSktLWbx4MatXr6aurg6z2czSpUtZtGgRGgoVb7kv2SfqctWJhRBiPpKAI8Q8d7lloUzYS9o3SsY/zr+NT6GPnz0V5fP5yGazGI1GCgsLqa2tpb29nWXLlmE2m6mrq6Opqem8Y9+X6hN1JdWJhRBivpGAI8Q8d7FloWw8TNo7TCYwScp9nJRvnKlCHclIkFQqhclkoqioiKqqKpqamlixYgUOh4OysjLa2tqw2WwXXPNSfaKutDqxEELMJxJwhJjn1tQXU+GwMBFMkE0lSHuHSQcnSU2cJDV1GjUVQ5cME1Z1GI1GHA4HFRUVLF68mBUrVlBRUYHNZqO1tRWn03nZe0mfKCHEjUICjhDznF6n8PVtS3jk7/6NlN9NevoMyfETZBMR1HgE0nFKC/PIz8/D6XRSX19PW1sb9fX1FBQU0NTURG1tba6ezZXcT46CCyEWOgk4QsxjmUyGwcFBDBOD3Ffq5X/v+Q2xoBctHkFNxTDo9ZSXOFhU7aKuro6GhgZaWlqw2WwsWbKEJUuWYDDI/+ZCiJuP/MsnxDX2YasNX4yqqpw5c4ajx46zo7efPe/uIuSZpMYYYSIbAJMOq70QV3kxtbW1LFq0iNbWVkpLS6mpqaGpqQmLRY50CyFuXhJwhLiGPmy14ffTNI3x8XGOHTvGr7r7+acXXiU4OYKajKPFg+gUhfIiO64SB7W1tdTU1NDS0oLL5aKysjI3eyOEEDc7CThCXCMfptrwxUxNTXH06FHcbjfP/usb/Oo3+1CTUdRYAAUNDCYwmvHqClnR0MamtR3U1tbicrloaWmhuLiYrKqxe9ArG4SFEDc9CThCXANXWm34Ys0qvV4vx44dY2Jigp6eHvbu3cvuo8NkwkFQMygGM4rRgj6/EL29FFOhi+NU8WcdK2hv++3JqKudPRJCiBuJBBwhroEPqjZ8sWaVgUCAY8eO4Xa7OXjwIN3d3Xi9XsYmpogHQ2eDjcWOzmpDby/HaCvFUOjCYCslXlyNtXYZTmcpcPWzR0IIcaORgCPENXClTSinwglCoRDHjx9nfHyc/v5+du/ezdTUFD6fj2g0iqo3oVjt6C02dPYyTLYSDA4nensZxuIqDI5yFEWHJ5IErm72SAghblQScIS4Bq6kCaWaSuA/c5ztI1EGBgbYsWMHExMT+P1+wuEwVquVsrIyMFoIRq0Y7CUY7GUY7OUYS6oxOJwoOv0F95zJ7JEQQtzoJOAIcQ28t9rw+2dS1FSCjH8MuxomeMbGv+7YwejoKIFAgGAwiNVqxel0kp+f/382DVcQOB4nZi7CWFqHodB1XrB5f+PLDzN7JIQQNwsJOEJcA3qdckGzSjWdIOMbJxueJh3y0KAb4+fPTBIIBAiFQpjNZpxOJwUFBVRXV1NVVYXT6WTRokW03OriW7/xo+gNH9j48kpmjz7M+4QQ4kYgAUeIa+Rcs8q/+EUvo2dOkQlNkw570bn7KdX8eDIJgsEgRqOR8vJy8vPzqaqqoqamhtLSUurq6lixYgUNDQ2YTCZqFl94KupijS8vN3sEF874CCHEzUACjhDXSCwWoyI7xZ93qrxrzLBzx0GmJ4ZRUzFCoRBGoxGn04nVaqWyspKamhpKSkqora1l9erVLFmyBLPZnLvelTa+vNjs0TkXm/ERQoibgaJp2sU+9N3wQqEQDoeDYDCI3W6f6+GIBSwaPbtpeHR0lMnJSd58800GBgYIhUKEQiEMBgPFxcXk5+dTVlZGfX09RUVF1NbW0tXVRWNjIyaT6arHIXVwhBA3gyv9/S0zOELMUDgcZmBggPHxcSYnJ3njjTc4efJkLtjo9XpcLhdWq5Xy8nLq6+ux2+1UV1ezbt06li5dek2CzTlXOuMjhBA3Awk4QnxIoVCIEydO4Ha7GR8f580332RwcJBIJEIwGESv11NRUYHFYqGsrIzFixdjs9morq5m/fr1NDU1YTQaZ2Vsep0iR8GFEIJ5GnB27NjBk08+SU9PD263mxdeeIEHHngg97qmaTz++OP8wz/8A36/n7Vr1/J3f/d3tLW1zeGoxY3O5/Nx8uRJJicnGRkZ4Y033uD06dO5YGM0GqmqqsJsNlNSUkJDQ0PuhNTGjRtZunQpBsO8/F9OCCFuOPPyX9toNEpHRwd/+Id/yEMPPXTB63/5l3/JX/3VX/GTn/yEpUuX8s1vfpM77riD48ePSydlccWyqnZFyzlTU1OcPHmS6elpzpw5wxtvvMHw8DCRSIRQKITJZKK6uhqLxZILNvn5+dTU1LB582YaGxvR6XRz8IRCCHHzmpcB5+677+buu+++6GuapvH973+fP//zP+djH/sYAD/96U9xOp0888wz/Mmf/MlFfy6ZTJJMJnNfh0Khaz9wsWB80IZcTdNwu92cPHkSv9/P8ePH2b59O+Pj40Sj0QuCTWlpKUuXLsVisVBbW8uWLVtoaGhAUWT/ixBCzIV5GXAuZ2hoiImJCe68887c98xmM1u2bGHXrl2XDDhPPPEEjz/++PUappjHLteY8nM/28d/u91FtT5MKBTiwIED7Ny5E4/HQzgcJhKJYLVac8GmvLycxsZGzGYzS5YsYdOmTdTX10uwEUKIObbgAs7ExAQATqfzvO87nU7OnDlzyZ/7+te/zle+8pXc16FQiJqamtkZpJi3LtWYUstmSAcnyQYn+e8/zvBAVZx9+/bi8/kIh8NEo1EKCgqoq6vLVSA+V5CvubmZW2+9lcrKyjl5JiGEEBdacAHnnPd/QtY07bKfms1m83lF1MTN6f2NKdV0kkxggmzIQyYWIjl+DP/0aZ7P08gmY8RiMRwORy7YVFdXU19fj9FoZNmyZWzZsuVsg0whhBDzyoILOC6XCzg7k1NR8dviZVNTUxfM6gjxfucaTqrJGGn/ONmIj3TIQ2r8OGnfOFoyhpqO4Y2Bs7SYsrIyrFYrixYtora2FrPZTEdHB5s2baKoqGiOn0YIIcSlLLiAU19fj8vl4rXXXqOzsxOAVCrF9u3b+c53vjPHoxPzmaZpGFJhkmPHSEf9pH2jpN0nyYQ8aKkoajIOaGDKw1VRiqu0kCVLllBRUYHNZmPVqlVs2LCBvLy8uX4UIYQQH2BeBpxIJMLJkydzXw8NDXHgwAGKi4upra3ly1/+Mt/+9rdpbGyksbGRb3/72+Tl5fHJT35yDkctrrUrPcb9gdfJZjkzPMIrOw8wODyKdqaH6NgA2VgQNRlBTSVRdAqKOR/FYMZqs7FpQxcuZzmlpaWsWbOGVatWXdOqw0IIIWbXvAw4+/bt43d+53dyX5/bHPyZz3yGn/zkJ/zn//yficfjfP7zn88V+vv3f/93qYFzA7kWfZUSiQSnT5/mf7/Vw49/vQ/P8AnS3jGyiQhqMgKZFOj0KJYCdAYThnwHhkIX969pYvWKpaxbt47W1tZZqzoshBBi9kizTWm2Oe9c6hj3ubmbpz618rIhJxAIMDQ0xKlTp3hl536e+/UuMiEPmWQULRkFNQM6Axgt6AxmDLYSDA4nRcXFfGLrSv74oW0sXrwYvV4/a88ohBBiZqTZpliQLnWMG0DjbMh5/OV+7mh1nbdcpaoqbrebwcFBTp48yeDgIP39R3m99zjJaAQ1FQNNBb0JxWJDb7KgLyijoMzFvWtaaGlayifuuZWa6iqpYSOEEDcACThiXnn/Me730wB3MEH3kI/1S0pIJpOcOXOGo0ePcurUKY4fP87w8DA+nw+3x0fMH/g/wcaIYrKgNxegt5dicJRjKChBV1LD733iY9yzpun6PaQQQohZJwFHzCvnjnF/kFNjExh8p+jr62NoaIihoSHGx8cJBoOEw2GSySSxRBpFbwCTFb3VjsFehqGoEkNBMcbSGkxFlegsBaRNBbP8VEIIIa43CTjiql2r004A5TbLJV/T1CzZsJfU1BDdv97PnvAUIyMjeL1eQqEQsVgs128sLy8Pq72IcNSAweHCVFqNIb8YQ0k1BocTnfG3RR8vd08hhBALkwQccVWuxWmn91pTX0yFw8JEMJHbh6MmY6R9o6QmT5EJeTCmAhwNa4TDIeLxOLFYjEQigV6vx263Y7fbcTqdNDU38/+dyBDR287O3NjLUHS/3TisAC7H2UAmhBDixiKnqOQU1Yxd7Wmny133cz/bRybkITlxirRvlGzERzbqQ40GqLXryTMqhMNh0uk0JpMJu92OzWajtraWzs5OampqqKioYDhdwDdeHUVRlPPGebVjFEIIMTfkFJWYVTM97fRBQqEQVu8J7jUc5pf9h4n5fajhabKpGAZFpdyskI5F8WWz5OfnU1JSgsPhoKWlhc7OTsrKyqiqqqK+vp7CwkIAisucF8wyua5ilkkIIcT8JwFHzMiHPe10Odlslv7+fvbv38/g4CB+v58pt5sl2hQ+XZRkXpaUkiabSmJSzOTZbOTl5eFyuVi5ciXt7e3YbLbz+kW9113tFdzR6rpm+4SEEELMfxJwxIxc6Wmny71vbGyMnp4ejhw5gt/vx+PxMDo6SjweR9M0EokE8XCATCZzdtNwoZ38/HyWLl3K2rVrqaiooLi4ONefTKfTXfJeep3ygUFLCCHEjUMCjpiRKz15NDAZZvegNzdj4vf72b9/P4cOHcLj8RAIBBgbG8Pr9QJnZ3MCgQDxeBy9Xk9+fj5ms5mSkhJWrlzJypUrycvLo7q6mrq6OhwOx2w+phBCiAVKAo6YkYuddrqYv31rkL9+tQ9HYoKt5XEKtTChUIiJiQkmJydRVRVN00gmkwQCZ2drrFYrhYWFWK1WFi1axPr166mrq8Nut7No0SKqqqqkP5QQQojLkoAjZkSvU3js/lYeeboXBS4IOWoyRjrgJu05TSY4RTAZ50cxP00OlQLD2Xdrmobf7yeZTKIoCnl5ebkTUZ2dnXR0dFBUVERFRQV1dXWUlMgSkxBCiCsjAUfM2F3tFTz1qZW5E0pqIkImNEXKc4ZsyEM2GUONB8lGg6iJKOjgREShucxCPB7PHfG22WxYrVaqqqpYvXo1S5YsoaioiNraWqqrqy/YNCyEEEJ8EAk4YsZUVWW1y8hT95Tx8vZuuvtOsr3/FGoihBoLkk1EAQ1N01C0LGoqTSKTZtqYosRhw263k5eXR3NzM52dnbhcLiorK6mrq6O4WIrvCSGEmDkJOOJDSafTeDwe3G43J0+eZHJykunpaTweD9Onhkm6x0DNAhqoWVRFgWwaLZ0GvQHFnEee3cGiRdUsX76c5uZmnE4ndXV1srdGCCHENSMBR3ygWCzG5OQkk5OTnD59Gq/Xmws209PTBINBkskkiVgKMinQ6QAFLZtF0VQwGMBiQmcwYShycc8Dt3P7mnaqqqqoqamRk1BCCCGuOQk44gKqquL1epmammJqagqPx4PX68XtdjM1NYXP5yMUCpFIJFBVFVVVMRqNlBfZGfRESfyfhpcYTSiAzmLDVL4IU2kdruoavvzZj1Hhcl62bo0QQghxNSTgCADi8fh5gSYSieD1ehkdHcXj8eD3+wmHw7kifAAGgwGHw4GqqsTjcSKRMLXFeZyYzKAzGDEUujCWL8ZUUoPRUY6hoIQn/2AtVZXSHkEIIcTskoBzk8pkMrklpnOB5lzIGRkZwePxEAwGiUQixGIxdDodiqJgNBqx2+0UFBQQiUTw+Xy5WjY6nY7m+gpWrNnA3kgREXMpelsJOnPeVXUYF0IIIT4sCTg3CU3TCAQCeDye3IyMpmnEYjFGR0cZGRnJLT3FYrFcJWGdTkd+fj4FBQVUVFQQjUZz+2/OzeRYLBaWLl1KV1cXnZ2dVFdX4ygsYu9pv/R+EkIIMSck4NygNE0jFAoxPT2N1+vF6/WSyWTQNA2fz8fw8DBjY2MEAgESiQTRaJRUKoVer8dkMlFYWEhhYSG1tbVkMhmGh4c5dOgQ2WwWAJ1OR01NDcuXL2fLli3U19dTXl5+3r4a6f0khBBirkjAuUFomkY4HD4v0KTTaeDs0e7x8XFOnz7N1NRUboYmnU6TTCYxGAxYrVaKi4tzocZoNHLixAneffdd4vF4bomqsLCQlpYW7rjjDtra2nA6nRgM8p+REEKI+UV+My1QqqoSCATwer34fD78fn8u0GSzWXw+H0NDQ3g8Hnw+H+l0OhdoUqkU+fn5OBwOCgsLsdlsVFVVkZ+fz5EjR3jnnXcIBAK5JSqLxcLixYvZunUrmzZtwuVySb0aIYQQ85oEnAUilUrh9/vx+Xz4fD4CgQCqqgJnA43f72d8fJzJyUm8Xi/xeBxFUc7Wp0kkAHA4HFRUVOBwOMjLy8PlcmGxWDh69Ci/+c1vmJiYyIUas9lMZWUlGzdu5M4772TRokUSaoQQQiwYEnDmIVVVCYVC+P1+/H4/gUCAaDSaez2ZTDI5OUkgEMgtSSWTSTRNI5PJEI/HyWQy5OfnU11dTVFRERaLBYPBQHFxMVarlYGBAd58803Gx8dRFAWdTofRaKS4uJh169Zx//3309zcLMtPQgghFiT57TXHNE0jGo0SCAQIBoMEAoHzZmdUVcXv9+PxeIhGo4RCIQKBAJlMhlQqhaqqJJNJMpkMBoOB8vLy3CxNOp1Gp9NRVFSE0Wjk1KlTvPXWW4yPj6OqKnq9Hr1ej81mo7Ozk3vvvZd169ZJAT4hhBALngScWZZVNbqHfEyFE5QVmGkrNxMJh3KBJhgMkslkgLNhJhKJ4PF4CAQCxONxotEomUyGbDZLIpEglUqRzWZRFAWTyURBQQE1NTWUlJSgaRqRSASz2UxhYSGqqnL69Gl6e3uZmJjInZJSFAWHw8GyZcvYtm0bmzdvluUnIYQQNxQJOLMkk8nwwp4TfPvF/UxO+1FTMbRkjCKrgU+sqWFZRQHhcDi3nyaZTBKPxwFQFIVMJpM77aQoClarlfz8fAoLCykpKcHpdKIoCqlUKldd2Gq1YjAYOHPmDN3d3UxPT+dOQMHZPTjnTkDdeuutWK3WufwrEkIIIWaNBJxrrK+vj6mpKXb0j/DU24Oo2QxaJoWWjKMmI7hTcf7q6D5uXVrCopJ8LBYLcLZVQjweJ5FIkJeXh81mw263YzQaMRqNOByOXCAJh8NEo1Hy8vIwm81ks1lOnz6N2+1menqaWCyW21dTXFxMY2MjW7ZsYfPmzRQVFc3lX48QQghxXUjAucaGhoY4euwY//jKAeLhKBoaaBooCoreADoD6A10D/modliIx+MUFBRQWVmJ2WxGURQSiQQmkwmbzYbJZMoV4lNVFYPBgMlkwu/3MzIygtfrzYUaTdMwGo2UlpayZMkSNm/ezLp163A6nXP91yKEEEJcVxJwrrGBgQH2Hx0iHE+hmCwoCqApgIai6EBvQDEYSVkKKF7UQJXNQDgcRlEUzBYrkzENf1LFEE9RTThXYE9VVTweD7FYLHdUPBaLAWA2m3G5XNTW1rJ27VrWrVtHdXU1iiKtEYQQQtycJOBcY83NzfS7wxgcgM6AouhR9HowmCCbQUvFUJNxtHgYTyBCU0U1NpuNPSfGeGl7H6FIDBQFVJU8fYb1dTaqHWamp6cJBoNEo1EMBgNms5ma/7+9u41pq2z4AP5vCy0tT+kGDErDi8WxwWAoL5sZA7dExXuSGWN0zrmXZPqBPExBEgOKZhuPA90iWSLCUpN7X8wiH5w6d2NidQtIlt0QXnRhi8wMYfeAB6db25UBpb2eDws8D8JY99wd5/Tw/yX90Ku0/HslHP455zrnJCQgPj4e2dnZWLduHRITE3kGFBEREVhwAk6v1yPl4WRornghPOMQk2PwjrsAIaAK1UGt1UOjC4fwjMMQAly7dg2Xrt3EqZ+HIKY88Pm8gPBBeDy4Oe7EPwbHsCZWD/OyO+t1VqxYAYvFguzsbGRlZSEpKYmlhoiI6C9YcALs9u3b0Lh/h+H277jlVUOl1UOjN0JMjsM37saU83fAO4X/0Glg8EXA5RL4/sIgvBMeCO8kfLdvQUxNAFMeQK2BWhuGq5PhKMrJRU5ONjIzM1lqiIiI7oEFJ8AiIiKQbH0IL2xQ4+9nL0CMOeDzTAIQgFoNNVQQQiB/5QoInxdXRm7g5s0bwORtCK8HUKmhCgmFWm+CxhgNbWwytDEP4eldz2Bjygqpvx4REVFQYMEJsN7eXly+fBkatRpPrzSi7fIEXF41oAKg0sAYpkW+1YgYPTA0NIR/Df8JMeGGOlQHtSECIRGx0MY8hNDoRGjCl88sFL7unpT2ixEREQURFpwAMxgM0Ov1UKvVyI6ORl7Gw+gfdWH05k14x1wI9bgx6fxv/HE7FOHh4bAmPYRLKoHQ6ESERsVDo4+Y93NjjGGL/E2IiIiCFwtOgCUkJECj0cDtdmNkZAS//fYbHA4HvF4vDAYDDKYIREdHY/Xq1UhLS0PywyvxdP0/MewYX/Bzb7gnFukbEBERBT8WnADr7OzExYsX4Xa7odVqYTAYYLFYYLFYkJKSgtWrVyMxMRE6nW7mPe8VrcF/nuha8HP/6x+X8HRGHDRqXtuGiIjoXlhwAiwkJATh4eFISEhAQkICVq1ahVWrVmHFihV3PfNpebj2np877BhHe/+f2PBwVKAjExERKU5QF5yGhgYcOXIEw8PDSE9Px9GjR1FQUCBppqKiIoyNjcFqtcJoNPr1nlHXwoen7vfniIiIlrqgvZhKU1MTysrKUFVVhe7ubhQUFGDLli0YHByUNNfKlSuRmZnpd7kB/F9AzIXGRERE/gnaglNXV4dXX30Vr732GtLS0nD06FEkJCSgsbFR6mj3bb01EnGmMNxtdY0KQJwpDOutkYsZi4iIKGgFZcGZnJxEZ2cnCgsLZ40XFhbi3Llz875nYmICTqdz1kMuNGoV9m9dAwBzSs708/1b13CBMRERkZ+CsuBcv34dXq8XsbGxs8ZjY2MxMjIy73tqa2thMplmHgkJCYsR1W9/y4hD485smE2zD0OZTWFo3JmNv2XESZSMiIgo+AT1IuPpq/xOE0LMGZv29ttvo7y8fOa50+mUZcl5ao0Z7f1/YtQ1jhjjncNS3HNDRER0f4Ky4ERHR0Oj0czZWzM6Ojpnr840nU4369ozcqVRq3gqOBER0b8pKA9RabVa5OTkwG63zxq32+3Iy8uTKBURERHJRVDuwQGA8vJy7Nq1C7m5udiwYQNsNhsGBwdRXFwsdTQiIiKSWNAWnJdeegl//PEHqqurMTw8jIyMDDQ3NyMpKUnqaERERCQxlRBCSB1CCk6nEyaTCQ6HAxER89/Bm4iIiOTF3//fQbkGh4iIiGghLDhERESkOCw4REREpDgsOERERKQ4LDhERESkOEF7mvi/a/rkMTnddJOIiIgWNv1/+14ngS/ZguNyuQBAdvejIiIiontzuVwwmUx3fX3JXgfH5/NhaGgIRqPxrjfo/P+Yvonn1atXeX2de+Bc3R/Ol/84V/7jXPmPc+W/BzlXQgi4XC5YLBao1XdfabNk9+Co1WrEx8c/sM+PiIjgH4CfOFf3h/PlP86V/zhX/uNc+e9BzdVCe26mcZExERERKQ4LDhERESmO5sCBAwekDqE0Go0GmzdvRkjIkj0C6DfO1f3hfPmPc+U/zpX/OFf+k3quluwiYyIiIlIuHqIiIiIixWHBISIiIsVhwSEiIiLFYcEhIiIixWHBCbCGhgZYrVaEhYUhJycHP/74o9SRZKe2thbr1q2D0WhETEwMnnvuOfzyyy9SxwoKtbW1UKlUKCsrkzqKLF27dg07d+5EVFQUDAYDHn30UXR2dkodS3ampqbw7rvvwmq1Qq/XIzk5GdXV1fD5fFJHk4XW1lZs3boVFosFKpUKX3311azXhRA4cOAALBYL9Ho9Nm/ejN7eXonSSmuhufJ4PKioqMDatWsRHh4Oi8WC3bt3Y2hoaFGyseAEUFNTE8rKylBVVYXu7m4UFBRgy5YtGBwclDqarLS0tKCkpATnz5+H3W7H1NQUCgsL4Xa7pY4max0dHbDZbMjMzJQ6iizduHEDGzduRGhoKL799ltcvHgRH330EZYtWyZ1NNn58MMPcezYMdTX1+PSpUs4fPgwjhw5go8//ljqaLLgdrvxyCOPoL6+ft7XDx8+jLq6OtTX16OjowNmsxlPPfXUzD0Ol5KF5mpsbAxdXV1477330NXVhZMnT6Kvrw/PPvvs4oQTFDDr168XxcXFs8ZSU1NFZWWlRImCw+joqAAgWlpapI4iWy6XS6SkpAi73S42bdokSktLpY4kOxUVFYKeBjoAAAXUSURBVCI/P1/qGEGhqKhI7N27d9bY888/L3bu3ClRIvkCIL788suZ5z6fT5jNZvHBBx/MjI2PjwuTySSOHTsmRUTZ+Otczae9vV0AEAMDAw88D/fgBMjk5CQ6OztRWFg4a7ywsBDnzp2TKFVwcDgcAIDIyEiJk8hXSUkJioqK8OSTT0odRbZOnTqF3NxcvPjii4iJiUFWVhY+/fRTqWPJUn5+Pn744Qf09fUBAH766Se0tbXhmWeekTiZ/PX392NkZGTWtl6n02HTpk3c1vvB4XBApVItyp5VXooxQK5fvw6v14vY2NhZ47GxsRgZGZEolfwJIVBeXo78/HxkZGRIHUeWPv/8c3R1daGjo0PqKLJ25coVNDY2ory8HO+88w7a29vxxhtvQKfTYffu3VLHk5WKigo4HA6kpqZCo9HA6/Xi0KFDePnll6WOJnvT2/P5tvUDAwNSRAoa4+PjqKysxI4dOxblZqUsOAGmUqlmPRdCzBmj/7Vv3z78/PPPaGtrkzqKLF29ehWlpaX47rvvEBYWJnUcWfP5fMjNzUVNTQ0AICsrC729vWhsbGTB+YumpiZ89tlnOHHiBNLT09HT04OysjJYLBbs2bNH6nhBgdv6++PxeLB9+3b4fD40NDQsyu9kwQmQ6OhoaDSaOXtrRkdH5zR9uuP111/HqVOn0Nraivj4eKnjyFJnZydGR0eRk5MzM+b1etHa2or6+npMTExAo9FImFA+4uLisGbNmlljaWlp+OKLLyRKJF9vvfUWKisrsX37dgDA2rVrMTAwgNraWhacezCbzQDu7MmJi4ubGee2/u48Hg+2bduG/v5+nDlzZlH23gA8iypgtFotcnJyYLfbZ43b7Xbk5eVJlEqehBDYt28fTp48iTNnzsBqtUodSbaeeOIJXLhwAT09PTOP3NxcvPLKK+jp6WG5+T82btw453IDfX19SEpKkiiRfI2NjUGtnr3512g0PE3cD1arFWazeda2fnJyEi0tLdzWz2O63Fy+fBnff/89oqKiFu13cw9OAJWXl2PXrl3Izc3Fhg0bYLPZMDg4iOLiYqmjyUpJSQlOnDiBr7/+GkajcWavl8lkgl6vlzidvBiNxjlrk8LDwxEVFcU1S3/x5ptvIi8vDzU1Ndi2bRva29ths9lgs9mkjiY7W7duxaFDh5CYmIj09HR0d3ejrq4Oe/fulTqaLNy6dQu//vrrzPP+/n709PQgMjISiYmJKCsrQ01NDVJSUpCSkoKamhoYDAbs2LFDwtTSWGiuLBYLXnjhBXR1deH06dPwer0z2/vIyEhotdoHG+6Bn6e1xHzyySciKSlJaLVakZ2dzVOf5wFg3sfx48eljhYUeJr43X3zzTciIyND6HQ6kZqaKmw2m9SRZMnpdIrS0lKRmJgowsLCRHJysqiqqhITExNSR5OFs2fPzruN2rNnjxDizqni+/fvF2azWeh0OvH444+LCxcuSBtaIgvNVX9//12392fPnn3g2VRCCPFgKxQRERHR4uIaHCIiIlIcFhwiIiJSHBYcIiIiUhwWHCIiIlIcFhwiIiJSHBYcIiIiUhwWHCIiIlIcFhwiIiJSHBYcIiIiUhwWHCIiIlIcFhwiIiJSHBYcIlKMmpoaqFSqOY+6ujqpoxHRIuPNNolIMVwuF9xu98zz6upqNDc3o62tDfHx8RImI6LFFiJ1ACKiQDEajTAajQCAgwcPorm5GS0tLSw3REsQD1ERkeIcPHgQx48fR0tLC5KSkqSOQ0QSYMEhIkVhuSEigAWHiBSE5YaIpnENDhEpwvvvv4/6+nqcPn0aOp0OIyMjAIDly5dDp9NJnI6IFhvPoiKioCeEwLJly+B0Oue8dv78eTz22GMSpCIiKbHgEBERkeJwDQ4REREpDgsOERERKQ4LDhERESkOCw4REREpDgsOERERKQ4LDhERESkOCw4REREpDgsOERERKQ4LDhERESkOCw4REREpDgsOERERKc7/APcFRTZHNUlBAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Posterior distribution of w: 𝒩(xi=[2.71e+02, 1.71e+03, 1.20e+04], w=[[15.00, 72.34, 4.42e+02][72.34, 4.42e+02, 3.00e+03][4.42e+02, 3.00e+03, 2.20e+04]])\n", "\n" ] } ], "source": [ "# Build factorgraph\n", "\n", "fg = FactorGraph()\n", "@RV w ~ GaussianMeanVariance(constant(zeros(3)), constant(Ξ£, id=:Ξ£), id=:w) # p(w)\n", "for t=1:N\n", " x_t = Variable(id=:x_*t)\n", " d_t = Variable(id=:d_*t) # d=w'*x\n", " DotProduct(d_t, x_t, w) # p(f|w,x)\n", " @RV y_t ~ GaussianMeanVariance(d_t, constant(Οƒ2, id=:Οƒ2_*t), id=:y_*t) # p(y|d)\n", " placeholder(x_t, :x, index=t, dims=(3,))\n", " placeholder(y_t, :y, index=t);\n", "end\n", "\n", "# Build and run message passing algorithm\n", "algo = messagePassingAlgorithm(w)\n", "source_code = algorithmSourceCode(algo)\n", "eval(Meta.parse(source_code))\n", "data = Dict(:x => x_train, :y => y_train)\n", "w_posterior_dist = step!(data)[:w]\n", "\n", "# Plot result\n", "println(\"Posterior distribution of w: $(w_posterior_dist)\")\n", "scatter(z, y_train); xlabel(L\"z\"); ylabel(L\"f([1.0, z, z^2]) + \\epsilon\");\n", "z_test = collect(0:0.2:12)\n", "x_test = [[1.0; z; z^2] for z in z_test]\n", "for sample=1:10\n", " w = ForneyLab.sample(w_posterior_dist)\n", " f_est(x) = (w'*x)[1]\n", " plot(z_test, map(f_est, x_test), \"k-\", alpha=0.3);\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##
OPTIONAL SLIDES
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sum-Product Messages for the Equality Node\n", "\n", "- LetΒ΄s compute the SP messages for the **equality node** $f_=(x,y,z) = \\delta(z-x)\\delta(z-y)$: \n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$\\begin{align*}\n", "\\overrightarrow{\\mu}_{Z}(z) &= \\iint \\overrightarrow{\\mu}_{X}(x) \\overrightarrow{\\mu}_{Y}(y) \\,\\delta(z-x)\\delta(z-y) \\,\\mathrm{d}x \\mathrm{d}y \\\\\n", " &= \\overrightarrow{\\mu}_{X}(z) \\int \\overrightarrow{\\mu}_{Y}(y) \\,\\delta(z-y) \\,\\mathrm{d}y \\\\\n", " &= \\overrightarrow{\\mu}_{X}(z) \\overrightarrow{\\mu}_{Y}(z) \n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- By symmetry, this also implies (for the same equality node) that\n", "\n", "$$\\begin{align*}\n", "\\overleftarrow{\\mu}_{X}(x) &= \\overrightarrow{\\mu}_{Y}(x) \\overleftarrow{\\mu}_{Z}(x) \\quad \\text{and} \\\\\n", "\\overleftarrow{\\mu}_{Y}(y) &= \\overrightarrow{\\mu}_{X}(y) \\overleftarrow{\\mu}_{Z}(y)\\,.\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Let us now consider the case of Gaussian messages $\\overrightarrow{\\mu}_{X}(x) = \\mathcal{N}(x|\\overrightarrow{m}_X,\\overrightarrow{V}_X)$, $\\overrightarrow{\\mu}_{Y}(y) = \\mathcal{N}(y| \\overrightarrow{m}_Y,\\overrightarrow{V}_Y)$ and $\\overrightarrow{\\mu}_{Z}(z) = \\mathcal{N}(z|\\overrightarrow{m}_Z,\\overrightarrow{V}_Z)$. LetΒ΄s also define the precision matrices $\\overrightarrow{W}_X \\triangleq \\overrightarrow{V}_X^{-1}$ and similarly for $Y$ and $Z$. Then applying the SP update rule leads to multiplication of two Gaussian distributions, resulting in \n", "\n", "$$\\begin{align*}\n", "\\overrightarrow{W}_Z &= \\overrightarrow{W}_X + \\overrightarrow{W}_Y \\\\ \n", "\\overrightarrow{W}_Z \\overrightarrow{m}_z &= \\overrightarrow{W}_X \\overrightarrow{m}_X + \\overrightarrow{W}_Y \\overrightarrow{m}_Y\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- It follows that **message passing through an equality node is similar to applying Bayes rule**, i.e., fusion of two information sources. Does this make sense?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sum-Product Messages for the Addition Node\n", "\n", "- Next, consider an **addition node** $f_+(x,y,z) = \\delta(z-x-y)$: \n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$\\begin{align*}\n", "\\overrightarrow{\\mu}_{Z}(z) &= \\iint \\overrightarrow{\\mu}_{X}(x) \\overrightarrow{\\mu}_{Y}(y) \\,\\delta(z-x-y) \\,\\mathrm{d}x \\mathrm{d}y \\\\\n", " &= \\int \\overrightarrow{\\mu}_{X}(x) \\overrightarrow{\\mu}_{Y}(z-x) \\,\\mathrm{d}x \\,, \n", "\\end{align*}$$\n", "i.e., $\\overrightarrow{\\mu}_{Z}$ is the convolution of the messages $\\overrightarrow{\\mu}_{X}$ and $\\overrightarrow{\\mu}_{Y}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Of course, for Gaussian messages, these update rules evaluate to\n", "\n", "$$\\begin{align*}\n", "\\overrightarrow{m}_Z &= \\overrightarrow{m}_X + \\overrightarrow{m}_Y \\\\\n", "\\overrightarrow{V}_z &= \\overrightarrow{V}_X + \\overrightarrow{V}_Y \\,.\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "-
Exercise: For the same summation node, work out the SP update rule for the backward message $\\overleftarrow{\\mu}_{X}(x)$ as a function of $\\overrightarrow{\\mu}_{Y}(y)$ and $\\overleftarrow{\\mu}_{Z}(z)$. And further refine the answer for Gaussian messages.
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sum-Product Messages for Multiplication Nodes\n", "- Next, let us consider a **multiplication** by a fixed (invertible matrix) gain $f_A(x,y) = \\delta(y-Ax)$\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$\\begin{align*}\n", "\\overrightarrow{\\mu}_{Y}(y) = \\int \\overrightarrow{\\mu}_{X}(x) \\,\\delta(y-Ax) \\,\\mathrm{d}x = \\overrightarrow{\\mu}_{X}(A^{-1}y) \\,.\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For a Gaussian message input message $\\overrightarrow{\\mu}_{X}(x) = \\mathcal{N}(x|\\overrightarrow{m}_{X},\\overrightarrow{V}_{X})$, the output message is also Gaussian with \n", "$$\\begin{align*}\n", "\\overrightarrow{m}_{Y} = A\\overrightarrow{m}_{X} \\,,\\,\\text{and}\\,\\,\n", "\\overrightarrow{V}_{Y} = A\\overrightarrow{V}_{X}A^T\n", "\\end{align*}$$\n", "since \n", "$$\\begin{align*}\n", "\\overrightarrow{\\mu}_{Y}(y) &= \\overrightarrow{\\mu}_{X}(A^{-1}y) \\\\\n", " &\\propto \\exp \\left( -\\frac{1}{2} \\left( A^{-1}y - \\overrightarrow{m}_{X}\\right)^T \\overrightarrow{V}_{X}^{-1} \\left( A^{-1}y - \\overrightarrow{m}_{X}\\right)\\right) \\\\\n", " &= \\exp \\big( -\\frac{1}{2} \\left( y - A\\overrightarrow{m}_{X}\\right)^T \\underbrace{A^{-T}\\overrightarrow{V}_{X}^{-1} A^{-1}}_{(A \\overrightarrow{V}_{X} A^T)^{-1}} \\left( y - A\\overrightarrow{m}_{X}\\right)\\big) \\\\\n", " &\\propto \\mathcal{N}(y| A\\overrightarrow{m}_{X},A\\overrightarrow{V}_{X}A^T) \\,.\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "-
Excercise: Proof that, for the same factor $\\delta(y-Ax)$ and Gaussian messages, the (backward) sum-product message $\\overleftarrow{\\mu}_{X}$ is given by \n", "$$\\begin{align*}\n", "\\overleftarrow{\\xi}_{X} &= A^T\\overleftarrow{\\xi}_{Y} \\\\\n", "\\overleftarrow{W}_{X} &= A^T\\overleftarrow{W}_{Y}A\n", "\\end{align*}$$\n", "where $\\overleftarrow{\\xi}_X \\triangleq \\overleftarrow{W}_X \\overleftarrow{m}_X$ and $\\overleftarrow{W}_{X} \\triangleq \\overleftarrow{V}_{X}^{-1}$ (and similarly for $Y$).
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "### Code example: Gaussian forward and backward messages for the Addition node\n", "\n", "Let's calculate the Gaussian forward and backward messages for the addition node in ForneyLab. \n", "

" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forward message on Z: 𝒩(m=3.00, v=2.00)\n" ] } ], "source": [ "# Forward message towards Z\n", "fg = FactorGraph()\n", "@RV x ~ GaussianMeanVariance(constant(1.0), constant(1.0), id=:x) \n", "@RV y ~ GaussianMeanVariance(constant(2.0), constant(1.0), id=:y)\n", "@RV z = x + y; z.id = :z\n", "\n", "q = PosteriorFactorization(fg)\n", "\n", "algo1 = messagePassingAlgorithm(z, id=:_forward_Z)\n", "source_code1 = algorithmSourceCode(algo1)\n", "eval(Meta.parse(source_code1))\n", "msg_forward_Z = step_forward_Z!(Dict())[:z]\n", "print(\"Forward message on Z: $(msg_forward_Z)\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Backward message on X: 𝒩(m=1.00, v=2.00)\n" ] } ], "source": [ "# Backward message towards X\n", "fg = FactorGraph()\n", "@RV x; x.id=:x\n", "@RV y ~ GaussianMeanVariance(constant(2.0), constant(1.0), id=:y)\n", "@RV z = x + y\n", "GaussianMeanVariance(z, constant(3.0), constant(1.0), id=:z) \n", "\n", "q = PosteriorFactorization(fg)\n", "\n", "algo2 = messagePassingAlgorithm(x, id=:_backward_X)\n", "source_code2 = algorithmSourceCode(algo2)\n", "eval(Meta.parse(source_code2))\n", "msg_backward_X = step_backward_X!(Dict())[:x]\n", "print(\"Backward message on X: $(msg_backward_X)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code Example: forward and backward messages for the Matrix Multiplication node\n", "\n", "In the same way we can also investigate the forward and backward messages for the matrix multiplication (\"gain\") node \n", "

" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forward message on Y: 𝒩(m=4.00, v=16.00)\n" ] } ], "source": [ "# Forward message towards Y\n", "fg = FactorGraph()\n", "@RV x ~ GaussianMeanVariance(1.0, 1.0)\n", "@RV y = 4.0 * x\n", "\n", "q = PosteriorFactorization(fg)\n", "#This is where the bugs live..\n", "algo3 = messagePassingAlgorithm(y, id=:_y_fwd)\n", "source_code3 = algorithmSourceCode(algo3)\n", "eval(Meta.parse(source_code3))\n", "msg_forward_Y = step_y_fwd!(Dict())[:y]\n", "print(\"Forward message on Y: $(msg_forward_Y)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " \n", " ### Example: Sum-Product Algorithm to infer a posterior\n", " \n", " - Consider a generative model \n", "$$p(x,y_1,y_2) = p(x)\\,p(y_1|x)\\,p(y_2|x) .$$ \n", " - This model expresses the assumption that $Y_1$ and $Y_2$ are independent measurements of $X$.\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ " \n", "- Assume that we are interested in the posterior for $X$ after observing $Y_1= \\hat y_1$ and $Y_2= \\hat y_2$. The posterior for $X$ can be inferred by applying the sum-product algorithm to the following graph:\n", "\n", "

" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " - (Note that) we usually draw terminal nodes for observed variables in the graph by smaller solid-black squares. This is just to help the visualization of the graph, since the computational rules are no different than for other nodes. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code for Sum-Product Algorithm to infer a posterior\n", "\n", "We'll use ForneyLab to build the above graph, and perform sum-product message passing to infer the posterior $p(x|y_1,y_2)$. We assume $p(y_1|x)$ and $p(y_2|x)$ to be Gaussian likelihoods with known variances:\n", "$$\\begin{align*}\n", " p(y_1\\,|\\,x) &= \\mathcal{N}(y_1\\,|\\,x, v_{y1}) \\\\\n", " p(y_2\\,|\\,x) &= \\mathcal{N}(y_2\\,|\\,x, v_{y2})\n", "\\end{align*}$$\n", "Under this model, the posterior is given by:\n", "$$\\begin{align*}\n", " p(x\\,|\\,y_1,y_2) &\\propto \\overbrace{p(y_1\\,|\\,x)\\,p(y_2\\,|\\,x)}^{\\text{likelihood}}\\,\\overbrace{p(x)}^{\\text{prior}} \\\\\n", " &=\\mathcal{N}(x\\,|\\,\\hat{y}_1, v_{y1})\\, \\mathcal{N}(x\\,|\\,\\hat{y}_2, v_{y2}) \\, \\mathcal{N}(x\\,|\\,m_x, v_x) \n", "\\end{align*}$$\n", "so we can validate the answer by solving the Gaussian multiplication manually." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sum-product message passing result: p(x|y1,y2) = 𝒩(1.1428571428571428,0.5714285714285714)\n", "Manual result: p(x|y1,y2) = 𝒩(1.1428571428571428, 0.5714285714285714)\n" ] } ], "source": [ "using ForneyLab \n", "\n", "# Data\n", "y1_hat = 1.0\n", "y2_hat = 2.0\n", "\n", "# Construct the factor graph\n", "fg = FactorGraph()\n", "@RV x ~ GaussianMeanVariance(constant(0.0), constant(4.0), id=:x) # Node p(x)\n", "@RV y1 ~ GaussianMeanVariance(x, constant(1.0)) # Node p(y1|x)\n", "@RV y2 ~ GaussianMeanVariance(x, constant(2.0)) # Node p(y2|x)\n", "Clamp(y1, y1_hat) # Terminal (clamp) node for y1\n", "Clamp(y2, y2_hat) # Terminal (clamp) node for y2\n", "# draw(fg) # draw the constructed factor graph\n", "\n", "# Perform sum-product message passing\n", "algo4 = messagePassingAlgorithm(x, id=:_x)\n", "source_code4 = algorithmSourceCode(algo4)\n", "eval(Meta.parse(source_code4))\n", "x_marginal = step_x!(Dict())[:x]\n", "println(\"Sum-product message passing result: p(x|y1,y2) = 𝒩($(mean(x_marginal)),$(var(x_marginal)))\")\n", "\n", "# Calculate mean and variance of p(x|y1,y2) manually by multiplying 3 Gaussians (see lesson 4 for details)\n", "v = 1 / (1/4 + 1/1 + 1/2)\n", "m = v * (0/4 + y1_hat/1.0 + y2_hat/2.0)\n", "println(\"Manual result: p(x|y1,y2) = 𝒩($(m), $(v))\")" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Backward message on X: 𝒩(xi=8.00, w=16.00)\n" ] } ], "source": [ "# Backward message towards X\n", "fg = FactorGraph()\n", "x = Variable(id=:x)\n", "@RV y = constant(4.0) * x\n", "GaussianMeanVariance(y, constant(2.0), constant(1.0))\n", "\n", "q = PosteriorFactorization(fg)\n", "\n", "algo5 = messagePassingAlgorithm(x, id=:_x_fwd2)\n", "source_code5 = algorithmSourceCode(algo5)\n", "eval(Meta.parse(source_code5))\n", "msg_backward_X = step_x_fwd2!(Dict())[:x]\n", "print(\"Backward message on X: $(msg_backward_X)\")" ] }, { "cell_type": "code", "execution_count": 16, "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 display(\"text/html\", read(f, String)) end\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.5.2", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.2" } }, "nbformat": 4, "nbformat_minor": 4 }