binned_spiking_activity_human_epilepsy/python_code/analysis_PlosCB/mr_epilepsy_fullrec.py

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


def mr_analysis_confint(data_path, result_path, outputfilename,
                        recordings='all', binnings=[], ksteps=(1, 400),
                        fit_method='offset'):

    if recordings == 'all':
        recordinglist = [dI for dI in os.listdir(data_path) if
                         os.path.isdir(os.path.join(data_path, dI))]
    else:
        recordinglist = recordings

    focus_list = '../Data/patients.txt'
    f = open(focus_list, 'r')
    foci = [line.rstrip('\n').split('\t') for line in f.readlines()]

    # -----------------------------------------------------------------------#
    # Run through all recordings and binnings and perform MR estimation
    # -----------------------------------------------------------------------#
    for recording in recordinglist:
        outputfile = '{}{}/{}'.format(data_path, recording, outputfilename)

        if not os.path.isfile(outputfile):
            print("No binning data for recording ", recording)
            print(outputfile)
            continue

        data = pd.read_pickle(outputfile)
        patient = int(data['patient'].unique()[0])

        # read the location of the focus
        focus = 'NA'
        for i, idf in enumerate(foci[1:]):
            if int(idf[0]) == int(patient):
                focus = idf[1]
        if focus == 'NA':
            print('Error: patient {} not in foci list.'.format(patient))
            continue

        existing_binnings = data['binning'].tolist()
        if not binnings:
            binnings = existing_binnings

        for binning in binnings:

            binning_focus = mr_utilities.binning_focus(binning, patient)
            dt = data['dt'].unique()[0]
            activity = data[data['binning'] == binning]['activity'].tolist()[0]

            if mr_utilities.count_nonzero(activity) < 1000:
                continue

            # split into trials for bootstrap confidence intervals
            ntrials = 30
            len_trial = int(np.floor(len(activity) / ntrials))
            unusable_activity = len(activity) - ntrials * len_trial
            activity = activity[unusable_activity:]
            trials = np.reshape(activity, (ntrials, len_trial))

            # apply MR estimator
            prepared = mre.input_handler(trials)
            rk = mre.coefficients(prepared, dt=dt, dtunit='ms',
                                  steps=ksteps, method='sm', knownmean=None)
            try:
                maxfev = 200 * (len(rk.coefficients)+1)
                fit_res = mre.fit(rk, fitfunc=mr_utilities.string_to_func(
                                    fit_method), maxfev=maxfev)

                if not math.isnan(fit_res.tau):
                    inconsistencies = mr_utilities.consistency_checks(
                                            activity, rk, fit_res)
                    mr_savefile(result_path, patient, recording,
                                focus, binning, binning_focus, ksteps,
                                fit_method, rk, fit_res, inconsistencies)

            except RuntimeError:
                print('RUNTIME-ERROR', recording, binning)
                pass


def mr_savefile(result_path, patient, recording, focus, binning,
                binning_focus, ksteps, fit_method, rk, fitres,
                inconsistencies):

    datafile_patient = result_path + 'mr_results.pkl'
    m_quantiles = fitres.mrequantiles
    tau_quantiles = fitres.tauquantiles
    rk_stderr = rk.stderrs

    save_res = dict()
    save_res['patient'] = patient
    save_res['recording'] = recording
    save_res['focus'] = focus
    save_res['binning_focus'] = binning_focus
    save_res['binning'] = binning
    save_res['kmin'] = ksteps[0]
    save_res['kmax'] = ksteps[1]
    save_res['rk'] = [list(rk.coefficients)]
    save_res['rk_steps'] = [list(rk.steps)]
    save_res['dt'] = rk.dt
    save_res['fit_method'] = fit_method
    save_res['fit_params'] = [list(fitres.popt)]
    save_res['m'] = fitres.mre
    save_res['ssres'] = fitres.ssres
    save_res['ssres_weighted'] = fitres.ssres*len(fitres.popt)
    save_res['tau'] = fitres.tau
    save_res['inconsistencies'] = inconsistencies
    save_res['rk_stderr'] = rk_stderr
    save_res['m_quantiles'] = m_quantiles
    save_res['tau_quantiles'] = tau_quantiles

    if os.path.isfile(datafile_patient):
        mr_results_patient = pd.read_pickle(datafile_patient)
        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_patient)