4-ImplementationBP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation - Belief Propagation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy  as np\n",
    "import igraph as ig\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc(\"text\", usetex=True)\n",
    "\n",
    "%run ./2-ImplementationFactor.ipynb\n",
    "%run ./3-ImplementationPGM.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 BP\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a class ```belief_propagation``` and store messages in ```self.__msg ``` since they can be reused in further calculations. Also there is ```self.__pgm``` field for a factor graph.\n",
    "\n",
    "Note, that we normalize each message in order to avoid numerical underflow/overflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class belief_propagation():\n",
    "    def __init__(self, pgm):\n",
    "        if type(pgm) is not factor_graph:\n",
    "            raise Exception('PGM is not a factor graph')\n",
    "        if not (pgm.is_connected() and not pgm.is_loop()):\n",
    "            raise Exception('PGM is not a tree')\n",
    "        \n",
    "        self.__msg = {}\n",
    "        self.__pgm = pgm\n",
    "    \n",
    "    def belief(self, v_name):\n",
    "        incoming_messages = []\n",
    "        for f_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(v_name)]['name']:\n",
    "            incoming_messages.append(self.get_factor2variable_msg(f_name_neighbor, v_name))\n",
    "        return self.__normalize_msg(joint_distribution(incoming_messages))\n",
    "    \n",
    "    # ----------------------- Variable to factor ------------\n",
    "    def get_variable2factor_msg(self, v_name, f_name):\n",
    "        key = (v_name, f_name)\n",
    "        if key not in self.__msg:\n",
    "            self.__msg[key] = self.__compute_variable2factor_msg(v_name, f_name)\n",
    "        return self.__msg[key]\n",
    "    \n",
    "    def __compute_variable2factor_msg(self, v_name, f_name):\n",
    "        incoming_messages = []\n",
    "        for f_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(v_name)]['name']:\n",
    "            if f_name_neighbor != f_name:\n",
    "                incoming_messages.append(self.get_factor2variable_msg(f_name_neighbor, v_name))\n",
    "        \n",
    "        if not incoming_messages:\n",
    "            # if the variable does not have its own distribution\n",
    "            return factor([v_name], np.array([1.]*self.__pgm.get_graph().vs.find(name=v_name)['rank']))\n",
    "        else:\n",
    "            # Since all messages have the same dimension (1, order of v_name) the expression after\n",
    "            # ```return``` is equivalent to ```factor(v_name, np.prod(incoming_messages))```\n",
    "            return self.__normalize_msg(joint_distribution(incoming_messages))\n",
    "    \n",
    "    # ----------------------- Factor to variable ------------\n",
    "    def get_factor2variable_msg(self, f_name, v_name):\n",
    "        key = (f_name, v_name)\n",
    "        if key not in self.__msg:\n",
    "            self.__msg[key] = self.__compute_factor2variable_msg(f_name, v_name)\n",
    "        return self.__msg[key]\n",
    "    \n",
    "    def __compute_factor2variable_msg(self, f_name, v_name):\n",
    "        incoming_messages = [self.__pgm.get_graph().vs.find(f_name)['factor_']]\n",
    "        marginalization_variables = []\n",
    "        for v_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(f_name)]['name']:\n",
    "            if v_name_neighbor != v_name:\n",
    "                incoming_messages.append(self.get_variable2factor_msg(v_name_neighbor, f_name))\n",
    "                marginalization_variables.append(v_name_neighbor)\n",
    "        return self.__normalize_msg(factor_marginalization(\n",
    "            joint_distribution(incoming_messages),\n",
    "            marginalization_variables\n",
    "        ))\n",
    "    \n",
    "    # ----------------------- Other -------------------------\n",
    "    def __normalize_msg(self, message):\n",
    "        return factor(message.get_variables(), message.get_distribution()/np.sum(message.get_distribution()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "Let's take a look at the following factor graph."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./img/4/1.1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrf = string2factor_graph('f1(a,b)f2(b,c,d)f3(c)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = factor(['a', 'b'],      np.array([[2,3],[6,4]]))\n",
    "f2 = factor(['b', 'd', 'c'], np.array([[[7,2,3],[1,5,2]],[[8,3,9],[6,4,2]]]))\n",
    "f3 = factor(['c'],           np.array([5, 1, 9]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrf.change_factor_distribution('f1', f1)\n",
    "mrf.change_factor_distribution('f2', f2)\n",
    "mrf.change_factor_distribution('f3', f3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the algorithm to compute marginal distribution of ```b```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.37398374, 0.62601626])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bp = belief_propagation(mrf)\n",
    "bp.belief('b').get_distribution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Loopy BP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class loopy_belief_propagation():\n",
    "    def __init__(self, pgm):\n",
    "        if type(pgm) is not factor_graph:\n",
    "            raise Exception('PGM is not a factor graph')\n",
    "        if not pgm.is_connected():\n",
    "            raise Exception('PGM is not connected')\n",
    "        if len(pgm.get_graph().es) - 1 == len(pgm.get_graph().vs):\n",
    "            raise Exception('PGM is a tree')\n",
    "        \n",
    "        self.__t       = 0\n",
    "        self.__msg     = {}\n",
    "        self.__msg_new = {}\n",
    "        self.__pgm     = pgm\n",
    "        \n",
    "        # Initialization of messages\n",
    "        for edge in self.__pgm.get_graph().es:\n",
    "            start_index, end_index = edge.tuple[0], edge.tuple[1]\n",
    "            start_name, end_name = self.__pgm.get_graph().vs[start_index]['name'], self.__pgm.get_graph().vs[end_index]['name']\n",
    "            \n",
    "            if self.__pgm.get_graph().vs[start_index]['is_factor']:\n",
    "                self.__msg[(start_name, end_name)] = factor([end_name],   np.array([1.]*self.__pgm.get_graph().vs[end_index]['rank']))\n",
    "            else:\n",
    "                self.__msg[(start_name, end_name)] = factor([start_name], np.array([1.]*self.__pgm.get_graph().vs[start_index]['rank']))\n",
    "            self.__msg[(end_name, start_name)] = self.__msg[(start_name, end_name)]\n",
    "            \n",
    "            self.__msg_new[(start_name, end_name)] = 0\n",
    "            self.__msg_new[(end_name, start_name)] = 0\n",
    "    \n",
    "    def belief(self, v_name, num_iter):\n",
    "        if self.__t > num_iter:\n",
    "            raise Exception('Invalid number of iterations. Current number: ' + str(self.__t))\n",
    "        elif self.__t < num_iter:\n",
    "            self.__loop(num_iter)\n",
    "        \n",
    "        incoming_messages = []\n",
    "        for f_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(v_name)]['name']:\n",
    "            incoming_messages.append(self.get_factor2variable_msg(f_name_neighbor, v_name))\n",
    "        return self.__normalize_msg(joint_distribution(incoming_messages))\n",
    "    \n",
    "    # ----------------------- Variable to factor ------------\n",
    "    def get_variable2factor_msg(self, v_name, f_name):\n",
    "        return self.__msg[(v_name, f_name)]\n",
    "    \n",
    "    def __compute_variable2factor_msg(self, v_name, f_name):\n",
    "        incoming_messages = []\n",
    "        for f_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(v_name)]['name']:\n",
    "            if f_name_neighbor != f_name:\n",
    "                incoming_messages.append(self.get_factor2variable_msg(f_name_neighbor, v_name))\n",
    "        \n",
    "        if not incoming_messages:\n",
    "            return factor([v_name], np.array([1]*self.__pgm.get_graph().vs.find(name=v_name)['rank']))\n",
    "        else:\n",
    "            return self.__normalize_msg(joint_distribution(incoming_messages))\n",
    "    \n",
    "    # ----------------------- Factor to variable ------------\n",
    "    def get_factor2variable_msg(self, f_name, v_name):\n",
    "        return self.__msg[(f_name, v_name)]\n",
    "    \n",
    "    def __compute_factor2variable_msg(self, f_name, v_name):\n",
    "        incoming_messages = [self.__pgm.get_graph().vs.find(f_name)['factor_']]\n",
    "        marginalization_variables = []\n",
    "        for v_name_neighbor in self.__pgm.get_graph().vs[self.__pgm.get_graph().neighbors(f_name)]['name']:\n",
    "            if v_name_neighbor != v_name:\n",
    "                incoming_messages.append(self.get_variable2factor_msg(v_name_neighbor, f_name))\n",
    "                marginalization_variables.append(v_name_neighbor)\n",
    "        return self.__normalize_msg(factor_marginalization(\n",
    "            joint_distribution(incoming_messages),\n",
    "            marginalization_variables\n",
    "        ))\n",
    "    \n",
    "    # ----------------------- Other -------------------------\n",
    "    def __loop(self, num_iter):\n",
    "        # Message updating\n",
    "        while self.__t < num_iter:\n",
    "            for edge in self.__pgm.get_graph().es:\n",
    "                start_index, end_index = edge.tuple[0], edge.tuple[1]\n",
    "                start_name, end_name   = self.__pgm.get_graph().vs[start_index]['name'], self.__pgm.get_graph().vs[end_index]['name']\n",
    "\n",
    "                if self.__pgm.get_graph().vs[start_index]['is_factor']:\n",
    "                    self.__msg_new[(start_name, end_name)] = self.__compute_factor2variable_msg(start_name, end_name)\n",
    "                    self.__msg_new[(end_name, start_name)] = self.__compute_variable2factor_msg(end_name, start_name)\n",
    "                else:\n",
    "                    self.__msg_new[(start_name, end_name)] = self.__compute_variable2factor_msg(start_name, end_name)\n",
    "                    self.__msg_new[(end_name, start_name)] = self.__compute_factor2variable_msg(end_name, start_name)\n",
    "            self.__msg.update(self.__msg_new)\n",
    "            self.__t += 1\n",
    "    \n",
    "    def __normalize_msg(self, message):\n",
    "        return factor(message.get_variables(), message.get_distribution()/np.sum(message.get_distribution()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "Let's consider the following factor graph with a loop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](./img/4/1.2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrf = string2factor_graph('f1(a,b)f2(a,c)f3(b,c)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = factor(['a', 'b'],  np.array([[2,3],[6,4]]))\n",
    "f2 = factor(['a', 'c'],  np.array([[7,2,3],[1,5,2]]))\n",
    "f3 = factor(['b', 'c'],  np.array([[7,9,3],[6,4,2]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mrf.change_factor_distribution('f1', f1)\n",
    "mrf.change_factor_distribution('f2', f2)\n",
    "mrf.change_factor_distribution('f3', f3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, let's compute exact marginal distribution of `b` and normalize it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.63451777, 0.36548223])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exact = factor_marginalization(joint_distribution([f1, f2, f3]), ['a', 'c']).get_distribution()\n",
    "exact = exact / np.sum(exact)\n",
    "exact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Secondly, we run the loopy belief propagation algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "lbp = loopy_belief_propagation(mrf)\n",
    "tol = []\n",
    "\n",
    "for i in range(15):\n",
    "    tol.append(np.linalg.norm(lbp.belief('b', i).get_distribution() - exact))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xa1b876a90>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgEAAAF6CAYAAACA64v+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhV930n/vfnLlrRgkASAswiCWNsDFgI747HWE76S+M2TcBu+0t/M8/UiKSZp+k0GTDpr52n/XWKsdt0SdIYOe0z02lngsHpNDQrwqkXjBchMF4Agy4Isy9CC1rvvefz++MekYsQkq50zvnee8/79Tx6kI6ke95HQud8zvd8F1FVEBERkf8ETAcgIiIiM1gEEBER+RSLACIiIp9iEUBERORTLAKIiIh8ikUAERGRT4VMB/DSzJkzdcGCBaZjEBEReWbfvn2XVLV8tM/5qghYsGABWlpaTMcgIiLyjIi03+xzfBxARETkUywCiIiIfIpFABERkU+xCCAiIvIpFgFEREQ+xSKAiIjIp3w1RJCIiIDu7m5cuHAB0WjUdBSagnA4jIqKChQXF0/6NVgEEBH5SHd3N86fP485c+YgPz8fImI6Ek2CqqK/vx+nT58GgEkXAnwcQETkIxcuXMCcOXNQUFDAAiCDiQgKCgowZ84cXLhwYdKv44siQEQeF5Gmrq4u01GIiIyKRqPIz883HWNKOjs7R30/0zhxHPn5+VN6rOOLIkBVd6pqY0lJiekoRETGZXoLwObNmwEkLpxNTU3Xtjc3N6Ompua6bV6ZzL5vdhypmOrv0hdFABERZb+GhgaUlpaioaEhrfb97LPPep5nolgETNLVviiqf+ef8YU/fdl0FCIiskUiEVRXV6fFviORCJ599llUV1ejqakJra2tRnKNhUXAJE0rCCMeEhy92Gs6ChERAWhtbUV9fX3a7Lu6uhpr1qzB5s2bUV1djbq6OiPZxsIhglOQH1N0IGY6BhERIfFcvrS0FM3NzWhtbUVdXZ1njwZG23d1dTV27NiBTZs2IRKJXNueTlgETEF5bhgnNQ7LshAIsFGFiMikXbt2Ye3atWhoaEBDQwNqamrQ1taG1tZWrFu3Dvv27fN83xs2bEBTUxMaGxtd2/dUsAiYggUzCvHx1T4c+7gLt86fbjoOEZGvtbS0YNeuXQASPe4jkQgAXLsrN7FvACkVAJ2dnWhubkZHRwd27dqFLVu2uJqdRcAU3LmgFK+934fX3j3HIoCIMtbv/d5PcODAOSP7XrFiFv7qr35pyq/T2tp6Q6e84Y+bm5uxatWqcV+jqalp3NaCjRs33nBRHmvfqXrxxRfR2dmJDRs2AAC2bNmCrVu3Tuq1JoJFwBTcd0cl/vb9M9h/7LLpKEREvtbR0XFdx7xt27Zh/fr1AHDtWXwkEkFzczOeeOIJlJaW3vAak22yH2vfqUrO0NbWhpqamkm9zkSxCJiCe5ZWQv/JwtHzPaajEBFNmhN34qYl33kPN6kP39Xv2rULa9asufa5jo6OUYsAN/Y9FZFIBFu2bJny64yFRcAU5ISDCPTHcBZx01GIiHyturoaNTU12LFjByKRCHbv3n3tc5FIBDt27MCGDRuuNbN7te/JevbZZ7F9+3YH0o2NRcAUFSOAblimYxAR+d5oF/hIJIKGhgY0NjZi/fr12LJli6OtAGPte7J27Nhx7bFAc3Ozq8McOa5timYX5cIqCKNvgOtyExGlm9bWVqxduxalpaUoLS29rtd+OmptbcXGjRvx6KOPoqamxvW8bAmYoltnFeHQuU68+d55rF4113QcIiJKMtwXAIDrz9edUFdXh7a2Ns/2x5aAKapbNAMA8OaHk1/PmYiIJseNpn0TTB2HqKqRHZtQX1+vLS0tjr7m8dPdeOSbr+H+gnz8rz9a7ehrExE57dChQ1iyZInpGOSg8X6nIrJPVUddVIGPA6Zo4ZxiaH8M7QNcSIiIiDILiwAH5A5ZuMQRAkRElGHYJ8ABM3JCGMjhj5KIiDILr1wOWFBWAMkP4eQ5zhxIRESZg0WAA+6Yl+jV+er+M4aTEBGNz08dwrPdVH+XLAIccM+SCgBA61EuJERE6S0cDqO/v990DHJIf38/wuHwpL+fRYAD7ls2C2opjpztNh2FiGhMFRUVOH36NPr6+tgikMFUFX19fTh9+jQqKiom/TocHeCAaQVhSG8UZ4QLCRFReisuLgYAnDlzBtEopzvPZOFwGJWVldd+p5PBIsAhRRB0WRwmSETpr7i4eEoXDsoefBzgkFmFuYgXhDAUZWsAERFlBhYBDllUOQ0SCqCFawgQEVGGYBHgkBU1iYWE9r7PIoCIiDIDiwCHPLhsFgDgYHuH4SREREQTwyLAIbctKIUOxnH8Up/pKERERBPC0QEOCQQCCA/GcZELCRERUYZgS4CDykJB9AVNpyAiIpoYFgEOmleaDykM40IHHwkQEVH6YxHgoCVzEwsJvbL/rOEkRERE42MR4KC7b5sJAGg5ctFwEiIiovGxCHDQgyuqoKo4fJoLCRERUfrj6AAHTS/Og/TGcEo4QoCIiNIfiwCHFVpAJ7h+ABERpT8+DnBYZUEYsbwgLK4oSEREaY5FgMNqyqdBcoI48NFl01GIiIjGxCLAYcsWlgEA9hw8ZzgJERHR2FgEOOyB4YWEjnMhISIiSm8ZWwSISIOI7DKdY6QVt86ADsXRdvGq6ShERERjytgiQFWbTWcYTSAQQGggjvN9UdNRiIiIxuRJESAidSKyXUQaRmwvFZENIrLG/rfOizxuKw0G0Zux5RUREfmF6/MEJF34q0f59HYA61U1Yn/tLhFZq6qdbudy09ziPFyKDeFK9wCmF+eZjkNERDQq14uA4WZ7Ebmup5yIlAKoHi4AbBEADQB2iEjjTV6vya2sTrltTjHePXkZrx84i8c/sdB0HCIiolGZnDGwHsDIO/5OAI8B2JEJF/ubqV9cjm0nL+Ptw5dYBBARUdoy+eS6FMDIcXSXAZRN5JtFZA2AahFptFsV0sbDd1UBAA6dyuinGkRElOVMrx0woQv+aFR1B4Ad432d/VihEQDmzZs32d2lpKKsANobxck+riFARETpy2RLQCcSrQHJZuDG1oEpUdUmVa1X1fry8nInX3pMBXGgI8YigIiI0pfJIqAFN7YElAJIuwmAJqM8L4xoLhcSIiKi9GWsCLCHAbaISPLQwXoAaTkJUKoWziyA5AZx+AT7BRARUXpyvQiwJwragMQFfqP9/rC1ANYMTxYEYF2mzxEwbNn8RCPH61xIiIiI0pQX8wS0AmgF8Owon+scbXs2uG9pBb753ikcaOOSwkRElJ44ua1L6m+vgMYsHD3PhYSIiCg9mR4imLVywkEE+2I4JxwhQERE6ckXLQEi8riINHV1dXm635JAAD1QT/dJREQ0Ub4oAlR1p6o2lpSUeLrf2UW50MIwrnJZYSIiSkO+KAJMWVxVDAkI9nKEABERpSEWAS6qWzQDAPDWoQuGkxAREd2IRYCLPnHXbADAByezYuoDIiLKMhwd4KJ5s4qg/TGc6O8zHYWIiOgGLAJcljdk4TK4fgAREaUfPg5w2czcEAZz+GMmIqL044urk6l5AgBgflkhJD+E46e7Pd83ERHRWHxRBJiaJwAAls4vBQC8euCM5/smIiIaiy+KAJPuvb0CANB6lAsJERFRemER4LJ776yExhUfnesxHYWIiOg6HB3gsoK8MAJ9UZwBFxIiIqL0wpYADxRB0M1hgkRElGZYBHigalourPwQhqJsDSAiovTBIsADiyqLIKEA3nr/vOkoRERE17AI8MBdtYmFhPZ+wCKAiIjSB4sADzy0fBYA4L0TXEiIiIjShy+KAJMzBgJA7S0l0IEYTlzuNbJ/IiKi0fiiCDA5YyAABAIB5AxZuDgYNbJ/IiKi0fiiCEgHZeEQ+kNiOgYREdE1LAI8Mr80H1IQxpmLfCRARETpgUWAR26/xV5IaP9Zw0mIiIgSWAR45O4l5QCAfR9dNJyEiIgogUWARx5aUQW1FEfOdJuOQkREBIALCHmmqDAH0hfDKS4kREREaYJFgIemWUAnFxIiIqI0wccBHqosyEE8P4RYjIUAERGZxyLAQ7UV0yDhAPYfYedAIiIyzxdFgOlpg4ctry4DALzxHhcSIiIi83xRBJieNnjYA8sSCwm9e+KK0RxERESAT4qAdHFnbRl0KI7IxaumoxAREXF0gJcCgQDCA3Fc4AgBIiJKA2wJ8FhpMIjeoOkURERELAI8d0tpPqQwjMtdA6ajEBGRz7EI8NiSOcUAgNf2nzGchIiI/I5FgMfqFycWEnrnyCXDSYiIyO9YBHjsoRVVAIBDp8zOWUBERMTRAR4rn54P7Y3iYy4kREREhrEIMKAgDlxhEUBERIbxcYABFXlhRPOCsCzOF0BEROawCDCgurwQkhPEB5EO01GIiMjHfFEEpMsCQsOWLUgsJPT6u1xIiIiIzJlwESAixSLyvIgcE5G37W3rRORr7sVzRrosIDTs/jsrAADvRi4bTkJERH6WSkvAdgDHVLV2+PtU9QUA690Ils1W3lYBjVo4doELCRERkTmpjA5Ypaqfst/XpO0zHMzjC6FQAMH+GM5xhAARERmUSktAi4j8dvIGEfk8gBZnI/lDaSCAq77okUFEROkqlZaAJ5AoBJ4GUCYi2wA0AFjpSrIsN6c4D5djUfT0DqGoMMd0HCIi8qEJ34uqaqfdH+CLAJ4B0AygWlVPuJQtqy2eXQwJCPa8e850FCIi8qlURgdsBgBV3a2qz9mdAquHt1NqVi6aCQB469AFw0mIiMivUnkqvWbkBlXdP9p2Gt/DdbMBAB9+3Gk4CRER+dW4fQKS7vRHu+tfCaDM8VQ+MLu8ENoXQ3s/pw4mIiIzJtIS0GG/yShvuwDUu5Yuy+XFLFweipmOQUREPjVuS4CqPgcAIrJKVZ92P5J/lOeG8bFyrgAiIjIjldEBT7gZxI8WzCiA5IVw7CT7BRARkfdSXTvgOyLyjogcTXr7qZsBs9md86cDAF5996zhJERE5Eeprh0wA4k5AgSJ+QJeBhBxIZcv3Ht7YiGh/ce4pDAREXkvlRkD61V1BgCIyHpV3Q1gt4i840607HfP0kro9ywcPddtOgoREflQKi0BV5Le7xSR5SJSDKDO4UyOE5HHRaSpq6vLdJTr5OWGEOiL4czVQdNRiIjIh1IpAraKyOfs958B8HMAxwHsdjyVw1R1p6o2lpSUmI5yg2IE0HPdooxERETemPDjgOGhgvb7rSKyEIm1A/a7kswnZhflojMexcBgDHm5qTydISIimppURgf81G7+BwCoahcLgKlbNKsIEgxg73tcSIiIiLyVyuOALgBr3QriV3W1MwAAb37IhYSIiMhbqbQ/bwbQLCI1SCwjfI2qvuxoKh95aEUV8Hob3m/nhEFEROStVIqAF5AYIfCE/TZMASxyMpSf1MwtgQ7EcGKgz3QUIiLymVQ6BnKhIJfkDlq4JFxNkIiIvJVKnwByyYycEAbC/FUQEZG3eOVJA/PLCiD5IZw6f9V0FCIi8hEWAWng9ltKAQCv7j9jOAkREfkJi4A0cM+ScgDAvqOXDCchIiI/SXUp4edF5JiIvG1vWyciX3Mvnj88sLwKaimOnOkxHYWIiHwk1aWEj6lq7fD3qeoLANa7EcxPphWEIX0xnO4ZMB2FiIh8JJV5Alap6qfs95NXvJnhYB7fKrKALnCYIBEReSeVloAWEfnt5A0i8nkALc5G8qdZhTmIF4QQi7EQICIib6RSBDwBYJOIHAVQLSLbADQBaHQlmc/UVhZBQgG8wzUEiIjIIxMuAlS10+4P8EUAzyCxfkC1qp5wKZuvrKguAwDsff+84SREROQXE+4TICKbAbyjqt8HsNve9l9EpExVN7kV0C8eXF4FvH0CB09cMR2FiIh8IpXHAY12AZDsBfBxgCOWLCyFDsZx/FKv6ShEROQTqYwOkJEbVLVTRG7YTqkLBAIID8ZxAXHTUYiIyCdSaQloFpGvJm8QkaeQAaMDRORxEWnq6uoyHWVMZaEg+oKsqYiIyBupFAHrAHxJRI6KyDYROQZgI4A17kRzjqruVNXGkpIS01HGdEtpPqQwjItX+k1HISIiH0hldECXPTrgS0jc/a9X1UWq2u1aOp9ZMjdRpHAhISIi8kLKCwiparOqPqequ90I5GerFicWEmo5woWEiIjIfal0DISILABQN3L7KKMGaBI+cVcV9Ecf4NDp9O67QERE2SGVeQKeR2I4YARAZ9KnFACLAAdML86D9MVwimsIEBGRB1JpCVgLoEZVj7sVhoCCOHCFwwSJiMgDqfQJOM4CwH2V+WHE8oKwLLYGEBGRu1JpCdgqIj8FsB2JRwLXqOrLjqbyseryaTjeeRUHj17GCrujIBERkRtSKQLWAygF8PSI7QpgkWOJfG75wjLs3n8Vew6eZxFARESumnARoKr1bgahhPvvrMQ39p/Eu8cvm45CRERZLuV5AshdK26dCY3GcewCFxIiIiJ3TbgIEJFiEXleRI6JyNv2tnUi8jX34vlPKBRAqD+O831DpqMQEVGWS6UlYDuAY/bUwQEAUNUXkOgrQA4qDQbQyzYaIiJyWSqXmlWq+uf2+5q0fYaDeQjA3OI8aGEIXT2DpqMQEVEWS6UIaBGR307eICKfRwYsJZxpFs8ugYjgtXfPmo5CRERZLJUi4AkAm0TkKIBqEdkGoAmJqYTJQfWLZwIA3jl00XASIiLKZqkMEewEUCsiDQDuQmL9gEZV5Wo3DvvEiiqg+Qg+PMUfLRERuSeVBYSOquoiVW0G0Jy0fQESLQIKYK2qdjsd0m9mzSyE9kVxso9TBxMRkXtSWjtARFaLyGYR+VzS9u0AWpEoDLY4ms7H8mOKjljMdAwiIspiqUwb3AlgIxIX+y+KSL2qfh3ASgCPqmq33V+AHFCeG8ZJWLAsC4EAxwsSEZHzUrm6PKqqn1LV51T1k0gsLQwA4CMA5y2cWQjJDeKjk+wXQERE7kj1ccByABCRu+x/i0d8jTgVzO/unD8dAPA6hwkSEZFLUikCGgG8JCJxJDoCHkeiL0BERLbZywxvdSGjL923tBIAsP9Yh+EkRESUrVIZItgKoHa0z4nIfwHQaU8jTA64+44KaMzC0fM9pqMQEVGWSmWI4E9xkyGAqvqco6kIOeEggn0xnJW46ShERJSlUnkc0IWkzoDkvuJAAFevW6aBiIjIOakMEdwMoFlEapA0WRAAqOrLjqYiAMDsabm4Eo+hbyCKgryw6ThERJRlUikCXgBwBYk1BJ5I2q4AFjkZihIWVxXhw7OdeOPgOTTcfYvpOERElGUm/DhAVetVtXaUt7QvAETkcRFp6urKrDH3dYsSCwm99SEXEiIiIuf5Yio6Vd2pqo0lJSWmo6TkE3dVAQDeP3nFcBIiIspGEy4CRKRYRJ4XkWMi8ra9bZ2IfM29eP42v6oY2h9De0ef6ShERJSFUmkJ2A7gmKrWDn+fPS/AejeCUULukIVLg1xIiIiInJdKx8BVqvop+/3kcWszHMxDI8zMDeG0cklhIiJyXiotAS0i8tvJG0Tk8wBanI1EyeaXFUDyQ2g/yzWaiIjIWakUAU8A2GQvF1wtItuQWEOg0ZVkBABYOi+xkNCr+7mQEBEROSuVIYKddn+ALwF4BokJg6pV9YRL2QjAPbeXAwBaj14ynISIiLJNKqMDVgOAqjar6nOq+oKqZtbA+wx0/7JZ0LjiyFkuJERERM5K5XFAk4hcFpHviMhy1xLRdQrywgj0RXHm6qDpKERElGVSeRxQC2AVEgsJvSQiR0Vks4gscCkb2aZB0G1xhAARETkrlSGCUNUIgKcBPC0iCwFsBNAGIOhCNrJVFebiiBXDUDSOnHD2/6j7BqL447/fh3Od/Z7uNxgI4OnfWIZb50/3dL9ERKakVAQAgIisAPAkgDUABMBzToei6y2qLMJHF7vw9gcX8OCKKtNxXGNZFv7k71vx3/d/DBTlGMlw9jtv4cfP/JKRfRMReW3CRYCIPIPEhX86ErMHPqGq+90KRr9wV20ZfnixC3vfP5+1RcDf7zyEzT86jGhJDgTAuiVV+L8/Wetphob/bzeO6fhfR0SULVJpCSgFsF5Vd7sVhkb34PIqYO9xvNeefQsJNb/9Mb76P1rRVRSChgWfKS/BX/zxvcjLTbmRaspqivJwyIrhSvcAphfneb5/IiKvTfhMq6pfdDMI3dyt80qgg3EcH+w1HcUxH0Y68MVv7UV7SIFcQX1ODp7fuBrl0/ONZXr4jkoc/vAsXmxuw/rP3WEsBxGRV1K63bJHAtSN3K6q33coD40iEAggZzCOi4ibjjJlFzr6sP4v96C1tx8IBrAgJtj6uw/htgXmO+P95mO1eP6DM/jZgTMsAojIF1LpE/A8ElMERwB0Jn1KAbAIcFlZKIRzGVwEDAzG8J+/uRc/OtkBKQihbEjxF/9+BVavmms62jXzq4oR7IniUA9XbSQif0ilJWAtgBpVPe5WGLq5edPzcX5wEOcu9WLWzELTcSbMsixs/ocD+O477dCiHORGLXz9wRr8h8/cZjraqObm5aAdcQwMxoz0SyAi8lIqMwYeZwFgzu1zSwAArx7InIWE/vHHR7Dky/+CFw4nMv/HW2fh8Ld/NW0LAAC4r3YGJCeIna+3m45CROS6VIqArSLyUxF5SkRWJ7+5lo6uWbUksZBQy5H0X0jolX2nUfef/gX/7yvHMJATwCfLpuGDv/gM/ug/rkQgkMp/Oe89uboGAPCvb500nISIyH2ptHeuR2KY4NMjtiuARY4lolE9tLwKuvN9HDmTvms2fdR+BV/85htoEwvIFawIhrD1D/9dRj2+qLutHLgaxcGr3aajEBG5LpUhgvVuBqGxlRTlQnpjOCXp1znwctcAvvSN1/FWdx8QDuCWqOA7X74Pd9bONB1tUioCAZwXhWVZad9yQUQ0FakOEVwBYBOAaiRGCfyZqr7rRjC6UaEFdCJ9FhIaisbxtW+/iX9puwgpDKNkyMJzT6zAp+6bZzralKycNx0/vtSNV1rP4JH69Bm9QETktAnf5ojIowBeBtAC4BkA+wC8LCKPuJSNRqgsyEEsP4hYzGwhYFkWnvvHA7jtKz/AD851IhxXfP3ehTj47c9mfAEAAJ99YD4A4PuvsXMgEWW3VFoCngGwMnmEgIg0A9gG9gnwRG1FISIdV3Hgo0uov73CSIYXm4/hv37/ffQXh4GA4AvV5fiTp+qzqtn8sXvmQr+3H++c6DAdhYjIVakUATfMEaCqrSIyw+FMdBPLF87Azzqu4o33znteBLzx7ln87nffwaXCIDQsWF1SiG/+0QOYVhD2NIcXAoEASmPAuTTsf0FE5KRUioAWEfk1Vf3n4Q0i8hQSjwfIAw8sq8Rz+9rxnVfa8JKHQ9j6Y3GczxEgV7BUgtj69Ycxt3KaZ/s3YemsIuzp7cd7xy5lbAdHIqLxpDpj4D4R+SISnQJXAVgIYKUbwehGyxbNQE7XEPpyAzhheTi1rQCzo4Jvf+nexBA6H/jlVXOx59+OYtvLERYBRJS1Uhki2AWgVkQ+j0QBsFVVX3AtGd0gEAjgo+/8mukYvvDZf7cQm3Ydxp4jF01HISJyzU2LAHs4IFT1QPJ2VX0JwEv21zwFoGXk1xBluoK8MAr64ziZwYs2ERGNZ6wu3S8AGG+CoOMAtjgXhyh9LC4rRKwojFPnr5qOQkTkirGKgBoAL471zaq6G+MXCkQZ6dFlVRARbGtuMx2FiMgVYw7uVtWJTKAuDmUhSiu//lgtNK54+f3MWbmRiCgVYxUBLePNBmjPIhhxNhJReiifno9wbxTHOvtNRyEicsVYRcCzAHaIyPzRPikiC5B4XLDB+VhE6WHhtDwMFIbQ0ztkOgoRkeNuWgSoajMSnf6Oi8g2EXlKRFbb/24D0AbgGVV92auwRF576LYKSCiAF3ezXwARZZ/x+gQ8C6AWief+TwNotv8VAPWq+pzrCYkM+s2GGgDAz1pPG05CROS8cScLUtUIgCc8yEKUdmrnlUJ6hvBBd9R0FCIix2XP0m9ELpmTE0ZPbgBDUU4cRETZhUUA0TjurZkByQ3iR3vaTUchInIUiwCicXz+EwsBADvf/NhwEiIiZ6WyiiCRL92ztALaG8WB3k7TUYiIHJWRLQEiUioia0SkUUS2i0i16UyUvQKBAMoRwCVRWJZlOg4RkWMysghAYrRCtao2AdgFYKPhPJTl6m4pgRSE8cbB86ajEBE5xpMiQETq7Dv2hhHbS0Vkg31Xv0FE6ibyeqraZM9hACQWOuJMLuSqx++dBwB46RXOkk1E2cP1PgFJF/7Rmuy3A1hvz0UAEdklImtVNZWHr9WqypYActX/df986PcP4u3jHaajEBE5xvUiwJ5+GCJy3dlTREqRuIAn31pFADQgsWZB401erynpNTao6lrnUxNdLxQKoDiqOIOY6ShERI4xOTqgHsDIO/5OAI8B2JF8sR+NiKwB0GS/3zBcbBC5ZWlFEfb2D+DwiSu4bcF003GIiKbMZMfAUgAj21YvAygb7xvtvgNbAOwWkTaM/qiByFG/VD8HAPC9ZnZBIaLsYHp0wLgX/NGoaquq1qjqSvvfm7Ya2MMIW0Sk5eLFi5NPSr73+UdqoDELrx+5YDoKEZEjTBYBnUi0BiSbgRtbB6bEHklQr6r15eXlTr40+cy0gjDy+mI40TtoOgoRkSNMFgEtuLEloBSJcf9EaenW0gJEp4VxoaPPdBQioikzVgTYwwBbRsz2Vw+AHfwoba2+swoSEPzvnx0zHYWIaMpcLwLsiYI2IHGB32i/P2wtgDXDkwUBWJfiHAFEnvqNx2qglmL3wbOmoxARTZkX8wS0AmgF8Owon+scbTtRupo1sxChq1F8hKjpKEREU2Z6dABRxplfkIv+/BCu9rEQIKLM5osiQEQeF5Gmrq4u01EoCzy0uBwSDuCf/43rCBBRZvNFEaCqO1W1saSkxHQUygJPPloDAPhRyynDSYiIpsbktMFEGen26jKgZwjv9/BxABFlNl+0BBA5rSoUQndYEItZpqMQEU0aiwCiSbh7YSbZ2QkAABSZSURBVBkkL4SfvfWx6ShERJPGIoBoEtY8vBAA8C972g0nISKaPBYBRJPwwPJZ0L4o9n3Mua2IKHOxCCCahEAggBmW4CLYJ4CIMheLAKJJWj6nBFIYxjsfnDcdhYhoUnxRBHCyIHLD4/fOAwC8+HNOGkREmckXRQAnCyI3fObB+dDBON5su2w6ChHRpHCyIKJJygkHMW3QwimJm45CRDQpvmgJIHLL7eXToEU5aDvFR01ElHlYBBBNwafqZgMAvtfcZjgJEVHqWAQQTcHa1TXQmIVXPuQIASLKPOwTQDQFJUW5yO2N4bjETEchIkoZWwKIpqimJB9DhWFc7howHYWIKCUsAoim6JGllZCgYNuuY6ajEBGlxBdFACcLIjf9ekMtVBW73j1rOgoRUUp8UQRwsiBy07xZRQj2RHH4cq/pKEREKfFFEUDktnn5OejLC2BgkB0EiShzsAggcsD9i2ZCcoL4P6+cMB2FiGjCWAQQOeDJ1TUAgB++/bHhJEREE8d5AogcsPzWmcDVKA72dJuOQkQ0YSwCiBxSGQziHCxYloVAgI1syc5c7EX72R5P9xkIAMsXzUReLk9zRDfDvw4ih9TPm44fXuzCyy2n0XD3LabjpI3X9p/BF/6hBZIb9HzfGrOQ0xfD7Lwc3Dm3BI+sqMKn7p2HaQVhz7MQpSMWAUQO+bWH5uOH3z+If36tnUVAkt/9u3eAHMEvl5cgFBTP9huNWzh6rgcnEcMJjaP9Qhf+9Wdd+P0ffYhwXwyVOWHcOacYDy+rwqfvn4eSolzPsmWqoWgc+49cxLvHOjA45O0S2qrAYDSOoZiFoWgcA9E4ojEr8XEsjqGYIhq3EI1ZiFqKWNxCNG4hZum1t7iliKv9vgKWKuIALCgsBSwBFPabACre/X9Ntv7uedj0/9R5si8WAUQOWV0/B/pPrWhpv2I6Str4ux8cwpVpIdSFc/Dtrz5oLEcsZmHPu2fxs5bT2Be5jHbEcQpxnL7cg5/8vAdPNx9BsDeKynAIS6qK8YmllfjlB+ajfHq+scymnDzXg7c+uICDbR04dq4bp6704/JgDH0BQPODkGBmPOpSVdhX+sSbArAUooCo/S8SveNFgaAIwgACkng/aP9rwtyZhZ7ti0UAkUMCgQCmx4Hz8PYOKV3FYhY2/+QwNCho+pPVRrOEQgE8vHIOHl4559o2y7Lw1vsX8NO3T6Gl7TKOI4YzsHC28ypefv0q/uurxxDojaI8GMRtlUV48I5KfOaB+Zhd7t0J2g19A1HsO3QRrUcu4fCpLpy41IvzVwfRbVmI5gQgeddfFlTjCEMxQwKoDIaxoLwQi+eUoDDf+8tHbjiI3Jwg8nOCyMsJIS83iPzcEPJzE+8X5iXez88LISfs/eOnTMQigMhBy6qK8WpPHw4cuYgVi8tNxzHq61vfRqw4B5+rKk3LO+pAIID7ls3CfctmXdtmWRYOHr2MH+79GG8fvYQIYjgPCxd6+vDqm8fxZ28eh/QMoSwQxK3lhXhwSSU+8+A8zK8qNngk17MsC5FT3Xj7wwt4/8QVHDvbg9Nd/egYiqM/AGhBCBL4xR2uxiwE1EIBBHMCIcwtzMets4txZ3UZ7r69AnMrpxk8GnKbL4oAEXkcwOO1tbWmo1CW+/SqW/Dqy0fw4s8jvi4CrnQPYNvhcwjFFc/96X2m40xYIBDAisXlN/zuPmjrwI/2nsSbH13EMcRwGRb29g1g7752PLevHbgaxXQIasoKUJTvbafDuKU42zWAC71D6FEL8bwgJOf6u2BVCzmqKJcgqkJhLKyYhtvnlWLlbeVYtmgG75p9zBdFgKruBLCzvr5+nekslN0++/ACbPzJh9jz0SXTUYz64jdehxSG8dX6+QiFMuMZ8ljuqCnDHTVl1207drITP3zjJPYeuYiPrsbQAUXL0BAwNOR5PrUsBGBhmgjKA0HcUlyAxXNKsKymDPfcUZmWLTGUHnxRBBB5JS83hMIBCx/7uF/A4RNX8GZ3H0oGLfzOmqWm47imdl4pvjKvFF9J2nahow89fVFPcwREML9qGuemoElhEUDksNtmFGLf0CBOnuvBvFlFpuN4bt1f7wFCAXzj399lOornKsoKUFE2/tcRpQuWjkQOa1heBRHB/951zHQUz/1oTztOhoGFlnCuBKIMwCKAyGFPPlYLjVv4tw/Om47iuQ3/tB+IWvjuVx4wHYWIJoCPA4gcNqMkDzm9MbRpzHQUT/319w7ianEYDxTko3Zeqek4RDQBbAkgckF1UR4GC0Po6hk0HcUTQ9E4/vrVNuBqFFsNzgxIRKlhEUDkgodvr4SEAti2u810FE/857/ZC6s4B19YNhtFhTmm4xDRBLEIIHLBk4/WAAB27T9jOIn7LnT04V/bLyOnawh/8lS96ThElAL2CSByQc3cEgS6h/AhvB0zbsJTf/4apCCEP/xELceqE2UY/sUSuWRObhhXcwMYimbvxEGthy/i3cEhzOiN4bc+vdh0HCJKEYsAIpfcVzsTkhvEztdOmI7imi/97ZuAAN9ad4/pKEQ0Cb4oAkTkcRFp6urqMh2FfOSJR6oBADvf/NhwEnds392G8wUBLAmGrluJj4gyhy+KAFXdqaqNJSUlpqOQj9TfXgHtjeLgmewrPi3Lwh++dBDaH8N3f/8h03GIaJLYMZDIRRUI4EJAYVlWVnWa+7P/sR8DxTn45PRpXG+eKINlz1mJKA2tnFcKKQjhtf1nTUdxTN9AFH/XchLSM4Rv/h6nBybKZCwCiFz02QcWAABeyqLOgV/+yz3Qohw03jMfeblsTCTKZCwCiFzUcPdc6EAM7xzvMB3FEe1nu/Hy+S7kd0ex8QsrTMchoiliGU/kolAogJKo4iyyYzGhp/5yD5AbxH/71duzqo8DkV/xr5jIZUsri4GiHHzQltmtAa/tP4OPrBhmDyg+Zw9/JKLMxiKAyGWfXjUHAPC9DF9M6Hf/7h1AFc9/+T7TUYjIISwCiFz2uUeqoVELez66YDrKpP39zkO4Mi2Eu/LzsPzWmabjEJFD2CeAyGUFeWHk98fQbjrIJMViFv7sx4ehIcELf7zadBwichBbAog8sHh6AWLTwjhzsdd0lJT9QdPbiBXn4Fery1E+Pd90HCJyEIsAIg88umw2JCD43q5jpqOkpKtnENsOnUOwewh//uV7TcchIoexCCDywK8/VgO1FC+/d850lJSs/8brQGEYv//IIuSEg6bjEJHD2CeAyAMVZQUIX43iKKKmo0zY4RNXsLerF8WDFr68dqnpOETkArYEEHlkQWEuBgpC6OkdMh1lQhr/5g0gFMA3fmul6ShE5BJfFAEi8riINHV1Zd+SrpQ5HrqtAhIK4KWfR0xHGddP955Ee0ixICZ47N5bTMchIpf4oghQ1Z2q2lhSUmI6CvnYbzTUAAB+su+04STj+9r/bAWiFl74yv2moxCRi9gngMgjt86fDukZwgfd6d0v4G+2vYee4jDuz8/DrfOnm45DRC7yRUsAUbqYHQ6hOzeAWMwyHWVUQ9E4/uqVY8DVKJ7//QdNxyEil7EIIPLQPdUzILlB/PiN9Jw/8KvfehNWcQ5+c2kVSopyTcchIpexCCDy0JqHFwIAfrD3pOEkN7rQ0YcfHL+EcPcQ/rRxlek4ROQBFgFEHrr3zkpobxStp9JvpMq6v3gNUhDCH3x6CQIBnhqI/IB/6UQeCgQCmIkALsGCZaVPv4ADRy7iwMAQyq7G8B8+c5vpOETkERYBRB5bMacEUhjGW++nz9LC6//2TUCAb62723QUIvIQiwAij/3KffMAADteSY9Jg156uQ3n8gSLAyHcv7zKdBwi8hCLACKP/fID86GDcbzZ1mE6CizLwh/sOAgMxvFdDgkk8h1OFkTksVAogOJBC2ckbjoKNv/DAQwU5+DR0mmYN6vIdBwi8hhbAogMuL1iGrQoBx+1XzGWoW8giu++0w7pGcK3fo/TAxP5EYsAIgM+tXIOAOB7u831C/hPf/UGtCgH6+6ej4K8sLEcRGQOiwAiA9auroHGLLx26LyR/Z8814Pd5zqR1z2Ep39rhZEMRGQe+wQQGVBUmIO8vhiOI2Zk/09943UgN4j/9itLOTEQkY/xr5/IkEUlBYhOC+NCR5+n+33j3bM4YsUwa0Dx+dU1nu6biNILWwKIDHlk6Sy8/94prP/LPbhtdrFn+/3Be2eB/CC2/s69nu2TiNITiwAiQ37jk7X4m33t2I8h7G+/5N2Oi8O4OzcXKxaXe7dPIkpLLAKIDJldXohdX3kI7Weverrf/LwQHlzBmQGJyCdFgIg8DuDx2tpa01GIrnPr/Om4df500zGIyKd80TFQVXeqamNJSYnpKERERGnDF0UAERER3YhFABERkU+xCCAiIvIpFgFEREQ+xSKAiIjIp1gEEBER+RSLACIiIp9iEUBERORTLAKIiIh8ikUAERGRT7EIICIi8ilRVdMZPCMiFwG0O/yyMwF4uA6sMTzO7MLjzC48zuzi9HHOV9VR1w73VRHgBhFpUdV60zncxuPMLjzO7MLjzC5eHicfBxAREfkUiwAiIiKfYhEwdU2mA3iEx5ldeJzZhceZXTw7TvYJICIi8im2BBAREfkUiwAiIiKfCpkOkIlEpBRAI4AIgGoAzaraajaV80SkDkCD/eEqAFtVtdlgJNeJSAOAUlXdYTqLG5J+pxEAZaqadc9Yk/4+O+1NkWz4f2v/7jZhxN9htp2PxjjOrDof3ew4R3yN6+cjFgGTsx3AelWNAICI7BKRtaraOc73ZZoGVX0WuHaiOS4ij2byCWYs9jFuBbDFdBY3DJ90VHWt/fE+ezxytv0+G4f/3wKAiGyxjzNj/z7tiwGQuMiPlDXno3GOM2vOR+Mc5/DXeHI+4uOAFNm/mOrhPzhbBL+oULNCUpUKALBPKC3IsuMc4QkAGXtnMQEvANiY9HFGnkAn4MkRH1/GGCfbTKCqzfbdYkfy9mw7H41xnFl1PrrZcY7gyfmIRUDq6vGLZsZhnQAeM5DFNfbFYe2IzdW48dizgl2ZZ20BYF8s6lQ1IiJ1IlKdiXeKExSxWzmqRaQawIwsLXYAno+y8v+wl+cjFgGpK8WN1dtlAGUGsrhqxPO4aiSO8UVzidxhXyBLR9xNZZt6JC6Oa2A/OxaRrYYzucJ+3BEB0IbE89aN43xLJuP5KMt4fT5iETA5WfcHNgFbkWg+zsbKuyFbOwImKcUvOo112ifUarsoyCp20/EuJO6Gq0Vkn+FIbuP5KLt4ej5iEZC6TiROqMlmYOxnOxlNRDYA2JiNTar2BSPrjmsUEQCdI06aEWRZs7F9h/ikqjbZz11rkGgBycrOnuD5KKuYOB9xdEDqWnBj5V2KxJ1H1rHvFK8NORKRuiz74ysDUC8iwx83ACgTEWTZ8LkIbrxYANn3TLUOwDsjtq1Dlo74AM9HPB9NEYuAFKlqp4i02B2rhp/Z1OP6XtdZwe6c0pn0B1eNxLFmzR/dyPG5IvIYgF1ZVgAM/7/dMeKkWY/EBTKbNCNxwU9uTq1HYhhd1uH5iOejqeLaAZOQbZNzjMb+A2sb5VMrs+1Yh9nNjJuQ+L1uzrZ+Avb/201I/F5rkDi5ZN2IiKRJZYZbOToy/XeZdEybkLj73zVizHxWnI9udpzZdj4a6/eZ9DWenI9YBBAREfkUOwYSERH5FIsAIiIin2IRQERE5FMsAoiIiHyKRQAREZFPsQggcoGIbBcRFZHGEdsbRGS0oU5O7bctaZlSz4nIBhG5IiJXbvL5BntImxEm9p/0M8nolQwpO7EIIHJPBFk4acvN2MXHk6o6XVWn3+TLtsPs8q+e7394/HeWL1BFGYpFAJF7dgDXpjr1gzqMs/ypXSAM/1w2uDmn/2ivn7x/r9jFUYuX+ySaKBYBRO65jMQUtptMB/HIDNMB0tRjyNK5/CnzsQggcpE953f1zZ7Tj3yGLyL7klsO7I8b7X+viMhWEam2P1YRGW1O/MdEZNfIz4tIqb39iv1v6Yj91Nl9GW56wRrxGvuGs9t33BsAbLCPadTn38PHa+fakvT1WyaT0b7bbxv59WO8fltS5lGP5WY/95v9TCbgWguJiKyxX7duCq9H5BgWAUTu24zJ9w0oBbBeVVcCWInEHPHb7Y+nA1gzSoFRqqqP2Z+vs+cgB4B99vdOR+LONLmpvBTAbnv72jHy7AOw1X6NtQC224sSbQTwLIAmVa0Z7/m3qq61v/5Z++uHfz6pZozY3z8dibUCtozz+uMeS9K+1ib/3KfQsa9eVVtFZI39KOIx3LjyH5ERLAKI3NcEoGEKF5GtwLWOZREA2+yPO5FYQW3k625P+vxGAOuHC4Xh1cjszmoji4cmVW2yv+8GdgtFZPiZup1nMxx63DGZjCOe729DYlW5iexrIseyPelzESTu6FM9pjoAERFpSNrX2mxcuIkyE5cSJnKZvdxrExJ3qZNpVk6+q+4c8XHHBL63DIkLWJmI7BvxWsnGe269asS+h1//yXG+b6JSzmhfZDchceeeSpE1kWMZ+XOfjAYkfkd1IlKa6asZUvZhEUDkjS1ILIU63oXW6THs1UhczDqRWGJ2rKb+8bQBWH+T13dCShnt5/+7ATxqN7c3YOJFltvHMuwxAFtUtdnuj9DKoYKUTvg4gMgD9ol/B27sG9AJ+w7Wvqt1YkKZtfbrVeMXrQ8vIvFIIrnTYePo3z66pE6OjUl5NyHRjD4Zl2E3sduvlWrG4Z/V8EV1ZPEw8vWvmeqx2J0KJ/K7qk9q+u8E0GFyMieikVgEEHlnM268yG8FsNVuAm9Aohf5ZJueAfuu3369fQB2JD1DX4lE/4Arkpi1cDKtDisBrLVnBNyOROe51klm3QGg3s7yZKoZ7f22ALiS9Agh+S77utd3+Fi2YJzOnnZLRfL8AFvxi8cDRGlBVNV0BiKijGK3VlQPzwZIlKnYEkBElLpVsGeEJMpkbAkgIiLyKbYEEBER+RSLACIiIp9iEUBERORTLAKIiIh8ikUAERGRT7EIICIi8ikWAURERD71/wNuBlh/uyi1ZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.semilogy(tol, label=r\"$\\|b_k - b^*\\|_2$\", color='navy')\n",
    "plt.xlabel(r\"Number of iteration, $k$\", fontsize=15)\n",
    "plt.ylabel(r\"Convergence rate\", fontsize=15)\n",
    "plt.xticks(fontsize = 15)\n",
    "plt.yticks(fontsize = 15)\n",
    "plt.legend(loc=\"best\", fontsize=15)\n",
    "plt.plot(tol)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}