|
| 1 | +"""Auditable Bayesian multi-objective acquisition and portfolio sensitivity.""" |
| 2 | +from __future__ import annotations |
| 3 | + |
| 4 | +import hashlib |
| 5 | +import json |
| 6 | +import pathlib |
| 7 | +from typing import Any, Iterable |
| 8 | + |
| 9 | +from motif.acquisition import rank_portfolio |
| 10 | + |
| 11 | + |
| 12 | +_SOURCE_ROOT = pathlib.Path(__file__).resolve().parent |
| 13 | +ACQUISITION_SOURCE_DIGESTS = { |
| 14 | + name: hashlib.sha256((_SOURCE_ROOT / name).read_bytes()).hexdigest() |
| 15 | + for name in ("acquisition.py", "advanced_acquisition.py") |
| 16 | +} |
| 17 | + |
| 18 | + |
| 19 | +def _digest(value: Any) -> str: |
| 20 | + return "sha256:" + hashlib.sha256(json.dumps( |
| 21 | + value, sort_keys=True, separators=(",", ":"), allow_nan=False).encode()).hexdigest() |
| 22 | + |
| 23 | + |
| 24 | +def botorch_qehvi(train_features, train_objectives, candidate_features, *, |
| 25 | + reference_point: Iterable[float], mc_samples: int = 128, |
| 26 | + seed: int = 0) -> dict[str, Any]: |
| 27 | + """Fit an exact multitask GP and score each candidate by q=1 noisy EHVI. |
| 28 | +
|
| 29 | + Objective columns must already be oriented so larger is better. Hard constraints |
| 30 | + are intentionally outside this function and remain exact Dirac gates. |
| 31 | + """ |
| 32 | + import torch |
| 33 | + from botorch.acquisition.multi_objective import qLogExpectedHypervolumeImprovement |
| 34 | + from botorch.fit import fit_gpytorch_mll |
| 35 | + from botorch.models import SingleTaskGP |
| 36 | + from botorch.sampling.normal import SobolQMCNormalSampler |
| 37 | + from botorch.utils.multi_objective.box_decompositions.non_dominated import ( |
| 38 | + NondominatedPartitioning, |
| 39 | + ) |
| 40 | + from gpytorch.mlls import ExactMarginalLogLikelihood |
| 41 | + |
| 42 | + torch.manual_seed(seed) |
| 43 | + x = torch.as_tensor(train_features, dtype=torch.double) |
| 44 | + y = torch.as_tensor(train_objectives, dtype=torch.double) |
| 45 | + candidates = torch.as_tensor(candidate_features, dtype=torch.double) |
| 46 | + ref = torch.as_tensor(list(reference_point), dtype=torch.double) |
| 47 | + if x.ndim != 2 or y.ndim != 2 or candidates.ndim != 2: |
| 48 | + raise ValueError("qEHVI requires 2D train/candidate matrices") |
| 49 | + if len(x) != len(y) or x.shape[1] != candidates.shape[1] or y.shape[1] != len(ref): |
| 50 | + raise ValueError("qEHVI matrix dimensions are inconsistent") |
| 51 | + if len(x) < 3 or y.shape[1] < 2: |
| 52 | + raise ValueError("qEHVI requires >=3 observations and >=2 objectives") |
| 53 | + model = SingleTaskGP(x, y) |
| 54 | + fit_gpytorch_mll(ExactMarginalLogLikelihood(model.likelihood, model)) |
| 55 | + partitioning = NondominatedPartitioning(ref_point=ref, Y=y) |
| 56 | + sampler = SobolQMCNormalSampler(sample_shape=torch.Size([mc_samples]), seed=seed) |
| 57 | + acquisition = qLogExpectedHypervolumeImprovement( |
| 58 | + model=model, ref_point=ref.tolist(), partitioning=partitioning, sampler=sampler) |
| 59 | + with torch.no_grad(): |
| 60 | + log_scores = acquisition(candidates.unsqueeze(1)).cpu() |
| 61 | + scores = torch.exp(log_scores).tolist() |
| 62 | + posterior = model.posterior(candidates) |
| 63 | + mean = posterior.mean.cpu().tolist() |
| 64 | + variance = posterior.variance.cpu().tolist() |
| 65 | + release = { |
| 66 | + "schema_version": "1.0", "kind": "botorch_qlogehvi", |
| 67 | + "botorch_version": __import__("botorch").__version__, |
| 68 | + "reference_point": ref.tolist(), "mc_samples": mc_samples, "seed": seed, |
| 69 | + "scores": [float(value) for value in scores], |
| 70 | + "posterior_mean": mean, "posterior_variance": variance, |
| 71 | + "policy": "q=1 independent candidate scoring; batch diversity applied separately", |
| 72 | + } |
| 73 | + release["digest"] = _digest(release) |
| 74 | + return release |
| 75 | + |
| 76 | + |
| 77 | +def information_value(posterior_variance: Iterable[Iterable[float]], |
| 78 | + costs: Iterable[float]) -> list[float]: |
| 79 | + """Transparent variance-reduction-per-cost proxy used beside, never inside, EHVI.""" |
| 80 | + import numpy as np |
| 81 | + |
| 82 | + variance = np.asarray(list(posterior_variance), dtype=float) |
| 83 | + cost = np.asarray(list(costs), dtype=float) |
| 84 | + if variance.ndim != 2 or len(variance) != len(cost) or (cost <= 0).any(): |
| 85 | + raise ValueError("VOI requires objective variances and strictly positive costs") |
| 86 | + return (variance.sum(axis=1) / cost).astype(float).tolist() |
| 87 | + |
| 88 | + |
| 89 | +def greedy_diversity(similarity, *, order: Iterable[int]) -> list[float]: |
| 90 | + """Score novelty against the already selected batch using a disclosed greedy order.""" |
| 91 | + import numpy as np |
| 92 | + |
| 93 | + matrix = np.asarray(similarity, dtype=float) |
| 94 | + indices = list(order) |
| 95 | + if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]: |
| 96 | + raise ValueError("diversity requires a square similarity matrix") |
| 97 | + scores = [0.0] * len(matrix) |
| 98 | + selected: list[int] = [] |
| 99 | + for index in indices: |
| 100 | + scores[index] = 1.0 if not selected else 1.0 - float(matrix[index, selected].max()) |
| 101 | + selected.append(index) |
| 102 | + return scores |
| 103 | + |
| 104 | + |
| 105 | +def selection_sensitivity(candidates: Iterable[dict[str, Any]], *, |
| 106 | + objectives: Iterable[dict[str, str]], |
| 107 | + hard_constraints: Iterable[dict[str, Any]], capacity: int, |
| 108 | + relative_perturbations: Iterable[float] = (-0.1, 0.1)) -> dict[str, Any]: |
| 109 | + """Re-run the disclosed policy over objective perturbations and capacity ±1.""" |
| 110 | + rows = [dict(row) for row in candidates] |
| 111 | + axes = list(objectives) |
| 112 | + baseline = rank_portfolio(rows, objectives=axes, hard_constraints=hard_constraints, |
| 113 | + capacity=capacity) |
| 114 | + baseline_ids = {row["proposal_id"] for row in baseline["selected"]} |
| 115 | + scenarios = [] |
| 116 | + capacities = sorted({max(0, capacity - 1), capacity, capacity + 1}) |
| 117 | + for scenario_capacity in capacities: |
| 118 | + result = rank_portfolio(rows, objectives=axes, hard_constraints=hard_constraints, |
| 119 | + capacity=scenario_capacity) |
| 120 | + selected = {row["proposal_id"] for row in result["selected"]} |
| 121 | + union = baseline_ids | selected |
| 122 | + scenarios.append({"kind": "capacity", "value": scenario_capacity, |
| 123 | + "selected_ids": sorted(selected), |
| 124 | + "jaccard": len(baseline_ids & selected) / len(union) if union else 1.0}) |
| 125 | + for axis in axes: |
| 126 | + for perturbation in relative_perturbations: |
| 127 | + changed = [] |
| 128 | + for row in rows: |
| 129 | + copy = {**row, "objectives": dict(row["objectives"])} |
| 130 | + value = float(copy["objectives"][axis["key"]]) |
| 131 | + copy["objectives"][axis["key"]] = value * (1.0 + perturbation) |
| 132 | + changed.append(copy) |
| 133 | + result = rank_portfolio(changed, objectives=axes, |
| 134 | + hard_constraints=hard_constraints, capacity=capacity) |
| 135 | + selected = {row["proposal_id"] for row in result["selected"]} |
| 136 | + union = baseline_ids | selected |
| 137 | + scenarios.append({"kind": "objective", "axis": axis["key"], |
| 138 | + "relative_perturbation": perturbation, |
| 139 | + "selected_ids": sorted(selected), |
| 140 | + "jaccard": len(baseline_ids & selected) / len(union) if union else 1.0}) |
| 141 | + minimum = min((row["jaccard"] for row in scenarios), default=1.0) |
| 142 | + report = {"schema_version": "1.0", "kind": "portfolio_sensitivity", |
| 143 | + "baseline_selected_ids": sorted(baseline_ids), "scenarios": scenarios, |
| 144 | + "minimum_jaccard": minimum, |
| 145 | + "unstable": minimum < 0.5} |
| 146 | + report["digest"] = _digest(report) |
| 147 | + return report |
| 148 | + |
| 149 | + |
| 150 | +__all__ = [ |
| 151 | + "botorch_qehvi", "greedy_diversity", "information_value", "selection_sensitivity", |
| 152 | +] |
0 commit comments