{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# meta-data does not work yet in VScode\n", "# https://github.com/microsoft/vscode-jupyter/issues/1121\n", "\n", "{\n", " \"tags\": [\n", " \"hide-cell\"\n", " ]\n", "}\n", "\n", "\n", "### Install necessary libraries\n", "\n", "try:\n", " import jax\n", "except:\n", " # For cuda version, see https://github.com/google/jax#installation\n", " %pip install --upgrade \"jax[cpu]\" \n", " import jax\n", "\n", "try:\n", " import distrax\n", "except:\n", " %pip install --upgrade distrax\n", " import distrax\n", "\n", "try:\n", " import jsl\n", "except:\n", " %pip install git+https://github.com/probml/jsl\n", " import jsl\n", "\n", "try:\n", " import rich\n", "except:\n", " %pip install rich\n", " import rich\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "{\n", " \"tags\": [\n", " \"hide-cell\"\n", " ]\n", "}\n", "\n", "\n", "### Import standard libraries\n", "\n", "import abc\n", "from dataclasses import dataclass\n", "import functools\n", "import itertools\n", "\n", "from typing import Any, Callable, NamedTuple, Optional, Union, Tuple\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "\n", "import jax\n", "import jax.numpy as jnp\n", "from jax import lax, vmap, jit, grad\n", "from jax.scipy.special import logit\n", "from jax.nn import softmax\n", "from functools import partial\n", "from jax.random import PRNGKey, split\n", "\n", "import inspect\n", "import inspect as py_inspect\n", "import rich\n", "from rich import inspect as r_inspect\n", "from rich import print as r_print\n", "\n", "def print_source(fname):\n", " r_print(py_inspect.getsource(fname))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{math}\n", "\n", "\\newcommand\\floor[1]{\\lfloor#1\\rfloor}\n", "\n", "\\newcommand{\\real}{\\mathbb{R}}\n", "\n", "% Numbers\n", "\\newcommand{\\vzero}{\\boldsymbol{0}}\n", "\\newcommand{\\vone}{\\boldsymbol{1}}\n", "\n", "% Greek https://www.latex-tutorial.com/symbols/greek-alphabet/\n", "\\newcommand{\\valpha}{\\boldsymbol{\\alpha}}\n", "\\newcommand{\\vbeta}{\\boldsymbol{\\beta}}\n", "\\newcommand{\\vchi}{\\boldsymbol{\\chi}}\n", "\\newcommand{\\vdelta}{\\boldsymbol{\\delta}}\n", "\\newcommand{\\vDelta}{\\boldsymbol{\\Delta}}\n", "\\newcommand{\\vepsilon}{\\boldsymbol{\\epsilon}}\n", "\\newcommand{\\vzeta}{\\boldsymbol{\\zeta}}\n", "\\newcommand{\\vXi}{\\boldsymbol{\\Xi}}\n", "\\newcommand{\\vell}{\\boldsymbol{\\ell}}\n", "\\newcommand{\\veta}{\\boldsymbol{\\eta}}\n", "%\\newcommand{\\vEta}{\\boldsymbol{\\Eta}}\n", "\\newcommand{\\vgamma}{\\boldsymbol{\\gamma}}\n", "\\newcommand{\\vGamma}{\\boldsymbol{\\Gamma}}\n", "\\newcommand{\\vmu}{\\boldsymbol{\\mu}}\n", "\\newcommand{\\vmut}{\\boldsymbol{\\tilde{\\mu}}}\n", "\\newcommand{\\vnu}{\\boldsymbol{\\nu}}\n", "\\newcommand{\\vkappa}{\\boldsymbol{\\kappa}}\n", "\\newcommand{\\vlambda}{\\boldsymbol{\\lambda}}\n", "\\newcommand{\\vLambda}{\\boldsymbol{\\Lambda}}\n", "\\newcommand{\\vLambdaBar}{\\overline{\\vLambda}}\n", "%\\newcommand{\\vnu}{\\boldsymbol{\\nu}}\n", "\\newcommand{\\vomega}{\\boldsymbol{\\omega}}\n", "\\newcommand{\\vOmega}{\\boldsymbol{\\Omega}}\n", "\\newcommand{\\vphi}{\\boldsymbol{\\phi}}\n", "\\newcommand{\\vvarphi}{\\boldsymbol{\\varphi}}\n", "\\newcommand{\\vPhi}{\\boldsymbol{\\Phi}}\n", "\\newcommand{\\vpi}{\\boldsymbol{\\pi}}\n", "\\newcommand{\\vPi}{\\boldsymbol{\\Pi}}\n", "\\newcommand{\\vpsi}{\\boldsymbol{\\psi}}\n", "\\newcommand{\\vPsi}{\\boldsymbol{\\Psi}}\n", "\\newcommand{\\vrho}{\\boldsymbol{\\rho}}\n", "\\newcommand{\\vtheta}{\\boldsymbol{\\theta}}\n", "\\newcommand{\\vthetat}{\\boldsymbol{\\tilde{\\theta}}}\n", "\\newcommand{\\vTheta}{\\boldsymbol{\\Theta}}\n", "\\newcommand{\\vsigma}{\\boldsymbol{\\sigma}}\n", "\\newcommand{\\vSigma}{\\boldsymbol{\\Sigma}}\n", "\\newcommand{\\vSigmat}{\\boldsymbol{\\tilde{\\Sigma}}}\n", "\\newcommand{\\vsigmoid}{\\vsigma}\n", "\\newcommand{\\vtau}{\\boldsymbol{\\tau}}\n", "\\newcommand{\\vxi}{\\boldsymbol{\\xi}}\n", "\n", "\n", "% Lower Roman (Vectors)\n", "\\newcommand{\\va}{\\mathbf{a}}\n", "\\newcommand{\\vb}{\\mathbf{b}}\n", "\\newcommand{\\vBt}{\\mathbf{\\tilde{B}}}\n", "\\newcommand{\\vc}{\\mathbf{c}}\n", "\\newcommand{\\vct}{\\mathbf{\\tilde{c}}}\n", "\\newcommand{\\vd}{\\mathbf{d}}\n", "\\newcommand{\\ve}{\\mathbf{e}}\n", "\\newcommand{\\vf}{\\mathbf{f}}\n", "\\newcommand{\\vg}{\\mathbf{g}}\n", "\\newcommand{\\vh}{\\mathbf{h}}\n", "%\\newcommand{\\myvh}{\\mathbf{h}}\n", "\\newcommand{\\vi}{\\mathbf{i}}\n", "\\newcommand{\\vj}{\\mathbf{j}}\n", "\\newcommand{\\vk}{\\mathbf{k}}\n", "\\newcommand{\\vl}{\\mathbf{l}}\n", "\\newcommand{\\vm}{\\mathbf{m}}\n", "\\newcommand{\\vn}{\\mathbf{n}}\n", "\\newcommand{\\vo}{\\mathbf{o}}\n", "\\newcommand{\\vp}{\\mathbf{p}}\n", "\\newcommand{\\vq}{\\mathbf{q}}\n", "\\newcommand{\\vr}{\\mathbf{r}}\n", "\\newcommand{\\vs}{\\mathbf{s}}\n", "\\newcommand{\\vt}{\\mathbf{t}}\n", "\\newcommand{\\vu}{\\mathbf{u}}\n", "\\newcommand{\\vv}{\\mathbf{v}}\n", "\\newcommand{\\vw}{\\mathbf{w}}\n", "\\newcommand{\\vws}{\\vw_s}\n", "\\newcommand{\\vwt}{\\mathbf{\\tilde{w}}}\n", "\\newcommand{\\vWt}{\\mathbf{\\tilde{W}}}\n", "\\newcommand{\\vwh}{\\hat{\\vw}}\n", "\\newcommand{\\vx}{\\mathbf{x}}\n", "%\\newcommand{\\vx}{\\mathbf{x}}\n", "\\newcommand{\\vxt}{\\mathbf{\\tilde{x}}}\n", "\\newcommand{\\vy}{\\mathbf{y}}\n", "\\newcommand{\\vyt}{\\mathbf{\\tilde{y}}}\n", "\\newcommand{\\vz}{\\mathbf{z}}\n", "%\\newcommand{\\vzt}{\\mathbf{\\tilde{z}}}\n", "\n", "\n", "% Upper Roman (Matrices)\n", "\\newcommand{\\vA}{\\mathbf{A}}\n", "\\newcommand{\\vB}{\\mathbf{B}}\n", "\\newcommand{\\vC}{\\mathbf{C}}\n", "\\newcommand{\\vD}{\\mathbf{D}}\n", "\\newcommand{\\vE}{\\mathbf{E}}\n", "\\newcommand{\\vF}{\\mathbf{F}}\n", "\\newcommand{\\vG}{\\mathbf{G}}\n", "\\newcommand{\\vH}{\\mathbf{H}}\n", "\\newcommand{\\vI}{\\mathbf{I}}\n", "\\newcommand{\\vJ}{\\mathbf{J}}\n", "\\newcommand{\\vK}{\\mathbf{K}}\n", "\\newcommand{\\vL}{\\mathbf{L}}\n", "\\newcommand{\\vM}{\\mathbf{M}}\n", "\\newcommand{\\vMt}{\\mathbf{\\tilde{M}}}\n", "\\newcommand{\\vN}{\\mathbf{N}}\n", "\\newcommand{\\vO}{\\mathbf{O}}\n", "\\newcommand{\\vP}{\\mathbf{P}}\n", "\\newcommand{\\vQ}{\\mathbf{Q}}\n", "\\newcommand{\\vR}{\\mathbf{R}}\n", "\\newcommand{\\vS}{\\mathbf{S}}\n", "\\newcommand{\\vT}{\\mathbf{T}}\n", "\\newcommand{\\vU}{\\mathbf{U}}\n", "\\newcommand{\\vV}{\\mathbf{V}}\n", "\\newcommand{\\vW}{\\mathbf{W}}\n", "\\newcommand{\\vX}{\\mathbf{X}}\n", "%\\newcommand{\\vXs}{\\vX_{\\vs}}\n", "\\newcommand{\\vXs}{\\vX_{s}}\n", "\\newcommand{\\vXt}{\\mathbf{\\tilde{X}}}\n", "\\newcommand{\\vY}{\\mathbf{Y}}\n", "\\newcommand{\\vZ}{\\mathbf{Z}}\n", "\\newcommand{\\vZt}{\\mathbf{\\tilde{Z}}}\n", "\\newcommand{\\vzt}{\\mathbf{\\tilde{z}}}\n", "\n", "\n", "%%%%\n", "\\newcommand{\\hidden}{\\vz}\n", "\\newcommand{\\hid}{\\hidden}\n", "\\newcommand{\\observed}{\\vy}\n", "\\newcommand{\\obs}{\\observed}\n", "\\newcommand{\\inputs}{\\vu}\n", "\\newcommand{\\input}{\\inputs}\n", "\n", "\\newcommand{\\hmmTrans}{\\vA}\n", "\\newcommand{\\hmmObs}{\\vB}\n", "\\newcommand{\\hmmInit}{\\vpi}\n", "\\newcommand{\\hmmhid}{\\hidden}\n", "\\newcommand{\\hmmobs}{\\obs}\n", "\n", "\\newcommand{\\ldsDyn}{\\vA}\n", "\\newcommand{\\ldsObs}{\\vC}\n", "\\newcommand{\\ldsDynIn}{\\vB}\n", "\\newcommand{\\ldsObsIn}{\\vD}\n", "\\newcommand{\\ldsDynNoise}{\\vQ}\n", "\\newcommand{\\ldsObsNoise}{\\vR}\n", "\n", "\\newcommand{\\ssmDynFn}{f}\n", "\\newcommand{\\ssmObsFn}{h}\n", "\n", "\n", "%%%\n", "\\newcommand{\\gauss}{\\mathcal{N}}\n", "\n", "\\newcommand{\\diag}{\\mathrm{diag}}\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(sec:hmm-intro)=\n", "# Hidden Markov Models\n", "\n", "In this section, we discuss the\n", "hidden Markov model or HMM,\n", "which is a state space model in which the hidden states\n", "are discrete, so $\\hmmhid_t \\in \\{1,\\ldots, K\\}$.\n", "The observations may be discrete,\n", "$\\hmmobs_t \\in \\{1,\\ldots, C\\}$,\n", "or continuous,\n", "$\\hmmobs_t \\in \\real^D$,\n", "or some combination,\n", "as we illustrate below.\n", "More details can be found in e.g., \n", "{cite}`Rabiner89,Fraser08,Cappe05`.\n", "For an interactive introduction,\n", "see https://nipunbatra.github.io/hmm/." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(sec:casino)=\n", "### Example: Casino HMM\n", "\n", "To illustrate HMMs with categorical observation model,\n", "we consider the \"Ocassionally dishonest casino\" model from {cite}`Durbin98`.\n", "There are 2 hidden states, representing whether the dice being used in the casino is fair or loaded.\n", "Each state defines a distribution over the 6 possible observations.\n", "\n", "The transition model is denoted by\n", "```{math}\n", "p(z_t=j|z_{t-1}=i) = \\hmmTrans_{ij}\n", "```\n", "Here the $i$'th row of $\\vA$ corresponds to the outgoing distribution from state $i$.\n", "This is a row stochastic matrix,\n", "meaning each row sums to one.\n", "We can visualize\n", "the non-zero entries in the transition matrix by creating a state transition diagram,\n", "as shown in \n", "{numref}`fig:casino`.\n", "\n", "```{figure} /figures/casino.png\n", ":scale: 50%\n", ":name: fig:casino\n", "\n", "Illustration of the casino HMM.\n", "```\n", "\n", "The observation model\n", "$p(\\obs_t|\\hidden_t=j)$ has the form\n", "```{math}\n", "p(\\obs_t=k|\\hidden_t=j) = \\hmmObs_{jk} \n", "```\n", "This is represented by the histograms associated with each\n", "state in {numref}`casino-fig`.\n", "\n", "Finally,\n", "the initial state distribution is denoted by\n", "```{math}\n", "p(z_1=j) = \\hmmInit_j\n", "```\n", "\n", "Collectively we denote all the parameters by $\\vtheta=(\\hmmTrans, \\hmmObs, \\hmmInit)$.\n", "\n", "Now let us implement this model in code." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# state transition matrix\n", "A = np.array([\n", " [0.95, 0.05],\n", " [0.10, 0.90]\n", "])\n", "\n", "# observation matrix\n", "B = np.array([\n", " [1/6, 1/6, 1/6, 1/6, 1/6, 1/6], # fair die\n", " [1/10, 1/10, 1/10, 1/10, 1/10, 5/10] # loaded die\n", "])\n", "\n", "pi = np.array([0.5, 0.5])\n", "\n", "(nstates, nobs) = np.shape(B)\n" ] }, { "cell_type": "code", "execution_count": 4, "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" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "import distrax\n", "from distrax import HMM\n", "\n", "\n", "hmm = HMM(trans_dist=distrax.Categorical(probs=A),\n", " init_dist=distrax.Categorical(probs=pi),\n", " obs_dist=distrax.Categorical(probs=B))\n", "\n", "print(hmm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Let's sample from the model. We will generate a sequence of latent states, $\\hid_{1:T}$,\n", "which we then convert to a sequence of observations, $\\obs_{1:T}$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Printing sample observed/latent...\n", "x: 633665342652353616444236412331351246651613325161656366246242\n", "z: 222222211111111111111111111111111111111222111111112222211111\n" ] } ], "source": [ "\n", "\n", "\n", "seed = 314\n", "n_samples = 300\n", "z_hist, x_hist = hmm.sample(seed=PRNGKey(seed), seq_len=n_samples)\n", "\n", "z_hist_str = \"\".join((np.array(z_hist) + 1).astype(str))[:60]\n", "x_hist_str = \"\".join((np.array(x_hist) + 1).astype(str))[:60]\n", "\n", "print(\"Printing sample observed/latent...\")\n", "print(f\"x: {x_hist_str}\")\n", "print(f\"z: {z_hist_str}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below is the source code for the sampling algorithm.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
  def sample(self,\n",
       "             *,\n",
       "             seed: chex.PRNGKey,\n",
       "             seq_len: chex.Array) -> Tuple:\n",
       "    \"\"\"Sample from this HMM.\n",
       "\n",
       "    Samples an observation of given length according to this\n",
       "    Hidden Markov Model and gives the sequence of the hidden states\n",
       "    as well as the observation.\n",
       "\n",
       "    Args:\n",
       "      seed: Random key of shape (2,) and dtype uint32.\n",
       "      seq_len: The length of the observation sequence.\n",
       "\n",
       "    Returns:\n",
       "      Tuple of hidden state sequence, and observation sequence.\n",
       "    \"\"\"\n",
       "    rng_key, rng_init = jax.random.split(seed)\n",
       "    initial_state = self._init_dist.sample(seed=rng_init)\n",
       "\n",
       "    def draw_state(prev_state, key):\n",
       "      state = self._trans_dist.sample(seed=key)\n",
       "      return state, state\n",
       "\n",
       "    rng_state, rng_obs = jax.random.split(rng_key)\n",
       "    keys = jax.random.split(rng_state, seq_len - 1)\n",
       "    _, states = jax.lax.scan(draw_state, initial_state, keys)\n",
       "    states = jnp.append(initial_state, states)\n",
       "\n",
       "    def draw_obs(state, key):\n",
       "      return self._obs_dist.sample(seed=key)\n",
       "\n",
       "    keys = jax.random.split(rng_obs, seq_len)\n",
       "    obs_seq = jax.vmap(draw_obs, in_axes=(0, 0))(states, keys)\n",
       "\n",
       "    return states, obs_seq\n",
       "\n",
       "
\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print_source(hmm.sample)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us check correctness by computing empirical pairwise statistics\n", "\n", "We will compute the number of i->j latent state transitions, and check that it is close to the true \n", "A[i,j] transition probabilites." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[225. 92.]\n", " [ 92. 90.]]\n", "[[0.7097792 0.29022083]\n", " [0.50549453 0.4945055 ]]\n" ] } ], "source": [ "\n", "\n", "import collections\n", "def compute_counts(state_seq, nstates):\n", " wseq = np.array(state_seq)\n", " word_pairs = [pair for pair in zip(wseq[:-1], wseq[1:])]\n", " counter_pairs = collections.Counter(word_pairs)\n", " counts = np.zeros((nstates, nstates))\n", " for (k,v) in counter_pairs.items():\n", " counts[k[0], k[1]] = v\n", " return counts\n", "\n", "\n", "def normalize(u, axis=0, eps=1e-15):\n", " u = jnp.where(u == 0, 0, jnp.where(u < eps, eps, u))\n", " c = u.sum(axis=axis)\n", " c = jnp.where(c == 0, 1, c)\n", " return u / c, c\n", "\n", "def normalize_counts(counts):\n", " ncounts = vmap(lambda v: normalize(v)[0], in_axes=0)(counts)\n", " return ncounts\n", "\n", "init_dist = jnp.array([1.0, 0.0])\n", "trans_mat = jnp.array([[0.7, 0.3], [0.5, 0.5]])\n", "obs_mat = jnp.eye(2)\n", "\n", "hmm = HMM(trans_dist=distrax.Categorical(probs=trans_mat),\n", " init_dist=distrax.Categorical(probs=init_dist),\n", " obs_dist=distrax.Categorical(probs=obs_mat))\n", "\n", "rng_key = jax.random.PRNGKey(0)\n", "seq_len = 500\n", "state_seq, _ = hmm.sample(seed=PRNGKey(seed), seq_len=seq_len)\n", "\n", "counts = compute_counts(state_seq, nstates=2)\n", "print(counts)\n", "\n", "trans_mat_empirical = normalize_counts(counts)\n", "print(trans_mat_empirical)\n", "\n", "assert jnp.allclose(trans_mat, trans_mat_empirical, atol=1e-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our primary goal will be to infer the latent state from the observations,\n", "so we can detect if the casino is being dishonest or not. This will\n", "affect how we choose to gamble our money.\n", "We discuss various ways to perform this inference below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(sec:lillypad)=\n", "## Example: Lillypad HMM\n", "\n", "\n", "If $\\obs_t$ is continuous, it is common to use a Gaussian\n", "observation model:\n", "```{math}\n", "p(\\obs_t|\\hidden_t=j) = \\gauss(\\obs_t|\\vmu_j,\\vSigma_j)\n", "```\n", "This is sometimes called a Gaussian HMM.\n", "\n", "As a simple example, suppose we have an HMM with 3 hidden states,\n", "each of which generates a 2d Gaussian.\n", "We can represent these Gaussian distributions are 2d ellipses,\n", "as we show below.\n", "We call these ``lilly pads'', because of their shape.\n", "We can imagine a frog hopping from one lilly pad to another.\n", "(This analogy is due to the late Sam Roweis.)\n", "The frog will stay on a pad for a while (corresponding to remaining in the same\n", "discrete state $\\hidden_t$), and then jump to a new pad\n", "(corresponding to a transition to a new state).\n", "The data we see are just the 2d points (e.g., water droplets)\n", "coming from near the pad that the frog is currently on.\n", "Thus this model is like a Gaussian mixture model,\n", "in that it generates clusters of observations,\n", "except now there is temporal correlation between the data points.\n", "\n", "Let us now illustrate this model in code.\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# Let us create the model\n", "\n", "initial_probs = jnp.array([0.3, 0.2, 0.5])\n", "\n", "# transition matrix\n", "A = jnp.array([\n", "[0.3, 0.4, 0.3],\n", "[0.1, 0.6, 0.3],\n", "[0.2, 0.3, 0.5]\n", "])\n", "\n", "# Observation model\n", "mu_collection = jnp.array([\n", "[0.3, 0.3],\n", "[0.8, 0.5],\n", "[0.3, 0.8]\n", "])\n", "\n", "S1 = jnp.array([[1.1, 0], [0, 0.3]])\n", "S2 = jnp.array([[0.3, -0.5], [-0.5, 1.3]])\n", "S3 = jnp.array([[0.8, 0.4], [0.4, 0.5]])\n", "cov_collection = jnp.array([S1, S2, S3]) / 60\n", "\n", "\n", "import tensorflow_probability as tfp\n", "\n", "if False:\n", " hmm = HMM(trans_dist=distrax.Categorical(probs=A),\n", " init_dist=distrax.Categorical(probs=initial_probs),\n", " obs_dist=distrax.MultivariateNormalFullCovariance(\n", " loc=mu_collection, covariance_matrix=cov_collection))\n", "else:\n", " hmm = HMM(trans_dist=distrax.Categorical(probs=A),\n", " init_dist=distrax.Categorical(probs=initial_probs),\n", " obs_dist=distrax.as_distribution(\n", " tfp.substrates.jax.distributions.MultivariateNormalFullCovariance(loc=mu_collection,\n", " covariance_matrix=cov_collection)))\n", "\n", "print(hmm)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(50,)\n", "(50, 2)\n" ] } ], "source": [ "\n", "n_samples, seed = 50, 10\n", "samples_state, samples_obs = hmm.sample(seed=PRNGKey(seed), seq_len=n_samples)\n", "\n", "print(samples_state.shape)\n", "print(samples_obs.shape)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACpOElEQVR4nOydd3gc5bm375ntTb33LtmWe5U7zXQIAUIgCenASc8hJyGEnIQcyJdOGgkppAIhISS0AKG7997Ue1+V7X1mvj9WO5YsyZZt2djO3telS9rd2ZnZ1e5v3vd5n+f3CIqiECdOnDhxLnzEd/sE4sSJEyfO9BAX9Dhx4sS5SIgLepw4ceJcJMQFPU6cOHEuEuKCHidOnDgXCdp368BpaWlKUVHRu3X4OHHixLkg2b1794CiKOkTPfauCXpRURG7du16tw4fJ06cOBckgiC0TfZYPOQSJ06cOBcJcUGPEydOnIuEuKDHiRMnzkVCXNDjxIkT5yIhLuhx4sSJc5EQF/Q4ceLEuUiIC3qcOHHiXCTEBT1OnDhxLhLigh4nTpw4FwlxQY8TJ06ci4S4oMeJEyfORUJc0OPEiRPnIiEu6HHixIlzkRAX9Dhx4sS5SIgLepw4ceJcJMQFPU6cOHEuEuKCHidOnDgXCXFBjxMnTpyLhLigx4kTJ85FQlzQ48SJE+ciIS7oceLEiXOREBf0OHHixLlIiAt6nDhx4lwkxAU9Tpw4cS4S4oIeJ06cOBcJcUGPEydOnIuEuKDHiRMnzkWC9t0+gThxpgs5EEByupAcw0jDDmSvB9njQfJ4UAIB5GAQJRBEiURAlkGRURQFQaNF0GoRdFoEgxHRakG0WNDYbGhSUtCmpqJJSUW0mBEE4d1+mXHiTEpc0OOc90huN+GuLsLd3YR7eojY7UT67UTsdqTBQSLDw0hDQyjB4Fk9D9FiQZeTjTYrG31+HvqyMgxlZRjKy9EmJ5/VY8eJMxXigh7nvED2egl1dhJqbibY1EyotZVQWxuhtjZkl+vdPj0geo7BhkaCDY14j3tMX1iIadFCzAsXYalZhi47+105xzj/2ZxU0AVB+B1wHdCvKEr1BI8LwE+AawAf8BFFUfZM94nGuTiQXC4CtbUE6+oJtYyId3MzEbv9zHeu06FJSECbkoyYmIjGloBotUZDKCYzgkGPaDAgaLUgakAcCZ9IEko4ghKJoAT8SB4PsteH5HIiDQ0TGRxAGjzxDCB28XE++w8AjLNmYbvicmyXX46hrOzMX1ucOFNgKiP0PwA/B/40yeNXA+UjP0uBX478jvMfjKIoRPrtBOvrCNTWEjhwkMDhw4S7u095X4LBgC4nZ+QnG21GJtqMDLTp6WjTUtEkJ6NJTjmrMW5FUZAcDiI9PYS7uwm1thJsaBgZsTeghMNjtg8cPkzg8GHsP/4JprlzSbr9/dguuYRgfT2C2YxxxgwEMZ6TEGd6OamgK4qyQRCEohNsciPwJ0VRFGCbIAhJgiBkK4rSM10nGef8R3I48O3di3//fgJHjhA4chRpYGDqO9Dp0OfkoCsswFBWjr64CH1hIfrCIrQZ6e/6YqQgCGiTk9EmJ2OcOXPMY3IwSODQIXw7d+HbsR3vjp0QiaiP+/fvx79/Pz1aLZqkpGgsPiuLnO99D11mxrl+KXEuYqYjhp4LdIy63Tly3zhBFwThLuAugIKCgmk4dJx3A0WSCDY1ETh4EP/+A/j27CbU2DSl5wo6HfqyMoyVlRjKy9GXFGMoKUGXmxsNhVyAiAYD5oULMS9cCPfcjeRy4Vm/Affrr+N++22Ijd4jEaSBAWSPBzkQoP/73yP3Bz94d08+zkXFdHyDJho6KRNtqCjKr4FfAyxatGjCbeKcfyjhMP6Dh/Dt2oVv5078e/Yge49fFhyPaDZjqKzEUFGBsXoWptmzMZSWIuh05+Cs3z00CQkkXn8diddfR2RwkP7vfx/nCy9GUyUBJRBACgTwbt6C7PMhms3v8hnHuViYDkHvBPJH3c4DTj1QGue8QRXw3bvw796Db+fOkwu4Votx5kzM8+djnDMb48yZ6AsL/+PjxNrUVLRXXIG0dSuKz4/G7VZHQNLwMP0/eoTMr93/roeU4lwcTIegvwB8RhCEp4kuhjrj8fMLj1BrK54NG/Fu3Ypvx46TCrg2PR3j3DmYqmdjmjcP09w5iCbTOTrbCwOPx0NdXR3dLhd5Oj16i4guI4NIRweEQgAMP/EE2vR00u6+610+2zgXA1NJW/wLsBZIEwShE/gGoANQFOUx4GWiKYuNRNMWP3q2TjbO9CG53Xg3b8G7ZQve7dsIt7WfcHtdTg7mZcswL1qEecniaMz7JKPKK6+8knXr1nHvvfdO56mf9/j9furq6ujs7ERRFGzJyeR+/3sEfv4ooY52dJmZoCiEOzsBsD/yCIbyMmyXXvoun3mcCx0hmpxy7lm0aJGya9eud+XY/6mE2tpwv/EmnvXr8e3ZMyYT43i02dlYli7FvGgh5kWL0BUWnnJYICMjg6eeeorLL7/8tM/5a1/7Gk899RSDg4MYjUZWr17Nj370o/NyUT0YDNLQ0EBbWxuyLKPRaKisrKS4uBhRFKOpjwMDCEYjosFA+8c/gW/nTiB6wSx99RUEvf5dfhVxzncEQditKMqiCR+LC/rFixKJ4Nu1C/dbb+HdsJFQa+uk2womE5aaGqyrVmKpqTktAR9NV1cXeXl52O120tLSTns/tbW1ZGdnk5iYiM/n44EHHmDbtm1s2bLltPc53YTDYRobG2lpaUGSJACysrKorq7GdIIwVGR4mOZrrkUaHo4+55vfIPn97z8n5xznwuVEgn5h5onFmRQ5FMK3bRuu117D88abSA7HpNsaZ87EsmY1lpoaTPPmIU7j6HDv3r3k5uaqYv7cc8/xxS9+ke9+97u8733vm/J+qqqq1L8VRUEURerq6qbtPM+ESCRCS0sLTU1NhEdSE81mM9XV1WRmZp70+drkZFI/8Qn6v/99AJzPvxAX9DhnRFzQLwIkjwfPW2/hfuttvBs3TrqgKRiNWGpqsF12KZbVq9FlnL2ilj179jBv3jyCwSBf+tKX2LJlC6+//jplZWV86lOf4qmnnpr0uffddx/33Xefevupp57iv/7rv3C5XGi1Wn70ox+dtfOeCrIs09raSkNDA6GRxU1RFCktLaW8vByNRjPlfSXe9B76f/ADUBT8+/cjORxokpLO0pnHudiJh1wuUORAAM+GDbhefgXP229P6jOizcjAdsUVWC+5BPOihYhG4zk5v5tuugmtVktjYyM1NTU88sgjGAyGM9pnb28vjz/+OCtWrGDt2rXTc6KngKIodHR0UF9fj9/vV+9PTU1lzpw5WK3W09pv83vfS/DIUQAKn3wiWqAUJ84kxEMuFwlKKIRnyxZc/3oZz5tvIvt8E26ny8vDdvnl2K5ch2nu3HclF3zv3r0Eg0EqKyv5xS9+MS37zMrK4pOf/CQlJSW0t7eTkpIyLfs9GYqi0NPTQ21tLd5Rsx+DwcDMmTPJy8s7o/3rMrNUQY8MDp7RvuL8ZxMX9PMcRVHw79uH66V/4frXvyaNiRuqqki46ipsl12KvqzsXS1UGR4epq2tjYMHD3L77bfzwAMP8NBDD6mP33PPPTzxxBOTPv/+++/n/vvvn/CxSCSC1+ulu7v7nAh6X18ftbW1uI6z8C0qKqKqqgrddFS9vkuz5DgXH3FBP08JdXbh/Mc/cD733KQOhfqiImxXX0XiNddgKC8/x2c4OXv37iUtLY3q6mpeeeUVampqyMvL45577gHgscce47HHHjvpfmRZ5he/+AXve9/7yMjIoLOzk89+9rOqmJ5NBgcHOXr0KMMjGSgxEhMTmTNnDknTGOcO9/aqf5/NdY04Fz9xQT+PkAMB3K+9hvO55/Fu3TrhyE2bnU3itdeQcM01GGbMOC9Lxvfu3cuCBQsAyMvL4+WXX2bt2rVkZWXxnve855T29fLLL/Otb30Lr9dLUlISa9eu5Y033kB7loy8HA4HtbW12O129Ho9RqORYDCIRqOhqqqKoqKiaX3PJZeLYH199IYgoC8qmrZ9x/nPI74oeh4QqK3F8czfcb744oTdecTERBLWXUHCtddiXrLkP94f5Wzgdrupra2lt7cXvV5PUlISLpeLQCBAbm4uM2fOxHgWFpRd/36Nrs9/Hog2xSh+9u/Tfow4FxfxRdHzEMnjxfXKyzie+TuBAwfGbyAIWFauJOmWW7BdsjZeQXiW8Hq91NfX09nZiVarpbCwEL/fT39/PxaLhWXLlpGenn7Wju966UX1b8vqVWftOHH+M4gL+jkmUF/P8BNP4nzpJZQJslR0+fkk3XwziTfecEH1pVQkCd/u3YQ7OtEXF2Oa9+5k10yVQCBAfX097e3tiKJISUkJWq2W5uZmZFmmsrKSsrIyxLP4GkLt7bjfelu9nXj9DWftWHH+M4gL+jlACYdxv/46w0/9Bd9EYSadjoQrriDpfbdekCEV2e+n63++HI0FSxJoNJjmziXn4YfOu5lFKBRSy/QBCgsLSU9Pp76+HqfTSXp6OrNnz8ZisZz1cxn4xS+j7xdgXrYMQ0nxWT9mnIubuKCfRSSnE8czzzD0xJNERmUyxNCXlZJ8660k3HAD2uTkcY8rskywoQEAQ3n5eSv0zhdfJFhbi5iQgCAII6mWe3G9/gaJ117zbp8eEE13bG5upqmpCUmSyMvLo7i4mLa2Nnbu3InRaGThwoXk5OSck/Px7tiB87nn1Nvpn/7UOTlunIubuKCfBUKdnQz94Y84nn0WZVRFIQAaDbbLLiPlQx/EtGjRpBkTobY2uu+/f6TQJNrPMvvbD2MoPv9Gcd6NGxF0OmRFweN2o9Vq0SkK3s2bxwi65HIR7uxEm5094QXsbCBJEq2trTQ2NhIKhcjOzqayshKn08n27dsJhUKUlJRQWVl51jJnxp2Ty0XP1x5Qb1svuwzz4sXn5NhxLm7igj6N+A8fZujx3+F69VW13VgMTWoqye9/P0nve99JGwMrskzP1/+XiH0Ajc0GRCsIex74OoV//tN5N1LXpGeg1NYhAIoiEw6HkHx+hj0ehg8fJi8vD+nFFxl+4kkQBBRJIum9N5F6111n7bXIsqyW6QcCAdLT06mqqkKr1XLgwAEGBwdJTk5m2bJlJCQknJVzmIze/3uIcEe0Da9otZL1v18/p8ePc/ESF/QzRFEUfNu3M/jr3+CdwNLVUFFByofvJOG66xCn6GUSamsj3NeHaLXi9/tRAJPVSmRwkFBzM4aysml+FWdG8i034928GQIBREFEEw6js1gIX34Zra2tdL/+OlnPPovWakNvMiHKMo6/P4uhomLamzooikJXVxd1dXX4fD5SUlJYsGABSUlJNDQ00NTUhEajYc6cORQUFJzzPH7Hs//A9eKxzJbsbz0YbXhxluh1BtjROkS/K0BYUohI0YFGdpKJwlQzhSlm0m2G87KeIc6pExf000SRZTxvvcXAr35N4ODBcY9blteQ8rGPY1mx/JS/LIJGAyj4/X5C4RBGw0j+s6LAKTj5nSuMM2eS/a1vMfCrX6FrbMSfmkrJN/4X87x5hEIhWt9Zj1+WCYRDBMJRd0IzAs6X/jWtgt7b20ttbS1ut5uEhASWLl1KRkYGfX19vP322/j9fvLz85kxY8YZG4WdDr7du+n55jfV24k33kjCNdO/xtDY7+Z3m1vZ2jRIy8DJm3ln2AysKEtjcVEKS4pTKE23xAX+AiUu6KeIIkm4Xn2Vwcd+pS5YqogiCVddRcrHP4Zp1qzTPoY2L49AUhJSRyfGpCQMBgOy240uP/+8rSS0LF2CZekSuru7adq9m4LcXBzd3fT39xNyOjlePgOBAI6ebgb376ewsPCMSuntdju1tbU4HA6sVisLFy4kOzubQCDAzp076e3txWazsXz5clJTU8/odZ4u/kOH6bjnv2DEN91QWTntoZZAWOLRtxt5bH0TYWnqBYP97iD/3NvFP/d2AVCeYeX6uTlcNyebkvTTc5CM8+4QrxSdIoos43rlFQYe/QWh5uYxjwkGA0k3v5eUj34UfX7+GR1HlmX27t1L35EjlPz7NbT9/QDoiwrJ/uY30Z2jLIzTweVy0dHRQfOo90en05E1NITxN79FZ7Mh6nSEAgECw8MMvOdGfDNmANHGEIWFheTl5U25InNoaIja2loGBwcxmUxUVlaSl5eHoig0NzdTP1JSX1FRQUlJyVnNKT8RodZWWu/4ANLQEACalBSK/vY39Hm503aMQ11OPvPUHloHx9Y26LUiCwuSqcyyYdCK6DQiYVmmc8hP25CX1gEfnuDkrQgXFSbzwWWFXFWdhVF3/s0O/xOJt6A7AxRZxv36Gwz8/GcEGxrHPCaazSTfcTspH/kI2jNosxZDlmV2795Nb28vs2bNori4mEhPDygK2pyc824aHIlEGBgYoL+/n76+PgKBwJjHV6xYQXJyMoIgMPz00wz98U8gCCBLJN5yK/IN13PkyBHcbveY56Wnp5Ofn09WVtaEzSJcLlfUb6W1FVtHJzl5uRRdcw1am42hoSEOHDiA2+2eUhu4s02wpYX2D3+EyMiFWZOYSMEf/4BxGs3FGvrcvPcXW3CPEuYFBUl8aV0lCwqTTyjEkqxwqMvJjpYhdrYOsbFhAH9YGrddslnH+xbl8/FVxWTYzo2nfpyJiQv6aeLdsoX+H/6IwOHDY+4XrVZS7vwQyR/60LSl30mSxK5du+jv72f27NkUnaehFa/XS19fH/39/QwODiLLMlqtlvT0dDIyMsjIyKCpqYnW1lauvvrqMaNiyeMh3N2NLisLzajMkoGBAerr6xk8zgtcq9WSk5NDfn4+KSkpeL1e6urq6OrqwtreTtZzz6MTBBAEBJ0Oz0c/Qrtej8lkorq6mqysrHP2vkxEoLaWjk/eRcRuB6Idowp+9zvMC+ZP2zG8wQjX/nSjOjK3GbTcd00Vty8uQBRPfQDgC0V482g/L+zv5u3afiLyWH0w6kRuW5TPJ1aVkJ9inpbXEOfUiAv6KeI/fBj7D380LmtFNJtJ/vCdpH7kI2gSE6fteJIksXPnTux2O3Pnzj2vOtrLsszg4KAq4rEGD1arlczMTDIyMkhJSRkj3L29vezcufOUY9aDg4PU19czMDBwwu3KiorQf+0BkCREo5FQKETA6UA2GDH++BEqqqpOqQ3c2cC7Ywed//UptR2gYDKR/8tfYlm2dFqP88PX6vjZW9GZo0mn4Zl7aqjOnZ7PZr87wDO7OnlqeztdjrH1FBpR4H2L8vni5eVkJMRH7OeSuDnXFAn39WH/0SM4n39+zP2CwUDyBz5A6ic/Me0FMZFIhB07djA0NMT8+fPPuPvNyZAcDiS3G11Ozkg2zXj8fr8q4AMDA0iShEajIS0tjZKSEjIyMjCbJx+dpaamIggCAwMDpyToqamp1NTUMDQ0RH19PfaRke3xuA4dItHnQ5+QgNfnJRwOozWZMSkyhWbzuy7mzpf+Rc/Xvqa2BRStVvJ/+YtpLx7qcwX49YZj6xUP3jhr2sQcIMNm5NOXlHHPmlLePNrHz95q5GCXE4iGav6yo53n9nbx8ZXFfHJ1CYmmaWj2EeeMiAs6US+Swcd/x+Djj4+t7BRFkm5+L2mf/jS6szB9D4fDbN++HYfDwfz588nNnb5FsuORg0H6f/gjPO+8A4KAJiGBzPu/inn+fGRZZnh4WBXxWEzbbDZTUFBARkYGqampUxZKnU5HYmIiAwMDVFZWnvK52mw2kpKSxgl6QUEBOp0O+4EDGP1+/IoMgoBWq8VitiB7PIjnwINlMhRJwv7IIwz+9nH1Pm16Ovm//S3GyoppP97fd3cSjETzymfnJnLLgrMzGNCIAutmZXHFzEw2NQ7w6NuNbGuOLvD6wxI/f7uRp3a085WrKrl1Yf5phXriTA//0YKuKAru11+n7zvfIdLdM+Yx66WXknHvf2MoLT0rxw6Hw2zbtg2Xy8WiRYvOerx38Pe/x/PWW4g2G4IoEvG4af/KfXi+eh/2QIBIJIIoiqSmpqoiPlHTYyUUwvX6G3i3bkWbnkbiDTdMaEeQlpZGU1MTkUhkyiX1kiTR0tJCY2Mj4XCY3NxcKisriUQiqjOiVqulaO5cqKggVFeLbDQSCYXwuN3o5s4jaDbzbowTJZeLrnu/hHfjRvU+fXEx+b96DP1ZCqE9v69L/fvjK4vPupAKgsCq8nRWlqWxoWGA//fyUWp7oxf/IW+Irzx7kKe2t/PgjdXMy086q+cSZ2L+Y2PoodZWeh96GO+mTWPuN1RVkXnfV7AsW3b2jh0KsXXrVjweD4sWLSLzLFYKxmi6/gYEQUAWRXx+H5IkIfoDuG+4HttVV5GZmUlaWtoJxVeRZbq//BX8+/eBqEGRJUStjpzvfgfTnDljtrXb7Wzbtk0t7jkRsizT1tZGQ0MDwWCQzMxMqqqqxpXku1wu6uvr6enpQfT5SHn136S2t0fNwGbPpmvVSmSDgeTkZPLy8sjNzZ2enp8nIdjYSOenP0OorU29z7pmDTk/+L5q3TDdDHqCLHzoDQD0GpG9/3sFFsO5HZ/JssKLB7r53qt142LsH1xWwFeuqsJmjIdhppt4DH0USijE4OOPM/DLx1BCIfV+TXIy6V/8Akk33zxpbHk6CAaDbN26Fa/Xy5IlS85q84QxSBHQ6VEUJRoTFzUYDQaKZs4kce7cKe3Cv3cv/gMHEG0JKIAoCMgeDwO/fIz8X/5izLaxhdKBgYFJBV1RFDo7O6mrq8Pv95OamsqiRYtIsljw792LJxTCNH++KooJCQksWrSITZs2MQwMvPcmhhWFwsJCZlRVMQPo7Oyks7OTgwcPcvjwYTIzM8nPzyc9PX3a89AVRcHxzDP0ffv/oYxK2Uy9+27SP/fZs/o52t/pUP+uzk0452IOIIoCN87LZd3MLH7xTiO/2tBMaCQE9MS2dt482s/DN1VzadXZH7DEiTKlT4EgCFcBPwE0wG8VRfnOcY8nAk8ABSP7/IGiKL+f5nM9Y3x79tLzv18n1Nh07E5BIPn295P+uc+hmcbGvxMRCATYunUrfr+fpUuXkjYNuetTxXrppbhfeRVNYiIGvYGQz4us0WJeMvWsi1BrK0okgiAIuEda5Rn0emhqGretRqMhOTl5wowVRVHo6emhrq4Oj8dDUlISc+fOJT09nWBzM22fvAvJ6wVFQdBoyPz6A1hragDo7u5meHiY8vJy8vLyaGhooKWjg9bOTgoLCykrK6O0tBSn06mKe09PDwaDgdzcXPLz86fFjEtyu+n91v+N8WURTCZyvv0wCVdffcb7PxltowqIqrLPrbnY8Zj0Gu5dF42ff/PFw7xVG82573EG+NgfdvGeeTk8eEM1ieb4aP1sc1JBFwRBAzwKXAF0AjsFQXhBUZQjozb7NHBEUZTrBUFIB+oEQXhSUZTQBLs850geL/ZHHmH4qafGNF42zppF1oMPYqo+/TL9qeL3+9m6dSvBYJBly5aRkpJy1o85mrR77iHc1UXg8BH0oogkiPRceSU5JuOUY8663FwErRZFUdDr9QSCAQIuJ0paOh0dHeTl5Y0pfkpLS6Ouro5wOKyGPvr7+6mtrcXpdGKz2Vi8eLG6fqAoCr0PPYzkdqujcjkYpO+hhzE/8zeCgsCBAwdISkqioqICURSZP38+FRUVNDY20traSltbGwUFBZSVlTFr1ixmzJhBf38/nZ2dtLa20tzcTEJCAvn5+eTm5p6Wp4t361a67/9atOhrBEN5GTk/+OFZWfyciF7XsRlBTuL5kTZYkGrm8Q8v4qUDPXzzhcMMeqNf/+f2dbOteYjv3zqHVeXnaEb6H8pURuhLgEZFUZoBBEF4GrgRGC3oCmATot9mKzAETF5PfA7xbttO9/1fHbPoKZjNZHzh8yR/4ANndVocw+fzsWXLFiKRCMuWLSP5HHmBj0ZjtZL7ox8Ram1FdjqhsIjeXTvZtWsXq1atmlKs2bxoEYbSUoL19ej0eiLBELICzisup3PfPhoaGqisrCRnpKo1JuiDg4Po9XqOHj3K0NAQZrNZzeoZfQGQBgYId3YiWq2EQiGCwSBarRZNOIR7334OKTKKorBgwYIx4ROLxcLcuXMpLy+nsbGR9vZ22tvbycvLo7y8nKysLLKysgiFQnR1ddHZ2cnhw4c5cuQIGRkZ5Ofnk5mZedKQjOz10v/Ijxl+4okx9yfecjNZDzyAeBaaSE+GJ3Ds63U+pQsKgsD1c3NYWZbGt146ovrD9LoCfOjxHdxZU8hXr56BSR+3ETgbTEXQc4GOUbc7gePn6T8HXgC6ARtwm6IoMu8ist9P/yOPMPynP4+537J6Fdnf+Aa6s5giOBqv18uWLVuQJImamhoSp7Eg6VQRBGFMRsqiRYvYsmULe/bsYcmSJSe1FhC0WnJ/+AOczz2PZ8sWDMlJ1OXnY50zh8VFRdTW1rJnzx4aGhqoqqpSY+c7d+4EwGg0MmfOHPLz8ycUT8F4zFVSVhQkWUIKyYiBANv37yOYl4fBYMDr9WIwGMYt4JrNZubMmTNG2GMzh/LyciwWC8XFxRQXF+N2u9WQTF9fHzqdTg3JTGQU5tm4id5vfINwd7d6nyYpiaxvfoOEq66a0vs/ncTSFSHq13K+kWzR88ht87hyVhb3//MgQyOj9T9tbWNz4wA/ef/8ac2ZjxPlpFkugiDcClypKMonRm5/CFiiKMpnR21zC7AC+G+gFHgdmKsoiuu4fd0F3AVQUFCwsG1UVsB0Eqirp+sLXyA00jcSol++zK/dT8J1150zTxS3283WrVtRFIWamppz3khhKrS2tnLw4EEqKipOK2e8ra2NAwcOUF1dTVFREd3d3dTX1+PxeMZsN3PmTIqKik6ay97zrf/Du2EDitmEx+vFKMmQmUHjHXeAKCKKIrIsIwgCSUlJpKWlkZaWRnJy8rh9BwIB1YZAURRyc3MpLy8fk46pKAoDAwN0dHTQ09ODLMtYrVby8/PJy8tD6/PR/93v4nz+hTH7tq5ZQ/bDD02Lh8/p8NV/HOAvO6LjrIdvquYDSwvflfOYCnZ3kK/+4yBvHO1T79OKAl++qpJPrCyJ562fImea5dIJjLYQzCM6Eh/NR4HvKNGrQ6MgCC1AFbBj9EaKovwa+DVE0xandvpTR1EUHH97hr5vf1ut0gOwrFlN9rf+76SdgqYTl8vF1q1bEUWR5cuXT5jTfT5QVFSEw+Ggvr6epKSkU06hLCwspLe3lyNHjpCenk5ycjKJiYnjBN1oNE6pMCnzS/di1+txv/0W2mCIQGkpw9dfj9lqZfXq1YiiyNDQEIODgwwMDNDY2EhDQwOiKJKSkkJaWhqpqakkJSVhNBqZNWsWZWVlqrB3dnaSk5NDRUUFNpsNQRBIT08nPT2dcDhMT08PHR0dHD1yhM4/P0H6v/6FOOq1aBITybz/qyTccMO7apZm0B57L/2h8WZa5xPpNgO/uXMhf9vVwYMvHsEXkojICt9+uZatTYP84Na5pFrPvT/9xchURuhaoB64DOgCdgJ3KIpyeNQ2vwT6FEX5piAImcAeoiP0SU05pjsPXfb76fnaA7hefvnYuZvNZH71PpJuueWcfvmcTidbt25Fo9GwfPnyc9JB/kyQJInNmzfj8/lYtWrVKZ9vMBjktddeU2+LokhxcTElJSXU19cTm4mlpaVRWVk5pQVhORSi7uhRGtvbAVi5ciVJCQl41q/Hs34DmsQEEq+/Hk1xsSruAwMDuEaybzQaDampqeoIPiEhgXA4rAp7JBIhOzubioqKcTOnQF093Q8+SHDPnrHvU00NmV+7n/TzoGPUj9+o58dvRP34P3tpGfeuO/XZ1btB64CXz/91H/s7HOp9WQlGfvWhhcyNFyNNiTMaoSuKEhEE4TPAv4mmLf5OUZTDgiDcM/L4Y8D/AX8QBOEgIABfOZGYTzfh7m46PvMZgkeOqvcZysvI/fGPz1ql52QMDw+zfft2dDodNTU1J/Q8OV/QaDQsXryYDRs2sHPnTlauXDkuPh3u74/a+GZkjLk4hsPhMf7nAJdddpnqaT579mw6OzuRJAm3283mzZvJyMigsrLyhE0tRL0eeeQczGYzSUlJ9D30MJ4NG0YseGVc/36NrG9+g8yaGnVmEQqFxgj8kSPRtXudTqcK/JIlS7Db7bS0tNDT00NWVhbl5eXYNBoGfv5zhp54EqRjo14xK4vgh++kKyuLlqNHMbe1qSGZd+v/m2rRq3/3u4In2PL8oijNwjN31/D9f9fym43RkGivK8Ctv9rKw++p5tZFZ9ZP4D+dC75S1LdzJ52f/4LaPAAg6X3vI/P+r57TrAOINlzYvn07BoOBmpqad9WH+3Sw2+1s376d7OxsFi5cCEBkYIDe/3uIwNHoxVJfXEz2N/4XISOD5uZmtbw/Ly8Pv9/P0NAQy5cvHzMK37lzJy6Xi7Vr19La2kpjYyOhUIisrCwqKysnXFvwer2sX78eSZJIT09nfmoqHZ/+DKLFgiIICAIo/gCapCQKn3xi0hlYIBAYI/A+XzR/W6/Xk5iYiNvtJuD1krB9B2lvvIE42ptdqyXlQx8i/TOfRrRYkCRJDcnE8utTU1PJy8sjJydnyhYH08GbR/v4+B+j35+VZWk88YnpdXE8F7xV28cXnt6Ha1TGzgeWFvC/188cE1KKM5aL1j53+C9/offhb0Nk5AOh1ZL1wAMkv/+2aTjDU2NgYIAdO3ZgMpmoqamZcted843GxkaOHj3KrFmzKCkpofOznyNw9CjiiOjKbjeR5GTaPvoRQpEIWVlZVFVVYbPZiEQirF+/HoA1a9aoAhdbeL3sssswm81EIpExF4OcnBwqKyvVdQZZltm8eTNer5eMjAx6enpYLooM/vgnYDbjcbtBEDDodWgDQUpfehFxihdPn893TODtdsT9B0h96SUMfX1jttMtWED+tx6ctCG33++ns7OTjo4OvF4vGo2GrKws8vPzSUtLO+shvoY+N1c8sgGAvGQTm74yvc22zxWtA17u+vMu6vuOrVPML0jilx9YSNZ5kl9/vnHRlf4risLAzx9l4NFH1fs0qank/eTHmBdN+DrPKna7nR07dmCxWKipqXlXGhBPF2VlZTgcDo4cOYLZ78d35AiyyYQUCBAKhUCREbu7SfZ6qbj88jFhE61Wy/z589myZQuHDx9m7oilQKwidmBggIKCArRaLRUVFRQXF9PU1ERzczM9PT3k5eVRUVFBW1sbDoeDRYsWYTQa6erqwmU2g6Lg83iiDosaDUGPh6DeQF1zMyVlZVN6381mM2azmTSHk/6//AXf1m1jHg8nJTF47TV45s7l6NGjFPp8VFZWjtu3yWSivLyc8vJyhoeH6ejooLu7m66uLoxGI3l5eeTn55+1xfDCVAsaUUCSFbocfnyhCGb9hfd1Lkqz8M9PreArzx7gpQPRWpG97Q5ufHQTv/7Qonhc/RS54D4BiiTR++C3cPztb+p9xlmzyPv5z9BlZ5/z8+nr62PXrl1YrVZqamrQ6/Unf9IIciCAb8cOZI8H09y55yw3fiIURcHj8eBwONCOVIPu3bWLnEAAGSUat4bobwGkUIi+vj5kWSYpKUnNK09JSaGsrIyGhgYyMzPJysrCarViNBpVQY+h0+moqqqipKSExsZGWlpa6OiIpuJlZGSQnZ2NoigYjUb6LRaM+flw9ChmqxVRlpF1eny33EJTSwvNra1qheiJQl3B5hbsP/sp7ldeHXO/aDaTctcn0b/3vZicTmpra4FoWmZsUbeoqIj09HRSU1PHFGIlJyeTnJxMdXU1vb29dHR00NTURGNjI0lJSWpV6nQahem1IsVpFhr7PSgKHO1xs7Dw3BesTQcWg5af3T6feflJ/L9XapFkhT5XkFsf28rDN8Xj6qfCBRVyUSIRur96/xj/DMvKleT97KdTnnJPJz09PezevZvExESWLVt2Sl/YUGcXXf/9RSSnCxQZEEj9xMdJft/7zt4Jj6AoCj6fD4fDof44nU6kkYVArVaLRqMhGAiQ85vfoB0aRrBYomGkgB/JYGTwS/ficLtRFAWtVktaWpqa/mcymdi0aRN+v5+1a9diMBjYs2cPAwMDrFu3btLzcrlcasgGoLi4mPLy8qhfS0sLRCJU9PZiPnwETWICSTffgnnBfLxeL42NjXR2dqIoCnl5eZSVlY0ZHYfa2xl49Bc4X3wR5FE1bxoNSTffTNpnPo3uOBOxSCTCoUOH1IvMaBITE9UMmpSUlHHx82AwqBYuuVwuRFGcdqOwL/51n1qJ+Y3rZ/LRFeNtjC80NjUM8Kknd4+Jq39pXQWfvqTsvOup+25xUcTQY9atrpdeUu9LuOF6sh96CPEURsXTRVdXF3v37iU5OZmlS5ee8oJY171fwn/ggNpbU4lEkP1+Cv/we3Q5OdN6rn6/f5x4h8NhIJpimJiYSFJSkvrj9/upq6tjeHgY7cAgWX/7G7aRi47GaiX74YcxVlYQDocZGBjAbrdjt9vVBUez2YzRaGRoaIiUlBRWrFhBe3s7+/fvZ+3atdgmsZTduXMn/f39LFy4kL6+Pjo6OhBFUb3QAFx//fUnfJ3Nzc20tbUhSRLZ2dkUm80En3oK5z+fG5O5AmBbt470L3wBQ8mJhVCWZTo6Oqivrx/TCHuqRU6jjcJCodC0GYX9cUsr33ghmj187exsHv3AgjGPK4rCG0f7+OvODnwhiStmZnL7koITNo0+H2gb9HL3n3erXusQXSz95g2z0GnOv6rYc80FL+iKotD38LfHeGgk3f5+sr7+dYRptkSdCh0dHezbt4/U1FSWLFlyymKuSBJN665ETEggIkmqEIp+P47rrsW/YAGCIKgjktjfE90+/rFIby/a9evR9/URzC/AtXgR8iR55TabjYSEBHW06HQ61TxuiKYzSpIEsoxlYIDszEwoLkbUaic8tt/vx263jysqGk1SUhKlpaXjnt/W1kZvby+ZmZnq4z6fjwMHDowR9BUrVqDT6U74noRCIVo2bsT/xJPY9u1DkMe6UFhWrCD9858b5+F+MmRZprOzk4aGBnw+H1arVfWQHxwcxOFwoCgKoiiSnJysCnwsJCXLsmoUFgtXnYlR2JFuF9f8NNpQI82qZ8f9l4+punxqexuPrW9GKwpoRAF/WGJRYTKP3DbvvB/tugNh7nliN5sbjzUOX1GWyi/uWPgf79p4wQt6/yM/ZvBXv1JvJ73/NrK+8Y135UMZK3VPT09n8eLFp9W/UlEUmm+4EUEQUDQagqEQKAqKz0f4Ix8mPGsWsiwjSZL6e/Tfsd/H/++0Q0Nk/+EPCMEQikaDKElIVivdH/8Y8vmQDx+JoLfbo12GzpJBmaGjg+R33sFy8BDCce+Pr7SUoSvXESguPuUL5uj7FUVRm2WPxmazqe37JsNkMpGUlEQkEpmwZ6perycrKwuNRnPS85IVhdueasYRiF7wHrupiMoMM4IgIMnwib83IiCg0wiIGhGNqMETjPDYhxZSlXX+2VAcTzAi8eW/H+D5fccK00vSLTzx8aXkJF1YKcHTyQWd5eJ84YUxYm676qroyPxdEPOWlhYOHTpERkYGixcvPu04qCAIJN1yM8N/fgLRbMao1yN7PGgyMii8444JQ0jhcBin0zkmbBIb2UN0NJ2xZy86SUKTkqI2oMDjYcbQMNqamnEXBLfbTW9v75jjpKSkRAVh1HaBQGDMKPl0MLa2kvbP5xDDIZAVgnl52G9+L/Jxax8ajQZRFNWR9slIT09Hp9WiPXwE7YsvIh46NH6juXPpX70aV8GxxTWLxTLGtvf4n4nuH32fzWZTR9wxYmIey8EPBoPjhN/v9+P3j+3uM5pQKET7SHVs7P04/tijKbcJ7AxEvwvP7Wzh2sLo494wuH0CJg2ER7bV6/Ug6uhxBC4IQTdoNfz4tnmUplv50ev1ADTbvdz62FZ+95HFVGadnW5QFzLntaAH6uroeeDr6m3LmtXkfv9758Ty9niampo4cuQIWVlZLFy48IwXtVI++EFAwPH3v6O4XJgWLiTj859D1OuRJGmMeDscjjHCEKucLCoqIikpicTERLRaLe3/+CdhWwKKVovbExUXMRzGu28vYs0yrFYrVqsVnU6nTvs1Gg3FxcWUlZWdcFHX6XSyadMmdc0g1vloKjOJzro6zH//e9RF0WBAK2owd3WR8u9/M/Ce94w5Tmwfo1FDP8cjSfhfew3jO+vRdx9vLwTeGVUMXXYZwcJCBEGIdlgaCb94PB4aGxuj/4uRRU1RFNFoNOpFZaL7Rv/WaDSUlZUhiiJ9fX00NTUhyzJut5uqqioKCgrGfE6CwaBa4DQ4OKj+T/V6vepBoygKw8PDqlGYyWRSq1JjtQ2jhV0s6GHnX/YBUBew8JN1y0feGpm/dO/CHYhg0omEQmF8gQAhKUwC42cXp8qAJ8i/D/Uy4AmyuCiFZSWpZ8VkSxAEPndZOSXpFr74132EpWia5s2/3MLvPrKYJcXntq/A+c55G3KRAwFab3s/wbo6APRlpRQ9/TSad8HkqqGhgdraWnJycpg/f/60tjKTJAmX04nT5VLF2+PxqF9Yo9E4ZsEyMTFx0tTI/h/9COfLr6BNTCQciRAOh5GcTobWrMFdM3GP1JihldVqxWazYbFYJg0jdXZ2snfvXkpLS5k5c+aUX6P7rbdp/9a3CI1U/+m0OvQ6HYrXS/3nPktxWRk5OTnqAuvQqKrfjIwMMjIySE9Pp76+nu7ubi5duhTXP5/D8eSTKMcVBCGKaC9Zi/aWW6CoaMILjiRJdHZ2jjvPmIfN8dufKaIootPpxl0QgsHgmFnWVMjJyVFTIDUaDUFJYe2Pt+EPRy9Uf797GQuLorOsrU0DfO2fhwhLsTUEhaXpCmvS/WRkZDBnzpzTqmZu6HPzmb/sxR+KEJYUQhGZqiwb37tlDrnJZy+0t77ezqee2I13xIzMqBP51YcWsabiP6tpxgUZQ+/9v4cYfvJJIOqTXfyPf5w0G+FsUFdXR319PXl5ecybd2aLSbHR2+jRt8vlUsVbr9ePE+9TqTgNd3fT/l//hez1IWq1KJEI2vR0TA8/xI7Dh0++gxHMZjM2m00d0cf+1ul0HDp0iJaWFhYuXEjOFLNx3G+/Td//+w5eRUGO2eTLMmIwRPv/fIklNTVqKl84HOatt95SrQHcbrc6ktUODJK0eTOJu3YhjMo2AVD0epyLFhG8ch1lK1aMa54xGX19fTQ0NDA8PIzBYKC0tJTCwsIxC93HXxQm+z367+7u7jEXJoDc3FwURZl0NjNRXH4qPNkgsK0/+loXpyvcWRFdmBVFkcGgyN4BhaAkMDdTR2WKFqfTOeb5+fn5Y2YjJ5uV3P9iA4d7PIiiQI8riCzLKECmzcj/vaeay2acvR6iR7pd3Pm7HQx4ov41eo3IT2+fz1XVWWftmOcbF5yg+3bvpu2DH1LbxWV98xskv//95/L0ADh69CiNjY0UFBQwZ86cUxLz0YU6MQF3Op3qlF+n040R7qSkpDPyfpE8XgYeewz3q69G4/FpaSTceAOORYto6usjHA6rozuPx0NfX58qOAaDQRVvvV6P1+vF4/Hg8XjU841tZ7FY1OfNmjWLnJyck150JI+H1tvejxSJ4JUi6LU6JJcT74yZDLznRvX9yMrKUnPJV6xYQUpKCookMfja6ww++STS7t3jFjpJSMBy661kfeyjDEci1NXV4XQ6sVqtVFZWkp2dPaX/2+DgIA0NDdjtdnQ6ndoI41QKxSaiv7+f+vp6hoeHMRqNlJWVUVBQcNLFdEVRcDqd9Pf309/fz/Dw8KTb9oZNfG9XiLCsYDOI/O32ElDkE15wPB7PuDWKSUNbx/HgbgG9CPYAgIBGFJFkhRSrngSjjhc/u/Kspka2DHj54G+30+WIrkUIAvzvdRdHHv5UuKAEXZFlWt57M8GRSj3rmjXkPfbLc74IevjwYZqbmykqKqK6uvqkx/d6vScs1ImJduz3dFvqdv3Pl/Hv3YtotaIoCmGnA3dRMb03v5eMjAyqqqrGdUsKBoP09fXR29uL3W5HlmV0Oh0ZGRlkZWWRnp5OKBTC4/HgdrtVkT9eXLRa7ZgRfWxUbzabj6U07ttH7/89hH94mHAwSKCwgMwHHiC/shK73U5PT8+YMEiu2Uzyvn1Ir75KuGN8eIT8fLyXXEJfVSWyTodWq1XdFCVJoqurC7fbTUJCApWVleoC6MlwOBw0NDTQ29uLRqOhqKiIkpKSM/bmsdvt1NfXMzQ0hMFgoKysjMLCwilnScUEPhaDH50hIyvwjV0CjlD0vf7zxxaxqmJqo+Suri4OHTpEOBymrKyM8vLycYvix18Q7nr6ML2uEHZPCAWQ5KiGpJo1mPU6fv6BBczJP7ux7S6Hnw/8Zhuto5plf/26mXx85cUv6heUoLtef52uz34OiHZRL335X2e9pF/yePBu3YoSDGFatJDa/n7a2tooKSlh1qzxDaRPtVDHYrGc1QtSuKuLto9+DNFqRZIkfH4fsiShC0fI/M2vSZ9CJyJJkrDb7fT29tLX10coFEIURdLS0sjKyiIzM3OMqHV3d7N7924g2uTC6/XidrsJjmosIooiFotFFXuLyUSopYW69nYiSUlceeWV6gjYbrezbcsWbE3NpO/Zg7Bnz7j8cQQB7ZLFtM+ezZwPf1htSjE4OKgWN8XCFkajMbqGMHJRTUpKoqqqivT0qcVb3W43jY2NdHV1IQgC+fn5lJWVnbFd7uDgIPX19QwMDKDX6yktLaWoqOiUaxlkWcbhcDAwMEBnZye/P+BjY2/0M7YiU+H9ZYoq0CfbdygU4siRI3R0dGCxWJgzZ47qvzMRbx7t43+fP0y3w89I/g0AWgEsOvifeQKl2SlkZmaSmZk5ZvASlmReP9LLm0f7STLruXlBHjNzTi/jZtAT5JN/2sWedod63/9772xuX1Iw+ZMuAi4oQW/94Afx74oKRcrHPkbml//nrJ5HoK6O7i9/GdkfAEUhGA7Td/llZN98MzNmzCAQCIzLOIlNVUVRxGazjRFvq9U6rYumU3oNtbV0ffG/Ec1m/H4/wVBUVDWhMMpXvkzavHmkpaVNuXBFURSGhoZUcY+JZFJSktpw2WazqS6KlZWVVFREu92Hw2F1JD96VD9ZfLi0tBRxYIDBvz1Dwq5daB2OcdtIZjOuxYtxr1hOYkUFdrudvLw85s+fP25bn8+nivvAwIB6oR3NsmXLpizsPp+PxsZGOjo61DZ2ZWVlk1a7TpWhoSHq6+ux2+3o9XpKSkooLi4+bQve7c0D3Pbr7QCYNAoPLVYY3Yc5NzdXzYqa7PNpt9s5cOAAPp+PgoICZs6cOWnm0193tnP/Pw8hyQoaAQQUZAWsOnj8hkz1Ag9gtVrJzMwkIyODH67vZnPTABpBRFJkREHgwRtmsaby9LqJ+UIR7nx8B7vajs0av3JVFfesKTnvi6dOlwtG0EOdXTRdfnn0hkZD2dtvjfPXmE4URaH9ox8j3NODxmrF5/MRCvjRSDLBB7+JIxJRS70FQcBms40ZfY+usnw3kUMhWm+5FUWSEI1GZFkm5PUiCQKdn/ss4ZGRbkJCguq3kpKSMuXpfixfvbe3F8eI4FosFjIzM+nt7cXn87FkyZITtq+LWeI6HA5SU1MZ6u3FcvAgCTt2YG5qnvA5wuzZWN9zIynXXot/xIu8p6dH/Z+kp6eTk5NDVlbWhLFuRVFwOBzY7Xb6+vrUc4+Rm5tLRUXFlBwRY/1JY7YCsaYYJ2rSMRWGh4dpaGhQG1XHhP1UjbxkWeGSH75D20gI4js3zaLK4JjQhyYpKYns7GzS0tJITEwcI3ySJFFfX09TUxN6vZ7q6uoJF79fPtjDd145SjAiM+QJAgomDWhFeHCJhmVLFmGxWOjr66Ovr4/BwUE63DKP14rYjBp0Oj06nY5AWCbRrOPv99SctgC7AmHu+M02DnUdq3K+a3UJX7266qIU9QtG0If+9Gf6vv1tACyrVlHwm1+f1XOIDA3Retv7EW02gqEQgUB0kUUMBPHccQemZUvHLFyeTlXoucKzbRt9D34LJRIBAQStjuxvPYhp4UKcTqc6ah0eHkaW5TE9ONPT08d9sScjEAiocfeBgYExi6YnagTd3d3N7p07KQsESTx4APdrryO7XOO2UxISCNQsY3DBAgKjmmTodDo1Pn+8SAmCQGpqqiruk81EwuEw/f397DmutRxAQUEB6enppKWlnXAhNBQK0dLSQktLC+FwmPT0dMrLy0lNTZ30OVPB6XRSX19Pb28vWq1WbeF3Kouyv3ynie++Gl17mpmdwL8+txJBEAgGg7S0tNDQ0DDuORqNRnWQTEtLU/usOp1O9u/fj9PpJDMzk9mzZ49ZtN/X4eALT+/FatDidrsRRRF/KILFqOOBpQa8Hg8zZ86kpKQEiBqd/W1rAz9d34ZBiC7YarVazGYzzkCEN/57zRktpLoCYe760y62NR/LLPryVZV8au273y5wurlgBL3j05/B8+abAGQ+8AApH/zAWT0H2e+n+T03IZpMyEQ/dKIoIgaD5P7wh5iqx8fPz2ciQ0P4tm8HQcC8dCnaCcrrJUliaGhIFfiYd4tOp1PFPS0tbUqLtpFIhP7+flpaWsak6MVsczMzM9Hr9Tj37aP2N7/Btmcv4gQhFUUQsK5cQdItt2K7ZC3CiIgFAoFxoZvj4/QTodPpKC8vJzc3d9LFzEgkwpEjR1Rr3NEkJSWpM5nk5OQJZ2GRSITW1laam5sJBoOqbfCpNtk+HpfLRX19PT09PWi1WnVRdirhsiFviOXfeZPASE767z+6mEuOC2W4XC46OjrGtQ2MEStyirlI9vf3U1dXhyAIzJgxg8KRIi1ZVvjUk3s43O1ECgVQBAFZgRsKFb5022UcPHhQ9bifO3cuoihyuNvJp5/cg9WgIRyJ4PP5UAQtaUk2/vmp5Wc8mg6EJT73l728duRYbcID187gE6tKJn3OlVdeybp167j33nvP6NjnkgtG0BvWXkJkpBS95MUXMJSXn/Xz6P/xj3G9+BKixQKiiOx2oy8uJv/Xv3pXjL/ONaOrF+12u1qWbjabxwj8yUaKfX197NixA4j6lUgtLVgPHCDx4CE0x9kLxAinpuJavBjXooVULFtG2RSbL4/OVS8tLVWFfrIiHUEQKC8vJyMjQ82nH72vlpYWGhsb1QVUg8FAKBRCURQ0Go0qcLHnj0aSJDo6OmhsbMTv95OQkEB5efmU0yUnw+1209DQQFdXFxqNhsLCQsqm0MTjmy8c5g9bWgGYnZvIC59ZMeF5TGQUFiNmJAag0Ruoc+k40ushwwRrSpOoWTQPm82GPyTx9M52/r6ljkSjhjtXlhHuPMScOXMoKCigoaGBuro6kpKSWLx4MQaDgf/5+wG2twyiE0UCoRDhUIi7FqfwyWuXn/Z7NZpgROLDv9sxZqT+pXUVfObSibUkIyODp556istjod7T4Ctf+QovvfQSHR0dWK1Wrr32Wr773e9OqRn66XBBCLoiy9TOnqNanFbu34d4Djr/KKEQg3/4I84XXkAJBbGsXEX6Zz6N9iz9M853vF6vOnofHBxUFxUTExNVcZ8s/l6/fj19zz9P+tFalJaWiQ+QmIjhssvoKCokXFzM6jVrOHjwIL29vaxcuXJcauVkxBwvV69erT5HkiR1Ma63t5fuCewAYoyujrVarRgMBjo7O2lpaUGWZbKyskhNTcXj8YzLnomN3tPT09ULnSzLdHV10djYiMfjwWKxUFZWRl5e3hmts3g8HlXYBUFQhX2ymUevM8Ca779NMBIV5L/etYylJScOB4VCIbq7u+no6MDhcCAIQrRiWG/kp7s8tA5HL26CABYt3DVDITvZwrJlywih5dHnN9Hnlbl++Wy09jqQJdauXYsgCPT29rJ37161EbnFlsjLB3t442gfyWYdS9IVGGqjrKyMGTNmnPb7NBp3IMzH/rCTna3HFkq/enUVd68Z2zC+q6uLvLw87Hb7CbN6Tsb999/PrbfeSnV1NQ6HgzvvvBO9Xs/zzz9/2vs8EReEoMvBIHVz50VvaDRUHTp4zhc0oh/ai28R5XSJLSrGRu8Txd+TBgaQt2zB89bbhJqaJt6PwYBm+XJcc2bTm50NIxeDxMREKisrSUxMZOPGjeh0OlavXj0lAQyFQrz22muUlZVRVVV1wm29Xq+6oHkidDodOp1uzEg/PT2duXPnRtseDgzQ398/JnsmMTGRjIwM9UIXE7GGhgacTicmk4nS0tIpFROd7DU0NDTQ2dmJIAgn7M50/z8P8tT2qMHXjfNy+Mn7x2cDTYbH46Gjo4POzk42dwR4qV0k0aRFq9VFfWa8QRamSlxbCMNB+O1RAV8kmrhoMZspTtHznoxhli87tkjudrt5dcN2GgaCzKkq5apFFWhH+ZofPHiQ1tZWtY/tdOALRfjkn3aNsd/96e3zuWHusQXel156iXvuuUetf3juuef44he/yHe/+13edwaNZv71r39xxx13jKvInS4uCEEHqFuwEHnky1S+cQPaKaaWxTk3RCIRBvv7Gdi6lcBbb6PbvRvdBDFxAFmrxTdzJqFFi1j8X/egH0nzi315IVqQFIlEVKvYSCQyaTriRGzbtg2fz8ell069QbLf71ezZUbH/WPFUaIo4vF4JozTxwy0LBZLNN9/pOH08PCwGp5JTU1VR+9+v5+GhgaGhobU1MSioqIzakU3Oo0SIC8vj/Ly8jH58Ye6nFz3s01AtDR+032XkGE7tcIoRVH46jO7eaduAL0gqQ6TYRnSbQYeXJvKd16tZd+AgFkbew4EFZH3FErMSYWrr74arVbLE9va+PWGJkLBYLTxSKKRxz+xknTbMbOx3bt309PTw/z588nLyzvt92c0gbDEnY/vYEdr9P8sClFRv25OVNS/9a1vsWPHDp599lm+9KUvsWXLFv76179SVlbGpz71KZ566qlJ933fffdx3333TfjYf//3f7Nz5042btw4La/jeC4YQe++76u4XnkFJRgk7xePYjuFL2qcs4fkdOLZsAH3W2/h3bIVeZKRh6zT4Ssvxzt7Np7qWShGI4sXL1arNHt7e9m5cycFBQXMnTsXWZYZGBhQ891HdwOaNWsWWVlZJyzkiXnTr1mz5rQ6/wSDQVXcBwcHURQFs9mspvRptVr6+/snzA6BaGw+1p0pGAwSCASIRI61TouFZzQaDcPDwzidTnQ6nbrQeSa2An6/n8bGRtrb29W2e+Xl5epi9k2/2MzekYKbu9eU8NWrTz2c8cctLfx2YwtiJIAgCFitVtzBCCvK0vj2TbO5+ZebGXB4kSPHLAQCssjcFIkbi0Zu6xJ4ZHcAqzFqyBYIBBj2Blmca+SnH1k9JmS1fft2BgcHWbJkCRnTlK487A1xy2NbaLJHQ2YaUeDnt8/n6tnZ3HTTTWi1WhobG6mpqeGRRx454wbvzz77LB/5yEdYv349CxYsOPkTToMLRtD7f/hDBn/zWwASrruO3B98/904tf94FFkmWFeHZ8NGPBs24N+7d2wfzlGICQnYLrkE25XrUObMYdDt5tBxnuSJiYmYzWZ6enqw2WwThlVixUxbtmwZc39CQoKaNXN8zncwGOS1116joqKCyilUw56IUChEb28vPT09ajqm0WgkJyeH7Oxs9Ho99fX1dHVFe3haLBYyMjLUTByv1ztmYfFkCIJAcXExpaWlZ2QrEMuPb21tVQufysvL2dzm4e4/Rwv0LHoNb//P2lMepdvdAW752Tt4ghLJNjNhOSqIj35gATOyE/jvv+5jd/swYiRAJBJd+/JH4LY5yVTpo6PiA24Lz9b5sGijr1mv16MIAh5fkG+vjF7wYxfjSCTC5s2b8Xq91NTUkDxNTVDs7iC3/2Ybjf3RTlp6rcjf7q7hPavmEQwGqays5J133jnj4zzzzDPcfffdPPvss1xyySVnvL/JuGAE3X/wIK23jsSuRDGa6VJaOsGz40w3Ebsd75YteNavx7tlK9IkoRQAbWYm1tWrsa1bh2XpEjXNEKKx3vXr16sLlcc7DgLq4mp6ejoJCQlj1i2cTicbN24kMTGR3Nxcent7GRoaQlEUjEajWqmampqKKIps3ryZcDjM2rVrp+29CIfD9PX10dPTQ39/P7IsYzAYyM7Oxmq1qrMKnU5HaWkpxcXFaDQafD6fmmI5+vfJDK8MBgPLly+fUoHTZASDQVXYJUkiOyeH+95x0WiPhjDfvzif79x8ai33Ghsb2bzvKIdD6TQ5ZUrSLHx0RbFaqn+oy8ln/7KHYDCEEgkTksGshU/PUrCMRJUi2XN45K0mzDqBcDhMKBhCEgRsRiNfmisTDoeZP38+2SP2HsFgkE2bNhGJRFixYsUZvSej6XcHeP+vt9E8MlJP0oTY/+33cvDgQW6//XZuvPFGHnroIXX7e+65hydGtbw8nvvvv5/7779fvf373/+ee++9lxdffJEVK1ZMyzlPxgUj6ADtn7wL70jsyVBRQeFTT6GxTq+RVRyIDA/j37sX3/bteLdsJThJWAEAQcBYXY3tssuwXnIJhorySVPhNm3ahM/nY+3atRgMBrZv364aScUWvOx2u1oWHst7jom82WxW/ecXLFhAbm4uoVBojImYJElotVoyMjLw+/0MDw9z6aWXTrvhGURHjaPFXZIk9Ho9JpNJtT7W6/WUl5dParYVCARUgXc4HBN6sY+muLiYlJQU1Z/+VLJkQqGQKuwH7BF+eST6XEGA5z+9gjl5SVPaz/DwMJs3byYrK4tFiybUDiAq6j959SCNvQ4um1PEstQg/sFj2UW+CPyuNYFgRMZq0BIMhRn2+rksR+Z9C3Jwu9243W4qKiqoqKhAEAS8Xi+bN29GFEVWrlw5pRmMOxDmzaP99LsDzMtPZlFh8riGG+2DPm54dBMOXxh/236GXvgeTe3diP4hampq+NrXvsY999wzpfdnND/96U958MEHefXVV1m8ePEpP/9UuaAEPXDkCC23vR9GsggsK1aQ97OfIp4PPTEvYCSnE9/u3fh27MS7ZQvB+voTbq9JScGybBnWNauxrFyJdgqVkEeOHKGpqWlM3DzmJw+MMeMKBAJq9szAwIAaP7dYLKSmpqpt2C6//PIxmRySJKkj5N7e3jEWsLNnzyYzM/OMbIhPhCRJ9Pf309PTQ19f35h4OUSzZCbqVDQZscrV2traEza6MJlMJCYmjvOnP5HvSygUorm5mXtfaOLwyCTp0ooUfvexmimd14YNG1AUhTVr1px0ETe2lnHFFVdgNBrVRigxenzw5oCVLo+CWa/h1gW5LEn20d7WhtFoRKPR4PV6ycrKYv78+Wi1Uc/2LVu2YDKZ1Mbgk9Hl8HPPn3fjCoSJSApaUWB1ZToPXj9rnKjvbB3iw7/bQe+mZ/C37GXhPT/gn59aQW9rPWvXruXxxx/nPcd10ToZgiCg1WrHxd9P1Cz9TLigBB3A8c/n6PnqV9XbhspK8n76E/SFhefq9C5oFFkm2NhI4OBB/AcP4t+3P9r56UT/a50O87x5WFaswHrJWgwjo6WpYrfb2bZtG4WFhcyZE53ax2LisUXDlJQUli1bNuF+Y/nesfz30WJZVlam+s+MFspYu7bNmzeP2VdiYqIamjmdxdKpIMuyavvb29s7zgSsurqaoqKiU3oPQ6EQ+/bto+/4LkyTYDQaxwh87O/RwnKkc5hrfh5dlxBQ+MFlyVy6aOYJ49OxjJMVK1ZMKY7d1dXFnj17uOSSS9QQSezzMOb1SbCiZilZmdEFz4GBAfbv34/P51Obb9tstmi+usXCwMAA27dvJykpiWXLlk2a9vm1fx5kQ72dJLMeUFAUcAcj/ODWuSwuGl9PsrVpkI/8foeaq7+sJIUnPr50TCrl+cwZC7ogCFcBPwE0wG8VRfnOBNusBX4M6IABRVHWnGifJ+tYZP/Zzxl49FH1tmixkPbpT5PywQ+MidnGiTa3CBw5jH/vPvx79uDbt2/STBQVnQ7TrFmYFizAUrMM88KFpz0LCoVCvPPOO2oeuUajIRAIsGHDBrRaLatWraK3t5d9+/ZNqX1dzBp2z549YxoqazQaUlJS1LTAmO9IY2MjR48eZdmyZTidTnp7e1XPdrPZrNoQpKamnpU6A1mWGRwcpLu7e0yDZ4jaIMyfP/+UUhUjkQjt7e3U1dWNmwWMRhAE1SZ49HYx35uY0P/vGz1sbol+HhZnCNxZLpGenk5FRcW4asbYaHvGjBlTrtyNVQmXz13Cc4eH2N/ppDTdwi1z0uk4vHPc9unp6cyZMwez2UwkEuHo0aNqKmvs/BctWkRaWppq0xwL/Yz+//l8PpxOJ3f86RCyFEEZWZS2WCy4gjIfXVE0adn/W7V9fPyPu9Qxzs0L8vj+LXPOSl/U6eaMBF0QBA1QD1wBdAI7gdsVRTkyapskYAtwlaIo7YIgZCiK0j/R/mKcTNABhp95hr5v/R/KqNGPrqCAtLvvIuH66xH/A4U9MjCA/9AhgnX1BOtqCRw+QugkBTMAaDQYZ87EvHgx5iWLsSxZMm1hrB07dmC321m1ahUJCQnIsszWrVtxOp2sXLlSHSXHctBPpX3dzp076e3tpaqqSrUpGB1/T09PV+PuowtTJmveEcuYSU9PP22r2hOhKAqDg4Ps27dvzMUIYN68eWRlZU1Z3GVZVm0FfD4fWq2WpKQkdDrdmOKmhIQEbDYbJpNpXMepYDBIqxt+eODY6PPeOTJFI+6/er2eBQsWkJ6ejtvtZsOGDaSmprJ06dIpX/wGBwf59ztb+GObFV9YwajTEIzI6DQCD6zNwNMdLTgbbSkAUFFRQVlZGRqNZsxoHaIXq5i5V0tLC4cOHUKn01FYWKjaWcde/y8OCzhCImajFiky0pBDY+C+a2Zy/dzJP2c/e7OBH75+LPT44ZpCvnnDrPO+uPBMBb0G+KaiKFeO3P4qgKIo/2/UNp8CchRFeWCqJzUVQYdo5kv3V+4jdJyZkCYlhcSb3kPCVVdhnEJHoQsNJRQi1N5OsLGJQF0twaO1BOrqiPT0TOn5mpQUTAvmY6qejXF2Naa5887K4nLME320mMa6PcUWNWPIssyWLVtwuVysWrVqSp7iwWCQ9evXYzAYWLVqFaIojom/2+32MUVARUVFqntgTDgjkciY5h3hcPiEzTumC0VROHLkyDgjrJPZ/k60n+7ubhoaGnC73VgsFkpLS7HZbGpzj1gVb2wWE2usbTAY8Hg8fP6Zw2xsiRqxFScIfL5aQjPJV2bmzJmkpqaeNE4fw+l08q2nN7B10EiqzaQWi7kCYS6pzGCZvkOdQWRmZjI0NDSmIcyCBQvIzs6ecLQO41vjJSYmjukAtq83yP++cASdRkAnCgy6vJg0Mr/9wBwqiidvdiHLCvf/8yBP7zzm3nn/NVXctfr8zqw7U0G/hejI+xMjtz8ELFUU5TOjtvkx0VDLLMAG/ERRlD9NsK+7gLsACgoKFp6sFDuGEg4z9OSTDPzysQlDCdqsLGyXXoplxXJMCxZM6DJ4PqIoCpG+PkKtbQSbGgm1thFubyfU2kqos1P1tTkpGg2GsjKMs6sxz1+AacF89KcYvz0dXC4XGzduJC0tjaVLlwLH4qnFxcVUV1ePe87xoZipjFZjU/rJ/D7cbjdbtmwZs0AqCILqPzPaNXF0846YlztAcnKyOno/0+YVx6MoCl1dXWMWCmPnmJqaSnZ2NllZWSe9qCiKoja1djgcGI1Gtal1bGYQu8jFFuSMRiNpaWmEDIl88MlawlL0+/7lKyv4wIIMnE4nR44cmbARSOz5EzUMHx2n93q93PObt2nxG0i0GHA5XdHqX62BgjQLX1+dytGjR0lISMDlcpGfn49b0vL6vhZ0IlQkRdMdjUYjkiRNeC4VFRX09vbicrmorq6muHhsq7n1df38fksr/a4giwqTWGxzgm+YuXPnUlAwuahLssLnn97LSweiAyWtKPDXu2tYWHj+asiZCvqtwJXHCfoSRVE+O2qbnwOLgMsAE7AVuFZRlElTKaY6Qh+N5PEw/NRfGP7LX044UtWXlWKaMxdDeXn0p6wUbXo6wjn2M5eDQSJ2O5H+fiL9dsK9PYS7ugl3dxPu7CTU3o5y3LT8ZAgGA8aZMzHOmIGhogLjrJkYyssRz8II80RIksTGjRsJhUKsWbMGg8GA2+1Wc8hramomzfSILZZmZmaOi4tOxoEDB2hra2P58uUTeo+73W7eeecdZs2aRVJSkpo9c3xZfixFMhZ/d7lcamhmdPOO2KJqcnLytF0YZVmms7OT+vr6ceEYgJSUFLKzs8nOzj5ppo7dbqehoYHBwcEJbQX8fv+Yzk2hUIh/d8BL7dH/iVknsuHLl5BmM6oXYZ1Op65/QPQiZ7FY1PDN8XH6mLjrdDqe2trEGz160hPNeL0+wuEwvghcVZnM19+7kLfeeovExMTojKvFy786BBQlmk6pE+HOCoW8kQnkzJkzSUpKorOzU12TMBqNLFy4kKamJnp7e08atpMkiV27dtHf339Sj5hQROa2X29VK2vTrAaeuaeG4rTzM136XIRc7gOMiqJ8c+T248CriqI8M9l+T0fQYyiRCJ6NG3H/+zU8b7+NNBUTHK0WXWYmuuxsNGlpaFNT0aSmoElKQmOxIFosiGYzgk4HWi2CVhf9tCkyKApKJIISDCIHgyiBALLXi+R2I3u8SE4n0vBw9LfTiTQ0hDQ0pPrSnC663Fz0JSUYKsoxVlVhqKzEUFwcPcd3mVg8PNbOLRwOs3HjRiKRCKtXrz7paDMWF62qqqJ8CjbJkUiEDRs2IMvypKl0b7/9NkajkZqaY6l5sZ6jsRBNbORqMBjG5L+bTCYCgYA6ch8cHESWZfR6/Zi4+3Q0OZFlmfb2dhoaGggEAhiNRhITE/H7/ao/fayrUHZ29gnz64eGhmhsbKSvr29S/3RFUXC5XPT09vGJvzfR6Y7Gsa/Ig9tnmvB6vWi1Wq688koAOjs7aWhowOfzkZCQQEVFBVlZWZP604dCIYIS/LZWYCAgqMe06ODuGQoJo6JK3jD88ICATgCNGBVqbyiCVYzwmVkK7jA4Ijpq5lZQVpjH4dZe+tvqEcN+RFGkurqajo4OjnY7yS2tYkllHmnWicv1ZVlmz5499PT0jGmTOBEdQz6u/elGXIHoRSs3ycQz99SQk3R2UmDPhDMVdC3RRdHLgC6ii6J3KIpyeNQ2M4CfA1cCemAH8H5FUQ6N32OUMxH00SiRCL5du/Fs3IB/1278hw/DCTIDzjfExET0RYUYiorRl5aiLyhAX1iAvqgI8SzlU58pMU+W0U20d+7cSV9fH8uXL5+yD/TevXvp7Oxk6dKlU/LuiKUo5ubmTmjgVVtbS2NjI+vWrZs0Nh0IBMaMXGPxd6vVqgp8bAYQi7v39/cTDofV7j6jm3ecCZIk0dbWRkNDA6FQiMzMTPLz8/F6vfT09KgzhoSEBNWCYLLKSZfLRUNDA93d3YiiSGFhIaWlpeNG+i8f7OFTT0Y7Nuk1AvfPk0gdufYaDAY1RJWWlqbOArxeLzabjYqKinE+7zFHzk2bNhGSYP8gtHsEsswK81JRK0Zj1DrgmWYRg3hMd6w2K+6AxOJMkY1tXjQCeMKAIJBmNSKIIqvzdCyyDBOS4Z9dFhoGoz2AjUYjd68p4wPLJk5pVhSF/fv309HRcdIMq52tQ3zo8e1qg5C5+Uk8c3cNeu35lc44HWmL1xBNSdQAv1MU5WFBEO4BUBTlsZFt/gf4KCATTW388Yn2OV2Cfjyyz4f/4CGCdXUEGxoINjYSam1FGh4++ZOnG60WbVoa2owMtOnp6DIz0OXkHPspLLxg4v0xAoEA69evx2g0qouUsbTBU7U/lSSJTZs24ff7WbVq1ZQqPWOFSosWLVLLxWM4nU42bNhw0rjpaFwulzp6HxwcRJIkBEFQuxalpaWRlJQ0Ju4eCETNqpKTk9XQzJlUqUYiEVpaWmhqaiIcDpOdnU1lZSUajUb1l4lZKNhsNnXkPlGOvdfrpbGxkc7OTtW0q6ysTL0QyLLC9T/fxOHu6ExgdorCj2+egV6vVy90sbWIhIQE0tLS1FlOLF88IyMDi8WC0+nE6XSOCcXEFipNJhO1tbUYDAZSU1NVb/p2D/yhTsAgRifAALIC3oiAQafBpFXwBSUGA4AAiTpIsugIo+OLa/PZcKCJDZ0RbProeogMCDojj39kCRWZE699KIrC4cOHaWlpobCwkNmzZ08aRltfb+fjf9hJRI7q4idWFvPAdSdOsz3XXHCFRWcD2e8n3NtLpKeHyOAg0tAQkYFBJJcL2euN/vh80dBKJAzhkQ+pKIIoIogigtGIYNAj6g2IVuvIjwWNLQFNSnI0fJOQiDYlGU1yMuJxPiUXOoqisG3bNoaHh1m9ejVWqxW73c727dvJzs5m4cKFp7xPn8/Hhg0bMJlMrFy58qQhjVizaa/Xy9q1a8eFdt58802sVqu6SHsqyLLM8PCwOnp3OBzjbHHT0tJQFEUV91iIxGq1kp2dTWZmJklJSaf1fw+HwzQ3N9Pc3KxaCVdUVGCxWAgEAmNsfxVFwWKxqOJ+vHGZ3++nqamJ9vZ2JEkiJyeHsrIyEhMT2d02zM2/PGaC9r1b5vC+RfnAsfBMf38/ra2tYxwwJ6KgoIDk5GT2799Pfn4+8+bNw+12Y7fbOXxYncSrKYuyAr86ImAPabDqNQSCIUJytHGGLwImbdRn3T+SD6AXIc0IQUmgOjeBfj+4vf6ow6MAAgLeCNy1toK71544dFdbW0tDQwO5ubnMmzdv0jWe32xo5uGXj6q3f/S+ubx3wfRY+k4HcUGPMy3ERuKxEbDf72fDhg0YDAZWrlx52nnd/f39bN++ndzc3ClZjno8HjZs2KBWno7myJEjtLS0sG7dujPyHYdj8ffYyDXWtWh0/N1iseBwONTO9jETsVjcPS0t7ZQ7FsX8WGLdk/Lz86moqFDDJ8FgcIwzpKIomEwmVdxHL+TGGkS3tLQQiUTIyMggPz+f+/6+jw3d0e++Safh+bsXYRGCOByOCUfeJ6KkpERNzTQajepFYPTf11xzDRqNhnfeeYfeITevdAgcdQhoBYUlGQqOoMChYTBpxgq6QYQMiwZfSCbPquCRNEREAwatiC8Ywi+BxxfmmlI937nzkpN+BmM+QVlZWSxcuHDC/42iKHz8j7t4qzZaSqMRBX73kcWsqTg/+jPEBT3OGROLk8Yq9mIjZY/Hw+rVq8/YGCv2RZsoJW0iYvnvs2fPpqioSL1/eHiYTZs2TWujhBixzJFYiCYWmrBaraSnp5OYmKh6zdjtdiKRCFqtdkzc/VQuMsFgkMbGRjUvu6CggPLy8jGzknA4rIp7rIDKaDSq4h7rohQOh2ltbVXDOkEJvr9foM8fFf7F6Qp3ViiIoqjmecdyva1Wqyp8gUCAvr4+6uvrJxy9i6LIzJkzycrKwmQyqf+nmL9PLKPmeFrd8Mf66GKpqNXR742gIJBmEtAi448o3Fys4InAvzsEEERc4aibhYJCsl7h7vk2Pn79qpNeQGPnlJ6ezuLFiyecFTp9YW779VZqe6NFbElmHf/63Cpyz4NF0rigxzkjJsowiaURjjbiOhMURWHXrl309fVRU1MzYWri8cQaIsTCP7H9vPHGG2pj4rOFoihqaMFutzM0NKTG35OTk1UhDQQC9Pf3EwwG1bzz2Oj9RM07RhPrfNTe3o4gCBQVFU3YMDpm9tXd3a3a/kK0G5PJZEKr1aqWCACNTvjJoRE3RuAvH53HkvLsCQVxdK/ZgYEBdfQeK/oZCkIgApmmaPZKSkqKanO8Z88evF4vJSUldHd3qxeC1NRUtSAKwJFQym+3dOIPhggqIhqtliSTDkmGKysTuDwn2tj6pRaZ17sEUAR0WoGcJDOyFEEKBfjBVZksmUIqbEdHB/v37yc5OZklS5ZMeKHtdwW44eeb6XVFz3dBQRJ/u7vmXfd8iQt6nDNi3759dHR0qDngsQbN09nYFzjl1MfYAq3ZbGbFihWqEB06dIi2tjauvPLKs1LePxGj4+92ux2n04miKGi1WlJSUtDpdEiSpKb7QXTRMbaoOpXm2D6fj/r6ejo7OxFFkZKSEkpLS9Hpov0+vV6vGjIZHBw8YU/LFStW4Ha7+dwzR9g/EBXUsjQj//r8Ggw6rRpu6u/vx263qwVYZrOZ9PR0MjIySE1NxReBrz93iB2NvciShEEDt5YqFE9SmxWbxXg8HrXT1IYNG9RzzcrJw48ee0cztsQkskpnUpBmGzHeii6k/3tvE//3ajOiIiESDXEZjUaGPX4+WR5gfkUhQkohzQMe8pJNzM8fb6UL0NPTw549e7DZbCxbtmzCrKWdrUO8/9fbkEYWST97aRn3rjuzZipnSlzQ45w2sSlyeXk5VVVVOJ1ONm3adELnxDMhVpyUkJDA8uXLTzp97unpYdeuXer5QdRbZMuWLafkGTPdhMNhBgYG1PBLLP5uNBoxmUxIkqT2JY3FwGMj91jzjslwu93s2bNHXZA9Ho1GQ0JCgppxYrPZcLlc7N+/X91Gq9WSmZlJxJjMB588SmikgnRdnsL1hcqY7WLrBbE1g9F844VDvHW0H50SRpIlBK2BUCjIF2ZJmCa5lsZi6waDgcsuu0xtiqLRaNRmIhkZGXR2dmIymVi6dOmY4zbbPXzsDzux6DVIsoRG1KAg4A9LfOeyVH7yThutPj1anRZREJiZncAP3zcPk358aKW/v5+dO3diNpupqamZcBDx87ca+MFrx2okH7ltLjfNf/cWSU8k6Odm+BLngsTn83HgwAGSk5OprKwkFAqxa9cuDAYDCxcuPCsZPDabjfnz57Nr1y4OHz7M7NmzT7h9dnY2+fn5NDY2kpmZqYY7DAYDPT0975qg63Q6NY4N0fdytP97LP6u1+uRZZlAIEBrayutra3odDoyMjLUYqZAIKAuVMZG4BN1QbLZbMybN4/ExMRx/5uGhgZEUaSmpoZIJEJPT89IFWYX1+TDc63RC8hrnQIlNoVZKSd3igxFZN6utZNg1BEIRpDDMkh+wpLAkCaZ5YXRUv/h41KGrVYrgUCAYDDIyy+/TPpIM3hJkli1ahX79u2jvb0di8WCz+dj06ZNLFmyRLXyLUm3srAwmR0tQ5j0GiRJIRCOcNvifIZ1Vpo8XeiUECajFp1Oy4EuJ//c28kdS8fnqmdkZLBs2TJ27NjB5s2bqampGRcKu2dNKdtbhtjYMADAfc8epDzDRnXuyWdV55r4CD3OhCiKohpprVmzBpPJpMasV6xYMS5Nbro5evQojY2NzJs3j/z8/BNuG4lEWL9+PYIgsHr1arRaLQcOHKCzs5Mrr7xyWqo7p5Pj4++xqtSpkJKSMm7B0uVyUVtbS39/PwaDgbKyMoqKitRRfmwRsLKykqSkJDWMEgv9yAr88ohArSN6EUg1wv3zZPSa6IWpuLiY4uLicSGJUETm8h+tx2bQIssSoVAIrU5LUBL43+tncmlVJhBdqG5oaFB93ouKisjNzVV97G02m+qgCZCTk8Pw8PA46+QFCxao6zX+kMSft7Xy6qFe9FqR66rTuazEykOvNrK704NOiYAACbYEAhGZknQrv7lz8s5LDoeDbdu2odFoqKmpGVfA5Q6EuekXW9S+pMVpFl75/CqMunP/2YqHXOKcMrECnphjYuz2nDlzKDwHjUZiOe9DQ0OsXLnypDHmoaEhNm/eTEFBAXPnzlUbLEzXou10oyiK2pJueHiYqRrVJSUlqRkzxxcWDQ0NUVdXx8DAAEajkfLyaKvAAwcOAMdywUVRJDU1VXVktNlsDHqCXPqDd3COlL5fna9wTcFYbSgtLaWkpGRMWCIWckk06RAEAX9IQhDgH59ajs04dmQ/OsNFq9WqC6vXXHMNwWCQN998E4jOWkYbrY2mrKyMtLQ0XC6X2r5udN/Wl9pg94BIglGDRhPtIuQORKgpTT1pT1WXy8W2bdtQFIVly5aN+8w12T3c+PPNeILR8/7U2lK+fFXVCfd5NogLepxTImaeFSuxj7kdxgpHzhWhUIgNGzYAsHr16pOW2sdG9YsXLyYjI4PXXntNDRu8m4wW71jYxOVyqSKk0WjG2MGazWbVInhgYECNvx9PrHlHVlaWmlUTDAY5evQoHR0dY7bV6/Xk5+er1gYTxeif2NbGA89F3Tq0osAPr87GGhoa1x4vKyuLWbNmYTabcfrDfP25QxzodCAIAhaDhm/dWM2CgokroGMeQKNH5TqdjksvvZSGhgZaWlq4/PLLcblcNDc3q/1oJ8JgMGCz2cb82AMi9zy1DwEw6zUEwjJhWeYn75/PvPykSfcVw+v1snXrViKRCEuWLBlnY/HU9nbu/+dBIFrp+psPLeLymZkn3e90Ehf0OFMmHA6zfv16RFFk9erVBINBNm7cqGaSnOvwhcPhYPPmzVNahJVlmY0bNxIIBFi7di1Hjx6lt7eXdevWnXJxz+kSC6eMjnmfSLyTkpKwWCwnfF2Txd9PxPGjXIvFQmVlJTk5OZMeS5IV3vvLLezvcABQlGrmxc+uRA766OnpobGxkeP1ImbO1uXw4wlEKE23nDCtLxgM8tZbb5GWlkZVVRXvvPPOSV/L8X7oMWILu7EF29jMYXfbEI+83kDLgJecJCOfuaSMNZUn9wqK4ff72bp1K4FAgMWLF6sxfohaJ9zx221sa45aMWhEqMpKYE1FOh9ZUUSG7ey7nsYFPc6U2bVrF729vaxcuRKbzaZ6raxevXrKedPTzamkSca67qSnp1NYWMiOHTumbP51qowW79joe7R4a7XaMdkmUxHvqRyzq6uLI0eOjGnsMRmFhYVkZGRQV1eHy+XCZrNRVVU1aRiqdcDLdT/bpIYVLp+RyW/uPLYA7na7aW5uHtdqLzc3l6qqqkk/I+FwWA2PxEJAJ6K0tJSMjIwxPVJlWaapqYna2lp1u9FdkGw2myru0XaD4mm3lAsGg2zbtg2Px8PChQvHvF8DniBX/2QDdvfIwrZGICPBSJrVwJ8/vmRcqGm6iQt6nCnR3t7O/v371X6Sp+qGeDaJTdUnMuU6nubmZg4fPkx1dTW1tbXk5OQwd+7cMzr+ROLtdDrHFO8c30nnTMU7RigUUkfo/f39amGOxWIhNTUVrVaLLMsMDQ1NmMqo0+koLS1VLwYej4fExESqqqom/L++sL+bz/3lWDOOR+9YwLVzxr/nsdDc8RqSlZWF1WpFlmVVxCfzhJkxYwZHjx7FarWyZMkSNmzYoMbWYwVUxztGxrzvY8yZM0ftShVbYBZFcUz/2YTT8FUKh8Ns27YNp9PJ/Pnzx3Tf+uBvt7OpcUC9nW4zYNZp+O8rKrhxfu5Eu5s24oIe56TE/FGSk5NZtmwZbW1tambEiXykzxWx9nVut1udPUzGaBMxm82Gz+dj3bp1U/5Cy7I8JuZ9LsU7dnyHw6Fmo8RsdHU6HWlpaWRkZJCWljbhaDgQCPD6669P+VjJyclUVVWRlpY25v4HnjvIE9uio/BUi55/f3H1GN/xWHFQLC3xRIu6iYmJ5OTkYLPZSEhIwGg0qpWaixYtwuv1cvToUdasWYPVauVf//oXgBomy8vLo7y8fMzrVRRFtWwGVCtnSZIYGhpSM4hiF7hY/9njwzMnIxKJsGPHDgYHB8ckBFz3s430OvwMeI91V8pMMHDnsiI+fenUmmufLnFBj3NCZFlm06ZN+Hw+1q5di9/vZ8uWLarXxfniGHkq7esCgQDvvPOO2s6spqZmnGjBePGOxbzPlXjHmKi0PmbjG8tGmYqLY6yfayy7x+fzYbfb6ejoGJcPPprk5GRmzpypLgI6/WGufGSDWva+vCiBr61Ow+eNivjohVJRFLFarZjNZhwOx4Sj8ZgbZXZ2NomJiSiKwvr165FlmRUrVvDmm2+Sk5PDvHnzVFfEFStW0NXVRXt7u2oFXF5ePqbIKObNHzuPyy67bIxYT9Z/9vjwzInWhkZ3P5o5cyaFhYV89g+b2NHhwRMRCI6E9zWiwM/vmM/V1SeeQZ4pcUGPc0KOHDlCU1MTixcvJjk5mQ0bNqiLomfqWDjdnEr7utFpckVFRcyaNUsV79Ex7+PFe3TM22w2n5ULWiQSGSM0sUwWk8mkCnhaWtopvf8xcZusn6uiKDidTnWBczJimTabm4d5ZM+xhdXLchU+NNs6LrPk+AucLMt0dXXR2Nio5rqPxmw2k52djUajob6+ntmzZ+PxeGhra+Oyyy5DURTefPNNysrKqKqqIhAI0NjYSFtbG4qikJubS3l5uZorfvys5ES1Cy6Xa1z+/1TCM7Iss3fvXtXXfTAo8KcmA94I9LmDRN1w4OrqLH7xgQVndRAUF/Q4kxLL1y4qKqK6upqtW7ficDhYuXLlhA0Uzgem0r4uFr+NpT3C2AW0cynecExMY3Hw0b1OY6X1scYRp0NsRmI2m1m5cuWkWT0xzxe3243D4aClpWXCDJIYr/QYebn5mKj/4gMLuGb21EagMd/4hoYGnE4nGo0Gi8WCTqcbY8oFUF1dzaFDh1QLhx07duBwOLj88svV1xIMBmlqaqK1tRVJklRht9lsqoFc7MKYkZHB3LlzTxhaOZXwjCRJHDlyRHW+TE5OZub8JTyzq51/bq2j1nFsv2fbPz0u6P+hSLJEUAoSlsOE5TAROUJYChNRIsiKTCAYYOeunWh1WubNm0drWyudnZ3RLIjMLAQENIIGURTRClq0ohadqEMratFr9OhFPVpRO6kIXnnllaxbt45777132l/bnj176OrqYtmyZaSmpk6YbXJ89WViYiKlpaVnXbxjxJwWY4IRC/8kJiaqo/Dk5OQzTqmMVfU6nU7WrFmDxWJBURR8Pt+Y4ptYL9DR74vZbFZj23q9nu7u7jGhGVmJ9go9OBR9r9Kset7877Ukmk9t5tbf309DQwNDQ0Po9XoKCgowGo00NjaOC9EsWbIEgB07dkzoxxMMBmlublYvRtnZ2VRUVGC1Wjlw4ICag6/RaKiurp5y96rJwjOjie2rsaWd1wasHB2CYMCPXwJXUBl5jwy89aU1JJylbJe4oF9gKIqCJ+xhKDDEUGAIZ9CJM+jEEXTgDrlxh9x4wh48IQ/eiBd/2I8v4sMf8eOP+AlEAgSlIJIy+chruhAQMGgMGLXG6I/GiFlnxqw189fb/8p1X7+OOcvnYNVZseltJBoSoz/6RJKNyaQYU0gyJKERJ45hPv300zz66KPs378fn89HKBTC7XYzODg4piNODJ1ONy7m7XK5iH3Wrr/++rP2XkiSNKYhRqxwxmg0jhnxnWkv0uOJuWHabDYSExNV4R498jaZTONCJTabbdLYcSzDw+Fw4I/AQ3sEXOGoqM9NVfjCIgszZsw45cbZQ0NDNDQ00N/fj06no6ioiIGBgXHxfZ1Op14Ar7322gkveqFQSBX2SCRCVlYWFRUV9PX1UVdXp26Xnp7O3Llzx2XLnIxYD4DRGimKIsnJyfxpv5OdvRKJpqhoe0MRBgKC6sr44ZpCHrxxfNhrOogL+nmCL+yj39eP3W+n39fPoH+QwcAgg/5BHEEHg/5BBgIDDPoHCcvhk+/wPCY8HKbui3VU/awKre3EHnACAsnGZNJMaaQaU1WhTzGm0LS5ieBwENkn8/P/93Oef/75cWlyEI2bpqSkTBqyePHFFwGmPQUzFpPt7+9naGhoTGn96JjsdBAIBMaNuI8XQqPROKFwn66NcCgUorGxkX/saOK3tcdE9fpCmXWjogqlpaXk5ORMaAw2EU6nk4aGBnp6etT7MjIy8Pv9uN1ucnNz6erqUh/Ly8sjOzubjIyMceIeDodpaWmhubmZcDhMZmYmJpNJjblDNMQWW9CcCk6nk7179+J2uyksLFSdRu12O339/Xz5HQ96EWJp7oIgRG1/PcfCU499cCFXVU+/7URc0M8RYTlMt6ebdlc77e52Ot2ddHm66PH20O3pxhWa2O70bBEbPes0OkRBJBAJEJEjCAjo0JGkTcJqseLzRi1cbVYbgiigKAoKCrIiqz8ROUJEiUTDNnKYsBQmJIWIKBO3KXPtc9H9x26qHol6Xbh2u+j5Sw9Z78siccmpu9R5jnpo/X4r838/n0xTJjm2HPIT8kkWkvF2e5mZPZN1i9dh1k1c2BIzqIJoKOh0R8nBYJCBgQE1lHJ81kRGRgYpKSlnVFEbDAbHCbfb7VZHrBAtezcYDGrcd+nSpSQnJ5+1RexAIMD//GUHL9YdM9H6WBXMTx1vKpaUlERBQQFpaWknXRPweDw0NjaOsypYsmQJCQkJvPHGG8CxEXusOjQm7qPf5+O7Mo3GYDAQDAZJS0tj7ty5kxZAybJMQ0MDDQ0NGAwG5s6dO24AIMsKa3/wNgZRQZYkwuEwiqIg6E0kmPUc6or+T5LMOl774upprx6NC/o0IskS3Z5uWlwttLnaaHO10eHuoNPdSbene1KBO1VMWpM6Sk0yJKmhCpvehk1nw6a3YdVbsWgt0RCHzoxJa8KkNWHUGDFoDGp8OyJH+Ni/P4bdZ8ems+HyuPBLfq6ruI5LNJfQ29s75S5BxxORI4SkEAEpQCASUMM+P/3uTzmw5wBf/tmX+dXDv6J2by23P3Q7ukwdz3/3eWrfqlUvHMd/BtOvTSf9urH9G2OCXv27E09j003p5NvyybPlUWAroCixiKKEIrKN2ax/cz0QtdxdtGhy573RxAp2RjeugNPPax5NLHx0/M/okn2dTqfGuEePuPV6PTt27MBut0/JvGw6CEsyd/x6KzvbHAAYNfCr9+RTVZBJT0/POFGOUVBQoGbtTHYhdTqdYxawAdasWUNDQwN2u53LL7+c4eFhuru76e3tJRQKodFoyMjIUJtzx2YhkUhEFfbR72Vubi59fX0oisKMGTM45NLz523tOHwhlpem8sGFmXQ3HcXpdJKXl0d1dfWkF8hvvnCIN0dMySRJYtjlwWA08cePL+UDv91BlyPqFHnVrCwe+9CpN08/EXFBP03cITe1Q7UcHTxK3XAdjY5Gmh3NBKQTd0KfDJ2oI8OcQbopnXRzuhpiSDGlqOKdZkojzZSGSTt9vQv39e/jgc0PYNFZ8Pv9hIIhjGYjvrCPu013M2/WPEpLS6fteAA33XQTWq2WxsZGFixYwJe+/CUODByg19VLcjiZZCFq3qTT6bAmWhHMAhFjhLA2jFfxMhQYUkNSdp+dIzuPsPUbW08q6JMhCiLp2nTSSCNTk8nK8pWsqlxFpjlzXIjA4/GoYZTBwUG1tdzo1LaphhZgbNn76J/Ri25arXZcmCQhIWFcm7kYo6thp9KDdbpw+sLc8Ogm2gajeejFNoV754mUlUa7J4VCIXp7e2lsbJzQc8Zms6kLwsfPZGJNyEcTc2UcnYqoKAqDg4P09PTQ29tLIBBAFEXS09PJzs4mKytL7RDV1tY2Zq0lLS0NQRB4+bCdV7u02MwG9FoNw54AZiHCvQu0LJw/96QOnU5fmK/84wC1PS4EIOD38qVLi7h5ZTWbGgb44OPb1W2f+sRSlpeNr4E4XeKCPgUG/AMctB+kdqhW/en2dp/yfjLMGRQlFJFvyyfflk+uLZccSw451hxSjCmIwrnvR7ilawvf3vFtLDoLPq9PnSL68fNQxUOsXrp6Wo4jyzIulwuHw8Hy5csJBoPk5OTwwEMP8LTnaRyKA0EUEEWR6wqu45NzPznlNL133nmHyy+/HIfPQbe3mx5PD12eLjrcHXS4O6IzJVcHEqe2EJygT6AiqYICUwEZSgZJgSSsYSuiIGKxWFQBT0tLO2kcOhKJqNWTo4V7dBaHRqOZULhPZYQf6xqVkZFxVvumTsau1iFuG9WW7apSM9dketDroxYDxcXFaLVagsEgfX199PT00N/fP+G+YimbaWlp2Gw23n77bURRVNMPdTodnkCYWgdkFpRSU5nLrJxjeeKKojA8PExPTw89PT34/X4EQSAtLY2cnByysrLQaDQ0NTWNWSj9yVED3kAYnRhd6JQkiYig5cH3zGZd9dRK9xVFoWXAizsQoWXfZspLorUOAPf+bT/P7ukEYHZuIs9/esVp+8ocT1zQj0OSJRodjezp38Pe/r0csB+gy9N18ieOkGpMpSSphMKEwrHibc2dNIb7buIKufjAvz6gphsGAgEcfgeZmkzuq7qP6urqU44pjxbvWLqg2+1WrWLvuOMOnnjiCR588EFyluagv0FPsik6Kt/56E7a32nHoDEgMP5Dfv/993P//fePuS8m6DGfj4kYdg7z0oaXCBqDWPOtdLg7aHW10upspcPdgcLUPusWrYU5aXOYnzWfhRkLqU6rHvN/jZW9Hz/iPr56cqLFSZPJdEbpkrF8a0mSWLNmzbRnzEyVx9Y38Z1XjplkfWZ1IWvS/fT29qLX69UmG7ERuCRJ2O12ent76erqmrChx+jMFoi+h5Xzl/HJP2zHHZKRZNBrRa6alck3bpo/YeaLw+FQxd3r9aqNubOzs0lPT+ett94C4Ju7BCy66HqRIAqYTCZ8YfjC5RXcvPDUc8jXr1+P0Whk6dKlAHQ7/Fzyg3cIRqKv8/5rqrhr9fTMgv/jBT0shzk8cJgdvTvY2buTgwMH8YYn9pgejVbQUpJUwoyUGcxInUFFcgWlSaWkGFNO+tzzjfUd6/nBrh9E88/9AfSKnk8VfgrBKaDVaqmoqBjT5WY0kiThcrnGWMLGxBui8eTRaYL79u3jzjvvVEdFq69eTcmNJRReXogiKyCAV/ZyWeJlXJZ8GaIootFoJvyB6MVj165dfOxjH6O+vh6NRqN2sT9++76+Pvbt20dRUdGY9nVDriH2tu9lY+1GeiO9dEld9Eg9BJWTOxaKgkihqZAyYxlFYhHZkWyMYnREHSt7P164z1ae+549e+ju7mb58uXjvLrPJYqi8Jm/7OVfB45lqfz+I4uZn6WntrYWu92OwWCgvLycwsLCMZ+r0aPqvr6+Sf3eAZ5uFGj26dHKIWQF/JJAUFK4c5aRD66ecUI7YJfLRXd3Nz09PWrFakpKCsPDw/y+Fto8YNJEZ016g4GgLPCbOxdRkTm5T9Bk7Nq1C5fLxaWXXqre96PX6vjpW9GKXK0o8MJnVjIz58wznv7jBF1RFNpcbWzp3sKW7i3s7N2JL+I74XMMGgOzUmcxM3UmVSlVVKVUUZJYgk5zfpW+nwnOoJPndz7PUP8Q7136XgpyC/B4PBw6dAi73Y7VamXGjBkYDIYxRTonEm+z2TxmxOrxePj973/Pjh07ePDBBwH44eEfsvnhzcy9ey45S3NAAU/Ew/Vp11OTUIMkSciyrDZOjv3Isszrr7/OT37yk3Gv5Te/+Q2ZmWfWWEBRFByyA7fJTXekm2ZvMy2BFrzKiS/2AgLltnKW5yxnbdFa5mTMQSee/c9JzEb4RBWy55JgROLjf9ilug4mm3W8/PlVZCeaGBoaora2lsHBQUwmExUVFeTl5U04YHC73fT29tLb26sakcV4eI+AKICkwGAAEARkJSrEV+fLXFNuoaysbNJ9jz5GbOTucrmw++F3dQKBCEiARoCVeXq+fE012dnZp3wxrq2tpbGxkWuvvVZ9bliSueWxraq//IKCJP5+z/IzDr38Rwh6WA6zu283b7e/zfrO9ScNoaSb0pmXMY/5GfOZnzGfyuTKi0q8J6Kzs5O9e/dSWlrKzJkz1ZG3w+Ggvr5+3CKWXq9XhVuv16uNBkaL9+jnjB6tWq1W9WdD/wYe3f8oJq0JrajFG/Zi0Bh4/MrHSTZO3NkmhqIoEwr98ffFyrh7e3vHpazFiHWbP9nxhuQh2iPttEZaaYu00S9NHP+NYRAMVOormW2aTbWlGpsuWrBzopnHZI9Ndr/f72fTpk1TavRxLhnyhrj6Jxvoc0VnOtW5CTz1yWVqleTAwAC1tbUMDw9jNpuprKwkNzd30vMPBAL09fWpnuk/OSjgDoMjBJIc7RIkK5CoFxAFhc9XKyQZUHupFhYWnjRl1Ov10tPTw56DRznqAE8Yim2QZ4nu32q1UlFRoY7+FUXhcLcLuzvIjOwEshLHr3fELraXXnrpmHWhxn4PV/9kA2EpqrM/vHXuaYV0RnPRCnpQCrK1eyuvtb7GOx3v4A67J90225LNkqwlLMlewqLMRWRbTv0qfCEzPDzMpk2bgGiRRmzhLvb/NxgMWCwWhoaGxjzPYrEQCATGVB3q9XpVrEeL92TxYUVR+HvD33n66NN4wh5Kk0r5/ILPU5lSeUavKRAIqOmEdrtdvbjEytgHBqIjx4yMDAKBwIRl74FAAFmW0ev1LF26VK0mHH3xcAQcHBw4yN6Bvewd3Euju3HSeLyISKWlkgXWBcw1zcUqWie8IE21KfRE6HS6KV8ITufiEXtsqt+P7c2D3P6bbYyskbK4KJk/fWwpJv0xYe3v76e2than04nVaqWysvKEI+FYGuPeAXiuVcAZFBCE6LsuAhkmCMnw3mKF6uMiTwUFBcycOROdToc7EOato/10O/zMyU9iWUkqmlEj5J6eHibTIYvFQm5hCT/aMkBdr1u9mHxoWSEfX1k85txjPW0nKlz73qu1/OKdpui5pZh589416E7Q1elknLGgC4JwFfATQAP8VlGU70yy3WJgG3Cboih/P9E+T1fQJVliR+8OXmx6kbc63po0Fm7WmlmSvYTlOctZkbOCfFv+f4yAS5I0xsfbbrePGZmO/iLH4tCRSASfzzdhFaZer6eqqkoV79NdiIsVKOk1p/f82Cg8VtQzulM8RKv1jEYjoVBonOFUrGLTZrPR399Pd3c3l19+OSaTSW2eMZnF7vE4g0629WxTQ3q93t4JtxMFkUWZi7iu5DrWFa3Dojs2cjuVmcf+/fuBaMFOWlrahNuc6PmnO2g7lQvBG80+frptUH3uFeWJfP2KAvXzFXtOf3+/2srOaDQyc+bMSePge/fupaOjk90D8HSzBlkGg0bBplXQCBCU4cMVCoWThLydQXii3YwvLCADGkFgUVEK37159pg2eT6fj23btuH1etVajKGhIRRF4d8dsLVfQ6rVgN6gR5ajZf6P3rGA6txjuf/BYJDXXnttwjRSdyDMqu+9jcMXnTk+fFM1H1h6+o3Wz0jQBUHQAPXAFUAnsBO4XVGUIxNs9zoQAH433YLe5eni2fpneaHpBfp8fRNuk23J5pL8S1ibv5ZFmYsu+hAKjBXvWMzb4/FM+UssitH0vIlCJS6Xi0OHDuFwOEhOTqa6upqkpKSz+4KOI1Za39HRMU7Aj2eisvfh4WEOHz6suvhBdMr91ltvMWvWLEpKSpAkSe2Us3bt2lOqtlQUhUZHI+90vMPbHW9zcODgxOemMXJ54eXcXH4zCzMXTnlwERtBxsJkp0NsRjAV8T+Tx19vl/hHy7HXdUeZTM0pLHMcP/MIBoNqrv5bXbDZrsOoFRAFAW9IIsMgcdcMhclC0v9sEdg/CGZt9HOu0+kIyAIP3zSbNZVjTywUCrFz506GhoaYNWsWubm59Pb28vGnj+INRtAKsZRTK8O+MB9aVsjda8Zmrbzyyivk5+dPaF386NuNfP/f0bRJm0HLv7+4mpyk06s1OZGgT8XgYQnQqChK88jOngZuBI4ct91ngWeBaUuMlRWZzV2b+WvdX9nQuWHCaW6+LZ91heu4ougKZqbMvKhH4Wcq3snJyeNCJSfKxkhOTmblypV0dnZy9OhRNm7cSH5+vrpwejaYrGflaDQaDUlJSeOqJycS4uTkZDweDw0NDSQmJpKdnY3FYiEhIYGenh5KSkrQaDTMnz+fTZs2cfDgQRYsWDDl8xUEgfLkcsqTy/nknE/S5+3jjfY3eK31Nfb271U/swEpwEvNL/FS80uUJZVxW+VtXFdyHVa9ddJ9+/1+9u/fT1JSknoxOh1EUYw6Zp6mn8tUuRoQn9nH33dH16/+1qxhxYKZrCxNnvCiEIlEaGtrUzNQwuGwum4jSRIGgwG73Q7A2hywaMNss4sEwrAwVeGSnMnFHKDZBfqRgbgsywSDQXwR+NubuzAMH/PaiVXe1tTUsGfPHg4fPozP52PWrFkkJXbidwQwaEAUFBhJtDXrx8fpLRbLhP7vAB9dUcTfdnXQNujDHYzw8L+O8ugHpv45mypTGaHfAlylKMonRm5/CFiqKMpnRm2TCzwFXAo8Drw00QhdEIS7gLsACgoKFk7WtiokhXi+6Xn+cOgPtLvHf7FTjClcU3wN15Ved9GKeCQSGZPnPVXxjo1SrVYrPp+Pvr4+0tPTWbZs2RmfT319PS0tLYiiSEVFBcXFxadt/RqrnnQ6nbS1tZ1w9B3r7B4T7lO9mIxuX7dq1SqsViv19fXU1dVxxRVXqEU9DQ0N1NbWsmDBgjH9I0+Xfl8/r7S8wvNNz9Mw3DDucYvOwq0Vt3LnzDtJN4+1Ohh9zqtXrz5tn/RzjS8U4aZHt1DXF/1/GrQiv//oYpaXTh7KkmWZjo4O6uvrCQQCpKWlUVVVRXJyMuFwmDfffJNwOEx+fj7z5s0b87zjZwrBYJD6+noGBwd5vFagyxu1KAAQRIGQouW26gTmJwXVz1zMxiEtLY20tDRaWlpoaGzm3/0WdvZKOPxhRCFqixtbF3jyE8vGLY7u2bOH4eFhLrvssglf5/bmQW779Tb19tN3LWNZyanbbZxpyOVW4MrjBH2JoiifHbXNM8APFUXZJgjCH5hE0EczUcjFF/bxt7q/8acjf8Lut497zoqcFdxaeSur81afkzSx0yUsh9nZu5Neby8liSXMSZ9zwgrR48Xb4XBMeqWH6KhwsjBJbBTm8/nYsGEDJpOJlStXnpFZ1Gi8Xi+HDx+mr68Pi8VCdXW1uggUkSN0e7pJ0CeQZExSX9vxBTixUddE6PV6SkpKyM/PPy1/lMmIta/T6XSsWrUKv9/PO++8w+zZsykqKvr/7Z13eJTXmbfvM129SwihghDqWFSBEEVGBgyOQzaOiY1L1mUdJyHOJus4LomTfCku2c0S2/HaOCaO119c4m7HxibU0EyREAgQIAkkgQQS6m1Gmpmzf4w0kkBl1IvPzTUXmnnfeefMuWZ+c85znvN7AEf4ZM+ePTQ0NLB06dJ+2632hJSSk1Unefv023xc+DHN1uYux/UaPWti1nB38t2Eezm2t7eXYRuqH5eR5GKtmZtf3EtJleN9ept0fPT9RUQG9P6j1L5VPz8/H4vFQnBwMPHx8U6r5O5KzPXGgbOV/PD1LFpaLOjbYu5uOlifJPEyakhISECv1zs90C0WC1LCZbsb24osHKuU+Bo1tAgDlY2OPPi4EE/+39dmMCfy6uysU6dOcebMGVavXt3jQOcHb2TzwRHHDvS5kX787f60fg9IByvoacAvpJQr2+4/AiClfKLTOWfBueUvEGgC7pNSvt/TdTsLutVu5b389/hj9h+pNFd2Oc/L4MXXY77O2ri1RHi7ZlQ/mtS31PPgzgc5X38eq92KTqNjTsgcHk97HJ1G1y/x1ul0XcS6c5ikt5GxzWZjz549NDU1sXjx4mEZ3ZWXl5Obm0tjY6MjiyTQzLPHnqWhpQGb3UayezIr3VbSau7dBjgkJMTpzDfcIYHKykr27dtHSEgI8+bNY/v27ZhMJtLS0pznNDY2snPnTmex7KGe/TW0NPBR4Ue8nvc6Z2vPdjmmEzrWxq3l5ik3czLrJBEREaSkpAzp648UxZVN3PziXmc6Y/wkL978dprTP7w3bDYbZ8+eJT8/32mH214MuvNaiCt8UVjJy7vPUlrTTHyQkVluVXhrr/5MJiYm4uMXwMPvHiO3rJHKZht2u8RNBwFuGjy9vKhuauGWeeE8kNl90fT2tOBrr73WWR7vSi7UNJPxu+3ONMa//tv8Xmcv3TFYQdfhWBTNBC7gWBRdJ6W8urqA4/xX6McIfV/pPp4++DT5NV1rHAa5BXFn4p3cHHdzl+yAsc4rua/wxqk38DY4doRZrVbqW+q5NeRWpsqp3Yq30Wi8aqTt5eU14BFqe8GD1NTUQW++uZLOeei1tbUUFhZSb6/npbqX0KBBL/ROn5hUYypL3ZY6n6vVars4FI5GGKHd1CohIQGr1Up+fj4rVqzokrlTXFxMTk6Oc9F0OLBLOztKdvDysZc5evlol2Nuwo2VPit5bNVjmAxDa706kuSU1HDzC/tosTlSNGdF+PLaPfPxMLr2w93uc15QUNDF8mH16tWDnnHW1taSk5PjdM8E2HcJPisR+LobuNRkp8Vqxy4lge5agnw8qG5qZV1qBN+9Nqbba7anBvf1vXvk3aO8fsDhTJkRF8Qrd6X2q+2DWhSVUlqFEOuBz3CkLW6SUh4XQtzfdvyFfrWmDZu08djux/iw4MMuj4e4h3DfNffxtZivDTi9bTTZW7oXo9aItEvqG+pBQots4cClA0QFRRESEnKVeA+lh3VRURElJSXExsYOSsztdnuPfiXtg4B2v5J8cz42bBiFI7YtNAKj3cjRlqOsCV3j9An39fUddLm1wRIdHU1NTQ15eXnEx8c76152LlMWERHBpUuXOHnypDN2P9RohIZlEcu4NvxaDl06xB+P/JHDlw4D0Cybeb/mffI+y+PX6b8edL7+aJES7ssTX5/Bf/zNkXaZXVzDXa8c5NW7UzHp+xZkvV7vXKvJz893Frb+5JNPyMzM7NHT3BV8fHxYssRhStfc3MyxY8fIPVmOVjgyXgwSLGgQCCx2DS1WiU4jWJ7Y83eqfYDSm5UBwH1LpjkFfcepCnadrmBJbFCvz3GVUdtY5D3NW0Y+3pGL6aZz494Z93JH4h1Dah070jz6z0fJqcjBQ+9Ba0srQgga7Y18M+6b3D3j7mF97ZqaGvbs2UNgYCCpqakuhQvsdnuX2pPtG44aGxudwt0es++u2rtGo+HDgg95/sjzeGo9MVvMjvUCDWh0Gt7/2vvD+p4Hgs1mY/fu3TQ3N2Oz2QgMDHSaKrVjsVjYuXMnRqORxYsXD/sPkZSS1w++zsZTG6m0d4QddRod30n5Dvck39Njmb6xzl/2nuPnH3ZM6L+aMpk/3DKz3+Gs0tJSDh92/OgJIYiIiCA2NnbI1loefucoewsuY8CGXdppsGqoamzFw6gjxNvIj5bH9VmBaPPmzUyePJlrrrmm1/N+9NYR3s1yZAPFhnjyyQOLu+TG98Zg0xaHBau9Ywq1auoqfjz3x1et9I9H1sat5Uj5EcxWM0aDkabWJkw6E6umrhry16psruTFoy+yr3QfJq2JJHsS6Z7pzJ49+6ovi6tFg9uFOzQ01Cncnp6evQpa6qRUXhQvYhM2Z4HiupY6boi4Ycjf81Cg1WqZO3cu//znPym0FPLW2bd4pe4VloYvZU3MGtx0bs5qNQcOHODUqVMkJCQMa5tqamrwLvfmybgnyXPP47kjz2GxWbDarTyb/SxZl7J4aslT+BiHv5DFUPOthVGYW2080ebO+GFOKcFeRh67IaFfoj558mSKi4s5WFhBQYMOXWkJ15wpYU58FDExMYNOpV07L5x9BZUInQF3nQZNqx2jTscv1ySSHhOIUdf3D6qHh0efI3SAh1fFszn3Ik0tNk5fauDdrAusnRc+qPbDKI7Qo5KiZMTPIvhZ2s9YHrl8VNowXOwr3cfLuS9T1lBGrF8s96fcP+TTZqvdyre3fJsLDRfw1HtS31hPs7WZ2xJv446EO7otqNBdtffON09PzwHHJrcVb2ND1gakdJSumxE4g8fTHh+TdsLtvHfsPTYc2YBA4G5yx66xk+ifyH9m/KczKyknJ4fi4mIWLlw4oIpOrtDa2squXbuQUrJ06VL0ej1na8/y090/7RJfD/MM45llzxDr1/2i3FhGSsmj7+Xy+oGONORvL43m4evj+yXqz3x+glf3FiI0GowGI9Jm5eaoVmL9NEydOpVp06YNylJ468lLPLc9n4p6C6E+Jn54XWy/ilNkZWVRVVXFdddd1/d72XqG3285DcDUQA/+8aOlXWwJemLMerls3bN1XI44xgKHLx3m53t+jofB4bViMVuwYcOOnQe8H3AKUufdk+0bcTqnNw4lTa1NFNQU4Gv0Jdx78KON4URKye2f3E51YzX2Fjs6vQ53d3caWxv57eLfkhLkyC5p9x+32+1kZGQMS78dPnyYsrIy0tPT8fPrSIez2q08f+R5Xjr2kvMxP6Mfm1ZuIsav+4W5sUyrzc76v2bx2fGOnd7fXxbDf6xwbbBTVtvMLRv3o7G3YG1pxcvLC4tN4mXU8tM0d8pKS9HpdERHRxMdHT3gtSkpJS02Owat63427bTvb7jhhhv6DNPVm1tZ+OQ26s2OaMULt8/m+uTQPl+jN0Ef1RUqJeYDp9ZS69yFaLPaEBqBQWtAaiUJMxJIT0/n+uuvZ/ny5SxYsICkpCTCw8Px9fUdtvRAd707M4JmjHkxB7BKKxXNFXi6eaLT67Db7U5nvbKGDo9vnU7HrFmzMJvN5ObmDnk7ioqKKC0tdW6k6YxOo+OB2Q+wIWMD7jrHTKfaUs29n997VcrjeECv1fDsrbO5LqFjYfHZbfm8caDnXcGdybtYj0aAu8lhAtfa2oqbXktVk5WYhBlkZGQQFBTE6dOn2bp1K2fOnOm1IEpPCCEw6rQDSll1dWEUwMuk544FHeuIr+13rR96Y3RTDhQDJjnQ4RdhtVvx8HRsZZdGSWJwIrFTY/H39x+2CvATAZ3QEe4VTrO12blJS0pH8Y2pPl3Nlfz8/Jg+fTolJSWUlZX1cMX+U19fT25uLkFBQb3WdM2MzGTjio1OUa80V/Lv2/+dptbePf7HIgadhudvm01GXMd62SPvHeOvX/QtZiFeJqR0CK6npycGowGrXWLUaXA3OGqyzp07lyVLluDv709eXh5bt26lsLDwKrO24aI/gg5w24JIp33B7vzLnLvs2vN6Qgn6OCXYPZhvJX2LZmszNZYa6i31uOvc+eGcH45208YFQgjWz3K4V9RZ6mhobaC+pZ70yendxqinT5+Or68vR48e7dNT3RVsNhuHDx9Gr9cza9asPkeDKUEpvLD8BUxaR0ZHYW0hv9r/q0G3YzQw6DQ8t242SW3Ve6SER987xqbdvc86EkK9SA7zoaa5FasdzK12Gi1Wbl8QiUHXIWU+Pj6kpqayaNEifHx8OH78ONu2bePcuXODsi12hf4KepivGxlxHXa7z23P7+XsvhnXfugKKK4rJqs8Cw+9B2mhab2aPSmu5lztOTaf3UyVpYpFYYtIn5zeY3pgQ0MDu3btIiAg4Ko0x/7Svti6YMECgoJcz+76sOBDHtv9mPP+a6tfc8b7xxs1TS3cuekAR893bO753Teu4ea5PYfsmlqsvLa/iC0nLuFu0LF27hRWz+i9tkFlZSV5eXlUVVX1WT1pKPjss8+YNGmSy7t89xdWckubx4tWI9jxYAbh/j0nE4zZRVEl6Irxxrlz5zh27FgXD5j+0p5PHRMTM6B0yB/v/DGbz20GIC00jY0rNg6oHWOBOnMrd//5IIeKqgGHoP1+bQprZg69f01FRQV5eXnU1NTg4eFBXFxcrzVJB8ru3bvRaDQsXLjQ5eese2k/ewsc+w/uWBDJr752tQVvO2N2UVShGG9ERUURHBzMiRMnejVQ64mmpiZycnLw8/MjLm5gqazrZ61HKxyziH1l+6ho6tnsbKzjbdKz6a55JIQ6wi82u+QHbxzh5T7CLwMhKCiIxYsXk5qailarJSsri507d1JWVjbgIiDd4enp6XLIpZ3vZnRkLb2bdZ4GS/8Xc0EJukLRb1JSUtBoNGRnZ/crJmu32zl8+DBCCGbPnj3gKX+kdyRzQuY47+84v2NA1xkreJv0/OWueUwP7ggX/urjE/xPW9m2oSYkJIQlS5Ywd+5cpJQcOnSIXbt2OQ3ABkt3ZRv7Ij0mwPn+G1tsfHCk95rIPaEEXaHoJyaTiZSUFGpqajhz5mqf855on+6npKQMyocEID0s3fn3qapTg7rWWCDY28Tf7k9jbidb2qc25/HqvnPD8npCCEJDQ8nIyGDWrFnYbDYOHDjA7t27e7V3doX+Loy2t2fd/A4/ofezlaArFCNGaGgo4eHhnDlzhurq6j7PLy8vp6CggMjISEJD+9480hftvulAjzVNxxu+7gb+9575pHUq+vD4B8d5ZuuZIQ2JdEYIwZQpU8jIyCAlJQWz2cz+/fvZu3cvlZWVfV+gGwYi6AA3pkx27hQ9eK6a0prmPp5xNUrQFYoBkpycjMlkIjs7u9cNLGazmezsbLy8vEhKShqS1+4scIKhXdQbTdwMWl761lxmRfg6H/v9ltP87INc7PbhS+DQaDRERESwbNkyZsyYQWNjI3v37mX//v0u/WB3ZqCCHuhpZOG0jh+zLSf6HwJSgq5QDBCdTsfs2bOdVZy6Q0pJdnY2NpuNOXPmDFnlqM6lGUM8htbzfrTxNOr4y92ppMd0iNtr+4t56J2jWG3Dm0eu0WiIiopi2bJlJCUlUVtby+7duzlw4EAX7/Te0Ol0GI3Gfgs6wIqkDjfHj4+W9vv5StAVikHg7+9PTEwMxcXF3S6q5efnc/nyZZKTk4fUV333hd3OvxMDEofsumMFb5OeP/9rKmtmTnY+9vbh83z/9WzMrcO/61Or1RIdHU1mZibx8fFUVVWxa9cuDh061Gv923Z6KxjdGyuTQrqEXQoq+ncNJegKxSCJi4vD29ubnJwcLBaL8/GqqipOnTpFWFhYlwIag6WwppDs8mzAEW5ZMmXJkF17LGHQafj92pnc0slW9tPci9yycT81TS0j0gadTsf06dPJzMwkNjaWiooKduzYQVZWVq8j8IGkLgIEe5nIjO/YOfpu1vl+PV8JukIxSDQaDbNnz6appYlndj7D+q3reXjnw7y17y3c3NyYMWPGkL2WlJJnsp/BLh2hh/SwdALd+leTcjyh1Qie+PoM7l3U4a9zpKSGb764nwsDWDQcKHq9nri4ODIzM4mJieHixYts376dnJwcmpuvboeHhwcWi2VA5mA3zZni/PuTYxf7tSCsBF2hGAI8PT3ZotnCR5c+oqCqgIMXDvJ61evUTqodUpO0d868w9birc777X40ExkhBD/9SiI/vzGR9k2dpy7Vs+a5PWQV92/BcrAYDAYSEhLIzMxk6tSpnD9/nm3btnHs2LEuHj/u7u7UWqC4vP/tWxobhIfBsdZy9nIj+eWuh12UoCsUQ8DJqpOcM5/DS+eFtEj0dj0eJg/ePPumczQ9WA6UHeCJL55w3v+XmH8hKWBosmbGA3elT+W/185Er3Wo+uUGC7ds3M/fjw6dA6arGI1GkpKSyMzMJDw8nKKiIrZu3crx48c5e6mWhz8pYkOu4I6/HGH9X7OobLD0fdE2THoti6d3+PvsPO16XrwSdIViCGj3UPfw8EAIgU6vw9PkSZ2lDrN18O6Mhy8d5ntbv0eL3RE7jvWL5ZH5jwz6uuONr80K47V75uPn7pj1tFjtrH89i7/sPTcq7TGZTFxzzTUsW7aMsLAwCgrP8u1NuzlzsRaTBtx0kFNSw+MfdJ8F1ROdi0a3e7y4ghJ0hWIImOoz1VFwRICXt6OAttlmJsg9aNBFzz89+ynf+cd3MNscPwzB7sFsyNgwroupD4b50QG8/710ooMc+d5Sws8/PM7jH+TSYh3etMaecHd3Z+bMmUxJmku9TYvWbkUIkHaJj5ue3NJaLtW5/sO+INrf+fehc1Uu5+CPqqBXNg9sJ5ZCMdaI9o0mY0oG9S31NLQ2UGepw263s37m+gG7+bXaWvndwd/x0K6HaLY6Ft4C3QLZtHLTuKgKNZxEBnjw9v0LmRnu63zs1X1F3P6nL6hqHJkMmO4wGN0wmYx4e3thMBjQthWWFjhK8LnK1EAPAj0dtVHrzFbOVrqWMTNqgl7RXMGN79/IRwUfDdu2XoViJHlw3oM8Ov9RFoUt4ivTvsIzy54hNTR1QNc6VXWKdZ+s49UTrzofi/KOYtPKTUR6R/byzC8P/h4G/vpv87lhRoeVwoFzVdz47G6OnXdtE9BQEzfJCy+jHrNV4u7ujtFgpN5iZYqfG2G+rs+ohBAkh3WU6My94Nr7GTU/dLepbjLmFw7LyGvDr+WR1EcI9Ry8x4VCMZ4xW838+fif2Xh0I1Z7R8pbRngGv130W7wMQ7c5aaIgpeSFnYU8/Vke7XJm0Gn49Zpk1s4b+ZnMidI6Hno7h8YWxwYofw8Dv1+bQmSAR7+u8/TmPJ5vc5xcf20MD6502C335oc+PNWCXcCgMTj/3l6ynb2le7kz8U7umXEPHvr+vXGFYrwjpeTTs5+yIWsDZY0dWRtGrZEHZj3A7Ym3oxFqyas7hBB8J2MasSGe/PubR6g3W2mx2nnonaMcvVDD419J6lKibrhJnOzNu99N50RZHTqNICHU27n7sz/EdLITdjV1cdQ+IdN8p7E2dq3zvsVm4aVjL7H63dX874n/HZcFcBWK/iKlZNf5Xaz7+zp+8s+fdBHzGYEzeOvGt7gz6U4l5i6QmRDCR+sXERfSMYt5bX8x617aT1ntyG1CAscMYWa4L8lhPgMSc4DIgA6L5fM1runhqJegO3zpME8ffJoTlSe6HPc1+nJbwm3cGn8rPkafHq6iUIxPrHYrW4q28PKxlzlV3dXP3M/ox/dmfo+bYm9Cpxm1SfS4panFykNvH+XjTvnpvu56/uvmFDITxo+R2cVaMwuecGwiC/AwcPhny4EhqCkqhLge+AOgBf4kpXzyiuO3AT9pu9sAfEdKmdPbNTvXFLVLO38v/Dt/yPoDl5q6Ghy56dz46rSvsjZubbfV2BWK8US1uZr389/nzVNvcqGhaxEDg8bA7Ym3c++Me1WsfJBIKXlxVyFPb86jc8bf3elTeej6OEz6oXG9HE5abXamP/YpABoB+b9ZjUYjBifoQggtcBpYDpwHDgK3SilPdDpnIXBSSlkthFgF/EJK2WtZ9O6KRFtsFt478x6vHH/lqg87wOzg2ayNW0tmRCYmnanXdisUYwUpJUcqjvD26bfZfHazc3NQO246N26afhP/mvSvE84Kd7Q5cLaKB17P5mKnHPC4EC/++5szSZzsPYotc43kn3/mrC+a8/gKfNz1gxb0NBwCvbLt/iMAUsonejjfD8iVUvZatrs7QW+n1d7K5rOb2ZS7ifya/KuOe+m9WBG1gq9O+yozg2eq+KJiTHK+/jwfF37MRwUfdfEvb8fH6MOt8beyLn4dfia/bq6gGAqqGlt48G85bMsrdz6m1wp+tDyO+5ZEDzjGPRLM/fUWLjc4BgAHHssk2Ms0aEH/BnC9lPLetvt3APOllN26AgkhHgTi28+/4th9wH0AERERc4qKinp9bSklhy4d4s1Tb7K1aCtWebVzWbB7MMsjl7MicgUpQSloNWN/KqWYuJTUlfB50edsKdrC8crut3snByTzzfhvcn3U9WqmOUJIKXnti2J+8/cTmFs7NvjMjfTj92tnEhEwuBqvw0Xqb/5Beb3DB2b/I5lM8hm8oN8MrLxC0FOllN/v5txrgeeBRVLKXreB9jZC746KpgrePfMuHxR8QEl9Sbfn+Jv8yQjPIGNKBgsmL/jSbo1WjBx2aef45eNsL9nO9pLt3c4oATz1nqyMWsk3Yr9BcmDyCLdS0U5BRQM/evMIOZ02HrkbtPx4ZRx3pkWNudF6yi8/p7a5FYAjjy/H190wMiEXIcQ1wHvAKinl6b4a2l9Bb6c9HvlB/gdsLd5KjaWm2/MMGgOzQ2aTPjmdhWELme47fcBbsBWKzlxuvsze0r3subCH/WX7qTJXdXueTqMjLTSNr0R/hWURy9RofIzQarPzx+35PLstH1unFdM5kX48ddMMYoLHxoK01WYn9qefYpcgBJz61SoMOs2gBV2HY1E0E7iAY1F0nZTyeKdzIoBtwJ1Syr2uNHaggt6ZVnsrhy4e4vOiz9lWvK3HLxY4Ru/zJ81nXug85oTMIco7SsXeFS5xufkyR8qPcODiAQ6UHaCgtqDHcw0aA2mT01geuZyM8AyVcjuGOXq+hgf/lsPpSx2bdvRawbeXTOO7107D3TC6KaPnq5tY9NR2AIK8jBx87DpgaNIWVwMbcKQtbpJS/kYIcT+AlPIFIcSfgJuA9qC4tacXbGcoBL0zNruNY5ePsb1kOztKdlBYW9jr+d4Gb1KCUpgdMpuUoBSSA5NViEaBzW6joLaAnIocsi9lk12ezfmG3suA+Rn9WDxlMcvCl5E2OQ13/diMxyquxmK18fz2Ap7fkU+rrUMLJ/uYeOyGRFbPmDRqM/vNuRe5/7XDAKRG+fPW/WnAEAj6cDDUgn4lFxsvsq90H7sv7OaLi19Qa+nd3EYrtMT6xZIYkEiCfwLxAfFM952uvpwTGKvdSnFdMSerTpJXlcfJqpPkXs6lsbV3Zzu9Rk9KUArpYemkT04nzj9OzfbGOacu1vPwu0fJLq7p8nhadAAPr4onpZOr40jxy4+O8+c95wC4b0k0j65OAL6kgt4Zu7RzpvoMX5R9wcFLB8kpz6Ha0ndpKIEgwjvCIfD+8Uz3m06sXywh7iEqHj/OqLXUUlBTQH5NvkO8K09ypuYMFlvflWQMGgPJgcnMDpnN/ND5pASlqNncBMRul7x9+DxPbc6j8goL3lXJk3jo+nimBo6Mz5TVZiftyW1UtGW4vHp3qrPoxZde0K9ESklRXRHZ5dkcqThCTnlOr3HRK/HUezLVZyqR3pFEeUcR4R3BFM8phHuF42vyHb6GK3rFbDVzvv485xvOU1JfwrnacxTVFXG29izlzeV9X6CNQLdAUoJSnCG5RP9E9NqhqwuqGNvUNrey4R+neXVfUZdFU51G8NWZk7l/6TRiQ4Z34fTNg8X85J1jgCN+vu/hZei0jlmgEnQXqLXUcrzyOHlVec5bUV1Rv+tB+hp9ifCKINw7nHCvcMI8w5jsMZlQz1AmuU9SwjAI7NJOlbmK0oZSShtLKW0opaS+hJK6Eorri7nYeNFRNagfBLsFE+sf65yFJQYkEuYZpmZgCgorGvivLae7rVm6KnkS9y+dNiyhmOrGFlZs2OUcnf/wulh+cN1053El6APEbDWTX5PPicoTnK4+zZnqM5ypOUN9S/2Ar+lv8ifEPYQg9yAC3QIJMAUQ4BaAn9EPP5MfQW6Ox72N3l+KuKyUkmZrM5ebL1PRXEGVuYpqczWV5koqmyupMldR3lROeVM5Fc0VXTzC+4NeoyfaJ5oYvxim+053rpWoGZWiL3JKanjy0zz2FV69tSY5zJt1qZGsmTkZD+Pgs2LMrTbu+ctB9uQ7XivE28j2BzO6ZNwoQR9CpJRUNFdQVFfEubpznKs912Wa314qbLBohRZfoy9+JofQ+xh88DH64G3wxsvghafBE0+9J+56dzz0Hrjr3HHTuWHSmRz/a00YtAb0Gv2wjDallFillVZbKxabBbPVTLOtmWZrM02tTY6btclZkq2hpYG6ljrqLHXUWGqobaml2lxNlbnKpTi2K2iEhjDPMOfMqD0kFukdSZhXGHqNmh0pBs7homr+Z0cB/zh56apjHgYtK5ImcWNKKItiggbkv37oXBWPvnesSxrlxjvmsCJpUpfzlKCPEHZpp7ypnJL6EorriimuL6asoYzSxlLKGsq4bL7c7xDOUGDQGNBpdOi1enRCh1ajRSscN43QIIRA0CH67WELu7Rjs9uwScfNarc6by32llF5L94GbyZ7TibUI5RQj1DCvcKJ8I5wirhBa+j7IgrFIDheWsvLu8/y96NlWLopSu1t0pE6NYB5UX6kxwT2WuDC3GrjSEkN72Vd4M1DXXfA/8fyWL6fOf2q5yhBHyNY7VZHaKGpgvLmciqbK52hhc5hhsrmSupbBx7WGW8YtUb8Tf4EuQXh7+aPv8nfGYpq/789TKWySxRjhZqmFt7JusD//6KIwoqeU111GsEUPzciAjwI8DDQarNjtUmqmlo4UlJDyxU/Cu4GLY+siueOtKhur6cEfRzSYmtxhiRqW2qptThudS11jjBGSwP1rfU0tzbTZG2isbURs9WM2Wam2dqM2Wqmxd4y4JizK2iEBqPWiEFrcIZ53HRuuOncnKEgD70HXnpHiMjL4IWP0ccZPvIz+RFgCsBN56YWIRXjFiklx0vr+CinlI+PlnGhZuBh1+sSQvjlmqReC0orQf8SY5d2WmwOYW+1t2K1W50hFLvdjk3akG3/kIBw5N8LBEKILqEZvVaPTqNDJ3QYtAZVTUehuAIpJQUVDRw6V83+wkr2F1Z18WLvjuhAD+ZH+7MqOdSZa94bY7JItGJk0AiNMoVSKEYIIQQxwV7EBHtxS2oE4CiJV1zVRFFlE/VmK3qtwKDVYNJrSZrsTbD30H0/laArFArFMOJu0BE/yZv4ScNfIWniJzorFArFlwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwQl6AqFQjFBUIKuUCgUEwSXBF0Icb0Q4pQQIl8I8XA3x4UQ4pm240eFELOHvqkKhUKh6I0+BV0IoQX+CKwCEoFbhRCJV5y2CpjedrsP+J8hbqdCoVAo+sCVEXoqkC+lLJRStgBvAGuuOGcN8Kp0sB/wFUKEDnFbFQqFQtELrhSJDgNKOt0/D8x34ZwwoKzzSUKI+3CM4AEsQojcfrV2YhMIXB7tRowRVF90RfVHV77s/RHZ0wFXBF1085gcwDlIKTcCGwGEEIeklHNdeP0vBao/OlB90RXVH11R/dEzroRczgPhne5PAUoHcI5CoVAohhFXBP0gMF0IMVUIYQBuAT684pwPgTvbsl0WALVSyrIrL6RQKBSK4aPPkIuU0iqEWA98BmiBTVLK40KI+9uOvwB8AqwG8oEm4C4XXnvjgFs9MVH90YHqi66o/uiK6o8eEFJeFepWKBQKxThE7RRVKBSKCYISdIVCoZggDLugK9uADlzoi9va+uCoEGKvECJlNNo5UvTVH53OmyeEsAkhvjGS7RtpXOkPIUSGEOKIEOK4EGLnSLdxpHDhu+IjhPhICJHT1heurNtNfKSUw3bDsYhaAEQDBiAHSLzinNXApzhy2RcAXwxnm0br5mJfLAT82v5eNVH7wtX+6HTeNhwL798Y7XaP8ufDFzgBRLTdDx7tdo9iXzwKPNX2dxBQBRhGu+2jfRvuEbqyDeigz76QUu6VUla33d2PI59/ouLKZwPg+8A7QPlINm4UcKU/1gHvSimLAaSUE7VPXOkLCXgJIQTgiUPQrSPbzLHHcAt6T5YA/T1nItDf93kPjpnLRKXP/hBChAH/Arwwgu0aLVz5fMQCfkKIHUKIw0KIO0esdSOLK33xHJCAYwPjMeAHUkr7yDRv7OLK1v/BMGS2ARMAl9+nEOJaHIK+aFhbNLq40h8bgJ9IKW2OgdiExpX+0AFzgEzADdgnhNgvpTw93I0bYVzpi5XAEWAZMA3YIoT4p5SybpjbNqYZbkFXtgEduPQ+hRDXAH8CVkkpK0eobaOBK/0xF3ijTcwDgdVCCKuU8v0RaeHI4up35bKUshFoFELsAlKAiSborvTFXcCT0hFEzxdCnAXigQMj08SxyXCHXJRtQAd99oUQIgJ4F7hjAo66rqTP/pBSTpVSRkkpo4C3ge9OUDEH174rHwCLhRA6IYQ7DtfTkyPczpHAlb4oxjFTQQgRAsQBhSPayjHIsI7Q5fDZBow7XOyLx4EA4Pm2UalVTlBXORf740uDK/0hpTwphNgMHAXswJ+klBPOgtrFz8avgFeEEMdwhGh+IqX8MlvqAmrrv0KhUEwY1E5RhUKhmCAoQVcoFIoJghJ0hUKhmCAoQVcoFIoJghJ0hUKhmCAoQVcoFIoJghJ0hUKhmCD8H1BCkqunfjj9AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "filenames": { "image/png": "/Users/kpmurphy/github/ssm-book/_build/jupyter_execute/chapters/ssm/hmm_17_0.png" }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\n", "# Let's plot the observed data in 2d\n", "xmin, xmax = 0, 1\n", "ymin, ymax = 0, 1.2\n", "colors = [\"tab:green\", \"tab:blue\", \"tab:red\"]\n", "\n", "def plot_2dhmm(hmm, samples_obs, samples_state, colors, ax, xmin, xmax, ymin, ymax, step=1e-2):\n", " obs_dist = hmm.obs_dist\n", " color_sample = [colors[i] for i in samples_state]\n", "\n", " xs = jnp.arange(xmin, xmax, step)\n", " ys = jnp.arange(ymin, ymax, step)\n", "\n", " v_prob = vmap(lambda x, y: obs_dist.prob(jnp.array([x, y])), in_axes=(None, 0))\n", " z = vmap(v_prob, in_axes=(0, None))(xs, ys)\n", "\n", " grid = np.mgrid[xmin:xmax:step, ymin:ymax:step]\n", "\n", " for k, color in enumerate(colors):\n", " ax.contour(*grid, z[:, :, k], levels=[1], colors=color, linewidths=3)\n", " ax.text(*(obs_dist.mean()[k] + 0.13), f\"$k$={k + 1}\", fontsize=13, horizontalalignment=\"right\")\n", "\n", " ax.plot(*samples_obs.T, c=\"black\", alpha=0.3, zorder=1)\n", " ax.scatter(*samples_obs.T, c=color_sample, s=30, zorder=2, alpha=0.8)\n", "\n", " return ax, color_sample\n", "\n", "\n", "fig, ax = plt.subplots()\n", "_, color_sample = plot_2dhmm(hmm, samples_obs, samples_state, colors, ax, xmin, xmax, ymin, ymax)\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtzElEQVR4nO3deXxU5dn/8c81S5JJ2CGEJayKFRcQSUFFRcUFrRSteyvVtkr10bq0j621ffSxfbraX2tbtUKVqq2gVUBBUbSKIipCgiDIGlkDgYQkkD2zXb8/ZoiTzIRMYCBwuN6vV15k7rPc932W75y555AjqooxxhjncrV3A4wxxhxaFvTGGONwFvTGGONwFvTGGONwFvTGGONwnvZuQCI9evTQgQMHtnczjDHmqFFQULBbVbMTTTsig37gwIHk5+e3dzOMMeaoISJbWppmQzfGGONwFvTGGONwFvTGGONwFvTGGONwjgn6cF0dtQUFNGzYQOzf79FQiLrPPqPus8/QUKjJMg2FhdQWFBCuqzvczW1RYFcJtUuXEigpae+mHFbBiopIv7dvP+h1NXzxRWS/1tY2KQ8UF1O7dCnBsrKk1hP2+6ld9in1a9c2PabCYeo+/5zaTz9FA4Gk1hUoie7XXYduv/q3baM2P5/Q3r0HtZ7Gc2blyrhz5kgU2Lkzsl9LS5uUhyorqc3Px7+l6XeUum+/rlnDsfK3vlq960ZE+gHPAb2AMDBVVf/cbB4B/gxcBtQCN6vqsui08dFpbuApVf1tSnsAVPz7JXb95jeI242GQnhz+9J/yhQCO3dS9IMfoPUNkXamp5P7lz/j7dOHbd+/DX9RUeMyOfffT9frrk1105Kmfj/b7/8p1e+8g3i9qN9Px0svpc///RLxetutXYeaqrLrN79lz4svIGnpqN9P5qhR5D76J1xZWW1aV2DXLrZN/j7+LVsQjwcNBul53310ufoqtv/ov6n54AMkLQ1taKDzFRPp9dBDiNudcF17589n5wM/a2yjp3s3+k2ZgjY0sO32/yJUWYmIgNtN30d+T4exYxP3LxBgxwM/o2r+/Ejdfj8dL76YPr/5dcr2a6iykqI7f0DdihWRYycQoNvNN5N9z92RNrZBzZIlbL/7HtTvB0AyMsh97K9kjhiRkramUtjvZ8d9P6Z6wQIkPR1taKDThMvp9fDDlD05hbK//71xe2Sccgr9Hn+MmqVLKf7pA6Aa2a9du9Lvyb+RPmRIe3fnkJLW3tFEpDfQW1WXiUhHoAC4QlVXx8xzGfADIkE/Gvizqo4WETewHrgIKAKWAjfELptIXl6eJnt7Zd3y5Wy5+Ttoff2XhS4X3n65BEt3o82u6vD5SMvJwb91K4TDX/YzI4P+06aReXr7HNC7fvd7KqZPRxsamrSp23e/Q8+77mqXNh0O5dNnUPLII2jMpypJS6PjxRfT9w+PtGldGydeQUNhIcRchYovg6wzz6Lmww/jtm2PO++gxy23xK2nYeNGNn3jqqbHlAju7t3RQIBwsytmycjguHmv4+3TJ25dJX/6E+XPPtdkXZKRQbdJN9LzRz9qU/9asu2226n+8EOI+XQhPh+9f/EwnSdMSHo9wYoKCsddGHfOuLKyOP69Bbg7dkxJe1Nl569+xZ6XXmq8kIPItu1wwflUL3ivyTGF14vv1FOpX7266X4F3N26MeT99476CyoRKVDVvETTWh26UdXifVfnqloFrAH6NpttIvCcRiwGukTfIEYBhaq6UVX9wAvReVOm/J//anICAxAOEyje2eSE/3JaKDI8EBPyANrQQPk//5nKpiVNVal48cW4fmh9PRXPT2+XNh0u5c880/SEJPLppuqttwg3OyH3p2HDhsibd7N9rnX1VL/3XuJt+1zi/b3nxX/HD8moEq6qij/WiAx17Jk1O+G6KqbPiAsWra+nYsYLrXUpKaE9e6hpFvIAWldH2dPT2rSuynnz4s4LiAxVVc2ff1DtTDVVZc9LLzcJeYhs26r5b8UdUwQC1C1fjgaD8etqaKDmo48OZXPbXZv+w5SIDARGAJ80m9QX2Bbzuihalqh8dAvrngxMBujfv3/SbQruLoUEn0rC4TASCtH8g6s2+Al7PPHvcKrsWreOT+fOBcDr9TJ+/Pik23FQVOMPzKhQVRVzo2067O06DMKVlYnLQyHeeOUVwtHhm9b6HayoiAzDJZim4XDccQDgr6hIuG2Du3cnvEgIhUKgGn/sBAIUFhSwOMGxE66pSdjeUE0Nc+fMgejQyoHu11B1NXg8cUEPUL2zuE3HTqhiT8I3snBDAys/+og96ekH1daUCgYTthUib7yJ9rcCkiDoA4EA+e8uoCp6LB4R/UuxpINeRDoAM4F7VLX52dnidm2hPL5QdSowFSJDN8m2q8P551O34rO4qyaXSGTcrlm5ZGTgCoXiGiEZGQy8+mpGRj/qxp4gh5q4XGScfDL1q1bFTcscMYIJMR+/D2e7DofM0aOpevvtuCtJb04Ol113XeMYc2v9zjjp5IRXa5Kejjsri1B5ebMJQsczRifcth3GnkvVu+/Gvfm6RBCXK+4LSsnM5JQbv0XHcePi2uobNoy65cvj2uUbdioTvv71uLrbytu7Ny6fj1DzCwW3mx4XXsSwNhw7WaNHUfb00/H9Tk8n7+ab8Q0bdlBtTSXxekk/4QQa1q2Lm+bNyYl84d7seHBF36ia988jwlm33kpabmSg4kjoX6olddeNiHiJhPzzqjorwSxFQL+Y17nAjv2Up0zXa67Bk5ODRHciRMYne9x9F1lnn434fE3KO5x5Jj3uvadpeXo6np496XJt+30Z2+vB/4m0ad+Xg243kplJr5890G5tOhx6/vDeyJeu0fFRFUEyMuj98MNt+iLR3SGLnj+8F8nIaCyT9HQ82T3o/ZtfR8r3bVuPB1dmJjk//nHCdXUaP560QYOarsvno9uNN9J54sS4Yypj6FA6nHdewnXl/PznSGZmY93qckX268//J+m+7Y+43fT63/+NtHXf9vJ6cXfsSPadd7RpXb68PLJGjYo/Z849pzHkjyS9Hnqo6Tnj8eDKyqL3b3+Du3NnJC0tUh49pno9/DDpxx8ft1+7XHtNY8g7VTJ33QjwNLBGVf/YwmxzgDtF5AUiQzN7VbVYREqBISIyCNgOXA98MzVNj3BlZTFo5stUTJ9B1Tv/wdO1K10nTaLDmDFoKMTeuXPZO3MWitLlG1fR+esTELcb34knUv7PfxIsL6fjuHF0/eY3cXdo210eqeQbNoxBs2ZS/vQ0di5eTK8xZ9H9e98jrQ3DWEejtAEDGDx3DuX/+Ae1yz6lxOth+IMPknHiiW1eV7dvf5v0E06g/NnnKP2ikP5XfoOuN34Ld8eODHr5JcqenkbDhg34hg+j23e+2+LJLWlpDJz+PBX/fomqefPYXVvLV+6+iw4XXABA1pizGr9T6TRhAl2uuqrFu3d8p5zM4NmzKHv6aepXr2FXVianP/wwaSn8o32dLr4Ib+9/UvaPaZSuXEWfSy6h+03fxpOd8O9btUhEyH38MfbOmcueWbMoKy/n+O9PptPll6esramUefoIBs18mbKnnmbXkiX0Oudsun/3u6T168fgOa9S/uxz1Hz8Md6+fej+ne/gGz6cTpdczJ6XXqbytdeQzEy63nA9HS+8sL27cuhp9Dajln6As4kMt3wGLI/+XAbcBtwWnUeAx4EvgJVAXszylxG58+YL4Get1aeqjBw5UtvbnDlzjsi627Ndh0NL/Wtrvw9kOx2OulNVRyrrbuv8R+IxmMo2HYn9SwaQry1kaqtX9Kq6iMRj7bHzKJDwc6KqzgPmtVaPMcaYQ8Mx/zPWGGNMYhb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcBb0xhjjcMk8SnAacDlQoqqnJJh+H/CtmPUNBbJVtVxENgNVQAgIqmpeqhpujDEmOclc0T8DjG9poqo+oqqnqeppwE+B91W1PGaW86PTLeSNMaYdtBr0qroQKG9tvqgbgBkH1SJjjDEplbIxehHJJHLlPzOmWIG3RKRARCa3svxkEckXkfzS0tJUNcsYY455qfwydgLwYbNhmzGqejpwKXCHiJzb0sKqOlVV81Q1Lzs7O4XNMsaYY1sqg/56mg3bqOqO6L8lwGxgVArrM8YYk4SUBL2IdAbGAq/GlGWJSMd9vwMXA6tSUZ8xxpjkJXN75QzgPKCHiBQBDwFeAFV9MjrblcBbqloTs2gOMFtE9tUzXVXfTF3TjTHGJKPVoFfVG5KY5xkit2HGlm0Ehh9ow4wxxqSG/c9YY4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxOAt6Y4xxuFaDXkSmiUiJiCR8DKCInCcie0VkefTnwZhp40VknYgUisj9qWy4McaY5CRzRf8MML6VeT5Q1dOiP78AEBE38DhwKXAScIOInHQwjTXGGNN2rQa9qi4Eyg9g3aOAQlXdqKp+4AVg4gGsxxhjzEFI1Rj9mSKyQkTeEJGTo2V9gW0x8xRFyxISkckiki8i+aWlpSlqljHGmFQE/TJggKoOB/4KvBItlwTzaksrUdWpqpqnqnnZ2dkpaJYxxhhIQdCraqWqVkd/nwd4RaQHkSv4fjGz5gI7DrY+Y4wxbXPQQS8ivUREor+Piq6zDFgKDBGRQSKSBlwPzDnY+owxxrSNp7UZRGQGcB7QQ0SKgIcAL4CqPglcDdwuIkGgDrheVRUIisidwHzADUxT1c8PSS+MMca0qNWgV9UbWpn+GPBYC9PmAfMOrGnGGGNSwf5nrDHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOJwFvTHGOFyrQS8i00SkRERWtTD9WyLyWfTnIxEZHjNts4isFJHlIpKfyoYbY4xJTjJX9M8A4/czfRMwVlWHAb8Epjabfr6qnqaqeQfWRGOMMQcjmSdMLRSRgfuZ/lHMy8VEHgJujDHmCJHqMfrvAW/EvFbgLREpEJHJ+1tQRCaLSL6I5JeWlqa4WcYYc+xq9Yo+WSJyPpGgPzumeIyq7hCRnsDbIrJWVRcmWl5VpxId9snLy9NUtcsYY451KbmiF5FhwFPARFUt21euqjui/5YAs4FRqajPGGNM8g466EWkPzALmKSq62PKs0Sk477fgYuBhHfuGGOMOXRaHboRkRnAeUAPESkCHgK8AKr6JPAg0B14QkQAgtE7bHKA2dEyDzBdVd88BH0wxhizH8ncdXNDK9NvAW5JUL4RGB6/hDHGmMPJ/mesMcY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4nAW9McY4XMqeGesU4bDy4Re7WVQMfbZUcHr/LkQfnsKq7XtZvm0PfbpkcO6QbDzuA3ufbAiGWLC2lLKaBkYN7MaQnI4H3N6de+t5f30JaR4X44bm0CnDC0BVfYB315ZQ5w9x7gnZ9OniAyAQCvP+ulJ2VdUzol9XTurTCQBVJX9LBet2VjGoRxZnDu6OyxXpd2FJNUs2ldM108v5J/Ykw+sGoLSqgQXrSnCLMG5oT7pkpgFQ6w/yzpoSqhuCnH18D/p1ywQgFFYWbihle0Udp/btzPB+XQ643xU1ft5ZW0JYlQtO7EmPDukHvK7Pd0T269ZyuDQUPuD9eiDW76pi6eZyNpXBJcEwaZ5I3SWV9by3rhSPWxg3NIfOvsh+rW4I8s6aXdT5Q5xzQjZ9o/u1JarKsq0VrCmuongPfC2sjft1Y2k1izeW0yXTywUx+7Uldf4Q76zdxYc7YVhZDQO6ZwGR/bqocDfbyms5uU8nTuv35TnzWdEePivaS9+uPs4dko07Wve28loWFe6mQ7qHcUN7kpnW/lFUWR/g3TUlLNkFX91bT6/OGUDknHlvXSklVfWMHNCVE3t9ec4s3VzB+l1VDM6OnDP7+l1YUsUnm8rpnpXO+Sdmk+6JbNuSqnreWxvdryfm0Dkzsl9rGoK8s7aEmmbnTKok84SpacDlQImqnpJgugB/Bi4DaoGbVXVZdNr46DQ38JSq/jaFbU+50qoGrp3yESVVDfj9wtynP+HkPp2YdvNX+cGMT/lkYzmK4nG56Jjh4aXbziS3a9t2yPpdVVw35WP8ISUUDoPCZcN684erhzeegMmauvAL/t9b63G5BBfwwKxVPHHj6XhdLib/Mx8BwgphVe4aN4TLh/Xm2ikfU90QJBSOPH/93CHZPHL1MG76x1LW76oirIpbhN5dfLw4+Qwemb+OVz7dDgJul+B1uZh+6xks31bBw3NX43YJAvzsFeUP1wynd+cMbpq2FAXC4TBhhZvPGsjNYwZyzZMfU1HrJxRWBGHkgK48fXNe40mQrDnLt3Pfy581hsb/vLKKX0w8meu+2r9N6wmEwtz+rwI+LCxDUTQkvPK7d3np+2fRv3tqT7TmwmHl3n8vZ/7nOwHQsDD7N+/w4vfPYOH63fzuzbW4RHAJPDB7JY/dcDqZ6W5uebbpfr197HHcc9EJCeuo9QeZ9PQS1hRXEtZI/9740/u8OPkM/vyfDbxUUBSpwwUel4vnbxnNKX07J1zXp1sr+Pa0JYRV8QeE2X9ayI1nDGDyOYO4ZspiymoaGvfr8H6dmTopjzumLyN/cwWK4nYJXXxpvHTbmTz38Rb+8eEmXAIulyAIz3znq+QN7HaoNner3ltXwu3/WoZLIBAUZj6ygHsvOoGLT8rhuimLqQsECUbPmXEn5vCrK09h0tNL+KK0uvGcye2ayYzJo/nV62t4fWUxAB6X4HW7eGHymXyyqYxfvb6m8bh9QFfy5+tOo0tmGt99dikQOS7CCreeM4j/vuTElPUvmbfRZ4DHgOdamH4pMCT6Mxr4GzBaRNzA48BFQBGwVETmqOrqg230ofLjl1ewrbwuukOFgD/EZ0V7ueXZfFYU7aE+EI7OGabWH+SuGZ8y67/GJL1+VeXW5/KpqA00KX9j5U7OPr4H3zg9N+l1fb5jL398ez0NwXCT8tv/WYDLJdT6Q03KH3u3kJcLiiitaiB6vALwwYZSbn0un9U7KvGHvlzXlrIabnl2Ket2VlPfpI4QN0/7hMr6YFzdP3pxOeleN9UNwSblzy3ewvvrSyneW9/4BgOQv7mcJ9/byN0XDkm63yVV9dz38mdxdT/06uecdVzbroSe/Wgziwp3x+xXobSqgTumL2PuD85Oej0HYuayIt76fFeTuv01fm6atpSymoa4/t0xvQCP2xW3X6cs3Mg5J/Rg5ID4kPzjW+tZtX1vzLqErWW13PJsPhtKqpvVEeK7zyxl8U/HxV1wBENhvvdsPlX1wcb1EA4zY8lWPt5YxvaKWkIxx9TyrXu49bl8Vmzb0+TYqQ/U851nlrK1rDauf997Np+lP7uw8RPN4VTdEOS/nl9GXWDfto3079H/rGf6J1vZXdOAxvTv3bUl3PpcPut2VuKP6fjG0mpueTafNTurYvYrQIib/rGEPTX+uH7f9cKneN1uahqa7tenF23m7CHZnDG4e0r62OpWVdWFQPl+ZpkIPKcRi4EuItIbGAUUqupGVfUDL0TnPSLVB0J8sGF347v2Pg3BMEs3lzfbcZErqlXbKymrbki6jg0l1ZRWxc9fFwjx/Cdb29TeVz7djr/ZQQNEr6Q1rrw+EGJrWQ3NJ9UFwuRvqWgS8gCBkPLptr3UBpoegAAVtQECCepQiDuQIfKRf92uqiYhD1AfDPNiftv6/eaqnUiCDz4h1carqGTNWLI14X5dv6uKksr6Nq2rrZ7/ZGtMsEQosHNvXQv7VQiGEuzXYIh/5xclrGPmsqK4/REMKyuK9sa9YUBk+OCz7XvjyvO3VOAPxs9f6w+xZkclzZtVv++caVZ3KKys31kV1+/ItDCfbCpL2I9DbcHakoTHVEMgzPY9dU1CHiLna8GWiiYhDxAIK59u3UNdgm27u6qBQCh+v4IkLK8PhHgpf1tburFfqRgY6wvEtqgoWpaofHRLKxGRycBkgP792/YRPBXCqsSfRhGqkSv8uGXCIea9+RZdosPDXq+X8ePHt1iHPxhOeEABlOwuZ+7cuY3raU19IBx3AAIEQ6FoS5tWpPvrRwvlkaUS9TtMKEEdwVCYyMVggnW1UEdVTV1jv6H1vvuDYcIJzpdgKMxnq1Yzt3J1UusB4t7c9gmHQ8yb/zbdM5Jr04FIFOYQ2UeRKU23VTgUc7XZZH74YtNW5s7dAjRta6I3huhScesBCIWCLHj/A7ZFR2/2rSsSRIkPXG1hXS0day3VHQwGWfTRYvas/bLu/Z1LqeQPRoZQm1MUDbftnGlpe4TD4YT7dX/na32CN8QDlYqgT35vJtyc0QmqU4GpAHl5eS3Od6hkpnk4tW9nVmzb06SRHpcwOLsDm3fXxL2D9+/egRuvGtv4BUxsYCUytHcn0j3xH9MyvC5uPu8rTDhncNLtvfTUXsxcVhR3ZeZ2uQCBZiHm87rITPNQVuNvUp7ucdG/WyYbS6ubXJmJwICuWeyqqqeu2VWvL91LMByOuxr2uF24RAg2C7HMNDedfV6K9za9Sva6hW/kDWDChLivflo0bmgOf5i/Lq483evmzivHNH65nIwJw/rw9KJNcVe9vbtkcdPV5zXu10PhihF92Li7Om4bdvR58Qc17qrX5XIhLWzb2y4bwYUn5cTVceFJOcxdsaPJp1QB+nXLpLTKH1eHx+Pljusvihs+yRvQLRpsTWWmuemamcb2PXVNyr0uYVB2BzbtriHQ7JzJ6eSjsj4Q/4nC5eau6y8iKz0SSa2dS6l07gnZcZ/kAXxeNxled9xQa7rHRb9uPjaWNv2E7BIY0C2L4sr6uP2ale4hGCZum7v3s18nDO97kD2LaVsK1lEE9It5nQvs2E/5EeuRq4fRMcNDhjeyWTLT3PTsmM7USXnkdsskMy3ypWGGx0VWuoc/XXdam8LA7RL+cv0IfF5348mUmebmxF4d+dYZA9rU1jMHd+eyU3vjS3MjgAslw+vip5cN5aEJJ5HhceF2RU7szDQ3Fw7NYcqkkWSluUmPqXtg9yymTBpJj47pjf3zed10zvAy9aY8TunbubHc6xYyvC7++s0RfHNUf3xed+QLNYm8Wd15wRB+f/UwMjwuPNFx3sw0N2cM7s7USXl0SPeQEVN3ny4+7m3hi8SWDOqRxW1jjyPD64p+elB8XjeTzhjQppAHuP284+gfs1/TPS6y0tw8en3b9uuB+PaZAzmhZ8fGutM8LnxeN09883QmntYHn9eNxGzb+8Z/hV9MPCVu2449IZsLTuyZsI6fXnYi2c32ayeflymT8hieG79fH73utIRj5L40N3+4ZhgZXhde95d15w3oytRJI+POmZzOGUydlEffLr4m50yHdA9TJo3kzMHdG8s9rkjdv/3GsMaQP9yyO6bzwGVDyfA2PWcuPaU3T944ksxm58xx2R2YMmkk3Tuk44veqeSLXsxM/XYeQ3t3+nK/uqP79caRXHV6btx+vfeiE/jNN04lw9t0v559fA8uTvDmfaBEE33+bz6TyEDgtRbuuvkacCeRu25GA39R1VEi4gHWA+OA7cBS4Juq+nlr9eXl5Wl+fn5b+pEye2sDzFxWxBel1QzP7cKE4X3wpblpCIZ4Y+VOlm4up3+3TK4amRt3S9/cuXOZMGFCq3UU763jpfwilqxcx40XjuTCoTlJ39IXW4eq8smmct5ctZOiLZu475qxfKVX5FbNwpIqZi3bTo0/yCUn9eLM4yK3fpVU1TOzYDvbK2oZNbg740/uRZrHRa0/yJzlO/hs+15O6NmBK0/PpbPPSyis/GfNLhZt2E12hzSuyuvXeEvfsq0VvLaiGLcLJp7Wt/GOjc27a3i5oIjK+gDjhuZwzvE9cLmE8ho/MwuK2FRWQ96Arlx2au9Wb+lraduu2r6XV5dvZ0PhRu6+cgwj+ndNaj3N+YNh3lhVzJJN5eR2zeTqkblkdzzwWzUTtbWlacFQmLdW7+Kjwt3kdM7gmpH96NU5o/G2vTdXFZPmcTHxtL4M7R15E/uitJpZBUVU+4NcfFIvzjqu+37flOr8Ieas2M6Kor0Mye7AN07PpXNmZL8uWFvC++tL6d4hjatH5rZ6B9nWslpeLthGRW2AC4b2ZOyQbFwuYU+tn5cLiti0u4bT+kXOmQyvm/pAiHkri8nfUsHA7plcPbIf3bLSCIeVDwp3886aXXTK8HL1yFwG9shKehu2dZsna/2uKmYtK6I+EOaSk3txxuBukXOmsp6XC4rYsaeOM47rziUn98LrdlHTEOSV5dv5fEclJ+Z05MrT+9Ixw0swFG48Z3p2yuCavFx6d/ahqhRsqWDeymI8bhdXnNa38eJkY2k1M5dtp6o+wEUn5XD28T3afLEhIgWqmpdwoqru9weYARQDASJX6d8DbgNui04XInfXfAGsBPJilr2MSNh/Afystbr2/YwcOVKPRnPmzDmk8+9vmQNZ19HkaOr3/tp0JLb3SHQ4ziWnAfK1hUxt9bOSqt7QynQF7mhh2jxgXmt1GGOMOXTsTyAYY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDJRX0IjJeRNaJSKGI3J9g+n0isjz6s0pEQiLSLTpts4isjE5rn+cDGmPMMazVJ0yJiJvIowIvIvIowaUiMkdVV++bR1UfAR6Jzj8BuFdVy2NWc76q7k5py40xxiQlmSv6UUChqm5UVT/wAjBxP/PfQOQ5s8YYY44AyQR9X2BbzOuiaFkcEckExgMzY4oVeEtECkRkckuViMhkEckXkfzS0tIkmmWMMSYZyQS9JCjTFuadAHzYbNhmjKqeDlwK3CEi5yZaUFWnqmqequZlZ2cn0SxjjDHJSCboi4B+Ma9zgR0tzHs9zYZtVHVH9N8SYDaRoSBjjDGHSTJBvxQYIiKDRCSNSJjPaT6TiHQGxgKvxpRliUjHfb8DFwOrUtFwY4wxyWn1rhtVDYrIncB8wA1MU9XPReS26PQno7NeCbylqjUxi+cAs0VkX13TVfXNVHbAGGPM/rUa9ACqOg+Y16zsyWavnwGeaVa2ERh+UC00xhhzUOx/xhpjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMNZ0BtjjMMlFfQiMl5E1olIoYjcn2D6eSKyV0SWR38eTHZZY4wxh1arDx4RETfwOHARkefHLhWROaq6utmsH6jq5Qe4rDHGmEMkmSv6UUChqm5UVT/wAjAxyfUfzLLGGGNSIJmg7wtsi3ldFC1r7kwRWSEib4jIyW1cFhGZLCL5IpJfWlqaRLOMMcYkI5mglwRl2uz1MmCAqg4H/gq80oZlI4WqU1U1T1XzsrOzk2iWMcaYZCQT9EVAv5jXucCO2BlUtVJVq6O/zwO8ItIjmWWNMcYcWskE/VJgiIgMEpE04HpgTuwMItJLRCT6+6joesuSWdYYY8yh1epdN6oaFJE7gfmAG5imqp+LyG3R6U8CVwO3i0gQqAOuV1UFEi57iPpijDEmgVaDHhqHY+Y1K3sy5vfHgMeSXdYYY8zhY/8z1hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHC6poBeR8SKyTkQKReT+BNO/JSKfRX8+EpHhMdM2i8hKEVkuIvmpbLwxxpjWtfqEKRFxA48DFxF52PdSEZmjqqtjZtsEjFXVChG5FJgKjI6Zfr6q7k5hu40xxiQpmSv6UUChqm5UVT/wAjAxdgZV/UhVK6IvFwO5qW2mMcaYA5VM0PcFtsW8LoqWteR7wBsxrxV4S0QKRGRySwuJyGQRyReR/NLS0iSaZYwxJhnJPBxcEpRpwhlFzicS9GfHFI9R1R0i0hN4W0TWqurCuBWqTiUy5ENeXl7C9RtjjGm7ZK7oi4B+Ma9zgR3NZxKRYcBTwERVLdtXrqo7ov+WALOJDAUZY4w5TJIJ+qXAEBEZJCJpwPXAnNgZRKQ/MAuYpKrrY8qzRKTjvt+Bi4FVqWq8McaY1rU6dKOqQRG5E5gPuIFpqvq5iNwWnf4k8CDQHXhCRACCqpoH5ACzo2UeYLqqvnlIemKMMSahZMboUdV5wLxmZU/G/H4LcEuC5TYCw5uXG2OMOXzsf8YaY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDWdAbY4zDJfX36E3LVJWPd3zMgm0LKK4t5pS9pzCo8yAA1pav5fWNr+MP+bl44MWc3vN0og9haZNAOMA7W99hSfESKmorGF07mp6ZPVPdlSNOWMMs2r6IhUUL2VW7i+GVw+nfqX97N6tFO2t28mrhqyypXoJ3s5dx/cfhdXnbu1lHldVlq3l94+sU1hSSW5LLiJ4jANhRvYNXC19ld91uzuxzJuf1Ow+Py+IrWUltKREZD/yZyBOmnlLV3zabLtHplwG1wM2quiyZZY9mYQ1zz4J7WFy8mLpgHS5cLJ67mJ+P/jll9WX8bcXfCIQCKMrswtlcPvhyHjzzwTbVURuo5aY3bmJr1VZqg7V48PDh7A95fNzjfLXXVw9Rz9pfMBzkznfuZFnJssZt+/Gcj/m/Mf/H+EHj27t5cT4p/oQ737mTsIbxh/2s+nAV01ZO49lLn8Xn8bV3844KU1ZM4amVT+EP+VGUZW8t48ohV3JO33P44Xs/JKQhAuEAr218jSFdhzDtkmmkudPau9lHhVaHbkTEDTwOXAqcBNwgIic1m+1SYEj0ZzLwtzYse9RasG1BY8gDhAnTEGrgl4t/yRPLn6Ah1ECYMIpSF6xj7hdzWV6yvE11PL/meTZVbqI2WAtAkCB1wTp+svAnhDWc6i4dMeZvnt8Y8vDltn3wowcby44UYQ3zk4U/oT5Ujz/sB6A2WMvGvRuZsWZGO7fu6LCtaht/X/l36kP1jedMfaieWetncd/C+6gP1RMIB4DItl1Xvo6Z62e2c6uPHsmM0Y8CClV1o6r6gReAic3mmQg8pxGLgS4i0jvJZY9ab256s8XQUTSurCHUwH+2/qdNdby+6XUaQg1x5dWBajbt3dSmdR1N5m2cl3DbusXNsl3L2qFFLSvcU5iwrQ2hBl7f9Ho7tOjos7BoYcLyhnAD/pA/rrw+VM9rm1471M1yjGSGbvoC22JeFwGjk5inb5LLAiAik4l8GqB//yN3HDZWujsdQeJCPRxOfKUtCFu/2Mrc4rkAeL2tj9+muRJ/NA0EAyx6bxGr3auTXtfRJN2TnrA8EAiQ/0k+5d5y4Mjod5orrcVPV7WVtcydO7fx9ZHQ3iNRmjsNV4LrTkFaPJ+qK6obt61t1/1LJugTfXvY/HK1pXmSWTZSqDoVmAqQl5eXcJ4jzZVDrmT+5vnUh+qblKd70glqkGAo2KTc6/Zyz/h7GNx5cNJ1XPuVa/n90t83uWIUhAGdB3DTFTcdXAeOYFcNuYpF2xfFXSn70n3cccUdR9QXcQM6DaB3Vm82V25u8qbv8/iYPGoyE4ZMaMfWHR3G9R/H75b8Lq7c6/LSJb0LJXUlTcp9Hh+3n3X7Efl9zZEomaGbIqBfzOtcYEeS8ySz7FFrZM5Ibjr5JtLcafjcPjI9mWR5snhs3GP8/tzfk+HOINOTic/jI82dxo/yftSmkAe48vgrGZs7lnR3OhnuDLK8WXTL6Maj5z96aDp1hDirz1lc+5VrI/32ZJDlyaKDtwOPjXvsiAp5ABHh0QsepVtGN7K8WWS4M0h3p3NBvwuYeLxjRioPqW4Z3fj12b9ucs6ku9P5yaif8MSFT9AlvUuTbXvpwEu5ZOAl7d3so4ao7v/iWUQ8wHpgHLAdWAp8U1U/j5nna8CdRO66GQ38RVVHJbNsInl5eZqfn3/AnTrciquL+WjHR2R5szg391wyvZkAVPoreX/b+wTDQc7JPYcevh4HXMeGig0sL11Oti+bMX3HHDO37W2r2saS4iV0SOvA2NyxZHgy2rtJLQqEAizavojd9bsZkT2C47se395NOursbdjLwqKFhDTEOX3PobuvOwD+kJ8Pij6goqGCkTkjG29hNl8SkQJVzUs4rbWgj67gMuBRIrdITlPVX4nIbQCq+mT09srHgPFEbq/8jqrmt7Rsa/UdbUFvjDHt7aCD/nCzoDfGmLbZX9Dbn0AwxhiHs6A3xhiHs6A3xhiHs6A3xhiHOyK/jBWRUmDLAS7eA9idwuYcLazfxxbr97ElmX4PUNXsRBOOyKA/GCKS39I3z05m/T62WL+PLQfbbxu6McYYh7OgN8YYh3Ni0E9t7wa0E+v3scX6fWw5qH47bozeGGNMU068ojfGGBPDgt4YYxzOMUEvIuNFZJ2IFIrI/e3dnkNJRKaJSImIrIop6yYib4vIhui/XduzjakmIv1EZIGIrBGRz0Xk7mi50/udISJLRGRFtN8PR8sd3e99RMQtIp+KyGvR18dKvzeLyEoRWS4i+/4S8AH33RFB7/SHkCfwDJE/CR3rfuAdVR0CvBN97SRB4EeqOhQ4A7gjuo+d3u8G4AJVHQ6cBowXkTNwfr/3uRtYE/P6WOk3wPmqelrM/fMH3HdHBD0Ofwh5c6q6EChvVjwReDb6+7PAFYezTYeaqhar6rLo71VETv6+OL/fqqrV0Zfe6I/i8H4DiEgu8DXgqZhix/d7Pw64704J+pYeTn4syVHVYoiEItCzndtzyIjIQGAE8AnHQL+jwxfLgRLgbVU9JvpN5IFFPwZinw5+LPQbIm/mb4lIgYhMjpYdcN+PrIdvHrikH0Jujm4i0gGYCdyjqpWRh5s5m6qGgNNEpAswW0ROaecmHXIicjlQoqoFInJeOzenPYxR1R0i0hN4W0TWHszKnHJF7+iHkCdpl4j0Boj+W9LO7Uk5EfESCfnnVXVWtNjx/d5HVfcA7xH5fsbp/R4DfF1ENhMZir1ARP6F8/sNgKruiP5bAswmMjx9wH13StAvBYaIyCARSQOuB+a0c5sOtznATdHfbwJebce2pFz0ucRPA2tU9Y8xk5ze7+zolTwi4gMuBNbi8H6r6k9VNVdVBxI5n99V1RtxeL8BRCRLRDru+x24GFjFQfTdMf8z9kAeQn60EpEZwHlE/nTpLuAh4BXg30B/YCtwjao2/8L2qCUiZwMfACv5csz2ASLj9E7u9zAiX7y5iVyY/VtVfyEi3XFwv2NFh27+W1UvPxb6LSKDiVzFQ2R4fbqq/upg+u6YoDfGGJOYU4ZujDHGtMCC3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHM6C3hhjHO7/A2a5ycuhtDvHAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "filenames": { "image/png": "/Users/kpmurphy/github/ssm-book/_build/jupyter_execute/chapters/ssm/hmm_18_1.png" }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Let's plot the hidden state sequence\n", "\n", "fig, ax = plt.subplots()\n", "ax.step(range(n_samples), samples_state, where=\"post\", c=\"black\", linewidth=1, alpha=0.3)\n", "ax.scatter(range(n_samples), samples_state, c=color_sample, zorder=3)\n" ] } ], "metadata": { "interpreter": { "hash": "6407c60499271029b671b4ff687c4ed4626355c45fd34c44476827f4be42c4d7" }, "kernelspec": { "display_name": "Python 3.9.2 ('spyder-dev')", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }