diff --git a/modelsusceptcontrasts.py b/modelsusceptcontrasts.py index d5ced2c..0ed5c51 100644 --- a/modelsusceptcontrasts.py +++ b/modelsusceptcontrasts.py @@ -3,8 +3,8 @@ import matplotlib.pyplot as plt from scipy.stats import pearsonr, linregress, gaussian_kde from thunderlab.tabledata import TableData from pathlib import Path -from spectral import whitenoise, diag_projection, peak_size -from plotstyle import plot_style, labels_params, significance_str +from plotstyle import plot_style, labels_params +from plotstyle import noise_files, plot_chi2, significance_str model_cells = ['2017-07-18-ai-invivo-1', # strong triangle @@ -16,101 +16,55 @@ data_path = Path('data') sims_path = data_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, data_file, rate): - fcutoff = 300 - data = np.load(data_file) - n = data['n'] - alpha = data['alpha'] +def load_chi2(file_path, cell_name, contrast=None, n=None): + files, nums = noise_files(sims_path, cell_name, contrast) + idx = -1 if n is None else nums.index(n) + data = np.load(files[idx]) + n = data['nsegs'] + fcutoff = data['fcutoff'] + contrast = 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)) - ax.set_visible(True) - ax.set_aspect('equal') - i0 = np.argmin(freqs < -fcutoff) - 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 = fac*ten - ten *= delta - break - pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax, - rasterized=True) - if 'noise_frac' in data: - ax.set_title('$c$=0\\,\\%', fontsize='medium') - else: - ax.set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium') - 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') - dfreqs, diag = diag_projection(freqs, chi2, 2*fcutoff) - sinorm, sirel, sif = peak_size(dfreqs, diag, rate, median=False) - ax.text(0.95, 0.88, f'SI($r$)={sinorm:.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) - if alpha == 0.1: - cb = fig.colorbar(pc, cax=cax, label=r'$|\chi_2|$ [Hz]') - else: - cb = fig.colorbar(pc, cax=cax) - cb.outline.set_color('none') - cb.outline.set_linewidth(0) - cax.set_yticks_delta(ten) + return freqs, chi2, fcutoff, contrast, n + - def plot_chi2_contrasts(axs, s, cell_name): - d = sims_path / f'baseline-{cell_name}.npz' + d = sims_path / f'{cell_name}-baseline.npz' data = np.load(d) rate = float(data['rate']) cv = float(data['cv']) print(f' {cell_name}: r={rate:3.0f}Hz, CV={cv:4.2f}') - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-split-{cell_name}-*.npz'), 1) - plot_chi2(axs[0], s, files[-1], rate) - for k, alphastr in enumerate(['010', '030', '100']): - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-noisen-{cell_name}-{alphastr}-*.npz'), 2) - plot_chi2(axs[k + 1], s, files[-1], rate) + freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name) + cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate) + cax.set_ylabel('') + axs[0].set_title(r'$c$=0\,\%', fontsize='medium') + for k, alpha in enumerate([0.01, 0.03, 0.1]): + freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name, + alpha) + cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate) + if alpha < 0.1: + cax.set_ylabel('') + axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium') def plot_si_cv(ax, s, data, alpha, cells): data = data[data['contrast'] == alpha, :] - r, p = pearsonr(data['cvbase'], data['dnli']) - l = linregress(data['cvbase'], data['dnli']) + r, p = pearsonr(data['cvbase'], data['dsinorm']) + l = linregress(data['cvbase'], data['dsinorm']) x = np.linspace(0, 1, 10) ax.set_visible(True) ax.set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium') ax.axhline(1, **s.lsLine) ax.plot(x, l.slope*x + l.intercept, **s.lsGrid) mask = data['triangle'] > 0.5 - ax.plot(data[mask, 'cvbase'], data[mask, 'dnli'], + ax.plot(data[mask, 'cvbase'], data[mask, 'dsinorm'], clip_on=False, zorder=30, label='strong', **s.psA1m) mask = data['border'] > 0.5 - ax.plot(data[mask, 'cvbase'], data[mask, 'dnli'], + ax.plot(data[mask, 'cvbase'], data[mask, 'dsinorm'], zorder=20, label='weak', **s.psA2m) - ax.plot(data['cvbase'], data['dnli'], clip_on=False, + ax.plot(data['cvbase'], data['dsinorm'], clip_on=False, zorder=10, label='none', **s.psB1m) for cell_name in cells: @@ -120,16 +74,16 @@ def plot_si_cv(ax, s, data, alpha, cells): color = s.psA2m['color'] elif data[mask, 'triangle']: color = s.psA1m['color'] - ax.plot(data[mask, 'cvbase'], data[mask, 'dnli'], + ax.plot(data[mask, 'cvbase'], data[mask, 'dsinorm'], zorder=40, marker='o', ms=s.psB1m['markersize'], mfc=color, mec='k', mew=0.8) - ax.set_ylim(0, 8) ax.set_xlim(0, 1) + ax.set_ylim(0, 9) + ax.set_yticks_delta(3) ax.set_minor_yticks_delta(1) ax.set_xlabel('CV$_{\\rm base}$') - ax.set_ylabel('SI') - ax.set_yticks_delta(4) + ax.set_ylabel('SI($r$)') ax.text(1, 0.9, f'$R={r:.2f}$', transform=ax.transAxes, ha='right', fontsize='small') ax.text(1, 0.75, significance_str(p), transform=ax.transAxes, @@ -139,13 +93,13 @@ def plot_si_cv(ax, s, data, alpha, cells): title='triangle', handlelength=0.5, handletextpad=0.5, labelspacing=0.2) - kde = gaussian_kde(data['dnli'], 0.15/np.std(data['dnli'], ddof=1)) - nli = np.linspace(0, 8, 100) - pdf = kde(nli) + kde = gaussian_kde(data['dsinorm'], 0.15/np.std(data['dsinorm'], ddof=1)) + si = np.linspace(0, 8, 100) + pdf = kde(si) dax = ax.inset_axes([1.04, 0, 0.3, 1]) dax.show_spines('') - dax.fill_betweenx(nli, pdf, **s.fsB1a) - dax.plot(pdf, nli, clip_on=False, **s.lsB1m) + dax.fill_betweenx(si, pdf, **s.fsB1a) + dax.plot(pdf, si, clip_on=False, **s.lsB1m) def plot_summary_contrasts(axs, s, cells): @@ -154,7 +108,7 @@ def plot_summary_contrasts(axs, s, cells): plot_si_cv(axs[0], s, data, 0, cells) print('noise split:') cdata = data[data['contrast'] == 0, :] - nli_split = cdata['dnli'] + nli_split = cdata['dsinorm'] print(f' mean SI = {np.mean(nli_split):.2f}, stdev = {np.std(nli_split):.2f}') n = np.sum(nli_split > nli_thresh) print(f' {n} cells ({100*n/len(nli_split):.1f}%) have SI > {nli_thresh:.1f}:') @@ -163,10 +117,10 @@ def plot_summary_contrasts(axs, s, cells): print(f' triangle cells have SI >= {np.min(nli_split[cdata["triangle"] > 0.5]):.2f}') print() for i, a in enumerate([0.01, 0.03, 0.1]): - plot_nli_cv(axs[1 + i], s, data, a, cells) + plot_si_cv(axs[1 + i], s, data, a, cells) print(f'contrast {100*a:2g}%:') cdata = data[data['contrast'] == a, :] - nli = cdata['dnli'] + nli = cdata['dsinorm'] r, p = pearsonr(nli_split, nli) print(f' correlation with split: r={r:.2f}, p={p:.1e}') print(f' mean SI = {np.mean(nli):.2f}, stdev = {np.std(nli):.2f}') diff --git a/modelsusceptlown.py b/modelsusceptlown.py index b463123..5861852 100644 --- a/modelsusceptlown.py +++ b/modelsusceptlown.py @@ -3,8 +3,9 @@ import matplotlib.pyplot as plt from scipy.stats import pearsonr, linregress, gaussian_kde from thunderlab.tabledata import TableData from pathlib import Path -from spectral import whitenoise, diag_projection, peak_size -from plotstyle import plot_style, labels_params, significance_str +from plotstyle import plot_style, labels_params +from plotstyle import noise_files, plot_chi2, significance_str +from modelsusceptcontrasts import load_chi2 model_cell = '2012-12-21-ak-invivo-1' @@ -13,103 +14,44 @@ data_path = Path('data') sims_path = data_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, data_file, rate): - fcutoff = 300 - data = np.load(data_file) - n = data['n'] - alpha = data['alpha'] - freqs = data['freqs'] - pss = data['pss'] - prss = data['prss'] - chi2 = np.abs(prss)*0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1)) - ax.set_visible(True) - ax.set_aspect('equal') - i0 = np.argmin(freqs < -fcutoff) - 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 = fac*ten - ten *= delta - break - pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax, - rasterized=True) - ns = f'$N={n}$' if n <= 100 else f'$N=10^{np.log10(n):.0f}$' - if 'noise_frac' in data: - ax.set_title(f'$c$=0\\,\\%, {ns}', fontsize='medium') - else: - ax.set_title(f'$c$={100*alpha:g}\\,\\%, {ns}', fontsize='medium') - 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') - dfreqs, diag = diag_projection(freqs, chi2, 2*fcutoff) - sinorm, sirel, sif = peak_size(dfreqs, diag, rate, median=False) - ax.text(0.95, 0.88, f'SI($r$)={sinorm:.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) - if alpha == 0.1: - cb = fig.colorbar(pc, cax=cax, label=r'$|\chi_2|$ [Hz]') - else: - cb = fig.colorbar(pc, cax=cax) - cb.outline.set_color('none') - cb.outline.set_linewidth(0) - cax.set_yticks_delta(ten) - - -def plot_chi2_contrasts(axs, s, cell_name, n=None): - d = sims_path / f'baseline-{cell_name}.npz' +def plot_chi2_contrasts(axs, s, cell_name, nsegs=None): + d = sims_path / f'{cell_name}-baseline.npz' data = np.load(d) rate = float(data['rate']) cv = float(data['cv']) print(f' {cell_name}: r={rate:3.0f}Hz, CV={cv:4.2f}') - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-split-{cell_name}-*.npz'), 1) - idx = -1 if n is None else nums.index(n) - plot_chi2(axs[0], s, files[idx], rate) - for k, alphastr in enumerate(['010', '030', '100']): - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-noisen-{cell_name}-{alphastr}-*.npz'), 2) - idx = -1 if n is None else nums.index(n) - plot_chi2(axs[k + 1], s, files[idx], rate) - + freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name, + None, nsegs) + ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$' + cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate) + cax.set_ylabel('') + axs[0].set_title(f'$c$=0\\,\\%, {ns}', fontsize='medium') + for k, alpha in enumerate([0.01, 0.03, 0.1]): + freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name, + alpha, nsegs) + ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$' + cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate) + if alpha < 0.1: + cax.set_ylabel('') + axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%, {ns}', + fontsize='medium') -def plot_nli_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): + +def plot_si_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): datax = data[data['contrast'] == alphax, :] datay = data[data['contrast'] == alphay, :] - nlix = datax['dnli'] - nliy = datay['dnli100'] - nfp = np.sum((nliy > ythresh) & (nlix < xthresh)) - ntp = np.sum((nliy > ythresh) & (nlix > xthresh)) - ntn = np.sum((nliy < ythresh) & (nlix < xthresh)) - nfn = np.sum((nliy < ythresh) & (nlix > xthresh)) - print(f' {ntp:2d} ({100*ntp/len(nlix):2.0f}%) true positive') - print(f' {nfp:2d} ({100*nfp/len(nlix):2.0f}%) false positive') - print(f' {ntn:2d} ({100*ntn/len(nlix):2.0f}%) true negative') - print(f' {nfn:2d} ({100*nfn/len(nlix):2.0f}%) false negative') - r, p = pearsonr(nlix, nliy) - l = linregress(nlix, nliy) + six = datax['dsinorm'] + siy = datay['dsinorm100'] + nfp = np.sum((siy > ythresh) & (six < xthresh)) + ntp = np.sum((siy > ythresh) & (six > xthresh)) + ntn = np.sum((siy < ythresh) & (six < xthresh)) + nfn = np.sum((siy < ythresh) & (six > xthresh)) + print(f' {ntp:2d} ({100*ntp/len(six):2.0f}%) true positive') + print(f' {nfp:2d} ({100*nfp/len(six):2.0f}%) false positive') + print(f' {ntn:2d} ({100*ntn/len(six):2.0f}%) true negative') + print(f' {nfn:2d} ({100*nfn/len(six):2.0f}%) false negative') + r, p = pearsonr(six, siy) + l = linregress(six, siy) x = np.linspace(0, 10, 10) ax.set_visible(True) ax.set_title(f'$c$={100*alphay:g}\\,\\%', fontsize='medium') @@ -119,10 +61,10 @@ def plot_nli_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): ax.axvline(xthresh, 0, 0.5, **s.lsLine) if alphax == 0: mask = datax['triangle'] > 0.5 - ax.plot(nlix[mask], nliy[mask], zorder=30, label='strong', **s.psA1m) + ax.plot(six[mask], siy[mask], zorder=30, label='strong', **s.psA1m) mask = datax['border'] > 0.5 - ax.plot(nliy[mask], nliy[mask], zorder=20, label='weak', **s.psA2m) - ax.plot(nlix, nliy, zorder=10, label='none', **s.psB1m) + ax.plot(siy[mask], siy[mask], zorder=20, label='weak', **s.psA2m) + ax.plot(six, siy, zorder=10, label='none', **s.psB1m) # mark cell: mask = datax['cell'] == cell_name color = s.psB1m['color'] @@ -131,7 +73,7 @@ def plot_nli_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): color = s.psA2m['color'] elif datax[mask, 'triangle']: color = s.psA1m['color'] - ax.plot(nlix[mask], nliy[mask], zorder=40, marker='o', + ax.plot(six[mask], siy[mask], zorder=40, marker='o', ms=s.psB1m['markersize'], mfc=color, mec='k', mew=0.8) box = dict(boxstyle='square,pad=0.1', fc='white', ec='none') @@ -159,13 +101,13 @@ def plot_nli_diags(ax, s, data, alphax, alphay, xthresh, ythresh, cell_name): title='triangle', handlelength=0.5, handletextpad=0.5, labelspacing=0.2) - kde = gaussian_kde(nliy, 0.15/np.std(nliy, ddof=1)) - nli = np.linspace(0, 8, 100) - pdf = kde(nli) + kde = gaussian_kde(siy, 0.15/np.std(siy, ddof=1)) + si = np.linspace(0, 8, 100) + pdf = kde(si) dax = ax.inset_axes([1.04, 0, 0.3, 1]) dax.show_spines('') - dax.fill_betweenx(nli, pdf, **s.fsB1a) - dax.plot(pdf, nli, clip_on=False, **s.lsB1m) + dax.fill_betweenx(si, pdf, **s.fsB1a) + dax.plot(pdf, si, clip_on=False, **s.lsB1m) def plot_summary_contrasts(axs, s, xthresh, ythresh, cell_name): @@ -173,7 +115,7 @@ def plot_summary_contrasts(axs, s, xthresh, ythresh, cell_name): data = TableData(data_path / 'Apteronotus_leptorhynchus-Punit-models.csv') for i, a in enumerate([0.01, 0.03, 0.1]): print(f'contrast {100*a:2g}%:') - plot_nli_diags(axs[1 + i], s, data, a, a, xthresh, ythresh, cell_name) + plot_si_diags(axs[1 + i], s, data, a, a, xthresh, ythresh, cell_name) print() @@ -182,7 +124,7 @@ def plot_summary_diags(axs, s, xthresh, ythresh, cell_name): data = TableData(data_path / 'Apteronotus_leptorhynchus-Punit-models.csv') for i, a in enumerate([0.01, 0.03, 0.1]): print(f'contrast {100*a:2g}%:') - plot_nli_diags(axs[1 + i], s, data, 0, a, xthresh, ythresh, cell_name) + plot_si_diags(axs[1 + i], s, data, 0, a, xthresh, ythresh, cell_name) if __name__ == '__main__': @@ -211,4 +153,3 @@ if __name__ == '__main__': fig.tag(axs, xoffs=-4.5, yoffs=1.8) axs[1, 0].set_visible(False) fig.savefig() - print() diff --git a/modelsusceptovern.py b/modelsusceptovern.py index d79c916..e713106 100644 --- a/modelsusceptovern.py +++ b/modelsusceptovern.py @@ -3,78 +3,26 @@ from scipy.stats import linregress import matplotlib.pyplot as plt from pathlib import Path from plotstyle import plot_style, labels_params +from plotstyle import noise_files, plot_chi2 +from modelsusceptcontrasts import load_chi2 data_path = Path('data') sims_path = data_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, data_file): - data = np.load(data_file) - n = data['n'] - alpha = data['alpha'] - freqs = data['freqs'] - pss = data['pss'] - prss = data['prss'] - chi2 = np.abs(prss)/0.5/(pss.reshape(1, -1)*pss.reshape(-1, 1)) - ax.set_visible(True) - ax.set_aspect('equal') - i0 = np.argmin(freqs < -300) - i0 = np.argmin(freqs < 0) - i1 = np.argmax(freqs > 300) - 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 = fac*ten - ten *= delta - break - pc = ax.pcolormesh(freqs, freqs, chi2, vmin=0, vmax=vmax, - rasterized=True) - ax.set_title(f'$N=10^{np.log10(n):.0f}$', fontsize='medium') - ax.set_xlim(0, 300) - ax.set_ylim(0, 300) - ax.set_xticks_delta(300) - ax.set_minor_xticks(3) - ax.set_yticks_delta(300) - ax.set_minor_yticks(3) - ax.set_xlabel('$f_1$', 'Hz') - ax.set_ylabel('$f_2$', 'Hz') - 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_yticks_delta(ten) - - def plot_overn(ax, s, files, nmax=1e6, title=False): ns = [] stats = [] for fname in files: data = np.load(fname) - if not 'n' in data: + if not 'nsegs' in data: return - n = data['n'] + n = data['nsegs'] if nmax is not None and n > nmax: continue - alpha = data['alpha'] + noise_frac = data['noise_frac'] + alpha = data['contrast'] freqs = data['freqs'] pss = data['pss'] prss = data['prss'] @@ -84,7 +32,7 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): i1 = np.argmax(freqs > 300) if i1 == 0: i1 = len(freqs) - chi2 = chi2[i0:i1, i0:i1] + chi2 = 1e-4*chi2[i0:i1, i0:i1] # Hz/%^2 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) @@ -94,12 +42,14 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): 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.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') if title: - if 'noise_frac' in data: + if noise_frac < 1: ax.set_title('$c$=0\\,\\%', fontsize='medium') else: ax.set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium') @@ -107,14 +57,12 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): ax.set_xscale('log') ax.set_yscale('log') ax.set_yticks_log(numticks=3) + ax.set_ylim(1e-1, 3e3) + ax.set_minor_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(5e0, 1e4) - ax.set_minor_yticks_log(numticks=5) ax.set_xticks_log(numticks=3) ax.set_minor_xticks_log(numticks=6) ax.set_xlabel('segments') @@ -126,10 +74,15 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): def plot_chi2_overn(axs, s, cell_name): print(cell_name) - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-split-{cell_name}-*.npz'), 1) - for k, n in enumerate([1e1, 1e2, 1e3, 1e6]): - plot_chi2(axs[k], s, files[nums.index(int(n))]) + files, nums = noise_files(sims_path, cell_name) + for k, nsegs in enumerate([1e1, 1e2, 1e3, 1e6]): + freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name, + None, nsegs) + ns = f'$N={n}$' if n < 1000 else f'$N=10^{np.log10(n):.0f}$' + cax = plot_chi2(axs[k], s, freqs, chi2, fcutoff) + if k < len(axs) - 2: + cax.set_ylabel('') + axs[k].set_title(ns, fontsize='medium') plot_overn(axs[-1], s, files) @@ -149,12 +102,10 @@ if __name__ == '__main__': for k in range(len(cells)): plot_chi2_overn(axs[k], s, cells[k]) cell_name = cells[0] - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-split-{cell_name}-*.npz'), 1) + files, nums = noise_files(sims_path, cell_name) plot_overn(axs[-1, 0], s, files, 1e7, True) - for k, alphastr in enumerate(['010', '030', '100']): - files, nums = sort_files(cell_name, - sims_path.glob(f'chi2-noisen-{cell_name}-{alphastr}-*.npz'), 2) + for k, alpha in enumerate([0.01, 0.03, 0.1]): + files, nums = noise_files(sims_path, cell_name, alpha) plot_overn(axs[-1, k + 1], s, files, 1e7, True) for k in range(4): fig.common_yticks(axs[k, :4]) diff --git a/noisesplit.py b/noisesplit.py index b5872a4..19552c0 100644 --- a/noisesplit.py +++ b/noisesplit.py @@ -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) diff --git a/plotstyle.py b/plotstyle.py index b226199..2b64b70 100644 --- a/plotstyle.py +++ b/plotstyle.py @@ -1,5 +1,7 @@ +import numpy as np import matplotlib as mpl import plottools.plottools as pt +from spectral import diag_projection, peak_size from plottools.spines import spines_params from plottools.labels import labels_params from plottools.colors import lighter, darker @@ -16,6 +18,66 @@ def significance_str(p): return '$p<0.001$' +def noise_files(data_path, cell_name, alpha=None): + if alpha is None: + file_pattern = f'{cell_name}-chi2-split-*.npz' + else: + file_pattern = f'{cell_name}-chi2-noise-{1000*alpha:03.0f}-*.npz' + files = sorted(data_path.glob(file_pattern), key=lambda x: x.stem) + 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, fcutoff, rate=None): + ax.set_visible(True) + 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 = 1e-4*chi2[i0:i1, i0:i1] # Hz/%^2 + vquantile = 0.996 + vmax = np.quantile(chi2, vquantile) + 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 + vmax = fac*ten + ten *= 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 rate is not None: + 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 = ax.get_figure().colorbar(pc, cax=cax) + cb.outline.set_color('none') + cb.outline.set_linewidth(0) + cax.set_ylabel(r'$|\chi_2|$', r'Hz/\%$^2$') + cax.set_yticks_delta(ten) + return cax + + def plot_style(): palette = pt.palettes['muted'] lwthick = 1.6