"""Exact fixed-magnetisation sector propagation for U(1)-symmetric spin models.

​

simulate(instance, times) builds the Hamiltonian from the instance's Pauli

terms restricted to the total-Z sector of the initial product state and

propagates it with a Chebyshev expansion. No offline training, no reference

data, no lookup. When the sector is too large to store (e.g. the 6 x 6

scaling instance), it falls back to overlapping exact cluster patches.

"""

import itertools

import math

​

import numpy as np

from scipy.sparse import csr_matrix, identity

from scipy.special import jv

​

PRIORS = {

"family": "exact-sector-chebyshev",

"parameters": {

"max_sector_dim": 1000000,

"propagator": "Chebyshev expansion of exp(-iH dt) between requested times, Gershgorin spectral bounds, Bessel tail < 1e-16",

"fallback": "overlapping exact clusters (up to 4 x 4) with total-Z drift correction",

"fallback_cluster": [4, 4],

"offline_training": "none",

},

"notes": (

"Adaptive/hybrid: if the fixed total-Z sector of the initial product "

"state has dimension <= max_sector_dim, the dynamics are computed "

"exactly in that sector (no approximation beyond floating point and "

"the Chebyshev truncation, Bessel coefficients below 1e-16). Otherwise (6 x 6 "

"scaling instance) each site's trace is taken from an exactly "

"simulated open-boundary cluster of at most 4 x 4 sites in which it is "

"most interior; this ignores correlations beyond the cluster and is "

"only a sanity-level approximation. The fallback subtracts the mean "

"total-Z drift per time and clips to [-1, 1]. No offline training, "

"preprocessing or use of reference traces."

),

}

​

​

def _parse_terms(instance):

terms = []

for coef, ops in instance["terms"]:

ops = {int(k): str(v).upper() for k, v in ops.items()}

c = complex(coef)

if c != 0:

terms.append((c, ops))

return terms

​

​

def _initial_bits(instance):

n = int(instance["n_sites"])

st = instance.get("initial_state", {})

bits = st.get("bits")

if bits is None:

typ = str(st.get("type", "neel")).lower()

geo = instance.get("geometry", {})

if typ == "neel" and "Ly" in geo:

ly = int(geo["Ly"])

bits = [((i // ly) + (i % ly)) % 2 for i in range(n)]

elif typ == "neel":

bits = [i % 2 for i in range(n)]

else:

bits = [0] * n

return [int(b) for b in bits]

​

​

def _sector_states(n, n_up_bits):

"""All n-bit integers with popcount n_up_bits (bit i = site i, 1 = down)."""

m = math.comb(n, n_up_bits)

states = np.empty(m, dtype=np.int64)

for k, combo in enumerate(itertools.combinations(range(n), n_up_bits)):

v = 0

for b in combo:

v |= 1 << b

states[k] = v

return np.sort(states)

​

​

def _sector_states_fast(n, k):

if math.comb(n, k) < 200000:

return _sector_states(n, k)

# Gosper-free vectorised construction via popcount filtering in blocks.

out = []

block = 1 << 22

total = 1 << n

for start in range(0, total, block):

x = np.arange(start, min(start + block, total), dtype=np.int64)

pc = np.zeros_like(x)

y = x.copy()

while True:

nz = y != 0

if not nz.any():

break

pc += y & 1

y >>= 1

out.append(x[pc == k])

return np.concatenate(out)

​

​

def _build_sector_hamiltonian(n, terms, states):

dim = states.size

index = {}

use_dict = n > 26

if not use_dict:

lookup = np.full(1 << n, -1, dtype=np.int64)

lookup[states] = np.arange(dim, dtype=np.int64)

else:

index = {int(s): i for i, s in enumerate(states)}

rows, cols, vals = [], [], []

diag = np.zeros(dim, dtype=np.complex128)

cols_all = np.arange(dim, dtype=np.int64)

for c, ops in terms:

flip = 0

amp = np.full(dim, c, dtype=np.complex128)

for site, p in ops.items():

b = (states >> site) & 1

sign = 1 - 2 * b

if p == "Z":

amp *= sign

elif p == "X":

flip |= 1 << site

elif p == "Y":

flip |= 1 << site

amp *= 1j * sign

elif p != "I":

raise ValueError("unknown Pauli " + p)

if flip == 0:

diag += amp

continue

tgt = states ^ flip

if not use_dict:

ti = lookup[tgt]

else:

ti = np.array([index.get(int(t), -1) for t in tgt], dtype=np.int64)

ok = ti >= 0

rows.append(ti[ok])

cols.append(cols_all[ok])

vals.append(amp[ok])

rows.append(cols_all)

cols.append(cols_all)

vals.append(diag)

r = np.concatenate(rows)

cc = np.concatenate(cols)

v = np.concatenate(vals)

H = csr_matrix((v, (r, cc)), shape=(dim, dim))

H.sum_duplicates()

H.eliminate_zeros()

if np.abs(H.data.imag).max(initial=0.0) < 1e-14:

H = csr_matrix((H.data.real, H.indices, H.indptr), shape=H.shape)

return H

​

​

def _is_u1(terms):

"""True if every term combination conserves total Z (checked pairwise XX/YY)."""

pend = {}

for c, ops in terms:

letters = set(ops.values()) - {"I", "Z"}

if not letters:

continue

flips = tuple(sorted(s for s, p in ops.items() if p in "XY"))

if len(flips) != 2:

return False

key = (flips, tuple(sorted((s, p) for s, p in ops.items() if p == "Z")))

a, b = flips

pa, pb = ops[a], ops[b]

pend.setdefault(key, {})

pend[key][(pa, pb)] = pend[key].get((pa, pb), 0) + c

for d in pend.values():

if abs(d.get(("X", "X"), 0) - d.get(("Y", "Y"), 0)) > 1e-12:

return False

if abs(d.get(("X", "Y"), 0) + d.get(("Y", "X"), 0)) > 1e-12:

return False

return True

​

​

def _evolve_sector(n, terms, bits, times):

k = int(sum(bits))

states = _sector_states_fast(n, k)

H = _build_sector_hamiltonian(n, terms, states)

init = 0

for i, b in enumerate(bits):

if b:

init |= 1 << i

psi0 = np.zeros(states.size, dtype=np.complex128)

psi0[int(np.searchsorted(states, init))] = 1.0

zsign = np.empty((n, states.size), dtype=np.float64)

for i in range(n):

zsign[i] = 1.0 - 2.0 * ((states >> i) & 1)

times = np.asarray(times, dtype=float)

out = np.empty((times.size, n))

prop = _ChebyshevPropagator(H)

order = np.argsort(times, kind="stable")

psi = np.column_stack([psi0.real, psi0.imag])

tprev = 0.0

for idx in order:

t = float(times[idx])

if t != tprev:

psi = prop.step(psi, t - tprev)

tprev = t

out[idx] = zsign @ (psi[:, 0] ** 2 + psi[:, 1] ** 2)