167 lines
5.8 KiB
Python
167 lines
5.8 KiB
Python
import numpy as np
|
|
from scipy.stats import linregress
|
|
import matplotlib.pyplot as plt
|
|
from pathlib import Path
|
|
from plotstyle import plot_style, labels_params
|
|
|
|
|
|
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']
|
|
dt_fix = 1 # 0.0005
|
|
prss = np.abs(data['prss'])/dt_fix*0.5/np.sqrt(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]
|
|
prss = prss[i0:i1, i0:i1]
|
|
vmax = np.quantile(prss, 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, prss, vmin=0, vmax=vmax,
|
|
cmap='viridis', 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:
|
|
return
|
|
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], '0.5', lw=1, zorder=50, label='99.8\\%')
|
|
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[:, 8], '0.0')
|
|
if title:
|
|
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(1e1, nmax)
|
|
ax.set_xscale('log')
|
|
ax.set_yscale('log')
|
|
ax.set_yticks_log(numticks=3)
|
|
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')
|
|
ax.set_ylabel('$|\\chi_2|$ [Hz]')
|
|
if alpha == 0.10:
|
|
ax.legend(loc='upper left', bbox_to_anchor=(1.4, 1.3),
|
|
markerfirst=False, title='$|\\chi_2|$ percentiles')
|
|
|
|
|
|
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))])
|
|
plot_overn(axs[-1], s, files)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
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
|
|
s = plot_style()
|
|
fig, axs = plt.subplots(6, 6, cmsize=(s.plot_width, 0.9*s.plot_width),
|
|
width_ratios=[1, 1, 1, 1, 0, 1],
|
|
height_ratios=[1, 1, 1, 1, 0, 1])
|
|
fig.subplots_adjust(leftm=8, rightm=0.5, topm=2, bottomm=4,
|
|
wspace=1, hspace=0.8)
|
|
for ax in axs.flat:
|
|
ax.set_visible(False)
|
|
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)
|
|
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)
|
|
plot_overn(axs[-1, k + 1], s, files, 1e7, True)
|
|
for k in range(4):
|
|
fig.common_yticks(axs[k, :4])
|
|
fig.common_xticks(axs[:4, k])
|
|
fig.common_xticks(axs[:4, -1])
|
|
fig.align_ylabels(axs[:4, -1], dist=12)
|
|
fig.common_yticks(axs[-1, :4])
|
|
fig.tag(axs, xoffs=-2.5, yoffs=1.8)
|
|
fig.savefig()
|