API reference

Every public function and class in the four subpackages and the pygfmc.protocols module, with its docstring. The entries are generated from the source code, so they always match the installed version.

`pygfmc.protocols`

The interfaces of a Hamiltonian or off-diagonal observable, a diagonal observable and a guiding wavefunction, as typing.Protocol classes. A class satisfies one by having the members; it does not inherit from it. They are also available directly from pygfmc.

`Operator`()	The rate protocol of a Hamiltonian or of an off-diagonal observable.
`DiagonalOperator`()	An operator measured through its diagonal only (`diag_operators`).
`Guide`()	A guiding wavefunction `psi_T`.

class pygfmc.protocols.Operator

The rate protocol of a Hamiltonian or of an off-diagonal observable.

The operator is a sum H = sum_i h_i of M_total terms, each of which produces at most one new basis state from a given one. Term 0 is the diagonal (all diagonal operators are absorbed into it). An operator passed as hamiltonian or in off_diag_operators must be a JAX pytree, because it is passed into jitted functions.

M_total: int

one diagonal term plus one per possible transition.

Type:: Number of terms

get_rates(config, return_configs=False)

Matrix elements of all terms acting on config.

Returns (rates, new_configs): rates has shape (M_total,) with the diagonal element first. new_configs is the (M_total, N_sites) array of successor states (row 0 is config itself) when return_configs=True, else None.

generate_state(config, idx)

The basis state produced by term idx acting on config.

Required when M_total > 1. The propagators use it to avoid materialising all successors.

move(config, idx)

Optional. The sites changed by term idx and their increments, (sites, deltas) of a fixed shape ((2,) each for the built-in classes; zero increments for idx == 0), so that generate_state(config, idx) is config with deltas added at sites. Needed only by a guide in the incremental mode (see Guide), which updates its per-walker state from them.

class pygfmc.protocols.DiagonalOperator

An operator measured through its diagonal only (diag_operators).

The built-in Hamiltonians provide diagonal_weight as well, so they can be measured as diagonal observables.

diagonal_weight(config)

The diagonal element for config: a scalar, or a vector of values for a vector-valued observable.

class pygfmc.protocols.Guide

A guiding wavefunction psi_T.

The engine needs only the ratios psi_T(x') / psi_T(x) to the M_total successors of the current configuration, so a guide implements at least one of the three methods below; the engine uses the first one it finds, in this order. build_cache is optional. A guide that is passed into a run must be a JAX pytree.

The class declares all four methods so that each is documented once; a guide that has only some of them is accepted at run time.

log_diff_from_operator(config, op, cache=None, state=None)

log psi_T(x'_k) - log psi_T(x) for every term k of op, shape (M_total,), computed from the operator’s structure without materialising the successors. The fastest path. state is the walker’s state of the incremental mode (passed as state= only when incremental is set; see init_state()).

log_diff_from_configs(config, new_configs, cache=None)

The same log-ratios from the materialised successors new_configs of shape (M_total, N_sites).

log_psi(config, cache=None)

log psi_T(config), a scalar. The last fallback: the engine then evaluates it on every successor.

build_cache(op)

Optional. Precompute operator-specific data once per run; the result is passed as cache to log_diff_from_operator. The configuration-based fallbacks are called without a cache.

incremental: bool

the engine keeps a per-walker state from init_state(), updates it with update_state() after every move (the Hamiltonian must provide Operator.move()), carries it through reconfiguration, passes it as state= to log_diff_from_operator(), and rebuilds it every 100 chain steps so that rounding cannot accumulate. The built-in Jastrow guides keep their weighted neighbour sums this way (constructor flag incremental). A run in this mode is not bit-identical to the same run without it.

Type:: Optional. True switches on the incremental mode

init_state(config, cache=None)

Optional (incremental mode). The per-walker state of config, a pytree of arrays.

update_state(state, config, sites, deltas, cache=None)

Optional (incremental mode). state after the move that adds deltas to sites of config (the result of Operator.move()); log_diff_from_operator then receives it as state=.

`pygfmc.walkers`

The walker population, its initializers and reconfiguration.

`Walker`(config, log_weight, e_loc)	Walker population (or a single walker) as a JAX pytree.
`prepare_walkers_random`(N_walkers, N_sites, ...)	Initialize walkers by drawing every site independently from `states`.
`prepare_walkers_sz`(N_walkers, N_sites, rng_key)	Initialize walkers in a fixed-charge sector drawn from `states`.
`prepare_walkers_from_configs`(configs)	Initialize walkers at user-supplied configurations.
`reconfiguration`(walkers, rng_key)	Systematic resampling of the walker population.

class pygfmc.walkers.Walker(config, log_weight, e_loc)

Walker population (or a single walker) as a JAX pytree.

config

configurations, shape (N_walkers, N_sites) for a batch ((N_sites,) for a single walker). Spin-1/2 models use the +-1 encoding.

Type:: jax.Array

log_weight

accumulated log weight since the last reconfiguration, shape (N_walkers,).

Type:: jax.Array

e_loc

local energy e_L(x) of the current configuration, shape (N_walkers,).

Type:: jax.Array

save(path, rng_key=None, group='walkers', **attrs)

Save walker arrays + arbitrary metadata to an HDF5 file.

Walker arrays (config, log_weight, e_loc) are stored as HDF5 datasets; the JAX PRNG key (if any) is also a dataset. config is stored in the smallest lossless dtype (int8 for ±1 / small-integer configs, ~4x smaller than float32; original dtype otherwise) — see _compact_config(). load casts it back to the compute dtype. Any additional keyword arguments — typically avg_log_weight=..., but anything scalar — are stored as group attributes (same convention as JastrowSpinHalf.save()).

Parameters:

path (str) – HDF5 file path (created or overwritten).
rng_key – JAX PRNG key to checkpoint alongside the walkers (a legacy PRNGKey or a typed jax.random.key; stored as raw key data and loaded back as a legacy key). Only Threefry keys with the default Threefry PRNG are supported; other implementations raise before writing. Pass the key run_gfmc returned (stats['parameters/rng_key_next']): the key carries the chain position, so it is the whole resume state. Walker.load gives it back, and passing it as rng_key continues the stream instead of replaying it.
group (str) – HDF5 group name to write under.
**attrs – extra scalar metadata stored as group attributes.

classmethod load(path, group='walkers', dtype=None, config_dtype=None)

Load a Walker (and optional RNG key / avg_log_weight) from an HDF5 checkpoint.

Parameters:

path (str) – path to the HDF5 file written by save().
group (str) – HDF5 group name to read from.
dtype – JAX dtype for weights and local-energy estimators. None (the default) uses JAX’s default float – float32, or float64 when jax_enable_x64 is set – exactly as the prepare_walkers_* initializers do, so a loaded ensemble has the same width as a freshly prepared one.
config_dtype – Optional JAX dtype for configurations. Defaults to dtype.

Returns:

(walker, rng_key, avg_log_weight) where rng_key and avg_log_weight are None when not saved.

pygfmc.walkers.prepare_walkers_random(N_walkers, N_sites, rng_key, states=(-1.0, 1.0))

Initialize walkers by drawing every site independently from states.

Each site of each walker is an independent uniform draw, so the total magnetization/occupation fluctuates binomially about N_sites * mean(states). Use this for Hamiltonians that do NOT conserve the corresponding charge (TFIM, toric code with h_x); use prepare_walkers_sz() to start inside a fixed sector.

Parameters:

N_walkers (int) – number of walkers.
N_sites (int) – number of sites.
rng_key (Array) – JAX random key.
states (Sequence[float]) – local values to draw from – (-1, 1) for the +-1 spin-1/2 encoding (the default), (-0.5, 0.5) for physical S^z, (-1, 0, 1) for S=1, (0, ..., n_max) for bosons. Drawn with equal probability.

Note there is no dtype argument: the configuration dtype is taken from states itself, so an integer ladder such as (-1, 0, 1) yields integer configurations and (-0.5, 0.5) yields floating-point ones. The weights and local energies use JAX’s default float dtype (float32, or float64 when jax_enable_x64 is set); they never inherit an integer dtype.

Returns:: Walker batch with zero log-weights.
Return type:: Walker

pygfmc.walkers.prepare_walkers_sz(N_walkers, N_sites, rng_key, sz=0, states=(-1.0, 1.0))

Initialize walkers in a fixed-charge sector drawn from states.

sz is simply the SUM over all sites of the local values, each of which is one of states:

sz = sum_i sigma_i, sigma_i in states

so its meaning follows the encoding rather than being fixed to spin-1/2. With the default states=(-1.0, 1.0) it is the sum of the +-1 config values, i.e. twice the physical S^z_tot, and the familiar constraint that N_sites + sz be even is just the two-state case of the general reachability condition. With states=(-0.5, 0.5) the same sz is S^z_tot itself; with (-1, 0, 1) it is S^z_tot for S=1; with (0, ..., n_max) it is the total boson number.

Every walker is an independent random permutation of one base multiset, so the sum is EXACT for every walker – which is what makes this usable to pin a conserved sector: for a Hamiltonian that conserves the charge, walkers started in the wrong sector can never leave it. A request outside the reachable range, or one the ladder cannot represent, raises rather than silently returning a config with the wrong sum.

Parameters:

N_walkers (int) – number of walkers.
N_sites (int) – number of sites.
rng_key (Array) – JAX random key.
sz (float) – target SUM over all sites of the local values. Default 0.
states (Sequence[float]) – local values available per site; default (-1, 1).

As in prepare_walkers_random() there is no dtype argument: the configuration dtype follows states, while the weights and local energies use JAX’s default float dtype.

Returns:: Walker batch with zero log-weights.
Return type:: Walker

pygfmc.walkers.prepare_walkers_from_configs(configs)

Initialize walkers at user-supplied configurations.

Like prepare_walkers_random() and prepare_walkers_sz(), this takes no dtype: the configurations keep the dtype of configs, so an int8 input stays int8. Cast the input yourself – jnp.asarray(configs, dtype=...) – to get something else. The weights and local energies use JAX’s default float dtype.

Useful whenever a specific sector must be seeded exactly, e.g. the all-ones plaquette sector of the toric code (prepare_walkers_from_configs(jnp.ones((N_walkers, N_sites)))), or to start from configurations prepared elsewhere.

Parameters:: configs (Sequence | ndarray | Array) – (N_walkers, N_sites) array, one configuration per walker.
Returns:: Walker batch with zero log-weights.
Return type:: Walker

pygfmc.walkers.reconfiguration(walkers, rng_key)

Systematic resampling of the walker population.

Resamples walkers proportional to their weights using the systematic scheme (one shared random offset). After resampling all weights are reset to one (log-weight = 0).

Parameters:

walkers (Walker) – Walker batch to resample.
rng_key (Array) – JAX PRNGKey.

Returns:

(walkers, log_G_n, E_n, indices). log_G_n is the logarithm of the mean walker weight G_n = (1/N) Σ_i w_i (the step’s branch weight, kept in log space because its product over a window spans hundreds of orders of magnitude), E_n the weight-averaged local energy Σ_i w_i e_L(x_i) / Σ_i w_i, and indices[i] the parent of new walker i (the parent map is needed to gather the new configurations, so it is always available and always returned).

Return type:

tuple[Walker, Array, Array, Array]

`pygfmc.engine`

The GFMC engine: propagation, estimators and HDF5 storage. These functions are also available directly from pygfmc.

`run`(walkers, n_steps, hamiltonian, ...[, ...])	Auto-dispatch GFMC: run_gfmc_pmap on >1 devices, run_gfmc otherwise.
`run_gfmc`(walkers, n_steps, hamiltonian, ...)	Run the GFMC production loop on a single device.
`statistics`(walkers, hamiltonian, step_size, ...)	Run branching steps and collect diagnostic statistics.
`evaluate_observables`(raw[, p_array, n_blocks])	Re-evaluate the energy and the observables from raw per-step series.
`write_to_HDF5`(filename, stats[, walkers, ...])	Append one GFMC run to an HDF5 file as a self-contained group.
`run_from_HDF5`(filename, hamiltonian[, ...])	Continue a stored chain in one call, and return what `run_gfmc` returns.
`load_raw_from_HDF5`(filename)	Read a file's raw series back in the layout `run` returns them in.
`check_run_compatible`(filename[, stats])	Validate a planned run against a stored file before walking.
`resume_from_HDF5`(filename[, hamiltonian, ...])	Assemble the `run_gfmc` arguments that continue a stored chain.

pygfmc.engine.run(walkers, n_steps, hamiltonian, step_size, guide, p_array, *, method='cont_time', rng_key=None, lam=None, **kwargs)

Auto-dispatch GFMC: run_gfmc_pmap on >1 devices, run_gfmc otherwise.

The recommended entry point. Drop-in replacement for pygfmc.engine.driver.run_gfmc() — same positional arguments and keyword arguments (plus n_devices, default jax.device_count()), same return tuple. lam is needed only by the discrete methods. When more than one device is used it dispatches to run_gfmc_pmap(), which shards walkers across devices; otherwise it calls run_gfmc unchanged.

Runs with diag_operators fall back to the single-device backend. Note that the two backends differ in a few details (see run_gfmc_pmap()).

pygfmc.engine.run_gfmc(walkers, n_steps, hamiltonian, step_size, guide, p_array, *, method='cont_time', rng_key=None, lam=None, diag_operators=None, off_diag_operators=None, return_raw=True, avg_log_weight=None, warmup_steps=None, verbosity=7, number_of_outputs=100, off_diag_stride=1, off_diag_batch_size=None, _step_offset=0)

Run the GFMC production loop on a single device.

Sequence:

Warm-up: a 4-stage step-size ramp of warmup_steps total steps (skipped when warmup_steps=0). It equilibrates the ensemble and, unless avg_log_weight is supplied, also measures the normalisation shift avg_log_weight.
Production scan of n_steps branching steps (evolve -> reconfigure -> diagonal ancestry / off-diagonal forward walking), printing progress every n_steps // number_of_outputs steps.
Windowed estimation over p_array and an end-of-run pygfmc.engine.estimators.statistics() pass (sampler health, offspring distribution, ancestor survival at m = p/2).

Parameters:

walkers (Walker) – Initial Walker batch.
n_steps (int) – Number of production branching steps.
hamiltonian (Operator) – Any object implementing the Hamiltonian protocol (M_total, get_rates, generate_state) — e.g. SpinHalf, Spin, BoseHubbard, ToricCode, or PlaquetteModel.
step_size (int | float) – Per-Markov-step size. For the discrete methods this is the integer count of imaginary-time sub-steps (must be integer-valued); for cont_time it is the imaginary-time length $beta$ (any positive float).
guide (Guide | None) – Guiding wavefunction or None.
p_array (Sequence[int] | ndarray) – Window sizes for the mixed/pure estimator (deduplicated and sorted). The projection depth per side is tau = (p/2)·step_size. When an off-diagonal operator has an off-diagonal transition (M_total > 1) all p values must be even.
method (str) – “single_step”, “multi_step”, or “cont_time” (default “cont_time”). This choice sets the meaning of step_size (integer sub-step count for the discrete methods, imaginary time beta for cont_time) and of avg_log_weight (log b per sub-step vs. minus the mean local energy per unit imaginary time); under cont_time lam has no effect (it is the lambda -> infinity limit).
rng_key (Array | None) – JAX PRNGKey or typed Threefry key. None (default) draws a fresh one from secrets.randbits(32) (printed at verbosity >= 1, stored in stats['parameters/rng_key']). The default PRNG must be Threefry; unsupported key implementations raise before running. THIS IS THE CHAIN POSITION, not just a seed: the production stream starts here and advances key_{t+1} = split(key_t) once per step, with no step index anywhere. Pass back the key a previous run returned (stats['parameters/rng_key_next'], which Walker.save(rng_key=...) checkpoints and Walker.load returns) and the continuation is bit-identical to one uninterrupted run – that one array is the whole resume state. Everything else needing randomness draws from tagged side-streams, so it cannot shift the production sequence.
lam (float | None) – GFMC shift λ, required by the discrete methods, where it must satisfy λ ≥ max_x H_xx over the reachable configurations (ValueError when omitted). Unused by cont_time, the default method, so it defaults to None.
diag_operators (Mapping[str, DiagonalOperator] | None) – dict name → operator of diagonal operators, or None. Each needs diagonal_weight(config) returning a scalar or an (n_vals,) vector. Measured with the ancestry-traced (forward-walking) pure estimator: the value at the ancestor m = p/2 steps back, weighted like the energy.
off_diag_operators (Mapping[str, Operator] | None) – dict name → operator of off-diagonal operators, or None. Each must implement the operator protocol (pygfmc.Operator), with non-negative rates. Measured via the forward-walking pure estimator (one shadow chain of depth max(p)/2 per operator per measurement). Every p must be even when an operator has M_total > 1; a purely diagonal operator (M_total == 1) accepts odd p as well.
return_raw (bool) – If True (default), append the raw dict to the return tuple (see Returns). The raw series can then be stored and re-estimated over the whole concatenated chain (evaluate_observables), where windows span run boundaries, so every off-diagonal measurement is kept. With False the per-call windowed numbers are final, and the off-diagonal measurements are restricted to the band of steps usable by EVERY p (the first max(p) steps become a lead-in and the last max(p)/2 a tail), so each p rests on the identical set of measurements; n_steps > 1.5*max(p) is then REQUIRED (ValueError otherwise).
avg_log_weight (float | None) – Pre-computed normalisation shift (e.g. from stats['avg_log_weight'] of a previous run). When supplied it is kept fixed; the warm-up ramp still runs (unless warmup_steps=0) but only equilibrates the walkers. When None it is measured by the warm-up (see warmup_steps).
warmup_steps (int | None) – Total warm-up steps split across 4 ramp stages (step_size/8 → step_size/4 → step_size/2 → step_size; integer division, floored at 1, for the discrete methods). Default None: 1000 steps (DEFAULT_WARMUP_STEPS) when no avg_log_weight is given, which the warm-up then measures, and 0 when one is – right when continuing from an equilibrated ensemble with stats['avg_log_weight'] of the previous run. A fresh start of a large system may need more, e.g. warmup_steps=10000. An explicit warmup_steps=0 without avg_log_weight raises ValueError.
verbosity (int) –
Integer 0–10 controlling printed output:
- 0 — silent
- 1 — step/it/N_config progress (no energy)
- 3 — progress + Energy[p]
- 8 — progress + Energy[p] + first diag/off-diag observable (re-estimated over the whole series at every progress line, so it costs host time) + warmup lines
- 10 — full output (all observables)
The end-of-run statistics report prints for any verbosity >= 1. Default 7.
number_of_outputs (int) – How many progress lines to print during the production run (default 100). Actual cadence is n_steps // number_of_outputs.
off_diag_stride (int) –
Measure the off-diagonal operators only every off_diag_stride-th production step (default 1 = every step). Each off-diagonal measurement launches a fresh shadow Markov chain m_max = max(p)/2 steps deep, so it dominates the per-step cost by ~m_max× — but consecutive main-chain steps are strongly correlated, so per-step measurement mostly averages redundant values. A stride ≈ the observable’s autocorrelation time keeps nearly all the statistical information at ~1/stride of the cost. Only the measured steps are stored (raw['off_diag'] has one row per measurement, raw['off_diag_steps'] gives each row’s step). Energy and diagonal observables are unaffected — they are cheap and stay per-step.

The stride counts BRANCHING STEPS, so choose it against step_size: decide the imaginary time between measurements and convert:
```
off_diag_stride = round(tau_between/step_size)
```
since ensemble decorrelation is a process in imaginary time. A stride copied between runs with different step_size silently changes meaning. A useful rule of thumb is tau_between ≈ tau_max/4 (i.e. stride ≈ m_max/4), bounded below by the observable’s autocorrelation time.
off_diag_batch_size (int | None) – optional walker-chunk size for APPLYING each off-diagonal operator (not for the shadow chain, whose cost is independent of the operator’s width). None = apply to all walkers at once. An int cuts the peak (N_walkers, n_off, max_aff) intermediate by off_diag_batch_size/N_walkers and is BIT-IDENTICAL: the operator application is per-walker independent and each walker’s PRNG key comes from its absolute index. Must divide the walker count. Needed only for a WIDE operator (e.g. a pair sum over all sites). Costs sequential passes, so leave it None unless memory demands it.
_step_offset (int) – Internal. The number of steps the chain already holds; run_from_HDF5 passes it so that the fixed refresh schedule of a guide in the incremental mode counts absolute chain steps and a split run reproduces a contiguous one.

Returns:

(walkers, energy, obs, stats, raw) by default (return_raw=True); with return_raw=False the 4-tuple (walkers, energy, obs, stats).

walkers: Walker batch after the production run (the production-final state; the extra lead-in and statistics steps are not applied to it).
energy: dict p → float energy estimates.
obs: dict name → dict p → value operator estimates (a float, or an (n_vals,) array for a vector-valued diagonal operator); empty when no operators are provided.

stats is a dict of run diagnostics — the end-of-run statistics() dict (health stats, offspring_k*, ancestor_m, ancestor_fraction, effective_sample_size) plus:

stats['avg_log_weight']: the normalisation shift used in production; pass it back to reuse it.
stats['it_per_sec']: production throughput.
stats['runtime']: wall-clock seconds for the production loop.
stats['total_wall_time']: complete call time, including warm-up, final estimation, and diagnostics.
stats['number_of_devices']: 1.
stats['autocorrelation/...']: integrated autocorrelation times of the energy, log weights and every operator series, in branching steps.
stats['parameters/...']: the run’s inputs (n_steps, n_walkers, step_size, lam, method, p_array, m_array, off_diag_stride, operator names and content hashes, dtypes, …) and the RNG stream: parameters/rng_key (where production started) and parameters/rng_key_next (where it ended — pass it as rng_key to continue).

raw holds the per-step GFMC values and the metadata needed to re-estimate them (evaluate_observables(raw) needs no further arguments):

raw['log_weights']: (n_steps,) log branch weights log G_n.
raw['energy']: (n_steps,) per-step energies E_n.
raw['diag']: dict name -> (n_steps, n_p) for a scalar operator or (n_steps, n_vals, n_p) for a vector-valued one; None without diagonal operators. Operators of different widths may share a run.
raw['additional_diag']: same layout, for the m_max steps FOLLOWING this run, traced with this run’s ancestry; the lead-in a continuation substitutes for its own first rows (see run_from_HDF5).
raw['off_diag']: dict name -> (n_meas, n_p) forward-walk values, one row per MEASURED step; None without off-diagonal operators.
raw['off_diag_norm']: dict name -> (n_meas, n_p) log shadow-chain weight products (same rows).
raw['off_diag_steps']: (n_meas,) 0-based production step (within this call) at which each row was measured.
raw['p_array'], raw['m_array'], raw['diag_operators'], raw['off_diag_operators']: the window grid and operator names; the last axis of every operator series is positional in p_array.

Return type:

tuple

pygfmc.engine.statistics(walkers, hamiltonian, step_size, guide, method, rng_key, avg_log_weight=None, lam=None, M=10, quiet=False, m_array=None, cache=None)

Run branching steps and collect diagnostic statistics.

One evolution pass provides three families of diagnostics:

Sampler health (first M steps): per-step min/max of the diagonal weight, survival probability p_d, diagonal sub-step count k, unique-configuration counts before/after each reconfiguration, and the changed/overlap mixing rates.
Offspring distribution (every step): from each reconfiguration’s parent map, the fraction of walkers with exactly k children for k = 0..3 plus an overflow bin (>= 4); fractions are in [0, 1] and sum to 1. fraction[0] is the per-step lineage-death rate (= the duplicate rate sum_{k>=2}, since resampling conserves N) and scales ~linearly with step_size.
Ancestor survival (optional, when m_array is given): each walker’s step-0 ancestor label is carried through the resampling indices; at every m in m_array the fraction of distinct surviving ancestors and the Kish effective sample size N^2 / sum_i c_i^2 (c_i = number of descendants of ancestor i) are recorded — the quantity that bounds the pure estimator’s effective sample size at projection depth m = p/2. The same ancestry gives the lineage correlation of the local energy: the correlation between a walker’s e_L at depth m and its own step-0 ancestor’s, which measures in branching steps how long the population remembers its state. The run is extended to max(M, max(m_array)) steps; steps beyond M skip the (expensive) sampler-health diagnostics.

Start from a stored, equilibrated population for meaningful ancestor curves (with near-uniform weights systematic resampling is the identity map, so the first steps of a fresh population show fraction 1.0).

Parameters:

walkers (Walker) – Walker batch to use as the starting point.
hamiltonian (Operator) – Hamiltonian object.
step_size (int | float) – Per-Markov-step size; forwarded to evolve_walkers (integer n_bra for the discrete methods, float beta for cont_time).
guide (Guide | None) – Guiding wavefunction or None.
method (str) – “single_step”, “multi_step”, or “cont_time”.
rng_key (Array) – JAX PRNGKey.
avg_log_weight (float | None) – Normalisation shift per sub-step (discrete methods) or per unit imaginary time (cont_time); evolve_walkers subtracts it * step_size. None (default) uses estimate_average_log_branch_weight() on the given walkers. Reconfiguration normalises the weights, so every diagnostic here is invariant to this shift up to float rounding; it only keeps the log weights near zero. Pass stats['avg_log_weight'] of a run to reproduce that run’s arithmetic exactly.
lam (float | None) – GFMC shift λ; required by the discrete methods, unused by cont_time.
M (int) – Number of full-diagnostic branching steps (default 10).
quiet (bool) – If True, suppress printed output.
m_array (Sequence[int] | None) – Optional list of ancestor-survival depths (in steps; tau = m * step_size). None (default) disables the lineage tracking and keeps the run at M steps.
cache (Any) – Optional guide.build_cache(hamiltonian) cache; computed internally when None.

The diagnostics run on the full population on one device. A population that only fits when sharded over several devices can exhaust memory here; pass a slice of it (jax.tree_util.tree_map(lambda x: x[:n], walkers)) in that case, noting that effective_sample_size then describes the slice.

Returns:

walkers: Walker batch after max(M, max(m_array or [0])) steps.
stats: dict with the keys below.
n_unique_pre (M,): unique configs before each reconfiguration.
n_unique_post (M,): unique configs after each reconfiguration.
unique_fraction (M,): n_unique_post / N_walkers – the walker diversity rho plotted against the step size in the paper’s convergence figure. Around 0.95 is a healthy value.
survival_rate (n_total,): 1 - offspring_k0, i.e. the fraction of walkers that produce at least one offspring. NOTE this is not the same as unique_fraction: it counts surviving parent slots, not distinct configurations, so two survivors carrying the same config count twice. The two agree only when the configuration space is large compared to N_walkers (they differ by a factor of ~3 for an L=8 chain with 256 walkers).
diag_min (M,): per-step minimum diagonal weight.
diag_max (M,): per-step maximum diagonal weight.
p_d_min / p_d_max / p_d_mean / hold_steps_mean / hold_steps_max (M,): discrete methods only – the stay probability and the sampled geometric hold length in sub-steps.
escape_rate_min / escape_rate_max / escape_rate_mean / hold_time_mean / hold_time_max / hops_per_step_mean / hops_per_step_max (M,): cont_time only – the guide-transformed escape rate Gamma, the sampled exponential hold time, and Gamma*step_size = the average number of configuration updates per branching step (see the block that sets them).
changed_rate (M,): per-step fraction of walkers whose config changed from step 0.
overlap_rate (M,): per-step fraction of step-t configs present in the step-0 population.
H_diag_max_observed / H_diag_min_observed (1,): scalar extremes of the diagonal matrix element H_diag(x) over all walkers and all M health steps (diag_min/max are the raw rates[0] extremes; the max is the number lambda tuning needs: choose lam >= it + margin).
offspring_k<j> (n_total,): fraction of walkers with exactly j children per step, j = 0..3.
offspring_k4plus (n_total,): fraction with >= 4 children.
ancestor_m / ancestor_fraction / effective_sample_size (len(m_array),): only when m_array is given. ancestor_fraction is a FRACTION in (0,1]: #{c_i > 0} / N. effective_sample_size is a COUNT in [1, N]: Kish’s ESS, n_eff = N^2 / sum_i c_i^2 (equivalently (sum_i c_i)^2 / sum_i c_i^2, since sum_i c_i = N), = 1 / P(two random walkers share an ancestor). At m = 0 it is N.
ancestor_correlation (len(m_array),): only when m_array is given. Pearson correlation between the local energy of each walker at depth m and that of its step-0 ancestor (every walker paired with its own ancestor, so lineages with many descendants weigh more). 1 at m = 0, decaying towards 0 over the population’s decorrelation time; the per-walker counterpart of the run’s autocorrelation/energy_estimator. NaN when every walker has the same local energy.

Return type:

(walkers, stats)

pygfmc.engine.evaluate_observables(raw, p_array=None, n_blocks=None)

Re-evaluate the energy and the observables from raw per-step series.

Takes the raw dict run/run_gfmc return, or the one load_raw_from_HDF5() reads back from a file; both carry the p grid (raw['p_array'], raw['m_array']) and the operator names alongside the series, so the call needs nothing else. A raw assembled by hand – an analysis script slicing the series – may carry neither, and then p_array is required: the operator arrays are indexed by POSITION in the grid and nothing in them says what it was. Operators of such a raw are named by their mapping keys, or positionally (diag_0, off_diag_0, …) for dense arrays.

The windowed values a run returns are fixed at the grid that run used. The raw series are not, and because a file holds one continuous chain they can be re-estimated over the whole of it – across run boundaries, which is legitimate here precisely because the boundaries are joins in one chain rather than seams between independent ones.

The energy takes any p. The observables do not: their ancestry depth and forward-walking shadow chains are fixed while the chain is walking, so they are only ever available at the grid they were measured on, and asking for other values returns the energy there and the observables at that grid.

Passing n_blocks additionally returns statistical errors from a blocked jackknife: the chain is cut into n_blocks contiguous blocks, each is deleted in turn and the whole estimator recomputed on the rest, and the error is the inflated scatter of those replicas, sqrt((nb-1)/nb * sum (x_b - mean)^2). Recomputing rather than averaging per-block values is what makes it correct for a ratio estimator like sum G*E / sum G, and it reuses _estimate_observables so the replicas inherit its window bookkeeping unchanged. If the sample count is not divisible by n_blocks, trailing samples are discarded with a warning, from BOTH the point estimates and the error calculation.

The block length is what matters, and it must exceed both the autocorrelation time and the window length: deleting one step from a 20000-step chain with tau_int = 9 underestimated the error 4-fold, and the estimate only settled once blocks were longer than ~4p. Choose n_blocks so that the block length n_steps / n_blocks clears max(2*tau_int, p) comfortably, and keep at least ~20 blocks or the error on the error takes over. Two effects push the same way here: a replica is the surviving steps concatenated, so windows spanning the cut combine weights from either side of it, which understates the error by ~p/blen (measured 19 % at blen = 2p, under 5 % at blen = 10p). Clearing the block-length condition suppresses that along with the autocorrelation, and at the plateau it is below the ~1/sqrt(2*n_blocks) uncertainty of the error estimate itself. The reliable way to choose it is to scan n_blocks and take the value where the error stops growing: on the chain above it rose 8.4e-3, 8.6e-3, 8.9e-3 for block lengths 40, 80, 160 and settled there. tau_int for each series is already in the file under run_<n>/statistics/autocorrelation/.

Note that this is the statistical error only. It says nothing about the systematic bias of a finite p, which is the larger danger at deep projection: an over-projected estimate can be far off while carrying a small error bar, because the surviving lineages agree with each other. Independent runs across several p remain the only check for that.

Parameters:

raw (Mapping[str, Any]) – raw dict from run/run_gfmc, or from load_raw_from_HDF5(). energy and log_weights are required; diag/off_diag and their companions are used when present.
p_array (Sequence[int] | ndarray | None) – Window sizes for the energy. Defaults to the grid the run used (raw['p_array']), in which case energy and observables share one grid; any other grid re-windows the energy alone. For a hand-built raw without p_array it is also the grid the operator arrays are taken to sit on.
n_blocks (int | None) – Number of jackknife blocks. None skips the error estimate and returns the point estimates alone.

Returns:

(energy, obs) in run_gfmc’s own form – dict p -> energy and dict name -> dict p -> value – or (energy, obs, energy_err, obs_err) with the same key structure when n_blocks is given. obs is empty when raw holds no operator series.

Raises:

ValueError – when raw holds no energy or log-weight series, when no p grid is available from either p_array or raw, when the operator arrays do not match the grid they are said to be on, or when n_blocks leaves fewer than 2 blocks.

Return type:

tuple

pygfmc.engine.write_to_HDF5(filename, stats, walkers=None, energy=None, obs=None, raw=None, walker_filename=None, verbose=True)

Append one GFMC run to an HDF5 file as a self-contained group.

Creates the file if it does not exist, otherwise appends the next /run_<n> (0-based). filename and stats are required; every other argument is optional. MUST_MATCH parameters are validated against the file’s /parameters group (and the RNG stream against the stored runs: a replayed key raises, a non-continuing key warns) before anything is written, and the write is skipped entirely if validation fails.

Parameters:

filename (str) – Path to the HDF5 file.
stats (Mapping[str, Any]) – stats dict from run_gfmc. Required. Its parameters/ entries carry everything that identifies the run, so without them a file has nothing to validate an append against and nothing to rebuild a continuation from — the two things the format exists for. The fixed ones go to /parameters, the per-run ones to run_<n>/parameters, and everything measured (autocorrelation/*, it_per_sec, runtime, total_wall_time, the sampler-health and lineage series, avg_log_weight) to run_<n>/statistics.
walkers (Walker | None) – Final Walker batch. Stored as the top-level /current_walkers checkpoint, overwriting the previous one — only the newest ensemble can continue the chain.
energy (Mapping[int, float] | None) – dict p → energy, as returned by run_gfmc.
obs (Mapping[str, Mapping[int, Any]] | None) – dict name → dict p → value.
raw (Mapping[str, Any] | None) – raw dict from run_gfmc. Appended to the top-level /raw series. Storing it is what makes the run re-estimable afterwards (via load_raw_from_HDF5() and evaluate_observables()); without it only the windowed values for this run’s p_array survive.
walker_filename (str | None) – Write the ensemble to this separate file instead of embedding it, and record the path. Useful when the ensemble dwarfs the measurements or lives on faster storage.
verbose (bool) – Print one line naming the group written (default True); False keeps the call silent.

Returns:

The index n of the group just written.

Return type:

int

pygfmc.engine.run_from_HDF5(filename, hamiltonian, guide=None, diag_operators=None, off_diag_operators=None, n_steps=None, walker_filename=None, **run_kwargs)

Continue a stored chain in one call, and return what run_gfmc returns.

Everything the continuation needs already sits in the file — the fixed parameters in /parameters, the ensemble in /current_walkers, the stream in the newest run’s rng_key_next. This reads all of it, verifies the objects handed in against the stored fingerprints, and runs:

walkers, energy, obs, stats, raw = run_from_HDF5(
    "out.h5", H, guide, diag_ops, off_diag_ops)
write_to_HDF5("out.h5", stats, walkers, energy, obs, raw)

No warmup is performed. The stored ensemble is already equilibrated and its avg_log_weight is passed through, so re-warming it would discard that equilibration and leave a discontinuity mid-chain.

The diagonal observables come out bit-identical to an uninterrupted run. A fresh ancestry map cannot trace back past its own start, so the first p/2 rows would otherwise be zeros; the previous run recorded exactly those values by walking m_max steps further while it still had the history, and they are substituted here. The off-diagonal observables need no such repair — their shadow chains only look forward.

lam, step_size, p_array, method and off_diag_stride come from the file rather than from arguments, which is what makes the continuation the same calculation instead of a similar one. Anything in MAY_DIFFER may still be passed through run_kwargs — verbosity and number_of_outputs typically are.

Parameters:

filename (str) – File written by write_to_HDF5().
hamiltonian (Operator) – The Hamiltonian. Compared against the stored class name and fingerprint.
guide (Guide | None) – Guiding wavefunction, compared the same way.
diag_operators (Mapping[str, DiagonalOperator] | None) – Diagonal operator dict, compared by name, order and fingerprint.
off_diag_operators (Mapping[str, Operator] | None) – As diag_operators.
n_steps (int | None) – Steps to run. Defaults to the newest stored run’s n_steps, so a self-looping campaign repeats equal blocks without restating the length.
walker_filename (str | None) – Read the ensemble from this file instead of the embedded checkpoint.
**run_kwargs – Passed to run_gfmc, overriding what was read from the file. Raises on any MUST_MATCH name, since overriding one of those is precisely what the file is there to prevent.

Returns:

(walkers, energy, obs, stats, raw) — exactly run’s tuple (a 4-tuple without raw if return_raw=False is passed).

Raises:

ValueError – on a parameter or fingerprint mismatch, an RNG stream replay, a missing ensemble, or a MUST_MATCH override.

Return type:

tuple

pygfmc.engine.load_raw_from_HDF5(filename)

Read a file’s raw series back in the layout run returns them in.

A file written by write_to_HDF5() holds ONE continuous chain: every run appends to the same /raw group, so what comes back here is the whole chain and not any single run’s slice. That is what makes re-estimating over it legitimate – the run boundaries are joins in one chain rather than seams between independent ones. run_<n>’s raw_step_begin/raw_step_end attributes record which rows each run contributed, for a caller that does want one run’s slice.

The keys are the ones run puts in raw – log_weights, energy, diag and off_diag/off_diag_norm (mappings keyed by operator name when the names are recorded) and off_diag_steps, so the result goes straight into evaluate_observables(). Keys the file does not hold are absent rather than None.

The metadata p_array, m_array, diag_operators and off_diag_operators are read from the file’s /parameters group, as in the raw that run returns. They matter because the operator arrays are indexed by POSITION in the p grid – off_diag[name][:, j] is the value at window p_array[j], and nothing in the array itself says so.

Parameters:: filename (str) – File written by write_to_HDF5().
Returns:: dict of the raw series, plus the metadata named above.
Raises:: ValueError – when the file stores no raw series, i.e. it was written without raw=... and there is nothing to re-estimate from.
Return type:: dict[str, Any]

pygfmc.engine.check_run_compatible(filename, stats=None, **parameters)

Validate a planned run against a stored file before walking.

write_to_HDF5 performs the same check, but only once the run has finished — by which point a mismatch has already cost the wall-clock. Call this first with the parameters you are about to use and the job stops immediately instead:

check_run_compatible(path, n_walkers=N, lam=lam, step_size=8,
                     method="multi_step", p_array=[4, 8, 16],
                     hamiltonian=H, guide=guide,
                     diag_operators=diag_ops,
                     off_diag_operators=off_ops)

hamiltonian, guide and the two operator dicts may be passed as the objects themselves, in which case their class names, operator names and content fingerprints (see pytree_fingerprint()) are all derived and compared — so a retrained guide or an edited coupling is caught even though every name still agrees. Names are checked, not only counts — two runs measuring the same number of different operators produce datasets of identical shape whose rows mean different things, which no later shape check can catch.

Accepted either as keyword arguments (bare names, no parameters/ prefix) or as a stats dict from a previous run_gfmc. Unknown and MAY_DIFFER names are ignored, so passing a whole parameter set is fine. A missing file is compatible with anything.

Raises:: ValueError – on a MUST_MATCH mismatch or an RNG stream replay.

pygfmc.engine.resume_from_HDF5(filename, hamiltonian=None, guide=None, diag_operators=None, off_diag_operators=None, walker_filename=None, dtype=None, config_dtype=None)

Assemble the run_gfmc arguments that continue a stored chain.

run_from_HDF5() does the whole continuation in one call; use this lower-level function when you want the pieces, e.g. to pass the stored ensemble and key to a run with different operators.

Reads the newest run’s embedded or external checkpoint, validates its walker hash when recorded, and returns it with that run’s parameters. An explicit walker_filename overrides the stored source but must still match the newest ensemble. Returns keyword arguments:

kw = resume_from_HDF5("out.h5")
walkers, E, obs, stats, raw = run_gfmc(
    kw.pop("walkers"), n_steps, H, step_size, guide, p_array,
    method=method, lam=lam, **kw)      # the stored run's settings
write_to_HDF5("out.h5", stats, walkers, E, obs, raw)

The returned rng_key is the stored production-final key, so the step keys continue past everything already walked, and warmup_steps=0 because the stored ensemble is already equilibrated — re-warming it would discard that equilibration and put a discontinuity in the middle of the chain, which is exactly what a single-chain-per-file layout is meant to prevent. The stored avg_log_weight is passed along for the same reason.

Pass the objects the continuation will use and they are verified against the stored run before anything is returned, so a mismatch costs nothing:

kw = resume_from_HDF5("out.h5", hamiltonian=H, guide=guide,
                     diag_operators=diag_ops)

This is the natural place for that check. write_to_HDF5 performs it as well, but only once the run has finished. Omitting the objects leaves the contents unverified and warns as much — the file records what the earlier runs used, and only the caller holds what this one is about to use.

Parameters:

filename (str) – File written by write_to_HDF5().
hamiltonian (Operator | None) – Hamiltonian the continuation will use. Verified by class name and content fingerprint against the stored run when given.
guide (Guide | None) – Guiding wavefunction, verified the same way.
diag_operators (Mapping[str, DiagonalOperator] | None) – Diagonal operator dict, verified by name, order and fingerprint.
off_diag_operators (Mapping[str, Operator] | None) – As diag_operators.
walker_filename (str | None) – Read the ensemble from this file instead of the embedded checkpoint (defaults to the path recorded at write time, if any).
dtype (Any) – JAX dtype for weights/local energies. Default: the weight_dtype recorded in the file.
config_dtype (Any) – JAX dtype for configurations. Default: the recorded config_dtype.

Returns:

walkers, avg_log_weight, rng_key (the stored production-final key, so the continuation resumes the stream) and warmup_steps.

Return type:

dict of run_gfmc keyword arguments

Building blocks

The functions the drivers are made of, for custom production loops and diagnostics. They are imported from their modules, e.g. from pygfmc.engine.evolution import evolve_walkers.

`evolve_walkers`(walkers, hamiltonian, ...[, ...])	Evolve the full walker batch (vmapped over walkers).
`evolve_carried`(walkers, carried, ...[, lam, ...])	`evolve_walkers()` with the carried pair `(rates, state)` in and out: `(walkers, carried)`.
`forward_walk_operator`(walkers, op, h_main, ...)	Forward-walking pure estimator for a single (off-diagonal) operator.
`require_lam`(method, lam)	Return `lam` validated for `method`: a float for the discrete methods, which need the shift λ, and `None` for `cont_time`, which does not use it.
`require_operator_protocol`(op[, role])	Raise `TypeError` unless `op` implements `pygfmc.Operator`.
`require_guide_protocol`(guide)	Raise `TypeError` unless `guide` is None or implements `pygfmc.Guide`, i.e. at least one of `log_diff_from_operator`, `log_diff_from_configs` and `log_psi`.

pygfmc.engine.evolution.evolve_walkers(walkers, hamiltonian, step_size, guide, method, rng_key, avg_log_weight, lam=None, cache=None, rates=None, state=None, return_rates=False)

Evolve the full walker batch (vmapped over walkers).

Parameters:

walkers (Walker) – Walker batch, shape (N_walkers, …)
hamiltonian (Operator) – Hamiltonian
step_size (int | float) – Per-Markov-step size. Interpreted as the number of imaginary-time sub-steps for the discrete methods single_step and multi_step (must be a positive integer-valued float), and as the imaginary-time length $beta$ for the continuous method cont_time (any positive float).
guide (Guide | None) – guiding wavefunction or None
method (str) – "single_step", "multi_step", or "cont_time".
rng_key (Array) – JAX PRNGKey
avg_log_weight (float) – normalisation shift subtracted from log_weight (per substep for the discrete methods, per unit imaginary time for cont_time).
lam (float | None) – GFMC shift λ. Required by the discrete methods, unused by cont_time (the default None).
cache (Any) – Optional guide.build_cache(hamiltonian) result, forwarded to the guide’s ratio routine.
rates (Array | None) – Optional (N_walkers, M_total) guide-transformed rates of the current configurations, as returned by the previous step with return_rates; the first rate evaluation of every walker is then skipped.
state (Any) – Optional per-walker guide state of the incremental mode (see pygfmc.Guide), updated after every move.
return_rates (bool) – Also return the rates and the guide state of the final configurations.

Returns:

Updated Walker batch with new config, log_weight, e_loc; with return_rates the tuple (walkers, rates, state).

Return type:

Any

pygfmc.engine.evolution.evolve_carried(walkers, carried, hamiltonian, step_size, guide, method, rng_key, avg_log_weight, lam=None, cache=None): evolve_walkers() with the carried pair (rates, state) in and out: (walkers, carried). With nothing carried the plain propagator runs and the final rates are not materialised (see _current_rates()); an entry that was not carried in stays None, so a scan carry keeps its pytree structure. The caller gathers the result with the resampling indices (gather_carried()).

pygfmc.engine.evolution.forward_walk_operator(walkers, op, h_main, step_size, guide, method, rng_key, avg_log_weight, m_array, lam=None, cache_h=None, cache_op=None, shadow_key=None, op_batch_size=None)

Forward-walking pure estimator for a single (off-diagonal) operator.

op may be any object implementing the operator protocol (pygfmc.Operator: M_total, get_rates, and generate_state when M_total > 1) — e.g. SpinHalf, PlaquetteModel, or ToricCode instances used as observables — whose guide-transformed rates are all non-negative. run_gfmc calls it at every measured production step (every off_diag_stride-th step), after reconfiguration. The drivers require even p (m = p // 2 shadow steps per window) for every operator with M_total > 1; a purely diagonal operator (M_total == 1) may use odd p, and at m = 0 its value is the plain mean of the diagonal elements (log_norm = 0).

Apply op to the current walker ensemble:

each walker x_i picks a new config x’_i proportional to the operator transition weights w(x_i, x’)

vals[i] = Σ_{x’} w(x_i, x’) (total operator weight, ≥ 0)

Propagate the operator-applied ensemble along a shadow chain for max(m_array) steps with h_main. After every step the shadow ensemble is reconfigured, vals is remapped through the parent map, and the log of the reconfiguration weight G_t is accumulated.

Record, at each m in m_array, mean(vals) and the accumulated log weight product.

Parameters:

walkers (Walker) – post-reconfiguration Walker batch (uniform weights)
op (Operator) – operator to measure (see above)
h_main (Operator) – main Hamiltonian (governs forward evolution)
step_size (int | float) – per-Markov-step size (integer sub-step count for the discrete methods, imaginary time beta for cont_time); forwarded to evolve_walkers.
guide (Guide | None) – guiding wavefunction or None
method (str) – “single_step”, “multi_step”, or “cont_time” (forwarded to evolve_walkers for the shadow chain)
rng_key (Array) – JAX PRNGKey for the operator application (and for the shadow chain when shadow_key is None)
avg_log_weight (float) – log branch-weight normalisation scalar
m_array (Array) – list of shadow-chain depths (one per p value, m = p // 2; m = 0 reads the un-propagated ensemble)
lam (float | None) – GFMC shift for h_main; required by the discrete methods, unused by cont_time
cache_h (Any) – optional guide.build_cache(h_main) result
cache_op (Any) – optional guide.build_cache(op) result
shadow_key (Array | None) – optional PRNG key to start the shadow chain from. run_gfmc passes the main chain’s NEXT production key: the shadow runs the same key recurrence, so from the (near-identical) operator-applied ensemble it replays the main chain’s future draws, and the population-level weight noise largely cancels between the estimator’s numerator and denominator (correlated sampling). Keys are state-independent, so expectations are unchanged. None derives an independent stream from rng_key (same expectation, larger variance at deep windows).
op_batch_size (int | None) – optional walker-chunk size for the OPERATOR APPLICATION (step 1). None (default) applies the operator to all walkers at once. An int splits the walker batch into chunks of that size via jax.lax.map, cutting the peak (N_walkers, n_off, max_aff) guide intermediate by op_batch_size/N_walkers. Bit-identical — the per-walker work and its PRNG key are unchanged, so no float reassociation is introduced. Must divide N_walkers. Use it when a wide operator (e.g. a site sum with many terms) would otherwise not fit in memory; it costs sequential passes.

Returns:

pair of jnp.ndarray, each shape (len(m_array),), ordered by sorted m.

val[j] = mean(vals after m_array[j] shadow steps) — the
operator measurement alone, linear and SIGNED.

log_norm[j] = log of the shadow chain’s accumulated weight
product after m_array[j] steps. It drifts roughly linearly with m and can reach -100 and beyond at deep windows, which is why it is kept in log space.

The two channels are kept separate so that the large, compensating factors of the estimator cancel in the exponent before any exp() (a single linear product underflows float32 at deep windows), while the sign lives only in val (a diagonal operator routed through this pipeline has a signed rates[0]). _estimate_observables consumes the pair at measurement step k with m = p/2 as

numerator term = G_main^{-m} · G_aux^{+m} · val denominator term = G_main^{-m} · G_main^{+m}

i.e. log_norm = log G_aux^{+m} is a numerator FACTOR: the ket half of the window is supplied by the shadow chain because that half carries the operator, while the denominator is main chain only. G_aux does not cancel — its ratio to G_main^{+m} is what converts the mixed estimator into the pure one.

Return type:

(val, log_norm)

pygfmc.engine.evolution.require_lam(method, lam)

Return lam validated for method: a float for the discrete methods, which need the shift λ, and None for cont_time, which does not use it. Raises ValueError when a discrete method is called without a shift.

pygfmc.engine.evolution.require_operator_protocol(op, role='hamiltonian')

Raise TypeError unless op implements pygfmc.Operator.

generate_state is required for any operator that can actually produce a transition (M_total > 1). It is NOT required at M_total == 1: such an operator is purely diagonal, the propagators short-circuit it, and there is no successor state to generate – a measurement-only identity is the canonical example.

Parameters:

op (Operator) – Hamiltonian or off-diagonal observable to check.
role (str) – How to name it in the error message.

pygfmc.engine.evolution.require_guide_protocol(guide)

Raise TypeError unless guide is None or implements pygfmc.Guide, i.e. at least one of log_diff_from_operator, log_diff_from_configs and log_psi.

`estimate_average_log_branch_weight`(walkers, ...)	Returns a scalar estimate of the average log branch weight across walkers.
`check_off_diag_operator_positivity`(walkers, ...)	Check off-diagonal operator weights are non-negative.
`count_unique_configs`(walkers)	Number of distinct configurations in the walker population (on the host).
`autocorrelation`(series[, c])	Integrated autocorrelation time of a per-step series.

pygfmc.engine.estimators.estimate_average_log_branch_weight(walkers, hamiltonian, guide, method='cont_time', lam=None, cache=None)

Returns a scalar estimate of the average log branch weight across walkers.

Used as the starting value of avg_log_weight before the warm-up ramp.

For the discrete methods (single_step, multi_step; lam required):: returns mean( log(lam - sum(rates)) ) — the average log branch weight per step.
For cont_time (the default; lam unused):: returns mean( -sum(rates) ) — the average negative local energy per unit imaginary time, which is the correct per-unit-time normalisation for the continuous-time propagator (avg_log_weight * beta is subtracted from each walker’s accumulated log-weight). Using the discrete formula here would give a value ~log(lam - E_diag) ≪ -E_L, causing float32 overflow in the display and a bad initial estimate that the warmup ramp would have to correct.

pygfmc.engine.estimators.check_off_diag_operator_positivity(walkers, off_diag_operators, guide, caches=None, verbose=True)

Check off-diagonal operator weights are non-negative.

M_total == 1 (diagonal only): always passes, any sign allowed. M_total > 1: all guide-transformed rates must be >= 0 (the estimator samples the successor by inverse CDF over the rates).

Parameters:

walkers (Walker) – Walker batch; its first 100 configurations are tested.
off_diag_operators (Mapping[str, Operator]) – dict name -> operator.
guide (Guide | None) – guiding wavefunction or None.
caches (Mapping[str, Any] | None) – optional dict name -> guide.build_cache(op).
verbose (bool) – print a one-line report when the check passes (the drivers pass verbosity >= 1).

Raises:

ValueError – if any rate is negative

pygfmc.engine.estimators.count_unique_configs(walkers)

Number of distinct configurations in the walker population (on the host).

pygfmc.engine.estimators.autocorrelation(series, c=5.0)

Integrated autocorrelation time of a per-step series.

Returns (tau_int, window). tau_int converts steps into independent samples – N_eff = N / (2 tau_int) – and inflates the error on a mean by sqrt(2 tau_int).

rho(t) is estimated by FFT in O(N log N). The self-consistent window is the smallest M with M >= c * (0.5 + sum(abs(rho[1:M+1]))). The absolute correlation envelope prevents premature termination for oscillating correlations; the returned time uses the SIGNED sum. Anticorrelated series can legitimately have 0 < tau_int < 0.5.

tau_int is NaN for insufficient/constant/nonfinite data, when no window is found, or when the truncated signed sum is nonpositive. These cases do not support a reliable effective-sample-count estimate.

NOTE this is the tau for statistics and budgeting. The initial transient decays with the SLOWEST MODE instead, which is longer, so it is the wrong number for choosing a warm-up length.

`run_gfmc_pmap`(walkers, n_steps, hamiltonian, ...)	Multi-device GFMC main loop.
`reconfiguration_global`(walkers, rng_key[, ...])	Reconfiguration that resamples globally across all devices.
`shard_walkers`(walkers[, n_devices])	Reshape walker arrays for pmap: (N, ...) → (n_devices, N//n_devices, ...).

pygfmc.engine.multi_gpu.run_gfmc_pmap(walkers, n_steps, hamiltonian, step_size, guide, p_array, *, method='cont_time', rng_key=None, lam=None, diag_operators=None, off_diag_operators=None, avg_log_weight=None, warmup_steps=None, verbosity=7, number_of_outputs=100, n_devices=None, return_raw=True, off_diag_stride=1, off_diag_batch_size=None, _step_offset=0)

Multi-device GFMC main loop. Energy + off-diagonal operators.

Args mirror pygfmc.engine.driver.run_gfmc(), with the same defaults, so that run() behaves the same whatever the device count. In addition:

Parameters:: n_devices (int | None) – number of devices to shard walkers across (default: all visible devices). The walker count must be divisible by it.

Differences from run_gfmc:

diag_operators raise NotImplementedError (ancestry tracing is single-device only).
Off-diagonal measurement gating (off_diag_stride > 1 and the equal windows implied by return_raw=False) works as in run_gfmc: the same steps are measured, unmeasured steps skip the shadow chain, and raw['off_diag_steps'] gives each row’s step.
Progress output: a progress line is printed from the first chunk on (run_gfmc starts once a full progress window has been walked), and the per-operator lines of verbosity >= 8 are not printed.

The warm-up follows run_gfmc: the ramp runs whenever warmup_steps > 0; a supplied avg_log_weight is kept fixed (the ramp then only equilibrates the walkers). The default warmup_steps=None runs 1000 warm-up steps without an avg_log_weight and none with one; an explicit warmup_steps=0 without avg_log_weight raises. stats carries the same parameters/* entries as run_gfmc, including parameters/rng_key_next.

Off-diagonal operators use the forward-walking pure estimator with global reconfiguration of the shadow ensemble.

rng_key is the production stream’s starting position, as in run_gfmc: keys advance key_{t+1} = split(key_t) in the scan carry and the warm-up, positivity check and statistics draw from side-streams, so the shadow chain replays the main chain’s own future draws (correlated sampling) and passing stats['parameters/rng_key_next'] back as rng_key (with the returned walkers, warmup_steps=0 and the returned avg_log_weight) continues the chain bit-identically on the same device count. The chain is NOT the single-device chain at equal seed: the evolve key carries a fold_in(device_id) and reconfiguration is global over all shards, so sharding changes the trajectory — as a different seed would, not as a bias.

stats carries the same pygfmc.engine.estimators.statistics() fields as run_gfmc — sampler health, the offspring distribution and ancestor survival / effective_sample_size — so a consumer of run() sees one contract whatever the device count. That pass runs on the UNSHARDED ensemble and therefore on a single device: a population that only fits when sharded can OOM there, after production has finished.

Returns the same tuple shape as pygfmc.engine.driver.run_gfmc(): four values when return_raw=False – stats always comes back – and five (adding raw) when return_raw=True.

pygfmc.engine.multi_gpu.reconfiguration_global(walkers, rng_key, axis_name='devices', aux=None)

Reconfiguration that resamples globally across all devices.

Each device contributes its walker shard. Weights are gathered globally, a single systematic resampling is performed, and each device picks the slice of resampled walkers it owns.

The rng_key must be identical on all devices so they agree on the random offset of the systematic resample.

If aux is given (a pytree of arrays of shape (N_per_device, …)), it is also gathered, remapped through the resampled indices, and the per-device slice is returned as a fourth element. Useful for carrying per-walker auxiliary state through reconfiguration (forward-walking values, the carried rates).

Must be called inside jax.pmap with the matching axis_name.

Returns:: (walkers_shard, log_G_n, E_n[, aux_shard]) with log_G_n the log of the global mean weight and E_n the global weighted local energy, as in pygfmc.walkers.reconfiguration().
Return type:: tuple

pygfmc.engine.multi_gpu.shard_walkers(walkers, n_devices=None)

Reshape walker arrays for pmap: (N, …) → (n_devices, N//n_devices, …).

Each leaf array gains a leading device axis. N must be divisible by n_devices.

`pygfmc.hamiltonian`

Hamiltonians, lattices, measurement operators and exact diagonalization.

`SpinHalf`()	Spin-1/2 XYZ model with longitudinal and transverse fields.
`Spin`()	Spin-S Hamiltonian; see the module docstring for H and the rate layout.
`BoseHubbard`()	Bose-Hubbard model; see the module docstring for H and the rate layout.
`ToricCode`()	Toric code Hamiltonian on an Lx × Ly torus.
`PlaquetteModel`()	Plaquette ring-exchange / Rokhsar–Kivelson Hamiltonian.
`StructureFactor`(positions, k[, channel, dtype])	The k-resolved correlator `AB(0) + AB(k)` as ONE `SpinHalf`.
`SzSzStructureFactor`(positions, k_vectors[, ...])	$S^{zz}(\vec k)$ for one or many wavevectors, from positions.
`is_valid_k`(k_vectors, lattice_vectors[, atol])	Whether each wavevector lies in the reciprocal lattice of the system.
`Lattice`(name, N, edges[, J_x, J_y, J_z, extras])	A finite lattice: a site count, an edge list and optional metadata.
`chain`(n_sites[, periodic])	Nearest-neighbour chain of `n_sites` sites.
`pyrochlore`(nx, ny, nz[, periodic])	Pyrochlore lattice built from explicit FCC unit-cell edge templates (O(N)).
`kagome`(Lx, Ly[, periodic])	Kagome lattice of `Lx` x `Ly` unit cells (O(N)).
`exact_diagonalization`(hamiltonian[, ...])	Lowest eigenpairs of `hamiltonian`, with observables on each.

Hamiltonians

class pygfmc.hamiltonian.SpinHalf

Spin-1/2 XYZ model with longitudinal and transverse fields.

H = Σ_<ij> [J_x S^x_i S^x_j + J_y S^y_i S^y_j + J_z S^z_i S^z_j]
− Σ_i [h_z S^z_i + h_x S^x_i]

The fields enter with a MINUS sign (the same convention as ToricCode), so a positive h_x is the sign-problem-free one.

The exchange term is rewritten per edge as J_pm (S+_i S-_j + S-_i S+_j) + J_pmpm (S+_i S+_j + S-_i S-_j) + J_z S^z_i S^z_j with J_pm = (J_x + J_y)/4 and J_pmpm = (J_x - J_y)/4. Configurations are ±1 floats, σ_i = 2 S^z_i. All couplings may be uniform or per-edge / per-site. Build instances with create().

Layout of the rates returned by get_rates():

rates[0]                 : diagonal H_{σσ}
rates[1 : 1+n_hx]        : single-site flips, −h_x/2 each (n_hx = N if h_x != 0)
rates[1+n_hx : M_total]  : one slot per edge, J_pm if the two spins are
                           anti-parallel, J_pmpm if parallel
                           (n_jpm = n_edges if J_pm or J_pmpm != 0)

Sign-problem-free convention: J_pm ≤ 0, J_pmpm ≤ 0 and h_x ≥ 0; J_z and h_z are unrestricted. create enforces all three by raising (never by swapping couplings). An off-diagonal measurement operator needs the opposite sign (rates ≥ 0) and is built with sign_check=False: for +Σ_i S^x_i that is create([], N=N, h_x=-1.0, sign_check=False), whose rates −h_x/2 are non-negative.

dtype: object: Alias for field number 0

edge_i: Array: Alias for field number 1

edge_j: Array: Alias for field number 2

N: int: Alias for field number 3

J_z: object: Alias for field number 4

J_pm: object: Alias for field number 5

J_pmpm: object: Alias for field number 6

h_z: object: Alias for field number 7

h_x_half: object: Alias for field number 8

n_hx_static: int: Alias for field number 9

n_jpm_static: int: Alias for field number 10

has_jz: bool: Alias for field number 11

uniform_jz: bool: Alias for field number 12

has_hz: bool: Alias for field number 13

property edges: Reconstruct (n_edges, 2) from edge_i/edge_j (Python-side only).

property n_hx: Number of transverse-field off-diagonal elements (static Python int).

property n_jpm: Number of edge off-diagonal elements (static Python int).

property n_off: Total number of off-diagonal elements (static Python int).

property M_total

n_off + 1.

Type:: Total number of transitions including diagonal (static Python int)

classmethod create(edges, J_x=0, J_y=0, J_z=0, h_z=None, h_x=None, N=None, *, dtype=jnp.float32, sign_check=True)

Build a spin-1/2 XYZ Hamiltonian.

Parameters:

edges – (n_edges, 2) site index pairs (may be empty).
J_x – scalar or (n_edges,) exchange couplings. They enter only through J_pm = (J_x+J_y)/4 and J_pmpm = (J_x-J_y)/4; J_pmpm > 0 raises when sign_check is True (pass J_x and J_y the other way round yourself if that is a symmetry of your model – it is not when h_x != 0).
J_y – scalar or (n_edges,) exchange couplings. They enter only through J_pm = (J_x+J_y)/4 and J_pmpm = (J_x-J_y)/4; J_pmpm > 0 raises when sign_check is True (pass J_x and J_y the other way round yourself if that is a symmetry of your model – it is not when h_x != 0).
J_z – scalar or (n_edges,) S^z S^z coupling (any sign).
h_z – scalar, (N,) or None; longitudinal field, entering as −h_z Σ S^z_i (any sign).
h_x – scalar, (N,) or None; transverse field, entering as −h_x Σ S^x_i. Must be ≥ 0 when sign_check is True.
N – number of sites. Inferred from edges and any per-site array when None; give it explicitly for an operator without two-body terms (edges=[]).
dtype – JAX dtype of the couplings.
sign_check – enforce the sign-problem-free convention: raise on J_pm > 0, J_pmpm > 0 or h_x < 0. Pass False for measurement operators, which need J_pm >= 0 and are never propagated.

diagonal_weight(config)

Return just the diagonal transition weight (scalar).

Equivalent to get_rates(config)[0][0] but skips the (M_total,) rates allocation and the off-diagonal branches — useful for diagonal-operator measurement in hot loops.

get_rates(config, return_configs=False)

Return transition rates, and optionally the full connected-configs tensor.

The default skips the (M_total, N) connected-configs tensor to avoid an O(n_walkers * M_total * N) allocation inside vmapped evolution loops; use generate_state(config, idx) to materialize the single selected config. Pass return_configs=True to obtain the full tensor.

Returns:

rates: shape (M_total,) — rates[0] = diagonal, rates[1:] = off-diagonal.
new_configs: shape (M_total, N) when return_configs=True (new_configs[0] = config unchanged, new_configs[1:] = connected configs for each off-diagonal element), else None.

Return type:

(rates, new_configs)

generate_state(config, idx)

Generate new config by applying off-diagonal transition.

idx == 0: return config unchanged (diagonal) idx >= 1: apply the (idx-1)-th off-diagonal flip

property local_hilbert_space

Per-site values, ordered so that position i is basis digit i.

Matches this class’s own basis convention: digit = (1 - sigma)/2, so sigma = +1 is digit 0.

class pygfmc.hamiltonian.Spin

Spin-S Hamiltonian; see the module docstring for H and the rate layout.

Build instances with create().

dtype: object: Alias for field number 0

S: float: Alias for field number 1

N: int: Alias for field number 2

edge_i: Array: Alias for field number 3

edge_j: Array: Alias for field number 4

J_z: object: Alias for field number 5

J_pm: object: Alias for field number 6

J_pmpm: object: Alias for field number 7

h_z: object: Alias for field number 8

h_pm: object: Alias for field number 9

n_jpm_off: int: Alias for field number 10

n_jpmpm_off: int: Alias for field number 11

n_hpm_off: int: Alias for field number 12

has_jz: bool: Alias for field number 13

has_jpm: bool: Alias for field number 14

has_jpmpm: bool: Alias for field number 15

has_hz: bool: Alias for field number 16

has_hpm: bool: Alias for field number 17

uniform_jz: bool: Alias for field number 18

uniform_jpm: bool: Alias for field number 19

uniform_jpmpm: bool: Alias for field number 20

uniform_hz: bool: Alias for field number 21

uniform_hpm: bool: Alias for field number 22

property edges: (n_edges, 2) numpy array of site pairs.

property n_edges: Number of edges (static Python int).

property n_off: Number of off-diagonal transitions (static Python int).

property M_total

n_off + 1.

Type:: Total number of transitions including the diagonal

property d: Local Hilbert-space dimension = 2S + 1.

classmethod create(edges, S=0.5, J_z=0.0, J_pm=0.0, J_pmpm=0.0, h_z=None, h_pm=None, N=None, *, dtype=jnp.float32, sign_check=True)

Build a spin-S Hamiltonian.

Parameters:

edges – (n_edges, 2) site index pairs.
S – spin per site (half-integer ≥ 0); 2S must be a non-negative int.
J_z – scalar or (n_edges,) Sz Sz coupling. Any sign (diagonal).
J_pm – scalar or (n_edges,) S+S- + h.c. coupling. ≤ 0 (SPF).
J_pmpm – scalar or (n_edges,) S+S+ + h.c. coupling. ≤ 0 (SPF).
h_z – scalar, (N,), or None. Any sign (diagonal); enters H as -h_z * sum_i S^z_i.
h_pm – scalar, (N,), or None. ≥ 0 (SPF). ONE amplitude multiplying both S^+ and S^-, i.e. the transverse field -h_pm * sum_i (S^+_i + S^-_i) = -2 h_pm * sum_i S^x_i. A single shared amplitude keeps H Hermitian.
N – number of sites. Inferred from edges and any per-site array when None; give it explicitly for an operator without two-body terms (edges=[]).
dtype – JAX dtype for amplitudes.
sign_check – enforce SPF conventions if True.

diagonal_weight(config): Scalar H_{σσ}: J_z Σ Sz_i Sz_j + h_z Σ Sz_i.

get_rates(config, return_configs=False)

Return (rates, new_configs|None).

rates has shape (M_total,). See module docstring for layout.

generate_state(config, idx)

Return the configuration after applying transition idx.

Every transition changes at most two sites, so the target sites and their increments are selected as scalars and the configuration is updated once, instead of building one full copy per transition block.

move(config, idx): The sites changed by transition idx and their increments, (sites, deltas) of shape (2,) each (zero increments for the diagonal idx == 0); generate_state is config with deltas added at sites. Lets a guide update per-walker state after a move (see pygfmc.Guide).

property local_hilbert_space

Per-site values, ordered so that position i is basis digit i.

Matches this class’s own basis convention: digit = m + S, so m = -S is digit 0.

class pygfmc.hamiltonian.BoseHubbard

Bose-Hubbard model; see the module docstring for H and the rate layout.

Build instances with create().

dtype: object: Alias for field number 0

edge_i: Array: Alias for field number 1

edge_j: Array: Alias for field number 2

N: int: Alias for field number 3

n_max: int: Alias for field number 4

t: object: Alias for field number 5

eta: object: Alias for field number 6

U: object: Alias for field number 7

mu: object: Alias for field number 8

h_pm: object: Alias for field number 9

n_hop_off: int: Alias for field number 10

n_eta_off: int: Alias for field number 11

n_hpm_off: int: Alias for field number 12

has_t: bool: Alias for field number 13

has_eta: bool: Alias for field number 14

has_U: bool: Alias for field number 15

has_mu: bool: Alias for field number 16

has_h_pm: bool: Alias for field number 17

uniform_t: bool: Alias for field number 18

uniform_eta: bool: Alias for field number 19

uniform_U: bool: Alias for field number 20

uniform_mu: bool: Alias for field number 21

uniform_hpm: bool: Alias for field number 22

property edges: (n_edges, 2) numpy array of site pairs.

property n_edges: Number of edges (static Python int).

property n_off: Number of off-diagonal transitions (static Python int).

property M_total

n_off + 1.

Type:: Total number of transitions including the diagonal

classmethod create(edges, t=0.0, eta=0.0, U=0.0, mu=0.0, h_pm=None, N=None, n_max=4, *, dtype=jnp.float32, sign_check=True)

Build a Bose-Hubbard Hamiltonian.

Parameters:

edges – (n_edges, 2) array of (i, j) site index pairs.
t – scalar or (n_edges,) hopping amplitude. Must be ≥ 0 for sign-problem-free GFMC (the Hamiltonian carries the explicit minus sign in front of t).
eta – scalar or (n_edges,) nearest-neighbour pairing amplitude η, coupling the term -η (b†_i b†_j + b_i b_j) on each edge (same explicit minus-sign convention as t). Must be ≥ 0 for sign-problem-free GFMC.
U – scalar or (N,) onsite repulsion. Any sign (diagonal).
mu – scalar or (N,) chemical potential. Any sign (diagonal).
h_pm – scalar or (N,) amplitude multiplying BOTH b^dagger_i and b_i, i.e. -h_pm * sum_i (b^dagger_i + b_i) (minus sign as for the fields of SpinHalf and Spin; stored negated). None to disable. One shared amplitude keeps H Hermitian. Must be ≥ 0 for sign-problem-free GFMC.
N – number of sites. Inferred from edges and any per-site array when None; give it explicitly for an operator without two-body terms (edges=[]).
n_max – hard cutoff on per-site occupation (≥ 1).
dtype – JAX dtype for amplitudes.
sign_check – if True, enforce the sign-problem-free convention.

diagonal_weight(config): Scalar H_{σσ}: (U/2) Σ n(n-1) − μ Σ n.

get_rates(config, return_configs=False)

Return (rates, new_configs|None).

rates has shape (M_total,) with rates[0] the diagonal and the off-diagonal layout described in the module docstring. Rates for moves that would exceed n_max (or go negative) are masked to 0.

new_configs (when return_configs=True) has shape (M_total, N) with new_configs[0] == config.

generate_state(config, idx)

Return the configuration after applying transition idx.

idx == 0 returns config unchanged. Off-diagonal indices index into the layout documented at the module level.

Every transition changes at most two sites, so the target sites and their increments are selected as scalars and the configuration is updated once (the same scheme as Spin.generate_state()), instead of building one full candidate configuration per block.

move(config, idx): The sites changed by transition idx and their increments, (sites, deltas) of shape (2,) each (zero increments for the diagonal idx == 0); generate_state is config with deltas added at sites. Lets a guide update per-walker state after a move (see pygfmc.Guide).

property local_hilbert_space

Per-site values, ordered so that position i is basis digit i.

Matches this class’s own basis convention: the occupation is itself the digit.

class pygfmc.hamiltonian.ToricCode

Toric code Hamiltonian on an Lx × Ly torus.

H = −A_s Σ_v ∏_{i∈star(v)} σ^x_i
− B_p Σ_p ∏_{i∈plaq(p)} σ^z_i − h_z Σ_i σ^z_i − h_x Σ_i σ^x_i

Sites are the 2·Lx·Ly edges of the underlying square lattice:: horizontal edge h(x,y): index x + Lx·y (x∈[0,Lx), y∈[0,Ly)) vertical edge v(x,y): index Lx·Ly + x + Lx·y (x∈[0,Lx), y∈[0,Ly))

Spin values: config[i] ∈ {+1, −1} = σ^z eigenvalue.

Star (vertex) at (x,y) acts on edges: h(x,y), h(x−1,y), v(x,y), v(x,y−1). Plaquette (face) at (x,y) acts on edges: h(x,y), h(x,y+1), v(x,y), v(x+1,y).

Sign-problem-free convention: A_s ≥ 0, h_x ≥ 0 (off-diagonal rates −A_s, −h_x ≤ 0). No restriction on B_p or h_z (both are diagonal).

Choosing lam: it must exceed the largest diagonal element H_{σσ} visited. lam = |B_p|·Lx·Ly + |h_z|·N, the worst case, is always safe. Equilibrated walkers sit near the bottom of the spectrum, so a smaller value often suffices; check the diagonal statistics.

Couplings may be scalars or per-object arrays (see create()).

Layout of rates returned by get_rates:: rates[0] : diagonal H_{σσ} rates[1 : 1+n_stars_off] : star moves (−A_s each; active for all σ) rates[1+n_stars_off : M_total] : single-site X moves (−h_x each)

dtype: object: Alias for field number 0

stars: Array: Alias for field number 1

plaquettes: Array: Alias for field number 2

N: int: Alias for field number 3

Lx: int: Alias for field number 4

Ly: int: Alias for field number 5

A_s: object: Alias for field number 6

B_p: object: Alias for field number 7

h_z: object: Alias for field number 8

h_x: object: Alias for field number 9

n_stars: int: Alias for field number 10

n_plaq: int: Alias for field number 11

n_stars_off: int: Alias for field number 12

n_hx_off: int: Alias for field number 13

has_As: bool: Alias for field number 14

has_Bp: bool: Alias for field number 15

has_hz: bool: Alias for field number 16

has_hx: bool: Alias for field number 17

uniform_As: bool: Alias for field number 18

uniform_Bp: bool: Alias for field number 19

uniform_hz: bool: Alias for field number 20

uniform_hx: bool: Alias for field number 21

property n_off

Number of off-diagonal transitions (static Python int).

Required for a ToricCode instance to be usable as an off_diag_operator in run_gfmc (which gates on op.n_off > 0); the star flips and the single-site X moves are the off-diagonal ones.

property M_total: Total transitions including diagonal (static Python int).

classmethod create(Lx, Ly, A_s=1.0, B_p=1.0, h_x=0.0, h_z=0.0, *, dtype=jnp.float32, sign_check=True)

Build the toric code Hamiltonian on an Lx × Ly torus.

Parameters:

Lx – lattice dimensions
Ly – lattice dimensions
A_s – star (vertex) coefficient; scalar or (n_stars,) = (Lx*Ly,) per star. ≥ 0 for sign-problem-free GFMC.
B_p – plaquette coefficient; scalar or (n_plaq,) = (Lx*Ly,) per plaquette (diagonal, any sign).
h_x – σ^x field; scalar or (N,) per site. ≥ 0 for sign-problem-free GFMC.
h_z – σ^z field; scalar or (N,) per site (diagonal, any sign).
dtype – JAX dtype (default float32)
sign_check – raise ValueError if any A_s < 0 or h_x < 0 (default True)

Every coupling may be uniform (a scalar) or spatially varying (an array of the matching length), so disorder and boundary-modified stabilizers can be expressed without subclassing.

diagonal_weight(config): Scalar H_{σσ}: −B_p Σ_p ∏_{i∈p} σ^z_i − h_z Σ_i σ^z_i.

get_rates(config, return_configs=False)

Rates (M_total,) and optionally the (M_total, N) connected-configs tensor.

Off-diagonal rates are (−A_s, −h_x) independent of config, because σ^x is always active (no flippability constraint). With per-star / per-site couplings they are config-independent but differ from slot to slot.

generate_state(config, idx)

Apply off-diagonal transition idx to config.

idx == 0 : return config unchanged (diagonal) 1 <= idx <= n_stars_off : flip the 4 sites of star (idx−1) n_stars_off < idx : flip site (idx − 1 − n_stars_off)

property local_hilbert_space

Per-site values, ordered so that position i is basis digit i.

Matches this class’s own basis convention: digit = (1 - sigma)/2, so sigma = +1 is digit 0.

class pygfmc.hamiltonian.PlaquetteModel

Plaquette ring-exchange / Rokhsar–Kivelson Hamiltonian.

H = −g · Σ_p (flip term on plaquette p) + mu · Σ_p |flippable_p⟩⟨flippable_p|

Each plaquette contributes

one off-diagonal element (the flip S⁺S⁻S⁺S⁻ … around the loop), active only when the plaquette is flippable — i.e. its spins strictly alternate as (+,−,+,−,…) or (−,+,−,+,…) in the order supplied to create().

one diagonal element (the RK counting potential μ_p) added when the plaquette is flippable.

Layout of the rates returned by get_rates:: rates[0] : diagonal, Σ_p μ_p · 1[flippable_p] rates[1 : N_p + 1] : one entry per plaquette, −g_p if flippable else 0

The site count N is inferred from the plaquettes (max index + 1); configurations are ±1 floats. Plaquettes must have EVEN length (create() raises otherwise): an odd loop cannot alternate all the way round, so its “flip” would change Sz_tot and is no ring exchange.

Fields:

plaquettes int array (N_plaquette, plaquette_size) — site indices: per plaquette, listed consecutively around the loop.
g flip coupling, shape (1,) uniform or (N_plaquette,).: Must be ≥ 0 (sign-problem-free convention): the off-diagonal matrix element is −g ≤ 0. The check can be disabled with sign_check=False when constructing measurement operators — note that an off-diagonal measurement operator needs rates ≥ 0, i.e. g ≤ 0, the opposite restriction.
mu RK counting coupling (diagonal, any sign), shape: (1,) uniform or (N_plaquette,).

dtype: object: Alias for field number 0

plaquettes: Array: Alias for field number 1

N: int: Alias for field number 2

g: Array: Alias for field number 3

mu: Array: Alias for field number 4

signs: Array: Alias for field number 5

n_plaquette: int: Alias for field number 6

plaquette_size: int: Alias for field number 7

has_g: bool: Alias for field number 8

has_mu: bool: Alias for field number 9

uniform_g: bool: Alias for field number 10

uniform_mu: bool: Alias for field number 11

property M_total

N_plaquette + 1.

Type:: Total transitions incl. diagonal (static Python int)

property n_off: Number of off-diagonal elements (static Python int).

classmethod create(plaquettes, g=None, mu=None, *, dtype=jnp.float32, sign_check=True)

Build a plaquette ring-exchange / RK Hamiltonian.

Parameters:

plaquettes – (N_plaquette, plaquette_size) site indices, listed consecutively around each loop. N = max index + 1, and every site from 0 to N-1 must appear in at least one plaquette (a gap raises ValueError). The plaquette length must be EVEN (odd raises, whatever sign_check).
g – scalar, (N_plaquette,) or None; flip amplitude in -g Σ flip. Must be ≥ 0 when sign_check.
mu – scalar, (N_plaquette,) or None; RK potential on flippable plaquettes (any sign).
dtype – JAX dtype of the couplings.
sign_check – enforce g ≥ 0. Pass False for an off-diagonal measurement operator, which instead needs g ≤ 0.

property local_hilbert_space

Per-site values, ordered so that position i is basis digit i.

Matches this class’s own basis convention: spin-1/2 degrees of freedom encoded as +-1.

diagonal_weight(config): Scalar diagonal: Σ_p μ_p · 1[flippable_p]. Matches get_rates(c)[0][0].

get_rates(config, return_configs=False): Return rates (M_total,) and optionally the (M_total, N) connected-configs tensor.

generate_state(config, idx)

Apply off-diagonal transition.

idx == 0 : return config unchanged (diagonal). idx >= 1 : flip all sites of plaquette (idx - 1).

Measurement operators

pygfmc.hamiltonian.StructureFactor(positions, k, channel='SpSm', dtype=jnp.float32)

The k-resolved correlator AB(0) + AB(k) as ONE SpinHalf.

O = sum_{i<j} [1 + cos(k.(r_i - r_j))] (A_iB_j + h.c.) O = AB(0) + AB(k) AB(q) = sum_{i<j} cos(q.(r_i - r_j)) (A_iB_j + h.c.)

Returns a SpinHalf with M_total = 1 + n_kept, where n_kept is the number of pairs i < j with nonzero weight, for use as an off-diagonal measurement operator (off_diag_operators of run_gfmc). Combined with the separately measured k = 0 operator it gives the second-moment correlation length xi^2 = [S(0)/S(k) - 1] / khat^2.

Positions, as in SzSzStructureFactor. The displacement of a pair is r_i - r_j of the positions as given; no periodic wrap is applied and no image is chosen. For a k in the reciprocal lattice of the simulation cell (k.L = 2*pi*integer for every periodic vector L of the system, see is_valid_k()) every periodic image gives the same cos(k.dr), so the result does not depend on which image the positions encode. k is not checked here.

``channel=”SzSz”`` is diagonal and warns. A = B = S^z gives A_iB_j + h.c. = 2 S^z_iS^z_j, so the operator is sum_{i<j} 2[1 + cos(k.dr)] S^z_iS^z_j, which is exactly SzSzStructureFactor(positions, [0, k]) summed over its two wavevectors, i.e. the same AB(0) + AB(k) structure as the other channels. The result has M_total = 1 and goes in diag_operators. A UserWarning points out that SzSzStructureFactor measures the same correlator for many k at once in O(N n_k) and needs neither the 1 + cos shift nor the k = 0 subtraction, since a diagonal operator has no sign constraint.

Parameters:

positions – (N, d) site coordinates; row i is the position of the site at index i in config, so N must match the Hamiltonian’s N. d is arbitrary; it only has to match k.
k – (d,) wavevector, in the same units as positions. One operator measures one wavevector.
channel – one of the four below; case- and space-insensitive.
dtype – JAX dtype of the couplings.

channel	per-pair operator
`"SpSm"`	`S+S- + S-S+` = `2(S^xS^x + S^yS^y)` (transverse)
`"SpSp"`	`S+S+ + S-S-`
`"SpSm+SpSp"`	`4 S^xS^x`
`"SzSz"`	`2 S^zS^z` (diagonal, warns)

Why these off-diagonal three. get_rates emits one off-diagonal slot per edge, where(anti_parallel, J_pm, J_pmpm), so per edge only PM = S+S- + S-S+ and PP = S+S+ + S-S- exist, and an off-diagonal measurement needs every rate >= 0. The reachable set is therefore the non-negative cone a*PM + b*PP (entered via J_x = 2(a+b)W, J_y = 2(a-b)W). Since PM + PP = 4 S^xS^x and PM - PP = 4 S^yS^y, YY needs a negative coefficient: build it as the "SpSm" operator minus the "SpSp" one. S.S = PM/2 + ZZ mixes the diagonal and off-diagonal paths: measure "SpSm" (off-diagonal) and "SzSz" (diagonal) separately, then add. Mixed species (XY, XZ, …) are not in the span.

The weight is ``1 + cos``, not ``cos``. Off-diagonal measurement samples a successor by inverse CDF over the rates, so mixed-sign weights are impossible; the +1 shift makes every weight non-negative at the price of the additive AB(0), which the caller removes by subtracting the separately measured k = 0 operator. Pairs whose weight vanishes (cos(k.dr) = -1, to within 1e-12, since an image other than the minimum one leaves a rounding residue) are dropped, so M_total does not reflect the pair count.

Memory. The operator has one transition slot per pair, about N^2/2, and applying it to the ensemble materialises (n_walkers, M_total) rates (times N for a guide that needs the successor configurations). Pass off_diag_batch_size to run_gfmc to apply it in walker chunks; the result is bit-identical.

Sign conditions. create’s J_pm <= 0 / J_pmpm <= 0 are propagator (sign-problem) conditions, whereas a measurement operator needs rates >= 0, i.e. J_pm >= 0 on every channel and J_pmpm > 0 on the two with b > 0. Every channel therefore builds the operator with sign_check=False.

Normalisation. For "SpSm" and "SzSz", AB(q) above equals the ordered sum Re sum_{i!=j} exp(-i q.r_ij) A_iB_j, with no factor 2, because the adjoint of A_iB_j is its i<->j partner (S+_jS-_i, resp. S^z_jS^z_i): listing each pair once and adding +h.c. reproduces the ordered sum. For "SzSz", AB(q) is therefore the S^{zz}(q) of SzSzStructureFactor. The identity with the ordered sum fails for "SpSp", whose ordered sum is not Hermitian; the definition of AB(q) above holds for every channel.

class pygfmc.hamiltonian.SzSzStructureFactor(positions, k_vectors, dtype=jnp.float32)

$S^{zz}(\vec k)$ for one or many wavevectors, from positions.

The observable is the pair sum

\[S^{zz}(\vec k) = \sum_{i<j} 2\cos(\vec k\cdot\vec{dr}_{ij}) S^z_i S^z_j ,\]

evaluated through the equivalent single sum

\[\rho_{\vec k} = \sum_i e^{i\vec k\cdot\vec r_i}\sigma^z_i , \qquad |\rho_{\vec k}|^2 = \sum_{ij}\cos(\vec k\cdot(\vec r_i-\vec r_j)) \sigma^z_i\sigma^z_j ,\]

so that (with $\sigma^z=\pm1$, $S^z=\sigma^z/2$) diagonal_weight returns $(|\rho_{\vec k}|^2 - N)/4$. The i == j terms are always dropped: they are the constant $N/4$ (sigma^2 = 1), carry no k dependence, and would shift every wavevector by the same amount.

Why the single sum. It is an exact algebraic rewriting, holding configuration by configuration, and costs O(N n_k) per walker instead of the pair sum’s O(N^2); measuring many k in one operator is therefore nearly free compared with a run that measures nothing.

It is exact only for ``k`` in the reciprocal lattice of the simulation cell, i.e. k . L = 2 pi x integer for every supercell vector L, which is exactly the set of wavevectors a finite periodic cell has. There every periodic image of a pair carries the same phase, so the result agrees with a pair sum built from minimum-image or absolute displacements. Off that grid the minimum-image pair sum is not a structure factor (it can even be negative, whereas $|\rho_k|^2 \geq 0$).

Supplying valid wavevectors is the caller’s job; nothing here checks it (use is_valid_k()). The positions alone cannot establish validity, because folded coordinates do not record the boundary identification (x = 0..7 is equally an L=8 ring, a half-empty L=16 cell or an open chain). Generate the grid rather than typing it:

cell = lat.extras['supercell_dims'][:, None] * lat.extras['primitive_vectors']
ks = 2 * np.pi * np.asarray(ms) @ np.linalg.inv(cell).T   # ms: integer rows

(pyrochlore ships both positions and that cell in extras.)

``diagonal_weight`` returns one value per wavevector, which the GFMC pipeline carries as a trailing axis: raw['diag'][name] is (n_steps, n_k, n_p), obs[name][p] is an (n_k,) array and the HDF5 dataset run_<n>/observables/<name> is (n_p, n_k). The k axis carries no labels – store op.k_vectors yourself, e.g. as stats['parameters/k_vectors'] before write_to_HDF5 (it is not a MUST_MATCH key, so it cannot break a later append). stats['autocorrelation/diag/<name>'] is resolved per wavevector, shape (n_k, n_p). The operator may share a run with diagonal operators of any other width; memory of its ancestry buffer grows linearly in n_k.

Parameters:

positions – Site coordinates, shape (N, d). Row i is the position of the site at index i in config.
k_vectors – One wavevector (d,) or several (n_k, d).
dtype – Dtype of the stored phase tables (default jnp.float32). They are always built in float64, since cos/sin of a large k.r is where the precision is lost, not the contraction.

Example

>>> ks = 2 * np.pi * np.arange(L) / L                  # allowed grid
>>> sf = SzSzStructureFactor(positions, ks[:, None])
>>> *_, obs, stats, raw = run_gfmc(..., diag_operators={"Szz": sf})
>>> obs["Szz"][p]                                      # (n_k,)

diagonal_weight(config): (n_k,) for a single (N,) config, (n, n_k) for a batch.

pygfmc.hamiltonian.is_valid_k(k_vectors, lattice_vectors, atol=1e-8)

Whether each wavevector lies in the reciprocal lattice of the system.

A wavevector is valid when k . L_j = 2 pi x integer for every periodic vector L_j of the system. Exactly then every periodic image of a pair carries the same phase: shifting a site to a periodic image, r -> r + sum_j m_j L_j with integers m_j, changes k . (r_i - r_j) by a multiple of 2 pi and leaves cos(k . (r_i - r_j)) unchanged. The direct displacements r_i - r_j of the positions as given therefore suffice, without choosing periodic images, so SzSzStructureFactor and StructureFactor() do not depend on which image the positions encode.

Parameters:

k_vectors – One wavevector (d,) or several (n_k, d).
lattice_vectors – (d, d), rows the periodic vectors of the whole system, e.g. supercell_dims[:, None] * primitive_vectors or kagome’s extras['torus']. They are the vectors the boundary identifies, not the primitive vectors of the unit cell.
atol – Absolute tolerance on k . L_j / (2 pi) being an integer.

Returns:

bool for a single (d,) wavevector, a boolean (n_k,) array otherwise.

Example

>>> L = np.array([[8.0]])                        # L = 8 ring
>>> is_valid_k(2 * np.pi * np.arange(3)[:, None] / 8, L)
array([ True,  True,  True])
>>> is_valid_k([1.0], L)
False

Lattices

class pygfmc.hamiltonian.Lattice(name, N, edges, J_x=None, J_y=None, J_z=None, extras=<factory>)

A finite lattice: a site count, an edge list and optional metadata.

name

lattice family, e.g. "chain", "square", "pyrochlore".

Type:: str

N

number of sites.

Type:: int

edges

list of (i, j) site pairs.

Type:: list[tuple[int, int]]

J_x, J_y, J_z: optional suggested couplings (scalar or per-edge list); only chain() sets them, and nothing in the package reads them automatically.

extras

factory-specific data. The pyrochlore factories store hexagons, type_per_hexagon, positions, supercell_dims, primitive_vectors and sublattice_offsets; kagome() stores hexagons, hexagon_centers, positions, supercell_dims, primitive_vectors, sublattice_offsets and torus.

Type:: dict

pygfmc.hamiltonian.chain(n_sites, periodic=True)

Nearest-neighbour chain of n_sites sites.

The closing bond (n_sites-1, 0) is added when periodic and n_sites > 2. Sets the suggested couplings J_x = J_y = -1, J_z = 1.

pygfmc.hamiltonian.pyrochlore(nx, ny, nz, periodic=True)

Pyrochlore lattice built from explicit FCC unit-cell edge templates (O(N)).

4 sites per FCC unit cell. N = 4*nx*ny*nz, n_edges = N*3 (periodic). Cell index: z + y*nz + x*ny*nz (z changes fastest). Periodic boundaries require nx, ny, nz >= 2. Positions use a conventional cubic cell of side 1 (NN distance 1/(2√2)); extras carries hexagons, positions and the supercell (see Lattice). With open boundaries only the hexagons lying completely inside the lattice are listed.

12 edge templates per cell:

Within A-tetrahedron (offset 0,0,0):
    (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
Inter-cell / B-tetrahedron connections:
    src=1 → tar=0, offset(+x, 0,  0)
    src=2 → tar=0, offset( 0,+y,  0)
    src=3 → tar=0, offset( 0, 0, +z)
    src=2 → tar=1, offset(-x,+y,  0)
    src=3 → tar=1, offset(-x, 0, +z)
    src=3 → tar=2, offset( 0,-y, +z)

pygfmc.hamiltonian.kagome(Lx, Ly, periodic=True)

Kagome lattice of Lx x Ly unit cells (O(N)).

A triangular Bravais lattice with primitive vectors a1 = (2, 0), a2 = (1, sqrt(3)) and three sites per cell at 0, a1/2 and a2/2 (nearest-neighbour distance 1), so N = 3*Lx*Ly and n_edges = 2N (periodic). Site index: 3*(x + Lx*y) + b for cell (x, y) and sublattice b. Periodic boundaries require Lx, Ly >= 2.

extras carries positions (N, 2), primitive_vectors, supercell_dims, sublattice_offsets, the torus vectors torus (rows Lx*a1 and Ly*a2), and one hexagon per cell: hexagons (n_hex, 6), the site indices in counter-clockwise order starting at angle 0, and hexagon_centers (n_hex, 2). With open boundaries only the hexagons fully inside the lattice are listed.

Exact diagonalization

pygfmc.hamiltonian.exact_diagonalization(hamiltonian, local_hilbert_space=None, N=None, diag_operators=None, off_diag_operators=None, number_of_states=1, return_states=False)

Lowest eigenpairs of hamiltonian, with observables on each.

Parameters:

hamiltonian – Any object implementing the Hamiltonian protocol (M_total, get_rates, generate_state). Its matrix is rebuilt from get_rates alone.
local_hilbert_space – The per-site values, in basis-digit order, e.g. (1, -1) for spin-1/2, (-1, 0, 1) for S = 1, (0, 1, 2) for bosons with n_max = 2. The order fixes the basis indexing of the returned states (energies do not depend on it). Defaults to hamiltonian.local_hilbert_space, which the pre-implemented classes provide.
N – Number of sites; the basis has len(local_hilbert_space) ** N states. Defaults to hamiltonian.N, so it only has to be given for a class that exposes no N – a user-defined Hamiltonian implementing just the protocol.
diag_operators – dict name → operator, measured as Σ_x |ψ_x|² O(x) via diagonal_weight.
off_diag_operators – dict name → operator, measured as ψᵀ O ψ with O built from its own get_rates.
number_of_states – How many lowest eigenpairs to solve for.
return_states – Also return the eigenvectors, shape (dim, number_of_states). Needed by ExactGroundState, which turns them into guide amplitudes.

Returns:

(energies, diag_obs, off_diag_obs), or that tuple with states appended when return_states is True. energies has shape (number_of_states,); the two dicts mirror the input dicts key for key, each value carrying one entry per state along the leading axis.

Note

The routine works in float64 throughout, whatever precision the caller is running in. The one residual it cannot remove is the Hamiltonian’s own dtype: a dtype=jnp.float32 model stores float32 couplings, and the spectrum inherits that ~1e-7 relative error; build the model with dtype=jnp.float64 when it matters.

Raises:: ValueError – on a basis larger than 2**24 states, an N or local_hilbert_space contradicting what the Hamiltonian reports (N, d or the local values), a local_hilbert_space with repeated or non-finite values, number_of_states outside [1, dim], a transition of nonzero amplitude that leaves the basis, a non-Hermitian or non-finite operator matrix, or an estimated allocation above _MAX_BYTES.

`pygfmc.wavefunction`

Guiding wavefunctions and their training. train_with_netket needs the optional extra netket: pip install "pygfmc[netket]", or from source pip install -e "./codes[netket]" in the root of the cloned repository (see Package and Installation).

`JastrowSpinHalf`(N, pairs, beta[, h, dtype])	Jastrow guiding wavefunction with sparse two-body couplings.
`ClusterCPS`(groups[, amplitudes, init, seed, ...])	Cluster-CPS guide composed from one or more patch groups.
`ClusterMLP`(groups[, params, seed, dtype, N])	Per-cluster MLP guide for SpinHalf, with one shared tanh MLP per cluster group.
`ExactGroundState`(hamiltonian, *[, log_amp, ...])	Guiding wavefunction with amplitudes from the exact ground state.
`JastrowBoseHubbard`(N, pairs, beta[, h, d, ...])	Jastrow guide for boson occupations.
`JastrowSpin`(N, pairs, beta[, h, d, dtype, ...])	Jastrow guide for spin-S configurations.
`train_with_netket`(hamiltonian, guide[, ...])	Train `guide` on `hamiltonian` with NetKet VMC and return a new guide.
`make_hilbert`(local_hilbert_space, N[, total])	A NetKet Hilbert space whose local states are exactly these values.

class pygfmc.wavefunction.JastrowSpinHalf(N, pairs, beta, h=None, dtype=jnp.float32)

Jastrow guiding wavefunction with sparse two-body couplings.

ψ_T(x) = exp(Σ_{(i,j) ∈ pairs} J_ij x_i x_j + Σ_i h_i x_i)

The constructor is Hamiltonian-agnostic: pass N (number of sites), a (n_pairs, 2) pair list and the couplings beta. No reference to a SpinHalf operator is kept on the guide. The coupling follows the convention log ψ = -Σ_<ij> β_ij x_i x_j + Σ_i h_i x_i and is stored internally as J = -beta (exposed as the attribute beta).

Parameters:

N – int, total number of sites in the system.
pairs – (n_pairs, 2) int array of site index pairs (i, j). Order within a pair does not matter (canonicalised to (min, max) internally).
beta – required; a scalar (one tied coupling for every pair) or an array of length n_pairs (one coupling per pair).
h – one-body fields — a scalar (one tied field on every site), an (N,) float array (one per site), or None (no field).
dtype – JAX dtype for amplitudes.

The guide keeps the parameters in the form they were given, and train_with_netket() trains exactly those: a scalar beta stays a single coupling, h=None stays without a field, a scalar h stays a single field.

Operators: log_diff_from_operator() and build_cache() index their output by the SpinHalf transition layout, so they accept a SpinHalf (Hamiltonian or measurement operator, e.g. StructureFactor()) or any operator exposing n_hx_static, n_jpm_static, edge_i and edge_j; any other operator (ToricCode, PlaquetteModel, …) raises TypeError. log_psi() and log_diff_from_configs() take configurations of any spin-1/2 operator.

property site_J: Couplings J = -beta per site and partner slot, (N, max_deg).

property h: One-body field per site, (N,) (zeros when constructed with h=None).

property tied_beta: True when beta was a single coupling shared by all pairs.

property beta

Per-pair coupling in the paper’s convention, beta = -J.

Read from the stored couplings (a tied scalar is broadcast to every pair), so it follows NetKet training or a tree_map and survives a pytree round trip / jit. The order is the order of pairs given to the constructor.

save(path, group='jastrow_spin_half', **attrs)

Write the Jastrow couplings to an HDF5 file.

Stores the pair list (in constructor order), per-pair J, the one-body field h, and N, so load() rebuilds the same guide with beta in the same order. Any extra keyword arguments are written as scalar HDF5 attributes on the group, e.g. E_vmc=..., E_vmc_err=....

Parameters:

path – HDF5 file path (created or overwritten).
group – HDF5 group name to write under.
**attrs – extra metadata stored as group attributes.

classmethod load(path, hamiltonian=None, group='jastrow_spin_half', dtype=jnp.float32)

Load a Jastrow saved with save().

Parameters:

path – HDF5 file path written by save().
hamiltonian – optional SpinHalf. If provided, N is taken as the largest of saved_N, hamiltonian.N and the highest edge index + 1, and h is zero-padded when the new system is larger; otherwise N=saved_N.
group – HDF5 group name to read from.
dtype – JAX dtype for the rebuilt guide.

classmethod from_jax(dtype, N, site_other, site_J, site_valid_mask, h, pair_site, pair_slot, site_pair)

Rebuild a guide from its internal arrays (used by pytree unflatten).

Bypasses __init__. pair_site/pair_slot locate each pair’s coupling in site_J (they define beta).

log_psi(config, cache=None)

log ψ_T(config) = Σ_pairs J_ij x_i x_j + Σ_i h_i x_i.

Accepts a single (N,) configuration or a (n, N) batch, matching the guide protocol’s third routine. cache is accepted for protocol parity and unused: the sum is already O(N·max_deg) and there is nothing operator-specific to cache.

Evaluated through site_other/site_J, which store each pair from both endpoints – hence the 1/2 on the pair term.

GFMC itself never calls this; it uses the ratio routines (log_diff_from_operator/log_diff_from_configs). It exists for diagnostics and for testing those ratios against an explicit amplitude.

log_diff_from_configs(config, new_configs, cache=None)

log ψ(x′) − log ψ(x) for each row of new_configs.

Returns shape (new_configs.shape[0],). Generic: works for any new_configs regardless of how they were produced.

cache is accepted for protocol parity and unused.

log_diff_from_operator(config, op, cache=None)

Return log-amplitude ratios indexed by a SpinHalf operator’s layout.

Layout matches op.get_rates of a SpinHalf: [diag (=0), op.n_hx site flips at sites 0..N-1, op.n_jpm op-edge flips]. No successor configurations are materialised.

cache=None (default) rebuilds the per-edge external-partner mask on the fly via _edge_log_diff_for_edges. Pass a JastrowOpCache (obtained from self.build_cache(op)) to use the cached fast path instead.

op must expose the SpinHalf layout (n_hx_static, n_jpm_static, edge_i, edge_j); any other operator raises TypeError.

build_cache(op)

Build and return the external-partner cache for op’s edges.

Returns a JastrowOpCache pytree of shape-(n_edges_op, 2, max_deg) (edge_other, edge_J) arrays — the result of joining self.site_other / self.site_J with op.edge_i / op.edge_j and masking out the i↔j bond endpoints themselves. Pass the returned cache to log_diff_from_operator(..., cache=...).

Results are memoised by the content of op’s edge layout and of the couplings J (which edge_J bakes in), so repeated calls with the same operator (e.g. once from run_gfmc and again from statistics()) skip the rebuild, while a guide whose couplings changed (e.g. a trained copy) gets its own tables.

op must expose the SpinHalf layout (n_hx_static, n_jpm_static, edge_i, edge_j); any other operator raises TypeError.

class pygfmc.wavefunction.ClusterCPS(groups, amplitudes=None, init='zero', seed=None, dtype=jnp.float32)

Cluster-CPS guide composed from one or more patch groups.

Parameters:

groups – list of patch-group specs (see the module docstring), or a single (n_patches, m) cluster array for a one-group guide.
amplitudes – optional LIST of amplitude tables, one per group. Each table has shape (n_types, M_g) (or (M_g,) if the group has one table) after that group’s symmetry and parity reduction. If omitted, amplitudes are initialised according to init. When groups is a bare cluster array a single table may be passed without the list.
init – 'zero' (default; the uniform state log psi == 0) or 'gaussian' (std 0.01) when amplitudes is omitted.
seed – RNG seed for Gaussian initialisation; None draws one, kept as guide.seed (None on a guide rebuilt from its pytree, e.g. after training).
dtype – JAX dtype for amplitude arrays.

Attributes: groups (tuple of ClusterCPSGroup), n_groups, N (largest site index + 1), n_params.

Operators: log_diff_from_operator() and build_cache() index their output by the SpinHalf transition layout, so they accept a SpinHalf (Hamiltonian or measurement operator) or any operator exposing n_hx_static, n_jpm_static, edge_i and edge_j; any other operator (ToricCode, PlaquetteModel, …) raises TypeError. log_psi() and log_diff_from_configs() take configurations of any spin-1/2 operator.

property n_params: Total number of amplitude-table entries over all groups.

log_psi(configs, cache=None)

Evaluate log psi for one config (N,) or a batch (..., N).

cache is accepted for protocol parity and ignored.

log_diff_from_configs(config, new_configs, cache=None): log psi(sigma') - log psi(sigma) for each row of new_configs.

build_cache(op)

Build affected-patch caches for SpinHalf-like site and edge flips.

op must expose the SpinHalf layout (n_hx_static, n_jpm_static, edge_i, edge_j and optionally N); any other operator raises TypeError. Returns a ClusterCPSCache. The cache bakes in the current amplitudes, so rebuild it after changing them. Not memoised.

log_diff_from_operator(config, op, cache=None)

Log-amplitude ratios in the SpinHalf transition layout.

Layout matches op.get_rates: [diag (=0), n_hx site flips, n_jpm edge flips]. With a ClusterCPSCache from build_cache() each affected patch index is updated by index arithmetic (idx_new = idx_old ± 2^slot) without materialising any successor configuration; with cache=None the successors are built via op.get_rates(config, return_configs=True). op must expose the SpinHalf layout either way; any other operator raises TypeError.

save(path, group='cps', **attrs)

Write the guide to path (overwritten) under group.

Each cluster group goes to a subgroup group_<i> (clusters, amplitudes, cluster_type, orbit tables, optional automorphisms); n_groups, N, n_params and any extra **attrs are written as attributes of group.

classmethod load(path, group='cps', dtype=jnp.float32): Load a guide written by save() from HDF5 group group.

class pygfmc.wavefunction.ClusterMLP(groups, params=None, seed=None, dtype=jnp.float32, N=None)

Per-cluster MLP guide for SpinHalf, with one shared tanh MLP per cluster group.

The log-amplitude is summed over groups g and their clusters c:

log ψ(x) = Σ_g Σ_{c ∈ clusters_g} MLP_g( x[clusters[g][c]] )

Parameters:

groups – list of [clusters, hidden] specs, one per group: clusters is an (n_clusters_g, k_g) int array (all clusters of a group share the size k_g) and hidden a list of hidden-layer widths ([] or None for a linear model; a bare int is rejected).
params – optional list, one entry per group, each a list of (W_i, b_i) pairs for layers i = 0 … n_layers-1. The last layer is the output (no activation). If omitted, weights are drawn from a Gaussian of std 0.01 and the biases are zero; if given, the layer shapes must match groups.
seed – RNG seed for the initialisation; None draws one, kept as guide.seed (None on a guide rebuilt from its pytree, e.g. after training).
dtype – JAX dtype.
N – number of sites; defaults to the largest site index + 1 (pass it if trailing sites belong to no cluster).

Operators: log_diff_from_operator() and build_cache() index their output by the SpinHalf transition layout, so they accept a SpinHalf (Hamiltonian or measurement operator) or any operator exposing n_hx_static, n_jpm_static, edge_i and edge_j; any other operator (ToricCode, PlaquetteModel, …) raises TypeError. log_psi() and log_diff_from_configs() take configurations of any spin-1/2 operator.

property n_params: Total number of weights and biases over all groups.

log_psi(config, cache=None)

log ψ(x) = Σ_g Σ_{c ∈ clusters_g} MLP_g(x[clusters[g][c]]).

config may be a single (N,) config or a batch (n_configs, N); the return is a scalar or (n_configs,) accordingly (same convention as ClusterCPS.log_psi()). cache is accepted for protocol parity but unused — the cluster cache only accelerates log_diff_from_operator().

log_diff_from_configs(config, new_configs, cache=None)

log ψ(x') − log ψ(x) for each row of new_configs.

Brute-force path: runs full log_psi per row. Memory: O(M·N) for new_configs. Use log_diff_from_operator(..., cache=cache) for the memory- and compute-light path when M can be huge.

log_diff_from_operator(config, op, cache=None)

Log-amplitude ratios indexed by op’s transition layout.

Layout matches op.get_rates of a SpinHalf:: [diag (=0), n_hx site flips at sites 0..N-1, n_jpm op-edge flips]

With cache=None: falls back to log_diff_from_configs(config, op.get_rates(config, return_configs=True)[1]). With a ClusterMLPCache from self.build_cache(op): uses the affected-cluster lookup to recompute MLPs only for clusters whose x_c actually changed under each transition. op must expose the SpinHalf layout (n_hx_static, n_jpm_static, edge_i, edge_j) either way; any other operator raises TypeError.

build_cache(op)

Build a ClusterMLPCache for fast log_diff_from_operator(op).

For each group g, computes:

site_to_clusters[g] of shape (N, max_aff_g) — for each site, the cluster indices of group g containing that site (padded with -1 to a uniform max length).
pair_to_clusters[g] of shape (n_jpm, max_aff_pair_g) — for each op-edge (edge_i[e], edge_j[e]), the union of cluster indices containing either endpoint.

Memoised by the byte content of op’s edge layout. op must expose the SpinHalf layout (n_hx_static, n_jpm_static, edge_i, edge_j); any other operator raises TypeError.

save(path, group='cluster_mlp', **attrs): Write clusters + per-group variable-depth MLP params to HDF5.

classmethod load(path, group='cluster_mlp', dtype=jnp.float32): Load a guide written by save().

class pygfmc.wavefunction.ExactGroundState(hamiltonian, *, log_amp=None, dtype=jnp.float32, local_hilbert_space=None, N=None)

Guiding wavefunction with amplitudes from the exact ground state.

Parameters:

hamiltonian – any Hamiltonian implementing the protocol (M_total, get_rates, generate_state), e.g. a SpinHalf / Spin / BoseHubbard / ToricCode instance.
log_amp – optional precomputed log |ψ_0|, shape (d^N,). If None, ψ_0 is computed with pygfmc.hamiltonian.exact.exact_diagonalization() and the ground-state energy is stored as E0 (NaN otherwise). The diagonalization returns ONE lowest eigenvector. When the ground state is degenerate – e.g. across sectors the Hamiltonian does not connect, such as the flippability sectors of PlaquetteModel or the B_p sectors of ToricCode – that vector is an arbitrary element of the degenerate subspace and can vanish on part of the basis (an entire sector), where |ψ_0| is floored at 1e-30 before the log. For a run inside one sector, supply that sector’s ground state through log_amp (e.g. from exact_diagonalization(..., number_of_states=k, return_states=True) restricted or combined as needed). No degeneracy check is made.
dtype – JAX dtype for the stored log_amp.
local_hilbert_space – Per-site values, e.g. (-1, 1) for spin-1/2. Defaults to hamiltonian.local_hilbert_space.
N – Number of sites. Defaults to hamiltonian.N.

basis_index(config): Index of config in the amplitude vector.

log_diff_from_configs(config, new_configs=None, cache=None)

Log amplitude ratios for every transition, shape (M_op,).

When new_configs is given (rows = successor configs of some operator, as produced by op.get_rates(config, return_configs=True)), the ratios follow that operator’s transition layout — this is the path GFMC takes (pygfmc.engine.evolution.transition_rates always passes the successors for this guide). When new_configs is None the successors are generated from the bound Hamiltonian: x_k_new = hamiltonian.generate_state(x, k), and element 0 is pinned to exactly 0.0 (diagonal, x unchanged).

For Hamiltonians with bounded local Hilbert space (e.g. Spin S = 1, BoseHubbard) successor configs may be out of range for transitions whose rate is zero (e.g. S+ acting on m = +S). Those produce out-of-range basis indices, which we clip into the valid range so the lookup is well-defined. The corresponding rate is 0, so the looked-up amplitude is multiplied by zero and never affects the dynamics. cache is accepted for protocol parity and unused.

save(path, group='exact_ground_state', **attrs)

Write log_amp plus metadata to an HDF5 file.

Datasets: log_amp and local_hilbert_space. Group attributes: N, family (e.g. "SpinHalf"), E0, plus any user keyword arguments. load() checks family and N against the Hamiltonian it is rebound to.

classmethod load(path, hamiltonian, group='exact_ground_state', dtype=jnp.float32)

Load a guide saved with save() and bind it to hamiltonian.

Unlike the other guides’ load, hamiltonian is required: the guide keeps a reference to it (for M_total and generate_state). Raises ValueError if the saved family or N does not match the Hamiltonian.

class pygfmc.wavefunction.JastrowBoseHubbard(N, pairs, beta, h=None, d=None, dtype=jnp.float32, incremental=False)

Jastrow guide for boson occupations.

\[\log\psi_T(x) = -\sum_{(i,j)} \beta_{ij} x_i x_j + \sum_i h_i x_i - \sum_i d_i x_i^2\]

Parameters:

N – number of sites.
pairs – (n_pairs, 2) site pairs of the two-body term, e.g. the Hamiltonian edges. Order within a pair does not matter.
beta – a scalar (one coupling shared by all pairs) or an array of length n_pairs. beta > 0 makes occupied neighbours less likely.
h – one-body field (fugacity): a scalar (shared), an (N,) array, or None (no field). It sets the mean density when the Hamiltonian does not conserve the particle number; at fixed particle number it is a constant and can be left out.
d – on-site (Gutzwiller) term: a scalar (shared), an (N,) array, or None (no term). d > 0 suppresses multiple occupancy.
dtype – JAX dtype of the amplitudes.
incremental – carry the weighted neighbour sums of every walker through the run and update them after each move instead of rebuilding them from all pairs at every rate evaluation (rebuilt every 100 chain steps against rounding drift). Not bit-identical to the default; the physics is unchanged.

The guide keeps every parameter in the form it was given, and train_with_netket() trains exactly those: a scalar stays a single tied value, None stays absent.

property beta: Coupling of every pair, (n_pairs,), in constructor order.

property h: One-body field per site, (N,) (zeros when constructed with h=None).

property d: On-site term per site, (N,) (zeros when constructed with d=None).

log_psi(config, cache=None): log ψ_T of one (N,) configuration or an (n, N) batch.

build_cache(op): The walker-independent tables of log_diff_from_operator() for a BoseHubbard operator, as a JastrowBoseHubbardOpCache (None for any other operator). Memoised per operator layout and parameter values.

init_state(config, cache=None): The per-walker state of the incremental mode: the weighted neighbour sums of config (weight -beta per pair, or 1 for a tied beta), shape (N,).

update_state(state, config, sites, deltas, cache=None): state after the move that adds deltas to sites (shape (2,) each, as pygfmc.hamiltonian.BoseHubbard.move() returns): every neighbour j of a moved site s gets w_sj Δ_s. Needs the operator’s build_cache() result, which holds the neighbour tables.

log_diff_from_operator(config, op, cache=None, state=None)

log ψ_T(x') - log ψ_T(x) for every transition of op.

Fast path for BoseHubbard operators (entry 0, the diagonal, is 0), from the tables of cache when it is the operator’s build_cache() result and built on the fly otherwise. state are the walker’s neighbour sums from init_state() / update_state() (rebuilt here when None). Any other operator falls back to its successor configurations.

save(path, group='jastrow_bose_hubbard', **attrs)

Write the guide to an HDF5 file (created or overwritten).

Stores the pairs, the parameters and the form of each one, so load() rebuilds the same parametrisation. Extra keyword arguments are stored as scalar group attributes.

classmethod load(path, group='jastrow_bose_hubbard', dtype=jnp.float32): Load a guide written by save().

class pygfmc.wavefunction.JastrowSpin(N, pairs, beta, h=None, d=None, dtype=jnp.float32, incremental=False)

Jastrow guide for spin-S configurations.

\[\log\psi_T(x) = -\sum_{(i,j)} \beta_{ij} x_i x_j + \sum_i h_i x_i - \sum_i d_i x_i^2\]

Parameters:

N – number of sites.
pairs – (n_pairs, 2) site pairs of the two-body term, e.g. the Hamiltonian edges.
beta – a scalar (one coupling shared by all pairs) or an array of length n_pairs. beta < 0 favours aligned neighbours.
h – one-body field: a scalar (shared), an (N,) array, or None (no field).
d – single-ion term: a scalar (shared), an (N,) array, or None (no term). d > 0 favours small |x_i|.
dtype – JAX dtype of the amplitudes.
incremental – carry the weighted neighbour sums W_i = Σ_j β_ij x_j of every walker through the run and update them after each move (two sites at most) instead of rebuilding them from all pairs at every rate evaluation. The sums are rebuilt every 100 chain steps, so float rounding cannot accumulate. The run is then no longer bit-identical to the default, which rebuilds the sums every time; the physics is unchanged.

The guide keeps every parameter in the form it was given, and train_with_netket() trains exactly those: a scalar stays a single tied value, None stays absent.

property beta: Coupling of every pair, (n_pairs,), in constructor order.

property h: One-body field per site, (N,) (zeros when constructed with h=None).

property d: Single-ion term per site, (N,) (zeros when constructed with d=None).

log_psi(config, cache=None): log ψ_T of one (N,) configuration or an (n, N) batch.

build_cache(op): The walker-independent tables of log_diff_from_operator() for a Spin operator, as a JastrowSpinOpCache (None for any other operator). Memoised per operator layout and parameter values, so repeated calls with the same operator skip the rebuild.

init_state(config, cache=None): The per-walker state of the incremental mode: the weighted neighbour sums W_i = Σ_j β_ij x_j of config, shape (N,).

update_state(state, config, sites, deltas, cache=None): state after the move that adds deltas to sites (shape (2,) each, as pygfmc.hamiltonian.Spin.move() returns): every neighbour j of a moved site s gets β_sj Δ_s. Needs the operator’s build_cache() result, which holds the neighbour tables.

log_diff_from_operator(config, op, cache=None, state=None)

log ψ_T(x') - log ψ_T(x) for every transition of op.

Fast path for Spin operators (entry 0, the diagonal, is 0), from the tables of cache when it is the operator’s build_cache() result and built on the fly otherwise. state are the walker’s neighbour sums from init_state() / update_state() (rebuilt here when None). Any other operator falls back to its successor configurations.

save(path, group='jastrow_spin', **attrs)

Write the guide to an HDF5 file (created or overwritten).

Stores the pairs, per-pair beta, h and d, and the form of each parameter, so load() rebuilds the same parametrisation. Extra keyword arguments are stored as scalar group attributes.

classmethod load(path, group='jastrow_spin', dtype=jnp.float32): Load a guide written by save().

pygfmc.wavefunction.train_with_netket(hamiltonian, guide, local_hilbert_space=None, N=None, sampler=None, initial_states=None, n_iter=100, learning_rate=0.02, n_samples=1024, n_chains=16, n_discard_per_chain=8, sr=False, diag_shift=1e-3, qgt=None, seed=0, show_progress=False, return_history=False)

Train guide on hamiltonian with NetKet VMC and return a new guide.

Needs only guide.log_psi plus the guide’s pytree registration, and hamiltonian.get_rates. Nothing is mutated: a new guide of the same class comes back, with parameters cast to their original dtypes.

Parameters:

hamiltonian – Any object exposing M_total and get_rates.
guide – Any registered-pytree guide exposing log_psi.
local_hilbert_space – Per-site values. Defaults to hamiltonian.local_hilbert_space.
N – Number of sites. Defaults to hamiltonian.N.
sampler – A ready-made NetKet sampler. It carries its own Hilbert space, so local_hilbert_space/N are then ignored. The default is MetropolisExchange, which swaps the values of two sites and so conserves their sum (S^z_tot for spins, the particle number for bosons); it is refused, with a ValueError, for a Hamiltonian that changes this sum. It samples one sector: that of initial_states when given, else the balanced one (S^z_tot = 0, or one particle per site), and swaps along the Hamiltonian’s edges (all site pairs if it has none). For a Hamiltonian that does not conserve this sum, e.g. a transverse field or a pair term, pass nk.sampler.MetropolisLocal(hilbert), which changes single sites. For local spaces with more than two values, an explicit ergodic sampler is required: swaps preserve the full value histogram, not just the conserved sum. Use suitable transfer proposals for a constrained higher-spin/boson space.
initial_states – (n, N) configurations seeding the Metropolis chains, resampled to the chain count. This is the argument that changes results rather than just cost: chains started from a random product state can take very long to random-walk into a constrained manifold such as the spin-ice 2-in-2-out states, and they must lie in the sector of interest to begin with.
n_iter – VMC settings.
learning_rate – VMC settings.
n_samples – VMC settings.
n_chains – VMC settings.
n_discard_per_chain – VMC settings.
sr – Stochastic reconfiguration, its regularisation and the quantum-geometric-tensor implementation. sr defaults to False because NetKet’s default on-the-fly QGT can fail with recent JAX versions (JAX registers Tracer as a jax.Array, so plum cannot disambiguate NetKet’s two is_scalar overloads and SR dies with AmbiguousLookupError inside qgt_onthefly_logic). With NetKet 3.21, qgt=nk.optimizer.qgt.QGTJacobianDense works, while QGTJacobianPyTree and the default fail. The optimiser is optax.sgd(learning_rate).
diag_shift – Stochastic reconfiguration, its regularisation and the quantum-geometric-tensor implementation. sr defaults to False because NetKet’s default on-the-fly QGT can fail with recent JAX versions (JAX registers Tracer as a jax.Array, so plum cannot disambiguate NetKet’s two is_scalar overloads and SR dies with AmbiguousLookupError inside qgt_onthefly_logic). With NetKet 3.21, qgt=nk.optimizer.qgt.QGTJacobianDense works, while QGTJacobianPyTree and the default fail. The optimiser is optax.sgd(learning_rate).
qgt – Stochastic reconfiguration, its regularisation and the quantum-geometric-tensor implementation. sr defaults to False because NetKet’s default on-the-fly QGT can fail with recent JAX versions (JAX registers Tracer as a jax.Array, so plum cannot disambiguate NetKet’s two is_scalar overloads and SR dies with AmbiguousLookupError inside qgt_onthefly_logic). With NetKet 3.21, qgt=nk.optimizer.qgt.QGTJacobianDense works, while QGTJacobianPyTree and the default fail. The optimiser is optax.sgd(learning_rate).
seed – Seeds parameter init and the sampler.
show_progress – Show NetKet’s progress bar.
return_history – Also return the per-iteration energy array.

Returns:

The trained guide, or (guide, history) when return_history.

pygfmc.wavefunction.make_hilbert(local_hilbert_space, N, total=None)

A NetKet Hilbert space whose local states are exactly these values.

Returns nk.hilbert.HomogeneousHilbert(StaticRange(...), N). Requires an evenly spaced ladder, which every convention in use is ((-1, 1), (-S, ..., S), (0, ..., n_max)). With total, the space is restricted to the configurations whose site values sum to total.

API reference

pygfmc.protocols

pygfmc.walkers

pygfmc.engine

Building blocks

pygfmc.hamiltonian

Hamiltonians

Measurement operators

Lattices

Exact diagonalization

pygfmc.wavefunction

`pygfmc.protocols`

`pygfmc.walkers`

`pygfmc.engine`

`pygfmc.hamiltonian`

`pygfmc.wavefunction`