1import itertools
  2import json
  3import numpy as np
  4
  5
  6def solve_allocation(R, cost, weight, eta, budget):
  7    """Exact multiple-choice knapsack; one expert per workload."""
  8    R = np.asarray(R, float)
  9    cost = np.asarray(cost, float)
 10    weight = np.asarray(weight, float)
 11    eta = np.asarray(eta, float)
 12    K, J = R.shape
 13    best = None
 14    best_obj = -np.inf
 15    for choice in itertools.product(range(K), repeat=J):
 16        choice = np.asarray(choice, int)
 17        total_cost = np.sum(weight * cost[choice, np.arange(J)])
 18        if total_cost <= budget + 1e-12:
 19            objective = np.sum(weight * (R[choice, np.arange(J)] - eta))
 20            # deterministic tie break by lexicographic choice
 21            key = tuple(choice.tolist())
 22            if objective > best_obj + 1e-12 or (abs(objective-best_obj) <= 1e-12 and (best is None or key < tuple(best.tolist()))):
 23                best_obj = objective
 24                best = choice.copy()
 25    if best is None:
 26        raise ValueError("No feasible allocation")
 27    return best, best_obj, np.sum(weight * cost[best, np.arange(J)])
 28
 29
 30def allocation_score(q, choice, weight):
 31    q = np.asarray(q)
 32    return float(np.sum(np.asarray(weight) * q[choice, np.arange(q.shape[1])]))
 33
 34
 35def run(seed=7):
 36    rng = np.random.default_rng(seed)
 37    # Each row is an expert and each column a workload cluster.
 38    J, K = 6, 4
 39    weight = rng.uniform(.5, 1.5, J)
 40    cost = rng.uniform(.7, 2.0, (K, J))
 41    R = rng.uniform(.45, .95, (K, J))
 42    true = np.clip(R + rng.normal(0, .08, (K, J)), 0, 1)
 43    budget = .60 * np.sum(weight * np.min(cost, axis=0)) + .40 * np.sum(weight * np.max(cost, axis=0))
 44
 45    nominal, nom_obj, nom_cost = solve_allocation(R, cost, weight, np.zeros(J), budget)
 46    rows = []
 47    for eta_level in np.linspace(0, .8, 9):
 48        eta = np.full(J, eta_level)
 49        pess, pess_obj, pess_cost = solve_allocation(R, cost, weight, eta, budget)
 50        rows.append({
 51            "eta": float(eta_level),
 52            "disputed_clusters": int(np.sum(nominal != pess)),
 53            "same_route": bool(np.array_equal(nominal, pess)),
 54            "nominal_true_score": allocation_score(true, nominal, weight),
 55            "pessimistic_true_score": allocation_score(true, pess, weight),
 56            "objective_shift": float(nom_obj - pess_obj),
 57        })
 58
 59    # Broad exact sweep: prediction is zero disagreements for every common eta.
 60    total = 0
 61    violations = 0
 62    for _ in range(40):
 63        jj = int(rng.integers(2, 5)); kk = int(rng.integers(2, 5))
 64        w = rng.uniform(.1, 2, jj)
 65        c = rng.uniform(.1, 3, (kk, jj))
 66        r = rng.uniform(-2, 2, (kk, jj))
 67        min_cost = np.sum(w * np.min(c, axis=0))
 68        b = min_cost + rng.uniform(0, 1) * np.sum(w * np.ptp(c, axis=0))
 69        a, _, _ = solve_allocation(r, c, w, np.zeros(jj), b)
 70        for e in [0, .1, 1, 100]:
 71            z, _, _ = solve_allocation(r, c, w, np.full(jj, e), b)
 72            total += 1
 73            violations += int(not np.array_equal(a, z))
 74
 75    # A small comparison with a noisy nominal router and an always-cheapest baseline.
 76    # Robust policy is mathematically identical to nominal here.
 77    true_nom = allocation_score(true, nominal, weight)
 78    cheapest = np.argmin(cost, axis=0)
 79    cheap_score = allocation_score(true, cheapest, weight)
 80    predicted_shift = float(np.sum(weight * np.full(J, .8)))
 81    shift_error = max(abs(x["objective_shift"] - float(np.sum(weight * np.full(J, x["eta"])))) for x in rows)
 82    checks = {
 83        "prediction_1_zero_disagreements": violations == 0,
 84        "prediction_2_shift_error": float(shift_error),
 85        "prediction_2_pass": shift_error < 1e-9,
 86        "prediction_3_identical_quality": abs(true_nom - allocation_score(true, pess, weight)) < 1e-12,
 87    }
 88    return {
 89        "seed": seed,
 90 "shape": [J, K], "budget": float(budget),
 91        "nominal_choice": nominal.tolist(), "nominal_cost": float(nom_cost),
 92        "cheapest_choice": cheapest.tolist(),
 93        "cheapest_cost": float(np.sum(weight * cost[cheapest, np.arange(J)])),
 94        "sweep": rows,
 95        "random_sweep_instances": total, "random_sweep_disagreements": violations,
 96        "nominal_true_weighted_quality": true_nom,
 97        "robust_true_weighted_quality": true_nom,
 98        "cheapest_true_weighted_quality": cheap_score,
 99        "checks": checks,
100    }
101
102if __name__ == "__main__":
103    print(json.dumps(run(), indent=2))