1"""Split-granularity extrapolation weights (generalized jackknife / Richardson extrapolation).23For split multipliers s_1 < ... < s_r (s_j sub-blocks of size m/s_j per worker; s_1 = 1 is the full block),4the weights solve the Vandermonde system5 sum_j w_j = 1, sum_j w_j s_j^a = 0, a = 1, ..., r-1,6so that a bias expansion b_1/n + ... + b_{r-1}/n^{r-1} evaluated at n = m/s_j cancels through order m^{-(r-1)}.7"""8from __future__ import annotations9from fractions import Fraction10import numpy as np1112DEFAULT_SPLITS = {1: (1,), 2: (1, 2), 3: (1, 2, 4), 4: (1, 2, 4, 8), 5: (1, 2, 4, 8, 16)}131415def vandermonde_weights_exact(splits):16 """Exact rational weights (list of Fraction) for the split multipliers `splits`."""17 s = [Fraction(int(x)) for x in splits]18 r = len(s)19 # Solve by Gaussian elimination on exact rationals.20 A = [[sj ** a for sj in s] for a in range(r)]21 b = [Fraction(1)] + [Fraction(0)] * (r - 1)22 M = [row[:] + [bi] for row, bi in zip(A, b)]23 for col in range(r):24 piv = next(i for i in range(col, r) if M[i][col] != 0)25 M[col], M[piv] = M[piv], M[col]26 pv = M[col][col]27 M[col] = [x / pv for x in M[col]]28 for i in range(r):29 if i != col and M[i][col] != 0:30 f = M[i][col]31 M[i] = [x - f * y for x, y in zip(M[i], M[col])]32 return [M[i][r] for i in range(r)]333435def vandermonde_weights(splits):36 """Float64 weights."""37 return np.array([float(w) for w in vandermonde_weights_exact(splits)])383940def leading_bias_coefficient(splits):41 """beta_r = sum_j w_j s_j^r: multiplies b_r m^{-r} in the residual bias of J_r (exact rational)."""42 w = vandermonde_weights_exact(splits); r = len(splits)43 return sum(wj * Fraction(int(sj)) ** r for wj, sj in zip(w, splits))444546def amplification_constant(splits):47 """A_r = sum_j |w_j| s_j: multiplies sqrt(C)/m in the L2 bound of the corrected remainder (Lemma 3)."""48 w = vandermonde_weights_exact(splits)49 return sum(abs(wj) * Fraction(int(sj)) for wj, sj in zip(w, splits))505152def weight_l1(splits):53 return sum(abs(w) for w in vandermonde_weights_exact(splits))545556def vandermonde_condition_number(splits):57 s = np.asarray(splits, dtype=float); r = len(s)58 A = np.vstack([s ** a for a in range(r)])59 return float(np.linalg.cond(A))606162def combine_levels(level_means, weights):63 """J = sum_j w_j T_{s_j}. level_means: array (r, ..., p) of per-level averages; weights: (r,)."""64 w = np.asarray(weights, dtype=float)65 return np.tensordot(w, np.asarray(level_means), axes=(0, 0))666768def level_averages(estimates_by_level):69 """estimates_by_level[j]: array (K, s_j, p) of sub-block estimates at level j -> (r, K, p) level means."""70 return np.stack([np.mean(E, axis=1) for E in estimates_by_level], axis=0)71