bpc_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basis profile curve identification to understand electrical stimulation effects in human brain networks\n",
    "\n",
    "Method developed in Matlab by Kai J. Miller, available at Miller KJ, Mueller KR, Hermes D. Basis profile curve identification to understand electrical stimulation effects in human brain networks. doi: https://doi.org/10.1101/2021.01.24.428020\n",
    "\n",
    "This project was supported by the National Institute Of Mental Health of the National Institutes of Health under Award Number R01MH122258. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health\n",
    "\n",
    "This Jupyter Notebook was written by Tal Pal Attia, Harvey Huang, Max van den Boom and Dora Hermes (2021)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview\n",
    "\n",
    "This notebook is a walk through of the methods presented in [Miller et al. (2021)](https://www.biorxiv.org/content/10.1101/2021.01.24.428020v1).\n",
    "\n",
    "> **Abstract (Miller et al. 2021)**\n",
    "<br /> Brain networks can be explored by delivering brief pulses of electrical current in one area while measuring voltage responses in other areas. We propose a convergent paradigm to study brain dynamics, focusing on a single brain site to observe the average effect of stimulating each of many other brain sites. Viewed in this manner, visually-apparent motifs in the temporal response shape emerge from adjacent stimulation sites. This work constructs and illustrates a data-driven approach to determine characteristic spatiotemporal structure in these response shapes, summarized by a set of unique “basis profile curves” (BPCs). Each BPC may be mapped back to underlying anatomy in a natural way, quantifying projection strength from each stimulation site using simple metrics. Our technique is demonstrated for an array of implanted brain surface electrodes in a human patient. This framework enables straightforward interpretation of single-pulse brain stimulation data, and can be applied generically to explore the diverse milieu of interactions that comprise the connectome.\n",
    "\n",
    "This notebook will walk you through the following five steps to compute and visualize Basis Profile Curves (BPCs) in an interactive manner:\n",
    "1. Getting started with the Python environment and packages\n",
    "2. Loading and checking the BIDS data and metadata\n",
    "3. Preprocessing the data\n",
    "4. Calculating BPCs \n",
    "5. Visualizing BPCs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Getting started with the Python environment and packages\n",
    "Before getting started, you need to have the necessary software installed and data downloaded. \n",
    "\n",
    "An iEEG dataset formatted according to the Brain Imaging Data Structure used this tutorial can be downloaded from OpenNeuro:\n",
    "https://openneuro.org/datasets/ds003708\n",
    "\n",
    "## 1.1 Setting up the environment\n",
    "This notebook requires the following Python packages to be installed: [bpc-requirements.txt](./bpc-requirements.txt)\n",
	"To ensure all requirements are installed run 'pip install -r bpc-requirements.txt'\n",
    "\n",
    "\n",
    "To easy setup an enviroment on a Mac or Linux use this file:  [bpc-dev-env.sh](./bpc-dev-env.sh) \\\n",
    "Run from terminal: <code>bash bpc-dev-env.sh</code> \\\n",
    "or from this notebook: <code>!bash bpc-dev-env.sh</code>\n",
    "\n",
    "You can run the following cell to setup your environment for this tutorial, after running, refresh the page and make sure the bpc kernel is selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating the enviroment for this toturial and refreshing the notebook\n",
    "# !bash bpc-dev-env.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Importing python packages\n",
    "\n",
    "We use several Python packages, and the following allow us to work with the BIDS dataset in this notebook:\n",
    " - [PyBIDS](https://bids-standard.github.io/pybids/#) ([Yarkoni et al., 2019](https://joss.theoj.org/papers/10.21105/joss.01294)), a Python library to centralize interactions with BIDS datasets. For more information on BIDS, see: [https://bids.neuroimaging.io](https://bids.neuroimaging.io).\n",
    " - [IeegPrep](https://pypi.org/project/ieegprep/), a python library to read, pre-process and epoch Intracranial Electroencephalography (iEEG) data that is structured according to the Brain Imaging Data Structure (BIDS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# packages importing\n",
    "import os\n",
    "\n",
    "# for scientific computing and data visualization \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ieegprep import IeegDataReader\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from sklearn.decomposition import NMF\n",
    "import itertools\n",
    "from tqdm import tqdm\n",
    "\n",
    "# for handling neuroimaging data\n",
    "import bids\n",
    "from nilearn import plotting\n",
    "\n",
    "# for this tutorial\n",
    "import functions.helperFunctions as helper\n",
    "import functions.pyBPCs as pyBPCs\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Loading and checking the BIDS data and metadata\n",
    "\n",
    "An example CCEP dataset is available on OpenNeuro ([link](https://openneuro.org/datasets/ds003708)). This dataset is formatted according to the Brain Imaging Data Structure ([BIDS](https://bids.neuroimaging.io/)) and contains one subject to work with in this tutorial. \n",
    "\n",
    "This dataset includes an electrocorticography (ECoG) dataset with single pulse stimulation, and accompanying metadata, such as electrode positions, channel information, stimulation events etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change this path to the full hardcoded BIDS data path to be analyzed\n",
    "BIDS_dataset_path = '' # '/full/path/to/Basis_profile_curve/data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Load BIDS metadata\n",
    "\n",
    "We use pyBIDS to initialize a BIDSLayout: this will index the files and metadata under the specified root folder (should be 1 subject, 1 session, 1 run)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the layout\n",
    "layout = bids.BIDSLayout(BIDS_dataset_path)\n",
    "\n",
    "# Print some basic information about the layout\n",
    "print(layout)\n",
    "print('Subjects in this BIDS layout:', layout.get_subjects())\n",
    "\n",
    "all_files = layout.get()\n",
    "df_layout = layout.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Select run to analyze\n",
    "\n",
    "User input regarding the subject, session and task to be analyzed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change this to make sure it fits your data\n",
    "bids_sub = '01' # The subject label\n",
    "bids_ses = 'ieeg01' # The session label\n",
    "bids_task = 'ccep' # The task name\n",
    "bids_run = '01' # The run name\n",
    "bids_space = 'MNI152NLin6Sym' # The electrodes space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can load filenames and metadata files needed for this analysis from the BIDSLayout:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve file names\n",
    "\n",
    "# pyBIDS does not yet support MEF3 data, we use the OS module.\n",
    "for (root, dirnames, _) in os.walk(BIDS_dataset_path):\n",
    "    for directory in dirnames:\n",
    "        if directory.endswith(\".mefd\"):\n",
    "             mefName = os.path.join(root, directory)\n",
    "print('We will load the following iEEG data:',mefName)\n",
    "                \n",
    "dataJson = layout.get(subject=bids_sub, session=bids_ses, run=bids_run, task=bids_task, suffix='ieeg', extension=\"json\")[0]\n",
    "# print('\\nSidecar JSON (*_ieeg.json):', dataJson)\n",
    "\n",
    "channels_tsv_name = layout.get(subject=bids_sub, session=bids_ses, run=bids_run, task=bids_task, suffix='channels', extension = \"tsv\")[0]\n",
    "# print('\\nChannels description (*_channels.tsv):', channels_tsv_name)\n",
    "\n",
    "events_tsv_name = layout.get(subject=bids_sub, session=bids_ses, run=bids_run, task=bids_task, suffix='events', extension = \"tsv\")[0]\n",
    "# print('\\nevents_tsv_name:', events_tsv_name)\n",
    "\n",
    "electrodes_tsv_name = layout.get(subject=bids_sub, session=bids_ses, suffix='electrodes', space = bids_space, extension = \"tsv\")[0]\n",
    "# print('\\nElectrode description (*_electrodes.tsv):', electrodes_tsv_name)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the pandas.DataFrame.head function to quickly test if our files have the expected data in them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_channels = channels_tsv_name.get_df()  # Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "# np.shape(df_channels)\n",
    "df_channels.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_events = events_tsv_name.get_df()  # Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "df_events.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_electrodes = electrodes_tsv_name.get_df()  # Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "df_electrodes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Render brain surface\n",
    "\n",
    "Render electrode labels in MNI space. This should open a separate browser window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_colors=[]\n",
    "\n",
    "xyz_list = df_electrodes[['x', 'y','z']].values.tolist()\n",
    "electrodeName_list = df_electrodes[['name']].name.values.tolist()\n",
    "\n",
    "# gray out bad channels\n",
    "for cords_name in electrodeName_list:\n",
    "    if df_channels[df_channels.name == cords_name].status.values[0] == 'good':\n",
    "        marker_colors.append('red')\n",
    "    else:\n",
    "        marker_colors.append('gray')\n",
    "\n",
    "view = plotting.view_markers(marker_coords=xyz_list, marker_labels=electrodeName_list, marker_size=7.5,\n",
    "                             marker_color=marker_colors, title='Electrodes in MNI space rendered on a standard brain') # Insert a 3d plot of markers in a brain\n",
    "view.open_in_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Preprocessing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the rendering with electrode labels and define the electrode of interest (e.g LMS2) and time window to be analyzed (e.g. -1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define electrode of interest\n",
    "el_interest = 'LMS2'\n",
    "# epoch_size for visualization\n",
    "epoch_limits = [-1,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Extracting relevant information from metadata\n",
    "Because we have the data in BIDS, we can automatically pull information essential for analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get electrode of interest from channels.tsv \n",
    "el_interest_nr = df_channels.loc[df_channels['name'] == el_interest]\n",
    "\n",
    "# Sampling rate from sidecar JSON (*_ieeg.json).\n",
    "dict_dataJson = dataJson.get_dict()\n",
    "srate = dict_dataJson[\"SamplingFrequency\"]\n",
    "\n",
    "# Construct a time vector\n",
    "tt = pyBPCs.create_time_vector(epoch_limits,srate)\n",
    "\n",
    "# Number of trials with electrical stimulation, good and not stimulating electrode of interest\n",
    "df_include_trials = df_events.loc[(df_events[\"trial_type\"]=='electrical_stimulation') &\n",
    "                               (df_events[\"status\"]=='good') &\n",
    "                               (~df_events[\"electrical_stimulation_site\"].str.contains(el_interest))]\n",
    "df_include_trials.head()\n",
    "\n",
    "# Filter channels to select good channels, only ECOG and SEEG\n",
    "df_good_channels = df_channels[(df_channels[\"status\"] == 'good') &\n",
    "                               ((df_channels[\"type\"] == 'ECOG') | (df_channels[\"type\"] == 'SEEG'))] \n",
    "# df_good_channels.head()\n",
    "\n",
    "reader = IeegDataReader(mefName, preload_data=False)\n",
    "\n",
    "print('Sampling frequency is', srate, 'Hz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Calculate convergent Matrix V\n",
    "\n",
    "**Start running the following cell because loading all the data and re-referencing to get CCEPs takes a few minutes.**\n",
    "\n",
    "In the meantime, we can discuss the convergent paradigm (A) used to calculate BPCs and contrast with the divergent paradigm (B) and the interpretation of the convergent paradign (E). \n",
    "     \n",
    "    Note: Rereferencing is tailored to these CCEP data, using the custom function ccep_CAR64blocks() which applies common average referencing on CCEP data. This excludes channels with high variance from 500-2000 ms after stimulation onset or high variance from 10-100 ms after stimulation onset (channels with stimulation artifact or large evoked responses). The function also assumes that noise is shared between blocks of 64 channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "We iterate through all trials with electrical stimulation, no artifacts and not stimulating the recorded electrode of interest.\n",
    "\"\"\"\n",
    "\n",
    "df_cceps = pd.DataFrame(columns = ['ccep_name_1', 'ccep_num_1', 'ccep_status_1', 'ccep_name_2', 'ccep_num_2', 'ccep_status_2'])\n",
    "V_pre = []  # single-trial stimulation-evoked voltage matrix\n",
    "\n",
    "for index, row in tqdm(df_include_trials.iloc[:].iterrows(), desc='Included trials', total=df_include_trials.shape[0]):\n",
    "    \n",
    "    df_data = pd.DataFrame(columns = ['channel', 'data']) # dataframe for mef3 data\n",
    "    samples = [row.onset * srate + element * srate for element in epoch_limits]  # get start and end samples\n",
    "    # Read mef3 data\n",
    "    for channel in df_electrodes.name:\n",
    "        data = reader.retrieve_sample_range_data(sample_start=samples[0], sample_end=samples[1], channels=channel)[0]  # Returns 1D numpy array with data from sample X to sample Y\n",
    "        curr_channel = pd.DataFrame({'channel': channel, 'data': [data]})\n",
    "        df_data = pd.concat([df_data, curr_channel], ignore_index=True)  # append to mef3 data dataframe\n",
    "    \n",
    "    # Get stimulated channels\n",
    "    trial_ccep_names = row.electrical_stimulation_site.split(\"-\")\n",
    "    curr_ccep = pd.DataFrame({ 'ccep_name_1': [trial_ccep_names[0]],\n",
    "                               'ccep_num_1': [df_channels.index[df_channels['name'] == trial_ccep_names[0]].tolist()[0]],\n",
    "                               'ccep_status_1': [df_channels['status'][df_channels.index[df_channels['name'] == trial_ccep_names[0]].tolist()[0]]],\n",
    "                               'ccep_name_2': [trial_ccep_names[1]],\n",
    "                               'ccep_num_2': [df_channels.index[df_channels['name'] == trial_ccep_names[1]].tolist()[0]],\n",
    "                               'ccep_status_2': [df_channels['status'][df_channels.index[df_channels['name'] == trial_ccep_names[1]].tolist()[0]]]})\n",
    "    df_cceps = pd.concat([df_cceps, curr_ccep], ignore_index=True)\n",
    "\n",
    "    # Delete current electrical stimulation site from dataFrame\n",
    "    curr_indices = df_good_channels[df_good_channels['name'].isin([trial_ccep_names[0], trial_ccep_names[1]])].index\n",
    "    df_incl_channels = df_good_channels[~df_good_channels.index.isin(curr_indices)]\n",
    "    \n",
    "    # Run the custom function ccep_CAR64blocks() in functions/pyBPCs.py\n",
    "    df_data_CAR = pyBPCs.ccep_CAR64blocks(df_data, tt, df_incl_channels)\n",
    "    \n",
    "    # append to single-trial stimulation-evoked voltage matrix\n",
    "    V_pre.append(df_data_CAR[df_data_CAR['channel'] == el_interest]['data'].values[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Plot convergent matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Calculate baseline subtraction and plot difference\n",
    "\"\"\"\n",
    "V_pre = np.array(V_pre)\n",
    "\n",
    "# calculate baseline\n",
    "baseline_V = np.mean(V_pre[:,np.logical_and(tt > -.5,tt < -0.05)],1)\n",
    "# make baseline into the same size\n",
    "baseline_Vrep = np.transpose(np.tile(baseline_V,(np.shape(tt)[0],1)))\n",
    "\n",
    "# subtract baseline from V_pre\n",
    "V_pre_baseSub = np.subtract(V_pre,baseline_Vrep)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 10), dpi=80)\n",
    "\n",
    "fig.suptitle('Convergent Matrix (V)',fontsize= 15)\n",
    "vmin=-250\n",
    "vmax=250\n",
    "im = ax1.imshow(V_pre, aspect='auto',extent=helper.extents(tt) + helper.extents(df_include_trials.index),vmin=vmin, vmax=vmax)\n",
    "ax1.set_xlabel('Time from stimulation (s)')\n",
    "ax1.set_ylabel('Stimulation trial index (k)')\n",
    "ax1.set_title('Convergent Matrix (V)')\n",
    "im = ax2.imshow(V_pre_baseSub, aspect='auto',extent=helper.extents(tt) + helper.extents(df_include_trials.index),vmin=vmin, vmax=vmax)\n",
    "ax2.set_xlabel('Time from stimulation (s)')\n",
    "ax2.set_ylabel('Stimulation trial index (k)')\n",
    "ax2.set_title('Convergent Matrix (V) after baseline correction')\n",
    "\n",
    "fig.subplots_adjust(right=0.8)\n",
    "cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "fig.colorbar(im, cax=cbar_ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Calculate BPCs\n",
    "\n",
    "Identifying basis profile curves (BPCs) that group characteristic shapes in the convergent CCEPs. \n",
    "\n",
    "The input for the BPC calculation is the convergent matrix (V_pre, numpy array), the time vector (tt, numpy array), and the pandas dataframe pair_types that we will generate in this section.\n",
    "\n",
    "## 4.1 Select time-frame for BPC extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select the epoch times to enter in the BPC analyses in seconds\n",
    "BPCs_epoch = [0.015, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Calculate the significance matrix\n",
    "\n",
    "To calculate the significance matrix, we project the unit-normalized stimulation trials into all other trials. We then calculate t-values across all subgroups of stimulation pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# V: single-trial stimulation-evoked voltage matrix (Dimensions of V are with K total stimulation events by T total timepoints)\n",
    "# tt_BPCs: BPCs time vector\n",
    "V, tt_BPCs = pyBPCs.bpcVoltage(V_pre_baseSub, tt, BPCs_epoch)  # take only Voltage at BPC times\n",
    "\n",
    "# You can compare results to taking the V_pre matrix before baseline correction, just to test\n",
    "# V, tt_BPCs = pyBPCs.bpcVoltage(V_pre, tt, BPCs_epoch)  # take only Voltage at BPC times\n",
    "\n",
    "# pair_types: structure with subgroups of stimulation pairs, with fields\n",
    "#    ccep_name_1, ccep_name_2: electrodes associated with the pair, cathode-anode ordering\n",
    "#    indices: indices of subgroup in V matrix (out of K total stims)\n",
    "\n",
    "pair_types = df_cceps.groupby(['ccep_name_1', 'ccep_num_1', 'ccep_name_2', 'ccep_num_2']).size().reset_index().rename(columns={0: 'count'})\n",
    "pair_types = pair_types.sort_values(by=['ccep_num_1']) \n",
    "\n",
    "indices_column = []\n",
    "for index, row in pair_types.iterrows():\n",
    "    indices_column.append(df_cceps.index[(df_cceps['ccep_name_1'] == row['ccep_name_1']) & (df_cceps['ccep_name_2'] == row['ccep_name_2'])].tolist())\n",
    "pair_types['indices'] = indices_column\n",
    "\n",
    "pair_types = pair_types.reset_index(drop=True)\n",
    "# pair_types.head()\n",
    "\n",
    "### calculate\n",
    "V0 = V/(np.ones((V.shape[0], 1)) * (np.sum(V ** 2, axis=0) ** 0.5))  # normalize (L2 norm) each trial\n",
    "P = V0.T @ V  # calculate internal projections\n",
    "\n",
    "tmat = pyBPCs.nativeNormalized(pair_types, P)\n",
    "\n",
    "# visualize the set of native normalized single stimulation cross-projections \n",
    "plt.figure(figsize=(10, 8), dpi=80)\n",
    "plt.imshow(tmat,aspect='auto') # Plot a matrix or an array as an image.\n",
    "plt.clim(0, 10)\n",
    "plt.ylabel('stimulation pair subgroup')\n",
    "plt.xlabel('stimulation pair subgroup')\n",
    "plt.title('Significance matrix '+r'$\\Xi$',fontsize= 15)\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 Iteratively decrease inner components of non-negative matrix factorization\n",
    "\n",
    "Using Non-Negative Matrix Factorization (NMF) on the significance matrix to cluster sites that produce similar measured responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Find two non-negative matrices (W, H) whose product approximates the factorize, \n",
    "non-negative and rescaled matrix (t0_scaled). \n",
    "At the end, the matrix H has size of number of clusters by stimulation pair sub-groups.\n",
    "\"\"\"\n",
    "\n",
    "#### Default settings and initialization\n",
    "\n",
    "# number of cluster dimensions\n",
    "cl_dim = 10\n",
    "    \n",
    "# number of iterations to re-run NNMF - because can get suboptimal factorization due to non-unique nature of NNMF (starts out with a random matrix)\n",
    "num_reruns = 20\n",
    "\n",
    "# convergence threshold: 1e-5\n",
    "conv_thresh = 0.00001\n",
    "\n",
    "t0 = tmat\n",
    "t0[t0 < 0] = 0 # factorize t-matrix, but first make non-negative\n",
    "\n",
    "# Globally rescale data to avoid potential overflow/underflow\n",
    "t0_scaled = t0 / (np.max(t0))    \n",
    "\n",
    "nnmf_xcorr_score = 100  # start off-diagonal penalty score out with high value\n",
    "\n",
    "WH_struct = pd.DataFrame(columns=['W', 'H', 'nnmf_xcorr_score'])  # saving structure\n",
    "\n",
    "#### Do the NMF\n",
    "\n",
    "#while nnmf_xcorr_score >. 5:\n",
    "while nnmf_xcorr_score > 1:\n",
    "  \n",
    "    cl_dim = cl_dim - 1 # reduce number of inner dimensions\n",
    "    \n",
    "    if cl_dim == 0:\n",
    "        break\n",
    "\n",
    "    # multiple run-throughs of  nnmf\n",
    "    tmp_mat_W = []\n",
    "    tmp_mat_H = []\n",
    "    tmp_err = []\n",
    "\n",
    "    for k in range(num_reruns):\n",
    "\n",
    "        model = NMF(n_components=cl_dim, init='random', solver='mu',tol=conv_thresh,max_iter=10000)\n",
    "        W = model.fit_transform(t0_scaled)\n",
    "        tmp_mat_W.insert(k, W)\n",
    "    \n",
    "        H = model.components_\n",
    "        tmp_mat_H.insert(k, H)\n",
    "\n",
    "        rec_err = model.reconstruction_err_\n",
    "        tmp_err.insert(k, rec_err)\n",
    "\n",
    "    # select factorization with smallest error\n",
    "    k_ind = np.argmin(tmp_err)    \n",
    "    W_min = tmp_mat_W[k_ind]\n",
    "    H_min = tmp_mat_H[k_ind]\n",
    "\n",
    "    # Normalize so rows of H have unit norm\n",
    "    Hnorms = np.sqrt(np.sum(H_min.T ** 2, axis=0))\n",
    "    Hnorms_ = (np.reshape(np.tile(Hnorms.T, t0_scaled.shape[1]),(-1, Hnorms.shape[0]))).T\n",
    "    H_min = H_min / Hnorms_\n",
    "\n",
    "    # Score for this decomposition - off diagonal element weights\n",
    "    HHT = H_min @ H_min.T\n",
    "    HHT_triu = np.triu(HHT, 1)\n",
    "    HHT_triu_reshaped = np.reshape(HHT_triu, (1,-1))\n",
    "    # nnmf_xcorr_score = np.max(HHT_triu_reshaped) # max of nonnegative matrix factorization component interdependencies\n",
    "    nnmf_xcorr_score = np.sum(HHT_triu_reshaped) # sum of off-diagonal nonnegative matrix factorization component interdependencies\n",
    "    print('Inner dimension: ', cl_dim, ', off diagonal score: ', nnmf_xcorr_score)\n",
    "\n",
    "    # save matrices and scores for plotting later\n",
    "    curr_WH = pd.DataFrame({'W': [W_min], 'H': [H_min], 'nnmf_xcorr_score': [nnmf_xcorr_score]})\n",
    "    WH_struct = pd.concat([WH_struct, curr_WH], ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Winner take all\n",
    "\n",
    "Each stimulation pair subgroup is assigned to one cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "winner-take-all on columns of H and then threshold by 1/(2*sqrt(length(pair_types))) \n",
    " -- since all equal would be 1/sqrt(length(pair_types))\n",
    "\"\"\"\n",
    "\n",
    "H0 = 0 * H_min\n",
    "\n",
    "for k in range(len(pair_types)):\n",
    "    k_ind = np.argmax(H_min, axis=0)[k]\n",
    "    H0[k_ind][k] = H_min[k_ind][k]\n",
    "\n",
    "H0_ = H0 > 1 / (2 * np.sqrt(len(pair_types)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.5 Identification of Basis Profile Curves \n",
    "\n",
    "BPCs are identified from the clustered groups (rows of H) using linear kernel PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "output will show the subgroup numbers clustered in each BPC\n",
    "\"\"\"\n",
    "\n",
    "B = []\n",
    "B_struct = pd.DataFrame(columns=['curve','pairs']) # saving structure\n",
    "\n",
    "for q in range(H0.shape[0]):\n",
    "    cl_pg = np.where(H0_[q]) # cluster pair groups\n",
    "    cl_inds = [] # cluster indices (e.g. concatenated trials from cluster pair groups)\n",
    "    for k in cl_pg:\n",
    "        for i in (pair_types['indices'][k]).values:\n",
    "            cl_inds.extend(i)\n",
    "\n",
    "    V_tmp = V[:, cl_inds]\n",
    "\n",
    "    E = pyBPCs.kpca(V_tmp)\n",
    "    B_tmp = E[:, 0] # basis vector is 1st PC\n",
    "\n",
    "    if np.mean(B_tmp.T @ V_tmp) < 0:\n",
    "        B_tmp = -B_tmp\n",
    "    \n",
    "    curr_B = pd.DataFrame({'curve': [B_tmp], 'pairs': [cl_pg[0]]})\n",
    "    B_struct = pd.concat([B_struct, curr_B], ignore_index=True)\n",
    "    \n",
    "# pairs not represented by any basis\n",
    "excluded_pairs = np.where(1 - (np.sum(H0_, axis=0)))\n",
    "\n",
    "\"\"\"\n",
    "Calculate statistics for each basis curve, eliminating those where there is no significant representation\n",
    "\"\"\"\n",
    "\n",
    "B_struct = pyBPCs.curvesStatistics(B_struct, V, B, pair_types)\n",
    "\n",
    "\"\"\"\n",
    "Calculate projection weights\n",
    "\"\"\"\n",
    "\n",
    "B_struct = pyBPCs.projectionWeights(B_struct)\n",
    "\n",
    "B_struct.pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Visualize BPCs\n",
    "\n",
    "## 5.1 Plot Calculated BPCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = iter(cm.tab10(np.linspace(0, 1, 10)))\n",
    "\n",
    "plt.figure(figsize=(10, 8), dpi=80)\n",
    "for q in range(len(B_struct)): #  cycle through basis curves\n",
    "    plt.plot(tt[tt_BPCs[0]:tt_BPCs[1]], B_struct.curve[q], color=next(colors), label=q)\n",
    "    \n",
    "plt.xlabel('Time from stimulation (s)')\n",
    "plt.ylabel('Normalized weight of BPCs')\n",
    "plt.title('Calculated BPCs',fontsize=15)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Spatial representation of the BPCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_colors = [] # marker_color array\n",
    "marker_sizes = [] # marker_size array\n",
    "\n",
    "# all electrodes\n",
    "xyz = df_electrodes[['x', 'y','z']]\n",
    "xyz_list = xyz.values.tolist()\n",
    "electrodeName_list = df_electrodes[['name']].name.values.tolist()\n",
    "for cords in xyz_list:\n",
    "    marker_colors.append('white')\n",
    "    marker_sizes.append(7.5)\n",
    "\n",
    "# el_interest\n",
    "xyz_el_interest = df_electrodes.loc[df_electrodes.name == el_interest_nr.name.values[0]][['x', 'y','z']]\n",
    "xyz_list.append(xyz_el_interest.values[0].tolist())\n",
    "marker_colors.append('black')\n",
    "marker_sizes.append(20)\n",
    "\n",
    "pair_types['interpolated_locs'] = ''\n",
    "\n",
    "for index, row in pair_types.iterrows():\n",
    "    xyz_1 = xyz.iloc[row['ccep_num_1']].values.tolist()\n",
    "    xyz_2 = xyz.iloc[row['ccep_num_2']].values.tolist()\n",
    "    pair_types.at[index, 'interpolated_locs'] = np.mean((xyz_1, xyz_2), axis=0).tolist()\n",
    "\n",
    "# non-significant stim pair sites    \n",
    "interpolated_locs_list = pair_types.iloc[excluded_pairs[0]].interpolated_locs.values.tolist()\n",
    "for cords in interpolated_locs_list:\n",
    "    marker_colors.append('gray')\n",
    "    marker_sizes.append(7.5)\n",
    "    \n",
    "xyz_all = [y for x in [xyz_list, interpolated_locs_list] for y in x]\n",
    "\n",
    "# plot BPCs, colored\n",
    "colors = iter(cm.tab10(np.linspace(0, 1, 10)))\n",
    "\n",
    "for q in range(B_struct.shape[0]):\n",
    "    # get electrodes for BPC\n",
    "    xyz_BPC_list = pair_types.interpolated_locs[B_struct.pairs[q]].tolist()\n",
    "    xyz_all = [y for x in [xyz_all, xyz_BPC_list] for y in x]\n",
    "    plotweights = B_struct.plotweights[q]\n",
    "    curr_color = next(colors)\n",
    "    for ind_cords in range(len(xyz_BPC_list)):\n",
    "        marker_colors.append(curr_color)\n",
    "        marker_sizes.append(B_struct.plotweights[q][ind_cords] * 25)\n",
    "\n",
    "# plot all\n",
    "view = plotting.view_markers(marker_coords=xyz_all,marker_size=marker_sizes, # marker_labels=electrodeName_list,\n",
    "                             marker_color=marker_colors, title='Spatial representation of BPCs rendered on an MNI brain') # Insert a 3d plot of markers in a brain\n",
    "view.open_in_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3 Optional parameters to change\n",
    "\n",
    "Once you have completed this, you can go back to Section 4.1 and change the time interval over which the BPCs are calculated (e.g. 0.2 - 1 sec) and look at the effects on the outputs.\n",
    "\n",
    "You can also select a different electrode in Section 3 to look at the various inputs into different regions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bpc",
   "language": "python",
   "name": "bpc"
  },
  "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.12"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}