part2/labs/11636785_11640758_lab2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Authors:\n",
    "Gabriele Bani 11640758\n",
    "\n",
    "Andrii Skliar 11636785"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "96a94c5280b3796c2e3bda8d652f9ee1",
     "grade": false,
     "grade_id": "cell-66ada6339bda77a6",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Lab 2: Inference in Graphical Models\n",
    "\n",
    "### Machine Learning 2 (2017/2018)\n",
    "\n",
    "* The lab exercises should be made in groups of two people.\n",
    "* The deadline is Thursday, 29.04, 23:59.\n",
    "* Assignment should be submitted through BlackBoard! Make sure to include your and your teammates' names with the submission.\n",
    "* Attach the .IPYNB (IPython Notebook) file containing your code and answers. Naming of the file should be \"studentid1\\_studentid2\\_lab#\", for example, the attached file should be \"12345\\_12346\\_lab1.ipynb\". Only use underscores (\"\\_\") to connect ids, otherwise the files cannot be parsed.\n",
    "\n",
    "Notes on implementation:\n",
    "\n",
    "* You should write your code and answers in an IPython Notebook: http://ipython.org/notebook.html. If you have problems, please ask.\n",
    "* Use __one cell__ for code and markdown answers only!\n",
    "    * Put all code in the cell with the ```# YOUR CODE HERE``` comment and overwrite the ```raise NotImplementedError()``` line.\n",
    "    * For theoretical questions, put your solution using LaTeX style formatting in the YOUR ANSWER HERE cell.\n",
    "* Among the first lines of your notebook should be \"%pylab inline\". This imports all required modules, and your plots will appear inline.\n",
    "* Large parts of you notebook will be graded automatically. Therefore it is important that your notebook can be run completely without errors and within a reasonable time limit. To test your notebook before submission, select Kernel -> Restart \\& Run All."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "2a0482dac028fae34551f9e75cc82068",
     "grade": false,
     "grade_id": "cell-7c4914def40aeb03",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### Introduction\n",
    "In this assignment, we will implement the sum-product and max-sum algorithms for factor graphs over discrete variables. The relevant theory is covered in chapter 8 of Bishop's PRML book, in particular section 8.4. Read this chapter carefuly before continuing!\n",
    "\n",
    "We will implement sum-product and max-sum and apply it to a poly-tree structured medical diagnosis example.\n",
    "\n",
    "For this assignment you should use numpy ndarrays (constructed with np.array, np.zeros, np.ones, etc.). We need n-dimensional arrays in order to store conditional distributions with more than one conditioning variable. If you want to perform matrix multiplication on arrays, use the np.dot function; all infix operators including *, +, -, work element-wise on arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "b85b4dd3e9c765f6ed0d81b98161f0cd",
     "grade": false,
     "grade_id": "cell-7f9d339abedb92d3",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Part 0: Doing the math (5 points)\n",
    "We start with a set of three Bernoulli distributed variables X, Y, Z. Calculate the marginals $p(Y=1)$ and $p(Y=0)$ using only the sum and product rule where,\n",
    "\n",
    "$$\n",
    "p(X=1) = 0.05 \\\\\\\\\n",
    "p(Z=1) = 0.2 \\\\\\\\\n",
    "$$\n",
    "$$\n",
    "p(Y = 1 | X = 1, Z = 1) = 0.99 \\\\\\\\\n",
    "p(Y = 1 | X = 1, Z = 0) = 0.9 \\\\\\\\\n",
    "p(Y = 1 | X = 0, Z = 1) = 0.7 \\\\\\\\\n",
    "p(Y = 1 | X = 0, Z = 0) = 0.0001 \\\\\\\\\n",
    "$$\n",
    "\n",
    "While implementing the message passing algorithms you should be able to use the results of this question as a guidance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "cafc6734ba3ef70e630191e714f3c0b3",
     "grade": true,
     "grade_id": "cell-653065eb00c64d46",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "\\begin{align*}\n",
    "p(Y=1) &= p(Y=1|X=1,Z=1)\\cdot p(X=1)\\cdot p(Z=1) + p(Y=1|X=1,Z=0)\\cdot p(X=1)\\cdot p(Z=0)\\\\\n",
    "&+ p(Y=1|X=0,Z=1)\\cdotp(X=0)\\cdotp(Z=1) + p(Y=1|X=0,Z=0)\\cdot p(X=0)\\cdot p(Z=0) \\\\\n",
    "&= p(Y=1|X=1,Z=1)\\cdot p(X=1)\\cdot p(Z=1) + p(Y=1|X=1,Z=0)\\cdot p(X=1)\\cdot (1 - p(Z=1))\\\\\n",
    "&+ p(Y=1|X=0,Z=1)\\cdot (1 - p(X=1)) \\cdot p(Z=1) + p(Y=1|X=0,Z=0)\\cdot (1 - p(X=1))\\cdot (1 - p(Z=1)) \\\\\n",
    "&= 0.99 \\cdot 0.05 \\cdot 0.2 + 0.9 \\cdot 0.05 \\cdot 0.8 + 0.7 \\cdot 0.95 \\cdot 0.2 + 0.0001 \\cdot 0.95 \\cdot 0.8 = 0.178976\\\\\n",
    "p(Y=0) &= (1 - p(Y=1)) = (1 - 0.178976) = 0.821024\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "joint -> y=0, x=0, z=0, 0.759924\n",
      "joint -> y=0, x=0, z=1, 0.057\n",
      "joint -> y=1, x=0, z=0, 7.6e-05\n",
      "joint -> y=1, x=0, z=1, 0.133\n",
      "joint -> y=0, x=1, z=0, 0.004\n",
      "joint -> y=0, x=1, z=1, 0.0001\n",
      "joint -> y=1, x=1, z=0, 0.036\n",
      "joint -> y=1, x=1, z=1, 0.0099\n",
      "marginal for y\n",
      "[0.821024 0.178976]\n",
      "marginal for x (given)\n",
      "[0.95 0.05]\n",
      "marginal for z (given)\n",
      "[0.8 0.2]\n"
     ]
    }
   ],
   "source": [
    "# Calculate marginals for all the variables\n",
    "# this is not required, we added it for double checking\n",
    "import numpy as np\n",
    "px = np.array([0.95, 0.05])\n",
    "pz = np.array([0.80, 0.20])\n",
    "m = np.array([\n",
    "    [\n",
    "        [0.9999, 0.3], \n",
    "        [0.1, 0.01]\n",
    "    ], \n",
    "    [\n",
    "        [0.0001, 0.7], \n",
    "        [0.9, 0.99]\n",
    "    ]\n",
    "])\n",
    "\n",
    "joint = np.zeros(m.shape)\n",
    "for x in range(2):\n",
    "    for y in range(2):\n",
    "        for z in range(2):\n",
    "            joint[x, y, z] = np.around(m[y, x, z] * px[x]*pz[z], 6)\n",
    "            print(f'joint -> y={y}, x={x}, z={z}, {joint[x, y, z]}')\n",
    "print('marginal for y')\n",
    "print(joint.sum(0).sum(1))\n",
    "print('marginal for x (given)')\n",
    "print(joint.sum(1).sum(1))\n",
    "print('marginal for z (given)')\n",
    "print(joint.sum(0).sum(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "28024e8ece48c8b90f5fa0fdd0c13a6e",
     "grade": false,
     "grade_id": "cell-9afa97d7b5de1808",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Part 1: The sum-product algorithm\n",
    "\n",
    "We will implement a data structure to store a factor graph and to facilitate computations on this graph. Recall that a factor graph consists of two types of nodes, factors and variables. Below you will find some classes for these node types to get you started. Carefully inspect this code and make sure you understand what it does. Step by step will update its functionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "1affe45cf22d1acc9bdb7a7994d07c6a",
     "grade": false,
     "grade_id": "cell-70843f7b3041a417",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "class Node(object):\n",
    "    \"\"\"\n",
    "    Base-class for Nodes in a factor graph. Only instantiate sub-classes of Node.\n",
    "    \"\"\"\n",
    "    def __init__(self, name):\n",
    "        # A name for this Node, for printing purposes\n",
    "        self.name = name\n",
    "        \n",
    "        # Neighbours in the graph, identified with their index in this list.\n",
    "        # i.e. self.neighbours contains neighbour 0 through len(self.neighbours) - 1.\n",
    "        self.neighbours = []\n",
    "        \n",
    "        # Reset the node-state (not the graph topology)\n",
    "        self.reset()\n",
    "        \n",
    "    def reset(self):\n",
    "        # Incomming messages; a dictionary mapping neighbours to messages.\n",
    "        # That is, it maps  Node -> np.ndarray.\n",
    "        self.in_msgs = {}\n",
    "        \n",
    "        # A set of neighbours for which this node has pending messages.\n",
    "        # We use a python set object so we don't have to worry about duplicates.\n",
    "        self.pending = set([])\n",
    "\n",
    "    def add_neighbour(self, nb):\n",
    "        self.neighbours.append(nb)\n",
    "\n",
    "    def send_sp_msg(self, other):\n",
    "        # To be implemented in subclass.\n",
    "        raise Exception('Method send_sp_msg not implemented in base-class Node')\n",
    "   \n",
    "    def send_ms_msg(self, other):\n",
    "        # To be implemented in subclass.\n",
    "        raise Exception('Method send_ms_msg not implemented in base-class Node')\n",
    "    \n",
    "    def receive_msg(self, other, msg):\n",
    "        # Store the incomming message, replacing previous messages from the same node\n",
    "        self.in_msgs[other] = msg\n",
    "\n",
    "        # TODO: add pending messages\n",
    "        # self.pending.update(...)\n",
    "    \n",
    "    def __str__(self):\n",
    "        # This is printed when using 'print node_instance'\n",
    "        return self.name\n",
    "\n",
    "\n",
    "class Variable(Node):\n",
    "    def __init__(self, name, num_states):\n",
    "        \"\"\"\n",
    "        Variable node constructor.\n",
    "        Args:\n",
    "            name: a name string for this node. Used for printing. \n",
    "            num_states: the number of states this variable can take.\n",
    "            Allowable states run from 0 through (num_states - 1).\n",
    "            For example, for a binary variable num_states=2,\n",
    "            and the allowable states are 0, 1.\n",
    "        \"\"\"\n",
    "        self.num_states = num_states\n",
    "        \n",
    "        # Call the base-class constructor\n",
    "        super(Variable, self).__init__(name)\n",
    "    \n",
    "    def set_observed(self, observed_state):\n",
    "        \"\"\"\n",
    "        Set this variable to an observed state.\n",
    "        Args:\n",
    "            observed_state: an integer value in [0, self.num_states - 1].\n",
    "        \"\"\"\n",
    "        # Observed state is represented as a 1-of-N variable\n",
    "        # Set all-but-one states to a very low probability;\n",
    "        # (a bit of a hack to avoid -inf values when taking logs)\n",
    "        self.observed_state[:] = 1e-10\n",
    "        self.observed_state[observed_state] = 1.0\n",
    "        \n",
    "    def set_latent(self):\n",
    "        \"\"\"\n",
    "        Erase an observed state for this variable and consider it latent again.\n",
    "        \"\"\"\n",
    "        # No state is preferred, so set all entries of observed_state to 1.0\n",
    "        # Using this representation we need not differentiate observed an latent\n",
    "        # variables when sending messages.\n",
    "        self.observed_state[:] = 1.0\n",
    "        \n",
    "    def reset(self):\n",
    "        super(Variable, self).reset()\n",
    "        self.observed_state = np.ones(self.num_states)\n",
    "        \n",
    "    def marginal(self):\n",
    "        \"\"\"\n",
    "        Compute the marginal distribution of this Variable.\n",
    "        It is assumed that message passing has completed when this function is called.\n",
    "        \"\"\"\n",
    "        # TODO: compute marginal\n",
    "        return None\n",
    "    \n",
    "    def unnormalized_log_marginal(self):\n",
    "        \"\"\"\n",
    "        Compute the unnormalized log marginal distribution of this Variable.\n",
    "        It is assumed that message passing has completed when this function is called.\n",
    "        \"\"\"\n",
    "        # TODO: compute unnormalized log marginal\n",
    "        return None\n",
    "    \n",
    "    def send_sp_msg(self, other):\n",
    "        # TODO: implement Variable -> Factor message for sum-product\n",
    "        pass\n",
    "   \n",
    "    def send_ms_msg(self, other):\n",
    "        # TODO: implement Variable -> Factor message for max-sum\n",
    "        pass\n",
    "\n",
    "class Factor(Node):\n",
    "    def __init__(self, name, f, neighbours):\n",
    "        \"\"\"\n",
    "        Factor node constructor.\n",
    "        Args:\n",
    "            name: a name string for this node. Used for printing\n",
    "            f: a numpy.ndarray with N axes, where N is the number of neighbours.\n",
    "               That is, the axes of f correspond to variables, and the index along that axes corresponds to a value of that variable.\n",
    "               Each axis of the array should have as many entries as the corresponding neighbour variable has states.\n",
    "            neighbours: a list of neighbouring Variables. Bi-directional connections are created.\n",
    "        \"\"\"\n",
    "        # Call the base-class constructor\n",
    "        super(Factor, self).__init__(name)\n",
    "        \n",
    "        f = np.array(f)  # For convenience, accept lists (of lists (of numbers or potentially lists of ...))\n",
    "        \n",
    "        assert len(neighbours) == f.ndim, 'Factor function f should accept as many arguments as this Factor node has neighbours'\n",
    "        \n",
    "        for nb_ind in range(len(neighbours)):\n",
    "            nb = neighbours[nb_ind]\n",
    "            assert f.shape[nb_ind] == nb.num_states, 'The range of the factor function f is invalid for input %i %s' % (nb_ind, nb.name)\n",
    "            self.add_neighbour(nb)\n",
    "            nb.add_neighbour(self)\n",
    "\n",
    "        self.f = f\n",
    "        \n",
    "    def send_sp_msg(self, other):\n",
    "        # TODO: implement Factor -> Variable message for sum-product\n",
    "        pass\n",
    "   \n",
    "    def send_ms_msg(self, other):\n",
    "        # TODO: implement Factor -> Variable message for max-sum\n",
    "        pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "24aa9afccd84106fdc0727236c1cb186",
     "grade": false,
     "grade_id": "cell-e0ada624a5583fce",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 1.1. Initialize the graph (5 points)\n",
    "The equations in Part 0 can be represented by a factor graph. Instantiate this graph by creating Variable and Factor instances and linking them according to the graph structure. \n",
    "To instantiate the factor graph, first create the Variable nodes and then create Factor nodes, passing a list of neighbour Variables to each Factor. To get you started, we initialize the variable node for $X$ and the factor node corresponding to the prior $p(X)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "e8e1e5033150190f7ba288f5a2b004b4",
     "grade": false,
     "grade_id": "cell-85011c26139bed2e",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "X = Variable(name='X', num_states=2)\n",
    "X_prior = Factor(name='p(X)',\n",
    "                 f=np.array([0.95, 0.05]),\n",
    "                 neighbours=[X])\n",
    "\n",
    "# Please stick to the naming convention used below, otherwise the test functionality throughout the lab won't work\n",
    "Z = Variable(name='Z', num_states=2)\n",
    "Z_prior = Factor(name='p(Z)',\n",
    "                 f=np.array([0.80, 0.20]),\n",
    "                 neighbours=[Z])\n",
    "                 \n",
    "Y = Variable(name='Y', num_states=2)\n",
    "Y_cond = Factor(name='p(Y | X, Z)',\n",
    "                 f=np.array([\n",
    "                 [[0.9999, 0.3], [0.1, 0.01]],\n",
    "                 [[0.0001, 0.7], [0.9, 0.99]],\n",
    "                 ]),\n",
    "                 neighbours=[Y, X, Z])  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "1d247840a0cad095db98f369d0531275",
     "grade": true,
     "grade_id": "cell-1bf4d184d405f181",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "assert Z_prior.f.shape == (2,)\n",
    "assert Y_cond.f.shape == (2, 2, 2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "09f3f672778ffb77e22ba8702e9a2911",
     "grade": false,
     "grade_id": "cell-af817385e2927300",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "We will be doing a lot of marginalizations, i.e. obtain the distribution of a *single* variable from a distribution over multiple variables. Let's first write a function `marginalize` for that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "d81924d48029a907c55cea627ae5b56f",
     "grade": false,
     "grade_id": "cell-74ecd6f85f99a2d5",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.57354891 0.42645109]\n",
      "[0.42323629 0.20899062 0.36777309]\n",
      "[0.23683897 0.16852547 0.3175838  0.27705176]\n"
     ]
    }
   ],
   "source": [
    "# To marginalize, we sum over all variables except for the one, that we are interested in\n",
    "def marginalize(P, dim):\n",
    "    n_dims = len(P.shape)\n",
    "    indexes = tuple(i for i in range(n_dims) if i != dim)\n",
    "    return np.sum(P, indexes)\n",
    "\n",
    "# Lets try it\n",
    "test_P = np.random.rand(2, 3, 4)\n",
    "test_P = test_P / test_P.sum()  # Normalize for proper distribution\n",
    "\n",
    "# Do the marginal distributions look like you expect?\n",
    "print (marginalize(test_P, 0))\n",
    "print (marginalize(test_P, 1))\n",
    "print (marginalize(test_P, 2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "83cf7c2b9bf87fc4a0268c36436195f6",
     "grade": false,
     "grade_id": "cell-56015dc94e6caaf4",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 1.2 Factor to variable messages (20 points)\n",
    "Write a method `send_sp_msg(self, other)` for the Factor class, that checks if all the information required to pass a message to `Variable` `other` is present, computes the message and sends it to `other`. \"Sending\" here simply means calling the `receive_msg` function of the receiving node (we will implement this later). The message itself should be represented as a numpy array (np.array) whose length is equal to the number of states of the variable.\n",
    "\n",
    "In the very end of 1.2 below we overwrite `send_sp_msg(self, other)` for the Factor class. In general, this is considered bad practise but in this lab it saves us from scrolling up and down all the time.\n",
    "\n",
    "You will implement a function `send_sp_msg` that sends a message from a factor to a variable for the max-sum algorith. This function implements Equation 8.66 from Bishop. The message should be a numpy array (np.array) whose length is equal to the number of states of the variable.\n",
    "\n",
    "It is a good idea to write a number of helper functions to implement these equations.\n",
    "\n",
    "Since it is always a good idea to include checks, you should first write a method `can_send_message` that checks whether a node `node` has all the information required to pass a message to `other`. This should work for both variable and factor nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "8b75b5269753ba1861f6c7dbd41fa1a5",
     "grade": false,
     "grade_id": "cell-4164822e95a6b68a",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# We check if all the messages for a receiver have been received\n",
    "# If there is at least one message missing - we return False\n",
    "def can_send_message(sender, receiver):\n",
    "    for neighbour in sender.neighbours:\n",
    "        if neighbour == receiver:\n",
    "            continue\n",
    "        if neighbour not in sender.in_msgs.keys():\n",
    "            return False\n",
    "    \n",
    "    return True\n",
    "\n",
    "\n",
    "\n",
    "# Do the results make sense?\n",
    "print (can_send_message(X, X_prior))\n",
    "print (can_send_message(X_prior, X))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "15c446225bb0473c69caa2219f4c80bd",
     "grade": false,
     "grade_id": "cell-816ded386be1b608",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "\n",
    "In Eq. 8.66, Bishop writes $f(x, x_1, ..., x_M)$, where $x$ corresponds to the variable that will receive the message. For now, assume this variable is the `index`-th neighbour of the factor. In order to ease implementation, it may be a good idea to write a function that rearanges the dimensions of `f` to comply with this notation, i.e. moves the dimension `index` to the front. Make sure to return a copy and keep all other dimensions in order! Use `np.moveaxis`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "296bbf4e2783a797fe753c30c7ad7877",
     "grade": false,
     "grade_id": "cell-702210e632c423bd",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def move_dimension_first(f, index):\n",
    "    moved = np.moveaxis(f, index, 0)\n",
    "    return moved\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should calculate the product of the incoming messages of all neighbours of the sending node except the receiving node. Therefore it may be useful to write a function `get_neighbour_messages` that gathers these messages in a list. If you want to not make things complicated, make sure the order of the messages in the list corresponds to the order of the variables in `neighbours`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "aec4aebe3b5cc9ee35ae118c1763ff26",
     "grade": false,
     "grade_id": "cell-9e7f7986263e1543",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# We get messages from all the neighbours except for the receiver\n",
    "# We do an additional check to ensure that the message from a specified neighbour is actually in the list of incoming messages\n",
    "def get_neighbour_messages(sender, receiver):\n",
    "    l = [sender.in_msgs[n] for n in sender.neighbours if n != receiver and n in sender.in_msgs]\n",
    "    return l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before marginalizing, we need to calculate $\\prod_{m\\in\\text{ne}(f_s)\\setminus x} \\mu_{x_m\\rightarrow f_s}(x_m)$ (Eq. 8.66) for all possible combinations of $x_1, ..., x_M$ (values of the neighbour nodes except the receiving node). An elegant and efficient way to calculate these is using the n-way outer product of vectors. This product takes n vectors $\\mathbf{x}^{(1)}, \\ldots, \\mathbf{x}^{(n)}$ and computes a $n$-dimensional tensor (ndarray) whose element $i_0,i_1,...,i_n$ is given by $\\prod_j \\mathbf{x}^{(j)}_{i_j}$. In python, this is realized as `np.multiply.reduce(np.ix_(*vectors))` for a python list `vectors` of 1D numpy arrays. Try to figure out how this statement works -- it contains some useful functional programming techniques. What should `vectors` be? Try to see the link between the result and Eq. 8.66."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "e5e5d3c0d6252b785451adf1c41bd5fc",
     "grade": false,
     "grade_id": "cell-67269e06a554035c",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def calc_other_neighbour_msg_prod(sender, receiver):\n",
    "    messages = get_neighbour_messages(sender, receiver)\n",
    "    prod = np.multiply.reduce(np.ix_(*messages))\n",
    "    \n",
    "    return prod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "00442c496ed4ebda7eff1672a64c9f59",
     "grade": false,
     "grade_id": "cell-d5366f0c524b30d0",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Following Eq. 8.66, before marginalizing, you should calculate the product of $f(x, x_1, ..., x_M)$ with $\\prod_{m\\in\\text{ne}(f_s)\\setminus x} \\mu_{x_m\\rightarrow f_s}(x_m)$ for all configurations of $x, x_1, ..., x_M$. Since the second part does not depend on $x$, its tensor representations are of different dimensions. You can overcome this problem by using a loop, but preferably use numpy broadcasting by first aligning the dimensions of the tensors. You can use `np.expand_dims` or `X[None, ...]` to insert one dimension at the front. Write a function `calculate_factor` that, given `f` (which is reordered such that $x$ corresponds to the first dimension) and the (outer) product of the other neighbour messages, computes $f(x, x_1, ..., x_M) \\prod_{m\\in\\text{ne}(f_s)\\setminus x} \\mu_{x_m\\rightarrow f_s}(x_m)$ for all configurations of $x, x_1, ..., x_M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "3ab59a95b589ff71451036c73b289138",
     "grade": false,
     "grade_id": "cell-d34359928f3a997d",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def calculate_factor(f_neighb_first, neighbour_msg_prod):\n",
    "    # Because of the broadcasting rules of numpy, we don't necessarily need to expand dimensions before multiplication\n",
    "    factor = f_neighb_first * neighbour_msg_prod\n",
    "    return factor\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "ef9e90a515d362d281f313955fbbc75f",
     "grade": false,
     "grade_id": "cell-d9b12fad1b574924",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Put all the pieces together to define a function `calc_sum_product_factor_to_variable_msg` that calculates Eq. 8.66."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "f3a2dc41ed909c8d0eb88b213ce66a99",
     "grade": false,
     "grade_id": "cell-ec6b421860d40ae0",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def calc_sum_product_factor_to_variable_msg(factor, variable):\n",
    "    var_index = factor.neighbours.index(variable)\n",
    "    f_neighb_first = move_dimension_first(factor.f, var_index)\n",
    "    neighbour_msg_prod = calc_other_neighbour_msg_prod(factor, variable)\n",
    "\n",
    "    f2 = calculate_factor(f_neighb_first, neighbour_msg_prod)\n",
    "    f = marginalize(f2, 0)\n",
    "    return f\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "c65b2bb803c0f8b53c5eb427e8bb2906",
     "grade": false,
     "grade_id": "cell-7d334ab3a40918b0",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Finally, we will define the send message function for you\n",
    "# Note, that we have edited the cell in order to return message \n",
    "# It is being printed during backward and forward pass, which greatly simplified the debugging procedure for us\n",
    "def factor_send_sp_msg(self, variable):\n",
    "    \n",
    "    assert isinstance(variable, Variable), \"Factor can only send messages to variable!\"\n",
    "    assert can_send_message(self, variable), \"Cannot send message!\"\n",
    "    \n",
    "    out_msg = calc_sum_product_factor_to_variable_msg(self, variable)\n",
    "    \n",
    "    # Send the message\n",
    "    variable.receive_msg(self, out_msg)\n",
    "    \n",
    "    # Remove the pending sign if present\n",
    "    self.pending.discard(variable)\n",
    "    return out_msg\n",
    "    \n",
    "Factor.send_sp_msg = factor_send_sp_msg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "4ddb724608cd7205b51cbadda1a3664f",
     "grade": true,
     "grade_id": "cell-429fba08a94a3bd2",
     "locked": true,
     "points": 20,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "# message from X_prior to X\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "\n",
    "X_prior.send_sp_msg(X)\n",
    "assert np.allclose(list(X.in_msgs.values()), [0.95, 0.05])\n",
    "\n",
    "# message from Z_prior to Z\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "\n",
    "Z_prior.send_sp_msg(Z)\n",
    "assert np.allclose(list(Z.in_msgs.values()), [0.8, 0.2])\n",
    "\n",
    "# message from Y_cond to Y\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "Y_cond.receive_msg(X, X_prior.f) # simulating that Y_cond received all necessary messages from X\n",
    "Y_cond.receive_msg(Z, Z_prior.f) # simulating that Y_cond received all necessary messages from Z\n",
    "Y_cond.send_sp_msg(Y)\n",
    "assert np.allclose(list(Y.in_msgs.values()), [0.821024, 0.178976])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "2fd488d52a30d4081cdbc910c280a279",
     "grade": false,
     "grade_id": "cell-2ac03f39f1b1281c",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 1.3 Variable to factor messages (10 points)\n",
    "\n",
    "Write a method `calc_sum_product_variable_to_factor_msg(variable, factor)` that computes the message to be sent to a neighbour variable by a factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "5c3c49af77b02653acaac2455d864132",
     "grade": false,
     "grade_id": "cell-93d0a0d839eb99f3",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def calc_sum_product_variable_to_factor_msg(variable, factor):\n",
    "    if len(variable.neighbours) == 1:\n",
    "        message = np.ones(variable.num_states)\n",
    "    else:\n",
    "        messages = get_neighbour_messages(variable, factor)\n",
    "        # we can do elementwise product because the messages to a variable have the same dimension\n",
    "        message = np.prod(messages, axis=0)\n",
    "\n",
    "    return message\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "bb0a564a418df48edc0a627ff3df1bc4",
     "grade": false,
     "grade_id": "cell-8ce3afc9581048ad",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Finally, we will define the send message function for you\n",
    "def variable_send_sp_msg(self, factor):\n",
    "    \n",
    "    assert isinstance(factor, Factor), \"Variable can only send messages to factor!\"\n",
    "    assert can_send_message(self, factor), \"Cannot send message!\"\n",
    "    \n",
    "    out_msg = calc_sum_product_variable_to_factor_msg(self, factor)\n",
    "    \n",
    "    # Send the message\n",
    "    factor.receive_msg(self, out_msg)\n",
    "    \n",
    "    # Remove the pending sign if present\n",
    "    self.pending.discard(factor)\n",
    "    return out_msg\n",
    "    \n",
    "Variable.send_sp_msg = variable_send_sp_msg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "80a738a13d17bc6014f56719bc935ca4",
     "grade": true,
     "grade_id": "cell-0cbe0bfa635450f3",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "# First message from X to Y_cond\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "X_prior.send_sp_msg(X) # simulating that X received all necessary messages\n",
    "X.send_sp_msg(Y_cond) \n",
    "assert np.allclose(list(Y_cond.in_msgs.values()), [0.95, 0.05])\n",
    "\n",
    "# Second message from Z to Y_cond\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Z_prior.send_sp_msg(Z) # simulating that Z received all necessary messages\n",
    "Z.send_sp_msg(Y_cond)\n",
    "assert np.allclose(list(Y_cond.in_msgs.values()), [[0.95, 0.05], [0.8, 0.2]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "72a4f47ec3cd02e4e32b7d68e9d0692b",
     "grade": false,
     "grade_id": "cell-6346fc68da031348",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### Testing a single forward pass and a single backward pass\n",
    "Before we go on we will make sure that messages can be passed in both directions the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d50a8fd57d73f85f1c8f81be1ba13f69",
     "grade": false,
     "grade_id": "cell-f44819044a283ec9",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "X_prior.reset()\n",
    "X.reset()\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "# Forward pass\n",
    "X_prior.send_sp_msg(X)\n",
    "Z_prior.send_sp_msg(Z)\n",
    "X.send_sp_msg(Y_cond)\n",
    "Z.send_sp_msg(Y_cond)\n",
    "Y_cond.send_sp_msg(Y)\n",
    "assert np.allclose(list(Y.in_msgs.values()), [0.821024, 0.178976])\n",
    "\n",
    "# Backward pass\n",
    "Y.send_sp_msg(Y_cond)\n",
    "Y_cond.send_sp_msg(X)\n",
    "Y_cond.send_sp_msg(Z)\n",
    "X.send_sp_msg(X_prior)\n",
    "Z.send_sp_msg(Z_prior)\n",
    "assert np.allclose(list(X.in_msgs.values()), [[0.95, 0.05],[1., 1.]])\n",
    "assert np.allclose(list(Z.in_msgs.values()), [[0.8, 0.2],[1., 1.]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "27ce20145a58d03708168b08fbbfeadd",
     "grade": false,
     "grade_id": "cell-f1b2c01d81e8d007",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 1.4 Compute marginal (10 points)\n",
    "Later in this assignment, we will implement message passing schemes to do inference. Once the message passing has completed, we will want to compute local marginals for each variable.\n",
    "Write the method `marginal` for the `Variable` class, that computes a marginal distribution over that node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "8002994abb8ad6ada1ad19a0b433c273",
     "grade": false,
     "grade_id": "cell-0f6434f5c5d4dc9c",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def marginal(self):\n",
    "    marginal = np.multiply.reduce(list(self.in_msgs.values()), 0)\n",
    "    \n",
    "    return marginal\n",
    "\n",
    "Variable.marginal = marginal\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "8dc0f9a18f2801cf7f3182732f45b52e",
     "grade": true,
     "grade_id": "cell-7f666a68bbbe81e5",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "# Simulate a single forward pass\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "X_prior.send_sp_msg(X)\n",
    "Z_prior.send_sp_msg(Z)\n",
    "X.send_sp_msg(Y_cond)\n",
    "Z.send_sp_msg(Y_cond)\n",
    "Y_cond.send_sp_msg(Y)\n",
    "\n",
    "assert np.allclose(X.marginal(), [0.95, 0.05])\n",
    "assert np.allclose(Z.marginal(), [0.8, 0.2])\n",
    "assert np.allclose(Y.marginal(), [0.821024, 0.178976])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "3d2d9dd5139b1ed6891107868305e57b",
     "grade": false,
     "grade_id": "cell-e38ca86d50f6c9c1",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 1.5 Receiving messages (10 points)\n",
    "In order to implement a message passing algorithms, we need some way to determine which nodes are ready to send messages to which neighbours. We make use of the concept of \"pending messages\", which is explained in Bishop (8.4.7): \n",
    "\"we will say that a (variable or factor)\n",
    "node a has a message pending on its link to a node b if node a has received any\n",
    "message on any of its other links since the last time it sent a message to b. Thus,\n",
    "when a node receives a message on one of its links, this creates pending messages\n",
    "on all of its other links.\"\n",
    "\n",
    "Before we say node a has a pending message for node b, we **must check that node a has received all messages needed to compute the message that is to be sent to b**.\n",
    "\n",
    "Modify the function `receive_msg`, so that it updates the self.pending variable as described above. The member self.pending is a set that is to be filled with Nodes to which self has pending messages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "d4256320e0a6fff7550f0a0c744ddb61",
     "grade": false,
     "grade_id": "cell-39b7b3a84cdf1e56",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# ANSWER 1.5\n",
    "def receive_msg(self, other, msg):\n",
    "    self.in_msgs[other] = msg\n",
    "    for n in self.neighbours:\n",
    "        if can_send_message(self, n):\n",
    "            self.pending.add(n)\n",
    "\n",
    "Node.receive_msg = receive_msg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "05683e9bfc7be13ffac9b88806845c3b",
     "grade": true,
     "grade_id": "cell-42c321c98bf9b6b1",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "assert X_prior.pending == set()\n",
    "\n",
    "X_prior.pending.add(X)\n",
    "assert str(list(X_prior.pending)[0]) == X.name\n",
    "\n",
    "X_prior.send_sp_msg(X)\n",
    "assert X_prior.pending == set()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.6 Inference Engine (10 points)\n",
    "Write a function `sum_product(node_list)` that runs the sum-product message passing algorithm on a tree-structured factor graph with given nodes. The input parameter `node_list` is a list of all Node instances in the graph, which is assumed to be ordered correctly. That is, the list starts with a leaf node, which can always send a message. Subsequent nodes in `node_list` should be capable of sending a message when the pending messages of preceding nodes in the list have been sent. The sum-product algorithm then proceeds by passing over the list from beginning to end, sending all pending messages at the nodes it encounters. Then, in reverse order, the algorithm traverses the list again and again sends all pending messages at each node as it is encountered. For this to work, we initialized pending messages for all the leaf nodes, e.g. `X_prior.pending.add(X)`, where `X_prior` is a Factor node corresponding the the prior, `X` is a `Variable` node and the only connection of `X_prior` goes to `X`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "0295fa5b1b0d1c2dd532935bd543056f",
     "grade": false,
     "grade_id": "cell-6a6fb2568d840c68",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def sum_product(node_list):\n",
    "    print('Forward pass')\n",
    "    for i, node in enumerate(node_list):\n",
    "        # Iterate over all neighbours, which haven't received all necessary messages in the current pass\n",
    "        # we are using set between neighbours and nodes further in the message passing chain\n",
    "        # to make sure that we don't send messages back to nodes that have \n",
    "        # already \"forwarded\" messages in a current pass\n",
    "        for neigh in set(node.neighbours) & set(node_list[i:]):\n",
    "            if can_send_message(node, neigh):\n",
    "                msg = node.send_sp_msg(neigh)\n",
    "                print(f'{node} -> {neigh}: {msg}')\n",
    "    \n",
    "    rev_node_list = node_list[::-1]\n",
    "    print('\\nBackward pass')\n",
    "    for i, node in enumerate(rev_node_list):\n",
    "        # Again, iterate over all neighbours, which haven't received all necessary messages in the current pass\n",
    "        for neigh in set(node.neighbours) & set(rev_node_list[i:]):\n",
    "            if neigh in node.pending:\n",
    "                msg = node.send_sp_msg(neigh)\n",
    "                print(f'{node} -> {neigh}: {msg}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "261e257a05c618b4e1c5d9568f7403d5",
     "grade": true,
     "grade_id": "cell-fb8315046f8fde42",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(X) -> X: [0.95 0.05]\n",
      "X -> p(Y | X, Z): [0.95 0.05]\n",
      "p(Z) -> Z: [0.8 0.2]\n",
      "Z -> p(Y | X, Z): [0.8 0.2]\n",
      "p(Y | X, Z) -> Y: [0.821024 0.178976]\n",
      "\n",
      "Backward pass\n",
      "Y -> p(Y | X, Z): [1. 1.]\n",
      "p(Y | X, Z) -> X: [1. 1.]\n",
      "p(Y | X, Z) -> Z: [1. 1.]\n",
      "Z -> p(Z): [1. 1.]\n",
      "X -> p(X): [1. 1.]\n"
     ]
    }
   ],
   "source": [
    "### Test test test\n",
    "nodes = [X_prior, X, Z_prior, Z, Y_cond, Y]\n",
    "for n in nodes:\n",
    "    n.reset()\n",
    "    \n",
    "X_prior.pending.add(X)\n",
    "Z_prior.pending.add(Z)\n",
    "Y.pending.add(Y_cond)\n",
    "\n",
    "sum_product(nodes)\n",
    "assert np.allclose(Y.marginal(), [0.821024, 0.178976])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.7 Observed variables and probabilistic queries (15 points)\n",
    "We will now use the inference engine to answer probabilistic queries. That is, we will set certain variables to observed values, and obtain the marginals over latent variables. We have already provided functions `set_observed` and `set_latent` that manage a member of Variable called `observed_state`. Modify the `calc_sum_product_variable_to_factor_msg` and `Variable.marginal` routines that you wrote before, to use `observed_state` so as to get the required marginals when some nodes are observed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "da29573eaddaf62947074fc00f714df1",
     "grade": false,
     "grade_id": "cell-55c50fb9c1836789",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def calc_sum_product_variable_to_factor_msg(variable, factor):\n",
    "    if len(variable.neighbours) == 1:\n",
    "        message = np.ones(variable.num_states)\n",
    "    else:\n",
    "        messages = get_neighbour_messages(variable, factor)\n",
    "        message = np.prod(messages, axis=0)\n",
    "\n",
    "    return message * variable.observed_state\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d27bd6da3b20252943b2828ceb5bb472",
     "grade": true,
     "grade_id": "cell-13801e018b6f9d27",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Test, test, test\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "Y_cond.reset()\n",
    "\n",
    "X_prior.send_sp_msg(X)\n",
    "X.set_observed(0)\n",
    "X.send_sp_msg(Y_cond)\n",
    "assert np.allclose(list(Y_cond.in_msgs.values()), [9.5e-01, 5.0e-12])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "55a9272099aabff3579cac4f41f08973",
     "grade": false,
     "grade_id": "cell-d0b91da808874e73",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def marginal(self):\n",
    "    marginal = np.multiply.reduce(list(self.in_msgs.values()), 0) * self.observed_state\n",
    "    \n",
    "    return marginal / marginal.sum()\n",
    "    \n",
    "\n",
    "Variable.marginal = marginal\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "fa342104c5f9d55b65ad4c202475a2c5",
     "grade": true,
     "grade_id": "cell-2f8e5d9052ed73ef",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test, test, test\n",
    "# Simulate a single forward pass\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "X.set_observed(0)\n",
    "Z.set_observed(0)\n",
    "\n",
    "X_prior.send_sp_msg(X)\n",
    "Z_prior.send_sp_msg(Z)\n",
    "X.send_sp_msg(Y_cond)\n",
    "Z.send_sp_msg(Y_cond)\n",
    "Y_cond.send_sp_msg(Y)\n",
    "\n",
    "assert np.allclose(Y.marginal(), [9.99900000e-01, 1.00000022e-04])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "97cce5eb6885d02431cf74e2affba069",
     "grade": false,
     "grade_id": "cell-2277868e4ec25b1d",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Part 2: The max-sum algorithm\n",
    "Next, we implement the max-sum algorithm as described in section 8.4.5 of Bishop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "34cda2c1859ee99af59b1e6f6609b5b6",
     "grade": false,
     "grade_id": "cell-f1dbc98349fb824f",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 2.1 Factor to variable messages (10 points)\n",
    "Implement the function `Factor.send_ms_msg` that sends Factor -> Variable messages for the max-sum algorithm. It is analogous to the `Factor.send_sp_msg` function you implemented before. Make sure it works for observed and unobserved nodes. Consider using a number of helper functions as seen in Part 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ms_calc_other_neighbour_msg_sum(sender, receiver):\n",
    "    messages = get_neighbour_messages(sender, receiver)\n",
    "    sum_var = np.add.reduce(np.ix_(*messages))\n",
    "    return sum_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def maximize(P, dim):\n",
    "    n_dims = len(P.shape)\n",
    "    indexes = tuple(i for i in range(n_dims) if i != dim)\n",
    "    return np.max(P, indexes) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ms_calculate_factor(f_neighb_first, neighbour_msg_prod):\n",
    "    factor = np.log(f_neighb_first) + neighbour_msg_prod\n",
    "    return factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ms_calc_max_sum_factor_to_variable_msg(factor, variable):\n",
    "    var_index = factor.neighbours.index(variable)\n",
    "    f_neighb_first = move_dimension_first(factor.f, var_index)\n",
    "\n",
    "    neighbour_msg_prod = ms_calc_other_neighbour_msg_sum(factor, variable) + np.log(variable.observed_state)\n",
    "    f = ms_calculate_factor(f_neighb_first, neighbour_msg_prod)\n",
    "    return maximize(f, 0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "c16471807291ab66dbe760bddc2f1480",
     "grade": false,
     "grade_id": "cell-389621edefd825a3",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def factor_send_ms_msg(self, variable):\n",
    "    \n",
    "    assert isinstance(variable, Variable), \"Factor can only send messages to variable!\"\n",
    "    assert can_send_message(self, variable), \"Cannot send message!\"\n",
    "\n",
    "    out_msg = ms_calc_max_sum_factor_to_variable_msg(self, variable)\n",
    "\n",
    "    # Send the message\n",
    "    variable.receive_msg(self, out_msg)\n",
    "\n",
    "    # Remove the pending sign if present\n",
    "    self.pending.discard(variable)\n",
    "    return out_msg\n",
    "    \n",
    "Factor.send_ms_msg = factor_send_ms_msg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "36080ee2c75dc501b57f32f527c65656",
     "grade": true,
     "grade_id": "cell-c6156e1505633bc4",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "# message from X_prior to X\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "\n",
    "X_prior.send_ms_msg(X)\n",
    "assert np.allclose(list(X.in_msgs.values()), [-0.05129329, -2.99573227])\n",
    "\n",
    "# message from Z_prior to Z\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "\n",
    "Z_prior.send_ms_msg(Z)\n",
    "assert np.allclose(list(Z.in_msgs.values()), [-0.22314355, -1.60943791])\n",
    "\n",
    "# message from Y_cond to Y\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "Y_cond.receive_msg(X, X_prior.f) # simulating that Y_cond received all necessary messages from X\n",
    "Y_cond.receive_msg(Z, Z_prior.f) # simulating that Y_cond received all necessary messages from Z\n",
    "Y_cond.send_ms_msg(Y)\n",
    "assert np.allclose(list(Y.in_msgs.values()), [1.74989999, 0.79332506])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "cfbc002bdf7c6ff9b40c06ef10403001",
     "grade": false,
     "grade_id": "cell-6eff0dd8439243ca",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 2.2 Variable to factor messages (10 points)\n",
    "Implement the `Variable.send_ms_msg` function that sends Variable -> Factor messages for the max-sum algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ms_calc_max_sum_variable_to_factor_msg(variable, factor):\n",
    "    if len(variable.neighbours) == 1:\n",
    "        message = np.zeros(variable.num_states) \n",
    "    else:\n",
    "        messages = get_neighbour_messages(variable, factor)\n",
    "        message = np.sum(messages, axis=0)\n",
    "\n",
    "    return message + np.log(variable.observed_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "a8e8c05dffcf10717b1e43299b96b627",
     "grade": false,
     "grade_id": "cell-183948f06650a7e0",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def variable_send_ms_msg(self, factor):\n",
    "    out_msg = ms_calc_max_sum_variable_to_factor_msg(self, factor)\n",
    "    factor.receive_msg(self, out_msg)\n",
    "    \n",
    "    # Remove the pending sign if present\n",
    "    self.pending.discard(factor)\n",
    "    return out_msg\n",
    "\n",
    "Variable.send_ms_msg = variable_send_ms_msg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "b597b2dc5896b79f1fb39862d19f2ea4",
     "grade": true,
     "grade_id": "cell-bbf08d58e9b309c9",
     "locked": true,
     "points": 10,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "# First message from X to Y_cond\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "X_prior.send_ms_msg(X) # simulating that X received all necessary messages\n",
    "X.send_ms_msg(Y_cond)\n",
    "assert np.allclose(list(Y_cond.in_msgs.values()), [-0.05129329, -2.99573227])\n",
    "\n",
    "# Second message from Z to Y_cond\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Z_prior.send_ms_msg(Z) # simulating that Z received all necessary messages\n",
    "Z.send_ms_msg(Y_cond)\n",
    "assert np.allclose(list(Y_cond.in_msgs.values()), [[-0.05129329, -2.99573227], [-0.22314355, -1.60943791]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d94749234e759cb16863184e99919a57",
     "grade": false,
     "grade_id": "cell-6e1b3e4ee8013fcb",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 2.3 Implement unnormalized log marginal (5 points)\n",
    "Write the method `unnormalized_log_marginal` for the `Variable` class, that computes a unnormalized log marginal distribution over that node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "2a82b636f2eb1326ea6c7b4cce15f428",
     "grade": false,
     "grade_id": "cell-96de1de3d9b1da94",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def unnormalized_log_marginal(self):\n",
    "    marginal = np.add.reduce(list(self.in_msgs.values()), 0) + np.log(self.observed_state)\n",
    "    return marginal\n",
    "\n",
    "Variable.unnormalized_log_marginal = unnormalized_log_marginal\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "977c4965c9797bd48b917039e495038a",
     "grade": true,
     "grade_id": "cell-d0c08fd4b0857dde",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "### Test test test\n",
    "# Simulate a single forward pass\n",
    "X_prior.reset()\n",
    "X.reset()\n",
    "Z_prior.reset()\n",
    "Z.reset()\n",
    "Y_cond.reset()\n",
    "Y.reset()\n",
    "\n",
    "X_prior.send_ms_msg(X)\n",
    "Z_prior.send_ms_msg(Z)\n",
    "X.send_ms_msg(Y_cond)\n",
    "Z.send_ms_msg(Y_cond)\n",
    "Y_cond.send_ms_msg(Y)\n",
    "\n",
    "assert np.allclose(X.unnormalized_log_marginal(), [-0.05129329, -2.99573227])\n",
    "assert np.allclose(Z.unnormalized_log_marginal(), [-0.22314355, -1.60943791])\n",
    "assert np.allclose(Y.unnormalized_log_marginal(), [-0.27453685, -2.01740615])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "811de1e611d918a437d3436251ae8f56",
     "grade": false,
     "grade_id": "cell-eb540987e2fb70d6",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 2.4 Find a MAP state (10 points)\n",
    "\n",
    "Using the same message passing schedule we used for sum-product, implement the max-sum algorithm. For simplicity, we will ignore issues relating to non-unique maxima. So there is no need to implement backtracking; the MAP state is obtained by a per-node maximization (eq. 8.98 in Bishop). Make sure your algorithm works with both latent and observed variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "17fb320cc9db70730f6e7a6fd0702903",
     "grade": false,
     "grade_id": "cell-f94559550902ac16",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def max_sum(node_list):\n",
    "    print('Forward pass')\n",
    "    for i, node in enumerate(node_list):\n",
    "        for neigh in set(node.neighbours) & set(node_list[i:]):\n",
    "            if can_send_message(node, neigh):\n",
    "                msg = node.send_ms_msg(neigh)\n",
    "                print(f'{node} -> {neigh}: {np.exp(msg)}')\n",
    "    \n",
    "    rev_node_list = node_list[::-1]\n",
    "    print('\\nBackward pass')\n",
    "    for i, node in enumerate(rev_node_list):\n",
    "        for neigh in set(node.neighbours) & set(rev_node_list[i:]):\n",
    "            if neigh in node.pending:\n",
    "                msg = node.send_ms_msg(neigh)\n",
    "                print(f'{node} -> {neigh}: {np.exp(msg)}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "49674d2f737491907099ead2f661ab25",
     "grade": true,
     "grade_id": "cell-ac2d26dcb27a8d75",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(X) -> X: [0.95 0.05]\n",
      "X -> p(Y | X, Z): [0.95 0.05]\n",
      "p(Z) -> Z: [0.8 0.2]\n",
      "Z -> p(Y | X, Z): [0.8 0.2]\n",
      "p(Y | X, Z) -> Y: [0.759924 0.133   ]\n",
      "\n",
      "Backward pass\n",
      "Y -> p(Y | X, Z): [1. 1.]\n",
      "p(Y | X, Z) -> X: [0.79992 0.72   ]\n",
      "p(Y | X, Z) -> Z: [0.949905 0.665   ]\n",
      "Z -> p(Z): [0.949905 0.665   ]\n",
      "X -> p(X): [0.79992 0.72   ]\n"
     ]
    }
   ],
   "source": [
    "### Test test test: unobserved\n",
    "nodes = [X_prior, X, Z_prior, Z, Y_cond, Y]\n",
    "for n in nodes:\n",
    "    n.reset()\n",
    "    \n",
    "X_prior.pending.add(X)\n",
    "Z_prior.pending.add(Z)\n",
    "Y.pending.add(Y_cond)\n",
    "\n",
    "max_sum(nodes)\n",
    "\n",
    "assert np.allclose(Y.unnormalized_log_marginal(), [-0.27453685, -2.01740615] )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d00995d3a4fb74ce8ed9481349e9044d",
     "grade": true,
     "grade_id": "cell-f06b8b4bec85b4dd",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(X) -> X: [0.95 0.05]\n",
      "X -> p(Y | X, Z): [0.95 0.05]\n",
      "p(Z) -> Z: [8.e-11 2.e-01]\n",
      "Z -> p(Y | X, Z): [8.e-21 2.e-01]\n",
      "p(Y | X, Z) -> Y: [0.057 0.133]\n",
      "\n",
      "Backward pass\n",
      "Y -> p(Y | X, Z): [1. 1.]\n",
      "p(Y | X, Z) -> X: [0.14  0.198]\n",
      "p(Y | X, Z) -> Z: [0.045  0.0495]\n",
      "Z -> p(Z): [4.50e-12 4.95e-02]\n",
      "X -> p(X): [0.14  0.198]\n"
     ]
    }
   ],
   "source": [
    "### Test test test: partiallY observed\n",
    "nodes = [X_prior, X, Z_prior, Z, Y_cond, Y]\n",
    "for n in nodes:\n",
    "    n.reset()\n",
    "    \n",
    "X_prior.pending.add(X)\n",
    "Z_prior.pending.add(Z)\n",
    "Y.pending.add(Y_cond)\n",
    "\n",
    "Z.set_observed(1)\n",
    "\n",
    "max_sum(nodes)\n",
    "\n",
    "assert np.allclose(Y.unnormalized_log_marginal(), [-2.86470401, -2.01740615])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "6e5b7cb4dc60f898345d08a07889f32c",
     "grade": false,
     "grade_id": "cell-f79ce013522649ed",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Given the max-marginals what do you have to do in order to find the global optimum? Why can we neglect the normalization constant here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "1bd8a1937fd5e81b184e36bc4c4b1b68",
     "grade": true,
     "grade_id": "cell-605e729ee910dbc8",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "\n",
    "In order to find the global optimum, as mentioned in Bishop, sec. 8.4.5., we would have to do following:\n",
    "- Keep track of which values of the variables gave rise to the maximum state of each variable\n",
    "- Once we know the most probable value of the ﬁnal node $x_N$,  simply follow the link back to ﬁnd the most probable state of node $x_{N−1}$ and so on back to the initial node $x_1$. This procedure is known as back-tracking.\n",
    "\n",
    "We can neglect the normalization constant, because we care for the maximum probability and $max (\\frac{\\mathbf{x}}{const}) = \\frac{max(\\mathbf{x})}{const}$, which means, that normalization constant won't change the final result, it will just change it proportionally.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "80bbe06f90fa8e5f278119c02cf5081d",
     "grade": false,
     "grade_id": "cell-d0592bafd25122c0",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### Part 3: Medical graph\n",
    "Now that we implemented the sum-product and max-sum algorithm. We will apply them to a poly-tree structured medical diagnosis example.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "9c99fffc963ea4db625082689e3b7c55",
     "grade": false,
     "grade_id": "cell-f0022d301754179a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 3.1 Initialize the graph (5 points)\n",
    "\n",
    "Convert the directed graphical model (\"Bayesian Network\") shown below to a factor graph. Instantiate this graph by creating Variable and Factor instances and linking them according to the graph structure. \n",
    "To instantiate the factor graph, first create the Variable nodes and then create Factor nodes, passing a list of neighbour Variables to each Factor.\n",
    "Use the following prior and conditional probabilities.\n",
    "\n",
    "$$\n",
    "p(\\verb+Influenza+) = 0.05 \\\\\\\\\n",
    "p(\\verb+Smokes+) = 0.2 \\\\\\\\\n",
    "$$\n",
    "$$\n",
    "p(\\verb+SoreThroat+ = 1 | \\verb+Influenza+ = 1) = 0.3 \\\\\\\\\n",
    "p(\\verb+SoreThroat+ = 1 | \\verb+Influenza+ = 0) = 0.001 \\\\\\\\\n",
    "p(\\verb+Fever+ = 1| \\verb+Influenza+ = 1) = 0.9 \\\\\\\\\n",
    "p(\\verb+Fever+ = 1| \\verb+Influenza+ = 0) 0.05 \\\\\\\\\n",
    "p(\\verb+Bronchitis+ = 1 | \\verb+Influenza+ = 1, \\verb+Smokes+ = 1) = 0.99 \\\\\\\\\n",
    "p(\\verb+Bronchitis+ = 1 | \\verb+Influenza+ = 1, \\verb+Smokes+ = 0) = 0.9 \\\\\\\\\n",
    "p(\\verb+Bronchitis+ = 1 | \\verb+Influenza+ = 0, \\verb+Smokes+ = 1) = 0.7 \\\\\\\\\n",
    "p(\\verb+Bronchitis+ = 1 | \\verb+Influenza+ = 0, \\verb+Smokes+ = 0) = 0.0001 \\\\\\\\\n",
    "p(\\verb+Coughing+ = 1| \\verb+Bronchitis+ = 1) = 0.8 \\\\\\\\\n",
    "p(\\verb+Coughing+ = 1| \\verb+Bronchitis+ = 0) = 0.07 \\\\\\\\\n",
    "p(\\verb+Wheezing+ = 1| \\verb+Bronchitis+ = 1) = 0.6 \\\\\\\\\n",
    "p(\\verb+Wheezing+ = 1| \\verb+Bronchitis+ = 0) = 0.001 \\\\\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we show the corresponding factor graph"
   ]
  },
  {
   "attachments": {
    "graph.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABWkAAAHfCAYAAAAiMVqkAAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAAsTAAALEwEAmpwYAAAAB3RJTUUH4gUJEyAgkjrJkAAAABl0RVh0Q29tbWVudABDcmVhdGVkIHdpdGggR0lNUFeBDhcAACAASURBVHja7N1/aFv3vf/xly+GK7iCCuYxl3nggEcUSSUKU6jDEubw9WjKHJJSH9XlJsy9S9lgCTeGBBxug45ZwIEEHHBGGB60Y4G2Pg7xiEtziYcNdUkucYlCLEUjGjGbxtVAA/3hPwxXQ98/bCmWLDm2Y/vox/MBoqklHZ3z1vl8zvu89Tmf05DNZrMCAADAuqXTaYXDYQIBYMP8fr9cLheBAAAABRoJAQAAwMaEw2EdPnyYQADYsKmpKXV0dBAIAABQgCItAADAJg0NDcnv9xMIAC8VDofV19dHIAAAQEkUaQEAADbJ7/czIg4AAADAK/sXQgAAAAAAAAAA9qFICwAAAAAAAAA2okgLAAAAAAAAADaiSAsAAAAAAAAANqJICwAAAAAAAAA2okgLAAAAAAAAADaiSAsAAAAAAAAANqJICwAAAAAAAAA2okgLAAAAAAAAADaiSAsAAAAAAAAANqJICwAAAAAAAAA2okgLAAAAAAAAADaiSAsAAAAAAAAANqJICwAAgKoRNb1q8JqKFvzVUrChQQ0NQVlb8BlWsEENDblHUNYWL7+Koi3T2yCvGWXHAwAA2GYUaQEAAFAdoqaMgagM05Rnmz7CCjYoaBHqJR6ZpqHogCHqtAAAANuLIi0AAACqgmUOKCpDhrFtnyDLkiSPQpGsstmsstlRGfUcdMOQoagGTCrXAAAA24kiLQAAAKrAcgHVMHagaOqV10PElywXxS1LlGkBAAC2D0VaAAAA2GDlPK9Lc58WzgNb/PKlIqHH493kZyzPZ9vw4nNeTGtQPOds7v/Xmod2rXlq13jOChasQ0NDgxpKzq+w3nUvfn25x8p1KY53gxoavGWnNPB6PHoxyhgAAADbgSItAAAAbGUFvRooKBBaChYVOC3LkuSREfRs8jMa5B2IrvrbTs4/awXLFGSt4Jo3JduSdfd45M3HtjjekhTVgLf0Mj1BQ578dwAAAIDtQJEWAAAANrJkWYZGs8VzwFoy80M7o4pEtvIzshpdnjPBMk1FlXsu99ml1ufVRE3vcgG0cD2y2YhCnuLt3ci65xQvN6ts7oXyKGStuNmaJ6RImdcWLrNIJCLuHwYAALA9KNICAADAVsboykKoodHlgmHUGl0uCkYUjUqbnyvWo1CksNhqmKGlomU0qsi2b2FUo1a05HpIHpnLf3uxvVux7paCy8NiPSFLpmdFfC3JKDflQqllerxLo3B3JFYAAAD1iSItAAAAbLR8Y6qCPy3fHGzLioIliru5wuOOyBWZl6YUKDtfbMnt3cy6R2V6l5fpCckyVyzACqrBO8CIWAAAgApDkRYAAACoIS/m+DU0GlkxzYGiMs3c6NpI4XQHWzStAwAAADaHIi0AAABsZGnV/agsa3kUaO5mV155PJIUUaTihoCWWKdopGhEbG79PQpFiuaC3epCqRXM3/yrcBqJpXWNLhdvTdPzknUu8Vz++wAAAMBWo0gLAAAAW1nB5Uvzl/7vxVyqRnB5FKhH3qVJURWtmElRc4XXqAbMlVVmS8FV0wl4ZJrG0mu9XpW8P1jUlLdhZRw2IWrK+6JCm7/B2Op1LrpJWdSUd60pECLRpee8XnnYXQEAALYFRVoAAADYyJBhWAoWz89aNNrTWJ64NlIxQ2lzhVctzfO6cv1zc+oWbGauaFpmXtpXnic2KtNYsYyCdcrF1aOgsRTT6IC38LNLrXNuyZFIwXcAAACArUeRFgAAALYyRrNFoz4NjRZf+r9cRIxWzlBayRhVtmi4qicUWfW3ldtZdkoDY/QVpzvITWWwNo8ZUSTkWfc6S1IkGlXJG7wBAABgyzRks9ksYQAAAFi/6elpHT58WFNTU+ro6CAgm2Ip2BCUVaogW+4dwQYFrfW/Hlv4PZUoSIN+AwAAbB1G0gIAAKAqGGZInlI3GsP2sSxZ8ihkUqAFAADYThRpAQAAUB08pqyQR5ZpKko0dkBUpmnJE7JkcscwAACAbdVICAAAAFAtPGZEWZM47FC0ZUayItwAAADbjyItAAAAbGBolFsjAAAAAJKY7gAAAAAAAAAAbEWRFgAAAAAAAABsRJEWAAAAAAAAAGxEkRYAAAAAAAAAbESRFgAAAAAAAABsRJEWAAAAAAAAAGxEkRYAAAAAAAAAbESRFgAAAAAAAABsRJEWAAAAAAAAAGxEkRYAAAAAAAAAbESRFgAAAAAAAABs1EgIAAAANiccDhMEAPQXAADglVGkBQAA2KS+vj6CAAAAAOCVNWSz2SxhAAAAWL90Os2oOJQUDofV19enoaEh+f1+AoJV/H6/XC4XgQAAAAUYSQsAALBBLpdLHR0dBAJl+f1+9hEAAACsGzcOAwAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAG1GkBQAAAAAAAAAbUaQFAAAAAAAAABtRpAUAAAAAAAAAGzUSgvq0sLCgVCqV//9MJqNkMqlMJvPS97pcLrlcrvz/NzU1yel0ElRsuXQ6rXA4TCAAbJjf7y84VgFbKZVKaWFhIf//yWRSi4uLkpQ/bhUfv5xOp5qamsrmUwD5EwDyJ9RLvlROvedLDdlsNsuuUxvJWDweVzqd1vz8vBKJhP7xj38onU4rlUopmUwqnU6vq1FsVmNjo5qbm9XY2KjW1tZ842pubta//du/qaWlpeBBYRcvMz09rcOHDxMIABs2NTWljo4OAoF151GpVErxeFzJZFLJZFJ///vflU6nlUgklEql8v9ezw/aG+FwONTc3Kzm5ma5XK78v1977TX5fD65XC653e6CExaA/AkA+RPIl2ovX2IkbRVZWFhQLBZTLBZTJBLR3NxcvjCbSqVKNoKmpqb8o7W1VS6XK19Mfe211/I7cGNjo5qamuRwOF66HsWF3nQ6rXQ6rX/+85/5YnCu4YbD4bIjdJ1Op9xut1paWuR2u7V79275fD653W4KuCgwNDQkv99PIAC8VDgcVl9fH4FASalUSuFwWOFwWI8fP1Yikcj/uF2cq6w8GWhpaZHf71dzc7P+9V//ddVVROu5qiiXL+VyulQqpb/97W/KZDL5k5pYLKbp6emSP6g3NzerqalJbrdbe/fuld/vl9/vV0tLC18syJ8AkD+BfKkG8iWKtBUok8koHo9rbm5Ojx8/ViwWUzgc1vz8fL5B5Aqtzc3NCgQCam1t1be+9S21tbXlC7LNzc0V1cgXFhYUj8eVSCSUTCb1pz/9SfPz84rH45qYmCho7K2trfL5fPL5fPJ6vfl/Nzayy9Yjv9/PL7oAgA2Zn5/X7OysHj9+nD/RSCQSBUl8W1ubAoGAfvazn6mlpSX/t+bmZlt/MF5cXFQsFlMqlcr/989//rPC4bAePHigsbGx/GtdLpcCgYD8fn/+ZMTtdpMzgfwJAEC+VGX5EtlbBUgmk5qdndXXX3+tyclJxWKxgnk72tra5Ha71d3drb1798rn86mtrW1do14rxcrRvKVkMpn8KOFYLKZvvvlGsVhMd+/eLShM+3w+dXR06M0331R7e3vZ5QEAgPqRyWTyJxZ/+MMf9ODBg4K599va2nTw4MGCERWV9GN2MYfDkR8B2dnZuer53JyjuREu4XBY165dy+dMTqdT7e3teuutt9Te3q5AIFBVeSMAACBfqsd8iSKtDdLptCYnJ/XHP/5Rs7OzCofD+Z2kvb1dvb292r17d75R1MOl/7kCrM/nK/j74uKi5ubmFA6HFYlEFA6Hdf36dV27dk2S5Ha7FQgE9NZbb6mzs7OiOxAAALB1EomExsfH9fXXX2tmZiY/6qOlpUWdnZ36wQ9+oEAgoEAgUHO5lMvlUkdHR8EoyVzONDs7q2+++UZ3797V5ORk/vUHDx7U//t//0/Hjx/nR24AAMiXyJcqMF+iSLtDksmkxsfHdevWLc3MzGhxcVFOp1M+n0/9/f364Q9/qEAgwE0hijgcjnyHkbOwsKBwOKzJyUl9/fXXmpiY0M2bN9XY2KhAIKB3331X3d3dnIAAAFBjUqmUPvvss3w+lclk1NTUJL/fr1AopPb29qq72mg7c6Z4PK4HDx7oD3/4g2ZmZjQxMaG+vj61t7frvffeU09PDz9wAwBAvkS+VCH5UkM2m82yC2+PeDyusbExffHFF5qZmZH04teKY8eO6ciRI1x6tkWmp6d169YtjY+P538ZCgQCOnbsmLq7u+V2uwlSlX6vhw8f5i6jAOg36tji4qLGx8f1u9/9TpOTk8pkMmppadHx48f13nvvqb29nflX1+nBgwe6deuWxsbGND8/L0nq6OjQT3/6U3V3d3PjVvpBAPQb9BvkS+RLNuZLFGm3oWGMjY3p888/z8+n2traqu7ubh07dow5wXbA7OxsvkHF43FJ0sGDB/Wzn/1MJ06coGMiWQBAv4EqkEgk9Otf/1o3b95UIpFQU1OTjh8/rnfffVednZ0cz1/RzMyMPv/8c42NjSmZTKq5uVk9PT36z//8T65Goh8EQL8B8iXYkC/9CyHfuoZhmqZef/11nTx5UuFwWKdPn9b9+/f17NkzXblyRQcPHqRAuwMCgYAGBwf17NkzPXr0SP39/UokEvrggw/0ve99T319fQV3KwQAAJVjbm5OhmFo165dunz5snw+nz799FM9f/5cIyMjOnLkCCccW+DgwYMaHh7W8+fPdfv2bbW3t+v69evatWuXjh49qnA4TJAAACBfIl/awXyJIu0ryhX/du3apYGBAbW3t8uyLD1//lxDQ0MMKbeZ3+/X4OCgnj9/rjt37hQ0qPfffz8/0hYAANhrfn5ehmFo3759Gh8fV29vr549e6Yvv/xSPT09XIq/TRwOh44fP67bt2/r6dOnOn36tKanp7Vv3z4ZhkGuBAAA+RL50g7lSxRpX6FhnDlzRnv27NHNmzfV29urhw8f6ssvv1R3dzeF2QrU1dWl27dv6+HDhzp16pTGxsa0Z88enT9/XqlUigABAGCDVCqlixcv6o033tD4+LhOnTqlp0+famRkRG1tbQRoB7W1tWl4eFhPnjzR2bNnNT4+ns+VkskkAQIAgHyJfGkb8yWKtBuUyWR0+fJlvfHGG7p+/bp6e3v15MkTjYyMFNwdDpXL7/frxo0bevbsmX7xi1/o2rVr2rVrly5duqRMJkOAAADYIWNjY9q3b58uXbqkzs5OPXr0SDdu3OBkw2atra0aGhrSkydPdPz4cV29elVvvPGGxsbGCA4AAORL2KZ8iSLtBiQSCe3fv18XLlyQz+fT/fv3NTw8LLfbTXCqtEENDw/r0aNH8vv9unjxovbv388cbAAA7IC+vj4ZhqGFhQX9/ve/1+3bt+Xz+QhMBXG73bIsS5ZlSZIMw9CZM2f4URsAAPIlbEO+RJF2nSYmJrRv3z4lEgl9+umnun//vtrb2wlMDfD5fPrqq69kWZYSiYQOHDjASBEAALZJOp3W0aNHde3aNXV3d+uvf/2rTpw4QWAqWHd3t54/f67u7m5dv35dR48eZaooAADIl7DF+RJF2nW4fPmy3nnnHbndbj18+FA9PT0EpUYb1KNHj9Te3i7DMGSaJkEBAGALpVIpHThwQBMTE/rVr36lTz/9lBtcVAmn0ynLsjQ4OKjJyUnt37+feWoBACBfwhbmSxRpX8I0TV24cEE9PT26d++eWltbCUoNa2lp0b1799Tb26uBgQEKtQAAbJGFhQW9/fbbisViGhkZ0UcffcSNVqtQf3+/RkZGlEgk9PbbbyudThMUAADIl7AF+RJF2jVcvXpVAwMD6u3t1ccffyyHw7GBd0dlehvU0FD0CFoVsnWWgg0l1q/sIygr/56grBr+3hsbG/Xxxx/r1KlTGhgY0KVLl2gMdS5qetXgNRWlndeAjW7f1sbDCtZnv1ruGOk1o3QwdeTkyZOanZ3VyMiITp06tUNtmDa5HXp7ezUyMqK5uTkZhkFAUIF52styoFpq4/blNqv70K3vp8ltQL5UieeR5EnblS9RpC3jwYMHunjxonp6ejQyMrKxXy6iprwNXg1ESx7J2KErMllZbWRkRCdOnNDAwIBmZmYIaf1m/jIGojJMUx7aec218x1d22CDguwUyzwyTUPRAUOcy9SHTz75ROPj4wqFQjaecNAmt/rE47/+6780OTmp69evExBUVp4G+lByG5AvoQrzJYq0ZVy4cEEul0s3btzY4NDyqExjYOmXXE9IkWxW2ZWPSKhCEghDo8Xrlh2V8dLn6suNGzfU0tKiDz74gDsZ12sCag4oKkOFP3zRzrHhPUlLN/v0KBQh3ku7oCFDUQ2YnOHVunQ6rfPnz8vn86m/v582WUNM01QgENDFixeZ9gAVlKe9LBeivdOHktuAfGnz55O09+3MlyjSljA+Pq7p6WkNDw/L5XJt8N0RRaPLByyrxC+6HlMRduiq4XQ6NTw8rHg8rvHxcQJSrwmoYRS1Wdo5NssrL0N98kmfYUiyLEad17hPPvlEqVRKw8PDG5w6ijZZDUZGRrSwsKDf/va3BAMVkqeBPpTcBuRLqM58iSJtCb/5zW/U1tam48eP7+jnRk3vqrk+Vs9pU3h574t5gIou97WC2z5PZvH6rl78+tZ1fdstlZ7/01vichJLwVVzIm1+7pSuri653W79+te/pnHUQDL/Yp8s3p9K7BvLCZbH46Wd13g7f7V9aa04FU/JsHJdy63nWtM4rPHcuveH9a578evXMyfVer+/5VM7j+fFSTZq1ueffy6/36+Ojo4qbZObbJdV2CY3w+/3KxAI6He/+x07O6o7T3tpm32xTqWb8vL7i+fH3WRfUDb/q5c+dN2xq6w+lNwGlZ0vbWc3sp72up39qJ196NblSxRpi2QyGU1PT6u7u3uTd9DLXW4T1YB3vTvE8iTjJSa3jA6Uu2FR+XmArGCZA/8WzpNpBVev71rzEpV+biPbbSlYcv7PpThv93xI3d3dmpmZ0eLiIo2kVk4DgsX708qiXy73tyR5ZAQ9tPM6aOc7FaftWofN7A9bsu4ej7yb/P48QUOefFtDLUqlUpqdndV7771Hm6yCNrlZ7777rubm5pRMJtnpUQF52na1WY+ChmfFZxe9dPlvHiOYv8rqVfoCu/KeSuhDNxu7SuhDyW1QLfmSPe11Z/rRSp0zez35UiPNoVA8Htfi4qL27t276WUYoxGFvLkOfama/2KfDCkSKbw8Omoay6/1KBSJyPTkn5DXO6BodECGGVTE9BQkKpZlaLTokuqo6V3eGYufi8r0ejUQtWSaURnmqyQzqz871wgs01TUKL78u9y6bnC7S8ROVlANQavoc5c+aymRW/25m7F3715lMhndvHlTbW1tNJQdEg6Ht+swUrRP5hL/le0jqkiEdl5P7Xy74jSaza7Yx4rX9dUvh9v8/rDe79hY3obV38mqKT/W/f0ViUQUlbGlczlvX/+BjXj+/LkymYxaW1tpk9XUJjco9/1OTEyQJ9V0/mTXsXXjedrKZRXkaFrdNjfSZj2mKWMgKMuyZMko7DuWl2Eut+2t7AuqJ6/ZGpuLXYX1oduQ29B/kC/Z0R9vZR8qaQf6UTv70PXlS8lkUs3NzaVflEWBqamprKTs1NTU1ixw1MhKKnp4sqFI7gWRbMiz9HdjdI33e0LZpbeMZg2Ve31uWSuXX7Cwpffml1XmeRnZ0TWfL7H8SCjrWfXe9azrerd75WeUehSv88u2ZXP7BQ97HlvWHtfaJ8u2tXXsQ7TzKmrnG31PuddvNk7r2Ya11rH4uc3sDxtd9/LtyFO8gA19f1vfV9Nf11M/vpNtciPtstrbJO2OdmenrczTXixr7X1/42121Cixjrn1y//xVfqCMnnTtuY2W9WH2pHbVFofuvW5zU7WQXjQb293H7p9/ehm+9DKqjcykrZIrpq9ZZdrGaNa+UPd0q95UQ0YpoIRU578DYjK3JXUMGTIkhWNKiKt+CWu1Otzy1q6BGOg7E+TxcvaqBITxHu88koqPT3PWuu6zu3O/7ppj1QqJWlpsmdGiOzsL7l9fX3bsOQS+13ZtkY7r5d2vjNx2mqvsj9sZt2jMr3LlxV5QrJW/mpdYd/f0NCQ/H4/HWkFjAz5j//4jx28DJ42aUebJE+ql/xpp2xlnvay0VQbb7OGYSzNh5sfSRmVaVrLq2lsQV9QJm+qiz70VWJX23kN+RP5kl398Vb3odr2ftTOPnR9+ZLL5Sr7Goq0Rdra2uRwOPT48WP19PRs/S4+OiqjIbi5QlDdetFgPaHIqsvBgzswEfTjx4/V2NioEydOcKdF0M5rtJ3j5V7MEWhotODyv8r7/qr6xgs1xOfzqbGxUfPz8wSjhttk7vvt6uoqf/keUDOJXq5gbGk0asrUqKxoZRcGUP15DfkT+RL9aPXLfb9r5UrcOKxIY2OjOjo6NDY2pkwmsw11iIgKp0/yyuNZ6uxLzi2+fNfSFxOYryW3LI9Ckayy2XKPSpibYyPb/WI0nlk8Z9OqeG6PsbExHTx4kAJt7aRkq/e7VW0tt49GFInSzuuhnVeeEvte2X1rB/YHK5ifgN8YLV7eJr+/3HPr2vdRjZqamhQIBPT555/XSbus8ja5Sbdu3ZLP56NAi+rI0175OGrIDC3Ni2uNRhUdtRSV5AmZK15TTflaneY229WHktugrvKlzbbX+uxH15MvUaQt4ec//7ni8bjGx8c3lVQEG8rdSc5S0DuwdLlFvtNecXe7oFfmyoNW1JQ3uPruduV5ZJqGloaDFy1r5TIr4pe/jWz3i0KPufKFuRsPbTQR2KCJiQnFYjH98pe/pHHUUvofXNkOLAVXtTWPvF5JiioaoZ3XejuvzGQnqgHTKr1v7fT+sOL7kjGqUaPcOm/w+4tEl57zehlxXsPee+89hcNhTU9P10G7rPI2uQnhcFizs7P66U9/ys6OCsnTNparbKbNeoJLN4SKWiGFrKgkj4ygp0rztTrMbbazDyW3QV3lS5tvr/XWj643X6JIW8Lx48fl9/t1/vx5LSwsbDKxaFBDQ/Ejt/MU3jXSY1oK5Q5a3hWvzxd6iubHWUv+IFO0rOJlVkJzXvd2vyj0RAe8ha8zjDK/muSGyRfHYWMNeHFxURcuXFBLS4u6urpoHDXDkGEsFVoL94vCX8tz8+FEylQAaee10c4rN9nR0nxoK7erVCy2fX+IyjRWLKNgnXLx3sz3J0WXb81tcH1oTTtx4oRcLpfOnz9fH+2yitvkZpw/f14Oh0O9vb3s7KioPG39H7eJNusJyvBIio5pLCrJYyjoqd58rb5ym+3tQ8ltUHf50mbba531o+vNlyjSlnHjxg0lk0l9+OGHG5z2wNBoNlvi17gXO1k2G1FhLcYjM5JVJLS6QOMJRZQtmB9nPW1kjSHfxmgFDQdf/3Z7zMiq13lCEWVHjTXiEFHoFX++/PDDDxWLxfTpp58y1UGtpf+jxe20xKToy0lYdNUQDdp5LbXzSk12ird7rVhs7/6Qu+TvJd/0Jr6/SDQq5vCrfU1NTRoaGtLs7KxM06yLdlmtbXKjTNPU5OSkBgcH1dTUxM6OCsnTNvd5G2uzL4p4UvmroaonX6un3GZ7+1ByG9RjvrS59lo//eiG8qUsygqFQllJ2XPnzhGMOtTf35+VlP3oo48Ihk2mpqaykrJTU1NbtMTRrCFlJSM7ut53GBt7PSrdxvcBbPN3YYxWeL+BrdLR0ZGVlL19+zbBqAFffvlltrGxMXvw4MHs//3f/xGQms6fKvcYTZ5GblMPuQ39BvkS6idfYiTtS6rd586d09WrV3XhwgUCUkcuXbqky5cv6+zZswqFQgSknkdymCF5yt34CsDmWZYseRQyGWpSL27fvi2/3y/DMDQxMUFAqtjk5KQMw1BbW5tu376txsZGggLyNIDcBuRLeMV8iSLtS1y5ckX9/f26fPmyjh49qnQ6TVBq2MLCggzD0MWLF3Xu3DlduXKFE4965zFlhTyyTLOO5g0DtltUpmnJE7JkcleNuuFyuXTnzh21trbqnXfe0c2bNwlKFRobG9M777yT/z6Z5gDkaQC5DciXsDX5EkXadRgcHNTw8LDu3r2rN954Qw8ePCAoNWh2dlb79+/X2NiYhoaGKNBiRf6/8TljAazZqpbmKuYspu60tLTo3r17crvdOnnypM6fP7/Buf9hp/Pnz+v9999XS0uLvvrqK7W1ca8+XgAAIABJREFUtREUkKcB5DYgX8IW5UsUadfp9OnTunPnjiTpwIEDOnPmDKNqa8TCwoL6+vp04MABpdNp3blzR2fPniUwNWnphl91dZMGAKgwra2tunfvnjo7O3X16lUdOHBAsViMwFSweDyuQ4cO6erVqzp48KDu3bun1tZWAgPyNAAgX8IW5ksUaTfgyJEjevr0qY4fP67r169rz549unnzJr9oVLHPPvtMe/bs0bVr19TT06MnT56oq6uLwACcAALYRs3Nzbp3755CoVD+SpZr166RU1Wg69ev68CBA5qZmVF/f7/u3bunlpYWAgOQ2wAgX8IW50sUaTfI6XTq9u3bsixLDodDJ0+e1KFDhzQ5OUlwqsjMzIwOHDig999/X42NjbIsS7///e+ZVw0AgB1kmqampqbU1NSUv6pldnaWwFSAcDisQ4cO6cyZM3I6nbp3754GBweZCgoAAPIlbFO+RJF2k7q7u/X06VMNDg4qFovpxz/+sfbt26exsTF+1ahgExMT2r9/vw4dOqRYLKZQKKQnT56ou7ub4AAAYIOOjg49ffpUv/rVrxSLxbR//379+Mc/VjgcJjg2mJub09tvv619+/ZpdnZWH330kZ48eaLOzk6CAwAA+RK2MV+iSPsKHA6H+vv79fz5cw0NDSmVSskwDO3atUuXL1/W/Pw8QaoAiURC165d065du3T06FElEgkNDg7q6dOnMk1TTqeTIAEAYHNO9dFHH+np06cKhUKamZnRvn37dPjwYU1OTvID+A6Ynp7On2xMT0+rv78/fzJIrgQAAPkStj9foki7BVwul86ePatnz55pZGRETqdTFy5c0Pe//329/fbb+uyzz7S4uEigdtDi4qLGxsb0zjvvaNeuXerr65PD4dDIyIiePn2q/v5+NTc3EygAACpIS0uLTNPU8+fPde7cOYXD4fzVSlevXuWmrVtsYWFB169fz5/gzczM6Be/+IWePXumwcFBbg4GAAD5EvnSDuZLFGm3kMPh0KlTp/T06VM9fPhQJ06c0OzsrN5//319+9vf1ocffqi7d+9SsN0mmUxG09PT6uvr0+uvvy7DMDQ5Oanu7m7dv39fT58+1alTp+RyuQgWAAAVrLm5WVeuXNH//u//anh4WJlMRufPn9e3v/1tnTx5knxqC/KlDz/8UK+//rrOnDmjhYWFgnhzYzAAAMiXyJd2Pl9i5v9tEggE9PHHH+dHdN66dUuffPKJfvvb38rpdKqjo0Pvvvuuurq6uFnVK0ilUpqcnNQXX3yhiYmJ/C9GR44c0b//+7/ryJEjxBcAgCrlcDh0+vRpnT59Wnfv3tXnn3+u8fFx3bx5U01NTTp+/LjeffdddXZ2ckOrl5ientatW7c0Pj6uRCIhp9Oprq4uvffee+rq6iJ+AACQL5Ev2Zwv8e3sQGM5ceKETpw4oXQ6rfHxcX3xxReanp7WxMSEJMntduvIkSN688035ff71dbWRsMpY35+XnNzc/r66681OTmpubk5LS4u5gvfP/nJT3T8+HGmMgAAoMYcOXJER44c0fDwsO7evatbt27ps88+K/gB/K233lJnZ6fcbnfdxysej2tmZkZ/+MMfNDMzo1QqJYfDoa6uLh07dowfsgEAIF8iX6qwfIlK4A5yuVzq7e1Vb2+vFhcXNTMzo//+7//WzMyMrl27ln+d0+lUZ2en3nzzTfl8Pvn9/rq87CyVSikcDiscDuubb77R5OSkUqlU/vlAIKDe3l795Cc/0cGDB5nGAACAOuB0OtXd3a3u7m5lMpn8CcjExET+B/CWlhYFAgH96Ec/yp+E1PoP4PPz85qYmND//M//aHZ2VvF4XJlMRk1NTers7NSxY8fU1dXFTcAAACBfIl+q0HyJIq1NHA6HOjs71dnZKUlKp9OamZnRH//4R83NzWlmZkbj4+P517tcrnzBdvfu3fL5fPL5fDUxAiKVSikWi2lubk6RSETxeFyzs7MFBdmV2/+jH/1I7e3tzJcGAECda2xsVFdXl7q6uiRJsVhM4+Pj+cQ7l0s5nU75/X75/X595zvfkdvtzl+9VG0SiYTC4bDm5ub097//XbFYrCBvampqkt/v13vvvacjR47I7/fL4XCwswAAQL5EvlTh+RJF2grhcrkKGo0khcNhzc7OKhKJ5HeumZmZgvc1NTXJ7XbL6XTK7XZr9+7d+X+3tbXZPro0k8kolUopkUhofn5eqVRKf/7znxWLxZROpxWPx5VMJktuk8/nk9frVXt7u/x+P1NAAACANbndbvX39+f/P5c7RSIRzc7O6re//W3BDTScTqd8Pl/+x+Ddu3fn8xCXy2Xb9EnxeFyxWEzJZFJ//vOfNTc3p1Qqpbm5OS0sLORf53A41NbWpq6uLv3gBz9Qe3u7AoEAOwIAACBfqsJ8iapXBcv9grFSboebm5vTn/70J83NzWl+fl7pdFp3794tuZyWlha5XC41NTWpublZDodDDodD3/nOdyQp/7dc43vZ6NxUKpXf4RcWFvK/RPz9739XOp1WKpVSMpnUwsKCEomEMpnMqmU4HA41NzfL5/Opp6dHu3fvVmtrq/x+P/PJYpVwOEwQANBfYMNyVx6tlLt65/HjxwqHw4rFYorFYiXzqKamJrlcLrW0tOT//d3vfjefX+V+QF6ZS5WyMh9KJpNaXFzUP/7xDy0sLORzqWQyqVQqVXAl0cqTo6ampvwP116vV4FAQG1tbYySBf0hAPoLkC/VSL7UkM1ms+yStSGZTCoej+dHqKbTaf3tb38rKJwuLi4qk8komUyWLJ6+isbGRjU3N+cLwk6nU83Nzfrud7+bb7DNzc1yu93cqALrMj09rcOHDxMIABs2NTWljo4OAoF1SyQSisfjSqVSisfj+R+fE4mEUqlU/t9bmT/l8iWXy1VwcvPaa6+pra1NLS0tamlpUWtrK1cUgfwJAPkTyJdqPF+iSFvH0um00um0pKVpCRKJxLrel2sYxQ0G2I59lF91AWyG3+/nhpLYNrkfvleesORORsLhsPr6+jQ0NFRwRdTKq5UcDoeampoovIL8CQD5E+oyXyqn3vMlirQAAADAFsmNYmQ0EgAAADbiXwgBAAAAAAAAANiHIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANiIIi0AAAAAAAAA2IgiLQAAAAAAAADYiCItAAAAAAAAANioMZ1OKxwOEwkAG+b3++VyuQgEttTCwoJSqVT+/zOZjJLJpDKZzEvf63K5CvbJpqYmOZ1OgloHyGdQKXL7IfsjyNcAAMBGNExNTWUPHz5MJABs2NTUlDo6OggE1pROpxWPx5VOpzU/P69EIqF//OMfSqfTSqVSSiaTSqfTSiaTWlxc3JZ1aGxsVHNzsxobG9Xa2iqn06mmpiY1Nzfr3/7t39TS0lLwoLBbfaanp0U+AwDkawAAVKvG3D+Ghobk9/uJCICXCofD6uvrIxDIW1hYUCwWUywWUyQS0dzcXL4wm0qlSo6CbWpqyj9aW1vlcrnyxdTXXntNTU1NSweqxkY1NTXJ4XC8dD2KC73pdFrpdFr//Oc/88XgVCqleDyucDhcdoSu0+mU2+1WS0uL3G63du/eLZ/PJ7fbTQG3wpHPAAD5GgAA1ShfpPX7/fzCCgBYUyaTUTwe19zcnB4/fqxYLKZwOKz5+fl8sTNXaG1ublYgEFBra6u+9a1vqa2tLV+QbW5urphtSqVSWlhYUDweVyKRUDKZ1J/+9CfNz88rHo9rYmKioJDb2toqn88nn88nr9eb/3djYyM7SAUgnwEAAABQjTijBACUlUwmNTs7q6+//lqTk5OKxWJaWFjIP9/W1ia3263u7m7t3btXPp9PbW1t6xr1WilWjuYtJZPJ5EcJx2IxffPNN4rFYrp7925BYdrn86mjo0Nvvvmm2tvbyy4PAAAAAIBiFGkBAHnpdFqTk5P64x//qNnZWYXD4Xwhsr29Xb29vdq9e7f8fr/8fn9dXPqfK8D6fL6Cvy8uLmpubk7hcFiRSEThcFjXr1/XtWvXJElut1uBQEBvvfWWOjs7K2r0MAAAAACgws49CQEA1LdkMqnx8XHdunVLMzMzWlxclNPplM/nU39/v374wx8qEAjk54jFEofDoUAgoEAgkP/bwsKCwuGwJicn9fXXX2tiYkI3b95UY2OjAoGA3n33XXV3dzPKFgAAAABQgCItANSheDyusbExffHFF5qZmZEktbS0qKenR8eOHdORI0eqasqCSuF0OnXw4EEdPHgw/7fp6WndunVL4+PjOn/+vM6fP69AIKBjx46pu7tbbrebwAEAAABAnfsXQgAA9WFxcVE3b97U0aNHtWfPHl24cEGJRELnzp3TV199pWfPnunjjz/W8ePHKdBuoY6ODg0PD+uvf/2rHj58qP7+fqXTaV28eFF79uzRoUOH9MknnxTcnAwAAAAAUF8o0gJAjUskEjJNU6+//rpOnjypcDis06dP6/79+3r27JmuXLmigwcPUpjdAYFAQIODg3r27JkePXqk/v5+JRIJffDBB/re976nvr4+JRIJAgUAAAAAdYYiLQDUqFzxb9euXRoYGFB7e7ssy9Lz5881NDSk9vZ2NTYy641d/H6/BgcH9fz5c925c0ft7e26fv26du3apffff1/xeJwgAQAAAECdoEgLADVmfn5eZ86c0Z49e3Tz5k319vbq4cOH+vLLL9Xd3U1htgJ1dXXp9u3bevjwoU6dOqWxsTHt2bNH58+fVyqVIkAAAAAAUOMo0gJAjchkMrp8+bLeeOMNXb9+Xb29vXry5IlGRkYUCAQIUBXw+/26ceOGnj17pl/84he6du2adu3apUuXLjFnLQAAAADUMIq0AFADEomE9u/frwsXLsjn8+n+/fsaHh6W2+0mOFWotbVVw8PDevTokfx+vy5evKj9+/crHA4THAAAAACoQRRpAaDKTUxMaN++fUokEvr00091//59tbe3E5ga4PP59NVXX8myLCUSCR04cEBjY2MEBgAAAABqDEVaAKhily9f1jvvvCO3262HDx+qp6eHoNSg7u5uPXr0SO3t7TIMQ6ZpEhQAAAAAqCEUaQGgSpmmqQsXLqinp0f37t1Ta2srQalhLS0tunfvnnp7ezUwMEChFgAAAABqCEXaKhE1vWrwmoqufkamt0ENDUWPoFVBa28p2NCghoagtnatlrbda0bZQVB3rl69qoGBAfX29urjjz+Ww+Gws4eq435oZzU2Nurjjz/WqVOnNDAwoEuXLtEYKu7YnNvXSj28qu5D1kba0Wba3Ha0040uc2vXwQqu/P6DsmqkL+LYQo4LAAC2HkXa6jgLlDEQlWGa8hT93dvg1UC05FlBHZwAeGSahqIDhshhUU8ePHigixcvqqenRyMjI2psbLS1f6rvfsgeIyMjOnHihAYGBjQzM0NAKunYvPabNOBtULCuG0X9FCmtYL1/1xxbyHEBAMBGUKSthiTfHFBUhgyj8ETPNAaWRu94Qopks8qufERCGzhprGKGIUNRDZicBaF+XLhwQS6XSzdu3LC3QEs/ZKsbN26opaVFH3zwgTKZDAGpiGNzwQFKo0VtYtTIvbfUlTE1d4Be3v5RGdv6nordQ2RZkuRRKJLbB2phuzi2kOMCAIDtQpG2WpJ8wyhK7COKRpeTf6vEKB6PqUhdnAwsnyBbFqP1UBfGx8c1PT2t4eFhuVwum9eGfshOTqdTw8PDisfjGh8fJyAVcWx+yRFrdLk9RKOKEMQ64ZWXX6s4tpDjAgCAdaBIa8OJ3YvL/Irn2ipx6dZyYubxeLf1s1/MeVW0Tmtcpxc1vavmCdvQ3FlRU97ce1fO6WcFNzT/mNfjeXHCDNS43/zmN2pra9Px48drrg+sqH5oE32RHbq6uuR2u/XrX/+axlG1x+biz185j2nhZ69vfy9cVvF7yu3Cq5f98vlz1152qWkNLAXz/188d2/5+VpLbffOX+6+3rgWb8PK7Sy3zmtNAbHGc+vuoza6T6w1r3LxupSaN9aOuZfJcQEAQPWjSGvnKWGweK6tlScvufNAS5JHRrB4HEHuEsul+e02erJS/NnRAa+CVlSmt2idrGCJpHT5ZgYlJgqLDpS7wVmJ5NW7dCmbMZpVNrI0UsIKlklW15h/zBM0lt5LBosal8lkND09re7ubpunOajtfkib7Ivs0t3drZmZGS0uLtJIbD82v2zZy8tZYwRu6XlMN7e/W8HV71m9/HLLjmrAW35fX9+yt+L7KL3dpb6XndtHdmbbX7YOm+mjtmTdPR55899BqXljX3XuZXJcclwAAOpTIyGw7TRQlmVoNH+5Vu5kw5JpRmWYHklRRda4HtIYjSiUTziXRhC8yOpCikTK3czEkhUJKZJdej5qeuUdiMoKepfn/ir8e9QaVXTFjVGiprH8mR6FIhGZnqKkNDogwwwqYnrKZc/LSWrh+6OmdzmhXxmXXMLs1UB0ZWxKiEQUlcEcmDsoHA4ThB30l7/8RYuLi9q7d2/FrFOt9UNb0hftsL179yqTyejmzZtqa2ujP7L52LxyWQXt4UWjUXbUWOfnv8r+vnpZuWKcZZqKGrk2VmbZisr0mutez1LLLrHxGs2OLsd09XZKxZd1vxhBaIy+mNP3RXxtKdG+dNtHs9kV+07xdr76peub76PW+70Zy9tQps9cOQVBqX5++bVr7wv1dWwhxwUAAOsyNTWVlZSdmprKYieMZg0pKylrjBY/ZWQlZeUJZSMFrzWyoy9d7PJ7Cx6ebChS6rM3+veVnx/Jhjxl1r/kNhQtJ7+exduUW27xOhQto2C5a60ntlOu3+Bhz6Ni++uq74detS+iPVZv+9jKY/OLZZV7rN5v1/j8De/v5dpSNpuNhLKegnV/ybLLbtt6ll0uTmvFr/i5F+unTR/nN5ojvGy917PtG93+zcRkI33URte9/P7pKV5A/v2lHluUl5HjbjqWnOcBAFBdGElr3/iz1XeENgwZsmQt31BkQ7+WG6NaOehhaWREVAOGqeCq0QblbmKxnptb5G7mUOaO1mtug6VgsNwIiNxyly5tGyg7fGQTscG2GRoakt/vJxA75C9/+Yt++tOfKplMVmi3Vu39UHX2RalUSpI0MjJS9yNp+/r6KujYXGqk6PKouWCDvKFIidF4pfbpze7vJdqSxyuvtOJy7Zcsu6z1LHsreGRaIVnLl44XjKhcczTldtqpbX9Z/7eZPmoz6x6VmZv2whOStXKfzY8a5dhCjgsAALYCRdoaZYyOymgIbq7gu91edqkYqobf71dHRweB2CGZTEY///nP9fjxY/X09NAP0Q9Jkh4/fqzGxkadOHFCDoeDhlLRPDIjo4o2BGUVXWaNciEzFcmaS/9eWRSMDsgb9K4xdQS2wov5XQ2NFhQfozLNpe/Cs+oHh+2dL5hjCwAAqFXcOMy+tHf1nVqX7xb94oYMXnk8khRRZKPDM6IRRbZlvXPrVOZOs6u2oSCt1ujyyVR0wFt0s4bccj0KRbLKZss9RlffbCW3rSU/E6gdjY2N6ujo0NjYmDKZTOWvcNX1Q6/YF9lkbGxMBw8epEBbDcfmgvYRXWf7eJX9/RWXXWmM0cL2Z1mqrdspldinVvWjO9hHWcH8zb+M0eLlvRh1ahYXJLet7yfHBQAAtY0irZ2ngsGVowwsBZczYY8RXB4V4JHXK0lRRSOrTySDDeXunGspmLs0cMuTOo+Chmd5/b0qyEGjpryrtqHECVbJJNYj0zS0dClY0XJXLr/UyIxIdGlbvV5GJKHm/fznP1c8Htf4+Hgl9GI11g+9Yl9kg4mJCcViMf3yl7+kcVTEsfllVl46vt528Yr7+7qXXdyWV6zrtlhPkdtSsKFEO6y54myuiBfVgGmV7kd3uo9asW/JGNXqAcsvCprmypXI3WCLYws5LgAA2DCKtLYxZBhLSWhDQ4Ma8olZ4YgEY3lSrEiZMxkrmHv/ykduWUV34N2qFNa0FMqdTHhXfG4+aS6as2y9SWz+JKBoucXLL85rl2+zbRhc8ojad/z4cfn9fp0/f14LCwsVsU411Q+9Ql+00xYXF3XhwgW1tLSoq6uLxlFBx+bllrFiObmH98Wd4zfQLl55f1/Xsovbcm5dtyPGpdpXuaJiiXaYLx4aqo0jf66Ip6UpHVbGpNQ2bnsfFZVprFhGwTrlvqsXBc3ogLfws0t+L7n2sP7iMTkuOS4AAPWGIq2dp4Kj2aKRCSVuMrKc6EZXDdcxNJrNquxUbMaostmItmdKLI/MSFaR0OqFe0IRZddzI4+iJLbBayq6HJOyl+gZoyWfi0Sj2vhNT4DqdePGDSWTSX344Yc2T3tQm/2QNtkX7bQPP/xQsVhMn376KVMdVMyxeT27bkiRDbeLLdjfN7xsj0KR7dnPjdGIQp719jGlX2uMZmtrPtoV/VHBd1tmG7e3j8pNZfCSPceMrNpv1lrnej+2kOMCAICXmpqaykrKTk1NZbETRrOGlJWM7Oh632Fs7PV1GU+D6Owk+g37hUKhrKTsuXPnCEYd6u/vz0rKfvTRRwRjS/oljs3kXdhWkVDWI2U9oQix2MEcl3wNAIDqwkjaahhcYYbkqZabeuw0y5Ilj0ImQwxQX0zT1Llz53T16lVduHCBgNSRS5cu6fLlyzp79qxCoRAB4dgMVLzoqKVoqZuMgRwXAADkUaStBh5TVsgjyzQrYh7ECkr5ZZqWPCFL5PyoR1euXFF/f78uX76so0ePKp1OE5QatrCwIMMwdPHiRZ07d05XrlxRY2MjgeHYDFR8vjZqRWWMjopyIzkuAAAoj7O7qjkXjChrEoeiqMiMZEVYUM8GBwf13e9+V319fXrjjTdkWZba29sJTI2ZnZ3VyZMnFYvFNDQ0pLNnzxIUjs1AVeVrIMcFAABrYyTtjlu6GUIl3HQGQG04ffq07ty5I0k6cOCAzpw5w6jaGrGwsKC+vj4dOHBA6XRad+7coUDLsRkAAABADaJICwA14MiRI3r69KmOHz+u69eva8+ePbp586YymQzBqVKfffaZ9uzZo2vXrqmnp0dPnjxRV1cXgQE2hAI8AAAAqgNFWgCoEU6nU7dv35ZlWXI4HDp58qQOHTqkyclJglNFZmZmdODAAb3//vtqbGyUZVn6/e9/r6amJoIDAAAAADWKIi0A1Jju7m49ffpUg4ODisVi+vGPf6x9+/ZpbGyMkbUVbGJiQvv379ehQ4cUi8UUCoX05MkTdXd3ExwAAAAAqHEUaQGgBjkcDvX39+v58+caGhpSKpWSYRjatWuXLl++rPn5eYJUARKJhK5du6Zdu3bp6NGjSiQSGhwc1NOnT2WappxOJ0ECAAAAgDpAkRYAapjL5dLZs2f17NkzjYyMyOl06sKFC/r+97+vt99+W5999pkWFxcJ1A5aXFzU2NiY3nnnHe3atUt9fX1yOBwaGfn/7d1/aBv3/cfxl7/oD0EFE0yjGlWZwxyiSBpRQCEO85jNDHWJSxxiGY+5zGEJTVlCa1ghgRTLpNBABslwSjpSSKCBtJZDXOrQjDrYUIe4xCUKlRSVqMSsHlOZxvSHYYap6PuHLcU/5NiyJZ8lPR8gmtry+fy5+3z0uvfdfe6yHj16pJMnT8put9NQAAAAAFBFTDQBAFQ+s9msI0eO6MiRI5qcnNR7772n4eFh3b59WxaLRZ2dnTp06JAaGxtlNptpsCJLp9MaHx/XJ598oqtXryqVSslisai9vV1vvPGG6uvraSQAAAAAqGIUaQGgyvh8Pl25ciV3ReeNGzd09epVffDBB7JYLGpsbNShQ4fIslkiAAAgAElEQVTU2trKw6o2IJlMamRkRLdu3dLw8LBSqZQkqaWlRb/73e/U0tJC+wIAAAAAJFGkBYCqZTab1dXVpa6uLqVSKQ0NDenWrVsaGxvT8PCwJMnpdKqlpUV79+6V1+tVXV2dTCY+OvKZmppSOBzW3bt3NTIyonA4rNnZ2Vzhe//+/Wpra2MqAwAAAADAMhxpAwBktVrV3d2t7u5uzc7Oanx8XH/72980Pj6uCxcu5N5nsVjU3NysvXv3yuPxyOv1yuFwVF17JZNJhUIhhUIhffXVVxoZGVEymcx93+fzqbu7W/v371dDQ4OsVis7GQAAAABgRRRpAQCLmM1mNTc3q7m5WZKUSqU0Pj6uO3fuKBwOa3x8XENDQ7n3W63WXMF2x44d8ng88ng8FXErfzKZVCwWUzgcViQSUTwe1+Tk5KKC7MK//9e//rXq6+ursnANAAAAAFg/irQAgGeyWq1qbW1Va2tr7muhUEiTk5OKRCIKhUK54u1CNptNTqdTFotFTqdTO3bsyP27rq7O8KtL0+m0ksmkpqenNTU1pWQyqW+//VaxWEypVErxeFyJRCLv3+TxeOR2u1VfXy+v18sUEAAAAACADTEtPOAGgLVgvIDX65XX6130tWQyqXA4rHA4rG+++UbhcFhTU1NKpVK6fft23uU4HA5ZrVbZbDbZ7XaZzWaZzWY9//zzkpT7mjQ31cJqV+cmk0nNzMxIkmZmZnJXvH7//fdKpVJKJpNKJBKamZnR9PS00un0smWYzWbZ7XZ5PB51dnZqx44dqq2tldfrZT5ZxicAYDwEAAAlkSvS9vT00BoAgHWz2WxqbGxUY2Pjsu8lEgnF4/HcFaqpVEr/+Mc/FhVOp6amlE6nlUgk8hZPN/RhZzLJbrfnCsK1tbVqbm7WCy+8IKvVKofDIbvdLqfTWRHTNFQz8gwAAACAclTzn//8J8NZVgDr4fV6eSASSiKVSimVSkmam5Zgenp6TT9ntVoX7ZM2m00Wi4UGrZJ9hjyDrSAUCqmnp0fnz59fdscBQF4DAAArqclkMhmaAQAAANi4sbExNTU1aXR0NO+dBQAAAEA+/0cTAAAAAAAAAIBxKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgSjSAgAAAAAAAICBKNICAAAAAAAAgIEo0gIAAAAAAACAgUw0QXWamZlRMpnM/X86nVYikVA6nV71Z61Wq6xWa+7/bTabLBYLjYqiS6VSCoVCNASAgnm93kWfVcBm5afs59bSzy/yE8hPAMhPIC89W7XnpZpMJpNhF6qMMBaPx5VKpTQ1NaXp6Wn9+9//ViqVUjKZVCKRUCqVUiKR0OzsbEnWwWQyyW63y2Qyqba2VhaLRTabTXa7Xc8995wcDseiFwcmWM3Y2JiamppoCAAFGx0dVWNjIw0B8hPITwBAfsIGJJPJRXkpkUiQl0qEK2nLyMzMjGKxmGKxmCKRiMLhcK6jJJPJvGclbDZb7lVbWyur1ZrbuX/0ox/JZrPldnibzSaz2bzqeizteKlUSqlUSj/88EOuc2Y7cSgUWvGMicVikdPplMPhkNPp1I4dO+TxeOR0OjkAwSLnz5+X1+ulIQCsKhQKqaenh4YA+QnkJ/ITAPIT1ih74joUCumbb77JZadsPsnHbrfLarWSl4qIIu0WlE6nFY/HFQ6H9fDhQ8ViMYVCIU1NTeV2vuyOb7fb5fP5VFtbqx//+Meqq6vLdRC73b5l/qZkMqmZmRnF43FNT08rkUjom2++0dTUlOLxuIaHhxd1rNraWnk8Hnk8Hrnd7ty/TSZ22Wrk9Xo5owsAID+Rn0B+AgBsQPbk9cJi7NJip9lszl2l2tDQILvdrhdeeEG1tbXkpRIjsW0BiURCk5OTunv3rkZGRhSLxTQzM5P7fl1dnZxOp9rb27Vr1y55PB7V1dWt6SzEVrHwapSVDqyyZ2pisZi++uorxWIx3b59e9GBlcfjUWNjo/bu3av6+voVlwcAAMhP5CfyEwAA1SwcDiscDuvOnTuamJhQPB7PXalqMpnkdDpVX1+vXbt2yel05l7ldIK3kvISRVoDpFIpjYyM6M6dO5qcnFQoFMrtGPX19eru7taOHTvk9Xrl9Xqr4ta1bIfweDyLvj47O6twOKxQKKRIJKJQKKSLFy/qwoULkiSn0ymfz6eXXnpJzc3NW+psDgAAID+RnwAAwGaJx+MaGhrSl19+qVAopHg8Lmnu6tiGhga1tLRo165d8nq9ZXfyuhryEkXaTZJIJDQ0NKQbN25ofHxcs7Ozslgs8ng8OnnypH75y1/K5/Pl5uyAcgOJz+eTz+fLfW1mZkahUEgjIyO6e/euhoeHde3aNZlMJvl8Ph06dEjt7e1cJQIAAPmJ/ER+AgCgYqXTaY2Pj+vWrVsaHh5WLBaTNDdfrNfr1RtvvCGfzyePx8Pc9WWQl2oymUyGTVMa8Xhcg4ODunXrlsbHxyVJDodDzc3NOnDggFpaWqrirMVmGBsb040bNzQ0NKTp6WlJks/n04EDB9Te3i6n00kjlel2bWpq4imjABg3yE/kJ/ITGAcBMG5Ac4XZkZER3bhxQ4ODg0qlUpKkhoaGXFZaetUoyiMvUaQtstnZWQ0ODurjjz/OzW9RW1ur9vZ2HThwQD6fjwOLEpucnMwNVtlL+xsaGvSHP/xBXV1dPDyDsACAcQPkJ5CfGAcBgHGjrMTjcf3lL3/JFQuzUxgcOnRIbW1tTF9UAXnp/2ji4pienlYgENBPf/pTvfrqqwqFQjp+/Lju3bunx48f69y5c2poaOAAYxP4fD69++67evz4sR48eKCTJ09qenpahw8f1osvvqienp7c2Q8AAEB+AvkJAICtaHZ2VkNDQ9qzZ4+2b9+u999/Xx6PR1euXNF3332nzz//XMeOHaNAWyF5iSJtEQ4uDh8+rG3btqmvr0/19fUKBoN68uSJzp8/r/r6eq48MJDX69W7776rJ0+e6NNPP1V9fb0uXryobdu26be//W3uzAcAACA/gfwEAMBWMDs7q4sXL2r79u06ePCgUqmU3n33Xf3rX//SZ599pu7ububkr8C8RJF2naampnTixAnt3LlT165dU3d3t+7fv6/PPvtM7e3tHFhsQa2trbp586bu37+vI0eOaHBwUDt37tRbb72lZDJJAwEAQH4C+QkAAMPMzs7qwoUL2r59u06cOKG6ujpdv35djx490smTJ2W1WmmkCs5LFGkLlE6ndfbsWf3iF7/QxYsX1d3dra+//lqXL19e9EQ4bF1er1eXLl3S48ePdezYMV24cEHbtm3TO++8o3Q6TQMBAEB+AvkJAIBNlS3s9fT0qLa2Vp999plGR0fV2dnJiewqyUsUaQswPT2tPXv26NSpU/J4PLp37576+/t58m2Zqq2tVX9/vx48eCCv16u3335be/bsUSgUonEAACA/gfwEAEDJpVIpHTx4UH6/X+l0WsFgUKOjo2ppaaFxqiwvUaRdo+HhYe3evVvT09O6fv267t27p/r6ehqmAng8Hn3xxRcKBoOanp7Wvn37NDg4SMMAAEB+AvkJAICSmZyc1L59+zQ8PKzTp0/r8ePHTAFVxXmJIu0anD17VgcPHpTT6dT9+/fV2dlJo1Sg9vZ2PXjwQPX19fL7/QoEAjQKAADkJ5CfAAAoumvXrqmpqUnpdFr379/XmTNnZDabaZgqzksUaVcRCAR06tQpdXZ26vPPP1dtbS2NUsEcDoc+//xzdXd3q6+vjwMNAADITyA/AQBQVBcuXNDhw4dVX1+ve/fuyev10ijkJYq0z/LnP/9ZfX196u7u1pUrV4p/RiMakLumRjVLX+6Aolu6ZYLqqKlRTU2HgiV5v7FMJpOuXLmiI0eOqK+vT++88w6dAQu6rXv1PlqWfbuc+/V61z3fy61AtFz3zs3YhhvZ7lEF3DVyl28Dw5D8VKn9dTP67OaO0+QnrD07zX0erNSHgx0L+nlHcOWclVtmeR1rlHeGIkcZMyaToSrNtWvX9NZbb6m5uVmffvqpbDbbRkfaIo6r5TamltfnwGp5iSLtCiYmJvT222+rs7NTly9fLvJ8IPMdyN2Xv2AT7ZO7rDpEZbp8+bK6urrU19en8fFxGgRSNCB/X1T+QEAu+nalbmT1uWvUUbUbqZQBx6VAwK9on18cY5Cf6K+Ve1BCfsLq2cklt3uuDwcHosv26eDCnTkYXLZvRweCcznL7V4hj4EcVWljMhmqksRiMR09elQNDQ26efNmkS4GZFwtNyvlJYq0Kzh16pSsVqsuXbpU9AOMYIdbffP9xj+QUSaz5BXprbCO4ddAJqNMZkD+MlvzS5cuyeFw6PDhw0qn03SMKhcM9Ckqv/x++nY59+vlf8PT14A/u623+h0NRm3DDW53v19+RdUX4FQF+Yn+Wvo+a9w4TX7CatnJP/+FaDSypFIQUUSSXO1qd0lSRJElHTwSjS5aBhmKHFUVYzIZqmKcOHEidzVlMe/WZlwtv3E1X16iSJvH0NCQxsbG1N/fL6vVWuSk0jF/ZnFuJxrItwe5AoqU/Yd2ZbBYLOrv71c8HtfQ0BANUt2HGXNnIP3+/H2Tvl0ZcXtgfvtEo4rQHCUJUH6/8p7BB/mJ/kp+QhVlJ7dr7sT1ks+Dp1dzdajP79Lyq8KyV4S55HLTwuQoMhTKLyuNjIyov7+/+PP1M65WRF6iSJvHX//6V9XV1amtra34UWX+OnNXb2BdhZpowL1s7p/lc9M86zaLlb63ZG6hAuYiWbpOy29vybecxV9bfRkbW8eNaG1tldPp1HvvvUfnqLADh6f7TnYOn+wrz/40/2HnWuGTayN9m369+f16Y/vLwnmdFq/jerbl2top37KfPe9b4dswqI5F/87XJ/Jvm3x/90rbz+1yadltVyA/VXl/LU2f3Xh/JT+hZNnJ1SF/niu6Fl7N5Zq7d3fxVWG5K8L86nCtty9l169jeV9Y6c1rfu+z5mt9dt+s9hy1tjG5nHPUyttlreMyGapyslJXV1fxF16icbUkY+qa37+WMbW8j02X5iWKtEuk02mNjY2pvb29BPOoRRWZ2/vl73AV/LMBd43cfcs/SaJ97g0+kGjhB0puoXKvYTKhYMfydQp2FDYP0dqWsf51LIb29naNj49rdnaWTlKJhxwLpilYaX+bK8Ku1HfX27fp10b26/zrPb8+K10xvWJbrG9brq2dVlp2VH3u/GGhGNtwI9t6xe0rydXhl0tPT2yA/FTt/XUz+2yh/ZX8hNJlp3zzJy65mivbtxdcFbbavIlr7UvBjhWKB8GOZQfihbx39SKKS+519tXKzVHrz8PVlqPIUJWRlZqbm0uUlYo/rpZiTC3quFohY+rCvGSiqywWj8c1OzurXbt2lWDpEc2dxHDLXWCNNhrwzwchl3ojEQVcuW/I7e5TNNonf6BDkUDhM15GA4H5HdGvgdyt2GsJ60EFgwt/5umOHgwEFPUH1jD/5tqWsf51LI5du3YpnU7r2rVrqquro6NsklAotBlxcsk+mN2vggoEovIH5m4JiUSK37fp18b26+yZ0mX8A8oM+Nf8t61/W661nVZYtqIKuANF3IZzPzPX/sv/xpXbI9tsS6f5CKqj5hlbMhJRVP6SzdO8OeMHsv7+97+XMD9Vcn8tZZ8NLrtKb939lfxEfip6dpqf+zAYnL+iy5Xnai63XC7NfT0q+V2rzZu41n7qXjBF1cJ+E1XA7VZf9OnfUch7n/bLzPKiQ8dcwbo3GJBr1ZRTXTkqGnCrY115uNxyVDDv3XDrGpdLmKHIT6XPSnv37i3Z7yjuuFr8MVUFv38tY2r5H5tm81IsFpMyWGR0dDQjKTM6OlqCpQ9k/FJG8mcGCvq5SKbXpYykjD/fDw74M5IycvVmIqv+nqXfyy7blemNLP21vRlX3uVkl7HWn8m3PoUsYz3rWJr9gpcxr9L2xzz9qqA+td6+Tb8uXr8utP2fbvuVXsu3yTP2l3Vvy0LaaYVll2wbrnVff7p+WvM2WO9nIeN1dY7XldpfN6PPFqO/0h/pj6XKTsv380iva24ZCxY+4J/7fa7eyDOWWaT8sXBZrt5MpKD3rt5mrtxCSvx3FJSj1vOZXMxxudAxuZxzVL7vFTouly5DMV6X+7FtMcfVUo2p63n/WsbU8q85LaxDciXtEna7XZKUSCRKsPTlZy3WJnuV3gpPlff75VdQwfkJ2gs7o/aMKwBdbrmlZ9xuvZ6fWc8yNrKOxZFMJiVJly9f5kqQTb4SpKenp8S/JU+/KrhPradv06+N7tfKe5XD/Bncjhq5eyN5rtrIt73Wuy0LaacVll2ybbhWLgWCvQq6++aXveCqGlevIpFAya6UXc358+fl9XoZSDfx6pDf//73JcpPldxfN7PPbn5/JT9Van4qRnbS/PyJfeqbz0/KczXXoqvCcleErXR7ayH9NKo+d436Vlq3aFSRgt6b72+OKpC9nd7Vq+Ca74yqphy1kTxMjiI/lWdWmpqaKt0vKeq4WuwxVXJtaFyt3DE1m5esVivTHSxVV1cns9mshw8fqrOzs+iDsNs9t8MFB6IKBFw0eJl4+PChTCaTurq6ZDabaRDQtyt4WwYiA4rWdCgYHFA0YFyRsXyaLKBIJjD379ztR5qfw8n9jFseS8vr9aqxsZHts0nS6bRee+21EuUn+mu59lfyE1brwx1+l/r6ogoOBBTN93Tx+aeVR4NBBVwRRSW5/B1l0defztvr10CkWsYnxuVyz1Hkp9JnpS+//FLd3d2Mq4ypy/KS0+nkwWFLmUwmNTY2anBwUOl0uujLz57BiPYFCpjXYv4qvZWe5Jid9HnZ2Y/FT/Wb+8XzZ0qWLXst7zWK8es4ODiohoYGDjAqUp5+taxPPWMfXHffpl9v+XWcP4Nbum1ZhGVvNf4BZTIZZbJX1gTzzL+26pVQID/RX7dMfyU/oYTZSdLTJ40Hg/mfLr7gaeXB4NyC3O6NHJpn18ul3khmvg/kew3IX9B7l1UTcg+k8Q+sZY73as1RpRyTK3BcJkNVRFYaGRkpaVba3HG10HFyneNqhY+pC/MSRdo8XnvtNcXjcQ0NDZVk4J07GTZ3K0PeB8VFA3Ivmpx47mzI3L7pViC65L3zC3l69iO7g0XVF1j0bFJ15G6lWLrspU+3zPdeoxi7jsPDw4rFYvrjH/9I56jUQ42OJfvVsj719EmZ0Uix+jb9emuu48LbaNYagAvdlutpp5WeVlyqieyffVC9aP+rWfI3r1bsiUSf+WRukJ/or6Xqs+vor+QnlDI7SU+fNB5d6bNhwXLWNZ3I8n4aCPjn80ee/rAorxXy3vxjydN8SI7a/DG5AnMUGapistK1a9dK90s2dVwtdJxcx7ha4WPq0rxEkTaPtrY2NTQ06MSJE5qZmSl+nxmIqNel3IdFTc2SV56dwBUIzv/M3Nwdy967aE6O7I6vuVsmcsvuUDDbYRctO7DgyXXPfq9h5RyD1nFmZkY9PT2qq6tTW1sbnaMi+eX3L9mv5r++cNqC7JWykWckrkL7Nv3a6HVc+HuzL/fTJwAH134bTWHbstB2Cq6wX2XXtdj9Ic/f8cynmi59b83TW/XybMvo/OO+/X4/ww/5if666X22sP5KfkKps9PTk9Ba8bNh0deKcQVh7iA/T39YmtcKeW+28Odf8LVFmW21z9PqzFGlHJMrLUeRoSojKzU2NurEiRMlnJt2k8fVQsfJgt5f2WNqvrxEkXYF586dUyqV0tGjR0twKbpLgUhGmUhv/g8uV68iy26bmfuZSK8rz9sjyiydk8M/sGz+GldvZIU5bfJM+u7qVWRgKw3+xqzj66+/rqmpKV25ckUmE1M4V+yhxkBmyRm5PPvb/OAcfeblIIX2bfr1llxHV68imYgKOx4ocFuqGMt2qTdS/Nt9Fp5sWNs2zP9+/0Am774ZiUa18SuhUJ35if66sT5beH8lP2EzslOH/+kJaVe+SsH8/IlS8eZN9A+sME1BNm8t+F4h73368Bly1NrH5VKOyZWVo8hQleHSpUtKp9M6fPiwZmdnS9EJN31cLWycLOT9lT2m5s1LGazo3LlzGUmZ7u7uzP/+97/qa4ABf0ZSRq7eTKQK1/HIkSMZSZkzZ87QGQwyOjqakZQZHR0txc6T8UsZyZ8ZWPPuVtj76debtY6Fb0sYtnHntpV/oEzHDZCfQH6q5vxUpdmporMeGYoMRX4ywocffpgxmUyZlpaWzH//+18apAqPn1fKS1xJ+wx/+tOf1Nvbq6tXr+rw4cPGPAhjUwTV4Q4svgS9KHMAlec6ptNpHT16VB988IF6e3t1+vRpOgPmzq0FeuUqlwcP0K+xJXfLoIJyqTfAJSDkJ1Qa8hPKPztxDAcyFDZHV1eXzp07p5GREb3yyitKpVI0SpWMqavlJe4/WkUgEJDZbNapU6c0NTWl69evy+FwVN4fGu2Tu6Zv+dc3OAdQua3j9PS0Xn31VY2Njen06dMKBAJ0AizY1wIK9gblDgQU9Qe2fvClX2Nr7ZAKBIJy9RZ6azrIT9jqyE+omOzEMRzIUNgkb775pmw2m44ePao9e/YoGAzK6/XSMBU8pq4lL3El7RqcPHlSn376qcLhsHbv3q3BwcEK+wvzzL0hzc0FEtkqYar06zg4OKh9+/ZpYmJCwWBQZ86cYedHnmONyBbqF/RrlFXvmZsTjqML8hMqCvkJlZOdOIYDGQqbq6urS1988YXS6bT27NmjQCDAHUgVOqauNS9RpF2j1tZWff3113I4HPL7/dq3b59CoVCF7ZAZZRa+BvxVsY7hcFhNTU3y+/2y2Wy6d++e2tvb2emrYhB+xoTl9GvWEQD5CeQnVFl2IusB2Fw+n08PHjxQa2ur+vr6tG3bNk5sV9CYWmheokhbAIfDofv37+vMmTO5q0JOnDiheDxO45Shqakp9fT0aPfu3ZqcnNSZM2d07949bjEAOGgEQH4C+QkgQwHYFFarVTdv3lQwGJTJZJLf79evfvUrjY2N0ThVlpco0hbIZDLp9OnT+vrrr3X8+HFdvXpVO3fu1NGjRzU5OUkDlYFQKKTXX39d27dv18WLF/Xmm2/qyZMnOn36tMxmM3KmuW4AAAXASURBVA0EAAD5CeQnAAA2VXt7ux49eqTz589rampKTU1Nampq0kcffcQ0CFWSlyjSrlNtba36+/v16NEjdXV16erVq9qzZ49efvllDQ4O0oG2oOHhYR08eFB79uzRBx98kBsAz507J5vNRgMBAEB+AvkJAADDmM1mvfnmm3r8+LH6+/sVj8f129/+Vjt37tTZs2eVSqVopArOSxRpN8jhcOjKlSt68uSJent7NTExIb/fr23btqmnp0cTExMccBgoFArp1KlT2rZtm1555RVNTEzo2LFjevLkia5fv666ujoaCQAA8hPITwAAbBlms1nHjx/X48ePdfPmTVmtVp06dUo/+clP9Morr+jq1atKJpM0VIXlJYq0RTzYCAQC+uc//6kPP/xQTqdTFy9e1L59+7R9+3a99dZbHHBsYkd5++23tX37du3evVtnz57NHQx+99136u/vl8PhoKEAACA/gfwEAMCWZTab1dbWpvv37+vRo0c6duyYJicndfjwYb344ot6+eWX9f7771OwrZC8VJPJZDI0c2nE43ENDg7qk08+0cTERO5gpKWlRYcOHVJzc7NMJhMNVQTj4+P6+OOPNTw8rKmpKUlzT0ncv3+/Ojs75XQ6aaQyNDY2pqamJo2OjqqxsZEGAcC4QX4iP5GfwDgIgHGjqqXTad2+fVs3btzQ0NBQbgqExsZG7d+/X62trXyGl2leoki7SRKJhAYHB3Xjxg1NTExodnZWFotFXq9Xv/nNb9TQ0CCfzyer1UpjrWJmZkbhcFgjIyO6e/euJicnlUwmZTKZ5PP5dODAAbW3t3MrHmEBAOMGDUJ+AvmJcRAAGDcqVjqd1tjYmG7duqXbt28rFotJmjvB7fV6tX//fvl8Pnk8Hh70WQZ5icsQNondbtfx48d1/PhxpVIp3b59W3fu3NHk5KT6+vpy78sebOzYsUNer1cej0cWi6Vq2212dlbhcFjhcFiRSEShUEjj4+OanZ2VJNXV1am5uVkvvfSSmpubuQ0PAADyE/mJ/AQAQFUwmUxqbm5Wc3Ozzp8/r1gspqGhIX355ZcKhUIaHh6WNDdtQmNjo7xer9xut7xer5xOZ1XfnbQV8xJFWgNYrVZ1dnaqs7NT0txVIhMTE7p7965GRkb0/vvv53aK7I7h8Xjk8Xi0a9cuOZ3OiutM6XRa8Xhc4XBYsVhMDx8+VCwWUywWy81DZzKZ5HQ6deTIEe3du1f19fVc7QEAAPmJ/ER+AgAAkpxOp06ePJn7/1AopHA4rDt37mhiYkIjIyOLMoLH45HT6ZTb7c79u66ujrxkUF6iSLsF2O12tbW1qa2tLbcDxWKxXDU/HA4rFAppaGjo6YYzmeRwOGSz2eR0OvXzn/9cVqtVTqdTVqtVdXV1stlsW+ZvTCaTmpmZUTweVyKRUCKR0Lfffpv7/3g8vujAyuFwyOPxqKWlJXdgxeX5AACA/ER+AgAAa+P1euX1etXV1SVp7nb+WCymUCikSCSiWCym8fFxffTRR7mfMZvNstvtstvtcjqdcjgcev755+V0OmWxWMhLJUSRdgvKns3weDyLvp7tTEsr/yMjI7p27dqKBzBWqzX3X6vVKpPJpBdeeEEWiyXXsUwmk+x2+5rOlmQ7QFYqlVIqldIPP/ygZDK56JVKpZRIJPI+ldliseTO0rS2tupnP/sZtygir1AoRCMAYLwA+Yn8BMZDAIwX2ACLxSKfzyefz7fo66lUatHJ7lgspng8ruHhYSWTybzLcjgcslgs5KUi4sFhFSLboRaePfj+++81MzOjRCKR24Gnp6fz7sDFOjjKdjyHwyGr1SqbzSa73a7nnntODodj0YsDCawmO4E9ABSKB1+A/ATyEwCQn7BxyWRyUV5KJpPkpRKhSFulZmZmFp0NSafTmp6eXtPPZs+QZNlsNg4YULKDZ87qAlgPr9e76LMKID+B/AQA5CeQl7ay/wetO0EkrQIc8gAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![graph.png](attachment:graph.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "4335e3cbfd46c7579a4b7d368c5353cc",
     "grade": true,
     "grade_id": "cell-216d36a73771fdb7",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "Influenza = Variable(name='Influenza', num_states=2)\n",
    "Influenza_prior = Factor(name='p(Influenza)',\n",
    "                 f=np.array([0.95, 0.05]),\n",
    "                 neighbours=[Influenza])\n",
    "\n",
    "Smoke = Variable(name='Smoke', num_states=2)\n",
    "Smoke_prior = Factor(name='p(Smoke)',\n",
    "                 f=np.array([0.80, 0.20]),\n",
    "                 neighbours=[Smoke])\n",
    "                 \n",
    "Bronchitis = Variable(name='Bronchitis', num_states=2)\n",
    "Bronchitis_cond = Factor(name='p(Bronchitis | Influenza, Smoke)',\n",
    "                 f=np.array([\n",
    "                 [[0.9999, 0.3], [0.1, 0.01]],\n",
    "                 [[0.0001, 0.7], [0.9, 0.99]],\n",
    "                 ]),\n",
    "                 neighbours=[Bronchitis, Influenza, Smoke])  \n",
    "\n",
    "SoreThroat = Variable(name='SoreThroat', num_states=2)\n",
    "SoreThroat_cond = Factor(name='p(SoreThroat | Influenza)',\n",
    "                 f=np.array([\n",
    "                     [0.999, 0.70], [0.001, 0.30],\n",
    "                 ]),\n",
    "                 neighbours=[SoreThroat, Influenza]) \n",
    "\n",
    "Fever = Variable(name='Fever', num_states=2)\n",
    "Fever_cond = Factor(name='p(Fever | Influenza)',\n",
    "                 f=np.array([\n",
    "                     [0.95, 0.10], [0.05, 0.90],\n",
    "                 ]),\n",
    "                 neighbours=[Fever, Influenza]) \n",
    "\n",
    "Coughing = Variable(name='Coughing', num_states=2)\n",
    "Coughing_cond = Factor(name='p(Coughing | Bronchitis)',\n",
    "                 f=np.array([\n",
    "                     [0.93, 0.20], [0.07, 0.80],\n",
    "                 ]),\n",
    "                 neighbours=[Coughing, Bronchitis]) \n",
    "\n",
    "Wheezing = Variable(name='Wheezing', num_states=2)\n",
    "Wheezing_cond = Factor(name='p(Wheezing | Bronchitis)',\n",
    "                 f=np.array([\n",
    "                     [0.999, 0.40], [0.001, 0.60],\n",
    "                 ]),\n",
    "                 neighbours=[Wheezing, Bronchitis]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "8e443821bd73340b3e8bc3eb58db013c",
     "grade": false,
     "grade_id": "cell-65ee0e0d85c57f32",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 3.2. Sum-product algorithm (10 points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "48d2bd6a0a7dba5d29de2c939ca5c67f",
     "grade": false,
     "grade_id": "cell-fed678d3f9f5ce50",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "3.2.1 Run the sum-product algorithm on an unobserved graph. Print the marginal for every variable node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "a7fadef738bc786fd10febd2b25f15f5",
     "grade": true,
     "grade_id": "cell-07681153d48af573",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(Influenza) -> Influenza: [0.95 0.05]\n",
      "SoreThroat -> p(SoreThroat | Influenza): [1. 1.]\n",
      "p(Smoke) -> Smoke: [0.8 0.2]\n",
      "Coughing -> p(Coughing | Bronchitis): [1. 1.]\n",
      "Wheezing -> p(Wheezing | Bronchitis): [1. 1.]\n",
      "Fever -> p(Fever | Influenza): [1. 1.]\n",
      "p(SoreThroat | Influenza) -> Influenza: [1. 1.]\n",
      "p(Fever | Influenza) -> Influenza: [1. 1.]\n",
      "p(Coughing | Bronchitis) -> Bronchitis: [1. 1.]\n",
      "p(Wheezing | Bronchitis) -> Bronchitis: [1. 1.]\n",
      "Smoke -> p(Bronchitis | Influenza, Smoke): [0.8 0.2]\n",
      "Bronchitis -> p(Bronchitis | Influenza, Smoke): [1. 1.]\n",
      "Influenza -> p(Bronchitis | Influenza, Smoke): [0.95 0.05]\n",
      "\n",
      "Backward pass\n",
      "p(Bronchitis | Influenza, Smoke) -> Influenza: [1. 1.]\n",
      "p(Bronchitis | Influenza, Smoke) -> Bronchitis: [0.821024 0.178976]\n",
      "p(Bronchitis | Influenza, Smoke) -> Smoke: [1. 1.]\n",
      "Influenza -> p(SoreThroat | Influenza): [0.95 0.05]\n",
      "Influenza -> p(Influenza): [1. 1.]\n",
      "Influenza -> p(Fever | Influenza): [0.95 0.05]\n",
      "Bronchitis -> p(Wheezing | Bronchitis): [0.821024 0.178976]\n",
      "Bronchitis -> p(Coughing | Bronchitis): [0.821024 0.178976]\n",
      "Smoke -> p(Smoke): [1. 1.]\n",
      "p(Wheezing | Bronchitis) -> Wheezing: [0.89179338 0.10820662]\n",
      "p(Coughing | Bronchitis) -> Coughing: [0.79934752 0.20065248]\n",
      "p(Fever | Influenza) -> Fever: [0.9075 0.0925]\n",
      "p(SoreThroat | Influenza) -> SoreThroat: [0.98405 0.01595]\n",
      "\n",
      "\n",
      "Marginal of SoreThroat is [0.98405 0.01595]\n",
      "Marginal of Coughing is [0.79934752 0.20065248]\n",
      "Marginal of Wheezing is [0.89179338 0.10820662]\n",
      "Marginal of Fever is [0.9075 0.0925]\n",
      "Marginal of Smoke is [0.8 0.2]\n",
      "Marginal of Bronchitis is [0.821024 0.178976]\n",
      "Marginal of Influenza is [0.95 0.05]\n"
     ]
    }
   ],
   "source": [
    "# maybe I just ordered it in a bad way..?\n",
    "nodes = [Influenza_prior, SoreThroat, Smoke_prior, Coughing,\n",
    "         Wheezing, Fever, SoreThroat_cond, Fever_cond, Coughing_cond, Wheezing_cond, \n",
    "         Smoke, Bronchitis, Influenza, Bronchitis_cond]\n",
    "\n",
    "leaves = [Influenza_prior, SoreThroat, Fever, Smoke_prior, Coughing, Wheezing]\n",
    "\n",
    "for n in nodes:\n",
    "    n.reset()\n",
    "\n",
    "for l in leaves:\n",
    "    for n in l.neighbours:\n",
    "        l.pending.add(n)\n",
    "\n",
    "sum_product(nodes)\n",
    "\n",
    "print('\\n')\n",
    "for n in nodes:\n",
    "    if type(n) == Variable:\n",
    "        print(f'Marginal of {n} is {n.marginal()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "7ab62562db323a454298e365f0989326",
     "grade": false,
     "grade_id": "cell-9cf9ad87ef338231",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "3.2.2 Rerun the sum-product algorithm on an partially observed graph, where the variable 'Influenza' is set to 0. Print the marginal for every variable node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0bf552eaf0812c495b3cb3c79f2fc431",
     "grade": true,
     "grade_id": "cell-83214d3cc28cb7db",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(Influenza) -> Influenza: [0.95 0.05]\n",
      "SoreThroat -> p(SoreThroat | Influenza): [1. 1.]\n",
      "p(Smoke) -> Smoke: [0.8 0.2]\n",
      "Coughing -> p(Coughing | Bronchitis): [1. 1.]\n",
      "Wheezing -> p(Wheezing | Bronchitis): [1. 1.]\n",
      "Fever -> p(Fever | Influenza): [1. 1.]\n",
      "p(SoreThroat | Influenza) -> Influenza: [1. 1.]\n",
      "p(Fever | Influenza) -> Influenza: [1. 1.]\n",
      "p(Coughing | Bronchitis) -> Bronchitis: [1. 1.]\n",
      "p(Wheezing | Bronchitis) -> Bronchitis: [1. 1.]\n",
      "Smoke -> p(Bronchitis | Influenza, Smoke): [0.8 0.2]\n",
      "Bronchitis -> p(Bronchitis | Influenza, Smoke): [1. 1.]\n",
      "Influenza -> p(Bronchitis | Influenza, Smoke): [9.5e-01 5.0e-12]\n",
      "\n",
      "Backward pass\n",
      "p(Bronchitis | Influenza, Smoke) -> Influenza: [1. 1.]\n",
      "p(Bronchitis | Influenza, Smoke) -> Bronchitis: [0.816924 0.133076]\n",
      "p(Bronchitis | Influenza, Smoke) -> Smoke: [0.95 0.95]\n",
      "Influenza -> p(SoreThroat | Influenza): [9.5e-01 5.0e-12]\n",
      "Influenza -> p(Influenza): [1.e+00 1.e-10]\n",
      "Influenza -> p(Fever | Influenza): [9.5e-01 5.0e-12]\n",
      "Bronchitis -> p(Wheezing | Bronchitis): [0.816924 0.133076]\n",
      "Bronchitis -> p(Coughing | Bronchitis): [0.816924 0.133076]\n",
      "Smoke -> p(Smoke): [0.95 0.95]\n",
      "p(Wheezing | Bronchitis) -> Wheezing: [0.86933748 0.08066252]\n",
      "p(Coughing | Bronchitis) -> Coughing: [0.78635452 0.16364548]\n",
      "p(Fever | Influenza) -> Fever: [0.9025 0.0475]\n",
      "p(SoreThroat | Influenza) -> SoreThroat: [0.94905 0.00095]\n",
      "\n",
      "\n",
      "Marginal of SoreThroat is [0.999 0.001]\n",
      "Marginal of Coughing is [0.8277416 0.1722584]\n",
      "Marginal of Wheezing is [0.91509208 0.08490792]\n",
      "Marginal of Fever is [0.95 0.05]\n",
      "Marginal of Smoke is [0.8 0.2]\n",
      "Marginal of Bronchitis is [0.85992 0.14008]\n",
      "Marginal of Influenza is [1.00000000e+00 5.26315789e-12]\n"
     ]
    }
   ],
   "source": [
    "for n in nodes:\n",
    "    n.reset()\n",
    "\n",
    "for l in leaves:\n",
    "    for n in l.neighbours:\n",
    "        l.pending.add(n)\n",
    "\n",
    "Influenza.set_observed(0)\n",
    "\n",
    "sum_product(nodes)\n",
    "\n",
    "print('\\n')\n",
    "for n in nodes:\n",
    "    if type(n) == Variable:\n",
    "        print(f'Marginal of {n} is {n.marginal()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d5c513d1f3fcdd0e153834c8cd7013e7",
     "grade": false,
     "grade_id": "cell-c2dbbc1def846b5e",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "### 3.3. max-sum algorithm (10 points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "e3397f5416e8ef7a287fd0e1844d5ec6",
     "grade": false,
     "grade_id": "cell-2aeb9d6f8bb0d6d9",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "3.3.1 Run the max_sum algorithm on an unobserved graph. Print the marginal for every variable node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "350a548fee4372022d9d7b51805ae71a",
     "grade": true,
     "grade_id": "cell-e5088e48afe006b9",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(Influenza) -> Influenza: [0.95 0.05]\n",
      "SoreThroat -> p(SoreThroat | Influenza): [1. 1.]\n",
      "p(Smoke) -> Smoke: [0.8 0.2]\n",
      "Coughing -> p(Coughing | Bronchitis): [1. 1.]\n",
      "Wheezing -> p(Wheezing | Bronchitis): [1. 1.]\n",
      "Fever -> p(Fever | Influenza): [1. 1.]\n",
      "p(SoreThroat | Influenza) -> Influenza: [0.999 0.7  ]\n",
      "p(Fever | Influenza) -> Influenza: [0.95 0.9 ]\n",
      "p(Coughing | Bronchitis) -> Bronchitis: [0.93 0.8 ]\n",
      "p(Wheezing | Bronchitis) -> Bronchitis: [0.999 0.6  ]\n",
      "Smoke -> p(Bronchitis | Influenza, Smoke): [0.8 0.2]\n",
      "Bronchitis -> p(Bronchitis | Influenza, Smoke): [0.92907 0.48   ]\n",
      "Influenza -> p(Bronchitis | Influenza, Smoke): [0.9015975 0.0315   ]\n",
      "\n",
      "Backward pass\n",
      "p(Bronchitis | Influenza, Smoke) -> Influenza: [0.74318167 0.3456    ]\n",
      "p(Bronchitis | Influenza, Smoke) -> Bronchitis: [0.72120587 0.12622365]\n",
      "p(Bronchitis | Influenza, Smoke) -> Smoke: [0.83756342 0.30293676]\n",
      "Influenza -> p(SoreThroat | Influenza): [0.67072146 0.015552  ]\n",
      "Influenza -> p(Influenza): [0.70531657 0.217728  ]\n",
      "Influenza -> p(Fever | Influenza): [0.70531657 0.012096  ]\n",
      "Bronchitis -> p(Wheezing | Bronchitis): [0.67072146 0.10097892]\n",
      "Bronchitis -> p(Coughing | Bronchitis): [0.72048467 0.07573419]\n",
      "Smoke -> p(Smoke): [0.83756342 0.30293676]\n",
      "p(Wheezing | Bronchitis) -> Wheezing: [0.67005074 0.06058735]\n",
      "p(Coughing | Bronchitis) -> Coughing: [0.67005074 0.06058735]\n",
      "p(Fever | Influenza) -> Fever: [0.67005074 0.03526583]\n",
      "p(SoreThroat | Influenza) -> SoreThroat: [0.67005074 0.0046656 ]\n",
      "\n",
      "\n",
      "Unnormalized log-marginal of SoreThroat is [-0.40040184 -5.36753884]\n",
      "Unnormalized log-marginal of Coughing is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Wheezing is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Fever is [-0.40040184 -3.34484082]\n",
      "Unnormalized log-marginal of Smoke is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Bronchitis is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Influenza is [-0.40040184 -4.52024098]\n"
     ]
    }
   ],
   "source": [
    "for n in nodes:\n",
    "    n.reset()\n",
    "\n",
    "for l in leaves:\n",
    "    for n in l.neighbours:\n",
    "        l.pending.add(n)\n",
    "\n",
    "max_sum(nodes)\n",
    "\n",
    "print('\\n')\n",
    "for n in nodes:\n",
    "    if type(n) == Variable:\n",
    "        print(f'Unnormalized log-marginal of {n} is {n.unnormalized_log_marginal()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "083775bc40177f2573543fd168e86e09",
     "grade": false,
     "grade_id": "cell-7e7771be2511ebd2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "3.3.2 Rerun the max_sum algorithm on an partially observed graph, where the variable 'Influenza' is set to 0. Print the marginal for every variable node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "b467a95a6dbe003435c78b970ed0e318",
     "grade": true,
     "grade_id": "cell-f74e50a37f4e31b9",
     "locked": false,
     "points": 5,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward pass\n",
      "p(Influenza) -> Influenza: [9.5e-01 5.0e-12]\n",
      "SoreThroat -> p(SoreThroat | Influenza): [1. 1.]\n",
      "p(Smoke) -> Smoke: [0.8 0.2]\n",
      "Coughing -> p(Coughing | Bronchitis): [1. 1.]\n",
      "Wheezing -> p(Wheezing | Bronchitis): [1. 1.]\n",
      "Fever -> p(Fever | Influenza): [1. 1.]\n",
      "p(SoreThroat | Influenza) -> Influenza: [0.999 0.7  ]\n",
      "p(Fever | Influenza) -> Influenza: [0.95 0.1 ]\n",
      "p(Coughing | Bronchitis) -> Bronchitis: [0.93 0.8 ]\n",
      "p(Wheezing | Bronchitis) -> Bronchitis: [0.999 0.6  ]\n",
      "Smoke -> p(Bronchitis | Influenza, Smoke): [0.8 0.2]\n",
      "Bronchitis -> p(Bronchitis | Influenza, Smoke): [0.92907 0.48   ]\n",
      "Influenza -> p(Bronchitis | Influenza, Smoke): [9.015975e-01 3.500000e-23]\n",
      "\n",
      "Backward pass\n",
      "p(Bronchitis | Influenza, Smoke) -> Influenza: [0.74318167 0.3456    ]\n",
      "p(Bronchitis | Influenza, Smoke) -> Bronchitis: [0.72120587 0.12622365]\n",
      "p(Bronchitis | Influenza, Smoke) -> Smoke: [0.83756342 0.30293676]\n",
      "Influenza -> p(SoreThroat | Influenza): [6.70721461e-01 1.72800000e-23]\n",
      "Influenza -> p(Influenza): [7.05316568e-01 2.41920000e-12]\n",
      "Influenza -> p(Fever | Influenza): [7.05316568e-01 1.20960000e-22]\n",
      "Bronchitis -> p(Wheezing | Bronchitis): [0.67072146 0.10097892]\n",
      "Bronchitis -> p(Coughing | Bronchitis): [0.72048467 0.07573419]\n",
      "Smoke -> p(Smoke): [0.83756342 0.30293676]\n",
      "p(Wheezing | Bronchitis) -> Wheezing: [0.67005074 0.06058735]\n",
      "p(Coughing | Bronchitis) -> Coughing: [0.67005074 0.06058735]\n",
      "p(Fever | Influenza) -> Fever: [0.67005074 0.03526583]\n",
      "p(SoreThroat | Influenza) -> SoreThroat: [0.67005074 0.00067072]\n",
      "\n",
      "\n",
      "Unnormalized log-marginal of SoreThroat is [-0.40040184 -7.30715662]\n",
      "Unnormalized log-marginal of Coughing is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Wheezing is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Fever is [-0.40040184 -3.34484082]\n",
      "Unnormalized log-marginal of Smoke is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Bronchitis is [-0.40040184 -2.80366912]\n",
      "Unnormalized log-marginal of Influenza is [ -0.40040184 -52.76916741]\n"
     ]
    }
   ],
   "source": [
    "for n in nodes:\n",
    "    n.reset()\n",
    "\n",
    "for l in leaves:\n",
    "    for n in l.neighbours:\n",
    "        l.pending.add(n)\n",
    "\n",
    "Influenza.set_observed(0)\n",
    "\n",
    "max_sum(nodes)\n",
    "\n",
    "print('\\n')\n",
    "for n in nodes:\n",
    "    if type(n) == Variable:\n",
    "        print(f'Unnormalized log-marginal of {n} is {n.unnormalized_log_marginal()}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ml2labs]",
   "language": "python",
   "name": "conda-env-ml2labs-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}