1import json
 2import numpy as np
 3from cross_ratio_lattice import *
 4
 5
 6def stats(x):
 7    x = np.asarray(x, float)
 8    return {'max': float(np.nanmax(x)), 'median': float(np.nanmedian(x)), 'mean': float(np.nanmean(x))}
 9
10rng = np.random.default_rng(123)
11# Direct formula and inverse formula on independent nonsingular triples.
12a = rng.normal(size=10000) + 1j*rng.normal(size=10000)
13b = rng.normal(size=10000) + 1j*rng.normal(size=10000)
14c = rng.normal(size=10000) + 1j*rng.normal(size=10000)
15d, bad = complete_d(a, b, c)
16finite = np.isfinite(d)
17formula_res = np.abs(cross_ratio(a[finite], b[finite], c[finite], d[finite]) + 1)
18c2, bad2 = complete_c(a, b, d)
19reverse_err = np.abs(c2[finite] - c[finite])
20results = {'formula': {'singular': int(bad.sum()), 'residual': stats(formula_res),
21                       'inverse_abs_error': stats(reverse_err)}}
22
23for dtype in (np.complex128, np.complex64):
24    z, singular = generate_lattice(18, seed=8, dtype=dtype)
25    # Cast noise to the lattice precision: do not accidentally promote complex64.
26    real_dtype = np.float32 if dtype == np.complex64 else np.float64
27    noise = (rng.normal(size=z.shape[0]).astype(real_dtype) +
28             1j * rng.normal(size=z.shape[0]).astype(real_dtype)).astype(dtype)
29    amp = dtype.type(1e-4) if hasattr(dtype, 'type') else 1e-4
30    top = (z[0, :] + amp * noise).astype(dtype)
31    left = (z[:, 0] + amp * noise).astype(dtype)
32    left[0] = top[0]
33    constrained = constrained_from_boundary(top, left)
34    additive = additive_from_boundary(top, left)
35    results[str(dtype)] = {
36        'generation_singular': singular,
37        'exact_residual': stats(plaquette_residuals(z)),
38        'constrained_perturbed_residual': stats(plaquette_residuals(constrained)),
39        'additive_control_residual': stats(plaquette_residuals(additive)),
40        'constrained_finite': bool(np.isfinite(constrained).all())}
41
42# Completion followed by inverse completion measures the claimed reversibility.
43for dtype in (np.complex128, np.complex64):
44    rd = np.float32 if dtype == np.complex64 else np.float64
45    aa = (rng.normal(size=5000).astype(rd) + 1j*rng.normal(size=5000).astype(rd)).astype(dtype)
46    bb = (rng.normal(size=5000).astype(rd) + 1j*rng.normal(size=5000).astype(rd)).astype(dtype)
47    cc = (rng.normal(size=5000).astype(rd) + 1j*rng.normal(size=5000).astype(rd)).astype(dtype)
48    dd, bad = complete_d(aa, bb, cc)
49    cb, _ = complete_c(aa, bb, dd)
50    ok = np.isfinite(dd)
51    results['roundtrip_' + str(dtype)] = {
52        'finite': int(ok.sum()),
53        'relative_error': float(np.linalg.norm((cb-cc)[ok]) / np.linalg.norm(cc[ok]))}
54
55with open('results.json', 'w') as f:
56    json.dump(results, f, indent=2)
57print(json.dumps(results, indent=2))