Completed and tested calibration portion of pipeline

Author: sammytheindi

conf/base/catalog.yml

@@ -98,8 +98,16 @@ prefixed_channels:
 #       compress: 3
 #   path: data/03_primary/CC01/ccCenterOut/downsampled_signals
 #   filename_suffix: ".lz4"
-#   layer: intermediate
+#   layer: primary
 
+calibration_statistics_pkl:
+  type: PartitionedDataSet
+  dataset: 
+    type: pickle.PickleDataSet
+    backend: pickle
+  path: data/03_primary/CC01/ccCenterOut/calibration_statistics
+  filename_suffix: ".pkl"
+  layer: primary
 
 
 

conf/base/parameters.yml

@@ -1,7 +1,10 @@
 patient_id: 'CC01'
+current_experiment: 'center_out'
+current_calibration: 'calibration'
 bci_states:
   CC01:
     center_out: ['StimulusCode', 'cursorX', 'cursorY', 'Baseline', 'ResultCode']
+    calibration: ['StimulusCode']
 gain: 0.25
 data_preprocessing:
   channel_labelling:
@@ -34,4 +37,9 @@ sessions:
         '20221027': ['S08', 'S09']
         '20221104': ['S09', 'S10']
         '20221108': ['S08', 'S09']
-        '20221111': ['S07', 'S08']
+        '20221111': ['S07', 'S08']
+    calibration:
+      '20221027': ['S01', 'S02', 'S03']
+      '20221104': ['S01', 'S02']
+      '20221108': ['S01', 'S02']
+      '20221111': ['S01', 'S02']

notebooks/Data Ingestion Notebook.ipynb

@@ -6,16 +6,25 @@
 from scripts.convert_bci_to_hdf5 import convert_bcistream, convert_dat
 
 
-def generate_center_out_hdf5_dataset(bcistreams, selected_sessions, patient_id):
-    sessions_dict = selected_sessions[patient_id]['center_out']
+def generate_hdf5_dataset(bcistreams, selected_sessions, patient_id, current_experiment, current_run_type):
+    sessions_dict = selected_sessions[patient_id][current_experiment]
+    calibration_dict = selected_sessions[patient_id]['calibration']
 
     for partition_key, partition_load_func in bcistreams.items():
-        
         continue_loop = True
         for paradigm_key, date_session_dict in sessions_dict.items():
             for date_key, sessions_list in date_session_dict.items():
-                if date_key in partition_key and partition_key.split('_')[-1] in sessions_list:
-                    continue_loop = False
+                if current_run_type == 'calibration':
+                    if date_key in partition_key:
+                        calibration_sessions_list = calibration_dict.get(date_key, [])
+                        if len(calibration_sessions_list):
+                            if partition_key.split('_')[-1] in calibration_sessions_list:
+                                continue_loop = False
+                        else:
+                            continue_loop = False
+                else:
+                    if date_key in partition_key and partition_key.split('_')[-1] in sessions_list:
+                        continue_loop = False
 
         if continue_loop:
             continue
@@ -29,26 +38,3 @@ def generate_center_out_hdf5_dataset(bcistreams, selected_sessions, patient_id):
         convert_dat(partition_data, h5filename=filename, add_everything=True)
 
     return {}
-
-def generate_calibration_hdf5_dataset(bcistreams, selected_sessions, patient_id):
-    sessions_dict = selected_sessions[patient_id]['center_out']
-
-    for partition_key, partition_load_func in bcistreams.items():
-        continue_loop = True
-        for paradigm_key, date_session_dict in sessions_dict.items():
-            for date_key, sessions_list in date_session_dict.items():
-                if date_key in partition_key:
-                    continue_loop = False
-
-        if continue_loop:
-            continue
-
-        partition_data = partition_load_func()
-
-        filename = partition_data.filename.replace(".dat", ".hdf5")
-
-        print(filename)
-
-        convert_dat(partition_data, h5filename=filename, add_everything=True)
-
-    return {}

src/decoding_pipeline/pipelines/data_generation/nodes.py

@@ -3,22 +3,22 @@
 generated using Kedro 0.18.3
 """
 
-from .nodes import generate_center_out_hdf5_dataset, generate_calibration_hdf5_dataset
+from .nodes import generate_hdf5_dataset
 
 from kedro.pipeline import Pipeline, node
 
 def create_pipeline(**kwargs) -> Pipeline:
     return Pipeline([
         node(
-            func=generate_center_out_hdf5_dataset,
-            inputs=["center_out_dat", "params:sessions", "params:patient_id"],
+            func=generate_hdf5_dataset,
+            inputs=["center_out_dat", "params:sessions", "params:patient_id", "params:current_experiment", "params:current_experiment"],
             outputs="center_out_hdf5",
             name="convert_center_out_dat_to_hdf5_node",
         ),
         node(
-            func=generate_calibration_hdf5_dataset,
-            inputs=["calibration_dat", "params:sessions", "params:patient_id"],
+            func=generate_hdf5_dataset,
+            inputs=["calibration_dat", "params:sessions", "params:patient_id", "params:current_experiment", "params:current_calibration"],
             outputs="calibration_hdf5",
-            name="generate_calibration_dat_to_hdf5_node"
+            name="convert_calibration_dat_to_hdf5_node"
         )
     ])

src/decoding_pipeline/pipelines/data_generation/pipeline.py

@@ -94,14 +94,14 @@ def prefix_single_channel_info(channels, patient_id, grid_split):
                 'ch_suffix_order': ch_suffix_order 
             }
 
-def extract_bci_data(h5_data, selected_channels, electrode_labels, states, patient_id, gain):
+def extract_bci_data(h5_data, selected_channels, electrode_labels, states, patient_id, gain, current_experiment):
     eeglabels = electrode_labels['eeglabels']
     auxlabels = electrode_labels['auxlabels']
 
     ch_include = selected_channels['ch_include']
     ch_exclude = selected_channels['ch_exclude']
 
-    selected_states = states[patient_id]['center_out']
+    selected_states = states[patient_id][current_experiment]
 
     num_channels = len(eeglabels)
 
@@ -166,8 +166,8 @@ def extract_bci_data(h5_data, selected_channels, electrode_labels, states, patie
 
     return save_dict
 
-def plot_bci_states(partitioned_data, states, patient_id):
-    state_names = states[patient_id]['center_out']
+def plot_bci_states(partitioned_data, states, patient_id, current_experiment):
+    state_names = states[patient_id][current_experiment]
 
     save_dict = {}
     for partition_key, partition_load_func in partitioned_data.items():

src/decoding_pipeline/pipelines/data_processing/nodes.py

@@ -49,26 +49,26 @@ def create_pipeline(**kwargs) -> Pipeline:
         ),
         node(
             func=extract_bci_data,
-            inputs=["center_out_hdf5", "selected_channels", "electrode_labels", "params:bci_states", "params:patient_id", "params:gain"],
+            inputs=["center_out_hdf5", "selected_channels", "electrode_labels", "params:bci_states", "params:patient_id", "params:gain", "params:current_experiment"],
             outputs="center_out_extracted_pkl",
             name="extract_bci_data_node"
         ),
         node(
             func=extract_bci_data,
-            inputs=["calibration_hdf5", "selected_channels", "electrode_labels", "params:bci_states", "params:patient_id", "params:gain"],
+            inputs=["calibration_hdf5", "selected_channels", "electrode_labels", "params:bci_states", "params:patient_id", "params:gain", "params:current_calibration"],
             outputs="calibration_extracted_pkl",
             name="extract_calibration_data_node"
         ),
     ],
     namespace="data_extraction",
     inputs=set(["calibration_hdf5", "center_out_hdf5", "selected_channels"]),
     outputs=set(["center_out_extracted_pkl", "calibration_extracted_pkl"]),
-    parameters={"params:patient_id": "params:patient_id", "params:gain": "params:gain", "params:bci_states": "params:bci_states"})
+    parameters={"params:patient_id": "params:patient_id", "params:gain": "params:gain", "params:bci_states": "params:bci_states", "params:current_experiment": "params:current_experiment", "params:current_calibration": "params:current_calibration"})
 
     dataset_metrics_pipeline = pipeline([
         node(
             func=plot_bci_states,
-            inputs=["center_out_extracted_pkl", "params:bci_states", "params:patient_id"],
+            inputs=["center_out_extracted_pkl", "params:bci_states", "params:patient_id", "params:current_experiment"],
             outputs="state_plots",
             name="plot_bci_states_node"
 
@@ -77,7 +77,7 @@ def create_pipeline(**kwargs) -> Pipeline:
     namespace="dataset_metrics",
     inputs=set(["center_out_extracted_pkl"]),
     outputs="state_plots",
-    parameters={"params:patient_id": "params:patient_id", "params:bci_states": "params:bci_states"})
+    parameters={"params:patient_id": "params:patient_id", "params:bci_states": "params:bci_states", "params:current_experiment": "params:current_experiment"})
 
     # return channel_labelling_pipeline + data_extraction_pipeline
 
@@ -86,5 +86,5 @@ def create_pipeline(**kwargs) -> Pipeline:
         namespace="data_preprocessing",
         inputs=set(["calibration_hdf5", "center_out_hdf5"]),
         outputs={"prefixed_channels": "prefixed_channels", "center_out_extracted_pkl": "center_out_extracted_pkl", "calibration_extracted_pkl": "calibration_extracted_pkl", "state_plots": "state_plots", "selected_channels": "selected_channels"},
-        parameters={"params:patient_id": "params:patient_id", "params:gain": "params:gain", "params:bci_states": "params:bci_states"}
+        parameters={"params:patient_id": "params:patient_id", "params:gain": "params:gain", "params:bci_states": "params:bci_states", "params:current_experiment": "params:current_experiment", "params:current_calibration": "params:current_calibration"}
     )

src/decoding_pipeline/pipelines/data_processing/pipeline.py

@@ -3,6 +3,7 @@
 generated using Kedro 0.18.3
 """
 import copy
+from turtle import update
 import numpy as np
 import scipy.signal as signal
 
@@ -229,12 +230,14 @@ def downsample_data_to_spectrogram(partitioned_sxx_data, partitioned_signal_data
         save_dict[partition_sxx_key] = create_closure_func(_downsample_data_to_spectrogram, sxx_data_func, signal_data_func)
     return save_dict
 
-def plot_downsampled_signals(partitioned_sxx_data, partitioned_signal_data):
+def plot_downsampled_signals(partitioned_sxx_data, partitioned_signal_data, partitioned_sxx_std_data):
     save_dict = {}
     for partition_sxx_key, sxx_data_func in partitioned_sxx_data.items():
         signal_data_func = partitioned_signal_data[partition_sxx_key]
+        sxx_std_func = partitioned_sxx_std_data[partition_sxx_key]
 
         sxx_data_dict = sxx_data_func()
+        sxx_std_data_dict = sxx_std_func()
         signal_data_dict = signal_data_func()
 
         stimuli = signal_data_dict['stimuli']
@@ -245,7 +248,9 @@ def plot_downsampled_signals(partitioned_sxx_data, partitioned_signal_data):
         t = sxx_data_dict['t']
         sxx = sxx_data_dict['sxx']
 
-        fig, (ax, ax1, ax2) = plt.subplots(3, figsize=(20,10))
+        sxx_std = sxx_std_data_dict['sxx']
+
+        fig, (ax, ax1, ax2, ax3) = plt.subplots(4, figsize=(20,10))
 
         ax.plot(t_seconds, stimuli[:, 0], color='k', linewidth=1)
         ax.margins(x=0)
@@ -268,11 +273,117 @@ def plot_downsampled_signals(partitioned_sxx_data, partitioned_signal_data):
 
         ax2.set_ylabel('Frequency (Hz)')
         ax2.set_ylim([0, 140])
-        ax2.set_xlabel('Time (s)')
 
-        ax.set_title('Downsampled States, Signals and Spectrogram')
+        ax3.pcolormesh(
+            t,
+            f,
+            sxx_std[:,:,0],
+            # norm=mpl.colors.PowerNorm(gamma=1.0 / 5),
+            cmap="seismic",
+            vmin=-3, 
+            vmax=3
+        #     cmap="YlGnBu"
+        )
+
+        ax3.set_ylabel('Frequency (Hz)')
+        ax3.set_ylim([0, 140])
+
+
+        ax3.set_xlabel('Time (s)')
+
+        ax.set_title('Downsampled States, Signals, Raw Spectrogram and Standardized Spectrogram')
 
         save_dict[f"{partition_sxx_key}.png"] = fig
 
         plt.close()
-    return save_dict
+    return save_dict
+
+def extract_calibration_statistics(partitioned_calibration_sxx, partitioned_calibration_data, selected_sessions, patient_id):
+    calibration_dict = selected_sessions[patient_id]['calibration']
+
+    # Find all dates that have no session data. By default, all sessions will be included for calibration
+    updated_dict = {}
+    for partition_key in list(partitioned_calibration_sxx.keys()):
+        # TODO: Switch the array based splitting to regex based splitting
+        date = partition_key.split('_')[-2]
+        session = partition_key.split('_')[-1]
+        if date in list(calibration_dict.keys()):
+            continue
+        else:
+            sessions_list = updated_dict.setdefault(date, [])
+            sessions_list.append(session)
+            updated_dict[date] = sessions_list
+    
+    calibration_dict.update(updated_dict)
+
+    stimuli = None
+    save_dict = {}
+    for date_key, sessions_list in calibration_dict.items():
+        intermed_list = []
+        for partition_sxx_key, sxx_data_func in partitioned_calibration_sxx.items():
+            continue_loop = True
+            if date_key in partition_sxx_key and partition_sxx_key.split('_')[-1] in sessions_list:
+                all_stimuli=partitioned_calibration_data[partition_sxx_key]()['stimuli']
+                stimuli=all_stimuli[:, 0]
+
+                # Check to make sure calibration stimuli is not all zeros
+                assert 1 in stimuli, "Stimuli contains all zeros, think about possibly excluding this calibration session"
+
+                continue_loop = False
+            
+            if continue_loop:
+                continue
+
+            data_dict = sxx_data_func()
+
+            sxx = data_dict['sxx']
+
+            sxx = sxx[:, stimuli == 1, :]
+
+            mean = np.mean(sxx, axis=1)[:,np.newaxis,:]
+            std = np.std(sxx, axis=1)[:,np.newaxis,:]
+            sxx_len = sxx.shape[1]
+
+            intermed_list.append({
+                'mean': mean,
+                'std': std,
+                'count': sxx_len
+            })
+        
+        total_len = np.sum([x['count'] for x in intermed_list])
+        fractions_len = [x['count']/total_len for x in intermed_list]
+
+        global_mean = sum([x['mean']*frac for x,frac in zip(intermed_list, fractions_len)])
+        global_std = sum([x['std']*frac for x,frac in zip(intermed_list, fractions_len)])
+
+        save_dict[f'Calibration_statistics_{date_key}'] = {
+            'mean': global_mean,
+            'std': global_std 
+        }
+    
+    return save_dict
+
+def _standardize_spectrogram(sxx_func, stats_func):
+    sxx_dict = sxx_func()
+    stats_dict = stats_func()
+
+    sxx = sxx_dict['sxx']
+    mean_sxx = stats_dict['mean']
+    std_sxx = stats_dict['std']
+
+    sxx_dict['sxx'] = (sxx - mean_sxx)/std_sxx
+
+    return sxx_dict
+
+def standardize_spectrograms(partitioned_sxx, partitioned_statistics):
+    
+    save_dict = {}
+    for partition_sxx_key, sxx_data_func in partitioned_sxx.items():
+        date = partition_sxx_key.split('_')[-2]
+        statistics_key = list(filter(lambda x: x.split('_')[-1] == date, partitioned_statistics.keys()))[0]
+
+        stats_func = partitioned_statistics[statistics_key]
+
+        save_dict[partition_sxx_key] = create_closure_func(_standardize_spectrogram, sxx_data_func, stats_func)
+
+    return save_dict