1import json
 2import numpy as np
 3from db_transport import energy, explicit_step, positivity_dt_bound, dissipation
 4
 5
 6def make_graph(n=18, seed=7):
 7    rng = np.random.default_rng(seed)
 8    a = rng.uniform(.15, 1.0, (n, n)); c = (a + a.T) / 2
 9    c *= (rng.random((n, n)) < .28)
10    c = np.triu(c, 1); c += c.T
11    pi = rng.uniform(.4, 1.6, n); pi /= pi.sum()
12    rho = rng.uniform(.15, 2.0, n); rho *= 1.0 / rho.sum()
13    return rho, pi, c
14
15
16def stability_sweep(rho, pi, c):
17    bound = positivity_dt_bound(rho, pi, c)
18    rows = []
19    for factor in [.25, .5, .99, 1.01, 3., 10., 30., 100., 300., 1000.]:
20        x = explicit_step(rho, pi, c, factor * bound)
21        row = {'factor': factor, 'min_mass': float(x.min()),
22               'mass_error': float(abs(x.sum() - rho.sum()))}
23        row['energy_change'] = (energy(x, pi) - energy(rho, pi)) if np.all(x > 0) else None
24        rows.append(row)
25    first_negative = next((r['factor'] for r in rows if r['min_mass'] < 0), None)
26    return bound, rows, first_negative
27
28
29def dissipation_sweep(rho, pi, c):
30    d = dissipation(rho, pi, c); rows = []
31    for dt in [1e-5, 3e-5, 1e-4, 3e-4, 1e-3]:
32        x = explicit_step(rho, pi, c, dt)
33        measured = (energy(rho, pi) - energy(x, pi)) / dt
34        rows.append({'dt': dt, 'measured': measured, 'predicted': d,
35                     'relative_error': abs(measured-d)/d})
36    return rows
37
38
39def conductance_sweep(rho, pi, c):
40    base = dissipation(rho, pi, c); rows = []
41    for scale in [.25, .5, 1., 2., 4.]:
42        actual = dissipation(rho, pi, c * scale)
43        rows.append({'scale': scale, 'dissipation': actual, 'predicted': base * scale,
44                     'relative_error': abs(actual-base*scale)/actual})
45    return rows
46
47
48def classification_utility(seed=11):
49    rng = np.random.default_rng(seed); n = 80
50    y = np.repeat(np.arange(2), n//2)
51    c = np.full((n, n), .12/(n-1)); same = y[:, None] == y[None, :]
52    c[same] = .88/(same.sum(1)[0]-1); np.fill_diagonal(c, 0); c = (c+c.T)/2
53    scores = np.full((n, 2), .15); scores[np.arange(n), y] = .85
54    pi = np.full((n, 2), 1/n); rho = scores / scores.sum(0, keepdims=True)
55    # Baseline is an unconstrained residual with a deliberately large multiplier.
56    x = scores.copy(); baseline_min = 1e9
57    # Transport uses the explicit positivity bound, as required by the idea.
58    dt_transport = .9 * positivity_dt_bound(rho, pi, c)
59    tr = rho.copy(); transport_min = 1e9
60    for _ in range(20):
61        x = x + 2.0 * (c @ x - x)
62        tr = explicit_step(tr, pi, c, dt_transport)
63        baseline_min = min(baseline_min, float(x.min())); transport_min = min(transport_min, float(tr.min()))
64    return {'baseline_accuracy': float((np.argmax(x, 1) == y).mean()),
65            'transport_accuracy': float((np.argmax(tr / pi, 1) == y).mean()),
66            'baseline_min_activation': baseline_min, 'transport_min_mass': transport_min,
67            'transport_dt_over_bound': .9}
68
69
70def main():
71    rho, pi, c = make_graph(); bound, stability, first_negative = stability_sweep(rho, pi, c)
72    diss = dissipation_sweep(rho, pi, c); conduct = conductance_sweep(rho, pi, c)
73    result = {'seed': 7, 'n': len(rho), 'dt_bound': bound,
74      'initial_energy': energy(rho, pi),
75      'stability_prediction': 'nonnegative for dt <= bound; sufficiently larger dt can fail',
76      'stability': stability, 'first_negative_factor_tested': first_negative,
77      'dissipation_prediction': '-dF/dt equals edge dissipation with O(dt) Euler error',
78      'dissipation': diss,
79      'conductance_prediction': 'dissipation scales linearly with global conductance',
80      'conductance': conduct, 'classification': classification_utility(),
81      'max_mass_error': max(r['mass_error'] for r in stability),
82      'max_conductance_relative_error': max(r['relative_error'] for r in conduct),
83      'smallest_dt_dissipation_relative_error': diss[0]['relative_error']}
84    print(json.dumps(result, indent=2))
85
86if __name__ == '__main__': main()