hmm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hidden Markov Models\n",
    "\n",
    "For users already familiar with the interface, the [API docs](../api/index_markov_hmm.rst).\n",
    "\n",
    "Hidden Markov models (HMM) are a type of Markov model where the underlying Markov process $X_t$ is _hidden_ and there is an _observable_ process $Y_t$ which depends on $X_t$. A nice introduction into HMM theory and related algorithms can be found in <cite data-footcite=\"rabiner1989tutorial\">(Rabiner, 1989)</cite>. When fitting such a model, the goal is to parameterize a Markov model over $X_t$ by observing $Y_t$. Here the dependency of $Y_t$ on $X_t$ must only depend on the current time step, i.e.,\n",
    "\n",
    "$$ \\mathbb{P}(Y_t \\in A \\mid X_1 = x_1, \\ldots, X_n = x_n) = \\mathbb{P}(Y_t \\in A \\mid X_n = x_n), $$\n",
    "\n",
    "where $A$ is some measurable set in observation space.\n",
    "\n",
    "The states of process $X_i$ are often called \"hidden states\" or \"macro states\", whereas $\\mathbb{P}(Y_t\\in A \\mid X_n = x_n)$ is referred to as output probability or emission probability. If the emission probability is a discrete probability distribution, the corresponding states are often called \"observable\" or \"micro\" states.\n",
    "\n",
    "In particular, this means that here - as opposed to in ordinary [Markov state models](mlmsm.ipynb) - one hidden state maps to a probability distribution over the observable space in which $Y_t$ takes its values.\n",
    "\n",
    "HMMs are generally more expressive as MSMs, as each MSM can be written in terms of an HMM with discrete output probabilities: The transition matrix can be used as hidden transition matrix and the emission probabilities reduce to the identity, i.e., there are as many hidden states as there are observable states and each hidden state maps to a delta distribution. HMMs relax this by allowing for other distributions too.\n",
    "\n",
    "Due to their more general nature, HMMs can recover from situations and produce good results when MSMs would not, for example in some cases of bad state space discretization. However, HMMs are harder to fit, computationally more expensive, and the estimation likes to get stuck in local optima.\n",
    "\n",
    "In deeptime, a hidden Markov model is determined by\n",
    "\n",
    "- a hidden transition model $P\\in\\mathbb{R}^{n\\times n}$ (of type [MarkovStateModel](../api/generated/deeptime.markov.msm.MarkovStateModel.rst#deeptime.markov.msm.MarkovStateModel)) describing the hidden process,\n",
    "- an [output model](../api/generated/deeptime.markov.hmm.OutputModel.rst#deeptime.markov.hmm.OutputModel), which holds information on the emission probabilities,\n",
    "- and an initial distribution $\\pi\\in\\mathbb{R}^n$ over the hidden states, which expresses a prior belief over the average most likely distribution of hidden states.\n",
    "\n",
    "Creating an [HMM](../api/generated/deeptime.markov.hmm.HiddenMarkovModel.rst#deeptime.markov.hmm.HiddenMarkovModel) is as simple as importing deeptime and fixing at least the transition model and the output model.\n",
    "\n",
    "Here $P\\in\\mathbb{R}^{2\\times 2}$ is the hidden transition matrix and $E\\in\\mathbb{R}^{2\\times 3}$ are the emission probabilities mapping from two hidden states to two three-state output probability distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:11.016869Z",
     "iopub.status.busy": "2020-12-02T13:38:11.014808Z",
     "iopub.status.idle": "2020-12-02T13:38:11.953939Z",
     "shell.execute_reply": "2020-12-02T13:38:11.952858Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import deeptime.markov.hmm as hmm\n",
    "\n",
    "P = np.array([[0.95, 0.05], [0.3, 0.7]])\n",
    "E = np.array([\n",
    "    [0.1, 0.1, 0.8], [0.5, 0.5, 0.]\n",
    "])\n",
    "\n",
    "ground_truth = hmm.HiddenMarkovModel(P, E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This hidden Markov model has two hidden (macro) states and three output (micro) states. The initial distribution defaults to a uniform distribution, which prefers no hidden state over the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:11.959144Z",
     "iopub.status.busy": "2020-12-02T13:38:11.958677Z",
     "iopub.status.idle": "2020-12-02T13:38:11.962449Z",
     "shell.execute_reply": "2020-12-02T13:38:11.961764Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# hidden states: 2\n",
      "# observation states: 3\n"
     ]
    }
   ],
   "source": [
    "print(\"# hidden states:\", ground_truth.n_hidden_states)\n",
    "print(\"# observation states:\", ground_truth.n_observation_states)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now generate an hidden and corresponding observation sequence from the hmm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:11.973447Z",
     "iopub.status.busy": "2020-12-02T13:38:11.972815Z",
     "iopub.status.idle": "2020-12-02T13:38:11.977425Z",
     "shell.execute_reply": "2020-12-02T13:38:11.976837Z"
    }
   },
   "outputs": [],
   "source": [
    "(hidden_trajectory, observation_trajectory) = ground_truth.simulate(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizing the hidden trajectory, one can observe that state 0 is preferred. The observation trajectory jumps between three states. If the hidden state is $0$, the observation trajectory is most of the time in micro state $2$, otherwise it is in state $0$ or $1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:11.984760Z",
     "iopub.status.busy": "2020-12-02T13:38:11.984280Z",
     "iopub.status.idle": "2020-12-02T13:38:12.840847Z",
     "shell.execute_reply": "2020-12-02T13:38:12.840416Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHxCAYAAACS8O5DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACHVElEQVR4nO2deZwcV3Xvv6eXmR4tM5Isr7JkyWAWY8AY2ewEwmYTwCxJMCFsITgkEPBLWENeIHl5IXm88LJAcBzisIRgCKtJDIYQgoEY8L7hTd6wvEiWLWtG0nTPdPd5f1RVd3V1rd3Vm+Z8P5/+zHR1LffWrbrn/s6991xRVQzDMAzDmCwKo06AYRiGYRjZMQNuGIZhGBOIGXDDMAzDmEDMgBuGYRjGBGIG3DAMwzAmEDPghmEYhjGBDMyAi8j5IrJbRK6P+F1E5G9EZIeIXCsipwwqLYZhGIZxqFEa4Lk/BXwM+EzE72cAJ7ifpwCfcP/GsnHjRt26dWs+KTQMwzCMMeeKK67Yo6qHB7cPzICr6iUisjVmlzOBz6gTSebHIrJORI5W1fvizrt161Yuv/zyPJOaC82mUm8qIlAu9ufYWKo3ASgXBRHJI3nGBLHcaBIXX6lYEIoFey7GDe+97ZdDoXzH7RlWVZYbnQnKu3715znvc4vIXWHbB6nAk9gE3O37vtPdFmvAxxFV5Xkf/T537DkAwJ++/CR+/anH9XSuD190I39/ye0APOfRh/OpN52WWzqN8eerV+3kf3zhmth9ZislLnnPc1m3amogaXj131/Ky04+htc+pbdn+FDiPV+6hvWrpnj/ix8bu9/Hv7eDj1x8cy7XXD1V5Lu//xyOmqvkcr4k3nnBVWzZsIrff+GjcznfN665l3dccFWsAV87XeJ7734OG9dM53LNJH7jU5fxvZsf6Nj2uqcex/96+Um5nP8/fraLsz97OU03zz9+//OGUn6jNOBhzZPQIheRs4GzAbZs2TLINPVEvancsecAzzphIz+94yF27N7f87lu2bXAUbMVNq6d4tZdvZ/HmEx27N5PQYisTG/ZtcDXr76X3Qu1gRnwq+9+mBOOXDOQc08a1+7cl8rI3LJrgfWryvzms47v63p37DnAl67Yyb37FodmwK/buY+Faj238+3YvR9VeNcLHxWqQnfs3s9Xr7qH+/dVh2bAb9m1n5M2zXLGSUcDcMFlP+fW3Qu5nf+2B/bTVDjn+SdQLhZYUxmOaR2lAd8JbPZ9Pxa4N2xHVT0POA9g+/btYxe8vbrcAODZJxzOTfcvUKs3+jhXk80bZnjUkWv51vX355VEY0KoLjeZKRd523MfGfr7d2/cxdevvrf1zOWNqlKrN6ku5+MOnnSce5F8r6vLDY5YW4kst7RcetuDfOmKnQMr3zCqy41cr1etN5gqFXj7L54Q+vv3b3mAr151T1/1ZFZq9QZPOPbwVvn8+PYHc220eO/L25/7SEp9dqFmYZTTyC4EXu+ORn8qsC+p/3tcqbl9X9PlAtOlArU+Kr9avcF0qch0qdg6r7FyqNUbTJeLkb9Pl4rufoN5Nrzz2rPnUFtupLoXtXqT6XL/1al3jmHe/1q9mev1astNpkvR98L7rZ96MivBNOVdv9bqDUoFGarxhgEqcBH5PPAcYKOI7AQ+CJQBVPVc4CLgxcAO4CDwpkGlZdB4rddKqUilXKTapwLfsLpApVwYaivcGA+qy00qMZVfxa3gB/VseJWqPXsO1QwKvFKKbnilxTtHbYIVeK3eoBLTCPV+66eezEo1kKZKuZDrPa4uN2PzPCgGOQr9NQm/K/C2QV1/mAxKgdebSr3RHHqrzhgdjpJLocAHpF48t6YpcIcsCnzNdP/VqSnw/Gk0nRHog1bgcXkeFKmuKCJHisg/isg33e8nisibB5u0ycFrvU7npMCny4WW0rKKdGVRXY6vCFoKfEDqpWoKvIP0CrzZalz1Q0udDun+1xtN6k3NvQ98nBS41ygNKvBc8zwiBZ62yfAp4GLgGPf7LcA5A0jPRJKvAm+6CtwM+ErEFPj4UG80aTQ1pQJv5NMHPuT3fhBjHsZNgXvXWbEKHNioql8EmgCqWgesie6SZx94bblBpVwYekvcGA+cvtTRKfBWhW7PHdV6em9EbbmZTx/4kN/7WoY8ZjlnGgU+7EbKoBV4XMN7UKQ14AdE5DDcedreqPGBpWrCGIgCNxf6imTUCtyr1Oy5azdiavUmGheVhAEo8CGp00GUd1I3kPfbsBop7S7OTgXujTHKg1Ep8LSjLn4PZ9rXI0TkR8DhwK8MLFUTht9FM10q9PwyNJvKUqPpnsdrpZoSWknUlhtMr40ObjHohp0p8Db+e+y8l9ENqyS3cVpKBaEgw1enS24jJY/wn7V6k7UxgUxG1U3gLz/vPVrKaZCwI7zG14DfAPwC8GicCGo3Y0uRtvAPkqiUiz23LP2unvZ0IVNCK4kk9+Og1Yt33qop8I57nDRILWngVlpEpK86JCv+6yQ9e1nOGXeeUrFAqSBDz2PF5yGplNr1ax4BDWvLjYFFRowjrRG+VFXrqnqDql6vqsvApYNM2CSRlwL3GgIdCtyU0IqiluB+FBGm+njGEq9vCryF/x7HecLCpin1Qz91SFY68piTWEijRkeRx04Fnq+HcywVuIgchbPAyIyIPIl2/PJZYNWA0zYxVHNS4J7a7lDgpoRWFNV6s0MphFEpDS7Ij1eh2XMXUKcxxi1smlI/DFOB+xtq1XqDOSfWVn/nTOGNqJSLQ+sebJePT4Hn7OFM8joMiiQX+ouAN+LEKf+ob/sC8AcDStPEEabAe+lPMgVuOAo8viKYLg8uzK5XoTUsiFBqBR42TakfhqlOqwNQ4NUU4wGmS4WhdQ9Wl0MUeM5jjMZSgavqp4FPi8irVPXLQ0rTxNHuYyn6XDPZ+5NMgRupFHjOYSD9dCqyJmtWsAEP9oFH7ncIKfBczplagQ/Lhb5yFTgAqvplEfkl4HFAxbf9TwaVsEmiVm9SEGcEqX+EZdYCNQW+svEChyQq8AEudNOpyBq5hAedVEyBZ0dVUynwqQF2AwWJVeA5pWFUCjxtKNVzgVcDv4vTD/4rwHEDTNdE4cUvF5G+BkcE55P7txmHPu3BNmkGAA12MRN/elYqaQd4hQ2S6gengTZ8BZ7HNZcaXh02um6gIC1h5FPgedevea1Gl5W0V3y6qr4e2Kuqfww8jc61vFc0Thxc51Z60xN6ac2GueItEtvKwV/+cTgu1kEpcL/beGU/e9WU7uWwaUr9MF0eYv+wz4Dlcc1qSm/EIAdiBml1TfoaWHlGvFt2PWd5ROLLStonbtH9e1BEjgGWgW2DSdLk4Slw6G96QnAwnHOela2CVhKmwMeL0SrwYcUJz1eBt9XuylHgfs/psEl7xX8TkXXAR4ArgTuBCwaUpokjTIH30pr1D4Zph1Rc2SpoJTEOCrxmCrxF2gFeeSvwQQ5SDFLLWYHXWmo3WYEPLY8hXoE8FXja93YQpB2h8n9UtQZ8WUT+DWcgW3VwyZosBqHARWSog1mM0TMOCrxqCryFKfAezjeGCrxabzDl1qmt6w9CgY/rIDZ8UddUtaaq+7BIbC0GocC9vytdBa0kTIGPF53TyIarwIe9Ghnk2wc+bgo8mJ48xxiNrQK3SGzpGIQC9/6udBW0khgXBS4CqqbAa/V092KSFXh1ueHL4zAVeGFoMS6cleI605OrAs95GmEWskRi+0vaBtwisfmo1Zut+bJeK7yXUejBdWuHGezAGD3twTCjC4JRqzdYO11ivlpf8c+eNw9+IeFehAUK6YdKeXANtCC1ejPX8k7fB14crgIPlE2+LvR0jZZBYJHYcqC63OSw1a4Cd1vhvUQ18lwxU8W2Al/pbsyVRNi6xWEM8rmoLTeZW1Vmvlpf8c9edbnJ6qkSi0uNBBd6/gp8uaE0mkqx0P/ynnFUlxvMzuRX3tUxVODVend44tYYo1xc6OPfB36siMyKwydF5EoReeFAUzZBOC6attGF3hX4VLFAwX1pp8vmQl9JtD0w6VZyUtUBpKHB3Ey5Iz0rFSckaCHR4xE2Takf2mu+D74BVas3WTVVZKqYT12T1p08XSq24u0Pmqi12vPqomxH0BzfeeC/oarzwAuBI4A3AX8+sFRNGM4gibbbG3obHFFdbnRUApWSDWJbSbQVePIIXtV21Ks8qdWbzFZcA77Cnz1vDfAkj0feCqyfgbBZcUKAJucxy/kgzUDM4cW5iAprndcg4fYaFuOrwD0/zouBf1LVa3zbog8SOV1EbhaRHSLyvpDfnyMi+0TkavfzR+mTPj6EKvAeHkzvZfIwBb6ySBsQYpBBfqrLpsA9vPc6SanVQqYp9UPea1XHUXXXn8+rrknfDTS8SJNeHrvSkFOeR6nA084Dv0JEvo0Tfe39IrIWiM25iBSBjwMvAHYCl4nIhar6s8CuP1DVl2RM91hR9Snw6T5az7XlRueKOaUiDx9czieRxthTSzkdxe/l8dRybmlwB2SK2DQy771OUmph05T6Ie+VsuKo1ZvMlItM5+TtG1cFPjvT/Z7k5eGsjVCBpzXgbwZOBm5X1YMichiOGz2O04Adqno7gIhcAJwJBA34xONX4KVigVJBel7MxN9SNAW+ssgyjQzyWT0qiLcsok1hdN7r1dMlplIo8DxHIOe9VnUc1eUG62bKpsD7YOz7wFW1iTOV7A9F5P/iLG5ybcJhm4C7fd93utuCPE1ErhGRb4rI49KkZ5xoNJXlhnYFyu8pkEtgTVnrA19ZeHNypxLW4K4M0MVac9cjtyBC7T7wpHtRDZmm1A/DVuCVcjG3aV1pG6Fj0QeeU/069n3gIvLnwDtx1PPPgHeIyIeTDgvZFhw2eyVwnKo+Efhb4GsR1z9bRC4XkcsfeOCBNEkeGlGB8k2BG1nxyj+pL7Wfbpp0aXAV+JBWxBpXsvSB56m+8l6rOg4n7fn1gdfqDUoFoZTQCG17GYZhwFe4AscZvPYCVT1fVc8HTgd+KeGYnXQuOXoscK9/B1WdV9X97v8XAWUR2Rg8kaqep6rbVXX74YcfnjLJwyFs2kSv7sdgRTA9xGAHxuipLaczBO1BTvlWfq1lEVtTp1b2s1ert/vAYw14xDSlXhnmSoTV5SbTLQWeTyjVNPei3Qgdhgs9ahpZPs+4F7GvXBzsnP0wsjx163z/z6XY/zLgBBHZJiJTwFnAhf4dROQoceWGiJzmpufBDGkaOcH45d7/vU0j63TFDTPYgTF60rpi22vO51v5+UOCOtOKVvaz5xi3QmLAj2q9kThoKwt5rpSVRM03Cr2X4FNd50t5LwbVCA0jKk2VnNZdry43qJSKuc1CyELaQWx/BlwlIt/DcY0/G3h/3AGqWheRtwMXA0XgfFW9QUTe6v5+LvDLwG+LSB1nzfGzdBDRKQZImAJPGvQSea4QBb7kBuwYxcNhDJe0rthBVX6tAUjlQm7qZJJpu5cPYQVe9zdSDj0FrqpDUeCjWAscUhhwESngTBl7KnAqjgF/r6ren3Ss6xa/KLDtXN//HwM+ljHNY8UgFbh/oMcoVroxhktqBV4eTOXXmgJUKuamTiYZJ4Z2kUpSIJd6I3SaUq8MS4GrKku+boJ8FHi6uqoyJAXuBTsKmyWQtwIfBYkGXFWbIvJ2Vf0iARe4Mfg+cO8aZsAPfVIr8AENAAoq8INL9VzPP0k0m8pSo5lqgNekKnB/4KD8FLgT1CaJYSnwuCh5eY0xGqUCT3vV74jIu0Rks4hs8D4DTdmEELYWbKXc24MRpcDzaBkb48/IFbhvUY6VrsD9AUmSphvl3Qee51rVcbRXDjt0FXgtxEPqkdcYo7FW4C6/4f59m2+bAsfnm5zJI2zeY18KvByuwI1Dn1q9waqp5FdyUArcPyVypfeBt6cGHcoK3F/e+Shwb1BcEq0FW4bUSIlS4HmMMRrrPnAAVd026IRMKmGhA3tZr1lV3WkrYX3gK7ciXUnU6k02rE6vwPN+Lqo+RbbSYxAEFbi3clbY/Oa8x6gMy4BXAwq8Vm/0bcyq9SbrUowHqAyoERokToHnNcbIv5jVsEkbyOVtIrLO9329iPzOwFI1QYSFDuxlZZ+lRhNVQhX4SnZlriSqKeeBe5Ha8n4uggp8JT93/vd6utWVFX4/okJ19kqea1XHEVTgTYXlRn+TgNIq8HJREBm8Ak/qA4f+PZzVemPs+8DfoqoPe19UdS/wloGkaMLwz5316EW9RLnind9Mga8EgpH4ougn3n4c/squ12iChwqdc+KjI6N5nrM8Dbhz3cF7QKrLIXnss8wdd3JyI9RrpAw6zkW7kRLSB55T/Zp3F0oW0l61ID6/irvS2NRgkjRZtAex9beOd3sEcKcr3vlt5SqhlYQXFSsNvcbbj8PvbnQGYq7c587/XldiFHh7JHe+LtRhxKJvl7cvj32WeW25kXpltl4H+2ZLj9dN0J2mvOrXas6L2WQh7SC2i4Evisi5OIPX3gp8a2CpmiByU+AR09Gca6xcJbSSiIrZHMYgFHItoMCXGk2aTaVQWHlBhNIq8LSLd2RlGGMQBqHAqxkGdA0j2l/1EFfgaQ34e4HfAn4bJ5DLt4FPDipRk4R/7qyHN+hludGknBDU3yNssIUp8JVFlvn+w1DgzrYmM1MrLwZBqAIPud9p13DPyjBWIvQr8OlcFXj6Z3jQ4iRure686te04WMHQdpR6E3gE+7H8BHad+0b3ZjWgIcNtjAFvnLwBw5JwyD6SP3Psv/ZW4kGvJZSnU6yAg/1MuTSBz66ZzgsPc61VrACF5ETgA8DJwIVb7uq2jxwd9Slf+qFPxTimul0Tg5T4CubsOmIcUwPoI/UH5RopT971ZTqNCyQUx4MQ4Gn9TKkpd5oUm9qJgU+zDyGXd/Zp/8+8FEp8LTNhn/CUd914LnAZ4DPDipRk0TYCNRe5nFaH/jKxh84JA2DUuAFgVJBVvyzZwq8j/OtIAXudJXq2I9Cn1HV7wKiqnep6oeAXxxcsiaHYPQ08M8vTP9gJLnijUObniq/ASjwStlZFnFYoS7HlbZHpK1Ow+5F3DSlfhhGJLwwBd5Pebfnzqf0Ig3RyxA6DzyHPC/FNBCGQdpBbFV3VbJb3SVC7wGOGFyyJoew+NW9uKPCXHHtQC4rUwWtJFrln8H9+PDicq5p8HuThrXYxLjiN0Zx72E1ZppSPwwjFn3akfZZz5cmnr+3396Dw8pjiAs9h/o1zkU/DNJe9RxgFfAO4MnArwOvH1CaJoqwFaR6cUeFPWjFglAuyopVQSuJcVDgzmCcYkc6VuqzF1ypy7+tc7/JVeAd0wZzKO84d3UYTh4HH8ilVJDQELijyHPepFXgitPnfRzgBbr9B+AJg0jUJBGmwHuZkhE1GGYYg1mM0dOLAs99OdF6o/UsD2tN6nHF73ptdSfEKfCcFdgwFHi13mCqWKBQkFzKO6sanS5nDzmdlWrM1MxDQYGnNeCfA94NXAeszCZ5BINU4DCcwSzG6Bk7Be6pzhU6Cr1WbzLlzi6JU2rtwYcDUOBDiFIW7DI5FBV41ACzlaTAH1DVCweakgllkH3gMJyBHsboyTodyVm/eQgKfIWOQq/6QoLGKbVBKfC81qqOwx8CNE8FnrYRWhmxAs9jjNGkKPAPisgnge8CNW+jqn5lIKmaIGr1BnOB5fNMgRtZiVu3OAxT4IPFvyhHa+WsMAWeceR1WvJaqzqOUAXeR3mPpwKPDrKSxxijSVHgbwIeg9P/7eVWATPgMQo80zzwiMEWldLKXlRipZA1kMug+sC9wEM2jaztjRAR5z2MWcxkEH3g3vkHFSTE73EpFwsUC/0Zs6xqtFIusFQfbLz96nL8QiP91q/+cLSjIK0Bf6KqPn6gKZlQqnF94BkUUjUiHJ+jwFemG3MlETdfNYzpUoF6U6k3mqEjbHuhttzksNWdCnyldt/4vREQPeDKc6FP5VQGrev51qoelAHvymOpP5d2LwocYKnRpFIYUB4Tlnp1uir6caGPVoGnfep+LCInDjQlE0pYHNxeBkeEBYQBz1W6MlXQSiJz5TeAaV7OM+i5VE2BB6Mihr2HcdOU+mEYkfD85e1dsy93cg+NUOe4wT1jXqjr6DTko8DT9vvnTdqrPhO4WkRuFpFrReQ6Ebk26SAROd09ZoeIvC/kdxGRv3F/v1ZETsmagVHjRa/y08v0hOpyM3LN2pU6kGglkd39mP80L+cZNAUO3YOfot7DuEFS/TCMWPQ1X3l71+xrQFcP3UDOcYN7xqoJXRB5KfC00z/zJq0L/fSsJxaRIvBx4AXATuAyEblQVX/m2+0M4AT38xSceOtPyXqtURLmookb9BJ7HlPgK5bs7sdBKPD2SlKFgjBVXLkDKLMo8EHEwR6WAl+/eqrjmrko8AxTIZ3jBqzA107HpGGyFXja5UTv6uHcpwE7VPV2ABG5ADgT8BvwM4HPqKriuOnXicjRqnpfD9fLzIP7a9x8/0LPxyvhhtcb9HL7ngP89449qc5138OL4eH+ykUeXlxKfR5jMtmxez+Q3v3oqYpLb3uQo+cqCXun4+BSvUNJTJcL3JnhGT6UePDAElsPW936XikXuX++2nUv7t67OFAF/tM7H+KBhVrC3r3x0MEljp6b6bjmffsWey7vW3c5z3CWYEQAP77jQXbuPdjTNZPYt7gcWz6VcoHdC93lmhbPfoy7Au+FTcDdvu876VbXYftsAjoMuIicDZwNsGXLltwSeOXPH+Ytn7m87/NsWFXu3rZ6in+/9j7+/dr0bZGnP+Kwrm3rV02xa77Gr33yJ32l0Rh/1q0qpx6Nu36Vo5x+/1+vyTUNG1a3n+XDVk/xrRvu51s33J/rNSaFJ21e3/p/w+op/uvmB0Lfw5M2zeZ+ba8cPvDV63M/t59nPvJw3zWn+O/bHuyrrlkzXaJcTPkMu3l8z5cSe2P7YoPPy9D126opvnvT7r7yPF0qUJkajQIXR/wO4MQivwK8SFV/0/3+OuA0Vf1d3z7/DnxYVX/ofv8u8B5VvSLqvNu3b9fLL+/f6ALsPbDELbt6V+AApaLw+E3rmAoop3seXmTnQ9lalY88Yg2Hrel09xxcqnP9PfMMqpyM8eGYdTNs3rAq1b7NpnLdPfty7aMuFITHb5prKZb79i3y8wcHo4wmgROPmWVtxTEycXXFto2rOWI2Hy+Ih6pyw73zHKjVcz1vkMdtmmtNHXxwf63lCeqVo+dm2HJYumdY1XmGF5cGO87i8cfOsWoqXKs+fHCpLy8swJGzFbZuXJ28Yx+IyBWquj24fZAKfCew2ff9WODeHvYZGOtXT/GU47tVbx5sWjfDpnUzyTsmsGqqxGnbNuSQIuNQolAQnrh53UCvcfTcTIeLdSUzyLoiDBHhpE1zQ7sewGFrprsExCAREZ5w7LqhXS+MdauGW655M0jdfxlwgohsE5Ep4CwgGI71QuD17mj0pwL7htX/bRiGYRiTzMAUuKrW3bXDLwaKwPmqeoOIvNX9/VzgIuDFwA7gIE7Et1iuuOKKPSLSy6C6KDYCh8ooHcvLeGJ5GU8sL+PJoZKXPPNxXNjGgfWBTwoicnlY38IkYnkZTywv44nlZTw5VPIyjHyMZuicYRiGYRh9YQbcMAzDMCYQM+Bw3qgTkCOWl/HE8jKeWF7Gk0MlLwPPx4rvAzcMwzCMScQUuGEYhmFMICvagCetljbOiMhmEfmeiNwoIjeIyDvd7R8SkXtE5Gr38+JRpzUNInKnu8rd1SJyubttg4h8R0Rudf+uTzrPKBGRR/vu+9UiMi8i50xKmYjI+SKyW0Su922LLAMReb/77twsIi8aTarDicjLR0TkJnflw6+KyDp3+1YRWfSVz7kjS3gIEXmJfKYmsFy+4MvHnSJytbt93Mslqg4e3jujqivygzM3/TbgeGAKuAY4cdTpypD+o4FT3P/XArcAJwIfAt416vT1kJ87gY2Bbf8HeJ/7//uAvxh1OjPkpwjcjzN/cyLKBHg2cApwfVIZuM/aNcA0sM19l4qjzkNCXl4IlNz//8KXl63+/cbtE5GX0GdqEssl8PtfAn80IeUSVQcP7Z1ZyQq8tVqaqi4B3mppE4Gq3qeqV7r/LwA34iwEcyhxJvBp9/9PAy8fXVIy8zzgNu1tJb+RoKqXAA8FNkeVwZnABapaU9U7cIIxnTaMdKYhLC+q+m1V9YKL/xgndPPYE1EuUUxcuXiIiAC/Cnx+qInqkZg6eGjvzEo24FEroU0cIrIVeBLgLanzdtdNeP64u519KPBtEblCnNXnAI5UN7Su+/eIkaUuO2fRWRFNYplAdBlM+vvzG8A3fd+3ichVIvJ9EXnWqBKVkbBnapLL5VnALlW91bdtIsolUAcP7Z1ZyQY8bM27iRuSLyJrgC8D56jqPPAJ4BHAyTjLsv7l6FKXiWeo6inAGcDbROTZo05Qr4gT+/9lwL+6mya1TOKY2PdHRD4A1IHPuZvuA7ao6pOA3wP+RUTyXyM0X6KeqYktF+A1dDZ6J6JcQurgyF1DtvVVNivZgI90JbQ8EJEyzoPzOVX9CoCq7lLVhqo2gX9gjNxncajqve7f3cBXcdK9S0SOBnD/7h5dCjNxBnClqu6CyS0Tl6gymMj3R0TeALwEeK26HZOuS/NB9/8rcPomHzW6VCYT80xNarmUgFcCX/C2TUK5hNXBDPGdWckGPM1qaWOL21/0j8CNqvpR3/ajfbu9Arg+eOy4ISKrRWSt9z/OYKPrccrjDe5ubwC+PpoUZqZDSUximfiIKoMLgbNEZFpEtgEnAD8dQfpSIyKnA+8FXqaqB33bDxeRovv/8Th5uX00qUxHzDM1ceXi8nzgJlXd6W0Y93KJqoMZ5jsz6pF8o/zgrIR2C07L7gOjTk/GtD8Tx/1yLXC1+3kx8FngOnf7hcDRo05rirwcjzM68xrgBq8sgMOA7wK3un83jDqtKfKyCngQmPNtm4gywWl03Acs46iFN8eVAfAB9925GThj1OlPkZcdOH2Q3vtyrrvvq9zn7hrgSuClo05/irxEPlOTVi7u9k8Bbw3sO+7lElUHD+2dsUhshmEYhjGBrGQXumEYhmFMLGbADcMwDGMCMQNuGIZhGBOIGXDDMAzDmEDMgBuGYRjGBGIG3DAMwzAmEDPghmEYhjGBmAE3DMMwjAnEDLhhGIZhTCBmwA3DMAxjAjEDbhiGYRgTiBlwwzAMw5hAzIAbhmEYxgRiBtwwDMMwJhAz4IZhGIYxgZgBNwzDMIwJxAy4YRiGYUwgZsANwzAMYwIxA24YhmEYE4gZcMMwDMOYQMyAG4ZhGMYEYgbcMAzDMCYQM+CGYRiGMYGYATcMwzCMCcQMuGEYhmFMIGbADcMwDGMCMQNuGIZhGBOIGXDDMAzDmEDMgBuGYRjGBGIG3DAMwzAmEDPghmEYhjGBmAE3DMMwjAnEDLhhGIZhTCBmwA3DMAxjAjEDbhiGYRgTiBlwwzAMw5hAzIAbhmEYxgRiBtwwDMMwJhAz4IZhGIYxgZgBNwzDMIwJxAy4YRiGYUwgZsANwzAMYwIxA24YhmEYE0hp1AkQkc3AZ4CjgCZwnqr+ddT+Gzdu1K1btw4pdYZhGIYxWq644oo9qnp4cPvIDThQB35fVa8UkbXAFSLyHVX9WdjOW7du5fLLLx9uCg3DMAxjRIjIXWHbR27AVfU+4D73/wURuRHYBIQa8DzZs7/GdTv3ceq2DayZdm7FrbsWOOHIta19rr77YR4+uBR5jmJBOHXrBirlYsf2O/Yc4K4HDwBwynHrma2UAbj34UVu2bUQeq7HHDXLUXOVjm3z1WWuvGsvAMcdtpptG1cDUF1u8MBCjc0bVgFQqze4/M69LDeaXeddt2qKkzevA0BVuezOvRxcqrd+nykX2b51A8WCAHDdzn08eKAWmecgIsIpW9ax1s3jzr0H2bF7f6pjTzhyLZvWzbTyetXPH0ZVE487Ym2FE4+ZBaDeaPLTOx9iqd6d99XTJbYftx4RJ29X/nwv84vLkectFQqcum090yWnPO9+6CAb10wzM+V8v3PPAe50yzXI4zfNcdiaaQAeWKhxw737EvPh59j1q3jkEWtS7dtoKpfd+RDV5QZrKyVO2eLkUVW58ud7WajWqZSLnOorV48b7t3HAws1SoUC27eubz27wbze9eAB7tjTmdcnHruO9aunAOf9EWjl2WPf4jK15QZHzHY+ywdqdS6/ay+qyrHrZ3jkEc57Vqs3uOyOvdSb3eWXlmBeb39gP1s2rKJUdHoJb7p/nvv3VSOPFxGetGVd6z0Nsnu+ys/umwfgsUfPcmQgb4NCVbnirr3sr9UT9y0XnfL0nt1bdy1wz8OLscc8afN65lZ15nn3QpWf3evkNaxOWqguc6DW6NoerJOSaLrP8OJyo+u3YJ102wP72XrY6q5n+ecPHuSI2enWM3zbA/u5+6GDnXncsp65mfBy9Xhwf43r7ul8Xx9x+JrEvNTqDe7fV+W4w1bH7jcoRm7A/YjIVuBJwE8C288GzgbYsmVLbte78q69nP3ZK/jG25/J44+d4/p79vGSv/0hX/mdp3PKlvXc9eABXv7xHyWe5w9/6bH85rOO79h21nmXsmveMYJveNpx/PGZJwHwO5+7kqvvfjj0PKdt28AXf+tpHdv+8uKb+fSlTuPr6LkKl77/eQB85tI7+dv/3MG1H3whIsKXrtjJB756fWQav//u53DcYau5/K69/OrfX9r1+z+/+Sk884SN7F6o8tKP/TAxz0F+85nb+MOXnOj8/+nLuen+8EZKkCccO8eFb38m0JnXJIoF4ZoPvpA10yW+87Nd/Pbnrozc92tvewYnb17Hjt37eeXf/Xfiuf/XmY/jdU/bCsBLP/ZD3vKs43nbcx8JwGs/+ZPISvElTziaj/3aKQC878vX8t2bdqfKi8dspcS1H3pRqn0vueUB3vSpy1rfv3XOs3jMUbPcdP8Cr/pEu3z/6Y2n8tzHHNH6vm9xmZf+7Q9pum2kPznzcbzezevLPvZDfuMZ2/jd550AwK/9Q3deX/mkTXz01ScD8HtfvIapYoFPvmF7xz5/9u83ctOuBb7+tmd0bP+7/9rBx793GwBrp0tc98dOXr942d38z6/fkCrfcXh53XtgiRf+v0v4y199ImeevInqcoOX/e2PWApp3Pp549O38qGXPS70t9//12v4wa17AHjq8Ru44Oynhe6XN9fu3Mcvn9v9vkbxZ694PL/2lC2oKmd+/EccXOo2jn7OOnUzf/6qJ3Rse9e/XssltzwAhNdJ/+87t3LJrQ/wH7/3Cx3bP/Xfd/Lx/9zBNR98IYWAoQ3j0tsf5LWf/Enk716d9MBCjRf+v0v4m7OexC894ejW70v1Jqf/9SW850WP5o3P2AbAr557KQ8e6BRcr3/acfyJW/9G8Ydfu55vXn9/x7bHHLWWb53z7NjjvnjZ3fzZRTdx1R+9oEvEDYOxMeAisgb4MnCOqs77f1PV84DzALZv354sz1Litcr2uYrsXreyuu/hKmyhZYD/+GWP4/HHzoWe46zzfszuhU612mwquxdq/Or2Y/nBrXvYe7Ct+B4+uMSzTtjI/3jBozqO+X/fuaWr5Qiw9+AyR89VeMq2DXzrhvYDtmf/EgvVOkuNJtOlIg+71/jC2U+lXGqPTbxu5z4+eOEN7F6ocdxhq9k176iQvz7rZDZvWMXOvYu84/NXsdf1Muxzz/P7L3gUzzhhY+z98/jdf7mKXb57sHuhxhknHcVbnn18zFHwif+6jWt3Ptz6vmu+xpYNq/irs06OPe47P9vFJ/7rNvZX66yZLrXu7ydfv50Na6Za+9255wC/98VrWnne7f7905ef1FLvfppN5ZfPvbR1vnqjycMHl1vPhaqya77Kq045ltc+tbMh+Qdfua5VBgB7Dy7xxGPn+GCEQQjyxcvu5oLL7qbR1C6VEYaXp3OefwJ/9R+3svfAcuu6AO983gn89Xdvbe3n8eD+Gk3tPq7RVPYeXOZeV6WqKrsXqrzylE38+lOPA+ADX72eXQvt89378CLTpe5xsHv219iz0O3B2XtwmbmZMi9+/NF8/qc/p95oUioWWvf7S299WqqKP8j9+6r8zueubD/Di8vUm8q9DztpffjgMkuNJr/9nEfwghOPDD3H7/zzlbGetocPLnPKlnWsmiq1nodh4JXfX/7KE9l2eLTKW643efV5P27dg1q9ycGlBq99yhZe9eRjQ4/5H1+4urW/n93zVZ6ybQNTpUJonbRnf409+7vL996HF1mo1Vmo1RMVL7SfVa8u8gjWSQ8fXKLRVO7b13nfH15c4uBSgz37nf1Ulb0Hl/iVJx/La57ivJ+/+y9XddS/cWl57NGz/O9XOIb+7753G9fd83DicQ/sX2JxucHDB5c5am6FGnARKeMY78+p6leGdV3PdeQZ8Ki/p2xZH2nA182UW0bPY6FWRxUedeRarrl7H7V6uxVcqzc5eq7CKVvWdxxz3GGruP6ebpdrrd5gtlJm84ZV1OpNVBURoea6nWp1x4DXlhuIOC1mz10MUHIrRC+NXp6eevxhHDlb4XDX/Vlz3c/e30cftbYrjVFsXDvdOq+qsm9xmeMPX514/LaNq/nBrQ+0vu9bXObI2enE4+544ICb1kbH3ycft77l3gVaeQuW55OPW89jj+424OAo+/Z5mx3HHVxqUG8qjzpyTVca16+a6irno2a7yzmKn9z+EOCoCs+FHYeXpu3HbXCv15nmU7du6NgveNwTj11HuShd99DrXlhcbrDcUE44ov0cHD1XYbfPgO9bXA414LV6s5WOju3LTdZMl9h62KrWfqVigVq9QakgbHfTnBWvYg8+w8FyP+mYucjyWFsphaa5nacGm9atYcOaKW68bz5yv7zx0n7atg2x7lyv2yl4D7ZtjH4P182UQ/O8b3GZx2+aY7pciKyTasvhx4HzDKUx4N45nrR5PVsOa+ctqk4KPsves+o9u/Wm0lSnLvXyvLZSatWVsWmpN9m4Zqp13KZ1FS6/K7lLx7v2vsXlri6FYTDyaWTiWJt/BG5U1Y8O89pBBR71N+5hnJspRz5YszNlKuVCx0tSqzdDXS1zM2Xmq/Wu/l9n/wKVchFVWG4EXtTl9kM+XSp0GO+4PHrbvQo4WJFncQf578GBpQaNpqZ6gedmylSXmx0vQZrjvLQFX/Bgmr0G2nyG8qyUCh331H9c3PFpyzk6T53lkMS+xWVKBWH96nJHWr20r19dpliQSAM+O1OmUip2HReX1+Czvm9xuev8Xh7C8lGrN5h2n+VgmsMaAmmpuH2+7UZt+3mKykvXOcrFBAPuvIfePUgzTiMP/OUVh4gwXSpkeo+dhn+4IZ6bKcfWSbV6o2t78H4n4d3v6XJn2UfVSVHPcrAe8MYAOOeOL9dWWpabHcdVyuH3Juy4sLQNi5EbcOAZwOuAXxSRq93Pi4dxYe+Fnq+2W47Q/SDOzkQ7KmZDDHjruEqZ6VKRqq8FWF1uhFZWs5UyjaZ2DVZx9i+2jqm6D7N3Tv9f/wPoP28wT1PFQut80+4LXnUfRO9vlgp1tlLqundRg4GCxwXTlua41r0I3IOpQJrXTJUQCSvP6GtMl4td97grbyHHpy3n6Dx1lkMS+xaX3Qaid1x3xT1bKUU+m3MzJabLhfY9jDB6/md/tlJqeXKqyw2W6k0WqnUazc6KvLrcDK38qm4l2VV+9UZf/YeeAajWO5/h7nKLfo+nS4WO8utOu/N+zVbK1Jua2LecF/OLy4g4YwaSmPY3PlO8x9PlQqvcPZYbjut9dqYcWyc1fWLCI6sB9+53JVBvRdVJUc9ysB6o+BoElYRybaXFbVy20lBy7k1SQy2qcTEsRm7AVfWHqiqq+gRVPdn9XDSMa8+Ui5SL0vXgBf+ujTEqYQrc3+KfDlFmYYY2qJQ79i+3DW5QHfr/hjYMAueddyt+T6lHtXane1TgXgWfRkl7aZsPpC2J6ZZabee9XJSuvuNCQZitlDvKs1gQVse4qKd7VOCh5VzOYsCzK/C5mbLvuO6KO8k7NB2nwEPKcW6mzEKtTrOpHeddqAaf2QZLjSbNZlC5OY2arvLrU4FPtxS4V27ZFXiw/IJ45Rn1ng6KfYvLrJ0upRob4KjN9O9xmAL336u4Osl/jeCxk6vAfQY84PGMOy4sbcNi5AZ8lIh0VvDzVael6Tcoayul2EFFjpspvJKcm/EUuFPIy40mjaZ2tBD953GODbZ2XdUSUFphf8NUTLEgrJ0u+Qx4nTmfEmmroc7Wblgao5ibKTPvuhXTVJb+48B5+BtNTT34Jag6q8uNrlZ8MG3edeZ8jZeoc7eVnKvA3eci1oUeqsCzuNCzKfD5ar1DgdcCSrpSLrZcoMHjvDyEKfCWN8q3n8fsTBlVWKjWO6bihT2zQFfFWVt2u4NKgfLL2N0QpFgQykXxeU6aHXlJ13VSTFTg3j31n3PQzFfrXdO8oqiUC93vcUzDqBKiwMMMeFT5Bp9Vf72ZBu9+BxtvUXVS8LxeuoJ14XQPCrwW8AIFPZ6ReQiMHRk2uRpwcfh1Efkj9/sWETktz2vkTYd6DPxNMxgjVoGv8hR45wCjbAq80anAQ5S39zdKxcyGGDGPqL6zLMZnbqbtVkzbZ+cd56XJU3FpDHiYNyJK7QbLN+n8jgLvHsSW1DgZRwUe1b1TKRdcV3a3Ap9fXO5Q2EEF7p0j2BfuJ/gs+bc7jdFgmhtd3R9Z8avJKAUe50mLU+Cq2nq/2kZteAo8zTsB7j0YigLvLl//+5FFgU+FjNuJqpMSFfhyd/06XS6GxofoSkuIAvefM+64sLQNi7wV+N8BTwNe435fAD6e8zVyJWjcgn+TXp7ZmXJXP6D/Jaj4XpKwPhr/efzHetSWm1RKxVDV6fwer8C9dMTlyT9go1cF7p17PsbIxR2XRbm37oWvrzqqwRHMe1LDYjpEgTeayoGlRof7OSxN3v7NprJUb0Z6BWLzlFaBu+UY9Vy0FHiw0jvYLv+KX4F7aVfYv1Tvy4BHKbTqgBR4Ky8BBe5viCd50uIU+FKjiSojUeBZDHgvCjzYyPI/43F1EnR6WKrLzZa7OUsfeFT6wuqkfQFPQFcfeL27fs3SB+5/BistL0A6BX6oGPCnqOrbgCqAqu4FpuIPGS1hCtzvekujwKGzH9Df15pVgQcr3DwUuN/NH5anjtZuy62VTYF75/auk8bt53fRZTHgvSrw+Wqyiz5MgUO7cRI1oMi5h840Py9gSCYF3sMo9LmZEqWCUJDO56AgzvTBKO+QfwZC8DkCx8iHqVb/oM+8FXg/feBOXroVuOdNSONJi1Pg7fd2NH3go1bgYXUSdBq3uOchCue9DU9fWJ3kddMFrxOvwOPHNoDTQF9uaLgCTzi2FuHeHxZ5G/BlESkCCiAih+MsUDK2xLnQsxjw4APs9bVWfH3gYX00rfOsCq8YqgEFXgsopqQ+8LA8dhlwf8u9NSWrNwW+b3GZgjgjwJOYDVF0aQx/UIHXYvrAHTdye2xDUnmG9YFD26jNVsqhA4q8aX5LjWbk6NrY62YYhe65K71nzD8C3nsORNoGPFjptRV4set58vYJU63+ZzSqwlbVRAU+HaLAswyaDGM6RIH7vQnJDbf2vQjSMgyTqMBj3uOwUeh+D1pcnQTdDdyw/+OIa7iF1UnOu9V9zfZ+3fVrcHZIaDpCptytVAX+N8BXgSNE5H8DPwQ+nPM1ciVUgfv6AZOmNQWnQnn/e9s7FHhIC9EjOOXJI58+cGcQW7OpzFe73cgVf8vdTeNUMcM0soAhXhth5IKUiwVWTRU7DEKWaWRpFPjsTKljgJ1XLnHnjlLgjgs+/Hh/+USNro29bgYF7s219+6Vfw66/zmYnSm33P/+fHjHRSnw+cXl0GffPyUxqsL2hyuNUuCVEAUe5+pNQ5gCh86GV+zx5ULLUATxD7Za4z4/w+wDT/NOQIQCj2lEVkpFlhsa2v03O1OKrZMgHwUeJTrC6qSo6wT38zecg/EZQtMRMuUuqwI/JAy4qn4OeA+O0b4PeLmqfjHPa+SN10+4uOTMa92weqqz5Z6gCOMUODgPhfeS1EJaiB7BKU/Q6dppTZUJuOP935PcyF6EuDAF7q/8SgVpLQKRhqACT6sY/GnL5kJPHwhkbqbcarmnU2KF1qCXWqCCijve36iIa6jFXdc7PongvQoaL++6yc9mMdSAR+U1WM5rpktMFQsd5w8O5PPjNS7a076c+7uUhwIvdTeU4/LSfbwz2Cls3q/fhV4sCGtD5tcPAm+ufZoBodDbPHCgY5DXvsX2IMe4OgnCFfiG1VOpGzfetMKotIU2yHznbkViC+wXVOD1plKPiYMfOv0s5aDSqChxwyLvUeifVdWbVPXjqvoxVb1RRD6b5zXyZnamRFPhXjcc4+b1zspYDyzUqNWbyS70EDeTfz5zO+pUwzewJNnV7R3jncNTLUF3vP973Hmry81W/OIwBe53P2YdUOSfz92rAfemhGTpA08TCMQ73737FlNFiPMPRvMrMq/fN9KA+waThQ2mSSLYLRBHcI52cACXd92WwT0Y9Wx2D2IDJ69h3Q2rpoqU3OhuznTEcsu74xGcSufhuNadcuoOvNK/Aq/4Knz/PUwqN//xEK64/AMDIXzmySDI0qgFr/unezBj5P4hbuLgvYqqk6C7gQtO/Zl+EFu8Ag92ifiv4/+/S4GXOxW4s0+0AQ8NAJNyFLp37CFhwIHH+b+4/eFPzvkaueI9rF7Qfi/esPc9qfUbNlfSP1iqQ5nFKHDvXP455f5W9HS5PwXuz1OSAs86oGjttONq68WAe1Od9i06EeLSGL1CQZgqFlIrcIjOe5AOt7LfqKVV4PXBK/D5ajYF7u0fnGufVYH7+9X3uX3ksxHPbPB/L051mALPOuUujOkIl2ur3BI8acFgMH78ChzaYY8HTZYZHRAehChegXe7icMMeGT5hijwzRtWpb43vSjwzvgDEQq81KnAg2ntTkf/CjwYC2RY5GLAReT9IrIAPEFE5kVkwf2+G/h6HtcYFMEKfkvAgPcziA06lVVWBe4PytFuLTdb05S8787f6KlUswlGrF8F7ne19aLAW32uM/FBVvwEg5AkKfC05elX4N1GLXoUe2WYCrzVT5lSgbv7B+faRynwuHJse0x88bIjFJo/L35FGFRFcd6jtPgHcPnTsDulJy0uFn1wkNNYK/CAZy7OgLe9eukVeLVDgXcb8GPXr+oaOBlFPwp8udFsje0I7hemwJOC9Pj39Z8jaVCp/36nnUGSJ7kYcFX9sKquBT6iqrOqutb9HKaq78/jGoPCqwR/HjDgP09Z4QfDsfpHCEN2Bd7hropQ4GEDhdIo8Kg89avA/WlPGw417Li5mFjVQYJBSJIUeNryDCpwkXR561eBe4MGB9kHHnpcQIGvW+VriIWoVi9uwj6fAe90oYcrcP/AvqliAZHBKvAN7qp0aT1pcbHog/3J42rAOwclJo9lCVOn84FGalSdBIGuCtcjs35V98DJKLIocK88/QONwelz7xpjNCQF7gX4CaZtmOQ9iO39IrJeRE4TkWd7nzyvkTfBCj6rAfe7FaF7NS6/sgrro/HjDyrjHePt71fgnf2MzcTAIYkGPNDa7WVAkd+12usgtizHOX2e6RX4z1NW5JVye9BL1b2nczNlds1XWWpEKzn/NL+4gD1RFArCVKl7Wk8YLdeqa2CDU6g8YxiMNR80CJVywR2U1GxN6Vk3U2b3fLRqDZZXXB9pqAIvFVuRtqr1JvVGk3pTMzV2wuiYQuVW+MWCpH6PDwUF7sX99qbyJXnSgmNJvGv635GoOgm6uyq858Gf9jiqMQ3vYJ0UtTTwEWu9pUcdD6cXAyEuj93p6E2Be4bfS8Mo5oLnPYjtN4FLgIuBP3b/fijPa+RNu4J3B7G1DPhix+9xhEVzi1XgMWrR737yt/xLxUJrrerOkb6NxMAhwTx2udBzUuD3z9dYbqRbStR/3MGlBnv21zId16mU0yjwdOXZqaQbrQUsko4PnUaW0Sj5+zDj2LfYOde+U4G3F8vxxiZ0KXDP8PvUiT9U6N17nbyGNXYSDXgKBd5Oc/tZztLYCSOowCvlArOVUoZyj66wg+7osTXgHc9g8nscXNYVumMlRNVJEO56Dxs4GUUw/nhn2jrrpJmpYseaDi0DPltppcs7n78bLiyP3enoTYG3DLibholX4MA7gVOBu1T1ucCTgAdyvkauBPtIj5qrUCxI2/WWMG/Y2ad7Na7uUehNX2UVpcBLLDeUxUAfrLe/t1Z1sJIMm//Yed52HosFYVVgNa6gK7WXynR2puRzV6Z3hXv3d+fexUyu94q78pLjxoquCLxIYmldqcHy8hR4+3lIUOD19lrYWe9j0prUHs4AsvZc+85of43WdQuBhWxafee++eOAa0idezjbkdfucpydKfHg/iUWlxsdfeDeymNR08iC3Qre/NxeGztBgnPhu8st/pnMosBnZ8os1Zuxqi4PvPJam6IOgpBnN6UC9/LnDXL0P+NRdVLw/y4DnsKYxcWuCKuTwoRSW4E3Q8/XFlBxhrj7fU0zqNQ7zkvDoWDAq6paBRCRaVW9CXh0ztfIldVTJQoC+2t1Vk8VKRedlru3Bm7aRTmi+xnbLpyo1Xf85/GfI7i/t1Z1h2uy3gyNQBR23v3uCOTgQDH/gLBaxlW0/Nfw7lkmBb6qM21pcdZvduIve6Obw/BWY9tfq6daVzlYXp4CT8qbfwWl6nI/CjzdIDZ/OvzR/hxvRPu6c6vink2ve6fZcr0n5bXj91VOhe3FTYDoaWTBZ9SLkNVLd0MY3vm86Wppy81/vJPmZAUeFSM8b7y59mljMviNVZr16INu4rBR71F1EnQ2drxnMsu9iYseGVYnzYbUs57x9J6l4Pla0ztjp5F1v6+lYoFSQWK7tDzjfigZ8J0isg74GvAdEfk6cG/O18iVQkFaD91c4K9n0JPwT7UIm+ID7RZisI8meB5oT0kLqpPpUAXe6FI3QbyIZ/5r+PFau/5Vl7IyG/LSpyGsskiD5zJNE3HKS1tUGNSO85b9bkjnXqTJ23SIAs96H/3dAnEE3ZzT5c7BS/6GnH+UeHCu/XRAgU+XionlEfzdU2utKT2JCrzQ+punAp8uFWiqM13N60bI8kzGuUyDccWzqMx+CA4oSyLYjZO0wlswz2Eu+6g6CbpHh89WfPHTE6ZVJdU1YXXS3Eypa7nbIz0Xeg4KPGxd8iwKPLjs6jBI7+tMgaq+wv33QyLyPWAO+Gae1xgEczNlHj64HGnI0xwfpXL80xj8caqjzuM/R1CdeHG6oxR4nIrx+prDPAr+Fz+uVRxHr4a41+Mq5QJ79tdTxXyemylzz8OLqc7vXykruAZ0XBorIQo86330TwOKI6jA/crdi53vT6//2fTPtferTm/6WRYDPjtTplZqts597PpkBd7qDip3KvB+FzPpmMa37HQj+LuKsnSdBGl3UXVOzxv0oKXggLIkuu9Bkgu9U4HHGfBUCnxVO3560r3xVniLW8zEuUbTV2/CHXsOdJz/8IACDzYEW/VvCgUe7IL0B8aJO+6Q6QP3R11T1e+r6oXA+XleYxAEDfZs4G+a4/0rH/mPDRsolJSOfQE106nA24q7WJBUCjwsj36CfWe9DmIL+z/LcVkqqywKPEuDLFyBt9u546DAu1zoZX9/YbcC9xtw/1x7f/x1T4EnqdZgOUc9s86Ay/QKvN/lRLvLrZ2XNJ60uNHK1XqDom9K1vAUeNapld3PbhzBfv/gIEeIrpOKBemIQeHNWljjdkkm3ZukQDNhdVLwWZ4uFVoeoHa/f1BFdwYNCk9Lfwp89XSptabDsMnbhT5xkdigu4LvRYH746f7+1rDFHhSOqIUuLdWtbd9bqacWoEHvQt+gn1nw1TgvbrevWlDrf6rBAWe9vxJCjyqkRFU4FPFQqoFXTrOkVqB1zvSMV0qsNRostxwxgR0K3D/amxtg1BJUOBhee2YYlTpXnLS/2zWQtTawBR4V7kVspV7ggL3h3odlgHPPrWyfQ9qPSjwYPef//9gnTQ3U26VqXfc7IzTRbW2kjxKv70yY4ICD7yH/sHCczPlrvo12JAPhu0NT0t4Y8K/MmHccZVA42KY5GLAQyKxzU9KJDboNm69GHAgdMnJgShwd7tXSWZT4N0t+rwV+NqUqydB56juQfWB96vAvePWTncur+mnc5pfb1Px0ihw1e71rb3y84yovzHjjdwNBhjqzGt3H3iUau1S4KvCn9nZSilbH/gAFHjQoxZ7fIIC96dvXA14VgUeXAEvjQvdX/dEDX5LY8yC3RJBohS4F/HMuzdd9Wu52wg714tX4FOlQlfX5lTCoFL/2IiJNuAhkdhmJyUSG0Qb7rQuXc/F6s2P9btcOyrJhOAKayvBl6XTteO5Sr3t3nSWNNOWPEMZNg3KP5glzfSTuPPHGbkwKuVie3RvBsMfnIYU730odfyNw+9yy2oIWtP86r0Fw/GmxsVRXW6y1Gh2PmPu/fMG9vgrxdlK52ps/jxU/JVfQIFH5TXY4Oqu4J0IYKunS4l94N5ce38eesUfy9xrQGV5j7Mo8LAlhAeBJwbS0jmVMcM0spg+8Kg6adanwIPHpTLgCQ23sDrJ/6x5BrwSqF+7FLivUROZlkD5eqRW4O4Ut4k14D7+TURWA4jIr4vIR0XkuJyvkTst4xZ44dO2fuNW4+qYXpSgzLwpT213ZLdqcVxF/lZwg7BpEEHS9IHvrzlR5PpR4Fn6sbvSlrDghJ/pLhdsct7TVeRtl5t/OlKa41vT/FJM4Qk9vtSOJhZFWCXrld++lgIPV4tRCtyL3+5vrEQ9+14ZzZSLTJUKrJ4qUiy0Qwl70bWC3oQwBe5cN98+8EX3fZguJ+fFz1SsAu9skJWKTr4HWWEv1ZutufZpaZdnM9Uz2I6I1zbEU8VCx3FRdZLfw9KLAW9H5otX4P46Kaye7RiI6YuB4BEM2xtGLeBh8Uia1un3/nnjoIZN3gb8E8BBEXkizrrgdwGfiTtARM4Xkd0icn3OaUlNHn3gEL4Mo4gTIjOqhRjEH6wg6NqZLhU6FLfTD5Uc4S0pT95xbQPQuwHPUuH0c2zLBZtzH3icAk8aUNSa5hfiykuDP5pYFGH9lF3lF9JfO19dZr4a3rjMosC9QUrefiLCbKXkm2bk9FcGvQnBrg5v6lt+Ctw53pv37Vfgacq9WBDKRYlQ4N3GcNAu01Y5Z2nUljoVeJpn0D9Qy4v1H3QlB+uk6VLBHcPQacA9z0QaY5ZWgXvn9oIMOdvqrWc5SYF7jZQkBR72/HmhaeOOc9JWOGQMeF2dmHtnAn+tqn8NrE045lPA6TmnIxN5GfComN6ea7UaGCEcda6Wuyrg2mkP/Gm3gv3f4wfItV+uIGEvS1ayei0601amWBBWT6W/bqVcRBUWap7bOHkeeJY+cE+BZxkM5bncaj2uruWP5x1FmAL3Kq2w8vP2e/hgd+MyGL/dP3I7Kq9e3IRgP2mYAvfnxYtTXS62x4b4FXi/i5kEvRCVjAocnGcoajnRoKFxjNrg5v1mDaMK/gGz7rOb4hmc9jW09kWMeg/WSW0Pi3tcYH16x50cf2+SGm5hosI/YNLrDvIr8OAMjPa54geHRq2lUCkV4o/zef9mUwzcGwS5zgMHFkTk/cDrgGe5o9Bjn0BVvUREtuacjkzkZcAvuWUPuxdqnLZtQ8fv3ktSW25y2Or4l6rjZQm4dtoDf/z9UCkV+KpkF3qYgkuL52rr1YDPhbT848jiNehXgXv9gEnH+wcZ9qTA/ZXpwWUuvuF+6s3OZRlv2bXQlRav4o5T4F+96h6aGq7c/QOFvPjpcXmdizHgnrELehO8QXL+KWx+Bd63Cz3EC5H1PfYvCuMnzB0dpbh2L1TZPV/jpE1zADy4v8Z3fraLZszqmoetmeJFjzsKcO7Tv197Hzt27weyT60EMj2D/pjjUYPmwuokf/x9z1j7G3/7Fpf4l5/8vOtcG1ZPcfpJRyVOHwyrk7y0XXzD/SxUnSA3XWOMQhot/jC7t+5aoFIutta8gHgFvhSnwH1T4eZmyhxYarDcaKYK/pUXeRvwVwO/BvyGqt4vIluAj/R7UhE5GzgbYMuWLf2erotHHrGGqWKBbRtXA7B142qmSgUeecSaVMevmS5xxNpp/v26+wBa5/FoTXkK6aMJMjdT5rYHnJfXv7azc562Ap8qFZhxV8464KnQmErwhCPWMlUqcFwgbdA2fvN9KHCAE45cw6OOTHfPOo9b23J9psVr2KRJ8/EbnfJNU57+aSlL7rzSYkF45BFreNSR8c4k/zS/nhS4q1pVlS9c/nP+7KKbQvebKhY4Zt1M+7ql6Huxaf0MU6UCn/+pU5luPaxd/t5+i0uNluEtFIQTjogvxxOOWMuRs9Ot7/4BPJ6xC3oTup7lUrEjnnhegVz892DjmmnWryqnfianYxR40JDOzZS568GDXfv+3fdu49+vu4/LPvB8AD7933fyN/+5I/HaP3zvczl2/Sp+cMsefu+L1wBQkPbqiGnw7q+3ImIqBe5b9Wvf4nJr1S8/YXVSxRd/f9/icsesheM3rma5ofzBV68LveYP3vPcxHIPq5OOnK2weqrIBZfd7Vzn8NXdY4wSFPjvffEaNq2b4dzXtWc3967A/Qa81ErvYSH3cFDkHYntfhH5HHCqiLwE+KmqxvaBpzzvecB5ANu3b09eKT4jjz5qLTf9r9NbU782rZvhpj85PfU8XhHh++9+LvNVZw548CVoTXlK0Qfe1dothStwz40F7dHHcZXgSZvmIvNUKfWvwAG+/NtP7+m4956ePVx+XL9vkC2Hreoo3zi8QS/te+rcm2+f82ySHAR+Bb56dfZXy2uULDWaPLh/ialSgR+857ld+62aKnZM1YtT4BvXTHPFHz6fg0sNysVCa+1icEL6FqTdDeEd9613xuf1H17fGdphbqbMPe4KZkkKvJ3X9rMr0l4PvVfCnoeZqSJX/OELEsvNn6YoBe5F/PKI6gN/YKHGg/trNJtKoSDsObDEhtVTfPOdzwq95g9u3cO7/vUaHjqwxLHrV/HggRoAX3/bM9i6cXXGcSHBwYzZFfgjD+9u7ITVSZ0KvFO5/+qpm3nuY46gqZ1V9Y92OI2Thw4sJSvwkDppzXSJn3zg+Ryo1SkWhI1uPTtVKrC4VO+KgdCRR/d6DyzUuuqKaAWe0Hdeb6+57p9OObEGXER+FUdx/xcgwN+KyLtV9Ut5XmcQBCv3rEE4ZqaKzET04XpqpJZGga/q7E8MKnBPcXsDhcBxt6YJHBL1u/eie31Zvc7JzeIC7/c4f97936PI0hibLhXa98J9sdMcXyk7I5N7VeB+NeFVil6s5/jj4u/F2ko5dG6+iDhpDhyXlNdgec0FFHglpQL30jwdMgc3K8HnIUu5+dMUpsCXQqZkRRnwfYvLraBOXr/oulXR5egp7GAo5uMPX50pngK0B+KlfSegU4FHxV4Pq5MqvsbOfLU75GuwwQOdeU2rwIN10prpEmsCCxJVSoXWeIQkBb5vcblrdbdqPTzEdCWp79w3LTg4nXJY5O1C/wBwqqruBhCRw4H/AMbegA+SLArcm55R9S0w0T6P51aqdyhwL6xgr+SlwIdJFgWe/dztKUJZuhPy6AMHOgJVpCFOgSdes1To+x56xsxbeMLpm0ynwPctLvfd/w35DMT0L8vqJywoyuxMmUW3m8W/aIg/Uthspdw1cDBIsOLft+iu9Z6wal5kHnzPbpry9BR4s6mhhhjC66TpUpHlhtJodgcIisKf114UeBTTZd/7Gjqf21HS3tS8oJHtR4F7x43KgOddUxc84+3yYNI1ROTzwKXAo0Vkp4i8Oec0jRzvJQmL1RvEP+2nFqLAvd/8Cny+utxXJKu2O7O/PvBh4s+7ozzye5Qr5ULrXmQxahV30EvUYJrE432BNZwgHukq8VYfeA/lVykX+y732Zky9aazZnTNU+C+/nyIVuDz1f4an63zlTvvQS/njFLg3mwEP1EVdlBJJwVj8QeBau2fcUBnRx58z25aBV6rN1mo1VENH/AXVif5p29508+SmO0w4Cn7wFPkpeN9DZ3P7eQxWDYe0X3gjsez3gg34mEK3Ot6GxZ5K/BvicjFwOfd768GLoo7QFVfk3Maxo7pUpG9B5ZZaqSbBw7OYIigayeouPNS4F19ZxOmwPNObz8KvBUUpU8FPl8NH1AUftzoFbh3/VYfuK8/3/NARSnwfpcShXZ/fr8KfO+Bpa7t3mwEP/48+93F3vXnfcbCP3AwSLcCz7aEaFceMipw75kNLsLkJ6xOmg40Nk/KqMAb7rD86GlkGRR4Qp6nSwX2HlzqmCnhX/MhToF7+4etye5X4LOVCVfg4jQZ/wb4e+AJwBOB81T1vXldY1KZLvlbiOkU+L7F5a4Hq7PSK3R+72MebavvrPUSjL8CH6wB782o+YPL9BqJDTr7wNPQ7W7M0ujIVuGH0fnMdg6wbK+S1tmt4K+gk7xSaXDGLvjy0lMgnfB5+FHTyKBz3etmU1modqq8JBf6TLlIuSitPtyk/dPkIctz4IW0jZt3HlYntRub6Z/Vihu9b95V4N4AsDCy1En+PIcq6UAeobPcouLG++Osh+Ffc93fyBkmudV8bgCXr6nqV1T191T1f6jqV/M6/yRTSeij8eN/WYKunUpHpVfs/N6n0e1UneOvwP3zRPN2+fvLK4shaE3zi3DJpTke3D7wg+kr8i4FninN8ZVfGlrP7MFlqm43kVfBe4OAqoHxH5WcFbh3zn4aoWGx6OuNJvWmdt2b2RAX+v6lemu+tzcmYL4ar6hFpGNAXJaGWxjTGZ9dT4GnNeDes+2V3/7aMgeX0od89fJaTdHITVsndbyvEQrc72WATkMbtQKjf0ppGP7jvDUdJlaBu/xYRE7N+ZwTz3SpwII3LSmhkkyjwBeq9Q4FvlCt5xDJypfGCVLg3r3I+9zevcjqQj+w1EC1NzXrHbO45PRJpjbgvnvh/57u2GJPx/npVuDtBWq8PuVglCzvGVuo1nNrMPrz0ss5wxR41LrVcyGKyxsxDc692F+r02hqYjn6Q5Wm7U+Oy0OW97hnBe6ee/d8LfK4MDwDXkvRyE1bJyW9r9MhCtz73xt42asCDwZN8j8DwyDvPvDnAr8lIncBB3CmkqmqPiHn60wUHSo6rQv94HJXy7DSoWCKnYomBwWeNo3jQOc9zV+Be2QdxNZPmrxj9uyvoZo+CpcXb3/JrVAyRbQLGSSZlU6F5g1yansTgK6BfZUQY94v/Z4zTIFHjZYOG8QWdNHGGcXgufJS4J2DXnNW4L46yTv37oXsBny+uszq6VIqBZ4mL0nva0uBV7vLKi4me5ICry03WLeqHVchamrhIMnbgJ+R8/kOCTpUdMpBbPsW610tvE4FU+issPo0ulnSOA50pjd/Be6RVYH3kybvmF3zVSD74i5LEUoi/rjugWVZ8Z7Zhw4suStHtRV4NUGBe2nPg34bod4YBj9Ro6Vbo8cPhrtl9y22DXhSQ2xupsyD+5dC12zPnodi6P/R+xcSFXhYndRS4AvZntW5mTK7F6psWD2d6I1MWycl7edNBwt6SCDaw+I/VyYFPuEu9KOBh1T1LlW9C3gIOCrna0wc0xkUeLlYYJW7VGFwEYXpQSpwL4iHb8GJcaaj1T2GCrzX9cChrWqyuFL9fXHZrtk9tSsrXvz0B9x0hynwaowCz8uD4j9nL5HdvDEM6osgFrVQ0HSp2NHnDkE1XvcZ8Hid5AV8ObjUoN7UTGuAd+chowL3vD4LNYoFYVVIMKqwOqmlwOe9ZzWdFpytlDoGO6ZJW1KdlOTh9ML2PhzS3REXi386QYEHPaSed2GY5G3APwHs930/4G5b0WRVt3MzZfbsr3UdWwko8OkBKHD/ghPjzLAUeLYR3fko8KxuSf+xWZ+DPBR4wV3Ixkt3lj5wf9r7xTtn3Ojm+OMLNJWOBWTi5isHFZf3/+Frp9m3uNxS5Gld6Gld7vF5yK7AwXnm4hYUCtZJ3rl39eBC9wY7plXgSXVSGgUOTgPTC7+6r7X8rSlwP6K+5quqNsnfTT9xZK3Y52bKLTdq1MPpD+SS9rxxeI2DfhsCw6JUdBYZgfxd/tMZVUzYcb2kqVWZ9uBC956FrNftN80ec6vaz2ynAm/SaCrLDY18D/LyoHh56VXRt5embCuu1oIVIc9BlAHfvH7GNeD11n5xzM2UWagu8/DBPAx4tka9d692zVcTI8a1y9evwJ1tab1FczNlFmp1qkuNxBk5aeukKJHT2lZq5/Gw1VMtbwL4GmihAWC8RmhEH3hgauTsIWDAbxeRd4hI2f28E7g952tMHFkHN83OlH3uyHCF5J9rm/a8cXiVV7+u+GHiVQB5D7rLqmLa6emvD9Yrwwf6UOBZr9tvmj3mfM9sZx94o1VJRr0HeSvwXs/Xji7WVlwtF2vIcxB0mXpRAY9eN9Nas9rbL465mTJNhfv2LabaPz4P2brVvHv1wEIt1gjPdpRvW4FnfVZnZ8qoOgM1kxV4ujqp81mKV+DeUrheubW6SEKeGe+81QgFHlzzYHamzEK13gpSMwzyNuBvBZ4O3APsBJ6CuwzoSmYQCny6XGytnJX2vHFMmgKHdqt5kAo805SsnBR4L4PYRq7AA8+sP5CL50aPeg/y6gMfiAKvJynwduhML/ztOp9LvFiQxLjmXjn//KGDHd97y4OTTkk5lqUXBe4v313zVWfef8pnx7vGrvnqQBR42L5+BT7rGvBBKHAvbwtD7AfPeznR3cBZeZ7zUCCrAvcWh+8+trMC9FbO8sfk7RVT4P7zOvcgzQpvYcf1mqZSsUCpIBxYcqJUhQ0oiqJfBV4uSqtLoheCz2xLvbiBbbztHoWCMFUsuKFW8ym/Sp8KfDpGgYcZqNmZMjfet9D67oVB9fdpz1ZKiWNKZnM04K3BjCnHsnj36kBCMJZg+XoG78BSo2Nt+CS8axxYauSmwKcTvDleuXp5nK92rq7mXCObAl9uNLvWXPdPLfRPLxskuRhwEXmPqv4fEflboMt/oKrvyOM6k0ovCjxsf09xq3aqzzRRjdKm0RR4+x5kHhCWg5qdLhWouxVNpvncvgE/ma5X7u24IMFnNkmBe9+XGr2t3BZGu9z6U+D+BU3a88DT9YF7BrzRVO7bt5hpla67c1Tgae+p/17NxYwkjyrf4G9JRNVtoWlLmRdvv6mIGAjTASM7Xy237nWsAvcadCEKvBbimRnFimR5KfAb3b+XE2LAVzpZpxf5p5H4j+1Q3D7FtW+xfzekv+U+KfSqOpNoK7mMxjCH/uRKuZiohqKO6+W6eXkxOhbdSaHAvf0WavXcnrl+FXgrcEfdP4gtRoFXyuyv1ak3nMUuvJXE/C7xLAb85w8dRISu9aqz5SHbe+xXnrGrpgXqpJ4N+Cr/eZJc6Nn6wKNc8pWAkZ2vlrk+qMBjXO9ho9CrIdPP2tH5hrciWS4GXFW/4f77M+APgK2+cyvwmTyuM6lkV+DtYulWLa7iDqhPU+B5ntdTpRmNYce0vt4VOMDajAa8dwWezz30V/BZFLiThpwVeK8u9BgFHtUHDs4SkhtWTzG/uMzm9TM+Rb3I9q3rE6/rzaG++6FF1kyXMnXbdOUho/doOsQAhRGskzwxUas3M81b73xOUk4jS6nAo965TgVe6nCht6cJdh9bLgoiCQo8pCEziQrc45+BdwPXAdEroa8wso647Wyldj5YLcUdUE79TsXpdRDUKBl0H3jW83YG1uldgUN2N2rPCjynwYv+9FbKxVZ/fnW5EapW/NfOW4H36o2KU+Bh5/RX2BtWT7UUuOeNWFxupJpeNefb/9g1Mz2lvZWHjGNZguo0irA6yYuj3qsLPa0CT6qTkp59//bZmTKzi3UOLjVYbjRjFbiIUCkVQ/vA4xT4JBvwB1T1wpzPOfH4lUGafs2OfqJyULW4D7UpcN89yDfNvfYLT+egwL3lCbMa8J4VeI/dBUHC+jY9hRYVLKP9LE+AAo8I5AJOCFVVbS0F6r8XacpxzXSJYkFSLXySmIesCjxkEFYYYXWSd0+yRAxcNVWkVBDqbsjd+LRlVOAR5+vqA59pLxwTN0jRu3aoAg/xKrXC606wAf+giHwS+C5Q8zaq6ldyvs5EkdXIdlaGgX7DUiHib14KfIIMeILrrPfz9q/Aex8J7SnwbK9mr+UXfI56JbSCdxcHaRvBYB94vgq8PSaix+6L1ij0doVdW24gEh6a1VOl/jCovRhwEWG2UmJvhiVkI/OQsTxTK/CQOsm7X1nS7C2f+uCBpQwKPMGAZ1DgczNl5n3lFtdA87aH9YGHud69td0n2YC/CXgMUKbtQldgRRvwlqswZcUS52YKxrzu1XUapN85tKMgr7x3n7e3RpF3XLEglHsI5Qlt13vPCjxzLPR87mHQhQ5OXqrLTZ+7MfAs56zA+22EtqLHBRR4lOfM7zL1B22ZzWjAvf1yMeAZY+L7n/E4JR1WJ3nllzXNngFPVOAp66RKVgVe9a2eFzE+o3VuNz5+kLAIfcG13YdB3gb8iar6+JzPOfFkdRXOjkCB56XEhkleee8+b2+GxZvm1889bCvwjAa8bwU+ABe620eapMDzKr9+u4HaA+98CrzejExflAH3FndRzWbA/X97Jet7nHY0eVid1IsC958r9TSylAo8TkV7BA14rd6IjZ2fRYFD59ruwyDv2vrHInJizueceLIO1umcsjEcBd7ralajpO3ZGIwCz+ra9Ubm9nMP+1Xg2Vcjy1eB+93NzpTHRqICz6v8+p0K2Z761q6wnRWnwtMXZcALBWm9w1lCjGbZP4qWFyK1AvcZt1UxBjykTupHgfvTGkXaOinJw+nf7p/m5ynwuPMnKfDgszFsBZ63AX8mcLWI3Cwi14rIdSJybc7XmDiyKjr/PEtT4NHkNQAr8rw9GBZ/HPCeru1WJlmXlBy1AvfmLvvdzYe6Ap8uFZgqFkLjnnt/05ajZ8CzDAiLSpP/bxLebAERWDMV7ZANq5O8+5w1zbkr8IQxRv7zVMrF1vXnXQUed/6sCjwYH3/Q5O1CPz3n8x0SOH2ikkkZOH1iS13hLYOt0qQBHGnJ2nc2DoybAveOHYUCr4xYgZeKBdZMlyj54m9XAgo8aFjHTYGX3RXu0ipwEWmtQNVa+7viGWKnas3qQu/XgPfiSauUixQLkjj/PFgn9a7AS6nSmLZOSlLgXtje2UCjKq0CP1DrDsxSi1Hgd+w5EJvePMk7FvpdWY8RkdOBvwaKwCdV9c/zTNO4MF0qZlIGczNlDi51u24iFXjfsdBNgXeddyQKvDdVM2oFDs4zW2+2jd90uRg70nfcFLh3jrQKHBxjFLb2d9Y+7VH1gXv7rk5YcAW666Re+8DnhqzAvd/8rvvpUqHVB5503EMHMvSBVybbhZ4JESkCHwfOAE4EXnOo9qE70amyKfCwBys4B7Pfua/t9A3GGA6SQTU6+pmO5I9C1tu1exzE1uO9yCuOADiNjuBUupq7nGhYnOq8y689/a/3ZzjoMk2q4D2X6fzickcY1NEZ8Ozvsd+4xRGsk3p9VlsGPHF+d7q8pHmOpsudeZybKTO/WKdWb7ZiL0SlIdyFHh6hb84dxNYc0pKiebvQs3IasENVbwcQkQuAM3FCsh5SOCs0ZVPgYRXRdKnYMWqy3+hT7fQNxh09SAY18K6f6Uj+lZp6u7aramIGFEVd1/83LXlOH5ybKdHwKfBKucj981V+fNuDsas95VV+lRwas5Vykcvu3MtHLr4JgNsfOMAjDl8Tuf/cTJkb7p1nf7XOWl8YVGcxmvRxzfMy4L28x5VyMbUBD66OOFUsZK4z2go8rWs8nVKPe46mS8UuA37Fz/eyVG+yYXX0ymHT5QK756ut58Hjirv2OtcM6QNvKuxfqmcex9ILozbgm4C7fd+9NcQ7EJGzcdcV37Jly3BSljOnbl3PY4+eTb//tg2hL/8Tjp3jKcdvaH1//KY5nrh5HTN9VoLHrlvF8RtX8+ij1vZ1nmFy4jGzPPrItWxck345wzRMlwo8acs6HnfMXOZjn3zc+r4M0uM3zfGkLetiBxSF8YjDV3PcYat4xBHRxiaM6VKBJx+3npM2pX82o3jKtsPYvdCK38RJx8zyzevu44Z75zue2dbvbl7zUuBHzVV4xOGrM71nQR6/aY7/vGk3t+5aaG172cnHRO7/pC3r+cGte3jowBLPPGFja/spW9Zz90OLqeOaP37THMdvXM3Ww1b1nHaAVVMlnrh5HSdleHZP3bqBTeuTQ7gG66QnHLuOXfO1TKvmATx+0zq2bVzNsQnX9OqkxxwVX56lYoHTtm7gpE3ReT5t2wYedWS7bnvycev50hU7AfiFRx0eedzjjpnjG9fcy99///au3x5z1NquNdfnZspMFQvsrw7HgIvq6BYPE5FfAV6kqr/pfn8dcJqq/m7UMdu3b9fLL798WEk0DMMwjFQ0m4oImRs1SYjIFaq6Pbh91Ap8J7DZ9/1Y4N4RpcUwDMMweqafleR6YdQKvATcAjwPuAe4DPg1Vb0h5pgHgMyj3WPYCOzJ8XyjxPIynlhexhPLy3hyqOQlz3wcp6pdvv6RKnBVrYvI24GLcaaRnR9nvN1jojssekBELg9zTUwilpfxxPIynlhexpNDJS/DyMeoXeio6kXARaNOh2EYhmFMEpMzZ8gwDMMwjBZmwOG8UScgRywv44nlZTyxvIwnh0peBp6PkQ5iMwzDMAyjN0yBG4ZhGMYEYgbcMAzDMCaQFW3AReR0d+3yHSLyvlGnJwsisllEviciN4rIDSLyTnf7h0TkHhG52v28eNRpTYOI3OmuH3+1iFzubtsgIt8RkVvdv+tHnc44ROTRvvt+tYjMi8g5k1ImInK+iOwWket92yLLQETe7747N4vIi0aT6nAi8vIREblJRK4Vka+KyDp3+1YRWfSVz7kjS3gIEXmJfKYmsFy+4MvHnSJytbt93Mslqg4e3jujqivygzPv/DbgeGAKuAY4cdTpypD+o4FT3P/X4gTEORH4EPCuUaevh/zcCWwMbPs/wPvc/98H/MWo05khP0XgfuC4SSkT4NnAKcD1SWXgPmvXANPANvddKo46Dwl5eSFQcv//C19etvr3G7dPRF5Cn6lJLJfA738J/NGElEtUHTy0d2YlK/DWSmiqugR4K6FNBKp6n6pe6f6/ANyIszjMocSZwKfd/z8NvHx0ScnM84DbVDXPqIEDRVUvAR4KbI4qgzOBC1S1pqp3ADtw3qmxICwvqvptVa27X3+ME7p57Ikolygmrlw8xAkg/qvA54eaqB6JqYOH9s6sZAMethLaRBpAEdkKPAn4ibvp7a6b8Pxxdzv7UODbInKFu/ocwJGqeh84LwtwxMhSl52z6KyIJrFMILoMJv39+Q3gm77v20TkKhH5vog8a1SJykjYMzXJ5fIsYJeq3urbNhHlEqiDh/bOrGQDHhZ1fuLm1InIGuDLwDmqOg98AngEcDJwH45LahJ4hqqeApwBvE1Enj3qBPWKiEwBLwP+1d00qWUSx8S+PyLyAaAOfM7ddB+wRVWfBPwe8C8i0v/6qoMl6pma2HIBXkNno3ciyiWkDo7cNWRbX2Wzkg34xK+EJiJlnAfnc6r6FQBV3aWqDVVtAv/AGLnP4lDVe92/u4Gv4qR7l4gcDeD+3T26FGbiDOBKVd0Fk1smLlFlMJHvj4i8AXgJ8Fp1OyZdl+aD7v9X4PRNPmp0qUwm5pma1HIpAa8EvuBtm4RyCauDGeI7s5IN+GXACSKyzVVMZwEXjjhNqXH7i/4RuFFVP+rbfrRvt1cA1wePHTdEZLWIrPX+xxlsdD1OebzB3e0NwNdHk8LMdCiJSSwTH1FlcCFwlohMi8g24ATgpyNIX2pE5HTgvcDLVPWgb/vhIlJ0/z8eJy+3jyaV6Yh5piauXFyeD9ykqju9DeNeLlF1MMN8Z0Y9km+UH+DFOCMHbwM+MOr0ZEz7M3HcL9cCV7ufFwOfBa5zt18IHD3qtKbIy/E4ozOvAW7wygI4DPgucKv7d8Oo05oiL6uAB4E537aJKBOcRsd9wDKOWnhzXBkAH3DfnZuBM0ad/hR52YHTB+m9L+e6+77Kfe6uAa4EXjrq9KfIS+QzNWnl4m7/FPDWwL7jXi5RdfDQ3hkLpWoYhmEYE8hKdqEbhmEYxsRiBtwwDMMwJhAz4IZhGIYxgZgBNwzDMIwJxAy4YRiGYUwgZsANYwUjIutE5Hfc/48RkS+NOk2GYaTDppEZxgrGjeH8b6p60qjTYhhGNkqjToBhGCPlz4FHuGsw3wo8VlVPEpE34qyiVAROwom1PQW8DqgBL1bVh0TkEcDHgcOBg8BbVPWmYWfCMFYi5kI3jJXN+3CWPT0ZeHfgt5OAX8OJs/2/gYPqLCxxKfB6d5/zgN9V1ScD7wL+bhiJNgzDFLhhGNF8T511jhdEZB/wDXf7dcAT3FWYng78qxMWGoDp4SfTMFYmZsANw4ii5vu/6fvexKk7CsDDrno3DGPImAvdMFY2C8DaXg5UZ+3jO0TkV8BZnUlEnphn4gzDiMYMuGGsYNRZb/lHInI98JEeTvFa4M0i4q0kd2ae6TMMIxqbRmYYhmEYE4gpcMMwDMOYQMyAG4ZhGMYEYgbcMAzDMCYQM+CGYRiGMYGYATcMwzCMCcQMuGEYhmFMIGbADcMwDGMCMQNuGIZhGBOIGXDDMAzDmEDMgBuGYRjGBGIG3DAMwzAmEDPghmEYhjGBmAE3DMMwjAnEDLhhGIZhTCBmwI1QRORTIvKnvu83iMhzYvaP/T2nNN0pIs8f5DXyuFbSsf57Fdx3HO5zyDU/LCLnDPOaxngiIj8VkceNOh2GgxlwIxWq+jhV/S8IN1D+34144u7VsO+ziKwXERWR/SJyUETu9RtrETkceD3w9+6+CyLyC4FzfFZEviIiknCtaRH5AxG5RkQeEJE9vs8b8spTnojIh6IaLyLyCBE5ICJH+7a91r2Hm1Oce+LuB/B/gT8ZdSIMh9KoE2AYw0ZESqpaH3U6xoSTgQdU9QgAETkT+JqIfElVdwJvBC5S1UVgUUTOA/4H8H13//8JnAg8S1U16iIiMg18D7gBeJWq7hhcloaDqt4mIv8GnAO8V0SeBnwMeJGq3h137ATfjwuBc0XkaFW9b9SJWemYAp9ARGSzq3geEJEHReRj7vbHish/icjDrqv1Zb5j7hSRd4nItSKyT0S+ICIV3+9PEpErXYX1BaASuOadIvJ8EfkssAX4hqva3uP/PYd0vE9EbnPT8TMReUXKe3KniLzfPWaviPxT4Lx3ish7ReRa4ICIlOLS6XJqzPmS0hl3bKSLPeN9PkZEvuw+B3eIyDt853mviNzjpu9mEXlexK07GbjM9/0n7t8p9+8ZuMba5aPAi1z1+SvA2cBLVfVgxPk93gtco6pvSWOsROR4Efk3V43uE5HvJB0zIv4C+C0ROQn4CvBWVf1piuMm8n6oahW4AnjhKK5vBFBV+0zQBygC1wD/D1iNY2ifCZSBHcAf4FS+vwgsAI92j7sT+ClwDLABuBGnssHd/y4cZVUGfhlYBv7Ud907gecH/w/+3k863N9/xf2tALwaOAAcHXXdwPWvBza75/1RSPqvdn+fSZnOuPMlpTMpLaH3MsN9LuBUpH/kpv944HbgRcCjgbuBY9xjtgKPiLhvnwE+6P6/DvhH4HJA3G0PAKcGjjkf+Jb72ykpn9ubgWMzPOc/BN6G87xXgGeM6H37EHBOwj7fdsv/jzKcdyLvh5uWvwE+Oqrr26f9MQU+eZyGYzjeraoHVLWqqj8EngqsAf5cVZdU9T+BfwNe4zv2b1T1XlV9CPgGjvrCPbYM/JWqLqvql+hUZVnoJx2o6r+6vzVV9QvArW6e0/AxVb3bPe//DlzTu+7d6riD06Qz8nwp0pmUln45FThcVf/ETf/twD8AZwENYBo4UUTKqnqnqt4WcZ6TgXeLyEM4DQLFUdSeO3wdTsPGz0dxGgpvU9Ur/T+IyHNFZGvIdbYA17vejrDPmwL7PwLHWBXdZ/xHSTfEvf7JIvLbafbNAxEp4NzvJo4aD/4+0vsxIBZwngtjxJgBnzw2A3dpdx/uMcDdqtr0bbsL2OT7fr/v/4M4Bsw79h5fpe0d2wv9pAMReb2IXO1VZMBJwMaU1/b3O97lpiXq9zTpjDxfinQmpaVfjgOO8Vf6ON6EI9VxyZ6Dox53i8gFItJ1fXH6YR8LPEZVN6jqI1T1N7Wzb3MvsDZw6BRQw3EZB/kNIGww2z3AE1R1XcTnnwL7vxY4E7hXRP5RRDbE3g0XVb1aVT+RtJ9rePPgL3GM2a04aQ4y0vvRCynuzVrg4UFd30iPGfDJ425gi4gEByDeC2wOvHxbcCqKJO4DNol0jCLeErN/5GClftIhIsfhqMi3A4ep6jocV3Ts6GYf/pG/W9y0+PGnO006Q8+XMp1JaUlD3H2+G7gjUOmvVdUXA6jqv6jqM3EMvRKiDnEaHQfUGawWxbXAowLbnghcH2xEumMIXgr8k4i8PnDMF4D3x1ynA1X9T1V9Hs4AuSfiDKZDRJ4jIheKMwbkahF5uYhcLM70psPc3/+vOHxcRH4gIt8TkcPd374hIl8D3iQi/ywi3xeRfxeR9WnT5svvbwGvAF6Oc3/f7X+HxuR+HOnm/xIR+aKIFHO4N4/F6cYzRowZ8MnjpzgG989FZLWIVETkGTiDjw4A7xGRsjhzhV8KXJDinJcCdeAd4gzueiXxbutdOH2uYfSTjtU4xuYBANeNeFKK4zzeJiLHuurkD3AqySjSpDPqfGnSmSUtUcTd558C8+IMVptxK+aTRORUEXm0iPyiq7CrwCKOmzfIk3BGQcdxEfALgW0n44wnCPJvwFWq+hxV/Uzgtz8Dni4ifyUiR8ZdUEReKSInuMZwLbA+cL0ZVX0l8HfAG1T1RTjldqZvn5cCTVV9lqo+F3jQ3b4eeCWwD9ipqr/gHvu7cWkKSePz3Ty9VFV3AV/C8Uz40zAO92Mv8AJVfTZOI/K59HFv3GfqycC4DipcUZgBnzBUtYHzAj4S+DmwE3i1qi4BL8MZNbwH52V+varelOKcSzgv7htxXvhXE+4e9fgw8Ieu6/ZdIefqNR0/w3FJXopjvB6PMwAsLf+CM6Dodvfzp1E7pkxn6PlSpjN1WmKIu8/ec3AycIebh08Cczj933/ubrsfOAKnERHkZBzPQRyfAV4sIjO+bU8k3IA/EmdwVheqegB4BrAE/FTi+3yfiTPyfQGnAfHn7hgFj2vdv/cG/vcrxcfiGz3v6yr5qfv/I2mP87gMOCEs3WGIyGNwDNvrVPU69/wNnLEB7/XtOg73YwPwryLyfeCXcLqI+rk3LwP+S1V78SgZOeONNDWMiUZE7gR+U1X/Y9RpOdQQkT8DdqvqXyXs93Jga9J+fablOcBLVPVdIvISYLuqfkhEzgKOxRlB/xLgEpwR/O9wjysAz/Yd+8s4o+vfKyKvA7apaleAEhH5EPBwL3kak/uhwEFV/YSI/C3OQMWH6PHeiMhPgDeralLDzxgCFsjFMIxYVDVMvYdxC/CnIrJVVc8ZYJLS8A3gdBH5IY7KfXXg968BrxSRS3C6UsIGoPXLONyP7wKfFZHTcbpSoI97o6pPGXiKjdSYAjcOCUyBG4NARN4KLKjq50adFsMIYgbcMAzDMCYQG8RmGIZhGBOI9YFPEBs3rNOtm/KOB2JkxrxWY8FdN9wy6iQYwB6ae1T18LzPu1lKWo0NhdC6/sWqenre158EzIBPEFs3HcNPvxqcTmoMnbotZDYO/M5jXzDqJBjA37O/16iNsdRQfoXVift9goW0kRoPOcyAG4ZhGGOHAKX4JeYdVrBDzAy4YRiGMZYU0gRRNgNuGIZhGOOFjbKOxwy4YRiGMXYIks6FvoIxA24YhmGMHUJKF/oKxgy4YRiGMZaYCz0eM+CGYRjG+CEg5kKPxQy4YRiGMXY408hGnYrxxgy4YRiGMZaYCz0eM+CGYRjG2OEMYjMJHocZcMMwDGPsMBd6MmbADcMwjLHEXOjxmAE3DMMwxpICJsHjsAaOYRiGMXZ4LvSkT+J5RE4XkZtFZIeIvC/k93eLyNXu53oRaYjIhgFkKXfMgBuGYRhjh4gTiS3pE38OKQIfB84ATgReIyIn+vdR1Y+o6smqejLwfuD7qvrQQDKVM2bADcMwjLGkgCR+EjgN2KGqt6vqEnABcGbM/q8BPp9T8geOGXDDMAxj7MjgQt8oIpf7Pmf7TrMJuNv3fae7rft6IquA04EvDyZH+WOD2AzDMIyxJOViJntUdXvEb2FniFpB/KXAjybFfQ5mwA3DMIwxRNK5yJPYCWz2fT8WuDdi37OYIPc5mAvdMAzDGFP6HcQGXAacICLbRGQKx0hfGNxJROaAXwC+nnceBokpcMMwDGPsEKDYpwBX1bqIvB24GCgC56vqDSLyVvf3c91dXwF8W1UP9HfF4WIG3DAMwxhL8gjkoqoXARcFtp0b+P4p4FN9X2zImAE3DMMwxg5J5yJf0ZgBNwzDMMaS4qgTMOaYATcMwzDGDltONBkz4IZhGMZYYuY7HjPghmEYxlhi85zjMQNuGIZhjCViLvRYzIAbhmEYY4dgLvQkzIAbhmEYY0m/gVwOdcyAG4ZhGGOJmAaPxQy4YRiGMXaYCz0ZM+CGYRjGWGIGPB4z4IZhGMYYIhRtFHosZsANwzCMscNc6MmYATcMwzDGD3EWNDGiMQNuGIZhjCVF0+CxmAE3DMMwxg5zoSdjoWYNwzCMsUQk+ZN8DjldRG4WkR0i8r6IfZ4jIleLyA0i8v288zEoTIEbhmEYY0mhTw0uIkXg48ALgJ3AZSJyoar+zLfPOuDvgNNV9ecickRfFx0ipsANwzCMscNZDzz5k8BpwA5VvV1Vl4ALgDMD+/wa8BVV/TmAqu7OOSsDwwy4YRiGMZZIik8Cm4C7fd93utv8PApYLyL/JSJXiMjr+074kDAXumEYhjGWpHShbxSRy33fz1PV89z/w06gge8l4MnA84AZ4FIR+bGq3pI1vcPGDLhhGIYxlqScB75HVbdH/LYT2Oz7fixwb8g+e1T1AHBARC4BngiMvQE3F/oISTM60jAMYyUiOAYq6ZPAZcAJIrJNRKaAs4ALA/t8HXiWiJREZBXwFODGfHIxWEyBj4g0oyMNwzBWMv3OA1fVuoi8HbgYKALnq+oNIvJW9/dzVfVGEfkWcC3QBD6pqtf3eemhYAZ8dLRGRwKIiDc60gy4YRgGUMghlqqqXgRcFNh2buD7R4CP9H2xIWMu9NGRZnQkInK2iFwuIpc/8NDeoSXOMAxjlKQZgb7SI7WZAR8daUZHoqrnqep2Vd1++Ib1Q0iWYRjGGCCCpPisZMyFPjrSjI40DMNYsRSLK9tAJ2EKfHSkGR1pGIaxIhHyiYV+KGMKfEREjY4ccbIMwzDGAzPQiZgBHyFhoyMNwzAMh0KKYOcrGTPghmEYxliy0gepJWEG3DAMwxg7vD5wIxoz4IZhGMb4IVA0F3osZsANwzCMMcTmeSdhBtwwDMMYOwQQm+gcixlwwzAMY/wQG8SWhBlwwzAMYyyxaWTxmAE3DMMwxhIT4PGYATcMwzDGDiGf5UQPZcyAG4ZhGOOHmAs9CTPghmEYxlhiAjweG6RvGIZhjB15rUYmIqeLyM0iskNE3hfy+3NEZJ+IXO1+/mgA2RkIpsANwzCM8UOEQp/rgYtIEfg48AJgJ3CZiFyoqj8L7PoDVX1JXxcbAabADcMwjLGkIJL4SeA0YIeq3q6qS8AFwJkDT/iQMANuGIZhjB0ZXOgbReRy3+ds32k2AXf7vu90twV5mohcIyLfFJHHDSxTOWMu9BwQkUcBnwCOVNWTROQJwMtU9U9HnDTDMIzJJP0o9D2quj36LF1o4PuVwHGqul9EXgx8DTghdTpHiCnwfPgH4P3AMoCqXgucNdIUGYZhTDgikvhJYCew2ff9WOBe/w6qOq+q+93/LwLKIrIxz3wMCjPg+bBKVX8a2FYfSUoMwzAOEXIYhX4ZcIKIbBORKRxhdWHnNeQocVsCInIajl18MP/c5I+50PNhj4g8Atc1IyK/DNw32iQZhmFMLiL0PQpdVesi8nbgYqAInK+qN4jIW93fzwV+GfhtEakDi8BZqhp0s48lZsDz4W3AecBjROQe4A7gtaNNkmEYxiSTz3rgrlv8osC2c33/fwz4WN8XGgFmwPNBVfX5IrIaKKjqgohsG3WiDMMwJhoLpRqL9YHnw5cBVPWAqi642740wvQYhmFMPnmEYjuEMQXeByLyGOBxwJyIvNL30yxQGU2qDMMwDgEEpGgaMw4z4P3xaOAlwDrgpb7tC8BbRpEgwzCMQwNT2EmYAe8DVf068HUReZqqXjrq9BiGYRwqiIBYH3gsZsDz4SoReRuOO73lOlfV3xhdkgzDMCYbc6HHY3cnHz4LHAW8CPg+TrSfhdgjDMMwjHhsEFssZsDz4ZGq+j+BA6r6aeCXgMePOE2GYRiTi4gzjSzps4IxF3o+LLt/HxaRk4D7ga2jS45hGMbkYy70eMyA58N5IrIe+EOcOLtrgP852iQZhmFMMN56okYkZsDz4buquhe4BDgewCKxGYZh9IeYAI/Fbk8+fDlkm0ViMwzD6BURpFhI/KxkTIH3gUViMwzDGCDmQo/FDHh/WCQ2wzCMQbHCR5knYQa8DywSm2EYxmBwpnmbAY9jZXcg5McrRGRWRMoi8l0R2SMivz7qRBmGYUwuAsVC8mcFs7Jznx8vVNV5HHf6TuBRwLtHmyTDMIzJRkQSPynOcbqI3CwiO0TkfTH7nSoiDRH55VwzMUDMgOdD2f37YuDzqvrQKBNjGIYx8Qh9R2ITkSLwceAM4ETgNSJyYsR+fwFcnH9GBocZ8Hz4hojcBGwHvisihwPVEafJMAxjoslhGtlpwA5VvV1Vl4ALgDND9vtdnOnAu/PNwWAxA54Dqvo+4GnAdlVdBg4S/pAYhmEYaUizkInjQt8oIpf7Pmf7zrIJuNv3fae7zXcZ2QS8Ajh30FnKGxuFnhNuJDbv/wPAgREmxzAMY+JJuR74HlXdHnWKkG0a+P5XwHtVtTFpo97NgBuGYRjjh5DHKPOdwGbf92OBewP7bAcucI33RuDFIlJX1a/1e/FBYwbcMAzDGEtyUMSXASe4a1PcA5wF/Jp/B1VtrVshIp8C/m0SjDeYAc8NEXkZ8Gz36/dV9Rt5X6O58+csvvsdeZ/WyMjtVwYb8MYoOPfAzlEnwQD+fvW6AZ25//W+VbUuIm/HGV1eBM5X1RtE5K3u7xPX7+3HDHgOiMiHcUY7fs7d9A4Rebqqvn+EyTIMw5hcBCj0P85aVS8CLgpsCzXcqvrGvi84RMyA58MvASerahNARD4NXAWYATcMw+iVCRtUNmzMgOfHOsAL4DI3wnQYhmEcAkguCvxQxgx4PvwZcJWIfA/H8fNsTH0bhmH0hynwWMyA94mIFIAm8FTgVBwD/l5VvX+kCTMMw5hkcuoDP5QxA94nqtoUkber6heBC0edHsMwjEMDc6EnYXcnH74jIu8Skc0issH7jDpRhmEYE026UKorFlPg+fAb7t+3+bYpcPwI0mIYhjH5mAs9ETPgOeCP5GMYhmHkgbnQk7C7kwMi8jYRWef7vl5EfmeESTIMw5h8zIUeixnwfHiLqj7sfXFXJnvL6JJjGIYx4Xgu9KTPCsZc6PlQEBFRVQUQkSIwNeI0GYZhTCyCICvcQCdhBjwfLga+KCLn4gxeeyvwrdEmyTAMY8JZ4S7yJMyA58N7gd8CfhvH8fNt4JMjTZFhGMYkY6PQEzEDngPuIiafcD+GYRhG39go9CTMgOeAiJwAfBg4Eah421XV5oEbhmH0irnQY7HmTT78E476rgPPBT4DfHakKTIMw5hkBJtGloAZ8HyYUdXvAqKqd6nqh4BfHHGaDMMwJhiBYjH5k3QWkdNF5GYR2SEi7wv5/UwRuVZErhaRy0XkmQPJzgAwF3o+VN1VyW4VkbcD9wBHjDhNhmEYk02fCtud0vtx4AXATuAyEblQVX/m2+27wIWqqiLyBOCLwGP6uvCQMAWeD+cAq4B3AE8GXge8YZQJMgzDmGjycaGfBuxQ1dtVdQm4ADjTv4Oq7vdieACrcaYCTwSmwHNAVS9z/90PvGmUaTEMwzg0kFQu8gQ2AXf7vu8EntJ1JZFX4AxEPgL4pX4vOizMgOeAiGwHPgAch++equoTRpYowzCMSSedC32jiFzu+36eqp7nnSFk/y6FrapfBb4qIs8G/hfw/KxJHQVmwPPhc8C7geuA5ojTYhiGMfl4LvRk9qjq9ojfdgKbfd+PBe6NOpGqXiIijxCRjaq6J3VaR4QZ8Hx4QFUvHHUiDMMwDh1ycaFfBpwgIttwBhefBfxax1VEHgnc5g5iOwVnHYsH+73wMDADng8fFJFP4oxmrHkbVfUro0uSYRjGhNPnKHRVrbszgy4GisD5qnqDiLzV/f1c4FXA60VkGVgEXu0b1DbWmAHPhzfhTDso03ahK2AG3DAMoxcEkP4nSqnqRcBFgW3n+v7/C+Av+r7QCDADng9PVNXHjzoRhmEYhw65uNAPaWweeD78WEROHHUiDMMwDikslGospsDz4ZnAG0TkDpw+cAHUppEZhmH0iNhqZEmYAc+H00edAMMwjEOOFa6wkzADngOqeteo02AYhnHIUbA+8DjMgBuGYRjjh7nQEzEDbhiGYYwn5kKPxQy4YRiGMX6ITSNLwgy4YRiGMZ7kEMjlUMYMuGEYhjGemAs9FjPghmEYxvhhLvREzIAbhmEY44kp8FjMgBuGYRjjifWBx2IG3DAMwxg/zIWeiBlwwzAMYzwxF3osZsANwzCMMUTMhZ6AGXDDMAxj/BCgYAo8DmveGIZhGONJoZj8SUBETheRm0Vkh4i8L+T314rIte7nv0XkiQPJywAwAz5CROR8EdktItePOi2GYRjjhbuYSdIn7gwiReDjwBnAicBrROTEwG53AL+gqk8A/hdw3gAyMxDMgI+WT2FriRuGYXQjOIPYkj7xnAbsUNXbVXUJuAA407+Dqv63qu51v/4YODbvrAwK6wMfIap6iYhsHXU6DMMwxg9Jux74RhG53Pf9PFX1VPQm4G7fbzuBp8Sc683ANzMlc4SYATcMwzDGk3Trge9R1e0Rv4VJdA3dUeS5OAb8mekSN3rMgI85InI2cDbA5pmpEafGMAxjSHgu9P7YCWz2fT8WuLfrUiJPAD4JnKGqD/Z70WFhfeBjjqqep6rbVXX7xunyqJNjGIYxJCSPUeiXASeIyDYRmQLOAi7suIrIFuArwOtU9ZaBZGVAmAI3DMMwxpM+A7moal1E3g5cDBSB81X1BhF5q/v7ucAfAYcBfyeO4q/HuOTHCjPgI0REPg88B2cQxk7gg6r6j6NNlWEYxhggkksgF1W9CLgosO1c3/+/Cfxm3xcaAWbAR4iqvmbUaTAMwxhb0o1CX7GYATcMwzDGEIuFnoQZcMMwDGMsEVuNLBYz4IZhGMb4IZgCT8AMuGEYhjGGCBStDzwOM+CGYRjGeGIKPBYz4IZhGMb4kU8ktkMaM+CGYRjGGJJ6MZMVixlwwzAMYzwxF3osZsANwzCM8SOnSGyHMmbADcMwjPHEXOixmAE3DMMwxhCLxJaEGXDDMAxjPLFR6LGYATcMwzDGDxEomomKw+6OYRiGMZZYLPR4zIAbhmEY44n1gcdid8cwDMMYP7xIbEmfpNOInC4iN4vIDhF5X8jvjxGRS0WkJiLvGkRWBoUpcMMwDGMM6T8Sm4gUgY8DLwB2ApeJyIWq+jPfbg8B7wBe3tfFRoApcMMwDGM8KRSSP/GcBuxQ1dtVdQm4ADjTv4Oq7lbVy4DlwWRicJgBNwzDMMaPNO5zx4W+UUQu933O9p1lE3C37/tOd9shgbnQDcMwjPEknQt9j6puj/gtrJNce0/QeGEG3DAMwxhT+p5GthPY7Pt+LHBvvycdF8yFbhiGYYwhqV3ocVwGnCAi20RkCjgLuHDgSR8SpsANwzCM8aTPeeCqWheRtwMXA0XgfFW9QUTe6v5+rogcBVwOzAJNETkHOFFV5/u6+BAwA24YhmGMH9488D5R1YuAiwLbzvX9fz+Oa33iMANuGIZhjCcWSTUWM+CGYRjGGGLLiSZhBtwwDMMYT2wxk1jMgBuGYRhjihnwOMyAG4ZhGOOJKfBYzIAbhmEY44dYH3gSZsANwzCM8cQUeCxmwA3DMIwxxQx4HGbADcMwjLFEkpcLXdGYATcMwzDGEMEUeDxmwA3DMIzxxPrAYzEDbhiGYYwfgo1CT8AMuGEYhjGemAKPxQy4YRiGMZ6Y/Y7FDLhhGIYxhlgglyTMgBuGYRjjR07rgR/KWPPGMAzDGFMkxSfhDCKni8jNIrJDRN4X8ruIyN+4v18rIqfkmoUBYgbcMAzDGE9Ekj+xh0sR+DhwBnAi8BoROTGw2xnACe7nbOAT+WdkMJgBNwzDMMYQtw886RPPacAOVb1dVZeAC4AzA/ucCXxGHX4MrBORo/PPT/6YATcMwzDGkz4VOLAJuNv3fae7Les+Y4kNYpsgrnr4wJ41X7v0rlGno082AntGnQjjECiH1etGnYI8mPxygOMGcdIrrrr6Ylm9bmOKXSsicrnv+3mqep77f5iF18D3NPuMJWbAJwhVPXzUaegXEblcVbePOh0rHSuH8cDKIRpVPT2H0+wENvu+Hwvc28M+Y4m50A3DMIxDlcuAE0Rkm4hMAWcBFwb2uRB4vTsa/anAPlW9b9gJ7QVT4IZhGMYhiarWReTtwMVAEThfVW8Qkbe6v58LXAS8GNgBHATeNKr0ZkVUJ8LVbxwiiMjZvv4pY0RYOYwHVg5GP5gBNwzDMIwJxPrADcMwDGMCMQNuDIWkcIbGcBCR80Vkt4hcP+q0rFREZLOIfE9EbhSRG0TknaNOkzGZmAvdGDhuOMNbgBfgTNm4DHiNqv5spAlbgYjIs4H9OJGnThp1elYibpSvo1X1ShFZC1wBvNzeByMrpsCNYZAmnKExBFT1EuChUadjJaOq96nqle7/C8CNTEjkL2O8MANuDIOJDVVoGINERLYCTwJ+MuKkGBOIGXBjGExsqELDGBQisgb4MnCOqs6POj3G5GEG3BgGExuq0DAGgYiUcYz351T1K6NOjzGZmAE3hkGacIaGsSIQEQH+EbhRVT866vQYk4sZcGPgqGod8MIZ3gh8UVVvGG2qViYi8nngUuDRIrJTRN486jStQJ4BvA74RRG52v28eNSJMiYPm0ZmGIZhGBOIKXDDMAzDmEDMgBuGYRjGBGIG3DAMwzAmEDPghmEYhjGBmAE3DMMwjAnEDLhhTBAi8rJhreYmIueIyKq89jMMI19sGplhHKK4AUNEVZs9Hn8nsF1V9+Sxn2EY+WIK3DDGABHZKiI3icgnReR6EfmciDxfRH4kIreKyGnufm8UkY+5/x8pIl8VkWvcz9Pd89woIn8HXAlsFpGPuOe8TkReHXLt1SLy7+45rheRV4vIO4BjgO+JyPfc/T4hIpe7a1j/sbstbL8XisilInKliPyrG/PbMIycMQVuGGOAuyrVDpyVqW7ACT97DfBm4GXAm1T15SLyRhy1+3YR+QJwqar+lbvm+hpgPXA78HRV/bGIvAp4K3A6sNE971NU9T7ftV8FnK6qb3G/z6nqvqCyFpENqvqQe63vAu9Q1Wv9+4nIRuArwBmqekBE3gtMq+qfDOzmGcYKxRS4YYwPd6jqda7L+wbgu+q0sK8Dtobs/4vAJwBUtaGq+9ztd6nqj93/nwl83v19F/B94NTAea4Dni8ifyEiz/KdJ8ivisiVwFXA44ATQ/Z5qrv9RyJyNfAG4LikjBuGkZ3SqBNgGEaLmu//pu97k2zv6gHf/2FLuXagqreIyJOBFwMfFpFvBxWziGwD3gWcqqp7ReRTQCXkdAJ8R1VfkyG9hmH0gClww5hcvgv8NoCIFEVkNmSfS4BXu78fDjwb+Kl/BxE5Bjioqv8M/F/gFPenBWCt+/8sTsNgn4gcCZzhO4V/vx8DzxCRR7rnXiUij+ovm4ZhhGEK3DAml3cC57krijVwjPl9gX2+CjwNpz9dgfeo6v2BfR4PfEREmsCyex6A84Bvish9qvpcEbkKx7V/O/Aj3/HB/d4IfF5Ept3f/xC4pf/sGobhxwaxGYZhGMYEYi50wzAMw5hAzIAbhmEYxgRiBtwwDMMwJhAz4IZhGIYxgZgBNwzDMIwJxAy4YRiGYUwgZsANwzAMYwIxA24YhmEYE8j/B4tsqV+YbomIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x504 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "f = plt.figure(constrained_layout=False, figsize=(7, 7))\n",
    "gs = f.add_gridspec(nrows=4, ncols=1)\n",
    "ax1 = f.add_subplot(gs[0, 0])\n",
    "ax2 = f.add_subplot(gs[1, 0])\n",
    "ax3 = f.add_subplot(gs[2:, 0])\n",
    "\n",
    "ax1.plot(hidden_trajectory[:200])\n",
    "ax1.set_ylabel(\"macrostate\")\n",
    "ax2.plot(observation_trajectory[:200])\n",
    "ax2.set_xlabel(\"time\")\n",
    "ax2.set_ylabel(\"microstate\")\n",
    "\n",
    "probabilities = np.zeros((2, 3))\n",
    "for i_hidden in range(2):\n",
    "    for i_obs in range(3):\n",
    "        # indices where macro state is i_hidden\n",
    "        hidden_ix = np.argwhere(hidden_trajectory == i_hidden).flatten()\n",
    "        subselected_obs = observation_trajectory[hidden_ix]\n",
    "        # indices where micro state is i_obs given i_hidden\n",
    "        obs_ix = np.argwhere(subselected_obs == i_obs).flatten()\n",
    "        probabilities[i_hidden, i_obs] = len(obs_ix) / len(hidden_ix)\n",
    "\n",
    "cb = ax3.matshow(probabilities, cmap=\"Reds\")\n",
    "ax3.xaxis.set_ticks_position('bottom')\n",
    "ax3.set_xlabel(\"micro state\")\n",
    "ax3.set_ylabel(\"macro state\")\n",
    "ax3.set_title(r\"conditional probabilities $P(Y_t\\in s_\\mathrm{micro}\\mid X_t\\in s_\\mathrm{macro})$\")\n",
    "f.colorbar(cb, ax=ax3)\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the output probabilities in this case map to discrete values, this is not the only option. We implement\n",
    "\n",
    "- the [discrete output model](../api/generated/deeptime.markov.hmm.DiscreteOutputModel.rst#deeptime.markov.hmm.DiscreteOutputModel), which maps to a discrete observable space and is parameterized by a row-stochastic output probability matrix $E\\in\\mathbb{R}^{n\\times m}$, where $n$ is the number hidden and $m$ the number of observable states,\n",
    "- and the [Gaussian output model](../api/generated/deeptime.markov.hmm.GaussianOutputModel.rst#deeptime.markov.hmm.GaussianOutputModel), which maps to a one-dimensional output space with $m$ Gaussians, each parameterized by mean and standard deviation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation\n",
    "\n",
    "The estimation of an HMM given sequences of observations can be performed using the Baum-Welch algorithm <cite data-footcite=\"baum1967inequality\">(Baum, 1967)</cite>, which belongs to the family of expectation-maximization algorithms.\n",
    "\n",
    "It crucially requires an initial guess for emission probabilities and hidden state transition matrix, as the optimization is prone to getting stuck in local optima. After optimization, the likelihood for a trained model is available and can be compared to models with different initial guesses.\n",
    "\n",
    "Going back to above example, we can try to fit an HMM based on the observation trajectory for different initial guesses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:12.849633Z",
     "iopub.status.busy": "2020-12-02T13:38:12.844601Z",
     "iopub.status.idle": "2020-12-02T13:38:12.866780Z",
     "shell.execute_reply": "2020-12-02T13:38:12.866257Z"
    }
   },
   "outputs": [],
   "source": [
    "hmm_est_real = hmm.MaximumLikelihoodHMM(\n",
    "    lagtime=1, initial_model=ground_truth\n",
    ").fit(observation_trajectory).fetch_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:12.873087Z",
     "iopub.status.busy": "2020-12-02T13:38:12.872586Z",
     "iopub.status.idle": "2020-12-02T13:38:12.874814Z",
     "shell.execute_reply": "2020-12-02T13:38:12.874478Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute deviation: 0.034\n"
     ]
    }
   ],
   "source": [
    "print(\"Maximum absolute deviation: {:.3f}\".format(\n",
    "    np.max(np.abs(hmm_est_real.transition_model.transition_matrix - P))\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When initialized with the ground truth, the deviation is not large from the initial model. Deeptime furthermore offers some initialization heuristics should the ground truth be unknown.\n",
    "\n",
    "For discrete output models:\n",
    "\n",
    "1. [init.discrete.metastable_from_data](../api/generated/deeptime.markov.hmm.init.discrete.metastable_from_data.rst#deeptime.markov.hmm.init.discrete.metastable_from_data) or [init.discrete.metastable_from_msm](../api/generated/deeptime.markov.hmm.init.discrete.metastable_from_msm.rst#deeptime.markov.hmm.init.discrete.metastable_from_msm): A reversible [MSM](mlmsm.ipynb) is estimated from given discrete trajectories. This is optional in case the MSM is already available. Then, the estimated transition matrix $P\\in\\mathbb{R}^{n\\times n}$ is coarse-grained with PCCA+ <cite data-footcite=\"roblitz2013fuzzy\">(Röblitz, 2013)</cite>, yielding $m$ metastable/macro states and a membership matrix $M\\in\\mathbb{R}^{n\\times m}$. The transition matrix is projected into hidden space via\n",
    "    $$ P_\\mathrm{coarse} = (M^\\top M)^{-1}M^\\top PM. $$\n",
    "    This procedure is described in greater detail in <cite data-footcite=\"noe2013projected\">(Noé, 2013)</cite>.\n",
    "2. [init.discrete.random_guess](../api/generated/deeptime.markov.hmm.init.discrete.random_guess.rst#deeptime.markov.hmm.init.discrete.random_guess): Initializes hidden transition matrix and initial distribution uniformly, the emission probabilities are drawn from the uniform distribution $\\mathcal{U}(0, 1)$ and normalized to form a row-stochastic matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a random initial model and use it as starting point for the maximum-likelihood estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:12.880890Z",
     "iopub.status.busy": "2020-12-02T13:38:12.880429Z",
     "iopub.status.idle": "2020-12-02T13:38:12.978167Z",
     "shell.execute_reply": "2020-12-02T13:38:12.977643Z"
    }
   },
   "outputs": [],
   "source": [
    "init_hmm_random = hmm.init.discrete.random_guess(n_observation_states=3, n_hidden_states=2)\n",
    "hmm_est_random = hmm.MaximumLikelihoodHMM(\n",
    "    lagtime=1, initial_model=init_hmm_random\n",
    ").fit(observation_trajectory).fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also here, we can compute the maximum deviation from ground truth hidden transition matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:12.981523Z",
     "iopub.status.busy": "2020-12-02T13:38:12.981118Z",
     "iopub.status.idle": "2020-12-02T13:38:12.983565Z",
     "shell.execute_reply": "2020-12-02T13:38:12.983158Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute deviation: 0.258\n"
     ]
    }
   ],
   "source": [
    "print(\"Maximum absolute deviation: {:.3f}\".format(\n",
    "    np.max(np.abs(hmm_est_random.transition_model.transition_matrix - P))\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And evaluate the likelihood of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1652.9913478370843"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hmm_est_random.likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as compute the likelihood of the generated data given our original model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1660.3598109385"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ground_truth.compute_observation_likelihood(observation_trajectory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian sampling\n",
    "\n",
    "Bayesian sampling of HMMs can be used to obtain confidences. In particular, Gibbs sampling can be performed around a reference HMM, as introduced in <cite data-footcite=\"chodera2011bayesian\">(Chodera, 2011)</cite>. In deeptime, a [BayesianHMM](../api/generated/deeptime.markov.hmm.BayesianHMM.rst#deeptime.markov.hmm.BayesianHMM) estimator is provided.\n",
    "\n",
    "The easiest method to obtain a BayesianHMM estimator from data is the [default()](../api/generated/deeptime.markov.hmm.BayesianHMM.rst#deeptime.markov.hmm.BayesianHMM.default) static method. It estimates an ML-MSM on data, coarse-grains it to the requested number of hidden states and uses it as starting point for an ML-HMM estimation. The resulting model is used as start sample for the BayesianHMM estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:12.986923Z",
     "iopub.status.busy": "2020-12-02T13:38:12.986472Z",
     "iopub.status.idle": "2020-12-02T13:38:13.071283Z",
     "shell.execute_reply": "2020-12-02T13:38:13.070701Z"
    }
   },
   "outputs": [],
   "source": [
    "estimator = hmm.BayesianHMM.default(observation_trajectory, n_hidden_states=2, lagtime=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum-likelihood estimation is done for fewer iterations and with less precision than the default settings, as a rough estimate suffices to start sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:13.075053Z",
     "iopub.status.busy": "2020-12-02T13:38:13.074611Z",
     "iopub.status.idle": "2020-12-02T13:38:13.076807Z",
     "shell.execute_reply": "2020-12-02T13:38:13.077199Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum absolute deviation: 0.282\n"
     ]
    }
   ],
   "source": [
    "print(\"Maximum absolute deviation: {:.3f}\".format(\n",
    "    np.max(np.abs(estimator.initial_hmm.transition_model.transition_matrix - P))\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling `fit` on the estimator starts sampling and yields a [BayesianHMMPosterior](../api/generated/deeptime.markov.hmm.BayesianHMMPosterior.rst#deeptime.markov.hmm.BayesianHMMPosterior) instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:13.081759Z",
     "iopub.status.busy": "2020-12-02T13:38:13.081353Z",
     "iopub.status.idle": "2020-12-02T13:38:17.894723Z",
     "shell.execute_reply": "2020-12-02T13:38:17.894297Z"
    }
   },
   "outputs": [],
   "source": [
    "posterior = estimator.fit(\n",
    "    observation_trajectory, \n",
    "    n_burn_in=100,  # number of sampling steps to discard at the beginning \n",
    "    n_thin=10  # take every 10th sample\n",
    ").fetch_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any quantity that can be evaluated on the hidden Markov state model object can also be accessed through the posterior but equipped with uncertainties. For this, call `gather_stats` with the name of the method or quantity that should be evaluated. If the method requires arguments, they can also be provided. If it is a method of a lower-level object, it can be separated by slashes `/`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:17.898119Z",
     "iopub.status.busy": "2020-12-02T13:38:17.897724Z",
     "iopub.status.idle": "2020-12-02T13:38:17.902143Z",
     "shell.execute_reply": "2020-12-02T13:38:17.901708Z"
    }
   },
   "outputs": [],
   "source": [
    "stats = posterior.gather_stats(\"transition_model/transition_matrix\", confidence=0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check on the confidence of, e.g., the transition matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:17.926502Z",
     "iopub.status.busy": "2020-12-02T13:38:17.924522Z",
     "iopub.status.idle": "2020-12-02T13:38:18.070038Z",
     "shell.execute_reply": "2020-12-02T13:38:18.070532Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATcklEQVR4nO3db2yVZ5rf8e8FIZn1zHam2tAyjcNx0KQFtgJKXZeq/6hohzBCgkqMlKy1I42isaicaiuBRKZIu1SLhfqiVUWW2WAtEV3paKOVd0rTIUOmSE3nxZZtnIXxmWw2Iy87BpqRYLKULXVFFPbqCzsjY47tA5zj43P7+5GOdJ77uX2e686t/PLk+RuZiSSp861odwGSpOYw0CWpEAa6JBXCQJekQhjoklQIA12SCrFgoEfEaxFxPSJ+OMf6iIjjETEeEWMRsbX5ZUqSFtLIHvpp4Ll51u8Cnp3+DAC/+ehlSZIe1GMLdcjM70dEzzxd9gC/nVN3KF2IiC9ExBcz8yfz/e6TTz6ZPT3z/awkabZ33333p5m5ut66BQO9AU8BV2csX5tuuy/QI2KAqb141q5dy+joaBM2L0nLR0RMzLWuGSdFo05b3ecJZOZwZvZmZu/q1XX/AyNJekjNCPRrwNMzlruBD5vwu5KkB9CMQH8D+Nr01S7bgFsLHT+XJDXfgsfQI+J3gO3AkxFxDfg1YBVAZr4KvAl8BRgHJoGvt6pYSdLcGrnK5YUF1icw2LSKJEkPxTtFJanFqtUqPT09rFixgp6eHqrVaku204zLFiVJc6hWqwwMDDA5OQnAxMQEAwMDAPT39zd1W9GuNxb19vam16FLKt1nPvMZ7ty5c197pVLhxz/+8QP/XkS8m5m99dZ5yEWSWqhemANcuXKl6dsy0CWphSqVSt32tWvXNn1bBroktdDQ0BBdXV33tHV1dTE0NNT0bRnoktRC/f39DA8PU6lUiAgqlQrDw8NNPyEKnhSVpI7iSVFJWgYMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJapKTJ0+yZs0atmzZwrp16zh9+vSibt9Al6QmGRsb48iRI1y6dImRkREOHDiwqNs30CWpSWq1Ghs2bACgu7ubu3fvLur2DXRJapJarcb69evJTI4fP87u3bsXdfu+U1SSmuDq1avcvn2bnTt3smrVKvr6+jhx4gRnzpzh7NmzXL9+ncHBQSYnJ+9Z/vKXv9y0Ggx0SWqCsbExduzYwblz5+5p37t3L3v37uXmzZscPHiQU6dO3bPczED3kIskNUGtVmPz5s1zrj969CiDg4NzLjeDgS5JTVCr1di0adN97ZnJoUOH2LVrF1u3br1vuZk85CJJTVCtVuu2v/LKK5w/f55bt24xPj7Oxx9/fM/y/v37m1aDr6CTpA7iK+gkaRkw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhGgr0iHguIj6IiPGIeLnO+s9HxH+JiB9ExHsR8fXmlypJms+CgR4RK4ETwC5gI/BCRGyc1W0Q+KPM3AxsB/5dRDze5FolSfNoZA+9DxjPzMuZ+THwOrBnVp8Efj4iAvgc8GfAJ02tVJI0r0YC/Sng6ozla9NtM/0GsAH4EKgBv5KZfzH7hyJiICJGI2L0xo0bD1myJKmeRgI96rTNfgDMTuAS8NeALcBvRMRfuu+PMoczszcze1evXv2ApUqS5tNIoF8Dnp6x3M3UnvhMXwe+nVPGgT8F1jenRElSIxoJ9HeAZyPimekTnc8Db8zqcwXYARARfxX4G8DlZhYqSZrfgs9Dz8xPIuIl4C1gJfBaZr4XEfun178K/DpwOiJqTB2iOZSZP21h3ZKkWRp6wUVmvgm8Oavt1RnfPwSa92I8SdID805RSSrEsg30kydPsmbNGrZs2cK6des4ffp0u0uSpEeybAN9bGyMI0eOcOnSJUZGRjhw4EC7S5KkR7JsA71Wq7FhwwYAuru7uXv3bpsrkqRHs6wDff369WQmx48fZ/fu3e0uSZIeybIM9KtXr3L79m127txJX18fN2/e5MSJE1y+fJkXX3yRffv2Ady3LElLWUOXLZZmbGyMHTt2cO7cuXvaP//5z3Pq1KmfBfi6devuWZakpWxZ7qHXajU2b97c7jIkqamWbaBv2rSp3WVIUlMty0CvVqv09/ff1/7RRx+xf/9+Ll68yLFjx+5blqSlLDJnPwl3cfT29ubo6Ghbti1JnSoi3s3M3nrrluUeuiSVyECXpEIY6JJUCANdkgphoEtSIQx0SSpEkYFerVbp6elhxYoV9PT0UK1W212SJLVccc9yqVarDAwMMDk5CcDExAQDAwMAdW8mkqRSdOSNRdu3b59z3YULF7hz58597U888QTbtm2r+zdvv/32Q9UhSYttWd1YVC/M52uXpFJ05CGX+faoe3p6mJiYuK+9Uqm4Jy6paMXtoQ8NDdHV1XVPW1dXF0NDQ22qSJIWR3GB3t/fz/DwMJVKhYigUqkwPDzsCVFJxevIk6KStFwtq5OikrRcGeiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQjQU6BHxXER8EBHjEfHyHH22R8SliHgvIv57c8uUJC1kwactRsRK4ATwz4BrwDsR8UZm/tGMPl8AvgU8l5lXIuKvtKheSdIcGtlD7wPGM/NyZn4MvA7smdXnl4BvZ+YVgMy83twyJUkLaSTQnwKuzli+Nt02018H/nJEvB0R70bE15pVoCSpMY284CLqtM1+RONjwN8GdgA/B/yPiLiQmT+654ciBoABgLVr1z54tZKkOTWyh34NeHrGcjfwYZ0+5zLz/2bmT4HvA5tn/1BmDmdmb2b2rl69+mFrliTV0UigvwM8GxHPRMTjwPPAG7P6/GfgH0bEYxHRBfxd4P3mlipJms+Ch1wy85OIeAl4C1gJvJaZ70XE/un1r2bm+xFxDhgD/gL4rcz8YSsLlyTdyzcWSVIH8Y1FkrQMGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklSIhgI9Ip6LiA8iYjwiXp6n39+JiLsRsa95JUqSGrFgoEfESuAEsAvYCLwQERvn6PdvgbeaXaQkaWGN7KH3AeOZeTkzPwZeB/bU6fcvgd8DrjexPklSgxoJ9KeAqzOWr023/UxEPAX8c+DV5pWm5eDkyZOsWbOGLVu2sG7dOk6fPt3ukqSO1UigR522nLX8H4BDmXl33h+KGIiI0YgYvXHjRoMlqmRjY2McOXKES5cuMTIywoEDB9pdktSxGgn0a8DTM5a7gQ9n9ekFXo+IHwP7gG9FxN7ZP5SZw5nZm5m9q1evfriKVZRarcaGDRsA6O7u5u7defcJJM3jsQb6vAM8GxHPAP8LeB74pZkdMvOZT79HxGngO5l5pnllqlS1Wo3169eTmRw/fpzdu3e3uySpYy24h56ZnwAvMXX1yvvA72bmexGxPyL2t7pAlevq1avcvn2bnTt30tfXx82bNzlx4gSXL1/mxRdfZN++qatfz5w5wze+8Q327NnD9773vTZXLS1dkTn7cPji6O3tzdHR0bZsW0vD2bNneeWVVzh37lzd9fv27WNkZORnyzdv3uTgwYOcOnVqsUqUlpyIeDcze+ut805RtU2tVmPz5s0N9z969CiDg4MtrEjqbAa62qZWq7Fp06YF+2Umhw4dYteuXWzdunURKpM6UyMnRaWWqFardds/+ugjDh8+zMWLFzl27Bif/exnOX/+PLdu3WJ8fJz9+z11I9XjMXRJ6iAeQ5ekZcBAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIGulqlWq/T09LBixQp6enrmvO5cUnN4Y5FaolqtMjAwwOTkJAATExMMDAwA0N/f387SpGK5h66WOHz48M/C/FOTk5McPny4TRVJ5TPQ1RJXrlx5oHZJj85AV0usXbv2gdolPToDXS0xNDREV1fXPW1dXV0MDQ21qSKpfAa6WqK/v5/h4WEqlQoRQaVSYXh42BOiUgv5tEVJ6iA+bVGSlgEDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEI0FOgR8VxEfBAR4xHxcp31/RExNv35/YjY3PxSJUnzWTDQI2IlcALYBWwEXoiIjbO6/SnwjzNzE/DrwHCzC5Ukza+RPfQ+YDwzL2fmx8DrwJ6ZHTLz9zPz5vTiBaC7uWVKkhbSSKA/BVydsXxtum0uLwLffZSiJEkP7rEG+kSdtrpvlo6If8JUoP+DOdYPAAMAa9eubbBESVIjGtlDvwY8PWO5G/hwdqeI2AT8FrAnMz+q90OZOZyZvZnZu3r16oepV5I0h0YC/R3g2Yh4JiIeB54H3pjZISLWAt8Gfjkzf9T8MiVJC1nwkEtmfhIRLwFvASuB1zLzvYjYP73+VeBXgV8AvhURAJ9kZm/rypYkzRaZdQ+Ht1xvb2+Ojo62ZduS1Kki4t25dpi9U1SSCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS4tcydPnmTNmjVs2bKFdevWcfr06XaXpIdkoEvL3NjYGEeOHOHSpUuMjIxw4MCBdpekh2SgS8tcrVZjw4YNAHR3d3P37t02V6SHZaBLy1ytVmP9+vVkJsePH2f37t3tLkkPqZF3ikoq1NWrV7l9+zY7d+5k1apV9PX1ceLECc6cOcPZs2e5fv06g4ODfOlLX2JoaIhbt24xMjLS7rI1BwNdWsbGxsbYsWMH586du6d979697N27l5s3b3Lw4EFOnTrFqVOn2LdvX5sqVSM85CItY7Vajc2bN8+5/ujRowwODi5iRXoUBrq0jNVqNTZt2nRfe2Zy6NAhdu3axdatW9tQmR6Gh1ykZaxardZtf+WVVzh//jy3bt1ifHycr371qxw+fJiLFy9y7NgxvvnNby5ypWqEL4mWpA7iS6IlaRkw0CWpEAa6JBXCQJekQhjoklQIA10qWLVapaenhxUrVtDT0zPnZYoqg9ehS4WqVqsMDAwwOTkJwMTEBAMDAwD09/e3szS1iNehSx1s+/btc667cOECd+7cua/9iSeeYNu2bXX/5u23325SZWoVr0OXlqF6YT5fuzqfh1ykDjbfHnVPTw8TExP3tVcqFffEC+UeulSooaEhurq67mnr6upiaGioTRWp1Qx0qVD9/f0MDw9TqVSICCqVCsPDw54QLZgnRSWpg3hSVJKWgYYCPSKei4gPImI8Il6usz4i4vj0+rGI8In4krTIFgz0iFgJnAB2ARuBFyJi46xuu4Bnpz8DwG82uU5J0gIa2UPvA8Yz83Jmfgy8DuyZ1WcP8Ns55QLwhYj4YpNrlSTNo5FAfwq4OmP52nTbg/YhIgYiYjQiRm/cuPGgtUqS5tHIjUVRp232pTGN9CEzh4FhgIi4ERH33/XQXE8CP23xNlrNMSwNnT6GTq8fHMOnKnOtaCTQrwFPz1juBj58iD73yMzVDWz7kUTE6FyX93QKx7A0dPoYOr1+cAyNaOSQyzvAsxHxTEQ8DjwPvDGrzxvA16avdtkG3MrMnzS5VknSPBbcQ8/MTyLiJeAtYCXwWma+FxH7p9e/CrwJfAUYByaBr7euZElSPQ09nCsz32QqtGe2vTrjewKDzS2tKYbbXUATOIalodPH0On1g2NYUNtu/ZckNZe3/ktSIQx0SSpEEYFewrNmGhjD9oi4FRGXpj+/2o465xIRr0XE9Yj44RzrO2EOFhrDUp+DpyPiv0XE+xHxXkT8Sp0+S3oeGhzDUp+Hz0TE/4yIH0yP4d/U6dOaecjMjv4wdeXNnwDrgMeBHwAbZ/X5CvBdpm6A2gb8QbvrfogxbAe+0+5a5xnDPwK2Aj+cY/2SnoMGx7DU5+CLwNbp7z8P/KgD/11oZAxLfR4C+Nz091XAHwDbFmMeSthDL+FZM42MYUnLzO8DfzZPl6U+B42MYUnLzJ9k5h9Of/8/wPvc/wiOJT0PDY5hSZv+Z3t7enHV9Gf21SctmYcSAr1pz5ppo0br+3vT/xv33Yj4xcUprWmW+hw0qiPmICJ6gL/F1N7hTB0zD/OMAZb4PETEyoi4BFwH/mtmLso8lPCS6KY9a6aNGqnvD4FKZt6OiK8AZ5h6XHGnWOpz0IiOmIOI+Bzwe8C/ysw/n726zp8suXlYYAxLfh4y8y6wJSK+APyniPibmTnz3ExL5qGEPfSWPGtmkS1YX2b++af/G5dTN3qtiognF6/ER7bU52BBnTAHEbGKqSCsZua363RZ8vOw0Bg6YR4+lZn/G3gbeG7WqpbMQwmBXsKzZhYcQ0SsiYiY/t7H1Nx9tOiVPrylPgcLWupzMF3bKeD9zPz3c3Rb0vPQyBg6YB5WT++ZExE/B/xT4I9ndWvJPHT8IZcs4FkzDY5hH/AvIuIT4P8Bz+f06fKlICJ+h6mrD56MiGvArzF1Mqgj5gAaGsOSngPg7wO/DNSmj98C/GtgLXTMPDQyhqU+D18E/mNMve1tBfC7mfmdxcgkb/2XpEKUcMhFkoSBLknFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgrx/wFatWs7sNIgFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [0, 1, 2, 3]\n",
    "y = stats.mean.flatten()\n",
    "errors = np.concatenate((stats.mean - stats.L, -stats.mean + stats.R)).T\n",
    "f, ax = plt.subplots()\n",
    "ax.errorbar(x, y, xerr=errors, fmt='o', color=\"k\")\n",
    "ax.annotate(r\"$P_{11}$\", (x[0]+.05, y[0]+.05))\n",
    "ax.annotate(r\"$P_{12}$\", (x[1]+.05, y[1]+.05))\n",
    "ax.annotate(r\"$P_{21}$\", (x[2]+.05, y[2]+.05))\n",
    "ax.annotate(r\"$P_{22}$\", (x[3]-.1, y[3]-.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, it is possible to instantiate a Bayesian HMM estimator directly by using a previously estimated HMM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:18.073776Z",
     "iopub.status.busy": "2020-12-02T13:38:18.073263Z",
     "iopub.status.idle": "2020-12-02T13:38:18.075233Z",
     "shell.execute_reply": "2020-12-02T13:38:18.074807Z"
    }
   },
   "outputs": [],
   "source": [
    "estimator = hmm.BayesianHMM(initial_hmm=posterior.prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And perform sampling as above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-02T13:38:18.080310Z",
     "iopub.status.busy": "2020-12-02T13:38:18.079690Z",
     "iopub.status.idle": "2020-12-02T13:38:21.810058Z",
     "shell.execute_reply": "2020-12-02T13:38:21.810376Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASdElEQVR4nO3df2zc913H8dfbaVrwNjJEDR5xc1drHXZAsQnGBPHLKFAnJVKMlErdrE1M0U5GLhpSIiXDYgQtp4g/QMidp+Y0Vx7SiQp5JZS1c0skwv6YMuos3t2y0slkOA6dlDQzHp5Ronhv/vC1ss9n+5vke767zz0fkqX7fj4f+/v+5hO/8sn3+737mrsLAFD7GipdAAAgHgQ6AASCQAeAQBDoABAIAh0AAvFQpXb86KOPejKZrNTuAaAmXbp06R13byrVV7FATyaTmpycrNTuAaAmmdnMen2ccgGAQBDoABAIAh0AAkGgA0AgCHQACASBDgBlls1mlUwm1dDQoGQyqWw2W5b9VOy2RQCoB9lsVqlUSouLi5KkmZkZpVIpSVJ/f3+s+2KFDgBl0tPTo6NHj74X5u9aXFzU0NBQ7Psj0AGgjG7fvl2y/dq1a7Hvi0AHgDK5cOGCEolEyb5du3bFvj8CHQDKKJ1Oq7GxcVVbY2Oj0ul07Psi0AGgjPr7+5XJZJRIJGRmSiQSymQysV8QlSSr1DNFu7q6nA/nAoB7Y2aX3L2rVB8rdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABICZnz55Vc3OzOjs71draqrGxsS3dP4EOADHJ5XI6deqUpqamND4+rmPHjm3p/gl0AIhJPp9Xe3u7JKmlpUVLS0tbun8CHQBiks/n1dbWJnfX8PCwDh06tKX75yHRABCD2dlZLSwsqLe3V9u3b1d3d7dGRkZ07tw5vfLKK7px44YGBwe1uLi4avvJJ5+MrQYCHQBikMvltH//fk1MTKxq7+vrU19fn+bm5nT8+HGNjo6u2o4z0DnlAgAxyOfz6ujoWLf/9OnTGhwcXHc7DgQ6AMQgn89rz549a9rdXSdOnNDBgwe1d+/eNdtx4pQLAMQgm82WbH/uued0/vx5zc/Pa3p6Wnfu3Fm1PTAwEFsNPFMUAGoIzxQFgDpAoANAIAh0AAgEgQ4AgSDQASAQkQLdzA6Y2VtmNm1mJ0v07zCzfzazb5nZFTP7ZPylAgA2smmgm9k2SSOSDkraLemjZra7aNigpO+4e4ekHkl/bWYPx1wrAGADUVbo3ZKm3f2qu9+R9KKkw0VjXNIHzMwkvV/SDyTdjbVSAMCGogT6TkmzK7avF9pW+rykdklvS8pL+rS7/ziWCgEAkUQJdCvRVvz20l5JU5J+XlKnpM+b2U+t+UFmKTObNLPJmzdv3mOpAICNRAn065IeW7HdouWV+EqflPSSL5uW9D1JbcU/yN0z7t7l7l1NTU33WzMAoIQogf6GpCfM7PHChc5nJL1cNOaapP2SZGY/J+kXJF2Ns1AAwMY2/bRFd79rZs9Kek3SNkkvuPsVMxso9D8v6XOSxswsr+VTNCfc/Z0y1g0AKBLp43Pd/VVJrxa1Pb/i9duS4nvsBgDgnvFOUQAIBIEOAIEg0AEgEAQ6AASibgP97Nmzam5uVmdnp1pbWzU2NlbpkgDggdRtoOdyOZ06dUpTU1MaHx/XsWPHKl0SADyQug30fD6v9vZ2SVJLS4uWlpYqXBEAPJi6DvS2tja5u4aHh3Xo0KFKlwQAD6QuA312dlYLCwvq7e1Vd3e35ubmNDIyoqtXr+ro0aM6cuSIJK3ZBoBqFumdoqHJ5XLav3+/JiYmVrXv2LFDo6Oj7wV4a2vrqm0AqGZ1uULP5/Pq6OiodBkAEKu6DfQ9e/ZUugwAiFVdBno2m1V/f/+a9lu3bmlgYECXL1/WmTNn1mwDQDUz9+KHD22Nrq4un5ycrMi+AaBWmdkld+8q1VeXK3QACBGBDgCBINABIBAEOgAEgkAHgEAEGejZbFbJZFINDQ1KJpPKZrOVLgkAyi64t/5ns1mlUiktLi5KkmZmZpRKpSSp5L3nABCKmrwPvaenZ92+ixcv6vbt22vaH3nkEe3bt6/k91y4cOG+6gCArVZX96GXCvON2gEgFDV5ymWjFXUymdTMzMya9kQiwUocQNCCW6Gn02k1NjauamtsbFQ6na5QRQCwNYIL9P7+fmUyGSUSCZmZEomEMpkMF0QBBK8mL4oCQL2qq4uiAFCvCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAIRKdDN7ICZvWVm02Z2cp0xPWY2ZWZXzOzf4i0TALCZTT+cy8y2SRqR9PuSrkt6w8xedvfvrBjzQUlfkHTA3a+Z2c+WqV4AwDqirNC7JU27+1V3vyPpRUmHi8Z8TNJL7n5Nktz9RrxlAgA2EyXQd0qaXbF9vdC20kck/bSZXTCzS2b2iVI/yMxSZjZpZpM3b968v4oBACVFCXQr0Vb8iV4PSfoVSX8gqVfSn5vZR9Z8k3vG3bvcvaupqemeiwUArC/KAy6uS3psxXaLpLdLjHnH3X8k6Udm9jVJHZK+G0uVAIBNRVmhvyHpCTN73MwelvSMpJeLxvyTpN8ys4fMrFHSr0l6M95SAQAb2XSF7u53zexZSa9J2ibpBXe/YmYDhf7n3f1NM5uQlJP0Y0lfdPdvl7NwAMBqPOACAGoID7gAgDpAoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIRKRAN7MDZvaWmU2b2ckNxv2qmS2Z2ZH4SgQARLFpoJvZNkkjkg5K2i3po2a2e51xfyXptbiLRJjOnj2r5uZmdXZ2qrW1VWNjY5UuCahpUVbo3ZKm3f2qu9+R9KKkwyXG/YmkL0u6EWN9CFgul9OpU6c0NTWl8fFxHTt2rNIlATUtSqDvlDS7Yvt6oe09ZrZT0h9Ken6jH2RmKTObNLPJmzdv3mutCEw+n1d7e7skqaWlRUtLSxWuCKhtUQLdSrR50fbfSjrh7hv+Rrp7xt273L2rqakpYokIVT6fV1tbm9xdw8PDOnToUKVLAmpalEC/LumxFdstkt4uGtMl6UUz+y9JRyR9wcz64igQYZqdndXCwoJ6e3vV3d2tubk5jYyM6OrVqzp69KiOHFm+rn7u3Dl96lOf0uHDh/X6669XuGqguj0UYcwbkp4ws8cl/bekZyR9bOUAd3/83ddmNibpK+5+Lr4yEZpcLqf9+/drYmJiVfuOHTs0Ojr6XqD39fWpr69Pc3NzOn78uJ588slKlAvUhE0D3d3vmtmzWr57ZZukF9z9ipkNFPo3PG8OlJLP59XR0RF5/OnTpzU4OFjGioDaF2WFLnd/VdKrRW0lg9zd/+jBy0Lo8vm8nnrqqU3HubtOnjypgwcPau/evVtQGVC7IgU6ELdsNluy/datWxoaGtLly5d15swZve9979P58+c1Pz+v6elpDQwMbHGlQO0w9+IbVrZGV1eXT05OVmTfAFCrzOySu3eV6uOzXAAgEAQ6AASCQAeAQBDoABAIAh0AAkGgo2yy2aySyaQaGhqUTCbXvVURQDy4Dx1lkc1mlUqltLi4KEmamZlRKpWSJPX391eyNCBYrNBRFkNDQ++F+bsWFxc1NDRUoYqA8BHoKItr167dUzuAB0egoyx27dp1T+0AHhyBjrJIp9NqbGxc1dbY2Kh0Ol2hioDwEegoi/7+fmUyGSUSCZmZEomEMpkMF0SBMuLDuQCghvDhXABQBwh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBiBToZnbAzN4ys2kzO1miv9/McoWvr5tZR/ylAgA2smmgm9k2SSOSDkraLemjZra7aNj3JP2Ou++R9DlJmbgLBQBsLMoKvVvStLtfdfc7kl6UdHjlAHf/urvPFTYvSmqJt0wAwGaiBPpOSbMrtq8X2tZzVNJXH6QoAMC9eyjCGCvR5iUHmv2ulgP9N9fpT0lKSdKuXbsilggAiCLKCv26pMdWbLdIert4kJntkfRFSYfd/VapH+TuGXfvcveupqam+6kXALCOKIH+hqQnzOxxM3tY0jOSXl45wMx2SXpJ0sfd/bvxlwkA2Mymp1zc/a6ZPSvpNUnbJL3g7lfMbKDQ/7ykz0r6GUlfMDNJuuvuXeUrGwBQzNxLng4vu66uLp+cnKzIvgGgVpnZpfUWzLxTFAACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHSgzp09e1bNzc3q7OxUa2urxsbGKl0S7hOBDtS5XC6nU6dOaWpqSuPj4zp27FilS8J9ItCBOpfP59Xe3i5Jamlp0dLSUoUrwv0i0IE6l8/n1dbWJnfX8PCwDh06VOmScJ+iPIIOQKBmZ2e1sLCg3t5ebd++Xd3d3RoZGdG5c+f0yiuv6MaNGxocHNSHP/xhpdNpzc/Pa3x8vNJlYx0EOlDHcrmc9u/fr4mJiVXtfX196uvr09zcnI4fP67R0VGNjo7qyJEjFaoUUXDKBahj+XxeHR0d6/afPn1ag4ODW1gRHgSBDtSxfD6vPXv2rGl3d504cUIHDx7U3r17K1AZ7genXIA6ls1mS7Y/99xzOn/+vObn5zU9Pa2nn35aQ0NDunz5ss6cOaPPfOYzW1wpouCZogBQQ3imKADUAQIdAAJBoANAIAh0AAgEgQ4AgSDQgYBls1klk0k1NDQomUyue5siwsB96ECgstmsUqmUFhcXJUkzMzNKpVKSpP7+/kqWhjLhPnSghvX09Kzbd/HiRd2+fXtN+yOPPKJ9+/aV/J4LFy7EVBnKhfvQgTpUKsw3akft45QLUMM2WlEnk0nNzMysaU8kEqzEA8UKHQhUOp1WY2PjqrbGxkal0+kKVYRyI9CBQPX39yuTySiRSMjMlEgklMlkuCAaMC6KAkAN4aIoANQBAh0AAkGgA0AgCHQACASBDgCBqNhdLmZ2U9Ladz3E61FJ75R5H+XGMVRerdcvcQzVIo5jSLh7U6mOigX6VjCzyfVu76kVHEPl1Xr9EsdQLcp9DJxyAYBAEOgAEIjQAz1T6QJiwDFUXq3XL3EM1aKsxxD0OXQAqCehr9ABoG4Q6AAQiCAC3cwOmNlbZjZtZidL9JuZDRf6c2a2txJ1ridC/T1mNm9mU4Wvz1aizo2Y2QtmdsPMvr1Of1XPgRTpGKp6HszsMTP7VzN708yumNmnS4yp6nmIeAxVOw9m9hNm9u9m9q1C/X9ZYkz55sDda/pL0jZJ/ympVdLDkr4laXfRmKckfVWSSdon6RuVrvse6++R9JVK17rJcfy2pL2Svr1Of9XOwT0cQ1XPg6QPSdpbeP0BSd+tpd+FeziGqp2Hwp/r+wuvt0v6hqR9WzUHIazQuyVNu/tVd78j6UVJh4vGHJb0d77soqQPmtmHtrrQdUSpv+q5+9ck/WCDIdU8B5IiHUNVc/fvu/s3C6//V9KbknYWDavqeYh4DFWr8Oe6UNjcXvgqvvOkbHMQQqDvlDS7Yvu61v4FiDKmUqLW9uuF/8Z91cx+cWtKi1U1z8G9qIl5MLOkpF/W8gpxpZqZhw2OQarieTCzbWY2JemGpH9x9y2bgxAeEm0l2or/RYwyplKi1PZNLX9+w4KZPSXpnKQnyl1YzKp5DqKqiXkws/dL+rKkP3X3HxZ3l/iWqpuHTY6hqufB3ZckdZrZByX9o5n9kruvvC5TtjkIYYV+XdJjK7ZbJL19H2MqZdPa3P2H7/43zt1flbTdzB7duhJjUc1zEEktzIOZbddyEGbd/aUSQ6p+HjY7hlqYB0ly9/+RdEHSgaKuss1BCIH+hqQnzOxxM3tY0jOSXi4a87KkTxSuLu+TNO/u39/qQtexaf1m1mxmVnjdreV5u7XllT6Yap6DSKp9Hgq1jUp6093/Zp1hVT0PUY6hmufBzJoKK3OZ2U9K+j1J/1E0rGxzUPOnXNz9rpk9K+k1Ld8x8oK7XzGzgUL/85Je1fKV5WlJi5I+Wal6i0Ws/4ikPzazu5L+T9IzXrhcXi3M7O+1fPfBo2Z2XdJfaPmCUNXPwbsiHEO1z8NvSPq4pHzhHK4k/ZmkXVLNzEOUY6jmefiQpC+Z2TYt/0PzD+7+la3KI976DwCBCOGUCwBABDoABINAB4BAEOgAEAgCHQACQaADQCAIdAAIxP8D7Z/bpk/EFgwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "posterior = estimator.fit(\n",
    "    observation_trajectory, \n",
    "    n_burn_in=100,  # number of sampling steps to discard at the beginning \n",
    "    n_thin=10  # take every 10th sample\n",
    ").fetch_model()\n",
    "stats = posterior.gather_stats(\"transition_model/transition_matrix\", confidence=0.95)\n",
    "\n",
    "x = [0, 1, 2, 3]\n",
    "y = stats.mean.flatten()\n",
    "errors = np.concatenate((stats.mean - stats.L, -stats.mean + stats.R)).T\n",
    "f, ax = plt.subplots()\n",
    "ax.errorbar(x, y, xerr=errors, fmt='o', color=\"k\")\n",
    "ax.annotate(r\"$P_{11}$\", (x[0]+.05, y[0]+.05))\n",
    "ax.annotate(r\"$P_{12}$\", (x[1]+.05, y[1]+.05))\n",
    "ax.annotate(r\"$P_{21}$\", (x[2]+.05, y[2]+.05))\n",
    "ax.annotate(r\"$P_{22}$\", (x[3]-.1, y[3]-.1));"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. rubric:: References\n",
    "\n",
    ".. footbibliography::"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.9"
  },
  "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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}