1import json
  2import math
  3from dataclasses import dataclass
  4from pathlib import Path
  5import numpy as np
  6
  7
  8@dataclass
  9class EventCell:
 10    theta: float = 1.0
 11    reset: float = 0.0
 12    leak: float = 1.0
 13
 14    def run(self, x, events):
 15        """Process (normalized_time, signed_mass) events in stable time order."""
 16        u = self.leak * float(x)
 17        spikes = []
 18        # Python's stable sort preserves input order for tied timestamps.
 19        for tau, delta, *rest in sorted(events, key=lambda e: e[0]):
 20            u += float(delta)
 21            fired = int(u >= self.theta)
 22            spikes.append(fired)
 23            if fired:
 24                u = self.reset
 25        return int(any(spikes)), u, spikes
 26
 27    def aggregate(self, x, events):
 28        """Baseline: one signed current and one threshold/reset per macro-step."""
 29        u = self.leak * float(x) + sum(float(e[1]) for e in events)
 30        fired = int(u >= self.theta)
 31        if fired:
 32            u = self.reset
 33        return fired, u
 34
 35
 36def order_condition(x, a, b, theta=1.0):
 37    return (x + a - b < theta) and (x + a >= theta)
 38
 39
 40def analytic_fraction(a, b, theta=1.0, x_low=0.0, x_high=1.0):
 41    # x must lie in [theta-a, theta-a+b), intersected with the sampling interval.
 42    lo = max(x_low, theta - a)
 43    hi = min(x_high, theta - a + b)
 44    return max(0.0, hi - lo) / (x_high - x_low)
 45
 46
 47def toy_verification(rng):
 48    theta = 1.0
 49    n = 500_000
 50    x = rng.uniform(0.0, theta, n)
 51
 52    # Prediction 1: unsafe/order-sensitive fraction is the interval length.
 53    a = 0.65
 54    b_grid = np.array([0.05, 0.15, 0.30, 0.50, 0.80])
 55    p_obs, p_pred = [], []
 56    for b in b_grid:
 57        p_obs.append(np.mean((x + a - b < theta) & (x + a >= theta)))
 58        p_pred.append(analytic_fraction(a, b, theta))
 59
 60    # Prediction 2: for fixed b, the fraction grows linearly with a until
 61    # clipping at the available x interval.
 62    b = 0.30
 63    a_grid = np.array([0.10, 0.30, 0.50, 0.70, 0.90, 1.10])
 64    pa_obs, pa_pred = [], []
 65    for aa in a_grid:
 66        pa_obs.append(np.mean((x + aa - b < theta) & (x + aa >= theta)))
 67        pa_pred.append(analytic_fraction(aa, b, theta))
 68
 69    # Prediction 3: conditioned on the strict margin interval, random E/I
 70    # ordering gives a 1/2 expected firing discrepancy.
 71    m = 0.10
 72    x0 = theta - a + b / 2.0  # strictly inside condition for all listed b>=.15
 73    b0 = 0.50
 74    # Explicitly use many random orderings of two events.
 75    trials = 200_000
 76    fires = []
 77    cell = EventCell(theta=theta)
 78    for _ in range(trials):
 79        if rng.random() < 0.5:
 80            ev = [(0.1, a), (0.9, -b0)]
 81        else:
 82            ev = [(0.1, -b0), (0.9, a)]
 83        fires.append(cell.run(x0, ev)[0])
 84    random_order_firing = float(np.mean(fires))
 85    # E-first is one, I-first is zero; baseline is always zero here.
 86    discrepancy_from_half = abs(random_order_firing - 0.5)
 87
 88    return {
 89        "prediction_1_interval_fraction": {
 90            "parameter": "b", "a": a,
 91            "rows": [{"b": float(bb), "predicted": float(pp), "observed": float(po),
 92                      "abs_error": float(abs(pp-po))}
 93                     for bb, pp, po in zip(b_grid, p_pred, p_obs)]
 94        },
 95        "prediction_2_interval_scaling": {
 96            "parameter": "a", "b": b,
 97            "rows": [{"a": float(aa), "predicted": float(pp), "observed": float(po),
 98                      "abs_error": float(abs(pp-po))}
 99                     for aa, pp, po in zip(a_grid, pa_pred, pa_obs)]
100        },
101        "prediction_3_random_order_half": {
102            "predicted": 0.5, "observed": random_order_firing,
103            "abs_error": discrepancy_from_half, "x": x0, "a": a, "b": b0,
104            "margin": min(theta-(x0+a-b0), x0+a-theta)
105        }
106    }
107
108
109def mini_experiment(rng):
110    # Each example has one excitatory and one inhibitory arrival. The reference
111    # is causal processing; the baseline collapses both into one signed pulse.
112    cell = EventCell(theta=1.0, reset=0.0)
113    n = 100_000
114    x = rng.uniform(0.0, 1.0, n)
115    a = rng.uniform(0.05, 0.95, n)
116    b = rng.uniform(0.05, 0.95, n)
117    # Generate both causal orderings, equally often.
118    order = rng.integers(0, 2, n)
119    micro = np.empty(n, dtype=np.int8)
120    base = np.empty(n, dtype=np.int8)
121    unsafe = np.zeros(n, dtype=bool)
122    margins = np.full(n, np.nan)
123    for i in range(n):
124        if order[i] == 0:  # E first
125            events = [(0.1, a[i]), (0.9, -b[i])]
126        else:
127            events = [(0.1, -b[i]), (0.9, a[i])]
128        micro[i] = cell.run(x[i], events)[0]
129        base[i] = cell.aggregate(x[i], events)[0]
130        unsafe[i] = order_condition(x[i], a[i], b[i])
131        if unsafe[i]:
132            margins[i] = min(1.0-(x[i]+a[i]-b[i]), x[i]+a[i]-1.0)
133
134    # In this setup aggregate is exactly the net-current model; compare it to
135    # the causal reference and expose the predicted 1/2 split on unsafe cases.
136    disagreement = float(np.mean(micro != base))
137    unsafe_rate = float(np.mean(unsafe))
138    unsafe_disagreement = float(np.mean((micro != base)[unsafe])) if unsafe.any() else 0.0
139    micro_rate = float(np.mean(micro))
140    base_rate = float(np.mean(base))
141    return {
142        "n": n,
143        "unsafe_fraction": unsafe_rate,
144        "mean_margin_on_unsafe": float(np.nanmean(margins)),
145        "baseline_spike_rate": base_rate,
146        "micro_event_spike_rate": micro_rate,
147        "overall_disagreement": disagreement,
148        "disagreement_conditioned_on_unsafe": unsafe_disagreement,
149        "predicted_conditioned_disagreement": 0.5,
150        "order_balance": float(np.mean(order == 0))
151    }
152
153
154def main():
155    rng = np.random.default_rng(2411)
156    toy = toy_verification(rng)
157    mini = mini_experiment(rng)
158    result = {"seed": 2411, "toy_verification": toy, "mini_experiment": mini}
159    Path("results.json").write_text(json.dumps(result, indent=2))
160    print(json.dumps(result, indent=2))
161
162
163if __name__ == "__main__":
164    main()