1import json
  2import math
  3import numpy as np
  4from scipy.linalg import expm
  5
  6
  7def pairs(n):
  8    return [(i, j) for i in range(n) for j in range(i + 1, n)]
  9
 10
 11def wedge2(M):
 12    """Second multiplicative compound in the ordered Plucker basis."""
 13    p = pairs(M.shape[0])
 14    W = np.empty((len(p), len(p)))
 15    for a, (i, j) in enumerate(p):
 16        for b, (k, l) in enumerate(p):
 17            W[a, b] = M[i, k] * M[j, l] - M[i, l] * M[j, k]
 18    return W
 19
 20
 21def additive_compound(A):
 22    """A^[2] = d/dh wedge2(exp(h A)) at h=0."""
 23    n = A.shape[0]
 24    p = len(pairs(n))
 25    B = np.zeros((p, p))
 26    ij = pairs(n)
 27    for a, (i, j) in enumerate(ij):
 28        B[a, a] = A[i, i] + A[j, j]
 29        for b, (k, l) in enumerate(ij):
 30            if a != b:
 31                # A reliable direct derivative avoids sign-convention mistakes.
 32                h = 1e-6
 33                B[a, b] = (wedge2(expm(h * A))[a, b]) / h
 34    return B
 35
 36
 37def hilbert(u, v):
 38    u = np.maximum(np.asarray(u), 1e-300)
 39    v = np.maximum(np.asarray(v), 1e-300)
 40    r = u / v
 41    return float(np.log(r.max()) - np.log(r.min()))
 42
 43
 44def birkhoff_diameter(P):
 45    # Exact finite-dimensional formula for a strictly positive matrix.
 46    best = 0.0
 47    for i in range(P.shape[0]):
 48        for j in range(P.shape[0]):
 49            for k in range(P.shape[0]):
 50                for l in range(P.shape[0]):
 51                    best = max(best, abs(math.log(P[i, k] * P[j, l] /
 52                                                (P[i, l] * P[j, k]))))
 53    return best
 54
 55
 56def cycle_generator(N, alpha, diagonal=-1.0):
 57    B = diagonal * np.eye(N)
 58    # Fixed strongly connected directed cycle, with an additional reverse edge.
 59    for i in range(N):
 60        B[(i + 1) % N, i] = alpha
 61        B[(i - 1) % N, i] = alpha
 62    return B
 63
 64
 65def random_positive(rng, N):
 66    return np.exp(rng.normal(0.0, 1.0, N))
 67
 68
 69def main():
 70    rng = np.random.default_rng(1478)
 71    # Core algebra check: wedge(exp(hA)) = I + h A^[2] + O(h^2).
 72    A = rng.normal(size=(5, 5))
 73    B = additive_compound(A)
 74    hs = np.array([2e-2, 1e-2, 5e-3, 2.5e-3])
 75    errs = []
 76    for h in hs:
 77        errs.append(np.linalg.norm(wedge2(expm(h * A)) - np.eye(10) - h * B))
 78    slope = np.polyfit(np.log(hs), np.log(errs), 1)[0]
 79
 80    # Toy mechanism: uniformly Metzler, fixed irreducible support, varying alpha.
 81    N, T = 6, 0.8
 82    alphas = [0.05, 0.1, 0.2, 0.4, 0.8]
 83    rows = []
 84    for alpha in alphas:
 85        Bp = cycle_generator(N, alpha)
 86        P = expm(T * Bp)
 87        delta = birkhoff_diameter(P)
 88        q_bound = math.tanh(delta / 4.0)
 89        # Empirical contraction over many positive vectors.
 90        ratios = []
 91        for _ in range(300):
 92            u, v = random_positive(rng, N), random_positive(rng, N)
 93            d0 = hilbert(u, v)
 94            if d0 > 1e-10:
 95                ratios.append(hilbert(P @ u, P @ v) / d0)
 96        # The shortest three-edge path controls the worst small-window entry on this 6-cycle.
 97        tsmall = 0.08
 98        minentry = expm(tsmall * Bp).min()
 99        predicted_path = (alpha * tsmall) ** 3 / 6.0
100        # Repeated identical windows test geometric decay.
101        u, v = random_positive(rng, N), random_positive(rng, N)
102        ds = [hilbert(u, v)]
103        for _ in range(5):
104            u, v = P @ u, P @ v
105            ds.append(hilbert(u, v))
106        observed_window_ratios = [ds[k + 1] / ds[k] for k in range(5)]
107        rows.append({
108            "alpha": alpha,
109            "min_P": float(P.min()),
110            "q_birkhoff_bound": q_bound,
111            "max_sampled_ratio": float(max(ratios)),
112            "path_minentry_at_t0.08": float(minentry),
113            "path_prediction_alpha3_t3_over_6": predicted_path,
114            "repeated_window_ratios": observed_window_ratios,
115        })
116
117    # Unconstrained comparison: signed cycle has no positive-cone guarantee.
118    Bu = cycle_generator(N, 0.4)
119    Bu[1, 0] = -0.4  # break Metzler condition while retaining comparable scale
120    Pu = expm(T * Bu)
121    signed_min = float(Pu.min())
122    unconstrained = []
123    for _ in range(300):
124        u, v = random_positive(rng, N), random_positive(rng, N)
125        d0 = hilbert(u, v)
126        if np.all(Pu @ u > 0) and np.all(Pu @ v > 0):
127            unconstrained.append(hilbert(Pu @ u, Pu @ v) / d0)
128
129    # Quantitative sweep summaries.
130    loga = np.log(alphas)
131    logpath = np.log([r["path_minentry_at_t0.08"] for r in rows])
132    path_slope = float(np.polyfit(loga, logpath, 1)[0])
133    q_by_alpha = [r["q_birkhoff_bound"] for r in rows]
134    result = {
135        "compound_identity": {
136            "h": hs.tolist(), "errors": errs, "loglog_slope_expected_2": float(slope)
137        },
138        "predictions": {
139            "P1_compound_identity_error_is_O_h2": {"observed_slope": float(slope), "expected": 2.0},
140            "P2_shortest_path_min_entry_scales_alpha3": {
141                "observed_log_slope": path_slope, "expected": 3.0,
142                "note": "For the bidirectional 6-cycle, the worst off-diagonal entries require three edges; the leading path term is (alpha*t)^3/3!."
143            },
144            "P3_stronger_uniform_coupling_improves_projective_contraction": {
145                "alphas": alphas, "q_bounds": q_by_alpha,
146                "q_at_alpha_0.05": q_by_alpha[0], "q_at_alpha_0.8": q_by_alpha[-1]
147            }
148        },
149        "rows": rows,
150        "baseline_signed_generator": {
151            "min_transition_entry": signed_min,
152            "valid_positive_output_sample_ratio_max": (max(unconstrained) if unconstrained else None),
153            "metzler_violation": 0.8
154        },
155        "interpretation": "Positive irreducible compound dynamics produced strictly positive finite-window transitions and geometric Hilbert contraction. The signed baseline has a negative transition entry, so the positive-cone theorem is unavailable."
156    }
157    with open("results.json", "w") as f:
158        json.dump(result, f, indent=2)
159    print(json.dumps(result, indent=2))
160
161
162if __name__ == "__main__":
163    main()