import os
osjoin = os.path.join # an alias for convenience
from pdb import set_trace
import json
from hetero.dynamics.li import activation_func
from hetero import utils
# from hetero.dataloader import ds_loader, spike_loader
from tqdm import tqdm
import numpy as np
import pandas as pd
from scipy.signal import find_peaks, periodogram, correlation_lags, correlate
from scipy.sparse import load_npz
from scipy.linalg import svd
import matplotlib as mpl
from matplotlib import ticker
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.colors import (LogNorm, Normalize, to_rgb, ListedColormap,
BoundaryNorm)#, LinearSegmentedColormap
from mpl_toolkits.axes_grid1 import make_axes_locatable, ImageGrid
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import matplotlib.patheffects as patheffects
from cycler import cycler
# import matplotlib.animation as animation
# from matplotlib.offsetbox import AnchoredText
# import matplotlib.cm as cm
# import scienceplots
# plt.style.use(['latex-sans', 'std-colors'])
# plt.style.use(['scienceplots'])
# plt.style.use(['no-latex'])
plt.rcParams['savefig.bbox'] = 'tight'
plt.rcParams['savefig.dpi'] = 150
plt.rcParams.update({'text.latex.preamble': r'\usepackage{amsfonts} \usepackage{amsmath} \usepackage{sfmath}'})
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif"
})
import warnings
warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)
mpl.colormaps.unregister("hcmap")
# blues = mpl.colormaps['Blues'].resampled(5)
hetero_colors = [
"#ae3832",
"#9fbfe2",
"#5d95d0",
"#336090",
]
# hetero_colors = blues(np.linspace(0, 1, 5))
# hetero_colors[0] = np.array([146./255, 19./255, 12./255, 1]) #
# hetero_colors[0] = np.array([1, 0, 0., 1])
hcmap = mpl.colors.ListedColormap(hetero_colors, name="hcmap")
mpl.colormaps.register(cmap=hcmap)
# mpl.colormaps.unregister("my_mix")
# blues = mpl.colormaps['Blues'].resampled(1000)
# grays = mpl.colormaps['Greys'].resampled(1000)
# bcolors = blues(np.linspace(0, 0.9, int((1 - np.exp(-1))*1000)))
# gcolors = grays(np.linspace(0, 0.65, int(np.exp(-1)*1000)))
# mixed_colores = np.vstack([gcolors[::-1], bcolors])
# my_mixed = mpl.colors.ListedColormap(mixed_colores, name="my_mix")
# mpl.colormaps.register(cmap=my_mixed)
EPS = 1e-12
NET_ID_MAP = {
'0400': 'IN', # Instanteneous - None
'1400': 'FN', # Fast - None
'1410': 'FL', # Fast - Low
'1420': 'FM', # Fast - Moderate
'1430': 'FH', # Fast - High
'1440': 'FvH', # Fast - Very High
'2400': 'MN', # Moderate - None
'2410': 'ML', # Moderate - Low
'2420': 'MM', # Moderate - Moderate
'2430': 'MH', # Moderate - High
'2440': 'MvH', # Moderate - High
'3400': 'SN', # Slow - None
'3410': 'SL', # Slow - Low
'3420': 'SM', # Slow - Moderate
'3430': 'SH', # Slow - High
'3440': 'SvH', # Slow - High
'4400': '-N', # Fixed
}
HETERO_LABEL_MAP = {
'4':'Very high',
'3':'High',
'2':'Medium',
'1':'Low',
'0':'None',
}
TAU_MAP = {
'0': '0',
'1': r'$\tau_0 / \sqrt{10}$',
'2': r'$\tau_0$',
'3': r'$\tau_0\sqrt{10}$',
'4': r'$\infty$',
}
STD_CMAP = ['Blues', 'Greens', 'Oranges', 'Reds', 'Purples', 'Greys']
CAPACITY_PARAMS = ['k', 'mono_degs'] # test size, powers
FORECAST_PARAMS = ['k', 'l', 'm'] # test size, forecast horizon (delta), memory
MFS = 18 # maximum font size
## UTILS
[docs]
def polish(fig, axs, lw=1, tight=True):
if type(axs) == type({}):
axs = list(axs.values())
fig.align_xlabels(axs)
fig.align_ylabels(axs)
fig.align_titles(axs)
axs = np.array(axs)
for ax in axs.ravel():
simpleaxis(ax, lw)
if tight:
plt.tight_layout()
[docs]
def simpleaxis(ax, lw=2):
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(lw)
ax.xaxis.set_tick_params(width=lw)
ax.yaxis.set_tick_params(width=lw)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
[docs]
def get_cax(ax, width=3., dist=.05):
"""
Returns a human-friendly colorbar axis from a given ax, to be used with
matplotlib's colorbar method.
"""
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
cax = inset_axes(ax, width=f"{width}%", height="100%",
loc='lower left', bbox_to_anchor=(1 + dist, 0., 1, 1),
bbox_transform= ax.transAxes, borderpad=0)
return cax
[docs]
def aggregate_scores(res_path, sims=[], sim_pattern=''):
## COLLECT SIM FOLDERS AND NET_IDS
if len(sims) == 0:
sims = sorted([f.path.split('/')[-1] for f in os.scandir(res_path)
if (f.is_dir() and sim_pattern in f.path)])
net_ids = sorted([f.name for f in os.scandir(osjoin(res_path, sims[0])) if f.is_dir()])
## COLLECT ALL SCORES IN EACH SIM FOLDER
# FILL ALL CONTAINERS
scr_mu = [] # score mean
scr_std= [] # score std
ros = [] # readout setup
Ns = [] # N neurons
ps = [] # sparisty
fs = [] # E-fraction
sigmas = [] # wdist sigma_0
distros = [] # heterog distribution
Ntrns = [] # num samples in the train set
n_timesteps = []
Js = [] # recurrent gains
Jns = [] # noise gain
Jus = [] # input gain
dyns = [] # dynamic
nets = [] # network id
simss = [] # sim name
rlzs = []
for sim in sims:
# print(sim)
md = utils.read_metadata(osjoin(res_path, sim))
N = md['net']['N']
p = md['net']['p']
f = md['net']['f']
sigma = md['net']['sig0']
Ntrn = md['Ntrn']
n_timestep = md['n_timesteps']
dim = md['dim']
task = list(md['task'].keys())[0]
distro = md['heterog']['distro']
net_ids = sorted([f.name for f in os.scandir(osjoin(res_path, sim)) if f.is_dir()])
for net_id in net_ids:
# print(net)
scrs = sorted([f.name for f in os.scandir(osjoin(res_path, sim, net_id)) if 'tst' in f.name])
# scrs = [scr for scr in scrs if '_euler.' not in scr] # exclude the imprecise euler scores
for scr in scrs:
print(sim, net_id, scr)
scr_df = pd.read_csv(osjoin(res_path, sim, net_id, scr))
# if task == 'tayloramus':
# else:
# scr_val.append(scr_df.set_index(['digit','speaker']).mean(axis=1).values)
rlz = int([x.split('rlz')[-1] for x in scr.split('-') if 'rlz' in x][0])
rlzs.append(rlz)
scr_mu.append(scr_df.set_index(['dim','delta','deg']).mean(axis=1).values)
scr_std.append(scr_df.set_index(['dim','delta','deg']).std(axis=1).values)
Ns.append(N)
ps.append(p)
fs.append(f)
sigmas.append(sigma)
distros.append(distro)
Ntrns.append(Ntrn)
n_timesteps.append(n_timestep)
_Jn, _J, _Ju = get_Js(scr)
Js.append(_J)
Jus.append(_Ju)
Jns.append(_Jn)
dyns.append('LIF' if 'LIF' in scr else 'LI')
nets.append(net_id)
simss.append(sim)
if 'LIF' in scr:
ro = 'd10' if 'd10' in scr else 'd50'
else:
ro = 'possig1'
ros.append(ro)
## BUILD THE MEAN AND STD DATAFRAMES
# NOTE: tasks, alongside other parameters are as columns
df0 = pd.DataFrame({'sim':simss, 'rlz':rlzs, 'N': Ns, 'p': ps, 'f': fs, 'sigma': sigmas,
'trn': Ntrns, 'n_timestep': n_timesteps,
'distro': distros, 'dyn': dyns, 'ro': ros,
'J': Js, 'Jn': Jns, 'Ju': Jus, 'net': nets})
df0['h'] = df0.net.map(lambda x: int(x[1]))
df0['n'] = df0.net.map(lambda x: int(x[0]))
# we use tasks as index. Let's extract the tasks param combination
# (dim, delta, deg) from a score file with maximal amount of tasks
# (e.g., lorenz, but not one-dimensional inputs like narma or sine)
scr_df = pd.read_csv(osjoin(res_path, simss[0], net_id, scr))
task_combi = scr_df.set_index(['dim', 'delta', 'deg']).index
df_mu = pd.DataFrame(scr_mu, columns=task_combi)
df_std = pd.DataFrame(scr_std, columns=task_combi)
df_mu = pd.concat([df0, df_mu], axis=1)
df_std = pd.concat([df0, df_std], axis=1)
## METLING THE DATAFRAME BASED ON NETWORK SETUP PARAM (IDS)
## AND DIFFERENT TASKS (VALUES)
tcols = df_mu.columns[17:].values.tolist() # 17 is the first column of tasks...
df = df_mu.melt(id_vars = [col for col in df_mu.columns if col not in tcols],
value_vars=tcols, var_name='task', value_name='score').reset_index(drop=True)
df['score_std'] = df_std.melt(id_vars = [col for col in df_mu.columns if col not in tcols],
var_name='task', value_name='score_std').reset_index(drop=True)['score_std']
## REMOVE TASKS THAT DON'T EXIST FOR ONE-DIM INPUTS
# we find valid_idx based on the mean score, rather than all
# columns. Because if only one realization or one trial is
# run, then std would be nan and dropna will remove all rows.
df = df.reset_index(drop=True)
valid_idxs = df['score'].dropna().index
df = df.loc[valid_idxs].reset_index(drop=True) # to remove those task that don't exist in 1-dimensional tasks
## ADD COMPLEXITY AND COMPLEIXTY TIERS
bins = [-1e-4, 1./3, 2./3, 1.0+1e-4]
labels = ['easy','intermediate','difficult']
df['comp'] = -1.
df = df.set_index('task')
for sim in df.sim.unique():
cmplx = pd.read_csv(osjoin(res_path, sim, 'complexities.csv'), index_col=0)
cmplx['task'] = [(int(row.dim), int(row.delta), int(row.deg)) for _, row in cmplx.iterrows()]
cmplx = cmplx.set_index('task')
df.loc[df.sim == sim, 'comp'] = cmplx.loc[df[df.sim == sim].index, 'complexity']
df['tier'] = pd.cut(df.comp, bins=bins, labels=labels)
df= df.reset_index()
return df
# def bilinear_color_interpolation(x, y):
# # c00 = np.array([.9, .9, .9]) # white at (0, 0)
# # c10 = np.array(to_rgb('dodgerblue')) # Blue at (1,0)
# # c01 = np.array(to_rgb('forestgreen')) # Green at (0, 1)
# # c11 = np.array(to_rgb('crimson')) # Red at (1, 1)
# c00 = np.array(to_rgb('lightsteelblue')) # white at (0, 0)
# c10 = np.array(to_rgb('tab:blue')) # Blue at (0, 1)
# c01 = np.array(to_rgb('tab:orange')) # Green at (1, 0)
# c11 = np.array(to_rgb('k')) # Red at (1, 1)
# if len(x)==1:
# color = (1-x)*(1-y) * c00 + x*(1-y) * c10 + (1-x)*y * c01 + x*y * c11
# else:
# locs = np.array([(1-x)*(1-y), x*(1-y), (1-x)*y, x*y]).T
# color = locs @ np.array([c00, c10, c01, c11])
# return color
[docs]
def append_subplot(fig, axes, position='bottom'):
# Get the current layout dimensions
nrows = max(ax.get_subplotspec().rowspan.stop for ax in axes)
ncols = max(ax.get_subplotspec().colspan.stop for ax in axes)
# Create a new GridSpec with an additional row or column
if position in ['bottom', 'top']:
gs = GridSpec(nrows + 1, ncols, figure=fig)
elif position in ['left', 'right']:
gs = GridSpec(nrows, ncols + 1, figure=fig)
else:
raise ValueError("Position must be 'bottom', 'top', 'left', or 'right'")
# Adjust the positions of the existing subplots
for ax in axes:
pos = ax.get_subplotspec()
if position == 'bottom':
ax.set_subplotspec(gs[pos.rowspan.start:pos.rowspan.stop, pos.colspan.start:pos.colspan.stop])
elif position == 'top':
ax.set_subplotspec(gs[pos.rowspan.start + 1:pos.rowspan.stop + 1, pos.colspan.start:pos.colspan.stop])
elif position == 'left':
ax.set_subplotspec(gs[pos.rowspan.start:pos.rowspan.stop, pos.colspan.start + 1:pos.colspan.stop + 1])
elif position == 'right':
ax.set_subplotspec(gs[pos.rowspan.start:pos.rowspan.stop, pos.colspan.start:pos.colspan.stop])
# Add the new subplot
if position == 'bottom':
new_ax = fig.add_subplot(gs[nrows, :])
elif position == 'top':
new_ax = fig.add_subplot(gs[0, :])
elif position == 'left':
new_ax = fig.add_subplot(gs[:, 0])
elif position == 'right':
new_ax = fig.add_subplot(gs[:, ncols])
return new_ax
[docs]
def plot_activity(t, v, activation, res_path, name, **kwargs):
fig, axs = plt.subplots(1,1,)
for idx, v_ in enumerate(v):
plt.plot(t , activation_func(v_, activation), **kwargs)
# plt.legend()
# plt.ylim(-2,2)
polish(fig, axs)
plt.savefig(osjoin(res_path, f'activity_{name}.png'), dpi=150,)
plt.close()
[docs]
def plot_weights(t, w, res_path, name, **kwargs):
fig, axs = plt.subplots(1,1,)
for idx, w_ in enumerate(w):
plt.plot(t , w_, **kwargs)
polish(fig, axs)
plt.savefig(osjoin(res_path, f'activity_{name}.png'), dpi=150,)
plt.close()
# def plot_score_setup(data_path, res_path, net_id, suffix, make_axs = True,
# # normalize_l=None, feature_size=None,
# test=True):
# if test:
# set_name = 'test'
# else:
# set_name = 'train'
# df = pd.read_csv(osjoin(data_path, net_id, f'scores{suffix}.csv'))#, converters={'pows': pd.eval})
# ks = np.sort(df.k.unique())
# ls = np.sort(df.l.unique())
# # if normalize_l is not None:
# # ls = ls*1. / normalize_l
# # if feature_size is not None:
# # ks = ks * feature_size
# stats = df.groupby(['l','k']).agg(['mean', 'std'])
# df_mean, df_std = utils.extract_mean_std_from_groupped(stats)
# dims = utils.get_dim(df_mean)
# if make_axs:
# fig, axs = plt.subplots(1, dims, figsize=(5*dims,4), sharey=True, sharex=True)
# if dims == 1:
# axs = np.array([axs])
# else:
# fig, axs = None, None
# return df_mean, df_std, dims, ks, ls, fig, axs, set_name
# def plot_score_forecast_line_by_l(data_path, res_path, net_id,
# normalize_l=None, feature_size=None,
# test=True, **kwargs):
# df_mean, df_std, dims, ks, ls, fig, axs, set_name = plot_score_setup(
# data_path, res_path, net_id, suffix= '_f',
# # normalize_l=normalize_l, feature_size=feature_size
# )
# # normalize delta and adjust the xlabel
# if normalize_l is not None:
# ls = ls/(1.*normalize_l)
# xlabel = r'$\boldsymbol{\Delta} / \boldsymbol{\tau_0}$'
# else:
# xlabel = r'$\boldsymbol{\Delta}$'
# for dim in range(dims):
# ax = axs[dim]
# # this is different from line_by_k
# _mean = utils.get_score_sub_df(df_mean, set_name+f'_dim_{dim+1}').unstack(level=1)
# _std = utils.get_score_sub_df(df_std, set_name+f'_dim_{dim+1}').unstack(level=1)
# for k_id, k in enumerate(ks):
# ax.plot(ls, _mean.iloc[:,k_id], label = f'k={k}', color=f'C{k_id}')
# ax.fill_between(ls,
# _mean.iloc[:,k_id]-_std.iloc[:,k_id],
# _mean.iloc[:,k_id]+_std.iloc[:,k_id],
# color=f'C{k_id}', alpha=0.2)
# for dim in range(dims):
# ax = axs[dim]
# ax.legend(ncols=1, loc='lower left')
# ax.set_title(r'$\textbf{dim = }$'+fr'$\textbf{dim+1}$')
# ax.set_xscale('log')
# ax.set_xlabel(xlabel)
# simpleaxis(ax, None)
# axs[0].set_ylabel(r'$\boldsymbol{R^2 (\Delta)}$')
# axs[0].set_ylim(0,1)
# plt.savefig(osjoin(res_path, f'score_l_{set_name}_{net_id}.svg'),
# bbox_inches='tight')
# plt.close()
# def plot_score_forecast_line_by_k(data_path, res_path, net_id,
# base_stride=None, feature_size=None,
# test=True, **kwargs):
# df_mean, df_std, dims, ks, ls, fig, axs, set_name = plot_score_setup(
# data_path, res_path, net_id, suffix = '_f',
# # base_stride = base_stride, feature_size = feature_size,
# )
# # Create a custom color cycler
# colormap = plt.cm.Blues
# num_colors = len(ls)
# colors_from_colormap = [colormap(i) for i in np.linspace(0, 1, num_colors)]
# custom_cycler = cycler(color=colors_from_colormap)
# for dim in range(dims):
# ax = axs[dim]
# ax.set_prop_cycle(custom_cycler)
# # this is different from line_by_l
# _mean = utils.get_score_sub_df(df_mean,f'test_dim_{dim+1}').unstack(level=0)
# _std = utils.get_score_sub_df(df_std,f'test_dim_{dim+1}').unstack(level=0)
# for l_id, l in enumerate(ls):
# ax.plot(ks, _mean.iloc[:,l_id],)
# ax.fill_between(ks,
# _mean.iloc[:,l_id]-_std.iloc[:,l_id],
# _mean.iloc[:,l_id]+_std.iloc[:,l_id],
# alpha=0.05,)
# # colorbar
# cax = make_axes_locatable(axs[-1]).append_axes("right", size="2.5%", pad=0.1)
# custom_ticks = [0,1]
# custom_tick_labels = ['imm.', 'far']
# # Create the colorbar with custom settings
# cbar = plt.colorbar(plt.cm.ScalarMappable(norm=Normalize(), cmap=colormap),
# cax=cax, ticks=custom_ticks,)
# cbar.ax.set_yticklabels(custom_tick_labels) # Set custom tick labels
# cbar.set_label('Future Horizon',)# pad=-10) # Set custom colorbar label
# for dim in range(dims):
# ax = axs[dim]
# ax.set_xlabel(r'$\textbf{Training set size (normalized)}$')
# ax.set_title(r'$\textbf{dim = }$'+fr'$\textbf{dim+1}$')
# ax.set_xlim(ks[0], ks[-1])
# simpleaxis(ax, None)
# axs[0].set_ylabel(r'$\boldsymbol{R^2}$')
# axs[0].set_ylim(0,1)
# plt.savefig(osjoin(res_path, f'score_k_{set_name}_{net_id}.svg'),
# bbox_inches='tight')
# plt.close()
# def plot_score_bar_on_sweep(data_path, net_id_extractor, v=None, w=None, hv=None, hw=None,
# normalize_l=None):
# cbar_labels = ['fast', 'slow']
# v_str, w_str = 'x', 'x'
# hv_str , hw_str = '', ''
# sort_by_heter = False
# pad = -15
# cbar_y = 0.55
# # s-sweep (var=0)
# if (v is not None) and (w is None):
# hv, hw = 0, 0 # overwrites the hv, hw arguments
# v_str = str(v)
# sweep_type = 's'
# cbar_title = r'$\boldsymbol{\mathbb E[\tau_s]}$'
# cmap = plt.cm.Greens_r
# # n-sweep (var=0)
# elif (v is None) and (w is not None):
# hv, hw = 0, 0 # overwrites the hv, hw arguments
# w_str = str(w)
# sweep_type = 'n'
# cbar_title = r'$\boldsymbol{\mathbb E[\tau_n]}$'
# cmap = plt.cm.Blues_r
# # h-sweep
# else:
# pad = -15
# cbar_labels = ['hom.', 'het.']
# v_str = str(v)
# w_str = str(w)
# hv_str , hw_str = 'x', 'x'
# # hn-sweep (var w = constant)
# if (hw is not None) and (hv is None):
# hw_str = str(hw)
# sweep_type = 'hn'
# cbar_title = r'$\boldsymbol{Var[\tau_n]}$'
# cmap = plt.cm.Reds_r
# # hs-sweep (var v = constant)
# elif (hw is None) and (hv is not None):
# hv_str = str(hv)
# sweep_type = 'hs'
# cbar_title = r'$\boldsymbol{Var[\tau_s]}$'
# cmap = plt.cm.Wistia_r
# else:
# print('Warning: You have requested all heterogenous combinations. It will be a mess!')
# sort_by_heter = True
# hv_str = str(hv_str)
# hw_str = str(hw_str)
# sweep_type = 'h'
# cbar_title = r'$\boldsymbol{Var[\tau_s]}$'+'\n + \n'+ r'$\boldsymbol{Var[\tau_n]}$'
# cbar_y = .75
# cmap = plt.cm.viridis
# net_ids = net_id_extractor(v,w,hv,hw)
# _df = pd.read_csv(osjoin(data_path, '0000', 'scores.csv'))
# ls = np.array(sorted(_df.l.unique()))
# dims = utils.get_dim(_df)
# if normalize_l is not None:
# ls = ls/(1.*normalize_l)
# xlabel = r'$\boldsymbol{\Delta} / \boldsymbol{\tau_0}$'
# else:
# xlabel = r'$\boldsymbol{\Delta}$'
# fig, axs = plt.subplots(dims, 1, figsize=(4.5*dims, 5), sharex=True, sharey=True)
# for ax in axs.ravel():
# ax.set_prop_cycle(cycler('color', cmap(np.linspace(0,0.7,len(net_ids)))))
# multiplier = 0
# width = 0.6/len(net_ids) # the width of the bars
# x = np.arange(len(ls)) - 0.3
# for net_id in net_ids:
# _, df_mean, df_std, _, ks, ls = utils.get_scores_of_max_k(data_path, net_id)
# offset = width * multiplier
# for dim in range(dims):
# score = df_mean[f'score_test_dim_{dim+1}'].values
# rects = axs[dim].bar(x + offset, score, width, label=net_id)
# # axs[dim].bar_label(rects, rotation=90, fontsize=5)
# multiplier += 1
# axs[dim].set_xticks(x, ls, rotation=90);
# axs[dim].set_ylim(0,1)
# xlims = x[0]-0.75, x[-1]+0.75
# axs[dim].set_xlim(*xlims)
# # colorbar
# for ax in axs:
# cax = make_axes_locatable(ax).append_axes("right", size="1%", pad=0.1)
# cbar = plt.colorbar(plt.cm.ScalarMappable(norm=Normalize(), cmap=cmap), cax=cax, ticks=[0,1])
# cbar.set_label(cbar_title, labelpad=pad, )# y=cbar_y, # rotation=0,) #
# cbar.ax.set_yticklabels(cbar_labels)
# # ylabel
# for ax in axs:
# ax.set_ylabel(r'$\boldsymbol{R^2 (\Delta)}$')
# ax.set_ylim(0,1)
# ax.set_xlabel(xlabel)
# # xlabel
# # for ax in axs:
# # ax.set_xlabel(xlabel)
# # titles
# axs[0].set_title(r'$\textbf{Periodic }$')
# axs[1].set_title(r'$\textbf{"Mildly" chaotic}$')
# axs[2].set_title(r'$\textbf{"Strongly" chaotic}$')
# for ax in axs.ravel():
# ax.hlines(0, *xlims, color='k', linestyle='--')
# # ax.legend()
# polish(fig, axs)
# plt.savefig(osjoin(data_path, f'sweep_{sweep_type}-{v_str}{w_str}{hv_str}{hw_str}.svg'))
# plt.close()
# def plot_score_on_sweep(data_path, net_id_extractor, suffix,
# v=None, w=None, hv=None, hw=None, test_only=True, titles=[],
# normalize_l=None,):
# # figh, figw = 4, 4
# nrows = 2
# if test_only:
# nrows = 1
# # figh, figw = 2, 8
# cbar_labels = ['fast', 'slow']
# v_str, w_str = 'x', 'x'
# hv_str , hw_str = '', ''
# sort_by_heter = False
# pad = -15
# cbar_y = 0.55
# # s-sweep (var=0)
# if (v is not None) and (w is None):
# hv, hw = 0, 0 # overwrites the hv, hw arguments
# v_str = str(v)
# sweep_type = 's'
# cbar_title = r'$\mathbb E[\boldsymbol{\tau_s}]$'
# cmap = plt.cm.Greens_r
# # n-sweep (var=0)
# elif (v is None) and (w is not None):
# hv, hw = 0, 0 # overwrites the hv, hw arguments
# w_str = str(w)
# sweep_type = 'n'
# cbar_title = r'$\mathbb E[\boldsymbol{\tau_n}]$'
# cmap = plt.cm.Blues_r
# # h-sweep
# else:
# pad = -15
# cbar_labels = ['hom.', 'het.']
# v_str = str(v)
# w_str = str(w)
# hv_str , hw_str = 'x', 'x'
# # hn-sweep (var w = constant)
# if (hw is not None) and (hv is None):
# hw_str = str(hw)
# sweep_type = 'hn'
# cbar_title = r'$Var[\boldsymbol{\tau_n}]$'
# cmap = plt.cm.Reds_r
# # hs-sweep (var v = constant)
# elif (hw is None) and (hv is not None):
# hv_str = str(hv)
# sweep_type = 'hs'
# cbar_title = r'$Var[\boldsymbol{\tau_s}]$'
# cmap = plt.cm.Wistia_r
# else:
# print('Warning: You have requested all heterogenous combinations. It will be a mess!')
# sort_by_heter = True
# hv_str = str(hv_str)
# hw_str = str(hw_str)
# sweep_type = 'h'
# cbar_title = r'$Var[\boldsymbol{\tau_s}]$'+'\n + \n'+ r'$Var[\boldsymbol{\tau_n}]$'
# cbar_y = .65
# cmap = plt.cm.viridis
# net_ids = net_id_extractor(v, w, hv, hw)
# for i, net_id in enumerate(net_ids):
# _, df_mean, df_std, dims, ks, ls = utils.get_scores_of_max_k(data_path, net_id, suffix=suffix)
# if normalize_l is not None:
# ls = ls/(1.*normalize_l)
# xlabel = r'$\boldsymbol{\Delta / \tau_0}$'
# else:
# xlabel = r'$\boldsymbol{\Delta}$'
# # setting up the figure
# if i == 0:
# fig, axs = plt.subplots(nrows, dims, figsize=(3.5*dims, 3), sharey=True)
# if type(axs) == np.ndarray:
# axs = np.array(axs).reshape(-1, dims)
# else:
# axs = np.array([axs]).reshape(-1, dims)
# #else:
# for ax in axs.ravel():
# ax.set_prop_cycle(cycler('color', cmap(np.linspace(0,0.7,len(net_ids)))))
# for dim in range(dims):
# axs[0,dim].semilogx(ls, df_mean[f'score_test_dim_{dim+1}'], label=net_id,)
# if not test_only:
# axs[1,dim].semilogx(ls, df_mean[f'score_train_dim_{dim+1}'], label=net_id, )
# if len(net_ids): # plot if there are any nets available
# # colorbar
# for ax in axs[:,-1]:
# cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.1)
# cbar = plt.colorbar(plt.cm.ScalarMappable(norm=Normalize(), cmap=cmap), cax=cax, ticks=[0,1])
# cbar.set_label(cbar_title, labelpad=pad, )# y=cbar_y, labelpad=pad, )#rotation=0,) #
# cbar.ax.set_yticklabels(cbar_labels)
# # ylabel
# for ax in axs[:,0]:
# ax.set_ylabel(r'$\boldsymbol{R^2 (\Delta)}$')
# ax.set_ylim(0,1)
# # xlabel
# for ax in axs[-1,:]:
# ax.set_xlabel(xlabel)
# # titles
# for i, title in enumerate(titles):
# axs[i,0].set_title(title)
# #axs[i,0].set_title(r'$\textbf{Periodic }$')
# #axs[0,1].set_title(r'$\textbf{"Mildly" chaotic}$')
# #axs[0,2].set_title(r'$\textbf{"Strongly" chaotic}$')
# for ax in axs.ravel():
# ax.hlines(0, ls[0], ls[-1], color='k', linestyle='--')
# ax.legend(loc='upper right')
# polish(fig, axs)
# plt.savefig(osjoin(data_path, f'sweep_{sweep_type}-{v_str}{w_str}{hv_str}{hw_str}.svg'))
# plt.close('all')
# def plot_score_forecast(data_path, res_path, net_id, base_stride, feature_size=None, **kwargs):
# """
# The vertical dash line specifies the minimal number of samples necessary to cover a full
# cycle of a purely periodic input with frequency equal to the dominant frequency of the
# input. Simply put, fof training sizes below that level, even a purely periodic singal
# may not be predicted perfectly.
# The horizontal line marks the second peak of autocorrelation peak, which is an optimistic
# proxy of the autocorrelation time (it's actually larger than that), which is defined by
# the first zero-crossing time of the autocorrelation of the input.
# """
# df = pd.read_csv(osjoin(data_path, net_id, 'scores.csv'))#, converters={'pows': pd.eval})
# ks = np.sort(df.k.unique())
# ls = np.sort(df.l.unique())
# stats = df.groupby(['k','l']).agg(['mean', 'std'])
# df_mean, df_std = utils.extract_mean_std_from_groupped(stats)
# # a min-max that is the same for all dims
# # mean_min = df_mean[get_score_cols(df_mean)].min().min()
# # mean_max = df_mean[get_score_cols(df_mean)].max().max()
# mean_min, mean_max = 0., 1.
# std_max = df_std[utils.get_score_cols(df_std)].max().max()
# std_min = df_std[utils.get_score_cols(df_std)].min().min()
# if std_min == 0:
# std_min += 1e-12
# grid_col_mapping = {'train': 0, 'test': 1}
# for dim in range(1, utils.get_dim(df_mean)+1):
# _df_mean = utils.get_dim_subdf(df_mean, dim)
# _df_std = utils.get_dim_subdf(df_std, dim)
# fig = plt.figure(figsize=(5, 5))
# grid = ImageGrid(fig, 111, # as in plt.subplot(111)
# nrows_ncols=(2,2),
# axes_pad=0.15,
# share_all=False,
# cbar_location="right",
# cbar_mode="edge",
# cbar_size="7%",
# cbar_pad=0.15,
# aspect=False
# )
# axs = np.array(grid.axes_row).reshape(2,2)
# for mode in ['train', 'test']:
# grid_col = grid_col_mapping[mode]
# get_cols = eval(f'get_{mode}_cols')
# # extract aggregated values from dataframes
# means = _df_mean[get_cols(_df_mean)].values.reshape(len(ks), len(ls))
# stds = _df_std[get_cols(_df_std)].values.reshape(len(ks), len(ls))
# mm = axs[0,grid_col].pcolormesh(ks, ls/base_stride, means.T, cmap='coolwarm_r',
# # norm=LogNorm(vmin=mean_min+1e-12, vmax=mean_max)
# vmin = mean_min, vmax = mean_max
# )
# ms = axs[1,grid_col].pcolormesh(ks, ls/base_stride, stds.T, cmap='PiYG_r',
# norm=LogNorm(vmin=std_min, vmax=std_max)
# # vmin = std_min, vmax = std_max)
# )
# if grid_col==1:
# axs[0,grid_col].cax.colorbar(mm, label= 'mean of score')
# axs[1,grid_col].cax.colorbar(ms, label= 'std of score')
# for ax in axs.ravel():
# ax.set_yscale('log')
# xlims = ax.get_xlim()
# ax.hlines(1, xlims[0], xlims[1], linestyle='--', color='k')
# if feature_size is not None:
# ylims = ax.get_ylim()
# kmin = base_stride/feature_size
# ax.vlines(kmin, ylims[0], ylims[1], linestyle='--', color='k')
# for ax in axs[1,:]:
# ax.set_xlabel(r'training size ($k/k_0$)')
# for ax in axs[:,0]:
# ax.set_ylabel(r'horizon ($l/l_0$)')
# axs[0,0].set_title('Train')
# axs[0,1].set_title('Test')
# plt.suptitle(net_id+ f' - dim {dim}')
# plt.savefig(osjoin(res_path, f'score_{dim}_{net_id}.png'), dpi=100,
# bbox_inches='tight')
# plt.close()
# def plot_score_capacity(data_path, res_path, net_id, **kwargs):
# df = pd.read_csv(osjoin(data_path, net_id, 'scores.csv'))#, converters={'pows': pd.eval})
# df['degree'] = df['pows'].apply(lambda row: sum(eval(row)))
# stats = df.groupby(['k', 'pows', 'degree']).aggregate(['mean','std'])
# df_mean, df_std = utils.extract_mean_std_from_groupped(stats)
# df_mean = df_mean.reset_index().sort_values('degree')
# df_std = df_std.reset_index().sort_values('degree')
# # a min-max that is the same for all dims
# # mean_min = df_mean[get_score_cols(df_mean)].min().min()
# # mean_max = df_mean[get_score_cols(df_mean)].max().max()
# mean_min, mean_max = 0., 1.
# std_max = df_std[utils.get_score_cols(df_std)].max().max()
# std_min = df_std[utils.get_score_cols(df_std)].min().min()
# if std_min == 0:
# std_min += 1e-12
# d2c = {
# 0: 'Reds',
# 1: 'Greens',
# 2: 'Oranges',
# 3: 'Purples',
# 4: 'Blues',
# 5: 'Greys',
# 6: 'copper_r',
# }
# for dim in range(1, utils.get_dim(df_mean)+1):
# fig, axs = plt.subplots(2,2, figsize=(20,2), sharex=True, sharey=True, layout='constrained')
# for mode_id, mode in enumerate(['train', 'test']):
# _mean = df_mean.pivot(index='k', columns = 'pows', values=f'score_{mode}_dim_{dim}')
# _std = df_std.pivot(index='k', columns = 'pows', values=f'score_{mode}_dim_{dim}')
# _mean.loc['deg',:] = [int(sum(eval(col))) for col in _mean.columns]
# _std.loc['deg',:] = [int(sum(eval(col))) for col in _std.columns]
# # to sort out the columns by the overal degree
# _mean = _mean.loc[:, df_mean[df_mean.k==df_mean.k[0]].pows]
# _std = _std.loc[:, df_std[df_std.k==df_std.k[0]].pows]
# for order in range(int(_mean.iloc[-1,:].max())+1):
# _m = _mean.copy()
# _s = _std.copy()
# _m.loc[:, _m.iloc[-1,:]!=order] = np.nan
# _s.loc[:, _s.iloc[-1,:]!=order] = np.nan
# axs[0, mode_id].imshow(_m.iloc[:-1, :], cmap = d2c[order], origin='lower',
# vmin = mean_min, vmax = mean_max)
# axs[1, mode_id].imshow(_s.iloc[:-1, :], cmap = d2c[order], origin='lower',
# vmin = std_min, vmax = std_max)
# for ax in axs[:,0]:
# ax.set_yticks(range(len(_s.index[:-1])), [str(k) for k in _s.index[:-1]])
# ax.set_ylabel(r'k')
# for ax in axs[-1,:]:
# ax.set_xticks(range(len(_s.columns)), _s.columns, rotation=90)
# ax.set_xlabel(r'Test polynomial')
# axs[0,0].set_title('Train')
# axs[0,1].set_title('Test')
# plt.suptitle(net_id+ f' - dim {dim}')
# plt.savefig(osjoin(res_path, f'score_heatmap_{dim}_{net_id}.png'), dpi=100, bbox_inches='tight')
# plt.close()
# def plot_on_heter_grid(data_dict, res_path, title=None,
# vlims=None, cmap='viridis', norm=None,
# ):
# """
# Given a data dictionary composed of `net_id` as key and data as an array,
# plots the data as a heat-map on the heterogeneity grid.
# """
# net_ids = sorted(list(data_dict.keys()))
# nvar = max([int(net_id[0]) for net_id in net_ids]) - 1
# data, net_labels = utils.group_and_create_arrays(data_dict, nvar)
# # tick labels
# if nvar == 1:
# var_labels = ['homog', 'heter']
# elif nvar == 2:
# var_labels = ['none', 'low', 'high']
# elif nvar == 3:
# var_labels = ['none', 'low', 'med', 'high']
# elif nvar == 4:
# var_labels = ['none', 'low', 'mid-low', 'mid-high', 'high']
# elif nvar == 5:
# var_labels = ['none', 'low', 'mid-low', 'mid', 'mid-high', 'high']
# # colors range
# if vlims is None:
# vals = np.concatenate([v.ravel() for v in data.values()])
# vlims = [min(vals), max(vals)]
# # color normalization
# if norm is None:
# norm = Normalize
# cmap_ticks = np.linspace(vlims[0], vlims[1], 5)
# elif norm == 'log':
# norm = LogNorm
# cmap_ticks = np.logspace(np.log10(vlims[0]), np.log10(vlims[1]), 5)
# else:
# raise ValueError(f'Specified color norm {norm} not recognized.')
# fig, axs = plt.subplots(nvar + 2, nvar + 2,
# sharex='col', sharey='row', figsize=(6.5, 6),
# height_ratios=[1]+nvar*[nvar+1]+[1],
# width_ratios=[1]+nvar*[nvar+1]+[1],
# )
# for w_id in range(nvar + 2):
# for v_id in range(nvar + 2):
# ax = axs[nvar+1- w_id, v_id] # ax selection
# im = ax.pcolormesh(data[f'{v_id}{w_id}'],
# norm = norm(vlims[0], vlims[1]),
# cmap = cmap)
# # homog. labels
# ax.annotate(f'{v_id}{w_id}', (0.5, 0.5), transform=ax.transAxes,
# ha='center', va='center', fontsize=.65*MFS*3/nvar, color='white',
# path_effects=[patheffects.withStroke(foreground='black', capstyle="round", linewidth=.5,)])
# # heterog. labels
# if net_labels is not None:
# labels = net_labels[f'{v_id}{w_id}'] # heter labels
# for wh_id in range(labels.shape[0]):
# for vh_id in range(labels.shape[1]):
# if (wh_id or vh_id):
# ax.annotate(labels[wh_id, vh_id], (0.5+vh_id, .5+wh_id),
# transform=ax.transAxes, ha='center', va='center',
# fontsize=.45*MFS*3/nvar, color='lightgray',
# path_effects=[patheffects.withStroke(foreground='black',
# capstyle="round",
# linewidth=.5,)])
# # tick labels
# if v_id == 0:
# nw = nvar + 1
# if w_id==0 or w_id==nvar+1:
# nw = 1
# ax.set_yticks(np.arange(nw)+0.5, var_labels[:nw])
# if w_id == 0:
# nv = nvar + 1
# if v_id==0 or v_id==nvar+1:
# nv = 1
# ax.set_xticks(np.arange(nv)+0.5, var_labels[:nv], rotation = 60)
# polish(fig, axs, 1)
# # the following must be done after tight_layout otherwise, things become dirty!
# fig.text(0.5, -0.01, r'$Var[\boldsymbol{\tau_n}]$', ha='center', fontsize=.85*MFS, color='k')
# fig.text(-0.01, 0.5, r'$Var[\boldsymbol{\tau_s}]$', ha='center', fontsize=.85*MFS, color='k', rotation=90)
# cbar = fig.colorbar(im, ax=axs.ravel().tolist(), use_gridspec=True, pad=0.03)
# # ticks format is different for scores and horizon index
# if cmap_ticks[-1]>1:
# cbar.ax.set_yticks(cmap_ticks, ['{:.2e}'.format(tick) for tick in cmap_ticks]);
# else:
# cbar.ax.set_yticks(cmap_ticks, ['{:.2}'.format(tick) for tick in cmap_ticks]);
# if title:
# fig.suptitle(title.replace('_',' '), y=1.02, fontsize=MFS)
# plt.savefig(osjoin(res_path, f'{title}.svg'))
# plt.close('all')
# def plot_on_half_grid(data_dict, res_path, title=None,
# vlims=None, cmap='viridis', norm=None,
# ):
# max_mu = max([eval(net_id[0]) for net_id in data_dict.keys()])
# max_std= max([eval(net_id[2]) for net_id in data_dict.keys()])
# # nvar = max_mu - 1
# scores = np.full((max_mu+1, max_std+1), np.nan)
# if vlims is None:
# vlims = [min(data_dict.values()), max(data_dict.values())]
# # color normalization
# if norm is None:
# norm = Normalize
# cmap_ticks = np.linspace(vlims[0], vlims[1], 5)
# elif norm == 'log':
# norm = LogNorm
# cmap_ticks = np.logspace(np.log10(vlims[0]), np.log10(vlims[1]), 5)
# else:
# raise ValueError(f'Specified color norm {norm} not recognized.')
# for net_id, data in data_dict.items():
# i, j = eval(net_id[0]), eval(net_id[2])
# scores[i,j] = data
# fig, ax = plt.subplots(1,1, figsize=(3.5, 3))
# im = ax.pcolormesh(scores.T, norm = norm(vlims[0], vlims[1]), cmap = cmap)
# polish(fig, [ax], 1.5)
# # the following must be done after tight_layout otherwise, things become dirty!
# fig.text(0.5, -0.01, r'$E[\boldsymbol{\tau_n}]$', ha='center', fontsize=.85*MFS, color='k')
# fig.text(-0.01, 0.5, r'$Var[\boldsymbol{\tau_n}]$', ha='center', fontsize=.85*MFS, color='k', rotation=90)
# cbar = fig.colorbar(im, ax=ax, use_gridspec=True, pad=0.03)
# # ticks format is different for scores and horizon index
# if cmap_ticks[-1]>1:
# cbar.ax.set_yticks(cmap_ticks, ['{:.2e}'.format(tick) for tick in cmap_ticks]);
# else:
# cbar.ax.set_yticks(cmap_ticks, ['{:.2}'.format(tick) for tick in cmap_ticks]);
# if title:
# fig.suptitle(title.replace('_',' '), y=1.02, fontsize=MFS)
# plt.savefig(osjoin(res_path, f'{title}.svg'))
# plt.close('all')
# def plot_ipc_k(data_path, res_path, net_id, max_capacity, test=True, **kwargs):
# ipc = pd.read_csv(osjoin(data_path, net_id, 'scores_ipc.csv'))
# ipc['degree'] = ipc['mono_degs'].apply(lambda row: sum(eval(row)))
# ks = sorted(ipc.k.unique())
# stats = ipc.groupby(['k', 'mono_degs', 'degree']).aggregate(['mean','std'])
# df_mean, df_std = utils.extract_mean_std_from_groupped(stats)
# df_mean = df_mean.reset_index().sort_values(['degree', 'k'])
# df_std = df_std.reset_index().sort_values(['degree', 'k'])
# fig, axs = plt.subplots(2,1,figsize=(9,7), sharex=True)
# if test:
# col_name = 'score_test_dim_1'
# else:
# col_name = 'score_train_dim_1'
# df_mean.groupby(['k', 'degree']).agg('sum').reset_index().\
# pivot(index='k',columns='degree', values=col_name).\
# plot(kind='bar', stacked=True, ax=axs[0])
# df_std.groupby(['k', 'degree']).agg('sum').reset_index().\
# pivot(index='k',columns='degree', values=col_name).\
# plot(kind='bar', stacked=True, ax=axs[1])
# for ax in axs[::-1]:
# ax.hlines(max_capacity, -1,len(ks), linestyle='--', color='k')
# ylims = ax.get_ylim()
# ax.set_ylim(0, max(ylims[1], max_capacity*1.3))
# axs[0].set_ylabel('Mean norm. IPC')
# axs[1].set_ylabel('Std. of norm. IPC');
# plt.savefig(osjoin(res_path, f'IPC_{net_id}.svg'))
# plt.close('all')
[docs]
def plot_score_map(res_path):
root, sim = os.path.split(res_path)
df = utils.aggregate_scores(root, sim,)
score_names = set(df.scr)
# set_trace()
hns = sorted(set(df.hn))
degs = sorted(set(df.degree))
deltas = sorted(set(df.delta))
ns = sorted(set(df.n))
u = np.load(osjoin(res_path,'u.npy')).astype(np.float64)
dims = len(u)
for scr in score_names:
for dim in range(dims):
fig, axs = plt.subplots(len(degs), len(ns),
sharex=True, sharey=True,
figsize=(15/5*len(ns), 7/3*len(degs)))
axs = axs.reshape(len(degs), len(ns))
lags = correlation_lags(len(u[0]), len(u[0]),)# mode='same' )
idxs = slice(np.where(lags == min(deltas))[0][0], np.where(lags == max(deltas))[0][0])
lags = lags[idxs]
for row, deg in enumerate(degs):
corr = correlate(u[dim]**deg, u[dim]**deg,)# mode='same' )
corr /= np.max(corr)
corr = abs(corr[idxs])
for col, n in enumerate(ns):
ax = axs[row, col]
sub_df = df[(df.n == n) * (df.degree == deg) * (df.scr == scr)]
scores = np.ones((len(deltas), len(hns)))*np.nan
for hn_idx, hn in enumerate(sorted(set(sub_df.hn))):
scores[:, hn_idx] = sub_df[sub_df.hn==hn][f'score_test_dim_{dim+1}'].values
g = ax.pcolormesh(deltas, hns, scores.T,
vmin=0, vmax=1, cmap='Blues')
ax2 = ax.twinx() # instantiate a second Axes that shares the same x-axis
ax2.plot(lags, corr, lw=1, color='r', ls='--')
ax2.plot(lags, np.zeros_like(lags), color='r', ls=':', lw=.5)
ax2.tick_params(axis='y', labelcolor='r')
if n == max(ns):
ax2.set_ylabel(r'$|C_{uu}(\Delta)|$', color='r', fontsize=14) # we already handled the x-label with ax1
else:
ax2.set_yticklabels([])
for ax in axs[:,0]:
ax.set_ylabel('Hetero level', fontsize=14)
ax.set_yticks(hns, ['None', 'Low', 'Medium', 'High', 'v-High', 'vv-High'][:len(hns)])
for ax in axs[-1, :]:
ax.set_xticks([deltas[0], 0, deltas[-1]],
[r'-2$\tau_0$', '0', r'2$\tau_0$'])
ax.set_xlabel(r'$\Delta$')
for n_idx, n in enumerate(ns):
# if n==max(ns):
# n_ = r'$\infty$'
# else:
# n_ = n
# axs[0, n_idx].set_title(rf'E[$\tau]$={n_}', fontsize=16)
n_ = TAU_MAP[str(int(n))]
axs[0, n_idx].set_title(rf'E[$\tau]$={n_}', fontsize=16)
cb = fig.colorbar(g, ax=axs, location='top', orientation='horizontal',
aspect=40, pad=.1).set_label(label='(exp-transformed) score',
size=20,weight='bold')
polish(fig, axs, tight=False)
plt.savefig(osjoin(res_path, f'scoremap_dim{dim+1}_{scr}.svg'))
plt.close('all')
# def plot_deviation(net_ids, data_path, t_lyap = None, xlims=None, ylims=None,):
# fig, axs = plt.subplots(5,1, figsize=(12,6), sharex=True, sharey=False)
# for net_id in net_ids:
# try:
# y1 = np.load(osjoin(data_path, net_id, 'y_autonom.npy'))#.astype(float)
# except FileNotFoundError:
# print(f'No autonomous dynamic is simulated. Halting autonom plot.')
# break
# y0 = np.load(osjoin(data_path, net_id, 'y.npy'))#.astype(float)
# t = np.load(osjoin(data_path, 't.npy'))#.astype(float)
# N = len(y0)
# L_tst = y1.shape[1]
# L_trn = y0.shape[1] - L_tst
# dy = y0[:,L_trn:] - y1
# # dy/= y0[:,-task.L_tst:] # normalizing
# col = eval(net_id[0])
# axs[col].plot(t[L_trn:], np.linalg.norm(dy, axis=0)/N,
# # np.linalg.norm(dy, axis=0)/np.linalg.norm(y0[:,task.L_trn:], axis=0),
# label=NET_ID_MAP[net_id],
# # label=net_id,
# lw=1)
# # print(net_id, np.linalg.norm(dy[:,0]), np.linalg.norm(dy[:,1]))
# for ax in axs:
# ax.set_xlabel('time [s]')
# ax.set_ylabel('mean deviation norm')
# if xlims is not None:
# ax.set_xlim(*xlims)
# else:
# ax.set_xlim(t[L_trn], t[-1])
# if ylims is not None:
# ax.set_ylim(*ylims)
# if t_lyap is not None:
# ax.vlines(t[L_trn] + t_lyap, *ax.get_ylim(),
# ls='--', color='k', label='Lyapunov time')
# ax.set_xlim(t[L_trn], t[L_trn] + 2*t_lyap)
# ax.legend(ncols=4, loc='upper right')
# fig.suptitle('Autonomous dynamic deviation')
# polish(fig, [axs], tight=True)
# fig.savefig(osjoin(data_path, f'deviation.svg'))
# plt.close('all')
[docs]
def plot_score_score(res_path, rlz_id=None):
"""Plots, for each mean-timescale level ``n`` (first digit of the net_id),
the trial-averaged score of every task of a heterogeneous network
(``net_id`` with a non-zero second digit ``h``) against the score of the
exact same task for the homogeneous network of the same group (the
``net_id`` with ``h=0``). Each ``h`` level is colored by ``hetero_colors``.
The score of the homogeneous network is on the x-axis and the score of the
heterogeneous network on the y-axis, so points above the identity line
indicate that heterogeneity improves performance.
The trn_/tst_ score files of every network are read, and the per-task score
is obtained by averaging across all trials. When ``rlz_id`` is given, only
the score files of that realization (e.g. ``LI-rlz000``) are read. One
figure per score file (state ``LI-rlz000`` and readout such as ``possig1``)
is saved into ``res_path``, named ``score_score_{trn|tst}{scr_name}.png``
to keep the ``rlzXXX`` naming convention alive.
Parameters
----------
res_path : str
The path to a simulation directory, containing one subdirectory per
network ({n}{h} with h=0 the homogeneous network) with the trn_/tst_
score files.
rlz_id : int, optional
The realization to plot. By default all realizations found are plotted.
"""
net_ids = sorted([f.name for f in os.scandir(res_path) if f.is_dir()])
# derive the score-file suffixes (e.g. "_LI-rlz000-possig1") from the
# first network that has any tst_ file, and restrict to the given rlz
scr_names = set()
for net_id in net_ids:
for f in os.scandir(osjoin(res_path, net_id)):
if 'tst_' in f.name:
scr_names.add(f.name.split('.csv')[0].split('tst')[1])
scr_names = sorted(scr_names)
if rlz_id is not None:
scr_names = [s for s in scr_names if f'-rlz{rlz_id:03d}' in s]
if len(scr_names) == 0:
print(f'No score files found under {res_path}.')
return
for scr_name in scr_names:
for set_name in ['trn', 'tst']:
tasks = []
for net_id in net_ids:
path = osjoin(res_path, net_id, f'{set_name}{scr_name}.csv')
if not os.path.exists(path):
continue
df = pd.read_csv(path)
# average over the trial columns (0, 1, ..., n_trials-1)
df = df.set_index(['dim', 'delta', 'deg']).mean(axis=1)
df = df.rename('score').reset_index()
df['net_id'] = net_id
df['n'] = net_id[0]
df['h'] = net_id[1]
tasks.append(df)
if len(tasks) == 0:
continue
df = pd.concat(tasks)
n_group = sorted(set(df['n']))
if len(n_group) == 0:
continue
fig, axs = plt.subplots(
1, len(n_group), sharex=True, sharey=True,
figsize=(3 * len(n_group), 3), layout='tight')
if len(n_group) == 1:
axs = np.array([axs])
for ax, n in zip(axs, n_group):
hom = df[(df.n == n) & (df.h == '0')]
hom = hom.set_index(['dim', 'delta', 'deg'])['score']
if len(hom) == 0:
print(f'No homogeneous (h=0) network found for n={n}.')
continue
for h in sorted(set(df[(df.n == n) & ~df.h.isna()].h)):
het = df[(df.n == n) & (df.h == h)]
het = het.set_index(['dim', 'delta', 'deg'])['score']
idx = hom.index.intersection(het.index)
x = hom.loc[idx].values
y = hom.loc[idx].values if h == '0' else het.loc[idx].values
ax.scatter(x, y, s=8, c=hetero_colors[int(h)],
label=HETERO_LABEL_MAP[str(h)],
alpha=0.75, edgecolors='none')
ax.plot([0, 1], [0, 1], ls='--', lw=1, color='k', zorder=0)
ax.set_title(r'$\mathbb{E}[\tau_s]$ = ' + TAU_MAP[str(n)])
ax.set_xlim(-0.15, 1.05)
ax.set_ylim(-0.15, 1.05)
ax.set_aspect('equal')
ax.legend(loc='lower right', fontsize=8,
title='Heterog. level')
axs[0].set_ylabel(f'heterogeneous score')
for ax in axs:
ax.set_xlabel('homogeneous score')
polish(fig, axs, tight=True)
fig.savefig(osjoin(res_path, f'score_score_{set_name}{scr_name}.png'))
plt.close('all')
[docs]
def plot_angle_vs_score(res_path, rlz_id=None):
"""Plots, for every task, the minimal angle between the network state and
the task against the task's benchmark score, for all networks.
Networks are split into groups by their mean-timescale level ``n`` (the
first digit of the ``net_id``, i.e. the mean tau, per ``TAU_MAP``); each
group gets one panel. Within each panel the different heterogeneity
levels ``h`` (the second digit of the ``net_id``) are plotted with
different colors, using the same color scheme as the other plots of this
module (``hetero_colors``, labeled by ``HETERO_LABEL_MAP``). The top and
bottom rows use the train and test splits respectively, so every task
appears once per split.
The task angle is read from the ``angle_{state}.csv`` files produced by
:meth:`hetero.sim.Simulation.compute_minimal_angles` and the score from
the paired ``trn_/tst_{state}-{readout}.csv`` files, averaged across
trials. When ``rlz_id`` is given, only that realization (e.g. ``0`` ->
``LI-rlz000``) is plotted.
One figure per score file is saved into ``res_path``, each coloring the
scatter with a different scheme: by the heterogeneity/color level
``h`` (the default, ``angle_vs_score{scr_name}.png``), and by the task
properties ``dim``, ``delta`` and ``deg`` (``angle_vs_score{scr_name}-
{dim,delta,deg}.png``). The task-property figures reuse the colormaps of
:func:`plot_udf_vs_complexity`: ``dim`` and ``deg`` are categorical
(``Set2`` and ``Accent``), ``delta`` is continuous (``coolwarm``).
Parameters
----------
res_path : str
The path to a simulation directory, containing one subdirectory per
network ({n}{h}) with the angle and trn_/tst_ score files of each
benchmarked realization.
rlz_id : int, optional
The realization to plot. By default all realizations found are
plotted.
"""
net_ids = sorted([f.name for f in os.scandir(res_path) if f.is_dir()])
# angle files: (net_id, state token) -> path; the token is e.g. '_LI-rlz000'
ang_files = {}
for net_id in net_ids:
for f in os.scandir(osjoin(res_path, net_id)):
if f.name.startswith('angle_') and f.name.endswith('.csv'):
tok = f.name[len('angle'):-len('.csv')]
ang_files[(net_id, tok)] = f.path
# score-file suffixes (e.g. '_LI-rlz000-possig1') found on any tst_ file
scr_names = set()
for net_id in net_ids:
for f in os.scandir(osjoin(res_path, net_id)):
if 'tst_' in f.name:
scr_names.add(f.name.split('.csv')[0].split('tst')[1])
scr_names = sorted(scr_names)
if rlz_id is not None:
scr_names = [s for s in scr_names if f'-rlz{rlz_id:03d}' in s]
if len(scr_names) == 0:
print(f'No score files found under {res_path}.')
return
n_levels = sorted(set(int(net_id[0]) for net_id in net_ids))
h_levels = sorted(set(int(net_id[1]) for net_id in net_ids))
for scr_name in scr_names:
# the state token shared by this score file and an angle file
tok = next((t for (_, t) in ang_files if scr_name.startswith(t)), None)
if tok is None:
print(f'No angle files matching {scr_name}; skipping.')
continue
# one figure per coloring scheme: default 'h' (per-network colors +
# discrete colorbar on the heterogeneity levels) and the task
# properties dim / delta / deg (per-task colors + colorbar, same
# colormaps as plot_udf_vs_complexity).
schemes = [('h', 'h', 'hmap', 'heterogeneity level'),
('dim', 'cat', 'Set2', 'input dim'),
('delta', 'cont', 'coolwarm', r'$\Delta$'),
('deg', 'cat', 'Accent', 'degree')]
for prop, kind, cmap_name, prop_label in schemes:
# categorical cmap/norm for dim & deg (one band per distinct value)
cmap = norm = None
uniq = None
if kind == 'cat':
ref = next((ang_files[(net_id, tok)] for net_id in net_ids
if (net_id, tok) in ang_files), None)
if ref is None:
continue
uniq = np.unique(pd.read_csv(ref)[prop].to_numpy())
base = plt.get_cmap(cmap_name).colors
cmap = ListedColormap([base[k % len(base)]
for k in range(len(uniq))])
norm = BoundaryNorm([uniq[0] - 0.5]
+ list((uniq[:-1] + uniq[1:]) / 2)
+ [uniq[-1] + 0.5],
ncolors=len(uniq), clip=True)
elif kind == 'cont':
cmap = cmap_name
elif kind == 'h':
# discrete heterogeneity-level colorbar using hetero_colors,
# mirroring the categorical (dim/deg) construction
uniq = np.array(h_levels)
cmap = ListedColormap(
[hetero_colors[int(h) % len(hetero_colors)] for h in h_levels])
norm = BoundaryNorm([uniq[0] - 0.5]
+ list((uniq[:-1] + uniq[1:]) / 2)
+ [uniq[-1] + 0.5],
ncolors=len(uniq), clip=True)
fig, axs = plt.subplots(
2, len(n_levels), sharex=True, sharey='row',
figsize=(1.8 * len(n_levels), 3.), squeeze=False)
if len(n_levels) == 1:
fig.subplots_adjust(left=0.12, right=0.99, top=0.92,
bottom=0.12)
sc = None # last scatter (for the colorbar)
for row, set_name in enumerate(['trn', 'tst']):
for col, n in enumerate(n_levels):
ax = axs[row, col]
for net_id in net_ids:
if int(net_id[0]) != n:
continue
ang_path = ang_files.get((net_id, tok))
scr_path = osjoin(res_path, net_id, f'{set_name}{scr_name}.csv')
if ang_path is None or not os.path.exists(scr_path):
continue
ang = pd.read_csv(ang_path)
scr = pd.read_csv(scr_path)
# align tasks by the template columns, average over trials
scr = scr.set_index(['dim', 'delta', 'deg']).mean(axis=1)
scr = scr.reindex(ang.set_index(['dim', 'delta', 'deg']).index)
if prop == 'h':
cc = np.full(len(ang), int(net_id[1])) # heterogeneity level
else:
cc = ang[prop].to_numpy()
sc = ax.scatter(ang[f'angle_{set_name}'],
scr.to_numpy(), s=4, c=cc, cmap=cmap,
norm=norm, alpha=0.75,
edgecolors='none')
ax.axhline(0.0, ls=':', lw=1, color='gray', zorder=0)
ax.set_xlim(0, np.pi / 2)
ax.set_xticks([0, np.pi/4, np.pi/2])
ax.set_xticklabels([r'$0$', r'$\pi/4$', r'$\pi/2$'])
if row==0:
ax.set_title(r'$\mathbb{E}[\tau_s]$ = ' + TAU_MAP[str(n)])
# one colorbar hugging this row, placed *after* the panels so the
# layout math sees it
if prop == 'h':
cb = fig.colorbar(sc, ax=axs[row, :],
label='Heterogeneity level',
ticks=h_levels, pad=0.02)
cb.ax.set_yticklabels(
[f'{HETERO_LABEL_MAP[str(int(h))][:4]}' for h in h_levels],
rotation=90, fontsize=6, va='center')
elif kind == 'cat':
fig.colorbar(sc, ax=axs[row, :], label=prop_label,
ticks=uniq, pad=0.02,
format='%d' if np.all(uniq == np.round(uniq))
else None)
else:
fig.colorbar(sc, ax=axs[row, :], label=prop_label,
pad=0.02)
axs[0, 0].set_ylabel('train score')
axs[1, 0].set_ylabel('test score')
for ax in axs[-1]:
ax.set_xlabel('min angle (rad)')
polish(fig, axs, tight=False)
suffix = '' if prop == 'h' else f'-{prop}'
fig.savefig(osjoin(res_path, f'angles{scr_name}{suffix}.png'))
plt.close('all')
[docs]
def plot_overlap_cdf(res_path, rlz_id=None):
"""Plots the cumulative distribution of the minimal angles of every task,
grouped into panels by the mean-timescale level ``n`` (the first digit of
the ``net_id``, i.e. the mean tau, per ``TAU_MAP``). Within each panel the
different heterogeneity levels ``h`` (the second digit of the ``net_id``)
are drawn in the module's ``hetero_colors``, labeled by
``HETERO_LABEL_MAP``.
The angles are read from the ``angle_{state}.csv`` files produced by
:meth:`hetero.sim.Simulation.compute_minimal_angles`. Each panel pools the
train and test angles of every subtask of every network sharing that ``n``,
held separately for each ``h``, and draws the fraction of tasks whose angle
lies below each value (an empirical CDF): a curve shifting toward zero
means the tasks are, on average, better aligned with the network state. One
figure per angle file (state token) is saved into ``res_path`` as
``overlap_cdf{state}.png``. When ``rlz_id`` is given, only that realization
(e.g. ``0`` -> ``-rlz000``) is plotted.
Parameters
----------
res_path : str
The path to a simulation directory, containing one subdirectory per
network ({n}{h}) with the angle files of each benchmarked realization.
rlz_id : int, optional
The realization to plot. By default all realizations found are
plotted.
"""
net_ids = sorted([f.name for f in os.scandir(res_path) if f.is_dir()])
# angle files: (net_id, state token) -> path; the token is e.g. '_LI-rlz000'
ang_files = {}
for net_id in net_ids:
for f in os.scandir(osjoin(res_path, net_id)):
if f.name.startswith('angle_') and f.name.endswith('.csv'):
tok = f.name[len('angle'):-len('.csv')]
ang_files[(net_id, tok)] = f.path
if len(ang_files) == 0:
print(f'No angle files found under {res_path}.')
return
tokens = sorted({t for _, t in ang_files})
if rlz_id is not None:
tokens = [t for t in tokens if f'-rlz{rlz_id:03d}' in t]
if len(tokens) == 0:
print(f'No angle files for realization {rlz_id} under {res_path}.')
return
n_levels = sorted(set(int(net_id[0]) for net_id in net_ids))
h_levels = sorted(set(int(net_id[1]) for net_id in net_ids))
for tok in tokens:
ang_tok = {(n_id, t): p for (n_id, t), p in ang_files.items()
if t == tok}
fig, axs = plt.subplots(1, len(n_levels), sharex=True, sharey=True,
squeeze=False,
figsize=(2.6 * len(n_levels), 2))
for col, n in enumerate(n_levels):
ax = axs[0, col]
for h in h_levels:
pooled = []
for net_id in net_ids:
if int(net_id[0]) != n or int(net_id[1]) != h:
continue
p = ang_tok.get((net_id, tok))
if p is None:
continue
ang = pd.read_csv(p)[['angle_trn', 'angle_tst']]
# ang = np.cos(np.pi/2 - ang) # the overlap to state is important
ang = np.cos(ang) # the overlap to state is important
pooled.append(ang.to_numpy().ravel())
if len(pooled) == 0:
continue
pooled = np.concatenate(pooled)
pooled = pooled[~np.isnan(pooled)]
xs = np.sort(pooled)
ys = np.arange(1, len(xs) + 1) / len(xs)
ax.step(xs, ys, where='post', lw=1.6,
color=hetero_colors[h],
label=HETERO_LABEL_MAP[str(h)] if col == 0 else None)
ax.set_xlim(0, 1.02)
ax.set_ylim(0, 1.02)
# ax.set_xticks([0, np.pi/4, np.pi/2])
# ax.set_xticklabels([r'$0$', r'$\pi/4$', r'$\pi/2$'])
ax.set_title(r'$\mathbb{E}[\tau_s]$ = ' + TAU_MAP[str(n)])
axs[0][0].set_ylabel('cumulative fraction of tasks')
for ax in axs[-1]:
ax.set_xlabel(r'Overlap (cos($\textbf{X}, y_k)$)')
handles = [mpl.lines.Line2D([], [], color=hetero_colors[h], lw=2,
label=HETERO_LABEL_MAP[str(h)])
for h in h_levels]
axs[0][-1].legend(handles=handles, loc='lower right', fontsize=7)
polish(fig, axs, tight=True)
fig.savefig(osjoin(res_path, f'overlap_cdf{tok}.png'))
plt.close('all')
[docs]
def plot_svals_erank(res_path, rlz_id=None):
"""Plots the singular-value spectrum of every network's activity and the
cumulative distribution of its effective rank, grouped into panels by the
mean-timescale level ``n`` (the first digit of the ``net_id``, i.e. the
mean tau, per ``TAU_MAP``), with the heterogeneity levels ``h`` (the second
digit of the ``net_id``) drawn in the module's ``hetero_colors`` and labeled
by ``HETERO_LABEL_MAP`` -- the exact grouping of
:func:`plot_overlap_cdf`.
The spectra and effective ranks are read from the ``svals_{state}.csv`` and
``erank_{state}.csv`` files produced by
:meth:`hetero.sim.Simulation.compute_minimal_angles`. Two figures per state
token are saved into ``res_path``:
- ``svals{tok}.png`` -- the mean (solid line) +/- one standard deviation
(shaded band) of the singular-value spectra pooled over the train and
test slices and over every network sharing that ``n`` and ``h``. Each
spectrum is normalized by its own largest singular value so the curves
isolate the *shape* of the spectrum (how many comparable modes the
activity has) rather than the overall activity scale; ``x`` runs over
the modes in descending order, rescaled to the unit fraction ``[0, 1]``.
A curve hugging 1 over a long stretch means the activity is effectively
low-dimensional (few dominant modes); a fast decay to 0 means the
variance is spread over many modes.
- ``erank{tok}.png`` -- the empirical CDF of the effective rank (the Roy &
Vetterli spectral-entropy rank, see
:meth:`hetero.sim.Simulation._effective_rank`) pooled over the train and
test slices and over every network of each (``n``, ``h``) pair, i.e. the
fraction of states whose effective rank lies below each value.
When ``rlz_id`` is given, only that realization (e.g. ``0`` ->
``-rlz000``) is plotted.
Parameters
----------
res_path : str
The path to a simulation directory, containing one subdirectory per
network ({n}{h}) with the svals/erank files of each benchmarked
realization.
rlz_id : int, optional
The realization to plot. By default all realizations found are
plotted.
"""
net_ids = sorted([f.name for f in os.scandir(res_path) if f.is_dir()])
# (net_id, state token) -> path; the token is e.g. '_LI-rlz000'
svals_files, erank_files = {}, {}
for net_id in net_ids:
for f in os.scandir(osjoin(res_path, net_id)):
if not f.name.endswith('.csv'):
continue
if f.name.startswith('svals_'):
svals_files[(net_id, f.name[len('svals'):-len('.csv')])] = f.path
elif f.name.startswith('erank_'):
erank_files[(net_id, f.name[len('erank'):-len('.csv')])] = f.path
if len(svals_files) == 0 and len(erank_files) == 0:
print(f'No svals/erank files found under {res_path}.')
return
tokens = sorted({t for _, t in svals_files} | {t for _, t in erank_files})
if rlz_id is not None:
tokens = [t for t in tokens if f'-rlz{rlz_id:03d}' in t]
if len(tokens) == 0:
print(f'No svals/erank files for realization {rlz_id} under {res_path}.')
return
n_levels = sorted(set(int(net_id[0]) for net_id in net_ids))
h_levels = sorted(set(int(net_id[1]) for net_id in net_ids))
handles = [mpl.lines.Line2D([], [], color=hetero_colors[h], lw=2,
label=HETERO_LABEL_MAP[str(h)])
for h in h_levels]
# common x-grid for the spectra: mode index rescaled to the unit fraction
grid = np.linspace(0.0, 1.0, 200)
for tok in tokens:
# ------------------------------------------------------------------
# singular-value spectra: mean +/- std, each spectrum normalized by
# its own largest singular value (shape, not scale)
# ------------------------------------------------------------------
fig, axs = plt.subplots(1, len(n_levels), sharex=True, sharey=True,
squeeze=False,
figsize=(2.6 * len(n_levels), 2.2))
for col, n in enumerate(n_levels):
ax = axs[0, col]
for h in h_levels:
spectra = []
for net_id in net_ids:
if int(net_id[0]) != n or int(net_id[1]) != h:
continue
p = svals_files.get((net_id, tok))
if p is None:
continue
sv = pd.read_csv(p)
for c in ('s_trn', 's_tst'):
s = sv[c].to_numpy()
s = s[~np.isnan(s)]
if s.size == 0 or not np.any(s > 0):
continue
s = s / s.max()
xs = np.linspace(0.0, 1.0, s.size)
spectra.append(np.interp(grid, xs, s))
if len(spectra) == 0:
continue
spectra = np.vstack(spectra)
mu, sd = spectra.mean(axis=0), spectra.std(axis=0)
cc = hetero_colors[h]
ax.plot(grid, mu, color=cc, lw=1.6,
label=HETERO_LABEL_MAP[str(h)] if col == 0 else None)
ax.fill_between(grid, mu - sd, mu + sd, color=cc,
alpha=0.25, lw=0)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1.02)
ax.set_title(r'$\mathbb{E}[\tau_s]$ = ' + TAU_MAP[str(n)])
axs[0][0].set_ylabel(r'singular value ($s/s_{max}$)')
for ax in axs[-1]:
ax.set_xlabel('fraction of modes')
axs[0][-1].legend(handles=handles, loc='upper right', fontsize=7)
polish(fig, axs, tight=True)
fig.savefig(osjoin(res_path, f'svals{tok}.png'))
plt.close('all')
# ------------------------------------------------------------------
# effective rank: empirical CDF pooled over slices and networks
# ------------------------------------------------------------------
fig, axs = plt.subplots(1, len(n_levels), sharex=True, sharey=True,
squeeze=False,
figsize=(2.6 * len(n_levels), 2))
for col, n in enumerate(n_levels):
ax = axs[0, col]
for h in h_levels:
pooled = []
for net_id in net_ids:
if int(net_id[0]) != n or int(net_id[1]) != h:
continue
p = erank_files.get((net_id, tok))
if p is None:
continue
er = pd.read_csv(p)['erank'].to_numpy()
pooled.append(er[~np.isnan(er)])
if len(pooled) == 0:
continue
pooled = np.concatenate(pooled)
if pooled.size == 0:
continue
xs = np.sort(pooled)
ys = np.arange(1, len(xs) + 1) / len(xs)
ax.step(xs, ys, where='post', lw=1.6, color=hetero_colors[h],
label=HETERO_LABEL_MAP[str(h)] if col == 0 else None)
ax.set_xlim(0, None)
ax.set_ylim(0, 1.02)
ax.set_title(r'$\mathbb{E}[\tau_s]$ = ' + TAU_MAP[str(n)])
axs[0][0].set_ylabel('cumulative fraction of states')
for ax in axs[-1]:
ax.set_xlabel('effective rank of the state')
axs[0][-1].legend(handles=handles, loc='lower right', fontsize=7)
polish(fig, axs, tight=True)
fig.savefig(osjoin(res_path, f'erank{tok}.png'))
plt.close('all')
# @utils.timer
# def plot_svd_grid(net_ids, res_path, N, suffix='', conv_params=None, activation=None):
# # nv, nw = utils.get_ntau_from_net_ids(net_ids, True, True)
# # cmaps = sorted([cmap for cmap in plt.cm._colormaps() if cmap[-1]=='s'])
# # cmaps.remove('Greys') # Grays = Greys
# # cmaps = [eval('plt.cm.'+cmap) for cmap in cmaps]
# # cmaps = [cmap(np.linspace(0.2, 1 , nv+2)) for cmap in cmaps]
# # fig, axs = plt.subplots(nw+2, nv+2, figsize=(4.*(nv+2), 2.*(nw+2)),
# # gridspec_kw={'width_ratios':[0.25, 1, 1, 1, 0.25]}
# # )
# # axs = axs.reshape(nw+2, nv+2)
# nv, nw, nhv, nhw = utils.get_net_levels_from_net_ids(net_ids)
# fig, axs = plt.subplots(nw, nv, figsize=(3.5*nv, 3.5*nw), sharey=True,
# # gridspec_kw={'width_ratios':[0.25] + (nv-2)*[1] + [0.25]}
# )
# axs = axs.reshape(nw, nv)
# blues = mpl.colormaps['Blues'].resampled(nhv)
# newcolors = blues(np.linspace(0, 1, nhv))
# newcolors[0] = np.array([1, 0, 0., 1])
# for ax in axs.ravel():
# ax.plot([1, N], [1, 1/N], ':', color='C1', label=r'$1/n$')
# for net_id in net_ids:
# col, row = utils.get_net_ordinal_idx(net_id, net_ids)
# # color_id = int(net_id[2])
# # cmap_id = int(net_id[3])
# ax = axs[row, col]
# # cmap = cmaps[cmap_id % len(cmaps)]
# # color = cmap[color_id]
# try:
# if suffix != '':
# print(conv_params, suffix)
# _, _, _X = spike_loader(data_path=res_path, net_id=net_id, spikes_name=f'spikes{suffix}',
# convolve = True, conv_params=conv_params,
# add_bias = True, N_neurons=N, )
# else:
# _, _, _X = ds_loader(path_stim=res_path, path_states=osjoin(res_path, net_id),
# pass_through_activation=True, activation=activation,
# add_bias=True, N_neurons=N,)
# # print()
# except FileNotFoundError as e:
# print(f'{net_id} is not integrated. Skipped in SVD.')
# continue
# else:
# print(f'{net_id} compuiting SVD...')
# _sig = svd(_X-_X.mean(axis=1)[:,np.newaxis],
# full_matrices=False, compute_uv=False)
# # _sig = np.arange(1, N+1)[::-1]
# # according to gast2024: (sum lambda) **2 / sum (lambda **2)
# # and note that s**2 = lambda
# eff_dim = sum(_sig**2)**2 / (sum(_sig**4))
# ax.semilogy(1-np.cumsum(_sig**2/sum(_sig**2)), color = newcolors[eval(net_id[2])],
# label=rf"$d_{{ {NET_ID_MAP[net_id][1:]} }}=$" +f'{eff_dim:.2f}'
# # label=fr"$d_{{{net_id}}}=$" +f'{eff_dim:.2f}',
# )
# # ax.semilogy(_sig/_sig[0], label=NET_ID_MAP[net_id], color = color)
# # ax.text(0.45, .95, r'$d_{eff}=$'
# # verticalalignment='center', horizontalalignment='center',
# # transform=ax.transAxes, fontsize=MFS*.5, )
# avg_labels = ['Instant.', 'Fast', 'Medium', 'Slow', 'Fixed']
# for ax, label in zip(axs[-1,:], avg_labels):
# ax.set_xlabel(r'$n$', fontsize=MFS*.8)
# ax.set_title(f'{label} neurons', fontsize=MFS*.8)
# for ax in axs[:,0]:
# # ax.set_ylabel(r'$\sigma_n$ / $\sigma_1$', fontsize=MFS*.8)
# ax.set_ylabel(r'Frac. Res. Var.', fontsize=MFS*.8)
# for ax in axs.ravel():
# ax.hlines(1e-2, 0, N, linestyle='--', color='k', lw=1)
# # ax.legend(loc='upper right')
# ax.legend(loc='lower left', fontsize=MFS*.6, ncols=1)
# ax.set_ylim(1e-8, 1.01)
# # ax.set_xlim(0, N/5)
# ax.set_xscale('log')
# polish(fig, [axs], 1., tight=True)
# plt.savefig(osjoin(res_path, f'svds{suffix}.svg'))
# plt.close('all')
# def plot_busy_svds(net_ids, res_path, N=None, also_hidden=False):
# # we only keep heterogenous networks to plot
# # ntau = utils.get_ntau_from_net_ids(net_ids)
# # hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
# # fig, axs = plt.subplots(1, ntau, figsize=(10, 2.))
# # for i, net_id in enumerate(hnet_ids):
# # print(net_id)
# # _y = np.load(osjoin(res_path, net_id, 'y.npy')).astype(float)
# # _sig = np.linalg.svd(_y-_y.mean(axis=1)[:,np.newaxis], compute_uv=False)
# # _sig /= _sig[0]
# # axs[i // (ntau+1)].semilogy(_sig, label='net_id') # NET_ID_MAP[net_id])
# # axs[i // (ntau+1)].legend(ncol=2)
# ntau = utils.get_ntau_from_net_ids(net_ids)
# sqrtn = int(np.sqrt(len(net_ids)))
# print(ntau, sqrtn)
# # hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
# fig, axs = plt.subplots(ntau+2, ntau+2, figsize=(3.*(ntau+2), 3*(ntau+2)),
# gridspec_kw={'width_ratios':[0.25, 1, 1, 1, 0.25]})
# cmaps = [plt.cm.Reds, plt.cm.Blues, plt.cm.Greens, plt.cm.Purples, plt.cm.Grays]
# cmaps = [cmap(np.linspace(0.2, 1 , ntau+2)) for cmap in cmaps]
# for net_id in net_ids:
# row_id, col_id = int(net_id[1]), int(net_id[0])
# color = cmaps[int(net_id[3])][int(net_id[2])]
# # alpha = np.clip(.2 + 0.8 * float(net_id[2])/3, 0, 1)
# print(net_id, (row_id, col_id))
# ax = axs[row_id, col_id]
# _y = np.load(osjoin(res_path, net_id, 'y.npy')).astype(float)
# _sig = np.linalg.svd(_y-_y.mean(axis=1)[:,np.newaxis], compute_uv=False)
# _sig /= _sig[0]
# ax.loglog(_sig, label=net_id, color = color, lw=1)#, alpha=alpha)
# for ax in axs[-1,:]:
# ax.set_xlabel(r'$n$', fontsize=MFS*.8)
# for ax in axs[:,0]:
# ax.set_ylabel(r'$\sigma_n$ / $\sigma_1$', fontsize=MFS*.8)
# for ax in axs.ravel():
# ax.set_ylim(1e-3, 1)
# # ax.set_ylim(-1,5)
# if ax in axs[:,0] or ax in axs[:,-1]:
# ax.legend(ncol=1, fontsize=6, )
# else:
# ax.legend(ncol=4, fontsize=6, )
# # ax.set_ylim(-1,5)
# polish(fig, [axs], 1., tight=False)
# plt.savefig(osjoin(res_path, 'svds_busy_all.svg'))
# plt.close('all')
# if also_hidden:
# assert N is not None
# fig, axs = plt.subplots(ntau+2, ntau+2, figsize=(3.*(ntau+2), 3*(ntau+2)),
# gridspec_kw={'width_ratios':[0.25, 1, 1, 1, 0.25]})
# cmaps = [plt.cm.Reds, plt.cm.Blues, plt.cm.Greens, plt.cm.Purples, plt.cm.Grays]
# cmaps = [cmap(np.linspace(0.2, 1 , ntau+2)) for cmap in cmaps]
# for net_id in net_ids:
# row_id, col_id = int(net_id[1]), int(net_id[0])
# color = cmaps[int(net_id[3])][int(net_id[2])]
# print(net_id, (row_id, col_id))
# ax = axs[row_id, col_id]
# _y = np.load(osjoin(res_path, net_id, 'y.npy')).astype(float)[:N,:]
# _sig = np.linalg.svd(_y-_y.mean(axis=1)[:,np.newaxis], compute_uv=False)
# _sig /= _sig[0]
# ax.loglog(_sig, label=net_id, color = color, lw=1)#, alpha=alpha)
# for ax in axs[-1,:]:
# ax.set_xlabel(r'$n$', fontsize=MFS*.8)
# for ax in axs[:,0]:
# ax.set_ylabel(r'$\sigma_n$ / $\sigma_1$', fontsize=MFS*.8)
# for ax in axs.ravel():
# ax.set_ylim(1e-3, 1)
# # ax.set_ylim(-1,5)
# if ax in axs[:,0] or ax in axs[:,-1]:
# ax.legend(ncol=1, fontsize=6, )
# else:
# ax.legend(ncol=4, fontsize=6, )
# polish(fig, [axs], 1., tight=False)
# plt.savefig(osjoin(res_path, 'svds_busy_v.svg'))
# plt.close('all')
[docs]
def plot_svds(net_ids, res_path, N, activation='possig', conv_params=None, suffix=''):
# we only keep heterogenous networks to plot
# ntau = utils.get_ntau_from_net_ids(net_ids)
# hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
# fig, axs = plt.subplots(1, ntau, figsize=(10, 2.))
# for i, net_id in enumerate(hnet_ids):
# print(net_id)
# _y = np.load(osjoin(res_path, net_id, 'y.npy')).astype(float)
# _sig = np.linalg.svd(_y-_y.mean(axis=1)[:,np.newaxis], compute_uv=False)
# _sig /= _sig[0]
# axs[i // (ntau+1)].semilogy(_sig, label='net_id') # NET_ID_MAP[net_id])
# axs[i // (ntau+1)].legend(ncol=2)
ntau = utils.get_ntau_from_net_ids(net_ids)
sqrtn = int(np.sqrt(len(net_ids)))
print(ntau, sqrtn)
# hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
fig, axs = plt.subplots(sqrtn, ntau+2, figsize=(3.*(ntau+2), 2*sqrtn))
row_mapping = {
'00': 0,
'10': 1,
'11': 2,
'12': 3,
'13': 4,
'20': 5,
'21': 6,
'22': 7,
'23': 8,
'30': 9,
'31': 10,
'32': 11,
'33': 12,
'40': 13,
}
cmaps = [plt.cm.Blues, plt.cm.Greens, plt.cm.Reds, plt.cm.Grays]
for net_id in net_ids:
row_base, col_id = int(net_id[1]), int(net_id[0])
cmap_id = int(net_id[3])
color_id = int(net_id[2])
row_id = row_mapping[net_id[1]+net_id[3]]
print(net_id, (row_id, col_id))
ax = axs[row_id, col_id]
# cmap = cmaps[cmap_id](np.linspace(0.2,1 , ntau+1))
color = cmaps[color_id]
# for row_id, mu_w in range(ntau+1):
# for col_id, mu_v in range(ntau+1):
# for i, net_id in enumerate(net_ids):
if suffix != '':
_, _, _X = spike_loader(data_path=res_path, net_id=net_id, spikes_name=f'spikes{suffix}',
convolve = True, conv_params=conv_params,
add_bias = True, N_neurons=N, )
else:
_, _, _X = ds_loader(path_stim=res_path, path_states=osjoin(res_path, net_id),
pass_through_activation=True, activation=activation,
add_bias=True, N_neurons=N,)
# print(net_id)
#_y = np.load(osjoin(res_path, net_id, 'y.npy')).astype(float)[:20,:]
_sig = svd(_X-_X.mean(axis=1)[:,np.newaxis], compute_uv=False, full_matrices=False)
# _sig /= _sig[0]
# # axs[i // sqrtn, i%sqrtn].semilogy(_sig, label=net_id)
# # axs[i // sqrtn, i%sqrtn].legend(ncol=2)
ax.plot(np.cumsum(_sig**2/sum(_sig**2)), label=net_id, )#color = color)
# ax.loglog(_sig, label=net_id, color = color)
# ax.text(color_id, col_id, net_id, rotation=90, color=color, fontsize=8)
for ax in axs[-1,:]:
ax.set_xlabel(r'$n$', fontsize=MFS*.8)
for ax in axs[:,0]:
# ax.set_ylabel(r'$\sigma_n$ / $\sigma_1$', fontsize=MFS*.8)
ax.set_ylabel(r'Frac. Expl. Var.', fontsize=MFS*.8)
for ax in axs.ravel():
ax.set_ylim(1e-3, 1.05)
ax.legend(fontsize=6)
# ax.set_ylim(-1,5)
polish(fig, [axs], 1., tight=False)
plt.savefig(osjoin(res_path, 'svds{suffix}.svg'))
plt.close('all')
[docs]
def plot_connectivity(data_path, rlz_id=None):
metadata = utils.read_metadata(data_path)
N = metadata['net']['N'] # network size
p = metadata['net']['p'] # sparsity
f = metadata['net']['f'] # fraction of excitatory neurons
w_ff, w, w_fb = utils.load_weights(data_path)
A_adj = utils.load_topology(data_path)
wmax = max(abs(w))
posts, pres = A_adj.nonzero()
Ne = int(N*f) # number of excitatory neurons
pres_e = pres<Ne
we = w[pres_e]
wi =w[~pres_e]
fig, axs = plt.subplots(2,1, figsize=(4,4), height_ratios=[.8, .2],)
axs[0].imshow(A_adj._with_data(w).todense(), cmap='seismic', vmin=-wmax, vmax=wmax, aspect='equal')
axs[0].set_ylabel('post')
axs[0].set_xlabel('pre')
# axs[0].set_axis_off()
axs[1].hist(we, bins=100, density=False, color='r', histtype='step', label='E');
axs[1].hist(wi, bins=100, density=False, color='b', histtype='step', label='I');
axs[1].legend(loc='upper left', fontsize=8)
axs[1].set_xlabel('weight value')
axs[1].set_ylabel('count')
polish(fig, axs)
# keep the top (square, aspect='equal') and bottom panels on the same width
fig.canvas.draw()
top = axs[0].get_position()
bot = axs[1].get_position()
axs[1].set_position([top.x0, bot.y0, top.width, bot.height])
if rlz_id is None:
name = 'connectivity.png'
else:
name = f'connectivity-rlz{rlz_id:03d}.png'
plt.savefig(osjoin(data_path, name))
plt.close('all')
[docs]
def plot_tau_dist(data_path, bins=50, rlz_id=None):
"""Plots the distribution of the neuronal time constants of every network
ID of the realization into the result path. Histograms are ``step`` type
with raw counts (``density=False``) on a log ``tau`` axis, colored by
the heterogeneity (``h``) level -- the second digit of the ``net_id`` --
following ``hetero_colors`` and labeled by ``HETERO_LABEL_MAP``. For the
homogeneous (``h=0``) networks, whose time constants are all equal, a vertical line at the constant value is
drawn instead of a histogram. The ``n`` levels (first digit, i.e. the mean
``E[tau]`` level per ``TAU_MAP``) are differentiated by line width. A
seaborn-style legend built from dummy handles reports the ``h`` colors and
the ``n`` widths. The x/y ranges are fixed across realizations: x in
[5e-3, 50] and y in [N/1000, N/2] with ``N`` the number of neurons per
network."""
net_ids = sorted([f.name for f in os.scandir(data_path) if f.is_dir()])
n_levels = sorted(set(int(net_id[0]) for net_id in net_ids))
h_levels = sorted(set(int(net_id[1]) for net_id in net_ids))
taus = {net_id: np.load(osjoin(data_path, net_id, 'taus.npy'))
for net_id in net_ids}
N = len(next(iter(taus.values())))
x_lo, x_hi = 5e-3, 50
bins = np.logspace(np.log10(x_lo), np.log10(x_hi), bins)
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
for i, n in enumerate(n_levels):
lw = 1 + i
for net_id in net_ids:
if int(net_id[0]) != n:
continue
color = hetero_colors[int(net_id[1])]
if int(net_id[1]) == 0:
ax.axvline(taus[net_id][0], color=color, lw=lw)
else:
ax.hist(taus[net_id], bins=bins, density=False,
histtype='step', lw=lw, color=color)
# seaborn-style legend: dummy handles, colored lines for the h levels,
# black lines of increasing width for the E[tau] (n) levels per TAU_MAP
h_handles = [mpl.lines.Line2D([], [], color=hetero_colors[h], lw=2,
label=HETERO_LABEL_MAP[str(h)])
for h in h_levels]
n_handles = [mpl.lines.Line2D([], [], color='k', lw=1 + i,
label=TAU_MAP[str(n)])
for i, n in enumerate(n_levels)]
ax.legend(handles=h_handles + n_handles, fontsize=9, loc='upper right',
frameon=True, borderpad=0.4, labelspacing=0.3)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(x_lo, x_hi)
ax.set_ylim(N / 1000, N / 2)
ax.set_xlabel(r'$\tau$', fontsize=16)
ax.set_ylabel('Counts', fontsize=11)
polish(fig, ax)
if rlz_id is None:
name = 'tau_dist.png'
else:
name = f'tau_dist-rlz{rlz_id:03d}.png'
plt.savefig(osjoin(data_path, name))
plt.close('all')
# def plot_w_norm(net_ids, res_path, suffix=''):
# # we only keep heterogenous networks to plot
# # ntau = utils.get_ntau_from_net_ids(net_ids)
# # hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
# # fig, axs = plt.subplots(1, ntau, figsize=(10, 2.))
# nv, nw, nhn, nhw = utils.get_net_levels_from_net_ids(net_ids)
# hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(nv-1)]]
# fig, axs = plt.subplots(1, nv-2, figsize=(10, 2.))
# for i, net_id in enumerate(hnet_ids):
# print(net_id)
# _w = np.load(osjoin(res_path, net_id, f'W_tayloramus{suffix}.npy'))
# _w = np.linalg.norm(_w, axis=-1) # norm across components #[:,:-1] # the last one is bias term
# _w = np.linalg.norm(_w, axis=-1) # norm across neurons
# # axs[i // (ntau+1)].semilogy(_w, lw=1, label=net_id,)#NET_ID_MAP[net_id], )
# # axs[i // (ntau+1)].legend(ncol=2)
# axs[i // nhn].semilogy(_w, lw=1, label=NET_ID_MAP[net_id], )
# axs[i // nhn].legend(ncol=2)
# for ax in axs:
# ax.set_xlabel('training id', fontsize=MFS*.8)
# axs[0].set_ylabel(r'$|W|$', fontsize=MFS*.8)
# polish(fig, [axs], 1., tight=False)
# plt.savefig(osjoin(res_path, f'w_norm{suffix}.svg'))
# plt.close('all')
[docs]
def plot_states_map(net_ids, data_path, activation='possig', fe=0.8):
# we only keep heterogenous networks to plot
ntau = utils.get_ntau_from_net_ids(net_ids)
hnet_ids = [net_id for net_id in net_ids if net_id[0] not in ['0', str(ntau+1)]]
t = np.load(osjoin(data_path, 't.npy'))
dt = t[-1] - t[-2]
u = np.load(osjoin(data_path,'u.npy'))
w_in = np.load(osjoin(data_path,'w_in.npy'))
W = load_npz(osjoin(data_path, 'A_adj.npz'))
W = W._with_data(np.load(osjoin(data_path, 'w0.npy')))
J = np.load(osjoin(data_path,'J.npy'))
delays = np.load(osjoin(data_path, 'delays.npy'))
N = W.shape[0]
Ne = int(N*fe)
#set_trace()
if len(delays.shape)>0:
print('Synaptic delays are different. Shifting u_ext for each neuron.')
delays = np.rint(delays/dt).astype(int)
u = np.array([np.roll(u[i], delays[i]) for i in range(u.shape[0])])
Iext = np.einsum('ik,kt->ti', w_in, u)
# for hists
fig1, axs1 = plt.subplots(int(len(hnet_ids)//ntau), ntau,
figsize=(len(hnet_ids)//ntau*3.5, (ntau)*3.5),
sharey=True, sharex=True)
# for crosscorr
fig2, axs2 = plt.subplots(int(len(hnet_ids)//ntau), ntau,
figsize=(len(hnet_ids)//ntau*2, (ntau)*3.5),
sharey=True, sharex=True)
# for autocorr
fig3, axs3 = plt.subplots(int(len(hnet_ids)//ntau), ntau,
figsize=(len(hnet_ids)//ntau*3.5, (ntau)*3.5),
sharey=True, sharex=True)
counts_ext, bins_ext = np.histogram(Iext.ravel(), bins=50, density=True,)
for i, net_id in enumerate(hnet_ids):
print(net_id)
loc = i%(ntau+1), i//(ntau+1) # ax locator
for ax in [axs1, axs2, axs3]:
ax[loc].set_title(net_id)
# reading tau to sort by tau
taus = np.load(osjoin(data_path, net_id, 'taus_v.npy'))
sidx = np.argsort(taus)
v = np.load(osjoin(data_path, net_id, 'y.npy')).astype(float)
f = activation_func(v, activation=activation)
Isyn = (W @ (J*f)).T
# histograms
counts_syn, bins_syn = np.histogram(Isyn.ravel(), bins=50, density=True,)
counts_e, bins_e = np.histogram(v[:Ne,:].ravel(), bins=50, density=True,)
counts_i, bins_i = np.histogram(v[Ne:,:].ravel(), bins=50, density=True,)
l1 = axs1[loc].stairs(counts_ext/max(counts_ext), bins_ext, alpha=0.8, label='External')
l2 = axs1[loc].stairs(counts_syn/max(counts_syn), bins_syn, alpha=0.8, label='Synaptic')
l3 = axs1[loc].stairs(-counts_e/max(counts_e), bins_e, alpha=0.8, label='Voltage E')
l4 = axs1[loc].stairs(-counts_i/max(counts_i), bins_i, alpha=0.8, label='Voltage I')
curr_legend = axs1[loc].legend(handles=[l1,l2], loc='upper right')
axs1[loc].add_artist(curr_legend)
axs1[loc].legend(handles=[l3,l4], loc='lower right')
axs1[loc].set_ylabel('Normalized density')
# cross-corrs
Q = np.hstack((Isyn[sidx], Iext[sidx], v.T[sidx]))
C = np.corrcoef(Q.T)
axs2[loc].imshow(C, vmin=-1, vmax=1, cmap='Blues')
for i in range(3):
# axs2[loc].vlines(C.shape[0]/3*i + C.shape[0]/3*fe, 0, C.shape[0], ls='--', lw=.5, color='grey')
# axs2[loc].hlines(C.shape[0]/3*i + C.shape[0]/3*fe, 0, C.shape[0], ls='--', lw=.5, color='grey')
axs2[loc].vlines(C.shape[0]/3*i, 0, C.shape[0], ls='--', color='k', lw=1)
axs2[loc].hlines(C.shape[0]/3*i, 0, C.shape[0], ls='--', color='k', lw=1)
axs2[loc].set_xlim(0, C.shape[0])
axs2[loc].set_ylim(0, C.shape[0])
every = 50
# power spectrum
# samples = np.random.choice(Isyn.shape[1], 10, replace=False)
freqxx, Psyn = periodogram(Isyn.T, fs=1./dt, )
freqxx, Pext = periodogram(Iext.T, fs=1./dt, )
freqxx, Uxx = periodogram(u, fs=1./dt, )
freqxx, Vxx = periodogram(v, fs=1./dt, )
# Psyn = (Psyn.T / np.max(Psyn, axis=1)).T
# Pext = (Pext.T / np.max(Pext, axis=1)).T
# Uxx = (Uxx.T / Uxx.max(axis=1)).T
# Vxx = (Vxx.T / Vxx.max(axis=1)).T
Psyn = Psyn.mean(axis=0)
Pext = Pext.mean(axis=0)
Vxx = Vxx.mean(axis=0)
Psyn /= Psyn.max()
Pext /= Pext.max()
Vxx /= Vxx.max()
fmean = np.mean(1/taus)
fstd = np.std(1./taus)
if fstd>1e-6:
axs3[loc].fill_between([fmean-fstd, fmean+fstd], [1,1],
color='r', alpha=0.1, label='$\\pm std$')
axs3[loc].vlines(fmean, 0, 1, ls='-.', lw=1, color='r',label=r'$\langle \tau^{-1} \rangle$')
# for dim in range(Uxx.shape[0]):
# axs3[loc].semilogy(freqxx[::every], Uxx[dim][::every], lw=.5,
# label=f'input {dim+1}', color=f'C{dim+2}')
shifter = 0
for name, psd, color in zip(['ext', 'syn','V'],
[Pext, Psyn, Vxx],
['k', 'C0', 'C1']):
# axs3[loc].semilogy(freqxx[::every], psd.mean(axis=0)[::every],
# lw=1, color=color)
axs3[loc].semilogy(freqxx[::every], psd[::every],
lw=1, color=color)
sub_freq = freqxx[freqxx>fmean]
# sub_psd = psd.mean(axis=0)[freqxx>fmean]
sub_psd = psd[freqxx>fmean]
p_opt = utils.linear_fit(sub_freq, sub_psd,
logtranfrom_x=True, logtranfrom_y=True,
bounds=((-np.inf, -np.inf), (1, np.inf)))
freqs = [sub_freq[0], sub_freq[-1]]
ys = p_opt[0]*np.log(freqs)+p_opt[1]
ys = np.exp(ys - 1.5*shifter) # shift down for visualization
axs3[loc].semilogy(freqs, ys, ls='--', color=color)
axs3[loc].plot([], [], lw=0.5, c=color, label=name + f'({p_opt[0]:.3f})')
shifter += 1
# axs3[loc].semilogy(freqxx[::every], Pext.mean(axis=0)[::every], lw=1, color='k')
# axs3[loc].semilogy(freqxx[::every], Psyn.mean(axis=0)[::every], lw=1, color='C1')
# axs3[loc].semilogy(freqxx[::every], Vxx.mean(axis=0)[::every], lw=1, color='C0')
# axs3[loc].fill_between(freqxx[::every],
# Pext.mean(axis=0)[::every] + Pext.std(axis=0)[::every],
# Pext.mean(axis=0)[::every] - Pext.std(axis=0)[::every],
# lw=1, color='k', alpha=0.1)
# axs3[loc].fill_between(freqxx[::every],
# Psyn.mean(axis=0)[::every] + Psyn.std(axis=0)[::every],
# Psyn.mean(axis=0)[::every] - Psyn.std(axis=0)[::every],
# lw=1, color='C1', alpha=0.1)
# axs3[loc].fill_between(freqxx[::every],
# Vxx.mean(axis=0)[::every] + Vxx.std(axis=0)[::every],
# Vxx.mean(axis=0)[::every] - Vxx.std(axis=0)[::every],
# lw=1, color='C0', alpha=0.1, )
# axs3[loc].plot([], [], lw=0.5, c='k', label='ext')
# axs3[loc].plot([], [], lw=0.5, c='C1', label='syn')
# axs3[loc].plot([], [], lw=0.5, c='C0', label='V')
#tf = 1/np.sqrt(1+(np.outer(taus, freqxx))**2)
# axs3[loc].semilogy(freqxx, tf.mean(axis=0), color='r', ls='--')
# axs3[loc].set_xlim(0,.3)
axs3[loc].set_xscale('log')
axs3[loc].set_ylim(1e-15, 2)
axs3[loc].legend()
axs3[loc].set_xlabel('Freq (Hz)')
axs3[loc].set_ylabel('Normalized PSD')
polish(fig1, axs1, 1.)
polish(fig2, axs2, 1.)
polish(fig3, axs3, 1.)
fig1.savefig(osjoin(data_path, f'state_hist.svg'))
fig2.savefig(osjoin(data_path, f'state_ccorr.svg'), dpi=300)
fig3.savefig(osjoin(data_path, f'state_acorr.svg'))
plt.close('all')
[docs]
def plot_complexity_map(data_path):
"""Creates and saves a complexity map per input dimension without recomputing
the tasks.
Rather than re-synthesizing the target timeseries and recomputing
information-theoretic measures for every subtask, this function simply loads
the previously calculated complexities from
``{data_path}/complexities.csv``, the file that
:meth:`hetero.sim.Simulation.compute_complexity` (via
``utils.compute_complexity``) has already saved during initialization.
Parameters
----------
task : hetero.tasks.TemporalTask
The task object (its ``name`` is used in the output file name).
data_path : str
Path to the simulation directory containing ``complexities.csv``.
Notes
-----
- ``complexities.csv`` contains one row per *flattened* subtask with columns
``dim, delta, deg, complexity``.
- A single figure is written to ``{task.name}_complexity.svg`` with one row
per input dimension: each row is a pcolormesh of the complexity with
``delta`` on the x-axis and ``degree`` on the y-axis.
- Tasks that do not vary along both a degree axis and a delta axis (e.g.
``Taylor`` or ``Nostradamus`` alone) have no 2D map to draw and are
skipped with a warning.
"""
path_csv = osjoin(data_path, 'complexities.csv')
if not os.path.exists(path_csv):
print(f'plot_complexity_map: {path_csv} does not exist -- skipping.')
return
df = pd.read_csv(path_csv, index_col=0)
if not {'dim', 'delta', 'deg'}.issubset(df.columns):
print('plot_complexity_map: complexities.csv lacks the columns '
'(dim, delta, deg). Skipping.')
return
dims = np.sort(df['dim'].unique())
if len(dims) == 0:
print('plot_complexity_map: no dimensions found in complexities.csv -- skipping.')
return
fig, axs = plt.subplots(len(dims), 1, figsize=(4, 2.2 * len(dims)),
sharex=True, squeeze=False)
for ax, dim_id in zip(axs[:, 0], dims):
sub = df[df['dim'] == dim_id].dropna(subset=['complexity'])
piv = sub.pivot(index='deg', columns='delta', values='complexity')
piv = piv.sort_index(axis=0).sort_index(axis=1)
m = ax.pcolormesh(piv.columns, piv.index, piv.values, edgecolor='face')
ax.set_ylabel('degree')
ax.set_title(f'input dim {int(dim_id)}', fontsize=9)
fig.colorbar(m, ax=axs, location='right', label='Complexity')
axs[-1, 0].set_xlabel(r'$\Delta$')
fig.savefig(osjoin(data_path, f'complexity.png'))
plt.close('all')
[docs]
def plot_udf_vs_complexity(data_path):
"""Plots the input-based baseline scores against the precomputed
complexities, per subtask.
Loads ``{data_path}/complexities.csv`` and ``{data_path}/udf.csv``, merges
them on ``(dim, delta, deg)`` and draws, for each subtask, the R2 scores of
the linear input-to-target decoder (``udf``) against the associated
complexity (``1 - |cos(u, y)|``). A single figure is saved as
``{task.name}_udf_vs_complexity.svg``.
The layout is a grid of subtask properties (rows) vs. data splits (columns):
each row colors the points by a different property of the subtask (input
dimension, delta shift ``\\Delta``, degree), the two columns are the train
and test scores, and a colorbar at the right of every row identifies the
coloring.
The udf scores (higher = easier) quantify how much of each target subtask is
*linearly reachable from the input alone*; the complexity (higher = harder)
is the static cosine-based difficulty. The dashed ``complexity = 1 - score``
line is drawn on every panel as the reference of perfect agreement.
For a linear (OLS) input-to-target decoder and a zero-mean, univariate
input, the two measures are translatable: ``score = |cos|^2 = (1 -
complexity)^2``. This mapping is drawn as a dash-dotted curve on every
panel.
Parameters
----------
task : hetero.tasks.TemporalTask
The task object (its ``name`` is used in the output file name).
data_path : str
Path to the simulation directory containing ``complexities.csv`` and
``udf.csv``.
"""
path_cmp = osjoin(data_path, 'complexities.csv')
path_udf = osjoin(data_path, 'udf.csv')
if not (os.path.exists(path_cmp) and os.path.exists(path_udf)):
print('plot_udf_vs_complexity: complexities.csv and/or udf.csv '
'missing -- skipping.')
return
cmp = pd.read_csv(path_cmp, index_col=0)
udf = pd.read_csv(path_udf, index_col=0)
m = cmp.merge(udf, on=['dim', 'delta', 'deg'], how='inner')
if len(m) == 0:
print('plot_udf_vs_complexity: no common subtasks found -- skipping.')
return
fig, axs = plt.subplots(3, 2, figsize=(6, 6), sharex=True, sharey=True, squeeze=False)
x_curve = np.linspace(0, 1, 300)
# rows: one coloring property per row -- dim, delta, degree.
# dim and deg are discrete labels and get a categorical colormap (one band
# per distinct value); delta is a continuous shift and stays diverging.
rows = [('dim', 'input dim', 'cat', 'Set2'),
('delta', r'$\Delta$', 'cont', 'coolwarm'),
('deg', 'degree', 'cat', 'Accent')]
cols = [('trn', 'train'), ('tst', 'test')]
for i, (prop, prop_label, kind, cmap_name) in enumerate(rows):
cvals = np.asarray(m[prop])
uniq = None
norm = None
if kind == 'cat':
uniq = np.unique(cvals)
base = plt.get_cmap(cmap_name).colors
cmap = ListedColormap([base[k % len(base)]
for k in range(len(uniq))])
norm = BoundaryNorm([uniq[0] - 0.5]
+ list((uniq[:-1] + uniq[1:]) / 2)
+ [uniq[-1] + 0.5],
ncolors=len(uniq), clip=True)
else:
cmap = cmap_name
for j, (key, col_title) in enumerate(cols):
ax = axs[i, j]
sc = ax.scatter(m['complexity'], m[key], s=5,
c=cvals, cmap=cmap, norm=norm, rasterized=True)
ax.plot(x_curve, (1 - x_curve)**2, ls='-.', c='k', lw=1.5,
label=r'(1-complexity)$^2 = \cos(u_k , y_{\theta})^2$'
)
if i == 0:
ax.set_title(col_title)
if i == len(rows) - 1:
ax.set_xlabel('complexity (= 1 - $| \\cos(u_k , y_{\\theta}) |$)')
if j == 0:
ax.set_ylabel('score (R$^2$ of $\\textbf{u} \\rightarrow y_{\\theta}$)')
if i == 0 and j == 0:
ax.legend(loc='lower left', fontsize=6.5)
# one colorbar per row (dim / delta / degree coloring)
if uniq is not None:
fig.colorbar(sc, ax=axs[i, :], label=prop_label, ticks=uniq,
format='%d' if np.all(uniq == np.round(uniq)) else None)
else:
fig.colorbar(sc, ax=axs[i, :], label=prop_label)
axs[0, 0].set_ylim(-0.15, 1.1)
# fig.suptitle(f'{task.name} | udf vs complexity')
polish(fig, axs, tight=False)
fig.savefig(osjoin(data_path, f'score_v_comp.svg'))
plt.close('all')
[docs]
def plot_difficulty_map_mi(task, u, data_path):
"""Information-theoretic difficulty map, kept alongside
:func:`plot_complexity_map` for comparison.
This version *recomputes* the target timeseries of every subtask and measures
the difficulty of the input-to-target mapping with mutual information
(``utils.compute_complexity_mi``): for each subtask and input dimension the
possibility, realizability and feasibility of the mapping are evaluated, and
a single figure is saved as ``{task.name}_feasibility.svg`` with one row per
input dimension and three columns (Possibility, Realizability, Feasibility).
Parameters
----------
task : hetero.tasks.TemporalTask
The task object. Works with both un-flattened subtasks (``[delta, deg]``
or ``[deg]``) and flattened ones (``[dim, delta, deg]``).
u : ndarray of shape (n_dims, n_timesteps)
The (compound) input timeseries.
data_path : str
Path to the simulation directory in which the figures are saved.
.. note::
Contrary to :func:`plot_complexity_map`, this function does *not* reuse
the precomputed ``complexities.csv`` and is therefore the expensive
reference against which the cheap (loaded) map is compared.
"""
u = u.astype(float)
n_dims = u.shape[0]
records = []
for subtask in task.subtasks:
# un-flattened subtasks are [degree] or [delta, degree]; flattened
# ones have the leading dim id.
if len(subtask) == 3:
dims = [int(subtask[0])]
delta, deg = subtask[1], subtask[2]
elif len(subtask) == 2:
dims = range(n_dims)
delta, deg = subtask[0], subtask[1]
else: # single degree, e.g., Taylor without flattening
dims = range(n_dims)
delta, deg = 0, subtask[0]
for dim_id in dims:
subu = u[dim_id, :]
suby = task.make_subtarget(u, delta, deg)[dim_id]
P, R, F = utils.compute_complexity_mi(subu, suby, return_prf=True)
records.append((int(dim_id), delta, deg, P, R, F))
df = pd.DataFrame(records, columns=['dim', 'delta', 'deg', 'P', 'R', 'F'])
dims = np.sort(df['dim'].unique())
n_dims = len(dims)
grid = [('P', 'Possibility'), ('R', 'Realizability'), ('F', 'Feasibility')]
fig, axs = plt.subplots(n_dims, 3, figsize=(8, 1.5 * n_dims),
sharex=True, sharey='row', squeeze=False)
for row, dim_id in enumerate(dims):
sub = df[df['dim'] == dim_id]
for col, (key, label) in enumerate(grid):
ax = axs[row, col]
piv = sub.pivot(index='deg', columns='delta', values=key)
piv = piv.sort_index(axis=0).sort_index(axis=1)
m = ax.pcolormesh(piv.columns, piv.index, piv.values, edgecolor='face')
plt.colorbar(m, ax=ax, location='top', label=label, shrink=0.8)
if row == n_dims - 1:
ax.set_xlabel(r'$\Delta$')
axs[row, 0].set_ylabel(f'dim {int(dim_id)}\ndegree')
fig.suptitle(task.name)
fig.savefig(osjoin(data_path, f'{task.name}_feasibility.svg'))
plt.close('all')
# def plot_prediction_trace(sim, net_id, subtask_ids=None, subtasks=None):
# assert (subtask_ids is not None) or (subtasks is not None), 'no subtask is specified.'
# task = sim.tasks[0]
# if subtask_ids is None:
# subtask_ids = []
# for subtask in subtasks:
# subtask_ids.append(task.get_subtask_index(subtask))
# t,u,x = ds_loader(path_stim=sim.data_path,
# path_states=osjoin(sim.data_path, net_id),
# N_neurons=sim.N, add_bias = True,
# pass_through_activation = True, activation=sim.activation,
# )
# W = np.load(osjoin(sim.res_path, net_id, f'W_{task.name}.npy'))
# Ntrials = int(len(W)/len(task.subtasks)) # number ot trials per subtask
# for subtask_id in subtask_ids:
# print(f'Task id {subtask_id} with params {task.subtasks[subtask_id]} is being plotted')
# y = task.make_target(t,u,x, *task.subtasks[subtask_id])
# w = W[Ntrials * subtask_id]
# yhat = x @ w
# fig, axs = plt.subplots(w.shape[1],1, sharex=True)
# for cidx in range(w.shape[1]):
# ax = axs[cidx]
# ax.plot(t[task.L_trn:], y[task.L_trn:, cidx], label='GT')
# ax.plot(t[task.L_trn:], yhat[task.L_trn:, cidx], label='Pred.')
# ax.set_xlim(t[task.L_trn], t[task.L_trn]+500)
# ax.set_xlabel('time')
# ax.set_ylabel(f'comp. {cidx+1}')
# axs[0].legend(ncols=2, loc='upper left')
# polish(fig, axs,1)
# plt.savefig(osjoin(sim.res_path, f'pred_{net_id}_subtask_id_{subtask_id}.svg'),
# bbox_inches='tight', dpi=100)
# plt.close('all')
# def plot_mismatch(task, data_path, activation, delta=[0], degree= 1, v=1, w=4,
# suffix='', conv_params=None, spikes_name=''):
# # net_ids = [f'{v}4{i}0' for i in range(5)]
# net_ids = sorted([d for d in os.listdir(osjoin(data_path)) if os.path.isdir(osjoin(data_path,d))])
# net_ids = [net_id for net_id in net_ids if net_id[0]==f'{v}']
# nhet = len(set([net_id[2] for net_id in net_ids]))
# colors = plt.cm.Blues(np.linspace(.2,.9, nhet))
# N = task.N
# s = slice(task.L_trn, None, 10)
# for i, net_id in enumerate(net_ids):
# if suffix != '':
# t, u, X = spike_loader(data_path=data_path, net_id=net_id, spikes_name=spikes_name,#f'spikes{suffix}',
# convolve = True, conv_params=conv_params,
# add_bias = True, N_neurons=N, )
# else:
# t, u, X = ds_loader(path_stim=data_path, path_states=osjoin(data_path, net_id),
# pass_through_activation=True, activation=activation,
# add_bias=True, N_neurons=N,)
# W = np.load(osjoin(data_path, net_id, f'W_tayloramus{suffix}.npy'))
# ntrials = int(round(W.shape[0]/len(task.subtasks)))
# subtask = [delta, degree]
# subtask_id = task.get_subtask_index(subtask)
# y = task.make_target(t, u, X, *subtask)
# w = W[subtask_id*ntrials]
# yhat = X@w
# if i == 0:
# ndim = W.shape[-1]
# fig, axs = plt.subplots(ndim, 1, figsize=(9, 2.5*ndim), sharex=True)
# axs = np.array(axs).reshape(ndim)
# for dim in range(ndim):
# ax = axs[dim]
# ax.plot(t[s], yhat[s, dim], label=NET_ID_MAP[net_id][1:],
# lw=1.5, color = colors[i])
# # plt.savefig(f'mismatch_{i+1}.svg', bbox_inches='tight')
# # axs[0, 0].set_ylabel(viz.NET_ID_MAP[net_id])
# # axs[i, dim].set_ylim(-2,2)
# for dim in range(ndim):
# axs[dim].plot(t[s], y[s, dim], label='ground truth',
# lw=1.5, color='k', ls='--' , zorder=99)
# axs[dim].set_xlabel('time')
# axs[dim].set_title(f'dim {dim+1}')
# axs[dim].legend(ncol=1, bbox_to_anchor=(1, 1.1),
# title='level of\n heterog.', frameon=False)
# polish(fig, axs, 1.)
# plt.savefig(osjoin(data_path,f'mismatch{suffix}.svg'), bbox_inches='tight')
# plt.close()
# # axs[0, 0].legend(loc='upper right')
# #for _dim, ax in enumerate(axs[-1,:]):
# # ax.plot(t[s],y[s, _dim], lw=2, color='k', ls=':')#,label='GT')
# # ax.set_xlabel('time')
# # ax.legend(ncol=1, bbox_to_anchor=(1, 1.1),
# # title='level of\n heterog.', frameon=False)
# # for dim in range(sim.stim_dim):
# # axs[0, dim].set_title(f'dim {dim+1}')
# # axs[0].set_ylabel('disposition')
# # for ax in axs.ravel():
# # plt.savefig(f'mismatch_{name}.svg', bbox_inches='tight')
# # plt.close()
[docs]
def plot_uvx(u, v, x, res_path, final_chunk=1000, state_name=None, rlz_id=None):
"""
Plots the input ``u``, voltage ``v``, and activation ``x`` arrays in a 3x2
grid: the last ``final_chunk`` timesteps in the first column (heatmap when
the array has more than one feature, a line otherwise) and the flattened
distribution of each array as a histogram in the second column.
Arrays must be shaped (n_features, timesteps).
"""
cmaps = ['rainbow', 'PRGn', 'Grays']
labels = ['u', 'v', 'x']
fig, axs = plt.subplots(3, 2, layout='tight', figsize=(8, 6))
for i, (q, label, cmap) in enumerate(zip([u, v, x], labels, cmaps)):
q = np.asarray(q)
n_features = q.shape[0]
_q = q[:, -final_chunk:]
ax = axs[i, 0]
if n_features > 1:
mx = np.quantile(abs(_q), .95)
im = ax.pcolormesh(_q, cmap=cmap, vmin=-mx, vmax=mx, rasterized=True)
fig.colorbar(im, ax=ax).set_label(label)
else:
ax.plot(_q[0], rasterized=True)
ax.set_ylabel(label)
axs[i, 1].hist(q.ravel(), bins=25, density=True, alpha=0.75, rasterized=True)
axs[i, 1].set_ylabel('density')
axs[2, 0].set_xlabel('time step')
axs[2, 1].set_xlabel('value')
polish(fig, axs, 1)
if state_name is not None:
name = f'uvx-{state_name}.png'
elif rlz_id is None:
name = 'uvx.png'
else:
name = f'uvx-rlz{rlz_id:03d}.png'
plt.savefig(osjoin(res_path, name), bbox_inches='tight', dpi=150)
plt.close('all')