examples_detection.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import rippl_AI\n",
    "import importlib\n",
    "importlib.reload(rippl_AI)\n",
    "import aux_fcn\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib qt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic detection example\n",
    "In this section, a use example of the predict and detection functions are provided"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data already exists. Moving on.\n"
     ]
    }
   ],
   "source": [
    "from figshare.figshare import Figshare\n",
    "fshare = Figshare()\n",
    "\n",
    "article_ids = [14959449] \n",
    "sess=['Dlx1']                                  \n",
    "for id,s in zip(article_ids,sess):\n",
    "    datapath = os.path.join('Downloaded_data', f'{s}')\n",
    "    if os.path.isdir(datapath):\n",
    "        print(\"Data already exists. Moving on.\")\n",
    "    else:\n",
    "        print(\"Downloading data... Please wait, this might take up some time\")        # Can take up to 10 minutes\n",
    "        fshare.retrieve_files_from_article(id,directory=datapath)\n",
    "        print(\"Data downloaded!\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Channel map value:  [0, 1, 2, 3, 4, 5, 6, 7]\n",
      "Downloaded_data\\Dlx1\\figshare_14959449/lfp_Dlx1-2021-02-12_12-46-54.dat\n",
      "fileStart  0\n",
      "fileStop  490242048\n",
      "nSamples  245121024\n",
      "nSamplesPerChannel  30640128\n",
      "nSamplesPerChunk  10000\n",
      "size data  30640128\n",
      "Sampling frequency:  30000\n",
      "Shape of the original data (30640128, 8)\n",
      "[[1 1 2]\n",
      " [2 1 2]\n",
      " [3 1 2]\n",
      " [4 1 2]\n",
      " [5 1 2]\n",
      " [6 1 2]\n",
      " [7 1 2]\n",
      " [8 1 2]]\n"
     ]
    }
   ],
   "source": [
    "path=os.path.join('Downloaded_data','Dlx1','figshare_14959449')\n",
    "\n",
    "sf, expName, ref_channels, dead_channels = aux_fcn.load_info(path)\n",
    "channels_map = aux_fcn.load_channels_map(path)\n",
    "\n",
    "# Reformat channels into correct values\n",
    "channels, shanks, ref_channels = aux_fcn.reformat_channels(channels_map, ref_channels)\n",
    "# Read .dat\n",
    "print('Channel map value: ',channels)\n",
    "LFP = aux_fcn.load_raw_data(path, expName, channels, verbose=True)\n",
    "print('Sampling frequency: ', sf)\n",
    "print('Shape of the original data',LFP.shape)\n",
    "print(channels_map)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The predict function takes care of normalizing and subsampling your data\n",
    "# If no architecture or model is specified, the best CNN1D will be used\n",
    "prob,LFP_norm=rippl_AI.predict(LFP,sf) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An interactive GUI will be displayed, choose your deired threshold\n",
    "det_ind=rippl_AI.get_intervals(prob,LFP_norm=LFP_norm) \n",
    "print(f\"{det_ind.shape[0]} events where detected\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **manual event curation GUI** is included in the toolbox. You can select and discard events. As a result of the curation a .txt file will be generated with the events' times.\n",
    "For a detailed description of its features, parameters and functionalities, check the comments of *manual_curation* in *aux_fcn.py*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aux_fcn.manual_curation(det_ind,LFP_norm,file_path='C:\\\\rippl-AI\\\\Curated_txt')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get performances after detection\n",
    "Every model predict, get_intervals is used automatically and the performance metric is ploted"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data download"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from figshare.figshare import Figshare\n",
    "fshare = Figshare()\n",
    "\n",
    "article_ids = [14959449] \n",
    "sess=['Dlx1']                                    \n",
    "for id,s in zip(article_ids,sess):\n",
    "    datapath = os.path.join('Downloaded_data', f'{s}')\n",
    "    if os.path.isdir(datapath):\n",
    "        print(\"Data already exists. Moving on.\")\n",
    "    else:\n",
    "        print(\"Downloading data... Please wait, this might take up some time\")        # Can take up to 10 minutes\n",
    "        fshare.retrieve_files_from_article(id,directory=datapath)\n",
    "        print(\"Data downloaded!\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path=os.path.join('Downloaded_data','Dlx1','figshare_14959449')\n",
    "\n",
    "sf, expName, ref_channels, dead_channels = aux_fcn.load_info(path)\n",
    "\n",
    "channels_map = aux_fcn.load_channels_map(path)\n",
    "# Now the ground truth (GT) tagged events is loaded \n",
    "ripples=aux_fcn.load_ripples(path)/sf\n",
    "# Reformat channels into correct values\n",
    "channels, shanks, ref_channels = aux_fcn.reformat_channels(channels_map, ref_channels)\n",
    "# Read .dat\n",
    "print('Channel map value: ',channels)\n",
    "LFP = aux_fcn.load_raw_data(path, expName, channels, verbose=True)\n",
    "print('Sampling frequency: ', sf)\n",
    "print('Shape of the original data',LFP.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Two loops going over every possible model\n",
    "architectures=['XGBOOST','SVM','LSTM','CNN1D','CNN2D']\n",
    "SWR_prob=[[None]*5]*5\n",
    "for i,architecture in enumerate(architectures):\n",
    "    print(i,architecture)\n",
    "    for n in range(1,6):\n",
    "        # Make sure the selected model expected number of channels is the same as the channels array passed to the predict fcn\n",
    "        # In this case, we are manually setting the channel array to 3 \n",
    "        if architecture=='CNN2D' and n>=3:\n",
    "            channels=[0,3,7]\n",
    "        else:\n",
    "            channels=[0,1,2,3,4,5,6,7]\n",
    "        SWR_prob[i][n-1],_=rippl_AI.predict(LFP,sf,arch=architecture,model_number=n,channels=channels)\n",
    "\n",
    "# SWR_prob contains the output of each model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "th_arr=np.linspace(0.1,1,10)\n",
    "fig,axs=plt.subplots(5,5,figsize=(10,10),sharex='all',sharey='all')\n",
    "for i in range(5):\n",
    "    for j in range(5):\n",
    "        F1_arr=np.zeros(shape=(len(th_arr)))\n",
    "        for k,th in enumerate(th_arr):\n",
    "            det_ind=rippl_AI.get_intervals(SWR_prob[i][j],threshold=th)\n",
    "            #print(ripples)\n",
    "            _,_,F1_arr[k],_,_,_=aux_fcn.get_performance(det_ind,ripples)\n",
    "        axs[i,j].plot(th_arr,F1_arr)\n",
    "    axs[i,0].set_title(architectures[i])\n",
    "\n",
    "axs[0,0].set_xlabel('Threshold')\n",
    "axs[0,0].set_ylabel('F1')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detecting with less than 8 channels\n",
    "Detectors need several channels for optimal performance. We found out that 8 channels have enough information to assure good performance. But what happens if we don't have 8? We have seen that interpolating the missing channels also works. In this section, we will show how to use the interpolation function we have created for this purpose, inside the `aux_fcn` package."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data download"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from figshare.figshare.figshare import Figshare\n",
    "fshare = Figshare()\n",
    "\n",
    "article_ids = [14959449] \n",
    "sess=['Dlx1']                                    \n",
    "for id,s in zip(article_ids,sess):\n",
    "    datapath = os.path.join('Downloaded_data', f'{s}')\n",
    "    if os.path.isdir(datapath):\n",
    "        print(\"Data already exists. Moving on.\")\n",
    "    else:\n",
    "        print(\"Downloading data... Please wait, this might take up some time\")        # Can take up to 10 minutes\n",
    "        fshare.retrieve_files_from_article(id,directory=datapath)\n",
    "        print(\"Data downloaded!\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data load\n",
    "To ilustrate how 'interpolate_channels' can be used to extract the desired number of channels, we will be simulating two cases using the DLx1 session:\n",
    "1. We are using a recording probe that extracts 4 channels, when we need 8.\n",
    "2. Some channels are dead or have to much noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path=os.path.join('Downloaded_data','Dlx1','figshare_14959449')\n",
    "\n",
    "sf, expName, ref_channels, dead_channels = aux_fcn.load_info(path)\n",
    "\n",
    "channels_map = aux_fcn.load_channels_map(path)\n",
    "channels, shanks, ref_channels = aux_fcn.reformat_channels(channels_map, ref_channels)\n",
    "LFP = aux_fcn.load_raw_data(path, expName, channels, verbose=False)\n",
    "print('Sampling frequency: ', sf)\n",
    "print('Shape of the original data',LFP.shape)\n",
    "LFP_linear=LFP[:,[0,2,4,6]]\n",
    "print('Shape of the 4 channels simulated data: ',LFP_linear.shape)\n",
    "LFP[:,[2,5]]=0\n",
    "LFP_dead=LFP\n",
    "print('Sample of the simulated dead LFP: ',LFP_dead[0])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After interpolation, the data is ready to use in prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define channels\n",
    "channels_interpolation = [0,-1,1,-1,2,-1,-1,3]\n",
    "\n",
    "# Make interpolation\n",
    "LFP_interpolated = aux_fcn.interpolate_channels(LFP_linear, channels_interpolation)\n",
    "print('Shape of the interpolated LFP: ',LFP_interpolated.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " # Define channels\n",
    "channels_interpolation = [0,1,-1,3,4,-1,6,7]\n",
    "\n",
    "# Make interpolation\n",
    "LFP_interpolated = aux_fcn.interpolate_channels(LFP_dead, channels_interpolation)\n",
    "print('Value of the 1st sample of the interpolated LFP: ',LFP_interpolated[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using an ensemble model\n",
    "In this section, we will show how to use an ensemble model that combines the output of the best models of each architecture. This model has better performance and more stability than the individual models. In this case, only the best ensemble model will be provided.\n",
    "\n",
    "First, the output of the 5 selected models needs to reshaped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5 outputs are generated\n",
    "architectures=['XGBOOST','SVM','LSTM','CNN1D','CNN2D']\n",
    "output=[]\n",
    "for architecture in architectures:\n",
    "    channels=[0,1,2,3,4,5,6,7]\n",
    "    SWR_prob,_=rippl_AI.predict(LFP,sf,arch=architecture,model_number=1,channels=channels)\n",
    "    output.append(SWR_prob)\n",
    "ens_input=np.array(output).transpose()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generating ensemble model output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "prob_ens=rippl_AI.predict_ens(ens_input)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots()\n",
    "th_arr=np.linspace(0.1,1,10)\n",
    "F1_arr=np.zeros(shape=(len(th_arr)))\n",
    "for k,th in enumerate(th_arr):\n",
    "    det_ind=rippl_AI.get_intervals(prob_ens,threshold=th)\n",
    "    _,_,F1_arr[k],_,_,_=aux_fcn.get_performance(det_ind,ripples)\n",
    "ax.plot(th_arr,F1_arr)\n",
    "ax.set_title('Ensemble model')\n",
    "ax.set_ylim(-0.05,0.8)\n",
    "ax.set_xlabel('Threshold')\n",
    "ax.set_ylabel('F1')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PublicBCG_d",
   "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.9.15"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}