1import json, time
 2import numpy as np
 3from scipy.optimize import minimize
 4from pdas_filter import solve_active_set, kkt_residual
 5
 6np.set_printoptions(precision=5, suppress=True)
 7
 8# Strictly convex QP: min 1/2 ||v-v_nom||^2, A v >= b.
 9def make_problem(v_nom, alpha=1.0, h=None, A=None):
10    n = len(v_nom)
11    if A is None:
12        A = np.vstack([np.eye(n), -np.eye(n)])
13    if h is None:
14        h = -np.ones(A.shape[0])
15    return np.eye(n), -np.asarray(v_nom), np.asarray(A, float), -alpha*np.asarray(h)
16
17
18def projected_gradient(v_nom, A, b, steps=250, lr=.12):
19    # Generic approximate baseline: gradient step followed by sequential projections.
20    v = np.zeros_like(v_nom)
21    for _ in range(steps):
22        v -= lr * (v - v_nom)
23        for a, bj in zip(A, b):
24            gap = bj - a @ v
25            if gap > 0:
26                v += gap * a / (a @ a)
27    return v
28
29# Prediction 1: for one constraint v >= alpha*q, activity changes at alpha=v_nom/q.
30v_nom = np.array([0.37]); q = 0.8
31alphas = np.linspace(0., 1.2, 25)
32active = []
33for alpha in alphas:
34    H,c,A,b = make_problem(v_nom, alpha, np.array([-q]), np.array([[1.]]))
35    _,_,act,_,_ = solve_active_set(H,c,A,b)
36    active.append(bool(act))
37pred_boundary = v_nom[0]/q
38changes = np.where(np.diff(np.array(active, dtype=int)) != 0)[0]
39obs_boundary = alphas[changes[0]+1] if len(changes) else np.nan
40
41# Prediction 2: after activation, mu = alpha*q-v_nom, slope dmu/dalpha=q.
42alpha_fit = np.linspace(pred_boundary + .08, 1.6, 12)
43mus = []
44for alpha in alpha_fit:
45    H,c,A,b = make_problem(v_nom, alpha, np.array([-q]), np.array([[1.]]))
46    _,mu,_,_,_ = solve_active_set(H,c,A,b)
47    mus.append(mu[0])
48slope = np.polyfit(alpha_fit, mus, 1)[0]
49
50# Prediction 3: smoothly moving nominal velocities preserve active sets and warm starts
51# converge in one solve iteration after the first transition.
52A = np.array([[1., 0.], [0., 1.], [-1., 0.]])
53h = np.array([-0.42, -0.25, 1.0])
54trajectory = np.column_stack([np.linspace(.1, .8, 40), np.linspace(.1, .35, 40)])
55prev = None; iterations = []; changes_count = 0; prior_act = None; max_kkt = 0.
56for vn in trajectory:
57    H,c,A2,b = make_problem(vn, 1., h, A)
58    v,mu,act,it,_ = solve_active_set(H,c,A2,b,active_init=prev)
59    iterations.append(it)
60    changes_count += int(prior_act is not None and act != prior_act)
61    prior_act = act; prev = act
62    max_kkt = max(max_kkt, kkt_residual(H,c,A2,b,v,mu))
63
64# Same small collection of changing QPs: exact PDAS versus projected-gradient approximation.
65rng = np.random.default_rng(2265)
66A4 = np.eye(4)
67h4 = -np.ones(4)*.55
68pdas_t = []; pg_t = []; pdas_err = []; pg_err = []; violations = []
69for vn in rng.normal(0, .5, size=(60,4)):
70    H,c,A5,b = make_problem(vn, 1., h4, A4)
71    t=time.perf_counter(); v,mu,act,it,_=solve_active_set(H,c,A5,b); pdas_t.append(time.perf_counter()-t)
72    t=time.perf_counter(); vp=projected_gradient(vn,A5,b); pg_t.append(time.perf_counter()-t)
73    pdas_err.append(np.linalg.norm(v-vn)**2/2); pg_err.append(np.linalg.norm(vp-vn)**2/2)
74    violations.append(max(0., np.max(b-A5@vp)))
75
76result = {
77 "prediction_boundary": {"predicted_alpha": float(pred_boundary), "observed_grid_alpha": float(obs_boundary), "grid_error": float(abs(obs_boundary-pred_boundary))},
78 "prediction_dual_slope": {"predicted_dmu_dalpha": float(q), "observed_slope": float(slope), "relative_error": float(abs(slope-q)/q)},
79 "prediction_warm_start": {"mean_iterations": float(np.mean(iterations)), "median_iterations": float(np.median(iterations)), "fraction_one_iteration": float(np.mean(np.array(iterations)==1)), "active_set_changes": int(changes_count), "max_kkt_residual": float(max_kkt)},
80 "baseline_comparison": {"pdas_mean_objective": float(np.mean(pdas_err)), "projected_gradient_mean_objective": float(np.mean(pg_err)), "pdas_mean_ms": float(1000*np.mean(pdas_t)), "projected_gradient_mean_ms": float(1000*np.mean(pg_t)), "projected_gradient_max_violation": float(max(violations))},
81 "all_iterations": iterations
82}
83print(json.dumps(result, indent=2))
84with open("results.json", "w") as f: json.dump(result, f, indent=2)