binned_spiking_activity_human_epilepsy/python_code/analysis_PlosCB/mr_epilepsy_timeresolved.py

import numpy as np
import pandas as pd
import mrestimator as mre
import os
import time
import mr_utilities
import math


def analyse_timeresolved_mr(data_path, result_path, recordings,
                            binnings='all', ksteps=(1,400),
                            windowsize=20000, windowstep=500,
                            fit_method="offset"):

    for recording in recordings:
        outputfile = '{}{}/pkl/activity_SU_4ms.pkl'.format(data_path, recording)

        if os.path.isfile(outputfile):
            data = pd.read_pickle(outputfile)
        else:
            print('Error: no datafile for recording {}'.format(recording))
            continue

        existing_binnings = data['binning'].unique().tolist()
        patient = int(data['patient'].unique()[0])

        if binnings == 'all':
            binnings = existing_binnings

        for binning in binnings:

            dt = data['dt'].unique()[0]
            activity_list = data[data['binning'] == binning]['activity'].tolist()

            if len(activity_list) == 0:
                continue
            if len(activity_list) > 1:
                print('Error: Multiple activity entries for binning.')
                data = data[~data['activity'].apply(tuple).duplicated()]
                activity_list = data[data['binning'] == binning][
                    'activity'].tolist()
                activity = activity_list[0]
            else:
                activity = activity_list[0]

            t_onset = 600 # in seconds, given 10min recordings
            # the recording is significantly less than 10min long,
            # it is unclear whether seizure onset is indeed at the end of
            # the recording, therefore better skip
            if len(activity) < 0.8 * t_onset * 10**3 / dt:
                continue

            save_path_recording = '{}{}/'.format(result_path, recording)

            timeresolved_mr(activity, windowsize, windowstep,
                            ksteps[0], ksteps[1], fit_method,
                            save_path_recording, patient, recording,
                            binning, dt)
            print('break to cool cpu...')
            time.sleep(30)


def timeresolved_mr(activity, windowsize, windowstep, kmin, kmax,
                    fit_method, save_path_recording, patient, recording,
                    binning, dt):

    """ performs MR estimation in sliding windows of provided activity and
    calls function to store results in pkl file
    """

    for t in range(0, len(activity) - windowsize, windowstep):
        windowed_data = activity[t:t + windowsize]
        prepared = mre.input_handler(windowed_data)
        rk = mre.coefficients(prepared, dt=dt, dtunit='ms', steps=(kmin, kmax))
        try:
            fit_res = mre.fit(rk, fitfunc=mr_utilities.string_to_func(
                fit_method))
            t0 = t * dt

            mt_savefile(save_path_recording, patient, recording, binning,
                        dt, kmin, kmax, fit_method, windowsize, windowstep,
                        t0=t0, activity=windowed_data, rk=rk, fitres=fit_res)
        except (RuntimeError, ValueError):
            continue


def mt_savefile(result_path, patient, recording, binning, dt, kmin, kmax,
                fit_method, windowsize, windowstep, t0, activity, rk,
                fitres):

    """ save results of m(t) estimation performed by function
    timeresolved_mr to pkl file.
    """

    datafile = '{}mt_results_{}_{}_kmin{}_kmax{}_winsize{}_winstep{}_dt{}' \
               '.pkl'.format(result_path, binning, fit_method, kmin, kmax,
                             windowsize, windowstep, dt)

    if math.isnan(fitres.tau):
        inconsistencies = ['not_converged']
    else:
        inconsistencies = mr_utilities.consistency_checks(activity, rk, fitres)
    try:
        fano = np.var(activity)/np.mean(activity)
        variance = np.var(activity)
    except ZeroDivisionError:
        fano = np.nan
        variance = np.nan

    save_res = dict()
    save_res['patient'] = patient
    save_res['recording'] = recording
    save_res['binning'] = binning
    save_res['kmin'] = kmin
    save_res['kmax'] = kmax
    save_res['windowsize'] = windowsize
    save_res['windowstep'] = windowstep
    save_res['t0'] = t0
    save_res['rk'] = [list(rk.coefficients)]
    save_res['rk_steps'] = [list(rk.steps)]
    save_res['dt'] = dt
    save_res['fit_method'] = fit_method
    save_res['fit_params'] = [list(fitres.popt)]
    save_res['m'] = fitres.mre
    save_res['tau'] = fitres.tau
    save_res['inconsistencies'] = inconsistencies
    save_res['fano'] = fano
    save_res['variance'] = variance
    save_res['offset'] = fitres.popt[2]

    if os.path.isfile(datafile):
        mr_results_patient = pd.read_pickle(datafile)
        mr_results_patient = mr_results_patient.append(save_res,
                                                       ignore_index=True)
    else:
        mr_results_patient = pd.DataFrame.from_records([save_res])
    mr_results_patient.to_pickle(datafile)


def mt_x_parts(x, outputfilename, data_path, result_path, recordings,
               binnings='all', ksteps=(1, 400), fit_method ="offset"):

    mx_allresults = None
    for recording in recordings:
        outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)

        if os.path.isfile(outputfile):
            data = pd.read_pickle(outputfile)
        else:
            print(' File not found {}'.format(outputfile))
            continue

        existing_binnings = data['binning'].unique().tolist()
        patient = int(data['patient'].unique()[0])

        kmin = ksteps[0]
        kmax = ksteps[1]

        if binnings == 'all':
            binnings = existing_binnings

        for binning in binnings:
            dt = data['dt'].unique()[0]
            activity_list = data[data['binning'] == binning]['activity'].tolist()
            activity = activity_list[0]
            partlength = int(np.floor(len(activity)/x))
            tau = []
            m = []
            binning_focus = mr_utilities.binning_focus(binning, patient)

            for i in range(x):
                activity_part = activity[i*partlength:(i+1)*partlength]
                nspikes = np.sum(np.array(activity_part))
                n_nonzero = mr_utilities.count_nonzero(activity_part)
                prepared = mre.input_handler(activity_part)
                rk = mre.coefficients(prepared, dt=dt, dtunit='ms',
                                      steps=(kmin, kmax))
                try:
                    fit_res = mre.fit(rk,
                                      fitfunc=mr_utilities.string_to_func(
                                          fit_method))

                    if not math.isnan(fit_res.tau):
                        tau.append(fit_res.tau)
                        m.append(fit_res.mre)

                        inconsistencies = mr_utilities.consistency_checks(
                            activity_part, rk, fit_res)
                        fano = np.var(activity_part)/np.mean(activity_part)
                        variance = np.var(activity_part)
                        mx_result = {'patient': patient,
                                      'recording': recording,
                                      'binning': binning,
                                      'binning_focus': binning_focus,
                                      'partno': i+1,
                                      'm': fit_res.mre,
                                      'tau': fit_res.tau,
                                      'fano': fano,
                                      'variance': variance,
                                      'offset': fit_res.popt[2],
                                      'rk': [list(rk.coefficients)],
                                      'fit_params': [list(fit_res.popt)],
                                      'nspikes': nspikes,
                                      'n_nonzero': n_nonzero,
                                      'inconsistencies': inconsistencies,
                                      'kmin': kmin,
                                      'kmax': kmax,
                                      'fit_method': fit_method}

                        if mx_allresults is None:
                            mx_allresults = pd.DataFrame.from_records(
                                [mx_result])

                        else:
                            mx_allresults = mx_allresults.append(
                                mx_result, ignore_index=True)

                except (RuntimeError, ValueError):
                    continue

    datafile_mx = '{}mt_{}_results_allrecords_kmin{}_kmax{}_{}.pkl'.format(
        result_path, x, ksteps[0], ksteps[1], fit_method)

    if mx_allresults is not None:
        mx_allresults.to_pickle(datafile_mx)