Local calibration and streamed prediction¶

This 2.0.0 offline example composes spatial-temporal out-of-fold beta calibration with AdaSTEM. It uses a small synthetic fixture, two ensembles and two workers, and compares native streaming/matrix and memmap predictions. Install the monitoring extra for notebook 10 and the tuning extra to construct the optional XGBoost recipe shown here.

The figures and metrics below are newly executed API demonstrations. They do not establish a best calibration method or validate copula inference. The older examples are retained separately; see the execution-status overview.

In [1]:

import sys
from pathlib import Path
sys.path.insert(0, str(Path.cwd()))
import tempfile
import time
import platform
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.special import expit
from sklearn.ensemble import RandomForestClassifier
from IPython.display import Image, display
import stemflow
from stemflow.model.AdaSTEM import AdaSTEMClassifier
from stemflow.model_selection import ST_train_test_split
from stemflow.lazyloading import training_databases
from stemflow.calibration import (BetaCalibrator, CrossValidatedCalibratedClassifier,
    SpatioTemporalKFold, calibration_scores, plot_calibration_diagnostics)
from stemflow.tuning import local_beta_xgboost
print("stemflow", stemflow.__version__, "Python", platform.python_version())
from matplotlib import font_manager
try:
    font_manager.findfont("Arial", fallback_to_default=False)
    plt.rcParams["font.family"] = "Arial"
    print("Plot font: Arial")
except ValueError:
    print("Arial is unavailable; install it to reproduce the reference typography. Using the local Matplotlib default.")
plt.rcParams.update({"font.size": 12, "axes.labelsize": 13,
    "xtick.labelsize": 11, "ytick.labelsize": 11, "legend.fontsize": 10,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": False, "pdf.fonttype": 42, "svg.fonttype": "none"})

import sys from pathlib import Path sys.path.insert(0, str(Path.cwd())) import tempfile import time import platform import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy.special import expit from sklearn.ensemble import RandomForestClassifier from IPython.display import Image, display import stemflow from stemflow.model.AdaSTEM import AdaSTEMClassifier from stemflow.model_selection import ST_train_test_split from stemflow.lazyloading import training_databases from stemflow.calibration import (BetaCalibrator, CrossValidatedCalibratedClassifier, SpatioTemporalKFold, calibration_scores, plot_calibration_diagnostics) from stemflow.tuning import local_beta_xgboost print("stemflow", stemflow.__version__, "Python", platform.python_version()) from matplotlib import font_manager try: font_manager.findfont("Arial", fallback_to_default=False) plt.rcParams["font.family"] = "Arial" print("Plot font: Arial") except ValueError: print("Arial is unavailable; install it to reproduce the reference typography. Using the local Matplotlib default.") plt.rcParams.update({"font.size": 12, "axes.labelsize": 13, "xtick.labelsize": 11, "ytick.labelsize": 11, "legend.fontsize": 10, "axes.spines.top": False, "axes.spines.right": False, "axes.grid": False, "pdf.fonttype": 42, "svg.fonttype": "none"})

stemflow 2.0.0 Python 3.13.15
Plot font: Arial

In [2]:

rng = np.random.default_rng(42)
X = pd.DataFrame({"longitude": rng.uniform(-100, -70, 2500),
    "latitude": rng.uniform(25, 50, 2500), "DOY": rng.integers(121, 201, 2500),
    "environment": rng.normal(size=2500)})
truth = expit(-2.5 + X.environment + .8*np.sin((X.DOY-121)/80*np.pi))
y = pd.Series(rng.binomial(1, truth), index=X.index, name="occ")
X_train, X_test, y_train, y_test = ST_train_test_split(X, y.to_frame(),
    Spatio_blocks_count=10, Temporal_blocks_count=5, test_size=.25, random_state=42)
print(f"Training rows/events: {len(X_train)}/{int(y_train.occ.sum())}")
print(f"Held-out rows/events: {len(X_test)}/{int(y_test.occ.sum())}")

rng = np.random.default_rng(42) X = pd.DataFrame({"longitude": rng.uniform(-100, -70, 2500), "latitude": rng.uniform(25, 50, 2500), "DOY": rng.integers(121, 201, 2500), "environment": rng.normal(size=2500)}) truth = expit(-2.5 + X.environment + .8*np.sin((X.DOY-121)/80*np.pi)) y = pd.Series(rng.binomial(1, truth), index=X.index, name="occ") X_train, X_test, y_train, y_test = ST_train_test_split(X, y.to_frame(), Spatio_blocks_count=10, Temporal_blocks_count=5, test_size=.25, random_state=42) print(f"Training rows/events: {len(X_train)}/{int(y_train.occ.sum())}") print(f"Held-out rows/events: {len(X_test)}/{int(y_test.occ.sum())}")

Training rows/events: 1896/309
Held-out rows/events: 604/103

Each stixel receives three retained fold models and one identity-regularized beta map. Longitude/latitude metadata are supplied automatically to the CV wrapper and excluded from forest predictors; DOY remains a predictor. The generic wrapper can also be fitted directly outside AdaSTEM. Set sample_weights_for_classifier=False for this unweighted strategy.

In [3]:

base = CrossValidatedCalibratedClassifier(
    RandomForestClassifier(n_estimators=15, max_depth=5, min_samples_leaf=10,
        n_jobs=1, random_state=42),
    BetaCalibrator(ridge=100), SpatioTemporalKFold(
        n_splits=3, Spatio_blocks_count=5, Temporal_blocks_count=5))
with tempfile.TemporaryDirectory(prefix="stemflow_example_") as scratch:
    model = AdaSTEMClassifier(base_model=base, ensemble_fold=2,
        min_ensemble_required=1, points_lower_threshold=200,
        grid_len_upper_threshold=25, grid_len_lower_threshold=5,
        temporal_start=121, temporal_end=201, temporal_step=40,
        temporal_bin_interval=40, temporal_bin_start_jitter=0,
        sample_weights_for_classifier=False, ensemble_bootstrap=True,
        lazy_loading=True, lazy_loading_dir=str(Path(scratch)/"models"),
        n_jobs=2, random_state=42, save_gridding_plot=False, verbosity=0)
    started = time.perf_counter()
    with training_databases(X_train, y_train, directory=scratch) as (xp, yp):
        model.fit(xp, yp, temporal_window_prequery=True,
            training_data_backend="memmap", inner_max_num_threads=1, persist_prediction_state=False)
    fit_seconds = time.perf_counter()-started
    started = time.perf_counter()
    probability, std = model.predict_proba(X_test, return_std=True,
        temporal_window_prequery=True)
    prediction_seconds = time.perf_counter()-started
    execution = model.prediction_execution_info_.copy()
    mapped, mapped_std = model.predict_proba(X_test, return_std=True, n_jobs=2,
        temporal_window_prequery=True, prediction_backend="memmap")
    np.testing.assert_allclose(probability, mapped, atol=1e-12, equal_nan=True)
    np.testing.assert_allclose(std, mapped_std, atol=1e-12, equal_nan=True)
    matrix, matrix_std = model.predict_proba(X_test, return_std=True,
        temporal_window_prequery=True, streaming=False)
    np.testing.assert_allclose(probability, matrix, atol=1e-12, equal_nan=True)
    np.testing.assert_allclose(std, matrix_std, atol=1e-12, equal_nan=True)
    diagnostics = model.calibration_diagnostics(X_test, y_test, probability[:, 1])
    plots = Path(scratch)/"plots"
    plot_calibration_diagnostics(diagnostics, plots, max_curves=100)
    display(Image(filename=str(plots/"calibration_diagnostics.png")))
    display(Image(filename=str(plots/"stixel_calibration_overlay.png")))
    del model
print(f"Measured fit: {fit_seconds:.3f}s; streamed prediction: {prediction_seconds:.3f}s")
print(execution)
print("Native streaming/matrix and memmap predictions agree; temporary files removed.")

base = CrossValidatedCalibratedClassifier( RandomForestClassifier(n_estimators=15, max_depth=5, min_samples_leaf=10, n_jobs=1, random_state=42), BetaCalibrator(ridge=100), SpatioTemporalKFold( n_splits=3, Spatio_blocks_count=5, Temporal_blocks_count=5)) with tempfile.TemporaryDirectory(prefix="stemflow_example_") as scratch: model = AdaSTEMClassifier(base_model=base, ensemble_fold=2, min_ensemble_required=1, points_lower_threshold=200, grid_len_upper_threshold=25, grid_len_lower_threshold=5, temporal_start=121, temporal_end=201, temporal_step=40, temporal_bin_interval=40, temporal_bin_start_jitter=0, sample_weights_for_classifier=False, ensemble_bootstrap=True, lazy_loading=True, lazy_loading_dir=str(Path(scratch)/"models"), n_jobs=2, random_state=42, save_gridding_plot=False, verbosity=0) started = time.perf_counter() with training_databases(X_train, y_train, directory=scratch) as (xp, yp): model.fit(xp, yp, temporal_window_prequery=True, training_data_backend="memmap", inner_max_num_threads=1, persist_prediction_state=False) fit_seconds = time.perf_counter()-started started = time.perf_counter() probability, std = model.predict_proba(X_test, return_std=True, temporal_window_prequery=True) prediction_seconds = time.perf_counter()-started execution = model.prediction_execution_info_.copy() mapped, mapped_std = model.predict_proba(X_test, return_std=True, n_jobs=2, temporal_window_prequery=True, prediction_backend="memmap") np.testing.assert_allclose(probability, mapped, atol=1e-12, equal_nan=True) np.testing.assert_allclose(std, mapped_std, atol=1e-12, equal_nan=True) matrix, matrix_std = model.predict_proba(X_test, return_std=True, temporal_window_prequery=True, streaming=False) np.testing.assert_allclose(probability, matrix, atol=1e-12, equal_nan=True) np.testing.assert_allclose(std, matrix_std, atol=1e-12, equal_nan=True) diagnostics = model.calibration_diagnostics(X_test, y_test, probability[:, 1]) plots = Path(scratch)/"plots" plot_calibration_diagnostics(diagnostics, plots, max_curves=100) display(Image(filename=str(plots/"calibration_diagnostics.png"))) display(Image(filename=str(plots/"stixel_calibration_overlay.png"))) del model print(f"Measured fit: {fit_seconds:.3f}s; streamed prediction: {prediction_seconds:.3f}s") print(execution) print("Native streaming/matrix and memmap predictions agree; temporary files removed.")

Measured fit: 2.525s; streamed prediction: 1.485s
{'aggregation_backend': 'streaming', 'accumulator_bytes': 12080, 'rows': 604, 'ensembles': 2}
Native streaming/matrix and memmap predictions agree; temporary files removed.

/private/tmp/stemflow_ci_fix_20261028/stemflow/utils/validation.py:415: UserWarning: The grid_len_upper_threshold is larger than the scale of longitude and latitude (x and y). Be sure if this is desired.
  warnings.warn(
/private/tmp/stemflow_ci_fix_20261028/stemflow/utils/validation.py:161: UserWarning: Input X is pandas dataframe; self.max_mem is not used.
  warnings.warn('Input X is pandas dataframe; self.max_mem is not used.')
/private/tmp/stemflow_ci_fix_20261028/stemflow/model/AdaSTEM.py:1061: UserWarning: There are 8.94% points (54 points) falling out of predictable range.
  warnings.warn(f"There are {100 * missing / len(mean):.3g}% points ({missing} points) falling out of predictable range.")
/private/tmp/stemflow_ci_fix_20261028/stemflow/utils/validation.py:161: UserWarning: Input X is pandas dataframe; self.max_mem is not used.
  warnings.warn('Input X is pandas dataframe; self.max_mem is not used.')
/private/tmp/stemflow_ci_fix_20261028/stemflow/lazyloading/shared_prediction.py:396: UserWarning: There are 8.94% points (54 points) falling out of predictable range.
  warnings.warn(f'There are {100 * missing / len(mean):.3g}% points ({missing} points) '
/private/tmp/stemflow_ci_fix_20261028/stemflow/utils/validation.py:161: UserWarning: Input X is pandas dataframe; self.max_mem is not used.
  warnings.warn('Input X is pandas dataframe; self.max_mem is not used.')
/private/tmp/stemflow_ci_fix_20261028/stemflow/model/AdaSTEM.py:1255: UserWarning: There are 0.08940397350993377% points (54 points) falling out of predictable range.
  warnings.warn(f"There are {nan_frac}% points ({nan_count} points) falling out of predictable range.")

In [4]:

overall = diagnostics["scores"].query("level == 'overall'")
display(overall[["n_valid", "detections", "coverage", "AUC", "Brier", "log_loss",
    "quantile_ECE_pp", "quantile_ECE_relative", "quantile_zero_event_mass"]])
display(diagnostics["summary"].query("level == 'stixel' and scope == 'all'"))
valid = np.isfinite(probability[:, 1])
print("Mean absolute probability error against synthetic truth:",
    np.mean(np.abs(probability[valid, 1]-truth.loc[X_test.index].to_numpy()[valid])))

overall = diagnostics["scores"].query("level == 'overall'") display(overall[["n_valid", "detections", "coverage", "AUC", "Brier", "log_loss", "quantile_ECE_pp", "quantile_ECE_relative", "quantile_zero_event_mass"]]) display(diagnostics["summary"].query("level == 'stixel' and scope == 'all'")) valid = np.isfinite(probability[:, 1]) print("Mean absolute probability error against synthetic truth:", np.mean(np.abs(probability[valid, 1]-truth.loc[X_test.index].to_numpy()[valid])))

Mean absolute probability error against synthetic truth: 0.05999128076839138

	n_valid	detections	coverage	AUC	Brier	log_loss	quantile_ECE_pp	quantile_ECE_relative	quantile_zero_event_mass
0	550	96	0.910596	0.677783	0.136654	0.439788	4.703917	1.382759	0.0

	scope	level	metric	n_groups	q25	median	q75	p90
39	all	stixel	ECE_pp	11	6.567155	7.738166	9.571750	10.519768
40	all	stixel	quantile_ECE_pp	11	9.624685	10.335133	10.891143	11.779608
41	all	stixel	ECE_relative	11	1.557367	1.703946	1.897413	73.944759
42	all	stixel	quantile_ECE_relative	11	5.817360	25.741909	102.287045	134.734561
43	all	stixel	quantile_relative_floor_1e8	11	10.479059	78.909127	582.069958	850.117607
44	all	stixel	quantile_zero_event_mass	11	0.125563	0.243243	0.377717	0.400000
45	all	stixel	bias_pp	11	-2.277472	1.782357	7.165384	7.723256
46	all	stixel	rate_ratio	11	0.889015	1.108215	1.586803	1.635023
47	all	stixel	AUC	11	0.627092	0.651838	0.679302	0.693510
48	all	stixel	average_precision	11	0.228585	0.369791	0.372404	0.402976
49	all	stixel	Brier	11	0.113479	0.125869	0.151032	0.166122
50	all	stixel	log_loss	11	0.378724	0.409223	0.490562	0.526650
51	all	stixel	coverage	11	1.000000	1.000000	1.000000	1.000000

The red overlay averages actual common probability intervals, including observed zeros before plotting. Zero observations break logarithmic lines. Training-supported stixels use training counts, not held-out labels. Sparse bins and overlapping stixels limit interpretation of relative ECE; support counts, expected detections, loss and coverage should accompany it. The regression wrapper calibrates point means only. Detailed profiling is opt-in through STEMFLOW_PROFILE_DIR, STEMFLOW_PROFILE_LINES and STEMFLOW_PROFILE_ALLOCATIONS; run those profiles separately from timing comparisons. See local calibration and tuning and parallel execution.

Use the XGBoost recipe instead¶

The fitted example above uses a fixed small forest with blocked calibration. The following constructs the XGBoost alternative without launching another fit. Set tuning=True to add nested search; trial/fold/tree counts are independent controls. A generic tuner with another estimator/search space can be wrapped in the same calibration class. These settings illustrate the API, not a universally optimal ridge or search configuration.

In [5]:

recipe = local_beta_xgboost(
    tuning=False, ridge=100, n_trials=4,
    n_estimators=20, tuning_n_estimators=8,
    n_splits=3, tuning_n_splits=2,
    Spatio_blocks_count=5, Temporal_blocks_count=5,
    max_trials=10, tuning_verbose=False,
)
print(type(recipe).__name__)
print("Calibration folds:", recipe.cv.n_splits)
print("Spatial/temporal block counts:", recipe.cv.Spatio_blocks_count, recipe.cv.Temporal_blocks_count)
print("Ridge:", recipe.calibrator.ridge)

recipe = local_beta_xgboost( tuning=False, ridge=100, n_trials=4, n_estimators=20, tuning_n_estimators=8, n_splits=3, tuning_n_splits=2, Spatio_blocks_count=5, Temporal_blocks_count=5, max_trials=10, tuning_verbose=False, ) print(type(recipe).__name__) print("Calibration folds:", recipe.cv.n_splits) print("Spatial/temporal block counts:", recipe.cv.Spatio_blocks_count, recipe.cv.Temporal_blocks_count) print("Ridge:", recipe.calibrator.ridge)

CrossValidatedCalibratedClassifier
Calibration folds: 3
Spatial/temporal block counts: 5 5
Ridge: 100