benda/nonlinearbaseline2025
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added correlation with baseline rate

Browse Source
This commit is contained in:
Jan Benda 2025-05-26 18:11:08 +02:00
parent 04578fd77e
commit d56c0c7c7e
4 changed files with 67 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
dataoverview.py
View File

@ -117,12 +117,12 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
ha='right', va='center', fontsize='small')
# statistics:
r, p = pearsonr(xdata, ydata)
ax.text(1, 0.9, f'$R={r:.2f}$ **', ha='right',
ax.text(1, 0.9, f'$R={r:.2f}$', ha='right',
transform=ax.transAxes, fontsize='small')
ax.text(1, 0.77, f'{significance_str(p)}', ha='right',
transform=ax.transAxes, fontsize='small')
#ax.text(1, 0.77, f'{significance_str(p)}', ha='right',
# transform=ax.transAxes, fontsize='small')
if 'cvbase' in xcol:
ax.text(1, 0.77, f'$n={data.rows()}$', ha='right',
ax.text(1, 0.64, f'$n={data.rows()}$', ha='right',
transform=ax.transAxes, fontsize='small')
print(f' correlation {xcol:<8s} - {ycol}: r={r:5.2f}, p={p:.2g}')
return cax
@ -154,6 +154,15 @@ def si_stats(title, data, sicol, si_thresh, nsegscol):
contrasts = 100*data['contrast']
print(' contrasts: ', ' '.join([f'{c:.2g}%' for c in np.unique(contrasts)]))
print(f' contrasts: {np.min(contrasts):.2g}% - {np.max(contrasts):.2g}%, median={np.median(contrasts):.2g}%, mean={np.mean(contrasts):.2g}%, std={np.std(contrasts):.2g}%')
cols = ['cvbase', 'respmod2', 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim']
for i in range(len(cols)):
for j in range(i + 1, len(cols)):
xcol = cols[i]
ycol = cols[j]
if xcol not in data or ycol not in data:
continue
r, p = pearsonr(data[xcol], data[ycol])
print(f' correlation {xcol:<11s} - {ycol:<11s}: r={r:5.2f}, p={p:.5f}')
def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
@ -204,6 +213,19 @@ def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
cax.set_ylabel('CV$_{\\rm base}$')
def plot_rate_si_punit(ax, data, ycol, si_thresh, color):
ax.set_xlabel('Baseline rate $r$', 'Hz')
ax.set_xlim(0, 700)
ax.set_xticks_delta(200)
ax.set_ylabel('SI($r$)')
ax.set_ylim(0, 6.5)
ax.set_yticks_delta(2)
examples = punit_examples if 'stimindex' in data else model_examples
cax = plot_corr(ax, data, 'ratebase', ycol, 'cvbase', 0, 1.5, 0.016,
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('CV$_{\\rm base}$')
def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_xlabel('CV$_{\\rm base}$')
ax.set_xlim(0, 0.2)
@ -248,6 +270,19 @@ def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
cax.set_ylabel('CV$_{\\rm base}$')
cax.set_yticks_delta(0.1)
def plot_rate_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_xlabel('Baseline rate $r$', 'Hz')
ax.set_xlim(50, 200)
ax.set_xticks_delta(50)
ax.set_ylabel('SI($r$)')
ax.set_ylim(0, 10)
ax.set_yticks_delta(2)
cax = plot_corr(ax, data, 'ratebase', ycol, 'cvbase', 0, 0.2, 0.06,
'coolwarm', color, si_thresh, *ampul_examples)
cax.set_ylabel('CV$_{\\rm base}$')
cax.set_yticks_delta(0.1)
if __name__ == '__main__':
punit_model = TableData(data_path /
@ -304,8 +339,10 @@ if __name__ == '__main__':
s.model_color2)
plot_rmod_si_punit(axs[0, 1], punit_model, 'dsinorm100', si_thresh,
s.model_color2)
plot_cvstim_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh,
s.model_color2)
#plot_cvstim_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh,
# s.model_color2)
plot_rate_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh,
s.model_color2)
print()
si_stats('P-unit data:', punit_data, 'sinorm' + si, si_thresh,
@ -316,8 +353,10 @@ if __name__ == '__main__':
s.punit_color2)
plot_rmod_si_punit(axs[1, 1], punit_data, 'sinorm' + si, si_thresh,
s.punit_color2)
plot_cvstim_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh,
s.punit_color2)
#plot_cvstim_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh,
# s.punit_color2)
plot_rate_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh,
s.punit_color2)
print()
si_stats('Ampullary data:', ampul_data, 'sinorm' + si, si_thresh,
@ -328,8 +367,10 @@ if __name__ == '__main__':
s.ampul_color2)
plot_rmod_si_ampul(axs[2, 1], ampul_data, 'sinorm' + si, si_thresh,
s.ampul_color2)
plot_cvstim_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh,
s.ampul_color2)
#plot_cvstim_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh,
# s.ampul_color2)
plot_rate_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh,
s.ampul_color2)
print()
fig.common_xticks(axs[:2, 0])

10
modelsusceptcontrasts.py
View File

@ -30,20 +30,21 @@ def load_chi2(file_path, cell_name, contrast=None, n=None):
return freqs, chi2, fcutoff, contrast, n
def plot_chi2_contrasts(axs, s, cell_name):
def plot_chi2_contrasts(axs, s, cell_name, vmax):
d = sims_path / f'{cell_name}-baseline.npz'
data = np.load(d)
rate = float(data['rate'])
cv = float(data['cv'])
print(f' {cell_name}: r={rate:3.0f}Hz, CV={cv:4.2f}')
freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name)
cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate)
cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate, vmax)
cax.set_ylabel('')
axs[0].set_title(r'$c$=0\,\%', fontsize='medium')
for k, alpha in enumerate([0.01, 0.03, 0.1]):
freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name,
alpha)
cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate)
cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate, vmax)
vmax /= 2
if alpha < 0.1:
cax.set_ylabel('')
axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium')
@ -142,8 +143,9 @@ if __name__ == '__main__':
fig.subplots_adjust(leftm=7, rightm=9, topm=2, bottomm=3.5,
wspace=1, hspace=0.5)
print('Example cells:')
vmax = [2, 6, 8, 40]
for k in range(len(model_cells)):
plot_chi2_contrasts(axs[k], s, model_cells[k])
plot_chi2_contrasts(axs[k], s, model_cells[k], vmax[k])
for k in range(4):
fig.common_yticks(axs[k, :])
fig.common_xticks(axs[:4, k])

14
modelsusceptlown.py
View File

@ -14,7 +14,7 @@ data_path = Path('data')
sims_path = data_path / 'simulations'
def plot_chi2_contrasts(axs, s, cell_name, nsegs=None):
def plot_chi2_contrasts(axs, s, cell_name, nsegs=None, vmax=None):
d = sims_path / f'{cell_name}-baseline.npz'
data = np.load(d)
rate = float(data['rate'])
@ -23,14 +23,18 @@ def plot_chi2_contrasts(axs, s, cell_name, nsegs=None):
freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name,
None, nsegs)
ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$'
cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate)
cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate, vmax)
cax.set_ylabel('')
axs[0].set_title(f'$c$=0\\,\\%, {ns}', fontsize='medium')
for k, alpha in enumerate([0.01, 0.03, 0.1]):
freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name,
alpha, nsegs)
ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$'
cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate)
cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate, vmax)
if n < 1000:
vmax /= 10
else:
vmax /= 4
if alpha < 0.1:
cax.set_ylabel('')
axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%, {ns}',
@ -139,8 +143,8 @@ if __name__ == '__main__':
for ax in axs.flat:
ax.set_visible(False)
print('Example cells:')
plot_chi2_contrasts(axs[0], s, model_cell)
plot_chi2_contrasts(axs[1], s, model_cell, nsmall)
plot_chi2_contrasts(axs[0], s, model_cell, None, 40)
plot_chi2_contrasts(axs[1], s, model_cell, nsmall, 600)
for k in range(2):
fig.common_yticks(axs[k, :])
for k in range(4):

2
modelsusceptovern.py
View File

@ -75,7 +75,7 @@ def plot_overn(ax, s, files, nmax=1e6, title=False):
def plot_chi2_overn(axs, s, cell_name):
print(cell_name)
files, nums = noise_files(sims_path, cell_name)
for k, nsegs in enumerate([1e1, 1e2, 1e3, 1e6]):
for k, nsegs in enumerate([1e2, 1e3, 1e4, 1e6]):
freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name,
None, nsegs)
ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$'