{
 "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",
    "$$\n",
    "p(x_2|x_3) = \\frac{p(x_2,x_3)}{p(x_3)} = \\frac{\\int \\cdots \\int p(x_1,x_2,x_3,x_4,x_5) \\, \\mathrm{d}x_1  \\mathrm{d}x_4 \\mathrm{d}x_5}{\\int \\cdots \\int 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}$$"
   ]
  },
  {
   "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": [
    "- <a id='factorization-helps'>Factorization helps</a>. 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",
    "$$\n",
    "p(x_2|x_3) = \\frac{\\int \\cdots \\int 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}{\\int \\cdots \\int 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",
    "$$\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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-example-1.png\" width=\"400px\"></p>"
   ]
  },
  {
   "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",
    " <p style=\"text-align:center;\"><img src=\"./figures/ffg-wEquality-node.png\" width=\"450px\"></p>\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- 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 }{ \\int\\cdots\\int f(x_1,x_2,x_3,x_4) \\,\\mathrm{d}x_1 \\mathrm{d}x_3 \\mathrm{d}x_4} \\\\\n",
    "  &= \\frac{\\int\\cdots\\int 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 }{ \\int\\cdots\\int 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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-example-prob-model.png\" width=\"400px\"></p>"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-message-passing-1.png\" width=\"600px\"><p style=\"text-align:center;\">"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-message-passing.png\" width=\"700px\"></p>"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-sum-product.png\" width=\"500px\"></p>\n",
    "\n",
    "- This is called the **Sum-Product Message** (SPM) update rule. (Look at the formula to understand why it's called the 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** (or alternatively, the **belief propagation** algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Just as a final note, inference by sum-product message passing is much like replacing the sum-of-products\n",
    "$$\n",
    "ac + ad + bc + b$$\n",
    "by the following product-of-sums:\n",
    "$$\n",
    "(a + b)(c + d) \\,.$$\n",
    "Which of these two computations is cheaper to execute?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <a id=\"sp-for-equality-node\">Sum-Product Messages for the Equality Node</a>\n",
    "\n",
    "- As an example, let´s evaluate the SP messages for the **equality node** $f_=(x,y,z) = \\delta(z-x)\\delta(z-y)$: \n",
    "\n",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-equality-node.png\" width=\"250px\"></p>\n",
    "\n",
    "$$\\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*}$$\n",
    "\n",
    "- By symmetry, this also implies (for the same equality node) that\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*}$$\n",
    "\n",
    "- 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 (see [Roweis notes](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Roweis-1999-gaussian-identities.pdf)), 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*}$$\n",
    "\n",
    "- 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": [
    "### 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 to represent observations, e.g., add a factor $f(y)=\\delta(y−3)$ to terminate the half-edge for variable $Y$  if  $y=3$  is observed.\n",
    "\n",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-observation-y-3.png\" width=\"600px\"></p>\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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-terminal-node-message.png\" width=\"250px\"></p>\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 more update rules for 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 [RxInfer](http://rxinfer.ml), 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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-bayesian-linear-regression.png\" width=\"500px\"></p>\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 [RxInfer](https://biaslab.github.io/rxinfer-website/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuWklEQVR4nO3df3RU9Z3/8ddNwEnCTsYGzcxEUIOyakxRAaEggmsNG2qzuvZUrcW667ceQKymnFaLbjfG1aRgF9ldWlasunaRr57TYpVdyTGntUEEmgimGNK1ijkaMdOsBCcBSQLJ/f7BN9MMSSCZX/fOvc/HOXNO587Nzdu0x3n18+P9MUzTNAUAAOAgGVYXAAAAkGgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DjjrC7AKv39/frkk0/k9XplGIbV5QAAgFEwTVNdXV0qKChQRsbI4zSuDTiffPKJJk+ebHUZAAAgBq2trZo0adKIn7s24Hi9Xkkn/kC5ubkWVwMAAEajs7NTkydPjnyPj8S1AWdgWio3N5eAAwBAmjnd8hIWGQMAAMch4AAAAMch4AAAAMch4AAAAMch4AAAAMch4AAAAMch4AAAAMch4AAAAMch4AAAAMdxbSdjAACQeH39pupbOtTe1a18b5ZmFeYpMyP1h1oTcAAAQELUNLWpckuz2sLdkWtBX5YqyopUWhxMaS1MUQEAgLjVNLVp2cY9UeFGkkLhbi3buEc1TW0prYeAAwAA4tLXb6pyS7PMYT4buFa5pVl9/cPdkRwEHAAAEJf6lo4hIzeDmZLawt2qb+lIWU0EHAAAEJf2rpHDTSz3JQIBBwAAxCXfm5XQ+xKBgAMAAOIyqzBPQV+WRtoMbujEbqpZhXkpq4mAAwAA4pKZYaiirEiShoScgfcVZUUp7YdDwAEAAHErLQ5q/eLpCviip6ECviytXzw95X1waPQHAAASorQ4qJKiAJ2MAQCAs2RmGJpzwUSry2CKCgAAOA8BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA59cAAAcIm+ftMWTfhSgYADAIAL1DS1qXJLs9rC3ZFrQV+WKsqKUn6MQiowRQUAgMPVNLVp2cY9UeFGkkLhbi3buEc1TW0WVZY8BBwAABysr99U5ZZmmcN8NnCtckuz+vqHuyN9EXAAAHCw+paOISM3g5mS2sLdqm/pSF1RKUDAAQDAwdq7Rg43sdyXLgg4AAA4WL43K6H3pQsCDgAADjarME9BX5ZG2gxu6MRuqlmFeaksK+kIOAAAOFhmhqGKsiJJGhJyBt5XlBU5rh8OAQcAAIcrLQ5q/eLpCviip6ECviytXzzdkX1waPQHAIALlBYHVVIUoJMxAABwlswMQ3MumGh1GSlBwAEAIE246SypeBFwAABIA247SypeLDIGAMDm3HiWVLwIOAAA2Jhbz5KKFwEHAAAbc+tZUvEi4AAAYGNuPUsqXgQcAABszK1nScWLgAMAgI259SypeBFwAACwsUSdJdXXb2rn/oN6ufGAdu4/6PhFyfTBAQDA5gbOkjq5D05glH1w3NhDxzBN09kRbgSdnZ3y+XwKh8PKzc21uhwAAE4rlk7GAz10Tv6yH/ipdDtsc7Tf37afoqqurpZhGCovL49cM01TDz/8sAoKCpSdna1rrrlG+/bts65IAABSYOAsqRsuP0dzLpg4qmkpt/bQsXXAaWho0IYNGzRt2rSo66tXr9aaNWu0bt06NTQ0KBAIqKSkRF1dXRZVCgCA/bi5h45tA87hw4f1zW9+U0899ZS+8IUvRK6bpqm1a9fqoYce0k033aTi4mI999xz+vzzz7Vp06YRn9fT06POzs6oFwAATubmHjq2DTjLly/X9ddfr+uuuy7qektLi0KhkBYuXBi55vF4tGDBAu3YsWPE51VXV8vn80VekydPTlrtAADYgZt76Ngy4LzwwgvavXu3qqurh3wWCoUkSX6/P+q63++PfDaclStXKhwOR16tra2JLRoAAJtxcw8d2wWc1tZW3XfffXr++eeVlTVyojSM6P+6TNMccm0wj8ej3NzcqBcAAE6WqB466ch2AWf37t1qb2/XjBkzNG7cOI0bN051dXX613/9V40bNy4ycnPyaE17e/uQUR0AANxuoIdOwBc9aBDwZaXdFvGxsF2jvy9/+ct65513oq79/d//vS6++GI98MADmjJligKBgGpra3XFFVdIknp7e1VXV6dVq1ZZUTIAALZWWhxUSVFgzD100pntAo7X61VxcXHUtQkTJmjixImR6+Xl5aqqqtLUqVM1depUVVVVKScnR7fddpsVJQMAYHsDPXTcwnYBZzTuv/9+HT16VHfffbcOHTqk2bNn67XXXpPX67W6NAAAYAMc1cBRDQAApI3Rfn+n5QgOAAB2F8u5UUgcAg4AAAnmxtO77cZ228QBAEhnA6d3n3wGVCjcrWUb96imqc2iytyFgAMAQIK4+fRuuyHgAACQIG4+vdtuCDgAACSIm0/vthsCDgAACeLm07vthoADAECCuPn0brsh4AAAkCCnOr1bOrEG54fXO/P0brsh4AAAkEAjnd494J/+u5mt4ilAwAEAIMFKi4P64fWXDPsZ/XBSg4ADAECC9fWb+qf//sOwn9EPJzUIOAAAJBj9cKxHwAEAIMHoh2M9Ag4AAAlGPxzrEXAAAEgw+uFYj4ADAECCnaofzsD7ijL64SQTAQcAgCQYqR9OwJel9Yunq7Q4aFFl7jDO6gIAAHCq0uKgSooCqm/pUHtXt/K9J6alGLlJPgIOAABJlJlhaM4FE60uw3UIOAAAJEhfv8lojU0QcAAASICapjZVbmmOavAX9GWpoqyI9TYWYJExAABxqmlq07KNe4Z0L+bcKesQcAAAiENfv6nKLc0a7lQpzp2yDgEHAIA4cO6UPRFwAACIA+dO2RMBBwCAOHDulD0RcAAAiAPnTtkTAQcA4Ep9/aZ27j+olxsPaOf+gzEvAubcKXuiDw4AwHUS3bNm4Nypk58ZoA+OZQzTNF25b62zs1M+n0/hcFi5ublWlwMASJGBnjUnf/kNjK/EcxAmnYyTb7Tf34zgAABc43Q9awyd6FlTUhSIKZhw7pR9sAYHAOAa9KxxDwIOAMA16FnjHgQcAIBr0LPGPQg4AADXOF3PGkk6M3u8+k2Ts6PSHAEHAOAap+pZM+Czo8f0zZ/9TvNW/YZTwNMYAQcA4CoDPWsCvlNPQ4XC3Vq2cQ8hJ00RcAAArlNaHNT2B67V8/9nts7MHj/sPQMTVJVbmpmuSkMEHACAK2VmGMrIMPTZ0WMj3sO28fRFoz8AgOsMdBzeOsrpJ7aNpx8CDgDAVYY7h+p02Daefgg4AADXGOkcqpEYOnFg5qzCvGSWhSRgDQ4AwBVOdQ7VcAa2kVeUFXFgZhpiBAcA4AqnO4fqZAFflirKimI+WRzWIuAAAFxhtAuFvzXnPC0qDmpWYR4jN2mMgAMAcIXRLhReVBzUnAsmJrkaJBsBBwDgWAPbwdu7unXWBI8CuVn6U2f3sOtwWFDsLAQcAIAjDbcd/Myc8TJ1IswMDjksKHYeAg4AwHFG2g4e/vxE12Jfznh99vmfOxizoNh5CDgAAEc51XbwgdGbrHEZev7bs/Xp4R7le7NYUOxABBwAgKOcbju4KSnU2aMMw9ANl5+TusKQUjT6AwA4ymi3g3O+lLMRcAAAjjLa7eCcL+VsBBwAgKPMKsxT0JelkVbUGJKCbAd3PAIOAMBRMjMMVZQVSdKQkMN2cPcg4AAAHKe0OKj1i6cr4Iuehgr4srR+8XS2g7sAu6gAAI5UWhxUSVEg0smY7eDuQsABADhWZobBuVIuxRQVAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHNsFnPXr12vatGnKzc1Vbm6u5syZo61bt0Y+N01TDz/8sAoKCpSdna1rrrlG+/bts7BiAABgN7YLOJMmTdKPfvQjvfXWW3rrrbd07bXX6oYbboiEmNWrV2vNmjVat26dGhoaFAgEVFJSoq6uLosrBwAAdmGYpmlaXcTp5OXl6fHHH9edd96pgoIClZeX64EHHpAk9fT0yO/3a9WqVVqyZMmon9nZ2Smfz6dwOKzc3NxklQ4AABJotN/fcY/gXHfddVFTSAP6+vrifbT6+vr0wgsv6MiRI5ozZ45aWloUCoW0cOHCyD0ej0cLFizQjh07Tvmsnp4edXZ2Rr0AAIAzxR1w3nrrLZ1//vmSpJaWlsj1p59+WrfffntMz3znnXf0F3/xF/J4PFq6dKleeuklFRUVKRQKSZL8fn/U/X6/P/LZSKqrq+Xz+SKvyZMnx1QbAACwv7gDTm9vr7xeryTpsssu0wcffCBJmjt3rn7961/H9MyLLrpIjY2N2rVrl5YtW6Y77rhDzc3Nkc8NI/qgNNM0h1w72cqVKxUOhyOv1tbWmGoDAAD2F/dhmxdeeKF+97vfyev16siRI/rss88kSV6vVx0dHTE984wzztCFF14oSZo5c6YaGhr0L//yL5F1N6FQSMHgn4+6b29vHzKqczKPxyOPxxNTPQAAIL3EPYJz991369vf/rYWLFigyy67TBs2bJAkvfHGG6cNHaNlmqZ6enpUWFioQCCg2trayGe9vb2qq6vT3LlzE/K7AABA+ot7BGfp0qU6++yz9d577+muu+7SrbfeqilTpqitrU333HPPmJ/34IMPatGiRZo8ebK6urr0wgsv6Le//a1qampkGIbKy8tVVVWlqVOnaurUqaqqqlJOTo5uu+22eP9RAACAQ8QdcCTpa1/7WuQ/b926VS+99JJ6e3t16623jvlZf/rTn3T77berra1NPp9P06ZNU01NjUpKSiRJ999/v44ePaq7775bhw4d0uzZs/Xaa69F1gEBAACkRR+cZKAPDgAA6SdlfXAAAADshoADAAAcJ6EB580331RPT08iHwkAADBmCQ04ixYt0oEDBxL5SABAmujrN7Vz/0G93HhAO/cfVF+/K5d4wiYSsotqgEvXKwOA69U0talyS7Pawt2Ra0FflirKilRaHDzFTwLJwRocAEBcapratGzjnqhwI0mhcLeWbdyjmqY2iyqDm8U1gvPzn/886v3x48e1efNm5efnR65961vfiudXAABsrK/fVOWWZg03fm9KMiRVbmlWSVFAmRmnPjMQSKS4As6zzz4b9f7YsWP6xS9+oezsbEknDsUk4ACAc9W3dAwZuRnMlNQW7lZ9S4fmXDAxdYXB9eIKOK+//nrUe6/Xq02bNmnKlClxFQUASA/tXSOHm1juAxKFNTgAgJjle7MSeh+QKAQcAEDMZhXmKejL0kirawyd2E01qzAvlWUBiQ04Dz74oPLy+B8xALhFZoahirIiSRoScgbeV5QVscAYKcdhmxy2CQBxow8OUmW0398JbfQHAHCn0uKgSooCqm/pUHtXt/K9J6alGLmBVQg4AICEyMww2AoO22CRMQAAcJyEjOAcO3ZMoVBIn3/+uc4++2wWGgMAAEvFPIJz+PBhPfnkk7rmmmvk8/l0/vnnq6ioSGeffbbOO+883XXXXWpoaEhkrQAAAKMSU8B54okndP755+upp57Stddeq82bN6uxsVHvvvuudu7cqYqKCh0/flwlJSUqLS3Ve++9l+i6AQAARhTTNvGvf/3r+sd//Ed98YtfPOV9PT09evrpp3XGGWfo29/+dsxFJgPbxAEASD+j/f6mDw4BBwCAtDHa7292UQEAAMeJaRfVK6+8MuafKSkpUXZ2diy/DgAAYExiCjg33njjmO43DEPvvfeepkyZEsuvAwAAGJOYp6hCoZD6+/tH9crJyUlkzQAAAKcUU8C54447xjTdtHjxYhbyAgCAlGEXFbuoAABIG0nbRXX06FEdOHBgyPV9+/aN9VEAAABJMaaA84tf/EJ/+Zd/qa985SuaNm2afve730U+u/322xNeHAAAQCzGFHAeffRR7dmzR7///e/1zDPP6M4779SmTZskSS6d6QIAADY0pm3ix44d09lnny1JmjlzprZt26abbrpJ77//vgzDSEqBAAAAYzWmEZz8/Hzt3bs38n7ixImqra3VH/7wh6jrAIA/6+s3tXP/Qb3ceEA79x9UXz8j3kCyjWkX1ccff6xx48YpEAgM+ezNN9/UVVddldDikoldVABSoaapTZVbmtUW7o5cC/qyVFFWpNLioIWVAekpKbuoJk2aNCTcrF27Vm1tbWkVbgAgFWqa2rRs456ocCNJoXC3lm3co5qmNosqA5wv7sM2V6xYoauuukoff/xx1PXe3l41NDTE+3gASEt9/aYqtzRruCHygWuVW5qZrgKSJCGniS9atEjz58+PCjmHDh3Sl770pUQ8HgDSTn1Lx5CRm8FMSW3hbtW3dKSuKMBFYjpsczDDMFRRUaH8/HzNnz9f27Zt06RJkySxdRyAe7V3jRxuYrkPwNjEHXAGVFRUSFIk5IwfP56t4wBcK9+bldD7AIxN3AFn8CjN4JDz4osvxvtoAEhbswrzFPRlKRTuHnYdjiEp4MvSrMK8VJcGuELcAeexxx7ThAkTIu8HQs71118f76MBIG1lZhiqKCvSso17ZEhRIWdgbLuirEiZGYx0A8mQtNPEH330Uf34xz/WZ599lozHx40+OABSgT44QGKN9vs7aQHH7gg4AFKlr99UfUuH2ru6le89MS3FyA0Qm9F+fydskTEAYHiZGYbmXDDR6jIAV0lIHxwAAAA7IeAAAADHSWrAycjI0LXXXqvdu3cn89cAAABESWrAeeaZZ7RgwQLde++9yfw1AAAAUdhFxS4qAADSxmi/v+MewVm7dq0++eSTeB8DAACQMHEHnBUrVujqq6+OOklcknp7e9XQ0BDv4wEAAMYsIWtwSktLNX/+/KiQc+jQIX3pS19KxOMBAADGJO5Gf4ZhqKKiQvn5+ZGTxCdNmiQp+iBOAACAVElYJ+PBJ4lv27ZN48ePl2HQihwAAKRe3AFn8CjN4JDz4osvxvtoAACAmMQdcB577DFNmDAh8n4g5Fx//fXxPhoAACAmcQeclStXDrlWUVGhzMxM/fjHP4738QAAAGNGoz8a/QEAkDZS1ugPAADAbgg4AADAcQg4AADAcZIacDIyMnTttddq9+7dyfw1AAAAUZIacJ555hktWLBA9957bzJ/DQAAQBR2UbGLCgCAtMEuKgAA4FpxNfr7+OOPtX79eu3YsUOhUEiGYcjv92vu3LlaunSpJk+enKg6AQAARi3mKart27dr0aJFmjx5shYuXCi/3y/TNNXe3q7a2lq1trZq69atuuqqqxJdc0IwRQUAQPoZ7fd3zAHnyiuv1Lx58/TEE08M+/l3v/tdbd++XQ0NDWN6bnV1tTZv3qz/+Z//UXZ2tubOnatVq1bpoosuitxjmqYqKyu1YcMGHTp0SLNnz9ZPfvITXXrppaP+PQQcAADST9LX4DQ1NWnp0qUjfr5kyRI1NTWN+bl1dXVavny5du3apdraWh0/flwLFy7UkSNHIvesXr1aa9as0bp169TQ0KBAIKCSkhJ1dXXF9M8CAACcJeY1OMFgUDt27IgaWRls586dCgaDY35uTU1N1Ptnn31W+fn52r17t+bPny/TNLV27Vo99NBDuummmyRJzz33nPx+vzZt2qQlS5aM/R8GAAA4SswB53vf+56WLl2q3bt3q6SkRH6/X4ZhKBQKqba2Vj/72c+0du3auAsMh8OSpLy8PElSS0uLQqGQFi5cGLnH4/FowYIF2rFjx4gBp6enRz09PZH3nZ2dcdcGAADsKeaAc/fdd2vixIl64okn9OSTT6qvr0+SlJmZqRkzZujnP/+5br755riKM01TK1as0Lx581RcXCxJCoVCkiS/3x91r9/v14cffjjis6qrq1VZWRlXPQAAID3EtU38lltu0S233KJjx47p008/lSSdddZZGj9+fEKKu+eee7R3715t3759yGeGYUS9N01zyLXBVq5cqRUrVkTed3Z2so0dAACHiivgDBg/fnxM621O5Tvf+Y5eeeUVbdu2TZMmTYpcDwQCkk6M5Az+ne3t7UNGdQbzeDzyeDwJrREAANhT0joZX3fddZoyZcqYf840Td1zzz3avHmzfvOb36iwsDDq88LCQgUCAdXW1kau9fb2qq6uTnPnzo27bgAAkP4SMoIznL/927+NTFuNxfLly7Vp0ya9/PLL8nq9kTU3Pp9P2dnZMgxD5eXlqqqq0tSpUzV16lRVVVUpJydHt912W6L/MQAAQBqy3WGbI62jefbZZ/V3f/d3kv7c6O/JJ5+MavQ3sBB5NGj0BwBA+kl6J+N0R8ABACD9WH6aeGtrq+68885kPR4AAGBESQs4HR0deu6555L1eAAu1ddvauf+g3q58YB27j+ovn5XDkIDOI2YFxm/8sorp/z8gw8+iPXRADCsmqY2VW5pVlu4O3Itb8J4PXpDsb4yrcDCygDYTcxrcDIyMmQYhk7144ZhRDoc2w1rcID0UtPUpmUb92ikf+MsmV+olV8pSmlNAFIv6WtwgsGgfvnLX6q/v3/Y1549e2J9NABE6es3VbmlecRwI0lPbmvRq3vbUlYTAHuLOeDMmDHjlCHmdKM7ADBa9S0dUdNSI/nhy02syQEgKY41ON///vd15MiRET+/8MIL9frrr8f6eACIaO86fbiRpINHelXf0qE5F0xMckUA7C7mgHP11Vef8vMJEyZowYIFsT4eACLyvVmjvne0YWhAX7+p+pYOtXd1K9+bpVmFecrMGPngXgDpIaaAs3fvXhUXFysjY3QzXPv27dNFF12kceOSdjIEAIcZHDzOmuDRF3LG69Dnx077c2MJQ8Ptygr6slRRVqTS4sQeIAwgtWLaRZWZmalQKKSzzz57VPfn5uaqsbExpsM3k4VdVIB9DRc8fNnjFD56/JQ/F/RlafsD145qBGakXVkDP7l+8XRCDmBDo/3+jmlIxTRN/fCHP1ROTs6o7u/t7Y3l1wBwoZGCR+dpwo0hqaKsaFTh5lS7ssz//6zKLc0qKQowXQWkqZgCzvz58/Xuu++O+v45c+YoOzs7ll8FwEVGEzx82eOUkZGhjiN//j9OY51WOt2uLFNSW7ibBctAGosp4Pz2t79NcBkAMLrg8dnR43r+/8xWRoYR88Lg0S5EHuuCZQD2wapfALYx2kDx6ZEe3XD5OWN69uBFy5929YzqZ8ayYBmAvcQUcD766COde+65o77/wIEDOuecsf3LCID7jDZQjDV4DLdoOcOQRuoJaEgK+E6MDAFITzF1Mr7yyit11113qb6+fsR7wuGwnnrqKRUXF2vz5s0xFwjAPWYV5inoy9JIk02GTqy3GUvwGFi0fPLU16nCjTT6BcsA7CmmEZwbbrhBXq9XpaWlGj9+vGbOnKmCggJlZWXp0KFDam5u1r59+zRz5kw9/vjjWrRoUaLrBuBAmRmGKsqKtGzjHhlS1GLjWILHaM6wOnkkJ0AfHMARYuqDc8YZZ6i1tVW5ubny+/26+eabdfDgQR09elRnnXWWrrjiCv31X/+1iouLk1FzQtAHB7CvRDXg27n/oL7x1K7T3vfD6y/RWV4PnYyBNJDUPjjnnHOO3n77bZWWlurw4cOqqqpSfn5+zMUCwGClxUGVFAXiPkJhtIuWz/J6xrxoGYC9xRRwvve97+lv/uZvNHPmTBmGoeeff17z5s1TcXEx/W4AJERmhhF3D5pkLVoGYH8xLTJevny53n77bX31q1+VaZr6yU9+ojlz5ig3N1eXXHKJbr31Vv3oRz/S1q1bE10vAIxaMhYtA0gPMa3BGezCCy/Url27NGHCBO3du1eNjY2RV1NTk7q6uhJVa0KxBgdwh4FdVNLwi5Y5cwpIL6P9/o474JyKaZoyDHsu1iPgAO7BqeGAcyR1kfFo2TXcAHCXRC1aBpA+OKoBgCskYtEygPQR0yJjAAAAOyPgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAxyHgAAAAx+E0ccBCff2m6ls61N7VrXxvlmYV5ikzw7C6LABIewQcwCI1TW2q3NKstnB35FrQl6WKsiKVFgctrAwA0h9TVIAFapratGzjnqhwI0mhcLeWbdyjmqY2iyoDAGcg4AAp1tdvqnJLs8xhPhu4VrmlWX39w90BABgNAg6QYvUtHUNGbgYzJbWFu1Xf0pG6ogDAYQg4QIq1d40cbmK5DwAwFAEHSLF8b1ZC7wMADEXAAVJsVmGegr4sjbQZ3NCJ3VSzCvNSWRYAOAoBB0ixzAxDFWVFkjQk5Ay8rygroh8OAMSBgANYoLQ4qPWLpyvgi56GCviytH7xdPrgAECcaPQHWKS0OKiSogCdjAEgCQg4gIUyMwzNuWCi1WUAgOMwRQUAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHoxqANNDXb3JmFQCMAQEHsLmapjZVbmlWW7g7ci3oy1JFWRGnjgPACJiiAmyspqlNyzbuiQo3khQKd2vZxj2qaWqzqDIAsDcCDmBTff2mKrc0yxzms4FrlVua1dc/3B0A4G4EHMCm6ls6hozcDGZKagt3q76lI3VFAUCasGXA2bZtm8rKylRQUCDDMPSrX/0q6nPTNPXwww+roKBA2dnZuuaaa7Rv3z5rigWSpL1r5HATy30A4Ca2DDhHjhzRZZddpnXr1g37+erVq7VmzRqtW7dODQ0NCgQCKikpUVdXV4orBZIn35uV0PsAwE1suYtq0aJFWrRo0bCfmaaptWvX6qGHHtJNN90kSXruuefk9/u1adMmLVmyJJWlAkkzqzBPQV+WQuHuYdfhSNKZOeM1qzAvpXUBQDqw5QjOqbS0tCgUCmnhwoWRax6PRwsWLNCOHTtG/Lmenh51dnZGvQA7y8wwVFFWNGK4kaTPPj+m2uZQymoCgHSRdgEnFDrxL3O/3x913e/3Rz4bTnV1tXw+X+Q1efLkpNYJjKSv39TO/Qf1cuMB7dx/8JS7oEqKAjozZ/yInxtiJxUADMeWU1SjYRjRXVxN0xxybbCVK1dqxYoVkfednZ2EHKTcWJv21bd06LPPj434vME7qeZcMDEZJQNAWkq7EZxAICBJQ0Zr2tvbh4zqDObxeJSbmxv1AlIplqZ97KQCgNikXcApLCxUIBBQbW1t5Fpvb6/q6uo0d+5cCysDRhZr0z52UgFAbGw5RXX48GG9//77kfctLS1qbGxUXl6ezj33XJWXl6uqqkpTp07V1KlTVVVVpZycHN12220WVg2MbCxN+wZPNZ1uJ5UhKeDLYicVAJzElgHnrbfe0l/91V9F3g+snbnjjjv0H//xH7r//vt19OhR3X333Tp06JBmz56t1157TV6v16qSgVOKdappYCfVso17ZEhRIWdgxVlFWREniwPASQzTNF25/aKzs1M+n0/hcJj1OEi6nfsP6htP7Trtff/3ri8Nu1iYE8UB4ITRfn/bcgQHcJp4p5pKi4MqKQqovqVD7V3dyveeuJeRGwAYHgEHSIFETDVlZhhsBQeAUUq7XVRAuiotDmr94ukK+KJ3PAV8WVq/eDpTTQCQQIzgACnEVBMApAYBBxhBX7+ZlCDCVBMAJB8BBxgGu5YAIL2xBgc4SSxHKgAA7IWAAwwS65EKAAB7IeAAg4zlSAUAgH0RcIBBOL0bAJyBgAMMwundAOAMBBxgkIEjFUbaDG7oxG4qTu8GAHsj4ACDDBypIGlIyOH0bgBIHwQcpERfv6md+w/q5cYD2rn/oK13IXGkAgCkPxr9IenSsWkeRyoAQHozTNO07/+VTqLOzk75fD6Fw2Hl5uZaXY5jDTTNO/l/ZAMxgRERAMBYjPb7mykqJI2dm+al05QZAGDsmKJC0oylaV4qD59MxykzAMDYMIKDpLFj0zzOmQIAdyDguEwqp2bs1jTPzlNmAIDEYorKRVI9NTPQNC8U7h42VBg6sfU6VU3z7DplBgBIPEZwXMKKqRm7Nc2z45QZACA5CDg2l4gpJSunZuzUNM9uU2YAgORhisrGEjWlZPXUjF2a5tltygwAkDyM4NhUIqeU7DA1k5lhaM4FE3XD5edozgUTLekIbLcpMwBA8hBwbCjRU0pMzfyZnabMAADJwxSVDSV6SompmWh2mTIDACQPAceGEj2lNDA1s2zjHhlSVMhx69TMwJQZAMCZmKKyoWRMKTE1AwBwE0ZwbChZU0pMzQAA3IKAY0PJnFJiagYA4AZMUdkUU0oAAMSOERwbY0oJAIDYEHBsjiklAADGjikqAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgODT6s0Bfv0l3YgAAkoiAk2I1TW2q3NKstnB35FrQl6WKsiLOlwIAIEGYokqhmqY2Ldu4JyrcSFIo3K1lG/eopqnNosoAAHAWAk6K9PWbqtzSLHOYzwauVW5pVl//cHcAAICxIOCkSH1Lx5CRm8FMSW3hbtW3dKSuKAAAHIqAkyLtXSOHm1juAwAAIyPgpEi+Nyuh9wEAgJERcFJkVmGegr4sjbQZ3NCJ3VSzCvNSWRYAAI5EwEmRzAxDFWVFkjQk5Ay8rygroh8OAAAJQMBJodLioNYvnq6AL3oaKuDL0vrF0+mDAwBAgtDoL8VKi4MqKQrQyRgAgCQi4CTQaI9gyMwwNOeCiRZUCACAOxBwEiTeIxg4nwoAgMQh4CTAwBEMJ/cgHjiC4XTrazifCgCAxGKRcZziPYKB86kAAEg8Ak6c4jmCgfOpAABIDgJOnOI5goHzqQAASA4CTpziOYKB86kAAEgOAk6c4jmCgfOpAABIDgJOnOI5goHzqQAASA4CTgLEegQD51MBAJAchmmartyi09nZKZ/Pp3A4rNzc3IQ8M9ZmffTBAQBgdEb7/U3ASWDAiQedjAEAOL3Rfn/TydgmOJ8KAIDESes1OD/96U9VWFiorKwszZgxQ2+88YbVJQEAABtI24Dz4osvqry8XA899JDefvttXX311Vq0aJE++ugjq0sDAAAWS9s1OLNnz9b06dO1fv36yLVLLrlEN954o6qrq0/783ZbgwMAAE5vtN/faTmC09vbq927d2vhwoVR1xcuXKgdO3YM+zM9PT3q7OyMegEAAGdKy4Dz6aefqq+vT36/P+q63+9XKBQa9meqq6vl8/kir8mTJ6eiVAAAYIG0DDgDDCN6G7VpmkOuDVi5cqXC4XDk1dramooSAQCABdJym/hZZ52lzMzMIaM17e3tQ0Z1Bng8Hnk8nlSUBwAALJaWIzhnnHGGZsyYodra2qjrtbW1mjt3rkVVAQAAu0jLERxJWrFihW6//XbNnDlTc+bM0YYNG/TRRx9p6dKlVpcGAAAslrYB55ZbbtHBgwf1yCOPqK2tTcXFxXr11Vd13nnnjernB3bHs5sKAID0MfC9fbouN2nbBydeH3/8MTupAABIU62trZo0adKIn7s24PT39+uTTz6R1+sdcecV4tfZ2anJkyertbWVhooW4O9vLf7+1uLvb61k/f1N01RXV5cKCgqUkTHyUuK0naKKV0ZGximTHxIrNzeXf8FYiL+/tfj7W4u/v7WS8ff3+XynvSctd1EBAACcCgEHAAA4DgEHSeXxeFRRUUGTRYvw97cWf39r8fe3ltV/f9cuMgYAAM7FCA4AAHAcAg4AAHAcAg4AAHAcAg4AAHAcAg4Srrq6WldeeaW8Xq/y8/N144036t1337W6LNeqrq6WYRgqLy+3uhTXOHDggBYvXqyJEycqJydHl19+uXbv3m11Wa5w/Phx/cM//IMKCwuVnZ2tKVOm6JFHHlF/f7/VpTnWtm3bVFZWpoKCAhmGoV/96ldRn5umqYcfflgFBQXKzs7WNddco3379iW9LgIOEq6urk7Lly/Xrl27VFtbq+PHj2vhwoU6cuSI1aW5TkNDgzZs2KBp06ZZXYprHDp0SFdddZXGjx+vrVu3qrm5Wf/8z/+sM8880+rSXGHVqlX693//d61bt05/+MMftHr1aj3++OP6t3/7N6tLc6wjR47osssu07p164b9fPXq1VqzZo3WrVunhoYGBQIBlZSUqKurK6l1sU0cSfe///u/ys/PV11dnebPn291Oa5x+PBhTZ8+XT/96U/16KOP6vLLL9fatWutLsvxfvCDH+jNN9/UG2+8YXUprvTVr35Vfr9fTz/9dOTa1772NeXk5Og///M/LazMHQzD0EsvvaQbb7xR0onRm4KCApWXl+uBBx6QJPX09Mjv92vVqlVasmRJ0mphBAdJFw6HJUl5eXkWV+Iuy5cv1/XXX6/rrrvO6lJc5ZVXXtHMmTP19a9/Xfn5+briiiv01FNPWV2Wa8ybN0+//vWv9cc//lGS9Pvf/17bt2/XV77yFYsrc6eWlhaFQiEtXLgwcs3j8WjBggXasWNHUn+3aw/bRGqYpqkVK1Zo3rx5Ki4utroc13jhhRe0e/duvfXWW1aX4joffPCB1q9frxUrVujBBx9UfX297r33Xnk8Hn3rW9+yujzHe+CBBxQOh3XxxRcrMzNTfX19euyxx/SNb3zD6tJcKRQKSZL8fn/Udb/frw8//DCpv5uAg6S65557tHfvXm3fvt3qUlyjtbVV9913n1577TVlZWVZXY7r9Pf3a+bMmaqqqpIkXXHFFdq3b5/Wr19PwEmBF198URs3btSmTZt06aWXqrGxUeXl5SooKNAdd9xhdXmuZRhG1HvTNIdcSzQCDpLmO9/5jl555RVt27ZNkyZNsroc19i9e7fa29s1Y8aMyLW+vj5t27ZN69atU09PjzIzMy2s0NmCwaCKioqirl1yySX65S9/aVFF7vL9739fP/jBD3TrrbdKkr74xS/qww8/VHV1NQHHAoFAQNKJkZxgMBi53t7ePmRUJ9FYg4OEM01T99xzjzZv3qzf/OY3KiwstLokV/nyl7+sd955R42NjZHXzJkz9c1vflONjY2EmyS76qqrhrRF+OMf/6jzzjvPoorc5fPPP1dGRvRXW2ZmJtvELVJYWKhAIKDa2trItd7eXtXV1Wnu3LlJ/d2M4CDhli9frk2bNunll1+W1+uNzMH6fD5lZ2dbXJ3zeb3eIeudJkyYoIkTJ7IOKgW++93vau7cuaqqqtLNN9+s+vp6bdiwQRs2bLC6NFcoKyvTY489pnPPPVeXXnqp3n77ba1Zs0Z33nmn1aU51uHDh/X+++9H3re0tKixsVF5eXk699xzVV5erqqqKk2dOlVTp05VVVWVcnJydNtttyW3MBNIMEnDvp599lmrS3OtBQsWmPfdd5/VZbjGli1bzOLiYtPj8ZgXX3yxuWHDBqtLco3Ozk7zvvvuM88991wzKyvLnDJlivnQQw+ZPT09VpfmWK+//vqw/86/4447TNM0zf7+frOiosIMBAKmx+Mx58+fb77zzjtJr4s+OAAAwHFYgwMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAPAMaqqqmQYxpDXmjVrrC4NQIpx2CYAx+jq6tKRI0ci7x955BG9+uqr2r59uyZNmmRhZQBSbZzVBQBAoni9Xnm9XklSZWWlXn31VdXV1RFuABdiigqA41RWVurZZ59VXV2dzjvvPKvLAWABAg4ARyHcAJAIOAAchHADYABrcAA4wqOPPqp169bpv/7rv+TxeBQKhSRJX/jCF+TxeCyuDkCqsYsKQNozTVNnnnmmOjs7h3y2a9cuzZ4924KqAFiJgAMAAByHNTgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBxCDgAAMBx/h8/O3w9lq8RugAAAABJRU5ErkJggg==",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using PyPlot, 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": [
    "Now build the factor graph in RxInfer, perform sum-product message passing and plot results (mean of posterior)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "using RxInfer, Random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior distribution of w: MvNormalWeightedMeanPrecision(\n",
      "xi: [318.19138919704756, 2358.5409231141653, 19126.67282976907]\n",
      "Λ: [15.000010000030002 81.69926298019712 580.4592797749924; 81.69926298019712 580.4592897750223 4589.625579040317; 580.4592797749924 4589.625579040317 38295.96922430482]\n",
      ")\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7rUlEQVR4nOzdeXicdb3//+fsmSQzk3WyJ03adN/3BWhZWihQQJayCLIcDwiiIkdRRK+Cl5YDHpHfkSMK56goAp6jrApIge6lzb6n2fc9mWQymWT2+/dHvx0JbSFtk06Svh/XNZfOPTP3/b5Hybz4rCpFURSEEEIIIaYRdagLEEIIIYQYbxJwhBBCCDHtSMARQgghxLQjAUcIIYQQ044EHCGEEEJMOxJwhBBCCDHtSMARQgghxLSjDXUBoRIIBGhvb8dkMqFSqUJdjhBCCCHGQFEUHA4HycnJqNWnbqc5bwNOe3s7aWlpoS5DCCGEEGegpaWF1NTUU75+3gYck8kEHPuCzGZziKsRQgghxFgMDg6SlpYW/B0/lfM24BzvljKbzRJwhBBCiCnmi4aXyCBjIYQQQkw7EnCEEEIIMe1IwBFCCCHEtCMBRwghhBDTjgQcIYQQQkw7EnCEEEIIMe1IwBFCCCHEtCMBRwghhBDTjgQcIYQQQkw7EnCEEEIIMe1IwBFCCCHEtCMBRwghhBDTjgQcIYQQQowrv99PT09PSGuQgCOEEEKIcVVSUsLhw4epq6sLWQ3akF1ZCCGEENNOXV0dLS0teL1eLBZLyOqQgCOEEEKIcdHT00NlZSXt7e3o9Xp8Pl/IapEuKiGEEEKcNafTSX5+Pn19fYyMjKDRaPB4PCGrR1pwhBBCCHFWfD4fubm52O12urq6AGhsbMTtdoesJmnBEUIIIcQZUxSFwsJCbDYbTU1NGAwGysrKOHr0KI2NjSGrSwKOEEIIIc5YdXU17e3t1NfXExkZSVlZGcPDw1itVuLj40NWlwQcIYQQQpyRjo4OqquraW5uxmg0cvToUfr7+4mLi+Oqq65i6dKlIatNxuAIIYQQ4rQ5HA4KCwvp6enB7/fT2tpKe3s7JpOJCy+8kMsuuyyk9UnAEUIIIcRp8Xg85OTkMDAwQG9vLy6Xi/r6esLCwli2bBk33XQTanVoO4mki0oIIYQQYxYIBMjPz6e/v5/m5mZUKhXl5eVotVpmz57NV77yFfR6fajLlIAjhBBCiLGrqKigq6uL+vp6IiIiKC4uRlEUUlNTufPOO4mKigLA6/WGtE4JOEIIIYQYk6amJurr60eFG5fLhdVq5bbbbiMjIwOA5uZmdu/ezcDAQMhqlYAjhBBCiC/U19dHaWkpbW1taLVaqqqqsNlsREdHs23bNpYtWwYc266hpKSEgYEBuru7Q1bvpAw4bW1t3H777cTGxhIeHs7SpUvJz88Pvq4oCo8//jjJyckYjUY2bdpEeXl5CCsWQgghpq/h4WHy8vLo7e3F6XTS0dFBa2trcMbUli1bABgcHCQvL4+qqioOHz5MWFhYyGqedAGnv7+fDRs2oNPpeO+996ioqODnP/95sE8P4Omnn+aZZ57hueeeIzc3l8TERDZv3ozD4Qhd4UIIIcQ0dHwbhv7+frq6unC5XFRVVWE0Glm6dCk333wzKpUKl8vFkSNHqKuro6ioiIGBgZA2Pky6aeJPPfUUaWlp/O53vwsemzFjRvC/K4rCs88+y2OPPcb1118PwEsvvURCQgKvvPIK991330nP63a7R+2JMTg4ODE3IIQQQkwTx7dh6O3tpaGhAZ1OR15eHnq9ntmzZ3P33XdjMBjw+Xzk5OTQ0NBATk4OPp+PhIQEEhMTQ1b7pGvBefvtt1m5ciU33XQTVquVZcuW8eKLLwZfb2hooLOzM9gcBmAwGNi4cSOHDh065XmffPJJLBZL8JGWljah9yGEEEJMddXV1bS1tVFXV0dkZCQFBQUApKamcs899xAVFYWiKBQUFFBfX8/BgwfxeDzExcWxZcsWFi5cGLLaJ13Aqa+v5/nnnyc7O5t//OMffO1rX+Ob3/wmf/jDHwDo7OwEICEhYdTnEhISgq+dzKOPPordbg8+WlpaJu4mhBBCiCmuvb2d6upqGhsbCQsLo6CgALfbjdVq5Stf+UqwoaCsrIza2lr27t2Ly+UiKiqKLVu2cNFFF6HT6UJW/6TrogoEAqxcuZKdO3cCsGzZMsrLy3n++ef5yle+EnyfSqUa9TlFUU449mkGgwGDwTAxRQshhBDTyMDAAIWFhbS3t6MoClVVVQwMDBAbG8sNN9zA4sWLgWONEpWVlXz88ccMDw8TGRnJ5s2bWb16NVVVVSxatChkKxpPuhacpKQk5s+fP+rYvHnzaG5uBgj25322taa7u/uEVh0hhBBCnB6Xy0Vubi59fX0MDAwEZ0yZzWYuu+wyLr74YuDY73BhYSF79uxhcHAQo9HIpZdeyoYNGygrK6O5uTmkg4wnXcDZsGEDVVVVo45VV1cHFw/KzMwkMTGRXbt2BV/3eDzs3buX9evXn9NahRBCiOnE7/eTk5ODzWajra0Nt9tNVVUVkZGRrFq1ihtvvBGVSkV/fz+5ubns37+fnp6e4FjYTZs2cfToUYaHh4mIiGD27Nkhu5dJ10X17W9/m/Xr17Nz5062b99OTk4OL7zwAi+88AJwrGvqoYceYufOnWRnZ5Odnc3OnTsJDw/ntttuC3H1QgghxNRVVFRET08P9fX1aLVacnJyCAsLY968edx9991otVqcTieHDx/mwIEDtLa2otPp2LBhA5deeikNDQ04HA7CwsJYu3ZtSIeGTLqAs2rVKt544w0effRRfvzjH5OZmcmzzz7Ll7/85eB7HnnkEUZGRnjggQfo7+9nzZo1fPDBB5hMphBWLoQQQkxd1dXVtLS0UFdXh9Fo5NChQ2i1WjIyMvjqV79KREQEHo+Hw4cPc+jQIWpra9Fqtaxdu5bLLruMjo4ObDYbOp2ONWvWEB4eHtL7USmKooS0ghAZHBzEYrFgt9sxm82hLkcIIYQImfb2dvLz86mtrUWtVpObm8vIyAhJSUk8+OCDzJ07F7/fzyeffMLu3bspLCxEq9WycuVKtm3bxsjICG1tbWg0GtauXUtMTMyE1TrW3+9JNwZHCCGEEOfO8RlTra2tKIpCeXk5Q0NDREdH8+Uvf5m5c+cGF/z75JNPKCoqQqVSsWjRIi6//HICgQBtbW2oVCpWrlw5oeHmdEjAEUIIIc5Tx2dM9fT0MDg4SHNzMx0dHVgsFq655hrWrVsHQEVFBQcPHuTIkSMoisL8+fO5/PLLMRqN1NfXA8eWdbFaraG8nVEk4AghhBDnoeMzpnp7e2lra8PhcFBXV4fZbGbjxo1s27YNOLbWzb59+zhw4AA+n49Zs2ZxxRVXEB8fH5z1vGDBAlJSUkJ5OyeQgCOEEEKcZ453OXV3d9PQ0ABAcXExERERLFu2jNtvvx2VSkVHRwcff/wxe/fuxe12k5mZydatW0lLS6OsrAyA7OxssrKyQnk7JyUBRwghhDjPVFVV0draSm1tLTqdjpycHIxGI9nZ2dx7773o9XpsNhsfffQRu3fvZmhoiPT0dK644gqys7MpLS0FICMjg7lz54b4bk5OAo4QQghxHmlpaaG6upq6ujoMBgNHjhxBp9ORkpLC/fffj9lsZmhoiI8++ohdu3YxMDBAUlISl19+OQsXLqSkpIRAIEBycjKLFi0K9e2ckgQcIYQQ4jzR19dHSUkJzc3NqFQq8vPz8fv9xMbGcu+995KcnIzb7Wb37t28//779PT0EB8fz9atW1mxYgWlpaX4/X6sVivLli373D0gQ23SLfQnhBBCiPHndDrJzc2lo6MDp9NJXV0dg4ODxMbGcscddzBv3jx8Ph/79+/nnXfeoaOjg5iYGC6//HLWrFlDWVkZXq+XmJgYVq5cGbJNNMdKAo4QQggxzXm9XnJycujp6aG7u5uenh7a29uJjo7mmmuu4YILLiAQCHD48GHeeOMNWltbsVgsXH755VxwwQVUVlbidrsxm82sXr0ajUYT6lv6QpM7fgkhhBDirAQCAfLy8uju7qa5uRmn00l1dTVms5mLLrqIa6+9FoCCggL+8pe/0NDQQEREBJs3b+aiiy6itraWkZERIiIiWLt2LTqdLsR3NDYScIQQQohprKysjPb2durq6ggEAhQVFREZGcny5cu56667UKvVVFZW8tprr1FTU0NYWBiXXnopmzZtorm5maGhoUmxeebpkoAjhBBCTFP19fXU19dTU1ODWq0mJyeH8PBw5syZw3333YdOp6OpqYk//vGPHD16FLVazcUXX8xll11GV1cXAwMD6PV61q1bF/LNM0+XBBwhhBBiGurs7KSsrIz6+no0Gg1HjhxBr9eTmprKgw8+iMlkorOzk9///veUlZWhKAoXX3wxl156KTabjd7eXrRaLWvWrCEyMjLUt3PaJOAIIYQQ08zAwAAFBQW0trbi8/nIz88HID4+ngcffBCr1YrNZuP3v/89BQUFuN1uLrjgAi677DJcLhddXV2o1WpWr15NVFRUaG/mDEnAEUIIIaaRkZERcnJy6OjowG63c/ToUYaHh4mKiuJf/uVfmDlzJg6Hgz/84Q8cPnyY4eFh1q9fz+WXX46iKMGdwVetWkVsbGyob+eMScARQgghpgmv18uRI0fo6uqio6ODtrY2urq6sFgs3HzzzaxatQqXy8Wf/vQn9u3bx9DQEKtWrWLr1q3o9frgAoDLly+fVDuDnwkJOEIIIcQ0EAgEyM/Pp6uri6amJux2O7W1tcH1bK644gq8Xi+vvvoqH374IXa7nSVLlnDVVVdhNpuDm24uXryY5OTkEN/N2ZOAI4QQQkwDZWVltLW1UVtbi9vtprS0FIvFwrp167j11lvx+/28/vrrvPfee9hsNhYsWMA111yD1WqlpqYGgAULFpCenh7iOxkfEnCEEEKIKa6uro76+nqqq6tRFIW8vDxMJhOLFy/m3nvvRa1W8/e//5033niD3t5eZs+ezTXXXENaWhpHjx4FYM6cOWRlZYX4TsaPBBwhhBBiCuvo6KC8vJza2loUReGTTz4hIiKCmTNn8uCDD2IwGPjoo4949dVX6erqIjMzky996UvMnDmTiooKALKyspg9e3aI72R8ScARQgghpqj+/n4KCwtpbGwM7jel1+tJSkriW9/6FhaLhU8++YSXXnqJrq4uUlNTueGGG5gzZw7l5eUoikJ6ejoLFiwI9a2MOwk4QgghxBTkdDrJycmhtbUVh8NBYWEhiqIQFxfHN7/5TZKSkigpKeHXv/41bW1tWK1WbrzxRhYsWBAMNykpKSxevDjUtzIhJOAIIYQQU4zH4+HIkSN0dHTQ1dVFdXU1TqeT6Oho7r33XrKzs6murubZZ5+lra2N6OhobrrpJpYsWUJ5eTmBQIDExESWLl2KSqUK9e1MCAk4QgghxBQSCATIzc2ls7OTlpYW2tvb6erqIjo6mltvvZVVq1bR2NjIz372M1paWggPD2f79u2sXLmSyspK/H4/8fHxrFixArV6+saA6XtnQgghxDSjKAqFhYXB6eADAwPU1dURHR3NlVdeyeWXX05XVxdPPfUUTU1N6PV6tm/fzpo1azh69Cg+n4+YmBhWrVo1rcMNSMARQgghpoyjR4/S1NRETU0NIyMjlJaWEh0dzQUXXMAtt9yCzWbjJz/5CXV1dajVam666SY2bNhAdXU1Ho+HqKgo1qxZg0ajCfWtTDgJOEIIIcQU0NTURFVVFVVVVbjdbvLz8zGbzSxbtox//dd/xel08tOf/jS4Fs4NN9zARRddRF1dHW63G7PZzNq1a9FqtaG+lXNCAo4QQggxyXV3d1NSUkJtbS0ej4ecnBxMJhNz587lG9/4Bn6/n5/+9KeUl5fj9Xr50pe+xMUXX0xjYyMjIyNERkaydu1adDpdqG/lnJGAI4QQQkxidrudvLw8GhoaGB4eJjc3F6PRSHp6Og8//DB6vZ4nn3ySkpIS3G431113HZdddhktLS04nU7Cw8NZt24dBoMh1LdyTknAEUIIISap4eFhjhw5QktLC/39/RQUFKDRaEhMTOTf/u3fMJvNPP300+Tn5+N0OrnyyivZsmULHR0dDA0NYTQaWbduHWFhYaG+lXNOAo4QQggxCXm9Xo4cOUJbWxtdXV1UVFTg9XqJjY3lW9/6FsnJyTz77LN88sknOBwOLr/8cq666iq6u7sZHBwkLCyMdevWER4eHupbCQkJOEIIIcQkEwgEyMnJoa2tjaamJurq6rDb7cTGxnLfffcxZ84cnnvuOfbu3cvg4CCXXHIJ11xzDTabDbvdjl6vZ926dURERIT6VkJGAo4QQggxiSiKQkFBAS0tLdTV1dHR0UFHRwcxMTHcdtttrF69mhdffJEPP/yQgYEBLrzwQq6//noGBwfp7+9Hp9Oxbt06IiMjQ30rISUBRwghhJhEKioqaGxspLq6mt7eXurq6oiJiWHbtm1cfvnlvPTSS/ztb3+jv7+fdevWsX37doaHh+nr60Or1bJu3TrMZnOobyPkJOAIIYQQk0R9fX1wrRu73U5lZSUxMTFcfPHF3HLLLfz5z3/mjTfeoL+/n5UrV3Lrrbficrno6elBo9Gwdu1aLBZLqG9jUpCAI4QQQkwC7e3tlJSUUFNTg8PhoLi4mKioKFavXs2//Mu/8MYbb/Daa69hs9lYvHgxX/7yl/F6vXR3d6PRaFizZg3R0dGhvo1JQwKOEEIIEWJ9fX0UFBRQX1+P3W4PrlK8cOFCHnzwQd5//31efvll+vr6mD9/PnfeeSeKotDV1YVarWbVqlXExsaG+jYmFQk4QgghRAg5HA5yc3Opr6+nr6+P/Px8IiMjmTVrFv/2b//GgQMH+N3vfkdPTw/Z2dncddddqNVqOjs7g+EmPj4+1Lcx6UjAEUIIIULE5XJx+PBhmpqa6O3tpaCgAIPBQGpqKo888gh5eXn85je/oaenh6ysLO6++270ej0dHR2o1WpWrlyJ1WoN9W1MSufHjltCCCHEJOP1ejl8+DAtLS20trZSUlKCWq0mMTGR7373uxw9epRf//rXdHV1kZaWxj333EN4eDjt7e2oVCpWrFhBQkJCqG9j0pIWHCGEEOIc8/v95OTk0NzcTGNjI5WVlXg8HqxWK9/+9rdpbW3lueeeo6Ojg+TkZO655x5MJlMw3CxfvpzExMRQ38akJgFHCCGEOIeOL+TX1NREbW0ttbW1DA4OEhcXxwMPPIDD4eC5556jvb2dpKQk/uVf/oWYmBja2tpQqVQsW7aM5OTkUN/GpDfpAs7jjz+OSqUa9fh0SlUUhccff5zk5GSMRiObNm2ivLw8hBULIYQQY1daWkpDQwPV1dU0NTXR3d1NXFwcd955JyqVil/+8pe0traSkJDA3XffTXx8PK2trQAsXbqUlJSUEN/B1DDpAg7AggULgktTd3R0UFpaGnzt6aef5plnnuG5554jNzeXxMRENm/ejMPhCGHFQgghxBerqamhurqayspKWltbaWlpIS4ujptvvpm4uDj+8z//k+bmZuLj47nrrrtITEykpaUFgCVLlpCamhriO5g6JuUgY61We9K+RUVRePbZZ3nssce4/vrrAXjppZdISEjglVde4b777jvlOd1uN263O/h8cHBw/AsXQgghTqG5uZnS0lKqqqro6Oigrq6O+Ph4rr76ambMmMGzzz5LU1MTMTEx3HnnnaSkpIwKN+np6SG+g6llUrbg1NTUkJycTGZmJrfccgv19fUANDQ00NnZyZYtW4LvNRgMbNy4kUOHDn3uOZ988kksFkvwkZaWNqH3IIQQQhzX1dVFYWEhNTU1dHZ2UlVVRWxsLJs2bWLJkiX853/+J/X19URHR3PHHXeQkZERDDeLFy+WcHMGJl3AWbNmDX/4wx/4xz/+wYsvvkhnZyfr16+nr6+Pzs5OgBOmxSUkJARfO5VHH30Uu90efBz/P44QQggxkWw2Gzk5OdTW1tLe3k5ZWRnR0dGsXbuWDRs28Mtf/pLa2lqioqK47bbbmDlzJs3NzcCxcJORkRHiO5iaJl0X1datW4P/fdGiRaxbt46ZM2fy0ksvsXbtWgBUKtWozyiKcsKxzzIYDBgMhvEvWAghhDgFh8NBTk4O9fX1tLW1UVpaSlRUFEuXLuXSSy/l+eefp6qqCovFwi233MKcOXOC4WbRokUSbs7CpGvB+ayIiAgWLVpETU1NcFzOZ1truru7ZbEjIYQQk8rIyAiHDx+mrq6OtrY2SkpKMJlMzJkzhyuvvJIXX3yRyspKLBYLN998M/Pnzx8VbmbMmBHaG5jiJn3AcbvdVFZWkpSURGZmJomJiezatSv4usfjYe/evaxfvz6EVQohhBD/5PF4+OSTT6ivr6e5uZmioiKMRiNZWVlce+21/P73v6e8vByLxcKNN97IokWLguFm4cKFEm7GwaTrovrOd77Dtm3bSE9Pp7u7m5/85CcMDg4G1wd46KGH2LlzJ9nZ2WRnZ7Nz507Cw8O57bbbQl26EEIIgc/n48iRIzQ2NtLQ0EBxcTE6nY7U1FS+9KUv8eqrr1JaWorZbOZLX/oSS5cupampCTgWbjIzM0N8B9PDpAs4ra2t3HrrrfT29hIfH8/atWs5fPhwsB/ykUceYWRkhAceeID+/n7WrFnDBx98gMlkCnHlQgghzneBQIC8vDwaGhqoqamhpKQElUpFcnIyN9xwA6+//jrFxcWYzWauu+46VqxYEQw3CxYskHAzjlSKoiihLiIUBgcHsVgs2O12zGZzqMsRQggxxSmKQmFhIZWVlZSXl1NWVobb7SYlJYUbbriBjz76iPz8fMxmM1dffTVr164d1S0l4WZsxvr7PenH4AghhBBTQXl5eXCV4srKSkZGRkhISODqq6/m448/pqCgALPZzNatW1mzZo2EmwkmAUcIIYQ4S8eDTUVFBUePHsXhcGC1WtmyZQuHDh0iPz8fk8nEli1bWLduXXAttkWLFkm4mSAScIQQQoiz0NDQQGlpKeXl5VRVVdHX10dcXBybNm2iqKiIvLw8TCYTl112GRs2bAhunClTwSfWpBtkLIQQQkwVbW1twXE3tbW1dHd3Ex8fz7p166iurh4Vbi644ALa2toAWaH4XJAWHCGEEOIMdHV1kZeXR1VVFdXV1bS2thIbG8uKFStobW0lNzcXs9kcDDft7e2AhJtzRVpwhBBCiNNks9k4cuQIVVVVVFVV0dTURFxcHAsXLqS3t5e8vDwsFgsbN25kw4YNwXAju4KfO9KCI4QQQpyGwcFBPvnkE2pqaqiqqqKhoYHY2FjmzJmD0+kMttxcdNFFbNiwgY6ODgCWLl0q4eYckhYcIYQQYoycTmcw3Bw9epSamhqio6OZMWMGXq832HJzPNx0d3ejUqlYunQpqampoS7/vCIBRwghhBgDl8vFJ598QnV1NRUVFVRXVxMdHU1qaioajSYYbi644ALWrl1LT08PKpWKZcuWkZKSEuryzzsScIQQQogv4PF4OHz4MDU1NVRUVFBVVYXZbMZqtWIwGILhZv369axbtw6bzYZKpWLFihUkJSWFuvzzkgQcIYQQ4nN4vV4OHz5MdXU1ZWVlVFZWEhkZSUxMDEajkcLCQqKjo1m3bh1r167FZrOhVqtZsWIFiYmJoS7/vCUBRwghhDgFv99PTk4ONTU1lJaWUlFRQXh4OBaLBbPZTGlp6ahwMzAwgFqtZuXKlSQkJIS6/POaBBwhhBDiJI7vDF5dXU1RcQmH8opQVFoI0zEjOpqysjJiYmJYt24dq1atCoabVatWYbVaQ13+eU8CjhBCCPEZx3cGr6qq4q2PD/HxoTw8AVBrw1CNGMmvOsjimals2LCBZcuWMTg4iEajYdWqVcTHx4e6fIGsgyOEEEKMoigKJSUlVFZW8saHB/nHvlw8flCpdRBuxm9rxaePoNybgNucytDQEFqtlrVr10q4mUQk4AghhBCfUlFRQVlZGQUFhXx0MBcFBdQaiLAQGOhAbbRgSJ6PPmEGL+2rRq05Fm5iYmJCXbr4FAk4QgghxP9TVVVFSUkJRUVF7M/Jx+3zo6g0qMLNKPYu1EYL+uQ56OLSUPxe+l0B9MlziY6ODnXp4jNkDI4QQggB1NfXU1xcHHyMuDyACrXRhOLoRR0ehSF5DtqYVAj4UGl0GJLnMow+1KWLk5CAI4QQ4rzX1NREfn4+xcXFFBYW4na70es0qIxmlOEBNOHR6JPnoI1JRqX4UWn0GFLmotYbsZrCQl2+OAnpohJCCHFea21tJTc3l+LiYvLz83G5XKhUKjJTEjF4h9BERKFPmYsuJgmVEkClNWBImYdGbyTJEsbqTBl7MxlJC44QQojzVnt7Ozk5OZSUlJCfn4/T6USlUhEbG4t9oJ/lc2dQ7LGitVhBUVDpwjAkz0WjMwCwY9t8NGpViO9CnIwEHCGEEOelrq4uDh8+THFxMbm5uTgcDgCio6NxOBwkJSWxdu1aLvKH8+e8NgZ8GgzJc1Fp9SRawtixbT5XLJR9piYrCThCCCHOOz09PXzyySeUlpaSk5OD3W4HICYmhpGRkWC4SU5OJl2jYdPiLDQJs+l3B7CajnVLScvN5CYBRwghxHnFZrPxySefUFJSwpEjRxgYGAAgKioKj8dDcnIya9asISkpCY1GQ0xMDKtXr0an04W2cHFaJOAIIYQ4bwwMDHDo0CGKi4s5cuQINpuNQCBAVFQUiqKQkpLC6tWrSUhIQKvVEhcXx6pVq9Bq5edyqpH/xYQQQkw7/oBCToONbocr2KXkHHJw8OBBSkpKOHz4ML29vQQCAcxmMxqNhpSUFFauXInVakWn05GQkMDKlStRq2XC8VQkAUcIIcS08n5ZB0+8U0GH3RU8Zg0LcHWCHd1AC4cOHaKnpwe/309kZCQGg4HU1FRWrlxJbGwsOp2O5ORkli1bJuFmCpOAI4QQYtp4v6yD+18uQPnUsYBnhKb6Cp75uJ7ZqnYCzmPdUhEREURERJCens6yZcuIjY1Fr9eTlpbGkiVLUKlkEPFUJgFHCCHEtOAPKDzxTsVnwo0LV1sFnq4G3G0VlLv6mRUfQWREBBaLhYyMDBYvXkxcXBx6vZ7MzEwWLFgg4WYakIAjhBBiWshpsI3qlgp4XbjbKvB0NuBqLcc/1AuKgl+lIyYmhhkzZrBgwQLi4+PR6XTMnj2bOXPmhPAOxHiSgCOEEGJa6HZ8Oty4cbdV4uluwNVaht/RC0oAlTYMozmarKws5s6dGxxQvGDBArKyskJYvRhvEnCEEEJMC8c3vVR8nmPhpquRkeYS/IN9EAiAPgx1ZBRZM7OYPXs2CQkJ6PV6Fi9eTHp6eoirF+NNhocLIYSYFlZnxpAQocbdWomnq4GR5mJ89h4I+EBnQBsZQ3RSBhuWLSQ5ORmDwcDy5csl3ExTEnCEEEJMCz6vh+tTnLi763E1F+Gzd6NS/KALQ2eOQxeTwpcuWEJaagoGg4FVq1aRnJwc6rLFBJEuKiGEEJPayRbt++w+UB6Ph0OHDqG3tzDTU0f5YA/4faAPR2exEp2cwbXrFrBl9QIMBgNr1qwhJiYmRHckzgUJOEIIISatky3al/SZnbyPh5vCwkI++OADhm1dpEbp0RpjiLYmM3f+PJbMTCcjI53w8HDWrFmDxWIJ1S2Jc0QCjhBCiEnpZIv2AXTaXdz/cgHP376cS+fEcejQIQoKCnj//fdpaWnB4/EQExNDZmYmCxcuJCoqioyMDCIiIli3bh0REREhuR9xbknAEUIIMemcbNG+4xRABex4swTDBgOFhQW8++67tLa24vV6iYuLY8aMGSxevBiz2Ux6ejoWi4W1a9cSFhZ2ju9EhMpZB5z9+/czb9484uLixqMeIYQQ4oRF+z4r4PfRWFHK68Neyg/vpq2tDY/HQ0JCApmZmSxatIjIyEgyMjKIjo5mzZo16PX6c3gHItTOehbVHXfcQUVFBQC33347fX19AAwNDREIBM729EIIIc5Dn16077MUvw9X21HcHTV8/MGxlhu3201CQgIzZ85k8ZIl9LhUdAVMdLr1rF6zVsLNeeisA05nZyczZswA4K233sJutwNw+PBh1q5de7anF0IIcR46vmjfZ/0z3FQzUnMYe08nbreb5ORksrOzMSRk8cKhNl4uH+G3hXZ+mudl08/38X5Zxzm+AxFqZx1w4uPjsdlsJxyfO3cuNTU1Z3t6IYQQ56HVmTEkWcL49GTwYLhpr2ak+hAM96NS/KSkpDB79mx0cem8VdaHUx+NzhKPPnEWKpU6OChZQs755awDzrZt23j55ZdPON7a2orH4znb0/Pkk0+iUql46KGHgscUReHxxx8nOTkZo9HIpk2bKC8vP+trCSGEmBw0ahU7ts0Hjg0o/me4Ocpw9UF8zn7MeshIT2fevHnMmTuX3Q1ONOY4dNHJ6K1ZwR3Bjw9UfuKdCvyBkw1bFtPRWQecxx9/nDfffJOf/OQnKMqx/+PU1dXxne98h1WrVp3VuXNzc3nhhRdYvHjxqONPP/00zzzzDM899xy5ubkkJiayefNmHA7HWV1PCCHE5HHFwiSev3051kgtrrZK3G2VDFcdRDUyQLRBzdzsmSxYsIDs7GzsfgNOrRldXAa62LQTzqUAHXYXOQ0n9jiI6emsA47VauXIkSMUFxfj9XrJzs5m9uzZVFVV8cwzz5zxeYeGhvjyl7/Miy++SHR0dPC4oig8++yzPPbYY1x//fUsXLiQl156ieHhYV555ZVTns/tdjM4ODjqIYQQYnK7dE4c/77BwHWJDmK784lRDRMfoWX+nFksWrSIjIwMrFYrWlMM+oSZ6KISP/d8nzd4WUwv47IXVWxsLP/3f/9HR0cH77zzDn/729+or69n+fLlZ3zOr3/961x11VVcdtllo443NDTQ2dnJli1bgscMBgMbN27k0KFDpzzfk08+icViCT7S0k5M+EIIISYPr9fLwYMH2bdvLzm738U3PIhOo2LWrFksW7aM1NRUkpOTSUxMZN3aNWhNX7xcyakGL4vpZ1wX+ouJieHKK6886/O89tpr5Ofnk5eXd8JrnZ2dACQkJIw6npCQQFNT0ynP+eijj/Lwww8Hnw8ODkrIEUKIScrj8XDw4EH27NnDm2++SV9fHyrVsXCzYsUKoqKiSE1NxWq1smbNGsyWKJI+aqPT7jrp4oAqINFybB8rcX6YdCsZt7S08K1vfYsPPvjgc1ecPD547DhFUU449mkGgwGDwTBudQohhJgYx8PNBx98wN/+9jf6+vrQarVkZ2ezevVqIiIimDFjBomJiaxdu5bIyEgAdmybz/0vFxwblPyp8x3/Zdixbf4Jm3SK6WtcuqjGU35+Pt3d3axYsQKtVotWq2Xv3r3853/+J1qtNthyc7wl57ju7u4TWnWEEEJMLW63m4MHD/Luu+/y1ltv0dfXh06nY/bs2axfv56IiAhmzpxJSkoKGzZsCIYb+Oeg5ETL6H85TrSE8fzty4Obc4rzw7i24Bw8eJCVK1eeVUvJpZdeSmlp6ahjd999N3PnzuV73/seWVlZJCYmsmvXLpYtWwYcS/t79+7lqaeeOqv6hRBCjD9/QCGnwUa3w4XVdKyb6GQtKS6XiwMHDvC3v/2N999/n4GBAYxGI7NmzeKCCy5ArVaTnZ1Namoqq1evRqfTnXCOKxYmsXl+4piuJ6a3cQ04W7dupaioiKysrDM+h8lkYuHChaOORUREEBsbGzz+0EMPsXPnTrKzs8nOzmbnzp2Eh4dz2223nVX9Qgghxtf7ZR088U7FqH2lkixh7Ng2f1SLisvlYv/+/fz1r39lz5499Pf3ExERwZw5c1i/fj1wbAHZGTNmsHz5cjQazSmvqVGrWDczduJuSkwJ4xpwjq+DM9EeeeQRRkZGeOCBB+jv72fNmjV88MEHmEymc3J9IYQQX+z9sg7uf7nghEG/x1cWPt5tNDIywt69e/nf//1f9u/fz8DAAGazmblz57J27VrUajVz584lOzubRYsWfe54SyGOUynjmEpMJhPFxcVn1YJzrgwODmKxWLDb7ZjN5lCXI4QQ04o/oHDBUx+fckfw47OaPvjGWvbt3cOrr77KJ598wuDgIFFRUcyfP5/Vq1ejVquZM2cOCxcuZPbs2ef2JsSkNNbf77NqwfnDH/4w6rnP5+P111/HarUGj33lK185m0sIIYSYgnIabKcMN3BsllNbzwD/9fJfyf/oHXJzc3E6nURHR7No0SKWLVuGTqdjzpw5LFu2jPT09HNXvJgWzirg/O53vxv13Ov18pe//AWj0Qgcm8otAUcIIc4/X7RicMAzwkh9AW8UFtBRV8HIyAixsbEsXbqUBQsWYDQag9PCExM/f3ViIU7mrALO7t27Rz03mUy88sorU6KLSgghxMT5vBWDA+5hhmtzGKrcy9GBJgxqhfj4eFasWEF2djZms5nZs2ezZs0aYmJkYT5xZibdQn9CCCGmvtWZMSRZwk5YWTjgGmKo+hOclfvwdTfiUQJEJVnZsGEDKSkpxMbGMmfOnFEL+AlxJibdQn9CCCGmPo1axY5t80cd8484cFTsYahsN76uegKKH3W4mZG4uSQmJpGUlMSiRYu48MILJdyIszauAecHP/iBNCcKIYQAji2699Blx2Y++ZwDOEp24azYh7+3kQAqNBHRGGcsxa2NwB8Rx9KlS1m/fr1sqyPGxbh2UT366KPjeTohhBBT3Iy4cHyOPhwluxipy8XX3w5qDdqIaMIylqDWh6GNSsKaOY9Vq1ahVkvHghgfMgZHCCHEhAkM9uAofI+R+ly8A12otHo0EVGEZSxDo9OjjUnGkDSbNSuXS7gR40oCjhBCiAlRU1PDwXdewVt/BO9Az7FwExlLWMYSNDo9uthUwlLmkpqWzupMGd4gxpcEHCGEEOOuuKSUnb/4Lw7u/Qicfai1BjQWK8b0hai0BnRx6RhT56OJiGLHtvmyGaYYd+MScLxeL52dnQwPDxMfHy8DjYUQ4jz2i5ff4clfvoCtKg//yBBqfRiG2GQiM+bjVenRxWcSlraAFGvMCZtuCjFezjjgDA0N8ac//YlXX32VnJwc3G538LXU1FS2bNnCvffey6pVq8alUCGEEJPfzhde46fP/Q+uxmICHhdqgxFtTAq6pGy8GiO3bt3Ipos2kBYfzerMGGm5ERPmjEZ0/eIXv2DGjBm8+OKLXHLJJbz++usUFRVRVVXFJ598wo4dO/D5fGzevJkrrriCmpqa8a5bCCHEJBIIBPjLX//Kz3/134zUFxDwuVAbwtDFZRCWNBt1uIWw1HnkexK5YfVM1s2MlXAjJtQZteAcOnSI3bt3s2jRopO+vnr1au655x5+/etf8z//8z/s3buX7OzssypUCCHE5OT1enn11Vf5/Wuv01+dh0IAtS4MQ+JstLEpaCKiMaTMR5+QSdeQj5wGG+tmxoa6bDHNnVHA+b//+78xvc9gMPDAAw+cySWEEEJMAcPDw/zxj3/k3XffJe+TA4CCWh+OIXke2qgEdKZjs6Z0MSnBz3zRRpxCjAeZRSWEEOIE/oBCToONbocLqynspONl+vv7eemll3j33XcpKChAp1GjMkRgSJ2PzhSHLjqJsIylaM1xoz73eRtxCjFezijgvP3226f9mc2bN2M0Gs/kckIIIc6h98s6eOKdCjrs/2xpSbKEjZrx1NbWxp/+9CfeeecdKioq0Ol0REVFoTHNxKM3o4tNIyxzGRqjOXgOFZBoCZM1b8Q5oVIURfnit412uqtNqlQqampqyMrKOt1LTZjBwUEsFgt2ux2z2fzFHxBCiPPA+2Ud3P9yAZ/9YTjedvP87cuZobXzv//7v7z99ts0NDRgMBiIiYnhoosuwqEK54OucMLSF6PSG0/6eZkWLs7GWH+/z7iLqrOzE6vVOqb3mkymM72MEEKIc8QfUHjinYoTwg1w7Jii8N1fvc7VUe28/dabdHZ2YjAYsFqtXHjhhSQkJLBu3Tq2x8xi5/u1o1qAEj/TAiTERDujgHPnnXeeVnfT7bffLq0kQggxyeU02EaFkk9TAn5cLWX01eXyh64i3E4HBoOB1NRUNmzYQGJiIps2bQpumHnlkvQvHMMjxEQ6oy6q6UC6qIQQYrS3itr41mtFJxwPeN24mopx1uXiqssjQuMj0hhGVlYWa9euJSkpia1btzJ//vxzX7Q474z19/u0F/obGRmhra3thOPl5eWneyohhBCTyMlmNwVcQ4zU5jB0dD+umiMoPhdOvxprahYbNmxgxowZ3HTTTRJuxKRzWgHnL3/5C7Nnz+bKK69k8eLFHDlyJPjaHXfcMe7FCSGEOHdWZ8aQZAkLDgj2Ofpw1uYwVLkPd0MhSsCHotaji8ukQZ2MJyKBW2+9lfT09JDWLcTJnFbA+clPfkJBQQHFxcX89re/5Z577uGVV14B4Dzt6RJCiGlDo1axY9uxlhivrY3h+lyGynbjbqtAUQKgNRCWNAtDQhb6hEwOKrOJipYp32JyOq1Bxl6vl/j4eABWrlzJvn37uP7666mtrUWlksFjQggx1W2eZ+Xby7Q8+fsihkr24OtvBxRUujD0SdnoY9MwpMzFmLGU7hGVbLsgJq3TasGxWq2UlJQEn8fGxrJr1y4qKytHHRdCCDH1uFwunvnj2/zst3+hJ/8DfP1tKAqoDJGEpc3HEJ+JMWslxszlqDTH/v1Ytl0Qk9VpteD88Y9/RKsd/RG9Xs+rr77Kgw8+OK6FCSGEOHfsdju/fPXvPP2Ht3BWHUIZcaAAmnATYakL0EUlEJ69Dl1c2qjPybYLYrI6rRac1NRUEhMTRx179tln6ejoYMOGDeNamBBCiHOjs7OT3Xv28twf/4/hyv0oI4MoKtBGxhKWvhh9XBoRCy8dFW5UHNu+QbZdEJPVaU8T/6yHH36YDRs20NraOuq4x+MhNzf3bE8vhBBiAtXW1rJ3336e/s1L9JQeIuAeApUanSWBsPSFGBJmErFoM1rTP8fZHB9xuWPbfFm8T0xa47Kb+NatW7nooovYt28fqampwLFdZteuXYvf7x+PSwghhBhHfr+f4uJi/vzhYX7z+1dxtB5F8XlArUEbnYwhcTaGlDlEzFqLSqsb9VnZdkFMBWcdcFQqFTt27MBqtZ4QcmTquBBCTD4ul4vc3Fze2JvHf734R3y2FgJ+N2q1Dl1cOvrEbIwzlhCWsRiVanRD/4+umsddGzKl5UZMeuPSggOwY8cOgGDI0el0MnVcCCEmmYGBAXJzc6moPMoLL/4PXlsnis+DWmtAb81Cn5BFxJz16BOyRn1OxbGWGwk3Yqo464Dz6VaaT4ecP//5z2d7aiGEEOOora2NoqIiioqK+ONr/8dwXzuK34dab8SQmI0+cRaR8zeitVhHfU7G3Iip6KwDzk9/+lMiIiKCz4+HnKuuuupsTy2EEGIcKIpCVVUVVVVV7N+/nw8++ICujk6UgB91mImw5Dnok2cTuWATmjDTCZ+XMTdiKjrrgPPoo4+ecGzHjh1oNBr+4z/+42xPL4QQ4iz4fD4KCwtpbm7mvffe49ChQ/T19aEEAmgjojAkzyMsbSHhczeg1upP+LyMuRFTlUo5T0cCj3W7dSGEmKqcTie5ubm0t7fz5ptvUlRURH9/P2q1mtTUNFr1KajSVmLIXHrCYOLjY24OfO8SCTdiUhnr7/e4DTIWQggxefT09JCfn09LSwtvvfUWlZWVDA4OotFoyM7OZu3atWSs3szP870AfPrfdGXMjZgOJOAIIcQ0U19fT3l5OVVVVbz//vvU1tYyNDSEXq9nwYIFbNiwgbvvvpuMjAxmzevgiXcq6LD/c08pGXMjpgMJOEIIMU34/X5KSkpobGykoKCAPXv20NjYyPDwMOHh4SxfvpyLL76Yu+66i6ioKACuWJjE5vmJ5DTY6Ha4sJqObb8gLTdiqpvQgKNWq9m0aRM/+9nPWLFixUReSgghzmvHF+/r6OjgwIEDHD58mNbWVjweDxaLhZUrV3L11Vdz6623otONXplYo1axbmbsKc4sxNQ0oQHnt7/9LU1NTXzzm9/k4MGDE3kpIYQ4b9lsNvLy8ujs7GTPnj0UFhbS0dGBz+cjPj6e1atX8+Uvf5ktW7aEulQhzhmZRSWzqIQQU1hTUxOlpaU0NTWxb98+ysrK6O3tJRAIkJqayvr167n//vtZuHBhqEsVYlyM9ff7rHcTf/bZZ2lvbz/b0wQ9//zzLF68GLPZjNlsZt26dbz33nvB1xVF4fHHHyc5ORmj0cimTZsoLy8ft+sLIcRUEAgEKC4uprCwkLKyMnbt2kVhYSHd3d0AZGdnc/XVV/PDH/5Qwo04L511wHn44Ye58MILaW1tHXXc4/GQm5t72udLTU3l3//938nLyyMvL49LLrmEa6+9Nhhinn76aZ555hmee+45cnNzSUxMZPPmzTgcjrO9FSGEmBJcLheHDh2iqqqKnJwc9u7dS3FxMX19feh0OhYtWsTNN9/M97//fZKSZCaUOD+ddReVWq3m/vvv57333hu1k3hXVxfJycn4/f6zLjImJoaf/exn3HPPPSQnJ/PQQw/xve99DwC3201CQgJPPfUU99133ynP4Xa7cbvdweeDg4OkpaVJF5UQYkr59HibgoICiouLqa2tZXh4mIiICJYuXcrtt9/ODTfcgEajCXW5Qoy7c9ZFpVKp2LFjB3feeScXXXTRqJacsx3e4/f7ee2113A6naxbt46GhgY6OztHDZQzGAxs3LiRQ4cOfe65nnzySSwWS/CRlpZ2VrUJIcS51tjYyMGDB6mtreXQoUPk5ORQVVWF0+kkOjqaDRs28Mgjj7B9+3YJN+K8N26zqD69k/i+ffvQ6XSoVGe2jkJpaSnr1q3D5XIRGRnJG2+8wfz584MhJiEhYdT7ExISaGpq+txzPvroozz88MPB58dbcIQQYrILBAKUlJTQ0NBATU0NR48eDc6U8nq9JCYmctFFF/HII4+QlZUV6nKFmBTOOuB8upXm0yHnz3/+8xmfc86cORQVFTEwMMBf//pX7rzzTvbu3Rt8/bPBSVGULwxTBoMBg8FwxjUJIUQoDA8Pk5eXR0dHB+Xl5dTW1lJYWEhvby8AM2bM4Oqrr+bb3/52cPE+IcQ4BJyf/vSnREREBJ8fDzlXXXXVGZ9Tr9cza9YsAFauXElubi7/3//3/wXH3XR2do4aONfd3X1Cq44QQkx1x/eT6urqory8nLq6OgoLC7Hb7eh0OmbPns0dd9zBXXfddcLifUKc78464Dz66KMnHNuxYwcajYb/+I//ONvTA8daaNxuN5mZmSQmJrJr1y6WLVsGHJuttXfvXp566qlxuZYQQoSaoijU1tZSWVlJS0sLdXV1VFZWUl5ejtPpxGg0snjxYr71rW+xZcuWMx4OIMR0NmErGf/whz/khz/84Wl/7gc/+AFbt24lLS0Nh8PBa6+9xp49e3j//fdRqVQ89NBD7Ny5k+zsbLKzs9m5cyfh4eHcdtttE3AXQghxbnm9XoqKimhpaaG6upr29vbgTKmRkRGioqJYu3Ytjz32GAsWLAh1uUJMWpNus82uri7uuOMOOjo6sFgsLF68mPfff5/NmzcD8MgjjzAyMsIDDzxAf38/a9as4YMPPsBkMoW4ciGEODv9A3Ze/tvHNDS3MdDZhN4zRF5eLi0tLfh8PhISEti6dSs/+MEPiIuLC3W5QkxqslWDrIMjhJgEXv6ogJ/88X16ujrwDnThc9gIdFWj9w1h0GpIS0vjy1/+Mt/4xjfQ6/WhLleIkBnr7/eka8ERQojzSSAQ4Ddv7mHHnz4+FmyG+vH2d+LtrkVxOfGqtaTPnMmPHvse1113nYy3EWKMznqhv889uVrNJZdcQn5+/kReRgghpqSRkRH27T/AL/5vD+6eJvwOG57uBrydR1FcQ6h0Bgwpc1GtvYtrrpVwI8TpmNCA89vf/paNGzfyzW9+cyIvI4QQU05PTw/79u3jo5wSuppqCIwM4Wo/ire7HsXjQhUWiXHWGqIvfwBHZBo5DbZQlyzElDKhXVR33XUX8M+1cYQQ4nynKAo1NTVUVlbS2trK0aoafCMOPG0V+Id6QQF1ZAyRCy7FvHIbat2xBUq7Ha4QVy7E1CJjcIQQ4hzxeDwUFBTQ1tZGbW0tQ0NDOPu7cTcWEhgZRKVWo4lOxLzqesLnrEet/mcju9UUFsLKhZh6zirgtLa28vzzz3Po0CE6OztRqVQkJCSwfv16vva1r8leT0II8f8MDAwEdwGvqanB7/dTXl5OaV4+imsQlVaHzjqTqAtvIyxpdvBzKiDREsbqzJjQFS/EFHTGAefAgQPBBfm2bNnCli1bUBSF7u5u3nzzTX75y1/y3nvvsWHDhvGsVwghppyGhgbKy8tpbW2lra0Nn8/HkSNHqKysxO12ExcdhdO6kOgLbkUTGRv83PEhxTu2zUejlgHGQpyOM14HZ9WqVVxwwQX84he/OOnr3/72tzlw4AC5ublnVeBEkXVwhBATzev1UlxcTHNzM3V1dTidTux2O/v376epqYlAIIDVauXuu+8m++IbePrDRjrs/xxrk2QJY8e2+VyxMOlzriLE+WWsv99nHHCMRiNFRUXMmTPnpK8fPXqUZcuWMTIyciann3AScIQQE8lut5Ofn09HRwf19fUANDc3s3//fnp6elCr1cycOZPvf//7XHPNNajVavwBhZwGG90OF1bTsW4pabkRYrQJX+gvKSmJQ4cOnTLgfPLJJ6N2/BZCiPNFU1MTpaWltLW10dHRAUBlZSWHDx/G4XBgMBhYs2YNP/7xj1m8eHHwcxq1inUzY091WiHEaTjjgPOd73yHr33ta+Tn57N582YSEhJQqVR0dnaya9cu/vu//5tnn312HEsVQojJzefzUVJSQmNjI/X19YyMjOD1ejly5AhlZWV4vV5MJhPXXHMNP/jBD+RfAoWYQGcccB544AFiY2P5xS9+wW9+8xv8fj8AGo2GFStW8Ic//IHt27ePW6FCCDGZDQ4OBruk6urqUKvV9Pf3s3fvXpqbm4Pjbb72ta9x7733YjQaQ12yENPauGy26fV66e3tBSAuLg6dTnfWhU00GYMjhBgvn+2S0mq1NDQ0sHfvXvr7+1Gr1cydO5fvf//7XHHFFWg0mlCXLMSUdU4329TpdNLUKoQ473y6S6qurg632w1AYWEhhw8fxuVyYTAY2LhxIz/60Y+YP39+iCsW4vwxYSsZX3bZZdTX1wdnDwghxHTy2VlSer2e4eFhDhw4QFVVFT6fD4vFwq233sq3v/1t4uLiQl2yEOeVCQs4X/rSl4LdVkIIMZ0cX7ivpaWF7u5uwsLC6Ojo4OOPP6a9vR2A9PR0/u3f/o0bb7yRsDDZZkGIc23CAs7Xv/71iTq1EEKEhNfrpaioKLhwn9/vR6PRcPToUT7++GOGhobQarUsX76cxx57jPXrR+8nJYQ4d2SzTSHEtDTei+bZbDYKCgpob2+nsbGRiIgIhoaGyMvLo7CwEI/Hg9FoZNu2bXz3u98lKytrHO9GCHG6JizgtLS0sGPHDn77299O1CWEEOKk3i/r4Il3KsZl2wNFUaitraWyspKWlhb6+vowGo10dnaye/dumpqaUBQFq9XKfffdx1133UVUVNQ435EQ4nSNyzTxkykuLmb58uXB9XEmG5kmLsT09H5ZB/e/XMBn/7Adb7t5/vblYw45LpeLwsJCWlpaqKurQ6PR4PV6aW5u5sMPP2RgYACNRsOCBQv4zne+w+bNm9Hr9eN6P0KI0SZ8mvjbb7/9ua/L7CkhxLnmDyg88U7FCeEGQOFYyHninQo2z0/8wu6qrq4uioqKaGtro6WlBbPZzODgICUlJXzyySd4PB7CwsK47LLL+O53v8uCBQtQqWTfKCEmizMOONdddx0qlYrPawCSf9iFEOdSToNtVLfUZylAh91FToPtlHs+BQIBKisrqa6upqmpCafTidlspqenhz179gQHF8fExHDXXXfx1a9+FavVOkF3JIQ4U2c8vD8pKYm//vWvBAKBkz4KCgrGs04hhPhC3Y5Th5uxvG9oaIgDBw5QUlJCeXk5iqKgUqmoq6vjL3/5CzU1NQQCAWbNmsVPfvITHnroIQk3QkxSZ9yCs2LFCgoKCrjuuutO+voXte4IIcR4s5rGtt7Myd7X1NREWVkZra2tdHV1ER0djc1mo6ysjEOHDuF2u9Hr9Vx44YU8/PDDrFq1SqaACzGJnXHA+e53v4vT6Tzl67NmzWL37t1nenohhDhtqzNjSLKE0Wl3nXQcjgpItBybMn6cx+OhuLiYpqYm6uvrURQFs9lMR0dHsEsqEAgQFRXFLbfcwr333kt6evo5uychxJmZsFlUk53MohJiejo+iwoYFXJONouqt7eXwsJC2traaG5uJiYmBofDQVtbG7t27cJms6FWq8nMzOSb3/wm27Ztw2QyndsbEkKMMtbf7zNqXy0pKSEQCIz5/eXl5fh8vjO5lBBCnJYrFibx/O3LSbSM7oZKtIQFw83xgcQHDhygvLyc1tZW4uLigov5vfXWW9hsNvR6PRs3buTnP/8527dvl3AjxBRyRi04Go2Gzs5O4uPjx/R+s9lMUVHRpFrZU1pwhJjeTrWSsdPppKCggJaWFurr64mMjCQQCNDT08PevXupr6/H7/cTHR3NTTfdxFe/+lUyMzNDfTtCiP9nQtfBURSFH/3oR4SHh4/p/R6P50wuI4QQZ0yjVp0wFby5uZnS0lKam5vp7u4mISEBm81Gc3MzH3/8Mf39/Wg0GmbOnMn999/Pddddh8ViCdEdCCHOxhkFnIsuuoiqqqoxv3/dunUYjcYzuZQQQpy14wOJGxsbqa+vR61WY7Va6ejooLCwkIKCAtxuN2FhYWzYsIEHHniAdevWodPpQl26EOIMnVHA2bNnzziXIYQQE6O7u5uioiJaWlpoaWkhISGB4eFhGhoa2L17Ny0tLSiKQlxcHNu3b+crX/kKM2fODHXZQoizJLuJCyGmJb/fT0VFBTU1NTQ0NOB2u0lLS6Ozs5Pa2lr279/P4OAgGo2G2bNnc99993HVVVfJRplCTBNnFHCam5tPax2ItrY2UlJSzuRSQghx2ux2e3AgcUNDA7GxsRgMBpqamjhy5AhlZWX4fD7Cw8PZuHEj9957L6tXr5YuKSGmkTOaJr5q1Sr+9V//lZycnFO+x2638+KLL7Jw4UJef/31My5QCCHGSlEUampq2LNnD6WlpTQ2NpKWlsbw8DDV1dW8+eabFBcX4/f7SUhI4Ktf/SqPP/44GzZskHAjxDRzRi041157LSaTiSuuuAKdTsfKlStJTk4mLCyM/v5+KioqKC8vZ+XKlfzsZz9j69at4123EEKM4nQ6KSwspLm5mfr6esLDw0lKSqKlpYXy8nJycnJwuVzodDrmzZvH3XffzVVXXSXLRAgxTZ3ROjh6vZ6WlhbMZjMJCQls376dvr4+RkZGiIuLY9myZVx++eUsXLhwImoeF7IOjhDTR2NjI2VlZTQ1NdHT00NqairDw8O0trZy8OBBGhoaCAQCmM1mLr30Uu666y5WrFiBVivDEIWYaiZ0HZyUlBQKCwu54oorGBoaYufOnbKjrhDinHO5XBQVFQWnf+t0OjIzM2lra6O2tpaDBw9it9vRaDTMmDGDm2++mZtuuonU1NRQly6EmGBnFHC+853vcM0117By5UpUKhV/+tOfuOCCC1i4cKGsdyOEOCfa2tooLi6mubmZjo4OUlJSUBSF6upq8vPzKS8vx+v1YjQaWbFiBXfeeSeXXHIJERERoS5dCHEOnPFmm+Xl5bz11lv88Ic/JCsri8bGRlQqFbNmzWLJkiUsXbqUJUuWTNrxN9JFJcTU5PF4KC0tpa6uLrj7d2pqKp2dnTQ2NnLw4EG6u7sBiI2N5corr+TLX/4yCxcuRK0+o3kVQohJZKy/32e9m/isWbM4fPgwERERlJSUUFRUFHyUlZXhcDjO5vQTRgKOEFNPZ2dnsNWmtbWVxMRE9Ho9zc3NlJSUUFxcHBxInJ2dzc0338x111035n3zhBCT3zkLOJ9HURRUKtVEnf6sSMARYurwer2Ul5dTU1NDfX09Pp+P9PR0bDYb9fX1HD58mLa2tuBA4g0bNnDrrbdy4YUXYjAYQl2+EGIcjfX3e0Lba88k3Dz55JOsWrUKk8mE1WrluuuuO2HfK0VRePzxx0lOTsZoNLJp0ybKy8vHq2whxCTS09PD7t27yc3NpaysjMjISNLT02lsbOTw4cO8//77tLS0oFKpmDFjBnfeeSc/+MEPuOyyyyTcCHEem3Qd0nv37uXrX/86hw8fZteuXfh8PrZs2YLT6Qy+5+mnn+aZZ57hueeeIzc3l8TERDZv3jxpu8OEEKfP5/NRUlLCnj17KCwspKuri+zsbDQaDUVFRXz44YccPHgQp9NJREQE69at49vf/jYPPvig7CUlhJjYLqrx0NPTg9VqZe/evVx00UUoikJycjIPPfQQ3/ve9wBwu90kJCTw1FNPcd99943pvNJFJcTk1dvbS1FREQ0NDbS2thIfH09MTAxNTU1UVlaSn5+P3W5HpVKRkJDAJZdcwq233sqSJUtkILEQ09yEroNzLtntdgBiYmIAaGhooLOzky1btgTfYzAY2LhxI4cOHTplwHG73bjd7uDzwcHBCaxaCHEmfD4flZWVHD16NLhBZnZ2NiMjIxQVFZGXl0dDQwNerxeDwcDcuXO57rrrZCCxEOIEkzrgKIrCww8/HFxjB47NogBISEgY9d6EhASamppOea4nn3ySJ554YuKKFUKclZ6enuCifa2trcTFxZGenk5LSwulpaUUFRXR398PQHx8PGvWrOHmm29mw4YN6PX6EFcvhJhsJnXAefDBBykpKeHAgQMnvPbZAcxfNGPr0Ucf5eGHHw4+HxwcJC0tbfyKFUKcEZ/PR0VFRbDVxuv1Mnv2bDweD0VFReTn51NfX4/b7SYsLIzMzEwuv/xybrzxRmbMmBHq8oUQk9SkDTjf+MY3ePvtt9m3b9+oZdUTExOBYy05SUlJwePd3d0ntOp8msFgkBkVQkwyPT09FBYW0tDQQFtbG1arlfj4eNrb24PratlsNhRFIS4ujqVLl3L99ddz2WWXER4eHuryhRCT2KQLOIqi8I1vfIM33niDPXv2kJmZOer1zMxMEhMT2bVrF8uWLQOOrWy6d+9ennrqqVCULIQ4TV6vN9hqU19fj9/vZ86cOfj9foqKiigsLKSurg63241OpyMjI4ONGzdyyy23MHv27Em7vpYQYvKYdAHn61//Oq+88gpvvfUWJpMpOObGYrFgNBpRqVQ89NBD7Ny5k+zsbLKzs9m5cyfh4eHcdtttIa5eiPObP6CQ02Cj2+HCagpjdWYMGvXoMNLZ2RmcIdXR0UFCQgJWq5X29nYKCgooLS0NttpERUWxYMECtm3bxtVXXy0zHoUQYzbpAs7zzz8PwKZNm0Yd/93vfsddd90FwCOPPMLIyAgPPPAA/f39rFmzhg8++ACTyXSOqxVCHPd+WQdPvFNBh90VPBYToeMn1y7kysXJuN1uysrKqKqqoqGhAYC5c+cSCATIzc2lsLCQ5ubm4FYLqamprFq1iltvvZWlS5fK9G8hxGmZ9OvgTBRZB0eI8fN+WQf3v1zAqf6Y3DLfyAUxw9TX19Pb20tSUhJWq5XW1lby8vIoLy9nYGCAQCCAxWJh9uzZbN26lWuvvZbY2Nhzei9CiMlt2qyDI4SY3PwBhSfeqThpuAl43Xh7GvnvslZa0zUsSIth3rx5KIrCkSNHKCgooLW1Ndhqk56ezvLly9m+fTsrV65Eq5U/UUKIMyN/PYQQZyWnwTaqWwqOTRbw2bvwdjfitXcScA2xT53AdZfOprWlhSNHjlBVVYXdbsfv9xMVFcWsWbO4+JJLmbHyYrr0JnKb7CcdwyOEEGMhAUcIcUaODyh+r6xj1PGAexh3dz2+/nb8gz2otAZ0cRkMB/z8398+pLvxKJ2dnaNabRYtWkTaiov5Q0s4Xa/XBc+VZAljx7b5XLEw6bOXF0KIzyVjcGQMjhCn7WQDihUlgM/WhqenCd9AJ4rXjcYchyrMhM/egaetikTVAGqfm0AgQFRUFJmZmWzcuJHYRRfxw3+0ndDNdbzt5vnbl0vIEUIAMgZHCDFBTjag2D8yiLur7lirzVA/akMEuvgMAgE/rsYivL2N+IftuMI1WCLDSUlJYe7cudx4442sWbuOS35x4KRjeBSOhZwn3qlg8/xE6a4SQoyZBBwhxJh9dkCx4vfh7Ws+1mpj7wJFQRuVhDosAm9vC+62SnyD3SgeNxoNJCckkpmZyYUXXsj1N9xIi8vAf+1tOGEMz6cpQIfdRU6DjXUzZUaVEGJsJOAIIcbs0wOKfY5ePF31+AY68A8Pogm3oDHHE3CPMFKTg7evCb/LieL3otYamJ+dwbqVS7jpppsYNM3gpj9Ufm6w+axux9jfK4QQEnCEEGPW7XAR8LjwdDfgtbXgH+xFpdagjUlFrTfi6arF3XYUv3MAxecGVBhN0SydN5NbrrmCm266icKeAN/4nDVzTsVqCpuIWxJCTFMScIQQYxIIBHD3tTFcn4d/oJOA14U2MgZ1ZCz+4X6GK3PxDnQQ8AwTE6YhPDqKBGsCq5ctZPv27axduxa1RssTv/34tMKNCki0HNv2QQghxkoCjhDiC/X19VFYWEjX0aOEDbbiVOnQxWWg0uhwtZTi6agl4BpC8bnRarTMzkwjNTWVjRs3cuONNxIfHw/AJ3V9p9UtdXxI8Y5t82WAsRDitEjAEUKcksfjobKykuLiYpqamlAUhe2XrODlon68tnZcLaX4Hb0oXjeK4ketD2fRnAzWrl7G9u3bWb58+ajViE93HE2irIMjhDhDEnCEECdQFIWWlhYKCgqoq6tjYGCA+Ph4MjIycLlcVBTkcLC2GN/IMAGfB7VGi9Ecy8oFWdx54zVce+21REdHn3DesY6jefDiWWyYFScrGQshzpgEHCHEKIODgxQWFlJZWUl7eztGo5F58+ZhsVjIzc3lwIED9Pf3kxKh4NLp0UdEEx8bx8b1q7n1lptZsGDBKXf+Xp0ZQ5IljE6766TjcI6Pt/n25tkSbIQQZ0UCjhACAJ/PR1VVFQUFBTQ1NeHz+UhLSyM5OZmOjg7efPNNGhsbGRkZwePxEB4eTlpqHCkpKWzdupUrr7wSk8n0udfQqFXs2Daf+18uQAWjQo6MtxFCjCcJOEII2tvbyc/Pp7q6moGBAWJiYsjIyECtVvPBBx+Qn5+P0+nE5XKh1WqxWq3ExcWxcuVKtm/fzqxZs1CpxhZKrliYxPO3Lz9hqwcZbyOEGE8ScIQ4jw0NDVFSUkJJSQnt7e3o9Xpmz55NXFwcxcXF7Nmzh76+PkZGRggEApjNZmJiYsjMzOSaa65h06ZNhIWd/vo0VyxMYvP8RHIabHQ7XFhNYTLeRggxriTgCHEe8vl8VFdXk5+fT1NTEx6Ph+TkZFJTU+nr6+MPf/gDDQ0NuN1uPB4PBoOBqKgo4uPj2bBhAzfddBNJSWfX0qJRq2TrBSHEhJGAI8R5pq2tjfz8fKqqqhgcHCQqKoq5c+ei1WrZtWsXhYWFo7qjoqKiiImJYc6cOSed+i2EEJOR/JUS4jzhcDgoLCyktLSUrq4u9Ho92dnZWK1WioqKgt1RHo+HQCBAREQE0dHRJCQkcNlll3HNNdcQFRUV6tsQQogxkYAjxDTn9XqpqqoiNzeXlpYWfD4fycnJpKen097ezksvvURjYyMulwufz4fBYMBkMhETE8PSpUu5+eabyc7OHvMgYiGEmAwk4AgxTR1frC8nJ4fa2lqcTifR0dFkZGSgKAp///vfKS0txel04vV60Wq1WCwWLBYLmZmZXHvttVx44YVnNIhYCCFCTQKOENOQzWYjPz+f8vJy+vr6CAsLY86cOURFRXHkyBFycnLo7+/H6/WiUqkwGo1YLBbi4+O58MILueGGG4L7RwkhxFQkAUeIacTlclFaWkpeXh6dnZ2oVCrS0o5tfFldXc0bb7xBe3s7brcbRVHQ6XREREQQFRXFggUL2L59OwsXLjzlSsRCCDFVSMARYgL5A8o5WeslEAhQX1/PoUOHaGxsxOv1BveO6u/v5y9/+Qt1dXU4nU4CgQAajQaDwYDFYiEtLY2rr76aSy65hPDw8HGvTQghQkECjhAT5P2yjhNW602agNV6Ozo6OHLkCFVVVQwNDREZGcmcOXMA2L17N+Xl5djtdnw+H2q1OjiI2Gq1sn79eq6//noSExPHrR4hhJgMJOAIMQHeL+vg/pcLTthQstPu4v6XC3j+9uVnHXLsdjv5+fkUFxdjs9nQ6XRkZWURHR1NXl4excXF9PT04PV60Wg06PV6jEYj0dHRLFy4kJtuuulzN8YUQoipTAKOEOPMH1B44p2Kk+6WrXBsU8kn3qlg8/zEM+qucrlclJWVceTIEbq6ugBISkoiNTWVyspK3nvvPTo7OxkeHkar1aLT6dDpdERFRZGWlsa2bdvYtGkTRqPxrO5TCCEmMwk4QoyznAbbqG6pz1KADruLnAbbaW1V4Pf7qa2t5dChQzQ3N+Pz+YiNjWXGjBm0tbXx17/+lba2NgYHB1GpVOj1egDMZjNWq5UNGzbwpS99CavVera3KIQQk54EHCHGWbfj1OHmTN4Hx7ZXOHToENXV1bhcLiIjI8nMzMThcPDee+/R2tpKX18fgUAAvV6Pz+cjLCyM2NhYlixZwg033MCcOXOkO0oIcd6QgCPEOLOaxrYw3ljeZ7PZOHz4MGVlZTgcDnQ6HbNmzUKj0XDgwAGampro7e3F6/ViMBjwer2o1WoSExPJysri2muvZe3atRgMhrO9LSGEmFIk4AgxzlZnxpBkCaPT7jrpOBwVkGg5NmX8VJxOJwUFBeTl5WGz2VCr1aSkpBAVFUVeXh6NjY309fUxNDREeHg4iqLg9/uJjY0lOTmZSy+9lCuvvBKz2Txh9ymEEJOZBBwhxplGrWLHtvnc/3IBKhgVco4PKd6xbf5JBxh7PB7Ky8s5dOhQcABxXFwcycnJlJSUsGfPHvr7++nv7yc8PByj0Yjb7cZsNhMfH8/q1au54YYbSE1NnfD7FEKIyUwCjhAT4IqFSTx/+/IT1sFJPMU6OIFAgNraWvbv309LSwt+vx+LxUJKSgq1tbXk5OQwODhIT09PcB0bh8NBREQE6enpzJs3jxtuuIEFCxag0WjO9e0KIcSko1IU5WSt6NPe4OAgFosFu90uzfhiwnzRSsaKotDc3Mz+/fupr6/H6/USHh5OSkoK7e3tlJeXMzQ0RG9vLwDh4eH09/ej1+uxWq1kZmZy9dVXs379ehlnI4Q4L4z191tacISYQBq16pRTwbu6uti/fz9Hjx7F7Xaj1+uZMWMG/f39fPjhhwwPD9PX14fX68VisdDf34/dbichIYGUlBQ2b97M5s2bJaALIcRJSMAR4hwbGBjg4MGDlJSUMDw8jEajISUlhZGREQ4cOIDT6WRgYICRkREsFguDg4N0d3cTGxtLYmIi69at49prr5XtFYQQ4nNIwBHiHHE6nRw5coTc3FyGhoZQqVQkJCTg8/mCx4aGhnA4HJjNZvx+Px0dHURHR5OamsqSJUu47rrrmD17NirV+G/YKYQQ04kEHCEmmMvlIj8/n8OHDzMwMABAVFQUiqJQVFSEw+HA5XIxODhIREQERqORjo4OIiIimDlzJvPmzeOaa65hyZIlaLXyj6wQQoyF/LUUYoJ4vV6Ki4s5ePBgcJBwZGQkKpWK6upq+vv76bIP0z9gJzIinOiICDo6OtDr9WRkZDBr1iyuuuoq1qxZQ1jY2BYPFEIIcYwEHCHGmd/vp6ysjP379wfXsgkLC0Oj0dDY2IjNZqOt30lZQweugAZFpUZprkev0bBgVhorF81ly5YtXHzxxURGRob4boQQYmqSgCPEOAkEAlRXV7N79246OjpQFAWdTodKpaKlpYX+/n4URaGho5fyzmECKj2BERsE/KjDLagsVupNC3h4+11su3BhqG9HCCGmNAk4QpwlRVGoqalhz549tLa2oigKKpUKlUpFa2srdrsdRVFwuVxoNFrq+v34h+0EvG7U4Ra0ZithKXMwZq1EZ4rl2QNd3LhhwUlXOhZCCDE2k3Jr4X379rFt2zaSk5NRqVS8+eabo15XFIXHH3+c5ORkjEYjmzZtory8PDTFivOWoijU1dXxP//zP7z88su0tLTg8/kIBAJ0dnZSVVXFwMAAXq8XlUpFeHg4Da3tDHW3gFaHLiGL8FlriFq/ncjFW9CaYlGADruLnAZbqG9PCCGmtEnZguN0OlmyZAl33303N9xwwwmvP/300zzzzDP8/ve/Z/bs2fzkJz9h8+bNVFVVYTKZQlCxON80NDSwZ88eGhoaUBQFj8dDIBBgYGCAwcFBFEVBURQMBgMajYa2tjbsdjteTRja+HQMCTMJn7kGbUwSKtWJ/57R7XCd5KpCCCHGalIGnK1bt7J169aTvqYoCs8++yyPPfYY119/PQAvvfQSCQkJvPLKK9x3333nslRxnmlqamL37t3U19cTCAQYGRnB6/UyPDyMw+EAQK1WExYWhsfjobW1lf7+fkwmE/PmzSM6NZt+WwL6uHRU6lPvGdXYO3yubkkIIaalSRlwPk9DQwOdnZ1s2bIleMxgMLBx40YOHTp0yoDjdrtxu93B54ODgxNeq5g+Ghsb2b17Nw0NDQQCARwOBx6PB7fbjdPpBECv1xMWFsbIiIuSo7X09fURGRHOvDlzWLpkCVu3bmXBwkUU/sdeOgfdn3u913KbefCSWTIORwghztCUCzidnZ0AJCQkjDqekJBAU1PTKT/35JNP8sQTT0xobWL6+XSw8fl82O12PB5PMNyoVCqMRiPh4eGMjIyQU1JBTVMnfpUGjSmOgcgsVHHrufPKq1m+LB2AW1en84sPaz73usfH4ZxqHyshhBCfb8oFnOM+u1T98Zkrp/Loo4/y8MMPB58PDg6SlpY2YfWJqa2hoYGPP/6YpqYm3G43/f39eL1e/H4/Ho8HjUZDeHg4kZGROJ1Ojh49SmNbF212L2pTPGHxmYRlrSAsaRZDGh3f+HMpWp2OKxYmMSMuYkw1yDgcIYQ4c1Mu4BzfYLCzs5OkpKTg8e7u7hNadT7NYDBgMBgmvD4xtdXW1rJ3716amppwOp309fXh9/sB8Hg8GAwGzGYzJpMJu91OaWkp/f39aLVa+tUWwjLTCctaTljSbFRaHQAKoAKeeKeCzfMTsZrGtirxWN8nhBDiRFMu4GRmZpKYmMiuXbtYtmwZcOyHZ+/evTz11FMhrk5MRYqicPToUfbt20dLSwuDg4P09fUFWwV9Ph/h4eGYzWYiIyMZGBggPz+fgYEBtFotaWlpxKXPodWeiCFxVjDYjLoG/+x2Wp0ZQ5IljE67C+Uk9aiAREsYqzNjJvrWhRBi2pqUAWdoaIja2trg84aGBoqKioiJiSE9PZ2HHnqInTt3kp2dTXZ2Njt37iQ8PJzbbrsthFWLqUZRFMrKyoLBxm6309fXh0ajQa1WoygKERERREREEB4eTk9PD4cPH8Zut6PVaklNTWXZsmVceeWVtGkS2P165Rdes9vhQqNWsWPbfO5/uQAVjAo5xztZd2ybLwOMhRDiLEzKgJOXl8fFF18cfH587Mydd97J73//ex555BFGRkZ44IEH6O/vZ82aNXzwwQeyBo4Yk0AgQFFREfv27aO1tRWbzUZ/fz86nQ6d7ljri8lkwmQyodfr6e7upqCgALvdjk6nCwabq6++mmXLlqHX6/mkrm9M1z7e7XTFwiSev305T7xTQYf9n2NtEi1h7Ng2nysWJp3qFEIIIcZApSjKyVrJp73BwUEsFgt2ux2z2RzqcsQ54PP5OHLkCAcPHgwGG4fDgU6nw2AwoNPpMJlMDHjUuAIKI3Yb/W2NDA050Ov1JCYmsnLlSq666ioWL16MXq8PntsfULjgqY+/sNvpwPcuGdUy4w8o5DTY6Ha4sJqOdUtJy40QQpzaWH+/J2ULjhDjaWRkhEOHDrF//346Ozux2WyMjIyg0+mIjIwkIiKCyMhI2hw+PiizYe/txGdrI+BxYTAYWDArnWs2b+Tqq69mwYIFaLUn/mNzpt1OGrVKpoILIcQEkIAjpq3BwUH27t3Lvn376O3tZWBgIDgTKiYmhsjISEwmE1qtltqeId45WIKvvwPF70Gl1aOLSkCXPJeWmatZevU1LFmS+rnXk24nIYSYPKSLSrqopp3u7m4+/PBDDh06FOyGOj4T6niwiYg4thbNsZ3Aa/n7wSI8Xi9qrR5NZDT6lHlEZK9FG52EWqU+affSqUi3kxBCTBzpohLnleM7e//973+npKQEu93O8PAwgUCAqKgoYmNjiYyMJCwsjEAggNfrpbq6mqamJhwjHnwqLbroOMIyFmOcuRqtOT64cOSnp3iPpTtJup2EECL0JOCIKc3v95Obm8vf//536urqGB4eZnh4GI1GQ0JCAtHR0URGRqLRaPB6vQwNDVFdXU17eztwbLZURHwUvaoZGGetRBcRfcprycrCQggxdUjAEVOS0+nkH//4B7t376arq4vh4WFGRkaIjIxk5syZmM1mjEYjiqLg8Xjo7++ntraWrq4utFot0dHRpKSksHXrVmJmr+be//3iNWxkZWEhhJg6JOCIkDudMSt1dXW8++675OXlMTQ0hNPpxOPxYLVamTNnDuHh4eh0OrxeL8PDw/T19dHQ0EB/fz8GgwGr1UpmZiZXXXUVl156KeHh4fgDCkn/aJCVhYUQYhqRgCNC6v2yjhNmHSV9ZtbR8PAwR44c4aOPPqK2than08nw8DBarZaUlBQSEhIICwtDo9EwNDTE0NBQMNg4nU6MRiOpqanMnj2bbdu2sXbdegpbHeyq7sdqGmF1ZoysLCyEENOMzKKSWVQh835ZB/e/XHBCq4kKUAJ+Hr8kAU1PNYcPH6arq4uhoSFcLhcxMTGkpaVhNpuDg4b7+/txuVz09PTQ1taG2+3GbDYTHR3NokWL2LZtG8uXL+eDiq5TBirgC8OWEEKI0Brr77cEHAk4YzLeU5+Pr/z76TChKAECw3a8tjbcrZVohzpZm6wLDhpOTEwkJSWFsLAwDAYDLpcLm83G0NAQPT099PT0oCgKFouF+Pj44HYKc+fOBT4/UAE8f/tyNs9PlCneQggxick0cTFuxtKNdLpyGmzB8/ldQ/gGu/F0N+LtbcY/0Inf7QS/j/7oZBbNnh2c5q0oCn19fXR1dWG32+np6cHhcKDRaIiNjSUlJYU1a9awZcsW0tLSgtfzBxSeeKfipGNsFI6FnCfeqWDz/ESZ4i2EENOABJzzxJm2wJyq1aPT7uL+lwt4/vblZxRyGjt68Pa14BvowtPfhrevjYCjl4DPDWo1aqMZbWQ0abOyyMpKxOVy0dzczMjICA6Hg76+PtxuNwaDgaSkJGbMmMGGDRu46KKLiI+PP+F6nw5UJ3O6a90IIYSY3CTgnAfOtAXmdFo9xhKWHA4H7e3ttLW1cbSwiqGjB48Fm+EBUKlQ68LQWBLRGE1owiJRUKFRji3INzQ0xODgYHDxPqPRSEpKCjNnzuSCCy5g9erVWCyWU157rGvYyFo3QggxPUjAmebOpgVmPFo9nE5nMNTYbDa6u7uprKykobEJVWcbSkCFOjwajTEctdGCWhdGwOfBZ+9Gr3gY7jXT7XDg9XrRarWYzWZSUlKYM2cOa9asYcGCBZhMpi/8Hsa6ho2sdSOEENODBJxp7GxbYM601WNwcJCOjg46OjpwOBw4HA6ampqoqKhgYGAAgPDwcNYsmsXhNg8qoxmVSo3fOYBnoJOAzwUeN7MSIhhyONDpdMTFxTFjxgzmzJnDsmXLyM7Oxmg0jvm7WJ0ZQ5IlTNa6EUKI84QEnGnsbFtgxtqaER9poL+/PxhqhoeH8Xq99PT0UF5eTldXFyMjIxgMBqKjo4mIiCAlJYXw8HASGrp4/3A5Qw47it+N4vWi16mZlxFLSnQkcXFxZGdnk52dzbx588jMzESv15/2d6FRq2StGyGEOI9IwJnGznbcyee1eihKAGV4kCjVMIN1BRyodKMoCna7nYaGBpqbmxkYGEClUqHRaIiLiyM+Pp7Y2FjUajWtra309PTgdDpZEu3BrteA1kJkRDgJURGkp6Uxa9Ys0tPTyc7OJiMjA41Gc1bfxxULk3j+9uUnjEdKlLVuhBBi2pGAM42d7biTz7Z6BPxe/M4B/M5+AsN2UAJcv2kmziEHHR0dNDU10dXVhcvlQqVSYTQaMZvNJCcnYzAYsNlsHD16FJvNhtvtxuv1otPpiImOJsGqw2QyMXPmTNLT00lISGDWrFkkJSWhVqvH7Tu5YmGSrHUjhBDnAQk409h4jDtZnx7Bjy6K5pm3c+nu7QsejzZquXpuFEZnJx991MTg4CBerxe1Wo3FYiExMRGTyYTL5RrVWuPxeNBoNERFRREZGYlWqyUuLo6ZM2cSFxdHYmJi8L9PFI1aJVPBhRBimpOAM42dybgTn89HT08P3d3ddHd343K5iAd+sjWTmq44eh3DKCODqIZ6sTW20OZwoCgKRqORpKQk4uPjg1sm1NbWMjg4yMjICAAmk4mUlBQMBkNwf6jU1NRRU77HMiNKCCGE+CIScKa5Lxp3cvmCRAYHB4OBxmaz8endOzQaDXq9nu7ubgL9rThbWrDb7Xg8HrRaLUlJScTExKAoCjabjdLSUvr7+xkZGcHn8xEeHk5qairR0dHBLqu0tDQsFgt6vZ4ZM2YwY8YMwsJkerYQQojxI3tRnSd7UX16JWOLHmZE+LH19dLT04PH4xn13oiICMLDw+np6aGlpYX6+noGBgYYGhoCID4+nujoaAwGAz09PQwMDNDf34/D4cDj8RAWFkZ8fDzx8fFYLBYiIyOJjo4mPj4erVZLZGQkWVlZpKamnvXAYSGEEOcX2YtKBHm9Xvr6+jCN9OLq78ExNETpp14/Psvp+EDguro6Ghsb6e/vx2634/P5iI6OJjs7G61WS19fH42NjTgcjmAXlE6nIz4+HqvVGgw2ZrMZi8USbJ2Jj48nKyuL+Ph4VCoZ1CuEEGLiSMCZhrxeLzabjd7eXvr6+rDb7aNeV6lUREVFERMTg1arxWazUV5eTktLCzabDbvdzsjICBaLhZkzZwJgs9mor69neHgYh8PB0NAQWq2WmJgYZs2ahdVqxWq1EhMTg8FgwGAwAMfCU0pKCllZWTK+RgghxDkjAWca8Hg82Gw2+vr6ThpoACIjI4mNjSUqKopAIEB7ezuHDx+mq6sLm83G4OAgQ0NDREREkJaWRiAQwGazUVVVhcvlYnh4GLvdjkqlIiYmhsWLFxMXF0dSUhKJiYlERkbi9XqDLTNGo5HMzEzS09PR6XTn+isRQghxnpOAMwU5nU5sNlvwYR90UNPlYGDES5RRR3aCCYvZRGxsLLGxsYSFhTE4OEhLSwt5eXn09vbS398f3LxSp9MFu5a6urooLy/H4/Hgdrux2+0oioLZbGbBggUkJCSQmJhIVlYWMTExjIyMBAcUq1Qq4uLiyMzMJCEhQbqhhBBChIwEnEnO5/Nht9vp7+8PPtxud/D1/CYbr+W0MODToA4zoTaaSI7X8PAl8Vg1KsrLy4MbXdpsNpxOZ3AhvtjYWMxmM11dXZSVleFyufD7/djtdgKBAFFRUcybN4+kpCSsVmtwZWGfzxeccQXHuqFSU1PJzMyUbighhBCTgsyimkSzqBRFYWhoKBhkBgYGcPy/dWY+Ta1WExUVRXmfn5/sakYVZkIJ+AkM2/EPDxBw9uMfGeRL8ywkhSsMDw/j8/kIBAKEh4fj9Xrp6uqiv78fj8eDz+fD6XTi9/uxWCykpaWRmpoaHBQ8d+5cNBoNnZ2dwVADx7q9ZsyYQWpqqnRDCSGEOCdkFtUkdzzMDAwMYLfbgw+/33/Ce41GI1FRUURHRxMdHU1UVBRen59vPv4WvmE7gd4W/G4nAfcQgZEhAm4nBAK8m9vFzStTUKsgEAgEQ43X68Xn8zE8PAyAxWJh7ty5pKWlBbuf5s2bh8Viobu7m4aGhuBUcpVKRWJiIjNmzJjQ1YaFEEKIsyEB5xzw+XzBKdWDg4PY7XYGBwdPGmaOb2Pw6TATFhYW3Miyp6eHo0ePcqiiieajZQRcTgIjDgKeEVRqDSgBFJ+PgM/FgNNNXtkgWu9IsKXG4/EEW4A+ve/TzJkzmTt3LgkJCfT399PU1ER5eXmwrrCwMDIyMkhPT5dF+YQQQkx6EnDG2cjISDDAHH84nc6Tvler1WKxWLBYLERFRWGxWIiIiAgOznU4jm1ieXy6t9frxev10t/fT1l5A97eZlRaA6Cg+L34Hb34vW4UjxPF40LxeWnz6jGHadHr9VgsFlJTU5kxYwZJSUlkZ2eTnZ2N1WrF4/HQ3NzMvn37Ri38Z7VaycjIkEHDQgghphQJOOOsoqKC9vb2E46HhYVhMplGLYD36TADx2ZHtbS00NvbS29vb3AwscvlCs56UqvVhIeHYwkPwz88eGzMjdcDXheK143ic6MoflRqDSqNgYTERBbNziI9PZ3U1FTmzZtHRkZGsHups7MzOLPq07WmpaWRkZGB0Wic4G9MCCGEGH8ScMZZdHQ0TqcTs9k86qHX60e97/gYnONr19hsNlwuV/A1p9PJ4OAgPp8PrVaLwWDA6/VSW1tLX18fIy4Xqv4ufCPDKJ5jn1NpDagMRrRhJjSmWKLjk/i3e69k7pzZJCcnExUVhUqlwm63U1FRQVtbG16vN1iTtNYIIYSYLiTgjLOsrCyysrJOOB4IBILTvY8Hmk93Bfn9foaGhvD7/fj9fnQ6HYFAgIaGBpqamo7t9+T0MOJy43U5UQW8xGpUtKNCbYpFHWZCGxmNxpKAPi4DfXQSv/yXjVy78lgtXq+XpqYmmpubRy0EaDQaSU9PJy0tTVprhBBCTBsScCbI8dWF+/v7sdlsDAwMEAgERr3H5/OhKAp+vx+Px4PH46GtrY36+no6OjqCXVTtfQM0d/bj8/lAp0etjyDcFMuqeWmsjU8gpz8MZ0QKOlMMaqOZpKhwbl2djqLR8/eco8Srhujp7gpeX61Wk5iYSHp6OnFxcdJaI4QQYtqRgDPOampqaGlpOenAYp1Oh0ajwe/3MzIyQk9PDz09PbS3t9Pa2hrskvL7/TidTtxuNyNehY4RUEfGoguLRBMWgcYUhz5hJhWxaTx3z0ZeXpEV3Cm8sXeYl/cf5anXPsLv6EPxe4gJ13PL6jQuXjSD9PR0UlJSTugyE0IIIaYTCTjjzOPxBMONyWQiPDw8uDpwS0sL7e3t9PT0BNekcbvdeL3e4H8qioJerycmJoaYmFj2NzoJi4tAF5uCPn4GWksCaqMJlUqNCnjqw0a2rchiRZqJV3fX8e9/2I3f/c9wpVJrcehi+F2dkbWbZpOZmRSib0YIIYQ4dyTgjLPU1FR8Ph89PT1UVlbS1NQU3ADz+N5Pbrcbj8eDoijBVYojIiJITU0lOjo6uJqwJjaDwzlDGM1xqDQn/k8VCPhpaW3j9299SJzGxZN/Kcbv9gAqNBFRaE3xqCMswTD0xDsVbJ6fiEYtXVJCCCGmNwk44+x//ud/KC4uxuFw4HK5GBoawuVy4Xa70Wg0qNVqAAwGAxaLhZiYGKKiorBarcyePZvs7Gzmzp2L2WzmraI2dDVFo86vKAoBl+PYmjcOG4rip745gM2oY8CvQxeXhNYUi0ozeusEBeiwu8hpsLFuZuw5+jaEEEKI0JCAM866u7spKytDURRUKhU6nQ6t9thCe3q9nqioKGJiYkhKSiItLY2srCzmzJlDSkpKMPwcZzX9c8XggHsYv6MP31Afiu+fm22qtAYWzpuLKiKasPqoL67P4Rq3exVCCCEmKwk448zn82E2m1Gr1Wg0GoxGIzExMaSmppKYmEhqaiqzZs0iLS3tCzf5XGA1EOXtpbO9nYB3JHhcpdaiiYxGa4ojOSGe7ZetIafBBtR+YX2fDk1CCCHEdCUBZ5zNmDGDnp4eUlJSSExMJCkpiczMTJKSkrBarV+46/bIyAjt7e20tbVht9u5LhOebxpBpVKjDo9CY4pFExGFWnWstefxaxagUatYnRlDkiWMTruLk20PrwISLWGszowZ/5sWQgghJpkpHXB+9atf8bOf/YyOjg4WLFjAs88+y4UXXhjSmq6++mqWLFlCQkICVqsVi8XyhevMuFwu2tvb6ejowGazBY+rVCq2rprLzLkLee5wL11DvuBriZYwdmybzxULj82K0qhV7Ng2n/tfLkAFo0LO8avv2DZfBhgLIYQ4L6iU49N4ppg///nP3HHHHfzqV79iw4YN/OY3v+G///u/qaioID09/Qs/Pzg4iMViwW63f2FX0UQ4VagBiImJISUlheTk5OB6Nf6AElzrxmo61hJzsrDyflkHT7xTQYf9n2Ntkj4ThoQQQoipaqy/31M24KxZs4bly5fz/PPPB4/NmzeP6667jieffPILPx+KgDMyMkJHR8cpQ01SUhLJycmEhZ3dOJmxhiEhhBBiqhnr7/eU7KLyeDzk5+fz/e9/f9TxLVu2cOjQoZN+xu12B7c+gGNf0LkwNDRER0cHnZ2dDAwMjHptPEPNp2nUKpkKLoQQ4rw2JQNOb28vfr+fhISEUccTEhLo7Ow86WeefPJJnnjiiXNRHna7PRhqHA7HqNdiY2NJTEwc91AjhBBCiH+akgHnuM8O3j2+9szJPProozz88MPB54ODg6SlpY17TXl5eXR0dASfq9Vq4uLiSEpKIiEhAYPBMO7XFEIIIcRoUzLgxMXFodFoTmit6e7uPqFV5ziDwXBOwkV0dDTd3d1YrdYxTw0XQgghxPiakgFHr9ezYsUKdu3axZe+9KXg8V27dnHttdeGsDLIyMhgxowZaDSakNYhhBBCnM+mZMABePjhh7njjjtYuXIl69at44UXXqC5uZmvfe1rIa1Lq52yX6kQQggxbUzZX+Obb76Zvr4+fvzjH9PR0cHChQt59913ycjICHVpQgghhAixKbsOztkK9UJ/QgghhDh9Y/39Vp/yFSGEEEKIKUoCjhBCCCGmHQk4QgghhJh2JOAIIYQQYtqRgCOEEEKIaUcCjhBCCCGmHQk4QgghhJh2JOAIIYQQYtqRgCOEEEKIaUcCjhBCCCGmHQk4QgghhJh2puxmm2fr+BZcg4ODIa5ECCGEEGN1/Hf7i7bSPG8DjsPhACAtLS3ElQghhBDidDkcDiwWyylfP293Ew8EArS3t2MymVCpVON23sHBQdLS0mhpaZFdysdAvq+xk+9q7OS7Gjv5rsZOvquxm8jvSlEUHA4HycnJqNWnHmlz3rbgqNVqUlNTJ+z8ZrNZ/gE4DfJ9jZ18V2Mn39XYyXc1dvJdjd1EfVef13JznAwyFkIIIcS0IwFHCCGEENOOBJxxZjAY2LFjBwaDIdSlTAnyfY2dfFdjJ9/V2Ml3NXbyXY3dZPiuzttBxkIIIYSYvqQFRwghhBDTjgQcIYQQQkw7EnCEEEIIMe1IwBFCCCHEtCMBZ5z96le/IjMzk7CwMFasWMH+/ftDXdKk8+STT7Jq1SpMJhNWq5XrrruOqqqqUJc1JTz55JOoVCoeeuihUJcyKbW1tXH77bcTGxtLeHg4S5cuJT8/P9RlTUo+n48f/vCHZGZmYjQaycrK4sc//jGBQCDUpYXcvn372LZtG8nJyahUKt58881RryuKwuOPP05ycjJGo5FNmzZRXl4emmJD7PO+K6/Xy/e+9z0WLVpEREQEycnJfOUrX6G9vf2c1CYBZxz9+c9/5qGHHuKxxx6jsLCQCy+8kK1bt9Lc3Bzq0iaVvXv38vWvf53Dhw+za9cufD4fW7Zswel0hrq0SS03N5cXXniBxYsXh7qUSam/v58NGzag0+l47733qKio4Oc//zlRUVGhLm1Seuqpp/j1r3/Nc889R2VlJU8//TQ/+9nP+OUvfxnq0kLO6XSyZMkSnnvuuZO+/vTTT/PMM8/w3HPPkZubS2JiIps3bw7ucXg++bzvanh4mIKCAn70ox9RUFDA66+/TnV1Nddcc825KU4R42b16tXK1772tVHH5s6dq3z/+98PUUVTQ3d3twIoe/fuDXUpk5bD4VCys7OVXbt2KRs3blS+9a1vhbqkSed73/uecsEFF4S6jCnjqquuUu65555Rx66//nrl9ttvD1FFkxOgvPHGG8HngUBASUxMVP793/89eMzlcikWi0X59a9/HYIKJ4/Pflcnk5OTowBKU1PThNcjLTjjxOPxkJ+fz5YtW0Yd37JlC4cOHQpRVVOD3W4HICYmJsSVTF5f//rXueqqq7jssstCXcqk9fbbb7Ny5UpuuukmrFYry5Yt48UXXwx1WZPWBRdcwEcffUR1dTUAxcXFHDhwgCuvvDLElU1uDQ0NdHZ2jvpbbzAY2Lhxo/ytHwO73Y5KpTonLavn7Wab4623txe/309CQsKo4wkJCXR2doaoqslPURQefvhhLrjgAhYuXBjqcial1157jfz8fPLy8kJdyqRWX1/P888/z8MPP8wPfvADcnJy+OY3v4nBYOArX/lKqMubdL73ve9ht9uZO3cuGo0Gv9/PT3/6U2699dZQlzapHf97frK/9U1NTaEoacpwuVx8//vf57bbbjsnm5VKwBlnKpVq1HNFUU44Jv7pwQcfpKSkhAMHDoS6lEmppaWFb33rW3zwwQeEhYWFupxJLRAIsHLlSnbu3AnAsmXLKC8v5/nnn5eAcxJ//vOfefnll3nllVdYsGABRUVFPPTQQyQnJ3PnnXeGurxJT/7Wnx6v18stt9xCIBDgV7/61Tm5pgSccRIXF4dGozmhtaa7u/uEpC+O+cY3vsHbb7/Nvn37SE1NDXU5k1J+fj7d3d2sWLEieMzv97Nv3z6ee+453G43Go0mhBVOHklJScyfP3/UsXnz5vHXv/41RBVNbt/97nf5/ve/zy233ALAokWLaGpq4sknn5SA8zkSExOBYy05SUlJwePyt/7UvF4v27dvp6GhgY8//victN6AzKIaN3q9nhUrVrBr165Rx3ft2sX69etDVNXkpCgKDz74IK+//joff/wxmZmZoS5p0rr00kspLS2lqKgo+Fi5ciVf/vKXKSoqknDzKRs2bDhhuYHq6moyMjJCVNHkNjw8jFo9+idAo9HINPEvkJmZSWJi4qi/9R6Ph71798rf+pM4Hm5qamr48MMPiY2NPWfXlhaccfTwww9zxx13sHLlStatW8cLL7xAc3MzX/va10Jd2qTy9a9/nVdeeYW33noLk8kUbPWyWCwYjcYQVze5mEymE8YmRUREEBsbK2OWPuPb3/4269evZ+fOnWzfvp2cnBxeeOEFXnjhhVCXNilt27aNn/70p6Snp7NgwQIKCwt55plnuOeee0JdWsgNDQ1RW1sbfN7Q0EBRURExMTGkp6fz0EMPsXPnTrKzs8nOzmbnzp2Eh4dz2223hbDq0Pi87yo5OZkbb7yRgoIC/va3v+H3+4N/72NiYtDr9RNb3ITP0zrP/Nd//ZeSkZGh6PV6Zfny5TL1+SSAkz5+97vfhbq0KUGmiZ/aO++8oyxcuFAxGAzK3LlzlRdeeCHUJU1ag4ODyre+9S0lPT1dCQsLU7KyspTHHntMcbvdoS4t5Hbv3n3Sv1F33nmnoijHporv2LFDSUxMVAwGg3LRRRcppaWloS06RD7vu2poaDjl3/vdu3dPeG0qRVGUiY1QQgghhBDnlozBEUIIIcS0IwFHCCGEENOOBBwhhBBCTDsScIQQQggx7UjAEUIIIcS0IwFHCCGEENOOBBwhhBBCTDsScIQQQggx7UjAEUIIIcS0IwFHCCGEENOOBBwhhBBCTDsScIQQ08bOnTtRqVQnPJ555plQlyaEOMdks00hxLThcDhwOp3B5z/+8Y959913OXDgAKmpqSGsTAhxrmlDXYAQQowXk8mEyWQC4IknnuDdd99l7969Em6EOA9JF5UQYtp54okn+N3vfsfevXvJyMgIdTlCiBCQgCOEmFYk3AghQAKOEGIakXAjhDhOxuCI/799O7ZiEIiBKLjZpVw1NEgJlEIRlEJMA7gA5+Z530yqROEPJKiwbVv2fc9xHBlj5LquJMmcM2OMl7cDfs0XFfD3nufJsiy57/trdp5n1nV9YSvgTQIHAKjjBgcAqCNwAIA6AgcAqCNwAIA6AgcAqCNwAIA6AgcAqCNwAIA6AgcAqCNwAIA6AgcAqPMBkdwE0ktiuOgAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Build model\n",
    "@model function linear_regression(N, Σ, σ2)\n",
    "\n",
    "    w ~ MvNormalMeanCovariance(constvar(zeros(3)),constvar(Σ))\n",
    "   \n",
    "    x = datavar(Vector{Float64}, N)\n",
    "    y = datavar(Float64, N)\n",
    "    \n",
    "    for i in 1:N\n",
    "        y[i] ~ NormalMeanVariance(dot(w , x[i]), σ2)\n",
    "    end\n",
    "    return w, x, y\n",
    "end\n",
    "# Run message passing algorithm \n",
    "results = inference(\n",
    "    model = linear_regression(length(x_train), Σ, σ2 ),\n",
    "    data  = (y = y_train, x = x_train),\n",
    "    returnvars   = (w = KeepLast()),\n",
    "    iterations   = 20\n",
    ");\n",
    "# Plot result\n",
    "w = results.posteriors[:w]\n",
    "println(\"Posterior distribution of w: $(w)\")\n",
    "scatter(z, y_train); xlabel(L\"z\"); ylabel(L\"f([1.0, z, z^2]) + ϵ\");\n",
    "z_test = collect(0:0.2:12)\n",
    "x_test = [[1.0; z; z^2] for z in z_test]\n",
    "for i=1:10\n",
    "    w_sample = rand(results.posteriors[:w])\n",
    "    f_est(x) = (w_sample'*x)[1]\n",
    "    plot(z_test, map(f_est, x_test), \"k-\", alpha=0.3);\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final thoughts: Modularity and Abstraction \n",
    "\n",
    "- The great Michael Jordan (no, not [this one](https://youtu.be/cuLprHh_BRg), but [this one](https://people.eecs.berkeley.edu/~jordan/)), wrote:   \n",
    "\n",
    "  > I basically know of two principles for treating complicated systems in simple ways: the first is the principle of **modularity** and the second is the principle of **abstraction**. I am an apologist for computational probability in machine learning because I believe that probability theory implements these two principles in deep and intriguing ways — namely through factorization and through averaging. Exploiting these two mechanisms as fully as possible seems to me to be the way forward in machine learning. — Michael Jordan, 1997 (quoted in [Fre98](https://mitpress.mit.edu/9780262062022/)).\n",
    "\n",
    "- Factor graphs realize these ideas nicely, both visually and computationally.\n",
    "\n",
    "- Visually, the modularity of conditional independencies in the model are displayed by the graph structure. Each node hides internal complexity and by closing-the-box, we can hierarchically move on to higher levels of abstraction. \n",
    "\n",
    "- Computationally, message passing-based inference uses the Distributive Law to avoid any unnecessary computations.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>OPTIONAL SLIDES</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###  <a id=\"sp-for-multiplication-node\">Sum-Product Messages for Multiplication Nodes</a>\n",
    "- Next, let us consider a **multiplication** by a fixed (invertible matrix) gain $f_A(x,y) = \\delta(y-Ax)$\n",
    "\n",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-gain-node.png\" width=\"250px\"></p>"
   ]
  },
  {
   "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 \\\\\n",
    "&= \\int  \\overrightarrow{\\mu}_{X}(x) \\,|A|^{-1}\\delta(x-A^{-1}y) \\,\\mathrm{d}x \\\\\n",
    "&= |A|^{-1}\\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) &= |A|^{-1}\\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": [
    "- <b>Exercise</b>: 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 RxInfer. \n",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-addition-node.png\" width=\"250px\"></p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward message on Z:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "NormalMeanVariance{Float64}(μ=3.0, v=2.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"Forward message on Z:\")\n",
    "@call_rule typeof(+)(:out, Marginalisation) (m_in1 = NormalMeanVariance(1.0, 1.0), m_in2 = NormalMeanVariance(2.0, 1.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Backward message on X:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "NormalMeanVariance{Float64}(μ=1.0, v=2.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"Backward message on X:\")\n",
    "@call_rule typeof(+)(:in1, Marginalisation) (m_out = NormalMeanVariance(3.0, 1.0), m_in2 = NormalMeanVariance(2.0, 1.0))"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-gain-node.png\" width=\"250px\"></p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward message on Y:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "NormalMeanVariance{Float64}(μ=4.0, v=16.0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"Forward message on Y:\")\n",
    "@call_rule typeof(*)(:out, Marginalisation) (m_A = PointMass(4.0), m_in = NormalMeanVariance(1.0, 1.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Backward message on X:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "NormalWeightedMeanPrecision{Float64}(xi=8.0, w=16.0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"Backward message on X:\")\n",
    "@call_rule typeof(*)(:in, Marginalisation) (m_out = NormalMeanVariance(2.0, 1.0), m_A = PointMass(4.0), meta = TinyCorrection())"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-observations.png\" width=\"375px\"></p>"
   ]
  },
  {
   "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",
    "<p style=\"text-align:center;\"><img src=\"./figures/ffg-observations-2.png\" width=\"650px\"></p>"
   ]
  },
  {
   "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 RxInfer 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": 12,
   "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": [
    "# Data\n",
    "y1_hat = 1.0\n",
    "y2_hat = 2.0\n",
    "\n",
    "# Construct the factor graph\n",
    "@model function my_model()\n",
    "\n",
    "    # `x` is the hidden states\n",
    "    x = randomvar()\n",
    "    # `y1` and `y2` are \"clamped\" observations\n",
    "    y1 = datavar(Float64,)\n",
    "    y2 = datavar(Float64,)\n",
    "    \n",
    "    x ~ NormalMeanVariance(constvar(0.0), constvar(4.0))\n",
    "    y1 ~ NormalMeanVariance(x, constvar(1))\n",
    "    y2 ~ NormalMeanVariance(x, constvar(2))\n",
    "    \n",
    "    return x\n",
    "end\n",
    "\n",
    "result = inference(model=my_model(), data=(y1=y1_hat, y2 = y2_hat,))\n",
    "println(\"Sum-product message passing result: p(x|y1,y2) = 𝒩($(mean(result.posteriors[:x])),$(var(result.posteriors[:x])))\")\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": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<!--\r\n",
       "This HTML file contains custom styles and some javascript.\r\n",
       "Include it a Jupyter notebook for improved rendering.\r\n",
       "-->\r\n",
       "\r\n",
       "<!-- Fonts -->\r\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\r\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\r\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\r\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\r\n",
       "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\r\n",
       "\r\n",
       "<!-- Custom style -->\r\n",
       "<style>\r\n",
       "\r\n",
       "@font-face {\r\n",
       "    font-family: \"Computer Modern\";\r\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\r\n",
       "}\r\n",
       "\r\n",
       "#notebook_panel { /* main background */\r\n",
       "    background: rgb(245,245,245);\r\n",
       "}\r\n",
       "\r\n",
       "div.container {\r\n",
       "    min-width: 960px;\r\n",
       "}\r\n",
       "\r\n",
       "div #notebook { /* centre the content */\r\n",
       "    background: #fff; /* white background for content */\r\n",
       "    margin: auto;\r\n",
       "    padding-left: 0em;\r\n",
       "}\r\n",
       "\r\n",
       "#notebook li { /* More space between bullet points */\r\n",
       "    margin-top:0.8em;\r\n",
       "}\r\n",
       "\r\n",
       "/* draw border around running cells */\r\n",
       "div.cell.border-box-sizing.code_cell.running {\r\n",
       "    border: 1px solid #111;\r\n",
       "}\r\n",
       "\r\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\r\n",
       "div.cell.code_cell {\r\n",
       "    background-color: rgb(256,256,256);\r\n",
       "    border-radius: 0px;\r\n",
       "    padding: 0.5em;\r\n",
       "    margin-left:1em;\r\n",
       "    margin-top: 1em;\r\n",
       "}\r\n",
       "\r\n",
       "div.text_cell_render{\r\n",
       "    font-family: 'Alegreya Sans' sans-serif;\r\n",
       "    line-height: 140%;\r\n",
       "    font-size: 125%;\r\n",
       "    font-weight: 400;\r\n",
       "    width:800px;\r\n",
       "    margin-left:auto;\r\n",
       "    margin-right:auto;\r\n",
       "}\r\n",
       "\r\n",
       "\r\n",
       "/* Formatting for header cells */\r\n",
       ".text_cell_render h1 {\r\n",
       "    font-family: 'Nixie One', serif;\r\n",
       "    font-style:regular;\r\n",
       "    font-weight: 400;\r\n",
       "    font-size: 45pt;\r\n",
       "    line-height: 100%;\r\n",
       "    color: rgb(0,51,102);\r\n",
       "    margin-bottom: 0.5em;\r\n",
       "    margin-top: 0.5em;\r\n",
       "    display: block;\r\n",
       "}\r\n",
       "\r\n",
       ".text_cell_render h2 {\r\n",
       "    font-family: 'Nixie One', serif;\r\n",
       "    font-weight: 400;\r\n",
       "    font-size: 30pt;\r\n",
       "    line-height: 100%;\r\n",
       "    color: rgb(0,51,102);\r\n",
       "    margin-bottom: 0.1em;\r\n",
       "    margin-top: 0.3em;\r\n",
       "    display: block;\r\n",
       "}\r\n",
       "\r\n",
       ".text_cell_render h3 {\r\n",
       "    font-family: 'Nixie One', serif;\r\n",
       "    margin-top:16px;\r\n",
       "    font-size: 22pt;\r\n",
       "    font-weight: 600;\r\n",
       "    margin-bottom: 3px;\r\n",
       "    font-style: regular;\r\n",
       "    color: rgb(102,102,0);\r\n",
       "}\r\n",
       "\r\n",
       ".text_cell_render h4 {    /*Use this for captions*/\r\n",
       "    font-family: 'Nixie One', serif;\r\n",
       "    font-size: 14pt;\r\n",
       "    text-align: center;\r\n",
       "    margin-top: 0em;\r\n",
       "    margin-bottom: 2em;\r\n",
       "    font-style: regular;\r\n",
       "}\r\n",
       "\r\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\r\n",
       "    font-family: 'Nixie One', sans-serif;\r\n",
       "    font-weight: 400;\r\n",
       "    font-size: 16pt;\r\n",
       "    color: rgb(163,0,0);\r\n",
       "    font-style: italic;\r\n",
       "    margin-bottom: .1em;\r\n",
       "    margin-top: 0.8em;\r\n",
       "    display: block;\r\n",
       "}\r\n",
       "\r\n",
       ".text_cell_render h6 { /*use this for copyright note*/\r\n",
       "    font-family: 'PT Mono', sans-serif;\r\n",
       "    font-weight: 300;\r\n",
       "    font-size: 9pt;\r\n",
       "    line-height: 100%;\r\n",
       "    color: grey;\r\n",
       "    margin-bottom: 1px;\r\n",
       "    margin-top: 1px;\r\n",
       "}\r\n",
       "\r\n",
       ".CodeMirror{\r\n",
       "    font-family: \"PT Mono\";\r\n",
       "    font-size: 90%;\r\n",
       "}\r\n",
       "\r\n",
       ".boxed { /* draw a border around a piece of text */\r\n",
       "  border: 1px solid blue ;\r\n",
       "}\r\n",
       "\r\n",
       "h4#CODE-EXAMPLE,\r\n",
       "h4#END-OF-CODE-EXAMPLE {\r\n",
       "    margin: 10px 0;\r\n",
       "    padding: 10px;\r\n",
       "    background-color: #d0f9ca !important;\r\n",
       "    border-top: #849f81 1px solid;\r\n",
       "    border-bottom: #849f81 1px solid;\r\n",
       "}\r\n",
       "\r\n",
       ".emphasis {\r\n",
       "    color: red;\r\n",
       "}\r\n",
       "\r\n",
       ".exercise {\r\n",
       "    color: green;\r\n",
       "}\r\n",
       "\r\n",
       ".proof {\r\n",
       "    color: blue;\r\n",
       "}\r\n",
       "\r\n",
       "code {\r\n",
       "  padding: 2px 4px !important;\r\n",
       "  font-size: 90% !important;\r\n",
       "  color: #222 !important;\r\n",
       "  background-color: #efefef !important;\r\n",
       "  border-radius: 2px !important;\r\n",
       "}\r\n",
       "\r\n",
       "/* This removes the actual style cells from the notebooks, but no in print mode\r\n",
       "   as they will be removed through some other method */\r\n",
       "@media not print {\r\n",
       "  .cell:nth-last-child(-n+2) {\r\n",
       "    display: none;\r\n",
       "  }\r\n",
       "}\r\n",
       "\r\n",
       "footer.hidden-print {\r\n",
       "    display: none !important;\r\n",
       "}\r\n",
       "    \r\n",
       "</style>\r\n",
       "\r\n",
       "<!-- MathJax styling -->\r\n",
       "<script>\r\n",
       "    MathJax.Hub.Config({\r\n",
       "                        CommonHTML: {\r\n",
       "                            scale: 200\r\n",
       "                        },\r\n",
       "                TeX: {\r\n",
       "                    extensions: [\"AMSmath.js\"],\r\n",
       "                    equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\r\n",
       "                },\r\n",
       "                tex2jax: {\r\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\r\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\r\n",
       "                },\r\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\r\n",
       "                \"HTML-CSS\": {\r\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\r\n",
       "                }\r\n",
       "        });\r\n",
       "</script>\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "open(\"../../styles/aipstyle.html\") do f display(\"text/html\", read(f, String)) end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Julia 1.8.2",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}