notebooks/factorial-hmm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a819aa7b-10e2-4169-9366-ab2c29a87c86",
   "metadata": {},
   "source": [
    "# Factorial HMM with Gaussian emissions\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/lindermanlab/ssm-jax-refactor/blob/main/notebooks/factorial-hmm.ipynb)\n",
    "\n",
    "This notebook showcases a simple factorial hidden Markov model where the latent state factors into multiple groups, $z_t = (z_{t1}, z_{t2}, \\ldots, z_{tG}) \\in [K_1] \\times [K_2] \\times \\cdots \\times [K_G]$ where $G$ denotes the number of groups, $K_g$ is the number of discrete states $z_{tg}$ can take on, and $[K_g]$ is shorthand for the set $\\{1,\\ldots, K_g\\}$. \n",
    "\n",
    "The factorial HMM assumes that each group transitions independently, but that they are all linked in the likelihood $p(x_t \\mid z_t)$. The joint probability is,\n",
    "\\begin{align}\n",
    "p(z_{1:T}, x_{1:T}) &= \n",
    "\\left[ \\prod_{g=1}^G p(z_{1g}) \\right] \n",
    "\\left[ \\prod_{t=2}^T \\prod_{g=1}^G p(z_{tg} \\mid z_{t-1,g}) \\right]\n",
    "\\left[ \\prod_{t=1}^T p(x_t \\mid z_{t}) \\right]\n",
    "\\end{align}\n",
    "Here, we consider a simple model where the emissions are scalar Gaussians and the likelihood is,\n",
    "\\begin{align}\n",
    "p(x_t \\mid z_t) &= \\mathcal{N}\\left(x_t \\,\\Big|\\, \\sum_{g=1}^G \\mu_{g,z_{tg}}, \\sigma^2 \\right) \n",
    "\\end{align}\n",
    "For example, suppose you observe the aggregate power usage of a home. At any point in time, the expected power reading may be the sum of power usages for each appliace that is on at that time. To model this setting, let $G$ denotes the number of appliances (groups), $z_{tg} \\in \\{0,1\\}$ indicate whether that appliance is on or off, and then let $\\mu_{g,1}$ denote the average usage when that appliace is on, and set $\\mu_{g,0} = 0$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a0acb76a-e986-44bb-a496-8cfc4a1cee39",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import ssm\n",
    "except:\n",
    "    !pip install git+https://github.com/lindermanlab/ssm-jax-refactor.git -qqq\n",
    "    import ssm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "33e5a806-8112-48ed-98b2-27c018da8a88",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as np\n",
    "import jax.random as jr\n",
    "\n",
    "from ssm.factorial_hmm.models import NormalFactorialHMM\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "139e213f-9ece-4644-b421-4f05563ff115",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "num_states = (2, 2, 2)\n",
    "factorial_hmm = NormalFactorialHMM(num_states=num_states, seed=jr.PRNGKey(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c0296932-6fd5-4cb8-bef2-3cfe420e92fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFzCAYAAAA9nXBaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABHu0lEQVR4nO3deXhjZ3k3/u8tyZIsS/Iuz2bPJDOTgSwkTCYhYQkQKBDWNqQECoWmC/CmbOXXly6UAi1cZelyheUFkl9LgZemLQ0pgUIaCCUBkpDMJEMWktlnbM9ieZdk7dL9/nF0ZNmWbNmWdKTj7+e6dI1HtqTn6EjPfZ7tfkRVQURERPbjsLoAREREVB8M8kRERDbFIE9ERGRTDPJEREQ2xSBPRERkUwzyRERENuWyugCr0dfXpzt27ICq4vTp0xgbG6v4t06nE7t27YLf71/yu0wmg6NHjyIej6/4mm1tbdi9ezfa29vXVfZ6Onv2LM6cObOmx7a3t2PXrl1wu92retzY2BhGR0fX9JrlbN26FZs2barZ861FLpfD8ePHEYlELC3HavT29mJoaAgOR/2v12dmZnD8+HHUatltf38/BgcHISI1eb71CofDGBkZsboYFTXb+9VM8vk8RkZGMDExUbfX6Onpwfbt2xvyXSt18uRJTExMrP2kq2rL3C6//HJVVU2n0/qhD31IAVS8dXZ26v3336/lnD59Wq+88splH2/eBgcH9Yknnij7PM0gm83qxz/+8aqOpdzt0ksv1eHh4VW9Zi6X089+9rNrfs1yt0996lOay+Xq9C5VZ3p6Wq+77rqaHle9b+94xzt0bm6uIe/PnXfeqW63u2Zlv/nmmzWZTDak7CvJ5/P6uc99TkXE8nNa6fae97ynad6vZhOPx/X3fu/36vr+v+Utb9FYLNbwYyvEvTXHTXbXExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZVMODvIh0isinROQZEZkSkUkRebpwX1ejy0NERGRXVrTk/x3ANICXqGqPqvYCeGnhvn+3oDxERES2ZEWQ36Gqn1bVc+YdqnpOVT8NYLsF5SEiIrIlK4L8KRH5kIgMmHeIyICI/AmAEQvKQ0REZEtWBPkbAfQCuK8wJj8F4CcAegC8yYLyEBER2ZKr0S+oqtMA/qRwIyIiojrhEjoiIiKbYpAnIiKyKQZ5IiIim7IsyIuIT0Q+IiK3Ff6/W0Rea1V5iIiI7MbKlvxXAaQAXF34/2kAn7CuOERERPZiZZDfqaqfAZABAFWNAxALy0NERGQrVgb5tIi0A1AAEJGdMFr2REREVAMNXydf4mMA7gYwKCLfBPACADdZWB4iIiJbsSzIq+o9InIAwFUwuunfr6oTVpWHiIjIbqycXX+vqk6q6n+p6vdUdUJE7rWqPERERHbT8Ja8iHgB+AD0iUg35ifbBQFsbXR5iIiI7MqK7vp3AfgAgC0ADmA+yEcAfMGC8hAREdmSFRvU3ALgFhF5r6p+vtGvT0REtFFYOfHu8yJyMYALAXhL7v+6VWUiIiKyE8uCvIh8FMBLYAT57wO4DsDPADDIExER1YCVyXBuAPAyAOdU9SYAlwLotLA8REREtmJlkE+oah5AVkSCAMIABi0sDxERka1YmfFuv4h0AbgNxiz7GIAHLSwPERGRrVg58e7mwo9fFpG7AQRV9XGrykNERGQ3lma8M39W1ZOq+jgz3hEREdUOM94RERHZFDPeERER2RQz3hEREdlUw8fkReQKEdlkBngRebuIfEdEPiciPY0uDxERkV1ZMfHuKwDSACAi1wD4FIwsd7MAbrWgPERERLZkxZi8U1WnCj/fCOBWVb0DwB0ictCC8hAREdmSFS15p4iYFxcvA/Djkt9ZmZyHiIjIVqwIqrcDuE9EJgAkAPwUAERkF4wueyIiIqoBK2bXf7KQ9GYzgHtUVQu/cgB4b6PLQ0REZFeWdI+r6kNl7jtsRVmIiIjsyspd6IiIiKiOGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbEpU1eoyVE1EogAOWV2OOuoDMGF1IerIzsdn52MDeHytjsfXuvaoamCtD3bVsiQNcEhV91ldiHoRkf08vtZk52MDeHytjsfXukRk/3oez+56IiIim2KQJyIisqlWC/K3Wl2AOuPxtS47HxvA42t1PL7Wta5ja6mJd0RERFS9VmvJExERUZUY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyqpTao6evr0x07dpT9XSKRwJEjR+ByubBr1y643e7GFq6G8vk8hoeHMTk5WfFv2trasHv3brS3t1f8m1wuhxMnTiAajWLnzp0IBoOYmprCiRMn0NPTg+3bt8PhWPk6L5PJ4OjRo4jH4xX/JhgM4rzzzoPLVZuPVDabxbFjxxCLxWryfLXmdDqxa9cu+P1+jI+PY3h42OoiAShfrv7+fgwODkJErC7eArFYDEePHkUul1v1Y0UEO3fuRGdnZ1OVq5GqqQNakVmXZzKZNT1++/bt6Ovrq3GprHPy5ElMTEys/curqi1zu/zyy7WSgwcP6pYtW/SSSy7RkydPVvy7VhCPx/Wmm25SABVvg4OD+uSTTy77PNPT0/qa17xGfT6f3nPPPaqq+s1vflNdLpe+7W1v07m5uarKMzo6qvv27Vu2PK985St1cnJy3cduCofDes011yz7mlbeuru79ec//7nm83n94he/qCJieZkqlevmm2/WRCJRs3NTKz/72c+0q6trTcfp9Xr1rrvuqku57r///jWXq5G3zZs368GDB+vyHljp8ccf161bt67pPXE4HHrbbbdZfQg1VYh7a46b7K4nIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbaqogLyI/sLoMREREduFq9AuKyN5KvwJwWQOLQkREZGsND/IAHgFwH4ygvlhXY4tCRERkX1YE+acBvEtVjyz+hYiMWFAeIiIiW7JiTP5jy7zuextYDiIiIltreEteVf9jmd/9ZwOLQkREZGtNNbueiIiIaodBnoiIyKYY5ImIiGzKsiAvIj4R+YiI3Fb4/24Rea1V5SEiIrIbK1vyXwWQAnB14f+nAXzCuuIQERHZi5VBfqeqfgZABgBUNY7yCXKIiIhoDawM8mkRaQegACAiO2G07ImIiKgGrMh4Z/oYgLsBDIrINwG8AMBNFpaHiIjIViwL8qp6j4gcAHAVjG7696vqhFXlISIishsrZ9ffq6qTqvpfqvo9VZ0QkXutKg8REZHdWLHVrBeAD0CfiHRjfrJdEMDWRpeHiIjIrqzorn8XgA8A2ALgAOaDfATAFywoDxERkS1ZsUHNLQBuEZH3qurnG/36REREG4WVE+8+LyIXA7gQgLfk/q9bVSYiIiI7sSzIi8hHAbwERpD/PoDrAPwMAIM8ERFRDViZDOcGAC8DcE5VbwJwKYBOC8tDRERkK1YG+YSq5gFkRSQIIAxg0MLyEBER2YqVGe/2i0gXgNtgzLKPAXjQwvIQERHZipUT724u/PhlEbkbQFBVH7eqPERERHZjacY782dVPamqjzPjHRERUe0w4x0REZFNMeMdERGRTTHjHRERkU01fExeRK4QkU1mgBeRt4vId0TkcyLS0+jyEBER2ZUVE+++AiANACJyDYBPwchyNwvgVgvKQ0REZEtWjMk7VXWq8PONAG5V1TsA3CEiBy0oDxERkS1Z0ZJ3ioh5cfEyAD8u+Z2VyXmIiIhsxYqgejuA+0RkAkACwE8BQER2weiyJyIiohqwYnb9JwtJbzYDuEdVtfArB4D3Nro8REREdmVJ97iqPlTmvsNWlIWIiMiurNyFjoiIiOqIQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsSlTV6jJUTUSiAA5ZXY466gMwYXUh6sjOx2fnYwN4fK2Ox9e69qhqYK0PdtWyJA1wSFX3WV2IehGR/Ty+1mTnYwN4fK2Ox9e6RGT/eh7P7noiIiKbYpAnIiKyqVYL8rdaXYA64/G1LjsfG8Dja3U8vta1rmNrqYl3REREVL1Wa8kTERFRlRjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbKqlNqjp6+vTHTt2WF0MW8lmszh27BhisVjxvsHBQYRCoZq/VjKZxJEjR5BOp9f1PH19fRgaGoKI1KhktZVOp3H06FEkEonifTt27EBvb29Vj5+ZmcHx48dRr0RVLpcLu3fvhs/nW/K7VCqFI0eOIJVK1fQ1nU4ndu/ejY6OjnU9Ty6Xw8mTJzEzM1ObgpXYvn07+vr6av68rUJVMTo6inA4vK7naWtrw+7du9He3r7kd2v5fDkcDuzatQuBwPxGbHNzczh69Ciy2eySv6/3ecxmszh+/Dii0WjxvnrVmQBw8uRJTExMrL2yU9WWuV1++eVKtRUOh/Waa65RAApARUS/8IUvaD6fr/lrPf3007pjx47ia6319s53vlMTiUTNy1crIyMjetlllxXL63Q69etf/3rVj7/zzjvV7Xav+32qdNu0aZM++uijZV/76NGjesEFF9T8Nbu7u/XBBx9c93s7Ozur119/fc3L53A49NZbb63L575VpFIpff/737/u93LLli36+OOPl32Nw4cP686dO1f1fIFAQH/yk58seJ5f/OIX2tfXV/Y83nbbbXV9nyYmJvTaa69dUGfecsstdfvsFOLemuMmu+uJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGh7kRSQoIn8jIt8Qkd9a9Lv/0+jyEBER2ZUVLfmvAhAAdwB4s4jcISKewu+usqA8REREtmRFkN+pqn+qqv+pqq8H8CiAH4tIrwVlISIisi2XBa/pERGHquYBQFU/KSKnAdwPwG9BeYiIiGzJipb8dwFcW3qHqv4zgP8PQNqC8hAREdlSw1vyqvqhCvffDWB3g4tDRERkW1xCR0REZFMM8kRERDbFIE9ERGRTlgV5EfGJyEdE5LbC/3eLyGutKg8REZHdWNmS/yqAFICrC/8/DeAT1hWHiIjIXqwM8jtV9TMAMgCgqnEYmfCIiIioBqwM8mkRaQegACAiO2G07ImIiKgGrMh4Z/oYgLsBDIrINwG8AMBNFpaHiIjIViwL8qp6j4gcgLEpjQB4v6pOWFUeIiIiu7Fydv29qjqpqv+lqt9T1QkRudeq8hAREdlNw1vyIuIF4APQJyLdmJ9sFwSwtdHlISIisisruuvfBeADALYAOID5IB8B8AULykNERGRLVmxQcwuAW0Tkvar6+Ua/PhER0UZh5cS7z4vIxQAuBOAtuf/rVpWJiIjITiwL8iLyUQAvgRHkvw/gOgA/A8AgT0REVANWJsO5AcDLAJxT1ZsAXAqg08LyEBER2YqVQT6hqnkAWREJAggDGLSwPERERLZiZca7/SLSBeA2GLPsYwAetLA8REREtmLlxLubCz9+WUTuBhBU1cetKg8REZHdWJrxzvxZVU+q6uPMeEdERFQ7zHhHRERkU8x4R0REZFPMeEdERGRTDR+TF5ErRGSTGeBF5O0i8h0R+ZyI9DS6PERERHZlxcS7rwBIA4CIXAPgUzCy3M0CuNWC8hAREdmSFWPyTlWdKvx8I4BbVfUOAHeIyEELykNERGRLVrTknSJiXly8DMCPS35nZXIeIiIiW7EiqN4O4D4RmQCQAPBTABCRXTC67ImIiKgGrJhd/8lC0pvNAO5RVS38ygHgvY0uDxERkV1Z0j2uqg+Vue+wFWUhIiKyKyt3oSMiIqI6YpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIimxJVtboMVRORKIBDVpejjvoATFhdiDqy8/HZ+dgAHl+r4/G1rj2qGljrg121LEkDHFLVfVYXol5EZD+PrzXZ+dgAHl+r4/G1LhHZv57Hs7ueiIjIphjkiYiIbKrVgvytVhegznh8rcvOxwbw+Fodj691revYWmriHREREVWv1VryREREVCUGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptqqQ1q+vr6dMeOHVYXg4iIqCEOHDgwoar9a318SwX5HTt2YP/+dW3IQ0RE1DJE5NR6Hs/ueiJaIprM4IFjExiZiltdlJpjKu+VqSoyufyKf5fN5ZFI55DK5hpQKlqLlmrJU33c+dgoHjk5jWt29+OFu/vg9zT+Y3FsPIYbv/IQLhjw44W7+/DW521HZ3tbw8thB9946BTGoyl88NcuAGBU2I+PzuLsbBLxdBYv3RNCd4e77GOzuTzecttD2H9qGqqA0yG48YpBvP9luzEQ9DbyMFb0Dz88jC/ddwxd7W3o9rnR5WtDl68NbpcTLodgai6NMzMJtDkduHBLEEFvGx4+OYmnz0YxEPBgR18H5lJZnJ5J4reeN1R8vzayb+0fwbcfPY2nzswikszC73GhP+DBRVuCuHhrJ/KqmIlncGQsiidORzARSwEAHAKc3+/HZYNd+MhrLkSnr/x3N57O4v//6Qls7vRiV8iPx0dn8aOnxxCOpCACXLGjB3/96xc38pBXLZHO4cP/+QT6/R5cNtiFlz4rBG+b0+piVcQgT/jaA6dwcGQG//KLYWwKevHgn10LEWloGR4bnsFELIWg14XP3H0IAsH/esnOhpbBLr77yzM4PZ0oBq2nzkTwhi/+vPj79127Cx98xZ6yjz07m8QjJ6fx+ku34DeeuxU/ORTGvzw8jAMnp/Hff3RNQ8pfrf2nptDjc+PFF/RjOp7GTDyDkxNxpHN5ZHJ5dPvcOK+vA8lsHj85FEYkmcXeoS783gvPw3g0hZOTcwi2t2F0OoH9J6esPpym8A8/PIy8Aq+9dAsGAl7MJNI4O5PEY8Mz+N7jZwEAbpcDO3p9ePEF/Rjq8cHtciCRzuKxkRn8x4FRvOxZIVx3yeayz3/foXH8/Q8PL7hvZ38HdoX8ODIWw7/tH8FfveGihtc/q/HE6Vl8+9HTEAFUgT986U7871c+y+piVcQgT4ilsvi1Cwdwfn8HvnLfcUQS2YpX4vUyOh2HCHD3B67Bvk/8EGdnEw19fTsZj6YwHk1BVSEiGC50uX/prXvxse8+hdMzyYqPDUeNltlv7N2Kl+4J4aXPCqHL58Yt9x5BMpNrqhZLOJLCZYNd+PQNz1nxb1UV+ULPxGI3ffVhTMTS9ShiS8nnFeOxFH7/RefjT161NGjNJjLwuBwVPwMTsRT2feJHxc9QOWMR47P3H+++GuPRFHYPBLAr5AcA3Hb/cXzy+09bUv+sRjhqHMNdf/hC/K9vHsDIVHPXVRyTJ8SSWfT43LhoSyeA+Q9xI41MJbAp6IXb5cBA0FusDGj1wpEk0rk8IoksgPmK9crzerC5s33Z8zte+F0o4Cnet73XBwA4PdNclVk4mkJ/STmXIyJlAzwA+L1tiKWytSxaS5qOp5HJ6YJzX6qzvW3Zi7wenxsuhyz73R2LptDmFOwd6sZ1l2wuBngACAWN17Wi/lmNsYhxEbOtux2bO71NX14GeUI0mYHf68JA4cttfogbaXQ6jm3d7QBQCPKNL4MdxFJZzKWNSVBm5RMuVKzdPjcGgh6El3lvzVZYKDA//j7UYwT54SaahJfM5DCbyFQMSKvh97gQTTLIlzv3q+FwCPoDnmVb8uFICv1+DxxlLrjMC7blHt8MwtEk3E4HunxtCAW8y36fmgGD/AaXyyvm0jn4PS6EChOrrGhFj04nMNhtBJNQ0IPxJv+iN6twybkzK8uxSLJYsYYCXowt0/IIR1JwOgS9JRPzBgtBfrSJgrw54cts/a1H0OtCNJlZ9/O0umKQX8d7GlopyEeTxXpm6WO9xb9pZuFICqGgByKCUHD5420GDPIb3FzaaMEEvK5iq6jRH9psLo9zkeSClnw4mkQ+z6VOq1V67szKcjyaKlasA0EPZuIZJDPllzyFo0n0+d0LWlr9fg88LkdTteTX2+os5fe4kMrmkc6uvGTMzswLxIF1vKehoHfBhebS10hV7H0ZMLvrm7xlPBZJFo8hFPAavWdNPNxjSZAXkd8UkadEJC8i+6woAxliyfkg3+FxIeBxNbwlf3Y2iVxesa3Qkh8IeJDJKabjnAy1WguCfGS+JV+slArBvlJPyXiZcW6HQ7Ctu72pJhiZx1btmPxy/F5j/nEzV9SN0IiW/Fg0WXEppt/jQnubs+lbxuFoqngMxQuTJi6zVS35JwFcD+B+i16fCsyxSL/HmM3ab0FX+ci00UIsbckD1swNaHVmK8oh8xVPaaU031tT/kIuHE2VbR0P9fiK56kZlJsguFYBr/HZ3+jj8uFIEgGva10rKEIBL6bm0mV7RZKZHGbimWJgXKxVur8Xt+QBLNt7YTVLgryqPq2qh6x4bVooljLGIs3WzECg8TPbR6eNFqI59lucG9DkY3PNaDyagtvlwGCPD+FoaknFutIFlBHkl1bCgz2+puuudwjQ66/NxDsAiKY29rh8pXO/GubnbDy29PM1XsUQSyjgaeqAmUjnEE1mFwx/AcaqgWbFMfkNbr4lb1R0oaCn4cF1dDoBhwCbOhd1gTXxl71ZmRW1WVkurlgHlplcmcsrJmMVgny3D9FkFrPx5giE49EUev2eisviViNQuMCNbfCW/Fikcld6tULLfHfN3qPlhgNCAW9TT7o1j2G+Z2wDt+RF5Eci8mSZ2xtW+TzvFJH9IrJ/fHy8XsXdsMwgHzRb8kFjSUgj83uPTsWxubMdbU7j49hv4VK+Vmd2JYYCXozHUksq1m5fG9qcUva9nYylkNfy49yDTbaMrhatTpMZ5Dd8d30N3lMz6JX7fJnzKJZrya+0BM9qY8VjMN6nYLsLbpejqS9M6pbxTlVfXqPnuRXArQCwb98+TreuMTMJiNldHwp4kMrmG5p1anQ6ga2F8XgA8Lic6Pa1Nf1SmmYUjqawq9+P/oAH9x9OlVRKRsUqYiyjK/fempVrf5lKeLDHOD8j03Fcsq2zXsWvWjiarMmkO2C+F2sjJ8RRVSPI16glP17m82X2HlUakzcfH0tlEU9n4XM3X0LWxS15EcFA0NPUybvYXb/BxZZ01zd+rerodLy4Rt7EhDhrE44kEQp60B/wIJrK4tSk0fIurVhDFRLijC8zu9psyTfLrnTLLcVaLfMCN7qBg3wkkUU6m1/3e9rb4YFDyrfkzWx33b7ymyMBpd3fzfndN49rwfcp4G3q3gerltD9hoiMArgawH+JyH9bUQ4yKjYRoMNtTrxrbFd5OrtwjbxppfW2tFQyk0MkmS2OyQPAk2dm4XIsrFhDgfItj/AyM9aDXmOHt2bors/lFZNz6ZqskQeMYwOwoRPijBWHddb3njqLWe/KfL6WyXZnsipXR7XC0STcLseCHTLZki9DVe9U1W2q6lHVAVV9pRXloEJKW7er+MVrdNa7s7MJ5BVLgvxAwMOW/CrNt8S9xfP45OlZhAILK9ZKewOstPZ8sNuHkWnr18pPzaWRy2tNst0BgMflgMshG3riXXjRWPN6VGrZLpftrvjYJs9fb/Ygle6Sx5Y8NbVYMlvsrgQafyVtLp/bVqa7fjyWQo5Z76pmBu7SlvypyfiSinUg6EUkmV2S9S4cTaGzvQ0eV/l10kM9vqborl+ux2EtRAR+r2tDj8nX8j01WrZluusjyWXH443Xb+7u+nA0ueQ9CgU9iCazSKTLZ5G0WvPNbKCGiqWyxfF4AOjwuOCvY9a7rz94Eo8Nz2A6nkaHx1VMXWtO7DINBD2FbtnyyVloqdJUr6UV0eJKqbgRSCSFoV5fyeOXVmCltvW044e/GkMml4fLIYgksjgzm8DZ2QTOzhrL9abn0sgrcPn2buwd6kYslcV4LIW9Q13FpDO1Os5aTbwDjBn2G7klX5yguc7uesCYuHlwZGbJ/eFoCs87r3fZx5qrP5q1ZTwWSWF3yc55wMKc+9t7O6wo1rIY5De4WCpbXEJkqtcGMdlcHh+96yl0tbdhW7cPR8ZiOD2TQMDjwqZFlUtxAmCEQb5a5hyGUNCD7sK2n9m8Lln7PFCSbKg0yBs57isHzh29HUjn8tj94R8Un3uxoNeFvALfeOjUgvvfd+0ufPAVe9Z8bKXGq1iKtVp+TxsiGzjIh6NJdLidCy741yoU8GAilkYmly8ui10p251JRNDvLz+m3wzGIkm8cFffgvtKez8Z5KnpRJLZBZNIgPplvZuIpaEK/PEr9+Ctz9sOwKhcsjmFy7lw5Ghh+lXrl2y1gnA0BZdD0OMzNpjp83twLrK0dV5pI5BwNIUrdvRUfP7XXboF6WweM/EMEpkcejvc2NzlxeZOY1/t/oAHbU4HcnnF02cj+OXoDLra3fj4d58qDsvUgplNraYteY+rmP1xI6rF8jmTeRE5EUthc6fRQ1dNtjtTf7A5E+KY2e4Wf+6WSzDVDBjkN7hYMoNtXYtntnvw6PB0zV+ruE625Ite6UvP/PWrFy5sLjM/idII8otb8vMJS+YrJXOd9HKB0+9x4R3P37FiOZwOwcVbO3HxVuPi7B9/drymWRTDkSSC68yxvpjf62ra1mMjjEeWP/erESoZDjKDfDXZ7kofPzxp/dyPxRavkTeVHm8z4sS7DW7xmDxQv6x3q9nlaj7r3cateFdrccYy8+f+Re93MetdSVCr1TrpcszPU63UstVpCnhdGzrj3XK7w61WuZbt/PrylV8jVGEJntXM+mvxkEOXrw1up6Np5xG0VEs+k1vdfs/xdBb3H57AszcHMNTjW7DsoZZlcohUzKGdyuZwZCyGgaC3qitlc1JTPcpaTjRZZky+JOsdYKRurEV55md/r/xFb3M60Od349HhGXz3l2fQ0+HG887rWdKtT8CBU1NIZfMYmYpjZ//8pCAzc93i/cHNrHfDk/FiZjGzUq1lF7gpFPDg50cnavZ84aix3rqW/J7WnHg3lzIuzrp8bWv+jqpqTZMLmRfxT52JYEtXO1LZPA6cMnoGq3mNUMCL6XgGpybnEElkceGW4JL6dWoujbsOnkaHx4XNne14/s7eZdff10Kl+kukkBugygaJqmImnkEmn4dTBD0d7iXnLpdXpLI5uGtQ37VUkH/mXBTv+sZ+3PySXbh0sGvFv/+PA6P4y+88BcD4cO3ZFMDOfj/efvV2nN/vX+HR8x4dnsat9x3H5FwKiUwOO/v9OL/PjydOz+DnRyfhdAguHezE5s52jEdTmImnoTASvRwbjyGTUzgEuPK8HuwZCGA8lsJkLI1kJodUNg+3ywG304GxaBKnpxO4fu82/O1vXrpiue4/PI6nz0bwrhfvrPpYSuXying6t2AJHTA/6e2qv7kXiUwOf/+mS3H93m1reo1S4WgKIkCfv3LGq1Ln9XXg/sPjuP+wsWdBf8CDGy7fhj9+xZ6abExiBycm5vDGLz1Y/P+LL+gv/jy/h/zSinVrdzt+8OQ5/ODJc2hvc6LDzHhYh0mOocKSvUQ6h3b3+rvYw9Ek9g5116Bk8/xeV0tmvLvulp9ieCoOb5sDmzvbsSnoxeZOLzZ1zjcqkpk8Do9F8cuRGcTTOfT63ej1e9Dnd8Pb5sTTZyNIZHIrToqrVm+HG26XA7fcewS33HukeH/A41o2253JLMeLP/sTAMAXfuu5eO1ztiz4mzsOjOKT33+6+P+v3nQFXronVIPSV3ZutnJa3oGgBz8+FMYNX3oAOVXk84qcKnJ5I6iHgl6c39eBiVgKDx6bxORcuvjYbl8b9mwKoM3pQCqTx3gshdHpODK52vSkVhXkReQqAJ8H8GwAbgBOAHOqGqxJKarUH/Dg4RNTeOTkIzjwFy9f8cr1xMQcfG4n/vzVz8ajp6ZxdDyGbzx0CqlsHn9z/SVVv+6/PjyMHx8KY+9QF3o6PHjkxBS+c/AMtnW34zf3bYOqcSFwNBxDKOBFT4cbDgEcInjxnn5cvKUTR8Ix/OCJs/j2mdMIBTzo7fCg0+eG2+lAJpdHOpvHZYPdcIrgqTORqsr1b4+M4MfPhPHOa85f01V8MW/9ou76q87vwasu2oRQ0IPbHx7G4bHYqp+7nHAkid4OT9Wt8X/8nStwejoBl0NwbDyGrz94Cl/6yTG87jlbcOGWhn70mtaJCePc/PWvX4xNQS/2bZ8Pfm/cuw0+txO9HUsr1r+94VL84oRR2UwWLjpzqnhOHfLSm120tVhipKrGKoAa9zgEvW1IZ/NIZXMV8wQ0m1Q2h+GpOF58QT8uGPDj7GwS52aT+MWJKYxFkgtWP/QHPLhssAtd7W3Fc34sHEMslcWegQD+4EXn4Teeu/4LeQBwOR24/Q+uwpmZBNwuBzwuh7H9cbevqtb2qy7ehImYkbPhI995qmyWxXORJHxuJ7717qvxms/9DCfG5/DS2izeqGh0OgG/x7VkojIA3HjFIL796Gk4HUavrtm76xCBiJH0a//JKfg9LlxzQT8u3toJt8uBdDaPw+eiOByOIpnJw+Ny4MLNQbziogF0+9zIZPN436fXV+5qW/JfAPBmAN8CsA/A2wFcsL6XXr1NQS9uvnY3/up7v8LUXHrFvaRHphIY7PbhbVdtx9uuMmZzv+GLP191Qo/hqTgu2dqJf33n1cX74uks2tucqwquH/y1ld+yP7/zCfz3k+eqer5wNIlEJoeJWHpN3axmkF/aXe/Fl3/7cgDAj341VrPxsXA0tarWQtDbhuBm4wu1eyCA/oAXD3zpAYxFk7gQDPKA8RkHgFdeNLCkFT7U66vYyzPU61uwfK6earnEKJrKIpnJ17zHobhJTTILj781grw5A/26izfhzVcOLfhdPq+YSWTgEGPoy+deXV21Xpdv78bl29fW29Llc+M91+4GAHzm7kNl53OY808u3BxEh9uJken6T9Q7PZPAtu72su/jjVcM4cYrhso8ap6qrukcvG/Vj1io6g5/VT0KwKmqOVX9KoBXrfO112R7oWI6VUWgHp2OL0myMtTjW3X+7ZGpBIZ6FlaIPndtxqkXCwU8mJxLVzX/wJzosdYPuJmre7kkJf2B2q2ZXynZykq4z/xSo9NxeFyOmo9R11ItlxgV06/WqGvZ1Io70S03kdXhMMZ6u3xudHjqU1c1Qn+FnB3hiLELoYhgW7evpks0KxmdTmDropVIq2HVOag2yMdFxA3goIh8RkT+aBWPramhKnfDUlWMTMWXpEsd6mnH6ZkEslVO4ktn8zgzmyjuwlVvZgtlpcBqTpYB1r4z2OId6MrpD9RuzerYOhPb9Df5UhUrjE5Xbl00i1BxpcT6z1txgmCNL2pacU/58ComsraqSolxxmPzdclgT3tD0i2PTseX7LHRCqoN1L9d+Nv3AJgDMAjg+noVajlmsD21wjrK6XgGc+nckpMy1ONDLq84O1tdq+L0TAKqWNKSr5dia3WFwBpLZZEo5B5f61VsdNFe8uXUKvtdNpfHZGx13fWLmfvM13LNdasbmV56Idts5pcYrf+8Lbcd7nr4WzHIFxPMNG8vznqFKiTGKV3Xv63bh9PTiZov+S01m8ggmsw2/XetnGqD/K+ralJVI6r6cVX9IIDX1rNglXjbnBgIelbscjev7Ba3wM3/V9tlb/5do4L8/AYNy1eIpRcB623JB5cJ8v3+6ocPljNZyGnev861uKFAbddct7rR6cSSIalmIyIV97BfLbPC769x6zXgMYasWqq7PpKCQ7Di3KRWZqyZX/i5SaRziKayJUG+HdFUFrOJ+mUsPF1oSG21cUv+HWXu+50almNVtvd0rJgRyWzdDi7prm/uIG+2dMdWaD2bFaZD1jMmb3bXVx6TN1tMk7F0xb+phlnegXW2OkJBz4rvzUYRTWYwE8+0ROui0h72qxWOpuBxOZa9MF2L+e761kltG44m0ef32Ho5aSjgQTydW3DxNb6oB8P8/NdzXH60UMfarrteRN4iIt8FcJ6I3FVy+x8AU40p4lKDVUyeMwPf4lbO5s52uBxSdZAfmYrD7XI0rEus1++BQ4DxFVvyxu+ftSlYnGG9Wmau7uW6682xz/V2tRYTSdSgJb/Se7NRVLqQbUYDwdrsuV064aqWzO9AS7XkV9hQyA7m5+HMf+cXJ24y6/h6jstX2hK7Fax0OfwAgLMA+gD8Xcn9UQCP16tQK9ne68MdjyaRzOQq5q8emYqjy9e2ZOa40yHY1t1efUt+Mo7B7va6Z1MyOR2CXn/5/ZhLmS3jy7d34/aHh5HL66qv6GPJLEQA3zI5wM2gvN5x+VqNHw4Eje67fF4bdk6a1XzF0/yti4GgFz+rQda7xal7a8WcfNpSY/KRFDZ12nfSHbBwIrKZwGzxZjeNaMmfnkmgvc2YE9Rqlg3yqnoKwCkAVy/3d6slIp8F8DoAaQDHANykqjPVPt7sOh+djmNXKFD2b0amExVbOIM9vqqv+oan4g3rqjcZgWzllrzH5cCFW4LI5hVnZxOrvsqMJI289csFy/5AdRMBV1KrtKmhgAfZvGIqnkafjcciq2F+hlshyPcHPIgms8U0ums1Hk0tSN1bK942J9xOR2sF+WiqLsmLmkmozETkxUsHO9vbEPS66rpW3pxZ38yrWCqpakxeRK4SkUdEJCYiaRHJiUh1adnK+yGAi1X1OQAOA/iz1TzYTOKxXGt8dGrpGvni46tcK28uw2t0kA8FvCu35AtddeaFzFq67GOpLAIr7B9tpqBdb0t+LJJCb4e7uL/0Wg2U7DO/0Y1OJ+BzO9FTJqNds6nVeatnF7Xf2zrbzWZzeUzO1adXo5mEyjQyxqMpOAtbKpvqvVbeXKraiqqtcb8A4C0AjgBoB/D7AL641hdV1XtU1bxkfgjAqvIpDq2wjC6f18JJKR+ch3p8mIlnVpyNORPPIJrKNmyNvMnskl5OuLDmvDgetYar2Fgyu+x4PGAsW+vyta17TH48mqzJzmGh4sREjsuPtFDrotqloctJZnKYTWTqFthaaZOaiVgaquuf49LsOtuXLr80Jhy6F/RA1nut/Oh0oiVn1gPNkfHudwH8YDUP6O1wo8PtrNgaH4+lkM7lMVjhpFSbUKfRM+tN/QEvJudSyybsCUeTGAh6sKWr3Zhhv4YPeLltZssJ1SDr3ViNdrkqjtGxJW8sn2uRiUDl9rBfrfnlc/UJ8q203Wxxf3abt+TNHd7GF3XXL/4MmC35eqyVjyaNBmErTroD6pjxTkR+JCJPlrm9oeRvPgwgC+CbyzzPO0Vkv4jsHx8fN+9bdly9OFZZITgPrjbINyjHtykU8EDVuFqvxJiA5EWb09h9ai1BPprMLJvSdr48658ZbV6UrBf3mTeoKkanWicDV3Fp6DrOW3jRhKta83taZye6+fS+9m7JA0tTa49Hl2bOHOxuRyKTW7C7W62cnmmdCa7lrCfj3RuXe4CqvlxVLy5z+w4AiMjvwEio81Zd5vJLVW9V1X2quq+/f34bzaEeX8Xu+uLyuUrd9VWM6VfzPPVSumtXOYl0DtHkwmQQI2sYj4qmVu6uB1DYK3ntQT6XV0zE0jWpnL1t5vDBxm7JRxJZRFOtk4Grs70NbpdjXT1CjWjJt0p3/UbIdmcKLap/yq2wqOcM+9HCfKf15K23UlVBvjDLPgDAY2a8K3Tfr4mIvArAhwC8XlXXNJCyvdeYPFfu+sCchFbpyivobUOXr62qrHl9fndxr+1GWSnX9+KuuqFVrBYoFUuuPPHOfJ3xWGrNXWGTcynk8lqzCVO1SqxST2ORJI6Nx5BI5+ry/JXyQDQrEUEo4MG5dXXXm7kW6hXk2xBtkYl3Zh2wEVaYhEpWG+XyisnY0u76antnTfm8sed7NeYT4bTGBfViy9bwYszo+SiMFryjcFcWwOdV9a/W8bpfAOAB8MPCpKGHVPXdq3mCoR4fUtk8/uDr+xEKevHea3dhc6dR4Y1Ox9Ef8FRcQ28+/n+eCeOvvvsrtLsdODkRx+i0kfjG53Yhlc3hV2ci2Bmq/XKdlazUkp9fQmJu0OBDOJrC1x44CZdTcDQcwzNno3jH87fjVRdvBgA8cnIKX7nvOL78tr3F/dxjqeySbWbL6Q94kM7mEUlk0bnMOtHvPX4GD5+YQiSRQV6N1ls6m8cDx4310eb5Wa9aJVZZr58eGceJiTm8/eodS373+i/8rHiR9sa92/B3b7q06uc9Nh7Dzf/3Ubicgg6PC5FEBuFoCgGvC3sGAtjU6cWJiTkArVXxbO704jsHz+DHz4Rx6bYufO13r1xVbodwtJDGtaN+E+9aZ0w+hZ4ON9wuS/YJa6h+vxfT8QzS2TxmEkZ67KUteaNuueXeI/j3/SPo83vwvPN6sKnTi0dOTuHps1Fs7/XhgoEADg7P4J5fncNcKoet3e3YVrj1B7xIZXKIJDMIR1I4O5tET4cbkWQGHpejuNKo1axUw/8RgBcAuEJVTwCAiJwP4Esi8keq+g9reVFV3bWWx5W65oJ+PH9nL05NxvGjp8PYHfLjphecB8BoAW9eIUnEb+4bxNceOInbHx5GKpvDUI8Pgz0+ZHOK6XgaHpcD+3b04Pq9W9db1FXr87shUnm5UXE8rvBBv2ywCw4BPnrXUwCA9jYnsvk8+gOeYpC/79A4fvT0GMLRFLZ0tSOZySGezqHLt/IH17xqHo8llw3yf/GfTyKZySEU8ELE2NQhn1dceV4v/uBF5+PFF/RXfOxq9Ac8OBaO1eS51uMbD57CgVPTS4J8Ip3DWCSFV1+yCadnkvjFiclVPe+BU9M4NBbFC3b1IpNTbOv2Ye/2bszGM3j6XAQPn5xCh9uFK3f0YJcFF6Fr9bHXX4T/eSaMQ2MxfPeXZ3D/4XG89Fmhqh8fjqTQW8c0rr1+N2YKwaTZg2e4RhNZW4HZczMRS2GqMOa+uCXf4XHhjXu34Ug4irlUFk+fjeLOx04DAFwOwXl9HXjg2ASSmTz8Hhde9uwQNnV6MTqdwOh0Aj/81RgmYka9H/C2IRTwFFc5PXMuiuds62yJVSzlrBTkfxvAr6lqMVWVqh4XkbcBuAfAmoJ8LWzv7cC//MFVUFVc8Bc/WJIsYcsKQf63r9qO375qO1QVubwWW7fNwOV0oLfDvUxLfmF3/TUX9OOpj78Kc+ks0tk8BoJe3PiVB5dMVjEeawR580KhmvHN0i1eKyUfSmZymIln8MevuADvuXZ3lUe6NgNBL8Zj1me9G4umMDmXXhIUJmLGe/uSPSEcDcfwzw+chKpWXUmY5+of33HFsr1RreaiLZ24aEun0btzdAK3Pzy8qiBvbC9av8BmzhmZnEvVrNepXsI1WpLaCkrXyk/HzSC/9NhLe8tUFcfG5zAWSeKywS50eFzI5vIYnopjS1d72e9VpayhiXQOjuYJD6u2UtHbSgO8SVXHATRFfj8RMfYcjiyafVnluJ2INFWANy2321o4moLLIeguaYW3u53o8xtL6pwOY9nJ4rWlwHwO6NUswSkuW4tV7iJfnGqynkIBDzKFHhcrmTn0Jxa9L+GSCWIhc6hjFd3AY5Ekgl6XrQJ8KbfLgRsu34Z7nwmvuNtiqXA0WecgP38x2+w2VEu+ZGfO8Uh1Ew5FBLtCfrxgV19xTpXL6cD5/f6K36tKPUTtbic8rtb9Lq4U3ZarRa2tYUv0B70LJmZMzaWKG6u0KmO3tQot+cJeysu1Yhdv0RguacmX/ltNUO6vovIrpq1twIYZ83MWrKuM83ld8l6airPA/Z75oY5VJO8JR1K2b6XdeMUgcnnFtw6MVv0YMwFUvdQqhXO95fOKiTr3ajQTs8E2HksVGxr1WmFhRyt1119aIX2tAGiaWqjf7ynOgJyMpYx9y1v8QzAQ8GL/yWl84nu/QpvLgZ39fjx7cwChgBfnIokVg0Ao6EU0mS1u4jO+OMgXWlDVrF0Pel3wuBwYjxkJekRkyVXv/FayjWnJA8Cdj53Gtx8dxchUAuFoEh0eF17+7AE8f2cvHA5BOptHJmfc0llFOpfH1FwKpybjiKdzuGhLEBduNnL/JzM5XLSls+qx2Ol4GtnC7NzFM/3NiihUyNcOLD/UsVi9W6zN4Px+P553Xg/+9ZFhBNvbcCwcQ8DrwmCPD0OFW2d7G6LJLFLZHPoDHiOw1fEicj5PenOu3FBVHB6L4fBYFNm82v4zYurtMOYoPXUmgkgiY+ternpYaYOalngnQ0EPHhueBlDaVdo01yBr8rzze/D9J87i9oeHkc7lkcktXO7xigsHln38fAvSGIM3u5TNFuVYmS7/SkQEoaAHX/35Cdz20+O4aEsQ33vvixb8zfxWsvWveMz0krfefxwelwNDPT70Bzw4M5MoTj5cjgjQ5nAgvSij4F+85tn4/RedX1UZxiJLe0lM4+YscL+nmFxlNa3DcDSFK3b0VP33reqtV23H+25/DB/5zyfhczuRzOSw0qqmel689/k9y054tdp3Hz+L993+WPH/F2yq7qKx1bmcDmwKevEvvxgGAFy0JWhxiVpLYxeA10ko4MHkXBqZXN423TnX792G6/caKf1zecWJiTkcOhfF5FwKkURmxQlL812PSXjaHMXK06zAqunyL/Xea3fjseFpnJiYwy9OTC3Z5tecJ9BTxUXDem3ubMe/v+tq+D0u7B7wL9j05th4DL8cmYHTIXA7HXC7HGhzGje3y4HO9jZs6zbmLRwei+LIWAwelwN/fucTODJW/Yz9BfMdFrfko0n0dBizwEOB1bUOVbVu26k2m9c9ZzO2dHqxuasdWzq9yOYVZ2YSGJ6KY3gqjkgii2C7C21OB8Zmk5iKp/HKizbVrTxtTgd6fO5l555Y6dC5CJwOwbfefTUGu30tX8etxldvugKnJuPwtjmxZ2BjXNzUik2CvNFqn4ilqp6Y0UqcDmMSyWqWS5VOIjKDsUjpmPzquoTftG8Qb9o3iDsOjOKh41M4PZNYsOWnmU+6UbPdrzyvfEt3Z7+/6q1IzdneAHDrT4+vapMf82LJUablN16SW9vvccHb5qi6dRhJGCskNkIFLiLYV9Jj0eYUbO/twPbeDsvKtN7sjvU0MpXAli4v9g51W12UhnvWpiCetYkt+LVovmnla1DaNW2Xlvx6FWekRlPFwL6z319sURqBaPVDGmZK4MWZpcYirT2OPNjtW12QL7yPO/v9SyZIlgZ5I9Nb9cl7xopZ3Vp7uKlVGXnSm3NMfngq3jIbElHzsEWQL221jhcyg230iRk9HW44HWJc+BRaJhdtCWIilkYur0ZQXsP4eXH/+kU5oo1li60bmAZ72nFmJrnszn+lxiIpdLa3YajHV74lX7K6I7RoOeNyFic6osaqxWZM9TI6zSBPq2ePIF+yV3U4mtzwrXjA6OI3E+qYvRsXbQkWA/x0PLOmmfChgAdul2NJS77Vx5EHu33I5RVnZ6sMxoVd9ULBhUsVVRXji3JrL/6blZ4XYJC3Siho7HhWbV7zRomns5iIpVtmrwJqHrYI8uYmDeFosrANIStIYD64hAvJVcxWwK/ORIq/Xy2HQ4xd70qCfDqbx9RcbXaZs0pxg4squ+zHCmu2QwEvJudSxR6A2UQGmZwuDPIBb7E3ZSWL9yWgxgoFPMjmFTOJ5tqoxtxdbbDC9tlEldgiyLc5HejpcBtd02sca7ajUMCL8cKYfCjoLQaOJ07PFn6/touhxePXZk9BLfaLt4p5AWRuK7kSM6tiKOiBKjARSxfvBxa+t/0BYyldNTvShSMpdLid8Dd450MyzM9laa5xefOimkGeVssWQR6Yz/AWjrZ+trtaMd8Tc4zYDDxPnTGD/NouhgZ72jE8OR/kww1cI18vm7u8cEh1LXljmVsSoYC3OORh5gkot+f5apbRbaSc5M2oOPTXZDPsi0GeY/K0SrYJ8v0BD05OzCGezrV0sKml/oAHkzFjy8RQcD7F6pOn195dDxjb9EaSWcwWujTHipPFWjc4tTkd2NzZXtV+1NNxo0veHJMH5rvZy63uCK0iDa+Zv4Cs0e9feD6bxch0Au1tzpbd7pSsY5sgHwp4cWzcSGbClrwhFPAgr8DpmQRCAQ+8bU50trfhXCRZ2Jd7bRVGcYZ9ISCOR1u/JQ8Ye1IvXjVQTjG7X2FMvvS+crv7hUqWeK5kI6S0bWbNmtp2ZCqObd3tLbvdKVnHNkG+vxDQzJ9pYWpf8z0xA0iv37Pm3feKk9QKQX4sUkjj2tHa7/tgj6+qlvxYSd7/Pr97QZKh8VjK2JO6ZEx9folnNd319d2EhZbnc7vg97iquiBrpJHpBMfjaU1sE+RDi5Ys0cL3wQwc5n3raS0unoluLlustFVjqxjs9iEcTSGZWX6CXOkOfi6nA70d8wlUzEQ4pS2ubp8bLoes2AUcS2URT+daegKjHSzewdFqqoqRqTgGu7l8jlbPkiAvIn8tIo+LyEERuUdEtqz3OUsDGrvrDYsTshj/ehf8fy0629sQ9LowUpiJbpfWp7kGeXSFLvvFEw0Hgp7ivIRySzgdDjHSpa4QOOwwgdEO+gOeqpc8NsJMPINYKsuWPK2JVS35z6rqc1T1MgDfA/CX631CM6BVu7PaRtC/aBkXsDTYr9Vgjw/DJd31dhhHNivR0RVm2IejqQXbXZZmtCtNaVuqqiAfbf0JjHbQv4oMhY1g9pgxyNNaWBLkVbV0j/oOAOtOL2XOYO7zN26TlGZnTrQD5gOHGYDW2yU81DO/Vn7cJsu+KqXsXWwsksRAyfGGAt75lnysfJAPBTwrjsnPT+hr/QumVtZsqW3NHjMun6O1sCzjhoh8EsDbAcwCeOl6n8+sGDnpbqH+gAeJTA7BduNUm8G4f51BebDHh3ufCSOdzWNyLm2LwGSm7D0WjuHx0RnEkllcsq0TAW8bwtEkDg7PIJNTHA3HsKlz/v0bCBpLFZOZHKbm0uj3L31v+wNePDo8gxMTc4gmjaWHDhFMx9MYi6QgAA4MTxfK0foXTK0sFPQgns7h2HgM2ZzC53aip8ONRCaHMzMJZHKKzZ1edLa3YWoujUgyg/P7/Gh3Gz07iXQOc+ks/B4XPC4HRAS5vGI2kcFMPI1gext6O9wL5m3MxjM4OTmHgNeFbp8bOVUkMzm0tzlxamoOAJjSltakbkFeRH4EoNzmzx9W1e+o6ocBfFhE/gzAewB8tMLzvBPAOwFgaGio4ut1eFzocDsZ5BcJBTxIpHPFCmVzIThtqkGQT2fzeM7H/xuqWNCybVUOh2BbVzv++YGT+OcHTgIwtucNBebH3E1Xnd9b/Hlbtw95Ba7+m3sBlL/Q3NzpxdRcGi/9258sW4aAx1W8ICNrmN+Nl/3dfVU/xuUQPGtzAPFUDicm56Ar9E1624wJmx0eJxKZXLG1XkmXrw0Bb1vV5SEy1a02UdWXV/mn3wTwfVQI8qp6K4BbAWDfvn3LfnVesKsPe7dvvL2Wl3PjFYMLsnddPtSNz9zwHLxkT/+6nvfVF2/CyFQcqooOjwuvvqTc9Vzr+dCr9uCJ07O4aEsnfG4nHhuewcnJOVy0JYjLt3cj4G1DLq84v39+z/Nff+5WtLkEP35mHE+MzuC5Q11Lnve3njeELl8b/B4XgoXKOq+KLp8boYAHCuDcbBLdHW1cC22xX7twAB973YXwtDkR8LoQT+UwFU+jvc2JTZ1euJ0OnJ1NYjaRQU9HGzo8Ljx9NoKDIzPY0tmO1126Bb1+N2KpLJLpHCAChxgTVjvb2zCbyOD0dALT8QzmUlm4nIK3XDmEnf1+xNNZTM9l4HIKvC7jAmA8msIFmwJWvy3UokRXuuSsx4uK7FbVI4Wf3wvgxap6w0qP27dvn+7fv7/u5SMiImoGInJAVfet9fFW9Qt+SkT2AMgDOAXg3RaVg4iIyLYsCfKq+kYrXpeIiGgjsaS7fq1EJArgkNXlqKM+ABNWF6KO7Hx8dj42gMfX6nh8rWuPqq55UkarTeM9tJ6xiWYnIvt5fK3JzscG8PhaHY+vdYnIuiai2SZ3PRERES3EIE9ERGRTrRbkb7W6AHXG42tddj42gMfX6nh8rWtdx9ZSE++IiIioeq3WkiciIqIqtUSQF5FXicghETkqIn9qdXnWS0QGReR/RORXIvKUiLy/cP/HROS0iBws3F5tdVnXSkROisgThePYX7ivR0R+KCJHCv+2ZA5iEdlTco4OikhERD7QyudPRP5JRMIi8mTJfWXPlxg+V/g+Pi4ie60reXUqHN9nReSZwjHcKSJdhft3iEii5Dx+2bKCV6HCsVX8LIrInxXO3SEReaU1pa5eheP7t5JjOykiBwv3t9S5A5aNB7X5/qlqU98AOAEcA3A+ADeAXwK40OpyrfOYNgPYW/g5AOAwgAsBfAzAH1tdvhod40kAfYvu+wyAPy38/KcAPm11OWtwnE4A5wBsb+XzB+AaAHsBPLnS+QLwagA/ACAArgLwC6vLv8bjewUAV+HnT5cc347Sv2v2W4VjK/tZLNQzvwTgAXBeoW51Wn0Mqz2+Rb//OwB/2YrnrlDmSvGgJt+/VmjJXwngqKoeV9U0gH8F8AaLy7QuqnpWVR8t/BwF8DSArdaWqiHeAOBrhZ+/BuDXrStKzbwMwDFVPWV1QdZDVe8HMLXo7krn6w0Avq6GhwB0icjmhhR0jcodn6reo6rZwn8fArCt4QWrgQrnrpI3APhXVU2p6gkAR2HUsU1rueMTYzenNwG4vaGFqqFl4kFNvn+tEOS3Ahgp+f8obBQQRWQHgOcC+EXhrvcUumD+qVW7swsUwD0ickCM7YIBYEBVzxZ+PgdgwJqi1dSbsbCCscv5AyqfLzt+J38XRuvIdJ6IPCYi94nIi6wq1DqV+yza7dy9CMCYFjY8K2jZc7coHtTk+9cKQd62RMQP4A4AH1DVCIAvAdgJ4DIAZ2F0Q7WqF6rqXgDXAfhDEbmm9Jdq9Du19NIOEXEDeD2AbxXustP5W8AO56sSEfkwgCyMba8B49wNqepzAXwQwL+ISNCq8q2RbT+Li7wFCy+yW/bclYkHRev5/rVCkD8NYLDk/9sK97U0EWmDcUK/qarfBgBVHVPVnKrmAdyGJu9GW46qni78GwZwJ4xjGTO7lQr/hq0rYU1cB+BRVR0D7HX+CiqdL9t8J0XkdwC8FsBbCxUpCl3Zk4WfD8AYt77AskKuwTKfRTudOxeA6wH8m3lfq567cvEANfr+tUKQfwTAbhE5r9ByejOAuywu07oUxpH+EcDTqvr3JfeXjqv8BoAnFz+2FYhIh4gEzJ9hTHB6EsZ5e0fhz94B4DvWlLBmFrQi7HL+SlQ6X3cBeHthlu9VAGZLuhVbhoi8CsCHALxeVeMl9/eLiLPw8/kAdgM4bk0p12aZz+JdAN4sIh4ROQ/GsT3c6PLVyMsBPKOqo+YdrXjuKsUD1Or7Z/XMwipnH74axozDYwA+bHV5anA8L4TR9fI4gIOF26sBfAPAE4X77wKw2eqyrvH4zocxg/eXAJ4yzxmAXgD3AjgC4EcAeqwu6zqOsQPAJIDOkvta9vzBuFg5CyADY4zv9yqdLxizer9Y+D4+AWCf1eVf4/EdhTG2aX4Hv1z42zcWPrcHATwK4HVWl38Nx1bxswjgw4VzdwjAdVaXfy3HV7j/nwG8e9HfttS5K5S5UjyoyfePGe+IiIhsqhW664mIiGgNGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5Ik2IBHpLdmp61zJjmUxEfk/VpePiGqDS+iINjgR+RiAmKr+rdVlIaLaYkueiIpE5CUi8r3Czx8Tka+JyE9F5JSIXC8inxGRJ0Tk7kIqTojI5YXNQA6IyH83+450RBsJgzwRLWcngGthbMTzfwH8j6peAiAB4DWFQP95ADeo6uUA/gnAJ60qLBEt5LK6AETU1H6gqhkReQKAE8DdhfufALADwB4AFwP4oZGCG04YKUiJqAkwyBPRclIAoKp5Ecno/CSePIz6QwA8papXW1VAIqqM3fVEtB6HAPSLyNWAsWWmiFxkcZmIqIBBnojWTFXTAG4A8GkR+SWMHbSeb2mhiKiIS+iIiIhsii15IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisqn/Bz6kQtlpzeHXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "states, data = factorial_hmm.sample(jr.PRNGKey(0), 10000)\n",
    "\n",
    "# Plot the states and data\n",
    "fig, axs = plt.subplots(len(num_states) + 1, 1, figsize=(8, 6), sharex=True)\n",
    "for i in range(len(num_states)):\n",
    "    axs[i].imshow(states[i][None, :], cmap=\"Greys\", aspect=\"auto\")\n",
    "    axs[i].set_ylabel(\"State {}\".format(i))\n",
    "    axs[i].set_yticks([])\n",
    "axs[-1].plot(data)\n",
    "axs[-1].set_ylabel(\"Data\")\n",
    "axs[-1].set_xlabel(\"Time\")\n",
    "axs[-1].set_xlim(0, 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee26e4ce-e166-45d8-a559-eef6f54c3950",
   "metadata": {},
   "source": [
    "# Fit the data with another factorial HMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "18ecdf14-ebcd-4704-8a22-ac8135788e68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing...\n",
      "Done.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/collinschlager/Code/ssm-jax-refactor/ssm/factorial_hmm/base.py:28: UserWarning: FactorialHMM.initialize() is not implemented!\n",
      "  warnings.warn(UserWarning(\"FactorialHMM.initialize() is not implemented!\"))\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78d34dadcbcf4cb090a04bac75061721",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAERCAYAAABVU/GxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdUUlEQVR4nO3df7xcdX3n8debexu0CPIrxhCIiRJ/oK6BvaKsroshdIG6Big/Wy2r2LS7YtVWK5Z9bK2rfYBV6foQqamA0VIUgUgqFIQIglZ+3GCAhBQTo9aEkKQREGUrzsxn/zjfuXPuvTOTM3PvuTO59/18POYx5/d8Tgbmc78/jyICMzOzTu3T6wDMzGzv5ARiZmZdcQIxM7OuOIGYmVlXnEDMzKwrTiBmZtaVGZdAJF0paaek9QWPP0vSI5I2SPqHsuMzM9tbaKaNA5H0JuAXwJci4lV7OHYRcC2wJCKekPSCiNg5FXGamfW7GVcCiYi7gJ/lt0l6iaRbJK2VdLekl6ddfwBcFhFPpHOdPMzMkhmXQFpYAbwnIv4j8AHgc2n7S4GXSvqupHskndSzCM3M+sxgrwPoNUnPA/4T8DVJ9c37pvdBYBFwPHA4cJekV0fEk1McpplZ35nxCYSsFPZkRCxusm8rcG9E/Br4kaQfkCWU+6cwPjOzvjTjq7Ai4udkyeFMAGVek3Z/naz0gaRDyaq0tvQgTDOzvjPjEoika4DvAS+TtFXS+cDvAedLehDYACxLh98K7Jb0CHAH8MGI2N2LuM3M+s2M68ZrZmaTY8aVQMzMbHLMqEb0Qw89NBYsWNDrMMzM9ipr1679t4iYPXb7jEogCxYsYHh4uNdhmJntVST9pNl2V2GZmVlXnEDMzKwrTiBmZtYVJxAzM+uKE4iZmXWl5wlE0kmSHpW0WdKFTfbvK+mraf+9khbk9n04bX9U0n+d0sDNzGa4niYQSQPAZcDJwFHAuZKOGnPY+cATEXEkcClwSTr3KOAc4JXAScDn0vXMzGwK9HocyLHA5ojYAiDpK2TzUD2SO2YZ8JG0fB3wWWXzri8DvhIRvyKbDHFzut73ygx46xPPcN3ardRqngLGzPYepx1zOAsP3W9Sr9nrBDIP+GlufSvwulbHRERF0lPAIWn7PWPOnTf2AyQtB5YDzJ8/f8IBXzu8lc+s2UTj0SFmZv3vmBcdNO0SSOkiYgXZEwcZGhqacLHh2UqNWYP78IOPnTzh2MzM9ma9bkTfBhyRWz88bWt6jKRB4PnA7oLnTrpqrcaAix9mZj1PIPcDiyQtlDSLrFF89ZhjVgPnpeUzgG9FNgf9auCc1EtrIdmTAu8rO+BKLRjcxwnEzKynVVipTeMCsgc3DQBXRsQGSR8FhiNiNXAF8OXUSP4zsiRDOu5asgb3CvDuiKiWHXO1FgwMOIGYmfW8DSQibgZuHrPtf+eW/x04s8W5Hwc+XmqAY1RdAjEzA3pfhbXXqdaCAScQMzMnkE5lbSD+ZzMz8y9hh1wCMTPLOIF0yL2wzMwyTiAdqtZqLoGYmeEE0rFK1VVYZmbgBNKxai0Y9DgQMzMnkE5VasGAe2GZmTmBdMoDCc3MMk4gHaq4Ed3MDHAC6ZhLIGZmGSeQDlU8kNDMDHAC6ZhLIGZmGSeQDmXjQPzPZmbmX8IOuQRiZpZxAulQpVbzA6XMzHAC6ZhLIGZmGSeQDrkXlplZxgmkQy6BmJllepJAJB0s6TZJm9L7QU2OWSzpe5I2SHpI0tm5fV+U9CNJ69Jr8VTF7rmwzMwyvfolvBBYExGLgDVpfaxngN+PiFcCJwF/I+nA3P4PRsTi9FpXdsB1LoGYmWV6lUCWASvT8krg1LEHRMQPImJTWn4M2AnMnqoAW6lUPReWmRn0LoHMiYjtaflxYE67gyUdC8wCfpjb/PFUtXWppH3bnLtc0rCk4V27dk04cJdAzMwypSUQSbdLWt/ktSx/XEQEEG2uMxf4MvCOiKilzR8GXg68FjgY+FCr8yNiRUQMRcTQ7NkTL8BUauFxIGZmwGBZF46Ipa32SdohaW5EbE8JYmeL4w4AbgIuioh7cteul15+Jekq4AOTGHpbLoGYmWV6VYW1GjgvLZ8H3Dj2AEmzgFXAlyLiujH75qZ3kbWfrC8z2LqIcC8sM7OkV7+EFwMnStoELE3rSBqS9IV0zFnAm4D/3qS77tWSHgYeBg4FPjYVQddSRZtLIGZmJVZhtRMRu4ETmmwfBt6Vlv8e+PsW5y8pNcAWKrWsCca9sMzMPBK9I9VUBHEJxMzMCaQjlZRAXAIxM3MC6Ui16hKImVmdE0gHRkogA/5nMzPzL2EH3AZiZtbgBNIB98IyM2twAumASyBmZg1OIB1wLywzswYnkA40SiD+ZzMz8y9hBypVl0DMzOqcQDrgNhAzswYnkA6M9MLy80DMzJxAOuESiJlZgxNIB9wLy8yswQmkA+6FZWbW4F/CDrgEYmbW4ATSgWpqRHcbiJmZE0hHPA7EzKyhZwlE0sGSbpO0Kb0f1OK4au6Z6Ktz2xdKulfSZklflTSr7JhH2kDcjdfMrKclkAuBNRGxCFiT1pv5fxGxOL3emtt+CXBpRBwJPAGcX264jTYQV2GZmfU2gSwDVqbllcCpRU+UJGAJcF0353erOtKI7po/M7Ne/hLOiYjtaflxYE6L454jaVjSPZJOTdsOAZ6MiEpa3wrMKy/UzEgvLLkEYmY2WObFJd0OvLDJrovyKxERkqLFZV4UEdskvRj4lqSHgac6iGE5sBxg/vz5RU9rquqpTMzMRpSaQCJiaat9knZImhsR2yXNBXa2uMa29L5F0p3A0cD1wIGSBlMp5HBgW4vzVwArAIaGhlolqULcBmJm1tDLKqzVwHlp+TzgxrEHSDpI0r5p+VDgDcAjERHAHcAZ7c6fbDUPJDQzG9HLBHIxcKKkTcDStI6kIUlfSMe8AhiW9CBZwrg4Ih5J+z4E/ImkzWRtIleUHbBLIGZmDaVWYbUTEbuBE5psHwbelZb/GXh1i/O3AMeWGeNYVZdAzMxGuD9qByqeTNHMbIR/CTvgEoiZWYMTSAfqc2G5DcTMzAmkI9VaDQn2cQIxM9tzApF0yFQEsjeo1MKlDzOzpEgJ5B5JX5N0SpqDasaq1sLtH2ZmSZEE8lKykdxvBzZJ+itJLy03rP6UlUBc62dmBgUSSGRui4hzgT8gG/V9n6RvSzqu9Aj7iEsgZmYNexxImNpA3kZWAtkBvIdsGpLFwNeAhSXG11cqtZrbQMzMkiIj0b8HfBk4NSK25rYPS/rbcsLqTy6BmJk1FKnQ/18R8X/yyUPSmQARcUlpkfWhStW9sMzM6ookkGaPmv3wZAeyN6jWws8CMTNLWlZhSToZOAWYJ+kzuV0HAJXmZ01v7oVlZtbQrg3kMWAYeCuwNrf9aeD9ZQbVr9wGYmbW0DKBRMSDwIOSrs49e3xGcy8sM7OGdlVY10bEWcD3mz2vPCL+Q6mR9SGXQMzMGtpVYb03vb9lKgLZG3guLDOzhnZVWNvT+0+mLpz+5hKImVlDuyqsp4FxVVeAyGY4OaC0qPpUNg7EvbDMzKDNOJCI2D8iDmjy2n+iyUPSwZJuk7QpvR/U5Jg3S1qXe/27pFPTvi9K+lFu3+KJxFOUSyBmZg0tE4ikA9L7wc1eE/zcC4E1EbEIWEOTwYoRcUdELI6IxcAS4Bngm7lDPljfHxHrJhhPIZVajUEPJDQzA9o3ov8DWQP6WrKqrPwvZwAvnsDnLgOOT8srgTuBD7U5/gzgnyLimQl85oS5BGJm1tCuCust6X1hRLw4vddfE0keAHPqjfTA48CcPRx/DnDNmG0fl/SQpEsl7dvqREnLJQ1LGt61a9cEQnYvLDOzvCKz8SLpdOCNZCWPuyPi6wXOuR14YZNdF+VXIiKajTPJXWcu8Grg1tzmD5MlnllkD7v6EPDRZudHxIp0DENDQy0/pwiXQMzMGoo8D+RzwJE0SgB/JOnEiHh3u/MiYmmba+6QNDcitqcEsbPNpc4CVkXEr3PXrpdefiXpKuADe7qPyeC5sMzMGoqUQJYAr4iIAJC0Etgwwc9dTfZkw4vT+41tjj2XMbP/5pKPgFOB9ROMpxCXQMzMGor8Ob0ZmJ9bPyJtm4iLgRMlbQKWpnUkDUn6Qv0gSQvS5317zPlXS3oYeBg4FPjYBOMpxHNhmZk1tBtI+I9kbR77Axsl3ZfWXwfcN5EPjYjdwAlNtg8D78qt/xiY1+S4JRP5/G5Vqy6BmJnVtavC+uSURbGXqNTC40DMzJJ2c2GNrTaa8dwGYmbWsMc2EEmvl3S/pF9IelZSVdLPpyK4fuNeWGZmDUV+DT9L1hNqE/BcsjaKy8oMql+5BGJm1lDoz+mI2AwMREQ1Iq4CTio3rP7kXlhmZg1FxoE8I2kWsE7SJ4DtFEw8041LIGZmDUUSwdvTcRcAvyQbl/E7ZQbVrzwXlplZwx5LIBHxk1QCWQDcADwaEc+WHVi/qdWCCBhwI7qZGVBsLqzfBv4W+CHZlO4LJf1hRPxT2cH1k0otm4fR40DMzDJF2kA+Bbw5NaQj6SXATcCMSiDVlEDcBmJmlilSH/N0PXkkW4CnS4qnb1VqNQC3gZiZJe3mwjo9LQ5Luhm4lmwurDOB+6cgtr7iEoiZ2WjtqrD+W255B/Bf0vIu4DmlRdSnRtpAnEDMzID2c2G9YyoD6XeNEoh7YZmZQbG5sA6XtErSzvS6XtLhUxFcP3EJxMxstCJ/Tl9F9gTBw9LrH9O2GaVadRuImVlekQQyOyKuiohKen0RmF1yXH1npBeWx4GYmQHFEshuSW+TNJBebwN2lx1Yv3EvLDOz0YokkHcCZwGPk02keAYw4xrY3QZiZjZa2wQiaQD4q4h4a0TMjogXRMSpEfGvE/1gSWdK2iCpJmmozXEnSXpU0mZJF+a2L5R0b9r+1TRfV2nqJZB95ARiZgZ7SCARUQVeVNKP83rgdOCuVgekBHYZcDJwFHCupKPS7kuASyPiSOAJ4PwSYhzhubDMzEYrMhfWFuC7klaTTecOQER8eiIfHBEbAdT+L/pjgc0RsSUd+xVgmaSNwBLgd9NxK4GPAJdPJKZ2qqkR3eNAzMwyRRLID9NrH2D/csMZZx7w09z6VuB1wCHAkxFRyW2f1+wCkpYDywHmz5/fdSCVqttAzMzyijwP5C8BJB2QrUbhiRQl3Q68sMmuiyLixsJRTkBErABWAAwNDUW316mGe2GZmeUVeR7IENnAwf3T+lPAOyNi7Z7OjYilE4xvG9kTEOsOT9t2AwdKGkylkPr20lTdC8vMbJQiFfpXAv8zIhZExALg3UzdSPT7gUWpx9Us4BxgdUQEcAdZl2KA84BSSzQVjwMxMxulSAKpRsTd9ZWI+A5QaXN8IZJOk7QVOA64SdKtafthafp4UuniAuBWYCNwbURsSJf4EPAnkjaTtYlcMdGY2qmOtIG4Ed3MDIo1on9b0ueBa8ieB3I2cKekYwAi4oFuPjgiVgGrmmx/DDglt34zcHOT47aQ9dKaEi6BmJmNViSBvCa9/8WY7UeTJZQlkxpRn6p6HIiZ2ShFemG9eSoC6XeVkXEgTiBmZlCsDcRwLywzs7GcQApyG4iZ2WhOIAU1SiD+JzMzg2IDCU9vsvkp4OGI2Dn5IfUnl0DMzEYr0gvrfLKxGnek9eOBtcBCSR+NiC+XFFtfqVbTEwmdQMzMgGIJZBB4RUTsAJA0B/gS2aSGdwEzIoGMlEDcjdfMDCjWBnJEPXkkO9O2nwG/Lies/uNeWGZmoxUpgdwp6RvA19L6GWnbfsCTZQXWb9wGYmY2WpEE8m6yJwe+Ma2vBK5PExrOmEGG7oVlZjZakZHoIek7wLNkU5fcl5LHjFIZeSZ6jwMxM+sTe/xzWtJZwH1kVVdnAfdKOqP9WdNPtVZjcB/t6RG8ZmYzRpEqrIuA19bHfEiaDdwOXFdmYP2mUgu3f5iZ5RSp0N9nzIDB3QXPm1aq1XAPLDOznCIlkFvSw56uSetn0+T5HNOdSyBmZqMVaUT/oKTfAd6QNq1ID4OaUaq1YHBgxhW8zMxaKlICISKuB64vOZa+5hKImdloLROIpKfJuu2O20XWu/eA0qLqQ/VeWGZmlmlZJxMR+0fEAU1e+080eUg6U9IGSTVJQy2OOULSHZIeSce+N7fvI5K2SVqXXqc0u8ZkcgnEzGy0QlVYJVhPNrr9822OqQB/GhEPSNofWCvptoh4JO2/NCI+WXagddWae2GZmeX1JIFExEag7aC8iNgObE/LT0vaCMwDHml5UolcAjEzG22v6FYkaQFwNHBvbvMFkh6SdKWkg9qcu1zSsKThXbt2dR1DNg5kr/jnMjObEqX9Ikq6XdL6Jq9lHV7neWQ9wN4XET9Pmy8HXgIsJiulfKrV+RGxIiKGImJo9uzZ3d0MLoGYmY1VWhVWRCyd6DUk/QZZ8rg6Im7IXXtH7pi/A74x0c/ak2qtxqAfJmVmNqJv62SUNZBcAWyMiE+P2Tc3t3oaWaN8qVwCMTMbrScJRNJpkraSPWv9pjRVCpIOk1SfJuUNwNuBJU26635C0sOSHiJ7Jsn7y47ZvbDMzEbrVS+sVcC46VAi4jHglLT8HbJBi83Of3upATbhEoiZ2Wh9W4XVb7ISiP+5zMzq/ItYkEsgZmajOYEU5LmwzMxGcwIpqFJ1CcTMLM8JpKDseSBOIGZmdU4gBVVrwYAb0c3MRvgXsaCKx4GYmY3iBFJQ1b2wzMxGcQIpqOJeWGZmoziBFOQSiJnZaE4gBbkNxMxsNCeQgqpV98IyM8vzL2JBFY8DMTMbxQmkILeBmJmN5gRSkHthmZmN5gRSQK0W1AL2kROImVmdE0gB1QgAl0DMzHKcQAqo1rIEMuBGdDOzEb16JvqZkjZIqkkaanPcj9Ozz9dJGs5tP1jSbZI2pfeDyoy3UnMJxMxsrF6VQNYDpwN3FTj2zRGxOCLyieZCYE1ELALWpPXSVKupBOJxIGZmI3ryixgRGyPi0QlcYhmwMi2vBE6dcFBtVGo1wCUQM7O8fv+TOoBvSloraXlu+5yI2J6WHwfmtLqApOWShiUN79q1q6sg6o3oHgdiZtYwWNaFJd0OvLDJrosi4saCl3ljRGyT9ALgNkn/EhGjqr0iIiRFqwtExApgBcDQ0FDL49qpug3EzGyc0hJIRCydhGtsS+87Ja0CjiVrN9khaW5EbJc0F9g50c9qp1J1CcTMbKy+rcKStJ+k/evLwG+RNb4DrAbOS8vnAUVLNF0ZKYG4G6+Z2YhedeM9TdJW4DjgJkm3pu2HSbo5HTYH+I6kB4H7gJsi4pa072LgREmbgKVpvTT1brzuhWVm1lBaFVY7EbEKWNVk+2PAKWl5C/CaFufvBk4oM8Y8t4GYmY3nP6kLqHfjdRuImVmDE0gBLoGYmY3nBFJAow3ECcTMrM4JpIBGCcT/XGZmdf5FLMDjQMzMxnMCKcDjQMzMxnMCKcC9sMzMxnMCKcC9sMzMxnMCKcC9sMzMxnMCKcC9sMzMxvMvYgEugZiZjecEUkDVTyQ0MxvHCaQAjwMxMxvPCaQAjwMxMxvPCaQAt4GYmY3nBFKAe2GZmY3nX8QCXAIxMxvPCaQA98IyMxvPCaQAl0DMzMbrSQKRdKakDZJqkoZaHPMySetyr59Lel/a9xFJ23L7Tikz3mrVc2GZmY012KPPXQ+cDny+1QER8SiwGEDSALANWJU75NKI+GSJMY5wCcTMbLyeJJCI2AggFf5BPgH4YUT8pLSg2qjWgoF91Em8ZmbT3t7SBnIOcM2YbRdIekjSlZIOanWipOWShiUN79q1q6sPr6QEYmZmDaUlEEm3S1rf5LWsw+vMAt4KfC23+XLgJWRVXNuBT7U6PyJWRMRQRAzNnj278xsh64Xl9g8zs9FKq8KKiKWTdKmTgQciYkfu2iPLkv4O+MYkfVZTLoGYmY23N1RhncuY6itJc3Orp5E1ypemWguXQMzMxuhVN97TJG0FjgNuknRr2n6YpJtzx+0HnAjcMOYSn5D0sKSHgDcD7y8z3qPmHsCJR80p8yPMzPY6iohexzBlhoaGYnh4uNdhmJntVSStjYhxY/b2hiosMzPrQ04gZmbWFScQMzPrihOImZl1xQnEzMy64gRiZmZdcQIxM7OuOIGYmVlXZtRAQkm7gG6nhD8U+LdJDGdvMRPveybeM8zM+56J9wyd3/eLImLcbLQzKoFMhKThZiMxp7uZeN8z8Z5hZt73TLxnmLz7dhWWmZl1xQnEzMy64gRS3IpeB9AjM/G+Z+I9w8y875l4zzBJ9+02EDMz64pLIGZm1hUnEDMz64oTSAGSTpL0qKTNki7sdTxlkHSEpDskPSJpg6T3pu0HS7pN0qb0flCvY51skgYkfV/SN9L6Qkn3pu/7q5Jm9TrGySbpQEnXSfoXSRslHTfdv2tJ70//ba+XdI2k50zH71rSlZJ2Slqf29b0u1XmM+n+H5J0TCef5QSyB5IGgMuAk4GjgHMlHdXbqEpRAf40Io4CXg+8O93nhcCaiFgErEnr0817gY259UuASyPiSOAJ4PyeRFWu/wvcEhEvB15Ddv/T9ruWNA/4Y2AoIl4FDADnMD2/6y8CJ43Z1uq7PRlYlF7Lgcs7+SAnkD07FtgcEVsi4lngK8CyHsc06SJie0Q8kJafJvtBmUd2ryvTYSuBU3sSYEkkHQ78NvCFtC5gCXBdOmQ63vPzgTcBVwBExLMR8STT/LsGBoHnShoEfhPYzjT8riPiLuBnYza3+m6XAV+KzD3AgZLmFv0sJ5A9mwf8NLe+NW2btiQtAI4G7gXmRMT2tOtxYE6v4irJ3wB/BtTS+iHAkxFRSevT8fteCOwCrkpVd1+QtB/T+LuOiG3AJ4F/JUscTwFrmf7fdV2r73ZCv29OIDaKpOcB1wPvi4if5/dF1ud72vT7lvQWYGdErO11LFNsEDgGuDwijgZ+yZjqqmn4XR9E9tf2QuAwYD/GV/PMCJP53TqB7Nk24Ijc+uFp27Qj6TfIksfVEXFD2ryjXqRN7zt7FV8J3gC8VdKPyaoml5C1DRyYqjlgen7fW4GtEXFvWr+OLKFM5+96KfCjiNgVEb8GbiD7/qf7d13X6rud0O+bE8ie3Q8sSr01ZpE1vK3ucUyTLtX9XwFsjIhP53atBs5Ly+cBN051bGWJiA9HxOERsYDse/1WRPwecAdwRjpsWt0zQEQ8DvxU0svSphOAR5jG3zVZ1dXrJf1m+m+9fs/T+rvOafXdrgZ+P/XGej3wVK6qa488Er0ASaeQ1ZUPAFdGxMd7G9Hkk/RG4G7gYRrtAX9O1g5yLTCfbCr8syJibAPdXk/S8cAHIuItkl5MViI5GPg+8LaI+FUPw5t0khaTdRyYBWwB3kH2B+W0/a4l/SVwNlmPw+8D7yKr759W37Wka4DjyaZs3wH8BfB1mny3KZl+lqw67xngHRExXPiznEDMzKwbrsIyM7OuOIGYmVlXnEDMzKwrTiBmZtYVJxAzM+uKE4hZFyT9c3pfIOl3J/naf97ss8z6jbvxmk1AfvxIB+cM5uZfarb/FxHxvEkIz6xULoGYdUHSL9LixcB/lrQuPW9iQNJfS7o/PV/hD9Pxx0u6W9JqshHQSPq6pLXpGRXL07aLyWaMXSfp6vxnpdHCf52eZ/GwpLNz174z93yPq9MAMbNSDe75EDNr40JyJZCUCJ6KiNdK2hf4rqRvpmOPAV4VET9K6+9Mo4GfC9wv6fqIuFDSBRGxuMlnnQ4sJnt+x6HpnLvSvqOBVwKPAd8lm+fpO5N9s2Z5LoGYTa7fIptbaB3ZNDCHkD2sB+C+XPIA+GNJDwL3kE1ot4j23ghcExHViNgBfBt4be7aWyOiBqwDFkzCvZi15RKI2eQS8J6IuHXUxqyt5Jdj1pcCx0XEM5LuBJ4zgc/Nz99Uxf9v2xRwCcRsYp4G9s+t3wr8jzQ1PpJemh7WNNbzgSdS8ng52WOE635dP3+Mu4GzUzvLbLKnCt43KXdh1gX/lWI2MQ8B1VQV9UWy54ksAB5IDdm7aP6Y1FuAP5K0EXiUrBqrbgXwkKQH0vTydauA44AHyR4I9GcR8XhKQGZTzt14zcysK67CMjOzrjiBmJlZV5xAzMysK04gZmbWFScQMzPrihOImZl1xQnEzMy68v8BevFVZEEr5q4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEGCAYAAACgt3iRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaI0lEQVR4nO3dfbRddX3n8ffHXGlUxADeRTEBko60GFgtOIeHjNihGGchWvNg8KGDwyidlPEJn0ZCnbWEmdUusA7YLh0lS9S0ZXgQRKJUU6Eg6oyRE0gk4TZDCliCAS5WlEgdGvKZP/bv6uFy7slJ9j335Nz7ea111j37tx/Od7uRD3v/9v5t2SYiIqKO5/W7gIiIGHwJk4iIqC1hEhERtSVMIiKitoRJRETUNtTvAvrlpS99qefPn9/vMiIiBsqGDRsetz08vn3Ghsn8+fNpNpv9LiMiYqBI+mG79lzmioiI2hImERFRW8IkIiJqS5hERERtCZOIiKgtYRIREbUlTCIioraESURE1JYwiYiI2voSJpLOkrRF0m5JjZb2kyRtLJ9Nkpa1zDtD0lZJ2yStamlfIGl9ab9W0gFTvT8RETNdv85MNgPLgTvatDdsHw+cAVwhaUjSLODTwOuAhcDbJC0s61wKXG775cBPgHOnoP6IiGjRlzCxPWJ7a5v2p2zvKpOzgbF3Cp8EbLN9v+2ngWuAJZIEnA5cX5ZbAyztafEREfEc+12fiaSTJW0B7gHOK+EyF3ioZbHtpe1Q4ImWABprj4iIKdSzUYMl3QL8eptZH7V900Tr2V4PHCvpFcAaSV+fxJpWAisBjjzyyMnabETEjNezMLG9uOb6I5J2AscBDwNHtMyeV9p+DMyRNFTOTsbaJ9rmamA1QKPR8ETLRUTE3tmvLnOVO7OGyvejgGOAB4E7gaPL/AOAtwJrbRu4DVhRNnEOMOFZT0RE9Ea/bg1eJmk7sAi4WdK6MutUYJOkjcCNwLtsP17OOt4DrANGgOtsbynrXAB8UNI2qj6UK6dwVyIiAlD1H/czT6PRcN60GBGxdyRtsN0Y375fXeaKiIjBlDCJiIjaEiYREVFbwiQiImpLmERERG0Jk4iIqC1hEhERtSVMIiKitoRJRETUljCJiIjaEiYREVFbwiQiImpLmERERG0Jk4iIqC1hEhERtSVMIiKitoRJRETUljCJiIjaEiYREVFbwiQiImrrS5hIOkvSFkm7JTVa2k+StLF8NklaVtqPkHSbpHvLeue3rHOIpG9Kuq/8Pbgf+xQRMZP168xkM7AcuKNNe8P28cAZwBWShoBdwIdsLwROAd4taWFZZxVwq+2jgVvLdERETKG+hIntEdtb27Q/ZXtXmZwNuLTvsH1X+f4kMALMLcstAdaU72uApT0sPSIi2tjv+kwknSxpC3APcF5LuIzNnw+cAKwvTYfZ3lG+PwIc1mHbKyU1JTVHR0cnv/iIiBmqZ2Ei6RZJm9t8lnRaz/Z628cCJwIXSprdss0DgRuA99v+WZt1TTmbmWDbq203bDeGh4f3ed8iIuLZhnq1YduLa64/ImkncBzQlPR8qiC5yvaXWxZ9VNLhtndIOhx4rM7vRkTE3tuvLnNJWlA63JF0FHAM8KAkAVcCI7YvG7faWuCc8v0c4KapqjciIir9ujV4maTtwCLgZknryqxTgU2SNgI3Au+y/TjwKuDtwOkttw6fWda5BHitpPuAxWU6IiKmkKpuhpmn0Wi42Wz2u4yIiIEiaYPtxvj2/eoyV0REDKaESURE1JYwiYiI2hImERFRW8IkIiJqS5hERERtCZOIiKgtYRIREbUlTCIioraESURE1JYwiYiI2hImERFRW8IkIiJqS5hERERtCZOIiKgtYRIREbUlTCIioraESURE1JYwiYiI2hImERFRW1/CRNJZkrZI2i2p0dJ+kqSN5bNJ0rJx682SdLekr7W0LZC0XtI2SddKOmAq9yUiIvp3ZrIZWA7c0aa9Yft44AzgCklDLfPPB0bGrXMpcLntlwM/Ac7tScURETGhvoSJ7RHbW9u0P2V7V5mcDXhsnqR5wOuBz7W0CTgduL40rQGW9qjsiIiYwH7XZyLpZElbgHuA81rC5ZPAR4DdLYsfCjzRssx2YG6Hba+U1JTUHB0dnfziIyJmqJ6FiaRbJG1u81nSaT3b620fC5wIXChptqQ3AI/Z3lCnJturbTdsN4aHh+tsKiIiWgzteZF9Y3txzfVHJO0EjgNeBbxR0plUl78OkvTXwNuBOZKGytnJPODhmqVHRMRe2uOZiaRDp6KQ8lsLxjrcJR0FHAM8aPtC2/NszwfeCvyd7bNtG7gNWFE2cQ5w01TVGxERlW4uc31P0pcknVk6vGuTtEzSdmARcLOkdWXWqcAmSRuBG4F32X58D5u7APigpG1UfShXTkaNERHRPVX/cd9hgSpAFgPvpOrHuA74ou3/2/vyeqfRaLjZbPa7jIiIgSJpg+3G+PY9npm48k3bbwP+E9WlpO9L+pakRT2oNSIiBsweO+BLn8nZVJ3djwLvBdYCxwNfAhb0sL6IiBgA3dzN9X+AvwKW2t7e0t6U9NnelBUREYOkmw74/2r7v7cGiaSzAGxf2rPKIiJiYHQTJqvatF042YVERMTgmvAyl6TXAWcCcyX9Rcusg4Bd7deKiIiZqFOfyY+AJvBGoHUYkyeBD/SyqIiIGCwThontTVQPEF7VMpBiRETEc3S6zHWd7TcDd0t6zpONtn+7p5VFRMTA6HSZ6/zy9w1TUUhERAyuTpe5dpS/P5y6ciIiYhB1usz1JC1vOmydRTXKykE9qyoiIgZKpzOTF09lIRERMbg6nZkcZPtnkg5pN9/2P/WurIiIGCSdOuD/F1Xn+waqy12t7zIx8Bs9rCsiIgZIp8tcbyh/MypwRER01NU74CUtp3oLooFv2/5KL4uKiIjB0s074P8ncB5wD7AZOE/Sp3tdWEREDI5uzkxOB17h8n5fSWuALT2tKiIiBko3Q9BvA45smT6itEVERAAdwkTSVyWtBV4MjEi6XdJtwEhp22eSzpK0RdJuSY2W9pMkbSyfTZKWtcybI+l6SX8vaWTs/fOSDpH0TUn3lb8H16ktIiL2XqfLXJ/o4e9uBpYDV7Rpb9jeJelwqlGLv1pGLf5z4Bu2V0g6AHhhWWcVcKvtSyStKtMX9LD2iIgYp9Otwd/q1Y/aHgGQNL79qZbJ2ZThXCS9BPhd4D+W5Z4Gni7LLQFOK9/XALeTMImImFLd3M11iqQ7Je2U9LSkZyT9rFcFSTpZ0haqu8fOK2clC4BR4AuS7pb0OUkvKqscNjYoJfAIcFiHba+U1JTUHB0d7dUuRETMON10wH8KeBtwH/AC4A+BPd4aLOkWSZvbfJZ0Ws/2etvHAicCF0qaTXUG9UrgM7ZPAH5Om3fTlzvO2g1OOTZ/te2G7cbw8PCediEiIrrU1UOLtrdJmmX7GcrZAXDhHtZZXKcw2yOSdgLHAduB7bbXl9nX86sweVTS4bZ3lH6Wx+r8bkRE7L1uzkyeKh3eGyV9XNIHulxvr0laIGmofD8KOAZ40PYjwEOSfqss+hrg3vJ9LXBO+X4OcFMvaouIiIl1EwpvL8u9h+ry0hHAm+r8qKRlkrYDi4CbJa0rs06luoNrI3Aj8C7bj5d57wWukvQD4HjgT0v7JcBrJd0HLC7TERExhVQebO+8UHVmcgxVf8TWcjfVQGs0Gm42m/0uIyJioEjaYLsxvn2PfSaSXg98FvgHqmHoF0j6I9tfn/wyIyJiEHXTAf8/gN+zvQ1A0r8CbgZmZJhc/NUt3Pujnt0ZHRHRUwtfdhAf+/1jJ3273fSZPDkWJMX9wJOTXklERAysTq/tXV6+NiX9DXAdVZ/JWcCdU1DbfqkXiR4RMeg6Xeb6/ZbvjwL/tnwfpRrqJCIiAug8Ntc7prKQiIgYXN2MzTVP0o2SHiufGyTNm4riIiJiMHTTAf8FqqfMX1Y+Xy1tERERQHdhMmz7C7Z3lc8XgYySGBERv9RNmPxY0tmSZpXP2cCPe11YREQMjm7C5J3Am6neFbIDWAGkcz4iIn6p4xPwkmYBf2r7jVNUT0REDKCOZybl/SVHlYEeIyIi2upmbK77ge9KWks1BD0Ati/rWVURETFQugmTfyif5wEv7m05ERExiPYYJrYvBpB0UDXpDPIYERHP0s0T8A1J9wA/AO6RtEnSv+59aRERMSi6ucz1earX534bQNKpVE/A/3YvC4uIiMHRzXMmz4wFCYDt7wC7eldSREQMmm7OTL4l6Qrgaqr3mbwFuF3SKwFs39XD+iIiYgB0c2byO8BvAh8DLgJeAZxA9TrfT+zLj0o6S9IWSbslNVraT5K0sXw2SVrWMu8DZZ3Nkq6WNLu0L5C0XtI2SdfmmZiIiKnXzd1cv9eD390MLAeuaNPesL1L0uHAJklfBQ4D3gcstP3Pkq4D3gp8EbgUuNz2NZI+C5wLfKYHNUdExAS6OTOZdLZHbG9t0/6U7bH+mNlUl9XGDAEvkDQEvBD4kSQBpwPXl2XWAEt7VnhERLTVlzDpRNLJkrYA9wDnlWHvH6a6pPaPVINN/tT23wKHAk+0BNB2YG6Hba+U1JTUHB0d7e2ORETMID0LE0m3lP6N8Z8lndazvd72scCJwIWSZks6GFgCLKB6QdeLylD4e8X2atsN243h4bySJSJisuyxz0TS8jbNPwXusf3YROvZXlynMNsjknYCx1GFyAO2R0tNXwb+DXAVMEfSUDk7mQc8XOd3IyJi73Vza/C5wCLgtjJ9GrABWCDpv9n+q8kqRtIC4KHSAX8UcAzwIDALOEXSC4F/Bl4DNG1b0m1U71i5BjgHuGmy6omIiO50c5lrCHiF7TfZfhOwkKpj/GTggn35UUnLJG2nCqmbJa0rs06luoNrI3Aj1ZP3j9teT9XJfhdVX8rzgNVlnQuAD0raRtWHcuW+1BQREftOtjsvIN1re2HLtIAtthdKutv2Cb0ushcajYabzWa/y4iIGCiSNthujG/v5jLX7ZK+BnypTK8obS8Cnpi8EiMiYlB1EybvpnrA8NQyvQa4wdUpTS8eaIyIiAHTzRPwlvQd4GmqvpLve0/XxiIiYkbp5n0mbwa+T3V5683Aekkrel1YREQMjm4uc30UOHHsmRJJw8At/GoIk4iImOG6uTX4eeMeTvxxl+tFRMQM0c2ZyTfKcyBXl+m3AH/Tu5IiImLQdNMB/18kvQl4VWlabfvG3pYVERGDpJszE2zfANzQ41oiImJATRgmkp7k2e8T+eUsqjuGD+pZVRERMVAmDBPbL57KQiIiYnDlrqyIiKgtYRIREbUlTCIioraESURE1JYwiYiI2hImERFRW8IkIiJqS5hERERtCZOIiKitL2Ei6SxJWyTtlvScF9NLOlLSTkkfbmk7Q9JWSdskrWppXyBpfWm/VtIBU7UfERFR6deZyWaq98rfMcH8y4Cvj01ImgV8GngdsBB4m6SFZfalwOW2Xw78BDi3V0VHRER7fQkT2yO2t7abJ2kp8ACwpaX5JGCb7fttPw1cAyyRJOB0fvXWxzXA0l7VHRER7e1XfSaSDgQuAC4eN2su8FDL9PbSdijwhO1d49on2v5KSU1JzdHR0ckrPCJihutZmEi6RdLmNp8lHVa7iOqS1c5e1GR7te2G7cbw8HAvfiIiYkbq6uVY+8L24n1Y7WRghaSPA3OA3ZJ+AWwAjmhZbh7wMNX76OdIGipnJ2PtERExhXoWJvvC9qvHvku6CNhp+1OShoCjJS2gCou3An9g25JuA1ZQ9aOcA9w09ZVHRMxs/bo1eJmk7cAi4GZJ6zotX8463gOsA0aA62yPddBfAHxQ0jaqPpQre1d5RES0I7vdm3mnv0aj4Waz2e8yIiIGiqQNtp/zfOB+dTdXREQMpoRJRETUljCJiIjaEiYREVFbwiQiImpLmERERG0Jk4iIqC1hEhERtSVMIiKitoRJRETUljCJiIjaEiYREVFbwiQiImpLmERERG0Jk4iIqC1hEhERtSVMIiKitoRJRETUljCJiIjaEiYREVFbX8JE0lmStkjaLek5L6aXdKSknZI+XKaPkHSbpHvLeue3LHuIpG9Kuq/8PXgq9yUiIvp3ZrIZWA7cMcH8y4Cvt0zvAj5keyFwCvBuSQvLvFXArbaPBm4t0xERMYX6Eia2R2xvbTdP0lLgAWBLy/I7bN9Vvj8JjABzy+wlwJryfQ2wtDdVR0TERParPhNJBwIXABd3WGY+cAKwvjQdZntH+f4IcFiHdVdKakpqjo6OTk7RERHRuzCRdIukzW0+SzqsdhFwue2dE2zzQOAG4P22fzZ+vm0DnmjjtlfbbthuDA8P790ORUTEhIZ6tWHbi/dhtZOBFZI+DswBdkv6he1PSXo+VZBcZfvLLes8Kulw2zskHQ48Vrv4iIjYKz0Lk31h+9Vj3yVdBOwsQSLgSmDE9mXjVlsLnANcUv7eNEXlRkRE0a9bg5dJ2g4sAm6WtG4Pq7wKeDtwuqSN5XNmmXcJ8FpJ9wGLy3REREwhVd0MM0+j0XCz2ex3GRERA0XSBtvPeT5wv7qbKyIiBlPCJCIiakuYREREbQmTiIioLWESERG1JUwiIqK2hElERNSWMImIiNoSJhERUVvCJCIiakuYREREbQmTiIioLWESERG1JUwiIqK2hElERNSWMImIiNoSJhERUVvCJCIiakuYREREbQmTiIiorS9hIuksSVsk7Zb0nBfTSzpS0k5JHx7XPkvS3ZK+1tK2QNJ6SdskXSvpgKnYh4iI+JV+nZlsBpYDd0ww/zLg623azwdGxrVdClxu++XAT4BzJ6vIiIjoTl/CxPaI7a3t5klaCjwAbBnXPg94PfC5ljYBpwPXl6Y1wNLJrzgiIjrZr/pMJB0IXABc3Gb2J4GPALtb2g4FnrC9q0xvB+Z22P5KSU1JzdHR0ckpOiIiGOrVhiXdAvx6m1kftX3TBKtdRHXJamd10vHLbb0BeMz2Bkmn7WtNtlcDq8s2RyX9cC9Wfynw+L7+9gDLfs8s2e+ZZ2/3/ah2jT0LE9uL92G1k4EVkj4OzAF2S/oF1dnGGyWdCcwGDpL018DbgTmShsrZyTzg4S7rG96bwiQ1bT/nZoHpLvs9s2S/Z57J2veehcm+sP3qse+SLgJ22v5UabqwtJ8GfNj22WX6NmAFcA1wDjDRWU9ERPRIv24NXiZpO7AIuFnSuhqbuwD4oKRtVH0oV05GjRER0b2+nJnYvhG4cQ/LXDRB++3A7S3T9wMnTV51E1o9Bb+xP8p+zyzZ75lnUvZdtidjOxERMYPtV7cGR0TEYEqYREREbQmTPZB0hqStZeyvVf2up1ckHSHpNkn3lnHTzi/th0j6pqT7yt+D+11rL4wf922mjPkmaY6k6yX9vaQRSYtmwjGX9IHyz/lmSVdLmj0dj7mkz0t6TNLmlra2x1eVvyj7/wNJr9yb30qYdCBpFvBp4HXAQuBtkhb2t6qe2QV8yPZC4BTg3WVfVwG32j4auLVMT0fjx32bKWO+/TnwDdvHAL9D9b/BtD7mkuYC7wMato8DZgFvZXoe8y8CZ4xrm+j4vg44unxWAp/Zmx9KmHR2ErDN9v22n6Z6lmVJn2vqCds7bN9Vvj9J9S+VuVT7u6YsNi3HPhs/7ttMGfNN0kuA36XcTm/7adtPMAOOOdWdrC+QNAS8ENjBNDzmtu8A/mlc80THdwnwl658j+qB8MO7/a2ESWdzgYdapjuO/TVdSJoPnACsBw6zvaPMegQ4rF919dAnefa4b3s15tsAWwCMAl8ol/g+J+lFTPNjbvth4BPAP1KFyE+BDcyMYw4TH99a/75LmMSzlME2bwDeb/tnrfNc3Uc+re4lbx33rd+19MEQ8ErgM7ZPAH7OuEta0/SYH0z1X+ELgJcBL+K5l4JmhMk8vgmTzh4GjmiZ7nrsr0Ek6flUQXKV7S+X5kfHTnXL38f6VV+PvIpq3LcHqS5jnk7VjzCnXAKB6XvctwPbba8v09dThct0P+aLgQdsj9r+F+DLVP8czIRjDhMf31r/vkuYdHYncHS5y+MAqk66tX2uqSdKP8GVwIjty1pmraUa8wym4dhnti+0Pc/2fKrj+3e2/z0wNuYbTMP9BrD9CPCQpN8qTa8B7mWaH3Oqy1unSHph+ed+bL+n/TEvJjq+a4H/UO7qOgX4acvlsD3KE/B7UEYq/iTVHR+ft/0n/a2oNySdCnwbuIdf9R38MVW/yXXAkcAPgTfbHt+hNy20DCL6Bkm/QXWmcghwN3C27f/Xx/J6QtLxVDceHADcD7yD6j8yp/Uxl3Qx8BaquxjvBv6Qqn9gWh1zSVcDp1ENM/8o8DHgK7Q5viVYP0V1ye8p4B22m13/VsIkIiLqymWuiIioLWESERG1JUwiIqK2hElERNSWMImIiNoSJhE1Sfrf5e98SX8wydv+43a/FbG/ya3BEZOk9TmVvVhnqGU8qHbzd9o+cBLKi+ipnJlE1CRpZ/l6CfBqSRvL+zJmSfozSXeW90P8UVn+NEnflrSW6slrJH1F0obyjo2Vpe0SqpFtN0q6qvW3ylPKf1bex3GPpLe0bPv2lneUXFUeRovoqaE9LxIRXVpFy5lJCYWf2j5R0q8B35X0t2XZVwLH2X6gTL+zPIX8AuBOSTfYXiXpPbaPb/Nby4Hjqd5B8tKyzh1l3gnAscCPgO9SjTv1ncne2YhWOTOJ6J1/RzXW0UaqYWkOpXrxEMD3W4IE4H2SNgHfoxps72g6OxW42vYzth8FvgWc2LLt7bZ3AxuB+ZOwLxEd5cwkoncEvNf2umc1Vn0rPx83vRhYZPspSbcDs2v8but4Us+Q/5/HFMiZScTkeRJ4ccv0OuA/l6H9kfSb5eVT470E+EkJkmOoXps85l/G1h/n28BbSr/MMNUbE78/KXsRsQ/yXywRk+cHwDPlctUXqd6LMh+4q3SCj9L+VbDfAM6TNAJspbrUNWY18ANJd5Wh8cfcCCwCNlG93Ogjth8pYRQx5XJrcERE1JbLXBERUVvCJCIiakuYREREbQmTiIioLWESERG1JUwiIqK2hElERNT2/wHvdUmNKu9fQgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_factorial_hmm = NormalFactorialHMM(num_states, seed=jr.PRNGKey(1))  # <-- different seed\n",
    "lps, _, posteriors = test_factorial_hmm.fit(data, tol=-1)\n",
    "\n",
    "plt.plot(lps)\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"log probability\")\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(np.arange(1, len(lps)), lps[1:])\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"log probability\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d0853d48-fd6e-45be-842c-24a3907519f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFzCAYAAAA9nXBaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABKPElEQVR4nO3deXhjZ3k3/u8tyZK8SN5GsmexZyazkQUSJhNIWAIEyg5pISWkaUNTWpaUrfwKXaAFWrga6HYFyBtI3pYmvDRtaUgJFNJAaBIgCclMEiYJyewzXmbG8m7J2qX798fRkWVbkmVb1pGOv5/r0jUeSUd6jo507vNs9yOqCiIiIrIfh9UFICIiorXBIE9ERGRTDPJEREQ2xSBPRERkUwzyRERENsUgT0REZFMuqwuwHCJScr5fV1cXtm7dimw2i+PHjyMcDteyaHXJ4/Fg165dcLvdGB4exsjICDZv3oze3l5Eo1EcOXIETU1N2LlzJ9xud3670dFRDAwMlH3tbdu2obu7O///RCKBo0ePIh6Pr6isZrlM6XQax44dQyQSWdHrrYT5eXk8niWfq6oYHBzE6Ogo+vv7EQgE8o8lk0kcPXoUsVgsf9/Cz6ucqakpHD9+HCuZ3up0OrFz5060tbXlj2MgEEBfXx9EBACQzWYxMDCA8fHxRdu3t7dj+/btcDqdVS1XtYkIduzYgfb29qq9ZigUwuDg4KLPq15EIhEcPXoUmUxmTV7f4XBg586d8Pl8q3qdbDaLU6dOYWJiYtFj5b5fplgshiNHjiCVSq2qHGvF6/Xmz6umsbExnDp1as3eU1VX/mVU1Ya5AdBSt2uvvVYjkYiOjo7qFVdcUfJ56+m2a9cuPXr0qCaTSf3EJz6hAPTGG2/UTCajTz75pPb29uqFF16oAwMDaspms3rLLbeow+Eo+bpOp1Nvv/12LXT8+HE999xzV1xWs1ymUCikl19+eU0/r927d+uxY8e0ErFYTD/4wQ+qiOjNN9+s2Ww2/9jg4KBedNFF8z6vO+64o6LXVVW9++671e12r2gfOjs79ec//7lms1m9+eabVUT0hhtu0Fgsln/92dlZve6664pu/5a3vEUnJyerXq5q37xer95zzz0Vf6ZLyWaz+uUvf1lFRD/0oQ9pPB6v2mtXy0MPPaTt7e1r9pn6fD594IEHVl3OSCSi11xzTdH3uPLKK3Vqaqrs9gcPHtTNmzdb/h0rdbvgggv05MmT+fJms1m97bbbyp4zV3vTVcRNNtcTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdlUzYO8iLSLyI0i8ryITIjIuIg8l7uvo9blISIisisravL/AWASwKtVtUtVuwG8Jnfff1hQHiIiIluyIshvU9UvqupZ8w5VPauqXwSw1YLyEBER2ZIVQf6UiHxSRHrMO0SkR0T+BMCgBeUhIiKyJSuC/NUAugE8mOuTnwDwAIAuAO+yoDxERES25Kr1G6rqJIA/yd2IiIhojXAKHRERkU0xyBMREdkUgzwREZFNWRbkRaRFRP5CRG7L/X+XiLzVqvIQERHZjZU1+W8ASAC4LPf/YQCft644RERE9mJlkN+hql8CkAIAVY0CEAvLQ0REZCtWBvmkiDQDUAAQkR0wavZERERUBTWfJ1/gswDuBdAnIt8C8HIA11tYHiIiIluxLMir6n0icgDApTCa6T+qqmNWlYeIiMhurBxdf7+qjqvqf6vq91V1TETut6o8REREdlPzmryIeAG0ANggIp2YG2znB7C51uUhIiKyKyua698P4GMANgE4gLkgPwPgqxaUh4iIyJasWKDmJgA3iciHVfUrtX5/IiKi9cLKgXdfEZELAJwHwFtw/x1WlYmIiMhOLAvyIvIZAK+GEeR/AOBNAH4GgEGeiIioCqxMhnMVgNcCOKuq1wO4EEC7heUhIiKyFSuDfExVswDSIuIHEALQZ2F5iIiIbMXKjHf7RaQDwG0wRtlHADxiYXmIiIhsxcqBdzfk/vyaiNwLwK+qB60qDxERkd1YmvHO/FtVT6rqQWa8IyIiqh5mvCMiIrIpZrwjIiKyKWa8IyIisqma98mLyCUi0msGeBG5TkS+KyJfFpGuWpeHiIjIrqwYePd1AEkAEJHLAdwII8vdNIBbLSgPERGRLVnRJ+9U1Ync31cDuFVV7wJwl4g8ZUF5iIiIbMmKmrxTRMyLi9cC+EnBY1Ym5yEiIrIVK4LqnQAeFJExADEAPwUAEdkJo8meiIiIqsCK0fVfyCW92QjgPlXV3EMOAB+udXmIiIjsypLmcVV9tMh9h60oCxERkV1ZuQodERERrSEGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbElW1ugwVE5EwgENWl2MNbQAwZnUh1pCd98/O+wZw/xod969x7VFV30o3dlWzJDVwSFX3WV2ItSIi+7l/jcnO+wZw/xod969xicj+1WzP5noiIiKbYpAnIiKyqUYL8rdaXYA1xv1rXHbeN4D71+i4f41rVfvWUAPviIiIqHKNVpMnIiKiCjHIExER2RSDPBERkU0xyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2VRDLVAjIszcU2XNzc3YuXMn3G53/r7R0VEMDAzUtBybN29Gb29v/v/pdBrHjh1DJBKpaTkandPpxM6dO9HW1pY/joFAAH19fRARAEA2m8XAwADGx8cXbd/e3o7t27fD6XQuemxqagrHjx9HPSTQEhHs2LED7e3tVXvNUCiEwcHBqr1eo3E4HNi5cyd8vhUveAbA+H6dOnUKExMTix4r9/0yxWIxHDlyBKlUalXlWCterxe7du2ad84cGxvDqVOn1uw9VVVWs3HD3AAob9W9XXjhhTowMKCmbDart9xyizocjpqW48Ybb9RMJpMvRygU0ssvv9zyz6fRbp2dnfrzn/9cs9ms3nzzzSoiesMNN2gsFst/trOzs3rdddcV3f4tb3mLTk5OajF33323ut1uy/cRgHq9Xr3nnnuKlnMlstmsfvnLX9ZcRWJd3nw+nz7wwAOr/iwjkYhec801Rd/jyiuv1KmpqbLbHzx4UDdv3mz551HqdsEFF+jJkyfnfXduu+22NT1n6iriJpvriYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIpuoqyIvID60uAxERkV24av2GIrK31EMALqphUYiIiGyt5kEewOMAHoQR1BfqqG1RiIiI7MuKIP8cgPer6pGFD4jIoAXlISIisiUr+uQ/W+Z9P1zDchAREdlazWvyqvqfZR77rxoWhYiIyNbqanQ9ERERVQ+DPBERkU0xyBMREdmUZUFeRFpE5C9E5Lbc/3eJyFutKg8REZHdWFmT/waABIDLcv8fBvB564pDRERkL1YG+R2q+iUAKQBQ1SiKJ8ghIiKiFbAyyCdFpBmAAoCI7IBRsyciIqIqsCLjnemzAO4F0Cci3wLwcgDXW1geIiIiW7EsyKvqfSJyAMClMJrpP6qqY1aVh4iIyG6sHF1/v6qOq+p/q+r3VXVMRO63qjxERER2Y8VSs14ALQA2iEgn5gbb+QFsrnV5iIiI7MqK5vr3A/gYgE0ADmAuyM8A+KoF5SEiIrIlKxaouQnATSLyYVX9Sq3fn4iIaL2wcuDdV0TkAgDnAfAW3H+HVWUiIiKyE8uCvIh8BsCrYQT5HwB4E4CfAWCQJyIiqgIrk+FcBeC1AM6q6vUALgTQbmF5iIiIbMXKIB9T1SyAtIj4AYQA9FlYHiIiIluxMuPdfhHpAHAbjFH2EQCPWFgeIiIiW7Fy4N0NuT+/JiL3AvCr6kGrykNERGQ3lma8M/9W1ZOqepAZ74iIiKqHGe+IiIhsihnviIiIbIoZ74iIiGyq5n3yInKJiPSaAV5ErhOR74rIl0Wkq9blISIisisrBt59HUASAETkcgA3wshyNw3gVgvKQ0REZEtW9Mk7VXUi9/fVAG5V1bsA3CUiT1lQHiIiIluyoibvFBHz4uK1AH5S8JiVyXmIiIhsxYqgeieAB0VkDEAMwE8BQER2wmiyJyIioiqwYnT9F3JJbzYCuE9VNfeQA8CHa10eIiIiu7KkeVxVHy1y32ErykJERGRXVq5CR0RERGuIQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsikGeiIjIphjkiYiIbIpBnoiIyKYY5ImIiGyKQZ6IiMimGOSJiIhsSlTV6jJUTETCAA5ZXY41tAHAmNWFWEN23j877xvA/Wt03L/GtUdVfSvd2FXNktTAIVXdZ3Uh1oqI7Of+NSY77xvA/Wt03L/GJSL7V7M9m+uJiIhsikGeiIjIphotyN9qdQHWGPevcdl53wDuX6Pj/jWuVe1bQw28IyIioso1Wk2eiIiIKsQgT0REZFMM8kRERDbFIE9ERGRTDPJEREQ2xSBPRERkUwzyRERENsUgT0REZFMNtUDNhg0bdNu2bVBVDA8PY2RkJP9YV1cXtm7dimw2i+PHjyMajWLnzp1oa2tb9DqpVApHjx5FNBpd8j2bmpqwa9cuNDc3V3VfyslmsxgYGMD4+HhF5Tpz5gxOnz696Dkejwe7du2C2+3Of16bN29Gb28votEojhw5gnQ6vWi7vr4+BIPB/P+LfV7btm1Dd3d3/v+JRAJHjx5FPB4vWebOzk5s3boVTqdz0WMTExM4ceJE/v/Nzc3YuXMn3G73oufOzs7iyJEjyGQyJd9rJczPy+PxLPlcVcXg4CBGR0fR39+PQCCQfyyZTOLo0aOIxWL5+xZ+XuVMTU3h+PHjqCRRVXd3N/r7++FwOEqWKxKJ4OjRo/nPy+VyYdeuXWhpaVn0eolEAkeOHEEikVj02MaNG7Fp06aK9mEhVcXAwADGxla3UJjb7c4fozNnzuDMmTOrer1iAoEA+vr6ICJVf+3VWHgcq31uUlUMDQ0hFAqt6nXKlavc98tknssLv9MLy7Vhwwb09fXln5PNZjE4OLjq71c5C2NMOBzOP7bwnFlNJ0+exNjY2Mq/jKraMLeLL75YVVWTyaR+8pOfVAD527XXXquRSERHR0f1iiuu0Pb2dn3ooYe0mOHhYX3JS14yb/tSt76+Pn366aeLvs5aiUajev311y9ZrmeeeUbT6bR+7nOfK/qcXbt26dGjRzWZTOonPvEJBaA33nijZjIZffLJJ7W3t3fRNiKiX/3qVzWbzebLMzQ0pPv27cs/x+l06u233z6vzMePH9dzzz23bJmvuuoqnZ6eLrrP3/rWt9TlcuWfe9FFF+nAwEDR5z7yyCPa2dlZ0fFbzm337t167Nixio5RLBbTD37wgyoievPNN8/7vAYHB/Wiiy6a93ndcccdFb2uqurdd9+tbre7ojK/5z3v0dnZ2bLl+tnPfjbv8+rt7dUnnnii6HsfPXpUd+/eXfS9PvOZz2g6na54PxZ+Xu973/tWfYy2bdumzz//vKZSKf30pz9d9e8AAP3Qhz6k8Xh8Rfu5lh566CFtb2/Pl3Pjxo361FNPVe31E4mEfvSjH13157dp0yY9ePBg0fc4fPiw7tixo+z211xzjYbD4bLleu9736vRaDT/nGg0qu9973vX5PtQWK5IJKJjY2N6xRVX5O8XEb3pppvmnQOqKRf3Vhw32VxPRERkUwzyRERENsUgT0REZFMM8kRERDbFIE9ERGRTDPJEREQ2xSBPRERkUwzyRERENsUgT0REZFMM8kRERDbFIE9ERGRTDPJEREQ2xSBPRERkUwzyRERENsUgT0REZFMM8kRERDbFIE9ERGRTDPJEREQ2VfMgLyJ+EfkbEfmmiPzWgsf+T63LQ0REZFdW1OS/AUAA3AXg3SJyl4h4co9dakF5iIiIbMmKIL9DVf9UVf9LVd8O4AkAPxGRbgvKQkREZFsuC97TIyIOVc0CgKp+QUSGATwEoM2C8hAREdmSFTX57wG4ovAOVf0XAP8fgKQF5SEiIrKlmtfkVfWTJe6/F8CuGheHiIjItjiFjoiIyKYY5ImIiGyKQZ6IiMimLAvyItIiIn8hIrfl/r9LRN5qVXmIiIjsxsqa/DcAJABclvv/MIDPW1ccIiIie7EyyO9Q1S8BSAGAqkZhZMIjIiKiKrAyyCdFpBmAAoCI7IBRsyciIqIqsCLjnemzAO4F0Cci3wLwcgDXW1geIiIiW7EsyKvqfSJyAMaiNALgo6o6ZlV5iIiI7MbK0fX3q+q4qv63qn5fVcdE5H6rykNERGQ3Na/Ji4gXQAuADSLSibnBdn4Am2tdHiIiIruyorn+/QA+BmATgAOYC/IzAL5qQXmIiIhsyYoFam4CcJOIfFhVv1Lr9yciIlovrBx49xURuQDAeQC8BfffYVWZiIiI7MSyIC8inwHwahhB/gcA3gTgZwAY5ImIiKrAymQ4VwF4LYCzqno9gAsBtFtYHiIiIluxMsjHVDULIC0ifgAhAH0WloeIiMhWrMx4t19EOgDcBmOUfQTAIxaWh4iIyFasHHh3Q+7Pr4nIvQD8qnrQqvIQERHZjaUZ78y/VfWkqh5kxjsiIqLqYcY7IiIim2LGOyIiIptixjsiIiKbqnmfvIhcIiK9ZoAXketE5Lsi8mUR6ap1eYiIiOzKioF3XweQBAARuRzAjTCy3E0DuNWC8hAREdmSFX3yTlWdyP19NYBbVfUuAHeJyFMWlIeIiMiWrKjJO0XEvLh4LYCfFDxmZXIeIiIiW7EiqN4J4EERGQMQA/BTABCRnTCa7ImIiKgKrBhd/4Vc0puNAO5TVc095ADw4VqXh4iIyK4saR5X1UeL3HfYirIQERHZlZWr0BEREdEaYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIim2KQJyIisikGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5IiIimxJVtboMFRORMIBDVpdjDW0AMGZ1IdaQnffPzvsGcP8aHfevce1RVd9KN3ZVsyQ1cEhV91ldiLUiIvu5f43JzvsGcP8aHfevcYnI/tVsz+Z6IiIim2KQJyIisqlGC/K3Wl2ANcb9a1x23jeA+9fouH+Na1X71lAD74iIiKhyjVaTJyIiogoxyBMREdkUgzwREZFNMcgTERHZFIM8ERGRTTHIExER2RSDPBERkU0xyBMREdlUQy1Qs2HDBt22bZvVxSAiIqqJAwcOjKlqYKXbN1SQ37ZtG/bvX9WCPERERA1DRE6tZns21xPRIuF4Cg8fG8PgRNTqolQdU3kvTVWRymSXfF46k0UsmUEinalBqWglGqomT2vj7ieH8PjJSVy+K4BX7NqANk/tvxbHRiO4+uuPYndPG16xawOufelWtDc31bwcdvDNR09hNJzAx39tNwDjhH1waBpnpuOIJtN4zZ4gOlvdRbdNZ7K45rZHsf/UJFQBp0Nw9SV9+Ohrd6HH763lbizpH390GLc8eAwdzU3obHGjo6UJHS1NcLuccDkEE7NJnJ6KocnpwHmb/PB7m/DYyXE8dyaMHp8H2za0YjaRxvBUHL/10v7857WefXv/IL7zxDCePT2NmXgabR4XAj4Pzt/kxwWb25FVxVQ0hSMjYTw9PIOxSAIA4BDgnEAbLurrwF+85Ty0txT/7UaTafzfn57AxnYvdgbbcHBoGj9+bgShmQREgEu2deGvf/2CWu7yssWSGXzqv55GoM2Di/o68JoXBOFtclpdrJIY5Am3P3wKTw1O4V9/MYBevxeP/NkVEJGaluHJgSmMRRLwe1340r2HIBB88NU7aloGu/jeL09jeDKWD1rPnp7BlTf/PP/4R67YiY+/fk/Rbc9Mx/H4yUm8/cJN+I0Xb8YDh0L418cGcODkJP7njy6vSfkrtf/UBLpa3HjV7gAmo0lMRVM4ORZFMpNFKpNFZ4sb2ze0Ip7O4oFDIczE09jb34H3vmI7RsMJnByfhb+5CUOTMew/OWH17tSFf/zRYWQVeOuFm9Dj82IqlsSZqTieHJjC9w+eAQC4XQ5s627Bq3YH0N/VArfLgVgyjScHp/CfB4bw2hcE8aYXbiz6+g8eGsU//OjwvPt2BFqxM9iGIyMR/Pv+QfzVlefX/PyzHE8PT+M7TwxDBFAF/vA1O/CJN7zA6mKVxCBPiCTS+LXzenBOoBVff/A4ZmLpklfia2VoMgoR4N6PXY59n/8RzkzHavr+djIaTmA0nICqQkQwkGtyv+Xavfjs957F8FS85LahsFEz+429m/GaPUG85gVBdLS4cdP9RxBPZeqqxhKaSeCivg588aoXLflcVUU21zKx0PXfeAxjkeRaFLGhZLOK0UgCv//Kc/Anb1wctKZjKXhcjpLfgbFIAvs+/+P8d6iYkRnju/efH7gMo+EEdvX4sDPYBgC47aHj+MIPnrPk/LMcobCxD/f84SvwwW8dwOBEfZ+r2CdPiMTT6Gpx4/xN7QDmvsS1NDgRQ6/fC7fLgR6/N38yoOULzcSRzGQxE0sDmDuxvmR7Fza2N5c9vqO5x4I+T/6+rd0tAIDhqfo6mYXCCQQKylmOiBQN8ADQ5m1CJJGuZtEa0mQ0iVRG5x37Qu3NTWUv8rpa3HA5pOxvdyScQJNTsLe/E2964cZ8gAeAoN94XyvOP8sxMmNcxGzpbMbGdm/dl5dBnhCOp9DmdaEn9+M2v8S1NDQZxZbOZgDIBfnal8EOIok0ZpPGICjz5BPKnVg7W9zo8XsQKvPZmrWwoG+u/72/ywjyA3U0CC+eymA6lioZkJajzeNCOM4gX+zYL4fDIQj4PGVr8qGZBAJtHjiKXHCZF2zltq8HoXAcbqcDHS1NCPq8ZX9P9YBBfp3LZBWzyQzaPC4EcwOrrKhFD03G0NdpBJOg34PROv+h16tQwbEzT5YjM/H8iTXo82KkTM0jNJOA0yHoLhiY15cL8kN1FOTNAV9m7W81/F4XwvHUql+n0eWD/Co+0+BSQT4cz59nFm/rzT+nnoVmEgj6PRARBP3l97ceMMivc7NJowbj87rytaJaf2nTmSzOzsTn1eRD4TiyWU51Wq7CY2eeLEfDifyJtcfvwVQ0hXiq+JSnUDiODW3ueTWtQJsHHpejrmryq611FmrzuJBIZ5FMLz1lzM7MC8SeVXymQb933oXm4vdIlGx96TGb6+u8ZjwyE8/vQ9DnNVrP6ri7x5IgLyK/KSLPikhWRPZZUQYyROJzQb7V44LP46p5Tf7MdByZrGJLribf4/MglVFMRjkYarnmBfmZuZp8/qSUC/alWkpGi/RzOxyCLZ3NdTXAyNy3Svvky2nzGuOP6/lEXQu1qMmPhOMlp2K2eVxobnLWfc04FE7k9yF/YVLHZbaqJv8MgHcAeMii96ccsy+yzWOMZg1Y0FQ+OGnUEAtr8oA1YwManVmLcsjciafwpDTXWlP8Qi4UThStHfd3teSPUz0oNkBwpXxe47u/3vvlQzNx+LyuVc2gCPq8mJhNFm0ViacymIqm8oFxoUZp/l5YkwdQtvXCapYEeVV9TlUPWfHeNF8kYfRFmrWZHl/tR7YPTRo1RLPvNz82oM775urRaDgBt8uBvq4WhMKJRSfWpS6gjCC/+CTc19VSd831DgG626oz8A4Awon13S9f6tgvh/k9G40s/n6NVtDFEvR56jpgxpIZhOPped1fgDFroF6xT36dm6vJGye6oN9T8+A6NBmDQ4De9gVNYHX8Y69X5onaPFkuPLH2lBlcmckqxiMlgnxnC8LxNKaj9REIR8MJdLd5Sk6LWw5f7gI3ss5r8iMzpZvSKxUs89s1W4/KdQcEfd66HnRr7sNcy9g6rsmLyI9F5JkityuX+TrvE5H9IrJ/dHR0rYq7bplB3m/W5P3GlJBa5vcemohiY3szmpzG1zFg4VS+Rmc2JQZ9XoxGEotOrJ0tTWhyStHPdjySQFaL93P31dk0umrUOk1mkF/3zfVV+EzNoFfs+2WOoyhXk19qCp7VRvL7YHxO/mYX3C5HXV+YrFnGO1V9XZVe51YAtwLAvn37ONy6yswkIGZzfdDnQSKdrWnWqaHJGDbn+uMBwONyorOlqe6n0tSjUDiBnYE2BHwePHQ4UXBSMk6sIsY0umKfrXlyDRQ5Cfd1GcdncDKKF25pX6viVywUjldl0B0w14q1nhPiqKoR5KtUkx8t8v0yW49K9cmb20cSaUSTabS46y8h68KavIigx++p6+RdbK5f5yKLmutrP1d1aDKanyNvYkKclQnNxBH0exDweRBOpHFq3Kh5F55YgyUS4oyWGV1t1uTrZVW6clOxlsu8wA2v4yA/E0sjmc6u+jPtbvXAIcVr8ma2u86W4osjAYXN3/X52zf3a97vyeet69YHq6bQ/YaIDAG4DMB/i8j/WFEOMk5sIkCr2xx4V9um8mR6/hx501LzbWmxeCqDmXg63ycPAM+cnobLMf/EGvQVr3mEyoxY93uNFd7qobk+k1WMzyarMkceMPYNwLpOiDOS79ZZ3WfqzGe9K/L9KpPtzmRVro5KhcJxuF2OeStksiZfhKrerapbVNWjqj2q+gYrykG5lLZuV/6HV+usd2emY8gqFgX5Hp+HNfllmquJe/PH8ZnhaQR980+spdYGWGrueV9nCwYnrZ8rPzGbRCarVcl2BwAelwMuh6zrgXehBX3Nq1GqZlsu211+2zrPX2+2IBWukseaPNW1SDydb64Ean8lbU6f21KkuX40kkCGWe8qZgbuwpr8qfHoohNrj9+LmXh6Uda7UDiB9uYmeFzF50n3d7XURXN9uRaHlRARtHld67pPvpqfqVGzLdJcPxMv2x9vvH99N9eHwvFFn1HQ70E4nkYsWTyLpNXqb2QD1VQkkc73xwNAq8eFtjXMenfHIyfx5MAUJqNJtHpc+dS15sAuU4/fk2uWLZ6chRYrTPVaeCJaeFLKLwQyk0B/d0vB9otPYIW2dDXjR78aQSqThcshmImlcXo6hjPTMZyZNqbrTc4mkVXg4q2d2NvfiUgijdFIAnv7O/JJZ6q1n9UaeAcYI+zXc00+P0Bzlc31gDFw86nBqUX3h8IJvHR7d9ltzdkf9VozHplJYFfBynnA/Jz7W7tbrShWWQzy61wkkc5PITKt1QIx6UwWn7nnWXQ0N2FLZwuOjEQwPBWDz+NC74KTS34A4AyDfKXMMQxBvweduWU/01ldNPe5pyDZUGGQN3Lclw6c27pbkcxksetTP8y/9kJ+rwtZBb756Kl593/kip34+Ov3rHjfCo1WMBVrudo8TZhZx0E+FI6j1e2cd8G/UkGfB2ORJFKZbH5a7FLZ7kwigkBb8T79ejAyE8crdm6Yd19h6yeDPNWdmXh63iASYO2y3o1FklAF/vgNe3DtS7cCME4u6YzC5ZzfczQ//ar1U7YaQSicgMsh6GoxFpjZ0ObB2ZnFtfNSC4GEwglcsq2r5Ou/7cJNSKazmIqmEEtl0N3qxsYOLza2G+tqB3weNDkdyGQVz52ZwS+HptDR7MbnvvdsvlumGsxsalWtyXtc+eyP61E1ps+ZzIvIsUgCG9uNFrpKst2ZAv76TIhjZrtb+L0rl2CqHjDIr3OReApbOhaObPfgiYHJqr9Xfp5swQ+91I+e+euXL5RbXGZuEKUR5BfW5OcSlsydlMx50uUCZ5vHhfe8bNuS5XA6BBdsbscFm42Ls3/62fGqZlEMzcThX2WO9YXavK66rT3WwuhM+WO/HMGC7iAzyFeS7a5w+4Fx68d+LLRwjrypcH/rEQferXML++SBtct6t5xVruay3q3fE+9yLcxYZv4dWPB557PeFQS1as2TLsb8PlVLNWudJp/Xta4z3pVbHW65itVs5+aXL/0ewRJT8Kxmnr8Wdjl0tDTB7XTU7TiChqrJpzLLW+85mkzjocNjOHejD/1dLfOmPVSzTA6Rkjm0E+kMjoxE0OP3VnSlbA5qWouyFhOOF+mTL8h6BxipG6tRnrnR30v/0JucDmxoc+OJgSl875en0dXqxku3dy1q1ifgwKkJJNJZDE5EsSMwNyjIzFy3cH1wM+vdwHg0n1nMPKlWswncFPR58POjY1V7vVDYmG9dTW2exhx4N5swLs46WppW/BtV1aomFzIv4p89PYNNHc1IpLM4cMpoGazkPYI+LyajKZwan8VMLI3zNvkXnV8nZpO456lhtHpc2NjejJft6C47/74aSp2/RHK5ASqskKgqpqIppLJZOEXQ1epedOwyWUUinYG7Cue7hgryz58N4/3f3I8bXr0TF/Z1LPn8/zwwhL/87rMAjC/Xnl4fdgTacN1lW3FOoG2Jrec8MTCJWx88jvHZBGKpDHYE2nDOhjY8PTyFnx8dh9MhuLCvHRvbmzEaTmAqmoTCSPRybDSCVEbhEOAl27uwp8eH0UgC45Ek4qkMEuks3C4H3E4HRsJxDE/G8I69W/B3v3nhkuV66PAonjszg/e/akfF+1Iok1VEk5l5U+iAuUFvl/7N/YilMviHd12Id+zdsqL3KBQKJyACbGgrnfGq0PYNrXjo8CgeOmysWRDweXDVxVvwx6/fU5WFSezgxNgs3nnLI/n/v2p3IP/33Bryi0+smzub8cNnzuKHz5xFc5MTrWbGwzUY5BjMTdmLJTNodq++iT0UjmNvf2cVSjanzetqyIx3b7rppxiYiMLb5MDG9mb0+r3Y2O5Fb/tcpSKeyuLwSBi/HJxCNJlBd5sb3W0ebGhzw9vkxHNnZhBLZZYcFFep7lY33C4Hbrr/CG66/0j+fp/HVTbbncksx6v+9gEAwFd/68V464s2zXvOXQeG8IUfPJf//zeuvwSv2ROsQulLOztdOi1vj9+DnxwK4apbHkZGFdmsIqOKTNYI6kG/F+dsaMVYJIFHjo1jfDaZ37azpQl7en1ocjqQSGUxGklgaDKKVKY6LakVBXkRuRTAVwCcC8ANwAlgVlX9VSlFhQI+Dx47MYHHTz6OA59+3ZJXrifGZtHiduLP33wunjg1iaOjEXzz0VNIpLP4m3e8sOL3/bfHBvCTQyHs7e9AV6sHj5+YwHefOo0tnc34zX1boGpcCBwNRRD0edHV6oZDAIcIXrUngAs2teNIKIIfPn0G3zk9jKDPg+5WD9pb3HA7HUhlskims7iorxNOETx7eqaicv3744P4yfMhvO/yc1Z0FZ/PW7+guf7Sc7rwxvN7EfR7cOdjAzg8Eln2axcTmomju9VTcW38n373EgxPxuByCI6NRnDHI6dwywPH8LYXbcJ5m2r61atbJ8aMY/PXv34Bev1e7Ns6F/zeuXcLWtxOdLcuPrH+3VUX4hcnjJPNeO6iM6OKF61BXnqzibYaU4xU1ZgFUOUWB7+3Ccl0Fol0pmSegHqTSGcwMBHFq3YHsLunDWem4zg7HccvTkxgZCY+b/ZDwOfBRX0d6Ghuyh/zY6EIIok09vT48Aev3I7fePHqL+QBwOV04M4/uBSnp2JwuxzwuBzG8sedLRXVtt94QS/GIkbOhr/47rNFsyyenYmjxe3Etz9wGd7y5Z/hxOgsXlOdyRslDU3G0OZxLRqoDABXX9KH7zwxDKfDaNU1W3cdIhAxkn7tPzmBNo8Ll+8O4ILN7XC7HEimszh8NozDoTDiqSw8LgfO2+jH68/vQWeLG6l0Fh/54urKXWlN/qsA3g3g2wD2AbgOwO7VvfXy9fq9uOGKXfir7/8KE7PJJdeSHpyIoa+zBb996Vb89qXGaO4rb/75shN6DExE8cLN7fi3912Wvy+aTKO5ybms4PrxX1v6I/vzu5/G/zxztqLXC4XjiKUyGIskV9TMagb5xc31Xnztdy4GAPz4VyNV6x8LhRPLqi34vU3wbzR+ULt6fAj4vHj4locxEo7jPDDIA8Z3HADecH7Polp4f3dLyVae/u6WedPn1lI1pxiFE2nEU9mqtzjkF6mJp+Fpa4wgb45Af9MFvXj3S/rnPZbNKqZiKTjE6PpqcS/vXLVaF2/txMVbV9ba0tHixoeu2AUA+NK9h4qO5zDHn5y30Y9WtxODk2s/UG94KoYtnc1FP8erL+nH1Zf0F9lqjqqu6Bh8ZNlbzFdxg7+qHgXgVNWMqn4DwBtX+d4rsjV3YjpVQaAemowuSrLS39Wy7PzbgxMx9HfNPyG2uKvTT71Q0OfB+GyyovEH5kCPlX7BzVzd5ZKUBHzVmzO/VLKVpXCd+cWGJqPwuBxV76OupmpOMcqnX61S07KpEVeiKzeQ1eEw+no7Wtxo9azNuaoWAiVydoRmjFUIRQRbOluqOkWzlKHJGDYvmIm0HFYdg0qDfFRE3ACeEpEvicgfLWPbquqvcDUsVcXgRHRRutT+rmYMT8WQrnAQXzKdxenpWH4VrrVm1lCWCqzmYBlg5SuDLVyBrpiAr3pzVkdWmdgmUOdTVawwNFm6dlEvgvmZEqs/bvkBglW+qGnENeVDyxjI2qhKJcYZjcydS/q6mmuSbnloMrpojY1GUGmg/p3ccz8EYBZAH4B3rFWhyjGD7akl5lFORlOYTWYWHZT+rhZksooz05XVKoanYlDFopr8WsnXVpcIrJFEGrFc7vGVXsWGF6wlX0y1st+lM1mMR5bXXL+Quc58NedcN7rBycUXsvVmborR6o9bueVwV6OtEYN8PsFM/bbirFawRGKcwnn9WzpbMDwZq/qU30LTsRTC8XTd/9aKqTTI/7qqxlV1RlU/p6ofB/DWtSxYKd4mJ3r8niWb3M0ru4U1cPP/lTbZm8+rVZCfW6Ch/Amx8CJgtTV5f5kgH2irvPugnPFcTvPAKufiBn3VnXPd6IYmY4u6pOqNiJRcw365zBN+oMq1V5/H6LJqqOb6mQQcgiXHJjUyY878/O9NLJlBOJEuCPLNCCfSmI6tXcbC4VxFarONa/LvKXLf71axHMuytat1yYxIZu22b1FzfX0HebOmO7JE7dk8YTpkNX3yZnN96T55s8Y0HkmWfE4lzPL2rLLWEfR7lvxs1otwPIWpaKohahel1rBfrlA4AY/LUfbCdCXmmusbJ7VtKBzHhjaPraeTBn0eRJOZeRdfowtaMMzv/1r2yw/lzrG2a64XkWtE5HsAtovIPQW3/wUwUZsiLtZXweA5M/AtrOVsbG+GyyEVB/nBiSjcLkfNmsS62zxwCDC6ZE3eePwFvf78COvlMnN1l2uuN/s+V9vUmk8kUYWa/FKfzXpR6kK2HvX4q7PmduGAq2oyfwMNVZNfYkEhO5gbhzP3m1+YuMk8x69lv3ypJbEbwVKXww8DOANgA4C/L7g/DODgWhVqKVu7W3DXE3HEU5mS+asHJ6LoaGlaNHLc6RBs6WyuvCY/HkVfZ/OaZ1MyOR2C7rbi6zEXMmvGF2/txJ2PDSCT1WVf0UfiaYgALWVygJtBebX98tXqP+zxG8132azW7JjUq7kTT/3XLnr8XvysClnvFqburRZz8GlD9cnPJNDbbt9Bd8D8gchmArOFi93UoiY/PBVDc5MxJqjRlA3yqnoKwCkAl5V73nKJyN8CeBuAJIBjAK5X1alKtzebzocmo9gZ9BV9zuBkrGQNp6+rpeKrvoGJaM2a6k1GIFu6Ju9xOXDeJj/SWcWZ6diyrzJn4kbe+nLBMuCrbCDgUqqVNjXo8yCdVUxEk9hg477ISpjf4UYI8gGfB+F4Op9Gd6VGw4l5qXurxdvkhNvpaKwgH06sSfKiehIsMhB54dTB9uYm+L2uNZ0rb46sr+dZLKVU1CcvIpeKyOMiEhGRpIhkRKSytGzF/QjABar6IgCHAfzZcjY2k3iUq40PTSyeI5/fvsK58uY0vFoH+aDPu3RNPtdUZ17IrKTJPpJIw7fE+tFmCtrV1uRHZhLobnXn15deqZ6CdebXu6HJGFrcTnQVyWhXb6p13NayibrN2zjLzaYzWYzPrk2rRj0JFqlkjIYTcOaWVDat9Vx5c6pqI6r0jPtVANcAOAKgGcDvA7h5pW+qqvepqnnJ/CiAZeVT7F9iGl02q7mDUjw493e1YCqaWnI05lQ0hXAiXbM58iazSbqcUG7Oeb4/agVXsZF4umx/PGBMW+toaVp1n/xoOF6VlcOC+YGJ7JcfbKDaRaVTQ8uJpzKYjqXWLLA10iI1Y5EkVFc/xqXetTcvnn5pDDh0z2uBXOu58kOTsYYcWQ/UR8a73wPww+Vs0N3qRqvbWbI2PhpJIJnJoq/EQak0oU6tR9abAj4vxmcTZRP2hMJx9Pg92NTRbIywX8EXvNgys8UEq5D1bqRKq1zl++hYkzemzzXIQKBia9gv19z0ubUJ8o203Gx+fXab1+TNFd5GFzTXL/wOmDX5tZgrH44bFcJGHHQHrGHGOxH5sYg8U+R2ZcFzPgUgDeBbZV7nfSKyX0T2j46OmveV7VfP91WWCM59yw3yNcrxbQr6PFA1rtZLMQYgedHkNFafWkmQD8dTZVPazpVn9SOjzYuS1eI68wZVxdBE42Tgyk8NXcVxCy0YcFVtbZ7GWYluLr2vvWvywOLU2qPhxZkz+zqbEUtl5q3uVi3DU40zwLWY1WS8e2e5DVT1dap6QZHbdwFARH4XRkKda7XM5Zeq3qqq+1R1XyAwt4xmf1dLyeb6/PS5Us31FfTpV/I6a6Vw1a5iYskMwvH5ySAGV9AfFU4s3VwPILdW8sqDfCarGIskq3Jy9jaZ3QfruyY/E0sjnGicDFztzU1wuxyrahGqRU2+UZrr10O2O1Nwwfmn2AyLtRxhP5Qb77SavPVWqijI50bZ+wB4zIx3ueb7FRGRNwL4JIC3q+qKOlK2dhuD54pdH5iD0Epdefm9Tehoaaooa96GNnd+re1aWSrX98Kmuv5lzBYoFIkvPfDOfJ/RSGLFTWHjswlkslq1AVPVSqyylkZm4jg2GkEsmVmT1y+VB6JeiQiCPg/Orqq53sy1sFZBvgnhBhl4Z54D1sMMk2DBbKNMVjEeWdxcX2nrrCmbNdZ8r8RcIpzGuKBeqOwZXowRPZ+BUYN35O5KA/iKqv7VKt73qwA8AH6UGzT0qKp+YDkv0N/VgkQ6iz+4Yz+Cfi8+fMVObGw3TnhDk1EEfJ6Sc+jN7f/3+RD+6nu/QrPbgZNjUQxNGolvWtwuJNIZ/Or0DHYEqz9dZylL1eTnppCYCzS0IBRO4PaHT8LlFBwNRfD8mTDe87KteOMFGwEAj5+cwNcfPI6v/fbe/HrukUR60TKzxQR8HiTTWczE0mgvM0/0+wdP47ETE5iJpZBVo/aWTGfx8HFjfrR5fFarWolVVuunR0ZxYmwW1122bdFjb//qz/IXae/cuwV//64LK37dY6MR3PD/noDLKWj1uDATSyEUTsDndWFPjw+97V6cGJsF0Fgnno3tXnz3qdP4yfMhXLilA7f/3kuWldshFM6lcW1du4F3jdMnn0BXqxtulyXrhNVUoM2LyWgKyXQWUzEjPfbimrxxbrnp/iP4j/2D2NDmwUu3d6G33YvHT07guTNhbO1uwe4eH54amMJ9vzqL2UQGmzubsSV3C/i8SKQymImnEJpJ4Mx0HF2tbszEU/C4HPmZRo1mqTP8HwF4OYBLVPUEAIjIOQBuEZE/UtV/XMmbqurOlWxX6PLdAbxsRzdOjUfx4+dC2BVsw/Uv3w7AqAFvXCJJxG/u68PtD5/EnY8NIJHOoL+rBX1dLUhnFJPRJDwuB/Zt68I79m5ebVGXbUObGyKlpxvl++NyX/SL+jrgEOAz9zwLAGhuciKdzSLg8+SD/IOHRvHj50YQCiewqaMZ8VQG0WQGHS1Lf3HNq+bRSLxskP/0fz2DeCqDoM8LEWNRh2xW8ZLt3fiDV56DV+0OlNx2OQI+D46FIlV5rdX45iOncODU5KIgH0tmMDKTwJtf2IvhqTh+cWJ8Wa974NQkDo2E8fKd3UhlFFs6W7B3ayemoyk8d3YGj52cQKvbhZds68JOCy5CV+qzbz8f//t8CIdGIvjeL0/jocOjeM0LghVvH5pJoHsN07h2t7kxlQsm9R48Q1UayNoIzJabsUgCE7k+94U1+VaPC+/cuwVHQmHMJtJ47kwYdz85DABwOQTbN7Ti4WNjiKeyaPO48Npzg+ht92JoMoahyRh+9KsRjEWM877P24Sgz5Of5fT82TBetKW9IWaxFLNUkP8dAL+mqvlUVap6XER+G8B9AFYU5Ktha3cr/vUPLoWqYvenf7goWcKmJYL871y6Fb9z6VaoKjJZzddu64HL6UB3q7tMTX5+c/3luwN49nNvxGwyjWQ6ix6/F1d//ZFFg1WMbY0gb14oVNK/WbjEa6nkQ/FUBlPRFP749bvxoSt2VbinK9Pj92I0Yn3Wu5FwAuOzyUVBYSxifLav3hPE0VAE//LwSahqxScJ81j903suKdsa1WjO39SO8ze1G607R8dw52MDywryxvKiaxfYzDEj47OJqrU6rZVQlaakNoLCufKTUTPIL973wtYyVcWx0VmMzMRxUV8HWj0upDNZDExEsamjuejvqlTW0FgyA0f9hIdlW6roTYUB3qSqowDqIr+fiBhrDs8sGH1ZYb+diNRVgDeVW20tFE7A5RB0FtTCm91ObGgzptQ5Hca0k4VzS4G5HNDLmYKTn7YWKd1EvjDV5FoK+jxI5VpcrGTm0B9b8LmECgaIBc2ujmU0A4/MxOH3umwV4Au5XQ5cdfEW3P98aMnVFguFwvE1DvJzF7P1bl3V5AtW5hydqWzAoYhgZ7ANL9+5IT+myuV04JxAW8nfVakWoma3Ex5X4/4Wl4pu5c6i1p5hCwT83nkDMyZmE/mFVRqVsdpaiZp8bi3lcrXYhUs0hgpq8oX/VhKUAxWc/PJpa2uwYMbcmAXrTsbZrC76LE35UeBtnrmujmUk7wnNJGxfS7v6kj5ksopvHxiqeBszAdRaqVYK57WWzSrG1rhVo56YFbbRSCJf0VirGRZ2tFRz/YUl0tcKgLo5CwXaPPkRkOORhLFueYN/CXp8Xuw/OYnPf/9XaHI5sCPQhnM3+hD0eXF2JrZkEAj6vQjH0/lFfEYXBvlcDaqSuet+rwselwOjESNBj4gsuuqdW0q2NjV5ALj7yWF854khDE7EEArH0epx4XXn9uBlO7rhcAiS6SxSGeOWTCuSmSwmZhM4NR5FNJnB+Zv8OG+jkfs/nsrg/E3tFffFTkaTSOdG5y4c6W+eiIK5fO1A+a6Ohda6xloPzgm04aXbu/Bvjw/A39yEY6EIfF4X+rpa0J+7tTc3IRxPI5HOIODzGIFtDS8i5/Kk1+fMDVXF4ZEIDo+Ekc6q7b8jpu5WY4zSs6dnMBNL2bqVay0stUBNQ3ySQb8HTw5MAihsKq2ba5AVeek5XfjB02dw52MDSGaySGXmT/d4/Xk9Zbefq0EaffBmk7JZoxwp0uRfiogg6PfgGz8/gdt+ehznb/Lj+x9+5bznzC0lu/YnHjO95K0PHYfH5UB/VwsCPg9OT8Xygw/LEQGaHA4kF2QU/PRbzsXvv/KcisowMrO4lcQ0ao4Cb/Pkk6ssp3YYCidwybauip/fqK69dCs+cueT+Iv/egYtbifiqQyWmtW0lhfvG9o8ZQe8Wu17B8/gI3c+mf//7t7KLhobncvpQK/fi3/9xQAA4PxNfotL1FhqOwF8jQR9HozPJpHKZG3TnPOOvVvwjr1GSv9MVnFibBaHzoYxPpvATCy15ICluabHODxNjvzJ0zyBVdLkX+jDV+zCkwOTODE2i1+cmFi0zK85TqCrgouG1drY3oz/eP9laPO4sKunbd6iN8dGI/jl4BScDoHb6YDb5UCT07i5XQ60NzdhS6cxbuHwSBhHRiLwuBz487ufxpGRykfszxvvsLAmH46jq9UYBR70La92qKprtpxqvXnbizZiU7sXGzuasandi3RWcXoqhoGJKAYmopiJpeFvdqHJ6cDIdBwT0STecH7vmpWnyelAV4u77NgTKx06OwOnQ/DtD1yGvs6Whj/HLcc3rr8Ep8aj8DY5sadnfVzcVItNgrxRax+LJCoemNFInA5jEMlypksVDiIyg7FIYZ/88pqE37WvD+/a14e7Dgzh0eMTGJ6KzVvy08wnXavR7i/ZXrymuyPQVvFSpOZobwC49afHl7XIj3mx5ChS8xstyK3d5nHB2+SouHY4EzNmSKyHE7iIYF9Bi0WTU7C1uxVbu1stK9NqszuupcGJGDZ1eLG3v9PqotTcC3r9eEEva/ArUX/DylegsGnaLjX51cqPSA0n8oF9R6AtX6M0AtHyuzTMlMALM0uNzDR2P3JfZ8vygnzuc9wRaFs0QLIwyBuZ3ipP3jOSz+rW2N1NjcrIk16fffIDE9GGWZCI6octgnxhrXU0lxlsvQ/M6Gp1w+kQ48InVzM5f5MfY5EkMlk1gvIK+s/z69cvyBFtTFts3MDU19WM01Pxsiv/FRqZSaC9uQn9XS3Fa/IFszuCC6YzlrMw0RHVVjUWY1orQ5MM8rR89gjyBWtVh8LxdV+LB4wmfjOhjtm6cf4mfz7AT0ZTKxoJH/R54HY5FtXkG70fua+zBZms4sx0hcE4t6pe0D9/qqKqYnRBbu2Fz1nqdQEGeasE/caKZ5XmNa+VaDKNsUiyYdYqoPphiyBvLtIQCsdzyxDyBAnMBZdQLrmKWQv41emZ/OPL5XCIsepdQZBPprOYmK3OKnNWyS9wUWGT/UhuznbQ58X4bCLfAjAdSyGV0flB3ufNt6YsZeG6BFRbQZ8H6axiKlZfC9WYq6v1lVg+m6gUWwT5JqcDXa1uo2l6hX3NdhT0eTGa65MP+r35wPH08HTu8ZVdDC3svzZbCqqxXrxVzAsgc1nJpZhZFYN+D1SBsUgyfz8w/7MN+IypdJWsSBeaSaDV7URbjVc+JMPcWJb66pc3L6oZ5Gm5bBHkgbkMb6Fw42e7qxbzMzH7iM3A8+xpM8iv7GKor6sZA+NzQT5Uwznya2VjhxcOqawmb0xziyPo8+a7PMw8AcXWPF/ONLr1lJO8HuW7/upshH0+yLNPnpbJNkE+4PPg5NgsoslMQwebagr4PBiPGEsmBv1zKVafGV55cz1gLNM7E09jOtekOZIfLNa4wanJ6cDG9uaK1qOejBpN8mafPDDXzF5sdkdwGWl4zfwFZI1A2/zjWS8GJ2NobnI27HKnZB3bBPmgz4tjo0YyE9bkDUGfB1kFhqdiCPo88DY50d7chLMz8dy63Cs7YeRH2OcC4mi48WvygLEm9cJZA8Xks/vl+uQL7yu2ul+wYIrnUtZDStt6Vq+pbQcnotjS2dywy52SdWwT5AO5gGb+TfNT+5qfiRlAuts8K159Lz9ILRfkR2ZyaVxbG/tz7+tqqagmP1KQ939Dm3tekqHRSMJYk7qgT31uimclzfVruwgLldfidqHN46rogqyWBidj7I+nFbFNkA8umLJE8z8HM3CY962mtrhwJLo5bbHUUo2Noq+zBaFwAvFU+QFyhSv4uZwOdLfOJVAxE+EU1rg6W9xwOWTJJuBIIo1oMtPQAxjtYOEKjlZTVQxORNHXyelztHyWBHkR+WsROSgiT4nIfSKyabWvWRjQ2FxvWJiQxfjXO+//K9He3AS/14XB3Eh0u9Q+zTnIQ0s02S8caNjj9+THJRSbwulwiJEudYnAYYcBjHYQ8HkqnvJYC1PRFCKJNGvytCJW1eT/VlVfpKoXAfg+gL9c7QuaAa3SldXWg8CCaVzA4mC/Un1dLRgoaK63Qz+yeRIdWmKEfSicmLfcZWFGu8KUtoUqCvLhxh/AaAeBZWQorAWzxYxBnlbCkiCvqoVr1LcCWHV6KXME84a22i2SUu/MgXbAXOAwA9Bqm4T7u+bmyo/aZNpXqZS9C43MxNFTsL9Bn3euJh8pHuSDPs+SffJzA/oa/4KpkdVbaluzxYzT52glLMu4ISJfAHAdgGkAr1nt65knRg66my/g8yCWysDfbBxqMxgHVhmU+7pacP/zISTTWYzPJm0RmMyUvcdCERwcmkIknsYLt7TD521CKBzHUwNTSGUUR0MR9LbPfX49fmOqYjyVwcRsEoG2xZ9twOfFEwNTODE2i3DcmHroEMFkNImRmQQEwIGByVw5Gv+CqZEF/R5EkxkcG40gnVG0uJ3oanUjlsrg9FQMqYxiY7sX7c1NmJhNYiaewjkb2tDsNlp2YskMZpNptHlc8LgcEBFksorpWApT0ST8zU3obnXPG7cxHU3h5PgsfF4XOlvcyKginsqgucmJUxOzAMCUtrQiaxbkReTHAIot/vwpVf2uqn4KwKdE5M8AfAjAZ0q8zvsAvA8A+vv7S75fq8eFVreTQX6BoM+DWDKTP6FszAWn3ioE+WQ6ixd97n+gink120blcAi2dDTjXx4+iX95+CQAY3neoG+uz9106Tnd+b+3dLYgq8Blf3M/gOIXmhvbvZiYTeI1f/dA2TL4PK78BRlZw/xtvPbvH6x4G5dD8IKNPkQTGZwYn4Uu0TbpbTIGbLZ6nIilMvnaeikdLU3weZsqLg+Rac3OJqr6ugqf+i0AP0CJIK+qtwK4FQD27dtX9qfz8p0bsHfr+ltruZyrL+mbl73r4v5OfOmqF+HVewKret03X9CLwYkoVBWtHhfe/MJi13ON55Nv3IOnh6dx/qZ2tLideHJgCifHZ3H+Jj8u3toJn7cJmazinMDcmue//uLNaHIJfvL8KJ4emsKL+zsWve5vvbQfHS1NaPO44M+drLOq6GhxI+jzQAGcnY6js7WJc6Et9mvn9eCzbzsPniYnfF4XookMJqJJNDc50dvuhdvpwJnpOKZjKXS1NqHV48JzZ2bw1OAUNrU3420XbkJ3mxuRRBrxZAYQgUOMAavtzU2YjqUwPBnDZDSF2UQaLqfgmpf0Y0egDdFkGpOzKbicAq/LuAAYDSewu9dn9cdCDUp0qUvOtXhTkV2qeiT394cBvEpVr1pqu3379un+/fvXvHxERET1QEQOqOq+lW5vVbvgjSKyB0AWwCkAH7CoHERERLZlSZBX1Xda8b5ERETriSXN9SslImEAh6wuxxraAGDM6kKsITvvn533DeD+NTruX+Pao6orHpTRaMN4D62mb6Leich+7l9jsvO+Ady/Rsf9a1wisqqBaLbJXU9ERETzMcgTERHZVKMF+VutLsAa4/41LjvvG8D9a3Tcv8a1qn1rqIF3REREVLlGq8kTERFRhRoiyIvIG0XkkIgcFZE/tbo8qyUifSLyvyLyKxF5VkQ+mrv/syIyLCJP5W5vtrqsKyUiJ0Xk6dx+7M/d1yUiPxKRI7l/GzIHsYjsKThGT4nIjIh8rJGPn4j8s4iEROSZgvuKHi8xfDn3ezwoInutK3llSuzf34rI87l9uFtEOnL3bxORWMFx/JplBa9AiX0r+V0UkT/LHbtDIvIGa0pduRL79+8F+3ZSRJ7K3d9Qxw4oGw+q8/tT1bq+AXACOAbgHABuAL8EcJ7V5VrlPm0EsDf3tw/AYQDnAfgsgD+2unxV2seTADYsuO9LAP409/efAvii1eWswn46AZwFsLWRjx+AywHsBfDMUscLwJsB/BCAALgUwC+sLv8K9+/1AFy5v79YsH/bCp9X77cS+1b0u5g7z/wSgAfA9ty51Wn1Pix3/xY8/vcA/rIRj12uzKXiQVV+f41Qk38JgKOqelxVkwD+DcCVFpdpVVT1jKo+kfs7DOA5AJutLVVNXAng9tzftwP4deuKUjWvBXBMVU9ZXZDVUNWHAEwsuLvU8boSwB1qeBRAh4hsrElBV6jY/qnqfaqazv33UQBbal6wKihx7Eq5EsC/qWpCVU8AOArjHFu3yu2fGKs5vQvAnTUtVBWViQdV+f01QpDfDGCw4P9DsFFAFJFtAF4M4Be5uz6Ua4L550Ztzs5RAPeJyAExlgsGgB5VPZP7+yyAHmuKVlXvxvwTjF2OH1D6eNnxN/l7MGpHpu0i8qSIPCgir7SqUKtU7Ltot2P3SgAjmlvwLKdhj92CeFCV318jBHnbEpE2AHcB+JiqzgC4BcAOABcBOAOjGapRvUJV9wJ4E4A/FJHLCx9Uo92poad2iIgbwNsBfDt3l52O3zx2OF6liMinAKRhLHsNGMeuX1VfDODjAP5VRPxWlW+FbPtdXOAazL/IbthjVyQe5K3m99cIQX4YQF/B/7fk7mtoItIE44B+S1W/AwCqOqKqGVXNArgNdd6MVo6qDuf+DQG4G8a+jJjNSrl/Q9aVsCreBOAJVR0B7HX8ckodL9v8JkXkdwG8FcC1uRMpck3Z47m/D8Dot95tWSFXoMx30U7HzgXgHQD+3byvUY9dsXiAKv3+GiHIPw5gl4hsz9Wc3g3gHovLtCq5fqR/AvCcqv5Dwf2F/Sq/AeCZhds2AhFpFRGf+TeMAU7PwDhu78k97T0AvmtNCatmXi3CLsevQKnjdQ+A63KjfC8FMF3QrNgwROSNAD4J4O2qGi24PyAiztzf5wDYBeC4NaVcmTLfxXsAvFtEPCKyHca+PVbr8lXJ6wA8r6pD5h2NeOxKxQNU6/dn9cjCCkcfvhnGiMNjAD5ldXmqsD+vgNH0chDAU7nbmwF8E8DTufvvAbDR6rKucP/OgTGC95cAnjWPGYBuAPcDOALgxwC6rC7rKvaxFcA4gPaC+xr2+MG4WDkDIAWjj++9pY4XjFG9N+d+j08D2Gd1+Ve4f0dh9G2av8Gv5Z77ztz39ikATwB4m9XlX8G+lfwuAvhU7tgdAvAmq8u/kv3L3f8vAD6w4LkNdexyZS4VD6ry+2PGOyIiIptqhOZ6IiIiWgEGeSIiIptikCciIrIpBnkiIiKbYpAnIiKyKQZ5onVIRLoLVuo6W7BiWURE/o/V5SOi6uAUOqJ1TkQ+CyCiqn9ndVmIqLpYkyeiPBF5tYh8P/f3Z0XkdhH5qYicEpF3iMiXRORpEbk3l4oTInJxbjGQAyLyP/W+Ih3ResIgT0Tl7ABwBYyFeP4fgP9V1RcCiAF4Sy7QfwXAVap6MYB/BvAFqwpLRPO5rC4AEdW1H6pqSkSeBuAEcG/u/qcBbAOwB8AFAH5kpOCGE0YKUiKqAwzyRFROAgBUNSsiKZ0bxJOFcf4QAM+q6mVWFZCISmNzPRGtxiEAARG5DDCWzBSR8y0uExHlMMgT0YqpahLAVQC+KCK/hLGC1sssLRQR5XEKHRERkU2xJk9ERGRTDPJEREQ2xSBPRERkUwzyRERENsUgT0REZFMM8kRERDbFIE9ERGRTDPJEREQ29f8D9+97YgXOHa8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute marginal probabilities for each state\n",
    "expected_states = posteriors.expected_states[0]\n",
    "marginals = []\n",
    "for j in range(len(num_states)):\n",
    "    axes = tuple(np.concatenate([np.arange(j), np.arange(j+1, len(num_states))]) + 1)\n",
    "    marginals.append(np.sum(expected_states, axis=axes))\n",
    "\n",
    "# Plot the states and data\n",
    "fig, axs = plt.subplots(len(num_states) + 1, 1, figsize=(8, 6), sharex=True)\n",
    "for i in range(len(num_states)):\n",
    "    axs[i].imshow(marginals[i][None, :, 1], cmap=\"Greys\", aspect=\"auto\", vmin=0, vmax=1)\n",
    "    axs[i].set_ylabel(\"State {}\".format(i))\n",
    "    axs[i].set_yticks([])\n",
    "    \n",
    "axs[-1].plot(data)\n",
    "axs[-1].set_ylabel(\"Data\")\n",
    "axs[-1].set_xlabel(\"Time\")\n",
    "axs[-1].set_xlim(0, 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a02f5d9d-865e-49c6-9f47-49b66cdeeb6c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[-1.5112934  -0.63249874]\n",
      "  [-2.3792994  -1.5005047 ]]\n",
      "\n",
      " [[-2.1733613  -1.2945665 ]\n",
      "  [-3.041367   -2.1625724 ]]]\n",
      "[[[-0.29943854 -0.7436352 ]\n",
      "  [-3.2376444  -3.681841  ]]\n",
      "\n",
      " [[ 1.3292058   0.88500917]\n",
      "  [-1.609      -2.0531967 ]]]\n"
     ]
    }
   ],
   "source": [
    "print(factorial_hmm._emissions._distribution.loc)\n",
    "print(test_factorial_hmm._emissions._distribution.loc)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}