import numpy as np import matplotlib.pyplot as plt from scipy.stats import pearsonr, linregress, gaussian_kde from thunderlab.tabledata import TableData from pathlib import Path 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 '2012-12-13-ao-invivo-1', # triangle '2012-12-20-ac-invivo-1', # weak border triangle '2013-01-08-ab-invivo-1'] # no triangle data_path = Path('data') sims_path = data_path / 'simulations' 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)) return freqs, chi2, fcutoff, contrast, n 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, 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, vmax) vmax /= 2 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['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, 'dsinorm'], clip_on=False, zorder=30, label='strong', **s.psA1m) mask = data['border'] > 0.5 ax.plot(data[mask, 'cvbase'], data[mask, 'dsinorm'], zorder=20, label='weak', **s.psA2m) ax.plot(data['cvbase'], data['dsinorm'], clip_on=False, zorder=10, label='none', **s.psB1m) for cell_name in cells: mask = data[:, 'cell'] == cell_name color = s.psB1m['color'] if data[mask, 'border']: color = s.psA2m['color'] elif data[mask, 'triangle']: color = s.psA1m['color'] 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_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($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, ha='right', fontsize='small') if alpha == 0: ax.legend(loc='upper left', bbox_to_anchor=(1.15, 1.05), title='triangle', handlelength=0.5, handletextpad=0.5, labelspacing=0.2) 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(si, pdf, **s.fsB1a) dax.plot(pdf, si, clip_on=False, **s.lsB1m) def plot_summary_contrasts(axs, s, cells): nli_thresh = 1.2 data = TableData(data_path / 'Apteronotus_leptorhynchus-Punit-models.csv') plot_si_cv(axs[0], s, data, 0, cells) print('noise split:') cdata = data[data['contrast'] == 0, :] nli_split = cdata['dsinorm'] print(f' mean SI = {np.mean(nli_split):.2f}, stdev = {np.std(nli_split):.2f}') print(f' triangle cells have SI >= {np.min(nli_split[cdata["triangle"] > 0.5]):.2f}') ntriangle = np.sum(cdata['triangle'] > 0.5) print(f' triangle cells: {ntriangle:2d} ({100*ntriangle/len(cdata):.1f}%)') nborder = np.sum(cdata['border'] > 0.5) print(f' border cells: {nborder:2d} ({100*nborder/len(cdata):.1f}%)') nother = len(cdata) - ntriangle - nborder print(f' other cells: {nother:2d} ({100*nother/len(cdata):.1f}%)') n = np.sum(nli_split > nli_thresh) print(f' {n} cells ({100*n/len(nli_split):.1f}%) have SI > {nli_thresh:.1f}:') for name, cv in cdata[nli_split > nli_thresh, ['cell', 'cvbase']].row_data(): print(f' {name:<22} CV={cv:4.2f}') print() for i, a in enumerate([0.01, 0.03, 0.1]): plot_si_cv(axs[1 + i], s, data, a, cells) print(f'contrast {100*a:2g}%:') cdata = data[data['contrast'] == a, :] 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}') n = np.sum(nli > nli_thresh) print(f' triangle cells have SI >= {np.min(nli[cdata["triangle"] > 0.5]):.2f}') print(f' {n} cells ({100*n/len(nli):.1f}%) have SI > {nli_thresh:.1f}:') for name, cv, respmod in cdata[nli > nli_thresh, ['cell', 'cvbase', 'respmod2']].row_data(): print(f' {name:<22} CV={cv:4.2f} ' f'response modulation={respmod:4.0f}Hz') print() print('overall lowest baseline CV:', ' '.join([f'{cv:.2f}' for cv in np.unique(data['cvbase'])[:5]])) if __name__ == '__main__': s = plot_style() #labels_params(xlabelloc='right', ylabelloc='top') fig, axs = plt.subplots(5, 4, cmsize=(s.plot_width, 0.95*s.plot_width), height_ratios=[3, 3, 3, 3, 0, 2]) 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], vmax[k]) for k in range(4): fig.common_yticks(axs[k, :]) fig.common_xticks(axs[:4, k]) print() plot_summary_contrasts(axs[4], s, model_cells) fig.common_yticks(axs[4, :]) fig.tag(axs, xoffs=-4.5, yoffs=1.8) fig.savefig() print()