vibespin.utils package

Submodules

utils.equilibration module

Equilibration algorithms for Monte Carlo simulations.

Provides adaptive and two-start convergence equilibration routines, together with the underlying relaxation-time estimation and quasi-steady-state detection helpers that support them. These algorithms were previously split between utils.system (equilibrate functions) and utils.analysis (relaxation diagnostic helpers), which required a lazy cross-import to avoid circular dependencies. Unifying them here eliminates that workaround.

utils.equilibration.adaptive_equilibrate(sim: _Sim, *, min_steps: int, probe_steps: int = 500, factor: float = 50.0, max_steps: int = 200000) → int[source]

Equilibrate a simulation adaptively, extending the burn-in until the probe window covers factor integrated autocorrelation times.

After the mandatory min_steps burn-in, the function repeatedly runs a probe_steps measurement via sim.run() and computes tau_int from the resulting magnetization series. If probe_steps >= factor * tau_int, the probe spans enough correlation times for the initial state to be forgotten and equilibration is declared complete. Otherwise the probe has advanced the system state and the loop continues. In the ordered phase the magnetization series has zero variance; this is treated as full equilibration and the function returns immediately.

Parameters:

sim (Any simulation object implementing equilibrate and run.)
min_steps (Mandatory burn-in passed to sim.equilibrate before probing.)
probe_steps (MC steps per probe run (default 500).)
factor (Required ratio probe_steps / tau_int (default 50.0).)
max_steps (Hard cap on total steps to prevent unbounded runtime near) – criticality (default 200 000).

Return type:

Total number of MC steps run (burn-in + probes).

Raises:

ValueError – If the adaptive-equilibration parameters are invalid.:

utils.equilibration.convergence_equilibrate(sim_random: _Sim, sim_ordered: _Sim, *, chunk_size: int = 500, max_steps: int = 200000, qs_sigma_threshold: float = 0.05, qs_min_steps: int = 1500, **kwargs: Any) → int[source]

Equilibrate two simulations (random- and ordered-start) until they converge.

Uses estimate_relaxation_time_two_start to detect convergence via the mutual cross-band criterion: each smoothed trajectory must enter and sustain a band defined by the other trajectory’s tail statistics. A sigma floor prevents false positives when one trace is nearly flat. This is more robust than one-start adaptive methods for complex energy landscapes.

Parameters:

sim_random (Simulation instance started from a random state.)
sim_ordered (Simulation instance started from an ordered state.)
chunk_size (Number of steps to run between convergence checks.)
max_steps (Hard cap on total steps.)
qs_sigma_threshold (Tail-std threshold used to detect a quasi-steady,) – non-converged stuck state and exit early.
qs_min_steps (Minimum accumulated steps before stuck detection is allowed) – to fire. Ensures tail statistics are computed from enough data to be reliable, preventing false positives when traces have not yet had time to settle.
**kwargs (Passed to estimate_relaxation_time_two_start (k, smooth_window,) – dwell_window, min_fraction_inside, sigma_floor, etc.).

Return type:

Total number of MC steps run per simulation.

utils.equilibration.convergence_equilibrate_with_status(sim_random: _Sim, sim_ordered: _Sim, *, chunk_size: int = 500, max_steps: int = 200000, qs_sigma_threshold: float = 0.05, qs_min_steps: int = 1500, qs_allow_stuck: bool = False, **kwargs: Any) → tuple[int, bool][source]

Equilibrate two simulations and report whether convergence was reached.

Uses estimate_relaxation_time_two_start to detect convergence via the mutual cross-band criterion: each smoothed trajectory must enter and sustain a band defined by the other trajectory’s tail statistics. A sigma floor prevents false positives when one trace is nearly flat. This variant returns both the total number of MC steps executed and a boolean convergence flag.

Parameters:

sim_random (Simulation instance started from a random state.)
sim_ordered (Simulation instance started from an ordered state.)
chunk_size (Number of steps to run between convergence checks.)
max_steps (Hard cap on total steps.)
qs_sigma_threshold (Tail-std threshold used to detect a quasi-steady,) – non-converged stuck state and exit early.
qs_min_steps (Minimum accumulated steps before stuck detection is allowed) – to fire. Ensures tail statistics are computed from enough data to be reliable, preventing false positives when traces have not yet had time to settle.
qs_allow_stuck (If True, detecting a quasi-steady stuck state (where the) – ordered trace is stable but the random trace is stranded) is treated as a successful equilibration. Useful for Ising sweeps where random-start domain-wall trapping is physically expected.
**kwargs (Passed to estimate_relaxation_time_two_start (k, smooth_window,) – dwell_window, min_fraction_inside, sigma_floor, etc.).

Return type:

Tuple (total_steps, converged).

utils.equilibration.estimate_relaxation_time_two_start(trace_random: ndarray, trace_ordered: ndarray, *, k: float = 1.0, smooth_window: int = 60, dwell_window: int = 60, min_fraction_inside: float = 0.85, sigma_floor: float = 0.02, skip_validation: bool = False) → int[source]

Estimate thermalization time from convergence of random- and ordered-start traces.

Compares a trajectory started from a random spin configuration against one started from a fully ordered configuration. Returns the first step at which both smoothed traces pass a mutual cross-band test sustained over a dwell window.

The criterion is asymmetric in the correct physical sense: the smoothed random-start trace must lie within k standard deviations of the ordered-start tail mean, and the ordered-start trace must simultaneously lie within k standard deviations of the random-start tail mean. Tail statistics are computed from the second half of each smoothed trace independently. A sigma_floor prevents band collapse when one trace is nearly variance-free (e.g. in the deep ordered phase near T=0).

Parameters:

trace_random (1-D magnetization trace from a random initial state.)
trace_ordered (1-D magnetization trace from an ordered initial state.)
k (Half-width of each convergence band in units of the respective tail) – standard deviation. Larger values are more permissive.
smooth_window (Window length (steps) for the moving-average smoother.)
dwell_window (Window length (steps) for the sustained-convergence test.)
min_fraction_inside (Fraction of steps within dwell_window that must) – satisfy the mutual cross-band condition to declare convergence.
sigma_floor (Minimum allowed standard deviation for each tail, preventing) – band collapse when a trajectory is nearly flat.
skip_validation (If True, skip np.asarray and finite-prefix checks.)

Returns:

Estimated relaxation time as a 0-indexed position in the original input
traces (i.e. a sweep count minus one). The first index at which both
smoothed traces have been mutually inside each other’s band for
dwell_window steps is mapped back to the corresponding raw-trace
position as hits[0] + smooth_window - 1. Returns n (full trace
length) if no convergence is detected, or 0 for very short traces.

utils.evolution_helpers module

Shared infrastructure for ordering evolution visualisation scripts.

Provides the snapshot loop and figure generation shared by the Ising, XY, and Clock ordering-evolution scripts.

utils.evolution_helpers.run_ordering_evolution(*, model_cls: type, model_kwargs: dict[str, Any], capture_vorticity: bool, title: str, size: int, temp: float, step_targets: list[int], output_dir: str, logger: Logger | None = None) → None[source]

Run an ordering evolution simulation and save the multi-panel figure.

Quenches the model to temp, advances to each step in step_targets, captures spin configurations and correlation functions, optionally captures vorticity maps, and writes the figure to output_dir.

Parameters:

model_cls (type) – Simulation class to instantiate (e.g. XYSimulation).
model_kwargs (dict[str, Any]) – Extra keyword arguments forwarded to the constructor beyond size, temp, and update='random'.
capture_vorticity (bool) – Whether to record vorticity maps at each snapshot. Pass True for vector-spin models (XY, Clock) and False for scalar models (Ising). Controls is_vector in the output figure.
title (str) – Figure-level suptitle.
size (int) – Linear lattice size L.
temp (float) – Quench temperature T.
step_targets (list[int]) – MC step counts at which to capture snapshots. Need not be sorted; the function sorts them internally.
output_dir (str) – Directory in which to write the output figure.
logger (logging.Logger, optional) – Logger to use; creates a module-level logger if not provided.

utils.exceptions module

Project-specific exception types for analysis and runtime failures.

exception utils.exceptions.NumericalAnalysisError[source]

Bases: VibeSpinError, RuntimeError

Base class for mathematically undefined or failed analysis results.

exception utils.exceptions.VibeSpinError[source]

Bases: Exception

Base class for project-specific exceptions.

exception utils.exceptions.ZeroVarianceAutocorrelationError[source]

Bases: NumericalAnalysisError

Raised when autocorrelation is undefined because the series variance is zero.

utils.kinetics_helpers module

Shared infrastructure for ordering kinetics analysis scripts.

Provides the simulation loop, power-law fitting, and result plotting shared by the Ising, XY, and Clock ordering-kinetics scripts. The compute_mean_intercept_length function is also placed here as a reusable domain-size estimator for scalar-spin models.

utils.kinetics_helpers.compute_mean_intercept_length(sim: Any) → float[source]

Estimate domain size via the stereological mean intercept length (MIL).

Counts domain-wall crossings along every row and column of the lattice, including the periodic wrap-around crossing, then divides the total line length by the crossing count to obtain the mean intercept length.

Parameters:: sim (Any) – A simulation instance with a spins attribute (an integer-sign array) and a size property giving the linear lattice dimension.
Returns:: Estimated domain size in lattice units. Returns float(size) when there are no domain walls (fully ordered state).
Return type:: float

utils.kinetics_helpers.run_ordering_kinetics(*, model_cls: type, model_kwargs: dict[str, Any], third_metric_fn: Callable[[Any], float], third_metric_label: str, title: str, left_title: str, right_title: str, size: int, temp: float, max_steps: int, samples: int, fit_min: int, output_dir: str, y_label: str = 'Characteristic Length Scale $L(t)$ (lattice units)', logger: Logger | None = None, npz_path: str | None = None, n_seeds: int = 1, base_seed: int = 42) → None[source]

Run an ordering kinetics simulation and save the two-panel kinetics figure.

Quenches the model to temp, steps forward to each of samples logarithmically spaced targets up to max_steps, measures growth scales and a model-specific third metric, fits power laws, and writes the figure to output_dir.

When n_seeds > 1, the simulation is repeated for an ensemble of seeds. The median trajectory is used for fitting and plotting, and per-seed arrays together with IQR-based error bars are saved to the NPZ file.

Parameters:

model_cls (type) – Simulation class to instantiate (e.g. IsingSimulation).
model_kwargs (dict[str, Any]) – Extra keyword arguments forwarded to the constructor beyond size, temp, and update='random'.
third_metric_fn (Callable[[Any], float]) – Called with the live simulation instance at each measurement step; returns a scalar third metric (e.g. vortex density or MIL).
third_metric_label (str) – Y-axis label for the third-metric panel.
title (str) – Figure-level suptitle.
left_title (str) – Title for the left (growth-scale) subplot.
right_title (str) – Title for the right (third-metric) subplot.
size (int) – Linear lattice size L.
temp (float) – Quench temperature T.
max_steps (int) – Total MC steps to run.
samples (int) – Number of logarithmically spaced measurement points.
fit_min (int) – Minimum MC step used for power-law fitting.
output_dir (str) – Directory in which to write the output figure.
y_label (str, optional) – Y-axis label for the growth-scale panel.
logger (logging.Logger, optional) – Logger to use; creates a module-level logger if not provided.
npz_path (str, optional) – If given, save the kinetics data arrays to this .npz file for notebook consumption.
n_seeds (int, optional) – Number of independent seeds to run. Defaults to 1 (legacy single-trajectory behavior).
base_seed (int, optional) – Starting seed; seeds are base_seed, base_seed+1, ....

utils.observables module

Physics-related utility functions for calculating thermodynamic observables and correlations.

utils.observables.calculate_entropy(*, temperatures: ndarray, specific_heat: ndarray, s_ref: float = 0.0) → ndarray[source]

Compute entropy by integrating specific heat over a temperature sweep.

Integrates C(T)/T from high temperature downward so that S(T_max) = s_ref and S(T) = s_ref - integral from T to T_max of C/T dT’. This avoids the C/T divergence near T = 0 and anchors the result at the high-temperature limit where the entropy is analytically known.

Parameters:

temperatures (1-D array of temperature values. Need not be sorted.)
specific_heat (Specific heat C_v at each temperature (same length).)
s_ref (Reference entropy at the highest temperature. Pass 0.0) – for relative entropy or np.log(q) (per site, in units of k_B) for absolute entropy of a q-state model.

Returns:

1-D array of entropy values, one per temperature, in the original
unsorted order of temperatures.

Raises:

ValueError – If arrays differ in length, contain fewer than 2 points,: or include non-positive temperatures.

utils.observables.calculate_thermodynamics(*, mags: ndarray, engs: ndarray, T: float, L: int) → tuple[float, float, float, float][source]

Calculate average magnetization, energy, susceptibility, and specific heat.

Parameters:

mags (Array of magnetization measurements.)
engs (Array of energy measurements.)
T (Temperature.)
L (Linear lattice size.)

Return type:

A tuple of (avg_mag, avg_eng, susceptibility, specific_heat).

Raises:

ValueError – If T is not positive or L is not a positive integer.:

utils.observables.compute_kinetics_metrics(*, sim: _Sim) → dict[str, float][source]

Calculate common kinetics metrics (R_sk, xi) for a simulation state.

Parameters:: sim (MonteCarloSimulation instance.)
Return type:: Dictionary containing ‘R_sk’ and ‘xi’.

utils.observables.get_averaged_correlation(*, sim: _Sim, total_steps: int, sample_interval: int) → tuple[ndarray, ndarray][source]

Run simulation and average the correlation function over multiple configurations.

Parameters:

sim (An instance of MonteCarloSimulation.)
total_steps (Total number of MC steps to run.)
sample_interval (Interval between correlation samples. Must be ≥ 1.)

Returns:

r (Radial distances.)
G_r_avg (Averaged correlation values.)

Raises:

ValueError – If sample_interval is less than 1 or total_steps is negative.:

utils.observables.pair_correlation_x(*, spins: ndarray) → tuple[ndarray, ndarray][source]

Compute the real-space spin-spin pair correlation G(r) along x.

Uses the Wiener-Khinchin theorem: the autocorrelation of each row is the inverse FFT of the row’s power spectrum. Results are averaged over all rows (y positions) and normalised so G(0) = 1.

Parameters:

spins ((N, N) or (N, N, 2) spin array.)

Returns:

r_vals (Lag distances r = 0 … N//2 in lattice units.)
G (Normalised pair correlation G(r) / G(0).)

utils.observables.radial_average_sk(*, spins: ndarray) → tuple[ndarray, ndarray][source]

Compute the circularly averaged structure factor $S(|k|)$.

Bins $S(k)$ by integer pixel radius from the DC centre of the shifted FFT, then averages within each annular bin.

Parameters:

spins ((N, N) or (N, N, 2) spin array.)

Returns:

k_vals (Wavevector magnitudes in units of 2π/N (reciprocal lattice).)
S_radial (Mean S(k) value for each annular bin.)

utils.plotting module

Domain plotting and file-output helpers for VibeSpin simulations.

Groups the filesystem output helper (save_plot(), ensure_results_dir()) together with the three domain-specific visualisation functions that were previously embedded in utils.system. Keeping visualisation logic here makes system_helpers a pure infrastructure module (logging, parallelism) with no matplotlib dependency.

utils.plotting.ensure_results_dir(*, directory: str = 'results') → str[source]

Ensure the results directory exists.

Parameters:: directory (Name of the directory to create.)
Return type:: The path to the directory.

utils.plotting.plot_ordering_evolution(*, targets: Sequence[int], snapshots: Sequence[ndarray], gr_data: Sequence[tuple[ndarray, ndarray]], vorticity_data: Sequence[ndarray] | None, title: str, filename: str, directory: str, is_vector: bool = False) → None[source]

Generate and save a multi-row figure showing order evolution over time.

Row 0: Spin configurations (binary for Ising, HSV for XY/Clock). Row 1: Circularly-averaged structure factor $S(|k|)$ or vorticity maps. Row 2: Radially averaged correlation functions G(r) with xi estimates.

Parameters:

targets (List of MC steps for each snapshot.)
snapshots (List of spin arrays.)
gr_data (List of (r, G) tuples.)
vorticity_data (Optional list of vorticity arrays.)
title (Figure-level suptitle.)
filename (Output filename.)
directory (Output directory.)
is_vector (Whether the spins are 2D vectors (True) or scalars (False).)

utils.plotting.plot_ordering_evolution_snapshots(*, targets: Sequence[int], snapshots: Sequence[ndarray], gr_data: Sequence[tuple[ndarray, ndarray]], sk_data: Sequence[tuple[ndarray, ndarray]], size: int, temp: float) → tuple[Figure, ndarray, list[float]][source]

Generate a 3-row diagnostic grid of ordering evolution snapshots.

Row 1 shows binary spin matrices. Row 2 shows the structure factor $S(|k|)$ on log-log axes. Row 3 shows the spatial correlation function G(r) with the 1/e crossing (correlation length xi) annotated.

Parameters:

targets (Target Monte Carlo sweep counts for each column.)
snapshots (Spin matrices, one per target.)
gr_data (Correlation data tuples (r, G(r)), one per target.)
sk_data (Structure factor tuples (k, S(|k|)), one per target.)
size (Lattice size L used in the suptitle.)
temp (Quench temperature used in the suptitle.)

Returns:

fig (The matplotlib Figure object.)
axes (The matplotlib Axes grid.)
xi_values (List of extracted correlation lengths (np.nan if undefined).)

utils.plotting.plot_ordering_kinetics(*, t: ndarray, R_sk: ndarray, R_xi: ndarray, third_metric: ndarray | None, third_metric_label: str | None, exponents: dict[str, float | None], prefactors: dict[str, float | None], fit_mask: ndarray, title: str, filename: str, directory: str, y_label: str = 'Characteristic Length Scale $L(t)$ (lattice units)', left_title: str = 'Domain Coarsening', right_title: str = 'Defect/Boundary Evolution') → None[source]

Generate and save a 2-panel figure showing ordering kinetics.

Left Panel: Growth of length scales R_sk and xi. Right Panel: Growth/Decay of a third metric (Vortex density or MIL).

Parameters:

t (Time array (Monte Carlo sweeps).)
R_sk (Domain size from structure factor.)
R_xi (Correlation length from G(r).)
third_metric (Optional third metric array.)
third_metric_label (Label for the third metric.)
exponents (Dict of fitted exponents.)
prefactors (Dict of fitted prefactors.)
fit_mask (Mask used for fitting.)
title (Figure-level suptitle.)
filename (Output filename.)
directory (Output directory.)
y_label (Y-axis label for the left panel.)
left_title (Title for the left subplot.)
right_title (Title for the right subplot.)

Generate and save a standardized temperature sweep plot.

Displays magnetization, energy, susceptibility, and specific heat as functions of temperature in a dedicated thermodynamics figure. When entropy or tau_int are provided a companion diagnostics figure is saved alongside the main figure.

Parameters:

temperatures (Array of temperature values (x-axis).)
avg_m (Average absolute magnetization per temperature point.)
avg_e (Average energy per temperature point.)
susc (Magnetic susceptibility per temperature point.)
spec_h (Specific heat per temperature point.)
title (Figure-level suptitle string (e.g. '2D Ising Model: Temperature Sweep (L=50)').)
filename (Output filename passed to save_plot() (e.g. ‘temperature_sweep.png’).)
directory (Output directory passed to save_plot() (e.g. ‘results/ising’).)
entropy (Optional entropy per temperature point.)
entropy_err (Optional symmetric uncertainty for entropy.)
entropy_ci_low (Optional lower confidence band for entropy.)
entropy_ci_high (Optional upper confidence band for entropy.)
tau_int (Optional integrated autocorrelation time per temperature point.) – Reveals critical slowing down as a peak near T_c.
tau_int_ci_low (Optional lower confidence band for tau_int.)
tau_int_ci_high (Optional upper confidence band for tau_int.)
tau_unstable_flag (Optional per-temperature mask flagging unstable) – tau_int intervals.
low_effective_sample_flag (Optional per-temperature mask identifying) – points with low effective sample size.
diagnostics_note (Optional text shown on diagnostics figure summarizing) – uncertainty quality metrics.
run_metadata_note (Optional concise metadata line for the diagnostics) – header (for example lattice size, seeds, confidence level, method).
quality_summary (Optional dictionary containing quality counts. Expected) – keys are total_points, well_conditioned_count, low_effective_count, unstable_interval_count, and undefined_count.
transition_temperatures (Optional mapping of transition-marker labels) – to temperatures rendered as vertical dashed guide lines.
transition_window (Optional (low, high) temperature window highlighted) – with a light background band.
entropy_reference (Optional (label, value) rendered as a horizontal) – reference on the entropy diagnostics panel.
annotate_peaks (If True, annotate peak temperatures for susceptibility,) – specific heat, and tau_int when present.
min_visible_rel_error (Minimum relative error bar size used for visibility) – when finite uncertainties are extremely small.
mark_invalid_uncertainty (If True, points with non-finite uncertainty are) – marked with x markers and a legend label.

utils.plotting.save_plot(*, filename: str, directory: str = 'results', tight_layout: bool = True) → None[source]

Save the current matplotlib plot to the results directory.

Parameters:

filename (Name of the output file (e.g., 'plot.png').)
directory (Output directory name.)
tight_layout (Whether to apply plt.tight_layout() before saving.)

utils.statistics module

Physics-related utility functions for calculating thermodynamic observables and correlations.

utils.statistics.blocking_error(*, time_series: ndarray, min_block_size: int = 2, max_block_size: int | None = None) → dict[str, float][source]: Estimate standard error using standard blocking/binning analysis.

utils.statistics.calculate_autocorr(*, time_series: ndarray, max_lag: int | None = None, window_constant: float = 6.0) → tuple[ndarray, float][source]

Compute the normalized autocorrelation function and integrated autocorrelation time.

Uses an FFT-based estimator for C(t): zero-pads to 2N, computes the power spectrum, IFFTs, and normalizes by variance and the number of pairs at each lag. Applies the Madras-Sokal automatic windowing rule to stop integration at the lag W where W = window_constant * tau_int, preventing the estimator from accumulating noise at large lags.

Parameters:

time_series (1-D array of sequential measurements (e.g. magnetization).)
max_lag (Maximum lag to compute. Defaults to len(time_series) // 2.)
window_constant (Multiplier for the Madras-Sokal window (default 6.0).)

Returns:

A tuple (C_t, tau_int) where C_t is the normalized autocorrelation
array truncated at the window endpoint, and tau_int is the integrated
autocorrelation time (scalar).

Raises:

ValueError – If time_series has fewer than 3 elements.:
ZeroVarianceAutocorrelationError – If time_series has zero variance.:

utils.statistics.estimate_effective_sample_size(*, time_series: ndarray, tau_int: float | None = None) → float[source]: Estimate effective sample size from integrated autocorrelation time.

utils.statistics.power_fit(*, t_arr: ndarray, y_arr: ndarray, mask: ndarray) → tuple[float, float] | tuple[None, None][source]

Fit a power law y = prefactor * t^exponent via log-log linear regression.

Parameters:

t_arr (Independent variable array (e.g., time or length scale).)
y_arr (Dependent variable array; only positive values are used in the fit.)
mask (Boolean mask selecting the subset of points to include.)

Returns:

exponent (Fitted power-law exponent, or None if fewer than 3 valid points.)
prefactor (Fitted prefactor, or None if fewer than 3 valid points.)

utils.statistics.summarize_asymmetric_replicate_uncertainty(*, samples: ndarray, confidence: float = 0.68) → dict[str, float][source]: Summarize 1-D replicate samples with asymmetric percentile intervals.

utils.statistics.summarize_derived_observable(*, magnetization_series: ndarray | None = None, energy_series: ndarray | None = None, temperature: float, L: int, observable: str, method: str = 'blocking', confidence: float = 0.68, bootstrap_resamples: int = 0) → dict[str, float][source]: Summarize derived observables (chi, Cv) with uncertainty estimates.

utils.statistics.summarize_entropy_observable(*, temperatures: ndarray, specific_heat_samples: ndarray, confidence: float = 0.68, s_ref: float = 0.0, specific_heat_err: ndarray | None = None, method: str = 'bootstrap', bootstrap_resamples: int = 400, rng_seed: int = 0) → dict[str, ndarray][source]

Summarize entropy from replicate specific-heat curves.

Each replicate column is converted to an entropy curve via calculate_entropy(), then uncertainty is computed either via replicate spread (blocking-like) or bootstrap over replicate curves. For a single finite replicate, uncertainty is reported as NaN unless specific_heat_err is provided, in which case entropy uncertainty is propagated from pointwise specific-heat errors.

utils.statistics.summarize_primary_observable(*, time_series: ndarray, confidence: float = 0.68) → dict[str, float][source]: Summarize a primary observable with blocking-based uncertainty.

utils.statistics.summarize_replicate_samples(*, samples: ndarray, confidence: float = 0.68) → dict[str, ndarray | float][source]: Summarize replicate samples using NaN-safe median and percentile bands.

utils.statistics.summarize_seed_ensemble(*, values: ndarray, within_seed_errors: ndarray, confidence: float = 0.68) → dict[str, float][source]

Aggregate per-seed estimates into one mean with hierarchical uncertainty.

The total variance combines between-seed and within-seed components: Var(mean) = Var_between / n_seeds + mean(Var_within) / n_seeds.

utils.sweep_helpers module

Shared infrastructure for per-temperature Monte Carlo sweep workers.

This module provides the two-layer decomposition used by all three model temperature sweeps (Ising, XY, Clock).

simulate_at_temperature owns equilibration and raw measurement. compute_thermo_observables owns statistical summarization. simulate_thermo_point chains both layers and is the worker function passed to parallel_sweep.

Post-processing helpers build_uncertainty_bundle and build_quality_flags aggregate per-seed results over a full temperature axis and are also shared here, eliminating identical copies that previously lived in each sweep script.

class utils.sweep_helpers.RawThermoData(temperature_index: int, seed_index: int, equilibrated_flag: float, equilibration_steps: float, temperature: float, size: int, mags_arr: np.ndarray | None, engs_arr: np.ndarray | None)[source]

Bases: NamedTuple

Raw simulation output from one (temperature, seed) equilibration + measurement.

Parameters:

temperature_index (int) – Index of this temperature in the sweep’s temperature array.
seed_index (int) – Replica index within the multi-seed ensemble for this temperature.
equilibrated_flag (float) – 1.0 if equilibration converged, 0.0 otherwise.
equilibration_steps (float) – Total equilibration steps taken (for diagnostics).
temperature (float) – Temperature at which the simulation was run.
size (int) – Linear lattice dimension L.
mags_arr (np.ndarray or None) – Magnetisation time series of length meas_steps. None when equilibration did not converge.
engs_arr (np.ndarray or None) – Energy per site time series of length meas_steps. None when equilibration did not converge.

engs_arr: ndarray | None: Alias for field number 7

equilibrated_flag: float: Alias for field number 2

equilibration_steps: float: Alias for field number 3

mags_arr: ndarray | None: Alias for field number 6

seed_index: int: Alias for field number 1

size: int: Alias for field number 5

temperature: float: Alias for field number 4

temperature_index: int: Alias for field number 0

class utils.sweep_helpers.ThermoPoint(temperature: float, size: int, meas_steps: int, eq_probe_steps: int, eq_max_steps: int, eq_qs_sigma_threshold: float, eq_qs_min_steps: int, qs_allow_stuck: bool, prefer_ordered_start: bool, temperature_index: int, seed_index: int, seed: int, model_cls: type, model_kwargs: dict[str, Any], confidence: float, derived_method: str, bootstrap_resamples: int)[source]

Bases: NamedTuple

Typed worker payload for one (temperature, seed) point in a sweep.

Parameters:

temperature (float) – Simulation temperature T.
size (int) – Linear lattice dimension L.
meas_steps (int) – Number of Monte Carlo sweeps to record after equilibration.
eq_probe_steps (int) – Chunk size passed to convergence_equilibrate_with_status.
eq_max_steps (int) – Hard cap on total equilibration steps.
eq_qs_sigma_threshold (float) – Tail-std threshold for quasi-steady stuck-state detection.
eq_qs_min_steps (int) – Minimum steps before stuck detection is allowed to fire.
qs_allow_stuck (bool) – If True, a quasi-steady stuck state counts as convergence (Ising low-T).
prefer_ordered_start (bool) – If True, switch to the ordered-start simulation when its magnetisation substantially exceeds the random-start value. Appropriate for Ising below T_c.
temperature_index (int) – Index of this temperature in the sweep’s temperature array.
seed_index (int) – Replica index within the multi-seed ensemble for this temperature.
seed (int) – RNG seed for both the random- and ordered-start simulations.
model_cls (type) – Simulation class to instantiate (e.g. IsingSimulation). Must be importable by its qualified name so it survives multiprocessing pickle.
model_kwargs (dict[str, Any]) – Extra keyword arguments forwarded to the model constructor beyond size, temp, init_state, and seed. Empty dict for Ising and XY; carries {'q': q, 'A': aniso} for Clock.
confidence (float) – Two-sided confidence level for uncertainty intervals (0 < c < 1).
derived_method (str) – Uncertainty method for susceptibility and specific heat ('blocking' or 'bootstrap').
bootstrap_resamples (int) – Bootstrap resamples used when derived_method='bootstrap'.

bootstrap_resamples: int: Alias for field number 16

confidence: float: Alias for field number 14

derived_method: str: Alias for field number 15

eq_max_steps: int: Alias for field number 4

eq_probe_steps: int: Alias for field number 3

eq_qs_min_steps: int: Alias for field number 6

eq_qs_sigma_threshold: float: Alias for field number 5

meas_steps: int: Alias for field number 2

model_cls: type: Alias for field number 12

model_kwargs: dict[str, Any]: Alias for field number 13

prefer_ordered_start: bool: Alias for field number 8

qs_allow_stuck: bool: Alias for field number 7

seed: int: Alias for field number 11

seed_index: int: Alias for field number 10

size: int: Alias for field number 1

temperature: float: Alias for field number 0

temperature_index: int: Alias for field number 9

utils.sweep_helpers.build_quality_flags(*, tau_int: ndarray, ci_low: ndarray, ci_high: ndarray, n_eff: ndarray, min_effective_samples: float, max_tau_relative_width: float) → dict[str, ndarray][source]

Build per-temperature diagnostic flags for uncertainty quality.

Parameters:

tau_int (np.ndarray) – Shape (n_temps,). Integrated autocorrelation time estimates.
ci_low (np.ndarray) – Shape (n_temps,). Lower CI bound on tau_int.
ci_high (np.ndarray) – Shape (n_temps,). Upper CI bound on tau_int.
n_eff (np.ndarray) – Shape (n_temps,). Effective sample size estimates.
min_effective_samples (float) – Threshold below which a point is flagged as low-N_eff.
max_tau_relative_width (float) – CI-width-to-tau ratio above which the tau estimate is considered unstable.

Returns:

Keys: 'undefined_autocorr_flag', 'low_effective_sample_flag', 'tau_interval_unstable_flag'. All uint8 arrays of length n_temps.

Return type:

dict[str, np.ndarray]

utils.sweep_helpers.build_uncertainty_bundle(*, values_by_seed: ndarray, errors_by_seed: ndarray, tau_by_seed: ndarray, n_eff_by_seed: ndarray, confidence: float) → dict[str, ndarray][source]

Aggregate per-seed results for one observable across a temperature axis.

For multi-seed sweeps, applies hierarchical seed-ensemble aggregation via summarize_seed_ensemble. Falls back to direct single-seed blocking results when only one seed is available.

Parameters:

values_by_seed (np.ndarray) – Shape (n_temps, n_seeds). Per-seed point estimates.
errors_by_seed (np.ndarray) – Shape (n_temps, n_seeds). Per-seed within-seed uncertainties.
tau_by_seed (np.ndarray) – Shape (n_temps, n_seeds). Per-seed integrated autocorrelation times.
n_eff_by_seed (np.ndarray) – Shape (n_temps, n_seeds). Per-seed effective sample sizes.
confidence (float) – Two-sided confidence level for the returned CI bounds.

Returns:

Keys: 'value', 'err', 'ci_low', 'ci_high', 'tau_int', 'n_eff', 'samples'. All arrays have length n_temps except 'samples' which has shape (n_temps, n_seeds).

Return type:

dict[str, np.ndarray]

utils.sweep_helpers.compute_thermo_observables(raw: RawThermoData, point: ThermoPoint) → dict[str, float][source]

Compute blocked/summarized thermodynamic observables from raw time series.

Takes the raw output of simulate_at_temperature and applies autocorrelation-aware blocking (primary observables) and variance-based summarization (derived observables) to produce the standard per-point result dictionary.

Parameters:

raw (RawThermoData) – Raw time series from simulate_at_temperature.
point (ThermoPoint) – Originating worker payload; provides temperature, size, and statistical configuration (confidence, derived_method, bootstrap_resamples).

Returns:

Keys: temperature_index, seed_index, equilibrated_flag, equilibration_steps, and — when equilibrated — avg_m_*, avg_e_*, susc_*, spec_h_* with suffixes _value, _err, _tau_int, _n_eff.

Return type:

dict[str, float]

utils.sweep_helpers.simulate_at_temperature(point: ThermoPoint) → RawThermoData[source]

Run equilibration and measurement for one (temperature, seed) point.

Creates a random-start and an ordered-start simulation with the same seed, runs two-start convergence equilibration, then records meas_steps sweeps on the active simulation. Returns raw magnetisation and energy time series for downstream statistical processing.

Parameters:: point (ThermoPoint) – Fully configured worker payload.
Returns:: Raw time series; mags_arr and engs_arr are None when equilibration did not converge.
Return type:: RawThermoData

utils.sweep_helpers.simulate_thermo_point(point: ThermoPoint) → dict[str, float][source]

Run a full (equilibration + measurement + summarization) point.

Thin wrapper over simulate_at_temperature followed by compute_thermo_observables. This is the function passed as worker_func to parallel_sweep in temperature-sweep scripts.

Parameters:: point (ThermoPoint) – Fully configured worker payload.
Returns:: Standard per-point result dictionary; see compute_thermo_observables for the full key set.
Return type:: dict[str, float]

utils.system module

Technical infrastructure for logging, parallel execution, and CLI argument parsing.

utils.system.parallel_sweep(*, worker_func: Callable, params: Iterable, num_processes: int | None = None) → list[source]

Run a parallel sweep over a set of parameters using a worker function. Uses multiprocessing.Pool and tqdm for progress tracking.

Parameters:

worker_func (Function to execute in parallel.)
params (Iterable of parameters to pass to the worker function.)
num_processes (Number of processes to use. Defaults to CPU count.)

Return type:

List of results from the worker function.

utils.system.parse_args_compat(parser: ArgumentParser) → Namespace[source]: Parse CLI args with compatibility for parser wrappers used by some runners.

utils.system.setup_logging(*, level: int = 20, log_file: str | None = None) → Logger[source]

Configure project-wide logging.

Parameters:

level (Logging level (e.g., logging.INFO).)
log_file (Optional path to a file to save logs to.)

Return type:

The configured logger instance.

Module contents

Public API for VibeSpin utility and analysis functions.

exception utils.NumericalAnalysisError[source]

Bases: VibeSpinError, RuntimeError

Base class for mathematically undefined or failed analysis results.

exception utils.VibeSpinError[source]

Bases: Exception

Base class for project-specific exceptions.

exception utils.ZeroVarianceAutocorrelationError[source]

Bases: NumericalAnalysisError

Raised when autocorrelation is undefined because the series variance is zero.

utils.blocking_error(*, time_series: ndarray, min_block_size: int = 2, max_block_size: int | None = None) → dict[str, float][source]: Estimate standard error using standard blocking/binning analysis.

utils.parse_args_compat(parser: ArgumentParser) → Namespace[source]: Parse CLI args with compatibility for parser wrappers used by some runners.

utils.summarize_derived_observable(*, magnetization_series: ndarray | None = None, energy_series: ndarray | None = None, temperature: float, L: int, observable: str, method: str = 'blocking', confidence: float = 0.68, bootstrap_resamples: int = 0) → dict[str, float][source]: Summarize derived observables (chi, Cv) with uncertainty estimates.

utils.summarize_primary_observable(*, time_series: ndarray, confidence: float = 0.68) → dict[str, float][source]: Summarize a primary observable with blocking-based uncertainty.

utils.summarize_seed_ensemble(*, values: ndarray, within_seed_errors: ndarray, confidence: float = 0.68) → dict[str, float][source]

Aggregate per-seed estimates into one mean with hierarchical uncertainty.

The total variance combines between-seed and within-seed components: Var(mean) = Var_between / n_seeds + mean(Var_within) / n_seeds.