diff --git a/dataoverview.py b/dataoverview.py index 284cf47..92831f6 100644 --- a/dataoverview.py +++ b/dataoverview.py @@ -117,12 +117,12 @@ def plot_corr(ax, data, xcol, ycol, zcol, zmin, zmax, xpdfmax, cmap, color, ha='right', va='center', fontsize='small') # statistics: r, p = pearsonr(xdata, ydata) - ax.text(1, 0.9, f'$R={r:.2f}$ **', ha='right', + ax.text(1, 0.9, f'$R={r:.2f}$', ha='right', + transform=ax.transAxes, fontsize='small') + ax.text(1, 0.77, f'{significance_str(p)}', ha='right', transform=ax.transAxes, fontsize='small') - #ax.text(1, 0.77, f'{significance_str(p)}', ha='right', - # transform=ax.transAxes, fontsize='small') if 'cvbase' in xcol: - ax.text(1, 0.77, f'$n={data.rows()}$', ha='right', + ax.text(1, 0.64, f'$n={data.rows()}$', ha='right', transform=ax.transAxes, fontsize='small') print(f' correlation {xcol:<8s} - {ycol}: r={r:5.2f}, p={p:.2g}') return cax @@ -154,6 +154,15 @@ def si_stats(title, data, sicol, si_thresh, nsegscol): 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}%') + cols = ['cvbase', 'respmod2', 'ratebase', 'vsbase', 'serialcorr1', 'burstfrac', 'ratestim', 'cvstim'] + for i in range(len(cols)): + for j in range(i + 1, len(cols)): + xcol = cols[i] + ycol = cols[j] + if xcol not in data or ycol not in data: + continue + r, p = pearsonr(data[xcol], data[ycol]) + print(f' correlation {xcol:<11s} - {ycol:<11s}: r={r:5.2f}, p={p:.5f}') def plot_cvbase_si_punit(ax, data, ycol, si_thresh, color): @@ -204,6 +213,19 @@ def plot_rmod_si_punit(ax, data, ycol, si_thresh, color): cax.set_ylabel('CV$_{\\rm base}$') +def plot_rate_si_punit(ax, data, ycol, si_thresh, color): + ax.set_xlabel('Baseline rate $r$', 'Hz') + ax.set_xlim(0, 700) + ax.set_xticks_delta(200) + ax.set_ylabel('SI($r$)') + ax.set_ylim(0, 6.5) + ax.set_yticks_delta(2) + examples = punit_examples if 'stimindex' in data else model_examples + cax = plot_corr(ax, data, 'ratebase', ycol, 'cvbase', 0, 1.5, 0.016, + 'coolwarm', color, si_thresh, *examples) + cax.set_ylabel('CV$_{\\rm base}$') + + def plot_cvbase_si_ampul(ax, data, ycol, si_thresh, color): ax.set_xlabel('CV$_{\\rm base}$') ax.set_xlim(0, 0.2) @@ -248,6 +270,19 @@ def plot_rmod_si_ampul(ax, data, ycol, si_thresh, color): cax.set_ylabel('CV$_{\\rm base}$') cax.set_yticks_delta(0.1) + +def plot_rate_si_ampul(ax, data, ycol, si_thresh, color): + ax.set_xlabel('Baseline rate $r$', 'Hz') + ax.set_xlim(50, 200) + ax.set_xticks_delta(50) + ax.set_ylabel('SI($r$)') + ax.set_ylim(0, 10) + ax.set_yticks_delta(2) + cax = plot_corr(ax, data, 'ratebase', ycol, 'cvbase', 0, 0.2, 0.06, + 'coolwarm', color, si_thresh, *ampul_examples) + cax.set_ylabel('CV$_{\\rm base}$') + cax.set_yticks_delta(0.1) + if __name__ == '__main__': punit_model = TableData(data_path / @@ -304,8 +339,10 @@ if __name__ == '__main__': s.model_color2) plot_rmod_si_punit(axs[0, 1], punit_model, 'dsinorm100', si_thresh, s.model_color2) - plot_cvstim_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh, - s.model_color2) + #plot_cvstim_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh, + # s.model_color2) + plot_rate_si_punit(axs[0, 2], punit_model, 'dsinorm100', si_thresh, + s.model_color2) print() si_stats('P-unit data:', punit_data, 'sinorm' + si, si_thresh, @@ -316,8 +353,10 @@ if __name__ == '__main__': s.punit_color2) plot_rmod_si_punit(axs[1, 1], punit_data, 'sinorm' + si, si_thresh, s.punit_color2) - plot_cvstim_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh, - s.punit_color2) + #plot_cvstim_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh, + # s.punit_color2) + plot_rate_si_punit(axs[1, 2], punit_data, 'sinorm' + si, si_thresh, + s.punit_color2) print() si_stats('Ampullary data:', ampul_data, 'sinorm' + si, si_thresh, @@ -328,8 +367,10 @@ if __name__ == '__main__': s.ampul_color2) plot_rmod_si_ampul(axs[2, 1], ampul_data, 'sinorm' + si, si_thresh, s.ampul_color2) - plot_cvstim_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh, - s.ampul_color2) + #plot_cvstim_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh, + # s.ampul_color2) + plot_rate_si_ampul(axs[2, 2], ampul_data, 'sinorm' + si, si_thresh, + s.ampul_color2) print() fig.common_xticks(axs[:2, 0]) diff --git a/modelsusceptcontrasts.py b/modelsusceptcontrasts.py index 12e4d03..892bf69 100644 --- a/modelsusceptcontrasts.py +++ b/modelsusceptcontrasts.py @@ -30,20 +30,21 @@ def load_chi2(file_path, cell_name, contrast=None, n=None): return freqs, chi2, fcutoff, contrast, n -def plot_chi2_contrasts(axs, s, cell_name): +def plot_chi2_contrasts(axs, s, cell_name, vmax): 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}') freqs, chi2, fcutoff, contrast, n = load_chi2(sims_path, cell_name) - cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate) + cax = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate, vmax) 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) + cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate, vmax) + vmax /= 2 if alpha < 0.1: cax.set_ylabel('') axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%', fontsize='medium') @@ -142,8 +143,9 @@ if __name__ == '__main__': fig.subplots_adjust(leftm=7, rightm=9, topm=2, bottomm=3.5, wspace=1, hspace=0.5) print('Example cells:') + vmax = [2, 6, 8, 40] for k in range(len(model_cells)): - plot_chi2_contrasts(axs[k], s, model_cells[k]) + plot_chi2_contrasts(axs[k], s, model_cells[k], vmax[k]) for k in range(4): fig.common_yticks(axs[k, :]) fig.common_xticks(axs[:4, k]) diff --git a/modelsusceptlown.py b/modelsusceptlown.py index 5145b07..5a2f32c 100644 --- a/modelsusceptlown.py +++ b/modelsusceptlown.py @@ -14,7 +14,7 @@ data_path = Path('data') sims_path = data_path / 'simulations' -def plot_chi2_contrasts(axs, s, cell_name, nsegs=None): +def plot_chi2_contrasts(axs, s, cell_name, nsegs=None, vmax=None): d = sims_path / f'{cell_name}-baseline.npz' data = np.load(d) rate = float(data['rate']) @@ -23,14 +23,18 @@ def plot_chi2_contrasts(axs, s, cell_name, nsegs=None): 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 = plot_chi2(axs[0], s, freqs, chi2, fcutoff, rate, vmax) 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) + cax = plot_chi2(axs[k + 1], s, freqs, chi2, fcutoff, rate, vmax) + if n < 1000: + vmax /= 10 + else: + vmax /= 4 if alpha < 0.1: cax.set_ylabel('') axs[k + 1].set_title(f'$c$={100*alpha:g}\\,\\%, {ns}', @@ -139,8 +143,8 @@ if __name__ == '__main__': for ax in axs.flat: ax.set_visible(False) print('Example cells:') - plot_chi2_contrasts(axs[0], s, model_cell) - plot_chi2_contrasts(axs[1], s, model_cell, nsmall) + plot_chi2_contrasts(axs[0], s, model_cell, None, 40) + plot_chi2_contrasts(axs[1], s, model_cell, nsmall, 600) for k in range(2): fig.common_yticks(axs[k, :]) for k in range(4): diff --git a/modelsusceptovern.py b/modelsusceptovern.py index e713106..4120421 100644 --- a/modelsusceptovern.py +++ b/modelsusceptovern.py @@ -75,7 +75,7 @@ def plot_overn(ax, s, files, nmax=1e6, title=False): def plot_chi2_overn(axs, s, cell_name): print(cell_name) files, nums = noise_files(sims_path, cell_name) - for k, nsegs in enumerate([1e1, 1e2, 1e3, 1e6]): + for k, nsegs in enumerate([1e2, 1e3, 1e4, 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}$'