12_synchro.processing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#default_exp synchro.processing\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# synchro.processing\n",
    "> Processing functions to align stimuli, detect frame timings and correct errors of the display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "import numpy as np\n",
    "import datetime\n",
    "import glob\n",
    "import os\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def get_thresholds(data):\n",
    "    \"\"\"Function that attempts to get the high and low thresholds. Not working very well\"\"\"\n",
    "    max_val = max(data[len(data)//2:len(data)//2 + 10000000]) #Looking for a max in a portion of the data, from the middle\n",
    "    high_thresh = max_val*3/4 # High threshold set at 3/4th of the max\n",
    "    low_thresh  = max_val*1/4\n",
    "    return low_thresh, high_thresh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading the data... 100%      \n"
     ]
    }
   ],
   "source": [
    "from theonerig.synchro.io import *\n",
    "from theonerig.core import *\n",
    "from theonerig.utils import *\n",
    "import matplotlib.pyplot as plt\n",
    "photodiode_data = load_adc_raw(\"./files/basic_synchro/photodiode_data\", sampling_rate=30000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Supposidly, get_thresholds should provide low and high threshold for the data, but the low_treshold is a sensitive value that should be checked manually in a record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def get_first_high(data, threshold):\n",
    "    if np.any(data>threshold):\n",
    "        return np.argmax(data>threshold)\n",
    "    else:\n",
    "        return -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`get_first_high` finds the idx of the first frame higher than the threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def detect_frames(data, low_threshold, high_threshold, increment, do_reverse=True, precision=.95):\n",
    "    \"\"\"Frame detection (or ON signal detection). Capable of finding frame times produced in a regular\n",
    "    fashion:\n",
    "        - data: raw data\n",
    "        - low_threshold: threshold used to detect begginning of each frame.\n",
    "        - high_threshold: threshold used to assign label to the frames, and used to detect the beggining of the reading frame.\n",
    "        - increment: Number of timepoints separating frames. Eg, for 30KHz recording of 60hz stimulus: 30000/60\n",
    "        - do_reverse: boolean to indicate if the reverse detection should be done after detecting the first frame.\n",
    "        - precision: Value that indicates how precise are the events recorded to accelerate the detection. \n",
    "        DLP is very stable (.95) whereas some camera triggers have more jitter (.6). Too low value (bellow 0.5) can\n",
    "        result in an overdetection of frames.\n",
    "        \n",
    "    \"\"\"\n",
    "    assert (precision>0) and (precision<=1)\n",
    "    frame_timepoints, frame_signals = [], []\n",
    "    increment      = int(increment)\n",
    "    safe_increment = int(increment*precision)\n",
    "\n",
    "    first_high = get_first_high(data, high_threshold)\n",
    "    if first_high == -1:\n",
    "        print(\"No high frame detected. Detection can't work.\")\n",
    "        return\n",
    "\n",
    "    frame_timepoints.append(first_high)\n",
    "    frame_signals.append(1)\n",
    "    \n",
    "    if do_reverse:\n",
    "        new_timepoints   = reverse_detection(data, frame_timepoints, low_threshold, increment, precision)\n",
    "        if len(new_timepoints)>1:\n",
    "            new_extrapolated = extend_timepoints(new_timepoints)\n",
    "        else:\n",
    "            new_extrapolated = []\n",
    "        frame_timepoints = new_extrapolated + new_timepoints + frame_timepoints\n",
    "        frame_signals    = [0]*(len(new_timepoints)+len(new_extrapolated)) + frame_signals\n",
    "\n",
    "    i = first_high + safe_increment\n",
    "    while i < len(data):\n",
    "        data_slice = data[i:i+increment//2+(increment-safe_increment)*2]\n",
    "        if np.any(data_slice>low_threshold):\n",
    "            i = i+np.argmax(data_slice>low_threshold)\n",
    "        else:\n",
    "            break #This frame sequence is over. Pass the next sequence through this function if there are frames left\n",
    "        frame_timepoints.append(i)\n",
    "        frame_signals.append(int(np.any(data_slice > high_threshold)))\n",
    "        i += safe_increment\n",
    "\n",
    "    frame_timepoints = np.array(frame_timepoints)\n",
    "    frame_signals    = np.array(frame_signals)\n",
    "    frame_timepoints = frame_timepoints - 3 # A slight shift of the timepoints \n",
    "                                            # to include the begginning of the peaks.\n",
    "        \n",
    "    error_check(frame_timepoints)\n",
    "\n",
    "    return frame_timepoints, frame_signals\n",
    "\n",
    "def reverse_detection(data, frame_timepoints, low_threshold, increment, precision=.95):\n",
    "    \"\"\"Detect frames in the left direction.\"\"\"\n",
    "    new_timepoints = []\n",
    "    new_signals = []\n",
    "\n",
    "    safe_increment = int(increment * (1+(1-precision)))\n",
    "\n",
    "    i = frame_timepoints[0]-safe_increment\n",
    "    while i>0:\n",
    "        data_slice = data[i:i+increment//2+(safe_increment-increment)*2]\n",
    "        if np.any(data_slice > low_threshold):\n",
    "            i = i+np.argmax(data_slice > low_threshold)\n",
    "        else:\n",
    "            break #No low threshold crossing found -> no more frames to detect\n",
    "        new_timepoints.append(i)\n",
    "        i-= safe_increment #We move backward of almost a frame\n",
    "\n",
    "    return new_timepoints[::-1]\n",
    "\n",
    "def extend_timepoints(frame_timepoints, n=10):\n",
    "    \"\"\"Extrapolates points to the left. Not really needed now except for the signals idx that would change\n",
    "    otherwise (and some starting index were set manually)\"\"\"\n",
    "    frame_timepoints = np.array(frame_timepoints)\n",
    "    typical_distance = int(np.mean(np.diff(frame_timepoints)))\n",
    "    extended_tp = [frame_timepoints[0]-(i+1)*typical_distance for i in range(n) if (frame_timepoints[0]-(i+1)*typical_distance)>0]\n",
    "    return extended_tp[::-1]\n",
    "\n",
    "def error_check(frame_tp):\n",
    "    \"\"\"Search error by looking at the time between each frame. \n",
    "    DLP is regular and odd time reveal misdetections.\"\"\"\n",
    "    deriv_frame_tp = np.diff(frame_tp)\n",
    "    error_len_th = np.mean(deriv_frame_tp)+np.std(deriv_frame_tp)*6\n",
    "\n",
    "    error_frames = np.abs(deriv_frame_tp)>error_len_th\n",
    "    if np.any(error_frames):\n",
    "        print(\"Error in timepoints detected in frames\", np.where(error_frames)[0], \n",
    "              \"at timepoint\", frame_tp[np.where(error_frames)[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "detect_frames do frame detection. Works for camera pulses and photodiode data emitted by a DLP. It does it by:\n",
    "* Finding the first frame higher than a threshold\n",
    "* Detecting the frames before if flag do_reverse is set to True\n",
    "* Detect frames\n",
    "* Assign a binary value of if it's higher than the high threshold\n",
    "* Do a quick check on the frames to spot weird results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1feb1c707f0>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2deZwV5Znvvw8QNOj1RgENAgrJ5ZMJOolLX4OTO0lGZpQYR5zP3NxgzIRJvHEm40ySmUkciJqMuUNiFveFhMSFKIGgMYG4oNCCuIKNyr50szdbNyBrSzfd/dw/TjWc5a2m65yqc6pOPd/P53zOOU9Vve/zVr3vr96t3hJVxTAMw0gHvSrtgGEYhlE+TPQNwzBShIm+YRhGijDRNwzDSBEm+oZhGCmiT6UdOBEDBgzQYcOGVdoNwzCMRLFkyZLdqjow3x570R82bBh1dXWVdsMwDCNRiMhml926dwzDMFKEib5hGEaKMNE3DMNIESb6hmEYKcJE3zAMI0WcUPRF5GERaRKRFY5t3xYRFZEBWbaJItIgImtF5Ios+8Uistzbdq+ISHjJ6AHTpsGwYdCrFwwYAKeeCiLHP717Z75PtG3AgEwYw4ZlwgwSx4AB8E//lPnuif3UU4/bevXKtZcr7O7OSzFxRhl2GtOTpHNVbemJKuwBAwq1JUxUtdsP8CngImBFnn0o8DywGRjg2UYCS4GTgOHAeqC3t20xcCkgwHPAZ08Ut6py8cUXa8k8/rhqv36qEO6nX79M2FHGYR/72Cd9n759j2tLkQB1qoWaesKavqouBPY6Nt0F3ARolm0sMENVW1V1I9AAXCIig4DTVPV1z5lfA9f0/NZUIjffDC0t4Yfb0pIJO8o4DMNIH21tx7UlZIrq0xeRq4Ftqro0b9NgYGvW/0bPNtj7nW/3C/8GEakTkbrm5uZiXMxly5bSwzhR2FHGYRhG+ohIUwKLvoj0A24Gvufa7LBpN3YnqjpFVWtUtWbgwIKniINzzjmlh3GisKOMwzCM9BGRphRT0/8wmf76pSKyCRgCvCUiHyRTgx+ate8QYLtnH+Kwl4dJk6Bfv/DD7dcvE3aUcRiGkT769j2uLWHj6ujP/wDDyBvIzdq2ieMDueeRO5C7geMDuW8Cozg+kHtlT+IOZSBXNTMocu65qiKq/furnnJK7sBJr16Z7xNt698/E8a55xYOtJwojv79Vb/+9cx3T+ynnHLcJpJrL1fY3Z2XYuKMMuw0pidJ56ra0hNV2P37lzyIq6qKz0CuZLb5IyLTgc8AA4BdwPdV9aGs7ZuAGlXd7f2/Gfgq0A58S1Wf8+w1wKPA+z3R/xc9UeRATU2N2oJrhmEYwRCRJapaU2Dvge5WFBN9wzCM4PiJvj2RaxiGkSJM9A3DMFKEib5hGEaKMNE3DMNIESb6hmEYKcJE3zAMI0WY6BuGYaQIE33DMIwUYaJvGIaRIkz0DcMwUoSJvmEYRoow0TcMw0gRJvqGYRgpwkTfMAwjRZjoG4ZhpIhUiX5npzJl4XpeW7+b1xp2c9OTS/nOE0u57Y8raTpwhJ+/tJ41Ow+wcF0z3/39cib8bhmTF6xnf8tR7py7jjU7DzBnxU7mrNjpG8eLa3bx69c35dhUlYde2cjL9bkveT/c2s6dL6xl696WHPvuQ608+upGWtrac+xvb3mXh17ZSGdn7jsQZr2zjRdW5vrU2t7BPfPqWbFtf459z6FWfvHSeva1tOXYV20/wN3z1nHkaEeOfc6KHcxYnPuCZlVlZt1Wljfmhn3wyFEemN9QkJ4te1q4/8V6Drfmpufl+mZmLN5C9jsdVJWn3mrkrS3v5ux75GgHd85dx9qdB3Psuw4c4c6569h7ODc9K7btZ/KC9Rzt6Myxz3pnG79/uzHHdrSjk8kL1rOscV+O/XBrO1Nf20TzwdYc+4bmQzwwv4H32nLPVde1z0/P1Nc28eamvTn7trS1c8+8+oJz1fhuCz+Zs4b97x3Nsb+5aS+TF6ynPS89f1y6nTkrduTY2to7uXveuoJrv/tQKz97fi27D+Wmp37XQR6Y30Bre/6131lw7Ts7lZlvbmXl9tywD7W285M5a9i853COfcueFu6au45Dedf+tfW7+eXCDQV5efriLTyfl5ePHO3gvtp6Nu3ODXvn/iM8ML+Bg0dyz9W6XQd55NWNBefq+ZU7mfXOthxbe0cnUxauL8hv+1uOctfcdTQdOJJjb2g6yL219QXlZPHGvUx35OVHX93IgrVNOfu+19bBXXPXsdGRnntr69nfkpuesEnVS1R+9OxqfrFwQyhhzbrxk3x86AdybM0HW/mfk+YBcO+1F3L1x88G4P4X6/nZC+sAWHHbFZx6Uh8Ahk145tixm27/3LHfJ7Jf94lzmPQ3fwpkhOarj2bOT+2/f5oPDzwVgMvuWMCG5sOBw/7gaSfzxndHA7C++RCj73gJgIfG1zD6o2cBMOmZVfzy5Y2Bw862HznawZ/cOgeAH4w9jy9fOgyAJ5c08u0nlgKw+LujOfO0kwE473tzOOyJbJA4Lxl2BjP/8VIAlmzey99Ofh2Amf9wKZcMPwOALz+8mIXrmktKz64DR/jED2sBuPsLF3DNhYMB+OXCDUx6djUAy/7zck47+X2Bw862X3PB2dw97kIgI5xf/OUiAJ79xp8z8uzTALjqvpdZse1A4LCHnvF+Xr7pMgA27znMp3+6AIApf3cxl5/3QSC3DK3/4ZX07iWB09PRqXz4u88C8A+f+hATr/wokBHlf3hsCQALv/MXnNM/887p87///LGbRpD0/NmH+/Obr40CMpWAq+57BYDpXxvFpR/uD8CN097imeU7Aoedbd//3lE+ftsLQG5enr54CxOfWg7A27f+Faef0jdw2KVgL1EBXl2/O7SwtuTV0AAOZNU4GnZlaqSqyqKNx2t5rXk1hHzyay0uFnothqMdnTz11vGaS1eNVFWPCX42+TVWFzuzajbrsmrV2eE9101LJ591uw467W9s2HPs92Lv/Bw52sETdVuP2fdmtUYOt3V/3vxYnFXD3rr3vazfx69fl+Bns3P/kQJbd2TXfNdmpfmdrcdbDy2tmTSs2n4gUNjZ1K45XmvcvOd4GnYeOJ62LsEPSvb5yW4NrM+69gvrj5ehTq/CGLTiuOfw8bCz80FD06Fjv7OvfX4roTuyWyuvrT8e9qasFsiWvZnf77V1HBP8Usi+nl2tq45OzYn/sNdqz28hdLG++ZDTHgWpEv1K8OM5a3m5vvBmc8+8euf+f/qfLxTYdu4/klML6GLSM6t5ellhpv37R950ht3VCsnnsdc3FdieXradr097q8B+97x1NL77XoF9x/5CG8Dldy0ssL20rtnp44TfLcu5QXbx3d8vd4b985fWF9ja2jv5yiOLC+wdnZojwF3kd1N0MepHtQU2VeWCHxRen/lrmo61tvLtLlG58t6XnXF+7deFYXR0Krc/t6bAvqxx37FaZDa/WbSlwAZwxwtrC2yqyn217ny4z9HFsHP/EVbvKLyhfPL2F51hdLXYslneuJ9LJhWe25a2dl5tKCwn+V2IXVz8/+YW2PYebuMjt8wpsD+9bDv//Ju3C+xfemiRM+zP3lN4fVTVWQbnr23i2l++UWD/r2dW8cel2wvsf/GzBc44u1rU5eCEoi8iD4tIk4isyLL9VETWiMgyEfm9iHwga9tEEWkQkbUickWW/WIRWe5tu1dEJPzkxI+pr21y2u+at67HYeT3NZ/I/pKj5upH47st3DprZYE9u5aSzd0+N6tLf+Qu+H5xunCJMvgLmUsMV2zfz/y1hemfvKCBRx3X4nP3vtKNp7m0d6pTDPP7zrv4xvRCoemOuat2Fdherm923txcaQf/G+R9LzYU2JoPtXLH3MJ8+PzKnVw/tfAG9LeTX3OGvd2nVfTkksYC2+qd7lbId55Y5sxzf32/+/rsyRvDAf+KxzzHeQVYsvldp911Y+vodLdm8seYuvjjUncLYkfAFmQU9KSm/ygwJs82FzhfVT8GrAMmAojISGAccJ53zIMi0ts7ZjJwAzDC++SHaVQAv8xcbaz2KZxxpzPKMTefoBdtKGxtAWzb5xbVMMju2jGi5YSir6oLgb15thdUtauj7Q1giPd7LDBDVVtVdSPQAFwiIoOA01T1dc10AP4auCasRPSUmI9Zp5aUNPpCwc5UvJEEXKEw+vS/Cjzn/R4MbM3a1ujZBnu/8+1OROQGEakTkbrm5p53VVSa/JtKZ6fynmPgxtU3CBRM6QTY19LGvS8Wdqms3L6ft7YUdoe4mqYAj7y60Wl3DSC92rDb2aWSP83zRLi6md5r62D+mly7kKlFbtpT2O3j13XiGnxtbe8omGIImdkV+X3D3d1n8qcRdnHoSM8HFIPiGgzdurfF2c2yr6WNdbt6XjP2a80FaeS95pgEIVAwHbGL/KmokJnueNOTy3oUn0DBtNViCdKlEsZsxqaDRwqmxYK73B/t6GTsA6+WHGcQShJ9EbkZaAemdZkcu2k3dieqOkVVa1S1ZuDAgaW4mBduaEGdGBGmOQTIj7U7D/J3DxUOQH5v1krnbAy/vmjXINTh1nZu++OqAvuzy3c4ByC/8qh7IPiCHxQOnoE7M7+5aS/jHy5Mzy1/WMG81bl9rCLuwUBBjk2zy2bn/iN82RH2fbUNzKwr7Ee+7ldvFMy2EqHgeYouxjsGgrfte48LHYOHzy7f4ewXf2bZDg46Zp24Bk4Pt7YzfOKzBfZJz6x25tm/vv8Vp6i4zgnABbcVDj63tXc6B6t/9fIGHnZUELqmhubjGpBvaDrEn/9kfoHdNfCMCC+u2ZUz48kzH5sGmU17RyeX3bGgwH6otZ2r7y8Uz8fe2OycHOAaZAX4F8c4jKoyweF7Q9PBgrEVEZwD1X6s3XmQpT5jWVHRp9gDRWQ8cBUwWo/fHhuBoVm7DQG2e/YhDnvVIsCtf1hxwv268KvV5D+gVQx+/cJ+rYIw2HPI3SrYtKdwKmlQXK0nyJ0KmI3fFMbvOQawAedsq5ccg8MAP32+cFYMwI2/KZz5BDhvEGsCjDeISM7UymxcrZ+29k7nzeexNzY7w7h/fuGAb1DyHzrqDgFnxcOP3YfanNOR76utd7Zo/MqgS9wB52y4+WubnIPSX3+88BonoaeyqJq+iIwB/gO4WlWzq1CzgXEicpKIDCczYLtYVXcAB0VklDdr58vArBJ9Ty1J6DcMQnWlJlqCior6NKjfC6EyESeOdkTXjN/tU4FxE//cfMKavohMBz4DDBCRRuD7ZGbrnATM9Qbh3lDVf1TVlSIyE1hFptvnRlXtqpZ9ncxMoPeTGQN4DsPABnLDIMqKgN/lseuWTE4o+qp6rcP8UDf7TwImOex1wPmBvAuZMOsCxeb3MApKJcpalHFWQjqiFEm/sSPTSCMO2BO5ReIu2JWdE+qnKbGfqRqSGNqU3PLiN9MljBkwkbZcUt4qMtGPiLCusav8+JUpv6KWVDGMtHURnzIYiEq4XYn8U23XJ04LW6ZK9ONw4oP6UBGPA9xoInWj8pfrGH4Doq7r6SdYYaQnRqckFCoh7n7XMi2kSvSN8uHbzHXagpX8tNR2gxAnIfPzpNpmnYVBJW56JvoREVjIKjE4W/4onSS1KR8W/v3AIYTtc5X9+pjjfimiPFdBSWq+NdFPIIEzW3wqgZES99q4kUyCFLck3AhM9EMkrqITp6a/i7DKSdzTGQYJ0JTYlgMjg4l+RIQ3e6f8Jajc4hmnJntSqcQ0RKOQJIxbmOjHnDRUmipRUCoxb9rEsxC/U2LnKjpSJfrlrDSHFVeUy/nGhSQUcN9nIMrqRbwe8gmD6ktPpT04MekS/RCLaLmWYTgYguiHcQMKmt64T6sM7J9P4Jsd7wCwPu14E1bLMqlv+0qV6JeTREzB9HEyStHqCOMR/QTUpsKgpc29hHQQAt+sAxzgfyXdW1zX3m/MKgl940FeQhOnFo2JfpE4l0eIIp4QQg0aQnuE7811vUg6tEHvcIIpO37JP+xYBz8J+N3X1+0q/T3FvmMAJYdcGSpxczPRj4iglzKpmTYcwkl9Gl7yHnU+2Xs42Csx00CQWrqtvWOUTJC8Ep9sFYwgSzZEHWfcqchaNbbgWo+JUzeOH6kS/TjcbCsy7z4OCe+GsJaEDpLOoM3qMM5gzC9DaojTQ3yV8CVdol9pBxKM73LOMRqYjcv1TUBlLzb4rgFk5zAyUiX6cSYuglUJkjBTI/a19Bi9eyDupyrtmOiHSLYwxKmmUolCmNRVHKMk0jzh94rGCKMMSpx8STMnFH0ReVhEmkRkRZbtDBGZKyL13vfpWdsmikiDiKwVkSuy7BeLyHJv272ShBGPEkhC8iryblrHeQn+4Ff8z20QbEmIANg6TSXTk5r+o8CYPNsEoFZVRwC13n9EZCQwDjjPO+ZBEentHTMZuAEY4X3yw0wUrkzWk4wXyovRA+4f+64JH8Iox0lezC0JPrrwfzF8kKmPCU18QGI5T19VFwJ788xjgane76nANVn2GaraqqobgQbgEhEZBJymqq9rZuTv11nHlI04z2IJMoOlEl0nSRUgPyLtaQkpm6VF+MpNEt7UFiXF9umfpao7ALzvMz37YGBr1n6Nnm2w9zvf7kREbhCROhGpa25uLtJFwwhxrf4YVxjiRpBpiNW2yqZ/yzI+CQp7INeVMu3G7kRVp6hqjarWDBw4MDTnrNjmEpfz0SvCAhGnwhYO8UlPJRbyq0SmDRJlEuoGxYr+Lq/LBu+7ybM3AkOz9hsCbPfsQxz2xOJceyemFzxWD6OEcJLCqEkmlaq7hxllp1jRnw2M936PB2Zl2ceJyEkiMpzMgO1irwvooIiM8mbtfDnrmKokVoUzPpofac3bdT+ptge/jHgTq3LvQ58T7SAi04HPAANEpBH4PnA7MFNErge2AJ8HUNWVIjITWAW0Azeqatf6sF8nMxPo/cBz3qe8VEnJ9astJyF5Ybxk2m8gzjnojTi3JKFwBiIB6QlyzhOQnEDEaUzohKKvqtf6bBrts/8kYJLDXgecH8g7I9LMUm0Fy4lvIlORel+inBnkO2UzQBi+N4h0X7ZQsCdyjarBpQdx0ohQnjsIIQwjPlSixWmiHxMCNX0rkFP8F1wrPWy/5MSoRRzsVYx+6QkaaQiX2f/dvtGd3CAhV+JZhKBpr7YbbapEP2oNiaIgufupSw+jUkTpi/tBtggjNIw8kvBAXbpEP8SqY9nEJITXMlZkPnVM8CuESU2PEW+SkK9SJfrVTpzyWxjjcGEUoEpM2YyyWyoJohIOET6wV4GSEqeHBE30IyK0jBWfvBIZURZCmwRSfvxa1IEWXLMLFBkm+nGn2mqNRa5OWuz+caph+RF/D93E6Ulvo+eY6BdJnGaWnIi0FE5XDTPJC645F6yqwAqecbop+Y/RxMnLeGOiHyJR6IK9kLtEErDqoWGUk1SJfmcZxc80pfyEMb3VKDMJuEDVtnxEqkQ/3t0cyWxuR0t1pTRoapJacQhlinDpQRg+pEr04/3gRM9Liv9TlvGn3Esl2AvaS6cS+Sqp3W/+T5fHp3SmSvTLSZymClam0FYg0viUq1CIcqnoSpBg13tMjLTdFxP9Iqlk4Qv+YvQE5MQQcHXf9UqD0iSEJN+wqolUiX7Uffo90dbgSyhEubRyPJ5MjFIMfLt3AsYZxmXwC6ISYhgXAQ5rnCNQODFJe6VIleiXk4oU5PJHGfselYo8fxbWw9gpF6fIqMADj3Eao0iV6KeklyNxxKlAhEGU+awSeTj4qq6lO1ldOcIfW0/fiIQ43eyizORBBj6r7D5TdSIZp+sTJ1/CwEQ/NlTZEyAxIaxWRCWe8Yj3FONoX6xjREdJoi8i/yoiK0VkhYhMF5GTReQMEZkrIvXe9+lZ+08UkQYRWSsiV5TufuU4UcYOq8AmtfxUonbkrOn77FuRQeyyx2iUmyR0VRYt+iIyGPgGUKOq5wO9gXHABKBWVUcAtd5/RGSkt/08YAzwoIj0Ls39YERdA4n3E79GNZAATemGRDtfNZTavdMHeL+I9AH6AduBscBUb/tU4Brv91hghqq2qupGoAG4pMT4Y0sSCmdcfEzC8gRx6bLwa6HExD3A/1y5rpv/FMyYZM6AJMHrokVfVbcBPwO2ADuA/ar6AnCWqu7w9tkBnOkdMhjYmhVEo2crQERuEJE6Ealrbm4u1kUjKYQ02Op+OCv+xbASr3SMy1kJKu6+e8clQQmglO6d08nU3ocDZwOniMiXujvEYXPWCVR1iqrWqGrNwIEDi3UxloSRN4M/kRvMXk34CmdCRSLtXYiVSH0oD+bFqLCV0r3zl8BGVW1W1aPAU8CfAbtEZBCA993k7d8IDM06fgiZ7qBE4myqZilJpHO1A9rjRLm1tso0PzUkoIEWCpXoxipF9LcAo0Skn2SGrEcDq4HZwHhvn/HALO/3bGCciJwkIsOBEcDiEuIPTJzutj3F6XLAZMSpdmgvhUkvsbluEb6eMwn0KfZAVV0kIk8CbwHtwNvAFOBUYKaIXE/mxvB5b/+VIjITWOXtf6OqdpTof6zIFtewMkpngHmISc2bkb4YPbR5+uUnqS/viIu2x4k4TeUsWvQBVPX7wPfzzK1kav2u/ScBk0qJs1qJUZ4oO2G1wKpt+YOkksQWdWgkoBzbE7llJmhxCONmEPcyOOXlDZGFneabadwIMmUzqSRhqqmJfooJ54YS7I7inMIV8KYUqNsjAasexsiVSKnEzDUnkU6yCBZ4JcbbUiX6Ma/wGiUS5PomQWejHeuILOhQCOpfEmrYcSFdoh+i6rvCiqIbJcqXd6QBE4N4E1aZidMMtbiTKtFPInHvj48TrnOVhNcl+j48l1AhC8PrOHW/RUnS5ukb3eD7mj6f/YOIu70jt+ckeWllIx4EEeYktCxN9IvkRFoS5aVPQi0oCZk/CNV234w0PUEqMNWVTRKBiX5MsMxfSJSzPZJwvgPVMH0SFIa42zt/q4tUiX4Sm+iughK0uyZOqS57K6Xqlu71E/c4XeWe47+0cpSRRhh0ArJVukQ/meUiMiqRP0OpvfuWrMILnIAyGOg9vpWodQctN36Vq2AtFz+731hZEq50PEiV6CeRcJZ1LT2MpFKJ8Y8ktiiN9GCib1Q1sar/pWBJjbQTq/zmg4l+RCShb68SRDppxLUgaUhv5TKMKKhEPkyV6JezkpR6TQnhBATtmon2hmJV7HxCmRmUgJJSbRWEVIl+tRC8n9pdOuMuY4FnKTn2933/bFEe9QxfIYv7CQ+JIC9G9yPSm0FKroMfJvpFUq61d8KgEq9u9KPsr0usslpaEgaJw1mGIYRADCepEv24ijIEE8OwykNFpmyG0e0TKL50qEeU6Qxtmmg4wcSaJGS3VIl+mJTr4rq6LGJ870oOCSicUVKJvvRQWgB+9phczySMUZQk+iLyARF5UkTWiMhqEblURM4QkbkiUu99n561/0QRaRCRtSJyRenuB8Xk0kg2cRKVJHQ1OYl0uCD+56TUmv49wBxV/RPg48BqYAJQq6ojgFrvPyIyEhgHnAeMAR4Ukd4lxu/PtGkwYECmCtD12bUrvPC/8IXcsEWg5uLj27/xDfdxAwa47aNGOc26dWuhcdUq577y6U+7w/7Yx9z2W291213827+57X37uu1XXum2B7gG8pGPuDcMH+4067PPFYbxztvuMC67zG3v435ttN4S4Fz9+Ha3/eyz3faLLnLbX3ih0Lat0b3vmDFue79+bvu3v+22u/jOd9z2U09122+4wW135eV5c937Dh7stp91ltv+I59z7uKrX3Xb/ZoLX7yu52FP+qHbfs65bvv55/v7MmBARsdCpmjRF5HTgE8BDwGoapuq7gPGAlO93aYC13i/xwIzVLVVVTcCDcAlxcbfLdOmwVe+Anv2RBJ80tDOzrLHKVEOoHR09HjXXkH98AlbQ+k/qEAtMMC5CotAbzDzuT6++acC6QkD8Vuaorv8uWdP5gYVsvCXUtP/ENAMPCIib4vIr0TkFOAsVd0B4H2f6e0/GMi+1Td6tvC5+WY4ejSSoLsjWxgq0Qj3yz7qt0hXXDpCK4BEOqMpvee1O/yEz02cuklK96XoPNHWltGzEClF9PsAFwGTVfVC4DBeV44PzndiO3cUuUFE6kSkrrm5ObhnW7YEP8aoSlZ88H847cEEKN3YTSy69Peo4hWynpUi+o1Ao6ou8v4/SeYmsEtEBgF4301Z+w/NOn4IsN0VsKpOUdUaVa0ZOHBgcM/OOcdpTmLmTXNtPCgaxto2pQdRESrht28LMoHlLNb46FmxFC36qroT2CoiXaNto4FVwGxgvGcbD8zyfs8GxonISSIyHBgBLC42/m6ZNAne975Cn8sooEH7tFNTTEJ5dN/vCeNknsWi+nsjijNKouxSi0vLLXQ/+vbN6FmIuKcq9Jx/AaaJSF9gA/AVMjeSmSJyPbAF+DyAqq4UkZlkbgztwI2qGs2ozHXeaPs3v2mDuYD2jm6SVJT4il6vCNMTo3PlEpDAN4IQ0hNUyCKtXMXo+rgINeX9+8M99xzXs5Aoacqmqr7jdcN8TFWvUdV3VXWPqo5W1RHe996s/Sep6odV9SOqWji/Lkyuuw527848htv1KaaryI/f/jY3bFV4883j2x+433mY7NntDm/RIrd98JDCMEaOdO/70ktOs77tM23xBz9w213ceafb3tbmtj/nc3k/+MGex1lf7zTrpk3u/f2mLbpYsMBt95kAoLcFOFf/cZPbvmOH2+53fcZ8tsAkQ4c6dgSZM8cdRkuL237HHW67i5/+1G0/eNBtnzzZbT/XMW3x8sudu8qOne4w/Kb8Tpzotrt4+GG33e+G+psAs2duucVt37LZbV+xwt+X3btDF3ywJ3LLjl++8ltcbOeBIxF6Uwl6Xmts74yyeyMYOw+81+N9D7znM3MsjBfilB5EIojTUFZ7R/nP+v/5xes8/obPjaJETPQjImieDSJw9U2HAoX98CsbnfY7567rcRh+a7t0RCjMP5mzJtD+R9t7/jxC0LVqpi92PFgUcN+G5mDXzcWO/e5KQJtP2jfuPlxynPfVNjjtflf+7nnuFtq2fYfAtAgAAAvLSURBVIU3zpfrfVq+AXlwwfoe79vpk2ePdrjP4b8/sbTHYd/3ovtcbfa5DpfftdBpX7xxL7f8wacVUCIm+mXm3RZ3d8gPn10dWZxPLPF5ijMAfi2Rp97a5rT//SNvOu3vtvT8+YlXGtyC4HejeX1D6eM3TQdbSw7Dj9rVTU77jAA3FD++9us6p93vhn+otb3HYR8MsC/A7kOln8Ote326pULAr6yNuDm6Hucv/sqn+7YCpEr0y/oiDJ+a5GV3uPvd396yr+QoWwPUdMOipS2YIARpGbS0ucf575rX8xaKH6u273faP/HD2pLD9uMhHwF+zKcZv6yx9Dzhh19tPAidEZanv3nwNaf9ntrS/T5wJFierTYk7m8Eqqmp0bo6dy0mKBf84AX2Bahpdsf9X7yQqz6Wu5bKssZ9XH3/q0BG82N+agNx8vt6ceRo+W8qhpFmNt3+uaKPFZElqlqTb6/qmv7O/UdC6dMshmoSfMAE3zCqhKoW/VE/quUvfrYAyHQrhFXLh0xf9pGjme4HVeUXL63n1lkrQwvfMAwjCkp9OCsRfGvG24EGrnrCi2ua+K9nVnHL50ayfNt+fvRcsJkmhmEYlSAVov+Hd5xL/JTM429s4fE3bHE3wzCSQ9V278xeGo3QG4ZhJJmqFf1vTPd5tN0wDCPFVK3oG4ZhGIWY6BuGYaQIE33DMIwUYaJvGIaRIkz0DcMwUoSJvmEYRoow0TcMw0gRJvqGYRgpomTRF5HeIvK2iDzt/T9DROaKSL33fXrWvhNFpEFE1orIFaXGbRiGYQQjjJr+N4HsV9FMAGpVdQRQ6/1HREYC44DzgDHAgyIS71fbG4ZhVBklib6IDAE+B/wqyzwWmOr9ngpck2WfoaqtqroRaAAuKSV+wzAMIxil1vTvBm4Cst+wcZaq7gDwvs/07IOB7JeBNnq2AkTkBhGpE5G65ubmEl00DMMwuiha9EXkKqBJVZf09BCHzfl+KVWdoqo1qlozcODAYl00DMMw8ihlPf1PAleLyJXAycBpIvI4sEtEBqnqDhEZBDR5+zcCQ7OOHwLY+seGYRhlpOiavqpOVNUhqjqMzADti6r6JWA2MN7bbTwwy/s9GxgnIieJyHBgBLC4aM8NwzCMwETx5qzbgZkicj2wBfg8gKquFJGZwCqgHbhRVTsiiN8wDMPwIRTRV9UFwALv9x5gtM9+k4BJYcRpGIZhBMeeyDUMw0gRJvqGYRgpwkTfMAwjRZjoG4ZhpAgTfcMwjBRhom8YhpEiTPQNwzBShIm+YRhGijDRNwzDSBEm+oZhGCnCRN8wDCNFmOgbhmGkCBN9wzCMFGGibxiGkSJM9A3DMFKEib5hGEaKMNE3DMNIESb6hmEYKcJE3zAMI0UULfoiMlRE5ovIahFZKSLf9OxniMhcEan3vk/POmaiiDSIyFoRuSKMBBiGYRg9p5Safjvw76r6UWAUcKOIjAQmALWqOgKo9f7jbRsHnAeMAR4Ukd6lOG8YhmEEo2jRV9UdqvqW9/sgsBoYDIwFpnq7TQWu8X6PBWaoaquqbgQagEuKjd8wDMMITih9+iIyDLgQWAScpao7IHNjAM70dhsMbM06rNGzucK7QUTqRKSuubk5DBcNwzAMQhB9ETkV+B3wLVU90N2uDpu6dlTVKapao6o1AwcOLNVFwzAMw6Mk0ReR95ER/Gmq+pRn3iUig7ztg4Amz94IDM06fAiwvZT4DcMwjGCUMntHgIeA1ap6Z9am2cB47/d4YFaWfZyInCQiw4ERwOJi4zcMwzCC06eEYz8J/B2wXETe8WzfBW4HZorI9cAW4PMAqrpSRGYCq8jM/LlRVTtKiN8wDMMISNGir6qv4O6nBxjtc8wkYFKxcRqGYRilYU/kGoZhpAgTfcMwjBRhom8YhpEiTPQNwzBShIm+YRhGijDRNwzDSBEm+oZhGCnCRN8wDCNFmOgbhmGkCBN9wzCMFGGibxiGkSJM9A3DMFKEib5hGEaKMNE3DMNIESb6hmEYKcJE3zAMI0WY6BuGYaQIE33DMIwUYaJvGIaRIsou+iIyRkTWikiDiEwod/yGYRhppqyiLyK9gQeAzwIjgWtFZGQ5fTAMw0gz5a7pXwI0qOoGVW0DZgBjy+yDYRhGImhr7ww9zHKL/mBga9b/Rs+Wg4jcICJ1IlLX3NxcNucMwzDihKKhh1lu0ReHrSBVqjpFVWtUtWbgwIFlcMswDCN+9BaXZJZGuUW/ERia9X8IsL3MPhiGYSSC3r2SL/pvAiNEZLiI9AXGAbOjiOjPRwyIIljDMIyyIRHU9PuEHmI3qGq7iPwz8DzQG3hYVVdGEddj138iimANwzASTVlFH0BVnwWeLXe8hmEYhj2RaxiGkSpM9A3DMFKEib5hGEaKMNE3DMNIESb6hmEYKcJE3zAMI0WY6BuGYaQIUQ1/QZ8wEZFmYHORhw8AdofoTrkwv8tLUv2G5PpufkfPuapasHhZ7EW/FESkTlVrKu1HUMzv8pJUvyG5vpvflcO6dwzDMFKEib5hGEaKqHbRn1JpB4rE/C4vSfUbkuu7+V0hqrpP3zAMw8il2mv6hmEYRhYm+oZhGCmiKkVfRMaIyFoRaRCRCRXyYaiIzBeR1SKyUkS+6dnPEJG5IlLvfZ+edcxEz+e1InJFlv1iEVnubbtXvNfpiMhJIvJbz75IRIaF6H9vEXlbRJ5Oit8i8gEReVJE1njn/dKE+P2vXh5ZISLTReTkuPotIg+LSJOIrMiylcVXERnvxVEvIuND8PunXl5ZJiK/F5EPxM3vSFDVqvqQeSPXeuBDQF9gKTCyAn4MAi7yfv83YB0wEvgJMMGzTwB+7P0e6fl6EjDcS0Nvb9ti4FIyL5Z/DvisZ/8n4Ofe73HAb0P0/9+A3wBPe/9j7zcwFfi/3u++wAfi7jcwGNgIvN/7PxP4+7j6DXwKuAhYkWWL3FfgDGCD93269/v0Ev2+HOjj/f5xHP2O4lOxiCNLUOaCPJ/1fyIwMQZ+zQL+ClgLDPJsg4C1Lj/JvFLyUm+fNVn2a4FfZO/j/e5D5klBCcHXIUAtcBnHRT/WfgOnkRFPybPH3e/BwFZPFPoAT3tiFFu/gWHkimfkvmbv4237BXBtKX7nbfsbYFoc/Q77U43dO12FqItGz1YxvKbehcAi4CxV3QHgfZ/p7ebn92Dvd7495xhVbQf2A/1DcPlu4CagM8sWd78/BDQDj3jdUr8SkVPi7reqbgN+BmwBdgD7VfWFuPudRzl8jbpcf5VMzT1pfgemGkXf9fr4is1LFZFTgd8B31LVA93t6rBpN/bujikaEbkKaFLVJT09xMeHsvpNpnZ1ETBZVS8EDpPpavAjFn57/d9jyXQjnA2cIiJf6u4QHx/Kfb57Qpi+RpYGEbkZaAemleBD2f0ulmoU/UZgaNb/IcD2SjgiIu8jI/jTVPUpz7xLRAZ52wcBTZ7dz+9G73e+PecYEekD/Hdgb4lufxK4WkQ2ATOAy0Tk8QT43Qg0quoi7/+TZG4Ccff7L4GNqtqsqkeBp4A/S4Df2ZTD10jKtTewehVwnXr9L0nwuxSqUfTfBEaIyHAR6UtmUGV2uZ3wRvUfAlar6p1Zm2YDXSP448n09XfZx3mzAIYDI4DFXnP5oIiM8sL8ct4xXWH9b+DFrIxbFKo6UVWHqOowMufuRVX9UgL83glsFZGPeKbRwKq4+02mW2eUiPTz4hsNrE6A39mUw9fngctF5HSvdXS5ZysaERkD/Adwtaq25KUntn6XTCUHFKL6AFeSmS2zHri5Qj78LzLNuGXAO97nSjL9fLVAvfd9RtYxN3s+r8WbFeDZa4AV3rb7Of4k9cnAE0ADmVkFHwo5DZ/h+EBu7P0GLgDqvHP+BzKzJZLg923AGi/Ox8jMGoml38B0MmMPR8nUYq8vl69k+t0bvM9XQvC7gUx/e1f5/Hnc/I7iY8swGIZhpIhq7N4xDMMwfDDRNwzDSBEm+oZhGCnCRN8wDCNFmOgbhmGkCBN9wzCMFGGibxiGkSL+P2o4YeXqUvkpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "frame_timepoints, frame_signals = detect_frames(photodiode_data, 200, 1000, increment=500)\n",
    "plt.figure()\n",
    "plt.plot(photodiode_data)\n",
    "plt.scatter(frame_timepoints, frame_signals*800+600, c=\"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def cluster_frame_signals(data, frame_timepoints, n_cluster=5):\n",
    "    \"\"\"Cluster the `frame_timepoints` in `n_cluster` categories depending on the area under the curve.\n",
    "        - data: raw data used to compute the AUC\n",
    "        - frame_timepoints: timepoints delimitating each frame\n",
    "        - n_cluster: Number of cluster for the frame signals\"\"\"\n",
    "    frame_aucs = np.fromiter(map(np.trapz, np.split(data, frame_timepoints)), float)\n",
    "    if frame_timepoints[0] != 0: #We need to remove the first part if it wasn't a full frame\n",
    "        frame_aucs = frame_aucs[1:]\n",
    "    frame_auc_sorted = np.sort(frame_aucs)\n",
    "    deriv = np.array(frame_auc_sorted[1:]-frame_auc_sorted[:-1])\n",
    "    deriv[:5]  = 0 #removing tails values that can show weird stuff\n",
    "    deriv[-5:] = 0\n",
    "    threshold_peak = np.std(deriv)*3\n",
    "    n          = n_cluster - 1\n",
    "    idx_gaps = np.zeros(n+3, dtype=\"int\")\n",
    "    tmp_deriv = deriv.copy()\n",
    "    zero_set_range = 10#int(len(deriv)*0.05) #Around the peaks, we set the values to 0 around\n",
    "    for i in range(n+3): #Detecting more peaks than needed and then taking them starting on the right\n",
    "        if tmp_deriv[np.argmax(tmp_deriv)] < threshold_peak:\n",
    "            if i<n_cluster-1:\n",
    "                print(\"Less transition in AUC detected than needed, results will be weird\")\n",
    "            break\n",
    "        idx_gaps[i] = np.argmax(tmp_deriv)\n",
    "        tmp_deriv[idx_gaps[i]-zero_set_range:idx_gaps[i]+zero_set_range] = 0\n",
    "    idx_gaps = np.sort(idx_gaps)\n",
    "    idx_gaps = idx_gaps[-(n_cluster-1):]\n",
    "    thresholds = np.zeros(n, dtype=\"float\")\n",
    "    for i, idx in enumerate(idx_gaps):\n",
    "        thresholds[i] = (frame_auc_sorted[idx+1] + frame_auc_sorted[idx])/2\n",
    "\n",
    "    return np.array([np.sum(auc>thresholds) for auc in frame_aucs], dtype=int)\n",
    "\n",
    "def cluster_by_epochs(data, frame_timepoints, frame_signals, epochs):\n",
    "    \"\"\"Does the same thing as `cluster_frame_signals`, but working on epochs around which the\n",
    "    number of cluster can differ. Useful when a record contains stimuli with different signals sizes.\"\"\"\n",
    "\n",
    "    frame_aucs = np.fromiter(map(np.trapz, np.split(data, frame_timepoints)), float)\n",
    "    if frame_timepoints[0] != 0: #We need to remove the first part if it wasn't a full frame\n",
    "        frame_aucs = frame_aucs[1:]\n",
    "\n",
    "    max_cluster = max([nclust-1 for (_,_,nclust) in epochs])\n",
    "    \n",
    "    for start,stop,n_cluster in epochs:\n",
    "        n          = n_cluster - 1\n",
    "        norm_clust = max_cluster/n\n",
    "        frame_auc_sorted = np.sort(frame_aucs[start:stop])\n",
    "        deriv = np.array(frame_auc_sorted[1:]-frame_auc_sorted[:-1])\n",
    "        deriv[:5]  = 0 #removing tails values that can show weird stuff\n",
    "        deriv[-5:] = 0\n",
    "        idx_gaps   = np.zeros(n, dtype=\"int\")\n",
    "        tmp_deriv = deriv.copy()\n",
    "        zero_set_range = 10\n",
    "        for i in range(n):\n",
    "            idx_gaps[i] = np.argmax(tmp_deriv)\n",
    "            tmp_deriv[idx_gaps[i]-zero_set_range:idx_gaps[i]+zero_set_range] = 0\n",
    "        idx_gaps = np.sort(idx_gaps)\n",
    "        thresholds = np.zeros(n, dtype=\"float\")\n",
    "        for i, idx in enumerate(idx_gaps):\n",
    "            thresholds[i] = (frame_auc_sorted[idx+1] + frame_auc_sorted[idx])/2\n",
    "\n",
    "        frame_signals[start:stop] = np.array([np.sum(auc>thresholds)*norm_clust for auc in frame_aucs[start:stop]], dtype=int)\n",
    "    return frame_signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def cluster_by_list(data, frame_timepoints, frame_signals, stim_list):\n",
    "    \"\"\"Assign the stimulus identity values from stim_list to the frames in data. stim_list contains only the \n",
    "        sequence of stimuli. Those need to be expanded. Opposed to cluster_frame_signals and cluster_by_epochs no \n",
    "        AUC operation is performed. \n",
    "        Input:\n",
    "            - data: raw data used to compute the stimulus times\n",
    "            - frame_timepoints: timepoints delimiting each frame\n",
    "            - frame_signals: binary 1-D numpy array indicating if high_threshold was passed in 'detect_frames'\n",
    "            - stim_list: 1-D numpy array containing the sequence of the stimuli presented\n",
    "        Output:\n",
    "            - frame_signals: [same size as frame_timepoints] stim_signals list containing the correct value from \n",
    "                stim_list at every entry\"\"\"\n",
    "    # Determine the stimulus on- and offsets and their location relative to data\n",
    "    stim_change = np.where(frame_signals[:-1] != frame_signals[1:])[0] \n",
    "    stim_change = stim_change + 1 # since I need to compare to [1:] all values are shifted by 1\n",
    "    #stim_idx = frame_timepoints[stim_change]\n",
    "    \n",
    "    # QDSpy currently is set to emit a short peak at the end to indicate the end of the stimulus presentation\n",
    "    # This peak needs to be ignored\n",
    "    epoch_end = stim_change[-2:] # return it for future analysis\n",
    "    stim_change = stim_change[:-2] #add it to the no stimulus category\n",
    "\n",
    "    # Split into on times & values vs off times & values\n",
    "    stim_ons = stim_change[0::2]\n",
    "    #stim_ons_idx = stim_idx[0::2]\n",
    "    stim_offs = stim_change[1::2]\n",
    "    #stim_offs_idx = stim_idx[1::2]\n",
    "    \n",
    "    # Replace the frame_signal entries with the stimulus codes\n",
    "    frame_signals[frame_signals == 0] = -1 # To avoid confusion with the '0' stimulus code\n",
    "    for i,stim_type in enumerate(stim_list):\n",
    "        frame_signals[stim_ons[i]:stim_offs[i]] = stim_type\n",
    "\n",
    "    return frame_signals, stim_ons, stim_offs, epoch_end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Frame signals are then refined using cluster_frame_signals of the signals to attribute them a value in a defined range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame_signals = cluster_frame_signals(photodiode_data, frame_timepoints, n_cluster=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1feb222a198>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO19ebwdRZn2896sJJCwBQwJEMCIgKMikXH53AYZEB3BGXVAHHHGmXz6oY7LqEF0ZAZxHJ1xEPewOKABRASJRFCIIKuEGyB7QkLWm/VmX+/Nvee83x/dfbqquqq6ejnn9On08/slt053LW9VV71bbcTMqFChQoUKhza62k1AhQoVKlRoPyphUKFChQoVKmFQoUKFChUqYVChQoUKFVAJgwoVKlSoAGBouwmIw7HHHsuTJk1qNxkVKlSo0FGYO3fuVmYe5xq/8MJg0qRJ6O7ubjcZFSpUqNBRIKI1SeJXbqIKFSpUqFAJgwoVKlSoUAmDChUqVKiAShhUqFChQgVUwqBChQoVKqASBhUqVKhQAZUwqFChQoUKOMSFwUCtjjvnrMW67fsBAE8s39oIz+/ZiVnzNwIA1u88gAcWbGwbnc3EE8u34pqZi8DM2LSrDzPnbQAz48DBGm57enWjPe6e24O12/a3l9gmYLBWxz3P9TTq+eDCTVi6aTcAYO6aHfhl9zoAwKqt+/DgwnL2gUUbduHOOWsBAOu278edc9aCmbGvfxDXzVqM9TsPYLBWx6/m9mDz7r42U5s/6nXGTY+vxKqt+8DM+M28DVizbR8A4KkVW/HYi70AgBfW7cTvF21qJ6lNxSEtDLpX78C0exbgK79eiL6BGj588zP4wI+fBgB88CdP48rbn0PfQA1X3bMAn5jxHLbu7W8zxfnjwzc/g/99ajV6dhzAl341H5++43ms2LIXf1i6Bf963yJ87w/LsW77fvzLL+fhazMXtpvc3PHc2p343F3z8KVfzUf/YA0f//lcfOr25wEAf/Ojp/CFu+fj4GAdH/3pHHz8589hT99AmynOH5dN/xOm3bMAm3b14crbn8O0exZg5dZ9+O2Cjbjx8VW4+fFVWLRhNz7/y3m49v7F7SY3dyzZtBtfn7UEX7x7Hjbs6sOn7ngen/3FCwCAD930DD5yyxwAwCU/eBJTfza3naQ2FYe0MBis1wEAq7ftQ63uXfKzydd8+ga8d7U6NzSD/sF6G6h0x8HBOnbsO5goTRd5f2t1bmjEB2t1DNS8uu7YP4D+wRoAYM6q7fkR2yQM1upIcmFTUM91O/bD7w5YvmWvFIfBWONbRYO1Yl8GxcwYrCXrp7v7BgEANWas3uppxLU6o8/v732DtUY7rezdlyO1zcGabfuwr3/QOf6A/0037OzDgF/nF9btbAptRcYhLQwCnkEAXIZ40W+F+/Qdz+Psax/KNU+C0E5EueadN5gZL7/6AVwzc1GCNELYoRcUuwcAX5+1BC+/+oGGcpMUulSu46MoeNu3H21o8y4IxjVRZ9Uzb8QKAyK6hYi2ENFC4dm3iWgpEc0nonuJ6Ejh3VVEtIKIlhHRBcLzc4hogf/uBioAZwk+PBEZGb34tOCyAA9m8Geq9dQxxqILw4C8W592P5IlqCeBjN+34NWWcNvTqwGEVm8SiN+3k+qsw9w1O5zjdnhVc4OLZfC/AC5Unj0E4FXM/GoALwK4CgCI6EwAlwI4y0/zQyIa4qf5EYCpACb7/9Q8Ww6p87eRjnbCRdgdKm1TBuswACEHXatEwsGG0PLtnO/bDMQKA2Z+DMB25dnvmTlwyv0JwEQ/fDGAO5m5n5lXAVgB4FwiGg9gDDM/zV5r3wbgkrwqkRYNywDl7uw2uGjDndI2aciU69n5bqIsMLnMDiUmKXoLDjXkMWfwDwAe8MMTAKwT3vX4zyb4YfV5eyFJA4foh8Z48NxEGsZQ9AGShWG5+otb2Qe6V2/Hsk17EqXJTF+MddjqLrC7Jau3kllAZRWMmYQBEV0NYBDAjOCRJhpbnpvynUpE3UTU3dvbm4VEZ5gmD4v24Wt1xpY9+a/1NtXT0wpzLy4W2/cdxDd+u0S7Mmbqbd24w18XL8JGJjNr65h0TshlkjkvvP/HT+OC6x+LPO9evR0XXv8Y+gZqxrRZ6fTmjfTPW4XfzNuAV1/zeyzo2RV598nbn8OtT63OpZykdSoYS8gNqYUBEV0B4D0ALudwlPUAOFGINhHABv/5RM1zLZh5OjNPYeYp48Y5X9STGEkHTCsZwYXXP4Z/vPXZyPPrZi3BudfNxs797ktIDw7W8U+3dWP5Zjctk5WatqPzXzNzEaY/thIPL9kSeff7xZtx1T0LIs9tdP7Tbd045arfatIEE8joGB/QNb9ZhKWb9mD55r3xkVNCbct29IHHl3uK4OKNUWFw//yN+FqCVWOAp2AEm8lEhA4CNxdBh3STxEglDIjoQgBfAvBeZha3pc4EcCkRjSCiU+BNFM9h5o0A9hDRG/xVRB8BcF9G2p3xT7d147Qve4zg4GAdL/V6g0haWlowN9HSTXu0jPDhJZsBALsOuJvPL6zbiYcWb8aX740yUBFF6uTBuva0SyRV6NoyFYrUSBaofXXXgQFMmjYLDy/enDgPAgmuwrwobD3e+q1H8LZvPxp53onzY82Ay9LSOwA8DeB0Iuohoo8B+D6AIwA8REQvENGPAYCZFwG4C8BiAA8CuJKZA1v2EwBugjep/BLCeYam46HFmxtM5ep7F+C8//4jtu87KK2fN/UBNoSLDmZu7JZV3SMrtuyNWBaRpaXSSqvW1zxgOknKVuMu3rAbK3vt2rP8fU2uQn38diOubfoGapjxzBowc8Mq/OGjK9Ll3YaKp2HMNrfuXmEj2m1Pr8akabOwt39Qu8/AJvOK5jrOC7F3IDPzZZrHN1viXwfgOs3zbgCvSkRdE/DUS9sAAPv6BxUGWJ4PzMz4xbPrMO2eBZj9+bdhsB6upQeAd37njzhh7Eg8ddV5UrpgADBk07kdTZNmaeTuA/Ku04tueBwAsPqb75ae79h3EPc+vx5//+ZJYXnkVs+id5NBwZL6798vw42Pr8Ixo4dj3BEjpHhPv7QN48eOxKRjRzeeqfMEnTgmDipzTDPnbcARI4biHa88Tnp+yxOrAAC9e/r1cyOWMjqvVdwQKww6FfN7duKwYUMw+fgjjHFEf3GZ1pgzA7OXem6RFVv2Njr+nNXhCuENu6KT0EWsXZIm/9nTq53ifeHu+Xh4yWacfdKRUqWLWH8bbAKTGdjmH02yt7+GsYd5tQvqeNmNfwIQFZRienEFUTvaJo1L6ldz10u/P32Hd86UWk+jJ6AECkFalPY4is/fNQ//8/CL1jhSZ3dZb58PaS1BQ8tnxsqt8efJRNbba/yoLXUXN9xE7nCNG7jIBmrhVLk3b2RwE4lLDwvUC1xpYWZ8d7Y3Fp5fm+7MnXYwwDRlDiQ8l0mdLxT7w6GG0goDFyQ9l6ZTwJC1qqQWjcv8SRFRd6xnt39Ugap5ulmHCYkqCDbsjF+OrI6BooyPJC5D1z5gEgBFW17cSpRaGOg+7PqdBwAEjEDwpZfMPAwGUJ29f0lRFG04rYbmugop6UqSInUBq5tI2kEsx1vtaCm6TKYWDeJnr7t2/IQftZP4QBKUVhi4+BulM0mccu2cXpDp0DpW2qYdbqIUEAepi4ZIyg+n1UQF4gRxQvqe50L/uTge9lqOdza6yopTbStczxtbu31/JJ7rIoKyorTCwAXqcsrYSB0CcUCoDNx83wF3/EF1dYlpx8ePzBV1SkUdoNbfRZCrrpNOPMgxqeBW+0BZXUAuKLUwsPUF1XR2Osu+Q/uJqCWLGpEZ6g7kzqh45iMYXOIUqClcfelZv594xHvRN525zhmI0M0Z2OpZpD6QJ0orDFwGSif7i5PAdUVUY6NXx1Y0aQJysw4LCufVRM75CWGH1XXNRpqiWkFeWa2H0goDF7icxilNpHZIHxDJdD1+OHJSabuX1Kaa9E4GsW28vSampaVCuAP7ADh52zAK0AdSIN2u5eaX0QkotTCwSfDo+uLyInFnhygoO6fzZzlCw3kHcoF6SpIll65WsDyFEh0grXQTpSmqLvUBlzIE67CD+nozUFph4LSaKIgLR+25QIzABrUqexJcDm7IMWP6FEjBCZIO5E6721eFc38kt7ki+dpLcWlqe/p+mhIXbYiecBpbTkIlopP7jA2lFQYukA6ocvCRlkFryOQyauXkoULmgYO1WNprCT9QdBGBgRSJSSYqom3gHMzewtQ1pt9t29uPdb6w++2CcEm1tNTYsudAVgqDsO2oj6I0TL4otTCwryZqHR1Z8dRLW3Gzf75Qs2BaPSQNinYYCARs2tWHM/71Qdzy5Gpr3J8a3vcPWi6BSbgUsUhsIJd7jg1gFGiuJKbsc77+MN7yrUescQbq+mMqyGAedooXIE+UWhjEQbexyiV+q/GhG5/Btfcvdo5v68j2Kuj9rW3lAwz07Ai0vo2J0unCIjxWGl1WaM+3Pa3xwrqduKt7nfSMwVi7bT8mTZuFR5ZZ7muIHLthnijXVa99J9e2BuJCEpdqllVMlFoYxH20cH0xOWkCnaktJB9SDEhcodOtYhcr0LquPD9SUuOSHzyJL949P/L8ubXeOUu/fl4+rTMrzRFXYcb8UtEghPsGajjn2ofwh6UJLudJQXW4uoyMgr/Tx4MJpRUGcZe3i9qOdc4gZ7ryxiqHc2Zc4NWfhLAPcWlqG1xrcQfuPb68F7+ZZ7xB1cvDQHhk96nTaqLiwNlNpBDtlk6+/FRkki2HPwG+bd9B/Mdvl+aefWKhV6ROkCNKKwySwsk8LFAnYAYeWboF7/ivRzFTYYbp6SyemygOf3fzHHzKP7M+KVTG1mm70JOsJnJJZ3SnSdNGBWqAGLh+qyJssCsCSi0M4j5mUQ8gc8Uy/yrDRevNy+kixzQ7ab+d1xZpkXyStP1tk3R3fRpdvghuIh2aQYduNZGtrLKOj9IKg7gBoN532mnLComALr+Sac5jUWE9y77N9U6z8ivpQXXS5iNLXu1uCyA5M0oSW5e365LLToWsCHaWdZgnSisMXFBURu+KLgrvLMgVLPqIy6sJqehE69CIVFXRuwkzZFh4SEqhOIdkOcalnC1RemFg/mySNeCoehaVKdosg2jNXFZNJU6SO7KfQMrasD2eS77th047t+6pSRBXjNNu4ZiHqysO6jEknTJX1gyUVhg48XdpkUy8FlBUxTEv0/1QGxRiH7HuNTmE3ETGPNpc7wI0ewPtFpLNQmmFgQvEQdU3kOwi7SIj68mN6uThs6u3ZycqIbIKuHQrSaJLKSPxC8SWxDaK3uecnE51QUW7J5CznhLg3gZ6d3GnLjdPi1hhQES3ENEWIlooPDuaiB4iouX+36OEd1cR0QoiWkZEFwjPzyGiBf67GyhuI0AOsDKEe+4Bf2maR9vzz+Hqmx6Nzy8fsoqDSZPC8KxZ2ihEhOsfXg4AODjYXoGZe/v/9gHgk5/0wkuXYvO9v827hKbCldlZB5rYB+6/XysA2nVybVY3kQvo3nvBH/+EF162DNsfeDgRXWWCi2XwvwAuVJ5NAzCbmScDmO3/BhGdCeBSAGf5aX5IREP8ND8CMBXAZP+fmmeuiBM1/MUvgrd7Gi/192PJnnhG1ynmoe0ydKkKa9aE4X//d2zzr8TkAui/ed5aZvxs1/47eNs2L3zwID69qKZN0wl3WiShS4oq9AG+5hr0+wK/SNUUe3CSMegUddo08NatXnhgAB9fOxqAp/i0fxS0FrHCgJkfA6D6CS4GcKsfvhXAJcLzO5m5n5lXAVgB4FwiGg9gDDM/zd7XvE1I0x4cOBBqPsww6R5FHfx5g/v6w/AhUmf09UnDff/QEY2w+Fy0iIrUNvaTNVNk2NdnzKOdjLGpJUf4gFCu8FM89bSsQiLtnMHxzLwRAPy/x/nPJwAQT9Pq8Z9N8MPqcy2IaCoRdRNRd29vb0oS4zsRJzQ+y9kF9Gi3FSQzuhTnKznQT6IbghlsKOZgTRAGBeoFsgWYc94c5k9E2NuX9U6MfJDEu+xiHTII7JDn4o279RmXCHlPIOta1aR2G5uUmacz8xRmnjJu3LiUhNg/sCgICGZGINOVipSmw7ZhLF1+IcjwvNloHdM1WYR6V1th+0AOdIljQhI0AKbds8ALt+doogaaoaRUSqGHtMJgs+/6gf83OD+3B8CJQryJADb4zydqnrcPh42UNAJTh5CZUud3A+PqycNGhuECcby0q4qcaiDUGQCY9MOhqHturW4iSzrj9x2ptIcmWoG6Riyc7qcYdZglfRjupPtP0iKtMJgJ4Ao/fAWA+4TnlxLRCCI6Bd5E8RzflbSHiN7gryL6iJCmabB2hm99C3zMMV54xEhgyBBz3ILDxhSsffjkk8PwV/9VeqVtunasKGlmoV/9KjiwPIcPBw8bJpTbWUjNrMQ+cM3XwjC3vw20R2PkzZX/4z/CPiB8/0i5KL51mBUuS0vvAPA0gNOJqIeIPgbgmwDOJ6LlAM73f4OZFwG4C8BiAA8CuJKZgyUanwBwE7xJ5ZcAPJBzXRS67e/5fe8Drr3Wi/v6KUBXvFzslE7gTObq1WH4onclT19gOH2ri94Nvv5/AAB05plGhUA+RjsH4loMK8lCH+B3vzs2r3ZryLlbrZdcAr7hBi98xhlyWYZd7EWaN8oTQ+MiMPNlhlfnGeJfB+A6zfNuAK9KRF1GxE8ge7AdUpYkv3ahGZ1Tm2db/MWt0chcT3ctKiNQ6ZWvLk2Rn5Bpu2ps21TnAqcxLbmCLMuwHZ53Okq7Azmu70SO6D1Eths6dXChbahdM8gCUmmj0rc152s6sDLpmf+dBNcqtLuqum/QDJqczjAq7MxRfiitMHBB0hu8isQI5N2ZqkZjJtTt0ncDI2wTe2hVuxdV63dB9DgKfdiGoh3VHUDU2JOwZNf6hNffqs9N8cuJUguD2Mtt/L+uHaxIq2xEJGFinaIVZoVrm7hYSu0+o8cFeXdNRng2Ubv7fbPL11rBlnLb3R7NQnmFQYxvgSF3gk7TAoiasLLCRyvOhEmCZk5aNhQCx3uwi8YImts2xaprYjSJ/IJ1gdxQXmHgAHmQJ4tfBKS5fMP9LPt0NOWFzBvnHPz/Ijp1EUGaSc6k1lC7kYfSY73T4tCYLoxFqYVB7GoiDnyFjpfbFLR3JJncctVstZN3Ba1/WpDlNisRRbMGdEjGL7O50DoFLkqAd+xGBaDEwiB+NVHyAVFUs9l6f3Gq/NqP9q9njw8XCYlOLXVYNssI+0FR65wXGkqh6/LikrZHaYXBoYA05vMzq/QX1UTXqWviJC4tPfI8XymPeFlSNBOmLmC/6jMeJmWpCHNIrnDtQ8ZoDkKzTCi1MHDV/p15aoH6gOfXd/GFy7h7bo82npx3mG8L7iCyIo/indaRKy6jIu8zaJWlUoS65gWrNRQsJLEd69JJUjAlSisMYo+jSDFJ2jFjo0luonb4ztMW6TYXAGNblWpS0cHqU597TcOR5C2tf0br1N06NC3EKK5C0AyUVhjkhaKahHkvLY3OOxSz3nmjqN/XBvGzG7uAdTVRZ9W5mUp5qlVXzSGl7TikhYFZ8h9aGoEK0z6DdlQ/rbxLSqtt92nRJpCz0rNq6z6HQqDXzNtQf7HIZDuQ5Qlxc7yE9BShEzQBpRUGSTpNpx5SZt5nYKaznmIVle1ZK5G2eLPby42xyncgF6sPpMGP//hSovglqLIVQfVsR3pUcwaHCFxNxTIPCrlunV/RxN8qcmKlSdC2Hy5uItnPn2bNlD5VOxSitHz46ZXbYuMwbHMoxe0DzUCphUEck+90n2Azj6MoYhskqe3Wvf2xcSSGl8FiajWS0pDmxE3TarUi1N8V63cccIpnOqjOGL+D2iAJSisMkjBK5wOqshDUQqSZFDPGK5B5nKT9Z85zu1U1rKdiGejiFAy2Lp6G5qLWMw/YT/I1PDenyEpOIVFaYQDEm7Qun1RmCuXsBCqKUMs8D8tL891Mu7iLNG/kvqlKWSnmmHeQf6HqnHNcZqF2jpfblBWlFQaxx1E4yv2yCQDzmmohzJ3LABLlm6aexWmWpkM7Z3AI1V9EdQfyIYpOn0vIcxldu9Dq1RuR4jqgD0h3Myvv0k4am/KLe94qpL3cJk1EF8WpTCi1MIjrDIkZYwf2AtWycXONcaGERrs2HYVxCtQYApp7FIV4uU3zynGixRBOm4f6vJhft/UorTBwOY7C+E5aV54TQS1EmrOF0giNZiLzQXXO/nTTc72W3In9wRWRwwo1rdMOwZhWGUh6253NOpSssJL2gdIKAxckPXukSH70NLjktSc4xStaZ89KTp7VKVIfSLuayMkaSk5OS9AMKzHp6sEi9YE8kUkYENFniWgRES0kojuIaCQRHU1EDxHRcv/vUUL8q4hoBREtI6ILspNvhytT6/TLbQAAM2aE4XPO0UY55vAR7ppdmyubeM5AqP+ogX5g8SK3dA5Hj3TMxkOxD7zpjdoow4e6D/ki1tVKklh/APz8C2HYpOBxmKdrnytiu+SB1MKAiCYA+DSAKcz8KgBDAFwKYBqA2cw8GcBs/zeI6Ez//VkALgTwQyIako18C30xDJ7h6i/WhwuFpcuAqVMbP7lHf0y13TUmh1n3ooVIvLRUqP+Yvj2gBx5IVoblKAKX5+2A9D1Xr5baAD3rtWnef85Ec36mvB2eNx8OvWDqVLxmwzIAwGs3LAXuvdcpZ5e6kuF5mZDVTTQUwGFENBTAKAAbAFwM4Fb//a0ALvHDFwO4k5n7mXkVgBUAzs1YfktR2D7w1JPA/v3aV6IlMHpE02Rv0+BsICj158FBp2SOXmVXKtqH+fOlNmCh5cYdMaIRJnSqm8OB5v378ZfL/wQAOGnnJmBgIFO+YjuJyoLr+V6dhtTCgJnXA/gvAGsBbASwi5l/D+B4Zt7ox9kI4Dg/yQQA64QsevxnERDRVCLqJqLu3t7etCSm7vRJzyppO/bssdY00AYPGz7EefVMUaqahgyG2+YhtX8MHxIOB+nES8k6LEjDQLFm9kdPIn351rUAgCFdQr0AHDbMRSlgbV0H6/WEVKZHq1o68aLC4nSBXJHFTXQUPG3/FAAnABhNRB+2JdE80zYrM09n5inMPGXcuHEpCbS/tm5Pl1aSxK+9bgck8o84Qn4nVJ6IMMTnGmnOHCpSnZsBsR3/+Z2TAQDHjB5ujt9sghJAElKjRkffG5zgo0cMNeQnC0AWnp9zsjf1d8zhIzQpW4Fkk0jESv2tE+qsLcHkIq4sgyjeCWAVM/cy8wCAewC8CcBmIhoPAP7fLX78HgAnCuknwnMrFQKiKS2iE+YM+E1vBkaNcoo7arheK7QtK3z9JI8RTDjysHQEZoAzCxDqT2Bg2DBttMgSWoERBG0zpIvMvKOgfQCvfrXUBkxomA55WDN/NmEsAKBWb10DmE+LMkCtv6EPSLlK80ZuJVbCIIq1AN5ARKPIa8XzACwBMBPAFX6cKwDc54dnAriUiEYQ0SkAJgOYk6H8WMRNmIqDZOpbTgUAHD9mhCVdQTvB6acDP5ke/p4YThKOVpj/6OF6rVCF2Aavn3Q0AGDEsDavRLZ90Olh/UeiDlx4YabspHhSuDh9QDrKetIkqQ0wYSIwfnwkzQjH1UTqAotgrAzUCuwmmj4dfJS/eHH0aOCSS+zxg3ISLrVtoTxsKbLMGTwD4G4AzwFY4Oc1HcA3AZxPRMsBnO//BjMvAnAXgMUAHgRwJTPXMlFvQdKViV1d8SmKpBCIpBAB/KEPhe+6n82Wt8FHXqT6R3D55Zj6Vk+gn/O21wFnnhW+s1h3phU0nba0lJmByy8PHzz5BDBmTOPnZeeeBAAYM9KsLatV092B3D44jOjLLwf+5V+88HsvBr/6NY1X9rPITG6i4s8b5Qk3NdEAZv4agK8pj/vhWQm6+NcBuC5Lma2AusxSF+4UyMzOVgO2/PLQdvM4ZiH44QZfeJ5odxO4Qj111W3SWJ/e5Xkz0IrjqVQL6NRxo7Gydx9GDO0y1rWyDDoQVtbHbtK+EzTjOLrkebT4SpiOPG67MIhBUnNfZQTS8w47pKwZFx3p26Z1LZBqJZkiBJPifa/1Fjiecqw8IS/Wu+jjIC1KKwySjA33y22K0wmsZwnlQGbjkDLNs1ZA29ZNICDIMcJMJaYiKgTF6QMitLeS6eJZ8xDjhbVOswqtCMhyPpUtbb2kpkFphQGA2B7s0gnaxQyTwPWMmiQnURRJ8GWFKzM3vcmqbeYN2woYfYJsZbg8bwbS2Dou31B9HrzL++a4TkNphUH8eUNuTKET1hfbBFhc3CTPW8sINN+PyFkzd3KHsX5jlZe+M2BrD+kIhaBGCT5iu7u7mbHHE5bEM6A9nRXmjZed0jeSorTCIEBRzfo8kFRzYrh1ZFYitqMJ07iJ8rJmIsK1zW2hwskYsLh2zHNlrAnJaYpQfxcksW4blgHk/SWdfpR9UpRaGMTegZxBS243VMZum+do9a1heYIo38lRPevTxDO+bH8nSOomUseBe3OyNn2rYCLTun/IILiNgjFF3YrEB/JEaYVB0OETM3zLk6L60aMT4Pp4cZvwxDBrnrcDnk83m4bmwjwIdoGapfwiQLcgQHweDevdJK2sPlt+xSGRm0jsBI1n0Qn1sqO0wsAFxg9sHCDNpSctmkGXdnVKURtAAYONc0buK0yKO5fgNmecL6XtZIaux4tLcXIg1zR3WFaUWhgkWRLntEa9QB2C2W3OIOr/dptYlX63mRG4X90Zhh02lFv8ymzRmJ1IaSqS0mCj38VV5rIfp5nwinTsA0KNTN0mciBfwu9bViuhtMIg9g5k2FxF+om0onUBE21mrdbcjU2TptLASUxh/nChgUA48Wj9wX1ZB3LZGIHp+3b68mIieUVauiXFnVv/NCitMAiQ9IOqmkKWvCrkC2crD4zTxh3ulC7p8sGidgHtXICB1neecXyzyckVrptC5TjmKz5VS0AzZWBdeVfUPpAVzT/MpY1IwwRs8YrWB5zcRGLY2V+ur2vRB0FyX7JbhYrWB5IurlInQoOjFsaPHWldZSNNOmtcRq1CkjLFuH9+6kJCQgQAACAASURBVNFuaXTLRlRXaRE+fJNRWssgMBGtAkEKm1xDBeMEDVjMXvbO5A/hdsmHmkcjekHq3YwVsp3oInAmLYPQAArz2SNwchWSe/UblgEpfCCh1djpKK0waAaK5kNNq/Vrnzuov0WovxOTZndmLjIC9Vnq8gsAkxB3nTB1mU9rFdJsM1EXj9jm11xmE4rQ95uNUgsD212+DNPkgOWIgg7sDyw4j5OQX8TOn5Yi1+9p3rBUrLaQLrVxiM/IPnfWzhawfadI3JzK0c7BNMLF6g95obTCINmmE7d3ReoCtqWljGwuFZNrrPBjIAWBbmcYZSoid7gtg3a0jCS3CGvfFa3+WRFpmxIpf1lQWmEQIOn5JF4aw/NO7xzMGHuY4W7gaNTIu3ZXnxxpsM4TKfUKV5J08JkdCaD27U7wi5O6zMcBtv6sPhf7t6MXspQotTCIYwpJP2rR3AUiTBqeSvEVbzoZgP1y+yIIPR0NydxcGco2lO+9K0Dj+HBfESdT7VwHjSLUjtqzMgFipV/ZUOaev+a5IVxWlFoYAK6SXu5ARWCGeUIcyEO6vE8+dEjyHZ3tapesOruNEerOpbGmL3DfcNkg6OpvlwWI2ZXUcuRgwEWsBmgknTpnUlR/cY4otTCI67dpBEVR4EyKxTw2/U4z4Zg3yLAa1oWZJVlJJMoCNmif7W6LtMhqHRUFqpsozzkTL67mmXPq8qC0wiDuPBvbmmqTAOiUDsKcbjmeKa9wYLWuBZrFf7OupS+SXHCeD7O9s8yJuV6l2QpkXjptUepM7iCnvUclQmmFgQvK9FHdtCU3jSnNfErzkHxqN9HcQmAZSPsMuPOUAIulJ/42WT+6/EJLS1NAu5CoMySjWbUGiyT4W4FMwoCIjiSiu4loKREtIaI3EtHRRPQQES33/x4lxL+KiFYQ0TIiuiA7+XaoE4Eu682t/b5AvcNdU0rXqbWmc1uqr/j406wOM7rKuHGVqSpyzBpmcZFmwjSxQtTuBshoIUTiObud9OEyIatl8F0ADzLzKwG8BsASANMAzGbmyQBm+79BRGcCuBTAWQAuBPBDIhqSsXwj4m9ATs5UWoYZM8LwpEmN4Iub92TK1r7X0uQvz1RkOsyYAfz+d1744x8HP/gggPj6p1kWGPwmKtggN/QBI1wvcxfCUQEo94G2tseMGcC993rhz30OuP9+p2RGmi0CULufIpJXkTpHc5BaGBDRGABvBXAzADDzQWbeCeBiALf60W4FcIkfvhjAnczcz8yrAKwAcG7a8l1h9RUapH1btcIZM4CpU8Pfa9Y0ggcO1nIrJvZKUE27tbT++w94v7f0Av/xDQBA30DdKYtETIyjjMD7bZg3agWHtPQBESYmb3gQTR/D8KXs29EH9u71yty6DXzNNY3y8/wELPrDAKk/GI+CL6lcyGIZnAqgF8BPieh5IrqJiEYDOJ6ZNwKA//c4P/4EAOuE9D3+s+bB8tWkzSZF+rhXXw3s3699xUhu1qcxb5mBervaRK0/AXSgr/Ez67dSl1/qsos+a3FjWPrAS717pd+mVVcqkt7c1la3iF//Onnsaf3Y44A+rw+s2rrPLQ+L+8c0gVwkNtAOZBEGQwG8DsCPmPlsAPvgu4QM0Bmy2vYnoqlE1E1E3b29vamI092BnMSPGKZpsctk7VqnaFaXT6SzJ/efBEyz5Uss/fr3jvammm45572Zs3Rpq6jrwJRXC2DpA7sPDGDbvn4AntvM7ZPEC0Almq+Ba8RiC/vAwpedBgD4zls+LNNg1Njzo021mjpuQUEKZBEGPQB6mPkZ//fd8ITDZiIaDwD+3y1C/BOF9BMBbNBlzMzTmXkKM08ZN25cBhLt0E+SslEAtGQgnHSS8ZUoTQ8Ohi6T/QcHUxVlvgAeqLt5ZPKHX//1Y7zvPuekP3NK1j9YC32/cNOEVVeAi5bcEk5g6QMMYPNuTxjc/MSqxvNZCzbK8ZytQLeILXUTaerPDsuItu49mMoCaig+BtegLX2ZkFoYMPMmAOuI6HT/0XkAFgOYCeAK/9kVAO7zwzMBXEpEI4joFACTAcxJW74TjXA3nQuzzPS664BRpisbgXXbPV/6DbNXNJ7f1d3jZM2oNbT7i6NMsiUtpKv/yJGtKFlCWzVBSx8QEbefRMfcEm3IC8KtboEc+kBaml1SdeomxDhkvensUwBmENFwACsB/D08AXMXEX0MwFoAHwAAZl5ERHfBExiDAK5k5vxmRBVofVJGTVhxERgYQUu6wOWXe38X+L9PPlkbbef+g9ZsyOFYt8F6lOF7P4Rbrlrd7/368xzhC375y8BL9mRpllWKv1195C1hjJY+IPbrGqv9Vv89bb5za//WWs76vHJFUP8nwrkiXPM1YLk9maurL5JOY/W4LjktEzItLWXmF3x3zquZ+RJm3sHM25j5PGae7P/dLsS/jplPY+bTmfmB7OQ70Wh5mehx6xAMBgBYvVobxbuVyQ0iwxObo25pG92blg2Oyy8Hjjsu/H3hhYmSR48oNsSzzyYkKjN3GPqASNWKLfJkcp4wLb9sGS6/HDT55WH57363NlrdcaWDzf/vZg3IacqIUu9AtvEEVSvSxYnLr53o3dNvfOfacY0aozF+6xog6a7jrPsEGBaG4Wg1dAJEH3nw2xg3mIMxuCDbBZGevYb5sjTfybbwoqyuIRGlFQbB2UR2Zhh2dtNEcVE7wept8tJDp0taEgxk0XRuRxOYBZIYtjA1Ywbm/Iy0FKgPJN0XEDxPrP12iJtEndzPKgTEuRXXpalp0LNjPz580zPY0zcQeffL7nWYu2ZH9kISorTCQIc0fuWWzxm0CJGVNCoj0NQ2S/3/7uZncP987eKx3OFMp6Getjzy6AP1OmPeup055JQcLvMClsepGWGtznjV136Hu7rXxUf2kfnocot7MIwTOgutimNGWgBg5rwN+Mgt3pqZ6x9ejidWbMUDCzdF4n3h7vn4mx89lUOJyVBqYWDbYdn27fY5wLY0VIyTqprCPEMeeHz5Vnzy9udTpTWtmHEW7lI4nvtFr3/MFzc9sRIX/+BJPPXS1hSp86XGdb9K1gUFBwZq2Ns/iH+bucg5jXx4oCGSo/aeeLMm1HGUKLkWn77jeTz2YnTf1PZ9BzFp2izc98L67IVkQGmFgetF4eJfIMYM7yTpYdji57KJzihA2lR9V591FqGn5t1MF8HSTd4ZSxt29sXETA+nuljSwFFjTos/vtjbPuvI5Cq0rM5SUuRIDLDS31V+29P6Y0dahdIKA8D/tikmv/KeMPzpk6swadqszGcL2YSRs8/cEMX1ELtmYl//IOas2h4f0YeRKqtAl6MlrVnatvj275Zi0rRZCi3RvG56fGWiNpDzi39nZIQtQlDcFbfMwcU/eBIAsGVPH1ZscT+EMU8WbRKIWS2LOCRZ+t4qlFYYxF3uYpsgU2+MzYrpj60EAOyI2Rsw+erf4vN3zctcnnG9eao89HnlhR37DmK1f97MZ37xAj74k6etK6VckZXJ530MyQ8eCTdK2G5o+PqsJfjgT542vk9CS+I2UAWlYe4sb7z5m3/AO7/zmDVOUpeg5wZO7hoKw2b2nzfPzusiqqworTAIYL/hKGFn0USfu2YHNuw84EiLd9DWpGmzsGLLXjAzrpu1uGEmDtQYv3qux5o+jp4GHPyttjmHVuko7/zOH/H2/3oUALBk424AQN+AbEGZfMdJj1IAoi5B021eJtdUp3gKjfRHhJ45vc570gztdaCmz9PlWqNU1r7y3OWek1YMiHZ3rVILg7hleBrlN7KkLK7v/82PnsJbvvWINY7YpX8zz1tRc98L67Fm237c+Pgq/OOt3fZCHGmJT+8+aa6L1wxGsG1f1FpK68oIl8O60ym6T9rF6Jds3I3BmuPx3MpvkWEWYQ+ACa7KLzOjf9BTBlqtMbuOe/XVQK2OSdNm4abHVyYvs0DfrMTCwN9n4MxI4l0rpg9XS3nec5DKthNYhPUcGlPYsVObBkIzmOS0X83HtfcvjjxvnDTrqu3lSZSat9GFmGMZ8CYP3/Xdx/GfDy7NMWdTaX7IoOnL/UGvFDWbdd34+Eqc/pUHsXWvbVNlvCYvWjZx6XXjgJU3tnof8C3Z7z4cc16GgKK4hkSUWBiUD82e1JLzdns2Z9V2XPyDJ6VTVAPM79kZcfkAwJ3PrpNO3IyDy8BJJ/QdW7CJPuKtez3L6IWcV9bYhLsUz56L/3/rtNdfP+9Zzpt2ySutstJgTm1xoRndq0nLiEMxJEOphYHN9+1pPvrO7uRDTEMPMzbvzrKcMLl/1EvFkef2CXR3l9BV98zHvHU7sXa7fOnIlj19eO/3n8QX7p7vlI8rXKyeNNYMKz+MR1snyHjh+l2YNX9j5LneBeeWZ6tcWba5BBV3PbsOU77+cFtWw+SxOsrle9jytbHyn/9pDb7kOAbaPR+V9dTSwkKnTbo0tipA8vhAjaMxGJjxjHdxR6ARBmU2D25+IjWWZjolUVvs6/csggU97hqvacLQfSWJQYhbzP1mDsD3fO8JAMC7X60/ZA3I111gmjRW34nPbJPD2jSaNp52z3zU2bsdb4jDKr5WwFZMtK9HlSU1nstx9zph+JVfLwQA/Of7X21IUxyXUaktA8Cu4bppmW5+w6Sopbg9xr6OXH5pulbOSSDaHK7GNIbnybKx5hXNWz9Ak2jZISMwM8K4JZZf+fUCPLJsi+aNrkyWBn+w98S1ndJo32m/gYkZmtPIcbbs6cOkabNw99weZ4a32F9RpqMntnxHN48tTnNGexQuK6VajVILg+hEkfypxZUkcjpTfu5lL+jZhZ/8MeYQ/oww0+PmYU2+Rd89ftDVk7SZOO/gIpzToBma6c//tBZ//9Nnk5fPaEwcz12zQzvv0gw4abmGH7r2IwOnX9nruQ7v6l6Xa7unE25mZcHFaoqLnxT7NXNpbfYSldhN5P+1m4tGz6gmlAx/9X3PRfB/33aavgT2lqQB2VcjWeOofi9jPL3gdD3u2Ia12/SXu6vYlHA+xXUpoDG9IZ3JjZC2HBt27g9PrexeE7/z2N2CMDE5N0s5Lm4iMPDwks0AwpU3LjAJGSVr5bfe5eOSR1RIuH13FzpVBEvMAdmKd92z1AyU2jJQ4dRBbK6YlKJhveED//hRz3Lo2eG+aU0EOb5rxHG0gOLeSeU4DIQ5q5Mdr8BQNpolSQj4q4T0A9n1lFLzapyMlkmm1CnLVJaJyi8taTTxktAv9oyXmngRTxa4Ngtb3sU9dy6fGds1+25ahdILA9cVJ6y8i0ufFQxgd194MYcLk3FdFmhl8o7pxfXW4aord2SdFMt/+Z5Lma1h1WLTiO3ktONWITFpO6v2X/I0KeMJhLq2c9IuZFpFF3lnUIpUS9G13lm9/6JjoJ2TyaUVBk5r05W/6nMv3BwXgZrXqq379BGzliOF3SuQ19pyBqcY1MnLTutHFi84kvNzsyDSlLl62z5rGdb0KSZJozSI+cHww+xLV+HyfcU4Kb2iCZHRglOsqWZeeHWjf3bZvJ5dueabFKUVBgFsywqNadpgy7vuQhbhvurGlF7fNi5n/sdB1HKbpe2YLCVruzhYV82+z+DZ1Tv8cpKntVkPyVfPkPMkqZOH1dhtkos9VxdmOsVBVfAChcDMK5rJEnp2hnNq7VxlVGphYO8nrb/cJrOJb5v8czBvE02qsfw3Nr571pmQZpVROqtBn3fWPqMml9xETtasnINrH9KzY1Xo6YWg3LatGzROl9uoyGgpyVklyyxt3yjKMtPSCoNGA1u1iFAjME4yZtRCTGjVFn/X+QM1jW74p3Jr5FBN8+mqbsJRTiOHTf7iVoCRghE4uHJsz6PZmdtWe6KrtT/pX1KL2J30bRnmDm7NQ3AbOqbPYvkyV5vOWoI4ZuEykJrJHxL70iO/zdSZr4p0cwFlZYxJNV4LKYnTJEnf6juQ1cUJmSfaJTeooS4mRadFKNLJnEDU6nNagqt8N22+rhPylsHV0RPIRDSEiJ4novv930cT0UNEtNz/e5QQ9yoiWkFEy4jogqxl2+ny/iZlElZ3c5o+PWNGmP5tb7VETL6SxDVe0rOWOPK/OX4s1Xv3gj7z2fC30B6pYBmUaVZkmYrIdee5WOdXTBbyTePvluGqc5uXysaHvd/R/mArAxCWHW/pBX/3u+GL22+PJxgWxmhQ1tJa7kmPo0gFsQ+ccoqUr2l1WauRh2XwzwCWCL+nAZjNzJMBzPZ/g4jOBHApgLMAXAjgh0Q0JIfyU0N0Eeie54KpU8Pw+nwvvHZx5bb7Okvetg3YKlz8LrZHi2CzhhqvHP1pqRiOUGdeu9YYLQ0fsLp5NO9cJ2bF33E1jmVgK5aDdoZnVPHHPx6TQENbUv+94zubBZjXXAIAuQ+sCfsAM6QG7NgJZCKaCODdAG4SHl8M4FY/fCuAS4TndzJzPzOvArACwLlZyo+D1RUEs4lt3rCUAvsNu2+b6Ct0OUjLelSH0m6uTMFQKEjMzNQechKFNkM8Qxp1zkNOYxOO8WWmglBnUQ9kVvcc5GcdSmks7qOk7rBUwrBWByFZH0gL0zd0EYIMc4MkX6mlwFDnIjnQsloG1wP4IgDxUJXjmXkjAPh/j/OfTwCwTojX4z+LgIimElE3EXX39vamIiztqaWRNAX6XGlpMTI5B2U4+6R5Wonn4vLJRpuoO9s1yeb0AbvLxazxm9OkoMG0vJhVBSk9mGQh6IqkWnKuFj2b2zrzijJlmVRB5o/TCwMieg+ALcw81zWJ5pm2WZl5OjNPYeYp48aNS0uiX0CywWPbup9OmIhh0j53zivVYNfTYtNCzStJ0jKi/Eap+Z6B7O2T5bkzlHZPeq5N6nkj5a8+jWItsvw3Pr0jbSlM4uTfydZXFStYk5ea2nWuKSnUNO2cM8hyUN2bAbyXiC4CMBLAGCL6OYDNRDSemTcS0XgAwdm+PQBOFNJPBLABTUTk4xhcCUn8qIlpGDXa+C5PzUd188TlHTl7JpJfTlB796hRrgljY5jbQ8/U7GnMcTILAFOdlYyT7JpPAuvEqImZJijJ62uW+F1DZDeRcx/QI4+d29I7B0Hj6nIyQqhz5ryahNSWATNfxcwTmXkSvInhPzDzhwHMBHCFH+0KAPf54ZkALiWiEUR0CoDJAOakpjwGATNMq0k1nkvhFF/uJz8JwxNkr1jWDUdZ4zqvJPHzSlL7xmquY44Bjj02fDF9ukNqlkozfsOsgprDTKLzS8ndNEYIdeaTTpJKyHpUh5FhGdNb8rbE1VkJ0fQG4fbylwNjjwxf/OhH5kyEvLKMCYbbmI6kCcIxypKUziWi2O9PPllIy+iSzqdqH5qxz+CbAM4nouUAzvd/g5kXAbgLwGIADwK4kpndz7PNAZnvOE3jhvjQh8Lwo48KeWXfAe20lBJuAzmat9uz2AE7ahToe98Lf19+uTsRMXTZ1tibfe7m/IzlZNXexDovCS++Z4s0aJqbyhHinIGsyabQyseNAz77mfD3pZc5pbOdyhvGkSfkk9JnWnWVO4Q+wCtWCOXL6FQ3UQPM/CiAR/3wNgDnGeJdB+C6PMp0pMu6KkLb2ZFS+7OUk19e7vGMbiIDMzVFK5IZq8K0gqiZCwUyC3DJ4snX3WFME7F6gryi8bIg777ikp3LpjvXgtTJdNPqO2MpBR4rLijvDmQXCevw8bJuPrJpVVl35rporzL97n7TrEIszYUfJlqM8Qxhe96y+0nPJDlV3m7lh+HVyqU/Tl3WWSEwO0PMZy3FW87aXDUbPPNAFiXZNu/jOm/Uij7ADHQV5DyK8goDH66HjIlapevqCbfyzeVRQl9hmpNNbcg6n+KUNpJXNpFqs/Ti4uhzjheWee5GFtP/71OrM5/ZL79zc4MGsJWtWtVp6m3eQJytFdlAmGmuy56ZMPbh1j/TLPs1ITqmO3TTWdFh+zQsvlcHkTFRZ9iB5tMno8+06S3WhOtAIEPYBam0fOc0bs/b7rNP+DwuL92RC+Jz9Z2U3rHj5OkSHay79bs0cFpcIUqJ2LjZys9yjleeKK0w0F3I7qKlRjpKxj5oWhefUUHW/cxETyRem+WekymvxE9q0bFjOXlC7XOii8Bp9YwjwWmEY6R/OmbidrlNMi530+MrnRokOqYNz215GBLlehKBsewcM8uI0gqDOCQxtxvPM5YjhmfO25B4gNjcRBFjMyZrglmzFleSRMpJ2AhpOvus+Rst+ZkGKGvjxNGgtZrY0jZZXRzWt8mZn3vemvgxH0fXpmnKFxdreL/jadu6J/ldwOmtJs3zyMMcLRMhq4O1ujlii1FuYRAx/S3uoKCzW6zDdkvxrD5+NoRtedjaAzBrfOYjtGMKBvDCup1OF5uY7p1whfOa/QJpb1G/uEsamRGb+oGx/gnoE5HpnH+wseGNNLsqAVZFMIySNH2afvLt3y2TbwVMnkVuKK0wSHLol8nU9N5l1ATtnDRZXknyFuKkO8I6+Tt9fLfNQ2oZTgMrB9O9FZaOLX3ylSSObiIHP4n1SBIDL057JEnW3fZmvmyybt3byYkPNFMhkOYMqgnkpkBlKqq5rzUP1d85CgbbQHYTXqplY3aN6NMnyLuN2jAzY8JRhwEAjhw1TH5nCguuLaswizAZfRo5Xo6NoWSV9BIi9fGQGIYeLV7uM9bFApqwjqoGCRmbaef+0DXkmpWz0LAsNtCmseTruqii01BaYZBNvuboH8w1L2D4kGyfzNmaCJikOigcO75s+ib/Gm849RgAwD+95VRjHKugdhDiKbwITYXr2URv9NsGAIYP9frDMaOHp7L0pHimNovJoHH0S8aW27ZPEAbsllue7jxVMDarH6h0Zll5lydKKwwC2DdahQyPwwS5Liu0m87J8zrn5KOk37FpUnp/k5rItjZLM6xcjiJIOimpphEzj66rT553mvKTWoRHjRqOUcOjd0KpFhRs7jlHPp/4DuQcXWgMYPTwodp3chpz/zZbCWIc1TpySJNR6KqoNp21AFbXB7t9vMwd3JJXcv+gsirDFIsdmWlCF1oQb0HPLtw/3ztwtnEgnTMzzg9Nm/BN4E5LnLVFK3RBYAlE8o0pizXPbEhSTd0Vs1nBDLzxNM8C+uibJinv9ILaZdWYtUyl/MTnHLnOU0QUgkTFNA2lFQZJmVQA21EEaTp73aJ1JDUP1bqMEBiDkx9Uyct5842Cv/r+E/jk7c/HxgvocHETSf5iaYCbBbo8eN0O/lOZhzGJRdBs3HUAk6bNwsOLN8cXGAPTLnSzcHNTSy0KszGetQ8FFpQf+KvvPYEzvvogAHFPj76H1xlOq8OkOMLzkcPcbsdNK4xsik8YjlfCnMuLKIX6cKtRWmGgg9FFEJkQyrHMHPM65vAR0u+J/iTr6ccf4VS+lfmpaSxuHxGBdWNzE7l08O2ivxihaZOHJptk4jwo31bsgp5dAIA7nzXfZ2wsL3EKu9XjujtYm5fyJHpQWxCOZrxg/S4cGPAOHW70ASWOaDFkcjVaKqbuvWnQbLUSFMtCSGNaWWRD0m+q0px1fi0vlFoYWHia/y5qbtoGWyrGbkmTdMfpkC59ZxkbWXETf5UekfuAkRiGauI6pY+HyvDCSUmboDFbI/mePCvnHScA7Xllo8tFi40r1+3MHkPbah4HfSDK5NIj2k7x9KRd7GAq35R3VgWvnl/XzBWlFQbajmgzpaXnepOQwegbqOHc6x7GI0u3RBNqizQLE5fBIp7Rorp2nDulRtBFB44c3UXDSlR4o5zkg1rOTV//dHMWCls0+aIdBWDy8tXjKPQ5uygktg2VJuZPlvxMeWijxy0tZVZcQG4fy8makMLmBNa+36ibmwByoSVJxGrOoEVIftmF/aP27NiPLXv6ce2sxdZ8ks5ZmDrEvHU7Jdq0YBuTMJdpI80yriXE1dO9+fWa613d6xxTp5NaOqHn6nLKQ8FzXU4atrPZGnKBdY5Afec/iPuGJlnQZXAfWajT0mLw7vrx3AS4rURd/WxigcHYsPMATrlqFhau36Wn2QLVgkp6PlWzUGphEDeoE2u/DLjqhboBkupiEjHMejeR8w1uCsM1LsuzTHqrCAe8mUklnRwX5xl6dhxwThOGLVqyo9vP9q3iGLMNNuuQDPG87x7QpeQnhZNr3IlroElgcpuJioJbU4Ut4Hpcuzo+ssBmeavfevbSLWAGbp+TfN4o76Po80JphYFuUis6kPTaqJOWGKctNQaI24c3TRzVJTeRfmozSd9y1eBdmabJMggFbfJBHaVNz73yWOHBjb/mvFSNM8tSysj3S6gKun8/Oayrp1qyOgaSMlqzQsBOq4nUFnU7XttMg+vQ1VoGbEkvmBNdMXNvOqhzBpWbqM2wXWRinjx1/3BJd+gbjyUwhL3f8bmzIR6BrIMlKaNLKTPDeCk4q9UaMFpHah565mHTnpOfs+OoEBj7AAvKTX5apcDTtGVqXWi6vmRQCLoEiyHNBkuRTktMLW3RepnHuzmNGbV6IAySc/LIjYeJc2gOSi0MxA4N2LTXNHk7xnOxMlLkpdKiakS6ycgkPtVGXIvGCdjWmOvjm4VGci1fFZRp5oe0zy35pJq0tkj0LgdOsGFnX8gwLMw7DfLYaR2n+NQtWrZEi0EAWPeapPoeZgvEVI4cP9TuUwkD5Xc7D6cTUVph4OSnNj0zdrb4JZuN8huuBPNgS3KyqpqXHCf6POnBdapW7MwXTCtgEtCh0mIr3MQk7PND5ngc81xNnwa2qslzBiTEC2Pu2H/Q6Joyu08MwtXybaOuJV2/MiRGtJ1crF1bXkF2Nz6+yhzPmolDObAIF8vYDfz+aYRBnktw80RphUEIsxlp7gQOucYwiMbAtkRzE1gOTF39LTE/mSkk1YLUd6ZllmaN33U/heWd0/OUGq4mrs3KYCD50eMZhUm9Hl+3uElaV5+76Z0tXpdpAjnYK+KoRKkKiemdCGmHKvF/AgAAIABJREFUv4XeNJ/AuB+AWRAG2fMtiGGQXhgQ0YlE9AgRLSGiRUT0z/7zo4noISJa7v89SkhzFRGtIKJlRHRBHhWwwdlstKRR39FvfuP9WLECmDTJmGfwgTNvMLGolTY3l2h1uE7YhU9DzhIRHr/4RRieNAm0datfhsWmdoBpMt8VmeccXPNgBh55xAs++CDY0gfEvHVlAvJ3MjGFWp1j5ynS9HXbOzmsWBkzZoQvJ00C9fdpaUh7JIxapjVeCqZvtg4V9YINffLZZ1H/xn8AALpuvBH45S+daA1QtzCFrK7kLMhiGQwC+DwznwHgDQCuJKIzAUwDMJuZJwOY7f+G/+5SAGcBuBDAD4nI7dCRNHDqiPFashResgT4ylfClGvWmIvXuYlSMLm6wrCMyMgMXbPjT30q/LFmDWjlS9a8olZLPC02LddmzdmYrokikyvEmNfChaDvftd7TmTtA3o6k7sIaiam5D/RlmkoP0KPFFaZoSbfeh2YOjV8sGYNsGu3ljbTWUOuGw9dZXNYRiqp4xZNjHf77eCd3v4f2rUT+MxnkpdbQKQWBsy8kZmf88N7ACwBMAHAxQBu9aPdCuASP3wxgDuZuZ+ZVwFYAeDctOW70Wj+bdMEjZ3qySdBB8zr3l1oyAKba0uKZwpbmX88IwAA7FfqX9ff4Zr8xEeRFudUhvQqw88GSaA/8giovz91ehUu14PW65CUm1YojxGBGPyqM7B/fzQyUownSzybEuG0JBxxJ0xF83ZWXA4OoE4e6+xiTXvEILrPoBh+olzmDIhoEoCzATwD4Hhm3gh4AgPAcX60CQDE7aQ9/jNdflOJqJuIunt7e9PRJJxt4wIXpsl79wrx7R9QV35ac1lnZchxUuRr0X7FyUN1IpEN3CvpjV122twEUxqtMMIwAkbg6qTetUuJlnARgJKvyzJV2TIIlzXKZXDke+omka235Zmei/RryA0WKUdZXLIznES9wuZKEeGqRNhzM40rc6pAGAzhmjVnbVq1D6Q4qqMZyCwMiOhwAL8C8Blm3m2LqnmmrTkzT2fmKcw8Zdy4cVlJNBZmHiBiWPh1xBGW1fkydLtU00wUMYvLNw15sL0LhZoPOw9Mdx+0vVKucyZRbT4+YZT5uZUl5eEgUKQ4Y8eCAk04h2WFIsxzBnVpCa/uchs7I0zetkFZLmi0hyF+ndmJaYsacz3FB02leBiEHiALJ1kgEur+x+pKUabppNV2I5MwIKJh8ATBDGa+x3+8mYjG++/HAwhOdOsBcKKQfCKADVnKd4FZq1SPY5Df6dLzm94MjPSOjY5jBCIDj41kAUPeTX3CkV75xx0xIhRmahpZmumLJpURymFjux02Sn7QZbhsxaINh8/FwQ9tOJLGpLE6lOelN6RxLAdvfwdoxHA/jpswsPnlXbKo1eXjoI8a7ZX/tleMs9ZVu9fEQpvRvyjEY+oCRsl9wHSEdVItV7R4XFO6WI1x75K6WzF8WMMyIGaw0h76MlkbLhKyrCYiADcDWMLM3xFezQRwhR++AsB9wvNLiWgEEZ0CYDKAOWnLj6cv+izz0QivPB30jevC3yefbCnfvukLUNaVOzBMZsZQfy3b6ycdbU2r5TEKw5QZnt46igzwG74b/jj5ZOBU8x3FXnp1P7deAIjln3xM/OBS87LNEySZK2nQpp8KAb/qLNBnPxs+sPSBuDIA4KJXvSw2fU1ZRRDUZ5hwH3bUy8WJGZCzG3L69DB88snAkUdqaTCmNyBiGYh5mBQ313kGBwXBtc340svAfp27xo4B/ud6c+YxtKh3mrdTTmSxDN4M4O8A/AURveD/uwjANwGcT0TLAZzv/wYzLwJwF4DFAB4EcCVzCodbQliPtRXDDto0M4D3XuyFTzkFWL06vnwlL/Eyc3MaWSuNO/Qubz+jVXP64AfD8OrVwLhjY/NyEnTC8yknH+XEQVwHuymNLZ21O/zFed7z884Dr1yZqEy1vDGHhXdRGM+nYkbcpfPRgwfNZcY99/KzvL/88jC8ejWooRknK1+FKACGDekyd4EUFqGc3GD1mYuR23bKFNQ/9WkAQNfnPgd84P0OZYaoM/CO0z3398feckph9hkMjY+iBzM/AbORe54hzXUArtO9awasxwrAxgjM2kbys4nCHJZv2YPxY0e6ZSBkRBppEJkglAiVBYVxktBqSrP0Ny6Njam6aJxDuwhHjhqGnfsHJI0XMA9enTZso8VWvqhxe3mYGVue9z6LeVj3GWjLdGNqanmu7jVXJSPusEIdDTrULbSYIFsQ5nFrQ9I9MoxwUr/mSKhq9Zw67nA8sqwXw4Z0YYLv+m03SrsDWeezj3Z2PeSJo5TajmaA3PrUGkw82tOiLjv3JIUWUzmyVtjw3abSiFzjuQ+lBj2WQo8foxeA4gAZ0kU49djR/nPHslMwQw15sc9VCyb5oWth+hc378EZ48cAAKa+1e5iC+BZBkFeIXWeq89UaHy+qvAxasIx74JsTN+NmTHuCOHKVoPQPUy46zjtEdbZbjeTf9vuL5+7ZgcA4PqHlyfOu64M3teeeKSxnFaitMIgFimERJKOpnPtqBdfv+bEsbH5OFkvaRmhrQ1cuKRDfAZL57fYtNIsJ3NG90mYNMZslg6Axsf1VsnIEbtXb8fiDfKiOjHGUy9ta7gKxYvex4w0G+lDhDMP9LWKsT6E9xHmaXHt6BQP3e+471ZnYPyRoUJgnjfSl6/SbaPF9M5mXTq1qZJm655+6beIHfsOYl//oEyLqXyY263VKLUwYMR1Fj2TSKtdPLViKx5evBmA/j4DQkgQM3DMaE9bOsZfHRKWH4atg8VABwvqq2mS1gqBSYiMxJokhQUVEQZ+2HbKpcrMtGHHwWU8zE15opY5xP+2Ou31/T9+Ghfd8LhSThiWGLvw4lhRc1bSHHv4CG1/smnvNqRhrE75MvDE8q2YNG0W1m4LN2J5x5uY0oiCgTHWn0M574zj9AmiOWhC6o+YHAKhZ3iug81yOfvah/DO7/zRmFcaq74VKK0wyDIpY2Q4gHb/QIAP3fQM/vG2bql8m+ks5mvWUFlyOWk1PEQZoOhW0OVdV7RCZ43Zov0BQP9gTUrPwv+2vMTNdc6mmhLFZWBF45iZvqn8gKHrNn8ZifOhWkmsiaOiXmfJ0tQxL12TuQlxMY3SJ43xZIhzBr96rgcA8Ozq7U5Wl9iENWacMf4IAMBJR5tXlLkoBDa48l8Tw2fEuzI37uqT00h1lgVg8PNlY0ZmcnNlRWmFQQDbdnVjGkfmY4POsBeXm9qW0YkYNrRLv5rIMtpZKCvi92z8Nm9Ac2UkIk0MxgvrduL0rzyIR5dtCd9b6qbuLUi8a1xiBOkGUegKcdOyGRwKA2fmE0Y85+SjGr/Vb2OqQk1RCETaw98WYlJaEyz0FRcw2NxXDWWqFrl4AmrwbsTQLmtfjdKriWd4VeeQAlc6g3RJILvGhLk/Ic6IYe1lx6UWBnHfK/wgakSDRiAMvri8dffCyhNkYbm2TnjM6OEGa8Q8wIHwaF23mmloYANtZi6J7tXbAQCPvbg1tjyxjEY4sKbqjhsCXbT6SPoI2bHx1P4RWgb1xNr36OFyHwjD5pxkC0TV3oPfMaawi4Jj0LJNfSOAaBno3VluZdZZtKjDF2eeMMaNZoTVFJm8DXW20CopbOJjdj4uQ5evyfXb7ruRSysMXM5Fqas9qfFcfGxgSjFdTdSQgqVjf/26CVJnbeTlqD2z4bkOXZpBCYRnvtiYvO6cGyMCoSE8IopvHyBqGegE2FtfYT6OxMqwDPGiecR/X1UGB21bM2xMi+aldxuq2qKpfNPxBapy4OIyYeGldRe+QEHdZoZC3h/RJQgG84Fy+ratG47qth/5bKLfDPO8i/oNxDRyHup5UXEwL4GVVwi2UxyUVhjooA5wo0/QpCFJA9RelqgtBWfJSLciGTQaFeJZ9rJlwpH6iOGGm6gum8GNsJ18YwR5gMRrf3ZBIv8I6in5kS3at5l92oqUaRbbw/TdTQeLuR+mFro/JBeBoJWqzFOE6Ywc3W+JziC9w7fU/Q4eqOnVeLJlEKTRM9zhQ+XNZLK7VJ8+Ur5BQVIVGhdYjz6R+oCiVCXk2kZFzkGAtwqlFgY29lCr69fs29K5aF46GvQd3FK+MpBCGaLXRmwahdpp640yLXMGjkLDJChdGZaJEYi02SZpjWdIWWiW0sPSBkK4JnBj6bmlDdW8Ao25pnzbMGxO76KF6vqwOG+k69ORfmPS2GMqKW6w1M37qC4WMzNnyb0q9gHzwgWxbZILENNNaeo71+MxTFAtIJGWNsuABkorDPRnE4Vh0ySO906fRvwd3xnCTh1qhaYyzFqHN7EapUVlxCYtVVVCGpooFAZuq4qYh2GAiMW7ruRSBZounXpGkJGRuBUZybsh9CwCeVC6djIMu64mEi010+F8VldhXd10FoZtFARnCEbONmqUb6ff3DbybxL6t243sih0o8qJohAgGs9Kp6ET2FxgavkuSo169Whiy0DsN4Jwrit9q53WQWmFASBrFypq9bCTqZNN6odrPBfjxJStm1Sr18VBoWgHRiarahECM28IJrls1/kIKY2qrRg0LrmcMKxquS7FqnSGbiKREZsd866Wmip0xPJZ6AMmqEzfxWqRy+eGZSCd2W9xhUjlK98i7ANs/U62/RBA/HcS28bEMAFRWRE0e4PVpWr56vfQ3hCYQBaYhI4xPctjypS3mp/zsmKhnACDNfGbmctoNUorDIjI6uusCQJA7Wzih44wgqDjJDCdu7R+1DBu1B8pCwPdyiTbPgUvHB3IzOFF3iojUM9YafwyPYcsRBnCpLVFVJoYuKpVBq9qbE4j5yvXM6kwMjF8NZ7UZhb3hZwZpDkDUxmmetbq+jQmKyl4p7NIRXqN31xBkqMhtBPlFi6rKkHS/hgDnVJWpv7g2AdtdTMpaAyOKDJxEGOICo43DkPh3k6UVhgM7SKjrxSQGaPKCMTfgzX5wzXCMeWLjE3vJjILI3Ug6Q69U2kwdc46yzFFhiliUGE4ZqarD4P1LgIVJjpN/t6aYhnYDyYLBKC5fDmvML9BVWMVwoOGDN0tA737Q1y1ZctKUgjEfGPThHTqvpvX/iZBY6YtUq5oBWsiJZmk1a2Cq4uSATB+Jykv1b1osW5NY8KqsCW8e8FbROCFB2o2+tsnEEorDIZ0kcTIVdTq3OikNs1cZQQNMzSGEeiW25m0QnlIykxGHNTmoylYEQxyPNPyVElbUdrKZDqLME3e2coxhll/xo36CU3M0KRJ6yBO6Ju+p+kKRuk7uWrMDGU1kdC2LNOkQ005OLHBvCRGElVougT3ZIDlm4WrW21CHzKTtNVU3GDZpRFaQVli3gEkK7iutyzsS0vFfOVxY0klxzONK4OCwqwuBDAWJKUf6p/GK7mJBDlXuYmahCFdpNGIRCYTjgTRVPPeQY4X5qAJ6SEyHHHysDFxZDEz6nW5EwabnOr1sOeok7+2OQeR+dRFRiRAvmVKcbNYBlxIs+w7NsE4H8LyGvWQLkVIsf4bqG4iOY0cls5A8t+pQl8c7CbLwHlpKRhdXVE3UZ1lt50Jzj5/RSDqLNLlW0Jh4HyBTEw1RTemftmt5XsIzz16oquJrJaFkm/oxjWn6R8UrX3zmLZNYpuEhpnO8GKqwXrdQJ/b92gWSisMhnaRVXPzBmIYFmFiBKIUjx0gQjyttiMIILUzyUyWhU1OCi1KfiL9Jq2ukUYRgCYLyAbV4jBdf2jKVx1QojAJ3kUmHA15DUqT8/G0A/6Evp+jak1IAtngsx9UlA0TrH3AlEahRb9KJ+w7O/YPyPQLFqXYn//8lKOFOIp1GOnrGsVFpVNhpnEr5wC5v7vsM3BdXiyeFBqxVEV6pEl8MzM391W7MmfKKxAGoptIdIGpY7rFsqC8wqArsAyEZ7Iv2tzZxcEvuk82CIdPxU4gC4xRz8zNjL3GcmfpapyFo2dKkr/BlndkUAplOkxSqolUbVw3kCMwaFsmF5jozovSJjMMmbEbBrgw9EUXQU3hjCJjTOsqDMuUd4SLCoXL5GyNGTt9Zv/LuetCMi1JBxVhFhQzfGhXpH8E+K/fL5PpFoSrcU8H6xUftW8ZSTV8T1nom9OLz//7oRcbYdfVWaZ+7+WtjycqknFliWlCN5HbPGTfYNMvgpRQWmEw1BcGItQlo4GGoH5MEyN47MVe2d9rwfqdBwAAK7bs1S7xi4bFwQvhXahRiLuJxQ69cuu+GM0nqEs9PI5CoXdA6KAPLNzYqGdkct3AwOsM6SABk7A0zptwaAEt2RTeB+BZORHJF6mDSL86WFVITC4oR7UMlL4S99yGFzftkeaNQteQaFmY+6oodJ5csU1bFxUq89KF1d/vO3tCmK9lPkJNL57iGu6WN8/HmCza6OQ6R56r6UXa3n/ORIkul7axb4ITaZO/h+pWjcNgjTFsSGAZiJsYzd9GtfaajdIKgyFdJK0WAZRBXTcNSrOL4M8mhJfR2PiA2FHue2G9fAS1/0pMLvquvd9yBwmYpCiYVPeJ6t9kMeznN1jTa8UqzZt390sCRNXyJNpEX7hmCaz626jhCVrlbxdsCrX8mnmC78BAqDmpO8pty4PDvEJpoO5nkNKrZj2CNG7CYOOuvoalWKuHeQvFR2kz0CLHUZisIY1JEw5oCDB+7GHSc1Fomsr1XDthH9AKvQjdENLr23aNcB+C69LS48eMkJ6brC65D4f0ROea9OMrsqjBoRvUmTHU3wUouon29A1K8xxS27R4Rrm8wkDjmjGF1TY3zRmcdPQoiZmaIEr+s04IBYhtk5FpsooFzcu6Ll4ZoOGZN+HzAYFh7j8om6BiPU8YO7IRb1A5p1ndNBUM/sG6/vhilQbTRPe+g/LNUAFqrKzpFt59/w8rGuElm/ZIVpOYRvweYj339g824qn+c2l5sSEsul9sOOuEMcKkNTfaVJrcVxUCiZmH5/t//vxXSELPpv2KmrnoQhOTSPNGhtV3cZOnYh8IlYMwXZoFBTPnbZBoNGUxWA8nZnv39EtuS5PgFgW/TaCbLBBxUUMXxXsJGnQOCSeQA9z3glxPEeJNeK1AeYVB0PCiVqcwspARyB/BzAjq+PXz6wF4KxK27zvYeNcnaKlimtedFN5vWmdgi39dXvQIChOd4SqbrXsPSu4bk8YsuklYyHuwVlfqqhcubzj1mAbzsK2yYQaeX7vTT1/HV369EACwa/8A/up7T2hKMc8NfGPWEqzaug8qvJ3i+jQi87ph9nJBk5XbRhQGA8JKkusfXi5PoNb07SEJ4RrjZ39a7cWvsyTEtu0Nr0IUMSgwYObQxaBahCKkSUaB4R0mHIHdvWaHcQLYpHhElQihTMn1EUIVlJt3i3NnwCL/ms9avY4FPbsaz4NlrNE5B307MyvaMIdxTK6YG2YvbyhL3/jtUolmMW+xD4jPr71/caOciBIjhM1uQ2Dd9tCKOSj0L9XDIE4g6+Zg1PJdN/vlhfIKA51loHTC4NVOdSWGKAAkBgF8T9BGX3ftQ43wJmFyWVzVIJa/p28Ad8/1boKqM2N3Xxiv279gW6VTdMXcMHt5g2bb6p9BwQX2+0WbQy0/ohWGYZP2rGqLotDb0xe2mxht1oKN2HdQduE0yjT4SHv39uOlXk8YvPvV44WJ3bqi4YUF/f2bT5FoC/Je2btPYQT6MGDW+FSXXCPMjLu6exq/p3z94UZYvBx97prtUvpeXwmo1cO8GYx5PTv9MPCnleF8gFhPsQ/YNFlR0Imb6GT65TSmDZYPLtyktS4B4Lan1zTCcl8HnvbrsHDDLnzxV/MBeO48cYyZ3JO1OmvrZ3OjLt64W7pKNMDa7fvl/i22jdoH/Hba0zdoUbDC8IGD8pi4dPqfGuHP/uKFRvi5teGYHqxzo38P1hi3+m047ogRmPHM2kb5WwWFwtUNmRdaLgyI6EIiWkZEK4hoWrPKGdIVNckeXLipEV6+Za9sIgsNH0z+AjKDUNe8i/hF97pGWBQMg3XGnFUeY/jhoy/htHGjAXja/1d9TRpA4+5kQO7sL27ag5W9ocYcCArbztzvCVryss17BJdP3eIHlQdlEM1zLYXvfvFsWM++Abe2ufHxleEPQ5mfePvLG+FhwuCu1fVaFAAtEwhgchMdrMnuMXG4iUL4eWEgi3Qu27THWObP/hQyycBiAuQ+eMPs5dKcwT3PeZbm9n0H8Zv5octAdOPVONQqV/XuazCvYw8fIbWNuHGtXme86GvmogDf2ycrPmI5uw+EjH37voONtlmq1LlXuAz+iRXhRUZiH5g1f6OUJhAMKmQrvI4/vtgbiSNa8SqdF71qfEPx06ULsFbQ3gcsm1FFLBUWMjwnKGu3PLnKmGbRhl2NsFgXtZ4Bzj7xSCxYH6a5W1A0aq5X6eUEauXGBiIaAuBFAOcD6AHwLIDLmHmxKc2UKVO4u7s7cVlX3TMfd8xZFx+xJJh83OHShqKiYvTwIZLVkBTjx46M3C9bobPwypcdEREwFaK48h2n4QsXvDJ1eiKay8xTXOO32jI4F8AKZl7JzAcB3Ang4mYUdCgJAgAdIQgAZBIEQPSi8Qqdh0oQuOEHj7zU0vJaLQwmABC5dI//rEKFChUqKBBdfM1Gq4WBzrkX8VMR0VQi6iai7t7eqA+xQoUKFQ4FDLXMi+WNVguDHgAnCr8nAtigRmLm6cw8hZmnjBtnvhC9QoUKFcqM4AiLVqDVwuBZAJOJ6BQiGg7gUgAzm1HQea88rhnZVqhQoUIp0VJhwMyDAD4J4HcAlgC4i5kXNaOsm64IJ9GPO2IErnzHaY3f93/q/zTCMz/5Zun5x98WxvvomyY1wo994R2N8Hc++JqwnI+E5Xzi7afhxKPDLf3f/Os/a4TFc1OmvStcIfDKlx3RCI8aPgRnjB/T+D3l5KO0aX50+esa4Q9OCfNV4/3mk2E9v37JqxrhD5wjp/nxh88B4G1ue/yLYT1v/6c/1+b1mhPDjXSHDRsiCd6rhPJv/YdzG+ExI4c2jvM4YexInCnU898vPqsRnv535zTCn3xHuNz0VH9JLuAdtnbxa09o/P7vD4TfQ8xLxWff+YpG+H/+Nkxz/d++thEee9gwXHauZ7x+9T1n4tdXhv1jzpfPa4TF8sVjED53flgGAHzqL8I6iN/gH4Q9El++SF4xcsLYkY3wV959RiMs9scbLjsbh48YCh0uEWj7z78J+6CY1zvPCL/ZBWcdj+9/6OzG7599LPxud/3fNzbC1//tazH2sGEAvL564Vkva7y77n1h3d5wangy6k8/+vpG+OqLwvLfMvlY6SykL1xwujb810Kc8UK7HDN6OL54YRhPHIfXCn3g/DOPb4Tfcfo4XPr60DEh9vUffCgcU2KZrzj+cPz168LfYhuK31Ms8/+9PeQhZ590JC76M6+dJh51WCMMAP/23jDNLR8N6X/Xq16G08aNxqeFvtMKtHRpaRqkXVpaoUKFCocyir60tEKFChUqFBCVMKhQoUKFCpUwqFChQoUKlTCoUKFChQqohEGFChUqVEAlDCpUqFChAiphUKFChQoVUAmDChUqVKiADth0RkS9ANbERtTjWABbY2N1DspUnzLVBShXfcpUF6Bc9UlSl5OZ2flwt8ILgywgou4kO/CKjjLVp0x1AcpVnzLVBShXfZpZl8pNVKFChQoVKmFQoUKFChXKLwymt5uAnFGm+pSpLkC56lOmugDlqk/T6lLqOYMKFSpUqOCGslsGFSpUqFDBAZUwqFChQoUK5RQGRHQhES0johVENK3d9JhARCcS0SNEtISIFhHRP/vPjyaih4houf/3KCHNVX69lhHRBcLzc4hogf/uBiJq3U3acp2GENHzRHR/CepyJBHdTURL/W/0xk6tDxF91u9jC4noDiIa2Ul1IaJbiGgLES0UnuVGPxGNIKJf+M+fIaJJbajPt/2+Np+I7iWiI4V3za8PM5fqH4AhAF4CcCqA4QDmATiz3XQZaB0P4HV++AgALwI4E8C3AEzzn08D8J9++Ey/PiMAnOLXc4j/bg6ANwIgAA8AeFeb6vQ5ALcDuN//3cl1uRXAP/rh4QCO7MT6AJgAYBWAw/zfdwH4aCfVBcBbAbwOwELhWW70A/h/AH7shy8F8Is21OcvAQz1w//Z6vq0fIC1oNO8EcDvhN9XAbiq3XQ50n4fgPMBLAMw3n82HsAyXV3g3SX9Rj/OUuH5ZQB+0gb6JwKYDeAvEAqDTq3LGHgMlJTnHVcfeMJgHYCjAQwFcL/PeDqqLgAmKcwzN/qDOH54KLxdvtSsuujqo7x7H4AZraxPGd1EQccP0OM/KzR8M+5sAM8AOJ6ZNwKA/ze4vdxUtwl+WH3ealwP4IsA6sKzTq3LqQB6AfzUd3vdRESj0YH1Yeb1AP4LwFoAGwHsYubfowProiBP+htpmHkQwC4AxzSN8nj8AzxNH2hRfcooDHQ+zEKvnyWiwwH8CsBnmHm3LarmGVuetwxE9B4AW5h5rmsSzbNC1MXHUHhm/I+Y+WwA++C5IkwobH18X/rF8FwMJwAYTUQftiXRPCtEXRyRhv7C1I2IrgYwCGBG8EgTLff6lFEY9AA4Ufg9EcCGNtESCyIaBk8QzGDme/zHm4lovP9+PIAt/nNT3Xr8sPq8lXgzgPcS0WoAdwL4CyL6OTqzLvDp6GHmZ/zfd8MTDp1Yn3cCWMXMvcw8AOAeAG9CZ9ZFRJ70N9IQ0VAAYwFsbxrlBhDRFQDeA+By9n08aFF9yigMngUwmYhOIaLh8CZPZraZJi38mf+bASxh5u8Ir2YCuMIPXwFvLiF4fqm/UuAUAJMBzPFN5D1E9AY/z48IaVoCZr6KmScy8yR4bf4HZv5wJ9YFAJh5E4B1RHS6/+g8AIvRmfVZC+ANRDTKp+E8AEvQmXURkSf9Yl7vh9d/W23BXQjgSwDey8yrf2NGAAAA20lEQVT7hVetqU+rJn9a+Q/ARfBW5rwE4Op202Oh8//AM93mA3jB/3cRPN/ebADL/b9HC2mu9uu1DMJKDgBTACz0330fTZ78iqnX2xFOIHdsXQC8FkC3/31+DeCoTq0PgH8DsNSn42fwVqZ0TF0A3AFvvmMAntb7sTzpBzASwC8BrIC3QufUNtRnBTw/f8ALftzK+lTHUVSoUKFChVK6iSpUqFChQkJUwqBChQoVKlTCoEKFChUqVMKgQoUKFSqgEgYVKlSoUAGVMKhQoUKFCqiEQYUKFSpUAPD/AQ4HD09qjkPvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(photodiode_data[120000:131800])\n",
    "plt.scatter(frame_timepoints[frame_timepoints>120000]-120000, frame_signals[frame_timepoints>120000]*200+200, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the frame detected, we can create our record master, often named reM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "265\n"
     ]
    }
   ],
   "source": [
    "ref_timepoints, ref_signals = extend_sync_timepoints(frame_timepoints, frame_signals, up_bound=len(photodiode_data))\n",
    "reM = RecordMaster([(ref_timepoints, ref_signals)])\n",
    "print(len(reM[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though the reM we just created is from a tiny portion of real data. From now one we will use a premade reM generated from the same dataset, in full."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Importing the record master\n",
      "\n"
     ]
    }
   ],
   "source": [
    "reM = import_record(\"./files/basic_synchro/reM_basic_synchro.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def parse_time(time_str, pattern=\"%y%m%d_%H%M%S\"):\n",
    "    \"\"\"Default parser of rhd timestamps. (serve as a template too)\"\"\"\n",
    "    return datetime.datetime.strptime(time_str, pattern)\n",
    "\n",
    "def get_position_estimate(stim_time, record_time, sampling_rate):\n",
    "    \"\"\"Estimate where in the record should a stimulus start, in sample points\"\"\"\n",
    "    if stim_time < record_time:\n",
    "        return -1\n",
    "    else:\n",
    "        return (stim_time - record_time).seconds * sampling_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-03-31 17:08:49\n"
     ]
    }
   ],
   "source": [
    "record_time = parse_time(\"200331_170849\") #Starting time of that example record found on the filename of the record\n",
    "print(record_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def match_starting_position(frame_timepoints, frame_signals, stim_signals, estimate_start, search_size=1000):\n",
    "    \"\"\"\n",
    "    Search the best matching index between frame_signals and stim_signals.\n",
    "    params:\n",
    "        - frame_timepoints: Indexes of the frames in the record\n",
    "        - frame_signals: Signals of the detected frames\n",
    "        - stim_signals: Expected stimulus signals\n",
    "        - estimate_start: Estimated start index\n",
    "        - search_size: Stimulus is searched in frame_signals[idx_estimate-search_size: idx_estimate+search_size]\n",
    "    return:\n",
    "        - best match for the starting position of the stimulus\n",
    "    \"\"\"\n",
    "    stim_matching_len = min(600, np.where(np.diff(stim_signals)!=0)[0][50]) #Way of getting the 50th change in the signals\n",
    "    idx_estimate = np.argmax(frame_timepoints>estimate_start)\n",
    "    search_slice = slice(max(0, idx_estimate-search_size-frame_signals.idx), min(idx_estimate+search_size-frame_signals.idx, len(frame_signals)))\n",
    "    return search_slice.start + np.argmax(np.correlate(frame_signals[search_slice], \n",
    "                                                       stim_signals[:stim_matching_len]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "match_starting_position seaks in the record the first frame of a stimulus. We can use functions from theonerig.synchro.extracting to find out the stimuli used in that record, and get their values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "flickering_bars_pr WARNING dt of frame #15864 was 50.315 m\n",
      "flickering_bars_pr WARNING dt of frame #19477 was 137.235 m\n"
     ]
    }
   ],
   "source": [
    "from theonerig.synchro.extracting import get_QDSpy_logs, unpack_stim_npy\n",
    "log = get_QDSpy_logs(\"./files/basic_synchro\")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "checkerboard             eed21bda540934a428e93897908d049e at 2020-03-31 17:09:26\n"
     ]
    }
   ],
   "source": [
    "print(log.stimuli[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(54039, 1, 18, 32) (54039,) None\n"
     ]
    }
   ],
   "source": [
    "#Unpacking the stimulus printed above\n",
    "unpacked_checkerboard = unpack_stim_npy(\"./files/basic_synchro/stimulus_data\", \"eed21bda540934a428e93897908d049e\")\n",
    "print(unpacked_checkerboard[0].shape, unpacked_checkerboard[1].shape, unpacked_checkerboard[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "get_position_estimate can approximately tell us where the stimulus should be to reduce the search time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimate position in sample points 1110000\n"
     ]
    }
   ],
   "source": [
    "estimate_start = get_position_estimate(log.stimuli[2].start_time, record_time, sampling_rate=30000)\n",
    "print(\"Estimate position in sample points\", estimate_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2235\n"
     ]
    }
   ],
   "source": [
    "stim_start_frame = match_starting_position(reM[\"main_tp\"][0], reM[\"signals\"][0], stim_signals=unpacked_checkerboard[1], estimate_start=estimate_start)\n",
    "print(stim_start_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def display_match(match_position, reference=None, recorded=None, corrected=None, len_line=50):\n",
    "    start, mid, end = 0, len(reference)//2, len(reference)-len_line\n",
    "    for line in [start, mid, end]:\n",
    "        if reference is not None:\n",
    "            print(\"REF [\"+str(line)+\"] \",\" \".join(map(str,map(int, reference[line:line+len_line]))))\n",
    "        if recorded is not None:\n",
    "            print(\"REC [\"+str(line)+\"] \",\" \".join(map(str,map(int, recorded[line+match_position:line+len_line+match_position]))))\n",
    "        if corrected is not None:\n",
    "            print(\"COR [\"+str(line)+\"] \",\" \".join(map(str,map(int, corrected[line:line+len_line]))))\n",
    "        print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the match we obtain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "REF [0]  0 0 0 0 0 0 4 0 4 4 4 0 4 0 4 4 4 0 4 0 4 0 0 0 0 0 0 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1\n",
      "REC [0]  0 0 0 0 0 0 4 0 4 4 4 0 4 0 4 4 4 0 4 0 4 0 0 0 0 0 0 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1\n",
      "\n",
      "REF [27019]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0\n",
      "REC [27019]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0\n",
      "\n",
      "REF [53989]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 4 4 4 4 4 0 0 0 0 0 0\n",
      "REC [53989]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 4 4 4 4 4 0 0 0 0 0 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "display_match(stim_start_frame, reference=unpacked_checkerboard[1], recorded=reM[\"signals\"][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a match!! But be sure to check it everytime, as mismatches occurs. Set then stim_start_frame manually"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def frame_error_correction(signals, unpacked, algo=\"nw\", **kwargs):\n",
    "    \"\"\"Correcting the display stimulus frame values. Shifts are first detected with one of\n",
    "    `shift_detection_conv` or `shift_detection_NW` and applied to the stimulus template. Then single frame\n",
    "    mismatch are detected and corrected.\n",
    "        - signals: true signal values recorded\n",
    "        - unpacked: stimulus tuple (inten,marker,shader)\n",
    "        - algo: algorithm for shift detection among [nw, conv]\n",
    "        - **kwargs: extra parameter for shift detection functions\n",
    "        \n",
    "        returns: stim_tuple_corrected, shift_log, (error_frames_idx, replacement_idx)\"\"\"\n",
    "    \n",
    "    if algo==\"no_shift\":\n",
    "        intensity, marker, shader  = unpacked[0].copy(), unpacked[1].copy(), unpacked[2]\n",
    "        if shader is not None:\n",
    "            shader = shader.copy()\n",
    "        error_frames, replacements = error_frame_matches(signals, marker, range_=5)\n",
    "        shift_log = []\n",
    "    else:\n",
    "        if algo==\"nw\":\n",
    "            shift_log = shift_detection_NW(signals.astype(int), unpacked[1].astype(int), **kwargs)\n",
    "        elif algo==\"conv\":\n",
    "            shift_log = shift_detection_conv(signals.astype(int), unpacked[1].astype(int), range_=5, **kwargs)\n",
    "        intensity, marker, shader = apply_shifts(unpacked, shift_log)\n",
    "        error_frames, replacements = error_frame_matches(signals, marker, range_=5)\n",
    "    if len(error_frames)>0:\n",
    "        intensity[error_frames]    = intensity[replacements]\n",
    "        marker[error_frames]       = marker[replacements]\n",
    "        if shader is not None:\n",
    "            shader[error_frames] = shader[replacements]\n",
    "    return (intensity, marker, shader), shift_log, list(zip(map(int,error_frames), map(int,replacements)))\n",
    "\n",
    "def error_frame_matches(signals, marker, range_):\n",
    "    \"\"\"Find the frames mismatching and finds in the record the closest frame with an identical signal value\"\"\"\n",
    "    error_frames = np.nonzero(signals!=marker)[0]\n",
    "    if len(error_frames)>0:\n",
    "        where_equal = [((np.where(marker[err_id-range_:err_id+(range_+1)] == signals[err_id])[0]) - range_) for err_id in error_frames]\n",
    "    #filtering out the frames where no match was found\n",
    "        tmp    = np.array([[wheq,err] for (wheq, err) in zip(where_equal, error_frames) if len(wheq)>0], dtype=\"object\")\n",
    "        if len(tmp)==0:\n",
    "            replacements = np.empty(shape=(0,), dtype=int)\n",
    "            error_frames = np.empty(shape=(0,), dtype=int)\n",
    "        else:\n",
    "            where_equal  = tmp[:,0]\n",
    "            error_frames = tmp[:,1]\n",
    "\n",
    "        #Choosing among the equal frame signals the one that is the closest\n",
    "            closest_equal = [wheq[(np.abs(wheq)).argmin()] for wheq in where_equal]\n",
    "            error_frames  = np.array(error_frames, dtype=int)\n",
    "            replacements  = error_frames + np.array(closest_equal, dtype=int)\n",
    "    else:\n",
    "        replacements = np.empty(shape=(0,), dtype=int)\n",
    "        error_frames = np.empty(shape=(0,), dtype=int)\n",
    "    \n",
    "    return error_frames, replacements\n",
    "\n",
    "def apply_shifts(unpacked, op_log):\n",
    "    \"\"\"Applies the shifts found by either shift_detection functions\"\"\"\n",
    "    inten, marker, shader = unpacked[0].copy(), unpacked[1].copy(), unpacked[2]\n",
    "    if shader is not None:\n",
    "        shader = shader.copy()\n",
    "    orig_len = len(marker)\n",
    "    for idx, op in op_log:\n",
    "        if op==\"ins\": #We insert a frame\n",
    "            marker = np.insert(marker, idx, marker[idx], axis=0)\n",
    "            inten  = np.insert(inten , idx, inten[idx], axis=0)\n",
    "            if shader is not None:\n",
    "                shader = np.insert(shader, idx, shader[idx], axis=0)\n",
    "        elif op==\"del\": #We concatenate without the deleted frame\n",
    "            marker = np.concatenate((marker[:idx],marker[idx+1:]))\n",
    "            inten  = np.concatenate((inten[:idx],inten[idx+1:]))\n",
    "            if shader is not None:\n",
    "                shader = np.concatenate((shader[:idx],shader[idx+1:]))\n",
    "    marker = marker[:orig_len]\n",
    "    inten  = inten[:orig_len]\n",
    "    if shader is not None:\n",
    "        shader = shader[:orig_len]\n",
    "    return (inten, marker, shader)\n",
    "\n",
    "def shift_detection_conv(signals, marker, range_):\n",
    "    \"\"\"Detect shifts with a convolution method. First look at how far the next closest frame are, and average\n",
    "    it over the record. When the average cross the -1 or 1 threshold, shift the reference accordingly.\"\"\"\n",
    "    marker = marker.copy()\n",
    "    shift_detected = True\n",
    "    shift_log = []\n",
    "    while shift_detected:\n",
    "        error_frames, replacements = error_frame_matches(signals, marker, range_)\n",
    "\n",
    "        all_shifts = np.zeros(len(marker))\n",
    "        all_shifts[error_frames] = replacements-error_frames\n",
    "        all_shifts_conv = np.convolve(all_shifts, [1/20]*20, mode=\"same\") #Averaging the shifts to find consistant shifts\n",
    "\n",
    "        shift_detected = np.any(np.abs(all_shifts_conv)>.5)\n",
    "        if shift_detected: #iF the -.5 threshold is crossed, we insert a \"fake\" frame in the reference and we repeat the operation\n",
    "            change_idx = np.argmax(np.abs(all_shifts_conv)>.5)\n",
    "            if all_shifts_conv[change_idx]>.5:#Need to delete frame in reference\n",
    "                #Need to refine index to make sure we delete a useless frame\n",
    "                start,stop = max(0,change_idx-2), min(len(marker),change_idx+2)\n",
    "                for i in range(start,stop):\n",
    "                    if marker[i] not in signals[start:stop]:\n",
    "                        change_idx = i\n",
    "                        break\n",
    "                shift_log.append([int(change_idx), \"del\"])\n",
    "                marker = np.concatenate((marker[:change_idx], marker[change_idx+1:], [0]))\n",
    "            else:#Need to insert frame in reference\n",
    "                shift_log.append([int(change_idx), \"ins\"])\n",
    "                #inserting a frame and excluding the last frame to keep the references the same length\n",
    "                marker     = np.insert(marker, change_idx, marker[change_idx], axis=0)[:-1] \n",
    "    return shift_log\n",
    "\n",
    "def shift_detection_NW(signals, marker, simmat_basis=[1,-1,-3,-3,-1], insdel=-10, rowside=20):\n",
    "    \"\"\"Memory optimized Needleman-Wunsch algorithm.\n",
    "    Instead of an N*N matrix, it uses a N*(side*2+1) matrix. Indexing goes slightly differently but\n",
    "    result is the same, with far less memory consumption and exection speed scaling better with\n",
    "    size of the sequences to align.\"\"\"\n",
    "    #Setting the similarity matrix\n",
    "    side = rowside\n",
    "    sim_mat = np.empty((len(marker), side*2+1), dtype=\"int32\")\n",
    "    #Setting the errors\n",
    "    insertion_v = insdel #insertions are commons not so high penalty\n",
    "    deletion_v  = insdel #deletions detection happens during periods of confusion but are temporary. High value\n",
    "    error_match = np.array(simmat_basis) #The value for a 0 matching with [0,1,2,3,4]\n",
    "    error_mat = np.empty((len(simmat_basis),len(simmat_basis)))\n",
    "    for i in range(len(simmat_basis)):\n",
    "        error_mat[i] = np.roll(error_match,i)\n",
    "                \n",
    "    #Filling the similarity matrix\n",
    "    sim_mat[0, side] = error_mat[marker[0], signals[0]]\n",
    "    #Initialization: Setting the score of the first few row and first few column cells\n",
    "    for j in range(side+1, side*2+1):\n",
    "        sim_mat[0,j] = sim_mat[0,side] + insertion_v*j\n",
    "    for i in range(1, side+1):\n",
    "        sim_mat[i,side-i] = sim_mat[0,side] + deletion_v*i\n",
    "          \n",
    "    #Corpus: if j is the first cell of the row, the insert score is set super low\n",
    "    #        if j is the last  cell of the row, the delete score is set super low\n",
    "    for i in range(1, sim_mat.shape[0]):\n",
    "        start = max(side-i+1, 0)\n",
    "        stop  = min(side*2+1, side+sim_mat.shape[0]-i)\n",
    "        for j in range(start, stop):\n",
    "            if j==0:#j==start and i>side:\n",
    "                insert = -99999\n",
    "                delete = sim_mat[i-1, j+1] + deletion_v\n",
    "            elif j==side*2:\n",
    "                delete = -99999\n",
    "                insert = sim_mat[i, j-1] + insertion_v\n",
    "            else:\n",
    "                insert = sim_mat[i, j-1] + insertion_v\n",
    "                delete = sim_mat[i-1, j+1] + deletion_v\n",
    "            match  = sim_mat[i-1, j] + error_mat[marker[i], signals[j+i-side]]\n",
    "            sim_mat[i,j] = max(insert,delete,match)\n",
    "            \n",
    "    #Reading the similarity matrix\n",
    "    #In general, it's the same, at the difference that when i decrement, must add 1 to j compared to usual.\n",
    "    i = len(marker)-1\n",
    "    j = side\n",
    "    shift_log = []\n",
    "    while (i > 0 or j>side-i):\n",
    "        if (i > 0 and j>side-i and sim_mat[i,j]==(sim_mat[i-1,j]+error_mat[marker[i], signals[j+i-side]])):\n",
    "            i -= 1\n",
    "        elif(i > 0 and sim_mat[i,j] == sim_mat[i-1,j+1] + deletion_v):\n",
    "            shift_log.insert(0,(j+i-side+1, \"del\")) #Insert the j value for deletion too because all shifts\n",
    "            i-=1                                        #are relative to the signals recorded, unlike normal NW\n",
    "            j+=1\n",
    "        else:\n",
    "            shift_log.insert(0,(j+i-side, \"ins\"))\n",
    "            j-=1\n",
    "                \n",
    "    return shift_log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We correct the stimulus values with frame_error_correction and it gives us back the changes it made to keep track of the errors made."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(53516, 'ins'), (53537, 'del')] 45\n"
     ]
    }
   ],
   "source": [
    "signals = reM[\"signals\"][0][stim_start_frame:stim_start_frame+len(unpacked_checkerboard[0])]\n",
    "corrected_checkerboard, shift_log, error_frames = frame_error_correction(signals, unpacked_checkerboard, algo=\"nw\")\n",
    "print(shift_log, len(error_frames))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def chop_stim_edges(first_frame, last_frame, stim_tuple, shift_log, frame_replacement):\n",
    "    \"\"\"Cut out the stimulus parts not containing actual stimulus, and change the idx values of `shift_log`\n",
    "    and `frame_replacement` to match the new indexing.\"\"\"\n",
    "    inten, marker, shader = stim_tuple\n",
    "    if last_frame<0: #Using negative indexing\n",
    "        last_frame = len(marker)+last_frame\n",
    "    inten = inten[first_frame:last_frame]\n",
    "    marker = marker[first_frame:last_frame]\n",
    "    if shader is not None:\n",
    "        shader = shader[first_frame:last_frame]\n",
    "    \n",
    "    shift_log = [(shift[0]-first_frame, shift[1]) for shift in shift_log if shift[0]<last_frame]\n",
    "    frame_replacement = [(fr[0]-first_frame, fr[1]-first_frame) for fr in frame_replacement if fr[0]<last_frame]\n",
    "    \n",
    "    return (inten, marker, shader), shift_log, frame_replacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def detect_calcium_frames(scanning_data, epoch_threshold=-8):\n",
    "    \"\"\"Detect the timing of the 2P frames, epoch by epoch over a record.\"\"\"\n",
    "    \n",
    "    #Finds the start of a stack recording\n",
    "    start_set = np.where((scanning_data[1:] > epoch_threshold) & (scanning_data[:-1] < epoch_threshold))[0]\n",
    "    #Finds the end of a stack recording\n",
    "    end_set   = np.where((scanning_data[1:] < epoch_threshold) & (scanning_data[:-1] > epoch_threshold))[0]\n",
    "    #Splits the records into the epochs\n",
    "    list_epoch = np.array_split(scanning_data, np.ravel(list(zip(start_set, end_set))))[1::2]\n",
    "\n",
    "    def detect_peak_sync(epoch):\n",
    "        #Finds the peaks in an epoch. Peaks have strong SNR so this works fine\n",
    "        return signal.find_peaks(epoch, prominence=2)[0]\n",
    "\n",
    "    return [arr + start_set[i] for i, arr in enumerate(list(map(detect_peak_sync, list_epoch)))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 00_core.ipynb.\n",
      "Converted 01_utils.ipynb.\n",
      "Converted 02_processing.ipynb.\n",
      "Converted 03_modelling.ipynb.\n",
      "Converted 04_plotting.ipynb.\n",
      "Converted 05_database.ipynb.\n",
      "Converted 06_eyetrack.ipynb.\n",
      "Converted 10_synchro.io.ipynb.\n",
      "Converted 11_synchro.extracting.ipynb.\n",
      "Converted 12_synchro.processing.ipynb.\n",
      "Converted 13_leddome.ipynb.\n",
      "Converted 99_testdata.ipynb.\n",
      "Converted index.ipynb.\n"
     ]
    }
   ],
   "source": [
    "#hide\n",
    "from nbdev.export import *\n",
    "notebook2script()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tor",
   "language": "python",
   "name": "tor"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}