Splitting Simulator.py into a simulator and feature extraction script

Simulator.py

+# Function to run simulations with the computational model, given a parameter set.
+
 from brian2 import *
 from numpy import *
 from scipy import signal
 import pickle
 import torch
 import sys
-from scipy.signal import find_peaks
-from scipy.fft import fft, fftfreq
-from scipy.stats import norm
-from itertools import combinations
 import numpy as np
 
-
 def MEAnetsimulate(parameter_set, simtime=165*second, transient=5*second):
     # Network parameters
     Nl = 10
@@ -176,186 +173,4 @@ def MEAnetsimulate(parameter_set, simtime=165*second, transient=5*second):
 
     APs = APs.reshape(len(APs) // 2, 2)
 
-    return APs, simtime, transient, fs
-
-
-def compute_features(APs, simtime, transient, fs):
-    numelectrodes = 12                              # number of electrodes
-    timbin = 25 * ms                                # timebin to compute network firing rate
-
-    # Initialize
-    APs_inbin = zeros((numelectrodes, int(floor((simtime - transient) / timbin))))
-    spikerate = zeros((numelectrodes, int(floor((simtime - transient) / timbin))))
-
-    # Delete transient
-    APs_wot = APs[APs[:, 1] > transient * fs / second, :]
-    APs_wot[:, 1] = APs_wot[:, 1] - transient * fs / second
-
-    # Calculate the network firing rate
-    for l in range(numelectrodes):
-        APt = APs_wot[APs_wot[:, 0] == l, :]         # go through one electrode first
-        APt = APt[:, 1]
-        binsize = timbin * fs / second
-        for k in range(int(floor((simtime - transient) / timbin))):
-            APs_inbin[l, k] = sum((APt > (k * binsize)) & (APt < ((k + 1) * binsize)))
-            spikerate[l, k] = APs_inbin[l, k] * second / timbin
-
-    APs_inbintot = sum(APs_inbin, axis=0)
-    spikeratetot = sum(spikerate, axis=0)
-
-    # Smoothen the spikerate by convolution with gaussian kernel 
-    width = 11                      
-    sigma = 3.0
-    x = np.arange(0, width, 1, float)
-    x = x - width // 2
-    kernel = norm.pdf(x, scale=sigma)
-    kernel /= np.sum(kernel)
-    spikeratesmooth = np.convolve(spikeratetot, kernel, mode='same')
-
-    # Detect fragmented bursts on smoothed spikerate
-    MB_th = (1/16) * max(spikeratetot)
-    peaks, ph = find_peaks(spikeratesmooth, height=MB_th, prominence=(1 / 20) * max(spikeratetot))
-
-    # Set parameters for burst detection
-    act_elec = sum(mean(spikerate, axis=1) > 0.02)      # calculate the number of active electrodes
-
-    start_th = 0.25 * max(spikeratetot)                 # spikerate threshold to start a burst
-    t_th = int((50 * ms) / timbin)                      # how long it has to surpass threshold for to start burst
-    e_th = 0.5 * act_elec                               # how many electrodes need to be active in the burst
-    stop_th = (1 / 50) * max(spikeratetot)              # threshold to end a burst
-
-    # Initialize burst detection
-    i = 0
-    NB_count = 0
-    max_NBs = 1000                                      # maximum amount of to be detected bursts
-    NBs = zeros((max_NBs, 4))
-
-    # Detect NBs
-    while (i + t_th) < len(spikeratetot):
-        if (all(spikeratetot[i:i + t_th] > start_th)) \
-                & (sum(sum(APs_inbin[:, i:i + t_th], axis=1) > t_th) > e_th):
-            NBs[NB_count, 2] = NBs[NB_count, 2] + sum(APs_inbintot[i:i + t_th])
-            NBs[NB_count, 0] = i
-            i = i + t_th
-            while any(spikeratetot[i:i + 2 * t_th] > stop_th):
-                NBs[NB_count, 2] = NBs[NB_count, 2] + APs_inbintot[i]
-                i = i + 1
-            NBs[NB_count, 3] = sum((peaks > NBs[NB_count, 0]) & (peaks < i))
-            NBs[NB_count, 1] = i
-            NB_count = NB_count + 1
-        else:
-            i = i + 1
-
-    NBs = NBs[0:NB_count, :]
-
-    MNBR = NB_count * 60 * second / (simtime - transient)
-    NBdurations = (array(NBs[:, 1]) - array(NBs[:, 0])) * timbin / second
-    MNBD = mean(NBdurations)
-    PSIB = sum(NBs[:, 2] / len(APs_wot)) * 100
-    MFR = len(APs_wot) / ((simtime - transient) / second) / numelectrodes
-    IBI = (array(NBs[1:, 0]) - array(NBs[0:-1, 1])) * timbin / second
-    CVIBI = np.std(IBI) / np.mean(IBI)
-    if NB_count == 0:
-        MNBD = 0.0
-        MNMBs = 0.0
-        NFBs = 0
-    else:
-        NFBs = sum(NBs[:, 3]) / NB_count
-
-    if NB_count < 2:
-        CVIBI = 0.0
-
-    # Calculate MAC metric as defined by Maheswaranathan
-    yf = fft(spikeratetot)
-    xf = fftfreq(len(spikeratetot), timbin / second)[:len(spikeratetot) // 2]
-    MAC = max(np.abs(yf[1:len(spikeratetot)])) / np.abs(yf[0])
-
-    # Calculate cross-correlation between binarized spike trains
-    # Binarize the spike trains 
-    all_combinations = list(combinations(list(arange(numelectrodes)), 2))
-    trans_timebin = 0.2 * second                         # timebin to transform spiketrains to binary
-    bin_timeseries = list(range(0, int((simtime - transient) / ms), int(trans_timebin / ms)))
-    binary_signal = zeros((numelectrodes, len(bin_timeseries)))
-    for i in range(numelectrodes):
-        signal = APs_inbin[i, :]
-        grosssignal = [sum(signal[x:x + int((trans_timebin / timbin))]) for x in
-                       range(0, len(signal), int(trans_timebin / timbin))]
-        binary_signal[i, :] = [1 if x > 0 else 0 for x in grosssignal]
-
-    # Calculate coefficients between every pair of electrodes
-    coefficients = []
-    N = len(binary_signal[0, :])
-    for i, j in all_combinations:
-            signal1 = binary_signal[i, :]
-            signal2 = binary_signal[j, :]
-            if (i != j) & (not list(signal1) == list(signal2)):
-                coefficients.append((N * sum(signal1 * signal2) - sum(signal1) * (sum(signal2)))
-                                    * ((N * sum(signal1 ** 2) - sum(signal1) ** 2) ** (-0.5))
-                                    * ((N * sum(signal2 ** 2) - sum(signal2) ** 2) ** (-0.5)))
-
-    mean_corr = mean(coefficients)
-    sd_corr = std(coefficients)
-
-    if not coefficients:
-        mean_corr = 1
-        sd_corr = 0
-        
-    # Compute continuous ISI arrays
-    time_vector = np.arange(0, (simtime - transient) / second, 1/fs)
-    isi_arrays = np.zeros((numelectrodes, len(time_vector)))
-
-    for electrode in range(numelectrodes):
-        # Extract spike times for the current electrode
-        electrode_spike_times = APs_wot[APs_wot[:, 0] == electrode, 1]
-
-        for i in range(len(electrode_spike_times) - 1):
-            spike1 = electrode_spike_times[i]
-            spike2 = electrode_spike_times[i + 1]
-            tisi = (spike2 - spike1) / fs
-
-            # Fill ISI values in the appropriate range
-            if i == 0:
-                isi_arrays[electrode, 0:spike1] = NaN
-            isi_arrays[electrode, spike1:spike2] = tisi
-            if (i + 1) == (len(electrode_spike_times) - 1):
-                isi_arrays[electrode, spike2:] = NaN
-
-    # Compute ISI measures
-    meanisi_array = np.nanmean(isi_arrays, axis=0)
-    mean_ISI = np.nanmean(meanisi_array)
-    sdmean_ISI = np.nanstd(meanisi_array)
-    sdtime_ISI = np.nanmean(np.nanstd(isi_arrays, axis=0))
-
-    # Calculate the ISI-distance and ISI correlations
-    ISI_distances = np.zeros(len(all_combinations))
-    isicoefficients = np.zeros(len(all_combinations))
-    N = len(isi_arrays[0, :])
-    j = 0
-    
-    # Iterate through the electrode combinations
-    for electrode1_key, electrode2_key in all_combinations:
-        # Get the ISI arrays for the selected electrodes
-        isit1_wn = isi_arrays[electrode1_key, :]
-        isit2_wn = isi_arrays[electrode2_key, :]
-        isit1 = isit1_wn[~isnan(isit1_wn)]
-        isit2_wn = isit2_wn[~isnan(isit1_wn)]
-        isit2 = isit2_wn[~isnan(isit2_wn)]
-        isit1 = isit1[~isnan(isit2_wn)]
-
-        isi_diff = isit1 / isit2
-        ISI_distances[j] = np.mean(np.where(isi_diff <= 1, abs(isi_diff - 1), -1 * (1 / isi_diff - 1)))
-
-        if (i != j) & (not list(isit1) == list(isit2)):
-            isicoefficients[j] = ((N * sum(isit1 * isit2) - sum(isit1) * (sum(isit2)))
-                                  * ((N * sum(isit1 ** 2) - sum(isit1) ** 2) ** (-0.5))
-                                  * ((N * sum(isit2 ** 2) - sum(isit2) ** 2) ** (-0.5)))
-
-        j += 1
-
-    mean_ISIcorr = mean(isicoefficients)
-    sd_ISIcorr = std(isicoefficients)
-    ISI_distance = np.mean(ISI_distances)
-
-    return [MFR, MNBR, MNBD, PSIB, NFBs, CVIBI, mean_corr, sd_corr, mean_ISIcorr, sd_ISIcorr, ISI_distance, mean_ISI, sdmean_ISI, sdtime_ISI, MAC]
-
-
+    return APs, simtime, transient, fs
\ No newline at end of file