Skip to content

stemflow.calibration

Composable calibration maps, OOF estimators and spatial-temporal diagnostics. See local calibration and tuning.

Public choices

Component Purpose
BetaCalibrator Fit an increasing binary probability map with an identity-centered ridge penalty.
AffineCalibrator Calibrate regression point means.
CrossValidatedCalibratedClassifier / Regressor Fit the map from out-of-fold scores and retain fold estimators.
SpatioTemporalKFold Block by space/time and audit fold support, retries and fallback status.
calibration_scores, calibration_bins Summarize held-out scores and bin counts/uncertainty references.
calibration_by_group, stixel_calibration_diagnostics Diagnose supplied groups or actual fitted planar stixel geometry.
plot_calibration_diagnostics Save multi-level diagnostics and stixel overlays from those tables.

For SpatioTemporalKFold, n_splits controls folds independently of spatial/temporal block counts. min_pos/min_neg are feasibility thresholds for each training and validation fold, not statistical reliability guarantees. rare_policy="allow" reports a fallback; "raise" rejects insufficient support. Inspect fold_support_, support_sufficient_, status_ and trials_.

Generic wrappers use shuffled folds when cv is omitted. Supply the spatial-temporal splitter explicitly for blocked calibration. Intentional all-zero/all-one AdaSTEM stixels keep their dummy behavior. Local calibration and relative ECE do not establish exact conditional probabilities or copula accuracy; examine support, coverage, proper losses and floor sensitivity together.

Standalone probability and point-prediction calibration maps.

AffineCalibrator

Bases: BaseEstimator

Ridge correction of a regression point mean; not interval calibration.

Source code in stemflow/calibration/maps.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
class AffineCalibrator(BaseEstimator):
    """Ridge correction of a regression point mean; not interval calibration."""

    def __init__(self, ridge=0.):
        self.ridge = ridge

    def fit(self, p, y, sample_weight=None):
        p, y = np.asarray(p, float).reshape(-1), np.asarray(y, float).reshape(-1)
        if not len(y) or len(p) != len(y) or not np.isfinite(p).all() or not np.isfinite(y).all():
            raise ValueError("Expected matching finite regression predictions and targets.")
        if not np.isfinite(self.ridge) or self.ridge < 0:
            raise ValueError("ridge must be finite and nonnegative.")
        x, w = np.column_stack([p, np.ones(len(p))]), _weights(sample_weight, len(p))
        center = np.array([1., 0.])
        self.theta_ = np.linalg.lstsq(x.T @ (w[:, None] * x) + self.ridge*np.eye(2),
            x.T @ (w*y) + self.ridge*center, rcond=None)[0]
        return self

    def predict(self, p):
        check_is_fitted(self, "theta_")
        return self.theta_[0] * np.asarray(p, float).reshape(-1) + self.theta_[1]

BetaCalibrator

Bases: BaseEstimator

Increasing beta map, regularized toward the identity probability map.

Minimizes SUM weighted binary log loss plus ridge / 2 * ((a - 1)**2 + (b - 1)**2 + c**2) for sigmoid(a*log(p) - b*log(1-p) + c). Ridge strength therefore depends on sample size and weight scale. The default matches the supplied pilot.

Source code in stemflow/calibration/maps.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
class BetaCalibrator(BaseEstimator):
    """Increasing beta map, regularized toward the identity probability map.

    Minimizes SUM weighted binary log loss plus
    ``ridge / 2 * ((a - 1)**2 + (b - 1)**2 + c**2)`` for
    ``sigmoid(a*log(p) - b*log(1-p) + c)``. Ridge strength therefore depends
    on sample size and weight scale. The default matches the supplied pilot.
    """

    def __init__(self, ridge=100., eps=1e-10):
        self.ridge, self.eps = ridge, eps

    def _design(self, p):
        p = np.clip(np.asarray(p, float).reshape(-1), self.eps, 1 - self.eps)
        return np.column_stack([np.log(p), -np.log1p(-p), np.ones(len(p))])

    def fit(self, p, y, sample_weight=None):
        p, y = np.asarray(p, float).reshape(-1), np.asarray(y, float).reshape(-1)
        if not len(y) or len(p) != len(y) or not np.isfinite(p).all():
            raise ValueError("Expected matching nonempty finite probabilities and labels.")
        if not np.isin(y, [0, 1]).all() or np.any((p < 0) | (p > 1)):
            raise ValueError("Expected binary labels and probabilities in [0,1].")
        if not np.isfinite(self.ridge) or self.ridge < 0 or not 0 < self.eps < .5:
            raise ValueError("ridge must be nonnegative and eps between 0 and 0.5.")
        w = _weights(sample_weight, len(y))
        x, center = self._design(p), np.array([1., 1., 0.])
        def objective(theta):
            z, delta = x @ theta, theta - center
            value = np.sum(w * (np.logaddexp(0, z) - y*z)) + self.ridge/2 * (delta @ delta)
            gradient = x.T @ (w * (expit(z) - y)) + self.ridge * delta
            return value, gradient
        bounds = [(0.02, 20.), (0.02, 20.), (-30., 30.)]
        result = minimize(objective, center, jac=True, method="L-BFGS-B", bounds=bounds,
            options={"maxiter": 250, "ftol": 1e-11, "gtol": 1e-7, "maxls": 50})
        if not result.success:
            result = minimize(objective, result.x, jac=True, method="SLSQP", bounds=bounds,
                              options={"maxiter": 300, "ftol": 1e-10})
            gradient = objective(result.x)[1].copy()
            for j, (lo, hi) in enumerate(bounds):
                if result.x[j] <= lo + 1e-7:
                    gradient[j] = min(gradient[j], 0.)
                if result.x[j] >= hi - 1e-7:
                    gradient[j] = max(gradient[j], 0.)
            if not result.success and np.max(np.abs(gradient)) > 1e-5:
                raise RuntimeError(f"Beta optimizer failed: {result.message}")
        if not np.isfinite(result.fun) or result.fun > objective(center)[0] + 1e-6:
            raise RuntimeError("Beta optimizer failed its objective check.")
        self.theta_, self.n_samples_ = result.x, len(y)
        self.penalty_ = float(self.ridge/2 * np.sum((self.theta_ - center)**2))
        self.objective_ = float(result.fun)
        return self

    def predict(self, p):
        check_is_fitted(self, "theta_")
        return expit(self._design(p) @ self.theta_)

Composable OOF calibration estimators, usable inside or outside AdaSTEM.

ConstantBinaryClassifier

Bases: ClassifierMixin, BaseEstimator

Smoothed fold fallback, distinct from AdaSTEM's intentional dummy zeros.

Source code in stemflow/calibration/estimators.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
class ConstantBinaryClassifier(ClassifierMixin, BaseEstimator):
    """Smoothed fold fallback, distinct from AdaSTEM's intentional dummy zeros."""

    def fit(self, X, y, sample_weight=None):
        y = np.asarray(y).reshape(-1)
        w = _weights(sample_weight, len(y))
        self.probability_ = float((w @ y + .5) / (w.sum() + 1))
        self.classes_ = np.array([0, 1])
        self.feature_importances_ = np.zeros(X.shape[1])
        return self

    def predict_proba(self, X):
        return np.tile([1 - self.probability_, self.probability_], (len(X), 1))

    def predict(self, X):
        return np.full(len(X), int(self.probability_ >= .5))

CrossValidatedCalibratedClassifier

Bases: ClassifierMixin, _CrossValidatedCalibration, BaseEstimator

Fit retained fold models and one map on pooled OOF probabilities.

Each fold probability is calibrated before averaging, matching the supplied local-beta strategy. Inside AdaSTEM this is a separate map per stixel. Calibrators follow fit(probability, y) / predict(probability). Use an explicit SpatioTemporalKFold for within-stixel blocked calibration. A fixed estimator requires no tuning; wrap it in SelfTuningClassifier only when nested hyperparameter selection is wanted.

Source code in stemflow/calibration/estimators.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
class CrossValidatedCalibratedClassifier(ClassifierMixin, _CrossValidatedCalibration, BaseEstimator):
    """Fit retained fold models and one map on pooled OOF probabilities.

    Each fold probability is calibrated before averaging, matching the supplied
    local-beta strategy. Inside AdaSTEM this is a separate map per stixel.
    Calibrators follow ``fit(probability, y)`` / ``predict(probability)``.
    Use an explicit SpatioTemporalKFold for within-stixel blocked calibration.
    A fixed estimator requires no tuning; wrap it in SelfTuningClassifier only
    when nested hyperparameter selection is wanted.
    """
    _classification = True

    def __init__(self, estimator, calibrator=None, cv=None, exclude_features=None):
        self.estimator, self.calibrator, self.cv = estimator, calibrator, cv
        self.exclude_features = exclude_features

    def predict_proba(self, X):
        p = self._predict_mean(X)
        return np.column_stack([1 - p, p])

    def predict(self, X):
        return (self._predict_mean(X) >= .5).astype(int)

CrossValidatedCalibratedRegressor

Bases: RegressorMixin, _CrossValidatedCalibration, BaseEstimator

OOF calibration of point means; supplies no prediction-interval guarantee.

Source code in stemflow/calibration/estimators.py 
163
164
165
166
167
168
169
170
171
172
class CrossValidatedCalibratedRegressor(RegressorMixin, _CrossValidatedCalibration, BaseEstimator):
    """OOF calibration of point means; supplies no prediction-interval guarantee."""
    _classification = False

    def __init__(self, estimator, calibrator=None, cv=None, exclude_features=None):
        self.estimator, self.calibrator, self.cv = estimator, calibrator, cv
        self.exclude_features = exclude_features

    def predict(self, X):
        return self._predict_mean(X)

Audited spatial-temporal folds for local calibration and nested tuning.

InsufficientSplitError

Bases: ValueError

There are no nonempty independent spatial-temporal folds.

Source code in stemflow/calibration/cross_validation.py 
10
11
class InsufficientSplitError(ValueError):
    """There are no nonempty independent spatial-temporal folds."""

SpatioTemporalKFold

Bases: BaseEstimator

Retry blocked folds, recording event support rather than hiding scarcity.

rare_policy='allow' retains the first nonempty blocked split when the requested class support cannot be met. No random row split is substituted. period folds circular-year copies onto their original seasonal day. Spatial and temporal block counts are configurable, independently of n_splits. They retain ST_KFold's linspace-cutpoint semantics and use the observed range of the data passed to each split, without a buffer.

Source code in stemflow/calibration/cross_validation.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
class SpatioTemporalKFold(BaseEstimator):
    """Retry blocked folds, recording event support rather than hiding scarcity.

    ``rare_policy='allow'`` retains the first nonempty blocked split when the
    requested class support cannot be met. No random row split is substituted.
    ``period`` folds circular-year copies onto their original seasonal day.
    Spatial and temporal block counts are configurable, independently of
    ``n_splits``. They retain ST_KFold's linspace-cutpoint semantics and use
    the observed range of the data passed to each split, without a buffer.
    """

    def __init__(self, n_splits=3, Spatio1="longitude", Spatio2="latitude",
                 Temporal1="DOY", Spatio_blocks_count=5, Temporal_blocks_count=5,
                 random_state=42, max_trials=10, min_pos=1, min_neg=1,
                 rare_policy="allow", period=366):
        self.n_splits = n_splits
        self.Spatio1, self.Spatio2, self.Temporal1 = Spatio1, Spatio2, Temporal1
        self.Spatio_blocks_count = Spatio_blocks_count
        self.Temporal_blocks_count = Temporal_blocks_count
        self.random_state, self.max_trials = random_state, max_trials
        self.min_pos, self.min_neg = min_pos, min_neg
        self.rare_policy, self.period = rare_policy, period

    @property
    def metadata_columns(self):
        return (self.Spatio1, self.Spatio2, self.Temporal1)

    @property
    def feature_exclusions(self):
        return (self.Spatio1, self.Spatio2)

    def get_n_splits(self, X=None, y=None, groups=None):
        return self.n_splits

    def split(self, X, y=None, groups=None):
        if not isinstance(X, pd.DataFrame):
            raise TypeError("Spatial-temporal cross-validation requires a DataFrame.")
        if self.n_splits < 2 or self.max_trials < 1 or self.min_pos < 0 or self.min_neg < 0:
            raise ValueError("Invalid fold count, trial count or event threshold.")
        if self.rare_policy not in {"allow", "raise"}:
            raise ValueError("rare_policy must be 'allow' or 'raise'.")
        geometry = X[list(self.metadata_columns)].copy()
        if self.period is not None:
            if self.period <= 0:
                raise ValueError("period must be positive or None.")
            geometry[self.Temporal1] = (geometry[self.Temporal1] - 1) % self.period + 1
        if not np.isfinite(geometry.to_numpy()).all():
            raise ValueError("Spatial-temporal coordinates must be finite.")
        labels = None if y is None else np.asarray(y).reshape(-1)
        if labels is not None and len(labels) != len(X):
            raise ValueError("X and y lengths differ.")
        binary = labels is not None and np.isin(labels, [0, 1]).all()
        first, selected, sufficient = None, None, False
        for trial in range(self.max_trials):
            seed = None if self.random_state is None else self.random_state + trial
            splitter = ST_KFold(Spatio1=self.Spatio1, Spatio2=self.Spatio2,
                Temporal1=self.Temporal1, Spatio_blocks_count=self.Spatio_blocks_count,
                Temporal_blocks_count=self.Temporal_blocks_count,
                n_splits=self.n_splits, random_state=seed)
            folds = [(np.asarray(a, int), np.asarray(b, int)) for a, b in splitter.split(geometry)]
            if not all(len(a) and len(b) for a, b in folds):
                continue
            if first is None:
                first = folds
            sufficient = not binary or all(
                np.sum(labels[v] == 1) >= self.min_pos and
                np.sum(labels[v] == 0) >= self.min_neg for fold in folds for v in fold)
            if sufficient:
                selected = folds
                break
        self.support_sufficient_ = sufficient
        self.trials_ = trial + 1
        selected = selected if selected is not None else first
        if selected is None or (not sufficient and self.rare_policy == "raise"):
            raise InsufficientSplitError("No blocked split meets the requested support.")
        self.status_ = "supported" if sufficient else "rare_fold_fallback"
        self.fold_support_ = [dict(n_train=len(a), n_validation=len(b),
            training_events=int(np.sum(labels[a] == 1)) if binary else None,
            validation_events=int(np.sum(labels[b] == 1)) if binary else None)
            for a, b in selected]
        yield from selected

Binary calibration metrics with explicit rare-bin support diagnostics.

Held-out diagnostics for final probabilities inside actual fitted stixels.

calibration_by_group(y, p, groups, n_bins=8, include_bins=True, common_bin_edges=None)

Score overlapping positional groups; units are not assumed independent.

Source code in stemflow/calibration/diagnostics.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def calibration_by_group(y, p, groups, n_bins=8, include_bins=True, common_bin_edges=None):
    """Score overlapping positional groups; units are not assumed independent."""
    y, p = np.asarray(y).reshape(-1), np.asarray(p).reshape(-1)
    if len(y) != len(p):
        raise ValueError("Prediction/outcome lengths differ.")
    scores, tables = [], []
    for metadata, positions in groups:
        values, bins = calibration_scores(y[positions], p[positions], n_bins)
        scores.append({**metadata, **values})
        if include_bins and len(bins):
            tables.append(bins.assign(**metadata))
            if common_bin_edges is not None:
                good = np.isfinite(y[positions]) & np.isfinite(p[positions])
                common = calibration_bins(y[positions][good], p[positions][good],
                    bin_edges=common_bin_edges).assign(strategy="common", **metadata)
                tables.append(common)
    return pd.DataFrame(scores), pd.concat(tables, ignore_index=True) if tables else pd.DataFrame()

plot_calibration_diagnostics(diagnostics, output_path, max_curves=200)

Save overall/stixel diagnostics and transparent log-scale overlays.

Source code in stemflow/calibration/diagnostics.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
def plot_calibration_diagnostics(diagnostics, output_path, max_curves=200):
    """Save overall/stixel diagnostics and transparent log-scale overlays."""
    from pathlib import Path
    import matplotlib.pyplot as plt
    from matplotlib.ticker import FuncFormatter
    output = Path(output_path)
    output.mkdir(parents=True, exist_ok=True)
    scores, bins = diagnostics["scores"], diagnostics["bins"]
    if bins.empty:
        raise ValueError("Plotting requires diagnostics computed with include_bins=True.")
    stixels = scores[(scores.level == "stixel") & (scores.n_valid > 0)]
    fig, axes = plt.subplots(2, 3, figsize=(15, 8))
    for strategy in ["uniform", "quantile"]:
        b = bins[(bins.level == "overall") & (bins.strategy == strategy)]
        axes[0, 0].plot(b.predicted, b.observed, "o-", label=strategy)
    axes[0, 0].plot([0, 1], [0, 1], "k--")
    axes[0, 0].set(xlabel="Predicted probability", ylabel="Observed frequency", title="Overall reliability")
    axes[0, 0].legend()
    axes[0, 1].scatter(stixels.predicted_rate, stixels.observed_rate, s=5, alpha=.25, rasterized=True)
    axes[0, 1].plot([0, 1], [0, 1], "k--")
    axes[0, 1].set(xlabel="Stixel mean prediction", ylabel="Stixel observed frequency", title="Actual fitted stixels")
    for axis, metric in zip([axes[0, 2], axes[1, 0]], ["quantile_ECE_pp", "quantile_ECE_relative"]):
        for label, group in [("All", stixels), ("Training-supported", stixels[stixels.training_supported])]:
            v = np.sort(group[metric].dropna())
            if len(v):
                axis.plot(v, np.arange(1, len(v)+1)/len(v), label=f"{label}: {len(v)}")
        axis.set(xlabel=metric, ylabel="Fraction of stixels", ylim=(0, 1))
        axis.legend(fontsize=8)
        if metric.endswith("relative"):
            axis.set_xscale("log")
    spatial = scores[scores.level == "spatial"]
    half = max(float(spatial.bias_pp.abs().max()), 1e-9)
    artist = axes[1, 1].scatter(spatial.longitude, spatial.latitude, c=spatial.bias_pp,
        cmap="coolwarm", vmin=-half, vmax=half, s=20)
    axes[1, 1].set(xlabel="Longitude", ylabel="Latitude", title="Spatial-bin bias")
    fig.colorbar(artist, ax=axes[1, 1], label="Predicted − observed (pp)")
    temporal = scores[scores.level == "temporal"].sort_values("DOY")
    axes[1, 2].plot(temporal.DOY, temporal.observed_rate, "o-", label="Observed")
    axes[1, 2].plot(temporal.DOY, temporal.predicted_rate, "o-", label="Predicted")
    axes[1, 2].set(xlabel="DOY", ylabel="Occurrence probability", title="Temporal-bin means")
    axes[1, 2].legend()
    fig.tight_layout()
    for ext in ["png", "pdf"]:
        fig.savefig(output/f"calibration_diagnostics.{ext}", dpi=180)
    plt.close(fig)
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    names = stixels.group.to_numpy()
    if len(names) > max_curves:
        names = np.random.default_rng(42).choice(names, max_curves, replace=False)
    for axis, strategy in zip(axes, ["uniform", "quantile"]):
        b = bins[(bins.level == "stixel") & (bins.strategy == strategy)]
        for _, curve in b[b.group.isin(names)].groupby("group", sort=False):
            good = (curve.predicted > 0) & (curve.observed > 0)
            axis.plot(curve.predicted.where(good), curve.observed.where(good), color="C0", alpha=.08)
        # Common probability intervals avoid averaging unrelated quantile-bin numbers.
        common = bins[(bins.level == "stixel") & (bins.strategy == "common")]
        mean = common.groupby("bin")[["predicted", "observed"]].mean()
        good = (mean.predicted > 0) & (mean.observed > 0)
        axis.plot(mean.predicted.where(good), mean.observed.where(good), color="crimson", lw=2,
                  label="Equal-stixel mean in common probability bins")
        axis.plot([1e-8, 1], [1e-8, 1], "k--", label="Perfect calibration")
        axis.set(xscale="log", yscale="log", xlim=(1e-8, 1), ylim=(1e-8, 1),
                 xlabel="Mean predicted probability", ylabel="Observed frequency", title=strategy)
        formatter = FuncFormatter(lambda value, _: f"{value:.8f}".rstrip("0").rstrip("."))
        axis.xaxis.set_major_formatter(formatter)
        axis.yaxis.set_major_formatter(formatter)
        axis.tick_params(axis="x", rotation=45)
        axis.legend(fontsize=8)
    fig.suptitle("Held-out final probabilities within stixels; zeros break log lines")
    fig.tight_layout()
    for ext in ["png", "pdf"]:
        fig.savefig(output/f"stixel_calibration_overlay.{ext}", dpi=180)
    plt.close(fig)

stixel_calibration_diagnostics(model, X, y, pred=None, n_bins=8, spatial_bin_width=5.0, temporal_bin_width=40.0, include_bins=True, support_min_events=20, support_min_samples=200)

Evaluate on caller-supplied held-out rows, using the fitted partition.

Returns overall, spatial, temporal, space-time and overlapping stixel scores. Group accuracy uses final ensemble probabilities, not in-sample fold scores. Relative ECE is floor-sensitive and accompanied by zero-event bin mass.

Source code in stemflow/calibration/diagnostics.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def stixel_calibration_diagnostics(model, X, y, pred=None, n_bins=8,
                                   spatial_bin_width=5., temporal_bin_width=40.,
                                   include_bins=True, support_min_events=20,
                                   support_min_samples=200):
    """Evaluate on caller-supplied held-out rows, using the fitted partition.

    Returns overall, spatial, temporal, space-time and overlapping stixel scores.
    Group accuracy uses final ensemble probabilities, not in-sample fold scores.
    Relative ECE is floor-sensitive and accompanied by zero-event bin mass.
    """
    from ..lazyloading.prediction import SpatialWindowIndex
    from ..model.static_func_AdaSTEM import transform_pred_set_to_STEM_quad
    if not isinstance(X, pd.DataFrame) or not X.index.is_unique:
        raise TypeError("Diagnostics require a DataFrame with unique observation IDs.")
    if isinstance(y, (pd.Series, pd.DataFrame)) and not y.index.equals(X.index):
        raise ValueError("Held-out X and y indices must match in order.")
    if spatial_bin_width <= 0 or temporal_bin_width <= 0:
        raise ValueError("Diagnostic bin widths must be positive.")
    if support_min_events < 0 or support_min_samples < 1:
        raise ValueError("Diagnostic support thresholds must be nonnegative/positive.")
    if "stixel_calibration_point_transformed_left_bound" not in model.ensemble_df:
        raise NotImplementedError("Actual-stixel diagnostics currently require planar AdaSTEM geometry.")
    coordinates = X[[model.Spatio1, model.Spatio2, model.Temporal1]]
    if not np.isfinite(coordinates.to_numpy()).all():
        raise ValueError("Diagnostic coordinates must be finite.")
    y = np.asarray(y).reshape(-1)
    if pred is None:
        pred = model.predict_proba(X, temporal_window_prequery=True)[:, 1]
    pred = np.asarray(pred).reshape(-1)
    if len(X) != len(y) or len(X) != len(pred):
        raise ValueError("Held-out X, y and predictions must have matching lengths.")
    spatial = list(zip(np.floor(X[model.Spatio1].to_numpy()/spatial_bin_width).astype(int),
                       np.floor(X[model.Spatio2].to_numpy()/spatial_bin_width).astype(int)))
    temporal = np.floor((X[model.Temporal1].to_numpy()-1)/temporal_bin_width).astype(int)
    def location(positions):
        return dict(longitude=X[model.Spatio1].iloc[positions].mean(),
                    latitude=X[model.Spatio2].iloc[positions].mean(),
                    DOY=X[model.Temporal1].iloc[positions].mean())
    def groups():
        yield dict(level="overall", group="all"), np.arange(len(X))
        keys = pd.DataFrame(dict(spatial=pd.Series(spatial), temporal=temporal))
        for level, cols in [("spatial", ["spatial"]), ("temporal", ["temporal"]),
                            ("spatiotemporal", ["spatial", "temporal"])]:
            for name, positions in keys.groupby(cols, sort=True).indices.items():
                yield dict(level=level, group=str(name), **location(positions)), positions
        columns = [model.Temporal1, model.Spatio1, model.Spatio2]
        for ensemble, grid in model.ensemble_df.groupby("ensemble_index", sort=True):
            for start, window_grid in grid.groupby(f"{model.Temporal1}_start", sort=True):
                positions = np.flatnonzero((X[model.Temporal1].to_numpy() >= start) &
                    (X[model.Temporal1].to_numpy() < start + model.temporal_bin_interval))
                coords = transform_pred_set_to_STEM_quad(model.Spatio1, model.Spatio2,
                    X.iloc[positions][columns], grid)
                lookup = SpatialWindowIndex(coords[f"{model.Spatio1}_new"].to_numpy(),
                                            coords[f"{model.Spatio2}_new"].to_numpy(), len(window_grid))
                for _, row in window_grid.iterrows():
                    name = str(row.unique_stixel_id)
                    if f"{name}_model" not in model.model_dict:
                        continue
                    estimator = model.model_dict[f"{name}_model"]
                    support = vars(estimator).get("_stemflow_training_support_", {})
                    selected = lookup.select(row)
                    yield dict(level="stixel", group=name, ensemble_index=ensemble,
                        DOY_start=start, partition_n=row.stixel_checklist_count,
                        training_n=support.get("n", np.nan),
                        training_detections=support.get("positives", np.nan),
                        **location(positions[selected])), positions[selected]
    scores, bins = calibration_by_group(y, pred, groups(), n_bins, include_bins,
                                      common_bin_edges=np.r_[0., np.geomspace(1e-8, 1, 65)])
    for name in ["training_n", "training_detections"]:
        if name not in scores:
            scores[name] = np.nan
    scores["supported"] = ((scores.n_valid >= support_min_samples) &
        (scores.detections >= support_min_events) & (scores.non_detections >= support_min_events))
    scores["training_supported"] = ((scores.level == "stixel") &
        (scores.training_detections >= support_min_events) &
        (scores.training_n-scores.training_detections >= support_min_events))
    metrics = ["ECE_pp", "quantile_ECE_pp", "ECE_relative", "quantile_ECE_relative",
               "quantile_relative_floor_1e8", "quantile_zero_event_mass", "bias_pp",
               "rate_ratio", "AUC", "average_precision", "Brier", "log_loss", "coverage"]
    records = []
    for scope, selected in [("all", scores), ("supported", scores[scores.supported]),
                            ("training_supported", scores[scores.training_supported])]:
        for level, rows in selected.groupby("level", sort=True):
            for metric in metrics:
                values = rows[metric].replace([np.inf, -np.inf], np.nan).dropna()
                records.append(dict(scope=scope, level=level, metric=metric, n_groups=len(values),
                    **dict(zip(["q25", "median", "q75", "p90"], values.quantile([.25, .5, .75, .9])))))
    return dict(scores=scores, bins=bins, summary=pd.DataFrame(records))