feat(models): SHAP top-k explainer for BBB predictions

src/models/bbb_model.py

@@ -125,3 +125,83 @@ def predict_with_proba(
         "label": label,
         "confidence": float(proba[label_idx]),
     }
+
+
+def explain_prediction(
+    model: RandomForestClassifier,
+    smiles: str,
+    top_k: int = 5,
+    n_bits: int = 2048,
+    radius: int = 2,
+) -> list[dict[str, object]]:
+    """Return the top-`top_k` feature attributions (SHAP values) for `smiles`.
+
+    Uses `shap.TreeExplainer` (exact for tree ensembles, no sampling). The
+    explanation is for the *predicted* class — i.e. SHAP values that pushed
+    the model toward whichever label was returned by `predict_with_proba`.
+
+    Reads fingerprint column names from `model._neurobridge_fp_cols` (set by
+    `train()`). Falls back to `fp_<index>` if the attribute is missing — useful
+    for models loaded from a joblib without the project-owned attribute.
+
+    Args:
+        model: Fitted classifier from `train()` or `load()`.
+        smiles: A SMILES string (validated via `is_valid_smiles`).
+        top_k: How many top features to return. Default 5 — matches the
+            jury-demo budget (more bars = noisier waterfall chart).
+        n_bits / radius: Must match training-time fingerprint settings.
+
+    Returns:
+        A list of `{"feature": "fp_<bit_idx>", "shap_value": float}` dicts,
+        sorted by `abs(shap_value)` descending.
+
+    Raises:
+        ValueError: if `smiles` cannot be parsed by RDKit.
+    """
+    import shap  # local import — heavy module, only loaded when needed
+
+    if not is_valid_smiles(smiles):
+        raise ValueError(f"invalid SMILES: {smiles!r}")
+    fp = compute_morgan_fingerprint(smiles, n_bits=n_bits, radius=radius)
+    X = fp.reshape(1, -1)
+
+    explainer = shap.TreeExplainer(model)
+    # uint8 fingerprints cause benign additivity violations in SHAP's
+    # reconstruction (base + sum != model output within tolerance); the
+    # default check produces false-positive errors on tree ensembles
+    # over quantized inputs, so we skip it.
+    shap_values = explainer.shap_values(X, check_additivity=False)
+    # `shap_values` shape varies by sklearn / shap versions:
+    #   - older: list of (1, n_features) arrays, one per class
+    #   - newer: ndarray of shape (1, n_features, n_classes) for binary RF
+    #   - or (1, n_features) when output already condensed
+    if isinstance(shap_values, list):
+        proba = model.predict_proba(X)[0]
+        label_idx = int(np.argmax(proba))
+        per_feature = shap_values[label_idx][0]
+    else:
+        arr = np.asarray(shap_values)
+        if arr.ndim == 3:
+            # (1, n_features, n_classes)
+            proba = model.predict_proba(X)[0]
+            label_idx = int(np.argmax(proba))
+            per_feature = arr[0, :, label_idx]
+        else:
+            # (1, n_features)
+            per_feature = arr[0]
+
+    fp_cols = (
+        list(model._neurobridge_fp_cols)
+        if hasattr(model, "_neurobridge_fp_cols")
+        else [f"fp_{i}" for i in range(len(per_feature))]
+    )
+
+    pairs = sorted(
+        zip(fp_cols, per_feature, strict=True),
+        key=lambda p: abs(p[1]),
+        reverse=True,
+    )
+    return [
+        {"feature": str(name), "shap_value": float(value)}
+        for name, value in pairs[:top_k]
+    ]

tests/models/test_bbb_model.py

@@ -89,3 +89,41 @@ class TestPredictWithProba:
         raw_proba = model.predict_proba(fp)[0]
         result = bbb_model.predict_with_proba(model, "CCO")
         assert abs(result["confidence"] - float(max(raw_proba))) < 1e-9
+
+
+class TestExplainPrediction:
+    def test_returns_top_k_features(self, trained_model_and_features):
+        model, _ = trained_model_and_features
+        attributions = bbb_model.explain_prediction(model, "CCO", top_k=5)
+        assert len(attributions) == 5
+        for a in attributions:
+            assert "feature" in a
+            assert "shap_value" in a
+            assert isinstance(a["shap_value"], float)
+
+    def test_features_sorted_by_absolute_shap_value_descending(
+        self, trained_model_and_features,
+    ):
+        model, _ = trained_model_and_features
+        attributions = bbb_model.explain_prediction(model, "CCO", top_k=10)
+        abs_vals = [abs(a["shap_value"]) for a in attributions]
+        assert abs_vals == sorted(abs_vals, reverse=True)
+
+    def test_features_named_fp_INDEX(self, trained_model_and_features):
+        model, _ = trained_model_and_features
+        attributions = bbb_model.explain_prediction(model, "CCO", top_k=3)
+        for a in attributions:
+            assert a["feature"].startswith("fp_")
+            int(a["feature"].split("_")[1])  # parses cleanly
+
+    def test_raises_on_invalid_smiles(self, trained_model_and_features):
+        model, _ = trained_model_and_features
+        with pytest.raises(ValueError):
+            bbb_model.explain_prediction(model, "still_not_a_smiles", top_k=5)
+
+    def test_deterministic_output(self, trained_model_and_features):
+        """AGENTS.md §4 rule 3: identical input → identical SHAP attributions."""
+        model, _ = trained_model_and_features
+        r1 = bbb_model.explain_prediction(model, "CCO", top_k=5)
+        r2 = bbb_model.explain_prediction(model, "CCO", top_k=5)
+        assert r1 == r2