def make_instance(family, size, seed):

import numpy as np

rng = np.random.default_rng(seed)

ham = family["hamiltonian"]; kind = ham["type"]

terms = []

if kind == "xxz_chain":

L = size; J, delta = ham["J"], ham["delta"]

for i in range(L - 1):

terms += [[J, {str(i): "X", str(i + 1): "X"}], [J, {str(i): "Y", str(i + 1): "Y"}], [J * delta, {str(i): "Z", str(i + 1): "Z"}]]

W = ham["field"]["W"]

for i in range(L): terms.append([float(rng.uniform(-W, W)), {str(i): "Z"}])

n_sites, geometry = L, {"type": "chain", "L": L}

elif kind == "tfim_chain_longitudinal":

L = size; J, hx = ham["J"], ham["hx"]; W = ham["hz"]["W"]

for i in range(L - 1): terms.append([-J, {str(i): "Z", str(i + 1): "Z"}])

for i in range(L): terms += [[-hx, {str(i): "X"}], [float(rng.uniform(-W, W)), {str(i): "Z"}]]

n_sites, geometry = L, {"type": "chain", "L": L}

elif kind == "heisenberg_2d":

Lx, Ly = size; J = ham["J"]; W = ham["field"]["W"]

idx = lambda x, y: x * Ly + y

for x in range(Lx):

for y in range(Ly):

for (x2, y2) in ((x + 1, y), (x, y + 1)):

if x2 < Lx and y2 < Ly:

for p in "XYZ": terms.append([J, {str(idx(x, y)): p, str(idx(x2, y2)): p}])

for s in range(Lx * Ly): terms.append([float(rng.uniform(-W, W)), {str(s): "Z"}])

n_sites, geometry = Lx * Ly, {"type": "square_open", "Lx": Lx, "Ly": Ly}

else:

raise ValueError(f"unknown hamiltonian type {kind}")

init = family["initial_state"]

if init["type"] == "random_product_z":

bits = (rng.uniform(0, 1, n_sites) >= init["p_up"]).astype(int).tolist()

elif init["type"] == "neel":

bits = [(x+y) % 2 for x in range(size[0]) for y in range(size[1])] if kind == "heisenberg_2d" else [i % 2 for i in range(n_sites)]

elif init["type"] == "domain_wall":

bits = [0] * (n_sites // 2) + [1] * (n_sites - n_sites // 2)

else:

raise ValueError("unknown initial state")

times = [float(x) for x in np.linspace(0.0, family["times"]["T"], family["times"]["num"])]

obs = [{"site": i, "pauli": family["observable"]["pauli"]} for i in range(n_sites)]

return {"family": family["family"], "n_sites": n_sites, "geometry": geometry, "terms": terms,

"initial_state": {"type": "product_z", "bits": bits}, "observables": obs, "times": times,

"conserved": family.get("conserved", [])}

​

def gate(trace, instance, tol_bounds, tol_cons, label):

import numpy as np

arr = np.asarray(trace, dtype=float)

nt, nobs = len(instance["times"]), len(instance["observables"])

if arr.shape != (nt, nobs): return f"{label}: output shape {arr.shape} differs from ({nt}, {nobs})"

if not np.all(np.isfinite(arr)): return f"{label}: output contains non-finite values"

if np.any(arr < -1 - tol_bounds) or np.any(arr > 1 + tol_bounds): return f"{label}: a Pauli expectation value lies outside [-1, 1]"

bits = np.array(instance["initial_state"]["bits"]); z0 = 1.0 - 2.0 * bits

if np.max(np.abs(arr[0] - z0)) > tol_cons: return f"{label}: the t = 0 row does not match the initial product state"

for c in instance.get("conserved", []):

if c == "total_z":

tot = arr.sum(axis=1)

if np.max(np.abs(tot - tot[0])) > tol_cons * nobs: return f"{label}: total magnetisation is not conserved (max drift {float(np.max(np.abs(tot - tot[0]))):.3g})"

return None

​