diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 0d08e261a..0a51ccf3c 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -9,3 +9,7 @@ updates: directory: "/" # Location of package manifests schedule: interval: "weekly" + - package-ecosystem: "uv" + directory: "/" # Location of pyproject.toml / uv.lock + schedule: + interval: "weekly" diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 000000000..0d1653aa7 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,96 @@ +# Changelog + +All notable changes to this project are documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), +and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [Unreleased] + +## [0.3.0] - 2026-08-04 + +### Added + +- `cyclops.monitor`: `Detector.detect_shift_by_subgroup()` runs the fitted + drift tester independently on each subgroup of a target dataset (e.g. + age band, sex, hospital site, defined via the existing `SliceSpec`), + instead of only testing the aggregate population - a model can look + stable overall while drifting badly for a specific subgroup, which + matters for health-equity-aware monitoring. Includes Bonferroni + correction across subgroups and a minimum-sample-size guard. +- `cyclops.monitor`: `DCTester.explain_shift()` explains a detected shift + using SHAP on the domain classifier trained by `tester_method="classifier"`, + returning features ranked by how strongly they indicate a sample belongs + to the shifted distribution. +- `cyclops.evaluate.metrics.experimental`: added calibration metrics - + `BinaryBrierScore`/`MulticlassBrierScore` and `BinaryCalibrationError` + (Expected Calibration Error / Maximum Calibration Error), plus their + functional counterparts. Discrimination metrics like AUROC don't tell + you whether a predicted probability can be trusted at face value, which + matters for clinical risk scores that are often acted on directly. + +### Fixed + +- `cyclops.monitor`: `errorfill()` crashed on the default `color=None` because + matplotlib removed `ax._get_lines.prop_cycler`. +- `cyclops.monitor`: `TSTester.test_shift()` mutated `p_val_threshold` in + place on every call, so the Bonferroni correction compounded across + repeated calls (e.g. in `Detector`'s sweep loops) instead of being + computed fresh each time; also fixed an `UnboundLocalError` when the + input isn't a plain `np.ndarray`. +- `cyclops.monitor`: `ContextMMDWrapper` was missing the + `preprocess_at_init` argument in its positional argument list to + alibi-detect's `ContextMMDDrift`, silently shifting every later + argument by one slot. `ContextMMDWrapper` and `LKWrapper` now pass + keyword arguments so future alibi-detect signature changes fail loudly + instead of silently misaligning. +- `cyclops.monitor`: `Reductor` raised `TypeError` when torchvision wasn't + installed and `transforms` was passed, because `isinstance(transforms, + Compose)` was called with `Compose is None`. +- `cyclops.monitor`: removed `plot_label_distribution`, an unreachable, + untested, and uncalled function with a use-before-assignment bug. +- `cyclops.report`: `ModelCardReport.export()` raised `IndexError` when no + `PerformanceMetric` had been logged. +- `cyclops.report`: `_process_metric_name()` raised `UnboundLocalError` + for any metric `type` not prefixed with `Binary`/`Multiclass`/ + `Multilabel` (e.g. a custom metric name). +- `cyclops.report`: `export()`'s default output filename was a static + `model_card.html`/`.json`, so repeated calls into the same + `output_dir` silently overwrote prior reports and broke trend/history + comparisons. The default filename is now timestamped per call. +- `cyclops.report`: `Citation.content` (raw BibTeX text) was rendered + with Jinja's `|safe` filter, bypassing autoescaping for no reason. +- `cyclops.utils`: `exchange_extension()` dropped the filename entirely + for paths with no existing extension (e.g. `"myfile"` -> `".csv"` + instead of `"myfile.csv"`). +- `cyclops.models`: `MLPModel` (the packaged `"mlp_pt"` model) could not + be constructed at all - an activation class was inserted into + `nn.Sequential` instead of an instance, the first hidden layer was + double-wrapped in a list, and the loop connecting hidden layers used + the wrong input dimension. The packaged `mlp_pt.yaml` config also + passed a nonexistent `layer_dim` argument left over from the RNN/GRU/ + LSTM configs. +- `cyclops.data`: `SliceSpec`'s datetime `day` component filter + (`filter_datetime`) silently matched on year instead of day of month. + +### Changed + +- `cyclops.monitor.utils`: removed ~470 lines of unused temporal-modeling + scaffolding (`Data`, `get_data`, `run_model`, `get_serving_data`, + `scale`, `daterange`, `get_obj_from_str`, `load_model`/`save_model`, + `print_metrics_binary`, `load_ckp`, `get_device`, `get_temporal_model`, + `Loader`, and a stray `__main__` demo block) that was neither exported, + imported elsewhere in the repo, nor tested. + +### CI / infra + +- Added a `uv` ecosystem entry to Dependabot so `pyproject.toml`/ + `uv.lock` dependencies get automated update PRs. +- Fixed the README's "integration tests" badge, which linked to a + workflow file that no longer exists; replaced with the (existing, + previously unlinked) unit tests badge. +- Expanded `CONTRIBUTING.md` with environment setup, test-running, and + repository layout sections. + +[Unreleased]: https://github.com/VectorInstitute/cyclops/compare/v0.3.0...HEAD +[0.3.0]: https://github.com/VectorInstitute/cyclops/compare/v0.2.12...v0.3.0 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5e27a2db8..fc3123dcd 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -5,14 +5,49 @@ Thanks for your interest in contributing to cyclops! To submit PRs, please fill out the PR template along with the PR. If the PR fixes an issue, don't forget to link the PR to the issue! -## Pre-commit hooks +## Setting up your environment -Once the python virtual environment is setup, you can run pre-commit hooks using: +cyclops uses [uv](https://docs.astral.sh/uv/getting-started/installation/) to +manage dependencies. Once uv is installed, set up a development environment +with the test dependency group and activate it: ```bash -pre-commit run --all-files +uv sync --group test +source .venv/bin/activate +``` + +Some modules have optional dependencies (e.g. `torch`, `xgboost`, `monai`, +`alibi-detect`) that aren't installed by default - see the table in +[README.md](README.md) for the full list of extras. To work on a module that +needs one of these, install the matching extra, e.g.: + +```bash +uv sync --group test --extra alibi-detect +``` + +Install the pre-commit hooks so code-style issues are caught before you push: + +```bash +pre-commit install +``` + +## Running tests + +Run the unit test suite with: + +```bash +python -m pytest -m "not integration_test" ``` +Tests marked `@pytest.mark.integration_test` require external +infrastructure (e.g. a live database via +[cycquery](https://github.com/VectorInstitute/cycquery)) that isn't +available in a plain checkout, which is why they're excluded above and not +run in CI. Only run them locally if you have that infrastructure set up. + +Pass `-k ` to scope a run to a subset of tests, and +`--cov=cyclops` to see coverage for the code you changed. + ## Coding guidelines For code style, we recommend the [PEP 8 style guide](https://peps.python.org/pep-0008/). @@ -24,3 +59,27 @@ analysis. Ruff checks various rules including [flake8](https://docs.astral.sh/ru Last but not the least, we use type hints in our code which is then checked using [mypy](https://mypy.readthedocs.io/en/stable/). + +You can run all pre-commit checks (ruff, ruff-format, mypy, notebook +stripping) against the whole repository at once with: + +```bash +pre-commit run --all-files +``` + +## Repository layout + +- `cyclops/data` - dataset construction, loading, and slicing +- `cyclops/models` - scikit-learn and PyTorch model wrappers and implementations +- `cyclops/tasks` - task formulations (e.g. binary/multi-label classification) tying data and models together +- `cyclops/evaluate` - metrics and fairness evaluation for clinical prediction tasks +- `cyclops/monitor` - dataset shift / drift detection for deployed models +- `cyclops/report` - model report card generation +- `cyclops/utils` - small shared utilities used across the other modules + +Each module has a corresponding test package under `tests/cyclops/`. + +## Code of Conduct + +Participation in this project is governed by our +[Code of Conduct](CODE_OF_CONDUCT.md). diff --git a/README.md b/README.md index e11c6c18c..cbff2d5b3 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![PyPI](https://img.shields.io/pypi/v/pycyclops)](https://pypi.org/project/pycyclops) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/pycyclops) [![code checks](https://github.com/VectorInstitute/cyclops/actions/workflows/code_checks.yml/badge.svg)](https://github.com/VectorInstitute/cyclops/actions/workflows/code_checks.yml) -[![integration tests](https://github.com/VectorInstitute/cyclops/actions/workflows/integration_tests.yml/badge.svg)](https://github.com/VectorInstitute/cyclops/actions/workflows/integration_tests.yml) +[![unit tests](https://github.com/VectorInstitute/cyclops/actions/workflows/unit_tests.yml/badge.svg)](https://github.com/VectorInstitute/cyclops/actions/workflows/unit_tests.yml) [![docs](https://github.com/VectorInstitute/cyclops/actions/workflows/docs.yml/badge.svg)](https://github.com/VectorInstitute/cyclops/actions/workflows/docs.yml) [![codecov](https://codecov.io/gh/VectorInstitute/cyclops/branch/main/graph/badge.svg)](https://codecov.io/gh/VectorInstitute/cyclops) [![docker](https://github.com/VectorInstitute/cyclops/actions/workflows/docker.yml/badge.svg)](https://hub.docker.com/r/vectorinstitute/cyclops) diff --git a/cyclops/data/slicer.py b/cyclops/data/slicer.py index 92e74563a..f0282bd1c 100644 --- a/cyclops/data/slicer.py +++ b/cyclops/data/slicer.py @@ -768,7 +768,7 @@ def _apply_mask( months = pc.month(example_values) mask = _apply_mask(months, month, mask) if day is not None: - days = pc.year(example_values) + days = pc.day(example_values) mask = _apply_mask(days, day, mask) if hour is not None: hours = pc.hour(example_values) diff --git a/cyclops/evaluate/metrics/experimental/__init__.py b/cyclops/evaluate/metrics/experimental/__init__.py index f80f1c431..7b7d427ba 100644 --- a/cyclops/evaluate/metrics/experimental/__init__.py +++ b/cyclops/evaluate/metrics/experimental/__init__.py @@ -15,6 +15,13 @@ MulticlassAveragePrecision, MultilabelAveragePrecision, ) +from cyclops.evaluate.metrics.experimental.brier_score import ( + BinaryBrierScore, + MulticlassBrierScore, +) +from cyclops.evaluate.metrics.experimental.calibration_error import ( + BinaryCalibrationError, +) from cyclops.evaluate.metrics.experimental.confusion_matrix import ( BinaryConfusionMatrix, MulticlassConfusionMatrix, @@ -95,6 +102,9 @@ "BinaryAveragePrecision", "MulticlassAveragePrecision", "MultilabelAveragePrecision", + "BinaryBrierScore", + "MulticlassBrierScore", + "BinaryCalibrationError", "BinaryConfusionMatrix", "MulticlassConfusionMatrix", "MultilabelConfusionMatrix", diff --git a/cyclops/evaluate/metrics/experimental/brier_score.py b/cyclops/evaluate/metrics/experimental/brier_score.py new file mode 100644 index 000000000..292c6e4c8 --- /dev/null +++ b/cyclops/evaluate/metrics/experimental/brier_score.py @@ -0,0 +1,163 @@ +"""Brier score metric.""" + +from typing import Any, Optional + +from cyclops.evaluate.metrics.experimental.functional.brier_score import ( + _binary_brier_score_compute, + _binary_brier_score_format_arrays, + _binary_brier_score_update, + _binary_brier_score_validate_args, + _binary_brier_score_validate_arrays, + _multiclass_brier_score_format_arrays, + _multiclass_brier_score_update, + _multiclass_brier_score_validate_args, + _multiclass_brier_score_validate_arrays, +) +from cyclops.evaluate.metrics.experimental.metric import Metric +from cyclops.evaluate.metrics.experimental.utils.types import Array + + +class BinaryBrierScore(Metric): + """Brier score for binary classification tasks. + + The Brier score is the mean squared error between predicted + probabilities and the (binary) target, and is a proper scoring rule + for probabilistic predictions - it rewards models whose predicted + probabilities are well-calibrated, not just well-ranked, which matters + for clinical risk scores that are acted on directly. + + Parameters + ---------- + ignore_index : int, optional, default=None + Values in the target array to ignore when computing the metric. + **kwargs : Any + Additional keyword arguments common to all metrics. + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental import BinaryBrierScore + >>> target = anp.asarray([0, 1, 1, 0]) + >>> preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + >>> metric = BinaryBrierScore() + >>> metric(target, preds) + Array(0.0375, dtype=float32) + + """ + + name: str = "Brier Score" + + def __init__(self, ignore_index: Optional[int] = None, **kwargs: Any) -> None: + super().__init__(**kwargs) + _binary_brier_score_validate_args(ignore_index=ignore_index) + self.ignore_index = ignore_index + + self.add_state_default_factory( + "sum_squared_error", + lambda xp: xp.asarray(0.0, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + self.add_state_default_factory( + "num_obs", + lambda xp: xp.asarray(0.0, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + + def _update_state(self, target: Array, preds: Array) -> None: + """Update the state of the metric.""" + xp = _binary_brier_score_validate_arrays( + target, + preds, + ignore_index=self.ignore_index, + ) + target, preds = _binary_brier_score_format_arrays( + target, + preds, + self.ignore_index, + xp=xp, + ) + sum_squared_error, num_obs = _binary_brier_score_update(target, preds) + self.sum_squared_error += sum_squared_error # type: ignore + self.num_obs += num_obs # type: ignore + + def _compute_metric(self) -> Array: + """Compute the binary Brier score.""" + return _binary_brier_score_compute( + self.sum_squared_error, # type: ignore + self.num_obs, # type: ignore + ) + + +class MulticlassBrierScore(Metric): + """Brier score for multiclass classification tasks. + + Computed as the mean squared error between the predicted probability + vector for each sample and the one-hot encoded target. + + Parameters + ---------- + num_classes : int + The number of classes in the classification task. + ignore_index : int, optional, default=None + Values in the target array to ignore when computing the metric. + **kwargs : Any + Additional keyword arguments common to all metrics. + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental import MulticlassBrierScore + >>> target = anp.asarray([0, 1, 2]) + >>> preds = anp.asarray( + ... [[0.7, 0.2, 0.1], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]], + ... ) + >>> metric = MulticlassBrierScore(num_classes=3) + >>> metric(target, preds) + Array(0.14666666, dtype=float32) + + """ + + name: str = "Brier Score" + + def __init__( + self, + num_classes: int, + ignore_index: Optional[int] = None, + **kwargs: Any, + ) -> None: + super().__init__(**kwargs) + _multiclass_brier_score_validate_args(num_classes, ignore_index=ignore_index) + self.num_classes = num_classes + self.ignore_index = ignore_index + + self.add_state_default_factory( + "sum_squared_error", + lambda xp: xp.asarray(0.0, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + self.add_state_default_factory( + "num_obs", + lambda xp: xp.asarray(0.0, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + + def _update_state(self, target: Array, preds: Array) -> None: + """Update the state of the metric.""" + xp = _multiclass_brier_score_validate_arrays(target, preds, self.num_classes) + target, preds = _multiclass_brier_score_format_arrays( + target, + preds, + self.ignore_index, + self.num_classes, + xp=xp, + ) + sum_squared_error, num_obs = _multiclass_brier_score_update(target, preds) + self.sum_squared_error += sum_squared_error # type: ignore + self.num_obs += num_obs # type: ignore + + def _compute_metric(self) -> Array: + """Compute the multiclass Brier score.""" + return _binary_brier_score_compute( + self.sum_squared_error, # type: ignore + self.num_obs, # type: ignore + ) diff --git a/cyclops/evaluate/metrics/experimental/calibration_error.py b/cyclops/evaluate/metrics/experimental/calibration_error.py new file mode 100644 index 000000000..c22ac8230 --- /dev/null +++ b/cyclops/evaluate/metrics/experimental/calibration_error.py @@ -0,0 +1,122 @@ +"""Calibration error metric.""" + +from typing import Any, Literal, Optional + +from cyclops.evaluate.metrics.experimental.functional.brier_score import ( + _binary_brier_score_format_arrays, +) +from cyclops.evaluate.metrics.experimental.functional.calibration_error import ( + _binary_calibration_error_compute, + _binary_calibration_error_update, + _binary_calibration_error_validate_args, + _binary_calibration_error_validate_arrays, +) +from cyclops.evaluate.metrics.experimental.metric import Metric +from cyclops.evaluate.metrics.experimental.utils.types import Array + + +class BinaryCalibrationError(Metric): + """Calibration error for binary classification tasks. + + Groups predicted probabilities into equal-width bins and measures, + within each bin, the gap between the average predicted probability + (confidence) and the observed event rate (accuracy). The default + `"l1"` norm gives the Expected Calibration Error (ECE), the most + commonly reported calibration metric. + + A well-calibrated clinical risk model should have a low calibration + error: among patients given, say, a 30% predicted risk, roughly 30% + should actually experience the event. This matters even for models + with good discrimination (e.g. high AUROC), since discrimination + alone doesn't guarantee predicted probabilities can be trusted at + face value - which is often how clinical risk scores are actually + used. + + Parameters + ---------- + n_bins : int, optional, default=15 + Number of equal-width bins to group predicted probabilities into. + norm : {'l1', 'l2', 'max'}, optional, default='l1' + Norm used to aggregate the per-bin calibration gaps. `'l1'` gives + the Expected Calibration Error (ECE), `'max'` gives the Maximum + Calibration Error (MCE). + ignore_index : int, optional, default=None + Values in the target array to ignore when computing the metric. + **kwargs : Any + Additional keyword arguments common to all metrics. + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental import BinaryCalibrationError + >>> target = anp.asarray([0, 1, 1, 0]) + >>> preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + >>> metric = BinaryCalibrationError(n_bins=2) + >>> metric(target, preds) + Array(0.17499998, dtype=float32) + + """ + + name: str = "Calibration Error" + + def __init__( + self, + n_bins: int = 15, + norm: Literal["l1", "l2", "max"] = "l1", + ignore_index: Optional[int] = None, + **kwargs: Any, + ) -> None: + super().__init__(**kwargs) + _binary_calibration_error_validate_args( + n_bins=n_bins, + norm=norm, + ignore_index=ignore_index, + ) + self.n_bins = n_bins + self.norm = norm + self.ignore_index = ignore_index + + self.add_state_default_factory( + "bin_confidence_sums", + lambda xp: xp.zeros(n_bins, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + self.add_state_default_factory( + "bin_correct_sums", + lambda xp: xp.zeros(n_bins, dtype=xp.float32, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + self.add_state_default_factory( + "bin_counts", + lambda xp: xp.zeros(n_bins, dtype=xp.int64, device=self.device), # type: ignore + dist_reduce_fn="sum", + ) + + def _update_state(self, target: Array, preds: Array) -> None: + """Update the state of the metric.""" + xp = _binary_calibration_error_validate_arrays( + target, + preds, + ignore_index=self.ignore_index, + ) + target, preds = _binary_brier_score_format_arrays( + target, + preds, + self.ignore_index, + xp=xp, + ) + bin_confidence_sums, bin_correct_sums, bin_counts = ( + _binary_calibration_error_update(target, preds, self.n_bins, xp=xp) + ) + self.bin_confidence_sums += bin_confidence_sums # type: ignore + self.bin_correct_sums += bin_correct_sums # type: ignore + self.bin_counts += bin_counts # type: ignore + + def _compute_metric(self) -> Array: + """Compute the binary calibration error.""" + return _binary_calibration_error_compute( + self.bin_confidence_sums, # type: ignore + self.bin_correct_sums, # type: ignore + self.bin_counts, # type: ignore + norm=self.norm, + ) diff --git a/cyclops/evaluate/metrics/experimental/functional/__init__.py b/cyclops/evaluate/metrics/experimental/functional/__init__.py index 09251c4bb..fe4e31e56 100644 --- a/cyclops/evaluate/metrics/experimental/functional/__init__.py +++ b/cyclops/evaluate/metrics/experimental/functional/__init__.py @@ -15,6 +15,13 @@ multiclass_average_precision, multilabel_average_precision, ) +from cyclops.evaluate.metrics.experimental.functional.brier_score import ( + binary_brier_score, + multiclass_brier_score, +) +from cyclops.evaluate.metrics.experimental.functional.calibration_error import ( + binary_calibration_error, +) from cyclops.evaluate.metrics.experimental.functional.confusion_matrix import ( binary_confusion_matrix, multiclass_confusion_matrix, @@ -95,6 +102,9 @@ "binary_average_precision", "multiclass_average_precision", "multilabel_average_precision", + "binary_brier_score", + "multiclass_brier_score", + "binary_calibration_error", "binary_confusion_matrix", "multiclass_confusion_matrix", "multilabel_confusion_matrix", diff --git a/cyclops/evaluate/metrics/experimental/functional/brier_score.py b/cyclops/evaluate/metrics/experimental/functional/brier_score.py new file mode 100644 index 000000000..fec7f8083 --- /dev/null +++ b/cyclops/evaluate/metrics/experimental/functional/brier_score.py @@ -0,0 +1,282 @@ +"""Functional interface for the Brier score metric.""" + +from types import ModuleType +from typing import Optional, Tuple, Union + +import array_api_compat as apc + +from cyclops.evaluate.metrics.experimental.functional._stat_scores import ( + _binary_stat_scores_validate_args, + _binary_stat_scores_validate_arrays, +) +from cyclops.evaluate.metrics.experimental.utils.ops import ( + _to_one_hot, + flatten, + remove_ignore_index, + sigmoid, + squeeze_all, + to_int, +) +from cyclops.evaluate.metrics.experimental.utils.types import Array +from cyclops.evaluate.metrics.experimental.utils.validation import ( + _basic_input_array_checks, + is_floating_point, +) + + +def _binary_brier_score_validate_args(ignore_index: Optional[int] = None) -> None: + """Validate arguments for binary Brier score computation.""" + _binary_stat_scores_validate_args(threshold=0.5, ignore_index=ignore_index) + + +def _binary_brier_score_validate_arrays( + target: Array, + preds: Array, + ignore_index: Optional[int] = None, +) -> ModuleType: + """Validate `target` and `preds` for binary Brier score computation.""" + return _binary_stat_scores_validate_arrays(target, preds, ignore_index=ignore_index) + + +def _binary_brier_score_format_arrays( + target: Array, + preds: Array, + ignore_index: Optional[int], + *, + xp: ModuleType, +) -> Tuple[Array, Array]: + """Format `target` and `preds` for binary Brier score computation. + + Unlike the stat-scores formatting used for other binary classification + metrics, `preds` is kept as a continuous probability (not thresholded + into a hard label), since the Brier score is a proper scoring rule + computed directly on predicted probabilities. + """ + target = flatten(target) + preds = flatten(preds) + + if ignore_index is not None: + target, preds = remove_ignore_index(target, preds, ignore_index=ignore_index) + + if not is_floating_point(preds): + preds = xp.astype(preds, xp.float32) + elif not xp.all(to_int(preds >= 0) * to_int(preds <= 1)): # preds are logits + preds = sigmoid(preds) + + return xp.astype(target, preds.dtype), preds + + +def _binary_brier_score_update(target: Array, preds: Array) -> Tuple[Array, int]: + """Update and return variables required to compute the binary Brier score.""" + xp = apc.array_namespace(target, preds) + diff = preds - target + sum_squared_error = xp.sum(diff * diff, dtype=xp.float32) + return sum_squared_error, target.shape[0] + + +def _binary_brier_score_compute( + sum_squared_error: Array, + num_obs: Union[int, Array], +) -> Array: + """Compute the binary Brier score from the accumulated state.""" + return squeeze_all(sum_squared_error / num_obs) + + +def binary_brier_score( + target: Array, + preds: Array, + ignore_index: Optional[int] = None, +) -> Array: + """Compute the Brier score for binary classification tasks. + + The Brier score is the mean squared error between predicted + probabilities and the (binary) target, and is a proper scoring rule + for probabilistic predictions - it rewards models whose predicted + probabilities are well-calibrated, not just well-ranked, which matters + for clinical risk scores that are acted on directly (e.g. a 30% + predicted mortality risk should correspond to an observed 30% event + rate). + + Parameters + ---------- + target : Array + Ground truth binary labels (0 or 1). + preds : Array + Predicted probabilities (or logits, which are converted to + probabilities via the sigmoid function) of the positive class. + ignore_index : int, optional, default=None + Values in `target` to ignore when computing the metric. + + Returns + ------- + Array + The Brier score, in the range [0, 1] (lower is better). + + Raises + ------ + TypeError + If `target` or `preds` is not an array object that is compatible + with the Python array API standard. + ValueError + If `target` or `preds` is empty, not a numeric array, or if the + shape of `target` and `preds` are not the same. + RuntimeError + If `target` contains values other than 0, 1 (and `ignore_index`, + if specified). + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental.functional import ( + ... binary_brier_score, + ... ) + >>> target = anp.asarray([0, 1, 1, 0]) + >>> preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + >>> binary_brier_score(target, preds) + Array(0.0375, dtype=float32) + + """ + _binary_brier_score_validate_args(ignore_index=ignore_index) + xp = _binary_brier_score_validate_arrays(target, preds, ignore_index=ignore_index) + target, preds = _binary_brier_score_format_arrays( + target, + preds, + ignore_index, + xp=xp, + ) + sum_squared_error, num_obs = _binary_brier_score_update(target, preds) + return _binary_brier_score_compute(sum_squared_error, num_obs) + + +def _multiclass_brier_score_validate_args( + num_classes: int, + ignore_index: Optional[int] = None, +) -> None: + """Validate arguments for multiclass Brier score computation.""" + if not isinstance(num_classes, int) or num_classes < 2: + raise ValueError( + f"Expected argument `num_classes` to be an integer larger than 1, " + f"but got {num_classes}", + ) + if ignore_index is not None and not isinstance(ignore_index, int): + raise ValueError( + "Expected argument `ignore_index` to either be `None` or an integer, " + f"but got {ignore_index}", + ) + + +def _multiclass_brier_score_validate_arrays( + target: Array, + preds: Array, + num_classes: int, +) -> ModuleType: + """Validate `target` and `preds` for multiclass Brier score computation.""" + _basic_input_array_checks(target, preds) + xp = apc.array_namespace(target, preds) + + if not (preds.ndim == target.ndim + 1 and is_floating_point(preds)): + raise ValueError( + "Expected `preds` to be a floating point array with one more " + "dimension than `target`, containing predicted probabilities for " + f"each of the {num_classes} classes. Got `preds` with shape " + f"{preds.shape} and `target` with shape {target.shape}.", + ) + if preds.shape[-1] != num_classes: + raise ValueError( + "Expected the last dimension of `preds` to be equal to " + f"`num_classes` ({num_classes}), but got {preds.shape[-1]}.", + ) + return xp # type: ignore[no-any-return] + + +def _multiclass_brier_score_format_arrays( + target: Array, + preds: Array, + ignore_index: Optional[int], + num_classes: int, + *, + xp: ModuleType, +) -> Tuple[Array, Array]: + """Format `target` and `preds` for multiclass Brier score computation.""" + target = flatten(target) + preds = xp.reshape(preds, (-1, num_classes)) + + if ignore_index is not None: + target, preds = remove_ignore_index(target, preds, ignore_index=ignore_index) + + target = _to_one_hot(to_int(target), num_classes=num_classes) + return xp.astype(target, xp.float32), preds + + +def _multiclass_brier_score_update(target: Array, preds: Array) -> Tuple[Array, int]: + """Update and return variables required to compute the multiclass Brier score.""" + xp = apc.array_namespace(target, preds) + diff = preds - target + sum_squared_error = xp.sum(xp.sum(diff * diff, axis=-1), dtype=xp.float32) + return sum_squared_error, target.shape[0] + + +def multiclass_brier_score( + target: Array, + preds: Array, + num_classes: int, + ignore_index: Optional[int] = None, +) -> Array: + """Compute the Brier score for multiclass classification tasks. + + Computed as the mean squared error between the predicted probability + vector for each sample and the one-hot encoded target. + + Parameters + ---------- + target : Array + Ground truth class labels, shape `(N, ...)`. + preds : Array + Predicted probabilities for each class, shape `(N, C, ...)`. Rows + are expected to sum to 1. + num_classes : int + Number of classes. + ignore_index : int, optional, default=None + Values in `target` to ignore when computing the metric. + + Returns + ------- + Array + The (multiclass) Brier score, in the range [0, 2] (lower is + better). + + Raises + ------ + TypeError + If `target` or `preds` is not an array object that is compatible + with the Python array API standard. + ValueError + If `num_classes` is not an integer larger than 1, if `preds` does + not have one more dimension than `target`, or if the size of the + last dimension of `preds` is not equal to `num_classes`. + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental.functional import ( + ... multiclass_brier_score, + ... ) + >>> target = anp.asarray([0, 1, 2]) + >>> preds = anp.asarray( + ... [[0.7, 0.2, 0.1], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]], + ... ) + >>> multiclass_brier_score(target, preds, num_classes=3) + Array(0.14666666, dtype=float32) + + """ + _multiclass_brier_score_validate_args(num_classes, ignore_index=ignore_index) + xp = _multiclass_brier_score_validate_arrays(target, preds, num_classes) + target, preds = _multiclass_brier_score_format_arrays( + target, + preds, + ignore_index, + num_classes, + xp=xp, + ) + sum_squared_error, num_obs = _multiclass_brier_score_update(target, preds) + return _binary_brier_score_compute(sum_squared_error, num_obs) diff --git a/cyclops/evaluate/metrics/experimental/functional/calibration_error.py b/cyclops/evaluate/metrics/experimental/functional/calibration_error.py new file mode 100644 index 000000000..80902f2da --- /dev/null +++ b/cyclops/evaluate/metrics/experimental/functional/calibration_error.py @@ -0,0 +1,213 @@ +"""Functional interface for the (binary) calibration error metric.""" + +from types import ModuleType +from typing import Literal, Optional, Tuple + +import array_api_compat as apc + +from cyclops.evaluate.metrics.experimental.functional._stat_scores import ( + _binary_stat_scores_validate_args, + _binary_stat_scores_validate_arrays, +) +from cyclops.evaluate.metrics.experimental.functional.brier_score import ( + _binary_brier_score_format_arrays, +) +from cyclops.evaluate.metrics.experimental.utils.ops import ( + bincount, + safe_divide, + to_int, +) +from cyclops.evaluate.metrics.experimental.utils.types import Array + + +_ALLOWED_NORMS = ("l1", "l2", "max") + + +def _binary_calibration_error_validate_args( + n_bins: int = 15, + norm: Literal["l1", "l2", "max"] = "l1", + ignore_index: Optional[int] = None, +) -> None: + """Validate arguments for binary calibration error computation.""" + if not isinstance(n_bins, int) or n_bins < 1: + raise ValueError( + f"Expected argument `n_bins` to be a positive integer, but got {n_bins}", + ) + if norm not in _ALLOWED_NORMS: + raise ValueError( + f"Expected argument `norm` to be one of {_ALLOWED_NORMS}, but got {norm}", + ) + _binary_stat_scores_validate_args(threshold=0.5, ignore_index=ignore_index) + + +def _binary_calibration_error_validate_arrays( + target: Array, + preds: Array, + ignore_index: Optional[int] = None, +) -> ModuleType: + """Validate `target` and `preds` for binary calibration error computation.""" + return _binary_stat_scores_validate_arrays(target, preds, ignore_index=ignore_index) + + +def _binary_calibration_error_update( + target: Array, + preds: Array, + n_bins: int, + *, + xp: ModuleType, +) -> Tuple[Array, Array, Array]: + """Compute per-bin confidence sum, correctness sum, and count.""" + bin_ids = to_int(xp.floor(preds * n_bins)) + # `preds == 1.0` falls in its own out-of-range bin; fold it into the last one + bin_ids = xp.where( + bin_ids >= n_bins, + xp.asarray(n_bins - 1, dtype=bin_ids.dtype, device=apc.device(bin_ids)), + bin_ids, + ) + + bin_confidence_sums = xp.astype( + bincount(bin_ids, weights=preds, minlength=n_bins), + xp.float32, + ) + bin_correct_sums = xp.astype( + bincount(bin_ids, weights=target, minlength=n_bins), + xp.float32, + ) + bin_counts = bincount(bin_ids, minlength=n_bins) + return bin_confidence_sums, bin_correct_sums, bin_counts + + +def _binary_calibration_error_compute( + bin_confidence_sums: Array, + bin_correct_sums: Array, + bin_counts: Array, + norm: Literal["l1", "l2", "max"] = "l1", +) -> Array: + """Compute the binary calibration error from the accumulated per-bin state.""" + xp = apc.array_namespace(bin_confidence_sums, bin_correct_sums, bin_counts) + bin_counts = xp.astype(bin_counts, xp.float32) + bin_confidence_sums = xp.astype(bin_confidence_sums, xp.float32) + bin_correct_sums = xp.astype(bin_correct_sums, xp.float32) + + avg_confidence = safe_divide(bin_confidence_sums, bin_counts) + avg_accuracy = safe_divide(bin_correct_sums, bin_counts) + gaps = xp.abs(avg_confidence - avg_accuracy) + + if norm == "max": + return xp.astype(xp.max(gaps), xp.float32) # type: ignore[no-any-return] + + bin_weights = safe_divide( + bin_counts, + xp.sum(bin_counts, dtype=xp.float32), + ) + if norm == "l2": + return xp.astype( # type: ignore[no-any-return] + xp.sqrt(xp.sum((gaps**2) * bin_weights, dtype=xp.float32)), + xp.float32, + ) + # l1, i.e. the "expected calibration error" (ECE) + return xp.sum(gaps * bin_weights, dtype=xp.float32) # type: ignore[no-any-return] + + +def binary_calibration_error( + target: Array, + preds: Array, + n_bins: int = 15, + norm: Literal["l1", "l2", "max"] = "l1", + ignore_index: Optional[int] = None, +) -> Array: + """Compute the calibration error for binary classification tasks. + + Groups predicted probabilities into `n_bins` equal-width bins and + measures, within each bin, the gap between the average predicted + probability (confidence) and the observed event rate (accuracy). The + `"l1"` norm (the default) gives the Expected Calibration Error (ECE), + the most commonly reported calibration metric. + + A well-calibrated clinical risk model should have a low calibration + error: among patients given, say, a 30% predicted risk, roughly 30% + should actually experience the event. This matters even for models + with good discrimination (e.g. high AUROC), since discrimination alone + doesn't guarantee that predicted probabilities can be trusted at face + value - which is often how clinical risk scores are actually used. + + Parameters + ---------- + target : Array + Ground truth binary labels (0 or 1). + preds : Array + Predicted probabilities (or logits, which are converted to + probabilities via the sigmoid function) of the positive class. + n_bins : int, optional, default=15 + Number of equal-width bins to group predicted probabilities into. + norm : {'l1', 'l2', 'max'}, optional, default='l1' + Norm used to aggregate the per-bin calibration gaps: + + - `'l1'`: the (sample-size-)weighted average absolute gap, i.e. + the Expected Calibration Error (ECE). + - `'l2'`: the (sample-size-)weighted root mean square gap. + - `'max'`: the largest gap across bins, i.e. the Maximum + Calibration Error (MCE). + ignore_index : int, optional, default=None + Values in `target` to ignore when computing the metric. + + Returns + ------- + Array + The calibration error, in the range [0, 1] (lower is better). + + Raises + ------ + TypeError + If `target` or `preds` is not an array object that is compatible + with the Python array API standard. + ValueError + If `n_bins` is not a positive integer, if `norm` is not one of + `'l1'`, `'l2'`, `'max'`, or if `target` or `preds` is empty, not a + numeric array, or not the same shape. + RuntimeError + If `target` contains values other than 0, 1 (and `ignore_index`, + if specified). + + Examples + -------- + >>> import numpy.array_api as anp + >>> from cyclops.evaluate.metrics.experimental.functional import ( + ... binary_calibration_error, + ... ) + >>> target = anp.asarray([0, 1, 1, 0]) + >>> preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + >>> binary_calibration_error(target, preds, n_bins=2) + Array(0.17499998, dtype=float32) + + """ + _binary_calibration_error_validate_args( + n_bins=n_bins, + norm=norm, + ignore_index=ignore_index, + ) + xp = _binary_calibration_error_validate_arrays( + target, + preds, + ignore_index=ignore_index, + ) + target, preds = _binary_brier_score_format_arrays( + target, + preds, + ignore_index, + xp=xp, + ) + bin_confidence_sums, bin_correct_sums, bin_counts = ( + _binary_calibration_error_update( + target, + preds, + n_bins, + xp=xp, + ) + ) + return _binary_calibration_error_compute( + bin_confidence_sums, + bin_correct_sums, + bin_counts, + norm=norm, + ) diff --git a/cyclops/models/configs/mlp_pt.yaml b/cyclops/models/configs/mlp_pt.yaml index b23d54d8b..d54cd460d 100644 --- a/cyclops/models/configs/mlp_pt.yaml +++ b/cyclops/models/configs/mlp_pt.yaml @@ -1,6 +1,5 @@ model__input_dim: null model__hidden_dims: [256, 256, 256, 256] -model__layer_dim: 2 model__output_dim: 1 model__activation: 'ReLU' criterion: 'BCEWithLogitsLoss' diff --git a/cyclops/models/neural_nets/mlp.py b/cyclops/models/neural_nets/mlp.py index 25eb550e8..f48a27460 100644 --- a/cyclops/models/neural_nets/mlp.py +++ b/cyclops/models/neural_nets/mlp.py @@ -50,13 +50,17 @@ def __init__( self.input_dim = input_dim self.hidden_dims = hidden_dims self.output_dim = output_dim - self.activation = get_module("activation", activation) + self.activation = ( + get_module("activation", activation)() + if isinstance(activation, str) + else activation + ) - layers = [self._layer(input_dim, hidden_dims[0], self.activation)] + layers = self._layer(input_dim, hidden_dims[0], self.activation) for i in range(len(hidden_dims) - 1): layers.extend( self._layer( - self.hidden_dims[i] if i > 0 else input_dim, + self.hidden_dims[i], self.hidden_dims[i + 1], activation, ), diff --git a/cyclops/monitor/detector.py b/cyclops/monitor/detector.py index 3ea8c5574..ba03b0d98 100644 --- a/cyclops/monitor/detector.py +++ b/cyclops/monitor/detector.py @@ -7,6 +7,7 @@ from datasets import concatenate_datasets from datasets.arrow_dataset import Dataset +from cyclops.data.slicer import SliceSpec from cyclops.monitor.reductor import Reductor from cyclops.monitor.tester import DCTester, TSTester from cyclops.monitor.utils import get_args @@ -192,6 +193,133 @@ def _detect_shift_sample(self, ds_target: Dataset) -> Dict[str, Any]: "shift_detected": shift_detected, } + def detect_shift_by_subgroup( + self, + ds_target: Dataset, + slice_spec: SliceSpec, + correction: str = "bonferroni", + min_sample_size: int = 30, + batched: bool = True, + batch_size: int = 1000, + num_proc: int = 1, + ) -> Dict[str, Dict[str, Any]]: + """Detect distribution shift independently within each subgroup. + + A model can look stable when tested against the whole target + population while drifting badly for a specific clinically or + socially relevant subgroup (e.g. an age band, sex, or hospital + site) - an aggregate test can mask this. This method runs the + already-fit tester separately on each subgroup of `ds_target` + defined by `slice_spec`, so that subgroup-level shift can be + detected and reported on its own, which is useful for + health-equity-aware monitoring of deployed models. + + Parameters + ---------- + ds_target : Dataset + Target dataset to test for shift, split into subgroups. + slice_spec : SliceSpec + Specification of the subgroups (slices) of `ds_target` to test + independently. See :class:`cyclops.data.slicer.SliceSpec`. + correction : str, optional + Multiple-testing correction applied to the p-value threshold + across all subgroups tested, to control the false-positive + rate that testing many subgroups simultaneously would + otherwise inflate. One of "bonferroni" or "none". Default is + "bonferroni". + min_sample_size : int, optional + Minimum number of samples required in a subgroup for the + shift test to be run. Subgroups with fewer samples than this + are still returned (with their sample size), but with + `p_val`/`distance`/`shift_detected` set to None, since a + statistical test on too few samples is unreliable. Default + is 30. + batched : bool, optional + Whether to filter the dataset in batches. Default is True. + batch_size : int, optional + Batch size to use when filtering. Default is 1000. + num_proc : int, optional + Number of processes to use when filtering. Default is 1. + + Returns + ------- + dict + Dictionary mapping each subgroup's slice name to a dictionary + with keys `p_val`, `distance`, `shift_detected`, and + `sample_size`. + + Examples + -------- + >>> import numpy as np + >>> from datasets import Dataset + >>> from cyclops.data.slicer import SliceSpec + >>> from cyclops.monitor.detector import Detector + >>> from cyclops.monitor.reductor import Reductor + >>> from cyclops.monitor.tester import TSTester + >>> np.random.seed(0) + >>> ds_source = Dataset.from_dict( + ... { + ... "feature_0": np.random.rand(100), + ... "sex": ["M", "F"] * 50, + ... }, + ... ) + >>> ds_target = Dataset.from_dict( + ... { + ... "feature_0": np.random.rand(100), + ... "sex": ["M", "F"] * 50, + ... }, + ... ) + >>> reductor = Reductor("nored", feature_columns=["feature_0"]) + >>> tester = TSTester("mmd") + >>> detector = Detector("sensitivity_test", reductor, tester) + >>> detector.fit(ds_source) + >>> slice_spec = SliceSpec( + ... spec_list=[{"sex": {"value": "M"}}, {"sex": {"value": "F"}}], + ... ) + >>> results = detector.detect_shift_by_subgroup(ds_target, slice_spec) + + """ + if correction not in ("bonferroni", "none"): + raise ValueError( + f"Unknown correction method: {correction}. " + "Must be one of 'bonferroni', 'none'.", + ) + + slices = slice_spec.get_slices() + base_threshold = self.tester.p_val_threshold + threshold = ( + base_threshold / len(slices) + if correction == "bonferroni" + else base_threshold + ) + + results: Dict[str, Dict[str, Any]] = {} + for slice_name, slice_fn in slices.items(): + ds_subgroup = ds_target.filter( + slice_fn, + batched=batched, + batch_size=batch_size, + num_proc=num_proc, + ) + sample_size = ds_subgroup.shape[0] + if sample_size < min_sample_size: + results[slice_name] = { + "p_val": None, + "distance": None, + "shift_detected": None, + "sample_size": sample_size, + } + continue + + drift_result = self._detect_shift_sample(ds_subgroup) + results[slice_name] = { + "p_val": drift_result["p_val"], + "distance": drift_result["distance"], + "shift_detected": 1 if drift_result["p_val"] < threshold else 0, + "sample_size": sample_size, + } + return results + def sensitivity_test( self, ds_source: Dataset, diff --git a/cyclops/monitor/explainer.py b/cyclops/monitor/explainer.py index 1ea93bf1a..057522d9b 100644 --- a/cyclops/monitor/explainer.py +++ b/cyclops/monitor/explainer.py @@ -14,10 +14,20 @@ if TYPE_CHECKING: import shap else: - shap = import_optional_module( - "shap", - error="warn", - ) + # imported lazily (see _ensure_shap_imported) rather than at module load, + # since shap depends on a third-party package also named `slicer`, which + # can collide with this repo's own cyclops/data/slicer.py under some + # import mechanisms (e.g. doctest's per-file `sys.path` handling) if shap + # were imported merely by importing this module. + shap = None + + +def _ensure_shap_imported() -> Any: + """Import shap on first use and cache it at module scope.""" + global shap # noqa: PLW0603 + if shap is None: + shap = import_optional_module("shap", error="warn") + return shap class Explainer: @@ -38,6 +48,7 @@ def __init__( data: Optional[Any] = None, explainer_type: Optional[str] = None, ) -> None: + _ensure_shap_imported() self.model = model self.data = data self.explainer_type = explainer_type @@ -51,6 +62,8 @@ def get_explainer(self) -> Any: explainer = shap.DeepExplainer(self.model, self.data) elif self.explainer_type == "gradient": explainer = shap.GradientExplainer(self.model, self.data) + elif self.data is not None: + explainer = shap.Explainer(self.model, self.data) else: explainer = shap.Explainer(self.model) return explainer diff --git a/cyclops/monitor/plotter.py b/cyclops/monitor/plotter.py index 79fdbea17..cd0d6a6d5 100644 --- a/cyclops/monitor/plotter.py +++ b/cyclops/monitor/plotter.py @@ -78,7 +78,7 @@ def errorfill( """Create custom error fill.""" ax = ax if ax is not None else plt.gca() if color is None: - color = next(ax._get_lines.prop_cycler)["color"] + color = ax._get_lines.get_next_color() if np.isscalar(yerr) or len(yerr) == len(y): ymin = y - yerr ymax = y + yerr @@ -187,120 +187,6 @@ def set_bars_color(bars: mpl.container.BarContainer, color: str) -> None: bar_item.set_color(color) -def plot_label_distribution( - X: pd.DataFrame, - y: pd.DataFrame, - label: str, - features: List[str], -) -> None: - """Set color attribute for bars in bar plots. - - Parameters - ---------- - bars: mpl.container.BarContainer - Bars. - X: pd.DataFrame - Feature values. - y: pd.DataFrame - Label outcome values. - label: str - Column name of outcome variable. - features: list - Names of features to plot. - - """ - data = pd.concat([X, y], axis=1) - data_pos = data.loc[data[label] == 1] - data_neg = data.loc[data[label] == 0] - _, axs = plt.subplots(2, 2, figsize=(30, 15), tight_layout=True) - - # Across age. - age = None - ages = data[age] - ages_pos = data_pos[age] - ages_neg = data_neg[age] - print( - f"Mean Age: Outcome present: {np.array(ages_pos).mean()}, \ - No outcome: {np.array(ages_neg).mean()}", - ) - - (_, bins, _) = axs[0][0].hist(ages, bins=50, alpha=0.5, color="g") - axs[0][0].hist(ages_pos, bins=bins, alpha=0.5, color="r") - setup_plot( - axs[0][0], - "Age distribution", - "Age", - "Num. of encounters", - ["All", "Outcome present"], - ) - - # Across sex. - sex = None - sex = list(data[sex].unique()) - sex_counts = list(data[sex].value_counts()) - sex_counts_pos = list(data_pos[sex].value_counts()) - - sex_bars = axs[0][1].bar(sex, sex_counts, alpha=0.5) - set_bars_color(sex_bars, "g") - sex_bars_pos = axs[0][1].bar(sex, sex_counts_pos, alpha=0.5) - set_bars_color(sex_bars_pos, "r") - setup_plot( - axs[0][1], - "Sex distribution", - "Sex", - "Num. of encounters", - ["All", "Outcome present"], - ) - - # Across features. - len_features = len(features) - width = 0.04 - x = np.arange(0, len([0, 1])) - - for i, feature in enumerate(features): - feature_counts = list(data[feature].value_counts()) - feature_counts_pos = list(data_pos[feature].value_counts()) - if len(feature_counts) == 1: - feature_counts.append(0) - icd_counts_pos: List[int] = [] - if len(icd_counts_pos) == 1: - feature_counts_pos.append(0) - position = x + (width * (1 - len_features) / 2) + i * width - feature_bars = axs[1][0].bar(position, feature_counts, width=width, alpha=0.5) - set_bars_color(feature_bars, "g") - feature_bars_pos = axs[1][0].bar( - position, - feature_counts_pos, - width=width, - alpha=0.5, - ) - set_bars_color(feature_bars_pos, "r") - - setup_plot( - axs[1][0], - "Feature distribution", - "Feature", - "Num. of encounters", - ["All", "Outcome present"], - ) - - # Across labels. - label_counts = y.value_counts().to_dict().values() - labels = data[label].value_counts().to_dict().keys() - - label_bars = axs[1][1].bar(labels, label_counts, alpha=0.5) - set_bars_color(label_bars, "g") - setup_plot( - axs[1][1], - "Outcome distribution", - "Outcome", - "Num. of encounters", - ["All"], - ) - - plt.show() - - def plot_drift_experiment( results: dict[str, dict[str, np.ndarray[float, np.dtype[np.float64]]]], plot_distance=False, diff --git a/cyclops/monitor/reductor.py b/cyclops/monitor/reductor.py index d306da2b1..0b2bce9dd 100644 --- a/cyclops/monitor/reductor.py +++ b/cyclops/monitor/reductor.py @@ -71,7 +71,7 @@ def __init__( dr_method: str, batch_size: int = 32, num_workers: int = 0, - device: str = None, + device: Optional[str] = None, transforms: Optional[Union[Callable, Compose]] = None, feature_columns: Optional[Union[str, List[str]]] = None, **kwargs: Any, @@ -80,7 +80,7 @@ def __init__( self.batch_size = batch_size self.num_workers = num_workers self.device = device - if isinstance(transforms, Compose): + if Compose is not None and isinstance(transforms, Compose): self.transforms = partial(apply_transforms, transforms=transforms) else: self.transforms = transforms @@ -120,7 +120,7 @@ def __init__( else: self.model = wrap_model(self.model) - def load_model(self, output_path: str = None) -> None: + def load_model(self, output_path: Optional[str] = None) -> None: """Load pre-trained model from path. For scikit-learn models, a pickle is loaded from disk. For the torch models, the diff --git a/cyclops/monitor/tester.py b/cyclops/monitor/tester.py index ac4b1d7c5..bc2e6cf47 100644 --- a/cyclops/monitor/tester.py +++ b/cyclops/monitor/tester.py @@ -165,6 +165,7 @@ def __init__( self.tester_method = tester_method self.method: Any = None self.p_val_threshold = p_val_threshold + self._base_p_val_threshold = p_val_threshold # dict where the key is the string of each test_method # and the value is the class of the test_method @@ -256,6 +257,7 @@ def test_shift( Tuple[float, float] p-value and distance between reference and target datasets """ + num_features = None if isinstance(X_t, np.ndarray): X_t = X_t.astype("float32") num_features = X_t.shape[1] @@ -287,8 +289,10 @@ def test_shift( p_val = p_val[idx] dist = dist[idx] - if self.tester_method in ["ks", "chi2", "fet", "tabular"]: - self.p_val_threshold = self.p_val_threshold / num_features + if self.tester_method in ["ks", "chi2", "fet", "tabular"] and num_features: + # Bonferroni-correct relative to the original threshold each call, + # so repeated calls (e.g. in Detector's loops) don't compound. + self.p_val_threshold = self._base_p_val_threshold / num_features return p_val, dist @@ -408,6 +412,7 @@ def __init__( self.p_val_threshold = p_val_threshold self.method_args = kwargs self.tester: Any = None + self.X_s: Any = None self.tester_methods = { "spot_the_diff": SpotTheDiffDrift, @@ -431,6 +436,7 @@ def fit( """Initialize test method to source data.""" if isinstance(X_s, np.ndarray): X_s = X_s.astype("float32") + self.X_s = X_s if self.tester_method == "spot_the_diff": if not isinstance(X_s, np.ndarray): @@ -475,6 +481,109 @@ def test_shift( dist = preds["data"]["distance"] return p_val, dist + def explain_shift( + self, + X_t: np.ndarray[float, np.dtype[np.float64]], + feature_names: Optional[List[str]] = None, + **explainer_kwargs: Any, + ) -> Dict[str, float]: + """Explain which features are driving a detected shift. + + Only supported for ``tester_method="classifier"``: that test trains + a classifier to discriminate reference (source) from test (target) + samples, which SHAP can then explain directly - the features that + most strongly indicate a sample belongs to the target distribution + are the ones most responsible for the detected drift. Must be + called after :meth:`fit` and :meth:`test_shift`. + + Parameters + ---------- + X_t : np.ndarray + Target data to compute SHAP values for (the same data, or data + from the same distribution, passed to :meth:`test_shift`). + feature_names : list of str, optional + Names for each feature/column of `X_t`, used as keys in the + returned dictionary. Defaults to stringified column indices. + **explainer_kwargs : Any + Additional keyword arguments passed to + :class:`cyclops.monitor.explainer.Explainer`. + + Returns + ------- + dict + Dictionary mapping each feature name to its mean absolute SHAP + value, sorted by descending importance (most drift-responsible + feature first). + + Examples + -------- + >>> from cyclops.monitor.tester import DCTester + >>> import numpy as np + >>> np.random.seed(0) + >>> X_s = np.random.normal(0, 1, (100, 10)) + >>> X_t = np.random.normal(1, 1, (100, 10)) + >>> from sklearn.linear_model import LogisticRegression + >>> model = LogisticRegression() + >>> tester = DCTester("classifier", model=model) + >>> tester.fit(X_s) + >>> p_val, dist = tester.test_shift(X_t) + >>> importances = tester.explain_shift(X_t) # doctest: +SKIP + + """ + # imported lazily: cyclops.monitor.explainer eagerly imports shap, and + # importing shap at module load time (i.e. every time cyclops.monitor + # is imported) is both unnecessary for users who never call + # explain_shift() and can collide with this repo's own + # cyclops/data/slicer.py under some import mechanisms (e.g. doctest's + # per-file `sys.path` handling), since shap depends on a third-party + # package also named `slicer`. + from cyclops.monitor.explainer import Explainer # noqa: PLC0415 + + if self.tester_method != "classifier": + raise ValueError( + 'explain_shift() is only supported for tester_method="classifier" ' + f"(got {self.tester_method!r}); the other domain-classifier " + "methods don't expose a single trained model to explain.", + ) + if self.tester is None: + raise ValueError("Must call fit() and test_shift() before explain_shift().") + + try: + trained_model = self.tester._detector.model # noqa: SLF001 + except AttributeError as exc: + raise RuntimeError( + "Could not access the trained classifier from the underlying " + "alibi-detect ClassifierDrift detector; explain_shift() may be " + "incompatible with the installed alibi-detect version.", + ) from exc + predict_fn = getattr( + trained_model, + "predict_proba", + getattr(trained_model, "predict", trained_model), + ) + # cap background data size, since SHAP's model-agnostic explainers scale + # poorly with the number of background samples + background = self.X_s[:100] if self.X_s is not None else None + + if isinstance(X_t, np.ndarray): + X_t = X_t.astype("float32") + explainer = Explainer(predict_fn, data=background, **explainer_kwargs) + shap_values = np.asarray(explainer.get_shap_values(X_t).values) + if shap_values.ndim == 3: # (samples, features, classes/outputs) + shap_values = shap_values.mean(axis=-1) + importances = np.abs(shap_values).mean(axis=0) + + if feature_names is None: + feature_names = [str(i) for i in range(len(importances))] + + return dict( + sorted( + zip(feature_names, importances.tolist()), + key=lambda item: item[1], + reverse=True, + ), + ) + class ContextMMDWrapper: """Wrapper for ContextMMDDrift.""" @@ -488,6 +597,7 @@ def __init__( backend: str = "tensorflow", p_val: float = 0.05, preprocess_x_ref: bool = False, + preprocess_at_init: bool = True, update_ref: Optional[Dict[str, int]] = None, preprocess_fn: Optional[Callable[..., Any]] = None, x_kernel: Optional[Callable[..., Any]] = None, @@ -505,25 +615,26 @@ def __init__( c_source = context_generator.transform(ds_source) - args = [ - backend, - p_val, - preprocess_x_ref, - update_ref, - preprocess_fn, - x_kernel, - c_kernel, - n_permutations, - prop_c_held, - n_folds, - batch_size, - device, - input_shape, - data_type, - verbose, - ] - - self.tester = ContextMMDDrift(X_s, c_source, *args) + self.tester = ContextMMDDrift( + X_s, + c_source, + backend=backend, + p_val=p_val, + x_ref_preprocessed=preprocess_x_ref, + preprocess_at_init=preprocess_at_init, + update_ref=update_ref, + preprocess_fn=preprocess_fn, + x_kernel=x_kernel, + c_kernel=c_kernel, + n_permutations=n_permutations, + prop_c_held=prop_c_held, + n_folds=n_folds, + batch_size=batch_size, + device=device, + input_shape=input_shape, + data_type=data_type, + verbose=verbose, + ) def predict( self, @@ -584,35 +695,36 @@ def __init__( kernel_b = GaussianRBF(trainable=True) if kernel_b is None else kernel_b kernel = DeepKernel(self.proj, kernel_a, kernel_b, eps) - args = [ - backend, - p_val, - x_ref_preprocessed, - preprocess_at_init, - update_x_ref, - preprocess_fn, - n_permutations, - batch_size_permutations, - var_reg, - reg_loss_fn, - train_size, - retrain_from_scratch, - optimizer, - learning_rate, - batch_size, - batch_size_predict, - preprocess_batch_fn, - epochs, - num_workers, - verbose, - train_kwargs, - device, - dataset, - dataloader, - input_shape, - data_type, - ] - self.tester = LearnedKernelDrift(X_s, kernel, *args) + self.tester = LearnedKernelDrift( + X_s, + kernel, + backend=backend, + p_val=p_val, + x_ref_preprocessed=x_ref_preprocessed, + preprocess_at_init=preprocess_at_init, + update_x_ref=update_x_ref, + preprocess_fn=preprocess_fn, + n_permutations=n_permutations, + batch_size_permutations=batch_size_permutations, + var_reg=var_reg, + reg_loss_fn=reg_loss_fn, + train_size=train_size, + retrain_from_scratch=retrain_from_scratch, + optimizer=optimizer, + learning_rate=learning_rate, + batch_size=batch_size, + batch_size_predict=batch_size_predict, + preprocess_batch_fn=preprocess_batch_fn, + epochs=epochs, + num_workers=num_workers, + verbose=verbose, + train_kwargs=train_kwargs, + device=device, + dataset=dataset, + dataloader=dataloader, + input_shape=input_shape, + data_type=data_type, + ) def predict( self, diff --git a/cyclops/monitor/utils.py b/cyclops/monitor/utils.py index d990c5096..63c0c355a 100644 --- a/cyclops/monitor/utils.py +++ b/cyclops/monitor/utils.py @@ -1,244 +1,17 @@ """Utilities for the drift detector module.""" -import datetime -import importlib import inspect -import pickle -from datetime import timedelta -from itertools import cycle -from shutil import get_terminal_size -from threading import Thread -from time import sleep -from typing import TYPE_CHECKING, Any, Dict, Generator, List, Optional, Tuple +from typing import TYPE_CHECKING, Any, Dict, Optional -import numpy as np -import pandas as pd -from sklearn import metrics -from sklearn.preprocessing import StandardScaler - -from cyclops.models.neural_nets.gru import GRUModel -from cyclops.models.neural_nets.lstm import LSTMModel -from cyclops.models.neural_nets.rnn import RNNModel -from cyclops.models.wrappers import SKModel from cyclops.utils.optional import import_optional_module if TYPE_CHECKING: import torch from torch import nn - from torch.optim import Optimizer - from torch.utils.data import DataLoader, TensorDataset - from torch.utils.data import Dataset as TorchDataset else: torch = import_optional_module("torch", error="warn") nn = import_optional_module("torch.nn", error="warn") - Optimizer = import_optional_module( - "torch.optim", - attribute="Optimizer", - error="warn", - ) - DataLoader = import_optional_module( - "torch.utils.data", - attribute="DataLoader", - error="warn", - ) - TensorDataset = import_optional_module( - "torch.utils.data", - attribute="TensorDataset", - error="warn", - ) - TorchDataset = import_optional_module( - "torch.utils.data", - attribute="Dataset", - error="warn", - ) - - -def print_metrics_binary( - y_test_labels: Any, - y_pred_values: Any, - y_pred_labels: Any, - verbose: int = 1, -) -> Dict[str, Any]: - """Print metrics for binary classification.""" - conf_matrix = metrics.confusion_matrix(y_test_labels, y_pred_labels) - if verbose: - print("confusion matrix:") - print(conf_matrix) - conf_matrix = conf_matrix.astype(np.float32) - tn, fp, fn, tp = conf_matrix.ravel() - acc = (tn + tp) / np.sum(conf_matrix) - prec0 = tn / (tn + fn) - prec1 = tp / (tp + fp) - rec0 = tn / (tn + fp) - rec1 = tp / (tp + fn) - - auroc = metrics.roc_auc_score(y_test_labels, y_pred_values) - - (precisions, recalls, _) = metrics.precision_recall_curve( - y_test_labels, - y_pred_values, - ) - auprc = metrics.auc(recalls, precisions) - minpse = np.max([min(x, y) for (x, y) in zip(precisions, recalls)]) - - if verbose: - print(f"accuracy: {acc}") - print(f"precision class 0: {prec0}") - print(f"precision class 1: {prec1}") - print(f"recall class 0: {rec0}") - print(f"recall class 1: {rec1}") - print(f"AUC of ROC: {auroc}") - print(f"AUC of PRC: {auprc}") - print(f"min(+P, Se): {minpse}") - - return { - "acc": acc, - "prec0": prec0, - "prec1": prec1, - "rec0": rec0, - "rec1": rec1, - "auroc": auroc, - "auprc": auprc, - "minpse": minpse, - } - - -def load_ckp( - checkpoint_fpath: str, - model: nn.Module, -) -> Tuple[nn.Module, Optimizer, int]: - """Load checkpoint.""" - checkpoint = torch.load(checkpoint_fpath) # type: ignore - model.load_state_dict(checkpoint["model"]) - optimizer = checkpoint["optimizer"] - return model, optimizer, checkpoint["n_epochs"] - - -def get_device() -> torch.device: - """Get device.""" - return torch.device("cuda") if torch.cuda.is_available() else torch.device("cpu") - - -def get_temporal_model(model: str, model_params: Dict[str, Any]) -> nn.Module: - """Get temporal model. - - Parameters - ---------- - model: string - String with model name (e.g. rnn, lstm, gru). - - """ - models = {"rnn": RNNModel, "lstm": LSTMModel, "gru": GRUModel} - return models[model.lower()](**model_params) - - -class Data(TorchDataset[Tuple[torch.Tensor, torch.Tensor]]): - """Data class.""" - - def __init__(self, inputs: pd.DataFrame, target: pd.DataFrame) -> None: - """Initialize Data class.""" - self.inputs = inputs - self.target = target - - def __getitem__(self, idx: int) -> Tuple[Any, Any]: - """Get item for iterator. - - Parameters - ---------- - idx: int - Index of sample to fetch from dataset. - - Returns - ------- - tuple - Input and target. - - """ - return self.inputs[idx], self.target[idx] - - def __len__(self) -> int: - """Return size of dataset, i.e. no. of samples. - - Returns - ------- - int - Size of dataset. - - """ - return len(self.target) - - def dim(self) -> Any: - """Get dataset dimensions (no. of features). - - Returns - ------- - int - Number of features. - - """ - return self.inputs.size(dim=1) - - def to_loader( - self, - batch_size: int, - num_workers: int = 0, - shuffle: bool = False, - pin_memory: bool = True, - ) -> DataLoader[Any]: - """Create dataloader. - - Returns - ------- - DataLoader with input data - - """ - return DataLoader( - TensorDataset(self.inputs, self.target), - batch_size=batch_size, - num_workers=num_workers, - shuffle=shuffle, - pin_memory=pin_memory, - ) - - -def get_data(X: np.ndarray[float, np.dtype[np.float64]], y: List[int]) -> Data: - """Convert pandas dataframe to dataset. - - Parameters - ---------- - X: numpy matrix - Data containing features in the form of [samples, timesteps, features]. - y: list - List of labels. - - """ - inputs = torch.tensor(X, dtype=torch.float32) - target = torch.tensor(y, dtype=torch.float32) - return Data(inputs, target) - - -def run_model( - model_name: str, - X: pd.DataFrame, - y: pd.DataFrame, - X_val: pd.DataFrame, - y_val: pd.DataFrame, -) -> SKModel: - """Choose and run a model on the data and return the best model.""" - if model_name == "mlp": - model = SKModel("mlp", save_path="./mlp.pkl") - model.fit(X, y, X_val, y_val) - elif model_name == "lr": - model = SKModel("lr", save_path="./lr.pkl") - model.fit(X, y, X_val, y_val) - elif model_name == "rf": - model = SKModel("rf", save_path="./rf.pkl") - model.fit(X, y, X_val, y_val) - elif model_name == "xgb": - model = SKModel("xgb", save_path="./xgb.pkl") - model.fit(X, y, X_val, y_val) - return model def get_args(obj: Any, kwargs: Dict[str, Any]) -> Dict[str, Any]: @@ -267,242 +40,6 @@ def get_args(obj: Any, kwargs: Dict[str, Any]) -> Dict[str, Any]: return args -def get_obj_from_str(string: str, reload: bool = False) -> Any: - """Get object from string.""" - module, cls = string.rsplit(".", 1) - if reload: - module_imp = importlib.import_module(module) - importlib.reload(module_imp) - return getattr(importlib.import_module(module, package=None), cls) - - -def load_model(model_path: str) -> Any: - """Load pre-trained model from path. - - Loads pre-trained model from specified model path. - For scikit-learn models, a pickle is loaded from disk. - For the pytorch models, the "state_dict" is loaded from disk. - - Returns - ------- - model - loaded pre-trained model - - """ - file_type = model_path.split(".")[-1] - if file_type in ("pkl", "pickle"): - with open(model_path, "rb") as file: - model = pickle.load(file) - elif file_type == "pt": - model = torch.load(model_path) # type: ignore - return model - - -def save_model(model: Any, output_path: str) -> None: - """Save the model to disk. - - For scikit-learn models, a pickle is saved to disk. - For the pytorch models, the "state_dict" is saved to disk. - - Parameters - ---------- - output_path: String - path to save the model to - - """ - file_type = output_path.split(".")[-1] - if file_type in ("pkl", "pickle"): - with open(output_path, "wb") as file: - pickle.dump(model, file) - elif file_type == "pt": - torch.save(model.state_dict(), output_path) - - -def scale(x: pd.DataFrame) -> pd.DataFrame: - """Scale columns of temporal dataframe. - - Returns - ------- - model: torch.nn.Module - feed forward neural network model. - - """ - numerical_cols = [ - col for col in x if not np.isin(x[col].dropna().unique(), [0, 1]).all() - ] - - for col in numerical_cols: - scaler = StandardScaler().fit(x[col].values.reshape(-1, 1)) - x[col] = pd.Series( - np.squeeze(scaler.transform(x[col].values.reshape(-1, 1))), - index=x[col].index, - ) - - return x - - -def daterange( - start_date: datetime.date, - end_date: datetime.date, - stride: int, - window: int, -) -> Generator[datetime.date, None, None]: - """Output a range of dates. - - Outputs a range of dates after applying a shift of - a given stride and window adjustment. - - Returns - ------- - datetime.date - range of dates after stride and window adjustment. - - """ - for date in range(int((end_date - start_date).days)): - if start_date + timedelta(date * stride + window) < end_date: - yield start_date + timedelta(date * stride) - - -def get_serving_data( - X: pd.DataFrame, - y: pd.DataFrame, - admin_data: pd.DataFrame, - start_date: datetime.date, - end_date: datetime.date, - stride: int = 1, - window: int = 1, - ids_to_exclude: Optional[List[str]] = None, - encounter_id: str = "encounter_id", - admit_timestamp: str = "admit_timestamp", -) -> Dict[str, Any]: - """Transform a static set of patient encounters with timestamps into serving data. - - Transforms a static set of patient encounters with timestamps into - serving data that ranges from a given start date and goes until - a given end date with a constant window and stride length. - - Returns - ------- - dictionary - dictionary containing keys timestamp, X and y - - """ - X_target_stream = [] - y_target_stream = [] - timestamps = [] - - admit_df = admin_data[[encounter_id, admit_timestamp]].sort_values( - by=admit_timestamp, - ) - for single_date in daterange(start_date, end_date, stride, window): - if single_date.month == 1 and single_date.day == 1: - print( - single_date.strftime("%Y-%m-%d"), - "-", - (single_date + timedelta(days=window)).strftime("%Y-%m-%d"), - ) - encounters_inwindow = admit_df.loc[ - ( - (single_date + timedelta(days=window)).strftime("%Y-%m-%d") - > admit_df[admit_timestamp].dt.strftime("%Y-%m-%d") - ) - & ( - admit_df[admit_timestamp].dt.strftime("%Y-%m-%d") - >= single_date.strftime("%Y-%m-%d") - ), - encounter_id, - ].unique() - if ids_to_exclude is not None: - encounters_inwindow = [ - x for x in encounters_inwindow if x not in ids_to_exclude - ] - encounter_ids = X.index.get_level_values(0).unique() - X_inwindow = X.loc[X.index.get_level_values(0).isin(encounters_inwindow)] - y_inwindow = pd.DataFrame(y[np.in1d(encounter_ids, encounters_inwindow)]) - if not X_inwindow.empty: - X_target_stream.append(X_inwindow) - y_target_stream.append(y_inwindow) - timestamps.append( - (single_date + timedelta(days=window)).strftime("%Y-%m-%d"), - ) - return {"timestamps": timestamps, "X": X_target_stream, "y": y_target_stream} - - -def reshape_2d_to_3d(data: pd.DataFrame, num_timesteps: int) -> pd.DataFrame: - """Reshape 2D data to 3D data.""" - data = data.unstack() - num_encounters = data.shape[0] - return data.values.reshape((num_encounters, num_timesteps, -1)) - - -# from https://stackoverflow.com/a/66558182 -class Loader: - """Loaing animation.""" - - def __init__( - self, - desc: str = "Loading...", - end: str = "Done!", - timeout: float = 0.1, - ) -> None: - """Loader-like context manager. - - Parameters - ---------- - desc (str, optional): The loader's description. Defaults to "Loading...". - end (str, optional): Final print. Defaults to "Done!". - timeout (float, optional): Sleep time between prints. Defaults to 0.1. - - """ - self.desc = desc - self.end = end - self.timeout = timeout - - self._thread = Thread(target=self._animate, daemon=True) - self.steps = ["⢿", "⣻", "⣽", "⣾", "⣷", "⣯", "⣟", "⡿"] - self.done = False - - def start(self) -> "Loader": - """Start the loader.""" - self._thread.start() - return self - - def _animate(self) -> None: - """Animate the loader.""" - for cycle_itr in cycle(self.steps): - if self.done: - break - print(f"\r{self.desc} {cycle_itr}", flush=True, end="") - sleep(self.timeout) - - def __enter__(self) -> None: - """Start the thread.""" - self.start() - - def stop(self) -> None: - """Stop the loader.""" - self.done = True - cols = get_terminal_size((80, 20)).columns - print("\r" + " " * cols, end="", flush=True) - print(f"\r{self.end}", flush=True) - - def __exit__(self, exc_type: Any, exc_value: Any, exc_traceback: Any) -> None: - """Stop the thread.""" - # handle exceptions with those variables ^ - self.stop() - - -if __name__ == "__main__": - with Loader("Loading with context manager..."): - for _i in range(10): - sleep(0.25) - - loader = Loader("Loading with object...", "That was fast!", 0.05).start() - for _i in range(10): - sleep(0.25) - loader.stop() - - class DCELoss(torch.nn.Module): """Disagreement Cross Entropy Loss.""" diff --git a/cyclops/report/report.py b/cyclops/report/report.py index df5f9e6cc..81a1c5f43 100644 --- a/cyclops/report/report.py +++ b/cyclops/report/report.py @@ -970,8 +970,9 @@ def log_performance_metrics( results: Dict[str, Any], metric_descriptions: Dict[str, str], pass_fail_thresholds: Union[float, Dict[str, float]] = 0.7, - pass_fail_threshold_fn: Callable[[float, float], bool] = lambda x, - threshold: bool(x >= threshold), + pass_fail_threshold_fn: Callable[[float, float], bool] = lambda x, threshold: ( + bool(x >= threshold) + ), ) -> None: """ Log all performance metrics to the model card report. @@ -1135,16 +1136,23 @@ def export( today_now = synthetic_timestamp else: today_now = dt_datetime.now().strftime("%Y-%m-%d %H:%M:%S") + # filesystem-safe timestamp for the default output filename, so that + # repeated export() calls into the same output_dir don't silently + # overwrite one another and lose trend/history data. + filename_timestamp = today_now.replace(" ", "_").replace(":", "-") current_report_metrics: Union[ List[List[PerformanceMetric]], List[PerformanceMetric] ] = [] sweep_metrics(self._model_card, current_report_metrics) - current_report_metrics_set = ( - current_report_metrics[0] - if isinstance(current_report_metrics[0], list) - else [current_report_metrics[0]] - ) + if len(current_report_metrics) == 0: + current_report_metrics_set: List[PerformanceMetric] = [] + else: + current_report_metrics_set = ( + current_report_metrics[0] + if isinstance(current_report_metrics[0], list) + else [current_report_metrics[0]] + ) report_paths = glob.glob( os.path.join( @@ -1160,7 +1168,11 @@ def export( latest_report = ModelCard.model_validate_json(f_handle.read()) latest_report_metric_cards: List[List[MetricCard]] = [] sweep_metric_cards(latest_report, latest_report_metric_cards) - latest_report_metric_cards_set = latest_report_metric_cards[0] + latest_report_metric_cards_set = ( + latest_report_metric_cards[0] + if len(latest_report_metric_cards) != 0 + else None + ) else: latest_report_metric_cards_set = None # check if overview section exists @@ -1206,7 +1218,7 @@ def export( report_path = os.path.join( self.output_dir, "cyclops_report", - output_filename or "model_card.html", + output_filename or f"model_card_{filename_timestamp}.html", ) self._write_file(report_path, content) if save_json: diff --git a/cyclops/report/templates/model_report/macros.jinja b/cyclops/report/templates/model_report/macros.jinja index f75ab204b..bf9befd79 100644 --- a/cyclops/report/templates/model_report/macros.jinja +++ b/cyclops/report/templates/model_report/macros.jinja @@ -15,7 +15,7 @@
  • {# {% for name, value in values %} #} {# {% if value %} #} - {{ values.content | safe }} + {{ values.content }} {# {% endif %} #} {# {% endfor %} #}
    • diff --git a/cyclops/report/utils.py b/cyclops/report/utils.py index 7f20e82d9..4f5d2c77f 100644 --- a/cyclops/report/utils.py +++ b/cyclops/report/utils.py @@ -734,6 +734,8 @@ def _process_metric_name( "Multilabel", ): name = metric["type"][10:] + else: + name = metric["type"] for key, value in _METRIC_NAMES_DISPLAY_MAP.items(): name = name.replace(key, value) else: @@ -1049,13 +1051,18 @@ def create_metric_card_plot( return GraphicsCollection(description="plot", collection=[graphic]) -def regex_replace(string: str, find: str, replace: str) -> str: +def regex_replace(string: Any, find: str, replace: str) -> Any: """Replace a regex pattern with a string.""" + if not isinstance(string, str): + # e.g. Jinja's Undefined when a template indexes into an empty list + return string return sub(find, replace, string) -def regex_search(string: str, find: str) -> List[Any]: +def regex_search(string: Any, find: str) -> List[Any]: """Search a regex pattern in a string and return the match.""" + if not isinstance(string, str): + return [] return findall(r"\((.*?)\)", string) diff --git a/cyclops/utils/file.py b/cyclops/utils/file.py index ab312acf3..39e7d76d8 100644 --- a/cyclops/utils/file.py +++ b/cyclops/utils/file.py @@ -55,7 +55,8 @@ def exchange_extension(file_path: str, new_ext: str) -> str: # Remove a leading dot new_ext = new_ext.strip(".") _, old_ext = os.path.splitext(file_path) - return file_path[: -len(old_ext)] + "." + new_ext + stem = file_path[: -len(old_ext)] if old_ext else file_path + return stem + "." + new_ext def process_file_save_path( diff --git a/docs/source/monitoring.rst b/docs/source/monitoring.rst index 591b5bd00..526086549 100644 --- a/docs/source/monitoring.rst +++ b/docs/source/monitoring.rst @@ -1,6 +1,17 @@ Monitoring ========== +.. note:: + + This page covers tracking a model's *logged performance metrics* over time + through repeated report card evaluations. To proactively test whether the + data a deployed model is seeing has statistically drifted from its + training/reference distribution - before a performance drop is even + observed - see the :doc:`drift detection API ` + (:mod:`cyclops.monitor`), which implements two-sample statistical tests, + the Detectron harmful-covariate-shift test, and clinically meaningful + shift simulators (e.g. by age, sex, hospital type, or time). + After initial evaluation and model report generation, how can we monitor model performance over time? diff --git a/docs/source/reference/api/cyclops.monitor.rst b/docs/source/reference/api/cyclops.monitor.rst index d513a61ca..d24fba32f 100644 --- a/docs/source/reference/api/cyclops.monitor.rst +++ b/docs/source/reference/api/cyclops.monitor.rst @@ -14,3 +14,7 @@ cyclops.monitor clinical_applicator synthetic_applicator + detector + reductor + tester + explainer diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 0105f629e..291043ec5 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -5,3 +5,4 @@ Tutorials :maxdepth: 3 tutorials_use_cases + tutorials_monitor diff --git a/pyproject.toml b/pyproject.toml index a4ec9e067..a77e47c55 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pycyclops" -version = "0.2.12" +version = "0.3.0" description = "Framework for healthcare ML implementation" authors = [{ name = "Vector AI Engineering", email = "cyclops@vectorinstitute.ai" }] license = { text = "Apache-2.0" } diff --git a/tests/cyclops/data/test_slicer.py b/tests/cyclops/data/test_slicer.py index 225f7e491..22ec08a99 100644 --- a/tests/cyclops/data/test_slicer.py +++ b/tests/cyclops/data/test_slicer.py @@ -6,6 +6,7 @@ import cycquery.ops as qo import numpy as np import pandas as pd +import pyarrow as pa import pytest from cycquery import OMOPQuerier from datasets import Dataset @@ -297,6 +298,27 @@ def test_filter_datetime( assert result.all() +def test_filter_datetime_day(): + """Test that filter_datetime's `day` argument filters on day, not year. + + Regression test: `day` filtering used `pc.year(...)` instead of + `pc.day(...)`, so it silently matched on year instead of day of + month. Self-contained (no database), unlike the other filter_datetime + tests in this file, which all require a live Synthea database and + are therefore excluded from CI. + """ + dates = pd.to_datetime( + ["2020-01-05", "2020-02-14", "2020-03-14", "2020-04-21"], + ) + table = pa.table({"visit_date": dates}) + + result = filter_datetime(table, column_name="visit_date", day=14) + assert result == [False, True, True, False] + + result = filter_datetime(table, column_name="visit_date", day=[5, 21]) + assert result == [True, False, False, True] + + @pytest.mark.integration_test() @pytest.mark.parametrize( ("column_name", "contains", "negate", "keep_nulls"), diff --git a/tests/cyclops/evaluate/metrics/experimental/test_brier_score.py b/tests/cyclops/evaluate/metrics/experimental/test_brier_score.py new file mode 100644 index 000000000..8ad6dd5d5 --- /dev/null +++ b/tests/cyclops/evaluate/metrics/experimental/test_brier_score.py @@ -0,0 +1,164 @@ +"""Tests for the Brier score metric.""" + +import array_api_compat.torch +import numpy as np +import numpy.array_api as anp +import pytest +import torch +from sklearn.metrics import brier_score_loss + +from cyclops.evaluate.metrics.experimental import ( + BinaryBrierScore, + MulticlassBrierScore, +) +from cyclops.evaluate.metrics.experimental.functional import ( + binary_brier_score, + multiclass_brier_score, +) + + +@pytest.mark.parametrize("xp", [anp, array_api_compat.torch]) +def test_binary_brier_score_matches_sklearn(xp): + """Binary Brier score must match sklearn's brier_score_loss.""" + target_list = [0, 1, 1, 0, 1, 0, 0, 1] + preds_list = [0.1, 0.9, 0.8, 0.3, 0.4, 0.2, 0.6, 0.7] + expected = brier_score_loss(target_list, preds_list) + + target = xp.asarray(target_list) + preds = xp.asarray(preds_list) + result = binary_brier_score(target, preds) + + assert float(result) == pytest.approx(expected, abs=1e-5) + + +def test_binary_brier_score_perfect_predictions(): + """Brier score for perfect predictions must be 0.""" + target = anp.asarray([0, 1, 0, 1]) + preds = anp.asarray([0.0, 1.0, 0.0, 1.0]) + assert float(binary_brier_score(target, preds)) == pytest.approx(0.0) + + +def test_binary_brier_score_worst_predictions(): + """Brier score for maximally wrong predictions must be 1.""" + target = anp.asarray([0, 1, 0, 1]) + preds = anp.asarray([1.0, 0.0, 1.0, 0.0]) + assert float(binary_brier_score(target, preds)) == pytest.approx(1.0) + + +def test_binary_brier_score_from_logits(): + """Logits (values outside [0, 1]) must be converted via sigmoid.""" + target = anp.asarray([0, 1, 1, 0]) + logits_np = np.asarray([-3.0, 3.0, 2.0, -1.0]) + logits = anp.asarray(logits_np) + + result = float(binary_brier_score(target, logits)) + expected = brier_score_loss([0, 1, 1, 0], 1 / (1 + np.exp(-logits_np))) + assert result == pytest.approx(expected, abs=1e-5) + + +def test_binary_brier_score_ignore_index(): + """Values matching ignore_index must be excluded.""" + target = anp.asarray([0, 1, 1, -1]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.99]) + result = binary_brier_score(target, preds, ignore_index=-1) + expected = brier_score_loss([0, 1, 1], [0.1, 0.9, 0.8]) + assert float(result) == pytest.approx(expected, abs=1e-5) + + +def test_binary_brier_score_invalid_target_raises(): + """Non-binary target values must raise.""" + target = anp.asarray([0, 1, 2]) + preds = anp.asarray([0.1, 0.9, 0.8]) + with pytest.raises(RuntimeError): + binary_brier_score(target, preds) + + +class TestBinaryBrierScoreClass: + """Tests for the BinaryBrierScore metric class.""" + + def test_single_call(self): + """Test single-call usage matches the functional API.""" + target = anp.asarray([0, 1, 1, 0]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + metric = BinaryBrierScore() + assert float(metric(target, preds)) == pytest.approx( + float(binary_brier_score(target, preds)), + ) + + def test_streaming_matches_batch(self): + """Accumulating over multiple updates must match a single batch call.""" + target = [0, 1, 1, 0, 1, 0] + preds = [0.1, 0.9, 0.8, 0.3, 0.4, 0.2] + + batch_result = float( + binary_brier_score(anp.asarray(target), anp.asarray(preds)), + ) + + metric = BinaryBrierScore() + for t, p in zip([target[:3], target[3:]], [preds[:3], preds[3:]]): + metric.update(anp.asarray(t), anp.asarray(p)) + streaming_result = float(metric.compute()) + + assert streaming_result == pytest.approx(batch_result, abs=1e-5) + + def test_torch_backend(self): + """Test the metric works with a torch backend.""" + target = torch.tensor([0, 1, 1, 0]) + preds = torch.tensor([0.1, 0.9, 0.8, 0.3]) + metric = BinaryBrierScore() + result = metric(target, preds) + assert isinstance(result, torch.Tensor) + assert float(result) == pytest.approx(0.0375, abs=1e-4) + + +def test_multiclass_brier_score_matches_manual_computation(): + """Multiclass Brier score must match a manually one-hot-encoded MSE.""" + target = anp.asarray([0, 1, 2]) + preds_np = np.asarray([[0.7, 0.2, 0.1], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]]) + preds = anp.asarray(preds_np) + + one_hot = np.asarray([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float) + expected = np.mean(np.sum((preds_np - one_hot) ** 2, axis=1)) + + result = multiclass_brier_score(target, preds, num_classes=3) + assert float(result) == pytest.approx(expected, abs=1e-5) + + +def test_multiclass_brier_score_perfect_predictions(): + """Multiclass Brier score for perfect one-hot predictions must be 0.""" + target = anp.asarray([0, 1, 2]) + preds = anp.asarray([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + assert float(multiclass_brier_score(target, preds, num_classes=3)) == pytest.approx( + 0.0, + ) + + +def test_multiclass_brier_score_invalid_num_classes(): + """num_classes < 2 must raise.""" + target = anp.asarray([0, 1]) + preds = anp.asarray([[1.0], [1.0]]) + with pytest.raises(ValueError, match="num_classes"): + multiclass_brier_score(target, preds, num_classes=1) + + +def test_multiclass_brier_score_wrong_preds_shape(): + """Preds without one more dimension than target must raise.""" + target = anp.asarray([0, 1, 2]) + preds = anp.asarray([0.1, 0.9, 0.8]) + with pytest.raises(ValueError, match="preds"): + multiclass_brier_score(target, preds, num_classes=3) + + +class TestMulticlassBrierScoreClass: + """Tests for the MulticlassBrierScore metric class.""" + + def test_single_call(self): + """Test single-call usage matches the functional API.""" + target = anp.asarray([0, 1, 2]) + preds = anp.asarray( + [[0.7, 0.2, 0.1], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]], + ) + metric = MulticlassBrierScore(num_classes=3) + assert float(metric(target, preds)) == pytest.approx( + float(multiclass_brier_score(target, preds, num_classes=3)), + ) diff --git a/tests/cyclops/evaluate/metrics/experimental/test_calibration_error.py b/tests/cyclops/evaluate/metrics/experimental/test_calibration_error.py new file mode 100644 index 000000000..18ed8ba4d --- /dev/null +++ b/tests/cyclops/evaluate/metrics/experimental/test_calibration_error.py @@ -0,0 +1,134 @@ +"""Tests for the (binary) calibration error metric.""" + +import array_api_compat.torch +import numpy.array_api as anp +import pytest +import torch + +from cyclops.evaluate.metrics.experimental import BinaryCalibrationError +from cyclops.evaluate.metrics.experimental.functional import binary_calibration_error + + +@pytest.mark.parametrize("xp", [anp, array_api_compat.torch]) +def test_binary_calibration_error_two_bins(xp): + """Test binary calibration error against a hand-computed example. + + target = [0, 1, 1, 0], preds = [0.1, 0.9, 0.8, 0.3], n_bins=2. + Bin [0, 0.5): preds=[0.1, 0.3], target=[0, 0] -> conf=0.2, acc=0.0, gap=0.2 + Bin [0.5, 1]: preds=[0.9, 0.8], target=[1, 1] -> conf=0.85, acc=1.0, gap=0.15 + ECE = 0.5 * 0.2 + 0.5 * 0.15 = 0.175 + """ + target = xp.asarray([0, 1, 1, 0]) + preds = xp.asarray([0.1, 0.9, 0.8, 0.3]) + result = binary_calibration_error(target, preds, n_bins=2) + assert float(result) == pytest.approx(0.175, abs=1e-4) + + +def test_binary_calibration_error_perfect_calibration(): + """A perfectly calibrated model must have (near) zero ECE.""" + target = anp.asarray([0, 0, 0, 0, 0, 1, 1, 1, 1, 1]) + preds = anp.asarray([0.0] * 5 + [1.0] * 5) + assert float(binary_calibration_error(target, preds, n_bins=2)) == pytest.approx( + 0.0, + ) + + +def test_binary_calibration_error_max_norm(): + """The 'max' norm must return the largest per-bin gap (MCE).""" + target = anp.asarray([0, 1, 1, 0]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + result = binary_calibration_error(target, preds, n_bins=2, norm="max") + assert float(result) == pytest.approx(0.2, abs=1e-4) + + +def test_binary_calibration_error_norms_ordering(): + """For a fixed input, max-norm gap must be >= l2 gap >= l1 (ECE) gap.""" + target = anp.asarray([0, 1, 1, 0, 1, 0, 0, 1]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.3, 0.4, 0.2, 0.6, 0.7]) + ece = float(binary_calibration_error(target, preds, n_bins=4, norm="l1")) + l2 = float(binary_calibration_error(target, preds, n_bins=4, norm="l2")) + mce = float(binary_calibration_error(target, preds, n_bins=4, norm="max")) + assert ece <= l2 <= mce + + +def test_binary_calibration_error_ignore_index(): + """Values matching ignore_index must be excluded from binning.""" + target = anp.asarray([0, 1, 1, -1]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.99]) + without_ignored = float( + binary_calibration_error( + anp.asarray([0, 1, 1]), + anp.asarray([0.1, 0.9, 0.8]), + n_bins=2, + ), + ) + with_ignored = float( + binary_calibration_error(target, preds, n_bins=2, ignore_index=-1), + ) + assert with_ignored == pytest.approx(without_ignored, abs=1e-5) + + +@pytest.mark.parametrize( + ("kwargs", "match"), + [ + ({"n_bins": 0}, "n_bins"), + ({"n_bins": -1}, "n_bins"), + ({"norm": "l3"}, "norm"), + ], +) +def test_binary_calibration_error_invalid_args(kwargs, match): + """Invalid n_bins or norm arguments must raise ValueError.""" + target = anp.asarray([0, 1]) + preds = anp.asarray([0.1, 0.9]) + with pytest.raises(ValueError, match=match): + binary_calibration_error(target, preds, **kwargs) + + +def test_binary_calibration_error_invalid_target_raises(): + """Non-binary target values must raise.""" + target = anp.asarray([0, 1, 2]) + preds = anp.asarray([0.1, 0.9, 0.8]) + with pytest.raises(RuntimeError): + binary_calibration_error(target, preds) + + +class TestBinaryCalibrationErrorClass: + """Tests for the BinaryCalibrationError metric class.""" + + def test_single_call(self): + """Test single-call usage matches the functional API.""" + target = anp.asarray([0, 1, 1, 0]) + preds = anp.asarray([0.1, 0.9, 0.8, 0.3]) + metric = BinaryCalibrationError(n_bins=2) + assert float(metric(target, preds)) == pytest.approx( + float(binary_calibration_error(target, preds, n_bins=2)), + ) + + def test_streaming_matches_batch(self): + """Accumulating bin counts over multiple updates must match a batch call.""" + target = [0, 1, 1, 0, 1, 0] + preds = [0.1, 0.9, 0.8, 0.3, 0.4, 0.2] + + batch_result = float( + binary_calibration_error( + anp.asarray(target), + anp.asarray(preds), + n_bins=4, + ), + ) + + metric = BinaryCalibrationError(n_bins=4) + for t, p in zip([target[:3], target[3:]], [preds[:3], preds[3:]]): + metric.update(anp.asarray(t), anp.asarray(p)) + streaming_result = float(metric.compute()) + + assert streaming_result == pytest.approx(batch_result, abs=1e-5) + + def test_torch_backend(self): + """Test the metric works with a torch backend.""" + target = torch.tensor([0, 1, 1, 0]) + preds = torch.tensor([0.1, 0.9, 0.8, 0.3]) + metric = BinaryCalibrationError(n_bins=2) + result = metric(target, preds) + assert isinstance(result, torch.Tensor) + assert float(result) == pytest.approx(0.175, abs=1e-4) diff --git a/tests/cyclops/evaluate/test_evaluator.py b/tests/cyclops/evaluate/test_evaluator.py new file mode 100644 index 000000000..4850740e4 --- /dev/null +++ b/tests/cyclops/evaluate/test_evaluator.py @@ -0,0 +1,188 @@ +"""Integration tests for the top-level evaluate() function. + +These exercise cyclops.evaluate.evaluator.evaluate() end-to-end against a +small in-memory dataset - previously the main public entry point for +evaluating models had zero test coverage. +""" + +import pytest +from datasets import Dataset, DatasetDict +from datasets.splits import Split + +from cyclops.data.slicer import SliceSpec +from cyclops.evaluate.evaluator import evaluate +from cyclops.evaluate.fairness.config import FairnessConfig +from cyclops.evaluate.metrics.experimental import BinaryAccuracy, BinaryPrecision +from cyclops.evaluate.metrics.experimental.metric_dict import MetricDict + + +@pytest.fixture +def classification_dataset() -> Dataset: + """Create a small synthetic binary classification dataset.""" + data = { + "target": [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1], + "prediction": [0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0], + "group": ["A"] * 6 + ["B"] * 6, + } + return Dataset.from_dict(data) + + +def test_evaluate_basic(classification_dataset): + """evaluate() with a single metric and no slicing computes an overall result.""" + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + ) + + assert "model_for_prediction" in results + overall = results["model_for_prediction"]["overall"] + assert 0 <= float(overall["BinaryAccuracy"]) <= 1 + assert overall["sample_size"] == classification_dataset.num_rows + + +def test_evaluate_with_slice_spec(classification_dataset): + """evaluate() with a slice_spec computes per-slice results.""" + metrics = MetricDict([BinaryAccuracy()]) + slice_spec = SliceSpec( + spec_list=[{"group": {"value": "A"}}, {"group": {"value": "B"}}], + ) + + results = evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + slice_spec=slice_spec, + ) + + model_results = results["model_for_prediction"] + assert set(model_results.keys()) == {"group:A", "group:B", "overall"} + assert model_results["group:A"]["sample_size"] == 6 + assert model_results["group:B"]["sample_size"] == 6 + assert model_results["overall"]["sample_size"] == 12 + + +def test_evaluate_multiple_prediction_columns(classification_dataset): + """evaluate() with multiple prediction columns computes results per model.""" + dataset = classification_dataset.add_column( + "prediction_2", + [1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0], + ) + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate( + dataset=dataset, + metrics=metrics, + target_columns="target", + prediction_columns=["prediction", "prediction_2"], + ) + + assert set(results.keys()) == {"model_for_prediction", "model_for_prediction_2"} + + +def test_evaluate_empty_slice_raises(classification_dataset): + """An empty slice must raise when raise_on_empty_slice=True.""" + metrics = MetricDict([BinaryAccuracy()]) + slice_spec = SliceSpec( + spec_list=[{"group": {"value": "nonexistent"}}], + include_overall=False, + ) + with pytest.raises(RuntimeError, match="empty"): + evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + slice_spec=slice_spec, + raise_on_empty_slice=True, + ) + + +def test_evaluate_empty_slice_warns_and_returns_nan(classification_dataset): + """An empty slice must warn and produce NaN metric values by default.""" + metrics = MetricDict([BinaryAccuracy()]) + slice_spec = SliceSpec( + spec_list=[{"group": {"value": "nonexistent"}}], + include_overall=False, + ) + with pytest.warns(RuntimeWarning, match="empty"): + results = evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + slice_spec=slice_spec, + raise_on_empty_slice=False, + ) + + slice_result = results["model_for_prediction"]["group:nonexistent"] + assert slice_result["BinaryAccuracy"] != slice_result["BinaryAccuracy"] # NaN + + +def test_evaluate_missing_required_column_raises(classification_dataset): + """A missing target/prediction column must raise ValueError.""" + metrics = MetricDict([BinaryAccuracy()]) + with pytest.raises(ValueError, match="missing_column"): + evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="missing_column", + prediction_columns="prediction", + ) + + +def test_evaluate_dataset_dict_without_split_uses_choose_split( + classification_dataset, +): + """A DatasetDict with split=None must fall back to choose_split(), not error.""" + dataset_dict = DatasetDict({"test": classification_dataset}) + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate( + dataset=dataset_dict, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + ) + assert "model_for_prediction" in results + + +def test_evaluate_dataset_dict_split_all_raises(classification_dataset): + """A DatasetDict with split=Split.ALL must raise ValueError.""" + dataset_dict = DatasetDict({"test": classification_dataset}) + metrics = MetricDict([BinaryAccuracy()]) + with pytest.raises(ValueError, match="Split.ALL"): + evaluate( + dataset=dataset_dict, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + split=Split.ALL, + ) + + +def test_evaluate_with_fairness_config(classification_dataset): + """evaluate() with a fairness_config populates a "fairness" results key.""" + metrics = MetricDict([BinaryAccuracy(), BinaryPrecision()]) + fairness_config = FairnessConfig( + metrics=metrics, + dataset=classification_dataset, # overridden by evaluate() with the real dataset + groups="group", + target_columns="target", + ) + + results = evaluate( + dataset=classification_dataset, + metrics=metrics, + target_columns="target", + prediction_columns="prediction", + fairness_config=fairness_config, + ) + + assert "fairness" in results + fairness_results = results["fairness"] + assert set(fairness_results.keys()) == {"group:A", "group:B", "overall"} + for group_result in fairness_results.values(): + assert "BinaryAccuracy" in group_result + assert "BinaryAccuracy Parity" in group_result diff --git a/tests/cyclops/evaluate/test_fairness_evaluator.py b/tests/cyclops/evaluate/test_fairness_evaluator.py new file mode 100644 index 000000000..06d947bc8 --- /dev/null +++ b/tests/cyclops/evaluate/test_fairness_evaluator.py @@ -0,0 +1,124 @@ +"""Integration tests for cyclops.evaluate.fairness.evaluator.evaluate_fairness(). + +evaluate_fairness() (989 lines) previously had zero test coverage despite +being the module's main entry point for fairness/subgroup evaluation. +""" + +import pytest +from datasets import Dataset + +from cyclops.evaluate.fairness.evaluator import evaluate_fairness +from cyclops.evaluate.metrics.experimental import BinaryAccuracy +from cyclops.evaluate.metrics.experimental.metric_dict import MetricDict + + +@pytest.fixture +def classification_dataset() -> Dataset: + """Create a small synthetic binary classification dataset.""" + data = { + "target": [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1], + "prediction": [0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0], + "group": ["A"] * 6 + ["B"] * 6, + "age": [20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75], + } + return Dataset.from_dict(data) + + +def test_evaluate_fairness_basic(classification_dataset): + """Basic categorical group fairness evaluation.""" + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate_fairness( + metrics=metrics, + dataset=classification_dataset, + groups="group", + target_columns="target", + prediction_columns="prediction", + ) + + assert set(results.keys()) == {"group:A", "group:B", "overall"} + for slice_result in results.values(): + assert "BinaryAccuracy" in slice_result + assert "BinaryAccuracy Parity" in slice_result + assert 0 <= float(slice_result["BinaryAccuracy"]) <= 1 + + # parity relative to the overall metric value must be 1.0 for "overall" itself + assert float(results["overall"]["BinaryAccuracy Parity"]) == pytest.approx(1.0) + + +def test_evaluate_fairness_group_base_values(classification_dataset): + """Parity must be computed relative to an explicit group_base_values.""" + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate_fairness( + metrics=metrics, + dataset=classification_dataset, + groups="group", + target_columns="target", + prediction_columns="prediction", + group_base_values={"group": "A"}, + ) + + accuracy_a = float(results["group:A"]["BinaryAccuracy"]) + parity_a = float(results["group:A"]["BinaryAccuracy Parity"]) + # base group's parity relative to itself must be 1.0 + assert parity_a == pytest.approx(1.0) + assert accuracy_a > 0 + + +def test_evaluate_fairness_group_bins_continuous(classification_dataset): + """Continuous groups must be bucketed via group_bins.""" + metrics = MetricDict([BinaryAccuracy()]) + results = evaluate_fairness( + metrics=metrics, + dataset=classification_dataset, + groups="age", + target_columns="target", + prediction_columns="prediction", + group_bins={"age": 3}, + ) + + assert "overall" in results + # binning into 3 groups should yield multiple non-overall slice keys + assert len(results) > 2 + for slice_result in results.values(): + if slice_result["sample_size"] > 0: + assert "BinaryAccuracy" in slice_result + + +def test_evaluate_fairness_invalid_dataset_type(): + """A non-Dataset `dataset` argument must raise TypeError.""" + metrics = MetricDict([BinaryAccuracy()]) + with pytest.raises(TypeError, match="Dataset"): + evaluate_fairness( + metrics=metrics, + dataset="not a dataset", # type: ignore[arg-type] + groups="group", + target_columns="target", + prediction_columns="prediction", + ) + + +def test_evaluate_fairness_missing_group_column_raises(classification_dataset): + """A missing group column must raise ValueError.""" + metrics = MetricDict([BinaryAccuracy()]) + with pytest.raises(ValueError, match="missing_group"): + evaluate_fairness( + metrics=metrics, + dataset=classification_dataset, + groups="missing_group", + target_columns="target", + prediction_columns="prediction", + ) + + +def test_evaluate_fairness_invalid_array_lib(classification_dataset): + """An unsupported array_lib must raise NotImplementedError.""" + metrics = MetricDict([BinaryAccuracy()]) + with pytest.raises(NotImplementedError): + evaluate_fairness( + metrics=metrics, + dataset=classification_dataset, + groups="group", + target_columns="target", + prediction_columns="prediction", + array_lib="not_a_real_lib", # type: ignore[arg-type] + ) diff --git a/tests/cyclops/models/neural_nets/__init__.py b/tests/cyclops/models/neural_nets/__init__.py new file mode 100644 index 000000000..bccfe9fbf --- /dev/null +++ b/tests/cyclops/models/neural_nets/__init__.py @@ -0,0 +1 @@ +"""Cyclops models neural_nets test package.""" diff --git a/tests/cyclops/models/neural_nets/test_mlp.py b/tests/cyclops/models/neural_nets/test_mlp.py new file mode 100644 index 000000000..e6e0d52ef --- /dev/null +++ b/tests/cyclops/models/neural_nets/test_mlp.py @@ -0,0 +1,58 @@ +"""Tests for the MLPModel.""" + +import torch + +from cyclops.models.catalog import create_model +from cyclops.models.neural_nets.mlp import MLPModel + + +def test_mlp_model_forward_pass(): + """MLPModel must be constructible and runnable with default arguments. + + Regression test: get_module("activation", activation) returns the + activation *class* (e.g. torch.nn.ReLU), not an instance, and the + first hidden layer was wrapped in an extra list + (`layers = [self._layer(...)]` instead of `self._layer(...)`), both + of which made `nn.Sequential(*layers)` raise a TypeError. + """ + model = MLPModel(input_dim=10) + output = model(torch.randn(4, 10)) + assert output.shape == (4, 1) + + +def test_mlp_model_multiple_hidden_layers(): + """Hidden-to-hidden layer dimensions must chain correctly. + + Regression test: the loop connecting hidden layers used `input_dim` + instead of `hidden_dims[i]` for the first hidden-to-hidden + connection, causing a shape mismatch whenever hidden_dims[0] != + input_dim. + """ + model = MLPModel(input_dim=10, hidden_dims=(32, 16, 8), output_dim=2) + linear_layers = [m for m in model.model if isinstance(m, torch.nn.Linear)] + dims = [(layer.in_features, layer.out_features) for layer in linear_layers] + assert dims == [(10, 32), (32, 16), (16, 8), (8, 2)] + + output = model(torch.randn(4, 10)) + assert output.shape == (4, 2) + + +def test_mlp_model_accepts_module_instance_as_activation(): + """Activation may be passed as an already-instantiated nn.Module.""" + model = MLPModel(input_dim=10, activation=torch.nn.Tanh()) + assert isinstance(model.activation, torch.nn.Tanh) + output = model(torch.randn(4, 10)) + assert output.shape == (4, 1) + + +def test_mlp_pt_config_initializes(): + """The packaged mlp_pt config must initialize without error. + + Regression test: configs/mlp_pt.yaml set model__layer_dim, a + leftover from the RNN/GRU/LSTM configs, which MLPModel.__init__ + doesn't accept. + """ + wrapped_model = create_model("mlp_pt", model__input_dim=10) + wrapped_model.initialize() + output = wrapped_model.model_(torch.randn(4, 10)) + assert output.shape == (4, 1) diff --git a/tests/cyclops/monitor/test_detector.py b/tests/cyclops/monitor/test_detector.py index f78bd98cb..0c87967c7 100644 --- a/tests/cyclops/monitor/test_detector.py +++ b/tests/cyclops/monitor/test_detector.py @@ -7,6 +7,7 @@ synthetic_nih_dataset, ) +from cyclops.data.slicer import SliceSpec from cyclops.monitor.detector import Detector from cyclops.monitor.reductor import Reductor from cyclops.monitor.tester import TSTester @@ -51,3 +52,73 @@ def test_detector_pca_mmd(source_target): ds_source, ds_target = source_target results = detector.detect_shift(ds_source, ds_target) assert results["p_val"].shape == (2, 3) + + +def test_detector_detect_shift_by_subgroup(source_target): + """Test Detector.detect_shift_by_subgroup.""" + reductor = Reductor( + "pca", + n_components=2, + feature_columns=[f"feature_{i}" for i in range(10)], + ) + tester = TSTester("mmd") + detector = Detector("sensitivity_test", reductor, tester) + ds_source, ds_target = source_target + detector.fit(ds_source) + + slice_spec = SliceSpec( + spec_list=[{"mortality": {"value": 0}}, {"mortality": {"value": 1}}], + include_overall=False, + ) + results = detector.detect_shift_by_subgroup(ds_target, slice_spec) + + assert set(results.keys()) == set(slice_spec.get_slices().keys()) + for subgroup_result in results.values(): + assert subgroup_result["sample_size"] > 0 + assert 0 <= subgroup_result["p_val"] <= 1 + assert subgroup_result["shift_detected"] in (0, 1) + + +def test_detector_detect_shift_by_subgroup_small_subgroup_skipped(source_target): + """Subgroups below min_sample_size must be skipped, not tested.""" + reductor = Reductor( + "pca", + n_components=2, + feature_columns=[f"feature_{i}" for i in range(10)], + ) + tester = TSTester("mmd") + detector = Detector("sensitivity_test", reductor, tester) + ds_source, ds_target = source_target + detector.fit(ds_source) + + slice_spec = SliceSpec( + spec_list=[{"mortality": {"value": 0}}], + include_overall=False, + ) + results = detector.detect_shift_by_subgroup( + ds_target, + slice_spec, + min_sample_size=10_000, + ) + + (subgroup_result,) = results.values() + assert subgroup_result["p_val"] is None + assert subgroup_result["shift_detected"] is None + assert subgroup_result["sample_size"] < 10_000 + + +def test_detector_detect_shift_by_subgroup_invalid_correction(source_target): + """An unknown correction method must raise a clear error.""" + reductor = Reductor( + "pca", + n_components=2, + feature_columns=[f"feature_{i}" for i in range(10)], + ) + tester = TSTester("mmd") + detector = Detector("sensitivity_test", reductor, tester) + ds_source, ds_target = source_target + detector.fit(ds_source) + + slice_spec = SliceSpec(spec_list=[{"mortality": {"value": 0}}]) + with pytest.raises(ValueError, match="correction"): + detector.detect_shift_by_subgroup(ds_target, slice_spec, correction="invalid") diff --git a/tests/cyclops/monitor/test_tester.py b/tests/cyclops/monitor/test_tester.py index 482b3140d..3790c67a7 100644 --- a/tests/cyclops/monitor/test_tester.py +++ b/tests/cyclops/monitor/test_tester.py @@ -88,3 +88,40 @@ def test_dctester(source_target, generic_source_target, method): tester.fit(X_source) p_val = tester.test_shift(X_target)[0] assert 0 <= p_val <= 1 + + +def test_dctester_explain_shift(source_target): + """Test DCTester.explain_shift for the classifier tester method.""" + X_source, X_target = source_target + model = RandomForestClassifier() + tester = DCTester("classifier", model=model) + tester.fit(X_source) + tester.test_shift(X_target) + + feature_names = [f"feature_{i}" for i in range(X_source.shape[1])] + importances = tester.explain_shift(X_target, feature_names=feature_names) + + assert set(importances.keys()) == set(feature_names) + assert all(value >= 0 for value in importances.values()) + # sorted by descending importance + values = list(importances.values()) + assert values == sorted(values, reverse=True) + + +def test_dctester_explain_shift_unsupported_method(source_target): + """explain_shift must raise a clear error for non-classifier methods.""" + X_source, X_target = source_target + tester = DCTester("spot_the_diff") + tester.fit(X_source) + tester.test_shift(X_target) + + with pytest.raises(ValueError, match="classifier"): + tester.explain_shift(X_target) + + +def test_dctester_explain_shift_before_fit(): + """explain_shift must raise a clear error if called before fit/test_shift.""" + model = RandomForestClassifier() + tester = DCTester("classifier", model=model) + with pytest.raises(ValueError, match="fit"): + tester.explain_shift(np.random.rand(10, 10)) diff --git a/tests/cyclops/report/test_report.py b/tests/cyclops/report/test_report.py index 7ffc85783..2fc5e762a 100644 --- a/tests/cyclops/report/test_report.py +++ b/tests/cyclops/report/test_report.py @@ -1,5 +1,6 @@ """Test cyclops report module model report.""" +import os from unittest import TestCase import numpy as np @@ -371,6 +372,47 @@ def test_export(self): assert isinstance(report_path, str) +def test_export_default_filename_is_timestamped_per_call(tmp_path): + """Repeated export() calls without output_filename must not overwrite each other. + + Regression test: the default output filename used to be the static + "model_card.html"/"model_card.json", so every export() call into the + same output_dir silently overwrote the previous report, defeating the + trend/history comparison the export() docstring promises. + """ + report = ModelCardReport(str(tmp_path)) + report.log_owner(name="John Doe") + + path_1 = report.export( + interactive=False, + save_json=True, + synthetic_timestamp="2024-01-01 00:00:00", + ) + path_2 = report.export( + interactive=False, + save_json=True, + synthetic_timestamp="2024-01-02 00:00:00", + ) + + assert path_1 != path_2 + assert os.path.exists(path_1) + assert os.path.exists(path_2) + + +def test_export_with_no_performance_metrics(tmp_path): + """Test that export() does not crash when no PerformanceMetric was logged. + + Regression test: previously raised IndexError because + `current_report_metrics[0]` was indexed unconditionally on a + possibly-empty list. + """ + report = ModelCardReport(str(tmp_path)) + report.log_owner(name="John Doe") + + report_path = report.export(interactive=False, save_json=False) + assert isinstance(report_path, str) + + def test_log_performance_metrics(): """Test log_performance_metrics.""" report = ModelCardReport() diff --git a/tests/cyclops/report/test_utils.py b/tests/cyclops/report/test_utils.py index 200a022bd..bbab8eaec 100644 --- a/tests/cyclops/report/test_utils.py +++ b/tests/cyclops/report/test_utils.py @@ -21,6 +21,7 @@ QuantitativeAnalysis, ) from cyclops.report.utils import ( + _process_metric_name, create_metric_card_plot, create_metric_cards, extract_performance_metrics, @@ -305,6 +306,22 @@ def model_card(): return model_card +def test_process_metric_name_with_recognized_prefix(): + """Test _process_metric_name strips known Binary/Multiclass/Multilabel prefixes.""" + assert _process_metric_name({"type": "BinaryAccuracy"}) == "Accuracy" + assert _process_metric_name({"type": "MulticlassPrecision"}) == "Precision" + assert _process_metric_name({"type": "MultilabelRecall"}) == "Recall" + + +def test_process_metric_name_with_unrecognized_prefix(): + """A metric type without a Binary/Multiclass/Multilabel prefix must not crash. + + Regression test: previously raised UnboundLocalError because `name` was + only assigned inside the prefix-matching branches. + """ + assert _process_metric_name({"type": "CustomMetric"}) == "CustomMetric" + + def test_sweep_tests(model_card): """Test sweep_tests function.""" tests = [] diff --git a/tests/cyclops/utils/test_file.py b/tests/cyclops/utils/test_file.py index c85ef3912..57d15e719 100644 --- a/tests/cyclops/utils/test_file.py +++ b/tests/cyclops/utils/test_file.py @@ -156,6 +156,17 @@ def test_exchange_extension(): assert exchange_extension("/tmp/file.txt", "csv") == "/tmp/file.csv" +def test_exchange_extension_no_existing_extension(): + """Test exchange_extension fn on a path with no existing extension. + + Regression test: os.path.splitext returns "" for old_ext on an + extensionless path, and `file_path[:-len(old_ext)]` evaluated to + `file_path[:-0]` == `file_path[:0]` == "", silently dropping the + filename instead of appending the new extension. + """ + assert exchange_extension("/tmp/myfile", "csv") == "/tmp/myfile.csv" + + def test_process_file_save_path(): """Test process_file_save_path fn.""" with pytest.raises(ValueError): diff --git a/tests/cyclops/utils/test_index.py b/tests/cyclops/utils/test_index.py index f189fd1d4..5c3862fe6 100644 --- a/tests/cyclops/utils/test_index.py +++ b/tests/cyclops/utils/test_index.py @@ -15,7 +15,7 @@ def test_index_axis(): indices = index_axis(4, 2, (10, 20, 30)) assert indices[0] == slice(None, None, None) - assert indices[0] == slice(None, None, None) + assert indices[1] == slice(None, None, None) assert indices[2] == 4