1import json
  2import numpy as np
  3from scipy.optimize import minimize_scalar
  4
  5SEED = 3148
  6rng = np.random.default_rng(SEED)
  7
  8# A bounded activation model makes the infinite-domain construction numerically
  9# concrete while preserving the paper's local asymptotics.
 10XMAX = 8.0
 11RATE = 1.0
 12SCALE = 0.5
 13
 14
 15def sample_x(n=400000):
 16    u = rng.random(n)
 17    # inverse CDF of Exp(1) conditional on x <= XMAX
 18    return -np.log(1.0 - u * (1.0 - np.exp(-RATE * XMAX))) / RATE
 19
 20
 21def density(x):
 22    z = 1.0 - np.exp(-RATE * XMAX)
 23    return RATE * np.exp(-RATE * x) / z
 24
 25
 26def utility(x):
 27    return np.log1p(np.maximum(x, 0.0) / SCALE) / np.log1p(XMAX / SCALE)
 28
 29
 30def utility_prime(x):
 31    return 1.0 / ((SCALE + x) * np.log1p(XMAX / SCALE))
 32
 33
 34def weighted_thresholds(n, grid_n=200001):
 35    grid = np.linspace(0.0, XMAX, grid_n)
 36    g = np.sqrt(density(grid) * utility_prime(grid))
 37    # cumulative trapezoidal integration
 38    cdf = np.concatenate([[0.0], np.cumsum((g[:-1] + g[1:]) * np.diff(grid) * 0.5)])
 39    cdf /= cdf[-1]
 40    targets = np.arange(1, n) / n
 41    internal = np.interp(targets, cdf, grid)
 42    return np.r_[0.0, internal, XMAX]
 43
 44
 45def uniform_thresholds(n):
 46    return np.linspace(0.0, XMAX, n + 1)
 47
 48
 49def left_edge_quantize(x, edges):
 50    # edges has n+1 entries; the last edge is the finite support endpoint.
 51    idx = np.searchsorted(edges, x, side="right") - 1
 52    return edges[np.clip(idx, 0, len(edges) - 2)]
 53
 54
 55def lloyd_max(x, n, iterations=80):
 56    # Standard MSE scalar quantizer, with nearest-centroid reconstruction.
 57    cent = np.quantile(x, (np.arange(n) + 0.5) / n)
 58    for _ in range(iterations):
 59        boundaries = (cent[:-1] + cent[1:]) * 0.5
 60        labels = np.searchsorted(boundaries, x)
 61        new = cent.copy()
 62        for k in range(n):
 63            vals = x[labels == k]
 64            if len(vals):
 65                new[k] = vals.mean()
 66        if np.max(np.abs(new - cent)) < 1e-10:
 67            break
 68        cent = new
 69    return cent
 70
 71
 72def lloyd_quantize(x, cent):
 73    boundaries = (cent[:-1] + cent[1:]) * 0.5
 74    return cent[np.searchsorted(boundaries, x)]
 75
 76
 77def utility_constant():
 78    grid = np.linspace(0.0, XMAX, 400001)
 79    g = np.sqrt(density(grid) * utility_prime(grid))
 80    z = np.trapz(g, grid)
 81    return 0.5 * z * z
 82
 83
 84def evaluate(x, q, ideal_utility=None):
 85    if ideal_utility is None:
 86        ideal_utility = utility(x)
 87    gap = ideal_utility - utility(q)
 88    return {
 89        "utility_gap": float(np.mean(gap)),
 90        "n_gap": None,
 91        "mse": float(np.mean((x - q) ** 2)),
 92        "upward_fraction": float(np.mean(q > x + 1e-12)),
 93        "mean_abs_error": float(np.mean(np.abs(x - q))),
 94    }
 95
 96
 97def main():
 98    x = sample_x()
 99    ideal = utility(x)
100    c_theory = utility_constant()
101    rows = []
102    scaling = []
103    for n in [4, 8, 16, 32, 64]:
104        ew = weighted_thresholds(n)
105        eu = uniform_thresholds(n)
106        qw = left_edge_quantize(x, ew)
107        qu = left_edge_quantize(x, eu)
108        cent = lloyd_max(x, n)
109        ql = lloyd_quantize(x, cent)
110        rw = evaluate(x, qw, ideal)
111        ru = evaluate(x, qu, ideal)
112        rl = evaluate(x, ql, ideal)
113        rw["n_gap"] = n * rw["utility_gap"]
114        ru["n_gap"] = n * ru["utility_gap"]
115        rl["n_gap"] = n * rl["utility_gap"]
116        rows.append({"n": n, "weighted_left": rw, "uniform_left": ru, "mse_lloyd": rl})
117        scaling.append({"n": n, "weighted_n_gap": rw["n_gap"], "uniform_n_gap": ru["n_gap"], "theory_Cw": c_theory})
118
119    # Direct numerical check of the local allocation rule: weighted bins have
120    # nearly equal integrals of sqrt(p Q'), unlike uniform bins.
121    e = weighted_thresholds(32)
122    grid = np.linspace(0, XMAX, 300001)
123    g = np.sqrt(density(grid) * utility_prime(grid))
124    gcdf = np.concatenate([[0], np.cumsum((g[:-1] + g[1:]) * np.diff(grid) / 2)])
125    bin_mass = np.diff(np.interp(e, grid, gcdf))
126    allocation_check = {
127        "weighted_g_mass_cv": float(np.std(bin_mass) / np.mean(bin_mass)),
128        "weighted_g_mass_min": float(bin_mass.min()),
129        "weighted_g_mass_max": float(bin_mass.max()),
130    }
131
132    result = {
133        "seed": SEED,
134        "distribution": "Exp(rate=1) truncated to [0,8]",
135        "utility": "log1p(x/0.5)/log1p(8/0.5)",
136        "theory_Cw": c_theory,
137        "results": rows,
138        "scaling": scaling,
139        "allocation_check": allocation_check,
140        "interpretation": {
141            "claimed_signal": "weighted left-edge utility gap scales approximately as Cw/n and is lower than simple uniform left-edge",
142            "conservative_rule": "weighted and uniform left-edge have zero upward rounding by construction",
143            "mse_control": "Lloyd-Max minimizes squared error, not the monotone left-edge utility gap",
144        },
145    }
146    with open("results.json", "w") as f:
147        json.dump(result, f, indent=2)
148    print(json.dumps(result, indent=2))
149
150
151if __name__ == "__main__":
152    main()