1import json
  2import math
  3from pathlib import Path
  4import numpy as np
  5
  6
  7def paper_jacobian(t):
  8    t1, t2, t3, t4, t5, t6, t7, t8, t9 = np.asarray(t, dtype=float)
  9    return np.array([
 10        [-t1+t6+t7, 0, 0, 0, 0, t1-t6+t7, t1+t6-t7, 0, 0],
 11        [0, -t2+t3+t9, t2-t3+t9, 0, 0, 0, 0, 0, t2+t3-t9],
 12        [0, 0, 0, -t4+t5+t8, t4-t5+t8, 0, 0, t4+t5-t8, 0],
 13        [0, 0, 0, 0, 0, 0, -t7+t8+t9, t7-t8+t9, t7+t8-t9],
 14    ], dtype=float)
 15
 16
 17def numerical_rank(a, delta=1e-7):
 18    singular = np.linalg.svd(np.asarray(a), compute_uv=False)
 19    return int(np.sum(singular > delta)), singular
 20
 21
 22def triangle_area(lengths):
 23    a, b, c = map(float, lengths)
 24    s = (a + b + c) / 2.0
 25    return math.sqrt(max(s * (s-a) * (s-b) * (s-c), 1e-14))
 26
 27
 28def triangle_area_gradient(lengths):
 29    # Exact gradient of Heron's area with respect to (a,b,c).
 30    a, b, c = map(float, lengths)
 31    area = triangle_area(lengths)
 32    return np.array([
 33        a * (b*b + c*c - a*a),
 34        b * (a*a + c*c - b*b),
 35        c * (a*a + b*b - c*c),
 36    ]) / (8.0 * area)
 37
 38
 39def finite_difference_check():
 40    lengths = np.array([1.1, 1.3, 1.6])
 41    eps = 1e-6
 42    analytic = triangle_area_gradient(lengths)
 43    numeric = np.array([
 44        (triangle_area(lengths + eps*np.eye(3)[i]) -
 45         triangle_area(lengths - eps*np.eye(3)[i])) / (2*eps)
 46        for i in range(3)
 47    ])
 48    return {
 49        "analytic": analytic.tolist(), "finite_difference": numeric.tolist(),
 50        "max_abs_error": float(np.max(np.abs(analytic - numeric))),
 51    }
 52
 53
 54def make_candidates(seed=7, n_edges=30, n_candidates=180):
 55    rng = np.random.default_rng(seed)
 56    rows, lengths = [], []
 57    for _ in range(n_candidates):
 58        center = int(rng.integers(0, n_edges))
 59        tri = np.array([(center + int(rng.integers(-3, 4))) % n_edges,
 60                        int(rng.integers(0, n_edges)), int(rng.integers(0, n_edges))])
 61        while len(set(tri.tolist())) < 3:
 62            tri = rng.choice(n_edges, size=3, replace=False)
 63        pts = rng.normal(size=(3, 2))
 64        ls = np.array([np.linalg.norm(pts[1]-pts[2]),
 65                       np.linalg.norm(pts[0]-pts[2]),
 66                       np.linalg.norm(pts[0]-pts[1])])
 67        row = np.zeros(n_edges)
 68        row[tri] = triangle_area_gradient(ls)
 69        rows.append(row)
 70        lengths.append(ls)
 71    return np.asarray(rows), np.asarray(lengths)
 72
 73
 74def logdet_information(rows, lam=1e-3):
 75    gram = lam * np.eye(rows.shape[1]) + rows.T @ rows
 76    sign, value = np.linalg.slogdet(gram)
 77    return float(value) if sign > 0 else -np.inf
 78
 79
 80def greedy_select(rows, budget, lam=1e-3):
 81    selected, remaining = [], list(range(len(rows)))
 82    info = lam * np.eye(rows.shape[1])
 83    rank_trace = []
 84    for _ in range(min(budget, len(rows))):
 85        inv = np.linalg.inv(info)
 86        gains = [math.log1p(float(rows[i] @ inv @ rows[i])) for i in remaining]
 87        pos = int(np.argmax(gains))
 88        idx = remaining.pop(pos)
 89        selected.append(idx)
 90        info += np.outer(rows[idx], rows[idx])
 91        rank_trace.append(numerical_rank(rows[selected])[0])
 92    return selected, rank_trace
 93
 94
 95def random_select(rows, budget, rng):
 96    return rng.choice(len(rows), size=budget, replace=False).tolist()
 97
 98
 99def mean_cosine(rows):
100    z = rows / np.maximum(np.linalg.norm(rows, axis=1, keepdims=True), 1e-12)
101    vals = [abs(float(z[i] @ z[j])) for i in range(len(z)) for j in range(i)]
102    return float(np.mean(vals)) if vals else 0.0
103
104
105def metrics(rows, selected, lam=1e-3):
106    chosen = rows[selected]
107    rank, singular = numerical_rank(chosen)
108    return {
109        "rank": rank,
110        "logdet_information": logdet_information(chosen, lam),
111        "mean_pairwise_cosine_abs": mean_cosine(chosen),
112        "smallest_singular_value": float(singular[-1]) if len(singular) else 0.0,
113    }
114
115
116def main():
117    # Explicit witness from the paper: columns 1,2,4,7 are 1-indexed.
118    witness = paper_jacobian([2, 1, 0, 1, 0, 0, 1, 0, 0])
119    witness_rank, witness_singular = numerical_rank(witness)
120    witness_minor_det = float(np.linalg.det(witness[:, [0, 1, 3, 6]]))
121
122    rows, lengths = make_candidates()
123    rng = np.random.default_rng(123)
124    fractions = [0.10, 0.25, 0.50]
125    sweep = {}
126    for fraction in fractions:
127        budget = int(round(fraction * len(rows)))
128        selected, rank_trace = greedy_select(rows, budget)
129        random_metrics = [metrics(rows, random_select(rows, budget, rng)) for _ in range(40)]
130        keys = ["rank", "logdet_information", "mean_pairwise_cosine_abs", "smallest_singular_value"]
131        sweep[str(int(100*fraction)) + "%"] = {
132            "budget": budget,
133            "greedy": metrics(rows, selected),
134            "greedy_rank_trace": rank_trace,
135            "random_mean": {k: float(np.mean([x[k] for x in random_metrics])) for k in keys},
136            "random_std": {k: float(np.std([x[k] for x in random_metrics])) for k in keys},
137            "message_fraction": budget / len(rows),
138        }
139
140    output = {
141        "finite_difference_gradient_check": finite_difference_check(),
142        "paper_witness": {
143            "rank": witness_rank,
144            "singular_values": witness_singular.tolist(),
145            "selected_minor_determinant": witness_minor_det,
146        },
147        "candidate_pool": {"candidates": len(rows), "global_edges": rows.shape[1]},
148        "all_candidates": metrics(rows, list(range(len(rows)))),
149        "budget_sweep": sweep,
150        "interpretation": "Greedy selects rank-increasing rows and maximizes regularized global edge information; random is the equal-budget control.",
151    }
152    Path("results.json").write_text(json.dumps(output, indent=2))
153    print(json.dumps(output, indent=2))
154
155
156if __name__ == "__main__":
157    main()