perf+robustness: O(log n) conformal p-value, sklearn NotFittedError, unattainable-alpha warning (#21)

- ConformalAnomalyDetector.pvalue now uses jnp.searchsorted on the
  sorted calibration scores (O(log n)) instead of an O(n) masked sum;
  bit-identical results, vmaps cleanly through pvalue_batch.
- HDClassifier/HDAnomalyDetector now call check_is_fitted so calling
  predict/predict_proba/pvalue/score_samples/decision_function before
  fit raises sklearn.exceptions.NotFittedError.
- HDAnomalyDetector.fit warns when alpha is below the conformal
  resolution floor 1/(n_cal+1), where no point can ever be flagged.
- Align sklearn anomaly threshold (p <= alpha) with the core
  ConformalAnomalyDetector decision rule.

Squash-merge of cursor/anomaly-sklearn-robustness-795a (PR #21).

Co-authored-by: Rajdeep Singh <rlogger@users.noreply.github.com>

Author: Cursor Agent

bayes_hdc/anomaly.py

@@ -481,7 +481,12 @@ def pvalue(self, query_hv: jax.Array) -> jax.Array:
         """
         score = self.score(query_hv)
         n = self.n_calibration
-        ge_count = jnp.sum((self.calibration_scores >= score).astype(jnp.float32))
+        # ``calibration_scores`` is stored sorted ascending (see ``fit``), so the
+        # number of calibration scores >= ``score`` is an O(log n) binary search
+        # rather than an O(n) scan. ``side="left"`` returns the first index whose
+        # value is >= ``score``; everything from there to the end is the >= count.
+        # This is exact (identical to the masked sum) and vmaps cleanly.
+        ge_count = n - jnp.searchsorted(self.calibration_scores, score, side="left")
         return (1.0 + ge_count) / (n + 1.0)
 
     def pvalue_batch(self, queries: jax.Array) -> jax.Array:

bayes_hdc/sklearn.py

@@ -34,6 +34,7 @@
 
 from __future__ import annotations
 
+import warnings
 from typing import Any
 
 import jax
@@ -42,6 +43,7 @@
 
 try:
     from sklearn.base import BaseEstimator, ClassifierMixin, OutlierMixin
+    from sklearn.utils.validation import check_is_fitted
 except ImportError as exc:  # pragma: no cover - exercised only without sklearn
     raise ImportError(
         "bayes_hdc.sklearn requires scikit-learn. Install it with "
@@ -137,10 +139,12 @@ def _encode(self, X: Any) -> jax.Array:
         return self.encoder_.encode_batch(jnp.asarray(_as_f32(X)))
 
     def predict(self, X: Any) -> np.ndarray:
+        check_is_fitted(self, "classifier_")
         idx = np.asarray(self.classifier_.predict(self._encode(X)))
         return self.classes_[idx]
 
     def predict_proba(self, X: Any) -> np.ndarray:
+        check_is_fitted(self, "classifier_")
         return np.asarray(self.classifier_.predict_proba(self._encode(X)))
 
 
@@ -216,19 +220,38 @@ def fit(self, X: Any, y: Any = None) -> HDAnomalyDetector:
             k_neighbors=self.k_neighbors,
         ).fit(normal_hvs)
         self.detector_ = ConformalAnomalyDetector.create(scorer).fit(cal_hvs)
+
+        # Conformal p-values live on the grid {1/(n_cal+1), ..., 1}, so the
+        # smallest attainable p-value is 1/(n_cal+1). If ``alpha`` is below that
+        # floor, ``predict`` can never flag *any* point as an outlier (no matter
+        # how extreme) and the detector silently returns all-inliers. Warn loudly
+        # rather than letting users mistake "0 detections" for "0 anomalies".
+        n_cal = int(self.detector_.n_calibration)
+        floor = 1.0 / (n_cal + 1.0)
+        if self.alpha < floor:
+            warnings.warn(
+                f"alpha={self.alpha:g} is below the conformal resolution floor "
+                f"1/(n_calibration+1)={floor:g} (n_calibration={n_cal}). No point "
+                f"can be flagged as an outlier at this alpha. Increase the amount "
+                f"of fit data, raise calibration_fraction, or use a larger alpha "
+                f"(alpha >= {floor:g}).",
+                UserWarning,
+                stacklevel=2,
+            )
         return self
 
     def _encode(self, X: Any) -> jax.Array:
         return self.encoder_.encode_batch(jnp.asarray(_as_f32(X)))
 
     def pvalue(self, X: Any) -> np.ndarray:
         """Split-conformal p-values; small = anomalous."""
+        check_is_fitted(self, "detector_")
         return np.asarray(self.detector_.pvalue_batch(self._encode(X)))
 
     def predict(self, X: Any) -> np.ndarray:
         """+1 for inliers, -1 for outliers (scikit-learn convention)."""
         pvals = self.pvalue(X)
-        return np.where(pvals < self.alpha, -1, 1).astype(int)
+        return np.where(pvals <= self.alpha, -1, 1).astype(int)
 
     def score_samples(self, X: Any) -> np.ndarray:
         """Higher = more normal (scikit-learn convention).

tests/test_sklearn.py

@@ -5,12 +5,15 @@
 
 from __future__ import annotations
 
+import warnings
+
 import numpy as np
 import pytest
 
 pytest.importorskip("sklearn")
 
 from sklearn.datasets import make_blobs, make_classification  # noqa: E402
+from sklearn.exceptions import NotFittedError  # noqa: E402
 from sklearn.model_selection import cross_val_score  # noqa: E402
 from sklearn.pipeline import make_pipeline  # noqa: E402
 from sklearn.preprocessing import StandardScaler  # noqa: E402
@@ -172,3 +175,36 @@ def test_anomaly_get_params_roundtrip():
     assert det.get_params()["alpha"] == 0.05
     det.set_params(alpha=0.01)
     assert det.get_params()["alpha"] == 0.01
+
+
+def test_anomaly_warns_when_alpha_below_resolution_floor():
+    # n=20, calibration_fraction=0.3 -> n_cal=6 -> floor = 1/7 ~= 0.143.
+    # alpha=0.05 is below the floor, so no point can ever be flagged.
+    X_norm, _ = make_blobs(n_samples=20, n_features=4, centers=1, random_state=0)
+    rng = np.random.default_rng(0)
+    X_out = rng.normal(loc=15.0, size=(10, 4))
+    with pytest.warns(UserWarning, match="conformal resolution floor"):
+        det = HDAnomalyDetector(alpha=0.05, dimensions=2000, random_state=0).fit(X_norm)
+    # The warning's premise: even extreme outliers cannot be flagged at this alpha.
+    assert (det.predict(X_out) == -1).sum() == 0
+
+
+def test_anomaly_no_warning_when_alpha_attainable():
+    X_norm, _ = make_blobs(n_samples=300, n_features=4, centers=1, random_state=0)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")  # any UserWarning would fail the test
+        HDAnomalyDetector(alpha=0.05, dimensions=2000, random_state=0).fit(X_norm)
+
+
+def test_estimators_raise_not_fitted_before_fit():
+    clf = HDClassifier(dimensions=1000)
+    det = HDAnomalyDetector(dimensions=1000)
+    X = np.zeros((3, 4), dtype=np.float32)
+    with pytest.raises(NotFittedError):
+        clf.predict(X)
+    with pytest.raises(NotFittedError):
+        clf.predict_proba(X)
+    with pytest.raises(NotFittedError):
+        det.predict(X)
+    with pytest.raises(NotFittedError):
+        det.pvalue(X)