diff --git a/modelsusceptovern.py b/modelsusceptovern.py index 40a2021..aeaae02 100644 --- a/modelsusceptovern.py +++ b/modelsusceptovern.py @@ -93,10 +93,10 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): ns = ns[indx] stats = stats[indx] ax.set_visible(True) - ax.plot(ns, stats[:, 7], '0.5', lw=1, zorder=50, label='99.8\\%') + ax.plot(ns, stats[:, 7], zorder=50, label='99.8\\%', **s.lsMax) ax.fill_between(ns, stats[:, 2], stats[:, 6], fc='0.85', zorder=40, label='5--95\\%') ax.fill_between(ns, stats[:, 3], stats[:, 5], fc='0.5', zorder=45, label='25-75\\%') - ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsSpine) + ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsMedian) #ax.plot(ns, stats[:, 8], '0.0') if title: if 'noise_frac' in data: diff --git a/noisesplit.py b/noisesplit.py index ae397b0..b3fca7e 100644 --- a/noisesplit.py +++ b/noisesplit.py @@ -1,6 +1,7 @@ import numpy as np import matplotlib.pyplot as plt from pathlib import Path +from spectral import whitenoise from plotstyle import plot_style @@ -33,9 +34,11 @@ def plot_chi2(ax, s, freqs, chi2, nsegs): 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 = fac*ten - ten *= delta + 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, 300) @@ -54,11 +57,60 @@ def plot_chi2(ax, s, freqs, chi2, nsegs): cb.outline.set_linewidth(0) cax.set_ylabel(r'$|\chi_2|$ [Hz]') cax.set_yticks_delta(ten) + return cax -def plot_chi2_contrast(ax1, ax2, s, cell_name, contrast, nsmall, nlarge): - data_files = sims_path.glob(f'chi2-noisen-{cell_name}-{1000*contrast:03.0f}-*.npz') - files, nums = sort_files(cell_name, data_files, 2) +def plot_overn(ax, s, files, nmax=1e6): + ns = [] + stats = [] + for fname in files: + data = np.load(fname) + n = data['n'] + if nmax is not None and n > nmax: + continue + alpha = data['alpha'] + freqs = data['freqs'] + pss = data['pss'] + dt_fix = 1 # 0.0005 + chi2 = np.abs(data['prss'])/dt_fix*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1)) + ns.append(n) + i0 = np.argmin(freqs < 0) + i1 = np.argmax(freqs > 300) + if i1 == 0: + i1 = len(freqs) + chi2 = chi2[i0:i1, i0:i1] + stats.append(np.quantile(chi2, [0, 0.001, 0.05, 0.25, 0.5, + 0.75, 0.95, 0.998, 1.0])) + ns = np.array(ns) + stats = np.array(stats) + indx = np.argsort(ns) + ns = ns[indx] + stats = stats[indx] + ax.set_visible(True) + ax.plot(ns, stats[:, 7], zorder=50, label='99.8\\%', **s.lsMax) + ax.fill_between(ns, stats[:, 2], stats[:, 6], fc='0.85', zorder=40, label='5--95\\%') + ax.fill_between(ns, stats[:, 3], stats[:, 5], fc='0.5', zorder=45, label='25-75\\%') + ax.plot(ns, stats[:, 4], zorder=50, label='median', **s.lsMedian) + #ax.plot(ns, stats[:, 8], '0.0') + ax.set_xlim(1e2, nmax) + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_yticks_log(numticks=5) + if nmax > 1e6: + ax.set_ylim(3e-1, 5e3) + ax.set_minor_yticks_log(numticks=5) + ax.set_xticks_log(numticks=4) + ax.set_minor_xticks_log(numticks=8) + else: + ax.set_ylim(1e0, 1.3e3) + #ax.set_minor_yticks_log(numticks=5) + ax.set_xticks_log(numticks=5) + #ax.set_minor_xticks_log(numticks=6) + ax.set_xlabel('segments') + ax.set_ylabel('$|\\chi_2|$ [Hz]') + + +def plot_chi2_contrast(ax1, ax2, s, files, nums, nsmall, nlarge): for ax, n in zip([ax1, ax2], [nsmall, nlarge]): i = nums.index(n) data = np.load(files[i]) @@ -67,12 +119,11 @@ def plot_chi2_contrast(ax1, ax2, s, cell_name, contrast, nsmall, nlarge): freqs = data['freqs'] pss = data['pss'] chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1)) - plot_chi2(ax, s, freqs, chi2, n) + cax = plot_chi2(ax, s, freqs, chi2, n) + cax.set_ylabel('') -def plot_chi2_split(ax1, ax2, s, cell_name, nsmall, nlarge): - data_files = sims_path.glob(f'chi2-split-{cell_name}-*.npz') - files, nums = sort_files(cell_name, data_files, 1) +def plot_chi2_split(ax1, ax2, s, files, nums, nsmall, nlarge): for ax, n in zip([ax1, ax2], [nsmall, nlarge]): i = nums.index(n) data = np.load(files[i]) @@ -82,11 +133,17 @@ def plot_chi2_split(ax1, ax2, s, cell_name, nsmall, nlarge): freqs = data['freqs'] pss = data['pss'] chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1)) - plot_chi2(ax, s, freqs, chi2, n) + cax = plot_chi2(ax, s, freqs, chi2, n) + cax.set_ylabel('') return alpha, noise_frac def plot_chi2_data(ax, s, cell_name, run): + data_file = data_path / f'{cell_name}-baseline.npz' + data = np.load(data_file) + 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 = np.load(data_file) n = data['n'] @@ -94,57 +151,138 @@ def plot_chi2_data(ax, s, cell_name, run): freqs = data['freqs'] pss = data['pss'] chi2 = np.abs(data['prss'])*0.5/np.sqrt(pss.reshape(1, -1)*pss.reshape(-1, 1)) - print(f'Measured cell {data_file.name} at {100*alpha:.1f}% contrast') + print(f'Measured cell {data_file.name} at {100*alpha:.1f}% contrast: r={ratebase:3.0f}Hz, CV={cvbase:4.2f}') plot_chi2(ax, s, freqs, chi2, n) - return alpha + return alpha, ratebase, eodf -def plot_noise_split(ax, contrast, noise_contrast, noise_frac): - axr, axs, axn = ax.subplots(3, 1, hspace=0.1) +def plot_ram(ax, contrast, eodf, wtime, wnoise): tmax = 50 + am = 1 + contrast*wnoise + eod = np.sin(2*np.pi*eodf*wtime)*am + + ax.show_spines('l') + ax.plot(1e3*wtime, eod, clip_on=False, **s.lsEOD) + ax.plot(1e3*wtime, +am, clip_on=False, **s.lsAM) + ax.plot(1e3*wtime, -am, clip_on=False, **s.lsAM) + ax.set_xlim(0, tmax) + ax.set_ylim(-1.3, 1.3) + ax.set_yticks_delta(1) + ax.set_ylabel('EOD') + ax.text(1, 1, f'RAM ($c={100*contrast:.0f}$\\,\\%)', ha='right', + transform=ax.transAxes, color=s.lsAM['color']) + +def plot_noise_split(ax, contrast, noise_contrast, noise_frac, + wtime, wnoise): + axr, axs, axn = ax.subplots(3, 1, hspace=0.2) + cmax = 26 + cdelta = 20 + tmax = 50 + axr.show_spines('l') + axr.axhline(0, **s.lsGrid) + axr.plot(1e3*wtime, 100*contrast*wnoise, clip_on=False, **s.lsAM) axr.set_xlim(0, tmax) - axr.set_ylim(-8, 8) - axr.set_yticks_delta(6) + axr.set_ylim(-cmax, cmax) + axr.set_yticks_delta(cdelta) axr.set_ylabel('\\%') + axr.text(1, 1, f'RAM ($c={100*contrast:.0f}$\\,\\%)', ha='right', + transform=axr.transAxes, color=s.lsAM['color']) axs.show_spines('l') + axs.axhline(0, **s.lsGrid) + axs.plot(1e3*wtime, 100*noise_contrast*wnoise, clip_on=False, **s.lsAMsplit) axs.set_xlim(0, tmax) - axs.set_ylim(-8, 8) - axs.set_yticks_delta(6) + axs.set_ylim(-cmax, cmax) + axs.set_yticks_delta(cdelta) axs.set_ylabel('\\%') - - axn.set_ylim(-6, 6) + if noise_contrast > 0: + axs.text(1, 1, f'$s_{{\\xi}}(t)$ ($c={100*noise_contrast:.0f}$\\,\\%)', + ha='right', transform=axs.transAxes, + color=s.lsAMsplit['color']) + + ntime = np.linspace(0, 1e-3*tmax, 800) + rng = np.random.default_rng(45432) + nnoise = rng.normal(size=len(ntime)) + axn.show_spines('l') + axn.axhline(0, **s.lsGrid) + axn.plot(1e3*ntime, noise_frac*nnoise, clip_on=False, **s.lsNoise) + axn.set_ylim(-2, 2) axn.set_xlim(0, tmax) - axn.set_yticks_delta(6) + axn.set_yticks_delta(5) axn.set_yticks_blank() - axn.set_xticks_delta(25) - axn.set_xlabel('Time', 'ms') + #axn.set_xticks_delta(25) + #axn.set_xlabel('Time', 'ms') + y = 0.8 if noise_frac < 1 else 1.2 + axn.text(1, y, f'Intrinsic noise (${100*noise_frac:.0f}$\\,\\%)', + ha='right', transform=axn.transAxes, color=s.lsNoise['color']) + if noise_frac < 1: + axn.xscalebar(1, 0, 10, 'ms', ha='right') + + return axr if __name__ == '__main__': - cell_name = '2012-07-03-ak-invivo-1' + #cell_name = ['2012-07-03-ak-invivo-1', 0] + cell_name = ['2017-07-18-ai-invivo-1', 1] # Take this! at 3% model, 5% data nsmall = 100 nlarge = 1000000 contrast = 0.03 + + wdt = 0.0001 + wnoise = whitenoise(0, 300, wdt, 0.05, rng=np.random.default_rng(51234)) + wtime = np.arange(len(wnoise))*wdt s = plot_style() - fig, axs = plt.subplots(2, 4, cmsize=(s.plot_width, 0.4*s.plot_width), - width_ratios=[1, 0, 1, 1, 1]) - fig.subplots_adjust(leftm=7, rightm=8, topm=2, bottomm=3.5, - wspace=0.4, hspace=0.6) - axs[1, 0].set_visible(False) - data_contrast = plot_chi2_data(axs[0, 0], s, cell_name[:13], 0) - plot_noise_split(axs[0, 1], data_contrast, 0, 1) - plot_chi2_contrast(axs[0, 2], axs[0, 3], s, cell_name, contrast, nsmall, nlarge) - noise_contrast, noise_frac = plot_chi2_split(axs[1, 2], axs[1, 3], s, - cell_name, nsmall, nlarge) - plot_noise_split(axs[1, 1], contrast, noise_contrast, noise_frac) + fig, axs = plt.subplots(3, 4, cmsize=(s.plot_width, 0.7*s.plot_width), + width_ratios=[1, 0, 1, 1, 0.15, 1]) + fig.subplots_adjust(leftm=8, rightm=1.5, topm=3, bottomm=4, + wspace=0.25, hspace=0.8) + axs[0, 2].set_visible(False) + axs[0, 3].set_visible(False) + + # data: + axs[0, 1].text(-2.42, 1.2, 'P-unit data', fontsize='large', + transform=axs[0, 1].transAxes, color=s.punit_color1) + data_contrast, ratebase, eodf = plot_chi2_data(axs[0, 1], s, cell_name[0], + cell_name[1]) + plot_ram(axs[0, 0], data_contrast, eodf, wtime, wnoise) + axs[0, 1].text(-1.5, 1.2, f'$r={ratebase:.0f}$\\,Hz', + transform=axs[0, 1].transAxes, fontsize='large') + + # model: + data_files = sims_path.glob(f'chi2-noisen-{cell_name[0]}-{1000*contrast:03.0f}-*.npz') + files, nums = sort_files(cell_name[0], data_files, 2) + axs[1, 1].text(-2.42, 1.2, 'P-unit model', fontsize='large', + transform=axs[1, 1].transAxes, color=s.model_color1) + plot_chi2_contrast(axs[1, 1], axs[1, 2], s, files, nums, nsmall, nlarge) + axr1 = plot_noise_split(axs[1, 0], contrast, 0, 1, wtime, wnoise) + plot_overn(axs[1, 3], s, files, nmax=1e6) + axs[1, 3].legend(loc='lower center', bbox_to_anchor=(0.5, 1.1), + markerfirst=False, title='$|\\chi_2|$ percentiles') + + # model noise split: + data_files = sims_path.glob(f'chi2-split-{cell_name[0]}-*.npz') + files, nums = sort_files(cell_name[0], data_files, 1) + axs[2, 1].text(-2.42, 1.2, 'P-unit model', fontsize='large', + transform=axs[2, 1].transAxes, color=s.model_color1) + axs[2, 1].text(-1.5, 1.2, f'(noise split)', fontsize='large', + transform=axs[2, 1].transAxes) + noise_contrast, noise_frac = plot_chi2_split(axs[2, 1], axs[2, 2], s, + files, nums, nsmall, nlarge) + axr2 = plot_noise_split(axs[2, 0], 0, noise_contrast, noise_frac, + wtime, wnoise) + plot_overn(axs[2, 3], s, files, nmax=1e6) + + fig.common_xticks(axs[:, 1]) fig.common_xticks(axs[:, 2]) fig.common_xticks(axs[:, 3]) - fig.common_yticks(axs[0, 2:]) - fig.common_yticks(axs[1, 2:]) - #fig.tag(axs, xoffs=-4.5, yoffs=1.8) + fig.common_yticks(axs[1, 1:3]) + fig.common_yticks(axs[2, 1:3]) + fig.tag([axs[0, :2], + [axr1] + axs[1, 1:].tolist(), + [axr2] + axs[2, 1:].tolist()], + xoffs=-4.5, yoffs=2) fig.savefig() print() diff --git a/plotstyle.py b/plotstyle.py index 3fb2a95..b67f0ad 100644 --- a/plotstyle.py +++ b/plotstyle.py @@ -55,7 +55,11 @@ def plot_style(): pt.make_line_styles(ns, 'ls', 'Line', '', palette['black'], '-', lwthin) pt.make_line_styles(ns, 'ls', 'EOD', '', palette['gray'], '-', lwthin) - pt.make_line_styles(ns, 'ls', 'AM', '', palette['red'], '-', lwthick) + pt.make_line_styles(ns, 'ls', 'AM', '', palette['red'], '-', lwmid) + pt.make_line_styles(ns, 'ls', 'AMsplit', '', palette['orange'], '-', lwmid) + pt.make_line_styles(ns, 'ls', 'Noise', '', palette['gray'], '-', lwmid) + pt.make_line_styles(ns, 'ls', 'Median', '', palette['black'], '-', lwthick) + pt.make_line_styles(ns, 'ls', 'Max', '', palette['black'], '-', lwmid) ns.lsStim = dict(color='gray', lw=ns.lwmid) ns.lsRaster = dict(color='black', lw=ns.lwthin) diff --git a/punitexamplecell.py b/punitexamplecell.py index 3a1b6af..2b037c3 100644 --- a/punitexamplecell.py +++ b/punitexamplecell.py @@ -14,9 +14,9 @@ run2 = 1 example_cells = [ ['2021-06-18-ae-invivo-1', 3], # 98Hz, 1%, ok ['2012-03-30-ah', 2], # 177Hz, 2.5%, 2.0, nice - ##['2012-07-03-ak', 0], # 120Hz, 2.5%, 1.8, broader - ##['2012-12-20-ac', 0], # 213Hz, 2.5%, 2.1, ok - #['2017-07-18-ai-invivo-1', 1], # 78Hz, 5%, 2.3, weak + ##['2012-07-03-ak', 0], # 120Hz, 2.5%, 1.8, broader, the one model cell, nice triangle up to 1%! + ##['2012-12-20-ac', 0], # 213Hz, 2.5%, 2.1, ok, model cell, weak triangle up to 1%! + #['2017-07-18-ai-invivo-1', 1], # 78Hz, 5%, 2.3, weak, nice model cell with clear triangle up to 10%! ##['2019-06-28-ae', 0], # 477Hz, 10%, 2.6, weak ##['2020-10-27-aa-invivo-1', 4], # 259Hz, 0.5%, 2.0, ok ##['2020-10-27-ae-invivo-1', 4], # 375Hz, 0.5%, 4.3, nice, additional low freq line