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)