1import json
 2import numpy as np
 3from adaptive_cur import cur_factors, residual_leverage_refresh, rel_error
 4
 5
 6def one(seed):
 7    rng = np.random.default_rng(seed)
 8    m, n, r, k = 64, 80, 8, 12
 9    U, _ = np.linalg.qr(rng.standard_normal((m, r)))
10    V, _ = np.linalg.qr(rng.standard_normal((n, r)))
11    W0 = (U * np.linspace(3.0, .5, r)) @ V.T
12    rows, cols = np.arange(k), np.arange(k)
13    stale = cur_factors(W0, rows, cols)
14    W1 = W0.copy()
15    # New structure is deliberately outside the original index sets.
16    ir = np.arange(40, 52); ic = np.arange(55, 67)
17    W1[np.ix_(ir, ic)] += .8 * rng.standard_normal((len(ir), len(ic)))
18    stale_e = rel_error(W1, stale)
19    af, ar, ac, _ = residual_leverage_refresh(W1, rows, cols, probes=16, retain=.25,
20                                               rng=np.random.default_rng(seed + 100))
21    adaptive_e = rel_error(W1, af)
22    rf = cur_factors(W1, rng.choice(m, k, replace=False), rng.choice(n, k, replace=False))
23    random_e = rel_error(W1, rf)
24    return stale_e, adaptive_e, random_e
25
26vals = np.array([one(s) for s in range(20)])
27print(json.dumps({
28    'seeds': 20,
29    'mean_relative_errors': dict(stale=float(vals[:,0].mean()), adaptive=float(vals[:,1].mean()), random=float(vals[:,2].mean())),
30    'median_relative_errors': dict(stale=float(np.median(vals[:,0])), adaptive=float(np.median(vals[:,1])), random=float(np.median(vals[:,2]))),
31    'adaptive_beats_stale': int(np.sum(vals[:,1] < vals[:,0])),
32    'adaptive_beats_random': int(np.sum(vals[:,1] < vals[:,2])),
33    'per_seed': vals.tolist()
34}, indent=2))