1import json
  2import numpy as np
  3from scipy.linalg import expm
  4
  5SEED = 2080
  6A = np.diag([0.4, 0.4])
  7C1 = np.array([[1.0, 0.0]])
  8C2 = np.array([[0.0, 1.0]])
  9I = np.eye(2)
 10
 11
 12def obs(A, C):
 13    n = A.shape[0]
 14    return np.vstack([C @ np.linalg.matrix_power(A, k) for k in range(n)])
 15
 16
 17def matrices(gain, tau1=0.5, tau2=0.5):
 18    L1 = np.array([[gain], [0.0]])
 19    L2 = np.array([[0.0], [gain]])
 20    F1 = A - L1 @ C1
 21    F2 = A - L2 @ C2
 22    P1, P2 = expm(F1 * tau1), expm(F2 * tau2)
 23    return F1, F2, P1, P2, P2 @ P1
 24
 25
 26def spectral_radius(M):
 27    return float(np.max(np.abs(np.linalg.eigvals(M))))
 28
 29
 30def exact_boundary():
 31    # For this diagonal complementary system, rho(Phi)=exp((.4-g) tau)
 32    # and the predicted boundary is gain=.8 for equal dwell time.
 33    gains = np.linspace(0.0, 1.2, 241)
 34    values = np.array([spectral_radius(matrices(g)[-1]) for g in gains])
 35    crossing = gains[np.argmin(np.abs(values - 1.0))]
 36    return crossing, gains, values
 37
 38
 39def empirical_decay(gain, tau=0.5, cycles=25):
 40    Phi = matrices(gain, tau, tau)[-1]
 41    e = np.array([1.0, -0.7])
 42    norms = []
 43    for _ in range(cycles):
 44        norms.append(np.linalg.norm(e))
 45        e = Phi @ e
 46    norms = np.maximum(np.asarray(norms), 1e-300)
 47    slope = np.polyfit(np.arange(cycles), np.log(norms), 1)[0]
 48    return slope, norms[-1]
 49
 50
 51def noisy_switched(gain, tau=0.5, noise=0.15, cycles=80, trials=200):
 52    # Discrete exact observer with measurement noise. For each dwell, use
 53    # midpoint-equivalent correction z-C h at the start of the dwell.
 54    rng = np.random.default_rng(SEED + 9)
 55    F1, F2, P1, P2, _ = matrices(gain, tau, tau)
 56    L1 = np.array([[gain], [0.0]])
 57    L2 = np.array([[0.0], [gain]])
 58    # Integral exp(F(t-s)) L ds, computed robustly by augmented exponential.
 59    def noise_map(F, L):
 60        aug = np.zeros((3, 3)); aug[:2, :2] = F; aug[:2, 2:] = L
 61        return expm(aug * tau)[:2, 2:]
 62    Q1, Q2 = noise_map(F1, L1), noise_map(F2, L2)
 63    vals = []
 64    for _ in range(trials):
 65        e = np.array([1.0, -0.7])
 66        for _ in range(cycles):
 67            e = P1 @ e + Q1 @ rng.normal(0, noise, size=(1,))
 68            e = P2 @ e + Q2 @ rng.normal(0, noise, size=(1,))
 69        vals.append(np.dot(e, e))
 70    return float(np.mean(vals))
 71
 72
 73def residual_fusion(gain, tau=0.5, noise=0.15, cycles=80, trials=200):
 74    # Baseline: each observed coordinate is independently corrected, but its
 75    # residual is held at zero during the missing-channel dwell (ordinary
 76    # modality-specific residual fusion, no cycle-aware gain design).
 77    rng = np.random.default_rng(SEED + 9)
 78    F1, F2, P1, P2, _ = matrices(gain, tau, tau)
 79    L1 = np.array([[gain], [0.0]]); L2 = np.array([[0.0], [gain]])
 80    def noise_map(F, L):
 81        aug = np.zeros((3, 3)); aug[:2, :2] = F; aug[:2, 2:] = L
 82        return expm(aug * tau)[:2, 2:]
 83    Q1, Q2 = noise_map(F1, L1), noise_map(F2, L2)
 84    vals = []
 85    for _ in range(trials):
 86        e = np.array([1.0, -0.7])
 87        for _ in range(cycles):
 88            # baseline uses only currently available channel but resets the
 89            # unobserved residual branch; this is equivalent here to no
 90            # complementary correction during that half-cycle.
 91            e = P1 @ e + Q1 @ rng.normal(0, noise, size=(1,))
 92            e = expm(A * tau) @ e + Q2 @ rng.normal(0, noise, size=(1,))
 93        vals.append(np.dot(e, e))
 94    return float(np.mean(vals))
 95
 96
 97def main():
 98    O1, O2 = obs(A, C1), obs(A, C2)
 99    Ostack = np.vstack([O1, O2])
100    crossing, gains, rhos = exact_boundary()
101    # Prediction 1: equal-dwell stability boundary at gain=1.2.
102    boundary = {"predicted": 0.8, "observed_grid": float(crossing),
103                "rho_at_0_7": float(spectral_radius(matrices(.7)[-1])),
104                "rho_at_0_9": float(spectral_radius(matrices(.9)[-1]))}
105    # Prediction 2: log decay slope equals log spectral radius per cycle.
106    decay = []
107    for g in (0.9, 1.0, 1.2):
108        slope, final = empirical_decay(g)
109        predicted = np.log(spectral_radius(matrices(g)[-1]))
110        decay.append({"gain": g, "predicted_log_slope": float(predicted),
111                      "observed_log_slope": float(slope), "final_norm": float(final)})
112    # Prediction 3: stable switched cycle can tolerate an individually expanding
113    # channel; scan dwell ratios while keeping gain=1.2.
114    dwell = []
115    for t1, t2 in ((0.2,0.8),(0.3,0.7),(0.4,0.6),(0.5,0.5),(0.6,0.4),(0.7,0.3),(0.8,0.2)):
116        F1,F2,P1,P2,Phi = matrices(1.2,t1,t2)
117        dwell.append({"tau1":t1,"tau2":t2,"rho_cycle":spectral_radius(Phi),
118                      "rho_channel1":spectral_radius(P1),"rho_channel2":spectral_radius(P2)})
119    gain = 1.2
120    comparison = {"switched_mse": noisy_switched(gain),
121                  "residual_fusion_mse": residual_fusion(gain),
122                  "gain": gain}
123    out = {"observability_rank_channel1": int(np.linalg.matrix_rank(O1)),
124           "observability_rank_channel2": int(np.linalg.matrix_rank(O2)),
125           "observability_rank_stacked": int(np.linalg.matrix_rank(Ostack)),
126           "predictions": {"boundary": boundary, "decay": decay, "dwell_sweep": dwell},
127           "comparison": comparison}
128    with open("results.json", "w") as f: json.dump(out, f, indent=2)
129    print(json.dumps(out, indent=2))
130
131if __name__ == "__main__": main()