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

updated model figures to new analysis with the right units

Browse Source
This commit is contained in:
Jan Benda
2025-05-25 12:29:46 +02:00
parent e785d51b18
commit e87d63c46b
5 changed files with 219 additions and 351 deletions
Download Patch File Download Diff File Expand all files Collapse all files

132
noisesplit.py
View File

@@ -1,8 +1,8 @@
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from spectral import whitenoise, diag_projection, peak_size
from plotstyle import plot_style
from spectral import whitenoise
from plotstyle import plot_style, noise_files, plot_chi2
example_cell = ['2017-07-18-ai-invivo-1', 1]
@@ -13,78 +13,21 @@ data_path = base_path / 'cells'
sims_path = base_path / 'simulations'
def sort_files(cell_name, all_files, n):
files = [fn for fn in all_files if '-'.join(fn.stem.split('-')[2:-n]) == cell_name]
if len(files) == 0:
return None, 0
nums = [int(fn.stem.split('-')[-1]) for fn in files]
idxs = np.argsort(nums)
files = [files[i] for i in idxs]
nums = [nums[i] for i in idxs]
return files, nums
def plot_chi2(ax, s, freqs, chi2, nsegs, rate):
fcutoff = 300
ax.set_aspect('equal')
i0 = np.argmin(freqs < 0)
i1 = np.argmax(freqs > fcutoff)
if i1 == 0:
i1 = len(freqs)
freqs = freqs[i0:i1]
chi2 = chi2[i0:i1, i0:i1]
vmax = np.quantile(chi2, 0.996)
ten = 10**np.floor(np.log10(vmax))
for fac, delta in zip([1, 2, 3, 4, 6, 8, 10],
[0.5, 1, 1, 2, 3, 4, 5]):
if fac*ten >= vmax:
vmax = prev_fac*ten
ten *= prev_delta
break
prev_fac = fac
prev_delta = delta
pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax,
rasterized=True)
ax.set_xlim(0, fcutoff)
ax.set_ylim(0, fcutoff)
ax.set_xticks_delta(100)
ax.set_yticks_delta(100)
ax.set_xlabel('$f_1$', 'Hz')
ax.set_ylabel('$f_2$', 'Hz')
if nsegs < 10000:
ax.text(1, 1.1, f'$N={nsegs}$',
ha='right', transform=ax.transAxes)
else:
ax.text(1, 1.1, f'$N=10^{{{np.log10(nsegs):.0f}}}$',
ha='right', transform=ax.transAxes)
dfreqs, diag = diag_projection(freqs, chi2, 2*fcutoff)
nli, nlirel, nlif = peak_size(dfreqs, diag, rate, median=False)
ax.text(0.95, 0.88, f'SI($r$)={nli:.1f}', ha='right', zorder=50,
color='white', fontsize='medium', transform=ax.transAxes)
cax = ax.inset_axes([1.04, 0, 0.05, 1])
cax.set_spines_outward('lrbt', 0)
cb = fig.colorbar(pc, cax=cax)
cb.outline.set_color('none')
cb.outline.set_linewidth(0)
cax.set_ylabel(r'$|\chi_2|$ [Hz]')
cax.set_yticks_delta(ten)
return cax
def plot_overn(ax, s, files, nmax=1e6):
fcutoff = 300
ns = []
stats = []
for fname in files:
data = np.load(fname)
n = data['n']
fcutoff = data['fcutoff']
n = data['nsegs']
if nmax is not None and n > nmax:
continue
alpha = data['alpha']
alpha = data['contrast']
freqs = data['freqs']
pss = data['pss']
prss = data['prss']
chi2 = np.abs(prss)/0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1))
chi2 *= 1e-4 # Hz/%^2
ns.append(n)
i0 = np.argmin(freqs < 0)
i1 = np.argmax(freqs > fcutoff)
@@ -114,26 +57,33 @@ def plot_overn(ax, s, files, nmax=1e6):
ax.set_xticks_log(numticks=4)
ax.set_minor_xticks_log(numticks=8)
else:
ax.set_ylim(4e0, 1.3e3)
ax.set_ylim(1e-1, 2e3)
#ax.set_minor_yticks_log(numticks=5)
ax.set_minor_yticks_off()
ax.set_xticks_log(numticks=5)
ax.set_xticks_log(numticks=6)
#ax.set_minor_xticks_log(numticks=6)
ax.set_xlabel('segments')
ax.set_ylabel('$|\\chi_2|$ [Hz]')
ax.set_ylabel(r'$|\chi_2|$', r'Hz/\%$^2$')
def plot_chi2_contrast(ax1, ax2, s, files, nums, nsmall, nlarge, rate):
for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
i = nums.index(n)
data = np.load(files[i])
n = data['n']
alpha = data['alpha']
nsegs = data['nsegs']
fcutoff = data['fcutoff']
alpha = data['contrast']
freqs = data['freqs']
pss = data['pss']
prss = data['prss']
if nsegs < 10000:
ax.text(1, 1.1, f'$N={nsegs}$',
ha='right', transform=ax.transAxes)
else:
ax.text(1, 1.1, f'$N=10^{{{np.log10(nsegs):.0f}}}$',
ha='right', transform=ax.transAxes)
chi2 = np.abs(prss)*0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1))
cax = plot_chi2(ax, s, freqs, chi2, n, rate)
cax = plot_chi2(ax, s, freqs, chi2, fcutoff, rate)
cax.set_ylabel('')
print(f'Modeled cell {"-".join(files[i].name.split("-")[2:-2])} at {100*alpha:4.1f}% contrast: noise_frac={1:3.1f}, nsegs={n}')
print()
@@ -143,14 +93,21 @@ def plot_chi2_split(ax1, ax2, s, files, nums, nsmall, nlarge, rate):
for ax, n in zip([ax1, ax2], [nsmall, nlarge]):
i = nums.index(n)
data = np.load(files[i])
n = data['n']
alpha = data['alpha']
nsegs = data['nsegs']
fcutoff = data['fcutoff']
alpha = data['contrast']
noise_frac = data['noise_frac']
freqs = data['freqs']
pss = data['pss']
prss = data['prss']
chi2 = np.abs(prss)*0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1))
cax = plot_chi2(ax, s, freqs, chi2, n, rate)
if nsegs < 10000:
ax.text(1, 1.1, f'$N={nsegs}$',
ha='right', transform=ax.transAxes)
else:
ax.text(1, 1.1, f'$N=10^{{{np.log10(nsegs):.0f}}}$',
ha='right', transform=ax.transAxes)
cax = plot_chi2(ax, s, freqs, chi2, fcutoff, rate)
cax.set_ylabel('')
print(f'Modeled cell {"-".join(files[i].name.split("-")[2:-1])} at {100*alpha:4.1f}% contrast: noise_frac={noise_frac:3.1f}, nsegs={n}')
print()
@@ -163,9 +120,10 @@ def plot_chi2_data(ax, s, cell_name, run):
eodf = float(data['eodf'])
ratebase = float(data['ratebase/Hz'])
cvbase = float(data['cvbase'])
data_file = data_path / f'{cell_name}-spectral-s{run:02d}.npz'
data_file = data_path / f'{cell_name}-spectral-all-s{run:02d}.npz'
data = np.load(data_file)
n = data['n']
nsegs = data['n']
fcutoff = data['fcutoff']
nfft = data['nfft']
deltat = data['deltat']
alpha = data['alpha']
@@ -173,9 +131,11 @@ def plot_chi2_data(ax, s, cell_name, run):
pss = data['pss']
prss = data['prss']
chi2 = np.abs(prss)*0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1))
print(f'Measured cell {"-".join(data_file.name.split("-")[:-2])} at {100*alpha:4.1f}% contrast: r={ratebase:3.0f}Hz, CV={cvbase:4.2f}, dt={1000*deltat:4.2f}ms, nfft={nfft}, win={1000*deltat*nfft:6.1f}ms, nsegs={n}')
print(f'Measured cell {"-".join(data_file.name.split("-")[:-2])} at {100*alpha:4.1f}% contrast: r={ratebase:3.0f}Hz, CV={cvbase:4.2f}, dt={1000*deltat:4.2f}ms, nfft={nfft}, win={1000*deltat*nfft:6.1f}ms, nsegs={nsegs}')
print()
plot_chi2(ax, s, freqs, chi2, n, ratebase)
ax.text(1, 1.1, f'$N={nsegs}$',
ha='right', transform=ax.transAxes)
plot_chi2(ax, s, freqs, chi2, fcutoff, ratebase)
return alpha, ratebase, eodf
@@ -277,34 +237,34 @@ if __name__ == '__main__':
# model 5%:
axss = axs[1]
data_files = sims_path.glob(f'chi2-noisen-{model_cell}-{1000*data_contrast:03.0f}-*.npz')
files, nums = sort_files(model_cell, data_files, 2)
files, nums = noise_files(sims_path, model_cell, data_contrast)
axss[1].text(xt, yt, 'P-unit model', fontsize='large',
transform=axs[1, 1].transAxes, color=s.model_color1)
plot_chi2_contrast(axss[1], axss[2], s, files, nums, nsmall, nlarge, ratebase)
plot_chi2_contrast(axss[1], axss[2], s, files, nums, nsmall, nlarge,
ratebase)
axr1 = plot_noise_split(axss[0], data_contrast, 0, 1, wtime, wnoise)
plot_overn(axss[3], s, files, nmax=1e6)
axss[3].legend(loc='lower center', bbox_to_anchor=(0.5, 1.2),
markerfirst=False, title='$|\\chi_2|$ percentiles')
# model 1%:
axss = axs[2]
data_files = sims_path.glob(f'chi2-noisen-{model_cell}-{1000*contrast:03.0f}-*.npz')
files, nums = sort_files(model_cell, data_files, 2)
plot_chi2_contrast(axss[1], axss[2], s, files, nums, nsmall, nlarge, ratebase)
files, nums = noise_files(sims_path, model_cell, contrast)
plot_chi2_contrast(axss[1], axss[2], s, files, nums, nsmall, nlarge,
ratebase)
axr2 = plot_noise_split(axss[0], contrast, 0, 1, wtime, wnoise)
plot_overn(axss[3], s, files, nmax=1e6)
# model noise split:
axss = axs[3]
data_files = sims_path.glob(f'chi2-split-{model_cell}-*.npz')
files, nums = sort_files(model_cell, data_files, 1)
files, nums = noise_files(sims_path, model_cell)
axss[1].text(xt, yt, 'P-unit model', fontsize='large',
transform=axss[1].transAxes, color=s.model_color1)
axss[1].text(xt + 0.9, yt, f'(noise split)', fontsize='large',
transform=axss[1].transAxes)
noise_contrast, noise_frac = plot_chi2_split(axss[1], axss[2], s,
files, nums, nsmall, nlarge, ratebase)
files, nums, nsmall, nlarge,
ratebase)
axr3 = plot_noise_split(axss[0], 0, noise_contrast, noise_frac,
wtime, wnoise)
plot_overn(axss[3], s, files, nmax=1e6)