1import json
  2import math
  3import numpy as np
  4from scipy.linalg import expm
  5
  6SEED = 1468
  7rng = np.random.default_rng(SEED)
  8LAMBDA = 1e-8
  9
 10
 11def symvec(S):
 12    return np.array([S[0, 0], S[1, 1], S[2, 2], math.sqrt(2)*S[0, 1],
 13                     math.sqrt(2)*S[0, 2], math.sqrt(2)*S[1, 2]], dtype=float)
 14
 15
 16def symmat(c):
 17    S = np.zeros((3, 3), dtype=float)
 18    S[0, 0], S[1, 1], S[2, 2] = c[:3]
 19    S[0, 1] = S[1, 0] = c[3]/math.sqrt(2)
 20    S[0, 2] = S[2, 0] = c[4]/math.sqrt(2)
 21    S[1, 2] = S[2, 1] = c[5]/math.sqrt(2)
 22    return S
 23
 24
 25def axes_icosa():
 26    p = (1 + math.sqrt(5)) / 2
 27    raw = [(0,1,p), (1,p,0), (p,0,1), (0,1,-p),
 28           (1,-p,0), (p,0,-1)]
 29    return np.asarray(raw, dtype=float) / math.sqrt(1 + p*p)
 30
 31
 32def axes_coordinate():
 33    return np.eye(3)
 34
 35
 36def axes_random(n=6):
 37    x = rng.normal(size=(n, 3))
 38    return x / np.linalg.norm(x, axis=1, keepdims=True)
 39
 40
 41def frame(v, F):
 42    w = v @ F.T
 43    # Each row is the Frobenius-orthonormal vectorization of ww^T.
 44    A = np.stack([symvec(np.outer(q, q)) for q in w], axis=0)
 45    return A
 46
 47
 48def sl_deformation(anisotropy, rotation=True):
 49    H = np.diag([anisotropy, -anisotropy/2, -anisotropy/2])
 50    if rotation:
 51        Q, _ = np.linalg.qr(rng.normal(size=(3, 3)))
 52        if np.linalg.det(Q) < 0: Q[:, 0] *= -1
 53        H = Q @ H @ Q.T
 54    F = expm(H)
 55    # exp(trace-free H) has determinant one up to floating point precision.
 56    return F
 57
 58
 59def random_tracefree():
 60    S = rng.normal(size=(3, 3)); S = (S + S.T)/2
 61    return S - np.trace(S)*np.eye(3)/3
 62
 63
 64def reconstruct(v, F, S, noise_std=0.0, ridge=LAMBDA):
 65    A = frame(v, F)
 66    y = A @ symvec(S) + rng.normal(0, noise_std, size=len(v))
 67    c = np.linalg.solve(A.T @ A + ridge*np.eye(6), A.T @ y)
 68    # This is the stated optional trace-free postprocessing.
 69    Sh = symmat(c)
 70    Sh -= np.trace(Sh)*np.eye(3)/3
 71    return Sh, A
 72
 73
 74def rank_and_svals(v, F):
 75    A = frame(v, F)
 76    return np.linalg.matrix_rank(A, tol=1e-10), np.linalg.svd(A, compute_uv=False, full_matrices=False)
 77
 78
 79def main():
 80    ico = axes_icosa(); coord = axes_coordinate()
 81    # A fixed random frame is used as a secondary practical comparator.
 82    rand6 = axes_random()
 83    results = {"seed": SEED, "lambda": LAMBDA}
 84
 85    # Prediction 1: at identity, the six icosa projectors span all Sym(3),
 86    # while coordinate projectors have rank exactly three.
 87    r_i, sv_i = rank_and_svals(ico, np.eye(3))
 88    r_c, sv_c = rank_and_svals(coord, np.eye(3))
 89    results["prediction_1_identity_rank"] = {
 90        "predicted": "ico rank=6; coordinate rank=3",
 91        "observed": {"ico_rank": r_i, "coordinate_rank": r_c,
 92                     "ico_sigma_min": float(sv_i[-1]),
 93                     "coordinate_smallest_singular_value": 0.0}}
 94
 95    # Prediction 2: congruence by any invertible F preserves rank. Sweep SL(3)
 96    # deformation strengths and record ranks/minimum singular values.
 97    strengths = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]
 98    rank_sweep = []
 99    for a in strengths:
100        F = sl_deformation(a)
101        ri, si = rank_and_svals(ico, F)
102        rc, sc = rank_and_svals(coord, F)
103        rank_sweep.append({"anisotropy": a, "ico_rank": ri,
104                           "coordinate_rank": rc,
105                           "ico_sigma_min": float(si[-1]),
106                           "coordinate_smallest_nonzero_singular_value": float(sc[-1]),
107                           "coordinate_has_zero_singular_values": True})
108    results["prediction_2_rank_preservation"] = {
109        "predicted": "ico rank remains 6 for every invertible F; coordinate remains rank 3",
110        "observed": rank_sweep}
111
112    # Prediction 3: at fixed measurement noise, reconstruction error tracks
113    # inverse conditioning and rises under strong anisotropy; compare frames.
114    noise = 1e-3
115    trials = 100
116    error_sweep = []
117    # reset independent RNG sequence is not needed; all results deterministic.
118    for a in strengths:
119        row = {"anisotropy": a}
120        for name, v in [("icosahedral", ico), ("random6", rand6), ("coordinate", coord)]:
121            errs, invmins, conds = [], [], []
122            for _ in range(trials):
123                F = sl_deformation(a)
124                S = random_tracefree()
125                Sh, A = reconstruct(v, F, S, noise_std=noise)
126                errs.append(np.linalg.norm(Sh-S) / np.linalg.norm(S))
127                ss = np.linalg.svd(A, compute_uv=False)
128                invmins.append(1.0 / max(ss[-1], 1e-30))
129                conds.append(float(ss[0] / max(ss[-1], 1e-30)))
130            row[name] = {"relative_error_mean": float(np.mean(errs)),
131                         "relative_error_std": float(np.std(errs)),
132                         "inv_sigma_min_mean": float(np.mean(invmins)),
133                         "condition_mean": float(np.mean(conds))}
134        error_sweep.append(row)
135    results["prediction_3_noise_conditioning"] = {
136        "predicted": "noise error increases with inverse smallest singular value; coordinate is non-identifiable",
137        "noise_std": noise, "trials": trials, "observed": error_sweep}
138
139    # Exact noiseless reconstruction sanity check at random SL deformations.
140    exact = {}
141    for name, v in [("icosahedral", ico), ("random6", rand6)]:
142        es = []
143        for _ in range(100):
144            F = sl_deformation(rng.uniform(0, 3)); S = random_tracefree()
145            Sh, _ = reconstruct(v, F, S, noise_std=0.0, ridge=0.0)
146            es.append(np.linalg.norm(Sh-S)/np.linalg.norm(S))
147        exact[name] = float(np.max(es))
148    results["exact_reconstruction_max_relative_error"] = exact
149
150    # Compact verdict-relevant aggregates.
151    ico_err = [x["icosahedral"]["relative_error_mean"] for x in error_sweep]
152    ico_inv = [x["icosahedral"]["inv_sigma_min_mean"] for x in error_sweep]
153    results["quantitative_predictions"] = {
154        "P1_identity_rank": {"predicted_icosa_rank": 6, "observed_icosa_rank": int(r_i),
155                             "predicted_coordinate_rank": 3, "observed_coordinate_rank": int(r_c)},
156        "P2_SL3_rank_invariance": {"predicted_icosa_rank_at_all_strengths": 6,
157                                    "observed_icosa_ranks": [int(x["ico_rank"]) for x in rank_sweep],
158                                    "predicted_coordinate_rank_at_all_strengths": 3,
159                                    "observed_coordinate_ranks": [int(x["coordinate_rank"]) for x in rank_sweep]},
160        "P3_noise_amplification": {"prediction": "error rises as inverse smallest singular value rises",
161                                    "strengths": strengths,
162                                    "ico_error_at_strength_0_and_3": [ico_err[0], ico_err[-1]],
163                                    "ico_inverse_sigma_at_strength_0_and_3": [ico_inv[0], ico_inv[-1]],
164                                    "error_growth_factor": ico_err[-1] / ico_err[0],
165                                    "inverse_sigma_growth_factor": ico_inv[-1] / ico_inv[0]}}
166    results["summary"] = {
167        "identity_ico_full_rank": bool(r_i == 6),
168        "identity_coordinate_deficient": bool(r_c < 6),
169        "all_ico_swept_full_rank": bool(all(x["ico_rank"] == 6 for x in rank_sweep)),
170        "all_coordinate_swept_deficient": bool(all(x["coordinate_rank"] < 6 for x in rank_sweep)),
171        "exact_ico_error": exact["icosahedral"],
172        "exact_random6_error": exact["random6"]}
173    with open("results.json", "w") as f:
174        json.dump(results, f, indent=2, default=lambda x: x.item() if isinstance(x, np.generic) else x)
175    print(json.dumps(results, indent=2, default=lambda x: x.item() if isinstance(x, np.generic) else x))
176
177if __name__ == "__main__":
178    main()