1import json
  2import math
  3import os
  4import random
  5import numpy as np
  6from sklearn.metrics import f1_score
  7
  8SEED = 2057
  9rng = np.random.default_rng(SEED)
 10
 11
 12def group_prox(B, lr_lam):
 13    # B shape [outputs, regulators]; one regulator is one regulator-to-state group.
 14    norms = np.linalg.norm(B, axis=0, keepdims=True)
 15    return B * np.maximum(0.0, 1.0 - lr_lam / (norms + 1e-12))
 16
 17
 18def fit_group(X, Y, lam=0.02, epochs=500, lr=None):
 19    # Minimize mean squared multivariate regression loss + group lasso.
 20    X = np.asarray(X); Y = np.asarray(Y)
 21    n, p = X.shape
 22    if lr is None:
 23        lr = 0.8 / (np.linalg.norm(X, 2) ** 2 / max(n, 1) + 1e-8)
 24    W = np.zeros((p, Y.shape[1]))
 25    # Work with rows as regulator groups.
 26    for _ in range(epochs):
 27        G = (X.T @ (X @ W - Y)) / n
 28        W -= lr * G
 29        W = group_prox(W.T, lr * lam).T
 30    return W
 31
 32
 33def make_data(n=5000, d=8, q=12, dt_choices=(0.05, 0.1, 0.2), sigma=0.05,
 34              rho=0.0, seed=0):
 35    rg = np.random.default_rng(seed)
 36    active = np.array([0, 3, 7])
 37    B = np.zeros((q, d))
 38    # correlated regulator features, with a controllable active Gram condition.
 39    z = rg.normal(size=(n, 1))
 40    R = rg.normal(size=(n, q))
 41    for j in active[1:]:
 42        R[:, j] = rho * R[:, 0] + math.sqrt(max(1-rho*rho, 0)) * R[:, j]
 43    B[active] = rg.normal(0, 0.65, size=(len(active), d))
 44    dt = rg.choice(dt_choices, size=n)
 45    # latent increments from the integral equation, noisy observations at endpoints.
 46    clean = (dt[:, None] * (R @ B))
 47    eps0 = rg.normal(0, sigma, size=(n, d))
 48    eps1 = rg.normal(0, sigma, size=(n, d))
 49    y_int = clean + eps1 - eps0
 50    # Derivative target is exactly the same noisy increment divided by dt.
 51    y_fd = y_int / dt[:, None]
 52    return R, dt, y_int, y_fd, B, active
 53
 54
 55def noise_scaling():
 56    # Directly test the two claims: Var(delta x noise)=2 sigma^2 and
 57    # Var((delta x)/dt noise)=2 sigma^2/dt^2.
 58    sigma = 0.07
 59    out = []
 60    for dt in [0.025, 0.05, 0.1, 0.2, 0.4]:
 61        rg = np.random.default_rng(1000 + int(dt * 10000))
 62        e = rg.normal(0, sigma, size=(250000, 1))
 63        f = rg.normal(0, sigma, size=(250000, 1))
 64        inc = (f-e).ravel()
 65        fd = inc / dt
 66        vi, vf = float(np.var(inc)), float(np.var(fd))
 67        out.append({"dt": dt, "integral_var": vi, "integral_pred": 2*sigma*sigma,
 68                    "fd_var": vf, "fd_pred": 2*sigma*sigma/(dt*dt),
 69                    "fd_over_integral": vf/vi, "ratio_pred": 1/(dt*dt)})
 70    return out
 71
 72
 73def support_sweep():
 74    # Prediction: as rho -> 1, lambda_min(active Gram) -> 0 and support recovery degrades.
 75    rows = []
 76    for rho in [0.0, 0.5, 0.8, 0.95, 0.99]:
 77        R, dt, yi, yf, B, active = make_data(n=7000, sigma=0.045, rho=rho, seed=20+int(rho*100))
 78        # Integral formulation uses integrated features dt*r and increments.
 79        Xi = R * dt[:, None]
 80        W = fit_group(Xi, yi, lam=0.018, epochs=350)
 81        # FD formulation uses r and derivative targets; same equal-weight objective.
 82        Wfd = fit_group(R, yf, lam=0.018, epochs=350)
 83        gram = (R[:, active].T @ R[:, active]) / len(R)
 84        eig = float(np.linalg.eigvalsh(gram).min())
 85        def score(w):
 86            pred = np.linalg.norm(w, axis=1) > 0.06
 87            truth = np.zeros(w.shape[0], dtype=bool); truth[active] = True
 88            return float(f1_score(truth, pred)), int(pred.sum())
 89        fi, ni = score(W); ff, nf = score(Wfd)
 90        rows.append({"rho": rho, "lambda_min": eig, "integral_f1": fi,
 91                     "fd_f1": ff, "integral_edges": ni, "fd_edges": nf})
 92    return rows
 93
 94
 95def method_comparison():
 96    # Irregular intervals: integrated loss naturally weights observations by physical
 97    # increment size, while FD noise is amplified on the short intervals.
 98    R, dt, yi, yf, B, active = make_data(n=12000, sigma=0.06,
 99                                          dt_choices=(0.025, 0.05, 0.1, 0.2), rho=0.5, seed=77)
100    Xi = R * dt[:, None]
101    Wi = fit_group(Xi, yi, lam=0.0008, epochs=450)
102    Wf = fit_group(R, yf, lam=0.02, epochs=450)
103    # Evaluate physical increment prediction on fresh noiseless transitions.
104    R2, dt2, y2, yf2, _, _ = make_data(n=4000, sigma=0.0,
105                                         dt_choices=(0.025, 0.05, 0.1, 0.2), rho=0.5, seed=91)
106    true_inc = dt2[:, None] * (R2 @ B)
107    pi = dt2[:, None] * (R2 @ Wi)
108    pf = dt2[:, None] * (R2 @ Wf)
109    def f1(w):
110        pred = np.linalg.norm(w, axis=1) > 0.06
111        truth = np.zeros(w.shape[0], dtype=bool); truth[active] = True
112        return float(f1_score(truth, pred))
113    return {"integral_clean_increment_mse": float(np.mean((pi-true_inc)**2)),
114            "finite_difference_clean_increment_mse": float(np.mean((pf-true_inc)**2)),
115            "integral_support_f1": f1(Wi), "finite_difference_support_f1": f1(Wf),
116            "active_edges": active.tolist()}
117
118
119if __name__ == "__main__":
120    result = {"seed": SEED,
121              "prediction_1_noise_scaling": noise_scaling(),
122              "prediction_2_conditioning_support": support_sweep(),
123              "mini_experiment": method_comparison()}
124    with open("results.json", "w") as f:
125        json.dump(result, f, indent=2)
126    print(json.dumps(result, indent=2))