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