Listing: Listing

#!/usr/bin/env python3
"""Em03-H-v8 tools. New, statically revised, unexecuted artifact."""
import argparse
import contextlib
import datetime
from fractions import Fraction as F
import hashlib
import itertools
import json
import math
import os
import errno
from pathlib import Path
import platform
import shutil
import sys
import tempfile
import time
import traceback

np = None
TOL = 2e-10
PTOL = 1e-12
DEN = 256*65521
TIMES = (0., .25, .5, .75, 1.)
ANGLES = (0., 0., math.pi/4, math.pi/2, math.pi/2)
CONTRACT = {"t0":1, "T":1, "Omega":1,
            "v":1, "w":0, "packet":"fixed_f"}
TRANSIENT = {errno.EINTR, errno.EAGAIN, errno.EBUSY, errno.EIO}

def need(ok, message):
    if not ok:
        raise AssertionError(message)

def sha(path):
    return hashlib.sha256(Path(path).read_bytes()).hexdigest()

def atomic_json(path, data):
    path = Path(path)
    payload = json.dumps(data, indent=2, allow_nan=False)
    for attempt in range(3):
        tmp = None
        try:
            fd, name = tempfile.mkstemp(
                prefix="."+path.name+".", suffix=".tmp", dir=path.parent)
            tmp = Path(name)
            with os.fdopen(fd, "w", encoding="utf-8") as stream:
                stream.write(payload)
                stream.flush()
                os.fsync(stream.fileno())
            os.replace(tmp, path)
            return
        except OSError as exc:
            if attempt == 2 or exc.errno not in TRANSIENT:
                raise
            time.sleep(.05*(attempt+1))
        finally:
            if tmp is not None and tmp.exists():
                try:
                    tmp.unlink()
                except OSError:
                    pass

class Recorder:
    def __init__(self, root, metadata):
        self.root = Path(root)
        self.meta = metadata
        self.records = {}
        self.sequence = 0

    def record(self, name, value):
        self.sequence += 1
        self.records[name] = value
        line = {"sequence":self.sequence, "name":name, "value":value}
        # A journal failure is fatal: do not call the inspection complete.
        with (self.root/"checks.jsonl").open("a", encoding="utf-8") as stream:
            stream.write(json.dumps(line, allow_nan=False)+"\n")
            stream.flush()
            os.fsync(stream.fileno())
        if self.sequence % 25 == 0:
            self.snapshot()

    def snapshot(self):
        atomic_json(self.root/"inspection.json",
                    {"records":self.records, "validated":False,
                     "last_sequence":self.sequence})

    def exact(self, name, ok):
        value = {"passed":bool(ok), "exact_predicate":True}
        self.record(name, value)
        need(ok, "Exact check failed: "+name)

    def call(self, name, function):
        value = {"passed":False}
        try:
            result = function()
            value["passed"] = True
            return result
        except BaseException as exc:
            value["failure"] = repr(exc)
            raise
        finally:
            self.record(name, value)

    def check(self, name, actual, expected, tol=TOL):
        value = {"passed":False, "tolerance":float(tol)}
        try:
            a, b = np.asarray(actual), np.asarray(expected)
            value["actual_shape"] = list(a.shape)
            value["expected_shape"] = list(b.shape)
            need(a.shape == b.shape, "Shape: "+name)
            need(np.all(np.isfinite(a)) and np.all(np.isfinite(b)),
                 "Nonfinite: "+name)
            discrepancy = 0. if a.size == 0 else float(np.max(np.abs(a-b)))
            value["observed_max_abs"] = discrepancy
            value["passed"] = discrepancy <= tol
            need(value["passed"], "Tolerance: "+name)
        finally:
            self.record(name, value)

def execute(args, action, source, numerical=True,
            success="COMPLETED_NOT_VALIDATED"):
    global np
    root = Path(args.out)
    root.mkdir(parents=True, exist_ok=False)
    start = time.monotonic()
    meta = {"status":"RUNNING", "mode":getattr(args,"mode","defect"),
            "argv":sys.argv, "validated":False,
            "python":sys.version, "platform":platform.platform(),
            "start_utc":datetime.datetime.now(datetime.timezone.utc).isoformat()}
    r = Recorder(root, meta)
    primary = None
    try:
        atomic_json(root/"manifest.json", meta)
        source = Path(source)
        shutil.copyfile(source, root/source.name)
        meta["source_name"] = source.name
        meta["source_sha256"] = sha(source)
        atomic_json(root/"manifest.json", meta)
        with (root/"stdout.txt").open("w") as out, \
             (root/"stderr.txt").open("w") as err, \
             contextlib.redirect_stdout(out), contextlib.redirect_stderr(err):
            if numerical:
                import numpy as numpy
                np = numpy
                meta["numpy"] = np.__version__
                atomic_json(root/"numpy_configuration.json",
                            np.show_config(mode="dicts"))
            result = action(args, root, r)
            if result is not None:
                print(json.dumps(result, indent=2, allow_nan=False))
        meta["status"] = success
    except BaseException as exc:
        primary = exc
        meta["status"] = "FAILED"
        meta["failure"] = repr(exc)
        try:
            (root/"traceback.txt").write_text(traceback.format_exc())
        except OSError as io:
            meta["traceback_write_failure"] = repr(io)
            print(traceback.format_exc(), file=sys.__stderr__)
        raise
    finally:
        meta["elapsed_seconds"] = time.monotonic()-start
        meta["end_utc"] = datetime.datetime.now(datetime.timezone.utc).isoformat()
        errors = []
        try:
            r.snapshot()
        except BaseException as exc:
            errors.append({"inspection_snapshot":repr(exc)})
        hashes = {}
        try:
            for p in sorted(root.iterdir()):
                if p.is_file() and p.name != "manifest.json":
                    hashes[p.name] = sha(p)
        except BaseException as exc:
            errors.append({"artifact_hashing":repr(exc)})
        meta["files_sha256"] = hashes
        if errors:
            meta["status"] = "FAILED"
            meta["finalization_errors"] = errors
        try:
            atomic_json(root/"manifest.json", meta)
        except BaseException as exc:
            errors.append({"manifest_write":repr(exc)})
        if errors:
            print("Incomplete persistence: "+repr(errors), file=sys.__stderr__)
            if primary is None:
                raise OSError("Run finalization incomplete")

def load(path):
    with np.load(path, allow_pickle=False) as data:
        return {key:data[key].copy() for key in data.files}

def archive_hashes(src, r):
    src = Path(src)
    manifest = json.loads((src/"manifest.json").read_text())
    r.record("input_manifest_sha256", sha(src/"manifest.json"))
    table = manifest.get("files_sha256")
    need(isinstance(table,dict) and table, "Missing archive hash table")
    for name, expected in table.items():
        r.exact("safe_name_"+name, name not in (".","..") and
                Path(name).name == name)
        r.exact("hash_format_"+name, isinstance(expected,str) and
                len(expected) == 64 and
                all(x in "0123456789abcdef" for x in expected))
        path = src/name
        r.exact("file_exists_"+name, path.is_file())
        observed = sha(path)
        r.record("hash_values_"+name, {"reported":expected,"observed":observed})
        r.exact("hash_match_"+name, observed == expected)
    r.exact("input_completed",
            manifest.get("status") == "COMPLETED_NOT_VALIDATED")
    return manifest

def coefficients():
    rows, a = [], 1
    for q in itertools.combinations(range(1,17),4):
        a = (75*a+74) % 65521
        rows.append({"quartet":list(q),"numerator":2*a-65521})
    return {"denominator":64*65521,"rows":rows}

def tensor_assembly(data):
    eye = np.eye(2,dtype=complex)
    z = np.diag([1,-1]).astype(complex)
    x = np.array([[0,1],[1,0]],complex)
    y = np.array([[0,-1j],[1j,0]],complex)
    generators = []
    for j in range(8):
        for middle in (x,y):
            matrix = np.array([[1]],complex)
            for k in range(8):
                matrix = np.kron(matrix,z if k < j else
                                 middle if k == j else eye)
            perm = np.argmax(np.abs(matrix),axis=0)
            generators.append((perm,matrix[perm,np.arange(256)]))
    re = np.zeros((256,256),dtype=np.int64)
    im = np.zeros_like(re)
    cols = np.arange(256)
    for row in data["rows"]:
        p, value = cols.copy(), np.ones(256,complex)
        for a in reversed(row["quartet"]):
            perm, sign = generators[a-1]
            value *= sign[p]
            p = perm[p]
        re[p,cols] += row["numerator"]*value.real.astype(np.int64)
        im[p,cols] += row["numerator"]*value.imag.astype(np.int64)
    return re,im

def creator_transfer(v,t,angle):
    c = 0. if angle == math.pi/2 else math.cos(angle)
    s = 0. if angle == 0 else math.sin(angle)
    out = np.zeros((v.shape[1],256,256),complex)
    for p in range(256):
        modes = [j for j in range(8) if p & (1 << (7-j))]
        terms = {0:1.+0j}
        for j in reversed(modes):
            newer = {}
            for word,value in terms.items():
                for index,coefficient in ((j,c),(8+j,-1j*s)):
                    bit = 1 << (15-index)
                    need(not (word & bit),"Double occupation")
                    before = (word >> (16-index)).bit_count()
                    target = word | bit
                    newer[target] = newer.get(target,0j)+ \
                        value*coefficient*((-1)**before)
            terms = newer
        free = np.exp(-1j*t*len(modes))
        for word,value in terms.items():
            boundary,bath = word >> 8, word & 255
            sign = (-1)**(boundary.bit_count()*bath.bit_count())
            out[:,bath,boundary] += free*sign*value*v[p]
    return out

def objects(z):
    d = z.shape[0]
    u = np.array([[z[i]@z[j].conj().T for j in range(d)]
                  for i in range(d)])
    c = np.array([[z[i].T@z[j].conj() for j in range(d)]
                  for i in range(d)])
    ju = np.block([[u[i,j]/d for j in range(d)] for i in range(d)])
    jc = np.block([[c[i,j]/d for j in range(d)] for i in range(d)])
    return u,c,ju,jc

def audit_choi(name,j,r):
    r.check(name+"_Hermitian",j,j.conj().T)
    value = float(np.linalg.eigvalsh((j+j.conj().T)/2)[0])
    r.record(name+"_minimum_eigenvalue",{"value":value,"rigorous":False})
    r.exact(name+"_exploratory_PSD",value >= -TOL)

def matrix_json(a):
    return [[[float(z.real),float(z.imag)] for z in row] for row in a]

def inspect_p8(args,root,r):
    src = Path(args.archive)
    manifest = archive_hashes(src,r)
    r.exact("p8_mode",manifest.get("mode") == "p8")
    results = manifest["results"]
    r.exact("stdout_matches_manifest",
            json.loads((src/"stdout.txt").read_text()) == results)
    data = json.loads((src/"coefficients.json").read_text())
    r.exact("frozen_coefficients",data == coefficients())
    re,im = tensor_assembly(data)
    saved = load(src/"integer_hamiltonian.npz")
    r.exact("Hamiltonian_denominator",int(saved["denominator"]) == DEN)
    r.exact("real_assembly",np.array_equal(re,saved["real"]))
    r.exact("imag_assembly",np.array_equal(im,saved["imag"]))
    r.exact("real_symmetry",np.array_equal(re,re.T))
    r.exact("imag_antisymmetry",np.array_equal(im,-im.T))
    even = np.array([p for p in range(256) if p.bit_count()%2 == 0])
    odd = np.array([p for p in range(256) if p.bit_count()%2 == 1])
    r.exact("even_order",np.array_equal(even,saved["even_words"]))
    r.exact("odd_order",np.array_equal(odd,saved["odd_words"]))
    r.exact("parity_blocks",not np.any(re[np.ix_(even,odd)]) and
            not np.any(im[np.ix_(even,odd)]))
    r.exact("trace",sum(int(x) for x in np.diag(re)) == 0 and
            sum(int(x) for x in np.diag(im)) == 0)
    squares = sum(int(x)**2 for x in re.flat)+sum(int(x)**2 for x in im.flat)
    r.exact("Frobenius_identity",squares ==
            256*sum(row["numerator"]**2 for row in data["rows"]))
    h = (re.astype(float)+1j*im.astype(float))/DEN
    spectra = load(src/"spectra_exploratory.npz")
    ee = np.linalg.eigvalsh(h[np.ix_(even,even)])
    eo = np.linalg.eigvalsh(h[np.ix_(odd,odd)])
    r.check("even_spectrum",spectra["even_unshifted"],ee)
    r.check("odd_spectrum",spectra["odd_unshifted"],eo)
    emin = min(ee[0],eo[0])
    for name,value in (("global_minimum_float",emin),
                       ("even_minimum_float",ee[0]),
                       ("odd_minimum_float",eo[0]),
                       ("parity_minimum_difference_float",ee[0]-eo[0])):
        r.check(name,value,results[name])
    r.check("selected_unshifted",ee[:2],results["selected_even_unshifted"])
    r.check("selected_gaps",[ee[1]-ee[0],ee[2]-ee[1]],
            results["selected_neighbor_gaps"])
    code = load(src/"code_and_quench_exploratory.npz")
    v,energy = code["V"],code["energies"]
    r.exact("code_shape",v.shape == (256,2))
    r.check("orthonormality",v.conj().T@v,np.eye(2))
    r.check("odd_code",v[odd],np.zeros((128,2)))
    r.check("shifted_energies",energy,ee[:2]-emin)
    residual = h@v-v*(energy+emin)[None,:]
    r.check("eigenpair_residual",residual,np.zeros_like(v))
    r.record("residual_norms_float",
             [float(np.linalg.norm(residual[:,i])) for i in range(2)])
    anchors = []
    for i in range(2):
        nz = np.flatnonzero(np.abs(v[:,i]) > PTOL)
        r.exact("anchor_exists_%d"%i,len(nz) > 0)
        p = int(nz[0]); anchors.append(p)
        r.check("anchor_imag_%d"%i,v[p,i].imag,0.)
        r.exact("anchor_positive_%d"%i,v[p,i].real > 0)
    r.exact("phase_summary",anchors == results["phase_surrogate_anchors"])
    pop = np.array([p.bit_count() for p in range(256)],float)
    nu = v.conj().T@(pop[:,None]*v)
    mu = v.conj().T@(pop[:,None]**2*v)
    lam = mu-nu@nu
    leak = pop[:,None]*v-v@nu
    for name,value in (("nu",nu),("mu",mu),("Lambda",lam),
                       ("direct_leakage",leak)):
        r.check("saved_"+name,code[name],value)
    r.check("direct_Gram",leak.conj().T@leak,lam)
    r.check("E_code",np.linalg.eigvalsh((nu+nu.conj().T)/2)[-1],
            results["E_code_float"])
    r.check("L_N",np.linalg.norm(leak,2),results["L_N_float"])
    r.record("compressed_matrices",{"nu":matrix_json(nu),
        "mu":matrix_json(mu),"Lambda":matrix_json(lam)})
    for m in (0,2,4,6,8):
        block = v[pop == m]
        occ = load(src/("occupation_%d.npz"%m))
        r.check("occupation_block_%d"%m,occ["projected_code"],block)
        r.check("occupation_Gram_%d"%m,occ["gram"],block.conj().T@block)
        r.check("occupation_singular_%d"%m,np.linalg.svd(block,compute_uv=False),
                results["occupation_data"][str(m)]["singular_values_float"])
    dep = load(src/"departure_exploratory.npz")
    r.check("departure_grid",dep["times"],np.arange(33)/128,tol=0)
    for k,t in enumerate(dep["times"]):
        f = v.conj().T@(np.exp(-1j*t*pop)[:,None]*v)
        r.check("departure_F_%d"%k,dep["F"][k],f)
        gram = f.conj().T@f
        p = 1-np.linalg.eigvalsh((gram+gram.conj().T)/2)[0]
        r.check("departure_value_%d"%k,dep["departure"][k],p)
    r.check("departure_1_8",dep["departure"][16],results["departure_1_over_8"])
    r.check("departure_1_4",dep["departure"][32],results["departure_1_over_4"])
    bases = load(src/"full_bases_exploratory.npz")
    for name,inds,values in (("even",even,spectra["even_unshifted"]),
                            ("odd",odd,spectra["odd_unshifted"])):
        a = bases[name+"_vectors"]
        r.check(name+"_basis_orthonormal",a.conj().T@a,np.eye(128))
        r.check(name+"_basis_residual",h[np.ix_(inds,inds)]@a-a*values[None,:],
                np.zeros_like(a))
    for i in range(2):
        col = bases["even_vectors"][:,i].copy()
        nz = np.flatnonzero(np.abs(col) > PTOL)
        need(len(nz) > 0,"Basis phase anchor")
        col *= np.exp(-1j*np.angle(col[nz[0]]))
        r.check("selected_basis_link_%d"%i,v[even,i],col)
    stored = v*np.exp(-1j*energy)[None,:]
    for k,(t,angle) in enumerate(zip(TIMES,ANGLES)):
        saved = load(src/("channel_%d.npz"%k))
        r.check("frozen_time_%d"%k,saved["time"],t,tol=0)
        r.check("frozen_area_%d"%k,saved["theta"],angle,tol=0)
        z = saved["Z"]
        r.check("creator_reconstruction_%d"%k,z,creator_transfer(stored,t,angle))
        r.check("isometry_%d"%k,np.einsum("irb,jrb->ij",z.conj(),z),np.eye(2))
        u,c,ju,jc = objects(z)
        for name,value in (("units",u),("complement",c),
                           ("Choi_accessible",ju),("Choi_complement",jc)):
            r.check(name+"_%d"%k,saved[name],value)
        for i in range(2):
            for j in range(2):
                for name,obj in (("R",u),("E",c)):
                    r.check("adjoint_%s_%d_%d_%d"%(name,k,i,j),
                            obj[i,j].conj().T,obj[j,i])
                    r.check("trace_%s_%d_%d_%d"%(name,k,i,j),
                            np.trace(obj[i,j]),float(i == j))
        audit_choi("Choi_R_%d"%k,saved["Choi_accessible"],r)
        audit_choi("Choi_E_%d"%k,saved["Choi_complement"],r)
        if k == 4:
            lt = np.exp(-1j*(1+math.pi/2)*pop)[:,None]*stored
            for i in range(2):
                for j in range(2):
                    r.check("endpoint_R_%d_%d"%(i,j),u[i,j],
                            np.outer(lt[:,i],lt[:,j].conj()))
                    target = np.zeros((256,256),complex)
                    target[0,0] = float(i == j)
                    r.check("endpoint_E_%d_%d"%(i,j),c[i,j],target)
    return {"inspection":"completed", "validated":False,
            "exact_code_selection_certified":False,
            "exact_ranks_certified":False}

def g(t):
    return 2/(2+1j*t)

def bump(t):
    t = np.asarray(t,float)
    out = np.zeros_like(t)
    good = (t > .25) & (t < .75)
    out[good] = np.exp(-1/((t[good]-.25)*(.75-t[good])))
    return out

def pulse_b0(panels):
    x = np.arange(panels+1)/panels
    y = bump(x)
    return float((y[0]/2+y[1:-1].sum()+y[-1]/2)/panels)

def volterra(n,b0,v=1,coupled=True):
    t,h = np.arange(n+1)/n,1/n
    kap = math.pi*bump(t)/(2*b0) if coupled else np.zeros(n+1)
    a,b = np.zeros(n+1,complex),np.zeros(n+1,complex)
    a[0] = g(1)
    for i in range(1,n+1):
        weights = np.ones(i); weights[0] = .5
        ra = g(t[i]+1)-1j*h*np.sum(
            weights*g(t[i]-t[:i])*(v*a[:i]+kap[:i]*b[:i]))
        rb = -1j*h*np.sum(
            weights*np.exp(-1j*(t[i]-t[:i]))*kap[:i]*a[:i])
        matrix = np.array([[1+1j*h*v/2,1j*h*kap[i]/2],
                           [1j*h*kap[i]/2,1]],complex)
        a[i],b[i] = np.linalg.solve(matrix,[ra,rb])
    return t,a,b

def spectral(ne,nt,b0):
    x,weights = np.polynomial.laguerre.laggauss(ne)
    need(np.all(weights > 0) and np.all(np.isfinite(weights)),
         "Laguerre weights")
    f,omega = np.sqrt(weights),x/2
    initial = np.r_[f*np.exp(-1j*omega),0j]
    state,dt = initial.copy(),1/nt
    def rhs(t,y):
        a = np.vdot(f,y[:-1])
        k = float(math.pi*bump(np.array([t]))[0]/(2*b0))
        return -1j*np.r_[omega*y[:-1]+f*a+k*f*y[-1],y[-1]+k*a]
    for j in range(nt):
        t = j*dt
        k1 = rhs(t,state)
        k2 = rhs(t+dt/2,state+dt*k1/2)
        k3 = rhs(t+dt/2,state+dt*k2/2)
        k4 = rhs(t+dt,state+dt*k3)
        state += dt*(k1+2*k2+2*k3+k4)/6
    return x,weights,initial,state

def ratio(value):
    n,d = float(value).as_integer_ratio()
    return [str(n),str(d)]

def encode(array):
    return [[ratio(z.real),ratio(z.imag)] for z in array]

def candidate(args,root,r):
    meshes = args.meshes
    need(meshes == sorted(set(meshes)) and
         all(n >= 8 and n & (n-1) == 0 for n in meshes),
         "Increasing power-of-two meshes")
    panels = max(32768,8*max(meshes))
    b0 = pulse_b0(panels)
    need(math.isfinite(b0) and b0 > 0,"Pulse normalization")
    report = {"contract":CONTRACT,"B0_float":b0,"pulse_panels":panels,
              "volterra":[],"spectral":[],"validated":False}
    previous = None
    for n in meshes:
        t,a,b = volterra(n,b0)
        np.savez_compressed(root/("volterra_%d.npz"%n),times=t,A=a,B=b)
        atomic_json(root/("nodes_%d.json"%n),{
            "contract":CONTRACT,"interpolation":"continuous_piecewise_linear",
            "node_representation":"exact_binary64_dyadic_ratios",
            "A":encode(a),"B":encode(b)})
        q = b[-1]
        row = {"n":n,"q":[float(q.real),float(q.imag)],
               "eta":float(abs(q)**2),
               "successive_difference":None if previous is None else
                                       float(abs(q-previous))}
        report["volterra"].append(row)
        r.record("Volterra_%d"%n,row)
        atomic_json(root/"candidate_report.json",report)
        previous = q
    for ne,nt in ((32,2048),(64,4096),(96,8192)):
        x,weights,initial,state = spectral(ne,nt,b0)
        name = "spectral_%d_%d.npz"%(ne,nt)
        np.savez_compressed(root/name,nodes=x,weights=weights,
            initial=initial,final=state,energy_nodes=ne,time_steps=nt,
            B0=b0,t0=1,T=1,v=1,w_control=0)
        q = state[-1]
        row = {"file":name,"energy_nodes":ne,"time_steps":nt,
               "q":[float(q.real),float(q.imag)],
               "finite_system_norm":float(np.vdot(state,state).real),
               "weight_sum":float(weights.sum()),
               "difference_from_finest_volterra":float(abs(q-previous))}
        report["spectral"].append(row)
        r.record(name,row)
        atomic_json(root/"candidate_report.json",report)
    return report

def node_link(values,array):
    need(array.dtype == np.complex128 and array.ndim == 1,
         "Require one-dimensional complex128 nodes")
    need(len(values) == len(array),"Node lengths")
    for item,z in zip(values,array):
        need(len(item) == 2,"Complex node shape")
        for pair,value in zip(item,(z.real,z.imag)):
            need(len(pair) == 2 and int(pair[1]) > 0,"Rational shape")
            need(F(int(pair[0]),int(pair[1])) ==
                 F(*float(value).as_integer_ratio()),"Exact JSON/NPZ mismatch")
    return True

def inspect_candidate(args,root,r):
    src = Path(args.archive)
    archive_hashes(src,r)
    report = json.loads((src/"candidate_report.json").read_text())
    r.exact("candidate_contract",report["contract"] == CONTRACT)
    b0 = report["B0_float"]
    r.check("pulse_replay",b0,pulse_b0(report["pulse_panels"]))
    t,a,b = volterra(128,b0,coupled=False)
    r.check("zero_coupling_B",b,np.zeros_like(b))
    t,a,b = volterra(128,b0,v=0,coupled=False)
    r.check("free_A",a,g(t+1)); r.check("free_B",b,np.zeros_like(b))
    rows = report["volterra"]
    need(rows,"No Volterra results")
    sizes = [row["n"] for row in rows]
    r.exact("mesh_order",sizes == sorted(set(sizes)))
    finest = None
    for row in rows:
        n = row["n"]
        r.exact("power_two_%d"%n,n >= 8 and n & (n-1) == 0)
        path = src/("volterra_%d.npz"%n)
        saved = load(path)
        t,a,b = saved["times"],saved["A"],saved["B"]
        r.check("grid_%d"%n,t,np.arange(n+1)/n,tol=0)
        r.exact("node_shapes_%d"%n,a.shape == (n+1,) and b.shape == (n+1,))
        nodes_path = src/("nodes_%d.json"%n)
        nodes = json.loads(nodes_path.read_text())
        r.exact("JSON_contract_%d"%n,nodes["contract"] == CONTRACT)
        r.exact("JSON_interpolant_%d"%n,
                nodes["interpolation"] == "continuous_piecewise_linear")
        r.call("A_exact_link_%d"%n,lambda:node_link(nodes["A"],a))
        r.call("B_exact_link_%d"%n,lambda:node_link(nodes["B"],b))
        r.record("file_link_%d"%n,{"npz_sha256":sha(path),
            "json_sha256":sha(nodes_path),"exact_node_equality":True})
        tt,aa,bb = volterra(n,b0)
        r.check("A_solver_replay_%d"%n,a,aa)
        r.check("B_solver_replay_%d"%n,b,bb)
        q = b[-1]; finest = q
        r.check("q_summary_%d"%n,q,complex(*row["q"]))
        r.check("eta_summary_%d"%n,abs(q)**2,row["eta"])
        r.check("initial_A_%d"%n,a[0],(4-2j)/5)
        r.check("initial_B_%d"%n,b[0],0)
    for row in report["spectral"]:
        saved = load(src/row["file"])
        ne,nt = row["energy_nodes"],row["time_steps"]
        for key,value in (("energy_nodes",ne),("time_steps",nt),
                          ("t0",1),("T",1),("v",1),("w_control",0)):
            r.check(key+"_"+row["file"],saved[key],value,tol=0)
        r.check("B0_"+row["file"],saved["B0"],b0,tol=0)
        x,w,initial,final = spectral(ne,nt,b0)
        for key,value in (("nodes",x),("weights",w),
                          ("initial",initial),("final",final)):
            r.check(key+"_"+row["file"],saved[key],value)
        r.check("q_"+row["file"],saved["final"][-1],complex(*row["q"]))
        norm = float(np.vdot(saved["final"],saved["final"]).real)
        r.check("norm_summary_"+row["file"],norm,row["finite_system_norm"])
        r.check("finite_system_norm_"+row["file"],norm,1.,tol=1e-7)
        r.check("difference_"+row["file"],abs(saved["final"][-1]-finest),
                row["difference_from_finest_volterra"])
    return {"inspection":"completed","JSON_NPZ_link_checked":True,
            "solver_replay_checked":True,"continuum_validated":False}

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("mode",
        choices=["inspect-p8","candidate","inspect-candidate"])
    parser.add_argument("--archive")
    parser.add_argument("--out",required=True)
    parser.add_argument("--meshes",nargs="+",type=int,default=[256,512])
    args = parser.parse_args()
    if args.mode != "candidate":
        need(args.archive is not None,"--archive required")
    action = {"inspect-p8":inspect_p8,"candidate":candidate,
              "inspect-candidate":inspect_candidate}[args.mode]
    execute(args,action,__file__)

if __name__ == "__main__":
    main()