UMAP / tSNE neuro / behavioral states

analysis/Behavior/MoSeq Ratio MX GLM to W1 - W4.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 1,
    "id": "4001a375",
    "metadata": {},
    "outputs": [],
@@ -15,7 +15,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 2,
    "id": "570920e2",
    "metadata": {},
    "outputs": [],
@@ -40,7 +40,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 3,
    "id": "26ba4176",
    "metadata": {},
    "outputs": [
@@ -68,7 +68,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 4,
    "id": "4daf16cf",
    "metadata": {
     "scrolled": true
@@ -80,7 +80,7 @@
        "['009266_hippoSIT_2023-04-17_17-04-17', '009266_hippoSIT_2023-04-18_10-10-37']"
       ]
      },
-     "execution_count": 36,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -103,7 +103,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 5,
    "id": "c3501bc7",
    "metadata": {},
    "outputs": [],
@@ -135,7 +135,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 6,
    "id": "d0dd41e8",
    "metadata": {
     "scrolled": true
@@ -145,15 +145,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "W1 - original: 0.09, 0.000; shuffled: -0.06, 0.013\n",
-      "W2 - original: 0.29, 0.000; shuffled: -0.00, 0.848\n",
-      "W3 - original: 0.32, 0.000; shuffled: 0.01, 0.709\n",
-      "W4 - original: 0.43, 0.000; shuffled: 0.02, 0.530\n",
+      "W1 - original: 0.09, 0.000; shuffled: 0.03, 0.246\n",
+      "W2 - original: 0.32, 0.000; shuffled: -0.03, 0.223\n",
+      "W3 - original: 0.35, 0.000; shuffled: 0.00, 0.932\n",
+      "W4 - original: 0.44, 0.000; shuffled: 0.04, 0.146\n",
       "009266_hippoSIT_2023-04-17_17-04-17 done\n",
-      "W1 - original: 0.18, 0.000; shuffled: 0.01, 0.658\n",
-      "W2 - original: 0.28, 0.000; shuffled: 0.01, 0.661\n",
-      "W3 - original: 0.26, 0.000; shuffled: -0.02, 0.363\n",
-      "W4 - original: 0.54, 0.000; shuffled: -0.00, 0.989\n",
+      "W1 - original: 0.14, 0.000; shuffled: 0.00, 0.976\n",
+      "W2 - original: 0.32, 0.000; shuffled: 0.02, 0.423\n",
+      "W3 - original: 0.24, 0.000; shuffled: 0.01, 0.707\n",
+      "W4 - original: 0.53, 0.000; shuffled: 0.01, 0.669\n",
       "009266_hippoSIT_2023-04-18_10-10-37 done\n"
      ]
     }
@@ -165,7 +165,8 @@
     "    animal = session.split('_')[0]\n",
     "    s_path           = os.path.join(source, animal, session)\n",
     "    meta_file        = os.path.join(source, animal, session, 'meta.h5')\n",
-    "    moseq_class_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_class_500.h5')\n",
+    "    moseq_file       = os.path.join(source, animal, session, 'MoSeq.h5')\n",
+    "    moseq_class_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')\n",
     "    moseq_fit_file   = os.path.join(source, animal, session, 'analysis', 'MoSeq_W1-W4_GLM.h5')\n",
     "    \n",
     "    with h5py.File(meta_file, 'r') as f:\n",
@@ -175,24 +176,36 @@
     "    with h5py.File(moseq_class_file, 'r') as f:\n",
     "        syl_ratio_mx = np.array(f['syl_ratio_mx'])\n",
     "        idxs_srm_tl  = np.array(f['idxs_srm_tl'])\n",
+    "    with h5py.File(moseq_file, 'r') as f:\n",
+    "        moseq = np.array(f['moseq'])\n",
     "        \n",
+    "    width = 70  # 100 points ~= 1 sec with at 100Hz\n",
+    "    kernel = signal.gaussian(width, std=(width) / 7.2)\n",
+    "    dx = np.sqrt(np.square(np.diff(moseq[:, 3])) + np.square(np.diff(moseq[:, 4])))\n",
+    "    dt = np.diff(moseq[:, 0])\n",
+    "    speed = np.concatenate([dx/dt, [dx[-1]/dt[-1]]])\n",
+    "    speed_smooth = np.convolve(speed, kernel, 'same') / kernel.sum()\n",
+    "    \n",
     "    for j, phase in enumerate([1, 2, 3, 4]):\n",
     "        w_pca = activity_at_phase(s_path, phase, do_pca=True)\n",
     "        w_int = np.interp(tl[idxs_srm_tl][:, 0], events[:, 0], w_pca)\n",
+    "\n",
+    "        pred = syl_ratio_mx\n",
+    "        pred = np.column_stack([syl_ratio_mx, speed_smooth[idxs_srm_tl]])\n",
     "        \n",
     "        # original\n",
-    "        corr, pval, params, pvalues = get_GLM_and_prediction(syl_ratio_mx, w_int)\n",
+    "        corr, pval, params, pvalues = get_GLM_and_prediction(pred, w_int)\n",
     "        \n",
     "        # diff train / split\n",
     "        corr_mx_chun = np.zeros([iter_count, 2])\n",
     "        for k in range(iter_count):\n",
-    "            corr_s, pval_s, _, _ = get_GLM_and_prediction(syl_ratio_mx, w_int)\n",
+    "            corr_s, pval_s, _, _ = get_GLM_and_prediction(pred, w_int)\n",
     "            corr_mx_chun[k] = (corr_s, pval_s)\n",
     "            \n",
     "        # shuffled\n",
     "        corr_mx_shuf = np.zeros([iter_count, 2])  # coeff, pval for each shuffle\n",
     "        for k in range(iter_count):\n",
-    "            syl_ratio_mx_s = syl_ratio_mx.copy()\n",
+    "            syl_ratio_mx_s = pred.copy()\n",
     "            np.random.shuffle(syl_ratio_mx_s)\n",
     "            corr_s, pval_s, _, _ = get_GLM_and_prediction(syl_ratio_mx_s, w_int)\n",
     "            corr_mx_shuf[k] = (corr_s, pval_s)\n",
@@ -217,28 +230,6 @@
     "    print(session + ' done')"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 51,
-   "id": "b9011b94",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "array([385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397,\n",
-       "       398, 399, 400, 401, 402, 403])"
-      ]
-     },
-     "execution_count": 51,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "np.arange(385, 404)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 45,
@@ -321,6 +312,22 @@
    "metadata": {},
    "outputs": [],
    "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b26aafe",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0bd036ed",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {

analysis/Behavior/MoSeq UMAP BRG TGT SIL success.ipynb

@@ -38,7 +38,11 @@
    "outputs": [
     {
      "data": {
-      "application/javascript": "IPython.OutputArea.prototype._should_scroll = function(lines) {\n    return false;\n}\n",
+      "application/javascript": [
+       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
+       "    return false;\n",
+       "}\n"
+      ],
       "text/plain": [
        "<IPython.core.display.Javascript object>"
       ]
@@ -175,9 +179,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pysit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "pysit"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -189,7 +193,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.8.10"
   }
  },
  "nbformat": 4,

analysis/Behavior/MoSeq UMAP pellet chasing.ipynb

@@ -39,7 +39,11 @@
    "outputs": [
     {
      "data": {
-      "application/javascript": "IPython.OutputArea.prototype._should_scroll = function(lines) {\n    return false;\n}\n",
+      "application/javascript": [
+       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
+       "    return false;\n",
+       "}\n"
+      ],
       "text/plain": [
        "<IPython.core.display.Javascript object>"
       ]
@@ -225,9 +229,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pysit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "pysit"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -239,7 +243,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.8.10"
   }
  },
  "nbformat": 4,

analysis/Behavior/MoSeq W4 Separate by success - fail.ipynb

@@ -40,7 +40,11 @@
    "outputs": [
     {
      "data": {
-      "application/javascript": "IPython.OutputArea.prototype._should_scroll = function(lines) {\n    return false;\n}\n",
+      "application/javascript": [
+       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
+       "    return false;\n",
+       "}\n"
+      ],
       "text/plain": [
        "<IPython.core.display.Javascript object>"
       ]
@@ -238,9 +242,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pysit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "pysit"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -252,7 +256,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.8.10"
   }
  },
  "nbformat": 4,

analysis/Behavior/behavior.py

@@ -71,6 +71,7 @@ def get_idxs_in_patches(fit, patches, extent, bin_count=100):
 
 
 def get_idxs_behav_state(source, session, idxs_tl_sample, fit_type='tSNE', fit_parm=70, sigma=0.3, margin=10, bin_count=100):
+    # returns idxs to timeline!
     animal = session.split('_')[0]
     meta_file        = os.path.join(source, animal, session, 'meta.h5')
     moseq_class_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')
@@ -85,7 +86,36 @@ def get_idxs_behav_state(source, session, idxs_tl_sample, fit_type='tSNE', fit_p
     idxs_state = np.array([i for i, x in enumerate(idxs_srm_tl) if x in idxs_tl_sample], dtype=np.int32)
     
     extent = get_extent(fit, margin=margin)
-    behav_map     = density_map(fit[idxs_state], extent, sigma=sigma, bin_count=bin_count)
-    state_patches = get_field_patches(behav_map, 0.2)
-    return get_idxs_in_patches(fit, state_patches, extent, bin_count=bin_count)
+    behav_map      = density_map(fit[idxs_state], extent, sigma=sigma, bin_count=bin_count)
+    state_patches  = get_field_patches(behav_map, 0.2)
+    idxs_srm_state = get_idxs_in_patches(fit, state_patches, extent, bin_count=bin_count)
+    
+    # convert to timeline idxs
+    bins_to_fill = int((idxs_srm_tl[1] - idxs_srm_tl[0])/2)
+    idxs_res = []
+    for idx in idxs_srm_tl[idxs_srm_state]:
+        idxs_res += list(range(idx - bins_to_fill, idx + bins_to_fill))
+    idxs_res = np.array(idxs_res, dtype=np.int32)
+    idxs_res = idxs_res[idxs_res > 0]
+    idxs_res = idxs_res[idxs_res < len(tl) - 1]
+        
+    return idxs_res
+    
+    
+def get_idxs_neuro_state(source, session, idxs_ev_sample, fit_type='tSNE', fit_parm=70, sigma=0.3, margin=10, bin_count=100):
+    # returns idxs in sound events space
+    animal = session.split('_')[0]
+    meta_file        = os.path.join(source, animal, session, 'meta.h5')
+    umap_file = os.path.join(source, animal, session, 'analysis', 'W1-W4_tSNE_UMAP.h5')
+
+    with h5py.File(meta_file, 'r') as f:
+        tl = np.array(f['processed']['timeline'])
+        tgt_mx = np.array(f['processed']['target_matrix'])
+    with h5py.File(umap_file, 'r') as f:
+        fit = np.array(f[fit_type][str(fit_parm)])  # already in event sampling
+
+    extent = get_extent(fit, margin=margin)
+    selected_map  = density_map(fit[idxs_ev_sample], extent, sigma=sigma, bin_count=bin_count)
+    state_patches = get_field_patches(selected_map, 0.3)
 
+    return get_idxs_in_patches(fit, state_patches, extent, bin_count=bin_count)

analysis/Behavior/dPCA.ipynb

@@ -0,0 +1,78 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "bca0a270",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from numpy import *\n",
+    "from numpy.random import rand, randn, randint\n",
+    "from dPCA import dPCA"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "1fef3742",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# number of neurons, time-points and stimuli\n",
+    "N,T,S = 100, 250, 6\n",
+    "\n",
+    "# noise-level and number of trials in each condition\n",
+    "noise, n_samples = 0.2, 10\n",
+    "\n",
+    "# build two latent factors\n",
+    "zt = (arange(T)/float(T))\n",
+    "zs = (arange(S)/float(S))\n",
+    "\n",
+    "# build trial-by trial data\n",
+    "trialR = noise*randn(n_samples,N,S,T)\n",
+    "trialR += randn(N)[None,:,None,None]*zt[None,None,None,:]\n",
+    "trialR += randn(N)[None,:,None,None]*zs[None,None,:,None]\n",
+    "\n",
+    "# trial-average data\n",
+    "R = mean(trialR,0)\n",
+    "\n",
+    "# center data\n",
+    "R -= mean(R.reshape((N,-1)),1)[:,None,None]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eaa5ef01",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a dataset for dPCA\n",
+    "# n, r, s, t = neuron, run/rest, state, time\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pysit",
+   "language": "python",
+   "name": "pysit"
+  },
+  "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

analysis/Population/GLM Speed + State => Pop. Act..ipynb

@@ -0,0 +1,310 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "f5ab8420",
+   "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(), '..', '..', '..', 'pplSIT', 'workflow', 'utils'))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "5a7026b5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import pandas as pd\n",
+    "from session.sessions import selected_009266, selected_008229, selected_009265\n",
+    "from imports import *\n",
+    "from scipy import stats\n",
+    "from scipy import signal\n",
+    "from sklearn.manifold import TSNE\n",
+    "from sklearn import decomposition\n",
+    "from sklearn.model_selection import train_test_split\n",
+    "from statsmodels.formula.api import ols, glm\n",
+    "\n",
+    "from Behavior.behavior import get_extent, get_idxs_behav_state\n",
+    "from population import unit_response_matrix, activity_at_phase"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "3d645f04",
+   "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": 12,
+   "id": "72133ead",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_GLM_and_prediction(syl_ratio_mx, pop_at_phase, test_size=0.33, glm_min_pval=0.95):\n",
+    "\n",
+    "    # separate train / test\n",
+    "    X_train, X_test, y_train, y_test = train_test_split(syl_ratio_mx, pop_at_phase, test_size=test_size)\n",
+    "    \n",
+    "    # train glm to get contributions of each syllable\n",
+    "    data = np.column_stack([y_train, X_train])\n",
+    "    columns = ['state'] + [\"x%d\" % x for x in range(X_train.shape[1])]\n",
+    "    AM_df = pd.DataFrame(data, columns=columns)\n",
+    "\n",
+    "    model = glm('state ~ ' + ' + '.join(columns[1:]), data=AM_df).fit()\n",
+    "    #model.summary()\n",
+    "    \n",
+    "    glm_coeffs = dict([(i, coef) for i, coef in enumerate(model.params[1:]) if model.pvalues[1:][i] < glm_min_pval])\n",
+    "    #if len(glm_coeffs) == 0:\n",
+    "    #    return 0, 0, model.params, model.pvalues\n",
+    "    target_fit = np.zeros(len(y_test))\n",
+    "    for idx, coef in glm_coeffs.items():\n",
+    "        target_fit += coef * X_test[:, idx]\n",
+    "\n",
+    "    corr, pval = stats.pearsonr(target_fit, y_test)\n",
+    "    \n",
+    "    return corr, pval, model.params, model.pvalues"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "id": "851f1f03",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "['009265_hippoSIT_2023-02-24_09-53-26',\n",
+       " '009265_hippoSIT_2023-02-24_17-22-46',\n",
+       " '009265_hippoSIT_2023-02-28_09-16-50',\n",
+       " '009265_hippoSIT_2023-02-28_13-16-10',\n",
+       " '009265_hippoSIT_2023-02-28_20-45-04',\n",
+       " '009265_hippoSIT_2023-03-01_10-46-12',\n",
+       " '009265_hippoSIT_2023-03-02_09-32-54',\n",
+       " '009265_hippoSIT_2023-03-02_16-27-42',\n",
+       " '009265_hippoSIT_2023-03-02_20-11-35',\n",
+       " '009265_hippoSIT_2023-03-03_09-37-07',\n",
+       " '009265_hippoSIT_2023-03-03_16-00-47',\n",
+       " '009265_hippoSIT_2023-03-04_11-12-04',\n",
+       " '009265_hippoSIT_2023-03-05_11-52-17',\n",
+       " '009265_hippoSIT_2023-03-05_18-31-32',\n",
+       " '009265_hippoSIT_2023-03-08_18-10-07',\n",
+       " '009265_hippoSIT_2023-03-09_20-03-08',\n",
+       " '009265_hippoSIT_2023-03-10_09-57-34',\n",
+       " '009265_hippoSIT_2023-04-13_09-54-39',\n",
+       " '009265_hippoSIT_2023-04-20_11-39-02']"
+      ]
+     },
+     "execution_count": 41,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "source = '/home/sobolev/nevermind_ag-grothe/AG_Pecka/data/processed/'\n",
+    "source = '/home/sobolev/nevermind/AG_Pecka/data/processed/'\n",
+    "\n",
+    "sessions = [s for s in selected_009266.keys()]\n",
+    "sessions = [s for s in selected_009265.keys()]\n",
+    "#sessions = [s for s in selected_008229.keys()]\n",
+    "sessions.sort()\n",
+    "\n",
+    "selected = sessions[:]\n",
+    "try:\n",
+    "    selected.remove('009266_hippoSIT_2023-04-20_15-24-14')\n",
+    "except:\n",
+    "    pass\n",
+    "try:\n",
+    "    selected.remove('009265_hippoSIT_2023-02-27_10-18-32')\n",
+    "    selected.remove('009265_hippoSIT_2023-02-27_15-33-46')\n",
+    "except:\n",
+    "    pass\n",
+    "selected"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "id": "643e902f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2023-02-24_09 - Pred: 0.10, pval: 0.000; Coeffs: -0.01, -0.11\n",
+      "2023-02-24_17 - Pred: 0.19, pval: 0.000; Coeffs: 0.16, -0.28\n",
+      "2023-02-28_09 - Pred: 0.04, pval: 0.024; Coeffs: 0.01, -0.02\n",
+      "2023-02-28_13 - Pred: 0.07, pval: 0.000; Coeffs: 0.42, 0.04\n",
+      "2023-02-28_20 - Pred: 0.18, pval: 0.000; Coeffs: 0.55, -0.07\n",
+      "2023-03-01_10 - Pred: 0.10, pval: 0.000; Coeffs: 0.16, -0.13\n",
+      "2023-03-02_09 - Pred: 0.02, pval: 0.232; Coeffs: 0.06, 0.00\n",
+      "2023-03-02_16 - Pred: 0.21, pval: 0.000; Coeffs: 0.77, -0.15\n",
+      "2023-03-02_20 - Pred: 0.34, pval: 0.000; Coeffs: 0.89, -0.37\n",
+      "2023-03-03_09 - Pred: 0.44, pval: 0.000; Coeffs: 0.61, -0.53\n",
+      "2023-03-03_16 - Pred: 0.27, pval: 0.000; Coeffs: 0.68, -0.28\n",
+      "2023-03-04_11 - Pred: 0.27, pval: 0.000; Coeffs: 0.39, -0.29\n",
+      "2023-03-05_11 - Pred: 0.26, pval: 0.000; Coeffs: 0.39, -0.20\n",
+      "2023-03-05_18 - Pred: 0.07, pval: 0.007; Coeffs: 0.14, -0.17\n",
+      "2023-03-08_18 - Pred: 0.01, pval: 0.688; Coeffs: 0.19, 0.07\n",
+      "2023-03-09_20 - Pred: 0.20, pval: 0.000; Coeffs: 0.48, -0.45\n",
+      "2023-03-10_09 - Pred: 0.20, pval: 0.000; Coeffs: 0.44, -0.18\n",
+      "2023-04-13_09 - Pred: 0.40, pval: 0.000; Coeffs: 0.62, -0.47\n",
+      "2023-04-20_11 - Pred: 0.31, pval: 0.000; Coeffs: 0.57, -0.40\n"
+     ]
+    }
+   ],
+   "source": [
+    "#session = selected[0]\n",
+    "\n",
+    "for session in selected:\n",
+    "    ft = 'tSNE'\n",
+    "    fp = 70\n",
+    "    animal = session.split('_')[0]\n",
+    "    session_path     = os.path.join(source, animal, session)\n",
+    "    meta_file        = os.path.join(source, animal, session, 'meta.h5')\n",
+    "    moseq_file       = os.path.join(source, animal, session, 'MoSeq.h5')\n",
+    "    moseq_class_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')\n",
+    "\n",
+    "    with h5py.File(meta_file, 'r') as f:\n",
+    "        tl = np.array(f['processed']['timeline'])\n",
+    "        tgt_mx = np.array(f['processed']['target_matrix'])\n",
+    "        events = np.array(f['processed']['sound_events'])\n",
+    "    with h5py.File(moseq_class_file, 'r') as f:\n",
+    "        idxs_srm_tl = np.array(f['idxs_srm_tl'])\n",
+    "        fit = np.array(f[ft][str(fp)])\n",
+    "    with h5py.File(moseq_file, 'r') as f:\n",
+    "        moseq = np.array(f['moseq'])\n",
+    "\n",
+    "    # take only silence\n",
+    "    idxs_sil_ev = np.where(events[:, 1] == 0)[0]\n",
+    "    idxs_sil_ev = np.where( (events[:, 1] == 0) | (events[:, 1] == 1) )[0]\n",
+    "    #idxs_sil_ev = np.arange(len(events))\n",
+    "\n",
+    "    # 1. Speed\n",
+    "    width = 70  # 100 points ~= 1 sec with at 100Hz\n",
+    "    kernel = signal.gaussian(width, std=(width) / 7.2)\n",
+    "    dx = np.sqrt(np.square(np.diff(moseq[:, 3])) + np.square(np.diff(moseq[:, 4])))\n",
+    "    dt = np.diff(moseq[:, 0])\n",
+    "    speed = np.concatenate([dx/dt, [dx[-1]/dt[-1]]])\n",
+    "    speed_smooth = np.convolve(speed, kernel, 'same') / kernel.sum()\n",
+    "    speed_ev = speed_smooth[events[:, 2].astype(np.int32)]\n",
+    "    speed_sil_ev = speed_ev[idxs_sil_ev]\n",
+    "\n",
+    "    # 2. Population activity \n",
+    "    bins, unit_mx = unit_response_matrix(session_path, [1, 2], times_to_event=[])\n",
+    "\n",
+    "    # z-score\n",
+    "    unit_act_matrix = unit_mx.T\n",
+    "    for u, unit_data in enumerate(unit_act_matrix):\n",
+    "        unit_act_matrix[u] = stats.zscore(unit_data)\n",
+    "    unit_mx = unit_act_matrix.T\n",
+    "\n",
+    "    # take first PC from PCA on multiple units\n",
+    "    pca = decomposition.PCA(n_components=3)\n",
+    "    pca.fit(unit_mx)\n",
+    "    X = pca.transform(unit_mx)\n",
+    "    pop_act = X[:, 0]  # PC1 score\n",
+    "\n",
+    "    # filter SILENCE\n",
+    "    pop_act = pop_act[idxs_sil_ev]\n",
+    "\n",
+    "    # smooth\n",
+    "    k_width = 30\n",
+    "    kernel  = signal.gaussian(k_width, std=(k_width) / 7.2)\n",
+    "    pop_act = np.convolve(pop_act, kernel, 'same') / kernel.sum()\n",
+    "\n",
+    "    # filter slow oscillations\n",
+    "    sos = signal.butter(10, 0.001, fs=4, analog=False, btype='highpass', output='sos')\n",
+    "    pop_act = signal.sosfiltfilt(sos, pop_act)\n",
+    "\n",
+    "    # just to test - W4 only\n",
+    "    w4 = activity_at_phase(session_path, phase=4, do_pca=True)\n",
+    "    w4 = w4[idxs_sil_ev]\n",
+    "\n",
+    "\n",
+    "    # 3. Behavioral state\n",
+    "    idxs_tl_tgt_succ = []\n",
+    "    for tgt_rec in tgt_mx[tgt_mx[:, 4] == 1]:\n",
+    "        idxs_tl_tgt_succ += list(range(tgt_rec[2], tgt_rec[3]))\n",
+    "    idxs_tl_tgt_succ = np.array(idxs_tl_tgt_succ)\n",
+    "    idxs_AL_tl = get_idxs_behav_state(source, session, idxs_tl_tgt_succ, fit_type='tSNE', fit_parm=70, sigma=0.3, margin=10, bin_count=100)\n",
+    "    idxs_AL_ev = np.array([k for k, x in enumerate(events) if x[2] in idxs_AL_tl], dtype=np.int32)\n",
+    "\n",
+    "    bs = np.zeros(len(events))\n",
+    "    bs[idxs_AL_ev] = 1\n",
+    "    bs = bs[idxs_sil_ev]\n",
+    "\n",
+    "    # try smth random\n",
+    "    rnd = np.random.rand(len(bs))\n",
+    "\n",
+    "    # 4. Compute GLM and prediction\n",
+    "    corr, pval, params, pvalues = get_GLM_and_prediction(np.column_stack([speed_sil_ev, bs, rnd]), pop_act)\n",
+    "\n",
+    "    print(\"%s - Pred: %.2f, pval: %.3f; Coeffs: %.2f, %.2f\" % (session[-19:-6], corr, pval, params[1], params[2]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e5ba407a",
+   "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
+}

analysis/Population/LR Speed + Pop. Act. => Behav State.ipynb

@@ -0,0 +1,437 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "28ab8ba6",
+   "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(), '..', '..', '..', 'pplSIT', 'workflow', 'utils'))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "5869aeac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import pandas as pd\n",
+    "from session.sessions import selected_009266, selected_008229, selected_009265\n",
+    "from imports import *\n",
+    "from scipy import stats\n",
+    "from scipy import signal\n",
+    "from sklearn.manifold import TSNE\n",
+    "from sklearn import decomposition\n",
+    "from sklearn.model_selection import train_test_split\n",
+    "from sklearn.linear_model import LogisticRegression\n",
+    "from sklearn.metrics import classification_report, confusion_matrix\n",
+    "from statsmodels.formula.api import ols, glm\n",
+    "from regression import get_GLM_and_prediction\n",
+    "\n",
+    "from Behavior.behavior import get_extent, get_idxs_behav_state\n",
+    "from population import unit_response_matrix, activity_at_phase"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "90eed9ee",
+   "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": 15,
+   "id": "0e120760",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "['009266_hippoSIT_2023-04-17_17-04-17',\n",
+       " '009266_hippoSIT_2023-04-18_10-10-37',\n",
+       " '009266_hippoSIT_2023-04-18_17-03-10',\n",
+       " '009266_hippoSIT_2023-04-19_10-33-51',\n",
+       " '009266_hippoSIT_2023-04-20_08-57-39',\n",
+       " '009266_hippoSIT_2023-04-21_08-43-00',\n",
+       " '009266_hippoSIT_2023-04-21_13-12-31',\n",
+       " '009266_hippoSIT_2023-04-24_10-08-11',\n",
+       " '009266_hippoSIT_2023-04-24_16-56-55',\n",
+       " '009266_hippoSIT_2023-04-26_08-20-17',\n",
+       " '009266_hippoSIT_2023-05-02_12-22-14',\n",
+       " '009266_hippoSIT_2023-05-04_19-47-15',\n",
+       " '009266_hippoSIT_2023-05-22_09-27-22',\n",
+       " '009266_hippoSIT_2023-05-23_09-18-05',\n",
+       " '009266_hippoSIT_2023-05-25_15-55-57',\n",
+       " '009266_hippoSIT_2023-06-14_08-21-23',\n",
+       " '009266_hippoSIT_2023-06-19_08-58-35']"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "source = '/home/sobolev/nevermind_ag-grothe/AG_Pecka/data/processed/'\n",
+    "source = '/home/sobolev/nevermind/AG_Pecka/data/processed/'\n",
+    "\n",
+    "sessions = [s for s in selected_009266.keys()]\n",
+    "#sessions = [s for s in selected_009265.keys()]\n",
+    "#sessions = [s for s in selected_008229.keys()]\n",
+    "sessions.sort()\n",
+    "\n",
+    "selected = sessions[:]\n",
+    "try:\n",
+    "    selected.remove('009266_hippoSIT_2023-04-20_15-24-14')\n",
+    "except:\n",
+    "    pass\n",
+    "try:\n",
+    "    selected.remove('009265_hippoSIT_2023-02-27_10-18-32')\n",
+    "    selected.remove('009265_hippoSIT_2023-02-27_15-33-46')\n",
+    "except:\n",
+    "    pass\n",
+    "selected"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "id": "442f3fdb",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2023-04-17_17 - Scores: (AL: 0.08)\n",
+      "0.74: Speed + W4\n",
+      "0.73: Speed + W4 + RND\n",
+      "0.68: Speed only\n",
+      "0.69: W4 only\n",
+      "0.47: RND\n",
+      "2023-04-18_10 - Scores: (AL: 0.22)\n",
+      "0.77: Speed + W4\n",
+      "0.77: Speed + W4 + RND\n",
+      "0.73: Speed only\n",
+      "0.71: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-18_17 - Scores: (AL: 0.24)\n",
+      "0.69: Speed + W4\n",
+      "0.69: Speed + W4 + RND\n",
+      "0.68: Speed only\n",
+      "0.59: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-19_10 - Scores: (AL: 0.48)\n",
+      "0.65: Speed + W4\n",
+      "0.63: Speed + W4 + RND\n",
+      "0.62: Speed only\n",
+      "0.60: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-20_08 - Scores: (AL: 0.34)\n",
+      "0.72: Speed + W4\n",
+      "0.71: Speed + W4 + RND\n",
+      "0.68: Speed only\n",
+      "0.69: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-21_08 - Scores: (AL: 0.41)\n",
+      "0.65: Speed + W4\n",
+      "0.63: Speed + W4 + RND\n",
+      "0.63: Speed only\n",
+      "0.58: W4 only\n",
+      "0.51: RND\n",
+      "2023-04-21_13 - Scores: (AL: 0.47)\n",
+      "0.57: Speed + W4\n",
+      "0.56: Speed + W4 + RND\n",
+      "0.57: Speed only\n",
+      "0.56: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-24_10 - Scores: (AL: 0.38)\n",
+      "0.68: Speed + W4\n",
+      "0.70: Speed + W4 + RND\n",
+      "0.66: Speed only\n",
+      "0.65: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-24_16 - Scores: (AL: 0.34)\n",
+      "0.64: Speed + W4\n",
+      "0.63: Speed + W4 + RND\n",
+      "0.62: Speed only\n",
+      "0.55: W4 only\n",
+      "0.49: RND\n",
+      "2023-04-26_08 - Scores: (AL: 0.30)\n",
+      "0.76: Speed + W4\n",
+      "0.77: Speed + W4 + RND\n",
+      "0.72: Speed only\n",
+      "0.75: W4 only\n",
+      "0.49: RND\n",
+      "2023-05-02_12 - Scores: (AL: 0.31)\n",
+      "0.70: Speed + W4\n",
+      "0.71: Speed + W4 + RND\n",
+      "0.67: Speed only\n",
+      "0.62: W4 only\n",
+      "0.51: RND\n",
+      "2023-05-04_19 - Scores: (AL: 0.19)\n",
+      "0.72: Speed + W4\n",
+      "0.77: Speed + W4 + RND\n",
+      "0.74: Speed only\n",
+      "0.70: W4 only\n",
+      "0.50: RND\n",
+      "2023-05-22_09 - Scores: (AL: 0.41)\n",
+      "0.67: Speed + W4\n",
+      "0.67: Speed + W4 + RND\n",
+      "0.67: Speed only\n",
+      "0.60: W4 only\n",
+      "0.50: RND\n",
+      "2023-05-23_09 - Scores: (AL: 0.30)\n",
+      "0.72: Speed + W4\n",
+      "0.69: Speed + W4 + RND\n",
+      "0.68: Speed only\n",
+      "0.69: W4 only\n",
+      "0.49: RND\n",
+      "2023-05-25_15 - Scores: (AL: 0.39)\n",
+      "0.68: Speed + W4\n",
+      "0.68: Speed + W4 + RND\n",
+      "0.65: Speed only\n",
+      "0.63: W4 only\n",
+      "0.47: RND\n",
+      "2023-06-14_08 - Scores: (AL: 0.32)\n",
+      "0.70: Speed + W4\n",
+      "0.69: Speed + W4 + RND\n",
+      "0.70: Speed only\n",
+      "0.56: W4 only\n",
+      "0.49: RND\n",
+      "2023-06-19_08 - Scores: (AL: 0.29)\n",
+      "0.68: Speed + W4\n",
+      "0.66: Speed + W4 + RND\n",
+      "0.65: Speed only\n",
+      "0.62: W4 only\n",
+      "0.49: RND\n"
+     ]
+    }
+   ],
+   "source": [
+    "#session = selected[0]\n",
+    "\n",
+    "conf_mxs = {}\n",
+    "scores   = {}\n",
+    "for session in selected[:]:\n",
+    "    ft = 'tSNE'\n",
+    "    fp = 70\n",
+    "    animal = session.split('_')[0]\n",
+    "    session_path     = os.path.join(source, animal, session)\n",
+    "    meta_file        = os.path.join(source, animal, session, 'meta.h5')\n",
+    "    moseq_file       = os.path.join(source, animal, session, 'MoSeq.h5')\n",
+    "    moseq_class_file = os.path.join(source, animal, session, 'analysis', 'MoSeq_tSNE_UMAP.h5')\n",
+    "\n",
+    "    with h5py.File(meta_file, 'r') as f:\n",
+    "        tl = np.array(f['processed']['timeline'])\n",
+    "        tgt_mx = np.array(f['processed']['target_matrix'])\n",
+    "        events = np.array(f['processed']['sound_events'])\n",
+    "    with h5py.File(moseq_class_file, 'r') as f:\n",
+    "        idxs_srm_tl = np.array(f['idxs_srm_tl'])\n",
+    "        fit = np.array(f[ft][str(fp)])\n",
+    "    with h5py.File(moseq_file, 'r') as f:\n",
+    "        moseq = np.array(f['moseq'])\n",
+    "\n",
+    "    # events filter\n",
+    "    idxs_filt_ev = np.where(events[:, 1] == 0)[0]  # silence filter\n",
+    "    idxs_filt_ev = np.where( (events[:, 1] == 0) | (events[:, 1] == 1) )[0]  # non-target filter\n",
+    "    #idxs_filt_ev = np.arange(len(events))  # no filter\n",
+    "\n",
+    "    # 1. Speed\n",
+    "    width = 70  # 100 points ~= 1 sec with at 100Hz\n",
+    "    kernel = signal.gaussian(width, std=(width) / 7.2)\n",
+    "    dx = np.sqrt(np.square(np.diff(moseq[:, 3])) + np.square(np.diff(moseq[:, 4])))\n",
+    "    dt = np.diff(moseq[:, 0])\n",
+    "    speed = np.concatenate([dx/dt, [dx[-1]/dt[-1]]])\n",
+    "    speed_smooth = np.convolve(speed, kernel, 'same') / kernel.sum()\n",
+    "    speed_ev = speed_smooth[events[:, 2].astype(np.int32)]\n",
+    "    speed_filt_ev = speed_ev[idxs_filt_ev]\n",
+    "\n",
+    "    # 2. Population activity \n",
+    "    bins, unit_mx = unit_response_matrix(session_path, [1, 2], times_to_event=[])\n",
+    "\n",
+    "    # z-score\n",
+    "    unit_act_matrix = unit_mx.T\n",
+    "    for u, unit_data in enumerate(unit_act_matrix):\n",
+    "        unit_act_matrix[u] = stats.zscore(unit_data)\n",
+    "    unit_mx = unit_act_matrix.T\n",
+    "\n",
+    "    # take first PC from PCA on multiple units\n",
+    "    pca = decomposition.PCA(n_components=3)\n",
+    "    pca.fit(unit_mx)\n",
+    "    X = pca.transform(unit_mx)\n",
+    "    pop_act = X[:, 0]  # PC1 score\n",
+    "\n",
+    "    # filter SILENCE\n",
+    "    pop_act_filt = pop_act[idxs_filt_ev]\n",
+    "\n",
+    "    # smooth\n",
+    "    k_width = 30\n",
+    "    kernel  = signal.gaussian(k_width, std=(k_width) / 7.2)\n",
+    "    pop_act_filt = np.convolve(pop_act_filt, kernel, 'same') / kernel.sum()\n",
+    "\n",
+    "    # filter slow oscillations\n",
+    "    sos = signal.butter(10, 0.001, fs=4, analog=False, btype='highpass', output='sos')\n",
+    "    pop_act_filt = signal.sosfiltfilt(sos, pop_act_filt)\n",
+    "\n",
+    "    # just to test - W4 only\n",
+    "    w4 = activity_at_phase(session_path, phase=4, do_pca=True)\n",
+    "    w4_filt_ev = w4[idxs_filt_ev]\n",
+    "\n",
+    "\n",
+    "    # 3. Behavioral state\n",
+    "    idxs_tl_tgt_succ = []\n",
+    "    for tgt_rec in tgt_mx[tgt_mx[:, 4] == 1]:\n",
+    "        idxs_tl_tgt_succ += list(range(tgt_rec[2], tgt_rec[3]))\n",
+    "    idxs_tl_tgt_succ = np.array(idxs_tl_tgt_succ)\n",
+    "    idxs_AL_tl = get_idxs_behav_state(source, session, idxs_tl_tgt_succ, fit_type='tSNE', fit_parm=70, sigma=0.3, margin=10, bin_count=100)\n",
+    "    idxs_AL_ev = np.array([k for k, x in enumerate(events) if x[2] in idxs_AL_tl], dtype=np.int32)\n",
+    "\n",
+    "    bs = np.zeros(len(events))\n",
+    "    bs[idxs_AL_ev] = 1\n",
+    "    bs_filt_ev = bs[idxs_filt_ev]\n",
+    "\n",
+    "    # try smth random to show that random staff doesn't affect the score\n",
+    "    rnd = np.random.rand(len(bs_filt_ev))\n",
+    "\n",
+    "    \n",
+    "    # 4. Compute LR and prediction\n",
+    "    combs = [\n",
+    "        np.column_stack([speed_filt_ev, w4_filt_ev]),\n",
+    "        np.column_stack([speed_filt_ev, w4_filt_ev, rnd]),\n",
+    "        np.column_stack([speed_filt_ev, rnd]),\n",
+    "        np.column_stack([w4_filt_ev, rnd]),\n",
+    "        np.column_stack([rnd])\n",
+    "    ]\n",
+    "    labels = ['Speed + W4', 'Speed + W4 + RND', 'Speed only', 'W4 only', 'RND']\n",
+    "\n",
+    "    # 5a. Logistic regression\n",
+    "    # try to compensate for bias\n",
+    "    bs_high = np.where(bs_filt_ev > 0)[0]\n",
+    "    bs_low  = np.where(bs_filt_ev == 0)[0]\n",
+    "    if len(bs_high) != len(bs_low):\n",
+    "        idx_rem = 0 if len(bs_high) < len(bs_low) else 1\n",
+    "        idxs_exclude = []\n",
+    "        while len(idxs_exclude) < int(np.abs(len(bs_high) - len(bs_low))):\n",
+    "            idx = np.random.randint(len(bs_filt_ev))\n",
+    "            if bs_filt_ev[idx] == idx_rem and not idx in idxs_exclude:\n",
+    "                idxs_exclude.append(idx)\n",
+    "    idxs_include = np.array([x for x in range(len(bs_filt_ev)) if x not in idxs_exclude])\n",
+    "\n",
+    "    print(\"%s - Scores: (AL: %.2f)\" % (session[-19:-6], (len(idxs_include)/2)/len(bs_filt_ev)))\n",
+    "    \n",
+    "    for j, lr_data in enumerate(combs):\n",
+    "        X_train, X_test, y_train, y_test = train_test_split(lr_data[idxs_include], bs_filt_ev[idxs_include], test_size=0.33)\n",
+    "        model = LogisticRegression(solver='liblinear', random_state=0).fit(X_train, y_train)\n",
+    "        print(\"%.2f: %s\" % (model.score(X_test, y_test), labels[j]))\n",
+    "        #print(confusion_matrix(y_test, model.predict(X_test)))\n",
+    "        \n",
+    "        if j == 0:\n",
+    "            scores[session]   = model.score(X_test, y_test)\n",
+    "            conf_mxs[session] = confusion_matrix(y_test, model.predict(X_test))\n",
+    "       \n",
+    "    # 5b. Use SVMs\n",
+    "    \n",
+    "    # 5c. Use GLMs\n",
+    "    #for j, lr_data in enumerate(combs):\n",
+    "    #    corr, pval, params, pvalues = get_GLM_and_prediction(lr_data, bs_filt_ev)\n",
+    "    #    print(\"%.2f, %.3f: %s\" % (corr, pval, labels[j]))\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "id": "f55c9daa",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2023-06-19_08: score 0.74, FP/FN: 1.84\n",
+      "2023-06-19_08: score 0.77, FP/FN: 1.62\n",
+      "2023-06-19_08: score 0.69, FP/FN: 2.34\n",
+      "2023-06-19_08: score 0.65, FP/FN: 1.48\n",
+      "2023-06-19_08: score 0.72, FP/FN: 1.71\n",
+      "2023-06-19_08: score 0.65, FP/FN: 2.10\n",
+      "2023-06-19_08: score 0.57, FP/FN: 2.49\n",
+      "2023-06-19_08: score 0.68, FP/FN: 1.21\n",
+      "2023-06-19_08: score 0.64, FP/FN: 3.42\n",
+      "2023-06-19_08: score 0.76, FP/FN: 1.28\n",
+      "2023-06-19_08: score 0.70, FP/FN: 1.64\n",
+      "2023-06-19_08: score 0.72, FP/FN: 1.47\n",
+      "2023-06-19_08: score 0.67, FP/FN: 2.15\n",
+      "2023-06-19_08: score 0.72, FP/FN: 1.79\n",
+      "2023-06-19_08: score 0.68, FP/FN: 1.38\n",
+      "2023-06-19_08: score 0.70, FP/FN: 4.27\n",
+      "2023-06-19_08: score 0.68, FP/FN: 1.52\n"
+     ]
+    }
+   ],
+   "source": [
+    "for k, v in conf_mxs.items():\n",
+    "    print(\"%s: score %.2f, FP/FN: %.2f\" % (session[-19:-6], scores[k], v[0][1]/v[1][0]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4be67fb1",
+   "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
+}

analysis/PSTH/UMAP tSNE on W1 - W4.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "id": "4001a375",
    "metadata": {},
    "outputs": [],
@@ -15,22 +15,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "id": "570920e2",
    "metadata": {},
    "outputs": [
     {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "2024-03-20 14:22:56.896776: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA\n",
-      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
-      "2024-03-20 14:22:57.022611: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory\n",
-      "2024-03-20 14:22:57.022631: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n",
-      "2024-03-20 14:22:57.044694: E tensorflow/stream_executor/cuda/cuda_blas.cc:2981] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered\n",
-      "2024-03-20 14:22:57.630050: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer.so.7'; dlerror: libnvinfer.so.7: cannot open shared object file: No such file or directory\n",
-      "2024-03-20 14:22:57.630115: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libnvinfer_plugin.so.7'; dlerror: libnvinfer_plugin.so.7: cannot open shared object file: No such file or directory\n",
-      "2024-03-20 14:22:57.630122: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Cannot dlopen some TensorRT libraries. If you would like to use Nvidia GPU with TensorRT, please make sure the missing libraries mentioned above are installed properly.\n"
+     "ename": "ModuleNotFoundError",
+     "evalue": "No module named 'behavior'",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
+      "\u001b[0;32m/tmp/ipykernel_1646657/2689562968.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mdecomposition\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0msklearn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcluster\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mDBSCAN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mbehavior\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mget_extent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mpopulation\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0munit_response_matrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivity_at_phase\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'behavior'"
      ]
     }
    ],

analysis/Target/Speed after successful target.ipynb

@@ -156,7 +156,6 @@
     "    dx = np.sqrt(np.square(np.diff(moseq[:, 3])) + np.square(np.diff(moseq[:, 4])))\n",
     "    dt = np.diff(moseq[:, 0])\n",
     "    speed = np.concatenate([dx/dt, [dx[-1]/dt[-1]]])\n",
-    "\n",
     "    speed_smooth = np.convolve(speed, kernel, 'same') / kernel.sum()\n",
     "\n",
     "    ax = axes[i]\n",