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

finished new figure captions

Browse Source
This commit is contained in:
Jan Benda
2025-05-23 11:25:53 +02:00
parent 3e07684093
commit 585e5ec27b
6 changed files with 116 additions and 156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

153
dataoverview.py
View File

@@ -9,16 +9,27 @@ from plotstyle import plot_style, lighter, significance_str
data_path = Path('data')
from noisesplit import model_cell as model_split_example
from modelsusceptcontrasts import model_cells as model_contrast_examples
from modelsusceptlown import model_cell as model_lown_example
from punitexamplecell import example_cell as punit_example
from punitexamplecell import example_cells as punit_examples
from noisesplit import example_cell as punit_split_example
from ampullaryexamplecell import example_cell as ampul_example
from ampullaryexamplecell import example_cells as ampul_examples
model_examples = ([[model_lown_example, 0.01],
[model_lown_example, 0.03],
[model_lown_example, 0.1]],
[[model_split_example, 0.01]],
[[m, a] for m in model_contrast_examples for a in [0.01, 0.03, 0.1]])
punit_examples = (punit_example, [punit_split_example], punit_examples)
ampul_examples = (ampul_example, [], ampul_examples)
def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
nli_thresh, example=[], examples=[]):
ax.axhline(nli_thresh, color='k', ls=':', lw=0.5)
si_thresh, example=[], split_example=[], examples=[]):
ax.axhline(si_thresh, color='k', ls=':', lw=0.5)
xmax = ax.get_xlim()[1]
ymax = ax.get_ylim()[1]
mask = (data[xcol] < xmax) & (data[ycol] < ymax)
@@ -31,13 +42,38 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=6, marker='^', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
for cell, run in split_example:
mask = (data['cell'] == cell) & (data['stimindex'] == run)
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=6, marker='s', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
for cell, run in examples:
mask = (data['cell'] == cell) & (data['stimindex'] == run)
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=5, marker='o', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
else:
for cell, alpha in example:
mask = (data['cell'] == cell) & (data['contrast'] == alpha)
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=6, marker='^', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
for cell, alpha in split_example:
mask = (data['cell'] == cell) & (data['contrast'] == alpha)
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=6, marker='s', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
for cell, alpha in examples:
mask = (data['cell'] == cell) & (data['contrast'] == alpha)
ax.scatter(data[mask, xcol], data[mask, ycol], c=data[mask, zcol],
s=5, marker='o', linewidth=0.5, edgecolors='black',
clip_on=False, cmap=cmap, vmin=zmin, vmax=zmax,
zorder=20)
# color bar:
fig = ax.get_figure()
cax = ax.inset_axes([1.3, 0, 0.04, 1])
@@ -66,9 +102,9 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
yax.set_xlim(left=0)
# threshold:
if 'cvbase' in xcol:
ax.text(xmax, 0.4*ymax, f'{100*np.sum(data[ycol] > nli_thresh)/len(data):.0f}\\%',
ax.text(xmax, 0.4*ymax, f'{100*np.sum(data[ycol] > si_thresh)/len(data):.0f}\\%',
ha='right', va='bottom', fontsize='small')
ax.text(xmax, 0.3, f'{100*np.sum(data[ycol] < nli_thresh)/len(data):.0f}\\%',
ax.text(xmax, 0.3, f'{100*np.sum(data[ycol] < si_thresh)/len(data):.0f}\\%',
ha='right', va='center', fontsize='small')
# statistics:
r, p = pearsonr(data[xcol], data[ycol])
@@ -83,76 +119,80 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color,
return cax
def nli_stats(title, data, column, nli_thresh):
def si_stats(title, data, column, si_thresh):
print(title)
print(f' nli threshold: {nli_thresh:.1f}')
cells = np.unique(data['cell'])
ncells = len(cells)
nrecs = len(data)
print(f' cells: {ncells}')
print(f' recordings: {nrecs}')
hcells = np.unique(data[data(column) > nli_thresh, 'cell'])
print(f' high nli cells: n={len(hcells):3d}, {100*len(hcells)/ncells:4.1f}%')
print(f' high nli recordings: n={np.sum(data(column) > nli_thresh):3d}, '
f'{100*np.sum(data(column) > nli_thresh)/nrecs:4.1f}%')
print(f' cells: {ncells}')
print(f' recordings: {nrecs}')
print(f' SI threshold: {si_thresh:.1f}')
hcells = np.unique(data[data(column) > si_thresh, 'cell'])
print(f' high SI cells: n={len(hcells):3d}, {100*len(hcells)/ncells:4.1f}%')
print(f' high SI recordings: n={np.sum(data(column) > si_thresh):3d}, '
f'{100*np.sum(data(column) > si_thresh)/nrecs:4.1f}%')
nsegs = data['nsegs']
print(f' number of segments: {np.min(nsegs):4.0f} - {np.max(nsegs):4.0f}, median={np.median(nsegs):4.0f}, mean={np.mean(nsegs):4.0f}, std={np.std(nsegs):4.0f}')
nrecs = []
for cell in cells:
nrecs.append(len(data[data["cell"] == cell, :]))
print(f' number of recordings per cell: {np.min(nrecs):4.0f} - {np.max(nrecs):4.0f}, median={np.median(nrecs):4.0f}, mean={np.mean(nrecs):4.0f}, std={np.std(nrecs):4.0f}')
fcutoff = data['fcutoff']
print(' cutoff frequencies:', ' '.join([f'{f:3.0f}Hz' for f in np.unique(fcutoff)]))
print(' cutoff frequencies:', ' '.join([f'{np.sum(fcutoff == f):3d}' for f in np.unique(fcutoff)]))
print(f' cutoff frequencies: {np.min(fcutoff):.0f}Hz - {np.max(fcutoff):.0f}Hz, median={np.median(fcutoff):.0f}Hz, mean={np.mean(fcutoff):.0f}Hz, std={np.std(fcutoff):.0f}Hz')
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}%')
def plot_cvbase_nli_punit(ax, data, ycol, nli_thresh, color):
def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color):
ax.set_xlabel('CV$_{\\rm base}$')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 1.5)
ax.set_ylim(0, 7.2)
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, nli_thresh,
punit_example, punit_examples)
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('Response mod.', 'Hz')
def plot_cvstim_nli_punit(ax, data, ycol, nli_thresh, color):
def plot_cvstim_si_punit(ax, data, ycol, si_thresh, color):
ax.set_xlabel('CV$_{\\rm stim}$')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 1.6)
ax.set_ylim(0, 7.2)
ax.set_xticks_delta(0.5)
ax.set_yticks_delta(2)
examples = punit_examples if 'stimindex' in data else model_examples
#cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 250, 2,
# 'coolwarm', color, nli_thresh,
# punit_example, punit_examples)
# '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, nli_thresh,
punit_example, punit_examples)
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('CV$_{\\rm base}$')
#cax = plot_corr(ax, data, 'cvstim', ycol, 'ratebase', 50, 450, 2,
# 'coolwarm', color, nli_thresh,
# punit_example, punit_examples)
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('$r$', 'Hz')
#cax = plot_corr(ax, data, 'cvstim', ycol, 'serialcorr1', -0.6, 0, 2,
# 'coolwarm', color, nli_thresh,
# punit_example, punit_examples)
# 'coolwarm', color, si_thresh, *examples)
#cax.set_ylabel('$\\rho_1$')
def plot_mod_nli_punit(ax, data, ycol, nli_thresh, color):
def plot_rmod_si_punit(ax, data, ycol, si_thresh, color):
ax.set_xlabel('Response modulation', 'Hz')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 250)
ax.set_ylim(0, 7.2)
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, nli_thresh,
punit_example, punit_examples)
'coolwarm', color, si_thresh, *examples)
cax.set_ylabel('CV$_{\\rm base}$')
def plot_cvbase_nli_ampul(ax, data, ycol, nli_thresh, color):
def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_xlabel('CV$_{\\rm base}$')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 0.2)
@@ -160,12 +200,11 @@ def plot_cvbase_nli_ampul(ax, data, ycol, nli_thresh, color):
ax.set_xticks_delta(0.1)
ax.set_yticks_delta(5)
cax = plot_corr(ax, data, 'cvbase', ycol, 'respmod2', 0, 80, 20,
'coolwarm', color, nli_thresh,
ampul_example, ampul_examples)
'coolwarm', color, si_thresh, *ampul_examples)
cax.set_ylabel('Response mod.', 'Hz')
def plot_cvstim_nli_ampul(ax, data, ycol, nli_thresh, color):
def plot_cvstim_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_xlabel('CV$_{\\rm stim}$')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 0.85)
@@ -173,22 +212,19 @@ def plot_cvstim_nli_ampul(ax, data, ycol, nli_thresh, color):
ax.set_xticks_delta(0.2)
ax.set_yticks_delta(5)
#cax = plot_corr(ax, data, 'cvstim', ycol, 'respmod2', 0, 80, 6,
# 'coolwarm', color, nli_thresh,
# ampul_example, ampul_examples)
# '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, nli_thresh,
ampul_example, ampul_examples)
'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, nli_thresh,
# ampul_example, ampul_examples)
# 'coolwarm', color, si_thresh, *ampul_examples)
#cax.set_ylabel('$r$', 'Hz')
#cax.set_yticks_delta(30)
def plot_mod_nli_ampul(ax, data, ycol, nli_thresh, color):
def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color):
ax.set_xlabel('Response modulation', 'Hz')
ax.set_ylabel('SI($r$)')
ax.set_xlim(0, 80)
@@ -196,8 +232,7 @@ def plot_mod_nli_ampul(ax, data, ycol, nli_thresh, color):
ax.set_xticks_delta(20)
ax.set_yticks_delta(5)
cax = plot_corr(ax, data, 'respmod2', ycol, 'cvbase', 0, 0.2, 0.06,
'coolwarm', color, nli_thresh,
ampul_example, ampul_examples)
'coolwarm', color, si_thresh, *ampul_examples)
cax.set_ylabel('CV$_{\\rm base}$')
cax.set_yticks_delta(0.1)
@@ -213,25 +248,25 @@ if __name__ == '__main__':
ampul_data = TableData(data_path /
'Apteronotus_leptorhynchus-Ampullary-data.csv',
sep=';')
nli_thresh = 1.8
si_thresh = 1.8
u, p = mannwhitneyu(punit_model['cvbase'], punit_data['cvbase'])
print('CV differs between P-unit models and data:')
print(f' U={u:g}, p={p:.2g}')
print(f' median model: {np.median(punit_model["cvbase"]):.2f}')
print(f' median data: {np.median(punit_data["cvbase"]):.2f}')
print(f' median data: {np.median(punit_data["cvbase"]):.2f}')
print()
u, p = mannwhitneyu(punit_model['respmod2'], punit_data['respmod2'])
print('Response modulation differs between P-unit models and data:')
print(f' U={u:g}, p={p:.2g}')
print(f' median model: {np.median(punit_model["respmod2"]):.2f}')
print(f' median data: {np.median(punit_data["respmod2"]):.2f}')
print(f' median data: {np.median(punit_data["respmod2"]):.2f}')
print()
u, p = mannwhitneyu(punit_model['dnli100'], punit_data['nli'])
print('NLI does not differ between P-unit models and data:')
print('SI does not differ between P-unit models and data:')
print(f' U={u:g}, p={p:.2g}')
print(f' median model: {np.median(punit_model["dnli100"]):.1f}')
print(f' median data: {np.median(punit_data["nli"]):.1f}')
print(f' median data: {np.median(punit_data["nli"]):.1f}')
print()
s = plot_style()
@@ -240,28 +275,28 @@ if __name__ == '__main__':
fig.subplots_adjust(leftm=6.5, rightm=13.5, topm=4.5, bottomm=4,
wspace=1.1, hspace=0.6)
nli_stats('P-unit model:', punit_model, 'dnli100', nli_thresh)
si_stats('P-unit model:', punit_model, 'dnli100', si_thresh)
axs[0, 0].text(0, 1.35, 'P-unit models',
transform=axs[0, 0].transAxes, color=s.model_color1)
plot_cvbase_nli_punit(axs[0, 0], punit_model, 'dnli100', nli_thresh, s.model_color2)
plot_mod_nli_punit(axs[0, 1], punit_model, 'dnli100', nli_thresh, s.model_color2)
plot_cvstim_nli_punit(axs[0, 2], punit_model, 'dnli100', nli_thresh, s.model_color2)
plot_cvbase_si_punit(axs[0, 0], punit_model, 'dnli100', si_thresh, s.model_color2)
plot_rmod_si_punit(axs[0, 1], punit_model, 'dnli100', si_thresh, s.model_color2)
plot_cvstim_si_punit(axs[0, 2], punit_model, 'dnli100', si_thresh, s.model_color2)
print()
nli_stats('P-unit data:', punit_data, 'nli', nli_thresh)
si_stats('P-unit data:', punit_data, 'nli', si_thresh)
axs[1, 0].text(0, 1.35, 'P-unit data',
transform=axs[1, 0].transAxes, color=s.punit_color1)
plot_cvbase_nli_punit(axs[1, 0], punit_data, 'nli', nli_thresh, s.punit_color2)
plot_mod_nli_punit(axs[1, 1], punit_data, 'nli', nli_thresh, s.punit_color2)
plot_cvstim_nli_punit(axs[1, 2], punit_data, 'nli', nli_thresh, s.punit_color2)
plot_cvbase_si_punit(axs[1, 0], punit_data, 'nli', si_thresh, s.punit_color2)
plot_rmod_si_punit(axs[1, 1], punit_data, 'nli', si_thresh, s.punit_color2)
plot_cvstim_si_punit(axs[1, 2], punit_data, 'nli', si_thresh, s.punit_color2)
print()
nli_stats('Ampullary data:', ampul_data, 'nli', nli_thresh)
si_stats('Ampullary data:', ampul_data, 'nli', si_thresh)
axs[2, 0].text(0, 1.35, 'Ampullary data',
transform=axs[2, 0].transAxes, color=s.ampul_color1)
plot_cvbase_nli_ampul(axs[2, 0], ampul_data, 'nli', nli_thresh, s.ampul_color2)
plot_mod_nli_ampul(axs[2, 1], ampul_data, 'nli', nli_thresh, s.ampul_color2)
plot_cvstim_nli_ampul(axs[2, 2], ampul_data, 'nli', nli_thresh, s.ampul_color2)
plot_cvbase_si_ampul(axs[2, 0], ampul_data, 'nli', si_thresh, s.ampul_color2)
plot_rmod_si_ampul(axs[2, 1], ampul_data, 'nli', si_thresh, s.ampul_color2)
plot_cvstim_si_ampul(axs[2, 2], ampul_data, 'nli', si_thresh, s.ampul_color2)
print()
fig.common_xticks(axs[:2, 0])
@@ -270,6 +305,6 @@ if __name__ == '__main__':
fig.common_yticks(axs[0, :])
fig.common_yticks(axs[1, :])
fig.common_yticks(axs[2, :])
fig.tag(xoffs=-3.5, yoffs=2)
fig.tag(axs, xoffs=-3.5, yoffs=2)
fig.savefig()
print()