{ "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 algorithms\n", "- Materials \n", " - Mandatory\n", " - These lecture notes \n", " - [Loeliger, 2007](./files/Loeliger-2007-The-factor-graph-approach-to-model-based-signal-processing.pdf), pp. 1295-1300 (until section IV)\n", " - Optional\n", " - [Video lecture](https://www.youtube.com/watch?v=Fv2YbVg9Frc&t=31) by Frederico Wadehn (ETH Zurich) (**highly recommended**)\n", " \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 generative model $p(x_1,x_2,x_3,x_4,x_5)$, the inference task $p(x_2|x_3)$ is given by \n", "\n", "\\begin{align*}\n", "p(x_2|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 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 generative distribution 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." ] }, { "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](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=910573&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F18%2F19638%2F00910573.pdf%3Farnumber%3D910573))\n", "\n", " 1. A **node** for every factor;\n", " 1. An **edge** (or **half-edge**) for every variable;\n", " 1. Node $g$ is connected to edge $x$ **iff** variable $x$ appears in factor $g$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Some FFG Terminology\n", "\n", "- $f$ is called the **global function** and $f_\\bullet$ are the **factors**. \n", "\n", "- A **configuration** is an assigment of values to all variables.\n", "\n", "- The **configuration space** is the set of all configurations, i.e., the domain of $f$\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}$ 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", " &\\triangleq g(x_1 \\mid x_2)\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "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", "" ] }, { "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 (i.e., by message passing). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume we wish to compute the marginal\n", "$$\n", "\\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)\n", "$$\n", "where $f$ is factorized as given by the following FFG\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Due to the factorization, we can decompose this sum by the **distributive law** as\n", "\\begin{align*}\n", "\\bar{f}(x_3) = & \\underbrace{ \\left( \\sum_{x_1,x_2} f_a(x_1)\\,f_b(x_2)\\,f_c(x_1,x_2,x_3)\\right) }_{\\overrightarrow{\\mu}_{X_3}(x_3)} \\\\\n", " & \\underbrace{ \\cdot\\left( \\sum_{x_4,x_5} f_d(x_4)\\,f_e(x_3,x_4,x_5) \\cdot \\underbrace{ \\left( \\sum_{x_6,x_7} f_f(x_5,x_6,x_7)\\,f_g(x_7)\\right) }_{\\overleftarrow{\\mu}_{X_5}(x_5)} \\right) }_{\\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": "subslide" } }, "source": [ "- Messages may flow in both directions on any edge (here on $X_3$). We often draw _directed edges_ in a 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). With directed edges, the FFG looks as follows: \n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Crucially, drawing arrows on edges is only meant as a notational convenience. Technically, an FFG is an undirected graph. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that $\\overleftarrow{\\mu}_{X_5}(x_5)$ is obtained by multiplying all enclosed factors ($f_f$, $f_g$) by the green dashed box, followed by marginalization over all enclosed variables ($x_6$, $x_7$). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This is the **Closing the Box**-rule, which is a general recipe for marginalization of hidden variables and leads to a new factor with outgoing (sum-product) message \n", "$$\\mu_{\\text{SP}} = \\sum_{ \\stackrel{ \\textrm{enclosed} }{ \\textrm{variables} } } \\;\\prod_{\\stackrel{ \\textrm{enclosed} }{ \\textrm{factors} }}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Evaluating the closing-the-box rule for individual nodes\n", " \n", "\n", "- First, closing a box around the **terminal nodes** leads to $\\overrightarrow{\\mu}_{X_1}(x_1) \\triangleq f_a(x_1)$, $\\overrightarrow{\\mu}_{X_2}(x_2) \\triangleq f_b(x_2)$ etc. \n", " - So, the message out of a terminal node is the factor itself.\n", " - (Exercise) Derive now that the message coming from the open end of a half-edge always equals $1$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The messages from **internal nodes** evaluate to:\n", "\\begin{align*}\n", "\\overrightarrow{\\mu}_{X_3}(x_3) &= \\sum_{x_1,x_2} f_a(x_1) \\,f_b(x_2) \\,f_c(x_1,x_2,x_3) \\\\\n", " &= \\sum_{x_1,x_2} \\overrightarrow{\\mu}_{X_1}(x_1) \\overrightarrow{\\mu}_{X_2}(x_2) \\,f_c(x_1,x_2,x_3) \\\\\n", "\\overleftarrow{\\mu}_{X_5}(x_5) &= \\sum_{x_6,x_7} f_f(x_5,x_6,x_7)\\,f_g(x_7) \\\\\n", " &= \\sum_{x_6,x_7} \\overrightarrow{\\mu}_{X_7}(x_7)\\, f_f(x_5,x_6,x_7) \\\\\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Crucially, all message update rules can be computed from information that is **locally available** at each node." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Sum-Product Algorithm\n", "\n", "- (**Sum-Product update rule**). This recursive pattern for computing messages applies generally and is called the **Sum-Product update rule**, which is really just a special case of the closing-the-box rule: For any node, the outgoing message is obtained by taking the product of all incoming messages and the node function, followed by summing out (marginalization) all incoming variables. What is left (the outgoing message) is a function of the outgoing variable only: \n", "\n", "$$\\boxed{\n", "\\overrightarrow{\\mu}_{Y}(y) = \\sum_{x_1,\\ldots,x_n} \\overrightarrow{\\mu}_{X_1}(x_1)\\cdots \\overrightarrow{\\mu}_{X_n}(x_n) \\,f(y,x_1,\\ldots,x_n) }\n", "$$\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- (**Sum-Product Theorem**). If the factor graph for a function $f$ has **no cycles**, 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 the Sum-Product Theorem:\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": "subslide" } }, "source": [ "- (**Sum-Product Algorithm**). If folows 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": "fragment" } }, "source": [ "- (Exercise) Verfiy for yourself that all maginals in a cycle-free graph (a tree) can be computed exactly by starting with messages at the terminals and working towards the root of the tree." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Processing Observations in a Factor Graph\n", "\n", " - Terminal nodes can be used describe **observed variables**, e.g., use a factor $$f_Y(y) = \\delta(y-3)$$ to terminate the edge for variable $Y$ if $y=3$ is observed.\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " \n", " ##### Example\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": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "source": [ "- (Exercise) Can you draw the messages that infer $p(x\\,|\\,y_1,y_2)$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### CODE EXAMPLE\n", "\n", "We'll use ForneyLab, a factor graph toolbox for Julia, 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": 5, "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", "eval(Meta.parse(sumProductAlgorithm(x, name=\"_algo1\"))) # Automatically derives a message passing schedule\n", "x_marginal = step_algo1!(Dict())[:x] # Execute algorithm and collect marginal distribution of 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": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example: SP 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}(\\overrightarrow{m}_X,\\overrightarrow{V}_X)$, $\\overrightarrow{\\mu}_{Y}(y) = \\mathcal{N}(\\overrightarrow{m}_Y,\\overrightarrow{V}_Y)$ and $\\overrightarrow{\\mu}_{Z}(z) = \\mathcal{N}(\\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": [ "### (OPTIONAL SLIDE) Example: SP Messages for the Addition Nodes\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 \\,,\\,\\text{and}\\,\\,\\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": [ "### (OPTIONAL SLIDE) Example: SP 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}(\\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 \\left( -\\frac{1}{2} \\left( y - A\\overrightarrow{m}_{X}\\right)^T A^{-T}\\overrightarrow{V}_{X}^{-1} A \\left( y - A\\overrightarrow{m}_{X}\\right)\\right) \\\\\n", " &\\propto \\mathcal{N}(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", "### (OPTIONAL SLIDE) \n", "#### CODE EXAMPLE\n", "\n", "Let's calculate the Gaussian forward and backward messages for the addition node in ForneyLab. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forward message on Z: 𝒩(m=3.00, v=2.00)\n", "Backward message on X: 𝒩(m=1.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", "eval(Meta.parse(sumProductAlgorithm(z, name=\"_z_fwd\")))\n", "msg_forward_Z = step_z_fwd!(Dict())[:z]\n", "print(\"Forward message on Z: $(msg_forward_Z)\")\n", "\n", "# Backward message towards X\n", "fg = FactorGraph()\n", "@RV x = Variable(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", "eval(Meta.parse(sumProductAlgorithm(x, name=\"_x_bwd\")))\n", "msg_backward_X = step_x_bwd!(Dict())[:x]\n", "print(\"Backward message on X:$(msg_backward_X)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### (OPTIONAL SLIDE) \n", "#### CODE EXAMPLE\n", "\n", "In the same way we can also investigate the forward and backward messages for the gain node " ] }, { "cell_type": "code", "execution_count": 3, "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(constant(1.0), constant(1.0), id=:x)\n", "@RV y = constant(4.0) * x; y.id = :y\n", "\n", "eval(Meta.parse(sumProductAlgorithm(y, name=\"_y_fwd\")))\n", "msg_forward_Y = step_y_fwd!(Dict())[:y]\n", "print(\"Forward message on Y: $(msg_forward_Y)\")" ] }, { "cell_type": "code", "execution_count": 4, "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", "eval(Meta.parse(sumProductAlgorithm(x, name=\"_x_fwd2\")))\n", "msg_backward_X = step_x_fwd2!(Dict())[:x]\n", "print(\"Backward message on X:$(msg_backward_X)\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Example: Bayesian Linear Regression\n", "\n", "- Recall: the goal of regression is to estimate an unknown function from a set of (noisy) function values." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Assume we want to estimate some function $f: \\mathbb{R}^D \\rightarrow \\mathbb{R}$ from data set $D = \\{(x_1,y_1), \\ldots, (x_N,y_N)\\}$, where $y_i = f(x_i) + \\epsilon_i$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- 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(D,w) &= \\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": [ "- We are interested in inferring the posterior $p(w|D)$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Here's the factor graph for this model\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### CODE EXAMPLE\n", "\n", "Let's build the factor graph in Julia (with the FFG toolbox ForneyLab)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAG0CAYAAAA2BP2yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3X9wVPW9//HXJkAC6WY1gWQ3l1yMOLdtCFD5pQGKiASDGKW2KlYsVmUqIkqR1mJrY/yK0eut6L3epkr9HRBnqhFoIbe0SihFBBKjhHQsxdwaZdNUwE2gZIFkv39wd8uShCT7I+fs2edjZmc8Zz85+8bMsC8+P20+n88nAAAAC0swugAAAIBoI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLG2B0AUbp6OjQwYMHZbfbZbPZjC4HAAD0gs/nU2trq7KyspSQ0Pt+m7gNPAcPHlR2drbRZQAAgBA0NjZq+PDhvW4ft4HHbrdLOv0/LDU11eBqAABAb7S0tCg7OzvwPd5bcRt4/MNYqampBB4AAGJMX6ejMGkZAABYHoEHAABYHoEHAABYHoEHAABYHoEHAABYHoEHAABYHoEHAABYnukDT2lpqWw2m5YuXRq45/V6tWTJEg0dOlQpKSm65ppr9OmnnxpYJQAAMDNTB57du3frueee05gxY4LuL126VBUVFVq3bp22b9+uo0eP6uqrr1Z7e7tBlQIAADMzbeA5evSobr75Zq1evVrnn39+4L7H49Hzzz+vn/3sZ5o5c6YuvvhilZeXa+/evfrd735nYMUAAMSX9g6f3j1wSOtrP9O7Bw6pvcNndEndMu3REosXL9acOXM0c+ZMPfLII4H71dXVOnnypGbNmhW4l5WVpby8PO3YsUNXXnlll8/zer3yer2B65aWlugVDwCAxVXWuVWysV5uT1vgnsuRrOKiXBXmuQysrGum7OFZt26dampqVFpa2um9pqYmDRo0KKjXR5IyMzPV1NTU7TNLS0vlcDgCL05KBwAgNJV1bi0qrwkKO5LU5GnTovIaVda5Daqse6YLPI2Njbr33ntVXl6u5OTkXv+cz+c750FiK1askMfjCbwaGxsjUS4AAHGlvcOnko316mrwyn+vZGO96Ya3TBd4qqur1dzcrPHjx2vAgAEaMGCAqqqq9J//+Z8aMGCAMjMzdeLECR05ciTo55qbm5WZmdntc5OSkgIno3NCOgAAodnVcLhTz86ZfJLcnjat2vJnU83rMV3gueKKK7R3717V1tYGXhMmTNDNN98c+O+BAwdqy5YtgZ9xu92qq6vT5MmTDawcAADra27tPuyc6Zl3/qKbVu/U1MffNsUQl+kmLdvtduXl5QXdS0lJUXp6euD+7bffrvvuu0/p6elKS0vT8uXLNXr0aM2cOdOIkgEAiBsZ9t5PN5H+Oa+nbP44Qyczmy7w9MaqVas0YMAA3XDDDTp+/LiuuOIKvfTSS0pMTDS6NAAA+kV7h0+7Gg6rubVNGfZkTcpJU2JC93NZI2VSTppcjmQ1edq6nMdzNp8km07P6ynIdfZLjV2x+Xw+cwyu9bOWlhY5HA55PB7m8wAAYorRS8L9q7Qk9Sr0+L228FLlj0wP67ND/f423RweAADQPTMsCS/Mc6ls/jg5HX0b3urt/J9oIPAAABAjzLQkvDDPpe33z9BrCy/V3ZeP7NXP9HX+TyQReAAAiBG9XRK+q+Fwv9STmGBT/sh0fb/gy3I5ktXd7BybTg+5TcpJ65e6ukLgAQAgRvR2SKi/h44SE2wqLsqVpE6hx39dXJRr2IRlicADAEDM6O2QkBFDR93N63E6kg1fki7F6LJ0AADiUU9Lwm06HTCMGjoqzHOpINdpyHL5nhB4AACIEf6ho0XlNbIpeEm4WYaO/PN6zIYhLQAAYojZh47Mih4eAABijJmHjsyKwAMAQAwy69CRWTGkBQAALI/AAwAALI/AAwAALI/AAwAALI/AAwAALI/AAwAALI/AAwAALI/AAwAALI+NBwEAMLn2Dh+7KoeJwAMAgIlV1rlVsrFebk9b4J7LkaziolzOzeoDhrQAADCpyjq3FpXXBIUdSWrytGlReY0q69wGVRZ7CDwAAJhQe4dPJRvr5eviPf+9ko31au/oqgXORuABAMCEdjUc7tSzcyafJLenTbsaDvdfUTGMwAMAgAk1t3YfdkJpF+8IPAAAmFCGPTmi7eIdgQcAABOalJMmlyNZ3S0+t+n0aq1JOWn9WVbMIvAAAGBCiQk2FRflSlKn0OO/Li7KZT+eXiLwAABgUoV5LpXNHyenI3jYyulIVtn8cezD0wdsPAgAgIkV5rlUkOtkp+UwEXgAADC5xASb8kemG11GTGNICwAAWB49PAAAmAgHhUYHgQcAAJPY9OFB/WR9nQ4fOxm4x0GhkcGQFgAAJlC6qV53rX0/KOxIp4+P4KDQ8Jky8JSVlWnMmDFKTU1Vamqq8vPztXnz5sD706dPl81mC3rNmzfPwIoBAAjdpg/denZbQ7fv+8RBoeEyZeAZPny4HnvsMe3Zs0d79uzRjBkzdO2112rfvn2BNgsXLpTb7Q68nn32WQMrBgAgNO0dPv1kfV2P7TgoNDymnMNTVFQUdL1y5UqVlZVp586dGjVqlCRpyJAhcjqdRpQHAEDE7Go4rMPHTvSqLQeFhs6UPTxnam9v17p163Ts2DHl5+cH7q9Zs0ZDhw7VqFGjtHz5crW2tp7zOV6vVy0tLUEvAACM1pcQw0GhoTNlD48k7d27V/n5+Wpra9OXvvQlVVRUKDf39JkiN998s3JycuR0OlVXV6cVK1bogw8+0JYtW7p9XmlpqUpKSvqrfAAAeqW3ISYtZSAHhYbB5vP5TDkD6sSJE/rkk0/0xRdf6I033tAvf/lLVVVVBULPmaqrqzVhwgRVV1dr3LhxXT7P6/XK6/UGrltaWpSdnS2Px6PU1NSo/TkAADiX9g6fpj7+ttyec/f0/PzbF+uqMVn9VJV5tbS0yOFw9Pn727RDWoMGDdJFF12kCRMmqLS0VGPHjtXTTz/dZdtx48Zp4MCB2r9/f7fPS0pKCqz68r8AADCa/1T0c20t+L1pOYSdMJk28JzN5/MF9dCcad++fTp58qRcLjZlAgDEHv+p6K6zTkVPTxmkn397nFZc1Xl0A31jyjk8DzzwgGbPnq3s7Gy1trZq3bp12rp1qyorK3XgwAGtWbNGV111lYYOHar6+nrdd999uvjiizVlyhSjSwcAoEddHR/BqejRZcrA87e//U233HKL3G63HA6HxowZo8rKShUUFKixsVG///3v9fTTT+vo0aPKzs7WnDlzVFxcrMTERKNLBwDgnCrr3CrZWB80Z+fM4yM4FT06TDtpOdpCnfQEAECoKuvcWlReo7O/eP19OGXzx3FmVg8sN2kZAAArae/wqWRjfaewIylwj+MjoofAAwBAP9jVcPicS8994viIaCLwAADQD3q7o3KT53iUK4lPBB4AAPpBb3dU/n+/+ZMq69xRrib+EHgAAOgHk3LS5HIkn3ODQUk6cuyEFpXXEHoijMADAEA/8O+o3BMmMEcHgQcAgH7i31E5LWXgOdsxgTnyCDwAAPSjwjyXHrx6VK/a9naiM3pG4AEAoJ85U3s3gbm3E53RMwIPAAD9rKcJzDadPm5iUk5af5ZlaQQeAAD62ZkTmM8OPf7r4qJcDg6NIAIPAAAG8E9gdjqCh62cjmTO1IoCU56WDgBAPCjMc6kg16ldDYfV3NqmDPvpYSx6diKPwAMAgIESE2zKH5ludBmWx5AWAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwvAFGFwAAQH9p7/BpV8NhNbe2KcOerEk5aUpMsBldFvoBgQcAEBcq69wq2Vgvt6ctcM/lSFZxUa4K81wGVob+YMohrbKyMo0ZM0apqalKTU1Vfn6+Nm/eHHjf6/VqyZIlGjp0qFJSUnTNNdfo008/NbBiAICZVda5tai8JijsSFKTp02LymtUWec2qDL0F1MGnuHDh+uxxx7Tnj17tGfPHs2YMUPXXnut9u3bJ0launSpKioqtG7dOm3fvl1Hjx7V1Vdfrfb2doMrBwCYTXuHTyUb6+Xr4j3/vZKN9Wrv6KoFrMLm8/li4jeclpamJ554Qt/61rc0bNgwvfrqq7rxxhslSQcPHlR2drY2bdqkK6+8slfPa2lpkcPhkMfjUWpqajRLBwAY6N0Dh3TT6p09tntt4aXKH5neDxUhHKF+f5uyh+dM7e3tWrdunY4dO6b8/HxVV1fr5MmTmjVrVqBNVlaW8vLytGPHDgMrBQCYUXNrW8+N+tAOscm0k5b37t2r/Px8tbW16Utf+pIqKiqUm5ur2tpaDRo0SOeff35Q+8zMTDU1NXX7PK/XK6/XG7huaWmJWu0AAPPIsCdHtB1ik2l7eL785S+rtrZWO3fu1KJFi7RgwQLV19d3297n88lm635pYWlpqRwOR+CVnZ0djbIBACYzKSdNLkeyuvuGsOn0aq1JOWn9WRb6mWkDz6BBg3TRRRdpwoQJKi0t1dixY/X000/L6XTqxIkTOnLkSFD75uZmZWZmdvu8FStWyOPxBF6NjY3R/iMAAEwgMcGm4qJcSeoUevzXxUW57MdjcaYNPGfz+Xzyer0aP368Bg4cqC1btgTec7vdqqur0+TJk7v9+aSkpMAyd/8LABAfCvNcKps/Tk5H8LCV05Gssvnj2IcnDphyDs8DDzyg2bNnKzs7W62trVq3bp22bt2qyspKORwO3X777brvvvuUnp6utLQ0LV++XKNHj9bMmTONLh0AYFKFeS4V5DrZaTlOmTLw/O1vf9Mtt9wit9sth8OhMWPGqLKyUgUFBZKkVatWacCAAbrhhht0/PhxXXHFFXrppZeUmJhocOUAADNLTLCx9DxOxcw+PJHGPjwAAMQey+7DAwAAEC4CDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsDwCDwAAsLwBRhcAALCW9g6fdjUcVnNrmzLsyZqUk6bEBJvRZSHOEXgAABFTWedWycZ6uT1tgXsuR7KKi3JVmOcysDLEO4a0AAARUVnn1qLymqCwI0lNnjYtKq9RZZ27y59r7/Dp3QOHtL72M7174JDaO3z9US7iDD08AICwtXf4VLKxXl1FFZ8km6SSjfUqyHUGDW/RI4T+Qg8PACBsuxoOd+rZOZNPktvTpl0NhwP3Qu0RAkJB4AEAhK25tfuw01W7nnqEpNM9QgxvIVIIPACAsGXYk/vULpQeISAcBB4AQNgm5aTJ5UhWd4vPbTo9N2dSTpqkvvcIAeEKO/DMnDlTmzdv7nS/vb093EcDAGJEYoJNxUW5ktQp9Pivi4tyAxOW+9ojBIQr7MCzZ88eXXDBBZKkhoaGwP3nn39et9xyS7iPBwDEiMI8l8rmj5PTERxSnI5klc0fF7Tqqq89QkC4wl6WfuLECdntdknS2LFjVVtbqwsvvFCTJ0/WQw89FO7jAQAxpDDPpYJcZ487Lft7hBaV18gmBU1e7qpHCAhX2IHnoosu0nvvvSe73a5jx47piy++kCTZ7XYdPsxkMwCIN4kJNuWPTO+xnb9H6Ox9eJzsw4MoCDvw3HXXXbrjjjs0YsQIjR07Vs8995x+8Ytf6A9/+IMyMzMjUSMAwKJ62yMEhCvswHPnnXdq2LBh2r9/vxYuXKh58+bpwgsvlNvt1t133x2JGgEAFtbbHiEgHDafzxfRXZ1OnTqliooKnThxQvPmzVNiYmIkHx8xLS0tcjgc8ng8Sk1NNbocAADQC6F+f0f8LK0BAwbo+uuvj/RjAQAAQsbGgwAAwPI4LR0A4kB7h4+JwYhrEQ08f/zjHzVhwgQlJSVF8rEAgDBU1rk7Lf12sfQbcSaiQ1qzZ8/WZ599FslHAgDCUFnn1qLymk4HdTZ52rSovEaVdW6DKgP6V0QDT4QXfAEAwtDe4VPJxnp19Tez/17Jxnq1d/B3N6zPlJOWS0tLNXHiRNntdmVkZGju3Ln66KOPgtpMnz5dNpst6DVv3jyDKgYA89nVcLhTz86ZfJLcnjbtamBXfFhfWHN4XnnllaDrU6dO6c0331RGRkbg3ne+850+P7eqqkqLFy/WxIkTderUKf34xz/WrFmzVF9fr5SUlEC7hQsX6uGHHw5cDx48OIQ/BQBYU3Nr92EnlHZALAsr8Lz44otB1ydPntSvfvWrQPCw2WwhBZ7KyspOn5ORkaHq6mpNmzYtcH/IkCFyOp0hVA4A1pdhT+65UR/aAbEsrMDzzjvvBF3b7XatXbtWF154YVhFnc3j8UiS0tLSgu6vWbNG5eXlyszM1OzZs1VcXBw4uf1sXq9XXq83cN3S0hLRGgHAbCblpMnlSFaTp63LeTw2nT6oc1JOWhfvAtZiyjk8Z/L5fFq2bJmmTp2qvLy8wP2bb75Zr732mrZu3aoHH3xQb7zxhq677rpun1NaWiqHwxF4ZWdn90f5AGCYxASbiotyJZ0ON2fyXxcX5bIfD+JCRM/Sstvt+uCDDyLaw7N48WL95je/0fbt2zV8+PBu21VXV2vChAmqrq7WuHHjOr3fVQ9PdnY2Z2kBsDz24YGVmOIsrQceeKDTsFM4lixZog0bNmjbtm3nDDuSNG7cOA0cOFD79+/vMvAkJSWxISKAuFSY51JBrjPsnZbZrRmxLKKBZ8WKFRF5js/n05IlS1RRUaGtW7cqJyenx5/Zt2+fTp48KZeLf60AwNkSE2zKH5ke8s/TS4RYF9EhrUi56667tHbtWq1fv15f/vKXA/cdDocGDx6sAwcOaM2aNbrqqqs0dOhQ1dfX67777tPgwYO1e/duJSYm9vgZoXaJAUC88e/WfPaXhb9vp2z+OEIP+k2o39+mDDw2W9ddpC+++KJuvfVWNTY2av78+aqrq9PRo0eVnZ2tOXPmqLi4uNdDagQeAOhZe4dPUx9/u9sNDP0rvbbfP4PhLfQLU8zhiZSeMlh2draqqqr6qRoAiF992a05nCEzINoiEnhOnjyppqYm/eMf/9CwYcMiOnEZAGAcdmuGVYS8D8/Ro0f17LPPavr06XI4HLrggguUm5urYcOGacSIEVq4cKF2794dyVoBAP2M3ZphFSEFnlWrVumCCy7Q6tWrNWPGDL355puqra3VRx99pHfffVfFxcU6deqUCgoKVFhYqP3790e6bgBAP/Dv1tzd7BybTq/WYrdmmF1Ik5avv/56/fSnP9Xo0aPP2c7r9er555/XoEGDdMcdd4RcZDQwaRkAese/SktS0EotVmnBCJZapdUfCDwA0HvswwOzsNQqLQCAuURqt2bAKCEFng0bNvT5ZwoKCjR48OBQPg4AYALh7tYMGCmkwDN37tw+tbfZbNq/f39EDxUFAADorZCXpTc1Namjo6NXryFDhkSyZgAAgD4JKfAsWLCgT8NT8+fPZ2IwAAAwDKu0WKUFAEDMCPX7u889PMePH9dnn33W6f6+ffv6+igAAIB+0afA86tf/Ur/9m//pquuukpjxozRe++9F3jvlltuiXhxAAAAkdCnwPPII4+opqZGH3zwgV544QXddtttWrt2raSeTzgHAAAwSp+WpZ88eVLDhg2TJE2YMEHbtm3Tddddp7/85S+y2dh8CgAAmFOfengyMjL04YcfBq7T09O1ZcsW/elPfwq6DwAAYCZ9CjyvvvqqMjIygu4NGjRIr732mqqqqiJaGAAAQKT0KfAMHz5cTqcz6N5TTz0lt9utKVOmRLQwAACASAl5p2W/ZcuWacqUKfr000+D7p84cUK7d+8O9/EAAABhCzvwSNLs2bM1bdq0oNBz5MgRXXrppZF4PAAAQFhCOjz0TDabTcXFxcrIyNC0adO0bds2DR8+XBJL1QEAgDmEHXj8iouLJSkQegYOHMhSdQAAYAphB54ze3HODD2vv/56uI8GAACIiLADz8qVK5WSkhK49oeeOXPmhPtoAACAiAg78KxYsaLTveLiYiUmJuo//uM/wn08AABA2Gy+OJ1ZHOrx8gAAwDihfn9HZFk6AACAmRF4AACA5RF4AACA5UU18CQkJGjGjBmqrq6O5scAAACcU1QDzwsvvKDLLrtM99xzTzQ/BgAA4JxYpcUqLQAAYoZhq7SeeuopHTx4MNzHAAAARE3YgWfZsmX6+te/HnRSuiSdOHFCu3fvDvfxAAAAYYvIHJ7CwkJNmzYtKPQcOXJEl156aUjPKy0t1cSJE2W325WRkaG5c+fqo48+Cmrj9Xq1ZMkSDR06VCkpKbrmmms6hS4AAAApAoHHZrOpuLhYCxYs6BR6Qp0eVFVVpcWLF2vnzp3asmWLTp06pVmzZunYsWOBNkuXLlVFRYXWrVun7du36+jRo7r66qvV3t4e7h8JAABYTNiTlhMTE+V2u5WRkaGSkhK9/PLL2rZtmwYOHKisrKyIBJC///3vysjIUFVVlaZNmyaPx6Nhw4bp1Vdf1Y033ihJOnjwoLKzs7Vp0yZdeeWVPT6TScsAAMQewyYtn5mXzuzp+eSTT8J9dIDH45EkpaWlSZKqq6t18uRJzZo1K9AmKytLeXl52rFjR8Q+FwD6qr3Dp3cPHNL62s/07oFDau+Iy4WwgOmEfVr6ypUrlZKSErguLi6WJM2ZMyfcR0s6HaiWLVumqVOnKi8vT5LU1NSkQYMG6fzzzw9qm5mZqaampi6f4/V65fV6A9ctLS0RqQ8A/Crr3CrZWC+3py1wz+VIVnFRrgrzXAZWBiDsHp4VK1YEBR7pdOi55557ZLfbw3287r77bn344Yd67bXXemzr8/lks9m6fK+0tFQOhyPwys7ODrs2APCrrHNrUXlNUNiRpCZPmxaV16iyzm1QZQCkKO60/JOf/ERffPFFWM9YsmSJNmzYoHfeeUfDhw8P3Hc6nTpx4oSOHDkS1L65uVmZmZldPmvFihXyeDyBV2NjY1i1AYBfe4dPJRvr1dXglf9eycZ6hrcAA5ny8FCfz6e7775bb775pt5++23l5OQEvT9+/HgNHDhQW7ZsCdxzu92qq6vT5MmTu3xmUlKSUlNTg14AEAm7Gg536tk5k0+S29OmXQ2H+68oAEHCnsMTDYsXL9batWu1fv162e32wLwch8OhwYMHy+Fw6Pbbb9d9992n9PR0paWlafny5Ro9erRmzpxpcPUA4k1za/dhJ5R2ACLPlIGnrKxMkjR9+vSg+y+++KJuvfVWSdKqVas0YMAA3XDDDTp+/LiuuOIKvfTSS0pMTOznagGYSXuHT7saDqu5tU0Z9mRNyklTYkLXc/siJcOeHNF2ACIvqoeHJiQkaPr06XriiSc0fvz4aH1MSNiHB7Aeo1ZJtXf4NPXxt9XkaetyHo9NktORrO33z4h6+AKszrB9eM7lhRde0GWXXaZ77rknmh8DAIaukkpMsKm4KFfS6XBzJv91cVEuYQcwUFR7eMyMHh7AOvw9LN1NHO6vHhb24QGiL9Tvb1PO4QGAvujLKqn8kelRq6Mwz6WCXGe/zyEC0LOwAs+nn36qsrIy7dixQ01NTbLZbMrMzNTkyZN15513srkfgH5hplVSiQm2qIYqAKEJOfBs375ds2fPVnZ2tmbNmqVZs2bJ5/OpublZb731lv7rv/5Lmzdv1pQpUyJZLwB0Eq1VUkas+AIQHSEHnu9///u64447tGrVqm7fX7p0qXbv3h1ycQDQG5Ny0uRyJPe4SmpSTlqvn8l8HMBaQl6lVVdXpzvvvLPb97/3ve+prq4u1McDQK9FepUU52IB1hNy4HG5XNqxY0e377/77rtyufhXEID+UZjnUtn8cXI6goetnI5klc0f1+teGc7FAqwp5CGt5cuX684771R1dbUKCgqUmZkpm82mpqYmbdmyRb/85S/11FNPRbJWADinSKySMsuKLwCRFXLgueuuu5Senq5Vq1bp2WefVXt7uyQpMTFR48eP1yuvvKIbbrghYoUCQG+Eu0rKTCu+AEROWMvSb7zxRt144406efKkPv/8c0nS0KFDNXDgwIgUBwD9jXOxAGuKyMaDAwcOZL4OAEuIxoovAMaL2llaM2fO1IUXXhitxwNAVHAuFmBNUQs83/jGN7RgwYJoPR4AoiZSK74AmAeHh3J4KIBusNMyYD4cHgoAEca5WIB1RG1Iq7GxUbfddlu0Hg8AANBrUQs8hw8f1ssvvxytxwNAn7R3+PTugUNaX/uZ3j1wiJ2SgTgT8pDWhg0bzvn+xx9/HOqjASCiOAgUQMiTlhMSEmSz2XSuH7fZbIEdmM2GSctAfPAfBHr231T+qcesugJiS6jf32EdHvrGG2+oo6Ojy1dNTU2ojwaAiOAgUAB+IQee8ePHnzPU9NT7AwDR1peDQAFYW8hzeH7wgx/o2LFj3b5/0UUX6Z133gn18QAQNg4CBeAXcuD5+te/fs73U1JSdNlll4X6eAAIGweBAvALaUjrww8/VEdHR6/b79u3T6dOnQrlowAgZP6DQLvbG9mm06u1OAgUsL6QAs/FF1+sQ4cO9bp9fn6+Pvnkk1A+CkAMMOseNxwECsAvpCEtn8+nBx98UEOGDOlV+xMnToTyMQBigNn3uPEfBHp2jU4T1Qgg+kLah2f69Omy2fr2L6K1a9fK5TLPXyzswwOEz+g9bvpyuCcHgQLWEOr3N6elE3iAkLR3+DT18be7XfZt0+lelO33z4hKsDB7zxKA6Oj3jQcBxDcj97jx9yyd/flNnjYtKq9RZZ074p8JILaFFHj6OgH5s88+C+VjAJiYUXvcRHr3ZLNOuAYQWSEFnokTJ2rhwoXatWtXt208Ho9Wr16tvLw8vfnmmyEXCMCcjNrjJpI9S5V1bk19/G3dtHqn7l1Xq5tW79TUx9+mhwiwoJBWaV177bWy2+0qLCzUwIHCXX1DAAATX0lEQVQDNWHCBGVlZSk5OVlHjhxRfX299u3bpwkTJuiJJ57Q7NmzI103AIP597hp8rR12dvin8Pj3+MmUpOGI9Wz1N2Ea/+wGIeKAtYSUuB56aWX1NjYqEceeUSZmZlyuVz6/PPPdfz4cQ0dOlQ333yzrrzySuXl5UW6XgAm4d/jZlF5jWxSUHA4e4+bSE4wjkTPUk/DYjadHhYryHWykguwiJACz7/8y7/o/fffV2FhoY4ePapHH31UGRkZka4NgMn1Zo+bSPek9LVnqSt9GRbLH5ne69oAmFdIc3iWL1+ua665RpMnT5bNZtOaNWu0e/duHT9+PCJFbdu2TUVFRcrKypLNZtNbb70V9P6tt94qm80W9Lr00ksj8tkA+qYwz6Xt98/Qawsv1dPzvqbXFl6q7ffPUGGeK+ITjKXI7J7MoaJA/Akp8CxevFjvv/++rr76avl8Pv33f/+38vPzlZqaqq9+9auaN2+eHnvsMW3evDmkoo4dO6axY8fqmWee6bZNYWGh3G534LVp06aQPgtA+BITbMofma5rv/Yvyh+ZHggb0Vq67u9ZcjqCh62cjuRe9RhxqCgQf0I+LX3UqFEaNWqUXnjhBe3cuVMpKSn68MMPVVtbq9raWq1fv14rV65Ua2trn589e/bsHic6JyUlyel0hlo+gH4QzZ6UwjyXCnKdIU2EjsSwGIDYEnLg8fvLX/4S+O9LLrlEl1xySeA6mps4b926VRkZGTrvvPN02WWXaeXKleecR+T1euX1egPXLS0tUasNwGnR7knx9yyF8nO9nXANwBqiutNyX8/b6q3Zs2drzZo1evvtt/Wzn/1Mu3fv1owZM4ICzdlKS0vlcDgCr+zs7KjUBuCf/D0p3f1NYNPp1VpG9KSEOywGILaY/iwtm82miooKzZ07t9s2brdbI0aM0Lp163Tdddd12aarHp7s7GzO0gKizL9KS+q6J8XocMGhokBsCfUsrbCHtMzA5XJpxIgR2r9/f7dtkpKSlJSU1I9VAZB6t3TdSKEOiwGILZYIPIcOHVJjY6NcLrqgATMKZ4IxAESCKQPP0aNHgyZDNzQ0qLa2VmlpaUpLS9NDDz2kb37zm3K5XPrf//1fPfDAAxo6dKi+8Y1vGFg1gHOhJwWAkUwZePbs2aPLL788cL1s2TJJ0oIFC1RWVqa9e/fqlVde0RdffCGXy6XLL79cr7/+uux2u1ElAwAAEzP9pOVoCXXSEwAAME6o399RXZYOAABgBgQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeabceBCAtXBAJwCjEXgARFVlnbvTwaEukxwcCiB+MKQFIGoq69xaVF4TFHYkqcnTpkXlNaqscxtUGYB4Q+ABEBXtHT6VbKxXV2fX+O+VbKxXe0dcnm4DoJ8ReABExa6Gw516ds7kk+T2tGlXw+H+KwpA3CLwAIiK5tbuw04o7QAgHAQeAFGRYU+OaDsACAeBB0BUTMpJk8uRrO4Wn9t0erXWpJy0/iwLQJwi8ACIisQEm4qLciWpU+jxXxcX5bIfD4B+QeABEDWFeS6VzR8npyN42MrpSFbZ/HHswwOg37DxIICoKsxzqSDXyU7LAAxF4AEQdYkJNuWPTDe6DABxjCEtAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQQeAABgeQOMLgCwmvYOn3Y1HFZza5sy7MmalJOmxASb0WUBQFwzZQ/Ptm3bVFRUpKysLNlsNr311ltB7/t8Pj300EPKysrS4MGDNX36dO3bt8+gaoF/qqxza+rjb+um1Tt177pa3bR6p6Y+/rYq69xGlwYAcc2UgefYsWMaO3asnnnmmS7f//d//3c9+eSTeuaZZ7R79245nU4VFBSotbW1nysF/qmyzq1F5TVye9qC7jd52rSovIbQAwAGsvl8Pp/RRZyLzWZTRUWF5s6dK+l0705WVpaWLl2q+++/X5Lk9XqVmZmpxx9/XN/73vd69dyWlhY5HA55PB6lpqZGrX7Eh/YOn6Y+/nansONnk+R0JGv7/TOiOrzFcBoAqwv1+zvm5vA0NDSoqalJs2bNCtxLSkrSZZddph07dnQbeLxer7xeb+C6paUl6rUi9oQaGHY1HO427EiST5Lb06ZdDYeVPzI9ghX/U2WdWyUb64PqcDmSVVyUq8I8V1Q+EwBiRcwFnqamJklSZmZm0P3MzEz99a9/7fbnSktLVVJSEtXaENvCCQzNrd2HnVDa9ZV/OO3s7lr/cFrZ/HGEHgBxzZRzeHrDZgv+V7fP5+t070wrVqyQx+MJvBobG6NdImJIuPNvMuzJvfqc3rbri/YOn0o21ncKO5IC90o21qu9w9Sj1wAQVTEXeJxOp6R/9vT4NTc3d+r1OVNSUpJSU1ODXoAUmcAwKSdNLkeyuovcNp3uLZqUkxZmtZ31ZTgNAOJVzAWenJwcOZ1ObdmyJXDvxIkTqqqq0uTJkw2sDLEqEoEhMcGm4qJcSeoUevzXxUW5UZlAbPRwGgDEAlMGnqNHj6q2tla1tbWSTk9Urq2t1SeffCKbzaalS5fq0UcfVUVFherq6nTrrbdqyJAh+va3v21w5YhFkQoMhXkulc0fJ6cjeNjK6UiO6hwaI4fTACBWmHLS8p49e3T55ZcHrpctWyZJWrBggV566SX98Ic/1PHjx3XXXXfpyJEjuuSSS/Tb3/5WdrvdqJIRwyIZGArzXCrIdfbr0nD/cFqTp63LYTn/kvhoDKcBQKww/T480cI+PPDz76HTU2CI9h464fBPupYU9GfwV8sqLQBWEer3tymHtID+ZOT8m0gxajgNAGIFPTz08OD/WGHjPnZaBmB1oX5/E3gIPDhDT4GBQAEAxoqboyWAaEpMsHV79IMVeoAAIF4xhwfoBU5CB4DYRuABesDRDQAQ+wg8QA84ugEAYh9zeIAeRGInZiY7A4CxCDxAD8LdiZnJzgBgPIa0gB6EcxI6k50BwBwIPEAPQt2JmcnOAGAeBB6gF0I5uoHJzgBgHszhAXqpryehR2KyMwAgMgg8QB+cayfms4U72RkAEDkMaQFREs5kZwBAZBF4gCgJdbIzACDyCDxAFIUy2RkAEHnM4QGirK+TnQEAkUfgAXoh3KMh+jLZGQAQeQQeoAccDQEAsY85PMA5cDQEAFgDgQfoBkdDAIB1EHiAbnA0BABYB4EH6AZHQwCAdRB4gG5wNAQAWAeBB+gGR0MAgHUQeEykvcOndw8c0vraz/TugUNMhjUYR0MAgHWwD49JsNeLOfmPhjj7d+PkdwMAMcXm8/nishuhpaVFDodDHo9Hqamphtbi3+vl7F+Ev9+AM5eMF+5OywCAyAj1+5seHoP1tNeLTaf3einIdfIFayCOhgCA2MYcHoOx1wsAANFH4DEYe70AABB9BB6DsdcLAADRR+AxGHu9AAAQfQQeg7HXCwAA0ReTgeehhx6SzWYLejmdTqPLCpl/rxenI3jYyulIZkk6AAARELPL0keNGqXf/e53gevExEQDqwlfYZ5LBblO9noBACAKYjbwDBgwIKZ7dbrCXi8AAERHTA5pSdL+/fuVlZWlnJwczZs3Tx9//PE523u9XrW0tAS9AABAfIjJwHPJJZfolVde0f/8z/9o9erVampq0uTJk3Xo0KFuf6a0tFQOhyPwys7O7seKAQCAkSxxltaxY8c0cuRI/fCHP9SyZcu6bOP1euX1egPXLS0tys7ONsVZWgAAoHfi+iytlJQUjR49Wvv37++2TVJSkpKSkvqxKgAAYBYxOaR1Nq/Xqz/96U9yuVi+DQAAOovJwLN8+XJVVVWpoaFB7733nr71rW+ppaVFCxYsMLo0AABgQjE5pPXpp5/qpptu0ueff65hw4bp0ksv1c6dOzVixAijSwMAACYUk4Fn3bp1RpcAAABiSEwOaQEAAPQFgQcAAFgegQcAAFgegQcAAFgegQcAAFgegQcAAFgegQcAAFgegQcAAFgegQcAAFheTO60bGbtHT7tajis5tY2ZdiTNSknTYkJNqPLAgAgrhF4Iqiyzq2SjfVye9oC91yOZBUX5aowj5PcAQAwCkNaEVJZ59ai8pqgsCNJTZ42LSqvUWWd26DKAAAAgScC2jt8KtlYL18X7/nvlWysV3tHVy0AAEC0EXgiYFfD4U49O2fySXJ72rSr4XD/FQUAAAIIPBHQ3Np92AmlHQAAiCwCTwRk2JMj2g4AAEQWgScCJuWkyeVIVneLz206vVprUk5af5YFAAD+D4EnAhITbCouypWkTqHHf11clMt+PAAAGITAEyGFeS6VzR8npyN42MrpSFbZ/HHswwMAgIHYeDCCCvNcKsh1stMyAAAmQ+CJsMQEm/JHphtdBgAAOANDWgAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPLidqdln88nSWppaTG4EgAA0Fv+723/93hvxW3gaW1tlSRlZ2cbXAkAAOir1tZWORyOXre3+foakSyio6NDBw8elN1ul83G4Z79raWlRdnZ2WpsbFRqaqrR5eAM/G7Mi9+NOfF76V8+n0+tra3KyspSQkLvZ+bEbQ9PQkKChg8fbnQZcS81NZW/IEyK34158bsxJ34v/acvPTt+TFoGAACWR+ABAACWl/jQQw89ZHQRiE+JiYmaPn26BgyI25FV0+J3Y178bsyJ34v5xe2kZQAAED8Y0gIAAJZH4AEAAJZH4AEAAJZH4AEAAJZH4EG/Ki0t1cSJE2W325WRkaG5c+fqo48+MrosnKW0tFQ2m01Lly41uhRI+uyzzzR//nylp6dryJAh+trXvqbq6mqjy4p7p06d0k9+8hPl5ORo8ODBuvDCC/Xwww+ro6PD6NLQBdbPoV9VVVVp8eLFmjhxok6dOqUf//jHmjVrlurr65WSkmJ0eZC0e/duPffccxozZozRpUDSkSNHNGXKFF1++eXavHmzMjIydODAAZ133nlGlxb3Hn/8cf3iF7/Qyy+/rFGjRmnPnj367ne/K4fDoXvvvdfo8nAWlqXDUH//+9+VkZGhqqoqTZs2zehy4t7Ro0c1btw4/fznP9cjjzyir33ta3rqqaeMLiuu/ehHP9If//hH/eEPfzC6FJzl6quvVmZmpp5//vnAvW9+85saMmSIXn31VQMrQ1cY0oKhPB6PJCktLc3gSiBJixcv1pw5czRz5kyjS8H/2bBhgyZMmKDrr79eGRkZuvjii7V69Wqjy4KkqVOn6ve//73+/Oc/S5I++OADbd++XVdddZXBlaErDGnBMD6fT8uWLdPUqVOVl5dndDlxb926daqpqdHu3buNLgVn+Pjjj1VWVqZly5bpgQce0K5du3TPPfcoKSlJ3/nOd4wuL67df//98ng8+spXvqLExES1t7dr5cqVuummm4wuDV0g8MAwd999tz788ENt377d6FLiXmNjo+6991799re/VXJystHl4AwdHR2aMGGCHn30UUnSxRdfrH379qmsrIzAY7DXX39d5eXlWrt2rUaNGqXa2lotXbpUWVlZWrBggdHl4SwEHhhiyZIl2rBhg7Zt26bhw4cbXU7cq66uVnNzs8aPHx+4197erm3btumZZ56R1+tVYmKigRXGL5fLpdzc3KB7X/3qV/XGG28YVBH8fvCDH+hHP/qR5s2bJ0kaPXq0/vrXv6q0tJTAY0IEHvQrn8+nJUuWqKKiQlu3blVOTo7RJUHSFVdcob179wbd++53v6uvfOUruv/++wk7BpoyZUqnrRv+/Oc/a8SIEQZVBL9//OMfSkgIngqbmJjIsnSTIvCgXy1evFhr167V+vXrZbfb1dTUJElyOBwaPHiwwdXFL7vd3mkeVUpKitLT05lfZbDvf//7mjx5sh599FHdcMMN2rVrl5577jk999xzRpcW94qKirRy5Ur967/+q0aNGqX3339fTz75pG677TajS0MXWJaOfmWz2bq8/+KLL+rWW2/t32JwTtOnT2dZukn8+te/1ooVK7R//37l5ORo2bJlWrhwodFlxb3W1lY9+OCDqqioUHNzs7KysnTTTTfppz/9qQYNGmR0eTgLgQcAAFge+/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAsIxHH31UNput0+vJJ580ujQABuMsLQCW0draqmPHjgWuH374YW3atEnbt2/X8OHDDawMgNEGGF0AAESK3W6X3W6XJJWUlGjTpk2qqqoi7ABgSAuA9ZSUlOjFF19UVVWVRowYYXQ5AEyAwAPAUgg7ALpC4AFgGYQdAN1hDg8AS3jkkUf0zDPP6Ne//rWSkpLU1NQkSTr//POVlJRkcHUAjMYqLQAxz+fz6bzzzlNLS0un93bu3KlLLrnEgKoAmAmBBwAAWB5zeAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOUReAAAgOX9f7JW0wlzmiuGAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject