{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(pymc_pytensor)=\n",
    "\n",
    "# PyMC and PyTensor\n",
    "\n",
    "**Authors:** [Ricardo Vieira](https://github.com/ricardoV94) and [Juan Orduz](https://juanitorduz.github.io/)\n",
    "\n",
    "In this notebook we want to give an introduction of how PyMC models translate to PyTensor graphs. The purpose is not to give a detailed description of all [`pytensor`](https://github.com/pymc-devs/pytensor)'s capabilities but rather focus on the main concepts to understand its connection with PyMC. For a more detailed description of the project please refer to the official documentation."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Prepare Notebook\n",
    "\n",
    "First import the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pytensor\n",
    "import pytensor.tensor as pt\n",
    "import scipy.stats\n",
    "\n",
    "import pymc as pm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to PyTensor\n",
    "\n",
    "We start by looking into `pytensor`. According to their documentation\n",
    "\n",
    "> PyTensor is a Python library that allows one to define, optimize, and efficiently evaluate mathematical expressions involving multi-dimensional arrays."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![pytensor logo](https://raw.githubusercontent.com/pymc-devs/pytensor/main/doc/images/PyTensor_RGB.svg)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple example\n",
    "\n",
    "To begin, we define some pytensor tensors and show how to perform some basic operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "x type: Scalar(float64, shape=())\n",
      "x name = x\n",
      "---\n",
      "y type: Vector(float64, shape=(?,))\n",
      "y name = y\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x = pt.scalar(name=\"x\")\n",
    "y = pt.vector(name=\"y\")\n",
    "\n",
    "print(\n",
    "    f\"\"\"\n",
    "x type: {x.type}\n",
    "x name = {x.name}\n",
    "---\n",
    "y type: {y.type}\n",
    "y name = {y.name}\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have defined the `x` and `y` tensors, we can create a new one by adding them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = x + y\n",
    "z.name = \"x + y\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the computation a bit more complex let us take the logarithm of the resulting tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = pt.log(z)\n",
    "w.name = \"log(x + y)\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the {func}`~pytensor.dprint` function to print the computational graph of any given tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log [id A] 'log(x + y)'\n",
      " └─ Add [id B] 'x + y'\n",
      "    ├─ ExpandDims{axis=0} [id C]\n",
      "    │  └─ x [id D]\n",
      "    └─ y [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(w)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this graph does not do any computation (yet!). It is simply defining the sequence of steps to be done. We can use {func}`~pytensor.function` to define a callable object so that we can push values trough the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = pytensor.function(inputs=[x, y], outputs=w)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the graph is compiled, we can push some concrete values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x=0, y=[1, np.e])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "Sometimes we just want to debug, we can use {meth}`~pytensor.graph.basic.Variable.eval` for that:\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.eval({x: 0, y: [1, np.e]})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can set intermediate values as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.eval({z: [1, np.e]})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTensor is clever!\n",
    "\n",
    "One of the most important features of `pytensor` is that it can automatically optimize the mathematical operations inside a graph. Let's consider a simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True_div [id A] 'a / b'\n",
      " ├─ a [id B]\n",
      " └─ b [id C]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = pt.scalar(name=\"a\")\n",
    "b = pt.scalar(name=\"b\")\n",
    "\n",
    "c = a / b\n",
    "c.name = \"a / b\"\n",
    "\n",
    "pytensor.dprint(c)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us multiply `b` times `c`. This should result in simply `a`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mul [id A] 'b * c'\n",
      " ├─ b [id B]\n",
      " └─ True_div [id C] 'a / b'\n",
      "    ├─ a [id D]\n",
      "    └─ b [id B]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = b * c\n",
    "d.name = \"b * c\"\n",
    "\n",
    "pytensor.dprint(d)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graph shows the full computation, but once we compile it the operation becomes the identity on `a` as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DeepCopyOp [id A] 0\n",
      " └─ a [id B]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = pytensor.function(inputs=[a, b], outputs=d)\n",
    "\n",
    "pytensor.dprint(g)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is in a PyTensor graph?\n",
    "\n",
    "The following diagram shows the basic structure of an `pytensor` graph.\n",
    "\n",
    "![pytensor graph](https://raw.githubusercontent.com/pymc-devs/pytensor/main/doc/tutorial/apply.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can can make these concepts more tangible by explicitly indicating them in the first example from the section above. Let us compute the graph components for the tensor `z`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "z type: Vector(float64, shape=(?,))\n",
      "z name = x + y\n",
      "z owner = Add(ExpandDims{axis=0}.0, y)\n",
      "z owner inputs = [ExpandDims{axis=0}.0, y]\n",
      "z owner op = Add\n",
      "z owner output = [x + y]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"\"\"\n",
    "z type: {z.type}\n",
    "z name = {z.name}\n",
    "z owner = {z.owner}\n",
    "z owner inputs = {z.owner.inputs}\n",
    "z owner op = {z.owner.op}\n",
    "z owner output = {z.owner.outputs}\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code snippet helps us understand these concepts by going through the computational graph of `w`. The actual code is not as important here, the focus is on the outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---\n",
      "Checking variable log(x + y) of type Vector(float64, shape=(?,))\n",
      " > Op is Log\n",
      " > Input 0 is x + y\n",
      "---\n",
      "Checking variable x + y of type Vector(float64, shape=(?,))\n",
      " > Op is Add\n",
      " > Input 0 is ExpandDims{axis=0}.0\n",
      " > Input 1 is y\n",
      "---\n",
      "Checking variable ExpandDims{axis=0}.0 of type Vector(float64, shape=(1,))\n",
      " > Op is ExpandDims{axis=0}\n",
      " > Input 0 is x\n",
      "---\n",
      "Checking variable y of type Vector(float64, shape=(?,))\n",
      " > y is a root variable\n",
      "---\n",
      "Checking variable x of type Scalar(float64, shape=())\n",
      " > x is a root variable\n"
     ]
    }
   ],
   "source": [
    "# start from the top\n",
    "stack = [w]\n",
    "\n",
    "while stack:\n",
    "    print(\"---\")\n",
    "    var = stack.pop(0)\n",
    "    print(f\"Checking variable {var} of type {var.type}\")\n",
    "    # check variable is not a root variable\n",
    "    if var.owner is not None:\n",
    "        print(f\" > Op is {var.owner.op}\")\n",
    "        # loop over the inputs\n",
    "        for i, input in enumerate(var.owner.inputs):\n",
    "            print(f\" > Input {i} is {input}\")\n",
    "            stack.append(input)\n",
    "    else:\n",
    "        print(f\" > {var} is a root variable\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this is very similar to the output of {func}`~pytensor.dprint` function introduced above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log [id A] 'log(x + y)'\n",
      " └─ Add [id B] 'x + y'\n",
      "    ├─ ExpandDims{axis=0} [id C]\n",
      "    │  └─ x [id D]\n",
      "    └─ y [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(w)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph manipulation 101\n",
    "\n",
    "Another interesting feature of PyTensor is the ability to manipulate the computational graph, something that is not possible with TensorFlow or PyTorch. Here we continue with the example above in order to illustrate the main idea around this technique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x, y]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get input tensors\n",
    "list(pytensor.graph.graph_inputs(graphs=[w]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a simple example, let's add an {func}`~pytensor.tensor.exp` before the {func}`~pytensor.tensor.log` (to get the identity function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "parent_of_w = w.owner.inputs[0]  # get z tensor\n",
    "new_parent_of_w = pt.exp(parent_of_w)  # modify the parent of w\n",
    "new_parent_of_w.name = \"exp(x + y)\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the graph of `w` has actually not changed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log [id A] 'log(x + y)'\n",
      " └─ Add [id B] 'x + y'\n",
      "    ├─ ExpandDims{axis=0} [id C]\n",
      "    │  └─ x [id D]\n",
      "    └─ y [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(w)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To modify the graph we need to use the {func}`~pytensor.clone_replace` function, which *returns a copy of the initial subgraph with the corresponding substitutions.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log [id A] 'log(exp(x + y))'\n",
      " └─ Exp [id B] 'exp(x + y)'\n",
      "    └─ Add [id C] 'x + y'\n",
      "       ├─ ExpandDims{axis=0} [id D]\n",
      "       │  └─ x [id E]\n",
      "       └─ y [id F]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_w = pytensor.clone_replace(output=[w], replace={parent_of_w: new_parent_of_w})[0]\n",
    "new_w.name = \"log(exp(x + y))\"\n",
    "pytensor.dprint(new_w)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can test the modified graph by passing some input to the new graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 2.71828183])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_w.eval({x: 0, y: [1, np.e]})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the new graph is just the identity function."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "Again, note that `pytensor` is clever enough to omit the `exp` and `log` once we compile the function.\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Add [id A] 'x + y' 1\n",
      " ├─ ExpandDims{axis=0} [id B] 0\n",
      " │  └─ x [id C]\n",
      " └─ y [id D]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = pytensor.function(inputs=[x, y], outputs=new_w)\n",
    "\n",
    "pytensor.dprint(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 2.71828183])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(x=0, y=[1, np.e])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyTensor RandomVariables\n",
    "\n",
    "Now that we have seen pytensor's basics we want to move in the direction of random variables.\n",
    "\n",
    "How do we generate random numbers in [`numpy`](https://github.com/numpy/numpy)? To illustrate it we can sample from a normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAIQCAYAAACBuKavAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCsklEQVR4nO3de3yP9f/H8ednZgeHbTazWexAckqSwyxyXI3kFEqJkcg5VmF9Q3zLVIrIIX1z6BudiJIiOX7LHKJSisgcog2xLZON7fr90W2fn08bZsbns7fH/Xa7bjef93Vd78/ruj7Xtqfr876uy2ZZliUAAADAQG7OLgAAAAC4Vgi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLvAdWaz2fTcc89d1/dMSUlR165dFRAQIJvNpqlTp17X98fVOXDggGw2m+bPn1+o9efPny+bzaYDBw7Y21q0aKEWLVoUSX2X889j/rnnnpPNZtOJEyeuy/uHh4erd+/e1+W9rpWrPQaAGxlhF8XSDz/8oK5duyosLExeXl666aabdPfdd2v69OnOLs0ljRgxQqtWrVJ8fLz++9//qk2bNs4uCcXQpk2b9Nxzzyk1NdXZpeThyrUBcC53ZxcAXKlNmzapZcuWCg0NVb9+/RQcHKzDhw9r8+bNeu211zR06FBnl+hy1q5dq44dO+qpp55ydilwEV988cUVr7Np0yaNHz9evXv3lp+fX4HX++uvv+Tufm3/3Fyqtj179sjNrXif2wkLC9Nff/2lkiVLOrsUoNgh7KLYeeGFF+Tr66tt27bl+aN27Ngx5xTl4o4dO1agcJKRkaHSpUtf+4KKubNnz8rDw6NYBygPD49r2n9OTo6ysrLk5eUlLy+va/pel+Pp6enU9y8KNpvN6fsRKK6K729q3LB+/fVX1a5dO9/wVqFCBYfX8+bNU6tWrVShQgV5enqqVq1amjVrVp71wsPDdd9992n9+vVq0KCBvL29VadOHa1fv16S9NFHH6lOnTry8vJS/fr19e233zqs37t3b5UpU0b79+9XTEyMSpcurZCQEE2YMEGWZV12m44cOaJHH31UQUFB8vT0VO3atTV37tw8y02fPl21a9dWqVKlVK5cOTVo0ECLFi26aL+5YzUty9KMGTNks9lks9kc5m3YsEGDBg1ShQoVVKlSJfu6M2fOVO3ateXp6amQkBANHjw4z1fELVq00K233qqdO3eqefPmKlWqlG6++WYtXrxYkrRhwwZFRkbK29tb1atX15dffnnZfZGVlaWxY8eqfv368vX1VenSpXXXXXdp3bp1l11X+v/P8quvvlKjRo3k5eWlKlWq6O23386z7P79+9WtWzf5+/urVKlSaty4sVasWOGwzPr162Wz2fTee+/p2Wef1U033aRSpUopPT3d/rkfOnRI9913n8qUKaObbrpJM2bMkPT3cJtWrVqpdOnSCgsLy/NZnTx5Uk899ZTq1KmjMmXKyMfHR23bttX3339foG3Nz65du9SqVSt5e3urUqVKev7555WTk5NnufzG7F7q+Hruuef09NNPS5IiIiLsx1LuOGCbzaYhQ4Zo4cKF9uNm5cqV9nn5jVM/ceKEHnjgAfn4+CggIEBPPPGEzp49a59/qXGqF/Z5udryG7N7JZ/9Bx98oBdeeEGVKlWSl5eXWrdurX379uWp6Z969+6t8PDwPO25Y5YvtHr1ajVt2lR+fn4qU6aMqlevrmeeeeaS+yL3+Dty5Ig6deqkMmXKKDAwUE899ZSys7Md+v/jjz/Us2dP+fj4yM/PT7Gxsfr+++8LNA4493fF119/rbi4OAUGBqp06dLq3Lmzjh8/7rDsxT7rf34GuX1+9dVXGjZsmAIDA+Xn56fHH39cWVlZSk1NVa9evVSuXDmVK1dOI0eOdPhdmrs/Jk+erClTpigsLEze3t5q3ry5fvzxR/ty8+bNk81my/M7W5ImTpyoEiVK6MiRI5fcfhR/nNlFsRMWFqbExET9+OOPuvXWWy+57KxZs1S7dm116NBB7u7uWr58uQYNGqScnBwNHjzYYdl9+/bp4Ycf1uOPP65HHnlEkydPVvv27TV79mw988wzGjRokCQpISFBDzzwQJ6vRrOzs9WmTRs1btxYL730klauXKlx48bp/PnzmjBhwkVrTElJUePGje1hITAwUJ9//rn69u2r9PR0DR8+XJL05ptvatiwYeratas9FOzcuVNbtmzRww8/nG/fzZo103//+1/17NlTd999t3r16pVnmUGDBikwMFBjx45VRkaGpL//GI8fP17R0dEaOHCg9uzZo1mzZmnbtm36+uuvHb5KPXXqlO677z51795d3bp106xZs9S9e3ctXLhQw4cP14ABA/Twww/r5ZdfVteuXXX48GGVLVv2ovsjPT1d//nPf/TQQw+pX79++vPPP/XWW28pJiZGW7du1e23337RdXPt27dPXbt2Vd++fRUbG6u5c+eqd+/eql+/vmrXrm3f73feeafOnDmjYcOGKSAgQAsWLFCHDh20ePFide7c2aHPf//73/Lw8NBTTz2lzMxM+5nR7OxstW3bVs2aNdNLL72khQsXasiQISpdurT+9a9/qUePHrr//vs1e/Zs9erVS1FRUYqIiJD0d+BatmyZunXrpoiICKWkpOiNN95Q8+bN9dNPPykkJOSy23qh5ORktWzZUufPn9fo0aNVunRpzZkzR97e3pdd93LH1/33369ffvlF7777rqZMmaLy5ctLkgIDA+19rF27Vh988IGGDBmi8uXL5xv0LvTAAw8oPDxcCQkJ2rx5s6ZNm6ZTp07l+x+TSylIbRe60s9+0qRJcnNz01NPPaW0tDS99NJL6tGjh7Zs2XJFdV7Mrl27dN999+m2227ThAkT5OnpqX379unrr7++7LrZ2dmKiYlRZGSkJk+erC+//FKvvPKKqlatqoEDB0r6+yx7+/bttXXrVg0cOFA1atTQxx9/rNjY2Cuqc+jQoSpXrpzGjRunAwcOaOrUqRoyZIjef//9Qm13bp/BwcEaP368Nm/erDlz5sjPz0+bNm1SaGioJk6cqM8++0wvv/yybr311jy/w95++239+eefGjx4sM6ePavXXntNrVq10g8//KCgoCB17dpVgwcP1sKFC1WvXj2HdRcuXKgWLVropptuKnT9KCYsoJj54osvrBIlSlglSpSwoqKirJEjR1qrVq2ysrKy8ix75syZPG0xMTFWlSpVHNrCwsIsSdamTZvsbatWrbIkWd7e3tbBgwft7W+88YYlyVq3bp29LTY21pJkDR061N6Wk5NjtWvXzvLw8LCOHz9ub5dkjRs3zv66b9++VsWKFa0TJ0441NS9e3fL19fXvg0dO3a0ateufZm9kz9J1uDBgx3a5s2bZ0mymjZtap0/f97efuzYMcvDw8O65557rOzsbHv766+/bkmy5s6da29r3ry5JclatGiRvW337t2WJMvNzc3avHmzvT13f86bN++StZ4/f97KzMx0aDt16pQVFBRkPfroo5fd1tzPcuPGjQ7b5OnpaT355JP2tuHDh1uSrP/973/2tj///NOKiIiwwsPD7du+bt06S5JVpUqVPMdT7uc+ceJEh1q9vb0tm81mvffee3n2y4Wf/dmzZx32sWVZVlJSkuXp6WlNmDDBoa0g+y53m7Zs2eKw7b6+vpYkKykpyd7evHlzq3nz5vbXBTm+Xn755Tz95Mr9zHft2pXvvAu3e9y4cZYkq0OHDg7LDRo0yJJkff/995ZlXXq7/9nnpWoLCwuzYmNj7a+v9LOvWbOmwzH52muvWZKsH374Ic97XSg2NtYKCwvL0567/bmmTJliSXL4PfFP+e2L3OPvwmPFsiyrXr16Vv369e2vlyxZYkmypk6dam/Lzs62WrVqVaDjKvd3RXR0tJWTk2NvHzFihFWiRAkrNTXV3vbPzyXXPz+D3D5jYmIc+oyKirJsNps1YMAAe9v58+etSpUqORyvufvD29vb+u233+ztW7ZssSRZI0aMsLc99NBDVkhIiMPP2o4dOwq07TADwxhQ7Nx9991KTExUhw4d9P333+ull15STEyMbrrpJn3yyScOy154RistLU0nTpxQ8+bNtX//fqWlpTksW6tWLUVFRdlfR0ZGSpJatWql0NDQPO379+/PU9uQIUPs/849U5uVlXXRr+8ty9KSJUvUvn17WZalEydO2KeYmBilpaVpx44dkiQ/Pz/99ttv2rZtW4H2U0H169dPJUqUsL/+8ssvlZWVpeHDhzucue7Xr598fHzyfNVbpkwZde/e3f66evXq8vPzU82aNe37Srr0frtQiRIl7GdNc3JydPLkSZ0/f14NGjSw74vLqVWrlu666y7768DAQFWvXt3hvT/77DM1atRITZs2ddiW/v3768CBA/rpp58c+oyNjb3oGdLHHnvM/m8/Pz9Vr15dpUuX1gMPPGBvz90vF9bg6elp38fZ2dn6448/7F9hF3RbL/TZZ5+pcePGatSokcO29+jR47LrFsXx1bx5c9WqVavAy//z25Xci0s/++yzQtdQEFf62ffp08dhjHPusXW5Y7mgcodkffzxx/kOObmcAQMGOLy+6667HGpbuXKlSpYsqX79+tnb3Nzc8uz/y+nfv7/D8Iu77rpL2dnZOnjw4BXXnKtv374OfUZGRsqyLPXt29feVqJECTVo0CDf/d2pUyeHM7ONGjVSZGSkwzHUq1cvHT161GEo1MKFC+Xt7a0uXboUunYUH4RdFEsNGzbURx99pFOnTmnr1q2Kj4/Xn3/+qa5duzr8ofr6668VHR2t0qVLy8/PT4GBgfZxcP8MuxcGWkny9fWVJFWuXDnf9lOnTjm0u7m5qUqVKg5tt9xyiyQ53N/0QsePH1dqaqrmzJmjwMBAh6lPnz6S/v+iu1GjRqlMmTJq1KiRqlWrpsGDBxfoa87Lyf1KPVfuH67q1as7tHt4eKhKlSp5/rBVqlQpz/hDX1/fAu+3/CxYsEC33XabvLy8FBAQoMDAQK1YsSLPZ3Yx//wsJalcuXIO733w4ME82yhJNWvWtM+/0D/3Uy4vL688X5f7+vpedL9cWENOTo6mTJmiatWqydPTU+XLl1dgYKB27txZ4G290MGDB1WtWrU87flt5z8VxfF1sX10Mf+stWrVqnJzc7voz0tRudLP/p/HU7ly5SQV7FguiAcffFBNmjTRY489pqCgIHXv3l0ffPBBgYJvfsdffsd6xYoVVapUKYflbr755iuq81rshyv5vZvf++R3vN9yyy0Ox9Ddd9+tihUrauHChZL+/rl799131bFjx0sOqYI5CLso1jw8PNSwYUNNnDhRs2bN0rlz5/Thhx9K+vtCttatW+vEiRN69dVXtWLFCq1evVojRoyQpDx/SC48u1mQdqsAF55dTm4NjzzyiFavXp3v1KRJE0l//yHes2eP3nvvPTVt2lRLlixR06ZNNW7cuKuqoSDjOS+lqPfbO++8o969e6tq1ap66623tHLlSq1evVqtWrUq8Fmva/GZXWw/Xc32T5w4UXFxcWrWrJneeecdrVq1SqtXr1bt2rULdYbvahTF8XW1x9I//3Pwz9e5/nnx1bVW2OOpoPV7e3tr48aN+vLLL9WzZ0/t3LlTDz74oO6+++7LbuvFarsWrubn6mLbcSU/P4X9+S1RooQefvhhLVmyRGfPntW6det09OhRPfLII4XqD8UPF6jBGA0aNJAk/f7775Kk5cuXKzMzU5988onD2YOCXtV/pXJycrR//3772VxJ+uWXXyTpohfqBAYGqmzZssrOzlZ0dPRl36N06dJ68MEH9eCDDyorK0v333+/XnjhBcXHxxfZbYnCwsIk/X1v0gvPVGdlZSkpKalAdV6NxYsXq0qVKvroo48cwsLVhvp/CgsL0549e/K079692z7/Wlu8eLFatmypt956y6E9NTXVfpHVlQgLC9PevXvztOe3nfm53PF1sfBWWHv37nU4G7xv3z7l5OTYf15yzxz+8y4g+X1tfiW1Xa/Pvly5cvk+5CK/+t3c3NS6dWu1bt1ar776qiZOnKh//etfWrdu3VX/zIWFhWndunU6c+aMw9ndgtxR4krlt81ZWVn238tFLb/j/ZdffsnzO7dXr1565ZVXtHz5cn3++ecKDAxUTEzMNakJroczuyh21q1bl+//8HPHaOV+PZl7ZuDCZdPS0jRv3rxrVtvrr79u/7dlWXr99ddVsmRJtW7dOt/lS5QooS5dumjJkiUOt8vJdeFtff744w+HeR4eHqpVq5Ysy9K5c+eKaAuk6OhoeXh4aNq0aQ777q233lJaWpratWtXZO+Vn/w+ty1btigxMbFI3+fee+/V1q1bHfrNyMjQnDlzFB4efkVjTwurRIkSeY7lDz/8sNC3Qrr33nu1efNmbd261d52/Phx+9e3l1KQ4yv3HsxF9ZSy3Fu05cp9AmLbtm0lST4+Pipfvrw2btzosNzMmTPz9HUltV2vz75q1apKS0vTzp077W2///67li5d6rDcyZMn86ybe9eRzMzMq64jJiZG586d05tvvmlvy8nJybP/i0LVqlXzfF5z5sy5Zmfjly1b5vDzsnXrVm3ZssV+DOW67bbbdNttt+k///mPlixZou7du1/zB53AdfBJo9gZOnSozpw5o86dO6tGjRrKysrSpk2b9P777ys8PNw+1vWee+6Rh4eH2rdvr8cff1ynT5/Wm2++qQoVKlyTswxeXl5auXKlYmNjFRkZqc8//1wrVqzQM888c9FbIEl/39Zo3bp1ioyMVL9+/VSrVi2dPHlSO3bs0Jdffmn/Q3jPPfcoODhYTZo0UVBQkH7++We9/vrrateuXZGOOwsMDFR8fLzGjx+vNm3aqEOHDtqzZ49mzpyphg0bXvOv/u677z599NFH6ty5s9q1a6ekpCTNnj1btWrV0unTp4vsfUaPHq13331Xbdu21bBhw+Tv768FCxYoKSlJS5YsuS4PjLjvvvs0YcIE9enTR3feead++OEHLVy4MM/Y74IaOXKk/XHQTzzxhP3WY2FhYQ6BKz8FOb7q168vSfrXv/6l7t27q2TJkmrfvn2hH0SSlJSkDh06qE2bNkpMTNQ777yjhx9+WHXr1rUv89hjj2nSpEl67LHH1KBBA23cuNH+jcmFrqS26/XZd+/eXaNGjVLnzp01bNgwnTlzRrNmzdItt9zicAHihAkTtHHjRrVr105hYWE6duyYZs6cqUqVKjlcRFdYnTp1UqNGjfTkk09q3759qlGjhj755BP775aiPGP/2GOPacCAAerSpYvuvvtuff/991q1alWhvqkoiJtvvllNmzbVwIEDlZmZqalTpyogIEAjR47Ms2yvXr3sT5FkCMONhbCLYmfy5Mn68MMP9dlnn2nOnDnKyspSaGioBg0apGeffdZ+ZXP16tW1ePFiPfvss3rqqacUHBysgQMHKjAwUI8++miR11WiRAmtXLlSAwcO1NNPP62yZctq3LhxGjt27CXXCwoK0tatWzVhwgR99NFHmjlzpgICAlS7dm29+OKL9uUef/xxLVy4UK+++qpOnz6tSpUqadiwYXr22WeLfFuee+45BQYG6vXXX9eIESPk7++v/v37a+LEidf8caW9e/dWcnKy3njjDa1atUq1atXSO++8ow8//ND+kI+iEBQUpE2bNmnUqFGaPn26zp49q9tuu03Lly+/5mevcz3zzDPKyMjQokWL9P777+uOO+7QihUrNHr06EL1V7FiRa1bt05Dhw7VpEmTFBAQoAEDBigkJMTh6vb8FOT4atiwof79739r9uzZWrlypXJycpSUlFTosPv+++9r7NixGj16tNzd3TVkyBC9/PLLDsuMHTtWx48f1+LFi/XBBx+obdu2+vzzz/M8QOZKarten31AQICWLl2quLg4jRw5UhEREUpISNDevXsdwm6HDh104MABzZ07VydOnFD58uXVvHlzjR8/3n7B1tUoUaKEVqxYoSeeeEILFiyQm5ubOnfurHHjxqlJkyZF+mS2fv36KSkpyT7e/q677tLq1asv+u3W1erVq5fc3Nw0depUHTt2TI0aNdLrr7+uihUr5lm2R48eGjVqlKpWrepwxxKYz2YVxVU2wA2ud+/eWrx4cZGeeQSAa2nZsmXq3LmzvvrqK/uFsMXFgQMHFBERoZdfftl+tvZyTpw4oYoVK2rs2LEaM2bMNa4QroQxuwAAGO6vv/5yeJ2dna3p06fLx8dHd9xxh5Oqur7mz5+v7Oxs9ezZ09ml4DpjGAMAAIYbOnSo/vrrL0VFRSkzM1MfffSRNm3apIkTJ171LeNc3dq1a/XTTz/phRdeUKdOnS77GGuYh7ALAIDhWrVqpVdeeUWffvqpzp49q5tvvlnTp093eOqjqSZMmKBNmzapSZMm9jt+4MbCmF0AAAAYizG7AAAAMBZhFwAAAMZizK7+fpLM0aNHVbZs2SJ/HCYAAACunmVZ+vPPPxUSEnJFD38h7Eo6evSoKleu7OwyAAAAcBmHDx9WpUqVCrw8YVeyPwrz8OHD8vHxcXI1AAAA+Kf09HRVrlzZntsKirCr/38uuI+PD2EXAADAhV3pkFMuUAMAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLGcGnY3btyo9u3bKyQkRDabTcuWLcuzzM8//6wOHTrI19dXpUuXVsOGDXXo0CH7/LNnz2rw4MEKCAhQmTJl1KVLF6WkpFzHrQAAAICrcmrYzcjIUN26dTVjxox85//6669q2rSpatSoofXr12vnzp0aM2aMvLy87MuMGDFCy5cv14cffqgNGzbo6NGjuv/++6/XJgAAAMCF2SzLspxdhPT3c46XLl2qTp062du6d++ukiVL6r///W++66SlpSkwMFCLFi1S165dJUm7d+9WzZo1lZiYqMaNGxfovdPT0+Xr66u0tDT5+Phc9bYAAACgaBU2r7nsmN2cnBytWLFCt9xyi2JiYlShQgVFRkY6DHXYvn27zp07p+joaHtbjRo1FBoaqsTERCdUDQAAAFfismH32LFjOn36tCZNmqQ2bdroiy++UOfOnXX//fdrw4YNkqTk5GR5eHjIz8/PYd2goCAlJydftO/MzEylp6c7TAAAADCPu7MLuJicnBxJUseOHTVixAhJ0u23365NmzZp9uzZat68eaH7TkhI0Pjx44ukTgAAALgulz2zW758ebm7u6tWrVoO7TVr1rTfjSE4OFhZWVlKTU11WCYlJUXBwcEX7Ts+Pl5paWn26fDhw0VePwAAAJzPZcOuh4eHGjZsqD179ji0//LLLwoLC5Mk1a9fXyVLltSaNWvs8/fs2aNDhw4pKirqon17enrKx8fHYQIAAIB5nDqM4fTp09q3b5/9dVJSkr777jv5+/srNDRUTz/9tB588EE1a9ZMLVu21MqVK7V8+XKtX79ekuTr66u+ffsqLi5O/v7+8vHx0dChQxUVFVXgOzEAAADAXE699dj69evVsmXLPO2xsbGaP3++JGnu3LlKSEjQb7/9purVq2v8+PHq2LGjfdmzZ8/qySef1LvvvqvMzEzFxMRo5syZlxzG8E/cegwAAMC1FTavucx9dp2JsAsAxU/46BXOLiFfBya1c3YJgJGMu88uAAAAcLUIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIzl7uwCAAD/L3z0CmeXkK8Dk9o5uwQAKBTO7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCynht2NGzeqffv2CgkJkc1m07Jlyy667IABA2Sz2TR16lSH9pMnT6pHjx7y8fGRn5+f+vbtq9OnT1/bwgEAAFAsODXsZmRkqG7dupoxY8Yll1u6dKk2b96skJCQPPN69OihXbt2afXq1fr000+1ceNG9e/f/1qVDAAAgGLE3Zlv3rZtW7Vt2/aSyxw5ckRDhw7VqlWr1K5dO4d5P//8s1auXKlt27apQYMGkqTp06fr3nvv1eTJk/MNxwCAKxc+eoWzSwCAQnHpMbs5OTnq2bOnnn76adWuXTvP/MTERPn5+dmDriRFR0fLzc1NW7ZsuWi/mZmZSk9Pd5gAAABgHpcOuy+++KLc3d01bNiwfOcnJyerQoUKDm3u7u7y9/dXcnLyRftNSEiQr6+vfapcuXKR1g0AAADX4LJhd/v27Xrttdc0f/582Wy2Iu07Pj5eaWlp9unw4cNF2j8AAABcg8uG3f/97386duyYQkND5e7uLnd3dx08eFBPPvmkwsPDJUnBwcE6duyYw3rnz5/XyZMnFRwcfNG+PT095ePj4zABAADAPE69QO1SevbsqejoaIe2mJgY9ezZU3369JEkRUVFKTU1Vdu3b1f9+vUlSWvXrlVOTo4iIyOve80AAABwLU4Nu6dPn9a+ffvsr5OSkvTdd9/J399foaGhCggIcFi+ZMmSCg4OVvXq1SVJNWvWVJs2bdSvXz/Nnj1b586d05AhQ9S9e3fuxAAAAADnDmP45ptvVK9ePdWrV0+SFBcXp3r16mns2LEF7mPhwoWqUaOGWrdurXvvvVdNmzbVnDlzrlXJAAAAKEacema3RYsWsiyrwMsfOHAgT5u/v78WLVpUhFUBAADAFC57gRoAAABwtQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLHcnV0AADhD+OgVzi4BAHAdcGYXAAAAxiLsAgAAwFiEXQAAABjLqWF348aNat++vUJCQmSz2bRs2TL7vHPnzmnUqFGqU6eOSpcurZCQEPXq1UtHjx516OPkyZPq0aOHfHx85Ofnp759++r06dPXeUsAAADgipwadjMyMlS3bl3NmDEjz7wzZ85ox44dGjNmjHbs2KGPPvpIe/bsUYcOHRyW69Gjh3bt2qXVq1fr008/1caNG9W/f//rtQkAAABwYTbLsixnFyFJNptNS5cuVadOnS66zLZt29SoUSMdPHhQoaGh+vnnn1WrVi1t27ZNDRo0kCStXLlS9957r3777TeFhIQU6L3T09Pl6+urtLQ0+fj4FMXmAHBx3I0B18qBSe2cXQJgpMLmtWI1ZjctLU02m01+fn6SpMTERPn5+dmDriRFR0fLzc1NW7ZsuWg/mZmZSk9Pd5gAAABgnmITds+ePatRo0bpoYcesqf55ORkVahQwWE5d3d3+fv7Kzk5+aJ9JSQkyNfX1z5Vrlz5mtYOAAAA5ygWYffcuXN64IEHZFmWZs2addX9xcfHKy0tzT4dPny4CKoEAACAq3H5J6jlBt2DBw9q7dq1DmM0goODdezYMYflz58/r5MnTyo4OPiifXp6esrT0/Oa1QwAAADX4NJndnOD7t69e/Xll18qICDAYX5UVJRSU1O1fft2e9vatWuVk5OjyMjI610uAAAAXIxTz+yePn1a+/bts79OSkrSd999J39/f1WsWFFdu3bVjh079Omnnyo7O9s+Dtff318eHh6qWbOm2rRpo379+mn27Nk6d+6chgwZou7duxf4TgwAAAAwl1PD7jfffKOWLVvaX8fFxUmSYmNj9dxzz+mTTz6RJN1+++0O661bt04tWrSQJC1cuFBDhgxR69at5ebmpi5dumjatGnXpX4AAAC4NqeG3RYtWuhSt/ktyC2A/f39tWjRoqIsCwAAAIZw6TG7AAAAwNUg7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwlruzCwAAwCTho1c4u4Q8Dkxq5+wSAKfhzC4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwllPD7saNG9W+fXuFhITIZrNp2bJlDvMty9LYsWNVsWJFeXt7Kzo6Wnv37nVY5uTJk+rRo4d8fHzk5+envn376vTp09dxKwAAAOCqnBp2MzIyVLduXc2YMSPf+S+99JKmTZum2bNna8uWLSpdurRiYmJ09uxZ+zI9evTQrl27tHr1an366afauHGj+vfvf702AQAAAC7M3Zlv3rZtW7Vt2zbfeZZlaerUqXr22WfVsWNHSdLbb7+toKAgLVu2TN27d9fPP/+slStXatu2bWrQoIEkafr06br33ns1efJkhYSEXLdtAQAAgOtx2TG7SUlJSk5OVnR0tL3N19dXkZGRSkxMlCQlJibKz8/PHnQlKTo6Wm5ubtqyZct1rxkAAACuxalndi8lOTlZkhQUFOTQHhQUZJ+XnJysChUqOMx3d3eXv7+/fZn8ZGZmKjMz0/46PT29qMoGAACAC3HZM7vXUkJCgnx9fe1T5cqVnV0SAAAArgGXDbvBwcGSpJSUFIf2lJQU+7zg4GAdO3bMYf758+d18uRJ+zL5iY+PV1pamn06fPhwEVcPAAAAV+CyYTciIkLBwcFas2aNvS09PV1btmxRVFSUJCkqKkqpqanavn27fZm1a9cqJydHkZGRF+3b09NTPj4+DhMAAADM49Qxu6dPn9a+ffvsr5OSkvTdd9/J399foaGhGj58uJ5//nlVq1ZNERERGjNmjEJCQtSpUydJUs2aNdWmTRv169dPs2fP1rlz5zRkyBB1796dOzEAAADAuWH3m2++UcuWLe2v4+LiJEmxsbGaP3++Ro4cqYyMDPXv31+pqalq2rSpVq5cKS8vL/s6Cxcu1JAhQ9S6dWu5ubmpS5cumjZt2nXfFgAAALgem2VZlrOLcLb09HT5+voqLS2NIQ3ADSJ89ApnlwBcNwcmtXN2CcBVK2xec9kxuwAAAMDVIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADBWocJuq1atlJqamqc9PT1drVq1utqaAAAAgCJRqLC7fv16ZWVl5Wk/e/as/ve//111UQAAAEBRcL+ShXfu3Gn/908//aTk5GT76+zsbK1cuVI33XRT0VUHAAAAXIUrCru33367bDabbDZbvsMVvL29NX369CIrDgAAALgaVxR2k5KSZFmWqlSpoq1btyowMNA+z8PDQxUqVFCJEiWKvEgAAACgMK4o7IaFhUmScnJyrkkxAAAAQFG6orB7ob1792rdunU6duxYnvA7duzYqy4MAAAAuFqFCrtvvvmmBg4cqPLlyys4OFg2m80+z2azEXYBAADgEgoVdp9//nm98MILGjVqVFHXAwAAABSZQt1n99SpU+rWrVtR1wIAAAAUqUKF3W7duumLL74o6loAAACAIlWoYQw333yzxowZo82bN6tOnToqWbKkw/xhw4YVSXEAAADA1bBZlmVd6UoREREX79Bm0/79+6+qqOstPT1dvr6+SktLk4+Pj7PLAXAdhI9e4ewSgOvmwKR2zi4BuGqFzWuFOrOblJRUmNUAAACA66rQ99kFgILiLCoAwFkKFXYfffTRS86fO3duoYoBAAAAilKhwu6pU6ccXp87d04//vijUlNT1apVqyIpDAAAALhahQq7S5cuzdOWk5OjgQMHqmrVqlddFAAAAFAUCnWf3Xw7cnNTXFycpkyZUlRdAgAAAFelyMKuJP366686f/58UXYJAAAAFFqhhjHExcU5vLYsS7///rtWrFih2NjYIikMAAAAuFqFOrP77bffOkw7d+6UJL3yyiuaOnVqkRWXnZ2tMWPGKCIiQt7e3qpatar+/e9/68LnYFiWpbFjx6pixYry9vZWdHS09u7dW2Q1AAAAoPgq1JnddevWFXUd+XrxxRc1a9YsLViwQLVr19Y333yjPn36yNfX1/5I4pdeeknTpk3TggULFBERoTFjxigmJkY//fSTvLy8rkudAAAAcE1X9VCJ48ePa8+ePZKk6tWrKzAwsEiKyrVp0yZ17NhR7dr9/ZjD8PBwvfvuu9q6daukv8/qTp06Vc8++6w6duwoSXr77bcVFBSkZcuWqXv37kVaDwAAAIqXQg1jyMjI0KOPPqqKFSuqWbNmatasmUJCQtS3b1+dOXOmyIq78847tWbNGv3yyy+SpO+//15fffWV2rZtK+nvxxYnJycrOjravo6vr68iIyOVmJh40X4zMzOVnp7uMAEAAMA8hQq7cXFx2rBhg5YvX67U1FSlpqbq448/1oYNG/Tkk08WWXGjR49W9+7dVaNGDZUsWVL16tXT8OHD1aNHD0lScnKyJCkoKMhhvaCgIPu8/CQkJMjX19c+Va5cuchqBgAAgOsoVNhdsmSJ3nrrLbVt21Y+Pj7y8fHRvffeqzfffFOLFy8usuI++OADLVy4UIsWLdKOHTu0YMECTZ48WQsWLLiqfuPj45WWlmafDh8+XEQVAwAAwJUUaszumTNn8pxNlaQKFSoU6TCGp59+2n52V5Lq1KmjgwcPKiEhQbGxsQoODpYkpaSkqGLFivb1UlJSdPvtt1+0X09PT3l6ehZZnQAAAHBNhQq7UVFRGjdunN5++237HQ/++usvjR8/XlFRUUVW3JkzZ+Tm5njyuUSJEsrJyZEkRUREKDg4WGvWrLGH2/T0dG3ZskUDBw4ssjoAACjOwkevcHYJ+TowqZ2zS8ANoFBhd+rUqWrTpo0qVaqkunXrSvr74jFPT0998cUXRVZc+/bt9cILLyg0NFS1a9fWt99+q1dffVWPPvqoJMlms2n48OF6/vnnVa1aNfutx0JCQtSpU6ciqwMAAADFU6HCbp06dbR3714tXLhQu3fvliQ99NBD6tGjh7y9vYusuOnTp2vMmDEaNGiQjh07ppCQED3++OMaO3asfZmRI0cqIyND/fv3V2pqqpo2baqVK1dyj10AAADIZl34OLICSkhIUFBQkP0Ma665c+fq+PHjGjVqVJEVeD2kp6fL19dXaWlp8vHxcXY5gHFc9StUAM7FMAZcicLmtULdjeGNN95QjRo18rTXrl1bs2fPLkyXAAAAQJErVNhNTk52uPtBrsDAQP3+++9XXRQAAABQFAoVditXrqyvv/46T/vXX3+tkJCQqy4KAAAAKAqFukCtX79+Gj58uM6dO6dWrVpJktasWaORI0cW6RPUAAAAgKtRqLD79NNP648//tCgQYOUlZUlSfLy8tKoUaMUHx9fpAUCAAAAhVWosGuz2fTiiy9qzJgx+vnnn+Xt7a1q1arxVDIAAAC4lEKF3VxlypRRw4YNi6oWAAAAoEgV6gI1AAAAoDgg7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAY7k7uwAARSd89ApnlwAAgEvhzC4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIzl8mH3yJEjeuSRRxQQECBvb2/VqVNH33zzjX2+ZVkaO3asKlasKG9vb0VHR2vv3r1OrBgAAACuwqXD7qlTp9SkSROVLFlSn3/+uX766Se98sorKleunH2Zl156SdOmTdPs2bO1ZcsWlS5dWjExMTp79qwTKwcAAIArcHd2AZfy4osvqnLlypo3b569LSIiwv5vy7I0depUPfvss+rYsaMk6e2331ZQUJCWLVum7t27X/eaAQAA4Dpc+szuJ598ogYNGqhbt26qUKGC6tWrpzfffNM+PykpScnJyYqOjra3+fr6KjIyUomJiRftNzMzU+np6Q4TAAAAzOPSYXf//v2aNWuWqlWrplWrVmngwIEaNmyYFixYIElKTk6WJAUFBTmsFxQUZJ+Xn4SEBPn6+tqnypUrX7uNAAAAgNO4dNjNycnRHXfcoYkTJ6pevXrq37+/+vXrp9mzZ19Vv/Hx8UpLS7NPhw8fLqKKAQAA4EpcOuxWrFhRtWrVcmirWbOmDh06JEkKDg6WJKWkpDgsk5KSYp+XH09PT/n4+DhMAAAAMI9Lh90mTZpoz549Dm2//PKLwsLCJP19sVpwcLDWrFljn5+enq4tW7YoKirqutYKAAAA1+PSd2MYMWKE7rzzTk2cOFEPPPCAtm7dqjlz5mjOnDmSJJvNpuHDh+v5559XtWrVFBERoTFjxigkJESdOnVybvEAAABwOpcOuw0bNtTSpUsVHx+vCRMmKCIiQlOnTlWPHj3sy4wcOVIZGRnq37+/UlNT1bRpU61cuVJeXl5OrBwAAACuwGZZluXsIpwtPT1dvr6+SktLY/wuirXw0SucXQIAFNiBSe2cXQKKkcLmNZceswsAAABcDcIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwVrEKu5MmTZLNZtPw4cPtbWfPntXgwYMVEBCgMmXKqEuXLkpJSXFekQAAAHAZxSbsbtu2TW+88YZuu+02h/YRI0Zo+fLl+vDDD7VhwwYdPXpU999/v5OqBAAAgCspFmH39OnT6tGjh958802VK1fO3p6Wlqa33npLr776qlq1aqX69etr3rx52rRpkzZv3uzEigEAAOAKikXYHTx4sNq1a6fo6GiH9u3bt+vcuXMO7TVq1FBoaKgSExOvd5kAAABwMe7OLuBy3nvvPe3YsUPbtm3LMy85OVkeHh7y8/NzaA8KClJycvJF+8zMzFRmZqb9dXp6epHVCwAAANfh0md2Dx8+rCeeeEILFy6Ul5dXkfWbkJAgX19f+1S5cuUi6xsAAACuw6XD7vbt23Xs2DHdcccdcnd3l7u7uzZs2KBp06bJ3d1dQUFBysrKUmpqqsN6KSkpCg4Ovmi/8fHxSktLs0+HDx++xlsCAAAAZ3DpYQytW7fWDz/84NDWp08f1ahRQ6NGjVLlypVVsmRJrVmzRl26dJEk7dmzR4cOHVJUVNRF+/X09JSnp+c1rR0AAADO59Jht2zZsrr11lsd2kqXLq2AgAB7e9++fRUXFyd/f3/5+Pho6NChioqKUuPGjZ1RMgAAAFyIS4fdgpgyZYrc3NzUpUsXZWZmKiYmRjNnznR2WQAAAHABNsuyLGcX4Wzp6eny9fVVWlqafHx8nF0OUGjho1c4uwQAKLADk9o5uwQUI4XNay59gRoAAABwNQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxnJ3dgEAAODGFD56hbNLyOPApHbOLgFFjDO7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwlsuH3YSEBDVs2FBly5ZVhQoV1KlTJ+3Zs8dhmbNnz2rw4MEKCAhQmTJl1KVLF6WkpDipYgAAALgKlw+7GzZs0ODBg7V582atXr1a586d0z333KOMjAz7MiNGjNDy5cv14YcfasOGDTp69Kjuv/9+J1YNAAAAV2CzLMtydhFX4vjx46pQoYI2bNigZs2aKS0tTYGBgVq0aJG6du0qSdq9e7dq1qypxMRENW7c+LJ9pqeny9fXV2lpafLx8bnWmwBDhI9e4ewSAABF7MCkds4uARdR2Lzm8md2/yktLU2S5O/vL0navn27zp07p+joaPsyNWrUUGhoqBITE51SIwAAAFyDu7MLuBI5OTkaPny4mjRpoltvvVWSlJycLA8PD/n5+TksGxQUpOTk5Hz7yczMVGZmpv11enr6NasZAAAAzlOswu7gwYP1448/6quvvrqqfhISEjR+/PgiqgoAAJjCVYeoMbyi8IrNMIYhQ4bo008/1bp161SpUiV7e3BwsLKyspSamuqwfEpKioKDg/PtKz4+Xmlpafbp8OHD17J0AAAAOInLh13LsjRkyBAtXbpUa9euVUREhMP8+vXrq2TJklqzZo29bc+ePTp06JCioqLy7dPT01M+Pj4OEwAAAMzj8sMYBg8erEWLFunjjz9W2bJl7eNwfX195e3tLV9fX/Xt21dxcXHy9/eXj4+Phg4dqqioqALdiQEAAADmcvmwO2vWLElSixYtHNrnzZun3r17S5KmTJkiNzc3denSRZmZmYqJidHMmTOvc6UAAABwNS4fdgtyG2AvLy/NmDFDM2bMuA4VAQAAoLhw+TG7AAAAQGERdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADCWu7MLAC4nfPQKZ5cAAACKKc7sAgAAwFiEXQAAABiLsAsAAABjEXYBAABgLMIuAAAAjEXYBQAAgLEIuwAAADAWYRcAAADGIuwCAADAWIRdAAAAGIuwCwAAAGMRdgEAAGAswi4AAACMRdgFAACAsQi7AAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAY7k7u4AbVfjoFc4uIV8HJrVzdgkAAOAfyA2Fx5ldAAAAGIuwCwAAAGMRdgEAAGAsY8LujBkzFB4eLi8vL0VGRmrr1q3OLgkAAABOZsQFau+//77i4uI0e/ZsRUZGaurUqYqJidGePXtUoUIFZ5dXrLjqAHgAAIDCMOLM7quvvqp+/fqpT58+qlWrlmbPnq1SpUpp7ty5zi4NAAAATlTsz+xmZWVp+/btio+Pt7e5ubkpOjpaiYmJ+a6TmZmpzMxM++u0tDRJUnp6+rUt9gI5mWeu23sBAABcC9czO+W+l2VZV7ResQ+7J06cUHZ2toKCghzag4KCtHv37nzXSUhI0Pjx4/O0V65c+ZrUCAAAYCLfqdf/Pf/44w/5+voWePliH3YLIz4+XnFxcfbXOTk5OnnypAICAmSz2ZxYWcGkp6ercuXKOnz4sHx8fJxdzg2H/e9c7H/nYv87F/vfudj/zpWWlqbQ0FD5+/tf0XrFPuyWL19eJUqUUEpKikN7SkqKgoOD813H09NTnp6eDm1+fn7XqsRrxsfHhx82J2L/Oxf737nY/87F/ncu9r9zubld2SVnxf4CNQ8PD9WvX19r1qyxt+Xk5GjNmjWKiopyYmUAAABwtmJ/ZleS4uLiFBsbqwYNGqhRo0aaOnWqMjIy1KdPH2eXBgAAACcyIuw++OCDOn78uMaOHavk5GTdfvvtWrlyZZ6L1kzh6empcePG5RmKgeuD/e9c7H/nYv87F/vfudj/zlXY/W+zrvT+DQAAAEAxUezH7AIAAAAXQ9gFAACAsQi7AAAAMBZhFwAAAMYi7BZzHTp0UGhoqLy8vFSxYkX17NlTR48edXZZN4QDBw6ob9++ioiIkLe3t6pWrapx48YpKyvL2aXdMF544QXdeeedKlWqVLF8MExxM2PGDIWHh8vLy0uRkZHaunWrs0u6YWzcuFHt27dXSEiIbDabli1b5uySbhgJCQlq2LChypYtqwoVKqhTp07as2ePs8u6YcyaNUu33Xab/UEeUVFR+vzzz6+oD8JuMdeyZUt98MEH2rNnj5YsWaJff/1VXbt2dXZZN4Tdu3crJydHb7zxhnbt2qUpU6Zo9uzZeuaZZ5xd2g0jKytL3bp108CBA51divHef/99xcXFady4cdqxY4fq1q2rmJgYHTt2zNml3RAyMjJUt25dzZgxw9ml3HA2bNigwYMHa/PmzVq9erXOnTune+65RxkZGc4u7YZQqVIlTZo0Sdu3b9c333yjVq1aqWPHjtq1a1eB++DWY4b55JNP1KlTJ2VmZqpkyZLOLueG8/LLL2vWrFnav3+/s0u5ocyfP1/Dhw9Xamqqs0sxVmRkpBo2bKjXX39d0t9PqqxcubKGDh2q0aNHO7m6G4vNZtPSpUvVqVMnZ5dyQzp+/LgqVKigDRs2qFmzZs4u54bk7++vl19+WX379i3Q8pzZNcjJkye1cOFC3XnnnQRdJ0lLS5O/v7+zywCKVFZWlrZv367o6Gh7m5ubm6Kjo5WYmOjEyoDrLy0tTZL4Xe8E2dnZeu+995SRkaGoqKgCr0fYNcCoUaNUunRpBQQE6NChQ/r444+dXdINad++fZo+fboef/xxZ5cCFKkTJ04oOzs7z1Mpg4KClJyc7KSqgOsvJydHw4cPV5MmTXTrrbc6u5wbxg8//KAyZcrI09NTAwYM0NKlS1WrVq0Cr0/YdUGjR4+WzWa75LR792778k8//bS+/fZbffHFFypRooR69eolRqcU3pXuf0k6cuSI2rRpo27duqlfv35OqtwMhdn/AHA9DB48WD/++KPee+89Z5dyQ6levbq+++47bdmyRQMHDlRsbKx++umnAq/PmF0XdPz4cf3xxx+XXKZKlSry8PDI0/7bb7+pcuXK2rRp0xWd4sf/u9L9f/ToUbVo0UKNGzfW/Pnz5ebG/yGvRmGOf8bsXltZWVkqVaqUFi9e7DBONDY2VqmpqXybdJ0xZtc5hgwZoo8//lgbN25URESEs8u5oUVHR6tq1ap64403CrS8+zWuB4UQGBiowMDAQq2bk5MjScrMzCzKkm4oV7L/jxw5opYtW6p+/fqaN28eQbcIXM3xj2vDw8ND9evX15o1a+wBKycnR2vWrNGQIUOcWxxwjVmWpaFDh2rp0qVav349QdcF5OTkXFHOIewWY1u2bNG2bdvUtGlTlStXTr/++qvGjBmjqlWrclb3Ojhy5IhatGihsLAwTZ48WcePH7fPCw4OdmJlN45Dhw7p5MmTOnTokLKzs/Xdd99Jkm6++WaVKVPGucUZJi4uTrGxsWrQoIEaNWqkqVOnKiMjQ3369HF2aTeE06dPa9++ffbXSUlJ+u677+Tv76/Q0FAnVma+wYMHa9GiRfr4449VtmxZ+zh1X19feXt7O7k688XHx6tt27YKDQ3Vn3/+qUWLFmn9+vVatWpVwTuxUGzt3LnTatmypeXv7295enpa4eHh1oABA6zffvvN2aXdEObNm2dJynfC9REbG5vv/l+3bp2zSzPS9OnTrdDQUMvDw8Nq1KiRtXnzZmeXdMNYt25dvsd6bGyss0sz3sV+z8+bN8/Zpd0QHn30USssLMzy8PCwAgMDrdatW1tffPHFFfXBmF0AAAAYiwGGAAAAMBZhFwAAAMYi7AIAAMBYhF0AAAAYi7ALAAAAYxF2AQAAYCzCLgAAAIxF2AUAAICxCLsAAAAwFmEXAAAAxiLsAgAAwFiEXQAAABjr/wBmKNGJv3OqRAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = np.random.normal(loc=0, scale=1, size=1_000)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.hist(a, color=\"C0\", bins=15)\n",
    "ax.set(title=\"Samples from a normal distribution using numpy\", ylabel=\"count\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to do it in PyTensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TensorType(float64, shape=())"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = pt.random.normal(loc=0, scale=1, name=\"y\")\n",
    "y.type"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we show the graph using {func}`~pytensor.dprint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normal_rv{\"(),()->()\"}.1 [id A] 'y'\n",
      " ├─ RNG(<Generator(PCG64) at 0x7FC6504CC740>) [id B]\n",
      " ├─ NoneConst{None} [id C]\n",
      " ├─ 0 [id D]\n",
      " └─ 1 [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(y)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inputs are always in the following order:\n",
    "1. `rng` shared variable\n",
    "2. `size`\n",
    "4. `arg1`\n",
    "5. `arg2`\n",
    "6. `argn`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We *could* sample by calling {meth}`~pytensor.graph.basic.Variable.eval`. on the random variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(0.67492335)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.eval()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note however that these samples are always the same!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0: 0.6749233482557402\n",
      "Sample 1: 0.6749233482557402\n",
      "Sample 2: 0.6749233482557402\n",
      "Sample 3: 0.6749233482557402\n",
      "Sample 4: 0.6749233482557402\n",
      "Sample 5: 0.6749233482557402\n",
      "Sample 6: 0.6749233482557402\n",
      "Sample 7: 0.6749233482557402\n",
      "Sample 8: 0.6749233482557402\n",
      "Sample 9: 0.6749233482557402\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(f\"Sample {i}: {y.eval()}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We always get the same samples! This has to do with the random seed step in the graph, i.e. `RandomGeneratorSharedVariable` (we will not go deeper into this subject here). We will show how to generate different samples with `pymc` below."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyMC\n",
    "\n",
    "![pymc logo](https://raw.githubusercontent.com/pymc-devs/pymc/main/docs/logos/PyMC.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do so, we start by defining a `pymc` normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normal_rv{\"(),()->()\"}.1 [id A]\n",
      " ├─ RNG(<Generator(PCG64) at 0x7FC650553D80>) [id B]\n",
      " ├─ NoneConst{None} [id C]\n",
      " ├─ 0 [id D]\n",
      " └─ 1 [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = pm.Normal.dist(mu=0, sigma=1)\n",
    "pytensor.dprint(x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that `x` is just a normal `RandomVariable` and which is the same as `y` except for the `rng`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try to generate samples by calling {meth}`~pytensor.graph.basic.Variable.eval` as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0: 0.3880069666747013\n",
      "Sample 1: 0.3880069666747013\n",
      "Sample 2: 0.3880069666747013\n",
      "Sample 3: 0.3880069666747013\n",
      "Sample 4: 0.3880069666747013\n",
      "Sample 5: 0.3880069666747013\n",
      "Sample 6: 0.3880069666747013\n",
      "Sample 7: 0.3880069666747013\n",
      "Sample 8: 0.3880069666747013\n",
      "Sample 9: 0.3880069666747013\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(f\"Sample {i}: {x.eval()}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before we get the same value for all iterations. The correct way to generate random samples is using {func}`~pymc.draw`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAIQCAYAAACbhEYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABC5klEQVR4nO3deVjU5f7/8dcgsggyiCJIIqJ53DIzV9LjSuGuuZRlhktabmW2GJ3j1iJlm2nuv1w6aZ6sXDu55FqJS1pZluaWmgZqyqCYoPL5/dHFfB1BBQQGb5+P65rrcu7PMu/P4vjynnvusVmWZQkAAAAwgIe7CwAAAADyC+EWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYoImw2m8aMGVOor5mUlKRu3bqpdOnSstlsmjBhQqG+Pm7Mb7/9JpvNpjlz5uRp+zlz5shms+m3335ztjVv3lzNmzfPl/qu58p7fsyYMbLZbDp58mShvH7FihXVu3fvQnmtgnKj9wBgIsItjPLjjz+qW7duioiIkI+Pj2677Tbde++9mjRpkrtLK5KefvpprVy5UnFxcfrPf/6j1q1bu7sk3IQ2bdqkMWPGKDk52d2lZFGUawNQMDzdXQCQXzZt2qQWLVqoQoUK6t+/v0JDQ3XkyBFt3rxZ7777roYOHeruEouctWvXqlOnTnr22WfdXQqKiFWrVuV6m02bNmns2LHq3bu3AgMDc7zdX3/9JU/Pgv1n6Fq17dmzRx4eN3cfT0REhP766y8VL17c3aUARQbhFsZ49dVXZbfbtW3btiz/iB0/ftw9RRVxx48fz1EYSU1NlZ+fX8EXdJM7f/68vLy8burA5OXlVaD7z8jIUHp6unx8fOTj41Ogr3U93t7ebn39/GCz2dx+HoGi5uZ9BwausH//ftWsWTPbsFa2bFmX57Nnz1bLli1VtmxZeXt7q0aNGpo6dWqW7SpWrKj27dtr/fr1qlevnnx9fVWrVi2tX79ekvTZZ5+pVq1a8vHxUd26dfXdd9+5bN+7d2/5+/vrwIEDiomJkZ+fn8LCwvTSSy/JsqzrHtPRo0fVt29fhYSEyNvbWzVr1tSsWbOyrDdp0iTVrFlTJUqUUKlSpVSvXj3Nnz//qvvNHGtpWZYmT54sm80mm83msmzDhg0aNGiQypYtq/Llyzu3nTJlimrWrClvb2+FhYVp8ODBWT7ybd68ue644w7t3LlTzZo1U4kSJXT77bfrk08+kSRt2LBBDRs2lK+vr6pWraovv/zyuuciPT1do0aNUt26dWW32+Xn56d//vOfWrdu3XW3lf7vWn799ddq0KCBfHx8VKlSJX3wwQdZ1j1w4IC6d++uoKAglShRQo0aNdLnn3/uss769etls9m0YMEC/fvf/9Ztt92mEiVKKCUlxXndDx8+rPbt28vf31+33XabJk+eLOnv4TMtW7aUn5+fIiIislyrU6dO6dlnn1WtWrXk7++vgIAAtWnTRj/88EOOjjU7u3btUsuWLeXr66vy5cvrlVdeUUZGRpb1shtze637a8yYMXruueckSZGRkc57KXMcr81m05AhQzRv3jznfbNixQrnsuzGmZ88eVIPPPCAAgICVLp0aT311FM6f/68c/m1xplevs/r1ZbdmNvcXPuPP/5Yr776qsqXLy8fHx+1atVK+/bty1LTlXr37q2KFStmac8cc3y51atXq0mTJgoMDJS/v7+qVq2qF1988ZrnIvP+O3r0qDp37ix/f38FBwfr2Wef1aVLl1z2/+eff6pXr14KCAhQYGCgYmNj9cMPP+RoHG/me8XGjRv1+OOPq3Tp0goICNCjjz6q06dPO9eLjY1VmTJldOHChSz7uO+++1S1alXn88z7ZeHChapRo4Z8fX0VFRWlH3/8UZI0ffp03X777fLx8VHz5s1dxotn2rJli9q2batSpUrJz89Pd955p959991rHgvMQs8tjBEREaGEhAT99NNPuuOOO6657tSpU1WzZk117NhRnp6eWrZsmQYNGqSMjAwNHjzYZd19+/bp4Ycf1uOPP65HHnlEb775pjp06KBp06bpxRdf1KBBgyRJ8fHxeuCBB7J81Hnp0iW1bt1ajRo10vjx47VixQqNHj1aFy9e1EsvvXTVGpOSktSoUSPnm31wcLC++OIL9evXTykpKRo2bJgkaebMmXryySfVrVs3ZwjYuXOntmzZoocffjjbfTdt2lT/+c9/1KtXL91777169NFHs6wzaNAgBQcHa9SoUUpNTZX09z++Y8eOVXR0tAYOHKg9e/Zo6tSp2rZtm7755huXj0ZPnz6t9u3bq0ePHurevbumTp2qHj16aN68eRo2bJieeOIJPfzww3rjjTfUrVs3HTlyRCVLlrzq+UhJSdH/+3//Tw899JD69++vM2fO6P3331dMTIy2bt2qu+6666rbZtq3b5+6deumfv36KTY2VrNmzVLv3r1Vt25d1axZ03ne77nnHp07d05PPvmkSpcurblz56pjx4765JNPdP/997vs8+WXX5aXl5eeffZZpaWlOXs+L126pDZt2qhp06YaP3685s2bpyFDhsjPz0//+te/1LNnT3Xp0kXTpk3To48+qqioKEVGRkr6O2AtXrxY3bt3V2RkpJKSkjR9+nQ1a9ZMP//8s8LCwq57rJdLTExUixYtdPHiRb3wwgvy8/PTjBkz5Ovre91tr3d/denSRb/++qs++ugjvfPOOypTpowkKTg42LmPtWvX6uOPP9aQIUNUpkyZbIPd5R544AFVrFhR8fHx2rx5syZOnKjTp09n+x+Ra8lJbZfL7bV/7bXX5OHhoWeffVYOh0Pjx49Xz549tWXLllzVeTW7du1S+/btdeedd+qll16St7e39u3bp2+++ea62166dEkxMTFq2LCh3nzzTX355Zd66623VLlyZQ0cOFDS373oHTp00NatWzVw4EBVq1ZNS5YsUWxsbK7qHDJkiAIDAzVmzBjne8KhQ4ec/wno1auXPvjgA61cuVLt27d3bpeYmKi1a9dq9OjRLvv76quvtHTpUud7cXx8vNq3b6/nn39eU6ZM0aBBg3T69GmNHz9effv21dq1a53brl69Wu3bt1e5cuX01FNPKTQ0VL/88ouWL1+up556KlfHhZuYBRhi1apVVrFixaxixYpZUVFR1vPPP2+tXLnSSk9Pz7LuuXPnsrTFxMRYlSpVcmmLiIiwJFmbNm1ytq1cudKSZPn6+lqHDh1ytk+fPt2SZK1bt87ZFhsba0myhg4d6mzLyMiw2rVrZ3l5eVknTpxwtkuyRo8e7Xzer18/q1y5ctbJkyddaurRo4dlt9udx9CpUyerZs2a1zk72ZNkDR482KVt9uzZliSrSZMm1sWLF53tx48ft7y8vKz77rvPunTpkrP9vffesyRZs2bNcrY1a9bMkmTNnz/f2bZ7925LkuXh4WFt3rzZ2Z55PmfPnn3NWi9evGilpaW5tJ0+fdoKCQmx+vbte91jzbyWGzdudDkmb29v65lnnnG2DRs2zJJkffXVV862M2fOWJGRkVbFihWdx75u3TpLklWpUqUs91PmdR83bpxLrb6+vpbNZrMWLFiQ5bxcfu3Pnz/vco4ty7IOHjxoeXt7Wy+99JJLW07OXeYxbdmyxeXY7Xa7Jck6ePCgs71Zs2ZWs2bNnM9zcn+98cYbWfaTKfOa79q1K9tllx/36NGjLUlWx44dXdYbNGiQJcn64YcfLMu69nFfuc9r1RYREWHFxsY6n+f22levXt3lnnz33XctSdaPP/6Y5bUuFxsba0VERGRpzzz+TO+8844lyeV94krZnYvM++/ye8WyLKtOnTpW3bp1nc8//fRTS5I1YcIEZ9ulS5esli1b5ui+ynyvqFu3rsv77Pjx4y1J1pIlS5z7LF++vPXggw+6bP/2229bNpvNOnDggLNNkuXt7e1yvTLfW0NDQ62UlBRne1xcnMu1vXjxohUZGWlFRERYp0+fdnmtjIyMax4LzMKwBBjj3nvvVUJCgjp27KgffvhB48ePV0xMjG677TYtXbrUZd3Le6wcDodOnjypZs2a6cCBA3I4HC7r1qhRQ1FRUc7nDRs2lCS1bNlSFSpUyNJ+4MCBLLUNGTLE+efMntj09PSrfhxvWZY+/fRTdejQQZZl6eTJk85HTEyMHA6HduzYIUkKDAzU77//rm3btuXoPOVU//79VaxYMefzL7/8Uunp6Ro2bJhLz3T//v0VEBCQ5aNbf39/9ejRw/m8atWqCgwMVPXq1Z3nSrr2ebtcsWLFnL2iGRkZOnXqlC5evKh69eo5z8X11KhRQ//85z+dz4ODg1W1alWX1/7f//6nBg0aqEmTJi7HMmDAAP3222/6+eefXfYZGxt71R7Qxx57zPnnwMBAVa1aVX5+fnrggQec7Znn5fIavL29nef40qVL+vPPP50fSef0WC/3v//9T40aNVKDBg1cjr1nz57X3TY/7q9mzZqpRo0aOV7/yk9PMr8M+r///S/PNeREbq99nz59XMYoZ95b17uXcypziNWSJUuyHUJyPU888YTL83/+858uta1YsULFixdX//79nW0eHh5Zzv/1DBgwwOVTm4EDB8rT09N5vTw8PNSzZ08tXbpUZ86cca43b9483XPPPc5PLDK1atXKpXc/8z2ia9euLp/uXPne8d133+ngwYMaNmxYluFpVw73gNkItzBK/fr19dlnn+n06dPaunWr4uLidObMGXXr1s3lH6ZvvvlG0dHR8vPzU2BgoIKDg53j2K4Mt5cHWEmy2+2SpPDw8GzbLx9rJv39xl6pUiWXtn/84x+SlO14MUk6ceKEkpOTNWPGDAUHB7s8+vTpI+n/viQ3YsQI+fv7q0GDBqpSpYoGDx6co48tr+fKf3AOHTokSS7j46S/v4BUqVIl5/JM5cuXz/IPit1uz/F5y87cuXN15513ysfHR6VLl1ZwcLA+//zzLNfsaq68lpJUqlQpl9c+dOhQlmOUpOrVqzuXX+7K85TJx8cny8ffdrv9qufl8hoyMjL0zjvvqEqVKvL29laZMmUUHBysnTt35vhYL3fo0CFVqVIlS3t2x3ml/Li/rnaOrubKWitXriwPD4+r/n3JL7m99lfeT6VKlZKUs3s5Jx588EE1btxYjz32mEJCQtSjRw99/PHHOQq62d1/2d3r5cqVU4kSJVzWu/3223NV55XXy9/fX+XKlXO5Xo8++qj++usvLVq0SNLfM1Vs375dvXr1yrK/vL7n7t+/X5KuOywN5iPcwkheXl6qX7++xo0bp6lTp+rChQtauHChpL/fAFu1aqWTJ0/q7bff1ueff67Vq1fr6aeflqQs/3Bc3nuZk3YrB18Uu57MGh555BGtXr0620fjxo0l/f0P7549e7RgwQI1adJEn376qZo0aZJlHFtu5WQ85rXk93n78MMP1bt3b1WuXFnvv/++VqxYodWrV6tly5Y57tUqiGt2tfN0I8c/btw4DR8+XE2bNtWHH36olStXavXq1apZs2aeevBuRH7cXzd6L135n4Gr9cJd+WWpgpbX+ymn9fv6+mrjxo368ssv1atXL+3cuVMPPvig7r333use69Vqc5caNWqobt26+vDDDyX9/ffZy8vL5VOMTO54z4VZ+EIZjFevXj1J0h9//CFJWrZsmdLS0rR06VKXHoKcfus+tzIyMnTgwAFnb60k/frrr5J01S/WBAcHq2TJkrp06ZKio6Ov+xp+fn568MEH9eCDDyo9PV1dunTRq6++qri4uHybJigiIkLS3z0ul/dEp6en6+DBgzmq80Z88sknqlSpkj777DOXcHCjIf5KERER2rNnT5b23bt3O5cXtE8++UQtWrTQ+++/79KenJzs/FJUbkRERGjv3r1Z2rM7zuxc7/7K74989+7d69Lbu2/fPmVkZDj/vmT2kF45S8eVPatS7j6OLqxrX6pUqWx/VCK7+j08PNSqVSu1atVKb7/9tsaNG6d//etfWrdu3Q3/nYuIiNC6det07tw5l97bnMz4cLm9e/eqRYsWzudnz57VH3/8obZt27qs9+ijj2r48OH6448/NH/+fLVr1855LfND5cqVJUk//fRTgb8foWij5xbGWLduXbb/g88c95X5cWPm//4vX9fhcGj27NkFVtt7773n/LNlWXrvvfdUvHhxtWrVKtv1ixUrpq5du+rTTz/VTz/9lGX5iRMnnH/+888/XZZ5eXmpRo0asiwr26l38io6OlpeXl6aOHGiy7l7//335XA41K5du3x7rexkd922bNmihISEfH2dtm3bauvWrS77TU1N1YwZM1SxYsVcjR3Nq2LFimW5lxcuXKijR4/maX9t27bV5s2btXXrVmfbiRMnNG/evOtum5P7K3MO5Pz6FbDMKdMyZf7CYJs2bSRJAQEBKlOmjDZu3Oiy3pQpU7LsKze1Fda1r1y5shwOh3bu3Ols++OPP5wf2Wc6depUlm0zZwVJS0u74TpiYmJ04cIFzZw509mWkZGR5fxfz4wZM1zea6ZOnaqLFy86r1emhx56SDabTU899ZQOHDigRx555MYO4Ap33323IiMjNWHChCzXm97dWws9tzDG0KFDde7cOd1///2qVq2a0tPTtWnTJv33v/9VxYoVnWNV77vvPnl5ealDhw56/PHHdfbsWc2cOVNly5Z19u7mJx8fH61YsUKxsbFq2LChvvjiC33++ed68cUXrzolkfT3NEPr1q1Tw4YN1b9/f9WoUUOnTp3Sjh079OWXXzr/4bvvvvsUGhqqxo0bKyQkRL/88ovee+89tWvX7ppTa+VWcHCw4uLiNHbsWLVu3VodO3bUnj17NGXKFNWvXz/f/6G6Uvv27fXZZ5/p/vvvV7t27XTw4EFNmzZNNWrU0NmzZ/PtdV544QV99NFHatOmjZ588kkFBQVp7ty5OnjwoD799NNC+YGG9u3b66WXXlKfPn10zz336Mcff9S8efOyjN3Oqeeff97588pPPfWUcyqwiIgIl4CVnZzcX3Xr1pUk/etf/1KPHj1UvHhxdejQIc8//HHw4EF17NhRrVu3VkJCgj788EM9/PDDql27tnOdxx57TK+99poee+wx1atXTxs3bnR+InK53NRWWNe+R48eGjFihO6//349+eSTOnfunKZOnap//OMfLl8YfOmll7Rx40a1a9dOEREROn78uKZMmaLy5cu7fOktrzp37qwGDRromWee0b59+1StWjUtXbrU+d6S017v9PR0tWrVyjkV4pQpU9SkSRN17NjRZb3g4GC1bt1aCxcuVGBgYL7/h9jDw0NTp05Vhw4ddNddd6lPnz4qV66cdu/erV27dmnlypX5+noowgp9fgaggHzxxRdW3759rWrVqln+/v6Wl5eXdfvtt1tDhw61kpKSXNZdunSpdeedd1o+Pj5WxYoVrddff92aNWtWlimDIiIirHbt2mV5LWUzhVbmlDxvvPGGsy02Ntby8/Oz9u/fb913331WiRIlrJCQEGv06NFZpnrSFVMYWZZlJSUlWYMHD7bCw8Ot4sWLW6GhoVarVq2sGTNmONeZPn261bRpU6t06dKWt7e3VblyZeu5556zHA7Hdc9ZdseROb3Ptm3bst3mvffes6pVq2YVL17cCgkJsQYOHJhl2p1mzZplO31Ubs7nlTIyMqxx48ZZERERlre3t1WnTh1r+fLlV51WKaevfeXUV5ZlWfv377e6detmBQYGWj4+PlaDBg2s5cuXu6yTOR3UwoULs+wz87pn91o5OS/nz5+3nnnmGatcuXKWr6+v1bhxYyshISFLrTmdCsyyLGvnzp1Ws2bNLB8fH+u2226zXn75Zev999+/7lRgOb2/Xn75Zeu2226zPDw8XPZ5rWt75T2fORXWzz//bHXr1s0qWbKkVapUKWvIkCHWX3/95bLtuXPnrH79+ll2u90qWbKk9cADD1jHjx/P9u/R1Wq7ciowy7qxa5+b67Fq1SrrjjvusLy8vKyqVataH374YZapwNasWWN16tTJCgsLs7y8vKywsDDroYcesn799ddrvubV7r8r929ZlnXixAnr4YcftkqWLGnZ7Xard+/e1jfffGNJcpmyLjuZ7xUbNmywBgwYYJUqVcry9/e3evbsaf3555/ZbvPxxx9bkqwBAwZkuzyn762WdfXr8PXXX1v33nuvVbJkScvPz8+68847rUmTJl3zWGAWm2XRVw8UlN69e+uTTz7J155FAChIixcv1v3336+vv/7a+cXV7MyZM0d9+vTRtm3bnN9tuJ4lS5aoc+fO2rhxo8u0fEB+YswtAAC3qL/++svl+aVLlzRp0iQFBATo7rvvzvfXmzlzpipVqpQvwyqAq2HMLQAAt6ihQ4fqr7/+UlRUlNLS0vTZZ59p06ZNGjdu3A1P4Xa5BQsWaOfOnfr888/17rvv8qMKKFCEWwAAblEtW7bUW2+9peXLl+v8+fO6/fbbNWnSJJdfVcwPDz30kPz9/dWvXz8NGjQoX/cNXIkxtwAAADAGY24BAABgDMItAAAAjMGYW/39iyzHjh1TyZIlGeQOAABQBFmWpTNnzigsLOyaP6pCuJV07NgxhYeHu7sMAAAAXMeRI0dUvnz5qy4n3ErOn5A8cuSIAgIC3FwNAAAArpSSkqLw8PDr/rQ84Vb/9/vZAQEBhFsAAIAi7HpDSPlCGQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjOHWcLtx40Z16NBBYWFhstlsWrx4cZZ1fvnlF3Xs2FF2u11+fn6qX7++Dh8+7Fx+/vx5DR48WKVLl5a/v7+6du2qpKSkQjwKAAAAFBVuDbepqamqXbu2Jk+enO3y/fv3q0mTJqpWrZrWr1+vnTt3auTIkfLx8XGu8/TTT2vZsmVauHChNmzYoGPHjqlLly6FdQgAAAAoQmyWZVnuLkL6+3eCFy1apM6dOzvbevTooeLFi+s///lPtts4HA4FBwdr/vz56tatmyRp9+7dql69uhISEtSoUaMcvXZKSorsdrscDocCAgJu+FgAAACQv3Ka14rsmNuMjAx9/vnn+sc//qGYmBiVLVtWDRs2dBm6sH37dl24cEHR0dHOtmrVqqlChQpKSEhwQ9UAAABwpyIbbo8fP66zZ8/qtddeU+vWrbVq1Srdf//96tKlizZs2CBJSkxMlJeXlwIDA122DQkJUWJi4lX3nZaWppSUFJcHAAAAbn6e7i7gajIyMiRJnTp10tNPPy1Juuuuu7Rp0yZNmzZNzZo1y/O+4+PjNXbs2HypEwAAAEVHke25LVOmjDw9PVWjRg2X9urVqztnSwgNDVV6erqSk5Nd1klKSlJoaOhV9x0XFyeHw+F8HDlyJN/rBwAAQOErsuHWy8tL9evX1549e1zaf/31V0VEREiS6tatq+LFi2vNmjXO5Xv27NHhw4cVFRV11X17e3srICDA5QEAAICbn1uHJZw9e1b79u1zPj948KC+//57BQUFqUKFCnruuef04IMPqmnTpmrRooVWrFihZcuWaf369ZIku92ufv36afjw4QoKClJAQICGDh2qqKioHM+UAAAAAHO4dSqw9evXq0WLFlnaY2NjNWfOHEnSrFmzFB8fr99//11Vq1bV2LFj1alTJ+e658+f1zPPPKOPPvpIaWlpiomJ0ZQpU645LOFKTAUGAABQtOU0rxWZeW7diXALoMgYY3d3Bdkb43B3BQBucTf9PLcAAABAbhFuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjeLq7AAAA8mSM3d0VZG+Mw90VALc0em4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGIN5bgEA11dU55QFgCvQcwsAAABjEG4BAABgDMItAAAAjOHWcLtx40Z16NBBYWFhstlsWrx48VXXfeKJJ2Sz2TRhwgSX9lOnTqlnz54KCAhQYGCg+vXrp7NnzxZs4QAAACiS3BpuU1NTVbt2bU2ePPma6y1atEibN29WWFhYlmU9e/bUrl27tHr1ai1fvlwbN27UgAEDCqpkAAAAFGFunS2hTZs2atOmzTXXOXr0qIYOHaqVK1eqXbt2Lst++eUXrVixQtu2bVO9evUkSZMmTVLbtm315ptvZhuGAQAAYK4iPeY2IyNDvXr10nPPPaeaNWtmWZ6QkKDAwEBnsJWk6OhoeXh4aMuWLVfdb1pamlJSUlweAAAAuPkV6XD7+uuvy9PTU08++WS2yxMTE1W2bFmXNk9PTwUFBSkxMfGq+42Pj5fdbnc+wsPD87VuAAAAuEeRDbfbt2/Xu+++qzlz5shms+XrvuPi4uRwOJyPI0eO5Ov+AQAA4B5FNtx+9dVXOn78uCpUqCBPT095enrq0KFDeuaZZ1SxYkVJUmhoqI4fP+6y3cWLF3Xq1CmFhoZedd/e3t4KCAhweQAAAODmV2R/frdXr16Kjo52aYuJiVGvXr3Up08fSVJUVJSSk5O1fft21a1bV5K0du1aZWRkqGHDhoVeMwAAANzLreH27Nmz2rdvn/P5wYMH9f333ysoKEgVKlRQ6dKlXdYvXry4QkNDVbVqVUlS9erV1bp1a/Xv31/Tpk3ThQsXNGTIEPXo0YOZEgAAAG5Bbh2W8O2336pOnTqqU6eOJGn48OGqU6eORo0aleN9zJs3T9WqVVOrVq3Utm1bNWnSRDNmzCiokgEAAFCEubXntnnz5rIsK8fr//bbb1nagoKCNH/+/HysCgAAADerIvuFMgAAACC3CLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYw63hduPGjerQoYPCwsJks9m0ePFi57ILFy5oxIgRqlWrlvz8/BQWFqZHH31Ux44dc9nHqVOn1LNnTwUEBCgwMFD9+vXT2bNnC/lIAAAAUBR4uvPFU1NTVbt2bfXt21ddunRxWXbu3Dnt2LFDI0eOVO3atXX69Gk99dRT6tixo7799lvnej179tQff/yh1atX68KFC+rTp48GDBig+fPnF/bhALiZjLG7uwIAQAGwWZZlubsISbLZbFq0aJE6d+581XW2bdumBg0a6NChQ6pQoYJ++eUX1ahRQ9u2bVO9evUkSStWrFDbtm31+++/KywsLEevnZKSIrvdLofDoYCAgPw4HABFHeEWBWWMw90VAEbKaV67qcbcOhwO2Ww2BQYGSpISEhIUGBjoDLaSFB0dLQ8PD23ZsuWq+0lLS1NKSorLAwAAADe/mybcnj9/XiNGjNBDDz3kTOuJiYkqW7asy3qenp4KCgpSYmLiVfcVHx8vu93ufISHhxdo7QAAACgcN0W4vXDhgh544AFZlqWpU6fe8P7i4uLkcDicjyNHjuRDlQAAAHA3t36hLCcyg+2hQ4e0du1alzEWoaGhOn78uMv6Fy9e1KlTpxQaGnrVfXp7e8vb27vAagYAAIB7FOme28xgu3fvXn355ZcqXbq0y/KoqCglJydr+/btzra1a9cqIyNDDRs2LOxyAQAA4GZu7bk9e/as9u3b53x+8OBBff/99woKClK5cuXUrVs37dixQ8uXL9elS5ec42iDgoLk5eWl6tWrq3Xr1urfv7+mTZumCxcuaMiQIerRo0eOZ0oAAACAOdw6Fdj69evVokWLLO2xsbEaM2aMIiMjs91u3bp1at68uaS/f8RhyJAhWrZsmTw8PNS1a1dNnDhR/v7+Oa6DqcCAWxBTgaGgMBUYUCBymtfc2nPbvHlzXStb5yR3BwUF8YMNAAAAkFTEx9wCAAAAuUG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGJ7uLgDALWCM3d0VAABuEfTcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMdwabjdu3KgOHTooLCxMNptNixcvdlluWZZGjRqlcuXKydfXV9HR0dq7d6/LOqdOnVLPnj0VEBCgwMBA9evXT2fPni3EowAAAEBR4dZwm5qaqtq1a2vy5MnZLh8/frwmTpyoadOmacuWLfLz81NMTIzOnz/vXKdnz57atWuXVq9ereXLl2vjxo0aMGBAYR0CAAAAihCbZVmWu4uQJJvNpkWLFqlz586S/u61DQsL0zPPPKNnn31WkuRwOBQSEqI5c+aoR48e+uWXX1SjRg1t27ZN9erVkyStWLFCbdu21e+//66wsLAcvXZKSorsdrscDocCAgIK5PiAW9oYu7srAArPGIe7KwCMlNO8VmTH3B48eFCJiYmKjo52ttntdjVs2FAJCQmSpISEBAUGBjqDrSRFR0fLw8NDW7ZsKfSaAQAA4F6e7i7gahITEyVJISEhLu0hISHOZYmJiSpbtqzLck9PTwUFBTnXyU5aWprS0tKcz1NSUvKrbAAAALhRkQ23BSk+Pl5jx451dxkAABMVxWE4DJXALaTIDksIDQ2VJCUlJbm0JyUlOZeFhobq+PHjLssvXryoU6dOOdfJTlxcnBwOh/Nx5MiRfK4eAAAA7lBkw21kZKRCQ0O1Zs0aZ1tKSoq2bNmiqKgoSVJUVJSSk5O1fft25zpr165VRkaGGjZseNV9e3t7KyAgwOUBAACAm59bhyWcPXtW+/btcz4/ePCgvv/+ewUFBalChQoaNmyYXnnlFVWpUkWRkZEaOXKkwsLCnDMqVK9eXa1bt1b//v01bdo0XbhwQUOGDFGPHj1yPFMCAAAAzOHWcPvtt9+qRYsWzufDhw+XJMXGxmrOnDl6/vnnlZqaqgEDBig5OVlNmjTRihUr5OPj49xm3rx5GjJkiFq1aiUPDw917dpVEydOLPRjAQAAgPsVmXlu3Yl5boECVhS/YAPcSvhCGQxw089zCwAAAOQW4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDHyFG5btmyp5OTkLO0pKSlq2bLljdYEAAAA5Emewu369euVnp6epf38+fP66quvbrgoAAAAIC88c7Pyzp07nX/++eeflZiY6Hx+6dIlrVixQrfddlv+VQcAAADkQq7C7V133SWbzSabzZbt8ANfX19NmjQp34oDAAAAciNX4fbgwYOyLEuVKlXS1q1bFRwc7Fzm5eWlsmXLqlixYvleJAAAAJATuQq3ERERkqSMjIwCKQYAAAC4EbkKt5fbu3ev1q1bp+PHj2cJu6NGjbrhwgAAAIDcylO4nTlzpgYOHKgyZcooNDRUNpvNucxmsxFuAQAA4BZ5CrevvPKKXn31VY0YMSK/6wEAAADyLE/z3J4+fVrdu3fP71oAAACAG5KncNu9e3etWrUqv2sBAAAAbkiehiXcfvvtGjlypDZv3qxatWqpePHiLsuffPLJfCkOAAAAyA2bZVlWbjeKjIy8+g5tNh04cOCGiipsKSkpstvtcjgcCggIcHc5gHnG2N1dAXBrG+NwdwXADctpXstTz+3BgwfzXBgAAABQUPI05hYAAAAoivLUc9u3b99rLp81a1aeigEAAABuRJ7C7enTp12eX7hwQT/99JOSk5PVsmXLfCkMAAAAyK08hdtFixZlacvIyNDAgQNVuXLlGy4KAAAAyIt8G3Pr4eGh4cOH65133smvXQIAAAC5kq9fKNu/f78uXryYn7sEAAAAcixPwxKGDx/u8tyyLP3xxx/6/PPPFRsbmy+FAQAAALmVp57b7777zuWxc+dOSdJbb72lCRMm5Ftxly5d0siRIxUZGSlfX19VrlxZL7/8si7/3QnLsjRq1CiVK1dOvr6+io6O1t69e/OtBgAAANw88tRzu27duvyuI1uvv/66pk6dqrlz56pmzZr69ttv1adPH9ntdudP/I4fP14TJ07U3LlzFRkZqZEjRyomJkY///yzfHx8CqVOAAAAFA15CreZTpw4oT179kiSqlatquDg4HwpKtOmTZvUqVMntWvXTpJUsWJFffTRR9q6daukv3ttJ0yYoH//+9/q1KmTJOmDDz5QSEiIFi9erB49euRrPQAAACja8jQsITU1VX379lW5cuXUtGlTNW3aVGFhYerXr5/OnTuXb8Xdc889WrNmjX799VdJ0g8//KCvv/5abdq0kfT3zwAnJiYqOjrauY3dblfDhg2VkJBw1f2mpaUpJSXF5QEAAICbX57C7fDhw7VhwwYtW7ZMycnJSk5O1pIlS7RhwwY988wz+VbcCy+8oB49eqhatWoqXry46tSpo2HDhqlnz56SpMTERElSSEiIy3YhISHOZdmJj4+X3W53PsLDw/OtZgAAALhPnsLtp59+qvfff19t2rRRQECAAgIC1LZtW82cOVOffPJJvhX38ccfa968eZo/f7527NihuXPn6s0339TcuXNvaL9xcXFyOBzOx5EjR/KpYgAAALhTnsbcnjt3LktvqSSVLVs2X4clPPfcc87eW0mqVauWDh06pPj4eMXGxio0NFSSlJSUpHLlyjm3S0pK0l133XXV/Xp7e8vb2zvf6gQAAEDRkKee26ioKI0ePVrnz593tv31118aO3asoqKi8q24c+fOycPDtcRixYopIyNDkhQZGanQ0FCtWbPGuTwlJUVbtmzJ1zoAAABwc8hTz+2ECRPUunVrlS9fXrVr15b095e9vL29tWrVqnwrrkOHDnr11VdVoUIF1axZU999953efvtt9e3bV5Jks9k0bNgwvfLKK6pSpYpzKrCwsDB17tw53+oAAADAzSFP4bZWrVrau3ev5s2bp927d0uSHnroIfXs2VO+vr75VtykSZM0cuRIDRo0SMePH1dYWJgef/xxjRo1yrnO888/r9TUVA0YMEDJyclq0qSJVqxYwRy3AAAAtyCbdfnPfeVQfHy8QkJCnD2omWbNmqUTJ05oxIgR+VZgYUhJSZHdbpfD4VBAQIC7ywHMM8bu7gqAW9sYh7srAG5YTvNansbcTp8+XdWqVcvSXrNmTU2bNi0vuwQAAABuWJ7CbWJiosvsBJmCg4P1xx9/3HBRAAAAQF7kKdyGh4frm2++ydL+zTffKCws7IaLAgAAAPIiT18o69+/v4YNG6YLFy6oZcuWkqQ1a9bo+eefz9dfKAMAAAByI0/h9rnnntOff/6pQYMGKT09XZLk4+OjESNGKC4uLl8LBAAAAHIqT7MlZDp79qx++eUX+fr6qkqVKjftr34xWwJQwJgtAXAvZkuAAXKa1/LUc5vJ399f9evXv5FdAAAAAPkmT18oAwAAAIoiwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxvB0dwEAAKCAjbG7u4LsjXG4uwIYqMj33B49elSPPPKISpcuLV9fX9WqVUvffvutc7llWRo1apTKlSsnX19fRUdHa+/evW6sGAAAAO5SpMPt6dOn1bhxYxUvXlxffPGFfv75Z7311lsqVaqUc53x48dr4sSJmjZtmrZs2SI/Pz/FxMTo/PnzbqwcAAAA7lCkhyW8/vrrCg8P1+zZs51tkZGRzj9blqUJEybo3//+tzp16iRJ+uCDDxQSEqLFixerR48ehV4zAAAA3KdI99wuXbpU9erVU/fu3VW2bFnVqVNHM2fOdC4/ePCgEhMTFR0d7Wyz2+1q2LChEhISrrrftLQ0paSkuDwAAABw8yvS4fbAgQOaOnWqqlSpopUrV2rgwIF68sknNXfuXElSYmKiJCkkJMRlu5CQEOey7MTHx8tutzsf4eHhBXcQAAAAKDRFOtxmZGTo7rvv1rhx41SnTh0NGDBA/fv317Rp025ov3FxcXI4HM7HkSNH8qliAAAAuFORDrflypVTjRo1XNqqV6+uw4cPS5JCQ0MlSUlJSS7rJCUlOZdlx9vbWwEBAS4PAAAA3PyKdLht3Lix9uzZ49L266+/KiIiQtLfXy4LDQ3VmjVrnMtTUlK0ZcsWRUVFFWqtAAAAcL8iPVvC008/rXvuuUfjxo3TAw88oK1bt2rGjBmaMWOGJMlms2nYsGF65ZVXVKVKFUVGRmrkyJEKCwtT586d3Vs8AAAACl2RDrf169fXokWLFBcXp5deekmRkZGaMGGCevbs6Vzn+eefV2pqqgYMGKDk5GQ1adJEK1askI+PjxsrBwAAgDvYLMuy3F2Eu6WkpMhut8vhcDD+FigIRfWnPwG4Fz+/i1zIaV4r0mNuAQAAgNwg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGJ7uLgBAPhpjd3cFAAC4FT23AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMMZNFW5fe+012Ww2DRs2zNl2/vx5DR48WKVLl5a/v7+6du2qpKQk9xUJAAAAt7lpwu22bds0ffp03XnnnS7tTz/9tJYtW6aFCxdqw4YNOnbsmLp06eKmKgEAAOBON0W4PXv2rHr27KmZM2eqVKlSznaHw6H3339fb7/9tlq2bKm6detq9uzZ2rRpkzZv3uzGigEAAOAON0W4HTx4sNq1a6fo6GiX9u3bt+vChQsu7dWqVVOFChWUkJBQ2GUCAADAzTzdXcD1LFiwQDt27NC2bduyLEtMTJSXl5cCAwNd2kNCQpSYmHjVfaalpSktLc35PCUlJd/qBQAAgPsU6Z7bI0eO6KmnntK8efPk4+OTb/uNj4+X3W53PsLDw/Nt3wAAAHCfIh1ut2/fruPHj+vuu++Wp6enPD09tWHDBk2cOFGenp4KCQlRenq6kpOTXbZLSkpSaGjoVfcbFxcnh8PhfBw5cqSAjwQAAACFoUgPS2jVqpV+/PFHl7Y+ffqoWrVqGjFihMLDw1W8eHGtWbNGXbt2lSTt2bNHhw8fVlRU1FX36+3tLW9v7wKtHQAAAIWvSIfbkiVL6o477nBp8/PzU+nSpZ3t/fr10/DhwxUUFKSAgAANHTpUUVFRatSokTtKBgAAgBsV6XCbE++88448PDzUtWtXpaWlKSYmRlOmTHF3WQAAAHADm2VZlruLcLeUlBTZ7XY5HA4FBAS4uxwg78bY3V0BAOTcGIe7K8BNJKd5rUh/oQwAAADIDcItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAY3i6uwAAAHCLGmN3dwVZjXG4uwLcIHpuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABj8CMOQF4VxcnHAQC4xdFzCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjFHkw218fLzq16+vkiVLqmzZsurcubP27Nnjss758+c1ePBglS5dWv7+/uratauSkpLcVDEAAADcpciH2w0bNmjw4MHavHmzVq9erQsXLui+++5Tamqqc52nn35ay5Yt08KFC7VhwwYdO3ZMXbp0cWPVAAAAcAebZVmWu4vIjRMnTqhs2bLasGGDmjZtKofDoeDgYM2fP1/dunWTJO3evVvVq1dXQkKCGjVqdN19pqSkyG63y+FwKCAgoKAPAaYYY3d3BQCA/DbG4e4KcBU5zWtFvuf2Sg7H3zddUFCQJGn79u26cOGCoqOjnetUq1ZNFSpUUEJCgltqBAAAgHt4uruA3MjIyNCwYcPUuHFj3XHHHZKkxMREeXl5KTAw0GXdkJAQJSYmZruftLQ0paWlOZ+npKQUWM0AAAAoPDdVz+3gwYP1008/acGCBTe0n/j4eNntducjPDw8nyoEAACAO9004XbIkCFavny51q1bp/LlyzvbQ0NDlZ6eruTkZJf1k5KSFBoamu2+4uLi5HA4nI8jR44UZOkAAAAoJEV+WIJlWRo6dKgWLVqk9evXKzIy0mV53bp1Vbx4ca1Zs0Zdu3aVJO3Zs0eHDx9WVFRUtvv09vaWt7d3gdcOAABuMkX1y8J80S3Hiny4HTx4sObPn68lS5aoZMmSznG0drtdvr6+stvt6tevn4YPH66goCAFBARo6NChioqKytFMCQAAADBHkQ+3U6dOlSQ1b97cpX327Nnq3bu3JOmdd96Rh4eHunbtqrS0NMXExGjKlCmFXCkAAADc7aab57YgMM8t8qSofnQFADAPwxLMnecWAAAAuBrCLQAAAIxR5MfcAnz8DwAAcoqeWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxPN1dAAAAAK5jjN3dFWRvjMPdFWRBzy0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBie7i7gljXG7u4KsjfG4e4KAAAA8oyeWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDGYLQGuiuosDgAAADlgTM/t5MmTVbFiRfn4+Khhw4baunWru0sCAABAITMi3P73v//V8OHDNXr0aO3YsUO1a9dWTEyMjh8/7u7SAAAAUIiMCLdvv/22+vfvrz59+qhGjRqaNm2aSpQooVmzZrm7NAAAABSim37MbXp6urZv3664uDhnm4eHh6Kjo5WQkJDtNmlpaUpLS3M+dzj+/lWulJSUgi3WpQir8F4LAACgIBRidsrMaZZ17Qx104fbkydP6tKlSwoJCXFpDwkJ0e7du7PdJj4+XmPHjs3SHh4eXiA1AgAAGOm1wv8i+pkzZ2S3X/11b/pwmxdxcXEaPny483lGRoZOnTql0qVLy2azubGy3ElJSVF4eLiOHDmigIAAd5dzy+C8Fz7OuXtw3t2D8+4enPfCl9tzblmWzpw5o7CwsGuud9OH2zJlyqhYsWJKSkpyaU9KSlJoaGi223h7e8vb29ulLTAwsKBKLHABAQH8RXQDznvh45y7B+fdPTjv7sF5L3y5OefX6rHNdNN/oczLy0t169bVmjVrnG0ZGRlas2aNoqKi3FgZAAAACttN33MrScOHD1dsbKzq1aunBg0aaMKECUpNTVWfPn3cXRoAAAAKkRHh9sEHH9SJEyc0atQoJSYm6q677tKKFSuyfMnMNN7e3ho9enSWIRYoWJz3wsc5dw/Ou3tw3t2D8174Cuqc26zrzacAAAAA3CRu+jG3AAAAQCbCLQAAAIxBuAUAAIAxCLcAAAAwBuHWEB07dlSFChXk4+OjcuXKqVevXjp27Ji7yzLab7/9pn79+ikyMlK+vr6qXLmyRo8erfT0dHeXZrxXX31V99xzj0qUKHFT/wBLUTd58mRVrFhRPj4+atiwobZu3erukoy2ceNGdejQQWFhYbLZbFq8eLG7SzJefHy86tevr5IlS6ps2bLq3Lmz9uzZ4+6yjDd16lTdeeedzh9viIqK0hdffJFv+yfcGqJFixb6+OOPtWfPHn366afav3+/unXr5u6yjLZ7925lZGRo+vTp2rVrl9555x1NmzZNL774ortLM156erq6d++ugQMHursUY/33v//V8OHDNXr0aO3YsUO1a9dWTEyMjh8/7u7SjJWamqratWtr8uTJ7i7llrFhwwYNHjxYmzdv1urVq3XhwgXdd999Sk1NdXdpRitfvrxee+01bd++Xd9++61atmypTp06adeuXfmyf6YCM9TSpUvVuXNnpaWlqXjx4u4u55bxxhtvaOrUqTpw4IC7S7klzJkzR8OGDVNycrK7SzFOw4YNVb9+fb333nuS/v7lx/DwcA0dOlQvvPCCm6szn81m06JFi9S5c2d3l3JLOXHihMqWLasNGzaoadOm7i7nlhIUFKQ33nhD/fr1u+F90XNroFOnTmnevHm65557CLaFzOFwKCgoyN1lADckPT1d27dvV3R0tLPNw8ND0dHRSkhIcGNlQMFyOBySxPt4Ibp06ZIWLFig1NRURUVF5cs+CbcGGTFihPz8/FS6dGkdPnxYS5YscXdJt5R9+/Zp0qRJevzxx91dCnBDTp48qUuXLmX5lceQkBAlJia6qSqgYGVkZGjYsGFq3Lix7rjjDneXY7wff/xR/v7+8vb21hNPPKFFixapRo0a+bJvwm0R9sILL8hms13zsXv3buf6zz33nL777jutWrVKxYoV06OPPipGneRebs+7JB09elStW7dW9+7d1b9/fzdVfnPLy3kHgPwyePBg/fTTT1qwYIG7S7klVK1aVd9//722bNmigQMHKjY2Vj///HO+7Jsxt0XYiRMn9Oeff15znUqVKsnLyytL+++//67w8HBt2rQp37r5bxW5Pe/Hjh1T8+bN1ahRI82ZM0ceHvyfMS/ycr8z5rZgpKenq0SJEvrkk09cxnzGxsYqOTmZT4UKAWNuC9eQIUO0ZMkSbdy4UZGRke4u55YUHR2typUra/r06Te8L898qAcFJDg4WMHBwXnaNiMjQ5KUlpaWnyXdEnJz3o8ePaoWLVqobt26mj17NsH2BtzI/Y785eXlpbp162rNmjXOcJWRkaE1a9ZoyJAh7i0OyEeWZWno0KFatGiR1q9fT7B1o4yMjHzLLIRbA2zZskXbtm1TkyZNVKpUKe3fv18jR45U5cqV6bUtQEePHlXz5s0VERGhN998UydOnHAuCw0NdWNl5jt8+LBOnTqlw4cP69KlS/r+++8lSbfffrv8/f3dW5whhg8frtjYWNWrV08NGjTQhAkTlJqaqj59+ri7NGOdPXtW+/btcz4/ePCgvv/+ewUFBalChQpurMxcgwcP1vz587VkyRKVLFnSOabcbrfL19fXzdWZKy4uTm3atFGFChV05swZzZ8/X+vXr9fKlSvz5wUs3PR27txptWjRwgoKCrK8vb2tihUrWk888YT1+++/u7s0o82ePduSlO0DBSs2Njbb875u3Tp3l2aUSZMmWRUqVLC8vLysBg0aWJs3b3Z3SUZbt25dtvd1bGysu0sz1tXew2fPnu3u0ozWt29fKyIiwvLy8rKCg4OtVq1aWatWrcq3/TPmFgAAAMZggCAAAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxvj/wScvqoNwoR4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.hist(pm.draw(x, draws=1_000), color=\"C1\", bins=15)\n",
    "ax.set(title=\"Samples from a normal distribution using pymc\", ylabel=\"count\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay! We learned how to sample from a `pymc` distribution!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is going on behind the scenes?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now look into how this is done inside a {class}`~pymc.Model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normal_rv{\"(),()->()\"}.1 [id A] 'z'\n",
      " ├─ RNG(<Generator(PCG64) at 0x7FC6482D1FC0>) [id B]\n",
      " ├─ NoneConst{None} [id C]\n",
      " ├─ [0 0] [id D]\n",
      " └─ [1 2] [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    z = pm.Normal(name=\"z\", mu=np.array([0, 0]), sigma=np.array([1, 2]))\n",
    "\n",
    "pytensor.dprint(z)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are just creating random variables like we saw before, but now registering them in a `pymc` model. To extract the list of random variables we can simply do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[z ~ Normal(<constant>, <constant>)]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.basic_RVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normal_rv{\"(),()->()\"}.1 [id A] 'z'\n",
      " ├─ RNG(<Generator(PCG64) at 0x7FC6482D1FC0>) [id B]\n",
      " ├─ NoneConst{None} [id C]\n",
      " ├─ [0 0] [id D]\n",
      " └─ [1 2] [id E]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(model.basic_RVs[0])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try to sample via {meth}`~pytensor.graph.basic.Variable.eval` as above and it is no surprise that we are getting the same samples at each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0: [-0.78480847  2.20329511]\n",
      "Sample 1: [-0.78480847  2.20329511]\n",
      "Sample 2: [-0.78480847  2.20329511]\n",
      "Sample 3: [-0.78480847  2.20329511]\n",
      "Sample 4: [-0.78480847  2.20329511]\n",
      "Sample 5: [-0.78480847  2.20329511]\n",
      "Sample 6: [-0.78480847  2.20329511]\n",
      "Sample 7: [-0.78480847  2.20329511]\n",
      "Sample 8: [-0.78480847  2.20329511]\n",
      "Sample 9: [-0.78480847  2.20329511]\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(f\"Sample {i}: {z.eval()}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the correct way of sampling is via {func}`~pymc.draw`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample 0: [-1.10363734 -4.33735303]\n",
      "Sample 1: [ 0.69639479 -0.81137315]\n",
      "Sample 2: [ 1.25238284 -0.0119145 ]\n",
      "Sample 3: [ 1.21683809 -3.08878544]\n",
      "Sample 4: [1.63496743 2.58329782]\n",
      "Sample 5: [0.4128748  3.29810689]\n",
      "Sample 6: [1.76074607 3.33727713]\n",
      "Sample 7: [ 0.92855273 -0.14005723]\n",
      "Sample 8: [ 2.04166261 -1.25987621]\n",
      "Sample 9: [-0.24230627 -2.91013171]\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(f\"Sample {i}: {pm.draw(z)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAKqCAYAAACTnV4oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1TklEQVR4nO3deZjVdd34/9eZGWYYVoEARZBFuy3UsguQXL4GSi4/l6xE79JCMnIB16zEu0LtTiy9WkTF5SrsvtLbNZcyF/SG6s41SdNMuklJwlQQHZBltvP5/eHt3E0gQvF+nxl4PK7rXDqHM+/X+5w558xzPnNmplQURREAAJBYVaU3AADAtkF4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngC25xSqRTnn39+pbexUSeccEIMGzas0tsA2KKEJ/APefrpp+Poo4+OoUOHRteuXWPHHXeMj370ozFr1qxKb63DWbx4cZRKpbj00ks3+O/nn39+lEqlWL58+T8159lnn43zzz8/Fi9e/E+tA5CK8AQ220MPPRSjR4+Op556KqZMmRKXX355fP7zn4+qqqr4/ve/X+ntbRWuvfbaWLhw4Wa9z7PPPhsXXHCB8AQ6rJpKbwDofL75zW9G79694/HHH4/tttuu3b+9+uqrldnUVqZLly6V3sJmW7duXdTW1kZVlWMawIZ5dgA225/+9KfYbbfd1ovOiIgBAwa0e3vOnDlxwAEHxIABA6Kuri5GjhwZs2fPXu/9hg0bFocffnjMnz8/Ro8eHfX19bHHHnvE/PnzIyLiJz/5Seyxxx7RtWvXGDVqVPz2t79t9/4nnHBC9OjRI55//vk4+OCDo3v37jFo0KC48MILoyiKd71OS5cujc997nMxcODAqKuri9122y1++MMfrne5WbNmxW677RbdunWLPn36xOjRo+OGG2541/U314Ze43njjTfGqFGjomfPntGrV6/YY4892o4wX3fddTFx4sSIiBg/fnyUSqUolUptt19ExJVXXhm77bZb1NXVxaBBg2Lq1KnxxhtvrDf7iiuuiBEjRkR9fX3stdde8atf/SrGjRsX48aNa7vM/Pnzo1QqxY033hhf/epXY8cdd4xu3brFypUrY8WKFXHOOefEHnvsET169IhevXrFoYceGk899VS7OW+vcfPNN8cFF1wQO+64Y/Ts2TOOPvroaGhoiMbGxjjzzDNjwIAB0aNHj5g8eXI0NjZukdsXqAxHPIHNNnTo0Hj44YfjmWeeid13332jl509e3bstttuceSRR0ZNTU389Kc/jVNPPTXK5XJMnTq13WUXLVoUn/70p+Okk06K448/Pi699NI44ogj4qqrrorzzjsvTj311IiImDlzZhxzzDGxcOHCdkfXWltb45BDDokPf/jD8e1vfzvuvffemDFjRrS0tMSFF174jnt85ZVX4sMf/nCUSqWYNm1a9O/fP+6555448cQTY+XKlXHmmWdGxFvf/j799NPj6KOPjjPOOCPWrVsXv/vd7+LRRx+NT3/60+96u61Zs2aDr+Ncs2bNu77v3Llz41Of+lQceOCB8a1vfSsiIv7whz/Er3/96zjjjDNi//33j9NPPz0uu+yyOO+88+L9739/RETbf88///y44IILYsKECXHKKafEwoULY/bs2fH444/Hr3/967YjrLNnz45p06bF//t//y/OOuusWLx4cRx11FHRp0+fGDx48Hr7+sY3vhG1tbVxzjnnRGNjY9TW1sazzz4bd9xxR0ycODGGDx8er7zySlx99dXxkY98JJ599tkYNGhQuzVmzpwZ9fX1ce6558aiRYti1qxZ0aVLl6iqqorXX389zj///HjkkUfiuuuui+HDh8fXv/71d729gA6qANhM999/f1FdXV1UV1cXe++9d/HlL3+5uO+++4qmpqb1LrtmzZr1zjv44IOLESNGtDtv6NChRUQUDz30UNt59913XxERRX19ffHnP/+57fyrr766iIhi3rx5bedNmjSpiIjitNNOazuvXC4Xhx12WFFbW1ssW7as7fyIKGbMmNH29oknnljssMMOxfLly9vt6V//9V+L3r17t12Hj33sY8Vuu+32LrfO+l544YUiIt719Ld7nDRpUjF06NC2t88444yiV69eRUtLyzvOueWWW9a7XYqiKF599dWitra2OOigg4rW1ta28y+//PIiIoof/vCHRVEURWNjY9GvX79izJgxRXNzc9vlrrvuuiIiio985CNt582bN6+IiGLEiBHrfYzXrVvXbs7bt0FdXV1x4YUXrrfG7rvv3u6+86lPfaoolUrFoYce2m6Nvffeu91tAnQ+vtUObLaPfvSj8fDDD8eRRx4ZTz31VHz729+Ogw8+OHbccce466672l22vr6+7f8bGhpi+fLl8ZGPfCSef/75aGhoaHfZkSNHxt5779329tixYyMi4oADDoiddtppvfOff/759fY2bdq0tv9/+whmU1NTPPDAAxu8LkVRxG233RZHHHFEFEURy5cvbzsdfPDB0dDQEAsWLIiIiO222y7+8pe/xOOPP75Jt9Pf+8IXvhBz585d7/SZz3zmXd93u+22i9WrV8fcuXM3e+4DDzwQTU1NceaZZ7Y7QjxlypTo1atX3H333RER8Zvf/CZee+21mDJlStTU/N83xI477rjo06fPBteeNGlSu49xRERdXV3bnNbW1njttdeiR48eseuuu7bdln/rs5/9bLvXtI4dOzaKoojPfe5z7S43duzYWLJkSbS0tGzmLQB0FL7VDvxDxowZEz/5yU+iqakpnnrqqbj99tvju9/9bhx99NHx5JNPxsiRIyMi4te//nXMmDEjHn744fW+pdzQ0BC9e/due/tv4zIi2v5tyJAhGzz/9ddfb3d+VVVVjBgxot15//Iv/xIR8Y4/6b1s2bJ444034pprrolrrrlmg5d5+wemvvKVr8QDDzwQe+21V+yyyy5x0EEHxac//enYd999N/h+f++9731vTJgwYb3z//u///td3/fUU0+Nm2++OQ499NDYcccd46CDDopjjjkmDjnkkHd93z//+c8REbHrrru2O7+2tjZGjBjR9u9v/3eXXXZpd7mampp3/J2iw4cPX++8crkc3//+9+PKK6+MF154IVpbW9v+rV+/futdfnM+7uVyORoaGja4DtDxOeIJ/FNqa2tjzJgxcdFFF8Xs2bOjubk5brnlloh464eQDjzwwFi+fHl85zvfibvvvjvmzp0bZ511VkS8FSh/q7q6eoMz3un8YhN+aOjdvL2H448/foNHI+fOndsWlu9///tj4cKFceONN8Z+++0Xt912W+y3334xY8aMf3of72bAgAHx5JNPxl133RVHHnlkzJs3Lw499NCYNGlS8tkb8/dHOyMiLrroojj77LNj//33jx//+Mdx3333xdy5c2O33XZb72MeUZmPO1AZjngCW8zo0aMjIuKvf/1rRET89Kc/jcbGxrjrrrvaHdWaN29ekvnlcjmef/75tqOcERF//OMfIyLe8Yhd//79o2fPntHa2rrBo5F/r3v37nHsscfGscceG01NTfGJT3wivvnNb8b06dOja9euW+R6vJPa2to44ogj4ogjjohyuRynnnpqXH311fG1r30tdtlllyiVSht8v6FDh0ZExMKFC9sdEW5qaooXXnih7Xq/fblFixbF+PHj2y7X0tISixcvjg984AObtM9bb701xo8fHz/4wQ/anf/GG2/Ee97znk2/wsBWxxFPYLPNmzdvg0edfv7zn0fE/31L9+0jVn972YaGhpgzZ06yvV1++eVt/18URVx++eXRpUuXOPDAAzd4+erq6vjkJz8Zt912WzzzzDPr/fuyZcva/v+1115r92+1tbUxcuTIKIoimpubt9A12LC/n11VVdUWgm//iqHu3btHRKz3K5ImTJgQtbW1cdlll7X7WPzgBz+IhoaGOOywwyLirS8c+vXrF9dee22711Fef/31672sYWOqq6vXu3/ccsstsXTp0k1eA9g6OeIJbLbTTjst1qxZEx//+Mfjfe97XzQ1NcVDDz0UN910UwwbNiwmT54cEREHHXRQ21G6k046Kd5888249tprY8CAAW1HRbekrl27xr333huTJk2KsWPHxj333BN33313nHfeedG/f/93fL+LL7445s2bF2PHjo0pU6bEyJEjY8WKFbFgwYJ44IEHYsWKFW3XZ/vtt4999903Bg4cGH/4wx/i8ssvj8MOOyx69uy5xa/P3/r85z8fK1asiAMOOCAGDx4cf/7zn2PWrFmx5557tv3KpD333DOqq6vjW9/6VjQ0NERdXV3b71CdPn16XHDBBXHIIYfEkUceGQsXLowrr7wyxowZE8cff3xEvBXS559/fpx22mlxwAEHxDHHHBOLFy+O6667Lnbeeed3PKL69w4//PC48MILY/LkybHPPvvE008/Hddff/16r78Ftj2OeAKb7dJLL43x48fHz3/+8zj77LPj7LPPjsceeyxOPfXUePTRR9t+sfyuu+4at956a5RKpTjnnHPiqquuii984QtxxhlnJNlXdXV13HvvvfHyyy/Hl770pXj88cdjxowZ8Y1vfGOj7zdw4MB47LHHYvLkyfGTn/wkpk2bFt///vdjxYoVbb8zMyLa4vk73/lOTJ06Ne644444/fTT48c//nGS6/O3jj/++OjatWtceeWVceqpp8aPfvSjOPbYY+Oee+5p+wny7bffPq666qp49dVX48QTT4xPfepT8eyzz0bEW7/H8/LLL48XX3wxzjrrrLj55pvjC1/4Qtx///3tfqJ82rRpcdlll8WLL74Y55xzTvzqV7+Ku+66K7bbbrtNfinBeeedF1/84hfjvvvuizPOOCMWLFgQd99993o/LARse0qFV2kDW4ETTjghbr311njzzTcrvZWtTrlcjv79+8cnPvGJuPbaayu9HaATc8QTgDbr1q1b7/WZ//Ef/xErVqxo9yczAf4RXuMJQJtHHnkkzjrrrJg4cWL069cvFixYED/4wQ9i9913b/tb8AD/KOEJQJthw4bFkCFD4rLLLosVK1ZE375947Of/WxcfPHFUVtbW+ntAZ2c13gCAJCF13gCAJCF8AQAIIsO/RrPcrkcL730UvTs2XOTf3ExAAD5FEURq1atikGDBrX9XuF30qHD86WXXvILhwEAOoElS5bE4MGDN3qZDh2eb/8Juv3i/4ua6PIul4ZtTCnxK2WKctr1I9Jfh4gs16NUk/75qWhJ+7fgtxpbyX3K43sT5bgevKuWaI7/jp9v0p8O7tDh+fa312uiS9SUhCe0k/xJfSv5xJThepQyPD8VXm20abaS+5TH96YSnh3C//5+pE15WaQfLgIAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZFFT6Q3AVqmU4Wu6opx0+ara2qTrR0QU5SL9jJa0t1NERKmqlHxG8lsqw322VF2dfEbR0px8Bpso8XMUnZMjngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsnDc+nSpXH88cdHv379or6+PvbYY4/4zW9+k3osAAAdTNJfp/T666/HvvvuG+PHj4977rkn+vfvH//zP/8Tffr0STkWAIAOKGl4futb34ohQ4bEnDlz2s4bPnx4ypEAAHRQSb/Vftddd8Xo0aNj4sSJMWDAgPjQhz4U1157bcqRAAB0UEnD8/nnn4/Zs2fHe9/73rjvvvvilFNOidNPPz1+9KMfbfDyjY2NsXLlynYnAAC2Dkm/1V4ul2P06NFx0UUXRUTEhz70oXjmmWfiqquuikmTJq13+ZkzZ8YFF1yQcksAAFRI0iOeO+ywQ4wcObLdee9///vjxRdf3ODlp0+fHg0NDW2nJUuWpNweAAAZJT3iue+++8bChQvbnffHP/4xhg4dusHL19XVRV1dXcotAQBQIUmPeJ511lnxyCOPxEUXXRSLFi2KG264Ia655pqYOnVqyrEAAHRAScNzzJgxcfvtt8d//ud/xu677x7f+MY34nvf+14cd9xxKccCANABJf1We0TE4YcfHocffnjqMQAAdHD+VjsAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyS/wJ52BZVden8D61yc0vyGVVd65LPKGX4WJTXNSafUb3ddmkHtLamXT8iipb096lSVSn5jHJTU/IZUZTTz4AKcMQTAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIIuaSm+AzqOqtjb5jHJTU/IZOWS5HqW0XzfWvKdv0vUjIspvrk4+I4fqPr3TD2luST8jsaoe3ZPPaH29IfmMUk2X5DOKluak6+d4Pi/KRfoZra3JZ0RRTj9jG+KIJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCxqKr0BtpxSTZek65ebmpKuHxERpfRfC5Wqq5PPKFqak8+oqu+adP1iXWPS9SMiqvr1TT6jeHN18hmlPtslnxE1ie+3K99Mu34mVT26V3oLW0ZVKenyOR4XRWtr8hlRlNPPYItyxBMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAssoXnxRdfHKVSKc4888xcIwEA6ECyhOfjjz8eV199dXzgAx/IMQ4AgA4oeXi++eabcdxxx8W1114bffr0ST0OAIAOKnl4Tp06NQ477LCYMGHCu162sbExVq5c2e4EAMDWoSbl4jfeeGMsWLAgHn/88U26/MyZM+OCCy5IuSUAACok2RHPJUuWxBlnnBHXX399dO3adZPeZ/r06dHQ0NB2WrJkSartAQCQWakoiiLFwnfccUd8/OMfj+rq6rbzWltbo1QqRVVVVTQ2Nrb7tw1ZuXJl9O7dO8bFx6Km1CXFNrcupbSvnCi9y8ersyhVldIPqUn6zYSIiKjarnfaAWmeGtqrq0s+onVg4tspIqpWrks+o7l/9+QzUqv9n79mGJLhc0VzS/oZXdI+h5SXvZZ0/YiIoqk5+Ywcz+flpqbkMzq7lqI55sed0dDQEL169droZZPdsw888MB4+umn2503efLkeN/73hdf+cpX3jU6AQDYuiQLz549e8buu+/e7rzu3btHv3791jsfAICtn79cBABAFulfiPY35s+fn3McAAAdiCOeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIIuaSm+AzqNobU0+o6pLhrtkdXXyEVXduyWfUR7YJ+n6VStWJV0/IqJx+HuSz8hh7bD0H+9S4odfl7XltAMiYu2+w5LPqGtoST6jy+uNyWdUv7Q86fql2tqk60dElNeuSz4jIv3zeZQyHKMr0j/+OgpHPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGRRU+kN0HlUda1LPqNUk+EuWV2dfkavnulnJNY8pF/yGY19uySfsWpw+q+va1cVyWes7V9Kun51Y/rbqdur6W+nmtUtyWe0dkv/PFXdrT7p+qXm5qTrR+T5nFE0p/94lzJ8zihaUw8oJx6w6RzxBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMiiptIbYAsqykmXL9Wkv7sUjY3JZ1Rt1zv5jGhN+7GIiGjq1y3t+r3Tf7xXDk3/tW9L2pspIiIadm9JPqPmjbQfjx67vp50/YiIN57om3xGdWPX5DNKrUXyGdWNPdKuv2Zt0vUj8jzX5vic0fp6Q/IZ6T9/d0m7fhERm/g06IgnAABZCE8AALIQngAAZCE8AQDIIml4zpw5M8aMGRM9e/aMAQMGxFFHHRULFy5MORIAgA4qaXj+4he/iKlTp8YjjzwSc+fOjebm5jjooINi9erVKccCANABJf39HPfee2+7t6+77roYMGBAPPHEE7H//vunHA0AQAeT9TWeDQ1v/S6svn3T/y43AAA6lmy/QL5cLseZZ54Z++67b+y+++4bvExjY2M0/s0vg125cmWu7QEAkFi2I55Tp06NZ555Jm688cZ3vMzMmTOjd+/ebachQ4bk2h4AAIllCc9p06bFz372s5g3b14MHjz4HS83ffr0aGhoaDstWbIkx/YAAMgg6bfai6KI0047LW6//faYP39+DB8+fKOXr6uri7q6upRbAgCgQpKG59SpU+OGG26IO++8M3r27Bkvv/xyRET07t076uvrU44GAKCDSfqt9tmzZ0dDQ0OMGzcudthhh7bTTTfdlHIsAAAdUPJvtQMAQIS/1Q4AQCbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFkn/chF/o9T5G79oaan0FraM2trkI4pe3ZLPKNelvU+11JeSrh8R8eawcvIZRU36v6CW4+Hd2j3tbTWo18qk60dE/GFg7+QzlldVJ5+xw0PNyWe01qW9HtX1XZOuHxERK1clH1Gsa0w+I4r0z1OpFa2tadcvNn39zl9DAAB0CsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIIuaSm9gW1Gqrk4/o0vaD2epri7p+hERpdouyWdE1/TXo1yb/qHV2Cvt142rty8lXT8iom7Qm+ln1DYnn/Hb0Tcln/GpFw5Iuv4bTV2Trh8REV1bk49o7pn+eMqK96d/nnrPU+Wk6xfd03+8S2vXJZ8RRZF8RFVtbfIZ5eaWpOunbpBSUY7YxKvgiCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFjWV3sC2olRVSj6jvK4x6fql5pak60dEVPfpnXxGrF6TfETrjtsln9HUK+3XjdVp704REVFUl5PP2L7nquQzHlxbnXxGbVXax9++/Z5Pun5ExMrGrslnLHt1YPIZUaQfsWaH2qTr176W/nkwatNeh4iIUk36jGlN/Lk1IiKK1M+F6Z+jNpUjngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyyBKeV1xxRQwbNiy6du0aY8eOjcceeyzHWAAAOpDk4XnTTTfF2WefHTNmzIgFCxbEBz/4wTj44IPj1VdfTT0aAIAOJHl4fuc734kpU6bE5MmTY+TIkXHVVVdFt27d4oc//GHq0QAAdCBJw7OpqSmeeOKJmDBhwv8NrKqKCRMmxMMPP7ze5RsbG2PlypXtTgAAbB2Shufy5cujtbU1Bg4c2O78gQMHxssvv7ze5WfOnBm9e/duOw0ZMiTl9gAAyKhD/VT79OnTo6Ghoe20ZMmSSm8JAIAtpCbl4u95z3uiuro6XnnllXbnv/LKK7H99tuvd/m6urqoq6tLuSUAACok6RHP2traGDVqVDz44INt55XL5XjwwQdj7733TjkaAIAOJukRz4iIs88+OyZNmhSjR4+OvfbaK773ve/F6tWrY/LkyZu+SKnqrVMqRTnd2hmVqquTrl9V3zXp+hERkeGId9Gze/oZVaXkM2pXFUnXX/ahpMtHRERrQ33yGW92W5t8xsLGHZLP+GCvtC89en5t/6TrR0Q0tiT/lBPVjclHRGuGb8x1+2tT2gHNrWnXj4hoSnwdIqJoTP8Br+qS/n5bbm5Jun7RmvbjXRSbvn7yW/PYY4+NZcuWxde//vV4+eWXY88994x77713vR84AgBg65Y+4yNi2rRpMW3atByjAADooDrUT7UDALD1Ep4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWNZXewCYpyhFRrvQu/inlpqbkM0o1XZKuX7S0JF0/IqKUfEJE0TX93b6pd/oZNWvTPiaqmquTrh8R0bo6/YxlDT2Sz1g3OO1jLyKisZx2Ro/qxqTrR0T0qluXfMaa1uQjosubGWY0pL2tmvt3T7p+RESXv76afEZRLpLPKFXl+MyUWJG4oTZjfUc8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALKoqfQGOoTS1tHfpapS2gGtrWnXj4jyiteTz6jqWpd8Rt2K9DPWDkg7o/+CctL1IyIadq5OPqPbsLXJZ1z+yIHJZ/Tuvyrp+n3q099Oi5/bIfmMuvR3qeiyukg/pEg7o8uy1UnXj4go1XdNP6OxMfmM8rr0M5JL3jlVEZt4l906igsAgA5PeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkUVPpDXQIRTn9jFL6xi/KRdL1S13S311K1dXJZxQrXk8+o7a5OfmMlm7bJ12/XFtKun5ERPW69I+LN57qn3xGqUf655A3V/RJuv7axr5J14+IqF+XfER0WZV+Rs8Xm9IPKSV+/P3lr2nXj4jo2SP9jNbW5COyfF5KfT1Sd85mrO+IJwAAWQhPAACySBaeixcvjhNPPDGGDx8e9fX1sfPOO8eMGTOiqSnDtygAAOhwkr1o77nnnotyuRxXX3117LLLLvHMM8/ElClTYvXq1XHppZemGgsAQAeVLDwPOeSQOOSQQ9reHjFiRCxcuDBmz54tPAEAtkFZf6q9oaEh+vZ955+cbGxsjMbGxra3V65cmWNbAABkkO2HixYtWhSzZs2Kk0466R0vM3PmzOjdu3fbaciQIbm2BwBAYpsdnueee26USqWNnp577rl277N06dI45JBDYuLEiTFlypR3XHv69OnR0NDQdlqyZMnmXyMAADqkzf5W+xe/+MU44YQTNnqZESNGtP3/Sy+9FOPHj4999tknrrnmmo2+X11dXdTV1W3ulgAA6AQ2Ozz79+8f/ftv2l8CWbp0aYwfPz5GjRoVc+bMiaoqvzYUAGBbleyHi5YuXRrjxo2LoUOHxqWXXhrLli1r+7ftt0/7p/4AAOh4koXn3LlzY9GiRbFo0aIYPHhwu38rirR/UxwAgI4n2fe+TzjhhCiKYoMnAAC2PV50CQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWdRUegMdQilDfxfl9CNa065f6pL+7lK0Jr4SERGNjclHlHr2SD6jdmVz0vWLDA+LASuaks9YObw++YyIUvIJbw5OO6Pr8qTLvzXj9fTPgz2WrEk+o2pt2sdeRERp+RtJ1y9yPJ+vWZt8Rrm5JfmMHJ+/tyWOeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkUVPpDXQIRbnSO9gyEl+P8rrGpOtHRFR1rUs+I6qr08+oSf/Qqnnp9aTrF726JV0/IqKlV/qPd1VLkXxG/fKW5DPqGtIeJ6h9oznp+hERXd5Yl3xGNLcmH1HUd0k+o5R4/WJt+o9Fjs8ZpQzP50XLVtIIHYQjngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALKoqfQGthkljb8pyusak8+o7tE9+YzWP/8l+YxSbZek61fVbJ90/YiILg1vJp/R5+W0t1Mu5Z5dk65fampJun5ERJTTjyitWp1+xpr0z+fFmrVpB7S2pl1/K1JVX598Rnlt4o93B6KGAADIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyyhGdjY2PsueeeUSqV4sknn8wxEgCADiZLeH75y1+OQYMG5RgFAEAHlTw877nnnrj//vvj0ksvTT0KAIAOrCbl4q+88kpMmTIl7rjjjujWrVvKUQAAdHDJwrMoijjhhBPi5JNPjtGjR8fixYvf9X0aGxujsbGx7e2VK1em2h4AAJlt9rfazz333CiVShs9PffcczFr1qxYtWpVTJ8+fZPXnjlzZvTu3bvtNGTIkM3dHgAAHVSpKIpic95h2bJl8dprr230MiNGjIhjjjkmfvrTn0apVGo7v7W1Naqrq+O4446LH/3oR+u934aOeA4ZMiTGxceiptRlc7bZ8ZT85qqOorpH9+QzymvXJZ9Rqk37mKgatH3S9SMiork5/Ywunfy543+Ve3ZNun6pqSXp+hERUU4/orRqdfoh1emfz4uGVWnXX7Mm6foREeXm9PepUnV1+hldkr4qMSIiymvXJp+RUkvRHPPjzmhoaIhevXpt9LKbfWv2798/+vfv/66Xu+yyy+Lf//3f295+6aWX4uCDD46bbropxo4du8H3qauri7q6us3dEgAAnUCyjN9pp53avd2jR4+IiNh5551j8ODBqcYCANBB+f4vAABZpH/hwv8aNmxYbObLSfPJ8frLIsOLm9gkrW+mfw1YVX3a1+NFRBRNaV8f2bp4SdL1I9K/TjUiomqHgclnxJr0r8+qSv162MT3p2xaW5OPKDK8hrtoSfv6yKKc4fNxhs97RfoPdxQtW8ljo4NwxBMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWNZXeQIdQlCu9gy2jlPjriK3ldspwPYqm5uQzSlWlpOsX5SLp+hERpdra5DPKf30l+YwsH+/aLonXT/+xKNasST4jqqvTz8igvHZtpbfQOWwtn5e2IY54AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyKKm0htgCyrKld4B/6tobc0wI+36VV3r0g6IiNaGlclnlKqrk8/Iobx2XdL1S03NSdePiChVlZLPKJpb0s/I8Piuqq1Nun45w+2Uhc97nY4jngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZ1FR6A5BdKcPXW0U5/YzE16O8dm3S9SMiSjVdks8oWluTzyhVVyefkfo+VaS/maJoyfC4yCHDc0i5qSntgBzPg7AB7nkAAGSRNDzvvvvuGDt2bNTX10efPn3iqKOOSjkOAIAOLNm32m+77baYMmVKXHTRRXHAAQdES0tLPPPMM6nGAQDQwSUJz5aWljjjjDPikksuiRNPPLHt/JEjR6YYBwBAJ5DkW+0LFiyIpUuXRlVVVXzoQx+KHXbYIQ499NB3PeLZ2NgYK1eubHcCAGDrkCQ8n3/++YiIOP/88+OrX/1q/OxnP4s+ffrEuHHjYsWKFe/4fjNnzozevXu3nYYMGZJiewAAVMBmhee5554bpVJpo6fnnnsuyuW3fmXGv/3bv8UnP/nJGDVqVMyZMydKpVLccsst77j+9OnTo6Ghoe20ZMmSf+7aAQDQYWzWazy/+MUvxgknnLDRy4wYMSL++te/RkT713TW1dXFiBEj4sUXX3zH962rq4u6urrN2RIAAJ3EZoVn//79o3///u96uVGjRkVdXV0sXLgw9ttvv4iIaG5ujsWLF8fQoUP/sZ0CANCpJfmp9l69esXJJ58cM2bMiCFDhsTQoUPjkksuiYiIiRMnphgJAEAHl+z3eF5yySVRU1MTn/nMZ2Lt2rUxduzY+K//+q/o06dPqpEAAHRgpaIoikpv4p2sXLkyevfuHePiY1FTSv83ndlG+FvtmybDdfC32jdd0dKcdsDW8rjIYWu4rbaWv9W+tdynOrmWojnmx53R0NAQvXr12uhlt5J7HgAAHZ3wBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALGoqvQHIrihXegdbxlZwPUpVpeQzitbkI6JoaU4/hI5jK3jsbRXXgU7JEU8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZ1FR6A7BVKmX4mq4op5+RWLmpqdJb4G1bwf0pIjz2oINzxBMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALIQnAABZCE8AALIQngAAZCE8AQDIQngCAJCF8AQAIAvhCQBAFsITAIAshCcAAFkITwAAshCeAABkITwBAMhCeAIAkEVNpTewMUVRRERESzRHFBXeDGyWDF/TFeX0M6DT8diD3FqiOSL+r9s2pkOH56pVqyIi4r/j5xXeCWwmXyhBZXjsQcWsWrUqevfuvdHLlIpNydMKKZfL8dJLL0XPnj2jVCpVejudzsqVK2PIkCGxZMmS6NWrV6W302m5HbcMt+OW47bcMtyOW4bbccvozLdjURSxatWqGDRoUFRVbfy7Dh36iGdVVVUMHjy40tvo9Hr16tXp7sQdkdtxy3A7bjluyy3D7bhluB23jM56O77bkc63+eEiAACyEJ4AAGQhPLdidXV1MWPGjKirq6v0Vjo1t+OW4XbcctyWW4bbcctwO24Z28rt2KF/uAgAgK2HI54AAGQhPAEAyEJ4AgCQhfAEACAL4bmNOPLII2OnnXaKrl27xg477BCf+cxn4qWXXqr0tjqVxYsXx4knnhjDhw+P+vr62HnnnWPGjBnR1NRU6a11St/85jdjn332iW7dusV2221X6e10GldccUUMGzYsunbtGmPHjo3HHnus0lvqdH75y1/GEUccEYMGDYpSqRR33HFHpbfUKc2cOTPGjBkTPXv2jAEDBsRRRx0VCxcurPS2Op3Zs2fHBz7wgbZfHL/33nvHPffcU+ltJSM8txHjx4+Pm2++ORYuXBi33XZb/OlPf4qjjz660tvqVJ577rkol8tx9dVXx+9///v47ne/G1dddVWcd955ld5ap9TU1BQTJ06MU045pdJb6TRuuummOPvss2PGjBmxYMGC+OAHPxgHH3xwvPrqq5XeWqeyevXq+OAHPxhXXHFFpbfSqf3iF7+IqVOnxiOPPBJz586N5ubmOOigg2L16tWV3lqnMnjw4Lj44ovjiSeeiN/85jdxwAEHxMc+9rH4/e9/X+mtJeHXKW2j7rrrrjjqqKOisbExunTpUuntdFqXXHJJzJ49O55//vlKb6XTuu666+LMM8+MN954o9Jb6fDGjh0bY8aMicsvvzwiIsrlcgwZMiROO+20OPfccyu8u86pVCrF7bffHkcddVSlt9LpLVu2LAYMGBC/+MUvYv/996/0djq1vn37xiWXXBInnnhipbeyxTniuQ1asWJFXH/99bHPPvuIzn9SQ0ND9O3bt9LbYBvQ1NQUTzzxREyYMKHtvKqqqpgwYUI8/PDDFdwZvKWhoSEiwnPiP6G1tTVuvPHGWL16dey9996V3k4SwnMb8pWvfCW6d+8e/fr1ixdffDHuvPPOSm+pU1u0aFHMmjUrTjrppEpvhW3A8uXLo7W1NQYOHNju/IEDB8bLL79coV3BW8rlcpx55pmx7777xu67717p7XQ6Tz/9dPTo0SPq6uri5JNPjttvvz1GjhxZ6W0lITw7sXPPPTdKpdJGT88991zb5b/0pS/Fb3/727j//vujuro6PvvZz4ZXWmz+7RgRsXTp0jjkkENi4sSJMWXKlArtvOP5R25LoPObOnVqPPPMM3HjjTdWeiud0q677hpPPvlkPProo3HKKafEpEmT4tlnn630tpLwGs9ObNmyZfHaa69t9DIjRoyI2tra9c7/y1/+EkOGDImHHnpoqz2cv6k293Z86aWXYty4cfHhD384rrvuuqiq8vXb2/6R+6TXeG6apqam6NatW9x6663tXo84adKkeOONN3wH4x/kNZ7/vGnTpsWdd94Zv/zlL2P48OGV3s5WYcKECbHzzjvH1VdfXemtbHE1ld4A/7j+/ftH//79/6H3LZfLERHR2Ni4JbfUKW3O7bh06dIYP358jBo1KubMmSM6/84/c59k42pra2PUqFHx4IMPtkVSuVyOBx98MKZNm1bZzbFNKooiTjvttLj99ttj/vz5onMLKpfLW+3nZ+G5DXj00Ufj8ccfj/322y/69OkTf/rTn+JrX/ta7Lzzztv80c7NsXTp0hg3blwMHTo0Lr300li2bFnbv22//fYV3Fnn9OKLL8aKFSvixRdfjNbW1njyyScjImKXXXaJHj16VHZzHdTZZ58dkyZNitGjR8dee+0V3/ve92L16tUxefLkSm+tU3nzzTdj0aJFbW+/8MIL8eSTT0bfvn1jp512quDOOpepU6fGDTfcEHfeeWf07Nmz7bXGvXv3jvr6+grvrvOYPn16HHroobHTTjvFqlWr4oYbboj58+fHfffdV+mtpVGw1fvd735XjB8/vujbt29RV1dXDBs2rDj55JOLv/zlL5XeWqcyZ86cIiI2eGLzTZo0aYO35bx58yq9tQ5t1qxZxU477VTU1tYWe+21V/HII49Uekudzrx58zZ435s0aVKlt9apvNPz4Zw5cyq9tU7lc5/7XDF06NCitra26N+/f3HggQcW999/f6W3lYzXeAIAkIUXqAEAkIXwBAAgC+EJAEAWwhMAgCyEJwAAWQhPAACyEJ4AAGQhPAEAyEJ4AgCQhfAEACAL4QkAQBbCEwCALP5/oFHUtHwYyCQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "z_draws = pm.draw(vars=z, draws=10_000)\n",
    "ax.hist2d(x=z_draws[:, 0], y=z_draws[:, 1], bins=25)\n",
    "ax.set(title=\"Samples Histogram\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Enough with Random Variables, I want to see some (log)probabilities!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall we have defined the following model above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\n",
       "            \\begin{array}{rcl}\n",
       "            \\text{z} &\\sim & \\operatorname{Normal}(\\text{<constant>},~\\text{<constant>})\n",
       "            \\end{array}\n",
       "            $$"
      ],
      "text/plain": [
       "z ~ Normal(<constant>, <constant>)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pymc` is able to convert `RandomVariable`s to their respective probability functions. One simple way is to use {func}`~pymc.logp`, which takes as first input a RandomVariable, and as second input the value at which the logp is evaluated (we will discuss this in more detail later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_value = pt.vector(name=\"z\")\n",
    "z_logp = pm.logp(rv=z, value=z_value)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`z_logp` contains the PyTensor graph that represents the log-probability of the normal random variable `z`, evaluated at `z_value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check{sigma > 0} [id A] 'z_logprob'\n",
      " ├─ Sub [id B]\n",
      " │  ├─ Sub [id C]\n",
      " │  │  ├─ Mul [id D]\n",
      " │  │  │  ├─ ExpandDims{axis=0} [id E]\n",
      " │  │  │  │  └─ -0.5 [id F]\n",
      " │  │  │  └─ Pow [id G]\n",
      " │  │  │     ├─ True_div [id H]\n",
      " │  │  │     │  ├─ Sub [id I]\n",
      " │  │  │     │  │  ├─ z [id J]\n",
      " │  │  │     │  │  └─ [0 0] [id K]\n",
      " │  │  │     │  └─ [1 2] [id L]\n",
      " │  │  │     └─ ExpandDims{axis=0} [id M]\n",
      " │  │  │        └─ 2 [id N]\n",
      " │  │  └─ ExpandDims{axis=0} [id O]\n",
      " │  │     └─ Log [id P]\n",
      " │  │        └─ Sqrt [id Q]\n",
      " │  │           └─ 6.283185307179586 [id R]\n",
      " │  └─ Log [id S]\n",
      " │     └─ [1 2] [id L]\n",
      " └─ All{axes=None} [id T]\n",
      "    └─ MakeVector{dtype='bool'} [id U]\n",
      "       └─ All{axes=None} [id V]\n",
      "          └─ Gt [id W]\n",
      "             ├─ [1 2] [id L]\n",
      "             └─ ExpandDims{axis=0} [id X]\n",
      "                └─ 0 [id Y]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(z_logp)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "There is also a handy `pymc` function to compute the log cumulative probability of a random variable {func}`~pymc.logcdf`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that, as explained at the beginning, there has been no computation yet. The actual computation is performed after compiling and passing the input. For illustration purposes alone, we will again use the handy {meth}`~pytensor.graph.basic.Variable.eval` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.91893853, -1.61208571])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z_logp.eval({z_value: [0, 0]})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is nothing else than evaluating the log probability of a normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.91893853, -1.61208571])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.norm.logpdf(x=np.array([0, 0]), loc=np.array([0, 0]), scale=np.array([1, 2]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pymc` models provide some helpful routines to facilitating the conversion of `RandomVariable`s to probability functions. {meth}`~pymc.Model.logp`, for instance can be used to extract the joint probability of all variables in the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check{sigma > 0} [id A] 'z_logprob'\n",
      " ├─ Sub [id B]\n",
      " │  ├─ Sub [id C]\n",
      " │  │  ├─ Mul [id D]\n",
      " │  │  │  ├─ ExpandDims{axis=0} [id E]\n",
      " │  │  │  │  └─ -0.5 [id F]\n",
      " │  │  │  └─ Pow [id G]\n",
      " │  │  │     ├─ True_div [id H]\n",
      " │  │  │     │  ├─ Sub [id I]\n",
      " │  │  │     │  │  ├─ z [id J]\n",
      " │  │  │     │  │  └─ [0 0] [id K]\n",
      " │  │  │     │  └─ [1 2] [id L]\n",
      " │  │  │     └─ ExpandDims{axis=0} [id M]\n",
      " │  │  │        └─ 2 [id N]\n",
      " │  │  └─ ExpandDims{axis=0} [id O]\n",
      " │  │     └─ Log [id P]\n",
      " │  │        └─ Sqrt [id Q]\n",
      " │  │           └─ 6.283185307179586 [id R]\n",
      " │  └─ Log [id S]\n",
      " │     └─ [1 2] [id L]\n",
      " └─ All{axes=None} [id T]\n",
      "    └─ MakeVector{dtype='bool'} [id U]\n",
      "       └─ All{axes=None} [id V]\n",
      "          └─ Gt [id W]\n",
      "             ├─ [1 2] [id L]\n",
      "             └─ ExpandDims{axis=0} [id X]\n",
      "                └─ 0 [id Y]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<ipykernel.iostream.OutStream at 0x7fc6cd274a90>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pytensor.dprint(model.logp(sum=False))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we only have one variable, this function is equivalent to what we obtained by manually calling `pm.logp` before. We can also use a helper {meth}`~pymc.Model.compile_logp` to return an already compiled PyTensor function of the model logp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "logp_function = model.compile_logp(sum=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function expects a \"point\" dictionary as input. We could create it ourselves, but just to illustrate another useful {class}`~pymc.Model` method, let's call {meth}`~pymc.Model.initial_point`, which returns the point that most samplers use when deciding where to start sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'z': array([0., 0.])}"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point = model.initial_point()\n",
    "point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([-0.91893853, -1.61208571])]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logp_function(point)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What are value variables and why are they important?\n",
    "\n",
    "As we have seen above, a logp graph does not have random variables. Instead it's defined in terms of input (value) variables. When we want to sample, each random variable (RV) is replaced by a logp function evaluated at the respective input (value) variable. Let's see how this works through some examples. RV and value variables can be observed in these [`scipy`](https://github.com/scipy/scipy) operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<scipy.stats._distn_infrastructure.rv_continuous_frozen at 0x7fc64808be10>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv = scipy.stats.norm(0, 1)\n",
    "\n",
    "# Equivalent to rv = pm.Normal(\"rv\", 0, 1)\n",
    "scipy.stats.norm(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.00721939, -0.60656542, -0.28202337])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Equivalent to rv_draw = pm.draw(rv, 3)\n",
    "rv.rvs(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.7001885332046727"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Equivalent to rv_logp = pm.logp(rv, 1.25)\n",
    "rv.logpdf(1.25)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's look at how these value variables behave in a slightly more complex model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model_2:\n",
    "    mu = pm.Normal(name=\"mu\", mu=0, sigma=2)\n",
    "    sigma = pm.HalfNormal(name=\"sigma\", sigma=3)\n",
    "    x = pm.Normal(name=\"x\", mu=mu, sigma=sigma)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each model RV is related to a \"value variable\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{mu ~ Normal(0, 2): mu,\n",
       " sigma ~ HalfNormal(0, 3): sigma_log__,\n",
       " x ~ Normal(mu, sigma): x}"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_2.rvs_to_values"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that for sigma the associated value is in the *log* scale as in practice we require unbounded values for NUTS sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[mu, sigma_log__, x]"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_2.value_vars"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know how to extract the model variables, we can compute the element-wise log-probability of the model for specific values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -1.61208572, -11.32440366,   9.08106147])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# extract values as pytensor.tensor.var.TensorVariable\n",
    "mu_value = model_2.rvs_to_values[mu]\n",
    "sigma_log_value = model_2.rvs_to_values[sigma]\n",
    "x_value = model_2.rvs_to_values[x]\n",
    "# element-wise log-probability of the model (we do not take te sum)\n",
    "logp_graph = pt.stack(model_2.logp(sum=False))\n",
    "# evaluate by passing concrete values\n",
    "logp_graph.eval({mu_value: 0, sigma_log_value: -10, x_value: 0})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This equivalent to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "mu_value -> -1.612085713764618\n",
      "sigma_log_value -> -11.324403641427345\n",
      "x_value -> 9.081061466795328\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"\"\"\n",
    "mu_value -> {scipy.stats.norm.logpdf(x=0, loc=0, scale=2)}\n",
    "sigma_log_value -> {- 10 + scipy.stats.halfnorm.logpdf(x=np.exp(-10), loc=0, scale=3)}\n",
    "x_value -> {scipy.stats.norm.logpdf(x=0, loc=0, scale=np.exp(-10))}\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{Note}\n",
    "For `sigma_log_value` we add the $-10$ term for the `scipy` and `pytensor` to match because of the jacobian.\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we already saw, we can also use the method {meth}`~pymc.Model.compile_logp` to obtain a compiled pytensor function of the model logp, which takes a dictionary of `{value variable name : value}` as inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array(-1.61208572), array(-11.32440366), array(9.08106147)]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_2.compile_logp(sum=False)({\"mu\": 0, \"sigma_log__\": -10, \"x\": 0})"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The {class}`~pymc.Model` class also has methods to extract the gradient ({meth}`~pymc.Model.dlogp`) and the hessian ({meth}`~pymc.Model.d2logp`) of the logp."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to go deeper into the internals of `pytensor` RandomVariables and `pymc` distributions please take a look into the [distribution developer guide](implementing-a-distribution)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Tue Jun 25 2024\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.11.8\n",
      "IPython version      : 8.22.2\n",
      "\n",
      "pytensor: 2.20.0+3.g66439d283.dirty\n",
      "\n",
      "pytensor  : 2.20.0+3.g66439d283.dirty\n",
      "pymc      : 5.15.0+1.g58927d608\n",
      "scipy     : 1.12.0\n",
      "numpy     : 1.26.4\n",
      "matplotlib: 3.8.3\n",
      "\n",
      "Watermark: 2.4.3\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w -p pytensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "dIYxBNT_5HgW",
    "JhmIBByY6T9h",
    "k7spOYvvTdZL",
    "qIMN2gHGzKof",
    "aMM0sJy6z7Ur",
    "WY6bUgqZztC3",
    "58gng2f17_P7",
    "C8Us4nEyhRdu",
    "wkZR0gDWRAgK",
    "wPw9kCvASOeJ",
    "WdZcUfvLUkwK",
    "EHH1hP0uzaWG",
    "I_Ph4o7ZW_XD",
    "jOI-E76fZ2bG",
    "Rq2W2fmobmFg"
   ],
   "name": "PyMC intro.ipynb",
   "provenance": []
  },
  "hide_input": false,
  "kernelspec": {
   "display_name": "pymc",
   "language": "python",
   "name": "pymc"
  },
  "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.11.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "867ba48c05011db76db56a12fb95ccd32f7ac276df8f4ae698e0d475911a6ba0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}