Simulation

Simulation, the object that owns a parameter tree and a result directory. It names the run, enumerates the (mean-timescale x heterogeneity) network grid, integrates each realization, benchmarks it with a linear readout, computes the complexity / angle / effective-rank measures, and renders the figures.

class hetero.sim.Simulation(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)

Bases: object

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.

generate_sim_name()

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.

generate_state_name(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.

generate_readout_name(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.

generate_optim_name(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.

get_state_path(net_id, state_name, dyn=None)

state_files_exist(state_name, dyn=None)

score_files_exist(score_name, mode='tst')

generate_score_template(state_name, readout_params, optim_params, ntrial='')

get_score_paths_for_template(net_id, template)

get_score_paths(net_id, state_name, readout_params, optim_params, ntrial)

scan_unbenchmarked(state_name, readout_params, optim_params, verbose=False)

benchmark_files_exist(state_name, readout_params, optim_params)

progress_count(skip_benchmark=False)

setup_seed(seed)

read_metadata()

Checks if there’s already an available metadata. If so, loads it, otherwise initializes the previous metadata with an empty dictionary.

update_metadata()

Updates the metadata file with the current simulation parameters, and saves it on the disk.

mark_integrated(rlz_id, state_name, dyn=None)

mark_benchmarked(rlz_id, state_name, readout_params, optim_params)

get_param_from_net_id(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.

generate_net_ids()

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.

sync_rng(state_name, verbose=True)

save_rng_state(state_name)

needs_new_realization(rlz)

rlz_suffix(rlz_id=0)

setup_integration_interval()

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.

synthesize_timescales()

Generates and saves the heterogeneous timescales for all networks based on the heterogeneity parameters.

setup_task()

Initializes the task object on the specified task parameters, and further creates all subtasks.

setup_stim_n_target()

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:` au_{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.

compute_complexity()

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.

compute_udf(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 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.

synthesize_noise()

Generates and saves the white noise time series for the network.

synthesize_topology()

Creates the network topology and weights.

plot_connectivity(rlz_id=0)

Plots the connectivity matrix and the recurrent weight distribution of the current realization into the result path.

plot_taus(rlz_id=0)

Plots the distribution of the neuronal time constants of all networks for the current realization into the result path.

plot_scores(rlz_id=0)

Plots the heterogeneous-vs-homogeneous score of every task for the given realization after it has been benchmarked.

plot_angle_vs_score(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 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.

plot_overlap_cdf(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 hetero.viz.plot_overlap_cdf()).

Meant to be called once all networks are benchmarked; when rlz_id is None all realizations found are plotted.

plot_svals_erank(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 hetero.viz.plot_svals_erank()).

Meant to be called once all networks are benchmarked; when rlz_id is None all realizations found are plotted.

plot_uvx(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.

plot_complexity_map()

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.

plot_udf_vs_complexity()

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.

make_task_basis(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 benchmark(), this method is called a single time (before the per-network loop) and its output is fed to 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 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 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).

compute_minimal_angles(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 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 _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 make_task_basis() and prepared exactly once in 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 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 numpy.arcsin() where cos^2 >= 0.5 (small angles) and 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 _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].

get_result_path()

Creates the result path.

prepare_dirs()

Creates one subdirectory for each network in the result path.

run_dynamic(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.

benchmark(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 compute_minimal_angles()).

Warning

The rate-based dynamic is parallelized via joblib. Depending on theoperating system, one might need to change the parallelizing backend.

get_n_integrated(n_realizations, state_name, dyn)

returns the number of fully integrated realizations

get_n_benchmarked(n_realizations, state_name, optim_params, readout_params)

returns the number of realizations that are fully benchmarked

setup_runner_log(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.

update_runner(running_state_name)

replot_all(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 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.

replot_net_figs(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 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.

integrate_and_benchmark(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 (plot_overlap_cdf()) and the singular-value/effective-rank figures (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.