S21_sorter_state.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "Vackr_erxAp1",
   "metadata": {},
   "source": [
    "# Modeling the World State"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wcN7yQNq-aLh",
   "metadata": {
    "colab_type": "text",
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "no-tex"
    ]
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/gtbook/robotics/blob/main/S21_sorter_state.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "SfD_sWF1YbYG",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "9699f781-e0d4-4fad-a0a0-7ea90ff91392",
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install -U -q gtbook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "zxdBVtluL_Mh",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gtsam\n",
    "\n",
    "import plotly.express as px\n",
    "try:\n",
    "    import google.colab\n",
    "except:\n",
    "    import plotly.io as pio\n",
    "    pio.renderers.default = \"png\"\n",
    "\n",
    "import gtbook\n",
    "from gtbook.discrete import Variables\n",
    "from gtbook.display import show\n",
    "\n",
    "def pretty(obj): \n",
    "    return gtbook.display.pretty(obj, VARIABLES)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "iaryKoy0AdjG",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "```{index} state; discrete state\n",
    "```\n",
    "\n",
    "> The physical properties of a piece of trash comprise all of the information needed by the robot."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3e0174a",
   "metadata": {
    "tags": [
     "no-pdf"
    ]
   },
   "source": [
    "<img src=\"Figures2/S21-Trash_sorting_robot_with_gripper-09.jpg\" alt=\"Splash image with robot approaching a trash bin\" width=\"40%\" align=center style=\"vertical-align:middle;margin:10px 0px\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2-E-_s0tCYOi",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For our simple trash sorting robot, the only thing that matters at a given moment is the category\n",
    "of the item of trash on the conveyor belt. Remember that items of trash are presented individually to the robot,\n",
    "so there is no clutter, and no circumstance in which multiple pieces of trash are simultaneously in the workspace.\n",
    "Therefore, it is natural to define the world state explicitly in terms of the category of the current\n",
    "item of trash. We consider five possible categories:\n",
    "- cardboard\n",
    "- paper\n",
    "- can\n",
    "- scrap metal\n",
    "- bottle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eeoyuiaQX1p",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "For simplicity, we assume here that there will never be a piece of trash that does not belong to one of these categories.  We do not, however, assume that the category of an item can be reliably determined with 100% accuracy. \n",
    "Instead, we use probability theory to quantify the uncertainty associated to an object's categorization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "J7aya4jqdDFQ",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Using Probability to Model Uncertainty\n",
    "\n",
    "> Probability theory provides a rigorous methodology for reasoning about uncertainty."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "uowcFQT9MCy5",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will use probability theory to model uncertainty.\n",
    "While a comprehensive coverage of probability theory is beyond the scope of this book,\n",
    "we introduce key concepts and methods throughout the text, as needed,\n",
    "to deal with various kinds of uncertainty that occur in robotics applications.\n",
    "Rigorous introductions can be found in many textbooks, including the book  [\"Introduction to Probability for Data Science\"](https://probability4datascience.com/index.html) {cite:p}`Chan23book_prob4ds` (which is available online)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "JUVRQKCigGOW",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```{index} sample space, outcome\n",
    "```\n",
    "The starting point for reasoning with uncertainty is to define the set of outcomes that might occur.\n",
    "The set of all possible outcomes is called the **sample space**, often denoted by $\\Omega.$\n",
    "In our example, when an item of trash arrives on the conveyor belt,\n",
    "there are five possible outcomes,\n",
    "\\begin{equation}\n",
    "\\Omega = \\{ \\rm{cardboard, paper, can, scrap \\; metal, bottle}\\}.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "D4bXMY8c_7Bn",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "```{index} event, probability distribution\n",
    "```\n",
    "## Probability Distributions\n",
    "A subset of the sample space $\\Omega$ is called an **event**.  A **probability distribution**, $P$, assigns a probability $0 \\leq P(A) \\leq 1$ to each event $A \\subseteq \\Omega$, with $P(\\emptyset) = 0$ and $P(\\Omega)=1$. \n",
    "In addition, for disjoint events, $A_i \\cap A_j = \\emptyset$, we have\n",
    "$P(A_i \\cup A_j) = P(A_i) + P(A_j)$.\n",
    "Using this property, it is a simple matter to compute the probability for any $A \\subseteq \\Omega$\n",
    "if we are provided with the probabilities of the individual outcomes.\n",
    "Further, since $P(\\Omega)=1$, it follows immediately that \n",
    "\\begin{equation}\n",
    "P(\\Omega) = \\sum_{\\omega \\in \\Omega} P(\\{\\omega\\}) = 1\n",
    "\\end{equation}\n",
    "i.e., that the probabilities of the individual outcomes sum to unity.\n",
    "As a slight abuse of notation, for singleton events, we will often write $P(\\omega)$ rather than $P(\\{\\omega\\})$\n",
    "to simplify notation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eoKrdorUJOEL",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "In robotics applications, the probability assigned to an outcome reflects our certainty in that outcome.\n",
    "These probabilities can change based on the arrival of new evidence.\n",
    "In robotics, this can occur when the robot acts in the world, or based on sensor data.\n",
    "How evidence affects the propagation of probability values is a recurring topic in this book."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "IQOJDtsbjpPg",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "```{index} prior probability distribution\n",
    "```\n",
    "\n",
    "## Prior Probability Distributions\n",
    "> A *prior* describes our beliefs before any sensor data is obtained.\n",
    "\n",
    "Once we have enumerated the set of possible outcomes, we confront a fundamental question: *Where do the\n",
    "probability values come from?* In this section we explicitly consider the notion of prior knowledge that is available in a particular application. High-quality \"priors\" can make a big difference in performance, especially when measurements are few or unreliable.\n",
    "\n",
    "In some cases, we merely assume that all outcomes are equally likely, for example, when rolling a die or tossing coin.\n",
    "In such cases, the probability of any outcome is merely $P(\\omega) = 1/N$ for each $\\omega \\in \\Omega$, where $N =| \\Omega |$.\n",
    "This leads to $P(\\mathrm{heads}) = P(\\mathrm{tails}) = 0.5$ when tossing a fair coin, \n",
    "where $\\Omega = \\{ \\mathrm{heads, tails} \\}$.\n",
    "\n",
    "In other cases, we can estimate probabilities using data.\n",
    "Suppose, for example, that the owner of the trash-sorting facility has told us (or we have kept statistics over time) that for every 1000 pieces of trash, the observed category counts are approximately as follows:\n",
    "- cardboard: 200\n",
    "- paper: 300\n",
    "- can: 250\n",
    "- scrap metal: 200\n",
    "- bottle: 50"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "J4OeG1SU0QFD",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It is common to assume that outcomes occur in proportion to their probability. There are a number of\n",
    "technical conditions that underlie this assumption, such as the condition that outcomes are independent,\n",
    "but we will not address these here. Thus, from the above observed frequencies, we might estimate that the probability of seeing a piece of cardboard in the work cell is given by\n",
    "\\begin{equation}\n",
    "P(\\mathrm{cardboard}) \\approx 200/1000 = 0.2\n",
    "\\end{equation}\n",
    "Using the same logic, we can do this for all categories, yielding:\n",
    "\n",
    "| *Category (C)* | *P(C)* |\n",
    "|----------------|:------:|\n",
    "| cardboard      |   0.20 |\n",
    "| paper          |   0.30 |\n",
    "| can            |   0.25 |\n",
    "| scrap metal    |   0.20 |\n",
    "| bottle         |   0.05 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cGAke7fQFxLJ",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "```{index} prior\n",
    "```\n",
    "We call this type of probabilistic knowledge about the state of the world, in the absence of any other information, a **prior**, because it represents our belief *before* any evidence (e.g., sensor data) has been acquired.\n",
    "\n",
    "<!-- Here, a word about notation is in order.  We use upper case letters to denote random quantities. For example, we use $C$ to denote the random outcome of a particular\n",
    "category of trash arriving to the sorting station.  When you see an uppercase letter, it indicates that we do not know what value this outcome will take, but\n",
    "that we do know a probability distribution for the values that could occur.\n",
    "We use lower case letters, possibly subscripted, to denote specific outcomes. For example, we might define $c_1 =$ cardboard,\n",
    "$c_2 =$ paper, etc.  Thus, we might write $P(C = c_1)$ or merely $P(c_1)$ to denote the prior probability that an item of trash is cardboard. -->"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "_GZ72xtYko-4",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Probability Distributions in Python\n",
    "> We represent probability distributions using the `DiscreteDistribution` class in GTSAM.\n",
    "\n",
    "```{index} GTSAM\n",
    "```\n",
    "\n",
    "```{index} pair: structure from motion; SfM\n",
    "```\n",
    "\n",
    "```{index} pair: simultaneous localization and mapping; SLAM\n",
    "```\n",
    "\n",
    "The GTSAM library (GTSAM stands for “Georgia Tech Smoothing and\n",
    "Mapping”) is a BSD-licensed C++ library based on factor graphs which we make heavy use of in the code examples.\n",
    "It provides state of the art solutions to important problems in robotics,\n",
    "such as the Simultaneous Localization and Mapping (SLAM) and Structure from Motion (SfM) problems,\n",
    "but can also be used to model and solve both\n",
    "simpler and more complex estimation problems. More information is\n",
    "available at <http://gtsam.org>.\n",
    "\n",
    "GTSAM also provides both a MATLAB and a python interface,\n",
    "enabling rapid prototype development, visualization, and user interaction.\n",
    "The python library can be imported directly into a Google colab via\n",
    "`import gtsam`. A large subset of the GTSAM functionality can be\n",
    "accessed through wrapped classes from within python. \n",
    "To not interrupt the flow of the book too much, we do not always fully explain the code throughout the text, but rather include a \"GTSAM 101\" subsection at the end that elaborates on the types and functions we used. \n",
    "\n",
    "The code below illustrates the use of GTSAM. First we create a `Variables` data structure that will be used to obtain more informative output from other code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "A6pbb7RuPpfX",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "VARIABLES = Variables()\n",
    "categories = [\"cardboard\", \"paper\", \"can\", \"scrap metal\", \"bottle\"]\n",
    "Category = VARIABLES.discrete(\"Category\", categories)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "TQVU6l0xWrkf",
   "metadata": {},
   "source": [
    "Conceptually, the `Variables` class keeps track of the names of variables and what values each variable can take on. For example in the above, we need the variable `Category`, and it can take on the values `cardboard`, `paper`, `can`, `scrap metal`, and `bottle`. We do this so that later when we print, it can show us a nicely rendered outputs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vDiuY2aG0D4k",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can now create a prior probability distribution $P(Category)$ on the category using a `DiscreteDistribution` constructor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eqNBywYH4MHq",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "category_prior = gtsam.DiscreteDistribution(Category, \"200/300/250/200/50\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7pw0jLvfKvoT",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The constructor automatically normalizes the numbers given to it to a proper probability distribution, i.e., it makes the probabilities sum to one. It is rendered in notebook as a table below, where we can verify this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "xv15iFdGujC1",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 175
    },
    "outputId": "c8706115-a1e8-430e-f99b-e31d82e37350",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<p>  <i>P(Category):</i></p>\n",
       "<div>\n",
       "<table class='DecisionTreeFactor'>\n",
       "  <thead>\n",
       "    <tr><th>Category</th><th>value</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>cardboard</th><td>0.2</td></tr>\n",
       "    <tr><th>paper</th><td>0.3</td></tr>\n",
       "    <tr><th>can</th><td>0.25</td></tr>\n",
       "    <tr><th>scrap metal</th><td>0.2</td></tr>\n",
       "    <tr><th>bottle</th><td>0.05</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<gtbook.display.pretty at 0x10556d460>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pretty(category_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hd-wCIHE9pLx",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can evaluate the prior for any category value, e.g., \"can\", but that function takes integer indices, not strings. Hence, we use the built-in python function `index` to obtain that integer (2 in this case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "iYf7Ryv6w6Lp",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "a9817da7-938b-49a5-cfd1-9c1debc9f65e",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P('can') = 0.25\n"
     ]
    }
   ],
   "source": [
    "index = categories.index('can')  # we still have to use an integer value\n",
    "P_can = category_prior(index)\n",
    "print(f\"P('can') = {P_can}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "C0CivmClcs4c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can also recover all values in the probability distribution at once, using the `pmf` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "IxEaCuBC8Q9x",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "786197c4-d1c7-4fd1-87e8-9a4326a10996",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.2, 0.3, 0.25, 0.2, 0.05]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PMF = category_prior.pmf()\n",
    "PMF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "FQJWZ6ePO5Pm",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "```{index} pair: probability mass function; PMF\n",
    "```\n",
    "\n",
    "Here \"PMF\" is short for \"probability mass function\", which we define more precisely in Section 2.2. Note that the ordering of the array was fixed when we defined `categories` above.  It is the programmer's responsibility\n",
    "to maintain consistency when using arrays to store values associated to a collection of variables.\n",
    "\n",
    "We can display probability distributions in various ways, including as a bar graph, as shown in Figure [2.1](#fig:discrete-distribution)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4_GcI74bqHSr",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 542
    },
    "outputId": "495376af-aae7-4c01-d35d-4a958aa67866",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4Xu2dd5wV1fn/n6UKShNEjV1iLNEYkxijJrGL3WBLwF6ixF5AVCT2gohEsaCgYokKGruARGzEHkWjiTHR+NVYEAQFFJC2v9eMv92wsLBz9jPPzOzc9/0POM9zzrw/M2ffO3fupaq6urraeEEAAhCAAAQgAAEIQKCkBKoQ3pImy2FBAAIQgAAEIAABCMQEEF5OBAhAAAIQgAAEIACBUhNAeEsdLwcHAQhAAAIQgAAEIIDwcg5AAAIQgAAEIAABCJSaAMJb6ng5OAhAAAIQgAAEIAABhJdzAAIQgAAEIAABCECg1AQQ3lLHy8FBAAIQgAAEIAABCCC8nAMQgAAEIAABCEAAAqUmgPCWOl4ODgIQgAAEIAABCEAA4eUcgAAEIAABCEAAAhAoNQGEt9TxcnAQgAAEIAABCEAAAggv5wAEIAABCEAAAhCAQKkJILyljpeDgwAEIAABCEAAAhBAeDkHIAABCEAAAhCAAARKTQDhLXW8HBwEIAABCEAAAhCAAMLLOQABCEAAAhCAAAQgUGoCCG+p4+XgIAABCEAAAhCAAAQQXs4BCEAAAhCAAAQgAIFSE0B4Sx0vBwcBCEAAAhCAAAQggPByDkAAAhCAAAQgAAEIlJoAwlvqeDk4CEAAAhCAAAQgAAGEl3MAAhCAAAQgAAEIQKDUBBDeUsfLwUEAAhCAAAQgAAEIILycAxCAAAQgAAEIQAACpSaA8JY6Xg4OAhCAAAQgAAEIQADh5RyAAAQgAAEIQAACECg1AYS31PFycBCAAAQgAAEIQAACCC/nAAQgAAEIQAACEIBAqQkgvKWOl4ODAAQgAAEIQAACEEB4OQcgAAEIQAACEIAABEpNAOEtdbwcHAQgAAEIQAACEIAAwss5AAEIQAACEIAABCBQagIIb6nj5eAgAAEIQAACEIAABBBezgEIQAACEIAABCAAgVITQHhLHS8HBwEIQAACEIAABCCA8HIOQAACEIAABCAAAQiUmgDCW+p4OTgIQAACEIAABCAAAYSXcwACEIAABCAAAQhAoNQEEN5Sx8vBQQACEIAABCAAAQggvJwDEIAABCAAAQhAAAKlJoDwljpeDg4CEIAABCAAAQhAAOHlHIAABCAAAQhAAAIQKDUBhLfU8XJwEIAABCAAAQhAAAIIL+cABCAAAQhAAAIQgECpCSC8pY6Xg4MABCAAAQhAAAIQQHg5ByAAAQhAAAIQgAAESk0A4S11vBwcBCAAAQhAAAIQgADCyzkAAQhAAAIQgAAEIFBqAghvqePl4CAAAQhAAAIQgAAEEF7OAQhAAAIQgAAEIACBUhNAeEsdLwcHAQhAAAIQgAAEIIDwcg5AAAIQgAAEIAABCJSaAMJb6ng5OAhAAAIQgAAEIAABhJdzAAIQgAAEIAABCECg1AQQ3lLHy8FBAAIQgAAEIAABCCC8nAMQgAAEIAABCEAAAqUmgPCWOl4ODgIQgAAEIAABCEAA4eUcgAAEIAABCEAAAhAoNQGEt9TxcnAQgAAEIAABCEAAAgiveA58Mm2O2IFyCEAAAhCAAAQgsHwC3+ncBkQCAYRXgBeVIrwiQMohAAEIQAACEGiQAMLbIKLlDkB4NX4Ir8iPcghAAAIQgAAEGiaA8DbMaHkjEF6NH8Ir8qMcAhCAAAQgAIGGCSC8DTNCeDVGy63mkQZHuLSGAAQgAAEIQCAmgPBqJwJ3eDV+3OEV+VEOAQhAAAIQgEDDBBDehhlxh1djxB1eR360hgAEIAABCECgYQIIb8OMEF6NEcLryI/WEIAABCAAAQg0TADhbZgRwqsxQngd+dEaAhCAAAQgAIGGCSC8DTNCeDVGCK8jP1pDAAIQgAAEINAwAYS3YUYIr8YI4XXkR2sIQAACEIAABBomgPA2zAjh1RghvI78aA0BCEAAAhCAQMMEEN6GGSG89RCYN2++fTHjK+vapaNVVVUtl+KChQvt8+kzrHpRtXXt0smaN29WO57v4dVOQKohAAEIQAACEGiYAMLbMCOEdzEC1dXVdsPtD9t1tz4Q/+3KHdvZtZeeaptv0q1eTqMeetIuHHJ77b+tukonu+bik23TDdeL/w7h1U5AqiEAAQhAAAIQaJgAwtswI4R3MQKT3vq3HXLiJXbH0HNss43Wt2tuvt8em/CCPTHqKmvWbOk7vY+Mf946dljJfvyDDS2609vnguttwYKFdsuQfgivdu5RDQEIQAACEIBAQgIIb0JQyxhWcf/T2uBho+3tdz+wEVf2jZFM+fxL2+GAU+2+4RfYxhus0yDNPhfeYIsWVdtV5x+P8DZIiwEQgAAEIAABCKRBAOHVKFac8EbC2qnDStb/lENryX1/+yPs+stOs+223nyZNB8e/5w9+ZdJ9q///NeuOv8E2+i7ayO82rlHNQQgAAEIQAACCQkgvAlBcYf3WwLH9r3SNuy2tp3R+6BaJFvu3tvO73OE7bnTz5ZJ8w/D77NX//Yvm/L5F3bRmUfbT7fYKB47a84CLQGqi0mg2mxhdXUx19bEV9W8nkeHmvghsXwIQAAC7gTatWnhPkeZJ6jIO7zRB9XOOfmQ2lyT3OGtGXzjHY/YnX8abxMfHPqt8M6eX+bzo2KP7bOpi2zsE2YzZy7/GzwqFlAjD7x9e7Pdd6myVbs0sgFlEIAABCqUQLu2LSv0yNM57IoT3ugZ3nfe+9BuGtQnJhj6DO/4Z/5qp513rb0x4WZr0bw539KQznlYuC7Tp1fZHXdV2bTp//sKusItsgkuqHPnajus1yLr1Im7500wPpYMAQjkSIBHGjT4FSe8//uWhv622cbr29Uj7rMxE16s/ZaGkaPH2YSJr8Xf4hC9rh/5oG37081sw25r2bQvZlr0DHCb1q34lgbtvCt8NcLrExHC68OVrhCAQPkJILxaxhUnvNH38F576wM27PaHY3Jt26xgNw06w7bYdIP4z4Ouv8dGP/K0vTJ2WPzn/pePsAfH/aWWcjTu8v7H2pqrrxL/Hd/Dq52ARa1GeH2SQXh9uNIVAhAoPwGEV8u44oS3Btfcb+bZ9C9m2mpdO9f7/buLY43+V7Yp0760ldq2ib+Td/EXwqudgEWtRnh9kkF4fbjSFQIQKD8BhFfLuGKFV8P2v2qENy2SxeqD8PrkgfD6cKUrBCBQfgIIr5Yxwqvx45EGkV9RyxFen2QQXh+udIUABMpPAOHVMkZ4NX4Ir8ivqOUIr08yCK8PV7pCAALlJ4DwahkjvBo/hFfkV9RyhNcnGYTXhytdIQCB8hNAeLWMEV6NH8Ir8itqOcLrkwzC68OVrhCAQPkJILxaxgivxg/hFfkVtRzh9UkG4fXhSlcIQKD8BBBeLWOEV+OH8Ir8ilqO8Pokg/D6cKUrBCBQfgIIr5YxwqvxQ3hFfkUtR3h9kkF4fbjSFQIQKD8BhFfLGOHV+CG8Ir+iliO8PskgvD5c6QoBCJSfAMKrZYzwavwQXpFfUcsRXp9kEF4frnSFAATKTwDh1TJGeDV+CK/Ir6jlCK9PMgivD1e6QgAC5SeA8GoZI7waP4RX5FfUcoTXJxmE14crXSEAgfITQHi1jBFejR/CK/IrajnC65MMwuvDla4QgED5CSC8WsYIr8YP4RX5FbUc4fVJBuH14UpXCECg/AQQXi1jhFfjh/CK/IpajvD6JIPw+nClKwQgUH4CCK+WMcKr8UN4RX5FLUd4fZJBeH240hUCECg/AYRXyxjh1fghvCK/opYjvD7JILw+XOkKAQiUnwDCq2WM8Gr8EF6RX1HLEV6fZBBeH650hQAEyk8A4dUyRng1fgivyK+o5QivTzIIrw9XukIAAuUngPBqGSO8Gj+EV+RX1HKE1ycZhNeHK10hAIHyE0B4tYwRXo0fwivyK2o5wuuTDMLrw5WuEIBA+QkgvFrGCK/GD+EV+RW1HOH1SQbh9eFKVwhAoPwEEF4tY4RX44fwivyKWo7w+iSD8PpwpSsEIFB+AgivljHCq/FDeEV+RS1HeH2SQXh9uNIVAhAoPwGEV8sY4dX4Ibwiv6KWI7w+ySC8PlzpCgEIlJ8AwqtljPBq/BBekV9RyxFen2QQXh+udIUABMpPAOHVMkZ4NX4Ir8ivqOUIr08yCK8PV7pCAALlJ4DwahkjvBo/hFfkV9RyhNcnGYTXhytdIQCB8hNAeLWMEV6NH8Ir8itqOcLrkwzC68OVrhCAQPkJILxaxgivxg/hFfkVtRzh9UkG4fXhSlcIQKD8BBBeLWOEV+OH8Ir8ilqO8Pokg/D6cKUrBCBQfgIIr5YxwqvxQ3hFfkUtR3h9kkF4fbjSFQIQKD8BhFfLGOHV+CG8Ir+iliO8PskgvD5c6QoBCJSfAMKrZYzwavwQXpFfUcsRXp9kEF4frnSFAATKTwDh1TJGeDV+CK/Ir6jlCK9PMgivD1e6QgAC5SeA8GoZI7waP4RX5FfUcoTXJxmE14crXSEAgfITQHi1jBFejR/CK/IrajnC65MMwuvDla4QgED5CSC8WsYIr8YP4RX5FbUc4fVJBuH14UpXCECg/AQQXi1jhFfjh/CK/IpajvD6JIPw+nClKwQgUH4CCK+WMcKr8UN4RX5FLUd4fZJBeH240hUCECg/AYRXyxjh1fghvCK/opYjvD7JILw+XOkKAQiUnwDCq2WM8Gr8EF6RX1HLEV6fZBBeH650hQAEyk8A4dUyRng1fgivyK+o5QivTzIIrw9XukIAAuUngPBqGSO8Gj+EV+RX1HKE1ycZhNeHK10hAIHyE0B4tYwRXo0fwivyK2o5wuuTDMLrw5WuEIBA+QkgvFrGCK/GD+EV+RW1HOH1SQbh9eFKVwhAoPwEEF4t49IK76JF1TZl2hfWZeUO1qJ58wYpzZj1tX3zzXzr2qVjg2MXH/DJtDlB4xncNAggvD45Ibw+XOkKAQiUnwDCq2VcSuF95oU3rM+FN9jsOXNjOuedfrgdtM8O9ZL6fPoMO+zkS+2Djz6L/73bOt+x3x68l+296zbxnydMfM1OHnDNUrWvjR9urVu15A6vdv4Vthrh9YkG4fXhSlcIQKD8BBBeLePSCe+cufPslz1OthOP6mEH77ezPf3863bKgKH2+N2DbM3VV1mK1pTPv7QHx020fbpvayu2WcHuuG+83TpqnD37wDXWZoVW9sTEV+3sS4fbfcMvqFO79hpdraqqCuHVzr/CViO8PtEgvD5c6QoBCJSfAMKrZVw64Y3u7h5/9hCbNH64tWrVMqazxyH9Yvk9eL9dGqT10adTrXvPvnbH0HPsR5t9LxbeCwaPtIkPDq23lkcaGkTaJAcgvD6xIbw+XOkKAQiUnwDCq2VcOuEd/cjTNnLUWBtz58BaMif1v9rWXWt1O6P3QQ3SemDsRDt34M2x4K7csV0svNEd4n27b2utW7eyn2y+oXXffsva54IR3gaRNskBCK9PbAivD1e6QgAC5SeA8GoZl054R9z1mI176uU6jyBEz/Ou1LaNnd/niOXS+vf7H1mv4y+2ww/sHj8SEb3e/Of79vjTL1uHdivaJ59Ns9EPP2W9euxk/U85NP73WXMWaAlQXUgCk6csshEjq23a9GaFXF9TXVQkvMccXmWrda1K9xCqzRZWV6fbk24xgebNUs4KrhCAQKMItGvTolF1FH1LoHTC29g7vB9P/twOPekS2/KHG9mlZ/3WmjevX3TuH/OsDbjiFntjws3xXd5Zs+dzLpWQwOQp1TbiNoQ37WhrhXfpx+mlqT6busjGPmE2cyZyJoFcorh9e7Pdd6myVbuk2ZVeEIBAYwi0a/vtY5q8GkegdMJb8wzv638eYS1bfvvbUPRM7mEH7rrMZ3jfff9jO/K0y23Hn//IBpx22HK/xmziS29a736D7dXHb7IVWrfiQ2uNO+8KX8UjDT4ReT3SQF5NKy+f1dIVAuUmwCMNWr6lE97Zc76xLXc/zvqd0NN61fMtDbO+mm1HnjbQju65h+2+41b2znv/tf2OHmB77vQzO+no/axZs2/v7LZt09o6dWhndz0wwTbstpZt8r11bcasr6zvhcOsZYvmdsuQfvE4nuHVTsCiViNQPskgvD5cvbp65eW1XvpCoMwEEF4t3dIJb4TjyecmWfRBtZrXuaceaj1/tVP8xxkzv7Zt9jnBav5u7JMvxd/Zu+Qr+h7ey8851q66cbTdfPeY2n/+wSbdbNCA3rVfcYbwaidgUasRXp9kvASKvJpWXj6rpSsEyk0A4dXyLaXwRkgWLlxkk6dOt66dO9Y+2tBYVHO/mWdTp31p7VZsax07rFSnDcLbWKrFrkOgfPJBeH24enX1ystrvfSFQJkJILxauqUVXg1L8mqENzmrpjQS4fVJy0ugyKtp5eWzWrpCoNwEEF4tX4RX48czvCK/opYjUD7JILw+XL26euXltV76QqDMBBBeLV2EV+OH8Ir8ilqO8Pok4yVQ5NW08vJZLV0hUG4CCK+WL8Kr8UN4RX5FLUegfJJBeH24enX1ystrvfSFQJkJILxaugivxg/hFfkVtRzh9UnGS6DIq2nl5bNaukKg3AQQXi1fhFfjh/CK/IpajkD5JIPw+nD16uqVl9d66QuBMhNAeLV0EV6NH8Ir8itqOcLrk4yXQJFX08rLZ7V0hUC5CSC8Wr4Ir8YP4RX5FbUcgfJJBuH14erV1Ssvr/XSFwJlJoDwaukivBo/hFfkV9RyhNcnGS+BIq+mlZfPaukKgXITQHi1fBFejR/CK/IrajkC5ZMMwuvD1aurV15e66UvBMpMAOHV0kV4NX4Ir8ivqOUIr08yXgJFXk0rL5/V0hUC5SaA8Gr5IrwaP4RX5FfUcgTKJxmE14erV1evvLzWS18IlJkAwquli/Bq/BBekV9RyxFen2S8BIq8mlZePqulKwTKTQDh1fJFeDV+CK/Ir6jlCJRPMgivD1evrl55ea2XvhAoMwGEV0sX4dX4Ibwiv6KWI7w+yXgJFHk1rbx8VktXCJSbAMKr5YvwavwQXpFfUcsRKJ9kEF4frl5dvfLyWi99IVBmAgivli7Cq/FDeEV+RS1HeH2S8RIo8mpaefmslq4QKDcBhFfLF+HV+CG8Ir+iliNQPskgvD5cvbp65eW1XvpCoMwEEF4tXYRX44fwivyKWo7w+iTjJVDk1bTy8lktXSFQbgIIr5YvwqvxQ3hFfkUtR6B8kkF4fbh6dfXKy2u99IVAmQkgvFq6CK/GD+EV+RW1HOH1ScZLoMiraeXls1q6QqDcBBBeLV+EV+OH8Ir8ilqOQPkkg/D6cPXq6pWX13rpC4EyE0B4tXQRXo0fwivyK2o5wuuTjJdAkVfTystntXSFQLkJILxavgivxg/hFfkVtRyB8kkG4fXh6tXVKy+v9dIXAmUmgPBq6SK8Gj+EV+RX1HKE1ycZL4Eir6aVl89q6QqBchNAeLV8EV6NH8Ir8itqOQLlkwzC68PVq6tXXl7rpS8EykwA4dXSRXg1fgivyK+o5QivTzJeAkVeTSsvn9XSFQLlJoDwavkivBo/hFfkV9RyBMonGYTXh6tXV6+8vNZLXwiUmQDCq6WL8Gr8EF6RX1HLEV6fZLwEiryaVl4+q6UrBMpNAOHV8kV4NX4Ir8ivqOUIlE8yCK8PV6+uXnl5rZe+ECgzAYRXSxfh1fghvCK/opYjvD7JeAkUeTWtvHxWS1cIlJsAwqvli/Bq/BBekV9RyxEon2QQXh+uXl298vJaL30hUGYCCK+WLsKr8UN4RX5FLUd4fZLxEijyalp5+ayWrhAoNwGEV8sX4dX4Ibwiv6KWI1A+ySC8Ply9unrl5bVe+kKgzAQQXi1dhFfjh/CK/IpajvD6JOMlUOTVtPLyWS1dIVBuAgivli/Cq/FDeEV+RS1HoHySQXh9uHp19crLa730hUCZCSC8WroIr8YP4RX5FbUc4fVJxkugyKtp5eWzWrpCoNwEEF4tX4RX44fwivyKWo5A+SSD8Ppw9erqlZfXeukLgTITQHi1dBFejR/CK/IrajnC65OMl0CRV9PKy2e1dIVAuQkgvFq+CK/GD+EV+RW1HIHySQbh9eHq1dUrL6/10hcCZSaA8GrpIrwaP4RX5FfUcoTXJxkvgSKvppWXz2rpCoFyE0B4tXwRXo0fwivyK2o5AuWTDMLrw9Wrq1deXuulLwTKTADh1dJFeDV+CK/Ir6jlCK9PMl4CRV5NKy+f1dIVAuUmgPBq+SK8Gj+EV+RX1HIEyicZhNeHq1dXr7y81ktfCJSZAMKrpYvwavwQXpFfUcsRXp9kvASKvJpWXj6rpSsEyk0A4dXyRXg1fgivyK+o5QiUTzIIrw9Xr65eeXmtl74QKDMBhFdLt2KFd968+fbFjK+sa5eOVlVVtVyK8xcstM+nfWkrd2pvrVu1rDP2k2lztASoLiQBhNcnFi+BIq+mlZfPaukKgXITQHi1fCtOeKurq+2G2x+26259ICa3csd2du2lp9rmm3Srl+TwPz5qfxh+X+2/dd9+Szvv9COsQ/sV479DeLUTsKjVCJRPMgivD1evrl55ea2XvhAoMwGEV0u34oR30lv/tkNOvMTuGHqObbbR+nbNzffbYxNesCdGXWXNmi19p/feR5+2tb7T1Tbf5Lv230+m2NGnD7Sje+5pR/x6N4RXO/cKXY3w+sTjJVDk1bTy8lktXSFQbgIIr5ZvxQnv4GGj7e13P7ARV/aNyU35/Evb4YBT7b7hF9jGG6zTIM0BV9xiH3861W4Z0g/hbZBW0x2AQPlkh/D6cPXq6pWX13rpC4EyE0B4tXQrTnj7XHiDdeqwkvU/5dBact/f/gi7/rLTbLutN18uzehZ3u49+9ieO21tZ/Q+COHVzr1CVyO8PvF4CRR5Na28fFZLVwiUmwDCq+VbccJ7bN8rbcNua9cKa4Rvy9172/l9jrA9d/rZcmmed+WtNmbCS/bYHZfHH3aLXrPmLEiWQLXZwurqZGMZFUSgeT2PogQ1qGfw5CmLbMTIaps2vZnaivrFCETCe8zhVbZa1+V/UDQUGnmFEks23isvYz9MFkAjRnnsh41YBiUOBNq1aeHQtXJaVpzwRnd4ow+qnXPyIbUpJ7nDe/3IB+26kQ/aPcPOs802Wq+2dtbs+YnOls+mLrKxT5jNnJnuD/pEk5d4UPv2ZrvvUmWrdkn3ICdPqbYRtyG86VI1qxWoVdLtTF7p8qzp5pUX+6FPXl77oc9q6RpKoF3but8SFVpf6eMrTnijZ3jfee9Du2lQnzj7hp7hXbSo2gYPG2WjH3nabrv6LNvke+vWOWeSfksDb7n6XGq8Re7D1asreXmR9elLXj5cvbp65eW1XvqGEeCRhjBeS46uOOH937c09LfNNl7frh5xn42Z8GLttzSMHD3OJkx8Lf4Wh+h17sCb7YGxE23YwDNs/XVWr+W36iqdrEXz5om/lgzh1U7UZVV7bfDkRV4+BJpWV64v8mpaBMq9WoRXy7fihDf6Ht5rb33Aht3+cEyubZsV7KZBZ9gWm24Q/3nQ9ffEd3NfGTss/nP3nn3to0+nLkV5zJ0DbZ01V0V4tfNPruYHsoww0wbklSlueTLykhFm2sArr0wPgsmWSQDh1U6OihPeGlxzv5ln07+Yaat17Vzv9+8mxcojDUlJ+Yzz2uC5w0tePgSaVleuL/JqWgTKvVqEV8u3YoVXw/a/aoQ3LZKN68MP5MZxy6uKvPIi37h5yatx3PKq8sorr+Nh3roEEF7tjEB4NX480iDyU8u9Nnju8KrJ1F9PXj5cvbqSlxdZn75eefmslq6hBBDeUGJ1xyO8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfnl6+2wAACAASURBVGq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY4Ir8YP4RX5qeVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo0jwqvxQ3hFfmq51waP8KrJILw+BLPtyvWVLW91Nq+81HVRnw4BhFfjiPBq/BBekZ9a7rXBI7xqMgivD8Fsu3J9Zctbnc0rL3Vd1KdDAOHVOCK8Gj+EV+Snlntt8AivmgzC60Mw265cX9nyVmfzyktdF/XpEEB4NY6pCu/IUeNs3bVWs59vtZm1aN5cW1kTqf5k2pxEK0WgEmEKHuS1wZNXcBSJCsgrEabCDCKvwkSRaCFeeSWanEHuBBBeDXGqwnvBVbfZ6IefslVX6WSHH7Sb/ar7z61D+xW1FRa8GuHNNyCvDR7h9cmVvHy4enUlLy+yPn298vJZLV1DCSC8ocTqjk9VeKPWb779H7vnoSftwXF/iWc6aJ8d7Df77mgbdltLW2lBqxHefIPx2uARXp9cycuHq1dX8vIi69PXKy+f1dI1lADCG0rMWXhr2k//cpY9NO4vdsefxttnU7+wLX+4kR26/6623Tabl+pxB4RXOwHVaq8NHuFVk6m/nrx8uHp1JS8vsj59vfLyWS1dQwkgvKHEMhLeGTO/tofHP2e3jhobC2/bNivY7DlzbeWO7az3YfvawfvtrK28INUIb75BeG3wCK9PruTlw9WrK3l5kfXp65WXz2rpGkoA4Q0l5iy8b73zvo166Cm7f8yz8Uw7bruF9eqxs231o03snfc+tDvuG28vvvYPe/LeIdrKC1KN8OYbhNcGj/D65EpePly9upKXF1mfvl55+ayWrqEEEN5QYo7CW/OhtehubnQH98C9t7c1Vuuy1ApnzPraOrQrx4fZEF7tBFSrvTZ4hFdNpv568vLh6tWVvLzI+vT1ystntXQNJYDwhhJzFN4bbn/I1lxtFdtlu5/YCq1baSsTqxctqrYp076wLit3SPzM8IKFC61ZVTNr1qwq8ewIb2JULgO9NniE1yUuIy8frl5dycuLrE9fr7x8VkvXUAIIbygxR+HVlpJe9TMvvGF9LrwhfmY4ep13+uHxt0Us7zVn7jz79XHn27GH7G177bJ17dAJE1+zkwdcs1Tpa+OHW+tWLfmPJ9KLrVGdvDZ4hLdRcTRYRF4NIirUAPIqVBwNLsYrrwYnZkAmBBBeDXPqX0umLUevjsT1lz1OthOP6hE/VvH086/bKQOG2uN3D7I1V1+l3gmuHDbKbr1nbPxvA/sfV0d4n5j4qp196XC7b/gFdWrXXqOrVVVVIbx6ZFIHrw0e4ZViWWYxeflw9epKXl5kffp65eWzWrqGEkB4Q4nVHV864Y3u7h5/9hCbNH64tWrVMj7aPQ7pF8vvwfvtUi+tL2d8ZXPnzbNex19kpx970FLCe8HgkTbxwaH11vJIg3YCqtVeGzzCqyZTfz15+XD16kpeXmR9+nrl5bNauoYSQHhDiZVceEc/8rSNHDXWxtw5sPZIT+p/ta271up2Ru+Dlkure8++dtJR+y0lvNEd4n27b2utW7eyn2y+oXXffsva54IRXu0EVKu9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBrH0t3hHXHXYzbuqZfrPIIQPc+7Uts2dn6fI4KF981/vm+PP/1y/K0Sn3w2Lf6vk3v12Mn6n3Jo3GvWnAWJEpg8ZZGNGFlt06Y3SzSeQckIRBv8MYdX2Wpdk3/QMEln8kpCKXwMeYUzy7OCvPKkHz63V17hK6HCg0C7Ni082lZMz9IJb9p3eJc8E6LvFx5wxS32xoSb47u8s2bPT3SyTJ5SbSNuQ3gTwQoYVLvB1/94dkCnukPJq9HolltIXj5cvbqSlxdZn75eefmslq6hBNq1/fYxTV6NI1A64a15hvf1P4+wli2//W0oelThsAN3XeYzvDXo6nukYUmsE19603r3G2yvPn5T/NVrPNLQuBMvrSqvt/B4pCGthOr2IS8frl5dycuLrE9fr7x8VkvXUAI80hBKrO740gnv7Dnf2Ja7H2f9Tuhpver5loZZX822I08baEf33MN233GrmEb0/bvVi6ptr8POtt6H7WN77bx1rSzf9cAE27DbWrbJ99a1GbO+sr4XDrOWLZrbLUP6xbUIr3YCqtVeGzzCqyZTfz15+XD16kpeXmR9+nrl5bNauoYSQHhDiZVceKPDe/K5SRZ9UK3mde6ph1rPX+0U/3HGzK9tm31OsMX/7vTzr4+f01389ejtl9l6a69uV9042m6+e0ztP/1gk242aEDv2q84Q3i1E1Ct9trgEV41GYTXh2C2Xbm+suWtzuaVl7ou6tMhgPBqHEt3h7cGx8KFi2zy1OnWtXPH2ru1jUU195t5NnXal9ZuxbbWscNKddogvI2lmk6d1waP8KaTz5JdyMuHq1dX8vIi69PXKy+f1dI1lADCG0qs7vjSCq+GJXk1wpuclcdIrw0e4fVIy/ivhX2wunXl+nJD69LYKy+XxdI0mADCG4ysTgHCq/HjGV6Rn1rutcEjvGoy9deTlw9Xr67k5UXWp69XXj6rpWsoAYQ3lBh3eDViS1RzhzdVnMHNvDZ4hDc4ikQF5JUIU2EGkVdhoki0EK+8Ek3OIHcCCK+GmDu8Gj/u8Ir81HKvDR7hVZPhDq8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejSPCq/FDeEV+arnXBo/wqskgvD4Es+3K9ZUtb3U2r7zUdVGfDgGEV+OI8Gr8EF6Rn1rutcEjvGoyCK8PwWy7cn1ly1udzSsvdV3Up0MA4dU4IrwaP4RX5KeWe23wCK+aDMLrQzDbrlxf2fJWZ/PKS10X9ekQQHg1jgivxg/hFfmp5V4bPMKrJoPw+hDMtivXV7a81dm88lLXRX06BBBejWPFCu+8efPtixlfWdcuHa2qqqpBitXV1bZw0SJr0bx5nbGfTJvTYG00AIFKhCl4kNcGT17BUSQqIK9EmAoziLwKE0WihXjllWhyBrkTQHg1xBUnvJG43nD7w3bdrQ/E5Fbu2M6uvfRU23yTbssl+cj4523I8HvtyXuHILzaOZdqtdcGj/CmGlNtM/Ly4erVlby8yPr09crLZ7V0DSWA8IYSqzu+4oR30lv/tkNOvMTuGHqObbbR+nbNzffbYxNesCdGXWXNmi19p/fDjz+z3/a50j76dKqtukonhFc731Kv9trgEd7Uo4obkpcPV6+u5OVF1qevV14+q6VrKAGEN5RYhQvv4GGj7e13P7ARV/aNSUz5/Evb4YBT7b7hF9jGG6yzFM0FCxfa59Nn2JN/mWQj7noU4dXOt9SrvTZ4hDf1qBBeH6SuXbm+XPGm3twrr9QXSsNGEUB4G4Wttqji7vD2ufAG69RhJet/yqG1EL6//RF2/WWn2XZbb75MmmOffMkG3XAPwqudb6lXe23wCG/qUSG8Pkhdu3J9ueJNvblXXqkvlIaNIoDwNgpb5QrvsX2vtA27rW1n9D6oFsKWu/e28/scYXvu9LNg4Z01Z0GiBCZPWWQjRlbbtOnNEo1nUDIC0QZ/zOFVtlrXhj94mKzjt6PIK4RW8rHklZxVEUaSVxFSSL4Gr7ySr4CRngTatWnh2b70vSvyDm/0QbVzTj6kNlzlDu+s2fMTnSSTp1TbiNsQ3kSwAgbVbvCrBBQlGEpeCSA1Ygh5NQJajiXklSP8RkztlVcjlkKJA4F2bVs6dK2clhUnvNEzvO+896HdNKhPnHJDz/DWnAo80lDMi8LrLTweafDJm7x8uHp1JS8vsj59vfLyWS1dQwnwSEMosbrjK054//ctDf1ts43Xt6tH3GdjJrxY+y0NI0ePswkTX4u/xSF6RV9jtmDBQhv31Mvx15I9ftcgq2pWVft9vHwPr3YCqtVeGzzCqyZTfz15+XD16kpeXmR9+nrl5bNauoYSQHhDiVW48EYCe+2tD9iw2x+OSbRts4LdNOgM22LTDeI/D7r+Hhv9yNP2ythh8Z/fff9j2/fI/nWo7b3rNnb5OcfGf4fwaiegWu21wSO8ajIIrw/BbLtyfWXLW53NKy91XdSnQwDh1ThW3B3eGlxzv5ln07+Yaat17Vzv9+8mxYrwJiXlM85rg0d4ycuHQNPqyvVFXk2LQLlXi/Bq+Vas8GrY/leN8KZFsnF9+IHcOG55VZFXXuQbNy95NY5bXlVeeeV1PMxblwDCq50RCK/Gj0caRH5qudcGzx1eNZn668nLh6tXV/LyIuvT1ysvn9XSNZQAwhtKrO54hFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBpHhFfjh/CK/NRyrw0e4VWTQXh9CGbblesrW97qbF55qeuiPh0CCK/GEeHV+CG8Ij+13GuDR3jVZBBeH4LZduX6ypa3OptXXuq6qE+HAMKrcUR4NX4Ir8hPLffa4BFeNRmE14dgtl25vrLlrc7mlZe6LurTIYDwahwRXo0fwivyU8u9NniEV00G4fUhmG1Xrq9seauzeeWlrov6dAggvBrH0grvokXVNmXaF9Zl5Q7WonnzBimFjq9p+Mm0OQ32jgYgUIkwBQ/y2uDJKziKRAXklQhTYQaRV2GiSLQQr7wSTc4gdwIIr4a4lML7zAtvWJ8Lb7DZc+bGdM47/XA7aJ8dlklqeeMnTHzNTh5wzVK1r40fbq1bteQOr3b+ydVeGzzCK0dTbwPy8uHq1ZW8vMj69PXKy2e1dA0lgPCGEqs7vnTCO2fuPPtlj5PtxKN62MH77WxPP/+6nTJgqD1+9yBbc/VVlqLV0PgnJr5qZ1863O4bfkGd2rXX6GpVVVUIr3b+ydVeGzzCK0eD8PogzLQr11emuOXJvPKSF0aDVAggvBrG0glvdLf2+LOH2KTxw61Vq5YxnT0O6RfL78H77bIUrYbGR8J7weCRNvHBofWS5pEG7QRUq702eIRXTab+evLy4erVlby8yPr09crLZ7V0DSWA8IYSK/kd3tGPPG0jR421MXcOrD3Sk/pfbeuutbqd0fugpWg1ND4S3ugO8b7dt7XWrVvZTzbf0Lpvv2Xtc8EIr3YCqtVeGzzCqyaD8PoQzLYr11e2vNXZvPJS10V9OgQQXo1j6e7wjrjrMRv31Mt1HkGInuddqW0bO7/PEUvRamj8m/983x5/+mXr0G5F++SzaTb64aesV4+drP8ph8a9Zs1ZkCiByVMW2YiR1TZterNE4xmUjEC0wR9zeJWt1rUqWUHCUeSVEFTgMPIKBJbzcPLKOYDA6b3ysmqzhdXVgatheBICzZsl/9nVrk2LJC0ZswwCpRPehu7YLskhdPz9Y561AVfcYm9MuDm+yztr9vxEJ9fkKdU24jaENxGsgEG1G/zSj2cHdFl6KHlJ+JZZTF4+XL26kpcXWZ++Xnl9NnWRjX3CbObM5HLmc4Tl6tq+vdnuu1TZql2SHVe7tt8+psmrcQRKJ7w1z+S+/ucR1rLlt78Nde/Z1w47cNflPsObdPzEl9603v0G26uP32QrtG7Fh9Yad96lVuX1Fh6PNKQWUZ1G5OXD1asreXmR9elLXj5cvbqG5sUjDVoSpRPe2XO+sS13P876ndDTetXzLQ2zvpptR5420I7uuYftvuNW1tD4ux6YYBt2W8s2+d66NmPWV9b3wmHWskVzu2VIv5g8z/BqJ6BaHbphJJ0P4U1KKmwceYXxyns0eeWdQNj85BXGK+/RoXkhvFpipRPeCMeTz02y6INqNa9zTz3Uev5qp/iPM2Z+bdvsc4It/nfLG3/VjaPt5rvH1Pb6wSbdbNCA3rVfcYbwaiegWh26YSSdD+FNSipsHHmF8cp7NHnlnUDY/OQVxivv0aF5IbxaYqUU3gjJwoWLbPLU6da1c8faRxuWh2p54+d+M8+mTvvS2q3Y1jp2WKlOG4RXOwHV6tANI+l8CG9SUmHjyCuMV96jySvvBMLmJ68wXnmPDs0L4dUSK63waliSVyO8yVl5jAzdMJKuAeFNSipsHHmF8cp7NHnlnUDY/OQVxivv0aF5IbxaYgivxo9neEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcC7jTq+gAAHjlJREFUYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwaP4RX5KeWh24YSedDeJOSChtHXmG88h5NXnknEDY/eYXxynt0aF4Ir5YYwqvxQ3hFfmp56IaRdD6ENympsHHkFcYr79HklXcCYfOTVxivvEeH5oXwaokhvBo/hFfkp5aHbhhJ50N4k5IKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjCq/FDeEV+annohpF0PoQ3KamwceQVxivv0eSVdwJh85NXGK+8R4fmhfBqiSG8Gj+EV+SnloduGEnnQ3iTkgobR15hvPIeTV55JxA2P3mF8cp7dGheCK+WGMKr8UN4RX5qeeiGkXQ+hDcpqbBx5BXGK+/R5JV3AmHzk1cYr7xHh+aF8GqJIbwJ+c36arYtWLjQOnVoV6fik2lzEnVAoBJhCh4UumEknYC8kpIKG0deYbzyHk1eeScQNj95hfHKe3RoXgivlhjC2wC/2XPmWr+Lb7Qnn5sUj/zBJt1s6MUnW5eVO8R/Rni1E1CtDt0wks6H8CYlFTaOvMJ45T2avPJOIGx+8grjlffo0LwQXi0xhLcBfiPueszufeRpu2Nof2uzQiv73VlDbL21V7eLzjwK4dXOvVSqQzeMpJMivElJhY0jrzBeeY8mr7wTCJufvMJ45T06NC+EV0sM4W2A3wG/Pc+6b7+l/fbgveKRjz/9sp1+/vX21lO3WlVVFXd4tfNPrg7dMJJOiPAmJRU2jrzCeOU9mrzyTiBsfvIK45X36NC8EF4tMYS3AX5b7t7bLu53dCy90esf//o/O/DY8+35R66zDu1WRHi180+uDt0wkk6I8CYlFTaOvMJ45T2avPJOIGx+8grjlffo0LwQXi0xhHc5/Kqrq23THY606y87zbbbevN45Hv/97Htc0R/e2LUYFt91c6J6X/48QJ7aNxCmzmzKnENAxsm0L692b57Nre1V2ve8OCAEeQVACtgKHkFwCrAUPIqQAgBSyCvAFgFGOqVVwEOrZBLQHgbiCW6w3vJWcfYrtv9JB655B3eQqbKoiAAAQhAAAIQgAAEagkgvA2cDNEzvLvt8FM7ptee8cgln+HlXIIABCAAAQhAAAIQKDYBhLeBfIb/8VG779Fn4m9paNumtfXud1Wdb2kodrysDgIQgAAEIAABCEAA4W3gHPh69lzrc+EN9uyLb8QjN91wPRt6ySnWtUvHJn/2/Os/H1n0H2r8+Affsy9nfGUvvPp3233HrTI/rkWLqm3cUy/btj/dNP4gIC8IQAACEKhLIHp38Sebb2SdO7UvNZosfhZVCstSnyiNODiENyG0GbO+tvnzF9T+hxMJywo97PrbHrJ/vvuBXXPRyfbm2/+x3/zuwtqvW8ty4RHXH+5yjN03/ALbeIN1spyauSAAAQjkQqDfJTfGj8ptsN6aieb//vZH2O3XnBPfoCjzS/lZ9MwLb8Q/y048qkctopvvHmNrrt7Fum//09q/qxSWZT5PGnNsCG9jqJWkBuEtSZAcBgQg0OQIRNJ165Cz7KdbbJRo7ZUiaYrw/vH+J+J3C+8Yek4t05MHXGMbfXcdO/7wfRHeRGdaeQchvE0g29fe/JcNuek+++e7H8a/qR56wK623x6/tOgOwfOvvGXTv5xl3db5jp1wZI/a7wu+/Nq7bO01VrUZs76y51/5u/X81U62/Tab2xXX3WOPTXjRVmjd0tq2WcE2/O5ade7wRr8ZR/+z3Kyv5thxh+5d+2G96A73FdfdbeOf+au1W6mNHbDX9nbsIXtZi+bN7a9vvGMXXnWbfTplekxzh21+aP1PPTR+POHd9z+2/pePsLNO6mV33Dfepnz+pd15bX974a9/t8uG/tHe++CT+L9r/ts/3qvYO7wRo7Muvcl2+eVPbNTDT8bsI7Y1/9lJxO3WUWPts6lf2Mod28VZ/u7wfeP/+OSR8c/bU8+/biu2XSHe6KN/P/fUQ+0XW/0gzuKTyZ/HnF987W3b/Pvd7MC9tq89R3oef1E8z8SX3rS3//1B/H3T0XnEq/EEoq8y/NNjz9of7/+zffTp57bRd9e204870NZba3X73dlD4ushen1/w3Xt7JMOtg27rRX/Ocoium6i6+uDjz6z3+y7ox1/xK/i/92xUl9zv5lng4eNis/rud/Mj8/f/icfEn+G4tPPptmgG0bZK6+/bS1btrCdf/FjO+fkQ6y+fW/aFzMaff0syT7q37xZM3vvg4/j62brn3zfzjqhlw2/61F78i+TYnk9+ej9a3Nd1vV31Y2j7ds7j6tYx/YrWY89fmG7bf/T5Z4jlSa8oT+LPv70czvkxIvjn4fRo4fR65ADdrELr7o9/nn3nVW72Abrrxnvc4uznDN3nl094j577IkXrFOHdvbrfXew/fbYrqKvvbLuOQhvwZP98OPPbPeD+8WCu98ev7D/++9ke/3v79oFfY6Mf6h+d701rXPH9vb0C6/bkJvutecfvs46tF8x/i+Qo+eOo7dxoh8Um220vj3y5+ftmRdetxOO6GHfXW8NG3b7w9ayZfM6wrvnTj+zvXfdxl589R82cvQ4e/zuQfGmfOZFw2LhPv24g2z6lzPtsqF32am/PcAO3m9ne+ud9+3f//kofhxhztxv7LxBt9r22/wwHlvz2/qqq3Sy/ff4pa2wQuv4K95263Wm7dt921icJ0+Zbn0vuqFihbeGUQ37lya9bbfeM9bG/nFg/EtLJEEtWjS3tb6ziv334yl20rnX1H439MhR42zQDfdY78P2sR9s3M1GP/JU/MvDxAeH2vwFC23fI86xH37/u/EvSe9/ODnmPP6eK22N1brEm370Oni/Xew7q3WOz5XVu65c8Cui2MuLfgGJfnmJpGfrH29iz//179a+3Yq2584/swfGTrQfbbqBtWrV0m65e4z958NP43M+FuDtj4h/2eh92L7xh2P7XjTMrjr/+NpfXIp91D6ri/5b99tGj7NrLz3VmjdvZk89N8l+9qNN4vN53yP7W9cunezonntY9BmA6MPF0S/S9e17n0+f0ajrp76jivpHv+DX/BJz3pW32kefTo1/OY3k9/Z7H7d2K7W1y885drnX3+w5c+1XR55rZ57Q0zbZYB1brevK8XnS0DlSSY80hP4s6rH7L2zITaPtpdfetgGnHRbHt8Zqq9iZFw+ztdfoatG/r7Rim/jn1OLCe/6VI+Nf+E877sD4JsIFg0fa7w7bN/45yKtcBBDegud57S0PxHf9nn3gmvhiXPy1cOEie+e9D2MRje6cDr3lfht143nxb7fRxhzdPYqkNHrVPCd70ZlHxfIcvRp6pGGPQ/rFG3n0v8xF30c8aMDvbI+dvv1QW3Sn46XX/mEP3HJx/Oep07601978t035/ItY0Nq3a2vXXXpqrfC+PGZYfBcyet14xyN255/G1x5TpT/DW99beDXso006ekX/4ck//vWBTZ3+ZSzDxxy8lx1+YHeLhPcvr7xpI67sG4+LzoMdDjjVxtw50D6dMs2OPv0Ku+3qs2vZR5v7vrv93Hr12Cne9IcNPMN+sdVmBb8Kms7yDjnxElvzO6vEwrPkK7qT9Le337P/+/BTe/Of78dy8/enR9YKbyRsW2y6Qfzn6N2bLp06WN/jf9N0Dj7llUZ7X/RL+jUXn2zfW3/N2v0venfomD6D4nN8nTVXrTPrkvtezT825vpZsnfUK+r/o802qH335Q/D77N/v/9RvNdFr6een2S/v+KW+BfOF1/7R4PX35KPNDR0jlSS8L711K21mSf9WRT6SMMm31vXfrLbsdb/lENti02/G2d4/5hn7bPPv4hvBPEqFwGEt+B5Rj/4otfA/sfVWWn07RHRV6RFsrvjz7ew1bt2ju9y3H39gPgRgSU35v9+MiW+q/rI7ZfZ+muvHvdqSHhPP/+6+C2eQ/bfxfY67Ow6P2Ae/fMLdsFVt9krY4fZ2Cdfir/J4kebfc823mBti779IXoLKZKp+mTu3IE327x58+2KAb3jdSC8S39gsIZ9dKci+uUieqxhx223sHXWWs3GTHjRDt1/VzvyN7svJbwRz2//O+yjLDpHBlxxS61E1ZxAO2y7RXxnLBLexSWr4JdCk1hexP6sE3vZ/nt++0tlzSt6lOHI0y6P7/5t+cON7Jt58+PHUZYlvJdcfYctWLjIzjv98CZx3B6LjB6R6n/ZcIve8Ygev+r5qx3jO+DRfhNdE9Hes+RryX0v+vfGXj+Lf8ipZp4l+9905yP2xj/eqxXeGsmNco3EqaHrb3HhTXKOVKrwJv1ZFCq80SNg0c+26K7vCq3/9/hQ9C1MV51/gsdpTc8cCSC8OcJPMvWVw0bZsy+8YQ/fdmmd4RMmvmbRw/jPP3Jd7Vd5RQKzLOGN3t7+4c5Hx3cCo7fekghv955947uBB/fY2bbZ54R4U48eVYhe0d2XMU++GEvwPoefY7vtuFXthwJuuWeMvTzp7WUK7233Pm5/fuavsWwhvFbvLwU7HniaHbj39nbQ3jvYL3ucbLcM6WdbbbFxzKt3v8G21Rab1Cu8H0/+3Hb9TR8b+YezbPacb+JfRF549Lr4WeslXwhvkiswbEyPo861rX60SSy9i78GXnd3/LbpzYPPjN+ejySp1/EXIbwJ8EbP6778+j/t4j/cYWef1MtW7tjeTjjnD/bM/Vcv9a05SwrptC9mNvr6iX4xaUioo5sM0SNmNXd4Fxfe6BsDGrr+br7qzPgxjeiV5BypVOFN+rPorgcmxDcEan62RFzjD611Wzt+Jr7mVfNIQ/Ro3zZ7n2D33nS+RXd7eZWbAMJb8HyjZ2mPPuMK+/1ph9neu24bv00dfVAt+iqb6O/vv/kiW22VleMPokV3hZYlvDUXfvQYRL8TetqMmV/Hd2ijZzcX/1qysX+8Iv6O4eitxOjt75pHJKK3aldacQU77/Qj7IsZs+y0866zXbfb0s7ofZBF/xZ9GOD0Yw+Mn2eL6jp1XGmZwhvdld7/mN/Hj0hEH/KI7hZHz6FW6teS1dwFjx4P6dq5o90/9lkbPGx0nO3qq3a2rfc6Pv6gRcQ7en4w+iEafeK45g7vQ4//xW68oo99M2+eXTfyQXvu5Tdt/D2Dbd78+bbzQWfEz66d+tv94zP9ldffsfkLFsQf8kF407/4r7v1AbvnoSft0rN/G/9i+dfX37GvZs+xd979MP5w4Q2Xn2YLFiyMc1reIw3c4bX4MwrRnbfoHavo3Yrol4m+v/uN/ezHm9iuv+lre+2ydfysZfR8e/TsbPT41pLCO/Or2Y2+fur7wOCS/ZcnvNEHfZd3/R112kDbcouN7Jhee9ns2XPjx7waOkcqSXgb87Mo+oD3cWdeFX/+IfrFMvpAYJRRtG9G358fnUfRXd3Fn+GNcohuCEXvOHZZuUP8mOCrf/tX/MgYr3IRQHibQJ7Rh8cGXX9P7UqjDyhFHzyL3ub587N/jf8+erv7yecm2T03/N4223j9eOOPvq+x5r9EjsZMeuvfdmzfwRZ9YCJ6izD67XaVzh2+Fd5/vm+/6X1B/PfRv0evSIwP+/8X/fsffmqnDBgaf6tC9Iru9EbPKUZv0T73ylt21iU3xp+Ojb/5odta8d9HP9xr+i7+PFb0IZPogwTRW5M1vZ5+/nX704gL40+1V9qrRnijjThiGL0Wf9Y6+jR39Knu6BV9sCl6Ozz6poYjfr1b/EhD9MtCzSv6gOGgAb1jSajJPPqWjOiT/9EryifKbadf/AjhdTjRomcwL/7D7fbguL/U8h547nGxuJ3U/+r4Lm/0ip6bjj7lv7xHGhYuqo5/0a3UV/ROUfSLX815G33Y9YK+R8bvVkTP8fYfOCL+5pLoFT1OFX0VVX37nnL9LMl+yf5LCm/0+MWJ51xd+7hFtOcu6/qL3qU7f/Ct8TUfifv+e23X4DkSHWN0rGV+KT+LFixcaCee84f42opefx13k02eMi3+WRk9ahc9Ix/d/Y2Et4ZldA6dP3hk7X8uFdVF31AUffCUV7kIILxNJM/ozmz09lzH9ivGn/KueUWfQI5+k42etU3yin6TjT5YFt0VjuqWfEX//tnU6da5U4d6v5Yl+lBU69Ytl/of0aKNJnrrcbWuna1li6XfPq9vbdHao68UqvT/Xa1GeN+YcHN85z26K7FkNtGdiehu1ZLfolDzobUbLjvNZn09J757Ud+r5j9Oif6XpiU//JjkvGFMGIHoGfUvZ34d/69Yi2cZfU1Vxw7t4m9i4NUwgWhfmTZ9pnVeuX29j+VEe2L07GXNB2KX1VG9fhpe6fJHLOv6i/b16B2zxa9LzpFvWTb2Z1FUG/Fu1bJlnZ9h0bkSfRPGsn4+RV+DF+2/yzrX1HOA+vwJILz5Z8AKKpyA8kXrS35LQ4Wj5PAhEESA6ycIF4Mh0KQJILxNOj4WXwYC0Vua0fcj13wFWcgxRW/TRW/J8dViIdQYC4FvCXD9cCZAoHIIILyVkzVHCgEIQAACEIAABCqSAMJbkbFz0BCAAAQgAAEIQKByCCC8lZM1RwoBCEAAAhCAAAQqkgDCW5Gxc9AQgAAEIAABCECgcgggvJWTNUcKAQhAAAIQgAAEKpIAwluRsXPQEIAABCAAAQhAoHIIILyVkzVHCgEIQAACEIAABCqSAMJbkbFz0BCAAAQgAAEIQKByCCC8lZM1RwoBCEAAAhCAAAQqkgDCW5Gxc9AQgAAEIAABCECgcgggvJWTNUcKAQhAAAIQgAAEKpIAwluRsXPQEIAABCAAAQhAoHIIILyVkzVHCgEIQAACEIAABCqSAMJbkbFz0BCAAAQgAAEIQKByCCC8lZM1RwoBCEAAAhCAAAQqkgDCW5Gxc9AQgAAEIAABCECgcgggvJWTNUcKAQjkRGDoLffby5P+aZecdYytvUbXeBXvvPdfu2jI7fbrfXawvXfdJqeVMS0EIACByiCA8FZGzhwlBCCQI4HPp8+wHkeda6uusrLddd25Nn/BQjvw2POsc6cOdsuQftayRfMcV8fUEIAABMpPAOEtf8YcIQQgUAACr735Lzv0pEvt4P12sRmzvrLnX3nLHrjlYuuycocCrI4lQAACECg3AYS33PlydBCAQIEI3Hbv43bFdXfHKxp143m26YbrFWh1LAUCEIBAeQkgvOXNliODAAQKRmDiS3+z3v2uilc15s6Bts6aqxZshSwHAhCAQDkJILzlzJWjggAECkbgo0+nWo+jBthuO/zUXv3bO9aieXO7Z9h51rZN64KtlOVAAAIQKB8BhLd8mXJEEIBAwQjMmTvPeh1/oTVv3jz+0Np/PvzU9j/m9/G3M1x+zrEFWy3LgQAEIFA+Aghv+TLliCAAgYIROO/KW+2+R5+xsX8caGuv8e1jDPc89GT8tWQX9DnSDthru4KtmOVAAAIQKBcBhLdceXI0EIAABCAAAQhAAAJLEEB4OSUgAAEIQAACEIAABEpNAOEtdbwcHAQgAAEIQAACEIAAwss5AAEIQAACEIAABCBQagIIb6nj5eAgAAEIQAACEIAABBBezgEIQAACEIAABCAAgVITQHhLHS8HBwEIQAACEIAABCCA8HIOQAACEIAABCAAAQiUmgDCW+p4OTgIQAACEIAABCAAAYSXcwACEIAABCAAAQhAoNQEEN5Sx8vBQQACEIAABCAAAQggvJwDEIAABCAAAQhAAAKlJoDwljpeDg4CEIAABCAAAQhAAOHlHIAABCAAAQhAAAIQKDUBhLfU8XJwEIAABCAAAQhAAAIIL+cABCAAAQhAAAIQgECpCSC8pY6Xg4MABCAAAQhAAAIQQHg5ByAAAQhAAAIQgAAESk0A4S11vBwcBCAAAQhAAAIQgADCyzkAAQhAAAIQgAAEIFBqAghvqePl4CAAAQhAAAIQgAAEEF7OAQhAAAIQgAAEIACBUhNAeEsdLwcHAQhAAAIQgAAEIIDwcg5AAAIQgAAEIAABCJSaAMJb6ng5OAhAAAIQgAAEIAABhJdzAAIQgAAEIAABCECg1AQQ3lLHy8FBAAIQgAAEIAABCCC8nAMQgAAEIAABCEAAAqUmgPCWOl4ODgIQgAAEIAABCEAA4eUcgAAEIAABCEAAAhAoNQGEt9TxcnAQgAAEIAABCEAAAggv5wAEIAABCEAAAhCAQKkJILyljpeDgwAEIAABCEAAAhBAeDkHIAABCEAAAhCAAARKTQDhLXW8HBwEIAABCEAAAhCAAMLLOQABCEAAAhCAAAQgUGoCCG+p4+XgIAABCEAAAhCAAAQQXs4BCEAAAhCAAAQgAIFSE0B4Sx0vBwcBCEAAAhCAAAQggPByDkAAAhCAAAQgAAEIlJrA/wN9NMsD8JxBlgAAAABJRU5ErkJggg=="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#| caption: A discrete probability distribution as a bar graph.\n",
    "#| label: fig:discrete-distribution\n",
    "px.bar(y=PMF, x=categories)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rEvq8irysNEQ",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simulation by Sampling\n",
    "\n",
    "> We can simulate our trash sorting cell by sampling from the prior."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "JOy1xtoe7ear",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```{index} simulation, sampling\n",
    "```\n",
    "\n",
    "Suppose we wish to simulate our trash sorting system such that the behavior of the simulation matches, \n",
    "in some statistical sense, the behavior of the actual system.\n",
    "In our case, this amounts to generating samples from the probability distribution\n",
    "on trash categories.\n",
    "In particular, we would like to generate a sequence of categories,\n",
    "$\\omega_1, \\omega_2, \\dots, \\omega_n$ such that \n",
    "$\\omega_i = \\mathrm{cardboard}$ approximately 25% of the time,\n",
    "$\\omega_i = \\mathrm{paper}$ approximately 20% of the time, etc.\n",
    "How can we write a computer program to do this?\n",
    "\n",
    "```{index} uniform probability distribution\n",
    "```\n",
    "\n",
    "While most programming libraries do not include functions to generate samples from an arbitrary distribution,\n",
    "almost all include a random number generator that will generate a random number from the unit interval.\n",
    "We denote by $U(a,b)$ the *uniform probability distribution* on the interval $[a,b]$.\n",
    "In numpy, the function `np.random.rand()` generates a sample $x \\sim U(0,1)$. We will often use the symbol $\\sim$ to indicate a variable is sampled from a given distribution.\n",
    "Below we explain how we can use this to generate a sample from an arbitrary probability distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dnm1_zFnL19p",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Cumulative Distribution Function\n",
    "```{index} pair: cumulative distribution function; CDF\n",
    "```\n",
    "```{index} random variable\n",
    "```\n",
    "\n",
    "We begin by introducing the Cumulative Distribution Function (CDF) for a random variable $X$.\n",
    "We will more carefully introduce the notion of a discrete random variable in Section 2.2, but for now it\n",
    "is sufficient to know that a discrete random variable takes values from a countable set,\n",
    "each of which is assigned a probability value.\n",
    "For a random variable $X$, the CDF for $X$ is denoted by $F_X$, and is defined as\n",
    "\\begin{equation}\n",
    "F_X(\\alpha) = P(X \\leq \\alpha)\n",
    "\\end{equation}\n",
    "It follows immediately that $0 \\leq F_X(\\alpha) \\leq 1$,\n",
    "since $F_X(\\alpha)$ is itself a probability.\n",
    "In the case of discrete random variables,\n",
    "say $X \\in \\{ x_0, \\dots x_{n-1}\\}$, we can compute the CDF\n",
    "$F_X(\\alpha)$ by summing the probabilities assigned\n",
    "to all $x_i \\leq \\alpha$\n",
    "\\begin{equation}\n",
    "F_X(\\alpha) = \\sum_{x_i \\leq \\alpha} P(x_i) = \\sum_{i=0}^{k-1} P(x_i)\n",
    "\\end{equation}\n",
    "in which the rightmost summation follows if we choose $k$ such that $x_{k-1} \\leq \\alpha < x_k$.\n",
    "The terminology *Cumulative Distribution Function* stems from the fact that $F_X(\\alpha)$\n",
    "is the accumulated probability assigned to all outcomes less than or equal to $\\alpha$,\n",
    "which is apparent in these summation expressions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "qq755xF4HFdi",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "But what does this have to do with generating samples from our distribution on categories?\n",
    "The idea is simple: we can generate $x\\sim U(0,1)$, and a CDF takes on values in the\n",
    "interval $[0,1]$.  For a discrete random variable $X \\in \\{ x_0, \\dots x_{n-1}\\}$\n",
    "the probability that our sample $x$ corresponds to category $k$\n",
    "is exactly equal to $F_X(x_k) - F_X(x_{k-1})$, and we define $F_X(x_k)=0$ for $k < 0$.\n",
    "\n",
    "To see this, we impose an ordering on our categories,\n",
    "\\begin{align*}\n",
    "c_0 &= \\mathrm{cardboard} \\\\\n",
    "c_1 &= \\mathrm{paper} \\\\\n",
    "c_2 &= \\mathrm{can} \\\\\n",
    "c_3 &= \\mathrm{scrap \\; metal} \\\\\n",
    "c_4 &= \\mathrm{bottle}\n",
    "\\end{align*}\n",
    "and we define the random variable $X \\in \\{ 0,1,2,3,4\\}$ to be the index of the chosen category.\n",
    "The CDF for $X$ is given by:\n",
    "\n",
    "| $k$ | $x_k$ | $F_X$($x_k$) |\n",
    "|----|----------------|:------:|\n",
    "| 0 | cardboard    |   0.20 |\n",
    "| 1 | paper    |   0.50 |\n",
    "| 2 | can    |   0.75 |\n",
    "| 3 | scrap metal    |   0.95 |\n",
    "| 4 | bottle    |   1.00 |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "HqOZN-mR5hyZ",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example\n",
    "Some numpy code to generate the CDF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "A3fhqoo_rwWN",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "outputId": "b310dfef-21f8-48ba-e5d8-088c94cdb978",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.2  0.5  0.75 0.95 1.  ]\n"
     ]
    }
   ],
   "source": [
    "CDF = np.cumsum(PMF)\n",
    "print(CDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "o3hsxeV8n_Pr",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "\n",
    "Now, suppose we generate a random sample $x \\sim U(0,1)$,\n",
    "and use this $x$ to choose category $k$ such that $F_X(x_{k-1}) < x \\leq F_X(x_k)$.\n",
    "For example, we choose category 4 if $0.95 < x \\leq 1.0$.\n",
    "In this case, what is the probability of choosing category 4?\n",
    "The answer follows from the fact that,\n",
    "for the uniform distribution on the unit interval, $P(X \\in [a,b]) = b-a$.\n",
    "Therefore, the probability that our sample lies in the interval $[0.95,1.0]$ is\n",
    "$0.05$, which happens to be exactly the prior probability assigned to category 4!\n",
    "Likewise, if our sample $x$ satisfies $0.2 < x \\leq 0.5$,\n",
    "we choose category 1, and the probability that our sample lies in the interval\n",
    "$[0.2,0.5]$ is $0.3$, which is, as expected, exactly the prior probability assigned\n",
    "to category 1.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "KlIIihLa1t5c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The code to accomplish sampling is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "v63pYNenOwbY",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 542
    },
    "outputId": "55ffca22-626f-4563-a60b-d525bfa56484",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def sample():\n",
    "    u = np.random.rand()\n",
    "    for category in range(5):\n",
    "        if u<float(CDF[category]):\n",
    "            return category"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "Yst9AzWddrJy",
   "metadata": {},
   "source": [
    "If we simulate our system by generating a large number of samples, we expect that the frequencies of each category will approximate their prior probabilities. To see this, the code below generates a histogram for 1000 generated samples. While the histogram values do not exactly match the prior probabilities, we can see that the approximation is quite good."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72y0OJpL8eJq",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "```{index} histogram\n",
    "```\n",
    "\n",
    "We can use the plotly library to plot a set of 1000 samples generated using this method, as a *histogram*. The result is shown in Figure [2.2](#fig:discrete-samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "lnHzswHfuJhU",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4Xu3de7xcdXno/yfZEEhCmpCEYLhEKK20FLE9NfwOnlMDXrgIIiBBAZEAASMoyM2kBEyAIKCFABFoBMGCoAgVNShiFYIKRXoUsbyoxeIBKyp3MJAEctm/1wwmx5jA7KxnZq25vPdfbvZ8L+v9XbUfxrVnD+rv7+8PXwQIECBAgAABAgS6VGCQ4O3Sk3VZBAgQIECAAAECdQHB60YgQIAAAQIECBDoagHB29XH6+IIECBAgAABAgQEr3uAAAECBAgQIECgqwUEb1cfr4sjQIAAAQIECBAQvO4BAgQIECBAgACBrhYQvF19vC6OAAECBAgQIEBA8LoHCBAgQIAAAQIEulpA8Hb18bo4AgQIECBAgAABweseIECAAAECBAgQ6GoBwdvVx+viCBAgQIAAAQIEBK97gAABAgQIECBAoKsFBG9XH6+LI0CAAAECBAgQELzuAQIECBAgQIAAga4WELxdfbwujgABAgQIECBAQPC6BwgQIECAAAECBLpaQPB29fG6OAIECBAgQIAAAcHrHiBAgAABAgQIEOhqAcHb1cfr4ggQIECAAAECBASve4AAAQIECBAgQKCrBQRvVx+viyNAgAABAgQIEBC87gECBAgQIECAAIGuFhC8XX28Lo4AAQIECBAgQEDwugcIECBAgAABAgS6WkDwdvXxujgCBAgQIECAAAHB6x4gQIAAAQIECBDoagHB29XH6+IIECBAgAABAgQEr3uAAAECBAgQIECgqwUEb1cfr4sjQIAAAQIECBAQvO4BAgQIECBAgACBrhYQvF19vC6OAAECBAgQIEBA8LoHCBAgQIAAAQIEulpA8Hb18bo4AgQIECBAgAABweseIECAAAECBAgQ6GoBwdvVx+viCBAgQIAAAQIEBK97gAABAgQIECBAoKsFBG9XH6+LI0CAAAECBAgQELzuAQIECBAgQIAAga4WELxdfbwujgABAgQIECBAQPC6BwgQIECAAAECBLpaQPB29fG6OAIECBAgQIAAAcHrHiBAgAABAgQIEOhqAcHb1cfr4ggQIECAAAECBASve4AAAQIECBAgQKCrBQRvVx+viyNAgAABAgQIEBC87gECBAgQIECAAIGuFhC8XX28Lo4AAQIECBAgQEDwugcIECBAgAABAgS6WkDwdvXxujgCBAgQIECAAAHB6x4gQIAAAQIECBDoagHB29XH6+IIECBAgAABAgQEr3uAAAECBAgQIECgqwUEb1cfr4sjQIAAAQIECBAQvO4BAgQIECBAgACBrhYQvF19vC6OAAECBAgQIEBA8LoHCBAgQIAAAQIEulpA8Hb18bo4AgQIECBAgAABweseIECAAAECBAgQ6GoBwdvVx+viCBAgQIAAAQIEBG/yHvj100uSMxhOgAABAgQIEHhtgS3GDEWUEBC8CbzaUMGbBDScAAECBAgQaCggeBsSveYLBG/OT/Am/QwnQIAAAQIEGgsI3sZGr/UKwZvzE7xJP8MJECBAgACBxgKCt7GR4M0ZveZojzS0ENfUBAgQIECAQF1A8OZuBO/w5vy8w5v0M5wAAQIECBBoLCB4Gxt5hzdn5B3eFvqZmgABAgQIEGgsIHgbGwnenJHgbaGfqQkQIECAAIHGAoK3sZHgzRkJ3hb6mZoAAQIECBBoLCB4GxsJ3pyR4G2hn6kJECBAgACBxgKCt7GR4M0ZCd4W+pmaAAECBAgQaCwgeBsbCd6IeOnlZfHk08/FsKEbx+hRI9YyWbmyP554+tkYO3pkbNDXt9bPF72wOJavWBGbjlxzrI8ly92ARhMgQIAAAQKNBQRvY6OeD97Tz/9c3Hzr91c7/I83viHmzTk+Ro3cpP7P7vzX++OUsy6PxUuW1r+fddLhcdC+u9X/c+2fTZ8zP26/67769zvtsF19bC2Ma1+CN3cDGk2AAAECBAg0FhC8jY16PnjnX7sg/vfOb4w3bLd1/Obxp+LQ4+bEYQfuHsd84N2xZOnL8db9j4+PHLl/HHrAO2Lh3T+JE86YF7d98dOx1fjN4srrvxE3LlgY186bGUM3HhIfnjE3tp0wPs7++JGCN3fvGU2AAAECBAgMUEDwDhDqVV7Wc394Ytmy5fG2ySfGR488oP4ubu3d3WP/fm7c9+0rYsiQDetM7/rA9Hr8HnrAO+PAo2fFHrtOjKMP3af+s9sW3hsnzb4sHrjj6hg0aJB3eHP3n9EECBAgQIDAAAQE7wCQXuMlPRO8L7+8LK760q1x5z33x2ZjRsYnZxwdmwwfGl9esDA+f8Ot8c0vnL+a6aMzL45tth4fJ087KCbuNS3mTD+qHr21rwcfeiQmHzM77l5waYwcMVzw5u4/owkQIECAAIEBCAjeASAJ3qg/ujDzvCvjZ//1aIwbu2mce9oxMX7c6PojC9+649646YozVzPVnufdZNjQmHXy4bHjbkfEZeeeGJN2eVP95w8/8ljsO2VmfOeGC2L85mNi8UvLcydgNAECHSswqGN33n0b7+++S3JFBNYQGLbRBkQSAj3zDu8qo/7+/jj61H+I1202uv7O7UDe4T1nxtTYfdKb1/kO77OLXk7wG0qAQCcLiKz2OT3/8tE+Z2EnrRHYdMSQ1kzcI7P2XPDWzvWTl3whfvHL38SV/3Dq6md4f/IvV8aGG77yb097HHxqfHDy7quf4d1zt51j6iF713/mGd7O/r+MT5zl35A76QTP+oT/BaWTzsteCRBonYBHGnK2XR+8L7y4JD77hQWx/15/F1ttMa7+DO7Ukz9dD9gPHfbuWLzkpZi414di+nEHxyHr+JSGK667JW665c76pzQMG7pRTJt+oU9pyN1zlY4WvJXyr/figne9yQwgQKBLBQRv7mC7PnhfXLw0Dj/h3PiPnz+6Wmq/Pf93fOKkw2Oj338qQ+0zdmu/qLbq6/SPHRYH7/f2+re18bVner93z/3173fcftuYd84JMW7sqPr3Poc3dwOWPVrwli2eW0/w5vyMJkCgewQEb+4suz54V/HUwvXpZ5+PsaNH1d+p/eOvFStWxm+ffCbGjRm1+tGGP3zN84tejNpHmq36gxOrfiZ4czdg2aMFb9niufUEb87PaAIEukdA8ObOsmeCN8f06qMFb6tkWzOv4G2Na6tmFbytkjUvAQKdJiB4cycmeHN+HmlI+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOUvDm/ARv0q/s4YK3bPHceoI352c0AQLdIyB4c2cpeHN+gjfpV/ZwwVu2eG49wZvzM5oAge4RELy5sxS8OT/Bm/Qre7jgLVs8t57gzfkZTYBA9wgI3txZCt6cn+BN+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOUvDm/ARv0q/s4YK3bPHceoI352c0AQLdIyB4c2cpeHN+gjfpV/ZwwVu2eG49wZvzM5oAge4RELy5sxS8OT/Bm/Qre7jgLVs8t57gzfkZTYBA9wgI3txZCt6cn+BN+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOUvDm/ARv0q/s4YK3bPHceoI352c0AQLdIyB4c2cpeHN+gjfpV/ZwwVu2eG49wZvzM5oAge4RELy5sxS8OT/Bm/Qre7jgLVs8t57gzfkZTYBA9wgI3txZCt6cn+BN+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOUvDm/ARv0q/s4YK3bPHceoI352c0AQLdIyB4c2cpeHN+gjfpV/ZwwVu2eG49wZvzM5oAge4RELy5sxS8OT/Bm/Qre7jgLVs8t57gzfkZTYBA9wgI3txZCt6cn+BN+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOUvDm/ARv0q/s4YK3bPHceoI352c0AQLdIyB4c2cpeHN+gjfpV/ZwwVu2eG49wZvzM5oAge4RELy5sxS8OT/Bm/Qre7jgLVs8t57gzfkZTYBA9wgI3txZCt6cn+BN+pU9XPCWLZ5bT/Dm/IwmQKB7BARv7iwFb85P8Cb9yh4ueMsWz60neHN+RhMg0D0Cgjd3loI35yd4k35lDxe8ZYvn1hO8OT+jCRDoHgHBmztLwZvzE7xJv7KHC96yxXPrCd6cn9EECHSPgODNnaXgzfkJ3qRf2cMFb9niufUEb87PaAIEukdA8ObOsmeCd8nSl+PZ534Xrxs3JgYPHrTeaoteWBzLV6yITUeOWGPsr59est5zGVCdgOCtzr7IyoK3iJoxBAh0o4DgzZ1qTwTvR2deHLffdV9davSoEbHfnn8XJ087aLXcvoefFg8/+us1JI+bsl8cO2W/WLxkaUyfM3/1+J122C7mzTk+xo4eWX+94M3dgGWPFrxli+fWE7w5P6MJEOgeAcGbO8ueCN7PXHVz7L7rxJiw5bi450cPxnGnXRRfuvwT8ca//NO6Xi14937HLrHnbjuv1hw5YniMGrlJXHn9N+LGBQvj2nkzY+jGQ+LDM+bGthPGx9kfP1Lw5u69SkYL3krYCy8qeAvTGUiAQJcJCN7cgfZE8P4x0dsmnxjvf8/b4pgPvHt18E55355xwLveupbmgUfPij12nRhHH7pP/We3Lbw3Tpp9WTxwx9UxaNAg7/Dm7r/SRwve0slTCwreFJ/BBAh0kYDgzR1mzwXvo796PN71gelx2bknxqRd3rQ6eIcPHxrbvX6L2GLzMbHPO3eJCVtuXv/ZxL2mxZzpR9Wjt/b14EOPxORjZsfdCy6N2rvAHmnI3YBljxa8ZYvn1hO8OT+jCRDoHgHBmzvLngreFxcvjQ98ZE5sMnxYfP6iGdHXN7iud+nVN8fgvsHR3x9x+w9+HLUo/ucrz4yttxgXO+52xBpx/PAjj8W+U2bGd264IMZvPiaef3FZ7gSMLlXg5NNKXc5iSYELPpmcoMXD1//XX1u8oR6evr+Hr92l94bAyOEb9saFtugqeyZ4a5/ScMIZl8Rvn3gmrrnktPrzuev6WrZseexxyKlx2Ht3jyPev1f9Hd5zZkyN3Se9eZ3v8L6wRPC26N5sybQfm9GSaU3aIoGLzmvRxE2aVmQ1CbIJ0/iXjyYgmqKtBTYZKngzB9QTwfu7FxbH8adfEkuWvBTzP3Xyq8buKsj3fejMmPSWv45jD39P1J7hrf0y29RD9q7/2DO8mdut+rEeaaj+DNZnBx5pWB8tryVAoJsFPNKQO92uD97FS16K9087s/4ZunPP/EhsMnxoXWzw4MExftzo+OVjj9c/cqwWtWM2HRm33XFvTD9nfv1d4L/d6Q1xxXW3xE233Fn/lIZhQzeKadMv9CkNuXuu0tGCt1J+i3e5gH9B6fIDdnmVCgjeHH/XB+/jTz4btU9l+OOv2ufxfv+r8+rBO+Vj50Xtdau+ph93cHxw8h71b2vP/Z5y1uXxvXvur3+/4/bbxrxzTohxY0fVv/dLa7kbsOzRgrdscev1koDg7aXTdq1lCwjenHjXB+9AePr7++OZ5xbV/8hE7RfRNujrW2vY84tejNrzvav+4MSqFwjegQi3z2sEb/uchZ10n4Dg7b4zdUXtIyB4c2cheHN+3uFN+pU9XPCWLW69XhIQvL102q61bAHBmxMXvDk/wZv0K3u44C1b3Hq9JCB4e+m0XWvZAoI3Jy54c36CN+lX9nDBW7a49XpJQPD20mm71rIFBG9OXPDm/ARv0q/s4YK3bHHr9ZKA4O2l03atZQsI3py44M35Cd6kX9nDBW/Z4tbrJQHB20un7VrLFhC8OXHBm/MTvEm/socL3rLFrddLAoK3l07btZYtIHhz4oI35yd4k35lDxe8ZYtbr5cEBG8vnbZrLVtA8ObEBW/OT/Am/coeLnjLFrdeLwkI3l46bddatoDgzYkL3pyf4E36lT1c8JYtbr1eEhC8vXTarrVsAcGbExe8OT/Bm/Qre7jgLVvcer0kIHh76bRda9kCgjcnLnhzfoI36Vf2cMFbtrj1eklA8PbSabvWsgUEb05c8Ob8BG/Sr+zhgrdscev1koDg7aXTdq1lCwjenLjgzfkJ3qRf2cMFb9ni1uslAcHbS6ftWssWELw5ccGb8xO8Sb+yhwvessWt10sCgreXTtu1li0geHPigjfnJ3iTfmUPF7xli1uvlwQEby+dtmstW0Dw5sQFb85P8Cb9yh4ueMsWt14vCQjeXjpt11q2gODNiQvenJ/gTfqVPVzwli1uvV4SELy9dNqutWwBwZsTF7w5P8Gb9Ct7uOAtW9x6vSQgeHvptF1r2QKCNycueHN+gjfpV/ZwwVu2uPV6SUDw9tJpu9ayBQRvTlzw5vwEb9Kv7OGCt2xx6/WSgODtpdN2rWULCN6cuODN+QnepF/ZwwVv2eLW6yUBwdtLp+1ayxYQvDlxwZvzE7xJv7KHC96yxa3XSwKCt5dO27WWLSB4c+KCN+cneJN+ZQ8XvGWLW6+XBARvL522ay1bQPDmxAVvzk/wJv3KHi54yxa3Xi8JCN5eOm3XWraA4M2JC96cn+BN+pU9XPCWLW69XhIQvL102q61bAHBmxMXvDk/wZv0K3u44C1b3Hq9JCB4e+m0XWvZAoI3Jy54c36CN+lX9nDBW7a49XpJQPD20mm71rIFBG9OXPDm/ARv0q/s4YK3bHHr9ZKA4O2l03atZQsI3py44M35Cd6kX9nDBW/Z4tbrJQHB20un7VrLFhC8OXHBm/MTvEm/socL3rLFrddLAoK3l07btZYtIHhz4oI35yd4k35lDxe8ZYtbr5cEBG8vnbZrLVtA8ObEBW/OT/Am/coeLnjLFrdeLwkI3l46bddatoDgzYkL3pyf4E36lT1c8JYtbr1eEhC8vXTarrVsAcGbExe8OT/Bm/Qre7jgLVvcer0kIHh76bRda9kCgjcnLnhzfoI36Vf2cMFbtrj1eklA8PbSabvWsgUEb05c8Ob8BG/Sr+zhgrdscev1koDg7aXTdq1lCwjenLjgzfkJ3qRf2cMFb9ni1uslAcHbS6ftWssWELw5ccGb8xO8Sb+yhwvessWt10sCgreXTtu1li0geHPigjfnJ3iTfmUPF7xli1uvlwQEby+dtmstW0Dw5sQFb85P8Cb9yh4ueMsWt14vCQjeXjpt11q2gODNiQvenJ/gTfqVPVzwli1uvV4SELy9dNqutWwBwZsTF7w5P8Gb9Ct7uOAtW9x6vSQgeHvptF1r2QKCNycueHN+gjfpV/ZwwVu2uPV6SUDw9tJpu9ayBQRvTlzw5vwEb9Kv7OGCt2xx6/WSgODtpdN2rWULCN6cuODN+QnepF/ZwwVv2eLW6yUBwdtLp+1ayxYQvDlxwZvzE7xJv7KHC96yxa3XSwKCt5dO27WWLSB4c+I9E7xLlr4czz73u3jduDExePCgtdRWruyPJ55+NsaOHhkb9PWt9fNFLyyO5StWxKYjR6zxs18/vSR3AkaXKiB4S+W2WI8JCN4eO3CXW6qA4M1x90TwfnTmxXH7XffVpUaPGhH77fl3cfK0g1bL3fmv98cpZ10ei5csrf+zWScdHgftu1v9P9f+2fQ581eP32mH7WLenOPrYVz7Ery5G7Ds0YK3bHHrESDQrgL+BaVdT2bd+xK8ufPqieD9zFU3x+67TowJW46Le370YBx32kXxpcs/EW/8yz+N2ju/b93/+PjIkfvHoQe8Ixbe/ZM44Yx5cdsXPx1bjd8srrz+G3HjgoVx7byZMXTjIfHhGXNj2wnj4+yPHyl4c/deJaMFbyXsFiVAoA0FBG8bHsprbEnw5s6rJ4L3j4neNvnEeP973hbHfODdUXt399i/nxv3ffuKGDJkw/pL3/WB6fX4PfSAd8aBR8+KPXadGEcfuk/9Z7ctvDdOmn1ZPHDH1TFo0CDv8Obuv9JHC97SyS1IgECbCgjeNj2YV9mW4M2dV88F76O/erwetJede2JM2uVN8eUFC+PzN9wa3/zC+asla49AbLP1+PpjDxP3mhZzph9Vj97a14MPPRKTj5kddy+4NEaOGC54c/df6aMFb+nkFiRAoE0FBG+bHozgbcnB9FTwvrh4aXzgI3Nik+HD4vMXzYi+vsH1Rxa+dce9cdMVZ64Grj3Pu8mwoTHr5MNjx92OWB3HtRc8/Mhjse+UmfGdGy6I8ZuPiad+91JLDqbXJ1371wqbIzL9jMHNmcgsBAgQ6HCB889e2eFX0Jnb7y+47bF/slHBkYbVBHomeGvP6p5wxiXx2yeeiWsuOS1GjdykfgcM5B3ec2ZMjd0nvXmd7/C+9PIKd1ILBIr+F0KjrXzkVP8F38jIzwkQ6A2Bz3zaGwBVnHTRN3Q2GrL2J0hVsf9OXbMngvd3LyyO40+/JJYseSnmf+rk1bFbO7RVz/D+5F+ujA033KB+jnscfGp8cPLuq5/h3XO3nWPqIXvXf+YZ3k691V/Zt0caOvv87J4AgeYJeKSheZZlzOQZ3pxy1wfv4iUvxfunnVn/DN25Z34kNhk+tC42ePDgGD9udNR+PnGvD8X04w6OQ9bxKQ1XXHdL3HTLnfVPaRg2dKOYNv1Cn9KQu+cqHS14K+W3OAECbSQgeNvoMAawFcE7AKTXeEnXB+/jTz4btU9l+OOv2ufxfv+r8+r/uPYZvbVfVFv1dfrHDouD93t7/dvac7+1Z3q/d8/99e933H7bmHfOCTFu7Kj69z6HN3cDlj1a8JYtbj0CBNpVQPC268mse1+CN3deXR+8A+VZsWJl/PbJZ2LcmFGrH234w7HPL3oxli1bvvoPTqz6meAdqHB7vE7wtsc52AUBAtULCN7qz2B9diB410dr7dcK3pxfTD1hWXIGwwkQIECAQPkCgrd888yKgjej10Of0pBjevXRgrdVsuYlQIAAgVYKCN5W6jZ/bsGbM/UOb87PO7xJP8MJECBAoBoBwVuNe9FVBW9RuVfGCd6cn+BN+hlOgAABAtUICN5q3IuuKniLygnenNzvR3ukoSmMJiFAgACBkgUEb8ngyeUEbw7QO7w5P+/wJv0MJ0CAAIFqBARvNe5FVxW8ReW8w5uT8w5vU/xMQoAAAQLVCAjeatyLrip4i8oJ3pyc4G2Kn0kIECBAoBoBwVuNe9FVBW9ROcGbkxO8TfEzCQECBAhUIyB4q3EvuqrgLSoneHNygrcpfiYhQIAAgWoEBG817kVXFbxF5QRvTk7wNsXPJAQIECBQjYDgrca96KqCt6ic4M3JCd6m+JmEAAECBKoRELzVuBddVfAWlRO8OTnB2xQ/kxAgQIBANQKCtxr3oqsK3qJygjcnJ3ib4mcSAgQIEKhGQPBW4150VcFbVE7w5uQEb1P8TEKAAAEC1QgI3mrci64qeIvKCd6cnOBtip9JCBAgQKAaAcFbjXvRVQVvUTnBm5MTvE3xMwkBAgQIVCMgeKtxL7qq4C0qJ3hzcoK3KX4mIUCAAIFqBARvNe5FVxW8ReUEb05O8DbFzyQECBAgUI2A4K3GveiqgreonODNyQnepviZhAABAgSqERC81bgXXVXwFpUTvDk5wdsUP5MQIECAQDUCgrca96KrCt6icoI3Jyd4m+JnEgIECBCoRkDwVuNedFXBW1RO8ObkBG9T/ExCgAABAtUICN5q3IuuKniLygnenJzgbYqfSQgQIECgGgHBW4170VUFb1G5koL33vt+FiP/ZHhsv93Wa+z0yaefi3t+9GDs9fb/Lzbo68tdRYWjp56wrMLVLU2AAAECBIoJCN5iblWNErw5+UH9/f39uSlee/RHZ14cO2y/TXz4g+9Z44W//u1T8c73nxK3XHNubDthfCu30NK5BW9LeU1OgAABAi0SELwtgm3RtII3B1tZ8D740CMx+ZjZcet158eELTfPXUWFowVvhfiWJkCAAIHCAoK3MF0lAwVvjr1lwTvjk5+N555fFD/66c9j9KgRse2E163e6csvL48f3vcf8Zd//vq46Yozc1dQ8WjBW/EBWJ4AAQIECgkI3kJslQ0SvDn6lgXvGZ+6Kp5f9ELc9+8/jxGbDIs/23bL1TvdeMiQmPg3fxGT/udfx7ixo3JXUPFowVvxAVieAAECBAoJCN5CbJUNErw5+pYF76pt3Xzr9+N1m42OXd78V7mdtulowdumB2NbBAgQIPCaAoK3s24QwZs7r5YH76rtrVzZH0uWvrTWbocP2zh3BRWPFrwVH4DlCRAgQKCQgOAtxFbZIMGbo2958D7x1HMx/9qvx7fv/Ld45rlFa+327gWXxsgRw3NXUeFowVshvn83Oi8AACAASURBVKUJECBAoLCA4C1MV8lAwZtjb3nwfvKSL8R1X/lOHHfE/rHl68bGBhus+Zm7u7/1zbHhhhvkrqLC0YK3QnxLEyBAgEBhAcFbmK6SgYI3x97y4P27/T4ak9+9axx/1HtzO23T0YK3TQ/GtggQIEDgNQUEb2fdIII3d14tD95p0y+IrbcYFzNPOCy30zYdLXjb9GBsiwABAgQEbxfdA4I3d5gtD967/u2B+NgnPlP/AxNjR4/M7bYNRwveNjwUWyJAgACBhgLe4W1I1FYvELy542h58J5y1uVx6+0/fNVd+qW13AEaTYAAAQIEiggI3iJq1Y0RvDn7lgfvd7//4/jvXz/xqrs8eP+3x0ZDNsxdRYWjvcNbIb6lCRAgQKCwgOAtTFfJQMGbY2958Oa21/6jBW/7n5EdEiBAgMDaAoK3s+4KwZs7r5YHb39//2vucNCgQbkrqHi04K34ACxPgAABAoUEBG8htsoGCd4cfcuD9/gzLonaYw2v9uUZ3twBGk2AAAECBIoICN4iatWNEbw5+5YH7/d/+NP4zeNPr7XLeVd9Jf5q+23ikrOPjyGe4c2dotEECBAgQGA9BQTveoJV/HLBmzuAlgfvq23v5lu/H5+85Lq462vzBG/uDI0mQIAAAQLrLSB415us0gGCN8dfWfD+8rHHY69Dp8eNn50dO7xhm9xVVDjaM7wV4luaAAECBAoLCN7CdJUMFLw59kqCd+XK/vjS126Pcy6+Nm697lMxYctxuauocLTgrRDf0gQIECBQWEDwFqarZKDgzbG3PHjP+NRVsfDu+9bY5TPPLap/v8euE+PC2cflrqDi0YK34gOwPAECBAgUEhC8hdgqGyR4c/QtD95vfPeeeOSXv1ljl8OHDY23TNwx3vCnW+V23wajBW8bHIItECBAgMB6Cwje9SardIDgzfG3PHhz22vu6NqjFLXPBe7rG7zeEy96YXEsX7EiNh05Yo2xgne9KQ0gQIAAgTYQELxtcAjrsQXBux5Y63hpKcH70svL4pvfvSf+8+H/jiVLX4qtxm8We+62c2y9RXnP7tZCd/YFn68TnHnKEWtQ7Hv4afHwo79e458dN2W/OHbKfrF4ydKYPmd+3H7XK49l7LTDdjFvzvExdvTI+veCN3cDGk2AAAEC1QgI3mrci64qeIvKvTKu5cH71DPPx6HHzYlf/ebJ+oLDhm5cj8jaV+353dpzvK3+um3hvTHnomuj9uzwgftMWmfw7v2OXeoRvupr5IjhMWrkJnHl9d+IGxcsjGvnzYyhGw+JD8+YG9tOGB9nf/xIwdvqgzM/AQIECLRMQPC2jLYlEwveHGvLg7f2S2vfuuPeuOzcj9XfHd1oyIbxi1/+Ji74xxti4d0/if/zrc/WQ7KVX4uXvBS/e+HFmPvZG2PjjYasM3invG/POOBdb11rGwcePase5Ucfuk/9Z7V4Pmn2ZfHAHVdH7c8ie4e3lSdnbgIECBBolYDgbZVsa+YVvDnXlgfv2yafGPu8Y5c46UMHrbHTn/3XL+O9Uz8RN8yfFTtuv23uKgY4+qy518SKFSvWGbzDhw+N7V6/RWyx+ZjY5527xIQtN6/POnGvaTFn+lGr34l+8KFHYvIxs2PVn0QWvAPE9zICBAgQaCuBObOWt9V+bOa1BV43eiiihEDLg3f/I0+PN+3wZzH7lClrbPPe+34WR5x4XlsE76VX3xyD+wZHf3/E7T/4cTz6q8fjn688s/6M8Y67HRGXnXtiTNrlTfX9P/zIY7HvlJnxnRsuiPGbj/EOb+LmM5QAAQIEqhOYP3eD6ha38noL9A0etN5jDPh/Ai0P3gvnfzk+98Vv1oN357/+y/pzsT/66UMx/5qvx68ffypuv+mi2HCDvlLO5NXe4f3DxZctWx57HHJqHPbe3eOI9+9Vf4f3nBlTY/dJb66/zDu8pRyVRQgQIECgxQIeaWgxcJOn90hDDrTlwbtk6ctxwhmXxF3/9sAaOx09akRcfPbx8T/e+Oe5K1iP0QMJ3tp07/vQmTHpLX8dxx7+nqg9w1v7Zbaph+xdX8kzvOsB7qUECBAg0LYCgrdtj2adGxO8ufNqefCu2t59D/w8fv6LX0XtF8i22mKzeMub/6r+iQ1lfK1YsTJWrlwZcy6+NpYvXxGzT54SfX19MXjwoPjlY4/XP3KsFrVjNh0Zt91xb0w/Z35cc8lp8bc7vSGuuO6WuOmWO+uf0jBs6EYxbfqFPqWhjEOzBgECBAi0VEDwtpS36ZML3hxpy4P3P37+aNx6+w/jwH12jQlb/r/P3f3sFxbEZmNGxf57/V3uCgYw+stfvyPOvPCf1nhl7WPFap/KUAveKR87Lx5/8tnVP59+3MHxwcl71L9/cfHSOOWsy+N799xf/772C3bzzjkhxo0dVf/eL60N4AC8hAABAgTaTkDwtt2RvOaGBG/uvFoevDPPu7L+3OtNV5y1xl84u/7m78Y5F19byseSNSKq/VGK2mf01j4fuPaLaBv0rf1M8fOLXoza872r/uDEqjkFbyNdPydAgACBdhQQvO14Kq++J8GbO6+WB2/tr5jtu8f/Wv0M7KrtPvn0c7Hrez8WX/nc2bH9dlvnrqLC0YK3QnxLEyBAgEBhAcFbmK6SgYI3x97y4K39AtgO228Ts046fI2d1j6p4YPHfzIWXHNu/OmE8bmrqHC04K0Q39IECBAgUFhA8Bamq2Sg4M2xtzx4z7/0i3HNjbfF9ZedUX/+ta9vcDzx1HPxiU9/Ln7005/H3V//TGy4Yed+FqDgzd2ARhMgQIBANQKCtxr3oqsK3qJyr4xrefA+/7sXY/+jTq//UljtUxm2Gj82HvrFr+qLn3faMfHu3d+Su4KKRwveig/A8gQIECBQSEDwFmKrbJDgzdG3PHhr26v9MtgNX7sj/v1n/zeWLH0pttn6dfU/N/xX22+T230bjBa8bXAItkCAAAEC6y0geNebrNIBgjfHX0rw5rbY3qMFb3ufj90RIECAwLoFBG9n3RmCN3degjfn53N4k36GEyBAgEA1AoK3GveiqwreonKvjBO8OT/Bm/QznAABAgSqERC81bgXXVXwFpUTvDm534/2SENTGE1CgAABAiULCN6SwZPLCd4coHd4c37e4U36GU6AAAEC1QgI3mrci64qeIvKeYc3J+cd3qb4mYQAAQIEqhEQvNW4F11V8BaVE7w5OcHbFD+TECBAgEA1AoK3GveiqwreonKCNycneJviZxICBAgQqEZA8FbjXnRVwVtUTvDm5ARvU/xMQoAAAQLVCAjeatyLrip4i8oJ3pyc4G2Kn0kIECBAoBoBwVuNe9FVBW9ROcGbkxO8TfEzCQECBAhUIyB4q3EvuqrgLSoneHNygrcpfiYhQIAAgWoEBG817kVXFbxF5QRvTk7wNsXPJAQIECBQjYDgrca96KqCt6ic4M3JCd6m+JmEAAECBKoRELzVuBddVfAWlRO8OTnB2xQ/kxAgQIBANQKCtxr3oqsK3qJygjcnJ3ib4mcSAgQIEKhGQPBW4150VcFbVE7w5uQEb1P8TEKAAAEC1QgI3mrci64qeIvKCd6cnOBtip9JCBAgQKAaAcFbjXvRVQVvUTnBm5MTvE3xMwkBAgQIVCMgeKtxL7qq4C0qJ3hzcoK3KX4mIUCAAIFqBARvNe5FVxW8ReUEb05O8DbFzyQECBAgUI2A4K3GveiqgreonODNyQnepviZhAABAgSqERC81bgXXVXwFpUTvDk5wdsUP5MQIECAQDUCgrca96KrCt6icoI3Jyd4m+JnEgIECBCoRkDwVuNedFXBW1RO8ObkBG9T/ExCgAABAtUICN5q3IuuKniLygnenJzgbYqfSQgQIECgGgHBW4170VUFb1E5wZuTE7xN8TMJAQIECFQjIHircS+6quAtKid4c3KCtyl+JiFAgACBagQEbzXuRVcVvEXlBG9OTvA2xc8kBAgQIFCNgOCtxr3oqoK3qJzgzckJ3qb4mYQAAQIEqhEQvNW4F11V8BaVE7w5OcHbFD+TECBAgEA1AoK3GveiqwreonKCNycneJviZxICBAgQqEZA8FbjXnRVwVtUTvDm5ARvU/xMQoAAAQLVCAjeatyLrip4i8oJ3pyc4G2Kn0kIECBAoBoBwVuNe9FVBW9ROcGbkxO8TfEzCQECBAhUIyB4q3EvuqrgLSoneHNygrcpfiYhQIAAgWoEBG817kVXFbxF5QRvTk7wNsXPJAQIECBQjYDgrca96KqCt6ic4M3JCd6m+JmEAAECBKoRELzVuBddVfAWlRO8OTnB2xQ/kxAgQIBANQKCtxr3oqsK3qJygjcnJ3ib4mcSAgQIEKhGQPBW4150VcFbVE7w5uQEb1P8TEKAAAEC1QgI3mrci64qeIvK9WDwrlzZH/39/dHXN3gttdrPnnj62Rg7emRs0Ne31s8XvbA4lq9YEZuOHLHGz6aesCx3AkYTIECAAIEKBARvBeiJJQVvAi8iBvXXCrAHvmqXOfuCz9ev9MxTjljjiu/81/vjlLMuj8VLltb/+ayTDo+D9t2t/p9r/2z6nPlx+1331b/faYftYt6c4+thXPsSvD1w87hEAgQIdKGA4O2sQxW8ufPqieC9beG9Meeia+OZ5xbFgftMWiN4lyx9Od66//HxkSP3j0MPeEcsvPsnccIZ8+K2L346thq/WVx5/TfixgUL49p5M2PoxkPiwzPmxrYTxsfZHz9S8ObuPaMJECBAoEIBwVshfoGlBW8BtD8Y0hPBu3jJS/G7F16MuZ+9MTbeaMgawVt7d/fYv58b9337ihgyZMM6zbs+ML0ev4ce8M448OhZsceuE+PoQ/ep/6wWzyfNviweuOPqGDRokHd4c/ef0QQIECBQkYDgrQi+4LKCtyDc74f1RPCuIjpr7jWxYsWKNYL3ywsWxudvuDW++YXzV0t+dObFsc3W4+PkaQfFxL2mxZzpR9Wjt/b14EOPxORjZsfdCy6NkSOGC97c/Wc0AQIECFQkIHgrgi+4rOAtCCd4XxGoPbLwrTvujZuuOHO1ZO153k2GDY1ZJx8eO+52RFx27okxaZc31X/+8COPxb5TZsZ3brggxm8+RvDm7j+jCRAgQKAigflzN6ho5R5ftuBvTvX1DepxuNzle4d3AO/wnjNjauw+6c3e4c3da0YTIECAQBsJzJm1vI1200NbKditr9t0aA8hNf9Sez54Vz3D+5N/uTI23PCVf9vd4+BT44OTd1/9DO+eu+0cUw/Zu/4zz/A2/yY0IwECBAiUL+CRhvLNMyt6pCGj1yMfS7ZixcpYuXJlzLn42li+fEXMPnlK9PX1xeDBg6L2C20T9/pQTD/u4DhkHZ/ScMV1t8RNt9xZ/5SGYUM3imnTL/QpDbl7zmgCBAgQaAMBwdsGh7AeWxC864G1jpf2xDu8X/76HXHmhf+0xuXXPlbsgHe9tf7Pap+xW/tFtVVfp3/ssDh4v7fXv31x8dL6Z/R+757769/vuP22Me+cE2Lc2FH1730Ob+4GNJoAAQIEqhEQvNW4F11V8BaVe2VcTwTvQIhq7wL/9slnYtyYUasfbfjDcc8vejGWLVu++g9OrPqZ4B2IrtcQIECAQLsJCN52O5HX3o/gzZ2X4M35eYc36Wc4AQIECFQjIHircS+6quAtKucd3pzc70d7h7cpjCYhQIAAgZIFBG/J4MnlBG8O0Du8OT/v8Cb9DCdAgACBagQEbzXuRVcVvEXlvMObk/MOb1P8TEKAAAEC1QgI3mrci64qeIvKCd6cnOBtip9JCBAgQKAaAcFbjXvRVQVvUTnBm5MTvE3xMwkBAgQIVCMgeKtxL7qq4C0qJ3hzcoK3KX4mIUCAAIFqBARvNe5FVxW8ReUEb05O8DbFzyQECBAgUI2A4K3GveiqgreonODNyQnepviZhAABAgSqERC81bgXXVXwFpUTvDk5wdsUP5MQIECAQDUCgrca96KrCt6icoI3Jyd4m+JnEgIECBCoRkDwVuNedFXBW1RO8ObkBG9T/ExCgAABAtUICN5q3IuuKniLygnenJzgbYqfSQgQIECgGgHBW4170VUFb1E5wZuTE7xN8TMJAQIECFQjIHircS+6quAtKid4c3KCtyl+JiFAgACBagQEbzXuRVcVvEXlBG9OTvA2xc8kBAgQIFCNgOCtxr3oqoK3qJzgzckJ3qb4mYQAAQIEqhEQvNW4F11V8BaVE7w5OcHbFD+TECBAgEA1AoK3GveiqwreonKCNycneJviZxICBAgQqEZA8FbjXnRVwVtUTvDm5ARvU/xMQoAAAQLVCAjeatyLrip4i8oJ3pyc4G2Kn0kIECBAoBoBwVuNe9FVBW9ROcGbkxO8TfEzCQECBAhUIyB4q3EvuqrgLSoneHNygrcpfiYhQIAAgWoEBG817kVXFbxF5QRvTk7wNsXPJAQIECBQjYDgrca96KqCt6ic4M3JCd6m+JmEAAECBKoRELzVuBddVfAWlRO8OTnB2xQ/kxAgQIBANQKCtxr3oqsK3qJygjcnJ3ib4mcSAgQIEKhGQPBW4150VcFbVE7w5uQEb1P8TEKAAAEC1QgI3mrci64qeIvKCd6cnOBtip9JCBAgQKAaAcFbjXvRVQVvUTnBm5MTvE3xMwkBAgQIVCMgeKtxL7qq4C0qJ3hzcoK3KX4mIUCAAIFqBARvNe5FVxW8ReUEb05O8DbFzyQECBAgUI2A4K3GveiqgreonODNyQnepviZhAABAgSqERC81bgXXVXwFpUTvDk5wdsUP5MQIECAQDUCgrca96KrCt6icoI3Jyd4m+JnEgIECBCoRkDwVuNedFXBW1RO8ObkBG9T/ExCgAABAtUICN5q3IuuKniLygnenJzgbYqfSQgQIECgGgHBW4170VUFb1E5wZuTE7xN8TMJAQIECFQjIHircS+6quAtKid4c3KCtyl+JiFAgACBagQEbzXuRVcVvEXlBG9OTvA2xc8kBAgQIFCNgOCtxr3oqoK3qJzgzckJ3qb4mYQAAQIEqhEQvNW4F11V8BaVE7w5OcHbFD+TECBAgEA1AoK3GveiqwreonKCNycneJviZxICBAgQqEZA8FbjXnRVwVtUTvDm5ARvU/xMQoAAAQLVCAjeatyLrip4i8oJ3pyc4G2Kn0kIECBAoBoBwVuNe9FVBW9ROcGbkxO8TfEzCQECBAhUIyB4q3EvuqrgLSoneNdLbtELi2P5ihWx6cgRa4ybesKy9ZrHiwkQIECAQDsICN52OIWB70HwDtxqXa8c1N/f35+bovNH73v4afHwo79e40KOm7JfHDtlv1i8ZGlMnzM/br/rvvrPd9phu5g35/gYO3pk/XvB2/nn7woIECDQiwKCt7NOXfDmzkvwRkQtePd+xy6x5247r9YcOWJ4jBq5SVx5/TfixgUL49p5M2PoxkPiwzPmxrYTxsfZHz9S8ObuPaMJECBAoEIBwVshfoGlBW8BtD8YInh/H7xT3rdnHPCut66leeDRs2KPXSfG0YfuU//ZbQvvjZNmXxYP3HF1DBo0yDu8ufvPaAIECBCoSEDwVgRfcFnBWxDu98ME7++Dd/jwobHd67eILTYfE/u8c5eYsOXmdaKJe02LOdOPqkdv7evBhx6JycfMjrsXXBq1d4E90pC7AY0mQIAAAQIEGgtcefGGjV/kFa8qIHgj4tKrb47BfYOj9jTz7T/4cTz6q8fjn688M7beYlzsuNsRcdm5J8akXd5UR3z4kcdi3ykz4zs3XBDjNx8jeP0fFwECBAgQINByAcGbIxa8f+S3bNny2OOQU+Ow9+4eR7x/r/o7vOfMmBq7T3qzd3hz95rRBAgQIECAQEEBwVsQ7vfDBO86/N73oTNj0lv+Oo49/D1Re4a39stsUw/Zu/5Kz/DmbjijCRAgQIAAgfUXELzrb/aHI3o+eH/52OP1jxyrRe2YTUfGbXfcG9PPmR/XXHJa/O1Ob4grrrslbrrlzvqnNAwbulFMm36hT2nI3XNGEyBAgAABAuspIHjXE+yPXi54H3s8pnzsvHj8yWdX00w/7uD44OQ96t+/uHhpnHLW5fG9e+6vf7/j9tvGvHNOiHFjR9W/90truRvQaAIECBAgQKCxgOBtbPRar+j54K3h1P72xjPPLar/kYnaL6Jt0Ne3ltnzi16M2vO9q/7gxKoXCN7cDWg0AQIECBAg0FhA8DY2Erw5o9ccLXhbiGtqAgQIECBAoC4geHM3gnd4c34eaUj6GU6AAAECBAg0FhC8jY28w5sz8g5vC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5owEbwv9TE2AAAECBAg0FhC8jY0Eb85I8LbQz9QECBAgQIBAYwHB29hI8OaMBG8L/UxNgAABAgQINBYQvI2NBG/OSPC20M/UBAgQIECAQGMBwdvYSPDmjARvC/1MTYAAAQIECDQWELyNjQRvzkjwttDP1AQIECBAgEBjAcHb2Ejw5ozqoxe9sDiWr1gRm44cscZsU09Y1oTZTUGAAAECBAgQeHUBwZu7Owb19/f356bo7tGLlyyN6XPmx+133Ve/0J122C7mzTk+xo4eWf9e8Hb3+bs6AgQIECDQDgKCN3cKgreB35XXfyNuXLAwrp03M4ZuPCQ+PGNubDthfJz98SMFb+7eM5oAAQIECBAYoIDgHSDUq7xM8DbwO/DoWbHHrhPj6EP3qb/ytoX3xkmzL4sH7rg6Bg0a5B3e3P1nNAECBAgQIDAAAcE7AKTXeIngbeA3ca9pMWf6UfXorX09+NAjMfmY2XH3gktj5Ijhgjd3/xlNgAABAgQIDEBA8A4ASfAWQ6o93rzjbkfEZeeeGJN2eVN9kocfeSz2nTIzvnPDBTF+8zGCtxitUQQIECBAgMB6CAje9cBax0u9wzuAd3jPmTE1dp/05nW+w5vjN5oAAQIECBAgQKDVAoK3gXDtGd49d9s5ph6yd/2Vf/wMb6sPyPwECBAgQIAAAQI5AcHbwO+K626Jm265s/4pDcOGbhTTpl+4xqc05PiNJkCAAAECBAgQaLWA4G0g/OLipXHKWZfH9+65v/7KHbffNuadc0KMGzuq1WdT2fyv9kc2KtuQhV9ToPas+YqVK2ODvj5SbS6wbPmKeOrp52L0pn8SGw3ZsM13a3u1Pzb01DPPR//K/hg3dtPo6xsMhQCBDhUQvAM8uOcXvRjLli1f/QcnBjiso17W6I9sdNTF9NBmF3z77ph7xY1x+41ze+iqO+9Sa/9r0UVX3LR647VPfpl10pQY+SfDO+9iemDHN3zt9jhr7jWrr3TzzTaNS+YcX3/Tw1d7C8z97I1R+wz9f73lsviTTYa192btrjQBwVsadfsv1OiPbLT/FfTWDn/52ONx9Cn/EL/6zZNR+3/Ggre9z//GWxbG1luMizft8Gfx379+Io466fw46uC9Y8r79mzvjffo7mr/Ijlq5CbxtzttX/+z8qeceVksX74irpo7vUdFOuOyb771+3H6+Z+rb1bwdsaZlbVLwVuWdAes0+iPbHTAJfTUFlf9z623/+C+uPL6WwRvh53+GZ+6Kh77zZMCqkPOrfZo28qV/XHh7GM7ZMe9t81/+8nP4ti/vyjOOvWI+qOIgrf37oHXumLB635YLdDoj2ygak+BW2//YXz68i8J3vY8nnXuqvYs7x4HnxJ7v32XOHnaQR20897b6te/fVfU/qXyoV/8d1w4+7j4iz+b0HsIHXDFj/7q8ai9aXPRWR+JzcduGu85Yqbg7YBzK3OLgrdM7TZeayB/ZKONt9/TWxO8nXf8s/7h6vjmd38Y37j2vK7+BdjOO5m1d1x77vpHP30onnjq2Tj740fFzn/zF91wWV11Dc//7sU46EOz4/CD9oxD9n97/Nf/fUzwXB/EdgAAByxJREFUdtUJN+diBG9zHLtilto7vP7IRucdpeDtrDO77PNfjUs//9X40j/Oijf+hV+A6pTTm3/tgvjCP387vv/VeZ2y5Z7Z56rPx//g5D1iUEQ88/yiqD2D/b73vC0m7zMp/vLPX98zFi701QUEr7tjtYA/stGZN4Pg7Yxzqz3/ecE/3hBfXrAw/uniGbHDG7bpjI3bZV3g23f+nzhx1mfi/u9+zkcAttk98fAjj8V3f/Dj1buqfZTcdV/5TnzosHfH3m//n7HdNlu22Y5tpwoBwVuFepuu6Y9stOnBvMq2ao+h1H5r/Ft33Fv/WLLbrv90DBo8yP8zbtNjrP3meO03yP/x/JPjT18/fvUua5+w4TOU2+/Qau/E/6+d3xjbb7d1PP3s7+q/BDV0oyF+ybD9jmqtHXmkoQMOqYItCt4K0Nt1yV78IxvtehYD2deq/1L/w9e+e/e3xHmnHTOQ4V5TssAeB59a/wi5P/765hfOj9dvtXnJu7FcI4GZ510ZX/3WD1a/7G92/PM4b+YxsdX4zRoN9fOKBQRvxQfQpssL3jY9mCq31Qt/ZKNKX2sTINAZAi+/vCyeePq52GTY0Ppn8voiQKBzBQRv556dnRMgQIAAAQIECAxAQPAOAMlLCBAgQIAAAQIEOldA8Hbu2dk5AQIECBAgQIDAAAQE7wCQvIQAAQIECBAgQKBzBQRv556dnRMgQIAAAQIECAxAQPAOAMlLCBAgQIAAAQIEOldA8Hbu2dk5AQIECBAgQIDAAAQE7wCQvIQAAQIECBAgQKBzBQRv556dnRMgQIAAAQIECAxAQPAOAMlLCBAgQIAAAQIEOldA8Hbu2dk5AQIECBAgQIDAAAQE7wCQvIQAAQIECBAgQKBzBQRv556dnRMgQIAAAQIECAxAQPAOAMlLCBAgQIAAAQIEOldA8Hbu2dk5AQIECBAgQIDAAAQE7wCQvIQAAQIECBAgQKBzBQRv556dnRMgQIAAAQIECAxAQPAOAMlLCBAgQIAAAQIEOldA8Hbu2dk5AQIdIjDvqq/Evff9LM6ZMTUmbDmuvuv/fPi/4+y518T79t0t3r37WzrkSmyTAAECnSkgeDvz3OyaAIEOEnjqmedj/yNPj803Gx3XX3p6LFu+IiYfMyvGbDoyrpo7PTbcoK+DrsZWCRAg0HkCgrfzzsyOCRDoQIEf//tDcdhHPxmHHvDOeH7RC3H3vz0QN181J8aOHtmBV2PLBAgQ6CwBwdtZ52W3BAh0sMA/3XhbfOrSL9av4Ib5s2LH7bft4KuxdQIECHSOgODtnLOyUwIEOlzg+z/8aUybfmH9Kr75hfPj9Vtt3uFXZPsECBDoDAHB2xnnZJcECHS4wK9+82Tsf+QZseduO8ePfvqfsUFfX3zpH2fFsKEbdfiV2T4BAgTaX0Dwtv8Z2SEBAh0usGTpy3HIsWdFX19f/ZfWfvHL38R7p36i/ukM5512TIdfne0TIECg/QUEb/ufkR0SINDhArP+4eq46ZY749brzo8JW77yGMOXvnZ7/WPJzjzliDhwn0kdfoW2T4AAgfYWELztfT52R4AAAQIECBAgkBQQvElAwwkQIECAAAECBNpbQPC29/nYHQECBAgQIECAQFJA8CYBDSdAgAABAgQIEGhvAcHb3udjdwQIECBAgAABAkkBwZsENJwAAQIECBAgQKC9BQRve5+P3REgQIAAAQIECCQFBG8S0HACBAgQIECAAIH2FhC87X0+dkeAAAECBAgQIJAUELxJQMMJECBAgAABAgTaW0Dwtvf52B0BAgQIECBAgEBSQPAmAQ0nQIAAAQIECBBobwHB297nY3cECBAgQIAAAQJJAcGbBDScAAECBAgQIECgvQUEb3ufj90RIECAAAECBAgkBQRvEtBwAgQIECBAgACB9hYQvO19PnZHgAABAgQIECCQFBC8SUDDCRAgQIAAAQIE2ltA8Lb3+dgdAQIECBAgQIBAUkDwJgENJ0CAAAECBAgQaG8Bwdve52N3BAgQIECAAAECSQHBmwQ0nAABAgQIECBAoL0FBG97n4/dESBAgAABAgQIJAUEbxLQcAIECBAgQIAAgfYWELztfT52R4AAAQIECBAgkBQQvElAwwkQIECAAAECBNpbQPC29/nYHQECBAgQIECAQFJA8CYBDSdAgAABAgQIEGhvAcHb3udjdwQIECBAgAABAkkBwZsENJwAAQIECBAgQKC9BQRve5+P3REgQIAAAQIECCQFBG8S0HACBAgQIECAAIH2FhC87X0+dkeAAAECBAgQIJAUELxJQMMJECBAgAABAgTaW0Dwtvf52B0BAgQIECBAgEBSQPAmAQ0nQIAAAQIECBBobwHB297nY3cECBAgQIAAAQJJAcGbBDScAAECBAgQIECgvQUEb3ufj90RIECAAAECBAgkBQRvEtBwAgQIECBAgACB9hb4/wH+0hB+A2NsWwAAAABJRU5ErkJggg=="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#| caption: A histogram of samples from the discrete distribution.\n",
    "#| label: fig:discrete-samples\n",
    "px.histogram(x=[sample() for i in range(1000)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "EKOtYSxcJgYP",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Of course, GTSAM has all this machinery built-in, and Figure [2.3](#fig:discrete-samples-gtsam) is generated using GTSAM's `DiscreteDistribution.sample` method. Notice that even for a 1000 samples the histograms can look noticeably different because of the variability in sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "r4w6AyHQXUwE",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAH0CAYAAADfWf7fAAAgAElEQVR4Xu3deXxV1bnw8ScJBMLQhMFgQFBKKy1F7b0V3ou3FXBgEERAUAGVgIAoCjIZCiqDQcABlBQUwaGiVoTWtqCIVxkcqNKrqNeP9WLxghUtIpOBJECG97M3JjVMJ1nPOWuds/cvfzXmPHvt/V277a+nOydJZWVlZcIXAggggAACCCCAAAIBFUgieAO6s1wWAggggAACCCCAgC9A8HIjIIAAAggggAACCARagOAN9PZycQgggAACCCCAAAIEL/cAAggggAACCCCAQKAFCN5Aby8XhwACCCCAAAIIIEDwcg8ggAACCCCAAAIIBFqA4A309nJxCCCAAAIIIIAAAgQv9wACCCCAAAIIIIBAoAUI3kBvLxeHAAIIIIAAAgggQPByDyCAAAIIIIAAAggEWoDgDfT2cnEIIIAAAggggAACBC/3AAIIIIAAAggggECgBQjeQG8vF4cAAggggAACCCBA8HIPIIAAAggggAACCARagOAN9PZycQgggAACCCCAAAIEL/cAAggggAACCCCAQKAFCN5Aby8XhwACCCCAAAIIIEDwcg8ggAACCCCAAAIIBFqA4A309nJxCCCAAAIIIIAAAgQv9wACCCCAAAIIIIBAoAUI3kBvLxeHAAIIIIAAAgggQPByDyCAAAIIIIAAAggEWoDgDfT2cnEIIIAAAggggAACBC/3AAIIIIAAAggggECgBQjeQG8vF4cAAggggAACCCBA8HIPIIAAAggggAACCARagOAN9PZycQgggAACCCCAAAIEL/cAAggggAACCCCAQKAFCN5Aby8XhwACCCCAAAIIIEDwcg8ggAACCCCAAAIIBFqA4A309nJxCCCAAAIIIIAAAgQv9wACCCCAAAIIIIBAoAUI3kBvLxeHAAIIIIAAAgggQPByDyCAAAIIIIAAAggEWoDgDfT2cnEIIIAAAggggAACBC/3AAIIIIAAAggggECgBQjeQG8vF4cAAggggAACCCBA8HIPIIAAAggggAACCARagOAN9PZycQgggAACCCCAAAIEL/cAAggggAACCCCAQKAFCN5Aby8XhwACCCCAAAIIIEDwcg8ggAACCCCAAAIIBFqA4A309nJxCCCAAAIIIIAAAgQv9wACCCCAAAIIIIBAoAUI3kBvLxeHAAIIIIAAAgggQPByDyCAAAIIIIAAAggEWoDgDfT2cnEIIIAAAggggAACBC/3AAIIIIAAAggggECgBQjeQG8vF4cAAggggAACCCBA8HIPIIAAAggggAACCARagOAN9PZycQgggAACCCCAAAIEL/cAAggggAACCCCAQKAFCN5Aby8XhwACCCCAAAIIIEDwcg8ggAACCCCAAAIIBFqA4A309nJxCCCAAAIIIIAAAgQv9wACCCCAAAIIIIBAoAUI3kBvLxeHAAIIIIAAAgggQPByDyCAAAIIIIAAAggEWoDgDfT2cnEIIIAAAggggAACBK/yHvhyd6HyCIwjgAACCCCAAAKnFmjaKA0ihUBogvfQ4SOya/c+qZNWWxpm1K82Wf6BAikuKZEG6ZVnCd5qUzKAAAIIIIAAAtUUIHirCXbMy0MRvHfMeUxeWP1GxaX/+zlnS17uaMlIr+f/s16DJ8vW7V9WohmV3Vtuzu4tBYVFkpO7SNa+tdn/+bltWvmzjRum+98TvLobkGkEEEAAAQQQiCxA8EY2OtUrQhG8i5aulF+2P0fObtVcvtr5jQwalSvX9esiI669vCJ4e1zSQbp1bl9hlV6/rh/ES559UZavXC9L86ZIWu1UuWnSPGnZIkvuvn0owau795hGAAEEEEAAgSoKELxVhDrJy0IRvN+/9iNHiuWi/mPl1qF95apenSuCN/vqbtL3sguPY+o3fKp07dROhg/q6f9szfpNMm7aQvlo3ROSlJTEO7y6+49pBBBAAAEEEKiCAMFbBaRTvCQ0wXv48BF5/LnVsuHtD+S0Rulyz6ThUq/u0QfAvUca6tZNk1ZnNpWmTRpJz0s7SItmTfyftes+UnJzbvCj1/v6eMs26T9immxcuUC8d4F5pEF3AzKNAAIIIIAAApEFCN7IRqd6RWiCt7DosEyZvUQ++ft2yWzcQGZNHiFZmQ19mwVPvCDJKclSViay9s33ZPsXO+X3S6ZL86aZ0rbzEFk4a6x07HCe/9qt23ZIr+wp8uqyBySrSSM5WFSs2wGmTyiQhAsCCCCAAAIBFCgzvKa6tWsYTjLmCYQmeMu3u6ysTIZPvF9OP62h/87tsV/eIw9dB06U667sIkOu6e6/wztz0jDp0vF8/6XHvsO778Bh7qQYCJj+B0IMToVDIoAAAgggEDUB0zd0MuqlRu0cwnig0AWvt8n3zH9aPvv8K1ly/8QT7vnVN06Xjhf8XG4efIV4z/B6v8w2bGAP/7U8wxvGf5twzQgggAACCLgV4JEGnX/gg/fAwUJ59OmV0qf7r+SMppn+O7TDxt/nB+yN110un+/Y6X/kmBe1jRqky5p1myRn5iJ5av5k+cW5Z8viZ1bJilUb/E9pqJNWS0bmzOVTGnT3HNMIIIAAAgggUE0BgreaYMe8PPDBe7CgSAaPmSV/+3R7xaX37vZLuWvcYKmVWtMP3uzbZsvOXXsrfp4zaoBc37+r/703P2HGw/L62x/437dt3VLyZo6RzMYZ/vf80pruBmQaAQQQQAABBCILELyRjU71isAHb/nFe+G6e+9+adwww3+n9vtf3nO9e/bl+39kwvtFtBopKceZ7c8/KN7zveV/cKL8BQSv7gZkGgEEEEAAAQQiCxC8kY0IXp3RKacJ3hjicmgEEEAAAQQQ8AUIXt2NEJp3eHVMJ58meGMly3ERQAABBBBAoFyA4NXdCwSvzo9neJV+jCOAAAIIIIBAZAGCN7LRqV5B8Or8CF6lH+MIIIAAAgggEFmA4I1sRPDqjE45zSMNMcTl0AgggAACCCDgCxC8uhuBd3h1frzDq/RjHAEEEEAAAQQiCxC8kY14h1dnxDu8MfSzfei7ZvC3yG2ba9abcVexZpxZBBBAIDACBK9uK3mHV+fHO7xKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr28vQBO+hw0dk1+59UiettjTMqH+cWmlpmXy9e680bpguNVJSjvt5/oECKS4pkQbplWe/3F2o2wGmrQoQvFa51YsRvGpCDoAAAgERIHh1GxmK4L1jzmPywuo3KqT+/ZyzJS93tGSk1/P/2Ya/fCATZjwsBYVF/vdTxw2Wq3p19v+1989ychfJ2rc2+9+f26aVP+uFsfdF8OpuQNvTBK9tcd16BK/Oj2kEEAiOAMGr28tQBO+ipSvll+3PkbNbNZevdn4jg0blynX9usiIay+XwqLDcmGf0XLL0D4yqO8lsn7j+zLmzjxZ87v75Iys02TJsy/K8pXrZWneFEmrnSo3TZonLVtkyd23DyV4dfeek2mC1wm78aIErzEdgwggEDABgle3oaEI3u8THTlSLBf1Hyu3Du3rv4vrvbt786/nyeZXFktqak3/pZddm+PH76C+l0q/4VOla6d2MnxQT/9na9ZvknHTFspH656QpKQk3uHV3X/Wpwle6+SqBQleFR/DCCAQIAGCV7eZoQnew4ePyOPPrZYNb38gpzVKl3smDZd6ddPk+ZXr5cllq+Wlp+dUSN465SE5q3mWjB95lbTrPlJyc27wo9f7+njLNuk/YppsXLlA0uvXla/3HdLtANNWBSZPTba6HovpBO6ZXqo7ANMIIIBAQAQyM2oF5ErcXEZogtd7dGHK7CXyyd+3S2bjBjJr8gjJymzoP7Lw8rpNsmLx9Iod8J7nrVcnTaaOHyxtOw+RhbPGSscO5/k/37pth/TKniKvLntAspo0kuIS/gvZza1rturIcSVmg0w5EXhk7vG/QOrkRFgUAQQQcCxQI4U3bDRbEJrgLUcqKyuT4RPvl9NPa+i/c1uVd3hnThomXTqef8J3ePmlNc3tZ3+WRxrsm2tW5JEGjR6zCCAQJAEeadDtZuiC1+O6Z/7T8tnnX8mS+ydWPMP7/n8tkZo1a/iaXQdMlOv7d6l4hrdb5/YybGAP/2c8w6u74VxPE7yud6B66xO81fPi1QggEFwBgle3t4EP3gMHC+XRp1dKn+6/kjOaZvrP4A4bf58fsDded7kUFB6Sdt1vlJxRA2TgCT6lYfEzq2TFqg3+pzTUSaslI3Pm8ikNunvO6TTB65S/2osTvNUmYwABBAIqQPDqNjbwwXuwoEgGj5klf/t0e4VU726/lLvGDZZa330qg/cZu94vqpV/3XHbdTKg98X+t96890zv629/4H/ftnVLyZs5RjIbZ/jf80iD7ga0PU3w2hbXrUfw6vyYRgCB4AgQvLq9DHzwlvN44bp7735p3DDDf6f22K+SklL55649ktkoo+LRhu+/Zn/+QfE+0qz8D06U/4zg1d2AtqcJXtviuvUIXp0f0wggEBwBgle3l6EJXh3TyacJ3ljJxua4BG9sXGN1VII3VrIcFwEEEk2A4NXtGMGr8+ORBqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXBK/Oj+BV+tkeJ3hti+vWI3h1fkwjgEBwBAhe3V4SvDo/glfpZ3uc4LUtrluP4NX5MY0AAsERIHh1e0nw6vwIXqWf7XGC17a4bj2CV+fHNAIIBEeA4NXtJcGr8yN4lX62xwle2+K69QhenR/TCCAQHAGCV7eXoQnewqLDsnfft3J6ZiNJTk6qtlr+gQIpLimRBun1K81+ubuw2sdiwJ0AwevO3mRlgtdEjRkEEAiiAMGr29VQBO+tUx6StW9t9qUaZtSX3t1+JeNHXlUh12vwZNm6/ctKkqOye8vN2b2loLBIcnIXVcyf26aV5OWOlsYN0/3XE7y6G9D2NMFrW1y3HsGr82MaAQSCI0Dw6vYyFMH7m8dfkC6d2kmLZpny9rsfy6jJD8pzD98l5/z0h76eF7w9Lukg3Tq3r9BMr19XMtLryZJnX5TlK9fL0rwpklY7VW6aNE9atsiSu28fSvDq7j0n0wSvE3bjRQleYzoGEUAgYAIEr25DQxG8xxJd1H+sXHPFRTLi2ssrgjf76m7S97ILj9PsN3yqdO3UToYP6un/bM36TTJu2kL5aN0TkpSUxDu8uvvP+jTBa51ctSDBq+JjGAEEAiRA8Oo2M3TBu/2LnXLZtTmycNZY6djhvIrgrVs3TVqd2VSaNmkkPS/tIC2aNfF/1q77SMnNucGPXu/r4y3bpP+IabJx5QLx3gXmkQbdDWh7muC1La5bj+DV+TGNAALBESB4dXsZquA9WFAk196SK/Xq1pEnH5wkKSnJvt6CJ16Q5JRkKSsTWfvme+JF8e+XTJfmTTOlbechleJ467Yd0it7iry67AHJatJIig6X6HaAaasCt0wstboei+kEfnPf0X+P8oUAAgiEXaB2akrYCVTXH5rg9T6lYcyd8+WfX++Rp+ZP9p/PPdHXkSPF0nXgRLnuyi4y5Jru/ju8MycNky4dz/dffuw7vHvyD6s2gGG7ArffUf1P6LB7hqz2fYF7c8sAQQABBBDwfum+fioOCoFQBO+3Bwpk9B3zpbDwkCy6d/xJY7fc8eobp0vHC34uNw++QrxneL1fZhs2sIf/Y57hVdxtcTDKIw1xsAnVOAUeaagGFi9FAIFAC/BIg257Ax+8BYWH5JqR0/3P0J03/RapVzfNF0tOTpaszIby+Y6d/keOeVHbqEG6rFm3SXJmLvLfBf7FuWfL4mdWyYpVG/xPaaiTVktG5szlUxp095zTaYLXKX+1Fyd4q03GAAIIBFSA4NVtbOCDd+euveJ9KsOxX97n8b7xxzw/eLNvmy3e68q/ckYNkOv7d/W/9Z77nTDjYXn97Q/879u2bil5M8dIZuMM/3t+aU13A9qeJnhti+vWI3h1fkwjgEBwBAhe3V4GPnirwlNWViZ79uX7f2TC+0W0GinHPxi+P/+geM/3lv/BifLjErxVEY6f1xC88bMXVTkTgrcqSrwGAQTCIEDw6naZ4NX58Q6v0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFdesRvDo/phFAIDgCBK9uLwlenR/Bq/SzPU7w2hbXrUfw6vyYRgCB4AgQvLq9JHh1fgSv0s/2OMFrW1y3HsGr82MaAQSCI0Dw6vaS4NX5EbxKP9vjBK9tcd16BK/Oj2kEEAiOAMGr20uCV+dH8Cr9bI8TvLbFWS9MAvwPlDDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YkTvDo/glfpZ3uc4LUtznphEiB4w7TbXKttAYJXJ07w6vwIXqWf7XGC17Y464VJgOAN025zrbYFCF6dOMGr8yN4lX62xwle2+KsFyYBgjdMu8212hYgeHXiBK/Oj+BV+tkeJ3hti7NemAQI3jDtNtdqW4Dg1YmHJngLiw7L3n3fyumZjSQ5Oek4tdLSMvl6915p3DBdaqSkHPfz/AMFUlxSIg3S61f62Ze7C3U7wLRVAYLXKjeLhUyA4A3ZhnO5VgUIXh13KIL31ikPydq3NvtSDTPqS+9uv5LxI6+qkNvwlw9kwoyHpaCwyP9nU8cNlqt6dfb/tffPcnIXVcyf26aV5OWO9sPY+yJ4dTeg7WmC17Y464VJgOAN025zrbYFCF6deCiC9zePvyBdOrWTFs0y5e13P5ZRkx+U5x6+S8756Q/Fe+f3wj6j5ZahfWRQ30tk/cb3ZcydebLmd/fJGVmnyZJnX5TlK9fL0rwpklY7VW6aNE9atsiSu28fSvDq7j0n0wSvE3YWDYkAwRuSjeYynQgQvDr2UATvsUQX9R8r11xxkYy49nLx3t29+dfzZPMriyU1tab/0suuzfHjd1DfS6Xf8KnStVM7GT6op/+zNes3ybhpC+WjdU9IUlIS7/Dq7j/r0wSvdXIWDJEAwRuizeZSrQsQvDry0AXv9i92+kG7cNZY6djhPHl+5Xp5ctlqeenpORWS3iMQZzXP8h97aNd9pOTm3OBHr/f18ZZt0n/ENNm4coGk169L8OruP+vTBK91chYMkQDBG6LN5lKtCxC8OvJQBe/BgiK59pZcqVe3jjz54CRJSUn2H1l4ed0mWbF4eoWk9zxvvTppMnX8YGnbeUhFHHsv2Lpth/TKniKvLntAspo0kt35h3U7wLRVgZw7jv+FRasnwGIIBFhgTm5ZgK+OS0PArUCj+qluTyDBVw9N8HrP6o65c7788+s98tT8yZKRXs/fuqq8wztz0jDp0vH8E77De+hwSYLfAuE6/VETS8N1wZDzWwAAACAASURBVFwtAhYFFtyXbHE1lkIgXAK1Uo//BKlwCeiuNhTB++2BAhl9x3wpLDwki+4dXxG7Hl35M7zv/9cSqVmzhq/ZdcBEub5/l4pneLt1bi/DBvbwf8YzvLobzvU0jzS43gHWD7IAjzQEeXe5NtcCPNKg24GYB++mzZ9I+g/qSutWzSud6a7d+/xPTOh+8f874efe6i7rX9MFhYfkmpHT/c/QnTf9FqlXN83/YXJysmRlNhTv5+263yg5owbIwBN8SsPiZ1bJilUb/E9pqJNWS0bmzOVTGqK1OQ6OQ/A6QGfJ0AgQvKHZai7UgQDBq0OPefB6vwDWpvVZctP1V1Q60y//+Y1ces0EWfXULD8gY/W1c9de8T6V4dgv7/N43/hjnv+Pvc/o9c6z/OuO266TAb0v9r/1nvv1nul9/e0P/O/btm4peTPHSGbjDP97Poc3VjsXm+MSvLFx5agIeAIEL/cBArETIHh1ts6Ct/zTDlY/M0daNGuiu4ooTJeUlMo/d+2RzEYZFY82fP+w+/MPypEjxRV/cKL8ZwRvFPAtHoLgtYjNUqETIHhDt+VcsEUBgleHHbPgnXTPo7Jvf768++Gn/l83a9ni9IozPXy4WN7Z/Df56Y/PrPTpCLpLcTNN8LpxN12V4DWVYw6ByAIEb2QjXoGAqQDBayp3dC5mwXvnvY/L/vwDsvl/PpX69erIj1o2qzjT2qmp0u7ffiId/+PnFY8G6C7D3TTB687eZGWC10SNGQSqJkDwVs2JVyFgIkDwmqj9ayZmwVu+xAur35DTT2soHc7/me5M43Sa4I3TjTnJaRG8ibVfnG1iCRC8ibVfnG1iCRC8uv2KefCWn15paZkUFh067mzr1qmtuwLH0wSv4w2o5vIEbzXBeDkCCARWgP+BklhbS/Dq9ivmwfv1N/tk0dI/yysb/ip79uUfd7blf6JXdxnupgled/YmKxO8JmrMIIBAEAUI3sTaVYJXt18xD9575j8tz/zhVRk1pI80O72x1KhR+S+FdLnw/BN+KoLusuxNE7z2rKOxEsEbDUWOgQACQRAgeBNrFwle3X7FPHh/1ftW6X95Jxl9w5W6M43TaYI3TjfmJKdF8CbWfnG2CCAQOwGCN3a2sTgywatTjXnwjsx5QJo3zZQpY67TnWmcThO8cboxBG9ibQxniwAC1gUIXuvkqgUJXhVf7D6WrPy03vrrR3LbXb8R7w9MNG6YrjvbOJwmeONwU05xSrzDm1j7xdkigEDsBAje2NnG4sgEr0415u/wen+Wd/Xad056lvzSmm4Dma6eAMFbPS9ejQACwRUgeBNrbwle3X7FPHhfe+M9+ceXX5/0LAf0uVhqpdbUXYXDad7hdYhvsDTBa4DGCAIIBFKA4E2sbSV4dfsV8+DVnV78Tw8bcyT+T5IzRAABBBBA4BgBgjexbgmCV7dfMQ/esrKyU55hUlKS7gocTxO8jjeA5RFAAAEEjAQIXiM2Z0MEr44+5sE7+s754j3WcLKvRH+Gl+DV3YBMI4AAAgi4ESB43bibrkrwmsodnYt58L7xzofy1c7dx51l3uN/kJ+1Pkvm3z1aUhP4GV6CV3cDMo0AAggg4EaA4HXjbroqwWsqZyl4T3Z6L6x+Q+6Z/4y89ac8gle3h0wjgAACCCBQbQGCt9pkTgcIXh1/zN/hPdnpfb5jp3QflCPLH50mbc4+S3cVDqd5h9chPksjgAACCBgLELzGdE4GCV4du5PgLS0tk+f+tFZmPrRUVj9zr7Rolqm7CofTBK9DfJZGAAEEEDAWIHiN6ZwMErw69pgH7533Pi7rN26udJZ79uX733ft1E7mThuluwLH0wSv4w1geQQQQAABIwGC14jN2RDBq6OPefC++Nrbsu3zryqdZd06aXJBu7Zy9g/P0J19HEwTvHGwCZwCAggggEC1BQjeapM5HSB4dfwxD17d6cX/NMEb/3vEGSKAAAIIHC9A8CbWXUHw6vbLSvAeOnxEXnrtbfnfrf+QwqJDckbWadKtc3tp3jRxn90tZyd4dTcg0wgggAACbgQIXjfupqsSvKZyR+diHrzf7Nkvg0blyhdf7fIXrJNWWwoKi/x/7T2/6z3Hm8hfBG8i7x7njgACCIRXgOBNrL0neHX7FfPg9X5p7eV1m2ThrNvk3DatpFZqTfns86/kgUeWyfqN78t/v/yopNVO1V2Fw2mC1yE+SyOAAAIIGAsQvMZ0TgYJXh17zIP3ov5jpeclHWTcjVdVOtNP/v65XDnsLlm2aKq0bd1SdxUOpwleh/gsjQACCCBgLEDwGtM5GSR4dewxD94+Q++Q89r8SKZNyK50pps2fyJDxs4meHX7xzQCCCCAAAJGAgSvEZuzIYJXRx/z4J276Hl57Hcv+cHb/uc/lYz0evLuh1tk0VN/li93fiNrVzwoNWuk6K7C4TTv8DrEZ2kEEEAAAWMBgteYzskgwatjj3nwFhYdljF3zpe3/vpRpTNtmFFfHrp7tPz7OT/WXYHjaYLX8QawPAIIIICAkQDBa8TmbIjg1dHHPHjLT2/zR5/Kp599IQWFh+SMpqfJBef/zP/EhkT/IngTfQc5fwQQQCCcAgRvYu07wavbr5gH798+3S6r174j/Xp2khbN/vW5u48+vVJOa5Qhfbr/SncFjqcJXscbwPIIIIAAAkYCBK8Rm7MhgldHH/PgnTJ7iXy8ZZusWDxDUlKSK8722Rdek5kPLeVjyXT7xzQCCCCAAAJGAgSvEZuzIYJXRx/z4O01eLL06vqfMmxgj0pnumv3Pul05W3yh8fultatmuuuwuE07/A6xGdpBBBAAAFjAYLXmM7JIMGrY4958F5943Rp0/osmTpucKUz9T6p4frR98jKp2bJD1tk6a7C4TTB6xCfpRFAAAEEjAUIXmM6J4MEr4495sE7Z8Hv5Knla+TZhXf6f2DCe6zh62/2yV33PSbvfvipbPzzb6RmzRq6q3A4TfA6xGdpBBBAAAFjAYLXmM7JIMGrY4958O7/9qD0ueEO2blrr/+pDGdkNZYtn33hn/XsySPk8i4X6K7A8TTB63gDWB4BBBBAwEiA4DViczZE8OroYx683ukVFBbJsj+tk//55P+ksOiQnNX8dP/PDf+s9Vm6s4+DaYI3DjaBU0AAAQQQqLYAwVttMqcDBK+O30rw6k4xvqcJ3vjeH84OAQQQQODEAgRvYt0ZBK9uvwhenZ8QvEpAxhFAAAEEnAgQvE7YjRcleI3p/EGCV+dH8Cr9GEcAAQQQcCNA8LpxN12V4DWVOzpH8Or8CF6lH+MIIIAAAm4ECF437qarErymciEM3tLSMikrK6v0F9+qypd/oECKS0qkQXr9SiM80lBVQV6HAAIIIBBPAgRvPO1G5HMheCMbneoVoXmH1wvdaQ886VtMnzCkkon31+C2bv+y0j8bld1bbs7u7X/CRE7uIln71mb/5+e2aSV5uaOlccN0/3uCV3cDMo0AAggg4EaA4HXjbroqwWsqF6J3eNes3yS5Dy6VPfvypV/PjicM3h6XdJBundtXaKbXrysZ6fVkybMvyvKV62Vp3hRJq50qN02aJy1bZMndtw8leHX3HtMIIIAAAg4FCF6H+AZLE7wGaN8bCcU7vAWFh+TbAwdl3qPLpXat1BMGb/bV3aTvZRcep9lv+FTp2qmdDB/U0/+ZF8/jpi2Uj9Y9IUlJSbzDq7v/mEYAAQQQcCRA8DqCN1yW4DWE+24sFMFbTjRj3lNSUlJywuCtWzdNWp3ZVJo2aSQ9L+0gLZo18cfadR8puTk3+NHrfX28ZZv0HzFNNq5cIN67wDzSoLsBmUYAAQQQcCNA8LpxN12V4DWVOzpH8IrIgidekOSUZCkrE1n75nuy/Yud8vsl06V500xp23mILJw1Vjp2OM8H27pth/TKniKvLntAspo0Inh19x/TCCCAAAKOBB6Zm+JoZZY1EaiRkmwyxsx3AgTvMbfCkSPF0nXgRLnuyi4y5Jru/ju8MycNky4dz+cdXv5tgwACCCAQGIF7ppcG5lrCcCGZGbXCcJkxu0aC9wS0V984XTpe8HO5efAV4j3D6/0y27CBPfxX8gxvzO5FDowAAgggYFGARxosYkdhKR5p0CGGInhLSkqltLRUch9aKsXFJTJtfLakpKRIcnKSfL5jp/+RY17UNmqQLmvWbZKcmYvkqfmT5Rfnni2Ln1klK1Zt8D+loU5aLRmZM5dPadDdc0wjgAACCMSBAMEbB5tQjVMgeKuBdYKXhiJ4n//zOpk+97eVLt/7WDHvUxm84M2+bbbs3LW34uc5owbI9f27+t8fLCiSCTMeltff/sD/vm3rlpI3c4xkNs7wv+eX1nQ3INMIIIAAAm4ECF437qarErymckfnQhG8kYi8P0rhfUav90cmvF9Eq5Fy/IP8+/MPivd8b/kfnCg/JsEbSZefI4AAAgjEowDBG4+7cvJzInh1+0Xw6vx4h1fpxzgCCCCAgBsBgteNu+mqBK+pHO/w6uS+m+Yd3qgwchAEEEAAAcsCBK9lcOVyBK8OkHd4dX68w6v0YxwBBBBAwI0AwevG3XRVgtdUjnd4dXK8wxsVPw6CAAIIIOBGgOB14266KsFrKkfw6uQI3qj4cRAEEEAAATcCBK8bd9NVCV5TOYJXJ0fwRsWPgyCAAAIIuBEgeN24m65K8JrKEbw6OYI3Kn4cBAEEEEDAjQDB68bddFWC11SO4NXJEbxR8eMgCCCAAAJuBAheN+6mqxK8pnIEr06O4I2KHwdBAAEEEHAjQPC6cTddleA1lSN4dXIEb1T8OAgCCCCAgBsBgteNu+mqBK+pHMGrkyN4o+LHQRBAAAEE3AgQvG7cTVcleE3lCF6dHMEbFT8OggACCCDgRoDgdeNuuirBaypH8OrkCN6o+HEQBBBAAAE3AgSvG3fTVQleUzmCVydH8EbFj4MggAACCLgRIHjduJuuSvCayhG8OjmCNyp+HAQBBBBAwI0AwevG3XRVgtdUjuDVyRG8UfHjIAgggAACbgQIXjfupqsSvKZyBK9OjuCNih8HQQABBBBwI0DwunE3XZXgNZUjeHVyBG9U/DgIAggggIAbAYLXjbvpqgSvqRzBq5MjeKPix0EQQAABBNwIELxu3E1XJXhN5QhenRzBGxU/DoIAAggg4EaA4HXjbroqwWsqR/Dq5AjeqPhxEAQQQAABNwIErxt301UJXlM5glcnR/BGxY+DIIAAAgi4ESB43bibrkrwmsoRvDo5gjcqfhwEAQQQQMCNAMHrxt10VYLXVI7g1ckRvFHx4yAIIIAAAm4ECF437qarErymcgSvTo7gjYofB0EAAQQQcCNA8LpxN12V4DWVI3h1cgRvVPw4CAIIIICAGwGC14276aoEr6kcwauTI3ij4sdBEEAAAQTcCBC8btxNVyV4TeUIXp0cwRsVPw6CAAIIIOBGgOB14266KsFrKkfw6uQI3qj4cRAEEEAAATcCBK8bd9NVCV5TOYJXJ0fwRsWPgyCAAAIIuBEgeN24m65K8JrKEbw6OYI3Kn4cBAEEEEDAjQDB68bddFWC11SO4NXJEbxR8eMgCCCAAAJuBAheN+6mqxK8pnIEr06O4I2KHwdBAAEEEHAjQPC6cTddleA1lSN4dXIEb1T8OAgCCCCAgBsBgteNu+mqBK+pHMGrkyN4o+LHQRBAAAEE3AgQvG7cTVcleE3lCF6dHMEbFT8OggACCCDgRoDgdeNuuirBaypH8OrkCN6o+HEQBBBAAAE3AgSvG3fTVQleUzmCVydH8EbFj4MggAACCLgRIHjduJuuSvCayhG8OjmCNyp+HAQBBBBAwI0AwevG3XRVgtdULoTBW1paJmVlZZKSknycmvezr3fvlcYN06VGSspxP88/UCDFJSXSIL1+pZ8NG3NEtwNMI4AAAggg4ECA4HWArliS4FXgiUhSmVeAIfjyLnPaA0/6Vzp9wpBKV7zhLx/IhBkPS0Fhkf/Pp44bLFf16uz/a++f5eQukrVvbfa/P7dNK8nLHe2HsfdF8Ibg5uESEUAAgQAKELyJtakEr26/QhG8a9ZvktwHl8qeffnSr2fHSsFbWHRYLuwzWm4Z2kcG9b1E1m98X8bcmSdrfnefnJF1mix59kVZvnK9LM2bImm1U+WmSfOkZYssufv2oQSv7t5jGgEEEEDAoQDB6xDfYGmC1wDteyOhCN6CwkPy7YGDMu/R5VK7Vmql4PXe3b351/Nk8yuLJTW1pk9z2bU5fvwO6nup9Bs+Vbp2aifDB/X0f+bF87hpC+WjdU9IUlIS7/Dq7j+mEUAAAQQcCRC8juANlyV4DeG+GwtF8JYTzZj3lJSUlFQK3udXrpcnl62Wl56eUyF565SH5KzmWTJ+5FXSrvtIyc25wY9e7+vjLduk/4hpsnHlAkmvX5fg1d1/TCOAAAIIOBIgeB3BGy5L8BrCEbxHBbxHFl5et0lWLJ5eIek9z1uvTppMHT9Y2nYeIgtnjZWOHc7zf7512w7plT1FXl32gGQ1aUTw6u4/phFAAAEEHAncPzMUv8LjSDf6y2bUS43+QUN0RN7hrcI7vDMnDZMuHc/nHd4Q/RuDS0UAAQSCLvDQnKSgX2Kgrq9u7RqBuh7bFxP64C1/hvf9/1oiNWsevZm6Dpgo1/fvUvEMb7fO7WXYwB7+z3iG1/YtynoIIIAAArEQ4JGGWKjG7pg80qCzDUXwlpSUSmlpqeQ+tFSKi0tk2vhsSUlJkeTkJPF+oa1d9xslZ9QAGXiCT2lY/MwqWbFqg/8pDXXSasnInLl8SoPunmMaAQQQQCAOBAjeONiEapwCwVsNrBO8NBTB+/yf18n0ub+tdPnex4r1vexC/595n7Hr/aJa+dcdt10nA3pf7H97sKDI/4ze19/+wP++beuWkjdzjGQ2zvC/53N4dTcg0wgggAACbgQIXjfupqsSvKZyR+dCEbxVIfLeBf7nrj2S2Sij4tGG78/tzz8oR44UV/zBifKfEbxV0eU1CCCAAALxJkDwxtuOnPp8CF7dfhG8Oj/e4VX6MY4AAggg4EaA4HXjbroqwWsqxzu8OrnvpnmHNyqMHAQBBBBAwLIAwWsZXLkcwasD5B1enR/v8Cr9GEcAAQQQcCNA8LpxN12V4DWV4x1enRzv8EbFj4MggAACCLgRIHjduJuuSvCayhG8OjmCNyp+HAQBBBBAwI0AwevG3XRVgtdUjuDVyRG8UfHjIAgggAACbgQIXjfupqsSvKZyBK9OjuCNih8HQQABBBBwI0DwunE3XZXgNZUjeHVyBG9U/DgIAggggIAbAYLXjbvpqgSvqRzBq5MjeKPix0EQQAABBNwIELxu3E1XJXhN5QhenRzBGxU/DoIAAggg4EaA4HXjbroqwWsqR/Dq5AjeqPhxEAQQQAABNwIErxt301UJXlM5glcnR/BGxY+DIIAAAgi4ESB43bibrkrwmsoRvDo5gjcqfhwEAQQQQMCNAMHrxt10VYLXVI7g1ckRvFHx4yAIIIAAAm4ECF437qarErymcgSvTo7gjYofB0EAAQQQcCNA8LpxN12V4DWVI3h1cgRvVPw4CAIIIICAGwGC14276aoEr6kcwauTI3ij4sdBEEAAAQTcCBC8btxNVyV4TeUIXp0cwRsVPw6CAAIIIOBGgOB14266KsFrKkfw6uQI3qj4cRAEEEAAATcCBK8bd9NVCV5TOYJXJ0fwRsWPgyCAAAIIuBEgeN24m65K8JrKEbw6OYI3Kn4cBAEEEEDAjQDB68bddFWC11SO4NXJEbxR8eMgCCCAAAJuBAheN+6mqxK8pnIEr06O4I2KHwdBAAEEEHAjQPC6cTddleA1lSN4dXIEb1T8OAgCCCCAgBsBgteNu+mqBK+pHMGrkyN4o+LHQRBAAAEE3AgQvG7cTVcleE3lCF6dHMEbFT8OggACCCDgRoDgdeNuuirBaypH8OrkCN6o+HEQBBBAAAE3AgSvG3fTVQleUzmCVydH8EbFj4MggAACCLgRIHjduJuuSvCayhG8OjmCNyp+HAQBBBBAwI0AwevG3XRVgtdUjuDVyRG8UfHjIAgggAACbgQIXjfupqsSvKZyBK9OjuCNih8HQQABBBBwI0DwunE3XZXgNZUjeHVyBG9U/DgIAggggIAbAYLXjbvpqgSvqRzBq5MjeKPix0EQQAABBNwIELxu3E1XJXhN5QhenRzBGxU/DoIAAggg4EaA4HXjbroqwWsqR/Dq5AjeqPhxEAQQQAABNwIErxt301UJXlM5glcnR/BGxY+DIIAAAgi4ESB43bibrkrwmsoRvDo5gjcqfhwEAQQQQMCNAMHrxt10VYLXVI7g1ckRvFHx4yAIIIAAAm4ECF437qarErymcgSvTo7gjYofB0EAAQQQcCNA8LpxN12V4DWVI3h1cgRvVPw4CAIIIICAGwGC14276aoEr6kcwVstufwDBVJcUiIN0utXmhs25ki1jsOLEUAAAQQQiAcBgjcedqHq50DwVt3qRK9MKisrK9MdIvGnew2eLFu3f1npQkZl95abs3tLQWGR5OQukrVvbfZ/fm6bVpKXO1oaN0z3vyd4E3//uQIEEEAgjAIEb2LtOsGr2y+CV0S84O1xSQfp1rl9hWZ6/bqSkV5Pljz7oixfuV6W5k2RtNqpctOkedKyRZbcfftQgld37zGNAAIIIOBQgOB1iG+wNMFrgPa9EYL3u+DNvrqb9L3swuM0+w2fKl07tZPhg3r6P1uzfpOMm7ZQPlr3hCQlJfEOr+7+YxoBBBBAwJEAwesI3nBZgtcQ7rsxgve74K1bN01andlUmjZpJD0v7SAtmjXxidp1Hym5OTf40et9fbxlm/QfMU02rlwg3rvAPNKguwGZRgABBBBwI0DwunE3XZXgNZU7OkfwisiCJ16Q5JRk8Z5mXvvme7L9i53y+yXTpXnTTGnbeYgsnDVWOnY4zwfbum2H9MqeIq8ue0CymjQieHX3H9MIIIAAAo4EFj9Y45Qrh/4XfBzty8mWTU5KirMzSqzTIXiP2a8jR4ql68CJct2VXWTINd39d3hnThomXTqezzu8iXVvc7YIIIAAAqcQiPQOL3kVX7dPVqO0+DqhBDsbgvcEG3b1jdOl4wU/l5sHXyHeM7zeL7MNG9jDfyXP8CbYHc7pIoAAAgicUCBS8MIWXwI80qDbj9AH7+c7dvofOeZFbaMG6bJm3SbJmblInpo/WX5x7tmy+JlVsmLVBv9TGuqk1ZKROXP5lAbdPcc0AggggEAcCBC8cbAJ1TgFgrcaWCd4KcG7Y6dk3zZbdu7aW8GTM2qAXN+/q//9wYIimTDjYXn97Q/879u2bil5M8dIZuMM/3t+aU13AzKNAAIIIOBGgOB14266KsFrKnd0LvTB6yF4f3tjz758/49MeL+IViMl5TjV/fkHxXu+t/wPTpS/gODV3YBMI4AAAggggEBkgSUP1Yz8Il5xUgGCV3lzELxKQMYRQAABBBBAIKIAwRuR6JQvIHh1fjzSoPRjHAEEEEAAAQQiCxC8kY1O9QqCV+dH8Cr9GEcAAQQQQACByAIEb2QjgldndMppHmmIIS6HRgABBBBAAAFfgODV3Qi8w6vz4x1epR/jCCCAAAIIIBBZgOCNbMQ7vDoj3uGNoR+HRgABBBBAAIHIAgRvZCOCV2dE8MbQj0MjgAACCCCAQGQBgjeyEcGrMyJ4Y+jHoRFAAAEEEEAgsgDBG9mI4NUZEbwx9OPQCCCAAAIIIBBZgOCNbETw6owI3hj6cWgEEEAAAQQQiCxA8EY2Inh1RgRvDP04NAIIIIAAAghEFiB4IxsRvDojgjeGfhwaAQQQQAABBCILELyRjQhenRHBG0M/Do0AAggggAACkQUI3shGBK/OiOCNoR+HRgABBBBAAIHIAgRvZCOCV2dE8MbQj0MjgAACCCCAQGQBgjeyEcGrMyJ4Y+jHoRFAAAEEEEAgsgDBG9mI4NUZEbwx9OPQCCCAAAIIIBBZgOCNbETw6owI3hj6cWgEEEAAAQQQiCxA8EY2Inh1RgRvDP04NAIIIIAAAghEFiB4IxsRvDojgjeGfhwaAQQQQAABBCILELyRjQhenRHBG0M/Do0AAggggAACkQUI3shGBK/OiOCNoR+HRgABBBBAAIHIAgRvZCOCV2dE8MbQj0MjgAACCCCAQGQBgjeyEcGrMyJ4Y+jHoRFAAAEEEEAgsgDBG9mI4NUZEbwx9OPQCCCAAAIIIBBZgOCNbETw6owI3hj6cWgEEEAAAQQQiCxA8EY2Inh1RgRvDP04NAIIIIAAAghEFiB4IxsRvDojgjeGfhwaAQQQQAABBCILELyRjQhenRHBG0M/Do0AAggggAACkQUI3shGBK/OiOCNoR+HRgABBBBAAIHIAgRvZCOCV2dE8MbQj0MjgAACCCCAQGQBgjeyEcGrMyJ4Y+jHoRFAAAEEEEAgsgDBG9mI4NUZEbwx9OPQCCCAAAIIIBBZgOCNbETw6owI3hj6cWgEEEAAAQQQiCxA8EY2Inh1RgRvDP04NAIIIIAAAghEFiB4IxsRvDojgjeGfhwaAQQQQAABBCILELyRjQhenRHBG0M/Do0AAggggAACkQUI3shGBK/OiOCNoR+HRgABBBBAAIHIAgRvZCOCV2dE8MbQj0MjgAACCCCAQGQBgjeyEcGrM/Kn8w8USHFJiTRIr1/paMPGHInC0TkEAggggAACCCBwcgGCV3d3JJWVlZXpDhHs6YLCIsnJXSRr39rsX+i5bVpJXu5oadww3f+e4A32/nN1CCCAAAIIxIMAwavbBYI3gt+SZ1+U5SvXy9K8KZJWO1VumjRPWrbIkrtvH0rw6u49phFAAAEEEECgigIEbxWhTvIygjeCX7/hU6Vrp3YyfFBP/5Vr1m+ScdMWykfrnpCkpCTe4dXdf0wjgAACCCCAQBUECN4qIJ3ijpJKGwAACp9JREFUJQRvBL923UdKbs4NfvR6Xx9v2Sb9R0yTjSsXSHr9ugSv7v5jGgEEEEAAAQSqIEDwVgGJ4DVD8h5vbtt5iCycNVY6djjPP8jWbTukV/YUeXXZA5LVpBHBa0bLFAIIIIAAAghUQ4DgrQbWCV7KO7xVeId35qRh0qXj+Sd8h1fHzzQCCCCAAAIIIIBArAUI3gjC3jO83Tq3l2EDe/ivPPYZ3lhvEMdHAAEEEEAAAQQQ0AkQvBH8Fj+zSlas2uB/SkOdtFoyMmdupU9p0PEzjQACCCCAAAIIIBBrAYI3gvDBgiKZMONhef3tD/xXtm3dUvJmjpHMxhmx3htnxz/ZH9lwdkIsfEoB71nzktJSqZGSglScCxwpLpFvdu+Thg1+ILVSa8b52XJ63h8b+mbPfikrLZPMxg0kJSUZFAQQSFABgreKG7c//6AcOVJc8QcnqjiWUC+L9Ec2EupiQnSyK1/ZKPMWL5e1y+eF6KoT71K9/7fowcUrKk7c++SXqeOyJf0HdRPvYkJwxsv+tFZmzHuq4kqbnNZA5ueO9t/04Cu+BeY9uly8z9D/y6qF8oN6deL7ZDk7awIErzXq+F8o0h/ZiP8rCNcZfr5jpwyfcL988dUu8f7LmOCN7/1fvmq9NG+aKee1+ZH848uv5YZxc+SGAT0k++pu8X3iIT07739IZqTXk1+c29r/s/ITpi+U4uISeXxeTkhFEuOyX1j9htwx5zH/ZAnexNgzW2dJ8NqSToB1Iv2RjQS4hFCdYvn/3br2zc2y5NlVBG+C7f6d9z4uO77aRUAlyL55j7aVlpbJ3Gk3J8gZh+80//r+J3Lzrx+UGROH+I8iErzhuwdOdcUEL/dDhUCkP7IBVXwKrF77jtz38HMEb3xuzwnPynuWt+uACdLj4g4yfuRVCXTm4TvVP7/ylnj/o3LLZ/+QudNGyU9+1CJ8CAlwxdu/2CnemzYPzrhFmjRuIFcMmULwJsC+2TxFgtemdhyvVZU/shHHpx/qUyN4E2/7p97/hLz02jvy4tLZgf4F2MTbmePP2Hvu+t0Pt8jX3+yVu2+/Qdr/20+CcFmBuob93x6Uq26cJoOv6iYD+1wsf/+/HQRvoHY4OhdD8EbHMRBH8d7h5Y9sJN5WEryJtWcLn/yjLHjyj/LcI1PlnJ/wC1CJsnuLlq6Up3//irzxx7xEOeXQnGf55+Nf37+rJInInv354j2DffUVF0n/nh3lpz8+MzQWXOjJBQhe7o4KAf7IRmLeDARvYuyb9/znA48sk+dXrpffPjRJ2px9VmKcOGfpC7yy4b9l7NTfyAevPcZHAMbZPbF12w557c33Ks7K+yi5Z/7wqtx43eXS4+L/kFZnNYuzM+Z0XAgQvC7U43RN/shGnG7MSU7LewzF+63xl9dt8j+WbM2z90lSchL/ZRyn2+j95rj3G+SPzBkvPzwzq+IsvU/Y4DOU42/TvHfi/7P9OdK6VXPZvfdb/5eg0mql8kuG8bdVx50RjzQkwCY5OEWC1wF6vC4Zxj+yEa97UZXzKv8P9e+/9vIuF8jsySOqMs5rLAt0HTDR/wi5Y79eenqOnHlGE8tnw3KRBKbMXiJ/fPnNipf9W9sfy+wpI+SMrNMijfJzxwIEr+MNiNPlCd443RiXpxWGP7Lh0pe1EUAgMQQOHz4iX+/eJ/XqpPmfycsXAggkrgDBm7h7x5kjgAACCCCAAAIIVEGA4K0CEi9BAAEEEEAAAQQQSFwBgjdx944zRwABBBBAAAEEEKiCAMFbBSReggACCCCAAAIIIJC4AgRv4u4dZ44AAggggAACCCBQBQGCtwpIvAQBBBBAAAEEEEAgcQUI3sTdO84cAQQQQAABBBBAoAoCBG8VkHgJAggggAACCCCAQOIKELyJu3ecOQIIIIAAAggggEAVBAjeKiDxEgQQQAABBBBAAIHEFSB4E3fvOHMEEEAAAQQQQACBKggQvFVA4iUIIIAAAggggAACiStA8Cbu3nHmCCCAAAIIIIAAAlUQIHirgMRLEEAAAQQQQAABBBJXgOBN3L3jzBFAAAEEEEAAAQSqIEDwVgGJlyCAAAIIIIAAAggkrgDBm7h7x5kjgAACCCCAAAIIVEGA4K0CEi9BAAEEEEAAAQQQSFwBgjdx944zRwCBBBHIe/wPsmnzJzJz0jBp0SzTP+v/3foPuXveU3J1r85yeZcLEuRKOE0EEEAgMQUI3sTcN84aAQQSSOCbPfulz9A7pMlpDeXZBXfIkeIS6T9iqjRqkC6Pz8uRmjVSEuhqOFUEEEAg8QQI3sTbM84YAQQSUOC9/9ki1916jwzqe6nszz8gG//6kbzweK40bpiegFfDKSOAAAKJJUDwJtZ+cbYIIJDAAr9dvkbuXfA7/wqWLZoqbVu3TOCr4dQRQACBxBEgeBNnrzhTBBBIcIE33vlQRubM9a/ipafnyJlnNEnwK+L0EUAAgcQQIHgTY584SwQQSHCBL77aJX2G3indOreXdz/8X6mRkiLPPTJV6qTVSvAr4/QRQACB+BcgeON/jzhDBBBIcIHCosMy8OYZkpKS4v/S2meffyVXDrvL/3SG2ZNHJPjVcfoIIIBA/AsQvPG/R5whAggkuMDU+5+QFas2yOpn5kiLZkcfY3juT2v9jyWbPmGI9OvZMcGvkNNHAAEE4luA4I3v/eHsEEAAAQQQQAABBJQCBK8SkHEEEEAAAQQQQACB+BYgeON7fzg7BBBAAAEEEEAAAaUAwasEZBwBBBBAAAEEEEAgvgUI3vjeH84OAQQQQAABBBBAQClA8CoBGUcAAQQQQAABBBCIbwGCN773h7NDAAEEEEAAAQQQUAoQvEpAxhFAAAEEEEAAAQTiW4Dgje/94ewQQAABBBBAAAEElAIErxKQcQQQQAABBBBAAIH4FiB443t/ODsEEEAAAQQQQAABpQDBqwRkHAEEEEAAAQQQQCC+BQje+N4fzg4BBBBAAAEEEEBAKUDwKgEZRwABBBBAAAEEEIhvAYI3vveHs0MAAQQQQAABBBBQChC8SkDGEUAAAQQQQAABBOJbgOCN7/3h7BBAAAEEEEAAAQSUAgSvEpBxBBBAAAEEEEAAgfgWIHjje384OwQQQAABBBBAAAGlAMGrBGQcAQQQQAABBBBAIL4FCN743h/ODgEEEEAAAQQQQEApQPAqARlHAAEEEEAAAQQQiG8Bgje+94ezQwABBBBAAAEEEFAKELxKQMYRQAABBBBAAAEE4luA4I3v/eHsEEAAAQQQQAABBJQCBK8SkHEEEEAAAQQQQACB+BYgeON7fzg7BBBAAAEEEEAAAaUAwasEZBwBBBBAAAEEEEAgvgUI3vjeH84OAQQQQAABBBBAQClA8CoBGUcAAQQQQAABBBCIbwGCN773h7NDAAEEEEAAAQQQUAoQvEpAxhFAAAEEEEAAAQTiW4Dgje/94ewQQAABBBBAAAEElAIErxKQcQQQQAABBBBAAIH4FiB443t/ODsEEEAAAQQQQAABpQDBqwRkHAEEEEAAAQQQQCC+BQje+N4fzg4BBBBAAAEEEEBAKUDwKgEZRwABBBBAAAEEEIhvAYI3vveHs0MAAQQQQAABBBBQChC8SkDGEUAAAQQQQAABBOJb4P8DIbRYjeCmf38AAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#| caption: The same histogram produced by GTSAM's `DiscreteDistribution.sample` method.\n",
    "#| label: fig:discrete-samples-gtsam\n",
    "px.histogram(x=[category_prior.sample() for i in range(1000)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "_gm7mnKY68HZ",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## GTSAM 101\n",
    "\n",
    "> The GTSAM concepts used in this section, explained."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "IKt2DaIm1Brr",
   "metadata": {},
   "source": [
    "Above we created an instance of the `DiscreteDistribution` class. As with any GTSAM class, you can type\n",
    "```python\n",
    "help(gtsam.DiscreteDistribution)\n",
    "```\n",
    "{raw:tex}`\\noindent`\n",
    "to get documentation on its constructors and methods. In particular, we called the constructor\n",
    "```python\n",
    " __init__(self: gtsam.DiscreteDistribution, key: Tuple[int, int], spec: str) -> None\n",
    " ```\n",
    "{raw:tex}`\\noindent`\n",
    "which expects two arguments (besides `self`, which you can ignore):\n",
    "- `key`: Many GTSAM objects take a *key* to indicate which variable is involved. In the case of a DiscreteDistribution, the key is actually a tuple of ints:\n",
    "    - the first int is a 64-bit identifier for the variable;\n",
    "    - the second int is the *cardinality* of the variable.\n",
    "- `spec`: The `DiscreteDistribution` class specifies a PMF (remember: probability mass function) which is given as a string of numbers, separated by `/`.\n",
    "\n",
    "{raw:tex}`\\noindent`\n",
    "Let's look at an example below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "rVVfEO9DDv1T",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<p>  <i>P(42):</i></p>\n",
       "<div>\n",
       "<table class='DecisionTreeFactor'>\n",
       "  <thead>\n",
       "    <tr><th>42</th><th>value</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>0</th><td>0.4</td></tr>\n",
       "    <tr><th>1</th><td>0.1</td></tr>\n",
       "    <tr><th>2</th><td>0.5</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/markdown": [
       " *P(42):*\n",
       "\n",
       "|42|value|\n",
       "|:-:|:-:|\n",
       "|0|0.4|\n",
       "|1|0.1|\n",
       "|2|0.5|\n"
      ],
      "text/plain": [
       "Discrete Prior\n",
       " P( 42 ):\n",
       " Choice(42) \n",
       " 0 Leaf  0.4\n",
       " 1 Leaf  0.1\n",
       " 2 Leaf  0.5\n"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior = gtsam.DiscreteDistribution((42, 3), \"0.4/0.1/0.5\")\n",
    "prior\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sKXIMu0AsCvv",
   "metadata": {},
   "source": [
    "```{index} decision tree\n",
    "```\n",
    "As you can see, this is a PMF on the variable with id $42$, and it indeed has probabilities (that add up to one) for values `0..2`. Internally, GTSAM *actually* represents a PMF as a small **decision tree**, which you can reveal using `show`, as illustrated in [Figure 2.4](#fig:discrete-decision-tree)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "orGvLDaO2ZPv",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 9.0.0 (0)\n",
       " -->\n",
       "<!-- Title: G Pages: 1 -->\n",
       "<svg width=\"206pt\" height=\"109pt\"\n",
       " viewBox=\"0.00 0.00 206.00 109.46\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 105.46)\">\n",
       "<title>G</title>\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-105.46 202,-105.46 202,4 -4,4\"/>\n",
       "<!-- 0x1071efbf0 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>0x1071efbf0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-81.33\" rx=\"20.13\" ry=\"20.13\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-76.28\" font-family=\"Times,serif\" font-size=\"14.00\">42</text>\n",
       "</g>\n",
       "<!-- 0x12771bba0 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>0x12771bba0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"54,-25.2 0,-25.2 0,0 54,0 54,-25.2\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-7.55\" font-family=\"Times,serif\" font-size=\"14.00\"> 0.4</text>\n",
       "</g>\n",
       "<!-- 0x1071efbf0&#45;&gt;0x12771bba0 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>0x1071efbf0&#45;&gt;0x12771bba0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M84.43,-66.83C73.95,-57.12 59.67,-43.87 47.96,-33.03\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"50.45,-30.56 40.73,-26.33 45.69,-35.69 50.45,-30.56\"/>\n",
       "</g>\n",
       "<!-- 0x1277324d0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>0x1277324d0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"126,-25.2 72,-25.2 72,0 126,0 126,-25.2\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-7.55\" font-family=\"Times,serif\" font-size=\"14.00\"> 0.1</text>\n",
       "</g>\n",
       "<!-- 0x1071efbf0&#45;&gt;0x1277324d0 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>0x1071efbf0&#45;&gt;0x1277324d0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M99,-61.03C99,-53.45 99,-44.69 99,-36.77\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"102.5,-36.84 99,-26.84 95.5,-36.84 102.5,-36.84\"/>\n",
       "</g>\n",
       "<!-- 0x127715cf0 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>0x127715cf0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"198,-25.2 144,-25.2 144,0 198,0 198,-25.2\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-7.55\" font-family=\"Times,serif\" font-size=\"14.00\"> 0.5</text>\n",
       "</g>\n",
       "<!-- 0x1071efbf0&#45;&gt;0x127715cf0 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>0x1071efbf0&#45;&gt;0x127715cf0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M113.57,-66.83C124.05,-57.12 138.33,-43.87 150.04,-33.03\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"152.31,-35.69 157.27,-26.33 147.55,-30.56 152.31,-35.69\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<gtbook.display.show at 0x10556d340>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#| caption: The decision tree representation for a discrete probability distribution.\n",
    "#| label: fig:discrete-decision-tree\n",
    "show(prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "TDrC8HJbO4kt",
   "metadata": {},
   "source": [
    "Of course, it would be much nicer if we could print out these PMFs in a more readable format, that shows us a name for each variable as well as a pretty name for each. This is where the `Variables` class comes to the rescue. We actually defined a global variable at the top of this notebook, like so:\n",
    "\n",
    "```python\n",
    "VARIABLES = Variables()\n",
    "```\n",
    "\n",
    "{raw:tex}`\\noindent`\n",
    "which then allows us to give a name to a variable. It will also pick a unique ID for our variable. We can do this with the `discrete` method, which takes a name and a set of value names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "dVD_hUuwjbC8",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "key = (1, 3)\n"
     ]
    }
   ],
   "source": [
    "T = VARIABLES.discrete(\"TresCommas\", [\"one\", \"two\", \"three\"])\n",
    "print(f\"key = {T}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lX5GargYuy8L",
   "metadata": {},
   "source": [
    "As you can see, the id for the variable is 1 (`Category`, defined above, took id 0), and the cardinality was inferred to be three from the length of the list given as second argument.\n",
    "\n",
    "The `Variables` class can tell us about any discrete variable so defined, using the two methods `name` and `domain`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "AH4nq4S_CgoJ",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "key = (1, 3), name = TresCommas, domain = ['one', 'two', 'three']\n"
     ]
    }
   ],
   "source": [
    "print(f\"key = {T}, name = {VARIABLES.name(T)}, domain = {VARIABLES.domain(T)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fs7ewWz4PLP_",
   "metadata": {},
   "source": [
    "Conceptually, the `Variables` class keeps track of the names of variables and what values each variable can take on. We do this so that later when we print, it can show us a beautiful outputs. Of course, this only works if we actually use the key returned to use by `Variables.discrete`. We first define a prior using a string \"spec\" as shown in [Figure 2.5](#fig:unnormalized-spec). This also illustrates once again that the `spec` string does not need to contain normalized probabilities. The constructor will do the normalization for us! Note that the function `pretty(*)` above is just a shortcut for `gtbook.discrete.pretty(*, VARIABLES)`, and is also defined alongside `VARIABLES`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "Sy1PHI-doIJC",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<p>  <i>P(TresCommas):</i></p>\n",
       "<div>\n",
       "<table class='DecisionTreeFactor'>\n",
       "  <thead>\n",
       "    <tr><th>TresCommas</th><th>value</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>one</th><td>0.25</td></tr>\n",
       "    <tr><th>two</th><td>0.5</td></tr>\n",
       "    <tr><th>three</th><td>0.25</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<gtbook.display.pretty at 0x14234bc40>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#| caption: Creating a `DiscreteDistribution` using an unnormalized \"spec\".\n",
    "#| label: fig:unnormalized-spec\n",
    "prior_on_tres_commas = gtsam.DiscreteDistribution(T, \"2/4/2\")\n",
    "pretty(prior_on_tres_commas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "JHRmJg92Lb3p",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Finally, let us also look at the other `DiscreteDistribution` methods we called above:\n",
    "\n",
    "- `sample(self: gtsam.DiscreteConditional) -> int`: The `sample` method generates samples according to the PMF, returning the integer index of the sampled value, in $0\\dots cardinality-1$.\n",
    "- `pmf(self: gtsam.DiscreteDistribution) -> List[float]`: The `pmf` method will simply return all probability values, in order.\n",
    "- `__call__(self: gtsam.DiscreteDistribution, arg0: int) -> float`: The call operator: when given an integer value, will return just the one corresponding probability value.\n",
    "\n",
    "{raw:tex}`\\noindent`\n",
    "We have illustrated the use of all of these in the text, already."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3kXhNl6xYeo",
   "metadata": {},
   "source": [
    "Finally, we can also inspect the value of `VARIABLES`, as it has an HTML representation. The table shown in [Figure 2.6](#fig:nice-Variables-value) lists that we have defined two variables so far, `Category` and `TresCommas`, with respectively six and three possible values, enumerated in the `Domain` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "zMCqGQDEV7nn",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table class='Variables'>\n",
       "  <thead>\n",
       "    <tr><th>Variable</th><th>Domain</th></tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr><th>Category</th><td>cardboard, paper, can, scrap metal, bottle</td></tr>\n",
       "    <tr><th>TresCommas</th><td>one, two, three</td></tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "<gtbook.discrete.Variables at 0x10556d9a0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#| caption: A value of the `Variables` renders to a nice HTML output.\n",
    "#| label: fig:nice-Variables-value\n",
    "VARIABLES"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "S21_sorter_state.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "gtbook",
   "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.9.19"
  },
  "latex_metadata": {
   "affiliation": "Georgia Institute of Technology",
   "author": "Frank Dellaert and Seth Hutchinson",
   "title": "Introduction to Robotics"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}