Merge pull request #31 from neuropsychology/ecg_eventrelated-enhancement

Ecg eventrelated enhancement - 0.1.99

neurokit/__init__.py

@@ -1,7 +1,7 @@
 """
 neurokit module.
 """
-__version__ = "0.1.93"
+__version__ = "0.1.99"
 
 from .miscellaneous import *
 from .statistics import *

neurokit/bio/bio_ecg.py

@@ -821,38 +821,45 @@ def ecg_EventRelated(epoch, event_length=1, window_post=0):
     References
     -----------
     """
+    def compute_features(variable, prefix, response):
+        """
+        Internal function to compute features and avoid spaguetti code.
+        """
+        response[prefix + "_Baseline"] = epoch[variable][0]
+        response[prefix + "_Min"] = epoch[variable][0:window_end].min()
+        response[prefix + "_MinDiff"] = response[prefix + "_Min"] - response[prefix + "_Baseline"]
+        response[prefix + "_MinTime"] = epoch[variable][0:window_end].idxmin()
+        response[prefix + "_Max"] = epoch[variable][0:window_end].max()
+        response[prefix + "_MaxDiff"] = response[prefix + "_Max"] - response[prefix + "_Baseline"]
+        response[prefix + "_MaxTime"] = epoch[variable][0:window_end].idxmax()
+        response[prefix + "_Mean"] = epoch[variable][0:window_end].mean()
+        response[prefix + "_MeanDiff"] = response[prefix + "_Mean"] - response[prefix + "_Baseline"]
+
+        return(response)
+
     # Initialization
     ECG_Response = {}
     window_end = event_length + window_post
 
     # Heart Rate
     # =============
     if "Heart_Rate" in epoch.columns:
-        ECG_Response["Heart_Rate_Baseline"] = epoch["Heart_Rate"].ix[0]
-        ECG_Response["Heart_Rate_Min"] = epoch["Heart_Rate"].ix[0:window_end].min()
-        ECG_Response["Heart_Rate_MinDiff"] = ECG_Response["Heart_Rate_Min"] - ECG_Response["Heart_Rate_Baseline"]
-        ECG_Response["Heart_Rate_MinTime"] = epoch["Heart_Rate"].ix[0:window_end].idxmin()
-        ECG_Response["Heart_Rate_Max"] = epoch["Heart_Rate"].ix[0:window_end].max()
-        ECG_Response["Heart_Rate_MaxDiff"] = ECG_Response["Heart_Rate_Max"] - ECG_Response["Heart_Rate_Baseline"]
-        ECG_Response["Heart_Rate_MaxTime"] = epoch["Heart_Rate"].ix[0:window_end].idxmax()
-        ECG_Response["Heart_Rate_Mean"] = epoch["Heart_Rate"].ix[0:window_end].mean()
-        ECG_Response["Heart_Rate_MeanDiff"] = ECG_Response["Heart_Rate_Mean"] - ECG_Response["Heart_Rate_Baseline"]
-
-
+        ECG_Response = compute_features("Heart_Rate", "ECG_Heart_Rate", ECG_Response)
+#
     # Cardiac Phase
     # =============
     if "ECG_Systole" in epoch.columns:
-        ECG_Response["ECG_Phase_Systole"] = epoch["ECG_Systole"].ix[0]
+        ECG_Response["ECG_Phase_Systole"] = epoch["ECG_Systole"][0]
 
         # Identify beginning and end
         systole_beg = np.nan
         systole_end = np.nan
         for i in epoch[0:window_end].index:
-            if epoch["ECG_Systole"].ix[i] != ECG_Response["ECG_Phase_Systole"]:
+            if epoch["ECG_Systole"][i] != ECG_Response["ECG_Phase_Systole"]:
                 systole_end = i
                 break
         for i in epoch[:0].index[::-1]:
-            if epoch["ECG_Systole"].ix[i] != ECG_Response["ECG_Phase_Systole"]:
+            if epoch["ECG_Systole"][i] != ECG_Response["ECG_Phase_Systole"]:
                 systole_beg = i
                 break
 
@@ -863,39 +870,63 @@ def ecg_EventRelated(epoch, event_length=1, window_post=0):
     # RR Interval
     # ==================
     if "ECG_RR_Interval" in epoch.columns:
-        ECG_Response["ECG_RRi_Baseline"] = epoch["ECG_RR_Interval"].ix[0]
-        ECG_Response["ECG_RRi_Min"] = epoch["ECG_RR_Interval"].ix[0:window_end].min()
-        ECG_Response["ECG_RRi_MinDiff"] = ECG_Response["ECG_RRi_Min"] - ECG_Response["ECG_RRi_Baseline"]
-        ECG_Response["ECG_RRi_MinTime"] = epoch["ECG_RR_Interval"].ix[0:window_end].idxmin()
-        ECG_Response["ECG_RRi_Max"] = epoch["ECG_RR_Interval"].ix[0:window_end].max()
-        ECG_Response["ECG_RRi_MaxDiff"] = ECG_Response["ECG_RRi_Max"] - ECG_Response["ECG_RRi_Baseline"]
-        ECG_Response["ECG_RRi_MaxTime"] = epoch["ECG_RR_Interval"].ix[0:window_end].idxmax()
-        ECG_Response["ECG_RRi_Mean"] = epoch["ECG_RR_Interval"].ix[0:window_end].mean()
-        ECG_Response["ECG_RRi_MeanDiff"] = ECG_Response["ECG_RRi_Mean"] - ECG_Response["ECG_RRi_Baseline"]
+        ECG_Response = compute_features("ECG_RR_Interval", "ECG_RRi", ECG_Response)
+
 
     # RSA
     # ==========
     if "RSA" in epoch.columns:
-        ECG_Response["ECG_RSA_Baseline"] = epoch["RSA"].ix[0]
-        ECG_Response["ECG_RSA_Min"] = epoch["RSA"].ix[0:window_end].min()
-        ECG_Response["ECG_RSA_MinDiff"] = ECG_Response["ECG_RSA_Min"] - ECG_Response["ECG_RSA_Baseline"]
-        ECG_Response["ECG_RSA_MinTime"] = epoch["RSA"].ix[0:window_end].idxmin()
-        ECG_Response["ECG_RSA_Max"] = epoch["RSA"].ix[0:window_end].max()
-        ECG_Response["ECG_RSA_MaxDiff"] = ECG_Response["ECG_RSA_Max"] - ECG_Response["ECG_RSA_Baseline"]
-        ECG_Response["ECG_RSA_MaxTime"] = epoch["RSA"].ix[0:window_end].idxmax()
-        ECG_Response["ECG_RSA_Mean"] = epoch["RSA"].ix[0:window_end].mean()
-        ECG_Response["ECG_RSA_MeanDiff"] = ECG_Response["ECG_RSA_Mean"] - ECG_Response["ECG_RSA_Baseline"]
+        ECG_Response = compute_features("RSA", "ECG_RSA", ECG_Response)
 
     # HRV
     # ====
     if "ECG_R_Peaks" in epoch.columns:
-        rpeaks = epoch[epoch["ECG_R_Peaks"]==1].ix[0:event_length].index*1000
+        rpeaks = epoch[epoch["ECG_R_Peaks"]==1][0:event_length].index*1000
         hrv = ecg_hrv(rpeaks, sampling_rate=1000, hrv_features=["time"])
 
+        # HRV time domain feature computation
         for key in hrv:
             if isinstance(hrv[key], float):  # Avoid storing series or dataframes
                 ECG_Response["ECG_HRV_" + key] = hrv[key]
 
+        # Computation for baseline
+        if epoch.index[0] > -4:  # Sanity check
+            print("NeuroKit Warning: ecg_EventRelated(): your epoch starts less than 4 seconds before stimulus onset. That's too short to compute HRV baseline features.")
+        else:
+            rpeaks = epoch[epoch["ECG_R_Peaks"]==1][:0].index*1000
+            hrv = ecg_hrv(rpeaks, sampling_rate=1000, hrv_features=["time"])
+
+            for key in hrv:
+                if isinstance(hrv[key], float):  # Avoid storing series or dataframes
+                    ECG_Response["ECG_HRV_" + key + "_Baseline"] = hrv[key]
+
+            # Compute differences between features and baseline
+            keys = [key for key in ECG_Response.keys() if '_Baseline' in key]  # Find keys
+            keys = [key for key in keys if 'ECG_HRV_' in key]
+            keys = [s.replace('_Baseline', '') for s in keys]  # Remove baseline part
+            for key in keys:
+                try:
+                    ECG_Response[key + "_Diff"] = ECG_Response[key] - ECG_Response[key + "_Baseline"]
+                except KeyError:
+                    ECG_Response[key + "_Diff"] = np.nan
+
+
+
+    if "ECG_HRV_VHF" in epoch.columns:
+        ECG_Response = compute_features("ECG_HRV_VHF", "ECG_HRV_VHF", ECG_Response)
+
+    if "ECG_HRV_HF" in epoch.columns:
+        ECG_Response = compute_features("ECG_HRV_HF", "ECG_HRV_HF", ECG_Response)
+
+    if "ECG_HRV_LF" in epoch.columns:
+        ECG_Response = compute_features("ECG_HRV_LF", "ECG_HRV_LF", ECG_Response)
+
+    if "ECG_HRV_VLF" in epoch.columns:
+        ECG_Response = compute_features("ECG_HRV_VLF", "ECG_HRV_VLF", ECG_Response)
+
+
+
+
     return(ECG_Response)
 
 

neurokit/eeg/eeg_complexity.py

@@ -8,6 +8,23 @@
 import re
 
 
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+#def eeg_complexity(epochs):
+#    """
+#    """
+#    df = epochs.to_data_frame(index=["epoch", "time", "condition"])
+#
+#    for index, epoch in enumerate(epochs):
+#        df = epoch
+
+
 
 # ==============================================================================
 # ==============================================================================

neurokit/eeg/eeg_data.py

@@ -314,18 +314,24 @@ def eeg_create_events(onsets, conditions=None):
     """
     event_id = {}
 
+    # Sanity check
+    if len(conditions) != len(onsets):
+        print("NeuroKit Warning: eeg_create_events(): conditions parameter of different length than onsets. Aborting.")
+        return()
+
     if conditions is None:
         conditions = ["Event"] * len(onsets)
 
 
+
     event_names = list(set(conditions))
 #    event_index = [1, 2, 3, 4, 5, 32, 64, 128]
     event_index = list(range(len(event_names)))
     for i in enumerate(event_names):
         conditions = [event_index[i[0]] if x==i[1] else x for x in conditions]
         event_id[i[1]] = event_index[i[0]]
 
-    events = np.array([onsets, [0]*len(onsets), conditions]).T
+    events = np.array([onsets, [0]*len(onsets), conditions])
     return(events, event_id)
 
 
@@ -421,8 +427,105 @@ def eeg_add_events(raw, events_channel, conditions=None, treshold="auto", cut="h
 
 
 
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+def eeg_epochs_to_dict(epochs, include="all", exclude=None, hemisphere="both", include_central=True):
+    """
+    Convert mne.Epochs object to Python dict.
+    """
+    data = {}
+    for index, epoch in enumerate(epochs.get_data()):
+        epoch = pd.DataFrame(epoch.T)
+        epoch.columns = epochs.ch_names
+
+        selection = eeg_select_sensor_area(include=include, exclude=exclude, hemisphere=hemisphere, include_central=include_central)
+
+        data[index] = epoch[selection]
+    return()
+
+
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+# ==============================================================================
+def eeg_select_sensor_area(include="all", exclude=None, hemisphere="both", include_central=True):
+    """
+    Returns list of electrodes names (according to a 10-20 EEG montage). This function is probably not very flexibile. Looking for help to improve it.
+
+    Parameters
+    ----------
+    include : str
+        Sensor area to include.
+    exclude : str or None
+        Sensor area to exclude.
+    hemisphere : str
+        Select both hemispheres? "both", "left" or "right".
+    include_central : bool
+        if `hemisphere != "both"`, select the central line?
+
+    Returns
+    ----------
+    sensors : list
+        List of sensors corresponding to the selected area.
+
+    Example
+    ----------
+    >>> import neurokit as nk
+    >>> nk.eeg_select_sensor_area(include="F", exclude="C")
+
+
+    Notes
+    ----------
+    *Authors*
+
+    - Dominique Makowski (https://github.com/DominiqueMakowski)
+
+    References
+    ------------
+    - None
+    """
+    sensors = ['AF3', 'AF4', 'AF7', 'AF8', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'CP1', 'CP2', 'CP3', 'CP4', 'CP5', 'CP6', 'CPz', 'Cz', 'F1', 'F2', 'F3', 'F4', 'F5', 'F6', 'F7', 'F8', 'FC1', 'FC2', 'FC3', 'FC4', 'FC5', 'FC6', 'Fp1', 'Fp2', 'FT10', 'FT7', 'FT8', 'FT9', 'O1', 'O2', 'Oz', 'P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'PO3', 'PO4', 'PO7', 'PO8', 'POz', 'Pz', 'FCz', 'T7', 'T8', 'TP10', 'TP7', 'TP8', 'TP9', 'AFz']
+
+    if include != "all":
+        sensors = [s for s in sensors if include in s]
+
+    if exclude != None:
+        if isinstance(exclude, str):
+            exclude = [exclude]
+        for to_exclude in exclude:
+            sensors = [s for s in sensors if to_exclude not in s]
+
+    if hemisphere != "both":
+        if include_central == False:
+            if hemisphere == "left":
+                sensors = [s for s in sensors if "1" in s or "3" in s or "5" in s or "7" in s or "9" in s]
+            if hemisphere == "right":
+                sensors = [s for s in sensors if "2" in s or "4" in s or "6" in s or "8" in s or "10" in s]
+        else:
+            if hemisphere == "left":
+                sensors = [s for s in sensors if "1" in s or "3" in s or "5" in s or "7" in s or "9" in s or "z" in s]
+            if hemisphere == "right":
+                sensors = [s for s in sensors if "2" in s or "4" in s or "6" in s or "8" in s or "10" in s or "z" in s]
+
+
+    return(sensors)
+
 
 
+#==============================================================================
+#==============================================================================
+#==============================================================================
+#==============================================================================
 # ==============================================================================
 # ==============================================================================
 # ==============================================================================

neurokit/eeg/eeg_erp.py

@@ -1,7 +1,7 @@
 """
 ERP analysis EEG submodule.
 """
-from .eeg_preprocessing import eeg_select_sensor_area
+from .eeg_data import eeg_select_sensor_area
 import numpy as np
 import mne
 import matplotlib

neurokit/eeg/eeg_preprocessing.py

@@ -315,72 +315,3 @@ def eeg_epoching(raw, events, event_id, tmin=-0.2, tmax=1, eog_reject=600e-6, pr
         epochs.pick_types(meg=False, eeg=True)
     return(epochs)
 
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-# ==============================================================================
-def eeg_select_sensor_area(include="all", exclude=None, hemisphere="both", include_central=True):
-    """
-    Returns list of electrodes names (according to a 10-20 EEG montage). This function is probably not very flexibile. Looking for help to improve it.
-
-    Parameters
-    ----------
-    include : str
-        Sensor area to include.
-    exclude : str or None
-        Sensor area to exclude.
-    hemisphere : str
-        Select both hemispheres? "both", "left" or "right".
-    include_central : bool
-        if `hemisphere != "both"`, select the central line?
-
-    Returns
-    ----------
-    sensors : list
-        List of sensors corresponding to the selected area.
-
-    Example
-    ----------
-    >>> import neurokit as nk
-    >>> nk.eeg_select_sensor_area(include="F", exclude="C")
-
-
-    Notes
-    ----------
-    *Authors*
-
-    - Dominique Makowski (https://github.com/DominiqueMakowski)
-
-    References
-    ------------
-    - None
-    """
-    sensors = ['AF3', 'AF4', 'AF7', 'AF8', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'CP1', 'CP2', 'CP3', 'CP4', 'CP5', 'CP6', 'CPz', 'Cz', 'F1', 'F2', 'F3', 'F4', 'F5', 'F6', 'F7', 'F8', 'FC1', 'FC2', 'FC3', 'FC4', 'FC5', 'FC6', 'Fp1', 'Fp2', 'FT10', 'FT7', 'FT8', 'FT9', 'O1', 'O2', 'Oz', 'P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'PO3', 'PO4', 'PO7', 'PO8', 'POz', 'Pz', 'FCz', 'T7', 'T8', 'TP10', 'TP7', 'TP8', 'TP9', 'AFz']
-
-    if include != "all":
-        sensors = [s for s in sensors if include in s]
-
-    if exclude != None:
-        if isinstance(exclude, str):
-            exclude = [exclude]
-        for to_exclude in exclude:
-            sensors = [s for s in sensors if to_exclude not in s]
-
-    if hemisphere != "both":
-        if include_central == False:
-            if hemisphere == "left":
-                sensors = [s for s in sensors if "1" in s or "3" in s or "5" in s or "7" in s or "9" in s]
-            if hemisphere == "right":
-                sensors = [s for s in sensors if "2" in s or "4" in s or "6" in s or "8" in s or "10" in s]
-        else:
-            if hemisphere == "left":
-                sensors = [s for s in sensors if "1" in s or "3" in s or "5" in s or "7" in s or "9" in s or "z" in s]
-            if hemisphere == "right":
-                sensors = [s for s in sensors if "2" in s or "4" in s or "6" in s or "8" in s or "10" in s or "z" in s]
-
-
-    return(sensors)

neurokit/eeg/eeg_time_frequency.py

@@ -1,7 +1,7 @@
 """
 Time-frequency submodule.
 """
-from .eeg_preprocessing import eeg_select_sensor_area
+from .eeg_data import eeg_select_sensor_area
 from ..miscellaneous import Time
 
 import numpy as np