1import json
  2import math
  3import numpy as np
  4from pathlib import Path
  5
  6# Adjoint-weak fractional residual MVP.
  7# The Grünwald matrix is used only to make the discrete integration-by-parts
  8# identity explicit; weak measurements use its exact transpose.
  9
 10def grunwald_matrix(n, alpha, T=1.0):
 11    h = T / (n - 1)
 12    # left-sided GL derivative: (D f)_i = h^-alpha sum_{k<=i} (-1)^k C(alpha,k) f_{i-k}
 13    w = np.empty(n)
 14    w[0] = 1.0
 15    for k in range(1, n):
 16        w[k] = w[k-1] * (-(alpha-k+1)/k)
 17    A = np.zeros((n, n))
 18    for i in range(n):
 19        A[i, :i+1] = h**(-alpha) * w[:i+1][::-1]
 20    return A
 21
 22def trapezoid_weights(n, T=1.0):
 23    h = T/(n-1)
 24    q = np.full(n, h)
 25    q[[0,-1]] = h/2
 26    return q
 27
 28def gaussian_phi(t, center=.62, width=.16):
 29    return np.exp(-0.5*((t-center)/width)**2)
 30
 31def one_check(n, alpha, sigma, rng):
 32    hq = trapezoid_weights(n)
 33    A = grunwald_matrix(n, alpha)
 34    phi = gaussian_phi(np.linspace(0,1,n))
 35    # Integral phi^T W D f = f^T D^T W phi.
 36    lhs = phi @ (hq * (A @ rng.normal(size=n)))
 37    qweak = A.T @ (hq * phi)
 38    rhs = qweak @ rng.normal(size=n) # independent draw, used only for norm check below
 39    # exact coefficient vector acting on sampled f
 40    q = qweak
 41    return float(np.dot(q,q)), float(np.dot(hq*phi, hq*phi)), A, q, phi, hq
 42
 43def noisy_variance_sweep():
 44    rng = np.random.default_rng(11)
 45    n, alpha = 257, .8
 46    A = grunwald_matrix(n, alpha)
 47    t = np.linspace(0,1,n); hq = trapezoid_weights(n)
 48    phi = gaussian_phi(t)
 49    qweak = A.T @ (hq*phi)
 50    # A pointwise derivative at an interior location, versus one weak scalar.
 51    i = n//2
 52    qstrong = A[i].copy()
 53    rows=[]
 54    for sigma in [.01,.03,.1,.3]:
 55        vals_w=[]; vals_s=[]
 56        for _ in range(5000):
 57            e=rng.normal(0,sigma,n)
 58            vals_w.append(qweak@e); vals_s.append(qstrong@e)
 59        vw=np.var(vals_w, ddof=1); vs=np.var(vals_s, ddof=1)
 60        rows.append({'sigma':sigma,'weak_var':vw,'strong_var':vs,
 61                     'weak_over_sigma2':vw/sigma**2,'strong_over_sigma2':vs/sigma**2})
 62    return rows
 63
 64def resolution_sweep():
 65    alpha=.8; rng=np.random.default_rng(12); out=[]
 66    for n in [65,129,257,513]:
 67        t=np.linspace(0,1,n); hq=trapezoid_weights(n)
 68        A=grunwald_matrix(n,alpha); phi=gaussian_phi(t)
 69        qw=A.T@(hq*phi); qs=A[n//2]
 70        out.append({'n':n,'weak_norm2':float(qw@qw),'strong_norm2':float(qs@qs)})
 71    # consecutive ratios are the directly testable scaling predictions
 72    for i in range(1,len(out)):
 73        out[i]['weak_ratio_prev']=out[i]['weak_norm2']/out[i-1]['weak_norm2']
 74        out[i]['strong_ratio_prev']=out[i]['strong_norm2']/out[i-1]['strong_norm2']
 75    return out
 76
 77def coefficient_fit(seed, sigma, n=129, alpha=.8, trials=200):
 78    rng=np.random.default_rng(seed); t=np.linspace(0,1,n); hq=trapezoid_weights(n)
 79    A=grunwald_matrix(n,alpha); phi=gaussian_phi(t)
 80    # Synthetic relation D^alpha u = c u on a fixed smooth trajectory. We estimate
 81    # c from noisy sampled u; weak form transfers D to phi.
 82    u=np.exp(-.8*t)*(1+.2*np.sin(2*np.pi*t)); ctrue=.7
 83    # use a generated target equal to the discrete operator, avoiding continuum mismatch
 84    target=A@u
 85    qweak=A.T@(hq*phi); qid=hq*phi
 86    # Compare noise-induced error to each estimator's noiseless value; this
 87    # isolates robustness without claiming the arbitrary trajectory obeys a
 88    # pointwise constant-coefficient fractional ODE.
 89    weak0=(qweak@u)/(qid@u)
 90    inds=np.arange(n//2-3,n//2+4)
 91    strong0=np.sum(u[inds]*(A[inds]@u))/np.sum(u[inds]**2)
 92    csw=[]; css=[]
 93    for _ in range(trials):
 94        un=u+rng.normal(0,sigma,n)
 95        # weak scalar equation: int phi D u = c int phi u
 96        den=qid@un; csw.append((qweak@un)/den)
 97        # strong local equation, average a central neighborhood for fair scalar noise
 98        inds=np.arange(n//2-3,n//2+4)
 99        d=A[inds]@un
100        css.append(np.sum(un[inds]*d)/np.sum(un[inds]**2))
101    return {'sigma':sigma,'weak_rmse':float(np.sqrt(np.mean((np.array(csw)-weak0)**2))),
102            'strong_rmse':float(np.sqrt(np.mean((np.array(css)-strong0)**2))),
103            'weak_noiseless':float(weak0),'strong_noiseless':float(strong0)}
104
105def identity_check():
106    rng=np.random.default_rng(4); n=97; alpha=.65
107    A=grunwald_matrix(n,alpha); q=trapezoid_weights(n); f=rng.normal(size=n); phi=rng.normal(size=n)
108    lhs=phi@(q*(A@f)); rhs=f@(A.T@(q*phi))
109    return {'absolute_error':float(abs(lhs-rhs)), 'relative_error':float(abs(lhs-rhs)/(1+abs(lhs)))}
110
111def main():
112    result={'identity':identity_check(), 'noise_sweep':noisy_variance_sweep(),
113            'resolution_sweep':resolution_sweep(),
114            'coefficient_fit':[coefficient_fit(100+i,s) for i,s in enumerate([.01,.03,.1,.3])]}
115    # Predictions from the mechanism: variance is sigma^2 ||q||^2; weak q norm
116    # scales approximately h^(1/2) for a fixed smooth localized kernel. The
117    # GL point-row norm is dominated by summable near-diagonal weights, so its
118    # squared norm scales as h^(-2 alpha), giving ratio 2^(2 alpha).
119    r=result['resolution_sweep']
120    result['predictions']={
121      'noise_linearity':'weak variance / sigma^2 should be constant',
122      'weak_resolution_ratio_predicted':0.5,
123      'strong_resolution_ratio_predicted':2**(2*.8),
124      'observed_weak_ratios':[x.get('weak_ratio_prev') for x in r[1:]],
125      'observed_strong_ratios':[x.get('strong_ratio_prev') for x in r[1:]],
126      'identity_tolerance':'machine precision'}
127    Path('results.json').write_text(json.dumps(result,indent=2))
128    print(json.dumps(result,indent=2))
129
130if __name__=='__main__': main()