2023-06-30_dataset2_laurent_decoding.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f1fec7f6-4c47-4c18-90ec-11e430dab517",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9008a877",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 2267760\n",
      "-rw-r--r--  1 laurentperrinet  staff   211M Jun 14 10:42 H_Ready-Set_bin1ms_alignedReady.mat\n",
      "-rw-r--r--  1 laurentperrinet  staff   407M Jun 14 10:55 G_Ready-Set_bin1ms_alignedReady.mat\n",
      "-rw-r--r--@ 1 laurentperrinet  staff   2.0K Jun 15 12:29 README.md\n",
      "-rw-r--r--@ 1 laurentperrinet  staff   183M Jun 27 15:52 G_Ready-Set_bin1ms_alignedReady_session=0.npz\n",
      "-rw-r--r--@ 1 laurentperrinet  staff   212M Jun 27 15:52 G_Ready-Set_bin1ms_alignedReady_session=1.npz\n",
      "-rw-r--r--@ 1 laurentperrinet  staff    17M Jun 27 15:52 H_Ready-Set_bin1ms_alignedReady_session=0.npz\n",
      "-rw-r--r--@ 1 laurentperrinet  staff    68M Jun 27 15:52 H_Ready-Set_bin1ms_alignedReady_session=1.npz\n"
     ]
    }
   ],
   "source": [
    "%ls -ltrh /Users/laurentperrinet/quantic/grants/Polychronies/RSG-spikepatterns/Data/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "16476f59",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1200, 107, 1498)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_session = 0\n",
    "monkey = 'G'\n",
    "from pathlib import Path\n",
    "path = Path.home()/'quantic'/'grants'/'polychronies'/'RSG-spikepatterns'/'Data'/f'{monkey}_Ready-Set_bin1ms_alignedReady_session={i_session}.npz'\n",
    "\n",
    "dictdata = np.load(path)\n",
    "data, id_neuron, id_short, t_s, t_p = dictdata['data'], dictdata['id_neuron'], dictdata['id_short'], dictdata['t_s'], dictdata['t_p']\n",
    "data = data.astype(float)\n",
    "data[data==255] = np.nan\n",
    "N_time, N_neuron, N_trial = data.shape\n",
    "\n",
    "N_time, N_neuron, N_trial"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "62b0b690",
   "metadata": {},
   "source": [
    "### computing the firing rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "aa340cf9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 100,  140,  180,  220,  260,  300,  340,  380,  420,  460,  500,\n",
       "         540,  580,  620,  660,  700,  740,  780,  820,  860,  900,  940,\n",
       "         980, 1020, 1060, 1100, 1140, 1180]),\n",
       " 28)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_s_max = int(t_s.max())\n",
    "DeltaT = 100 # number of time bins on which we average\n",
    "dt = 40 # temporal distance between points on which we compute the firing rate\n",
    "\n",
    "timepoints = np.arange(DeltaT, t_s_max+1, dt)\n",
    "N_timepoint = len(timepoints)\n",
    "timepoints, N_timepoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5ec9a4e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fr = np.empty((N_timepoint, N_neuron, N_trial)) * np.nan\n",
    "for i_timepoint, timepoint in enumerate(timepoints):\n",
    "    fr[i_timepoint, :, :] = np.nansum(data[(timepoint-DeltaT):timepoint, :, :], axis=0) / (DeltaT/1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a102fd8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i_trial in range(N_trial):\n",
    "#     print(i_trial, torch.isnan(torch.tensor(fr[:, :, i_trial].astype(float))).sum())\n",
    "#     break\n",
    "# torch.isnan(torch.tensor(fr.astype(float))).sum(), torch.isnan(torch.tensor(data.astype(float))).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cb35aebf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(28, 107, 1498)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f825ab78",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([783, 418,  15, ..., 539, 858, 573])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratio_train = 80/100\n",
    "N_trial_train = int(N_trial*ratio_train)\n",
    "N_trial_test = N_trial - N_trial_train\n",
    "idx_trial = np.random.permutation(N_trial)\n",
    "idx_trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5b4b0c83",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x122acb510>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAAE0CAYAAACLqfDEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnL0lEQVR4nO3dfZCV5Xk4/mvZlYXV3SVg2GUVAlpBiy9NJRLUWlMZkGac+DJptaSDxq9O28UKtE0kjabmjWqm0dESqZ1W2/mGmDgTtNrEhqLi1xGIklLrxixq8QdRFxvSfQHrirvn90emG1fBcy/7POec3f18Zs7onnPtfV/nOddzP8+5ePacqkKhUAgAAACAHIwrdwIAAADA6KXxAAAAAORG4wEAAADIjcYDAAAAkBuNBwAAACA3Gg8AAABAbjQeAAAAgNxoPAAAAAC50XgAAAAAclNT7gTerb+/P1599dWor6+PqqqqcqcDAAAAHEKhUIienp5oaWmJceMOf11Dbo2HtWvXxte+9rXo6OiIM844I+68884466yziv7eq6++GtOnT88rLQAAACBDe/bsieOPP/6wj+fSePj2t78dq1atinXr1sX8+fPj9ttvj8WLF0d7e3tMnTr1fX+3vr4+IiLOjd+OmjjqsHHVp5yUlEvf8y+kJ56BLPNKGavUzy/VSM49S29cNK9oTN1Dz2Q2X6XuF1S21LopZiTX1X/+5UeS4k644enM5iz1OpnlfJWYe6pKrdO9rfOLxrQ8+vOiMeXY7pW6TUttJJ/7VGLuzZuOSYrb/lrxf7A87lPPDzcdjlCW63epVer+mqKU70HejoPxZHxv4H384eTSePj6178e11xzTVx11VUREbFu3br453/+5/j7v//7uOGGG973d//3zytq4qioqXqfxkN1bVIuVe8zRh6yzCtlrFI/v1QjOfcs1Rw1oXhMhtuhUvcLKltq3RQzkutq3ITi+2pE6ffXLLdplvNVYu6pKrVOq2sTjhcVut0rdZuW2kg+96nE3McfMz4prrqueO5Zrt0MTZbrd6lV6v6aoqTvQQq/+E+xj0nI/MMl33rrrdi+fXssXLjwl5OMGxcLFy6MLVu2vCe+t7c3uru7B90AAACA0SHzxsPPfvaz6Ovri6ampkH3NzU1RUdHx3vi16xZE42NjQM3n+8AAAAAo0fZv05z9erV0dXVNXDbs2dPuVMCAAAAMpL5Zzwce+yxUV1dHXv37h10/969e6O5ufk98bW1tVFbO3L/9gcAAAA4vMyveBg/fnyceeaZsWnTpoH7+vv7Y9OmTbFgwYKspwMAAAAqWC7farFq1apYtmxZzJs3L84666y4/fbb48CBAwPfcpGFvrb2zMaqVCnPsXrunMzGylLP7ElFY+ra8s9jNEl5rcfCfkH2UurmjUuKf9XfSN6nJ8zoKfmcI3mdfGXRlKIxzRnmnrKtIiLqd3YWjanUtbTh5b6Sz1mMY8rQZLWWRkTUbdg23HRGvCe2zE2Ka3miUDQm9Xw5hf2iPFKPAyn7Tpb1UKlSjoelPurk0nj43d/93fiv//qvuOmmm6KjoyN+7dd+LR555JH3fOAkAAAAMLrl0niIiFi+fHksX748r+EBAACAEaDs32oBAAAAjF4aDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMhNTbkTYHTqnlldNKauBHmUW92GbSWdr3runKS4vrb2nDNhJEmtm1KNE5FeoylzpoxV82Rj0nyjXZavYZbqd3aWO4XcpTzHVxZNKRrT3JZBMlABWp4oJMUlnWslrm09sycVn88+NiQpx+CUY0/qcaAvKWr0q8RadsUDAAAAkBuNBwAAACA3Gg8AAABAbjQeAAAAgNxoPAAAAAC50XgAAAAAcqPxAAAAAORG4wEAAADITU25Exht+trak+Kq587JbCzGhpR6SKkrRoc3LplfNKZuw7bM5ksaK7H+emZPKj5fW9JQmTnuB/vSAjNcu+t3dhYfK2mkNFmuIQ0vZ5lZdkbycTMl94bZxfd7KlvKfh+R7b6fNF8F7jupx7BMz6krcB9LXZcr8TWMSDtfSToeZvj8KnVbZanU5xgpXPEAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5Kam3Akcqeq5c5Li+trac87kyLyyaErRmOMi7TlWooaX+8qdwpiUUlcREc1tOSdC7rpnVheNqStBHqNF6rEi9diT5ZxZSck9OafZ84eZzS9lmlcFGunnK6QbybVc6tyznK9n9qSiMXWJ5z31OzuLxjjDrWwjeT/MUiU+R1c8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3NSUO4GxquHlvnKnkKv6nZ1FY0b3Fshe9dw55U6BClLqNSSl/vra2pPGqo+EWk6s99Q5yVb3zOqiMXUlyAMgy3POSjymVGJOQ1G3YVvRmJ5L5hcfpy2LbH7hjYT5ItJyr1RZnrdlJfMrHv7iL/4iqqqqBt1OPvnkrKcBAAAARoBcrniYO3du/Ou//usvJ6lxYQUAAACMRbl0BGpqaqK5uTmPoQEAAIARJJcPl3zhhReipaUlTjjhhFi6dGns3r37sLG9vb3R3d096AYAAACMDpk3HubPnx/33ntvPPLII3HXXXfFrl274jd+4zeip6fnkPFr1qyJxsbGgdv06dOzTgkAAAAok8wbD0uWLIlPfvKTcfrpp8fixYvje9/7XnR2dsZ3vvOdQ8avXr06urq6Bm579uzJOiUAAACgTHL/1MdJkybF7Nmz48UXXzzk47W1tVFbW5t3GgAAAEAZ5PIZD++0f//+eOmll2LatGl5TwUAAABUmMwbD3/6p38amzdvjpdffjmeeuqpuOSSS6K6ujquuOKKrKcCAAAAKlzmf2rx05/+NK644orYt29ffPCDH4xzzz03tm7dGh/84AeznmpEq9/ZWe4UGKbquXMyGaevrT0p7v/7YvHddfL/7RtuOowQKWtIajX0zJ40nFQG1LWlxaXUfFb7V9ZS99espG6HUueVIjX3lPpLra2RrFJrHvJQiWsW5ZPlOU2WUtZltZwu88bDfffdl/WQAAAAwAiV+2c8AAAAAGOXxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByU1PuBI5UX1t7uVMYlpGeP6V31nG7i8bs3dmSNFbfcJOh7Eq9hnTPrC4aU1eCPDhyPbMnFY2pa0sb67gf7BteMkNUPXdOUlwlHltTc0p5jnUbtg03nSEZydu9UlXqNvUapqvU15BfSDnW1e/szGys1ONmqVVi/bniAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5Kam3AlAMdVz5xSN6WtrL/lYKV5ZNKVozHFRPKeIiB99u/hYzW1PJY0F71S3YVvRmO6VZ5cgk1/Kcj9MkbI2RJQ+ryx1z6wuGlOXOFam22H2/OzGGsF6Zk8qGlPXln8e5CvldY4o/Wtd6vMjfmEsbPeU55jl+pd0rNuQuE1H8PGplO+fCn29Ec8XH8cVDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuasqdAIdXPXdO0Zi+tvYSZDJ69MyeVDSmri1trJRtv///fDRtsAQNL/dlNhbkIWXNirBujQZZHp/qNmwrHpRYWykq9dhav7OzeFCF5k667pnVSXF1OecxEqQeU7I8t0thHxualNenUqWsy5V6dl7q/SKFKx4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuasqdQN6q584pGtPX1l6CTAZLySvLccrxHIk4Znd2vb36nZ1FY/oym41Kl+XaltV6RGVreLm0K0Q5jk8pc/bMnlQ0pq4tg2SGKCWvFFnm7twhe6XeD8eC7pnVRWPqSpAHRy7Lc9zm254aVi6D5rQGRkTx49PbB9+MeL74OEN+V/TEE0/ERRddFC0tLVFVVRUPPPDAoMcLhULcdNNNMW3atJg4cWIsXLgwXnjhhaFOAwAAAIwCQ248HDhwIM4444xYu3btIR+/9dZb44477oh169bFtm3b4uijj47FixfHm2++OexkAQAAgJFlyH9qsWTJkliyZMkhHysUCnH77bfH5z//+fjEJz4RERH/+I//GE1NTfHAAw/E5ZdfPrxsAQAAgBEl0w+X3LVrV3R0dMTChQsH7mtsbIz58+fHli1bDvk7vb290d3dPegGAAAAjA6ZNh46OjoiIqKpqWnQ/U1NTQOPvduaNWuisbFx4DZ9+vQsUwIAAADKqOxfp7l69ero6uoauO3Zs6fcKQEAAAAZybTx0NzcHBERe/fuHXT/3r17Bx57t9ra2mhoaBh0AwAAAEaHTBsPs2bNiubm5ti0adPAfd3d3bFt27ZYsGBBllMBAAAAI8CQv9Vi//798eKLLw78vGvXrtixY0dMnjw5ZsyYEStWrIgvf/nLcdJJJ8WsWbPixhtvjJaWlrj44ouzzBsAAAAYAYbceHjmmWfiYx/72MDPq1atioiIZcuWxb333huf+cxn4sCBA3HttddGZ2dnnHvuufHII4/EhAkTssuaXFTPnZPZWH1t7RU5VvfM6qIxdZnNFtHwcl/RmCyfH+QlpZbJXpbrQ/3OzuLzJY6VklfqMSXLY89IlvL69MyelHseHLmU/aI+0uq91CvuSD4X2T+jv6Tzpa5ZKdt0JG/3VClrW4ost3uWY40FxV7Dt/t6k8YZcuPh/PPPj0KhcNjHq6qq4otf/GJ88YtfHOrQAAAAwChT9m+1AAAAAEYvjQcAAAAgNxoPAAAAQG40HgAAAIDcaDwAAAAAudF4AAAAAHKj8QAAAADkpqbcCXB4fW3t5U7hiL1xyfyiMXUbtpUgE8hf9dw5SXFZ7tMlH2t28X16LEh5rcuxdqfW4Eidr2f2pKS4lONK96Kzi4+TNFu2Up8jjAbZrqUfHV4y71DqtS1FOc4xSi1l/avf2Zl7HqNJ98zqojGpx7pir8/bB9+MeL74OK54AAAAAHKj8QAAAADkRuMBAAAAyI3GAwAAAJAbjQcAAAAgNxoPAAAAQG40HgAAAIDcaDwAAAAAudF4AAAAAHJTU+4EjlT13DnlTmFYemZPKhpT15Z/HnnpnlldNKauBHm82/4Z/WWYFUa++p2dRWP62tozmy91jc9yzhSlXruzPNZlua2StsOGbZnNVx+J9ZAQc9wP9hUPKkP9pexjKVK2AeVT6jVrLDhmt39HHQtSjjsREd2Lzi4a03zbU8PMZmwpdnx6u683aRx7KgAAAJAbjQcAAAAgNxoPAAAAQG40HgAAAIDcaDwAAAAAudF4AAAAAHKj8QAAAADkRuMBAAAAyE1NuRM4Uj2zJyXF1W3Ylm8iR6h+Z2fxoLlziob0tbUPP5kcxmqYPT+zsfiF6hLXA4wlyfvOCF7bKnUNSckrSynnD0nHaHiXLGvZ8bz0a0NE2nbvWHl20ljNbcPN5hdGei0k5Z9wbE1dl7tnTikak1pbI3nbN7zcV+4U3sMVDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDc15U7gSHXPrE6Kq8s5Dw6tbsO2cqdwSC1PFEo6X/3OzqIxfYljvbJoStGY5rbEwSAH1XPnJMX1tbVnElMOlbq2pWyv1NenEmVZD0nrchnqrxJfwyz3aXinLOu94eXUMynGgubbnioeNIKPh+VQbH/tKxxMGmfIVzw88cQTcdFFF0VLS0tUVVXFAw88MOjxK6+8MqqqqgbdLrzwwqFOAwAAAIwCQ248HDhwIM4444xYu3btYWMuvPDCeO211wZu3/rWt4aVJAAAADAyDflPLZYsWRJLlix535ja2tpobm4+4qQAAACA0SGXD5d8/PHHY+rUqTFnzpz4wz/8w9i3b18e0wAAAAAVLvMPl7zwwgvj0ksvjVmzZsVLL70Un/vc52LJkiWxZcuWqK5+7wdC9vb2Rm9v78DP3d3dWacEAAAAlEnmjYfLL7984P9PO+20OP300+PEE0+Mxx9/PC644IL3xK9ZsyZuvvnmrNMAAAAAKkAuf2rxTieccEIce+yx8eKLLx7y8dWrV0dXV9fAbc+ePXmnBAAAAJRI5lc8vNtPf/rT2LdvX0ybNu2Qj9fW1kZtbW3eaQAAAABlMOTGw/79+wddvbBr167YsWNHTJ48OSZPnhw333xzXHbZZdHc3BwvvfRSfOYzn4lf+ZVficWLF2eaOAAAAFD5htx4eOaZZ+JjH/vYwM+rVq2KiIhly5bFXXfdFc8++2z8wz/8Q3R2dkZLS0ssWrQovvSlL7mq4V16Zk8qGlO/s7NoTPXcOUnz9bW1J8WNdnUbtpV0viy3+/4Z/ZmNxciXsu+n1l/qOkLpWbt5p5Rzh7q20s6X9ZwwVCnny32JY6UcD4/7Qdq39aXOSZosj4eOreUx5MbD+eefH4VC4bCP/8u//MuwEgIAAABGj9w/XBIAAAAYuzQeAAAAgNxoPAAAAAC50XgAAAAAcqPxAAAAAORG4wEAAADIjcYDAAAAkJuacidwpBpe7kuKq547p2hMX1v7cNMZsvqdnSWfk2zrodS1dcxufUJ+qWf2pKIxdW355/FOKTlFlD6vkeyNS+YnxdVt2JZzJoOlHMPSjtJp62TKeps6VoqOlWcnxTXf9lQm80Wkv9YwVqQeU5xTj3yvnldVNOZXNmQ3X6mPKeWQ5XE6K97JAAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALmpKXcCY1XP7ElFY+p3duaeR16q584pGtPX1l6CTAZL2e51baUfK0XzbU9lNxiZKkctJ60PCfthqpTnWB9p8/UNN5khqtT1KEX3zOqkuPqE55jlmvXKoilFY45LrYcSb/uU+X797+qTxnr1tuFmMzR1G7aVdkIoo9Tz4CzXkKR10n7ICFDq9ykpXPEAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByU1PuBI5U98zqpLj6nTkncoTqd3YWjemZPaloTN2GbcNPZoiq584pGpOUe1sGyZRRymvYl38ajFGvLJpSNOa4H+xLGquvrb1oTMp+n5JTRETzCN/3S6nh5bRVJHXbF1OXyShD88Yl84vGpKy3Wbpnxv9Lilscv5bZnEnPMWE/TNmfKR+vT7rUde24yG6/KPW5XcqxdSzUzDG7S/tv4WNhm1YiVzwAAAAAudF4AAAAAHKj8QAAAADkRuMBAAAAyI3GAwAAAJAbjQcAAAAgNxoPAAAAQG40HgAAAIDcaDwAAAAAuakZSvCaNWviu9/9bvzkJz+JiRMnxtlnnx233HJLzJkzZyDmzTffjD/5kz+J++67L3p7e2Px4sXxjW98I5qamjJPPkVfW3tJ56ueO6d4UES8smhK0ZjjfrCvaExf0myl9/NP7S8aU7ehBIm8S/3OzqIxqdu0Z/akojF1bYmDwRA1vFzavT+l3sleypoVEdE9s/gxJcua2T+jv2hM8vF39vzMxko5BqeMddXu30iaL6InMa64lLzeuKT4tsryuNM9szopri67KWFAyjqTtVK/bxgLUt8blVLKWhoRUbdhW86ZjC1DuuJh8+bN0draGlu3bo2NGzfGwYMHY9GiRXHgwIGBmJUrV8ZDDz0U999/f2zevDleffXVuPTSSzNPHAAAAKh8Q7ri4ZFHHhn087333htTp06N7du3x3nnnRddXV3xd3/3d7F+/fr4rd/6rYiIuOeee+KUU06JrVu3xkc/+tHsMgcAAAAq3rA+46GrqysiIiZPnhwREdu3b4+DBw/GwoULB2JOPvnkmDFjRmzZsmU4UwEAAAAj0JCueHin/v7+WLFiRZxzzjlx6qmnRkRER0dHjB8/PiZNmjQotqmpKTo6Og45Tm9vb/T29g783N3dfaQpAQAAABXmiK94aG1tjeeeey7uu+++YSWwZs2aaGxsHLhNnz59WOMBAAAAleOIGg/Lly+Phx9+OB577LE4/vjjB+5vbm6Ot956Kzo7OwfF7927N5qbmw851urVq6Orq2vgtmfPniNJCQAAAKhAQ2o8FAqFWL58eWzYsCEeffTRmDVr1qDHzzzzzDjqqKNi06ZNA/e1t7fH7t27Y8GCBYccs7a2NhoaGgbdAAAAgNFhSJ/x0NraGuvXr48HH3ww6uvrBz63obGxMSZOnBiNjY1x9dVXx6pVq2Ly5MnR0NAQ1113XSxYsMA3WgAAAMAYNKTGw1133RUREeeff/6g+++555648sorIyLitttui3HjxsVll10Wvb29sXjx4vjGN76RSbLvdNwP9iXF9WU+czYaXi6eWV9bewkyGbqUvGqePLsEmZRX/c7OojGVWn+MDaVeQyp1Xa7UtTRFau4Ns+cXjanbsG246Qw4ZvewvhRrkCzzysqLt/xqUlxdVF7uWUo5V4G8pK4zY2GNH+1Szh+sRkNTie9ThtR4KBQKRWMmTJgQa9eujbVr1x5xUgAAAMDokN0/WQAAAAC8i8YDAAAAkBuNBwAAACA3Gg8AAABAbjQeAAAAgNxoPAAAAAC50XgAAAAAclNT7gSOVF9be7lTOKTUvOpjTvGxhptMGR33g31FY8rx/Cq1bmCo6jZsK+l89Ts7SzofQ1Pqemi+7amSzldqpd6eqSp1vx/J5ytZco6RrYaXVdZokLJfVM8t/r4oS2NhbavE9cgVDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuqgqFQqHcSbxTd3d3NDY2xvnxiaipOqrc6QAAAACH8HbhYDweD0ZXV1c0NDQcNs4VDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Ayp8bBmzZr4yEc+EvX19TF16tS4+OKLo729fVDM+eefH1VVVYNuf/AHf5Bp0gAAAMDIMKTGw+bNm6O1tTW2bt0aGzdujIMHD8aiRYviwIEDg+KuueaaeO211wZut956a6ZJAwAAACNDzVCCH3nkkUE/33vvvTF16tTYvn17nHfeeQP319XVRXNzczYZAgAAACPWsD7joaurKyIiJk+ePOj+b37zm3HsscfGqaeeGqtXr4433nhjONMAAAAAI9SQrnh4p/7+/lixYkWcc845ceqppw7c/3u/93vxoQ99KFpaWuLZZ5+Nz372s9He3h7f/e53DzlOb29v9Pb2Dvzc3d19pCkBAAAAFeaIGw+tra3x3HPPxZNPPjno/muvvXbg/0877bSYNm1aXHDBBfHSSy/FiSee+J5x1qxZEzfffPORpgEAAABUsCP6U4vly5fHww8/HI899lgcf/zx7xs7f/78iIh48cUXD/n46tWro6ura+C2Z8+eI0kJAAAAqEBDuuKhUCjEddddFxs2bIjHH388Zs2aVfR3duzYERER06ZNO+TjtbW1UVtbO5Q0AAAAgBFiSI2H1tbWWL9+fTz44INRX18fHR0dERHR2NgYEydOjJdeeinWr18fv/3bvx1TpkyJZ599NlauXBnnnXdenH766bk8AQAAAKByVRUKhUJycFXVIe+/55574sorr4w9e/bEpz71qXjuuefiwIEDMX369Ljkkkvi85//fDQ0NCTN0d3dHY2NjXF+fCJqqo5KTQ0AAAAoobcLB+PxeDC6urre9z3/kP/U4v1Mnz49Nm/ePJQhAQAAgFHsiD5cEgAAACCFxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5EbjAQAAAMiNxgMAAACQG40HAAAAIDcaDwAAAEBuNB4AAACA3Gg8AAAAALnReAAAAAByo/EAAAAA5Kam3Am8W6FQiIiIt+NgRKHMyQAAAACH9HYcjIhfvo8/nIprPPT09ERExJPxvTJnAgAAABTT09MTjY2Nh328qlCsNVFi/f398eqrr0Z9fX1UVVVFRER3d3dMnz499uzZEw0NDWXOEPKl3hlr1DxjiXpnLFHvjDVjseYLhUL09PRES0tLjBt3+E9yqLgrHsaNGxfHH3/8IR9raGgYMy8gqHfGGjXPWKLeGUvUO2PNWKv597vS4X/5cEkAAAAgNxoPAAAAQG5GROOhtrY2vvCFL0RtbW25U4HcqXfGGjXPWKLeGUvUO2ONmj+8ivtwSQAAAGD0GBFXPAAAAAAjk8YDAAAAkBuNBwAAACA3Gg8AAABAbkZE42Ht2rUxc+bMmDBhQsyfPz9++MMfljslGLY1a9bERz7ykaivr4+pU6fGxRdfHO3t7YNi3nzzzWhtbY0pU6bEMcccE5dddlns3bu3TBlDdv7yL/8yqqqqYsWKFQP3qXdGk1deeSU+9alPxZQpU2LixIlx2mmnxTPPPDPweKFQiJtuuimmTZsWEydOjIULF8YLL7xQxozhyPX19cWNN94Ys2bNiokTJ8aJJ54YX/rSl+Kdn2Gv5hmpnnjiibjooouipaUlqqqq4oEHHhj0eEpt//znP4+lS5dGQ0NDTJo0Ka6++urYv39/CZ9F+VV84+Hb3/52rFq1Kr7whS/Ej370ozjjjDNi8eLF8frrr5c7NRiWzZs3R2tra2zdujU2btwYBw8ejEWLFsWBAwcGYlauXBkPPfRQ3H///bF58+Z49dVX49JLLy1j1jB8Tz/9dPzN3/xNnH766YPuV++MFv/93/8d55xzThx11FHx/e9/P3784x/HX/3VX8UHPvCBgZhbb7017rjjjli3bl1s27Ytjj766Fi8eHG8+eabZcwcjswtt9wSd911V/z1X/91PP/883HLLbfErbfeGnfeeedAjJpnpDpw4ECcccYZsXbt2kM+nlLbS5cujba2tti4cWM8/PDD8cQTT8S1115bqqdQGQoV7qyzziq0trYO/NzX11doaWkprFmzpoxZQfZef/31QkQUNm/eXCgUCoXOzs7CUUcdVbj//vsHYp5//vlCRBS2bNlSrjRhWHp6egonnXRSYePGjYXf/M3fLFx//fWFQkG9M7p89rOfLZx77rmHfby/v7/Q3Nxc+NrXvjZwX2dnZ6G2trbwrW99qxQpQqY+/vGPFz796U8Puu/SSy8tLF26tFAoqHlGj4gobNiwYeDnlNr+8Y9/XIiIwtNPPz0Q8/3vf79QVVVVeOWVV0qWe7lV9BUPb731Vmzfvj0WLlw4cN+4ceNi4cKFsWXLljJmBtnr6uqKiIjJkydHRMT27dvj4MGDg+r/5JNPjhkzZqh/RqzW1tb4+Mc/PqiuI9Q7o8s//dM/xbx58+KTn/xkTJ06NT784Q/H3/7t3w48vmvXrujo6BhU742NjTF//nz1zoh09tlnx6ZNm2Lnzp0REfHv//7v8eSTT8aSJUsiQs0zeqXU9pYtW2LSpEkxb968gZiFCxfGuHHjYtu2bSXPuVxqyp3A+/nZz34WfX190dTUNOj+pqam+MlPflKmrCB7/f39sWLFijjnnHPi1FNPjYiIjo6OGD9+fEyaNGlQbFNTU3R0dJQhSxie++67L370ox/F008//Z7H1DujyX/+53/GXXfdFatWrYrPfe5z8fTTT8cf//Efx/jx42PZsmUDNX2o8xv1zkh0ww03RHd3d5x88slRXV0dfX198ZWvfCWWLl0aEaHmGbVSarujoyOmTp066PGampqYPHnymKr/im48wFjR2toazz33XDz55JPlTgVysWfPnrj++utj48aNMWHChHKnA7nq7++PefPmxVe/+tWIiPjwhz8czz33XKxbty6WLVtW5uwge9/5znfim9/8Zqxfvz7mzp0bO3bsiBUrVkRLS4uaByKiwj9c8thjj43q6ur3fKr53r17o7m5uUxZQbaWL18eDz/8cDz22GNx/PHHD9zf3Nwcb731VnR2dg6KV/+MRNu3b4/XX389fv3Xfz1qamqipqYmNm/eHHfccUfU1NREU1OTemfUmDZtWvzqr/7qoPtOOeWU2L17d0TEQE07v2G0+LM/+7O44YYb4vLLL4/TTjstfv/3fz9WrlwZa9asiQg1z+iVUtvNzc3v+WKEt99+O37+85+Pqfqv6MbD+PHj48wzz4xNmzYN3Nff3x+bNm2KBQsWlDEzGL5CoRDLly+PDRs2xKOPPhqzZs0a9PiZZ54ZRx111KD6b29vj927d6t/RpwLLrgg/uM//iN27NgxcJs3b14sXbp04P/VO6PFOeec856vR965c2d86EMfioiIWbNmRXNz86B67+7ujm3btql3RqQ33ngjxo0b/Laiuro6+vv7I0LNM3ql1PaCBQuis7Mztm/fPhDz6KOPRn9/f8yfP7/kOZdLxf+pxapVq2LZsmUxb968OOuss+L222+PAwcOxFVXXVXu1GBYWltbY/369fHggw9GfX39wN94NTY2xsSJE6OxsTGuvvrqWLVqVUyePDkaGhriuuuuiwULFsRHP/rRMmcPQ1NfXz/w+SX/6+ijj44pU6YM3K/eGS1WrlwZZ599dnz1q1+N3/md34kf/vCHcffdd8fdd98dERFVVVWxYsWK+PKXvxwnnXRSzJo1K2688cZoaWmJiy++uLzJwxG46KKL4itf+UrMmDEj5s6dG//2b/8WX//61+PTn/50RKh5Rrb9+/fHiy++OPDzrl27YseOHTF58uSYMWNG0do+5ZRT4sILL4xrrrkm1q1bFwcPHozly5fH5ZdfHi0tLWV6VmVQ7q/VSHHnnXcWZsyYURg/fnzhrLPOKmzdurXcKcGwRcQhb/fcc89AzP/8z/8U/uiP/qjwgQ98oFBXV1e45JJLCq+99lr5koYMvfPrNAsF9c7o8tBDDxVOPfXUQm1tbeHkk08u3H333YMe7+/vL9x4442FpqamQm1tbeGCCy4otLe3lylbGJ7u7u7C9ddfX5gxY0ZhwoQJhRNOOKHw53/+54Xe3t6BGDXPSPXYY48d8px92bJlhUIhrbb37dtXuOKKKwrHHHNMoaGhoXDVVVcVenp6yvBsyqeqUCgUytTzAAAAAEa5iv6MBwAAAGBk03gAAAAAcqPxAAAAAORG4wEAAADIjcYDAAAAkBuNBwAAACA3Gg8AAABAbjQeAAAAgNxoPAAAAAC50XgAAAAAcqPxAAAAAORG4wEAAADIzf8PLndzT5r89qsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1300x1300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 13))\n",
    "ax.imshow(fr[:, :, 102], cmap='viridis')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e9714952",
   "metadata": {},
   "source": [
    "### logistic regression model with NaNs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4c568a01",
   "metadata": {},
   "source": [
    "### inference model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6abe4441",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "if torch.backends.mps.is_available():\n",
    "    device = torch.device('mps')\n",
    "elif torch.cuda.is_available():\n",
    "    device = torch.device('cuda')\n",
    "else:\n",
    "    device = 'cpu'\n",
    "\n",
    "device = 'cpu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ff55e8fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/homebrew/lib/python3.11/site-packages/torch/__init__.py:615: UserWarning: torch.set_default_tensor_type() is deprecated as of PyTorch 2.1, please use torch.set_default_dtype() and torch.set_default_device() as alternatives. (Triggered internally at /Users/runner/work/pytorch/pytorch/pytorch/torch/csrc/tensor/python_tensor.cpp:453.)\n",
      "  _C._set_default_tensor_type(t)\n"
     ]
    }
   ],
   "source": [
    "import torch.nn.functional as F\n",
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "torch.set_default_tensor_type(\"torch.DoubleTensor\") # -> torch.tensor([1.2, 3]).dtype = torch.float64\n",
    "# torch.set_default_tensor_type(\"torch.FloatTensor\") # -> torch.tensor([1.2, 3]).dtype = torch.float\n",
    "# see https://sebastianraschka.com/faq/docs/pytorch-crossentropy.html#pytorch-loss-input-confusion-cheatsheet\n",
    "criterion = torch.nn.BCELoss(reduction=\"mean\") # loss divided by output size\n",
    "#criterion = torch.nn.NLLLoss(reduction=\"mean\") # loss divided by output size\n",
    "\n",
    "\n",
    "class LogisticRegressionModel(torch.nn.Module):\n",
    "    def __init__(self, N_neuron, N_timepoint, logit_0=-2., logit_wt=0.1):\n",
    "        super(LogisticRegressionModel, self).__init__() \n",
    "        self.bias = torch.nn.Parameter(logit_0*torch.ones(1, N_timepoint))\n",
    "        self.weight = torch.nn.Parameter(logit_wt * torch.randn((N_neuron, N_timepoint)))\n",
    "        # self.nl = torch.nn.Softmax(N_timepoint)\n",
    "        # self.nl = torch.nn.Sigmoid()\n",
    "\n",
    "    def forward(self, X_input):\n",
    "        # X_input.shape = N_batch, N_neuron\n",
    "        X_input[torch.isnan(X_input)] = 0.\n",
    "        # return self.nl(X_input @ self.weight +  self.bias)\n",
    "        # print(X_input.dtype, self.bias.dtype, self.weight.dtype)\n",
    "        return F.softmax(torch.tensordot(X_input, self.weight, dims=1) +  self.bias, dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "52fe6e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# F.softmax?\n",
    "# torch.tensor?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8b94dc1f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([28, 107]), torch.Size([28, 28]))"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logistic_model = LogisticRegressionModel(N_neuron, N_timepoint)\n",
    "i_trial = 3\n",
    "X = torch.tensor(fr[:, :, i_trial], dtype=torch.float64)\n",
    "outputs = logistic_model(X)\n",
    "X.shape, outputs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fb0d8d56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor([480., 490., 440., 370., 390., 520., 610., 490., 410., 440., 420., 460.,\n",
       "         480., 530., 500., 430., 480., 470., 410., 480., 540., 350., 100.,   0.,\n",
       "           0.,   0.,   0.,   0.]),\n",
       " torch.Size([28]))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.sum(axis=1), (X.sum(axis=1)>0).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0403d75d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2dd002a10>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAAE0CAYAAACLqfDEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAm/UlEQVR4nO3dfZBV9X0/8M+yCwur+xAw7LIRIhpBiw9NJRLUWlMZkGaciE5aLXbQ+NNpu1iRtomk0dQ8Uc00Oloitb9W25kQE2eCVpvYUFQYJ0CUlFo3ZlELAxEXG8w+gHXF3fP7I79s3Yjes8s592H39Zq5E/feD9/v5577veec+87Zu1VJkiQBAAAAkINxpW4AAAAAGL0EDwAAAEBuBA8AAABAbgQPAAAAQG4EDwAAAEBuBA8AAABAbgQPAAAAQG4EDwAAAEBuBA8AAABAbmpK3cCvGhgYiH379kV9fX1UVVWVuh0AAADgCJIkid7e3mhtbY1x4979uobcgoc1a9bEV7/61ejs7Iwzzzwz7r777jj77LML/rt9+/bF9OnT82oLAAAAyNDevXvj+OOPf9fHcwkevvWtb8XKlStj7dq1MW/evLjzzjtj0aJF0dHREVOnTn3Pf1tfXx8REefF70RNjH/XuupTT86s3/7nX8hsrNcvnpuqrv7F7szmTKP3Q40Fa+oeeSaz+dJsh7TzZTlWsWXZe5o1n+VazlIl916u9rfNK1jTvGZbxc5H9v7rrz5SsObEm57ObL5yPU6XqzTbq9jH8jRrJiLbdQO/lPacundGdcGatMcn5yvppTkviEi37bPc7o49xfdWHI6n4ruDn+PfTS7Bw9e+9rW49tpr4+qrr46IiLVr18a//Mu/xD/8wz/ETTfd9J7/9pe/XlET46Om6j2Ch+razPqteo95hqtm/MR0ddVvZDZnqvlS9PVe27uU8xW79yxl2XuaNZ/lWs5SJfderqpri/u+KPZ8ZG/cxCKvmTI9TperNNur2MfDNGsm6znhl9KeU1fXFg4exsK5VrGlOS+ISLfts9zujj0lkPzifwp9TULmXy755ptvxvbt22PBggX/O8m4cbFgwYLYsmXLO+r7+vqip6dnyA0AAAAYHTIPHn72s59Ff39/NDc3D7m/ubk5Ojs731G/evXqaGxsHLz5fgcAAAAYPUr+5zRXrVoV3d3dg7e9e/eWuiUAAAAgI5l/x8Nxxx0X1dXVsX///iH379+/P1paWt5RX1tbG7W12f0uDgAAAFA+Mr/iYcKECXHWWWfFxo0bB+8bGBiIjRs3xvz587OeDgAAAChjufxVi5UrV8ayZcti7ty5cfbZZ8edd94Zhw4dGvwrF0T0zmoqWFO/syuz+dKM1Z/ZbBF167P7s3rF7j1Lldw75a1hd3FXzlvnpfgTwHfk3wcjd+yekv925YhVz5mdqq6/vSPnTvKTpveehecUrKnLopn/L8s1MxZeQ9J7fUnhP8WY/jy46WhaGSLN+Xlde2bTUSJpXueIyn6t0+xz0+5vC42V9PdFPF94nFyCh9/7vd+L//7v/45bbrklOjs749d//dfjsccee8cXTgIAAACjWy7BQ0TE8uXLY/ny5XkNDwAAAFSAyr3uEgAAACh7ggcAAAAgN4IHAAAAIDeCBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3NaVuYKT62ztK3cIR1e/symys3llNmc1XrtsrjdHe++tL5mU2X117ZkMxhlTPmV3UcSr5PV3JDs4YKHULR5RmPWS1RivdB75/oGBNf4bztdzxgwxHg+FJe6yoj8L7h96U51o9J1QXrKlLNdLo17A7u71Nsc8L0rzOEZX9Wqf5HFnszw2ueAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAclNT6gZGm/72jlR1ry+ZV7CmfmdXwZreWU2p5qtrT1WWmeo5swvWpN1WWY5VbGl6TyvNeujPbLZslevrU8mKvR7e2FOfournGc5I1lo3J6VuYcTsQ2B0qVu/LbOx0p4Lp9Gwu1zPpAqr5PPlLKV5jg2zCn8OI3uueAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHJTU+oGxqr6nV1FHac/RU31nNnpxmrvSFWXld5ZTQVr6trTjZXmOWb5/NL03nNCdaqx6tYXd7uPBcVeD1nK8n0BeUizRiOs0+Eo9v4oy/OCct2XUvnSnEd94PsHUo2Vdr/FMM5fc+6DI6tbvy2zsQrtv/uTw6nGccUDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJCbmlI3MFb1t3cUrKmeM7sInYwtvbOaCtbUtWc332tXHixYU/NUY3YTMmbU7+wqXJRyH5Jmf3TsnnNSjUX52nd+VcGaD60vQiMjkOXxMM16L4WsnmOWz+/lhVNS1bWkOG6mfX7l+vqkkeY5VvLzK1cHZwwUrEm93WfNK1gyFtZyGmm2e7lKdQ4VEf35tlExCq35pL8v4vnC42R+xcNf/uVfRlVV1ZDbKaeckvU0AAAAQAXI5YqHOXPmxL/927/97yQ1LqwAAACAsSiXRKCmpiZaWlryGBoAAACoILl8ueQLL7wQra2tceKJJ8bSpUtjz54971rb19cXPT09Q24AAADA6JB58DBv3ry4//7747HHHot77rkndu3aFb/5m78Zvb29R6xfvXp1NDY2Dt6mT5+edUsAAABAiWQePCxevDg++clPxhlnnBGLFi2K7373u9HV1RXf/va3j1i/atWq6O7uHrzt3bs365YAAACAEsn9Wx+bmppi1qxZ8eKLLx7x8dra2qitrc27DQAAAKAEcvmOh7c7ePBgvPTSSzFt2rS8pwIAAADKTObBw5/92Z/Fpk2bYvfu3fGDH/wglixZEtXV1XHFFVdkPRUAAABQ5jL/VYuf/vSnccUVV8SBAwfi/e9/f5x33nmxdevWeP/735/1VJRI9ZzZBWv62zsym69+Z1fh+VKO1XNCdcGaupRjpdkOaRycMZDJOIwtad5jry+Zl2qsuvbCNcVep2nfX1nua0a7Y/fkfpEjQCppjk9167elGmvijCN/gf1I7Du/qmDN7J2ZTZcpx8PSKPbnorTKsa/Mg4cHHngg6yEBAACACuX//gAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByU1PqBvJWPWd2wZr+9o4idJKP3llNqerq2vPtI09Zvj4Nu/szGyuNN/bUF6xp3ZwUoRPGop4TqlPV1WU0XyXvS8eCD3z/QMGaLPeQ9Tu7UtWlmdPa+oVib4diHzMrnXVaGjVPNZa6hTHp2D3l+f9fp/nsNxa8vHBKwZqWlJ8PC+3b+pPDqcYpzxUDAAAAjAqCBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3NaVugHfX395RsKY+Zqcb62ibeftYKfqqnlO4rzTjlLN0/X+0YEXd+m1H3wwcQcPu7N75rZuTgjVp3vcR2b73s9rXlKJ3fmEsHC/SKMfn2HNCdaq6uhQ15fj8KJ36nV0Fa9IewVru+MFR9fJ2E2f0ZjYWpdE7qymzseraC9eU674ty3PArLjiAQAAAMiN4AEAAADIjeABAAAAyI3gAQAAAMiN4AEAAADIjeABAAAAyI3gAQAAAMiN4AEAAADIjeABAAAAyE1NqRvIW397R6lbGLHqObML1vTOako1Vn0UHquSt1Va9Tu7Ctb0ZzjfsXtke5ROzwnVqerqUtSkee9kKe3+KM1+kl9Ic7yoa8+/D8amtO/VcjwXqeTe00rzHMv1+RW795cXTklV1zLK96cNu7M8Yy4sy+N92nOa4j7DbBX7M08aPhUBAAAAuRE8AAAAALkRPAAAAAC5ETwAAAAAuRE8AAAAALkRPAAAAAC5ETwAAAAAuRE8AAAAALmpKXUDHJ36nV2p6npnNRWsqWtPN2f1nNkFa/rbO9INNsp94PsHCtb0pxzLds/eaN9eDbvTrq7CUu1D1m/LbL5iG+1rISLd8SK7FcNokOYYFjH6102a/V9EuvOoNMfyiNG/T+r4P+8rWPOhG4vQyK9on/+NgjVnPvXHReiEkUp1rBvl76+I8nyOw77iYfPmzXHxxRdHa2trVFVVxUMPPTTk8SRJ4pZbbolp06bFpEmTYsGCBfHCCy9k1S8AAABQQYYdPBw6dCjOPPPMWLNmzREfv/322+Ouu+6KtWvXxrZt2+KYY46JRYsWxRtvvHHUzQIAAACVZdi/arF48eJYvHjxER9LkiTuvPPO+NznPhef+MQnIiLin/7pn6K5uTkeeuihuPzyy4+uWwAAAKCiZPrlkrt27YrOzs5YsGDB4H2NjY0xb9682LJlyxH/TV9fX/T09Ay5AQAAAKNDpsFDZ2dnREQ0NzcPub+5uXnwsV+1evXqaGxsHLxNnz49y5YAAACAEir5n9NctWpVdHd3D9727t1b6pYAAACAjGQaPLS0tERExP79+4fcv3///sHHflVtbW00NDQMuQEAAACjQ6bBw8yZM6OlpSU2btw4eF9PT09s27Yt5s+fn+VUAAAAQAUY9l+1OHjwYLz44ouDP+/atSt27NgRkydPjhkzZsSKFSviS1/6Upx88skxc+bMuPnmm6O1tTUuueSSLPsGAAAAKsCwg4dnnnkmPvaxjw3+vHLlyoiIWLZsWdx///3x6U9/Og4dOhTXXXdddHV1xXnnnRePPfZYTJw4Mbuuh6F6zuyCNf3tHUXoZKg0fWVp3/lVBWtm7yxuT6VQ7Nc6y/l6ZzUVrKlrz2w6xpBi74+yVIr9dzGlfW0qeTtUcu+kl+YYFuE4ViqV/D7M8vzo6j2/WbCmYXd/usFGufqdXanqir21Xl44pWBNi/3MsBQ6F0n6+yKeLzzOsIOHCy64IJIkedfHq6qq4gtf+EJ84QtfGO7QAAAAwChT8r9qAQAAAIxeggcAAAAgN4IHAAAAIDeCBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3NaVugOI4dk/lZkzVc2ZnNlZ/e0dmY6Xx+pJ5BWvq1m9LNVb9zq6CNf2pRqJU0qzltGs0y/dF76ymzMaCXyr2/jZLad9flfwc03h54ZRUdS3thWt6TqhONVZdqqriSnP8jXAMHo7Z//fnBWvSbs8sz482b5lTsOZD67emHI0spT1Xadhd+NUeC/v4LD+DZKVyP40CAAAAZU/wAAAAAORG8AAAAADkRvAAAAAA5EbwAAAAAORG8AAAAADkRvAAAAAA5EbwAAAAAORG8AAAAADkpqbUDeStd1ZTwZq69vz7+FX97R0Fa6rnzM5svobd/QVr0vQUkW1fo92+86sK1nxofREaGUXSrL+0a3m06zmhOlVdmv3Da1ceLFhTV4K1bD2kV67HQ8aGNPuZcmUfAu+U5pgSUfzjSt36bYWLyvSzTJbnNGnOAetSjVR4zv7kcKpxXPEAAAAA5EbwAAAAAORG8AAAAADkRvAAAAAA5EbwAAAAAORG8AAAAADkRvAAAAAA5EbwAAAAAOSmptQNjFXVc2YXrOlv78hknIiI+p1dhedLNVJE76ymgjV17SkHSyHNdgBGLs3+4TWHi5Ko5P1f2uNTOT7HcuypFBp2pz0zIC1rK720+5Bij0VppDlXiYgIr3VERBycMVDqFt7BFQ8AAABAbgQPAAAAQG4EDwAAAEBuBA8AAABAbgQPAAAAQG4EDwAAAEBuBA8AAABAbgQPAAAAQG5qSt0AR6e/vSNVXfWc2ZnNWb+zq2BNf2azZSvNdki7TdOY/X9/Xni+lGO9vHBKwZqW9pSDVbAsX59iK3bvDbvTra40fb2x56NH2w5jTO+splR1dWNgv5VGsY9PafScUJ2qri7nPhh90uwf0pxvZq11c1L0OStV2tcnq88EWX7mSXt8qo/i7pezHOvYPeV3fcGwO9q8eXNcfPHF0draGlVVVfHQQw8Nefyqq66KqqqqIbeLLrooq34BAACACjLs4OHQoUNx5plnxpo1a9615qKLLopXXnll8PbNb37zqJoEAAAAKtOwf9Vi8eLFsXjx4vesqa2tjZaWlhE3BQAAAIwOufzyx5NPPhlTp06N2bNnxx/90R/FgQMH8pgGAAAAKHOZf7nkRRddFJdeemnMnDkzXnrppfjsZz8bixcvji1btkR19Tu/pKivry/6+voGf+7p6cm6JQAAAKBEMg8eLr/88sH/Pv300+OMM86Ik046KZ588sm48MIL31G/evXquPXWW7NuAwAAACgDuf+djRNPPDGOO+64ePHFF4/4+KpVq6K7u3vwtnfv3rxbAgAAAIok8yseftVPf/rTOHDgQEybNu2Ij9fW1kZtbW3ebQAAAAAlMOzg4eDBg0OuXti1a1fs2LEjJk+eHJMnT45bb701LrvssmhpaYmXXnopPv3pT8eHPvShWLRoUaaNAwAAAOVv2MHDM888Ex/72McGf165cmVERCxbtizuueeeePbZZ+Mf//Efo6urK1pbW2PhwoXxxS9+sWRXNdTv7CpY059/G+/QO6upYE19zC5Y09/ekUE35a16ju0QEXFwxkCpW+AoVfJaPnZP7r+ZNyLlur2I6DnhnV8ofSR1OfcxEmneqxGjf/057vxCluvB2vqFNOfnpVCunxtgNBh28HDBBRdEkiTv+vi//uu/HlVDAAAAwOhRnv8XFgAAADAqCB4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNwIHgAAAIDc1JS6gZGqnjO71C2MWf3tHaVuYcQqufdj98gJK13vrKaCNXXt6cZ6eeGUo2vm7XOmqGnY3Z/ZfJRG/c6ugjVZvsqVvGbSvFcj0r9fiahbv63ULeTOuWm20h7nPvD9AwVr0p7/vb5kXsEa7/tfSLufrI/C74tin5+nOR5GVPbnhlTviyL08XY+yQAAAAC5ETwAAAAAuRE8AAAAALkRPAAAAAC5ETwAAAAAuRE8AAAAALkRPAAAAAC5ETwAAAAAuRE8AAAAALmpKXUDeeud1VSwpq49/z5+Vc8J1QVr6tZ3FKypnjM71Xz97YXHKlfl+hoWW8sdPyh1C5SRgzMGCta0bk6K0AmVotjHgfqdXanq+vNtY0QqufcsZbkPGQvnK6SX5twuzXEuwpoplbT7yXJkzZSGKx4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNzUlLqBsaphd3+pW6gY9Tu7CtaU69bsb+8odQsjVj1ndsGaUjy/LPsq9nPMci1PnNGbourYlKMV1nNCdcGausxmIw/Ffu+kVY77mnLddxd7W9Wt35bZWJWsd1ZTqrpU+/gyXVvFlmZbTbwy3ceUYr8vXl8yL1XdaH//vLxwSqq6ljt+kHMnQ6V5rb2GpeGKBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3ggcAAAAgN4IHAAAAIDeCBwAAACA3ggcAAAAgN4IHAAAAIDc1wylevXp1fOc734mf/OQnMWnSpDjnnHPitttui9mzZw/WvPHGG/Gnf/qn8cADD0RfX18sWrQovv71r0dzc3Omjfe3d6Sqq4/ZBWv6j7aZt6meU3i+iIj6nV2Fi1KMlXY7pOkr7VjFVq59wXAVey2n2s9ERO+SeQVr3jqvu/BAd6SajhJ5eeGUgjUt7UVohDGpd1ZTqrq6MlyD+86vSlU3e2fOjTBiWZ6fp13LlXzunUbD7iw/QZG1NOu02PvbYV3xsGnTpmhra4utW7fGhg0b4vDhw7Fw4cI4dOjQYM2NN94YjzzySDz44IOxadOm2LdvX1x66aWZNw4AAACUv2Fd8fDYY48N+fn++++PqVOnxvbt2+P888+P7u7u+Pu///tYt25d/PZv/3ZERNx3331x6qmnxtatW+OjH/1odp0DAAAAZe+ovuOhu/sXl95Onjw5IiK2b98ehw8fjgULFgzWnHLKKTFjxozYsmXL0UwFAAAAVKBhXfHwdgMDA7FixYo499xz47TTTouIiM7OzpgwYUI0NTUNqW1ubo7Ozs4jjtPX1xd9fX2DP/f09Iy0JQAAAKDMjPiKh7a2tnjuuefigQceOKoGVq9eHY2NjYO36dOnH9V4AAAAQPkYUfCwfPnyePTRR+OJJ56I448/fvD+lpaWePPNN6Orq2tI/f79+6OlpeWIY61atSq6u7sHb3v37h1JSwAAAEAZGlbwkCRJLF++PNavXx+PP/54zJw5c8jjZ511VowfPz42btw4eF9HR0fs2bMn5s+ff8Qxa2tro6GhYcgNAAAAGB2G9R0PbW1tsW7dunj44Yejvr5+8HsbGhsbY9KkSdHY2BjXXHNNrFy5MiZPnhwNDQ1x/fXXx/z58/1FCwAAABiDhhU83HPPPRERccEFFwy5/7777ourrroqIiLuuOOOGDduXFx22WXR19cXixYtiq9//euZNDsS/e0dJZv7vaTpq3rO7CJ0AmNXue4f0qh5qrFgTX/7D9INNmteJvPBSFTy+zBLo/28oH5nV6q6/nzbGJFj9xzVH4HjCNKs95qnzkk5VuFjXZbvnbRrebSr5Pf0WJDm9Sn2azOs4CFJkoI1EydOjDVr1sSaNWtG3BQAAAAwOohwAQAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNzUlLqB0aa/vaPULRxRufZF+SrXNVOufRXbB75/oGBNfwXPR/ayfA29D0vDdi+NNO+dCK9P1lJv9zQ1KV+b6jmzMxtrtKvk7VC3flupW8hdOb4+rngAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByI3gAAAAAciN4AAAAAHIjeAAAAAByU5UkSVLqJt6up6cnGhsb44L4RNRUjS91OwAAAMARvJUcjifj4eju7o6GhoZ3rXPFAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkJthBQ+rV6+Oj3zkI1FfXx9Tp06NSy65JDo6OobUXHDBBVFVVTXk9od/+IeZNg0AAABUhmEFD5s2bYq2trbYunVrbNiwIQ4fPhwLFy6MQ4cODam79tpr45VXXhm83X777Zk2DQAAAFSGmuEUP/bYY0N+vv/++2Pq1Kmxffv2OP/88wfvr6uri5aWlmw6BAAAACrWUX3HQ3d3d0RETJ48ecj93/jGN+K4446L0047LVatWhWvv/760UwDAAAAVKhhXfHwdgMDA7FixYo499xz47TTThu8//d///fjgx/8YLS2tsazzz4bn/nMZ6KjoyO+853vHHGcvr6+6OvrG/y5p6dnpC0BAAAAZWbEwUNbW1s899xz8dRTTw25/7rrrhv879NPPz2mTZsWF154Ybz00ktx0kknvWOc1atXx6233jrSNgAAAIAyNqJftVi+fHk8+uij8cQTT8Txxx//nrXz5s2LiIgXX3zxiI+vWrUquru7B2979+4dSUsAAABAGRrWFQ9JksT1118f69evjyeffDJmzpxZ8N/s2LEjIiKmTZt2xMdra2ujtrZ2OG0AAAAAFWJYwUNbW1usW7cuHn744aivr4/Ozs6IiGhsbIxJkybFSy+9FOvWrYvf+Z3fiSlTpsSzzz4bN954Y5x//vlxxhln5PIEAAAAgPJVlSRJkrq4quqI9993331x1VVXxd69e+PKK6+M5557Lg4dOhTTp0+PJUuWxOc+97loaGhINUdPT080NjbGBfGJqKkan7Y1AAAAoIjeSg7Hk/FwdHd3v+dn/mH/qsV7mT59emzatGk4QwIAAACj2Ii+XBIAAAAgDcEDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQG8EDAAAAkBvBAwAAAJAbwQMAAACQm5pSN/CrkiSJiIi34nBEUuJmAAAAgCN6Kw5HxP9+jn83ZRc89Pb2RkTEU/HdEncCAAAAFNLb2xuNjY3v+nhVUiiaKLKBgYHYt29f1NfXR1VVVURE9PT0xPTp02Pv3r3R0NBQ4g4hX9Y7Y401z1hivTOWWO+MNWNxzSdJEr29vdHa2hrjxr37NzmU3RUP48aNi+OPP/6IjzU0NIyZFxCsd8Yaa56xxHpnLLHeGWvG2pp/rysdfsmXSwIAAAC5ETwAAAAAuamI4KG2tjY+//nPR21tbalbgdxZ74w11jxjifXOWGK9M9ZY8++u7L5cEgAAABg9KuKKBwAAAKAyCR4AAACA3AgeAAAAgNwIHgAAAIDcVETwsGbNmjjhhBNi4sSJMW/evPjhD39Y6pbgqK1evTo+8pGPRH19fUydOjUuueSS6OjoGFLzxhtvRFtbW0yZMiWOPfbYuOyyy2L//v0l6hiy81d/9VdRVVUVK1asGLzPemc0efnll+PKK6+MKVOmxKRJk+L000+PZ555ZvDxJEnilltuiWnTpsWkSZNiwYIF8cILL5SwYxi5/v7+uPnmm2PmzJkxadKkOOmkk+KLX/xivP077K15KtXmzZvj4osvjtbW1qiqqoqHHnpoyONp1vZrr70WS5cujYaGhmhqaoprrrkmDh48WMRnUXplHzx861vfipUrV8bnP//5+NGPfhRnnnlmLFq0KF599dVStwZHZdOmTdHW1hZbt26NDRs2xOHDh2PhwoVx6NChwZobb7wxHnnkkXjwwQdj06ZNsW/fvrj00ktL2DUcvaeffjr+9m//Ns4444wh91vvjBY///nP49xzz43x48fH9773vfjxj38cf/3Xfx3ve9/7Bmtuv/32uOuuu2Lt2rWxbdu2OOaYY2LRokXxxhtvlLBzGJnbbrst7rnnnvibv/mbeP755+O2226L22+/Pe6+++7BGmueSnXo0KE488wzY82aNUd8PM3aXrp0abS3t8eGDRvi0Ucfjc2bN8d1111XrKdQHpIyd/bZZydtbW2DP/f39yetra3J6tWrS9gVZO/VV19NIiLZtGlTkiRJ0tXVlYwfPz558MEHB2uef/75JCKSLVu2lKpNOCq9vb3JySefnGzYsCH5rd/6reSGG25IksR6Z3T5zGc+k5x33nnv+vjAwEDS0tKSfPWrXx28r6urK6mtrU2++c1vFqNFyNTHP/7x5FOf+tSQ+y699NJk6dKlSZJY84weEZGsX79+8Oc0a/vHP/5xEhHJ008/PVjzve99L6mqqkpefvnlovVeamV9xcObb74Z27dvjwULFgzeN27cuFiwYEFs2bKlhJ1B9rq7uyMiYvLkyRERsX379jh8+PCQ9X/KKafEjBkzrH8qVltbW3z84x8fsq4jrHdGl3/+53+OuXPnxic/+cmYOnVqfPjDH46/+7u/G3x8165d0dnZOWS9NzY2xrx586x3KtI555wTGzdujJ07d0ZExH/8x3/EU089FYsXL44Ia57RK83a3rJlSzQ1NcXcuXMHaxYsWBDjxo2Lbdu2Fb3nUqkpdQPv5Wc/+1n09/dHc3PzkPubm5vjJz/5SYm6guwNDAzEihUr4txzz43TTjstIiI6OztjwoQJ0dTUNKS2ubk5Ojs7S9AlHJ0HHnggfvSjH8XTTz/9jsesd0aT//qv/4p77rknVq5cGZ/97Gfj6aefjj/5kz+JCRMmxLJlywbX9JHOb6x3KtFNN90UPT09ccopp0R1dXX09/fHl7/85Vi6dGlEhDXPqJVmbXd2dsbUqVOHPF5TUxOTJ08eU+u/rIMHGCva2triueeei6eeeqrUrUAu9u7dGzfccENs2LAhJk6cWOp2IFcDAwMxd+7c+MpXvhIRER/+8Ifjueeei7Vr18ayZctK3B1k79vf/nZ84xvfiHXr1sWcOXNix44dsWLFimhtbbXmgYgo8y+XPO6446K6uvod32q+f//+aGlpKVFXkK3ly5fHo48+Gk888UQcf/zxg/e3tLTEm2++GV1dXUPqrX8q0fbt2+PVV1+N3/iN34iampqoqamJTZs2xV133RU1NTXR3NxsvTNqTJs2LX7t135tyH2nnnpq7NmzJyJicE07v2G0+PM///O46aab4vLLL4/TTz89/uAP/iBuvPHGWL16dURY84xeadZ2S0vLO/4wwltvvRWvvfbamFr/ZR08TJgwIc4666zYuHHj4H0DAwOxcePGmD9/fgk7g6OXJEksX7481q9fH48//njMnDlzyONnnXVWjB8/fsj67+joiD179lj/VJwLL7ww/vM//zN27NgxeJs7d24sXbp08L+td0aLc8899x1/Hnnnzp3xwQ9+MCIiZs6cGS0tLUPWe09PT2zbts16pyK9/vrrMW7c0I8V1dXVMTAwEBHWPKNXmrU9f/786Orqiu3btw/WPP744zEwMBDz5s0res+lUva/arFy5cpYtmxZzJ07N84+++y4884749ChQ3H11VeXujU4Km1tbbFu3bp4+OGHo76+fvB3vBobG2PSpEnR2NgY11xzTaxcuTImT54cDQ0Ncf3118f8+fPjox/9aIm7h+Gpr68f/P6SXzrmmGNiypQpg/db74wWN954Y5xzzjnxla98JX73d383fvjDH8a9994b9957b0REVFVVxYoVK+JLX/pSnHzyyTFz5sy4+eabo7W1NS655JLSNg8jcPHFF8eXv/zlmDFjRsyZMyf+/d//Pb72ta/Fpz71qYiw5qlsBw8ejBdffHHw5127dsWOHTti8uTJMWPGjIJr+9RTT42LLroorr322li7dm0cPnw4li9fHpdffnm0traW6FmVQKn/rEYad999dzJjxoxkwoQJydlnn51s3bq11C3BUYuII97uu+++wZr/+Z//Sf74j/84ed/73pfU1dUlS5YsSV555ZXSNQ0Zevuf00wS653R5ZFHHklOO+20pLa2NjnllFOSe++9d8jjAwMDyc0335w0NzcntbW1yYUXXph0dHSUqFs4Oj09PckNN9yQzJgxI5k4cWJy4oknJn/xF3+R9PX1DdZY81SqJ5544ojn7MuWLUuSJN3aPnDgQHLFFVckxx57bNLQ0JBcffXVSW9vbwmeTelUJUmSlCjzAAAAAEa5sv6OBwAAAKCyCR4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNwIHgAAAIDcCB4AAACA3AgeAAAAgNz8P3NIX/vNhpnDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1300x1300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 13))\n",
    "ax.imshow(X.numpy(), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bf3eed67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(28,\n",
       " array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "        17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]),\n",
       " array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "        17, 18, 19, 20, 21, 22]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_timepoint, np.arange(N_timepoint),  np.arange(N_timepoint)[X.sum(axis=1)>0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cdc9e161",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "        17, 18, 19, 20, 21, 22]),\n",
       " torch.Size([23, 107]),\n",
       " torch.Size([23, 28]))"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = torch.tensor(fr[:, :, i_trial], dtype=torch.float64)\n",
    "time_labels = np.arange(N_timepoint)[X.sum(axis=1)>0]\n",
    "X = X[time_labels, :]\n",
    "outputs = logistic_model(X)\n",
    "time_labels, X.shape, outputs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "98e008d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,\n",
       "         0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],\n",
       "        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "         0., 0., 0., 0., 1., 0., 0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.one_hot(torch.tensor(time_labels), num_classes=N_timepoint)*1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "928a2cc9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ddaba010>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAANnCAYAAACWLJu0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAr4klEQVR4nO3df5DV9X3v8fcuKysaWLMCu2wEBH+mUfBeoyvVeLVyBdJxYqRz1foHWq+ZKjhXibWhU39gc4ep00kdU6pz57bSzJT8mqk45t4xk5KIkxGw4nCtvVcihBQsP1Sm7IYlLD/2e//IuO1REFf2zWf37OMxcybsOV/3vP745is8/XK2oaqqKgAAAAASNJYeAAAAANQv4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0jSVHvBBfX19sWPHjhg7dmw0NDSUngMAAAAcRVVV8ctf/jI6OjqisfHY9zUMufCwY8eOmDx5cukZAAAAwMewffv2OOuss475+pALD2PHjo2IiKvii9EUpxReA1DWlj/7fOkJNc554NXSExjG3vuvl5eeUGP8/3yl9AQYNE3TppSeUOPw1m2lJ8CgaTz9tNITavT17C89od/hOBQ/jf/d/+f4Yxly4eH9v17RFKdEU4PwAIxsjWNOLT2hhusyJ2LUaOczZGlqbC49oZb/f1FHGhtGl55Qo6/hUOkJ/6b69f8c72MSfLgkAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0qSFh+XLl8fZZ58dp556anR2dsYrr7yS9VYAAADAEJUSHr773e/G4sWL45FHHonXXnstZs6cGXPmzIl33nkn4+0AAACAISolPHzjG9+Iu+66K+644474jd/4jXj66afjtNNOi7/+67/OeDsAAABgiBr08HDw4MHYsGFDzJ49+9/epLExZs+eHWvXrh3stwMAAACGsKbB/obvvfdeHDlyJNra2mqeb2trizfffPNDx/f29kZvb2//193d3YM9CQAAACik+E+1WLZsWbS0tPQ/Jk+eXHoSAAAAMEgGPTyMHz8+Ro0aFbt37655fvfu3dHe3v6h45csWRJdXV39j+3btw/2JAAAAKCQQQ8Po0ePjksvvTRWr17d/1xfX1+sXr06Zs2a9aHjm5ubY9y4cTUPAAAAoD4M+mc8REQsXrw4FixYEJ///Ofj8ssvjyeeeCJ6enrijjvuyHg7AAAAYIhKCQ8333xzvPvuu/Hwww/Hrl274pJLLokXXnjhQx84CQAAANS3lPAQEbFo0aJYtGhR1rcHAAAAhoHiP9UCAAAAqF/CAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANI0lR4A9W7UedNLT6hx5K2fl57AAJy3aH3pCTBo9n6ur/SEGhNKD2BAen/7stITajT/r38oPaHG4Z//ovSEIa1xxoWlJ9Toe/3N0hMYgL6entIThj13PAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkaaiqqio94t/r7u6OlpaWuCa+FE0Np5SeAyfshzs2lp5QY07HJaUnAAAAdeBwdShejOeiq6srxo0bd8zj3PEAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkKap9ACod3M6Lik9AWBI+Nkzl5aeUOP8OzaUnsAAbHv4N0tPqDHlsZdLTwAYNtzxAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJCmqfQAAGBk+IfZT5aeUOO2uLL0BAZgymMvl57AMNZ4+umlJ9To6+kpPQFOKnc8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKRpKj0AABgZPt04pvQEYITq27+/9AQY0dzxAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAmkEPD48++mg0NDTUPC688MLBfhsAAABgGGjK+Kaf+9zn4u///u//7U2aUt4GAAAAGOJSikBTU1O0t7dnfGsAAABgGEn5jIe33norOjo6Yvr06XHbbbfFtm3bMt4GAAAAGOIG/Y6Hzs7OWLFiRVxwwQWxc+fOWLp0aXzhC1+IN954I8aOHfuh43t7e6O3t7f/6+7u7sGeBAAAABQy6OFh3rx5/b+eMWNGdHZ2xtSpU+N73/te3HnnnR86ftmyZbF06dLBngEAAAAMAek/TvOMM86I888/PzZv3nzU15csWRJdXV39j+3bt2dPAgAAAE6S9PCwb9++2LJlS0yaNOmorzc3N8e4ceNqHgAAAEB9GPTw8MADD8SaNWviF7/4Rbz88svx5S9/OUaNGhW33nrrYL8VAAAAMMQN+mc8vP3223HrrbfGnj17YsKECXHVVVfFunXrYsKECYP9VgAAAMAQN+jh4Tvf+c5gf0sAAABgmEr/jAcAAABg5BIeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkKap9AAAYGR4fv+40hOAEWrUZ88rPaHGkf/7s9IT4KRyxwMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAmqbSAwCAkeGG07pLT6jxVOkBwElz5P+9VXoCjGjueAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABI01R6wHDQOOPC0hNq9L3+ZukJADBgvz318tITPuBw6QFQt5rOnlJ6Qo3Dv9hWekKNKetPLz2hxrbOntITqHPueAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABI01R6wHDQsL+39AQAjqahofSCWlVVesGQVh0+XHoCw9gvb76i9IQaY7+7rvSEIW37/LNKT6gxccP40hNqbL/6n0pPYBhraBo6f4xvqKqIj/Gvd3c8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKRpKj1gODiyeWvpCQxjDaeMLj2hRnXoYOkJMHiqqvQC4CRpeW5j6Qk1+koPGOI+s/y10hNqfOPNH5eeUOO+s3+z9ASGserw4dIT+lXVx9vijgcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIMODy89NJLccMNN0RHR0c0NDTEqlWral6vqioefvjhmDRpUowZMyZmz54db7311mDtBQAAAIaRAYeHnp6emDlzZixfvvyorz/++OPx5JNPxtNPPx3r16+P008/PebMmRMHDhw44bEAAADA8NI00H9g3rx5MW/evKO+VlVVPPHEE/HHf/zH8aUvfSkiIr71rW9FW1tbrFq1Km655ZYTWwsAAAAMK4P6GQ9bt26NXbt2xezZs/ufa2lpic7Ozli7du1gvhUAAAAwDAz4joePsmvXroiIaGtrq3m+ra2t/7UP6u3tjd7e3v6vu7u7B3MSAAAAUFDxn2qxbNmyaGlp6X9Mnjy59CQAAABgkAxqeGhvb4+IiN27d9c8v3v37v7XPmjJkiXR1dXV/9i+fftgTgIAAAAKGtTwMG3atGhvb4/Vq1f3P9fd3R3r16+PWbNmHfWfaW5ujnHjxtU8AAAAgPow4M942LdvX2zevLn/661bt8bGjRujtbU1pkyZEvfdd198/etfj/POOy+mTZsWDz30UHR0dMSNN944mLsBAACAYWDA4eHVV1+Na6+9tv/rxYsXR0TEggULYsWKFfHggw9GT09PfOUrX4m9e/fGVVddFS+88EKceuqpg7caAAAAGBYGHB6uueaaqKrqmK83NDTEY489Fo899tgJDQMAAACGv+I/1QIAAACoX8IDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0jSVHgD1rjp8qPQEABj2frDl5dITanzxM/+x9IQhre/AgdITatx39m+WngAjmjseAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANI0lR4wLDSOKr2gVt+R0gsYiKoqvaDGr750eekJNcY890rpCcBJMmrChNITahx5993SExiA81bdXXpCjfNifekJDMTlF5deUOuVfyy9AE4qdzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApGkqPWBY6DtSegEMmjHPbyg9ARihvvnqs6Un1Lhn6lWlJzAA5y1cX3oCw9ion20rPaGGP10w0rjjAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACBNU+kBwEnWd6T0AmCEOueUT5WeAIxQR/Z2lZ4AI5o7HgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSNJUeAPXun793cekJNabe/EbpCbWqqvQC4CSZ8Wf3lJ5QY1K8XHoCAIwI7ngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACDNgMPDSy+9FDfccEN0dHREQ0NDrFq1qub122+/PRoaGmoec+fOHay9AAAAwDAy4PDQ09MTM2fOjOXLlx/zmLlz58bOnTv7H9/+9rdPaCQAAAAwPDUN9B+YN29ezJs37yOPaW5ujvb29k88CgAAAKgPKZ/x8OKLL8bEiRPjggsuiLvvvjv27NlzzGN7e3uju7u75gEAAADUh0EPD3Pnzo1vfetbsXr16vjTP/3TWLNmTcybNy+OHDly1OOXLVsWLS0t/Y/JkycP9iQAAACgkAH/VYvjueWWW/p/ffHFF8eMGTPinHPOiRdffDGuu+66Dx2/ZMmSWLx4cf/X3d3d4gMAAADUifQfpzl9+vQYP358bN68+aivNzc3x7hx42oeAAAAQH1IDw9vv/127NmzJyZNmpT9VgAAAMAQM+C/arFv376auxe2bt0aGzdujNbW1mhtbY2lS5fG/Pnzo729PbZs2RIPPvhgnHvuuTFnzpxBHQ4AAAAMfQMOD6+++mpce+21/V+///kMCxYsiKeeeipef/31+Ju/+ZvYu3dvdHR0xPXXXx9/8id/Es3NzYO3GgAAABgWBhwerrnmmqiq6piv//CHPzyhQQAAAED9SP+MBwAAAGDkEh4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQpqn0AKh3U//LP5aeUGPU+DNLT6hx5L09pScAJ8mz/+3x0hNq3PONq0pPAIARwR0PAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmaSg84poaGXz+GgqoqvQAGzddeWV16Qo3/Pv2S0hOAk+SeqVeVnlCj8fTTS08Y0vp6ekpPqDVUfl/4Pr8/HF4aR5VeUKvvSOkFcFK54wEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgTVPpAcdUVRFRlV4Bdeex3/u90hNqjIrXSk8ARqi+np7SExiAS187UnpCjQ3/wX+/G1aqvtILYERzxQQAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABI01R6ANS7URMmlJ5Q68XXSi9gABpPO630hBp9+/eXnlCjoWlo/WusOny49IQhbVTbxNITahzZ/U7pCTWapp9dekKNwz//RekJNTZe++nSEz6gq/SAWlfMKL2g1rrXSy+o0dQxqfSEGof/ZUfpCQzAPz82q/SEGlMfXlt6woC54wEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgTVPpAVDvfn7vuaUn1Jj68LulJzAAffv3l54wpFWHD5eewAD87MHppSfUOOer75SeUOPwz39ResKQdmRvV+kJQ9sr/1R6wZBW7f9V6QkMY1MfXlt6wrDnjgcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIMKDwsW7YsLrvsshg7dmxMnDgxbrzxxti0aVPNMQcOHIiFCxfGmWeeGZ/61Kdi/vz5sXv37kEdDQAAAAwPAwoPa9asiYULF8a6deviRz/6URw6dCiuv/766Onp6T/m/vvvj+effz6+//3vx5o1a2LHjh1x0003DfpwAAAAYOhrGsjBL7zwQs3XK1asiIkTJ8aGDRvi6quvjq6urvirv/qrWLlyZfzWb/1WREQ888wz8dnPfjbWrVsXV1xxxeAtBwAAAIa8E/qMh66uroiIaG1tjYiIDRs2xKFDh2L27Nn9x1x44YUxZcqUWLt27VG/R29vb3R3d9c8AAAAgPrwicNDX19f3HfffXHllVfGRRddFBERu3btitGjR8cZZ5xRc2xbW1vs2rXrqN9n2bJl0dLS0v+YPHnyJ50EAAAADDGfODwsXLgw3njjjfjOd75zQgOWLFkSXV1d/Y/t27ef0PcDAAAAho4BfcbD+xYtWhQ/+MEP4qWXXoqzzjqr//n29vY4ePBg7N27t+auh927d0d7e/tRv1dzc3M0Nzd/khkAAADAEDegOx6qqopFixbFs88+Gz/+8Y9j2rRpNa9feumlccopp8Tq1av7n9u0aVNs27YtZs2aNTiLAQAAgGFjQHc8LFy4MFauXBnPPfdcjB07tv9zG1paWmLMmDHR0tISd955ZyxevDhaW1tj3Lhxce+998asWbP8RAsAAAAYgQYUHp566qmIiLjmmmtqnn/mmWfi9ttvj4iIP//zP4/GxsaYP39+9Pb2xpw5c+Iv//IvB2UsAAAAMLwMKDxUVXXcY0499dRYvnx5LF++/BOPAgAAAOrDJ/6pFgAAAADHIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgTVPpAVDvvjDn9dITamx7uPQCYKQ6/3+8V3pCjSOlB8BgqvpKLxjaJrSWXlDrX/+19AIGoqGh9IJaVVV6wYC54wEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgTVPpAVDv/uU/jyo9ocbWb88sPaHGtFv/T+kJQ1vj0Dp/ou9I6QUMY0c2bS49AepXVZVeMKS9858mlp5Q48yfbSk9YWjz+5+6444HAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApBEeAAAAgDTCAwAAAJBGeAAAAADSCA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaYQHAAAAII3wAAAAAKQRHgAAAIA0wgMAAACQRngAAAAA0ggPAAAAQBrhAQAAAEgjPAAAAABphAcAAAAgjfAAAAAApGkqPeCDqqqKiIjDcSiiKjwGBkFVHSw9oUbf/gOlJ9Q4XB0qPWFoq/pKL6hVHSm9AAAG7MhBv/8ZVvz+Z9g4HL8+l9//c/yxNFTHO+Ike/vtt2Py5MmlZwAAAAAfw/bt2+Oss8465utDLjz09fXFjh07YuzYsdHQ0PCJv093d3dMnjw5tm/fHuPGjRvEhXDyOZ+pJ85n6onzmXrifKaeOJ9Pjqqq4pe//GV0dHREY+OxP8lhyP1Vi8bGxo8sJQM1btw4Jxp1w/lMPXE+U0+cz9QT5zP1xPmcr6Wl5bjH+HBJAAAAII3wAAAAAKSp2/DQ3NwcjzzySDQ3N5eeAifM+Uw9cT5TT5zP1BPnM/XE+Ty0DLkPlwQAAADqR93e8QAAAACUJzwAAAAAaYQHAAAAII3wAAAAAKSp2/CwfPnyOPvss+PUU0+Nzs7OeOWVV0pPggF79NFHo6GhoeZx4YUXlp4FH8tLL70UN9xwQ3R0dERDQ0OsWrWq5vWqquLhhx+OSZMmxZgxY2L27Nnx1ltvlRkLx3G88/n222//0PV67ty5ZcbCR1i2bFlcdtllMXbs2Jg4cWLceOONsWnTpppjDhw4EAsXLowzzzwzPvWpT8X8+fNj9+7dhRbDsX2c8/maa6750PX593//9wstHrnqMjx897vfjcWLF8cjjzwSr732WsycOTPmzJkT77zzTulpMGCf+9znYufOnf2Pn/70p6UnwcfS09MTM2fOjOXLlx/19ccffzyefPLJePrpp2P9+vVx+umnx5w5c+LAgQMneSkc3/HO54iIuXPn1lyvv/3tb5/EhfDxrFmzJhYuXBjr1q2LH/3oR3Ho0KG4/vrro6enp/+Y+++/P55//vn4/ve/H2vWrIkdO3bETTfdVHA1HN3HOZ8jIu66666a6/Pjjz9eaPHIVZc/TrOzszMuu+yy+Iu/+IuIiOjr64vJkyfHvffeG1/72tcKr4OP79FHH41Vq1bFxo0bS0+BE9LQ0BDPPvts3HjjjRHx67sdOjo64qtf/Wo88MADERHR1dUVbW1tsWLFirjlllsKroWP9sHzOeLXdzzs3bv3Q3dCwFD37rvvxsSJE2PNmjVx9dVXR1dXV0yYMCFWrlwZv/M7vxMREW+++WZ89rOfjbVr18YVV1xReDEc2wfP54hf3/FwySWXxBNPPFF23AhXd3c8HDx4MDZs2BCzZ8/uf66xsTFmz54da9euLbgMPpm33norOjo6Yvr06XHbbbfFtm3bSk+CE7Z169bYtWtXzbW6paUlOjs7XasZtl588cWYOHFiXHDBBXH33XfHnj17Sk+C4+rq6oqIiNbW1oiI2LBhQxw6dKjm+nzhhRfGlClTXJ8Z8j54Pr/vb//2b2P8+PFx0UUXxZIlS2L//v0l5o1oTaUHDLb33nsvjhw5Em1tbTXPt7W1xZtvvlloFXwynZ2dsWLFirjgggti586dsXTp0vjCF74Qb7zxRowdO7b0PPjEdu3aFRFx1Gv1+6/BcDJ37ty46aabYtq0abFly5b4oz/6o5g3b16sXbs2Ro0aVXoeHFVfX1/cd999ceWVV8ZFF10UEb++Po8ePTrOOOOMmmNdnxnqjnY+R0T87u/+bkydOjU6Ojri9ddfjz/8wz+MTZs2xd/93d8VXDvy1F14gHoyb968/l/PmDEjOjs7Y+rUqfG9730v7rzzzoLLAPj3/v1fD7r44otjxowZcc4558SLL74Y1113XcFlcGwLFy6MN954w+dHUReOdT5/5Stf6f/1xRdfHJMmTYrrrrsutmzZEuecc87Jnjli1d1ftRg/fnyMGjXqQ5+8u3v37mhvby+0CgbHGWecEeeff35s3ry59BQ4Ie9fj12rqVfTp0+P8ePHu14zZC1atCh+8IMfxE9+8pM466yz+p9vb2+PgwcPxt69e2uOd31mKDvW+Xw0nZ2dERGuzydZ3YWH0aNHx6WXXhqrV6/uf66vry9Wr14ds2bNKrgMTty+fftiy5YtMWnSpNJT4IRMmzYt2tvba67V3d3dsX79etdq6sLbb78de/bscb1myKmqKhYtWhTPPvts/PjHP45p06bVvH7ppZfGKaecUnN93rRpU2zbts31mSHneOfz0bz/oe2uzydXXf5Vi8WLF8eCBQvi85//fFx++eXxxBNPRE9PT9xxxx2lp8GAPPDAA3HDDTfE1KlTY8eOHfHII4/EqFGj4tZbby09DY5r3759Nf81YevWrbFx48ZobW2NKVOmxH333Rdf//rX47zzzotp06bFQw89FB0dHTU/KQCGio86n1tbW2Pp0qUxf/78aG9vjy1btsSDDz4Y5557bsyZM6fgaviwhQsXxsqVK+O5556LsWPH9n9uQ0tLS4wZMyZaWlrizjvvjMWLF0dra2uMGzcu7r333pg1a5afaMGQc7zzecuWLbFy5cr44he/GGeeeWa8/vrrcf/998fVV18dM2bMKLx+hKnq1De/+c1qypQp1ejRo6vLL7+8WrduXelJMGA333xzNWnSpGr06NHVZz7zmermm2+uNm/eXHoWfCw/+clPqoj40GPBggVVVVVVX19f9dBDD1VtbW1Vc3Nzdd1111WbNm0qOxqO4aPO5/3791fXX399NWHChOqUU06ppk6dWt11113Vrl27Ss+GDznaeRwR1TPPPNN/zK9+9avqnnvuqT796U9Xp512WvXlL3+52rlzZ7nRcAzHO5+3bdtWXX311VVra2vV3NxcnXvuudUf/MEfVF1dXWWHj0ANVVVVJzN0AAAAACNH3X3GAwAAADB0CA8AAABAGuEBAAAASCM8AAAAAGmEBwAAACCN8AAAAACkER4AAACANMIDAAAAkEZ4AAAAANIIDwAAAEAa4QEAAABIIzwAAAAAaf4/RBM8mTSVjX4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1300x1300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 13))\n",
    "ax.imshow(outputs.detach().numpy(), cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "98dce17f",
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 0.00005\n",
    "logit_0 = -2.\n",
    "logit_wt = 0.1\n",
    "beta1, beta2 = 0.995, 0.999\n",
    "num_epochs = 2 ** 10 + 1\n",
    "amsgrad = True  # gives similar results\n",
    "amsgrad = False # gives similar results\n",
    "\n",
    "def fit_data(fr, logit_0=logit_0, logit_wt=0.1,\n",
    "                learning_rate=learning_rate,\n",
    "                num_epochs=num_epochs,\n",
    "                beta1=beta1, beta2=beta2, amsgrad=amsgrad,\n",
    "                verbose=False, \n",
    "        ):\n",
    "\n",
    "    logistic_model = LogisticRegressionModel(N_neuron, N_timepoint, logit_0=logit_0, logit_wt=logit_wt)\n",
    "    # logistic_model = logistic_model.to(device)\n",
    "    logistic_model.train()\n",
    "    optimizer = torch.optim.Adam(\n",
    "        logistic_model.parameters(), lr=learning_rate, betas=(beta1, beta2), amsgrad=amsgrad\n",
    "    )\n",
    "    for epoch in range(int(num_epochs)):\n",
    "        logistic_model.train()\n",
    "        losses = []\n",
    "        for i_trial in idx_trial[range(N_trial_train)]:\n",
    "            X = torch.tensor(fr[:, :, i_trial], dtype=torch.float64)\n",
    "            time_labels = np.arange(N_timepoint)[X.sum(axis=1)>0]\n",
    "            X = X[time_labels, :]\n",
    "            labels = F.one_hot(torch.tensor(time_labels), num_classes=N_timepoint)*1.\n",
    "            # X, labels = X.to(device), labels.to(device)\n",
    "            outputs = logistic_model(X)\n",
    "            loss = criterion(outputs, labels)\n",
    "\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            losses.append(loss.item())\n",
    "\n",
    "        if verbose and (epoch % (num_epochs // 32) == 0):\n",
    "            print(f\"Iteration: {epoch} - Loss: {np.mean(losses):.5f}\")\n",
    "\n",
    "    return logistic_model, np.mean(losses)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "36eb6b45",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration: 0 - Loss: 0.65640\n",
      "Iteration: 32 - Loss: 0.13167\n",
      "Iteration: 64 - Loss: 0.13023\n",
      "Iteration: 96 - Loss: 0.13010\n",
      "Iteration: 128 - Loss: 0.13001\n",
      "Iteration: 160 - Loss: 0.12995\n",
      "Iteration: 192 - Loss: 0.12990\n",
      "Iteration: 224 - Loss: 0.12986\n",
      "Iteration: 256 - Loss: 0.12983\n",
      "Iteration: 288 - Loss: 0.12981\n",
      "Iteration: 320 - Loss: 0.12979\n",
      "Iteration: 352 - Loss: 0.12978\n",
      "Iteration: 384 - Loss: 0.12977\n",
      "Iteration: 416 - Loss: 0.12976\n",
      "Iteration: 448 - Loss: 0.12975\n",
      "Iteration: 480 - Loss: 0.12974\n",
      "Iteration: 512 - Loss: 0.12974\n",
      "Iteration: 544 - Loss: 0.12973\n",
      "Iteration: 576 - Loss: 0.12973\n",
      "Iteration: 608 - Loss: 0.12973\n",
      "Iteration: 640 - Loss: 0.12972\n",
      "Iteration: 672 - Loss: 0.12972\n",
      "Iteration: 704 - Loss: 0.12972\n",
      "Iteration: 736 - Loss: 0.12972\n",
      "Iteration: 768 - Loss: 0.12972\n",
      "Iteration: 800 - Loss: 0.12972\n",
      "Iteration: 832 - Loss: 0.12972\n",
      "Iteration: 864 - Loss: 0.12971\n",
      "Iteration: 896 - Loss: 0.12971\n",
      "Iteration: 928 - Loss: 0.12971\n",
      "Iteration: 960 - Loss: 0.12971\n",
      "Iteration: 992 - Loss: 0.12971\n",
      "Iteration: 1024 - Loss: 0.12971\n"
     ]
    }
   ],
   "source": [
    "\n",
    "logistic_model, loss = fit_data(fr, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "63d1dd5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((107, 28), (2996,))"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logistic_model.weight.detach().numpy().shape, logistic_model.weight.detach().numpy().ravel().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "0bb4ce15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((28, 28),\n",
       " (28, 28),\n",
       " array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "        1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    counts = np.zeros((N_timepoint, N_timepoint))\n",
    "    outputs = np.zeros((N_timepoint, N_timepoint))\n",
    "    for i_trial in idx_trial[range(N_trial_train, N_trial)]:\n",
    "        X = torch.tensor(fr[:, :, i_trial], dtype=torch.float64)\n",
    "        time_labels = np.arange(N_timepoint)[X.sum(axis=1)>0]\n",
    "        counts[time_labels, :] += 1\n",
    "        X = X[time_labels, :]\n",
    "        output = logistic_model(X).numpy()\n",
    "        outputs[time_labels, :] += output\n",
    "\n",
    "    outputs = outputs/counts\n",
    "\n",
    "outputs.shape, counts.shape, outputs.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b7757287",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.001887879487640713, 0.21387404554392864)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "outputs.min(), outputs.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2252651e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(60.0, 300.0)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts.min(), counts.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "e3654260",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(tensor(0),)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.isnan(torch.tensor(fr[:, :, i_trial].astype(float))).sum(), "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "47780c09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2e9018e10>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABBgAAAQUCAYAAADTI6oSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtl0lEQVR4nO3df4xX9Z3v8fcAMmKdGTsiDLMOFLVqt6jNtUqJLbGBCDQx9UcTtd1EjbGpiybIunbZ1F/ZJmT9o2u6y+pfK9uktN0mVVOz68ZiGdIs2GivIf4hKyy94EWoNcuMjHVE53v/6HW2UwEdXjN+h+HxSE7CnO9nvudtcnKSeXq+59vSaDQaBQAAABCY0uwBAAAAgOOfwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGLTmj3AHxsaGqq9e/dWW1tbtbS0NHscAAAAOGE1Go164403qru7u6ZMOfo9ChMuMOzdu7d6enqaPQYAAADw/+3Zs6fOPPPMo66ZcIGhra2tqqr+z68+Ue2n+gQHAAAANEv/waGa979+Pfy3+tFMuMDw3sci2k+dUu1tAgMAAAA024d5hIG/4AEAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgNm6BYd26dfWJT3yiTj755Fq4cGH98pe/HK9DAQAAAE02LoHhRz/6Ua1evbruu++++tWvflUXXXRRLVu2rH7zm9+Mx+EAAACAJhuXwPCd73ynbr311rr55pvrT//0T+uRRx6pU045pf7pn/5pPA4HAAAANNmYB4a33367nn/++Vq6dOn/HGTKlFq6dGlt2bLlfesHBwerv79/xAYAAAAcX8Y8MPz2t7+td999t2bPnj1i/+zZs2vfvn3vW7927drq6OgY3np6esZ6JAAAAGCcNf1bJNasWVN9fX3D2549e5o9EgAAADBK08b6DWfOnFlTp06t/fv3j9i/f//+6urqet/61tbWam1tHesxAAAAgI/QmN/BMH369Lr44otr48aNw/uGhoZq48aNtWjRorE+HAAAADABjPkdDFVVq1evrhtvvLE++9nP1qWXXloPPfRQDQwM1M033zwehwMAAACabFwCw3XXXVevvfZa3XvvvbVv3776zGc+U0899dT7HvwIAAAATA4tjUaj0ewh/lB/f391dHTUf//nWdXe1vRnUAIAAMAJq/+Nofr4uf9VfX191d7eftS1/oIHAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABAbMwDw/33318tLS0jtvPPP3+sDwMAAABMINPG400//elP189+9rP/Oci0cTkMAAAAMEGMy1/+06ZNq66urvF4awAAAGACGpdnMLz88svV3d1dZ511Vn3ta1+r3bt3H3Ht4OBg9ff3j9gAAACA48uYB4aFCxfW+vXr66mnnqqHH364du3aVV/4whfqjTfeOOz6tWvXVkdHx/DW09Mz1iMBAAAA46yl0Wg0xvMABw4cqHnz5tV3vvOduuWWW973+uDgYA0ODg7/3N/fXz09PfXf/3lWtbf5kgsAAABolv43hurj5/5X9fX1VXt7+1HXjvvTF0877bQ699xza8eOHYd9vbW1tVpbW8d7DAAAAGAcjfstAgcPHqydO3fWnDlzxvtQAAAAQJOMeWC46667qre3t37961/Xf/zHf9TVV19dU6dOrRtuuGGsDwUAAABMEGP+EYlXXnmlbrjhhnr99dfrjDPOqM9//vO1devWOuOMM8b6UAAAAMAEMeaB4Yc//OFYvyUAAAAwwfmaBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMRGHRg2b95cV155ZXV3d1dLS0s9/vjjI15vNBp177331pw5c2rGjBm1dOnSevnll8dqXgAAAGACGnVgGBgYqIsuuqjWrVt32NcffPDB+u53v1uPPPJIPfvss/Wxj32sli1bVm+99VY8LAAAADAxTRvtL6xYsaJWrFhx2NcajUY99NBD9a1vfau+/OUvV1XV9773vZo9e3Y9/vjjdf3112fTAgAAABPSmD6DYdeuXbVv375aunTp8L6Ojo5auHBhbdmy5bC/Mzg4WP39/SM2AAAA4PgypoFh3759VVU1e/bsEftnz549/NofW7t2bXV0dAxvPT09YzkSAAAA8BFo+rdIrFmzpvr6+oa3PXv2NHskAAAAYJTGNDB0dXVVVdX+/ftH7N+/f//wa3+stbW12tvbR2wAAADA8WVMA8P8+fOrq6urNm7cOLyvv7+/nn322Vq0aNFYHgoAAACYQEb9LRIHDx6sHTt2DP+8a9eueuGFF6qzs7Pmzp1bq1atqm9/+9v1yU9+subPn1/33HNPdXd311VXXTWWcwMAAAATyKgDw3PPPVdf/OIXh39evXp1VVXdeOONtX79+rr77rtrYGCgvv71r9eBAwfq85//fD311FN18sknj93UAAAAwITS0mg0Gs0e4g/19/dXR0dH/fd/nlXtbU1/BiUAAACcsPrfGKqPn/tf1dfX94HPTPQXPAAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiE1r9gBHsul3U+pj06Y2ewwAAAA4YQ387sOvdQcDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACITWv2AEfys/4F1Tp0UrPHAAAAgBPW4MFDVbXjQ611BwMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACA2KgDw+bNm+vKK6+s7u7uamlpqccff3zE6zfddFO1tLSM2JYvXz5W8wIAAAAT0KgDw8DAQF100UW1bt26I65Zvnx5vfrqq8PbD37wg2hIAAAAYGKbNtpfWLFiRa1YseKoa1pbW6urq+tDvd/g4GANDg4O/9zf3z/akQAAAIAmG5dnMGzatKlmzZpV5513Xt122231+uuvH3Ht2rVrq6OjY3jr6ekZj5EAAACAcTTmgWH58uX1ve99rzZu3Fh/+7d/W729vbVixYp69913D7t+zZo11dfXN7zt2bNnrEcCAAAAxtmoPyLxQa6//vrhf19wwQV14YUX1tlnn12bNm2qJUuWvG99a2trtba2jvUYAAAAwEdo3L+m8qyzzqqZM2fWjh07xvtQAAAAQJOMe2B45ZVX6vXXX685c+aM96EAAACAJhn1RyQOHjw44m6EXbt21QsvvFCdnZ3V2dlZDzzwQF177bXV1dVVO3furLvvvrvOOeecWrZs2ZgODgAAAEwcow4Mzz33XH3xi18c/nn16tVVVXXjjTfWww8/XNu2bat//ud/rgMHDlR3d3ddccUV9Td/8zeeswAAAACT2KgDw+WXX16NRuOIr//7v/97NBAAAABw/Bn3ZzAAAAAAk5/AAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMSmNXuAI3nut3Nr2u9amz0GAAAAnLDeGRj80GvdwQAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGLTmj3Akfzf3afXlBknN3sMAAAAOGEN/e6tD73WHQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAAAxgQEAAACICQwAAABATGAAAAAAYgIDAAAAEBMYAAAAgJjAAAAAAMQEBgAAACAmMAAAAACxac0e4Ei6eqfWtJOmNnsMAAAAOGG9c2hqvfIh17qDAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACA2KgCw9q1a+uSSy6ptra2mjVrVl111VW1ffv2EWveeuutWrlyZZ1++ul16qmn1rXXXlv79+8f06EBAACAiWVUgaG3t7dWrlxZW7duraeffroOHTpUV1xxRQ0MDAyvufPOO+unP/1p/fjHP67e3t7au3dvXXPNNWM+OAAAADBxtDQajcax/vJrr71Ws2bNqt7e3lq8eHH19fXVGWecURs2bKivfOUrVVX10ksv1ac+9anasmVLfe5zn/vA9+zv76+Ojo767DXfrmknnXysowEAAAChdw69Vc/95FvV19dX7e3tR10bPYOhr6+vqqo6Ozurqur555+vQ4cO1dKlS4fXnH/++TV37tzasmXLYd9jcHCw+vv7R2wAAADA8eWYA8PQ0FCtWrWqLrvsslqwYEFVVe3bt6+mT59ep5122oi1s2fPrn379h32fdauXVsdHR3DW09Pz7GOBAAAADTJMQeGlStX1osvvlg//OEPowHWrFlTfX19w9uePXui9wMAAAA+etOO5Zduv/32evLJJ2vz5s115plnDu/v6uqqt99+uw4cODDiLob9+/dXV1fXYd+rtbW1Wltbj2UMAAAAYIIY1R0MjUajbr/99nrsscfqmWeeqfnz5494/eKLL66TTjqpNm7cOLxv+/bttXv37lq0aNHYTAwAAABMOKO6g2HlypW1YcOGeuKJJ6qtrW34uQodHR01Y8aM6ujoqFtuuaVWr15dnZ2d1d7eXnfccUctWrToQ32DBAAAAHB8GlVgePjhh6uq6vLLLx+x/9FHH62bbrqpqqr+7u/+rqZMmVLXXnttDQ4O1rJly+of//Efx2RYAAAAYGIaVWBoNBofuObkk0+udevW1bp16455KAAAAOD4cszfIgEAAADwHoEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiAkMAAAAQExgAAAAAGICAwAAABATGAAAAICYwAAAAADEBAYAAAAgJjAAAAAAMYEBAAAAiE1r9gBHctr/3l/TprQ2ewwAAAA4Yb0zNPih17qDAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAIDat2QMcyTu7dle1nNTsMQAAAOCE9U7j0Ide6w4GAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiowoMa9eurUsuuaTa2tpq1qxZddVVV9X27dtHrLn88surpaVlxPaNb3xjTIcGAAAAJpZRBYbe3t5auXJlbd26tZ5++uk6dOhQXXHFFTUwMDBi3a233lqvvvrq8Pbggw+O6dAAAADAxDJtNIufeuqpET+vX7++Zs2aVc8//3wtXrx4eP8pp5xSXV1dYzMhAAAAMOFFz2Do6+urqqrOzs4R+7///e/XzJkza8GCBbVmzZp68803j/geg4OD1d/fP2IDAAAAji+juoPhDw0NDdWqVavqsssuqwULFgzv/+pXv1rz5s2r7u7u2rZtW33zm9+s7du3109+8pPDvs/atWvrgQceONYxAAAAgAmgpdFoNI7lF2+77bb6t3/7t/rFL35RZ5555hHXPfPMM7VkyZLasWNHnX322e97fXBwsAYHB4d/7u/vr56enrq8vlzTWk46ltEAAACAMfBO41Btqieqr6+v2tvbj7r2mO5guP322+vJJ5+szZs3HzUuVFUtXLiwquqIgaG1tbVaW1uPZQwAAABgghhVYGg0GnXHHXfUY489Vps2bar58+d/4O+88MILVVU1Z86cYxoQAAAAmPhGFRhWrlxZGzZsqCeeeKLa2tpq3759VVXV0dFRM2bMqJ07d9aGDRvqS1/6Up1++um1bdu2uvPOO2vx4sV14YUXjst/AAAAANB8o3oGQ0tLy2H3P/roo3XTTTfVnj176s/+7M/qxRdfrIGBgerp6amrr766vvWtb33gZzXe09/fXx0dHZ7BAAAAAE02bs9g+KAW0dPTU729vaN5SwAAAGASmNLsAQAAAIDjn8AAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACA2LRmD/DHGo1GVVW9U4eqGk0eBgAAAE5g79Shqvqfv9WPZsIFhjfeeKOqqn5R/9rkSQAAAICq3/+t3tHRcdQ1LY0PkyE+QkNDQ7V3795qa2urlpaWY36f/v7+6unpqT179lR7e/sYTggfLecyk4VzmcnAecxk4VxmsnAuj79Go1FvvPFGdXd315QpR3/KwoS7g2HKlCl15plnjtn7tbe3O9GYFJzLTBbOZSYD5zGThXOZycK5PL4+6M6F93jIIwAAABATGAAAAIDYpA0Mra2tdd9991Vra2uzR4GIc5nJwrnMZOA8ZrJwLjNZOJcnlgn3kEcAAADg+DNp72AAAAAAPjoCAwAAABATGAAAAICYwAAAAADEBAYAAAAgNmkDw7p16+oTn/hEnXzyybVw4cL65S9/2eyRYFTuv//+amlpGbGdf/75zR4Ljmrz5s115ZVXVnd3d7W0tNTjjz8+4vVGo1H33ntvzZkzp2bMmFFLly6tl19+uTnDwlF80Ll80003ve8avXz58uYMC0ewdu3auuSSS6qtra1mzZpVV111VW3fvn3EmrfeeqtWrlxZp59+ep166ql17bXX1v79+5s0MRzehzmXL7/88vddl7/xjW80aeIT16QMDD/60Y9q9erVdd9999WvfvWruuiii2rZsmX1m9/8ptmjwah8+tOfrldffXV4+8UvftHskeCoBgYG6qKLLqp169Yd9vUHH3ywvvvd79YjjzxSzz77bH3sYx+rZcuW1VtvvfURTwpH90HnclXV8uXLR1yjf/CDH3yEE8IH6+3trZUrV9bWrVvr6aefrkOHDtUVV1xRAwMDw2vuvPPO+ulPf1o//vGPq7e3t/bu3VvXXHNNE6eG9/sw53JV1a233jriuvzggw82aeITV0uj0Wg0e4ixtnDhwrrkkkvqH/7hH6qqamhoqHp6euqOO+6ov/qrv2rydPDh3H///fX444/XCy+80OxR4Ji0tLTUY489VldddVVV/f7uhe7u7vqLv/iLuuuuu6qqqq+vr2bPnl3r16+v66+/vonTwpH98blc9fs7GA4cOPC+OxtgInvttddq1qxZ1dvbW4sXL66+vr4644wzasOGDfWVr3ylqqpeeuml+tSnPlVbtmypz33uc02eGA7vj8/lqt/fwfCZz3ymHnrooeYOd4KbdHcwvP322/X888/X0qVLh/dNmTKlli5dWlu2bGniZDB6L7/8cnV3d9dZZ51VX/va12r37t3NHgmO2a5du2rfvn0jrs8dHR21cOFC12eOS5s2bapZs2bVeeedV7fddlu9/vrrzR4Jjqqvr6+qqjo7O6uq6vnnn69Dhw6NuC6ff/75NXfuXNdlJrQ/Ppff8/3vf79mzpxZCxYsqDVr1tSbb77ZjPFOaNOaPcBY++1vf1vvvvtuzZ49e8T+2bNn10svvdSkqWD0Fi5cWOvXr6/zzjuvXn311XrggQfqC1/4Qr344ovV1tbW7PFg1Pbt21dVddjr83uvwfFi+fLldc0119T8+fNr586d9dd//de1YsWK2rJlS02dOrXZ48H7DA0N1apVq+qyyy6rBQsWVNXvr8vTp0+v0047bcRa12UmssOdy1VVX/3qV2vevHnV3d1d27Ztq29+85u1ffv2+slPftLEaU88ky4wwGSxYsWK4X9feOGFtXDhwpo3b179y7/8S91yyy1NnAyAP/xIzwUXXFAXXnhhnX322bVp06ZasmRJEyeDw1u5cmW9+OKLnufEce9I5/LXv/714X9fcMEFNWfOnFqyZEnt3Lmzzj777I96zBPWpPuIxMyZM2vq1Knve/rt/v37q6urq0lTQe60006rc889t3bs2NHsUeCYvHcNdn1mMjrrrLNq5syZrtFMSLfffns9+eST9fOf/7zOPPPM4f1dXV319ttv14EDB0asd11mojrSuXw4CxcurKpyXf6ITbrAMH369Lr44otr48aNw/uGhoZq48aNtWjRoiZOBpmDBw/Wzp07a86cOc0eBY7J/Pnzq6ura8T1ub+/v5599lnXZ457r7zySr3++uuu0UwojUajbr/99nrsscfqmWeeqfnz5494/eKLL66TTjppxHV5+/bttXv3btdlJpQPOpcP570Hpbsuf7Qm5UckVq9eXTfeeGN99rOfrUsvvbQeeuihGhgYqJtvvrnZo8GHdtddd9WVV15Z8+bNq71799Z9991XU6dOrRtuuKHZo8ERHTx4cMT/Kdi1a1e98MIL1dnZWXPnzq1Vq1bVt7/97frkJz9Z8+fPr3vuuae6u7tHPJ0fJoKjncudnZ31wAMP1LXXXltdXV21c+fOuvvuu+ucc86pZcuWNXFqGGnlypW1YcOGeuKJJ6qtrW34uQodHR01Y8aM6ujoqFtuuaVWr15dnZ2d1d7eXnfccUctWrTIN0gwoXzQubxz587asGFDfelLX6rTTz+9tm3bVnfeeWctXry4LrzwwiZPf4JpTFJ///d/35g7d25j+vTpjUsvvbSxdevWZo8Eo3Ldddc15syZ05g+fXrjT/7kTxrXXXddY8eOHc0eC47q5z//eaOq3rfdeOONjUaj0RgaGmrcc889jdmzZzdaW1sbS5YsaWzfvr25Q8NhHO1cfvPNNxtXXHFF44wzzmicdNJJjXnz5jVuvfXWxr59+5o9NoxwuHO4qhqPPvro8Jrf/e53jT//8z9vfPzjH2+ccsopjauvvrrx6quvNm9oOIwPOpd3797dWLx4caOzs7PR2traOOeccxp/+Zd/2ejr62vu4Ceglkaj0fgogwYAAAAw+Uy6ZzAAAAAAHz2BAQAAAIgJDAAAAEBMYAAAAABiAgMAAAAQExgAAACAmMAAAAAAxAQGAAAAICYwAAAAADGBAQAAAIgJDAAAAEDs/wGTxx0PVnxNIQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1300x1300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 13))\n",
    "ax.imshow(counts, cmap='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "65f764f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'true time')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABCsAAAQoCAYAAAAjX1T/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMFklEQVR4nO3de5TfBX3n//dnvnNJIMlAgNwkhAAKKLeKGHNQLiWHm4eWwu5BWpfLoWBtgkBKQVYRaKl0PbuCdlF61hV0j6jtCtjy22VrEUJZIRYwsOxqgDQISAIIJMMkmcnM9/v5/UH54kdyG8nM503m8Thnzsl855vkxef7+V7yZC5FWZZlAAAAACTRUfcAAAAAgF8lVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkEpn3QNGW6vViueffz4mT54cRVHUPQcAAADGrbIs47XXXotZs2ZFR8fmP39ih48Vzz//fMyePbvuGQAAAMC/evbZZ2PPPffc7Md3+FgxefLkiIj4+SN7x5RJOb7q5fO/PKDuCW3fuX9+3RMq9ngo12e/THpusO4JFV2vbqh7wpuGhuteUFEM59oTPpNrs8otFPRalGXdC1IrMh2fZrPuBVWtRMcmIt3jcqrHwUzncULlxo11T6jK9DyR7Nwph4bqnpBauTHR8Un2HFEO5zk2w+VQ3B//X/vf6puzw8eKN770Y8qkjpgyOccDX89AV90T2jomTqh7QkWjO9ELm4jo7Ey2p9Gqe8KbWo26F1QUyfakepGeTNnI8VjcluyFaDZFphdbZbJYEYmOTUSuf+BF5HocdD/fojLRTRURuc7lZOdOmel+lVCqc7nIdu7UveDXlLHVb9OQ6JEAAAAAQKwAAAAAkhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUnlHxIobb7wx9t5775gwYULMmzcvfvzjH9c9CQAAABgl6WPFd7/73Vi8eHFcddVV8cgjj8Shhx4aJ5xwQrz44ot1TwMAAABGQfpY8cUvfjHOP//8OPfcc+O9731v3HTTTbHTTjvF17/+9bqnAQAAAKMgdazYuHFjPPzww7FgwYL2ZR0dHbFgwYJ44IEHalwGAAAAjJbOugdsyS9/+ctoNpsxffr0yuXTp0+Pn/3sZ5v8PYODgzE4ONh+v6+vb1Q3AgAAANtX6s+s+E1cd9110dvb236bPXt23ZMAAACAEUgdK3bfffdoNBrxwgsvVC5/4YUXYsaMGZv8PVdccUWsXbu2/fbss8+OxVQAAABgO0kdK7q7u+Pwww+Pu+++u31Zq9WKu+++O+bPn7/J39PT0xNTpkypvAEAAADvHKm/Z0VExOLFi+Pss8+OD3zgA/HBD34wbrjhhli3bl2ce+65dU8DAAAARkH6WHHGGWfESy+9FJ/73Odi9erVcdhhh8Vdd931lm+6CQAAAOwY0seKiIhFixbFokWL6p4BAAAAjIHU37MCAAAAGH/ECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUumse8BYuebF90XP+q66Z0RExON9s+qe0Fbu1Kx7QsXgLt11T6joWZvjnHlDx1BP3RPaOjYka52NZHtaZd0Lqlqtuhe8qSjqXpBbmevcKRPdt4pkxyYi0f0qIqLRqHtBVab7erZzp5nr9Vc6mW6vbLdVpufziFy3FTucPK9AAAAAAEKsAAAAAJIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFLprHvAWLnt/nnRMWFC3TMiIqLsLOue0FZMHq57QkX/7DzHJiKiaDXqnlBV5DiHIyK6Xsv18NFYP1T3hIpiqFX3hIpiqFn3hDfJ5FtUDOZ6XI6Oou4FbWWZ7Dki2Z4o8txWERHRTPS4k+y2Kpu5niOizLWn3JhoT6bzOCKilejYRL7H5Wgl28Pb4iUjAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpdNY9YKzs8XBEo6vuFa8bnpinEa19T3fdEyqaew7UPaHi1SlJTpp/1ZyQ5y670+o853FExIRXG3VPqOhcN1z3hIqOoWbdE97UKuteUFE0W3VPqCiauY5PtPIcn2I40XkcETGU635ebhyqe0LVcK7jk0oz17lcJtuT6fiUyZ4joky2J5ky2WuMVMpEx2Ybt+T61wYAAAAw7okVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJBKZ90DxkrPq8PR2Tlc94yIiGh152lEwzvlOgXW7NJV94SKrl0H655QsX5mo+4JbR1Dec7jiIjGxrLuCb8m132rsz/P8ekYbtU9oapV1L2gqsxzW0VERCvRnnTHJtm5XCbb05HovtVs1r2gosx2Ljs+m5fpPI6IcijRscko2+Mgb0uuf20AAAAA455YAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpdNY9YKx09Q1FZ2ej7hkREdHqybEjIqLn1Vy9auIvcp2SA4O5jk+jrHvBmwZ7615QVTRz3VY9a4q6J1RMKPOcPJ3rm3VPqMh15kSUXXmeIyIiikTnThS57lfRyHVbRatV9wK2UbIzOcqhuhdUFdnu65lke9xJpsz1EiOXd+DByfYaDQAAABjnxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFJJHSuuvvrqKIqi8nbAAQfUPQsAAAAYRZ11D9ia973vffGP//iP7fc7O9NPBgAAAN6G9P/y7+zsjBkzZtQ9AwAAABgjqb8MJCLiySefjFmzZsU+++wTf/AHfxDPPPPMFq8/ODgYfX19lTcAAADgnSN1rJg3b17ccsstcdddd8VXv/rVWLlyZXzkIx+J1157bbO/57rrrove3t722+zZs8dwMQAAAPB2FWVZlnWP2FZr1qyJOXPmxBe/+MU477zzNnmdwcHBGBwcbL/f19cXs2fPjqM/+Jno7JwwVlO3qNXTqHtCW/+7uuueUNG3T65+NrB7q+4JFY2Bou4JbV19ebZEREx4JddDWc+aXHsmvDJc94S2zvXNuidUdAwl27NhqO4JFcXGPOdOsWFw61caS8O5zp1yKNe5k8pwnvM4IiKayc6dgWT3LTarbOZ6bZpNmey+lUorz7EZLofi3vh+rF27NqZMmbLZ66X/nhW/apdddon3vOc98dRTT232Oj09PdHT0zOGqwAAAIDtKdf/xt6K/v7+WLFiRcycObPuKQAAAMAoSR0rLr300liyZEk8/fTT8aMf/Sh+7/d+LxqNRpx55pl1TwMAAABGSeovA3nuuefizDPPjJdffjn22GOP+PCHPxwPPvhg7LHHHnVPAwAAAEZJ6ljxne98p+4JAAAAwBhL/WUgAAAAwPgjVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKp11DxgrHc1WdESr7hmv6yjqXtA2YU2z7gkVzV/kOTYREY2BXD2v2V33gjeVXXUvqBrcJde5UyR5uHlDY6hR94S2spHrtsp2P+9slXVPqEh1dIZzPWcVZa7bqiiTvaxr5XkgLDvyPAZGREQz17kcXcme1IeG6l7QVia7n6dT5rmfR0QUif6dVSZ7Pn8nSvUaBAAAAECsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUOuseMFaKoVYUrWbdMyIioijqXvCmotlV94SKrvVl3ROqss1JdI/dODnRiRwRwxPrXlC1sTfX8YkyT5vu6sl1bLr7c+3paCa6o0dEkehxsNg4XPeEqmar7gVV2fZ05HncicjxGrCtleiOFRHRSnbuwI6gdL96uzI9iwAAAACIFQAAAEAuYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQSmfdA8ZMWb7+lkDRyrEjIqIx0Kx7QkX3a3UvqOocKOqeUNHqzLOnGM7VOotWnmMTETE8se4FVRun5Dk+ra48WyIiijwPyRER0RhM9tTczHOAig2NuidUdeU5NhERRZLXOW2tVt0L2opEWyIiylau+3nRketxOdOZXDRzvVYuUx2dfMpE/87i7cv1rw0AAABg3BMrAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACCVzroHjJlmGRFl3Ste15FkR0QUzTxbIiKK4Vx7OgdbdU9Iq2O4UfeEX5NrT6uzqHtCRXNC3QveVOa6qaJzQ67banhirv+PULTyvFToGO6ue0JFx/q6FyTXbNa9oK1o5Xp9EUNDdS+oKnI97rAFHbmesyLP3ZwdkEcmAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEils+4BY6UY3BhFo6h7RkREFM1G3RPelOSYvKHVpZ9tSUezVfeEtnJDrnOnc2Kuc6cxmOv4tLrqXvCmMtdNFUM757qtOgdy7SkS3WCNwUTPnxFRDOd6GVUMNeueUJXo3ImOZPerzlznTjk8UPeEqo5E5042w8N1L6hqJHtcjjyPg2Ur2Xlc5jk22yrZEQQAAADGO7ECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAglVpjxX333RennHJKzJo1K4qiiDvuuKPy8bIs43Of+1zMnDkzJk6cGAsWLIgnn3yynrEAAADAmKg1Vqxbty4OPfTQuPHGGzf58S984Qvx5S9/OW666aZYunRp7LzzznHCCSfEwMDAGC8FAAAAxkpnnX/5SSedFCeddNImP1aWZdxwww3x2c9+Nn73d383IiK++c1vxvTp0+OOO+6Ij33sY2M5FQAAABgjab9nxcqVK2P16tWxYMGC9mW9vb0xb968eOCBBzb7+wYHB6Ovr6/yBgAAALxzpI0Vq1evjoiI6dOnVy6fPn16+2Obct1110Vvb2/7bfbs2aO6EwAAANi+0saK39QVV1wRa9eubb89++yzdU8CAAAARiBtrJgxY0ZERLzwwguVy1944YX2xzalp6cnpkyZUnkDAAAA3jnSxoq5c+fGjBkz4u67725f1tfXF0uXLo358+fXuAwAAAAYTbX+NJD+/v546qmn2u+vXLkyli1bFlOnTo299torLr744rj22mvj3e9+d8ydOzeuvPLKmDVrVpx66qn1jQYAAABGVa2x4qGHHopjjz22/f7ixYsjIuLss8+OW265JS677LJYt25dXHDBBbFmzZr48Ic/HHfddVdMmDChrskAAADAKCvKsizrHjGa+vr6ore3N47b96LobPTUPed1nY26F7S1dk5yTP7V0JRce7LpaLbqntA2PCHPeRwRMTC11vb6Fht2z/VVdkOT6l6QV1d/3QuqJr6c534eEdG1Ps+enleH6p5Q0ViXa0/HusG6J1Qles4qBjfWPaEq2Z5yYKDuCRVlonMnWom2REQ5PFz3hNyazboXtJWtZP/MbuU5NsPlUNwb34+1a9du8XtM5no1DQAAAIx7YgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApNJZ94CxUrRaUUSr7hkREVEWRd0T0ipaZd0TKopmtj05zuGIiEay87hrQ55jExHR7Mt1fIpmnj2t7roXVJWNuhdUNbvz3FYREWWR5/9rdAwme9mS62EnioHhuidUpDqTO/KcxxERMXFC3QuqhobqXlBRJHrciUaiLRER6zfUvSC11lCex8GiI9WjYEQjzwuwoiwituFhJ9m9DwAAABjvxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASKWz7gFjpWw0omw06p7xus4kOyKi1ZVnS0REUZZ1T6joGByue0JF0cpzfMpGrtbZsTHPsYmI6BzItadMdFcfLoq6J1S0uupeUNXsqXvBrynz3F7NCbkedxqDufaUE5K9rBuoe8CvSPacVWwYrHtCRdGV64GwHBqqe8KbhnO9Foxkz6HRbNa9oKLoSHZ8EimH89yvynLbtuR65AYAAADGPbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFLprHvAmOlsRDQada+IiIiyI08j6hhq1j2homy26p5Q0bEx2fEpirontGU7dxoDee5XERFdjTy31etyHZ9MhhPdryIioqx7QFWZ46kzIiLKZKdxsyfXoM7OXHtiYlfdC9qKVq7XF5mezyMiioGNdU+oSHV0uvKcxxER0dpQ94KKTP+uiYiIMtGTaCvRlogouvM8oRdlETG49eslO7sAAACA8U6sAAAAAFIRKwAAAIBUxAoAAAAgld8oVqxZsya+9rWvxRVXXBGvvPJKREQ88sgj8Ytf/GK7jgMAAADGnxH/NJDHHnssFixYEL29vfH000/H+eefH1OnTo3bbrstnnnmmfjmN785GjsBAACAcWLEn1mxePHiOOecc+LJJ5+MCRMmtC8/+eST47777tuu4wAAAIDxZ8Sx4p//+Z/jE5/4xFsuf9e73hWrV6/eLqMAAACA8WvEsaKnpyf6+vrecvkTTzwRe+yxx3YZBQAAAIxfI44Vv/M7vxN/9md/FkNDQxERURRFPPPMM3H55ZfH6aefvt0HAgAAAOPLiGPFf/pP/yn6+/tj2rRpsWHDhjj66KNjv/32i8mTJ8df/MVfjMZGAAAAYBwZ8U8D6e3tjR/84Adx//33x2OPPRb9/f3x/ve/PxYsWDAa+wAAAIBxZsSx4g0f/vCH48Mf/vD23AIAAADwm8WKf/7nf4577rknXnzxxWi1WpWPffGLX9wuwwAAAIDxacSx4vOf/3x89rOfjf333z+mT58eRVG0P/arvwYAAAD4TYw4VnzpS1+Kr3/963HOOeeMwhwAAABgvBvxTwPp6OiII488cjS2AAAAAIw8VlxyySVx4403jsYWAAAAgJF/Gcill14aH/3oR2PfffeN9773vdHV1VX5+G233bbdxgEAAADjz4hjxac+9am455574thjj43ddtvNN9UEAAAAtqsRx4pvfOMb8b3vfS8++tGPjsYeAAAAYJwb8fesmDp1auy7776jsQUAAABg5LHi6quvjquuuirWr18/GnsAAACAcW7EXwby5S9/OVasWBHTp0+Pvffe+y3fYPORRx7ZbuMAAACA8WfEseLUU08dhRkAAAAArxtxrLjqqqtGYwcAAABARPwG37MCAAAAYDRt02dWTJ06NZ544onYfffdY9ddd42iKDZ73VdeeWW7jQMAAADGn22KFddff31Mnjy5/estxQoAAACAt2ObYsXZZ5/d/vU555wzWlsAAAAARv49KxqNRrz44otvufzll1+ORqOxXUYBAAAA49eIY0VZlpu8fHBwMLq7u9/2IAAAAGB82+YfXfrlL385IiKKooivfe1rMWnSpPbHms1m3HfffXHAAQds/4UAAADAuLLNseL666+PiNc/s+Kmm26qfMlHd3d37L333nHTTTdt/4UAAADAuLLNsWLlypUREXHsscfGbbfdFrvuuuuojRoNrZ6uaDW66p6RT7Kf7FJs+quMalM2ch2faNU94E1lsnOno5no4EREMTzir7IbVY2Nee5crW1+5hkbZa5TOd3jYCbpHpNz3c2j1ZXre4cVRZ7H5aIj141VbObLqmvTmevciVaecyeazboXVHXkehxMdDePiIhU96xkt1U5NFz3hLay3LYtI37JeM8994x4DAAAAMC2ypWZAQAAgHFPrAAAAABSESsAAACAVMQKAAAAIJXfKFb80z/9U3z84x+P+fPnxy9+8YuIiPhv/+2/xf33379dxwEAAADjz4hjxfe+97044YQTYuLEifGTn/wkBgcHIyJi7dq18fnPf367DwQAAADGlxHHimuvvTZuuumm+C//5b9EV1dX+/IjjzwyHnnkke06DgAAABh/Rhwrli9fHkcdddRbLu/t7Y01a9Zsj00AAADAODbiWDFjxox46qmn3nL5/fffH/vss892GQUAAACMXyOOFeeff35cdNFFsXTp0iiKIp5//vn41re+FZdeeml88pOfHI2NAAAAwDjSOdLf8OlPfzparVYcd9xxsX79+jjqqKOip6cnLr300rjwwgtHYyMAAAAwjow4VhRFEZ/5zGfiT//0T+Opp56K/v7+eO973xuTJk0ajX0AAADAODPiWPGG7u7ueO9737s9twAAAACMPFYce+yxURTFZj/+wx/+8G0NAgAAAMa3EceKww47rPL+0NBQLFu2LB5//PE4++yzt9cuAAAAYJwacay4/vrrN3n51VdfHf39/W97EAAAADC+jfhHl27Oxz/+8fj617++vf44AAAAYJzabrHigQceiAkTJmyvPw4AAAAYp0b8ZSCnnXZa5f2yLGPVqlXx0EMPxZVXXrndhgEAAADj04hjRW9vb+X9jo6O2H///ePP/uzP4vjjj99uwwAAAIDxaUSxotlsxrnnnhsHH3xw7LrrrqO1CQAAABjHRvQ9KxqNRhx//PGxZs2aUZoDAAAAjHcj/gabBx10UPzLv/zLaGwBAAAAGHmsuPbaa+PSSy+NO++8M1atWhV9fX2VNwAAAIC3Y8TfYPPkk0+OiIjf+Z3fiaIo2peXZRlFUUSz2dx+6wAAAIBxZ8Sx4p577hmNHQAAAAAR8RvEirlz58bs2bMrn1UR8fpnVjz77LPbbRgAAAAwPo34e1bMnTs3Xnrppbdc/sorr8TcuXO3yygAAABg/BpxrHjje1P8uv7+/pgwYcJ2GQUAAACMX9v8ZSCLFy+OiIiiKOLKK6+MnXbaqf2xZrMZS5cujcMOO2y7DwQAAADGl22OFT/5yU8i4vXPrPg//+f/RHd3d/tj3d3dceihh8all166/RcCAAAA48o2x4o3fgrIueeeG1/60pdiypQpozYKAAAAGL9G/NNAbr755tHYAQAAABARv8E32AQAAAAYTSP+zIp3qtaERrQ6c/zndgy16p6QVyvXsSmaZd0TqhIdn8aGXMemVea4f7+h0UjWghPNKbMdm7f+gKtadQ7kum8N9+Q5QNkek5vduc7lzs48t1VERHTkOT5lT6PuCVUdPXUvqNqwse4FFUWZ6L7e3VX3gqpMxyYiotmse0FFkWlPstuqSPJv4YiIjrIjYt02XG/0pwAAAABsO7ECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAglVpjxX333RennHJKzJo1K4qiiDvuuKPy8XPOOSeKoqi8nXjiifWMBQAAAMZErbFi3bp1ceihh8aNN9642euceOKJsWrVqvbbt7/97TFcCAAAAIy1zjr/8pNOOilOOumkLV6np6cnZsyYMUaLAAAAgLql/54V9957b0ybNi3233//+OQnPxkvv/zyFq8/ODgYfX19lTcAAADgnSN1rDjxxBPjm9/8Ztx9993xH/7Df4glS5bESSedFM1mc7O/57rrrove3t722+zZs8dwMQAAAPB21fplIFvzsY99rP3rgw8+OA455JDYd9994957743jjjtuk7/niiuuiMWLF7ff7+vrEywAAADgHST1Z1b8un322Sd23333eOqppzZ7nZ6enpgyZUrlDQAAAHjneEfFiueeey5efvnlmDlzZt1TAAAAgFFS65eB9Pf3Vz5LYuXKlbFs2bKYOnVqTJ06Na655po4/fTTY8aMGbFixYq47LLLYr/99osTTjihxtUAAADAaKo1Vjz00ENx7LHHtt9/43tNnH322fHVr341HnvssfjGN74Ra9asiVmzZsXxxx8ff/7nfx49PT11TQYAAABGWa2x4phjjomyLDf78f/1v/7XGK4BAAAAMnhHfc8KAAAAYMcnVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKp11Dxgz5b++UVXmOijFULPuCRXFcKvuCVWJbq+yUfeCX9PKc2wiIqKoe0BeRbK7Vcdw3Qt+TbJTuTFU94I3lcn+F0u2c3l4Qq4H5sbGPAeooyvXsenYkOiOFRHRleufBGWi1zvFcLLXpt3ddU+oGkp2Lic6PuXGjXVPqCiH87zgKcttu18le9oHAAAAxjuxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABS6ax7wLhUlnUvaCuGW3VPqCiGmnVPqGrmOj5RFHUveFOj7gFVRSvP/SoioiPZfascrHtBXmVHrm7f0cx1LjcTHZ7OgVzHptWV6DE58j0OtjrzHJ+yI8+WiIjmpJ66J1Q0Xhuoe0JVI8+LjNaEXLdVx6t9dU+o6uqqe0HV8HDdC9qKZOdO0bFT3RPaOlobI7bhtWmilyAAAAAAYgUAAACQjFgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQSmfdA8ZK2SiibBR1z3jdUN0D3lQMNeueUDWcbE82nY26F7wpW+oskty/39As615QUSQ6dSLXoUm3pzGQa1DZkee+VQznOjap7lcRUeQ6PNHK8rorIlrduW6sbE+hZVeu4xOtugck1tNd94KqVrIbq5HoXG7lelAuBwbqntBWlhu36XrZHisBAACAcU6sAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUOuseMFaK8vW3DMqOou4JeZVJbqSsMh2fVt0DqorBobonVDRauQ5Q0cpz7hStRt0TKoYn5ur2jYFm3RMqmj15nrO6+ofrnlAxMKG77gkVnetyHZ+hSXleZhaZnj8jYv3MCXVPqJj0dK5zJ9NrjOFdeuqeUNHzVH/dE6q68tzPIyJiOM+5XHR11T2hasYedS9oK5qDEa9t/Xq5XqEBAAAA455YAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpdNY9YKwM7NYdnV3ddc+IiIjGYFn3hLbujqLuCRWNVqvuCRXF4FDdEyrKoeG6J7QVRa5zp2w06p5Q1ZGrBTcnJnq4z/MQGBER3a81655Q0erMdd/q7svzuNPqznW/mvDKxronVJTJzp3OgTz3rVZnrnNnp+cH6p5QMTBtp7onVHT153nc6XplQ90TKsopO9c9oaLoW1f3hIpiyuS6J7SVfa/VPaGi9cwv6p7Q1iq37d9YuR65AQAAgHFPrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIJVaY8V1110XRxxxREyePDmmTZsWp556aixfvrxynYGBgVi4cGHstttuMWnSpDj99NPjhRdeqGkxAAAAMNpqjRVLliyJhQsXxoMPPhg/+MEPYmhoKI4//vhYt25d+zqXXHJJ/P3f/3387d/+bSxZsiSef/75OO2002pcDQAAAIymzjr/8rvuuqvy/i233BLTpk2Lhx9+OI466qhYu3Zt/Nf/+l/j1ltvjd/+7d+OiIibb745DjzwwHjwwQfjQx/6UB2zAQAAgFGU6ntWrF27NiIipk6dGhERDz/8cAwNDcWCBQva1znggANir732igceeKCWjQAAAMDoqvUzK35Vq9WKiy++OI488sg46KCDIiJi9erV0d3dHbvsskvlutOnT4/Vq1dv8s8ZHByMwcHB9vt9fX2jthkAAADY/tJ8ZsXChQvj8ccfj+985ztv68+57rrrore3t/02e/bs7bQQAAAAGAspYsWiRYvizjvvjHvuuSf23HPP9uUzZsyIjRs3xpo1ayrXf+GFF2LGjBmb/LOuuOKKWLt2bfvt2WefHc3pAAAAwHZWa6woyzIWLVoUt99+e/zwhz+MuXPnVj5++OGHR1dXV9x9993ty5YvXx7PPPNMzJ8/f5N/Zk9PT0yZMqXyBgAAALxz1Po9KxYuXBi33nprfP/734/Jkye3vw9Fb29vTJw4MXp7e+O8886LxYsXx9SpU2PKlClx4YUXxvz58/0kEAAAANhB1RorvvrVr0ZExDHHHFO5/Oabb45zzjknIiKuv/766OjoiNNPPz0GBwfjhBNOiK985StjvBQAAAAYK7XGirIst3qdCRMmxI033hg33njjGCwCAAAA6pbiG2wCAAAAvEGsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAgFbECAAAASEWsAAAAAFIRKwAAAIBUxAoAAAAglc66B4xHxXCr7gltHQPDdU+oauY5NikVRd0L3pTstioyHZuIiMh132oMNOqe0FYmu6mGu3I9FXaty3XuNCfmOT4dg826J1QM75zn2EREdCY7d1o9eR53oqx7QNXGqd11T6joXJfrvpXJ8C4T6p5Q0fnqhronVE2aWPeCimLNa3VPaCu6uuqeUNExaee6J7R1tDZGDG7D9UZ/CgAAAMC2EysAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVDrrHjBWWl1FtLqKumdERERR1r3gTcW6gbonVBQbBuuekFrR3VX3hLayyHF/ahtu1r2gItnRiUZ/nvtWc6fuuidUlI1ct1bHcKvuCRWZ7llFM9ETaESUHcnOnY2Zbq1cGoO5js366TvXPaGi+5WNdU+oKMo89/X1sybWPaFi8kvr6p5QUXY16p5QUTQS7elMtCUihvedWfeEtuHhgYgHt349n1kBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqYgUAAACQilgBAAAApCJWAAAAAKmIFQAAAEAqnXUPGCsbduuIRneONtMYbNQ9oW3CcLPuCVXDw3UvqCqKuhdUbax7wJuKsqx7QkXZmed+9bpce8qOHI9/ERGR7G7VGGjVPaGi1Z3r3Olo5jk+rQm5jk1jY55jExExPKm77gkVRSvP80Qz2fP55J9vqHtCxdCUrronVHSuy/N6cOKqgbonVAxN3anuCRWNDUN1T6hoTZ1c94S2jtdy3c87n3q+7glvam3bP2oSvXoFAAAAECsAAACAZMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIJXOugeMlZ1eakVnV6vuGRER0fPLwbonvGkg0ZaIKIeH656QWlGWdU94U1HUvaAi15qIMtnx6ch07jRyHZvhSV11T6gaSHRbRcTwxDzHpzHYrHtCxfDERt0TKrr6cz2HNifkOT6d63Idm427dNc9oWLCSwN1T6hodeX5/6kD03rqnlAx4cVcr92biZ4jIiI6h/I8T5Q7Tah7QkVrt8l1T2hrNQciXtz69fI8EgAAAACEWAEAAAAkI1YAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqXTWPWCsDE8souwu6p4RERFdO+c57J2Tdqp7QkWxfqDuCRXl8HDdE/Iqy7oXMAJlR6I2nezc6X51sO4JFc2d8jxHREQ0Bpt1T2hr9jTqnlDRuT7PsYmI2Nib69zp2Jjnvt7qSvQYGBETXlhf94SK5qTuuidUdGzI8/prwou5niOGpnTVPaGiqz/PbRURMbzrxLontHW9tK7uCRUdK5+ve0JbR2vjtl1vlHcAAAAAjIhYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKmIFAAAAkIpYAQAAAKQiVgAAAACpiBUAAABAKrXGiuuuuy6OOOKImDx5ckybNi1OPfXUWL58eeU6xxxzTBRFUXn7oz/6o5oWAwAAAKOt1lixZMmSWLhwYTz44IPxgx/8IIaGhuL444+PdevWVa53/vnnx6pVq9pvX/jCF2paDAAAAIy2zjr/8rvuuqvy/i233BLTpk2Lhx9+OI466qj25TvttFPMmDFjrOcBAAAANUj1PSvWrl0bERFTp06tXP6tb30rdt999zjooIPiiiuuiPXr12/2zxgcHIy+vr7KGwAAAPDOUetnVvyqVqsVF198cRx55JFx0EEHtS///d///ZgzZ07MmjUrHnvssbj88stj+fLlcdttt23yz7nuuuvimmuuGavZAAAAwHaWJlYsXLgwHn/88bj//vsrl19wwQXtXx988MExc+bMOO6442LFihWx7777vuXPueKKK2Lx4sXt9/v6+mL27NmjNxwAAADYrlLEikWLFsWdd94Z9913X+y5555bvO68efMiIuKpp57aZKzo6emJnp6eUdkJAAAAjL5aY0VZlnHhhRfG7bffHvfee2/MnTt3q79n2bJlERExc+bMUV4HAAAA1KHWWLFw4cK49dZb4/vf/35Mnjw5Vq9eHRERvb29MXHixFixYkXceuutcfLJJ8duu+0Wjz32WFxyySVx1FFHxSGHHFLndAAAAGCU1BorvvrVr0ZExDHHHFO5/Oabb45zzjknuru74x//8R/jhhtuiHXr1sXs2bPj9NNPj89+9rM1rAUAAADGQu1fBrIls2fPjiVLlozRGgAAACCDjroHAAAAAPwqsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUhErAAAAgFTECgAAACAVsQIAAABIRawAAAAAUumse8BY6dxQRudwWfeMiIhoDLXqnvCmrmSnQLI9RUeyntdR1L2grUx2W0W22yqZYuNQ3RPelOymGu7trntCRWNdotsqIgb2mFj3hLau9cN1T6jYOKWr7gkVXa81655Q0ZyQ587eGMx1bAZm7FT3hIqdnu6re0JFc3JP3RPaBqbl2RIRsfPTr9U9oaI5Kdfx6frFmrontJU75To2xYw96p7wpuZgxCtbv1qeZxEAAACAECsAAACAZMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIJXOugeMlaJVRtEq657xumaSHRFRdibrVV25TsmiubHuCVVFUfeCtqLZqntC1dBw3QuqsjzevKFMtKeV7NyZMqHuBRUdfRvqnlDRmNJT94S2zrWDdU+oGNylq+4JFV19uY5Pq3ti3RPaOtYP1T2h4rWD8hybiIidH1tX94SKYkKe14N9c/JsiYiY/Nj6uidUNCfleY6IiIi1/XUvaCu6cz1H/OKE3eue0NYcHIj46davl+xfqgAAAMB4J1YAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqXTWPWC0lWUZERHDQwM1L3lTY3hj3RPe1Byse0FF0Ux0bCKiaOXaE0WivlgUdS+oarXqXlDVKuteUFXm2VM282yJiBgezvP8EBFRJHtcznR8OrIdm6FcL6OGm3luq4iI4aE8zxONZOdOc2Oy26qV7Pgkur2ag26rLcn0HBEREZleuyc6jyNynctvPAaWW3l9WpRbu8Y73HPPPRezZ8+uewYAAADwr5599tnYc889N/vxHT5WtFqteP7552Py5MlRvI3/E9zX1xezZ8+OZ599NqZMmbIdF8LYci6zI3Aes6NwLrOjcC6zo3Auj76yLOO1116LWbNmRUfH5j9zPNfnL46Cjo6OLdaakZoyZYqTlh2Cc5kdgfOYHYVzmR2Fc5kdhXN5dPX29m71Oom+AB4AAABArAAAAACSESu2UU9PT1x11VXR09NT9xR4W5zL7Aicx+wonMvsKJzL7Cicy3ns8N9gEwAAAHhn8ZkVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgV2+DGG2+MvffeOyZMmBDz5s2LH//4x3VPghG5+uqroyiKytsBBxxQ9yzYqvvuuy9OOeWUmDVrVhRFEXfccUfl42VZxuc+97mYOXNmTJw4MRYsWBBPPvlkPWNhC7Z2Lp9zzjlveZw+8cQT6xkLm3HdddfFEUccEZMnT45p06bFqaeeGsuXL69cZ2BgIBYuXBi77bZbTJo0KU4//fR44YUXaloMm7Yt5/IxxxzzlsflP/qjP6pp8fgkVmzFd7/73Vi8eHFcddVV8cgjj8Shhx4aJ5xwQrz44ot1T4MRed/73herVq1qv91///11T4KtWrduXRx66KFx4403bvLjX/jCF+LLX/5y3HTTTbF06dLYeeed44QTToiBgYExXgpbtrVzOSLixBNPrDxOf/vb3x7DhbB1S5YsiYULF8aDDz4YP/jBD2JoaCiOP/74WLduXfs6l1xySfz93/99/O3f/m0sWbIknn/++TjttNNqXA1vtS3nckTE+eefX3lc/sIXvlDT4vHJjy7dinnz5sURRxwR//k//+eIiGi1WjF79uy48MIL49Of/nTN62DbXH311XHHHXfEsmXL6p4Cv7GiKOL222+PU089NSJe/6yKWbNmxZ/8yZ/EpZdeGhERa9eujenTp8ctt9wSH/vYx2pcC5v36+dyxOufWbFmzZq3fMYFZPbSSy/FtGnTYsmSJXHUUUfF2rVrY4899ohbb701/s2/+TcREfGzn/0sDjzwwHjggQfiQx/6UM2LYdN+/VyOeP0zKw477LC44YYb6h03jvnMii3YuHFjPPzww7FgwYL2ZR0dHbFgwYJ44IEHalwGI/fkk0/GrFmzYp999ok/+IM/iGeeeabuSfC2rFy5MlavXl15jO7t7Y158+Z5jOYd6d57741p06bF/vvvH5/85Cfj5ZdfrnsSbNHatWsjImLq1KkREfHwww/H0NBQ5XH5gAMOiL322svjMqn9+rn8hm9961ux++67x0EHHRRXXHFFrF+/vo5541Zn3QMy++UvfxnNZjOmT59euXz69Onxs5/9rKZVMHLz5s2LW265Jfbff/9YtWpVXHPNNfGRj3wkHn/88Zg8eXLd8+A3snr16oiITT5Gv/ExeKc48cQT47TTTou5c+fGihUr4t//+38fJ510UjzwwAPRaDTqngdv0Wq14uKLL44jjzwyDjrooIh4/XG5u7s7dtlll8p1PS6T2abO5YiI3//93485c+bErFmz4rHHHovLL788li9fHrfddluNa8cXsQLGgZNOOqn960MOOSTmzZsXc+bMib/5m7+J8847r8ZlAERE5cuWDj744DjkkENi3333jXvvvTeOO+64GpfBpi1cuDAef/xx3wOLd7zNncsXXHBB+9cHH3xwzJw5M4477rhYsWJF7LvvvmM9c1zyZSBbsPvuu0ej0XjLdzB+4YUXYsaMGTWtgrdvl112ife85z3x1FNP1T0FfmNvPA57jGZHtM8++8Tuu+/ucZqUFi1aFHfeeWfcc889seeee7YvnzFjRmzcuDHWrFlTub7HZbLa3Lm8KfPmzYuI8Lg8hsSKLeju7o7DDz887r777vZlrVYr7r777pg/f36Ny+Dt6e/vjxUrVsTMmTPrngK/sblz58aMGTMqj9F9fX2xdOlSj9G84z333HPx8ssve5wmlbIsY9GiRXH77bfHD3/4w5g7d27l44cffnh0dXVVHpeXL18ezzzzjMdlUtnaubwpb3yjeo/LY8eXgWzF4sWL4+yzz44PfOAD8cEPfjBuuOGGWLduXZx77rl1T4Ntdumll8Ypp5wSc+bMieeffz6uuuqqaDQaceaZZ9Y9Dbaov7+/8n8wVq5cGcuWLYupU6fGXnvtFRdffHFce+218e53vzvmzp0bV155ZcyaNavyUxYggy2dy1OnTo1rrrkmTj/99JgxY0asWLEiLrvssthvv/3ihBNOqHE1VC1cuDBuvfXW+P73vx+TJ09ufx+K3t7emDhxYvT29sZ5550XixcvjqlTp8aUKVPiwgsvjPnz5/tJIKSytXN5xYoVceutt8bJJ58cu+22Wzz22GNxySWXxFFHHRWHHHJIzevHkZKt+qu/+qtyr732Kru7u8sPfvCD5YMPPlj3JBiRM844o5w5c2bZ3d1dvutd7yrPOOOM8qmnnqp7FmzVPffcU0bEW97OPvvssizLstVqlVdeeWU5ffr0sqenpzzuuOPK5cuX1zsaNmFL5/L69evL448/vtxjjz3Krq6ucs6cOeX5559frl69uu7ZULGpczgiyptvvrl9nQ0bNpR//Md/XO66667lTjvtVP7e7/1euWrVqvpGwyZs7Vx+5plnyqOOOqqcOnVq2dPTU+63337ln/7pn5Zr166td/g4U5RlWY5lHAEAAADYEt+zAgAAAEhFrAAAAABSESsAAACAVMQKAAAAIBWxAgAAAEhFrAAAAABSESsAAACAVMQKAGC72nvvveOGG25ov18URdxxxx1jvuPqq6+Oww47bMS/75hjjomLL754u+8BALZdZ90DAIAd26pVq2LXXXfdputeffXVcccdd8SyZctGd1RE3HvvvXHsscfGq6++Grvsskv78ttuuy26urpG/e8HADZPrAAA3mLjxo3R3d29Xf6sGTNmbJc/Z6xMnTq17gkAMO75MhAA2MEdc8wxsWjRoli0aFH09vbG7rvvHldeeWWUZdm+zt577x1//ud/HmeddVZMmTIlLrjggoiIuP/+++MjH/lITJw4MWbPnh2f+tSnYt26de3f9+KLL8Ypp5wSEydOjLlz58a3vvWtt/z9v/5lIM8991yceeaZMXXq1Nh5553jAx/4QCxdujRuueWWuOaaa+LRRx+NoiiiKIq45ZZbIiJizZo18Yd/+Iexxx57xJQpU+K3f/u349FHH638PX/5l38Z06dPj8mTJ8d5550XAwMDmz0mTz/9dBx77LEREbHrrrtGURRxzjnntI/Xr34ZyN577x3XXnttnHXWWTFp0qSYM2dO/N3f/V289NJL8bu/+7sxadKkOOSQQ+Khhx6q/B1bO3YAwOaJFQAwDnzjG9+Izs7O+PGPfxxf+tKX4otf/GJ87Wtfq1znP/7H/xiHHnpo/OQnP4krr7wyVqxYESeeeGKcfvrp8dhjj8V3v/vduP/++2PRokXt33POOefEs88+G/fcc0/89//+3+MrX/lKvPjii5vd0d/fH0cffXT84he/iL/7u7+LRx99NC677LJotVpxxhlnxJ/8yZ/E+973vli1alWsWrUqzjjjjIiI+Lf/9t/Giy++GP/zf/7PePjhh+P9739/HHfccfHKK69ERMTf/M3fxNVXXx2f//zn46GHHoqZM2fGV77ylc3umD17dnzve9+LiIjly5fHqlWr4ktf+tJmr3/99dfHkUceGT/5yU/iox/9aPy7f/fv4qyzzoqPf/zj8cgjj8S+++4bZ511VjsAbcuxAwC2oAQAdmhHH310eeCBB5atVqt92eWXX14eeOCB7ffnzJlTnnrqqZXfd95555UXXHBB5bJ/+qd/Kjs6OsoNGzaUy5cvLyOi/PGPf9z++E9/+tMyIsrrr7++fVlElLfffntZlmX513/91+XkyZPLl19+eZNbr7rqqvLQQw99y985ZcqUcmBgoHL5vvvuW/71X/91WZZlOX/+/PKP//iPKx+fN2/eW/6sX3XPPfeUEVG++uqrlcuPPvro8qKLLmq/P2fOnPLjH/94+/1Vq1aVEVFeeeWV7cseeOCBMiLKVatWlWW59WMHAGyZz6wAgHHgQx/6UBRF0X5//vz58eSTT0az2Wxf9oEPfKDyex599NG45ZZbYtKkSe23E044IVqtVqxcuTJ++tOfRmdnZxx++OHt33PAAQdUvlnlr1u2bFn81m/91oi+L8Sjjz4a/f39sdtuu1W2rFy5MlasWBERET/96U9j3rx5ld83f/78bf47tuaQQw5p/3r69OkREXHwwQe/5bI3Pqtka8cOANgy32ATAIiIiJ133rnyfn9/f3ziE5+IT33qU2+57l577RVPPPHEiP+OiRMnjvj39Pf3x8yZM+Pee+99y8e2FEa2p1/96SBvRJ9NXdZqtSJi68cOANgysQIAxoGlS5dW3n/wwQfj3e9+dzQajc3+nve///3x//7f/4v99ttvkx8/4IADYnh4OB5++OE44ogjIuL17/+wZs2azf6ZhxxySHzta1+LV155ZZOfXdHd3V35bI83dqxevTo6Oztj77333uSfe+CBB8bSpUvjrLPOqvw3bskbP+3k1/++7WFrxw4A2DJfBgIA48AzzzwTixcvjuXLl8e3v/3t+Ku/+qu46KKLtvh7Lr/88vjRj34UixYtimXLlsWTTz4Z3//+99vfJHL//fePE088MT7xiU/E0qVL4+GHH44//MM/3OJnT5x55pkxY8aMOPXUU+N//+//Hf/yL/8S3/ve9+KBBx6IiNd/8sbKlStj2bJl8ctf/jIGBwdjwYIFMX/+/Dj11FPjH/7hH+Lpp5+OH/3oR/GZz3ym/RM4Lrroovj6178eN998czzxxBNx1VVXxf/9v/93i/99c+bMiaIo4s4774yXXnop+vv7R3JIt2hrxw4A2DKxAgDGgbPOOis2bNgQH/zgB2PhwoVx0UUXtX886eYccsghsWTJknjiiSfiIx/5SPzWb/1WfO5zn4tZs2a1r3PzzTfHrFmz4uijj47TTjstLrjggpg2bdpm/8zu7u74h3/4h5g2bVqcfPLJcfDBB8df/uVftj/D4/TTT48TTzwxjj322Nhjjz3i29/+dhRFEf/jf/yPOOqoo+Lcc8+N97znPfGxj30sfv7zn7e/V8QZZ5wRV155ZVx22WVx+OGHx89//vP45Cc/ucX/vne9611xzTXXxKc//emYPn36dg0J23LsAIDNK8ryV37IOgCwwznmmGPisMMOixtuuKHuKQAA28RnVgAAAACpiBUAAABAKr4MBAAAAEjFZ1YAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqYgVAAAAQCpiBQAAAJCKWAEAAACkIlYAAAAAqfz/+iyI9ymP+VUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1300x1300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(13, 13))\n",
    "ax.imshow(outputs, cmap='viridis')\n",
    "ax.set_xlabel('predicted time')\n",
    "ax.set_ylabel('true time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "0fa287de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.21387405, 0.17206091, 0.1111422 , 0.06354874, 0.04119461],\n",
       "       [0.17045611, 0.15384128, 0.11962737, 0.07720544, 0.05201893],\n",
       "       [0.10863931, 0.11159411, 0.12346836, 0.10238995, 0.07511837],\n",
       "       [0.06167122, 0.06953499, 0.10203358, 0.12670563, 0.1141167 ],\n",
       "       [0.04795824, 0.0526353 , 0.07805046, 0.11621552, 0.13023582]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "outputs[:5, :5]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "59d44e8f",
   "metadata": {},
   "source": [
    "## scan parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "a8378946",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning_rate=5.000e-06:loss=1.297e-01\n",
      "learning_rate=1.581e-05:loss=1.295e-01\n",
      "learning_rate=5.000e-05:loss=1.297e-01\n",
      "learning_rate=1.581e-04:loss=1.307e-01\n",
      "learning_rate=5.000e-04:loss=1.340e-01\n"
     ]
    }
   ],
   "source": [
    "N_scan = 5\n",
    "for learning_rate_ in learning_rate*np.logspace(-1, 1, N_scan, base=10):\n",
    "    logistic_model, loss = fit_data(fr, learning_rate=learning_rate_, verbose=False)\n",
    "    print(f'learning_rate={learning_rate_:.3e}:{loss=:.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "a2b59fd2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logit_0=-2.000e-01:loss=1.297e-01\n",
      "logit_0=-6.325e-01:loss=1.297e-01\n",
      "logit_0=-2.000e+00:loss=1.297e-01\n",
      "logit_0=-6.325e+00:loss=1.297e-01\n",
      "logit_0=-2.000e+01:loss=1.297e-01\n"
     ]
    }
   ],
   "source": [
    "for logit_0_ in logit_0*np.logspace(-1, 1, N_scan, base=10):\n",
    "    logistic_model, loss = fit_data(fr, logit_0=logit_0_, verbose=False)\n",
    "    print(f'logit_0={logit_0_:.3e}:{loss=:.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "ad355197",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logit_wt=1.000e-03:loss=1.297e-01\n",
      "logit_wt=3.162e-03:loss=1.297e-01\n",
      "logit_wt=1.000e-02:loss=1.297e-01\n",
      "logit_wt=3.162e-02:loss=1.297e-01\n",
      "logit_wt=1.000e-01:loss=1.297e-01\n"
     ]
    }
   ],
   "source": [
    "for logit_wt_ in logit_wt*np.logspace(-1, 1, N_scan, base=10):\n",
    "    logistic_model, loss = fit_data(fr, logit_wt=logit_wt_, verbose=False)\n",
    "    print(f'logit_wt={logit_wt_:.3e}:{loss=:.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "ab94a08c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta1=9.900e-01:loss=1.298e-01\n",
      "beta1=9.982e-01:loss=1.295e-01\n",
      "beta1=9.997e-01:loss=1.293e-01\n",
      "beta1=9.999e-01:loss=1.532e-01\n",
      "beta1=1.000e+00:loss=1.532e-01\n"
     ]
    }
   ],
   "source": [
    "for beta1_ in 1-np.logspace(-2, -5, N_scan, base=10):\n",
    "    logistic_model, loss = fit_data(fr, beta1=beta1_, verbose=False)\n",
    "    print(f'beta1={beta1_:.3e}:{loss=:.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ce111043",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta2=9.000e-01:loss=2.284e-01\n",
      "beta2=9.822e-01:loss=1.298e-01\n",
      "beta2=9.968e-01:loss=1.297e-01\n",
      "beta2=9.994e-01:loss=1.297e-01\n",
      "beta2=9.999e-01:loss=1.297e-01\n"
     ]
    }
   ],
   "source": [
    "for beta2_ in 1-np.logspace(-1, -4, N_scan, base=10):\n",
    "    logistic_model, loss = fit_data(fr, beta2=beta2_, verbose=False)\n",
    "    print(f'beta2={beta2_:.3e}:{loss=:.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "c796488e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "amsgrad=True:loss=1.296e-01\n"
     ]
    }
   ],
   "source": [
    "logistic_model, loss = fit_data(fr, amsgrad=not(amsgrad), verbose=False)\n",
    "print(f'amsgrad={not(amsgrad)}:{loss=:.3e}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}