1"""Toy verification of the cross-degree anti-oscillation certificate.
  2
  3Dynamics: x[t+1] = 1{W x[t] >= R}.  For a bipartite graph, initializing one
  4side to one and the other to zero gives an exact period-two orbit whenever
  5minimum cross-gain >= R (with no within-side edges).  This script checks the
  6boundary, scaling, and threshold suppression predictions, then compares a
  7baseline threshold choice with a certificate-aware choice.
  8"""
  9import itertools
 10import json
 11import numpy as np
 12
 13SEED = 1629
 14
 15
 16def balanced_partitions(n):
 17    # Fix vertex 0 in side 0 to remove complement duplication.
 18    for bits in itertools.product((0, 1), repeat=n - 1):
 19        p = np.array((0,) + bits, dtype=np.int8)
 20        if 0 < p.sum() < n:
 21            yield p
 22
 23
 24def min_cross_gain(W, part):
 25    vals = []
 26    for i in range(len(part)):
 27        opposite = part != part[i]
 28        vals.append(float(np.abs(W[i, opposite]).sum()))
 29    return min(vals)
 30
 31
 32def exact_certificate(W):
 33    vals = [(min_cross_gain(W, p), p.copy()) for p in balanced_partitions(len(W))]
 34    return max(vals, key=lambda z: z[0])
 35
 36
 37def threshold_step(W, x, R):
 38    return (W @ x >= R).astype(np.int8)
 39
 40
 41def orbit(W, x0, R, steps=12):
 42    x = np.asarray(x0, dtype=np.int8).copy()
 43    out = [x.copy()]
 44    for _ in range(steps):
 45        x = threshold_step(W, x, R)
 46        out.append(x.copy())
 47    return np.asarray(out)
 48
 49
 50def is_period_two_from(out, start=0):
 51    return (len(out) >= start + 3 and
 52            np.array_equal(out[start], out[start + 2]) and
 53            not np.array_equal(out[start], out[start + 1]))
 54
 55
 56def make_bipartite(left, right, p, weight=1.0, seed=0):
 57    rng = np.random.default_rng(seed)
 58    n = left + right
 59    W = np.zeros((n, n), dtype=float)
 60    for i in range(left):
 61        for j in range(left, n):
 62            if rng.random() < p:
 63                W[i, j] = W[j, i] = weight
 64    # Ensure every vertex has at least one cross edge for useful sweeps.
 65    for i in range(n):
 66        if not np.any(W[i]):
 67            j = rng.integers(left, n) if i < left else rng.integers(0, left)
 68            W[i, j] = W[j, i] = weight
 69    return W
 70
 71
 72def boundary_sweep(W, x0, cert, thresholds):
 73    return [int(is_period_two_from(orbit(W, x0, R), 0)) for R in thresholds]
 74
 75
 76def main():
 77    # Prediction 1: exact transition is R <= c_*; test every integer R.
 78    W = make_bipartite(4, 4, p=0.72, weight=1.0, seed=1)
 79    c_star, part = exact_certificate(W)
 80    x0 = (part == 1).astype(np.int8)
 81    thresholds = list(range(0, int(c_star) + 3))
 82    observed = boundary_sweep(W, x0, c_star, thresholds)
 83    predicted = [int(1 <= R <= c_star) for R in thresholds]
 84    boundary_ok = observed == predicted
 85
 86    # Prediction 2: scaling all weights scales c_* and the boundary linearly.
 87    scales = [0.5, 1.0, 1.5, 2.0]
 88    scaling_rows = []
 89    for s in scales:
 90        Ws = s * W
 91        cs, ps = exact_certificate(Ws)
 92        xs = (ps == 1).astype(np.int8)
 93        # Probe integer and half-integer thresholds around the predicted c*.
 94        Rs = np.arange(max(0.0, cs - 1.0), cs + 1.01, 0.25)
 95        flags = [is_period_two_from(orbit(Ws, xs, float(R)), 0) for R in Rs]
 96        max_observed = max((float(R) for R, ok in zip(Rs, flags) if ok), default=-1.0)
 97        scaling_rows.append({"scale": s, "predicted_c_star": cs,
 98                             "observed_largest_period2_R": max_observed,
 99                             "boundary_error": abs(max_observed - cs)})
100    scaling_ok = all(row["boundary_error"] <= 0.26 for row in scaling_rows)
101
102    # Prediction 3: choosing R 10% above the certificate blocks this certified
103    # orbit, while choosing R at or below it preserves it.
104    suppress_rows = []
105    for s in [0.5, 1.0, 1.5, 2.0]:
106        Ws = s * W
107        cs, ps = exact_certificate(Ws)
108        xs = (ps == 1).astype(np.int8)
109        at = is_period_two_from(orbit(Ws, xs, cs), 0)
110        above = is_period_two_from(orbit(Ws, xs, 1.1 * cs), 0)
111        suppress_rows.append({"scale": s, "c_star": cs,
112                              "at_certificate": bool(at),
113                              "at_110_percent": bool(above)})
114    suppression_ok = all(r["at_certificate"] and not r["at_110_percent"] for r in suppress_rows)
115
116    # Baseline vs idea on random initial states. Baseline uses R=1, while
117    # certificate-aware selection uses R=1.1*c_* (minimum integer-like safe
118    # threshold is represented directly here because the toy permits reals).
119    rng = np.random.default_rng(SEED)
120    n_trials = 300
121    base_flags, idea_flags = [], []
122    safe_R = 1.1 * c_star
123    for _ in range(n_trials):
124        z = rng.integers(0, 2, size=len(W), dtype=np.int8)
125        base_flags.append(is_period_two_from(orbit(W, z, 1.0), 2))
126        idea_flags.append(is_period_two_from(orbit(W, z, safe_R), 2))
127    baseline_rate = float(np.mean(base_flags))
128    idea_rate = float(np.mean(idea_flags))
129
130    result = {
131        "seed": SEED,
132        "graph_nodes": len(W),
133        "edges": int(np.count_nonzero(W) // 2),
134        "exact_c_star": c_star,
135        "certificate_partition": part.tolist(),
136        "prediction_1_boundary": {"thresholds": thresholds, "observed": observed,
137                                  "predicted": predicted, "confirmed": boundary_ok},
138        "prediction_2_scaling": scaling_rows,
139        "prediction_2_confirmed": scaling_ok,
140        "prediction_3_suppression": suppress_rows,
141        "prediction_3_confirmed": suppression_ok,
142        "baseline_vs_idea": {"baseline_R": 1.0, "idea_R": safe_R,
143                              "trials": n_trials,
144                              "period2_rate_baseline": baseline_rate,
145                              "period2_rate_idea": idea_rate},
146        "all_mechanism_predictions_confirmed": bool(boundary_ok and scaling_ok and suppression_ok)
147    }
148    with open("results.json", "w") as f:
149        json.dump(result, f, indent=2)
150    print(json.dumps(result, indent=2))
151
152
153if __name__ == "__main__":
154    main()