susceptibility1/nonlin_regime.py
2024-06-18 17:44:19 +02:00

453 lines
24 KiB
Python

import os
import sys
import numpy as np
import pandas as pd
from IPython import embed
from matplotlib import gridspec, pyplot as plt
from plotstyle import plot_style
from threefish.calc_time import extract_am
from threefish.core import find_code_vs_not, find_folder_name, find_project_data, info
from threefish.defaults import default_figsize, default_ticks_talks
from threefish.load import load_savedir, resave_small_files, save_visualization
from threefish.plot.limits import join_x, join_y, set_same_ylim
from threefish.RAM.calc_fft import log_calc_psd
from threefish.RAM.calc_model import chose_old_vs_new_model
from threefish.RAM.core import find_load_function
from threefish.RAM.plot_labels import label_deltaf1, label_deltaf2, label_diff, label_f_eod_name_core_rm, \
label_fbasename_small, label_sum, \
onebeat_cond, \
remove_yticks
from threefish.RAM.plot_subplots import colors_suscept_paper_dots, plt_spikes_ROC, recalc_fr_to_DF1
from threefish.RAM.values import val_cm_to_inch, vals_model_full
from threefish.twobeat.calc_model import calc_roc_amp_core_cocktail
from threefish.twobeat.colors import colors_susept, twobeat_cond
from threefish.twobeat.labels import f_stable_name, f_vary_name
from threefish.twobeat.reformat import c_dist_recalc_func, c_dist_recalc_here, dist_recalc_phaselockingchapter, \
find_dfs, \
get_frame_cell_params
from threefish.twobeat.subplots import plt_psd_saturation, plt_single_trace, plt_stim_saturation, plt_vmem_saturation, \
power_spectrum_name
from threefish.values import values_nfft_full_model, values_stimuluslength_model_full
def nonlin_regime0(yposs=[450, 450, 450], show=True):
results_diff = pd.DataFrame()
position_diff = 0
plot_style()
default_figsize(column=2, length=5.5)
########################################
# für das model_full, die Freuqnezen
DF1_frmult, DF2_frmult = vals_model_full(val=0.30833333333333335)
frame_cvs = resave_small_files('csv_model_data.csv', load_folder = 'calc_base')
frame_cell = frame_cvs[frame_cvs.cell == '2012-07-03-ak-invivo-1']
for d in range(len(DF1_frmult)):
DF2_frmult[d] = recalc_fr_to_DF1(DF2_frmult, d, frame_cell.fr_data.iloc[0])
DF1_frmult[d] = recalc_fr_to_DF1(DF1_frmult, d, frame_cell.fr_data.iloc[0])
freq1 = DF1_frmult[3]
freq2 = DF2_frmult[3]
# sachen die ich variieren will
###########################################
auci_wo = []
auci_w = []
cells_here = ["2013-01-08-aa-invivo-1"]
for cell_here in cells_here:
c_nrs_orig = [0.01,0.04, 0.1, 0.2] # 0.0002, 0.05, 0.5
redo = False # True
log = 'log' # 'log'
grid0 = gridspec.GridSpec(1, 1, bottom=0.13, top=0.94, left=0.11,
right=0.95, wspace=0.04) #
grid00 = gridspec.GridSpecFromSubplotSpec(2, 1,
wspace=0.04, hspace=0.75,
subplot_spec=grid0[0], height_ratios=[1.5,1],) # height_ratios=[2,1],
############################
# overview
grid_down = gridspec.GridSpecFromSubplotSpec(1, 1,
hspace=0.75,
wspace=0.1,
subplot_spec=grid00[1]) # 1.2hspace=0.4,wspace=0.2,len(chirps)
ax_u1, c_dist_recalc, c_nrs, color0, color01, color012, color01_2, color02, eodf, f_counter, frame_cell, frame_cell_orig, freq1, freq2, i, sampling = plt_amplitudes_for_contrasts(
c_nrs_orig, cell_here, freq1, freq2, grid_down, yposs)
##############################
# single examples
grid_up = gridspec.GridSpecFromSubplotSpec(1, len(c_nrs_orig),
hspace=0.75,
wspace=0.25,
subplot_spec=grid00[
0]) # height_ratios=[1, 1, 0.7], 1.2hspace=0.4,wspace=0.2,len(chirps)
start = 200 # 1000
mults_period = 3
xlim = [start, start + (mults_period * 1000 / np.min([np.abs(freq1), np.abs(freq2)]))]
axts = []
axps = []
axes = []
axss = []
p_arrays_all = []
model_fit = ''#'_old_fit_' # ''#'_old_fit_'#''#'_old_fit_'#''#'_old_fit_'#''###'_old_fit_'
model_cells, reshuffled = chose_old_vs_new_model(model_fit=model_fit)
axe_all = []
for c_nn, c_nr in enumerate(c_nrs):
#################################
# hier wurden die psds vorgespeichert weil die lange brauchen
a_f2s = [c_nrs_orig[c_nn]]
xlimp = (0, 300)
runs = 1
n = 1
dev = 0.001
dev_name = ['05']
a_fr = 1
a = 0
auci_w, auci_wo, ff_p, p_arrays_here, position_diff, results_diff = load_psds_for_nonlin_regime(a_f2s, a_fr,
auci_w,
auci_wo,
c_nn,
c_nrs_orig,
cell_here,
dev,
dev_name,
eodf,
freq1,
freq2,
model_cells,
n,
position_diff,
redo,
reshuffled,
results_diff,
xlimp)
# #
#
############################
# spikes werden jedes mal neu generiert weil das schnell geht
trials_nr = 10
stimulus_length_here = 1 # values_stimuluslength_model_full()
nfft_here = stimulus_length_here * 20000
a_f2s = [c_nrs_orig[c_nn]]
_, arrays_spikes, arrays_stim, _, _, _, _, arrays, names, _, _ = calc_roc_amp_core_cocktail([freq1 + eodf],
[freq2 + eodf],
auci_wo, auci_w,
results_diff,
a_f2s,
trials_nr,
nfft_here, a_fr,
stimulus_length_here,
model_cells,
position_diff,
dev, cell_here,
a_f1s=[
c_nrs_orig[
c_nn]],
dev_name=dev_name,
min_amps='_minamps_',
n=n,
reshuffled=reshuffled,
mean_choice='',
printing_cell=False)
##################################################################
# preprocessing der arrays
time = np.arange(0, len(arrays[a][0]) / sampling, 1 / sampling)
time = time * 1000
arrays_time = [arrays[3]] # [v_mems[1],v_mems[3]]#[1,2]#[1::]
arrays_here = [arrays[3]] # [arrays[1],arrays[3]]#arrays[1::]#
arrays_st = [arrays_stim[3]] #1:: [arrays_stim[1],arrays_stim[3]]#
arrays_sp = [arrays_spikes[3]] # [arrays_spikes[1],arrays_spikes[3]]#arrays_spikes[1::]
colors_array_here = ['grey', 'grey', 'grey'] # colors_array[1::]
p_arrays_all.append(p_arrays_here)
for a in range(len(arrays_here)):
grid_pt = gridspec.GridSpecFromSubplotSpec(3, 1,
hspace=0.3,
wspace=0.2,
subplot_spec=grid_up[a, c_nn], height_ratios = [1.2, 1,2]
)
########################################
# plot stimulus
plt_stim_nonlin(arrays_st, axe_all, axes, c_nn, c_nrs, colors_array_here, grid_pt, i, time, xlim)
#############################
# plot spikes
plt_spikes_nonlin(a, arrays_sp, axps, axss, c_nn, grid_pt, xlim)
f_counter += 1
####################
# plot psds
a = plt_psd_nonlin(a, arrays_here, axps, c_nrs, color0, color01, color012, color01_2, color02,
colors_array_here, ff_p, frame_cell_orig, freq1, freq2, log, p_arrays_all, xlimp)
axps[0].get_shared_y_axes().join(*axps)
axps[0].get_shared_x_axes().join(*axps)
set_same_ylim(axe_all)
set_same_ylim(axps)
fig = plt.gcf()
fig.tag([axes[0], axes[1], axes[2], axes[3], ax_u1], xoffs=-2.3, yoffs=1.4)
save_visualization(cell_here, show)
def plt_spikes_nonlin(a, arrays_sp, axps, axss, c_nn, grid_pt, xlim):
axs = plt.subplot(grid_pt[1])
plt_spikes_ROC(axs, 'grey', np.array(arrays_sp[a], dtype=object) * 1000, xlim, lw=1)
if c_nn == 0:
axs.xscalebar(0.15, -0.1, 30, 'ms', va='right', ha='bottom')
axss.append(axs)
axp = plt.subplot(grid_pt[-1])
axps.append(axp)
def plt_stim_nonlin(arrays_st, axe_all, axes, c_nn, c_nrs, colors_array_here, grid_pt, i, time, xlim):
axe = plt.subplot(grid_pt[0])
axes.append(axe)
am, time_am = extract_am(arrays_st[0], time, extract='', norm=False)
am = am - np.mean(am)
if c_nn == 0:
spines = 'l'
else:
spines = ''
plt_stim_saturation(0, [], [am * 100], axe, colors_array_here, i,
c_nn, ['Contrast [$\%$]'], time,
xlim=xlim, lw=1, spines=spines) # np.array(arrays_sp)*1000
beat_here = r'$\rm{Contrast}=%s$' % (int(np.round(c_nrs[c_nn]))) + '$\%$' # +'$'
title_name = beat_here # fish + '\n' + +c1+c2#twobeat_cond(big=True, double=True,cond=False)
axe.text(1, 1.1, title_name, va='bottom', ha='right',
transform=axe.transAxes)
axe_all.append(axe)
def plt_psd_nonlin(a, arrays_here, axps, c_nrs, color0, color01, color012, color01_2, color02, colors_array_here, ff_p,
frame_cell_orig, freq1, freq2, log, p_arrays_all, xlimp):
pps = []
for c_nn, c_nr in enumerate(c_nrs):
for a in range(len(arrays_here)):
axps_here = [[axps[0], axps[1], axps[2], axps[3]]] # [axps[3], axps[4], axps[5]
axp = axps_here[a][c_nn]
pp = log_calc_psd(log, p_arrays_all[c_nn][a][0],
np.nanmax(p_arrays_all))
pps.append(pp)
colors_peaks = [color01, color02, color012, color01_2, color0]
markeredgecolors = [color01, color02, color012, color01_2, color0]
fr = frame_cell_orig.fr.unique()[0]
labels = [label_deltaf1(), label_deltaf2(),
label_sum(), label_diff(), label_fbasename_small(), label_fbasename_small()]
freqs = [np.abs(freq1), np.abs(freq2), np.abs(freq1) + np.abs(freq2), np.abs(np.abs(freq1) - np.abs(freq2)),
fr]
plt_psd_saturation(pp, ff_p, a, axp, colors_array_here, freqs=freqs,
colors_peaks=colors_peaks, xlim=xlimp,
markeredgecolor=markeredgecolors, labels=labels, log=log, add_log=5)
if c_nn == 0:
axp.legend(ncol=5, loc=(-0, -1.1))
if log:
scalebar = False
if scalebar:
axp.show_spines('b')
if c_nn == 0:
axp.yscalebar(-0.05, 0.5, 20, 'dB', va='center', ha='left')
axp.set_ylim(-33, 5)
else:
axp.show_spines('lb')
if c_nn == 0:
axp.set_ylabel('dB') # , va='center', ha='left'
else:
remove_yticks(axp)
axp.set_ylim(-39, 5)
else:
axp.show_spines('lb')
if c_nn != 0:
remove_yticks(axp)
else:
axp.set_ylabel(power_spectrum_name())
axp.set_xlabel('Frequency [Hz]')
return a
def load_psds_for_nonlin_regime(a_f2s, a_fr, auci_w, auci_wo, c_nn, c_nrs_orig, cell_here, dev,
dev_name, eodf, freq1, freq2, model_cells, n,
position_diff, redo, reshuffled,
results_diff, xlimp):
trials_nr = 1
stimulus_length_here = 500 # 500 # values_stimuluslength_model_full()
nfft_here = 20 * 20000
save_dir = load_savedir(level=0).split('/')[0] + '_afe_' + str(a_f2s[0]) + '_nfft_' + str(
nfft_here) + '_len_' + str(stimulus_length_here)
name_psd = find_project_data() + save_dir + '_psd.npy'
name_psd_f = find_project_data() + save_dir + '_psdf.npy'
if ((not os.path.exists(name_psd)) | (redo == True)): # | (do == True)
auci_w, auci_wo, ff_p, p_arrays_here, position_diff, results_diff = develop_save_nonlin_psd_values(a_f2s, a_fr,
auci_w,
auci_wo,
c_nn,
c_nrs_orig,
cell_here,
dev,
dev_name,
eodf,
freq1, freq2,
model_cells,
n, name_psd,
name_psd_f,
nfft_here,
position_diff,
reshuffled,
results_diff,
stimulus_length_here,
trials_nr,
xlimp)
else:
ff_p = np.load(name_psd_f) # p_arrays_p
p_arrays_here = np.load(name_psd) # p_arrays_p
return auci_w, auci_wo, ff_p, p_arrays_here, position_diff, results_diff
def develop_save_nonlin_psd_values(a_f2s, a_fr, auci_w, auci_wo, c_nn, c_nrs_orig, cell_here, dev,
dev_name, eodf, freq1, freq2, model_cells, n, name_psd,
name_psd_f, nfft_here, position_diff, reshuffled,
results_diff, stimulus_length_here, trials_nr, xlimp):
# psd generierung
_, _, arrays_stim, results_diff, position_diff, auci_wo, auci_w, arrays, names, p_arrays_p, ff_p = calc_roc_amp_core_cocktail(
[freq1 + eodf], [freq2 + eodf], auci_wo, auci_w, results_diff, a_f2s, trials_nr, nfft_here, a_fr,
stimulus_length_here, model_cells, position_diff, dev, cell_here,a_f1s=[c_nrs_orig[c_nn]], dev_name=dev_name, min_amps='_minamps_', n=n, reshuffled=reshuffled, mean_choice='first')
p_arrays_here = [p_arrays_p[3]]
# embed()
for p in range(len(p_arrays_here)):
p_arrays_here[p][0] = p_arrays_here[p][0][ff_p < xlimp[1]]
ff_p = ff_p[ff_p < xlimp[1]]
np.save(name_psd, p_arrays_here)
np.save(name_psd_f, ff_p)
return auci_w, auci_wo, ff_p, p_arrays_here, position_diff, results_diff
def plt_amplitudes_for_contrasts(c_nrs_orig, cell_here, freq1, freq2, grid_down, yposs):
full_names = [
'calc_model_amp_freqs-_F1_0.22833333333333333Fr_F2_1Fr_af_coupled__C1Len_100_FirstC1_0.0001_LastC1_0.3_FRrelavtiv__start_0.0001_end_0.3_StimLen_100_nfft_20000_trialsnr_1_mult_minimum_1_power_1_minamps__dev_original_05_not_log__point_1_fft2_temporal']
for i, full_name in enumerate(full_names):
frame = load_nonlin_regime_values(full_name)
frame_cell_orig = frame[(frame.cell == cell_here)].copy()
if len(frame_cell_orig) > 0:
f_counter = 0
frame_cell_orig, df1s, df2s, f1s, f2s = find_dfs(frame_cell_orig)
eodf = frame_cell_orig.f0.unique()[0]
#######################################################################################
# übersicht
frame_cell = frame_cell_orig[
(frame_cell_orig.df1 == freq1) & (frame_cell_orig.df2 == freq2)]
if len(frame_cell) < 1:
freq1 = frame_cell_orig.iloc[(np.argmin(np.abs(frame_cell_orig.df1 - freq1)))].df1
freq2 = frame_cell_orig.iloc[(np.argmin(np.abs(frame_cell_orig.df2 - freq2)))].df2
frame_cell = frame_cell_orig[
(frame_cell_orig.df1 == freq1) & (frame_cell_orig.df2 == freq2)]
# embed()
labels, alpha, color01, color01_012, color02, color02_012, colors, colors_array, linestyles, scores, linewidths = colors_susept(
add='_mean_original', nr=4)
sampling = 20000
c_dist_recalc = dist_recalc_phaselockingchapter()
c_nrs = c_dist_recalc_func(frame_cell, c_nrs=c_nrs_orig, cell=cell_here,
c_dist_recalc=c_dist_recalc)
if not c_dist_recalc:
c_nrs = np.array(c_nrs) * 100
letters = ['A', 'B', 'C', 'D']
scores = ['c_B1_012_original', 'c_B2_012_original', 'c_B1+B2_012_original', 'c_B1-B2_012_original']
color01, color012, color01_2, color02, color0_burst, color0 = colors_suscept_paper_dots()
colors = [color01, color02, color012, color01_2]
linestyles = ['-', '-', '-', '-']
index = [0, 1, 2, 3]
ax_u1 = plt.subplot(grid_down[0, i])
labels = [label_deltaf1(), label_deltaf2(),
label_sum(), label_diff(), label_fbasename_small()]
# ax_u1.legend(ncol = 4, loc = (0,1.2))
plt_single_trace([], ax_u1, frame_cell_orig, freq1, freq2,
scores=np.array(scores)[index], labels=np.array(labels)[index],
colors=np.array(colors)[index],
linestyles=np.array(linestyles)[index],
linewidths=np.array(linewidths)[index],
alpha=np.array(alpha)[index],
thesum=False, B_replace='F', default_colors=False,
c_dist_recalc=c_dist_recalc)
# ax_us.append(ax_u1)
frame_cell = frame_cell_orig[
(frame_cell_orig.df1 == freq1) & (frame_cell_orig.df2 == freq2)]
c1 = c_dist_recalc_here(c_dist_recalc, frame_cell)
ax_u1.set_xlim(0, 25)
if i != 0:
ax_u1.set_ylabel('')
remove_yticks(ax_u1)
ax_u1.scatter(c_nrs, (np.array(yposs[i]) - 0) * np.ones(len(c_nrs)), color='black',
marker='v',
clip_on=False)
#
for c_nn, c_nr in enumerate(c_nrs):
ax_u1.text(c_nr, yposs[i][c_nn] + 30, letters[c_nn], color='black', ha='center',
va='top')
# ax_u1.plot([c_nr, c_nr], [0, 435], color='black', linewidth=0.8, clip_on=False)
ylim = ax_u1.get_ylim()
ax_u1.set_ylim(0, ylim[-1])
ax_u1.set_xlabel('Contrast [$\%$]')
ax_u1.set_ylabel('Amplitude [Hz]')
return ax_u1, c_dist_recalc, c_nrs, color0, color01, color012, color01_2, color02, eodf, f_counter, frame_cell, frame_cell_orig, freq1, freq2, i, sampling
def load_nonlin_regime_values(full_name):
version_comp, subfolder, mod_name_slash, mod_name, subfolder_path = find_code_vs_not(
)
if os.path.exists(find_folder_name('calc_cocktailparty') + '/' + full_name + '.csv') & (version_comp != 'public'):
frame = pd.read_csv(find_folder_name('calc_cocktailparty') + '/' + full_name + '.csv')
load_function = find_load_function()
frame.to_csv(find_project_data() + load_function + '-traces.csv')
else:
load_function = find_load_function()
frame = pd.read_csv(
find_project_data() + load_function + '-traces.csv') # resave_small_files(full_name + '.csv', load_folder='calc_cocktailparty')
return frame
if __name__ == '__main__':
sys.excepthook = info
nonlin_regime0(yposs = [[220, 220, 220, 220]],)#, [430,470], [200,200], [200,200, 200, 200], [200,200, 200, 200], [200,200, 200, 200]