API: Metrics & helpers

The metrics take probabilities of shape (B, C, *spatial) and integer labels of shape (B, *spatial), or floating class distributions with the same shape as the probabilities. mask selects valid voxels for either target kind; ignore_index applies to hard labels. See multi-rater calibration.

negative_log_likelihood

negative_log_likelihood(probs: Any, targets: Any, mask: Any | None = None, ignore_index: int = -100) -> float

Mean negative log-likelihood (cross-entropy) over valid voxels.

Zero probabilities are floored at 1e-12 by :func:fiducio.utils.safe_log, capping the NLL at ~27.63 instead of inf.

Source code in src/fiducio/metrics/calibration.py

def negative_log_likelihood(
    probs: Any,
    targets: Any,
    mask: Any | None = None,
    ignore_index: int = -100,
) -> float:
    """Mean negative log-likelihood (cross-entropy) over valid voxels.

    Zero probabilities are floored at ``1e-12`` by :func:`fiducio.utils.safe_log`,
    capping the NLL at ~27.63 instead of ``inf``.
    """
    p, y, valid = _voxels(probs, targets, mask, ignore_index)
    if y.shape == p.shape:
        loss = torch.zeros_like(_select(p[:, 0], valid), dtype=torch.float64)
        for c in range(p.shape[1]):
            loss -= _select(y[:, c], valid).double() * safe_log(_select(p[:, c], valid).double())
        return float(loss.mean().item())
    true_p = _select(_true_class_probs(p, y), valid).double()
    if true_p.numel() == 0:
        return float("nan")
    return float((-safe_log(true_p)).mean().item())

expected_calibration_error

expected_calibration_error(probs: Any, targets: Any, mask: Any | None = None, ignore_index: int = -100, n_bins: int = 15, *, class_index: int | None = None) -> float

Top-1 expected calibration error (ECE) with uniform binning.

Pass class_index for class-specific calibration against hard or soft targets. This measures ECE against the supplied targets, not mean rater ECE.

The confidence is the maximum predicted probability and the accuracy is whether the argmax matches the label. Bins partition [0, 1] uniformly and each bin's gap is weighted by its share of voxels — see :func:unweighted_calibration_error for the unweighted variant.

n_bins defaults to 15. The paper reports ECE with n_bins=50 (ACE with n_bins=15) and averages metrics per image before pooling; pass n_bins=50 and compute per case to follow that convention. Boundary-aware ECE is not included.

Source code in src/fiducio/metrics/calibration.py

def expected_calibration_error(
    probs: Any,
    targets: Any,
    mask: Any | None = None,
    ignore_index: int = -100,
    n_bins: int = 15,
    *,
    class_index: int | None = None,
) -> float:
    """Top-1 expected calibration error (ECE) with uniform binning.

    Pass ``class_index`` for class-specific calibration against hard or soft
    targets. This measures ECE against the supplied targets, not mean rater ECE.

    The confidence is the maximum predicted probability and the accuracy is
    whether the argmax matches the label. Bins partition ``[0, 1]`` uniformly
    and each bin's gap is weighted by its share of voxels — see
    :func:`unweighted_calibration_error` for the unweighted variant.

    ``n_bins`` defaults to 15. The paper reports ECE with ``n_bins=50`` (ACE with
    ``n_bins=15``) and averages metrics per image before pooling; pass
    ``n_bins=50`` and compute per case to follow that convention. Boundary-aware
    ECE is not included.
    """
    return reliability_curve(
        probs, targets, mask, ignore_index, n_bins, class_index=class_index
    ).ece

unweighted_calibration_error

unweighted_calibration_error(probs: Any, targets: Any, mask: Any | None = None, ignore_index: int = -100, n_bins: int = 15, *, class_index: int | None = None) -> float

Top-1 unweighted calibration error (ACE) with uniform binning.

class_index selects a class instead of the default top-1 confidence.

The unweighted mean of the per-bin |confidence - accuracy| gap over non-empty uniform bins. This is not the adaptive calibration error of other work, which uses data-dependent (equal-mass) bins.

Uses the same uniform confidence bins as :func:expected_calibration_error, but weights each non-empty bin equally instead of by its share of voxels; a confidence region visited by only a handful of voxels therefore counts as much as a densely populated one. n_bins=15 matches the paper's ACE; the paper computes ACE on the pooled test voxels rather than per image. Kept as average_calibration_error for backwards compatibility.

Source code in src/fiducio/metrics/calibration.py

def unweighted_calibration_error(
    probs: Any,
    targets: Any,
    mask: Any | None = None,
    ignore_index: int = -100,
    n_bins: int = 15,
    *,
    class_index: int | None = None,
) -> float:
    """Top-1 unweighted calibration error (ACE) with uniform binning.

    ``class_index`` selects a class instead of the default top-1 confidence.

    The unweighted mean of the per-bin ``|confidence - accuracy|`` gap over
    non-empty uniform bins. This is **not** the adaptive calibration error of
    other work, which uses data-dependent (equal-mass) bins.

    Uses the same uniform confidence bins as :func:`expected_calibration_error`,
    but weights each non-empty bin equally instead of by its share of voxels; a
    confidence region visited by only a handful of voxels therefore counts as
    much as a densely populated one. ``n_bins=15`` matches the paper's ACE; the
    paper computes ACE on the pooled test voxels rather than per image. Kept as
    ``average_calibration_error`` for backwards compatibility.
    """
    return reliability_curve(
        probs, targets, mask, ignore_index, n_bins, class_index=class_index
    ).ace

Alias

average_calibration_error is a backwards-compatible alias of unweighted_calibration_error — the same function, not a separate estimator.

brier_score

brier_score(probs: Any, targets: Any, mask: Any | None = None, ignore_index: int = -100) -> float

Mean sum of squared class residuals against labels or soft distributions.

Against an MHR, this differs from the mean per-rater Brier by the rater disagreement term. Binary two-channel Brier is twice positive-channel MSE. Residuals accumulate in float64 to retain errors near perfect predictions.

Source code in src/fiducio/metrics/calibration.py

def brier_score(
    probs: Any,
    targets: Any,
    mask: Any | None = None,
    ignore_index: int = -100,
) -> float:
    """Mean sum of squared class residuals against labels or soft distributions.

    Against an MHR, this differs from the mean per-rater Brier by the rater
    disagreement term. Binary two-channel Brier is twice positive-channel MSE.
    Residuals accumulate in float64 to retain errors near perfect predictions.
    """
    p, y, valid = _voxels(probs, targets, mask, ignore_index)
    if valid is not None and not bool(valid.any()):
        return float("nan")
    # Direct residuals avoid cancellation near perfect predictions, one channel at a time.
    squared = torch.zeros_like(_select(p[:, 0], valid), dtype=torch.float64)
    for c in range(p.shape[1]):
        target = y[:, c] if y.shape == p.shape else y == c
        residual = _select(p[:, c], valid).double() - _select(target, valid).double()
        squared += residual.square()
    return float(squared.mean().item())

reliability_curve

reliability_curve(probs: Any, targets: Any, mask: Any | None = None, ignore_index: int = -100, n_bins: int = 15, *, class_index: int | None = None) -> ReliabilityCurve

Compute top-1 or class-specific reliability with uniform binning.

The confidence is the maximum predicted probability and the accuracy is whether the argmax matches the label, or q[argmax(p)] for soft targets. With class_index=c, compare p_c with 1[y=c] or q_c. ECE against MHR differs from mean per-rater ECE; retain the individual votes to compute the latter. Useful both for reporting ECE/ACE and for drawing reliability diagrams (see :func:fiducio.plots.reliability_diagram).

Source code in src/fiducio/metrics/calibration.py

def reliability_curve(
    probs: Any,
    targets: Any,
    mask: Any | None = None,
    ignore_index: int = -100,
    n_bins: int = 15,
    *,
    class_index: int | None = None,
) -> ReliabilityCurve:
    """Compute top-1 or class-specific reliability with uniform binning.

    The confidence is the maximum predicted probability and the accuracy is
    whether the argmax matches the label, or ``q[argmax(p)]`` for soft targets.
    With ``class_index=c``, compare ``p_c`` with ``1[y=c]`` or ``q_c``.
    ECE against MHR differs from mean per-rater ECE; retain the individual votes
    to compute the latter. Useful both for reporting ECE/ACE and
    for drawing reliability diagrams (see
    :func:`fiducio.plots.reliability_diagram`).
    """
    try:  # operator.index accepts numpy integers
        n_bins = -1 if isinstance(n_bins, bool) else operator.index(n_bins)
    except TypeError:
        n_bins = -1
    if n_bins < 1:
        raise ValueError("n_bins must be an integer >= 1")
    p, y, valid = _voxels(probs, targets, mask, ignore_index)
    if class_index is not None:
        try:
            class_index = -1 if isinstance(class_index, bool) else operator.index(class_index)
        except TypeError:
            class_index = -1
        if not 0 <= class_index < p.shape[1]:
            raise ValueError(f"class_index must be an integer in [0, {p.shape[1] - 1}]")
    return _binned_curve(p, y, valid, n_bins, class_index)

ReliabilityCurve dataclass

ReliabilityCurve(bin_edges: Tensor, bin_confidence: Tensor, bin_accuracy: Tensor, bin_counts: Tensor, ece: float, ace: float)

Per-bin top-1 or class-specific reliability statistics.

Attributes:

Name	Type	Description
bin_edges	Tensor	(n_bins + 1,) bin boundaries in [0, 1].
bin_confidence	Tensor	(n_bins,) mean predicted confidence per bin (0 for empty bins).
bin_accuracy	Tensor	(n_bins,) mean agreement or class proportion per bin (0 for empty bins).
bin_counts	Tensor	(n_bins,) int64 number of voxels per bin. Other statistics are float64; all tensors reside on the predictions' device.
ece	float	Expected calibration error (count-weighted mean |confidence - accuracy| over bins).
ace	float	Unweighted calibration error (unweighted mean |confidence - accuracy| over non-empty bins). Unlike ece, a sparsely populated bin counts as much as a densely populated one.

Diagnostics

See reliability diagnostics.

diagnose

diagnose(logits: Any, labels: Any, *, input_type: str = 'logits', mask: Any | None = None, ignore_index: int = -100, device: Any | None = None, n_bins: int = 15, class_names: list[str] | None = None, image_ids: list[str] | None = None, background_class: int | None = None, dice_failure_threshold: float = 0.1, n_worst: int = 5, n_boot: int = 200, seed: int = 0, min_images: int = 20, min_class_fraction: float = 0.01, min_class_voxels: int = 1000) -> ReliabilityReport

Diagnose the reliability of a segmenter on a calibration set.

Parameters:

Name	Type	Description	Default
logits	Any	(B, C, *spatial) logits (or probabilities with input_type="probs"); B is the number of images or cases. Alternatively a torch.nn.Module: labels is then a DataLoader and each batch is (inputs, labels) or (inputs, labels, mask). inputs are whatever the model accepts (moved to device); the model output must have shape (B, C, *spatial). The model runs under torch.no_grad() and its training flag is unchanged (call model.eval() yourself if you want inference-mode dropout/batchnorm). Only per-image statistics are retained, so the loader's images are never materialised together.	required
labels	Any	(B, *spatial) integer labels, or (B, C, *spatial) floating class distributions (soft targets / MHR, same convention as the metrics). With a model, the DataLoader instead of labels.	required
mask	Any | None	Voxel selection, as in the metrics. For soft targets only mask applies. With the DataLoader path, pass per-batch masks in the batch instead.	None
ignore_index	Any | None	Voxel selection, as in the metrics. For soft targets only mask applies. With the DataLoader path, pass per-batch masks in the batch instead.	None
device	Any | None	Device for the DataLoader path's batches (default: the device of the model's first parameter, otherwise CPU). Ignored for tensor inputs.	None
class_names	list[str] | None	Optional display names (default: indices).	None
image_ids	list[str] | None	Optional display names (default: indices).	None
background_class	int | None	Class excluded from the per-image mean Dice, whose value near 1 would hide a failure on the classes that matter. None excludes nothing.	None
dice_failure_threshold	float	A class fails on an image when it is in the labels and its Dice is below this. An image fails when its mean Dice over the classes present is below it.	0.1
n_worst	int	Number of lowest-Dice images listed.	5
n_boot	int	Bootstrap resamples over images for the 95% intervals (0 disables them; an interval needs at least two resamples, so 1 yields NaN) and the seed of the resampling. Per-class ECE, routed ECE and Dice intervals reuse the same resampling.	200
seed	int	Bootstrap resamples over images for the 95% intervals (0 disables them; an interval needs at least two resamples, so 1 yields NaN) and the seed of the resampling. Per-class ECE, routed ECE and Dice intervals reuse the same resampling.	200
min_images	int	Thresholds of the warnings. They are heuristics, not statistical guarantees.	20
min_class_fraction	int	Thresholds of the warnings. They are heuristics, not statistical guarantees.	20
min_class_voxels	int	Thresholds of the warnings. They are heuristics, not statistical guarantees.	20

Notes

Soft targets. Per-class ECE compares p_c with the target mass q_c and the top-1/routed ECE compares the top-1 confidence with q[argmax(p)]; NLL and Brier follow :mod:fiducio.metrics.calibration. Hard-label Dice and the failure statistics are undefined: every Dice field, the image Dice matrix and the per-image failure counts are NaN/0, and image_dice is all NaN. Class balance still counts the expected label mass sum(q_c).

Boundary bins. report.boundary describes the first [0, 1/n_bins] and last [(n_bins-1)/n_bins, 1] top-1 confidence bins: their voxel mass, their mean confidence minus accuracy, and their exact contribution to ece. It is additive reporting only; ece keeps its definition.

The ECE interval is a percentile bootstrap. ECE has a finite-sample noise floor that resampling raises further, so for a well-calibrated model the interval can sit slightly above the estimate (it missed it in 15 of 100 simulated calibrated models of 24 images). Read it as an order of magnitude. The bias line uses the voxel-weighted confidence-accuracy gap, which has no such noise floor, but its intervals are slightly narrow: it flagged 9 of those 100 calibrated models, against 5 expected (reproducibility/diagnose_interval_check.py reruns the simulation).

Cost is 2C + 1 reliability curves per image (global, per-class and routed), plus a few passes for NLL, Brier and the confusion matrix. With a DataLoader the model forward pass dominates; the accumulator keeps O(n_images) summaries (O(n_images * C * n_bins)), not the voxels. Classes absent from both labels and prediction have an undefined Dice and are never counted as failures.

Source code in src/fiducio/diagnostics.py

def diagnose(
    logits: Any,
    labels: Any,
    *,
    input_type: str = "logits",
    mask: Any | None = None,
    ignore_index: int = -100,
    device: Any | None = None,
    n_bins: int = 15,
    class_names: list[str] | None = None,
    image_ids: list[str] | None = None,
    background_class: int | None = None,
    dice_failure_threshold: float = 0.1,
    n_worst: int = 5,
    n_boot: int = 200,
    seed: int = 0,
    min_images: int = 20,
    min_class_fraction: float = 0.01,
    min_class_voxels: int = 1000,
) -> ReliabilityReport:
    """Diagnose the reliability of a segmenter on a calibration set.

    Parameters
    ----------
    logits:
        ``(B, C, *spatial)`` logits (or probabilities with ``input_type="probs"``);
        ``B`` is the number of images or cases. Alternatively a ``torch.nn.Module``:
        ``labels`` is then a ``DataLoader`` and each batch is ``(inputs, labels)`` or
        ``(inputs, labels, mask)``. ``inputs`` are whatever the model accepts (moved
        to ``device``); the model output must have shape ``(B, C, *spatial)``. The
        model runs under ``torch.no_grad()`` and its ``training`` flag is unchanged
        (call ``model.eval()`` yourself if you want inference-mode dropout/batchnorm).
        Only per-image statistics are retained, so the loader's images are never
        materialised together.
    labels:
        ``(B, *spatial)`` integer labels, or ``(B, C, *spatial)`` floating class
        distributions (soft targets / MHR, same convention as the metrics). With a
        model, the DataLoader instead of labels.
    mask, ignore_index:
        Voxel selection, as in the metrics. For soft targets only ``mask`` applies.
        With the DataLoader path, pass per-batch masks in the batch instead.
    device:
        Device for the DataLoader path's batches (default: the device of the model's
        first parameter, otherwise CPU). Ignored for tensor inputs.
    class_names, image_ids:
        Optional display names (default: indices).
    background_class:
        Class excluded from the per-image mean Dice, whose value near 1 would hide a
        failure on the classes that matter. ``None`` excludes nothing.
    dice_failure_threshold:
        A class fails on an image when it is in the labels and its Dice is below this.
        An image fails when its mean Dice over the classes present is below it.
    n_worst:
        Number of lowest-Dice images listed.
    n_boot, seed:
        Bootstrap resamples over images for the 95% intervals (``0`` disables them;
        an interval needs at least two resamples, so ``1`` yields NaN) and the seed
        of the resampling. Per-class ECE, routed ECE and Dice intervals reuse the same
        resampling.
    min_images, min_class_fraction, min_class_voxels:
        Thresholds of the warnings. They are heuristics, not statistical guarantees.

    Notes
    -----
    Soft targets. Per-class ECE compares ``p_c`` with the target mass ``q_c`` and the
    top-1/routed ECE compares the top-1 confidence with ``q[argmax(p)]``; NLL and
    Brier follow :mod:`fiducio.metrics.calibration`. Hard-label Dice and the failure
    statistics are undefined: every Dice field, the image Dice matrix and the
    per-image failure counts are NaN/0, and ``image_dice`` is all NaN. Class balance
    still counts the expected label mass ``sum(q_c)``.

    Boundary bins. ``report.boundary`` describes the first ``[0, 1/n_bins]`` and last
    ``[(n_bins-1)/n_bins, 1]`` top-1 confidence bins: their voxel mass, their mean
    confidence minus accuracy, and their exact contribution to ``ece``. It is
    additive reporting only; ``ece`` keeps its definition.

    The ECE interval is a percentile bootstrap. ECE has a finite-sample noise floor that
    resampling raises further, so for a well-calibrated model the interval can sit
    slightly above the estimate (it missed it in 15 of 100 simulated calibrated models
    of 24 images). Read it as an order of magnitude. The *bias* line uses the
    voxel-weighted confidence-accuracy gap, which has no such noise floor, but its
    intervals are slightly narrow: it flagged 9 of those 100 calibrated models, against
    5 expected (``reproducibility/diagnose_interval_check.py`` reruns the simulation).

    Cost is ``2C + 1`` reliability curves per image (global, per-class and routed),
    plus a few passes for NLL, Brier and the confusion matrix. With a DataLoader the
    model forward pass dominates; the accumulator keeps ``O(n_images)`` summaries
    (``O(n_images * C * n_bins)``), not the voxels. Classes absent from both labels and
    prediction have an undefined Dice and are never counted as failures.
    """
    if isinstance(ignore_index, bool) or not isinstance(ignore_index, numbers.Integral):
        raise TypeError(f"ignore_index must be an int, got {ignore_index!r}")
    ignore_index = int(ignore_index)
    if input_type not in ("logits", "probs"):
        raise ValueError(f"input_type must be 'logits' or 'probs', got {input_type!r}")
    for name, value in (
        ("n_boot", n_boot),
        ("n_worst", n_worst),
        ("min_images", min_images),
        ("min_class_voxels", min_class_voxels),
    ):
        if isinstance(value, bool) or not isinstance(value, numbers.Integral) or value < 0:
            raise ValueError(f"{name} must be an integer >= 0")
    if isinstance(seed, bool) or not isinstance(seed, numbers.Integral):
        raise ValueError("seed must be an integer")
    seed = int(seed)
    if (
        isinstance(min_class_fraction, bool)
        or not isinstance(min_class_fraction, numbers.Real)
        or not 0.0 <= min_class_fraction <= 1.0
    ):
        raise ValueError("min_class_fraction must be a number in [0, 1]")
    if (
        isinstance(dice_failure_threshold, bool)
        or not isinstance(dice_failure_threshold, numbers.Real)
        or not 0.0 <= dice_failure_threshold <= 1.0
    ):
        raise ValueError("dice_failure_threshold must be a number in [0, 1]")

    acc: _Accumulator | None = None
    if isinstance(logits, torch.nn.Module):
        if labels is None:
            raise ValueError("labels must be a DataLoader when logits is a torch.nn.Module")
        if isinstance(labels, (torch.Tensor, np.ndarray)):
            raise ValueError(
                "labels must be a DataLoader when logits is a torch.nn.Module, not an array"
            )
        if mask is not None:
            raise ValueError("with a model and DataLoader, pass each mask in the batch")
        model = logits
        if device is None:
            first_param = next(iter(model.parameters()), None)
            device = first_param.device if first_param is not None else torch.device("cpu")
        batch_device = torch.device(device)
        for batch in labels:
            try:
                items = list(batch)
            except TypeError as exc:
                raise ValueError(
                    "each loader batch must be (inputs, labels) or (inputs, labels, mask)"
                ) from exc
            if not 2 <= len(items) <= 3:
                raise ValueError(
                    "each loader batch must be (inputs, labels) or (inputs, labels, mask)"
                )
            inputs = _move_to_device(items[0], batch_device)
            target = _move_to_device(items[1], batch_device)
            batch_mask = _move_to_device(items[2], batch_device) if len(items) == 3 else None
            with torch.no_grad():
                output = model(inputs)
            if not torch.is_tensor(output):
                raise TypeError("the model must return a tensor of shape (B, C, *spatial)")
            xb = output.detach()
            if not xb.is_floating_point():
                xb = xb.float()
            if xb.ndim < 3:
                raise ValueError(
                    "the model output must have shape (B, C, *spatial) with at least one "
                    f"spatial axis; got {tuple(xb.shape)}"
                )
            mb = None if batch_mask is None else to_tensor(batch_mask, device=xb.device).bool()
            yb, num_classes = _prepare_batch(
                xb,
                target,
                mb,
                input_type,
                ignore_index,
                None if acc is None else acc.num_classes,
            )
            if acc is None:
                acc = _Accumulator(num_classes, n_bins, input_type, ignore_index)
            acc.add(xb, yb, mb)
        if acc is None or acc.n_images == 0:
            raise ValueError("at least one image is required")
    else:
        x = to_tensor(logits).detach()
        if not x.is_floating_point():
            x = x.float()
        m = None if mask is None else to_tensor(mask, device=x.device).bool()
        if x.ndim < 3:
            raise ValueError(
                f"logits must have shape (B, C, *spatial) with at least one spatial axis; got {tuple(x.shape)}"
            )
        if x.shape[0] == 0:
            raise ValueError("at least one image is required")
        y, num_classes = _prepare_batch(x, labels, m, input_type, ignore_index)
        acc = _Accumulator(num_classes, n_bins, input_type, ignore_index)
        acc.add(x, y, m)

    num_classes, n_images = acc.num_classes, acc.n_images
    names = (
        _string_sequence(class_names, "class_names")
        if class_names is not None
        else [str(c) for c in range(num_classes)]
    )
    ids = (
        _string_sequence(image_ids, "image_ids")
        if image_ids is not None
        else [str(i) for i in range(n_images)]
    )
    if len(names) != num_classes or len(ids) != n_images:
        raise ValueError("class_names must have one name per class and image_ids one per image")
    if background_class is not None:
        if isinstance(background_class, bool) or not isinstance(background_class, numbers.Integral):
            raise ValueError(f"background_class must be an integer in [0, {num_classes - 1}]")
        background_class = int(background_class)
        if not 0 <= background_class < num_classes:
            raise ValueError(f"background_class must lie in [0, {num_classes - 1}]")

    # Everything below is small: the accumulator holds O(n_images) summaries.
    top_s, class_s, routed_s, conf = acc.stacked()
    nll_t = torch.tensor(acc.nll_sum, dtype=torch.float64)
    brier_t = torch.tensor(acc.brier_sum, dtype=torch.float64)
    nvox = torch.tensor(acc.n_valid, dtype=torch.float64)
    soft = bool(acc.soft)
    total_valid = float(nvox.sum())
    if total_valid == 0:
        raise ValueError("no valid voxels: every label is ignored or masked out")

    # --- global metrics and bootstrap over images -----------------------------------------
    pooled = top_s.sum(0)
    point = {
        "ece": float(_ece(pooled)),
        "nll": float(nll_t.sum() / nvox.sum()),
        "brier": float(brier_t.sum() / nvox.sum()),
        "gap": float(_gap(pooled)),
    }
    intervals = {k: NAN_CI for k in point}
    class_ece_ci: list[tuple[float, float]] = [NAN_CI] * num_classes
    routed_ece_ci: list[tuple[float, float]] = [NAN_CI] * num_classes
    boot_w: torch.Tensor | None = None
    if n_boot > 0 and n_images >= 2:
        g = torch.Generator().manual_seed(seed)
        draws = torch.randint(n_images, (n_boot, n_images), generator=g)
        w = torch.zeros(n_boot, n_images, dtype=torch.float64).scatter_add_(
            1, draws, torch.ones_like(draws, dtype=torch.float64)
        )
        boot_top = torch.einsum("rb,bkj->rkj", w, top_s)
        boot_n = w @ nvox
        intervals = {
            "ece": _ci(_ece(boot_top)),
            "nll": _ci((w @ nll_t) / boot_n),
            "brier": _ci((w @ brier_t) / boot_n),
            "gap": _ci(_gap(boot_top)),
        }
        boot_class = torch.einsum("rb,bckj->rckj", w, class_s)
        boot_routed = torch.einsum("rb,bckj->rckj", w, routed_s)
        class_ece_ci = [_ci(_ece(boot_class[:, c])) for c in range(num_classes)]
        routed_images = (routed_s[..., 0] > 0).sum(dim=(0, 2))  # (C,) images with routed voxels
        routed_ece_ci = [
            _ci(_ece(boot_routed[:, c])) if int(routed_images[c]) >= 2 else NAN_CI
            for c in range(num_classes)
        ]
        boot_w = w
    metrics = {k: (point[k], *intervals[k]) for k in point}

    # --- Dice, failures -------------------------------------------------------------------
    if soft:
        # Hard-label Dice is undefined for class distributions; do not invent a soft Dice.
        dice = torch.full((n_images, num_classes), float("nan"))
        mean_dice = torch.full((n_images,), float("nan"))
        pooled_dice = torch.full((num_classes,), float("nan"))
        label_totals = torch.stack(acc.label_counts).sum(0)  # (C,) expected mass
        pred_totals = routed_s[..., 0].sum(dim=(0, 2))  # (C,) argmax counts
    else:
        assert conf is not None
        tp = conf.diagonal(dim1=1, dim2=2)  # (B, C)
        lab, prd = conf.sum(2), conf.sum(1)
        dice = 2 * tp / (lab + prd)  # NaN where the class is in neither labels nor prediction
        present = lab > 0
        failed_cls = present & (dice < dice_failure_threshold)
        spurious = (lab == 0) & (prd > 0)
        considered = present.clone()
        if background_class is not None:
            considered[:, background_class] = False
        mean_dice = dice.masked_fill(~considered, float("nan")).nanmean(dim=1)  # (B,)
        failed_img = mean_dice < dice_failure_threshold
        pooled_conf = conf.sum(0)
        label_totals, pred_totals = pooled_conf.sum(1), pooled_conf.sum(0)
        pooled_dice = 2 * pooled_conf.diagonal() / (label_totals + pred_totals)

    dice_ci: list[tuple[float, float]] = [NAN_CI] * num_classes
    if boot_w is not None and not soft:
        for c in range(num_classes):
            if int(present[:, c].sum()) < 2:
                continue
            boot_dice = 2 * (boot_w @ tp[:, c]) / (boot_w @ (lab[:, c] + prd[:, c]))
            dice_ci[c] = _ci(boot_dice)

    routed_pool, class_pool = routed_s.sum(0), class_s.sum(0)
    classes = []
    for c in range(num_classes):
        class_ece = float(_ece(class_pool[c]))
        routed_ece = float(_ece(routed_pool[c]))
        if soft:
            label_vox: int | float = float(label_totals[c])
            pred_vox = int(pred_totals[c])
            dice_c = median = p10 = lowest = float("nan")
            n_present = n_failed = n_spurious = 0
        else:
            d_c = dice[:, c][present[:, c]]
            median, p10, lowest = _dice_summary(d_c)
            label_vox, pred_vox = int(label_totals[c]), int(pred_totals[c])
            dice_c = float(pooled_dice[c])
            n_present = int(present[:, c].sum())
            n_failed = int(failed_cls[:, c].sum())
            n_spurious = int(spurious[:, c].sum())
        classes.append(
            ClassStats(
                name=names[c],
                label_voxels=label_vox,
                label_fraction=label_vox / total_valid,
                pred_voxels=pred_vox,
                pred_fraction=pred_vox / total_valid,
                dice=dice_c,
                class_ece=class_ece,
                routed_ece=routed_ece,
                n_present=n_present,
                n_failed=n_failed,
                n_spurious=n_spurious,
                dice_median=median,
                dice_p10=p10,
                dice_min=lowest,
                class_ece_ci=class_ece_ci[c],
                routed_ece_ci=routed_ece_ci[c],
                dice_ci=dice_ci[c],
            )
        )

    if soft:
        failures = FailureStats(
            threshold=dice_failure_threshold,
            n_failed=0,
            ece_all=point["ece"],
            ece_without_failed=float("nan"),
            confidence_failed=float("nan"),
            confidence_rest=float("nan"),
        )
    else:

        def mean_conf(sel: torch.Tensor) -> float:
            s = top_s[sel].sum(0)
            return float(s[:, 1].sum() / s[:, 0].sum())

        rest = ~failed_img
        ranked = torch.argsort(mean_dice.nan_to_num(float("inf")))[: min(n_worst, n_images)]
        worst = []
        for i in (
            int(j) for j in ranked if float(mean_dice[j]) < 1.0
        ):  # NaN and perfect images excluded
            row = dice[i].masked_fill(~considered[i], float("inf"))
            c = int(row.argmin())
            worst.append((ids[i], float(mean_dice[i]), names[c], float(dice[i, c])))
        failures = FailureStats(
            threshold=dice_failure_threshold,
            n_failed=int(failed_img.sum()),
            ece_all=point["ece"],
            ece_without_failed=float(_ece(top_s[rest].sum(0))) if rest.any() else float("nan"),
            confidence_failed=mean_conf(failed_img) if failed_img.any() else float("nan"),
            confidence_rest=mean_conf(rest) if rest.any() else float("nan"),
            worst=worst,
        )

    # --- curve of the pooled top-1 statistics ---------------------------------------------
    counts = pooled[:, 0]
    bin_conf, bin_acc = pooled[:, 1] / counts.clamp_min(1), pooled[:, 2] / counts.clamp_min(1)
    nonempty = counts > 0
    curve = ReliabilityCurve(
        torch.linspace(0.0, 1.0, n_bins + 1).double(),
        bin_conf,
        bin_acc,
        counts.long(),
        point["ece"],
        float((bin_conf - bin_acc).abs()[nonempty].mean()),
    )

    report = ReliabilityReport(
        n_images=n_images,
        n_valid=int(total_valid),
        n_bins=n_bins,
        saturated_fraction=sum(acc.n_saturated) / total_valid,
        boundary=_boundary_stats(pooled, total_valid),
        metrics=metrics,
        classes=classes,
        failures=failures,
        curve=curve,
        image_dice=dice,
        warnings=[],
        soft_targets=soft,
    )
    report.warnings = _warnings(report, min_images, min_class_fraction, min_class_voxels)
    if background_class is not None and not soft and bool(mean_dice.isnan().all()):
        report.warnings.insert(
            0,
            "background_class leaves no class in the per-image Dice; "
            "image failure statistics are undefined",
        )
    return report

ReliabilityReport dataclass

ReliabilityReport(n_images: int, n_valid: int, n_bins: int, saturated_fraction: float, boundary: BoundaryStats, metrics: dict[str, tuple[float, float, float]], classes: list[ClassStats], failures: FailureStats, curve: ReliabilityCurve, image_dice: Tensor, warnings: list[str], soft_targets: bool = False)

Result of :func:diagnose. print(report) shows :meth:summary.

metrics maps ece, nll, brier and gap (mean confidence minus accuracy) to (value, ci_low, ci_high); the interval is NaN when it could not be computed. image_dice is the (n_images, n_classes) Dice matrix, NaN where the class is absent from both labels and prediction, or everywhere for soft targets (hard-label Dice is undefined there). boundary reports the first and last confidence bins. soft_targets records which target kind was diagnosed.

bias property

bias: str

overconfident/underconfident if the interval of the gap excludes 0.

summary

summary() -> str

Plain-text report for a terminal or a notebook cell.

Source code in src/fiducio/diagnostics.py

def summary(self) -> str:
    """Plain-text report for a terminal or a notebook cell."""
    m = self.metrics

    def est(key: str) -> str:
        return _fmt_interval(*m[key])

    lines = [
        f"Reliability report: {self.n_images} images, {self.n_valid} valid voxels, "
        f"{len(self.classes)} classes",
        f"  ECE   {est('ece')}",
        f"  NLL   {est('nll')}",
        f"  Brier {est('brier')}",
        f"  Bias  {self.bias} (confidence - accuracy = {est('gap')})",
        f"  Saturated voxels (confidence > {_SATURATION}): {self.saturated_fraction:.1%}",
        f"  First/last bin mass {self.boundary.low_bin_mass:.1%}/"
        f"{self.boundary.high_bin_mass:.1%}, gap "
        f"{_f(self.boundary.low_bin_gap)}/{_f(self.boundary.high_bin_gap)}, "
        f"ECE contribution {_f(self.boundary.ece_contribution)}",
        "",
    ]
    rows = [["class", "label %", "pred %", "Dice", "class ECE", "routed vox", "routed ECE"]] + [
        [
            c.name,
            f"{100 * c.label_fraction:.2f}",
            f"{100 * c.pred_fraction:.2f}",
            _f(c.dice),
            _f(c.class_ece),
            str(c.pred_voxels),
            _f(c.routed_ece),
        ]
        for c in self.classes
    ]
    lines += _table(rows) + [""]
    lines.append("  Per-class 95% intervals (bootstrap over images):")
    rows = [["class", "class ECE", "routed ECE", "Dice"]] + [
        [
            c.name,
            _fmt_interval(c.class_ece, *c.class_ece_ci),
            _fmt_interval(c.routed_ece, *c.routed_ece_ci),
            _fmt_interval(c.dice, *c.dice_ci),
        ]
        for c in self.classes
    ]
    lines += _table(rows) + [""]
    f = self.failures
    lines.append(f"Failed images (mean Dice < {f.threshold:g}): {f.n_failed}/{self.n_images}")
    lines.append(
        f"  ECE with all images {_f(f.ece_all)}, without the failed ones "
        f"{_f(f.ece_without_failed)}; mean confidence {_f(f.confidence_failed)} "
        f"on failed vs {_f(f.confidence_rest)} on the rest"
    )
    lines.append(
        f"  Per class, a failure is Dice < {f.threshold:g} on an image where it is present:"
    )
    rows = [["class", "present", "failed", "spurious", "Dice median", "p10", "min"]] + [
        [
            c.name,
            str(c.n_present),
            str(c.n_failed),
            str(c.n_spurious),
            _f(c.dice_median),
            _f(c.dice_p10),
            _f(c.dice_min),
        ]
        for c in self.classes
    ]
    lines += _table(rows)
    if f.worst:
        lines.append("  Lowest-Dice images:")
        lines += [
            f"    {name}: mean Dice {_f(d)}, worst class {cls} ({_f(cd)})"
            for name, d, cls, cd in f.worst
        ]
    lines += ["", "Warnings" if self.warnings else "No warnings"]
    lines += [f"  - {w}" for w in self.warnings]
    return "\n".join(lines)

to_dict

to_dict() -> dict[str, Any]

JSON-serializable dictionary (NaN becomes None).

Source code in src/fiducio/diagnostics.py

def to_dict(self) -> dict[str, Any]:
    """JSON-serializable dictionary (NaN becomes ``None``)."""
    out: dict[str, Any] = _clean(dataclasses.asdict(self) | {"bias": self.bias})
    return out

plot

plot() -> Any

Reliability, confidence histogram, per-class ECE and per-image Dice (2x2).

Source code in src/fiducio/diagnostics.py

def plot(self) -> Any:
    """Reliability, confidence histogram, per-class ECE and per-image Dice (2x2)."""
    plt = _require_matplotlib()
    fig, axes = plt.subplots(2, 2, figsize=(10, 9))
    _draw_reliability(
        axes[0, 0], self.curve, class_index=None, title="Reliability (top-1)", metric="ece"
    )

    edges = self.curve.bin_edges.cpu().numpy()
    counts = self.curve.bin_counts.cpu().numpy()
    ax = axes[0, 1]
    ax.bar(
        (edges[:-1] + edges[1:]) / 2, counts / max(counts.sum(), 1), width=edges[1] - edges[0]
    )
    ax.set_yscale("log")
    ax.set_xlabel("confidence")
    ax.set_ylabel("share of voxels")
    ax.set_title("Confidence histogram")

    names = [f"{c.name}\n{100 * c.label_fraction:.1f}%" for c in self.classes]
    ax = axes[1, 0]
    ax.bar(
        range(len(names)),
        [0.0 if math.isnan(c.class_ece) else c.class_ece for c in self.classes],
    )
    ax.set_xticks(range(len(names)), names)
    ax.set_ylabel("class ECE")
    ax.set_title("Per-class ECE (share of labels below the name)")

    ax = axes[1, 1]
    dice = self.image_dice.cpu().numpy()
    # Deterministic jitter so that points of one class do not overlap.
    jitter = ((torch.arange(dice.shape[0]) * 0.618) % 1 - 0.5).numpy() * 0.4
    for c in range(dice.shape[1]):
        ax.scatter(c + jitter, dice[:, c], s=8, alpha=0.6)
    ax.axhline(self.failures.threshold, color="red", linestyle="--", label="failure threshold")
    ax.set_xticks(range(len(names)), [c.name for c in self.classes])
    ax.set_ylim(-0.02, 1.02)
    ax.set_ylabel("Dice per image")
    ax.set_title("Per-image Dice")
    ax.legend(loc="lower right", fontsize="small")
    fig.tight_layout()
    return fig

Helpers

mean_human_response

mean_human_response(annotations: ArrayLike, *, num_classes: int, ignore_index: int = DEFAULT_IGNORE_INDEX) -> tuple[torch.Tensor, torch.Tensor]

Return empirical class proportions and an annotated-voxel mask.

annotations has shape (B, R, *spatial) with rater axis 1. Missing votes equal ignore_index. Unannotated voxels have zero mass; pass the returned mask to fit and metrics. Each valid voxel receives equal weight, regardless of how many available raters it has.

Source code in src/fiducio/utils/tensors.py

def mean_human_response(
    annotations: ArrayLike, *, num_classes: int, ignore_index: int = DEFAULT_IGNORE_INDEX
) -> tuple[torch.Tensor, torch.Tensor]:
    """Return empirical class proportions and an annotated-voxel mask.

    ``annotations`` has shape ``(B, R, *spatial)`` with rater axis 1.
    Missing votes equal ``ignore_index``. Unannotated voxels have zero mass;
    pass the returned mask to ``fit`` and metrics. Each valid voxel receives
    equal weight, regardless of how many available raters it has.
    """
    try:
        classes = -1 if isinstance(num_classes, bool) else operator.index(num_classes)
    except TypeError:
        classes = -1
    if classes < 2:
        raise ValueError("num_classes must be an integer >= 2")
    a = to_tensor(annotations).detach()
    if a.ndim < 2:
        raise ValueError("annotations must have shape (B, R, *spatial)")
    a = integer_targets(a, device=a.device)
    present = a != ignore_index
    if bool(((a < 0) | (a >= classes))[present].any()):
        raise ValueError(f"annotation labels must lie in [0, {classes - 1}]")
    counts = present.sum(dim=1)
    proportions = torch.stack(
        [((a == c) & present).sum(dim=1) for c in range(classes)], dim=1
    ).float()
    proportions /= counts.clamp_min(1).unsqueeze(1)
    return proportions, counts > 0

two_channel_from_binary

two_channel_from_binary(scores: ArrayLike, *, input_type: str = 'logits') -> torch.Tensor

Convert single-channel binary outputs to the two-channel form Fiducio uses.

Fiducio represents binary segmentation as two channels (C = 2). Use this helper to convert a single-channel sigmoid output, whose class axis at dimension 1 has size 1 (shape (B, 1, *spatial)), into (B, 2, *spatial). Class 1 is the positive class.

Parameters:

Name	Type	Description	Default
scores	ArrayLike	(B, 1, *spatial) score for the positive class.	required
input_type	str	"logits" — a sigmoid logit z; returns logits [0, z] whose softmax equals [1 - sigmoid(z), sigmoid(z)]. "probs" — a probability p; returns [1 - p, p].	'logits'

Source code in src/fiducio/utils/tensors.py

def two_channel_from_binary(scores: ArrayLike, *, input_type: str = "logits") -> torch.Tensor:
    """Convert single-channel binary outputs to the two-channel form Fiducio uses.

    Fiducio represents binary segmentation as two channels (``C = 2``). Use this
    helper to convert a single-channel sigmoid output, whose class axis at
    dimension 1 has size 1 (shape ``(B, 1, *spatial)``), into ``(B, 2, *spatial)``.
    Class 1 is the positive class.

    Parameters
    ----------
    scores:
        ``(B, 1, *spatial)`` score for the positive class.
    input_type:
        ``"logits"`` — a sigmoid logit ``z``; returns logits ``[0, z]`` whose
        softmax equals ``[1 - sigmoid(z), sigmoid(z)]``. ``"probs"`` — a
        probability ``p``; returns ``[1 - p, p]``.
    """
    x = to_tensor(scores, dtype=torch.float32)
    if x.ndim < 2 or x.shape[1] != 1:
        raise ValueError(
            "scores must have a singleton class axis at dimension 1, i.e. shape "
            f"(B, 1, *spatial); got {tuple(x.shape)}"
        )
    if input_type == "logits":
        return torch.cat([torch.zeros_like(x), x], dim=1)
    if input_type == "probs":
        return torch.cat([1.0 - x, x], dim=1)
    raise ValueError(f"input_type must be 'logits' or 'probs', got {input_type!r}")