diff --git a/examples/tutorials/model_explainability.ipynb b/examples/tutorials/model_explainability.ipynb index 78b0acc5c..81f1cd672 100644 --- a/examples/tutorials/model_explainability.ipynb +++ b/examples/tutorials/model_explainability.ipynb @@ -31,6 +31,8 @@ "import warnings\n", "from typing import cast\n", "\n", + "from openstef_meta.presets import EnsembleForecastingWorkflowConfig, create_ensemble_forecasting_workflow\n", + "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "from openstef_core.testing import configure_notebook_display, setup_notebook_logging\n", @@ -174,7 +176,7 @@ }, "outputs": [], "source": [ - "from openstef_models.explainability import ExplainableForecaster\n", + "from openstef_models.explainability import ExplainableForecaster, ModelContributionsPlotter\n", "from openstef_models.models.forecasting_model import ForecastingModel\n", "\n", "forecaster = cast(ForecastingModel, workflow.model).forecaster\n", @@ -317,6 +319,70 @@ "fig.show()" ] }, + { + "cell_type": "markdown", + "id": "5118a84f", + "metadata": {}, + "source": [ + "## Ensemble model weights\n", + "\n", + "For each forecast quantile, a classifier estimates which base model is likely\n", + "to perform best. Its predicted probabilities become nonnegative model weights\n", + "that sum to one at each timestep. The top panel shows the ensemble forecast\n", + "for the selected quantile; the stacked area below shows each model's weight.\n", + "The plot defaults to P50 when no quantile is specified." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "621e0f94", + "metadata": { + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "ensemble_workflow = create_ensemble_forecasting_workflow(\n", + " config=EnsembleForecastingWorkflowConfig(\n", + " model_id=\"example_gblinear_lgbm_ensemble\",\n", + " ensemble_type=\"learned_weights\",\n", + " base_models=[\"gblinear\", \"lgbm\"],\n", + " combiner_model=\"lgbm\",\n", + " quantiles=[Q(0.1), Q(0.5), Q(0.9)],\n", + " horizons=[LeadTime.from_string(\"PT36H\")],\n", + " mlflow_storage=None,\n", + " temperature_column=\"temperature_2m\",\n", + " relative_humidity_column=\"relative_humidity_2m\",\n", + " wind_speed_column=\"wind_speed_10m\",\n", + " radiation_column=\"shortwave_radiation\",\n", + " pressure_column=\"surface_pressure\",\n", + " )\n", + ")\n", + "ensemble_result = ensemble_workflow.fit(train_dataset)\n", + "ensemble_forecast, model_weights = ensemble_workflow.model.predict_with_contributions(\n", + " predict_dataset, forecast_start=train_end\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df7a1d98", + "metadata": {}, + "outputs": [], + "source": [ + "fig = ModelContributionsPlotter.plot_stacked_area(\n", + " model_weights,\n", + " forecast=ensemble_forecast,\n", + " quantile=Q(0.1),\n", + ")\n", + "fig.update_layout(height=650)\n", + "fig.show()\n", + "# fig.write_html(\"ensemble_model_weights.html\", include_plotlyjs=True)" + ] + }, { "cell_type": "markdown", "id": "dc95f759", @@ -337,9 +403,21 @@ "notebook_metadata_filter": "-jupytext.text_representation.jupytext_version" }, "kernelspec": { - "display_name": ".venv", + "display_name": "openstef (3.12.10)", "language": "python", "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" } }, "nbformat": 4, diff --git a/examples/tutorials/model_explainability.py b/examples/tutorials/model_explainability.py index 5ddf925ba..efb4e149d 100644 --- a/examples/tutorials/model_explainability.py +++ b/examples/tutorials/model_explainability.py @@ -8,7 +8,7 @@ # format_name: percent # format_version: '1.3' # kernelspec: -# display_name: .venv +# display_name: openstef (3.12.10) # language: python # name: python3 # --- @@ -23,6 +23,8 @@ import warnings from typing import cast +from openstef_meta.presets import EnsembleForecastingWorkflowConfig, create_ensemble_forecasting_workflow + warnings.filterwarnings("ignore") from openstef_core.testing import configure_notebook_display, setup_notebook_logging @@ -122,7 +124,7 @@ # tiles represent more influential features. # %% tags=["hide-input"] -from openstef_models.explainability import ExplainableForecaster +from openstef_models.explainability import ExplainableForecaster, ModelContributionsPlotter from openstef_models.models.forecasting_model import ForecastingModel forecaster = cast(ForecastingModel, workflow.model).forecaster @@ -193,6 +195,47 @@ fig.update_layout(title="Prediction decomposition (timestep 48)", height=500) fig.show() +# %% [markdown] +# ## Ensemble model weights +# +# For each forecast quantile, a classifier estimates which base model is likely +# to perform best. Its predicted probabilities become nonnegative model weights +# that sum to one at each timestep. The top panel shows the ensemble forecast +# for the selected quantile; the stacked area below shows each model's weight. +# The plot defaults to P50 when no quantile is specified. + +# %% tags=["hide-input"] +ensemble_workflow = create_ensemble_forecasting_workflow( + config=EnsembleForecastingWorkflowConfig( + model_id="example_gblinear_lgbm_ensemble", + ensemble_type="learned_weights", + base_models=["gblinear", "lgbm"], + combiner_model="lgbm", + quantiles=[Q(0.1), Q(0.5), Q(0.9)], + horizons=[LeadTime.from_string("PT36H")], + mlflow_storage=None, + temperature_column="temperature_2m", + relative_humidity_column="relative_humidity_2m", + wind_speed_column="wind_speed_10m", + radiation_column="shortwave_radiation", + pressure_column="surface_pressure", + ) +) +ensemble_result = ensemble_workflow.fit(train_dataset) +ensemble_forecast, model_weights = ensemble_workflow.model.predict_with_contributions( + predict_dataset, forecast_start=train_end +) + +# %% +fig = ModelContributionsPlotter.plot_stacked_area( + model_weights, + forecast=ensemble_forecast, + quantile=Q(0.1), +) +fig.update_layout(height=650) +fig.show() +# fig.write_html("ensemble_model_weights.html", include_plotlyjs=True) + # %% [markdown] # ## Next steps # diff --git a/packages/openstef-meta/src/openstef_meta/models/ensemble_forecasting_model.py b/packages/openstef-meta/src/openstef_meta/models/ensemble_forecasting_model.py index 01672b190..f35c13cde 100644 --- a/packages/openstef-meta/src/openstef_meta/models/ensemble_forecasting_model.py +++ b/packages/openstef-meta/src/openstef_meta/models/ensemble_forecasting_model.py @@ -673,5 +673,41 @@ def predict_contributions( original_data=data, ) + def predict_with_contributions( + self, + data: TimeSeriesDataset, + forecast_start: datetime | None = None, + ) -> tuple[ForecastDataset, TimeSeriesDataset]: + """Generate forecasts and per-model contributions in a single pass. + + Args: + data: Input time series dataset. + forecast_start: Optional start time for forecasts. + + Returns: + A tuple of ``(forecast, contributions)`` where *forecast* is the + ensemble prediction (a :class:`ForecastDataset`) and *contributions* + is the :class:`TimeSeriesDataset` returned by the configured combiner. + Its meaning depends on the combiner; for learned weights, columns + contain per-model selection weights. + + Raises: + NotFittedError: If the model has not been fitted yet. + """ + if not self.is_fitted: + raise NotFittedError(self.__class__.__name__) + + ensemble_predictions = self._predict_forecasters(data=data, forecast_start=forecast_start) + + forecast = self._predict_transform_combiner( + ensemble_dataset=ensemble_predictions, + original_data=data, + ) + contributions = self._predict_contributions_combiner( + ensemble_dataset=ensemble_predictions, + original_data=data, + ) + return forecast, contributions + __all__ = ["EnsembleForecastingModel", "EnsembleModelFitResult", "ModelFitResult"] diff --git a/packages/openstef-meta/tests/unit/models/test_ensemble_forecasting_model.py b/packages/openstef-meta/tests/unit/models/test_ensemble_forecasting_model.py index b4592f944..5b7b99a57 100644 --- a/packages/openstef-meta/tests/unit/models/test_ensemble_forecasting_model.py +++ b/packages/openstef-meta/tests/unit/models/test_ensemble_forecasting_model.py @@ -144,6 +144,23 @@ def test_forecasting_model__predict__raises_error_when_not_fitted( model.predict(data=sample_timeseries_dataset) +def test_predict_with_contributions_matches_separate_calls( + sample_timeseries_dataset: TimeSeriesDataset, model: EnsembleForecastingModel +) -> None: + """Combined prediction matches separate forecast and contribution calls.""" + model.fit(data=sample_timeseries_dataset) + forecast_start = datetime.fromisoformat("2025-01-01T12:00:00") + + forecast_combined, contributions_combined = model.predict_with_contributions( + sample_timeseries_dataset, forecast_start=forecast_start + ) + forecast_separate = model.predict(sample_timeseries_dataset, forecast_start=forecast_start) + contributions_separate = model.predict_contributions(sample_timeseries_dataset, forecast_start=forecast_start) + + pd.testing.assert_frame_equal(forecast_combined.data, forecast_separate.data) + pd.testing.assert_frame_equal(contributions_combined.data, contributions_separate.data) + + def test_forecasting_model__score__returns_metrics( sample_timeseries_dataset: TimeSeriesDataset, model: EnsembleForecastingModel ): diff --git a/packages/openstef-models/src/openstef_models/explainability/__init__.py b/packages/openstef-models/src/openstef_models/explainability/__init__.py index bfc562145..3170cba98 100644 --- a/packages/openstef-models/src/openstef_models/explainability/__init__.py +++ b/packages/openstef-models/src/openstef_models/explainability/__init__.py @@ -8,11 +8,12 @@ """ from .mixins import ContributionsMixin, ExplainableForecaster -from .plotters import ContributionsPlotter, FeatureImportancePlotter +from .plotters import ContributionsPlotter, FeatureImportancePlotter, ModelContributionsPlotter __all__ = [ "ContributionsMixin", "ContributionsPlotter", "ExplainableForecaster", "FeatureImportancePlotter", + "ModelContributionsPlotter", ] diff --git a/packages/openstef-models/src/openstef_models/explainability/plotters/__init__.py b/packages/openstef-models/src/openstef_models/explainability/plotters/__init__.py index 027485f79..a79cb5522 100644 --- a/packages/openstef-models/src/openstef_models/explainability/plotters/__init__.py +++ b/packages/openstef-models/src/openstef_models/explainability/plotters/__init__.py @@ -10,8 +10,10 @@ from .contributions_plotter import ContributionsPlotter from .feature_importance_plotter import FeatureImportancePlotter +from .model_contributions_plotter import ModelContributionsPlotter __all__ = [ "ContributionsPlotter", "FeatureImportancePlotter", + "ModelContributionsPlotter", ] diff --git a/packages/openstef-models/src/openstef_models/explainability/plotters/common.py b/packages/openstef-models/src/openstef_models/explainability/plotters/common.py new file mode 100644 index 000000000..bdca4faf3 --- /dev/null +++ b/packages/openstef-models/src/openstef_models/explainability/plotters/common.py @@ -0,0 +1,37 @@ +# SPDX-FileCopyrightText: 2026 Contributors to the OpenSTEF project +# +# SPDX-License-Identifier: MPL-2.0 + +"""Layout helpers for explainability plots.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import plotly.graph_objects as go +from plotly.subplots import make_subplots + +if TYPE_CHECKING: + import pandas as pd + + +class VerticalTwoPanelPlotLayout: + """Create vertically aligned subplots.""" + + @staticmethod + def make_two_panel_time_figure() -> go.Figure: + """Create the vertically aligned subplot layout for two subplots.""" + return make_subplots(rows=2, cols=1, shared_xaxes=True, row_heights=[0.2, 0.8], vertical_spacing=0.03) + + @staticmethod + def add_top_panel_line_plot( + fig: go.Figure, x: pd.Index, y: pd.Series, *, name: str = "Prediction", show_legend: bool = False + ) -> go.Figure: + """Add a lineplot. Can be used to add a top panel prediction line plot.""" + return fig.add_trace( + go.Scatter( + x=x, y=y, mode="lines", line={"color": "black", "width": 1.5}, name=name, showlegend=show_legend + ), + row=1, + col=1, + ) diff --git a/packages/openstef-models/src/openstef_models/explainability/plotters/contributions_plotter.py b/packages/openstef-models/src/openstef_models/explainability/plotters/contributions_plotter.py index 5da8f021d..b533f3fe5 100644 --- a/packages/openstef-models/src/openstef_models/explainability/plotters/contributions_plotter.py +++ b/packages/openstef-models/src/openstef_models/explainability/plotters/contributions_plotter.py @@ -9,9 +9,9 @@ from typing import TYPE_CHECKING, cast import plotly.graph_objects as go -from plotly.subplots import make_subplots from openstef_core.datasets import TimeSeriesDataset # noqa: TC001 +from openstef_models.explainability.plotters.common import VerticalTwoPanelPlotLayout if TYPE_CHECKING: import pandas as pd @@ -69,26 +69,9 @@ def plot_heatmap( if bias is not None: prediction += bias - fig = make_subplots( - rows=2, - cols=1, - shared_xaxes=True, - row_heights=[0.2, 0.8], - vertical_spacing=0.03, - ) + fig = VerticalTwoPanelPlotLayout.make_two_panel_time_figure() - fig.add_trace( - go.Scatter( - x=df.index, - y=prediction, - mode="lines", - name="Prediction", - line={"color": "black", "width": 1.5}, - showlegend=False, - ), - row=1, - col=1, - ) + VerticalTwoPanelPlotLayout.add_top_panel_line_plot(fig, df.index, prediction) fig.add_trace(heatmap, row=2, col=1) fig.update_layout( diff --git a/packages/openstef-models/src/openstef_models/explainability/plotters/model_contributions_plotter.py b/packages/openstef-models/src/openstef_models/explainability/plotters/model_contributions_plotter.py new file mode 100644 index 000000000..e66202595 --- /dev/null +++ b/packages/openstef-models/src/openstef_models/explainability/plotters/model_contributions_plotter.py @@ -0,0 +1,108 @@ +# SPDX-FileCopyrightText: 2026 Contributors to the OpenSTEF project +# +# SPDX-License-Identifier: MPL-2.0 + +"""Visualizations for ensemble model-level contributions.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import plotly.graph_objects as go + +from openstef_core.datasets.validated_datasets import ENSEMBLE_COLUMN_SEP +from openstef_core.types import Q, Quantile +from openstef_models.explainability.plotters.common import ( + VerticalTwoPanelPlotLayout, +) + +if TYPE_CHECKING: + from openstef_core.datasets import ForecastDataset, TimeSeriesDataset + + +class ModelContributionsPlotter: + """Visualize learned per-model weights alongside the ensemble forecast.""" + + @staticmethod + def plot_stacked_area( + contributions: TimeSeriesDataset, + forecast: ForecastDataset, + quantile: Quantile = Q(0.5), + target_column: str = "load", + ) -> go.Figure: + """Create a stacked-area plot of per-model weights with the ensemble prediction. + + The top panel shows the ensemble prediction for the requested quantile. + The bottom panel shows each base model's selection weight as a stacked + area (summing to 1.0), so you can see how the combiner distributes + weight across models at each timestep. + + Args: + contributions: Per-model selection weights from a learned-weights + ``EnsembleForecastingModel``. Use the result of + ``predict_contributions()`` or the second element of + ``predict_with_contributions()``. Columns contain model weights + for each forecast quantile. + forecast: Output of + ``EnsembleForecastingModel.predict()`` or the first element + of ``EnsembleForecastingModel.predict_with_contributions()``. + Must contain a column for the requested quantile. + quantile: Forecast quantile to show in both panels. Default ``Q(0.5)``. + target_column: Name of the target column to exclude from + contributions. Default ``"load"``. + + Returns: + Plotly Figure with a prediction line (top) and stacked-area model + weights (bottom). + + Raises: + ValueError: If no contribution columns match the requested quantile. + KeyError: If the forecast dataset has no column for the requested + quantile. + """ + df = contributions.data.drop(columns=[target_column], errors="ignore") + + available_columns: dict[tuple[str, str], str] = {} + for col in df.columns: + model, separator, model_quantile = col.rpartition(ENSEMBLE_COLUMN_SEP) + if separator: + available_columns[(model, model_quantile)] = col + + selected_columns = [ + col for (model, model_quantile), col in available_columns.items() if model_quantile == quantile.format() + ] + trace_names = [col[: -len(ENSEMBLE_COLUMN_SEP + quantile.format())] for col in selected_columns] + + if not selected_columns: + available = sorted({q for _, q in available_columns}) + msg = f"No columns found for quantile {quantile.format()!r}. Available quantiles: {available}." + raise ValueError(msg) + + selected = df[selected_columns] + + forecast_series = forecast.data[quantile.format()] + + fig = VerticalTwoPanelPlotLayout.make_two_panel_time_figure() + VerticalTwoPanelPlotLayout.add_top_panel_line_plot(fig, forecast_series.index, forecast_series) + + for col, name in zip(selected_columns, trace_names, strict=True): + fig.add_trace( + go.Scatter( + x=df.index, + y=selected[col], + mode="lines", + stackgroup="contrib", + name=name, + ), + row=2, + col=1, + ) + + fig.update_layout( + yaxis_title=f"Prediction ({quantile.format().removeprefix('quantile_')})", + yaxis2_title=f"Model weight ({quantile.format().removeprefix('quantile_')})", + xaxis2_title="Time", + margin={"t": 30, "r": 10, "b": 40, "l": 60}, + ) + + return fig diff --git a/packages/openstef-models/tests/unit/explainability/plotters/test_model_contributions_plotter.py b/packages/openstef-models/tests/unit/explainability/plotters/test_model_contributions_plotter.py new file mode 100644 index 000000000..0a9e43d71 --- /dev/null +++ b/packages/openstef-models/tests/unit/explainability/plotters/test_model_contributions_plotter.py @@ -0,0 +1,152 @@ +# SPDX-FileCopyrightText: 2026 Contributors to the OpenSTEF project +# +# SPDX-License-Identifier: MPL-2.0 + +"""Tests for ModelContributionsPlotter.""" + +from datetime import timedelta +from pathlib import Path + +import numpy as np +import pandas as pd +import plotly.graph_objects as go +import pytest + +from openstef_core.datasets import ForecastDataset, TimeSeriesDataset +from openstef_core.types import Q +from openstef_models.explainability import ModelContributionsPlotter + +SAMPLE_INTERVAL = timedelta(minutes=15) + + +@pytest.fixture +def contributions_dataset() -> TimeSeriesDataset: + """Two model weights across two quantiles, summing to 1.0 per row.""" + index = pd.date_range("2025-01-01", periods=8, freq="15min") + gb_p50 = np.linspace(0.6, 0.4, 8) + data = pd.DataFrame( + { + "gblinear__quantile_P50": gb_p50, + "lgbm__quantile_P50": 1.0 - gb_p50, + "gblinear__quantile_P90": np.full(8, 0.5), + "lgbm__quantile_P90": np.full(8, 0.5), + "load": np.linspace(16, 30, 8), + }, + index=index, + ) + return TimeSeriesDataset(data=data, sample_interval=SAMPLE_INTERVAL) + + +@pytest.fixture +def forecast_dataset() -> ForecastDataset: + """Ensemble forecast with P50 and P90 quantile columns.""" + index = pd.date_range("2025-01-01", periods=8, freq="15min") + data = pd.DataFrame( + { + "load": np.linspace(16, 30, 8), + "quantile_P10": np.linspace(10, 20, 8), + "quantile_P50": np.linspace(15, 28, 8), + "quantile_P90": np.linspace(18, 32, 8), + }, + index=index, + ) + return ForecastDataset(data=data, sample_interval=SAMPLE_INTERVAL) + + +def test_returns_figure( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + result = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset) + assert isinstance(result, go.Figure) + + +def test_has_forecast_and_stacked_weights_in_separate_panels( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset, quantile=Q(0.5)) + assert len(fig.data) == 3 + assert fig.data[0].xaxis == "x" + assert all(trace.xaxis == "x2" for trace in fig.data[1:]) + assert all(trace.stackgroup == "contrib" for trace in fig.data[1:]) + + +def test_selects_requested_quantile( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset, quantile=Q(0.5)) + assert [trace.name for trace in fig.data[1:]] == ["gblinear", "lgbm"] + assert all("P90" not in trace.name for trace in fig.data[1:] if trace.name is not None) + + +def test_selected_weights_and_forecast_use_the_same_quantile( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area( + contributions_dataset, + forecast=forecast_dataset, + quantile=Q(0.9), + ) + + assert [trace.name for trace in fig.data[1:]] == ["gblinear", "lgbm"] + np.testing.assert_array_almost_equal(fig.data[0].y, forecast_dataset.data["quantile_P90"].to_numpy()) + expected_weights = contributions_dataset.data["gblinear__quantile_P90"].to_numpy() + np.testing.assert_array_almost_equal(fig.data[1].y, expected_weights) + assert fig.layout.yaxis.title.text == "Prediction (P90)" + assert fig.layout.yaxis2.title.text == "Model weight (P90)" + + +def test_top_panel_shows_prediction( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset, quantile=Q(0.5)) + top_trace = fig.data[0] + expected = forecast_dataset.data["quantile_P50"].to_numpy() + np.testing.assert_array_almost_equal(top_trace.y, expected) + + +def test_weights_sum_to_one( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset, quantile=Q(0.5)) + weight_traces = fig.data[1:] + stacked = np.sum([np.array(t.y, dtype=float) for t in weight_traces], axis=0) + np.testing.assert_array_almost_equal(stacked, 1.0) + + +def test_missing_quantile_raises( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + with pytest.raises(ValueError, match="quantile_P25"): + ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast=forecast_dataset, quantile=Q(0.25)) + + +def test_default_quantile_is_p50( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, +) -> None: + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast_dataset) + + np.testing.assert_array_almost_equal(fig.data[0].y, forecast_dataset.data["quantile_P50"].to_numpy()) + np.testing.assert_array_almost_equal(fig.data[1].y, contributions_dataset.data["gblinear__quantile_P50"].to_numpy()) + + +def test_figure_exports_as_standalone_html( + contributions_dataset: TimeSeriesDataset, + forecast_dataset: ForecastDataset, + tmp_path: Path, +) -> None: + html_path = tmp_path / "model-contributions.html" + fig = ModelContributionsPlotter.plot_stacked_area(contributions_dataset, forecast_dataset) + + fig.write_html(html_path, include_plotlyjs=True) + + html = html_path.read_text() + assert "plotly.js" in html + assert "gblinear" in html