wai-quantum 0.3.40

//! Neural quantum states — `wai.quantum.nqs`.
//!
//! A neural network can stand in for a many-body wavefunction: it gives the
//! amplitude of each configuration, and variational Monte Carlo trains it
//! toward the ground state. Results from such networks are offered as
//! evidence about systems no exact method reaches. This module is the method
//! in its checkable form, deterministic and held to exact answers.
//!
//! - **Models.** [`SpinModel`] is a spin-½ Hamiltonian on any graph, with
//!   exchange, Ising, and field terms. Presets are Heisenberg, the frustrated
//!   J₁–J₂ model, and transverse-field Ising, on chains and square lattices.
//! - **The network.** [`Rbm`] is a restricted Boltzmann machine with `α·n`
//!   hidden units, its parameters real or complex ([`Amplitude`]).
//!   - A real network is positive. On a bipartite lattice it can carry
//!     Marshall's sign, which makes the Heisenberg ground state positive in
//!     what remains.
//!   - A complex network carries a phase, so training can find a sign
//!     structure nobody supplied. Frustrated models need that, because no
//!     sign rule is known for them.
//!   - Weights and biases start small and random. With the biases at zero,
//!     the network would be exactly even under a global spin flip.
//! - **The sector.** [`Symmetry`] projects the network,
//!   `Ψ(s) = Σ_g χ(g) ψ(g·s)`, over a lattice's translations (momentum 0 or
//!   π along each direction), on a square lattice its whole space group
//!   ([`Symmetry::square`]: the translations with the eight rotations and
//!   reflections, in one of their four one-dimensional representations,
//!   [`PointIrrep`]), on any lattice the translations with its half turn and
//!   axis reflections ([`Symmetry::rectangle`], [`RectIrrep`]), and,
//!   optionally, the flip of every spin (even or odd).
//!   [`Symmetry::check`] proves the elements a group and the characters
//!   multiplicative, without which the sum projects onto nothing in
//!   particular. Ratios and log-derivatives become weighted averages over the
//!   transformed configurations, with a cache that Metropolis moves keep up
//!   to date.
//! - **Variational Monte Carlo.** [`Vmc`] samples `|ψ|²` by Metropolis
//!   chains: single flips, or swaps of antiparallel spins at fixed
//!   magnetisation. Or it enumerates every configuration (at most 24 spins)
//!   for expectations with no Monte Carlo noise. With a guide
//!   ([`VmcConfig::guide`]) the chains sample `|ψ|²` with a floor under it
//!   and each sample is weighted back, so the configurations near a node,
//!   where the local energy is large, are visited. Each step is stochastic
//!   reconfiguration, `θ ← θ − η (S + λ)⁻¹ F`, solved in parameter space or
//!   in sample space, whichever is smaller; the two give the same step. When
//!   Cholesky loses a pivot to rounding, the step comes from the
//!   eigendecomposition instead, with eigenvalues below `10⁻¹²` of the
//!   largest dropped. A projected network can take its metric there:
//!   amplitudes that nearly cancel give log-derivatives of order `10⁸`.
//! - **Referees.** [`SpinModel::exact_ground_energy`] is Lanczos on the
//!   amplitudes of a magnetisation sector, or all `2ⁿ` (at most 22 spins).
//!   [`SpinModel::exact_ground_energy_in`] is the same within one symmetry
//!   sector, on the sector's own basis (one state per orbit of the group,
//!   up to 24 spins), so a projected network is held to its own sector's
//!   answer. It refuses a group that is not one and a model the group does
//!   not preserve ([`SpinModel::commutes_with`]). [`SpinModel::opsum`] gives the model to
//!   DMRG on its exact MPO, for lattices past exact diagonalisation.
//!
//! Chains draw from their own seeded generators, and every sum runs in a
//! fixed order, so a run is the same bits at any thread count.
//!
//! # Checked
//!
//! - **The parts.** The model's exact diagonalisation equals the chain
//!   referees. Its operator sum's MPO equals its dense matrix on a 3 × 3
//!   lattice with every kind of term. The network's amplitude ratios equal
//!   ratios of its amplitudes to `10⁻¹²`, Marshall sign included, and its
//!   log-derivatives agree with central differences. Both spaces take the
//!   same reconfiguration step.
//! - **Exact expectations.** On an 8-site Heisenberg chain, 200 steps with
//!   every configuration enumerated reach the exact energy to `10⁻⁷`.
//! - **Sampling.** Sampled energies agree with enumerated ones within their
//!   error bars, and steps repeat bit for bit at any thread count.
//! - **Complex networks.** A complex network with real parameters is the real
//!   one, to `10⁻¹³`. Its ratios equal ratios of its amplitudes. Its
//!   holomorphic log-derivatives agree with differences along a real step.
//!   Both spaces take the same Hermitian reconfiguration step.
//! - **The projection.** It equals the explicit sum over the group, ratios
//!   and log-derivatives included. It carries its sector: a translation
//!   multiplies `Ψ` by its character, a global flip by its parity, and on the
//!   4 × 4 square a quarter turn, an axis reflection and a diagonal one by
//!   their characters in each of the four representations. Through accepted
//!   moves, the cache stays equal to a fresh one.
//! - **The groups.** Every square space group on 3 × 3 and 4 × 4, at both
//!   momenta, in every representation, with and without the flip, closes and
//!   multiplies its characters. A character that does not multiply, a set
//!   that does not close, and a character other than `±1` are each refused.
//! - **Guided sampling.** On an 8-site ring, a projected network's guided
//!   estimates agree with enumeration: the energy within five standard
//!   errors, the variance within 20 %. Steps repeat bit for bit at any thread
//!   count.
//! - **Chains start where the sector has amplitude.** A configuration fixed
//!   by an element of character −1 has none, and no move from it has a
//!   ratio; a chain is redrawn until its start is not one. With 32 chains at
//!   momentum (π, π) in A₁ on the 4 × 4 lattice, where a quarter of the
//!   configurations are such, every chain starts with amplitude.
//! - **The start.** On the 8-site ring, a projected network within `±0.01`
//!   starts on the ferromagnet (`E = n/4`, variance below `10⁻³`) and is
//!   within `10⁻²` of it after 40 steps; from `±0.1` it reaches the ground
//!   state to `10⁻³` in 200.
//! - **Two sector referees agree.** The orbit basis gives the same energy,
//!   to `10⁻⁹`, as Lanczos on `P H P` re-projected after every product
//!   (kept in the tests): in all 24 sectors of the 4 × 3 J₁–J₂ torus (four
//!   representations, both momenta along x, no flip, an even one and an odd
//!   one), on a 10-site ring, and on an Ising ring in a transverse field,
//!   where the magnetisation is not conserved. With the identity alone it is
//!   the plain diagonalisation; with the flip alone one half holds the
//!   ground state and the other lies above. Every rectangle group on 4 × 3
//!   and 6 × 4 closes and multiplies its characters.
//! - **The sector referee.** On a 10-site ring the sector of momentum π and
//!   odd parity holds the ground state to `10⁻⁹`, and the sector of momentum
//!   0 and even parity lies above it. On the 4 × 4 lattice at `J₂ = J₁/2`,
//!   the sector even under every rotation and reflection at momentum zero
//!   holds the ground state to `10⁻⁹`; the one odd under a quarter turn lies
//!   above. Open boundaries under ring translations, a flip out of a nonzero
//!   magnetisation, and a `σᶻ` field under the flip are refused.
//! - **Rings in their sector.** With every configuration enumerated, `α = 1`,
//!   reaches each ring's exact energy to rounding:
//!
//!   | ring (sector)                          | steps | relative error |
//!   |----------------------------------------|-------|----------------|
//!   | 6 sites, Heisenberg (π, odd)            | 150   | `< 10⁻⁸`       |
//!   | 10 sites, Heisenberg (π)                | 300   | 7·10⁻¹⁴        |
//!   | 10 sites, Majumdar–Ghosh (π)            | 300   | 3·10⁻¹⁵        |
//!   | 10 sites, J₁–J₂ 0.3 (π)                 | 300   | 3·10⁻¹⁴        |
//!   | 12 sites, Heisenberg (0, even)          | 400   | 7·10⁻¹⁶        |
//!   | 12 sites, Majumdar–Ghosh (0, even)      | 400   | 4·10⁻¹⁶        |
//!   | 12 sites, J₁–J₂ 0.3 (0, even)           | 400   | 1·10⁻¹⁵        |
//! - **A sign nobody supplied.** Exact diagonalisation gives the
//!   Majumdar–Ghosh ring (`J₂ = J₁/2`) its exact `−3n/8`. On 8-site rings
//!   with no sign rule given, Heisenberg and Majumdar–Ghosh, a real network
//!   stays more than 10 % off. A complex one comes within `2·10⁻³`, with
//!   every configuration enumerated. On 12-site rings, 400 steps reach:
//!
//!   | ring (no sign given)  | relative error |
//!   |-----------------------|----------------|
//!   | Majumdar–Ghosh        | 2·10⁻⁶         |
//!   | J₁–J₂, `J₂ = 0.3`     | 3.5·10⁻⁴       |
//!   | Heisenberg            | 2·10⁻³         |
//! - **Against exact answers** (`nqs_ground` example; relative error of the
//!   final estimate, with its statistical error):
//!
//!   | model             | lattice      | reference               | α | steps | relative error          |
//!   |-------------------|--------------|-------------------------|---|-------|-------------------------|
//!   | Heisenberg        | chain of 12  | exact diagonalisation   | 2 | 300   | 7.7·10⁻⁵ ± 7·10⁻⁵       |
//!   | Heisenberg        | 4 × 4 open   | exact diagonalisation   | 2 | 600   | 4.6·10⁻⁴ ± 2·10⁻⁴       |
//!   | Ising, `h = 3.044`| 4 × 4 open   | exact diagonalisation   | 1 | 400   | −3·10⁻⁶ ± 1·10⁻⁵        |
//!   | Heisenberg        | 6 × 6 open   | DMRG, bond 256          | 2 | 800   | 3.6·10⁻³ ± 3·10⁻⁴       |
//!
//!   Each run used 1000 samples a step. The 6 × 6 DMRG reference is itself
//!   an upper bound (discarded weight `2·10⁻⁵`). The network was still
//!   descending when it stopped there.
//!
//! # Measured
//!
//! On 10-site rings every model above settled about 9 % high, on a state
//! with essentially zero variance: an eigenstate, but not the ground state.
//! The cause is a symmetry the start imposes.
//!
//! - With its biases at zero, the network is exactly even under flipping
//!   every spin. Training keeps it even, because the biases' gradients vanish
//!   on an even state.
//! - A singlet is even or odd under that flip as half the ring's length is
//!   even or odd. On 10-site rings the singlet ground state is therefore out
//!   of reach, and training finds the lowest even state instead.
//! - Starting the biases away from zero frees it. The same 10-site runs then
//!   reach `3·10⁻⁴` (Heisenberg), `1·10⁻⁶` (Majumdar–Ghosh) and `1.4·10⁻⁵`
//!   (J₁–J₂ 0.3) in 400 steps.
//! - On 12-site rings, where the ground state is even, the even start did
//!   better than drawn biases (`2·10⁻³` against `2·10⁻²` for Heisenberg).
//!   It was a sector chosen by accident.
//!
//! Biases are now drawn, and a sector is chosen on purpose with
//! [`Symmetry`]. That is the table of rings above, where projection takes
//! every case to rounding.
//!
//! The variance cannot catch the trap, since an excited eigenstate has none.
//! Only a referee can.
//!
//! A projected network has a trap of its own at the start. Parameters near
//! zero make it nearly the projection of a constant, the uniform
//! superposition, which at zero magnetisation is the ferromagnet's member:
//! `E = n/4` on a Heisenberg ring, an eigenstate, so its energy gradient
//! vanishes. On the 8-site ring (every configuration enumerated, 200 steps)
//! a `±0.01` start ends `2.8·10⁻²` above the ground state, `±0.1` ends
//! `2.6·10⁻⁴` above, and `±0.3` `5·10⁻⁷`. How long a small start stays depends
//! on the lattice and the group: on the 4 × 4 lattice the translations
//! escaped it within 60 steps and ended the same from either start, while
//! the space group from `±0.01` ended ten times further off than from `±0.1`
//! (the table under the hardest point below). Projected networks in the CLI
//! and `nqs_ground` now start at `±0.1`.
//!
//! A run's sampled variance is not its variance. The space-group network
//! started at `±0.01` reported `1.4·10⁻²` from its final 8000 samples; over
//! every configuration its variance is `0.55`. 89 % of that sits in two
//! configurations related by symmetry, each of probability `2.6·10⁻⁶`
//! (expected 0.02 times in 8000 samples), where the local energy is 310
//! away from the mean: amplitudes near a node of the network, where the
//! local energy is a large number over a small one.
//!
//! - It is not the chains. Thinning them tenfold reported `1.9·10⁻²`;
//!   64 000 samples happened on one of the two and reported `1.55`.
//! - The network started at `±0.1` is the same: 97 % of its variance
//!   (`0.075`) sits in configurations whose probabilities total
//!   `2.4·10⁻³`.
//! - The energy holds better, but leans. Across both networks and the three
//!   samplings, every sampled energy came out below the enumerated one, by
//!   at most `4·10⁻⁴` relative: the tail that goes unsampled lies above the
//!   mean. A sampled energy is not a bound.
//!
//! A guide brings the tail into view. Sampling `|ψ|²` with a floor under it
//! ([`VmcConfig::guide`]) and weighting each sample back, the `±0.1`
//! network's sampled variance came within 5 % of the enumerated one at
//! every floor from `0.01` to `1`, four seeds each; without the guide its
//! four seeds gave from `0.11` to `2.0` times it. The `±0.01` network's
//! variance sits in two configurations of 12 870, which 8000 samples cannot
//! visit reliably whatever they are drawn from: guided seeds still gave from
//! `0.1` to `2.3` times it.
//!
//! Training sees what its samples see. Unguided, the `±0.1` space-group
//! network stopped near `10⁻³`, with nearly all its remaining variance in
//! configurations whose probabilities total `2.4·10⁻³`. Guided at `0.1`,
//! the same network, start and steps reach the ground energy to
//! `1.7·10⁻¹²`. Projected networks in the CLI and `nqs_ground` now sample
//! with a guide of `0.1`. Unprojected ones gain nothing from it (on the
//! 6 × 6 open Heisenberg lattice, 800 steps at `α = 2`, it doubled the
//! error) and sample `|ψ|²`.
//!
//! The same trap caught the referee once. Lanczos stopped a restart once its
//! lowest Ritz pair had converged, which is safe from DMRG's warm starts.
//! From a cold start that pair can be an excited state, and the exact
//! diagonalisation of the 8-site Majumdar–Ghosh ring returned `−2.595`
//! instead of `−3`. The network's energy, a true upper bound, came out below
//! it. The referees now build every Krylov vector.
//!
//! # Honest boundaries
//!
//! - **The hardest point, on a small sector.** At `J₂ = 0.5` on the
//!   periodic 4 × 4 lattice, exact diagonalisation gives `−8.45792335`, and
//!   the sector referee places it in the sector even under every rotation
//!   and reflection at momentum zero. The next state in that sector is 13 %
//!   higher.
//!   - Unprojected at `α = 1`, 600 steps of 1000 samples leave every network
//!     about 3 % above it (sampled estimates): complex with Marshall's sign,
//!     complex without, and real with it. A guide does not change that
//!     (complex, `2.78·10⁻²` over every configuration).
//!   - Projected, a complex `α = 1` network after 600 steps of 1000 samples,
//!     its energy then taken over every configuration (`nqs_ground
//!     --exact-final`):
//!
//!     | projection, even parity    | start   | guide | relative error |
//!     |----------------------------|---------|-------|----------------|
//!     | translations               | `±0.01` | none  | `1.87·10⁻²`    |
//!     | translations               | `±0.1`  | none  | `1.90·10⁻²`    |
//!     | space group, A₁            | `±0.01` | none  | `2.09·10⁻²`    |
//!     | space group, A₁            | `±0.1`  | none  | `1.27·10⁻³`    |
//!     | space group, A₁            | `±0.3`  | none  | `1.09·10⁻³`    |
//!     | space group, A₁            | `±0.1`  | `0.1` | `1.7·10⁻¹²`    |
//!     | space group, A₁            | `±0.1`  | `1`   | `1.2·10⁻¹²`    |
//!
//!     Unguided from `±0.1`, 2000 steps reach `4.60·10⁻⁴`, and `α = 2`
//!     reaches `6.91·10⁻⁴` in 600. Guided, the variance over every
//!     configuration is below `4·10⁻¹¹`.
//!   - The sector is small: 107 configurations up to its symmetries, against
//!     the network's 288 parameters. Reaching its ground state to rounding
//!     says the training finds what the network can hold. It says nothing of
//!     a sector larger than the network, which the next point is.
//! - **A sector larger than the network.** On the periodic 6 × 4 lattice at
//!   `J₂ = J₁/2`, exact diagonalisation gives `−12.54059845`, in the sector
//!   even under its half turn and both reflections at momentum zero, with
//!   even parity: 15 578 configurations up to its symmetries, against a
//!   complex `α = 1` network's 624 parameters. After 600 steps of 1000
//!   samples, energies over every configuration:
//!
//!   | projection, even parity | guide | relative error |
//!   |-------------------------|-------|----------------|
//!   | translations            | `0.1` | `1.06·10⁻³`    |
//!   | rectangle group, A₁     | none  | `8.32·10⁻³`    |
//!   | rectangle group, A₁     | `0.1` | `5.57·10⁻⁴`    |
//!
//!   The guide is worth a factor of fifteen here and the reflections a
//!   factor of two. What remains is the network's.
//! - **One-dimensional representations only.** Characters are `±1`: momenta
//!   0 and π, and the four one-dimensional representations of the square's
//!   rotations and reflections. The two-dimensional one and every other
//!   momentum are out of reach. Projection multiplies the cost of every
//!   amplitude by the group's size, 256 on the 4 × 4 lattice with the flip;
//!   600 steps took about seven times as long as with the translations.
//! - **Error bars ignore autocorrelation and the tail.** The reported error
//!   is the standard error of uncorrelated samples, from the sampled
//!   variance. Chains thinned by two sweeps are correlated, and the sampled
//!   variance can miss nearly all of the variance (see Measured), so the true
//!   error is larger. A guide recovers a tail spread over many
//!   configurations, not one held by a few. Only enumeration
//!   (`--exact-final`) settles it, up to 24 spins.
//! - **Dense reconfiguration.** `S` (or its sample-space twin) is formed and
//!   factored in full each step. That bounds the network and sample sizes
//!   long before iterative solvers would.

use crate::quantum_dmrg::{C, Op, OpSum, Scalar, default_threads, lanczos};
use crate::repro::{atan2, exp, ln, sin_cos};

// ---------------------------------------------------------------------------
// Spin models
// ---------------------------------------------------------------------------

/// A spin-½ model on any graph, in Pauli matrices:
/// `H = Σ J⊥ (σˣσˣ + σʸσʸ) + Σ J_z σᶻσᶻ + Σ hₓ σˣ + Σ h_z σᶻ`.
/// Spin `+1` is `|↑⟩`, the first basis state.
#[derive(Clone, Debug, PartialEq)]
pub struct SpinModel {
    pub n: usize,
    /// `(i, j, J⊥)`: `J⊥ (σˣᵢσˣⱼ + σʸᵢσʸⱼ)`, which swaps antiparallel spins with
    /// amplitude `2J⊥`.
    pub exchange: Vec<(usize, usize, f64)>,
    /// `(i, j, J_z)`: `J_z σᶻᵢσᶻⱼ`.
    pub zz: Vec<(usize, usize, f64)>,
    /// `(i, hₓ)`: `hₓ σˣᵢ`.
    pub x: Vec<(usize, f64)>,
    /// `(i, h_z)`: `h_z σᶻᵢ`.
    pub z: Vec<(usize, f64)>,
}

/// The bonds of an open (or periodic) chain.
pub fn chain_edges(n: usize, periodic: bool) -> Vec<(usize, usize)> {
    let mut e: Vec<(usize, usize)> = (0..n.saturating_sub(1)).map(|i| (i, i + 1)).collect();
    if periodic && n > 2 {
        e.push((n - 1, 0));
    }
    e
}

/// The nearest-neighbour bonds of an `lx × ly` square lattice, site
/// `y·lx + x`, open or periodic in both directions.
pub fn square_edges(lx: usize, ly: usize, periodic: bool) -> Vec<(usize, usize)> {
    let mut e = Vec::new();
    for y in 0..ly {
        for x in 0..lx {
            let i = y * lx + x;
            if x + 1 < lx {
                e.push((i, i + 1));
            } else if periodic && lx > 2 {
                e.push((i, y * lx));
            }
            if y + 1 < ly {
                e.push((i, i + lx));
            } else if periodic && ly > 2 {
                e.push((i, x));
            }
        }
    }
    e
}

/// The next-nearest-neighbour bonds of an open (or periodic) chain.
pub fn chain_next_nearest(n: usize, periodic: bool) -> Vec<(usize, usize)> {
    let mut e: Vec<(usize, usize)> = (0..n.saturating_sub(2)).map(|i| (i, i + 2)).collect();
    if periodic && n > 4 {
        e.push((n - 2, 0));
        e.push((n - 1, 1));
    }
    e
}

/// The diagonal (next-nearest) bonds of an `lx × ly` square lattice, open or
/// periodic in both directions.
pub fn square_diagonals(lx: usize, ly: usize, periodic: bool) -> Vec<(usize, usize)> {
    let mut e = Vec::new();
    let wrap = |v: isize, l: usize| -> Option<usize> {
        if v >= 0 && (v as usize) < l {
            Some(v as usize)
        } else if periodic && l > 2 {
            Some(v.rem_euclid(l as isize) as usize)
        } else {
            None
        }
    };
    for y in 0..ly {
        for x in 0..lx {
            let i = y * lx + x;
            for dx in [1isize, -1] {
                if let (Some(xx), Some(yy)) = (wrap(x as isize + dx, lx), wrap(y as isize + 1, ly)) {
                    let j = yy * lx + xx;
                    if i != j && !e.contains(&(i, j)) && !e.contains(&(j, i)) {
                        e.push((i, j));
                    }
                }
            }
        }
    }
    e
}

/// The checkerboard sublattice of an `lx × ly` lattice (a chain is `ly = 1`).
pub fn checkerboard(lx: usize, ly: usize) -> Vec<bool> {
    (0..lx * ly).map(|i| (i % lx + i / lx).is_multiple_of(2)).collect()
}

impl SpinModel {
    /// The Heisenberg model `J Σ Sᵢ·Sⱼ` (with `S = σ/2`) on `edges`.
    pub fn heisenberg(n: usize, edges: &[(usize, usize)], j: f64) -> SpinModel {
        SpinModel { n, exchange: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), zz: edges.iter().map(|&(a, b)| (a, b, j / 4.0)).collect(), x: Vec::new(), z: Vec::new() }
    }

    /// The frustrated `J₁ Σ_nn Sᵢ·Sⱼ + J₂ Σ_nnn Sᵢ·Sⱼ`. With `J₂ > 0` no
    /// sign rule is known in general, and a real network cannot follow it.
    pub fn j1j2(n: usize, nn: &[(usize, usize)], nnn: &[(usize, usize)], j1: f64, j2: f64) -> SpinModel {
        let mut m = SpinModel::heisenberg(n, nn, j1);
        let second = SpinModel::heisenberg(n, nnn, j2);
        m.exchange.extend(second.exchange);
        m.zz.extend(second.zz);
        m
    }

    /// The transverse-field Ising model `−J Σ σᶻσᶻ − h Σ σˣ` on `edges`.
    pub fn ising(n: usize, edges: &[(usize, usize)], j: f64, h: f64) -> SpinModel {
        SpinModel { n, exchange: Vec::new(), zz: edges.iter().map(|&(a, b)| (a, b, -j)).collect(), x: (0..n).map(|i| (i, -h)).collect(), z: Vec::new() }
    }

    /// Whether `Σ σᶻ` is conserved (no transverse field).
    pub fn conserves_magnetisation(&self) -> bool {
        self.x.is_empty()
    }

    /// The model as an operator sum, for DMRG on its MPO.
    pub fn opsum(&self) -> OpSum {
        let mut s = OpSum::new(self.n);
        for &(i, j, c) in &self.exchange {
            // σˣσˣ + σʸσʸ = 2(σ⁺σ⁻ + σ⁻σ⁺).
            s.add(2.0 * c, &[(i, Op::plus()), (j, Op::minus())]);
            s.add(2.0 * c, &[(i, Op::minus()), (j, Op::plus())]);
        }
        for &(i, j, c) in &self.zz {
            s.add(c, &[(i, Op::z()), (j, Op::z())]);
        }
        for &(i, c) in &self.x {
            s.add(c, &[(i, Op::x())]);
        }
        for &(i, c) in &self.z {
            s.add(c, &[(i, Op::z())]);
        }
        s
    }

    /// The states with total `Σ σᶻ = magnetisation` (all `2ⁿ` when none), in
    /// ascending order, and each one's place among them.
    fn sector_states(&self, magnetisation: Option<i32>) -> (Vec<usize>, Vec<u32>) {
        let n = self.n;
        let states: Vec<usize> = (0..1usize << n).filter(|&x| magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m)).collect();
        let mut place = vec![u32::MAX; 1usize << n];
        for (k, &x) in states.iter().enumerate() {
            place[x] = k as u32;
        }
        (states, place)
    }

    /// `out = H inp` on the amplitudes of `states` alone, which `H` must map
    /// among themselves. Site 0 is the most significant bit of a state, and
    /// a set bit is `|↓⟩`.
    fn apply_sector(&self, states: &[usize], place: &[u32], inp: &[f64], out: &mut [f64]) {
        let n = self.n;
        let bit = |i: usize| 1usize << (n - 1 - i);
        let spin = |x: usize, i: usize| if x & bit(i) == 0 { 1.0 } else { -1.0 };
        out.iter_mut().for_each(|o| *o = 0.0);
        for (k, (&x, &amp)) in states.iter().zip(inp).enumerate() {
            if amp == 0.0 {
                continue;
            }
            let mut diag = 0.0;
            for &(i, j, c) in &self.zz {
                diag += c * spin(x, i) * spin(x, j);
            }
            for &(i, c) in &self.z {
                diag += c * spin(x, i);
            }
            out[k] += diag * amp;
            for &(i, j, c) in &self.exchange {
                if spin(x, i) != spin(x, j) {
                    out[place[x ^ bit(i) ^ bit(j)] as usize] += 2.0 * c * amp;
                }
            }
            for &(i, c) in &self.x {
                out[place[x ^ bit(i)] as usize] += c * amp;
            }
        }
    }

    /// The exact ground energy by Lanczos (at most 22 sites), restricted to
    /// total `Σ σᶻ = magnetisation` when given (the model must conserve it).
    pub fn exact_ground_energy(&self, magnetisation: Option<i32>) -> Option<f64> {
        let n = self.n;
        if n == 0 || n > 22 {
            return None;
        }
        if magnetisation.is_some() && !self.conserves_magnetisation() {
            return None;
        }
        let (states, place) = self.sector_states(magnetisation);
        if states.is_empty() {
            return None;
        }
        let start: Vec<f64> = states.iter().map(|&x| 1.0 + ((x as u64).wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect();
        let apply = |inp: &[f64], out: &mut [f64]| self.apply_sector(&states, &place, inp, out);
        Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0)
    }

    /// The projector referee the orbit basis replaced, kept to check it:
    /// Lanczos on `P H P`, `P = (1/|G|) Σ_g χ(g) g`, re-projected after
    /// every product, on the magnetisation sector's amplitudes (at most 22
    /// sites).
    #[cfg(test)]
    fn projected_ground_energy(&self, magnetisation: Option<i32>, sym: &Symmetry) -> Option<f64> {
        let n = self.n;
        if n == 0 || n > 22 || sym.check().is_err() || !self.commutes_with(sym) {
            return None;
        }
        if magnetisation.is_some() && !self.conserves_magnetisation() {
            return None;
        }
        if magnetisation.is_some_and(|m| m != 0) && sym.flips.iter().any(|&f| f) {
            return None;
        }
        let bit = |i: usize| 1usize << (n - 1 - i);
        let (states, place) = self.sector_states(magnetisation);
        // images[k·|G| + g]: the place of g·states[k], whose spin at i is the
        // state's at perms[g][i], flipped if the element flips.
        let order = sym.len();
        let mut images = Vec::with_capacity(states.len() * order);
        for &x in &states {
            for (p, &flip) in sym.perms.iter().zip(&sym.flips) {
                let mut y = 0usize;
                for (i, &pi) in p.iter().enumerate() {
                    if (x & bit(pi) != 0) != flip {
                        y |= bit(i);
                    }
                }
                images.push(place[y]);
            }
        }
        let project = |v: &[f64], out: &mut [f64]| {
            for (k, o) in out.iter_mut().enumerate() {
                let row = &images[k * order..(k + 1) * order];
                let mut acc = 0.0;
                for (&y, &c) in row.iter().zip(&sym.chars) {
                    acc += c * v[y as usize];
                }
                *o = acc / order as f64;
            }
        };
        let raw: Vec<f64> = states.iter().map(|&x| 1.0 + ((x as u64).wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect();
        let mut start = vec![0.0; states.len()];
        project(&raw, &mut start);
        let (before, after) = (raw.iter().map(|v| v * v).sum::<f64>(), start.iter().map(|v| v * v).sum::<f64>());
        if after <= 1e-20 * before {
            return None;
        }
        let apply = |inp: &[f64], out: &mut [f64]| {
            let mut h = vec![0.0; states.len()];
            self.apply_sector(&states, &place, inp, &mut h);
            project(&h, out);
        };
        Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0)
    }

    /// Whether every element of `sym` maps the model onto itself: each bond
    /// and field to one of the same strength, and under a flip of every spin
    /// no `σᶻ` field, which the flip reverses.
    pub fn commutes_with(&self, sym: &Symmetry) -> bool {
        let n = self.n;
        if sym.perms.iter().any(|p| p.len() != n) {
            return false;
        }
        let bonds = |terms: &[(usize, usize, f64)], p: &[usize]| {
            let mut b: Vec<(usize, usize, u64)> = terms.iter().map(|&(i, j, c)| (p[i].min(p[j]), p[i].max(p[j]), c.to_bits())).collect();
            b.sort_unstable();
            b
        };
        let fields = |terms: &[(usize, f64)], p: &[usize]| {
            let mut f: Vec<(usize, u64)> = terms.iter().map(|&(i, c)| (p[i], c.to_bits())).collect();
            f.sort_unstable();
            f
        };
        let identity: Vec<usize> = (0..n).collect();
        let (exchange, zz, x, z) = (bonds(&self.exchange, &identity), bonds(&self.zz, &identity), fields(&self.x, &identity), fields(&self.z, &identity));
        sym.perms.iter().zip(&sym.flips).all(|(p, &flip)| bonds(&self.exchange, p) == exchange && bonds(&self.zz, p) == zz && fields(&self.x, p) == x && fields(&self.z, p) == z && !(flip && !self.z.is_empty()))
    }

    /// The exact ground energy within the sector `sym` projects onto (at most
    /// 24 sites), and in the magnetisation sector when given, by Lanczos on
    /// the sector's own basis: for each orbit of the group, the state
    /// `Σ_g χ(g) |g·r⟩` of its least configuration `r`, normalised, unless an
    /// element of character −1 fixes `r` and the state vanishes. `None` when
    /// `sym` is not a group with characters, when the model is not invariant
    /// under it, when a flip would leave the magnetisation sector, or when
    /// the sector is empty.
    pub fn exact_ground_energy_in(&self, magnetisation: Option<i32>, sym: &Symmetry) -> Option<f64> {
        let n = self.n;
        if n == 0 || n > 24 || sym.check().is_err() || !self.commutes_with(sym) {
            return None;
        }
        if magnetisation.is_some() && !self.conserves_magnetisation() {
            return None;
        }
        if magnetisation.is_some_and(|m| m != 0) && sym.flips.iter().any(|&f| f) {
            return None;
        }
        let bit = |i: usize| 1usize << (n - 1 - i);
        // g·x: its spin at i is x's at perms[g][i], flipped if g flips.
        let image = |g: usize, x: usize| -> usize {
            let mut y = 0usize;
            for (i, &pi) in sym.perms[g].iter().enumerate() {
                if (x & bit(pi) != 0) != sym.flips[g] {
                    y |= bit(i);
                }
            }
            y
        };
        // The least configuration of x's orbit, and the character of an
        // element taking x to it.
        let least = |x: usize| -> (usize, f64) {
            let (mut r, mut chi) = (x, 1.0);
            for g in 0..sym.len() {
                let y = image(g, x);
                if y < r {
                    r = y;
                    chi = sym.chars[g];
                }
            }
            (r, chi)
        };
        // Representatives, ascending, each with Σ χ over its stabiliser:
        // |G| times that is its state's squared norm, and zero when an
        // element of character −1 fixes it.
        let mut reps: Vec<usize> = Vec::new();
        let mut weight: Vec<f64> = Vec::new();
        for x in 0..1usize << n {
            if magnetisation.is_some_and(|m| n as i32 - 2 * x.count_ones() as i32 != m) {
                continue;
            }
            let mut stab = 0.0;
            let mut least_here = true;
            for g in 0..sym.len() {
                let y = image(g, x);
                if y < x {
                    least_here = false;
                    break;
                }
                if y == x {
                    stab += sym.chars[g];
                }
            }
            if least_here && stab > 0.0 {
                reps.push(x);
                weight.push(stab);
            }
        }
        if reps.is_empty() {
            return None;
        }
        // H on the basis, column by column: H|r⟩ = Σ_x h_xr |x⟩, and the
        // state of x's orbit carries it with χ(k) √(S_r′/S_r), x = k·r′.
        let mut columns: Vec<Vec<(usize, f64)>> = Vec::with_capacity(reps.len());
        for (a, &r) in reps.iter().enumerate() {
            let spin = |i: usize| if r & bit(i) == 0 { 1.0 } else { -1.0 };
            let mut diag = 0.0;
            for &(i, j, c) in &self.zz {
                diag += c * spin(i) * spin(j);
            }
            for &(i, c) in &self.z {
                diag += c * spin(i);
            }
            let mut col = vec![(a, diag)];
            let mut reach = |x: usize, h: f64| {
                let (r2, chi) = least(x);
                if let Ok(b) = reps.binary_search(&r2) {
                    col.push((b, h * chi * (weight[b] / weight[a]).sqrt()));
                }
            };
            for &(i, j, c) in &self.exchange {
                if spin(i) != spin(j) {
                    reach(r ^ bit(i) ^ bit(j), 2.0 * c);
                }
            }
            for &(i, c) in &self.x {
                reach(r ^ bit(i), c);
            }
            col.sort_by_key(|&(b, _)| b);
            let mut merged: Vec<(usize, f64)> = Vec::with_capacity(col.len());
            for (b, v) in col {
                match merged.last_mut() {
                    Some(last) if last.0 == b => last.1 += v,
                    _ => merged.push((b, v)),
                }
            }
            columns.push(merged);
        }
        let apply = |inp: &[f64], out: &mut [f64]| {
            out.iter_mut().for_each(|o| *o = 0.0);
            for (col, &v) in columns.iter().zip(inp) {
                for &(b, h) in col {
                    out[b] += h * v;
                }
            }
        };
        let start: Vec<f64> = reps.iter().map(|&x| 1.0 + ((x as u64).wrapping_mul(2_654_435_761) % 1000) as f64 / 1000.0).collect();
        Some(lanczos(&apply, &start, 60, 30, 1e-9, false).0)
    }
}

// ---------------------------------------------------------------------------
// The restricted Boltzmann machine
// ---------------------------------------------------------------------------

/// The numbers a network's parameters are: real, or complex so the network
/// can carry a sign structure nobody supplied.
pub trait Amplitude: Scalar {
    fn exp(self) -> Self;
    /// `ln cosh`, stable for large real parts (shifted by `ln 2`).
    fn lncosh(self) -> Self;
    fn tanh(self) -> Self;
    /// `self / d` for a real `d`.
    fn div_re(self, d: f64) -> Self;
    fn div(self, o: Self) -> Self;
    /// An initial weight, uniform in `±scale` (each part, when complex).
    fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> Self;
}

impl Amplitude for f64 {
    fn exp(self) -> f64 {
        exp(self)
    }
    fn lncosh(self) -> f64 {
        let a = self.abs();
        a + ln(1.0 + exp(-2.0 * a)) - core::f64::consts::LN_2
    }
    fn tanh(self) -> f64 {
        let e = exp(-2.0 * self.abs());
        let t = (1.0 - e) / (1.0 + e);
        if self < 0.0 { -t } else { t }
    }
    fn div_re(self, d: f64) -> f64 {
        self / d
    }
    fn div(self, o: f64) -> f64 {
        self / o
    }
    fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> f64 {
        scale * (2.0 * uniform() - 1.0)
    }
}

fn c_ln(z: C) -> C {
    C { re: ln((z.re * z.re + z.im * z.im).sqrt()), im: atan2(z.im, z.re) }
}

fn c_div(a: C, b: C) -> C {
    let d = b.re * b.re + b.im * b.im;
    C { re: (a.re * b.re + a.im * b.im) / d, im: (a.im * b.re - a.re * b.im) / d }
}

impl Amplitude for C {
    fn exp(self) -> C {
        let m = exp(self.re);
        let (s, c) = sin_cos(self.im);
        C { re: m * c, im: m * s }
    }
    fn lncosh(self) -> C {
        // cosh is even; with Re z ≥ 0, ln cosh z = z + ln(1 + e^{−2z}) − ln 2.
        if self.re < 0.0 {
            return Scalar::scale(self, -1.0).lncosh();
        }
        let e = Scalar::scale(self, -2.0).exp();
        let l = c_ln(C { re: 1.0 + e.re, im: e.im });
        C { re: self.re + l.re - core::f64::consts::LN_2, im: self.im + l.im }
    }
    fn tanh(self) -> C {
        if self.re < 0.0 {
            return Scalar::scale(Scalar::scale(self, -1.0).tanh(), -1.0);
        }
        let e = Scalar::scale(self, -2.0).exp();
        c_div(C { re: 1.0 - e.re, im: -e.im }, C { re: 1.0 + e.re, im: e.im })
    }
    fn div_re(self, d: f64) -> C {
        C { re: self.re / d, im: self.im / d }
    }
    fn div(self, o: C) -> C {
        c_div(self, o)
    }
    fn draw(uniform: &mut dyn FnMut() -> f64, scale: f64) -> C {
        let re = scale * (2.0 * uniform() - 1.0);
        let im = scale * (2.0 * uniform() - 1.0);
        C { re, im }
    }
}

/// One of the four one-dimensional representations of the square's eight
/// symmetries, by its character on a rotation and on a reflection.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PointIrrep {
    /// Even under every rotation and reflection.
    A1,
    /// Even under rotations, odd under every reflection.
    A2,
    /// Odd under a quarter turn and a diagonal reflection, even under a
    /// reflection in an axis.
    B1,
    /// Odd under a quarter turn and a reflection in an axis, even under a
    /// diagonal reflection.
    B2,
}

impl PointIrrep {
    /// The character of `turns` quarter turns after a reflection in the
    /// vertical axis if `reflect`. An odd number of turns after it makes a
    /// diagonal reflection, an even number an axis one.
    pub fn character(self, turns: usize, reflect: bool) -> f64 {
        let rotation = if turns % 2 == 1 { -1.0 } else { 1.0 };
        let reflection = if reflect { -1.0 } else { 1.0 };
        match self {
            PointIrrep::A1 => 1.0,
            PointIrrep::A2 => reflection,
            PointIrrep::B1 => rotation,
            PointIrrep::B2 => rotation * reflection,
        }
    }
}

/// One of the four one-dimensional representations of a rectangle's
/// symmetries (the half turn, and the reflections in its two axes), by its
/// character on each.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum RectIrrep {
    /// Even under the half turn and both reflections.
    A1,
    /// Even under the half turn, odd under both reflections.
    A2,
    /// Odd under the half turn and the reflection `y → −y`, even under
    /// `x → −x`.
    B1,
    /// Odd under the half turn and the reflection `x → −x`, even under
    /// `y → −y`.
    B2,
}

impl RectIrrep {
    /// The character of reflecting `x → −x` if `flip_x` and `y → −y` if
    /// `flip_y` (both: the half turn).
    pub fn character(self, flip_x: bool, flip_y: bool) -> f64 {
        let sx = if flip_x { -1.0 } else { 1.0 };
        let sy = if flip_y { -1.0 } else { 1.0 };
        match self {
            RectIrrep::A1 => 1.0,
            RectIrrep::A2 => sx * sy,
            RectIrrep::B1 => sy,
            RectIrrep::B2 => sx,
        }
    }
}

/// A symmetry group of the lattice: site permutations, each possibly
/// followed by flipping every spin, each with a real character. The
/// projection `Ψ(s) = Σ_g χ(g) ψ(g·s)` picks the sector: a momentum of 0 or
/// π along each direction, on a square lattice a representation of its
/// rotations and reflections, and with the flip even or odd parity.
#[derive(Clone, Debug, PartialEq)]
pub struct Symmetry {
    /// `perms[g][i]`: the site whose spin lands on site `i`.
    pub perms: Vec<Vec<usize>>,
    /// Whether element `g` also flips every spin.
    pub flips: Vec<bool>,
    /// `±1`, the character of each element.
    pub chars: Vec<f64>,
    /// `inverse[g][j]`: where site `j`'s spin lands.
    inverse: Vec<Vec<usize>>,
}

impl Symmetry {
    pub fn new(perms: Vec<Vec<usize>>, chars: Vec<f64>) -> Symmetry {
        let flips = vec![false; perms.len()];
        Symmetry::with_flips(perms, flips, chars)
    }

    fn with_flips(perms: Vec<Vec<usize>>, flips: Vec<bool>, chars: Vec<f64>) -> Symmetry {
        assert_eq!(perms.len(), chars.len());
        assert_eq!(perms.len(), flips.len());
        let inverse = perms
            .iter()
            .map(|p| {
                let mut inv = vec![usize::MAX; p.len()];
                for (i, &j) in p.iter().enumerate() {
                    inv[j] = i;
                }
                assert!(inv.iter().all(|&x| x != usize::MAX), "not a permutation");
                inv
            })
            .collect();
        Symmetry { perms, flips, chars, inverse }
    }

    /// Add the flip of every spin to the group, projecting onto even or odd
    /// parity under it. A singlet's parity is `(−1)^{n/2}`.
    pub fn with_spin_flip(self, odd: bool) -> Symmetry {
        let mut perms = self.perms.clone();
        let mut flips = self.flips.clone();
        let mut chars = self.chars.clone();
        for g in 0..self.perms.len() {
            perms.push(self.perms[g].clone());
            flips.push(!self.flips[g]);
            chars.push(if odd { -self.chars[g] } else { self.chars[g] });
        }
        Symmetry::with_flips(perms, flips, chars)
    }

    /// The translations of a ring of `n` sites, at momentum 0 or π (π needs
    /// `n` even: going round the ring must multiply by one).
    pub fn ring(n: usize, momentum_pi: bool) -> Symmetry {
        assert!(!momentum_pi || n.is_multiple_of(2), "momentum π on a ring of odd length");
        let perms = (0..n).map(|t| (0..n).map(|i| (i + t) % n).collect()).collect();
        let chars = (0..n).map(|t| if momentum_pi && t % 2 == 1 { -1.0 } else { 1.0 }).collect();
        Symmetry::new(perms, chars)
    }

    /// The translations of an `lx × ly` torus (site `y·lx + x`), at momentum
    /// 0 or π along each direction (π along one needs its length even).
    pub fn torus(lx: usize, ly: usize, kx_pi: bool, ky_pi: bool) -> Symmetry {
        assert!(!kx_pi || lx.is_multiple_of(2), "momentum π along an odd length");
        assert!(!ky_pi || ly.is_multiple_of(2), "momentum π along an odd length");
        let mut perms = Vec::new();
        let mut chars = Vec::new();
        for ty in 0..ly {
            for tx in 0..lx {
                perms.push((0..lx * ly).map(|i| ((i / lx + ty) % ly) * lx + (i % lx + tx) % lx).collect());
                let odd = (kx_pi && tx % 2 == 1) ^ (ky_pi && ty % 2 == 1);
                chars.push(if odd { -1.0 } else { 1.0 });
            }
        }
        Symmetry::new(perms, chars)
    }

    /// The whole space group of an `l × l` torus: every translation, each
    /// composed with the eight symmetries of the square about a site (quarter
    /// turns, and reflections in the axes and the diagonals). Momentum 0 or
    /// `(π, π)`, both of which every rotation and reflection leaves alone, so
    /// a character of the eight extends to the group; `irrep` picks which.
    pub fn square(l: usize, momentum_pi: bool, irrep: PointIrrep) -> Symmetry {
        // At l = 2 a reflection and the identity move every site alike, and
        // a character that tells them apart cannot be one.
        assert!(l >= 3, "the square's symmetries are distinct from l = 3");
        assert!(!momentum_pi || l.is_multiple_of(2), "momentum π along an odd length");
        let mut perms = Vec::new();
        let mut chars = Vec::new();
        for ty in 0..l {
            for tx in 0..l {
                for turns in 0..4 {
                    for reflect in [false, true] {
                        perms.push(
                            (0..l * l)
                                .map(|i| {
                                    let (mut x, mut y) = (i % l, i / l);
                                    if reflect {
                                        x = (l - x) % l;
                                    }
                                    for _ in 0..turns {
                                        (x, y) = ((l - y) % l, x);
                                    }
                                    ((y + ty) % l) * l + (x + tx) % l
                                })
                                .collect(),
                        );
                        let momentum = if momentum_pi && (tx + ty) % 2 == 1 { -1.0 } else { 1.0 };
                        chars.push(momentum * irrep.character(turns, reflect));
                    }
                }
            }
        }
        Symmetry::new(perms, chars)
    }

    /// The space group of an `lx × ly` torus with its rectangle's symmetries:
    /// every translation, each composed with the identity, the reflections
    /// `x → −x` and `y → −y` about a site, or both (the half turn), in one of
    /// their four one-dimensional representations. Momentum 0 or π along
    /// each direction, all of which the reflections leave alone.
    pub fn rectangle(lx: usize, ly: usize, kx_pi: bool, ky_pi: bool, irrep: RectIrrep) -> Symmetry {
        // At length 2 a reflection moves every site as the identity does.
        assert!(lx >= 3 && ly >= 3, "a rectangle's reflections are distinct from length 3");
        assert!(!kx_pi || lx.is_multiple_of(2), "momentum π along an odd length");
        assert!(!ky_pi || ly.is_multiple_of(2), "momentum π along an odd length");
        let mut perms = Vec::new();
        let mut chars = Vec::new();
        for ty in 0..ly {
            for tx in 0..lx {
                for flip_x in [false, true] {
                    for flip_y in [false, true] {
                        perms.push(
                            (0..lx * ly)
                                .map(|i| {
                                    let (mut x, mut y) = (i % lx, i / lx);
                                    if flip_x {
                                        x = (lx - x) % lx;
                                    }
                                    if flip_y {
                                        y = (ly - y) % ly;
                                    }
                                    ((y + ty) % ly) * lx + (x + tx) % lx
                                })
                                .collect(),
                        );
                        let odd = (kx_pi && tx % 2 == 1) ^ (ky_pi && ty % 2 == 1);
                        chars.push(if odd { -1.0 } else { 1.0 } * irrep.character(flip_x, flip_y));
                    }
                }
            }
        }
        Symmetry::new(perms, chars)
    }

    /// Whether the elements form a group under composition and every
    /// character is `±1` and multiplies as the elements do: what makes
    /// `Σ_g χ(g) ψ(g·s)` a projection onto one sector, not onto a mixture or
    /// onto nothing.
    pub fn check(&self) -> Result<(), String> {
        if self.is_empty() {
            return Err("no elements".into());
        }
        if let Some(c) = self.chars.iter().find(|&&c| c != 1.0 && c != -1.0) {
            return Err(format!("character {c} is not ±1"));
        }
        let mut index: std::collections::BTreeMap<(&[usize], bool), usize> = std::collections::BTreeMap::new();
        for g in 0..self.len() {
            match index.get(&(self.perms[g].as_slice(), self.flips[g])) {
                Some(&h) if self.chars[h] != self.chars[g] => return Err(format!("elements {h} and {g} act alike with different characters")),
                Some(_) => {}
                None => {
                    index.insert((self.perms[g].as_slice(), self.flips[g]), g);
                }
            }
        }
        let n = self.perms[0].len();
        let identity: Vec<usize> = (0..n).collect();
        match index.get(&(identity.as_slice(), false)) {
            Some(&e) if self.chars[e] == 1.0 => {}
            Some(_) => return Err("the identity's character is not 1".into()),
            None => return Err("no identity".into()),
        }
        // (g·(h·s))ᵢ = s at p_h(p_g(i)), flipped by either.
        let mut q = vec![0usize; n];
        for g in 0..self.len() {
            for h in 0..self.len() {
                for (i, qi) in q.iter_mut().enumerate() {
                    *qi = self.perms[h][self.perms[g][i]];
                }
                let Some(&k) = index.get(&(q.as_slice(), self.flips[g] ^ self.flips[h])) else {
                    return Err(format!("elements {g} and {h} compose to one outside the set"));
                };
                if self.chars[k] != self.chars[g] * self.chars[h] {
                    return Err(format!("the characters of {g} and {h} do not multiply to their product's"));
                }
            }
        }
        Ok(())
    }

    pub fn len(&self) -> usize {
        self.perms.len()
    }

    pub fn is_empty(&self) -> bool {
        self.perms.is_empty()
    }
}

/// What a configuration's amplitude ratios and derivatives need: `θ(s)`,
/// or for a projected network `θ` of every transformed configuration and its
/// share `p_g = χ(g) ψ(s∘g) / Ψ(s)` of the amplitude.
#[derive(Clone, Debug)]
enum Cache<T> {
    Plain(Vec<T>),
    Projected { configs: Vec<Vec<i8>>, thetas: Vec<Vec<T>>, shares: Vec<T>, log_abs: f64 },
}

/// A restricted Boltzmann machine: `ψ(s) = sign(s) · exp(Σ aᵢsᵢ) ·
/// Π_j cosh(b_j + Σ W_jᵢ sᵢ)`, with `m = α·n` hidden units. Its parameters are
/// real (`Rbm<f64>`), so `ψ` is positive up to the optional sign, or complex
/// (`Rbm<C>`), so `ψ` carries a phase the training finds. On a bipartite
/// lattice the optional sign is Marshall's, `(−1)` to the number of up spins
/// on sublattice A, which makes the Heisenberg ground state positive in the
/// rest.
#[derive(Clone, Debug, PartialEq)]
pub struct Rbm<T = f64> {
    pub n: usize,
    pub m: usize,
    pub a: Vec<T>,
    pub b: Vec<T>,
    /// `W[j·n + i]`.
    pub w: Vec<T>,
    pub sublattice: Option<Vec<bool>>,
    /// Project onto a momentum: `Ψ(s) = Σ_g χ(g) ψ(s∘g)`.
    pub symmetry: Option<Symmetry>,
}

/// A deterministic generator (splitmix64).
#[derive(Clone, Debug, PartialEq)]
struct Rng(u64);

impl Rng {
    fn next(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
        z ^ (z >> 31)
    }

    /// Uniform on `[0, 1)`.
    fn uniform(&mut self) -> f64 {
        (self.next() >> 11) as f64 / 9_007_199_254_740_992.0
    }

    fn below(&mut self, k: usize) -> usize {
        (self.next() % k as u64) as usize
    }
}

impl Rbm<f64> {
    /// A real machine: `α·n` hidden units, weights uniform in `±scale`,
    /// biases zero.
    pub fn new(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm {
        Rbm::random(n, alpha, scale, seed)
    }
}

impl Rbm<C> {
    /// A complex machine: weights with real and imaginary parts uniform in
    /// `±scale`, biases zero.
    pub fn complex(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm<C> {
        Rbm::random(n, alpha, scale, seed)
    }
}

impl<T: Amplitude> Rbm<T> {
    fn random(n: usize, alpha: usize, scale: f64, seed: u64) -> Rbm<T> {
        let m = alpha.max(1) * n;
        let mut rng = Rng(seed);
        let mut u = || rng.uniform();
        let w = (0..m * n).map(|_| T::draw(&mut u, scale)).collect();
        // Biases drawn too: with both at zero the machine is exactly even
        // under flipping every spin, and training keeps it so, which shuts
        // it out of every odd state (singlets on rings whose half-length is
        // odd among them).
        let a = (0..n).map(|_| T::draw(&mut u, scale)).collect();
        let b = (0..m).map(|_| T::draw(&mut u, scale)).collect();
        Rbm { n, m, a, b, w, sublattice: None, symmetry: None }
    }

    /// Multiply by Marshall's sign for this sublattice.
    pub fn with_marshall_sign(mut self, sublattice: Vec<bool>) -> Rbm<T> {
        assert_eq!(sublattice.len(), self.n);
        self.sublattice = Some(sublattice);
        self
    }

    /// Project onto the momentum `symmetry` selects.
    pub fn with_symmetry(mut self, symmetry: Symmetry) -> Rbm<T> {
        assert!(symmetry.perms.iter().all(|p| p.len() == self.n));
        self.symmetry = Some(symmetry);
        self
    }

    /// The configuration's cache.
    fn cache(&self, s: &[i8]) -> Cache<T> {
        let Some(sym) = &self.symmetry else {
            return Cache::Plain(self.theta(s));
        };
        let configs: Vec<Vec<i8>> = sym.perms.iter().zip(&sym.flips).map(|(p, &f)| p.iter().map(|&j| if f { -s[j] } else { s[j] }).collect()).collect();
        let thetas: Vec<Vec<T>> = configs.iter().map(|t| self.theta(t)).collect();
        let logs: Vec<T> = configs.iter().map(|t| self.log_psi(t)).collect();
        let top = logs.iter().map(|l| l.re()).fold(f64::NEG_INFINITY, f64::max);
        // χ(g) sign(s∘g) e^{ln ψ(s∘g) − top}, summed in order.
        let terms: Vec<T> = logs.iter().zip(&configs).zip(&sym.chars).map(|((l, t), &c)| l.sub(T::ONE.scale(top)).exp().scale(c * self.sign(t))).collect();
        let total = terms.iter().fold(T::ZERO, |acc, &x| acc.add(x));
        let mag2 = total.norm2();
        // A sum below the rounding error of its own terms is zero: the
        // configuration is one the symmetry forbids, and what is left is
        // noise. Its log-derivative, noise over noise, would otherwise feed
        // the metric something of order one.
        let scale: f64 = terms.iter().map(|x| x.norm2().sqrt()).sum();
        let vanishes = mag2.sqrt() <= 1e-13 * scale;
        let log_abs = if vanishes { f64::NEG_INFINITY } else { top + 0.5 * ln(mag2) };
        let shares = if vanishes { vec![T::ZERO; terms.len()] } else { terms.iter().map(|&x| x.div(total)).collect() };
        Cache::Projected { configs, thetas, shares, log_abs }
    }

    /// `Ψ(s′)/Ψ(s)` where `s′` is `s` with `flips` flipped.
    fn ratio_cached(&self, s: &[i8], cache: &Cache<T>, flips: &[usize]) -> T {
        match cache {
            Cache::Plain(theta) => self.ratio(s, theta, flips),
            Cache::Projected { configs, thetas, shares, .. } => {
                let sym = self.symmetry.as_ref().expect("a projected cache has a symmetry");
                let mut r = T::ZERO;
                let mut moved = Vec::with_capacity(flips.len());
                for g in 0..configs.len() {
                    if shares[g].is_zero() {
                        continue;
                    }
                    moved.clear();
                    moved.extend(flips.iter().map(|&j| sym.inverse[g][j]));
                    r = r.add(shares[g].mul(self.ratio(&configs[g], &thetas[g], &moved)));
                }
                r
            }
        }
    }

    /// `∂ ln Ψ / ∂(a, b, W)`.
    fn derivatives_cached(&self, s: &[i8], cache: &Cache<T>, out: &mut [T]) {
        match cache {
            Cache::Plain(theta) => self.log_derivatives(s, theta, out),
            Cache::Projected { configs, thetas, shares, .. } => {
                out.iter_mut().for_each(|o| *o = T::ZERO);
                let mut one = vec![T::ZERO; out.len()];
                for g in 0..configs.len() {
                    if shares[g].is_zero() {
                        continue;
                    }
                    self.log_derivatives(&configs[g], &thetas[g], &mut one);
                    for (o, &d) in out.iter_mut().zip(&one) {
                        *o = o.add(shares[g].mul(d));
                    }
                }
            }
        }
    }

    /// `ln |Ψ(s)|` from the cache (up to a constant).
    fn log_abs_cached(&self, s: &[i8], cache: &Cache<T>) -> f64 {
        match cache {
            Cache::Plain(_) => self.log_amplitude(s),
            Cache::Projected { log_abs, .. } => *log_abs,
        }
    }

    /// Flip `flips` in `s` and bring the cache along, given the ratio `r`
    /// already computed for the move.
    fn accept(&self, s: &mut [i8], cache: &mut Cache<T>, flips: &[usize], r: T) {
        let n = self.n;
        match cache {
            Cache::Plain(theta) => {
                for &i in flips {
                    let ds = -2.0 * f64::from(s[i]);
                    for (j, t) in theta.iter_mut().enumerate() {
                        *t = t.add(self.w[j * n + i].scale(ds));
                    }
                    s[i] = -s[i];
                }
            }
            Cache::Projected { configs, thetas, shares, log_abs } => {
                let sym = self.symmetry.as_ref().expect("a projected cache has a symmetry");
                for g in 0..configs.len() {
                    let moved: Vec<usize> = flips.iter().map(|&j| sym.inverse[g][j]).collect();
                    let rg = self.ratio(&configs[g], &thetas[g], &moved);
                    shares[g] = shares[g].mul(rg).div(r);
                    for &i in &moved {
                        let ds = -2.0 * f64::from(configs[g][i]);
                        for (j, t) in thetas[g].iter_mut().enumerate() {
                            *t = t.add(self.w[j * n + i].scale(ds));
                        }
                        configs[g][i] = -configs[g][i];
                    }
                }
                *log_abs += 0.5 * ln(r.norm2());
                for &i in flips {
                    s[i] = -s[i];
                }
            }
        }
    }

    /// The number of variational parameters (complex ones counted once).
    pub fn params(&self) -> usize {
        self.n + self.m + self.n * self.m
    }

    fn theta(&self, s: &[i8]) -> Vec<T> {
        (0..self.m)
            .map(|j| {
                let row = &self.w[j * self.n..(j + 1) * self.n];
                let mut t = self.b[j];
                for (wi, &si) in row.iter().zip(s) {
                    t = t.add(wi.scale(f64::from(si)));
                }
                t
            })
            .collect()
    }

    /// `ln ψ(s)` up to a constant, without the sign: real part `ln |ψ|`,
    /// imaginary part the phase.
    pub fn log_psi(&self, s: &[i8]) -> T {
        let mut l = T::ZERO;
        for (a, &si) in self.a.iter().zip(s) {
            l = l.add(a.scale(f64::from(si)));
        }
        for t in self.theta(s) {
            l = l.add(t.lncosh());
        }
        l
    }

    /// `ln |ψ(s)|` up to a constant.
    pub fn log_amplitude(&self, s: &[i8]) -> f64 {
        self.log_psi(s).re()
    }

    /// The sign of `ψ(s)` from Marshall's rule (one without it).
    pub fn sign(&self, s: &[i8]) -> f64 {
        match &self.sublattice {
            Some(sub) => {
                let ups = s.iter().zip(sub).filter(|&(&si, &a)| a && si > 0).count();
                if ups.is_multiple_of(2) { 1.0 } else { -1.0 }
            }
            None => 1.0,
        }
    }

    /// `ψ(s′)/ψ(s)` where `s′` is `s` with `flips` flipped, given `θ(s)`.
    fn ratio(&self, s: &[i8], theta: &[T], flips: &[usize]) -> T {
        let mut l = T::ZERO;
        for &i in flips {
            l = l.sub(self.a[i].scale(2.0).scale(f64::from(s[i])));
        }
        for (j, &t) in theta.iter().enumerate() {
            let mut d = T::ZERO;
            for &i in flips {
                d = d.sub(self.w[j * self.n + i].scale(2.0).scale(f64::from(s[i])));
            }
            l = l.add(t.add(d).lncosh().sub(t.lncosh()));
        }
        let mut r = l.exp();
        if let Some(sub) = &self.sublattice
            && flips.iter().filter(|&&i| sub[i]).count() % 2 == 1
        {
            r = r.scale(-1.0);
        }
        r
    }

    /// `∂ ln ψ / ∂(a, b, W)` at `s`, given `θ(s)`.
    fn log_derivatives(&self, s: &[i8], theta: &[T], out: &mut [T]) {
        let (n, m) = (self.n, self.m);
        for (o, &si) in out[..n].iter_mut().zip(s) {
            *o = T::ONE.scale(f64::from(si));
        }
        for (j, &th) in theta.iter().enumerate() {
            let t = th.tanh();
            out[n + j] = t;
            for (o, &si) in out[n + m + j * n..n + m + (j + 1) * n].iter_mut().zip(s) {
                *o = t.scale(f64::from(si));
            }
        }
    }

    fn shift(&mut self, delta: &[T], step: f64) {
        let (n, m) = (self.n, self.m);
        for (x, d) in self.a.iter_mut().zip(&delta[..n]) {
            *x = x.sub(d.scale(step));
        }
        for (x, d) in self.b.iter_mut().zip(&delta[n..n + m]) {
            *x = x.sub(d.scale(step));
        }
        for (x, d) in self.w.iter_mut().zip(&delta[n + m..]) {
            *x = x.sub(d.scale(step));
        }
    }
}

/// `E_loc(s) = Σ_{s′} H_{ss′} ψ(s′)/ψ(s)`.
fn local_energy<T: Amplitude>(model: &SpinModel, rbm: &Rbm<T>, s: &[i8], cache: &Cache<T>) -> T {
    let mut e = T::ZERO;
    for &(i, j, c) in &model.zz {
        e = e.add(T::ONE.scale(c * f64::from(s[i]) * f64::from(s[j])));
    }
    for &(i, c) in &model.z {
        e = e.add(T::ONE.scale(c * f64::from(s[i])));
    }
    for &(i, j, c) in &model.exchange {
        if s[i] != s[j] {
            e = e.add(rbm.ratio_cached(s, cache, &[i, j]).scale(2.0 * c));
        }
    }
    for &(i, c) in &model.x {
        e = e.add(rbm.ratio_cached(s, cache, &[i]).scale(c));
    }
    e
}

// ---------------------------------------------------------------------------
// Variational Monte Carlo
// ---------------------------------------------------------------------------

/// Settings for [`Vmc`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct VmcConfig {
    /// Samples per iteration, over all chains.
    pub samples: usize,
    pub chains: usize,
    /// Sweeps (`n` moves each) discarded when a chain starts.
    pub burn: usize,
    /// Sweeps between kept samples.
    pub thin: usize,
    /// Enumerate every configuration (at most 24 spins) instead of sampling:
    /// exact expectations, no Monte Carlo noise.
    pub exact: bool,
    /// Keep `Σ σᶻ` at this value: moves swap antiparallel spins. The model
    /// must conserve it.
    pub magnetisation: Option<i32>,
    /// Learning rate of stochastic reconfiguration.
    pub lr: f64,
    /// Diagonal shift added to the metric `S`.
    pub shift: f64,
    pub seed: u64,
    /// Threads over chains and configurations; the result does not depend on
    /// it.
    pub threads: usize,
    /// Sample `|ψ|² + guide·|ψ̄|²` instead of `|ψ|²`, `|ψ̄|` the chains'
    /// typical amplitude when the batch starts, and weight each sample by
    /// `|ψ|²` over that. Near a node of the network the local energy is a
    /// large number over a small one, and `|ψ|²` visits there too rarely for
    /// the variance to be measured; the floor visits, the weight corrects.
    /// Zero samples `|ψ|²`.
    pub guide: f64,
}

impl Default for VmcConfig {
    fn default() -> VmcConfig {
        VmcConfig { samples: 2000, chains: 8, burn: 50, thin: 2, exact: false, magnetisation: None, lr: 0.05, shift: 1e-3, seed: 1, threads: default_threads(), guide: 0.0 }
    }
}

/// One iteration's estimate.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Estimate {
    /// `Re ⟨E_loc⟩`.
    pub energy: f64,
    /// `⟨|E_loc − E|²⟩`: zero exactly at an eigenstate.
    pub variance: f64,
    /// The statistical error of `energy`, from the variance and the sample
    /// count, ignoring autocorrelation (zero when exact).
    pub error: f64,
    /// Fraction of proposed Metropolis moves accepted (one when exact).
    pub acceptance: f64,
}

/// Variational Monte Carlo with stochastic reconfiguration on an RBM, real
/// or complex.
#[derive(Clone, Debug, PartialEq)]
pub struct Vmc<T = f64> {
    pub model: SpinModel,
    pub rbm: Rbm<T>,
    pub cfg: VmcConfig,
    chains: Vec<(Vec<i8>, Rng)>,
}

/// Samples: their weights (summing to one), local energies and
/// log-derivatives (`params` per row).
struct Batch<T> {
    weights: Vec<f64>,
    energies: Vec<T>,
    derivs: Vec<T>,
    acceptance: f64,
    /// Samples drawn from a guide, weighted back to `|ψ|²`.
    guided: bool,
}

/// The guide `|ψ|² + floor·e^{2c}`, by its logarithm.
#[derive(Clone, Copy, Debug)]
struct Guide {
    /// The chains' mean `ln |ψ|` when the batch starts.
    c: f64,
    ln_floor: f64,
}

impl Guide {
    /// `ln (e^{2(l−c)} + floor)`.
    fn ln_q(self, l: f64) -> f64 {
        let x = 2.0 * (l - self.c);
        let (hi, lo) = if x > self.ln_floor { (x, self.ln_floor) } else { (self.ln_floor, x) };
        hi + ln(1.0 + exp(lo - hi))
    }

    /// A sample's weight `|ψ|² / (|ψ|² + floor·e^{2c})`.
    fn weight(self, l: f64) -> f64 {
        exp(2.0 * (l - self.c) - self.ln_q(l))
    }
}

impl<T: Amplitude> Vmc<T> {
    pub fn new(model: SpinModel, rbm: Rbm<T>, cfg: VmcConfig) -> Vmc<T> {
        assert_eq!(model.n, rbm.n, "the model and the machine differ in size");
        if cfg.magnetisation.is_some() {
            assert!(model.conserves_magnetisation(), "the model does not conserve the magnetisation");
        }
        let n = model.n;
        let chains = (0..cfg.chains.max(1))
            .map(|c| {
                let mut rng = Rng(cfg.seed ^ (c as u64).wrapping_mul(0xa076_1d64_78bd_642f));
                let mut s = start_configuration(n, cfg.magnetisation, &mut rng);
                // A configuration the projection forbids has no amplitude, so
                // no move from it has a ratio, and a chain started there would
                // never leave. Draw again; an unprojected network never does.
                let mut draws = 1;
                while !rbm.log_abs_cached(&s, &rbm.cache(&s)).is_finite() {
                    assert!(draws < 10_000, "the projection vanishes on every start drawn: an empty sector?");
                    s = start_configuration(n, cfg.magnetisation, &mut rng);
                    draws += 1;
                }
                (s, rng)
            })
            .collect();
        let mut v = Vmc { model, rbm, cfg, chains };
        if !cfg.exact {
            for c in 0..v.chains.len() {
                let (mut s, mut rng) = v.chains[c].clone();
                metropolis(&v.rbm, &mut s, &mut rng, cfg.burn * n, cfg.magnetisation.is_some(), None);
                v.chains[c] = (s, rng);
            }
        }
        v
    }

    fn batch(&mut self) -> Batch<T> {
        if self.cfg.exact { self.enumerate() } else { self.sample() }
    }

    fn enumerate(&self) -> Batch<T> {
        let n = self.model.n;
        assert!(n <= 24, "exact enumeration takes at most 24 spins");
        let configs: Vec<usize> = (0..1usize << n).filter(|&x| self.cfg.magnetisation.is_none_or(|m| n as i32 - 2 * x.count_ones() as i32 == m)).collect();
        let p = self.rbm.params();
        let k = configs.len();
        let mut logs = vec![0.0; k];
        let mut energies = vec![T::ZERO; k];
        let mut derivs = vec![T::ZERO; k * p];
        let threads = self.cfg.threads.max(1).min(k.div_ceil(64)).max(1);
        let per = k.div_ceil(threads);
        let work = |range: core::ops::Range<usize>, logs: &mut [f64], energies: &mut [T], derivs: &mut [T]| {
            let start = range.start;
            for idx in range {
                let x = configs[idx];
                let s: Vec<i8> = (0..n).map(|i| if x >> (n - 1 - i) & 1 == 0 { 1 } else { -1 }).collect();
                let cache = self.rbm.cache(&s);
                let r = idx - start;
                logs[r] = self.rbm.log_abs_cached(&s, &cache);
                if logs[r] == f64::NEG_INFINITY {
                    // The projection vanishes here: no weight, nothing to add.
                    continue;
                }
                energies[r] = local_energy(&self.model, &self.rbm, &s, &cache);
                self.rbm.derivatives_cached(&s, &cache, &mut derivs[r * p..(r + 1) * p]);
            }
        };
        if threads <= 1 || cfg!(target_arch = "wasm32") {
            work(0..k, &mut logs, &mut energies, &mut derivs);
        } else {
            std::thread::scope(|scope| {
                let mut lrest: &mut [f64] = &mut logs;
                let mut erest: &mut [T] = &mut energies;
                let mut drest: &mut [T] = &mut derivs;
                for t in 0..threads {
                    let lo = t * per;
                    let hi = ((t + 1) * per).min(k);
                    if lo >= hi {
                        break;
                    }
                    let (l, lr) = lrest.split_at_mut(hi - lo);
                    let (e, er) = erest.split_at_mut(hi - lo);
                    let (d, dr) = drest.split_at_mut((hi - lo) * p);
                    lrest = lr;
                    erest = er;
                    drest = dr;
                    let work = &work;
                    scope.spawn(move || work(lo..hi, l, e, d));
                }
            });
        }
        // |ψ|², normalised, in a fixed order.
        let top = logs.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        let raw: Vec<f64> = logs.iter().map(|l| exp(2.0 * (l - top))).collect();
        let z: f64 = raw.iter().sum();
        Batch { weights: raw.iter().map(|r| r / z).collect(), energies, derivs, acceptance: 1.0, guided: false }
    }

    fn sample(&mut self) -> Batch<T> {
        let n = self.model.n;
        let p = self.rbm.params();
        let nc = self.chains.len();
        let per_chain = self.cfg.samples.div_ceil(nc).max(1);
        let exchange = self.cfg.magnetisation.is_some();
        let thin = self.cfg.thin.max(1) * n;
        let rbm = &self.rbm;
        let model = &self.model;
        let guide = (self.cfg.guide > 0.0).then(|| {
            let logs: Vec<f64> = self.chains.iter().map(|(s, _)| rbm.log_abs_cached(s, &rbm.cache(s))).filter(|l| l.is_finite()).collect();
            Guide { c: logs.iter().sum::<f64>() / logs.len().max(1) as f64, ln_floor: ln(self.cfg.guide) }
        });
        let run = |chain: &mut (Vec<i8>, Rng)| -> (Vec<T>, Vec<T>, Vec<f64>, f64) {
            let (s, rng) = chain;
            let mut energies = Vec::with_capacity(per_chain);
            let mut derivs = vec![T::ZERO; per_chain * p];
            let mut weights = Vec::with_capacity(if guide.is_some() { per_chain } else { 0 });
            let (mut accepted, mut proposed) = (0usize, 0usize);
            for k in 0..per_chain {
                let (a, m) = metropolis(rbm, s, rng, thin, exchange, guide);
                accepted += a;
                proposed += m;
                let cache = rbm.cache(s);
                if let Some(g) = guide {
                    weights.push(g.weight(rbm.log_abs_cached(s, &cache)));
                }
                energies.push(local_energy(model, rbm, s, &cache));
                rbm.derivatives_cached(s, &cache, &mut derivs[k * p..(k + 1) * p]);
            }
            (energies, derivs, weights, accepted as f64 / proposed.max(1) as f64)
        };
        let threads = self.cfg.threads.max(1).min(nc);
        let results: Vec<(Vec<T>, Vec<T>, Vec<f64>, f64)> = if threads <= 1 || cfg!(target_arch = "wasm32") {
            self.chains.iter_mut().map(run).collect()
        } else {
            let per = nc.div_ceil(threads);
            std::thread::scope(|scope| {
                let handles: Vec<_> = self
                    .chains
                    .chunks_mut(per)
                    .map(|group| {
                        let run = &run;
                        scope.spawn(move || group.iter_mut().map(run).collect::<Vec<_>>())
                    })
                    .collect();
                handles.into_iter().flat_map(|h| h.join().expect("a sampling thread panicked")).collect()
            })
        };
        let total = per_chain * nc;
        let mut energies = Vec::with_capacity(total);
        let mut derivs = Vec::with_capacity(total * p);
        let mut raw = Vec::with_capacity(if guide.is_some() { total } else { 0 });
        let mut acc = 0.0;
        for (e, d, w, a) in results {
            energies.extend(e);
            derivs.extend(d);
            raw.extend(w);
            acc += a;
        }
        let weights = if guide.is_some() {
            let z: f64 = raw.iter().sum();
            raw.iter().map(|w| w / z).collect()
        } else {
            vec![1.0 / total as f64; total]
        };
        Batch { weights, energies, derivs, acceptance: acc / nc as f64, guided: guide.is_some() }
    }

    /// Estimate the energy without changing the machine.
    pub fn estimate(&mut self) -> Estimate {
        let b = self.batch();
        summarise(&b, self.cfg.exact)
    }

    /// One step of stochastic reconfiguration: `θ ← θ − η (S + λ)⁻¹ F`, with
    /// `S` the covariance of the log-derivatives and `F` their covariance
    /// with the local energy (Hermitian ones, for a complex machine). It is
    /// solved in parameter space or in sample space, whichever is smaller;
    /// the two give the same step.
    pub fn step(&mut self) -> Estimate {
        let b = self.batch();
        let est = summarise(&b, self.cfg.exact);
        let delta = sr_direction(&b, self.rbm.params(), self.cfg.shift);
        self.rbm.shift(&delta, self.cfg.lr);
        est
    }
}

fn summarise<T: Amplitude>(b: &Batch<T>, exact: bool) -> Estimate {
    let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e.re()).sum();
    let et = T::ONE.scale(energy);
    let variance: f64 = b
        .weights
        .iter()
        .zip(&b.energies)
        .map(|(&w, e)| {
            let d = e.sub(et);
            d.scale(w).mul(d.conj()).re()
        })
        .sum();
    // Weighted samples: the variance of a weighted mean, Σ w² |E − Ē|².
    let error = if exact {
        0.0
    } else if b.guided {
        b.weights.iter().zip(&b.energies).map(|(&w, e)| w * w * e.sub(et).norm2()).sum::<f64>().sqrt()
    } else {
        (variance / b.energies.len() as f64).sqrt()
    };
    Estimate { energy, variance, error, acceptance: b.acceptance }
}

/// `(S + λ)⁻¹ F` from a batch.
fn sr_direction<T: Amplitude>(b: &Batch<T>, p: usize, shift: f64) -> Vec<T> {
    let k = b.weights.len();
    let energy: f64 = b.weights.iter().zip(&b.energies).map(|(w, e)| w * e.re()).sum();
    let et = T::ONE.scale(energy);
    let mut mean = vec![T::ZERO; p];
    for (r, &w) in b.weights.iter().enumerate() {
        for (m, &d) in mean.iter_mut().zip(&b.derivs[r * p..(r + 1) * p]) {
            *m = m.add(d.scale(w));
        }
    }
    // Rows √w (O − Ō) and √w (E − Ē).
    let mut o = vec![T::ZERO; k * p];
    let mut eps = vec![T::ZERO; k];
    for r in 0..k {
        let sw = b.weights[r].sqrt();
        for c in 0..p {
            o[r * p + c] = b.derivs[r * p + c].sub(mean[c]).scale(sw);
        }
        eps[r] = b.energies[r].sub(et).scale(sw);
    }
    if p <= k { sr_parameter_space(&o, &eps, k, p, shift) } else { sr_sample_space(&o, eps, k, p, shift) }
}

/// `(O†O + λ)⁻¹ O†ε`.
fn sr_parameter_space<T: Amplitude>(o: &[T], eps: &[T], k: usize, p: usize, shift: f64) -> Vec<T> {
    let mut s = vec![T::ZERO; p * p];
    for r in 0..k {
        let row = &o[r * p..(r + 1) * p];
        for i in 0..p {
            let oi = row[i].conj();
            if oi.is_zero() {
                continue;
            }
            for j in i..p {
                s[i * p + j] = s[i * p + j].add(oi.mul(row[j]));
            }
        }
    }
    for i in 0..p {
        s[i * p + i] = s[i * p + i].add(T::ONE.scale(shift));
        for j in 0..i {
            s[i * p + j] = s[j * p + i].conj();
        }
    }
    let mut f = vec![T::ZERO; p];
    for r in 0..k {
        for c in 0..p {
            f[c] = f[c].add(o[r * p + c].conj().mul(eps[r]));
        }
    }
    hermitian_solve(&mut s, &mut f, p);
    f
}

/// `O† (OO† + λ)⁻¹ ε`.
fn sr_sample_space<T: Amplitude>(o: &[T], eps: Vec<T>, k: usize, p: usize, shift: f64) -> Vec<T> {
    let mut t = vec![T::ZERO; k * k];
    for i in 0..k {
        for j in i..k {
            let mut acc = T::ZERO;
            for c in 0..p {
                acc = acc.add(o[i * p + c].mul(o[j * p + c].conj()));
            }
            t[i * k + j] = acc;
            t[j * k + i] = acc.conj();
        }
        t[i * k + i] = t[i * k + i].add(T::ONE.scale(shift));
    }
    let mut x = eps;
    hermitian_solve(&mut t, &mut x, k);
    let mut delta = vec![T::ZERO; p];
    for r in 0..k {
        for c in 0..p {
            delta[c] = delta[c].add(o[r * p + c].conj().mul(x[r]));
        }
    }
    delta
}

/// Solve `A x = b` for Hermitian positive-definite `A`, `b` overwritten by
/// `x`. Cholesky first. If a pivot is lost to rounding, the metric is
/// singular at the scale of its own entries (projected networks with nearly
/// cancelling amplitudes drive it there), and the step is taken from the
/// eigendecomposition instead, dropping eigenvalues below `10⁻¹²` of the
/// largest.
fn hermitian_solve<T: Amplitude>(a: &mut [T], b: &mut [T], n: usize) {
    let (a0, b0) = (a.to_vec(), b.to_vec());
    if cholesky_solve(a, b, n) {
        return;
    }
    let (vals, vecs) = T::eigh(&a0, n);
    let top = vals.iter().fold(0.0f64, |m, v| m.max(v.abs()));
    for x in b.iter_mut() {
        *x = T::ZERO;
    }
    for (i, &lambda) in vals.iter().enumerate() {
        if lambda <= 1e-12 * top {
            continue;
        }
        // (v_i† b₀ / λ_i) v_i
        let mut c = T::ZERO;
        for r in 0..n {
            c = c.add(vecs[r * n + i].conj().mul(b0[r]));
        }
        let c = c.div_re(lambda);
        for r in 0..n {
            b[r] = b[r].add(vecs[r * n + i].mul(c));
        }
    }
}

/// `A x = b` by Cholesky (`A` overwritten by `L`, `A = L L†`; `b` by `x`).
/// Returns false, leaving both unusable, if a pivot is lost: not positive
/// or below rounding against its diagonal entry.
#[allow(clippy::needless_range_loop)]
fn cholesky_solve<T: Amplitude>(a: &mut [T], b: &mut [T], n: usize) -> bool {
    for j in 0..n {
        let diag = a[j * n + j].re();
        let mut d = diag;
        for k in 0..j {
            d -= a[j * n + k].norm2();
        }
        if !d.is_finite() || d <= 1e-14 * diag.abs() {
            return false;
        }
        let d = d.max(1e-300).sqrt();
        a[j * n + j] = T::ONE.scale(d);
        for i in j + 1..n {
            let mut s = a[i * n + j];
            for k in 0..j {
                s = s.sub(a[i * n + k].mul(a[j * n + k].conj()));
            }
            a[i * n + j] = s.div_re(d);
        }
    }
    for i in 0..n {
        let mut s = b[i];
        for k in 0..i {
            s = s.sub(a[i * n + k].mul(b[k]));
        }
        b[i] = s.div_re(a[i * n + i].re());
    }
    for i in (0..n).rev() {
        let mut s = b[i];
        for k in i + 1..n {
            s = s.sub(a[k * n + i].conj().mul(b[k]));
        }
        b[i] = s.div_re(a[i * n + i].re());
    }
    b.iter().all(|x| x.re().is_finite() && x.norm2().is_finite())
}

/// A random configuration, at the given magnetisation if any.
fn start_configuration(n: usize, magnetisation: Option<i32>, rng: &mut Rng) -> Vec<i8> {
    match magnetisation {
        None => (0..n).map(|_| if rng.next() & 1 == 0 { 1 } else { -1 }).collect(),
        Some(m) => {
            let ups = ((n as i32 + m) / 2).clamp(0, n as i32) as usize;
            let mut s: Vec<i8> = (0..n).map(|i| if i < ups { 1 } else { -1 }).collect();
            // Fisher–Yates.
            for i in (1..n).rev() {
                let j = rng.below(i + 1);
                s.swap(i, j);
            }
            s
        }
    }
}

/// `moves` Metropolis moves on `s` (single flips, or swaps of antiparallel
/// spins), accepted with probability `|ψ(s′)/ψ(s)|²`, or with a guide its
/// `q(s′)/q(s)`. Returns the number accepted and the number proposed: a draw
/// of two parallel spins proposes nothing.
fn metropolis<T: Amplitude>(rbm: &Rbm<T>, s: &mut [i8], rng: &mut Rng, moves: usize, exchange: bool, guide: Option<Guide>) -> (usize, usize) {
    let n = s.len();
    let mut cache = rbm.cache(s);
    let (mut accepted, mut proposed) = (0, 0);
    let mut l = if guide.is_some() { rbm.log_abs_cached(s, &cache) } else { 0.0 };
    for _ in 0..moves {
        let flips: Vec<usize> = if exchange {
            let i = rng.below(n);
            let j = rng.below(n);
            if s[i] == s[j] {
                continue;
            }
            vec![i, j]
        } else {
            vec![rng.below(n)]
        };
        proposed += 1;
        let r = rbm.ratio_cached(s, &cache, &flips);
        let r2 = r.norm2();
        let p = match guide {
            None => r2,
            // A ratio below 10⁻¹⁰ in size is a move to where the projection
            // vanishes, to rounding: outside the guide, as it is outside |ψ|².
            Some(_) if r2 < 1e-20 => 0.0,
            Some(g) => exp(g.ln_q(l + 0.5 * ln(r2)) - g.ln_q(l)),
        };
        if rng.uniform() < p {
            rbm.accept(s, &mut cache, &flips, r);
            accepted += 1;
            if guide.is_some() {
                l += 0.5 * ln(r2);
            }
        }
    }
    (accepted, proposed)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum_dmrg::{Chain, exact_ground_energy};

    #[test]
    fn exact_diagonalisation_matches_the_chain_referees() {
        let n = 10;
        let e = chain_edges(n, false);
        let a = SpinModel::heisenberg(n, &e, 1.0).exact_ground_energy(Some(0)).unwrap();
        let b = exact_ground_energy(&Chain::heisenberg(n, 1.0, 1.0, 0.0)).unwrap();
        assert!((a - b).abs() < 1e-9, "{a} vs {b}");
        let a = SpinModel::ising(n, &e, 1.0, 0.7).exact_ground_energy(None).unwrap();
        let b = exact_ground_energy(&Chain::ising(n, 1.0, 0.7)).unwrap();
        assert!((a - b).abs() < 1e-9, "{a} vs {b}");
    }

    #[test]
    fn the_operator_sum_is_the_model() {
        // A 3 × 3 lattice with every kind of term.
        let n = 9;
        let mut m = SpinModel::heisenberg(n, &square_edges(3, 3, false), 0.8);
        m.x.push((4, 0.3));
        m.z.push((2, -0.45));
        m.zz.push((0, 8, 0.2));
        let dense = m.opsum().mpo().to_dense().unwrap();
        let dim = 1usize << n;
        let (states, place) = m.sector_states(None);
        let mut col = vec![0.0; dim];
        let mut out = vec![0.0; dim];
        for c in 0..dim {
            col.iter_mut().for_each(|x| *x = 0.0);
            col[c] = 1.0;
            m.apply_sector(&states, &place, &col, &mut out);
            for r in 0..dim {
                assert!((out[r] - dense[r * dim + c]).abs() < 1e-12, "({r}, {c})");
            }
        }
    }

    #[test]
    fn ratios_and_derivatives_are_the_machines() {
        let n = 7;
        let mut rbm = Rbm::new(n, 2, 0.4, 3).with_marshall_sign(checkerboard(n, 1));
        let mut rng = Rng(5);
        for v in rbm.a.iter_mut().chain(rbm.b.iter_mut()) {
            *v = 0.3 * (rng.uniform() - 0.5);
        }
        for _ in 0..20 {
            let s = start_configuration(n, None, &mut rng);
            let theta = rbm.theta(&s);
            for flips in [vec![rng.below(n)], vec![0, 3], vec![2, 6]] {
                let mut t = s.clone();
                for &i in &flips {
                    t[i] = -t[i];
                }
                let want = rbm.sign(&t) * rbm.sign(&s) * exp(rbm.log_amplitude(&t) - rbm.log_amplitude(&s));
                let got = rbm.ratio(&s, &theta, &flips);
                assert!((got - want).abs() < 1e-12 * want.abs().max(1.0), "{got} vs {want}");
            }
            // Derivatives by central differences.
            let p = rbm.params();
            let mut d = vec![0.0; p];
            rbm.log_derivatives(&s, &theta, &mut d);
            let h = 1e-6;
            for k in [0, n, n + 3, n + rbm.m + 5, p - 1] {
                let mut e = vec![0.0; p];
                e[k] = -h;
                let mut up = rbm.clone();
                up.shift(&e, 1.0);
                e[k] = h;
                let mut dn = rbm.clone();
                dn.shift(&e, 1.0);
                let fd = (up.log_amplitude(&s) - dn.log_amplitude(&s)) / (2.0 * h);
                assert!((fd - d[k]).abs() < 1e-7, "param {k}: {fd} vs {}", d[k]);
            }
        }
    }

    #[test]
    fn both_spaces_take_the_same_step() {
        let mut rng = Rng(9);
        for (k, p) in [(12usize, 5usize), (5, 12)] {
            let o: Vec<f64> = (0..k * p).map(|_| rng.uniform() - 0.5).collect();
            let eps: Vec<f64> = (0..k).map(|_| rng.uniform() - 0.5).collect();
            let a = sr_parameter_space(&o, &eps, k, p, 1e-3);
            let b = sr_sample_space(&o, eps.clone(), k, p, 1e-3);
            for (x, y) in a.iter().zip(&b) {
                assert!((x - y).abs() < 1e-9 * x.abs().max(1.0), "{x} vs {y}");
            }
        }
    }

    #[test]
    fn exact_vmc_finds_the_ground_state() {
        let n = 8;
        let model = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0);
        let exact = model.exact_ground_energy(Some(0)).unwrap();
        let rbm = Rbm::new(n, 2, 0.05, 1).with_marshall_sign(checkerboard(n, 1));
        let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-4, ..VmcConfig::default() };
        let mut vmc = Vmc::new(model, rbm, cfg);
        let mut last = vmc.step();
        for _ in 0..300 {
            let e = vmc.step();
            last = e;
        }
        let rel = (last.energy - exact) / exact.abs();
        assert!(last.energy >= exact - 1e-9, "variational: {} vs {exact}", last.energy);
        assert!(rel < 1e-3, "{} vs {exact} (relative {rel:e})", last.energy);
        assert!(last.variance < 1e-2, "{}", last.variance);
    }

    #[test]
    fn sampling_agrees_with_enumeration_and_repeats() {
        let n = 8;
        let model = SpinModel::ising(n, &chain_edges(n, false), 1.0, 1.0);
        let rbm = Rbm::new(n, 1, 0.3, 4);
        let exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, ..VmcConfig::default() }).estimate();
        let cfg = VmcConfig { samples: 8000, chains: 8, threads: 1, ..VmcConfig::default() };
        let sampled = Vmc::new(model.clone(), rbm.clone(), cfg).estimate();
        assert!((sampled.energy - exact.energy).abs() < 5.0 * sampled.error, "{sampled:?} vs {exact:?}");
        // The same bits at any thread count, and again.
        let mut one = Vmc::new(model.clone(), rbm.clone(), VmcConfig { samples: 400, threads: 1, ..VmcConfig::default() });
        let mut many = Vmc::new(model, rbm, VmcConfig { samples: 400, threads: 3, ..VmcConfig::default() });
        for _ in 0..3 {
            assert_eq!(one.step(), many.step());
        }
        assert_eq!(one.rbm, many.rbm);
        assert_eq!(one.chains, many.chains);
    }

    #[test]
    fn a_complex_machine_with_real_parameters_is_the_real_one() {
        let n = 6;
        let real = Rbm::new(n, 2, 0.4, 5);
        let cplx = Rbm { n: real.n, m: real.m, a: real.a.iter().map(|&re| C { re, im: 0.0 }).collect(), b: real.b.iter().map(|&re| C { re, im: 0.0 }).collect(), w: real.w.iter().map(|&re| C { re, im: 0.0 }).collect(), sublattice: None, symmetry: None };
        let mut rng = Rng(8);
        for _ in 0..10 {
            let s = start_configuration(n, None, &mut rng);
            let (tr, tc) = (real.theta(&s), cplx.theta(&s));
            assert!((real.log_amplitude(&s) - cplx.log_amplitude(&s)).abs() < 1e-14);
            assert!(cplx.log_psi(&s).im.abs() < 1e-14);
            for flips in [vec![1], vec![0, 4]] {
                let (r, c) = (real.ratio(&s, &tr, &flips), cplx.ratio(&s, &tc, &flips));
                assert!((r - c.re).abs() < 1e-13 * r.abs().max(1.0) && c.im.abs() < 1e-13);
            }
        }
    }

    #[test]
    fn complex_ratios_and_derivatives_are_the_machines() {
        let n = 6;
        let mut rbm = Rbm::complex(n, 2, 0.4, 6);
        let mut rng = Rng(12);
        for v in rbm.a.iter_mut().chain(rbm.b.iter_mut()) {
            *v = C { re: 0.3 * (rng.uniform() - 0.5), im: 0.3 * (rng.uniform() - 0.5) };
        }
        let p = rbm.params();
        for _ in 0..10 {
            let s = start_configuration(n, None, &mut rng);
            let theta = rbm.theta(&s);
            for flips in [vec![rng.below(n)], vec![1, 5]] {
                let mut t = s.clone();
                for &i in &flips {
                    t[i] = -t[i];
                }
                let want = rbm.log_psi(&t).sub(rbm.log_psi(&s)).exp();
                let got = rbm.ratio(&s, &theta, &flips);
                assert!(got.sub(want).norm2().sqrt() < 1e-12 * want.norm2().sqrt().max(1.0), "{got:?} vs {want:?}");
            }
            // Holomorphic: the derivative along a real step is the derivative.
            let mut d = vec![C::ZERO; p];
            rbm.log_derivatives(&s, &theta, &mut d);
            let h = 1e-6;
            for k in [0, n + 1, n + rbm.m + 3, p - 1] {
                let mut e = vec![C::ZERO; p];
                e[k] = C { re: -h, im: 0.0 };
                let mut up = rbm.clone();
                up.shift(&e, 1.0);
                e[k] = C { re: h, im: 0.0 };
                let mut dn = rbm.clone();
                dn.shift(&e, 1.0);
                let fd = up.log_psi(&s).sub(dn.log_psi(&s)).scale(0.5 / h);
                assert!(fd.sub(d[k]).norm2().sqrt() < 1e-7, "param {k}: {fd:?} vs {:?}", d[k]);
            }
        }
    }

    #[test]
    fn both_spaces_take_the_same_complex_step() {
        let mut rng = Rng(19);
        for (k, p) in [(12usize, 5usize), (5, 12)] {
            let o: Vec<C> = (0..k * p).map(|_| C { re: rng.uniform() - 0.5, im: rng.uniform() - 0.5 }).collect();
            let eps: Vec<C> = (0..k).map(|_| C { re: rng.uniform() - 0.5, im: rng.uniform() - 0.5 }).collect();
            let a = sr_parameter_space(&o, &eps, k, p, 1e-3);
            let b = sr_sample_space(&o, eps.clone(), k, p, 1e-3);
            for (x, y) in a.iter().zip(&b) {
                assert!(x.sub(*y).norm2().sqrt() < 1e-9 * x.norm2().sqrt().max(1.0), "{x:?} vs {y:?}");
            }
        }
    }

    #[test]
    fn exact_diagonalisation_solves_majumdar_ghosh() {
        // At J₂ = J₁/2 a ring's ground states are the two dimer coverings,
        // at exactly −3n/8.
        for n in [8usize, 12] {
            let m = SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5);
            let e = m.exact_ground_energy(Some(0)).unwrap();
            assert!((e + 3.0 * n as f64 / 8.0).abs() < 1e-9, "n={n}: {e}");
        }
    }

    #[test]
    fn a_complex_machine_learns_a_sign_a_real_one_cannot() {
        // The Heisenberg ring and the Majumdar–Ghosh ring, with no sign rule
        // given. A real network is positive and stays far off; a complex one
        // finds the signs.
        let n = 8;
        let models = [
            SpinModel::heisenberg(n, &chain_edges(n, true), 1.0),
            SpinModel::j1j2(n, &chain_edges(n, true), &chain_next_nearest(n, true), 1.0, 0.5),
        ];
        for model in models {
            let exact = model.exact_ground_energy(Some(0)).unwrap();
            let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-3, ..VmcConfig::default() };
            let mut real = Vmc::new(model.clone(), Rbm::new(n, 2, 0.05, 1), cfg);
            let mut cplx = Vmc::new(model, Rbm::complex(n, 2, 0.05, 3), cfg);
            let (mut er, mut ec) = (real.step(), cplx.step());
            for _ in 0..250 {
                er = real.step();
                ec = cplx.step();
            }
            assert!((er.energy - exact) / exact.abs() > 0.1, "real: {} vs {exact}", er.energy);
            assert!(ec.energy >= exact - 1e-9, "variational: {} vs {exact}", ec.energy);
            assert!((ec.energy - exact) / exact.abs() < 2e-3, "complex: {} vs {exact}", ec.energy);
        }
    }

    /// `Ψ(s)` by summing over the group.
    fn projected<T: Amplitude>(rbm: &Rbm<T>, s: &[i8]) -> T {
        let sym = rbm.symmetry.as_ref().unwrap();
        let mut total = T::ZERO;
        for g in 0..sym.len() {
            let t: Vec<i8> = sym.perms[g].iter().map(|&j| if sym.flips[g] { -s[j] } else { s[j] }).collect();
            total = total.add(rbm.log_psi(&t).exp().scale(sym.chars[g] * rbm.sign(&t)));
        }
        total
    }

    #[test]
    fn the_projection_is_the_group_sum() {
        let n = 6;
        for odd in [false, true] {
            let rbm = Rbm::complex(n, 1, 0.3, 4).with_symmetry(Symmetry::ring(n, true).with_spin_flip(odd));
            let mut rng = Rng(3);
            for _ in 0..12 {
                let s = start_configuration(n, Some(0), &mut rng);
                let base = projected(&rbm, &s);
                if base.norm2() < 1e-20 {
                    continue;
                }
                let cache = rbm.cache(&s);
                for flips in [vec![0usize, 1], vec![2, 5], vec![3]] {
                    let mut t = s.clone();
                    for &i in &flips {
                        t[i] = -t[i];
                    }
                    let want = projected(&rbm, &t).div(base);
                    let got = rbm.ratio_cached(&s, &cache, &flips);
                    assert!(got.sub(want).norm2().sqrt() < 1e-10 * want.norm2().sqrt().max(1.0), "{got:?} vs {want:?}");
                }
                // The sector: a translation by one site multiplies by −1
                // (momentum π), a global flip by the parity.
                let shifted: Vec<i8> = (0..n).map(|i| s[(i + 1) % n]).collect();
                let flipped: Vec<i8> = s.iter().map(|&x| -x).collect();
                assert!(projected(&rbm, &shifted).add(base).norm2().sqrt() < 1e-10 * base.norm2().sqrt());
                let sign = if odd { -1.0 } else { 1.0 };
                assert!(projected(&rbm, &flipped).sub(base.scale(sign)).norm2().sqrt() < 1e-10 * base.norm2().sqrt());
                // Derivatives of ln Ψ by central differences along a real step.
                let p = rbm.params();
                let mut d = vec![C::ZERO; p];
                rbm.derivatives_cached(&s, &cache, &mut d);
                let h = 1e-6;
                for k in [0, n + 1, p - 1] {
                    let mut e = vec![C::ZERO; p];
                    e[k] = C { re: -h, im: 0.0 };
                    let mut up = rbm.clone();
                    up.shift(&e, 1.0);
                    e[k] = C { re: h, im: 0.0 };
                    let mut dn = rbm.clone();
                    dn.shift(&e, 1.0);
                    let fd = projected(&up, &s).sub(projected(&dn, &s)).div(base).scale(0.5 / h);
                    assert!(fd.sub(d[k]).norm2().sqrt() < 1e-6, "param {k}: {fd:?} vs {:?}", d[k]);
                }
            }
        }
    }

    #[test]
    fn accepted_moves_keep_the_cache_exact() {
        let n = 8;
        let rbm = Rbm::complex(n, 1, 0.3, 9).with_symmetry(Symmetry::ring(n, false).with_spin_flip(false));
        let mut rng = Rng(21);
        let mut s = start_configuration(n, Some(0), &mut rng);
        let mut cache = rbm.cache(&s);
        for _ in 0..200 {
            let (i, j) = (rng.below(n), rng.below(n));
            if s[i] == s[j] {
                continue;
            }
            let r = rbm.ratio_cached(&s, &cache, &[i, j]);
            rbm.accept(&mut s, &mut cache, &[i, j], r);
        }
        let fresh = rbm.cache(&s);
        match (&cache, &fresh) {
            (Cache::Projected { configs: c1, shares: p1, log_abs: l1, thetas: t1 }, Cache::Projected { configs: c2, shares: p2, log_abs: l2, thetas: t2 }) => {
                assert_eq!(c1, c2);
                assert!((l1 - l2).abs() < 1e-9, "{l1} vs {l2}");
                for (a, b) in p1.iter().zip(p2) {
                    assert!(a.sub(*b).norm2().sqrt() < 1e-9, "{a:?} vs {b:?}");
                }
                for (a, b) in t1.iter().flatten().zip(t2.iter().flatten()) {
                    assert!(a.sub(*b).norm2().sqrt() < 1e-12);
                }
            }
            _ => unreachable!(),
        }
    }

    #[test]
    fn projection_reaches_the_ground_state_in_its_sector() {
        // A 6-site ring: half its length is odd, so the singlet ground state
        // has momentum π and odd parity under a global flip.
        let n = 6;
        let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0);
        let exact = model.exact_ground_energy(Some(0)).unwrap();
        let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.1, shift: 1e-3, ..VmcConfig::default() };
        let rbm = Rbm::complex(n, 1, 0.05, 3).with_symmetry(Symmetry::ring(n, true).with_spin_flip(true));
        let mut v = Vmc::new(model, rbm, cfg);
        let mut e = v.step();
        for _ in 0..150 {
            e = v.step();
        }
        assert!(e.energy >= exact - 1e-9, "variational: {} vs {exact}", e.energy);
        assert!((e.energy - exact) / exact.abs() < 1e-8, "{} vs {exact}", e.energy);
    }

    #[test]
    fn square_groups_are_groups_with_characters() {
        for l in [3, 4] {
            for momentum_pi in [false, true] {
                if momentum_pi && l % 2 == 1 {
                    continue;
                }
                for irrep in [PointIrrep::A1, PointIrrep::A2, PointIrrep::B1, PointIrrep::B2] {
                    let sym = Symmetry::square(l, momentum_pi, irrep);
                    assert_eq!(sym.len(), 8 * l * l);
                    assert_eq!(sym.check(), Ok(()), "{l} {momentum_pi} {irrep:?}");
                }
                for odd in [false, true] {
                    assert_eq!(Symmetry::square(l, momentum_pi, PointIrrep::B2).with_spin_flip(odd).check(), Ok(()));
                }
            }
        }
        for (lx, ly) in [(4, 3), (6, 4)] {
            for (kx, ky) in [(false, false), (true, false), (false, true), (true, true)] {
                if (kx && lx % 2 == 1) || (ky && ly % 2 == 1) {
                    continue;
                }
                for irrep in [RectIrrep::A1, RectIrrep::A2, RectIrrep::B1, RectIrrep::B2] {
                    let sym = Symmetry::rectangle(lx, ly, kx, ky, irrep);
                    assert_eq!(sym.len(), 4 * lx * ly);
                    assert_eq!(sym.check(), Ok(()), "{lx}×{ly} {kx} {ky} {irrep:?}");
                }
            }
        }
        // Characters that do not multiply, and elements that do not close.
        let ring = Symmetry::ring(4, false);
        assert!(Symmetry::new(ring.perms.clone(), vec![1.0, -1.0, 1.0, 1.0]).check().is_err());
        assert!(Symmetry::new(ring.perms[..2].to_vec(), vec![1.0, 1.0]).check().is_err());
        assert!(Symmetry::new(ring.perms.clone(), vec![1.0, 0.5, 1.0, 0.5]).check().is_err());
    }

    #[test]
    fn a_square_projection_carries_its_representation() {
        let (l, n) = (4, 16);
        let quarter_turn = |x: usize, y: usize| ((l - y) % l, x);
        let axis = |x: usize, y: usize| ((l - x) % l, y);
        let diagonal = |x: usize, y: usize| (y, x);
        let step = |x: usize, y: usize| ((x + 1) % l, y);
        for irrep in [PointIrrep::A1, PointIrrep::A2, PointIrrep::B1, PointIrrep::B2] {
            let rbm = Rbm::complex(n, 1, 0.3, 11).with_symmetry(Symmetry::square(l, true, irrep).with_spin_flip(false));
            let mut rng = Rng(5);
            let mut seen = 0;
            // A configuration some element with character −1 leaves alone
            // has no amplitude in the sector; on sixteen sites that is up to
            // two in five of them.
            for _ in 0..40 {
                if seen == 6 {
                    break;
                }
                let s = start_configuration(n, Some(0), &mut rng);
                let base = projected(&rbm, &s);
                if base.norm2() < 1e-16 {
                    continue;
                }
                seen += 1;
                let moved = |f: &dyn Fn(usize, usize) -> (usize, usize)| -> Vec<i8> {
                    (0..n)
                        .map(|i| {
                            let (x, y) = f(i % l, i / l);
                            s[y * l + x]
                        })
                        .collect()
                };
                // Momentum (π, π): one step along an axis multiplies by −1.
                let want = [(irrep.character(1, false), moved(&quarter_turn)), (irrep.character(0, true), moved(&axis)), (irrep.character(1, true), moved(&diagonal)), (-1.0, moved(&step))];
                for (chi, t) in want {
                    let got = projected(&rbm, &t);
                    assert!(got.sub(base.scale(chi)).norm2().sqrt() < 1e-10 * base.norm2().sqrt(), "{irrep:?}: {got:?} vs {chi} × {base:?}");
                }
            }
            assert_eq!(seen, 6, "{irrep:?}: too few configurations survive the projection");
        }
    }

    #[test]
    fn a_projected_start_near_zero_sits_on_the_ferromagnet() {
        // An 8-site ring projected onto momentum 0 and even parity, the
        // ground state's sector. Parameters within ±0.01 make the network
        // nearly the projection of a constant: the uniform superposition,
        // which is the ferromagnet's member at zero magnetisation, energy
        // n/4, an eigenstate. Reconfiguration barely moves it.
        let n = 8;
        let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0);
        let exact = model.exact_ground_energy(Some(0)).unwrap();
        let ferromagnet = n as f64 / 4.0;
        let cfg = VmcConfig { exact: true, magnetisation: Some(0), lr: 0.02, shift: 1e-3, ..VmcConfig::default() };
        let run = |init: f64, steps: usize| {
            let rbm = Rbm::complex(n, 1, init, 7).with_symmetry(Symmetry::ring(n, false).with_spin_flip(false));
            let mut v = Vmc::new(model.clone(), rbm, cfg);
            let first = v.estimate();
            for _ in 0..steps {
                v.step();
            }
            (first, v.estimate())
        };
        let (first, after) = run(0.01, 40);
        assert!((first.energy - ferromagnet).abs() < 1e-3 && first.variance < 1e-3, "{first:?}");
        assert!((after.energy - ferromagnet).abs() < 1e-2, "{after:?}");
        let (_, after) = run(0.1, 200);
        assert!((after.energy - exact) / exact.abs() < 1e-3, "{} vs {exact}", after.energy);
    }

    #[test]
    fn guided_sampling_agrees_with_enumeration() {
        // A projected network on an 8-site ring, momentum 0 and even parity.
        let n = 8;
        let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0);
        let rbm = Rbm::complex(n, 1, 0.3, 5).with_symmetry(Symmetry::ring(n, false).with_spin_flip(false));
        let exact = Vmc::new(model.clone(), rbm.clone(), VmcConfig { exact: true, magnetisation: Some(0), ..VmcConfig::default() }).estimate();
        let cfg = VmcConfig { samples: 4000, chains: 8, magnetisation: Some(0), guide: 0.5, threads: 1, ..VmcConfig::default() };
        let g = Vmc::new(model.clone(), rbm.clone(), cfg).estimate();
        assert!((g.energy - exact.energy).abs() < 5.0 * g.error, "{g:?} vs {exact:?}");
        assert!((g.variance - exact.variance).abs() < 0.2 * exact.variance, "{g:?} vs {exact:?}");
        // The same bits at any thread count.
        let cfg = |threads| VmcConfig { samples: 200, magnetisation: Some(0), guide: 0.5, threads, ..VmcConfig::default() };
        let mut one = Vmc::new(model.clone(), rbm.clone(), cfg(1));
        let mut many = Vmc::new(model, rbm, cfg(3));
        for _ in 0..2 {
            assert_eq!(one.step(), many.step());
        }
        assert_eq!(one.rbm, many.rbm);
    }

    #[test]
    fn no_chain_starts_where_the_projection_vanishes() {
        // At momentum (π, π), even under every rotation and reflection, about
        // a quarter of the 4 × 4 lattice's configurations are fixed by an
        // element of character −1 and have no amplitude. A chain there has
        // no move with a ratio and would never leave.
        let (l, n) = (4, 16);
        let model = SpinModel::heisenberg(n, &square_edges(l, l, true), 1.0);
        let rbm = Rbm::complex(n, 1, 0.1, 3).with_symmetry(Symmetry::square(l, true, PointIrrep::A1).with_spin_flip(false));
        let cfg = VmcConfig { chains: 32, burn: 1, magnetisation: Some(0), threads: 1, ..VmcConfig::default() };
        let v = Vmc::new(model, rbm, cfg);
        for (c, (s, _)) in v.chains.iter().enumerate() {
            assert!(v.rbm.log_abs_cached(s, &v.rbm.cache(s)).is_finite(), "chain {c} sits where the projection vanishes: {s:?}");
        }
    }

    #[test]
    fn sector_diagonalisation_splits_the_spectrum() {
        // A ring of ten: half its length is odd, so the singlet ground state
        // has momentum π and odd parity under the flip, and the sector of
        // momentum 0 and even parity starts higher.
        let n = 10;
        let model = SpinModel::heisenberg(n, &chain_edges(n, true), 1.0);
        let all = model.exact_ground_energy(Some(0)).unwrap();
        let ground = model.exact_ground_energy_in(Some(0), &Symmetry::ring(n, true).with_spin_flip(true)).unwrap();
        let other = model.exact_ground_energy_in(Some(0), &Symmetry::ring(n, false).with_spin_flip(false)).unwrap();
        assert!((ground - all).abs() < 1e-9, "{ground} vs {all}");
        assert!(other > all + 0.1, "{other} vs {all}");
        // Refusals: a model the group does not preserve, a flip out of the
        // magnetisation, and a field the flip reverses.
        let open = SpinModel::heisenberg(n, &chain_edges(n, false), 1.0);
        assert_eq!(open.exact_ground_energy_in(Some(0), &Symmetry::ring(n, false)), None);
        assert_eq!(model.exact_ground_energy_in(Some(2), &Symmetry::ring(n, false).with_spin_flip(false)), None);
        let mut field = model.clone();
        field.z = (0..n).map(|i| (i, 0.1)).collect();
        assert!(field.commutes_with(&Symmetry::ring(n, false)));
        assert!(!field.commutes_with(&Symmetry::ring(n, false).with_spin_flip(false)));
    }

    #[test]
    fn the_orbit_basis_agrees_with_the_projector() {
        let close = |a: Option<f64>, b: Option<f64>, what: &str| match (a, b) {
            (Some(x), Some(y)) => assert!((x - y).abs() < 1e-9 * y.abs().max(1.0), "{what}: {x} vs {y}"),
            (x, y) => assert_eq!(x.is_some(), y.is_some(), "{what}: {x:?} vs {y:?}"),
        };
        // Every sector of a 4 × 3 torus at J₂ = J₁/2: four representations,
        // both momenta along x, and no flip, an even one or an odd one.
        let (lx, ly, n) = (4, 3, 12);
        let model = SpinModel::j1j2(n, &square_edges(lx, ly, true), &square_diagonals(lx, ly, true), 1.0, 0.5);
        for kx in [false, true] {
            for irrep in [RectIrrep::A1, RectIrrep::A2, RectIrrep::B1, RectIrrep::B2] {
                let base = Symmetry::rectangle(lx, ly, kx, false, irrep);
                for sym in [base.clone(), base.clone().with_spin_flip(false), base.with_spin_flip(true)] {
                    close(model.exact_ground_energy_in(Some(0), &sym), model.projected_ground_energy(Some(0), &sym), &format!("4×3 {kx} {irrep:?} {}", sym.len()));
                }
            }
        }
        // A 10-site ring, and an Ising ring in a transverse field, where the
        // magnetisation is not conserved.
        let ring = SpinModel::heisenberg(10, &chain_edges(10, true), 1.0);
        let ising = SpinModel::ising(8, &chain_edges(8, true), 1.0, 0.7);
        for k in [false, true] {
            for odd in [false, true] {
                let sym = Symmetry::ring(10, k).with_spin_flip(odd);
                close(ring.exact_ground_energy_in(Some(0), &sym), ring.projected_ground_energy(Some(0), &sym), "ring");
                let sym = Symmetry::ring(8, k).with_spin_flip(odd);
                close(ising.exact_ground_energy_in(None, &sym), ising.projected_ground_energy(None, &sym), "ising");
            }
        }
        // The identity alone is no reduction, and the flip alone splits the
        // spectrum in two, one half holding the ground state.
        let identity = Symmetry::new(vec![(0..n).collect()], vec![1.0]);
        let all = model.exact_ground_energy(Some(0)).unwrap();
        close(model.exact_ground_energy_in(Some(0), &identity), Some(all), "identity");
        let even = model.exact_ground_energy_in(Some(0), &identity.clone().with_spin_flip(false)).unwrap();
        let odd = model.exact_ground_energy_in(Some(0), &identity.with_spin_flip(true)).unwrap();
        assert!((even.min(odd) - all).abs() < 1e-9 && even.max(odd) > all + 1e-6, "{even} {odd} {all}");
    }

    #[test]
    fn the_frustrated_square_ground_state_is_fully_symmetric() {
        // 4 × 4 at J₂ = J₁/2: the ground state is even under every rotation
        // and reflection, at momentum zero.
        let (l, n) = (4, 16);
        let model = SpinModel::j1j2(n, &square_edges(l, l, true), &square_diagonals(l, l, true), 1.0, 0.5);
        let all = model.exact_ground_energy(Some(0)).unwrap();
        let sector = |momentum_pi, irrep| model.exact_ground_energy_in(Some(0), &Symmetry::square(l, momentum_pi, irrep)).unwrap();
        let a1 = sector(false, PointIrrep::A1);
        assert!((a1 - all).abs() < 1e-9, "{a1} vs {all}");
        let b1 = sector(false, PointIrrep::B1);
        assert!(b1 > all + 0.1, "{b1} vs {all}");
    }
}