jsl/demos/fixed_lag_smoother.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "fixed_lag_smoother.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/probml/JSL/blob/main/jsl/demos/fixed_lag_smoother.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Fixed Lag Smoother\n",
        "\n",
        "This notebook achieves the following: \n",
        "\n",
        "1. implements JAXified versions of fixed lag smoother;\n",
        "\n",
        "2. demonstrates the correctness of the implementations by comparison of their full-lag smoothed posteriors with the results of JSL's implementation of `hmm_forwards_backwards_jax`;\n",
        "\n",
        "3. demonstrates the improved \"online\" performance of a version that uses a vectorized approach to compute the $\\beta$ values across its sliding window against the one that iteratively smooths backwards inside the window, as implemented in Kevin Murphy's HMM Toolbox for Matlab. \n",
        "\n",
        "The JAX State-Space Models Library (JSL) is available at:\n",
        "https://github.com/probml/JSL\n",
        "\n",
        "Kevin Murphy's Hidden Markov Model (HMM) Toolbox for Matlab is available at:\n",
        "https://www.cs.ubc.ca/~murphyk/Software/HMM/hmm.html\n",
        "\n",
        "Author: Peter G. Chang ([@peterchang0414](https://github.com/peterchang0414))"
      ],
      "metadata": {
        "id": "YNVwtTv8P0bZ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 0. Imports"
      ],
      "metadata": {
        "id": "Fob7r0qaR1Nk"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install -q flax"
      ],
      "metadata": {
        "id": "vR05NbgRYED1",
        "outputId": "bac3616b-b00b-4020-dc46-db858eacfacf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\u001b[?25l\r\u001b[K     |█▉                              | 10 kB 30.0 MB/s eta 0:00:01\r\u001b[K     |███▋                            | 20 kB 32.1 MB/s eta 0:00:01\r\u001b[K     |█████▍                          | 30 kB 36.4 MB/s eta 0:00:01\r\u001b[K     |███████▏                        | 40 kB 37.7 MB/s eta 0:00:01\r\u001b[K     |█████████                       | 51 kB 27.0 MB/s eta 0:00:01\r\u001b[K     |██████████▊                     | 61 kB 30.1 MB/s eta 0:00:01\r\u001b[K     |████████████▌                   | 71 kB 28.4 MB/s eta 0:00:01\r\u001b[K     |██████████████▎                 | 81 kB 28.6 MB/s eta 0:00:01\r\u001b[K     |████████████████                | 92 kB 30.6 MB/s eta 0:00:01\r\u001b[K     |█████████████████▉              | 102 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |███████████████████▋            | 112 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |█████████████████████▍          | 122 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |███████████████████████▏        | 133 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |█████████████████████████       | 143 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |██████████████████████████▊     | 153 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |████████████████████████████▌   | 163 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |██████████████████████████████▎ | 174 kB 29.2 MB/s eta 0:00:01\r\u001b[K     |████████████████████████████████| 184 kB 29.2 MB/s \n",
            "\u001b[K     |████████████████████████████████| 140 kB 56.8 MB/s \n",
            "\u001b[K     |████████████████████████████████| 72 kB 780 kB/s \n",
            "\u001b[?25h"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!git clone https://github.com/probml/JSL.git"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "N4P2DsduVbW5",
        "outputId": "77d339b4-f0d2-4aa9-a903-36efcefedd1c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Cloning into 'JSL'...\n",
            "remote: Enumerating objects: 1781, done.\u001b[K\n",
            "remote: Counting objects: 100% (1726/1726), done.\u001b[K\n",
            "remote: Compressing objects: 100% (1138/1138), done.\u001b[K\n",
            "remote: Total 1781 (delta 1126), reused 1124 (delta 567), pack-reused 55\u001b[K\n",
            "Receiving objects: 100% (1781/1781), 5.78 MiB | 29.44 MiB/s, done.\n",
            "Resolving deltas: 100% (1145/1145), done.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import sys\n",
        "sys.path.insert(0,'/content/JSL')\n",
        "from jsl.hmm.hmm_lib import HMMJax, hmm_forwards_backwards_jax, normalize"
      ],
      "metadata": {
        "id": "2cHM_Qq4VqCh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from functools import partial\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import jax.numpy as jnp\n",
        "from jax import vmap\n",
        "from jax import jit\n",
        "from jax.random import PRNGKey, split, choice, uniform"
      ],
      "metadata": {
        "id": "PNu3vU_vXxky"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, a quick aside. The current JSL implementation of `normalize` does not work as expected when applied to a rectangular matrix with `axis=1`, as demonstrated by the error raised below:"
      ],
      "metadata": {
        "id": "JwQLLuagdSvu"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "A = jnp.array([[1, 2], [2, 3], [3, 4]])\n",
        "\n",
        "_ = normalize(A, axis=0) # Works as expected\n",
        "# _ = normalize(A, axis=1) # Raises \"Incompatible shapes for broadcasting error\""
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "U2BmOvwfdUql",
        "outputId": "c722afce-03e4-4c51-89bf-75c6c1ce1ccc"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "The following implementation fixes the issue by adding `keepdims=True`."
      ],
      "metadata": {
        "id": "KKswVJtgd8AS"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "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, keepdims=True)\n",
        "    c = jnp.where(c == 0, 1, c)\n",
        "    return u / c, c\n",
        "\n",
        "_ = normalize(A, axis=1)     # Works as expected"
      ],
      "metadata": {
        "id": "vn0UidP4eGNx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1. Fixed Lag Smoother - Two Implementations"
      ],
      "metadata": {
        "id": "jbmIOFafYk7b"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QDJ7YcaZGSTI"
      },
      "outputs": [],
      "source": [
        "# Naive (un-vectorized) version\n",
        "@partial(jit, static_argnums=(1))\n",
        "def fixed_lag_smoother_iterative(params, win_len, alpha_win, obs_seq_win, obs, act=None):\n",
        "    '''\n",
        "    Computes the smoothed posterior for each state in the lagged window of\n",
        "    fixed size, win_len.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    params      : HMMJax\n",
        "        Hidden Markov Model (with action-dependent transition)\n",
        "\n",
        "    win_len     : int\n",
        "        Desired window length (>= 2)\n",
        "    \n",
        "    alpha_win   : array\n",
        "        Alpha values for the most recent win_len steps, excluding current step\n",
        "    \n",
        "    obs_seq_win : array\n",
        "        Observations for the most recent win_len steps, excluding current step\n",
        "    \n",
        "    obs         : int\n",
        "        New observation for the current step\n",
        "    \n",
        "    act         : array\n",
        "        (optional) Actions for the most recent win_len steps, including current step\n",
        "    \n",
        "    Returns\n",
        "    -------\n",
        "    * array(win_len, n_states)\n",
        "        Updated alpha values\n",
        "    \n",
        "    * array(win_len)\n",
        "        Updated observations for the past d steps\n",
        "    \n",
        "    * array(win_len, n_states)\n",
        "        Smoothed posteriors for the past d steps\n",
        "    '''\n",
        "    if len(alpha_win.shape) < 2:\n",
        "        alpha_win = jnp.expand_dims(alpha_win, axis=0)\n",
        "    curr_len = alpha_win.shape[0]\n",
        "    win_len = min(win_len, curr_len+1)\n",
        "    assert win_len >= 2, \"Must keep a window of length at least 2.\"\n",
        "\n",
        "    trans_mat, obs_mat = params.trans_mat, params.obs_mat\n",
        "    n_states, n_obs = obs_mat.shape\n",
        "    \n",
        "    # If trans_mat is independent of action, adjust shape\n",
        "    if len(trans_mat.shape) < 3:\n",
        "        trans_mat = jnp.expand_dims(trans_mat, axis=0)\n",
        "        act = None\n",
        "    if act is None:\n",
        "        act = jnp.zeros(shape=(curr_len+1,), dtype=jnp.int8)\n",
        "\n",
        "    # Shift window forward by 1\n",
        "    if curr_len == win_len:\n",
        "        alpha_win = alpha_win[1:]\n",
        "        obs_seq_win = obs_seq_win[1:]\n",
        "    new_alpha, _ = normalize(\n",
        "        obs_mat[:, obs] * (alpha_win[-1][:, None] * trans_mat[act[-1]]).sum(axis=0)\n",
        "    )\n",
        "    alpha_win = jnp.concatenate((alpha_win, new_alpha[None, :]), axis=0)\n",
        "    obs_seq_win = jnp.append(obs_seq_win, obs)\n",
        "\n",
        "    # Smooth backwards inside the window\n",
        "    beta_win = jnp.ones(shape=(win_len, n_states))\n",
        "    gamma_win = jnp.array(alpha_win)\n",
        "    for t in range(win_len-2, -1, -1):\n",
        "        new_beta, _ = normalize(\n",
        "            (beta_win[t+1,:] * obs_mat[:, obs_seq_win[t+1]] *\n",
        "             trans_mat[act[t]]).sum(axis=1)\n",
        "        )\n",
        "        beta_win = beta_win.at[t, :].set(new_beta)\n",
        "\n",
        "        new_gamma, _ = normalize(alpha_win[t, :]*beta_win[t, :])\n",
        "        gamma_win = gamma_win.at[t, :].set(new_gamma)\n",
        "    return alpha_win, obs_seq_win, beta_win, gamma_win"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Vectorized version\n",
        "@partial(jit, static_argnums=(1))\n",
        "def fixed_lag_smoother_vectorized(params, win_len, alpha_win, bmatrix_win, obs, act=None):\n",
        "    '''\n",
        "    Computes the smoothed posterior for each state in the lagged window of\n",
        "    fixed size, win_len.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    params      : HMMJax\n",
        "        Hidden Markov Model (with action-dependent transition)\n",
        "    \n",
        "    win_len     : int\n",
        "        Desired window length (>= 2)\n",
        "    \n",
        "    alpha_win   : array\n",
        "        Alpha values for the most recent win_len steps, excluding current step\n",
        "    \n",
        "    bmatrix_win : array\n",
        "        Beta transformations for the most recent win_len steps, excluding current step\n",
        "    \n",
        "    obs         : int\n",
        "        New observation for the current step\n",
        "    \n",
        "    act         : array\n",
        "        (optional) Actions for the most recent win_len steps, including current step\n",
        "    \n",
        "    Returns\n",
        "    -------\n",
        "    * array(win_len, n_states)\n",
        "        Updated alpha values\n",
        "    \n",
        "    * array(win_len, n_states)\n",
        "        Updated beta transformations\n",
        "    \n",
        "    * array(win_len, n_states)\n",
        "        Smoothed posteriors for the past d steps\n",
        "    '''\n",
        "    if len(alpha_win.shape) < 2:\n",
        "        alpha_win = jnp.expand_dims(alpha_win, axis=0)\n",
        "    curr_len = alpha_win.shape[0]\n",
        "    win_len = min(win_len, curr_len+1)\n",
        "    assert win_len >= 2, \"Must keep a window of length at least 2.\"\n",
        "\n",
        "    trans_mat, obs_mat = params.trans_mat, params.obs_mat\n",
        "    n_states, n_obs = obs_mat.shape\n",
        "    \n",
        "    # If trans_mat is independent of action, adjust shape\n",
        "    if len(trans_mat.shape) < 3:\n",
        "        trans_mat = jnp.expand_dims(trans_mat, axis=0)\n",
        "        act = None\n",
        "    if act is None:\n",
        "        act = jnp.zeros(shape=(curr_len+1,), dtype=jnp.int8)\n",
        "\n",
        "    # Shift window forward by 1\n",
        "    if curr_len == win_len:\n",
        "        alpha_win = alpha_win[1:]\n",
        "        bmatrix_win = bmatrix_win[1:]\n",
        "        \n",
        "    # Perform one forward operation\n",
        "    new_alpha, _ = normalize(\n",
        "        obs_mat[:, obs] * (alpha_win[-1][:, None] * trans_mat[act[-1]]).sum(axis=0)\n",
        "    )\n",
        "    alpha_win = jnp.concatenate((alpha_win, new_alpha[None, :]))\n",
        "\n",
        "    # Smooth inside the window in parallel\n",
        "    def update_bmatrix(bmatrix):\n",
        "        return (bmatrix @ trans_mat[act[-2]]) * obs_mat[:, obs]\n",
        "    bmatrix_win = vmap(update_bmatrix)(bmatrix_win)\n",
        "    bmatrix_win = jnp.concatenate((bmatrix_win, jnp.eye(n_states)[None, :]))\n",
        "    \n",
        "    # Compute beta values by row-summing bmatrices\n",
        "    def get_beta(bmatrix):\n",
        "        return normalize(bmatrix.sum(axis=1))[0]\n",
        "    beta_win = vmap(get_beta)(bmatrix_win)\n",
        "    \n",
        "    # Compute posterior values\n",
        "    gamma_win, _ = normalize(alpha_win * beta_win, axis=1)\n",
        "    return alpha_win, bmatrix_win, gamma_win"
      ],
      "metadata": {
        "id": "6YFXqHP00UmQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 2. Helper Functions"
      ],
      "metadata": {
        "id": "U_eoeO7IgQll"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We define some helper functions to run the two implementations on any given dataset."
      ],
      "metadata": {
        "id": "Muz8w0oYgSvd"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def get_fls_result_iterative(params, data, win_len, act=None):\n",
        "    assert data.size > 2, \"Complete observation set must be of size at least 2\"\n",
        "    prior, obs_mat = params.init_dist, params.obs_mat\n",
        "    alpha, _ = normalize(prior * obs_mat[:, data[0]])\n",
        "    obs_seq = jnp.array([data[0]])\n",
        "    for obs in data[1:]:\n",
        "        alpha, obs_seq, beta, gamma = fixed_lag_smoother_iterative(params, win_len, alpha, obs_seq, obs)\n",
        "    return alpha, beta, gamma\n",
        "\n",
        "def get_fls_result_vectorized(params, data, win_len, act=None):\n",
        "    assert data.size > 2, \"Complete observation set must be of size at least 2\"\n",
        "    prior, obs_mat = params.init_dist, params.obs_mat\n",
        "    n_states = obs_mat.shape[0]\n",
        "    alpha, _ = normalize(prior * obs_mat[:, data[0]])\n",
        "    bmatrix = jnp.eye(n_states)[None, :]\n",
        "    for obs in data[1:]:\n",
        "        alpha, bmatrix, gamma = fixed_lag_smoother_vectorized(params, win_len, alpha, bmatrix, obs)\n",
        "    return alpha, gamma"
      ],
      "metadata": {
        "id": "RH_qnaLtgZX5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 3. Correctness of Implementation"
      ],
      "metadata": {
        "id": "tEg7cAk0aNAH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, we randomly generate a Hidden Markov Model with state space of size 3, observation space of size 5.\n"
      ],
      "metadata": {
        "id": "WXfW503oaXXW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "num_states = 3\n",
        "num_obs = 5\n",
        "key = PRNGKey(0)\n",
        "\n",
        "# Generate transition matrix\n",
        "trans_mat = uniform(key, shape=(num_states, num_states))\n",
        "trans_mat, _ = normalize(trans_mat, axis=1)\n",
        "key, _ = split(key)\n",
        "\n",
        "# Generate observation likelihood matrix\n",
        "obs_mat = uniform(key, shape=(num_states, num_obs))\n",
        "obs_mat, _ = normalize(obs_mat, axis=1)\n",
        "key, _ = split(key)\n",
        "\n",
        "# Generate uniform random prior\n",
        "prior, _ = normalize(jnp.ones(num_states))\n",
        "\n",
        "# Generate HMM\n",
        "hmm = HMMJax(trans_mat=trans_mat, obs_mat=obs_mat, init_dist=prior)"
      ],
      "metadata": {
        "id": "36-3uX2_a3sC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then, we randomly generate an observed dataset of size 20."
      ],
      "metadata": {
        "id": "Fsi7mh1Ab7E3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "num_timesteps = 20\n",
        "\n",
        "# Generate random data\n",
        "data = choice(key, num_obs, shape=(num_timesteps,))"
      ],
      "metadata": {
        "id": "hDqPeGlLcLOS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now, we test correctness of the two implementations by comparison of their full-lag results with JSL's `hmm_forwards_backwards_jax` implementation."
      ],
      "metadata": {
        "id": "Ax8eZ9SEe9sM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "_, _, gamma_fb, _ = hmm_forwards_backwards_jax(hmm, data)\n",
        "*_, gamma_it = get_fls_result_iterative(hmm, data, win_len=data.size)\n",
        "*_, gamma_vec = get_fls_result_vectorized(hmm, data, win_len=data.size)"
      ],
      "metadata": {
        "id": "RjrgZNuQe0DB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "assert jnp.allclose(gamma_fb, gamma_it) and jnp.allclose(gamma_fb, gamma_vec)"
      ],
      "metadata": {
        "id": "EWNY_tBjfpeM"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 4. Timed Experiments"
      ],
      "metadata": {
        "id": "S2ZVMRwEkQML"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We compare the performance of the two implementations by timing them against different sizes of windows. We reuse the HMM from the previous section. We use total observation set of size 1000.\n"
      ],
      "metadata": {
        "id": "hhWfZkxckwFL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot timed results against different window sizes\n",
        "import numpy as np\n",
        "\n",
        "key = PRNGKey(0)\n",
        "num_timestep = 1000\n",
        "window_grid = np.insert(np.arange(10, 81, 10), 0, 2)\n",
        "data = choice(key, num_obs, shape=(num_timestep,))\n",
        "\n",
        "it_times, vec_times = [], []\n",
        "for window_size in window_grid:\n",
        "    print(\"window of size \", window_size)\n",
        "    it_time = %timeit -o get_fls_result_iterative(hmm, data, window_size)\n",
        "    vec_time = %timeit -o get_fls_result_vectorized(hmm, data, window_size)\n",
        "    it_times.append(it_time.best)\n",
        "    vec_times.append(vec_time.best)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "sx0V530-lGYg",
        "outputId": "a4fdae5d-ff21-4d83-923a-0cdb87f5b19b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "window of size  2\n",
            "The slowest run took 51.57 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 21.6 ms per loop\n",
            "The slowest run took 20.09 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 20.6 ms per loop\n",
            "window of size  10\n",
            "The slowest run took 243.44 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 26.6 ms per loop\n",
            "The slowest run took 94.51 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 21.3 ms per loop\n",
            "window of size  20\n",
            "The slowest run took 632.30 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 35.5 ms per loop\n",
            "The slowest run took 179.90 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 22.2 ms per loop\n",
            "window of size  30\n",
            "The slowest run took 1196.74 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 42.2 ms per loop\n",
            "The slowest run took 253.59 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 23.3 ms per loop\n",
            "window of size  40\n",
            "The slowest run took 1948.70 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 46.7 ms per loop\n",
            "The slowest run took 334.61 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 23.7 ms per loop\n",
            "window of size  50\n",
            "The slowest run took 2801.64 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 53.1 ms per loop\n",
            "The slowest run took 406.91 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 24.7 ms per loop\n",
            "window of size  60\n",
            "The slowest run took 3677.37 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 58.8 ms per loop\n",
            "The slowest run took 454.20 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 26.2 ms per loop\n",
            "window of size  70\n",
            "The slowest run took 4519.96 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 72.9 ms per loop\n",
            "The slowest run took 611.26 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 25.9 ms per loop\n",
            "window of size  80\n",
            "The slowest run took 5452.36 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 71.8 ms per loop\n",
            "The slowest run took 618.85 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "1 loop, best of 5: 26.4 ms per loop\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot the result\n",
        "plt.figure()\n",
        "plt.plot(window_grid, it_times, label='iterative')\n",
        "plt.plot(window_grid, vec_times, label='vectorized')\n",
        "plt.xlabel('window length')\n",
        "plt.ylabel('best time in s')\n",
        "plt.legend();\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "id": "QKMz5jOqvssT",
        "outputId": "d4338dec-28d8-467a-bd66-b77e1e29ed77"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUVfrA8e9LCiEkhBBCL0Gq1IChCSo2BBvWBXFdLCvLrq66bhH92XBdRde1t7UgViyoiGJhLQiKhVCk95bQSSCN9Ly/P+4NDGGAETKZkvfzPPPMLWdm3kyS8957zr3niKpijDHGVFUn0AEYY4wJTpYgjDHGeGUJwhhjjFeWIIwxxnhlCcIYY4xXkYEOoLo0btxYU1JSAh2GMcaElPnz5+9W1WRv+8ImQaSkpJCenh7oMIwxJqSIyKbD7bMmJmOMMV5ZgjDGGOOVJQhjjDFehU0fhDelpaVkZmZSVFQU6FDCTkxMDK1atSIqKirQoRhj/CSsE0RmZibx8fGkpKQgIoEOJ2yoKllZWWRmZtKuXbtAh2OM8ZOwbmIqKioiKSnJkkM1ExGSkpLszMyYMBfWCQKw5OAn9r0aE/7CuonJGGOqw+zVu1i1PY9mCTE0T4ihWUIMTeJjiI4M72NsSxB+dvLJJzN37lw2btzI3LlzGT16dLW99wMPPMAdd9xxyGcZY6rPj+uzuPqVn6moMnWOCDSOq0vzhBiaNjiQOJonxNCsQb396zFREYEJvBpIuEwYlJaWplXvpF6xYgUnnnhigCI62KxZs3jkkUf45JNPfH5NWVkZkZGHz+FxcXHk5+dXR3jHJJi+X2P8IbughHOfmENsdARvXT+AnMJStuUUsj2niG05RezIdZ6d9UJyi8oOeY/E2CiaJdSjWYO6NEuod1AicZbrEVc3cMfqIjJfVdO87bMzCD+rrMTHjx/PihUrSE1NZcyYMdx0002MHz+eWbNmUVxczA033MAf/vAHZs2axV133UViYiIrV65k9erVXHTRRWRkZFBUVMTNN9/M2LFjGT9+PIWFhaSmptKtWzfefPPN/Z81atQorrrqKs477zwArr76as4//3wuvvhir59pjDmUqvKPqb+QXVDCS2NOpplbsXduFn/Y1xQUl7E9t2h/AtmeU3jQ+uLMHLIKSg55XXzdyP3v36zBgcThmUwS6kXVeN9frUkQEz5exvKtudX6nl1bNOCeC7r5VHbixIkHnUG88MILJCQkMG/ePIqLixk0aBBDhw4FYMGCBSxdunT/JaSTJk2iUaNGFBYW0rdvXy699FImTpzI008/zaJFiw75rJEjR/Luu+9y3nnnUVJSwldffcVzzz3Hyy+/7PUz7VJVYw41ee5Gvlyxk3su6Er3lgk+vaZ+3UjaJ8fRPjnusGWKSsvZmVvsnIl4nIFszyliW24Rq3fsYmdeMVUbd2Ki6tA8oR7NGsTsTybN3YRyQnIcHZoc/jOPVa1JEMFm5syZLF68mKlTpwKQk5PDmjVriI6Opl+/fgdV2k8++SQffvghABkZGaxZs4akpKTDvvfw4cO5+eabKS4u5vPPP+fUU0+lXr16h/1MSxDGHGzplhwe/HQlZ53YhKtPTqnW946JiqBNUixtkmIPW6a0vIJdecUHNV95Nmf9vCGbHblFlLkdI+f1bM4zo/tUa5xQixKEr0f6NUVVeeqppzjnnHMO2j5r1izq169/0PqXX37JDz/8QGxsLEOGDDnq/QcxMTEMGTKEL774gnfeeYdRo0Yd8TONMQfkF5fx5ykLaVQ/mn9f1isgl3RHRdShRcN6tGhY77BlKiqU3QXFbM8pom6kfzrCw/sarSASHx9PXl7e/vVzzjmH5557jtLSUgBWr15NQUHBIa/LyckhMTGR2NhYVq5cyY8//rh/X1RU1P7XVzVy5EheeeUV5syZw7Bhw37VZxpTm9390VI2ZRXw+KhUEutHBzqcw6pTR2gSH0PPVg2P2C9yPGrNGUSg9ezZk4iICHr16sXVV1/NzTffzMaNG+nTpw+qSnJyMtOmTTvkdcOGDeP555/nxBNPpHPnzgwYMGD/vrFjx9KzZ0/69OnDm2++edDrhg4dylVXXcWIESOIjnb+yH//+9/79JnG1FYfLszkgwVbuPnMjgw44fDNuLWFXy9zFZFhwBNABPCSqk6ssr8u8BpwEpAFjFTVjSJyJfB3j6I9gT6qemiPrCvYL3MNR/b9mnCyYXcB5z85h24tE3jr9/2JjKgdDSxHuszVb9+AiEQAzwDDga7AFSLStUqx64A9qtoBeAx4CEBV31TVVFVNBa4CNhwpORhjzPEoLivnz1MWEBVZhydGpdaa5HA0/vwW+gFrVXW9qpYAbwMjqpQZAbzqLk8FzpRDe4SucF9rjDF+8fDnq1i6JZeHL+1J84TDdwzXNv5MEC2BDI/1THeb1zKqWgbkAFUb/kYCU7x9gIiMFZF0EUnftWtXtQRtjKldvl65g5e/28DVJ6cwtFuzQIcTVIL6PEpE+gP7VHWpt/2q+oKqpqlqWnJycg1HZ4wJddtzivjbe4s5sXkDxg/vEuhwgo4/E8QWoLXHeit3m9cyIhIJJOB0VlcaxWHOHowx5niUVyi3vLOQotJynh7dO6QH1fMXfyaIeUBHEWknItE4lf30KmWmA2Pc5cuAr9W9rEpE6gC/wfofjDF+8Ow3a/lxfTYTLux2xKExajO/JQi3T+FG4AtgBfCuqi4TkftE5EK32MtAkoisBW4Fxnu8xalAhqqu91eMoWLatGksX778V7/u+eef57XXXjvuz09JSWH37t3H/T7GBIt5G7N57MvVXJTagstOahXocIKWX2+UU9VPgU+rbLvbY7kIuPwwr50FDPC2r7aZNm0a559/Pl27Vr1K+PDKysoYN26cH6MyJjTt3VfCzVMW0rpRLPdf3MNmRzwCu5Paz8aPH0/r1q254YYbALj33nuJi4tDVXn33XcpLi7m4osvZsKECQC89tprPPLII4gIPXv25I9//CPTp0/n22+/5f777+f9998nLy+PcePGsW/fPtq3b8+kSZNITExkyJAhpKam8t1333HFFVeQl5dHXFwco0eP5txzz90f05IlS1i/fj2xsbGMGzeOzZs3A/D4448zaNAgsrKyuOKKK9iyZQsDBw4kXOYMMUZVue39xezKL+b9P54c0HkYQkHt+XY+Gw/bl1TvezbrAcMnHrHIyJEjueWWW/YniHfffZfbbruN77//np9//hlV5cILL2T27NkkJSVx//33M3fuXBo3bkx2djaNGjXiwgsv5Pzzz+eyyy4DnGE7nnrqKU477TTuvvtuJkyYwOOPPw5ASUkJlXeU33vvvQC0aNFi/7DgzzzzDN9++y1t27Zl9OjR/OUvf2Hw4MFs3ryZc845hxUrVjBhwgQGDx7M3XffzYwZM3j55Zer93szJkDe+HETXyzbwZ3nnUjPVg0DHU7Qqz0JIkB69+7Nzp072bp1K7t27SIxMZElS5Ywc+ZMevfuDUB+fj5r1qzhl19+4fLLL6dx48YANGrU6JD3y8nJYe/evZx22mkAjBkzhssvP9BKN3LkyMPG8v333/Piiy/y3XffAfDll18e1LeRm5tLfn4+s2fP5oMPPgDgvPPOIzEx8Ti/BWMCb8W2XP45YwWnd07m2kE2xL0vak+COMqRvj9dfvnlTJ06le3btzNy5Eg2bdrE7bfffshsbk899dRxf5bnUOGetm3bxnXXXcf06dOJi3Ou2KioqODHH38kJibmuD/XmGC2r6SMG99aQMN6UTxyeS/q1LF+B18E9Y1y4WLkyJG8/fbbTJ06lcsvv5xzzjmHSZMm7Z9PesuWLezcuZMzzjiD9957j6ws51aQ7Oxs4OChwhMSEkhMTGTOnDkAvP766/vPJg6ntLSUyy+/nIceeohOnTrt3z506NCDklJlM9Spp57KW2+9BcBnn33Gnj17quNrMCZgJkxfzvrdBTw+MpWkuLqBDidkWIKoAd26dSMvL4+WLVvSvHlzhg4dyujRoxk4cCA9evTgsssuIy8vj27duvF///d/nHbaafTq1Ytbb70VgFGjRvHvf/+b3r17s27dOl599VX+/ve/07NnTxYtWsTdd999xM+fO3cu6enp3HPPPaSmppKamsrWrVt58sknSU9Pp2fPnnTt2pXnn38egHvuuYfZs2fTrVs3PvjgA9q0aeP378gYf5n+y1beSc/ghiEdOLlD40CHE1L8Otx3TbLhvmuefb8m2G3O2se5T86hc7N43hk7wEZp9SIgw30bY0wglZRV8OcpC6gj2BDex6j2dFIbY2qV/8xcxS+ZOTz/2z60SowNdDghKexTarg0oQUb+15NMPt29S7+O3s9vx3QhmHdmwc6nJAV1gkiJiaGrKwsq8yqmaqSlZVll8eaoLQzr4i/vruILs3iufM834enMYcK6yamVq1akZmZiU0mVP1iYmJo1coGOTPBpaJCufWdX8gvLmPK9QNsCO/jFNYJIioqinbt7I5JY2qL52ev47u1u5l4SQ86No0PdDghL6ybmIwxtcf8TXv4z8zVnNezOSP7tj76C8xRWYIwxoS8nMJSbpqykOYJMTx4iQ3hXV3CuonJGBP+VJU7PljCjtwi3hs3kAYxUYEOKWzYGYQxJqRN+TmDGUu28bdzOtO7jY08XJ0sQRhjQtbqHXlM+HgZp3RszNhTTgh0OGHHEoQxJiQVlpRz41sLiI+J4tHfpNoQ3n5gfRDGmJD0zxnLWb0jn9eu7UdyvA3h7Q92BmGMCTmfLtnGWz9tZtxp7Tm1U3KgwwlbliCMMSElI3sft72/mNTWDfnr0E5Hf4E5ZpYgjDEho7S8gpvfXggKT13RmygbwtuvrA/CGBMyHvvfahZs3svTo3vTupEN4e1vln6NMSHhuzW7ee7bdYzq25rze7YIdDi1giUIY0zQ251fzF/eXUT75DjuuaBboMOpNayJyRgT1CoqlL+++wu5haW8fl0/6kXbEN41xc4gjDFB7aXv1vPt6l3cdX5XujRrEOhwahVLEMaYoPVLxl4e/nwVw7s348r+bQIdTq3j1wQhIsNEZJWIrBWR8V721xWRd9z9P4lIise+niLyg4gsE5ElImLzWxpTi+QVlfLnKQtp2iCGiZf0tCG8A8BvCUJEIoBngOFAV+AKEak6Qex1wB5V7QA8BjzkvjYSeAMYp6rdgCFAqb9iNcYEF1Xl/z5cypa9hTx5RSoJsTaEdyD48wyiH7BWVderagnwNjCiSpkRwKvu8lTgTHEOE4YCi1X1FwBVzVLVcj/GaowJIu/Nz2T6L1u59exOnNS2UaDDqbX8mSBaAhke65nuNq9lVLUMyAGSgE6AisgXIrJARP7h7QNEZKyIpItI+q5du6r9BzDG1Ly1O/O456NlnNw+iXGntQ90OLVasHZSRwKDgSvd54tF5MyqhVT1BVVNU9W05GQbsMuYUFdUWs6Nby0kNjqCx0amEmFDeAeUPxPEFsBz5vBW7javZdx+hwQgC+dsY7aq7lbVfcCnQB8/xmqMCQIPfLqCldvzeOQ3vWjawK5LCTR/Joh5QEcRaSci0cAoYHqVMtOBMe7yZcDXqqrAF0APEYl1E8dpwHI/xmqMCbAvlm3ntR828fvB7Ti9c5NAh2Pw453UqlomIjfiVPYRwCRVXSYi9wHpqjodeBl4XUTWAtk4SQRV3SMij+IkGQU+VdUZ/orVGBM4WfnFPPnVGt78aTM9Wibwj2FdAh2ScYlzwB760tLSND09PdBhGGN8VFhSzqTvN/DcrHUUlpYzqm9r/jq0M43qRwc6tFpFROarapq3fTYWkzGmRpVXKB8syOQ/M1ezPbeIs7s25bZhXejQJC7QoZkqLEEYY2rMt6t38aDbEd2rdUOeGJVK/xOSAh2WOQxLEMYYv1u2NYeJn61kzprdtGkUy9Oje3Nej+Y2fEaQswRhjPGbrXsLeWTmKj5cuIWEelHcdX5XfjugDXUjbcjuUGAJwhhT7XKLSnlu1jomfbcBBcaeegJ/GtKBhHo2plIosQRhjKk2JWUVvPnTJp78ag179pVySe+W3Dq0E60Sbf7oUGQJwhhz3FSVT5ds5+EvVrIpax+DOiRx+/AT6d4yIdChmeNgCcIYc1zmbczmXzNWsChjL52bxvPKNX0Z0inZOqDDgCUIY8wxWbcrn4c+W8nM5Tto2qAuD1/ak0tPamUD7IURSxDGmF9ld34xj3+5mik/ZxATWYe/De3EtYPbERtt1Um4sd+oMcYnhSXlvDRnPc9/u46isgpG92vDzWd1pHFc3UCHZvzEEoQx5ojKK5Sp8zN49H+r2ZFbzDndmvKPYV1on2xDY4S7oyYIEakPFKpqhYh0AroAn6mqzRFtTBhTVWat3sXET1eyakcevds05OnRfeibYlOA1ha+nEHMBk4RkURgJs4Q3CNxZnszxoShpVtyeODTFcxdl0XbpFievbIPw7s3syuTahlfEoSo6j4RuQ54VlUfFpFF/g7MGFPzMvfs4z8zV/Phwi0kxkZxzwVdubJ/W6Ijg3V2YuNPPiUIERmIc8ZwnbvNBlIxJozkFJby7DdreWXuRgT445D2/HFIexrE2NAYtZkvCeJm4HbgQ3dGuBOAb/wbljGmJhSXlfPGj5t56us15BSWcnHvlvx1aGdaNqwX6NBMEDhqglDV2Tj9EJXr64Gb/BmUMca/VJVPFm/j4S9WkpFdyCkdGzN+eBe6tbChMcwBdpmrMbXMT+uzeODTFfySmUOXZvG8em0/TuuUHOiwTBCyBGFMLbFlbyETpi9j5vIdNGsQw78v68klfWxoDHN4liCMCXNl5RW88v1GHvtyNarwt6GduG7wCdSLtmtNzJH5cqNcMnA9kOJZXlWv9V9YxpjqsChjL3d8sITl23I5o0sT7hvRzeZmMD7z5QziI2AO8CVQ7t9wjDHVIbeolH9/voo3ftpEk/i6PHdlH4bZjW7mV/IlQcSq6m1+j8QYc9xUlRlLtjHh4+Vk5RczZmAKfx3aiXi7n8EcA18SxCcicq6qfur3aIwxxywjex93TlvKt6t30b1lA14ek0bPVg0DHZYJYb7eKHeHiBQDpYAAqqoN/BqZMcYnpeUVvDhnPU9+tYYIEe4+vyu/G9iWyAgbHsMcH19ulIuviUCMMb/e/E3Z3PHBUlbtyOOcbk2598JuNE+wu6BN9ThsghCRLqq6UkT6eNuvqgv8F5Yx5khy9pUy8fOVTPl5My0SYnjxd2mc3bVpoMMyYeZIZxC3AmOB/3jZp8AZfonIGHNYqspHi7Zy/4zl7NlXyu8Ht+MvZ3eifl27pclUv8P+VanqWPf59GN9cxEZBjyBM/rrS6o6scr+usBrwElAFjBSVTeKSAqwAljlFv1RVccdaxzGhIMNuwu4a9pSvlu7m16tG/Lqtd1t7CTjV3477BCRCOAZ4GwgE5gnItNVdblHseuAParaQURGAQ/hTEYEsE5VU/0VnzGhorisnP9+u56nv1lL3Yg63DeiG1f2b2tDZBi/8+d5aT9grTv6KyLyNjAC8EwQI4B73eWpwNNid/IYs9+P67P4vw+XsG5XAef1bM7d53elaYOYQIdlagl/JoiWQIbHeibQ/3BlVLVMRHKAJHdfOxFZCOQCd6rqnKofICJjcfpJaNOmTfVGb0wAZReU8OCnK3hvfiatEuvxyjV9Ob1zk0CHZWoZX8ZiEpzZ5E5Q1ftEpA3QTFV/9mNc24A2qpolIicB00Skm6rmehZS1ReAFwDS0tLUj/EYUyNUlfcXbOFfM5aTV1TGuNPac/OZHW1gPRMQvpxBPAtU4Fy1dB+QB7wP9D3K67YArT3WW7nbvJXJFJFIIAHIUlUFigFUdb6IrAM6Aek+xGtMSFq7M587py3hx/XZnNQ2kX9d3J0uzex+VBM4viSI/qrax23uQVX3iEi0D6+bB3QUkXY4iWAUMLpKmenAGOAH4DLga1VVdwTZbFUtd6c47Qis9+1HMia0FJWW8+w3a3nu23XUi4rggYt7MKpva+pYJ7QJMF8SRKl7RZLC/uG/K472IrdP4UbgC5zLXCe5c1rfB6Sr6nTgZeB1EVkLZOMkEYBTgftEpNT9rHGqmv0rfzZjgt73a3dz57SlbNhdwIjUFtx5XleS4+sGOixjABCnNecIBUSuxLn0tA/wKs6R/p2q+p7/w/NdWlqapqdbC5QJDbvzi/nXjBV8uHALKUmx/POi7pzS0ab9NDVPROarapq3fb6MxfSmiMwHzsQZqO8iVV1RzTEaUytUVCjvpmfw4Gcr2VdSxp/P6MANp3cgJso6oU3w8fUy1x04kwZFAvVEpI+NxWTMr7N6Rx53fLCE9E176NeuEQ9c3J0OTWwsTBO8fLnM9Z/A1cA63H4IbCwmY3xWWFLOk1+v4cXZ64mLieThy3py+UmtbHY3E/R8OYP4DdBeVUv8HYwx4WbWqp3c9dFSMrILubRPK+44twtJcdYJbUKDLwliKdAQ2OnnWIwJGztzi7jvk+V8sngbJyTXZ8r1AxjYPunoLzQmiPiSIB4EForIUtyb1wBU9UK/RWVMiNpXUsbbP2fw2P9WU1xewV/O6sS4ISdQN9I6oU3o8SVBvIozyuoSfLj/wZjaKHPPPl7/YRNTft5MblEZJ7dP4v6LunNCclygQzPmmPmSIPap6pN+j8SYEKOq/Lwhm8lzN/LFsu2ICOd0a8o1g9qR1jbROqFNyPMlQcwRkQdxhsXwbGKyy1xNrVRUWs7Hv2xl8tyNLNuaS0K9KMae2p6rBralZUObD9qED18SRG/3eYDHNrvM1dQ6O3OLeOPHTbz502ayCkro1DSOBy7uwcW9W9poqyYs+XIn9TFPOWpMOFiUsZdXvt/AjMXbKFflzC5NuPrkdgzqkGTNSCasHTZBiMhvVfUNEbnV235VfdR/YRkTWKXlFXy2dDuvfL+BhZv3Elc3kqsGtmXMwBRSGtcPdHjG1IgjnUFU/hd4GwvAJucxYSm7oIQpP2/m9R82sT23iJSkWO65oCuXndSK+JioQIdnTI06bIJQ1f+6i1+q6vee+0RkkF+jMqaGrdiWy+TvNzJt0RaKyyoY3KEx/7q4O6d3bmLzMphay5dO6qdwhvo+2jZjQkp5hfLlih288v0GflyfTUxUHS49qRVXn5xCp6Y2iJ4xR+qDGAicDCRX6YdogDMBkDEhKaewlPfSM3j1h41kZBfSIiGG8cO7MKpvaxrG+jJZojG1w5HOIKKBOLeM5+FULs6kQcaElHW78nl17kamzs9kX0k5fVMSuX34iQzt2pTIiDqBDs+YoHOkPohvgW9FZLKqbqrBmIypNhUVyuw1u5g8dyOzVu0iOqIOF/RqwTWDUujeMiHQ4RkT1Hy5D8KSgwk5BcVlfLAgk8lzN7JuVwHJ8XX5y1mdGN2/jc35bIyPfJ1RzpiQkJG9j9d+2Mjb8zLIKyqjZ6sEHhvZi/N6tCA60pqRjPk1fJlRbpC3y1yrbjMmUFSVnzZk88r3G/jf8h2ICMO7N+OaQSn0aWOD5hlzrOwyVxOyikrLmb5oK6/M3ciKbbkkxkYx7jRn0LzmCTZonjHHyy5zNSGnqLSc/367nld/2Eh2QQmdm8Yz8ZIeXNS7JTFR9qdpTHWxy1xNSPl5Qzbj31/M+t0FnHViE64d1I6B7W3QPGP84Vdd5ioidYA4Vc2tqQCNAcgrKmXiZyt586fNtG5Ujzeu68/gjo0DHZYxYc2nOalFZBxQDswDGojIE6r6b/+GZozjqxU7uHPaUnbkFnHd4Hb8dWgnYqPtAjxj/M2X/7KuqporIlcCnwHjgfmAJQjjV1n5xUz4eDnTf9lK56bxPPfbk0ht3TDQYRlTa/iSIKJEJAq4CHhaVUtFxIb7Nn6jqny0aCsTPl5GfnEZfzmrE38c0t7uYzCmhvnyH/dfYCPO/BCzRaQtTkf1UYnIMBFZJSJrRWS8l/11ReQdd/9PIpJSZX8bEckXkb/58nkm9G3ZW8g1k+dxyzuLSGlcnxk3ncLNZ3W05GBMAPgy1MaTwJMemzaJyFGnIRWRCOAZ4GwgE5gnItNVdblHseuAParaQURGAQ8BIz32P4rTrGXCXEWF8sZPm3jos5UocM8FXfndwBQibC4GYwLGlzupmwIPAC1UdbiIdAUGAi8f5aX9gLWqut59n7eBEYBnghgB3OsuTwWeFhFRVRWRi4ANQMGv+HlMCFq7M5/x7y8mfdMeTunYmAcu7kHrRrGBDsuYWs+X8/bJwBdAC3d9NXCLD69rCWR4rGe627yWUdUyIAdIEpE44DZggg+fY0JUaXkFT3+9hnOfmMOanfn85/JevHZtP0sOxgQJXzqpG6vquyJyOzgVuYiU+zmue4HHVDX/SDdAichYYCxAmzZt/BySqU6LM/fyj6mLWbk9j/N6NufeC7rZKKvGBBlfEkSBiCQBCiAiA3CO9I9mC9DaY72Vu81bmUwRiQQSgCygP3CZiDwMNAQqRKRIVZ/2fLGqvgC8AJCWlmZXVoWAwpJyHvtyNS/NWU9yfF1e/F0aZ3dtGuiwjDFe+JIgbgWmA+1F5HsgGd+G2pgHdBSRdjiJYBQwukqZ6cAY4Af3Pb9WVQVOqSwgIvcC+VWTgwk9c9ft5vYPlrApax9X9GvD7ed2oUFMVKDDMsYchi9XMS0QkdOAzoAAq1S11IfXlYnIjTj9FxHAJFVdJiL3AemqOh2no/t1EVkLZOMkERNmcgpLefDTFbw9L4OUpFimXD+Age2TAh2WMeYoxDlgP0IBkRjgT8BgnGamOcDzqlrk//B8l5aWpunp6YEOw1TxxbLt3DVtKVkFJfz+lHb85axONuKqMUFEROarapq3fb40Mb0G5OHMAQFOM9HrwOXVE54JRzvzirh3+jI+XbKdE5s34OUxfenRyuaANiaU+JIguqtqV4/1b0Rk+WFLm1pNVZk6P5P7Z6ygsLScv5/TmbGnnkBUhN0JbUyo8SVBLBCRAar6I4CI9AesLcccIiN7H3d8uIQ5a3bTNyWRiZf2pH1yXKDDMsYcoyPNKLcEp88hCpgrIpvd9bbAypoJz4SC8gpl8tyNPPLFKuoI/POi7lzZrw11bJgMY0Lake3RGtkAABW4SURBVM4gzq+xKEzIWr0jj39MXcyijL2c0aUJ91/UnRYNbT5oY8LBkWaU21STgZjQUlxWzrPfrOPZWWuJj4niiVGpXNirhU39aUwYsWm5zK+2YPMebpu6mDU787kotQV3X9CNRvWjAx2WMaaaWYIwPisoLuORmauYPHcjzRvE8MrVfTm9S5NAh2WM8RNLEMYnc9bs4vYPlpC5p5DfDWzLP4Z1Ia6u/fkYE87sP9wc0d59JfzzkxW8vyCTE5Lr8964gfRNaRTosIwxNcAShDmsz5du485pS9m7r5QbT+/AjWd0sGEyjKlFLEGYQ6gqj/1vNU9+vZYeLRN47dr+dG3RINBhGWNqmCUIc5DisnJum7qYaYu28pu0Vvzr4h42TIYxtZQlCLPf3n0ljH19Pj9vyObv53TmT0Pa230NxtRiliAMABt3F3Dt5Hlk7inkySt6c2GvFkd/kTEmrFmCMMzflM31r81HVXnz+v52lZIxBrAEUet9sngrt777Cy0b1mPS1X1p17h+oEMyxgQJSxC1lKry3LfrePjzVfRNSeSFq9JItOEyjDEeLEHUQqXlFdw1bSlvz8vgwl4tePiynnZ/gzHmEJYgapncolJueHMBc9bs5sbTO3Dr2Z1s3gZjjFeWIGqRLXsLufaVeazblc/Dl/bkN31bBzokY0wQswRRSyzJzOHaV+dRVFLO5Gv6Mbhj40CHZIwJcpYgaoH/Ld/BTVMW0qh+NG/+qT+dmsYHOiRjTAiwBBHmXvl+A/d9spyeLRN4cUwaTeJjAh2SMSZEWIIIU+UVyj8/Wc7kuRsZ2rUpT4zqTb1ou1LJGOM7SxBhaF9JGTdNWcSXK3Zw3eB23HHuiUTYlUrGmF/JEkSY2ZlbxHWvprNsaw73jejG7wamBDokY0yIsgQRRlZtz+PayfPYs6+EF3+XxpknNg10SMaYEGYJIkzMXr2LG95cQL3oCN79w0C6t0wIdEjGmBDn15lgRGSYiKwSkbUiMt7L/roi8o67/ycRSXG39xORRe7jFxG52J9xhrq3f97MNZPn0TKxHtNuGGTJwRhTLfx2BiEiEcAzwNlAJjBPRKar6nKPYtcBe1S1g4iMAh4CRgJLgTRVLROR5sAvIvKxqpb5K95QVFGh/HvmKp6btY7TOiXz9OjexMdEBTosY0yY8OcZRD9graquV9US4G1gRJUyI4BX3eWpwJkiIqq6zyMZxADqxzhDUlFpOTe9vZDnZq1jdP82vDwmzZKDMaZa+bMPoiWQ4bGeCfQ/XBn3bCEHSAJ2i0h/YBLQFrjK29mDiIwFxgK0adOm2n+AYJWVX8zY1+czf9Mebh/ehbGnnmBTgxpjql3Qzkavqj+pajegL3C7iBxyC7CqvqCqaaqalpycXPNBBsD6Xflc8txclm7J4dkr+/CH02zeaGOMf/gzQWwBPIcLbeVu81pGRCKBBCDLs4CqrgDyge5+izRE/Lwhm0uem0t+URlTxg7g3B7NAx2SMSaM+TNBzAM6ikg7EYkGRgHTq5SZDoxxly8DvlZVdV8TCSAibYEuwEY/xhr0pi3cwm9f+olG9aP58E+D6NMmMdAhGWPCnN/6INw+hRuBL4AIYJKqLhOR+4B0VZ0OvAy8LiJrgWycJAIwGBgvIqVABfAnVd3tr1iDmary1NdrefR/q+nfrhH/veokGsba1KDGGP8T1fC4QCgtLU3T09MDHUa1Kimr4I4PlzB1fiaX9G7JxEt7Eh0ZtN1GxpgQJCLzVTXN2z67kzpI5RSW8sc35jN3XRa3nNWRm8/saJ3RxpgaZQkiCGVk7+OayfPYlFXAo7/pxSV9WgU6JGNMLWQJIsgsytjL71+dR0lZBa9d25+B7ZMCHZIxppayBBFEPl+6nVveWUhyfF3eHjuQDk3iAh2SMaYWswQRBFSVl+Zs4IHPVpDauiEv/i6NxnF1Ax2WMaaWswQRYGXlFdz78TLe+HEz5/ZoxqO/SSUmyqYGNcYEniWIANq6t5DxHyxh9upd/OG0E7jtnC7UsalBjTFBwhJEAJSUVTDp+w08+dUaKlR54OIejO5fewYbNOYQqpC/A8qKILIeRNaFqHoQEQ3hfnl3RQWUFkBxPpTkQ3Ge++zjekk+nDAEht5f7aFZgqhhP6zL4q6PlrJ2Zz5nd23K3ed3pXWj2ECHZUzNKM6DrLWwe63znLXGfV7nVHSHEIiMgagY57ny4e/1yLqHT0yqTiIrzoeSPI+K25d1bxV8AT7PaBAd5zzqVj7HQ4OWENfsWH8jR2QJoobszCvigRkrmLZoK60S6/HyGJsz2oSp8lLYs+ngBFCZEPK3exQUaNgGGneENgMhqQNExTqVb+WjtAjKCqGsGErdZ8/1kgLYt9tdr1K2ovT4fo79ycM9oykvO1Dha7nv77G/Qo93nusnQ6N2Byr4qhX+4daj6kOdmh1JwRKEn5WVV/DGj5v4z8zVFJdVcNMZHfjT6R2sI9qEtsomoay1sLvyLMB97NkIFR7Tt8QmOZV/h7Mgqb2TEJI6QGI758jdXyrK3URTmVyqJp7KZS+J55D1IqgT5VFxe1T4h12Pg4jQnsTLEoQfLdi8hzs/XMrybbmc0rEx943oTrvG9QMdljG+K85zmn88E8DuNW6TUN6BcpEx0Kg9NO0GXUdAkpsEktpDbKPAxF4nAqLrOw9zTCxB+EF2QQkPfbaSd9IzaNYghmev7MPw7s1sLCUTnMrLYO8m72cDeds8Cgo0bO1U/K37O8+NOzjPDVrVePOH8T9LENWookJ5Jz2Dhz5fSX5RGWNPPYGbzuxIXF37mk0AqULRXsjfCXnbnSagLPcsYPca2LPh4Caheo2cSr/9Gc4ZQOXZQKMT/NskZIKO1VzVZOmWHO6ctpRFGXvp164R91/UnU5N4wMdlglnpUVQsPNAxZ+/w1ne/7z9wHp5ycGvjajrVP5NusCJF7hnA24iCFSTkAk6liCOU05hKY/OXMXrP26iUf1oHv1NLy7u3dKak8yxqaiAwmy3kt8BeTuqVPwej6IcL28gTqdwfDOIawKNOznPcU3dRxNo2BYSWjlt9MYcgSWIY6SqTFu0hX/NWEl2QTFXDWjLrUM7k1AvtK9aMH5SUnCgos/bXqXC9zza3+n9Esqo2AOVfJMTnRuj9lf8zQ4s128c8lfOmOBhCeIYrN6Rx13TlvLThmx6tW7I5Gv60r1lQqDDMtVF1bmWv6zQucTR87F/2z6niad0X5Xt7mNf1sEVv7ebwKQO1G9yoHJv2gPiPY709x/1N3UumzSmhlmC+BUKist44qs1TPpuA3ExkTx4SQ9GprW28ZNqWkU5FOyColyngi4r8lJheyyXelTqlde9H7Ktsqy77OuNUAcR96aqGKeZJ64pNE91nr1V/LFJ1sxjgpolCB+oKp8t3c4/P1nOtpwiRqa15rbhXWhUPzrQoYWf4jzI3QZ5W6s8b4PcrQc6Y39NBR5R17n6JirWvTM21qnIo+o5FXZUPWf8n6gqj0O2eb7e430qtx1peAZjQpAliKPYsLuAuz9aypw1u+navAFPj+7DSW0TAx1W6Kkodyr2qpX+/op/m7PN8+arSjEJEN8CGjSHJl2d5/hmENPQS2UeWyUZ1LOjdGOOkSWIwygqLeeZb9by32/XUzeyDvde0JXfDmhLZITdDHSIolyPyt7bUf8296i/4uDX1Yl0OlgbNIfkLs519/HNoUELj+dmdiesMQFiCcKLr1bs4N6Pl5GRXchFqS2447wTaRJfC28QKi9zrrP32uTjkQS8dcAe7qi/clt8C2fQMrv71pigZQnCQ0b2PiZ8vJwvV+ygY5M4plw/gIHtkwIdVvUqLXI6eAt2QsFu5wqbgp2Qv+vA9srlfVkcMgyx51F/kxOhw5l21G9MmLIEARSXlfPSnA089fUaBGH88C5cO6gd0ZEhcHSr6nTsFriVev7OIy8X53p/n2h3GOL6yc4dtm0GOMvxTe2o35haqtYniOVbc7nxrQWs313A8O7NuOv8rrRoWC+wQVVUQOEejyP6w1X67v6yIu/vUy/xwHX2zXu6y24SqNxev7GzHG2TFhljDlbrE0RyfF3iYiKZfE1fhnRu4r8PKitxJjbJd5t2CrxU9PubdnYfPHhaJYlwKvfKSj6pY5UK32PZ7qg1xhwnSxDxdfnohkG/fuwkVae5Zn87vke7/v4j/d0HtnsdNwfnUszKyj2hFbTs7b3Cj2viXNZpzTvGmBpS6xMEcCA5lJc6HbMHVfCVlb9nInAr/vJi729Yr9GB9vxm3Q8+qq+fXKVpp77dXGWMCUp+TRAiMgx4AogAXlLViVX21wVeA04CsoCRqrpRRM4GJgLRQAnwd1X92i9BblsM7//eqfALs72XiYg+UOHXT3Yu26zanFNZ8ccmWdOOMSYs+C1BiEgE8AxwNpAJzBOR6aq63KPYdcAeVe0gIqOAh4CRwG7gAlXdKiLdgS+Aln4JNCYBkjtDyuCD2/frN3GfGztl7CjfGFPL+PMMoh+wVlXXA4jI28AIwDNBjADudZenAk+LiKjqQo8yy4B6IlJXVQ/TpnMcEtvCyNer/W2NMSbU+bPHsyWQ4bGeyaFnAfvLqGoZkANUvTPtUmCBt+QgImNFJF1E0nft2lVtgRtjjPFvgjhuItINp9npD972q+oLqpqmqmnJyck1G5wxxoQ5fyaILUBrj/VW7javZUQkEkjA6axGRFoBHwK/U9V1fozTGGOMF/5MEPOAjiLSTkSigVHA9CplpgNj3OXLgK9VVUWkITADGK+q3/sxRmOMMYfhtwTh9inciHMF0grgXVVdJiL3iciFbrGXgSQRWQvcCox3t98IdADuFpFF7sOPtzkbY4ypSlT16KVCQFpamqanpwc6DGOMCSkiMl9V07ztC+pOamOMMYFjCcIYY4xXYdPEJCK7gE1HKNIY5w7tYGSxHRuL7dhYbMcmXGNrq6pe7xMImwRxNCKSfrh2tkCz2I6NxXZsLLZjUxtjsyYmY4wxXlmCMMYY41VtShAvBDqAI7DYjo3FdmwstmNT62KrNX0Qxhhjfp3adAZhjDHmV7AEYYwxxquwTxAiMkxEVonIWhEZf/RX+D2eSSKyU0SWemxrJCL/E5E17nNiAOJqLSLfiMhyEVkmIjcHUWwxIvKziPzixjbB3d5ORH5yf7fvuINCBoSIRIjIQhH5JJhiE5GNIrLEHc8s3d0W8N+pG0dDEZkqIitFZIWIDAyG2ESks8cYcItEJFdEbgmG2Nz4/uL+HywVkSnu/4df/t7COkF4THs6HOgKXCEiXQMbFZOBYVW2jQe+UtWOwFccGLSwJpUBf1XVrsAA4Ab3uwqG2IqBM1S1F5AKDBORAThzhTymqh2APThT2AbKzTiDUlYKpthOV9VUj+vkg+F3Cs589Z+rahegF873F/DYVHWV+32lAicB+3CmHgh4bCLSErgJSFPV7kAEzkjZ/vl7U9WwfQADgS881m8Hbg+CuFKApR7rq4Dm7nJzYFUQxPgRznziQRUbEAssAPrj3Dka6e13XcMxtcKpMM4APgEkiGLbCDSusi3gv1OcuV824F4oE0yxVYlnKPB9sMTGgVk4G+FMGf0JcI6//t7C+gwC36Y9DQZNVXWbu7wdaBrIYEQkBegN/ESQxOY24SwCdgL/A9YBe9UZVh4C+7t9HPgHUOGuJxE8sSkwU0Tmi8hYd1sw/E7bAbuAV9ymuZdEpH6QxOZpFDDFXQ54bKq6BXgE2Axsw5mmeT5++nsL9wQRctQ5BAjYtcciEge8D9yiqrme+wIZm6qWq3PK3wroB3QJRBxVicj5wE5VnR/oWA5jsKr2wWlmvUFETvXcGcDfaSTQB3hOVXsDBVRpsgmC/4Vo4ELgvar7AhWb2+8xAifBtgDqc2iTdbUJ9wThy7SnwWCHiDQHcJ93BiIIEYnCSQ5vquoHwRRbJVXdC3yDcxrdUJypaiFwv9tBwIUishF4G6eZ6Ykgia3yiBNV3YnTjt6P4PidZgKZqvqTuz4VJ2EEQ2yVhgMLVHWHux4MsZ0FbFDVXapaCnyA8zfol7+3cE8Qvkx7Ggw8p14dg9P+X6NERHBm+Fuhqo8GWWzJ4kxDi4jUw+kbWYGTKC4LZGyqeruqtlLVFJy/r69V9cpgiE1E6otIfOUyTnv6UoLgd6qq24EMEensbjoTWB4MsXm4ggPNSxAcsW0GBohIrPs/W/m9+efvLZAdQDXUqXMusBqnzfr/giCeKThth6U4R1HX4bRZfwWsAb4EGgUgrsE4p8yLgUXu49wgia0nsNCNbSlwt7v9BOBnYC1OM0DdAP9uhwCfBEtsbgy/uI9llX//wfA7deNIBdLd3+s0IDGIYqsPZAEJHtuCJbYJwEr3f+F1oK6//t5sqA1jjDFehXsTkzHGmGNkCcIYY4xXliCMMcZ4ZQnCGGOMV5YgjDHGeGUJwtRKIvJp5b0VPpZP8RyBt5piyK/O93PfM1VEzvVYv1dE/lbdn2NqB0sQplZS1XPVuSs73KTi3L9izHGzBGHCjoj8XURucpcfE5Gv3eUzRORNd3mjiDR2zwxWiMiL7hj7M927tRGRk8SZg+IX4AaP948RkVfceRYWisjp7vYZItLTXV4oIne7y/eJyPU+xDxPRBbLgfkujhRbX7fsIhH5tzs3QDRwHzDS3T7SffuuIjJLRNZXfi/G+MIShAlHc4BT3OU0IM4dZ+oUYLaX8h2BZ1S1G7AXuNTd/grwZ3XmofB0A854bT1whmN4VURiKj9XRBJw5tcY5JY/3OcCICJD3Rj64ZwBnOQxqN6RYvuDOgMYluMEVALcDbyjznwG77hlu+AMCd0PuMf9Low5KksQJhzNx6lkG+BMNvQDTqI4BacSr2qDqi7yeG2K2z/RUFUrK/bXPcoPBt4AUNWVwCagk/vep+Ikhhk4iSkWaKeqq44Q71D3sRBnrosuOInhSLHFq+oP7va3jvRlADNUtVhVd+MMMBfoIbRNiIg8ehFjQouqlorIBuBqYC7OWD+nAx04eNa3SsUey+VAvWP86Hk4iWg9zpwVjYHrcSr2IxHgQVX970EbnXk5qiO2qu9h//fGJ3YGYcLVHOBvOE07c4BxwEL1cfAxtwN7r4gMdjddWeW9rwQQkU5AG5zZxUpwJqi6HOesxTOGI/kCuNadiwMRaSkiTY4SW56I9Hc3jfLYnQfEH/0nNOboLEGYcDUHZ1rIH9QZz78I781LR3IN8Iw7k514bH8WqCMiS4B3gKtVtfIofQ7OBEKF7nKro32uqs7EaSb6wX3PqRy9kr8OeNGNrT7OzGLgDPvctUontTHHxEZzNSYEiUicqua7y+Nx5kq+OcBhmTBjbZHGhKbzROR2nP/hTTj9LcZUKzuDMMYY45X1QRhjjPHKEoQxxhivLEEYY4zxyhKEMcYYryxBGGOM8er/ASdlrLC9wwxBAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Note that the speedup is even more pronounced if we remove jit from the equation by running each code once."
      ],
      "metadata": {
        "id": "dcvmVu_OGCib"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot timed results against different window sizes - don't use cache\n",
        "import time\n",
        "import numpy as np\n",
        "\n",
        "key = PRNGKey(0)\n",
        "num_timestep = 1000\n",
        "window_grid = np.arange(3, 94, 10)\n",
        "data = choice(key, num_obs, shape=(num_timestep,))\n",
        "\n",
        "it_times, vec_times = [], []\n",
        "for window_size in window_grid:\n",
        "    print(\"window of size \", window_size)\n",
        "    start = time.time()\n",
        "    _ = get_fls_result_iterative(hmm, data, window_size)\n",
        "    it_time = time.time()-start\n",
        "    it_times.append(it_time)\n",
        "    print(\"iterative took \", it_time, \" seconds.\")\n",
        "    start = time.time()\n",
        "    _ = get_fls_result_vectorized(hmm, data, window_size)\n",
        "    vec_time = time.time()-start\n",
        "    vec_times.append(vec_time)\n",
        "    print(\"vectorized took \", vec_time, \" seconds.\")"
      ],
      "metadata": {
        "id": "XP07QTs1Fd8z",
        "outputId": "0f265cdf-73ca-4ae6-c9ff-88360e43f830",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "window of size  3\n",
            "iterative took  1.3882253170013428  seconds.\n",
            "vectorized took  0.579291820526123  seconds.\n",
            "window of size  13\n",
            "iterative took  9.687594413757324  seconds.\n",
            "vectorized took  3.332158327102661  seconds.\n",
            "window of size  23\n",
            "iterative took  29.268041372299194  seconds.\n",
            "vectorized took  4.641597032546997  seconds.\n",
            "window of size  33\n",
            "iterative took  63.42604160308838  seconds.\n",
            "vectorized took  6.314378976821899  seconds.\n",
            "window of size  43\n",
            "iterative took  106.45702171325684  seconds.\n",
            "vectorized took  8.493563413619995  seconds.\n",
            "window of size  53\n",
            "iterative took  162.71045875549316  seconds.\n",
            "vectorized took  9.996003150939941  seconds.\n",
            "window of size  63\n",
            "iterative took  232.62644505500793  seconds.\n",
            "vectorized took  11.84960126876831  seconds.\n",
            "window of size  73\n",
            "iterative took  323.3961088657379  seconds.\n",
            "vectorized took  14.298079490661621  seconds.\n",
            "window of size  83\n",
            "iterative took  437.9278955459595  seconds.\n",
            "vectorized took  16.133213758468628  seconds.\n",
            "window of size  93\n",
            "iterative took  574.5705089569092  seconds.\n",
            "vectorized took  17.906277418136597  seconds.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot the result\n",
        "plt.figure()\n",
        "plt.plot(window_grid, it_times, label='iterative')\n",
        "plt.plot(window_grid, vec_times, label='vectorized')\n",
        "plt.xlabel('window length')\n",
        "plt.ylabel('compute time in s')\n",
        "plt.legend();\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "08q8pi-6GUKn",
        "outputId": "1e1b2044-19fd-4c66-a691-d5a4a5ba18a2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        }
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXgUVdb48e/JRlhDCIuQAEEFZDEJEDZxYVBxQ3BUBFwGEEXndUbUGWfw9f3hOjM4OorboCggOCogbijujggIKFvY9z1hDxASQtY+vz+6EjohIU1Ip5P0+TxPP111q7rqpGju6ap765aoKsYYYwxAkL8DMMYYU3VYUjDGGFPIkoIxxphClhSMMcYUsqRgjDGmkCUFY4wxhXyaFESkoYjMFpGNIrJBRHqLSCMR+U5Etjjvkc66IiKviMhWEVktIl19GZsxxpjT+fpM4WXga1W9CIgHNgBjgR9UtS3wgzMPcB3Q1nmNBib6ODZjjDHFiK9uXhORCCAJOF89diIim4C+qrpPRJoD81S1vYi86Ux/UHy90vbRuHFjjY2N9Un8xhhTUy1fvvywqjYpaVmID/fbBjgETBWReGA5MAZo5lHR7weaOdPRwB6Pzyc7ZUWSgoiMxn0mQatWrVi2bJnP/gBjjKmJRGRXact8efkoBOgKTFTVLsAJTl0qAsA5gzirUxVVnaSqiaqa2KRJiYnOGGNMOfkyKSQDyar6izM/G3eSOOBcNsJ5P+gsTwFaenw+xikzxhhTSXyWFFR1P7BHRNo7RVcC64E5wHCnbDjwmTM9B/id0wupF5B2pvYEY4wxFc+XbQoAfwTeE5EwYDswEncimiUio4BdwG3Oul8C1wNbgUxn3bOWm5tLcnIyWVlZ5xq7KSY8PJyYmBhCQ0P9HYoxxkd8mhRUNQlILGHRlSWsq8AD57rP5ORk6tevT2xsLCJyrpszDlUlNTWV5ORk2rRp4+9wjDE+UuPuaM7KyiIqKsoSQgUTEaKiouwMzJgarsYlBcASgo/YcTWm5quRScEYY2oql0t59ov17E7N9Mn2LSn4wCWXXALAzp07ef/99yt023//+99L3JcxJjC8tWA7by/cwaJth32yfUsKPrBo0SKgfEkhLy/vjMuLJ4WCfRljar41yWm88O0mru10HkO6tyz7A+VgScEH6tWrB8DYsWNZsGABCQkJvPTSS+Tn5/Poo4/SvXt34uLiePPNNwGYN28el112GQMHDqRjx44A3HTTTXTr1o1OnToxadKkwu2dPHmShIQE7rjjjiL7Gjp0KHPnzi2MYcSIEcyePbvUfRpjqpcT2Xk8OGMlUXVrMf6Wi33Wxufr+xT86qnP17F+7/EK3WbHFg144sZOXq07fvx4XnjhBb744gsAJk2aREREBEuXLiU7O5s+ffrQv39/AFasWMHatWsLu3tOmTKFRo0acfLkSbp3784tt9zC+PHjee2110hKSjptX0OGDGHWrFnccMMN5OTk8MMPPzBx4kQmT55c4j6tW6kx1cvTn69nZ+oJ3r+nFw3rhPlsPzU6KVQ13377LatXr2b27NkApKWlsWXLFsLCwujRo0eRivqVV17hk08+AWDPnj1s2bKFqKioUrd93XXXMWbMGLKzs/n666+5/PLLqV27dqn7tKRgTPUxd/U+Zi7bwwO/uYDeF5ReD1SEGp0UvP1FX1lUlVdffZVrrrmmSPm8efOoW7dukfnvv/+exYsXU6dOHfr27Vvm/QHh4eH07duXb775hpkzZzJ06NAz7tMYUz2kHDvJYx+vJr5lQx66qp3P92dtCj5Uv3590tPTC+evueYaJk6cSG5uLgCbN2/mxIkTp30uLS2NyMhI6tSpw8aNG1myZEnhstDQ0MLPFzdkyBCmTp3KggULuPbaa89qn8aYqiffpTw8I4l8l/LK0ARCg31fZdfoMwV/i4uLIzg4mPj4eEaMGMGYMWPYuXMnXbt2RVVp0qQJn3766Wmfu/baa3njjTfo0KED7du3p1evXoXLRo8eTVxcHF27duW9994r8rn+/ftz1113MWjQIMLC3Ncc77nnHq/2aYypev7941Z+3XmEF2+Lp3VU3bI/UAF89uS1ypCYmKjFH7KzYcMGOnTo4KeIaj47vsZUjuW7jnLbm4sZENecCUMSKrS3kYgsV9WSxqWzy0fGGFPVHM/KZcyMlTSPCOeZmzpX6hAzdvnIGGOqmHGfrmVfWhaz7utNg/DKHarezhSMMaYK+WRlMp8m7WXMlW3p1jqy0vdvScEYY6qIXakn+H+frqNHbCMe+M2FfonBkoIxxlQBufkuxsxIQgReGppAcJB/hqq3NgVjjKkCXv5+C0l7jvHa7V2Ibljbb3HYmUIV9Omnn7J+/fqz/twbb7zB9OnTz3n/sbGxHD7sm2F5jTGnW7I9ldfnbWVwtxgGxLXwayx2plAFffrppwwYMKBwxFRv5OXlcf/99/swKmOMLxzLzOHhmUnERtXlyYH+H5rHkkIFGzt2LC1btuSBBx4A4Mknn6RevXqoKrNmzSI7O5vf/va3PPXUUwBMnz6dF154AREhLi6O3//+98yZM4effvqJZ599lo8++oj09HTuv/9+MjMzueCCC5gyZQqRkZH07duXhIQEFi5cyLBhw0hPT6devXrcfvvtXH/99YUxrVmzhu3bt1OnTh3uv/9+du/eDcCECRPo06cPqampDBs2jJSUFHr37k11vqHRmOpEVRn70RoOpWfz8f9cQt1a/q+S/R+BL301FvavqdhtnncxXDe+1MVDhgzhoYceKkwKs2bN4q9//Ss///wzv/76K6rKwIEDmT9/PlFRUTz77LMsWrSIxo0bc+TIERo1asTAgQMZMGAAt956K+AeLuPVV1/liiuuYNy4cTz11FNMmDABgJycHAru6n7yyScBaNGiReHw2q+//jo//fQTrVu35vbbb+fhhx/m0ksvZffu3VxzzTVs2LCBp556iksvvZRx48Yxd+5cJk+eXLHHzBhToplL9/D1uv2Mve4i4mIa+jscoKYnBT/o0qULBw8eZO/evRw6dIjIyEjWrFnDt99+S5cuXQDIyMhgy5YtrFq1isGDB9O4cWMAGjVqdNr20tLSOHbsGFdccQUAw4cPZ/DgwYXLhwwZUmosP//8M2+99RYLFy4E4Pvvvy/SVnH8+HEyMjKYP38+H3/8MQA33HADkZGV3zfamECz9WAGT32+nksuiGL0Zef7O5xCNTspnOEXvS8NHjyY2bNns3//foYMGcKuXbt47LHHuO+++4qs9+qrr57zvjyH3Pa0b98+Ro0axZw5cwqfzuZyuViyZAnh4eHnvF9jTPll5+UzZsZKwkODePG2BIL81P20JNb7yAeGDBnCjBkzmD17NoMHD+aaa65hypQpZGRkAJCSksLBgwfp168fH374IampqQAcOXIEKDrkdkREBJGRkSxYsACAd999t/CsoTS5ubkMHjyY5557jnbtTo2/3r9//yKJqOAS0+WXX174LOmvvvqKo0ePVsRhMMaU4oVvNrFu73GeuyWO8yKq1o80Swo+0KlTJ9LT04mOjqZ58+b079+f22+/nd69e3PxxRdz6623kp6eTqdOnXj88ce54ooriI+P55FHHgHcz1t+/vnn6dKlC9u2bWPatGk8+uijxMXFkZSUxLhx4864/0WLFrFs2TKeeOIJEhISSEhIYO/evbzyyissW7aMuLg4OnbsyBtvvAHAE088wfz58+nUqRMff/wxrVq18vkxMiZQzd98iLcW7ODOXq3o3+k8f4dzGp8OnS0iO4F0IB/IU9VEEWkEzARigZ3Abap6VNzDAL4MXA9kAiNUdcWZtm9DZ1c+O77GlN/hjGyunbCAyDqhzPnDpdQOC/ZLHP4eOvs3qprgEcBY4AdVbQv84MwDXAe0dV6jgYmVEJsxxlQKVeUvs1dzPCuXV4Z18VtCKIs/Lh8NAqY509OAmzzKp6vbEqChiDT3Q3zGGFPhpi/exX83HuSx6y6iQ/MG/g6nVL5OCgp8KyLLRWS0U9ZMVfc50/uBZs50NLDH47PJTlkRIjJaRJaJyLJDhw6VvFO7+con7LgaUz4b9x/nb19uoG/7Joy4JNbf4ZyRr5PCparaFfeloQdE5HLPhequZc6qplHVSaqaqKqJTZo0OW15eHg4qampVoFVMFUlNTXVurMac5aycvMZ80ESDcJDeWFwfKU+Ra08fHqfgqqmOO8HReQToAdwQESaq+o+5/LQQWf1FKClx8djnLKzEhMTQ3JyMqWdRZjyCw8PJyYmxt9hGFOt/OPLDWw6kM47I7vTuF4tf4dTJp8lBRGpCwSparoz3R94GpgDDAfGO++fOR+ZA/xBRGYAPYE0j8tMXgsNDaVNmzYV8ScYY8w5+WHDAaYt3sWoS9vQt31Tf4fjFV+eKTQDPnFOlUKA91X1axFZCswSkVHALuA2Z/0vcXdH3Yq7S+pIH8ZmjDE+dfB4Fo/OXk2H5g34y7Xt/R2O13yWFFR1OxBfQnkqcGUJ5Qo84Kt4jDGmsrhcyp8+XEVmTh6vDkugVkjV7H5aEruj2RhjKtjkhTtYsOUw4wZ04sKm9f0dzlmxpGCMMRVobUoa//xmI9d0asawHi3L/kAVY0nBGGMqSGZOHg9+sJKourUYf3Ncle9+WpKaPXS2McZUoqc/X8+O1BO8d09PIuuG+TuccrEzBWOMqQBfrtnHjKV7+P0VF3DJBY39HU65WVIwxphztPfYScZ+tJr4mAgevrpd2R+owiwpGGPMOch3KQ/NTCLfpbw8tAuhwdW7WrU2BWOMOQcT523l1x1H+NfgeGIbl/x43Oqkeqc0Y4zxoxW7j/LS91sYGN+Cm7ueNqhztWRJwRhjyiE9K5cxM1bSPCKcZ3/buVp2Py2JXT4yxphyGPfZOlKOnuTD+3vTIDzU3+FUGDtTMMaYs/TJymQ+WZnCmCvb0a11I3+HU6EsKRhjzFnYnZrJ//t0Hd1jI3ngNxf4O5wKZ0nBGGO8lJvv4sEZKxGBl4YkEFLNu5+WxNoUjDHGSy9/v4WkPcd47fYuxETW8Xc4PlHz0pwxxvjA3NX7eO3HrQzuFsOAuBb+DsdnLCkYY0wZVuw+yiOzkujWOpJnburs73B8ypKCMcacwZ4jmYyevoxmDcKZdFc3wkOrz1PUysPaFIwxphRpJ3MZ+c5ScvJczBjdnah6tfwdks9ZUjDGmBLk5rv4n/eWsyv1BNPv7smFTev5O6RKYUnBGGOKUVX+75O1/Lw1lRcGx9P7gih/h1RprE3BGGOKeeOn7cxctoc/9ruQW7vF+DucSmVJwRhjPMxdvY/nvt7IwPgWPFLNH5hTHpYUjDHGUdD1NLF1JP+8Na7GjHx6NspMCiLSR0TqOtN3isiLItLa96EZY0zl2XMkk3unLeO8iHDeDICup6Xx5kxhIpApIvHAn4BtwHRvdyAiwSKyUkS+cObbiMgvIrJVRGaKSJhTXsuZ3+osjz3rv8YYY8qhoOtpnkuZMiIwup6WxpukkKeqCgwCXlPV14H6Z7GPMcAGj/nngJdU9ULgKDDKKR8FHHXKX3LWM8YYn/LsevrGnd24oElgdD0tjTdJIV1EHgPuBOaKSBDg1RMlRCQGuAF425kXoB8w21llGnCTMz3ImcdZfqUE4gU9Y0yl8ex6Ov7muIDqeloab5LCECAbGKWq+4EY4Hkvtz8B+AvgcuajgGOqmufMJwMFDzaNBvYAOMvTnPWLEJHRIrJMRJYdOnTIyzCMMeZ0E3/axsxle3iw34XcEmBdT0tTZlJQ1f2q+qKqLnDmd6tqmW0KIjIAOKiqyysgTs94JqlqoqomNmnSpCI3bYwJIHNX7+OfX29iYHwLHg7Arqel8eUdzX2AgSJyPRAONABeBhqKSIhzNhADpDjrpwAtgWQRCQEigFQfxmeMCVArdh/l4QDveloan92noKqPqWqMqsYCQ4H/quodwI/Arc5qw4HPnOk5zjzO8v86DdzGGFNhCrqeNo8IZ9LvEgO262lp/HHz2l+BR0RkK+42g8lO+WQgyil/BBjrh9iMMTVY2slcRkz9tbDraaO6Yf4Oqcop8/KRiPQBngRaO+sLoKp6vrc7UdV5wDxnejvQo4R1soDB3m7TGGPORk6ei9//Zzm7j2Ty7qieAd/1tDTetClMBh4GlgP5vg3HGGMqnqryf5+uYdG2VP41OJ5e51vX09J4kxTSVPUrn0dijDE+MvGnbcxalmxdT73gTVL4UUSeBz7Gfb8CAKq6wmdRGWNMBfli9V7++fUmBiVY11NveJMUejrviR5livvOZGOMqbKW7zrKI7NWkdg6kudusa6n3igzKajqbyojEGOMqUi7UzMZPd26np6tUpOCiNypqv8RkUdKWq6qL/ouLGOMKb+0zFxGvuPuejrVup6elTOdKdR13s9mRFRjjPGrnDwXv3/vVNfT863r6VkpNSmo6pvO+1OVF44xxpSfdT09d/Y4TmNMjfHveU7X0yvbWtfTcrKkYIypEb5YvZfnv3G6nl7V1t/hVFuWFIwx1V5B19Pusdb19FyVmRREpJmITBaRr5z5jiIyqqzPGWNMZfDsevrmXdb19Fx5c6bwDvAN0MKZ3ww85KuAjDHGW9b1tOJ5kxQaq+osnEdqOg/HsYHxjDF+5dn1dNJd3azraQXxZpiLEyIShXtoC0SkF+7nJxtjjF+oKo9/4u56+uJt8fS0rqcVxpuk8Ajup6JdICI/A0049eQ0Y4ypdP+et40Pl7u7nt7c1bqeViRvxj5aISJXAO1xP2Bnk6rm+jwyY4wpweerrOupL3nz5LVg4Hog1lm/v4jY2EfGmEq3fNdR/vShu+vpP2+1rqe+4M3lo8+BLGANTmOzMcZUtt2pmdw7fRktnK6ntUKs66kveJMUYlQ1zueRGGNMKQq6nrpUmWJdT33Kmy6pX4lIf59HYowxJUg7mcvd05ay+0gmb95pXU99zZszhSXAJyISBOTibmxWVW3g08iMMQHvYHoWv5v8K9sOZfDK0C7W9bQSeJMUXgR6A2tUVX0cjzHGALDnSCZ3Tv6Fg8ezmTy8O5e3a+LvkAKCN0lhD7DWEoIxprJs3H+c303+lew8F+/d25OurSL9HVLA8CYpbAfmOQPiZRcUWpdUY4wvLN91lJFTf6V2WDAf3t+bds3s4Y+VyZuksMN5hTkvY4zxiZ82H+L+d5fTrEEt3h3Vk5aN6vg7pIDjzR3N5Xocp4iEA/OBWs5+ZqvqEyLSBpgBRAHLgbtUNUdEagHTgW5AKjBEVXeWZ9/GmOrn81V7eWRWEm2b1mfa3T1oUr+Wv0MKSKV2SRWRCc775yIyp/jLi21nA/1UNR5IAK51BtN7DnhJVS8EjgIFz2YYBRx1yl9y1jPGBIB3l+ziwRkr6dIykhn39bKE4EdnOlN413l/oTwbdhqmM5zZUOelQD/gdqd8GvAkMBEY5EwDzAZeExGxBm5jai5V5fUft/LCt5u58qKmvH5HV3tIjp+VeqagqsudyQRV/cnzhfuXf5lEJFhEkoCDwHfANuCY80wGgGQg2pmOxt3TqeCZDWm4LzEV3+ZoEVkmIssOHTrkTRjGmCrI5VKenbuBF77dzG+7RPPGXd0sIVQB3tzRPLyEshHebFxV81U1AYgBegAXeR9aqducpKqJqprYpIn1WzamOsrLd/Ho7NVMXriDkX1i+dfgeEKD7ZHxVUGpl49EZBjuyzxtirUh1AeOnM1OVPWYiPyI+ya4hiIS4pwNxAApzmopQEsgWURCgAjcDc7GmBokKzefP7y/ku83HOCRq9vxx34X2minVciZ2hQWAfuAxsC/PMrTgdVlbVhEmgC5TkKoDVyNu/H4R9wP6ZmB+yzkM+cjc5z5xc7y/1p7gjE1y/GsXO6dtoxfdx7hmUGduKt3rL9DMsWUmhRUdRewC/ev+/JoDkxznscQBMxS1S9EZD0wQ0SeBVYCk531JwPvishW3GciQ8u5X2NMFXQ4I5vhU35l0/50JgxJYFBCdNkfMpXOm5vXykVVVwNdSijfjrt9oXh5FjDYV/EYY/wn+Wgmd03+lX1pJ3l7eCJ92zf1d0imFD5LCsYYA7DlQDp3Tf6VzJw83runJ91aN/J3SOYMvGruF5HaItLe18EYY2qWpD3HGPzmYvJVmXlfb0sI1UCZSUFEbgSSgK+d+QQv72g2xgSwhVsOc/tbS2gQHspH919Ch+b2CJbqwJszhSdxtwEcA1DVJKCND2MyxlRzX63Zx93vLKVVozrMvr83raJsYLvqwps2hVxVTSvWj9i6ihpjSvTBr7t5/JM1dGkVyZTh3YmoE+rvkMxZ8CYprBOR24FgEWkLPIj7HgZjjCli4rxtPPf1Rq5o14SJd3alTpj1ZaluvLl89EegE+5RT9/HPSbRGF8GZYypXlSVf3y5gee+3sjA+Ba89btESwjVlDf/ajeo6uPA4wUFIjIY+NBnURljqo28fBf/+8kaZi1L5q5erXlqYCeCgmzYiurKmzOFx7wsM8YEmKzcfB54fwWzliXz4JVteXqQJYTq7kwD4l0HXA9Ei8grHosaAHklf8oYEygysvMYPX0Zi7al8sSNHRnZxzol1gRnuny0F1gGDMT92MwC6cDDvgzKGFO1pWZkM/Kdpazbe5yXhsTz2y4x/g7JVJAzDYi3ClglIu95PBTHGBPg9h47yZ2TfyHl6Ekm3dWNKzs083dIpgJ509C8RUROuy9BVc/3QTzGmCps68EMfjf5F9Kz8ph+dw96nn/awxFNNedNUkj0mA7HPZKpDWBiTIBZnXyMEVOXEiQw475edGoR4e+QjA+U2ftIVVM9XimqOgG4oRJiM8ZUEYu2HWbYpCXUDg3mw/svsYRQg5V5piAiXT1mg3CfOdhdKcYEiG/W7eeP768ktnEdpt/dk/Miwv0dkvEhbyp3z0dx5gE7gdt8Eo0xpspQVd5ZtJNnvlhPXExD3hnZnYZ1wvwdlvGxMpOCqv6mMgIxxlQdh9Kz+cvsVfy46RBXXtSUV4Z1oW4tu0AQCLy5fBQFPAFcint01IXA06qa6uPYjDF+8OPGgzw6exXpWXk8PagTd/VqTbFRkk0N5k3qnwHMB25x5u8AZgJX+SooY0zly8rNZ/xXG3ln0U4uOq8+79/bi3bN6vs7LFPJvEkKzVX1GY/5Z0VkiK8CMsZUvg37jjNmxko2H8hg1KVtePSa9oSHBvs7LOMH3iSFb0VkKDDLmb8V+MZ3IRljKovLpUxdtJPnvtpIRJ1Qpt3dgyvaNfF3WMaPvEkK9wIPAf9x5oOAEyJyH6Cqag9eNaYaOng8iz/PXs38zYe4qkMznrvlYqLq1fJ3WMbPvOl9ZBcVjalhvlt/gL9+tJrMnDz+9tvO3N6jlTUmG8DLm9BEJA6I9VxfVT/2UUzGGB85mZPPs3PX894vu+nYvAGvDEvgwqb2u8+c4k2X1ClAHLAOcDnFCpwxKYhIS2A60MxZf5KqviwijXD3XorFuRFOVY+K+2fKy7if4ZAJjFDVFeX4m4wxJVibksaYGSvZdugE911+Po/0b0etEGtMNkV5c6bQS1U7lmPbecCfVHWFiNQHlovId8AI4AdVHS8iY4GxwF+B64C2zqsnMNF5N8acA5dLeXvhdp7/ZhON6obx3j096XNhY3+HZaoob5LCYhHpqKrrz2bDqroP2OdMp4vIBiAaGAT0dVabBszDnRQGAdNVVYElItJQRJo72zHGlMP+tCz+9GESP29N5ZpOzRh/cxyRdW2oClM6b5LCdNyJYT+QDQjuXkdx3u5ERGKBLsAvQDOPin4/7stL4E4Yezw+luyUFUkKIjIaGA3QqlUrb0MwJuB8vXY/Yz9eTXaui/E3X8yQ7i2tMdmUyZukMBm4C1jDqTYFr4lIPeAj4CFVPe75pVRVLekBPmeiqpOASQCJiYln9VljAkFmTh5Pf76eGUv3cHF0BC8PTeD8JvX8HZapJrxJCodUdU55Ni4iobgTwnsevZUOFFwWEpHmwEGnPAVo6fHxGKfMGOOl1cnHeGhGEjtST/A/fS/goavaERZS5mNTjCnkTVJYKSLvA5/jvnwElN0l1elNNBnYoKoveiyaAwwHxjvvn3mU/0FEZuBuYE6z9gRjvJPvUt6cv40Xv91Mk/q1+ODeXvSyR2WacvAmKdTGnQz6e5SV2SUV6INz2UlEkpyy/8WdDGaJyChgF6eezfAl7u6oW3F3SR3pzR9gTKDbe+wkD89M4pcdR7jh4ub8/bcXE1En1N9hmWrKmzuay1U5q+pC3I3SJbmyhPUVeKA8+zImUM1dvY/HPl5Nvkt5YXA8t3SNtsZkc068uXktBngV9y9/gAXAGFVN9mVgxpjSZWTn8eScdcxenkxCy4ZMGJJAbOO6/g7L1ADeXD6aCrwPDHbm73TKrvZVUMaY0q3cfZSHZiax50gmD/a7kD9e2ZbQYGtMNhXDm6TQRFWnesy/IyIP+SogY0zJ8l3Kv3/cyoQftnBeg3Bm3teb7rGN/B2WqWG8SQqpInIn8IEzPwywR3EaU4mSj2by8Mwklu48ysD4FjxzU2cialtjsql43iSFu3G3KbyEu9fRIqxnkDGV5rOkFP7vk7UoMGFIAjd1ifZ3SKYG86b30S5gYCXEYozxcDwrlyc+W8cnK1Po1jqSCUMSaNmojr/DMjVcma1TIjJNRBp6zEc6w2kbY3zkp82HuP7lBcxZtZeHr2rHzNG9LCGYSuHN5aM4VT1WMOM8+6CLD2MyJmBtPZjB3+au58dNh2jTuC6z7utNt9aR/g7LBBBvkkKQiESq6lEA5yE5Xj2xzRjjnbTMXCb8sJl3F++idmgwj1/fgeGXxNq4RabSeVO5/wv30NkfOvODgb/5LiRjAkdevov3f93Ni99t5vjJXIb2aMUjV7ejcb1a/g7NBChvGpqni8gyoJ9TdPPZPnDHGHO6+ZsP8cwX69lyMIPe50cx7saOdGjewN9hmQDn1WUgJwlYIjCmAmw/lMHf5m7gh40HadWoDm/e1Y3+HZvZmEWmSrC2AWMqSdrJXF75YQvTFu0kPDSYsdddxMg+sdQKCfZ3aMYUsqRgjI/l5buYsXQPL363maOZOQxJbMmf+i0FqxYAABGLSURBVLenSX1rNzBVjyUFY3zo562Hefrz9Ww6kE6PNo0YN6AjnaMj/B2WMaWypGCMD+w8fIJn527g+w0HiImszcQ7unJt5/Os3cBUeZYUjKlAx7Nyee2/W5n68w7CgoP4y7XtubtPG8JDrd3AVA+WFIypAPkuZebSPfzr200cycxhcLcY/ty/PU0bhPs7NGPOiiUFY87R4m2pPP3FejbsO0732EjeGdCDi2Os3cBUT5YUjCmn3amZ/O3L9Xyz7gDRDWvz2u1duOHi5tZuYKo1SwrGnKX0rFxe+3ErUxfuJCRY+HP/dtxz2fnWbmBqBEsKxngp36XMXr6H57/ZzOGMbG7pGsNfrm1PM2s3MDWIJQVjvPDLdne7wbq9x+nWOpLJwxOJb9mw7A8aU81YUjDmDPYcyeQfX23gyzX7aRERzivDunBjnLUbmJrLkoIxJcjIzuPfP27l7YU7CBbhkavbce9l51M7zNoNTM3ms6TgPLJzAHBQVTs7ZY2AmUAssBO4zXmSmwAvA9cDmcAIVV3hq9iMKU1uvotPVqTw/LebOJSezc1donn02vY0j6jt79CMqRS+PFN4B3gNmO5RNhb4QVXHi8hYZ/6vwHVAW+fVE5jovBtTKY5l5vDBr3uYvngn+9Ky6NKqIZPu6kaXVvYoTBNYfJYUVHW+iMQWKx4E9HWmpwHzcCeFQcB0VVVgiYg0FJHmqrrPV/EZA+5nIk/9eQcfrUgmK9fFJRdE8cygzvS7qClBQdZuYAJPZbcpNPOo6PcDzZzpaGCPx3rJTpklBVPhVJUFWw4zeeEOftp8iLCQIG5KaMHIPm3syWcm4PmtoVlVVUT0bD8nIqOB0QCtWrWq8LhMzXUyJ59PVqYw9ecdbDmYQeN6tXjk6nbc3rOVPRPZGEdlJ4UDBZeFRKQ5cNApTwFaeqwX45SdRlUnAZMAEhMTzzqpmMCzPy2Ld5fs5P1fdnM0M5eOzRvwr8HxDIhvbk89M6aYyk4Kc4DhwHjn/TOP8j+IyAzcDcxp1p5gztWqPceY8vMO5q7eR74qV3doxqhL29CjTSO7z8CYUviyS+oHuBuVG4tIMvAE7mQwS0RGAbuA25zVv8TdHXUr7i6pI30Vl6nZ8vJdfLv+AJMX7mD5rqPUqxXC73rHMuKSWFpF1fF3eMZUeb7sfTSslEVXlrCuAg/4KhZT86WdzGXm0t1MW7SLlGMnadmoNuMGdGRwYgz1w0P9HZ4x1Ybd0WyqtR2HT/DOzzv4cHkymTn59GzTiHE3duSqDs0Iti6lxpw1Swqm2lFVFm1LZcrCHfx300FCg4K4Mb4FI/vE0jnaHm5jzLmwpGCqjazcfD5LSmHKwp1sOpBOVN0wHuzXljt6taJpfRu+2piKYEnBVHkHj2fxnyW7+M8vuzlyIoeLzqvPP2+NY2B8C3uwjTEVzJKCqbLWpqQxZeEOPl+9lzyXcuVFTbn70jb0Pj/KupQa4yOWFEyVku9Svlt/gCk/7+DXHUeoExbMHT1bM/ySWNo0ruvv8Iyp8SwpmCrhwPEs5iTtZdrinSQfPUl0w9r83w0dGJzYkoja1qXUmMpiScH4zbHMHL5au5/PklL4ZccRVCGxdSSPX9+Bqzs2IyQ4yN8hGhNwLCmYSnUiO4/vNxxgTtJe5m85RG6+cn7jujzYry03xrfgwqb1/B2iMQHNkoLxuey8fH7adIg5q/by/YYDZOW6aB4Rzsg+bRgY34JOLRpYw7ExVYQlBeMT+S5l8bZU5qxK4eu1+zmelUdknVBu6RrDwPgWdI9tZA+xMaYKsqRgKoyqsmL3MT5ftZcvVu/jcEY2dcOCuabTedyY0IJLL2xMqLUTGFOlWVIw52zj/uN8lrSXz1ftJfnoScJCgujXvikDE1rQ76KmdoOZMdWIJQVTLrtST/D5qr3MWbWXzQcyCA4S+lzYmIeuakf/Ts1oYCOTGlMtWVIwXjtwPIsvVu9jzqq9rNpzDIDusZE8M6gT113c3B5paUwNYEnBnFHBvQRzkvayZEcqqtCxeQPGXncRN8a3ILphbX+HaIypQJYUzGlKupegTeO6/LFfWwbavQTG1GiWFAxQ9F6CHzYc5GRuPuc1CGfEJbEMjI+mc7TdS2BMILCkEKDSs3JZt/c4a1PSWJ2cxrxNBwvvJbi5a7TdS2BMgLKkEACOZ+WyLsWdANakpLE2JY3th08ULm8eEc5VHZrZvQTGGEsKNc3xrFzWOhX/GicR7CiWADpHR3BTl2gujo6gc3QETepbryFjjJslhWos7WQu65xf/wVnADtTMwuXt3ASwM1doukcE8HF0RHWbdQYc0aWFKqJtMxc1u4tmgB2eSSA6Ia16RzdgFu7xdA52p0AoiwBGOMdVXDlQX4uuHIhP88978p1yvJOvRdf7spz5j2nPT+b77HMY75we3nFtn+Gac8YL30EOt1U4YfCkkIV5E0CuDg6gtsSW9I5OoLOLRpYAjBVg6q70srPdr/nZUN+jkdZTrHynDOU5UBesXWKV6ZFKuxiFfRp6xWv7D3W0/xKPlACwaEQFApBIRAc4n4PCi06XbisYDoUQmu750Pr+CQySwp+diwzh7Upxwsr/zUpaew+cioBxESeSgAFbQCN6ob5MWLjEy6Xu2Iq8ovwTL9W88pZSZ6hYiyxAnWW5ec6FbNnZV9CmSu34o9NUCgEh52qHAsr0+BT08UrzpBa3q1X+F6sAg52Pue5v9Iqbs91Pbdb6nQoBFXdzhyWFCqYy6UcO5lLakY2hzNyOHIih9QTBdPZpGbkkHoih9SMbFJP5HAs89R/opaN3AlgaA8nAbSIINISQFEul1NR5RSrrHKcCq3gV2XeuZUXeeV7Me/NOqXMaz6oq3KPowR5VGhnqhidCi2kFoTVheBICAlzKmmPV0ElXGqZx3Th5z3LCqZrOcucir2KV6A1UZVKCiJyLfAyEAy8rarj/RwSqsrxk3mknsguUpmnZhSdLqj8j5zIwaWnb0cEGtYOJapeLRrVDaP9efWJqluLFk5bQLkTgOqpCqbwl6ZTSWl+sWWlzBd8rvB6p8f1UVf+6dc0i/yi9HL9fI/lLs9fpZ7r5Z6qnF2eFb5Hhe3L03wp/ovS8xV8+rwEF50PqXWG9YuXBZe8Tc95CT69gvb8tenNL+fTfqkWzFtFa0pWZZKCiAQDrwNXA8nAUhGZo6rrK3pfJ3PyOZieVXLlnpFJesYJjmdkkpGZwcnMTIJcOYSRRxi5hJFLLckljDwiwlw0qwUX1VIahrmIaOoiItRFvRAX9ULyqRuUT+3gfGpLHrUkl6DC66TZ7ldqNhzMhY157l+Kp1XurtIr8IJ1KSED+UNhhVpQuXlWRp4VX+ipyq9g/ZCwoqfehb8uCyq4MB+Vey5z4jImwFWZpAD0ALaq6nYAEZkBDAIqPCksnv0SrTe+TRPyiHYq+FpOhR8ixU7jyxoBOtt5FSdB7lPhkDAICT81XeQ9HGrVP/Wr0PPXY+Gv0CCP6YJlQSX8Wg0q9rngousV30fx9Tx/UXq+il8vLa2StyEwjKkRqlJSiAb2eMwnAz2LryQio4HRAK1atSrXjtq2iYWMeELDwgmpFU5oeG3CwmsTElb79Mr7TBX6aWXOK7iWu+I0xphqptrVXKo6CZgEkJiYWK5rJy173wq9b63QuIwxpiaoSq1NKUBLj/kYp8wYY0wlqUpJYSnQVkTaiEgYMBSY4+eYjDEmoFSZy0eqmicifwC+wd0ldYqqrvNzWMYYE1CqTFIAUNUvgS/9HYcxxgSqqnT5yBhjjJ9ZUjDGGFPIkoIxxphClhSMMcYUEtUqMnZOOYjIIWBXseLGwGE/hFNV2fE4xY5FUXY8igqk49FaVZuUtKBaJ4WSiMgyVU30dxxVhR2PU+xYFGXHoyg7Hm52+cgYY0whSwrGGGMK1cSkMMnfAVQxdjxOsWNRlB2Poux4UAPbFIwxxpRfTTxTMMYYU06WFIwxxhSqUUlBRK4VkU0islVExvo7nsokIi1F5EcRWS8i60RkjFPeSES+E5Etznukv2OtTCISLCIrReQLZ76NiPzifEdmOsO013gi0lBEZovIRhHZICK9A/m7ISIPO/9P1orIByISHqjfjeJqTFIQkWDgdeA6oCMwTEQ6+jeqSpUH/ElVOwK9gAecv38s8IOqtgV+cOYDyRhgg8f8c8BLqnohcBQY5ZeoKt/LwNeqehEQj/uYBOR3Q0SigQeBRFXtjHuo/qEE7nejiBqTFIAewFZV3a6qOcAMYJCfY6o0qrpPVVc40+m4/9NH4z4G05zVpgE3+SfCyiciMcANwNvOvAD9gNnOKgFxPEQkArgcmAygqjmqeowA/m7gfmxAbREJAeoA+wjA70ZJalJSiAb2eMwnO2UBR0RigS7AL0AzVd3nLNoPNPNTWP4wAfgL4HLmo4BjqprnzAfKd6QNcAiY6lxKe1tE6hKg3w1VTQFeAHbjTgZpwHIC87txmpqUFAwgIvWAj4CHVPW45zJ19z8OiD7IIjIAOKiqy/0dSxUQAnQFJqpqF+AExS4VBdh3IxL3WVIboAVQF7jWr0FVITUpKaQALT3mY5yygCEiobgTwnuq+rFTfEBEmjvLmwMH/RVfJesDDBSRnbgvJfbDfV29oXPJAALnO5IMJKvqL878bNxJIlC/G1cBO1T1kKrmAh/j/r4E4nfjNDUpKSwF2jo9CMJwNxzN8XNMlca5Xj4Z2KCqL3osmgMMd6aHA59Vdmz+oKqPqWqMqsbi/i78V1XvAH4EbnVWC4jjoar7gT0i0t4puhJYT4B+N3BfNuolInWc/zcFxyPgvhslqVF3NIvI9bivIwcDU1T1b34OqdKIyKXAAmANp66h/y/udoVZQCvcw4zfpqpH/BKkn4hIX+DPqjpARM7HfebQCFgJ3Kmq2f6MrzKISALuBvcwYDswEvePwoD8bojIU8AQ3L32VgL34G5DCLjvRnE1KikYY4w5NzXp8pExxphzZEnBGGNMIUsKxhhjCllSMMYYU8iSgjHGmEKWFEzAEJEvRaThWawfKyJrKziGjIrcnrPNBKc7dsH8kyLy54rejwkMlhRMwFDV652B4GqaBOD6MtcyxguWFEyNICKPisiDzvRLIvJfZ7qfiLznTO8UkcbOGcAGEXnLGVP/WxGp7azTTURWicgq4AGP7YeLyFQRWeMMKvcbp3yuiMQ50ytFZJwz/bSI3OtFzEtFZLVzMxVlxNbdWTdJRJ53ngUQBjwNDHHKhzib7ygi80Rke8FxMcYblhRMTbEAuMyZTgTqOWNBXQbML2H9tsDrqtoJOAbc4pRPBf6oqvHF1n8A97hxFwPDgGkiEl6wX2d46jzcY+hwhv0CICL9nRh64P6l301ELvcitvtUNQHIxx1QDjAOmKmqCao601n3IuAaZ/tPOMfCmDJZUjA1xXLcFWsDIBtYjDs5XIa74i5uh6omeXw21mlvaKiqBZX5ux7rXwr8B0BVN+IeFqKds+3LcSeDubiTUR2gjapuOkO8/Z3XSmAF7kq8bRmx1VfVxU75+2c6GMBcVc1W1cO4B7oLiGGxzbkLKXsVY6o+Vc0VkR3ACGARsBr4DXAhRZ+8VsBzTJt8oHY5d70Ud/LZDnwHNAbuxV2Zn4kA/1DVN4sUup+FURGxFd+G/V83XrEzBVOTLAD+jPuyzQLgfmClejnAl9MIfcwZXBDgjmLbvgNARNrhHkRuk3P5Zg8wGPfZiWcMZ/INcLfz/AtEJFpEmpYRW7qI9HSKhnosTgfql/0XGlM2SwqmJlkANAcWq+oBIIuSLx2dyUjgdRFJwv1rvsC/gSARWQPMBEZ4jKC5APcDfU460zFl7VdVv8V9CWixs83ZlF2xjwLecmKri/uJYeAe8rljsYZmY8rFRkk1ppoQkXqqmuFMjwWaq+oYP4dlahi7zmhM9XGDiDyG+//tLtztJ8ZUKDtTMMYYU8jaFIwxxhSypGCMMaaQJQVjjDGFLCkYY4wpZEnBGGNMof8PdtkCvNveCgIAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# # Plot timed results against different observation set sizes\n",
        "\n",
        "# %%capture\n",
        "# key = PRNGKey(0)\n",
        "# window_size = 10\n",
        "# timestep_grid = jnp.arange(10, 500, 20)\n",
        "\n",
        "# it_times, vec_times = [], []\n",
        "# for timestep in timestep_grid:\n",
        "#     print(\"timestep \", timestep)\n",
        "#     data = choice(key, num_obs, shape=(timestep,))\n",
        "#     key, _ = split(key)\n",
        "#     it_time = %timeit -o get_fls_result_iterative(hmm, data, window_size)\n",
        "#     vec_time = %timeit -o get_fls_result_vectorized(hmm, data, window_size)\n",
        "#     it_times.append(1000*it_time.best)\n",
        "#     vec_times.append(1000*vec_time.best)\n",
        "\n",
        "# # Plot the result\n",
        "# plt.figure()\n",
        "# plt.plot(timestep_grid, it_times, label='iterative')\n",
        "# plt.plot(timestep_grid, vec_times, label='vectorized')\n",
        "# plt.xlabel('size of observation set')\n",
        "# plt.ylabel('best time in ms')\n",
        "# plt.legend();\n",
        "# plt.show()"
      ],
      "metadata": {
        "id": "a8iYj2-Fm_Zg"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}