##from update_project import **#model_and_data import os import numpy as np import pandas as pd from IPython import embed from matplotlib import gridspec, pyplot as plt try: from plotstyle import plot_style, spines_params except: print('plotstyle not installed') from threefish.core_plot_subplots import plot_lowpass2, plt_single_square_modl, plt_time_arrays from threefish.core_plot_suscept import perc_model_full, plt_data_susept from threefish.core_filenames import overlap_cells from threefish.core_load import resave_small_files, save_visualization from threefish.core_reformat_RAM import get_flowchart_params, load_model_susept from threefish.core import find_folder_name from threefish.core_plot_labels import label_noise_name, nonlin_title, remove_xticks, remove_yticks, set_xlabel_arrow, \ set_ylabel_arrow, title_find_cell_add, xlabel_xpos_y_modelanddata import itertools as it from threefish.defaults import default_figsize, default_settings from threefish.core_plot_lims import join_x, join_y, set_clim_same, set_same_ylim #from plt_RAM import model_and_data, model_and_data_sheme, model_and_data_vertical2 def table_printen(table): print(table.keys()) for l in range(len(table)): list_here = np.array(table.iloc[l]) l1 = "& ".join(list_here) print(l1) def model_and_data2(eod_metrice=False, width=0.005, nffts=['whole'], powers=[1], cells=["2013-01-08-aa-invivo-1"], contrasts=[0], noises_added=[''], D_extraction_method=['additiv_cv_adapt_factor_scaled'], internal_noise=['RAM'], external_noise=['RAM'], level_extraction=[''], receiver_contrast=[1], dendrids=[''], ref_types=[''], adapt_types=[''], c_noises=[0.1], c_signal=[0.9], cut_offs1=[300]): # ['eRAM'] stimulus_length = 1 # 20#550 # 30 # 15#45#0.5#1.5 15 45 100 trials_nrs = [1] # [100, 500, 1000, 3000, 10000, 100000, 1000000] # 500 rep = 1000000 # 500000#0 good_data, remaining = overlap_cells() cells_all = [good_data[0]] plot_style() default_figsize(column=2, length=3.1) #2.9.254.75 0.75# bottom=0.07, top=0.94, grid = gridspec.GridSpec(2, 5, wspace=0.95, bottom=0.13, hspace=0.60, top=0.88, width_ratios=[2, 0, 2, 2, 2], left=0.09, right=0.93, ) #bottom=0.09, hspace=0.25, top=0.9, a = 0 maxs = [] mins = [] mats = [] ims = [] iternames = [D_extraction_method, external_noise, internal_noise, powers, nffts, dendrids, cut_offs1, trials_nrs, c_signal, c_noises, ref_types, adapt_types, noises_added, level_extraction, receiver_contrast, contrasts, ] nr = '2' # cell_contrasts = ["2013-01-08-aa-invivo-1"] # cells_triangl_contrast = np.concatenate([cells_all,cell_contrasts]) # cells_triangl_contrast = 1 # cell_contrasts = 1 rows = len(cells_all) # len(good_data)+len(cell_contrasts) perc = 'perc' lp = 2 label_model = r'Nonlinearity $\frac{1}{S}$' for all in it.product(*iternames): var_type, stim_type_afe, stim_type_noise, power, nfft, dendrid, cut_off1, trial_nrs, c_sig, c_noise, ref_type, adapt_type, noise_added, extract, a_fr, a_fe = all # print(trials_stim,stim_type_noise, power, nfft, a_fe,a_fr, dendrid, var_type, cut_off1,trial_nrs) fig = plt.figure() hs = 0.45 ################################# # data cells grid_data = gridspec.GridSpecFromSubplotSpec(1, 1, grid[0, 0], hspace=hs) fr_print = False nr = 1 ax_data, stack_spikes_all, eod_frs = plt_data_susept(fig, grid_data, cells_all, cell_type='p-unit', width=width, cbar_label=True, fr_print=fr_print, eod_metrice=eod_metrice, nr=nr, amp_given=1, xlabel=False, lp=lp, title=True) for ax_external in ax_data: ax_external.set_xticks_delta(100) set_ylabel_arrow(ax_external, xpos=xlabel_xpos_y_modelanddata(), ypos=0.87) set_xlabel_arrow(ax_external, ypos=ypos_x_modelanddata()) #embed() #plt.show() ################################## # model part trial_nr = 500000 cell = '2013-01-08-aa-invivo-1' cell = '2012-07-03-ak-invivo-1' print('cell' + str(cell)) cells_given = [cell] save_name_rev = find_folder_name( 'calc_model') + '/' + 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV_revQuadrant_' # for trial in trials:#.009 trial_nr = 1000000 #1000000 save_names = [ 'calc_RAM_model-2__nfft_whole_power_1_afe_0.025_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_afe_0.025_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', ] save_names = [ 'calc_RAM_model-2__nfft_whole_power_1_afe_0.023_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_afe_0.023_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', ] #calc_RAM_model-2__nfft_whole_power_1_afe_2.6_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV ########## # Erklärung # Ich habe hier 0.009 und nicht 0.25 weil das Modell einen Fehler hat # den Stimulus in den Daten habe ich überprüft der tatsächliche stimulus ist 2.3 Prozent # sollte aber 2.5 Prozent sein # Im Fall von 2.5 Prozent wäre das ein Fehler von 0.36 sonst von 0.39 # Hier werde ich nun mit dem Fehler von 0.36 verfahren # das bedeutet aber das sich den Stimulus zwar mit 0.009 ins Modell reintue später für die # Susceptiblitätsberechnung sollte ich ihn aber um den Faktor 0.36 teilen # oben habe ich einen bias factor weil die Zelle zu sensitiv gefittet ist, also passe ich das an dass die den # gleichen CV und feurrate hat, wie die Zelle in der Stimulation, deswegen ist dieser Bias faktor nur oben! bias_factors = [0.36, 0.36, 1, 1]#0.36 save_names = [ 'calc_RAM_model-2__nfft_whole_power_1_afe_0.009_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_afe_0.009_RAM_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', 'calc_RAM_model-2__nfft_whole_power_1_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_' + str( trial_nr) + '_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', ] #'calc_RAM_model-2__nfft_whole_power_1_afe_0.009_RAM_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_11_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', #'calc_RAM_model-2__nfft_whole_power_1_afe_0.009_RAM_RAM_additiv_cv_adapt_factor_scaled_cNoise_0.1_cSig_0.9_cutoff1_300_cutoff2_300no_sinz_length1_TrialsStim_500000_a_fr_1__trans1s__TrialsNr_1_fft_o_forward_fft_i_forward_Hz_mV', nrs_s = [3, 4, 8, 9] #, 10, 11 #embed() tr_name = trial_nr / 1000000 if tr_name == 1: tr_name = 1 #'$c=1\,\%$','$c=0\,\%$' c = 2.5 cs = ['$c=%.1f$' % c + '$\,\%$', '$c=0\,\%$'] titles = ['Model\n$N=11$', 'Model\n' + '$N=10^6$', 'Model\,(' + label_noise_name().lower() + ')' + '\n' + '$N=11$', 'Model\,(' + label_noise_name().lower() + ')' + '\n' + '$N=10^6$' ] #%#%s$' % (tr_name) + '\,million' #'Model\,('+noise_name().lower()+')' + '\n' + '$N=11$\n $c=1\,\%$',$N=%s $' % (tr_name) +'\,million' # 'Model\,('+noise_name().lower()+')' + '\n' + '$N=%s$' % (tr_name) + '\,million\n $c=1\,\%$ ' ax_model = [] for s, sav_name in enumerate(save_names): try: ax_external = plt.subplot(grid[nrs_s[s]]) except: print('vers something') embed() ax_model.append(ax_external) save_name = find_folder_name('calc_model') + '/' + sav_name cell_add, cells_save = title_find_cell_add(cells_given) perc = 'perc' path = save_name + '.pkl' # '../'+ # model = get_stack_one_quadrant(cell, cell_add, cells_save, path, save_name) # full_matrix = create_full_matrix2(np.array(model), np.array(stack_rev)) # stack_final = get_axis_on_full_matrix(full_matrix, model) # im = plt_RAM_perc(ax, perc, np.abs(model)) model = load_model_susept(path, cells_save, save_name.split(r'/')[-1] + cell_add) #embed() if len(model) > 0: add_nonlin_title, cbar, fig, stack_plot, im = plt_single_square_modl(ax_external, cell, model, perc, titles[s], width, eod_metrice=eod_metrice, titles_plot=True, resize=True, bias_factor=bias_factors[s], fr_print=fr_print, nr=nr) # if s in [1,3,5]: ims.append(im) mats.append(stack_plot) maxs.append(np.max(np.array(stack_plot))) mins.append(np.min(np.array(stack_plot))) col = 2 row = 2 ax_external.set_xticks_delta(100) ax_external.set_yticks_delta(100) # cbar[0].set_label(nonlin_title(add_nonlin_title)) # , labelpad=100 cbar.set_label(nonlin_title(' [' + add_nonlin_title), labelpad=lp) # rotation=270, if s in np.arange(col - 1, 100, col): # | (s == 0) remove_yticks(ax_external) else: set_ylabel_arrow(ax_external, xpos=xlabel_xpos_y_modelanddata(), ypos=0.87) if s >= row * col - col: set_xlabel_arrow(ax_external, ypos=ypos_x_modelanddata()) # ax.set_xlabel(F1_xlabel(), labelpad=20) else: remove_xticks(ax_external) if len(cells) > 1: a += 1 set_clim_same(ims, mats=mats, lim_type='up', nr_clim='perc', clims='', percnr=perc_model_full()) ################################################# # Flowcharts ax_ams, ax_external = plt_model_flowcharts(a_fr, ax_external, c, cs, grid, model, stimulus_length) set_same_ylim(ax_ams, up='up') axes = np.concatenate([ax_data, ax_model]) axes = [ax_ams[0], axes[1], axes[2], ax_ams[1], axes[3], axes[4], ] #ax_ams[2], axes[5], axes[6], #axd1 = plt.subplot(grid[1, 1]) #axd2 = plt.subplot(grid[2, 1]) #ax_data.extend([,]) #axd1.show_spines('') #axd2.show_spines('') #embed() #axes = [[ax_ams[0],ax_data[0],axes[2], axes[3]],[ax_ams[1],axd1,axes[4], axes[5]],[axd2,axd2, axes[6], axes[7]]] fig.tag([ax_data], xoffs=-3, yoffs=1.6) # ax_ams[3], fig.tag([axes[0:3]], xoffs=-3, yoffs=1.6) # ax_ams[3], #fig.tag([[axes[4]]], xoffs=-3, yoffs=1.6, minor_index=0) # ax_ams[3], fig.tag([axes[3:6]], xoffs=-3, yoffs=1.6) #, major_index = 1, minor_index = 2 ax_ams[3], #fig.tag([[axes[7]]], xoffs=-3, yoffs=1.6, major_index=2,minor_index=0) # ax_ams[3], #fig.tag([axes[8::]], xoffs=-3, yoffs=1.6, major_index=2, minor_index=2) # ax_ams[3], #fig.tag([axes[7::]], xoffs=-3, yoffs=1.6) # ax_ams[3], #fig.tag([ax_ams[0],ax_data[0],axes[2], axes[3]], xoffs=-3, yoffs=1.6)#ax_ams[3], save_visualization(pdf=True) def plt_model_flowcharts(a_fr, ax_external, c, cs, grid, stack, stimulus_length): var_types = ['', 'additiv_cv_adapt_factor_scaled'] # 'additiv_cv_adapt_factor_scaled', ##additiv_cv_adapt_factor_scaled a_fes = [c / 100, 0] # , 0.009 eod_fe = [750, 750] # , 750 ylim = [-0.5, 0.5] c_sigs = [0, 0.9] # , 0.9 grid_left = [[], grid[1, 0]] # , grid[2, 0] ax_ams = [] for g, grid_here in enumerate([grid[0, 2], grid[1, 2]]): # , grid[2, 0] grid_lowpass = gridspec.GridSpecFromSubplotSpec(4, 1, subplot_spec=grid_here, hspace=0.2, height_ratios=[1, 1, 1, 0.1]) models = resave_small_files("models_big_fit_d_right.csv", load_folder='calc_model_core') model_params = models[models['cell'] == '2012-07-03-ak-invivo-1'].iloc[0] cell = model_params.pop('cell') # .iloc[0]# Werte für das Paper nachschauen eod_fr = model_params['EODf'] # .iloc[0] deltat = model_params.pop("deltat") # .iloc[0] v_offset = model_params.pop("v_offset") # .iloc[0] eod_fr = stack.eod_fr.iloc[0] print(var_types[g] + ' a_fe ' + str(a_fes[g])) noise_final_c, spike_times, stimulus, stimulus_here, time, v_dent_output, v_mem_output, frame = get_flowchart_params( a_fes, a_fr, g, c_sigs[g], cell, deltat, eod_fr, model_params, stimulus_length, v_offset, var_types, eod_fe=eod_fe) if (len(np.unique(frame.RAM_afe)) > 1) & (len(np.unique(frame.RAM_noise)) > 1): grid_lowpass2 = gridspec.GridSpecFromSubplotSpec(4, 1, subplot_spec=grid_here, height_ratios=[1, 1, 1, 0.1], hspace=0.05) # if (np.unique(frame.RAM_afe) != 0):grid_left[g] ax_external = plt_time_arrays('red', grid_lowpass2, 1, frame.RAM_afe * 100, time=time, nr=0) # if (np.unique(frame.RAM_noise) != 0): remove_xticks(ax_external) ax_intrinsic = plt_time_arrays('purple', grid_lowpass2, 1, frame.RAM_noise * 100, time=time, nr=1) ax_intrinsic.text(-0.6, 0.5, '$\%$', rotation=90, va='center', transform=ax_intrinsic.transAxes) ax_intrinsic.show_spines('l') ax_external.show_spines('l') # ax_ams.append(axt_p2) # color_timeseries = 'black' # axt_p2.set_xlabel('Time [ms]') # axt_p2.text(-0.6, 0.5, '$\%$', rotation=90, va='center', transform=axt_p2.transAxes) # ax_ams.append(axt_p2) vers = 'all' elif len(np.unique(frame.RAM_afe)) > 1: color_timeseries = 'red' nr_plot = 0 print(str(g) + ' afevar ' + str(np.var(frame.RAM_afe)) + ' afenoise ' + str(np.var(frame.RAM_noise))) try: ax_external, ff, pp, ff_am, pp_am = plot_lowpass2([grid_lowpass[nr_plot]], time, (frame.RAM_afe + frame.RAM_noise) * 100, deltat, eod_fr, color1=color_timeseries, lw=1, extract=False) except: print('add up thing') embed() ax_external.show_spines('l') ax_intrinsic = plt.subplot(grid_lowpass[1]) ax_intrinsic.show_spines('l') ax_intrinsic.axhline(0, color='black', lw=0.5) ax_intrinsic.axhline(0, color='purple', lw=0.5) remove_xticks(ax_external) remove_xticks(ax_intrinsic) join_x([ax_intrinsic, ax_external]) join_y([ax_intrinsic, ax_external]) vers = 'first' elif len(np.unique(frame.RAM_noise)) > 1: color_timeseries = 'purple' nr_plot = 1 print(str(g) + ' afevar ' + str(np.var(frame.RAM_afe)) + ' afenoise ' + str(np.var(frame.RAM_noise))) try: ax_intrinsic, ff, pp, ff_am, pp_am = plot_lowpass2([grid_lowpass[nr_plot]], time, (frame.RAM_afe + frame.RAM_noise) * 100, deltat, eod_fr, color1=color_timeseries, lw=1, extract=False) except: print('add up thing') embed() ax_external = plt.subplot(grid_lowpass[0]) ax_external.show_spines('l') ax_intrinsic.show_spines('l') ax_external.axhline(0, color='black', lw=0.5) ax_external.axhline(0, color='red', lw=0.5) join_x([ax_intrinsic, ax_external]) join_y([ax_intrinsic, ax_external]) vers = 'second' ax_intrinsic.set_yticks_delta(6) ax_external.set_yticks_delta(6) ax_external.text(-0.6, 0.5, '$\%$', va='center', rotation=90, transform=ax_external.transAxes) ax_intrinsic.text(-0.6, 0.5, '$\%$', va='center', rotation=90, transform=ax_intrinsic.transAxes) remove_xticks(ax_intrinsic) # if (len(np.unique(frame.RAM_afe)) > 1) & (len(np.unique(frame.RAM_noise)) > 1): ax_external.set_xlabel('') # remove_yticks(ax) ax_ams.append(ax_external) remove_xticks(ax_external) ax_n, ff, pp, ff_am, pp_am = plot_lowpass2([grid_lowpass[2]], time, noise_final_c, deltat, eod_fr, extract=False, color1='grey', lw=1) remove_yticks(ax_n) if g == 1: # ax_n.set_xlabel('Time [ms]', labelpad=-0.5) ax_n.text(0.5, xlabel_xpos_y_modelanddata() * 3, 'Time [ms]', transform=ax_n.transAxes, ha='center') else: remove_xticks(ax_n) ax_n.set_ylim(ylim) if vers == 'first': ax_external.text(1, 1, 'RAM(' + cs[0] + ')', ha='right', color='red', transform=ax_external.transAxes) ax_n.text(start_pos_modeldata(), 1.1, noise_component_name(), ha='right', color='gray', transform=ax_n.transAxes) elif vers == 'second': ax_external.text(1, 1, 'RAM(' + cs[1] + ')', ha='right', color='red', transform=ax_external.transAxes) ax_intrinsic.text(start_pos_modeldata(), 1.1, signal_component_name(), ha='right', color='purple', transform=ax_intrinsic.transAxes) ax_n.text(start_pos_modeldata(), 0.9, noise_component_name(), ha='right', color='gray', transform=ax_n.transAxes) else: ax_n.text(start_pos_modeldata(), 0.9, noise_component_name(), ha='right', color='gray', transform=ax_n.transAxes) ax_external.text(1, 1, 'RAM', ha='right', color='red', transform=ax_external.transAxes) ax_intrinsic.text(start_pos_modeldata(), 1.1, signal_component_name(), ha='right', color='purple', transform=ax_intrinsic.transAxes) ax_external.tick_params(axis='y', which='major', labelsize=8.4) ax_intrinsic.tick_params(axis='y', which='major', labelsize=8.4) ax_n.tick_params(axis='y', which='major', labelsize=8.4) # xtick_labelsize = # embed() return ax_ams, ax_external def start_pos_modeldata(): return 1.03 def signal_component_name(): return r'$s_\xi(t)$' #r'$\xi_{signal}$'#signal noise' def noise_component_name(): #$\xi_{noise}$noise_name = return 'Intrinsic noise' ##r'$\xi_{noise}$'#'Noise component'#'intrinsic noise' def ypos_x_modelanddata(): return -0.45 if __name__ == '__main__': model = resave_small_files("models_big_fit_d_right.csv", load_folder='calc_model_core') cells = model.cell.unique() params = {'cells': cells} show = True # if show == False: # low CV: cells = ['2012-07-03-ak-invivo-1'] plot_style() default_settings(lw=0.5, column=2, length=3.35) #8.5 redo = False D_extraction_method = ['additiv_cv_adapt_factor_scaled'] # D_extraction_method = ['additiv_visual_d_4_scaled'] ########################## # hier printen wir die table Werte zum kopieren in den Text path = 'print_table_suscept-model_params_suscept_table.csv' if os.path.exists(path): table = pd.read_csv(path) table_printen(table) path = 'print_table_all-model_params_suscept_table.csv' if os.path.exists(path): table = pd.read_csv() print('model big') table_printen(table) #embed() ########################## #embed() print('hi') model_and_data2(eod_metrice=False, width=0.005, D_extraction_method=D_extraction_method) #r'$\frac{1}{mV^2S}$'