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

updated data and figures to new analysis with the right units

Browse Source
This commit is contained in:
Jan Benda
2025-05-25 23:12:11 +02:00
parent e87d63c46b
commit d043a8c47c
423 changed files with 607 additions and 923 deletions
Download Patch File Download Diff File Expand all files Collapse all files

211
ampullaryexamplecell.py
View File

@@ -3,8 +3,10 @@ import matplotlib.pyplot as plt
from pathlib import Path
from spectral import diag_projection, peak_size
from plotstyle import plot_style
from plotstyle import plot_chi2
from punitexamplecell import load_baseline, load_noise, load_spectra
from punitexamplecell import plot_chi2, plot_colorbar
from punitexamplecell import plot_response_spectrum, plot_response
from punitexamplecell import plot_gain, plot_diagonals, plot_isih_small
example_cell = [['2012-05-15-ac', 3],
@@ -12,7 +14,7 @@ example_cell = [['2012-05-15-ac', 3],
example_cells = [
['2010-11-26-an', 0],
['2011-10-25-ac', 0],
['2010-11-08-aa', 1],
['2011-02-18-ab', 1],
['2014-01-16-aj', 5],
]
@@ -26,150 +28,38 @@ def plot_isih(ax, s, rate, cv, isis, pdf):
ax.set_xlim(0, 12)
ax.set_xticks_delta(4)
ax.set_xlabel('ISI', 'ms')
ax.text(0, 1.08, 'Ampullary:', transform=ax.transAxes, color=s.cell_color1,
fontsize='large')
ax.text(0.95, 1.08, f'$r={rate:.0f}$Hz, CV$_{{\\rm base}}$={cv:.2f}',
ax.text(0.95, 1.08, f'CV$_{{\\rm base}}$={cv:.2f}, $r={rate:.0f}$Hz',
transform=ax.transAxes)
def plot_isih2(ax, s, rate, cv, isis, pdf):
ax.show_spines('b')
ax.fill_between(1000*isis, pdf, facecolor=s.cell_color1)
ax.set_xlim(0, 20)
#ax.set_xticks_delta(5)
#ax.set_xticks_blank()
#ax.set_xticks_fixed([0, 5, 10, 15, 20], ['0', '', '', '', '20\\,ms'])
ax.set_xticks_fixed([0, 5, 10, 15, 20], ['0', '5', '10', '15', '20\\,ms'])
ax.text(1, 1.1, f'CV$_{{\\rm base}}$={cv:.2f}', ha='right',
transform=ax.transAxes)
ax.text(1, 0.6, f'$r={rate:.0f}$Hz', ha='right', transform=ax.transAxes)
def plot_response_spectrum(ax, s, eodf, rate, freqs, prr):
rate_i = np.argmax(prr[freqs < 0.7*eodf])
eod_i = np.argmax(prr[freqs > 500]) + np.argmax(freqs > 500)
power_db = 10*np.log10(prr/np.max(prr))
ax.show_spines('b')
mask = (freqs > 30) & (freqs < 890)
#ax.plot(freqs[mask], power_db[mask], **s.lsC1)
#ax.plot(freqs[eod_i], power_db[eod_i] + 2, **s.psFEOD)
#ax.plot(freqs[rate_i], power_db[rate_i] + 2, **s.psF0)
#ax.set_ylim(-25, 5)
ax.plot(freqs[mask], 1e-3*prr[mask], **s.lsC1)
ax.plot(freqs[eod_i], 1e-3*prr[eod_i] + 0.4, **s.psFEOD)
ax.plot(freqs[rate_i], 1e-3*prr[rate_i] + 0.4, **s.psF0)
ax.set_ylim(0, 6)
ax.set_xlim(0, 900)
ax.set_xticks_delta(300)
ax.set_xlabel('$f$', 'Hz')
#ax.text(freqs[eod_i], power_db[eod_i] + 4, '$f_{\\rm EOD}$',
# ha='center')
#ax.text(freqs[rate_i], power_db[rate_i] + 4, '$r$',
# ha='center')
#ax.yscalebar(1.05, 0, 10, 'dB', ha='right')
ax.text(freqs[eod_i], 1e-3*prr[eod_i] + 0.8, '$f_{\\rm EOD}$',
ha='center')
ax.text(freqs[rate_i], 1e-3*prr[rate_i] + 0.8, '$r$',
ha='center')
ax.yscalebar(1.05, 0, 1, 'kHz', ha='right')
def plot_response(ax, s, eodf, time1, stimulus1, contrast1, spikes1, contrast2, spikes2):
t0 = 0.3
t1 = 0.4
maxtrials = 8
trials = np.arange(maxtrials)
ax.show_spines('')
ax.eventplot(spikes1[2:2+maxtrials], lineoffsets=trials - maxtrials + 1,
linelength=0.8, linewidths=1, color=s.cell_color1)
ax.eventplot(spikes2[2:2+maxtrials], lineoffsets=trials - 2*maxtrials,
linelength=0.8, linewidths=1, color=s.cell_color2)
am = contrast1*stimulus1
eod = np.sin(2*np.pi*eodf*time1) + am
ax.plot(time1, 4*eod + 7, **s.lsEOD)
ax.plot(time1, 4*(1 + am) + 7, **s.lsAM)
ax.set_xlim(t0, t1)
ax.set_ylim(-2*maxtrials - 0.5, 14)
ax.xscalebar(1, -0.05, 0.01, None, '10\\,ms', ha='right')
ax.text(t1 + 0.003, -0.5*maxtrials, f'${100*contrast1:.0f}$\\,\\%',
va='center', color=s.cell_color1)
ax.text(t1 + 0.003, -1.55*maxtrials, f'${100*contrast2:.0f}$\\,\\%',
va='center', color=s.cell_color2)
def plot_gain(ax, s, fbase, contrast1, freqs1, gain1,
contrast2, freqs2, gain2, fcutoff):
ax.axvline(fbase, **s.lsGrid)
ax.plot(freqs2, 1e-2*gain2, label=f'{100*contrast2:.0f}', **s.lsC2)
ax.plot(freqs1, 1e-2*gain1, label=f'{100*contrast1:.0f}', **s.lsC1)
ax.set_xlim(0, fcutoff)
ax.set_ylim(0, 12)
ax.set_xticks_delta(50)
ax.set_yticks_delta(4)
ax.set_xlabel('$f$', 'Hz')
ax.set_ylabel(r'$|\chi_1|$', r'Hz/\%')
ax.text(fbase, 12.5, '$r$', ha='center')
def plot_diagonals(ax, s, fbase, contrast1, freqs1, chi21, contrast2, freqs2, chi22, fcutoff):
diags = []
sis = []
sips = []
sifs = []
for contrast, freqs, chi2 in [[contrast1, freqs1, chi21], [contrast2, freqs2, chi22]]:
dfreqs, diag = diag_projection(freqs, chi2, 2*fcutoff)
diags.append([dfreqs, diag])
sinorm, sirel, sif = peak_size(dfreqs, diag, fbase, median=False)
sip = diag[np.argmin(np.abs(dfreqs - sif))]
sis.append(sinorm)
sips.append(sip)
sifs.append(sif)
print(f' SI at {100*contrast:.1f}% contrast: {sinorm:.2f}')
#ax.axvline(fbase, **s.lsGrid)
ax.plot(diags[1][0], 1e-4*diags[1][1], **s.lsC2)
ax.plot(diags[0][0], 1e-4*diags[0][1], **s.lsC1)
ax.plot(sifs[1], 1e-4*sips[1] + 0.3, clip_on=False, **s.psC2)
ax.plot(sifs[0], 1e-4*sips[0] + 0.3, clip_on=False, **s.psC1)
ax.set_xlim(0, 2*fcutoff)
ax.set_ylim(0, 17)
ax.set_xticks_delta(100)
ax.set_yticks_delta(5)
ax.set_xlabel('$f_1 + f_2$', 'Hz')
#ax.set_ylabel(r'$|\chi_2|$', r'Hz/\%$^2$')
ax.text(sifs[1] - 25, 1e-4*sips[1], f'{100*contrast2:.0f}\\%',
ha='right')
ax.text(sifs[1] + 35, 1e-4*sips[1], f'SI={sis[1]:.1f}')
ax.text(sifs[0] - 25, 1e-4*sips[0] + 0.5, f'{100*contrast1:.0f}\\%',
ha='right')
ax.text(sifs[0] + 35, 1e-4*sips[0] + 0.5, f'SI={sis[0]:.1f}')
#ax.text(fbase, 1.75, '$r$', ha='center')
if __name__ == '__main__':
"""
# find a nice example cell:
from thunderlab.tabledata import TableData
data = TableData('data/Apteronotus_leptorhynchus-Ampullary-data.csv')
data = data[(data['fcutoff'] > 140) & (data['fcutoff'] < 160), :]
data = data[(data['nli'] > 2) & (data['nli'] < 2.5), :]
data = data[(data['respmod2'] > 20) & (data['respmod2'] < 100), :]
data = data[(data['cvbase'] > 0.05) & (data['cvbase'] < 0.2), :]
data = data[(data['ratebase'] > 100) & (data['ratebase'] < 180), :]
data = data[(data['sinorm_nmax'] > 5) & (data['sinorm_nmax'] < 50), :]
data = data[(data['contrast'] > 0.04) & (data['contrast'] < 0.06), :]
data = data[(data['respmod2'] > 0) & (data['respmod2'] < 100), :]
data = data[(data['cvbase'] > 0) & (data['cvbase'] < 0.2), :]
data = data[(data['ratebase'] > 50) & (data['ratebase'] < 180), :]
for k in range(data.rows()):
print(f'{data[k, "cell"]:<22s} s{data[k, "stimindex"]:02.0f}: '
f'{100*data[k, "contrast"]:3g}%, {data[k, "respmod2"]:3.0f}Hz, '
f'nli={data[k, "nli"]:5.2f}')
f'{100*data[k, "contrast"]:3g}%, r={data[k, "ratebase"]:3.0f}Hz, '
f'CV={data[k, "cvbase"]:4.2f}, '
f'rmod={data[k, "respmod2"]:3.0f}Hz, '
f'SI={data[k, "sinorm_nmax"]:5.2f}')
print()
exit()
#exit()
"""
#mode = 'all'
mode = '100'
cell_name = example_cell[0][0]
print('Example Ampullary cell:', cell_name)
print('Example Ampullary cell:')
eodf, rate, cv, isis, pdf, freqs, prr = load_baseline(data_path, cell_name)
print(f' baseline firing rate: {rate:.0f}Hz')
print(f' baseline firing CV : {cv:.2f}')
print(f' {cell_name:<22s}: fbase={rate:3.0f}Hz, CV={cv:.2f}')
contrast1, time1, stimulus1, spikes1 = load_noise(data_path,
*example_cell[0])
contrast2, time2, stimulus2, spikes2 = load_noise(data_path,
@@ -190,56 +80,59 @@ if __name__ == '__main__':
s.psC1 = s.psA1
s.psC2 = s.psA2
fig, (ax1, ax2, ax3) = \
plt.subplots(3, 1, height_ratios=[3, 0, 3, 0.2, 4.7],
plt.subplots(3, 1, height_ratios=[3, 0, 3, 0.3, 4.7],
cmsize=(s.plot_width, 0.85*s.plot_width))
fig.subplots_adjust(leftm=8, rightm=9, topm=2, bottomm=4,
wspace=0.4, hspace=0.4)
axi, axp, axr = ax1.subplots(1, 3, width_ratios=[2, 3, 0, 10])
axg, axc1, axc2, axd = ax2.subplots(1, 4, wspace=0.4)
axg = axg.subplots(1, 1, width_ratios=[1, 0.1])
axd = axd.subplots(1, 1, width_ratios=[0.2, 1])
axs = ax3.subplots(2, 4, wspace=0.4, hspace=0.35, height_ratios=[1, 4])
fig.subplots_adjust(leftm=8, rightm=2, topm=2, bottomm=3.5,
wspace=0.4, hspace=0.42)
axi, axp, axr = ax1.subplots(1, 3, width_ratios=[2, 3, 0, 10, 0.2])
axg, axc1, axc2, axd = ax2.subplots(1, 4, wspace=0.2,
width_ratios=[3.5, 0.5, 4, 4, 0.8, 3.5])
axs = ax3.subplots(2, 4, wspace=0.4, hspace=0.35,
width_ratios=[1, 0.1, 1, 1, 1, 0.1],
height_ratios=[1, 4])
axi.text(0, 1.08, 'Ampullary:', transform=axi.transAxes,
color=s.cell_color1, fontsize='large')
plot_isih(axi, s, rate, cv, isis, pdf)
plot_response_spectrum(axp, s, eodf, rate, freqs, prr)
plot_response(axr, s, eodf, time1, stimulus1, contrast1, spikes1,
contrast2, spikes2)
contrast2, spikes2, am=False)
plot_gain(axg, s, rate, contrast1, freqs1, gain1,
contrast2, freqs2, gain2, fcutoff1)
pc = plot_chi2(axc1, s, contrast2, freqs2, chi22, fcutoff2, 10)
axc1.plot([0, fcutoff2], [0, fcutoff2], zorder=20, **s.lsDiag)
plot_gain(axg, s, contrast1, freqs1, gain1,
contrast2, freqs2, gain2, fcutoff1, ymax=12, dy=4)
axg.axvline(rate, **s.lsGrid)
axg.text(rate, 12.5, '$r$', ha='center')
axc = plot_chi2(axc1, s, freqs2, chi22, fcutoff2, None, 12)
axc.remove()
axc1.set_title(f'$c$={100*contrast2:g}\\,\\%',
fontsize='medium', color=s.cell_color2)
pc = plot_chi2(axc2, s, contrast1, freqs1, chi21, fcutoff1, 10)
plot_chi2(axc2, s, freqs1, chi21, fcutoff1, None, 12)
axc2.set_title(f'$c$={100*contrast1:g}\\,\\%',
fontsize='medium', color=s.cell_color1)
axc2.plot([0, fcutoff1], [0, fcutoff1], zorder=20, **s.lsDiag)
plot_colorbar(axc2, pc, 2)
plot_diagonals(axd, s, rate, contrast1, freqs1, chi21,
contrast2, freqs2, chi22, fcutoff1)
contrast2, freqs2, chi22, fcutoff1, ymax=17, toffs=1)
fig.common_yticks(axc1, axc2)
fig.tag([axi, axp, axr], xoffs=-3, yoffs=-1)
fig.tag([axi, axp, axr], xoffs=-3, yoffs=0)
fig.tag([axg, axc1, axc2, axd], xoffs=-3, yoffs=2)
print()
print('Additional example cells:')
axs[0, 0].text(0, 1.6, 'Ampullary cells:', transform=axs[0, 0].transAxes,
color=s.cell_color1, fontsize='large')
for k, (cell, run) in enumerate(example_cells):
eodf, rate, cv, isis, pdf, _, _ = load_baseline(data_path, cell)
fcutoff, contrast, freqs, gain, chi2 = load_spectra(data_path, mode,
cell, run)
dfreqs, diag = diag_projection(freqs, chi2, 2*fcutoff)
sinorm, sirel, sif = peak_size(dfreqs, diag, rate, median=False)
print(f' {cell:<22s}: run={run:2d}, fbase={rate:3.0f}Hz, CV={cv:.2f}, SI={sinorm:3.1f}')
plot_isih2(axs[0, k], s, rate, cv, isis, pdf)
pc = plot_chi2(axs[1, k], s, contrast, freqs, chi2, fcutoff, 20)
axs[1, k].text(0.95, 0.9, f'SI($r$)={sinorm:.1f}', ha='right', zorder=50,
color='white', fontsize='medium',
transform=axs[1, k].transAxes)
axs[0, 0].text(0, 1.6, 'Ampullary cells:', transform=axs[0, 0].transAxes,
color=s.cell_color1,
fontsize='large')
plot_colorbar(axs[1, -1], pc, 5)
print(f' {cell:<22s}: run={run:2d}, contrast={100*contrast:3.2g}%, '
f'fbase={rate:3.0f}Hz, CV={cv:.2f}')
plot_isih_small(axs[0, k], s, contrast, rate, cv, isis, pdf)
vmax = 15 if k > 0 else 3
axc = plot_chi2(axs[1, k], s, freqs, chi2, fcutoff, rate, vmax)
if k % 3 != 0:
axc.remove()
if k == 0:
axc.set_ylabel('')
fig.common_yticks(axs[1, :])
fig.tag([axs[0, :]], xoffs=-3, yoffs=1)