1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2844
  6rng = np.random.default_rng(SEED)
  7OUT = Path('results.json')
  8
  9
 10def relu(x):
 11    return np.maximum(x, 0.0)
 12
 13
 14def rollout(A, B, b, x0, U):
 15    xs = [x0.copy()]
 16    zs = []
 17    x = x0.copy()
 18    for u in U:
 19        z = A @ x + B @ u + b
 20        zs.append(z.copy())
 21        x = relu(z)
 22        xs.append(x.copy())
 23    return np.asarray(xs), np.asarray(zs)
 24
 25
 26def observability(A, C, T):
 27    blocks = []
 28    P = np.eye(A.shape[0])
 29    for _ in range(T + 1):
 30        blocks.append(C @ P)
 31        P = P @ A
 32    return np.vstack(blocks)
 33
 34
 35def rank_tol(M):
 36    s = np.linalg.svd(M, compute_uv=False)
 37    return int(np.sum(s > max(M.shape) * np.finfo(float).eps * s[0])) if s[0] else 0
 38
 39
 40def prediction_positive_regime():
 41    # Positive A, input, bias, and initial state imply every preactivation is positive.
 42    n, p, steps = 8, 3, 30
 43    A = 0.82 * np.eye(n) + 0.06 * np.ones((n, n)) / n
 44    B = 0.12 * np.ones((n, p))
 45    b = np.full(n, 0.03)
 46    x0 = np.full(n, 0.4)
 47    U = rng.uniform(0.0, 1.0, size=(steps, p))
 48    xr, z = rollout(A, B, b, x0, U)
 49    # Explicit affine recurrence, which must equal ReLU recurrence.
 50    xa = [x0.copy()]
 51    x = x0.copy()
 52    for u in U:
 53        x = A @ x + B @ u + b
 54        xa.append(x.copy())
 55    xa = np.asarray(xa)
 56    return {
 57        'predicted_max_difference': 0.0,
 58        'observed_max_difference': float(np.max(np.abs(xr - xa))),
 59        'negative_preactivation_fraction': float(np.mean(z < 0)),
 60        'pass': bool(np.max(np.abs(xr - xa)) < 1e-12 and np.all(z >= 0)),
 61    }
 62
 63
 64def prediction_rank_transition():
 65    # A diagonal system with distinct eigenvalues and one coordinate observed.
 66    # For C=[1,0,...], T+1 observations span k coordinates only when T+1 >= k.
 67    n = 6
 68    eig = np.array([0.92, 0.78, 0.64, 0.50, 0.36, 0.22])
 69    A = np.diag(eig)
 70    C = np.ones((1, n))  # all modes visible; Vandermonde rows become full rank at T=n-1
 71    rows = []
 72    for T in range(0, n + 3):
 73        Q = observability(A, C, T)
 74        s = np.linalg.svd(Q, compute_uv=False)
 75        rows.append({'T': T, 'rank': rank_tol(Q), 'sigma_min': float(s[-1])})
 76    observed_transition = next(r['T'] for r in rows if r['rank'] == n)
 77    predicted_transition = n - 1
 78    return {
 79        'predicted_full_rank_first_T': predicted_transition,
 80        'observed_full_rank_first_T': observed_transition,
 81        'sweep': rows,
 82        'pass': observed_transition == predicted_transition,
 83    }
 84
 85
 86def prediction_noise_scaling():
 87    # In the linear positive regime y=Qx0+noise. Least-squares error is governed by 1/sigma_min(Q).
 88    n = 5
 89    eig = np.array([0.95, 0.80, 0.65, 0.50, 0.35])
 90    A = np.diag(eig)
 91    C = np.ones((1, n))
 92    noise = 1e-3
 93    x_trials = rng.normal(size=(1500, n))
 94    result = []
 95    for T in range(n - 1, 13):
 96        Q = observability(A, C, T)
 97        s = np.linalg.svd(Q, compute_uv=False)
 98        # Keep initial conditions positive while preserving a linear identification problem.
 99        x = 0.5 + 0.15 * x_trials
100        y = x @ Q.T + rng.normal(scale=noise, size=(len(x), Q.shape[0]))
101        est = np.linalg.lstsq(Q, y.T, rcond=None)[0].T
102        rmse = float(np.sqrt(np.mean((est - x) ** 2)))
103        result.append({'T': T, 'sigma_min': float(s[-1]), 'rmse': rmse,
104                       'rmse_times_sigma': rmse * float(s[-1])})
105    # Compare the monotonic qualitative prediction and a fitted log slope.
106    sig = np.array([r['sigma_min'] for r in result])
107    err = np.array([r['rmse'] for r in result])
108    slope = float(np.polyfit(np.log(sig), np.log(err), 1)[0])
109    monotonic = bool(np.all(np.diff(sig) > 0) and np.all(np.diff(err) < 0))
110    return {
111        'predicted': 'sigma_min increases with T and reconstruction RMSE decreases approximately as 1/sigma_min',
112        'observed_log_error_vs_log_sigma_slope': slope,
113        'sweep': result,
114        'pass': monotonic and slope < -0.5,
115    }
116
117
118def positivity_penalty_sweep():
119    # A fourth practical check: shifting bias upward removes masks and penalty.
120    n, p, steps = 8, 2, 20
121    A = 0.65 * np.eye(n) + 0.03 * np.ones((n, n))
122    B = 0.08 * np.ones((n, p))
123    x0 = np.full(n, 0.2)
124    U = rng.uniform(-1, 1, size=(steps, p))
125    out = []
126    for bias in [-0.20, -0.10, 0.0, 0.10, 0.20]:
127        xs, z = rollout(A, B, np.full(n, bias), x0, U)
128        neg = np.maximum(-z, 0)
129        out.append({'bias': bias, 'negative_fraction': float(np.mean(z < 0)),
130                    'positivity_penalty': float(np.sum(neg * neg))})
131    return out
132
133
134def baseline_vs_positive():
135    # Same observable system dimensions and inputs; compare signed ReLU dynamics
136    # against a nonnegative system. The latter is the proposed regime.
137    n, m, steps = 8, 4, 12
138    C = np.eye(n)[:m]
139    x0 = np.full(n, 0.3)
140    U = rng.uniform(-0.05, 0.05, size=(steps, 2))
141    cases = {
142        'baseline_signed_relu': (0.55 * rng.normal(size=(n, n)) / np.sqrt(n),
143                                 0.05 * rng.normal(size=(n, 2)), np.zeros(n)),
144        'positive_regime_relu': (0.55 * (0.5 + rng.random((n, n))) / np.sqrt(n),
145                                 0.05 * (0.5 + rng.random((n, 2))), np.full(n, 0.12)),
146    }
147    out = {}
148    for name, (A, B, b) in cases.items():
149        xs, zs = rollout(A, B, b, x0, U)
150        # Finite differences test the linear perturbation claim along the realized path.
151        d = rng.normal(size=n); d *= 1e-6
152        xp, _ = rollout(A, B, b, x0 + d, U)
153        actual = (xp[1:] - xs[1:])
154        predicted = np.asarray([np.linalg.matrix_power(A, t + 1) @ d for t in range(steps)])
155        rel_err = np.linalg.norm(actual - predicted) / max(np.linalg.norm(actual), 1e-15)
156        Q = observability(A, C, steps - 1)
157        out[name] = {
158            'negative_preactivation_fraction': float(np.mean(zs < 0)),
159            'linearized_perturbation_relative_error': float(rel_err),
160            'sigma_min_Q': float(np.linalg.svd(Q, compute_uv=False)[-1]),
161            'observability_rank': rank_tol(Q),
162        }
163    return out
164
165
166def main():
167    results = {
168        'seed': SEED,
169        'positive_regime_exactness': prediction_positive_regime(),
170        'observability_rank_transition': prediction_rank_transition(),
171        'noise_reconstruction_scaling': prediction_noise_scaling(),
172        'positivity_bias_sweep': positivity_penalty_sweep(),
173        'baseline_vs_positive': baseline_vs_positive(),
174    }
175    results['mechanism_pass'] = all([
176        results['positive_regime_exactness']['pass'],
177        results['observability_rank_transition']['pass'],
178        results['noise_reconstruction_scaling']['pass'],
179    ])
180    OUT.write_text(json.dumps(results, indent=2))
181    print(json.dumps(results, indent=2))
182
183if __name__ == '__main__':
184    main()