1import json
 2import numpy as np
 3
 4OMEGA = 1.0 / 16.0
 5
 6
 7def heis(p, q):
 8    x, y, z = p
 9    X, Y, Z = q
10    return np.array([x + X, y + Y, z + Z + 0.5 * (x * Y - y * X)])
11
12
13def four(v):
14    return -OMEGA*v[0] + (0.5 + OMEGA)*(v[1] + v[2]) - OMEGA*v[3]
15
16
17def refine(seq, geometric=True):
18    """Odd/even refinement. Boundary inserted points use linear interpolation."""
19    n, d = seq.shape
20    out = np.empty((2*n-1, d), dtype=float)
21    out[::2] = seq
22    for j in range(n-1):
23        if 1 <= j <= n-3:
24            w = seq[j-1:j+3]
25            a, b = four(w[:, 0]), four(w[:, 1])
26            z = four(w[:, 2])
27            if geometric:
28                z += 0.5 * (seq[j, 0] * b - seq[j, 1] * a)
29            out[2*j+1] = [a, b, z]
30        else:
31            out[2*j+1] = 0.5 * (seq[j] + seq[j+1])
32    return out
33
34
35def path_z(x, y):
36    """Discrete signed area accumulated along a fine path, starting at zero."""
37    z = np.zeros(len(x))
38    z[1:] = np.cumsum(0.5 * (x[:-1]*y[1:] - y[:-1]*x[1:]))
39    return z
40
41
42def make_paths(seed=7, count=100, n=65):
43    rng = np.random.default_rng(seed)
44    t = np.linspace(0, 1, n)
45    paths = []
46    for _ in range(count):
47        # Smooth closed-ish controls with nontrivial winding and signed area.
48        phase = rng.uniform(0, 2*np.pi, 4)
49        ax, ay = rng.uniform(.7, 1.3, 2)
50        x = ax*np.cos(2*np.pi*t + phase[0]) + .18*np.sin(6*np.pi*t + phase[1])
51        y = ay*np.sin(2*np.pi*t + phase[2]) + .18*np.cos(4*np.pi*t + phase[3])
52        z = path_z(x, y)
53        paths.append(np.stack([x, y, z], axis=1))
54    return paths
55
56
57def main():
58    # Core algebra: associativity should hold to floating point precision.
59    rng = np.random.default_rng(123)
60    assoc = []
61    for _ in range(1000):
62        p, q, r = rng.normal(size=(3, 3))
63        assoc.append(np.max(np.abs(heis(heis(p,q),r)-heis(p,heis(q,r)))))
64    assoc_err = max(assoc)
65
66    # Reversal: reversing the traversal changes the signed polygonal area sign.
67    x = np.array([0., 1., 1., 0., 0.]); y = np.array([0., 0., 1., 1., 0.])
68    area = path_z(x, y)[-1]
69    reverse_area = path_z(x[::-1], y[::-1])[-1]
70
71    paths = make_paths()
72    errors = {"linear": [], "four_point": [], "heisenberg": []}
73    z_ranges = {k: [] for k in errors}
74    for fine in paths:
75        coarse = fine[::2]
76        # Linear control is standard interpolation; ordinary four-point is the direct baseline.
77        lin = np.empty_like(fine); lin[::2] = coarse
78        for j in range(len(coarse)-1): lin[2*j+1] = .5*(coarse[j]+coarse[j+1])
79        fp = refine(coarse, geometric=False)
80        hp = refine(coarse, geometric=True)
81        for name, pred in [("linear",lin),("four_point",fp),("heisenberg",hp)]:
82            errors[name].append(np.sqrt(np.mean((pred[:,2]-fine[:,2])**2)))
83            z_ranges[name].append(float(np.max(np.abs(pred[:,2]))))
84    summary = {
85        "associativity_max_abs_error": assoc_err,
86        "unit_square_area": area,
87        "reversed_unit_square_area": reverse_area,
88        "reversal_sum": area + reverse_area,
89        "rmse_mean": {k: float(np.mean(v)) for k,v in errors.items()},
90        "rmse_median": {k: float(np.median(v)) for k,v in errors.items()},
91        "max_abs_z_mean": {k: float(np.mean(v)) for k,v in z_ranges.items()},
92        "n_paths": len(paths), "fine_length": len(paths[0]), "coarse_length": len(paths[0][::2]),
93    }
94    print(json.dumps(summary, indent=2))
95
96if __name__ == '__main__':
97    main()