SNN_Demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### downgrade scipy if smm complained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import sys\n",
    "# !{sys.executable} -m pip install autograd\n",
    "\n",
    "# !{sys.executable} -m pip install seaborn\n",
    "\n",
    "# !{sys.executable} -m pip install ssm\n",
    "# import sys\n",
    "# !{sys.executable} -m pip install scipy==1.2.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import block_diag\n",
    "import autograd.numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import ssm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## all old stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading and dataset functions\n",
    "\n",
    "import pathlib\n",
    "import os.path\n",
    "import csv\n",
    "import numpy as np\n",
    "\n",
    "\"\"\" FILE NAMES \"\"\"\n",
    "# Quality of recordings for a given active neuron\n",
    "f_quality = \"clusters._phy_annotation.npy\"\n",
    "# Which neuron was detected by which channel:\n",
    "f_channels = \"clusters.peakChannel.npy\"\n",
    "# Active neurons indices:\n",
    "f_clusters = \"spikes.clusters.npy\"\n",
    "# Info about channels <-> brain areas:\n",
    "f_areas = \"channels.brainLocation.tsv\"\n",
    "\n",
    "# Time of spike:\n",
    "f_times = \"spikes.times.npy\"\n",
    "# Times of visual stimulations for each trial\n",
    "f_visual = \"trials.visualStim_times.npy\"\n",
    "# Times of auditory cue for each trial\n",
    "f_cue = \"trials.goCue_times.npy\"\n",
    "\n",
    "\n",
    "# Info about the trial time intervals:\n",
    "f_trials = \"trials.intervals.npy\"\n",
    "# Info about the trials used:\n",
    "f_trials_idx = \"trials.included.npy\"\n",
    "\n",
    "# Location of the project\n",
    "project_path = pathlib.Path().absolute()\n",
    "\n",
    "\n",
    "\"\"\" LOADING FUNCTIONS \"\"\"\n",
    "\n",
    "\n",
    "def load_visual_stim_times(recording_name):\n",
    "    # Returns a vector of length n_trials. Each value tells the time at which the stimulus was played for a given trial\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    visual_times_file = os.path.join(data_path, f_visual)\n",
    "    visual_times = np.load(visual_times_file)\n",
    "    visual_times = visual_times.reshape(len(visual_times))\n",
    "    return visual_times\n",
    "\n",
    "\n",
    "def load_cue_times(recording_name):\n",
    "    # Returns a vector of length n_trials. Each value tells the time at which the cue was played for a given trial\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    cue_times = os.path.join(data_path, f_cue)\n",
    "    cue_times = np.load(cue_times)\n",
    "    cue_times = cue_times.reshape(len(cue_times))\n",
    "    return cue_times\n",
    "\n",
    "\n",
    "def load_neuron_regions(recording_name):\n",
    "    # Returns a vector of length n_neurons. Each value tells the brain area of the given neuron\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "\n",
    "    # Load info about neuron's channel\n",
    "    channels_file = os.path.join(data_path, f_channels)\n",
    "    neurons_2_channels = np.load(channels_file)\n",
    "    neurons_2_channels = neurons_2_channels.reshape(len(neurons_2_channels))\n",
    "\n",
    "    # Load info about channel's regions\n",
    "    areas_file = os.path.join(data_path, f_areas)\n",
    "    reader = csv.reader(open(areas_file, \"r\"), delimiter=\"\\t\")\n",
    "    regions_2_channels = np.array(list(reader))[:, 3]\n",
    "    neurons_2_channels = neurons_2_channels.reshape(len(neurons_2_channels))\n",
    "    neurons_2_channels = neurons_2_channels.astype(int)\n",
    "\n",
    "    # Combine the information above to discover which region corresponds to each neuron\n",
    "    neuron_2_regions = regions_2_channels[neurons_2_channels]\n",
    "    return neuron_2_regions\n",
    "\n",
    "\n",
    "def load_neuron_grades(recording_name):\n",
    "    # Returns a vector of length n_neurons. Each value tells the grade (1=bad, 2=decent, 3=good) of the given neuron\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    grades_file = os.path.join(data_path, f_quality)\n",
    "    neuron_2_grades = np.load(grades_file)\n",
    "    neuron_2_grades = neuron_2_grades.reshape(len(neuron_2_grades))\n",
    "    neuron_2_grades = neuron_2_grades.astype(int)\n",
    "    return neuron_2_grades\n",
    "\n",
    "\n",
    "def load_trial_intervals(recording_name, idx=None):\n",
    "    # Returns a vector of length n_neurons. Each value tells the grade (1=bad, 2=decent, 3=good) of the given neuron\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    trials_file = os.path.join(data_path, f_trials)\n",
    "    trials_indices_file = os.path.join(data_path, f_trials_idx)\n",
    "\n",
    "    trials = np.load(trials_file)\n",
    "    trials_indices = np.load(trials_indices_file)\n",
    "    selected_trials = clean_trials(trials, trials_indices)\n",
    "\n",
    "    if idx is not None:\n",
    "        selected_trials = selected_trials[idx]\n",
    "\n",
    "    return selected_trials\n",
    "\n",
    "\n",
    "def load_spike_times(recording_name):\n",
    "    # Returns a vector of length n_spikes. Each value tells the time at which the spike was recorded (in seconds)\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    spikes_file = os.path.join(data_path, f_times)\n",
    "    spikes = np.load(spikes_file)\n",
    "    spikes = spikes.reshape(len(spikes))\n",
    "    return spikes\n",
    "\n",
    "\n",
    "def load_spikes_2_neurons(recording_name):\n",
    "    # Returns a vector of length n_spikes. Each value is the identifier of the neuron to which the spike belongs\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    neurons_file = os.path.join(data_path, f_clusters)\n",
    "    spikes_2_neurons = np.load(neurons_file)\n",
    "    spikes_2_neurons = spikes_2_neurons.reshape(len(spikes_2_neurons))\n",
    "    spikes_2_neurons = spikes_2_neurons.astype(int)\n",
    "    return spikes_2_neurons\n",
    "\n",
    "\n",
    "\"\"\" DATA-SET FUNCTIONS \"\"\"\n",
    "\n",
    "\n",
    "def clean_trials(trial_intervals, included_trials):\n",
    "    #   Remove trials that were not included\n",
    "    #   INPUT: trial_intervals  - onset and offset timestamps of a N trials\n",
    "    #          included_trials - boolean +1 included, 0 not included\n",
    "    #\n",
    "    idx = np.where(included_trials == 1)\n",
    "    idx = idx[0]\n",
    "    return trial_intervals[idx, :]\n",
    "\n",
    "\n",
    "def select_neurons_by_region(recording_name, region):\n",
    "    # Return a vector with the identifiers of all the neurons belonging to a given brain area\n",
    "    neuron_2_regions = load_neuron_regions(recording_name)\n",
    "    neuron_indices = np.where(np.char.equal(region, neuron_2_regions))[0]\n",
    "    return neuron_indices\n",
    "\n",
    "\n",
    "def select_neurons_by_grade(recording_name, minimum_grade):\n",
    "    # Return a vector with the identifiers of all the neurons with score >= minimum_grade\n",
    "    neuron_2_grades = load_neuron_grades(recording_name)\n",
    "    neuron_indices = np.where(neuron_2_grades > minimum_grade)[0]\n",
    "    return neuron_indices\n",
    "\n",
    "\n",
    "def generate_spike_counts(recording_name, brain_region, minimum_grade, dt, t0, tf, firing_rates=False):\n",
    "    # Selects the neuron belonging to a given brain region, and generates a data-set.\n",
    "    # Returns: a matrix [n_neurons * n_time_steps] describing how many times each neuron spiked at a given time step\n",
    "    # Arguments:\n",
    "    # brain_region: the brain region we want to analyze\n",
    "    # minimum_grade: the minimum score we accept to include a neuron in the data-set (1=bad, 2=decent, 3=good)\n",
    "    # t0 = the initial time step (seconds)\n",
    "    # tf = the final time step (seconds)\n",
    "    # dt = the length of each time step (seconds)\n",
    "    # firing_rates = if false, the returned values are spike counts. if true, it is firing rates\n",
    "\n",
    "    # select the neurons based on brain region and score\n",
    "    neuron_indices_1 = select_neurons_by_region(recording_name, brain_region)\n",
    "    neuron_indices_2 = select_neurons_by_grade(recording_name, minimum_grade)\n",
    "    neuron_indices = np.intersect1d(neuron_indices_1, neuron_indices_2)\n",
    "\n",
    "    # load spikes and spike mappings\n",
    "    spikes = load_spike_times(recording_name)\n",
    "    spikes_2_neurons = load_spikes_2_neurons(recording_name)\n",
    "\n",
    "    # generate the time steps sequence\n",
    "    time_steps = np.arange(t0, tf, dt)\n",
    "    n_time_bins = len(time_steps) - 1\n",
    "    time_bins = time_steps[1:]\n",
    "\n",
    "    n_neurons = len(neuron_indices)\n",
    "    spike_counts = np.zeros((n_neurons, n_time_bins))\n",
    "\n",
    "    for i, neuron_id in enumerate(neuron_indices):\n",
    "        # for each neuron, retrieve its spikes\n",
    "        neuron_spikes = spikes[spikes_2_neurons == neuron_id]\n",
    "\n",
    "        # generate an histogram of spike counts\n",
    "        [neuron_spike_counts, _] = np.histogram(neuron_spikes, time_steps)\n",
    "        if firing_rates:\n",
    "            neuron_spike_counts = neuron_spike_counts / dt\n",
    "        spike_counts[i, :] = neuron_spike_counts\n",
    "\n",
    "    return spike_counts, time_bins\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "# functions for selecting trials of choice\n",
    "\n",
    "import pathlib\n",
    "import os.path\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "\n",
    "\"\"\" FILE NAMES \"\"\"\n",
    "# Interval of trials [onset, offset]\n",
    "f_intervals = \"trials.intervals.npy\"\n",
    "# Which trials were included [boolean]\n",
    "f_included = \"trials.included.npy\"\n",
    "# Timing of visual stimulation onset\n",
    "f_visStim_times = \"trials.visualStim_times.npy\"\n",
    "# Timing of the auditory go cue\n",
    "f_cue_times = \"trials.goCue_times.npy\"\n",
    "# Response choice: -1 for right, +1 for left, 0 for no-go\n",
    "f_choice = \"trials.response_choice.npy\"\n",
    "# Response time\n",
    "f_response_times = \"trials.response_times.npy\"\n",
    "# Feedback time \n",
    "f_feedback_times = \"trials.feedback_times.npy\"\n",
    "# Feedback type:  -1 for white noise (miss), +1 for water reward (hit)\n",
    "f_feedback_types = \"trials.feedbackType.npy\"\n",
    "\n",
    "\n",
    "# Location of the project\n",
    "project_path = pathlib.Path().absolute()\n",
    "\n",
    "\n",
    "filenames = list([f_intervals, f_included, f_visStim_times, f_cue_times, f_choice, f_response_times, f_feedback_times, f_feedback_types])\n",
    "\n",
    "\n",
    "def load_trial_files (recording_name):\n",
    "    data_path = os.path.join(project_path, 'data', recording_name)\n",
    "    alldata = []\n",
    "    for filename in filenames:\n",
    "        file = os.path.join(data_path, filename)\n",
    "        data = np.load(file)\n",
    "        alldata.append(data)\n",
    "        \n",
    "    trials = {\n",
    "        'interval': alldata[0],\n",
    "        'included': alldata[1],\n",
    "        'visStim_times': alldata[2],\n",
    "        'cue_times': alldata[3],\n",
    "        'choice': alldata[4],\n",
    "        'response_time': alldata[5],\n",
    "        'feedback_times': alldata[6],\n",
    "        'feedback_types': alldata[7],\n",
    "    }\n",
    "    \n",
    "    return trials\n",
    "    \n",
    "\n",
    "def extract_clean_trials (recording_name):\n",
    "#   Remove trials that were not included\n",
    "#   INPUT: trial_intervals  - onset and offset timestamps of a N trials\n",
    "#          included_trials - boolean +1 included, 0 not included\n",
    "  trials = load_trial_files (recording_name)\n",
    "    \n",
    "  included_trials = trials['included']\n",
    "  idx = np.where(included_trials == 1 )\n",
    "  idx = idx[0]\n",
    "  \n",
    "  # Undate trials [dictionary]\n",
    "  trials['interval'] = trials['interval'][idx,:]\n",
    "  trials['included'] = trials['included'][idx,:]\n",
    "  trials['visStim_times'] = trials['visStim_times'][idx,:]\n",
    "  trials['cue_times'] = trials['cue_times'][idx,:]\n",
    "  trials['choice'] = trials['choice'][idx,:]\n",
    "  trials['response_time'] = trials['response_time'][idx,:]\n",
    "  trials['feedback_times'] = trials['feedback_times'][idx,:]\n",
    "  trials['feedback_types'] = trials['feedback_types'][idx,:]\n",
    "      \n",
    "  \n",
    "  return trials\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for mapping best colors scheme \n",
    "\n",
    "import numpy as np\n",
    "from sympy.utilities.iterables import multiset_permutations\n",
    "\n",
    "\n",
    "def map_states(n_states, sequence1, sequence2):\n",
    "    # Maps the best permutation of sequence 2 to sequence 1\n",
    "    all_states = range(n_states)\n",
    "\n",
    "    min_T = np.min([len(sequence1), len(sequence2)])\n",
    "    sequence1 = sequence1[:min_T]\n",
    "    sequence2 = sequence2[:min_T]\n",
    "\n",
    "    best_score = -np.inf\n",
    "    best_mapping = sequence2\n",
    "    for permutation in multiset_permutations(all_states):\n",
    "        permutation = np.array(permutation)\n",
    "        remapped_sequence2 = permutation[sequence2]\n",
    "        score = np.sum(remapped_sequence2 == sequence1)\n",
    "        if score > best_score:\n",
    "            best_score = score\n",
    "            best_mapping = remapped_sequence2\n",
    "    return best_mapping\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting functions\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def plot_psths(dataset, time_bins, title, visual_time, cue_time, feedback_time):\n",
    "    plt.imshow(dataset, cmap='binary', aspect='auto')\n",
    "    plt.title(title)\n",
    "    plt.ylabel('neurons')\n",
    "    plt.xlabel('time (ms)')\n",
    "\n",
    "    x_ticks = [0, len(time_bins)-1]\n",
    "    x_labels = (time_bins[x_ticks] - time_bins[0]) * 1000\n",
    "    plt.xticks(x_ticks, x_labels.astype(int))\n",
    "\n",
    "    visual_time_bin = np.where(visual_time < time_bins)[0][0]\n",
    "    plt.axvline(visual_time_bin, color='r')\n",
    "\n",
    "    cue_time_bin = np.where(cue_time < time_bins)[0][0]\n",
    "    plt.axvline(cue_time_bin, color='g')\n",
    "\n",
    "    feedback_time_bin = np.where(feedback_time < time_bins)[0][0]\n",
    "    plt.axvline(feedback_time_bin, color='k')\n",
    "\n",
    "    plt.legend(['visual', 'cue', 'feedback'])\n",
    "\n",
    "\n",
    "def add_states_2_psth(fig, states, colors, n_neurons):\n",
    "    # Create a Rectangle patch\n",
    "    previous_state_change = 0\n",
    "    state_changes = np.where(np.diff(states))[0]\n",
    "    for state_change in state_changes:\n",
    "\n",
    "        state = states[previous_state_change]\n",
    "        c = colors[state]\n",
    "\n",
    "        state_dt = state_change - previous_state_change\n",
    "        rect = patches.Rectangle((previous_state_change, 0), state_dt, n_neurons, edgecolor='none', facecolor=c, alpha=0.2)\n",
    "        fig.add_patch(rect)\n",
    "        previous_state_change = state_change\n",
    "\n",
    "    state = states[previous_state_change]\n",
    "    c = colors[state]\n",
    "\n",
    "    state_dt = len(states) - previous_state_change\n",
    "    rect = patches.Rectangle((previous_state_change, 0), state_dt, n_neurons, edgecolor='none', facecolor=c, alpha=0.2)\n",
    "    fig.add_patch(rect)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "# main\n",
    "\n",
    "# Experimental Parameters\n",
    "recording_name = 'Cori_2016-12-14'\n",
    "brain_region = 'ACA'\n",
    "neuron_min_score = 2\n",
    "\n",
    "# Model Parameters\n",
    "bin_dt = 0.01  # seconds\n",
    "pre_stim_dt = 0.5 # seconds\n",
    "post_resp_dt = 0.5 # seconds\n",
    "n_compoments = 5\n",
    "colors = ['r', 'g', 'b', 'y', 'm', 'c']\n",
    "\n",
    "# Choose trials\n",
    "trials = extract_clean_trials(recording_name)\n",
    "conditioned_trials = np.where(trials['choice'] == 0)[0]\n",
    "\n",
    "# Run the fit\n",
    "# plt.figure()\n",
    "n_trials = 8\n",
    "for i in range(7,n_trials):\n",
    "#     fig = plt.subplot(1, n_trials, i + 1)\n",
    "\n",
    "    # Load time pointers for the given trial\n",
    "    trial = conditioned_trials[i]\n",
    "    visual_time = trials['visStim_times'][trial]\n",
    "    cue_time = trials['cue_times'][trial]\n",
    "    feedback_time = trials['feedback_times'][trial]\n",
    "\n",
    "    # generate the spike count histograms\n",
    "    t0 = visual_time - pre_stim_dt\n",
    "    tf = feedback_time + post_resp_dt\n",
    "    [dataset, time_bins] = generate_spike_counts(recording_name, brain_region, neuron_min_score, bin_dt, t0, tf)\n",
    "    (n_neurons, n_bins) = dataset.shape\n",
    "\n",
    "#     # Create a hmm model\n",
    "#     model = hmm.GaussianHMM(n_components=n_compoments, n_iter=1000)\n",
    "#     model.fit(dataset.T)\n",
    "#     [logprob, states] = model.decode(dataset.T)\n",
    "\n",
    "    # Find the best mapping of the state sequences\n",
    "#     if i == 0:\n",
    "#         states_trial0 = states\n",
    "#     else:\n",
    "#         states = map_states(n_compoments, states_trial0, states)\n",
    "\n",
    "#     # Plot\n",
    "#     title = brain_region + ' trial#' + str(trial)\n",
    "#     plot_psths(dataset, time_bins, title, visual_time, cue_time, feedback_time)\n",
    "#     add_states_2_psth(fig, states, colors, n_neurons)\n",
    "\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## step by step smm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f7580bda320>"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAACYCAYAAAAfpakdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deXwc1ZUv8N+pKi1Wu2VJlizLtqzFlmUsyZI3WRIQIGY1SUwCJJCQAIE4EJzJAgOYNZCXBBIgy/vkJeMkzMvLyyNkSDLJAJOdJJOxAeMAxsaLJEuydsmyrH3r7vP+qJbUra7eu1XdrfP9fPxxq7Z77q2q29W3qu4lZoYQQojEopgdgBBCiMiTyl0IIRKQVO5CCJGApHIXQogEJJW7EEIkIKnchRAiAUWlcieiK4noBBE1ENH90UhDCCGEdxTp59yJSAVwEsBlANoAHARwIzO/G9GEhBBCeBWNK/dqAA3MfIqZJwH8DMCuKKQjhBDCCy0K21wJoNXl7zYA232tkEwpnApLFEIRQojENYT+M8ycYzQvGpV7QIhoN4DdAJCKNGynHWaFIoQQATv3iVqk9tuR+h+vmx0K/sgvtHibF43KvR1Avsvfq5zT3DDzPgD7ACCdsqSDGyFEXMj6+ZtgZsR6pRWNNveDAEqIqIiIkgHcAOA3UUhHCCHmnWN8HDwxYXYYfkX8yp2ZbUS0B8DvAKgAnmXmo5FORwghhHdRaXNn5pcBvByNbQshhPBP3lAVQogEJJW7EEIkIKnchRAiAUnlLoQQCUgq94VCUaGmp5sdhRCm0wryQSkpAS+vZmYCRGGlqVgsnmkSQc1YEtD6lJQMxWoNLs2glhZxSy0pQstd5WaHIYTpdv32EGx1ZQEv3/DP66EuCe/CqPeGjZi8qMJtmpqdjca7NwS0vqN6A7purvC/oIuI9woZinTKYul+QAghgvNHfuEQM281midX7kIIkYASonJvfbgOWnGh2WFAzV2G5q/Umh2GEEIkRuVe8MzbsDW3+l8wyuw9vSh+4ojZYQghgkTbKtBzV53ZYUSUaV3+RpJjZMTsEHTMcAwNmR2FECJI9E498possJsdSAQlxJW7ENHiuGgTJnZuMzuMhNazp86Ux3RHrt0OVOtPoDjGx2E/0+exTNfn60BacNfAXV+oAxTVeGZ1BUY/5HPsIq+0ogKcvTXwZt+EuHIXIlqSm88gKTkppCs6244taPoEI+tvKVCngIyfvIq+T9YAALJ/fBBss0U22DiVc2gEjrHxeU83/Vg/aGAYvvZC7sFRsD24vZ/7+ijgMF5Ha++DdWg8pOOJ+89h6VuLA+5HXip3IXywtYR+Lyfp3DhS65fA0m0D2fVT0tKtn9bsMP8R5FhBB942ZeAL+7sn/S5D//1W0Nv1tY6tvcNg6KLA2M8NAG8OBLy8VO5CRAkfOor8Q+7TUl80f2i2RKSUr4cyMgZbk9dR5xYcaXMXQsS9sdVW2HIDe5V/oZArdyFE3Et5+aDZIcQcuXIXQogEFHLlTkT5RPQKEb1LREeJ6HPO6V8ionYiesv5b2fkwhVCCBGIcJplbADuZuZ/EJEVwCEi+oNz3jeZ+anwwxNCCBGKkCt3Zu4E0On8PERExwCsjFRgQgghQheRNnciKgSwCcBrzkl7iOgwET1LRJmRSEMIIUTgwq7ciWgxgF8A+DwzDwL4HoA1AKqgX9k/7WW93UT0BhG9MYWJcMMQQoigDV+/HaMfDK07gFgX1qOQRJQEvWL/KTP/EgCYudtl/g8AvGi0LjPvA7AP0AfrCCcOIYQIxZLfHQMzw2F2IFEQztMyBOBHAI4x8zMu0/NcFvsgAOkDVwgRc0au2w77+oKE7ck1nCv38wF8HMA7RDTdmcIDAG4koioADKAZwKfDilAIIaIg/ehZ0Lkhnx2HxbNwnpb5OwCjIcFfDj0cIYSYH/Zj9WaHEFXyhqoQQiQgqdyFECIBSeUuhFiQtOW5UDMStyfJuK7c1eylXudRUnJC7zjhydfxEI31AkGaBjUzM+x0ohnjQtV7RTGmKovNDiNq4rpyb7h7HZTUVMN5vKkUHR8vm+eIhFmU1FQ03LMu6PVI01AfwnqBUgrz0XLHeVDS0tBwd+jp1N+zLuixPIVvWf/vEJT/Omx2GFFDzOa/P5ROWbyddkR+w0RADORPCBF7eu+oRcapKST9/g2zQwnZH/mFQ8y81WheXF+5+8K1lej8QuAjhYv51XNXHWzv3WJ2GGIBy/n+gZAq9v5bajF2TXUUIoosuXIX5pB9I2LA4I01UGyMxf/2mv+Fp5Hz9Z4YOH4X5JU7AJ+F/+FjXUB1xTwGI9zEwIkhwkdJyWj6mv4LufkrtVAsFpMjCs6SF/4B67+/GdxKzHFx/C7YOzS/+NCFoFMnEfu7SIjYxVOTWPutRtgBrPlOI+yjo2aHFBSemjQ7hKhZsJV7or96LMR8sXf3uP0fTyau2gbFznF9U9WbuGiW0YoLMXy93ucyJSWj71ML40bp4I010PJXRTWNvttmf0qP7aqGWlYa1fRE7FBL12L0Q559mZ/ZXRuVxy61vOUY+FhNxLcbDsvxHiw6GfyXEtdWwrbD/YEA0rSYqpvi48p9dAyWDn1AD7bbYW1N1H7c3C3umACPjUU1jfTWKbBNL89FnWOgwZGopidiBw2PYlGnZxt5+mkb2BH5Bksen8Di9thqBrE1tYS0ntY3DHUkya0feHZwTNVNif20jBBCJLC4f1pGzcxckM0Fysb1UNPTI7MtqxXKxvUR2VZYcVRtgJKWZnYYoJQU0CZ5g9kfpXx9xLrxUCwWKJXnRWRb4VLLSme6hXAVyeOCtpSBUlIisq1QxEXljuXZ6K3OMjuKeddTkwnKDj/fitUKe0UxemrneaxyIqhri9wmdZ2fASUzY37jMKBYF6PrQul7yJ/emkwgJzL92ihZmeiuNX/fA8CZbVlAXs7M34rVCrWkGNiw1vC40IoKQEnJQaXRdf4SKIvNezRUmmUWANpSht4t6cjed2Be01VSU9F872asfnz/vKYrRLBoSxla3r8EyeeA5d/yPF477q3D6p82w9beYUJ03vlqlpHKXQgh4pSvyj2sp2WIqBnAEAA7ABszbyWiLADPAyiEPobqh5m5P5x0hBAiYmKo+4BoikSb+yXMXOXy7XE/gD8xcwmAPzn/FkKImHD60Vo0PL0dtLXc7FCiKhrPue8CcLHz848B/AXAfVFIRwghgrb6Mf3eUyw0SUdTuFfuDOD3RHSIiHY7p+Uyc6fzcxeA3DDTEBFGmoamJ2LnTTqRWNr21ulPnsSotr210IoKzA4j6sK9cr+AmduJaBmAPxDRcdeZzMxEZPj16Pwy2A0AqTD/ueeFhG02rP2m3tmTEJFWsO84HAODZofhVcG+E7CfGzA7jKgLq3Jn5nbn/z1E9CsA1QC6iSiPmTuJKA+AYccNzLwPwD5Af1omnDhE8OKxkycRH+x9Z80Owaeu69dhajEh/6UzsL970uxwoibkZhkishCRdfozgMsBHAHwGwA3Oxe7GcCvww1SCCEiJXf/Oaz8yxDQ2Wt2KFEVTpt7LoC/E9HbAF4H8BIz/xbAEwAuI6J6AJc6/w6a44IqDN5Yg8Eba2C/ZLPPZceuqYaak+NzGWFs+MM1HoOMK1UbwHWVYW+bkpIxdIPvXgCVjevB51eFlY5iscz0GupVdcW8PB3huGjTguwqI5443j4GPvgO7P36E9pK+XoM3liDySsMHxc3PEfiQcjNMsx8CoBHDcDMfQDCfiOJWP+nb9RfMADY4WchYcjoiQFm/2UeIApkt0TgqQXjOztz5s/H0xGOeUpHRAwxu9c3c8Xp/pQ3VIUQIk7Ffa+QIjDq2qKoDLIg4pscFwuTVO4JpOX6PChLItNFsEgcclwsTNIsI4QQcUqaZYQQYoGRyl0IIRKQVO4LWONTNUGPLrNQdH+2Do4LN5kdhpgHE1dtQ9/t7n0tkaah8Wnf72i4UjMzcSrG+muSNvcFTElNhWN83OwwYhIlJYPtdsAhPfAkPEUFqSp4atJ9cpDnhxnnk7S5C0NSsXvHU5NSsS8UDvtMxa7mLkPb3jp9cpDnR6ydT1K5CyGEk6PvLAqfazM7jIiQyl0IIZzYZoOt+bTZYUSEVO5CCJGApHIXQogEJJW7EEIkIKnchRAiAUnlLoRY0BSLBWrGErPDiLiYqdwViyXyyyoqlDTfg28rqalBd4eqpKYCiup7mUBjJPIbYzhI00ApKRHdZjD7Kpx1ZtZNSwOIIpsOUcDLkqZFbSQeJS3N77EUiXVcUVKy1zeTw9lPbtsJ4ByJdJqhGrtoA3qu3xD0epE6LqKV/5ip3Bse3gg13X+3pEpqKhoe3RjQNtWSIpz+vO8h3Lpu3wx7bUVA25vWcedmYKuPg0FRUf9YRUAVkrZyBZruD39IO2/GL9uEMzf5HqYwWA2PbAz6gKz/0saQv2RO7a2EtiLP/4JEgZf78lycuj+w48h+fgW6botsGU5r21MFZWNww/Kd/lwV1PVrQk5z8NrNGN7l2bUCaZpefhHg9xxx0fjQRqhLsyKSbkgY+ghaQYrEcaGkpaHhkcCOw2CF3P0AEZUCeN5lUjGARwBkAPgUgOnRZx9g5pd9bUu6HxBCiOD56n4gnDFUTwCoAgAiUgG0A/gVgFsBfJOZnwp120IIIcITqWaZHQAambklQtsTEaaWrp3pM0MIAfR9qha2HVswecVWnL21duYcUaxWNH3Ns4fHjnvroJYF1oTW9LXaqLWlq+npaPqq/x4oI1W53wDgOZe/9xDRYSJ6logyA9nA6Uf1QnU1/r5qjFy33W0apaSg5TH/lRRtLUfvHYF1wTl0Qw0mr9wGAFDLStH9T8bbb3ugDmpOjsd0pWoDej7jso6iouXx2pm23/b76qAV5HtNX01Px+lH9PVbHq/1uME7+NEaTF26ZeZvrSAfHfd6xtj6UB3UpVnouasOSuV5AADbji0YuKkGjsZmFHzvqNcY5uK6SvTdZlx+px+t83oT2F9Z+OOtLLq+UAd1nXs7s7JxPXr26MsGclxMl4VaUozOLxov235/HbSVKwKKda7eO2tBm8r8LtfyWF1A9x/8HRfTHBduQv8t4XU3q1gsOP2ocZn03V4LrvW8L6RmLEHrw8brzC2LYI6L4eu3Y+LqbUHmIHg5PzuCsZwkTKaryP75YUAhsAY4hoex9qkTHsvn7zsKx4nGgLa99qkTcIyMRDpkAIB9aAhrn/aMb66wu/wlomQAHQDKmLmbiHIBnIF+m+LLAPKY+ZMG6+0GsBsAUpG25eKiO2BraQVc4lGXZgFEsJ/pc10RWuFq2Jp8/0hQLBYo6VbYOrv85kHNyQFsNtj7+6GkpkLJXgpbW7vHclpBPuztnWCbzT2ttDQomRmwtXfMLltUMBOjVpAPe0e3R5eisxtQoa1eCVvzabf13OKbmoT93IBeBEnJUJcvg63VvYMjrXA1bKfboeXlwtF/Do7RUf0Rr5QU2Lt7/JaDW0hWKxRLGmxd3Z7lUFSg979hcOz4Kwv/CRuXhbZqJRxn+tx63nNLK4DjYrosHOcGoC7L8Sg/IIB95YO2cgUc5wb8ntS+ys+Nn+NimpqeDkpbZLivAuaj/LTluXCMjMIxNOQ1Po915pRFMMeFmr0UsDtg7+8PPT8BUrOXAsyw950FpaR4PS5ila8290hU7rsA3MXMlxvMKwTwIjOX+9pG2rJ8Pv/s+dLFqhAioigpGd2f2opl/2u/2aFERbT7c78RLk0yROT6zNoHARzxtwF1nAEO4VkkETRKScHgRwMfYSastJKSMXBTDQZumk1v5Lrt+tWSydSMJRj6SPDloKSlYfDG2fUGPlYT8fcIxCxt5QqMXVMNJTV15lia/hfIM+ZstyPr+MQ8RDpr4KYaj+Ym+yWboZSvN16ByO0ciZSQn5YBACKyALgMwKddJn+diKqgN8s0z5lnvJ2hUYDMHxFqQXAwtPH5K+u5aanjDNjN/yJnuwPaeAhxOBxueUoaY8Ahx27U2O1QxxnMnsdtQK0ODju0Px+KUnDGjM4vZcIOmrIZLO19nXDJMHtCCBGnZJg9IYRYYKRyF0L4pGYvjWr/R9EW6qOt8U4qdyGET93XroO9qsTsMEKipKai4TMFZodhirBuqAohEl/2vxwwO4SQOcbHUfhg/MYfDrlyF0KIBCSVuwGtcLX+arQQQsQpqdwN2JpPo/DL8/tsrBlOPVnrddCGiKf1RK287CN86rhntk8kMyRa53pSuXsRSt8isWi6sywjJU8cn7d8ljx5HDwxv28KiviSv+8I+Gi9KWnT1nJ0XroMBd91f6G+5fH5uwCKNLmhmuBSDtYjJSUFRr32zEfHTGakJeKTfXDQtLTpWBNWtC2GbU4Ma55thy1OL/Skck9wZp4wQsQLx8iIYW+eRj1exgtplhFCiAQklbsQQiSgmKnc+fwqryPNzDfFYgFqfI9Irm5YB23VyqjGQVvL9QEmQqCWlfqNz3FBVcAjJEUC11XG3s0pIr0cIkhdWwR1bVFgyW8q0weliQVRKAtf5qY132Xheo6oucsi/qQO11YG1C1xtMRM5d6zOQ2UbHziq+npAZ8sM+tkLIFWXGg8b8M6n18kinUxzmz0Pf7h8LoM2FZFvl9y1z6f+yqsoBAr96HSDNjzsqDm5EDLXwUlLc1jmLqeLWkgVYV6XgkoJQVaUQHUzIBGRXSjFeR79NFOSclQyte7jTnZu8kCVK5zO+DVnJyAviTVdWvcxqRUl2ZBK1ytj5xVvh7qhnWz88pK3b605sbnVhakoGeLe78paklx0ONfquvWzPS/MlGQhfHCwCqp/vJ0ICsjqLQiyV9ZRA2RW1pq6Vqc2Rx4WahlpYCihpS0WroWSvl6nN2UiZHKFVDK12O8cjXObM7Qy0JRZ47b6WNp+hwB9FHB1NxlXrc/fQ6fqUoDWbwMR2lwPkZazFTuuf9zPxyjo8Yz8/PQvWN5UNvjwpXouTjPcN7pD2T77AjJ1tWN7H2+X1le9O+vA68eDiqmQLTsypr54sn61wMh39BJ++Vr4IPvYLIsH2cuyYejfA06L891W2b5t/eDbTa0XZUDNSsTve9ZAcfaVUGn1V+7ElPnrXabpqQvRss1WTh56+yXxbLv7kfzB6xQ8mbjmKhYjXO1/tPsvCwXlDd7QtlK89F3/gooOdlo2ZWF0+/PnpnX8oGlUFyeqe+7YCVspbNj2CrLl+llQQQtLxfLv+0+Sk/XpbmgVcbHjjcdV+RCydXH19X+dAhJfwzsPYmMnxyAvf5UUGlF0kxZAIDD7lEWUcPsllb7zmXI/d3pgMui5QNLoSwK7aq4/cplaPhEJjKOD4MVQsuuLIxnalj21050Xp4LxZKGk7foX84nbssEaUkz5wgADGxfhfEK4zGRtfxVaNmVBSgqcr53APa+s4bLKcuy0XFlruG8SImZ/tzrlnzQ65MdpGmglBSPu9lqevrMOorFAp6Y8Bjf1Nvy0+s4xsZDGt5PSUsDT9n8PieupqfDPjTkNl4mpaSAiNzGA/Wajku+KCkZlKS5fQm65mluWkpqqj6gwcZ16N5uxfIfv+MxDqa3dYyeSZ9Oa245zuQrKRnKolSw3Q7HyAiUtDSc2luJNd84OrN8y+O1WPOvHbA1tfhMa27eFKsVjpFRt31ldFx4i80zWIK2PBeNdxah4JEDnsdGWhpI0+CYmAj4+XyvxyARVKvVMK5gj1vFaoVjeNj/+KsG686k6aXcfR1LM1lxPQYjnS+j/ZucbHjRF/B+NlgPANo/WY7c10agvHEMPDXpdiwpFgua7tuIgkcOoOGZGpTc94+Z83zuee+6HiUlo+nRLSh8aPbCcG6+gj2H3SgqFEsaHENDM+vERX/uDQ+UzRQ8FNW9LbiiFG13eo6+Xv9g2cxP/PY7KoHKUo9lppGm4eTD7iPTt+6phLIhtN7uej9aifFLK/3+NGy4v2zmG3/a+KUb0feRTQGl03lbJRzVetxTF1ag5+Oz5aCkpqL+wdk8NTxQBjV79gp24JoqDL+/CnzwHeT9oQenP1vhGd9e93XOfngTxi73LGtKSUH9Q3pa9Q+UGTZbTF5UgZMPl6H9jkqQpoEnJ1H46Ot6uTv3J9lp5sDtv7YKIzvntLu6NJfVP1Q28/fpPRVQS4vdE6wsResel1iJcPLhMu9vwhLN7C9teS4aPqNX7NNpueap69YqnHy4zLAsvOnYXQne7Nluqy7LQeO9ZQZrAG2fqQSVrzOcB+iVwXS5A0DTP5dDWxH4rwrXc2Ta9HHhSrFaUf+Ay7G0twzqknSP7dnrytB1i76umpGBhr3G+fJWFrMJqm7HRfMXK6Ctdm+ec2wvR+ftxvfi3OqLIDQ8UIaTD5dh5cvd6N1kweQl+r013roBHbfr+9oxMjJzXLgiTUPXJ6tgr5nNM21Yi7a79PV4atKtYgeA5i9UQC2ZbVKevKgCvR9zOYfT0tzPYYP6Ypp63lq03lUBNT0dDQ8Yl7tbvLFy5e46ElPX5+uQe3AU9N9vmRiVf0MfqQErQPpzr5odCgCg6YlarH26Afbe3qDWa/5yLdb+sA22ltaIxKGWFOPEncugThCK9gbeI5+2aiUa7ljtcYJ4Q1vK0HbpEqx80r0pofGpGpQ8ctjjim9i5zYMr9Cw9IfG2z/1ZC1KnjguL1zFEKV8PU5/IAurvhqZ5qLpc6T11hKMLneg4Lc2JP3+Da/LNz5dg7UPvAmemEDb3jqsfuksHIePB5ze6UfqUPR8D+wnGiIRvoewrtyJ6Fki6iGiIy7TsojoD0RU7/w/0zmdiOg7RNRARIeJaHMoAS//1v6wK3aurUTnF0PvJ0JdmoVTT/ruPMz6/KtIf+5V0LYKdNzjPa2Wx+qgFbn3Ke24oApdnwssvr7bajFx9TYAeoXWfp/neo1P1WDN42/D3tuLjnvrQFvLPZbR8paj+SueeSp8+EBQFXvjN2o9rtyHPlIzM3C0vf4U1n7xVaz5H4fR+JT3gX/PfbwWI9dun/nb1tYeUMU+ct12DHysBnzoqEfFDgBr7nnV8Kd8yssHsfSHB6CeV4LWhzzLsPi+A0FV7A3P1Hj95TZ5xVac+bT/zufa76sDbfJ/Feaq9cE6txvIth1b0HuncVqNT9UY3l+a2LkNfbfPrqNYrWj8hvs2mr9cC61Ab1vuv6UWY7uq9WUrzzPsg2XucXFmdy0mr9zmsZySlmZ4XLTfV4eGZ2pge++WmWmOI8c9KvaGZ2p8PhAxfY6o2UvR9IRnuZzYuwb5v+zEknr9V8PcsnC15u5X0fj4ZqiZmVj1tf1uFXvHPXVI+avnvUDX42L14/sNK3a1pBinH3Evw8an9cHWw60vpvm9ciei9wAYBvB/mLncOe3rAM4y8xNEdD+ATGa+j4h2AvgsgJ0AtgP4NjNv97btaZEYQ7Xz7jqs+MsA+NDR6cBBquqzzQ8Aeu6qQ/bbo1D+7vllQknJgfW94ict0jTDdthA4gMwW4E47F7Xc42VNA1stxu22xnGEiTDcnGN0d+yftbxy2U9ragATTetRP6Xg7iycylD0jQ0fnUbiu8Nvs9vn3kjAkjxmzdf+yrgdXyk5TVGg3XmLut2rAR5DM6sww7j49AgLtI0gBSwbcpnefg9L13379w8OR/F5anJ2fiAmbLQ8pbP3Ifxlx5pmn5fZs69M9fl2x6oQ8EL3bCfbDRc37UMp9cLpr7wdeUeULMMERUCeNGlcj8B4GJm7iSiPAB/YeZSIvoX5+fn5i7na/upK/P5Pb0XhFXpKFYreGws6G0oFgt4ciphOgozQlvL0bPNipzvOQ9YRUXz49UBN3+EZTqth18NqgILZvvqYktY3SyoGUtgPzcQwaBEqPpuq0XmiTHDi61IO3trLZY0TcChEUbykpHxkwOzN4mHRwI6R9SyUrRfthTLv2V8caGmp8M+PDLzJWq/ZDMGC1KQ+b8jc+5F44ZqrkuF3QVg+pmelQBcf9+3Oaf5lHJmMqSKvW1v3czNM8fQUEjbcIyMJGzFPv6+akxctQ307iks/2Uj1JwcvanKYcfaH7Tp5ee8kum9o9b9udvqCpz7eOh92qtLs/SmKmdaUanYAcBhD7v/HKnYY8ey3zRAe3N+eobMebEBSYfqkfp6PZb+50l9IrN+PE0ft35wYwtWPe/98c3pbU1LeqMe2S8F3/5OW8rQf3Nw52PYT8uwfukf9JlLRLuJ6A0iemNycthj/vj7qsG1vp9UWP3SWXTfvsXrW5xaQT76bwm9giJNQ/dnvbdz2d67BVOXbvE632yL32yD5a1WOEZHYe/ugWNwEKt+ewYAYGtpxeqXzuo/gQEs/3MPuLMHAOC4cBOGiizI/nuHxza7/6lu5t/cF55oWwXGrtHbZR2Dw1j5ygB69tQZtuePXLs9qLbm/ptroRXOPktPW8ow+iHfLX59t9W6D45csxHj768OOE0x/+y9vYYdeEUtraEh2AcHYT/T5zHf330opWoDht5XCVtnV8BpOoaGPB94IEL3P/luT1dOdSDnv3w2gHgI9X3/biLKc2mW6XFObwfg+nT/Kuc0D8y8D8A+QG9znzvfcuocaHQcvq7FHYePIzu9Sn9W3SiNgUFkHg/9QGG7HTmHx7zOT2nVb74F/5T8/LC1u1fOPDEB+7snZ/52vTnk2iaY3HYWSWeSYWtq8dhm9uHZsuYx97JRu/phGZuCA3qbJr19AtmLPR+/BABr/QCUvkGf+9dV5rFh8MDsM/pq9zksnrLD4WOdrHdHwUOzFw5aZz8sQ2kxu7/E/Jq6fCuG85KQeWIEat8wRkqzkfri6wGvr/Seg5XZ5zEYEGa388qIvb8fCPIprlDb3L8BoM/lhmoWM99LRFcD2IPZG6rfYWa/l0qRuKEK6H2XaL1Dpr7xZ6bx91Uj7c9HvL/pGwX2izcjpbEHtlb9J6y2aiUmSnKhvvKPeYtBiEC5niNcV4nxZSmwnBqE0j+MyeIcKH990+wQg+Krzd3vlTsRPQfgYgDZRNQG4FEATwD4ORHdBqAFwIedi78MvWJvADAK4Nawow+CbXES1MEY65hqHk1kKEhTQ+tvI1RTVg3JKUmzE5KTMJmuYdG8RkisVt8AAAW5SURBVCGigWsroR5umLdmkvngeo7Q/rexCIDD+U9p9d/GHgjHRZug/O0tn/eZ1PR02DcURqULk2l+29yZ+UZmzmPmJGZexcw/YuY+Zt7BzCXMfCkzn3Uuy8x8FzOvYeYKZvb+dkAUJP3+DTiOBP6CwXxQy0qj0tslaZpbp1wAsOT/vurRvUC0pf7H67A3NM38bTvVjEW/DvynbSDUDetk/FWngMqCCMrG9R6T1ZLioN7qPFueBlrsvQM10jS3ju58UTOWQC0p9r9glGW92Q8ej+5wj72Vi/RHK32gdCvObgiuc7pgxUz3A4mq9eqlUHycIKFSrFa0Xh35XiljUduV2VBC7B0z0QRSFpScjNNXe/ZK2XNxLpAfeNcFS39wAPbuHq/zlbQ0tBqkY4RX56Hnouh2lBWI1p3ROR9d5X5nv9/3G2xt7ch6NrqPIsdk9wNCCCH8i4uOw4QQQkSOVO5CCJGAYqJZhoh6AYwAOGN2LFGQDclXPJF8xZeFnq8CZs4xmhETlTsAENEb3tqO4pnkK75IvuKL5Ms7aZYRQogEJJW7EEIkoFiq3PeZHUCUSL7ii+Qrvki+vIiZNnchhBCRE0tX7kIIISLE9MqdiK4kohPOcVfvNzuecBBRMxG9Q0RvEdEbzmmG483GMjPGzZ0PXvL1JSJqd+6zt5xDRU7P2+vM1wkiusKcqP0jonwieoWI3iWio0T0Oef0uN5nPvIV1/uMiFKJ6HUietuZr8ec04uI6DVn/M8TUbJzeorz7wbn/MKAEmJm0/4BUAE0AigGkAzgbQAbzIwpzPw0A8ieM+3rAO53fr4fwJNmxxlAPt4DYDOAI/7yAb0X0P8EQABqALxmdvxB5utLAO4xWHaD83hMAVDkPE5Vs/PgJV95ADY7P1sBnHTGH9f7zEe+4nqfOct9sfNzEoDXnPvh5wBucE7/PoA7nZ8/A+D7zs83AHg+kHTMvnKvBtDAzKeYeRLAzwDsMjmmSNsF4MfOzz8GcI2JsQSEmf8G4Oycyd7ysQv64OnMzK8CyHAO4BJzvOTLm10AfsbME8zcBL0b65gcxomZO5n5H87PQwCOQR/eMq73mY98eRMX+8xZ7tOjyCQ5/zGA9wJ4wTl97v6a3o8vANhBROQvHbMr95DGXI1hDOD3RHSIiHY7p3kbbzbeRHTc3Bizx9k88axLs1lc5sv5k30T9KvBhNlnc/IFxPk+IyKViN6CPordH6D/yjjHzNODk7nGPpMv5/wBAH67hDW7ck80FzDzZgBXAbiLiN7jOpP131Vx/3hSouTD6XsA1gCoAtAJ4GlzwwkdES0G8AsAn2dmt1HD43mfGeQr7vcZM9uZuQr6UKTVAALrGD8IZlfuAY+5Gg+Yud35fw+AX0Hfad3TP3nnjDcbb7zlI673ITN3O080B4AfYPZnfFzli4iSoFeAP2XmXzonx/0+M8pXouwzAGDmcwBeAVALvXlsemQf19hn8uWcvwSA54jec5hduR8EUOK8S5wM/WbBb0yOKSREZCEi6/RnAJcDOAI9Pzc7F7sZwK/NiTBs3vLxGwCfcD6BUQNgwKUpIObNaWv+IPR9Buj5usH5pEIRgBIAkR1iKkKc7a8/AnCMmZ9xmRXX+8xbvuJ9nxFRDhFlOD8vAnAZ9PsJrwC4zrnY3P01vR+vA/Bn5y8x32LgzvFO6HfBGwE8aHY8YeSjGPqd+rcBHJ3OC/S2sT8BqAfwR+iDiZser5+8PAf95+4U9La/27zlA/qd/+869987ALaaHX+Q+fqJM+7DzpMoz2X5B535OgHgKrPj95GvC6A3uRwG8Jbz385432c+8hXX+wzARgBvOuM/AuAR5/Ri6F9GDQD+DUCKc3qq8+8G5/ziQNKRN1SFECIBmd0sI4QQIgqkchdCiAQklbsQQiQgqdyFECIBSeUuhBAJSCp3IYRIQFK5CyFEApLKXQghEtD/BwpA/BvzJUNAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "109"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "310"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.float64"
      ]
     },
     "execution_count": 183,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(dataset[10][10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_states = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "109"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_neurons = dataset.shape[0]\n",
    "num_neurons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = dataset.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import block_diag\n",
    "import autograd.numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import ssm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create our HMM\n",
    "simple_hmm = ssm.HMM(num_states,\n",
    "                     num_neurons,\n",
    "                     observations=\"poisson\",\n",
    "                     transitions=\"standard\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "90ff9058dbcc4c2d87ccfabf4cd3c127",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "lls = simple_hmm.fit(dataset.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_ = simple_hmm.most_likely_states(dataset.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7580321dd8>]"
      ]
     },
     "execution_count": 192,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO2dfdQkVX3nv99qBlCGl4EZAYcZBgmJgAZknwMkkqi7SgbOrpOc9ZyFuIbk6JmNK3nb7J7AukdY/GNjXkzWSCRsnEXzAiaoceKCSJSE3SDIoIiAQUbAMBMij4xgeAnDdP32j6rqrqefeum+dW/dvrd/n3Oe83RXdXfd6qr763t/L99LEYGiKIoSL4nvBiiKoihuUUOvKIoSOWroFUVRIkcNvaIoSuSooVcURYmcg3w3oIr169fLli1bfDdDURQlGO65557visiGqn1zaei3bNmCXbt2+W6GoihKMJD8dt0+dd0oiqJEjhp6RVGUyFFDryiKEjlq6BVFUSJHDb2iKErktBp6kptI3kbyQZIPkPyliteQ5IdI7iZ5H8mzSvsuIflw/neJ7RNQFEVRmpkmvfIAgF8Vka+QPBzAPSRvFZEHS6+5AMAp+d85AD4C4BySRwO4AsASAMnfu1NEvmf1LBRFUZRaWg29iDwB4In88T+R/AaAjQDKhn4bgI9Lpnl8J8mjSB4P4I0AbhWRfQBA8lYAWwFcb/UsHPLC/iGuu+MxvLD/wKp96w8/BO8490SQ9NAypQ/2H0jxv//2UTz34urrX8dxR74MP33OZoet8ktVn1h32MH42R/dYtQXvrX8LD5z7z8AHSTTzz/9OLxm45HG7zfh+f0HcN0dj+Gf9w+tfebLDzkIP/+Gk619XsFMBVMktwB4HYC7JnZtBPB46fmefFvd9qrP3g5gOwBs3jw/neTLj+3DBz73dwCA8j1c3JNvPvVYvPKol3lomdIH9+15Gv/j5tXXv47ivrjgNcdh3WEHO2yZP+569KkVfaI45zf+0Ctw0vrDZv68j93xGD7+pW9P9f1WIQJ8a/k5XP32s9pfbJG7Ht2H3/jcQwCmuzemYf3aQ/waepJrAXwSwC+LyPdtN0RErgVwLQAsLS3NzWooLx1IAQB/eel5eO0J4xHDJ+/Zg1/986/hwHBumqo4YP8wu/43bD8X577qmNbX//Gd38Z/+4v78VKaum6aN17K7/nP/sJ5eM3GI/HZ+/4Bl/7pV3FgaHbOLw1TvOLwQ/Dl977Z6P1bf/d2vGR47C4UtqH4HuaZqbJuSK5BZuT/REQ+VfGSvQA2lZ6fkG+r2x4Mw3y4MvmLnSQr9ytxUtjrZMohW/G6iO08hunKPlGcs2lfGKYy9fdbRUIi9dAPi2N2aXtfTJN1QwAfBfANEflgzct2AviZPPvmXADP5L79WwCcT3IdyXUAzs+3BUOx1OIgWXkxRx1aDX3UpKPrP93ri9fFfF9M9omuP26prO5fszBIiNTD110cs0vb+2Ia183rAbwDwNdJ3ptv+68ANgOAiFwD4CYAFwLYDeB5AD+X79tH8v0A7s7fd1URmA2FYc2Ibnxzx9uhlfKMbrrOXLxuGPF9MZwYyRZ2zvTHLU2lk487oZ/vuzhmAHZ+qqyb/weg8VTybJv31OzbAWCHUevmgLoRXfErHnF/VlAavU5piYrXRTygH93zhaEf9wVDQy/SaVScJJ5dNwFYeq2MbSGtGdEV1zbmkZtSP6OrYxFiN+nESDbpOIsZSjc/t/ro21FD30JaM6JTH/1iMB61Tff6Rbgv0kkffcfZbSrSyf0xIL0Ev4tjTjvb84ka+hZaffQRd2ilPHqdNesm3vtimNr30XcZFZN+ZlB1GXnziBr6FupGdOqjXwxmzaxYhPuisKnFSH7Q8cetq49+kHAUS+mTuoy8eUQNfQvFzbsqvTKJP7tCKWeYTPf6RYjdTH4no0wj4zz67j56P1k34+PPO2roW5hMJSvoOl1VwkBmDLgtgktvMm5VDIJMTzkVmToGUkXiLY9+tviNTwJool8mU8kKuk5XlTAY1szo6uiaahgCo7hF4bopMo26uG46jIoH9PN91yVqzCNq6FuYTCUr6DpdVcJgMvDYRtdUwxCY/E66u26kkwKsP9eNpldGw2QqWUHX6aoSBpOBxza6phqGwChAXbhuRkViZictHSUQ/Lluxsefd9TQtzAWcNKCqUXENBgbteumSCvMrcd4FmP2eZmomXl7EvpxodbN9ucRNfQtFP21Lusm5g6tzO6HXYTYzeR3knQUcssKprqKmnn00Qdg6dXQt1A3oluE7AplbLBnFjWL+L6YTCvsWiTW1dCT9FowpT76CKgLuAw6TleVMDAtmIrYzq9KK+xaJNZZppj08n2P4jdq6MOnVo9+AXTHldmlaBchdjMpCzE6505ZN+btUZnidtTQt6B69IvNrFK0ixC7mcy6STpm3YQuU6w++ghIa3z0i6BposwuRbsIsZtJMa+utQNdC6YS0lvWDTl9/MYnrQuPkNwB4F8DeFJEXlOx/78AeHvp804FsCFfXeoxAP8EYAjggIgs2Wp4X6RSfTG7TleVMBjOKEW7CLGbNE+HLPpEMegx1qNPuxnLgcdgbAj+eWC6Ef11ALbW7RSR3xSRM0XkTACXA/ibieUC35TvD87IA/Wjja7TVSUMZtajX4DYzaSrJekYgBaRqdfkrcJnwVQI8gfAFIZeRG4HMO06rxcDuL5Ti+aMOmW9RSh1V1SPvoqhrJQssBGM7aZe6a9gKhA7b89HT/LlyEb+nyxtFgCfJ3kPye0t799OchfJXcvLy7aa1RmpUdZTH/1iMBl4bGMR7guZGMkOOsYlMvXKMAumQgjEAnaDsf8GwN9OuG3OE5GzAFwA4D0kf7zuzSJyrYgsicjShg0bLDarG3WjjWJTzCM3ZfZVhBYhdjMpWcCOs5hUwtWjj8lHPy0XYcJtIyJ78/9PAvg0gLMtHq8XhjU++lEAKuIOreQzuhkyKxYhdjM5Au9eMCUYdMqj91Mw1XWt2z6xYuhJHgngDQA+U9p2GMnDi8cAzgdwv43j9YlIdQ511+mqEgaz+o8XIXYzucZr1yIxGz56HwOuri6nPpkmvfJ6AG8EsJ7kHgBXAFgDACJyTf6ynwLweRF5rvTWYwF8Oh8JHQTgT0Xkc/aa3g91ynpdp6tKGKQ1P/R1LIKPflKyoGuRWJoG7KMPxHXTauhF5OIpXnMdsjTM8rZHAJxh2rB5oS7gsggdWpl9er4IsZvhxHfStUisa5oiSaQe6ha65v/3iVbGtpBK9eo3i6BpomQGexYjtAhLCcpEodDYjWn2ecOOa8YOEj/fd9f8/z4JpJn+SNPq0cYiaJoohRGawdAvhEzxZMHUeLsJkz8cs+KtMnbGQYBP1NC3MDlNLVgETRNldeCxjUWI3UymFXYtEusajGWeddN3ptNk4dg8o4a+hbrIetfpqhIGs2qlL0LsZrKIsGtf6KxH7+k777rWbZ+ooW8hTauDsV2nq0oY1M3o6liE2M2kmFfx0NR90lVKwNd33nWt2z5RQ9/CsKZqbxE0TZTZ/ceLELuZzJIhiYQd9ei75NF7+s5DyqNXQ99CXXqdum4WA9OCqagNfcUIvIsMwawB70l8fedd17rtEzX0LdQF47pOV5UwmNlHvwADgKraki5SwWlHzRhf33ldRt48ooa+hbqCqa7TVSUMZvUf05O/uE+qZjkJu6lXdslH9/WdD2UBZYpjpan6zZdqntIfs0rRFq+NeQBQpTY56LCcX9eVmnx957KgMsVR0lT95mtlG6U/6oLxdYxFzVy1yD9pRSVrYli0JCKZcGBHmWLAV9aNGvooaBptdJmuKmEwq9ZN8dqY74uqLJkkMZMKLmxzJ0PvKY9+OKPgnU/U0LfQtChCl+mqEgZ1dRR1FLGbmA39sEJtcpCYuTGL76mLj96XZLh01NHvEzX0LTR19CTxo7Gh9IfJ9Dz22E1VWqGpJnzxPXWREvBbMBWGpVdD30JT9VuiI/roMVnmLvbYTVVaYbbKk4mPPvvfJaipBVPtqKFvoakoYhB5h1aqA49tRO+6qUgrNJ3FFLOALvZyXKVu/hkmZPn//R7TlNZbmOQOkk+SrFwGkOQbST5D8t78732lfVtJPkRyN8nLbDa8L5oMva8lzJT+MCnPjz12U5VWaDroKX4cuqVXZv99jOhjSq+8DsDWltf8XxE5M/+7CgBIDgBcDeACAKcBuJjkaV0a64OmykjT6aoSDsN0dila01TDUKjyTZNmuk8yCsZaSK/0IFMcjY9eRG4HsM/gs88GsFtEHhGR/QBuALDN4HO8MmyojIw96KaYSdGaphqGQtU6uqbrttoY0Rfv7XvQZRK/8YUtH/2PkPwayZtJnp5v2wjg8dJr9uTbKiG5neQukruWl5ctNas7TdVv6qOPHxMp2oRxSyBU1RZksxiTz8rfb2NE37uPfrFkir8C4EQROQPA7wH4C5MPEZFrRWRJRJY2bNhgoVl2aJqemU5XlXAwUSg0Hd2GQmXBlGEAOrUQjFUffTudDb2IfF9Ens0f3wRgDcn1APYC2FR66Qn5tqCYXDatzEDz6KPHxNCTcRv6Kv0n01TjUcFUx6UEs3ZpHn0dnQ09yeOYf9Mkz84/8ykAdwM4heRJJA8GcBGAnV2P1zdNWjcDqusmdiYXwp6GQeSxm6yIcOU208pYK1k3qkffykFtLyB5PYA3AlhPcg+AKwCsAQARuQbA2wC8m+QBAC8AuEiyqMgBkpcCuAXAAMAOEXnAyVk4pOlXW1038VMVeGwj9thNpR694aBHLPjofa0Z23Wt2z5pNfQicnHL/g8D+HDNvpsA3GTWtPmgqfotdl+sMruoGRD/ACArmJoUNTMbUY9H9Obt8aVH33Wt2z7RytgWmlKoNL0yfowKpiIfAIis9qkPDOMSY1Gz8PToFyoYGztpg0Kd6XRVCYemhWfqME01DIWqlFMaDnrGWTcB6tEH5KNXQ99Ck4/edLqqhENTML6O2LVuqtyZA8MisSL33Yah97FmrBr6SJCGYJzpdFUJB1OZ4ph99GnNmrG+9Oh9LfZiEr/xhRr6FpoqI02nq0o4mEjRxu6jTyt89InhoMeGHv0460Z99HWooW9h2CqBEG+HVsz0TLIBgKMGzQHDmjVjTfrCSI8+yIKpbj9QfaKGvgVp0qMne9fAVvqlKRhfxyDpPwOkT6r6hLFMcRGM7bKU4CjrxvwzTDCJ3/gikGb6o61gSiUQ4qZqfdQ2BgsgU7yqYKrjmrE2KmN9ZN10mYn0iRr6Fpqq37JMg3g7tFIdeGwj9thN1eDHWNTMggTCqGCqbz16g7UKfKGGvoWm6jctmIqfqsBjG6aphqEgFXEL84Kp/P0BFkyZrFXgCzX0LTRVRsa+CLRSHXhsI3Y9+mFFWqFpANrmwiN9B8BN1irwhRr6FoYN6XWxF8YozcH4OkwzUEKhKq3QNAAduh59FzG2PlFD30JTep0WTMWPccFUxPdFWqNH361gqnt6pcoU16OGvoUq7e2CLNOg3/Yo/WIiRbsYMsUrtyUd14ztVDDlzdB3y//vEzX0LTQJFyWRy9EqZlK0XAgfvV09+i4jevXRt9Nq6EnuIPkkyftr9r+d5H0kv07yDpJnlPY9lm+/l+Qumw3vAxGpzDAo0MrY+DGVKY417bauTwwM41U29OgTDz764vrG5KO/DsDWhv2PAniDiLwWwPsBXDux/00icqaILJk10R+jFeprC6biLoxRmoPxdSQR3xd1faKrj96KemWPs6jiXENx3UyzwtTtJLc07L+j9PROZIuAR0Gbst6AcedLK2ZaN0nE0hh1fSIxrB2wUhnrYSnBYYQj+ll4J4CbS88FwOdJ3kNye9MbSW4nuYvkruXlZcvNMqMtUBR7vrTSHIyvYxDxOgUjV8tkeqXxiD5/vw0ffa+um5XHnndaR/TTQvJNyAz9eaXN54nIXpKvAHAryb8Tkdur3i8i1yJ3+ywtLc1FL2kLFJlmGijhYLKKUMwV03Uj8CQxM7QjF4gNPXoPrptABvR2RvQkfxjAHwLYJiJPFdtFZG/+/0kAnwZwto3j9cVoeta0lGCkHVppD8bXEfMAYDQCr/DRdymYCk2P3kb+f590NvQkNwP4FIB3iMg3S9sPI3l48RjA+QAqM3fmlTb/4UDXjI2atmB8HTGvJTx2Z67c3rlgKjA9+iIGE4qoWavrhuT1AN4IYD3JPQCuALAGAETkGgDvA3AMgN/PT/pAnmFzLIBP59sOAvCnIvI5B+fgjDRt/tVOknizKxTzZe5MUw1DQGpGssZ69BbWjPWhRz/+gervmF2YJuvm4pb97wLwrortjwA4Y/U7wqFNcEkLpuLGtGozZh99XZ+gYV8YzZot+Oj7HHQtetZNVIym7rqU4EJiWrVpmmoYAnV9wlim2KJ6pQ8ffShZN2roG2hT1ot55Ka0B+PriDnttq5PmLoxbaZX9jm7Ti24nPpEDX0DbYGiRAumosZ01BbzTK+uT5gGoIejrBvzNvkomDKN3/gikGb6YRofvQZj46UtGF9HzDLFtQVTidmIWixk3Yx89F7y6HVEHzxS448siHnkpph35phdenUuC1N9n6Hhj2kZklkwWH30taihb6Ct+o0Ra5oo7cH4OmLWo6/VusndmLMWTdnQowf6XwTIRmyhT9TQNzCU5tHGwLDsWwkD02XuTFMNQ2BYM5IdZ77M9nk29OiL4/epR19XODavqKFvQFqmZ7qUYNyYVm3GfF/U9QnTdVtNM5smSQzXrDWlrnBsXlFD30Bb1R4Np6tKGBj76COumK7rE6YyBLZ83X3HRepmNvOKGvoGWvXoPaR1Kf3RFoyvI2atm7o+YSpDYKNgCuhfd0rz6CNiGj368uuUuDCVoo1ZGqOuT5jKENgKavrLuuntkJ1QQ9/AKFBUZ+g9yKMq/WEqRRtz2m1dnzCVIbCl6973d75wMsUxMxYuqt7vQ2ND6Q9TrXTmboQYYzd1fcJUhkBEQHZPr+zdR68FU/EwjR599rremqT0SN0iG20Ur4/Qztf2CdN41VDEygLbSc+1C2PVTTX0wdMWKKL66KPGdJm74vUxZt7UyUIUBm/WvjBM7YyKBz2v9mY6CPCFGvoG2sqzR6MYNfRRYlq16WPFo76oc1mM1m2d8cdNRDpp0ZeP36sefYxrxpLcQfJJkpVLATLjQyR3k7yP5FmlfZeQfDj/u8RWw/ugbSk5H2tVKv3RFoyvw8eKR31R1ycGHYKxNkb0fa/TG6vr5joAWxv2XwDglPxvO4CPAADJo5EtPXgOsoXBryC5zrSxfdOWQjUaucXYo5XWYHwdPlY86otaPXrjgik77o+kb9dNYHn0rUsJAoCI3E5yS8NLtgH4uGRpBneSPIrk8cjWmr1VRPYBAMlbkf1gXN+l0U38z796GH/3j9+38lnfffZFAA2um/wi/9qN9+HQNQMrx/TFIQcluPzCU3HsEYeu2vf08/tx1WcfxAv7h9aOd+amo/Af3nDyqu1f+MZ3cOM9e6wdZ1oOPijBr219NY5ZezCu+MwDeOaFl/C95/cDMFOvBIBfvuFerAllUdEpqesTxcj2vZ++Hy8/uL4v/MAr1uJXz/+h0fM0z7rpyiAh7np0H979x/cAAN7xIyfiR09e3/2DawhNj34qQz8FGwE8Xnq+J99Wt30VJLcjmw1g8+bNxg35vS8+jCNetgbr1x5s/Bllztp8FE5af1jlvjM2HYnTX3kE9j79gpVj+WL/gRSPPfU83vTqV2Dbmasvz1cffxqf+spebDr6ZXiZhR+073z/RXzpkacqDf31X34ctz+8jC3HvLzzcablpaHg0e8+hx87ZQNOO/4I3HD343jlkYdi7aEH4YxNR+EHjz18ps8768R1OO34I/D3+55z1GK/VPWJ1248Eq/ZeASeeKa+Lzz17H7cfP8/4lfe/IMralBs5KKff9qxuO2hJ/Gt5WfxyPJzWHvIQU4N/dAw9dYXtgx9Z0TkWgDXAsDS0pLxHCwVwdvP2bxi1OCK0195JP7PL/6Y8+O45tHvPoc3/dZf1/qUi3zwD198Fs7YdFTn41258wF86ivVo3YRwQ8euxaf/YX+vtc933se533gNqQio5HalW89HeeffpzR5521eR1u+qXw74tZ+KHjDm+9Zh/6wsP44K3fRCqCBGNXjw33x+UXnorLLzwVAPD6X/+i81RLGwum9ImticdeAJtKz0/It9Vtd4KIIJVwfmXnhTYphzZxt9mPV5/zPBQ7HX/W9gBZ9lRoC0qERFXsIhX7Ac0kcZ8gYbtPuMaWod8J4Gfy7JtzATwjIk8AuAXA+STX5UHY8/NtTggtt3VeSFqCyuPl4+wcb5A0/aj0b+gL18FQxMqKR0o1I3dNSTc+TcV6iuKghyrZ1DBQ74upXDckr0cWWF1Pcg+yTJo1ACAi1wC4CcCFAHYDeB7Az+X79pF8P4C784+6qgjMuiC0AMm8ME4HrO4ctrW3m1LhRPo3suVFM0xXlVLaqUrBtFUZW6aPVEvT9YR9MW3WzcUt+wXAe2r27QCwY/amzY6tZckWjXFqXPV+29rbTYtnDx2M8Nrbk/1f6brptw2LQNXMMRWx77rpYeEX1aP3iK1lyRaNpGV1oLbCsVlp0g5PPfjoy4Vvo5FaIB04JAqDLqtcN3a/60EPaznb7hOuicrQ21qWbNFoU+Eca/7YOl69j96HoS9LFoSWNhcSdcFY2wMz9iCHYLtPuCYqQ68ZE2YMSlknVdjW3k4aYgIuOn4bZckCnRW6o0oyZGipYGryOK4lolWP3iO2liVbNMa+0+r9trW3m8rlh6n9jt/envzYpaybQPpvULBiQCEugrE9ZN2oHr1HbC1LtmgUPvr6rJvidZZ89A3a5WKpUnIWyq6r0MSqQqJq/QYX6bR9aNPb7hOuicrQ62jMjDZBKtuxj6JfV8UE5qVgSoOx9qnS6XdSMMUeCqYCiwdGZehFR2NGtK0OZNv4NUnappYWopipPaXzD63iMSSqXDeuCqac59EHNiCIytAPA/vy54Wps24su26qZhCZyJWVw0xNWQIitIrHkKgtmLIugdBDZWxgNTtR3c6hBUjmhXatG7vf63hkV30sH+mV5Mo8er2H7FP8eJbvs1Tsf9cJq+8tm4QWD4zK0IcWIJkX2lbKsq0hVEi0V7puHFRKTkMx3Q+tA4dEUhGMdeK66UECYaSJFMiAICpDr8FYM6p8p2WKTkNLd8toIemaPHofo+ksJS+8IFtIVLkIbenRTx7HecGU5T7hmkCaOR2hFTHME4OGlDTbgaemmEAqAh+LMiVJFswfBfQDGamFRGXBVCrW/dxNMti20GCsR0a/soF8+fNE0lA27kKPHpgfHz0wLrLROI87qmJBIvaNZeajd+26KY4Vxn0SmaHP/ofyKztPNCn+2c5EGTSIqImDvOppKITWVI/eHYVRLF/2oYj1DKc+fPShZWcF0szpUB+9OQlZ76N3lHVTJ4Hg4/oVWTca0HdHVVqtCxE79rHwSGAzv6gMvZavm9Pso89fY7lgqmrQ5SI4Nw3FKFCDse6oDMY6kil2PKAPznswlaEnuZXkQyR3k7ysYv/vkLw3//smyadL+4alfTttNn6SNDC/2TxBtksg2Ppak4pS+IJU7AfnpqHw0YcWZAuJpCoY66Rgyr1Mse0+4ZrWFaZIDgBcDeAtAPYAuJvkThF5sHiNiPxK6fW/AOB1pY94QUTOtNfkekaVsVHNU/qhSdpVJHOn2DLAzVk3foxssfxcaBWPITFayWtFHr2Lgin3PnrbfcI105jEswHsFpFHRGQ/gBsAbGt4/cUArrfRuFlRPXpzBg25x8PU7qhrlGZX46P34rrJVyXSgil3DCpiM6k4KpjqQaY4pHtkGkO/EcDjped78m2rIHkigJMAfLG0+VCSu0jeSfIn6w5Ccnv+ul3Ly8tTNGs1oQVI5gnmBUNVDC27U6rWDi1IPejRA+P0Ug3ou6NS1CzQginbfcI1tp0cFwG4UUSGpW0nisgSgJ8G8LskT656o4hcKyJLIrK0YcMGo4PraMycQdKsR2/TndKUR586WIhiGkauGw3oO6NKJdWJHn0Pa8a6yP93yTSGfi+ATaXnJ+TbqrgIE24bEdmb/38EwF9jpf/eKsM0rADJPNG0Ko/tlMexr7ZGj96DkS3SS9X9546qNWNd1E30okfvKQ3YlGkM/d0ATiF5EsmDkRnzVdkzJF8NYB2AL5W2rSN5SP54PYDXA3hw8r22KEakIf3SzgtNZeO2hcaaRNR8ad0U6aWF+0rvIfvUZd2EKGrmS3zPlNasGxE5QPJSALcAGADYISIPkLwKwC4RKYz+RQBukJXz/1MB/AHJFNmPyq+Xs3VsM9RptzFJ0qxHb9P4Nq1o5ULNcBqY++hDq3gMiapF6F246priTbZwkf/vklZDDwAichOAmya2vW/i+ZUV77sDwGs7tG8mivsnpAswLzStypOK3bhHUuGrHR/LX9aNqB69U6pliu2nKDbFm2xhu0+4JqpxS6oZE8Y0+ugtT6/rfPSS68H7LZjKnqvrxj7VC4/YX1Gsr6ybkOxMXIZeZYqNSZL6snGxrEdSNYXPjrNyf58khY8+sIrHkBiLmq2UKXaRXuk6j952n3BNVIZeJWbNSZokECz7I0eiZhO/LD51Zgpp29AqHkNiJGo2sfBIiHr0vuS0TYnK0GtqnDlN091haneWNK6MnTyOv2D6IOGoYErvHzdU6dG7kLwYJPWDFluoj94jWjBlTqvWjcU7pU6PXjxev2IUmHrSw18EKvXoHWRZJX2kV6b2dfRdElBT29HydXPag7GL4bpxob2iZFSl1brIR+9D1Mx2n3BNVIZey9fNKYKRVdieXg8qgnLZcfy53kZ69KkfCYZFoKpQzpUevfs1Y8PKzIrT0Ad0AeaFprJx20Jj45Hd6uOU9/dJsSpRaBWPIVFVGevC192UWGALX+J7psRl6LV83Zjmgim7KXBJjY/eZ4ylWJUotIrHkKjSo89UIC0fJ6meMdrEV2GfKVEZes2BNqdd1MyFeuWEj95jjKVYlSi0bIqQqNKjFwcSCE0SG7YILTsrKk30ySMAABMISURBVEM/EjXTjjozmdZN9T7bQmNVcrXA+Pp5U6/M14zV28cNrIjNuCiYqru/bOJLfM+UqAz9UNeMNaapmjC1nF5ZJVdbfu7j+hXnH1rFY0iMCqby+8yV5AVHLiK3rhtNr/TEWL3Sc0MCpEna1bbCYFUpfHacvC1eZYrV0Lti0kfvSvJiJLHh2kcf0H0SlUlUPXpzssrY6n3D1I0e/aQPdZR148l1M0zFehWwMmYy68ZV3UTd/WWTYbrYSwl6RbVuzCkKhqpILbsz6oJlXoOxeXqp7SpgZcxkMNZV3ct4bVqrH7uCKLNuSG4l+RDJ3SQvq9j/sySXSd6b/72rtO8Skg/nf5fYbPwkIz36gC7AvNDoukktrxk7Sn+bOI7HYPqoYEp99M6Y1KNPHcXUBn346C33Cde0LjxCcgDgagBvAbAHwN0kd1asFPUJEbl04r1HA7gCwBIAAXBP/t7vWWn9BKpHb05RMFSF7VznumBs0TF969GH1IFDYrJ+YvzDbvs41RIbNnGR/++Sab7iswHsFpFHRGQ/gBsAbJvy838CwK0isi837rcC2GrW1HZUj96comCoCrE8Ta0LlvkMxhZ6/KFVPIbEZP2EqyyrpIdgrO0+4ZppDP1GAI+Xnu/Jt03yb0neR/JGkptmfC9Ibie5i+Su5eXlKZq1Gp/peaFTFAxV4UqPfq4KpjheMzakDhwSIx99fp+JI9dN0oOPPrTsLFuTpr8EsEVEfhjZqP1js36AiFwrIksisrRhwwajRhR2KqQLMC80Kf7Zlu6tK2jxKUpXSECE1oFDYpzfnv0fOpqB18lg2yQ0OetpDP1eAJtKz0/It40QkadE5MX86R8C+BfTvtcmKlNsTmvBlAsf/ar0ynFb+oYk0jS8iseQILkiu8tVfx3JYDtMrwxNznoaQ383gFNInkTyYAAXAdhZfgHJ40tP3wrgG/njWwCcT3IdyXUAzs+3OUF99OYUBUNVWC+YqlAxLD+3HZybhkGSHT+0isfQKM8cXUlejGWwrX7sCkIrmGrNuhGRAyQvRWagBwB2iMgDJK8CsEtEdgL4RZJvBXAAwD4AP5u/dx/J9yP7sQCAq0Rkn4PzADAOpIVUyDAvsEHadZja/U7rgmVjUTqfBVNhdeDQSJLxkpXOgrH5D7XTrBvLfcI1rYYeAETkJgA3TWx7X+nx5QAur3nvDgA7OrRxajQH2pws66a6Y2QZBnaPBdSLmvnKuhmP6PUeckU5u8tVllV/WTfOPt46ATW1Hc2BNmeQNC0O7kaPfnVl7LgtfVOsSmS7ClhZSXlRkMJXb/vrHi8+71YCISQXcVyGXnOgjckKpqr3ZcUhDlw3NemVPq5hYYBCq3gMjSThKgkE+ytM9VUwFc59Epeh1xxoYwZJ/Yo8YnmmNM+um9AqHkMjKbkIXWlT9ZFHb7tPuCYqQz9MNTXOlEy9sqlgyt6xiktUq0fva+GRXI9eBwvuKLsIXWlTjeWQXa8w5ezjrROVoQ8tt3WeaF94xG5lbMJ6PXofP9aqR98PmUpo9jh1NIMb1KTv2iS0oH18hj6gL3+eSNiQR+/A+FWtUetTlI4jCYSwKh5DozygcFUw1ceasaEtIh+doQ/JbzZPFAVDVbjIZkoqCrS8yhTnvmOdFbqlXDDlSvJiXJBn9WNXEFqGX1SGPrQihnmiaoRdMHRQLVos9LHiOB4XjtGCqX4YJOPsLleSF7346AOroA6oqe2kaVhFDPNEIdNbhYsFswcVMQGfPvpihqGuG7ckpewuV5IXo6wuh0P60BaRj8osquvGnEKmtwoXAcqqLJ/xVN7qoaZsT/Z/mKbqunFI+bq7krwYiZo5z7oJ50aJytCHVsQwTwyaXDcOqgCTZPWIvji+jx/r4pgvDTW90iXl+8xV3cS4Mtbqx65AK2M9IqLKlaaM13FdbezFgXRvlVqmTz364pgvDVMdLDiEHKtKDl2tGduDHr2LPuGSqAx9aEUM80RTStrQQSZKlatoZOg9BWMB4MBQ3X8uGVRIINh21fXiugksOysqQ6959ObUrfqUbXPguqlQyyym2l5cN3lPOKA+eqesSK905Kob69G7LZgKyXsQn6HX0ZgR42XeVneO1EHaalU65zg4Z/VQU7cHyHz0OlhwR9nQu5K8GM9OrX7sClz0CZfEZehVedCYsdBYhaF3oL1d5aMfBec8ad0AwIFhqoMFhyRJWQIh3+Zo4RHXEgghpXJP1VSSW0k+RHI3ycsq9v8nkg+SvI/kF0ieWNo3JHlv/rdz8r02UeVBc9p99LZT4Kpkile2pU+K35aXtGDKKeWsG1eSF3Uy2DYJbZGj1hWmSA4AXA3gLQD2ALib5E4RebD0sq8CWBKR50m+G8BvAPh3+b4XRORMy+2uJA0s5WmeqCsbFxGHWTfzk0df3DcHhqm6bhySJKslEGz32aZ4kw1c9QmXTNOlzgawW0QeEZH9AG4AsK38AhG5TUSez5/eCeAEu82cjtACJPPEqGy8p2rVrHBm5TafWTfk2DjoLeSOFT56Z3r02X9XWTc+K7hNmcbQbwTweOn5nnxbHe8EcHPp+aEkd5G8k+RP1r2J5Pb8dbuWl5enaNZqhhJWgGSeqJN2dVWmXqV14yoLYxrKAwQdLLgjk77IHrscRADusm5c9QmXTLU4+LSQ/PcAlgC8obT5RBHZS/JVAL5I8usi8q3J94rItQCuBYClpSWjKyQiGGgfNaJu+bXRqMvB9HqVj97RQhRTtadkbEIaqYVGIQcNuHPVFT/UrmSKXfUJl0zzFe8FsKn0/IR82wpIvhnAewG8VUReLLaLyN78/yMA/hrA6zq0t5HQ9Cfmibrl11y5U+ZRj75A7yF3lH/gXS084lqP3qeL0ZRpDP3dAE4heRLJgwFcBGBF9gzJ1wH4A2RG/snS9nUkD8kfrwfwegDlIK5VtGDKnLqy8aKvuOiMc6VHn5RH9L0ffmGo8tFbr9EYyXlY/dgRrvqES1pdNyJygOSlAG4BMACwQ0QeIHkVgF0ishPAbwJYC+DP84v29yLyVgCnAvgDkimyH5Vfn8jWsUqa+snYiAHWjILGndHu8ZKKhU6GHkdK5WOqj94dSTIOwrv6YXcdjHXVJ1wylY9eRG4CcNPEtveVHr+55n13AHhtlwbOQiqCNWrpjRiXja/c7qqIacDV6ZXFU1969FWPFbsMWNKjdyR50VT8ZwOfhX2mRGUVQytimCeK38faYKwDzfC62YOP/lM+ZkD9NzjKsRlXkhejVFnXwdiAbE1Uhj4NrIhhnkhqRkGjFDgHWTeTAy7fa8ZWPVbsUl4r2Nls0XHBlKs+4ZK4DL3KFBtTVzY+zjCwfbzV8YDUUXBuGsrH1FoMdyQl6QtXkhfj1cJcZ904+XgnRGXoQ1v1ZZ6oGwW5TIGrmj34un5aMNUPZekLV3n0YzkPxwVTAQ0IojL0KlNsTt0oyF2Z+mpD73MxB/XR9wNLa8a6rNEof75t1EfvGTX05tR1DnHoo6+aPfi6fpp10w8DjmMzrhcecSdqlv0P6T6JzNDrtNuUOq2b0YLdDqbXVT56b64blUDohXJsZiR54UiP3rUEQkhaNwE1tZ00VT16U+rKxl0VMVWJmg1Tf0Y20aybXij/wI8kL6wL5jlOr4xUAiEYVKbYnCY9esCBHn1lMNajjz4pP9Z7yBWD0lrBzvTonbtu1NB7RQumzBnp0VeMsrP9LgqmVm7zqVWUrHDdeGnCQpCUgrGuRsbFx7mTQMj+h2RrojL0qcepf+gMaqa7rrS3B8lqvfBU/C3jp+mV/bCyYCrf5mAQkVD16MsE1NR2fE79Q2ckalYTjLWuMFgpgeCvWKl832jBlDtWFky5Kzyqur9s4apPuCQqQ68FU+YU39tqUbN8vwMp2dWpnOJtlKTB2H6oKphy0WeTivRdW7jqEy6JytCnEtav7DxRWzDlqnqxQo/e58Ix6qPvh/JI26XkRVVWly1c9QmXBNTUdlKPI8LQqSsbd1W9OKjoiD5F6VYsPKKW3hlJuWDKYd1Ltjata62bcO6TqMyiz2Be6NSqVzqUQFhVMCXibZSkSwn2Q0KsyLpx9Ztazu6xjas+4ZKpuhXJrSQfIrmb5GUV+w8h+Yl8/10kt5T2XZ5vf4jkT9hr+mqGqajrxpDiB3Iy5XFcBejAR1+hqzMfWTdemrAQDCYKplwZy6r7yxau+oRLWm9pkgMAVwO4AMBpAC4medrEy94J4Hsi8gMAfgfAB/L3noZsjdnTAWwF8Pv55zlBVALBmGIkXatH76RgauU2n3n0KoHQD0lSdt24S56o0lKyhas+4ZJplhI8G8BuEXkEAEjeAGAbVi7yvQ3AlfnjGwF8mNnQehuAG0TkRQCPktydf96X7DR/JUPVozemuGmv3PkAfuuWh0bbn98/zPdbPl4CLD/7It7ywb8ZbXvimX/GcUceavdAU0I19L2QENg/TPGWD/4Nlp990d2InsBn7t2LOx95yvpnu+oTLpnG0G8E8Hjp+R4A59S9Jl9M/BkAx+Tb75x478aqg5DcDmA7AGzevHmatq/iJ04/Fqe98gij9y46J29Yi4vP3oxnXti/at/rf+AYnL7xSKvH+6nXnYDvv3AAgvGw65Rj1+LHTtlg9TjTcuIxL8fbz9mMF/YPcc5JR3tpwyJwwWuOx7efeh6pCE45di1eu/EoJ8f5+TecjK/8/fecfDbgpk+4hG3VYyTfBmCriLwrf/4OAOeIyKWl19yfv2ZP/vxbyH4MrgRwp4j8cb79owBuFpEbm465tLQku3btMj4pRVGURYPkPSKyVLVvmrDTXgCbSs9PyLdVvobkQQCOBPDUlO9VFEVRHDKNob8bwCkkTyJ5MLLg6s6J1+wEcEn++G0AvijZVGEngIvyrJyTAJwC4Mt2mq4oiqJMQ6uPPve5XwrgFgADADtE5AGSVwHYJSI7AXwUwB/lwdZ9yH4MkL/uz5AFbg8AeI+IDB2di6IoilJBq4/eB+qjVxRFmY2uPnpFURQlYNTQK4qiRI4aekVRlMhRQ68oihI5cxmMJbkM4NuGb18P4LsWm+MDPYf5IYbz0HOYD1yfw4kiUllaPpeGvgskd9VFnkNBz2F+iOE89BzmA5/noK4bRVGUyFFDryiKEjkxGvprfTfAAnoO80MM56HnMB94O4fofPSKoijKSmIc0SuKoigl1NAriqJETjSGvm0B83mG5GMkv07yXpK78m1Hk7yV5MP5/3W+21mG5A6ST+aLzhTbKtvMjA/l1+Y+kmf5a/mYmnO4kuTe/FrcS/LC0r7eFrqfFpKbSN5G8kGSD5D8pXx7MNei4RxCuxaHkvwyya/l5/Hf8+0nkbwrb+8ncrl35PLtn8i330Vyi7PGiUjwf8jkk78F4FUADgbwNQCn+W7XDO1/DMD6iW2/AeCy/PFlAD7gu50T7ftxAGcBuL+tzQAuBHAzAAI4F8BdvtvfcA5XAvjPFa89Lb+vDgFwUn6/DebgHI4HcFb++HAA38zbGsy1aDiH0K4FAazNH68BcFf+Hf8ZgIvy7dcAeHf++D8CuCZ/fBGAT7hqWywj+tEC5iKyH0CxgHnIbAPwsfzxxwD8pMe2rEJEbke29kCZujZvA/BxybgTwFEkj++npfXUnEMdo4XuReRRAMVC914RkSdE5Cv5438C8A1k6zIHcy0azqGOeb0WIiLP5k/X5H8C4F8CKJZPnbwWxTW6EcC/It2slh6Loa9awLzpRpk3BMDnSd6TL5IOAMeKyBP5438EcKyfps1EXZtDuz6X5m6NHSWX2dyfQz71fx2ykWSQ12LiHIDArgXJAcl7ATwJ4FZks42nReRA/pJyW0fnke9/BsAxLtoVi6EPnfNE5CwAFwB4D8kfL++UbG4XVB5siG3O+QiAkwGcCeAJAL/ttznTQXItgE8C+GUR+X55XyjXouIcgrsWIjIUkTORrY99NoBXe24SgHgMfdCLkIvI3vz/kwA+jewG+U4xpc7/P+mvhVNT1+Zgro+IfCfvrCmA/4WxS2Buz4HkGmQG8k9E5FP55qCuRdU5hHgtCkTkaQC3AfgRZO6xYtnWcltH55HvPxLAUy7aE4uhn2YB87mE5GEkDy8eAzgfwP1YueD6JQA+46eFM1HX5p0AfibP+DgXwDMlt8JcMeGv/ilk1wKY04Xuc5/uRwF8Q0Q+WNoVzLWoO4cAr8UGkkflj18G4C3I4g23AXhb/rLJa1Fco7cB+GI++7KP70i1rT9k2QTfROYTe6/v9szQ7lchyyD4GoAHirYj89V9AcDDAP4KwNG+2zrR7uuRTadfQuZ3fGddm5FlI1ydX5uvA1jy3f6Gc/ijvI33IeuIx5de/978HB4CcIHv9udtOg+ZW+Y+APfmfxeGdC0aziG0a/HDAL6at/d+AO/Lt78K2Q/RbgB/DuCQfPuh+fPd+f5XuWqbSiAoiqJETiyuG0VRFKUGNfSKoiiRo4ZeURQlctTQK4qiRI4aekVRlMhRQ68oihI5augVRVEi5/8DC8CU+iQANGYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(states_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [],
   "source": [
    "post_trial = 7\n",
    "\n",
    "trial = conditioned_trials[post_trial]\n",
    "visual_time = trials['visStim_times'][trial]\n",
    "cue_time = trials['cue_times'][trial]\n",
    "feedback_time = trials['feedback_times'][trial]\n",
    "\n",
    "t0 = visual_time - pre_stim_dt\n",
    "tf = feedback_time + post_resp_dt\n",
    "[dataset, time_bins] = generate_spike_counts(recording_name, brain_region, neuron_min_score, bin_dt, t0, tf)\n",
    "(n_neurons, n_bins) = dataset.shape\n",
    "\n",
    "post_dataset = dataset\n",
    "post_dataset = post_dataset.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior = simple_hmm.filter(post_dataset.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f7580269c18>"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOy9d5gc1Zm+fZ9Ok/MojxhlCQlQBBQIIogFbBOMMckLrL3LOgnnXbzB9mL7W+z1j7W9a/BiG4OxCQJsokAmiCgJJZCEchxplCbn0Ol8f5yq7p7OaTpI574uXa2pqu6qDvXUW8/7nvcIKSUajUajOT2wZPsANBqNRpM5tOhrNBrNaYQWfY1GozmN0KKv0Wg0pxFa9DUajeY0wpatHdfW1soJEyZka/cajUaTl2zatKlFSjki2ednTfQnTJjAxo0bs7V7jUajyUuEEA2pPF/bOxqNRnMaoUVfo9FoTiO06Gs0Gs1phBZ9jUajOY3Qoq/RaDSnETFFXwjxsBCiSQjxcYT1QgjxSyHEPiHEViHEvPQfpkaj0WjSQTyR/iPAlVHWXwVMNf7dBTyY+mFpNBqNZjiIKfpSyneAtiibXAv8QSrWAZVCiDGxXre1xxn/UWo0Go0mLaTD0x8HHAn4u9FYFoIQ4i4hxEYhxMaTnX1p2LVGo9FoEiGjiVwp5UNSygVSygVeBB6vnsBFo9FoMkk6RP8oMD7g7zpjWVQk0NQ9kIbdazQajSZe0iH6LwC3G1U8C4FOKeXxeJ7Y2N6fht1rNBqNJl5iNlwTQjwBLAVqhRCNwPcBO4CU8tfASuBqYB/QB/xdvDtvbO/j3AnViR+1RqPRaJIipuhLKW+JsV4CX0lm541tOtLXaDSaTJK1Ebk2i9D2jkaj0WSYrIm+w2ahsUOXbWo0Gk0myZro260WHelrNBpNhslqpH+so1/X6ms0Gk0Gydp0iQ6rBZdH0tQ9wJiKomwdRuIc+xCkF8bNH7rcPQg7X4SiShgzF0pq1PLeFjj2ETh7YOa1IETk13b2we6VMPlSKM5wVZPbCa5e9T7cA1A+Dqz26M8Z7IHWfdDbrP4VVcH0q1I7BiTYCpJ/jXjZ/yaMOBPKY3YM8dPeAAffgap6qJ4EZWPBkoG4SUr1vdgLh39f8eL1QOt+6Dmp/rkH4JybwZqApBz+AMpGQdWEYTtMTShZE3271cIgqlY/ZdFvWAubfg/OXvXP41I/pBHTYMQMmHABOErif722A3DwXZh/x9Dlgz3w2PXQ3w6zb4Fl90LpSNj3Oqz8J2jb79+2Yrw6Wbsa/cvueBEmXhR5vx8+Bq/8E9gK4awb4Ny/h3EZaFradgAeWAzuALtt9i1w/a+jP++Rq+H4lqHLvrIeRkxP7jievgMGuuDOl6JfHIPxesBiDV3+9n+p7yT4fbTshcc+DdOuhFufHLru5A547kvwmYehZvLQdc99GRre8/9dVAW3Pw9jZsd/rMmw9Sn1+/rKOigfO3SdNO6UE/m8AvF6YdW/wFmfhvHnxf+8174Ha/936DKPExZ8Pv7XWPG3UFEHf//G0ONvb4DHb1Kf/9RlMGUZVITt7KJJgqzaO6Bq9VNCSnj5m7DrZSVezh7wumHfa+qH+fhn1WMivP8LePFu2L966PLNjxqCfytsewb+Zz784Vr44w1q/S1PKWFf9kN1Ap2xEK74Edz+ApSOgnf/X/T9NqyB0tFKcLc/B7+5BNbFEN500NmoBP+8f4RP/re64Gx9SoljJHpblODPvxO+8Br84ztgLYD1DyV3DFLC4XVKVA+sjr7tYLf6/F+4G345F/5zPLTsG7pNTzO8+zPY8gQceHvouvd/DkjY8wqcCOoY/sZ/wPGPYPX/N3T5kfXq2C6+Rwn9J+4HYYE37k3q7SZE0w4Y7IT3fj50uZTwxM3wxC1+8U+U/W/ABw/Ci19XF4B4OfYhjJypfttf/gDqzoV37zfu1uJgoFPdIRzdBAeDvp+3f6LO5WMfwotfg/+eCW/+KP5jyyV2vQyr/hVW3A4PLYVnPp/Y5zwMZE/0rYbop1qrf2C1Oimu+il8eS38/evwhVXw7T3wz4dg1NnKgkiEQ++rx9d/4P+C3E5Y879QfwFc/yB8eZ36oR/7EC77ntr39CtVJL/kbhUpfuZ3sHg5TLoYFn0VDrwFjZvC71NKOPKBuiv51M/hWzthwoXqBBjsSfLDiROPcaKedYOK1K78ibrbeOdnkZ9zeJ16nH2LusCNmQ1nfwY+egL6OxI/hp6T0G80c33rJ5FFzOtVd1vPfkFdGGuNu4rVQaLwwa+VJVJcA6t/7H+9zqOw5Sk4+7PgKIX3/tv/nCMbYM+rUHEGfPwsNO30r3vvv1Vkv3g5TFoK534BlnxN3eUd/iDx95sIPc3qcdMj0BUw2P3jZ9Xx7nkFdjyf3GuvexAsdmjarl4vXlr3w9h56rc9cgZc/M/QeQS2Phn7uQBtB/3/DwyGWvapC/W5X4BvbFfn2dS/gbW/Su53lU0Gu+Gpv4X1v4GT20FY1Wcc72c0TGRN9IWAEWUFqVfwrP0VlIxUghNMUZW6Rew4ErouEt0noXUvjJmjIr4dz6nlW5+C7mNw4TfU37VT4G//DPcchgu/FduHXvB3UFgJ790ffn3HYeg+ru4OAAor4LLvKyHc8Jv4jz8ZPG71aHr4pSOU+G9boU7ucBxeqyL7sXP9y867S+UFPvrT0G1dA+pfNE5uV48zr4Uj65R3Ho7tf4bGDeoi/88HlT2z6Muw/S8qdwLqZNvwGzjzk3DJv6iL6b431Lq1/6tyMpf+m3qP2/+sokpQF47iWvi7leqC8NZ/Gse2Q+Vazv8iFJT6j+Xcv4eSEfBW0F1BuultgrIx6g72/V8Y77EH/vpv6mI76iz467+DK8FzqXmPivQv+o56jdU/VtZoLAa7oecE1EzyL5tyufotvPOz+F7DtEJn36q+6yMb1N9v36cCjgu+oURi5JlwyXfB1QdbsiuWCXPkA5AeuPUpWL5J3RHXnauchyxewLI6XWJdVVFqtfpNO1Wkdd5dkUW3ok7ZF/He/jYYUf7V/wUjZ8GbP1QR4/u/gNHnwOTLkjvWgjI4/x9h10vQtCt0/REjWhx/vn/Z+HPVyfT+L4c32jcj/cDE7eK7weqIbEkdXquS2YGf+9g5cMYiZfF4PWpZ+yH4n3nw1G3Rj8GMqq+8Twnc2z8J3cY9qOyUUWcrwTV9/MXL1QXVtAA2PaLsgyXfgLm3q8h99Y+ht1WtO/tGlYxd9BUV5b7/C5XDOfAWXPhNqBwPC7+koufjW5UdZC9Rv7NAHCWw5OvqeQ1ror+/VOg5qX57s29Wuavuk+p76T4OV/9MfWadh2HN/yT2uh/8Wl24z/0CXPrv0H4QPvxj7OeZF8maKf5lQsBF/wQdDbDt6fhfY9m9Kjh77351cd32jDpPSkf6tx07V/3WNj6cvI2VDRrWqui+7lz1t8Wivq/eFn9AkQWyLPrFqUX66x5QUUG05FHlGeAZVNUl8dCwRp3gY+fC5d9XP86nPqeifzP6SJbzv6heO9BSMDm8DhxlMGrW0OVLvzv80b7XiMysDv+yslHqc93ypP8ENXH2Kj/fvCsJ5Px/VEK/96/QdQwevUY97nsdjm6OfAxNO1Teo3ys+pwb3odD7w3dZsNvlahcce/QxG1hhXrOvteUAK99QFljdfPB5oCLvwPHNqvv0dUHF3zdeI+jYe5t8NHjsOq76mJj/pYWfUW97svfUkK04O/CV1Qt+Lw67uAcQDrpaVZ3Xxd+S0XRK7+t7lhMa23iheoO6d37VYATD/3tykY5+0YoqYVpf6MCjrd/GvuOwbRLq4MS3dOvUhfkd37mv+hHou2g+rxLR6jzYvdKeP4r6g5r8d2h2y/4ArTsDv1N5DINa1QgFHh3OHaO+s2sfyg0n5Qhsh7pJ12r39uivNnZt/jLI8NRUace47V4Gt6HM85XUe/UK+CMxUrAqiepEysViquVeGx7WgljIIfXqcg+uAqlboGqXkhDtP/bbb9lw4kNoSvM23FLUDHX4rvVsuBo/+gmZTWcsSj0tWZ8UpV7vnu/SnL3tcHtz0FBOaz5ZeSDO7ldJQYB5t3uF1Jnr1rW364EafKl6l8w592lkuBP3a5suCVf96+bfYuq5jq8BqZ/QlkGge/R64ET2+Cib4PdqCQrqoRFy6FxvUrYLorQXspRDBd8Ew69G9mSSgWvVwUspaOUVXnOZ2HnCypCv/wH/u2W/VDZVq//IMILBbH5D+oCuPCL6m8hVG6q+xhs+F3057YaQUD1pKHLhYCLjSq2WPmBtgP+5593lwqGjm1Wn3O4i+tZn1Z3cxtjHFuu4BqAoxuhfnHoukv/Tb2Xld/Oyp1L1kXfrNVPmA2/UxH8wi9H367CaPXfGYfo97aqiNP8ooSAZf+hTvoLvxW+LDBRFn1Fvd7aX/mX9Xeo/Y4PEzkDLL0nLdH+b7f9loc/fjh0hSdMpA+qhn3+HSra7wyYIuHwOkCEL/Gz2pVd0LheXWhvfUolPhd8XtklwXcNoES3eZdf9O1FymdueB9+Nk1FgCv/SVk2y34Y/s05ipXgDHaqaHNKgA1ntSv7QljV9xhI9USYc6uyKubePnTdwi8qsZ13e2ipZCDz71RRa7TEd7L0tytfuMSwOy76jrqAXvbv6k7FpKpeFRBse1pdwKLhcavk4oQLYfTZ/uUTLlAX1Pfujx5gtO5TF3ZHcei6GZ9U4x/e++/ogta6X332oER+4ZfUZ70owvlsL4K5n1NjYbpPRH9/ucDRTco2rV8Suq64WrkIh9eqO5wMk3V7B5Loq+/sVbdHU69QtfjRqExA9A8bvmz9Bf5l48+Db+6COTE86XgpH6tuqT/8kzqhARo3AlLdYYSjboHy9tf8b+yEaBRcHhebT27G5Q1KtIXz9E0WfUVFkB8ElI42rFE2VFFl+B3N/ztVA3/L4zDB+NGf/0V11xB4sTNpP6QG94ya6V927t/DnSth5nXw8V9UUnnOrTD6rMhvcN7tcNZn4IofhtpwZ38GvrNPWT7BfOoX8KW1ygoKpKAMvroRrgqTXwjEXqje38G3VQ4gnfScVI+lI9RjzWT49l5lowWz6CsqYo7l7a//P3U+nP/F0HVL/wX6WqNH1G37Q6N8E4tFVTU17YC9r4XfZrBbJacD7aFL/hW+tlVZapFY8Hl1h7n5scjb5AoNawAR3gIFmPM5lWt6P8rd7zCR9UgfkqjVX/8b6GuBC78de9vCSuUTBnudUqr67cBKg4Y1KkcQPCCqbFRqXn4wi76iqlw2/l79fWSdikLHLYj8nCVfU+9524qkdimlxOV10efuY0frjqErvWb1jiP0iVUTlK216VF1snrcqnom0o8ZVCRz61NDbZjyMXDOTSpR2NsydHuzcmdkgOgLoS4Y1/0KvrMXbn5CJSyjYbWrMtnJl0Q+rnBYrKGCb1JYHntkMqho31EaOmApVXqb1GPpKP+ySCNzi6rUndnHz0b29j96XA3Gmv6J8KOnx5+rvrf3f+m31oJp3T80iRvMWTeoOwGz0igYs1wz8MJhscQecVwzWd01bnok9xO6h83AqCr8eqtN6cCRdWoMSLwMdqd8aFkV/XGVhugnUqs/2K1+TFMujxwZByKEsniCPf2jm+AP16gTwOTQeyrTPtxtAEafBZMuUXcrbqeyS0afNTThE4x5K772V0n94D3Sg0Q9L8TXNyP9YE/fZNFyZZtsfgxOfqwGwIXz82Ox+G4V0a8PsqmadgBCjZ4Oh6MEZlytBDhXKapUdxofPzvUCksVs0a/ZGT07UwWfkn9PtY9GLpux/PKKpu0VI0jiWRXXnyPCjDCeft9bcpqDB6tHIjNoQSt4T3jLjaItgg5gXiYcIEa5R5PWWi28LjV2I1Y58jcz6mgNNLFMRwH3krp0CDLol9otyZeq//Br9WPbum/xN7WpKJOlbQFcsK4DV//kDoZ+juUFxrOgxsOFn9VldxtfUpdgCL5+SZCwMKvKO97/xsJ787p8Y+UXH88KLKI5Omb1M1XP+B1D6qEJSQn+iOmqQhz/UOqz5BJ0w7l74bziPOJ87+oBDfQCkuVYHsnFpVnwKzr1Z3ZQKd/+e5X4ZkvqKDm5sejR9VnnK+CkjVhov1w5ZrhmHe7smre/3noOp/oT4z9foIxf6MBv+ec48QWdScfLokbSEGpsjF3vRw6ojwSe/+a8uFlVfQhwVr9gU7lV067Krw3G4nK8aG3uyd3qBLJcfPh+a8aAz+k34MebiZfpuyM1/5dVVHEc9dy1g2qQmXtAwnvzvTxHRYHHzZ9iCswUvKJfhQbY/FydeF877+VF5lsL5RFX1EX7Y+f8S87uWOotZOvVNXDrOuMcQJd6XnN3iYldIUR8ifhWHI3OLuVfej1wFv3qXYNo2bCrSvi60O19B5VNbQxKPEfqVwzmIIyOPcfYOdLoYLWtl/ZVQVl8b8nE1P0g/NSuYQ5ZiOW6IPKzVgd8duCwa1hkiAHRD+BWv21Dyjhv+S7ie2kYrxKmgZWJDTtUKV7n/m9iqJfvUcN1DEHUgw3QhgCaCRzY0X6oG6bz/sHFemf3BF7+wBM0Z8/aj4DngG2tgQkHL0uVVEUrTpp2lXqRO9rhfokonyT+sVq0Nv6h1RU7OpXInAqiD6odhuDXap5XjroaVbWTiI5pTGzYeLF6o7jj59WA4Fm3wx/90rk5HswZyxUr/H+L4belbXuV7+VeDpjnv9FJWjBpbptB5OzdsAfmOSyvdOwRp0rgdVVkSgdqb6bjx6Hnqbo2/a2xleQEoMcEP04a/X729VgrDM/lXhXQ1/ZphHtS6lEf9RMFZ1d+wAgVdRv1mlngrNvVBFPxfj4I+cFnwdbkfosEsBtJGsXjV2EQLD+RIDF43GqC140LAG16uPjuCuJhBDqwnVimxqF3LxbVQeNOkVEf9w8Vf217kF/e4tU6DkZv7UTyJK7lX14eB1c8z9w3YOJdZoFVd7a26zsB5O2/cpCipT4DqR0BMy5RVmYgXc+gTX6iZLr9o7Xq0Q/nijfZPFy9X6Cc13BnIxRihsnWWutbGLW6p/sGmBsZRTB3f2qiqAu+EbiOwks2xw5Q9X59rf7o8szP6lOjKokPMZUsBXADb+NvzMhqAqUObeqKpjLvjd0uHoUTE+/tqiW6dXT2XBiA1+a/SW10uOO7OcHMvdzKol79o3xH284zvksvPZ9Fe1PWaaWjZwV/Tn5xMIvqtG/+15LbX4BUPZOeRJW2uTL4BP/T91BRitzjcaEC6G8TlWMnWN85637Yls7gcy+Vdldu1eqiNbZqy5Gyfj54A9OPE6cbi/feWYLbb3hz5+xFUX856fPxmKJ7y7J5fHyo5d2cPN5Z3DmmCSLBpp3wUBHYrnB2qmqhfTmR9VYk0g2a6zxF3GS9Uh/dLlKKDV3D0bfsHGD8uDHzEl8J+aoXPPWqClMieC829Vw9kwz8SKYenliz1n4ZRUZJJAwNO0du8XOuaPPZUvTFgY9xmfuccY3+YWtQJWOplpF4yhRF5Adz6vJTKwFyUd+uci0q9RgrWA/PBl6mlVTt0QRQiUJkxV8UHd3Z39GNavraVZ3yK0HYidxAxl/nsoBbTNyOOZI9EQuHIH47B03R9r7eP6jYzS299M76B7yr7G9n6c2HuFYZ/xFIs9uauTRtQ28tuNkcscGRiUaqt1CIiz4grqrC7yrCubENjVxT4pkXfQri9WX2N4XI9ptXK+St8mMii0bo8oRzbJN0w8P7nOTL9ROgZnXwPrfxp0w9Im+1c55o8/D6XWypcmYAMXrii/STyfnfkGND9i2QlX1JDLjUq5jtcG8O9TgpOB2G4nga8EQZ7nmcHDOZ9WI4O1/UZ6zszt6uWYwQsDZN6iLe2+Lv2trGuwdl0e1Pf/O30znz19eMuTfj65TF7sjcZaDD7o9/M+bKuHc2hMjAI1GT5hxFfEwdZmyeaMFCie2DR1BnSQ5IPrqS+zsj5KYcfaqATzJJlktVjUS1hfp71RVMJmekjCdLPm6qp3f9Pu4NjerdewWO/NHzcciLH5f3+OK7emnm5rJp6a1YzLvdiV4mx5N/jX625TgJiog6WTULPX9bFvhb4eciOiDGiUtPapNeSrlmhAQ6Ttxe1Qe0BbGvhlvjPY/EufAzxUbjnC0ox+H1UJLBLsoLnqbVYCZSLUVKI2af4ca1R2ufNM1oPJfp4LoVxmi3x7tgz66WSX76hKYzi2YioCyzabt+Z84HDdPDbJZ+6u4WjM4verztVvslDnKOLP6TP8gLY8rvlGn6cZsVZzv30U4KsYpm+fDxxLL2QRiRo3J2Dvp5JzPKnt13+vq70StmVGz1MC7bc8q0S+ujd5uIRrm79Tr9kX6dluojI2pLMQi4mvxMuDy8L+r93HuhCpmj69ILdLvbVLfVzJzJ8+9XV0wwgVyzTvVhfNUEP3yQnVb394XJdJvNMSpLkqbgliYo3K9HnXFPBVKBC/4hvIB45iJJ9DTB5hZM5MDnUbU5XFmR/SnXK6mZ5zzuczvOxOc+3mj+uXF5J7vG5iVRXsHjAmKBHzwkLojNKvh4kUIFe0fXqOa6KWSvxli76hI3x5GYO1WC2Mqimhsix3p/+mDw5zsGuSby6ZTW1pAa08KkX6yORhQ7V7O/JQq0ghub20mcU8F0bdZLZQX2qLbO40bVPIoFTumok61jW3Zo1oBnAqiP/Fi1ff//V/E7F9u2jsO46QpsBb4G69546zeSTcWiypBjdYaO5+ZdClU1vt7LCWKOQdENu0dUOdO/RLl51dPTC7/ctan1WPrvsTtoUAsgfaOivRt1vDVOXVVRTHtnT6nmwff2sfiyTUsmlxDTamD1pTsnabU7swWfF5V/2x/bujy41tVIUsaKgyzLvoAVSWOyIlcKVVDolQHTVWOVxaROW3eqWApCKGi/bYDqsd6FIIjfbvF7qvdV3X6p1AiNVewWNT8CYfeVXeXiZIr9g74SzaTrbqpmazm1IU0RfpuXMbYHntE0S+Omch9/IPDtPQ4+dYVqltvTUkB7X3+C0rC9Lakdmc24UKomRra5fTENlWJlYxtFEROiH5lsSOyvdN+SDV/SlX0zVvSvauI2twr35jxSTU6MsbEF4GePoDNYgsQ/SxU75wuzP1b1bk1XEvpWPScNFowJOl/p5OZ16pBgYET0CSKOY91SqIfGunbreFlbHx1ESe7Bxh0h78Ldnu8/P79Q5w/sZr59cpFqC11IGUMuzkSUqoLdSoXaXPwYuMGf/dNr1c1OkyDtQO5IvpFdjoiRfo+Pz9Not+wRv3oMjnydjixWFXN+6F3/S1rw+Cr3rH6I32P9OCV3uwlck8HSmrVYLotT8YeZh+MOWNWOtt6J0tRFfzjO8kNjjSZfYual2JShNbX8RCmZNMWIfodX1WMlHCsI3yhw8qPT3C0o5+/v9B/EaopVR12W3uTSOYOdqmJnVLNwcy5bWj3zfaDalDkqST6VcV2OiJdWRs3qIkhUvXgzQFaXvepYe0EMvtW1Q/loz9F3MSM6h0WddLYDDvH7XUbdfpa9IeNRV81BtP9X2LPSzVqTDcjpqU2MK+4Gq57ILUcTkDvHV8iN4qnD3AkTDJXSslv3z3ApNoSLpvhF+maEnV+JJXM9bXBTvE7C+6+mcYkLsQp+kKIK4UQu4UQ+4QQ94RZf4YQYrUQ4kMhxFYhxNWJHISydyJ8yEfWq/LEVAfvOIqh2PixnWp14RXj1LD7jx6PmNA12zCYkf4Q0Y+n944meWomw4xPqIndE5nnuLcp+5U7uYavZNOF2xvL3ok8M9/6g21sbezk8xdMHNKmwYz0W5Ip2+xNk+jD0O6bJ7apSZZGpGCtBRBT9IUQVuBXwFXATOAWIURwqPxvwAop5VzgZiChbmCVxXa6B9yhyRNnn/Ky0tX50rR4TrVIH5TF03U0YuvV4ESuKfouryv+3jua5FnyNVWVkUj3zR4t+iEE2jtuY3BWhEh/VHkhdqsIW8Hz2/cOUlVs54Z5dUOW15amEOn7ZjlLw3cW2H3zwGoYMT32zGJxEk+kfx6wT0p5QErpBJ4Erg3aRgLmfV8FcCyRg6iKNCr3+EfKjgk3AXcymI3XToVyzWCmXwVF1RFFJVz1jm95vL13NMkz/jzV/GztA/F13/R6VSVIvDNmnS5YAuydGJG+1SIYW1kUYu8cbOnl9Z0n+duF9RQ5hrZ1KS+0Y7OI5Dz9dFdbmd03j25Km7UD8Yn+OCCwiXOjsSyQHwCfE0I0AiuB5eFeSAhxlxBioxBiY3Nzs2+5v/9OkOibSdxoc8cmQu00KKg4tZp7mdgKVGSw62XVdzuIwN47EM7T15H+sLPkbjURzY7nYm/ra8GgRX8IAZ5+tDYMJuPDzNfx8HsHsVst/O2iCSHbWyyC6hJHkpF+CyDUiON0UDsVphtOeYZFPx5uAR6RUtYBVwOPCSFCXltK+ZCUcoGUcsGIEf6rodl/J6SCp2mnaiubTD/xcFzwTfjHt5Jr2pYPzP2cEvAwk6ebnr5NKLEfau9koffO6ci0q1TAsTmOfjzmaNxcSuTmAmGqd8K1YTAZX11EY4C9I6XklY9PcMXMUYwoCz8Xdk1pAS3J2jvF1em9a77wm6oLbRqncY1H9I8CgeOu64xlgXwBWAEgpVwLFAJxX+6qjEg/pILH2ZPeGuWC0lMzyjcZNUuN0N30iLIHAnB5XdgtdoRR/mfaOyqRq6t3MoLFAmdeo8qGA+evDUey3RpPdcJV70QZsFRXVUxLj5M+p7LUdp3opqVnkIunRb6Y1pY6krd30m3H1S2AfzmqilnSRDyivwGYKoSYKIRwoBK1wcM/DwOXAQghzkSJfjNx4mu6Fhzpu/pPnXr6TGFOnh5kIZiibxJSvaNFPzNMv0rlqczmZZHwtWDQ9s4QLFZVnhxHGwbwl20eNSyed/aoz/XCqZFFvyZpe6dZjctIN2k+N2OKvpTSDXwVWAXsRFXpbBdC3CuEuMbY7FvAPwghtgBPAHdKKWPMf+inIlKk7+oHe3G8L6MB1eNkxAw1GXZA+abL4/L13YEgeydbvXdOR+rOVQn33a9G307bO5GxOsDr8rVhiOrpVw9tsfzu3sNEf2YAACAASURBVBamjSpldEXkSpjqkoLkOm1me+6DOInLfJJSrkQlaAOXfS/g/zuApE2nsgIbNougoz840u/TP/pEsVhh6Xfh6TvUbEWzbwJCI/2h9o7uvZMxLFaYeoVqB+JxR/Z/e5qUl5sLLRhyDatD2Tt4sVuFz7IMh3+AVj/9Tg/rD7Vx+8L6qC9fU+qg1+lhwOWh0J5A/s+cxD7HyYkRuUIIKovtodU72t5JjjOvUdn+t+9Tfj2x7B1dvZNRpl+p5mhuXB95GzNqzIUWDLmGxeazdyK1YDAZUVpAod1CY3sf6w+14XR7uTCKnw8BtfqJdNt09asupMNh76SZnBB9UBU8IdU7rj5t7ySDxQKX/KvqvrnlCUDZO3ZraKTvcg+q0kDt6WeOyZepaqndr0TepuekvsuNhBnpe2RUPx9UQGl223x3TzMOm4XzJkRv0V5TYvTfScTiyaMcTO6IflGY/js60k+eaVfCuPnw9k/B7Yxs75iTo2vRzxyF5TBhCeyJ4uu3HVTtNTShGKLv9npxRBiYFYjZV/+dvc2cN6E6ZEBWMDXJjMr19d3Roh83Ydsr60Ru8gihGn11HoHjW3B6neHtHbfRgVDX6WeWaVepCX3MicID6TqmOiuesSjzx5UPWG2+NgyxIn1QA7T2NvWw52QPF06Nbb/UJtN/pzfNo3GHkZwRfdVpM+DKKqVh7+hIP2nMH6C7P8Te8VXvmKKvPf3MMv1K9Rgu2j/0nnqccEHmjief8FXvxPb0QQ3QcrpVeWe0Uk2TmmQ8fd+4Ci36caMSuQEfssepZrrSop88gaMXI9o7pujr6p2MUjVBdU0M5+sfeldV7Yw6K+OHlRdY7b42DI4oo3FN6qqUW1BbWsCM0WUxty922CiyW5Pz9LW9Ez+VxQ4GXF4GXEZtucsYOq3tneSxGaJvePpmL30IjPRNT19H+hln+pVqdG5wr6RD76lh96dqu5BUsTp8bRii1eibjDdE/8KptUPaKEejpjTBAVq9zVBQnrZOmMNJzoh+la//juHrm7PB60g/eQIifafHGdbe8SVytaefec66QVVObf+zf1nnUVV1pa2dyFjshuhLbHEkcieNKGFibQnXzY0/MV5TWkBLovZOHpRrQg6Jvr/TpvFB+0RfR/pJYzUaSkWxd1y6eid7jD5bTeiz5Un/sob31aMW/chY7eBx4/Z6I86aFUhJgY3V314atd9OMLUljsTtnTywdiAHRd8f6Zv2jo70k8bmj/TdXndYe8ftNi6yWvSzw+yb4OhGNS0eaD8/Hgx7x+2REXvpp0pS9k4eJHEhh0S/Kri9srZ3Use0d9yDKtIPV73j0Z5+Vjn7RkDA1qfU39rPj43VbliW8Xn6yVBTWkBr7yBxtxDLtfmMo5Azoh8ykYpO5KZOsKcftveOcZHVnn52KB8Lky5Woq/9/Piw2sGrplcdtki/xIHLI+kaiGOWM49LTXqj7Z3E8EX6/TrSTxtBJZu2gKZq/kSutneyzjk3Q0cDvPv/1N9a9KNj2jteGZennwzmAK24fP0+o/pK2zuJUWi3Umi3+D19Z6961JF+8tiMRK5ZshnYWlmY9o4W/axz5ifBVgQbH9Z+fjwkWL2TDAkN0PLNjasj/YSpKnbQ3qsj/bRhRvYepxqRG2DhCCGwWWy4jS6c2tPPIgVlSviR2s+PB2NwlssTX/VOMiTUdC2PWjBAjol+RZGdjv7gOn0d6SeNEGAtQLoHQnrvgPL13caE6bqffpY552b1qK2d2JgN1+JorZwsZnvluObK7cmfDpsQ5yQqmaIqsL2yLtlMD1aHz7cPFn2bsOEyE7k60s8uky+F6x5UcyFoomO2VhbDV7JZVZJAp01fCwYd6SdMVUnARCpmpG/Top8SNoevqZojSNjtVruaRAW0p59tLBaYcysUlGb7SHIfs8vmMNo7dquFiiJ7fBOk9xqznBXE7uuTC+RUpF9R5Bg6OMtWqE4GTfJYC3yiHy7S93v6WvQ1eYLRZdMt4mutnCxxD9Dqya9ZznJK9M32ylJKhJ5AJT1Y7b4BWCGib7Hhkkakr+v0NfmC1aHq9HEPm70DUFtSQHM8idym7VAzZdiOI93kVBhdVezA7ZX0DLqNCVRKsn1I+Y+twFeWaQ+K5ofaO9rT1+QJvqq04RX9+ppi9pzsxuuNMirX1Q8nd8C4ecN2HOkmp0S/IrD/jp5AJT1YHZEjfWELEP2cuunTaCLjG3ToGrY2DAALJ9XQ0edi14nuyBud+Fh1Sh2rRT8phrRX1vZOerA6cJqibw1j7+hIX5NvGL9V4XUN2+AsgEWTawBYe6A18kbHNqtHHeknR1Vge2VXn67RTwdWh68sM1ydvvb0NXmHcVfqwI1jGBO5YyuLmFBTzNr9LZE3OroZSkdB2ZhhO450k1OiP6Snvo7004PNgcutKnTCJXLdXmOmMl29o8kXjEjfjntYI31Q0f4HB9pwe7zhNzi2WVk7eVK5Azkn+urL7Ow37R0d6aeM1YHLGHUbXKevqnc8KjGWRz9azWmO8Tu2Cc+wevoAiybX0j3oZvuxrtCVA13QsjevrB3INdEvMiL9Xp3ITRvWApze8JG+3WLHLT3az9fkF8ZdqR13XBOjp8LCSdVABF//+EeAzKskLuSY6NusFgrtFnqdbm3vpAubP9IPa+9Ij/bzNfmF8Xt14B623jsmI8sKmTqylLX7w4j+USOJO3busB5Duskp0QfVYnnA5dGJ3HRhjS76LunVfr4mvzDtHTzDOiLXZNHkGjYcasMV7Osf2wyV9VBSM+zHkE5yrji7yCf6OtJPC1YHLq8HsIb23rHYceeI6LtcLhobGxkYGMj2oeQUhYWF1NXVYbdn/zvKGQLsneHqvRPI4sk1/GFtA1sbO5hfX+1fcfRDqJs/7PtPN3GJvhDiSuAXgBX4rZTyvjDbfBb4ASCBLVLKW5M5oEK7lUGnCzyDOtJPB1aHUZZpjWDv5IboNzY2UlZWxoQJExA6qQyAlJLW1lYaGxuZOHFitg8ndzB+rw4xvCNyTc6fWIMQsGZfq1/0e1ug8zCc9w/Dvv90E/MTE0JYgV8BVwEzgVuEEDODtpkKfBdYIqWcBXw92QMqsFnwOPUEKmnDVuAbgBVW9PHmhKc/MDBATU2NFvwAhBDU1NTou59gAks2M9CQsarEwYzR5UOTuXnq50N8nv55wD4p5QEppRN4Erg2aJt/AH4lpWwHkFI2JXtAhXYr0jdVohb9lLHaVVkmYXrvWOyGp58b1Tta8EPRn0kYjN+xDU9G7B1QFs/GhnZlPYMxElfA2DkZ2X86iUf0xwFHAv5uNJYFMg2YJoR4XwixzrCDQhBC3CWE2CiE2Njc3Bx2Z0V2Kzj1rFlpw1rgF/2wkb7UfXc0+YURpDiGuctmILPHV+J0e2loNSZ3OroZaqflTQ/9QNL1idmAqcBS4BbgN0KIyuCNpJQPSSkXSCkXjBgRfpaZQrsF3NreSRtWB04jGApbp4/MmUg/2/z4xz9m1qxZnHPOOcyZM4cPPvgAgJ///Of09fXFfH682wXS1tbGsmXLmDp1KsuWLaO9vT2pYz+tsPgTuZmo3gGoMWbSau9zgpRw7MO8G5RlEo/oHwXGB/xdZywLpBF4QUrpklIeBPagLgIJU2i3BkyVqCP9lLE5cBkWQVh7B5kTnn62Wbt2LS+99BKbN29m69atvP7664wfr372wyn69913H5dddhl79+7lsssu4777QmokNMEMsXcyE+lX+joAO8HZq2bLGjE9I/tON/F8YhuAqUKIiUIIB3Az8ELQNs+honyEELUou+dAMgdUaLcidKSfPqwOXCjRt4mhNo6q3pE5Ub2TbY4fP05tbS0FBQUA1NbWMnbsWH75y19y7NgxLrnkEi655BIAvvSlL7FgwQJmzZrF97//fYCw2/31r39l0aJFzJs3jxtvvJGenp6Q/T7//PPccccdANxxxx0899xzmXi7+Y2ZyBXuYW/DYGJ2AG7vc+XdnLjBxDRzpZRuIcRXgVWoks2HpZTbhRD3AhullC8Y664QQuwAPMB3pJRR+pFGptBuwWNM76cj/TRgdeAS4LDYQ5KCNosNjwCvsOXUKL3/eHE7O8L1OkmBmWPL+f6nZkVcf8UVV3Dvvfcybdo0Lr/8cm666SYuvvhi7r77bu6//35Wr15NbW0toGyg6upqPB4Pl112GVu3bg3ZrqWlhR/96Ee8/vrrlJSU8JOf/IT777+f733ve0P2e/LkScaMUR0aR48ezcmTJ9P6vk9JrP4RuZmK9P2i74S+DrXwVBV9ACnlSmBl0LLvBfxfAt80/qVEod1Kv0dH+mnDVoBTCOyW0K/a9PjdVjunu6tfWlrKpk2bePfdd1m9ejU33XQT9913H3feeWfItitWrOChhx7C7XZz/PhxduzYwTnnnDNkm3Xr1rFjxw6WLFkCgNPpZNGiRVGPQQihq3XiIaBkM1OiX+SwUmCzqLk+eo1Wy8W1Gdl3usm5so1CuxWbe0AdmRb91LEqT98uQr9qm3EhcFttOSX60SLy4cRqtbJ06VKWLl3K2WefzaOPPhoi+gcPHuRnP/sZGzZsoKqqijvvvDNsHb2UkmXLlvHEE09E3eeoUaM4fvw4Y8aM4fjx44wcOTKdb+nUJMDTz1QiF1S0397rDLB38qv9gkku3dUDUGizYpfa3kkbVgfuCJG+KfquMOtON3bv3s3evXt9f3/00UfU19cDUFZWRne3mjKvq6uLkpISKioqOHnyJK+88orvOYHbLVy4kPfff599+/YB0Nvby549e0L2e8011/Doo48C8Oijj3LttcFDYDQhBDRcc2Qo0geVzG3vc0GfEemfyvZOJim0WyhCzfSkI/00YHXgAhzCGrLKtHe06ENPTw/Lly+no6MDm83GlClTeOihhwC46667uPLKKxk7diyrV69m7ty5zJgxg/Hjx/vsm3DbPfLII9xyyy0MDqrpKn/0ox8xbdq0Ifu95557+OxnP8vvfvc76uvrWbFiRebedL5izXzJJqhIv6PPqewdWxE4SjK273SSc2d7kcNKIeok0ZF+GrA5lKcfRvT99k7outON+fPns2bNmrDrli9fzvLly31/P/LII3Ftd+mll7Jhw4ao+62pqeGNN95I/IBPZyxWvMJqTKKSuUi/qsTO7hPdSvTzNMqHHLV3ioQTKay6lDAdGJ6+LZro60hfk2d4hT1jXTZNKosdKpHb15K3fj7koOgX2C0UMYi0Fekp/NKBtcBI5IZ+1T57J8w6jSaX8QpbRks2AaqK7XT0u5C9LXlbuQM5KPqFditFOPHYtJ+fFqx2VacfNdLX9o4mv/BYbFmp3vF4JbKnWds76aTQbqVQDOK1Fmb7UE4NbAU4EdgJPTm06GvyFY9p72TQ068sdgAS+lq1vZNOioxI361FPz1Es3eEWb2Tcz8DjSYqHmGjQLixZKgNAyh7p5hBLJ4Bbe+kk0LD09einyasdiX64SJ9I2fi1p6+Js/wCBsO4cnoPiuLHdSITvWHtnfSR6FdVe+4LVr004KtQHn6YUTfLtWjW0f6QHZaKz/99NPMmjULi8XCxo0bkzru0xE39oyLflWxnRrU4DtKdKSfNgptqk7fqUU/PZhtGMKJvvGoq3ey11r5rLPO4s9//jMXXXRRUsd9uuIWNhzCndF9VgVG+treSR+FDjUi1ykKsn0opwZGa+VwIx5saHvHJFutlc8880ymT8/PvuzZxIOVggxH+uVFdmpE/kf6OTcqp9BupZhBBrTopwejtXJY0ZfK33Hn2nCIV+6BE9vS+5qjz4arIk9Qkq3WyprkcAk7djIr+laLYKyjFyR5Lfo5F+IV2lTJ5oDQ9k5asNpVGwYZukrbO37M1soPPfQQI0aM4KabborYbmHFihXMmzePuXPnsn37dnbs2BGyTWBr5Tlz5vDoo4/S0NAwzO/i9MGNDXuG7R2AsbYeBkVh3vbdgRyM9O1WQRFOBtGRfloQQnn6MlT1zUjflWuRfpSIfDjJRmtlTXK4seIg86I/wtpNl7eC/K3dycFIXwBFDNKfUx3e85uIou/1AuDW7S6y1lpZkxwubNgz7OkD1Ipu2inP+H7TSc5F+nicWIWkX0f6aUFKiUsIHGFE32fvZPaQcpJstVb+y1/+wvLly2lubuYTn/gEc+bMYdWqVZl743mKCxv2LET6lbKTI7KcabE3zVlyT/RdquStT+pIPx24verEsEtvyDqbscwdppzzdCNbrZWvv/56rr/++sQP+DTHhRVbFkS/3NvJSc+YjO83neScvYNLzY/b69Winw5cXhXHh7d3zOqdMFlejSaHcUkbdplh0ZeSEncHTZ5SBt2Zt5bSRc6Kvo7004NP9L2hP1Iz+tf2jibfcKK6bGZ2pz3YvIO0ynLVVz9PyUHRV/ZOr1dPoJIOTNEP5+n77R2NJr9wyizYO71qbtw2ymnvc2Z232kkB0VfRfrdHh3ppwOnR/047Z7QqMhmRP9utL2jyS+c0opNZjja7msFoEWW096rI/30YUT63Z7cyzHnI2akb5Ohoi+8bmxS4tKir8kznDIL9k5vMwBtslxNkJ6n5KDoq0i/y6PtnXTg8piefphbYY8Tu5Q60tfkHU5pxZrpSH+IvaMj/fRhRPpdbi366cDn6XvCiL7XjU2CO4zffzqSjdbK3/nOd5gxYwbnnHMO119/PR0dHUkd++nGIFaseCFMgcKwYUT6rbJMe/ppRUf6acXpjezp43FiR+IitIb/dCNbrZWXLVvGxx9/zNatW5k2bRr/+Z//mdTxn244vcYUn54MRtx9rWAvBnuJtnfSiiH6nTrSTwt+eyfMyeFxYZMSd5iBW6cb2WqtfMUVV2CzqfzVwoULaWxszMTbzXsGpJHzC/e7Hi56W6C4lqpie17bO7mXLTXsnXaXnqw7Hfjq9MPZOx6XsndyzNP/yfqfsKttV1pfc0b1DP75vH+OuD4XWis//PDD3HTTTWl936cqTq8Rr2Yy0u9thpJaKi0OHemnFd+IXDsuj45AU8Xv6Yc5ObwuZe+Eqew53ch2a+Uf//jH2Gw2brvttnS9pVOaAZ+9k0Hx7WuBklqqSnSkn15cfbgtBUgsDLg82K25d13KJ3z2TjjR9zhz0t6JFpEPJ9lqrfzII4/w0ksv8cYbbyB0x9O4GJBZ8PR7W2HkLCqFg53HuzK33zQTl6IKIa4UQuwWQuwTQtwTZbsbhBBSCLEg6SNy9eOxqglUBly5JUb5iN/eCSf6qnrHlWOinw2y1Vr51Vdf5ac//SkvvPACxcXFw/b+TjUyHulL6bN3qorttPfmr70TM9IXQliBXwHLgEZggxDiBSnljqDtyoCvAR+kdESuvgDR17ZDqkQXfSd2wJXpeuccJFutlb/61a8yODjIsmXLAHWx+PWvf52hd52feLwSV6YjfWcPeAaV6OOgs9+F1yuxWPLvziwee+c8YJ+U8gCAEOJJ4Fog2Mj8IfAT4DspHZGrH6+tCNCinw58nr5rUEUrgfaB14UNf/vl05lstVY27wQ08ePyeHGZ0pWpSN8YmEXJCCqlA6+ErgEXlcX51y4mHntnHHAk4O9GY5kPIcQ8YLyU8uVoLySEuEsIsVEIsbG5uTn8RkNEX9sOqWL23rFJLwSLu8eFDaFFX5NXuL0SJxku2TRF3yjZBPI2mZtyllQIYQHuB74Va1sp5UNSygVSygUjRkSYZdLVB6bo53HP6lzBZ+8gQ6Mijws7wreNRpMPuNxe3GTY3ukzI/0aqozoPl9H5cYj+keB8QF/1xnLTMqAs4C3hBCHgIXAC0knc139YNf2TrrwT6ICuAeHrvQ4cyrSl7odRAj6MwnF5c2CvdO8Wz1WjKfSiPTztVY/HtHfAEwVQkwUQjiAm4EXzJVSyk4pZa2UcoKUcgKwDrhGSrkxqSNy9YFDVTH0O7Xop4pZsumQMjQq8rqxCUtOiH5hYSGtra1a5AKQUtLa2kphYWG2DyWncHskTnNEbqYi/f1vwsiZUDrSH+nnaXvlmIlcKaVbCPFVYBVgBR6WUm4XQtwLbJRSvhD9FRLE1Y8oVaI/4Naefqo4vU4E6ovDExrp27HkhL1TV1dHY2MjEXM9pymFhYXU1dVl+zByiqGJ3Az8dp19cHgtnHcXQN7bO3ENzpJSrgRWBi0LO55cSrk0pSNy9WMxIn1t76SOy+vCLqxq6vPgE8TjwiZyw96x2+1MnDgx24ehyQNcHhng6WdAeBvWqP1MvhSAskIbJQ4rB1t6h3/fw0DuDXd19WEpUKI/qEU/ZVweF3ZhXNtDPH0XNmHNCdHXaOLF7fVmtnpn/5tgLYD6xQBYLILzJlaz9kDr8O97GMg90Xf2YTU9fS36KePyunBYIiS9vC7s5Ianr9HEi9sjM2vv7H9TCb5RYAKwZEotB5p7OdEZ2oIj18kt0XcNgGcQa3EloOv004HL68IeSfQ9TmzCmhOevkYTL06PF7fMkL3TdQyad/qsHZNFk2sAWLO/ZXj3Pwzkluj3twFgLanBZhHa008DQ+ydENF3Y7foSF+TXwyN9IdZ9PevVo9Bon/m6HIqi+2s2Z9/Fk9uib4x2zzF1RTZrTrSTwNDIv0wdfp27elr8gy3J8DTDzdPRDrZ/yaUjIRRs4YstlgEiybVsHZ//pUZ55joq0if4hoK7FY9IjcNuLwu7FZjFrIwnr62dzT5htPjzUz1jtcLB1arKD9My+vFU2o52tFPQ2tiU2Rmm9wSfcPeoaiaQruFAT04K2WcHid2i9EUKkwbBpuw4ZGevItWNKcvGbN3TmxV7kOQtWOy2Ofr55fFk1uiH2DvFOpIPy0oe8eI9N1hRN+iIiZt8WjyBbfXi8uM9Ifzd7v/TfU4aWnY1ZNqSxhVXpB3ydwcE/129VikPf10EdXe8Th9fr+2eDT5gssjkViQFtvwRvpHPoARM6BsVNjVQggWT67NO18/t0S/vw0cZWBzKHtHV++kjMvjwmEtUH8Et2HwurEJLfqa/MKcO3vYRb/7OFSeEXWTxZNraO11svtk9/AdR5rJLdHva4XiKgAK7VY9OCsNqEjf9PSD2zA4sRmRvrZ3NPmC22NE1RbH8A7O6mlWlTtR8NXr78sfXz/HRL8NitWHWGDT9k46UKJvRPph2jCY1o+O9DX5gstr6ILVPnyiHzAnbjTqqoqprynOK18/x0S/FYqqASi0W3TvnTQwRPTDVe8YSV4d6WvyBZc7UPSHyd4Z6FB9fUqjR/oAs+sq2XVC2zvJ0e+P9FUiV4t+qjg9Tuy2CKLvdWl7R5N3uL2GvWMdRnunx2jxHcPeAagucdCRR1Mn5pbo97VDsRnpW3U//TTgK9m0OiJU7zh822myy4neE3lVBZItXJ4A0R+u321vk3osjTCtawBVxQ56Bt0480Svckf0PS4Y7Bxi7+iZs1JniOgH1ul7PSC92Kza3skF9nfs58pnr2RVw6psH0rO4zaqd8Rw2js9huiXxCH6Jcb0if35MalK7oh+v1GjPyTS1yNFU0WVbDqMSD8gkWvcFtu16OcEKw+uxCM9rDu2LtuHkvOYJZvDmsjtNSdCj23vmDNp5YvFkzui7+u74xd9KVWfDU3yOL1OFenbCoZGRcZtsU3bO1lHSsmqQyrC39y0OctHk/u4vBKH1YIIZ1mmi94mEBafHkXDP2eujvQTw2zBUOQXfdA99VNBSonb61bRvNU+1N7xRfoqyasj/eyxq20XDV0NTKyYyMHOg7QNtGX7kHIat8eLzSqGOZHbBMW1YLQpiUZlsbpbbteRfoL0+ztsgvL0Qc+TmwqmkCtPPyjSN04W7elnn1WHVmEVVr4x7xsAfNj0YZaPKLdxeSQ2iwCrbRjtnea4/HyAqhLT3tGRfmIENFsDKLSZkb4W/WQxLRuHxRFavWP834z0tb2THaSUvHroVRaOWciScUtwWBx8eFKLfjRcHi8OmyV8RVrYJ/TDw1fBkfXx76S3Oa7KHYAqHeknSZ+/rTJoeycdOH3CbgebI7ynr+2drLK9dTtHe47yNxP+BofVwVm1Z2lfPwZuj8RmscRfstmyBw6v8XfNjIeepriSuKDGFDlsFh3pJ0x/G9iKwJgUvcih7Z1UMaN3f8lmaPWOzejL45L5EaWcarx68FVsFhuXnqF6ts8fNZ+drTvpc+XXxByZxOU1PH1LnPZO+yHjsSH+nfQ2xzUaF1S3zapiO206kZsgfW1DMuXa3kmdENEPPEHMRK4xWtc1nI2rNGGRUrKqYRWLxy6moqACgLkj5+KWbra1bMvy0eUuLo/Ebk3A3jFFvyNO0R/sAVdfzL47gVQVO7S9kzBBol9g2Du602by+ETfag9Tp69OFpu1END2TjbY0ryFE70nuHLClb5ls0fORiC0xRMFt8eL3Ve9E8fvNtFIvzf+FgwmVcUObe8kTH+bz8+HwOod7ekni8/TN+v0h4zIVSeLz9OXWvQzzesNr2O32Fk6fqlvWbmjnGlV09h8Uot+JFw+Tz/OEbmm6HcfC509Lhym6Mdp74AalduuRT9B+lqHRPpFRqQ/qKdMTJpQeydM9Y5NRfra3sk8bze+zbmjz6XMUTZk+dyRc9nSvEXffUXA7TUj/QREX1hBeqGrMfb2CbRgMKkszp+mazkk+v4OmxBYvaNFP1lCRT+MvWPT9k42ONJ1hENdh7io7qKQdfNHzaff3c/u9t1ZOLLcx+XxBnj6MYTW64GOIzBunvo7HovHZ+/EL/pVxXY6+l150TYmN0Tf61G9d4bYO4anr5uuJY0ZvTusDqNkMzCRa9g79iJA2zuZ5p2j7wBw0bhQ0Z8zcg6Atngi4PJIY0SuPXbJZtcxtc3Ei9Xf8SRzkxJ9Bx6vpGsg98+j3BD9gU5ABkX6hqefJ+1Kc5FBI7IPX7IZNDhL2zsZ5e0jbzOxYiLjy8eHrBtdMppxpeP4qOmjLBxZ7uM20H5omgAAIABJREFUI31LHPaO6efXL1IlnvFE+j1NUFipAqU4qSzOn1G5cYm+EOJKIcRuIcQ+IcQ9YdZ/UwixQwixVQjxhhCiPqGjCGq2BrpkM1WklDy6/VGKbcXUl9cbbRgChN0cnGXTkX6m6XX1svHkRi6uuzjiNnVldZzoPZHBo8of/G0YHMqn90bRCFP0qydDRV2ckX5TQklc8I/KzYda/ZiiL4SwAr8CrgJmArcIIWYGbfYhsEBKeQ7wDPDThI4iqNkagMUicNgsunonSZ7b9xxrj6/lm/O/SU1RjZH0Ch2cJWwObMKmPf0Msu7YOlxeV1g/36SmsIbWgfyZbDuT+D19JbRRfX0ziVtRB5X1cXr6LQlZOxAY6ef+HXM8kf55wD4p5QEppRN4Erg2cAMp5WoppTmEcB1Ql9BR9IdG+gCFNouO9JPgZO9J/mvDf7Fg1AJunH6jWhjcWtk8Uax27Fa7tncyyDtH36HMXubz7sNRU1Sju21GwO2VQaIfJbpuP6QE32qHqvr4Iv2epoRFv9poupYPZZvxiP444EjA343Gskh8AXgl3AohxF1CiI1CiI3Nzc3+FWHsHTAmUtGinxBSSn647oe4vC7+Y/F/YBHGV2zeCpuDWcwTxWJXkb62dzKCV3p5t/FdFo9brHItEagprKHf3a/bMYRhSGtliB3pV01Q/6+sV0laZ2/0HaRg7+TDqNy0JnKFEJ8DFgD/FW69lPIhKeUCKeWCESMCrqRh7B3Qop8Mrxx8hbcb32b53OWcUX6Gf4XvBDHE3rRzrA5sFm3vZIqdbTtp7m+Oau2AivQBWvu1xROMvw1DnJG+KfrmY8fhyNu7B1VhSQKjcQHKC+1YxKmTyD0KBJYY1BnLhiCEuBz4V+AaKeVg8Pqo9LepTHzB0EEqRXar9vQT5Pfbf8+0qmncduZtQ1cYPXZ8vr55olht2C123Vo5Q7zT+A4CwQXjLoi6XU2hIfra1w/BFdiGASKXbQ72QF/L0Egfovv65jSJcbZVNrFYBBVF+TEqNx7R3wBMFUJMFEI4gJuBFwI3EELMBf4PJfhNCR+FORpXiCGLC+0WBtI4Irezz8XO411pe71cY0/7Hna17eLTUz+NNXjGHzMqMoeh+zx9HelnCiklqw+v5uwRZ1NdGH0aPh3pR8btNdowWGIkck3/3hfp1w9dHo7exEfjmgxH07WDnQe55rlrONh5MG2vGVP0pZRu4KvAKmAnsEJKuV0Ica8Q4hpjs/8CSoGnhRAfCSFeiPBy4elrC7F2QDVdS+fgrAfe2seNv16L15v7o+aS4aX9L2ETNq6aeFXoSqsZ6QeJvsWOzWLTkX4G+LDpQ3a27eRTkz4Vc1sd6UfG5fb6B2dBZHvHLNc0Rb9khGrfHi3S70m82ZpJZbE97fbO7z/+PQc7D7Lu+Lq0vaYtno2klCuBlUHLvhfw/8tTOor+9iEDs0wK7VY6+9MnRodae+kZdNPcM8io8sK0vW4u4PF6ePnAyywZtyR8FBni6QdU71jsOtLPAI9sf4TKgkqunXJtzG2rjSBIR/qhuLxeHGYbBogc6QeLvhBQeUZ8kX6C9g6oSP9450DCz4tEc18zLx14CYDdbelryZEbI3L7WqG4KmRxoc3CYBoTucc61BfS2H7qVUSsP7Gepv4mPjU5QhRpCxJ9j1ONUBRCR/oZ4GDnQd468hY3z7iZImNAXDTsFjsVBRU60g+DO7ANA0QX/YJyKArQllhlm0m0YDCpLHak1dN/fNfjuL1uJpRPYFfbrrS9bo6IflvYSL/IYU1rP/3jnf0ANLb3p+01c4WXDrxEmb1sSJveIZhRkdmKwePyLdOe/vDzhx1/wG6xc/P0m+N+Tk1hjY70g5BS+j39eOydqvqhucLKemiPUr3T0wz2EnCUJHxs1Wlsr9zn6uOp3U9xef3lXFR3Efs69qXtHM2+6EsZ0kvfZERpASc6B9LiwQ+4PLT0qC/kSNupFen3ufp4reE1rphwBQWmdx+Mz9N3+R+NRJiu3hleWvtbeWHfC1wz5RpfgjYeaor0qNxgXB6lBb6J0SFy9U5guaZJVT0MdipLORy9TUlZO6Ai/QGXNy1l5n/Z9xe6nd3cMesOZlTPYNAzSENXAtM9RiH7oj/YpWrGi0NFv762hEG3l5PdqftkJwK8tlMt0n/zyJv0u/v55KRPRt7IFxUZkX5vExSWAzrSH26e3P0kTq+T22fentDzTsVIv3OwM6Xnu72qhNvXewfCR/per0rYBot+rLLNBCZED6aqOD2jct1eN4/teIy5I+cye8RspldPB0ibxZN90feNxg2NgCbWqFusQy2pR+bHDGvHIk490X9x/4uMLRnLvFHzIm9kC6jekRIOvQf1i9UqLfrDRr+7nyd3PcnS8UuZWDExoefmWqT/wEcPsKV5S9LP39K8hYueuohtzcnP/2tG+jarReWkILyn33NCBTjhIn2I7Ov3tiQ8Gtf30uao3N7U7ppfP/w6R3uOcsesOwCYWDERu8WetvkVckf0w9g79TXFADS0xhg2HQdmEvescRUcOYUSuQ1dDaw9tpZPTf6Uv+VCOHyevhOadqiEldFjXNs7w8eqQ6voGOxIOMoHqC6sptfVy4A7fRUhybKnfQ8PbnmQv+z9S9Kv8fy+5/FKLztadyT9Gi6PivSHDM4KF+kHV+6YxIr0e5sSmhB9yEunqb3yk7uepK60jqV1SwF1fk6pnJK2Cp64SjaHlQjN1gDGVhZhtwoOtaYu0sc7VHS/oL6ax9YdwuOVWC0ixrNyn0e2P6IShDNiJAgDT5ADb6v/T1KiryP94ePZPc8yoXwCC0YtSPi5Zq1+20AbY0vHpvvQEuKVg6qdVrK+ssvj4q8NfwXgcHeURGoM3Eakbw8s2Vz/G9j716EbmqJeFXR3VVQJhRUq0h/sVpOsdB01Ho+pSsJk7Z2S1PvvHOw8yKaTm/javK8NGWA5vXo67za+m/TrBpJ90e82eoaHKZGyWgTjq4vTE+l39lNb6mDKyFJcHklT9wBjKmKXzuUyzX3NPL/vea6bch21RTGik0B75+Db/v7iZFf0T/SeYGTxyOh3KXnKvvZ9fNT8Ed+a/y2ESDzACByVm03Rl1Ky8oAapnO4KznBXnt8LZ2DnViFNenXAH+kb7MIKB8Lo8+BtgPqXzDjz1d1+cFU1sOG38GG34auKx0FE5YkdWzp8PT/vPfP2ISN66ZcN2T59KrpPLfvOVr6W5J+bZPsi37HYRAWnwAFM6GmhIMt6bF3xlQUUVelhP5IW3/ei/4fd/4Rj/Rw56w7Y29sJnKdvXDofTjnRt+qbNk7BzoPcMPzN3Dvknsjjy/IY/6878/YLDaumXJN7I3DkCujcrc0b+FY7zEmVUziQOcB+lx9FNuLE3qNlQdXUu4oZ+7IualF+t6ASL+gFL6YRPS79Ltw4C110aioU4/lY6FsjD84SoJKn6efnOi7PC5e2P8CF4+/OCSIS2cyN/vhVcdhKB/nF6Ug6muKaWjtS3nC4WMd/YytLGR8tfqx5vsArW5nNyt2r2BZ/bKh3TQjYZZsHvkAnN3+OUPJXqT/9O6ncUs321qST+zlKk6Pkxf3v8il4y+N2WcnErnSf+flAy9TYC3wBReJinafq483D7/JFROuYFLlJI50H8ETbbarKPg9/RSka8bVcPVP4YKvw9mfUQUNVRNSEnyAApuVYoc1aXvnzSNv0jbQxg1TbwhZN61qGpCekbm5IfrhbsEMJtaW0O/y0NydWOPOQKSUHOtQkf3YStV+Id8reJ7e8zQ9rh4+f9bn43uC6X/ufQ0QMNHf2jcbbRj63f08v/95APZ17MvovjPBG4ffoGOwgxumhZ7A8eIT/SxG+i6v8uIvrruYM2vOBBL39d9ufJt+dz9XT7yaM8rOwOV1cbLvZHLHY9o71tzMx1UVO5JO5D6751nGlIxh8djFIesqCioYWzL29BD9erNsM4VkbteAm16nh3GVRRTYrIwqL8jrAVqDnkEe2/EYC8csZGZN8MyVETDbMPQ2weizhyTOsxHpv3rwVbqd3UyumMz+jv0Z3XcmeHbPs4wrHcfCMQuTfo0CawGl9tKsRvofHP+AtoE2rp6kBBsS9/VXHlzJyOKRzBs5T83XTPIJYX8iN0dFP8lRuY3djaw9vpbrp1wf2iHXYHr19LSUbWZX9N1O6D4WVfQnGGWbh1JI5prtF8YYUf74quKcifS90puwdbXywEpa+lvij/LBH+mDr2rHJBFP3yvTM7/Bit0rmFwxmRum3UDbQFvWLYx0cqTrCB+c+IDrp1yfcoI6lVr9dHxXKw+spMxexoXjLqTYXszIopEJCXbnYCfvHX2PKydcidViZXyZmprjSPeRGM8Mjz+Rm/14NRzJtlf+y76/IBAhCdxApldP51DXoRSOTpHdT67rqJrCL4roj6sswmYRHEohmXvMKNccW6kSt3VVRSG1+g9+9KCvQiGTfPvtb7P8zeVxby+l5E87/8TUqqmJRZGBoj9x6ZBV8Ub6e9r3cP6fzufDpg/j328Ytrdu5+PWj7lx+o1MqZwCnFoWzzN7n8EiLFFP4HhJdlRu20Abl6y4hJcPvJz0vvvd/bxx+A2WTViGw/j9nFF+RkKi/3rD67i9bq6edDUAI4tHUmgtTKH0MyCRm4NUJmHvuLwuntv7HEvGLWFM6ZiI282ompGWC3l2Pzlz2rIoom+zWqirKqIhBXvHHJg1tsIU/WKOdw7gNqKGhq4GHtjyAI/ueDTpfSTDgc4DvNbwGmuOrWHQE1/OYtPJTexu381tM25LrAxQCNVrx2KH+kVDVsXr6f9h+x8Y8AykXC+8YvcKimxFXDP5GqZWTQVOHdHvd/fzzJ5nuOyMyxhVMirl10s20n9mzzO0DbTxTuM7Se/7r4f+Sp+7b8j8DPXl9XEncqWUPLX7KSZVTGJmtbIhLcJCXVld0hU8ZhuGnLV3iu0JR/qvHXqNpv4mbplxS9TtzAqeVMl50QeYUFvCwdYuVh1alZT3fKyjH5tFMKJMZefHVxfh8UpOdKmLwRO7ngBUOVSvK/Xy0Hh5fOfjgLrSxztK8fFdj1NRUOGLnBLCVgDjzwvpIGizqInRo9lMLf0trDyo7oQ2N21OfN8GXc4uXjn4CldNvIoyRxk1hTVUFlSyt31v0q+ZS7y4/0W6nF187szPpeX1qgurE470XR4XT+56EsDXNmHXiS62HOmI+zU8Xg+/2fYbpldN5/zR5/uW15fX0zbQRrezO+ZrmJPG3Hbm0AClvrw+6Vp9d2AbhhyksthB14ALTwJNIv+080/Ul9fHnEJzXOk4Su2lqR5iDoi+sKiSzShMqCnhsHM133772/x2W5gBFTE43jnAqPJC3wjcuiqVJzjS1k+vq5fn9j3HuNJxeKU3pd4iidA52KlqcuuUvx6PZXK85zhvHH6DG6beEFdP9hBmXQcLQvMANqOHSbQL6tO7n8bldXHhuAv5uOVjXJF6mMfghX0v0O/u57PTPwuAEIIplVNOiUjfK738cecfmVkzk7kj56blNWuKauhydiX0ea9qWEVzfzNLxi7haM9Rmvqa+KdntnLXYxvj7lj7yqFXaOhq4IuzvzhEsM3y4HhE+7Edj1HuKA8Zg3FG2RlJl206Pbkf6UtJ3JM/bWnewtaWrdx25m0x8z9CCOaMnJPyMWZf9KPU6JvU1xTjLVPThf3flv9LeIDC0Y5+xlX6RdIcoNXY3sfz+56n19XLDxb/AIuwsPn/b++8w6Oq0sf/OVMz6Y0ESCEhAQSSSFHEshYURRQQlbWBioCCZXVt237ruvts/e266+6qgAICIqDSxMJK70UIgYSEhFACSUgvk2Qyfc73j5kMCckkoUjU3M/zzJO5955777knZ977nve873vKLl6LvRBW5a/C7DDzwtAXSAhOIKOsY6G/PM+tvV1ITvYWTHjX7Zd8HlpPimVfk7k2p43lecu5OfZmJvabiNVpJaf6wvOnOFwOlhxdwpAeQxgcMdi7Pzk0mRO1Jy45FqOr2XN2D6eMp5g8cPJFReC2xYUGaEkpWZKzhITgBGYNmQXArsJ0soqNlNVZST/jI6VwM5wuJ3MPz6VfWD9GxY9qcaxPkNv7pqMJxaL6IjYXbmZS/0mtFJT44It323R8z236TVG51abOmWuX5CwhSBvEhKSOV1MDuL7X9R0X6oCuF/odmHYA1H7FqP3OMjFhKqF+ofxm528uSPMpMZq9njsAvUIMqAQUVptYlruM1MhURvYayYCwAZc8SdkZHC4HS3OXMqLnCAaED2BI1BAOVRxqd5LG7DCzMn8lt8ff3u5kz8XQpOn7EvrrTq2j2lLN5IGTvRrsofJDF3yfDac3UNxQzNSUqS32J4cm02BvuGjf7e8LS44uIdIQyZiEMZftmk2++tWW6k6VP1xxmOyqbCYPnMyg8EHo1XrWn9hL0/v0q8ySDq/xTcE3FNQVMDNtZivtMy44DoHoUNNfmrsUFao2c0J5XT8vwq7fIrXyRbK3ZC+7i3df9PntMai3O135tmMdp0soNZWy4fQG7u93f6cjnNvy4b9QfhBCP6tuPdKloZ/hbt4Y+QbHao4xN3Nup27hcklKjRav5w64F2DoGezHoapvKagr4NGBjwIwPHo4mRWZF2266Cybzmyi1FTqtfsOixpGrbWWAmOBz3O+Pvk1RquRR6969LLXpz3zjpSSJUeXkByazMheI4k0RBIXFHfBIyIpJR8e+ZCE4IRWq3slh7k9eHzZ9Z0uJ6eMp77XSeFOGk+ys3gnDw14CG0HI9cLwavpd9Ku/1HORwTpghiXNA6tWsvgiMEcqcok2E/DHQOjWHekpF0Tj0u6mJs5l+TQZO7o03rpa71aT6+AXpyu9+1902BrYFX+KkYnjKZnQM9Wxy/ERHQ+l8N75097/8QvdvziO8le2j86iKtjQ/h0f2GHI9flucuRSB4Z2P4EbnOSQpMutYpdKPSl7NBHH9wh3DtKNuCsT6O8VsVt8bcxPmk887LmkV2V3eFtKhus2J2S3iEtF0KPDfPnmHkdkYZI7upzFwDDoodhcVpamC6klGwt3Npp75rOsCRnCXFBcdwc646KbdKefY0ynC4nC7MXclX4VQyPHn7Z6tFEk3mnLaGaXpZObnVui8m4oVFDOVRxyGenrrfV8/r21/m25Fvvvn2l+zhafZQnBj/RSntsy22z0d7I/Kz5zNw4k5uW38T4NeOZnzW/3eeos9Vx+2e3eyfmm2Nz2th8ZnObdS6qL/I5l3Ow7CBlpo5HIEuPLkWn0jGp/6QOy14I50flSinZXrSdRntrb7azDWfZeGYjD/Z/0Ks5Xt3jamqdp7g+KZhxV/fu0MSz/vR6ThpP8kzaMz5tzPHB8e0K7DXH12Cym5gycEqbx6P8o9Cr9S2uUWmubNFffHGpaRgqzZUU1BVQa631Zg693Dx0bTx5ZfUcamfi3OwwsyJ/BaPiRhET2P6cZnMuh9mw64S+09ahjz64O6HJ3kCo8yZvVO7r175OmF8Ybx14q8PbFJ/no99EWGgtJvURJvWf5NXMvMK3mX1905lNvLD5BeYe7tzIoiPSy9I5VHGIxwY+5o286xPch3C/cJ9eMRvObKCgroDpqdMvm624Oe2Zd+YfmU+oPpR7+t7j3Tc0aijVlmqfvtbLcpex7tQ6Zm6cyfoCd8rbhUcWEuEX0WZitRB9CFGGqBZC/1/p/+Ltg29T2lDK2MSxpEWmsTR3KTZf66HiFrzljeXMy5zXarS24MgCXtzyYqsfupSSV7e9yvRvprfKYFjSUMK09dP49c5f+7wnuAXJ58c/556+91zQcoid4XxNf0/JHp7b9Bz/yfhPq7KLcxajQsUjA85pjtH6q0A4SYo1cvvAaHQaFUszvuWVra9QZ6trcb5Luph7eC59Q/oyus9on3XqE9yHgrqCNl+gdpfdO2+T2iO1zfNVQkVcUFyL0cLvdv+OGRtmcMp4qp3WgNMNOeijvuRiFf30snQAgrRBLMtd9p3MI427uhcGrZpPD/gOQFuVvwqj1cjkQZfHy+tC6FqhDx0K/VX5q0gITiApONWbYjlEH8K0lGnsL93PgdID7Z5f4lkm8fyMmlXqjUiXmonJ5yY2Iw2R9AnuQ3q5u2O4pIs5h+cAbhvlpS71BjD78Gwi/CK4v9/93n1CCIb0GNKmndwlXbyf+T6JIYncEd96uH058GXeyarIYmfxTp4Y/ESLybj2RiZmh5mPj37MiJ4jGBwxmFe3vcpbB95i19ldPDbwMZ9r+CaFJnnNO2WmMlbmr+SBfg+w5r41/Pb63/L80OeptlR73UbPp8HWwEc5HxETGEO5uZx1BetaHQOYmzm3hdfIjuIdZFdlY3FaWJTdMk5jXtY8HC4H35Z+6xUWbfHhkQ+xuWxMS53ms8zF4q/1x6AxUGWpQkrp7Y8rjq2gorHCW66isYIVx1YwLmlcizmf2mr3d13gGQL1Gm7t34PN5fNYf3p9q5HTxtMbOV57nGfSnvGZCgDcNvl6Wz211taa7Or81RQ3FDM9dXq7zxUfFE9hnVsoHqk8wvai7bikiw8yP/B5jsPlYF3pv9FF7KSoof2Xgy/Sy9IxaAy8MOwFjlYf5VDFhc9NtcdJ40l+vftlRg3W8cXhEhptrUfPVqeVBVkLGB49/DsZuXfE91ron6g9QUZ5Bg/0e4DEiEBOVZq8b+YH+j9AhF9Eh7b9pmjc5t47VeYqjpu3YDcOw25t6fc6LGoYGeUZuKSLrYVbyavJ48nBT2Kym7yC42LJKM9gX8k+pqZMbeXR0JRy9nxtc2vhVvJr8pmROqPdH+Kl4Mu8MydzDiH6kFZBI4khiYToQ9oU+qvzV1NtqebZIc/y/p3v85PYn7AweyEGjcHrptkWyWHJnDSexOlysuDIAqSULQTHyF4jSQ5NZknOkja1s+V5y6mz1fGPW/5Bv7B+LMxe6C23NHcpdbY6pqZM9QbEgVvLn31oNjGBMYxJGMMneZ94NepSUymrjq9ifNJ4IvwimH14dpv1rjRX8mnep9yTeI83r8zlpslXf1/pPjLKM3h80OM4XA4+zP7QW+bD7A9xuBzMSJ3R4tz0Uw5UzkgKGtwmy8F9q3D55RGsDWNJzhJKGtwTuy7pYk7mHBJDErkr4a526+Mrf06jvZHZh2czLGqY13Tpi/hgt9tmk2IVrAvmwf4P8tWpr3yOIFflr6LaXgTA/vI97V7fF+ll6QzpMYQJSRMI0gV5Y2UuF+9kvMPWwq3UBCyjwWpvc+J8df5qys3lzLp61mW9d2fpWqHfgY/+imMr0Kg0jEsaR58If+otDm+0m0FjYGrKVPaW7G3Xk+RsrQV/nZpgw7mlA5bnLccp7dirf8LB8+ybQ6OGYrQaOVl7kjmH5xAXFMeLw15kdJ/RfHz040vS9mcfmk24X3ibdt+h0a29YqSUzM2cS1xQXIuoyMtNW+ad7Mpsthdt5/FBjxOgbRnMpRIqhvQY0kro2112FmYvZGjUUIZHD8egMfD2bW8zdfBUXrv2NUL0IT7r0C+0H1anlYzyDK/GGht0bo0FIQRTBk0hryaPA2UtR3eN9kYWZy/mxpgbSYlM4YlBT5Bfk8/us7sx2U0szlnMzbE389Kwl0gKSWLO4Tm4pIsdxTs4UnWEGakzmDVkFhaHxRuV3RQP8vyQ55maMpV9JfvanLxelL0Im8vG02lPd6apL4qmqNzZh2YT5R/Fi8Ne5J6+9/BZ3mdUmiupNFfyWd5n3NP3HuKC47zn2Rwu9p6sIsZvIIfK3XMwB+uXIx2BjPT/fwC8c+gdALac2UJ+TT5Ppz3doXLRJPTP975ZmruUSnMlLw1/qUMzZHxwPDaXjS1ntrCtaBuPD3qc54Y8h1al5f3M91uVN9lNvHfoPSI0A3BaotlXeuHeN0arkfyafIZHD8df68/E5IlsOL2hU3M2naHAWMDG0xtJDEkkq3ovvWPyWpl4bE4b87LmMTRqKCN6jrgs971QulboB/X26aNfb6tn9fHVjO4zmghDBAN7uV2hJs/bxxeHz+J0SSb1n0SYPow5mXN83sadR9/g7YRNC1XfHHMr/cOT+MXKTHbknxsmNw233j74NkerjzIjdQYalYZn0p6hwd7AkqNLvGVP151mfcH6Tk3yHio/xJ6SPTw5+Mk23bOa3OuaC9JdZ3eRU5XD9NTpXsH8XdCWpt+kffnyFhoSNYSCuoIWroTrTq2jxFTCtJRzZg6tSsvL17zc4QRn02Tu7/f8Hqd0MiNtRqsyYxPHEqYP4930BaSfrvFq8p8d+4waaw0z02Z6y0UZoliYvZBlucswWo3MunoWKqFi5tUzOWE8wfrT65lzeA4xgTGMTxpP35C+jEkcw/Lc5eRV57EqfxX3Jd9Hr8Be/HTATwn3C2+l7VeZq/gk7xPuTrybhJCEdp/vUojwi+Bw+WEOlh9kWso0dGodM1JnYHPZWJy92OeL5+CZGhptTkb0GkqVpYpV+as4WH6ARM14duao+Gn/R/jixBfkVucyJ3MOCcEJ3J3QsXIRExSDWqhbeJsZrUYWZC3g1thbOxWY1uS2+ed9f3b3s4GPEmmIZFL/SXx18iuv6aeJ/xyYR5WlijP5t9NDczUZ5QfbnMxuj4zyDCTS+xt/+KqHcUkXnx779IKu44sPsz9Ep9Yx/875pESk4Ahdzf7CIo6VnYteXnN8DWWNZa2C3q4kXSf0HbZ2TTur8ldhspu8K8LfkBTBW5OuxuJw8sKyDG5/aytzthZxU9QD7CreRVZFy4U4pJSU1Vk4VWmiVzPPnTXH11BrreWp1CdZOmMkiZEBTF90wCv444LiiPCLYFvRNmICY7g36V7Anffijvg7WJKzhPyafN7c/SYT1kzglW2vcO/qe1mdv7pdl8I5h+cQpg/joQEPtXlcq9aSEplCRnkGVqeVvSV7+ffBf9MroBeTeTKZAAASgklEQVTj+n63q0p5bfrSXf+cqhy2Fm1lyqApBOraDvseFjUMODcycUkXC7IW0C+sX4dD+7ZockUrqCvg3r73erMxNuFySXblG1GZbuBAxS4enLeWsf/ZyeK9OSw48iHX9bzOG62oVWt5dOCj7C3ZyweZH3BTzE2kRKYAMLrPaPqG9OWPe/9IVmUW01OneyfyZ6bNxOKwMH399BbmJYPGwNTBrUeVi3IWYXFYvlMtH9yavsVpIcoQ5c3PnxCSwB1xd7EkZykf5yxjVOxdrcxLO/IrUKsEEwa6fbv/+u1fifaPZtbwyZTVWVm+oR86VQDPbXqB3OpcZqR1zoSoVWnpHdi7haY//8h8GuwNvDCsc8kDm+pabi5nyqApBOmCAHgq5Sk0Kg0fZH2AyyU5VFjL777cw8e5i3A1pPGrUXfx5zEPYnfZ2V+6v1P3aiK9LB2tSuudYI4LiuOW2FtYcWzFBXnn1VhqeHDtg7yT8Y5X8SgzlbH2xFomJk+kh38P3rzhTWyyAUP019z5r+1c88cNTHh3G//Y9x6JgYMYEX3xKbcvla5bLtHpW+g3eQBcE32NN3JTCMEDw2O5b2gM67NLeX/HSd7ZnI9LxBGY5M9jn7+EzhmHWi1QIWiwOt3uXWpBWOAAjtdEkBiSyOLsxaRFpjEsahhCCJbOGMmjH+xl+qID3JvWm+LaRoyuODBUYa28jdc/y6ZvZABx4f7cFPkIG89s5P6196MWGm7pOYG+QUPYVLKUN3a/wfuZ8+ih74PRbMfYaAcBIQYtgXoVmTW7mdz/WTJON1JlqkUtBP56NQE6Df4694o7/UNS+eTYYkZ+fD0OaUeg5vqg55m9tQB/nZpAvQZ/vYZAvRp/nca97dkfEaj3ppmQUrL3ZDULd59i94kqlk4fSWqsb9OKlO53/6+2/I0ATShVttP4qQLxN9/K8m/PoFIJQg1aQv11BOo1mO0OqkxRqIWWv+z5F4uzVuLCwgnjCR5O+CX//5s8jpXWExagIz7cnz4R/oQYtKhVArVKoFOrCPXXEeavJdjgzj9+ttZGmC6aWlsFsWIcy789g9nupLDazJlqE3ll9RRWm+kZNgJ1z2/oM2g5JeZG/p7nHmmcLXiSVz87TEKEP0IIqhqHoMaPRkcj/o138+6W4wT7aYgK9uPu2Cm8m/17ehh6kqC/hU1Hy6gy2ahqkMTpRnLGuocE3Sg+2WMi2O8kkYF6EgPvIES3gNe3/Ypov0QabU5ONKSTaLiRzZmw3+8MQX5agvw0BPlpCTFoiQjUEaTXXJJGZ2y0Y250z/+kBU1k5YEyzHYn245VsOfMIPQJ/wPg860DOZGzk5uSI4kN8ycqSM+mo+UMiw8lLeoqArQBmOwmXkt7jbED4vhsZgDvbTnOzpJbsEZ/hU5GsTU9luMn8gg2aFEJgRAgPL+9pu8IgQC0rih2FO5hyhfP4qdVc7BiD9dHjybjuB8rK4+iU6voHWqgV4gf4QE6nC6J0yVxSfDXqQnQB6BT6dGpdQwNGcfWvHJqGm1UNdhI1I1izfHP+fxwAVaHRKWrRGNwseC+N7kmpi82ZywGjYGdxTu5Je4WX03XivSydFIjU1s4E0weNJnp66ezJn8ND13VtkLWHCklf9z7R/Jq8sirycNkN/H6ta+zOGcxUkqvkjogfABPpTzFB1kfMKxXABa7oMJahVlWkZ0znpFHNzMmJZrEyECigvTuT7AfPYL0BOp9i+XOptFoD9FVoe/X9NbIA0t+D6N+0+rY1ye/5hc7fsF/R/23VSBPc0xWB5lFRpbmfMa3NSuRUuKSEolEo1Kh1QjUKkm11W2zizREUmmu5J+3/rOFS1q1ycYzHx3gVKWJ+HB/AkLzKZVbiTLP4FSFhbPGc0Ecuh7/Q2gasFWOQtqbFiKRaIKy0YbvRKjcE8c6jQopz/kVS2cQ5sIpIH0vyabSl6KP/hyXJQaHKRm1LQkh9VgdHadT1ardi8j3jQygqMZMbmk9Yf5aHC5JSu8Qls64zqfwefOrvXxa+AdQndN27NU3Yq+9rs3yTeijvkAdcM7N0mUPx1I0Ga1aQ2JkAEaznbK6zmtQuoitIJzYKm/37vPXqYkP9yc+3J+7U3tyb1pv3jv8X7YWbqVfaD/0Mpbyit7UVsdwqtJEuWeFNa1a4B+5G7RVmEvGe3O2eGqKX+9PcNSn4qhPaVEHP/8q9NGrcZX/lHpTcItjmuCD6CK2AxKEREgNjUWPIu2+F6XXqVWEBWjRaVRoVSo0aoFBpyHYT0OIQUtkoJ5BvYNJ6R1CUlQAx8sb2H+qmv0FNRwuqqWoxozaUIA2YjuW4kdAukclTe1R5/c5KhX0cj7A5txyDhXW0lwuvHpnf54f1Y/nNj3HydqTrL1vbYvgsUNFFby85RXUpuupr+5Peb2FzsgVTUg6uoht7g0pkM4ALGd/inSEolOrcLhcHV5HH/UlTms0DuO1LfartfUExH+En86GQSfQa1VMHvgojw9+3FvmhU0vcLz2OF/f/3WnXqqN9kZuWHYD01Kn8cLQc6MRKSVT1k2hvLGcryZ+1WFgXZNs+tnQn1FtqWbJ0SXcl3wf3xR8w6j4Ufz1J3/1lrU6rby45UXya/IRCIQQDA5P4c7I1/gis4QteeVY7K1/2/46NWH+OjRqgcajKJmsThqsDhqsDk7+5Z50KeU1HT60D7pQ6KvlgS8XwrCWARxSSh7+6mEa7Y18ft/nl7wIBbiHXtuKtrGlcAsA74x654I8YRptDs7WWjhba+ZsrRmnlIT56wj116LXqKmz2KlttFFvcZAQEUBabAihnhwcNSYbmcVGzlSZCAvQERmoJyJAh0uCyebAbHNisjow252YrE6cUhLvEd4xoQZUKoHD6cJkc9Joc2CyOjBZ3eeYPOfWWx0U15gpqDRRUGVCr1Hx2HV9GD+kN8u/PcObX+SwcOq13DogqtWz5ZytY9w7O5k4NIY3xw/GYnditrldGoUAtUrgcEqMZju1jXbqLXb89W6hFeSnxeZwUdFgpbLeit3pYlDvYAb0DEKvcbev2eaksKaReosDl3Rre1aHi9pGG7WN7muGBWjpFeLWCkMMWlQqgUq4BWZ4gO6CNGWzzYlKhff+TVgdToyN7pdQaZ2FsjoLOrWKiEAd4U3/l0Ad/rpzWpbTJWmwOKhosFJWZ/Gm4+4TEUDfHgFEBelxSWiwOqi32Km3ODwf93NVm2xUmWxUm6zYHC4cLonDKWm0O6kz26mz2CkzWjDZWiceiwk1MCQ+lMG9gxncO4S+kQHegCSVCnoE6ttsF7vTRWWDlbI6KzUmGyP7RmDQqTFajThcjg7jCBxOFyarE4lESpC4f5PuvyBx75RAo81JeZ2FsnorRrOd2FADST0CiQkzuM2r9VZKas3UNtpRewSYQNBoc3gFWNP/uOkTEaAn2NDx6Gh57nL+tO9PfDnxy055Te0+u5tnNjzD3NFzW6Uy2F60nec2PccfbvgDE/tN9HmNMlMZE9dOJDEkkUVjFqEWat4++DYLjiwAYOX4ld61bDuDlJLaRjvl9VbK6y1U1Fvd3+us1JptOD39xemS+OvVBPtpCdRreG3MVT9gob9zc6tVnPaX7uepb57ijevfuOzRjd0Rm8PFHf/chr9OzVc/+4nXBARuoXb/e7sorjWz8eVbvC8qhSuHyyU5VWXiSLGR4+UN9O0RwLUJ4d5MsAptU1hfyNhVY/nliF/y2MDHOiz/34z/Mj9rPrsf2d3KkUJKyUNfPkSDvYG1961t02lCSsmsTbNIL01nxfgV3heNlJKPcj6iwd7As0OevTwP1wFCiEsS+p1So4UQY4QQeUKI40KIX7ZxXC+E+MRzfJ8QIqFTdz/Ppm91WpmXNY9wv/DvfPKyu6DTqHj1rgHkltbz+aHiFscW7i7gcJGRN8YNVgR+F6FSCZJ6BDJhSAyv3DmAiUNjFYHfCeKC4ugT3Iddxbs6Vf5A6QEGRQxq03NOCMEzac9QWF/YImK7xlLDFye+4Le7fstdK+9iV/EuXr7m5RYjCyEEjw9+/IoJ/MtBhxO5Qgg18C4wGigC9gsh1kopm+fWnQbUSCmThRAPA38DOpgVEV4f/UZ7IyuOrWBR9iLKzeX8fPjP8dP4tX+6Qqe5N7UXH2w/yVvrjzE2tRcuKckrreet9XncNqAH49Iub9ZOBYUrwY29b2RV/iqsTmubkd4u6SK3OpfdZ3eTVZnV7ojgtvjbSA5N5oOsD4gJjOHTY5+yvmA9dpedEH0II3qOYNbVs5iQ3LkUyN9nOjTvCCGuB96UUt7l2f4VgJTyL83KfOMps0cIoQFKgR6ynYuHJhjkrW+7I/8qzBXU2eoY0XMET6c9zYieI7rMh/XHyq7jlTw2bx8BOrXXhhygU7P+5VtaRCsrKPxQ2FG0g2c3PUtcUBw6VeuRarWlmhqrO/hyYPhA/n7L39u1/687tY7Xt78OQKA2kHFJ45iQNIGBEQMvy9zi5eJSzTudcdmMAZpHShQB57t1eMtIKR1CCCMQAbTIKSCEeBp4GiAi3p++oX0BSIlM4YH+D1y21YYUWnNjciQvj+5PWZ07zXRsmIFh8WGKwFf4wXJdr+uY1H9SmzmAAFI0KYzoNYIbet9ApMG3h1UTd/a5k9yUXOKC4hibOLbTOe5/aHRG038QGCOlnO7ZngJcJ6V8vlmZI54yRZ7tE54yPlcSuOaaa+SBA+0nS1NQUFBQaMmVmMgtBpqHR8Z69rVZxmPeCQEubDVnBQUFBYXvnM4I/f1APyFEohBCBzwMrD2vzFrgCc/3B4HN7dnzFRQUFBS6hg5t+h4b/fPAN4AaWCClzBZC/AE4IKVcC8wHPhJCHAeqcb8YFBQUFBS+Z3Qq946U8mvg6/P2vdHsuwVQIqkUFBQUvud8f/yQFBQUFBS+cxShr6CgoNCNUIS+goKCQjdCEfoKCgoK3Yguy7IphKgH8rrk5t8/Ijkverkbo7TFOZS2OIfSFucYIKUMutiTu27lLMi7lKiyHxNCiANKW7hR2uIcSlucQ2mLcwghLimVgWLeUVBQUOhGKEJfQUFBoRvRlUL//S689/cNpS3OobTFOZS2OIfSFue4pLbosolcBQUFBYUrj2LeUVBQUOhGKEJfQUFBoRvRJUK/o4XWf6wIIeKEEFuEEDlCiGwhxIue/eFCiA1CiHzP37CuruuVQgihFkJkCCG+9GwnCiH2efrGJ5503j96hBChQogVQohcIcRRIcT13bVfCCF+7vl9HBFCLBNC+HWnfiGEWCCEKPcsTtW0r82+INz8x9MumUKIYR1d/4oL/WYLrd8NDAIeEUIMutL16CIcwCtSykHASOA5z7P/EtgkpewHbPJsdxdeBI422/4b8C8pZTJQA0zrklpdef4N/E9KeRVwNe426Xb9QggRA/wMuEZKmYI7nfvDdK9+sRAYc94+X33hbqCf5/M0MLuji3eFpj8COC6lPCmltAHLgR/+EvOdQEpZIqU86Plej/uHHYP7+Rd5ii0C7uuaGl5ZhBCxwD3APM+2AEYBKzxFukVbCCFCgJtxr0uBlNImpaylm/YL3EGjBs8qfP5ACd2oX0gpt+Nel6Q5vvrCBGCxdLMXCBVC9Grv+l0h9NtaaD2mC+rRpQghEoChwD4gWkpZ4jlUCkR3UbWuNG8DrwMuz3YEUCuldHi2u0vfSAQqgA89pq55QogAumG/kFIWA/8AzuAW9kYgne7ZL5rjqy9csDxVJnK7ACFEILASeElKWdf8mGeZyR+9H60Q4l6gXEqZ3tV1+R6gAYYBs6WUQwET55lyulG/CMOtvSYCvYEAWps6ujWX2he6Quh3ZqH1Hy1CCC1ugf+xlHKVZ3dZ05DM87e8q+p3BbkRGC+EKMBt4huF264d6hnWQ/fpG0VAkZRyn2d7Be6XQHfsF3cAp6SUFVJKO7AKd1/pjv2iOb76wgXL064Q+p1ZaP1HicdmPR84KqX8Z7NDzReWfwL4/ErX7UojpfyVlDJWSpmAuw9sllI+BmwBHvQU6y5tUQoUCiEGeHbdDuTQDfsFbrPOSCGEv+f30tQW3a5fnIevvrAWeNzjxTMSMDYzA7WNlPKKf4CxwDHgBPCbrqhDFz33TbiHZZnAIc9nLG5b9iYgH9gIhHd1Xa9wu9wKfOn53hf4FjgOfAbou7p+V6gNhgAHPH1jDRDWXfsF8HsgFzgCfATou1O/AJbhns+w4x4FTvPVFwCB2xvyBJCF2+up3esraRgUFBQUuhHKRK6CgoJCN0IR+goKCgrdCEXoKygoKHQjFKGvoKCg0I1QhL6CgoJCN0IR+goKCgrdCEXoKygoKHQj/g8GIAhiZCYivQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(num_states):\n",
    "    plt.plot(posterior[:,i], label=\"State %d\" % i)\n",
    "    plt.xlim(0,100)\n",
    "plt.legend()\n"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}