1import json, time
  2from pathlib import Path
  3import numpy as np
  4from scipy.spatial import ConvexHull
  5
  6SEED = 206
  7rng = np.random.default_rng(SEED)
  8N_DIR = 64
  9angles = np.linspace(0, 2*np.pi, N_DIR, endpoint=False)
 10DIRECTIONS = np.c_[np.cos(angles), np.sin(angles)]
 11
 12
 13def hull(points):
 14    p = np.asarray(points, float)
 15    if len(p) <= 2:
 16        return p
 17    return p[ConvexHull(p).vertices]
 18
 19
 20def random_polygon(n=12):
 21    # A compact, nondegenerate convex object represented by its hull vertices.
 22    return hull(rng.normal(size=(n, 2)) * np.array([1.0, .7]) + rng.normal(size=2))
 23
 24
 25def support_tuple(poly, dirs=DIRECTIONS):
 26    p = np.asarray(poly)
 27    # deterministic tie break is adequate away from measure-zero ties
 28    return p[np.argmax(p @ dirs.T, axis=0)]
 29
 30
 31def decode(x):
 32    return hull(np.asarray(x))
 33
 34
 35def tuple_add(x, y, alpha=1., beta=1.):
 36    return alpha * np.asarray(x) + beta * np.asarray(y)
 37
 38
 39def generic_minkowski(a, b):
 40    return hull((a[:, None, :] + b[None, :, :]).reshape(-1, 2))
 41
 42
 43def point_segment_distance(p, a, b):
 44    v = b-a
 45    den = np.dot(v, v)
 46    t = 0.0 if den == 0 else np.clip(np.dot(p-a, v)/den, 0, 1)
 47    return np.linalg.norm(p-(a+t*v))
 48
 49
 50def directed_poly_distance(a, b):
 51    if len(b) == 1:
 52        return max(np.linalg.norm(x-b[0]) for x in a)
 53    return max(min(point_segment_distance(x, b[i], b[(i+1)%len(b)]) for i in range(len(b))) for x in a)
 54
 55
 56def hausdorff(a, b):
 57    return max(directed_poly_distance(a,b), directed_poly_distance(b,a))
 58
 59
 60def tuple_distance(x, y):
 61    return np.max(np.linalg.norm(np.asarray(x)-np.asarray(y), axis=1))
 62
 63
 64def verify_math(trials=300):
 65    add_err, scale_err, lip_ratios = [], [], []
 66    for _ in range(trials):
 67        a, b = random_polygon(), random_polygon()
 68        xa, xb = support_tuple(a), support_tuple(b)
 69        # Support identity, checked by comparing every directional support point.
 70        lhs = support_tuple(generic_minkowski(a,b))
 71        rhs = xa + xb
 72        add_err.append(np.max(np.linalg.norm(lhs-rhs, axis=1)))
 73        lam = rng.uniform(0, 3)
 74        lhs_s = support_tuple(lam*a)
 75        scale_err.append(np.max(np.linalg.norm(lhs_s-lam*xa, axis=1)))
 76        # Perturb point tuples; convexification should not amplify Hausdorff error.
 77        noise = rng.normal(size=xa.shape)*.03
 78        dec0, dec1 = decode(xa), decode(xa+noise)
 79        lip_ratios.append(hausdorff(dec0,dec1)/(np.max(np.linalg.norm(noise,axis=1))+1e-12))
 80    return {"max_add_identity_error": float(max(add_err)),
 81            "max_scaling_identity_error": float(max(scale_err)),
 82            "max_convexification_ratio": float(max(lip_ratios)),
 83            "median_convexification_ratio": float(np.median(lip_ratios))}
 84
 85
 86def chain_trial(length, n=8):
 87    objs = [random_polygon(n) for _ in range(length)]
 88    # Perturb each latent tuple independently, then compare decoded chain output.
 89    tuples = [support_tuple(x) for x in objs]
 90    noisy = [x + rng.normal(size=x.shape)*0.01 for x in tuples]
 91    structured, structured_noisy = tuples[0], noisy[0]
 92    t0 = time.perf_counter()
 93    for x, xn in zip(tuples[1:], noisy[1:]):
 94        structured = tuple_add(structured, x)
 95        structured_noisy = tuple_add(structured_noisy, xn)
 96    structured_time = time.perf_counter()-t0
 97    # Generic control repeatedly forms all pairwise vertex sums and hulls.
 98    generic = objs[0]
 99    t0 = time.perf_counter()
100    for x in objs[1:]:
101        generic = generic_minkowski(generic, x)
102    generic_time = time.perf_counter()-t0
103    exact = generic
104    struct_poly = decode(structured)
105    noisy_poly = decode(structured_noisy)
106    return {"structured_error": hausdorff(struct_poly, exact),
107            "structured_perturbed_error": hausdorff(noisy_poly, struct_poly),
108            "structured_tuple_bound": tuple_distance(structured_noisy, structured),
109            "generic_vertices": int(len(generic)),
110            "structured_time": structured_time, "generic_time": generic_time}
111
112
113def main():
114    math = verify_math()
115    rows = [chain_trial(k) for k in range(2, 11)]
116    out = {"seed": SEED, "directions": N_DIR, "math": math, "chains": rows}
117    Path("results.json").write_text(json.dumps(out, indent=2))
118    print(json.dumps(out, indent=2))
119
120if __name__ == "__main__":
121    main()