1import json
  2import numpy as np
  3
  4SEED = 1099
  5rng = np.random.default_rng(SEED)
  6
  7
  8def kernel(x, y, ell=0.18, amp=1.0):
  9    x = np.asarray(x)[:, None]
 10    y = np.asarray(y)[None, :]
 11    return amp * np.exp(-0.5 * ((x-y)/ell)**2)
 12
 13
 14def posterior_var(train_x, eval_x, noise=0.04, ell=0.18):
 15    train_x = np.asarray(train_x, dtype=float)
 16    eval_x = np.asarray(eval_x, dtype=float)
 17    K = kernel(train_x, train_x, ell) + noise * np.eye(len(train_x))
 18    # solve rather than explicitly invert; this is the formula in the proposal
 19    cross = kernel(eval_x, train_x, ell)
 20    sol = np.linalg.solve(K, cross.T)
 21    v = np.diag(kernel(eval_x, eval_x, ell)) - np.sum(cross * sol.T, axis=1)
 22    return np.maximum(v, 0.0)
 23
 24
 25def reduction(train_x, action, eval_x, weights):
 26    before = posterior_var(train_x, eval_x)
 27    after = posterior_var(np.r_[train_x, action], eval_x)
 28    return float(np.dot(weights, before-after)), before, after
 29
 30
 31def make_candidates(n=9, batch=5):
 32    # Three overlapping clusters plus two isolated-ish alternatives. Every action
 33    # is a possible simulator trajectory/batch and has a spatial transition cost.
 34    centers = np.array([-0.72, -0.42, -0.10, 0.22, 0.55, 0.82])
 35    actions = []
 36    for c in centers:
 37        pts = np.linspace(c-0.10, c+0.10, batch)
 38        actions.append(np.clip(pts, -1, 1))
 39    return centers, actions
 40
 41
 42def tree_objective(train_x, actions, eval_x, weights, gamma=0.9, lam=0.0):
 43    # Exact shallow MCTS value: exhaustive tree is the small-MVP equivalent of
 44    # rollouts/UCB, and is deterministic for reproducibility.
 45    immediate = np.array([reduction(train_x, a, eval_x, weights)[0] for a in actions])
 46    best_by_root = []
 47    paths = []
 48    for i, a in enumerate(actions):
 49        best = -np.inf; bestj = None
 50        for j, b in enumerate(actions):
 51            r2 = reduction(np.r_[train_x, a], b, eval_x, weights)[0]
 52            cost = abs(i-j) / max(1, len(actions)-1)
 53            val = immediate[i] + gamma * (r2 - lam*cost)
 54            if val > best:
 55                best, bestj = val, j
 56        best_by_root.append(best)
 57        paths.append((i, bestj))
 58    root = int(np.argmax(best_by_root))
 59    return root, paths[root], float(best_by_root[root]), immediate
 60
 61
 62def overlap_sweep():
 63    # Fixed first batch; construct second batches with controlled point overlap.
 64    train = np.array([-1.0, 1.0])
 65    z = np.linspace(-1, 1, 401); w = np.ones(len(z))/len(z)
 66    a = np.linspace(-0.10, 0.10, 8)
 67    # Same candidate location but replace points with distant points to control overlap.
 68    far = np.linspace(0.62, 0.82, 8)
 69    rows=[]
 70    for frac in np.linspace(0, 1, 5):
 71        m = int(round(frac*len(a)))
 72        b = np.r_[a[:m], far[:len(a)-m]]
 73        # marginal reward after observing a
 74        r, _, _ = reduction(np.r_[train, a], b, z, w)
 75        rows.append((float(frac), r))
 76    return rows
 77
 78
 79def lambda_sweep():
 80    train = np.array([-1.0, 1.0])
 81    z = np.linspace(-1, 1, 401); w = np.ones(len(z))/len(z)
 82    centers, actions = make_candidates(); gamma = .9
 83    # Every fixed depth-2 path has affine value A - B*lambda. Precompute these
 84    # lines once, then find the exact first root-policy transition.
 85    lines = {}
 86    for i in range(len(actions)):
 87        r1 = reduction(train, actions[i], z, w)[0]
 88        for j in range(len(actions)):
 89            r2 = reduction(np.r_[train, actions[i]], actions[j], z, w)[0]
 90            lines[(i,j)] = (r1 + gamma*r2, gamma*abs(i-j)/5)
 91    def policy(lam):
 92        vals=[]
 93        paths=[]
 94        for i in range(len(actions)):
 95            candidates=[(lines[(i,j)][0]-lam*lines[(i,j)][1],j) for j in range(len(actions))]
 96            val,j=max(candidates)
 97            vals.append(val); paths.append((i,j))
 98        root=int(np.argmax(vals))
 99        return root,paths[root]
100    grid=np.linspace(0,2,401)
101    chosen=[policy(x)[0] for x in grid]
102    changes=np.where(np.diff(chosen)!=0)[0]
103    if len(changes):
104        idx=int(changes[0]); observed=(idx+0.5)*2/400
105        old_path=policy(max(0,observed-1e-7))[1]
106        new_path=policy(min(2,observed+1e-7))[1]
107        a0,a1=lines[old_path]; b0,b1=lines[new_path]
108        predicted=(a0-b0)/(a1-b1) if abs(a1-b1)>1e-12 else None
109        return {'observed_threshold':float(observed),'predicted_threshold':None if predicted is None else float(predicted),
110                'old_root':int(chosen[idx]),'new_root':int(chosen[idx+1]),'num_changes':int(len(changes))}
111    return {'observed_threshold':None,'predicted_threshold':None,'num_changes':0}
112
113def main():
114    # Core numerical formula sanity check: conditioning cannot increase variance.
115    train=np.array([-0.85,-0.25,0.35,0.9]); z=np.linspace(-1,1,301)
116    v0=posterior_var(train,z); v1=posterior_var(np.r_[train,0.0],z)
117    max_increase=float(np.max(v1-v0))
118    overlap=overlap_sweep()
119    # Check monotonic diminishing returns and report endpoint ratio.
120    overlap_monotone=all(overlap[i+1][1] <= overlap[i][1]+1e-10 for i in range(len(overlap)-1))
121    lam=lambda_sweep()
122    # Mini comparison: random, greedy, and planned tree, two executed batches.
123    train0=np.array([-1.0,1.0]); z=np.linspace(-1,1,401); w=np.ones(len(z))/len(z)
124    centers, actions=make_candidates(); gamma=.9
125    def run(method):
126        tr=train0.copy(); total=0.; roots=[]
127        for t in range(2):
128            if method=='random': i=int(rng.integers(len(actions)))
129            elif method=='greedy': i=int(np.argmax([reduction(tr,a,z,w)[0] for a in actions]))
130            else: i=tree_objective(tr,actions,z,w,gamma,lam=0.8)[0]
131            total += reduction(tr,actions[i],z,w)[0]; tr=np.r_[tr,actions[i]]; roots.append(i)
132        return total,roots, float(np.mean(posterior_var(tr,z)))
133    cmp={m:run(m) for m in ['random','greedy','tree']}
134    out={'seed':SEED,'math_check':{'max_variance_increase':max_increase,
135          'nonincrease_verified':max_increase < 1e-9},
136         'prediction_1':{'claim':'posterior variance never increases after conditioning',
137                         'predicted':'max increase <= numerical tolerance 1e-9', 'observed':max_increase},
138         'prediction_2':{'claim':'overlap makes later marginal reduction smaller',
139                         'overlap_fraction_and_reduction':overlap,'monotone_verified':overlap_monotone,
140                         'endpoint_reduction_ratio':overlap[-1][1]/max(overlap[0][1],1e-12)},
141         'prediction_3':{'claim':'tree root changes at the analytic reward/cost crossover',**lam},
142         'comparison':{k:{'cumulative_reduction':v[0],'actions':v[1],'final_mean_variance':v[2]} for k,v in cmp.items()}}
143    with open('results.json','w') as f: json.dump(out,f,indent=2)
144    print(json.dumps(out,indent=2))
145
146if __name__=='__main__': main()