susceptibility1/nonlin_regime.py

470 lines
25 KiB
Python

import os
import sys
import numpy as np
import pandas as pd
from matplotlib import gridspec, pyplot as plt
from plotstyle import plot_style
from plottools.colors import lighter, darker
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
from threefish.load import load_savedir, resave_small_files, save_visualization
from threefish.plot.limits import 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_fbasename_small, label_sum, \
remove_yticks
from threefish.RAM.plot_subplots import colors_suscept_paper_dots, plt_spikes_ROC, recalc_fr_to_DF1
from threefish.RAM.values import vals_model_full
from threefish.twobeat.calc_model import calc_threefish_model
from threefish.twobeat.colors import colors_susept
from threefish.twobeat.reformat import c_dist_recalc_func, c_dist_recalc_here, dist_recalc_phaselockingchapter, \
find_dfs
from threefish.twobeat.subplots import plt_psd_saturation, plt_single_trace, plt_stim_saturation, power_spectrum_name
def nonlin_regime0():
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
c_nrs_orig = [0.002, 0.01, 0.05, 0.1]
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, panel E
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 = [[220, 220, 220, 220]])
##############################
# 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 = []
xlimp = (0, 300)
n = 1
dev = 0.001
dev_name = ['05']
a_fr = 1
a = 0
for c_nn, c_nr in enumerate(c_nrs):
a_f2s = [c_nrs_orig[c_nn]]
#################################
# hier wurden die psds vorgespeichert weil die lange brauchen
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]]
# v_mems, arrays_spikes, arrays_stim, frame_diff, position_diff, arrays_convovled, names, p_array, ff
_, spikes_arrays_all, stimulus_arrays_all, _, _, convovled_arrays, threefish_names, _, _ = calc_threefish_model(
results_diff, position_diff, model_cells, cell_here, trials_nr=trials_nr,
stimulus_length=stimulus_length_here, freq1=[freq1 + eodf], freq2=[freq2 + eodf], a_f2s=a_f2s, a_f1s=[
c_nrs_orig[
c_nn]], a_fr=a_fr, dev=dev, dev_name=dev_name, n=n, reshuffled=reshuffled, nfft=nfft_here,
psd_min_amps='_minamps_', mean_choice='', printing_cell=False)#means_different
##################################################################
# plot panels A, B, C, D
time = np.arange(0, len(convovled_arrays[a][0]) / sampling, 1 / sampling)
time = time * 1000
# die arrays die rausgegeben werden sind für das threefish problem und hier nehmen wir nur den 012, also wo alle fische dabei sind
stimuli_arrays = [stimulus_arrays_all[3]] #1:: [stimulus_arrays_all[1],stimulus_arrays_all[3]]#
spike_arrays = [spikes_arrays_all[3]] # [spikes_arrays_all[1],spikes_arrays_all[3]]#spikes_arrays_all[1::]
p_arrays_all.append(p_arrays_here)
for a in range(len(stimuli_arrays)):
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(stimuli_arrays, axe_all, axes, c_nn, c_nrs, ['grey'], grid_pt, i, time, xlim)
#############################
# plot spikes
plt_spikes_nonlin(a, spike_arrays, axss, c_nn, grid_pt, xlim)
#############################
# preallocate psd ax
axp = plt.subplot(grid_pt[-1])
axps.append(axp)
f_counter += 1
####################
# plot psds, panels A, B, C, D, bottom
plt_psd_nonlin(axps, c_nrs, color0, color01, color012, color01_2, color02,
['grey'], ff_p, frame_cell_orig, freq1, freq2, log, p_arrays_all, xlimp)
###################
# limits
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=True)
def plt_spikes_nonlin(a, arrays_sp, 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)
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, [r'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]))) + r'\%' # +'$'
beat_here = f'$\\rm{{Contrast}}={c_nrs[c_nn]:.2g}\\%$'
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(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(p_arrays_all[c_nn])):
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]')
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, arrays, names, p_arrays_p, ff_p = calc_threefish_model(results_diff,
position_diff,
model_cells,
cell_here,
trials_nr=trials_nr,
stimulus_length=stimulus_length_here,
freq1=[
freq1 + eodf],
freq2=[
freq2 + eodf],
a_f2s=a_f2s,
a_f1s=[
c_nrs_orig[
c_nn]],
a_fr=a_fr,
dev=dev,
dev_name=dev_name,
n=n,
reshuffled=reshuffled,
nfft=nfft_here,
psd_min_amps='_minamps_',
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 = [[220, 220, 220, 220]]):
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_u1.set_ylim(0, 220)
cc = np.arange(0, 8, 0.1)
ax_u1.plot(cc, 1.5 + 33*cc**2, color=darker(colors[2], 0.9), lw=0.5, zorder=100)
ax_u1.plot(cc, 1.5 + 15*cc**2, color=darker(colors[3], 0.8), lw=0.5, zorder=100)
ax_u1.plot(cc, 1.5 + 6.5*cc, color=lighter(colors[0], 0.5), lw=0.5, zorder=80)
ax_u1.plot(cc, 25 + 125*cc, color=lighter(colors[1], 0.5), lw=0.5, zorder=80)
# 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(r'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()