From 004d8550b152b581d2c19dc00d731bf5623d33e5 Mon Sep 17 00:00:00 2001 From: Weinreb Date: Sat, 8 Oct 2022 16:28:13 -0400 Subject: [PATCH] updated notebook --- examples/keypoint_slds.ipynb | 224 ++++++++++++++++++++++++++++------- 1 file changed, 178 insertions(+), 46 deletions(-) diff --git a/examples/keypoint_slds.ipynb b/examples/keypoint_slds.ipynb index bc6189c..a4d2386 100644 --- a/examples/keypoint_slds.ipynb +++ b/examples/keypoint_slds.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "varying-morrison", "metadata": {}, "outputs": [], @@ -29,10 +29,20 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 4, "id": "expressed-christian", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2022-10-08 16:24:18.130080: W external/org_tensorflow/tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory\n", + "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n", + "2022-10-08 16:24:18.130112: W external/org_tensorflow/tensorflow/stream_executor/cuda/cuda_driver.cc:269] failed call to cuInit: UNKNOWN ERROR (303)\n" + ] + } + ], "source": [ "# load dictionary {session_name: ndarray (time,keypoints,2)}\n", "keypoint_data_dict = pickle.load(open('example_keypoint_coords.p','rb'))\n", @@ -44,6 +54,17 @@ "Y,mask = jnp.array(Y),jnp.array(mask)" ] }, + { + "cell_type": "code", + "execution_count": 5, + "id": "convenient-wednesday", + "metadata": {}, + "outputs": [], + "source": [ + "Y = Y[:2]\n", + "mask = mask[:2]" + ] + }, { "cell_type": "markdown", "id": "packed-shade", @@ -54,7 +75,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 6, "id": "contrary-future", "metadata": {}, "outputs": [], @@ -73,26 +94,26 @@ "trans_hypparams = {\n", " 'gamma': 1e3, \n", " 'alpha': 5.7, \n", - " 'kappa': 2e5,\n", + " 'kappa': 1e6,\n", " 'num_states':num_states}\n", "\n", "ar_hypparams = {\n", - " 'nu_0': latent_dim+200,\n", - " 'S_0': 10*jnp.eye(latent_dim),\n", + " 'nu_0': latent_dim+2,\n", + " 'S_0': .01*jnp.eye(latent_dim),\n", " 'M_0': jnp.pad(jnp.eye(latent_dim),((0,0),((nlags-1)*latent_dim,1))),\n", - " 'K_0': 0.1*jnp.eye(latent_dim*nlags+1),\n", + " 'K_0': 10*jnp.eye(latent_dim*nlags+1),\n", " 'num_states':num_states,\n", " 'nlags':nlags}\n", "\n", "obs_hypparams = {\n", - " 'sigmasq_0': .1,\n", - " 'nu_0': keypoint_dim+200,\n", - " 'nu_k': jnp.ones(num_keypoints)*5,\n", - " 'M_0': jnp.zeros((keypoint_dim*num_keypoints, latent_dim+1)),\n", - " 'K_0': jnp.eye(latent_dim+1)*200}\n", + " 'sigmasq_0': 10,\n", + " 'sigmasq_C': .1,\n", + " 'nu_sigma': 1e5,\n", + " 'nu_s': 5,\n", + " 's_0': 1}\n", "\n", "translation_hypparams = {\n", - " 'sigmasq_loc': 0.1\n", + " 'sigmasq_loc': 0.5\n", "}" ] }, @@ -106,7 +127,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 7, "id": "falling-overall", "metadata": {}, "outputs": [], @@ -116,18 +137,48 @@ "states = {}\n", "params = {}\n", "\n", - "states['v'] = initial_location(Y)\n", - "states['h'] = initial_heading(Y, posterior_keypoints, anterior_keypoints)\n", - "states['x'],params['Cd'] = initial_latents(latent_dim=latent_dim, **data, **states)\n", + "states['v'] = initial_location(**data)\n", + "states['h'] = initial_heading(posterior_keypoints, anterior_keypoints, **data)\n", + "states['x'],params['Cd'], pca_model = initial_latents(latent_dim=latent_dim, **data, **states)\n", "\n", "params['betas'],params['pi'] = initial_hdp_transitions(key, **trans_hypparams)\n", "params['Ab'],params['Q']= initial_ar_params(key, **ar_hypparams)\n", - "params['sigmasq'] = initial_variance(**data, **states, **params)\n", + "params['sigmasq'] = jnp.ones(Y.shape[-2])\n", "\n", - "states['z'] = resample_stateseqs(key, **data, **states, **params)\n", + "states['z'],_ = resample_stateseqs(key, **data, **states, **params)\n", "states['s'] = resample_scales(key, **data, **states, **params, **obs_hypparams)\n" ] }, + { + "cell_type": "code", + "execution_count": 8, + "id": "posted-lingerie", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKwAAACICAYAAAB6DAD4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAATUUlEQVR4nO2deXSdVbmHn1/npGM63NB0CLQChZZS20gVEaigtLjkLqVAQREQAbGggoJw9QqiLl14XYJLL5XJ2qVSJvG2FQooDSAIdp5LpdA5HRMS0iHN8N4/9pd4SE/O2Wl7cs6X7Gets8437L3P7zTv2d3f3u/7bpkZgUBc6JRtAYFAawgGG4gVwWADsSIYbCBWBIMNxIpgsIFYkTGDlfSopF2SVrVwX5J+KeltSSskjc+UlkD7IZM97Exgcor7U4ATo9f1wAMZ1BJoJ6Q1WEmFkh6R9Fx0fqqka9PVM7NXgPIURf4TmGWON4B+kgb7Cg90TLp4lJkJ/Bb4bnS+HngceOQoP3sIsCXhfGt0rax5QUnX43ph8vLyJgwbNuywxhoaGujUKXeG5LmmB+Klaf369XvMbNBhN8ws5QtYGL0vTbi2LF29qNzxwKoW7s0Dzko4/xtQkq7NCRMmWDIWLFiQ9Hq2yDU9ZvHSBCyyJH9/n5/bPkkDAAOQ9FGgsnW/o6RsAxK7yqHRtUCgRXyGBLcCc4CRkl4DBgFTj8FnzwFukjQbmAhUmtlhw4FAIJG0BmtmSySdA5wMCHjLzGrT1ZP0GHAuMFDSVuAuoGvU5gzgWeBC4G1gP3DNEX6HQEx5/2Atu/c3tKpOWoOVNB34g5mtjs4LJF1uZv+bqp6ZXZ7mvgHTWyM2EC8aGoxd79ewuXw/m/buY3P5/ujYvZfvO8Sw3p245EL/Nn2GBNeZ2a8bT8ysQtJ1QEqDDXQMaurq2VpxgM17nVFuKt/PlgSjrKn7dw/aSVDUL4/h/fO5YHQhw/v35OCud1v1eT4G21mSoh4RSZ2Bbq36lEDsOXConrU7qli9rZI1ZVW8u2cfm/fup6zqIIkxAHldOzO8fz7HD+zJOScNonhAPsP651M8oCdD+uXRrcsHn/NLS7fQGnwMdj7wuKTfROc3RNcC7ZTKA7Ws2V7F6u2VrN5exaptlWzYXU1DZJh987oyclBPJo4YwPD++Qzvn0/xgHyGD8hnUK/uSMqYNh+D/Q7OSG+Mzl8EHs6YokCbsqe6pskoGw100979TfcL+3RnTFFfpow5jtFD+jJmSF+K+vbIqFGmwmeWoAG3zh/W+mOMmbH3QAMvrN7B6qj3XLWtih1VB5vKDO+fz+iiPlxaMozRRX0YXdSXQb27Z1H14fjMEnwcuBsojsoL95A/IrPSAkfDoboG1pRVsWhjOUs2V7BoYwW73q8BFiPByEG9+OiI/owu6svoIX0YPbgvffO7Zlt2WnyGBI8AtwCLgfrMygkcKRX7DjnD3FTB4k0VLN/yXtMT+tCCPD42cgC9a/bwuXNLOGVwb/K7+fzpcw8f1ZVm9lzGlQS8MTPe2bOPxRudcS7aVM6G3fsA6NJJjB7Sly9MLKbk+AImFBdQ2KcHAKWlpUwoLsim9KPGx2AXSPoZ8CegpvGimS3JmKrABzhYW8+KrZUs2lTOkqgHrdjvFhv75nVlQnEBnx8/lJLiAsYO7Udet85ZVpw5fAx2YvReknDNgE+mqyhpMnA/0Bl42Mx+2ux+MfAozj+hHPiimW310NTu2f1+Dc+tKmPe8jKWbqmgtt7NKY0Y1JPzTyls6j1HDOxFp07ZeWLPBj6zBJOOpOFogeHXwKdwvq4LJc0xszUJxf4H58T9O0mfBH4CXHkkn9ceqNxfy/OrdzBn+XZe37CHBoNRx/Xm2rNGUFJcwPjiAvr37NhrNl4jb0mfAUYDPRqvmdk9aaqdAbxtZu9EbczGRRkkGuypOG8wgAXAn71UtyP21dTx17U7mbu8jJfX76K23igekM/0SR/is6cXcVJh72xLzCl8prVmAPnAJNyCwVTgnx5tJ4somNiszHLg87hhw+eA3pIGmNneZhqaIg4KCwspLS097MOqq6uTXs8WqfQcqjdW7qnnzbI6lu2u51A9FHQX5w3rzMTBXTi+D0hlbF9bxva1baMpW7RaUzKvbvtgZMCKZu+9gFc96k3FjVsbz68EftWsTBHuYW4pzmi3Av1StRvXiINDdfVW+tYu+9YTy2zM9+db8Xfm2fh7XrDvPbPS3nxnr9XXN7S5plygtREHPkOCA9H7fklFwF7AJ1gwbUSBmW3H9bBI6gVcbGbvebQdCxoajIUby5m7YjvPrtxB+b5D9O7ehQvGHMdFpxdx5sgBdOmcWzFWuY6Pwc6T1A/4GbAEN0Pg40uwEDhR0gk4Q50GXJFYQNJAoNzc8u+duBmD2FNWeYDH1tVwx+svsaPqID26duL8Uwq56PQizj5pED26tt9pp0zjM0vww+jwaUnzgB5mljamy8zqJN0EPI+b1nrUzFZLugfX3c/BRST8RJIBrxBzh+4Dh+r5zSsbmPHyBmrrGpg0agB3XjiK808ppGf3eK4s5Rot/itK+qSZvSTp80nuYWZ/Ste4mT2LC4VJvPb9hOOngKdaJzn3MDPmLN/OT59bR1nlQT4zdjDn9nuPSy4sSV850CpS/ezPAV4CPpvknuEeljo8SzdXcM+8NSzd/B5jhvTh/mkf5owT+ufc03h7oUWDNbO7JHUCnjOzJ9pQUywoqzzAvfPf4pml2xjUuzs/mzqWi8cP7VCrTtkg5cDKzBok3Q4Eg41IHKc2GEyfNJIbz/0QvcIYtU3w+Vf+q6Rv49IT7Wu8aGap8ma1Ow4bp542mDumjGJY//xsS+tQ+BjsZdF74hO8AR3GgbulcWqg7fGZ1jqhLYTkIs3HqfdOHcvUME7NKr7OL2NwjiqJzi+zMiUq2xw4VM+Dr7zDjJc3UG8Wxqk5hI/zy124Cf5TcXOqU4C/A+3OYMM4Nffx6TKmAqfj0m1eI6kQ+L1P4x4O3MOB3wH9ojJ3RIsNbU5dfQM3/mEJL67ZyeiiPtx32TgmjhiQDSmBFHg5v0TTW3WS+gC7+KBTS1I8Hbi/BzxhZg9IauzBj2/tlzgW/Ogva3lxzU7unDKKr3xiBJ3DODUn8THYRZHzy0O4yNlq4B8e9XwcuA3oEx33Bbb7yT62zPrHRma+vpFrzzqBG84ZmQ0JAU9krdgcWdLxQB8zW+FRdiow2cy+Ep1fCUw0s5sSygwGXgAKgJ7A+Wa2OElbiQ7cE2bPnn3Y51VXV9OrVy/v79LIit11/GJxDacP6szXx3en0zHKaHKkejJJnDRNmjRpsZkd7oyRzEnWPuhkPQfnFtgzXdlm9XwcuG8FvhUdfwzX+3ZK1e6xdOB+a0eVjfn+fJt83ytWfbC21fVTESdn6WySiZTxPwfOAtZIekrSVEk90lXCLyX8tUTLvmb2D9y02UCPto+aPdU1fHnmQvK6deaRq0qC+19MSGuwZvaymX0Nt7L1G+BS3INXOpocuCV1wzlwz2lWZjNwHoCkU3AGu9tf/pFxsLae62ctYk91DQ9fVUJRv7xMf2TgGOG7cJCHczO8DBiPm4pKifk5cH8LeEjSLbgHsKuj/w4yhplx+1MrWLL5PWZ8cTxjh/bL5McFjjE+CwdP4J745wO/Al42F9KSFkvvwL0G+HhrBB8t9/31X8xZvp3bJ5/M5DFhH7u44ZsM7nIzi30iuP9bto37//YvLpkwlBvD9FUs8XF+eb4thGSaxZvKue3JFUw8oT8//txpWUvIGzg6OkSM8Zby/Vw/azFDCvKY8cUJh+XZD8SHdv+XqzpYy5dnLqSuwXjkqhIKOnhuqriTKmp2fKqKFoN0m3X1DUz/wxLe3bOPWdeewYhBubXKE2g9qcawP4/ee+BSbS7HpYsfCyzCrUzlLGbG3XNX8+q/9nDvxWM5c2SbrEcEMkyLQwIzm2Qu1WYZMN7MSsxsAvBhYrCJ8W9f28jv39jMDeeM4NKPpHUuC8QEnzHsyWa2svHEzFYBp2RO0tHz0rqd/Ogva7hgdCHfuWBUtuUEjiE+87ArJD3Mv522vwCk9dbKFmvLqrj5j0s5tagPv7hsXIi/amf49LDXAKuBb0SvNXjuvC1psqS3JL0t6Y4k938haVn0Wi/pvVZoP4xdVQe5duZCevfoyiNXfSS2O6UEWsZn4eBglNT4WTN7y7dhn4gDM7slofzNuPHxEVFTb1w3axEV+2t58qsfa9o5JdC+SNvDSroIWEa0v6ykcZKae10loyniwMwOAY0RBy1xOfCYR7uH0dBgPLSihhXbKrl/2jjGDOl7JM0EYoDP/5l34YyvFMDMlkU5X9PhkzIeaNpN5gRc8rlk91OmjH96/SEW7aznspO70W33OkpL13nIyyztIj17G9BaTT4GW2tmlc3W3o+1C+A04KmWHGzM7EHgQYCSkhI799xzm+4t3VzB3Pmvc/bQLvz06vNzxkegtLSURJ25QHvQ5GOwqyVdAXSWdCLwdeB1j3o+EQeNTOMIkxl/eHgBD3xhPF13r8sZYw1kDp9ZgptxWx7V4MaYVcA3Per5RBwgaRQuCNEnEjcpU04bTJcwfdUh8Jkl2A98N3p54xlxAM6QZ2c60iDQPvCJODgJ+DYuwUVTeTNLu3VnuoiD6PxuP6mBgN8Y9klgBm7nmNhHHQTijY/B1pnZAxlXEgh44PPQNVfS1yQNltS/8ZVxZYFAEnx62Kui99sSrnWoDNyB3CFk4A7EiiPaWA7w2lguEDjWhI3lArEi5cZy0buX72sg0Bb45tb6DG55NnFTjns86qVMGR+VuRS4G9drLzezK5qXCQQa8VnpmgHkA5NwiwdTgX961EvrwB0509wJfNzMKiT9xxF9i0CHwWce9kwz+xJQYWY/wIV3n+RRz8eB+zrg12ZWAWBmPmk8Ax0Yr005ovf9koqAvYBP2j8fB+6TACS9hhs23G1m85s3lM6BG3LPOTnX9EA70ZQsLbd9MK37f+O2JboY2IHLU/BDj3o+KePnAc8AXXERB1uAfqnaPZYp4zNJrukxi5cmWkgZ77Nw8MPo8GlJ84AeZlbp8VvwceDeCrxpZrXAu5LWAyfifGkDgcNItXCQdMEguuezcNDkwI0z1Gm4zT0S+TMu+PC3kgbihgjveOgOdFBS9bDJFgwaSbtwYH4O3M8Dn5a0Bue6eJuZ7W3VNwh0KFItHBz1goGlTxlvuK2Pbj3azwp0DHzyEgyQ9EtJSyQtlnS/pLAJayAr+MzDzsZtRXQx7sl/N/B4JkUFAi3hMw87OGGmAOBHki7LlKBAIBU+PewLkqZJ6hS9LsU9LAUCbY6PwV4H/BGXl6AGN0S4QdL7kqoyKS4QaI7PwkHvthASCPjgM0twbbPzzpLuypykQKBlfIYE50l6NoqaHQO8AYReN5AVfIYEV0SzAiuBfcAVZvZaxpUFAknwGRKciEsV/zSwCbhSUr5P4x4p46+WtDshbfxXWvsFAh0Ln3nYucB0M/ubXD7LW3GOLaNTVfKJOIh43Mxuar30QEfEx2DPMLMqaFr7/7mkuT71iCIOACQ1Rhw0N9hAwJtU7oW3m9m9ZlYl6RIzezLh9tXAf6Vp2zdl/MWSzgbWA7eY2ZbmBRIjDoBqSck2BxkI7EmjqS3JNT0QL03FSUsn8+p2HSlLkh0nO2+hvk/EwQCge3R8A/BSunZTfF5SD/VsvXJNT3vRlOqhSy0cJztPRtqIAzPba2Y10enDwASPdgMdmFQGay0cJztPRtqU8ZISgxkvAtZ6tBvowKR66Do98hUQkJfgNyASEmq0hPlFHHw92gesDijHjY2PlAePom4myDU90A40KRpHBAKxwGdpNhDIGYLBBmJF7A1W0qOSdklalW0tAJJ6SPqnpOWSVkv6QbY1AUjaKGlltAS+KAf0nJywJL9MUpWkb6atF/cxbLToUA3MMrMxOaBHQE8zq5bUFfg78A0zeyPLujYCJWaWawsHjcv424CJZrYpVdnY97Bm9gpuhiEnMEd1dNo1esW7V8g85wEb0hkrtAODzUUiJ/dlwC7gRTN7M8uSwP1oXohC9a9PW7ptmYbbFjYtwWAzgJnVm9k43OreGZHje7Y5y8zGA1OA6dFQKutEi0oX4TYwTEsw2AxiZu8BC4DJWZaCmW2L3nfhMkaekV1FTUzB+abs9CkcDPYYI2mQpH7RcR7OH3hdljX1lNS78Rj4NJATsyq4ZIBewwFoBwYr6THc1vUnS9raPGgyCwwGFkhagfOneNHM5mVZUyHwd0nLcen+/2JJEke3NdGP51O0Ykei2E9rBToWse9hAx2LYLCBWBEMNhArgsEGYkUw2ECsCAabA0iqjzyWVkl6sjFRiaTjJM2WtCFaUn1Wks+mfu2WYLC5wQEzGxd5mx0Cvhp5fT0DlJrZSDObgNvmtDCbQrON1+bIgTblVWAsbm/fWjOb0XjDzJZDU/Dm40Af3N/wRjN7NQta25zQw+YQkrrg1tZXAmOAxS0UvQJ4PnKwOR1Y1hb6coFgsLlBXuSOuAjYDDySpvxC4BpJdwOnmdn7mZWXOwSDzQ0ax7DjzOxmc7ufr6aFxCKR0/rZOC/9mZK+1IZas0ow2NzlJaB7orO1pLGSPiGpGNhpZg/hMuaMz5bItiY4v+QAkqrNrFeS60XAfbie9iCwEfgmcCZwG1CLi2f7kpm920Zys0ow2ECsCEOCQKwIBhuIFcFgA7EiGGwgVgSDDcSKYLCBWBEMNhAr/h+fG0wgWoHkHAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(np.arange(latent_dim)+1,np.cumsum(pca_model.explained_variance_ratio_))\n", + "plt.xlabel('PCs')\n", + "plt.ylabel('Explained variance')\n", + "plt.yticks(np.arange(0.5,1.01,.1))\n", + "plt.xticks(range(1,latent_dim+2,2))\n", + "plt.gcf().set_size_inches((2.5,2))\n", + "plt.grid()\n", + "plt.tight_layout()" + ] + }, { "cell_type": "markdown", "id": "million-reynolds", @@ -138,19 +189,69 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "capable-grove", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "260319e619e84d6e872563d085e5d153", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/500 [00:00" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAADQCAYAAABGO9SNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtNUlEQVR4nO3deZwdVZn/8c+XsO+QRAazkABhCSABGgiiyE5AICDIPgR+OBlHEFAZJzgKDMgIKG6oaJTIvgkiGQ1iCERcCKQDkQgYiCFIYiCBsO8Jz++POjcUze3u252+t7q6v+/X67666lSdqqf6duo+OffUOYoIzMzMzMysMVYqOgAzMzMzs97ECbiZmZmZWQM5ATczMzMzayAn4GZmZmZmDeQE3MzMzMysgZyAm5mZmZk10MpFB9Bo/fr1iyFDhhQdhplZh82YMeO5iOhfdByN1N3v2bNnzwZgyy23LDgSM+tu2rpn97oEfMiQITQ3NxcdhplZh0l6qugYGq2737P33HNPAKZOnVpoHGbW/bR1z3YXFDMz6xBJoyTNljRH0rgq2z8raZakmZL+KGl4btvZqd5sSQc0NnIzs+7BCbiZmdVMUh/gh8CBwHDg2HyCnVwfEdtFxAjgEuDbqe5w4BhgG2AU8KN0PDOzXqXXdUExM7MVsgswJyLmAki6ERgNPFrZISJezu2/FhBpeTRwY0S8BTwpaU463n2NCLwe3PXEzDrDCbiZmXXEAODp3Pp8YNeWO0k6FfgisCqwd67utBZ1B1SpOxYYCzB48OAuCdrMrDtxFxQzM+tyEfHDiNgM+C/gqx2sOz4imiKiqX//7j3oy7e+9S2+9a1vFR2GmZWMW8BrNGTcbwCYd9EnC47EzKxQC4BBufWBqaw1NwKXd7JuXVXu621p757/61//GoCzzjqrS2Iys97BLeBmZtYR04FhkoZKWpXsocqJ+R0kDcutfhJ4Ii1PBI6RtJqkocAw4IEGxGxm1q24BdzMzGoWEUslnQbcCfQBJkTEI5LOB5ojYiJwmqR9gXeAF4Axqe4jkm4me2BzKXBqRCwr5ELMzArkBNzMzDokIiYBk1qUnZNbPqONuhcCF9YvOjOz7s8JuJmZWSetscYaRYdgZiXkBNzMzKyT7rjjjqJDMLMS8kOYZmZmZmYN5ATczMysky644AIuuOCCosMws5JxAm5mZtZJU6ZMYcqUKUWHYWYl4wTczMzMzKyBnICbmZmZmTWQE3AzMzMzswbyMIRmZmad1Ldv36JDMLMScgJuZmbWSbfeemvRIZhZCbkLipmZmZlZAzkBNzMz66Szzz6bs88+u+gwzKxk3AXFzMysk+67776iQzCzEnILuJmZmZlZAzkBNzMzMzNrICfgZmZmZmYN5D7gZmZmnTRw4MCiQzCzEqpbC7ikQZLukfSopEcknZHKN5Q0WdIT6ecGqVySvi9pjqSHJe2YO9aYtP8TksbkyneSNCvV+b4k1et6zMwMJI2SNDvdd8dV2f7FdN9/WNIUSZvkti2TNDO9JjY28vq49tprufbaa4sOw8xKpp5dUJYCX4qI4cBI4FRJw4FxwJSIGAZMSesABwLD0msscDlkCTtwLrArsAtwbiVpT/v8W67eqDpej5lZryapD/BDsvv1cODYdF/PewhoioiPALcAl+S2vRERI9Lr0IYEbWbWDdUtAY+IhRHxYFp+BXgMGACMBq5Ku10FHJaWRwNXR2YasL6kjYEDgMkRsSQiXgAmA6PStnUjYlpEBHB17lhmZtb1dgHmRMTciHgbuJHs3r1cRNwTEa+n1WlAj+6jceaZZ3LmmWcWHYaZlUxD+oBLGgLsANwPbBQRC9OmZ4CN0vIA4OlctfmprK3y+VXKzcysPqrdj3dtY/9TgDty66tLaib7hvSiiPhVl0fYYDNnziw6BDMrobon4JLWBm4FzoyIl/PdtCMiJEUDYhhL1q2FwYMH1/t0Zma9nqQTgCbgE7niTSJigaRNgbslzYqIv1ep63u2mfVodR2GUNIqZMn3dRHxy1T8bOo+Qvq5KJUvAAblqg9MZW2VD6xS/gERMT4imiKiqX///it2UWZmvVdr9+P3kbQv8N/AoRHxVqU8Ihakn3OBqWTfjH6A79lm1tPVcxQUAVcAj0XEt3ObJgKVkUzGALfnyk9Mo6GMBF5KXVXuBPaXtEF6+HJ/4M607WVJI9O5Tswdy8zMut50YJikoZJWBY4hu3cvJ2kH4CdkyfeiXPkGklZLy/2A3YFHGxa5mVk3Us8uKLsD/wrMkjQzlX0FuAi4WdIpwFPAUWnbJOAgYA7wOnAyQEQskXQB2Y0f4PyIWJKWPwdcCaxB1s8w39fQzMy6UEQslXQaWcNIH2BCRDwi6XygOSImAt8E1gZ+kboc/iONeLI18BNJ75I1/lwUEaVPwLfYYouiQzCzEqpbAh4RfwRaG5d7nyr7B3BqK8eaAEyoUt4MbLsCYZqZ9UqStouIWR2tFxGTyBpM8mXn5Jb3baXen4HtOnq+7m78+PFFh2BmJeSp6M3MeqcfSXpA0uckrVd0MGZmvYkTcDOzXigiPg4cT/ZQ5QxJ10var+CwSmfs2LGMHTu26DDMrGQaMg64mZl1PxHxhKSvAs3A94Ed0kPtX8mNXGVtePzxx4sOwcxKyC3gZma9kKSPSPoO2SzFewOHRMTWafk7hQZnZtbDuQXczKx3ugz4GVlr9xuVwoj4Z2oVNzOzOnELuJlZ73RbRFyTT74lnQEQEdcUF5aZWc/nBNzMrHc6sUrZSY0OouxGjBjBiBEjig7DzErGXVDMzHoRSccCxwFDJeVnsVwHWFK9lrXmu9/9btEhmFkJOQE3M+td/gwsBPoBl+bKXwEeLiQiM7NepuYEXNKaEfF6PYMxM7P6ioingKeA3YqOpSc44YQTALj22msLjsTMyqTdPuCSPirpUeBvaX17ST+qe2RmZtblJP0x/XxF0su51yuSXi46vrKZP38+8+fPLzoMMyuZWlrAvwMcAEwEiIi/SNqjrlGZmVldRMTH0s91io7FzKy3qmkUlIh4ukXRsjrEYmZmDSJpM0mrpeU9JZ0uaf2CwzIz6xVqScCflvRRICStIuksspnTzMysvG4FlknaHBgPDAKuLzYkM7PeoZYuKJ8FvgcMABYAvwNOrWdQZmZWd+9GxFJJhwOXRcRlkh4qOqiy2W03P8tqZh3XbgIeEc8BxzcgFjMza5x30pjgY4BDUtkqBcZTSt/4xjeKDsHMSqjdBFzS96sUvwQ0R8TtXR+SmZk1wMlk33BeGBFPShoKeAp6M7MGqKUP+OrACOCJ9PoIMBA4RdJ36xaZmZnVTUQ8GhGnR8QNaf3JiLi4vXqSRkmaLWmOpHFVtn9R0qOSHpY0RdImuW1jJD2RXmO69oqKccQRR3DEEUcUHYaZlUwtfcA/AuweEcsAJF0O/AH4GDCrjrGZmVmdSNodOA/YhOyzQEBExKZt1OkD/BDYD5gPTJc0MSIeze32ENAUEa9L+g/gEuBoSRsC5wJNQAAzUt0Xuv7qGuf5558vOgQzK6FaEvANgLXJup0ArAVsGBHLJL1Vt8jMzKyergC+AMyg9qFldwHmRMRcAEk3AqOB5Ql4RNyT238acEJaPgCYHBFLUt3JwCjghhW4hrobMu43bW5/Zu7zjNy0b4OiMbOeopYE/BJgpqSpZC0kewD/K2kt4K46xmZmZvXzUkTc0cE6A4D8vBDzgV3b2P8UoHKOanUHdPD8ZmY9Qi2joFwhaRJZywfAVyLin2n5P+sWmZmZ1dM9kr4J/BJY/m1mRDzYFQeXdAJZd5NPdKLuWGAswODBg7siHDOzbqWWFnCAN4GFZA9kbi5p84i4t35hmZlZnVVarptyZQHs3UadBWQT9lQMTGXvI2lf4L+BT0TEW7m6e7aoO7XaSSJiPNnkQDQ1NUUb8RRu9U22Z599tiw6DDMrmVqGIfwMcAbZzXImMBK4j7Zv0mZm1o1FxF6dqDYdGJaGLFwAHAMcl99B0g7AT4BREbEot+lOsu6LG6T1/YGzOxFDt7L+7sfyta99sugwzKxkahmG8AxgZ+CpdMPeAXixvUqSJkhaJOmvubLzJC2QNDO9DsptOzsNazVb0gG58qpDXkkaKun+VH6TpFVru2QzM5O0kaQrJN2R1odLOqWtOhGxFDiNLJl+DLg5Ih6RdL6kQ9Nu3yR7cP8X6T4/MdVdAlxAlsRPB86vPJBpZtbb1JKAvxkRbwJIWi0i/gbU8n3blWRPuLf0nYgYkV6T0nGHk7WkbJPq/EhSn9yQVwcCw4Fj074AF6djbQ68QPawj5mZ1eZKskT6w2n9ceDM9ipFxKSI2CIiNouIC1PZORFRSbT3jYiNcvf5Q3N1J0TE5un1866+oCI8e/O5HHjggUWHYWYlU0sCPl/S+sCvgMmSbgeeaq9S6iNea+vGaODGiHgrIp4E5pA99Ll8yKuIeBu4ERgtSWRdYG5J9a8CDqvxXGZmBv0i4mbgXVjeul3rcISWxNK3eOONN4oOw8xKppZRUA5Pi+dJugdYD/jtCpzzNEknAs3Al9IkDAPIxoutyA9PVW3Iq77Ai+kDo+X+ZmbWvtck9SV78BJJI3lvvgczM6ujdlvAJW0mabXKKjAEWLOT57sc2IxsavuFwKWdPE6HSBorqVlS8+LFixtxSjOz7u6LwERgM0l/Aq4GPl9sSGZmvUMtXVBuBZZJ2pxsWKhBwPWdOVlEPBsRyyLiXeCnvDe2eGtDW7VW/jywvqSVW5S3dt7xEdEUEU39+/fvTOhmZj1KGu/7E8BHgX8HtomIh4uNysysd6hlHPB3I2KppMOByyLiMkkPdeZkkjaOiIVp9XCgMkLKROB6Sd8meyBoGPAAWYv7B4a8iohI3WGOJOsXPga4vTMxmZn1JpI+1cqmLSQREb9saEAlt8Zmu3DwJ7cuOgwzK5laEvB3JB1LluQekspWaa+SpBvIJl3oJ2k+cC6wp6QRZH0O55G1upCGsboZeBRYCpwaEcvScSpDXvUBJkTEI+kU/wXcKOnrwEPAFTVci5lZb1e5j3+IrPX77rS+F/BnspkxrUbr7fopzjrL44CbWcfUkoCfDHwWuDAinkyt0de0Vykijq1S3GqSnIazurBK+SRgUpXyubzXhcXMzGoQEScDSPodMLzyraSkjcmGJjQzszprtw94RDwaEadHxA1p/cmIuLj+oZmZWR0NynUJBHgWGFxUMGX1zPXj2HPPPYsOw8xKppap6J8kDVOVFxGb1iUiMzNrhCmS7gRuSOtHA3cVGI+ZWa9RSxeUptzy6sCngQ3rE46ZmTVCRJyWHq7fIxWNj4jbiozJzKy3qGUinudbFH1X0gzgnPqEZGZmjZASbifdZmYNVksXlB1zqyuRtYjX0nJuZmZmZmYt1JJI52erXAo8CRxVn3DMzMzKY62tPs5Rh21bdBhmVjK1dEHZqxGBmJlZ40g6BPhNmpnYOmmdHT/J5z7nccDNrGNqmYrezMx6nqOBJyRdImmrooMpq3ffeZPXX3+96DDMrGScgJuZ9UIRcQKwA/B34EpJ90kaK2mdgkMrlUW/OI+DDjqo6DDMrGScgJuZ9VIR8TJwC3AjsDFwOPCgpM8XGpiZWQ/XbgIuaU1JX5P007Q+TNLB9Q/NzMzqRdJoSbcBU4FVgF0i4kBge+BLRcZmZtbT1dIC/nPgLWC3tL4A+HrdIjIzs0b4FPCdiNguIr4ZEYsAIuJ14JTWKkkaJWm2pDmSxlXZvoekByUtlXRki23LJM1Mr4ldfUFmZmVRSwK+WURcArwDy2/OqmtUZmZWb89ExL35AkkXA0TElGoVJPUBfggcCAwHjpU0vMVu/wBOAq6vcog3ImJEeh26gvGbmZVWLQn425LWAAJA0mZkLeJmZlZe+1UpO7CdOrsAcyJibkS8TdZ3fHR+h4iYFxEPA71ieMO1t9uXk046qegwzKxkaknAzwV+CwySdB0wBfhyXaMyM7O6kPQfkmYBW0l6OPd6Eni4neoDgKdz6/NTWa1Wl9QsaZqkw9qIcWzar3nx4sUdOHzjOQE3s86oZSKeyZIeBEaSdT05IyKeq3tkZmZWD9cDdwDfAPJ9uF+JiCV1PvcmEbFA0qbA3ZJmRcTfW+4UEeOB8QBNTU1R55hWyLLXX+K5556jX79+RYdiZiXSagu4pB0rL2ATYCHwT2BwKjMzs/KJiJgHnAq8knshacN26i4ABuXWB6ayWk+8IP2cSzb6yg611u2uFv/qGxx55JHt72hmltNWC/ilbWwLYO8ujsXMzOrveuBgYAbZvTz/UH0Am7ZRdzowTNJQssT7GOC4Wk4qaQPg9Yh4S1I/YHfgko6Hb2ZWfq0m4BGxVyMDMTOz+ouIg9PPoZ2ou1TSacCdQB9gQkQ8Iul8oDkiJkraGbgN2AA4RNL/RMQ2wNbATyS9S/bt60UR8WgXXZaZWam02wdc0urA54CPkbWO/AH4cUS8WefYzMysi7XXhTAiHmxn+yRgUouyc3LL08m6prSs92dguw4Fa2bWQ7WbgANXk/UPvCytHwdcA3y6XkGZmVnduHuhmVnBaknAt42I/EQL90jy14ZmZiXk7oVda50dDuI/jvO4BGbWMbUk4A9KGhkR0wAk7Qo01zcsMzOrB0l7R8Tdkj5VbXtE/LLRMZXZWlvvwdFHf7LoMMysZFpNwNNEDQGsAvxZ0j/S+ibA3xoTnpmZdbFPAHcDh1TZFoAT8A5Y+vJinn76aQYNGtT+zmZmSVst4AevyIElTUjHWBQR26ayDYGbgCHAPOCoiHhBkoDvAQcBrwMnVR4EkjQG+Go67Ncj4qpUvhNwJbAG2QNBZ0REt56wwcysaBFxbvp5ctGx9ATP/fpS/vXRK5k6dWrRoZhZibQ6EU9EPJV/AW+QtY5UXu25EhjVomwcMCUihpFNaV+Zhe1AYFh6jQUuh+UJ+7nArsAuwLlpLFnSPv+Wq9fyXGZm1gpJfSV9X9KDkmZI+p6kvkXHZWbWG7SagFdIOlTSE8CTwO/JWq7vaK9eRNwLtJzWeDRwVVq+CjgsV351ZKYB60vaGDgAmBwRSyLiBWAyMCptWzcipqVW76tzx6qrIeN+s/xlZlZiNwKLgSOAI9PyTYVGZGbWS7SbgAMXACOBx9PEDfsA0zp5vo0iYmFafgbYKC0PAJ7O7Tc/lbVVPr9KeVWSxkpqltS8ePHiToZuZtajbBwRF0TEk+n1dd67J5uZWR3VkoC/ExHPAytJWiki7gGaVvTEqeW6IX22I2J8RDRFRFP//v0bcUozs+7ud5KOkbRSeh1FNsOlmZnVWS3DEL4oaW3gXuA6SYuA1zp5vmclbRwRC1M3kkWpfAGQf4R8YCpbAOzZonxqKh9YZX8zM2uDpFfIGj8EnAlcmzatBLwKnFVMZOW07i6H86UxOxcdhpmVTC0t4KPJHsD8AvBb4O9UH76qFhOBMWl5DHB7rvxEZUYCL6WuKncC+0vaID18uT9wZ9r2sqSRaQSVE3PHMjOzVkTEOhGxbvq5UkSsnF4rRcS6RcdXNmtuviuHHNLZj0Qz663abQGPiHxr91Wt7tiCpBvIWq/7SZpPNprJRcDNkk4BngKOSrtPIhuCcA7ZMIQnp3MvkXQBMD3td35EVB7s/BzvDUN4BzU8GGpmZu9JDRvDgNUrZekBeqvRO8/PZ/bs2Wy55ZZFh2JmJdLWRDyVryk/sImsC3ebLSURcWwrm/apsm8Ap7ZynAnAhCrlzcC2bcVgZmbVSfoMcAZZF76ZZA/b3wfsXWBYpfP8nT9g+zt/wL8cd1G7+867yDNmmlmmrXHAK19Ttnyt468pzcxK7wxgZ+CpiNgL2AF4sdCIzMx6ibZawDdsq2KuK4iZmZXPmxHxpiQkrRYRf5PkfhRmZg3QVh/wGbz3pHxLAWxal4jMzKwR5ktaH/gVMFnSC2TP5piZWZ21moCnSXfMzKwHiojD0+J5ku4B1iMb6crMzOqslqnob5V0kKRahiw0M7OSkLSjpNOBjwDzI+LtGuqMkjRb0hxJ46ps30PSg5KWSjqyxbYxkp5IrzEt65bReh89hvU+ekzRYZhZydSSVF8OHA88Ieki9xE0Mys/SeeQDS3bF+gH/FzSV9up0wf4IXAgMBw4VtLwFrv9AzgJuL5F3Q3JhqPdFdgFODcNg1hqawwZwRpDRhQdhpmVTLsJeETcFRHHAzsC84C7JP1Z0smSVql3gGZmVhfHAztHxLkRcS7ZMIT/2k6dXYA5ETE3tZbfSDZZ23IRMS8iHgbebVH3AGByRCyJiBeAycCorriQIr397FzefnZu0WGYWcnU1K1EUl+yFo3PAA8B3yNLyCfXLTIzM6unf5KbgAdYDVjQTp0BwNO59fmprBY115U0VlKzpObFixfXePhiLJkyniVTxhcdhpmVTLszYUq6DdgSuAY4JE0DD3CTpOZ6BmdmZl1L0mVkI1m9BDwiaXJa3w94oMjYKiJiPDAeoKmpqdqEcGZmpdZuAg58PyLuqbYhIpq6OB4zM6uvSsPJDOC2XPnUGuouAAbl1gfSfqt5vu6eLerWck4zsx6nrYl4dgaeriTfkk4EjiAbJ/Y8T8RjZlY+EXFVZVnSqsAWaXV2RLzTTvXpwDBJQ8kS6mOA42o89Z3A/+YevNwfOLvmwM3MepC2+oD/BHgbsmGlgIuAq8m+tnSHNzOzEpO0J/AE2agmPwIeT/f6VkXEUuA0smT6MeDmiHhE0vmSDk3H3VnSfODTwE8kPZLqLgEuIEvipwPnuyHHzHqrtrqg9MndHI8GxkfErcCtkmbWPTIzM6unS4H9I2I2gKQtgBuAndqqFBGTgEktys7JLU8n615Sre4EYMKKhd29rL9HjxjO3MwarM0EXNLKqcVjH2BsjfXMzKz7W6WSfANExOMeWrbjVh+4ddEhmFkJtZVI3wD8XtJzwBvAHwAkbU7WDcXMzMprhqSfAdem9eN57wFNq9Gb8x8DnIibWce0moBHxIWSpgAbA7+LiMpQUCsBn29EcGZmVjefBU4FTk/rfyDrC24d8OK92TOt/3LcRQVHYmZl0mZXkoiYVqXs8fqFY2Zm9ZamlP9LRGwFfLvoeMzMepuaZsI0M7OeIyKWAbMlDS46FjOz3sgPU5qZ9U4bkM2E+QDwWqUwIg4tLiQzs97BCbiZWe/0taIDMDPrrZyAm5n1IpJWJ3sAc3NgFnBFGm7WOmHDfca2v5OZWQtOwM3MepergHfIRj05EBgOnFFoRCW26kabFh2CmZWQE3Azs95leERsByDpCuCBguMptTfmzQRgjSEjCo3DzMrFCbiZWe/yTmUhIpZKKjKW0nvpzzcCtSXgQ8b9pt195l30yRUNycxKoJBhCCXNkzRL0kxJzalsQ0mTJT2Rfm6QyiXp+5LmSHpY0o6544xJ+z8haUwR12JmVjLbS3o5vV4BPlJZlvRy0cGZmfUGRY4DvldEjIiIprQ+DpgSEcOAKWkdsj6Kw9JrLHA5ZAk7cC6wK7ALcG4laTczs+oiok9ErJte60TEyrnldYuOz8ysN+hOE/GMJns4iPTzsFz51ZGZBqwvaWPgAGByRCyJiBeAycCoBsdsZmZmZtYhRSXgAfxO0gxJlTGcNoqIhWn5GWCjtDwAeDpXd34qa638AySNldQsqXnx4sVddQ1mZmZmZh1W1EOYH4uIBZI+BEyW9Lf8xogISdFVJ4uI8cB4gKampi47rpmZ9W59Dzit6BDMrIQKaQGPiAXp5yLgNrI+3M+mriWkn4vS7guAQbnqA1NZa+VmZmYNsUrfgazSd2DRYZhZyTQ8AZe0lqR1KsvA/sBfgYlAZSSTMcDtaXkicGIaDWUk8FLqqnInsL+kDdLDl/unMjMzqxNJoyTNTiNTjauyfTVJN6Xt90saksqHSHojjX41U9KPGx58Hbw+535en3N/0WGYWckU0QVlI+C2NPbsysD1EfFbSdOBmyWdAjwFHJX2nwQcBMwBXgdOBoiIJZIuAKan/c6PiCWNuwwzs95FUh/gh8B+ZM/dTJc0MSIeze12CvBCRGwu6RjgYuDotO3vETGikTHX28sP3AbAmpvv2iXH81jhZr1DwxPwiJgLbF+l/HlgnyrlAZzayrEmABO6OkYzM6tqF2BOuo8j6UaykaryCfho4Ly0fAvwA3m2HzOz9+lOwxCamVn3VsvoU8v3iYilwEtA37RtqKSHJP1e0sfrHayZWXflqejNzKwRFgKDI+J5STsBv5K0TUR8YPbNNDztWIDBgwc3OEwzs/pzC7iZmdWqltGnlu8jaWVgPeD5iHgrdTUkImYAfwe2qHaSiBgfEU0R0dS/f/8uvgQzs+K5BdzMzGo1HRgmaShZon0McFyLfSojWt0HHAncneZ26A8siYhlkjYFhgFzGxd6ffQ7+EtFh2BmJeQE3MzMahIRSyWdRjbkax9gQkQ8Iul8oDkiJgJXANdImgMsIUvSAfYAzpf0DvAu8NmeMHLVyus2voXeI6WYlZ8TcDMzq1lETCIbHjZfdk5u+U3g01Xq3QrcWvcAG+y1x+4FYK2t9yg4EjMrEyfgZmZmnfTKQ9n/RZyAm1lH+CFMMzMzM7MGcgJuZmZmZtZATsDNzMzMzBrICbiZmZmZWQP5IUwzM7NO6n/Y2UWHUFUtQxWChys0K4oTcDMzs07qs+Z6RYdgZiXkLihmZmad9Oqsu3h11l1Fh2FmJeMWcDMzs06qJN9rb7dvwZF0jmfVNCuGW8DNzMzMzBrICbiZmZmZWQO5C4qZmZm1yt1UzLqeW8DNzMzMzBrILeBmZmad9KFPn1d0CKXhlnSz9zgBNzMz66SVVlm96BC6hVon/jGzjBNwMzOzTnrlwSzxXGdHt9x2Bc/gab2FE3AzM7NOeu1vfwCcgDdad2xx938KyqmorlGlT8AljQK+B/QBfhYRFxUckplZj9befVfSasDVwE7A88DRETEvbTsbOAVYBpweEXc2MHSzuumqRK6r/nPh/xB0b6VOwCX1AX4I7AfMB6ZLmhgRjxYbmZlZz1TjffcU4IWI2FzSMcDFwNGShgPHANsAHwbukrRFRCxr7FWYFaORLfddea6u+o+D/1PwnlIn4MAuwJyImAsg6UZgNOAE3MysPmq5744GzkvLtwA/kKRUfmNEvAU8KWlOOt59DYrdzDqhq5L5ruzj3x27IXVE2RPwAcDTufX5wK4FxWJm1hvUct9dvk9ELJX0EtA3lU9rUXdA/UI1szIqe3Jdi7In4DWRNBYYm1ZflTS7k4fqBzy3/LgXr2hk3c77rq+H6unX2NOvD3r3NW7S6ECK0EX37Ib+nTx18cFddaiy/n077sZy3A2kizsdd6v37LIn4AuAQbn1gansfSJiPDB+RU8mqTkimlb0ON1VT78+6PnX2NOvD3yN3UAt993KPvMlrQysR/YwZsPu2d38d9gqx91YjruxHPd7yj4V/XRgmKShklYle7hnYsExmZn1ZLXcdycCY9LykcDdERGp/BhJq0kaCgwDHmhQ3GZm3UapW8BT38LTgDvJhsOaEBGPFByWmVmP1dp9V9L5QHNETASuAK5JD1kuIUvSSfvdTPbA5lLgVI+AYma9UakTcICImARMatDpVrgbSzfX068Pev419vTrA19j4arddyPinNzym8CnW6l7IXBhXQPMdOvfYRscd2M57sZy3ImybwXNzMzMzKwRyt4H3MzMzMysVJyA10DSKEmzJc2RNK7oeLqCpEGS7pH0qKRHJJ2RyjeUNFnSE+nnBkXHuiIk9ZH0kKRfp/Whku5P7+VN6SGy0pK0vqRbJP1N0mOSdutJ76GkL6S/z79KukHS6mV/DyVNkLRI0l9zZVXfM2W+n671YUk7Fhd5OZT5fi1pnqRZkmZKai46ntZ05G+4O2kl7vMkLUi/85mSDioyxmrK+nndRtzd+neePmcekPSXFPf/pPIu/exxAt4OvTft8oHAcOBYZdMpl91S4EsRMRwYCZyarmscMCUihgFT0nqZnQE8llu/GPhORGwOvEA2ZXaZfQ/4bURsBWxPdq094j2UNAA4HWiKiG3JHvirTGte5vfwSmBUi7LW3rMDyUYKGUY2LvblDYqxlHrI/XqviBjRzYdqu5La/4a7kyv5YNyQ3U9GpFejninriLJ+XrcWN3Tv3/lbwN4RsT0wAhglaSRd/NnjBLx9y6ddjoi3gcq0y6UWEQsj4sG0/ApZ4jaA7NquSrtdBRxWSIBdQNJA4JPAz9K6gL3JpsaG8l/fesAeZCNOEBFvR8SL9KD3kOxB8TWUjSW9JrCQkr+HEXEv2cggea29Z6OBqyMzDVhf0sYNCbSceuT9urvp4N9wt9FK3N1eWT+v24i7W0v321fT6irpFXTxZ48T8PZVm3a52/8BdYSkIcAOwP3ARhGxMG16BtioqLi6wHeBLwPvpvW+wIsRsTStl/29HAosBn6eutn8TNJa9JD3MCIWAN8C/kGWeL8EzKBnvYcVrb1nPf7+08XK/vsK4HeSZiibDbRMynzfOS118ZrQ3bpxtFTWz+sWcUM3/50r6746E1gETAb+Thd/9jgB7+UkrQ3cCpwZES/nt6WJM0o5TI6kg4FFETGj6FjqaGVgR+DyiNgBeI0WX0GW/D3cgKyFZyjwYWAtqn993KOU+T2zFfaxiNiRrAvNqZL2KDqgzijZ3/DlwGZkXQ0WApcWGk0byvp5XSXubv87j4hlETGCbLbeXYCtuvocTsDbV9PUyWUkaRWyfxTXRcQvU/Gzla+4089FRcW3gnYHDpU0j+xr6L3J+kuvn7ozQPnfy/nA/IiotCjcQpaQ95T3cF/gyYhYHBHvAL8ke1970ntY0dp71mPvP3VS6t9X+taHiFgE3Eb2wV8WpbzvRMSzKdl6F/gp3fR3XtbP62pxl+V3DpC6dd4D7EYXf/Y4AW9fj5zuPvWHvgJ4LCK+nduUn0J6DHB7o2PrChFxdkQMjIghZO/Z3RFxPNk/pCPTbqW9PoCIeAZ4WtKWqWgfshkGe8R7SNb1ZKSkNdPfa+X6esx7mNPaezYRODGNhjISeCn3lbN9UGnv15LWkrROZRnYH/hr27W6lVLed1o8U3E43fB3XtbP69bi7u6/c0n9Ja2fltcA9iPrv96lnz2eiKcGaYic7/LetMuNmMWtriR9DPgDMIv3+kh/hax/1s3AYOAp4KiIKN1DK3mS9gTOioiDJW1K1iK+IfAQcEJEvFVgeCtE0giyh0xXBeYCJ5P9x7pHvIdp+KejyZ6mfwj4DFm/u9K+h5JuAPYE+gHPAucCv6LKe5Y+wH5A1vXmdeDkiOi2w9N1B2W9X6d7021pdWXg+u4ae0f+hgsKsapW4t6TrCtEAPOAf+9u/8kt6+d1G3EfSzf+nUv6CNlDln1In6cRcX5X5w9OwM3MzMzMGshdUMzMzMzMGsgJuJmZmZlZAzkBNzMzMzNrICfgZmZmZmYN5ATczMzMzKyBnIBbjyPpvyU9kqa5nSlp13b2nyepX1p+tZ19h0iqOmappKmSmjofece0FYuZWS3KeL/s6nufpJMkfTi3/jNJw7vo2IdJOict95d0v6SHJH28K46/AnF9S9LeRcbQ263c/i5m5SFpN+BgYMeIeCt9UKxacFjtkrRyRCwtOg4z6z3Ker/sDEl9ImJZK5tPIpsM5p8AEfGZLjz1l4FD0/I+wKxqx28nvnq4jGwWyrsbeE7LcQu49TQbA89VBsePiOci4p+S9pb0q8pOkvaTdFtrB5G0tqQpkh6UNEvS6NzmlSVdJ+kxSbdIWrNK/f0l3Zfq/0LS2lX2mSrpu5KagTMkHZJrHblL0kZpv/MkTUj7z5V0epVjbZrq7dyRX5aZ9Wplul/uJOkvkv4CnJorP0nSD3Lrv06TryHpVUmXpjq7STpH0nRJf5U0XpkjgSbguvQNwBr51nlJx6Zr+quki3PneVXShSmmaZX7dYuYtwDeiojnlE2adgkwOneeduNLx5kq6TuSmtPvcWdJv5T0hKSv5853gqQH0vF/IqlPel2ZjjlL0hfSe/0U0FfSv7T2vlp9OQG3nuZ3wCBJj0v6kaRPpPJ7gK0k9U/rJwMT2jjOm8DhEbEjsBdwaeVmCGwJ/CgitgZeBj6Xr5hakb4K7JvqNwNfbOU8q0ZEU0RcCvwRGBkRO5DNtvXl3H5bAQcAuwDnSlold74tgVuBkyJiehvXZGaWV6b75c+Bz0fE9h24vrWA+yNi+4j4I/CDiNg5IrYF1gAOjohb0jmPj4gREfFGLrYPAxcDe5PN3LizpMNyx56W4rkX+Lcq598deBAgImYC5wA35c7Tbny5Y70dEU3Aj8mmQD8V2BY4SVJfSVuTzRq8e0SMAJYBx6e4B0TEthGxXfo9VjyYYrQCOAG3HiUiXgV2AsYCi4GbJJ0U2ZSv1wAnSFof2A24o41DCfhfSQ8Dd5FNf15p4Xg6Iv6Ulq8FPtai7khgOPAnSTOBMcAmrZznptzyQOBOSbOA/wS2yW37TUS8FRHPAYtysfQnuxkfHxF/aeN6zMzepyz3yxTD+hFxbyq6psZLXEbWOFGxV/qWcRZZUr1N9WrL7QxMjYjFqYvgdcAeadvbwK/T8gxgSJX6G5P9Xrsivonp5yzgkYhYmL65mAsMIuveshMwPf0e9wE2Tds3lXSZpFFk/wmqWAR8GCuE+4Bbj5P60U0FpqYb2RjgSrL/+f8fWWvNL9rpc308WXK7U0S8I2kesHrlFC1P2WJdwOSIOLaGcF/LLV8GfDsiJqavUM/LbXsrt7yM9/7tvgT8g+xD7dEazmdmtlzJ7pfVLOX9jYmr55bfrPSrlrQ68COgKSKelnRei3076p30HxV4/z057w1gvTaO0ZH4Kp8B7/L+z4N307kFXBURZ7c8iaTtyb5B/SxwFPD/0qbVU4xWALeAW48iaUtJw3JFI4CnACLin2QP2XyV938NV816wKL0YbIX72+RGazs4SWA48i6juRNA3aXtHmKaa3UF7A96wEL0vKYGvaHrBXmcOBEScfVWMfMrDT3y4h4EXhRUqX1/Pjc5nnACEkrSRpE1k2vmkoy+5yyPuZH5ra9AqxTpc4DwCck9ZPUBzgW+H0rx6/mMWDzGvdtK75aTAGOlPQhAEkbStokdfFZKSJuJXsvd8zV2YLs4VMrgFvAradZG7gsfWW5FJhD9vVqxXVA/4h4rJ3jXAf8X2oRagb+lts2GzhV0gSyVufL8xUjYrGkk4AbJK2Wir8KPN7OOc8DfiHpBbIn04e2s3/lfK9JOhiYLOnViJjYbiUzs3LdL08GJkgKsr7rFX8CnkzHfozU57qliHhR0k/JEs5ngPzzMlcCP5b0Bll3m0qdhZLGkfWJF1lXwNvb+D20dC+pP3yutbyqduJrV0Q8KumrwO8krQS8Q9ZP/A3g56kM4GyA9BzR5mTvlxVA7fxNmPUoyp6Wfygirig6FjOz7sz3yxUn6XvA/0XEXUXHkifpcLLhJ79WdCy9lRNw6zUkzSDrc71fZdgtMzP7IN8vu4ay4Ql37W7fTEr6NFnf+xeLjqW3cgJuZmZmZtZAfgjTzMzMzKyBnICbmZmZmTWQE3AzMzMzswZyAm5mZmZm1kBOwM3MzMzMGsgJuJmZmZlZA/1/5v7y+11J8GgAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/tmp/ipykernel_23841/404238302.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_iters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mparams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'betas'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'pi'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresample_hdp_transitions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mstates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mtrans_hypparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mparams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Ab'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Q'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mresample_ar_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mstates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mar_hypparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0mstates\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'z'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresample_stateseqs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mstates\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], "source": [ - "num_iters = 20\n", + "num_iters = 500\n", "plot_iters = 10\n", "keys = jr.split(key,num_iters)\n", "\n", "for i in tqdm.trange(num_iters):\n", " params['betas'],params['pi'] = resample_hdp_transitions(keys[i], **data, **states, **params, **trans_hypparams)\n", " params['Ab'],params['Q']= resample_ar_params(keys[i], **data, **states, **params, **ar_hypparams)\n", - " states['z'] = resample_stateseqs(keys[i], **data, **states, **params)\n", + " states['z'],_ = resample_stateseqs(keys[i], **data, **states, **params)\n", " \n", " if i % plot_iters == 0:\n", " usage,durations = stateseq_stats(states['z'], mask)\n", @@ -174,6 +275,20 @@ "### Gibbs sampling (full model)" ] }, + { + "cell_type": "code", + "execution_count": 10, + "id": "challenging-century", + "metadata": {}, + "outputs": [], + "source": [ + "trans_hypparams = {\n", + " 'gamma': 1e3, \n", + " 'alpha': 100, \n", + " 'kappa': 1e6/50,\n", + " 'num_states':num_states}" + ] + }, { "cell_type": "code", "execution_count": null, @@ -183,7 +298,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "8cb5598062394de8a92d2b924ec66eb8", + "model_id": "261b3292da634bf182e7e61ee81d2473", "version_major": 2, "version_minor": 0 }, @@ -193,18 +308,6 @@ }, "metadata": {}, "output_type": "display_data" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAADQCAYAAABGO9SNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzuElEQVR4nO3deZxcVZ338c+XsIU1EGIGCZhgAoooERoIo4MsAwZFAiMii0PCg2ZmgBFUHifMo5BhGcGRAQFFMxIJiyyCaEaCMQQyqEMgC5GYxECEIB0DCYR9D/yeP+4puGmqu6rTVXWrur/v16tede+559z63arOrZNTZ1FEYGZmZmZmjbFB0QGYmZmZmfUlroCbmZmZmTWQK+BmZmZmZg3kCriZmZmZWQO5Am5mZmZm1kCugJuZmZmZNdCGRQfQaNttt10MHTq06DDMzLpt3rx5T0XEoKLjaCTfs9e1dOlSAHbdddeCIzGzSrq6Z/e5CvjQoUOZO3du0WGYmXWbpMeKjqHRfM9e1wEHHADArFmzCo3DzCrr6p7tLihmZmZmZg3kCriZmVVN0mhJSyUtkzShzPH9Jc2XtFbS0WWObyWpXdIVjYnYzKz59LkuKGZmtn4k9QO+BxwCtANzJE2NiMW5bH8GxgFndnKa84B76hlnb+auJ2a9g1vAzcysWvsAyyLikYh4HbgRGJPPEBHLI+JB4K2OhSXtBQwGft2IYM3MmpUr4GZmVq0dgMdz++0prSJJGwAX03nLeD7veElzJc1dvXr1egXaW33nO9/hO9/5TtFhmFkPuQtKFYZOuP3t7eUXfrrASMzMWtYpwLSIaJfUZcaImARMAmhra4t6BZS/t3em2e75v/zlLwE488yK/48xsybmCriZmVVrBbBjbn9ISqvGfsDfSDoF2ALYWNKLEfGugZxmZr2dK+BmZlatOcAIScPIKt7HAsdXUzAiTihtSxoHtLnybWZ9lfuAm5lZVSJiLXAaMB1YAtwcEYsknSvpCABJe0tqBz4H/FDSouIiNjNrTm4BNzOzqkXENGBah7Szc9tzyLqmdHWOq4Gr6xBer9e/f/+iQzCzGnAF3MzMrEXccccdRYdgZjXgLihmZmZmZg3kCriZmVmLOO+88zjvvPOKDsPMeqhuFXBJu0pakHs8L+kMSdtKmiHp4fS8TcovSZdJWibpQUl75s41NuV/WNLYXPpekhamMpep0uSyZmZmLWzmzJnMnDmz6DDMrIfqVgGPiKURMTIiRgJ7AS8DtwETgJkRMQKYmfYBDgNGpMd44EoASdsC5wD7ki2DfE6p0p7yfClXbnS9rsfMzMzMrBYa1QXlYOBPEfEYMAaYktKnAEem7THANZGZDQyQtD3wSWBGRKyJiGeAGcDodGyriJgdEQFckzuXmZmZmVlTalQF/FjghrQ9OCJWpu0ngMFpewfg8VyZ9pTWVXp7mXQzMzMzs6ZV9wq4pI2BI4CfdjyWWq6jATGMlzRX0tzVq1fX++XMzMzqYuDAgQwcOLDoMMyshxoxD/hhwPyIeDLtPylp+4hYmbqRrErpK4Adc+WGpLQVwAEd0mel9CFl8r9LREwCJgG0tbXVvcJvZmZWD7feemvRIZhZDTSiC8pxvNP9BGAqUJrJZCzwi1z6iWk2lFHAc6mrynTgUEnbpMGXhwLT07HnJY1Ks5+cmDuXmZmZmVlTqmsLuKTNgUOAf8glXwjcLOlk4DHgmJQ+DfgUsIxsxpSTACJijaTzgDkp37kRsSZtn0K2nHF/4I70MDMz65XOOussAL71rW8VHImZ9URdK+AR8RIwsEPa02SzonTMG8CpnZxnMjC5TPpcYPeaBGtmZtbk7r333qJDMLMa8EqYZmZmZmYN5Aq4mZlVTdJoSUvTCsQTyhzfX9J8SWslHZ1LHynpXkmL0mrHn29s5GZmzcMVcDMzq4qkfsD3yGa32g04TtJuHbL9GRgH/KRD+svAiRHxIbJViy+VNKCuAZuZNalGTENoZma9wz7Asoh4BEDSjWSrGC8uZYiI5enYW/mCEfFQbvsvklYBg4Bn6x51LzJkyJDKmcys6bkCbmZm1Sq3MvG+3T2JpH2AjYE/1SiuPuO6664rOgQzqwF3QTEzs4ZJC7BdC5wUEW91kserF5tZr+YKuJmZVauzFYurImkr4Hbg/0XE7M7yRcSkiGiLiLZBgwatd7C90RlnnMEZZ5xRdBhm1kPugmJmZtWaA4yQNIys4n0scHw1BSVtDNwGXBMRt9QvxN5twYIFRYdgZjXgCriZmVUlItZKOg2YDvQDJkfEIknnAnMjYqqkvckq2tsAn5H0b2nmk2OA/YGBksalU46LiAUNv5BuGDrh9op5ll/46QZEYma9iSvgZmZWtYiYBkzrkHZ2bnsOWdeUjuWuAzyC0MwM9wE3MzMzM2sot4CbmZm1iF122aXoEMysBuraAi5pgKRbJP1R0hJJ+0naVtIMSQ+n521SXkm6LC1v/KCkPXPnGZvyPyxpbC59L0kLU5nLJKme12Nm1ltI+nDRMVj3TZo0iUmTJhUdhpn1UL27oHwX+FVEfADYA1gCTABmRsQIYGbah2xp4xHpMR64EkDStsA5ZIs97AOcU6q0pzxfypUbXefrMTPrLb4v6X5Jp0jauuhgzMz6krpVwNMNfX/gKoCIeD0iniVbtnhKyjYFODJtjyGbnirS/LAD0oINnwRmRMSaiHgGmAGMTse2iojZERHANblz1c3QCbdXNSrezKyZRcTfACeQzes9T9JPJB1ScFhWwfjx4xk/fnzRYZhZD9WzD/gwYDXwY0l7APOA04HBEbEy5XkCGJy2yy1xvEOF9PYy6e8iaTxZqzo77bTT+l+RmVkvEhEPS/oGMBe4DPho6sr3rxHxs2Kjs3IeeuihokMwsxqoZxeUDYE9gSsj4qPAS7zT3QSA1HIddYyh9DpeVc3MLEfSRyRdQtY18CDgMxHxwbR9SaHBmZn1cvWsgLcD7RFxX9q/haxC/mTqPkJ6XpWOd7bEcVfpQ8qkm5lZZZcD84E9IuLUiJgPEBF/Ab5RaGRmZr1c3SrgEfEE8LikXVPSwcBiYCpQmslkLPCLtD0VODHNhjIKeC51VZkOHCppmzT48lBgejr2vKRR6SfTE3PnMjOzrt0WEddGxCulBEmnA0TEtcWFZWbW+9V7HvB/Bq6XtDHwCHASWaX/ZkknA4+RLU8M2cpqnwKWAS+nvETEGknnAXNSvnMjYk3aPgW4GugP3JEeZmZW2YnApR3SxpHNXmVNauTIkUWHYGY1UNcKeEQsANrKHDq4TN4ATu3kPJOByWXS5wK79yxKM7O+Q9JxwPHAMElTc4e2BNaUL2XN4tJLLy06BDOrAa+EaWbWt/wvsBLYDrg4l/4C8GAhEbW4aqemXX7hp+sciZm1iqor4JI2i4iX6xmMmZnVV0Q8Rtb9b7+iY7Hu+8IXvgDAddddV3AkZtYTFQdhSvprSYuBP6b9PSR9v+6RmZlZzUn6bXp+QdLzuccLkp6vovxoSUslLZM0oczx/SXNl7RW0tEdjo2V9HB6jO1Y1iprb2+nvb29ckYza2rVtIBfQrYa5VSAiPi9pP3rGpWZmdVFRHw8PW/Z3bKS+gHfAw4hm2p2jqSpEbE4l+3PZIM5z+xQdlvgHLJxQUG2+ubUtMKxmVmfUtU0hBHxeIekN+sQi5mZNYik90vaJG0fIOnLkgZUKLYPsCwiHomI14EbgTH5DBGxPCIeBN7qUPaTwIyIWJMq3TOA0bW4FjOzVlNNBfxxSX8NhKSNJJ1JtnKamZm1rluBNyUNByaRLXj2kwpldgDyDTLtKa0aPSlrZtarVNMF5R/J5oXdgWylyV/TyXSBZmbWMt6KiLWSjgIuj4jLJT1QdFAAksYD4wF22mmngqNpLvvt57GzZr1BxQp4RDwFnNCAWMzMrHHeSHOCjwU+k9I2qlBmBVlLecmQlFaNFcABHcrOKpcxIiaRtcrT1tYWVZ6/T/jWt75VdAhmVgMVK+CSLiuT/BwwNyK89LuZWWs6iewXzgsi4lFJw4BKS9DPAUakvCuAY8kW9anGdODfJW2T9g8Fzup+2GZmra+aPuCbAiOBh9PjI2QtFydLurRukZmZWd1ExOKI+HJE3JD2H42IiyqUWQucRlaZXgLcHBGLJJ0r6QgASXtLagc+B/xQ0qJUdg1wHlklfg5wbkqzbvjsZz/LZz/72aLDMLMeqqYP+EeAj0XEmwCSrgR+A3wcWFjH2MzMrE4kfQyYCLyP7LtAQETEzl2Vi4hpwLQOaWfntueQNdKUKzsZmNyjwPu4p59+uugQzKwGqqmAbwNsQdbtBGBzYNuIeFPSa3WLzMzM6ukq4CvAPDy1rJlZQ1XTBeXbwAJJP5Z0NfAA8B+SNgfu7KqgpOWSFkpaIGluSttW0oy0EtqMUn9AZS5Lq6s9KGnP3HnKrp4maa90/mWprLr/FpiZ9UnPRcQdEbEqIp4uPYoOysysL6hYAY+Iq4C/Bn4O3AZ8PCJ+FBEvRcT/reI1DoyIkRHRlvYnADMjYgQwM+0DHAaMSI/xwJWwzupp+5ItAnFObhDPlcCXcuW8qIOZWXXulvQfkvaTtGfpUXRQZmZ9QTVdUABeBVaSDcgcLml4RNyznq85hnemoppCNg3Vv6T0ayIigNmSBkjaPuWdURqsI2kGMFrSLGCriJid0q8BjgTuWM+4zMz6kn3Tc1suLYCDCojFqnTwwQcXHYKZ1UA10xB+ETidbFDNAmAUcC/V3aQD+LWkAH6Y5nYdHBEr0/EngMFpu7NV0rpKby+TXu4avKiDmVlORBxYdAzWfd/85jeLDsHMaqCaFvDTgb2B2RFxoKQPAP9e5fk/HhErJL0HmCHpj/mDERGpcl5XXtTBzGxdkgaT3cvfGxGHSdoN2C91O7Q6GDrh9op5ll/46QZEYmZFq2YQ5qsR8SqApE0i4o/ArtWcPCJWpOdVZP3H9wGeTF1LSM+rUvbOVljrKn1ImXQzM6vsarL5vN+b9h8CzigqGKvOYYcdxmGHHVZ0GGbWQ9VUwNslDSAbhDlD0i+AxyoVkrS5pC1L22Srnv0BmEq29DHpubSa5lTgxDQbyiiyEforyb4gDpW0TRp8eSgwPR17XtKoNPvJiblzmZlZ17aLiJuBt+DtRXY8HWGTe+WVV3jllVeKDsPMeqhiF5SIOCptTpR0N7A18Ksqzj0YuC3NDLgh8JOI+JWkOcDNkk4mq8gfk/JPAz4FLANeJlsmmYhYI6m0ehqsu3raKWStOP3JBl96AKaZWXVekjSQbKwOpYaPYkMyM+sbqhmE+X6gPSJeI1spbSiwGfB6V+Ui4hFgjzLpTwPvGsadZj85tZNzlV09LSLmArtXugYzM3uXr5L98vh+Sb8DBgFHFxuSmVnfUE0XlFuBNyUNJxvIuCPwk7pGZWZmdRUR84FPkK3z8A/AhyLiwWKjMjPrG6qZBeWtiFgr6Sjg8oi4XNID9Q7MzMxqT9LfdXJoF0lExM8aGpB1y+GHH150CGZWA9VUwN+QdBzZgMnPpLSN6heSmZnVUek+/h6y1u+70v6BwP8CroA3sTPPPLPoEMysBqqpgJ8E/CNwQUQ8KmkYcG19wzIzs3qIiJMAJP0a2K20MFqaFvbqAkMzM+szKvYBj4jFEfHliLgh7T8aERfVPzQzM6ujHXOrEgM8CVRcKljSaElLJS2TNKHM8U0k3ZSO3ydpaErfSNIUSQslLZF0Vs2upA854IADOOCAA4oOw8x6qJpZUB4lTVOVFxE71yUiMzNrhJmSpgM3pP3PA3d2VUBSP+B7wCFAOzBH0tSIWJzLdjLwTEQMl3QscFE69+eATSLiw5I2AxZLuiEiltf0qszMWkA1XVDactubkt1Et61POGZm1ggRcVoaXL9/SpoUEbdVKLYPsCxNM4ukG4ExQL4CPgaYmLZvAa5Ii6UFsLmkDcnWbngdeL4W12Jm1mqqWYjn6Q5Jl0qaB5xdn5DMzKwRUoW7UqU7bwfg8dx+O7BvZ3nSDFrPAQPJKuNjgJVka0l8Jbeo2jokjQfGA+y0U8VeMWZmLaeaLih75nY3IGsRr6bl3MzMrGQfsqXu3wtsA/xG0p2l1vS8iJhEtu4EbW1t7+oCaWbW6qqpSF+c214LPMo7y8ebmVnfsYJsMbaSISmtXJ721N1ka+Bp4HjgVxHxBrAqrb7ZBryrAm6dO+YYf/2a9QbVdEE5sBGBmJlZ40j6DHB7RLzVjWJzgBFpOtoVwLFkFeu8qWTrRtxLtrT9XRERkv4MHARcK2lzYBRwac+uovcZOuH2Cjnex/ILP92QWMysfqpZir5HJPWT9ICkX6b9YWlqqmVpqqqNU3rZqavSsbNS+lJJn8yldzkdlpmZderzwMOSvi3pA9UUiIi1wGnAdGAJcHNELJJ0rqQjUrargIGSlgFfBUr35u8BW0haRFaR/3FEPFjD6+kT3nrjVV5++eWiwzCzHmpEX+7TyW7UW6X9i4BLIuJGST8gm7LqSjqZukrSbmStLB8i6zt4p6Rd0rkqTYdlZmZlRMQXJG0FHAdcLSmAHwM3RMQLXZSbBkzrkHZ2bvtVstmyOpZ7sVy6dc+qn07kU/MuZ9asWUWHYmY9UNcWcElDgE8DP0r7IvsJ8paUZQpwZNoek/ZJxw9O+ccAN0bEaxHxKLCMbDDP29NhRcTrQGk6LDMzq0JEPE92v70R2B44Cpgv6Z8LDczMrJerWAGXtJmkb0r6r7Q/QtLhVZ7/UuDrQKmP4UDg2fQzJmQt1zuk7XWmrgJKU1eVm/Zqhy7SzcysAkljJN0GzAI2AvaJiMOAPYCvFRmbmVlvV00L+I+B14D90v4K4PxKhVIlfVVEzFv/8GpD0nhJcyXNXb16ddHhmJk1g78j6w744Yj4j4hYBRARL5N1CTQzszqppgL+/oj4NvAGvH1zVhXlPgYcIWk52c+bBwHfBQakqalg3Sms3p7eqsPUVZ1Ne1XNdFikmCdFRFtEtA0aNKiK0M3Mer0nIuKefIKkiwAiYmYxIZmZ9Q3VVMBfl9SfbBlhJL2frEW8SxFxVkQMiYihZIMo74qIE4C7yaamgmyqql+k7dLUVZCbuiqlH5tmSRkGjADuJzcdVppJ5diU18zMKjukTNphDY/CumWLD/8t48aNKzoMM+uhamZBOQf4FbCjpOvJWrbH9eA1/wW4UdL5wANkU1aRnq9NU1etIatQk6a4uhlYTLYQ0KkR8SaApNJ0WP2AyRGxqAdxmZn1epL+CTgFeL+k/DSAWwK/KyYqq1ZWAfc84GatrpqFeGZImk+2aIKA0yPiqe68SETMIhvoQ1p2eJ8yecpOXZWOXQBcUCb9XdNhmZlZl34C3AF8i3fm6AZ4ISLWFBOSVevNl5/jqaeeYrvttis6FDPrgU67oEjas/QA3gesBP4C7JTSzMys9URELAdOBV7IPZC0bYFxWRVW//xbHH300ZUzmllT66oF/OIujgXZoEozM2stPwEOB+aR3cvzg+oD2LmIoMzM+pJOK+ARcWAjAzEzs/qLiMPT87CiYzEz66sq9gGXtCnZgJ2Pk7WO/Ab4QeqzbWZmLaRSF8KImN+oWGz9zH7kaYZOuL3LPMsv9EBNs2ZWzSwo15D1D7w87R8PXEsnAybNzKypuXuhmVnBqqmA7x4Ru+X275a0uF4BmZlZ/bh7YWvb8qOfKjoEM6uBairg8yWNiojZAJL2BebWNywzM6sHSQdFxF2S/q7c8Yj4WYXyo8lWNe4H/CgiLuxwfBOyX073IlvN+PNp1hUkfQT4IbAV8Bawt7szds/mH9y/6BDMrAY6rYBLWkj2c+RGwP9K+nPafx/wx8aEZ2ZmNfYJ4C7gM2WOBdBpBVxSP+B7ZKtotgNzJE2NiPyvoicDz0TEcEnHAhcBn5e0IXAd8PcR8XtJA4E3anJFfcja51cDsOFWgwqOxMx6oqsW8MMbFoWZmTVERJyTnk9aj+L7AMvSgmpIuhEYQ7ZScckYYGLavgW4QpKAQ4EHI+L36fWfXq8L6OOe+mXWhf+vjr+wQk4za2adLsQTEY/lH8ArZK0jpYeZmbUoSQMlXSZpvqR5kr6bWqW7sgPweG6/PaWVzRMRa4HngIHALkBImp5e8+u1uRIzs9bTaQW8RNIRkh4GHgX+B1hOtoyxmZm1rhuB1cBngaPT9k11fL0NyaazPSE9HyXp4HIZJY2XNFfS3NWrV9cxJDOzYlSsgAPnAaOAh9LCDQcDs+salZmZ1dv2EXFeRDyaHucDgyuUWQHsmNsfktLK5kn9vrcmG4zZDtwTEU9FxMvANKDsnOQRMSki2iKibdAg93U2s96nmgr4G6mv3gaSNoiIu4G2SoUkbSrpfkm/l7RI0r+l9GGS7pO0TNJNkjZO6Zuk/WXp+NDcuc5K6UslfTKXPjqlLZM0obsXb2bWh/1a0rGSNkiPY4DpFcrMAUak+/jGwLHA1A55pgJj0/bRwF0REencH5a0WaqYf4J1+46bmfUZ1UxD+KykLYB7gOslrQJeqqLca8BBEfGipI2A30q6A/gqcElE3CjpB2Qj5q+k85Hzu5Hd5D8EvBe4U9Iu6TUqjcY3M7McSS+QjeMRcAbZzCSQNci8CJzZWdmIWCvpNLLKdD9gckQsknQuMDcipgJXAddKWgasIbt/ExHPSPpPskp8ANMiouvlHO1dttrnqKJDMLMaqKYCPgZ4FfgKWd+9rYFzKxVKLR4vpt2N0qO0ytrxKX0K2Wj5K+l85PwY4MaIeA14NN3U90n5Ko3GNzOznIjYsoflp5F1H8mnnZ3bfpVOVkqOiOt4p8Jv62Gz4fsWHYKZ1UDFCnhE5Fu7p3Tn5GnO2HnAcLLW6j8Bz6aR8bDuCPp1Rs5LKo2c34F1+5zny3QcjV/2ziRpPDAeYKeddurOJZiZ9VqStgFGAJuW0iLinuIiskreeLodgI0GDik4EjPria4W4in9TPmuQ2QN3FtVOnlEvAmMlDQAuA34wHrG2SMRMQmYBNDW1uYpFM2sz5P0ReB0soGUC8gG299L9iulNamnp18BeB5ws1bX1TzgW0bEVmUeW1ZT+e5wrmeBu4H9gAFpAA6sO4K+s5HznY26r2Y0vpmZlXc6sDfwWEQcCHwUeLbQiMzM+ohOK+CStu3qUenEkgallm8k9ScbLLmErCJ+dMo2FvhF2u5s5PxU4Ng0S8owsp9L76e60fhmZlbeq6m/NpI2iYg/ArsWHJOZWZ/QVR/webwzUr6jAHaucO7tgSmpH/gGwM0R8UtJi4EbJZ0PPEA2Yh46Hzm/SNLNZIMr1wKnpq4tlBuNX+mCzcwMgPbUSPJzYIakZ4DHCo3IambohMoTzCy/8NMNiMTMyum0Ap4W3VlvEfEg2U+aHdMf4Z1ZTPLpXY2cvwC4oEz6u0bjm5lZZRFRms9uoqS7ybr9/arAkMzM+oyKs6BIupWsdfpXEfFW/UMyM7NGkLQn2bLwAfwuIl4vOCSrYOu/PrboEMysBqpZCfNKsvm/H5Z0oST3ETQza3GSziabWnYgsB3wY0nfKDYqq6T/0JH0Hzqy6DDMrIeqmQf8TrLVJ7cGjkvbjwP/BVwXEW/UOUYzM6u9E4A9cgMxLySbjvD8IoOyrr3+5CMAbDy40jAsM2tm1bSAI2kgMA74ItnAye8CewIz6haZmZnV01/ILcADbIKncm16a2ZOYs3MSUWHYWY9VE0f8NvIpqa6FvhMRKxMh26SNLeewZmZWW1Jupysz/dzwCJJM9L+IWRTvJqZWZ1VrIADl0XE3eUORERbjeMxM7P6KjWczCNbobhkVuNDsSJVM1UheLpCs3roain6vYHHS5VvSScCnyWbJ3ZiRKxpTIhmZlYrETGltJ0WMdsl7S71mB4zs8boqg/4D4HXASTtD1wIXEP2s6U7oJmZtTBJBwAPA98Dvg88lO71ZmZWZ11VwPvlWrk/D0yKiFsj4pvA8PqHZmZmdXQxcGhEfCIi9gc+CVxSqZCk0ZKWSlomaUKZ45tIuikdv0/S0A7Hd5L0oqQza3UhfcmA/ccyYP+xRYdhZj3UZQVcUqmLysHAXblj1fQdNzOz5rVRRCwt7UTEQ8BGXRWQ1I+sxfwwYDfgOEm7dch2MvBMRAwnq9Bf1OH4fwJ39DD2PmvTIR9k0yEfLDoMM+uhrirSNwD/I+kp4BXgNwCShpN1QzEzs9Y1T9KPgOvS/gm8M0CzM/sAyyLiEQBJNwJjgMW5PGOAiWn7FuAKSYqIkHQk8CjwUk2uoA96tX0JgCvhZi2u0xbwiLgA+BpwNfDxiIhcmX+udGJJO0q6W9JiSYsknZ7St5U0Q9LD6XmblC5Jl6WfLR9MSySXzjU25X9Y0thc+l6SFqYyl0nS+rwJZmZ90D+SVZy/nB6LgX+qUGYH4PHcfntKK5snItaSNdgMlLQF8C/Av1UKTNJ4SXMlzV29enUVl9J3PHvPFJ69Z0rljGbW1LrsShIRs8ukPVTludcCX4uI+ZK2JGttmUG2oM/MiLgw9R+cQHZTPgwYkR77AlcC+0raFjgHaCObq3aepKkR8UzK8yXgPmAaMBr/tGlm1qXUleT3EfEBsi4hjTARuCQiXqzUVhIRk0iD/dva2qLLzGZmLaiqlTDXR0SsjIj5afsFYAlZy8gYoPTf9ynAkWl7DHBNZGYDAyRtTzYwaEZErEmV7hnA6HRsq4iYnVrnr8mdy8zMOhERbwJLJe3UzaIrgB1z+0N49+qZb+dJ44i2Bp4ma1j5tqTlwBnAv0o6rdvBm5n1Ag0ZTJlGwX+UrKV6cG41zSeAwWm7s582u0pvL5NuZmaVbUO2Eub95PpkR8QRXZSZA4yQNIyson0scHyHPFOBscC9wNHAXamR5G9KGSRNBF6MiCtqcB1mZi2n7hXw1O/vVuCMiHg+/9NjGpRT958XJY0HxgPstFN3G3zMzHqlb3a3QESsTa3W04F+wOSIWCTpXGBuREwFrgKulbQMWENWSbcWVs2KmV4t06x76loBl7QRWeX7+oj4WUp+UtL2EbEydSNZldI7+2lzBXBAh/RZKX1Imfzv4v6EZmYZSZuSDcAcDiwErkqDJasSEdPIxtzk087Obb8KfK7COSZ2I2TL2fbg8UWHYGY1ULc+4GlGkquAJRGRH+RT+nmS9PyLXPqJaTaUUcBzqavKdOBQSdukGVMOBaanY89LGpVe68TcuczMrLwpZIPaF5INfr+42HCsOzYevDMbD9656DDMrIfq2QL+MeDvgYWSFqS0fyVb0v5mSScDjwHHpGPTgE8By4CXgZMAImKNpPPI+h4CnJtbofMUsmkS+5PNfuIZUMzMurZbRHwYQNJVwP0Fx2Pd8MryBQD0Hzqy0DjMrGfqVgGPiN8Cnc01dXCZ/AGc2sm5JgOTy6TPBXbvQZhmZn3NG6WN1Ke7yFism5773xuB5quAu5+4Wfd4SXkzs75lD0nPp20B/dO+yNpCtiouNDOzvsEVcDOzPiQi+hUdg5lZX1e3QZhmZmZmZvZubgE3MzOzunM/cbN3uAJuZmbWIgZ+8rSiQzCzGnAF3MzMrEVsNHBI5Uxm1vTcB7wHhk64vaqf1MzMzGrh5WX38fKy+4oOw8x6yC3gZmZmLeL5+28DYLPh+xYciZn1hCvgZmZm1hSq/VXZgzWt1bkLipmZmZlZA7kCbmZmVZM0WtJSScskTShzfBNJN6Xj90kamtIPkTRP0sL0fFDDgzczaxKugJuZWVUk9QO+BxwG7AYcJ2m3DtlOBp6JiOHAJcBFKf0p4DMR8WFgLHBtY6I2M2s+desDLmkycDiwKiJ2T2nbAjcBQ4HlwDER8YwkAd8FPgW8DIyLiPmpzFjgG+m050fElJS+F3A10B+YBpweEVGv6zEzM/YBlkXEIwCSbgTGAItzecYAE9P2LcAVkhQRD+TyLAL6S9okIl6rf9i9x3aHf63oEJqCF/WxVlfPFvCrgdEd0iYAMyNiBDAz7UPWmjIiPcYDV8LbFfZzgH3JbvznSNomlbkS+FKuXMfXMjOz2toBeDy3357SyuaJiLXAc8DADnk+C8zvrPItabykuZLmrl69uiaB9xYbbjWIDbcaVHQYZtZDdWsBj4h7Sn3/csYAB6TtKcAs4F9S+jWpBXu2pAGStk95Z0TEGgBJM4DRkmYBW0XE7JR+DXAkcEe9rsfMzHpO0ofIuqUc2lmeiJgETAJoa2vzL5s5Ly25B4DNP7h/wZE0P7eSWzNrdB/wwRGxMm0/AQxO2521qnSV3l4mvSy3ppiZ1cQKYMfc/pCUVjaPpA2BrYGn0/4Q4DbgxIj4U92j7YVeeGAaLzwwregwzKyHCpsHPCJCUkNaNtyaYmZWE3OAEZKGkVW0jwWO75BnKtkgy3uBo4G70v1+AHA7MCEifte4kM0651ZyK0qjW8CfTF1LSM+rUnpnrSpdpQ8pk25mZnWS+nSfBkwHlgA3R8QiSedKOiJluwoYKGkZ8FXeGetzGjAcOFvSgvR4T4MvwcysKTS6BbzUMnJhev5FLv20NKJ+X+C5iFgpaTrw77mBl4cCZ0XEGknPSxoF3AecCFzeyAsxM+uLImIa2cxT+bSzc9uvAp8rU+584Py6B2hm1gLqOQ3hDWSDKLeT1E42m8mFwM2STgYeA45J2aeRTUG4jGwawpMAUkX7PLKfPQHOLQ3IBE7hnWkI78ADMM3MzKzG3E3F6qGes6Ac18mhg8vkDeDUTs4zGZhcJn0usHtPYqylcv9A/Q/SzMxqadCRZxUdgpVRTSW9Wq479A2FDcLsC0r/IP2PyczMaqHfZlsXHYKZ1YCXojczM2sRLy68kxcX3ll0GGbWQ24Bb5B8a7hbxs3MbH2UKt9bfPhvC47E6sV9zvsGt4AXbOiE22vad8zMzMzMmptbwJtEvhLu/9mamZlZZ9xK3vpcATczMzPrZar9db2vV9SL+s+MK+BNyH3EzczMrFm4xb32XAFvcuUGb+b5D97MrO94z+cmFh2C9TIeh1YMV8BbnGdXMTPrOzbYaNOiQzArqxkr8s1cH/IsKL1UfnYVz7RiZtY7vDD/dl6Y7/u5WatzC3gf0lkXFndzMTNrDS/98TcAbLmn78tmlTRz42PLt4BLGi1pqaRlkiYUHU9v1FlrulvWzfqeSvdcSZtIuikdv0/S0Nyxs1L6UkmfbGjgZmZNpKVbwCX1A74HHAK0A3MkTY2IxcVG1rd0p2W9s1b2SuW6k9ct9mb1UeU992TgmYgYLulY4CLg85J2A44FPgS8F7hT0i4R8WZjr8LMrHgtXQEH9gGWRcQjAJJuBMYAroD3cfWu5Hcnr/8jYb1INffcMcDEtH0LcIUkpfQbI+I14FFJy9L57m1Q7GZmTaPVK+A7AI/n9tuBfQuKxaymmu2Xhe7k7S2vUYt4etl/pKq5576dJyLWSnoOGJjSZ3cou0P9QjUza16KiKJjWG+SjgZGR8QX0/7fA/tGxGkd8o0HxqfdXYGl6/mS2wFPrWfZZudra129+fp687VB96/vfRExqF7BVFLNPVfSH1Ke9rT/J7JK+kRgdkRcl9KvAu6IiFvKvE4t7tmt+rfjuBvLcTdWX4u703t2q7eArwB2zO0PSWnriIhJwKSevpikuRHR1tPzNCNfW+vqzdfXm68NWvL6qrnnlvK0S9oQ2Bp4usqyQG3u2S343gKOu9Ecd2M57ne0+iwoc4ARkoZJ2phsgM/UgmMyM+utqrnnTgXGpu2jgbsi+6l1KnBsmiVlGDACuL9BcZuZNZWWbgFP/QtPA6YD/YDJEbGo4LDMzHqlzu65ks4F5kbEVOAq4No0yHINWSWdlO9msgGba4FTPQOKmfVVLV0BB4iIacC0Br1cj7uxNDFfW+vqzdfXm68NWvD6yt1zI+Ls3ParwOc6KXsBcEFdA3xHy723ieNuLMfdWI47aelBmGZmZmZmrabV+4CbmZmZmbUUV8Cr0JuWu5e0o6S7JS2WtEjS6Sl9W0kzJD2cnrcpOtaekNRP0gOSfpn2h6VlsZelZbI3LjrG9SFpgKRbJP1R0hJJ+/Wmz07SV9Lf5R8k3SBp01b97CRNlrQqTctXSiv7WSlzWbrGByXtWVzkra1V79eSlktaKGmBpLlFx9OV7vxtN5NO4p4oaUV63xdI+lSRMZbTqt/bXcTd1O95+t65X9LvU9z/ltJr+l3kCngFemfp5cOA3YDjlC2p3KrWAl+LiN2AUcCp6XomADMjYgQwM+23stOBJbn9i4BLImI48AzZctmt6LvAryLiA8AeZNfYKz47STsAXwbaImJ3skF+paXMW/GzuxoY3SGts8/qMLJZQUaQzX99ZYNi7FV6wf36wIgY2QLTtF1N9X/bzeRq3h03ZPeXkenRqDFl3dGq39udxQ3N/Z6/BhwUEXsAI4HRkkZR4+8iV8Are3vp5Yh4HSgtvdySImJlRMxP2y+QVeB2ILumKSnbFODIQgKsAUlDgE8DP0r7Ag4iWxYbWvT6JG0N7E82ywQR8XpEPEsv+uzIBob3VzZ/9GbASlr0s4uIe8hmAcnr7LMaA1wTmdnAAEnbNyTQ3qVX3a+bVTf/tptGJ3E3vVb93u4i7qaW7sMvpt2N0iOo8XeRK+CVlVt6uen/gKohaSjwUeA+YHBErEyHngAGFxVXDVwKfB14K+0PBJ6NiLVpv1U/w2HAauDHqXvNjyRtTi/57CJiBfAd4M9kFe/ngHn0js+upLPPqtfeZxqsld/HAH4taZ6ylUBbTSvfh05LXb8mN1s3jo5a9Xu7Q9zQ5O+5sm6sC4BVwAzgT9T4u8gV8D5K0hbArcAZEfF8/lhaNKMlp8eRdDiwKiLmFR1LHWwI7AlcGREfBV6iw0+OLf7ZbUPWojMMeC+wOeV/Lu4VWvmzsrr4eETsSdZ95lRJ+xcd0Ppqsb/tK4H3k3U1WAlcXGg0XWjV7+0ycTf9ex4Rb0bESLIVe/cBPlDr13AFvLKql09uFZI2IvvHcH1E/CwlP1n6yTs9ryoqvh76GHCEpOVkPz8fRNZvekDq1gCt+xm2A+0RUWpBuIWsQt5bPru/BR6NiNUR8QbwM7LPszd8diWdfVa97j5TkJZ9H9MvQETEKuA2si/9VtKS96GIeDJVtt4C/osmfd9b9Xu7XNyt8p4DpG6edwP7UePvIlfAK+tVy92n/tBXAUsi4j9zh/LLR48FftHo2GohIs6KiCERMZTss7orIk4g+wd0dMrWktcXEU8Aj0vaNSUdTLaqYK/47Mi6noyStFn6Oy1dX8t/djmdfVZTgRPTbCijgOdyPy1b9Vryfi1pc0lblraBQ4E/dF2q6bTkfajDWIujaML3vVW/tzuLu9nfc0mDJA1I2/2BQ8j6r9f0u8gL8VQhTZFzKe8svdyoldxqTtLHgd8AC3mnj/S/kvXLuhnYCXgMOCYiWm6wSp6kA4AzI+JwSTuTtYhvCzwAfCEiXiswvPUiaSTZ4NKNgUeAk8j+I90rPrs03dPnyUbPPwB8kayfXct9dpJuAA4AtgOeBM4Bfk6Zzyp9UV1B1uXmZeCkiGjqqeiaVSver9P96ba0uyHwk2aOuzt/2wWFWFYncR9A1hUigOXAPzTbf35b9Xu7i7iPo4nfc0kfIRtk2Y/0/RoR59a6HuEKuJmZmZlZA7kLipmZmZlZA7kCbmZmZmbWQK6Am5mZmZk1kCvgZmZmZmYN5Aq4mZmZmVkDuQJuvYqk/ydpUVridoGkfSvkXy5pu7T9YoW8QyWVna9U0ixJbesfefd0FYuZWbVa8Z5Z6/ufpHGS3pvb/5Gk3Wp07iMlnZ22B0m6T9IDkv6mFufvQVzfkXRQkTH0dRtWzmLWGiTtBxwO7BkRr6UviY0LDqsiSRtGxNqi4zCzvqVV75nrQ1K/iHizk8PjyBaD+QtARHyxhi/9deCItH0wsLDc+SvEVw+Xk61CeVcDX9Ny3AJuvcn2wFOlifEj4qmI+IukgyT9vJRJ0iGSbuvsJJK2kDRT0nxJCyWNyR3eUNL1kpZIukXSZmXKHyrp3lT+p5K2KJNnlqRLJc0FTpf0mVzLyJ2SBqd8EyVNTvkfkfTlMufaOZXbuztvlpn1ea10z9xL0u8l/R44NZc+TtIVuf1fpkXYkPSipItTmf0knS1pjqQ/SJqkzNFAG3B9+gWgf751XtJx6Zr+IOmi3Ou8KOmCFNPs0j27Q8y7AK9FxFPKFlH7NjAm9zoV40vnmSXpEklz0/u4t6SfSXpY0vm51/uCpPvT+X8oqV96XJ3OuVDSV9Jn/RgwUNJfdfa5Wn25Am69ya+BHSU9JOn7kj6R0u8GPiBpUNo/CZjcxXleBY6KiD2BA4GLSzdCYFfg+xHxQeB54JR8wdSC9A3gb1P5ucBXO3mdjSOiLSIuBn4LjIqIj5KttPX1XL4PAJ8E9gHOkbRR7vV2BW4FxkXEnC6uycyso1a6Z/4Y+OeI2KMb17c5cF9E7BERvwWuiIi9I2J3oD9weETckl7zhIgYGRGv5GJ7L3ARcBDZyo17Szoyd+7ZKZ57gC+Vef2PAfMBImIBcDZwU+51KsaXO9frEdEG/IBsCfRTgd2BcZIGSvog2SrCH4uIkcCbwAkp7h0iYveI+HB6H0vmpxitAK6AW68RES8CewHjgdXATZLGRbbc67XAFyQNAPYD7ujiVAL+XdKDwJ1kS6GXWjcej4jfpe3rgI93KDsK2A34naQFwFjgfZ28zk257SHAdEkLgf8LfCh37PaIeC0ingJW5WIZRHYjPiEift/F9ZiZvUur3DNTDAMi4p6UdG2Vl/gmWQNFyYHpl8aFZJXqD5Uv9ra9gVkRsTp1E7we2D8dex34ZdqeBwwtU357sve1FvFNTc8LgUURsTL9cvEIsCNZ95a9gDnpfTwY2Dkd31nS5ZJGk/0nqGQV8F6sEO4Dbr1K6kM3C5iVbmJjgavJ/tf/32QtNT+t0Of6BLLK7V4R8Yak5cCmpZfo+JId9gXMiIjjqgj3pdz25cB/RsTU9PPpxNyx13Lbb/LOv9vngD+TfaEtruL1zMzW0WL3zHLWsm5j4qa57VdL/aolbQp8H2iLiMclTeyQt7veSP9RgXXvy3mvAFt3cY7uxFf6HniLdb8T3kqvLWBKRJzV8UUk7UH2K+o/AscA/ycd2jTFaAVwC7j1GpJ2lTQilzQSeAwgIv5CNsDmG6z7E1w5WwOr0hfJgazbGrOTsoFLAMeTdR3Jmw18TNLwFNPmqR9gJVsDK9L22CryQ9YCcxRwoqTjqyxjZga0zj0zIp4FnpVUaj0/IXd4OTBS0gaSdiTrqldOqTL7lLI+5kfnjr0AbFmmzP3AJyRtJ6kfcBzwP52cv5wlwPAq83YVXzVmAkdLeg+ApG0lvS918dkgIm4l+yz3zJXZhWzwqRXALeDWm2wBXJ5+rlwLLCP7abXkemBQRCypcJ7rgf9OrUFzgT/mji0FTpU0mazV+cp8wYhYLWkccIOkTVLyN4CHKrzmROCnkp4hG5U+rEL+0uu9JOlwYIakFyNiasVCZmaZVrpnngRMlhRkfddLfgc8ms69hNTnuqOIeFbSf5FVOJ8A8mNmrgZ+IOkVsu42pTIrJU0g6xMvsu6Av+jifejoHlJ/+FxreVkV4qsoIhZL+gbwa0kbAG+Q9RN/BfhxSgM4CyCNJRpO9nlZAVThb8Ks11A2Uv6BiLiq6FjMzJqd75k9J+m7wH9HxJ1Fx5In6Siy6Se/WXQsfZUr4NYnSJpH1uf6kNKUW2ZmVp7vmbWhbHrCfZvt10lJnyPre/9s0bH0Va6Am5mZmZk1kAdhmpmZmZk1kCvgZmZmZmYN5Aq4mZmZmVkDuQJuZmZmZtZAroCbmZmZmTWQK+BmZmZmZg30/wHUQw7A7J9X8gAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" } ], "source": [ @@ -214,13 +317,13 @@ "\n", "for i in tqdm.trange(num_iters):\n", " params['Ab'],params['Q'] = resample_ar_params(keys[i], **data, **states, **params, **ar_hypparams)\n", - " params['betas'],params['pi'] = resample_hdp_transitions(keys[i], **data, **states, **params, **trans_hypparams)\n", - " params['Cd'],params['sigmasq'] = resample_obs_params(keys[i], **data, **states, **params, **obs_hypparams)\n", - " states['z'] = resample_stateseqs(key, **data, **states, **params)\n", - " states['x'] = resample_latents(key, **data, **states, **params)\n", - " states['s'] = resample_scales(key, **data, **states, **params, **obs_hypparams)\n", - " states['h'] = resample_heading(key, **data, **states, **params)\n", + " params['sigmasq'] = resample_obs_variance(keys[i], **data, **states, **params, **obs_hypparams)\n", + " params['betas'],params['pi'] = resample_hdp_transitions(keys[i], **data, **states, **params, **trans_hypparams) \n", + " states['z'] = resample_stateseqs(keys[i], **data, **states, **params)[0]\n", + " states['x'] = resample_latents(keys[i], **data, **states, **params)\n", + " states['h'] = resample_heading(keys[i], **data, **states, **params)\n", " states['v'] = resample_location(key, **data, **states, **params, **translation_hypparams)\n", + " states['s'] = resample_scales(keys[i], **data, **states, **params, **obs_hypparams)\n", " \n", " if i % plot_iters == 0:\n", " usage,durations = stateseq_stats(states['z'], mask)\n", @@ -248,9 +351,9 @@ ], "metadata": { "kernelspec": { - "display_name": "jax", + "display_name": "jax4", "language": "python", - "name": "jax" + "name": "jax4" }, "language_info": { "codemirror_mode": { @@ -262,7 +365,36 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.10" + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false } }, "nbformat": 4,