analysis/imports.py

@@ -18,4 +18,7 @@ from session.utils import *
 from session.adapters import H5NAMES, EPOCH_NAMES, COLORS
 from postprocessing.spatial import place_field_2D, map_stats, get_field_patches, get_positions_relative_to
 from postprocessing.spatial import bins2meters, cart2pol, pol2cart
-from postprocessing.spiketrain import instantaneous_rate
+from postprocessing.spiketrain import instantaneous_rate
+
+source = '/home/sobolev/nevermind/Andrey/data'
+report = '/home/sobolev/nevermind/Andrey/analysis/'

analysis/target.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 
-def build_tgt_matrix(tl, aeps_events, tgt_dur):
+def build_tgt_matrix(tl, trials, aeps_events):
     # 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
@@ -32,9 +32,9 @@ def build_tgt_matrix(tl, aeps_events, tgt_dur):
 
     # 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
+    for idx_tl_success_end in trials[trials[:, 5] == 1][:, 1]:
+        idx_succ = np.abs(tl_tgt_end_idxs - idx_tl_success_end).argmin()
+        tgt_results[idx_succ] = 1
 
     # tl_idx_start, tl_idx_end, aep_idx_start, aer_idx_end, success / miss
     return np.column_stack([

postprocessing/execute.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 5,
    "id": "c37dce82",
    "metadata": {
     "scrolled": true
@@ -39,7 +39,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 6,
    "id": "a8ea2991",
    "metadata": {},
    "outputs": [
@@ -75,7 +75,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 7,
    "id": "3891b792",
    "metadata": {},
    "outputs": [],
@@ -100,16 +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-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",
+    "#'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",
@@ -149,7 +149,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 8,
    "id": "24cd2f6e",
    "metadata": {},
    "outputs": [
@@ -157,16 +157,7 @@
      "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-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"
+      "session 009266_hippoSIT_2023-04-26_08-20-17 done\n"
      ]
     }
    ],

reporting/State-dependent neural dynamics in the gerbil Auditory Cortex/functions.ipynb

@@ -0,0 +1,204 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "0642d602",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys, os\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..'))\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..', '..'))\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..', '..'))\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..', '..', 'analysis'))\n",
+    "sys.path.append(os.path.join(os.getcwd(), '..', '..', 'session'))\n",
+    "\n",
+    "import numpy as np\n",
+    "from imports import *\n",
+    "from matplotlib.patches import ConnectionPatch\n",
+    "from scipy.stats import pearsonr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "3c41fad7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_spike_counts(spk_times, pulse_times, hw=0.25, bin_count=51):\n",
+    "    collected = []\n",
+    "    for t_pulse in pulse_times:\n",
+    "        selected = spk_times[(spk_times > t_pulse - hw) & (spk_times < t_pulse + hw)]\n",
+    "        collected += [x for x in selected - t_pulse]\n",
+    "    collected = np.array(collected)\n",
+    "\n",
+    "    bins = np.linspace(-hw, hw, bin_count)\n",
+    "    counts, _ = np.histogram(collected, bins=bins)\n",
+    "    counts = (counts / len(pulse_times))# * 1/((2. * hw)/float(bin_count - 1))\n",
+    "    \n",
+    "    return bins, counts"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "1027884b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_tgt_bgr_psth(example_units):\n",
+    "    unit_count = np.array([len(vals) for vals in example_units.values()]).sum()\n",
+    "    rows = int(np.ceil(unit_count/3))\n",
+    "    fig = plt.figure(figsize=(15, rows*4))\n",
+    "    count = 0\n",
+    "\n",
+    "    for session, unit_ids in example_units.items():\n",
+    "        # read AEP events\n",
+    "        animal    = session.split('_')[0]\n",
+    "        aeps_file = os.path.join(source, animal, session, 'AEPs.h5')\n",
+    "        with h5py.File(aeps_file, 'r') as f:\n",
+    "            aeps_events = np.array(f['aeps_events'])\n",
+    "\n",
+    "        # read single units\n",
+    "        spike_times = {}\n",
+    "        h5_file = os.path.join(source, animal, session, session + '.h5')\n",
+    "        with h5py.File(h5_file, 'r') as f:\n",
+    "            cfg = json.loads(f['processed'].attrs['parameters'])\n",
+    "            for unit_id in unit_ids:\n",
+    "                spike_times[unit_id] = np.array(f['units'][unit_id][H5NAMES.spike_times['name']])\n",
+    "\n",
+    "        for unit_id in unit_ids:\n",
+    "            bins, counts_bgr = get_spike_counts(spike_times[unit_id], aeps_events[aeps_events[:, 1] == 1][:, 0])\n",
+    "            bins, counts_tgt = get_spike_counts(spike_times[unit_id], aeps_events[aeps_events[:, 1] == 2][:, 0])\n",
+    "\n",
+    "            ax = fig.add_subplot(rows, 3, count+1)\n",
+    "            tgt_dur, bgr_dur = cfg['sound']['sounds']['target']['duration'], cfg['sound']['sounds']['background']['duration']\n",
+    "            label_tgt = \"Tgt: %.2f\" % tgt_dur\n",
+    "            label_bgr = \"Bgr: %.2f\" % bgr_dur\n",
+    "            ax.hist(bins[:-1], bins=bins, weights=counts_tgt, edgecolor='black', color='tab:orange', alpha=0.9, label=label_tgt)\n",
+    "            ax.hist(bins[:-1], bins=bins, weights=counts_bgr, edgecolor='black', color='black', alpha=0.5, label=label_bgr)\n",
+    "            ax.axvline(0, color='black')\n",
+    "            ax.axvline(tgt_dur, color='tab:orange', ls='--')\n",
+    "            ax.axvline(tgt_dur - 0.25, color='tab:orange', ls='--')\n",
+    "            ax.axvline(bgr_dur, color='black', ls='--', alpha=0.5)\n",
+    "            ax.axvline(bgr_dur - 0.25, color='black', ls='--', alpha=0.5)\n",
+    "            #ax.set_xlabel('Pulse onset, s', fontsize=14)\n",
+    "            ax.axvspan(0, 0.05, alpha=0.3, color='gray')\n",
+    "            ax.set_title(\"%s : %s\" % (session[21:], unit_id), fontsize=14)\n",
+    "            ax.legend(loc='upper right', prop={'size': 10})\n",
+    "            if count % 3 == 0:\n",
+    "                ax.set_ylabel(\"Firing Rate, Hz\", fontsize=14)\n",
+    "            count += 1\n",
+    "        \n",
+    "    return fig"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "d97b4e42",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_psth_by_metric(area, m_name, example_units):\n",
+    "    unit_count = np.array([len(vals) for vals in example_units.values()]).sum()\n",
+    "    rows = int(np.ceil(unit_count/3))\n",
+    "    fig = plt.figure(figsize=(15, rows*4))\n",
+    "    count = 0\n",
+    "\n",
+    "    for session, unit_ids in example_units.items():\n",
+    "        # read AEP events\n",
+    "        animal    = session.split('_')[0]\n",
+    "        aeps_file = os.path.join(source, animal, session, 'AEPs.h5')\n",
+    "        with h5py.File(aeps_file, 'r') as f:\n",
+    "            aeps_events = np.array(f['aeps_events'])\n",
+    "            aeps = np.array(f[area]['aeps'])\n",
+    "\n",
+    "        # TODO find better way. Remove outliers\n",
+    "        aeps[aeps > 5000]  =  5000\n",
+    "        aeps[aeps < -5000] = -5000\n",
+    "\n",
+    "        # read single units\n",
+    "        spike_times = {}\n",
+    "        h5_file = os.path.join(source, animal, session, session + '.h5')\n",
+    "        with h5py.File(h5_file, 'r') as f:\n",
+    "            for unit_id in unit_ids:\n",
+    "                spike_times[unit_id] = np.array(f['units'][unit_id][H5NAMES.spike_times['name']])\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",
+    "            grp = f[area]\n",
+    "            for metric_name in grp['raw']:\n",
+    "                AEP_metrics_raw[metric_name]  = np.array(grp['raw'][metric_name])\n",
+    "                AEP_metrics_norm[metric_name] = np.array(grp['norm'][metric_name])\n",
+    "                AEP_metrics_lims[metric_name] = [int(x) for x in grp['raw'][metric_name].attrs['limits'].split(',')]\n",
+    "\n",
+    "        # separate high / low AEP metric states\n",
+    "        predictor = AEP_metrics_norm[m_name]\n",
+    "        low_state_idxs  = np.where(predictor < predictor.mean())[0]\n",
+    "        high_state_idxs = np.where(predictor > predictor.mean())[0]\n",
+    "        aeps_low_mean  = aeps[low_state_idxs].mean(axis=0)\n",
+    "        aeps_high_mean = aeps[high_state_idxs].mean(axis=0)\n",
+    "\n",
+    "        for unit_id in unit_ids:\n",
+    "            bins, counts_low  = get_spike_counts(spike_times[unit_id], aeps_events[low_state_idxs][:, 0])\n",
+    "            bins, counts_high = get_spike_counts(spike_times[unit_id], aeps_events[high_state_idxs][:, 0])\n",
+    "\n",
+    "            vals_max = np.array([counts_high.max(), counts_low.max()]).max()\n",
+    "            aep_low_profile  = (1/10) * vals_max * (aeps_low_mean/500)\n",
+    "            aep_high_profile = (1/10) * vals_max * (aeps_high_mean/500)\n",
+    "\n",
+    "            ax = fig.add_subplot(rows, 3, count+1)\n",
+    "            ax.hist(bins[:-1], bins=bins, weights=counts_high, edgecolor='black', color='red', alpha=0.8, label='%s >' % m_name)\n",
+    "            ax.hist(bins[:-1], bins=bins, weights=counts_low, edgecolor='black', color='black', alpha=0.5, label='%s <' % m_name)\n",
+    "            for x_l, x_r in [(-0.25, -0.051), (0.0, 0.199)]:\n",
+    "                ax.plot(np.linspace(x_l, x_r, len(aeps_low_mean)),  aep_high_profile, color='red', lw=2)\n",
+    "                ax.plot(np.linspace(x_l, x_r, len(aeps_high_mean)), aep_low_profile, color='black', lw=2)\n",
+    "            ax.axvline(0, color='black', ls='--')\n",
+    "            #ax.set_xlabel('Pulse onset, s', fontsize=14)\n",
+    "            ax.axvspan(0, 0.05, alpha=0.3, color='gray')\n",
+    "            ax.set_title(\"%s : %s\" % (session[21:], unit_id), fontsize=14)\n",
+    "            ax.legend(loc='upper right', prop={'size': 10})\n",
+    "            if count % 3 == 0:\n",
+    "                ax.set_ylabel(\"Firing Rate, Hz\", fontsize=14)\n",
+    "            count += 1\n",
+    "\n",
+    "    return fig"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "04c3c93b",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

session/adapters.py

@@ -135,16 +135,16 @@ class DatProcessor:
         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):
+    def get_single_channel(self, channel_no, block_duration=1):  # block duration 1 sec
         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
+        while offset < samples_no / self.s_rate - block_duration:
+            block = self.read_block_from_dat(block_duration, offset)  # read in 1 sec blocks
+            raw_signal[self.s_rate*offset:self.s_rate*(offset + block_duration)] = block[:, channel_no]
+            offset += block_duration
 
         return raw_signal

session/sessions.py

@@ -0,0 +1,25 @@
+selected_009266 = [
+# frequency
+'009266_hippoSIT_2023-04-17_17-04-17',  # ch17, 20 + 55 correction, 5067 events. Showcase for N2 / N3 mod in target
+'009266_hippoSIT_2023-04-18_10-10-37',  # ch17, 10 + 55 correction, 5682 events
+'009266_hippoSIT_2023-04-18_17-03-10',  # ch17, 6 + 55 correction, 5494 events: FIXME very weird 1-2nd in target, find out
+'009266_hippoSIT_2023-04-19_10-33-51',  # ch17, 4 + 55 correction, 6424 events: very weird 1-2nd in target, find out
+'009266_hippoSIT_2023-04-24_16-56-55',  # ch17, 5 + 55* correction, 6165 events, frequency
+'009266_hippoSIT_2023-05-02_12-22-14',  # ch20, 10 + 55 correction, 5976 events, FIXME very weird 1-2nd in target, find out
+    
+# duration
+'009266_hippoSIT_2023-04-20_08-57-39',  # ch1, 1 + 55 correction, 6424 events. Showcase for N2 / N3 mod in target
+'009266_hippoSIT_2023-04-26_08-20-17',  # ch17, 12 + 55* correction, 6095 events, duration - showcase for N2 
+'009266_hippoSIT_2023-05-04_19-47-15',  # ch20, 2 + 55 correction, 5678 events, duration
+
+# COMA
+'009266_hippoSIT_2023-05-04_09-11-06',  # ch17, 5 + 55* correction, 4487 events, COMA session with baseline AEPs
+]
+
+selected_008229 = [
+    # Old PPC sessions
+    '008229_hippoSIT_2022-05-17_21-44-43',  # chs: 0, 31, 54, 56
+    '008229_hippoSIT_2022-05-16_20-36-44',  # chs: 0, 56
+    '008229_hippoSIT_2022-05-20_15-54-39',  # chs: 0, 56
+    '008229_hippoSIT_2022-05-18_14-36-18',  # chs: 0, 56
+]