{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4001a375",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "sys.path.append(os.path.join(os.getcwd(), '..', '..'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "570920e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from session.sessions import selected_009266, selected_008229, selected_009265, selected_57\n",
    "from imports import *\n",
    "from loading import load_session_data\n",
    "from target import get_spike_counts, build_tgt_matrix, build_silence_matrix\n",
    "from scipy import stats\n",
    "from scipy import signal\n",
    "import pandas as pd\n",
    "from functools import reduce\n",
    "from statsmodels.formula.api import ols, glm\n",
    "from postprocessing.spiketrain import instantaneous_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "26ba4176",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62751662",
   "metadata": {},
   "source": [
    "## 01 - Syllable ratio matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "dd986c29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['57_SIT_2023-12-18_14-07-34',\n",
       " '57_SIT_2023-12-22_14-08-07',\n",
       " '57_SIT_2023-12-22_14-43-58',\n",
       " '57_SIT_2023-12-22_17-37-18',\n",
       " '57_SIT_2023-12-23_14-21-01',\n",
       " '57_SIT_2023-12-28_16-43-28',\n",
       " '57_SIT_2023-12-29_11-06-26',\n",
       " '57_SIT_2023-12-29_11-40-14',\n",
       " '57_SIT_2023-12-29_12-11-46',\n",
       " '57_SIT_2024-01-02_14-43-18',\n",
       " '57_SIT_2024-01-02_16-38-05',\n",
       " '57_SIT_2024-01-02_17-10-09',\n",
       " '57_SIT_2024-01-03_19-22-18',\n",
       " '57_SIT_2024-01-03_19-54-59',\n",
       " '57_SIT_2024-01-04_14-16-22',\n",
       " '57_SIT_2024-01-04_14-52-59',\n",
       " '57_SIT_2024-01-05_14-35-49',\n",
       " '57_SIT_2024-01-05_15-08-34',\n",
       " '57_SIT_2024-01-06_16-52-40',\n",
       " '57_SIT_2024-01-06_17-25-35',\n",
       " '57_SIT_2024-01-07_19-23-28',\n",
       " '57_SIT_2024-01-08_15-51-26',\n",
       " '57_SIT_2024-01-12_13-23-02',\n",
       " '57_SIT_2024-01-15_13-45-22',\n",
       " '57_SIT_2024-01-15_14-34-48']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#source = '/home/sobolev/nevermind_ag-grothe/AG_Pecka/data/processed/'\n",
    "source = '/home/sobolev/nevermind/AG_Pecka/data/processed/'\n",
    "\n",
    "sessions = selected_57\n",
    "sessions.sort()\n",
    "sessions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "63c25725",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Session 009265_hippoSIT_2023-02-27_10-18-32 done\n",
      "Session 009265_hippoSIT_2023-02-27_15-33-46 done\n",
      "Session 009265_hippoSIT_2023-02-28_09-16-50 done\n",
      "Session 009265_hippoSIT_2023-02-28_13-16-10 done\n",
      "Session 009265_hippoSIT_2023-02-28_20-45-04 done\n",
      "Session 009265_hippoSIT_2023-03-01_10-46-12 done\n",
      "Session 009265_hippoSIT_2023-03-02_09-32-54 done\n",
      "Session 009265_hippoSIT_2023-03-02_16-27-42 done\n",
      "Session 009265_hippoSIT_2023-03-02_20-11-35 done\n",
      "Session 009265_hippoSIT_2023-03-03_09-37-07 done\n",
      "Session 009265_hippoSIT_2023-03-03_16-00-47 done\n",
      "Session 009265_hippoSIT_2023-03-04_11-12-04 done\n",
      "Session 009265_hippoSIT_2023-03-05_11-52-17 done\n",
      "Session 009265_hippoSIT_2023-03-05_18-31-32 done\n",
      "Session 009265_hippoSIT_2023-03-08_18-10-07 done\n",
      "Session 009265_hippoSIT_2023-03-09_20-03-08 done\n",
      "Session 009265_hippoSIT_2023-03-10_09-57-34 done\n",
      "Session 009265_hippoSIT_2023-04-13_09-54-39 done\n",
      "Session 009265_hippoSIT_2023-04-20_11-39-02 done\n"
     ]
    }
   ],
   "source": [
    "win_l   = 2    # seconds\n",
    "step    = 1    # seconds\n",
    "s_rate  = 100  # Hz\n",
    "syl_num = 10\n",
    "\n",
    "for i, session in enumerate(selected_sessions):\n",
    "    session_data = load_session_data(session, load_aeps=False)\n",
    "    moseq_file = session_data['moseq_file']\n",
    "    moseq      = session_data['moseq']\n",
    "    #animal = session.split('_')[0]\n",
    "    #moseq_file = os.path.join(source, animal, session, 'moseq.h5')  # can be _full.h5 for dark sessions\n",
    "    \n",
    "    # build syllable ratio matrix\n",
    "    idxs_srm_tl = np.arange(0, len(session_data['tl']), step*s_rate)\n",
    "    syl_ratio_mx = np.zeros([len(idxs_srm_tl), syl_num])\n",
    "    for k, idx in enumerate(idxs_srm_tl):\n",
    "        curr_syls = moseq['syllables reindexed'][idx:idx + win_l*s_rate]\n",
    "        for j in np.arange(syl_num):\n",
    "            syl_ratio_mx[k, j] = np.sum(curr_syls == j) / (win_l*s_rate)\n",
    "    \n",
    "    # roll 1 sec to match\n",
    "    syl_ratio_mx = np.roll(syl_ratio_mx, step, axis=0)\n",
    "    syl_ratio_mx[0] = syl_ratio_mx[1]\n",
    "    \n",
    "    with h5py.File(moseq_file, 'a') as f:\n",
    "        if 'syl_ratio_mx' in f:\n",
    "            del f['syl_ratio_mx']\n",
    "        if 'idxs_srm_tl' in f:\n",
    "            del f['idxs_srm_tl']\n",
    "        \n",
    "        id_h5 = f.create_dataset('idxs_srm_tl', data=idxs_srm_tl)\n",
    "        ds_h5 = f.create_dataset('syl_ratio_mx', data=syl_ratio_mx)\n",
    "        ds_h5.attrs['win_l']    = win_l\n",
    "        ds_h5.attrs['step']     = step\n",
    "        ds_h5.attrs['syl_num']  = syl_num\n",
    "       \n",
    "    print(\"Session %s done\" % session)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ea4c1611",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA50AAAD8CAYAAAD0bMAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABIRklEQVR4nO3deXxU1f3/8ddhIiEQlGXAyBckgqiJEaMFFEQkqCCbsooFFbTiV4Vv+QlthdaFfoWKbbUWW75oCxGqtCIo0WBZ3YJBkEAgIFVEAoQQthBEFkG5vz9mmMxMJslkMpOZSd7Px2MeucuZez53Mttnzjn3GMuyEBEREREREQmFeuEOQERERERERGovJZ0iIiIiIiISMko6RUREREREJGSUdIqIiIiIiEjIKOkUERERERGRkFHSKSIiIiIiIiGjpFNERGqcMWaqMWZrDdTzmjEmM9T11CbGmDHGmO/8KJdvjPlFTcQkIiLRTUmniIj4zZnEWW63w8aYTGPMVeGOrRwTgHvDHUSUeRNod36lgh8IOgOzaiwqERGJWko6RUSkqlYBlzhvvYE44J2wRlQOy7KOWZZVUp1jGGMuCFI4UcGyrFOWZR30o9why7JO1kRMIiIS3ZR0iohIVX1vWVaR87YR+BNwlTEm7nwBY8x/GWP+ZYw56rwtNcZ08D6QMeYeY8xOY8xxY8wSY4zdbV9nY8wKZ2vqt8aYNcaYrm77FxhjFnsdr54xZq8xZqJz3aN7rTEm1hjzkjHmgDHmtDHmM2NMd7f9PZ0tuP2MMeuNMWeAPr4eBGPMRcaY/zPG7Hcea7sxZoTb/iHGmDxjzPfOmH5jjDFu+/ONMU87YzzuLDPCGNPE+dh9Z4zZYYzp7SO+AcaYXGe9OcaYn3jFVlndQ4wxW4wxp4wxxcaYj40xFzv3ubrXGmPGAM8AV7u1bo9xi/8Xbse81BjzjvNcjhtj3jbGtHbbP9UYs7Wi/7mIiNROSjpFRCRgxpjGwAggz7KsU85tDYEPgdPALUBXYD+wyrnvvETnfQfjaDG9Dpjutr8x8A/gZqALkAu8b4xp7tz/OtDfGHOR231uwdEC+89yQv69s84HnfXlAcuMMZd4lXseeBK4Cljn47wN8L6zvgeAZGAicMa5/yfAW8DbwDXAZGAKMN7rUP8PWA9cDywE5gELnMdOBT4BXjfGNPC63x+BJ4BOwDdA5vnHtrK6jTEJwL+cdSUBPXA8zr68CbwAfElp6/abPh6PekAGcDGQ5ry1Apa4J7tU/j8XEZHayLIs3XTTTTfddPPrBrwG/AB857xZwB4gxa3Mg8AOwLhtswFHgLud61NxJKUXuZX5DfB1BXUbHMnrvc71GOAA8DO3Mn8HVnjFm+lcboQjKbzfK66dwDTnek/nOQ2t5HG4HTgHJJWz/w3gA69tU4ECt/V84J9u6/HOume6bUt0buvkFd8or/uVAA/5UzeOBNcC2pYT+xjgO6/7bvVRLh/4hdvj8SOQ6La/nfMxui3Q/7luuummm26146aWThERqapPcLTCpeJogVwNrDDGtHHu/wlwGXDc2UX0O+AY0BRo73ac3ZZlHXNbLwRanl8xxrQ0xrxijPnKGHMMOO7cfymAZVk/4Gh1G+UsHwsMxdEC6kt74ALg0/MbLMv6EViLo6XS3YZKHoPrgP2WZW0vZ3+Sez1Oa4D/MsZc6LZti1ss3wEncbS+nnfA+bclntZ63S+P0nOorO7NOMblbjXGLDbGPGqMaVHOefgrCSi0LCvfLa5vcPxP3R/bCv/nIiJSO8WEOwAREYk6Jy3L+vr8ijHmIRxJ5cPAUziGbuQC9/i4b7Hb8lmvfRaewz7m4eiu+TiOVrXvcSS49d3KvA6sNcb8F3CDc9/bVT0hZ93uTgRwjEDq8vUYnPVRNlg/EluWZf3oHCd6I44urj8DnjPG3GJZ1uYg1eNRp9tyZf9zERGphfRGLyIi1WXh6EZ5frzmRuBy4LBlWV973YrLPUpZ3YGXLctaalnWNhwtnR5jLy3LWg98DfwUR4tnhrPlz5edOLrX3nR+gzHGhmPM6RdViAtgE3CJMSapnP3b3etx6o6ji+vxKtbly43nF4wxjYAUZ51+1W05rLUs67c4pj4pxDHW0pczOLohV2Q70MoYk+gWVzsc4zqr+tiKiEgto5ZOERGpqljnxWjA0WV2PI5xhe85t70B/ALIMMY8jWPMZxvgLmC2ZVk7/KznK+BeY8w6HOMxf4/zQj1e3gAewjH+cUh5B7Ms64Qx5v+A540xh4FdOFpRL6bq802uxnGBocXGmMedsV4ONLIsawmOi+98boyZiuPCQJ2BScCvq1hPeZ40xhzCkSw+jeNxWeDcV2HdxpgbgduA5Ti6716H4/9TXnKYD7Q1xlyP43953LKs773KrMLRVfgNY8wE57aXcfwA8UF1TlRERKKfWjpFRKSqbsNxQZ/9OBKvzsBwy7I+ArAcczf2wHFV1beA/+DoKtsUOFqFeh7Ekczm4Lja6lwcCZC314ErcXTxXVHJMZ/AMQ40HUcX4I7AHZZl7a9CXFiWdQ7oi2Ps5Os4Wvr+jLPrr+WYSmY4jjGmW4EZzttfqlJPBSbjSC43Ah2AAZZlnfCz7mM4WkIzcVzw6QXgWcuyyhsLuxjH1XRXA4dwtCp7sCzLwvGjwiEcVy7+ECgCBjn3iYhIHWb0WSAiIhIdjDE9cSR0LSzLOhzeaERERPyjlk4REREREREJGSWdIiIiIiIiEjLqXisiIiIiIiIho5ZOERERERERCRklnSIiIiIiIhIyNTJPp91utxITEwE4s690zu76/xXvc5v3+oFvvvY43sXtLvfY5r1+flthYaFrvVWrVmXWgYDK+DqHUNVVUd3n6/eOx/0+5R0n0Pgq+18Fq4w/51nVx6K88zy4u3Se9pZtGwOU2Rasx8sf3nW7r5/fdnrrNtd6g5Sr2XZkm0eZq5tfXeFxyjtP7+P6OodA/+eV8ed56x7f+Ri9Y/Z1DoE8//15vPx5LHzFU1nMvs7Tn8fYn/+5d13+nGeoygT63h6s52Sw3v/94f24+/N+5n0f78cGqPDxOl/Gn9eIP88Lf+ry3has93Z//lf+8K7L/b3z6uaOx8J7WyDvi/78P/15PwvWZyyFm0oranUdUHPfW3x9hvkjWK9hf74XbD5+0rXt2sYNy6wDHD+e59rWuPE1AT23ff3P/XkvDdZ7nvfr01dd3s9/X6+Ryury5/O8qu8F1Xr++yGQ91dfn08VvZeC7/9xoK+9ij5Tz59DsB7TMs8VH+8pADk5OYcty2qBDzUyprNTp07Whg0bAHhhxADX9klvZgJQMDnLta31jJvLrLvf5/z9vI/jq8zUqVNd61OnTi2z7qtuX2W8t3nfx1cZf44TSBn3us/XX9E5lHecQOOr7H8VrDL+nGdVH4vyzvOvj5TOWz5udi+AMtuC9Xj5w7tu9/Xz27ZfleRaT/rPdq6Zd41HmbzReRUep7zz9D6ur3MI9H9eGX+et+7xnY/RO2Zf5xDI89+fx8ufx8JXPJXF7Os8/XmM/fmfe9flz3mGqkyg7+3Bek4G8rwI5LkNZR93f97PvO/j6/OzosfrfBl/XiP+PC/8qct7W7De2/35X/nDuy7398680Y6kwntbIO+L/vw//Xk/C9ZnLFMvKq1o6jHHnxr63uLrM8wfwXoN+/O9IOHDXNe2orTUMusAqz9o79p2a6+dAT23ff3P/XkvDdZ7nvfr01dd3s9/X6+Ryury5/O8qu8F1Xr++yGQ91dfn08VvZeC7/9xoK+9ij5T/c2x/C1T5rni4z0FwBiTY1lWJ3yokZZOEREREZFwe+j0reEOQaROUtJZRzXeviHcIYhEvYySs67lcRVsExEREfGlRVGPcIdQI5R0Sq1yfMnDpSszttdo3fr1VCS61JUPehERkXBT0ikiEWvEv94sXanCuC0REZFQenPX867lSfg/vlvC5+zZsxQUFHD69Oky+7qP+6Vrefv27T63VbR+ftvZv/7FY73zvU3LlHHf5msdoE+fPh7bvNeBMvfzdQ4/3Fl6IaXj27eXWfe3jHddDW58jtYbn+eCMyX4KyKTzkh/Iav7nIjUJb0+cn+nq9keBMESyl4QFV387vzFHKR2CWevGhEJTEFBAY0bNyYxMRFjjMe+ovqlKVFC+w4+t3mvb2l8iccxklo34dSPP7rW45KSOLj7W48yLdte6LHN1zr4e/Vaz/v5OoczBaVXna3funGZdX/LuNfV4tLGHLEdp4AnuOyzKfgrIpJOJXH+c0/IITKTchEJvtYN3K+Sd6zMukQX76TFY925TcQfGltefY98vKR0JS2VN6yhbnt31nQ4EiKnT5/2mXBK1RhjaN4ohkMXtavS/eqFKB4REREREZGIoYQzOByPY9Uey4ho6RSR8KgN3Sb9oQvGiIiEh+fnDNTmzxqRykyfPp1/zHudejYb9Uw95qT/jcsSksotf8MNN7Bp0yYAOnTowI4dO/yqZ/ny5WRnZ/PosAcZ/fOHeO8fb5cpc+TIEYYNG8bn6z/nvuEj+fO0FwI7KT9FbdLZoOnEcIcgIiIiIiJRKHHy0nL2fOXHtrJl8mf0r7C+z3PWk5mZycrMT4iNjeVI8REuujgWzlZ4t4BkZWWRlpbGmvVrualzV59lGjRowLPPPktu1ga2fflF8IPwErVJZ7BE+kWLapL33J2aAiT4Vn/Q3rV8a6+d5W6LdnpdhcfsrhNcy+PIC2lLtvex60qreW19bt89pfTrQB5ln0uhpPGHIhJK5y/25n5V1qKd/rUYBtPBQ0XY7XZiY2MBaN6sOS1bXcjif77H39NfYd7fFgCwcuVKXnzxRebMmVPlOjKWLuXl2a+wp2AfGRkZHNhfROP4xqzftIHMle97lG3UqBHdu3fnP+u3Vv/k/KAxnUHSoqiH6yYiIiIiInJez5t7sXfvXrqmXc8TT04k+7M1AHTv2oOvd37F4SOHAUhPT2fEiBEB1XFX//6sfDeDlJQU8vLyuPqqq1m/bA1vp79Z+Z1DrM63dIrUVjXZUiEiIhIN1LIu4dKoUTw5OTm8+9YyPl2bxcPjH6Dk1PP0SxvCsCEjWPTOm/y83aOsXbuWGTNmBFzPzl27aNfOcWXZkydP0Di+cbBOoVqUdEqVVHiZ/yi9xL+vboF1patgTfKey3Dq1KmudfdliTzeXS9FIpH3e4z3OtTe7tEiUsr7dX5+vfsPv6zxWBof3+OxbrOlcFPXm7mp680kXZXM4sVv0XlAZ3qP7sX4UeNp0aoJw4cPJybGM0WLORsPGOff8vUZNITio0exgOTkZPYX7qdzn5v40//+gV6t+wT57KpGSWcYqSuuSCmNIRYREYl+vj7P9RkPX+3aRcPYWC6qfzEAW7/Io23btgC0TGhJi4QWTJs2jVWrVgVcx/IlbzPlmamMf3wieXl5fLuvmEdGjw1K/NWlpFNEREREJMhqQ28wCZ7vTp7k0dGjOXKoGFtMDJe1vYx5r6dzgAMA9B/Wn8XHFpOUlERhYWGFx0pNTWVFxic+9+Vt+4LU1FRmzZrF/4z87wqPk5iYyLclxzhz9izvLV/KitUrufzCNoGdYCWUdEqtFo2/tlV8VVDHNhGJHOqyGXre79uR/j4eDabyJ481kWBxHzcLkTt29rOxVwCQ0L6Dx9VsE9p3ACizzXt9S0FJleq7/uqryc7O5uDub13b7PYLOXDYkXRu+mwTY8c6WiWP1j8KwLJNy7Db7Rw88S27vihNRHNzcz2O4y5z0UJsNhtz5szhTMHxCmPKz8/3KFO/deNK7xMoJZ21QDR+4YnGmEVqivsYSojMcZQa51m7HN8e+EUrRKpq+1VJruWk/2z3GH8LpWNwpfZy/x4I+i549613E9cwjvRZ6eEOJWSUdIqIiIhIraSr1Uo0WLh6IYBrDs/aSPN0ioiIiIiISMiopVNEJMzUVTVyjfiX24TamtpHREQkIEo6xaXf5p3hDkFCaHbXCa7lcUptok7i6QWu5fzwhSFSJd5dG6PlAiPRblzCYLe1Y2GLoy7Rc1sAdjRp7VruGMY4IpG614qIiIiIiEjIqKVTRMJCUx5IKPVNesy1nKeWfRGJANFwZXIJrenTp/OPea9Tz2ajnqnHnPS/Ed8+vtzyva/vzeaNm4H6XJbcymPalIosX76c7OxsHh32IKN//hDv/ePtMmVWrlzJ5MmT+f7EaerXr89zv3mWPiMGBHpqlVLSKeXSpMbBd+WK10pXepW/rTaqbJocj+cbhP05p+5pdU+Loh7hDiGiNWg6MdwhVNtr/Xa7lieFMQ5/tW7g/gXQ8T6kq7FGr7xde8IdgrhJ+Een0mVf+ytZL9N9dmrF3xXW5eaSmZnJ5uWvERtbn8PFRznTrA1HOepnxP7LysoiLS2NNevXclPnrj7L2O123nvvPeznGrPtP18w4N7B7BvhX1IbiFqTdGo8Yim1IImIv6JxrK9+EBMRKPt9py5///H+YVfzoUee/YcPY7fbiY2tD4C9WVNo1Ypl7yzjjb+9wcz5MwFHC+SMl2Ywc97MKteRsXQpL89+hT0F+8jIyODA/iIaxzdm/aYNZK5836PsddddB8CZguMkX5nEqdOn+P777zHVPM/yaEyniIhU2fHtM1w3EZHzWhT1cN1EpNRt3bqxd+9erug+iMemPMfHa3MA6NK9C7t27KL4cDEA6enpDBk5JKA67urfn5XvZpCSkkJeXh5XX3U165et4e30Nyu83zvvZ5B6TWpI5wmtNS2dUru4/0IH+pVORKKHet6ISDB59+5QF+voFN+wITk5OWS9M5cPsz9nxKOTmXE8hs4DOjNg+AAyF2VyzfhrWLt2LU+8+ETA9ezctYt27doB8O13JZy94ATF358ggcY+y3/x5XZ+/bunWfrGkoDr9IeSzjqi8fYNYatbXTyim7oyikhNirTx3SIiwWKz2ejZrRM9u3Ximqs6MG/xYjoP6MzgkYMZP2o8ic0TGT58ODExgaVofQYNofjoUSwgOTmZwn37uG3gnUx7+ikGte9QpnxBQQHDx45k7kuv0j6xXTXPrmJKOkUk6NTSIyIiUjvp6uCB+WrXLhrGxtKhkWM9d9uXtG3bFoCWCS1pkdCCadOmsWrVKs5xLqA6li95mynPTGX84xPJy8ujaHc+Y+4d5bNsSUkJ/fv3Z/qU39Kt840B1VcVUZt09vrIvUNBaH8FrcsD00VEREREpHq+O3mSR0ePpuRwETExNi5PbMOr8xdygAMA9B/Wn8XHFpOUlMS2w9sqPFZqaiorMj7xuS9v2xekpqYya9Ys7h82tNxj/OUvf+Hrr79m+kvPM/0lR6/ElR+uoglxAZ5hxaI26QyVcHZD9SXS4hGR2i8ar2gbbZLuCd1l6SW43OdW1KtBpKxIG2M6KSnLr3JF9zm+Yye070DRzh2u7QnObqje27zXtxSUeByvzBQqXq6/+mqys7OhcFPpRrudA4cdSeemzzYxduxYj/us2LgCu93OwRPfeszRmZuby8Hd3/qsJ3PRQmw2G3PmzPGI2duTTz7Jk08+yZmC465t9Vs29lgPJiWdIiIiIiIiYXL3rXcT1zCO9Fnp4Q4lZJR0ikhY6AJTIiIidUOktYZGmoWrFwKEdMqScFPSKSIiEUsXpRIREYl+SjrrKH2RExGRqgjnWN9AekaoN4VI7de6wQC3tWNhi6M8O5q09livbNxnbaaksxbyfAFCJL4I6wp96RFfdGESERERqUuUdPpB/dBFJBrk7doT7hBEREREylDSKXXO8SUPl67MCO0crzVFyYaISM1RbwURCcT06dNZMH8uNls96pl6vDJ3PvHt48st3/v63mzeuBmoz2XJrTymTanI8uXLyc7O5p47BzBu4iQWzJ1Tpsz69et5+OGHsc6ew7Isnnp8CsPHjgz01CqlpFNEROqkXh+5912J/B+gNBY/enj8uAkh/YFzXMJgtzUNp4k23mOla/J9KRrnZA7mGM7b1wxxLKzx2uG97meZvNEVP4abc9aTmZnJxmULiI2tz+Hio5xp1oajHPU3ZL9lZWWRlpbGus830KVTJ59lUlJS2LBhA+eKTrH/QBGd+3Rj8AN3Bz2W85R0ioiIiIiIhNChgwew2+3ExtYHwN6sKbRqxbJ3lvHG395g5vyZAKxcuZIZL81g5ryZVa4jY+lSXp79CnsK9pGRkcH+wkIax8ezKXczy1ev9ijbsGFDAM4Ap78/jTGmeidYiXohPbqIiIiI1AoPnb7VdZPo06DpRNdNal63Hmns3buXK7oP4rEpz/Hx2hwAunTvwq4duyg+XAxAeno6Q0YOCaiOu/r3Z+W7GaSkpJCXl0fSlVew8t0lzHt1ts/y69atI/XWLvzk9q785XcvERMTuvZItXRK0IXqiq3ex/Wn+1JtHL8pApB4eoFrOT98YUSEvkmPuZbzoqSLWKTTOPHIpe60ItGpYaN4cnJyyHpnLh9mf86IRycz43gMnQd0ZsDwAWQuyuSa8dewdu1annjxiYDr2blrF+3atQPg5MlTxMeXP2b0hhtuIHf1erbv+JKHHv9vBt43JGQtkko6RURERETqGP24VPNsNhs9u3WiZ7dOXHNVB+YtXkznAZ0ZPHIw40eNJ7F5IsOHDw+4xbHPoCEUHz2KBSQnJ1O4bx+3DbyTaU8/xaD2Hcq9X1KHK4lvFM/WrVvpmHBlgGdXMSWdIiIiIiIiIZS/cwdxp5rQoZFjPXfbl7Rt2xaAlgktaZHQgmnTprFq1SrOcS6gOpYveZspz0xl/OMTycvLo2h3PmPuHeWz7K5du2jTpg0Auwv28OXXX5GYmAinA6q6Uko6RUREpE5Qy46IhMvJEycYPXoCJYeLiImxcXliG16dv5ADHACg/7D+LD62mKSkJLYd3lbhsVJTU1mR8YnPfXnbviA1NZVZs2Zx/7Ch5R5jzZo1zJgxgxhs1KtXjz9PfxG73c6ZguOBn2QFlHRKRAjVOFAREak5dfUCJbp+gEj0Wdn9bQAS2nfg1Natru1xKSkAFO3c4dqW0L5DmfUtBSVVqi+5YyrZ2dlQuKl0o93OgcOOpHPTZ5sYO3asx31WbFyB3W7n4IlvPebozM3N5eDub33Wk7loITabjTlz5njEDFD8fZFr+b777uO+++4LWZLpTUlnhGlR1CPcIYhImEVja8zdU0o/TnQpHzlPF70REanc3bfeTVzDONJnpYc7lJBR0ikiIiIiIhImC1cvBCA2NjbMkYROlZJOY8zFwMtAGmAD1gA/tywrP/ihVZ/3PFKaV0rCzVc3YnUtlmhs2RQRERHxV1VbOv8ObAamArHA/wPeAG7y9wCTkrKqWKUEgx538UVdIuse/c9FRKJbr4/Gua1pDLFEhwqTTmPMNOB/Lcs649yUBAy2LOsH5/7ngM9CG6KISPXoA1pEREQkfCpr6YwDNhljHrUs6xPg38AyY8xi4AJgNLA0xDGKiIiISC3lcfVf0BWARWqhCpNOy7ImGWMWAK8aYzYCvwbuBm4D6gELgL+GPEqRAKkroYiA3gtEgqFv0mOu5Ty9kkSqbPr06SyYPxebrR71TD1emTuf+Pbx5ZbvfX1vNm/cDNTnsuRWHtOmVGT58uVkZ2dzz50DGDdxEgvmzim37J59e0nt1YUnH5/C5Gm/qeop+a3SMZ2WZeUYY7oAk3B0pZ1iWVb5M42KiIigbs0iUntkfXKfa/nWXmEMRILmaP87HX/9Ketj/QLvQv+p+HNuc856MjMz2bhsAbGx9TlcfJQzzdpw1K8IqiYrK4u0tDTWfb6BLp06VVj2V//7a/qk3R70GLz5dSEhy7J+BH5vjFkE/J8x5j7gEcuy9oU0OhGRGqKrW4uIiEioHDp4ALvdTmxsfQDszZpCq1Yse2cZb/ztDWbOnwnAypUrmfHSDGbOm1nlOjKWLuXl2a+wp2AfGRkZ7C8spHF8PJtyN7N89eoy5ZcsWUJim7Y0atiweifnh8ouJHQtjivWXgVsAR60LKuPMeZ+YI0x5gXLsv4S8ihriRZFPcIdgohEgcTTCzzW88MThkiNUIt4+dzfC/LDF4aIBEG3HmnM/+sLXNF9ELfdfAMj7uzNLUOvo0v3Lkz71TSKDxeDHdLT0xkycojHfVtesBPDOVpesBO4rtw67urfnzv79WPo/WPIzs7m5m5dSf+/WcTHl+3C+9133/H888+z9LW3+dMrVU9wq6peJfvnAllAZ+AtYDaAZVnzgS7ADcYYXb1WRERERESkHA0bxZOTk8Orv3+SFs2bMOLRybz22msYYxgwfACZizIpKSlh7dq1dL+1e8D17Ny1i3bt2gFw8uQpnwknwNSpU3n88ceJb1T+mNJgqqx77RXACMuyvjbG7MAxLycAlmUdAu4zxvQOYXwiIhXSBWJERAS8roKrK+BKBLLZbPTs1ome3TpxzVUdmLd4MZ0HdGbwyMGMHzWexOaJDB8+nJgYv0ZAltFn0BCKjx7FApKTkynct4/bBt7JtKefYlD7Dh5l161bx6JFi/jVj7+k5Ntj1DOG+ISLeHjQ6CCcaVmVndFHOK5c+y+gF/CpdwHLslaEIC4RkVqhJpNidcUTERGJTPk7dxB3qgkdGjnWc7d9Sdu2bQFomdCSFgktmDZtGqtWreIc5wKqY/mSt5nyzFTGPz6RvLw8inbnM+beUT7LZmVlAXCm4DjPvvg7GjWMZ/z48ZwpOB5Q3ZWpLOm8H/gNcBewGZgRkihERETCTGMLJdyS7vFvOgQRiT4nT5xg9OgJlBwuIibGxuWJbXh1/kIOcACA/sP6s/jYYpKSkth2eFuFx0pNTWVFxic+9+Vt+4LU1FRmzZrF/cMiZ8KRyubpPAr8ooZiERGpE/J27amRetT1WEQkNBo0nVjl+9TWi8TN7jrBtTwuij5tmi59F4CE9h04tXWra3tcSgoARTt3uLYltO9QZn1LQUmV6kvumEp2djYUbirdaLdz4LAj6dz02SbGjh3rcZ8VG1dgt9uhcC/f7SjtcJqbm8vB3d/6rCdz0UJsNhtz5szxiLkiT038dZXOJRCBdRgWiVD6ki3hptYyEZHopjk5pabdfevdxDWMI31WerhDCRklnSIiIiIiFYikuZz1A3vts3D1QgBiY2PDHEnoKOkUkTonkr48SNXpC1f5JiVlhTsEERGp5bacu4wDlkXf0wv87iaupFMkCulLt4iIiIhECyWdIiIiIlJnNd6+IdwhiNR61U46jTFzgWWWZS0MQjwSZXR5dwF4c9fzruVJ3BzGSEQkVIJx1eVxCYO9thyr9jFFRCTyBaOlsxdwjzHmEcuydI0vpxZFPcIdgtSwmky8amrKDRERkYro+46I/6ZPn86C+XOx2epRz9TjlbnziW8fX2753tf3ZvPGzdiB+A43eUybUpHly5eTnZ3NPXcOYNzESSyYO6dMmfz8fJKSkriiXQcAulzfmb/9o2y5YKl20mlZVqIxJg5IC0I8IiIiIiIiIbX4D3udS3u99nzgo7R3Ge916Di74ra3zTnryczMZOOyBcTG1udw8VHONGvDUY76G7LfsrKySEtLY93nG+jSqVO55dq3b8/ny/xLZKsrKGM6Lcs6BbwfjGOJiIiISPjpSt+lNIxEquvQwQPY7XZiY+sDYG/WFFq1Ytk7y3jjb28wc/5MAFauXMmMl2Ywc97MKteRsXQpL89+hT0F+8jIyGB/YSGN4+PZlLuZ5atXB/V8qqqeP4WMMbcYY25wWx9jjFljjHnFGFN+m7CIiIjUSpOSslw3CZ/E0wtcNxHxz0Wnvq/xOrv1SGPv3r1c0X0Qj015jo/X5gDQpXsXdu3YRfHhYgDS09MZMnJIQHXc1b8/K9/NICUlhby8PJKuvIKV7y5h3quzfZbftWsXXe7ozm3D+rJmXXZgJ+Ynf1s6XwKmAhhjrgReAeYA3YE/AI9WJwjPCwvoogIiIiLRSJ/nIiK+NWwUT05ODlnvzOXD7M8Z8ehkZhyPofOAzgwYPoDMRZlcM/4a1q5dyxMvPhFwPTt37aJdu3YAnDx5ivh43+2Dl1xyCXv27KHxqfps3LKJ4Q+NZNvtX9AAE3DdFfE36byc0ukAhwIrLct6zNn6uZhqJp0iItFAX6hFRCQSNWg6MdwhiB9sNhs9u3WiZ7dOXHNVB+YtXkznAZ0ZPHIw40eNJ7F5IsOHDycmJrARkH0GDaH46FEsIDk5mcJ9+7ht4J1Me/opBrXv4FE2NjaW2NhYzhQc5/qO19Gu7WV89dVXdEy4MghnWpa/Z3QOsDmXbwXecS4XAc2DHZSIiIiISE1QF3GpCfk7dxB3qgkdGjnWc7d9Sdu2bQFomdCSFgktmDZtGqtWreIc5wKqY/mSt5nyzFTGPz6RvLw8inbnM+beUT7LHjp0iGbNmgHwze5dfL1rp6OF9GRAVVfK36Tzc+ApY8xK4GbgYef2RGB/COKqlXp9NM5tbXvY4gDNrym1h6aPEanY3VNKP+rzKigXTeryRV3cx27mhy8MEamikydOMHr0BEoOFxETY+PyxDa8On8h9jOOK+E+MugOXjp2mqSkJLYd3lbhsVJTU1mR8YnPfXnbviA1NZVZs2Zx/7Ch5R7jk08+4emnnyYGG/Xq1ePl516iWbNmnDl5PPCTrIC/Sef/AxYAdwHTLcva6dw+HFgbgrhERKSWU3cwEREJlwfHtwQgLiWFU1u3urbHpaQAULRzh2tbQvsOZda3FJRUqb7kjqlkZ2dD4abSjXY7FDqSzjXrNzF27FiP+6zYuAK7s4z7HJ25ubkc3P2tz3oyFy3EZrMxZ84cj5i9DR06lKFDh3KmIDRJpje/kk7LsrYCHX3s+gXwY1AjEhERERERqSN+csdIGjWM44XZ88MdSshUaZSqMaYdkAxYwHbLsr4JSVQS1V7rt9u1PCmMcUSr2tgVTiRa6fUoIiKhlrPM2W0+NhaAH0//VxijCQ2/kk5jzIU4pkgZCq6RrcYYsxj4mWVZNdMuW0WRPuZjdtcJruVx+jojIuKXcI7j7bd5Z+WFJCwi/TNfRKQu87el8884utemAednDr0JmI1jDs+fBT0yEZEa1DfpMddynn4EEhEREQkaf5POO4FBlmW5X1P6I2PMwzimT1HSWUfoqrciIiJ1U6RdLVyt21XzhyanXMvjKigXLvp/1m71/CwXBxzxsb0YaBC8cERERERERKQ28Tfp/BR41hjT8PwGY0wj4LeUdrcVERERERERH6ZPn87VacPoeNvdpN5+D+vWrauwfN+uHTl8+DAA8R1u8rue5cuX88wzz3C0pISRD5bfIXXLli30uOtWUm/twvW33cjp06f9rqOq/O1e+ziwHNhnjNni3HYNcBLoE4rARES8qXu3iIiIBMOsZycH9Xgd38yscP/mnPVkZmaycdkCYmPrc7j4KGeatQEOBDUOgKysLNLS0lj3+Qa6dOrks8wPP/zAvffey9w/zqZj8jUcOXqECy64gB85G/R4oArzdBpjOgCjgKucm/8BvGFZ1qny7ykiIiIiIlK3HTp4ALvdTmxsfQDszZpCq1Z8sPA9Zs79F0vmvgjAypUree6FP/PS31+vch0ZS5fy8uxX2FOwj4yMDPYXFtI4Pp5NuZtZvnq1R9kVK1bQsWNHOiZfA0Dzps2x2Wz8WM3zLI/f83RalnUS+FuI4qh1en3kPkR7e8DH0ZyXIiLRL5Ja6aNx7tFJSVmVFwqjaHxMRaRmdeuRxvy/vsAV3Qdx2803MOLO3twy9DrSburMY7+ewaEjR2nRCtLT0xk0YlSVj18v5mIG3/Ugg+58gMEj+5Kdnc3N3bqS/n+ziI+PL1P+q6++whhD/1GDOFx8hOF3DuXX05+i5EfLVaZltc7YU7lJpzFmCPCeZVlnncvlsizr7SDGJBKwSLuyntQ9+vIpItHoodO3hjsEkVqtYaN4cnJyyHpnLh9mf86IRycz43gMY3pfy31D+/H64qU80OZ61q5dyy9+NzPgenZ+8zXt2rUD4OTJUz4TTnB0r12zZg2fZnxIw7g47rhnIDf06sY1l3cOuO6KVNTSuQhIAA46l8tjAbZgBiUiIiIiEmzjEga7rR0LWxxSN9lsNnp260TPbp245qoOzFu8mDG9r+WBEXcxcMwEGlx8OcOHDycmxu/OqB56D7yF4uJizvEjycnJFO7bx20D72Ta008xqH0Hj7KtW7emR48e2Js1B+COtN5s3Lix5pNOy7Lq+VqOZpr/R0QAMkpKB8lH4lxloaTeACIiIjUvf+cO4k41oUMjx3ruti9p27YtAK0SWtDq4hZMmzaNVatWBXwpnxXvfczkpyYxbsIj5OXlUbQ7nzH3+u6q26dPH37/+99z8tRJ6l9Qn0/Wfcqkyb8IsObK+ZVMGmN6GGPKJKjGGJsxpkfwwxIREREREakdTp44wejRo0nuOZSOt93NFzu+YerUqa79o4b0o02bNiQlJVV6rNTU1HL3bdm6mdTUVLKysujapUu55Zo2bcrEiRPpNqAnnfvcxHUp19K/f/+qnFKV+Nt2+yFwCY6utu6aOPepe62IRLVQtQDW5ZZFjW8VEZFI9dhTMwCIS0nh1Natru1xKSkAFO3c4dqW0L5DmfUtBSVVqi+5YyrZ2dlQuKl0o90OhXsBWLN+E2PHjvW4z7/XbsFubwKFe/lux6eu7bm5uRzc/a3Pet5/ZxU2m405c+Z4xOzLvffey90976rSeQTK36TT4Bi76a05cCJ44YiIiIiIiARHv807ATh79ocwR1K+n9wxkkYN43hh9vxwhxIyFSadxph3nYsW8Lox5nu33TYgBcgOUWx1kga41z3n3wxFxD+Jpxe4lvPDF4aIiIvnVHlQnenypO7JWeb8XIuNDW8g5Wh5gft31esCOkZlLZ1HnH8NcBQ45bbvDLAGzd0pIiIiIiIi5agw6bQs6wEAY0w+8EfLstSVVqSOiaRJ7aNBXR7DKSISaTxbIGu29VGfnyKl/BrTaVnWb0MdiIiIRKbj22eEO4QaoQsfiUg00I+bEo38nnnUGPMA8FPgUqC++z7LstoFOS4RcXL/Igzh/zKs+W5Fgk8tIiISrfxJgmd3neBaHufnNxl93yi15dxlruWOYYyjOvydp/OXwAtADpAILAG2As2AuSGKTURERESiWNI9hR43kbps+vTpXJ02jI633U3q7fewbt26Csv37dqRw4cPA3Djla39rmf58uU888wzHC0pYeSDP/NZ5o033iA1NZXOfW6ic5+baHDpReTm5vpdR1X529I5FnjYsqxFxpjxwF8sy/rGGPMU0DZk0YmISJ3Wd9AfPdbz/byfup+JiEhFjrx+1LmU5bXHex0KKCqz3sy70IyKW2M356wnMzOTjcsWEBtbn8PFRznTrA1woCph+yUrK4u0tDTWfb6BLp06+SwzatQoRo0axZmC42zdvo1hD/2U1NTUcuf/rC5/k87WwHrn8ingQufyP53bx/q6Uyhp7E3o6RdJCaaamhomlO8NmqpDREREAnHo4AHsdjuxsY5RivZmTaFVKz5Y+B4z5/6LJXNfBGDlypU898Kfeenvr1e5jiXvLWbmrD+xpyCfjIwM9hcW0jg+nk25m1m+enW593szYxF33zkssBPzk1/da4EiwO5c3g10dS5fjmMOTxERqcPydu1x3ST6JZ5e4LqJRKp+m3e6biKRrluPNPbu3csV3Qfx2JTn+HhtDgBpN3XmP1/nc+iIo+U1PT2dQSNGBVTHoIFDWf1+FikpKeTl5ZF05RWsfHcJ816dXeH93npvMSPuioyk8wPgTufyHOBFY8yHwJvA26EITEREREREpDZo2CienJwcXv39k7Ro3oQRj07mtddewxjDfUP78fripZSUlLB27Vq6p90ecD07v/madu0c13g9efIU8fHxFZZfv+lzGsY15OqrkgOu0x/+dq99GGeCalnWbGPMUeAmYDHwSohiExERCYi6QouISKSx2Wz07NaJnt06cc1VHZi3eDFjel/LAyPuYuCYCTS4+HKGDx9OTIzfE4x46D3wFoqLiznHjyQnJ1O4bx+3DbyTaU8/xaD2HXzeZ2FG6Fs5wY+k0xhzATAd+CuOrrVYlvUmjlbOkGjQdGKoDi11jLr6iUhtNimp7AUvvOl9UEQk/PJ37iDuVBM6NHKs5277krZtHddjbZXQglYXt2DatGmsWrWKswHWseK9j5n81CTGTXiEvLw8inbnM+be8rvqnjt3jrcyF7Hknwso/r6IBBoHWHPlKk06Lcs6a4x5DJgVsihEREJIFx4TkWigFnrR51XtdfLECUaPnkDJ4SJiYmxcntiGV+cvhDN7ARg1pB+Hji8hKSmJLQUlFR4rNTWVFRmf+Ny3ZetmUlNTmTVrFvcPG1rhcT755BNaJVxC20svDeicqsLfttvlQC80J6eIiIiISI1RIhoaze9tCkBcSgqntm51bY9LSQGgaOcO17aE9h3KrHsnhpXNopncMZXs7Gwo3FS60W6HQkfSuWb9JsaO9ZwQ5N9rt2C3N6GwoITPvixwbc/NzS13apP331mFzWZjzpw5HjH70rNnT5YufquSyIPD36RzNfA7Y0xHIAc44b7TsixdTEhERERERKSKfnLHSBo1jOOF2fPDHUrI+Jt0/sX59+c+9lmALTjh+KYxnrWfP+OSwmnEv9yGME+dGrY4Qi3QX1N1uXqpy/T8FxGR6shZ5uxaHxvr930aH3cfr58S3IBCwK+k07Isf6dWERHxi76oRwZdZEZEJHz0Hix1RWDX4xURERERqSFv7nretTyJm8MYiVRmXMJgry3HwhKHRJaoTTqj8ZehujIQPBr/NyKB8Pxg1YdqMKgFvHxJ9xSGOwQRqeX0HU5CJWqTTolukT6G0x8Pnb413CGIiNQ6Sq4F1LIpNSuuWaAzY4q/lHSKSESItC+ami9PREREgmn69OksmD8Xm60e9Uw9Xpk7nxva1C+3fN+uHdm8aSMQw41XtvaYNqUiy5cvJzs7m3vuHMC4iZNYMHdOmTJnz57loYceYv1nn/HDjz8wfNAgfvfHFwI9tUop6RSRoPOnK3mkJZkiUveEsyuh3gOlLomE19r2hj94bH9+0SLHwvm/LhkB1dOxktkNNuesJzMzk43LFhAbW5/DxUc506wNcCCg+iqSlZVFWloa6z7fQJdOnXyWeeutt/j+++/58P1MTp46xS139OPh8f9DQ9Ms6PFAFZJOY0wDYADQHnjFsqwSY0x74KhlWcUhiU5ERKKCWoZFRCKHphuMPIcOHsButxMb62jZtDdrCq1a8cHC95g5918smfsiACtXruS5F/7MS39/vcp1LHlvMTNn/Yk9BflkZGSwv7CQxvHxbMrdzPLVqz3KGmM4ceIEP/zwA6dPn6b+BRdw4YUX8sPx6p+rL35NhWKMuRzYDswGpgPnU+BHgd+HJjQRkciWeHqB6yYiItDro3Gum4iU6tYjjb1793JF90E8NuU5Pl6bA0DaTZ35z9f5HDpyFID09HQGjRgVUB2DBg5l9ftZpKSkkJeXR9KVV7Dy3SXMe3V2mbLDhg2jUaNGXNv1Jjr16MkjDz1Is2ahaeUEP5NO4CVgJXAxcMpt+7tAWpBjEpEINikpy3UTERERkco1bBRPTk4Or/7+SVo0b8KIRyfz2muvYYzhvqH9eH3xUkpKSli7di3d024PuJ6d33xNu3btADh58hTx8fE+y61fvx6bzUZu9hrWf/QBr8xJ55tvvgm43sr42722G3CjZVk/GmPct+8BWgU9KqmSuvLlX1MpiIhEF3W7LlVXpk0TkfLZbDZ6dutEz26duOaqDsxbvJgxva/lgRF3MXDMBBpcfDnDhw8nJiawy+70HngLxcXFnONHkpOTKdy3j9sG3sm0p59iUPsOHmUXLFjAHXfcwQUXXIC9eXM6/+R6NmzYQM8b7gjGqZZRlTO6wMe2S9HkdCIiIiIiIuXK37mDuFNN6NDIsZ677Uvatm0LQKuEFrS6uAXTpk1j1apVBDqBy4r3PmbyU5MYN+ER8vLyKNqdz5h7fXfVvfTSS/nggw+4vduNnDx5kpxNuUx56ukAa66cv0nnCmAi8DPnumWMuRD4LbA0FIGJiIiISPi0bjDAbU1tDBI8r/Xb7bE+KUxx1KSTJ04wevQESg4XERNj4/LENrw6fyGc2QvAqCH9OHR8CUlJSWwpKKnwWKmpqax93feFhrZs3UxqaiqzZs3i/mFDyz3GuHHjeOCBB7jljn5YlsU9w4bSsWNHDu7+NuBzrIi/SedE4ENjzJdAA+BN4HIc1/i9OySRiYiIiIj4SV2YpSqeGDYMgLiUFCjcVLqj1XUAFO3c4dqU0L6DRyLYsXWTShNDb8kdU8nOzvasy26HQkfSuWb9JsaOHetxn3+v3YLd3oTCghKPOTpzc3M5tXWrz3ref2cVNpuNOXPmeJyDt/j4eN56660KywSTX0mnZVmFxphU4KfA9TguQPQq8IZlWacqum9tMC5hsNuafumLdvpQEql7wjlHnIiIVF1tvWbJlnOXuZY7Ov/+5I6RNGoYxwuz5/t9nLhmVe+EmxD3XZXvEyx+j+l0JpdznTephTRRdd1TV76I15XzlKqpqz9A6fUQHnrcRaQ8OcucF12LjQ1vICFUbtJpjBni70Esy3o7OOGIiIiIiEio6erSUpMqaulc5OcxLMAWhFhEREREJAx00SARCaVyk07LsurVZCAiIhJcdaU7nz/nqV/0SwX6WETSY1hbx3qJiNRWSixFREREREQkZPxOOo0x1xtj5htjNjhv/zDGXB/K4EREREREatJDp2913USCafr06Qy+tSvDbr+Ju/vczLp16yos37drRw4fPgzAjVe29rue5cuX88wzz3C0pISRD/7MZ5kzZ87wwAMPkNZvALcOGEj2ZxXHUl1+Xb3WGDMKmA98ALzv3HwjsN4YM8ayLN+zkwZJr4/GeW3ZHsrqIkpd6R4nIiIiIlJTsg/e5Vj4wGvHf8qW3bbbc331Vz4O2HpnhfVtzllPZmYmb77/EfVjYzlafIQ29jjggL8h+y0rK4u0tDTWfb6BLp06+Szzt7/9DYAP38/k8JEjjHzwIQb9dGTQYznP3ylTpgNPWZb1O/eNxpgpwDQgpEmniIiIiEgoRNJ45UDN7jrBtTwuxJNAqUEkMIcOHsBut1PfOS1K02bNadWqCR8sfI+Zc//FkrkvArBy5Uqee+HPvPT3qqdXS95bzMxZf2JPQT4ZGRnsLyykcXw8m3I3s3z1ao+yX3zxBb169QLA3rw5F13YmA0bNpB48VXVPFPf/O1e2wJY6GP7W0DL4IUjIiIiIiJSu3TrkcbevXsZ2KMT0389iQ1rPwUg7abO/OfrfA4dOQpAeno6g0aMCqiOQQOHsvr9LFJSUsjLyyPpyitY+e4S5r06u0zZa6+9lnfffZcffviBPXv3smXrNvbu3Rv4CVbC36TzQ6Cnj+09gY+DFYyIiIiIiEht07BRPDk5OTz9/Es0bW7nV+Me5LXXXsMYw31D+/H64qWUlJSwdu1auqfdHnA9O7/5mnbt2gFw8uQp4uPjfZZ78MEHad26NXcMHsLT035Hp+uvw2YL3SyY5XavNcYMcVv9N/CcMaYT8Jlz243AEGBqyKKLYndPKX1oQ9vJQSJFv80V9+WXyOE5TrzujBEXCaZI75KoaVVEJNLYbDY6d+1O567d6XBVMosXL2JM72t5YMRdDBwzgQYXX87w4cOJifF3BKSn3gNvobi4mHP8SHJyMvsLCujdfyAvTplCQvsOHmVjYmL405/+RNHOxwAYOHwEV1xxRbXPsTwVndEiH9sedt7cvQzMClpEQaQkoHZzT+xByb2IiIjUjD80OeVa9r7cpYgv+Tt3EHeqCcS1AODLbXm0bdsWgFYJLWh1cQumTZvGqlWrOBtgHSve+5jJT01i3IRHyMvLo+Trr/nve+7xWfbkyZNYlgXAx2s+xRZjIzk5mYO7vw2w9oqVm3RalqU5PEUk6ugCByIiUh1v7nretTyJm0NaV6T3GJDgOXniBKNHT+DA4WJsNhttEtvx5j/S4YxjHOWoIf04dHwJSUlJbCkoqfBYqamp5L6f7nPflq2bSU1NZdasWTw2cGC5xzh48CB9+vTh3I8/cMnFF/PyH/8Q8Ln5I7C2WxEpI+mewnCHICJBph8xRERqp24tMwCIS0mBwk2lO1pdB0DRzh2uTQntO3gkgh1bN6k0MfSW3DGV7Oxsj/vZ7U2g0JF0rlm/ibFjx3rc599rt2C3N6GwoITPvixwbc/NzfWM2c3776zCZrMxZ84cTm3dWm48iYmJfPnllx7nGUp+J53GmKZAX+BSoL77Psuy/jfIcUWdYI0R0y9eEk76gi3BFMjzSePhRURqv0j7vrHl3GWu5Y5hqP8nd4ykUcM4Xpg9Pwy11wy/kk5jzI3AUuB7HNOn7AMuca7nA3U+6RQRkdpPSbFIdFEvJIkGOcucjU7OOTxrI39bOv8AvAFMAL4FegEngH8Cc0ITWng0aDox3CGIiIiIiIjUGv4mnR2Bn1mWZRljfgRiLcv6xhjzBLAAR0IqEpX6Dvqjazm/Cvc7vsTtQs4zNO2GSGVmd53gWh4XxHbCSBuWEOh7Sl2llijRc0Ck9vM36TzjtnwAaItj4OJ3QKtgB+VNb0ZSl0XauAcRERERqTuCMebV36RzI9AZ+Ar4CJhmjLkYuBfYEmDdIiIiIiJlTErKCncI1RZpvTBEwsnfpPM3QGPn8pPAfOBlHEnoAyGIq8rcL+4AusCDhF6kt8DXZAupLq4iUj3uX05BX1BFRGqj6dOnM3f+69jq1aNevXrMm/t3bmhTv9zyfbt2ZPOmjUAMN17Z2mPalIosX76c7OxsHu7dmweeeIKM2bPLlDly5AjDhg1j/fr1jBgymN9Nfca1b3PeJn7+i8c4+8P39OvXj5/98rcYY6p8vu78Sjoty9rgtnwIx9QpIiIiIiIiUeeyQz84Fj7M9dzha9rKPZWXKWqdWmF9m3PWk5mZyZvvf0T92FiOFh+hjT0Ox8jF4MrKyiItLY1Pc3Lodv31Pss0aNCAZ599lk8//IAvv/rKY9+vnpzIC8/NpO+gXvTr149PP1pF97TbqxWT3/N0ikQajXWUQKllVkRCSb2vRMTboYMHsNvt1HdOi9K0WXNatWrCBwvfY+bcf7Fk7osArFy5kude+DMv/f31Ktex5L3FzJz1J/YU5JORkcGBwkIaN2rE53l5LP3wQ4+yjRo1onv37uRkf+qx/cDBIr47fpxO13fGGMP999/P4qVLQ5d0GmPyAMufg1iWFY55VEVEpAIaTyQiIhIZuvVIY/5fX2Bgj07c2P0W+gwcQsfh/Um7qTOP/XoGh44cpUUrSE9PZ9CIUQHVMWjgUO4aMITBI/uSnZ1N2o038uaf/0zjRo38Psb+okIuuaT0OrGtW7fmYNH+gOJxV1FL56JqH11ERETESWNXJdz0Y5wAxDU7C4eqN0axqho2iicnJ4fX3v43n2dn8atxD1LvxPOM6X0t9w3tx+uLl/JAm+tZu3Ytv/jdzIDr2fnN17Rr1w6AU2dO0LJNfeBskM4icOUmnZZl/bYmA4lk/rxBqbte9Ij0brmRHp+IRA59gZZIE+kX2RMJJ5vNRueu3enctTsdrkpm8eJFjOl9LQ+MuIuBYybQ4OLLGT58ODExgY2A7D3wFoqLiznHjyQnJ7O/sIDU2+/h5Wm/4ubB1/l1jEsSWrF/f+nruKCggJYJlwQUjzuN6QwjJRciIiIiwRHK71WVTeESymRb3xdrh/ydO4g71QTiWgDw5bY82rZtC0CrhBa0urgF06ZNY9WqVQG3S65472MmPzWJcRMeIS8vj5OFX/LYmLurdIyLWyYQ37gxGzZ+Tt9LezF//nwG/LT6k5VoTKeIBJ0+IEVE6obZXSe4lsepr1fE+EOTU67lcWGMQ0qdPHGC0aMncOBwMTabjTaJ7XjzH+lwZi8Ao4b049DxJSQlJbGloKTCY6WmppL7frrPfVu2biY1NZVZs2Yx8d6KJxxJTEzkWEkJZ86eZdnKVaz64APsjVrz/LMv8PNfPMb//OJ7+vbtW+2LCIHGdIaMvnSLSF2nrp8iIhKpiq50tq21us4jyevYugkAp7ZudW2LS0kpU6ayxNBbcsdUsrOzPe5ntzeBQkfSuWb9JsaOHetxn3+v3YLd3oTCghKPOTpzc3OhcJNH2cbHHblH1vy52Gw25syZU6aMt/z8fIp2ls7/ktC+Awd3f0tqx+v5ZMVntGx7IUCVz9UXjekUEREREREJk5/cMZJGDeN4Yfb8cIcSMn6N6TTGvAT83bKsrZWVFREREREREf/kLHP2DHLO4Vkb+Xshoc7A/xhjcoC/A/+0LOt46MKqumB1Zx2XMNht7VhQjikiuqKhiIhIbeQ9lEJTI4kvfiWdlmXdZIy5EngQeAZ40RjzNjDHsqyPQxmgSG2n6XZqnsZcSzD1HfRH13J++MIQEZEKWViWhTE1Oz9nbWRZFpZ/15t1qVeFg39pWdYTQBvgHiAeWGGM2WGMmWyMaValmkVERERERGpAg2PfcOTIESyrasmSeLIsix9OfsvukqpN7BLIPJ0XABcCFwE2YA9wH/CkMeZhy7IWVHRnERERERGRmtR64/MUXNmXQ4cOQcnB0h3HtnPgaOkUM9uPxwFw9sAB17YLbLYyZdzXfW3zt4x3LECFxzkfn/f9vOP1VcZXXd8eKt129MwPHD9y2rV+5GSDMvFccDyOjbuO8fK6o1SF30mnMaYTju619wAngXnAQ5Zl7XLufxT4E6CkU0REpJo0Dlqk+jR1k5x3wZkSLrvsMsfK1BtLd0w9Rt/JS12r+TP6A7B98BDXtqT/bC9Txn3d1zZ/y3jHAlR4nPPxed/PO15fZXzV9cLUX7o2TXozk78+8oFrfdzsXr7jmed5Xv7w9+q1ecCVwHJgDLDUsqwfvYq9Bfy1yhFIrfdav92u5UlhjENEpLqOb58RlOPoi7CIiNQl/rZ0LgTmWpa1r7wClmUdpgpjREWkZuniOSIiEq36bd4Z7hBEpBr8vXrts+7rxpgYoIFlWd+FJCoRERGRavDVmlxXWpjrynmKSPSoMOk0xtwKNLcsa6HbtsnAVCDGGLMKuMeyrJJQBiml3LuqQt3qrqqpRWo3zetVMyqaTy2/hmOZ3XWCx/o4vbJFRET8EqzP75oaBldZS+dk4N/nV4wxXYDfAXOA7cAvgd84/4pEBf0CLBIewRoPKSIinjSERiJdZUnnNTgSz/OGA9mWZY0FMMbsBaahpFNERERqIX2ZD75ovDKzfrAWqZ7Kks4mgNuELtwEvO+2/jnwX0GOSUQXDJA6RV9mIoOSCxEBmJSUFe4QRGqdypLO/UB7YK8xJha4DnjKbX9j4PsQxRZV9GVFRETKox8WQk+JgogEU99Bf3Qt54cvjFqjsilO/g383hjTC3geOAG4v6t3BL4OUWwiIiIiIiIS5Spr6XwaeBtYBXwHjLYs64zb/geBlSGKTURE6hhdbEhERKT2qTDptCzrMNDDGHMR8J1lWT96FRmOIxkVERERERERKaOylk4ALMs6Vs724uCG4zAuYbDbms+qa4W6PManLp+7RCY9J0VERERCo7IxnSIiIiIiIiIB86ulU0Qkmnm3YqpVU0REpHLRODuDPuMjU1Qkne5PHtCXRhFveo2ICOh1LyISberK+3ZUJJ0igaorL2QRERERkUilMZ0iIiIiIiISMmrplIh09xTPp2ZeDdffd9AfXcv5NVx3TQpnS3BtbYWureclUldE4xg2ERFfIuk7iZJOkToikt54RKR20vuMgxLXmuH+A3Ueev6JRDIlnSISsZLuKQx3CCIiUgEl2CLiDyWdIlJt+tIhIhK5fLUAqlVQpHaJ9KFhSjq9hPJNOFRfzCuag/D8NolcSthEREQihxJyqUu8v4fO7jrBtTwuiFdV0dVrRUREREREJGSUdIqIiIiIiEjIKOkUERERERGRkNGYzgD8ockp1/K4ENdVV8YV1JXz9IfG6JbP12Oh546IRAO9V4lIuIXzfUgtnSIiIiIiIhIyaukUEREREakG9UoSqViNJ52aKyp66P8igdLrXGoLPW9FRKKbvpNEBnWvFRERERERkZBR0ikiIiIiIiIhozGdIiIiIuXw1Q3v+PYZYYlFRCRaKekMkoqmuciv8WgiS96uPTVWlx53Eamt9P4mIiLRSkmniIiIiIjoxy0JGY3pFBERERERkZBR0ikiIiIiIiIhU2u61/qalLfvoD96rEvFIr1LRU2ODRURCbdIf08W8Yev72d6bovUPWrpFBERERERkZCpNS2dIiISHJoOQkRERIJJSadUi68uMsHoNqOutCIiIiIitYOSThEREZEg8jWOMZAyIiI1LVS9nTSmU0REREREREJGLZ0iIiIiIhJ1dCXk6KGkU0JObwil+m3eGe4QooqeOyIiIiLRT91rRUREREREJGSUdIqIiIiIiEjIqHttlNFUIiIiUpuoG70ESt+Jopeu3lz3qKVTREREREREQkZJp4iIiIiIiISMkk4REREREREJGY3plCrR2JtS3o+Fr8dmUlJWTYYkIiIiEpX0nal2U9IpIiIiIiJhpYaN2k3da0VERERERCRk1NIpIiIiIiIho1ZMUdIpIiIiIiJSB9XUDwLqXisiIiIiIiIho5ZOEXFR9xcRERERCTZjWVboKzHmELA75BWJiIiIiIhIOLS1LKuFrx01knSKiIiIiIhI3aQxnSIiIiIiIhIySjpFREREREQkZJR0ioiIiIiISMgo6RQREREREZGQUdIpIiIiIiIiIaOkU0REREREREJGSaeIiIiIiIiEjJJOERERERERCRklnSIiIiIiIhIy/x9LZivSzXBR8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check syll ratio matrix\n",
    "session = selected_sessions[0]\n",
    "session_data = load_session_data(session, load_aeps=False)\n",
    "tl = session_data['tl']\n",
    "moseq_file   = session_data['moseq_file']\n",
    "with h5py.File(moseq_file, 'r') as f:\n",
    "    syl_ratio_mx  = np.array(f['syl_ratio_mx'])\n",
    "    idxs_srm_tl   = np.array(f['idxs_srm_tl'])\n",
    "\n",
    "t1, t2 = 0, 200\n",
    "\n",
    "#colors_bold = [plt.cm.tab20(2*i) for i in range(10)]\n",
    "#colors_tran = [plt.cm.tab20(2*i + 1) for i in range(10)]\n",
    "#colors = colors_bold + colors_tran\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(16, 4))\n",
    "\n",
    "# syllable ratios\n",
    "idxs_to_idxs = np.where((tl[idxs_srm_tl][:, 0] > t1) & (tl[idxs_srm_tl][:, 0] < t2))[0]\n",
    "idxs_sel = idxs_srm_tl[idxs_to_idxs]\n",
    "bottom = np.zeros(len(idxs_sel))\n",
    "for i, syl_ratio in enumerate(syl_ratio_mx[idxs_to_idxs].T):\n",
    "    ax.bar(tl[idxs_sel][:, 0], syl_ratio, 0.7, bottom=bottom, label='Syll. # %s' % str(i+1))\n",
    "    bottom += syl_ratio\n",
    "ax.set_xlim(t1, t2)\n",
    "ax.set_title('Behavior composition', fontsize=14)\n",
    "ax.set_ylabel('Syllable ratios, %', fontsize=14)\n",
    "ax.legend(loc='upper right')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00026926",
   "metadata": {},
   "source": [
    "## 02 - Unit corr with shuffled behavior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "15a9c3c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_corr_to_syll_ratios(i_rate_binned, syl_ratio_mx):\n",
    "    data = np.column_stack([i_rate_binned, syl_ratio_mx])\n",
    "    columns = ['state'] + [\"x%d\" % x for x in range(syl_ratio_mx.shape[1])]\n",
    "    syl_ratio_df = pd.DataFrame(data, columns=columns)\n",
    "\n",
    "    model = glm('state ~ ' + ' + '.join(columns[1:]), data=syl_ratio_df).fit()\n",
    "\n",
    "    glm_coeffs = dict([(i, coef) for i, coef in enumerate(model.params[1:]) if model.pvalues[1:][i] < 0.05])\n",
    "    if len(glm_coeffs) == 0:\n",
    "        return 0, 0, model.params, model.pvalues\n",
    "    behav_fit = np.zeros(len(syl_ratio_mx))\n",
    "    for idx, coef in glm_coeffs.items():\n",
    "        behav_fit += coef * syl_ratio_mx[:, idx]\n",
    "        \n",
    "    corr, pval = stats.pearsonr(behav_fit, i_rate_binned)\n",
    "    return corr, pval, model.params, model.pvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "49886b72",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['57_SIT_2024-01-04_14-16-22', '57_SIT_2024-01-04_14-52-59']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sessions = selected_57\n",
    "sessions.sort()\n",
    "selected = sessions[:3]\n",
    "selected = [\n",
    "    \"57_SIT_2024-01-04_14-16-22\",\n",
    "    \"57_SIT_2024-01-04_14-52-59\",\n",
    "]\n",
    "selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b84f375f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1-10; 1-11; 1-12; 1-13; 1-14; 1-15; 1-16; 1-17; 1-18; 1-19; 1-2; 1-3; 1-4; 1-5; 1-6; 1-7; 1-8; 1-9; 2-10; 2-11; 2-12; 2-13; 2-14; 2-15; 2-16; 2-17; 2-2; 2-3; 2-4; 2-5; 2-6; 2-7; 2-8; 2-9; Session 57_SIT_2024-01-04_14-16-22 done\n",
      "1-10; 1-11; 1-12; 1-13; 1-14; 1-15; 1-2; 1-3; 1-4; 1-5; 1-6; 1-7; 1-8; 1-9; 2-10; 2-11; 2-2; 2-3; 2-4; 2-5; 2-6; 2-7; 2-8; 2-9; Session 57_SIT_2024-01-04_14-52-59 done\n"
     ]
    }
   ],
   "source": [
    "iter_count = 100\n",
    "\n",
    "for session in selected:\n",
    "    animal = session.split('_')[0]\n",
    "    units_file = os.path.join(source, animal, session, 'units.h5')\n",
    "    moseq_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')\n",
    "    moseq_units_glm = os.path.join(source, animal, session, 'analysis', 'MoSeq_units_GLM.h5')\n",
    "    \n",
    "    # load syll ratio matrix/tl_idxs\n",
    "    with h5py.File(moseq_file, 'r') as f:\n",
    "        syl_ratio_mx  = np.array(f['syl_ratio_mx'])\n",
    "        idxs_srm_tl   = np.array(f['idxs_srm_tl'])\n",
    "\n",
    "    # special processing for dark sessions - remove dark\n",
    "#     if session_data['animal'].startswith('0082'):\n",
    "#         tl = session_data['tl']\n",
    "#         idxs_light   = np.where((tl[idxs_srm_tl][:, 0] < 600) | (tl[idxs_srm_tl][:, 0] > 1800))[0]\n",
    "#         idxs_srm_tl  = idxs_srm_tl[idxs_light]\n",
    "#         syl_ratio_mx = syl_ratio_mx[idxs_light]\n",
    "\n",
    "    # load units\n",
    "    unit_names, single_units, spike_times = [], {}, {}\n",
    "    with h5py.File(units_file, 'r') as f:\n",
    "        unit_names = [x for x in f]\n",
    "    with h5py.File(units_file, 'r') as f:\n",
    "        for unit_name in unit_names:\n",
    "            spike_times[unit_name] = np.array(f[unit_name][H5NAMES.spike_times['name']])\n",
    "            single_units[unit_name] = np.array(f[unit_name][H5NAMES.inst_rate['name']])\n",
    "            #single_units[unit_name] = instantaneous_rate(unit_times, tl[:, 0], k_width=50)\n",
    "                \n",
    "    # create an empty MoSeq units GLM file\n",
    "    with h5py.File(moseq_units_glm, 'a') as f:\n",
    "        if not 'units' in f:\n",
    "            f.create_group('units')\n",
    "\n",
    "    for name, i_rate in single_units.items():\n",
    "#         if not name == '1-4':\n",
    "#             continue        \n",
    "        \n",
    "        i_rate_binned = i_rate[idxs_srm_tl]\n",
    "\n",
    "        # compute original\n",
    "        corr, pval, params, pvalues = get_corr_to_syll_ratios(i_rate_binned, syl_ratio_mx)\n",
    "\n",
    "        # compute shuffled\n",
    "        corr_mx = np.zeros([iter_count, 2])  # coeff, pval for each shuffle\n",
    "        for i in range(iter_count):\n",
    "            srm = syl_ratio_mx.copy()\n",
    "            np.random.shuffle(srm)\n",
    "            corr_s, pval_s, _, _ = get_corr_to_syll_ratios(i_rate_binned, srm)\n",
    "            corr_mx[i] = (corr_s, pval_s)\n",
    "            del srm\n",
    "\n",
    "        with h5py.File(moseq_units_glm, 'a') as f:\n",
    "            if not name in f['units']:\n",
    "                f['units'].create_group(name)\n",
    "            unit_grp = f['units'][name]\n",
    "\n",
    "            if 'corr_glm_fit_orig' in unit_grp:\n",
    "                del unit_grp['corr_glm_fit_orig']\n",
    "            if 'corr_glm_fit_shuffled' in unit_grp:\n",
    "                del unit_grp['corr_glm_fit_shuffled']\n",
    "            if 'glm_fit_params' in unit_grp:\n",
    "                del unit_grp['glm_fit_params']\n",
    "            if 'glm_fit_pvalues' in unit_grp:\n",
    "                del unit_grp['glm_fit_pvalues']\n",
    "\n",
    "            unit_grp.create_dataset('corr_glm_fit_orig', data=np.array([corr, pval]))\n",
    "            unit_grp.create_dataset('corr_glm_fit_shuffled', data=corr_mx)\n",
    "            unit_grp.create_dataset('glm_fit_params', data=params)\n",
    "            unit_grp.create_dataset('glm_fit_pvalues', data=pvalues)\n",
    "\n",
    "        print(name + '; ', end='')\n",
    "        \n",
    "    print('Session %s done' % session)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a03160",
   "metadata": {},
   "source": [
    "## 03 - Unit / GLM predictions for train / split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "db8ca96f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1-10; 1-11; 1-12; 1-13; 1-14; 1-15; 1-16; 1-17; 1-18; 1-19; 1-2; 1-3; 1-4; 1-5; 1-6; 1-7; 1-8; 1-9; 2-10; 2-11; 2-12; 2-13; 2-14; 2-15; 2-16; 2-17; 2-2; 2-3; 2-4; 2-5; 2-6; 2-7; 2-8; 2-9; Session 57_SIT_2024-01-04_14-16-22 done\n",
      "1-10; 1-11; 1-12; 1-13; 1-14; 1-15; 1-2; 1-3; 1-4; 1-5; 1-6; 1-7; 1-8; 1-9; 2-10; 2-11; 2-2; 2-3; 2-4; 2-5; 2-6; 2-7; 2-8; 2-9; Session 57_SIT_2024-01-04_14-52-59 done\n"
     ]
    }
   ],
   "source": [
    "iter_count = 100\n",
    "\n",
    "for session in selected:\n",
    "    animal = session.split('_')[0]\n",
    "    units_file = os.path.join(source, animal, session, 'units.h5')\n",
    "    moseq_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')\n",
    "    moseq_units_glm = os.path.join(source, animal, session, 'analysis', 'MoSeq_units_GLM.h5')\n",
    "    \n",
    "    # load syll ratio matrix/tl_idxs\n",
    "    with h5py.File(moseq_file, 'r') as f:\n",
    "        syl_ratio_mx  = np.array(f['syl_ratio_mx'])\n",
    "        idxs_srm_tl   = np.array(f['idxs_srm_tl'])\n",
    "\n",
    "    # load units\n",
    "    unit_names, single_units, spike_times = [], {}, {}\n",
    "    with h5py.File(units_file, 'r') as f:\n",
    "        unit_names = [x for x in f]\n",
    "    with h5py.File(units_file, 'r') as f:\n",
    "        for unit_name in unit_names:\n",
    "            spike_times[unit_name] = np.array(f[unit_name][H5NAMES.spike_times['name']])\n",
    "            single_units[unit_name] = np.array(f[unit_name][H5NAMES.inst_rate['name']])\n",
    "            #single_units[unit_name] = instantaneous_rate(unit_times, tl[:, 0], k_width=50)\n",
    "        \n",
    "    # second do train GLM on 50/50% of chunk-shuffled data\n",
    "    chunk_size  = 5  # seconds, syll_ratio_mx is sampled 4 Hz\n",
    "    chunk_count = int(len(syl_ratio_mx)/chunk_size)\n",
    "    split_ratio = 0.5\n",
    "\n",
    "    for name, i_rate in single_units.items():\n",
    "#         if not name == '6-3':\n",
    "#             continue\n",
    "        \n",
    "        i_rate_binned = i_rate[idxs_srm_tl]\n",
    "\n",
    "        # split into chunks\n",
    "        syll_chunks   = np.array(np.array_split(syl_ratio_mx[:chunk_count*chunk_size], chunk_count))\n",
    "        i_rate_chunks = np.array(np.array_split(i_rate_binned[:chunk_count*chunk_size], chunk_count))\n",
    "\n",
    "        # compute chunked\n",
    "        corr_tr_mx = np.zeros([iter_count, 2])  # coeff, pval for each random chunking\n",
    "        glm_cf_mx = np.zeros([iter_count, syl_ratio_mx.shape[1] + 1])  # + intercept  !!\n",
    "        glm_pv_mx = np.zeros([iter_count, syl_ratio_mx.shape[1] + 1])  # + intercept  !!\n",
    "        for i in range(iter_count):\n",
    "            # random indices\n",
    "            idxs_ch_rand = np.random.choice(np.arange(chunk_count).astype(np.int32), chunk_count, replace=False)\n",
    "\n",
    "            # randomized behav / unit firing\n",
    "            i_rate_rand = i_rate_chunks[idxs_ch_rand].flatten()\n",
    "            syll_chunks_rand = syll_chunks[idxs_ch_rand]\n",
    "            syll_chunks_rand = syll_chunks_rand.reshape(-1, syll_chunks_rand.shape[-1])\n",
    "\n",
    "            # build train, split\n",
    "            idx_split = int(split_ratio*len(syll_chunks_rand))\n",
    "            data = np.column_stack([i_rate_rand[:idx_split], syll_chunks_rand[:idx_split]])\n",
    "            columns = ['state'] + [\"x%d\" % x for x in range(syll_chunks_rand.shape[1])]\n",
    "            syl_ratio_df = pd.DataFrame(data, columns=columns)\n",
    "\n",
    "            if (syl_ratio_df['state'] == 0).all():\n",
    "                continue\n",
    "\n",
    "            model = glm('state ~ ' + ' + '.join(columns[1:]), data=syl_ratio_df).fit()\n",
    "\n",
    "            glm_cf_mx[i] = model.params\n",
    "            glm_pv_mx[i] = model.pvalues\n",
    "\n",
    "            glm_coeffs = dict([(i, coef) for i, coef in enumerate(model.params[1:]) if model.pvalues[1:][i] < 0.95])\n",
    "            if not len(glm_coeffs) == 0:\n",
    "                behav_fit = np.zeros(len(syll_chunks_rand[-idx_split:]))\n",
    "                for idx, coef in glm_coeffs.items():\n",
    "                    behav_fit += coef * syll_chunks_rand[-idx_split:][:, idx]\n",
    "                if not (np.diff(behav_fit) == 0).all():\n",
    "                    corr_tr_mx[i] = stats.pearsonr(behav_fit, i_rate_rand[-idx_split:])\n",
    "\n",
    "        with h5py.File(moseq_units_glm, 'a') as f:\n",
    "            if not 'units' in f:\n",
    "                f.create_group('units')\n",
    "            if not name in f['units']:\n",
    "                f['units'].create_group(name)\n",
    "            unit_grp = f['units'][name]\n",
    "\n",
    "            if 'corr_glm_fit_train_test' in unit_grp:\n",
    "                del unit_grp['corr_glm_fit_train_test']\n",
    "            if 'glm_cf_mx' in unit_grp:\n",
    "                del unit_grp['glm_cf_mx']\n",
    "            if 'glm_pv_mx' in unit_grp:\n",
    "                del unit_grp['glm_pv_mx']\n",
    "\n",
    "            unit_grp.create_dataset('corr_glm_fit_train_test', data=corr_tr_mx)\n",
    "            unit_grp.create_dataset('glm_cf_mx', data=glm_cf_mx)\n",
    "            unit_grp.create_dataset('glm_pv_mx', data=glm_pv_mx)\n",
    "        print(name + '; ', end='')\n",
    "    \n",
    "    print('Session %s done' % session)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2334e37",
   "metadata": {},
   "source": [
    "## Test - shuffled VS original syll ratios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "a4feeb51",
   "metadata": {},
   "outputs": [],
   "source": [
    "unit = '6-3'\n",
    "\n",
    "session = list(selected_009265.keys())[10]\n",
    "\n",
    "session_data = load_session_data(session, load_aeps=False)\n",
    "moseq_file = session_data['moseq_file']\n",
    "\n",
    "with h5py.File(moseq_file, 'r') as f:\n",
    "    grp = f['units'][unit]\n",
    "    corr_glm_fit_orig       = np.array(grp['corr_glm_fit_orig'])\n",
    "    corr_glm_fit_shuffled   = np.array(grp['corr_glm_fit_shuffled'])\n",
    "    corr_glm_fit_train_test = np.array(grp['corr_glm_fit_train_test'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "99d14146",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Unit 6-3 (0.42)')"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAADTCAYAAACIhgYXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAddElEQVR4nO3deXwUdZ7/8ddHroBEwKCRQwlegCQSQkAjHkGZH+I5qCCsKNlVcWCcWdcL0dHBY0Zc8RoRFJ0Rxl1FEZlRf+gqah4zogiSDYIgcgUJiBgUkgABDJ/9oyqxKnQnnaS7qyGf5+PRj1RXfavq3ZXOJ1XfrqoWVcUYY6ocEXQAY0xisaJgjPGxomCM8bGiYIzxsaJgjPGxomCM8bGiYIzxsaJwiBGRmSLydtA56iIiLURktYicmwBZjhWR70Wka9BZDgVWFOJARPJFZGqI8XkiUl7Pxf07MLquZYfJ0UlEZrl/IBUislJEzqtjnudFZJ2I7HHn+7uI9IpgdWOBLar6D8+yOojISyKy0328JCLtI8nuzv+ciKiI3O4Zd7SIPC0iX7kZN4nIdBFJqWqjqtuAvwL3R7qupsyKwiFGVXeq6o76zuf+8S0EBLgY6AX8BthWx6yfA3lu+yHu/AtEpEUt6xLgt8Cfa0x6GcgCLnQfWcBLEea/ChgAbKkxqTPQBbgTyMApmOcCr9Ro9yJwjYgcHcn6mjRVtUeMH0A+MDXE+Dyg3PN8JvA2zt7AZuBHnDdzm5ptPMNa45EWJsMfgYVReC2nu+vpUUubbOAA0N4zrpc730DPuLPrWpbbrpu7PXoBRcDtdbS/yF3/UTXGrwduCPr9kOgP21NIPOcA6cBg4GpgGE6RCOXfgU9xCkcn97EpTNtfAp+JyKsisk1ECkXkZve/ekRE5EjgX4FvcP44a3sN69S/R5MDlAOfeMYtBHYBZ9WyzuY4//UfUtVVEUY9CtgL7K4xfjFQ6+GSscOHRFQK/EpVV6nqe8Ac4IJQDVV1J7AP2K2qW91HZZjlngiMx/lvOQR4CpgM/LquQCIy3u37KAeGAheo6t5aZunGwbv5xwHfq/sv282vOIcvx9WyrPuBElWdXldON2t74EHgeVX9qcbkLUBaJMtpyqwoJJ6VNf6wtwDHRmG5RwAFqjpRVf9XVV8E/oRbFETkbhEp9zxO8Mz730BfnP+yXwNzRKRNLetqDVQ0NrCI5OIcYl0fYfu2wFs4hxp3hmiyx81mamFFIT5KgXYhxrcHdtYYt7/GcyU6v6dvgZU1xq0Cqv74nwUyPY/q//TqdG6uUeeThKuAU4Era1lXCdChxritwDHewxV3+Fh3Wii5OIdE34rITyLyE85eyCMiUuxt6BaE+e7TS1Q1VFE6Gvi+ltwGKwrxshrICnH8nuVOa4x9QLMI2i0EetQYdyqwEUBVf1DVtZ5HzV3vKuI+WtWyrv8FeoiI9/31KdAWp2+hSg5wJP5+Bq9pOB2bmfiL1RN4DqlEJBl4F2c7XKSq4T7mTQcKasltsKIQL9NxjumfFpE+ItJDRP4DGAU82shlFwEDRCRNRDrW+EP0egI4U0TuEZGTRWQ4zseGz4RbsNtugoj0E5ETROQsnD6OvTifkoTzEZCE8wcNgNtJ+C7wnIjkiEgO8BzOJymr3fV1cc83GObOs01VV3gfOHtSWz3zJAPv4eyZ5AFHishx7qOl57W0Afq5GUwtrCjEgaqux/ns/BScN/BiYCQwXFXfaeTip+DsLazE2TU+IVQjVV2C8wnECGAF8AfgXpz/xuHsxdmFfwdYC7wKlAE5qhpulx9V3Q68AVxTY9K/AMuA/3Efy4BrPdNb4OzNhDrUCqcfcCZwGk5/x7eeh/dTjcuBb1T1n/VYdpMkns5gY6JGRHrj7DGcrKqlCZBnMfCkqr4cdJZEZ3sKJiZU9UvgdqB70FlE5FjgdQ4+y9GEYHsKxhgf21MwxvhYUTDG+DSP58o6duyoaWlpIaft2rWLI488Mp5xwrIsoZWWltKsWfhTIpo1a0ZSUlLMczRqm+zYAftrnh/2s9XffQdAj9RUZ0SLFtC+fWyyRFltWZYuXVqiqsdEspy4FoW0tDQ+//zzkNPy8/PJzc2NZ5ywLEto8+fP56STTgo7vbS0lP79+8c8R6O2yfTp0DX8vVZy777bWccf/+iMKC6GceNikyXKassiIhsjXY4dPhhjfKwoGGN84nr4YEyiqz5saMLqLAoicjzO/e1Sca7Ym6GqT7m3tXoV5/r0ImCEqv5Y3wD79++nuLiYdu3asWpVpPfQiK1Ey7Jhwwa6du1KixZh74BmTNREsqfwE3Cbqha4F58sFZH3cS4++UBVJ4vIXcBdwIT6BiguLiY5OZmUlBSOOuqo+s4eE2VlZSQnJwcdA3A67/bt20dxcTHduwd+cuBhb8q8eQDcPmxYwEmCU2dRUNWqi0tQ1TIRWYVzo8zLcS6WAZiFcx/CeheFiooK0tLSKC+v702NmwYRISUlhe+/D+Y2ABPfWF49nHFgP3/9dCPX5XQLJEs8vL1kCWBFIWIikoZzB57PgFS3YIBzk4zUMPOMxbndN6mpqeTn5/umt2vXjvLyciorKykrK6tX+FhJtCzl5eVUVFQctO3iIf2IPdXDSQK9kysoKioK2baysjIuGcvLyxu+no4d4cCBsJN3uKf951e16dgRallXo7JEWbSyRFwU3DvbzAVuUdVS7/1CVFVFJORFFKo6A5gBkJ2drTU/R121ahXJyckJtcseKstFF13Eyy+/TPtaTmS57777OPfccxk8eHC915mfn8+UKVN4+23/bQqqsiQlJdG3b996L7ex/HsK6/myPInr0kPvKRwO5ym0d9/XuUe4H8yVlMDw4bHJEmXRyhJRUXDv8T8X+G9VfcMd/Z2IdFLVb0WkE3V/f8Ahqeq21/Pnz6+z7QMPPBCHRMbEVp3nKbi3EPszsEpVH/dMehMY4w6PAf4e/Xjx8fjjj5Oenk56ejpPPvkkGzdupEePHlx33XWkp6ezadMm0tLSKCkpAeDBBx+kR48enH322YwaNYopU6YAkJeXx+uvvw44Z2/+/ve/Jysri4yMDL766isAFi9eTE5ODn379uWss85i9erG3o3NRFPrli1p3bJl3Q0PY5HsKQzEuTvOchEpdMfdjXN78NdE5Hqc+/yNiEagULs/I0aMYPz48ezevZuLLrrooOl5eXnk5eVRUlLCVVdd5ZtW1zHW0qVLefHFF/nss89QVc444wyys7NZs2YNs2bN4swzz/S1X7JkCXPnzmXZsmXs37+frKws+vXrF3LZHTt2pKCggGnTpjFlyhReeOEFevbsyT//+U+aN2/OggULuPvuu5k7d27tG8XEzTuTJgUdIXCRfPrwMc6NOkMJ+X0Eh5KPP/6YYcOGVV9IcsUVV/DJJ5/QrVu3gwoCwMKFC7n88stJSkoiKSmJSy+9NOyyr7jiCgD69evHG284R107d+5kzJgxrFmzBhFhfy0X5xgThIQ7o7G2/+xt2rSpdXrHjh2j1hMcjSvfWrVybnjcrFkzfvrJuTnyvffey6BBg5g3bx5FRUUJ00llHA/Ong3AvSNHBpwkOE3+2odzzjmHv/3tb+zevZtdu3Yxb948zjor7LeYMXDgQN566y0qKiooLy8/6NOCuuzcuZMuXboAMHPmzMZENzHwwRdf8MEXXwQdI1BNvihkZWWRl5fHgAEDOOOMM7jhhhtq/dixf//+XHbZZZx++ukMHTqUjIwM2rWL/ObDd955JxMnTqRv377Vew/GJJKEO3wIwq233sqtt95a/bysrIwVK1b42nhP2Ln99tuZNGkSu3fv5txzz63uaPT+5/e2z87Orj6sycnJ4euvv66e9tBDDwFOB6sdSphEYEWhAcaOHcvKlSupqKhgzJgxZGVlBR3JmKixotAAL79sXx1wuEpJkLNqg2RFwRiPuRMnBh0hcE2+o9EY42dFwRiPibNmMXHWrKBjBMoOH4zx+NSuRbE9hXC8F0BF6k9/+hO9evXimmuuYe/evQwePJjMzExeffVVcnNzw97e3phEknB7Ct7r96Ph4Ssyorq82kybNo0FCxbQtWtXFi1aBEBhYSEA06dPj1sOYxrD9hRwvlnn4osvpk+fPqSnp1dftfj0008fdOnzpEmTqi+VBkhPT6eoqIhf/epXrF+/nqFDh/LII48wevRolixZQmZmJuvWrfOt77333iMnJ4esrCyGDx9ut6IzCcWKAvDuu+/SuXNnli1bxooVK6rvnFR16fO4ceN8hSCUZ599ls6dO/PRRx8xYcIEXnjhBc455xwKCwt936pUUlLCQw89xIIFCygoKCA7O5vHH3+8liWbeOqakkLXlJSgYwQq4Q4fgpCRkcFtt93GhAkTuOSSS8jMzARCX/rcWIsWLWLlypUMHDgQgH379pGTkxOVZZvG+6/bbgs6QuCsKACnnnoqBQUFzJ8/n9/97necffbZQOhLn5s3b84Bz40/Kyoq6rUuVeUXv/gFr7zySpTSGxNddvgAbNmyhTZt2jB69GjuuOMOli1bFrZtWloaBQUFABQUFLBhw4Z6revMM89k4cKFrF27FnD6M7wXSJlg3fL889zy/PNBxwiUFQVg+fLlDBgwgMzMTO6//37uuOOOsG2vvPJKfvjhB3r37s3UqVM59dRT67WuY445hpkzZzJq1ChOP/10cnJyqjsxTfAKN2ygsJ6F/nCTcIcP8fwIscqQIUMYMmRI9fOysrKwlz63bt2a9957L+RyvPPUvBTae0eo888/nyXul44Yk2hsT8EY42NFwRjjk3CHD8YE6dTOnYOOEDgrCsZ4zLj55qAjBM4OH4wxPlYUjPEYO3UqY6dODTpGoKwoRFl+fj6XXHIJAG+++SaTJ08O23bHjh1MmzYtXtFMBL7esoWvt2wJOkagEq9P4aaboru8556LymIqKytp1qxZvea57LLLuOyyy8JOryoK48ePb2w8Y6LG9hRwTjrq2bMn11xzDb169eLaa69l9+7dpKWlMWHCBLKyspgzZ07YS57fffddevbsSVZWlu/CqZkzZ3Kz23H13XffMWzYMPr06UOfPn345JNPuOuuu1i3bh2ZmZm1nkVpTDxZUXCtXr2a8ePHs2rVKpKTk6t361NSUigoKGDw4MEhL3muqKjgxhtv5K233mLp0qVs3bo15PJ/+9vfct5557Fs2TIKCgro3bs3kydP5qSTTqKwsJBHH300ni/XmLDqLAoi8hcR2SYiKzzjJonIZhEpdB8Hfz/8Ieb444+vvpz56quv5uOPP64eBv8lz5mZmcyaNYuNGzfy1Vdf0b17d0455RREhNGjR4dc/ocffsi4ceMA56rL+nzVnImfzO7dyezePegYgYqkT2EmMBX4a43xT6hq7XceOYSISMjnVd8+He6S56rbrZnDw5M33hh0hMDVuaegqv8AfohDlkB98803fPrppwDMmTOn+p4KVcJd8tyzZ0+Kioqqb7kW7j4JF1xwQfV9GisrK9m5cyfJycmUlZXF6iUZ0yCN6VO4WUS+cA8vOkQtUUB69OjBM888Q69evdixY0f1rn6VcJc8JyUlMWPGDC6++GKysrI49thjQy7/qaee4qOPPiIjI4N+/fqxcuVKUlJSGDhwIOnp6dUdjVV3fTLBGP3YY4x+7LGgYwRKVLXuRiJpwNuqmu4+TwVKAAUeBDqp6r+FmXcsMBYgNTW13+zZs33T27Vrx8knn9ygj/yiZePGjYwYMYLPPvsMaNjHj7FSlWXt2rXs3Lkz7uvfvGNP9XBr3UfFASGlbauQbSsrK2nTpk3MM5WXl9O2bduGzfz999CyZdjJt9xzDwBP/uEPzoh9++CYY2KTJcpqyzJo0KClqpodyXIadJ6Cqn5XNSwizwNv19J2BjADIDs7W2t+3XpVb39ZWRnJAX25Z9u2bTniiCOq1x9klpqqsiQlJdG3b9+4r997y/2MA+v5sjyJ69K7hWxbWlpK//79Y54pPz+fmu+jiE2fDl27hp3c3u1Lyj3C3YkuKYHhw2OTJcqilaVBhw8i0snzdBiwIlzbQ0FaWhorVhzSL8GYqKlzT0FEXgFygY4iUgz8HsgVkUycw4ciIMqnIRpjglJnUVDVUSFG/zmaISLp12jKbPvET06PHkFHCFzg1z4kJSWxfft2WtbS+dOUqSrbt28nKSkp6ChNwsNjxgQdIXCBF4WuXbtSXFzMjh07EuaNX1FRkVBZ2rdvT9daOseMiabAi0KLFi3o3r07+fn5gfSuh2JZmq4rH34YgLkTJwacJDiBFwVjEsl2O8PUrpI0xvhZUTDG+FhRMMb4WJ+CMR4XnH560BECZ0XBGI97R44MOkLg7PDBGONjRcEYj6GTJjF00qSgYwTKDh+M8dizb1/QEQJnewrGGB/bUzDmmWd+Ht68ObgcCcL2FIwxPranYIzHJbXcj7GpsKJgjMftTfyLYMAOH4wxNVhRMMYjd/FichcvDjpGoKwoGGN8rCgYY3ysKBhjfKwoGGN87CNJYzxGHHdc0BECZ0XBGI/xJ5wQdITAWVEwxmN3ZSUAsf/u7MRlRcEYj4uWLgUgP9gYgbKORmOMjxUFY4yPFQVjjE+dRUFE/iIi20RkhWfc0SLyvoiscX92iG1MY0y8RLKnMBO4sMa4u4APVPUU4AP3uTGHvLwuXcjr0iXoGIGq89MHVf2HiKTVGH05kOsOz8LprJ0QzWDGBKGpFwQAUdW6GzlF4W1VTXef71DV9u6wAD9WPQ8x71hgLEBqamq/2bNnh1xHeXk5bdu2rf8riAHL8rPNO/ZUD7fWfVQcEFLatgrZtrKykjZtYv8Jf6O2yfffQ8uW/nHbtlUP7nS/dbrdSSc5I/btg1ruxhT078ertiyDBg1aqqrZkSyn0ecpqKqKSNjKoqozgBkA2dnZmpubG7Jdfn4+4abFm2X52cQ3llcPZxxYz5flSVyX3i1k29LSUvr37x/zTI3aJtOnQ9eu/nFz51YPVt1LIX/WLGdESQkMHx6bLFEWrSwN/fThOxHpBOD+3FZHe2PMIaKhReFNYIw7PAb4e3TiGGOCFslHkq8AnwI9RKRYRK4HJgO/EJE1wGD3uTHmMBDJpw+jwky6IMpZjDEJwC6IMsZj3PHHBx0hcFYUjPG4ulOnoCMEzoqCMR6b9jjnZTTl/QUrCsZ4XLvcOS8jP9gYgbKrJI0xPlYUjDE+VhSMMT5WFIwxPtbRaJqOm26ClSvhyCPDNrktLS1+eRKUFQVjPC499tigIwTOioIxHqt37QKgR8A5gmRFwRiPm778ErDzFIwxppoVBWOMjxUFY4yPFQVjjI91NBrj8buquzg3YVYUjPEYnJISdITAWVEwxqOwtBSAzGBjBMqKgjEet3z1FWDnKRhjTDUrCsYYHysKxhgfKwrGGB/raDTG44+nnBJ0hMBZUTDG46wOHYKOEDgrCsZ4fPLjjwCcFXCOIFlRMMbj7jVrADtPwRhjqllRMMb4NOrwQUSKgDKgEvhJVbOjEcoYE5xo9CkMUtWSKCzHGJMArKPRGI8ne/YMOkLgRFUbPrPIBuBHQIHnVHVGiDZjgbEAqamp/WbPnh1yWeXl5bRt27bBWaLJsvxs84491cOtdR8VB4SUtq1Ctq2srKRNmzYxz9TgbfLNN7BnDxwRQVda1fc/7NsHxxwT/SwxUFuWQYMGLY308L6xRaGLqm4WkWOB94HfqOo/wrXPzs7Wzz//POS0/Px8cnNzG5wlmizLzya+sbx6OOPAer4sb811Od1Cti0tLaV///4xz9TgbRLBN0Qt2L4dgMH33eeMKC6GceOinyUGassiIhEXhUYdPqjqZvfnNhGZBwwAwhYFYxLdQ+vWATA44BxBavBHkiJypIgkVw0D/w9YEa1gxphgNGZPIRWYJyJVy3lZVd+NSipjTGAaXBRUdT3QJ4pZjDEJwM5oNMb42HkKxng817t30BECZ0XBGI8etXxc2VRYUTDG461t2wC4NOAcQbKiYIzHY0VFQNMuCtbRaIzxsaJgjPGxwwdTb3/9dGPI8b/sbTc9PRzYnoIxxsf2FIzxeCkjI+gIgbOiYIzH8a1bBx0hcFYUjPF49dtvAbj6mWecEbt2QWHhzw2eey7+oeLMioIxHtM3bQLg6k6dAk4SHOtoNMb4WFEwxvhYUTDG+FhRMMb4WEejMR6vZ2YGHSFwVhSM8ejYsmXQEQJnhw/GeMzcvJmZmzcHHSNQVhSM8bCiYEXBGFODFQVjjI8VBWOMjxUFY4yPfSRpoubvhVt4Y9Pyg8Y/fMWhc4+C+f36BR0hcFYUjPFo06xZ0BECZ0XBHH5uuqnBs0775hsAxp9wQrTSHHKsT8EYj9e2buW1rVuDjhGoRhUFEblQRFaLyFoRuStaoYwxwWlwURCRZsAzwFDgNGCUiJwWrWDGmGA0Zk9hALBWVder6j5gNnB5dGIZY4LSmI7GLsAmz/Ni4IzGxTGHo4lvHNofUzY1oqoNm1HkKuBCVb3BfX4tcIaq3lyj3VhgrPu0B7A6zCI7AiUNChN9liW0RMmSKDng0MnSTVWPiWQhjdlT2Awc73ne1R3no6ozgBl1LUxEPlfV7EbkiRrLElqiZEmUHHB4ZmlMn8IS4BQR6S4iLYGRwJuNDWSMCVaD9xRU9ScRuRn4H6AZ8BdV/TJqyYwxgWjUGY2qOh+YH6UsdR5ixJFlCS1RsiRKDjgMszS4o9EYc3iy05yNMT5xLQoicrSIvC8ia9yfHUK0GSQihZ5HhYj80p02U0Q2eKZlxjKL267Ss743PeO7i8hn7iner7qdrTHJISKZIvKpiHwpIl+IyNWeaY3eJnWdri4irdzXuNZ9zWmeaRPd8atFZEh9192ALLeKyEp3O3wgIt0800L+rmKYJU9Evves8wbPtDHu73SNiIyJQ5YnPDm+FpEdnmn12y6qGrcH8J/AXe7wXcAjdbQ/GvgBaOM+nwlcFc8sQHmY8a8BI93hZ4FxscoBnAqc4g53Br4F2kdjm+B0Eq8DTgRaAsuA02q0GQ886w6PBF51h09z27cCurvLaRbjLIM874dxVVlq+13FMEseMDXM+3a9+7ODO9whlllqtP8NTsd/g7ZLvA8fLgdmucOzgF/W0f4q4B1V3Z0AWaqJiADnA683ZP765lDVr1V1jTu8BdgGRHQiSgQiOV3dm/F14AJ3G1wOzFbVvaq6AVjrLi9mWVT1I8/7YRHO+TGx0JjT+IcA76vqD6r6I/A+cGEcs4wCXmnoyuJdFFJV9Vt3eCuQWkf7kRz84v7g7jo+ISKt4pAlSUQ+F5FFVYcxQAqwQ1V/cp8X45z2HcscAIjIAJz/Fus8oxuzTUKdrl7ztVS3cV/zTpxtEMm80c7idT3wjud5qN9VrLNc6W7710Wk6mS+wLaLezjVHfjQM7pe2yXqN1kRkQXAcSEm3eN9oqoqImE/+hCRTkAGznkQVSbi/OG0xPn4ZQLwQIyzdFPVzSJyIvChiCzH+aOIWJS3yUvAGFU94I6u1zY5XIjIaCAbOM8z+qDflaquC72EqHgLeEVV94rITTh7U+fHcH2RGAm8rqqVnnH12i5RLwqqOjjcNBH5TkQ6qeq37ht8Wy2LGgHMU9X9nmVX/UfdKyIvArfHOouqbnZ/rheRfKAvMBdoLyLN3f+cIU/xjmYOETkK+P/APaq6yLPsem2TECI5Xb2qTbGINAfaAdsjnDfaWRCRwTgF9TxV3Vs1PszvqqFFoc4sqrrd8/QFnP6hqnlza8yb38AcEWXxGAn8ukbO+m2XaHXMRNhh8ij+TrX/rKXtImBQjXGd3J8CPAlMjmUWnE6iVu5wR2ANbgcPMAd/R+P4GOZoCXwA3BJiWqO2Cc4/hvU4u5xVnVi9a7T5Nf6Oxtfc4d74OxrX07iOxkiyVL2hT4n0dxXDLJ08w8OARe7w0cAGN1MHd/joWGZx2/UEinDPP2rodonaH3yELy7FfXOvARZUbSic3cAXPO3ScCrhETXm/xBYDqwA/gtoG8sswFnu+pa5P6/3zH8isBinc21O1YaPUY7RwH6g0PPIjNY2AS4Cvnb/2O5xxz0AXOYOJ7mvca37mk/0zHuPO99qYGgU3iN1ZVkAfOfZDm/W9buKYZaHgS/ddX4E9PTM+2/u9loL/Guss7jPJ1Hjn0JDtoud0WiM8bEzGo0xPlYUjDE+VhSMMT5WFIwxPlYUjDE+VhSMMT5WFIwxPlYUjDE+/we40Zr5OBRRaAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(4, 3))\n",
    "\n",
    "coeffs_shuf = corr_glm_fit_shuffled[corr_glm_fit_shuffled[:, 1] < 0.95][:, 0]\n",
    "coeffs_chun = corr_glm_fit_train_test[corr_glm_fit_train_test[:, 1] < 0.95][:, 0]\n",
    "bins = np.linspace(-1, 1, 50)\n",
    "\n",
    "ax.hist(coeffs_shuf, bins=bins, alpha=0.6, density=True, label='shuffle', color='tab:blue')\n",
    "ax.hist(coeffs_chun, bins=bins, alpha=0.6, density=True, label='predict.', color='red')\n",
    "ax.axvspan(np.percentile(coeffs_shuf, 5), np.percentile(coeffs_shuf, 95), alpha=0.3, color='gray')\n",
    "ax.axvspan(np.percentile(coeffs_chun, 5), np.percentile(coeffs_chun, 95), alpha=0.3, color='red')\n",
    "ax.axvline(corr_glm_fit_orig[0], color='black', ls='--', label='original')\n",
    "ax.set_xlim(-0.8, 0.8)\n",
    "ax.grid()\n",
    "ax.legend(loc='upper left')\n",
    "ax.set_title('Unit %s (%.2f)' % (unit, corr_glm_fit_orig[0]), fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df33355f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}