import os
osjoin = os.path.join # an alias for convenience
from pdb import set_trace
import time
import json
import numpy as np
import pandas as pd
from scipy.sparse import save_npz
from scipy.linalg import svd
from hetero.config import CFG
from hetero import utils
from hetero import tasks
from hetero import stimulus
from hetero import networks
from hetero import dataloader
from hetero import fitting
# from hetero.viz import plot_score_profile
from hetero.dynamics.li import LI
from hetero.dynamics.lif import LIF
"""
Heterogeneous Neural Network Simulation Framework
Main class that handles simulation of heterogeneous neural networks with:
- Configurable network architectures and parameters
- Multiple types of neural dynamics
- Various input stimulus types (chaotic, periodic, etc.)
- Different learning tasks and evaluation metrics
- Heterogeneous timescales for neurons and synapses
The framework supports:
1. Network Configuration
- Network size and connectivity
- Synaptic weights initialization
- Heterogeneous timescales
2. Input Processing
- Multiple stimulus types (Lorenz, Mackey-Glass, sinusoidal)
- Input preprocessing and scaling
- Noise addition
3. Task Definition
- Taylor series approximation
- Prediction tasks
- Temporal processing
4. Results Management
- Automatic naming based on parameters
- Results storage and loading
- Performance analysis
Parameters
----------
task_params : dict
Parameters defining the learning task
Example: {'tayloramus': {'delta': (-2,2,31), 'degree': (1, 6, 5)}}
stim_params : dict
Input stimulus configuration
Example: {
'lorenz-0': {},
'mackey_glass-3': {'tau': 1},
'sin-0': {'phase': 21, 'omega': 3}
}
heterog_params : dict
Parameters controlling network heterogeneity
Example: {'lognormal':{'n_means': 3, 'n_vars': 3, 'distance': np.sqrt(10)}}
net_params : dict
Network architecture parameters
Example: {
'N': 250, # Number of neurons
'p': 0.1, # Connection probability
'g': 4.0, # Global gain
'mue': 1., # Mean excitatory weight
'sig0': 1., # Weight variance
}
set_params : dict
Training/testing set parameters
Example: {
'n_trn': 20, # Training samples
'n_tst': 5, # Test samples
'n_trials': 3 # Number of trials
}
root : str, optional
Results directory path (default: './results/')
save_type : str, optional
Data type for saving results (default: 'float16')
base_tau : float, optional
Base timescale for dynamics (default: 1.0)
dt : float, optional
Integration timestep (default: 0.001)
seed : int, optional
Random seed (default: 2025)
"""
NCORES = os.cpu_count()
def _variance_keep(svals, var_frac):
"""Boolean mask over a singular-value spectrum keeping the modes whose
relative contribution to the total variance is at least ``var_frac``.
A mode ``k`` is kept when ``s_k^2 / sum_j s_j^2 >= var_frac``: the tail of
the spectrum below the fraction is numerically non-robust (its directions
are dominated by roundoff/noise) and is truncated away. A spectrum with
zero total variance keeps nothing.
"""
if svals.size == 0:
return np.zeros(0, dtype=bool)
total = float(np.sum(svals ** 2))
if total == 0:
return np.zeros(svals.size, dtype=bool)
return (svals ** 2) / total >= var_frac
[docs]
class Simulation:
"""An simulation environment that orchestrates performance quantification of
heteorgeneous networks through the followings steps:
1. sampling systematically controlled heterogeneous time constants
2. network architectures creation and weight initialization,
3. handling different stimulus types (choatic, periodic, and mixed),
4. parameterizing a task family
5. perparing the train and test datasets
Parameters
----------
task_params : dict
Task parameters (see hetero.tasks for details).
stim_params : dict
Stimulus parameters (see hetero.stimulus for details).
heterog_params : dict
Heterogeneity setting.
net_params : dict
Stimulus configuration (see hetero.networks for details).
set_params : dict
Dataset settings.
root : str, optional
Root path for saving the results, by default './results/'
save_type : str, optional
State data type. Default is 'float16' to save disk space.
base_tau : float, optional
The value assigned to the base period of the stimulus. It defiens the
scale of time, and all other temporal variables, including the neuronal
time constants. Default is 1.
dt : float, optional
Integration timestep in the unit of ``base_tau``. Default is .01, which
alongside a unit ``base_tau`` result in 100 samples per each base period.
seed : int, optional
Random seed for reproduction. Default is 2025.
skip_init_plot : bool, optional
If True, skip the two figures rendered during initialization (the
complexity map and the udf-vs-complexity scatter). Both read the
data already saved in the folder, so a parameter sweep -- where the
same folder is re-entered once per case -- would otherwise
re-render them from the previous cases after each single run.
Default is False. The underlying ``complexities.csv`` /
``udf.csv`` are written either way.
"""
def __init__(self,
task_params,
stim_params,
heterog_params,
net_params,
set_params,
uq_params,
root,
seed,
save_type='float16',
base_tau=1., dt=.01,
skip_init_plot=False,
# naming_pattern = 'NyfusJd'
):
self.stype = save_type
self.task_params = task_params
self.stim_params = stim_params
self.heterog_params = heterog_params
self.net_params = net_params
self.set_params = set_params
# self.uq_params = uq_params
self.n_trials = uq_params['n_trials']
self.n_realizations = uq_params['n_realizations']
# TODO: must adjust the load utilities accordingly.
self.rlz_naming = True # whether to add a suffix to random variables
self.metadata = dict(
net=self.net_params,
task=self.task_params,
stim=self.stim_params,
dataset=self.set_params,
heterog=self.heterog_params,
)
self.name = self.generate_sim_name()
self.root = root
self.res_path = self.get_result_path()
self.net_ids = self.generate_net_ids()
self.prepare_dirs()
print('='*27)
print(f'Simulation name: {self.name}.')
print(f'Result path : {self.res_path}')
print('='*27)
print()
# loads the previouse md if exists. Will be used internally
# by *setup* methods below to override arguments in favor of
# the previously set values.
self.read_metadata()
# will be overridden later by the previous md if exists
self.base_tau = base_tau
self.dt = dt
self.bs = int(self.base_tau/self.dt)
# handling rng and seed
self.setup_seed(seed)
self.rng = np.random.default_rng(self.seed)
print(f'rng state initialized @ {self.rng.bit_generator.state["state"]["state"]}')
# self.sync_rng() # syncs the state of rng with the prev rng
# retrieving the number of distinct realized systems, those
# that are integrated, and those that are benchmakred.
# self.num_realized = self.prev_md.get('num_realized', 0)
# self.num_integrated = self.prev_md.get('num_integrated', 0)
# self.num_benchmarked = self.prev_md.get('num_benchmarked', 0)
self.setup_integration_interval() # uses prev_md
self.setup_task() # uses -
self.setup_stim_n_target() # uses prev_md
self.compute_complexity() # uses -
# self.synthesize_noise() # uses num_realized
# self.synthesize_timescales() # uses num_realized
# self.synthesize_topology() # uses num_realized
# post-synthesis crafts
# self.update_num_realized() # updates only upon new realization
# self.save_rng_state() # overrides the old one
# integrate the old an new metadata
self.update_metadata()
# print(f'rng state after init @ {self.rng.bit_generator.state["state"]["state"]}')
# compute the input-based baseline scores (udf) and compare them with
# the precomputed complexities. Only the *figures* are skippable: the
# complexities and the udf scores are data, saved to disk here and
# reused by --replot, so they are always computed.
self.compute_udf()
if not skip_init_plot:
# the complexity map of each input dimension, from the
# complexities computed and saved above
self.plot_complexity_map()
self.plot_udf_vs_complexity()
[docs]
def generate_sim_name(self):
"""Creates a suitable name for the simulation based on the specified
configuration. Generally, adds a "key_value" token to the name if the
value of the "key" is different than the default configuration file.
If the `value` is not a scaler, but rather a dictionary or list, only
the "key" is appended to the name. The rational is that for configurations
that require more in-depth specification (e.g., stimulus or task
parameterization), a unique label (defined in the config file) should be
used.
Returns
-------
str
The generated simulation name.
"""
# TODO: change the naming utilities so that they surely add the key
# if the value is a dictionary.
param_dicts = dict(
net=self.metadata['net'],
task=self.metadata['task'],
dataset=self.metadata['dataset'],
heterog=self.metadata['heterog'],
stim=self.metadata['stim'],
)
# print('='*20)
# print('Hello from the sim_namer')
tokens = []
for param, param_dict in param_dicts.items():
default_cfg = CFG[param] if param in CFG.keys() else CFG
for key, val in param_dict.items():
# print(key, val)
# always encode the stimulus in the name
if param == 'stim':
token = f'-{key}'
tokens.append(token)
# print(f' I appended {token} for {key},{val}')
else:
if val != default_cfg[key]:
val_str = utils.canonicalize_value(val)
token = f"-{key}_{val_str}"
token = utils.sanitize_token(token)
tokens.append(token)
# print(f' I appended {token} for {key},{val}')
return ''.join(tokens)[1:] # drop the trailing -
[docs]
def generate_state_name(self, dyn, Ju, J, Jn, solver):
"""Generates a suitable name for the network state based on the
specified dynamic ``dyn`` and the gain values. This name is
used to save/load the network states.
Parameters
----------
dyn: str
The dynamic rule (e.g. 'LI' or 'LIF').
Ju: float
The feedforward gain.
J: float
The recurrent gain.
Jn: float
The noise gain (intensity).
solver: str
The integration solver.
Returns
-------
str
The generated state name.
"""
tokens = [dyn]
for key in ['Ju','J','Jn']:
val = eval(key)
if val != CFG[key]:
val_str = utils.canonicalize_value(val)
token = f"-{key}_{val_str}"
token = utils.sanitize_token(token)
tokens.append(token)
# if solver != CFG['solver']:
# tokens += f"-{solver}"
# we only add the realization id to the states if we
# need to generate more than one realizations
# if self.n_realizations>1:
# tokens.append(f'-rlz{self.num_realized:03d}')
# tokens.append(rlz_suffix)
return ''.join(tokens)
[docs]
def generate_readout_name(self, readout_dict):
"""Generates a suitable name for score files resulted from reading out
the states with readout parameters defined in ``readout_dict``. As before,
a "key_value" token is added to the name if the value of the "key" is
different from the default configuration.
Parameters
----------
readout_dict: dict
The readout dictionary specifying the followings items used to define
create a linear decoder (see ``fitting.Decoder`` for more details):
- readout method
- readout regularization
- readout bias term
and
- readout activation dictionary
to instruct the type of activation / convolution to be applied on the
states prior to reading out,
Returns
-------
str
The generated state name.
"""
tokens = []
for key, val in readout_dict.items():
if isinstance(val, dict): # readout_activation params
for k, v in val.items():
tokens.append(f'-{k}')
# if v != CFG[key][k]:
# v_str = utils.canonicalize_value(v)
# token = f"-{k.replace('readout_','')}_{v_str}"
# token = utils.sanitize_token(token)
# tokens.append(token)
else:
if val != CFG[key]:
val_str = utils.canonicalize_value(val)
token = f"-{key.replace('readout_','')}_{val_str}"
token = utils.sanitize_token(token)
tokens.append(token)
tokens = ''.join(tokens)
if len(tokens)>0:
tokens = tokens
return tokens
[docs]
def generate_optim_name(self, optim_dict):
"""Generates a suitable name for the optimization procedure given
``optim_dict``. A "key_value" token is added to the name if the
value of the "key" is different from the default configuration.
Parameters
----------
optim_dict: dict
The optimization dictionary specifying the followings items used
to define optimization algorithm (see ``fitting.Decoder`` for
more details):
- alg
- reg
- fit_bias
Returns
-------
str
The generated state name.
"""
tokens = []
for key, val in optim_dict.items():
if val != CFG[key]:
val_str = utils.canonicalize_value(val)
token = f"-{key}_{val_str}"
token = utils.sanitize_token(token)
tokens.append(token)
tokens = ''.join(tokens)
if len(tokens)>0:
tokens = tokens
return tokens
[docs]
def get_state_path(self, net_id, state_name, dyn=None):
ext = 'npy'
if dyn == 'LIF' or (state_name is not None and 'LIF' in state_name):
ext = 'pkl'
return osjoin(self.res_path, net_id, f'{state_name}.{ext}')
[docs]
def state_files_exist(self, state_name, dyn=None):
exists = [os.path.exists(self.get_state_path(net_id, state_name, dyn=dyn))
for net_id in self.net_ids]
if sum(exists) == len(self.net_ids):
return True
elif sum(exists) == 0:
return False
else:
return 'partial'
[docs]
def score_files_exist(self, score_name, mode='tst'):
return all(
os.path.exists(osjoin(self.res_path, net_id, f'{mode}_{score_name}.csv'))
for net_id in self.net_ids
)
[docs]
def generate_score_template(self, state_name, readout_params, optim_params, ntrial=''):
ro_name = self.generate_readout_name(readout_params)
optim_name = self.generate_optim_name(optim_params)
# return f'{state_name}{ro_name}{optim_name}-nsplit{ntrial}'
return f'{state_name}{ro_name}{optim_name}'
[docs]
def get_score_paths_for_template(self, net_id, template):
return (
osjoin(self.res_path, net_id, f'trn_{template}.csv'),
osjoin(self.res_path, net_id, f'tst_{template}.csv'),
)
[docs]
def get_score_paths(self, net_id, state_name, readout_params, optim_params, ntrial):
template = self.generate_score_template(
state_name=state_name,
readout_params=readout_params,
optim_params=optim_params,
ntrial=ntrial,
)
return self.get_score_paths_for_template(net_id, template)
[docs]
def scan_unbenchmarked(self, state_name, readout_params, optim_params, verbose=False):
net_id_dict = {}
for net_id in self.net_ids:
trials_to_go = []
for trial in self.n_trials:
template = self.generate_score_template(
state_name=state_name,
readout_params=readout_params,
optim_params=optim_params,
ntrial=trial,
)
path_trn_scr, path_tst_scr = self.get_score_paths_for_template(net_id, template)
if os.path.exists(path_tst_scr) and os.path.exists(path_trn_scr):
if verbose:
print(f'scores {template} for network {net_id} exist. Skipped fitting.')
else:
trials_to_go.append(trial)
if len(trials_to_go) > 0:
net_id_dict[net_id] = trials_to_go
return net_id_dict
[docs]
def benchmark_files_exist(self, state_name, readout_params, optim_params):
return len(
self.scan_unbenchmarked(
state_name=state_name,
readout_params=readout_params,
optim_params=optim_params,
verbose=False,
)
) == 0
[docs]
def progress_count(self, skip_benchmark=False):
return self.num_integrated if skip_benchmark else self.num_benchmarked
[docs]
def setup_seed(self, seed):
# overrides the seed from the previous md if exists
self.seed = self.prev_md.get('seed', seed)
# def update_num_realized(self, rlz_id):
# # increment num realization if it was really realized
# if self.needs_new_realization():
# self.num_realized = rlz_id
[docs]
def mark_integrated(self, rlz_id, state_name, dyn=None):
if not self.state_files_exist(state_name, dyn=dyn):
return False
self.num_integrated =rlz_id# max(self.num_integrated, self.num_realized)
self.update_metadata()
return True
[docs]
def mark_benchmarked(self, rlz_id, state_name, readout_params, optim_params):
if not self.benchmark_files_exist(state_name, readout_params, optim_params):
return False
self.num_benchmarked = rlz_id # max(self.num_benchmarked, self.num_realized)
self.update_metadata()
return True
[docs]
def get_param_from_net_id(self, net_id):
"""Given a network ID, returns the corresponding index for mean and std levels.
For example network '21' corresponds to 2nd index of mean and the 2nd index of
std levels (recall the we always have a std of 0).
Parameters
----------
net_id: str
The network ID string (e.g., '21').
Returns
-------
tuple
A tuple for (mean, std) indices.
"""
# let's first collect all ids (may start from 1 instead of 0)
mids= np.array(sorted(set([int(digits[0]) for digits in self.net_ids])))
sids = np.array(sorted(set([int(digits[1]) for digits in self.net_ids])))
mid, sid = int(net_id[0]), int(net_id[1])
return self.means[mid - mids.min()], self.stds[sid - sids.min()]
[docs]
def generate_net_ids(self):
"""Generates the network IDs based on the heterogeneity parameters. Network
IDS are ordinal numerical strings where the first and second digits respectively
corresponds to the mean and variance level in the time constants.
Returns
-------
list
A list of network ID strings.
.. note:: A zero-std level is always added to the each mean level to account for
the homogeneous network acting as the control.
"""
n_vars = self.heterog_params['n_vars'] + 1 # plus 1 for homogeneous case
n_means = self.heterog_params['n_means']
if n_means%2 == 0:
n_means += 1
means_idxs = np.arange(n_means) + 1 # we reserve mean=0 for instantenous neurons
vars_idxs = np.arange(n_vars)
return utils.make_net_ids(means_idxs, vars_idxs)
[docs]
def sync_rng(self, state_name, verbose=True):
msg ='No previous rng detected.'
path_pkl = osjoin(self.res_path, f"rng_{state_name}.pkl")
if os.path.exists(path_pkl):
utils.load_rng(self.rng, path_pkl)
msg ='Previous rng state loaded.'
if verbose:
print(msg)
[docs]
def save_rng_state(self, state_name):
path_pkl = osjoin(self.res_path, f"rng_{state_name}.pkl")
utils.save_rng(self.rng, path_pkl)
# def get_num_realized(self, verbose=True):
# num = self.prev_md.get('num_realized', 0)
# if verbose:
# print(f'Number of previously made realizations: {num}.')
# return num
[docs]
def needs_new_realization(self, rlz):
return rlz < self.n_realizations-1
# if self.rlz < self.n_realizations:
# return True
# else:
# return False
# if self.num_benchmarked < self.num_realized:
# return False
# elif self.num_realized < self.n_realizations:
# return True
# else:
# return False
[docs]
def rlz_suffix(self, rlz_id=0):
if not self.rlz_naming:
return ''
else:
return f'-rlz{rlz_id:03d}' # TODO: shouldn't it be - than _
[docs]
def setup_integration_interval(self):
"""Finds the (temporal) length of training and testing intervals based on the
number of required samples in the training and testing sets. These intervals
also depend on the task's specifics. If a task requires inputs from the past or
future, the integration horizon must be extended accordingly.
"""
# extract the task name and its parameters
tname, tparams = next(iter(self.task_params.items()))
# check if we have calculated the train and test lengths already
if 'Ntrn' in self.prev_md.keys():
self.Ntrn = self.prev_md['Ntrn']
self.Ntst = self.prev_md['Ntst']
self.Ltrn = self.prev_md['Ltrn']
self.Ltst = self.prev_md['Ltst']
else:
if tasks.needs_synthesis(tname):
# compute the number of samples in each sets
self.Ntst = int(self.set_params['n_tst'] * self.bs)
self.Ntrn = int(self.set_params['n_trn'] * self.bs)
if self.set_params['scale_by_size']:
self.Ntrn *= (1+self.net_params['N'])
# and then find the integration legnth for each set
if tname == 'taylor':
self.Ltrn = self.Ntrn
self.Ltst = self.Ntst
elif tname == 'tayloramus' or tname == 'nostradamus':
# for both of these tasks the past and future horizons have
# to be extended by `delt_max`
self.Ltrn = self.Ntrn + 2 * tparams['delt_max']*self.bs
self.Ltst = self.Ntst + 2 * tparams['delt_max']*self.bs
else:
raise ValueError(f'Task {tname} is not supported for synthesis yet.')
else:
# To be extended for other task types ...
pass
# and finally select the train and test indices
if tname == 'taylor':
self.Strn = slice(0, self.Ltrn)
self.Stst = slice(self.Ltrn, self.Ltrn + self.Ltst)
elif tname == 'tayloramus' or tname == 'nostradamus':
self.Strn = slice(0 + tparams['delt_max']*self.bs,
self.Ltrn - tparams['delt_max']*self.bs)
self.Stst = slice(self.Ltrn + tparams['delt_max']*self.bs,
self.Ltrn + self.Ltst - tparams['delt_max']*self.bs)
self.n_timesteps = self.Ltrn + self.Ltst # total number of timesteps required
print(f'Total number of timesteps: {self.n_timesteps}')
[docs]
def synthesize_timescales(self):
"""Generates and saves the heterogeneous timescales for all networks based on the
heterogeneity parameters.
"""
distribution = self.heterog_params['distro']
log_distance = self.heterog_params['log_distance']
net_ids = self.net_ids
n_means = len(set([net_id[0] for net_id in net_ids]))
n_vars = len(set([net_id[1] for net_id in net_ids])) - 1 # homog will be added manually
mean_range = n_means//2 * log_distance
var_range = n_vars//2 * log_distance
self.means = np.logspace(-mean_range, +mean_range, n_means) * self.base_tau
self.stds = np.logspace(-var_range, +var_range, n_vars) * self.base_tau
self.stds = np.concat(([0], self.stds)) # always add homogeneous case
# distribution = self.heterog_params.keys()[0]
print('generating timescales')
for net_id in self.net_ids:
path_tau = osjoin(self.res_path, net_id, f'taus.npy')
# if not os.path.exists(path_tau):
mu_, std_ = self.get_param_from_net_id(net_id)
taus = utils.sample_taus(self.net_params['N'],
mu=mu_, std=std_,
clip_below=self.dt/2,
distr=distribution,
rng = self.rng)
np.save(path_tau, taus)
[docs]
def setup_task(self):
"""Initializes the task object on the specified task parameters,
and further creates all subtasks."""
name, params = list(self.task_params.items())[0]
if name == 'tayloramus':
deg_min = params['deg_min']
deg_max = params['deg_max']
deg_n = params['deg_n']
delt_min = params['delt_min']
delt_max = params['delt_max']
delt_n = params['delt_n']
self.task = tasks.Tayloramus(
deg_min=deg_min, deg_max=deg_max,
delt_min=delt_min, delt_max=delt_max, # in unit of tau0
base_stride = self.bs, # necessary to convert them into steps
deg_n=deg_n, delt_n=delt_n, delt_spacing='log',
)
elif name == 'taylor':
deg_min = params['deg_min']
deg_max = params['deg_max']
deg_n = params['deg_n']
self.task = tasks.Taylor(
deg_min=deg_min, deg_max=deg_max, deg_n=deg_n,
)
elif name == 'nostradamus':
delt_min = params['delt_min']
delt_max = params['delt_max']
delt_n = params['delt_n']
self.task = tasks.Nostradamus(
delt_min=delt_min, delt_max=delt_max, # in unit of tau0
base_stride = self.bs, # necessary to convert them into steps
delt_n=delt_n, delt_spacing='log',
)
# elif name == 'ashd':
# self.task = tasks.ASHD(
# path=params['path'],
# ds_fac=params['ds_fac'],
# encoding=params['encoding'],
# rng=self.rng
# )
else:
raise ValueError(f'Task {name} is not a valid task name')
self.task.make_subtasks() # make tasks for all parameters
[docs]
def setup_stim_n_target(self):
"""Generates the (compound) stimulus and the target timeseries, and
updates the metadata with the stimulus information, if not created
already.
Different stimuli have distinct timescales. To standardize all cases
we scale time so that the most prominant timescale in the stimulus
(what we refer to as the base period) is equated to the predefined
value of ``base_tau``, which we set it to be 1. To do so, we first
synthesize a temporary stimulus of size 100k to estimate the base
period and time scaling factor. Then, we resynthesize the stimulus
for the required number of timesteps with an adjusted sampling rate
that allocates :math:`\tau_{base}/dt` samples within the base period.
Then, the target is generated based on the synthesized stimulus, and
finally, everything is saved on the disk.
.. note:: A pure constant (DC) input carries no characteristic
timescale. For such input, the temporal characterization
and the time rescaling described above are bypassed
(``dt_scaler`` remains 1), since no temporal scaling is
necessary for a constant input.
"""
# skip if already created
created_already = os.path.exists(osjoin(self.res_path, 't.npy')) and \
os.path.exists(osjoin(self.res_path, 'u.npy')) and \
os.path.exists(osjoin(self.res_path, 'y.npy'))
if not created_already:
print('generating stim and target anew')
if tasks.needs_synthesis(self.task.name):
# A pure constant (DC) input has no characteristic timescale,
# so its temporal characterization (base-lag estimation) and
# the subsequent time rescaling are bypassed. No rescaling
# means dt_scaler=1, i.e., the synthesis step remains `h`.
if stimulus.is_constant_stim(self.stim_params):
print(' constant (DC) input: skipping temporal '
'characterization and time rescaling')
base_lags = [1. for _ in self.stim_params] # placeholder (no timescale)
base = 1. # placeholder (nominal base_tau)
dt_scaler = 1. # no temporal rescaling
else:
print(' analyzing stim')
# here we first create a stimulus of size 100k and use it
# to compute the main timescale of that input
# shape : (Σ ndim x ntimes)
u = self.task.make_stim(n_timesteps=100000, **self.stim_params)
# The ``base_lags`` are the porminant timescale of each component
# of u. ``base`` is their average.
base_lags, base = utils.find_base_lags(u)
# from which the we can estimate how faster or slower we
# should sample to ensure about 100 samples per each period.
dt_scaler = base / self.bs
print(f' scaling sampling rate by {dt_scaler}')
# and now we actually generate the stim
# shape : (Σ ndim x ntimes)
print(' synthesizing stim')
u = self.task.make_stim(n_timesteps=self.n_timesteps,
dt_scaler=dt_scaler,
**self.stim_params)
# and we pretend that one period equals base_tau unit of
# time. Consequently, the dt is exactly what we set up
# above in the initialziation, i.e. base_tau/self.bs
t = np.linspace(0, self.n_timesteps*self.dt,
self.n_timesteps, endpoint= False)
print(' generating target')
y = self.task.make_target(u)
# if input is multidimensional, so is the task. But we don't
# really care about the different components. In fact, all components
# are treated as distinct subtasks. So, we flatten the target
print(' flattening target')
self.task.subtasks, y = self.task.flatten(y)
ytype = self.stype
else:
# To be extended for other task types ...
pass
# t, u, y = self.task.make_stim_target(self.n_timesteps, self.dt)
# ytype = np.uint16
# t, base_lag, base_lags = utils.rescale_time(t, u, target_base_freq=1./self.base_tau)
# np.save(osjoin(self.res_path, 't.npy'), t.astype(self.stype))
# np.save(osjoin(self.res_path, 'u.npy'), u.astype(self.stype))
# np.save(osjoin(self.res_path, 'y.npy'), y.astype(int))
# save on disk
np.save(osjoin(self.res_path, 't.npy'), t)
np.save(osjoin(self.res_path, 'u.npy'), u.astype(self.stype))
np.save(osjoin(self.res_path, 'y.npy'), y.astype(ytype))
# and create the stimulus statistic metadata
stim_stats = dict(dim=u.shape[0], base_lag=base, base_lags=base_lags)
# and clean up
del t, u, y
else:
# just load the previous stim metadata
stim_stats = dict(dim=self.prev_md['dim'],
base_lag=self.prev_md['base_lag'],
base_lags=self.prev_md['base_lags'])
# flattening subtasks if necessary
if tasks.needs_synthesis(self.task.name):
self.task.subtasks = self.task.flatten_subtasks(n_dims=self.prev_md['dim'])
# finally update metadata
self.metadata.update(stim_stats)
[docs]
def compute_complexity(self):
"""Computes subtasks' complexity and saves the values. Complexity is defined as
the cosine similarity between the a particular subtask and a reference subtask.
The reference subtask is the task with minimal complexity. Thus, it is a relative
measure, between zero (very easy task) to 1 (maximally difficult).
The computed complexity is saved on the disk as a csv file.
.. note:: We used the relay task as the reference task since it needs no shift
to the past or future, or non-linear processing.
"""
complexity_path = osjoin(self.res_path, 'complexities.csv')
if not os.path.exists(complexity_path):
print('computing complexity')
_, u = utils.load_inputs(self.res_path, load_noise=False)
y = utils.load_target(self.res_path)
utils.compute_complexity(u, y, self.task.subtasks, complexity_path)
[docs]
def compute_udf(self, ncpu=-1):
"""Computes the input-based baseline scores (udf) and saves ``udf.csv``.
The udf baseline fits a linear decoder directly from the (compound)
input ``U`` to the task target ``Y`` and records, per subtask, the train
and test R2 scores. The resulting dataframe (columns ``dim, delta, deg,
trn, tst``) is saved in the result directory as ``udf.csv``.
The scores quantify how much of each target subtask is *linearly
reachable from the input alone* -- i.e. how much the recurrent network
can add on top of the input -- and serve as the baseline against which
both the network readouts and the (cosine-based) complexities are
compared (see :meth:`plot_udf_vs_complexity`).
The udf baseline is *always* fitted via OLS.
Parameters
----------
ncpu: int, optional
Number of CPU-cores used for the readout. Default is -1 (all
available cores).
.. note:: Unlike the complexities (``complexities.csv``), the udf scores
are absolute (not relative) and are measured with the very same
decoder/readout that is used for the network states.
"""
udf_path = osjoin(self.res_path, 'udf.csv')
if os.path.exists(udf_path):
return
print('computing udf (input-to-target baseline)')
# the udf baseline is always fitted via OLS
optim_params = {'alg': 'ridge', 'reg': 0, 'fit_bias': True}
decoder = fitting.Decoder(optim_params=optim_params, ncpu=ncpu)
_, U = utils.load_inputs(self.res_path, load_noise=False)
Y = utils.load_target(self.res_path)
s_trn, s_tst = self.Strn, self.Stst # train and test indices
print(U.shape)
# we use all samples in the train to get maximal possible score
U_trn, Y_trn = U[:, s_trn], Y[:, s_trn]
U_tst, Y_tst = U[:, s_tst], Y[:, s_tst]
decoder.fit(U_trn, Y_trn)
uscore_trn = decoder.score(decoder.predict(U_trn), Y_trn)
uscore_tst = decoder.score(decoder.predict(U_tst), Y_tst)
udf = self.task.get_score_template()
udf['trn'] = uscore_trn
udf['tst'] = uscore_tst
udf.to_csv(udf_path)
del Y
[docs]
def synthesize_noise(self):
"""Generates and saves the white noise time series for the network."""
path_xi = osjoin(self.res_path, f'xi.npy')
# if not os.path.exists(path_xi):
print('generating noise anew')
xi = self.rng.normal(0, 1, size=(self.net_params['N'], self.n_timesteps))
np.save(path_xi, xi.astype(self.stype))
[docs]
def synthesize_topology(self):
"""Creates the network topology and weights."""
print('Generating network topology and weights anew.')
path_A_adj = osjoin(self.res_path, f'A_adj.npz')
path_w = osjoin(self.res_path, f'w.npy')
path_w_ff = osjoin(self.res_path, f'w_ff.npy')
path_w_fb = osjoin(self.res_path, f'w_fb.npy')
# done = os.path.exists(path_A_adj) and (os.path.exists(path_w)) and \
# os.path.exists(path_w_fb) and (os.path.exists(path_w_ff))
# if not done:
# if self.needs_new_realization(): # will override previous ones
# N_syn = int(round(N**2 * self.net_params['p']))
# A_adj = networks.create_topology(rng = self.rng, N=N, N_syn=N_syn,
# allow_autapse = self.net_params['autapse'],
# topology=self.net_params['topology'])
w_ff, w, w_fb, A_adj = networks.intialize_weights(
N_input=self.metadata['dim'],
N = self.net_params['N'],
N_output = 0,
p = self.net_params['p'],
f = self.net_params['f'],
mu_e=self.net_params['mue'],
sigma_0=self.net_params['sig0'],
return_adj=True,
rng = self.rng)
np.save(path_w_ff, w_ff)
np.save(path_w, w)
np.save(path_w_fb, w_fb)
save_npz(path_A_adj, A_adj)
[docs]
def plot_connectivity(self, rlz_id=0):
"""Plots the connectivity matrix and the recurrent weight distribution
of the current realization into the result path."""
from hetero.viz import plot_connectivity as plot_conn
plot_conn(self.res_path, rlz_id=rlz_id)
[docs]
def plot_taus(self, rlz_id=0):
"""Plots the distribution of the neuronal time constants of all networks
for the current realization into the result path."""
from hetero.viz import plot_tau_dist as plot_taus
plot_taus(self.res_path, rlz_id=rlz_id)
[docs]
def plot_scores(self, rlz_id=0):
"""Plots the heterogeneous-vs-homogeneous score of every task for the
given realization after it has been benchmarked."""
from hetero.viz import plot_score_score
plot_score_score(self.res_path, rlz_id=rlz_id)
[docs]
def plot_angle_vs_score(self, rlz_id=None):
"""Plots, for every task, the minimal angle between the network state
and the task against the task's benchmark score, grouped by the
mean-timescale (``n``) and heterogeneity (``h``) levels of the
networks (see :func:`hetero.viz.plot_angle_vs_score`).
Meant to be called once all realizations are benchmarked; when
``rlz_id`` is None all realizations found are plotted."""
from hetero.viz import plot_angle_vs_score as _plot_angle_vs_score
try:
_plot_angle_vs_score(self.res_path, rlz_id=rlz_id)
except Exception as e:
print(f'plot_angle_vs_score failed: {e}')
[docs]
def plot_overlap_cdf(self, rlz_id=None):
"""Plots the cumulative distribution of the minimal angles of every
task, grouped by the mean-timescale level (``n``) of the networks, with
the heterogeneity levels (``h``) distinguished by color (see
:func:`hetero.viz.plot_overlap_cdf`).
Meant to be called once all networks are benchmarked; when ``rlz_id``
is None all realizations found are plotted."""
from hetero.viz import plot_overlap_cdf as _plot_overlap_cdf
try:
_plot_overlap_cdf(self.res_path, rlz_id=rlz_id)
except Exception as e:
print(f'plot_overlap_cdf failed: {e}')
[docs]
def plot_svals_erank(self, rlz_id=None):
"""Plots the singular-value spectrum of every network state and the
cumulative distribution of its effective rank, grouped by the
mean-timescale level (``n``) of the networks with the heterogeneity
levels (``h``) distinguished by color (see
:func:`hetero.viz.plot_svals_erank`).
Meant to be called once all networks are benchmarked; when ``rlz_id``
is None all realizations found are plotted."""
from hetero.viz import plot_svals_erank as _plot_svals_erank
try:
_plot_svals_erank(self.res_path, rlz_id=rlz_id)
except Exception as e:
print(f'plot_svals_erank failed: {e}')
[docs]
def plot_uvx(self, u, v, x, net_id, state_name):
"""Plots the input ``u``, the raw states ``v`` (before the readout
non-linearity), and the activated states ``x`` of the network
``net_id`` into its subfolder of the result path. The file name
embeds ``state_name`` so distinct dynamics/realizations (e.g. ``LI-
...-rlz000`` vs ``LIF-...-rlz000``) do not overwrite each other."""
from hetero.viz import plot_uvx as plot_uvx_viz
plot_uvx_viz(u, v, x, osjoin(self.res_path, net_id), state_name=state_name)
[docs]
def plot_complexity_map(self):
"""Plots the complexity of every subtask, per input dimension, from the
*previously calculated* complexities (``complexities.csv``) and saves the
figure(s) into the result directory.
The complexities are not recomputed here: they were already estimated in
``compute_complexity`` during initialization and persisted on the disk.
"""
from hetero.viz import plot_complexity_map as _plot_complexity_map
try:
_plot_complexity_map(self.res_path)
except Exception as e:
print(f'plot_complexity_map failed: {e}')
[docs]
def plot_udf_vs_complexity(self):
"""Plots the input-based baseline scores (``udf.csv``) against the
precomputed complexities (``complexities.csv``), per subtask, and saves
the figure(s) into the result directory.
Neither of the two quantities is recomputed here: the complexities were
estimated in ``compute_complexity`` and the udf scores in
``compute_udf``, both during initialization.
"""
from hetero.viz import plot_udf_vs_complexity as _plot_udf_vs_complexity
try:
_plot_udf_vs_complexity(self.res_path)
except Exception as e:
print(f'plot_udf_vs_complexity failed: {e}')
[docs]
def make_task_basis(self, Y, var_frac=0.01):
"""Prepares the task targets once for the minimal-angle analysis.
The targets ``Y`` (``n_tasks x T``) are *fixed* across all networks,
whereas the network state ``X`` varies from one network to the next.
To avoid re-orthogonalizing the targets for every network in
:meth:`benchmark`, this method is called a single time (before the
per-network loop) and its output is fed to :meth:`compute_minimal_angles`.
Each target is first *centered* (its temporal mean over the slice is
removed) and then row-normalized. Since each task spans a 1-D direction
in the trajectory space, its orthonormal basis -- the very object that
``scipy.linalg.orth`` (or, per task, the ``QB`` step of
``scipy.linalg.subspace_angles``) would produce -- is simply its
unit-norm direction. Those unit directions are computed here, once:
this is the "already orthogonalized Y" that the per-network angle
computation consumes.
Centering removes each task's constant (DC) part, which any linear
readout with a bias can fit regardless of the state. The remaining
direction is the task's *fluctuation*, and the angle measures whether
that fluctuation lies in the state's fluctuation subspace (see
:meth:`compute_minimal_angles`). Tasks with no fluctuation at all (a
constant target, or one collapsed by slicing) have zero centered norm;
they are marked separately and assigned a zero angle.
The task *subspace* is truncated exactly like the state side (see
:meth:`compute_minimal_angles`): the centered task matrix is SVD'ed
and only the modes whose relative contribution to the total task
variance is at least ``var_frac`` are retained. Each task direction is
then rebuilt from its projection onto that retained (numerically
robust) subspace and renormalized inside it. A task whose projected
fluctuation keeps less than ``var_frac`` of its own variance is as
numerically fragile as a constant target -- its direction is dominated
by roundoff -- so it is marked ``zero`` and assigned a zero angle (its
fluctuation is treated as trivially representable, the same convention
as constant tasks).
Parameters
----------
Y : ndarray, shape (n_tasks, T)
The (NaN-free, already sliced) target matrix.
var_frac : float, optional
Minimum relative variance contribution (``s_k^2 / sum s_j^2``) for
a task-subspace mode to be retained. Default ``0.01`` (1%).
Returns
-------
dict
With the keys:
- ``Yn``: ndarray ``(n_tasks, T)``, the projected, row-wise
normalized targets (unit-norm task directions in the retained
task subspace).
- ``zero``: ndarray ``(n_tasks,)`` boolean mask of tasks with no
robust fluctuation (constant after centering, or fluctuation
nearly confined to the truncated task-subspace modes).
"""
Yc = np.asarray(Y, dtype=np.float64) # shape: (n_tasks, T)
Yc = Yc - Yc.mean(axis=1, keepdims=True) # remove each task's DC part
norms = np.linalg.norm(Yc, axis=1)
zero = norms == 0
# joint truncation of the task subspace, the same variance-fraction
# recipe as the state side (see compute_minimal_angles): keep the modes
# of the centered task matrix with relative variance >= var_frac; the
# rows of vhy are an orthonormal basis of the retained subspace (rows
# of Yc live in column space == row space of Yc.T, spanned by vhy).
_, sy, vhy = svd(Yc, full_matrices=False) # vhy: (min(m, T) x T)
keep = _variance_keep(sy, var_frac)
if keep.any():
Qy = vhy[keep] # (k x T) robust task basis
Yproj = Yc @ Qy.T # (m x k) coefficients
kept_norm = np.linalg.norm(Yproj, axis=1)
# tasks that lose almost all their fluctuation to the truncated
# modes are numerically fragile directions; route them to the
# constant-task convention (zero angle) instead of the roundoff
# that would survive repeated normalization.
zero = zero | ((kept_norm ** 2) < var_frac * (norms ** 2))
with np.errstate(divide='ignore', invalid='ignore'):
Yn = Yproj @ Qy / kept_norm[:, None] # (m x T) unit directions
else:
kept_norm = norms
Yn = Yc / norms[:, None]
Yn[zero] = 0.0 # guard the 0/0 division
return dict(Yn=Yn, zero=zero)
@staticmethod
def _effective_rank(svals, n_neurons):
"""Effective rank of an activity spectrum -- the normalized
participation ratio.
Given the singular values ``svals`` of a matrix, the effective rank
is the participation ratio of the spectrum, normalized by the number
of neurons:
erank = (sum_k s_k**2)**2 / (n_neurons * sum_k s_k**4)
It equals ``1/n_neurons`` when a single mode dominates the activity,
``1`` when all ``n_neurons`` modes carry equal variance (each neuron's
activity is an independent, equally-strong direction), and interpolates
smoothly in between -- unlike the numerical rank, it does not depend on
an arbitrary threshold. The all-zero (no fluctuation) spectrum is
assigned an effective rank of ``0``.
"""
s2 = svals ** 2
total = s2.sum()
if total == 0:
return 0.0
return float(total ** 2 / (n_neurons * (s2 ** 2).sum()))
[docs]
def compute_minimal_angles(self, net_id, state_name, X,
slices, task_bases, var_frac=0.01):
"""Computes the minimal angle between the network state ``x`` and every
task, and saves the per-task angles into the network's subfolder.
For each task (subtask) ``i`` the *minimal angle* is the angle between
the centered task direction ``y_i`` and the linear subspace spanned by
the centered neuron activities ``X`` -- equivalently, the angle between
``y_i`` and its projection onto the span of ``X``:
cos(angle_i) = ||P_{span(X)} y_i|| / ||y_i||
``y_i`` is centered (mean removed) and unit-normalized by
:meth:`make_task_basis`, and ``X`` is centered here, per slice, before
its basis is built (each neuron's temporal mean over the slice is
subtracted). Both sides therefore have the constant (DC/bias) component
projected out, so the angle measures the reachability of the task's
*fluctuation* by the state's fluctuation -- equivalently, the best
achievable |correlation| of a linear readout (bias or not) with the
task. A zero angle means the task is perfectly reachable by a linear
readout of the state; an angle of pi/2 means it is orthogonal
(unreachable).
The results are saved in the *same* network folder and with the *same*
naming scheme used for the states. Alongside the per-task angles, the
singular-value spectrum of the (centered, per-slice) activity and its
effective rank are saved. Three files are written, all in the network's
subfolder:
- ``angle_{state_name}.csv`` -- one row per subtask (the columns of
``task.get_score_template()``) followed by one ``angle_{trn|tst}``
column per slice.
- ``svals_{state_name}.csv`` -- one row per singular value of the
centered activity, with one ``s_{trn|tst}`` column per slice
(singular values returned by the SVD, in descending order; the
shorter slice's column is NaN-padded when the train and test slices
differ in length).
- ``erank_{state_name}.csv`` -- one row per slice (``mode``) with the
slice's effective rank, computed from the largest singular values
via :meth:`_effective_rank` (the normalized participation ratio,
rank, see there for the definition).
Parameters
----------
net_id : str
The network ID (its subfolder in the result path).
state_name : str
The state name, embedded in the output file names.
X : ndarray, shape (N, T)
The (already loaded) activated state of this network, uncentered.
slices : tuple of slice
The train and test slices (``(s_trn, s_tst)``). They delimit the
NaN-free portions of the trajectory.
task_bases : tuple of dict
The precomputed task basis per slice (``(prep_trn, prep_tst)``),
as returned by :meth:`make_task_basis` and prepared exactly once in
:meth:`benchmark`. Passing them here avoids re-orthogonalizing the
fixed targets for every network.
var_frac : float, optional
Minimum relative variance contribution (``s_k^2 / sum s_j^2``) for
a singular-value mode of the (centered) activity to be retained in
the state subspace. Modes below this fraction are numerically
non-robust and are truncated away; the task side is truncated with
the same threshold in :meth:`make_task_basis`. Default ``0.01``
(1%).
Notes
-----
- The raw (uncentered) state ``X`` is left untouched and used by the
decoder as-is; only a per-slice centered copy enters the angle
computation, so the decoder's own bias term is not perturbed.
- The per-task angle is computed the way ``scipy.linalg.subspace_angles``
does, specialized to the 1-D task directions: the cosine is the norm
of the projected direction (the SVD singular value of the
cross-product matrix), the sine is the norm of the projection
residual, and the angle is taken with :func:`numpy.arcsin` where
``cos^2 >= 0.5`` (small angles) and :func:`numpy.arccos` otherwise
(large angles). A bare ``arccos(cos)`` loses precision when the
cosine is very small (angle near pi/2), which this hybrid scheme
avoids.
- A *single* SVD of the centered state drives both the angle subspace
and the saved spectrum: its left singular vectors span the state
subspace and its singular values ``s`` are the spectrum of the
activity itself. The retained columns are the modes whose relative
contribution to the total variance is at least ``var_frac``
(``s_k^2 / sum s_j^2 >= var_frac``, default 1%) -- a fixed,
scale-free cutoff that drops the numerically non-robust tail of the
spectrum, replacing the ``eps * max(N, T) * s_max`` tolerance of
``scipy.linalg.orth`` (near-threshold modes whose orientation is
dominated by roundoff are what made the old cutoff fragile). The
(full, untruncated) spectrum/effective rank therefore still describe
the *centered* activity per slice (see :meth:`_effective_rank`).
- No file is written if all three ``angle_/svals_/erank_`` files
already exist for this network, and each missing file is written
individually, so reruns, resumed benchmarks, and replots are safe
(a network whose angles were computed before this feature simply
gets its ``svals_``/``erank_`` files backfilled).
- ``X`` is expected to be already loaded; this method performs no disk
I/O on the state itself.
- Angles are reported in radians, within [0, pi/2].
"""
angle_path = osjoin(self.res_path, net_id, f'angle_{state_name}.csv')
svals_path = osjoin(self.res_path, net_id, f'svals_{state_name}.csv')
erank_path = osjoin(self.res_path, net_id, f'erank_{state_name}.csv')
if (os.path.exists(angle_path) and os.path.exists(svals_path)
and os.path.exists(erank_path)):
return
df = self.task.get_score_template()
svals = {} # mode -> singular values of the centered activity
erank = {} # mode -> effective rank of the centered activity
for mode, slc, basis in zip(('trn', 'tst'), slices, task_bases):
Xc = X[:, slc].astype(np.float64) # (N, T)
Xc = Xc - Xc.mean(axis=1, keepdims=True) # center per neuron, per slice
# Single SVD of the centered state feeds both outputs: the retained
# left singular vectors u span the state subspace, while the
# singular values s are the spectrum of the activity itself, saved
# below together with the effective rank. The columns are kept by
# their relative contribution to the total variance (>= var_frac,
# default 1%) -- a fixed, scale-free cutoff that drops the
# numerically non-robust tail of the spectrum (see _variance_keep).
u, s, vh = svd(Xc.T, full_matrices=False) # (T x n), (n,), (n x N); n = min(N, T)
Qx = u[:, _variance_keep(s, var_frac)] # (T x rank) state-space basis
Yn = basis['Yn'] # (m x T) unit task directions (precomputed QB)
# Mirror scipy.linalg.subspace_angles, specialized to the 1-D task
# directions. The cosine comes from the norm of the projected
# direction (the cross-products), the sine from the norm of the
# projection residual, and the angle is taken with arcsin where
# cos^2 >= 0.5 (small angles) and with arccos otherwise (large
# angles). Unlike a bare arccos(cos), this stays accurate when the
# cosine is very small (angle close to pi/2).
P = Qx.T @ Yn.T # (n x m) cross-products (QA_H_QB)
cos = np.linalg.norm(P, axis=0) # svdvals of each (n x 1): |cos(angle)|
R = Yn.T - Qx @ P # (T x m) projection residual (QB - QA QA_H_QB)
sin = np.linalg.norm(R, axis=0) # svdvals of each (T x 1): |sin(angle)|
angles = np.where(cos ** 2 >= 0.5,
np.arcsin(np.clip(sin, -1.0, 1.0)),
np.arccos(np.clip(cos, -1.0, 1.0)))
angles[basis['zero']] = 0.0 # zero tasks are trivially representable
df[f'angle_{mode}'] = angles
# spectrum and effective rank of this slice's activity
svals[mode] = s
erank[mode] = self._effective_rank(s, X.shape[0])
# write the three CSVs, each only when missing (resume-safe)
if not os.path.exists(angle_path):
df.to_csv(angle_path, index=False)
if not os.path.exists(svals_path):
pd.DataFrame({f's_{mode}': pd.Series(v)
for mode, v in svals.items()}).to_csv(svals_path, index=False)
if not os.path.exists(erank_path):
pd.DataFrame({'mode': list(erank),
'erank': list(erank.values())}).to_csv(erank_path, index=False)
[docs]
def get_result_path(self):
"""Creates the result path."""
return osjoin(self.root, self.name)
[docs]
def prepare_dirs(self):
"""Creates one subdirectory for each network in the result path."""
if not os.path.exists(self.res_path):
# to avoid workers' race condition we pass exist_ok
os.makedirs(self.res_path, exist_ok=True)
for net_id in self.net_ids:
if not os.path.exists(osjoin(self.res_path, net_id)):
# to avoid workers' race condition we pass exist_ok
os.makedirs(osjoin(self.res_path, net_id), exist_ok=True)
[docs]
def run_dynamic(self, dyn, dyn_params,
Ju, J, Jn,
rlz=0, purge_old=True,
ncpu=-1, solver='euler',
verbose=True, **kwargs):
"""
Integrates the neural dynamic by creating an object that encapsulates
the dynamical rule ``dyn``, parameters ``dyn_params`` and provides a
method to integrate the dynamic.
Parameters
----------
dyn : str
The name of the neural dynamic to be integrated.
dyn_params: dict
The parameters of the dynamic
Ju: float
Input (feedforward) gain
J: float
Recurrent gain
Jn: float
Noise intensity
ncpu: int
Number of CPU-cores used for parallelizing integration. Default
is -1, which uses all available cores. For large networks, this
might fill up the machine RAM, and reduce the overall speed. If
so, reduce it to half of your available cores or less.
solver: str
Numberical integration method. For simulations with non-zero noise
the only valid option is ``euler`` that implements Euler-Maruyama
method.
Returns
-------
state_name: str
The state name
.. note:: For spiking networks, we leverage the Brian2 library, whereas the
rate-based neurons are implemented in pure Python.
.. warning:: The rate-based dynamic is parallelized via joblib. Depending
on theoperating system, one might need to change the parallelizing backend.
"""
# the following contain no info of the realizaiton in their name
state_name = self.generate_state_name(dyn, Ju, J, Jn, solver) # acts as a template
w_ff, w_rec, _ = utils.load_weights(self.res_path) # overridden on each rlz
_, u, xi = utils.load_inputs(self.res_path, load_noise=True) # xi is overridden on each rlz
rlz_suffix = self.rlz_suffix(rlz)
state_rlz = state_name+rlz_suffix
print(f'Running dynamic {state_name} | realization {rlz} --> state_rlz: {state_rlz}')
namespace = {
# 'net_ids': self.net_ids,
# 'taus': utils.load_taus(),
'N': self.metadata['net']['N'],
'A_adj': utils.load_topology(self.res_path),
'delay': self.metadata['net']['delay'],
'wu': w_ff, 'w': w_rec,
'J': J, 'Ju': Ju, 'Jn':Jn,
'u': u, 'xi': xi, # time-dependent
'dt': self.metadata['dt'],
'n_timesteps': self.metadata['n_timesteps'],
'y0': utils.load_ic(self.res_path),
'verbose': verbose, # report integration progress
}
dyn_class = LIF if dyn == 'LIF' else LI
dyner = dyn_class(namespace=namespace,
sim_path=self.res_path,
dyn_params=dyn_params)
# limit the ncpus to num net_ids and no more (matters for HPCs)
if ncpu == -1:
ncpu = min(NCORES, len(self.net_ids))
dyner.run(Ju=Ju, J=J, Jn=Jn, ncpu=ncpu, solver=solver,
state_name=state_rlz, stype=self.stype,
verbose=verbose)
if purge_old and rlz>0:
# print(self.num_realized)
# print(state_name)
for old_rlz in range(rlz):
old_state = state_rlz.replace(f'rlz{rlz:03d}', f'rlz{old_rlz:03d}')
for net_id in self.net_ids:
old_file_path = self.get_state_path(net_id, old_state, dyn=dyn)
utils.purge(old_file_path)
return state_name # the temaplate without rlz info
[docs]
def benchmark(self, state_name, readout_params, optim_params, n_cpu,
skip_uvx_plot=False, skip_svd=False, var_frac=0.01):
"""
Quantifies the performance (r2 score) of a network whose state
specified by ``state_name`` using the readout parameters
``readout_params`` for all subtasks.
The results are saved in each network subdirectory as csv files
with a name composed of the ``state_name`` and a uniquely generated
name for the provided parameters ``readout_params``. The train and
test performances are distinguished by the leading ``trn_`` or ``tst_``
strings in file names.
Parameters
----------
state_name : str
The state name to be used for readout.
readout_params: dict
Readout parameters.
optim_params: dict
Optimization parameters.
ncpu: int
Number of CPU-cores used for parallelizing the readout. Default
is -1, which uses all available cores. Only works for ridge and
OLS decoders.
Returns
-------
state_name: str
The state name
.. note:: For spiking networks, we leverage the Brian2 library, whereas the
rate-based neurons are implemented in pure Python.
.. note:: In addition to the readout scores, the minimal angle between
the network state ``x`` and every task -- defined as the angle between
each task direction and the linear subspace spanned by the state --
is computed once the state is loaded, and saved per network as
``angle_{state_name}.csv``, together with the singular-value
spectrum (``svals_{state_name}.csv``) and the effective rank
(``erank_{state_name}.csv``) of the centered activity (see
:meth:`compute_minimal_angles`).
.. warning:: The rate-based dynamic is parallelized via joblib. Depending
on theoperating system, one might need to change the parallelizing backend.
"""
print('starting benchmark')
# remember the optimization parameters
self.optim_params = optim_params
# setting up the decorder
decoder = fitting.Decoder(optim_params=optim_params, ncpu=n_cpu)
_, U = utils.load_inputs(self.res_path, load_noise=False)
# scan for those that are not yet fitted. This is better than
# loading all and then ignorning, because the states are HUGE!
net_id_dict = self.scan_unbenchmarked(
state_name=state_name,
readout_params=readout_params,
optim_params=optim_params,
verbose=True,
)
# exit if all are fitted
if len(net_id_dict.values())==0:
return 1
# loading targets
print('loading targets')
Y = utils.load_target(self.res_path)
# retrieving the train/test sets
s_trn, s_tst = self.Strn, self.Stst # train and test indices
# ------------------------------------------------------------------
# Minimal-angle analysis: prepare the task subspaces *once*.
# The targets are fixed across all networks (only X varies), so
# orthogonalizing them here -- rather than for every network -- avoids
# redundant work. The per-network angles are computed in a separate
# code block below, right after the scores, reusing the loaded X.
# ------------------------------------------------------------------
if not skip_svd:
task_bases = dict(
trn=self.make_task_basis(Y[:, s_trn], var_frac=var_frac),
tst=self.make_task_basis(Y[:, s_tst], var_frac=var_frac),
)
# a distinct rng with a fixed hardcoded seed for shuffling training
# data. The reason for the divergence from self.rng is that, sometimes
# the integration phase is done but the benchmark is skipped. Resuming
# to the benchmark, however, may skip all the rng calls during network
# synthesis. So, effectively, a different shuffling might get used.
# So we fix the shuffling scheme here.
bm_rng = np.random.default_rng(12345) # seed 1234 + 5 :)
for net_id, net_trials in net_id_dict.items():
# load state
tic = time.time()
V, X = dataloader.load_states(
sim_path = self.res_path, net_id = net_id, state_name = state_name,
activation_params = readout_params['readout_activation'],
n_neuron = self.metadata['net']['N'],
n_timesteps= self.metadata['n_timesteps'],
)
print(f'state {net_id} loaded in {time.time() - tic:.1f} seconds.')
print(f'Set size V:{V.shape}, X:{X.shape}, Y:{Y.shape}')
if not skip_uvx_plot:
self.plot_uvx(U, V, X, net_id, state_name)
# spare the test set
X_tst, Y_tst = X[:, s_tst], Y[:, s_tst]
# X_tst = X[:, self.Ntrn:]
# Y_tst = Y[:, self.Ntrn:]
# shuffle the train set
# We don't use the sklearn's KFold because it doesn't accept our rng
# trn_idxs = np.arange(self.Ntrn)
trn_idxs = np.arange(s_trn.start, s_trn.stop, 1)
bm_rng.shuffle(trn_idxs)
# kf = KFold(n_splits=3, shuffle=False, random_state=self.random_state)
# for i, (trial_idxs, _) in enumerate(kf.split(X_trn)):
df_trn = self.task.get_score_template()
df_tst = self.task.get_score_template()
# print('ho ho ho 0')
# print(fit_params)
# split into different trials and fit
for ntrial in net_trials:
path_trn_scr, path_tst_scr = self.get_score_paths(
net_id=net_id,
state_name=state_name,
readout_params=readout_params,
optim_params=optim_params,
ntrial=ntrial,
)
tic = time.time()
print(f' Train split into {ntrial} trials.')
for i, trn_idx in enumerate(np.array_split(trn_idxs, ntrial)):
trial_tic = time.time()
X_trn, Y_trn = X[:, trn_idx], Y[:, trn_idx]
# reporting shapes
if i==0:
print(f' train input shape: {X_trn.shape}, train output shape: {Y_trn.shape}')
if decoder.alg == 'logreg': # for classification, subtasks must be done individually
# print('pre shapes')
# print(X_trn.shape, Y_trn.shape, len(pd.unique(Y_trn[2,:])))
# if max_norm:
# norm = np.linalg.norm(X_trn, axis=0)
# xx, yy = [], []
# for sample_id in pd.unique(Y_trn[2,:]):
# idxs = Y_trn[2,:] == sample_id
# xx.append(X_trn[:, idxs][:, np.argmax(norm[idxs])])
# yy.append(Y_trn[:, idxs][:,0])
# elif accum:
# xx, yy = [], []
# sample_ids = pd.unique(Y_trn[2,:])
# for q, sample_id in enumerate(sample_ids):
# print(sample_id, f'{q}/{len(sample_ids)}')
# idxs = Y_trn[2,:] == sample_id
# xx.append(X_trn[:, idxs].sum(axis=1))
# yy.append(Y_trn[:, idxs][:,0])
# else:
# xx = X_trn
# yy = Y_trn
# X_trn = np.array(xx).T
# Y_trn = np.array(yy).T
# print('post shapes')
trn_scrs, tst_scrs = [], []
for subtask_id in range(len(self.task.subtasks)):
decoder.fit(X_trn, Y_trn[subtask_id, :])
trn_scrs.append(decoder.score(decoder.predict(X_trn), Y_trn[subtask_id, :]))
tst_scrs.append(decoder.score(decoder.predict(X_tst), Y_tst[subtask_id, :]))
df_trn[i] = trn_scrs
df_tst[i] = tst_scrs
else: # for regression, subtasks can be done all at once
decoder.fit(X_trn, Y_trn)
df_trn[i] = decoder.score( decoder.predict(X_trn), Y_trn )
df_tst[i] = decoder.score( decoder.predict(X_tst), Y_tst )
# w = fitting.fit(X_trn, Y_trn, method='sk_ridge', alpha=1e-4)
# trn_scrs.append(fitting.score(w@X_trn, Y_trn))
# tst_scrs.append(fitting.score(w@X_tst, Y_tst))
# writing the result of each fold into the dataframe
# df_trn[i] = trn_scrs
# df_tst[i] = tst_scrs
del X_trn, Y_trn#, xx, yy
print(f' trial number {i+1} completed in {time.time() - trial_tic:.1f} seconds.')
# # makes scores
# df_trn = pd.DataFrame(np.array(trn_scrs).T)
# df_tst = pd.DataFrame(np.array(tst_scrs).T)
# pd.concat([df, df_trn], axis=1).to_csv(path_trn_scr, index=False)
# pd.concat([df, df_tst], axis=1).to_csv(path_tst_scr, index=False)
df_trn.to_csv(path_trn_scr, index=False)
df_tst.to_csv(path_tst_scr, index=False)
print(f' {ntrial} splits of network {net_id} results fitted in {time.time() - tic:.1f} seconds.')
print()
# ------------------------------------------------------------------
# Minimal angle between the network state x and the tasks y.
# Kept in its own block, after the fitting/scores: x is already in
# memory from the load above (loaded only once), and the task
# subspaces were prepared only once before the loop.
# ------------------------------------------------------------------
if not skip_svd:
self.compute_minimal_angles(
net_id=net_id,
state_name=state_name,
X=X,
slices=(s_trn, s_tst),
task_bases=(task_bases['trn'], task_bases['tst']),
var_frac=var_frac,
)
[docs]
def get_n_integrated(self, n_realizations, state_name, dyn):
"""
returns the number of fully integrated realizations
"""
integrated = 0
for rlz in range(n_realizations):
state_name_rlz = state_name + self.rlz_suffix(rlz)
print(state_name_rlz)
# either there's nothing integrated, or some integrated and prunned. So
# we look for the last one that is integrated
if self.state_files_exist(state_name_rlz, dyn=dyn) == 'partial':
integrated = rlz # try to integrate the partial one
break
elif self.state_files_exist(state_name_rlz, dyn=dyn) == True:
integrated = rlz + 1 # one fully integrated case
break
else:
continue
return integrated
[docs]
def get_n_benchmarked(self, n_realizations, state_name, optim_params, readout_params):
"""
returns the number of realizations that are fully benchmarked
"""
benchmarked = 0
for rlz in range(n_realizations):
state_name_rlz = state_name + self.rlz_suffix(rlz)
score_name = self.generate_score_template(
state_name=state_name_rlz,
readout_params=readout_params,
optim_params=optim_params,
)
if not self.score_files_exist(score_name):
break
benchmarked+=1
return benchmarked
[docs]
def setup_runner_log(self, n_realizations, dyn, solver,
gain_params, optim_params, readout_params):
"""a simple running log for a given state, indicating the number of
realization that are integrated and benchmarked."""
running_state_name = self.generate_state_name(dyn=dyn, solver=solver, **gain_params, )
n_integrated = self.get_n_integrated(n_realizations, running_state_name, dyn)
n_benchmarked = self.get_n_benchmarked(n_realizations, running_state_name, optim_params, readout_params)
self.rlog = {
'required_relizations': n_realizations,
'integrated': n_integrated,
'benchmarked': n_benchmarked,
}
self.update_runner(running_state_name)
[docs]
def update_runner(self, running_state_name):
with open(osjoin(self.res_path, f'log_{running_state_name}.json'), 'w') as file:
json.dump(dict(self.rlog), file, indent=2, cls=utils.NpEncoder)
[docs]
def replot_all(self, n_realizations, dyn, dyn_params, gain_params,
readout_params, optim_params, solver, ncpu):
"""Re-renders every figure a normal run would have produced, without
redoing work whose results already sit on the disk.
The following figures are enforced:
- ``complexity.png`` / ``score_v_comp.svg`` (from ``complexities.csv``
and ``udf.csv``; computed in :meth:`Simulation.__init__`),
- per realization: ``connectivity-rlzXXX.png``, ``tau_dist-rlzXXX.png``,
per-network ``uvx-<state>.png`` and ``score_score_{trn|tst}*.png``,
- once all realizations are done: ``angles*.png``,
``overlap_cdf*.png``, ``svals*.png`` and ``erank*.png``.
For each figure the required inputs are loaded from disk when they
already exist, and computed only when missing:
- topology / noise / timescales (``w*.npy``, ``A_adj.npz``, ``xi.npy``,
``<net_id>/taus.npy``) when absent,
- network states (per-realization ``<state>-rlzXXX.*``) when absent,
- readout scores (``trn_/tst_*.csv``) and minimal angles
(``angle_*.csv``) when absent.
.. note:: The topology, noise, and timescale files are *shared* across
realizations (overwritten by each ``synthesize_*`` call), so to
reproduce the per-realization connectivity/tau figures faithfully
the deterministic synthesis chain is replayed realization after
realization from the simulation seed (``self.rng`` is consumed only
by those ``synthesize_*`` calls). The figures are therefore
recomputed from the *replayed* (identical) values rather than from
stale on-disk leftovers of a later realization.
"""
# the input-level summary figures: ensure their inputs, then redraw
self.compute_udf()
self.plot_complexity_map()
self.plot_udf_vs_complexity()
if n_realizations != self.n_realizations:
self.n_realizations = n_realizations
state_name = self.generate_state_name(dyn=dyn, solver=solver, **gain_params)
self.setup_runner_log(n_realizations, dyn, solver,
gain_params, optim_params, readout_params)
# Replay the deterministic synthesis chain from the seed so that every
# realization's topology/noise/taus (and hence its figures) are exact.
# In a normal run self.rng is consumed *only* by the synthesize_* calls,
# so resetting it and replaying in order reproduces the original values.
self.rng = np.random.default_rng(self.seed)
print(f'Replot mode: replaying per-realization synthesis from seed '
f'{self.seed} to redraw the per-realization figures faithfully.')
for rlz in range(n_realizations):
state_name_rlz = state_name + self.rlz_suffix(rlz)
print(f'\n[replot] realization {rlz} -> {state_name_rlz}')
# --- synthesis artifacts (overwritten per realization) -----------
self.synthesize_topology()
self.synthesize_noise()
self.synthesize_timescales()
self.save_rng_state(state_name)
# per-realization figures whose inputs are now on disk
self.plot_connectivity(rlz)
self.plot_taus(rlz)
# --- network states, integrated only when missing ----------------
if not self.state_files_exist(state_name_rlz, dyn=dyn):
print(f'[replot] states of {state_name_rlz} missing; integrating.')
self.run_dynamic(
**gain_params, verbose=False,
dyn=dyn, dyn_params=dyn_params,
solver=solver, ncpu=ncpu,
rlz=rlz, purge_old=True,
)
self.rlog['integrated'] = self.get_n_integrated(
n_realizations, state_name, dyn)
self.update_runner(state_name)
# --- per-network uvx traces & minimal angles (from disk states) ---
self.replot_net_figs(state_name_rlz, readout_params)
# --- readout scores, fitted only when missing ---------------------
if not self.benchmark_files_exist(state_name_rlz,
readout_params, optim_params):
print(f'[replot] scores of {state_name_rlz} missing; benchmarking.')
self.benchmark(
state_name=state_name_rlz,
readout_params=readout_params,
optim_params=optim_params,
n_cpu=ncpu,
skip_uvx_plot=True, # uvx already re-rendered above
)
self.rlog['benchmarked'] = self.get_n_benchmarked(
n_realizations, state_name, optim_params, readout_params)
self.update_runner(state_name)
# per-realization score figure
self.plot_scores(rlz)
# once every realization is settled: the summary angle figures
self.plot_angle_vs_score()
self.plot_overlap_cdf()
self.plot_svals_erank()
[docs]
def replot_net_figs(self, state_name, readout_params, var_frac=0.01):
"""Re-renders the per-network ``uvx-<state>.png`` traces from the saved
states of ``state_name`` and computes the minimal angle csv (and the
accompanying singular-value spectrum and effective-rank csvs) of every
network that does not have them yet.
This mirrors the per-network block of :meth:`benchmark` but *only*
re-derives the visualization-related quantities: the states are loaded
(from disk) and plotted, and missing ``angle_/svals_/erank_{state_name}.csv``
files are filled, without re-training any readout.
"""
_, U = utils.load_inputs(self.res_path, load_noise=False)
Y = utils.load_target(self.res_path)
s_trn, s_tst = self.Strn, self.Stst
# prepare the (fixed) task subspaces once for the angle computation
task_bases = dict(
trn=self.make_task_basis(Y[:, s_trn], var_frac=var_frac),
tst=self.make_task_basis(Y[:, s_tst], var_frac=var_frac),
)
for net_id in self.net_ids:
state_path = self.get_state_path(net_id, state_name)
if not os.path.exists(state_path):
print(f'[replot] state of network {net_id} ({state_name}) '
f'not found; skipping its uvx/angle figures.')
continue
V, X = dataloader.load_states(
sim_path=self.res_path, net_id=net_id, state_name=state_name,
activation_params=readout_params['readout_activation'],
n_neuron=self.metadata['net']['N'],
n_timesteps=self.metadata['n_timesteps'],
)
self.plot_uvx(U, V, X, net_id, state_name)
self.compute_minimal_angles(
net_id=net_id,
state_name=state_name,
X=X,
slices=(s_trn, s_tst),
task_bases=(task_bases['trn'], task_bases['tst']),
var_frac=var_frac,
)
[docs]
def integrate_and_benchmark(
self, n_realizations, skip_integration, skip_benchmark,
dyn, dyn_params, gain_params, readout_params, optim_params,
solver, ncpu, skip_connectivity_plot=False, skip_score_plot=False,
skip_tau_plot=False, skip_uvx_plot=False, skip_angle_plot=False,
skip_summary_plot=False,
skip_svd=False,
replot=False
):
"""
Integrates and benchmarks multiple realizations of the simulation.
This method encapsulates the entire realization loop, handling random
variable regeneration, integration, and benchmarking for multiple
network realizations. It manages RNG state persistence to ensure
reproducibility and resumability.
Parameters
----------
n_realizations : int
Total number of realizations to complete
skip_integration : bool
If True, skip the integration step
skip_benchmark : bool
If True, skip the benchmarking step
dyn : str
Neural dynamic type (e.g., 'LI', 'LIF')
dyn_params : dict
Parameters for the neural dynamic
gain_params : dict
Gain parameters (Ju, J, Jn)
readout_params : dict
Parameters for readout layer
optim_params : dict
Parameters for optimization
solver : str
Numerical solver (e.g., 'euler')
ncpu : int
Number of CPU cores to use
skip_angle_plot : bool
If True, skip the summary angle-vs-score plot produced once all
realizations are benchmarked.
skip_summary_plot : bool
If True, skip the two summary figures produced once all
realizations are benchmarked -- the angle CDF
(:meth:`plot_overlap_cdf`) and the singular-value/effective-rank
figures (:meth:`plot_svals_erank`). Both aggregate *every* state on
disk, so in a parameter sweep (where the same folder accumulates
one subfolder set per case) they are re-rendered from all the
previously computed cases after each single run. Turn them off to
integrate the whole sweep and render the summaries once at the end.
replot : bool
If True, re-render every figure of the simulation instead of
running the usual synthesis/integration/benchmark pipeline. Each
figure's inputs are loaded from the disk when already computed and
computed (synthesis, integration, benchmark) only when missing.
The plot-skip flags above are ignored in this mode.
Returns
-------
int
Number of realizations completed after this call
Notes
-----
This method ensures that random variables (timescales, topology, noise)
are regenerated for each new realization by calling the synthesize_*
methods in sequence. It also manages RNG state persistence by calling
sync_rng() at the start and save_rng_state() after regeneration.
"""
if replot:
return self.replot_all(
n_realizations=n_realizations,
dyn=dyn,
dyn_params=dyn_params,
gain_params=gain_params,
readout_params=readout_params,
optim_params=optim_params,
solver=solver,
ncpu=ncpu,
)
self.setup_runner_log(n_realizations, dyn, solver,
gain_params, optim_params, readout_params)
# Generate state name for this realization to idenitfy if this run
# is now or not
state_name = self.generate_state_name(
dyn=dyn,
solver=solver,
**gain_params
)
# print('In ib', state_name)
# Update n_realizations if different from current
if n_realizations != self.n_realizations:
self.n_realizations = n_realizations
n_integrated = self.rlog['integrated']
n_benchmarked = self.rlog['benchmarked']
print(n_integrated, n_benchmarked)
# check if some benchmarking is left
if n_integrated > n_benchmarked:
assert n_benchmarked == (n_integrated - 1), "What has happened??"
rlz = n_integrated - 1 # likely this is broken
state_name_rlz = state_name + self.rlz_suffix(rlz)
self.benchmark(
state_name=state_name_rlz, # rlz info cuz we must load it
readout_params=readout_params,
optim_params=optim_params,
n_cpu=ncpu, skip_svd=skip_svd,
skip_uvx_plot=skip_uvx_plot,
)
self.rlog['benchmarked'] += 1
self.update_runner(state_name)
if not skip_score_plot:
self.plot_scores(rlz)
# continue with the integration and benchmarking
for rlz in range(n_integrated, n_realizations):
state_name_rlz = state_name + self.rlz_suffix(rlz)
# rlz_id = iteration + 1
# Check if we need a new realization
# if not self.needs_new_realization():
# break
# Sync RNG to continue from saved state
self.sync_rng(state_name=state_name, verbose=True)
# Regenerate all random variables for new realization
# print(f'before topoloy @ {self.rng.bit_generator.state["state"]["state"]}')
self.synthesize_topology()
if not skip_connectivity_plot:
self.plot_connectivity(rlz)
# print(f'before noise @ {self.rng.bit_generator.state["state"]["state"]}')
self.synthesize_noise()
# print(f'before tau @ {self.rng.bit_generator.state["state"]["state"]}')
self.synthesize_timescales()
if not skip_tau_plot:
self.plot_taus(rlz)
# Save RNG state after regeneration
self.save_rng_state(state_name)
# print(f'after save @ {self.rng.bit_generator.state["state"]["state"]}')
# Integration
if not skip_integration:
self.run_dynamic(
**gain_params, verbose=False,
dyn=dyn, dyn_params=dyn_params,
solver=solver, ncpu=ncpu,
rlz = rlz, purge_old=True,
)
self.rlog['integrated'] += 1
self.update_runner(state_name)
# integrated = self.mark_integrated(rlz, state_name, dyn=dyn)
# if not integrated:
# print(f'Integration for {state_name} is incomplete.')
else:
print('Skipping integration.')
# integrated = True
# integrated = self.mark_integrated(iteration, state_name, dyn=dyn)
# Benchmarking
if not skip_benchmark:# and integrated:
self.benchmark(
state_name=state_name_rlz, # rlz info cuz we must load it
readout_params=readout_params,
optim_params=optim_params,
n_cpu=ncpu,
skip_uvx_plot=skip_uvx_plot,
)
self.rlog['benchmarked'] += 1
self.update_runner(state_name)
if not skip_score_plot:
self.plot_scores(rlz)
# benchmarked = self.mark_benchmarked(rlz, state_name, readout_params, optim_params)
# if not benchmarked:
# print(f'Benchmarking for {state_name} is incomplete.')
else:
print('Skipping benchmarking.')
# else:
# print(f'Skipping benchmarking because {state_name} is not fully integrated.')
# Update metadata with new realization count
# self.update_num_realized(rlz_id)
# self.update_metadata()
# Check progress
# num_completed = self.progress_count(skip_benchmark=skip_benchmark)
# if num_completed <= initial_realized:
# print('No additional realization was completed in this run.')
# break
# once all realizations are benchmarked: one summary figure of the
# minimal angle of every task vs its score, across all networks.
if not skip_angle_plot:
self.plot_angle_vs_score()
# ------------------------------------------------------------------
# All networks of this realization are benchmarked and their per-task
# angles are on disk: aggregate everything into the angle CDF figure
# and the singular-value/effective-rank summary figures. Both pool
# *all* the states of the folder, so they are skipped in bulk sweeps
# (skip_summary_plot) and rendered once, at the end, instead.
# ------------------------------------------------------------------
rlz_id = None
if self.rlz_naming and 'rlz' in state_name:
rlz_id = int(state_name.rsplit('rlz', 1)[-1])
if not skip_summary_plot:
self.plot_overlap_cdf(rlz_id=rlz_id)
self.plot_svals_erank(rlz_id=rlz_id)
# return self.num_realized
# def visualize(self):
# plot_score_profile()