Refactoring plot_horizon_forecast to handle preview and eval modes

Author: ogrisel

content/python_files/feature_engineering.py

@@ -339,7 +339,7 @@ def iqr(col, *, window_size: int):
 
 # %% [markdown]
 #
-# ## Remark lagged features engineering and system lag
+# ## Important remark about lagged features engineering and system lag
 #
 # When working with historical data, we often have access to all the past
 # measurements in the dataset. However, when we want to use the lagged features
@@ -475,7 +475,7 @@ def define_prediction_time_range(prediction_start_time, prediction_end_time):
 
 prediction_time = define_prediction_time_range(
     prediction_start_time, prediction_end_time
-)
+).skb.subsample(n=1000, how="head")
 prediction_time
 
 
@@ -668,10 +668,10 @@ def build_targets(prediction_time, electricity, horizons):
 altair.Chart(
     pl.concat(
         [
-            targets.skb.eval(),
+            targets.skb.preview(),
             hgbr_predictions.rename(
                 {target_column_name: predicted_target_column_name}
-            ).skb.eval(),
+            ).skb.preview(),
         ],
         how="horizontal",
     ).tail(24 * 7)
@@ -742,7 +742,14 @@ def build_targets(prediction_time, electricity, horizons):
 
 
 # %%
-def collect_cv_predictions(pipelines, cv_splitter, predictions, prediction_time):
+def collect_cv_predictions(
+    pipelines,
+    cv_splitter,
+    predictions,
+    prediction_time,
+    method_name="predict",
+    method_kwargs=None,
+):
     index_generator = cv_splitter.split(prediction_time.skb.eval())
 
     def splitter(X, y, index_generator):
@@ -752,9 +759,13 @@ def splitter(X, y, index_generator):
         return X[train_idx], X[test_idx], y[train_idx], y[test_idx]
 
     results = []
+    if method_kwargs is None:
+        method_kwargs = {}
+
     for (_, test_idx), pipeline in zip(
         cv_splitter.split(prediction_time.skb.eval()), pipelines
     ):
+        method = getattr(pipeline, method_name)
         split = predictions.skb.train_test_split(
             predictions.skb.get_data(),
             splitter=splitter,
@@ -765,7 +776,7 @@ def splitter(X, y, index_generator):
                 {
                     "prediction_time": prediction_time.skb.eval()[test_idx],
                     "load_mw": split["y_test"],
-                    "predicted_load_mw": pipeline.predict(split["test"]),
+                    "predicted_load_mw": method(split["test"], **method_kwargs),
                 }
             )
         )
@@ -940,10 +951,10 @@ def splitter(X, y, index_generator):
 altair.Chart(
     pl.concat(
         [
-            targets.skb.eval(),
+            targets.skb.preview(),
             predictions_ridge.rename(
                 {target_column_name: predicted_target_column_name}
-            ).skb.eval(),
+            ).skb.preview(),
         ],
         how="horizontal",
     ).tail(24 * 7)
@@ -1114,25 +1125,25 @@ def splitter(X, y, index_generator):
 )
 
 # %%
-plot_at_time = datetime.datetime(2025, 5, 24, 0, 0, tzinfo=datetime.timezone.utc)
+plot_at_time = datetime.datetime(2021, 4, 19, 0, 0, tzinfo=datetime.timezone.utc)
 historical_timedelta = datetime.timedelta(hours=24 * 5)
 plot_horizon_forecast(
     targets,
     named_predictions,
     plot_at_time,
     historical_timedelta,
     target_column_name_pattern,
-)
+).skb.preview()
 
 # %%
-plot_at_time = datetime.datetime(2025, 5, 25, 0, 0, tzinfo=datetime.timezone.utc)
+plot_at_time = datetime.datetime(2021, 4, 20, 0, 0, tzinfo=datetime.timezone.utc)
 plot_horizon_forecast(
     targets,
     named_predictions,
     plot_at_time,
     historical_timedelta,
     target_column_name_pattern,
-)
+).skb.preview()
 
 # %%
 from sklearn.metrics import r2_score
@@ -1208,6 +1219,7 @@ def scoring(regressor, X, y):
 # TODO: Exercise using RandomForestRegressor
 from sklearn.ensemble import RandomForestRegressor
 
+
 multioutput_predictions_rf = features_with_dropped_cols.skb.apply(
     RandomForestRegressor(min_samples_leaf=30, random_state=0, n_jobs=-1),
     y=targets.skb.drop(cols=["prediction_time", "load_mw"]).skb.mark_as_y(),
@@ -1603,7 +1615,7 @@ def binned_coverage(y_true_folds, y_quantile_low, y_quantile_high, n_bins=10):
 
 # %% [markdown]
 #
-# ## Reliability diagram for quantile regression
+# ## Reliability diagrams for quantile regression
 
 # %%
 plot_reliability_diagram(
@@ -1674,6 +1686,7 @@ def fit(self, X, y):
             strategy="quantile",
             subsample=200_000,
             encode="ordinal",
+            quantile_method="averaged_inverted_cdf",
             random_state=random_state,
         )
 
@@ -1716,3 +1729,39 @@ def predict_quantiles(self, X, quantiles=(0.05, 0.5, 0.95)):
 
     def predict(self, X):
         return self.predict_quantiles(X, self.quantile).ravel()
+
+
+# %%
+from sklearn.ensemble import HistGradientBoostingClassifier
+from threadpoolctl import threadpool_limits
+
+
+# with threadpool_limits(1):
+if True:
+    predictions_bqr = features_with_dropped_cols.skb.apply(
+        BinnedQuantileRegressor(
+            RandomForestClassifier(
+                n_jobs=-1, n_estimators=200, min_samples_leaf=5, random_state=0
+            ),
+            # HistGradientBoostingClassifier(random_state=0),
+            n_bins=30,
+        ),
+        y=target,
+    )
+
+# %%
+predictions_bqr
+
+# %%
+cv_results_bqr = predictions_bqr.skb.cross_validate(
+    cv=ts_cv_5,
+    scoring={
+        "d2_pinball": make_scorer(d2_pinball_score, alpha=0.5),
+        "MAPE": make_scorer(mean_absolute_percentage_error),
+    },
+    return_pipeline=True,
+    verbose=1,
+    n_jobs=-1,
+)
+cv_results_bqr
+# %%

content/python_files/tutorial_helpers.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import polars as pl
+import polars.selectors as cs
 import altair
 import skrub
 
@@ -483,7 +484,7 @@ def plot_binned_residuals(cv_predictions, by="hour"):
         color="independent"
     )
 
-
+@skrub.deferred
 def plot_horizon_forecast(
     targets,
     named_predictions,
@@ -511,25 +512,23 @@ def plot_horizon_forecast(
     altair.Chart
         A chart with the true target and the forecast values for different horizons.
     """
-    merged_data = targets.skb.select(cols=["prediction_time", "load_mw"]).skb.concat(
-        [named_predictions], axis=1
+    merged_data = pl.concat(
+        [targets.select(pl.col("prediction_time"), pl.col("load_mw")), named_predictions],
+        how="horizontal",
     )
     start_time = plot_at_time - historical_timedelta
     end_time = plot_at_time + datetime.timedelta(
-        hours=named_predictions.skb.eval().shape[1]
+        hours=named_predictions.shape[1]
     )
     true_values_past = merged_data.filter(
         pl.col("prediction_time").is_between(start_time, plot_at_time, closed="both")
     ).rename({"load_mw": "Past true load"})
     true_values_future = merged_data.filter(
-        pl.col("prediction_time").is_between(plot_at_time, end_time, closed="both")
+        pl.col("prediction_time").is_between(plot_at_time, end_time, closed="right")
     ).rename({"load_mw": "Future true load"})
     predicted_record = (
-        merged_data.skb.select(
-            cols=skrub.selectors.filter_names(str.startswith, "predict")
-        )
+        merged_data.select(cs.starts_with("predict"))
         .row(by_predicate=pl.col("prediction_time") == plot_at_time, named=True)
-        .skb.eval()
     )
     forecast_values = pl.DataFrame(
         {
@@ -541,15 +540,14 @@ def plot_horizon_forecast(
         }
         for horizon in range(1, len(predicted_record))
     )
-
     true_values_past_chart = (
-        altair.Chart(true_values_past.skb.eval())
+        altair.Chart(true_values_past)
         .transform_fold(["Past true load"])
         .mark_line(tooltip=True)
         .encode(x="prediction_time:T", y="Past true load:Q", color="key:N")
     )
     true_values_future_chart = (
-        altair.Chart(true_values_future.skb.eval())
+        altair.Chart(true_values_future)
         .transform_fold(["Future true load"])
         .mark_line(tooltip=True)
         .encode(x="prediction_time:T", y="Future true load:Q", color="key:N")
@@ -562,4 +560,4 @@ def plot_horizon_forecast(
     )
     return (
         true_values_past_chart + true_values_future_chart + forecast_values_chart
-    ).interactive()
+    ).interactive()