A few holdouts for many campaigns: calibrating propensity estimates with a hierarchical model of their bias

import os

for v_ in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v_] = "1"  # nutpie runs chains as threads: a BLAS pool per chain would fight over the cores

import logging
import sys
import time
import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import nutpie
import pandas as pd
import pymc as pm
from scipy import special, stats

from pymc_challenges import data

RANDOM_SEED = 103
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning, module="arviz")
warnings.filterwarnings("ignore", message=".*Potential.*")
sys.modules.pop("jax", None)  # arviz leaves a lazy jax stub that makes matplotlib import all of JAX
pd.set_option("display.width", 160)
T_START = time.time()

BLUE, ORANGE, AQUA, PURPLE, RED = "#2a78d6", "#eb6834", "#1baf7a", "#8c5ac8", "#c8384e"
INK, MUTED, LIGHT = "#222222", "#8a8984", "#d9d8d3"
TYPES = ["retargeting", "broad", "prospecting"]
TCOL = {"retargeting": RED, "broad": AQUA, "prospecting": BLUE}
print(f"PyMC {pm.__version__}, ArviZ {az.__version__}")

data.describe("e101_criteo_uplift")
raw = np.load(data.path("e101_criteo_uplift"))
X = raw["X"].astype(float)
T = raw["treatment"].astype(int)
D_real = raw["exposure"].astype(int)
VISIT = raw["visit"].astype(float)
N = len(T)
ctl, tst = T == 0, T == 1
REACH_REAL = D_real[tst].mean()


def normal_scores(A):
    """Rank-based normal scores per column (ties share a rank)."""
    return np.column_stack([stats.norm.ppf((stats.rankdata(a) - 0.5) / len(a)) for a in A.T])


def fit_map(Z, y, prior_sd=0.5, a0=0.0, a_sd=1.5, iters=60):
    """Posterior mode of a logistic regression with N(a0, a_sd) intercept and N(0, prior_sd) slopes (Newton)."""
    Z1 = np.column_stack([np.ones(len(Z)), Z])
    prec = np.r_[1 / a_sd ** 2, np.full(Z.shape[1], 1 / prior_sd ** 2)]
    mu0 = np.r_[a0, np.zeros(Z.shape[1])]
    b = mu0.copy()
    for _ in range(iters):
        p = special.expit(Z1 @ b)
        g = Z1.T @ (y - p) - prec * (b - mu0)
        H = (Z1 * (p * (1 - p))[:, None]).T @ Z1 + np.diag(prec)
        step = np.linalg.solve(H, g)
        step *= min(1.0, 2.0 / max(np.abs(step).max(), 1e-12))    # damped: never move a log-odds by more than 2
        b += step
        if np.abs(step).max() < 1e-8:
            break
    return b


# E102's construction: the analyst sees 11 features; u is what f8 knows about visiting beyond them
R = normal_scores(X)
Zall = np.column_stack([R, R ** 2])
Zall = (Zall - Zall.mean(0)) / Zall.std(0)
HIDDEN = 8
SEEN = [k for k in range(12) if k != HIDDEN]
Zseen = Zall[:, SEEN + [12 + k for k in SEEN]]
tu = tst & (D_real == 0)
a0_v = special.logit(VISIT[tu].mean())
shift = (np.column_stack([np.ones(N), Zall]) @ fit_map(Zall[tu], VISIT[tu], a0=a0_v)
         - np.column_stack([np.ones(N), Zseen]) @ fit_map(Zseen[tu], VISIT[tu], a0=a0_v))
Zp = Zseen[ctl]
Z1p = np.column_stack([np.ones(len(Zp)), Zp])
Y0 = VISIT[ctl]                                                     # the audience's no-ad visits
n = len(Y0)
resid = shift[ctl] - Z1p @ np.linalg.lstsq(Z1p, shift[ctl], rcond=None)[0]
u = stats.norm.ppf((stats.rankdata(resid) - 0.5) / n)               # hidden engagement, N(0, 1)
eta = Z1p @ fit_map(Zseen[tst], D_real[tst], a0=special.logit(REACH_REAL))
eta = (eta - eta.mean()) / 1.7                                      # the real campaign's visible targeting
print(f"audience: {n:,} control users, visit rate {Y0.mean():.4f}; hidden trait u: corr with visits "
      f"{np.corrcoef(u, Y0)[0, 1]:.3f}")

E_CAP = 0.99
AUD = 20_000
REACH_RANGE = {"retargeting": (0.04, 0.30), "broad": (0.10, 0.60), "prospecting": (0.15, 0.70)}
GAMMA_MEAN = {"retargeting": 0.7, "broad": 0.0, "prospecting": -0.7}
SCENARIOS = {"tight": 0.1, "moderate": 0.4, "loose": 0.8, "unpredictable": None}


def aipw(d, idx, y=None):
    """The analyst's AIPW effect on the exposed (points) and its influence-function sd, on audience rows idx."""
    Z, Z1 = Zp[idx], Z1p[idx]
    y = Y0[idx] if y is None else y
    d = d.astype(float)
    t, un = d == 1, d == 0
    e = np.minimum(special.expit(Z1 @ fit_map(Z, d, a0=special.logit(d.mean()))), E_CAP)
    w = e / (1 - e)
    m0 = special.expit(Z1 @ fit_map(Z[un], y[un], a0=special.logit(max(y[un].mean(), 1e-4))))
    mu_w = np.average(y[un] - m0[un], weights=w[un])
    att = np.mean(y[t] - m0[t]) - mu_w
    phi = (t * (y - m0 - (att + mu_w)) - un * w * (y - m0 - mu_w) * t.sum() / w[un].sum()) / t.mean()
    return 100 * att, 100 * phi.std() / np.sqrt(len(y))


def deliver(kind, scen, rng_):
    """One campaign: its audience, delivery rule and exposures (who would be shown the ad)."""
    idx = np.sort(rng_.choice(n, AUD, replace=False))
    reach = rng_.uniform(*REACH_RANGE[kind])
    scale = rng_.uniform(0.5, 1.5)                                   # how hard it targets the visible features
    if SCENARIOS[scen] is None:
        gamma, supp = rng_.normal(0.0, 0.7), 0.0
    else:
        gamma = GAMMA_MEAN[kind] + rng_.normal(0.0, SCENARIOS[scen])
        supp = rng_.uniform(0.05, 0.2) if kind == "broad" else 0.0
    L = scale * eta[idx] + gamma * u[idx] + rng_.normal(size=AUD)
    if supp > 0:
        L[u[idx] > np.quantile(u[idx], 1 - supp)] = -np.inf          # suppression list: most engaged never exposed
    d = L > np.quantile(L, 1 - reach)
    return idx, d, {"type": kind, "scenario": scen, "reach": reach, "gamma": gamma, "suppressed": supp}


def make_campaign(kind, scen, rng_):
    idx, d, info = deliver(kind, scen, rng_)
    bias, ps_sd = aipw(d, idx)                                       # placebo: true effect 0, so this IS the bias
    return {**info, "bias": bias, "ps_sd": ps_sd, "base": 100 * Y0[idx][d].mean(), "exposed": int(d.sum()),
            "rows": idx[d].astype(np.int32)}


t0 = time.time()
rb = np.random.default_rng(RANDOM_SEED + 1)
bank = [make_campaign(k, sc, rb) for sc in SCENARIOS for k in TYPES for _ in range(60)]
ROWS = [c.pop("rows") for c in bank]                                 # each campaign's would-be-exposed users
bank = pd.DataFrame(bank)
print(f"{len(bank)} campaigns delivered and analysed in {time.time() - t0:.0f} s")
tab = bank.groupby(["scenario", "type"], sort=False).agg(
    bias=("bias", "mean"), bias_sd=("bias", "std"), ps_sd=("ps_sd", "mean"), base=("base", "mean"),
    reach=("reach", "mean"), exposed=("exposed", "mean"))
print(tab.round(2).to_string())
within = bank.groupby(["scenario", "type"], sort=False).bias.std().groupby("scenario", sort=False).mean()
print("bias sd within type, averaged over types: " + ", ".join(f"{s} {v:.1f}" for s, v in within.items()))

J, NS = 40, 1500                                                     # campaigns per portfolio, respondents per arm
LIFT_MEAN = {"retargeting": 2.0, "broad": 3.0, "prospecting": 3.5}


def experiment(j_bank, theta, rng_):
    """A ghost-ad holdout survey from real rows: NS exposed and NS held-out would-be-exposed respondents."""
    rows = ROWS[j_bank]
    y_c = Y0[rng_.choice(rows, NS)]                                  # held out: their real no-ad visits
    y_t = Y0[rng_.choice(rows, NS)]
    p = theta / (100 - bank.base[j_bank])                            # the ad converts this share of non-visitors
    y_t = np.where(y_t == 1, 1.0, (rng_.random(NS) < p).astype(float))
    return 100 * (y_t.mean() - y_c.mean()), 100 * np.sqrt(y_t.var() / NS + y_c.var() / NS)


def new_portfolio(scen, rng_):
    """40 campaigns from the scenario's bank, fresh true effects, costs and experiments."""
    pick = rng_.choice(np.flatnonzero(bank.scenario == scen), J, replace=False)
    port = bank.loc[pick, ["type", "reach", "gamma", "suppressed", "bias", "ps_sd", "base"]].reset_index(drop=True)
    port["bank"] = pick
    mean = port.type.map(LIFT_MEAN).to_numpy() * (1.3 - 0.9 * port.reach.to_numpy())
    port["theta"] = rng_.gamma(2.25, mean / 2.25)                    # true effect on the exposed (points)
    port["cost"] = rng_.uniform(1.0, 6.0, J)                         # break-even effect (points)
    port["est"] = port.theta + port.bias                             # the propensity estimate
    ex = [experiment(j, th, rng_) for j, th in zip(pick, port.theta)]
    port["exp"], port["se"] = [e[0] for e in ex], [e[1] for e in ex]
    return port


def test_order(port, rule, rng_, susp_sd=1.0):
    """Order in which campaigns get a holdout; the first k are tested."""
    ty, rz = port.type.to_numpy(), port.reach.to_numpy()
    if rule == "random":
        return rng_.permutation(J)
    if rule in ("designed", "suspicious within type"):
        if rule == "designed":                                      # spread over reach within each type
            vdc = lambda i: (sum(((i >> b) & 1) / 2 ** (b + 1) for b in range(8)) + 0.5) % 1
            lists = []
            for t in TYPES:
                g = np.flatnonzero(ty == t)
                g = g[np.argsort(rz[g])]
                o = [g[min(int(vdc(i) * len(g)), len(g) - 1)] for i in range(len(g))]
                o = list(dict.fromkeys(o))
                lists.append(o + [x for x in g if x not in o])
        else:                                                       # the team's hunch: bias plus noise
            sig = port.bias.to_numpy() + rng_.normal(0, susp_sd, J)
            lists = [list(np.flatnonzero(ty == t)[np.argsort(-sig[ty == t])]) for t in TYPES]
        s0 = rng_.integers(3)
        lists = lists[s0:] + lists[:s0]
        out = []
        while any(lists):
            for li in lists:
                if li:
                    out.append(li.pop(0))
        return np.array(out)
    if rule == "top estimates":
        return np.argsort(-port.est.to_numpy())
    if rule == "suspicious":
        return np.argsort(-(port.bias.to_numpy() + rng_.normal(0, susp_sd, J)))
    if rule == "retargeting only":
        g = rng_.permutation(np.flatnonzero(ty == "retargeting"))
        return np.r_[g, rng_.permutation(np.setdiff1d(np.arange(J), g))]


rp = np.random.default_rng(RANDOM_SEED + 3)
port = new_portfolio("moderate", rp)
K_MAIN = 12
tested = np.zeros(J, bool)
tested[test_order(port, "designed", rp)[:K_MAIN]] = True
port["tested"] = tested
show = port.copy()
show["raw bias"] = np.where(tested, show.est - show.exp, np.nan)
print(f"one portfolio: {J} campaigns, {tested.sum()} with a holdout; experiment sd {port.se[tested].min():.2f}-"
      f"{port.se[tested].max():.2f} points")
print(show[tested].sort_values("type")[["type", "reach", "theta", "est", "exp", "se", "raw bias", "bias"]]
      .round(2).to_string(index=False))
print(port.groupby("type")[["theta", "est", "bias"]].mean().round(2).to_string())

def model_arrays(port):
    t_idx = pd.Index(TYPES).get_indexer(port.type)
    rz = ((port.reach - 0.35) / 0.2).to_numpy()
    return t_idx, rz


def latent_model(port, mask):
    """The model as written: every campaign's true effect is a latent variable."""
    t_idx, rz = model_arrays(port)
    with pm.Model(coords={"campaign": np.arange(J), "type": TYPES}) as m:
        a_type = pm.Normal("a_type", 3.0, 3.0, dims="type")         # mean true effect by type
        a_reach = pm.Normal("a_reach", 0.0, 1.5)
        s = pm.HalfNormal("s", 3.0)                                  # spread of true effects
        z = pm.Normal("z", 0.0, 1.0, dims="campaign")
        theta = pm.Deterministic("theta", a_type[t_idx] + a_reach * rz + s * z, dims="campaign")
        b_type = pm.Normal("b_type", 0.0, 4.0, dims="type")         # mean bias by type
        b_reach = pm.Normal("b_reach", 0.0, 2.0)
        tau = pm.HalfNormal("tau", 4.0)                              # spread of the bias around its prediction
        pm.Normal("est", theta + b_type[t_idx] + b_reach * rz, tau, observed=port.est.to_numpy(), dims="campaign")
        pm.Potential("experiments", (mask * pm.logp(pm.Normal.dist(theta, port.se.to_numpy()),
                                                    np.nan_to_num(port.exp.to_numpy()))).sum())
    return m


t0 = time.time()
with latent_model(port, port.tested.to_numpy().astype(float)):
    idata_lat = pm.sample(random_seed=RANDOM_SEED, progressbar=False, target_accept=0.9)
sm_ = az.summary(idata_lat, var_names=["s", "tau"], round_to=4)
print(f"latent model: {time.time() - t0:.0f} s, divergences {int(idata_lat.sample_stats['diverging'].sum())}; "
      + "; ".join(f"{v}: r_hat {sm_.loc[v, 'r_hat']:.2f}, bulk ESS {sm_.loc[v, 'ess_bulk']:.0f}" for v in ["s", "tau"]))
del idata_lat

def calibration_model(port, mask):
    """The same model with every campaign's true effect integrated out (two Gaussian factors per campaign)."""
    t_idx, rz = model_arrays(port)
    est, ex, se = port.est.to_numpy(), np.nan_to_num(port.exp.to_numpy()), port.se.to_numpy()
    with pm.Model(coords={"campaign": np.arange(J), "type": TYPES}) as m:
        mask_ = pm.Data("mask", mask.astype(float), dims="campaign")
        a_type = pm.Normal("a_type", 3.0, 3.0, dims="type")
        a_reach = pm.Normal("a_reach", 0.0, 1.5)
        s = pm.HalfNormal("s", 3.0)
        b_type = pm.Normal("b_type", 0.0, 4.0, dims="type")
        b_reach = pm.Normal("b_reach", 0.0, 2.0)
        tau = pm.HalfNormal("tau", 4.0)
        m_j = a_type[t_idx] + a_reach * rz                           # expected true effect
        mu_b = b_type[t_idx] + b_reach * rz                          # predicted bias
        v = s ** 2 + tau ** 2
        pm.Normal("est", m_j + mu_b, pm.math.sqrt(v), observed=est, dims="campaign")
        k_ = s ** 2 / v                                              # the experiment, given the propensity estimate
        pm.Potential("experiments", (mask_ * pm.logp(
            pm.Normal.dist(m_j + k_ * (est - m_j - mu_b), pm.math.sqrt(s ** 2 * tau ** 2 / v + se ** 2)), ex)).sum())
    return m


def theta_draws(fit, port, mask, rng_):
    """True effects drawn from their exact conditional given each posterior draw of the hyperparameters."""
    po = az.extract(fit, var_names=["a_type", "a_reach", "s", "b_type", "b_reach", "tau"])
    t_idx, rz = model_arrays(port)
    m_j = po["a_type"].values[t_idx] + po["a_reach"].values * rz[:, None]
    mu_b = po["b_type"].values[t_idx] + po["b_reach"].values * rz[:, None]
    s, tau = po["s"].values, po["tau"].values
    se2 = port.se.to_numpy() ** 2
    prec = 1 / s ** 2 + 1 / tau ** 2 + (mask / se2)[:, None]
    num = (m_j / s ** 2 + (port.est.to_numpy()[:, None] - mu_b) / tau ** 2
           + (mask * np.nan_to_num(port.exp.to_numpy()) / se2)[:, None])
    return num / prec + rng_.normal(size=prec.shape) / np.sqrt(prec), po


cal = calibration_model(port, port.tested.to_numpy())
with cal:
    prior = pm.sample_prior_predictive(1000, random_seed=RANDOM_SEED)
pr = prior.prior
b_new = (pr["b_type"].values[0, :, :, None] + pr["tau"].values[0, :, None, None] * rp.normal(size=(1000, 3, 1)))
th_prior = pr["a_type"].values[0] + pr["s"].values[0, :, None] * rp.normal(size=(1000, 3))
pp_est = prior.prior_predictive["est"].values.ravel()
print("prior predictive, propensity estimates (points): 5/50/95% "
      + " / ".join(f"{v:.1f}" for v in np.percentile(pp_est, [5, 50, 95])))
print("prior predictive, bias of a new campaign of one type: 5/50/95% "
      + " / ".join(f"{v:.1f}" for v in np.percentile(b_new, [5, 50, 95])))
print("prior, true effect of a campaign: 5/50/95% "
      + " / ".join(f"{v:.1f}" for v in np.percentile(th_prior, [5, 50, 95])))
print(f"observed propensity estimates in the portfolio: {port.est.min():.1f} to {port.est.max():.1f}")

SG = np.exp(np.linspace(np.log(0.3), np.log(8.0), 16))              # grid over s ...
TG = np.exp(np.linspace(np.log(0.2), np.log(10.0), 22))             # ... and tau, uniform in the logs
GS, GT = [g.ravel() for g in np.meshgrid(SG, TG, indexing="ij")]
LOGPRIOR = stats.halfnorm.logpdf(GS, scale=3) + stats.halfnorm.logpdf(GT, scale=4) + np.log(GS * GT)


def designs(port):
    """Prior means/sds of the linear terms for the two calibration models."""
    oh = np.column_stack([(port.type == t).to_numpy(float) for t in TYPES])
    rz = ((port.reach - 0.35) / 0.2).to_numpy()
    one = np.ones((len(port), 1))
    cov = np.column_stack([oh, rz])
    return {"calibrated, covariates": (cov, np.r_[3.0, 3.0, 3.0, 0.0], np.r_[3.0, 3.0, 3.0, 1.5], cov,
                                       np.r_[4.0, 4.0, 4.0, 2.0]),
            "calibrated, pooled": (one, np.r_[3.0], np.r_[3.0], one, np.r_[4.0])}


def grid_posterior(est, ex, se, tested, Xt, tm, tsd, Xb, bsd):
    """Exact posterior of every theta on the (s, tau) grid, with the linear terms integrated out.
    Returns grid weights (G,), conditional means and sds of theta (G, J), and of tau's mean."""
    Jn, pt_, pb = len(est), Xt.shape[1], Xb.shape[1]
    ti = np.flatnonzero(tested)
    k, q = len(ti), pt_ + pb + Jn                                    # latent w = (alpha, beta, delta)
    H = np.zeros((Jn + k, q))
    H[:Jn, :pt_], H[:Jn, pt_:pt_ + pb], H[:Jn, pt_ + pb:] = Xt, Xb, np.eye(Jn)
    H[Jn:, :pt_] = Xt[ti]
    H[Jn + np.arange(k), pt_ + pb + ti] = 1
    y = np.r_[est, ex[ti]]
    mu0 = np.r_[tm, np.zeros(pb + Jn)]
    G = len(GS)
    pv = np.empty((G, q))
    pv[:, :pt_], pv[:, pt_:pt_ + pb], pv[:, pt_ + pb:] = tsd ** 2, bsd ** 2, GS[:, None] ** 2
    nz = np.empty((G, Jn + k))
    nz[:, :Jn], nz[:, Jn:] = GT[:, None] ** 2, se[ti] ** 2
    PVH = pv[:, :, None] * H.T[None]
    S = H[None] @ PVH + nz[:, :, None] * np.eye(Jn + k)              # marginal covariance of the data
    r = y - H @ mu0
    a = np.linalg.solve(S, np.broadcast_to(r, (G, Jn + k))[..., None])[..., 0]
    lw = -0.5 * a @ r - np.log(np.diagonal(np.linalg.cholesky(S), axis1=1, axis2=2)).sum(1) + LOGPRIOR
    wg = np.exp(lw - lw.max())
    wg /= wg.sum()
    A = np.zeros((Jn, q))
    A[:, :pt_], A[:, pt_ + pb:] = Xt, np.eye(Jn)                     # theta = Xt alpha + delta
    APV = A[None] * pv[:, None, :]
    C = APV @ H.T[None]
    mean = (A @ mu0)[None] + np.einsum("gjn,gn->gj", C, a)
    var = np.einsum("gjq,jq->gj", APV, A) - np.einsum("gjn,gnj->gj", C, np.linalg.solve(S, C.transpose(0, 2, 1)))
    return wg, mean, np.sqrt(np.maximum(var, 1e-12))


Xt, tm, tsd, Xb, bsd = designs(port)["calibrated, covariates"]
t0 = time.time()
wg, mn, sd = grid_posterior(port.est.to_numpy(), port.exp.fillna(0).to_numpy(), port.se.to_numpy(),
                            port.tested.to_numpy(), Xt, tm, tsd, Xb, bsd)
print(f"grid posterior in {1000 * (time.time() - t0):.0f} ms")
g_mean = wg @ mn
g_sd = np.sqrt(wg @ (sd ** 2 + mn ** 2) - g_mean ** 2)
print(f"theta, grid against PyMC: largest gap in posterior means {np.max(np.abs(g_mean - TH_post.mean(1))):.2f} "
      f"points, in sds {np.max(np.abs(g_sd - TH_post.std(1))):.2f} (posterior sds {TH_post.std(1).min():.2f}-"
      f"{TH_post.std(1).max():.2f})")
print(f"tau: grid mean {wg @ GT:.2f}, PyMC mean {post['tau'].values.mean():.2f}")

KS = [0, 2, 4, 8, 16]


def evaluate(port, tested):
    """Decisions of every policy on one portfolio with the given holdouts, scored against the truth."""
    th, c, est = port.theta.to_numpy(), port.cost.to_numpy(), port.est.to_numpy()
    ex, se = port.exp.to_numpy(), port.se.to_numpy()
    k = tested.sum()
    dec = {"oracle": th > c, "raw propensity": est > c,
           "experiment only": np.where(tested, ex > c, (ex[tested].mean() if k else 3.0) > c),
           "average correction": est - (np.mean((est - ex)[tested]) if k else 0.0) > c}
    extra = {"raw propensity": {"cover90": np.mean((np.abs(est - th) < 1.645 * port.ps_sd)[~tested]),
                                "error": np.mean((est - th)[~tested])},
             "average correction": {"correction": np.mean((est - ex)[tested]) if k else 0.0}}
    for name, (Xt_, tm_, tsd_, Xb_, bsd_) in designs(port).items():
        wg_, mn_, sd_ = grid_posterior(est, ex, se, tested, Xt_, tm_, tsd_, Xb_, bsd_)
        m_ = wg_ @ mn_
        pit = wg_ @ stats.norm.cdf((th[None] - mn_) / sd_)
        dec[name] = m_ > c
        extra[name] = {"cover90": np.mean(((pit > 0.05) & (pit < 0.95))[~tested]),
                       "error": np.mean((m_ - th)[~tested]), "tau": wg_ @ GT}
    best = np.maximum(th - c, 0)
    rows = []
    for pol, d in dec.items():
        loss = best - np.where(d, th - c, 0)
        rows.append({"policy": pol, "regret": loss.mean(), "wrong": np.mean(d != dec["oracle"]),
                     "wrong untested": np.mean((d != dec["oracle"])[~tested]), "best": best.mean(),
                     "worth": np.mean(th > c),
                     **{f"regret {t}": loss[(port.type == t).to_numpy()].mean() for t in TYPES}, **extra.get(pol, {})})
    return rows


REPS = 30
t0 = time.time()
rs = np.random.default_rng(RANDOM_SEED + 4)
res = []
for scen in SCENARIOS:
    for rep in range(REPS):
        pf = new_portfolio(scen, rs)
        for rule in ["random", "designed"]:
            o = test_order(pf, rule, rs)
            for k in KS:
                if k == 0 and rule == "designed":
                    continue
                t_ = np.zeros(J, bool)
                t_[o[:k]] = True
                res += [{**x, "scenario": scen, "rule": rule, "k": k, "rep": rep} for x in evaluate(pf, t_)]
res = pd.DataFrame(res)
print(f"{REPS * len(SCENARIOS)} portfolios x 9 holdout plans in {time.time() - t0:.0f} s; "
      f"the oracle earns {res[res.policy == 'oracle'].best.mean():.2f} points per campaign; "
      f"{res[res.policy == 'oracle'].worth.mean():.0%} of campaigns are worth continuing")
POL = ["raw propensity", "average correction", "experiment only", "calibrated, pooled", "calibrated, covariates"]
rr = res[(res.rule == "random") & res.policy.isin(POL)]
print("regret per campaign (points of visits), random holdouts:")
print(rr.pivot_table(index=["scenario", "k"], columns="policy", values="regret", sort=False)[POL].round(2).to_string())