update app.py

.gitignore

@@ -0,0 +1,3 @@
+
+*.pyc
+*.zip

BrainPulse/README.md

@@ -0,0 +1,2 @@
+# BrainPulse
+ 

BrainPulse/dataset.py

@@ -0,0 +1,153 @@
+from mne import Epochs, pick_types, events_from_annotations, create_info, EpochsArray
+from mne.channels import make_standard_montage
+from mne.io import concatenate_raws, read_raw_edf, RawArray
+from mne.datasets import eegbci
+import scipy.io
+import numpy as np
+import os
+import mne
+
+
+
+def eegbci_data(tmin, tmax, subject, filter_range = None, run_list = None):
+
+    event_id = dict(ev=0)
+    runs = run_list  # open eyes vs closed eyes
+
+    # raw_fnames = eegbci.load_data(subject, runs, path='./datasets', update_path=True)
+    raw_fnames = eegbci.load_data(subject, runs, update_path=True, path='../mne_data')
+    raw = concatenate_raws([read_raw_edf(f, preload=True) for f in raw_fnames])
+    eegbci.standardize(raw)  # set channel names
+    montage = make_standard_montage('standard_1005')
+    raw.set_montage(montage)
+
+    # strip channel names of "." characters
+    raw.rename_channels(lambda x: x.strip('.'))
+    raw.set_eeg_reference(projection=True)
+    raw.apply_proj()
+# Apply band-pass filter
+    if filter_range != None:
+        raw.filter(filter_range[0], filter_range[1], fir_design='firwin', skip_by_annotation='edge')
+
+    events, _ = events_from_annotations(raw, event_id=dict(T0=0))
+
+    picks = pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False,
+                       exclude='bads')
+    # Read epochs
+    epochs = Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks,
+                    baseline=None, preload=True)
+    # print(epochs.get_data())
+    return epochs, raw
+
+
+
+
+def eegMCI_data(stim_type, subject, run_type, path_to_data_folder):
+    """
+    eegMCI_data does read 3d data from .mat file structured as ( n epochs, n channels, n samples ).
+
+    :param stim_type:  'AV','A','V'
+    :param subject: select a subject from dataset
+    :param run_type: chose run type of selected subject, possible types ( button.dat, test.dat_1, test.dat_2, train.dat )
+    :param path_to_data_folder: root path to data folder
+
+    :return: epochs, raw mne objects
+    """
+
+    stim_type_long = ['audiovisual','audio','visual']
+
+    if stim_type == 'AV':
+        stim_typeL = stim_type_long[0]
+    elif stim_type == 'A':
+        stim_typeL = stim_type_long[1]
+    else:
+        stim_typeL = stim_type_long[2]
+
+    targets = ['allNTARGETS','allTARGETS']
+    path = os.path.join(path_to_data_folder,stim_typeL,'s'+str(subject)+'_'+stim_type+'_'+run_type+'.mat')
+
+    mat = scipy.io.loadmat(path)
+
+    raw_no_targets = mat[targets[0]].transpose((0, 2, 1))
+    raw_targets = mat[targets[1]].transpose((0, 2, 1))
+
+    electrodes_l = np.concatenate(mat['electrodes'])
+    electrodes_l = [str(x).replace('[','').replace(']','').replace("'",'') for x in electrodes_l]
+    ch_types = ['eeg'] * 16
+    info = create_info(electrodes_l, ch_types=ch_types, sfreq=512)
+    info.set_montage('standard_1020')
+
+    epochs_no_target = EpochsArray(raw_no_targets, info, tmin=-0.2)
+    epochs_target = EpochsArray(raw_targets, info, tmin=-0.2)
+
+    erp_no_target = epochs_no_target.average()
+    erp_target = epochs_target.average()
+
+    epochs = np.array([erp_no_target.get_data(),erp_target.get_data()],dtype=np.object)
+    raw = [raw_no_targets,raw_targets]
+    print(electrodes_l)
+    return EpochsArray(epochs,info,tmin=-0.2), np.array(raw)
+
+
+
+
+def eegMCI_data_epochs(stim_type1,stim_type2, subject, run_type, path_to_data_folder):
+    """
+    eegMCI_data does read 3d data from .mat file structured as ( n epochs, n channels, n samples ).
+
+    :param stim_type1:  'AV','A','V'
+    :param stim_type2:  'AV','A','V'
+    :param subject: select a subject from dataset
+    :param run_type: chose run type of selected subject, possible types ( button.dat, test.dat_1, test.dat_2, train.dat )
+    :param path_to_data_folder: root path to data folder
+
+    :return: epochs, raw mne objects
+    """
+
+    stim_type_long = ['audiovisual','audio','visual']
+
+    if stim_type1 == 'AV':
+        stim_typeL1 = stim_type_long[0]
+    elif stim_type1 == 'A':
+        stim_typeL1 = stim_type_long[1]
+    else:
+        stim_typeL1 = stim_type_long[2]
+
+    if stim_type2 == 'AV':
+        stim_typeL2 = stim_type_long[0]
+    elif stim_type2 == 'A':
+        stim_typeL2 = stim_type_long[1]
+    else:
+        stim_typeL2 = stim_type_long[2]
+
+    targets = ['allNTARGETS','allTARGETS']
+    path1 = os.path.join(path_to_data_folder,stim_typeL1,'s'+str(subject)+'_'+stim_type1+'_'+run_type+'.mat')
+    path2 = os.path.join(path_to_data_folder,stim_typeL2,'s'+str(subject)+'_'+stim_type2+'_'+run_type+'.mat')
+
+    mat1 = scipy.io.loadmat(path1)
+    mat2 = scipy.io.loadmat(path2)
+
+    raw_no_targets1 = mat1[targets[0]].transpose((0, 2, 1))
+    raw_targets1 = mat1[targets[1]].transpose((0, 2, 1))
+
+    raw_no_targets2 = mat2[targets[0]].transpose((0, 2, 1))
+    raw_targets2 = mat2[targets[1]].transpose((0, 2, 1))
+
+    electrodes_l = np.concatenate(mat1['electrodes'])
+    electrodes_l = [str(x).replace('[','').replace(']','').replace("'",'') for x in electrodes_l]
+    ch_types = ['eeg'] * 16
+    info = create_info(electrodes_l, ch_types=ch_types, sfreq=512)
+    info.set_montage('standard_1020')
+
+    # epochs_no_target = EpochsArray(raw_no_targets1, info, tmin=-0.2)
+    epochs_target1 = EpochsArray(raw_targets1, info, tmin=-0.2, baseline=(-0.2,0))
+    epochs_target2 = EpochsArray(raw_targets2, info, tmin=-0.2, baseline=(-0.2,0))
+
+    # erp_no_target = epochs_no_target.average()
+    erp_target1 = epochs_target1.average()
+    erp_target2 = epochs_target2.average()
+
+    evoked = np.array([erp_target1.get_data(),erp_target2.get_data()],dtype=np.object)
+    epochs = [epochs_target1,epochs_target2]
+    print(electrodes_l)
+    return EpochsArray(evoked,info,tmin=-0.2), epochs

BrainPulse/dependency.py

@@ -0,0 +1,8 @@
+import pynndescent
+import numpy
+import pyrqa
+import pandas
+import sklearn
+import scipy
+import torch
+import umap

BrainPulse/distance_matrix.py

@@ -0,0 +1,104 @@
+from numba import jit
+from pynndescent.distances import wasserstein_1d, spearmanr, euclidean
+import numpy as np
+from pyrqa.computation import RPComputation
+from pyrqa.time_series import TimeSeries, EmbeddedSeries
+from pyrqa.settings import Settings
+from pyrqa.analysis_type import Classic
+from pyrqa.metric import EuclideanMetric
+# from pyrqa.metric import Sigmoid
+# from pyrqa.metric import Cosine
+from pyrqa.neighbourhood import Unthresholded
+import torch
+
+
+#
+# @jit(parallel=True)
+@jit(nopython=True)
+def EuclideanPyRQA_RP_stft_cpu(stft_):
+    result = np.zeros((stft_.shape[0], stft_.shape[0]))
+
+    for i, fft in enumerate(stft_):
+        for j, v in enumerate(stft_):
+            d = euclidean(fft, v)
+            result[i, j] = d
+    return result
+
+@jit(nopython=True)
+def wasserstein_squereform(stft_):
+    result = np.zeros((stft_.shape[0], stft_.shape[0]))
+
+    for i, fft in enumerate(stft_):
+        for j, v in enumerate(stft_):
+            d = wasserstein_1d(fft, v)
+            result[i, j] = d
+    return result
+
+@jit(nopython=True)
+def spearmanr_squereform(stft_):
+    result = np.zeros((stft_.shape[0], stft_.shape[0]))
+
+    for i, fft in enumerate(stft_):
+        for j, v in enumerate(stft_):
+            d = spearmanr(fft, v)
+            result[i, j] = d
+    return result
+
+@jit(nopython=True)
+def wasserstein_squereform_binary(stft_, eps_):
+    result = np.zeros((stft_.shape[0], stft_.shape[0]))
+
+    for i, fft in enumerate(stft_):
+        for j, v in enumerate(stft_):
+            d = wasserstein_1d(fft, v)
+            if d <= eps_:
+                dist = 1
+            else:
+                dist = 0
+
+            result[i, j] = dist
+    return result
+
+
+@jit(nopython=True)
+def wasserstein_1d_array(stft_):
+    result = np.zeros((stft_.shape[0] * stft_.shape[0]))
+    k = 0
+    for i, fft in enumerate(stft_):
+        for j, v in enumerate(stft_):
+            d = wasserstein_1d(fft, v)
+            result[k] = d
+            k += 1
+
+    return result
+
+def set_epsilon(matrix, eps):
+    return np.heaviside(eps - matrix, 0)
+
+
+def EuclideanPyRQA_RP(signal, embedding = 2,timedelay = 9):
+    time_series = TimeSeries(signal,
+                             embedding_dimension=embedding,
+                             time_delay=timedelay)
+    settings = Settings(time_series,
+                        analysis_type=Classic,
+                        neighbourhood=Unthresholded(),
+                        similarity_measure=EuclideanMetric,
+                        theiler_corrector=1)
+
+    computation = RPComputation.create(settings)
+    result = computation.run()
+    return result.recurrence_matrix
+
+
+def EuclideanPyRQA_RP_stft(signal, embedding = 2,timedelay = 9):
+    time_series = EmbeddedSeries(signal)
+    settings = Settings(time_series,
+                        analysis_type=Classic,
+                        neighbourhood=Unthresholded(),
+                        similarity_measure=EuclideanMetric,
+                        theiler_corrector=1)
+
+    computation = RPComputation.create(settings)
+    result = computation.run()
+    return result.recurrence_matrix

BrainPulse/event.py

@@ -0,0 +1,377 @@
+"""Event segmentation using a Hidden Markov Model
+
+Adapted from the brainiak package for this workshop.
+See https://brainiak.org/ for full documentation."""
+
+import numpy as np
+from scipy import stats
+import logging
+import copy
+from sklearn.base import BaseEstimator
+from sklearn.utils.validation import check_is_fitted, check_array
+from sklearn.exceptions import NotFittedError
+import itertools
+
+
+def masked_log(x):
+    y = np.empty(x.shape, dtype=x.dtype)
+    lim = x.shape[0]
+    for i in range(lim):
+        if x[i] <= 0:
+            y[i] = float('-inf')
+        else:
+            y[i] = np.log(x[i])
+    return y
+
+
+class EventSegment(BaseEstimator):
+
+    def _default_var_schedule(step):
+        return 4 * (0.98 ** (step - 1))
+
+    def __init__(self, n_events=2,
+                 step_var=_default_var_schedule,
+                 n_iter=500, event_chains=None,
+                 split_merge=False, split_merge_proposals=1):
+        self.n_events = n_events
+        self.step_var = step_var
+        self.n_iter = n_iter
+        self.split_merge = split_merge
+        self.split_merge_proposals = split_merge_proposals
+        if event_chains is None:
+            self.event_chains = np.zeros(n_events)
+        else:
+            self.event_chains = event_chains
+
+    def _fit_validate(self, X):
+        if len(np.unique(self.event_chains)) > 1:
+            raise RuntimeError("Cannot fit chains, use set_event_patterns")
+
+        # Copy X into a list and transpose
+        X = copy.deepcopy(X)
+        if type(X) is not list:
+            X = [X]
+        for i in range(len(X)):
+            X[i] = check_array(X[i])
+            X[i] = X[i].T
+
+        # Check that number of voxels is consistent across datasets
+        n_dim = X[0].shape[0]
+        for i in range(len(X)):
+            assert (X[i].shape[0] == n_dim)
+
+        # Double-check that data is z-scored in time
+        for i in range(len(X)):
+            X[i] = stats.zscore(X[i], axis=1, ddof=1)
+
+        return X
+
+    def fit(self, X, y=None):
+        seg = []
+        X = self._fit_validate(X)
+        n_train = len(X)
+        n_dim = X[0].shape[0]
+        self.classes_ = np.arange(self.n_events)
+
+        # Initialize variables for fitting
+        log_gamma = []
+        for i in range(n_train):
+            log_gamma.append(np.zeros((X[i].shape[1], self.n_events)))
+        step = 1
+        best_ll = float("-inf")
+        self.ll_ = np.empty((0, n_train))
+        while step <= self.n_iter:
+            iteration_var = self.step_var(step)
+
+            # Based on the current segmentation, compute the mean pattern
+            # for each event
+            seg_prob = [np.exp(lg) / np.sum(np.exp(lg), axis=0)
+                        for lg in log_gamma]
+            mean_pat = np.empty((n_train, n_dim, self.n_events))
+            for i in range(n_train):
+                mean_pat[i, :, :] = X[i].dot(seg_prob[i])
+            mean_pat = np.mean(mean_pat, axis=0)
+
+            # Based on the current mean patterns, compute the event
+            # segmentation
+            self.ll_ = np.append(self.ll_, np.empty((1, n_train)), axis=0)
+            for i in range(n_train):
+                logprob = self._logprob_obs(X[i], mean_pat, iteration_var)
+                log_gamma[i], self.ll_[-1, i] = self._forward_backward(logprob)
+
+            if step > 1 and self.split_merge:
+                curr_ll = np.mean(self.ll_[-1, :])
+                self.ll_[-1, :], log_gamma, mean_pat = \
+                    self._split_merge(X, log_gamma, iteration_var, curr_ll)
+
+            # If log-likelihood has started decreasing, undo last step and stop
+            if np.mean(self.ll_[-1, :]) < best_ll:
+                self.ll_ = self.ll_[:-1, :]
+                break
+
+            self.segments_ = [np.exp(lg) for lg in log_gamma]
+            self.event_var_ = iteration_var
+            self.event_pat_ = mean_pat
+            best_ll = np.mean(self.ll_[-1, :])
+            
+            seg.append(self.segments_[0].copy())
+
+            step += 1
+
+        return seg
+
+    def _logprob_obs(self, data, mean_pat, var):
+        n_vox = data.shape[0]
+        t = data.shape[1]
+
+        # z-score both data and mean patterns in space, so that Gaussians
+        # are measuring Pearson correlations and are insensitive to overall
+        # activity changes
+        data_z = stats.zscore(data, axis=0, ddof=1)
+        mean_pat_z = stats.zscore(mean_pat, axis=0, ddof=1)
+
+        logprob = np.empty((t, self.n_events))
+
+        if type(var) is not np.ndarray:
+            var = var * np.ones(self.n_events)
+
+        for k in range(self.n_events):
+            logprob[:, k] = -0.5 * n_vox * np.log(
+                2 * np.pi * var[k]) - 0.5 * np.sum(
+                (data_z.T - mean_pat_z[:, k]).T ** 2, axis=0) / var[k]
+
+        logprob /= n_vox
+        return logprob
+
+    def _forward_backward(self, logprob):
+        logprob = copy.copy(logprob)
+        t = logprob.shape[0]
+        logprob = np.hstack((logprob, float("-inf") * np.ones((t, 1))))
+
+        # Initialize variables
+        log_scale = np.zeros(t)
+        log_alpha = np.zeros((t, self.n_events + 1))
+        log_beta = np.zeros((t, self.n_events + 1))
+
+        # Set up transition matrix, with final sink state
+        self.p_start = np.zeros(self.n_events + 1)
+        self.p_end = np.zeros(self.n_events + 1)
+        self.P = np.zeros((self.n_events + 1, self.n_events + 1))
+        label_ind = np.unique(self.event_chains, return_inverse=True)[1]
+        n_chains = np.max(label_ind) + 1
+
+        # For each chain of events, link them together and then to sink state
+        for c in range(n_chains):
+            chain_ind = np.nonzero(label_ind == c)[0]
+            self.p_start[chain_ind[0]] = 1 / n_chains
+            self.p_end[chain_ind[-1]] = 1 / n_chains
+
+            p_trans = (len(chain_ind) - 1) / t
+            if p_trans >= 1:
+                raise ValueError('Too few timepoints')
+            for i in range(len(chain_ind)):
+                self.P[chain_ind[i], chain_ind[i]] = 1 - p_trans
+                if i < len(chain_ind) - 1:
+                    self.P[chain_ind[i], chain_ind[i+1]] = p_trans
+                else:
+                    self.P[chain_ind[i], -1] = p_trans
+        self.P[-1, -1] = 1
+
+        # Forward pass
+        for i in range(t):
+            if i == 0:
+                log_alpha[0, :] = self._log(self.p_start) + logprob[0, :]
+            else:
+                log_alpha[i, :] = self._log(np.exp(log_alpha[i - 1, :])
+                                            .dot(self.P)) + logprob[i, :]
+
+            log_scale[i] = np.logaddexp.reduce(log_alpha[i, :])
+            log_alpha[i] -= log_scale[i]
+
+        # Backward pass
+        log_beta[-1, :] = self._log(self.p_end) - log_scale[-1]
+        for i in reversed(range(t - 1)):
+            obs_weighted = log_beta[i + 1, :] + logprob[i + 1, :]
+            offset = np.max(obs_weighted)
+            log_beta[i, :] = offset + self._log(
+                np.exp(obs_weighted - offset).dot(self.P.T)) - log_scale[i]
+
+        # Combine and normalize
+        log_gamma = log_alpha + log_beta
+        log_gamma -= np.logaddexp.reduce(log_gamma, axis=1, keepdims=True)
+
+        ll = np.sum(log_scale[:(t - 1)]) + np.logaddexp.reduce(
+            log_alpha[-1, :] + log_scale[-1] + self._log(self.p_end))
+
+        log_gamma = log_gamma[:, :-1]
+
+        return log_gamma, ll
+
+    def _log(self, x):
+        xshape = x.shape
+        _x = x.flatten()
+        y = masked_log(_x)
+        return y.reshape(xshape)
+
+    def set_event_patterns(self, event_pat):
+        if event_pat.shape[1] != self.n_events:
+            raise ValueError(("Number of columns of event_pat must match "
+                              "number of events"))
+        self.event_pat_ = event_pat.copy()
+
+    def find_events(self, testing_data, var=None, scramble=False):
+        if var is None:
+            if not hasattr(self, 'event_var_'):
+                raise NotFittedError(("Event variance must be provided, if "
+                                      "not previously set by fit()"))
+            else:
+                var = self.event_var_
+
+        if not hasattr(self, 'event_pat_'):
+            raise NotFittedError(("The event patterns must first be set "
+                                  "by fit() or set_event_patterns()"))
+        if scramble:
+            mean_pat = self.event_pat_[:, np.random.permutation(self.n_events)]
+        else:
+            mean_pat = self.event_pat_
+
+        logprob = self._logprob_obs(testing_data.T, mean_pat, var)
+        lg, test_ll = self._forward_backward(logprob)
+        segments = np.exp(lg)
+
+        return segments, test_ll
+
+    def predict(self, X):
+        check_is_fitted(self, ["event_pat_", "event_var_"])
+        X = check_array(X)
+        segments, test_ll = self.find_events(X)
+        return np.argmax(segments, axis=1)
+
+    def calc_weighted_event_var(self, D, weights, event_pat):
+        Dz = stats.zscore(D, axis=1, ddof=1)
+        ev_var = np.empty(event_pat.shape[1])
+        for e in range(event_pat.shape[1]):
+            # Only compute variances for weights > 0.1% of max weight
+            nz = weights[:, e] > np.max(weights[:, e])/1000
+            sumsq = np.dot(weights[nz, e],
+                           np.sum(np.square(Dz[nz, :] -
+                                  event_pat[:, e]), axis=1))
+            ev_var[e] = sumsq/(np.sum(weights[nz, e]) -
+                               np.sum(np.square(weights[nz, e])) /
+                               np.sum(weights[nz, e]))
+        ev_var = ev_var / D.shape[1]
+        return ev_var
+
+    def model_prior(self, t):
+        lg, test_ll = self._forward_backward(np.zeros((t, self.n_events)))
+        segments = np.exp(lg)
+
+        return segments, test_ll
+
+    def _split_merge(self, X, log_gamma, iteration_var, curr_ll):
+        # Compute current probabilities and mean patterns
+        n_train = len(X)
+        n_dim = X[0].shape[0]
+
+        seg_prob = [np.exp(lg) / np.sum(np.exp(lg), axis=0)
+                    for lg in log_gamma]
+        mean_pat = np.empty((n_train, n_dim, self.n_events))
+        for i in range(n_train):
+            mean_pat[i, :, :] = X[i].dot(seg_prob[i])
+        mean_pat = np.mean(mean_pat, axis=0)
+
+        # For each event, merge its probability distribution
+        # with the next event, and also split its probability
+        # distribution at its median into two separate events.
+        # Use these new event probability distributions to compute
+        # merged and split event patterns.
+        merge_pat = np.empty((n_train, n_dim, self.n_events))
+        split_pat = np.empty((n_train, n_dim, 2 * self.n_events))
+        for i, sp in enumerate(seg_prob):  # Iterate over datasets
+            m_evprob = np.zeros((sp.shape[0], sp.shape[1]))
+            s_evprob = np.zeros((sp.shape[0], 2 * sp.shape[1]))
+            cs = np.cumsum(sp, axis=0)
+            for e in range(sp.shape[1]):
+                # Split distribution at midpoint and normalize each half
+                mid = np.where(cs[:, e] >= 0.5)[0][0]
+                cs_first = cs[mid, e] - sp[mid, e]
+                cs_second = 1 - cs_first
+                s_evprob[:mid, 2 * e] = sp[:mid, e] / cs_first
+                s_evprob[mid:, 2 * e + 1] = sp[mid:, e] / cs_second
+
+                # Merge distribution with next event distribution
+                m_evprob[:, e] = sp[:, e:(e + 2)].mean(1)
+
+            # Weight data by distribution to get event patterns
+            merge_pat[i, :, :] = X[i].dot(m_evprob)
+            split_pat[i, :, :] = X[i].dot(s_evprob)
+
+        # Average across datasets
+        merge_pat = np.mean(merge_pat, axis=0)
+        split_pat = np.mean(split_pat, axis=0)
+
+        # Correlate the current event patterns with the split and
+        # merged patterns
+        merge_corr = np.zeros(self.n_events)
+        split_corr = np.zeros(self.n_events)
+        for e in range(self.n_events):
+            split_corr[e] = np.corrcoef(mean_pat[:, e],
+                                        split_pat[:, (2 * e):(2 * e + 2)],
+                                        rowvar=False)[0, 1:3].max()
+            merge_corr[e] = np.corrcoef(merge_pat[:, e],
+                                        mean_pat[:, e:(e + 2)],
+                                        rowvar=False)[0, 1:3].min()
+        merge_corr = merge_corr[:-1]
+
+        # Find best merge/split candidates
+        # A high value of merge_corr indicates that a pair of events are
+        # very similar to their merged pattern, and are good candidates for
+        # being merged.
+        # A low value of split_corr indicates that an event's pattern is
+        # very dissimilar from the patterns in its first and second half,
+        # and is a good candidate for being split.
+        best_merge = np.flipud(np.argsort(merge_corr))
+        best_merge = best_merge[:self.split_merge_proposals]
+        best_split = np.argsort(split_corr)
+        best_split = best_split[:self.split_merge_proposals]
+
+        # For every pair of merge/split candidates, attempt the merge/split
+        # and measure the log-likelihood. If any are better than curr_ll,
+        # accept this best merge/split
+        mean_pat_last = mean_pat.copy()
+        return_ll = curr_ll
+        return_lg = copy.deepcopy(log_gamma)
+        return_mp = mean_pat.copy()
+        for m_e, s_e in itertools.product(best_merge, best_split):
+            if m_e == s_e or m_e+1 == s_e:
+                # Don't attempt to merge/split same event
+                continue
+
+            # Construct new set of patterns with merge/split
+            mean_pat_ms = np.delete(mean_pat_last, s_e, axis=1)
+            mean_pat_ms = np.insert(mean_pat_ms, [s_e, s_e],
+                                    split_pat[:, (2 * s_e):(2 * s_e + 2)],
+                                    axis=1)
+            mean_pat_ms = np.delete(mean_pat_ms,
+                                    [m_e + (s_e < m_e), m_e + (s_e < m_e) + 1],
+                                    axis=1)
+            mean_pat_ms = np.insert(mean_pat_ms, m_e + (s_e < m_e),
+                                    merge_pat[:, m_e], axis=1)
+
+            # Measure log-likelihood with these new patterns
+            ll_ms = np.zeros(n_train)
+            log_gamma_ms = list()
+            for i in range(n_train):
+                logprob = self._logprob_obs(X[i],
+                                            mean_pat_ms, iteration_var)
+                lg, ll_ms[i] = self._forward_backward(logprob)
+                log_gamma_ms.append(lg)
+
+            # If better than best ll so far, save to return to fit()
+            if ll_ms.mean() > return_ll:
+                return_mp = mean_pat_ms.copy()
+                return_ll = ll_ms
+                for i in range(n_train):
+                    return_lg[i] = log_gamma_ms[i].copy()
+
+        return return_ll, return_lg, return_mp

BrainPulse/features_space.py

@@ -0,0 +1,282 @@
+import glob
+import numpy as np
+import pandas as pd
+from sklearn import svm
+import sklearn.model_selection as model_selection
+from sklearn.metrics import accuracy_score
+from sklearn.metrics import confusion_matrix
+from sklearn.impute import SimpleImputer
+from sklearn.svm import SVC
+from sklearn.model_selection import StratifiedKFold
+from sklearn.feature_selection import RFECV
+from sklearn.ensemble import RandomForestClassifier
+import matplotlib.pyplot as plt
+
+# import rcr
+
+def save_features_as_csv():
+    return
+
+def load_features_csv_concat(folder_path):
+    df_list = []
+    for file in glob.glob(folder_path+"/*.csv"):
+        df_ = pd.read_csv(file)
+        df_list.append(df_)
+    df = pd.concat(df_list)
+    df = df.reset_index(drop=True)
+    return df
+
+def exclude_subject(df,exluded_subjects):
+    condition_string = ''
+    for ex_sub in exluded_subjects:
+        condition_string += "(df['Subject'] !='" +ex_sub+"') & "
+    evaluation_string = 'df['+condition_string[:len(condition_string)-2]+']'
+    df_ex = eval(evaluation_string)
+    return df_ex.reset_index(drop=True)
+
+def electrode_wise_dataframe(df, condition_list, id_vars = ['Subject', 'Task', 'Electrode']):
+    stats_frame = df[
+        ['Subject', 'Task', 'Electrode','Lentr', 'TT', 'L', 'RR', 'LAM', 'DET', 'V','Vmax', 'Ventr', 'W','Wentr']
+    ]
+
+    stats_frame.melt(id_vars=id_vars, var_name='RQA_feature', value_name='feature_value')
+    stats = stats_frame.pivot_table(index=['Subject', 'Task'], columns='Electrode',
+                                     values=['Lentr', 'TT', 'L', 'RR', 'LAM', 'DET', 'V','Vmax', 'Ventr', 'W', 'Wentr']).reset_index()
+
+    stats = stats.replace(condition_list[0], 0)
+    stats = stats.replace(condition_list[1], 1)
+    y = stats.Task.values
+    return stats, y
+
+
+def electrode_wise_dataframe_epochs(df, condition_list, id_vars = ['Subject', 'Task', 'Epoch_id','Electrode']):
+    stats_frame = df[
+        ['Subject', 'Task','Epoch_id','Electrode','Lentr', 'TT', 'L', 'RR', 'LAM', 'DET', 'V','Vmax', 'Ventr', 'W','Wentr']
+    ]
+
+    stats_frame.melt(id_vars=id_vars, var_name='RQA_feature', value_name='feature_value')
+    stats = stats_frame.pivot_table(index=['Subject', 'Task'], columns=['Electrode', 'Epoch_id'],
+                                    values=['Lentr', 'TT', 'L', 'RR', 'LAM', 'DET', 'V','Vmax', 'Ventr', 'W', 'Wentr']).reset_index()
+
+    stats = stats.replace(condition_list[0], 0)
+    stats = stats.replace(condition_list[1], 1)
+    y = stats.Task.values
+    return stats, y
+
+
+def select_features_clean_and_normalize(df,features=['Lentr', 'TT', 'L', 'LAM', 'DET','V', 'Ventr', 'W','Wentr']):
+
+    stats_data = df[features].values
+
+    #rcr
+    stats_data_cleaned=np.empty((stats_data.shape[0],stats_data.shape[1]))
+    stats_data_cleaned[:]=np.nan
+    # r = rcr.RCR(rcr.SS_MEDIAN_DL)
+    r = stats_data_cleaned
+
+    for ii in range(stats_data.shape[1]):
+        # fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,figsize=(16,8),dpi=200)
+        # ax1.hist(stats_data[:,ii])
+        # ax1.set_title('Raw')
+
+        r.performBulkRejection(stats_data[:,ii])
+        cleaned_data_indices = r.result.indices
+        stats_data_cleaned[cleaned_data_indices,ii]=stats_data[cleaned_data_indices,ii]
+
+        # ax2.hist(stats_data_cleaned[:,ii][~np.isnan(stats_data_cleaned[:,ii])])
+        # ax2.set_title('Cleaned')
+
+        # plt.savefig('Feature_nr_'+str(ii)+'jpg')
+        # plt.close()
+
+
+    df_stats_data_cleaned=pd.DataFrame(stats_data_cleaned)
+    # df_stats_data_cleaned=df_stats_data_cleaned.fillna(method='mean', axis=0)#+df_stats_data_cleaned.fillna(method='bfill', axis=0))/2
+    # df_stats_data_cleaned.interpolate(limit=5, inplace=True)
+
+    imputer = SimpleImputer(missing_values=np.nan, strategy='mean')
+    imputer = imputer.fit(df_stats_data_cleaned)
+
+    stats_data_cleaned = imputer.transform(df_stats_data_cleaned)
+
+    ####normalize#########
+    stats_data_normed=np.empty((stats_data.shape[0],stats_data.shape[1]))
+    for ii in range(stats_data.shape[1]):
+        stats_data_normed[:,ii] = (stats_data_cleaned[:,ii]-stats_data_cleaned[:,ii].min(axis=0))/ (stats_data_cleaned[:,ii].max(axis=0)-stats_data_cleaned[:,ii].min(axis=0)) #stats_data[:,ii]-stats_data[:,ii].mean(axis=0))/ stats_data[:,ii].std(axis=0)
+
+    return stats_data_normed
+
+
+def clasyfication_SVM(df,y,cv=10,type='linear'):
+
+
+    clf=svm.SVC(kernel=type)
+    skf = StratifiedKFold(n_splits=cv)
+    # run split() again to generate folds
+    folds = skf.split(df, y)
+    print('folds shape ', folds)
+    performance = np.zeros(skf.n_splits)
+    performance_open= np.zeros(skf.n_splits)
+    performance_closed= np.zeros(skf.n_splits)
+
+    for i, (train_idx, test_idx) in enumerate(folds):
+
+        X_train = df[train_idx,:]
+        y_train = y[train_idx]
+
+        X_test = df[test_idx,:]
+        y_test = y[test_idx]
+
+        # call fit (on train) and predict (on test)
+        model = clf.fit(X=X_train, y=y_train)
+        y_hat = model.predict(X=X_test)
+
+        # calculate accuracy
+        performance[i] = accuracy_score(y_test, y_hat)
+        cm = confusion_matrix(y_test, y_hat)
+        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
+        # class_acuracy = cm.diagonal()
+        class_acuracy = cm.diagonal()
+        performance_open[i]=class_acuracy[0]*100
+        performance_closed[i]=class_acuracy[1]*100
+
+    # calculate average accuracy
+    print('Mean performance: %.3f' % np.mean(performance*100))
+    print('Mean performance 1st class: %.3f' % np.mean(performance_open))
+    print('Mean performance 2nd class: %.3f' % np.mean(performance_closed))
+
+
+    lin = svm.SVC(kernel=type).fit(X_train, y_train)
+    lin_pred = lin.predict(X_test)
+
+    return lin, lin_pred
+
+def cross_validation(df,y,cv=10,title = 'cv job',type='linear'):
+
+    # Create the RFE object and compute a cross-validated score.
+    svc = SVC(kernel=type)
+    # The "accuracy" scoring shows the proportion of correct classifications
+
+    min_features_to_select = 4  # Minimum number of features to consider
+    rfecv = RFECV(
+        estimator=svc,
+        step=1,
+        cv=StratifiedKFold(n_splits=cv),
+        scoring="accuracy",
+        min_features_to_select=min_features_to_select,
+    )
+    rfecv.fit(df, y)
+
+    print("Optimal number of features : %d" % rfecv.n_features_)
+
+    # Plot number of features VS. cross-validation scores
+    plt.figure()
+    plt.xlabel("Number of features selected")
+    plt.ylabel("Cross validation score (accuracy)")
+    plt.plot(
+        range(min_features_to_select, len(rfecv.cv_results_['mean_test_score']) + min_features_to_select),
+        rfecv.cv_results_['mean_test_score'],
+    )
+    plt.title(title)
+    plt.show()
+    plt.savefig(title+' classification with feature selection_more_features'+str(rfecv.n_features_)+'_'+str(round(max(rfecv.cv_results_['mean_test_score'])*100,2))+'.png', dpi=150)
+    plt.close()
+
+    return rfecv.transform(df)
+
+
+
+def compute_binary_SVM(df,y,predict_on_all_data = False,type='linear'):
+
+    # stats_data = df[['TT', 'RR', 'DET', 'LAM', 'L', 'Lentr']].values
+    X_train, X_test, y_train, y_test = model_selection.train_test_split(df, y, train_size=0.80, test_size=0.20,
+                                                                        random_state=101)
+    global lin
+
+    if predict_on_all_data:
+        print('SVM prediction on all data')
+        lin = svm.SVC(kernel=type).fit(X_train, y_train)
+
+        lin_pred = lin.predict(df)
+
+        lin_accuracy = accuracy_score(y, lin_pred)
+
+        print('Accuracy (Linear Kernel): ', "%.2f" % (lin_accuracy * 100))
+
+        cm = confusion_matrix(y, lin_pred)
+        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
+        class_acuracy = cm.diagonal()
+        print('Accuracy (1st class): ', "%.2f" % (class_acuracy[0] * 100))
+        print('Accuracy (2nd class): ', "%.2f" % (class_acuracy[1] * 100))
+    else:
+        print('SVM prediction on test data')
+        lin = svm.SVC(kernel=type).fit(X_train, y_train)
+
+        lin_pred = lin.predict(X_test)
+
+        lin_accuracy = accuracy_score(y_test, lin_pred)
+
+        print('Accuracy (Linear Kernel): ', "%.2f" % (lin_accuracy * 100))
+        print('Y train:', y_train)
+        print('Y test:', y_test)
+        print('Y pred:', lin_pred)
+
+        cm = confusion_matrix(y_test, lin_pred)
+        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
+        class_acuracy = cm.diagonal()
+        print('Accuracy (1st class): ', "%.2f" % (class_acuracy[0] * 100))
+        print('Accuracy (2nd class): ', "%.2f" % (class_acuracy[1] * 100))
+
+    return lin, lin_pred
+
+
+
+def clasyfication_RFC(df,y,cv=10,max_depth=2):
+
+    clf = RandomForestClassifier(max_depth=max_depth, random_state=0)
+    skf = StratifiedKFold(n_splits=cv)
+    # run split() again to generate folds
+    folds = skf.split(df, y)
+
+    performance = np.zeros(skf.n_splits)
+    performance_open= np.zeros(skf.n_splits)
+    performance_closed= np.zeros(skf.n_splits)
+
+    for i, (train_idx, test_idx) in enumerate(folds):
+
+        X_train = df[train_idx,:]
+        y_train = y[train_idx]
+
+        X_test = df[test_idx,:]
+        y_test = y[test_idx]
+
+        # call fit (on train) and predict (on test)
+        model = clf.fit(X=X_train, y=y_train)
+        y_hat = model.predict(X=X_test)
+
+        # calculate accuracy
+        performance[i] = accuracy_score(y_test, y_hat)
+        cm = confusion_matrix(y_test, y_hat)
+        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
+        # class_acuracy = cm.diagonal()
+        class_acuracy = cm.diagonal()
+        performance_open[i]=class_acuracy[0]*100
+        performance_closed[i]=class_acuracy[1]*100
+
+    # calculate average accuracy
+    print('Mean performance: %.3f' % np.mean(performance*100))
+    print('Mean performance 1st class: %.3f' % np.mean(performance_open))
+    print('Mean performance 2nd class: %.3f' % np.mean(performance_closed))
+
+
+    lin = RandomForestClassifier(max_depth=max_depth, random_state=0)
+    lin.fit(X=X_train, y=y_train)
+    lin_pred = lin.predict(X_test)
+
+    return lin, lin_pred
+
+
+
+
+
+

BrainPulse/frequency_recurrence.py

@@ -0,0 +1,124 @@
+import matplotlib.pyplot as plt
+# %matplotlib inline
+# %load_ext autotime
+plt.style.use('classic')
+import numpy as np
+
+
+
+def normalize(tSignal):
+    # copy the data if needed, omit and rename function argument if desired
+    signal = np.copy(tSignal) # signal is in range [a;b]
+    signal -= np.min(signal) # signal is in range to [0;b-a]
+    signal /= np.max(signal) # signal is normalized to [0;1]
+    signal -= 0.5 # signal is in range [-0.5;0.5]
+    signal *=2 # signal is in range [-1;1]
+    return signal
+
+
+
+
+def get_max_freqs(stft, cut_freq, crop=False, norm=True):
+    if crop==False:
+        crop=len(stft)
+    id_max_stft = stft[:crop].argmax(axis=1)*(cut_freq/stft.shape[1])
+    if norm:
+        id_max_norm=id_max_stft*(1/cut_freq)
+        return id_max_norm
+    else: 
+        return id_max_stft
+
+
+def symmetrize(a):
+    return a + np.transpose(a, (1, 0, 2)) - np.diag(a.diagonal())
+
+def symmetrized(m):
+    
+    import numpy as np
+    
+    i_lower = np.tril_indices(m.shape[0], -1)
+    m[:,:,0][i_lower] = m[:,:,0].T[i_lower] 
+    m[:,:,1][i_lower] = m[:,:,1].T[i_lower]
+    m[:,:,2][i_lower] = m[:,:,2].T[i_lower] 
+
+    return m
+
+def calc_color(id_max_norm, unique_values,colors,crop):
+    if crop==False:
+        crop=len(id_max_norm)
+    tmp = np.zeros((id_max_norm.shape[0],id_max_norm.shape[0],3), dtype = 'float64')
+
+    for i, v1 in enumerate(id_max_norm):
+        for j, v2 in enumerate(id_max_norm):
+            tmp[i,j] = v1,v2,0
+
+    tmp_color = np.zeros((id_max_norm.shape[0],id_max_norm.shape[0],3), dtype = 'float64')
+    for i in range(crop):
+        for j in range(crop):
+            x = tmp[i,j][0]
+            y = tmp[i,j][1]
+
+            for k in unique_values:
+
+                if x == k:
+                    id = (unique_values).index(x)
+                    
+                    tmp_color[i,j][0] = colors[id][0]
+                if y == k:
+                    id = (unique_values).index(y)
+                    tmp_color[i,j][1] = colors[id][1]
+                    
+     
+    return symmetrized(tmp_color),symmetrized(tmp)
+
+def calc_color_raw(id_max_stft):
+
+    tmp = np.zeros((id_max_stft.shape[0],id_max_stft.shape[0],3), dtype = 'float64')
+
+    for i, v1 in enumerate(id_max_stft):
+        for j, v2 in enumerate(id_max_stft):
+            tmp[i,j] = v1,v2,0
+    return symmetrized(tmp)
+
+def get_unique_colors(cut_freq,unique_values):
+    from matplotlib.colors import LinearSegmentedColormap
+    
+    vmax=cut_freq
+    cmap = LinearSegmentedColormap.from_list('mycmap1', [(0 / vmax, "violet"),
+                                                        (4 / vmax, 'blue'),
+                                                        (8 / vmax, 'green'),
+                                                        (15 / vmax, 'yellow'),
+                                                        (30 / vmax, 'red'),
+                                                        (60 / vmax, 'black')
+                                                        ])
+
+    colors=[cmap(each) for each in unique_values]
+    
+    return cmap, colors
+
+
+def freqRP(stft, cut_freq, crop=False, norm=True):
+
+    id_max_stft=get_max_freqs(stft,cut_freq,crop,norm)
+
+    unique_values = np.unique(id_max_stft).tolist()
+
+    cmap,colors=get_unique_colors(cut_freq, unique_values)
+
+    color_matrix, raw_c_matrix = calc_color(id_max_stft, unique_values,colors,crop)
+
+    return color_matrix, raw_c_matrix
+
+
+def plot_freqRP_interactive(color_matrix, raw_c_matrix, filename='', save=False):
+    import plotly.express as px
+    import numpy as np
+    fig=px.imshow(color_matrix, origin='lower')
+    fig.update_traces(customdata=np.round((raw_c_matrix*60),2),
+        hovertemplate="First frequency: %{customdata[0]}<br>Second frequency: %{customdata[1]}<extra></extra>"
+    )
+    if save:
+        fig.write_html(filename)
+    else:
+        fig.show()
+

BrainPulse/matrix_open_binary.npy

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:db7cd1f9b9fea62c6ecda8ec53c0e3759d6902e82590770080ef5c24531a8c0b
+size 92236944

BrainPulse/model_SVM.py

@@ -0,0 +1,45 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+import seaborn as sns
+from sklearn import svm, datasets
+import sklearn.model_selection as model_selection
+from sklearn.metrics import accuracy_score
+from sklearn.metrics import f1_score
+from sklearn.metrics import confusion_matrix
+from sklearn.preprocessing import StandardScaler
+import umap
+import umap.plot
+from umap.parametric_umap import ParametricUMAP
+from sklearn.impute import SimpleImputer
+from sklearn.pipeline import make_pipeline
+from sklearn.preprocessing import QuantileTransformer
+from mne.datasets import eegbci
+
+def compute_Linear_Kernel(df):
+    stats_frame = df[
+        ['Subject', 'Task', 'Electrode', 'Lentr', 'TT', 'L', 'RR', 'LAM', 'DET']
+    ]
+
+    stats_frame.melt(id_vars=['Subject', 'Task', 'Electrode'], var_name='RQA_feature', value_name='feature_value')
+    stats1 = stats_frame.pivot_table(index=['Subject', 'Task'], columns='Electrode',
+                                     values=['TT', 'RR', 'DET', 'LAM', 'L', 'Lentr']).reset_index()
+
+    y = stats1.Task.values
+    stats_data = stats1[['TT', 'RR', 'DET', 'LAM', 'L', 'Lentr']].values
+    X_train, X_test, y_train, y_test = model_selection.train_test_split(stats_data, y, train_size=0.80, test_size=0.20,
+                                                                        random_state=101)
+
+    lin = svm.SVC(kernel='linear').fit(X_train, y_train)
+
+    lin_pred = lin.predict(X_test)
+
+    lin_accuracy = accuracy_score(y_test, lin_pred)
+
+    print('Accuracy (Linear Kernel): ', "%.2f" % (lin_accuracy * 100))
+
+    cm = confusion_matrix(y_test, lin_pred)
+    cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
+    class_acuracy = cm.diagonal()
+    print('Accuracy (open): ', "%.2f" % (class_acuracy[0] * 100))
+    print('Accuracy (close): ', "%.2f" % (class_acuracy[1] * 100))

BrainPulse/plot.py

@@ -0,0 +1,634 @@
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import animation, cm
+import matplotlib
+import seaborn as sns
+import seaborn as sns
+import umap
+import umap.plot
+import pandas as pd
+from .event import EventSegment
+from sklearn.impute import SimpleImputer
+from sklearn.pipeline import make_pipeline
+from sklearn.preprocessing import QuantileTransformer
+
+sns.set_style("whitegrid")
+
+# plt.rcParams["font.family"] = "cursive"
+# plt.rcParams.update({'font.sans-serif':'Times'})
+# plt.rcParams.update({'font.family':'sans-serif'})
+# plt.rcParams['font.size'] = 14
+import matplotlib.font_manager as font_manager
+font = font_manager.FontProperties(family='Times')
+
+def explainer_(chan, stft, cut_freq, s_rate):
+
+    fig, axs = plt.subplots(4, figsize=(10, 14), dpi=150)  # figsize=(12, 12),
+    time_crop = np.linspace(0, int(chan[:400].shape[0]), chan[:400].shape[0])
+
+    axs[0].plot(chan[:400],'k')  # np.linspace(0, int(chan[:400].shape[0]/s_rate), chan[:400].shape[0]),
+    axs[0].fill_betweenx(y=[-210, 125], x1=time_crop[0],
+                        x2=time_crop[240], color='white', alpha=0.9, edgecolor='red' )
+
+    axs[0].fill_betweenx(y=[-210, 130], x1=time_crop[2]+20,
+                         x2=time_crop[260], color='white', alpha=0.9, edgecolor='green')
+
+    axs[0].fill_betweenx(y=[-210, 135], x1=time_crop[2]+40,
+                         x2=time_crop[280], color='white', alpha=0.9, edgecolor='blue')
+
+    axs[0].annotate('$fft_{1}$', xy=(.25, 72), xycoords='data',
+                xytext=(0.05, 1.45), textcoords='axes fraction',
+                arrowprops=dict(arrowstyle="->",facecolor='black',color='black'),
+                horizontalalignment='right', verticalalignment='top',
+                )
+
+    axs[0].annotate('$fft_{2}$', xy=(23.35, 85), xycoords='data',
+                    xytext=(0.15, 1.45), textcoords='axes fraction',
+                    arrowprops=dict(arrowstyle="->",facecolor='black',color='black'),
+                    horizontalalignment='right', verticalalignment='top',
+                    )
+
+    axs[0].annotate('$fft_{3}$', xy=(43.45, 95), xycoords='data',
+                    xytext=(0.25, 1.45), textcoords='axes fraction',
+                    arrowprops=dict(arrowstyle="->",facecolor='black ',color='black'),
+                    horizontalalignment='right', verticalalignment='top',
+                    )
+    axs[0].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, chan[:400].shape[0], 5)))
+    axs[0].set_xticklabels(
+        [str(np.round(x, 1)) for x in np.linspace(0, int(chan[:400].shape[0] / s_rate), axs[0].get_xticks().shape[0])])
+    axs[0].set_ylabel('Amplitude (µV)', )
+    axs[0].set_xlabel('Time (s)', )
+    axs[0].set_title('(a)', )
+    axs[0].xaxis.grid()
+    axs[0].yaxis.grid()
+
+
+    axs[1].plot((stft[100]/stft.shape[1])**2, 'red',label='$fft_{1}$',marker="o",markersize=3)
+    axs[1].legend(prop=font)
+    axs[1].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+    axs[1].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+    axs[1].set_xlim([0, 100])
+    # axs[1].set_ylim([0, 250])
+    axs[1].set_ylabel('Power ($\mu V^{2}$)', )
+    axs[1].set_xlabel('Freq (Hz)', )
+    # axs[1].set_title('Frequency Domain ($fft_{1}$, $fft_{2}$, $fft_{3}$)', fontsize=10)
+    axs[1].set_title('(b)', )
+    axs[1].xaxis.grid()
+    axs[1].yaxis.grid()
+
+
+    axs[2].plot((stft[115]/stft.shape[1])**2, 'green', label='$fft_{2}$', marker="o", markersize=3)
+    axs[2].legend(prop=font)
+    axs[2].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+    axs[2].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+    axs[2].set_xlim([0, 100])
+    # axs[2].set_ylim([0, 250])
+    axs[2].set_ylabel('Power ($\mu V^{2}$)', )
+    axs[2].set_xlabel('Freq (Hz)', )
+    axs[2].set_title('(c)', )
+    axs[2].xaxis.grid()
+    axs[2].yaxis.grid()
+
+
+    axs[3].plot((stft[140]/stft.shape[1])**2, 'blue', label='$fft_{3}$', marker="o", markersize=3)
+    axs[3].legend(prop=font)
+    axs[3].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+    axs[3].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+    axs[3].set_xlim([0, 100])
+    axs[3].set_ylabel('Power ($\mu V^{2}$)', )
+    axs[3].set_xlabel('Freq (Hz)', )
+    axs[3].set_title('(d)', )
+    axs[3].xaxis.grid()
+    axs[3].yaxis.grid()
+
+    # plt.title('Frequency Domain ($fft_{1}$, $fft_{2}$, $fft_{3}$)', fontsize=10)
+    plt.tight_layout()
+
+    plt.savefig('fig_4.png')
+
+
+def stft_collections(matrix, matrix_binary, s_rate, stft, cut_freq, task, info_args, max_indx = None, min_indx = None):
+    fig = plt.figure(figsize=(14, 12), dpi=150)
+    grid = plt.GridSpec(6, 8, hspace=0.0, wspace=3.5)
+    spectrogram = fig.add_subplot(grid[0:3, 0:4])
+    rp_plot = fig.add_subplot(grid[0:3, 4:])
+    fft_vector = fig.add_subplot(grid[4:, :])
+
+    if max_indx != None and min_indx != None:
+        max_index = max_indx
+        min_index = min_indx
+    else:
+        max_array = np.max(stft, axis=1)
+        max_value_stft = np.max(max_array, axis=0)
+        max_index =  list(max_array).index(max_value_stft)
+
+        min_array = np.min(stft, axis=1)
+        min_value_stft = np.min(min_array, axis=0)
+        min_index = list(min_array).index(min_value_stft)
+
+
+
+    # ręczne ustawienie wskaźników
+    # max_index = int(1.52*s_rate)
+    # min_index = int(2.4*s_rate)
+
+
+    # top = np.triu(matrix)
+    # bottom = np.tril(matrix_binary)
+
+    # np.linspace(0, stft.shape[1], stft.shape[1]),        np.linspace(0, stft.shape[0], cut_freq),
+    rp_plot.imshow(matrix_binary, cmap='Greys', origin='lower')  # interpolation='none'
+    # axs[0].imshow(bottom, cmap='jet', origin='lower') #interpolation='none'
+    rp_plot.plot(max_index, max_index, 'orange', marker="o", markersize=9)
+    rp_plot.plot(min_index, min_index, 'red', marker="o", markersize=9)
+    # axs[0].set_yticks(axs[0].get_yticks()[1:len(axs[0].get_yticks())-1])
+    # axs[0].set_xticks(axs[0].get_xticks()[1:len(axs[0].get_xticks())-1])
+    rp_plot.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    rp_plot.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    rp_plot.set_xticklabels(
+        [str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, rp_plot.get_xticks().shape[0])])
+    rp_plot.set_yticklabels(
+        [str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, rp_plot.get_yticks().shape[0])])
+    rp_plot.set_xlabel('Time (s)', )
+    rp_plot.set_ylabel('Time (s)', )
+    rp_plot.set_title('(b) Recurrence Plot', )
+    rp_plot.xaxis.grid()
+    rp_plot.yaxis.grid()
+
+    spectrogram.pcolormesh(stft.T,cmap='viridis') #,vmax=max_value_stft
+    spectrogram.plot(max_index,2,'orange', marker="|", markersize=40)
+    spectrogram.plot(min_index,2,'red', marker="|", markersize=40)
+
+    spectrogram.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[0], 5)))
+    spectrogram.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, stft.shape[0] / s_rate, spectrogram.get_xticks().shape[0])])
+    spectrogram.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 5)))
+    spectrogram.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 5)])
+    spectrogram.set_ylabel('Freq (Hz)', )
+    spectrogram.set_xlabel('Time (s)', )
+    spectrogram.set_title('(a) Spectrogram', )
+    # spectrogram.xaxis.grid()
+    # spectrogram.yaxis.grid()
+    # fig.colorbar(im1, cax=spectrogram, orientation='vertical')
+
+
+    max_index_ = stft[max_index]/stft.shape[1]
+    min_index_ = stft[min_index]/stft.shape[1]
+    fft_vector.plot(max_index_**2,'orange',label='$fft_{t_{1}}$')#,marker="o", markersize=2
+    fft_vector.plot(min_index_**2,'red',label='$fft_{t_{2}}}$')#,marker="o", markersize=2
+    fft_vector.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+    fft_vector.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+    fft_vector.set_xlim([0,100])
+    fft_vector.set_ylabel('Power ($\mu V^{2}$)', )
+    fft_vector.set_xlabel('Freq (Hz)', )
+    fft_vector.set_title('(c) Frequency Domain', )
+    fft_vector.legend(prop=font)
+    fft_vector.xaxis.grid()
+    fft_vector.yaxis.grid()
+
+    # plt.suptitle( 'Condition: '+ task  + '\n' + 'epsilon {},  FFT window size {} '.format(
+    #              str(info_args['eps']), str(info_args['win_len'])) + '\n'
+    #              + 'Subject {}, electrode {}, n_fft {}'.format(str(info_args['selected_subject']),str(info_args['electrode_name']),str(info_args['n_fft'])), fontsize=12 ,ha='left',va='top')
+    plt.tight_layout()
+    plt.savefig('fig_5.png')
+    # axs[0].imshow(matrix_binary, cmap='cividis', origin='lower') #interpolation='none'
+    # # axs[0].imshow(bottom, cmap='jet', origin='lower') #interpolation='none'
+    # axs[0].plot(max_index,max_index,'orange',marker="o", markersize=7)
+    # axs[0].plot(min_index,min_index,'red',marker="o", markersize=7)
+    # # axs[0].set_yticks(axs[0].get_yticks()[1:len(axs[0].get_yticks())-1])
+    # # axs[0].set_xticks(axs[0].get_xticks()[1:len(axs[0].get_xticks())-1])
+    # axs[0].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    # axs[0].yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    # axs[0].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs[0].get_xticks().shape[0])])
+    # axs[0].set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs[0].get_yticks().shape[0])])
+    # axs[0].set_xlabel('Time (s)')
+    # axs[0].set_ylabel('Time (s)')
+    # axs[0].set_title('Recurrence Plot', fontsize=12)
+
+def diagnostic(matrix, matrix_binary, s_rate, stft, cut_freq, task, info_args):
+
+    fig, axs = plt.subplots(3,1, figsize=(7,12), gridspec_kw={'height_ratios':[6,2,1]},dpi=150)
+
+    # Set up the axes with gridspec
+
+
+    max_array = np.max(stft, axis=1)
+    max_value_stft = np.max(max_array, axis=0)
+    max_index =  list(max_array).index(max_value_stft)
+
+    min_array = np.min(stft, axis=1)
+    min_value_stft = np.min(min_array, axis=0)
+    min_index = list(min_array).index(min_value_stft)
+
+    # top = np.triu(matrix)
+    # bottom = np.tril(matrix_binary)
+
+    axs[0].imshow(matrix_binary, cmap='cividis', origin='lower') #interpolation='none'
+    # axs[0].imshow(bottom, cmap='jet', origin='lower') #interpolation='none'
+
+    axs[0].plot(max_index,max_index,'orange',marker="o", markersize=7)
+    axs[0].plot(min_index,min_index,'red',marker="o", markersize=7)
+
+    axs[0].set_yticks(axs[0].get_yticks()[1:len(axs[0].get_yticks())-1])
+    axs[0].set_xticks(axs[0].get_xticks()[1:len(axs[0].get_xticks())-1])
+    axs[0].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    axs[0].yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    axs[0].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs[0].get_xticks().shape[0])])
+    axs[0].set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs[0].get_yticks().shape[0])])
+    axs[0].set_xlabel('Time (s)')
+    axs[0].set_ylabel('Time (s)')
+    axs[0].set_title('Recurrence Plot', )
+
+
+    # np.linspace(0, stft.shape[1], stft.shape[1]),        np.linspace(0, stft.shape[0], cut_freq),
+
+
+    axs[1].pcolormesh(stft.T, shading='gouraud') #,vmax=max_value_stft
+    axs[1].plot(max_index,0,'orange', marker="o", markersize=7)
+    axs[1].plot(min_index,0,'red', marker="o", markersize=7)
+    axs[1].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    axs[1].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs[1].get_xticks().shape[0])])
+    axs[1].yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 5)))
+    axs[1].set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 5)])
+    axs[1].set_ylabel('Freq (Hz)')
+    axs[1].set_xlabel('Time (s)')
+    axs[1].set_title('Spectrogram', )
+
+    max_index_ = stft[max_index]/stft.shape[1]
+    min_index_ = stft[min_index]/stft.shape[1]
+    axs[2].plot(max_index_**2,'orange')#,marker="o", markersize=2
+    axs[2].plot(min_index_**2,'red')#,marker="o", markersize=2
+    axs[2].xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+    axs[2].set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+    axs[2].set_xlim([0,100])
+    axs[2].set_ylabel('Power (µV^2)')
+    axs[2].set_xlabel('Freq (Hz)')
+    axs[2].set_title('Frequency Domain',)
+
+    plt.suptitle( 'Condition: '+ task  + '\n' + 'epsilon {},  FFT window size {} '.format(
+                 str(info_args['eps']), str(info_args['win_len'])) + '\n'
+                 + 'Subject {}, electrode {}, n_fft {}'.format(str(info_args['selected_subject']),str(info_args['electrode_name']),str(info_args['n_fft'])),
+                 ha='left',va='top')
+    plt.tight_layout()
+
+def RecurrencePlot(matrix, matrix_binary, s_rate, stft, cut_freq, task, info_args):
+
+    fig, axs = plt.subplots( figsize=(12,12),dpi=200)
+
+    # Set up the axes with gridspec
+
+
+    max_array = np.max(stft, axis=1)
+    max_value_stft = np.max(max_array, axis=0)
+    max_index =  list(max_array).index(max_value_stft)
+
+    min_array = np.min(stft, axis=1)
+    min_value_stft = np.min(min_array, axis=0)
+    min_index = list(min_array).index(min_value_stft)
+
+    # top = np.triu(matrix)
+    # bottom = np.tril(matrix_binary)
+
+    axs.imshow(matrix_binary, cmap='cividis', origin='lower') #interpolation='none'
+    # axs[0].imshow(bottom, cmap='jet', origin='lower') #interpolation='none'
+
+    # axs[0].plot(max_index,max_index,'orange',marker="o", markersize=7)
+    # axs[0].plot(min_index,min_index,'red',marker="o", markersize=7)
+
+    # axs[0].set_yticks(axs[0].get_yticks()[1:len(axs[0].get_yticks())-1])
+    # axs[0].set_xticks(axs[0].get_xticks()[1:len(axs[0].get_xticks())-1])
+    axs.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    axs.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    axs.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs.get_xticks().shape[0])])
+    axs.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, axs.get_yticks().shape[0])])
+    axs.set_xlabel('Time (s)')
+    axs.set_ylabel('Time (s)')
+    axs.set_title('Recurrence Plot')
+
+
+    # np.linspace(0, stft.shape[1], stft.shape[1]),        np.linspace(0, stft.shape[0], cut_freq),
+
+def features_hists(df, features_list, condition, dpi = 200):
+    fig, axs = plt.subplots(len(features_list),figsize=(6, len(features_list)*3), dpi=dpi)
+    abc = ['(a)','(b)','(c)','(d)','(e)','(f)']
+
+    for i,ax in enumerate(axs):
+        sns.histplot(data=df, x=features_list[i], hue=condition, alpha=0.8, element="bars", fill=False, ax=ax, kde=True)
+        ax.containers[1].remove()
+        ax.containers[0].remove()
+        ax.xaxis.grid()
+        ax.yaxis.grid()
+        ax.set_title(abc[i])
+        # plt.grid(b=None)
+
+    plt.autoscale(enable=True, axis='both', tight=None)
+    fig.tight_layout()
+
+def features_per_subjects_violin(df, features_list, condition, dpi = 200):
+    fig, axs = plt.subplots(len(features_list),figsize=(14, len(features_list)*2), dpi=dpi,sharex='col')
+
+    for i,ax in enumerate(axs):
+        sns.violinplot(data=df, x=df.Subject, y=features_list[i], hue=condition, ax=ax, split=True,linewidth=0.2)
+        ax.legend(loc='lower right')
+
+    axs[len(features_list)-1].set_xticklabels(axs[len(features_list)-1].get_xticklabels(), rotation=90)
+    # axs.set_ylim([0,1])
+
+
+    plt.tick_params(axis='x', which='major', labelsize=16)
+    fig.tight_layout()
+
+def umap_on_condition(df,y, title,labels_name,features_list=['TT', 'RR', 'DET', 'LAM', 'L', 'Lentr'], random_state = 70, n_neighbors = 15, min_dist = 0.25, metric = "hamming", df_type=True):
+
+
+    fig, ax1 = plt.subplots(figsize=(8, 8), dpi=150)
+
+    if df_type:
+        stats_data = df
+    else:
+        stats_data = df[features_list].values
+
+    # Preprocess again
+    pipe = make_pipeline(SimpleImputer(strategy="mean"), QuantileTransformer())
+    X = pipe.fit_transform(stats_data.copy())
+
+    # Fit UMAP to processed data
+    manifold = umap.UMAP(random_state=random_state, n_neighbors=n_neighbors, min_dist=min_dist, metric=metric).fit(X, y)
+    # X_reduced_2 = manifold.transform(X)
+    umap.plot.points(manifold, labels=labels_name, ax=ax1, color_key=np.array(
+        [(0, 0.35, 0.73), (1, 0.83, 0)]))  # ,color_key=np.array([(1,0.83,0),(0,0.35,0.73)])
+    ax1.set_title(title)
+
+def umap_side_by_side_plot(df1, df2, features_list=['TT', 'RR', 'DET', 'LAM', 'L', 'Lentr'], random_state = 70, n_neighbors = 15, min_dist = 0.25, metric = "hamming"):
+
+    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,figsize=(16,8),dpi=150)
+
+    stats_data = df1[features_list].values
+    y = df1.Task.values
+
+    # Preprocess again
+    pipe = make_pipeline(SimpleImputer(strategy="mean"), QuantileTransformer())
+    X = pipe.fit_transform(stats_data.copy())
+
+    # Fit UMAP to processed data
+    manifold = umap.UMAP(random_state=random_state, n_neighbors=n_neighbors, min_dist=min_dist, metric=metric).fit(X, y)
+    # X_reduced_2 = manifold.transform(X)
+    umap.plot.points(manifold, labels=y, ax=ax1, color_key=np.array(
+        [(0, 0.35, 0.73), (1, 0.83, 0)]))  # ,color_key=np.array([(1,0.83,0),(0,0.35,0.73)])
+    ax1.set_xlabel('(a) STFT Condition 0 - open eyes, 1 - closed eyes')
+
+    stats_data = df2[features_list].values
+    y = df2.Task.values
+
+    # Preprocess again
+    pipe = make_pipeline(SimpleImputer(strategy="mean"), QuantileTransformer())
+    X = pipe.fit_transform(stats_data.copy())
+
+    # Fit UMAP to processed data
+    manifold = umap.UMAP(random_state=random_state, n_neighbors=n_neighbors, min_dist=min_dist, metric=metric).fit(X, y)
+    # X_reduced_2 = manifold.transform(X)
+    umap.plot.points(manifold, labels=y, ax=ax2, color_key=np.array(
+        [(0, 0.35, 0.73), (1, 0.83, 0)]))  # ,color_key=np.array([(1,0.83,0),(0,0.35,0.73)])
+    ax2.set_xlabel('(b) TDEMB Condition 0 - open eyes, 1 - closed eyes')
+
+    return
+
+def SVM_histogram(df, lin, lin_pred,title):
+    stats_data = df #[features_list].values
+    plt.figure(dpi=150)
+    all_cechy=np.dot(stats_data, lin.coef_.T)
+    df_all=pd.DataFrame({'vectors':all_cechy.ravel(), 'Task':lin_pred})
+
+
+    a = sns.histplot(data=df_all, x='vectors', hue='Task', alpha=0.8, element="bars", fill=False,kde=True, kde_kws={'bw_adjust':0.4},palette=np.array([(0.3,0.85,0),(0.8,0.0,0.44)]))
+    a.containers[1].remove()
+    a.containers[0].remove()
+    # a = sns.kdeplot(data=df_all, x='vectors', hue='Task', alpha=0.8, bw_adjust=0.4,palette=np.array([(0.3,0.85,0),(0.8,0.0,0.44)]))
+    plt.title(title)
+    plt.xlabel('All')
+    plt.grid(b=None)
+    plt.show()
+
+def f_importances(coef, names):
+    imp = coef
+    imp,names = zip(*sorted(zip(imp,names)))
+    plt.figure()
+    plt.barh(range(len(names)), imp, align='center')
+    plt.yticks(range(len(names)), names)
+    plt.show()
+
+def SVM_features_importance(lin):
+
+
+    lebel_ll = np.array([['TT']*int(64)+ ['RR']*int(64)+
+              ['DET']*int(64)+ ['LAM']*int(64)+
+              ['L']*int(64)+ ['L_entr']*int(64)])
+
+    e_long = "Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8	Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8	Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8	Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8	Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8	Af3	Af4	Af7	Af8	Afz	C1	C2	C3	C4	C5	C6	CZ	Cp1	Cp2	Cp3	Cp4	Cp5	Cp6	Cpz	F1	F2	F3	F4	F5	F6	F7	F8	Fc1	Fc2	Fc3	Fc4	Fc5	Fc6	Fcz	Fp1	Fp2	Fpz	Ft7	Ft8	Fz	Iz	O1	O2	OZ	P1	P2	P3	P4	P5	P6	P7	P8	Po3	Po4	Po7	Po8	Poz	Pz	T10	T7	T8	T9	Tp7	Tp8".replace('\t',',').split(",")
+    y_e_long = np.array(np.unique(e_long, return_inverse=True)[1].tolist())
+
+    df = pd.DataFrame({'feature':lebel_ll[0],
+                       'electrode':e_long,
+                       'coef':lin.coef_[0],
+                       })
+    # df = df[(df.coef.values >= 0.15) | (df.coef.values <= -0.15)]
+
+    f_importances(df.coef, df.feature)
+    f_importances(df.coef, df.electrode)
+
+
+    sns.set_theme(style='darkgrid', rc={'figure.dpi': 120},
+                  font_scale=1.7)
+    fig, ax = plt.subplots(figsize=(16, 10))
+    ax.set_title('Weight of features by electrodes')
+    sns.barplot(x='feature', y='coef', data=df, ax=ax,
+                ci=None,
+                hue='electrode')
+    ax.legend(bbox_to_anchor=(1, 1), title='electrode',prop={'size': 7})
+
+##### HIDDEN MARKOV MODEL
+
+def soft_bounds(T,seg):
+
+    # Identify soft boundaries at each step of fitting
+    bounds_anim = []
+    K = seg[0].shape[1]
+    for it in range(1,len(seg)):
+        sb = np.zeros((T,T))
+        for k in range(K-1):
+            p_change = np.diff(seg[it][:,(k+1):].sum(1))
+            sb[1:,1:] += np.outer(p_change, seg[it][1:,k:(k+2)].sum(1))
+        sb = np.maximum(sb,sb.T)
+        sb = sb/np.max(sb)
+        bounds_anim.append(sb)
+    return bounds_anim
+
+
+def fitting_animation(seg,matrix,s_rate,meta_tick,metastate_id, state_width,color_states_matrix):
+
+    bounds_anim = soft_bounds(matrix.shape[0],seg)
+
+    # Plot timepoint-timepoint correation matrix, with boundaries animated on top
+
+    fig = plt.figure(figsize=(18, 12), dpi=300)
+    grid = plt.GridSpec(4, 12, hspace=0.0, wspace=3.5)
+    ax1 = fig.add_subplot(grid[:, 0:4])
+    ax2 = fig.add_subplot(grid[:, 4:8])
+    ax3 = fig.add_subplot(grid[:, 8:])
+
+
+    # fig, axs = plt.subplots(2,figsize=(8,8), dpi=120)
+    datamat = matrix # np.corrcoef(D)
+    bk = cm.viridis((datamat-np.min(datamat))/(np.max(datamat)-np.min(datamat)))
+    im = ax1.imshow(bk, interpolation='none',origin='lower')
+    fg = cm.gray(1-(sum(bounds_anim)/len(bounds_anim)))
+    # im.set_array(np.minimum(np.maximum(bk + fg, 0), 1))
+    im.set_array(bk * fg)
+
+
+    ax1.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax1.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax1.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax1.get_xticks().shape[0])])
+    ax1.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax1.get_yticks().shape[0])])
+    ax1.set_xlabel('Time (s)')
+    ax1.set_ylabel('Time (s)')
+    ax1.set_title('Metastates plot over recurrence plot', fontsize=10)
+    ax1.scatter(meta_tick,meta_tick,s=2)
+
+
+
+    ax2.imshow(fg, interpolation='none',origin='lower')
+    ax2.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax2.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax2.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax2.get_xticks().shape[0])])
+    ax2.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax2.get_yticks().shape[0])])
+    ax2.set_xlabel('Time (s)')
+    ax2.set_ylabel('Time (s)')
+    ax2.set_title('Metastates plot', fontsize=10)
+    ax2.scatter(meta_tick,meta_tick,s=2)
+
+    text_kwargs = dict(ha='center', va='center', fontsize=4, color='0')
+    for i,mstate in enumerate(metastate_id):
+        # ax1.text(meta_tick[i]-35, meta_tick[i]+(state_width[i]/2)+45, 's'+str(mstate)+'| '+ str(int(((1/160)*state_width[i])*1000)) + 'ms', **text_kwargs)
+        ax2.annotate('s '+str(mstate)+'| '+ str(int(((1/160)*state_width[i])*1000)) + 'ms', xy=(meta_tick[i], meta_tick[i]+(state_width[i]/2)),
+                     xytext =(meta_tick[i], meta_tick[i]+(state_width[i]/2)+70),
+                     xycoords='data',
+                     textcoords='data',
+                     arrowprops=dict(arrowstyle="->",facecolor='blue'),
+                     horizontalalignment='right', verticalalignment='top', fontsize=5
+                     )
+
+
+    color_states = color_states_matrix
+    ax3.imshow(fg[:,:,:3]*color_states, interpolation='none',origin='lower')
+    ax3.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax3.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+    ax3.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax3.get_xticks().shape[0])])
+    ax3.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, ax3.get_yticks().shape[0])])
+    ax3.set_xlabel('Time (s)')
+    ax3.set_ylabel('Time (s)')
+    ax3.set_title('Metastates plot', fontsize=10)
+    ax3.scatter(meta_tick,meta_tick,s=2)
+
+    text_kwargs = dict(ha='center', va='center', fontsize=4, color='0')
+    for i,mstate in enumerate(metastate_id):
+        # ax1.text(meta_tick[i]-35, meta_tick[i]+(state_width[i]/2)+45, 's'+str(mstate)+'| '+ str(int(((1/160)*state_width[i])*1000)) + 'ms', **text_kwargs)
+        ax3.annotate('s '+str(mstate)+'| '+ str(int(((1/160)*state_width[i])*1000)) + 'ms', xy=(meta_tick[i], meta_tick[i]+(state_width[i]/2)),
+                     xytext =(meta_tick[i], meta_tick[i]+(state_width[i]/2)+70),
+                     xycoords='data',
+                     textcoords='data',
+                     arrowprops=dict(arrowstyle="->",facecolor='blue'),
+                     horizontalalignment='right', verticalalignment='top', fontsize=5
+                     )
+
+
+    # def animate_func(i):
+    #     fg = cm.Greys(1-bounds_anim[i])
+    #     im.set_array(np.minimum(np.maximum(bk + fg,0),1))
+    #     return [im]
+    #
+    # anim = animation.FuncAnimation(fig, animate_func,
+    #                                frames = len(bounds_anim), interval = 1)
+    #
+    #
+    plt.savefig('Metastate.png')
+    # plt.close("all")
+
+    return fig
+
+    # return HTML(anim.to_jshtml(default_mode='Once'))
+
+def fit_HMM(matrix,n_events):
+    return EventSegment(n_events=n_events).fit(matrix)
+
+def metastates(seg,matrix,s_rate,meta_tick,metastate_id, state_width,color_states_matrix):
+    fitting_animation(seg,matrix,s_rate,meta_tick,metastate_id, state_width,color_states_matrix)
+
+# def diagnostic(matrix, matrix_binary, s_rate, stft, cut_freq, task, info_args):
+#
+#     # fig, axs = plt.subplots(3,1, figsize=(4,8), gridspec_kw={'height_ratios':[6,2,1]},dpi=120)
+#
+#     # Set up the axes with gridspec
+#     fig = plt.figure(figsize=(6, 6),dpi=120)
+#     grid = plt.GridSpec(6, 6, hspace=1.0, wspace=1.0)
+#     spectrogram = fig.add_subplot(grid[0:3, 0:3])
+#     rp_plot = fig.add_subplot(grid[0:3, 3:])
+#     fft_vector = fig.add_subplot(grid[3:,:])
+#
+#     max_array = np.max(stft, axis=1)
+#     max_value_stft = np.max(max_array, axis=0)
+#     max_index =  list(max_array).index(max_value_stft)
+#
+#     min_array = np.min(stft, axis=1)
+#     min_value_stft = np.min(min_array, axis=0)
+#     min_index = list(min_array).index(min_value_stft)
+#
+#     # top = np.triu(matrix)
+#     # bottom = np.tril(matrix_binary)
+#
+#     rp_plot.imshow(matrix_binary, cmap='cividis', origin='lower') #interpolation='none'
+#     # axs[0].imshow(bottom, cmap='jet', origin='lower') #interpolation='none'
+#     rp_plot.plot(max_index,max_index,'orange',marker="o", markersize=7)
+#     rp_plot.plot(min_index,min_index,'red',marker="o", markersize=7)
+#     # axs[0].set_yticks(axs[0].get_yticks()[1:len(axs[0].get_yticks())-1])
+#     # axs[0].set_xticks(axs[0].get_xticks()[1:len(axs[0].get_xticks())-1])
+#     rp_plot.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+#     rp_plot.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+#     rp_plot.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, rp_plot.get_xticks().shape[0])])
+#     rp_plot.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, rp_plot.get_yticks().shape[0])])
+#     rp_plot.set_xlabel('Time (s)')
+#     rp_plot.set_ylabel('Time (s)')
+#     rp_plot.set_title('Recurrence Plot', fontsize=10)
+#
+#
+#     # np.linspace(0, stft.shape[1], stft.shape[1]),        np.linspace(0, stft.shape[0], cut_freq),
+#
+#
+#     spectrogram.pcolormesh(stft.T, shading='gouraud') #,vmax=max_value_stft
+#     spectrogram.plot(max_index,0,'orange', marker="o", markersize=7)
+#     spectrogram.plot(min_index,0,'red', marker="o", markersize=7)
+#     spectrogram.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, matrix.shape[0], 5)))
+#     spectrogram.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, matrix.shape[0] / s_rate, spectrogram.get_xticks().shape[0])])
+#     spectrogram.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 5)))
+#     spectrogram.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 5)])
+#     spectrogram.set_ylabel('Freq (Hz)')
+#     spectrogram.set_xlabel('Time (s)')
+#     spectrogram.set_title('Spectrogram', fontsize=10)
+#
+#     max_index_ = stft[max_index]/stft.shape[1]
+#     min_index_ = stft[min_index]/stft.shape[1]
+#     fft_vector.plot(max_index_**2,'orange')#,marker="o", markersize=2
+#     fft_vector.plot(min_index_**2,'red')#,marker="o", markersize=2
+#     fft_vector.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft.shape[1], 9)))
+#     fft_vector.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 9)])
+#     fft_vector.set_xlim([0,100])
+#     fft_vector.set_ylabel('Power (µV^2)')
+#     fft_vector.set_xlabel('Freq (Hz)')
+#     fft_vector.set_title('Frequency Domain', size=10)
+#
+#     plt.suptitle( 'Condition: '+ task  + '\n' + 'epsilon {},  FFT window size {} '.format(
+#                  str(info_args['eps']), str(info_args['win_len'])) + '\n'
+#                  + 'Subject {}, electrode {}, n_fft {}'.format(str(info_args['selected_subject']),str(info_args['electrode_name']),str(info_args['n_fft'])),
+#                  fontsize=8,ha='left',va='top')
+#     plt.tight_layout()
+

BrainPulse/recurrence_quantification_analysis.py

@@ -0,0 +1,305 @@
+from numba import jit
+import numpy as np
+from pyrqa.computation import RPComputation, RQAComputation
+from pyrqa.time_series import TimeSeries, EmbeddedSeries
+from pyrqa.settings import Settings
+from pyrqa.analysis_type import Classic
+from pyrqa.metric import EuclideanMetric
+# from pyrqa.metric import Sigmoid
+# from pyrqa.metric import Cosine
+from pyrqa.neighbourhood import Unthresholded,FixedRadius
+
+def get_results(recurrence_matrix,
+                minimum_diagonal_line_length,
+                minimum_vertical_line_length,
+                minimum_white_vertical_line_length):
+
+    number_of_vectors = recurrence_matrix.shape[0]
+    diagonal = diagonal_frequency_distribution(recurrence_matrix)
+    vertical = vertical_frequency_distribution(recurrence_matrix)
+    white = white_vertical_frequency_distribution(recurrence_matrix)
+
+    number_of_vert_lines = number_of_vertical_lines(vertical, minimum_vertical_line_length)
+    number_of_vert_lines_points = number_of_vertical_lines_points(vertical, minimum_vertical_line_length)
+
+    RR = recurrence_rate(recurrence_matrix)
+    DET = determinism(number_of_vectors, diagonal, minimum_diagonal_line_length)
+    L = average_diagonal_line_length(number_of_vectors, diagonal, minimum_diagonal_line_length)
+    Lmax = longest_diagonal_line_length(number_of_vectors, diagonal)
+    DIV = divergence(Lmax)
+    Lentr = entropy_diagonal_lines(number_of_vectors, diagonal, minimum_diagonal_line_length)
+    DET_RR = ratio_determinism_recurrence_rate(DET, RR)
+    LAM = laminarity(number_of_vectors, vertical, minimum_vertical_line_length)
+    V = average_vertical_line_length(number_of_vectors, vertical, minimum_vertical_line_length)
+    Vmax = longest_vertical_line_length(number_of_vectors, vertical)
+    Ventr = entropy_vertical_lines(number_of_vectors, vertical, minimum_vertical_line_length)
+    LAM_DET = laminarity_determinism(LAM, DET)
+    W = average_white_vertical_line_length(number_of_vectors, white, minimum_white_vertical_line_length)
+    Wmax = longest_white_vertical_line_length(number_of_vectors, white)
+    Wentr = entropy_white_vertical_lines(number_of_vectors, white, minimum_white_vertical_line_length)
+    TT = trapping_time(number_of_vert_lines_points, number_of_vert_lines)
+
+    return [RR, DET, L, Lmax, DIV, Lentr, DET_RR, LAM, V, Vmax, Ventr, LAM_DET, W, Wmax, Wentr, TT]
+
+
+@jit(nopython=True)
+def diagonal_frequency_distribution(recurrence_matrix):
+    # Calculating the number of states - N
+    number_of_vectors = recurrence_matrix.shape[0]
+    diagonal_frequency_distribution = np.zeros(number_of_vectors + 1)
+
+    # Calculating the diagonal frequency distribution - P(l)
+    for i in range(number_of_vectors - 1, -1, -1):
+        diagonal_line_length = 0
+        for j in range(0, number_of_vectors - i):
+            if recurrence_matrix[i + j, j] == 1:
+                diagonal_line_length += 1
+                if j == (number_of_vectors - i - 1):
+                    diagonal_frequency_distribution[diagonal_line_length] += 1.0
+            else:
+                if diagonal_line_length != 0:
+                    diagonal_frequency_distribution[diagonal_line_length] += 1.0
+                    diagonal_line_length = 0
+    for k in range(1, number_of_vectors):
+        diagonal_line_length = 0
+        for i in range(number_of_vectors - k):
+            j = i + k
+            if recurrence_matrix[i, j] == 1:
+                diagonal_line_length += 1
+                if j == (number_of_vectors - 1):
+                    diagonal_frequency_distribution[diagonal_line_length] += 1.0
+            else:
+                if diagonal_line_length != 0:
+                    diagonal_frequency_distribution[diagonal_line_length] += 1.0
+                    diagonal_line_length = 0
+
+    return diagonal_frequency_distribution
+
+
+@jit(nopython=True)
+def vertical_frequency_distribution(recurrence_matrix):
+    number_of_vectors = recurrence_matrix.shape[0]
+
+    # Calculating the vertical frequency distribution - P(v)
+    vertical_frequency_distribution = np.zeros(number_of_vectors + 1)
+    for i in range(number_of_vectors):
+        vertical_line_length = 0
+        for j in range(number_of_vectors):
+            if recurrence_matrix[i, j] == 1:
+                vertical_line_length += 1
+                if j == (number_of_vectors - 1):
+                    vertical_frequency_distribution[vertical_line_length] += 1.0
+            else:
+                if vertical_line_length != 0:
+                    vertical_frequency_distribution[vertical_line_length] += 1.0
+                    vertical_line_length = 0
+
+    return vertical_frequency_distribution
+
+
+@jit(nopython=True)
+def white_vertical_frequency_distribution(recurrence_matrix):
+    number_of_vectors = recurrence_matrix.shape[0]
+
+    # Calculating the white vertical frequency distribution - P(w)
+    white_vertical_frequency_distribution = np.zeros(number_of_vectors + 1)
+    for i in range(number_of_vectors):
+        white_vertical_line_length = 0
+        for j in range(number_of_vectors):
+            if recurrence_matrix[i, j] == 0:
+                white_vertical_line_length += 1
+                if j == (number_of_vectors - 1):
+                    white_vertical_frequency_distribution[white_vertical_line_length] += 1.0
+            else:
+                if white_vertical_line_length != 0:
+                    white_vertical_frequency_distribution[white_vertical_line_length] += 1.0
+                    white_vertical_line_length = 0
+
+    return white_vertical_frequency_distribution
+
+
+@jit(nopython=True)
+def recurrence_rate(recurrence_matrix):
+    # Calculating the recurrence rate - RR
+    number_of_vectors = recurrence_matrix.shape[0]
+    return np.float(np.sum(recurrence_matrix)) / np.power(number_of_vectors, 2)
+
+
+def determinism(number_of_vectors, diagonal_frequency_distribution_, minimum_diagonal_line_length):
+    # Calculating the determinism - DET
+    numerator = np.sum(
+        [l * diagonal_frequency_distribution_[l] for l in range(minimum_diagonal_line_length, number_of_vectors)])
+    denominator = np.sum([l * diagonal_frequency_distribution_[l] for l in range(1, number_of_vectors)])
+    return numerator / denominator
+
+
+def average_diagonal_line_length(number_of_vectors, diagonal_frequency_distribution_, minimum_diagonal_line_length):
+    # Calculating the average diagonal line length - L
+    numerator = np.sum(
+        [l * diagonal_frequency_distribution_[l] for l in range(minimum_diagonal_line_length, number_of_vectors)])
+    denominator = np.sum(
+        [diagonal_frequency_distribution_[l] for l in range(minimum_diagonal_line_length, number_of_vectors)])
+    return numerator / denominator
+
+
+@jit(nopython=True)
+def longest_diagonal_line_length(number_of_vectors, diagonal_frequency_distribution_):
+    # Calculating the longest diagonal line length - Lmax
+    for l in range(number_of_vectors - 1, 0, -1):
+        if diagonal_frequency_distribution_[l] != 0:
+            longest_diagonal_line_length = l
+            break
+    return longest_diagonal_line_length
+
+
+@jit(nopython=True)
+def divergence(longest_diagonal_line_length_):
+    # Calculating the  divergence - DIV
+    return 1. / longest_diagonal_line_length_
+
+
+@jit(nopython=True)
+def entropy_diagonal_lines(number_of_vectors, diagonal_frequency_distribution_, minimum_diagonal_line_length):
+    # Calculating the entropy diagonal lines - Lentr
+    sum_diagonal_frequency_distribution = np.float(
+        np.sum(diagonal_frequency_distribution_[minimum_diagonal_line_length:-1]))
+    entropy_diagonal_lines = 0
+    for l in range(minimum_diagonal_line_length, number_of_vectors):
+        if diagonal_frequency_distribution_[l] != 0:
+            entropy_diagonal_lines += (diagonal_frequency_distribution_[
+                                           l] / sum_diagonal_frequency_distribution) * np.log(
+                diagonal_frequency_distribution_[l] / sum_diagonal_frequency_distribution)
+    entropy_diagonal_lines *= -1
+    return entropy_diagonal_lines
+
+
+@jit(nopython=True)
+def ratio_determinism_recurrence_rate(determinism_, recurrence_rate_):
+    # Calculating the  divergence - DIV
+    return determinism_ / recurrence_rate_
+
+
+def laminarity(number_of_vectors, vertical_frequency_distribution_, minimum_vertical_line_length):
+    # Calculating the laminarity - LAM
+    numerator = np.sum(
+        [v * vertical_frequency_distribution_[v] for v in range(minimum_vertical_line_length, number_of_vectors + 1)])
+    denominator = np.sum([v * vertical_frequency_distribution_[v] for v in range(1, number_of_vectors + 1)])
+    return numerator / denominator
+
+
+def average_vertical_line_length(number_of_vectors, vertical_frequency_distribution_, minimum_vertical_line_length):
+    # Calculating the average vertical line length - V
+    numerator = np.sum(
+        [v * vertical_frequency_distribution_[v] for v in range(minimum_vertical_line_length, number_of_vectors + 1)])
+    denominator = np.sum(
+        [vertical_frequency_distribution_[v] for v in range(minimum_vertical_line_length, number_of_vectors + 1)])
+    return numerator / denominator
+
+
+@jit(nopython=True)
+def longest_vertical_line_length(number_of_vectors, vertical_frequency_distribution_):
+    # Calculating the longest vertical line length - Vmax
+    longest_vertical_line_length_ = 0
+    for v in range(number_of_vectors, 0, -1):
+        if vertical_frequency_distribution_[v] != 0:
+            longest_vertical_line_length_ = v
+            break
+    return longest_vertical_line_length_
+
+
+@jit(nopython=True)
+def entropy_vertical_lines(number_of_vectors, vertical_frequency_distribution_, minimum_vertical_line_length):
+    # Calculating the entropy vertical lines - Ventr
+    sum_vertical_frequency_distribution = np.float(
+        np.sum(vertical_frequency_distribution_[minimum_vertical_line_length:]))
+    entropy_vertical_lines_ = 0
+    for v in range(minimum_vertical_line_length, number_of_vectors + 1):
+        if vertical_frequency_distribution_[v] != 0:
+            entropy_vertical_lines_ += (vertical_frequency_distribution_[
+                                           v] / sum_vertical_frequency_distribution) * np.log(
+                vertical_frequency_distribution_[v] / sum_vertical_frequency_distribution)
+    entropy_vertical_lines_ *= -1
+    return entropy_vertical_lines_
+
+
+@jit(nopython=True)
+def laminarity_determinism(laminarity_, determinism_):
+    # Calculating the ratio laminarity_determinism - LAM/DET
+    return laminarity_ / determinism_
+
+
+def average_white_vertical_line_length(number_of_vectors, white_vertical_frequency_distribution_,
+                                       minimum_white_vertical_line_length):
+    # Calculating the average white vertical line length - W
+    numerator = np.sum([w * white_vertical_frequency_distribution_[w] for w in
+                        range(minimum_white_vertical_line_length, number_of_vectors + 1)])
+    denominator = np.sum([white_vertical_frequency_distribution_[w] for w in
+                          range(minimum_white_vertical_line_length, number_of_vectors + 1)])
+    return numerator / denominator
+
+
+@jit(nopython=True)
+def longest_white_vertical_line_length(number_of_vectors, white_vertical_frequency_distribution_):
+    # Calculating the longest white vertical line length - Wmax
+    longest_white_vertical_line_length_ = 0
+    for w in range(number_of_vectors, 0, -1):
+        if white_vertical_frequency_distribution_[w] != 0:
+            longest_white_vertical_line_length_ = w
+            break
+    return longest_white_vertical_line_length_
+
+
+@jit(nopython=True)
+def entropy_white_vertical_lines(number_of_vectors, white_vertical_frequency_distribution_,
+                                 minimum_white_vertical_line_length):
+    # Calculating the entropy white vertical lines - Wentr
+    sum_white_vertical_frequency_distribution = np.float(
+        np.sum(white_vertical_frequency_distribution_[minimum_white_vertical_line_length:]))
+    entropy_white_vertical_lines_ = 0
+    for w in range(minimum_white_vertical_line_length, number_of_vectors + 1):
+        if white_vertical_frequency_distribution_[w] != 0:
+            entropy_white_vertical_lines_ += (white_vertical_frequency_distribution_[
+                                                  w] / sum_white_vertical_frequency_distribution) * np.log(
+                white_vertical_frequency_distribution_[w] / sum_white_vertical_frequency_distribution)
+    entropy_white_vertical_lines_ *= -1
+    return entropy_white_vertical_lines_
+
+def number_of_vertical_lines(vertical_frequency_distribution_, minimum_vertical_line_length):
+    if minimum_vertical_line_length > 0:
+        return np.sum(vertical_frequency_distribution_[minimum_vertical_line_length - 1:])
+
+    return np.uint(0)
+
+
+def number_of_vertical_lines_points(vertical_frequency_distribution_, minimum_vertical_line_length):
+    if minimum_vertical_line_length > 0:
+        return np.sum(
+            ((np.arange(vertical_frequency_distribution_.size) + 1) * vertical_frequency_distribution_)[minimum_vertical_line_length - 1:])
+
+    return np.uint(0)
+
+@jit(nopython=True)
+def trapping_time(number_of_vertical_lines_points_, number_of_vertical_lines_):
+    """
+    Trapping time (TT).
+    """
+    try:
+        return np.float32(number_of_vertical_lines_points_ / number_of_vertical_lines_)
+    except:
+        return 0
+
+
+
+
+
+def return_pyRQA_results(signal, nbr):
+    time_series = EmbeddedSeries(signal)
+    settings = Settings(time_series,
+                        analysis_type=Classic,
+                        neighbourhood=FixedRadius(nbr),
+                        similarity_measure=EuclideanMetric,
+                        theiler_corrector=1)
+    computation = RQAComputation.create(settings,
+                                        verbose=True)
+    result = computation.run()
+    return result

BrainPulse/vector_space.py

@@ -0,0 +1,16 @@
+import time
+import torch
+
+def compute_stft(epoch_, n_fft, win_len, s_rate, cut_freq):
+    # stft_time = time.time()
+
+    signal_tensor = torch.tensor(epoch_, dtype=torch.float)
+    stft_tensor = torch.stft(signal_tensor,n_fft=n_fft, win_length=win_len, hop_length=1,return_complex=True,window=torch.hann_window(win_len))
+
+    sft = torch.abs(stft_tensor).numpy()
+
+    freq_to_take = (((n_fft/2)+1)*cut_freq) / ((s_rate/2)+1)
+
+    sft = sft[:int(freq_to_take),::]
+
+    return sft.T

Dockerfile

@@ -0,0 +1,21 @@
+# syntax=docker/dockerfile:1
+
+FROM python:3.8-slim-buster
+RUN apt-get update && apt-get install -y git
+
+COPY ./requirements.txt /app/requirements.txt
+
+WORKDIR /app
+
+RUN pip3 install --no-cache-dir -r requirements.txt
+
+EXPOSE 5000
+
+COPY . /app
+
+# ENTRYPOINT [ "python" ]
+# ENTRYPOINT [ "python" ]
+# ENTRYPOINT [ "streamlit run BrainPulseAPP.py" ]
+
+# CMD [ "BrainPulseAPP.py" ]
+CMD ["streamlit", "run", "BrainPulseAPP.py"]

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2022 Łukasz Furman
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -1,13 +1,2 @@
----
-title: NCU
-emoji: 📚
-colorFrom: indigo
-colorTo: yellow
-sdk: streamlit
-sdk_version: 1.10.0
-app_file: app.py
-pinned: false
-license: mit
----
-
-Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference
+# BrainPulse_webapp
+ 

app.py

@@ -0,0 +1,260 @@
+from cProfile import run
+import streamlit as st
+import matplotlib.pyplot as plt
+import matplotlib
+import numpy as np
+import pandas as pd
+import plotly.graph_objects as go
+from complexRadar import ComplexRadar
+import math
+from zipfile import ZipFile
+from glob import glob
+import os
+from BrainPulse import (dataset,
+                        vector_space,
+                        distance_matrix,
+                        recurrence_quantification_analysis,
+                        features_space,
+                        plot)
+
+# path
+path = "./mne_data"
+path2 = "./RPs"
+
+# Remove the specified
+# file path
+try:
+    os.remove(path)
+    print("% s removed successfully" % path)
+except:
+    pass
+
+path = "./mne_data"
+os.makedirs(path, exist_ok = True)
+path1 = "./RPs"
+os.makedirs(path1, exist_ok = True)
+
+def run_computation(t_start, t_end, selected_subject, fir_filter, electrode_name, cut_freq, win_len, n_fft, percentile, run_list, options):
+    
+    epochs, raw = dataset.eegbci_data(tmin=t_start, tmax=t_end,
+                             subject=selected_subject,
+                             filter_range=fir_filter,run_list=run_list)
+
+    s_rate = epochs.info['sfreq']
+
+    electrode_index = epochs.ch_names.index(electrode_name)
+
+    electrode_open = epochs.get_data()[0][electrode_index]
+    electrode_close = epochs.get_data()[1][electrode_index]
+
+    stft_open = vector_space.compute_stft(electrode_open,
+                                        n_fft=n_fft, win_len=win_len,
+                                        s_rate=epochs.info['sfreq'],
+                                        cut_freq=cut_freq)
+
+    stft_close = vector_space.compute_stft(electrode_close,
+                                        n_fft=n_fft, win_len=win_len,
+                                        s_rate=epochs.info['sfreq'],
+                                        cut_freq=cut_freq)
+    del raw
+    del electrode_open, electrode_close
+    # matrix_open = distance_matrix.EuclideanPyRQA_RP_stft(stft_open)
+    # matrix_close = distance_matrix.EuclideanPyRQA_RP_stft(stft_close)
+    matrix_open = distance_matrix.EuclideanPyRQA_RP_stft_cpu(stft_open)
+    matrix_close = distance_matrix.EuclideanPyRQA_RP_stft_cpu(stft_close)
+
+    nbr_open = np.percentile(matrix_open, percentile)
+    nbr_close = np.percentile(matrix_close, percentile)
+
+    matrix_open_binary = distance_matrix.set_epsilon(matrix_open,nbr_open)
+    matrix_close_binary = distance_matrix.set_epsilon(matrix_close,nbr_close)
+
+    del matrix_open, matrix_close
+    # matrix_open_to_plot = matrix_open_binary
+    # matrix_closed_to_plot = matrix_close_binary
+
+    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,figsize=(16,8),dpi=200)
+    ax1.imshow(matrix_open_binary, cmap='Greys', origin='lower') #cividis
+    ax1.set_xticks(np.linspace(0, matrix_open_binary.shape[0] , ax1.get_xticks().shape[0]))
+    ax1.set_yticks(np.linspace(0, matrix_open_binary.shape[0] , ax1.get_xticks().shape[0]))
+    ax1.set_xticklabels([str(np.around(x,decimals=0)) for x in np.linspace(0, matrix_open_binary.shape[0] / s_rate, ax1.get_xticks().shape[0])])
+    ax1.set_yticklabels([str(np.around(x, decimals=0)) for x in np.linspace(0, matrix_open_binary.shape[0] / s_rate, ax1.get_yticks().shape[0])])
+    ax1.set_title(options[0]+' window size = 240 samples, ε = '+str(np.round(nbr_open,4)))
+    ax1.set_xlabel('time (s)')
+    ax1.set_ylabel('time (s)')
+
+    ax2.imshow(matrix_close_binary, cmap='Greys', origin='lower')
+    ax2.set_xticks(np.linspace(0, matrix_close_binary.shape[0] , ax1.get_xticks().shape[0]))
+    ax2.set_yticks(np.linspace(0, matrix_close_binary.shape[0] , ax1.get_xticks().shape[0]))
+    ax2.set_xticklabels([str(np.around(x,decimals=0)) for x in np.linspace(0, matrix_close_binary.shape[0] / s_rate, ax1.get_xticks().shape[0])])
+    ax2.set_yticklabels([str(np.around(x, decimals=0)) for x in np.linspace(0, matrix_close_binary.shape[0] / s_rate, ax2.get_yticks().shape[0])])
+    ax2.set_title(options[1]+' window size = 240 samples, ε = '+str(np.round(nbr_close,4)))
+    ax2.set_xlabel('time (s)')
+    ax2.set_ylabel('time (s)')
+
+    return fig, matrix_open_binary, matrix_close_binary, epochs, stft_open, stft_close
+
+
+def plot_rqa(matrix_open_binary, matrix_close_binary, min_vert_line_len, min_diagonal_line_len, min_white_vert_line_len,options):
+
+    categories = ['RR', 'DET', 'L', 'Lmax', 'DIV', 'Lentr', 'DET_RR', 'LAM', 'V', 'Vmax', 'Ventr', 'LAM_DET', 'W', 'Wmax', 'Wentr', 'TT']
+
+    result_rqa_open = recurrence_quantification_analysis.get_results(matrix_open_binary,min_vert_line_len, min_diagonal_line_len, min_white_vert_line_len)
+    result_rqa_closed = recurrence_quantification_analysis.get_results(matrix_close_binary,min_vert_line_len, min_diagonal_line_len, min_white_vert_line_len) 
+
+    data = pd.DataFrame([result_rqa_open,result_rqa_closed], columns=categories)
+    
+    data = data.drop(['RR', 'DIV', 'Lmax'],axis=1)
+    # print(data)
+    min_max_per_variable = data.describe().T[['min', 'max']]
+    min_max_per_variable['min'] = min_max_per_variable['min'].apply(lambda x: int(x))
+    min_max_per_variable['max'] = min_max_per_variable['max'].apply(lambda x: math.ceil(x))
+    # print(min_max_per_variable)
+
+
+    
+    variables = data.columns
+    ranges = list(min_max_per_variable.itertuples(index=False, name=None))   
+
+    format_cfg = {
+        #'axes_args':{'facecolor':'#84A8CD'},
+        'rad_ln_args': {'visible':True, 'linestyle':'dotted'},
+        'angle_ln_args':{'linestyle':'dotted'},
+        'outer_ring': {'visible':True, 'linestyle':'dotted'},
+        'rgrid_tick_lbls_args': {'fontsize':6},
+        'theta_tick_lbls': {'fontsize':9, 'backgroundcolor':'#355C7D', 'color':'#FFFFFF'},
+        'theta_tick_lbls_pad':3
+    }
+
+
+    fig = plt.figure(figsize=(5,3),dpi=100)
+    radar = ComplexRadar(fig, variables, ranges,n_ring_levels=3 ,show_scales=True, format_cfg=format_cfg)
+
+
+    custom_colors = ['#F67280', '#6C5B7B', '#355C7D']
+    k=0
+    for g,c in zip(data.index, custom_colors):
+        # radar.plot(data.loc[g].values, label=f"condition {g}", color=c, marker='o')
+        radar.plot(data.loc[g].values, label=options[k], color=c, marker='o')
+        radar.fill(data.loc[g].values, alpha=0.5, color=c)
+        k+=1
+
+    radar.use_legend(loc='upper left', bbox_to_anchor=(-0.4, 1.1), fontsize = 'xx-small') #, bbox_to_anchor=(0.15, -0.25),ncol=radar.plot_counter
+
+    return fig
+
+def waterfall_spectrum(stft1, stft2, s_rate, cut_freq, options):
+
+    fig = plt.figure(figsize=(14, 12), dpi=150)
+    grid = plt.GridSpec(8, 8, hspace=0.0, wspace=3.5)
+    spectrogram1 = fig.add_subplot(grid[0:3, 0:4])
+    spectrogram2 = fig.add_subplot(grid[0:3, 4:])
+
+    spectrogram1.pcolormesh(stft1.T,cmap='viridis')
+    spectrogram1.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft1.shape[0], 5)))
+    spectrogram1.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, stft1.shape[0] / s_rate, spectrogram1.get_xticks().shape[0])])
+    spectrogram1.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft1.shape[1], 5)))
+    spectrogram1.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 5)])
+    spectrogram1.set_ylabel('Freq (Hz)', )
+    spectrogram1.set_xlabel('Time (s)', )
+    spectrogram1.set_title(options[0] + ' Spectrogram', )
+    
+    spectrogram2.pcolormesh(stft2.T,cmap='viridis')
+    spectrogram2.xaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft2.shape[0], 5)))
+    spectrogram2.set_xticklabels([str(np.round(x, 1)) for x in np.linspace(0, stft2.shape[0] / s_rate, spectrogram2.get_xticks().shape[0])])
+    spectrogram2.yaxis.set_major_locator(matplotlib.ticker.FixedLocator(np.linspace(0, stft2.shape[1], 5)))
+    spectrogram2.set_yticklabels([str(np.round(x, 1)) for x in np.linspace(0, cut_freq, 5)])
+    spectrogram2.set_ylabel('Freq (Hz)', )
+    spectrogram2.set_xlabel('Time (s)', )
+    spectrogram2.set_title(options[1] +' Spectrogram', )
+    return fig
+
+def save(matrix_open_binary, matrix_close_binary):
+    
+    file_name_open = './RPs/subject-'+str(selected_subject)+'_electrode-'+electrode_name+'_percentile-'+str(percentile)+'_run-open_binary.npy'
+    np.save(file_name_open, np.asarray(matrix_close_binary, dtype=np.ubyte))
+    file_name_close = './RPs/subject-'+str(selected_subject)+'_electrode-'+electrode_name+'_percentile-'+str(percentile)+'_run-close_binary.npy'
+    np.save(file_name_close, np.asarray(matrix_close_binary, dtype=np.ubyte))
+
+def download():
+
+    file_paths = glob('./RPs/*')
+
+    with ZipFile('download.zip','w') as zip:
+        for file in file_paths:
+            # writing each file one by one
+            zip.write(file)
+
+    return open('download.zip', 'rb')
+
+# ---------------Settings--------------------
+
+st.set_page_config(layout="wide")
+st.title('BrainPulse Playground')
+sidebar = st.sidebar
+
+selected_subject = sidebar.slider('Select Subject', 0, 100, 25)
+
+electrode_name = sidebar.selectbox(
+    'Select Electrode',
+    ('FC5', 'FC3', 'FC1', 'FCz', 'FC2', 'FC4', 'FC6', 'C5', 'C3', 'C1', 'Cz', 'C2', 'C4', 'C6', 'CP5', 'CP3', 'CP1', 'CPz', 'CP2', 'CP4', 'CP6', 'Fp1', 'Fpz', 'Fp2', 'AF7', 'AF3', 'AFz', 'AF4', 'AF8', 'F7', 'F5', 'F3', 'F1', 'Fz', 'F2', 'F4', 'F6', 'F8', 'FT7', 'FT8', 'T7', 'T8', 'T9', 'T10', 'TP7', 'TP8', 'P7', 'P5', 'P3', 'P1', 'Pz', 'P2', 'P4', 'P6', 'P8', 'PO7', 'PO3', 'POz', 'PO4', 'PO8', 'O1', 'Oz', 'O2', 'Iz'))
+
+t_start, t_end = sidebar.slider(
+    'Select a time range in seconds',
+    0, 60, (0, 30), step=1)
+
+f1, f2 = sidebar.slider(
+    'Select a FIR filter range',
+    0, 60, (2, 50), step=1)
+fir_filter = [f1, f2]
+
+cut_freq = f2
+
+win_len = sidebar.slider('FFT window size', 0, 512, 170, step=1)
+
+n_fft = sidebar.slider('numer of FFT bins', 0, 1024, 512, step=1)
+
+min_vert_line_len = sidebar.slider('Minimum  vertical line length', 0, 250, 2, step=1)
+
+min_diagonal_line_len = sidebar.slider('Minimum diagonal line length', 0, 250, 2, step=1)
+
+min_white_vert_line_len = sidebar.slider('Minimum white vertical line length', 0, 250, 2, step=1)
+
+percentile = sidebar.slider('Precentile', 0, 100, 24, step=1)
+
+sidebar.download_button('Download file', download(),file_name='archive.zip')
+
+# ---------------Plot RPs--------------------
+# runs_ = ['Baseline open eyes', 'Baseline closed eyes', 'Motor execution: left vs right hand', 'Motor imagery: left vs right hand',
+#     'Motor execution: hands vs feet', 'Motor imagery: hands vs feet']
+#
+# options = st.multiselect('Select two runs to compare', runs_, ['Baseline open eyes', 'Baseline closed eyes'])
+
+
+# run_list = []
+#
+# for v in options:
+#     run_list.append(runs_.index(v)+1)
+# if len(run_list) <= 1:
+#     run_list = [1,2]
+st.markdown('Baseline open eyes vs Baseline closed eyes')
+options = ['Baseline open eyes', 'Baseline closed eyes']
+run_list = [1,2]
+
+rp_plot, matrix_open_binary, matrix_close_binary, epochs, stft1, stft2  = run_computation(t_start, t_end, selected_subject, fir_filter, electrode_name, cut_freq, win_len, n_fft, percentile, run_list,options)
+st.write(rp_plot)
+
+# ---------------Plot Spectrum-------------------- 
+st.write(waterfall_spectrum(stft1, stft2, 160, cut_freq, options))
+
+# ---------------Save RPs--------------------
+if st.button('Save RPs as *.npy'):
+    save(matrix_open_binary, matrix_close_binary)
+
+# ---------------Plot Radar--------------------
+rqa_radar = plot_rqa(matrix_open_binary, matrix_close_binary, min_vert_line_len, min_diagonal_line_len, min_white_vert_line_len, options)
+st.write(rqa_radar)
+
+
+
+

complexRadar.py

@@ -0,0 +1,168 @@
+import numpy as np 
+import textwrap
+
+class ComplexRadar():
+    """
+    Create a complex radar chart with different scales for each variable
+    Parameters
+    ----------
+    fig : figure object
+        A matplotlib figure object to add the axes on
+    variables : list
+        A list of variables
+    ranges : list
+        A list of tuples (min, max) for each variable
+    n_ring_levels: int, defaults to 5
+        Number of ordinate or ring levels to draw
+    show_scales: bool, defaults to True
+        Indicates if we the ranges for each variable are plotted
+    format_cfg: dict, defaults to None
+        A dictionary with formatting configurations
+    """
+    def __init__(self, fig, variables, ranges, n_ring_levels=5, show_scales=True, format_cfg=None):
+        
+        # Default formatting
+        self.format_cfg = {
+            # Axes
+            # https://matplotlib.org/stable/api/figure_api.html
+            'axes_args': {},
+            # Tick labels on the scales
+            # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.rgrids.html
+            'rgrid_tick_lbls_args': {'fontsize':8},
+            # Radial (circle) lines
+            # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.grid.html
+            'rad_ln_args': {},
+            # Angle lines
+            # https://matplotlib.org/3.2.2/api/_as_gen/matplotlib.lines.Line2D.html#matplotlib.lines.Line2D
+            'angle_ln_args': {},
+            # Include last value (endpoint) on scale
+            'incl_endpoint':False,
+            # Variable labels (ThetaTickLabel)
+            'theta_tick_lbls':{'va':'top', 'ha':'center'},
+            'theta_tick_lbls_txt_wrap':15,
+            'theta_tick_lbls_brk_lng_wrds':False,
+            'theta_tick_lbls_pad':25,
+            # Outer ring
+            # https://matplotlib.org/stable/api/spines_api.html
+            'outer_ring':{'visible':True, 'color':'#d6d6d6'}
+        }
+        
+        if format_cfg is not None:
+            self.format_cfg = { k:(format_cfg[k]) if k in format_cfg.keys() else (self.format_cfg[k]) 
+                 for k in self.format_cfg.keys()}        
+        
+        
+        # Calculate angles and create for each variable an axes
+        # Consider here the trick with having the first axes element twice (len+1)
+        angles = np.arange(0, 360, 360./len(variables))
+        axes = [fig.add_axes([0.1,0.1,0.9,0.9], 
+                             polar=True,
+                             label = "axes{}".format(i),
+                             **self.format_cfg['axes_args']) for i in range(len(variables)+1)]
+        
+        # Ensure clockwise rotation (first variable at the top N)
+        for ax in axes:
+            ax.set_theta_zero_location('N')
+            ax.set_theta_direction(-1)
+            ax.set_axisbelow(True)
+        
+        # Writing the ranges on each axes
+        for i, ax in enumerate(axes):
+
+            # Here we do the trick by repeating the first iteration
+            j = 0 if (i==0 or i==1) else i-1
+            ax.set_ylim(*ranges[j])
+            # Set endpoint to True if you like to have values right before the last circle
+            grid = np.linspace(*ranges[j], num=n_ring_levels, 
+                               endpoint=self.format_cfg['incl_endpoint'])
+            gridlabel = ["{}".format(round(x,2)) for x in grid]
+            gridlabel[0] = "" # remove values from the center
+            lines, labels = ax.set_rgrids(grid, 
+                                          labels=gridlabel, 
+                                          angle=angles[j],
+                                          **self.format_cfg['rgrid_tick_lbls_args']
+                                         )
+            
+            ax.set_ylim(*ranges[j])
+            ax.spines["polar"].set_visible(False)
+            ax.grid(visible=False)
+            
+            if show_scales == False:
+                ax.set_yticklabels([])
+
+        # Set all axes except the first one unvisible
+        for ax in axes[1:]:
+            ax.patch.set_visible(False)
+            ax.xaxis.set_visible(False)
+            
+        # Setting the attributes
+        self.angle = np.deg2rad(np.r_[angles, angles[0]])
+        self.ranges = ranges
+        self.ax = axes[0]
+        self.ax1 = axes[1]
+        self.plot_counter = 0
+        
+        
+        # Draw (inner) circles and lines
+        self.ax.yaxis.grid(**self.format_cfg['rad_ln_args'])
+        # Draw outer circle
+        self.ax.spines['polar'].set(**self.format_cfg['outer_ring'])
+        # Draw angle lines
+        self.ax.xaxis.grid(**self.format_cfg['angle_ln_args'])
+
+        # ax1 is the duplicate of axes[0] (self.ax)
+        # Remove everything from ax1 except the plot itself
+        self.ax1.axis('off')
+        self.ax1.set_zorder(9)
+        
+        # Create the outer labels for each variable
+        l, text = self.ax.set_thetagrids(angles, labels=variables)
+        
+        # Beautify them
+        labels = [t.get_text() for t in self.ax.get_xticklabels()]
+        labels = ['\n'.join(textwrap.wrap(l, self.format_cfg['theta_tick_lbls_txt_wrap'], 
+                                          break_long_words=self.format_cfg['theta_tick_lbls_brk_lng_wrds'])) for l in labels]
+        self.ax.set_xticklabels(labels, **self.format_cfg['theta_tick_lbls'])
+        
+        for t,a in zip(self.ax.get_xticklabels(),angles):
+            if a == 0:
+                t.set_ha('center')
+            elif a > 0 and a < 180:
+                t.set_ha('left')
+            elif a == 180:
+                t.set_ha('center')
+            else:
+                t.set_ha('right')
+
+        self.ax.tick_params(axis='both', pad=self.format_cfg['theta_tick_lbls_pad'])
+
+
+    def _scale_data(self, data, ranges):
+        """Scales data[1:] to ranges[0]"""
+        for d, (y1, y2) in zip(data[1:], ranges[1:]):
+            assert (y1 <= d <= y2) or (y2 <= d <= y1)
+        x1, x2 = ranges[0]
+        d = data[0]
+        sdata = [d]
+        for d, (y1, y2) in zip(data[1:], ranges[1:]):
+            sdata.append((d-y1) / (y2-y1) * (x2 - x1) + x1)
+        return sdata
+        
+    def plot(self, data, *args, **kwargs):
+        """Plots a line"""
+        sdata = self._scale_data(data, self.ranges)
+        self.ax1.plot(self.angle, np.r_[sdata, sdata[0]], *args, **kwargs)
+        self.plot_counter = self.plot_counter+1
+    
+    def fill(self, data, *args, **kwargs):
+        """Plots an area"""
+        sdata = self._scale_data(data, self.ranges)
+        self.ax1.fill(self.angle, np.r_[sdata, sdata[0]], *args, **kwargs)
+        
+    def use_legend(self, *args, **kwargs):
+        """Shows a legend"""
+        self.ax1.legend(*args, **kwargs)
+    
+    def set_title(self, title, pad=25, **kwargs):
+        """Set a title"""
+        self.ax.set_title(title,pad=pad, **kwargs)

papaerspace_config.txt

@@ -0,0 +1,5 @@
+image: datalab108/cookiewriter64:test
+port: 5005
+resources:
+  replicas: 1
+  instanceType: P4000

requirements.txt

@@ -0,0 +1,10 @@
+torch
+torchvision
+torchaudio
+umap-learn
+umap-learn[plot]
+PyRQA
+pandas
+mne
+plotly
+streamlit

test.py

@@ -0,0 +1,5 @@
+import time
+
+while True:
+    print("Hello World")
+    time.sleep(10)