AEP metrics new staff

analysis/AEPs/Performance prediction.ipynb

@@ -0,0 +1,235 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "7dbf62e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys, os\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
+    "\n",
+    "import numpy as np\n",
+    "import h5py\n",
+    "import json\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "from scipy import signal\n",
+    "from scipy import stats\n",
+    "from target import build_tgt_matrix\n",
+    "import pandas as pd\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "220a45cd",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/javascript": [
+       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
+       "    return false;\n",
+       "}\n"
+      ],
+      "text/plain": [
+       "<IPython.core.display.Javascript object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "%%javascript\n",
+    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
+    "    return false;\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "id": "06dca5bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "source = '/home/sobolev/nevermind/Andrey/data'\n",
+    "report = '/home/sobolev/nevermind/Andrey/analysis/PSTH'\n",
+    "\n",
+    "selected_sessions = [\n",
+    "'009266_hippoSIT_2023-04-17_17-04-17',  # ch17, 20 + 55 correction, 5067 events. Showcase for N2 / N3 mod in target\n",
+    "'009266_hippoSIT_2023-04-18_10-10-37',  # ch17, 10 + 55 correction, 5682 events\n",
+    "'009266_hippoSIT_2023-04-18_17-03-10',  # ch17, 6 + 55 correction, 5494 events: FIXME very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-04-19_10-33-51',  # ch17, 4 + 55 correction, 6424 events: very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-04-20_08-57-39',  # ch1, 1 + 55 correction, 6424 events. Showcase for N2 / N3 mod in target\n",
+    "'009266_hippoSIT_2023-04-24_16-56-55',  # ch17, 5 + 55* correction, 6165 events, frequency\n",
+    "'009266_hippoSIT_2023-04-26_08-20-17',  # ch17, 12 + 55* correction, 6095 events, duration - showcase for N2 \n",
+    "'009266_hippoSIT_2023-05-02_12-22-14',  # ch20, 10 + 55 correction, 5976 events, FIXME very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-05-04_09-11-06',  # ch17, 5 + 55* correction, 4487 events, coma session with baseline AEPs\n",
+    "'009266_hippoSIT_2023-05-04_19-47-15',  # ch20, 2 + 55 correction, 5678 events, duration\n",
+    "]\n",
+    "\n",
+    "session = selected_sessions[0]\n",
+    "\n",
+    "animal      = session.split('_')[0]\n",
+    "sessionpath = os.path.join(source, animal, session)\n",
+    "aeps_file   = os.path.join(sessionpath, 'AEPs.h5')\n",
+    "h5name      = os.path.join(sessionpath, session + '.h5')\n",
+    "report_path = os.path.join(report, session)\n",
+    "if not os.path.exists(report_path):\n",
+    "    os.makedirs(report_path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 32,
+   "id": "b5a09f6a",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "((5067, 200), (73, 5))"
+      ]
+     },
+     "execution_count": 32,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "with h5py.File(h5name, 'r') as f:\n",
+    "    tl = np.array(f['processed']['timeline'])  # time, X, Y, speed, etc.\n",
+    "    trials = np.array(f['processed']['trial_idxs'])  # t_start_idx, t_end_idx, x_tgt, y_tgt, r_tgt, result\n",
+    "    cfg = json.loads(f['processed'].attrs['parameters'])\n",
+    "    \n",
+    "with h5py.File(aeps_file, 'r') as f:\n",
+    "    aeps = np.array(f['aeps'])\n",
+    "    aeps_events = np.array(f['aeps_events'])\n",
+    "    \n",
+    "# TODO find better way. Remove outliers\n",
+    "aeps[aeps > 5000]  =  5000\n",
+    "aeps[aeps < -5000] = -5000\n",
+    "\n",
+    "# load metrics\n",
+    "AEP_metrics_lims = {}\n",
+    "AEP_metrics_raw  = {}\n",
+    "AEP_metrics_norm = {}\n",
+    "with h5py.File(aeps_file, 'r') as f:\n",
+    "    for metric_name in f['raw']:\n",
+    "        AEP_metrics_raw[metric_name]  = np.array(f['raw'][metric_name])\n",
+    "        AEP_metrics_norm[metric_name] = np.array(f['norm'][metric_name])\n",
+    "        AEP_metrics_lims[metric_name] = [int(x) for x in f['raw'][metric_name].attrs['limits'].split(',')]\n",
+    "\n",
+    "tgt_dur = cfg['experiment']['target_duration']\n",
+    "tgt_matrix = build_tgt_matrix(tl, aeps_events, cfg['experiment']['target_duration'])\n",
+    "\n",
+    "aeps.shape, tgt_matrix.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b0e36f45",
+   "metadata": {},
+   "source": [
+    "## Is performance dependent on AEP states?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "id": "9294c349",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# separate states based on AEP metrics - metric mean before entering\n",
+    "def compute_state_idxs(metric_name, n_pulses=10):\n",
+    "    idxs_aeps_high, idxs_aeps_low = [], []  # indices to tgt_matrix\n",
+    "\n",
+    "    metric_mean = AEP_metrics_norm[metric_name].mean()\n",
+    "    for i, tgt_enter in enumerate(tgt_matrix):\n",
+    "        if tgt_enter[2] - n_pulses < 0:\n",
+    "            continue\n",
+    "        metric_inst = AEP_metrics_norm[metric_name][tgt_enter[2] - n_pulses:tgt_enter[2]].mean()\n",
+    "        if metric_inst > metric_mean:\n",
+    "            idxs_aeps_high.append(i)\n",
+    "        else:\n",
+    "            idxs_aeps_low.append(i)\n",
+    "            \n",
+    "    return idxs_aeps_low, idxs_aeps_high"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "id": "02ffa773",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'N1': (0.6164383561643836, 0.3835616438356164),\n",
+       " 'P1': (0.5342465753424658, 0.4657534246575342),\n",
+       " 'P2': (0.6164383561643836, 0.3835616438356164),\n",
+       " 'P3': (0.3835616438356164, 0.6164383561643836)}"
+      ]
+     },
+     "execution_count": 34,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# test if the animal will stay in the island depending on high / low AEP metrics\n",
+    "\n",
+    "metric_name = 'P1'\n",
+    "n_pulses = 8\n",
+    "predictions = {}\n",
+    "\n",
+    "for metric_name in AEP_metrics_norm.keys():\n",
+    "    idxs_aeps_low, idxs_aeps_high = compute_state_idxs(metric_name, n_pulses)\n",
+    "\n",
+    "    actual = tgt_matrix[:, 4]\n",
+    "    predicted = np.zeros(len(tgt_matrix))\n",
+    "    predicted[idxs_aeps_high] = 1\n",
+    "    high_low = (predicted == actual).sum() / len(tgt_matrix)\n",
+    "    \n",
+    "    predicted = np.zeros(len(tgt_matrix))\n",
+    "    predicted[idxs_aeps_low] = 1\n",
+    "    low_high = (predicted == actual).sum() / len(tgt_matrix)\n",
+    "    \n",
+    "    predictions[metric_name] = (high_low, low_high)\n",
+    "    \n",
+    "predictions"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

analysis/target.py

@@ -0,0 +1,46 @@
+import numpy as np
+
+
+def build_tgt_matrix(tl, aeps_events, tgt_dur):
+    # compute timeline / AEP indices of entrances / exist to the target
+    tl_tgt_start_idxs = []  # timeline indices of entrance in target
+    tl_tgt_end_idxs   = []  # timeline indices of exit from target
+
+    for i in range(len(tl) - 1):
+        if tl[i][6] < 2 and tl[i+1][6] == 2:
+            tl_tgt_start_idxs.append(i + 1)
+        if tl[i][6] == 2 and tl[i+1][6] < 2:
+            tl_tgt_end_idxs.append(i)
+
+    # ignore first/last target if not ended
+    if tl_tgt_start_idxs[-1] > tl_tgt_end_idxs[-1]:
+        tl_tgt_start_idxs = tl_tgt_start_idxs[:-1]
+    if tl_tgt_end_idxs[0] < tl_tgt_start_idxs[0]:
+        tl_tgt_end_idxs = tl_tgt_end_idxs[1:]
+    tl_tgt_start_idxs = np.array(tl_tgt_start_idxs)
+    tl_tgt_end_idxs   = np.array(tl_tgt_end_idxs)
+
+    aeps_tgt_start_idxs = []  # indices of first AEPs in target
+    aeps_tgt_end_idxs   = []  # indices of last AEPs in target
+
+    for idx in tl_tgt_start_idxs:
+        aeps_tgt_start_idxs.append(np.abs(aeps_events[:, 0] - tl[idx][0]).argmin())
+    for idx in tl_tgt_end_idxs:
+        aeps_tgt_end_idxs.append(np.abs(aeps_events[:, 0] - tl[idx][0]).argmin())
+    aeps_tgt_start_idxs = np.array(aeps_tgt_start_idxs)
+    aeps_tgt_end_idxs = np.array(aeps_tgt_end_idxs)
+
+    # successful / missed
+    tgt_results = np.zeros(len(tl_tgt_start_idxs))
+    succ_idxs = np.where((tl_tgt_end_idxs - tl_tgt_start_idxs > tgt_dur * 100 - 10) & \
+                         (tl_tgt_end_idxs - tl_tgt_start_idxs < tgt_dur * 100 + 10))[0]
+    tgt_results[succ_idxs] = 1
+
+    # tl_idx_start, tl_idx_end, aep_idx_start, aer_idx_end, success / miss
+    return np.column_stack([
+        tl_tgt_start_idxs,
+        tl_tgt_end_idxs,
+        aeps_tgt_start_idxs,
+        aeps_tgt_end_idxs,
+        tgt_results
+    ]).astype(np.int32)

playground.ipynb

@@ -875,7 +875,9 @@
    "cell_type": "code",
    "execution_count": 79,
    "id": "e93f7f95",
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [
     {
      "data": {
@@ -939,10 +941,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "275103f0",
    "metadata": {},
    "outputs": [],
+   "source": [
+    "import pandas\n",
+    "human_data = pandas.read_csv('https://content.labxchange.org/labs/datasets/CodeSciLab/HumanHeightWeightData.csv')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10d00841",
+   "metadata": {},
+   "outputs": [],
    "source": []
   }
  ],

postprocessing/execute.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 46,
+   "execution_count": 1,
    "id": "c37dce82",
    "metadata": {
     "scrolled": true
@@ -39,7 +39,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 47,
+   "execution_count": 2,
    "id": "a8ea2991",
    "metadata": {},
    "outputs": [
@@ -75,7 +75,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 5,
    "id": "3891b792",
    "metadata": {},
    "outputs": [],
@@ -100,8 +100,16 @@
     "#sessions = [filebase]\n",
     "#sessions = processed_008229\n",
     "sessions = [\n",
+    "'009266_hippoSIT_2023-04-17_17-04-17',  # ch17, 20 + 55 correction, 5067 events. Showcase for N2 / N3 mod in target\n",
     "'009266_hippoSIT_2023-04-18_10-10-37',  # ch17, 10 + 55 correction, 5682 events\n",
     "'009266_hippoSIT_2023-04-18_17-03-10',  # ch17, 6 + 55 correction, 5494 events: FIXME very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-04-19_10-33-51',  # ch17, 4 + 55 correction, 6424 events: very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-04-20_08-57-39',  # ch1, 1 + 55 correction, 6424 events. Showcase for N2 / N3 mod in target\n",
+    "'009266_hippoSIT_2023-04-24_16-56-55',  # ch17, 5 + 55* correction, 6165 events, frequency\n",
+    "'009266_hippoSIT_2023-04-26_08-20-17',  # ch17, 12 + 55* correction, 6095 events, duration - showcase for N2 \n",
+    "'009266_hippoSIT_2023-05-02_12-22-14',  # ch20, 10 + 55 correction, 5976 events, FIXME very weird 1-2nd in target, find out\n",
+    "'009266_hippoSIT_2023-05-04_09-11-06',  # ch17, 5 + 55* correction, 4487 events, coma session with baseline AEPs\n",
+    "'009266_hippoSIT_2023-05-04_19-47-15',\n",
     "]\n",
     "\n",
     "# FIXME move occupancy outside units\n",
@@ -141,7 +149,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 49,
+   "execution_count": 6,
    "id": "24cd2f6e",
    "metadata": {},
    "outputs": [
@@ -149,8 +157,16 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "session 009266_hippoSIT_2023-04-17_17-04-17 done\n",
       "session 009266_hippoSIT_2023-04-18_10-10-37 done\n",
-      "session 009266_hippoSIT_2023-04-18_17-03-10 done\n"
+      "session 009266_hippoSIT_2023-04-18_17-03-10 done\n",
+      "session 009266_hippoSIT_2023-04-19_10-33-51 done\n",
+      "session 009266_hippoSIT_2023-04-20_08-57-39 done\n",
+      "session 009266_hippoSIT_2023-04-24_16-56-55 done\n",
+      "session 009266_hippoSIT_2023-04-26_08-20-17 done\n",
+      "session 009266_hippoSIT_2023-05-02_12-22-14 done\n",
+      "session 009266_hippoSIT_2023-05-04_09-11-06 done\n",
+      "session 009266_hippoSIT_2023-05-04_19-47-15 done\n"
      ]
     }
    ],

postprocessing/pack.py

@@ -180,7 +180,7 @@ def pack(session_path):
             if curr_sound_idx + 1 >= len(sounds):
                 break
 
-            if temp_tl[i][0] > sounds[curr_sound_idx + 1][0]:
+            if temp_tl[i][0] > sounds[curr_sound_idx][0]:
                 curr_sound_idx += 1
             sound_tl[i] = sounds[curr_sound_idx][1]
 

session/adapters.py

@@ -1,6 +1,9 @@
 import numpy as np
 import os
 import h5py
+from scipy import signal
+from scipy.signal import butter, sosfilt
+
 
 COLORS  = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:gray'] 
 EPOCH_NAMES = ('Original', 'Conflict', 'Control', 'All')
@@ -98,4 +101,50 @@ def create_dataset(h5name, where, descriptor, dataset):
             
         ds = target_group.create_dataset(descriptor['name'], data=dataset)
         for i, dim in enumerate(descriptor['dims']):
-            ds.attrs['dim%s' % i] = dim
+            ds.attrs['dim%s' % i] = dim
+            
+            
+class DatProcessor:
+    
+    def __init__(self, dat_file):
+        # read this from XML file
+        self.s_rate = 30000
+        self.ch_no = 64
+        self.dat_file = dat_file
+        
+    @staticmethod
+    def butter_bandpass_filter(data, lowcut, highcut, fs, order=6):
+        def butter_bandpass(lowcut, highcut, fs, order=5):
+            nyq = 0.5 * fs
+            low = lowcut / nyq
+            high = highcut / nyq
+            sos = butter(order, [low, high], analog=False, btype='band', output='sos')
+            return sos        
+
+        sos = butter_bandpass(lowcut, highcut, fs, order=order)
+        y = sosfilt(sos, data)
+        return y
+    
+    def read_block_from_dat(self, duration, offset):
+        """
+        duration      in seconds
+        offset        in seconds
+        """
+        count = self.s_rate * self.ch_no * duration  # number of values to read
+        offset_in_bytes = offset * self.s_rate * self.ch_no * 2  # assuming int16 is 2 bytes
+        block = np.fromfile(self.dat_file, dtype=np.int16, count=int(count), offset=int(offset_in_bytes))
+        return block.reshape([int(self.s_rate * duration), self.ch_no])
+    
+    def get_single_channel(self, channel_no):
+        size = os.path.getsize(self.dat_file)
+        samples_no = size / (64 * 2)
+
+        raw_signal = np.zeros(int(samples_no))  # length in time: samples_no / sample_rate
+        offset = 0
+
+        while offset < samples_no / self.s_rate - 1:
+            block = self.read_block_from_dat(1, offset)  # read in 1 sec blocks
+            raw_signal[self.s_rate*offset:self.s_rate*(offset + 1)] = block[:, channel_no]
+            offset += 1
+
+        return raw_signal