Repository: neuromodulation/py_neuromodulation (language: python)

Issue:
Refactor bids stream
Main change: Generalize Bids Stream to Generic Stream class

I think for many problems data in in BIDS format won't be available. This leaves the main requirements of the ABC nm_stream to be:
 - channel names
 - channel types
 - bad channels
 - rereferencing
 - fs 
 - line_noise
and subsequently defintion of a "data" np array for offline processing.

Thus I renamed `nm_BidsStream` to be `nm_GenericStream`.

Also `scripts/pn_nm_batch_trd_data.py` shows an example of utilizing epoch based data into a "stream" of data. 
This has implications of epoch discontinuities, since concatenating might not be optimal in some cases.

The advantages are though that the general pipeline is followed one to one compared to the "online" pipeline, so all pipeline code, analysis, potting, decoding stays the same.



Discussion / hints:

Fix patch:
diff --git a/examples/example_BIDS.py b/examples/example_BIDS.py
index e6c3e4ba..2b42ec54 100644
--- a/examples/example_BIDS.py
+++ b/examples/example_BIDS.py
@@ -1,36 +1,59 @@
 import os
+
+import py_neuromodulation as nm
+import xgboost
+from py_neuromodulation import (
+    nm_analysis,
+    nm_decode,
+    nm_define_nmchannels,
+    nm_IO,
+)
 from sklearn import linear_model, metrics, model_selection
 from skopt import space as skopt_space
-import xgboost
-from py_neuromodulation import nm_BidsStream, nm_analysis, nm_decode
 
 
 def run_example_BIDS():
     """run the example BIDS path in pyneuromodulation/tests/data"""
+    sub = "testsub"
+    ses = "EphysMedOff"
+    task = "buttonpress"
+    run = 0
+    datatype = "ieeg"
 
-    RUN_NAME = "sub-testsub_ses-EphysMedOff_task-buttonpress_run-0"
+    RUN_NAME = f"sub-{sub}_ses-{ses}_task-{task}_run-{run}"
     PATH_RUN = os.path.join(
         os.path.abspath(os.path.join("examples", "data")),
-        "sub-testsub",
-        "ses-EphysMedOff",
-        "ieeg",
+        f"sub-{sub}",
+        f"ses-{ses}",
+        datatype,
         RUN_NAME,
     )
     PATH_BIDS = os.path.abspath(os.path.join("examples", "data"))
     PATH_OUT = os.path.abspath(os.path.join("examples", "data", "derivatives"))
 
-    # read default settings
-    nm_BIDS = nm_BidsStream.BidsStream(
-        PATH_RUN=PATH_RUN,
-        PATH_BIDS=PATH_BIDS,
-        PATH_OUT=PATH_OUT,
-        LIMIT_DATA=False,
-        VERBOSE=True,
+    raw, _, _, _ = nm_IO.read_BIDS_data(
+        PATH_RUN=PATH_RUN, BIDS_PATH=PATH_BIDS, datatype=datatype
+    )
+
+    nm_channels = nm_define_nmchannels.set_channels(
+        ch_names=raw.ch_names,
+        ch_types=raw.get_channel_types(),
+        reference="default",
+        bads=raw.info["bads"],
+        new_names="default",
+        used_types=("ecog", "dbs"),
         target_keywords=("SQUARED_ROTATION",),
-        used_types=("ecog", "seeg"),
     )
 
-    nm_BIDS.run_bids()
+    nm_BIDS = nm.BidsStream(
+        settings=None,
+        nm_channels=nm_channels,
+        path_out=PATH_OUT,
+        path_grids=None,
+        verbose=True,
+    )
+
+    nm_BIDS.run(filepath=PATH_RUN, bids_root=PATH_BIDS, datatype=datatype)
 
     # init analyzer
     feature_reader = nm_analysis.Feature_Reader(
@@ -41,14 +64,19 @@ def run_example_BIDS():
     # feature_reader.plot_cort_projection()
 
     # plot for a single channel
-    ch_used = feature_reader.nm_channels.query('(type=="ecog") and (used == 1)').iloc[
-        0
-    ]["name"]
+    ch_used = feature_reader.nm_channels.query(
+        '(type=="ecog") and (used == 1)'
+    ).iloc[0]["name"]
 
-    feature_used = "stft" if feature_reader.settings["methods"]["stft"] else "fft"
+    feature_used = (
+        "stft" if feature_reader.settings["methods"]["stft"] else "fft"
+    )
 
     feature_reader.plot_target_averaged_channel(
-        ch=ch_used, list_feature_keywords=[feature_used], epoch_len=4, threshold=0.5
+        ch=ch_used,
+        list_feature_keywords=[feature_used],
+        epoch_len=4,
+        threshold=0.5,
     )
 
     # model = linear_model.LogisticRegression(class_weight='balanced')
@@ -56,7 +84,9 @@ def run_example_BIDS():
 
     bay_opt_param_space = [
         skopt_space.Integer(1, 100, name="max_depth"),
-        skopt_space.Real(10**-5, 10**0, "log-uniform", name="learning_rate"),
+        skopt_space.Real(
+            10**-5, 10**0, "log-uniform", name="learning_rate"
+        ),
         skopt_space.Real(10**0, 10**1, "uniform", name="gamma"),
     ]
 
diff --git a/py_neuromodulation/__init__.py b/py_neuromodulation/__init__.py
index 16d2143a..cdb3369c 100644
--- a/py_neuromodulation/__init__.py
+++ b/py_neuromodulation/__init__.py
@@ -1,8 +1,9 @@
 from . import (
     nm_decode,
     nm_analysis,
-    nm_BidsStream,
-    nm_stream,
+    nm_stream_abc,
     nm_cohortwrapper,
     nm_across_patient_decoding,
+    nm_stream_offline,
 )
+from .nm_stream_offline import Stream, BidsStream
diff --git a/py_neuromodulation/nm_BidsStream.py b/py_neuromodulation/nm_BidsStream.py
deleted file mode 100644
index 79b0965e..00000000
--- a/py_neuromodulation/nm_BidsStream.py
+++ /dev/null
@@ -1,247 +0,0 @@
-from typing import Optional, Iterable
-import os
-import pathlib
-import mne
-import numpy as np
-import pandas as pd
-import timeit
-
-from py_neuromodulation import nm_IO, nm_generator, nm_stream
-
-
-class BidsStream(nm_stream.PNStream):
-
-    PATH_ANNOTATIONS: str
-    PATH_GRIDS: str
-    PATH_BIDS: str
-    PATH_RUN: str
-    LIMIT_DATA: bool
-    LIMIT_LOW: int
-    LIMIT_HIGH: int
-    used_types: bool
-    raw_arr: mne.io.RawArray
-    raw_arr_data: np.array
-    annot: mne.Annotations
-    VERBOSE: bool = False
-
-    def __init__(
-        self,
-        PATH_RUN,
-        PATH_SETTINGS: str = os.path.join(
-            pathlib.Path(__file__).parent.resolve(), "nm_settings.json"
-        ),
-        PATH_NM_CHANNELS: str = str(),
-        PATH_OUT: str = os.getcwd(),
-        PATH_GRIDS: str = pathlib.Path(__file__).parent.resolve(),
-        VERBOSE: bool = True,
-        PATH_ANNOTATIONS: str = str(),
-        PATH_BIDS: str = str(),
-        LIMIT_DATA: bool = False,
-        LIMIT_LOW: int = 0,
-        LIMIT_HIGH: int = 10000,
-        used_types: Optional[Iterable[str]] = ("ecog", "dbs", "seeg"),
-        target_keywords: Optional[Iterable[str]] = ("mov", "squared", "label"),
-    ) -> None:
-
-        super().__init__(
-            PATH_SETTINGS=PATH_SETTINGS,
-            PATH_NM_CHANNELS=PATH_NM_CHANNELS,
-            PATH_OUT=PATH_OUT,
-            PATH_GRIDS=PATH_GRIDS,
-            VERBOSE=VERBOSE,
-        )
-
-        self.PATH_RUN = PATH_RUN
-        self.PATH_ANNOTATIONS = PATH_ANNOTATIONS
-        self.PATH_BIDS = PATH_BIDS
-        self.used_types = used_types
-        self.LIMIT_DATA = LIMIT_DATA
-        self.LIMIT_LOW = LIMIT_LOW
-        self.LIMIT_HIGH = LIMIT_HIGH
-
-        self.raw_arr, self.raw_arr_data, fs, line_noise = nm_IO.read_BIDS_data(
-            self.PATH_RUN, self.PATH_BIDS
-        )
-
-        self.set_fs(fs)
-
-        self.set_linenoise(line_noise)
-
-        self.nm_channels = self._get_nm_channels(
-            PATH_NM_CHANNELS,
-            ch_names=self.raw_arr.ch_names,
-            ch_types=self.raw_arr.get_channel_types(),
-            bads=self.raw_arr.info["bads"],
-            used_types=used_types,
-            target_keywords=target_keywords,
-            reference="default",
-        )
-
-        if self.PATH_ANNOTATIONS:
-            self.annot, self.annot_data, self.raw_arr = nm_IO.get_annotations(
-                self.PATH_ANNOTATIONS, self.PATH_RUN, self.raw_arr
-            )
-        if self.LIMIT_DATA:
-            self.raw_arr_data = self.raw_arr_data[:, LIMIT_LOW:LIMIT_HIGH]
-
-        self.coord_list, self.coord_names = self._get_bids_coord_list(self.raw_arr)
-
-        if self.coord_list is None and True in [
-            self.settings["project_cortex"],
-            self.settings["project_subcortex"],
-        ]:
-            raise ValueError("no coordinates could be loaded from BIDS Dataset")
-        elif self.coord_list is not None:
-            self.coords = self._add_coordinates(self.coord_names, self.coord_list)
-            self.sess_right = self._get_sess_lat(self.coords)
-
-        self.gen = nm_generator.ieeg_raw_generator(
-            self.raw_arr_data, self.settings, self.fs
-        )
-
-    def run_bids(self) -> None:
-        """process BIDS recording, add labels and save features after finish"""
-
-        # init features, projection etc. here
-        # settings, nm_channels might have been changed by the user in e.g. the ipynb
-        self._set_run()
-
-        self.run()
-
-        self._add_labels()
-
-        folder_name = os.path.basename(self.PATH_RUN)
-
-        self.save_after_stream(folder_name)
-
-    def get_data(self) -> np.array:
-        return next(self.gen, None)
-
-    def run(self, predict: bool = False) -> None:
-        """BIDS specific fun function
-        Does not need to run in parallel
-        """
-        # Loop
-        idx = 0
-        while True:
-            data = self.get_data()
-            if data is None:
-                break
-            feature_series = self.run_analysis.process_data(data)
-
-            # Measuring timing
-            # number_repeat = 100
-            # val = timeit.timeit(
-            #    lambda: self.run_analysis.process_data(data),
-            #    number=number_repeat
-            # ) / number_repeat
-
-            feature_series = self._add_timestamp(feature_series, idx)
-
-            # concatenate data to feature_arr
-            if idx == 0:
-                self.feature_arr = pd.DataFrame([feature_series])
-                idx += 1
-            else:
-                self.feature_arr = self.feature_arr.append(
-                    feature_series, ignore_index=True
-                )
-
-            if predict is True:
-                prediction = self.model.predict(feature_series)
-
-    def _add_timestamp(self, feature_series: pd.Series, idx: int = None) -> pd.Series:
-        """time stamp is added in ms
-        Due to normalization run_analysis needs to keep track of the counted samples
-        Those are accessed here for time conversion"""
-
-        if idx == 0:
-            feature_series["time"] = self.run_analysis.offset
-        else:
-            # sampling frequency is taken here from run_analysis, since resampling might change it
-            feature_series["time"] = (
-                self.run_analysis.cnt_samples * 1000 / self.run_analysis.fs
-            )
-
-        if self.VERBOSE:
-            print(
-                str(np.round(feature_series["time"] / 1000, 2))
-                + " seconds of data processed"
-            )
-
-        return feature_series
-
-    def _add_labels(self):
-        """add resampled labels to feature dataframe if there are target channels"""
-        if self.nm_channels.target.sum() > 0:
-            self.feature_arr = nm_IO.add_labels(
-                self.feature_arr,
-                self.settings,
-                self.nm_channels,
-                self.raw_arr_data,
-                self.fs,
-            )
-        else:
-            pass
-
-    def _get_bids_coord_list(self, raw_arr: mne.io.RawArray):
-        if raw_arr.get_montage() is not None:
-            coord_list = np.array(
-                list(dict(raw_arr.get_montage().get_positions()["ch_pos"]).values())
-            ).tolist()
-            coord_names = np.array(
-                list(dict(raw_arr.get_montage().get_positions()["ch_pos"]).keys())
-            ).tolist()
-        else:
-            coord_list = None
-            coord_names = None
-
-        return coord_list, coord_names
-
-    def _add_coordinates(self, coord_names: list, coord_list: np.array):
-        """set coordinate information to settings from RawArray
-        The set coordinate positions are set as lists,
-        since np.arrays cannot be saved in json
-        Parameters
-        ----------
-        raw_arr : mne.io.RawArray
-        PATH_GRIDS : string, optional
-            absolute path to grid_cortex.tsv and grid_subcortex.tsv, by default: None
-        """
-
-        coords = {}
-
-        def left_coord(val, coord_region):
-            if coord_region.split("_")[1] == "left":
-                return val < 0
-            else:
-                return val > 0
-
-        for coord_region in [
-            coord_loc + "_" + lat
-            for coord_loc in ["cortex", "subcortex"]
-            for lat in ["left", "right"]
-        ]:
-
-            coords[coord_region] = {}
-
-            ch_type = "ECOG" if "cortex" == coord_region.split("_")[0] else "LFP"
-
-            coords[coord_region]["ch_names"] = [
-                coord_names[ch_idx]
-                for ch_idx, ch in enumerate(coord_list)
-                if left_coord(coord_list[ch_idx][0], coord_region)
-                and (ch_type in coord_names[ch_idx])
-            ]
-
-            # multiply by 1000 to get m instead of mm
-            coords[coord_region]["positions"] = 1000 * np.array(
-                [
-                    ch
-                    for ch_idx, ch in enumerate(coord_list)
-                    if left_coord(coord_list[ch_idx][0], coord_region)
-                    and (ch_type in coord_names[ch_idx])
-                ]
-            )
-
-        return coords
diff --git a/py_neuromodulation/nm_EpochStream.py b/py_neuromodulation/nm_EpochStream.py
index f0140874..6a01d81f 100644
--- a/py_neuromodulation/nm_EpochStream.py
+++ b/py_neuromodulation/nm_EpochStream.py
@@ -58,7 +58,9 @@ def run(
             )  # None will be returned if generator ran through
             if data is None:
                 break
-            gen = nm_generator.ieeg_raw_generator(data, self.settings, self.fs)
+            gen = nm_generator.ieeg_raw_generator(
+                data, self.settings, self.sfreq
+            )
 
             def get_data_within_epoch() -> np.array:
                 return next(gen, None)
@@ -69,7 +71,9 @@ def get_data_within_epoch() -> np.array:
                 if data_within_epoch is None:
                     break
 
-                feature_series = self.run_analysis.process_data(data_within_epoch)
+                feature_series = self.run_analysis.process_data(
+                    data_within_epoch
+                )
                 if idx_within_epoch == 0:
                     self.feature_arr = pd.DataFrame([feature_series])
                     idx_within_epoch += 1
@@ -79,7 +83,9 @@ def get_data_within_epoch() -> np.array:
                     )
             self.feature_arr_list.append(self.feature_arr)
 
-    def _add_timestamp(self, feature_series: pd.Series, idx: int = None) -> pd.Series:
+    def _add_timestamp(
+        self, feature_series: pd.Series, idx: int = None
+    ) -> pd.Series:
         # in case of epochs no timestamp is necessary
         return feature_series
 
diff --git a/py_neuromodulation/nm_IO.py b/py_neuromodulation/nm_IO.py
index a5e0bc14..501f8735 100644
--- a/py_neuromodulation/nm_IO.py
+++ b/py_neuromodulation/nm_IO.py
@@ -1,19 +1,21 @@
-import mne_bids
-import mne
-import numpy as np
-import os
 import json
-from scipy import io
-import pandas as pd
+import os
 from pathlib import Path
 
+import mne
+import mne_bids
+import numpy as np
+import pandas as pd
+from scipy import io
 
-def read_settings(PATH_SETTINGS: str) -> None:
-    with open(PATH_SETTINGS, encoding="utf-8") as json_file:
-        return json.load(json_file)
+_PathLike = str | os.PathLike
 
 
-def read_BIDS_data(PATH_RUN, BIDS_PATH, datatype: str = "ieeg"):
+def read_BIDS_data(
+    PATH_RUN: _PathLike | mne_bids.BIDSPath,
+    BIDS_PATH: _PathLike | None = None,
+    datatype: str = "ieeg",
+) -> tuple[mne.io.Raw, np.ndarray, int | float, int]:
     """Given a run path and bids data path, read the respective data
 
     Parameters
@@ -29,17 +31,19 @@ def read_BIDS_data(PATH_RUN, BIDS_PATH, datatype: str = "ieeg"):
     fs : int
     line_noise : int
     """
-    entities = mne_bids.get_entities_from_fname(PATH_RUN)
-
-    bids_path = mne_bids.BIDSPath(
-        subject=entities["subject"],
-        session=entities["session"],
-        task=entities["task"],
-        run=entities["run"],
-        acquisition=entities["acquisition"],
-        datatype=datatype,
-        root=BIDS_PATH,
-    )
+    if isinstance(PATH_RUN, mne_bids.BIDSPath):
+        bids_path = PATH_RUN
+    else:
+        entities = mne_bids.get_entities_from_fname(PATH_RUN)
+        bids_path = mne_bids.BIDSPath(
+            subject=entities["subject"],
+            session=entities["session"],
+            task=entities["task"],
+            run=entities["run"],
+            acquisition=entities["acquisition"],
+            datatype=datatype,
+            root=BIDS_PATH,
+        )
 
     raw_arr = mne_bids.read_raw_bids(bids_path)
 
@@ -61,16 +65,21 @@ def read_grid(PATH_GRIDS: str, grid_str: str):
         )
     else:
         grid = pd.read_csv(
-            os.path.join(PATH_GRIDS, "grid_" + grid_str.name.lower() + ".tsv"), sep="\t"
+            os.path.join(PATH_GRIDS, "grid_" + grid_str.name.lower() + ".tsv"),
+            sep="\t",
         )
     return grid
 
 
-def get_annotations(PATH_ANNOTATIONS: str, PATH_RUN: str, raw_arr: mne.io.RawArray):
+def get_annotations(
+    PATH_ANNOTATIONS: str, PATH_RUN: str, raw_arr: mne.io.RawArray
+):
 
     try:
         annot = mne.read_annotations(
-            os.path.join(PATH_ANNOTATIONS, os.path.basename(PATH_RUN)[:-5] + ".txt")
+            os.path.join(
+                PATH_ANNOTATIONS, os.path.basename(PATH_RUN)[:-5] + ".txt"
+            )
         )
         raw_arr.set_annotations(annot)
 
@@ -81,7 +90,9 @@ def get_annotations(PATH_ANNOTATIONS: str, PATH_RUN: str, raw_arr: mne.io.RawArr
         print(
             "expected location: "
             + str(
-                os.path.join(PATH_ANNOTATIONS, os.path.basename(PATH_RUN)[:-5] + ".txt")
+                os.path.join(
+                    PATH_ANNOTATIONS, os.path.basename(PATH_RUN)[:-5] + ".txt"
+                )
             )
         )
     return annot, annot_data, raw_arr
@@ -169,7 +180,9 @@ def add_labels(
             for idx, label_ch in enumerate(nm_channels.name[ind_label]):
                 df_[label_ch] = label_downsampled[idx, :]
         else:
-            print("label dimensions don't match, saving downsampled label extra")
+            print(
+                "label dimensions don't match, saving downsampled label extra"
+            )
     else:
         print("no target specified")
 
@@ -216,23 +229,33 @@ def save_features_and_settings(
 
 def save_settings(settings: dict, PATH_OUT: str, folder_name: str = None):
     if folder_name is not None:
-        PATH_OUT = os.path.join(PATH_OUT, folder_name, folder_name + "_SETTINGS.json")
+        PATH_OUT = os.path.join(
+            PATH_OUT, folder_name, folder_name + "_SETTINGS.json"
+        )
 
     with open(PATH_OUT, "w") as f:
         json.dump(settings, f, indent=4)
     print("settings.json saved to " + str(PATH_OUT))
 
 
-def save_nmchannels(nmchannels: pd.DataFrame, PATH_OUT: str, folder_name: str = None):
+def save_nmchannels(
+    nmchannels: pd.DataFrame, PATH_OUT: str, folder_name: str = None
+):
     if folder_name is not None:
-        PATH_OUT = os.path.join(PATH_OUT, folder_name, folder_name + "_nm_channels.csv")
+        PATH_OUT = os.path.join(
+            PATH_OUT, folder_name, folder_name + "_nm_channels.csv"
+        )
     nmchannels.to_csv(PATH_OUT)
     print("nm_channels.csv saved to " + str(PATH_OUT))
 
 
-def save_features(df_features: pd.DataFrame, PATH_OUT: str, folder_name: str = None):
+def save_features(
+    df_features: pd.DataFrame, PATH_OUT: str, folder_name: str = None
+):
     if folder_name is not None:
-        PATH_OUT = os.path.join(PATH_OUT, folder_name, folder_name + "_FEATURES.csv")
+        PATH_OUT = os.path.join(
+            PATH_OUT, folder_name, folder_name + "_FEATURES.csv"
+        )
     df_features.to_csv(PATH_OUT)
     print("FEATURES.csv saved to " + str(PATH_OUT))
 
@@ -266,7 +289,11 @@ def save_general_dict(
 
     with open(PATH_OUT, "w") as f:
         json.dump(
-            dict_, f, default=default_json_convert, indent=4, separators=(",", ": ")
+            dict_,
+            f,
+            default=default_json_convert,
+            indent=4,
+            separators=(",", ": "),
         )
     print(f"{str_add} saved to " + str(PATH_OUT))
 
@@ -292,3 +319,39 @@ def get_run_list_indir(PATH: str) -> list:
             if "FEATURES" in x:
                 f_files.append(os.path.basename(dirpath))
     return f_files
+
+
+def loadmat(filename):
+    """
+    this function should be called instead of direct spio.loadmat
+    as it cures the problem of not properly recovering python dictionaries
+    from mat files. It calls the function check keys to cure all entries
+    which are still mat-objects
+    """
+    data = io.loadmat(filename, struct_as_record=False, squeeze_me=True)
+    return _check_keys(data)
+
+
+def _check_keys(dict):
+    """
+    checks if entries in dictionary are mat-objects. If yes
+    todict is called to change them to nested dictionaries
+    """
+    for key in dict:
+        if isinstance(dict[key], io.matlab.mio5_params.mat_struct):
+            dict[key] = _todict(dict[key])
+    return dict
+
+
+def _todict(matobj):
+    """
+    A recursive function which constructs from matobjects nested dictionaries
+    """
+    dict = {}
+    for strg in matobj._fieldnames:
+        elem = matobj.__dict__[strg]
+        if isinstance(elem, io.matlab.mio5_params.mat_struct):
+            dict[strg] = _todict(elem)
+        else:
+            dict[strg] = elem
+    return dict
diff --git a/py_neuromodulation/nm_RealTimeClientStreamMultiprocessing.py b/py_neuromodulation/nm_RealTimeClientStreamMultiprocessing.py
index 0e8fbbcf..25f0ec05 100644
--- a/py_neuromodulation/nm_RealTimeClientStreamMultiprocessing.py
+++ b/py_neuromodulation/nm_RealTimeClientStreamMultiprocessing.py
@@ -1,30 +1,22 @@
 import multiprocessing
 import os
-import sys
 import pathlib
-from typing import MutableMapping
-import numpy as np
-import pandas as pd
+import sys
 import time
 import timeit
+
+import numpy as np
+import pandas as pd
 from pynput.keyboard import Key, Listener
 
-from py_neuromodulation import \
-    (nm_projection,
-    nm_rereference,
-    nm_run_analysis,
-    nm_features,
-    nm_resample,
-    nm_stream,
-    nm_test_settings,
-    FieldTrip)
+from py_neuromodulation import FieldTrip, nm_stream_abc
 
 
-class RealTimePyNeuro(nm_stream.PNStream):
+class RealTimePyNeuro(nm_stream_abc.PNStream):
 
-    queue_raw:multiprocessing.Queue = multiprocessing.Queue(1)
-    queue_features:multiprocessing.Queue = multiprocessing.Queue(1)
-    ftc:FieldTrip.Client = None
+    queue_raw: multiprocessing.Queue = multiprocessing.Queue(1)
+    queue_features: multiprocessing.Queue = multiprocessing.Queue(1)
+    ftc: FieldTrip.Client = None
     filename: str
 
     def __init__(
@@ -36,22 +28,22 @@ def __init__(
         PATH_GRIDS: str = pathlib.Path(__file__).parent.resolve(),
         VERBOSE: bool = True,
         fs: int = 128,
-        line_noise: int = 50
+        line_noise: int = 50,
     ) -> None:
 
-        super().__init__(PATH_SETTINGS=PATH_SETTINGS,
-            PATH_NM_CHANNELS=PATH_NM_CHANNELS,
-            PATH_OUT=PATH_OUT,
-            PATH_GRIDS=PATH_GRIDS,
-            VERBOSE=VERBOSE)
+        super().__init__(
+            settings=PATH_SETTINGS,
+            nm_channels=PATH_NM_CHANNELS,
+            path_out=PATH_OUT,
+            path_grids=PATH_GRIDS,
+            verbose=VERBOSE,
+        )
 
         self.set_fs(fs)
         self.set_linenoise(line_noise)
         self.filename = filename
 
-        self.nm_channels = self._get_nm_channels(
-            self.PATH_NM_CHANNELS
-        )
+        self.nm_channels = self._get_nm_channels(self.path_nm_channels)
 
         # leave out coordinate setting for now
 
@@ -68,11 +60,11 @@ def init_fieldtrip() -> FieldTrip.Client:
         """
         ftc = FieldTrip.Client()
         # Python FieldTripBuffer https://www.fieldtriptoolbox.org/development/realtime/buffer_python/
-        ftc.connect('localhost', 1972)   # might throw IOError
+        ftc.connect("localhost", 1972)  # might throw IOError
         H = ftc.getHeader()
 
         if H is None:
-            print('Failed to retrieve header!')
+            print("Failed to retrieve header!")
             sys.exit(1)
 
         print(H)
@@ -81,29 +73,29 @@ def init_fieldtrip() -> FieldTrip.Client:
         return ftc
 
     def run(self) -> None:
-        """Start get_data, calcFeatures and sendFeature processes
-        """
+        """Start get_data, calcFeatures and sendFeature processes"""
         processes = [
             multiprocessing.Process(
-                target=self.get_data,
+                target=self._get_data,
                 args=(
                     self.queue_raw,
                     self.ftc,
-                )
+                ),
             ),
             multiprocessing.Process(
-                target=self.calcFeatures, args=(
+                target=self.calcFeatures,
+                args=(
                     self.queue_raw,
                     self.queue_features,
-                )
+                ),
             ),
             multiprocessing.Process(
                 target=self.sendFeatures,
                 args=(
                     self.queue_features,
                     self.ftc,
-                )
-            )
+                ),
+            ),
         ]
 
         for p in processes:
@@ -112,15 +104,11 @@ def run(self) -> None:
         for p in processes:
             p.join()
 
-    def get_data(
-        self,
-        queue_raw: multiprocessing.Queue,
-        ftc: FieldTrip.Client
+    def _get_data(
+        self, queue_raw: multiprocessing.Queue, ftc: FieldTrip.Client
     ) -> np.array:
-
         def on_press(key):
-            print('{0} pressed'.format(
-                key))
+            print("{0} pressed".format(key))
 
         def on_release(key):
             if key == Key.esc:
@@ -129,26 +117,24 @@ def on_release(key):
                 queue_raw.put(None)
                 return False
 
-        listener = Listener(
-            on_press=on_press,
-            on_release=on_release)
+        listener = Listener(on_press=on_press, on_release=on_release)
         listener.start()
         last_batch = 0
 
         while listener.is_alive() is True:
             # read new data
-            print('Trying to read last sample...')
-            #index = H.nSamples - 1
+            print("Trying to read last sample...")
+            # index = H.nSamples - 1
             ieeg_batch = ftc.getData()
 
-            #number_repeat = 10
-            #val = timeit.timeit(
+            # number_repeat = 10
+            # val = timeit.timeit(
             #    lambda: ftc.getData(),
             #    number=number_repeat
-            #) / number_repeat
+            # ) / number_repeat
 
             ieeg_batch = ieeg_batch[-128:, :2].T  # take last, #1: data
-            # check if time stamp changed 
+            # check if time stamp changed
             if np.array_equal(ieeg_batch, last_batch) is False:
                 last_batch = ieeg_batch
                 queue_raw.put(ieeg_batch)
@@ -159,7 +145,7 @@ def on_release(key):
     def calcFeatures(
         self,
         queue_raw: multiprocessing.Queue,
-        queue_features: multiprocessing.Queue
+        queue_features: multiprocessing.Queue,
     ):
 
         FLAG_STOP = False
@@ -169,16 +155,13 @@ def calcFeatures(
                 queue_features.put(None)
                 FLAG_STOP = True
             else:
-                feature_series = \
-                    self.run_analysis.process_data(ieeg_batch)
+                feature_series = self.run_analysis.process_data(ieeg_batch)
                 feature_series = self._add_timestamp(feature_series)
                 print("calc features")
                 queue_features.put(feature_series)
 
     def sendFeatures(
-        self,
-        queue_features: multiprocessing.Queue,
-        ftc: FieldTrip.Client
+        self, queue_features: multiprocessing.Queue, ftc: FieldTrip.Client
     ):
 
         idx = 0
@@ -197,8 +180,7 @@ def sendFeatures(
                     idx += 1
                 else:
                     self.feature_arr = self.feature_arr.append(
-                        features,
-                        ignore_index=True
+                        features, ignore_index=True
                     )
 
                 print("length of features:" + str(len(self.feature_arr)))
@@ -206,22 +188,20 @@ def sendFeatures(
                 # channel names are 1. data 2. ident 3. timestamp
                 # put at the end of the buffer the calculated features of the data channel
                 to_send = np.zeros([H.nSamples, H.nChannels])
-                #to_send[-features.shape[0]:, IDENT_FEATURES] = np.array(features)
-                to_send[-features.shape[0]:, 1] = np.array(features)
-                #ftc.putData(to_send)
+                # to_send[-features.shape[0]:, IDENT_FEATURES] = np.array(features)
+                to_send[-features.shape[0] :, 1] = np.array(features)
+                # ftc.putData(to_send)
 
-                #number_repeat = 10
-                #val = timeit.timeit(
+                # number_repeat = 10
+                # val = timeit.timeit(
                 #    lambda: ftc.putData(to_send),
                 #    number=number_repeat
-                #) / number_repeat
+                # ) / number_repeat
 
-    def _add_timestamp(self, feature_series: pd.Series, idx: int = None) -> pd.Series:
+    def _add_timestamp(
+        self, feature_series: pd.Series, idx: int = None
+    ) -> pd.Series:
 
         feature_series["time"] = time.time()  # UNIX Timestamp
 
         return feature_series
-
-    def _add_coordinates(self):
-        """Lateron add here method for providing coordinates"""
-        pass
\ No newline at end of file
diff --git a/py_neuromodulation/nm_RealTimeClientStreamSingleProcess.py b/py_neuromodulation/nm_RealTimeClientStreamSingleProcess.py
index 610a8a74..4862fc44 100644
--- a/py_neuromodulation/nm_RealTimeClientStreamSingleProcess.py
+++ b/py_neuromodulation/nm_RealTimeClientStreamSingleProcess.py
@@ -8,29 +8,31 @@
 import time
 import timeit
 from pynput.keyboard import Key, Listener
-#import pylsl
+
+# import pylsl
 #  clone TMSI-Python-Interface from here: https://gitlab.com/tmsi/tmsi-python-interface
 
-from py_neuromodulation import \
-    (nm_projection,
+from py_neuromodulation import (
+    nm_projection,
     nm_rereference,
     nm_run_analysis,
     nm_features,
     nm_resample,
     nm_stream,
     nm_test_settings,
-    FieldTrip)
+    FieldTrip,
+)
 
 
 class RealTimePyNeuro(nm_stream.PNStream):
 
-    queue_raw:multiprocessing.Queue = multiprocessing.Queue(1)
-    queue_features:multiprocessing.Queue = multiprocessing.Queue(1)
-    ftc:FieldTrip.Client = None
+    queue_raw: multiprocessing.Queue = multiprocessing.Queue(1)
+    queue_features: multiprocessing.Queue = multiprocessing.Queue(1)
+    ftc: FieldTrip.Client = None
     filename: str
-    esc_key_received : bool = False
-    use_FieldTripClient : bool = True
-    feature_count_idx : int = 0
+    esc_key_received: bool = False
+    use_FieldTripClient: bool = True
+    feature_count_idx: int = 0
 
     def __init__(
         self,
@@ -42,23 +44,23 @@ def __init__(
         VERBOSE: bool = True,
         fs: int = 128,
         line_noise: int = 50,
-        use_FieldTripClient : bool = True
+        use_FieldTripClient: bool = True,
     ) -> None:
 
-        super().__init__(PATH_SETTINGS=PATH_SETTINGS,
-            PATH_NM_CHANNELS=PATH_NM_CHANNELS,
-            PATH_OUT=PATH_OUT,
-            PATH_GRIDS=PATH_GRIDS,
-            VERBOSE=VERBOSE)
+        super().__init__(
+            path_settings=PATH_SETTINGS,
+            path_nm_channels=PATH_NM_CHANNELS,
+            path_out=PATH_OUT,
+            path_grids=PATH_GRIDS,
+            verbose=VERBOSE,
+        )
 
         self.set_fs(fs)
         self.set_linenoise(line_noise)
         self.filename = filename
         self.use_FieldTripClient = use_FieldTripClient
 
-        self.nm_channels = self._get_nm_channels(
-            self.PATH_NM_CHANNELS
-        )
+        self.nm_channels = self._get_nm_channels(self.path_nm_channels)
 
         # leave out coordinate setting for now
 
@@ -86,29 +88,28 @@ def init_fieldtrip(PORT: int = 1972) -> FieldTrip.Client:
         """
         ftc = FieldTrip.Client()
         # Python FieldTripBuffer https://www.fieldtriptoolbox.org/development/realtime/buffer_python/
-        ftc.connect('localhost', port=PORT)   # might throw IOError
+        ftc.connect("localhost", port=PORT)  # might throw IOError
         H = ftc.getHeader()
 
         if H is None:
-            print('Failed to retrieve header!')
+            print("Failed to retrieve header!")
             sys.exit(1)
 
         print(H)
         print(H.labels)
 
         return ftc
-    
+
     @staticmethod
-    def init_lsl(wait_max : int = 10, buffer_size : int = 1000):
+    def init_lsl(wait_max: int = 10, buffer_size: int = 1000):
 
-        #streams = pylsl.resolve_streams(wait_time=min(0.1, wait_max))
-        #print("Stream found")
-        #return pylsl.StreamInlet(info=streams[0], max_buflen=buffer_size)
+        # streams = pylsl.resolve_streams(wait_time=min(0.1, wait_max))
+        # print("Stream found")
+        # return pylsl.StreamInlet(info=streams[0], max_buflen=buffer_size)
         ...
 
     def run(self) -> None:
-        """Start get_data, calcFeatures and sendFeature processes
-        """
+        """Start get_data, calcFeatures and sendFeature processes"""
 
         self._set_run()
 
@@ -117,8 +118,8 @@ def run(self) -> None:
         while self.listener.is_alive() is True:
             ieeg_batch = self.get_data_client()
             ieeg_batch = ieeg_batch[-128:, :2]  # take last, #1: data
-            #ieeg_batch = np.random.random([128, 2]).T  # channels, samples
-            # check if time stamp changed 
+            # ieeg_batch = np.random.random([128, 2]).T  # channels, samples
+            # check if time stamp changed
             if np.array_equal(ieeg_batch, last_batch) is False:
 
                 last_batch = ieeg_batch
@@ -133,12 +134,12 @@ def run(self) -> None:
 
         self.disconnect()
 
-    def get_data_lsl(self, max_samples: int = 1000, timeout : int = 5):
-        #samples, _ = self.lsl_client.pull_chunk(
+    def get_data_lsl(self, max_samples: int = 1000, timeout: int = 5):
+        # samples, _ = self.lsl_client.pull_chunk(
         #    max_samples=max_samples,
         #    timeout=timeout
-        #)
-        #return np.vstack(samples).T
+        # )
+        # return np.vstack(samples).T
         ...
 
     def send_data_lsl(self, features: pd.Series):
@@ -147,12 +148,12 @@ def send_data_lsl(self, features: pd.Series):
         for reference implementation
         """
 
-        #info = pylsl.StreamInfo(name='py_nm', type=None,
+        # info = pylsl.StreamInfo(name='py_nm', type=None,
         #                        channel_count=None,
         #                        nominal_srate=None,
         #                        channel_format='float32', source_id=None)
-        #outlet = pylsl.StreamOutlet(info)
-        #outlet.push_sample(features.to_numpy())
+        # outlet = pylsl.StreamOutlet(info)
+        # outlet.push_sample(features.to_numpy())
         ...
 
     def disconnect_lsl(self):
@@ -164,19 +165,15 @@ def get_data_FieldTripClient(self):
     def send_data_FieldTripClient(self, features: pd.Series):
 
         # H = self.ftc_send.getHeader()  # retrieving header is not necessary for sending
-        self.ftc_send.putData(np.random.random([1,2]))
-        #self.ftc_send.putData(np.expand_dims(np.array(features), axis=0))
+        self.ftc_send.putData(np.random.random([1, 2]))
+        # self.ftc_send.putData(np.expand_dims(np.array(features), axis=0))
 
     def disconnect_FieldTripClient(self):
         self.ftc.disconnect()
 
-    def init_keyboard_listener(
-        self
-    ):
-
+    def init_keyboard_listener(self):
         def on_press(key):
-            print('{0} pressed'.format(
-                key))
+            print("{0} pressed".format(key))
 
         def on_release(key):
             if key == Key.esc:
@@ -187,11 +184,9 @@ def on_release(key):
                 self._sendFeatures(None)
                 return False
 
-        self.listener = Listener(
-            on_press=on_press,
-            on_release=on_release)
+        self.listener = Listener(on_press=on_press, on_release=on_release)
         self.listener.start()
-    
+
     def _sendFeatures(self, features: pd.Series):
         if features is None:
             self.save_after_stream(self.filename)
@@ -203,8 +198,7 @@ def _sendFeatures(self, features: pd.Series):
                 self.feature_count_idx += 1
             else:
                 self.feature_arr = self.feature_arr.append(
-                    features,
-                    ignore_index=True
+                    features, ignore_index=True
                 )
 
             print("length of features:" + str(len(self.feature_arr)))
@@ -215,7 +209,9 @@ def _sendFeatures(self, features: pd.Series):
 
         return False
 
-    def _add_timestamp(self, feature_series: pd.Series, idx: int = None) -> pd.Series:
+    def _add_timestamp(
+        self, feature_series: pd.Series, idx: int = None
+    ) -> pd.Series:
 
         feature_series["time"] = time.time()  # UNIX Timestamp
 
diff --git a/py_neuromodulation/nm_analysis.py b/py_neuromodulation/nm_analysis.py
index bbf9d929..d5c9acd3 100644
--- a/py_neuromodulation/nm_analysis.py
+++ b/py_neuromodulation/nm_analysis.py
@@ -54,19 +54,21 @@ def __init__(
             self.feature_file = feature_file
 
         FILE_BASENAME = Path(self.feature_file).stem
-        PATH_READ_FILE = str(Path(self.feature_dir, FILE_BASENAME, FILE_BASENAME))
+        PATH_READ_FILE = str(
+            Path(self.feature_dir, FILE_BASENAME, FILE_BASENAME)
+        )
 
         self.settings = nm_IO.read_settings(PATH_READ_FILE)
         self.sidecar = nm_IO.read_sidecar(PATH_READ_FILE)
         self.fs = self.sidecar["fs"]
-        self.line_noise = self.sidecar["line_noise"]
         self.nm_channels = nm_IO.read_nm_channels(PATH_READ_FILE)
         self.feature_arr = nm_IO.read_features(PATH_READ_FILE)
 
         self.ch_names = self.nm_channels.new_name
         self.used_chs = list(
             self.nm_channels[
-                (self.nm_channels["target"] == 0) & (self.nm_channels["used"] == 1)
+                (self.nm_channels["target"] == 0)
+                & (self.nm_channels["used"] == 1)
             ]["new_name"]
         )
         self.ch_names_ECOG = self.nm_channels.query(
@@ -78,11 +80,16 @@ def __init__(
 
         self.label_name = self._get_target_ch()
         self.label = self.read_target_ch(
-            self.feature_arr, self.label_name, binarize=binarize_label, binarize_th=0.3
+            self.feature_arr,
+            self.label_name,
+            binarize=binarize_label,
+            binarize_th=0.3,
         )
 
     def _get_target_ch(self) -> str:
-        target_names = list(self.nm_channels[self.nm_channels["target"] == 1]["name"])
+        target_names = list(
+            self.nm_channels[self.nm_channels["target"] == 1]["name"]
+        )
         target_clean = [
             target_name
             for target_name in target_names
@@ -148,7 +155,9 @@ def filter_features(
 
         if list_feature_keywords is not None:
             feature_select = [
-                f for f in feature_select if any(x in f for x in list_feature_keywords)
+                f
+                for f in feature_select
+                if any(x in f for x in list_feature_keywords)
             ]
 
             if (
@@ -171,9 +180,13 @@ def set_target_ch(self, ch_name: str) -> None:
 
     def plot_cort_projection(self) -> None:
         if self.sidecar["sess_right"]:
-            ecog_strip = np.array(self.sidecar["coords"]["cortex_right"]["positions"])
+            ecog_strip = np.array(
+                self.sidecar["coords"]["cortex_right"]["positions"]
+            )
         else:
-            ecog_strip = np.array(self.sidecar["coords"]["cortex_left"]["positions"])
+            ecog_strip = np.array(
+                self.sidecar["coords"]["cortex_left"]["positions"]
+            )
         self.nmplotter.plot_cortex(
             grid_cortex=np.array(self.sidecar["grid_cortex"])
             if "grid_cortex" in self.sidecar
@@ -221,7 +234,9 @@ def plot_target_averaged_channel(
     ):
 
         filtered_df = self.feature_arr[
-            self.filter_features(self.feature_arr.columns, ch, list_feature_keywords)
+            self.filter_features(
+                self.feature_arr.columns, ch, list_feature_keywords
+            )
         ]
 
         data = np.expand_dims(np.array(filtered_df), axis=1)
@@ -286,7 +301,9 @@ def plot_subject_grid_ch_performance(
 
         for ch in ch_:
             if "grid" not in ch and "combined" not in ch:
-                ecog_coords_strip.append(performance_dict[subject_name][ch]["coord"])
+                ecog_coords_strip.append(
+                    performance_dict[subject_name][ch]["coord"]
+                )
                 ecog_strip_performance.append(
                     performance_dict[subject_name][ch]["performance_test"]
                 )
@@ -303,7 +320,9 @@ def plot_subject_grid_ch_performance(
             grid_cortex=np.array(self.sidecar["grid_cortex"])
             if "grid_cortex" in self.sidecar
             else None,
-            ecog_strip=ecog_coords_strip if len(ecog_coords_strip) > 0 else None,
+            ecog_strip=ecog_coords_strip
+            if len(ecog_coords_strip) > 0
+            else None,
             grid_color=grid_performance if len(grid_performance) > 0 else None,
             strip_color=ecog_strip_performance
             if len(ecog_strip_performance) > 0
@@ -353,13 +372,17 @@ def get_epochs(data, y_, epoch_len, sfreq, threshold=0):
 
         ind_mov = ind_mov[low_limit & up_limit]
 
-        epoch_ = np.zeros([ind_mov.shape[0], epoch_lim, data.shape[1], data.shape[2]])
+        epoch_ = np.zeros(
+            [ind_mov.shape[0], epoch_lim, data.shape[1], data.shape[2]]
+        )
 
         y_arr = np.zeros([ind_mov.shape[0], int(epoch_lim)])
 
         for idx, i in enumerate(ind_mov):
 
-            epoch_[idx, :, :, :] = data[i - epoch_lim // 2 : i + epoch_lim // 2, :, :]
+            epoch_[idx, :, :, :] = data[
+                i - epoch_lim // 2 : i + epoch_lim // 2, :, :
+            ]
 
             y_arr[idx, :] = y_[i - epoch_lim // 2 : i + epoch_lim // 2]
 
@@ -384,6 +407,8 @@ def set_decoder(
         time_stack_n_samples=5,
         use_nested_cv=False,
         VERBOSE=False,
+        undersampling=False,
+        oversampling=False,
     ):
         if decoder is not None:
             self.decoder = decoder
@@ -408,6 +433,8 @@ def set_decoder(
                 time_stack_n_samples=time_stack_n_samples,
                 VERBOSE=VERBOSE,
                 use_nested_cv=use_nested_cv,
+                undersampling=undersampling,
+                oversampling=oversampling,
             )
 
     def run_ML_model(
@@ -572,10 +599,14 @@ def set_score(
                     obj_write[key_set] = val
 
                 set_score(
-                    key_set="performance_test", key_get="score_test", take_mean=True
+                    key_set="performance_test",
+                    key_get="score_test",
+                    take_mean=True,
                 )
                 set_score(
-                    key_set="performance_train", key_get="score_train", take_mean=True
+                    key_set="performance_train",
+                    key_get="score_train",
+                    take_mean=True,
                 )
 
                 if "coef" in obj_read:
@@ -598,16 +629,24 @@ def set_score(
                         take_mean=True,
                     )
                     set_score(
-                        key_set="fprate_test", key_get="fprate_test", take_mean=True
+                        key_set="fprate_test",
+                        key_get="fprate_test",
+                        take_mean=True,
                     )
                     set_score(
-                        key_set="fprate_train", key_get="fprate_train", take_mean=True
+                        key_set="fprate_train",
+                        key_get="fprate_train",
+                        take_mean=True,
                     )
                     set_score(
-                        key_set="tprate_test", key_get="tprate_test", take_mean=True
+                        key_set="tprate_test",
+                        key_get="tprate_test",
+                        take_mean=True,
                     )
                     set_score(
-                        key_set="tprate_train", key_get="tprate_train", take_mean=True
+                        key_set="tprate_train",
+                        key_get="tprate_train",
+                        take_mean=True,
                     )
 
                 if read_bay_opt_params is True:
@@ -624,7 +663,10 @@ def set_score(
 
             read_ML_performances(obj_read, obj_write)
 
-            if len([key_ for key_ in obj_read.keys() if "InnerCV_" in key_]) > 0:
+            if (
+                len([key_ for key_ in obj_read.keys() if "InnerCV_" in key_])
+                > 0
+            ):
                 read_ML_performances(obj_read, obj_write, set_inner_CV_res=True)
 
         if read_channels:
@@ -635,21 +677,24 @@ def set_score(
 
                 performance_dict[subject_name][ch] = {}
 
-                if len(self.sidecar["coords"]) > 0:  # check if coords are empty
-
-                    coords_exist = False
-                    for cortex_loc in self.sidecar["coords"].keys():
-                        for ch_name_coord_idx, ch_name_coord in enumerate(
-                            self.sidecar["coords"][cortex_loc]["ch_names"]
-                        ):
-                            if ch.startswith(ch_name_coord):
-                                coords = self.sidecar["coords"][cortex_loc][
-                                    "positions"
-                                ][ch_name_coord_idx]
-                                coords_exist = True  # optimally break out of the two loops...
-                    if coords_exist is False:
-                        coords = None
-                    performance_dict[subject_name][ch]["coord"] = coords
+                if "coords" in self.sidecar:
+                    if (
+                        len(self.sidecar["coords"]) > 0
+                    ):  # check if coords are empty
+
+                        coords_exist = False
+                        for cortex_loc in self.sidecar["coords"].keys():
+                            for ch_name_coord_idx, ch_name_coord in enumerate(
+                                self.sidecar["coords"][cortex_loc]["ch_names"]
+                            ):
+                                if ch.startswith(ch_name_coord):
+                                    coords = self.sidecar["coords"][cortex_loc][
+                                        "positions"
+                                    ][ch_name_coord_idx]
+                                    coords_exist = True  # optimally break out of the two loops...
+                        if coords_exist is False:
+                            coords = None
+                        performance_dict[subject_name][ch]["coord"] = coords
                 write_CV_res_in_performance_dict(
                     ML_res.ch_ind_results[ch],
                     performance_dict[subject_name][ch],
@@ -686,9 +731,11 @@ def set_score(
                     gp_str = grid_type + str(grid_point)
 
                     performance_dict[subject_name][gp_str] = {}
-                    performance_dict[subject_name][gp_str]["coord"] = self.sidecar[
-                        "grid_" + project_settings.split("_")[1]
-                    ][grid_point]
+                    performance_dict[subject_name][gp_str][
+                        "coord"
+                    ] = self.sidecar["grid_" + project_settings.split("_")[1]][
+                        grid_point
+                    ]
 
                     if gp_str in ML_res.active_gridpoints:
                         write_CV_res_in_performance_dict(
diff --git a/py_neuromodulation/nm_cohortwrapper.py b/py_neuromodulation/nm_cohortwrapper.py
index a4ce42ef..c57fd8b0 100644
--- a/py_neuromodulation/nm_cohortwrapper.py
+++ b/py_neuromodulation/nm_cohortwrapper.py
@@ -23,7 +23,8 @@
 import nibabel as nib
 
 import py_neuromodulation
-from py_neuromodulation import nm_decode, nm_analysis, nm_BidsStream, nm_IO
+from py_neuromodulation import nm_decode, nm_analysis, nm_IO
+from py_neuromodulation import nm_stream_offline
 
 
 class CohortRunner:
@@ -86,7 +87,8 @@ def __init__(
         self.target_keywords = target_keywords
         self.get_movement_detection_rate = get_movement_detection_rate
         self.grid_cortex = pd.read_csv(
-            os.path.join(py_neuromodulation.__path__[0], "grid_cortex.tsv"), sep="\t"
+            os.path.join(py_neuromodulation.__path__[0], "grid_cortex.tsv"),
+            sep="\t",
         ).to_numpy()
 
     def init_decoder(self) -> nm_decode.Decoder:
@@ -117,7 +119,7 @@ def multiprocess_pipeline_run_wrapper(self, PATH_RUN):
                 break
 
         if self.run_bids:
-            nm_BIDS = nm_BidsStream.BidsStream(
+            nm_BIDS = nm_stream_offline.BidsStream(
                 PATH_RUN=PATH_RUN,
                 PATH_BIDS=PATH_BIDS,
                 PATH_OUT=PATH_OUT,
@@ -180,7 +182,7 @@ def multiprocess_pipeline_run_wrapper(self, PATH_RUN):
                 estimate_gridpoints=self.estimate_gridpoints,
                 estimate_all_channels_combined=self.estimate_all_channels_combined,
                 save_results=True,
-                output_name="XGBNOMP"
+                output_name="XGBNOMP",
             )
 
         if self.plot_grid_performances:
@@ -234,9 +236,9 @@ def read_cohort_results(self, feature_path, cohort):
                 read_mov_detection_rates=True,
             )
 
-            sub = feature_file[feature_file.find("sub-") : feature_file.find("_ses")][
-                4:
-            ]
+            sub = feature_file[
+                feature_file.find("sub-") : feature_file.find("_ses")
+            ][4:]
             if sub not in performance_out:
                 performance_out[sub] = {}
             performance_out[sub][feature_file] = performance_run[
@@ -254,7 +256,9 @@ def cohort_wrapper_read_cohort(self):
         """Read results for multiple cohorts"""
 
         for cohort in self.cohorts.keys():
-            self.read_cohort_results(os.path.join(self.outpath, cohort), cohort)
+            self.read_cohort_results(
+                os.path.join(self.outpath, cohort), cohort
+            )
 
     def read_all_channels(
         self,
@@ -321,13 +325,19 @@ def read_all_channels(
                 channel_all[cohort][subject_name][ch] = {}
             channel_all[cohort][subject_name][ch][feature_file] = {}
 
-            channel_all[cohort][subject_name][ch][feature_file]["data"] = data_to_read[
-                ch
-            ]
-            channel_all[cohort][subject_name][ch][feature_file]["feature_names"] = [
-                ch_[len(ch) + 1 :] for ch_ in decoder.features.columns if ch in ch_
+            channel_all[cohort][subject_name][ch][feature_file][
+                "data"
+            ] = data_to_read[ch]
+            channel_all[cohort][subject_name][ch][feature_file][
+                "feature_names"
+            ] = [
+                ch_[len(ch) + 1 :]
+                for ch_ in decoder.features.columns
+                if ch in ch_
             ]
-            channel_all[cohort][subject_name][ch][feature_file]["label"] = decoder.label
+            channel_all[cohort][subject_name][ch][feature_file][
+                "label"
+            ] = decoder.label
             channel_all[cohort][subject_name][ch][feature_file][
                 "label_name"
             ] = decoder.label_name
@@ -336,9 +346,11 @@ def read_all_channels(
             lat = "CON"  # Beijing is always contralateral
             # Pittsburgh Subjects
             if (
-                "LEFT" in decoder.label_name and "LEFT" in decoder.features.columns[1]
+                "LEFT" in decoder.label_name
+                and "LEFT" in decoder.features.columns[1]
             ) or (
-                "RIGHT" in decoder.label_name and "RIGHT" in decoder.features.columns[1]
+                "RIGHT" in decoder.label_name
+                and "RIGHT" in decoder.features.columns[1]
             ):
                 lat = "IPS"
 
@@ -372,9 +384,14 @@ def cohort_wrapper_read_all_grid_points(self, read_channels=True):
                 )
 
         if read_channels is True:
-            np.save(os.path.join(self.outpath, "channel_all.npy"), grid_point_all)
+            np.save(
+                os.path.join(self.outpath, "channel_all.npy"), grid_point_all
+            )
         else:
-            np.save(os.path.join(self.outpath, "grid_point_all.npy"), grid_point_all)
+            np.save(
+                os.path.join(self.outpath, "grid_point_all.npy"),
+                grid_point_all,
+            )
 
     @staticmethod
     def rewrite_grid_point_all(d, outpath):
@@ -397,6 +414,8 @@ def rewrite_grid_point_all(d, outpath):
                         if f not in p[gp][cohort][sub].keys():
                             p[gp][cohort][sub][f] = {}
                         for key_ in d[cohort][sub][gp][f].keys():
-                            p[gp][cohort][sub][f][key_] = d[cohort][sub][gp][f][key_]
+                            p[gp][cohort][sub][f][key_] = d[cohort][sub][gp][
+                                f
+                            ][key_]
 
         np.save(os.path.join(outpath, "grid_point_all_re.npy"), p)
diff --git a/py_neuromodulation/nm_decode.py b/py_neuromodulation/nm_decode.py
index a3155e51..b7efc364 100644
--- a/py_neuromodulation/nm_decode.py
+++ b/py_neuromodulation/nm_decode.py
@@ -8,6 +8,7 @@
 from scipy.ndimage import binary_dilation, binary_erosion
 from scipy.ndimage import label as label_ndimage
 from imblearn.over_sampling import RandomOverSampler
+from imblearn.under_sampling import RandomUnderSampler
 import pandas as pd
 import os
 import json
@@ -60,6 +61,7 @@ class Decoder:
     STACK_FEATURES_N_SAMPLES: bool
     time_stack_n_samples: int
     ros: RandomOverSampler = None
+    rus: RandomUnderSampler = None
     VERBOSE: bool = False
     ch_ind_data: dict = {}
     grid_point_ind_data: dict = {}
@@ -102,6 +104,8 @@ def __init__(
         bay_opt_param_space: list = [],
         VERBOSE: bool = False,
         fs: int = None,
+        undersampling: bool = False,
+        oversampling: bool = False,
     ) -> None:
         """Initialize here a feature file for processing
         Read settings.json nm_channels.csv and features.csv
@@ -149,6 +153,8 @@ def __init__(
         self.time_stack_n_samples = time_stack_n_samples
         self.bay_opt_param_space = bay_opt_param_space
         self.VERBOSE = VERBOSE
+        self.undersampling = undersampling
+        self.oversampling = oversampling
 
         self.ch_ind_data = {}
         self.grid_point_ind_data = {}
@@ -158,7 +164,10 @@ def __init__(
         self.gridpoint_ind_results = {}
         self.all_ch_results = {}
 
-        if type(self.model) is discriminant_analysis.LinearDiscriminantAnalysis:
+        if undersampling:
+            self.rus = RandomUnderSampler(random_state=0)
+
+        if oversampling:
             self.ros = RandomOverSampler(random_state=0)
 
     def set_data(self, features, label, label_name, used_chs):
@@ -485,8 +494,10 @@ def evalerror(preds, dtrain):
         else:
 
             # check for LDA; and apply rebalancing
-            if type(model) is discriminant_analysis.LinearDiscriminantAnalysis:
+            if self.oversampling:
                 X_train, y_train = self.ros.fit_resample(X_train, y_train)
+            if self.undersampling:
+                X_train, y_train = self.rus.fit_resample(X_train, y_train)
 
             if type(model) is xgboost.sklearn.XGBClassifier:
                 model.fit(X_train, y_train, eval_metric="logloss")  # to avoid warning
diff --git a/py_neuromodulation/nm_features.py b/py_neuromodulation/nm_features.py
index dbc7dc1b..40424d9f 100644
--- a/py_neuromodulation/nm_features.py
+++ b/py_neuromodulation/nm_features.py
@@ -15,7 +15,7 @@
 
 
 class Features:
-    def __init__(self, s, ch_names, fs, line_noise, verbose=None) -> None:
+    def __init__(self, s, ch_names, fs, verbose=None) -> None:
         """Initialize Feature module
 
         Parameters
@@ -31,7 +31,6 @@ def __init__(self, s, ch_names, fs, line_noise, verbose=None) -> None:
             self.fs = s["raw_resampling_settings"]["resample_freq_hz"]
         else:
             self.fs = fs
-        self.line_noise = line_noise
         self.fband_names = [value for value in s["frequency_ranges_hz"].keys()]
         self.f_ranges = [
             self.s["frequency_ranges_hz"][fband_name]
diff --git a/py_neuromodulation/nm_generator.py b/py_neuromodulation/nm_generator.py
index 11469b52..a6fadc82 100644
--- a/py_neuromodulation/nm_generator.py
+++ b/py_neuromodulation/nm_generator.py
@@ -1,7 +1,11 @@
-from numpy import ceil
+from typing import Iterator
 
+import numpy as np
 
-def ieeg_raw_generator(ieeg_raw, settings: dict, fs: int):
+
+def ieeg_raw_generator(
+    ieeg_raw, settings: dict, fs: int
+) -> Iterator[np.ndarray]:
     """
     This generator function mimics online data acquisition.
     The nm_channels.csv selected raw channels are iteratively sampled with fs.
@@ -15,7 +19,7 @@ def ieeg_raw_generator(ieeg_raw, settings: dict, fs: int):
 
     cnt_fsnew = 0
     offset_time = settings["segment_length_features_ms"]
-    offset_start = ceil(offset_time / 1000 * fs).astype(int)
+    offset_start = np.ceil(offset_time / 1000 * fs).astype(int)
     fs_new = settings["sampling_rate_features_hz"]
 
     for cnt in range(ieeg_raw.shape[1]):
diff --git a/py_neuromodulation/nm_run_analysis.py b/py_neuromodulation/nm_run_analysis.py
index e4a58823..31ba7ab4 100644
--- a/py_neuromodulation/nm_run_analysis.py
+++ b/py_neuromodulation/nm_run_analysis.py
@@ -18,9 +18,9 @@ def __init__(
         self,
         features: nm_features.Features,
         settings: dict,
-        reference: nm_rereference.RT_rereference,
-        projection: nm_projection.Projection,
-        resample: nm_resample.Resample,
+        reference: nm_rereference.RT_rereference | None,
+        projection: nm_projection.Projection | None,
+        resample: nm_resample.Resample | None,
         notch_filter: nm_filter.NotchFilter | None,
         verbose: bool,
         feature_idx: list,
@@ -67,7 +67,6 @@ def __init__(
         self.settings = settings
         self.fs_new = int(settings["sampling_rate_features_hz"])
         self.fs = features.fs
-        self.line_noise = features.line_noise
         self.feature_idx = feature_idx
         self.sample_add = int(self.fs / self.fs_new)
         self.verbose = verbose
diff --git a/py_neuromodulation/nm_stream.py b/py_neuromodulation/nm_stream.py
deleted file mode 100644
index 7a643765..00000000
--- a/py_neuromodulation/nm_stream.py
+++ /dev/null
@@ -1,465 +0,0 @@
-from abc import ABC, abstractmethod
-import pathlib
-from sklearn import base
-import _pickle as cPickle
-import numpy as np
-import pandas as pd
-import os
-from enum import Enum
-from typing import Tuple, Union
-
-from py_neuromodulation import (
-    nm_filter,
-    nm_projection,
-    nm_rereference,
-    nm_run_analysis,
-    nm_features,
-    nm_resample,
-    nm_define_nmchannels,
-    nm_IO,
-    nm_plots,
-    nm_test_settings,
-)
-
-
-class GRIDS(Enum):
-    """Definition of possible projection grid types"""
-
-    CORTEX = "cortex"
-    SUBCORTEX = "subcortex"
-
-
-class PNStream(ABC):
-
-    resample: nm_resample.Resample
-    features: nm_features.Features
-    run_analysis: nm_run_analysis.Run
-    rereference: nm_rereference.RT_rereference
-    notch_filter: nm_filter.NotchFilter
-    projection: nm_projection.Projection
-    settings: dict
-    nm_channels: pd.DataFrame
-    coords: dict = {}
-    fs: Union[int, float]
-    line_noise: float
-    VERBOSE: bool
-    PATH_SETTINGS: str
-    PATH_NM_CHANNELS: str = str()
-    PATH_OUT: str = str()
-    PATH_GRIDS: str = str()
-    feature_arr: pd.DataFrame = pd.DataFrame()
-    CH_NAMES_USED: list
-    CH_TYPES_USED: list
-    FEATURE_IDX: list
-    LABEL_IDX: list
-    grid_cortex: np.array
-    grid_subcortex: np.array
-    sess_right: bool = None
-    feature_add: pd.DataFrame
-    model: base.BaseEstimator
-
-    @abstractmethod
-    def __init__(
-        self,
-        PATH_SETTINGS=os.path.join(
-            pathlib.Path(__file__).parent.resolve(), "nm_settings.json"
-        ),
-        PATH_NM_CHANNELS: str = str(),
-        PATH_OUT: str = os.getcwd(),
-        PATH_GRIDS: str = pathlib.Path(__file__).parent.resolve(),
-        VERBOSE: bool = True,
-    ) -> None:
-
-        self.PATH_SETTINGS = PATH_SETTINGS
-        self.PATH_NM_CHANNELS = PATH_NM_CHANNELS
-        self.PATH_OUT = PATH_OUT
-        self.VERBOSE = VERBOSE
-
-        self.settings = nm_IO.read_settings(self.PATH_SETTINGS)
-
-        if True in [
-            self.settings["methods"]["project_cortex"],
-            self.settings["methods"]["project_subcortex"],
-        ]:
-            self.grid_cortex, self.grid_subcortex = self.get_grids(
-                self.settings, self.PATH_GRIDS, GRIDS
-            )
-
-    @abstractmethod
-    def _add_coordinates(self) -> None:
-        """Write self.coords either from bids or from separate file
-        This method is implemented differently
-        for BIDS and real time data anylsis
-        """
-        pass
-
-    @abstractmethod
-    def get_data(
-        self,
-    ) -> np.array:
-        """Get new data batch from acquisition device or from BIDS"""
-        pass
-
-    @abstractmethod
-    def run(
-        self,
-    ):
-        """In this function data is first acquied iteratively
-        1. self.get_data()
-        2. data processing is called:
-        self.run_analysis.process_data(data) to calculate features
-        3. optional postprocessing
-        e.g. plotting, ML estimation is done
-        """
-        pass
-
-    @abstractmethod
-    def _add_timestamp(
-        self, feature_series: pd.Series, idx: int = None
-    ) -> pd.Series:
-        """Add to feature_series "time" keyword
-        For Bids specify with fs_features, for real time analysis with current time stamp
-        """
-        pass
-
-    def load_model(self, model_name: str):
-        """Load sklearn model, that utilizes predict"""
-        with open(model_name, "rb") as fid:
-            self.model = cPickle.load(fid)
-
-    def _set_run(self):
-        """Initialize preprocessing, and feature estimation modules"""
-
-        nm_test_settings.test_settings(self.settings, self.nm_channels)
-
-        (
-            self.CH_NAMES_USED,
-            self.CH_TYPES_USED,
-            self.FEATURE_IDX,
-            self.LABEL_IDX,
-        ) = self._get_ch_info(self.nm_channels)
-
-        self.features = self._set_features(
-            self.settings,
-            self.CH_NAMES_USED,
-            self.fs,
-            self.line_noise,
-            self.VERBOSE,
-        )
-
-        self.resample = self._set_resampling(
-            settings=self.settings,
-            sfreq_old=self.fs,
-        )
-
-        self.rereference, self.nm_channels = self._set_rereference(
-            self.settings, self.nm_channels
-        )
-
-        self.notch_filter = self._set_notch_filter(
-            settings=self.settings,
-            fs=self.fs
-            if self.settings["methods"]["raw_resampling"] is False
-            else self.settings["raw_resampling_settings"]["resample_freq_hz"],
-            line_noise=self.line_noise,
-        )
-
-        self.projection = self._get_projection(self.settings, self.nm_channels)
-        if "cortex_left" in self.coords or "cortex_right" in self.coords:
-            if (
-                self.projection is not None
-                or len(self.coords["cortex_left"]["positions"])
-                or len(self.coords["cortex_right"]["positions"])
-            ):
-                self.sess_right = self._get_sess_lat(self.coords)
-            else:
-                self.sess_right = None
-        else:
-            self.sess_right = None
-
-        self.run_analysis = nm_run_analysis.Run(
-            features=self.features,
-            settings=self.settings,
-            reference=self.rereference,
-            projection=self.projection,
-            resample=self.resample,
-            notch_filter=self.notch_filter,
-            verbose=self.VERBOSE,
-            feature_idx=self.FEATURE_IDX,
-        )
-
-    def _set_features(
-        self,
-        settings: dict,
-        CH_NAMES_USED: list,
-        fs: Union[int, float],
-        line_noise: int,
-        VERBOSE: bool,
-    ) -> None:
-        """initialize feature class from settings"""
-        return nm_features.Features(
-            settings, CH_NAMES_USED, fs, line_noise, VERBOSE
-        )
-
-    def set_fs(self, fs: Union[int, float]) -> None:
-        self.fs = fs
-
-    def _set_rereference(
-        self, settings: dict, nm_channels: pd.DataFrame
-    ) -> tuple[nm_rereference.RT_rereference, pd.DataFrame]:
-        """Initialize nm_rereference and update nm_channels
-        nm_channels are updated if no rereferencing is specified
-
-        Parameters
-        ----------
-        settings : dict
-            [description]
-        nm_channels : pd.DataFrame
-            [description]
-
-        Returns
-        -------
-        Tuple
-            nm_rereference object, updated nm_channels DataFrame
-        """
-        if settings["methods"]["re_referencing"] is True:
-            rereference = nm_rereference.RT_rereference(nm_channels)
-        else:
-            rereference = None
-            # reset nm_channels from default values
-            nm_channels["rereference"] = None
-            nm_channels["new_name"] = nm_channels["name"]
-        return rereference, nm_channels
-
-    def _set_resampling(
-        self, settings: dict, sfreq_old: int | float
-    ) -> nm_resample.Resample | None:
-        """Initialize Resampling
-
-        Parameters
-        ----------
-        settings : dict
-        fs : int
-
-        Returns
-        -------
-        nm_resample.Resample
-        """
-        ### We can actually leave out "is True", because this is implied by the
-        ### if statement. Actually leaving out "is True" or ""== True" makes
-        ### the code faster
-        if settings["methods"]["raw_resampling"]:
-            resample = nm_resample.Resample(
-                sfreq_old=sfreq_old,
-                sfreq_new=settings["raw_resampling_settings"][
-                    "resample_freq_hz"
-                ],
-            )
-        else:
-            resample = None
-        return resample
-
-    def _set_notch_filter(
-        self,
-        settings: dict,
-        fs: Union[int, float],
-        line_noise: int,
-        notch_widths: int = 3,
-        trans_bandwidth: int = 15,
-    ) -> nm_filter.NotchFilter | None:
-
-        if settings["methods"]["notch_filter"]:
-            return nm_filter.NotchFilter(
-                sfreq=fs,
-                line_noise=line_noise,
-                notch_widths=notch_widths,
-                trans_bandwidth=trans_bandwidth,
-            )
-        return
-
-    def set_linenoise(self, line_noise: int) -> None:
-        self.line_noise = line_noise
-
-    @staticmethod
-    def get_grids(
-        settings: dict(), PATH_GRIDS: str, GRID_TYPE: GRIDS
-    ) -> Tuple:
-        """Read settings specified grids
-
-        Parameters
-        ----------
-        settings : dict
-        PATH_GRIDS : str
-        GRID_TYPE : GRIDS
-
-        Returns
-        -------
-        Tuple
-            grid_cortex, grid_subcortex,
-            might be None if not specified in settings
-        """
-
-        if settings["methods"]["project_cortex"] is True:
-            grid_cortex = nm_IO.read_grid(PATH_GRIDS, GRID_TYPE.CORTEX)
-        else:
-            grid_cortex = None
-        if settings["methods"]["project_subcortex"] is True:
-            grid_subcortex = nm_IO.read_grid(PATH_GRIDS, GRID_TYPE.SUBCORTEX)
-        else:
-            grid_subcortex = None
-        return grid_cortex, grid_subcortex
-
-    def _get_projection(
-        self, settings: dict, nm_channels: pd.DataFrame
-    ) -> nm_projection.Projection:
-        """Return projection of used coordinated and grids"""
-
-        if any(
-            (
-                settings["methods"]["project_cortex"],
-                settings["methods"]["project_subcortex"],
-            )
-        ):
-            projection = nm_projection.Projection(
-                settings=settings,
-                grid_cortex=self.grid_cortex,
-                grid_subcortex=self.grid_subcortex,
-                coords=self.coords,
-                nm_channels=nm_channels,
-                plot_projection=False,
-            )
-        else:
-            projection = None
-        return projection
-
-    @staticmethod
-    def _get_ch_info(nm_channels: pd.DataFrame):
-        """Get used feature and label info from nm_channels"""
-
-        CH_NAMES_USED = nm_channels[nm_channels["used"] == 1][
-            "new_name"
-        ].tolist()
-        CH_TYPES_USED = nm_channels[nm_channels["used"] == 1]["type"].tolist()
-
-        # used channels for feature estimation
-        FEATURE_IDX = np.where(nm_channels["used"] & ~nm_channels["target"])[
-            0
-        ].tolist()
-
-        # If multiple targets exist, select only the first
-        LABEL_IDX = np.where(nm_channels["target"] == 1)[0]
-
-        return CH_NAMES_USED, CH_TYPES_USED, FEATURE_IDX, LABEL_IDX
-
-    @staticmethod
-    def _get_nm_channels(PATH_NM_CHANNELS: str, **kwargs) -> None:
-        """Read nm_channels from path or specify via BIDS arguments.
-        Nexessary parameters are then
-        ch_names (list),
-        ch_types (list),
-        bads (list)
-        used_types (list)
-        target_keywords (list)
-        reference Union[list, str]
-        """
-
-        if PATH_NM_CHANNELS and os.path.isfile(PATH_NM_CHANNELS):
-            nm_channels = pd.read_csv(PATH_NM_CHANNELS)
-        elif (
-            len(
-                [
-                    1
-                    for l in [
-                        kwargs.get("ch_names", None),
-                        kwargs.get("ch_types", None),
-                    ]
-                    if len(l) == 0
-                ]
-            )
-            == 0
-        ):
-
-            nm_channels = nm_define_nmchannels.set_channels(
-                ch_names=kwargs.get("ch_names"),
-                ch_types=kwargs.get("ch_types"),
-                bads=kwargs.get("bads"),
-                used_types=kwargs.get("used_types"),
-                target_keywords=kwargs.get("target_keywords"),
-                reference=kwargs.get("reference"),
-            )
-        return nm_channels
-
-    @staticmethod
-    def _get_sess_lat(coords):
-        if len(coords["cortex_left"]["positions"]) == 0:
-            sess_right = True
-        elif len(coords["cortex_right"]["positions"]) == 0:
-            sess_right = False
-        return sess_right
-
-    def save_sidecar(self, folder_name: str):
-        """Save sidecar incuding fs, line_noise, coords, sess_right to
-        PATH_OUT and subfolder 'folder_name'"""
-
-        sidecar = {
-            "original_fs": self.fs,
-            "fs": self.run_analysis.fs,
-            "line_noise": self.line_noise,
-            "ch_names": self.features.ch_names,
-            "coords": self.coords,
-            "sess_right": self.sess_right,
-        }
-
-        if self.settings["methods"]["project_cortex"]:
-            sidecar["grid_cortex"] = self.grid_cortex
-            sidecar["proj_matrix_cortex"] = self.projection.proj_matrix_cortex
-        if self.settings["methods"]["project_subcortex"]:
-            sidecar["grid_subcortex"] = self.grid_subcortex
-            sidecar[
-                "proj_matrix_subcortex"
-            ] = self.projection.proj_matrix_subcortex
-
-        nm_IO.save_sidecar(sidecar, self.PATH_OUT, folder_name)
-
-    def save_settings(self, folder_name: str):
-        nm_IO.save_settings(self.settings, self.PATH_OUT, folder_name)
-
-    def save_nm_channels(self, folder_name: str):
-        nm_IO.save_nmchannels(self.nm_channels, self.PATH_OUT, folder_name)
-
-    def save_features(self, folder_name: str):
-        nm_IO.save_features(self.feature_arr, self.PATH_OUT, folder_name)
-
-    def save_after_stream(
-        self, folder_name: str, save_features: bool = True
-    ) -> None:
-        """Save features, settings, nm_channels and sidecar after run"""
-
-        # create derivate folder_name output folder if doesn't exist
-        if os.path.exists(os.path.join(self.PATH_OUT, folder_name)) is False:
-            os.makedirs(os.path.join(self.PATH_OUT, folder_name))
-
-        self.save_sidecar(folder_name)
-        if save_features is True:
-            self.save_features(folder_name)
-
-        self.save_settings(folder_name)
-
-        self.save_nm_channels(folder_name)
-
-    def plot_cortical_projection(self):
-        """plot projection of cortical grid electrodes on cortex"""
-
-        if hasattr(self, "features") is False:
-            self._set_run()
-
-        nmplotter = nm_plots.NM_Plot(
-            ecog_strip=self.projection.ecog_strip
-            if self.projection is not None
-            else None,
-            grid_cortex=self.projection.grid_cortex
-            if self.projection is not None
-            else None,
-            sess_right=self.sess_right,
-        )
-        nmplotter.plot_cortex(set_clim=False)
diff --git a/py_neuromodulation/nm_stream_abc.py b/py_neuromodulation/nm_stream_abc.py
new file mode 100644
index 00000000..50815b3e
--- /dev/null
+++ b/py_neuromodulation/nm_stream_abc.py
@@ -0,0 +1,484 @@
+"""Module that contains PNStream ABC."""
+import os
+import pathlib
+from abc import ABC, abstractmethod
+from enum import Enum
+
+import _pickle as cPickle
+import numpy as np
+import pandas as pd
+from sklearn import base
+
+from py_neuromodulation import (
+    nm_features,
+    nm_filter,
+    nm_IO,
+    nm_plots,
+    nm_projection,
+    nm_rereference,
+    nm_resample,
+    nm_run_analysis,
+    nm_test_settings,
+)
+
+_PathLike = str | os.PathLike
+
+
+class GRIDS(Enum):
+    """Definition of possible projection grid types"""
+
+    CORTEX = "cortex"
+    SUBCORTEX = "subcortex"
+
+
+class PNStream(ABC):
+
+    settings: dict | _PathLike
+    nm_channels: pd.DataFrame | _PathLike
+    run_analysis: nm_run_analysis.Run | None
+    rereference: nm_rereference.RT_rereference | None
+    resample: nm_resample.Resample | None
+    notch_filter: nm_filter.NotchFilter | None
+    features: nm_features.Features
+    coords: dict
+    projection: nm_projection.Projection | None
+    sfreq: int | None
+    path_out: _PathLike | None
+    path_grids: _PathLike | None
+    model: base.BaseEstimator
+    verbose: bool
+
+    def __init__(
+        self,
+        nm_channels: pd.DataFrame | _PathLike,
+        settings: dict | _PathLike | None = None,
+        path_out: _PathLike | None = None,
+        path_grids: _PathLike | None = None,
+        verbose: bool = True,
+        folder_name: str = "sub",
+    ) -> None:
+        if settings is None:
+            settings = (
+                pathlib.Path(__file__).parent.resolve() / "nm_settings.json"
+            )
+        if path_out is None:
+            path_out = os.get_cwd()
+        if path_grids is None:
+            path_grids = pathlib.Path(__file__).parent.resolve()
+
+        self.path_out = path_out
+        self.path_grids = path_grids
+        self.verbose = verbose
+        self.settings = self._get_settings(settings)
+        self.nm_channels = self._get_nm_channels(nm_channels)
+        self.run_analysis = None
+        self.coords = {}
+        self.sess_right = None
+        self.projection = None
+        self.sfreq = None
+        self.model = None
+        self.folder_name = folder_name
+
+    @abstractmethod
+    def run(self):
+        """In this function data is first acquired iteratively
+        1. self.get_data()
+        2. data processing is called:
+        self.run_analysis.process_data(data) to calculate features
+        3. optional postprocessing
+        e.g. plotting, ML estimation is done
+        """
+
+    @abstractmethod
+    def _add_timestamp(
+        self, feature_series: pd.Series, idx: int | None = None
+    ) -> pd.Series:
+        """Add to feature_series "time" keyword
+        For Bids specify with fs_features, for real time analysis with current time stamp
+        """
+
+    def load_model(self, model_name: _PathLike):
+        """Load sklearn model, that utilizes predict"""
+        with open(model_name, "rb") as fid:
+            self.model = cPickle.load(fid)
+
+    def _set_run(self):
+        """Initialize preprocessing, and feature estimation modules"""
+        nm_test_settings.test_settings(self.settings, self.nm_channels)
+
+        (CH_NAMES_USED, _, FEATURE_IDX, _) = self._get_ch_info(self.nm_channels)
+
+        self.features = self._set_features(
+            self.settings,
+            CH_NAMES_USED,
+            self.sfreq,
+            self.verbose,
+        )
+        self.rereference, self.nm_channels = self._set_rereference(
+            self.settings, self.nm_channels
+        )
+
+        self.resample = self._set_resampling(
+            settings=self.settings,
+            sfreq_old=self.sfreq,
+        )
+
+        self.notch_filter = self._set_notch_filter(
+            settings=self.settings,
+            sfreq=self.sfreq
+            if self.settings["methods"]["raw_resampling"] is False
+            else self.settings["raw_resampling_settings"]["resample_freq_hz"],
+        )
+
+        if self.coords:
+            self.projection = self._get_projection(
+                self.settings, self.nm_channels
+            )
+
+        self.run_analysis = nm_run_analysis.Run(
+            features=self.features,
+            settings=self.settings,
+            reference=self.rereference,
+            projection=self.projection,
+            resample=self.resample,
+            notch_filter=self.notch_filter,
+            verbose=self.verbose,
+            feature_idx=FEATURE_IDX,
+        )
+
+    def _set_features(
+        self,
+        settings: dict,
+        CH_NAMES_USED: list,
+        fs: int,
+        VERBOSE: bool,
+    ) -> nm_features.Features:
+        """initialize feature class from settings"""
+        return nm_features.Features(
+            s=settings, ch_names=CH_NAMES_USED, fs=fs, verbose=VERBOSE
+        )
+
+    def set_fs(self, fs: int | float) -> None:
+        self.sfreq = int(fs)
+
+    def _set_rereference(
+        self, settings: dict, nm_channels: pd.DataFrame
+    ) -> tuple[nm_rereference.RT_rereference | None, pd.DataFrame | None]:
+        """Initialize nm_rereference and update nm_channels
+        nm_channels are updated if no rereferencing is specified
+
+        Parameters
+        ----------
+        settings : dict
+            [description]
+        nm_channels : pd.DataFrame
+            [description]
+
+        Returns
+        -------
+        Tuple
+            nm_rereference object, updated nm_channels DataFrame
+        """
+        if settings["methods"]["re_referencing"] is True:
+            rereference = nm_rereference.RT_rereference(nm_channels)
+        else:
+            rereference = None
+            # reset nm_channels from default values
+            nm_channels["rereference"] = None
+            nm_channels["new_name"] = nm_channels["name"]
+        return rereference, nm_channels
+
+    def _set_resampling(
+        self, settings: dict, sfreq_old: int | float
+    ) -> nm_resample.Resample | None:
+        """Initialize Resampling
+
+        Parameters
+        ----------
+        settings : dict
+        fs : int
+
+        Returns
+        -------
+        nm_resample.Resample
+        """
+        if settings["methods"]["raw_resampling"]:
+            resample = nm_resample.Resample(
+                sfreq_old=sfreq_old,
+                sfreq_new=settings["raw_resampling_settings"][
+                    "resample_freq_hz"
+                ],
+            )
+        else:
+            resample = None
+        return resample
+
+    def _set_notch_filter(
+        self,
+        settings: dict,
+        sfreq: int | float,
+        line_noise: int | float = 50,
+    ) -> nm_filter.NotchFilter | None:
+        if settings["methods"]["notch_filter"]:
+            return nm_filter.NotchFilter(sfreq=sfreq, line_noise=line_noise)
+        return None
+
+    def set_linenoise(self, line_noise: int) -> None:
+        self.line_noise = line_noise
+
+    @staticmethod
+    def get_grids(settings: dict, PATH_GRIDS: str, GRID_TYPE: GRIDS) -> tuple:
+        """Read settings specified grids
+
+        Parameters
+        ----------
+        settings : dict
+        PATH_GRIDS : str
+        GRID_TYPE : GRIDS
+
+        Returns
+        -------
+        Tuple
+            grid_cortex, grid_subcortex,
+            might be None if not specified in settings
+        """
+
+        if settings["methods"]["project_cortex"] is True:
+            grid_cortex = nm_IO.read_grid(PATH_GRIDS, GRID_TYPE.CORTEX)
+        else:
+            grid_cortex = None
+        if settings["methods"]["project_subcortex"] is True:
+            grid_subcortex = nm_IO.read_grid(PATH_GRIDS, GRID_TYPE.SUBCORTEX)
+        else:
+            grid_subcortex = None
+        return grid_cortex, grid_subcortex
+
+    def _get_projection(
+        self, settings: dict, nm_channels: pd.DataFrame
+    ) -> nm_projection.Projection | None:
+        """Return projection of used coordinated and grids"""
+
+        if any(
+            (
+                settings["methods"]["project_cortex"],
+                settings["methods"]["project_subcortex"],
+            )
+        ):
+            grid_cortex, grid_subcortex = self.get_grids(
+                self.settings, self.path_grids, GRIDS
+            )
+            projection = nm_projection.Projection(
+                settings=settings,
+                grid_cortex=grid_cortex,
+                grid_subcortex=grid_subcortex,
+                coords=self.coords,
+                nm_channels=nm_channels,
+                plot_projection=False,
+            )
+            return projection
+        return
+
+    @staticmethod
+    def _get_ch_info(nm_channels: pd.DataFrame):
+        """Get used feature and label info from nm_channels"""
+
+        CH_NAMES_USED = nm_channels[nm_channels["used"] == 1][
+            "new_name"
+        ].tolist()
+        CH_TYPES_USED = nm_channels[nm_channels["used"] == 1]["type"].tolist()
+
+        # used channels for feature estimation
+        FEATURE_IDX = np.where(nm_channels["used"] & ~nm_channels["target"])[
+            0
+        ].tolist()
+
+        # If multiple targets exist, select only the first
+        LABEL_IDX = np.where(nm_channels["target"] == 1)[0]
+
+        return CH_NAMES_USED, CH_TYPES_USED, FEATURE_IDX, LABEL_IDX
+
+    @staticmethod
+    def _get_settings(settings: dict | _PathLike) -> dict:
+        if isinstance(settings, dict):
+            return settings
+        return nm_IO.read_settings(str(settings))
+
+    @staticmethod
+    def _get_nm_channels(
+        nm_channels: pd.DataFrame | _PathLike,
+    ) -> pd.DataFrame:
+        """Read nm_channels from path or specify via BIDS arguments.
+        Nexessary parameters are then
+        ch_names (list),
+        ch_types (list),
+        bads (list)
+        used_types (list)
+        target_keywords (list)
+        reference Union[list, str]
+        """
+        if isinstance(nm_channels, pd.DataFrame):
+            return nm_channels
+        if nm_channels:
+            if not os.path.isfile(nm_channels):
+                raise ValueError(
+                    "PATH_NM_CHANNELS is not a valid file. Got: "
+                    f"{nm_channels}"
+                )
+            return pd.read_csv(nm_channels)
+
+    @staticmethod
+    def _get_sess_lat(coords: dict) -> bool:
+        if len(coords["cortex_left"]["positions"]) == 0:
+            return True
+        elif len(coords["cortex_right"]["positions"]) == 0:
+            return False
+        raise ValueError(
+            "Either cortex_left or cortex_right positions must" " be provided."
+        )
+
+    def save_sidecar(self, folder_name: str):
+        """Save sidecar incuding fs, line_noise, coords, sess_right to
+        PATH_OUT and subfolder 'folder_name'"""
+
+        sidecar = {
+            "original_fs": self.sfreq,
+            "fs": self.run_analysis.fs,
+            "ch_names": self.features.ch_names,
+            "sess_right": self.sess_right,
+        }
+        if self.projection:
+            sidecar["coords"] = (self.projection.coords,)
+            if self.settings["methods"]["project_cortex"]:
+                sidecar["grid_cortex"] = self.projection.grid_cortex
+                sidecar[
+                    "proj_matrix_cortex"
+                ] = self.projection.proj_matrix_cortex
+            if self.settings["methods"]["project_subcortex"]:
+                sidecar["grid_subcortex"] = self.projection.grid_subcortex
+                sidecar[
+                    "proj_matrix_subcortex"
+                ] = self.projection.proj_matrix_subcortex
+
+        nm_IO.save_sidecar(sidecar, self.path_out, folder_name)
+
+    def save_settings(self, folder_name: _PathLike):
+        nm_IO.save_settings(self.settings, self.path_out, folder_name)
+
+    def save_nm_channels(self, folder_name: _PathLike):
+        nm_IO.save_nmchannels(self.nm_channels, self.path_out, folder_name)
+
+    def save_features(self, folder_name: _PathLike, feature_arr: pd.DataFrame):
+        nm_IO.save_features(feature_arr, self.path_out, folder_name)
+
+    def save_after_stream(
+        self, folder_name: _PathLike, feature_arr: pd.DataFrame | None = None
+    ) -> None:
+        """Save features, settings, nm_channels and sidecar after run"""
+
+        # create derivate folder_name output folder if doesn't exist
+        if os.path.exists(os.path.join(self.path_out, folder_name)) is False:
+            os.makedirs(os.path.join(self.path_out, folder_name))
+
+        self.save_sidecar(folder_name)
+
+        if feature_arr is not None:
+            self.save_features(folder_name, feature_arr)
+
+        self.save_settings(folder_name)
+
+        self.save_nm_channels(folder_name)
+
+    def plot_cortical_projection(self):
+        """plot projection of cortical grid electrodes on cortex"""
+
+        if hasattr(self, "features") is False:
+            self._set_run()
+
+        ecog_strip = None
+        if self.projection is not None:
+            ecog_strip = self.projection.ecog_strip
+
+        grid_cortex = None
+        if self.projection is not None:
+            grid_cortex = self.projection.grid_cortex
+
+        sess_right = None
+        if self.projection is not None:
+            sess_right = self.projection.sess_right
+
+        nmplotter = nm_plots.NM_Plot(
+            ecog_strip=ecog_strip,
+            grid_cortex=grid_cortex,
+            sess_right=sess_right,
+        )
+        nmplotter.plot_cortex(set_clim=False)
+
+    def _set_coords(
+        self, coord_names: list | None, coord_list: list | None
+    ) -> dict:
+        if not any(
+            (
+                self.settings["methods"]["project_cortex"],
+                self.settings["methods"]["project_subcortex"],
+            )
+        ):
+            return {}
+
+        if any((coord_list is None, coord_names is None)):
+            raise ValueError(
+                "No coordinates could be loaded. Please provide coord_list and"
+                f" coord_names. Got: {coord_list=}, {coord_names=}."
+            )
+
+        return self._add_coordinates(
+            coord_names=coord_names, coord_list=coord_list
+        )
+
+    @staticmethod
+    def _add_coordinates(coord_names: list, coord_list: list) -> dict:
+        """set coordinate information to settings from RawArray
+        The set coordinate positions are set as lists,
+        since np.arrays cannot be saved in json
+        Parameters
+        ----------
+        raw_arr : mne.io.RawArray
+        PATH_GRIDS : string, optional
+            absolute path to grid_cortex.tsv and grid_subcortex.tsv, by default: None
+        """
+
+        def left_coord(val: int | float, coord_region: str) -> bool:
+            if coord_region.split("_")[1] == "left":
+                return val < 0
+            return val > 0
+
+        coords = {}
+
+        for coord_region in [
+            coord_loc + "_" + lat
+            for coord_loc in ["cortex", "subcortex"]
+            for lat in ["left", "right"]
+        ]:
+
+            coords[coord_region] = {}
+
+            ch_type = (
+                "ECOG" if "cortex" == coord_region.split("_")[0] else "LFP"
+            )
+
+            coords[coord_region]["ch_names"] = [
+                coord_name
+                for coord_name, ch in zip(coord_names, coord_list)
+                if left_coord(ch[0], coord_region) and (ch_type in coord_name)
+            ]
+
+            # multiply by 1000 to get m instead of mm
+            coords[coord_region]["positions"] = 1000 * np.array(
+                [
+                    coord
+                    for coord, coord_name in zip(coord_list, coord_names)
+                    if left_coord(coord[0], coord_region)
+                    and (ch_type in coord_name)
+                ]
+            )
+
+        return coords
diff --git a/py_neuromodulation/nm_stream_offline.py b/py_neuromodulation/nm_stream_offline.py
new file mode 100644
index 00000000..03854c96
--- /dev/null
+++ b/py_neuromodulation/nm_stream_offline.py
@@ -0,0 +1,176 @@
+import os
+import timeit
+
+import mne
+import mne_bids
+import numpy as np
+import pandas as pd
+
+from py_neuromodulation import nm_generator, nm_IO, nm_stream_abc
+
+_PathLike = str | os.PathLike
+
+
+class _OfflineStream(nm_stream_abc.PNStream):
+    def _run_offline(
+        self,
+        data: np.ndarray,
+    ) -> None:
+        self._set_run()
+
+        generator = nm_generator.ieeg_raw_generator(
+            data, self.settings, self.sfreq
+        )
+        features = []
+        first_sample = True
+        while True:
+            data_batch = next(generator, None)
+            if data_batch is None:
+                break
+            feature_series = self.run_analysis.process_data(data_batch)
+
+            # Measuring timing
+            # number_repeat = 100
+            # val = timeit.timeit(
+            #    lambda: self.run_analysis.process_data(data),
+            #    number=number_repeat
+            # ) / number_repeat
+
+            feature_series = self._add_timestamp(feature_series, first_sample)
+            features.append(feature_series)
+
+            if self.model is not None:
+                prediction = self.model.predict(feature_series)
+
+            if first_sample:
+                first_sample = False
+
+        feature_df = pd.DataFrame(features)
+        feature_df = self._add_labels(features=feature_df, data=data)
+
+        self.save_after_stream(self.folder_name, feature_df)
+
+    def _handle_data(self, data: np.ndarray | pd.DataFrame) -> np.ndarray:
+        names_expected = self.nm_channels["name"].to_list()
+
+        if isinstance(data, np.ndarray):
+            if not len(names_expected) == data.shape[0]:
+                raise ValueError(
+                    "If data is passed as an array, the first dimension must"
+                    " match the number of channel names in `nm_channels`. Got:"
+                    f" Data columns: {data.shape[0]}, nm_channels.name:"
+                    f" {len(names_expected)}."
+                )
+            return data
+        names_data = data.columns.to_list()
+        if not (
+            len(names_expected) == len(names_data)
+            and sorted(names_expected) == sorted(names_data)
+        ):
+            raise ValueError(
+                "If data is passed as a DataFrame, the"
+                "columns must match the channel names in `nm_channels`. Got:"
+                f"Data columns: {names_data}, nm_channels.name: {names_data}."
+            )
+        return data.to_numpy()
+
+    def _add_timestamp(
+        self, feature_series: pd.Series, first_sample: bool
+    ) -> pd.Series:
+        """time stamp is added in ms
+        Due to normalization run_analysis needs to keep track of the counted samples
+        Those are accessed here for time conversion"""
+
+        if first_sample:
+            feature_series["time"] = self.run_analysis.offset
+        else:
+            # sampling frequency is taken from run_analysis, since resampling
+            # might change it
+            feature_series["time"] = (
+                self.run_analysis.cnt_samples * 1000 / self.run_analysis.fs
+            )
+
+        if self.verbose:
+            print(
+                str(np.round(feature_series["time"] / 1000, 2))
+                + " seconds of data processed"
+            )
+
+        return feature_series
+
+    def _add_labels(
+        self, features: pd.DataFrame, data: np.ndarray
+    ) -> pd.DataFrame:
+        """Add resampled labels to features if there are target channels."""
+        if self.nm_channels.target.sum() > 0:
+            features = nm_IO.add_labels(
+                df_=features,
+                settings=self.settings,
+                nm_channels=self.nm_channels,
+                raw_arr_data=data,
+                fs=self.sfreq,
+            )
+        return features
+
+
+class Stream(_OfflineStream):
+    def run(
+        self,
+        data: np.ndarray | pd.DataFrame,
+        sfreq: int | float,
+        coord_names: list | None = None,
+        coord_list: list | None = None,
+    ) -> None:
+        """BIDS specific fun function. Does not need to run in parallel."""
+        self.sfreq = int(sfreq)
+
+        data = self._handle_data(data)
+        self.coords = self._set_coords(
+            coord_names=coord_names, coord_list=coord_list
+        )
+        self._set_run()
+
+        self._run_offline(data)
+
+
+class BidsStream(_OfflineStream):
+    def run(
+        self,
+        filepath: _PathLike | mne_bids.BIDSPath,
+        datatype: str,
+        bids_root: _PathLike | None = None,
+    ) -> None:
+        """Run stream with BIDS file."""
+
+        self.folder_name = os.path.basename(filepath)
+
+        raw, data, sfreq, _ = nm_IO.read_BIDS_data(
+            PATH_RUN=filepath,
+            BIDS_PATH=bids_root,
+            datatype=datatype,
+        )
+        coord_list, coord_names = _get_coord_list(raw)
+        self.coords = self._set_coords(
+            coord_names=coord_names, coord_list=coord_list
+        )
+        self.sfreq = int(sfreq)
+
+        self._run_offline(data)
+
+
+def _get_coord_list(
+    raw: mne.io.BaseRaw,
+) -> tuple[list, list] | tuple[None, None]:
+    montage = raw.get_montage()
+    if montage is not None:
+        coord_list = np.array(
+            list(dict(montage.get_positions()["ch_pos"]).values())
+        ).tolist()
+        coord_names = np.array(
+            list(dict(montage.get_positions()["ch_pos"]).keys())
+        ).tolist()
+    else:
+        coord_list = None
+        coord_names = None
+
+    return coord_list, coord_names
diff --git a/pyproject.toml b/pyproject.toml
index 07de284a..a58c9c46 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,3 +1,6 @@
 [build-system]
 requires = ["setuptools", "wheel"]
-build-backend = "setuptools.build_meta"
\ No newline at end of file
+build-backend = "setuptools.build_meta"
+
+[tool.black]
+line-length = 80
diff --git a/scripts/offline_stream_example.py b/scripts/offline_stream_example.py
new file mode 100644
index 00000000..596a23ab
--- /dev/null
+++ b/scripts/offline_stream_example.py
@@ -0,0 +1,169 @@
+import os
+from multiprocessing import Pool
+import py_neuromodulation as nm
+import xgboost
+import numpy as np
+from sklearn.preprocessing import OneHotEncoder
+from sklearn.preprocessing import LabelEncoder
+
+import py_neuromodulation as nm
+from py_neuromodulation import (
+    nm_analysis,
+    nm_decode,
+    nm_define_nmchannels,
+    nm_IO,
+)
+from sklearn import linear_model, metrics, model_selection
+from skopt import space as skopt_space
+
+
+def set_settings(settings: dict):
+    for method in list(settings["methods"].keys()):
+        settings["methods"][method] = False
+
+    settings["methods"]["fft"] = True
+    settings["methods"]["fooof"] = True
+    settings["methods"]["return_raw"] = True
+    settings["methods"]["raw_hjorth"] = True
+    settings["methods"]["re_referencing"] = False
+    settings["methods"]["raw_normalization"] = False
+    settings["methods"]["feature_normalization"] = True
+
+    settings["fooof"]["periodic"]["center_frequency"] = False
+    settings["fooof"]["periodic"]["band_width"] = False
+    settings["fooof"]["periodic"]["height_over_ap"] = False
+
+    settings["methods"]["sharpwave_analysis"] = True
+
+    for key in list(
+        settings["sharpwave_analysis_settings"]["sharpwave_features"].keys()
+    ):
+        settings["sharpwave_analysis_settings"]["sharpwave_features"][
+            key
+        ] = True
+    settings["sharpwave_analysis_settings"]["sharpwave_features"][
+        "peak_left"
+    ] = False
+    settings["sharpwave_analysis_settings"]["sharpwave_features"][
+        "peak_right"
+    ] = False
+    settings["sharpwave_analysis_settings"][
+        "apply_estimator_between_peaks_and_troughs"
+    ] = True
+    settings["sharpwave_analysis_settings"]["filter_low_cutoff"] = 5
+    settings["sharpwave_analysis_settings"]["filter_high_cutoff"] = 40
+
+    settings["sharpwave_analysis_settings"]["estimator"]["max"].append("trough")
+    settings["sharpwave_analysis_settings"]["estimator"]["mean"].append("width")
+    settings["sharpwave_analysis_settings"]["estimator"]["max"].append(
+        "decay_time"
+    )
+    settings["sharpwave_analysis_settings"]["estimator"]["max"].append(
+        "rise_time"
+    )
+    settings["sharpwave_analysis_settings"]["estimator"]["max"].append(
+        "rise_steepness"
+    )
+    settings["sharpwave_analysis_settings"]["estimator"]["max"].append(
+        "decay_steepness"
+    )
+    settings["sharpwave_analysis_settings"]["estimator"]["mean"].append(
+        "slope_ratio"
+    )
+
+    return settings
+
+
+def preprocess_trd_data(dat):
+
+    labels = dat["labels"][~np.array(dat["bad"], dtype="bool")]
+    data = np.swapaxes(np.swapaxes(dat["data"], 0, 2), 1, 2)
+    data = data[~np.array(dat["bad"], dtype="bool"), :, :]
+    NUM_CH = data.shape[1]
+
+    label_encoder = LabelEncoder()
+    integer_encoded = label_encoder.fit_transform(labels)
+
+    enc = OneHotEncoder(handle_unknown="ignore", sparse=False)
+    label_arr = enc.fit_transform(integer_encoded.reshape(-1, 1))
+
+    label_arr_concat = np.concatenate(
+        (np.expand_dims(integer_encoded + 1, axis=1), label_arr), axis=1
+    )  # integer encoded + 1, since REST will be set to zero
+
+    label_arr_epochs = np.zeros([data.shape[0], 4, data.shape[2]])
+    label_arr_epochs_names = list(label_encoder.classes_)
+    label_arr_epochs_names.insert(0, "ALL")
+
+    integer_encoded_names = label_arr_epochs_names.copy()
+    integer_encoded_names[0] = "REST"
+
+    # label_arr_epochs columns: ALL, NTR, PLS, UNPLS
+    # integer_encoded_names: 0 - REST, 1 - NTR, 2 - PLS, 3 - UNPLS
+
+    arr_insert = np.repeat(label_arr_concat[:, :, np.newaxis], 1000, axis=2)
+    label_arr_epochs[:, :, 3500:4500] = arr_insert
+
+    data_comb = np.concatenate((data, label_arr_epochs), axis=1)
+    data_stream = np.concatenate(data_comb, axis=1)
+
+    ch_names = list(dat["ch_names"])
+    ch_names = ch_names + label_arr_epochs_names
+
+    ch_types = ["seeg" for _ in range(NUM_CH)]
+    ch_types = ch_types + ["misc" for _ in range(len(label_arr_epochs_names))]
+
+    return data_stream, ch_names, ch_types, dat["fsample"]
+
+
+def run_patient_offline_stream(f):
+
+    file_name = os.path.basename(f)[: -len("_edit.mat")]
+    dat = nm_IO.loadmat(os.path.join(PATH_DATA, f))["D"]
+
+    data_stream, ch_names, ch_types, fs = preprocess_trd_data(dat)
+
+    # add nm_channels
+
+    nm_channels = nm_define_nmchannels.set_channels(
+        ch_names=ch_names,
+        ch_types=ch_types,
+        reference=None,
+        bads=[],
+        new_names="default",
+        used_types=("seeg",),
+        target_keywords=None,
+    )
+    stream.nm_channels.loc[
+        stream.nm_channels.query('type == "misc"').index, "target"
+    ] = 1
+
+    PATH_OUT = r"C:\Users\ICN_admin\Documents\TRD Analysis\features_epochs_OUT"
+    stream = nm.Stream(
+        settings=None,
+        nm_channels=nm_channels,
+        path_out=PATH_OUT,
+        path_grids=None,
+        verbose=True,
+    )
+    stream.settings = set_settings(stream.settings)
+    stream.run(
+        data=data_stream,
+        sfreq=fs,
+    )
+
+
+PATH_DATA = r"C:\Users\ICN_admin\Documents\TRD Analysis"
+
+
+def main():
+
+    files = [f for f in os.listdir(PATH_DATA) if "_edit" in f]
+    for f in files:
+        run_patient_offline_stream(f)
+    # pool = Pool(processes=len(files))
+    # pool.map(run_patient_GenericStream, files)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/py_nm_batch_trd_data.py b/scripts/py_nm_batch_trd_data.py
deleted file mode 100644
index 6ebc3986..00000000
--- a/scripts/py_nm_batch_trd_data.py
+++ /dev/null
@@ -1,168 +0,0 @@
-import os
-
-import numpy as np
-import scipy.io as spio
-from sklearn.preprocessing import OneHotEncoder
-from sklearn.preprocessing import LabelEncoder
-from multiprocessing import Pool
-import pickle
-
-from yaml import load
-
-
-def loadmat(filename):
-    """
-    this function should be called instead of direct spio.loadmat
-    as it cures the problem of not properly recovering python dictionaries
-    from mat files. It calls the function check keys to cure all entries
-    which are still mat-objects
-    """
-    data = spio.loadmat(filename, struct_as_record=False, squeeze_me=True)
-    return _check_keys(data)
-
-
-def _check_keys(dict):
-    """
-    checks if entries in dictionary are mat-objects. If yes
-    todict is called to change them to nested dictionaries
-    """
-    for key in dict:
-        if isinstance(dict[key], spio.matlab.mio5_params.mat_struct):
-            dict[key] = _todict(dict[key])
-    return dict
-
-
-def _todict(matobj):
-    """
-    A recursive function which constructs from matobjects nested dictionaries
-    """
-    dict = {}
-    for strg in matobj._fieldnames:
-        elem = matobj.__dict__[strg]
-        if isinstance(elem, spio.matlab.mio5_params.mat_struct):
-            dict[strg] = _todict(elem)
-        else:
-            dict[strg] = elem
-    return dict
-
-
-from py_neuromodulation import nm_EpochStream
-
-
-def run_patient(f):
-    epoch_stream = nm_EpochStream.EpochStream()
-    # epoch_stream.read_epoch_data(
-    #    os.path.join("scripts", "preproc_ecog_v0.3_Jan22_rerefSegments_data.npy")
-    # )
-    file_name = os.path.basename(f)[: -len("_edit.mat")]
-
-    dat = loadmat(os.path.join(PATH_DATA, f))["D"]
-    labels = dat["labels"][~np.array(dat["bad"], dtype="bool")]
-    # bring data into shape (samples, channels, time)
-    epoch_stream.data = np.swapaxes(np.swapaxes(dat["data"], 0, 2), 1, 2)
-
-    # use only good trials
-    epoch_stream.data = epoch_stream.data[~np.array(dat["bad"], dtype="bool"), :, :]
-
-    epoch_stream.set_fs(dat["fsample"])
-    epoch_stream.set_linenoise(50)
-    NUM_CH = epoch_stream.data.shape[1]
-
-    epoch_stream.nm_channels = epoch_stream._get_nm_channels(
-        PATH_NM_CHANNELS=None,
-        ch_names=dat["ch_names"],
-        ch_types=["lfp" for _ in range(NUM_CH)],
-        bads=[],
-        used_types=["lfp" for _ in range(NUM_CH)],
-        target_keywords=[],
-        reference=None,
-    )
-
-    for method in list(epoch_stream.settings["methods"].keys()):
-        epoch_stream.settings["methods"][method] = False
-
-    epoch_stream.settings["methods"]["fft"] = True
-    epoch_stream.settings["methods"]["sharpwave_analysis"] = True
-
-    for key in list(
-        epoch_stream.settings["sharpwave_analysis_settings"][
-            "sharpwave_features"
-        ].keys()
-    ):
-        epoch_stream.settings["sharpwave_analysis_settings"]["sharpwave_features"][
-            key
-        ] = False
-    epoch_stream.settings["sharpwave_analysis_settings"][
-        "apply_estimator_between_peaks_and_troughs"
-    ] = True
-    epoch_stream.settings["sharpwave_analysis_settings"]["sharpwave_features"][
-        "prominence"
-    ] = True
-    epoch_stream.settings["sharpwave_analysis_settings"]["sharpwave_features"][
-        "sharpness"
-    ] = True
-    epoch_stream.settings["sharpwave_analysis_settings"]["filter_low_cutoff"] = 5
-    epoch_stream.settings["sharpwave_analysis_settings"]["filter_high_cutoff"] = 40
-
-    epoch_stream.run()
-
-    label_encoder = LabelEncoder()
-    integer_encoded = label_encoder.fit_transform(labels)
-
-    enc = OneHotEncoder(handle_unknown="ignore", sparse=False)
-    label_arr = enc.fit_transform(integer_encoded.reshape(-1, 1))
-
-    dict_out = {}
-    dict_out["label_" + label_encoder.classes_[0]] = label_arr[:, 0]
-    dict_out["label_" + label_encoder.classes_[1]] = label_arr[:, 1]
-    dict_out["label_" + label_encoder.classes_[2]] = label_arr[:, 2]
-
-    dict_out["label"] = labels
-    dict_out["feature_names"] = list(epoch_stream.feature_arr.columns)
-
-    # add labels to nm_channels
-    for label_name in ["label_PLS", "label_UNPLS", "label_NTR"]:
-        epoch_stream.nm_channels = epoch_stream.nm_channels.append(
-            {
-                "name": label_name,
-                "rereference": None,
-                "used": 0,
-                "target": 1,
-                "type": "misc",
-                "status": "good",
-                "new_name": label_name,
-            },
-            ignore_index=True,
-        )
-
-    epoch_stream.PATH_OUT = os.path.join(PATH_DATA, "features")
-
-    dict_out["features"] = np.stack(
-        [
-            np.array(epoch_stream.feature_arr_list[idx])
-            for idx in range(len(epoch_stream.feature_arr_list))
-        ]
-    )
-
-    with open(os.path.join(PATH_DATA, "features", file_name, "dict_out.p"), "wb") as fp:
-        pickle.dump(dict_out, fp)
-
-    epoch_stream.save_after_stream(file_name, save_features=False)
-
-    print(f"finished processing {file_name}")
-
-
-PATH_DATA = r"C:\Users\ICN_admin\Documents\TRD Analysis"
-
-
-def main():
-
-    files = [f for f in os.listdir(PATH_DATA) if "_edit" in f]
-    # for f in files:
-    #    run_patient(f)
-    pool = Pool(processes=len(files))
-    pool.map(run_patient, files)
-
-
-if __name__ == "__main__":
-    main()
diff --git a/scripts/py_nm_batch_trd_analysis.py b/scripts/trd_analysis_epochs.py
similarity index 100%
rename from scripts/py_nm_batch_trd_analysis.py
rename to scripts/trd_analysis_epochs.py
diff --git a/scripts/trd_analysis_stream_old.py b/scripts/trd_analysis_stream_old.py
new file mode 100644
index 00000000..d4335d1f
--- /dev/null
+++ b/scripts/trd_analysis_stream_old.py
@@ -0,0 +1,61 @@
+import os
+from sklearn import linear_model, metrics, model_selection, discriminant_analysis
+from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
+from skopt import space as skopt_space
+import xgboost
+from py_neuromodulation import nm_GenericStream, nm_analysis, nm_decode
+
+
+def main():
+
+    PATH_OUT = r"C:\Users\ICN_admin\Documents\TRD Analysis\features_epochs_nonorm"
+    RUN_NAME = "effspm8_JUN_EMO"
+    feature_reader = nm_analysis.Feature_Reader(
+        feature_dir=PATH_OUT, feature_file=RUN_NAME, binarize_label=False
+    )
+
+    model = linear_model.LogisticRegression()
+    model = xgboost.XGBClassifier()
+    feature_reader.decoder = nm_decode.Decoder(
+        features=feature_reader.feature_arr,
+        label=feature_reader.label,
+        label_name=feature_reader.label_name,
+        used_chs=feature_reader.used_chs,
+        model=model,
+        eval_method=metrics.balanced_accuracy_score,
+        cv_method=model_selection.KFold(n_splits=3, shuffle=True),
+        # cv_method="NonShuffledTrainTestSplit",
+        get_movement_detection_rate=False,
+        min_consequent_count=2,
+        TRAIN_VAL_SPLIT=False,
+        RUN_BAY_OPT=False,
+        bay_opt_param_space=None,
+        use_nested_cv=False,
+        fs=feature_reader.settings["sampling_rate_features_hz"],
+        oversampling=False,
+        undersampling=True,
+    )
+
+    performances = feature_reader.run_ML_model(
+        estimate_channels=True,
+        estimate_gridpoints=False,
+        estimate_all_channels_combined=True,
+        save_results=True,
+    )
+
+    class_labels = ["rest", "ntr", "pls", "unpls"]
+    cm = confusion_matrix(
+        feature_reader.decoder.ch_ind_results["Cg25R01"]["y_test"][0],
+        feature_reader.decoder.ch_ind_results["Cg25R01"]["y_test_pr"][0],
+        normalize="true"
+    )
+
+    disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=class_labels)
+    disp.plot()
+
+    # feature_reader.decoder.ch_ind_results["Cg25R01"]["y_test_prfeature_reader"
+    print("Analysis finished")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/trd_analysis_stream_plotting.py b/scripts/trd_analysis_stream_plotting.py
new file mode 100644
index 00000000..821b6c86
--- /dev/null
+++ b/scripts/trd_analysis_stream_plotting.py
@@ -0,0 +1,306 @@
+from matplotlib import pyplot as plt
+import seaborn as sb
+from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
+import pandas as pd
+import numpy as np
+import pickle
+from scipy import stats
+import os
+import sys
+
+from py_neuromodulation import nm_plots, nm_analysis
+
+sys.path.append(r"C:\Users\ICN_admin\Documents\icn\icn_stats")
+import icn_stats
+
+
+def plot_df_subjects(
+    df, PATH_OUT, x_col="sub", y_col="performance_test", hue="all combined"
+):
+    alpha_box = 0.4
+    plt.figure(figsize=(5, 3), dpi=300)
+    sb.boxplot(
+        x=x_col,
+        y=y_col,
+        hue=hue,
+        data=df,
+        palette="viridis",
+        showmeans=False,
+        boxprops=dict(alpha=alpha_box),
+        showcaps=True,
+        showbox=True,
+        showfliers=False,
+        notch=False,
+        whiskerprops={"linewidth": 2, "zorder": 10, "alpha": alpha_box},
+        capprops={"alpha": alpha_box},
+        medianprops=dict(linestyle="-.", linewidth=5, color="gray", alpha=alpha_box),
+    )
+
+    ax = sb.stripplot(
+        x=x_col,
+        y=y_col,
+        hue=hue,
+        data=df,
+        palette="viridis",
+        dodge=True,
+        s=5,
+    )
+
+    # When creating the legend, only use the first two elements
+    # to effectively remove the last two.
+    handles, labels = ax.get_legend_handles_labels()
+    l = plt.legend(
+        handles[0:2], labels[0:2], bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0
+    )
+    plt.title("subject specific channel performances")
+    plt.ylabel(y_col)
+    plt.xticks(rotation=90)
+    plt.savefig(
+        os.path.join(PATH_OUT, "all_sub_performances.png"),
+        bbox_inches="tight",
+    )
+
+
+def main():
+
+    PATH_OUT = r"C:\Users\ICN_admin\Documents\TRD Analysis\features_epochs_nonorm"
+    feature_reader = nm_analysis.Feature_Reader(
+        feature_dir=PATH_OUT, feature_file="effspm8_JUN_EMO", binarize_label=False
+    )
+
+    with open("dict_res_out_PLS_UNPLS_MI.p", "rb") as handle:
+        d = pickle.load(handle)
+
+    class_labels = ["rest", "ntr", "pls", "unpls"]
+
+    depression_scores = {
+        "BDI24_am": {
+            "JUN": 37,
+            "KOR": 36,
+            "MIC": 32,
+            "NIL": 37,
+            "OHL": 1,
+            "SCH": 31,
+            "THI": 12,
+            "WES": 52,
+        },
+        "HAMD24_24": {
+            "JUN": 33,
+            "KOR": 27,
+            "MIC": 22,
+            "NIL": 29,
+            "OHL": 9,
+            "SCH": 16,
+            "THI": 7,
+            "WES": 33,
+        },
+    }
+
+    df = pd.DataFrame()
+    for run in d.keys():
+        for ch in d[run].keys():
+
+            if ch.startswith("all_ch_"):
+                ALL_COMB = True
+            else:
+                ALL_COMB = False
+            df = df.append(
+                {
+                    "performance_test": d[run][ch]["performances"],
+                    "sub": run[len("effspm8_") : len("effspm8_") + 3],
+                    "ch": ch,
+                    "all_comb": ALL_COMB,
+                    "BDI": depression_scores["BDI24_am"][
+                        run[len("effspm8_") : len("effspm8_") + 3]
+                    ],
+                    "HAMD": depression_scores["HAMD24_24"][
+                        run[len("effspm8_") : len("effspm8_") + 3]
+                    ],
+                },
+                ignore_index=True,
+            )
+
+    print("plot mean performances")
+
+    # plot depression correlation
+    metric = "BDI"
+    df_plt = df.query("all_comb == False").groupby("sub").apply("mean")
+    df_plt = df.query("all_comb == True")
+    rho, p = icn_stats.permutationTestSpearmansRho(
+        df_plt["performance_test"],
+        df_plt[metric],
+        False,
+        "R^2",
+        5000,
+    )
+    sb.regplot(x="performance_test", y=metric, data=df_plt)
+    plt.title(f"Accuracy~{metric} p={np.round(p, 2)} rho={np.round(rho, 2)}")
+
+    cm = []
+    x_epoch = []
+    y_epoch = []
+    y_pr = []
+    y_te = []
+
+    mi_series = []
+
+    for run in d.keys():
+        for ch in d[run].keys():
+            cm.append(d[run][ch]["cm"])
+            if ch.startswith("all") is False:
+                x_epoch.append(np.squeeze(d[run][ch]["X_epoch"].mean(axis=0)))
+                y_epoch.append(d[run][ch]["y_epoch"].mean(axis=0))
+            else:
+                y_pr.append(d[run][ch]["y_test_pr"])
+                y_te.append(d[run][ch]["y_test"])
+        mi_series.append(d[run][ch]["mi_series_sorted"])
+
+    # time plot
+    y_pr_ = np.concatenate(y_pr)
+    y_te_ = np.concatenate(y_te)
+    y_te_[np.where(y_te_ != 3)[0]] = 0
+    X_, y_ = feature_reader.get_epochs(
+        np.expand_dims(y_pr_, axis=(1, 2)), y_te_, epoch_len=4, sfreq=10, threshold=0.1
+    )
+    y_unpls = np.squeeze(X_).mean(axis=0)
+
+    # 0 - REST, 1 - NTR, 2 - PLS, 3 - UNPLS
+    y_te_ = np.concatenate(y_te)
+    y_te_[np.where(y_te_ != 2)[0]] = 0
+    X_, y_ = feature_reader.get_epochs(
+        np.expand_dims(y_pr_, axis=(1, 2)), y_te_, epoch_len=4, sfreq=10, threshold=0.1
+    )
+    y_pls = np.squeeze(X_).mean(axis=0)
+
+    y_te_ = np.concatenate(y_te)
+    y_te_[np.where(y_te_ != 1)[0]] = 0
+    X_, y_ = feature_reader.get_epochs(
+        np.expand_dims(y_pr_, axis=(1, 2)), y_te_, epoch_len=4, sfreq=10, threshold=0.1
+    )
+    y_ntr = np.squeeze(X_).mean(axis=0)
+
+    y_te_ = np.concatenate(y_te)
+    y_te_[np.where(y_te_ == 0)[0]] = 5
+    y_te_[np.where(y_te_ != 5)[0]] = 0
+    X_, y_ = feature_reader.get_epochs(
+        np.expand_dims(y_pr_, axis=(1, 2)), y_te_, epoch_len=4, sfreq=10, threshold=0.1
+    )
+    y_rest = np.squeeze(X_).mean(axis=0)
+
+    plt.figure(figsize=(7, 4), dpi=300)
+    plt.plot(y_rest, linewidth=2, color="black", label="rest")
+    plt.plot(y_ntr, linewidth=2, color="gray", label="ntr")
+    plt.plot(y_pls, linewidth=2, color="blue", label="pls")
+    plt.plot(y_unpls, linewidth=2, color="green", label="unpls")
+    plt.xticks(np.arange(0, 40, 2), np.round(np.arange(-2, 2, 0.2), 2))
+    plt.xlabel("Time [s]")
+    plt.legend()
+    plt.ylabel("Prediction")
+    plt.title(
+        "XGBOOST mean predictions all channels\n 0 - REST, 1 - NTR, 2 - PLS, 3 - UNPLS"
+    )
+    plt.tight_layout()
+    plt.show()
+
+    print("plot")
+
+    # plot best
+    best_series = d["effspm8_SCH_EMO"]["Cg25L23"]["mi_series_sorted"]
+    best_series = best_series.reset_index()
+    best_series[best_series["index"].str.startswith("Cg25L23")].set_index(
+        "index"
+    ).plot.bar()
+    plt.ylabel("Mutual Information")
+    plt.title("Sorted Best subject ch Cg25L23 MI feature scores")
+    plt.tight_layout()
+
+    # plot mean all
+    df_mean_sorted = (
+        pd.concat(mi_series, axis=1).mean(axis=1).sort_values(ascending=False)
+    )
+    df_mean_sorted.iloc[:20].plot.bar()
+    plt.ylabel("Mutual Information")
+    plt.title("Sorted Mean subject MI feature scores")
+    plt.tight_layout()
+
+    disp = ConfusionMatrixDisplay(
+        confusion_matrix=np.array(cm).mean(axis=0), display_labels=class_labels
+    )
+    disp.plot()
+
+    # plot best confusion matrix
+    disp = ConfusionMatrixDisplay(
+        confusion_matrix=np.array(d["effspm8_SCH_EMO"]["Cg25L23"]["cm"]),
+        display_labels=class_labels,
+    )
+    disp.plot()
+
+    plot_df_subjects(
+        df, PATH_OUT, x_col="sub", y_col="performance_test", hue="all_comb"
+    )
+
+    # have a look at the mi rates scors: in total there are 23 features per channel
+    # take average series of all subjects and sort?
+
+    # plot best
+    x_best = stats.zscore(
+        np.nan_to_num(
+            np.squeeze(d["effspm8_SCH_EMO"]["Cg25L23"]["X_epoch"].mean(axis=0))
+        ),
+        axis=1,
+    )
+    y_best = d["effspm8_SCH_EMO"]["Cg25L23"]["y_epoch"]
+    plot_epoch(x_best.T, y_best, x_best.T)
+
+    # plot all
+    plot_epoch(x_epoch, y_epoch)
+
+    def plot_epoch(x_epoch, y_epoch, X_epoch_mean=None):
+        epoch_len = 4
+        feature_names = [f[8:] for f in feature_reader.feature_arr.columns[:23]]
+        label_name = "UNPLS"
+        X_epoch = np.array(x_epoch)
+        sfreq = 10
+        if X_epoch_mean is None:
+            X_epoch_mean = stats.zscore(
+                np.nan_to_num(np.nanmean(np.squeeze(X_epoch), axis=0)), axis=0
+            ).T
+        y_epoch = np.stack(np.array(y_epoch))
+        plt.figure(figsize=(6, 6))
+        plt.subplot(211)
+        plt.imshow(X_epoch_mean, aspect="auto")
+        plt.yticks(np.arange(0, len(feature_names), 1), feature_names)
+        plt.xticks(
+            np.arange(0, X_epoch.shape[1], 1),
+            np.round(np.arange(-epoch_len / 2, epoch_len / 2, 1 / sfreq), 2),
+            rotation=90,
+        )
+        plt.gca().invert_yaxis()
+        plt.xlabel("Time [s]")
+        str_title = "UNPLS aligned features mean all channels"
+        plt.title(str_title)
+
+        plt.subplot(212)
+        for i in range(y_epoch.shape[0]):
+            plt.plot(y_epoch[i, :], color="black", alpha=0.4)
+        plt.plot(
+            y_epoch.mean(axis=0),
+            color="black",
+            alpha=1,
+            linewidth=3.0,
+            label="mean target",
+        )
+        plt.legend()
+        plt.ylabel("target")
+        plt.title(label_name)
+        plt.xticks(
+            np.arange(0, X_epoch.shape[1], 1),
+            np.round(np.arange(-epoch_len / 2, epoch_len / 2, 1 / sfreq), 2),
+            rotation=90,
+        )
+        plt.xlabel("Time [s]")
+        plt.tight_layout()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/trd_run_features.py b/scripts/trd_run_features.py
new file mode 100644
index 00000000..0c956038
--- /dev/null
+++ b/scripts/trd_run_features.py
@@ -0,0 +1,254 @@
+import os
+from re import VERBOSE
+
+import numpy as np
+import scipy.io as spio
+from sklearn.preprocessing import OneHotEncoder
+from sklearn.preprocessing import LabelEncoder
+from multiprocessing import Pool
+import pickle
+
+from py_neuromodulation import nm_EpochStream, nm_GenericStream
+
+from yaml import load
+
+
+def loadmat(filename):
+    """
+    this function should be called instead of direct spio.loadmat
+    as it cures the problem of not properly recovering python dictionaries
+    from mat files. It calls the function check keys to cure all entries
+    which are still mat-objects
+    """
+    data = spio.loadmat(filename, struct_as_record=False, squeeze_me=True)
+    return _check_keys(data)
+
+
+def _check_keys(dict):
+    """
+    checks if entries in dictionary are mat-objects. If yes
+    todict is called to change them to nested dictionaries
+    """
+    for key in dict:
+        if isinstance(dict[key], spio.matlab.mio5_params.mat_struct):
+            dict[key] = _todict(dict[key])
+    return dict
+
+
+def _todict(matobj):
+    """
+    A recursive function which constructs from matobjects nested dictionaries
+    """
+    dict = {}
+    for strg in matobj._fieldnames:
+        elem = matobj.__dict__[strg]
+        if isinstance(elem, spio.matlab.mio5_params.mat_struct):
+            dict[strg] = _todict(elem)
+        else:
+            dict[strg] = elem
+    return dict
+
+
+def set_settings(settings: dict):
+    for method in list(settings["methods"].keys()):
+        settings["methods"][method] = False
+
+    settings["methods"]["fft"] = True
+    settings["methods"]["fooof"] = True
+    settings["methods"]["return_raw"] = True
+    settings["methods"]["raw_hjorth"] = True
+    settings["methods"]["re_referencing"] = False
+    settings["methods"]["raw_normalization"] = False
+    settings["methods"]["feature_normalization"] = True
+
+    settings["fooof"]["periodic"]["center_frequency"] = False
+    settings["fooof"]["periodic"]["band_width"] = False
+    settings["fooof"]["periodic"]["height_over_ap"] = False
+
+    settings["methods"]["sharpwave_analysis"] = True
+
+    for key in list(
+        settings["sharpwave_analysis_settings"]["sharpwave_features"].keys()
+    ):
+        settings["sharpwave_analysis_settings"]["sharpwave_features"][
+            key
+        ] = True
+    settings["sharpwave_analysis_settings"]["sharpwave_features"][
+        "peak_left"
+    ] = False
+    settings["sharpwave_analysis_settings"]["sharpwave_features"][
+        "peak_right"
+    ] = False
+    settings["sharpwave_analysis_settings"][
+        "apply_estimator_between_peaks_and_troughs"
+    ] = True
+    settings["sharpwave_analysis_settings"]["filter_low_cutoff"] = 5
+    settings["sharpwave_analysis_settings"]["filter_high_cutoff"] = 40
+    return settings
+
+
+def run_patient_GenericStream(f):
+
+    file_name = os.path.basename(f)[: -len("_edit.mat")]
+    dat = loadmat(os.path.join(PATH_DATA, f))["D"]
+    labels = dat["labels"][~np.array(dat["bad"], dtype="bool")]
+
+    data = np.swapaxes(np.swapaxes(dat["data"], 0, 2), 1, 2)
+    data = data[~np.array(dat["bad"], dtype="bool"), :, :]
+    NUM_CH = data.shape[1]
+
+    label_encoder = LabelEncoder()
+    integer_encoded = label_encoder.fit_transform(labels)
+
+    enc = OneHotEncoder(handle_unknown="ignore", sparse=False)
+    label_arr = enc.fit_transform(integer_encoded.reshape(-1, 1))
+
+    label_arr_concat = np.concatenate(
+        (np.expand_dims(integer_encoded + 1, axis=1), label_arr), axis=1
+    )  # integer encoded + 1, since REST will be set to zero
+
+    label_arr_epochs = np.zeros([data.shape[0], 4, data.shape[2]])
+    label_arr_epochs_names = list(label_encoder.classes_)
+    label_arr_epochs_names.insert(0, "ALL")
+
+    integer_encoded_names = label_arr_epochs_names.copy()
+    integer_encoded_names[0] = "REST"
+
+    # label_arr_epochs columns: ALL, NTR, PLS, UNPLS
+    # integer_encoded_names: 0 - REST, 1 - NTR, 2 - PLS, 3 - UNPLS
+
+    arr_insert = np.repeat(label_arr_concat[:, :, np.newaxis], 1000, axis=2)
+    label_arr_epochs[:, :, 3500:4500] = arr_insert
+
+    data_comb = np.concatenate((data, label_arr_epochs), axis=1)
+    data_stream = np.concatenate(data_comb, axis=1)
+
+    ch_names = list(dat["ch_names"])
+    ch_names = ch_names + label_arr_epochs_names
+
+    ch_types = ["seeg" for _ in range(NUM_CH)]
+    ch_types = ch_types + ["misc" for _ in range(len(label_arr_epochs_names))]
+
+    stream = nm_GenericStream.GenericStream(
+        READ_BIDS=False,
+        ch_names=ch_names,
+        ch_types=ch_types,
+        fs=dat["fsample"],
+        line_noise=50,
+        PATH_NM_CHANNELS=None,
+        bads=[],
+        reference=None,
+        PATH_OUT=os.path.join(PATH_DATA, "features_epochs_nonorm"),
+        VERBOSE=False,
+    )
+
+    stream.set_data(data_stream)
+    stream.settings = set_settings(stream.settings)
+    stream.nm_channels.loc[
+        stream.nm_channels.query('type == "misc"').index, "target"
+    ] = 1
+
+    stream.run()
+
+    stream.save_after_stream(file_name)
+
+
+def run_patient_EpochStream(f):
+    epoch_stream = nm_EpochStream.EpochStream(VERBOSE=False)
+    # epoch_stream.read_epoch_data(
+    #    os.path.join("scripts", "preproc_ecog_v0.3_Jan22_rerefSegments_data.npy")
+    # )
+    file_name = os.path.basename(f)[: -len("_edit.mat")]
+
+    dat = loadmat(os.path.join(PATH_DATA, f))["D"]
+    labels = dat["labels"][~np.array(dat["bad"], dtype="bool")]
+    # bring data into shape (samples, channels, time)
+    epoch_stream.data = np.swapaxes(np.swapaxes(dat["data"], 0, 2), 1, 2)
+
+    # use only good trials
+    epoch_stream.data = epoch_stream.data[
+        ~np.array(dat["bad"], dtype="bool"), :, :
+    ]
+
+    epoch_stream.set_fs(dat["fsample"])
+    epoch_stream.set_linenoise(50)
+    NUM_CH = epoch_stream.data.shape[1]
+
+    epoch_stream.nm_channels = epoch_stream._get_nm_channels(
+        PATH_NM_CHANNELS=None,
+        ch_names=dat["ch_names"],
+        ch_types=["lfp" for _ in range(NUM_CH)],
+        bads=[],
+        used_types=["lfp" for _ in range(NUM_CH)],
+        target_keywords=[],
+        reference=None,
+    )
+
+    epoch_stream.settings = set_settings(epoch_stream.settings)
+    epoch_stream.run()
+
+    label_encoder = LabelEncoder()
+    integer_encoded = label_encoder.fit_transform(labels)
+
+    enc = OneHotEncoder(handle_unknown="ignore", sparse=False)
+    label_arr = enc.fit_transform(integer_encoded.reshape(-1, 1))
+
+    dict_out = {}
+    dict_out["label_" + label_encoder.classes_[0]] = label_arr[:, 0]
+    dict_out["label_" + label_encoder.classes_[1]] = label_arr[:, 1]
+    dict_out["label_" + label_encoder.classes_[2]] = label_arr[:, 2]
+
+    dict_out["label"] = labels
+    dict_out["feature_names"] = list(epoch_stream.feature_arr.columns)
+
+    # add labels to nm_channels
+    for label_name in ["label_PLS", "label_UNPLS", "label_NTR"]:
+        epoch_stream.nm_channels = epoch_stream.nm_channels.append(
+            {
+                "name": label_name,
+                "rereference": None,
+                "used": 0,
+                "target": 1,
+                "type": "misc",
+                "status": "good",
+                "new_name": label_name,
+            },
+            ignore_index=True,
+        )
+
+    epoch_stream.path_out = os.path.join(PATH_DATA, "features_epochs_nonorm")
+
+    dict_out["features"] = np.stack(
+        [
+            np.array(epoch_stream.feature_arr_list[idx])
+            for idx in range(len(epoch_stream.feature_arr_list))
+        ]
+    )
+
+    with open(
+        os.path.join(
+            PATH_DATA, "features_epochs_nonorm", file_name, "dict_out.p"
+        ),
+        "wb",
+    ) as fp:
+        pickle.dump(dict_out, fp)
+
+    epoch_stream.save_after_stream(file_name, save_features=False)
+
+    print(f"finished processing {file_name}")
+
+
+PATH_DATA = r"C:\Users\ICN_admin\Documents\TRD Analysis"
+
+
+def main():
+
+    files = [f for f in os.listdir(PATH_DATA) if "_edit" in f]
+    # for f in files:
+    #    run_patient_GenericStream(f)
+    pool = Pool(processes=len(files))
+    pool.map(run_patient_GenericStream, files)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/setup.cfg b/setup.cfg
index d947c8ff..41145c94 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [pycodestyle]
-max-line-length = 120
+max-line-length = 80
 
 [metadata]
 description-file = README.md


Test patch:
diff --git a/tests/test_example.py b/tests/test_example.py
index a07ffa36..44a606f4 100644
--- a/tests/test_example.py
+++ b/tests/test_example.py
@@ -1,5 +1,8 @@
 from examples import example_BIDS
 
-# test is run by "python -m pytest tests/" from pyneuromodulation root dr 
+# test is run by "python -m pytest tests/" from pyneuromodulation root dr
 def test_example():
     example_BIDS.run_example_BIDS()
+
+
+test_example()