1import json, math
 2import numpy as np
 3from scipy.special import roots_hermitenorm, eval_hermitenorm
 4from hermite_bottleneck import multi_indices, heat_identity_check
 5
 6
 7def one_dim_coeffs(var, maxdeg, nodes=80):
 8    # E_{N(0,var)}[He_k(X)/sqrt(k!)] by Gaussian quadrature.
 9    x, wt = roots_hermitenorm(nodes)
10    wt = wt / math.sqrt(2 * math.pi)
11    y = math.sqrt(var) * x
12    return np.array([np.sum(wt * eval_hermitenorm(k, y)) / math.sqrt(math.factorial(k))
13                     for k in range(maxdeg + 1)])
14
15
16def exact_report(vars_, maxdeg=12):
17    d = len(vars_)
18    one = [one_dim_coeffs(v, maxdeg) for v in vars_]
19    allidx = multi_indices(d, maxdeg)
20    coeff = np.array([np.prod([one[j][a[j]] for j in range(d)]) for a in allidx])
21    # Parseval norm of the density ratio, computed independently in closed form.
22    S = np.diag(vars_)
23    norm2 = np.linalg.det(S) ** -1 * np.linalg.det(2*np.linalg.inv(S)-np.eye(d)) ** -.5
24    rows = []
25    for N in range(maxdeg + 1):
26        keep = np.array([sum(a) <= N for a in allidx])
27        err = math.sqrt(max(0., norm2 - float(np.sum(coeff[keep]**2)))) / math.sqrt(norm2)
28        rows.append((N, int(keep.sum()), err))
29    q = max(abs(v - 1) for v in vars_)
30    return q, norm2, rows
31
32
33def run():
34    rng = np.random.default_rng(11)
35    # The semigroup formula is checked independently by Monte Carlo.
36    heat = heat_identity_check(rng, n=12, nz=12)
37    cases = [[1.04, .97, 1.01, .99], [1.20, .82, 1.08, .94],
38             [1.55, .62, 1.18, .88]]
39    reports = []
40    for v in cases:
41        q, norm2, rows = exact_report(v)
42        reports.append({'variances': v, 'q': q, 'norm2': norm2,
43                        'curve': [{'degree': n, 'coefficients': k, 'relative_tail': e}
44                                  for n,k,e in rows]})
45    # Fit the observed geometric envelope on degrees 2..8; report ratios rather than
46    # asserting a universal constant (the statement allows C(Sigma,u,d)).
47    ratios = []
48    for r in reports:
49        q = r['q']
50        for row in r['curve']:
51            if row['relative_tail'] > 1e-12 and q > 0:
52                ratios.append(row['relative_tail'] / q**((row['degree']+1)/2))
53    # Practical degree choice using a conservative empirical C from all cases.
54    C = max(ratios)
55    target = .08
56    for r in reports:
57        q = r['q']
58        chosen = next((x['degree'] for x in r['curve'] if C*q**((x['degree']+1)/2) <= target), 12)
59        r['chosen_degree'] = chosen
60        r['chosen_coefficients'] = next(x['coefficients'] for x in r['curve'] if x['degree']==chosen)
61        r['chosen_tail'] = next(x['relative_tail'] for x in r['curve'] if x['degree']==chosen)
62        fixed = next(x for x in r['curve'] if x['degree']==4)
63        r['fixed_degree_4_tail'] = fixed['relative_tail']
64        r['fixed_degree_4_coefficients'] = fixed['coefficients']
65    return {'heat_relative_error': heat, 'empirical_envelope_C': C,
66            'target': target, 'reports': reports}
67
68if __name__ == '__main__':
69    print(json.dumps(run(), indent=2))