synthetic_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from data_util import precompute_ar_covariates\n",
    "from twarhmm import TWARHMM, LinearRegressionObservations\n",
    "import matplotlib.pyplot as plt\n",
    "from util import random_rotation\n",
    "import numpy as np\n",
    "from scipy.linalg import logm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [],
   "source": [
    "# we'll generate some 2D autoregressive (1 lag) data\n",
    "data_dim = 2\n",
    "num_lags = 1\n",
    "\n",
    "# dictionary setting the parameters of the model\n",
    "true_config =  dict(\n",
    "    num_discrete_states=2,\n",
    "    data_dim=data_dim,\n",
    "    covariates_dim=3,\n",
    "    tau_scale=0.7,\n",
    "    num_taus=3,\n",
    "    kappa=5000,\n",
    "    alpha=2,\n",
    "    covariance_reg=1e-4)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [],
   "source": [
    "# initialize model using the above dictionary\n",
    "true_model = TWARHMM(true_config)\n",
    "\n",
    "# set transition matrices manually (don't have to do this)\n",
    "z_stay = .99\n",
    "transition_z = z_stay*np.eye(true_config['num_discrete_states']) + (1-z_stay)/(true_config['num_discrete_states']-1)*(np.ones((true_config['num_discrete_states'],true_config['num_discrete_states']))-np.eye(true_config['num_discrete_states']))\n",
    "tau_stay = 0.95\n",
    "transition_tau = tau_stay*np.eye(true_config['num_taus']) + (1-tau_stay)/(true_config['num_taus']-1)*(np.ones((true_config['num_taus'],true_config['num_taus']))-np.eye(true_config['num_taus']))\n",
    "true_model.transitions.transition_matrix = (transition_z, transition_tau)\n",
    "\n",
    "# set weight matrices manually (note that each matrix has an additional column associated with the bias (constant) term)\n",
    "true_model.observations.weights[0,:,:-1] = logm(0.92*random_rotation(data_dim,-np.pi/20))\n",
    "true_model.observations.weights[1,:,:-1] = logm(0.92*random_rotation(data_dim,np.pi/20))\n",
    "true_model.observations.weights[0,:,-1] = np.array([1,0])\n",
    "true_model.observations.weights[1,:,-1] = np.array([-1,0])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [],
   "source": [
    "# set the number of timesteps of data you'd like to generate\n",
    "T = 30000\n",
    "# sample for T timesteps, and give bias=True because the weight matrices have the additional bias column\n",
    "z,x = true_model.sample(T, bias=True)\n",
    "# also grab some test data\n",
    "z2,x2 = true_model.sample(T, bias=True)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "outputs": [
    {
     "data": {
      "text/plain": "[<matplotlib.lines.Line2D at 0x7ffb396a8580>,\n <matplotlib.lines.Line2D at 0x7ffb38e174f0>]"
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABB+UlEQVR4nO3dd3iUVfbA8e+d9B7SSC/0AEFKQhdQUMEKir2vbXd1V91dXV33t1Vdde27ri6KHcu6rg1FepXeISQkkEAaKSSk95n7++OGJqSQeSczk7mf58kTMjN53ztMcnLnvOeeK6SUaJqmac7LZO8BaJqmadbRgVzTNM3J6UCuaZrm5HQg1zRNc3I6kGuapjk5d3ucNCwsTCYmJtrj1JqmaU5r27ZtR6WU4T++3S6BPDExka1bt9rj1JqmaU5LCHH4bLfr1IqmaZqT04Fc0zTNyelArmma5uR0INc0TXNyOpBrmqY5OR3INU3TnJwO5JqmaU5OB3LNOI1VsOUtqD5i75Fomkuxy4IgrZdproNN/4YfXoHGStj3Fdz2NQhh75FpmkvQM3LtdIfXw4Jrob6ia48vy4JXR8PyP0PcOJj0IOSugd2f2nacmqadoGfk2klSwpL/g8KtsOi3cM2bHT++7igsmAvSDD9ZAvHjwGJRfwwW/w4GXgy+IT0zdk1zYXpGrp2Ut0EF8b4psOc/kLGw/ce2NMDHN0JtCdz4iQriACYTXP4yNFTC0j/0xKg1zeXpQN5bmFvUjPrHaorh6IGuHeOHV8E3FO78FiJTYOHDZ0+xtDbBlz+Dgi1w9TyITT39/sjhMPEB2PGBmp1rmmZTOpD3Bg3H4LVx8OktYDGfvL2+At6aAf8cA/MugM1vqvz1upfVY7/4qZo5A5Tth6xFkHYPeAfB7NfVcRc+BJX5KmXS2gxb5qucePoXcNFfYOhVZx/T1N+CfySse8nGT17TNENy5EKIYOAtYDgggZ9IKTcYcWytE1LC17+AihyoOAgrn4bp/6cC7xf3qdTHlEdh/yL47jcnv69PIlQVQtEOuPkzWP8PcPeGsfeo+yNTYNpvYcWTqgrFwxc8fKC+HGLHwlX/gP4Xtj8uTz8YdQuse1H9IQiOs+l/g6a5MqMudr4CfC+lnCuE8AR8DTqu1pmt8yHjGzU7Lj8Aa5+HqPPgaBZkL4FLn1fB+cInoHiPSrVEjwK/MMhdC5/erGbtDcdg1K3q9uPO/w0kng+l++BoNtSWwsgbof/0rpUWjr5VjWfHh3DB47b7P9A0Fyfk2fKq53IAIYKAnUA/2cWDpaamSr2xhAGK98KbF0LS+XDTZ2BpgXcvg5J0aG2EYXPgmvkdB93STFVuWF0AD2yF0P7GjvGDOapE8aHdYHIz9tia5mKEENuklKk/vt2IHHkSUAa8I4TYIYR4SwjhZ8BxtY7UV8Bnt4NPMMx+Q1WLuHvBdR+AVwCEDoArXul85hwxBO5bDXcvNz6IA4y+Xf2ROLDc+GNrmgYYE8jdgdHA61LKUUAd8NiPHySEuFcIsVUIsbWsrMyA07qwplo1i67Mh7nvgP8pW/gFRsHPN8I9K1RA7wrfEIgZbZuxDr4UfMNg+3u2Ob6maYYE8gKgQEq5qe3r/6IC+2mklPOklKlSytTw8DP2DtW6qrUZ/nMrFG2HuW9D4qQzH+Mb0vUgbmvunjDyJnWxtabY3qNxLulfwPb3IWsxFO1Ur72mnYXVFzullMVCiHwhxGAp5X5gOrDP+qFpZ2htgv/dCwdXwFWvQfLl9h5R14y+Hda/Cjs/gvN/Ze/ROIfyg/DZHaff5tMHhs9VF5yjR+teNtoJRlWt/AJY0FaxkgPcadBxteMq8+A/t6uZ+MVPqtI+ZxE2ACKGqcVBOpB3TfZS9fnORWDygKp8yFyoFllteVMF9KvnnbyAbDHDokfV51nPqXdCmsswJJBLKXcCZ1xJ1QxyYBl8frf6Jb1+gfPMxE8VNQIOrrT3KJxH9mIIHQgJE9XXcWkw/Gq1gGvDP2HN31X67PKXQFrgy5/D7k/UY6uL4Lr3VN2/5hJ00yxHV7wXFlwH4YPh+g9tU1nSEyJTYNfHqhbdP8Leo3FsTbVwaB2MvffM+3yC4cLfq5YMP7ysrofUlqogfsETah3Awl/Bh3Phxo/BO7Br56wtVat7EyapCijNqehA7uiW/p+aed3xrXN3EoxMUZ+L98CA6fYdi6PLXQPmZtU9sj0z/gQNFbD2BfX1lEdg6qPq316B6lrKR9fDHQs7rt8vzVQz/N3/AXMTDJoJc95Q+XhzK2x8DTa+AX0SIH48xI2HmDGnV0ppdqcDuSM7sExd2LzkaecO4gB9h6vPOpB3LnsxeAZA/IT2HyOE6jLp7gMBfWHyKdceUuaqGfuXP1VBetKDZz/G5jdV2wZ3bxh1MwTHw4qn4N9T4MI/qAvUxbvV6t6WetXGwdLWOycwVpWsXvCEWougnV1Trbq2kXwleNpuwbsO5I7KYoYlf1A9UdLutvdorOcbAkFxKpBr7ZNSXejsP63zC5YmN7j0ubPfd94NKoCseErNssMHn35/STosfgIGzIA588AvVN2eeL66qP6/u8G/L1z3vgpCQkBzverNU7RDXXQ/sFz9++7l6o+JdrrWZtWcLmel6ix63fvqwr8N6GSYo9r5EZSmq7fQ7l72Ho0xIlN0IO9MSTpUF3acVukKIdSFUE8/1XLY3HryvpZG+Lyty+Wcf58M4qBaEt+3Bi57Ae7frLpbHi9z9PRV6xYmPqDWMNz2lWqi9slNqj+9dpLFov7fc1bChAegpgjmTYP0L21yOh3Ie9iGg+Vc8/p6co/Wtf+g5nrVdTA2DYbO7rGx2VxkCpRnq+ennV32YvV5wEXWH8s/Ai79OxRugxV/gbpydfvyv6hJwlWvnd4k7Ti/UPUu0Ce44+NHj1QlkIVbVdWMxWL9mJ1ZbZla9FZbBosfh73/VROxS56C+9aqd0Wf3a66iRpMp1YMZrZIymqaKKysJ9zfm/jQk3kxi0Xy52/SySyu4db5m/j8ZxPpG+h94r4DZbUMjPBH5K6G2mKY/VrvWvQRmaJK5UozIHaMvUfjmLKWQOQI1WrBCMOvgcxv1cbYP7wCYYNUZ8y0e2CQlbN+gOQrVLBa9icIiFJBqzf9zHbV8r+cvPB83Pj7YdJD6t/BcWpNwKbXVarLYDqQG+iN1Qd5Ycl+WsyqCaSHm+DL+ycxLDoIgEV7i8ksruG+Kf34cONhbpu/mf/cN4HM4mqe/i6DXQVVnD8wjNdDl+Lv4avylb3JicqV3TqQn019BRRshvN/bdwxhYBr3oJxP4XDP6hFWUFxqu2xUSY9pGaiG19T1TaznnOtEsb0L1QQHzZH/c5Ki7omNHTO6X/U3D3bv/BsJR3IDZJfUc+LS7JISwzh0pQoIgK8+P2Xe3n40518/cBkPNxMvLQsi4ER/jw6cwhTBoVz5ztbmP7iao7WNhEV5M095yexYFMepXmLqI4YjVu9JNTPgrtbL/mlCE5QpXE6T352Gd+oIDDkMmOPa3JTe6rGj7PNylohYOYz4OahKlvMTaqixhXaFpdmwpf3qzTonHl2W1GrA7lBXlyahRDwwnXnERWkVtR5upu4450t/H3xfobHBHKgtJbXbhqNm0kwaUAYr9wwkqe+y+CRSwZz1+QkvD3cuHO4B9HvFPHXoguZ//RyTALCA7y4Y2IS907ph5vJid+2CqEveHZk7+cQ0g+iRtp7JOdOCLjor6qUcc3fVbMvk7v6CIhSK3sjR8CI651vt6hVz6prCte8DW6nhMzGKrUxi6efqkixY1sEHcgNkF5UxZc7C7lvSv8TQRxg2uAIbpuQwPx1uYT5ezIkMoBZwyNP3D8rJYpZKafnQqPLNwIw86qbSLTEUlbdyK6CKp79PpPlGSW8cN15JIQ6cbv3yBTY/oEqr3SFGVtX1ZbCobUqreKsOWYh1KrTiKFQlgmWVlXPXnkYjuxSF/nSv4CfrnOe51iRC2ueU88l4nmY1tah29wK/70Ljh2C27+BwGi7DlMHcgM8syiTIB8PfjbtzOXzj89KZt2Bo+SU1fH0nBRMnc2oc1aCfyRpaRNJa/thl1Ly1c4i/u+rvcx6ZS1v35HG+H6hHR/HUUWmQEud+gWxUU2tU9r3lUqrDL/G3iOx3vCrz3771rdh4cOqiibWSVozrX5OvasYcBGsfhb6TYO4capB2YGlKoV0vB+OHfWS5Kv9rMs+ytrsozxwwQCCfDzOuN/H0423bkvlj1cM5aKhnSyasFggZxX0v+C0GYsQgtmjYlj80BR8Pd34YONhg59FD4ocoT4X77bvOBzN3s/VTDYi2d4jsZ3hc9Um3s6yyUhZluphk3a3KrMMjlf196ufVXvlTnoQUh2j0asO5FZ6dXk2McE+3Dohod3H9Av3585JSYjO3k4W71YLLPpdcNa7o4N9mDY4grVZZbSanbRmN3yIast6ZKe9R+I4qgogb0P7M9newjtQPcc9n0NTjb1H07lVf1MtECY/rMZ+9Vtqsdaqv6n1HdP/ZO8RnqADuRWySmrYfKiC2ycm4OVuQL43p63Na79p7T5k2uBwqhtb2Zlfaf357MHdU7013fAarHn+9BWHrir9C/V5WC8P5ACj71Cptb2f23skHSveC+n/g/E/O7loKi5NLbAaNkc1FnOgEkvHGYkT+mhTHp5uJuaOMegq/MGVagOGDvpWnD8gHDeTYOX+UmPOaQ/Xf6D6d6z4K7x7KZScZUOp6iJ1Iak3ammE9f9UwawsS32OHuW8LYrPRWwqhCerqhZHtvYF8ApS7QhOlXYXXPuuw/V61xc7u6mh2czn2wuYlRJJiJ8BZUfN9ert9dl6UJ8iyNeD0fHBrNpfxiOXOGnXOd8QuPYdVS+98Ffw+gT1B2zoVaqfx76v1cIYD1+47Ws1E+pNMhfCkidOv+3iJ+0zlp4mBIy5Hb5/TM16I4fbe0RnMreqzqPD5qh2vk5Az8i76ZvdRdQ0tnLzuPZz4+dk31dqVVwXWrxOGxxBelE1pTWNxpzbXlLmwi+2wsxnVQOnVX+DJb+H1gbVHtW/L3x0rVp00Zvkt/2Runc1zH4dpv0ORt9m71H1nBHXg5snbHvX3iM5u6Id0FStig6chJ6Rd9NHm/IYEOFPWqIBf7Gb61WaIXo0JE3r9OFTB4Xz98X7Wb2/jGtTnWxxxY/5R8D4n6qPmhK1KjA4Xt2Xci28fQl8eDX8ZLHzLSRpT8FmtTlD9Ej14Wp8QyDlOtj2Dpx3o+O1a8hdpT4nTrHrMM6FnpF3Q3pRFTvzK7l5XHznlShdsfE1dTX8kqe6dAFlWHQg4QFerMoqs/7cjiSg78kgDhCSBLf8TzXnXzBX5Za7Q0rVkS5/C2QstG/3xeZ6tbI1tpeli87VJU+qFZ+f/wQaq+09mtPlrFZlsn7Os1ZDB/JzVFTZwPOL9+PtYeLqUbHWH7CmBNa9DEMu7/LCAiEE0waFO3cZYldFDle9r8syYcM/zry/o6qX5np10eq5fvD8AJg/Qy2p3jrfduPtTNEOtUowbqz9xuAIfPrANfOhMl8tEpLS3iNSmushfxP0m2rvkZwTHci7KONINfd/tJ3zn1vJ6qwytQDI98wFQOds1dPQ2ggz/nxO3zZtcATVja3scNYyxHMxcIaqclnzAlTmnbx91bPwXJIKjqeyWFT+9dVRqr1obJrKw9/4idpx6fD6nhz96Qo2q8+uPiMH1cTrgra+3TsX2Hs0St4Gda2qgxJgR6QDeRfUNLZw05sbWZtVxt2Tk1jz6AU8cOFA6w9clqXKsNLuPufl6pMHhuFmEny7+4j143AGlzytPi9uq/bY9G/1R7C5TvW8OHWBybI/wjcPqg2D7/webv6PysEPngWJkyFvo/02QcjfohpjnW1DB1c0+Veq9ev3v4OGSnuPBnJXqwVrHe2X6oB0IO+C99Yf4lh9Cx/ePY7HL00mto9Bm6hmtrUtPXXj3C4K8vHg6lExLNh0mINltcaMx5EFx8GUX0PG1/Dtb1SviyGXw61fwLFc+K5tB/mNr6tNg1PvUhdIE370Cxk3Xu0+X57d889BSjUjj3XxtMqpTG4w82/QVKX+ONtbzmq1YM3TuRrT6UDeierGFuatyWFGcl9GxAYbe/DcNZ0uAOrIozOH4O3uxl++2Yd0lByjLU34BfRJgi1vqlncNfNVLnPKI7DrI/jqAfj+cRXgL/372TvsHZ9p5W3s2bGDWuBUV9b76uKtFZmiXrONr6nWsPZSX6G6NDpZfhx0IO/U2+tyqW5s5aEZBqRSTtXapIJJUvdLnMIDvHhwxkBWZ5WxItOJV3p2lYe3al406ha44SP1NcCUR9VMe8cHajZ1zVvtt8gN7Q++YfYJ5AVb1Gc9Iz/TlEdUEN80z35jOLQWkE6XHwcdyDtUVd/C/LW5zBwWyfCYIGMPXrBVXeRMsm47t9smJNI/3I+/LNxHU6vZoME5sLixatNg78CTt7m5q5Wikx+GGz/uePm0EBA/Xl3U6mn5m8HDT3U51E4XPRIGzYIN/7RfOWLOKvAMUOs5nIxhgVwI4SaE2CGEWGjUMe3trXU51DS18qDRs3FQaRVhgoRJVh3G093EH68YxuHyet5ed8iYsTmjwGi1CbBvSOePjR+v8uo1JTYf1mkKNkPM6NN3mdFOmvZbaKyEzXaYlbc2wf5FamLlhK+PkTPyB4EMA49nV0WVDcxfl8tlKVEkRwV2/g3nKncNRJ0HPsFWH2rKoHCmDQ7nzbU5NDS7wKzcWsfz5Pk9mF5prlO9RVy9frwj0aNg4CWqM2Z3F3911+5PoeYIjL2nZ89rEEMCuRAiFrgMeMuI4zmCP3+TjkVKHptlg8ZUzfUqX5poXVrlVPdfMICKumY+3ZLX+YNdXeQI1We6J/PkRTtAmnV+vDPjf6aqijJ78I29xawW5UWNbHcvAEdn1Iz8ZeBRoN3iXCHEvUKIrUKIrWVljr20fHlGCYvTS/jl9IHEhRhUaniq/I1gaYEk466OpyWGkJrQhzfX5tLS21d7WsvdU/U66ck8+ZFd6nOM8+Vfe1TSVAhO6NmGWhlfQ8VBdY3FWfYS/RGrA7kQ4nKgVEq5raPHSSnnSSlTpZSp4eHh1p7WZhqazfzhq3QGRvhz9+R+tjlJ7hq1D2D8eEMP+/ML+lNY2cBXO4sMPW6vFD8ejuxWfVx6Qsk+8ItQTcK09plMqhPkobVQftD255MS1r4IoQMg+Qrbn89GjJiRTwKuFEIcAj4BLhRCfGjAce3i1RXZFFY28OTs4Xi626ioJ3etmhF6+Rt62AsGRzAkMoA3Vh/EYnGBunJrxE9QqY7CDucfxinZC311tUqXjLoFhFvP7O15cLnaYnHSQ+2XrDoBqyOVlPJxKWWslDIRuAFYIaW8xeqR2UFtUyvz1+Zy9egYxtlql/rGKijablX9eHuEEPxsWn8OlNayNKOHKzKcTVyaqho6sNT257KYVdOvvg64iYIjCohU7RR2fgStzcYeu+4ofH43vDVDfXzxMwiIVj3SnZiuIz/F+gNHaTZbuNaordvO5tAPalm+DQI5wGUpUcSH+PLikiydK++Id5BqxLX9fdunVypy1JoBXT/edaNvV6tgsxYZd8zCbfDvqZDxDXj6g1eA6q556XPquokTMzSQSylXSSkvN/KYPWl1Vhl+nm6MSbDh9k4b/6V2vokbZ5PDu7uZ+P1lyewvqeGttbk2OUevMeEB9Q5p18e2PU/JXvW57zDbnqc3GTAdAmNhm0HplR0L4O2ZKgd/1xK47UvVp+fWL5w6N36cnpG3kVKyOquMCf3DbJcbP7ROXcSZ9BC4e9nmHMDFwyK5eGhfXlmeRV65HTdRcHRxaaoccOO/VPrDVkr2qTROuJPusWoPJje1t+fB5bD/e+uOVXcUvn5ATZ7uXa3Wb/QyOpC3yT1aR8GxBqYOtmFFzapn1Gw89U7bnaPNn64chpsQ/P6rva7RUKu7JtyvUh9ZVgaLjpSkq6qI471htK6Z+EvVUOuL+07vQ3+u8jerdOYFv+vayl8npAN5mzVt26ZNHWijQH7oh5Oz8Y56gRgkOtiHX188mDVZZby5Nof9xTXUNnWwm46rGnK52l5uw2u2O0dpuk6rdIeHN1z7nnq39Nmd3b/wWbBFlftGjzJ2fA5EB/I2q7PKSArzIz7UBguAAFb33Gz8uNsnJjIyLpinv8vkkpfXMPyPi7nznc16hn4qN3cY9zM4/AMUbjf++E01qn1thA7k3RLaH676JxRuhWV/6t4xCraomX0HEyizRbLhYDnbDh9zytJd5+sOYwONLWY25lRwXaoBe3CezeH1ahHQJX/rkdn4cW4mwSf3jie9qJqiygY25JTz0aY8VmeVMW2wXphywqhbYOVTald3o1delra1H9Iz8u4bNhsO3aP6lY++DSLO4VqDuRUKt1GdfAM+ZgsebqfPXfPK6/locx5f7iikuFr1d4kO8uayEVGMiA1Goq6fDYkMZHBkgHHPyWA6kANbDx2jocVsm/y4lLDszz0+Gz/O20NV4YxJ6MMlwyJZlVnKP1ccYOqgcISTLkc2nHcgDJgBWUvUFnAmA9+olqSrzzqQW2fqb9Uf2h0fwCVPdf37StOhpZ4/bvdDtO7mxetGnrirsLKBy/+xlrpmM9MGhfP7y5NpNUu+2VXEu+sP0WI+OTOPCPBi4+PTMZkc83dGB3JgTXYZnm4mxttiEVD2EtVb5bIXe3Q2fjae7ibum9qfP36dzqbcCts8X2c1eBbs+xKO7DR2Vl6SrnpcB8cbd0xX5B8Og2bCrk9g+h+7Xvedrza73mIeQMH2Qi5K7suslCgsFslv/rMLs0Wy5OEp9A8/ucp69qgYqhpaKK5qxM0E67KP8qdv9rEt7xhpiY55sdTlc+Rmi2TV/lLSkvrg62nw3zWLRc3G+ySpt4QO4Pq0OML8vXht5QF7D8WxDLgIEJC12Njjlu5TS/P1ux/rjb4N6o+eW4VRwRaq3UOo8oxkRGwQv/tiD6XVjbz9Qy4bcsr5wxVDTwvixwX5eDA4MoABEQHMTY3D093k0Budu2wgzymr5dnvM5n0zAqySmq5eGik8SfZ+1/11u7C34Obh/HH7wZvDzfuOT+JtdlH2Zlfae/hOA6/UNUr3MgyRCnVYiC9otMY/adDQJRKr3RV/mZ2yoFM6B/Gi9eNpL7ZzH0fbuO5xfuZkdyX61I7X8Xt7+XO1EHhfL+32GEvhLpkIN9bWMUlL69h3pochkYH8tpNo7l1fIKxJ2lthhVPqqvlw6429thWunl8AkE+Hjz17T72FlbpKpbjBs1UqZVqg2Ze1UVq5ajOjxvDzR1G3gQHlqn/287UHYVjufzQ2I/JA8MYEOHP7y5NZkdeJYHe7jxzTUqXrxNdmhJJcXUjOxx08uNygVxKyZPf7iPA24Mffnshb9+RxmUjooy/iLFzAVQehul/MvbimQH8vdz5zSWD2Z5XyeX/WMfUv6/ixaVZus580Ez1Odug9MqJC526WZZhRt2iFvfs/Eh9LaX6Y3k2bZtdb7cMZNKAMABuHZ/Ary4axL9vHUOYf9dXV09P7ounm4lFe7r/R76uqZVHPttFTpnxvX0cK8L0gOUZpWzMqeChGQOJDLLhSrvMhRA2SPWMcEC3jk9g6xMzePaaFBLD/Hh1eTYzXljNd3uOuO4MPSIZguKtXxJ+XO5qtRBFt681Tkg/tbPW1nfgf/fCi0PhmQRVcfRj+Zsx40ZZQDL9wvwAMJkEv5w+kDEJ53bRMtDbg/MHhrFob3G3fz8+2HiYz7YVUNnQ0q3v74hLBfIWs4WnF2XQL9yPG8fasIqgtVnVjveb5tAXufr4eXJ9Wjzv/2Qs//v5REL8PPn5gu3cOn8zK/eXYnbQfKDNCAGDZ6rd1FsarDuWuRV2/0fN8r2DDBme1ib1J1BdAAeWQ/w4tWho4cNq8dUpZP5mMkkkbWCMIaW2s1KiKKxsYFdBO+8AOlDX1Mq8NTlMGRTO6Hjjm/K5VCD/ZHMeOWV1PD4r+YyFAYYq3AYt9TZrVWsLo+P78PUDk/jD5UPJLK7mzne2MOW5lbyx+qBrBfRBl0Brg1rAZY2Dy6GuVOV0NWMNmwO/yoRHDsC178Ls16G6EJb/9eRjmuuxFG5nc+sAJg8MM+S0FyX3xcNNdJpeaTFbyCyuPu22DzYepqKumQenDzRkLD/mMoG8tqmVl5ZlM75fCDOSbbyqMXcNICBxsm3PYzB3NxM/mZzE+sem89pNo0kI9eWZRZm8sboHttxyFAmTwcMP9n9n3XF2LgDfMBh4sTHj0k4SAgKjTr7bjRsLY++BzfMgfwsU74U3L8CttZ4VllFM7G9MIA/y9WDSgDA+21bAtsMVZ31MTWMLd76zhZkvr+X5xfuRUp6YjZ8/MMxmLbJdJpCvP3CUirpmfnnhQNuvaMxta5XpY8O+5jbk6W7ishFRLLh7HJePiOKlpVnsctCr9Ybz8FbplX1fg7mbucz6Cti/CEZc5zBlp73e9D9AYDR8dju8eSE0HOOZsKcpi5hEeIBxLaMfmzUEPy83rvv3Rv6xPPu0d6sl1Y1c+8YGNuaUM3lAGP9ceYA/f7OP9zeo2fhDM2wzGwcXCuS7C6pwMwlG2SA/dZrmerWarN9U256nBwgheGpOCn0DvXnwkx2uU9WSci00VMDBFd37/r2fg7lZp1V6klcAXP6SKkvsfwGNd6/l7eJ+TB5gzGz8uCGRgXz7y/O5LCWKF5ZmMfPlNdzz/lYe+WwXc177gfyKet6+I40P7hrL3ZOTeHf9IZ5bnNk2G7fdqlCXWaK/q6CSwX0D8PG08Qar+RvB0uJU+fGOBPl48NL1I7lh3gae+GIPs0fFUFzVSEVdM0MiA0hNDCHIp5fNOvtPB+9g2POZypmfq50L1PqByBTDh6Z1YNAl8OtM8O/LtoPlNLdaTpQdGinQ24NXbhjJtMHhfLa1gPyKevY2tODn5c6821IZHqMubj9xWTKBPh78a9UBHr5okOHjOJVLBHIpJbsLqrg0xQarN38sZzWYPNQu7b3E2KQQHrhgAK+uOMBXO09fiGESMDwmiGeuHsHQ6EA7jdBg7p6q497u/0BzHXj6df17S/ZB0Q6Y+YzNhqd1IED9jq8/eBQ3kyAtyTazYCEEV4+O5erR7XdMFUKVOv50an/b7TrWxiUC+eHyeqoaWjgvNtj2J8tdA7Fp5/bL7wQemjGIiQPUNnhRQd4EeHuwp6CKjTnlLNiUx4Of7OCbX0zG28PG73isJKXs2jWSlGth27sq150yt+snSP9CbeuWcm23x6hZb8PBcs6LDcLfy/4hztZBHFwkR76roBKAEbYO5A2Vaol3L0mrnMpkEozvF8ro+D5EBfng7+XOhP6hPHzRIJ6/dgTZpbW8tCzL3sPs0Le7j5D65DK+68rqvPiJEBij0ivnIm+DSqn4Gf+WXuua2qZWdhVUGVat4gxcI5DnV+HtYWJQ3zO7nBnq8A9q+XAvuNB5LqYNjuDGsXG8uSaHbYeP2Xs4Z2gxW/jrwn3c/9F2yuuau1ZOaTLB8GtUX4/6s5eancHcotYQxI23bsCaVbbkVmC2SCb0d502za4RyAsqGR4dhLstFwGBWmnm7gMxqbY9jwP63aXJRAX58Mhnu2hotuGO9Ocoq6SGm97cyPx1udwxMZHfX5bM7oIqdre9S+tQyrVgaVV9yruieI9aCBY/zpoha1Zaf/Aonm4mm9VsO6JeH8hbzBbSi6o4Ly7YtidqrFYXx5Kv6HrT+14kwNuDv88dQc7ROl52gBRLSXUjj32+m5kvryGzuIZXbhjJn64cxnVpcfh4uPHRpi7syh6ZovrlZHzTtZO2bWJAnA7k9rQhp5zRCcEOf73GSL0+kGeV1NDYYmFErI37Xez8CJprYPxPbXseBzZxQBjXpcYyf10u+4trOv8GG9mUU84Fz6/i8+0F3DExiTWPXMBVI2MAVTp21chovtpZRHVjJwt+hICEiSpd0pVGSfkbITAWgmy096vWqcr6ZtKLql0qPw4uEMh3tzW4sWnFisUCm96A2LEQM8Z253ECj81KJsDbnd9/uccuTfgPl9fx0w+3ERXkzbJfTeUPVwylj9/p75BuHpdAQ4uZL3cUdn7A6NGqTWpFTsePkxLyNum0ip1tzKlASlwqPw4uEMh35VcS5ONBQqiv7U6SvRiO5br0bPy4ED9PHp+VzJZDx/jv9oIePXd1Ywt3vbcVCcy/PY2E0LOXgKbEBjEiNogFG/M6b0l6fP/Owu0dP66qAGqKdFrFzjYcPIqPh1vPlBo7EKsDuRAiTgixUgixTwiRLoR40IiBGWVXQRUjYoNs219l4+uqVC35Studw4nMHRNLakIf/vZdBsfqmm1+PiklhZUN/OKjHRw6WsfrN48hMazjOv6bx8Wzv6SGrZ1V2YQnqwvYRZ0E8vxN6rMO5Ha1/mA5aUkhPVK77UiMeLatwK+llEOB8cD9QgiH6KTf0Gwmq6SGkba80FmyTzXJSrtbN0hqYzKpHi3Vja28/UOuzc6zt7CKW97axKi/LmXSMytYnVXGX2cP79Lb6ivOiybAy52PN3dy0dPNXTVAK9zW8ePyN6muiXo3ILspq2kiu7SWiS6WVgEDVnZKKY8AR9r+XSOEyABigH3WHttae4uqMFukbRcCbf63mrGNucN253BCgyMDGBEbxIaD5TY5fk5ZLbe9vRl3k2DW8CiGRgcyKi74RJ+Lzvh6unPlyGg+317An64cRqB3B3+EY0arHWnMrSqwn03eRogd0/79ms2lF6nrYaNsXaHmgAx9/yGESARGAZvOct+9QoitQoitZWVlRp62Xav3l6l+C4k2qidtbYb0L2HoleBru85mzmpsUgi7CippbDG2rry4qpFb529GAJ/eN4G/XZ3CreMTuhzEj7s+LY7GFgtf7+xkI9/o0WqzibKMs9/fVAMle/VCIDs7UKr2whzUN8DOI+l5hgVyIYQ/8DnwkJSy+sf3SynnSSlTpZSp4eHhRp22Q8sySkhN6EOwr43qunNXQ2MlDLvaNsd3cmMTQ2gxS3Ya2Mu8sr6Z29/eTGV9M+/eOZakTnLhHUmJCWJIZACfbsnv+IEnLni2k14p3KZW9OqKFbvKKqkhzN/zjColV2BIIBdCeKCC+AIp5f+MOKa18ivqySyu4aKhfW13kvQvwCsI+l9gu3M4sdSEEISAzbldXOLeiW2HK7js1XXkHq1j3m2ppFi5NkAIwQ1pcewprGJf0Rlzj5NC+qm2tu1VruRtAoRqlqbZTXZpLQMibNyGw0EZUbUigPlAhpTyReuHZIxlGSUATE+2USBvbYaMhTDkMnA3bgeS3iTI14MhkYFWB3KzRfLPFdlc9++NmEzw6X3jDeszPXtUDJ7uJv6ztYNZuRAQPar9ypUDSyF6pN5k2Y6klBwoqXXJtAoYMyOfBNwKXCiE2Nn2cakBx7XK8oxSBkT4W/XWu0M5K6GpSm0Eq7VrXFII2w4fo8VsOefvLaxs4NXl2Ux7fiXPL8ni0pQovv3l+Ybu8hTs68klwyL5Ykdhx7n8mNGqQqm5/vTba0uhYCsMtvuPvEsrqW6ipqmVgS46IzeiamUdYONNMM9NdWMLG3PKuev8JNudJP0LNQPrN8125+gF0hJDeHf9IdKLqrtUBlrV0MLivcV8ubOQDTnlSAmTBoTyu1nJzBweaZP1ADekxfHNriIWpxefWMp/hpgxIM2qMdapufCsxYCEwbMMH5fWdVklqiXEgAjXnJH3ylqp1fvLaLVILrJZWqUJMr9VC4BcsEHWuUhLUrPnzbnlHQby5lYLT3+XwUeb8mg2W0gM9eWXFw5k7phY4kJsuCoXmNAvlPAAL1ZmlrYfyKPbLngWbT89kO9fBEFxun7czrLbKlYG2rpVtYPqlYF8WUYJIX6ettto+eAKaKqG4Tqt0pmIAG/6hfmxObeCe6f0P+tjjtY28fMPt7P5UAU3pMVx49h426/GPYXJJDgvNog9hVXtPygwCgKiTq9caWlQPwujblF5dM1uDpTWEOLnSZi/a16v6nWBvMVsYWVmKRcPi8TNZINfroZKtSTfpw8kudYGEt2VlhjC9+nFWCwS049ek31F1dzz/laO1jbxyg0j258R21hKTDDLM0upbWptf3uwpCnqAvexw9AnQW3r19qg0yoOILvEdStWoBc2zdpyqILqxlZmJEcYe2ApYe/n8NpYOLQWpjyql+R30dikEKoaWsgqPb21bUl1I7e9vQmzRfLZTyfYLYgDpMQGIiUdlyFe+H9q5r3wYfXzsP878PSHxMk9N1DtDFJKsktrXfZCJ/TCGfmnW/IJ8HJnyqAuLjo6dhg2vAYe3qpBUkQyRAw9mfuWUr19XveSCuBRI+GmT1U5mtYlY9t2Mt9wsJwhkYEAtJot/OLjHdQ1mfn6gfEMtHPZ2PFVoXsKq06M9wzBcTD9j7DoEdj1Cez/HgZM1+WndlZW00RVQ4sO5L1FaU0j3+05wi3jE/D17OSpNdfBupdh/attmwZIMLd16nPzUnXB0aMgdy2UpoN/JMx6TjXHMrnOziNGiO3jw+C+AfxtUSZmi+Qnk5J4YWkWm3MreOn68+wexEHl8vsGerG3ozw5QNpdakPmhQ9Ba6MuO3QAJy902v/nyF56VSD/ZHM+LWbJreMTOn5gTTG8dRFU5cHwuXDRX8A/Qm0eUJKuLmgVbFGNkkIHwOzX1eN0hUq3CCFYcM84Hvt8D09+m8HXu4rYXVDFjWPjmTPKcXbTSYkJ7nwvT5MbXPkqvHE+CBMMvLhHxqa1L7ut9NBVK1agFwXyFrOFBZsOM2VQOP3CO3lBF/8Oakvgjm9Pz2+GD1Yfw9t6p1gsajd1zWph/l68edsYPt2Sz18W7mNYdCB/vMIhuh2fkBITxPLMko4veIJKv816Fo4d0s3SHEB2aS1BPh6Eu2jFCvSiQL4kvYSS6iaentPJbPzAcnXRctrvOr9IpYO4oYQQ3DA2nouHReLlbnK4zXFPveDZbp78uLS7emZQWqeyS9SFzp4qV3VEvSZSvbfhEPEhvkwb3EG1SksDfPtrlS6Z/FCPjU07XYifJ34dzXjt5NQLnppzkFKSVVrj0mkV6CWBPONINZtzK7h1fELHteNrX1R7a172oq400M4QEeBNZKA3ezrLk2sOo7yumcr6Fpddmn+c0wdyKSXPfp+Jn6cb16Z2cOGs/CD88DKMuB766YU82tkNj+lkhafmUI7X/Q/SM3Ln9v3eYlbtL+PhiwZ1vIHEyqfA5A4X/bXnBqc5nZSYIHKO1lHb1GrvoWhdsCKzFC93E2MSbNSOw0k4dSCvbWrlz9/sY0hkAHdMTGz/gUd2qQuc438OATbcaEJzel1a4ak5BCklS/eVcP7A8M7XjfRyTh3IX1mWRXF1I0/NGY67WwdPZflfVW+USb/sucFpTklf8HQe6UXVFFY2cLEtdwFzEk4byDOLq3n7h0PckBbHmIQOSsUOrVM7uEx+WO/gonXq+ArPdB3IHd6SfSWYBEw3uq+SE3LaQP7e+kN4u5v47cwh7T9ISlj2Z9V+dOy9PTc4zan1D/cn52idvYehdWLpvhJSE0IIdeGFQMc5ZSCXUrJ6fxmTB4a1v2P28Zrxgs0w9bfg4dOzg9ScVmKYH4fKdSB3ZPkV9WQcqebiYTqtAk4ayLNLaymqamx/8U/xHpg3DbbOhwkPwOjbenR8mnPrF+ZHZX0LlfXN9h6K1o4l+9Tm6hfp/DjgpIF81f5SAKb+uFWtlLBlPrx5ITQcg1v+B5c8pbsVauckMVRt2J2r0ysOa+m+Ygb3DSAh1EabqzsZpwzkq7PKGNTXn+jgU9IlLY3w1QPw7a/Uzj0/26B6RWvaOUoMU8FBp1cc07G6ZjbnVui0yimcrviyrqmVLbnHeCDNT225Zm4Gixn2fQVHdqp8+NTHdMMrrdviQ3wxCcg9Wm/voWhnsSa7DIuEGbbaXN0JOV0g33CwnEDzMe458CjszDt5h08fuPETvX+iZjVPdxMxfXx0asVBbTlUgb+X+4maf80ZA3lGLu97PYt301G44zuIOk/lwN08dS5cM0xiqB+HdCB3SFsPHWNUfLBtNld3Uk6Vf5AtDVya/msGiXzEde9D4iTw8lelhTqIawZKClOBXEpp76Fop6hqaGF/SQ2pHS0CdEFOFchr/vcwYyx72TjiSRh4kb2Ho/ViiaF+1DS1Ul6nSxAdyY68Y0gJqYmu3STrx5wqtbIsaC6bW/y4f9od9h6K1sslHa9cOVpHmF456DC2HT6Gm0kwMi7Y3kNxKIbMyIUQM4UQ+4UQB4QQjxlxzLPJEXHsjphNXIivrU6hacDJEkR9wdOxbD10jOSoAIfcYcqerA7kQgg34DVgFjAUuFEIYZNddX9zyWAW/qKTfTY1zQCxfXxwNwldS+5AWswWduZX6vz4WRgxIx8LHJBS5kgpm4FPgKsMOO5ZmfSVaq0HeLiZiAvx5ZCuJXcYGUeqaWgx6/z4WRgRyGOA/FO+Lmi77TRCiHuFEFuFEFvLysoMOK2m2VZiqK/uguhAthw6BqBn5GfRY1UrUsp5UspUKWVqeHh459+gaXaWGObH4XJdgugoth2uICbYh8ggb3sPxeEYEcgLgbhTvo5tu03TnFpSmB/1zWZKa5rsPRSXJ6Vk66FjOq3SDiMC+RZgoBAiSQjhCdwAfG3AcTXNrnQXRMdRcKyB0pomUl18k+X2WB3IpZStwAPAYiAD+I+UMt3a42qavZ1aS67Z19bDFQAdb+vowgwpxpRSfgd8Z8SxNM1RRAf74OlmIleXINrdppwKgnw8GBIZYO+hOCSnWqKvaT3JzSSID/Ult0wHcnvbkFPO2KQQXX7cDh3INa0DAyP8yS6ttfcwXFpRZQOHy+sZ3y/U3kNxWDqQa1oHBkcGcKi8jvrmVnsPxWVtyi0HYHw/nR9vjw7kmtaBIZEBSAnZJXpWbi8bDpYT5ONBcmSgvYfisHQg17QODGkLHvuLa+w8Ete1MaeCcTo/3iEdyDWtA/Ehvvh4uJFRXG3vobikwsoG8ip0frwzOpBrWgdMJsGgvv56Rm4nGw+q/PiE/jqQd0QHck3rxJDIQDKLa3TPFTvYmFNOsK8Hg/vq+vGO6ECuaZ0YHBlARV0zZbW650pP25hbrvPjXaADuaZ14vhqQp1e6VkFx+rJr2hggs6Pd0oHck3rxOC2QJ55RAfynrShLT8+XufHO6UDuaZ1ItTfi/AALzL1jLxHfbmzkJhgHwZF6Px4Z3Qg17QuGBIZwP4SXYLYUw6U1vLDgXJuGhev8+NdoAO5pnXBkMgAskpqaTVb7D0Ul7Bg02E83ATXp8V1/mBNB3JN64rBkYE0t1o4VK43Y7a1+uZW/rutgEtTogjz97L3cJyCDuSa1gW6cqXnfL2ziJrGVm4dn2DvoTgNHcg1rQsGRPjjZhJk6qX6NiWl5P0NhxkSGcAYva1bl+lArmld4O3hRmKor65csbEd+ZXsO1LNrRMSEEJf5OwqHcg1rYuGRAXq1IoNNbaY+evCfQR4uTN7ZIy9h+NUdCDXtC5Kjgwgr6Ke2ia9yYTRpJQ88cVeduRV8tzcEfh5GbKdsMvQgVzTukj3Jredt9bm8vn2Ah6aMZBZKVH2Ho7T0YFc07roxFJ9fcHTUCsyS/jbogwuTYnklxcOtPdwnJIO5JrWRbF9fPD3ctc9Vwy0MrOUn364neSoQJ6/9jy9irObdCDXtC4SQqil+jq1YojF6cXc+8FWBvX158O7xuHrqfPi3aUDuaadgyFRAWQUV+tNJqz03Z4j3L9gO8Oig1hw93j6+Hnae0hOTQdyTTsHQyIDqWlspaiq0d5DcVqL9hzhFx/vYGRcMB/cNZYgHw97D8npWRXIhRB/F0JkCiF2CyG+EEIEGzQuTXNIyVHHe5PrC57tqahrprT67H/olqQX84uPd3BebBDv/mQsAd46iBvB2qTUUuBxKWWrEOJZ4HHgt9YPS9Mc06C+xytXapie3NfOo3E8Pxw4yl3vbaGxxUJMsA8j44OJ6+OLr6cbrRbJ66sOMDwmiPd+MhZ/XStuGKv+J6WUS075ciMw17rhaJpjC/D2ILaPj8ss1W9sMfPhxsNcOTKaiADvDh+7LlsF8aQwP+aOiWVHfiU78ypZml5Cc1v731Hxwbx7p56JG83IP4k/AT5t704hxL3AvQDx8fEGnlbTetaQyECXSa18tq2AJ7/N4O11ubx1expDowPP+ri12WXc/d5WksL8WHD3OEJ/1H621WyhvsVMgJe77qFiA53myIUQy4QQe8/ycdUpj3kCaAUWtHccKeU8KWWqlDI1PDzcmNFrmh0kRwWQc7SOxhazvYdiU1JKPtqUR1KYHxKY+8Z6lqQXn/G4lZml3NUWxD+6Z/wZQRzA3c1EoLeHDuI20mkgl1LOkFIOP8vHVwBCiDuAy4Gbpa7J0lzAkMhAzBbJgdJaew/FpnYVVJFxpJq7Jifx1f2TGNg3gPs+3Mbvv9zDsbpmAL7fe7IW/ON7xhOiywjtwqrUihBiJvAoMFVKqbdO0VzC4FM2mRgeE2Tn0djOR5sO4+vpxlUjownw9uDTe8fzzKJMPth4mG92HeGqkdEs2JTHebFBvHOnLiO0J2vryP8JBABLhRA7hRBvGDAmTXNoiaG+eLmbHLrnitkiabFif9HqxpYTwfr4hUlvDzf+dOUwFj14PikxQby/4TBpiX344K5xOojbmbVVKwOMGoimOQt3NxOD+gY4ROXKN7uKKKxs4OKhfekX7k9Ds5lPt+Tx5tpcvNxNfHrfBMIDOt/3MuNINRtzyrl8RDThAV58taOQhhYzN449szBhUN8APrhrLOlF1QyI8Mfbw80WT007B7qQU9O6YVDfANYdKLPrGBpbzPz2893UN5t5ZlEmAyL8qahrpqKumVHxwWQeqeG2tzfzyb3j250xl1Q38sKS/Xy2rQAp4e+L93P3+f1Ykl7M8JhARsQGn/X7hBC9Oq3kbHQg17RuSAz15fPtTTS2mO02I121v5T6ZjPPX3setY0tLM0oITHUj3un9GNsUgirs8q4+70t3PPeVt6/a+yJcTa2mFmbfZTv9xbz3Z4jmC2SuyYlcfl50cxbc5BXl2cD8NSc4XZ5Xtq504Fc07ohPtQXgLyK+hOrPXvawt1HCPHzZPbIaNzdTNwxKem0+6cOCufF60byy092MPPlNfh5udPYYqaospGGFjOB3u5ccV4UD1ww8MTz+dfNY9iZX8mKjBKuGR1rj6eldYMO5JrWDQmhfgAcLrdPIG9oNrM8o5Q5o2Nwd2u/ZuGK86IxWyT/3VaAl7sJbw83zh8YzvTkCMb3C8XjLN87Mi6YkXHBNhy9ZjQdyDWtGxJC1Az2cHmdXc6/cn8pDS1mLu/CtmizR8Uwe5TezLg3021sNa0bgn09CPR253C5fZZPfLv7CGH+noxNCrHL+TXHogO5pnWDEIKEUD8OV/R8IK9vbmVFZikzh0d2mFbRXIf+KdC0booP9SXPDqmVlZllNLSYuSwlusfPrTkmHcg1rZsSQnwpONZAqxUrKLvj2z1FhPl76bSKdoIO5JrWTQmhvrRaJEd6cNu342mVWcMjcdM7zmttdCDXtG6KD1EliId6ML2yen8ZjS0WZqVE9tg5NcenA7mmdVNi2PESxJ674Pl9ejF9fD0Ym6jTKtpJOpBrWjf1DfDG091EXg9VrjS1mlmRUcrFQ3W1inY6/dOgad1kMgniQ3x7bFHQDweOUtPUykydVtF+RAdyTbNCQohvj6VWvt9bTICXOxP7h/bI+TTnoQO5plkhPtSXvIp6bL3LYavZwtJ9JUxPjsDLXff/1k6nA7mmWSEx1I/6ZjNltU02Pc/m3AqO1bcwc7hOq2hn0oFc06xwop2tjdMri/YW4+PhxtRBETY9j+acdPdDTbPCyS6I9aQaWBJY29TK7Nd+oLiqEZOAumYzFw/ti4+nTqtoZ9KBXNOsENvHF5PA8OZZn2zO40BpLTeOjcPL3Q0pJTeOO3P/TE0DHcg1zSqe7iaignwMbZ7VYrbw9rpcxiaF8LerRxh2XK330jlyTbNSQqivoTPyb3cfoaiqkfum9DPsmFrvpgO5plkpMcyPAyW11Da1Wn0sKSX/XpPDgAh/LhisL2xqXaMDuaZZ6brUOGqaWnl91QGrj7XuwFEyjlRz7/n9MOnuhloX6UCuaVYaGRfMnFExvLk2l3wrUyzz1uQQHuDFVaP0phFa1+lArmkGeHTmYEwCnvk+s9vHyCyuZm32Ue6YmKhXb2rnRAdyTTNAVJAPP53an293H2HLoYpuHeO99Yfw9jBxsy4z1M6RIYFcCPFrIYQUQoQZcTxNc0b3TelPVJA3f/lmHxbLufVeqaxv5osdhcwZFUOwr6eNRqj1VlYHciFEHHAxkGf9cDTNefl4uvGbiwezp7CKJftKzul7P92ST2OLhdsnJtpmcFqvZsSM/CXgUcC27d80zQlcNTKaxFBfXl2e3eWOiK1mC+9vOMz4fiEMiQy08Qi13siqQC6EuAoolFLu6sJj7xVCbBVCbC0rK7PmtJrmsNzdTDxw4UD2HalmWUZpl75nWUYphZUN3DExycaj03qrTgO5EGKZEGLvWT6uAn4H/KErJ5JSzpNSpkopU8PDw60dt6Y5rNkjo4kP8eWV5VldmpW/uz6XmGAfZiTrBUBa93QayKWUM6SUw3/8AeQAScAuIcQhIBbYLoTQDZM1l+buZuKBCwawt7CaFZntz8rNFsk/V2SzMaeCWyck6H04tW7r9k+OlHKPlDJCSpkopUwECoDRUspiw0anaU5qzugY4kJ8eKWdXHleeT3X/XsDzy/J4vIRUdw+IbHnB6n1GnoKoGk24NE2K99dUMXi9NMrWPYWVjHrlTVkldTwyg0j+edNo3Wfcc0qhgXytpn5UaOOp2nO7prRsfQP9+O5xZm0mi2ASqc88cUefDzd+f6hKVw1MsbOo9R6Az0j1zQbcXcz8disZHLK6vh0az4An2zJY1dBFU9cNoSYYB87j1DrLXQg1zQbmpEcQVpiH15elk1+RT3Pfb+fcUkhzNYzcc1AOpBrmg0JIXhsVjJlNU3M+dd66ppaeXL2cITQLWo14+hArmk2NiahDzOHRXK0tom7JicxsG+AvYek9TJ6z05N6wF/vHIo/SP8+Pm0AfYeitYL6UCuaT0gKsiHRy4ZYu9haL2UTq1omqY5OR3INU3TnJwO5JqmaU5OB3JN0zQnpwO5pmmak9OBXNM0zcnpQK5pmubkdCDXNE1zcqKrG8QaelIhyoDD3fz2MMAV2+W64vN2xecMrvm8XfE5w7k/7wQp5Rl7ZdolkFtDCLFVSplq73H0NFd83q74nME1n7crPmcw7nnr1IqmaZqT04Fc0zTNyTljIJ9n7wHYiSs+b1d8zuCaz9sVnzMY9LydLkeuaZqmnc4ZZ+SapmnaKXQg1zRNc3JOFciFEDOFEPuFEAeEEI/Zezy2IISIE0KsFELsE0KkCyEebLs9RAixVAiR3fa5j73HajQhhJsQYocQYmHb10lCiE1tr/enQghPe4/RaEKIYCHEf4UQmUKIDCHEhN7+WgshHm772d4rhPhYCOHdG19rIcTbQohSIcTeU24762srlFfbnv9uIcToczmX0wRyIYQb8BowCxgK3CiEGGrfUdlEK/BrKeVQYDxwf9vzfAxYLqUcCCxv+7q3eRDIOOXrZ4GXpJQDgGPAXXYZlW29AnwvpRwCnId6/r32tRZCxAC/BFKllMMBN+AGeudr/S4w80e3tffazgIGtn3cC7x+LidymkAOjAUOSClzpJTNwCfAVXYek+GklEeklNvb/l2D+sWOQT3X99oe9h4w2y4DtBEhRCxwGfBW29cCuBD4b9tDeuNzDgKmAPMBpJTNUspKevlrjdpi0kcI4Q74Akfoha+1lHINUPGjm9t7ba8C3pfKRiBYCBHV1XM5UyCPAfJP+bqg7bZeSwiRCIwCNgF9pZRH2u4qBvraa1w28jLwKGBp+zoUqJRStrZ93Rtf7ySgDHinLaX0lhDCj178WkspC4HngTxUAK8CttH7X+vj2nttrYpvzhTIXYoQwh/4HHhISll96n1S1Yz2mrpRIcTlQKmUcpu9x9LD3IHRwOtSylFAHT9Ko/TC17oPavaZBEQDfpyZfnAJRr62zhTIC4G4U76Obbut1xFCeKCC+AIp5f/abi45/lar7XOpvcZnA5OAK4UQh1ApswtRuePgtrff0Dtf7wKgQEq5qe3r/6ICe29+rWcAuVLKMillC/A/1Ovf21/r49p7ba2Kb84UyLcAA9uubnuiLpB8becxGa4tNzwfyJBSvnjKXV8Dt7f9+3bgq54em61IKR+XUsZKKRNRr+sKKeXNwEpgbtvDetVzBpBSFgP5QojBbTdNB/bRi19rVEplvBDCt+1n/fhz7tWv9Snae22/Bm5rq14ZD1SdkoLpnJTSaT6AS4Es4CDwhL3HY6PnOBn1dms3sLPt41JUzng5kA0sA0LsPVYbPf9pwMK2f/cDNgMHgM8AL3uPzwbPdySwte31/hLo09tfa+DPQCawF/gA8OqNrzXwMeo6QAvq3ddd7b22gEBV5R0E9qCqerp8Lr1EX9M0zck5U2pF0zRNOwsdyDVN05ycDuSapmlOTgdyTdM0J6cDuaZpmpPTgVzTNM3J6UCuaZrm5P4fVNGdfzg8eFIAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# see what it looks like\n",
    "plt.plot(x[:100,:])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "outputs": [],
   "source": [
    "# get data ready to train with\n",
    "train_dataset = [dict(data=x)]\n",
    "test_dataset = [dict(data=x2)]\n",
    "# First compute the autoregression covariates\n",
    "precompute_ar_covariates(train_dataset, num_lags=num_lags, fit_intercept=True)\n",
    "precompute_ar_covariates(test_dataset, num_lags=num_lags, fit_intercept=True)\n",
    "\n",
    "# Then precompute the sufficient statistics\n",
    "LinearRegressionObservations.precompute_suff_stats(train_dataset)\n",
    "LinearRegressionObservations.precompute_suff_stats(test_dataset)\n",
    "\n",
    "covariates_dim = train_dataset[0]['covariates'].shape[1]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [],
   "source": [
    "# initialize model to train\n",
    "config =  dict(\n",
    "    num_discrete_states=2,\n",
    "    data_dim=data_dim,\n",
    "    covariates_dim=3,\n",
    "    tau_scale=0.7,\n",
    "    num_taus=3,\n",
    "    kappa=100,\n",
    "    alpha=0.5,\n",
    "    covariance_reg=1e-4)\n",
    "twarhmm = TWARHMM(config)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "outputs": [
    {
     "data": {
      "text/plain": "  0%|          | 0/20 [00:00<?, ?it/s]",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "c10ed6b8ab054749961d389e49350cfc"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "  0%|          | 0/1 [00:00<?, ?it/s]",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "2851bf5891ce42239fd6d26415b5e5db"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average iteration time:  0.11177268028259277\n"
     ]
    }
   ],
   "source": [
    "# train\n",
    "train_lls, test_lls, train_posteriors, test_posteriors, = \\\n",
    "        twarhmm.fit_stoch(train_dataset,\n",
    "                         test_dataset,\n",
    "                          num_epochs=20, compute_posteriors=True, wandb_log=False)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0.5, 1.0, 'inferred discrete states')"
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYAklEQVR4nO3dfbRddX3n8feHhGcFFCJKEgja6BQ7OjhZSEenUqUW0CFdnXYKM11qxxFdS0fbOiq2XQxVu6yto1OXjIoPteIUitQ6aSctVsWxMy2W4FMFpEZ8SCJIRHyoluI95zt/7J1wuNwkJ2efcHbOfb/Wuitn7/275/z23Wd/8r2//btnp6qQJB38Dpl1ByRJ02GgS9KcMNAlaU4Y6JI0Jwx0SZoTBrokzQkDfZlLclOSs8Zs+7gkn0nyvSQvPbA9W/L11yWpJCvHbP/eJK9rH//rJLce2B5Ks2WgL3NV9fiq+viYzV8JXFdVD62qtxzAbk1dVf1VVT3uwXzNJF9JcvYBeN5Lk7x/P9qflWT7tPuh/jHQtT9OAW6a5BuXqqrHrbT76GDuu+aXgb7MjVaRbeV3dZL3tcMqNyXZ0G77GPCTwFuT/EOSxyY5PMkbk3wtyTeSvD3JkW37s5JsT/KqJHcAv98+/zVJ3p/ku8Dzkhyb5N1Jbk+yI8nrkqxon2NF+/zfTHIb8Kx97MvpST7V9v2PgCNGtt2vSm37taNte2uSZ4y85q8l+VK77cYka9ttleTFSb4IfLFd9+x2GOrbSf46yRPa9VcAJwN/2v68XtmuP7Nt9+0kn93bcNdSfUxyDvBrwC+0z/vZtu0vJbmlbXtbkhe2648G/hw4qW3/D0lOSnJIkovb/byrPe4Pb7/niPYY3dX284YkJ+7jraQ+qCq/lvEX8BXg7PbxpcA9wHnACuD1wPUjbT8O/KeR5TcDm4CHAw8F/hR4fbvtLGABeANwOHBk+/w/BH6Gppg4EvgT4B3A0cAjgL8FXtg+x4uALwBr29e4Dihg5RL7cRjwVeBXgEOBn2tf63Uj/dnePn4csA04qV1eBzymffwK4O/aNgGeCBzfbivgL9u+HAmcDtwJPLn9eT23/Xkevvhn2y6vBu5qf76HAD/VLq9aYn/21sdLgfcvav8s4DFtn58G/AB40uJ9H2n/MuB6YE17fN4BXNlue2F7LI9q9+tfAsfM+r3q1xjn86w74NeM3wAPDPSPjGw7DfjHkeWP0wZ6Gxzf3xUy7bofB77cPj4LuBc4YmT7pcAnRpZPBP4JOHJk3YU04/QAHwNeNLLtmew50H8C+DqQkXV/zdKB/iNtEJ8NHLroeW4FNu7hZ1XA00eW3wa8donvf9rin227/CrgikXtrwWeu8Rr7a2PDwj0Jb7/Q8DLFu/7yPZbgGeMLD+K5j/AlcB/bH92T5j1+9Ov/ftyyEWL3THy+AfAEXsYL15FU8Hd2P5a/m3gL9r1u+ysqnsWfd+2kcen0FTTt488xztoKnWAkxa1/+pe+n0SsKPadNpb+6raCvwyTTDemeSqJCe1m9cCX9rL6yzu/8t39b3t/9q2L0s5Bfj5Re2fShOm+9PHB0hybpLrk3yrfd7zgBP2sh+nAH8y0o9bgAHNf7JX0PxHc1WSryf5nSSH7uW51BMGuib1TeAfgcdX1XHt17FV9ZCRNkt9lOfoum00FfoJI89xTFU9vt1+O01A7nLyXvpzO7A6ScZpX1V/WFVPpQm2ohka2tWnx+zldRb3/7dG+n5cVR1VVVcu0XZX+ysWtT+6qn57P/t4v+dNcjjwx8AbgROr6jhgM81vUUv1Y1dfzl3UlyOqakdV/bCqfrOqTgP+FfBs4Dl7+ZmoJwx0TaSqhsA7gTcneQRAktVJfno/nuN24MPAf0tyTHuh7jFJntY2uRp4aZI1SR4GXLyXp/sbmjH7lyY5NMnPAmcs1TDNfPqnt0F4D81/TMN287uA1yZZn8YTkhy/h9d8J/CiJE9u2x6d5FlJHtpu/wbw6JH27wf+TZKfbi++HtFerF2zn338BrAuya7z9zCacfCdwEKSc2mGpxhpf3ySY0fWvR34rSSntK+3KsnG9vFPJvnnaS5Of5dmKGaIes9AVxevArYC16eZtfIRmot5++M5NIF0M3A3cA33DUG8k+ZX/88CnwI+uKcnqap7gZ8Fngd8C/iFvbQ/HPhtmt8y7qAZ4nl1u+1NNP+RfJgmzN5NcwF0qdfcArwAeGvb963t6+/yeuA32mGN/1JV24CNNLNUdtJUya9g6fNwb338QPvvXUk+VVXfA17a9vtu4N/TXKze1c8vAFcCt7V9OQn4vbbNh5N8j+YC6ZPbb3kkzXH4Ls1QzP+hGYZRz+X+Q46SpIOVFbokzQkDXZLmhIEuSXPCQJekOTGzDxg64YQTat26dbN6eUk6KN14443frKpVS22bWaCvW7eOLVu2zOrlJemglGSPfzHtkIskzQkDXZLmhIEuSXPCQJekOWGgS9Kc2GegJ3lPkjuTfH4P25PkLUm2JvlckidNv5uSpH0Zp0J/L3DOXrafC6xvvy6iuYuLJOlBts956FX1iSTr9tJkI/C+9k4x1yc5Lsmj2s+6lvbLV775fT746R3gp4Bqjj3jR0/kiWuPm/rzTuMPi1Zz/9tybW/XPSDQk1xEU8Vz8sl7u/mMlqv3X/9V3vV/v8z97jskzZlHHHNEbwN9bFV1OXA5wIYNGyzB9AD3DoY87KhD+fQlz9x3Y0n3M41ZLju4/30f17TrpP22MCxWHGJ5Lk1iGoG+CXhOO9vlTOA7jp9rUoOBgS5Nap9DLkmuBM4CTkiyHfivwKEAVfV2mruLn0dzP8UfAL90oDqr+bcwLFYe4p9HSJMYZ5bLhfvYXsCLp9YjLWvDskKXJmUppF5pKnQDXZqEga5eGQyHVujShAx09cqCF0WliRno6pXBsFi5wkCXJmGgq1eaeei+LaVJeOaoVwZeFJUmZqCrVxa8KCpNzEBXr1ihS5Mz0NUrfpaLNDkDXb1ihS5NzkBXrzTz0H1bSpPwzFGvWKFLkzPQ1SsLwyEr/MMiaSIGunrFCl2anIGuXnGWizQ5A129YoUuTc5AV6/4WS7S5Dxz1CtW6NLkDHT1ysLAz3KRJmWgq1es0KXJGejqlYVhOQ9dmpCBrl6xQpcmZ6CrN6rKWS5SB5456o1hNf9aoUuTMdDVGwvDIYCzXKQJGejqjUFboluhS5Mx0NUbC22gW6FLkzHQ1RuDgRW61IWBrt7YXaGv8G0pTcIzR73hGLrUjYGu3nCWi9TNWIGe5JwktybZmuTiJbafnOS6JJ9O8rkk502/q5p3VuhSN/sM9CQrgMuAc4HTgAuTnLao2W8AV1fV6cAFwP+Ydkc1/5zlInUzToV+BrC1qm6rqnuBq4CNi9oUcEz7+Fjg69PropaL+yp0RwKlSYxz5qwGto0sb2/XjboU+MUk24HNwH9e6omSXJRkS5ItO3funKC7mmcLAyt0qYtplUIXAu+tqjXAecAVSR7w3FV1eVVtqKoNq1atmtJLa144hi51M06g7wDWjiyvadeNej5wNUBV/Q1wBHDCNDqo5WP3LBc/D12ayDiBfgOwPsmpSQ6juei5aVGbrwHPAEjyozSB7piK9osVutTNPgO9qhaAlwDXArfQzGa5KclrkpzfNns58IIknwWuBJ5XVXWgOq355CwXqZuV4zSqqs00FztH110y8vhm4CnT7ZqWG2e5SN145qg3rNClbgx09cagvSjqGLo0GQNdveE8dKkbA129sXsM3WmL0kQMdPXGgtMWpU4MdPXGYPdFUd+W0iQ8c9QbVuhSNwa6emPgDS6kTgx09YYVutSNga7e2DWGfoiBLk3EQFdv7JqHboUuTcZAV28Myz8skrow0NUbC344l9SJZ456Y+CHc0mdGOjqDcfQpW4MdPXGYDgkcZaLNCkDXb2xMCyrc6kDA129MRiW4+dSBwa6eqOp0H1LSpPy7FFvWKFL3Rjo6o2F4dAxdKkDA129YYUudWOgqzcWBs5ykbow0NUbg2GxwvuJShMz0NUbznKRuvHsUW84hi51Y6CrN5zlInVjoKs3rNClbgx09Yaf5SJ1Y6CrN6zQpW7GCvQk5yS5NcnWJBfvoc2/S3JzkpuS/OF0u6nloJmHbo0hTWrlvhokWQFcBvwUsB24Icmmqrp5pM164NXAU6rq7iSPOFAd1vyyQpe6GaccOgPYWlW3VdW9wFXAxkVtXgBcVlV3A1TVndPtppaDheGQlf5hkTSxcQJ9NbBtZHl7u27UY4HHJvl/Sa5Pcs5ST5TkoiRbkmzZuXPnZD3W3LJCl7qZ1oDlSmA9cBZwIfDOJMctblRVl1fVhqrasGrVqim9tOaFs1ykbsYJ9B3A2pHlNe26UduBTVX1w6r6MvD3NAEvjc0KXepmnEC/AVif5NQkhwEXAJsWtfkQTXVOkhNohmBum143tRz4WS5SN/s8e6pqAXgJcC1wC3B1Vd2U5DVJzm+bXQvcleRm4DrgFVV114HqtOaTFbrUzT6nLQJU1WZg86J1l4w8LuBX2y9pIn6Wi9SNv9+qNwYDK3SpCwNdvbEwLOehSx0Y6OoNx9Clbgx09YazXKRuPHvUG1boUjcGunrDWS5SNwa6esMKXerGQFdv+FkuUjcGunphOCyqYIUXRaWJefaoFxaGBeA8dKkDA129MGgD3TF0aXIGunphYTgEcAxd6sBAVy9YoUvdGejqhd1j6Aa6NDEDXb1wX4XuW1KalGePesEKXerOQFcvDAaOoUtdGejqhd2zXJyHLk3MQFcvOMtF6s5AVy84hi51Z6CrF5zlInXn2aNesEKXujPQ1QuD9qLoIQa6NDEDXb2wMLBCl7oy0NULznKRujPQ1QuDskKXujLQ1QsLVuhSZwa6emGwewzdt6Q0Kc8e9YIVutSdga5eGHhPUakzA129sOvDuazQpcmNFehJzklya5KtSS7eS7t/m6SSbJheF7UcDPxLUamzfQZ6khXAZcC5wGnAhUlOW6LdQ4GXAZ+cdic1/xxDl7obp0I/A9haVbdV1b3AVcDGJdq9FngDcM8U+6dl4r4K3VFAaVLjnD2rgW0jy9vbdbsleRKwtqr+996eKMlFSbYk2bJz58797qzmlxW61F3ncijJIcCbgJfvq21VXV5VG6pqw6pVq7q+tObIYNDeschAlyY2TqDvANaOLK9p1+3yUODHgI8n+QpwJrDJC6PaH7srdKctShMbJ9BvANYnOTXJYcAFwKZdG6vqO1V1QlWtq6p1wPXA+VW15YD0WHPJWS5Sd/sM9KpaAF4CXAvcAlxdVTcleU2S8w90B7U8OIYudbdynEZVtRnYvGjdJXtoe1b3bmm5cZaL1J1nj3phV4VugS5NzkBXLwyGQ1YeEhITXZqUga5eWBiW4+dSRwa6emEwKGe4SB0Z6OoFK3SpOwNdvTAYFitX+HaUuvAMUi9YoUvdGejqhV2zXCRNzkBXL1ihS90Z6OqFwdBZLlJXBrp6wQpd6s5AVy8089B9O0pdeAapF6zQpe4MdPXCYDhkpTe3kDox0NULVuhSdwa6esFZLlJ3Brp6wQpd6s5AVy80FbpvR6kLzyD1ghW61J2Brl7ws1yk7gx09cLCwApd6spAVy80n4duoEtdGOjqhcGwWOFFUakTzyD1woLz0KXODHT1wsBZLlJnBrp6YcFZLlJnBrp6wQpd6s5AVy84hi51Z6CrFwYDZ7lIXXkGqRcWnIcudTZWoCc5J8mtSbYmuXiJ7b+a5OYkn0vy0SSnTL+rmmeOoUvd7TPQk6wALgPOBU4DLkxy2qJmnwY2VNUTgGuA35l2RzXfnOUidTdOhX4GsLWqbquqe4GrgI2jDarquqr6Qbt4PbBmut3UPBsOi2FhhS51NE6grwa2jSxvb9ftyfOBP19qQ5KLkmxJsmXnzp3j91JzbVAFYIUudTTVi6JJfhHYAPzuUtur6vKq2lBVG1atWjXNl9ZBbDBsAt1ZLlI3K8doswNYO7K8pl13P0nOBn4deFpV/dN0uqflYGF3oM+4I9JBbpxT6AZgfZJTkxwGXABsGm2Q5HTgHcD5VXXn9LupeTYYWKFL07DPM6iqFoCXANcCtwBXV9VNSV6T5Py22e8CDwE+kOQzSTbt4emkB3AMXZqOcYZcqKrNwOZF6y4ZeXz2lPulZWRhOASc5SJ15e+4mrldF0Wt0KVuDHTN3MLuMXQDXerCQNfM7a7Q/SwXqRMDXTO34Dx0aSo8gzRzjqFL02Gga+ac5SJNh4GumbNCl6bDQNfM3TeGbqBLXRjomrn7KnTfjlIXnkGaOeehS9NhoGvmnIcuTYeBrplzlos0HQa6Zs5ZLtJ0GOiaOWe5SNNhoGvmnOUiTYdnkGbOCl2aDgNdMzdoL4o6hi51Y6Br5pyHLk2Hga6Zcx66NB0GumbOMXRpOgx0zZyzXKTp8AzSzFmhS9NhoGvmnOUiTYeBrpmzQpemw0DXzA0GfpaLNA0GumbOCl2aDgNdMzcYFisOCYmBLnVhoGvmFtpAl9SNga6ZGwyHjp9LU2Cga+as0KXpMNA1c4NhWaFLUzBWoCc5J8mtSbYmuXiJ7Ycn+aN2+yeTrJt6TzW3mgrd2kLqap9nUZIVwGXAucBpwIVJTlvU7PnA3VX1I8CbgTdMu6OaX4OBFbo0DSvHaHMGsLWqbgNIchWwEbh5pM1G4NL28TXAW5OkqmqKfQXg6hu28c6/um3aT6sZuuO793DMEYfOuhvSQW+cQF8NbBtZ3g48eU9tqmohyXeA44FvjjZKchFwEcDJJ588UYePO+pQ1p/4kIm+V/20/sSHcOajj591N6SD3jiBPjVVdTlwOcCGDRsmqt6f+fhH8szHP3Kq/ZKkeTDOlagdwNqR5TXtuiXbJFkJHAvcNY0OSpLGM06g3wCsT3JqksOAC4BNi9psAp7bPv454GMHYvxckrRn+xxyacfEXwJcC6wA3lNVNyV5DbClqjYB7wauSLIV+BZN6EuSHkRjjaFX1WZg86J1l4w8vgf4+el2TZK0P/xrDkmaEwa6JM0JA12S5oSBLklzIrOaXZhkJ/DVCb/9BBb9FeoysRz3eznuMyzP/V6O+wz7v9+nVNWqpTbMLNC7SLKlqjbMuh8PtuW438txn2F57vdy3GeY7n475CJJc8JAl6Q5cbAG+uWz7sCMLMf9Xo77DMtzv5fjPsMU9/ugHEOXJD3QwVqhS5IWMdAlaU4cdIG+rxtWz4Mka5Ncl+TmJDcleVm7/uFJ/jLJF9t/Hzbrvk5bkhVJPp3kz9rlU9sbj29tb0R+2Kz7OG1JjktyTZIvJLklyY8vk2P9K+37+/NJrkxyxLwd7yTvSXJnks+PrFvy2KbxlnbfP5fkSfv7egdVoI95w+p5sAC8vKpOA84EXtzu58XAR6tqPfDRdnnevAy4ZWT5DcCb2xuQ301zQ/J583vAX1TVPwOeSLP/c32sk6wGXgpsqKofo/lo7guYv+P9XuCcRev2dGzPBda3XxcBb9vfFzuoAp2RG1ZX1b3ArhtWz5Wqur2qPtU+/h7NCb6aZl//oG32B8DPzKSDB0iSNcCzgHe1ywGeTnPjcZjPfT4W+AmaewpQVfdW1beZ82PdWgkc2d7l7CjgdubseFfVJ2juETFqT8d2I/C+alwPHJfkUfvzegdboC91w+rVM+rLgyLJOuB04JPAiVV1e7vpDuDEWfXrAPnvwCuBYbt8PPDtqlpol+fxeJ8K7AR+vx1qeleSo5nzY11VO4A3Al+jCfLvADcy/8cb9nxsO+fbwRboy0qShwB/DPxyVX13dFt7i7+5mXOa5NnAnVV146z78iBbCTwJeFtVnQ58n0XDK/N2rAHaceONNP+hnQQczQOHJubetI/twRbo49ywei4kOZQmzP9nVX2wXf2NXb+Ctf/eOav+HQBPAc5P8hWaobSn04wtH9f+Sg7zeby3A9ur6pPt8jU0AT/PxxrgbODLVbWzqn4IfJDmPTDvxxv2fGw759vBFujj3LD6oNeOHb8buKWq3jSyafRm3M8F/teD3bcDpapeXVVrqmodzXH9WFX9B+A6mhuPw5ztM0BV3QFsS/K4dtUzgJuZ42Pd+hpwZpKj2vf7rv2e6+Pd2tOx3QQ8p53tcibwnZGhmfFU1UH1BZwH/D3wJeDXZ92fA7SPT6X5NexzwGfar/NoxpQ/CnwR+Ajw8Fn39QDt/1nAn7WPHw38LbAV+ABw+Kz7dwD2918AW9rj/SHgYcvhWAO/CXwB+DxwBXD4vB1v4EqaawQ/pPlt7Pl7OrZAaGbxfQn4O5oZQPv1ev7pvyTNiYNtyEWStAcGuiTNCQNdkuaEgS5Jc8JAl6Q5YaBL0pww0CVpTvx/9esO8n78MMYAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get inferred discrete states via floor division by num_taus\n",
    "inferred_states = test_posteriors[0].get_states()//config['num_taus']\n",
    "plt.plot(inferred_states[:100])\n",
    "plt.title('inferred discrete states')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0.5, 1.0, 'inferred time constants')"
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAibUlEQVR4nO3debhcdZ3n8ffnbsUqCSTakIVEjQsuLHMFfOARtBGDbYtjtxraUWS0M/pIu4wzNowzwmDPqGOPOrYopDWgjg0ioh2ZKKKCuIG5IKJsEiKSxGAuhFXgrt/545y6t3Kpm1upqnNO5dTn9Tz1pOps9Tt1Ut/63e/5nfNVRGBmZuXVU3QDzMwsWw70ZmYl50BvZlZyDvRmZiXnQG9mVnIO9GZmJedAb7OSdKukExtc9rmSbpb0qKT3ZNuyuu+/TFJI6mtw+aWSHpPUm3XbzIrmQG+ziogXRMS1DS7+QeCaiNg/Ij6TYbOaIukeSSdVX0fEvRGxX0RMFNmuLKQ/eM9u07Z26wfUOpMDvbXLocCtzaxYL4g4sJi1jwO9zaq2FyzpXEmXSfpymp65VdJgOu+HwMuBz6bpkOdIqkj6R0n3SvqjpAsk7Z0uf6KkLZL+XtJ9wEXp9i+X9H8lPQK8TdIBkr4oaZukrZL+oZpqkdSbbv9+SZuAv9jFfnwFWAp8O23fB2f2VCVdm27/Z+ky35Z0kKSvSnpE0gZJy2q2+TxJV0vaIelOSW/cxfsfKOkiSX+Q9KCkb9XM+1tJG9PtrJN0SM28kPROSXdJekjS+ZKUznu2pB9Jejj9DL6WTr8uXf1X6X68SdJ8SVdKGk7f/0pJi2ve51pJH5H00/TYfk/SgnR2dXsPpdt76WzvbR0sIvzwo+4DuAc4KX1+LvAk8GqgF/gocH3NstcC76h5/SlgHXAgsD/wbeCj6bwTgXHg40AF2Dvd/hjwOpIOyN7AN4ELgX2BpwO/AP5Duo13AncAS9L3uAYIoG+ufUlfL6tdPm3/RuBZwAHAbcBvgZOAPuDLwEXpsvsCm4Ez0nlHAvcDh83y3v8P+BowH+gHTkinvyJd76j0c/gn4Lqa9QK4EphH8kM1DKxM510CfCj9rPYCjp+x3rNrXh8E/BWwT3osvg58a8axuxt4Tvq5Xwt8rN7nNNd7+9GZD/fobXf8JCLWR5LX/gpweL2F0l7nauD9EbEjIh4F/iewqmaxSeCciBiJiCfSaT+PiG9FxCTwNJIflfdFxJ8iYjvJj0d1G28EPh0RmyNiB8kPT6suioi7I+Jh4DvA3RHx/YgYJwmOR6bLvQa4JyIuiojxiPgl8A3gDXU+i4OBU4B3RsSDETEWET9KZ78ZWBsRN0XECHA28NLavxxIAu5DEXEvyY/ZEen0MZJ02SER8WRE/GS2nYqIByLiGxHxeHos/gdwQp19/216LC6reZ96Gn5v6wwO9LY77qt5/jiw1yy59IUkvccb05TDQ8B30+lVwxHx5Iz1Ntc8P5Sk97utZhsXkvTsAQ6Zsfzvd3Nf6vljzfMn6rzer6Ztx1TblbbtzcCf1dnmEmBHRDxYZ94h1LQ7Ih4DHgAW1Swz8zOvtuGDgIBfpGm0fz/bTknaR9KFkn6fpsWuA+Zp5xFHs71PPQ2/t3UGn/CyLNxPEhhfEBFbZ1mm3m1Ta6dtBkaABWmPeqZtJEG0aukcbWrnbVo3Az+KiFc2uOyBkuZFxEMz5v2B5EcDAEn7kqRZZvvMpkTEfcDfpusdD3xf0nURsbHO4h8AngscExH3SToC+CVJsJ7zrVp8b+sA7tFb26Wpl38GPiXp6QCSFkl61W5sYxvwPeB/S3qapB5Jz5JUTTlcBrxH0mJJ84Gz5tjkH4Fn7vbO1Hcl8BxJb5HUnz5eIun5s+zHd4DPpSdF+yW9LJ19CXCGpCMkVUjSWzdExD1zNUDSG2pOqD5IEpAn09cz93V/kh/ehyQdCJyzG/s6nG53antzvLd1IAd6y8rfk5zcvD5NF3yfpFe5O94KDJCcGH0QuBw4OJ33z8BVwK+Am4Ar5tjWR4H/mqZa/tNutmMnaZ77ZJLzBX8gSXtUTyzX8xaSvPYdwHbgfel2vg/8N5L8/jaSE8Gr6m/iKV4C3CDpMZKT3u+NiE3pvHOBL6X7+kbg0yQnWe8HridJozUkIh4nyen/NN3esXO8t3UgRbjwiJlZmblHb2ZWcg70ZmYl50BvZlZyDvRmZiXXkePoFyxYEMuWLSu6GWZme4wbb7zx/ohYWG9eRwb6ZcuWMTQ0VHQzzMz2GJJmvTrcqRszs5JzoDczKzkHejOzknOgNzMrOQd6M7OSmzPQS1oi6RpJt6X3nn5vnWUk6TNpSbRbJB1VM+/0tBTaXZJOb/cOmJnZrjUyvHIc+EBE3CRpf5JiEldHxG01y5wCrEgfxwCfJynMUL0l6iDJrUxvlLRuliIMZmaWgTkDfXo/7W3p80cl3U5SAac20J8KfDmSW2FeL2leWkLtRODqtNQbkq4GVpLch7srXPqLe/nDQ0/MvaCVSk+PeNNLlnDwAXsX3RSz3btgKq1leSRww4xZi9i5rNuWdNps0+ttezVJnVGWLp2rWNCe4eHHxzjril8DoEZq+VhpRECvxN/9+Yqim2LWeKCXtB9JgYT3RcQj7W5IRKwB1gAMDg6W4ib5T4xNAPDR17+I044ux4+XNWbFh9ZPHX+zojU06kZSP0mQ/2pE1Kvks5Wd63cuTqfNNr0rjIwnX/RKnwc3dZuB3h5Gxl1dzzpDI6NuBHwRuD0iPjnLYuuAt6ajb44FHk5z+1cBJ6e1MueTlF+7qk1t73jVL/qAA33XqfT3Tv3QmxWtkdTNcSQ1L38t6eZ02n8BlgJExAXAeuDVJDVCHwfOSOftkPQRYEO63nnVE7PdYGQsCfSVvt6CW2J5q/T1TB1/s6I1MurmJ8AuTyWmo23ePcu8tcDaplq3h3PqpntV+py6sc7hCJSh6hfdgb77VPqcurHO4QiUoakefb9TN92m0t/DqHv01iEc6DM0naP3x9xtnLqxTuIIlKHRCQf6bpWkbhzorTM4AmVoqkfv1E3XSXr0ztFbZ3Cgz5BH3XSvSr+HV1rncATKkEfddC+nbqyTOAJlaDrQO3XTbZy6sU7iQJ+hkbEJJOjv9a0ru41H3VgncaDP0Mj4JJW+HuR7FHedSn+vc/TWMRzoM5QEeqdtulE1dZPcHcSsWA70GRoZn/CJ2C5V6ethMmB80oHeiucolKGRsUkq/f6Iu1H1Lznn6a0TOAplyKmb7lX9gR9xlSnrAA70GXLqpntVj7t79NYJHIUyVB11Y93HqRvrJHMWHpG0FngNsD0iXlhn/n8G3lyzvecDC9PqUvcAjwITwHhEDLar4XuCkTGnbrrVdI/eqRsrXiPdzYuBlbPNjIhPRMQREXEEcDbwoxnlAl+ezu+qIA9p6sYnY7vSdI7ePXor3pxRKCKuAxqt83oacElLLSoRp266l1M31knaFoUk7UPS8/9GzeQAvifpRkmr51h/taQhSUPDw8PtalahPOqmezl1Y52knd3NvwR+OiNtc3xEHAWcArxb0stmWzki1kTEYEQMLly4sI3NKs7ImEfddKupHr1TN9YB2hmFVjEjbRMRW9N/twPfBI5u4/t1vJFxXzDVrQY8vNI6SFuikKQDgBOAf62Ztq+k/avPgZOB37Tj/fYUTt10L6durJM0MrzyEuBEYIGkLcA5QD9ARFyQLvZvge9FxJ9qVn0G8M30zo19wL9ExHfb1/TO5wumutfUqBv36K0DzBnoI+K0Bpa5mGQYZu20TcDhzTZsTzcxGYxNhHv0XWo6R+8evRXP3c2MjFarSzlH35V8CwTrJI5CGXFh8O5WPe6jDvTWARyFMuJ6sd2tr7eH3h65R28dwYE+I9Xx0+7Rdy8XCLdO4SiUkanUjXP0XcsFwq1TOAplxKkbq/S5QLh1Bgf6jFR79ANO3XStSr9TN9YZHIUy4hy9OXVjncJRKCPTqRt/xN2q0tfrQG8dwVEoI9Pj6J2j71YedWOdwoE+IyO+MrbrVfp7fDLWOoKjUEacozenbqxTOAplZGTCwyu7nVM31ikc6DNSvWuhUzfdy6NurFM4CmXEo27MF0xZp3AUykg10A/0+iPuVr5gyjrFnFFI0lpJ2yXVLQMo6URJD0u6OX18uGbeSkl3Stoo6ax2NrzTVatLpRW2rAs5dWOdopHu5sXAyjmW+XFEHJE+zgOQ1AucD5wCHAacJumwVhq7JxkZm3Tapst51I11ijkjUURcB+xoYttHAxsjYlNEjAKXAqc2sZ090sj4JJV+j7jpZpW+HiYmg/EJB3srVru6nC+V9CtJ35H0gnTaImBzzTJb0ml1SVotaUjS0PDwcJuaVRwXBjcXCLdO0Y5IdBNwaEQcDvwT8K1mNhIRayJiMCIGFy5c2IZmFWtk3KmbbjdVINyB3grWciSKiEci4rH0+XqgX9ICYCuwpGbRxem0rpDk6J266WbTBcI98saK1XKgl/RnSoeWSDo63eYDwAZghaTlkgaAVcC6Vt9vTzEyPuGLpbpctRaBx9Jb0frmWkDSJcCJwAJJW4BzgH6AiLgA+GvgXZLGgSeAVRERwLikM4GrgF5gbUTcmsledCCnbsypG+sUcwb6iDhtjvmfBT47y7z1wPrmmrZnGxmfZN7e/UU3wwrk1I11Cnc5MzIy5lE33c6jbqxTOBJlZNTj6LveVOrGOXormAN9Rpyjt+rxH51w6saK5UiUEV8wZVOpG/forWCORBnxOHrzqBvrFA70GUnudeOPt5t51I11CkeiDExOBqMTztF3u+lA7x69FcuRKAOjrhdrMDXqyjl6K5oDfQaqX2z36LubUzfWKRyJMlD9YjtH3936ekSPnLqx4jkSZWC6MLhTN91MkqtMWUdwoM/AVI/eqZuuV+nvYWTMqRsrliNRBp50jt5SLhBuncCRKANTqRvf66brOXVjncCBPgNO3VhV0qN36saKNWckkrRW0nZJv5ll/psl3SLp15J+Junwmnn3pNNvljTUzoZ3smoPbsCBvuslOXr36K1YjUSii4GVu5j/O+CEiHgR8BFgzYz5L4+IIyJisLkm7nk8jt6qnLqxTtBIhanrJC3bxfyf1by8nqQIeFebTt04R9/tnLqxTtDuLufbge/UvA7ge5JulLR6VytKWi1pSNLQ8PBwm5uVr+lx9O7RdzuPurFOMGePvlGSXk4S6I+vmXx8RGyV9HTgakl3RMR19daPiDWkaZ/BwcFoV7uKMD3qxoG+21X6ep2jt8K1JRJJejHwBeDUiHigOj0itqb/bge+CRzdjvfrdNULZJy6sUq/UzdWvJYDvaSlwBXAWyLitzXT95W0f/U5cDJQd+RO2UzfvdI9+m7n1I11gjlTN5IuAU4EFkjaApwD9ANExAXAh4GDgM9JAhhPR9g8A/hmOq0P+JeI+G4G+9BxPOrGqgYc6K0DNDLq5rQ55r8DeEed6ZuAw5+6RvmNjE8y0NdD+iNnXSzJ0Tt1Y8VylzMDLgxuVU7dWCdwNMrAyLgLg1ui0tfL+GQwPuFgb8VxoM/AyJjrxVqiOsR21IHeCuRolIGR8QmPoTdg+oT8qNM3ViBHoww4dWNV1f8HztNbkRzoM5AEen+0VlMg3FfHWoEcjTIwMuZRN5aopvB8dawVydEoAyPjk64uZYBTN9YZHOgz4NSNVU2lbtyjtwI5GmXAF0xZlXP01gkcjTKQjKN36samC8Q7dWNFcqDPQJKj90drTt1YZ3A0yoBTN1Y1Hejdo7fiOBplwBdMWdVU6sY5eiuQA32bRQSjHnVjKadurBM4GrWZ68VaLadurBM0FI0krZW0XVLdUoBKfEbSRkm3SDqqZt7pku5KH6e3q+GdairQO3Vj+IIp6wyNdjsvBlbuYv4pwIr0sRr4PICkA0lKDx5DUhj8HEnzm23snqD6J7pTNwbQ3yskXGXKCjVnKUGAiLhO0rJdLHIq8OWICOB6SfMkHUxSa/bqiNgBIOlqkh+MS1pq9Szee+kvOeE5C3n9UYuz2PysHh8d59x1t/Lok+M8MeZAb9MkUenr4cpbtnHX9sea2saLF8/jXSc+q80ts27SUKBvwCJgc83rLem02aY/haTVJH8NsHTp0qYa8YPbt7NgvwqvP2ruZdvp9m2PcNnQFhbN25t9K728cNHTOHzJvHwbYR3rL198CL/a8hB3D+9+oB9+dIQf33W/A721pF2BvmURsQZYAzA4OBjNbCOpz5n/n8jVoXOffOPhHPPMg3J/f+tsn3jD4U2v+/Hv3sEXf/y7NrbGulG78gtbgSU1rxen02abnolKX08h45WnR9r4BKy1V6Wvh9GJSSYnm+r7mAHtC/TrgLemo2+OBR6OiG3AVcDJkuanJ2FPTqdlYqCvp5DRDT4Ba1mpjtpxzVlrRUOpG0mXkJxYXSBpC8lImn6AiLgAWA+8GtgIPA6ckc7bIekjwIZ0U+dVT8xmodLXW0zqZmpIpQO9tVft3S/38l+M1qRGR92cNsf8AN49y7y1wNrdb9ruq/QX1KMfc+rGsrFzhar+Yhtje6xSdUGLy9E7dWPZ8AVX1g6likxO3VjZ+F451g6likyVwk7G+rYHlo1qoH/Sd7+0FpQr0BeWo59ASi53N2snV6iydihXoC8wdVPp60FyoLf2curG2qFkgb64C6actrEs+DbH1g7lC/QFXTDlE7GWhalRN87RWwtKFZ0q/b2MFjSOfsCB3jIw4NSNtUGpolP1pmbJ9Vv5GXHpQMuIUzfWDqWKTpW+HiYDxnO+AVSSunGO3tpv+spYB3prXskCfTFD0UbGJ10j1jIxnaN36saaV6roNNX7yflLMTLm1I1lw6kba4dSRaeivhRO3VhWHOitHUoW6AtM3bhHbxmQxEBfTyGjyaw8ShWdirqKcHR80rcotswUVSLTyqNcgX4qR+8evZVHcmsP9+iteQ1FJ0krJd0paaOks+rM/5Skm9PHbyU9VDNvombeuja2/SmKS934yljLTlG39rDymLPClKRe4HzglcAWYIOkdRFxW3WZiHh/zfJ/BxxZs4knIuKItrV4F4pK3SSjbpy6sWwkd2V16saa10g39GhgY0RsiohR4FLg1F0sfxpwSTsat7uKui+Ix9Fblpy6sVY1Ep0WAZtrXm9Jpz2FpEOB5cAPaybvJWlI0vWSXjfbm0hanS43NDw83ECznqqIqwgnJ4PRCefoLTtF3azPyqPd0WkVcHlE1P6deWhEDAJ/A3xa0rPqrRgRayJiMCIGFy5c2NSbF5G6GZ1wdSnLVpKjd+rGmtdIoN8KLKl5vTidVs8qZqRtImJr+u8m4Fp2zt+3VREnY6tpIvfoLSuVfqdurDWNRKcNwApJyyUNkATzp4yekfQ8YD7w85pp8yVV0ucLgOOA22au2y5TPfocez/Vvx6co7esOHVjrZpz1E1EjEs6E7gK6AXWRsStks4DhiKiGvRXAZfGzvcIfj5woaRJkh+Vj9WO1mm3InL0LgxuWfMFU9aqOQM9QESsB9bPmPbhGa/PrbPez4AXtdC+3TLQW0SgT3v0Tt1YRip9vR5Hby0pVXTq6+2ht0e59n6edI7eMpaMo3egt+aVLjrlfRXhVOrG97qxjDh1Y60qZ6B36sZKxBdMWatKF52SL0Weo26curFsVdLbFOddC9nKo3TRKe985vQ4eqduLBuuG2utKl+gzz1H73H0lq2i7spq5VG66OTUjZVNUXdltfIoXXTK/2SsUzeWrekrvt2jt+aUL9D351tfs3q7BaduLCvVobtO3VizShed8h6K5tSNZc2pG2tV6aJT3heXVAN99fYLZu02Hejdo7fmlC46FXHBVKWvB0m5vad1l6Iqp1l5lDDQ53sDqKRebOk+Rusg0+Ponbqx5pQuQuVdSDmpF+sRN5Ydp26sVeUL9AWlbsyy4gumrFUNRShJKyXdKWmjpLPqzH+bpGFJN6ePd9TMO13SXenj9HY2vp4iRt040FuWiqicZuUyZ+ERSb3A+cArgS3ABknr6lSK+lpEnDlj3QOBc4BBIIAb03UfbEvr66j09TAxGYxPTNKXw0iYkbFJBnyxlGXIqRtrVSOR8GhgY0RsiohR4FLg1Aa3/yrg6ojYkQb3q4GVzTW1MXnfAMqpG8uaUzfWqkYi1CJgc83rLem0mf5K0i2SLpe0ZDfXRdJqSUOShoaHhxtoVn15fymcurGsedSNtapdEerbwLKIeDFJr/1Lu7uBiFgTEYMRMbhw4cKmG5L3VYQedWNZm6qF7HH01qRGAv1WYEnN68XptCkR8UBEjKQvvwD8m0bXbbep3k9OX4qRMaduLFs9PWKg13VjrXmNRKgNwApJyyUNAKuAdbULSDq45uVrgdvT51cBJ0uaL2k+cHI6LTN5p25GnbqxHLhurLVizlE3ETEu6UySAN0LrI2IWyWdBwxFxDrgPZJeC4wDO4C3pevukPQRkh8LgPMiYkcG+zGlkNSNR91YxvK+K6uVy5yBHiAi1gPrZ0z7cM3zs4GzZ1l3LbC2hTbulkJOxvoWxZYxFwi3VpQuQuWeo/fwSstB3ld8W7mULkI5dWNlNNDX4ytjrWmlC/QDOV5FGBE+GWu5qPQ7dWPNK12Ems7RZ9/7maou5Ry9ZcyjbqwVpYtQeRZSdmFwy4tz9NaK8gb6HL4U1R6WUzeWtbwL6li5lC5CVW9HkEvqZsyFwS0feRfUsXIpXYQqJHXje91Yxpy6sVaULtD39YgeOXVj5eILpqwVpYtQkqj09TI6kefJ2NJ9jNZhKh5Hby0oZYSq9OfzpZjO0Tt1Y9lKcvTu0Vtzyhnoc8pnTqVuPI7eMlZN3URE0U2xPVApI1Re+Uynbiwv1f9jeaQkrXxKGaHyuorQF0xZXlwg3FpRzkDf35PP8Moxj7qxfExdH+KLpqwJpYxQuadunKO3jOV9V1Yrl4YilKSVku6UtFHSWXXm/0dJt0m6RdIPJB1aM29C0s3pY93MdbPg1I2VjVM31oo5K0xJ6gXOB14JbAE2SFoXEbfVLPZLYDAiHpf0LuB/AW9K5z0REUe0t9m7Vunr4bGR8czfxxdMWV6m7srq1I01oZEIdTSwMSI2RcQocClwau0CEXFNRDyevrweWNzeZu6evG4A5XvdWF6mKqc5dWNNaCRCLQI217zekk6bzduB79S83kvSkKTrJb1utpUkrU6XGxoeHm6gWbPL6wZQI+OTDPT1ICnz97Lu5tSNtaKh4uCNkvTvgEHghJrJh0bEVknPBH4o6dcRcffMdSNiDbAGYHBwsKWrQvK8YMq9ectD3kXvrVwaiVJbgSU1rxen03Yi6STgQ8BrI2KkOj0itqb/bgKuBY5sob0NyXPUjU/EWh6m78rq1I3tvkYC/QZghaTlkgaAVcBOo2ckHQlcSBLkt9dMny+pkj5fABwH1J7EzUReN4AaGXO9WMvHXv1O3Vjz5kzdRMS4pDOBq4BeYG1E3CrpPGAoItYBnwD2A76e5qvvjYjXAs8HLpQ0SfKj8rEZo3UykdcNoEbGJzyG3nLh1I21oqEcfUSsB9bPmPbhmucnzbLez4AXtdLAZlT6ehmfDMYnJunrzS4QO3VjefEFU9aKUnZH87oBVBLoS/kRWofxOHprRSmjVF7lBEfGJhhwoLccDHh4pbWglFFqukC4e/RWDgNO3VgLShml8spnOkdveentEf29co/emlLKQJ/Xn7kedWN5yuvWHlY+pYxSeZ248jh6y1Ned2W18illlHLqxsoor1t7WPmUPNBn+6UY9b1uLEeV/l5GHeitCaWMUtOjbnLo0TtHbzlx6saaVcoolcc4+ohw6sZy5dSNNavUgT7LK2Or23bqxvLiUTfWrFJGqanUTYZfiul6saX8CK0D5VVQx8qnlFEqj1E3U2UE+526sXw4dWPNKnmgz7JH78Lglq+8CupY+ZQySuVx726nbixvHnVjzSpllOrvFVK2ZdemUjcedWM5qfT3+GSsNaWhQC9ppaQ7JW2UdFad+RVJX0vn3yBpWc28s9Ppd0p6VRvbvqv2Zp7PnErdeBy95cSpG2vWnFFKUi9wPnAKcBhwmqTDZiz2duDBiHg28Cng4+m6h5HUmH0BsBL4XLq9zGX9pXDqxvLm1I01q5FSgkcDGyNiE4CkS4FT2bnI96nAuenzy4HPKikeeypwaUSMAL+TtDHd3s/b0/zZVfp6uOKmLfx04/2ZbP/xUZ+MtXxV+np4cmySV37yR0U3xTIyf58BLnvnS9u+3UYC/SJgc83rLcAxsy2TFhN/GDgonX79jHUX1XsTSauB1QBLly5tpO279K4Tn8WGe3a0vJ1dOf7ZCzjs4AMyfQ+zqpUvPJhN9/+JyYiim2IZedpe/Zlst6Hi4HmIiDXAGoDBwcGW/yefcdxyzjhuecvtMusUhx3yND77N0cV3QzbAzWSd9gKLKl5vTidVncZSX3AAcADDa5rZmYZaiTQbwBWSFouaYDk5Oq6GcusA05Pn/818MOIiHT6qnRUznJgBfCL9jTdzMwaMWfqJs25nwlcBfQCayPiVknnAUMRsQ74IvCV9GTrDpIfA9LlLiM5cTsOvDsiPGzAzCxHig48sTM4OBhDQ0NFN8PMbI8h6caIGKw3z2MDzcxKzoHezKzkHOjNzErOgd7MrOQ68mSspGHg902uvgDI5r4Hnasb9xm6c7+7cZ+hO/d7d/f50IhYWG9GRwb6Vkgamu3Mc1l14z5Dd+53N+4zdOd+t3OfnboxMys5B3ozs5IrY6BfU3QDCtCN+wzdud/duM/Qnfvdtn0uXY7ezMx2VsYevZmZ1XCgNzMrudIE+rkKmJeFpCWSrpF0m6RbJb03nX6gpKsl3ZX+O7/otrabpF5Jv5R0Zfp6eVqMfmNanH6g6Da2m6R5ki6XdIek2yW9tOzHWtL70//bv5F0iaS9ynisJa2VtF3Sb2qm1T22Snwm3f9bJO1WBZpSBPoGC5iXxTjwgYg4DDgWeHe6r2cBP4iIFcAP0tdl817g9prXHwc+lRalf5CkSH3Z/B/guxHxPOBwkv0v7bGWtAh4DzAYES8kuTX6Ksp5rC8GVs6YNtuxPYWknscKkpKrn9+dNypFoKemgHlEjALVAualExHbIuKm9PmjJF/8RST7+6V0sS8BryukgRmRtBj4C+AL6WsBryApRg/l3OcDgJeR1HsgIkYj4iFKfqxJ6mTsnVar2wfYRgmPdURcR1K/o9Zsx/ZU4MuRuB6YJ+ngRt+rLIG+XgHzukXIy0TSMuBI4AbgGRGxLZ11H/CMotqVkU8DHwQm09cHAQ9FxHj6uozHfDkwDFyUpqy+IGlfSnysI2Ir8I/AvSQB/mHgRsp/rKtmO7YtxbiyBPquI2k/4BvA+yLikdp5aRnH0oyblfQaYHtE3Fh0W3LWBxwFfD4ijgT+xIw0TQmP9XyS3uty4BBgX56a3ugK7Ty2ZQn0XVWEXFI/SZD/akRckU7+Y/VPufTf7UW1LwPHAa+VdA9JWu4VJLnreemf91DOY74F2BIRN6SvLycJ/GU+1icBv4uI4YgYA64gOf5lP9ZVsx3blmJcWQJ9IwXMSyHNTX8RuD0iPlkzq7ZA++nAv+bdtqxExNkRsTgilpEc2x9GxJuBa0iK0UPJ9hkgIu4DNkt6bjrpz0nqL5f2WJOkbI6VtE/6f726z6U+1jVmO7brgLemo2+OBR6uSfHMLSJK8QBeDfwWuBv4UNHtyXA/jyf5c+4W4Ob08WqSnPUPgLuA7wMHFt3WjPb/RODK9PkzgV8AG4GvA5Wi25fB/h4BDKXH+1vA/LIfa+C/A3cAvwG+AlTKeKyBS0jOQ4yR/PX29tmOLSCSkYV3A78mGZXU8Hv5FghmZiVXltSNmZnNwoHezKzkHOjNzErOgd7MrOQc6M3MSs6B3sys5BzozcxK7v8DMHyzogVOU/UAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get inferred taus via modding out by num_taus\n",
    "inferred_taus = np.mod(test_posteriors[0].get_states(),config['num_taus'])\n",
    "plt.plot(inferred_taus[:100])\n",
    "plt.title('inferred time constants')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}