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

finished data overview

Browse Source
This commit is contained in:
Jan Benda
2025-06-01 23:09:09 +02:00
parent 6e0098a273
commit f2ee764f78
4 changed files with 232 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

206
dataoverview.py
View File

@@ -92,7 +92,7 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
cb.outline.set_linewidth(0)
# pdf x-axis:
kde = gaussian_kde(xdata, 0.02*xmax/np.std(xdata, ddof=1))
xx = np.linspace(0, ax.get_xlim()[1], 400)
xx = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 400)
pdf = kde(xx)
xax = ax.inset_axes([0, 1.05, 1, 0.2])
xax.show_spines('')
@@ -102,7 +102,7 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
xax.set_ylim(0, xpdfmax)
# pdf y-axis:
kde = gaussian_kde(ydata, 0.02*ymax/np.std(ydata, ddof=1))
xx = np.linspace(0, ax.get_ylim()[1], 400)
xx = np.linspace(ax.get_ylim()[0], ax.get_ylim()[1], 400)
pdf = kde(xx)
yax = ax.inset_axes([1.05, 0, 0.2, 1])
yax.show_spines('')
@@ -124,10 +124,106 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
if 'cvbase' in xcol:
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}')
print(f' correlation {xcol:<8s} - {ycol}: r={r:5.2f}, p={p:.2e}')
return cax
def plot_corr_contrast(ax, data, xcol, ycol, xpdfmax, color,
si_thresh, example=[], split_example=[], examples=[]):
xdata = data[xcol]
ydata = data[ycol]
ax.axhline(si_thresh, color='k', ls=':', lw=0.5)
xmax = ax.get_xlim()[1]
ymax = ax.get_ylim()[1]
mask = (xdata < xmax) & (ydata < ymax)
if 'stimindex' in data:
for cell, run in example + split_example + examples:
mask &= ~((data['cell'] == cell) & (data['stimindex'] == run))
else: # simulations
for cell, alpha in example + split_example + examples:
mask &= ~((data['cell'] == cell) & (data['contrast'] == alpha))
contrasts = [[r'1\,\%', 0, 0.015, s.psC1],
[r'3\,\%', 0.015, 0.035, s.psC3],
[r'5\,\%', 0.035, 0.07, s.psC5],
[r'10\,\%', 0.07, 0.15, s.psC10],
[r'20\,\%', 0.15, 0.25, s.psC20]]
for l, c0, c1, sc in contrasts:
cmask = (data['contrast'] > c0) & (data['contrast'] < c1)
label = f'{l} ({np.sum(cmask)})'
cmask = cmask & mask
#if np.sum(cmask) > 0:
# ax.plot(xdata[cmask], ydata[cmask], label=l, clip_on=False,
# alpha=0.6, zorder=np.random.randint(0, 50), **sc)
for i in np.where(cmask)[0]:
ax.plot(xdata[i], ydata[i], label=label, clip_on=False,
alpha=0.7, zorder=np.random.randint(0, 50), **sc)
label = None
elw = 0.3
for l, c0, c1, sc in contrasts:
cmask = (data['contrast'] >= c0) & (data['contrast'] < c1)
if 'stimindex' in data:
for cell, run in example:
mask = cmask & (data['cell'] == cell) & (data['stimindex'] == run)
ax.plot(xdata[mask], ydata[mask], zorder=50, ms=4, marker='^',
color=sc['color'], mew=elw, mec='black', clip_on=False)
for cell, run in split_example:
mask = cmask & (data['cell'] == cell) & (data['stimindex'] == run)
ax.plot(xdata[mask], ydata[mask], zorder=51, ms=4, marker='s',
color=sc['color'], mew=elw, mec='black', clip_on=False)
for cell, run in examples:
mask = cmask & (data['cell'] == cell) & (data['stimindex'] == run)
ax.plot(xdata[mask], ydata[mask], zorder=52, ms=4, marker='o',
color=sc['color'], mew=elw, mec='black', clip_on=False)
else: # simulations
for cell, alpha in example:
mask = cmask & (data['cell'] == cell) & (data['contrast'] == alpha)
ax.plot(xdata[mask], ydata[mask], zorder=50, ms=4, marker='^',
color=sc['color'], mew=elw, mec='black', clip_on=False)
for cell, alpha in split_example:
mask = cmask & (data['cell'] == cell) & (data['contrast'] == alpha)
ax.plot(xdata[mask], ydata[mask], zorder=51, ms=4, marker='s',
color=sc['color'], mew=elw, mec='black', clip_on=False)
for cell, alpha in examples:
mask = cmask & (data['cell'] == cell) & (data['contrast'] == alpha)
ax.plot(xdata[mask], ydata[mask], zorder=52, ms=4, marker='o',
color=sc['color'], mew=elw, mec='black', clip_on=False)
# pdf x-axis:
kde = gaussian_kde(xdata, 0.02*xmax/np.std(xdata, ddof=1))
xx = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 400)
pdf = kde(xx)
xax = ax.inset_axes([0, 1.05, 1, 0.2])
xax.show_spines('')
xax.fill_between(xx, pdf, facecolor=color, edgecolor='none')
#xax.plot(xx, np.zeros(len(xx)), clip_on=False, color=color, lw=0.5)
xax.set_ylim(bottom=0)
xax.set_ylim(0, xpdfmax)
# pdf y-axis:
kde = gaussian_kde(ydata, 0.02*ymax/np.std(ydata, ddof=1))
xx = np.linspace(ax.get_ylim()[0], ax.get_ylim()[1], 400)
pdf = kde(xx)
yax = ax.inset_axes([1.05, 0, 0.17, 1])
yax.show_spines('')
yax.fill_betweenx(xx, pdf, facecolor=color, edgecolor='none')
#yax.plot(np.zeros(len(xx)), xx, clip_on=False, color=color, lw=0.5)
yax.set_xlim(left=0)
# threshold:
if 'cvbase' in xcol:
ax.text(xmax, 0.4*ymax, f'{100*np.sum(ydata > si_thresh)/len(data):.0f}\\%',
ha='right', va='bottom', fontsize='small')
ax.text(xmax, 0.3, f'{100*np.sum(ydata < si_thresh)/len(data):.0f}\\%',
ha='right', va='center', fontsize='small')
# statistics:
r, p = pearsonr(xdata, ydata)
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')
if 'cvbase' in xcol:
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:.2e}')
def si_stats(title, data, sicol, si_thresh, nsegscol):
print(title)
sidata = data[sicol]
@@ -176,7 +272,7 @@ def si_stats(title, data, sicol, si_thresh, nsegscol):
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}')
print(f' correlation {xcol:<11s} - {ycol:<11s}: r={r:5.2f}, p={p:.2e}')
def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
@@ -187,9 +283,11 @@ def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
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, 'cvbase', ycol, 'respmod2', 0, 250, 3,
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('Response mod.', 'Hz')
#cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 250, 3,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('Response mod.', 'Hz')
plot_corr_contrast(ax, data, 'cvbase', ycol, 3, color,
si_thresh, *examples)
def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
@@ -203,15 +301,17 @@ def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
#cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 250, 2,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('Response mod.', 'Hz')
cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 1.5, 2,
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('CV$_{\\rm base}$')
#cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 1.5, 2,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('CV$_{\\rm base}$')
#cax = plot_corr(ax, data, 'cvstim', ycol, 'ratebase', 50, 450, 2,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('$r$', 'Hz')
#cax = plot_corr(ax, data, 'cvstim', ycol, 'serialcorr1', -0.6, 0, 2,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('$\\rho_1$')
plot_corr_contrast(ax, data, 'cvstim', ycol, 2, color,
si_thresh, *examples)
def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
@@ -222,9 +322,11 @@ def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
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, 'respmod2', ycol, 'cvbase', 0, 1.5, 0.016,
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('CV$_{\\rm base}$')
#cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 1.5, 0.016,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('CV$_{\\rm base}$')
plot_corr_contrast(ax, data, 'respmod2', ycol, 0.016, color,
si_thresh, *examples)
def plot_rate_si_punit(ax, data, ycol, si_thresh, color):
@@ -235,9 +337,13 @@ def plot_rate_si_punit(ax, data, ycol, si_thresh, color):
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}$')
#cax = plot_corr(ax, data, 'ratebase', ycol, 'cvbase', 0, 1.5, 0.016,
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('CV$_{\\rm base}$')
plot_corr_contrast(ax, data, 'ratebase', ycol, 0.016, color,
si_thresh, *examples)
ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.9, 1),
markerfirst=False, handletextpad=0.5)
def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
@@ -247,9 +353,11 @@ def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_ylabel('SI($r$)')
ax.set_ylim(0, 10)
ax.set_yticks_delta(2)
cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 80, 20,
'coolwarm', color, si_thresh, *ampul_examples)
cax.set_ylabel('Response mod.', 'Hz')
#cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 80, 20,
# 'coolwarm', color, si_thresh, *ampul_examples)
#cax.set_ylabel('Response mod.', 'Hz')
plot_corr_contrast(ax, data, 'cvbase', ycol, 20, color,
si_thresh, *ampul_examples)
def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
@@ -262,14 +370,16 @@ def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
#cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 80, 6,
# 'coolwarm', color, si_thresh, *ampul_examples)
#cax.set_ylabel('Response mod.', 'Hz')
cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 0.2, 6,
'coolwarm', color, si_thresh, *ampul_examples)
cax.set_ylabel('CV$_{\\rm base}$')
cax.set_yticks_delta(0.1)
#cax = plot_corr(ax, data, 'cvstim', ycol, 'cvbase', 0, 0.2, 6,
# 'coolwarm', color, si_thresh, *ampul_examples)
#cax.set_ylabel('CV$_{\\rm base}$')
#cax.set_yticks_delta(0.1)
#cax = plot_corr(ax, data, 'cvstim', ycol, 'ratebase', 90, 180, 6,
# 'coolwarm', color, si_thresh, *ampul_examples)
#cax.set_ylabel('$r$', 'Hz')
#cax.set_yticks_delta(30)
plot_corr_contrast(ax, data, 'cvstim', ycol, 6, color,
si_thresh, *ampul_examples)
def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
@@ -279,10 +389,12 @@ def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_ylabel('SI($r$)')
ax.set_ylim(0, 10)
ax.set_yticks_delta(2)
cax = plot_corr(ax, data, 'respmod2', 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)
#cax = plot_corr(ax, data, 'respmod2', 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)
plot_corr_contrast(ax, data, 'respmod2', ycol, 0.06, color,
si_thresh, *ampul_examples)
def plot_rate_si_ampul(ax, data, ycol, si_thresh, color):
@@ -292,10 +404,14 @@ def plot_rate_si_ampul(ax, data, ycol, si_thresh, color):
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)
#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)
plot_corr_contrast(ax, data, 'ratebase', ycol, 0.06, color,
si_thresh, *ampul_examples)
ax.legend(title='contrast', loc='upper right', bbox_to_anchor=(1.95, 1),
markerfirst=False, handletextpad=0.5)
if __name__ == '__main__':
@@ -327,6 +443,14 @@ if __name__ == '__main__':
print(f' CV model: min={np.min(cvmodel):4.2f} max={np.max(cvmodel):4.2f} median={np.median(cvmodel):4.2f}')
print(f' CV data: min={np.min(cvdata):4.2f} max={np.max(cvdata):.2f} median={np.median(cvdata):4.2f}')
print()
ratemodel = punit_model['ratebase']
ratedata = punit_data['ratebase']
u, p = mannwhitneyu(ratemodel, ratedata)
print('Baseline rate differs between P-unit models and data:')
print(f' U={u:g}, p={p:.2g}')
print(f' baseline rate model: min={np.min(ratemodel):3.0f}Hz max={np.max(ratemodel):3.0f}Hz median={np.median(ratemodel):3.0f}Hz')
print(f' baseline rate data: min={np.min(ratedata):3.0f}Hz max={np.max(ratedata):3.0f}Hz median={np.median(ratedata):3.0f}Hz')
print()
rmmodel = punit_model['respmod2']
rmdata = punit_data['respmod2']
u, p = mannwhitneyu(rmmodel, rmdata)
@@ -343,12 +467,28 @@ if __name__ == '__main__':
print(f' SI model: min={np.min(simodel):4.1f} max={np.max(simodel):4.1f} median={np.median(simodel):4.1f}')
print(f' SI data: min={np.min(sidata):4.1f} max={np.max(sidata):4.1f} median={np.median(sidata):4.1f}')
print()
sipunit = punit_data['sinorm' + si]
siampul = ampul_data['sinorm' + si]
u, p = mannwhitneyu(sipunit, siampul)
print('SI differs between P-units and Ampullaries:')
print(f' U={u:g}, p={p:.2g}')
print(f' SI P-units: min={np.min(sipunit):4.1f} max={np.max(sipunit):4.1f} median={np.median(sipunit):4.1f}')
print(f' SI Ampullary: min={np.min(siampul):4.1f} max={np.max(siampul):4.1f} median={np.median(siampul):4.1f}')
print()
cvpunit = punit_data['cvbase']
cvampul = ampul_data['cvbase']
u, p = mannwhitneyu(cvpunit, cvampul)
print('CV differs between P-units and Ampullaries:')
print(f' U={u:g}, p={p:.2g}')
print(f' CV P-units: min={np.min(cvpunit):4.2f} max={np.max(cvpunit):4.2f} median={np.median(cvpunit):4.2f}')
print(f' CV Ampullary: min={np.min(cvampul):4.2f} max={np.max(cvampul):.2f} median={np.median(cvampul):4.2f}')
print()
s = plot_style()
fig, axs = plt.subplots(3, 3, cmsize=(s.plot_width, 0.75*s.plot_width),
height_ratios=[1, 0, 1, 0.3, 1])
fig.subplots_adjust(leftm=6.5, rightm=13.5, topm=4.5, bottomm=4,
wspace=1.1, hspace=0.6)
fig.subplots_adjust(leftm=6.5, rightm=18, topm=4.5, bottomm=4,
wspace=0.6, hspace=0.6)
si_stats('P-unit model:', punit_model, 'dsinorm100', si_thresh,
'nsegs100')