1import json, random
  2from dataclasses import dataclass
  3import numpy as np
  4
  5@dataclass
  6class LayerReport:
  7    lower: np.ndarray
  8    upper: np.ndarray
  9    active: int
 10    inactive: int
 11    unstable: int
 12    retained: int
 13    pruned_contribution: int
 14
 15
 16def affine_interval(W, b, lo, hi):
 17    W, b = np.asarray(W, float), np.asarray(b, float)
 18    lo, hi = np.asarray(lo, float), np.asarray(hi, float)
 19    lower = np.maximum(W, 0) @ lo + np.minimum(W, 0) @ hi + b
 20    upper = np.maximum(W, 0) @ hi + np.minimum(W, 0) @ lo + b
 21    return lower, upper
 22
 23
 24def interval_propagate(hidden_weights, hidden_biases, xlo, xhi):
 25    """Interval bounds for hidden preactivations and post-ReLU variables."""
 26    lo, hi = np.asarray(xlo, float), np.asarray(xhi, float)
 27    out = []
 28    for W, b in zip(hidden_weights, hidden_biases):
 29        L, U = affine_interval(W, b, lo, hi)
 30        zlo, zhi = np.maximum(0, L), np.maximum(0, U)
 31        out.append((L, U, zlo, zhi))
 32        lo, hi = zlo, zhi
 33    return out
 34
 35
 36def prune_report(hidden_weights, hidden_biases, xlo, xhi, tau=0.0):
 37    """Classify hidden neurons; only unstable neurons need SDP activation variables."""
 38    bounds = interval_propagate(hidden_weights, hidden_biases, xlo, xhi)
 39    reports = []
 40    for l, (L, U, _, _) in enumerate(bounds):
 41        # Make categories disjoint at the degenerate boundary.
 42        active_mask = L >= 0
 43        inactive_mask = (~active_mask) & (U <= 0)
 44        unstable_mask = ~(active_mask | inactive_mask)
 45        if l + 1 < len(hidden_weights):
 46            outgoing = np.max(np.abs(hidden_weights[l + 1]), axis=0)
 47        else:
 48            outgoing = np.ones(len(L))
 49        magnitude = np.maximum(np.abs(L), np.abs(U))
 50        contribution = outgoing * magnitude
 51        keep = unstable_mask & (contribution > tau)
 52        reports.append(LayerReport(
 53            L, U, int(active_mask.sum()), int(inactive_mask.sum()),
 54            int(unstable_mask.sum()), int(keep.sum()),
 55            int((unstable_mask & ~keep).sum())))
 56    return bounds, reports
 57
 58
 59def forward(weights, biases, x):
 60    z = np.asarray(x, float)
 61    for k, (W, b) in enumerate(zip(weights, biases)):
 62        z = np.asarray(W) @ z + np.asarray(b)
 63        if k < len(weights) - 1:
 64            z = np.maximum(0, z)
 65    return z
 66
 67
 68def output_interval(weights, biases, xlo, xhi):
 69    hidden_bounds = interval_propagate(weights[:-1], biases[:-1], xlo, xhi)
 70    L, U = hidden_bounds[-1][2], hidden_bounds[-1][3]
 71    return affine_interval(weights[-1], biases[-1], L, U)
 72
 73
 74def make_net(seed=7):
 75    rng = np.random.default_rng(seed)
 76    return ([rng.normal(0, 1.0, (12, 2)), rng.normal(0, .8, (10, 12)), rng.normal(0, .7, (2, 10))],
 77            [rng.normal(0, .9, 12), rng.normal(0, .7, 10), rng.normal(0, .3, 2)])
 78
 79
 80def main():
 81    random.seed(11); np.random.seed(11)
 82    W, b = make_net(); x = np.array([.15, -.2]); eps = .3
 83    xlo, xhi = x - eps, x + eps
 84    hiddenW, hiddenb = W[:-1], b[:-1]
 85    print('CORE_CHECK')
 86    for tau in [0.0, 1e-6, 1e-2, 1e-1, 3e-1, 1.0]:
 87        _, reps = prune_report(hiddenW, hiddenb, xlo, xhi, tau)
 88        print(json.dumps({'tau': tau, 'unstable': sum(r.unstable for r in reps),
 89            'fixed_sign': sum(r.active + r.inactive for r in reps),
 90            'retained': sum(r.retained for r in reps),
 91            'contribution_pruned': sum(r.pruned_contribution for r in reps),
 92            'sdp_variable_proxy': 2 * sum(r.retained for r in reps)}))
 93    bounds, reps = prune_report(hiddenW, hiddenb, xlo, xhi, 0.0)
 94    max_fixed_err = 0.0; samples = []
 95    for xx in np.linspace(xlo[0], xhi[0], 81):
 96        for yy in np.linspace(xlo[1], xhi[1], 81):
 97            v = np.array([xx, yy]); z = v
 98            for l, (Wl, bl) in enumerate(zip(hiddenW, hiddenb)):
 99                a = Wl @ z + bl; L, U = bounds[l][0], bounds[l][1]
100                pred = np.where(L >= 0, a, np.where(U <= 0, 0, np.maximum(0, a)))
101                max_fixed_err = max(max_fixed_err, float(np.max(np.abs(pred - np.maximum(0, a)))))
102                z = np.maximum(0, a)
103            samples.append(forward(W, b, v))
104    samples = np.asarray(samples); outL, outU = output_interval(W, b, xlo, xhi)
105    violations = int(np.sum((samples < outL - 1e-9) | (samples > outU + 1e-9)))
106    print('SOUNDNESS', json.dumps({'grid_points': len(samples), 'max_fixed_substitution_error': max_fixed_err,
107        'output_interval_lower': outL.tolist(), 'output_interval_upper': outU.tolist(),
108        'sample_interval_violations': violations}))
109    print('SUMMARY', json.dumps({'baseline_hidden_sdp_proxy': 2 * sum(len(v) for v in hiddenb),
110        'exact_pruned_hidden_proxy': 2 * sum(r.retained for r in reps),
111        'fixed_sign_neurons': sum(r.active + r.inactive for r in reps)}))
112
113if __name__ == '__main__': main()