1import json
 2import numpy as np
 3from scipy.optimize import brentq
 4import interval_equilibrium as ie
 5
 6
 7def roots(a, lo=-2.0, hi=2.0):
 8    grid = np.linspace(lo, hi, 40001)
 9    vals = a * np.tanh(grid) - grid
10    ans = []
11    for x, y, fx, fy in zip(grid[:-1], grid[1:], vals[:-1], vals[1:]):
12        if fx == 0 or fx * fy < 0:
13            z = x if fx == 0 else brentq(lambda t: a*np.tanh(t)-t, x, y)
14            if not ans or abs(z - ans[-1]) > 1e-6:
15                ans.append(float(z))
16    return ans
17
18
19def inclusion_check(a, lo, hi, p0, prad, n=20000, seed=7):
20    rng = np.random.default_rng(seed)
21    lo, hi = np.asarray(lo, float), np.asarray(hi, float)
22    out = ie.krawczyk(lo, hi, a, p0, prad)
23    jl, jh = ie.jac_interval(lo, hi, a)
24    flo = a*np.tanh(lo) + (p0-prad) - hi
25    fhi = a*np.tanh(hi) + (p0+prad) - lo
26    ok = True
27    for _ in range(n):
28        z = rng.uniform(lo, hi)
29        p = rng.uniform(p0-prad, p0+prad)
30        f = a*np.tanh(z) + p - z
31        j = a/np.cosh(z)**2 - 1
32        ok = ok and np.all(f >= flo-1e-12) and np.all(f <= fhi+1e-12)
33        ok = ok and np.all(j >= jl-1e-12) and np.all(j <= jh+1e-12)
34    return {'samples': n, 'included': bool(ok),
35            'q': None if out is None else float(out[2])}
36
37
38def main():
39    contraction = []
40    for a in [.7, .9, .99, 1.0, 1.01, 1.2]:
41        out = ie.krawczyk(np.array([-.2, -.2]), np.array([.2, .2]), a)
42        contraction.append({'gain': a, 'q': None if out is None else float(out[2])})
43
44    domains = []
45    for a in [.7, 1.2, 1.6]:
46        r = ie.classify(a, lo=(-2., -2.), hi=(2., 2.), budget=10000)
47        domains.append({'gain': a, 'scalar_roots': len(roots(a)),
48                        'status': r.status, 'certified': r.certified,
49                        'excluded': r.excluded, 'inconclusive': r.inconclusive,
50                        'visited': r.visited})
51
52    uncertainty = []
53    for prad in [0., .1, .2, .3, .5]:
54        r = ie.classify(.7, prad=prad, lo=(-1., -1.), hi=(1., 1.), budget=10000)
55        uncertainty.append({'prad': prad, 'status': r.status,
56                            'visited': r.visited})
57
58    baseline = []
59    for a in [.7, .9, 1.2]:
60        z, steps = ie.fixed_point(a, np.array([.05, .05]), steps=300)
61        baseline.append({'gain': a, 'steps': steps,
62                         'residual': float(np.max(np.abs(a*np.tanh(z)+.05-z)))})
63
64    result = {'contraction_q': contraction, 'global_domains': domains,
65              'parameter_uncertainty': uncertainty, 'fixed_point_baseline': baseline,
66              'interval_inclusion': inclusion_check(.7, [-.4, -.3], [.6, .5], .1, .2)}
67    with open('mvp_results.json', 'w') as f:
68        json.dump(result, f, indent=2)
69    print(json.dumps(result, indent=2))
70
71
72if __name__ == '__main__':
73    main()