Learning improvements / testing.

.github/copilot-instructions.md

@@ -203,9 +203,11 @@ eventdisplay-ml-plot-classification-gamma-efficiency --help
 
 **Model Architecture**:
 - **Stereo reconstruction**: Multi-output regression (XGBoost)
-  - Targets: `[MCxoff, MCyoff, MCe0]` (x offset, y offset, log energy)
+  - Targets: `[Xoff_residual, Yoff_residual, E_residual]` (residuals relative to DispBDT)
+  - Residuals computed as: MC truth - DispBDT prediction
+  - During inference: final prediction = DispBDT baseline + predicted residual
   - Single model handles all telescope multiplicities (2-4+ telescopes)
-  - Features: Telescope-level arrays + event-level parameters
+  - Features: Telescope-level arrays + event-level parameters (including DispBDT results)
 
 - **Classification**: Binary classification (XGBoost)
   - Target: Gamma vs hadron (implicit in training data split)

docs/changes/53.feature.md

@@ -0,0 +1,9 @@
+Improvement for learning and testing:
+
+- add evaluation set to training and save it to joblib.
+
+Add a training evaluation script to plot the training and evaluation metrics over the training iterations:
+
+```console
+python src/eventdisplay_ml/scripts/plot_training_evaluation.py --model_file tmp_cta_testing/stereo_south/dispdir_bdt.joblib --output_file test_training_curves.png
+```

pyproject.toml

@@ -64,6 +64,7 @@ scripts.eventdisplay-ml-apply-xgb-classify = "eventdisplay_ml.scripts.apply_xgb_
 scripts.eventdisplay-ml-apply-xgb-stereo = "eventdisplay_ml.scripts.apply_xgb_stereo:main"
 scripts.eventdisplay-ml-plot-classification-performance-metrics = "eventdisplay_ml.scripts.plot_classification_performance_metrics:main"
 scripts.eventdisplay-ml-plot-classification-gamma-efficiency = "eventdisplay_ml.scripts.plot_classification_gamma_efficiency:main"
+scripts.eventdisplay-ml-plot-training-evaluation = "eventdisplay_ml.scripts.plot_training_evaluation:main"
 scripts.eventdisplay-ml-train-xgb-classify = "eventdisplay_ml.scripts.train_xgb_classify:main"
 scripts.eventdisplay-ml-train-xgb-stereo = "eventdisplay_ml.scripts.train_xgb_stereo:main"
 

src/eventdisplay_ml/config.py

@@ -239,5 +239,8 @@ def configure_apply(analysis_type):
     )
     model_configs["energy_bins_log10_tev"] = par.get("energy_bins_log10_tev", [])
     model_configs["zenith_bins_deg"] = par.get("zenith_bins_deg", [])
+    if analysis_type == "stereo_analysis":
+        model_configs["target_mean"] = par.get("target_mean")
+        model_configs["target_std"] = par.get("target_std")
 
     return model_configs

src/eventdisplay_ml/data_processing.py

@@ -871,11 +871,36 @@ def load_training_data(model_configs, file_list, analysis_type):
                     max_tel_per_type=model_configs.get("max_tel_per_type", None),
                     preview_rows=model_configs.get("preview_rows", 20),
                 )
+
+                # Filter out events with invalid energy reconstruction for stereo training
+                if analysis_type == "stereo_analysis":
+                    n_before_erec_filter = len(df_flat)
+                    valid_erec_mask = (df_flat["ErecS"] > 0) & np.isfinite(df_flat["ErecS"])
+                    df_flat = df_flat[valid_erec_mask]
+                    n_removed_erec = n_before_erec_filter - len(df_flat)
+                    if n_removed_erec > 0:
+                        _logger.info(
+                            f"Removed {n_removed_erec} events with ErecS <= 0 or NaN "
+                            f"(fraction removed: {n_removed_erec / n_before_erec_filter:.4f})"
+                        )
+
                 if analysis_type == "stereo_analysis":
+                    mc_xoff = _to_numpy_1d(df["MCxoff"], np.float32)[valid_erec_mask]
+                    mc_yoff = _to_numpy_1d(df["MCyoff"], np.float32)[valid_erec_mask]
+                    mc_e0 = _to_numpy_1d(df["MCe0"], np.float32)[valid_erec_mask]
+
+                    disp_xoff = df_flat["Xoff_weighted_bdt"].values
+                    disp_yoff = df_flat["Yoff_weighted_bdt"].values
+                    disp_erec = df_flat["ErecS"].values
+
+                    # Compute log energies (ErecS already filtered > 0)
+                    mc_e0_log = np.where(mc_e0 > 0, np.log10(mc_e0), np.nan)
+                    disp_erec_log = np.log10(disp_erec)  # Safe since already filtered > 0
+
                     new_cols = {
-                        "MCxoff": _to_numpy_1d(df["MCxoff"], np.float32),
-                        "MCyoff": _to_numpy_1d(df["MCyoff"], np.float32),
-                        "MCe0": np.log10(_to_numpy_1d(df["MCe0"], np.float32)),
+                        "Xoff_residual": mc_xoff - disp_xoff,
+                        "Yoff_residual": mc_yoff - disp_yoff,
+                        "E_residual": mc_e0_log - disp_erec_log,
                     }
                 elif analysis_type == "classification":
                     new_cols = {
@@ -887,6 +912,22 @@ def load_training_data(model_configs, file_list, analysis_type):
                 for col_name, values in new_cols.items():
                     df_flat[col_name] = values
 
+                # Filter out events with NaN in residuals (can't train on these)
+                if analysis_type == "stereo_analysis":
+                    n_before_nan_filter = len(df_flat)
+                    valid_mask = (
+                        np.isfinite(df_flat["Xoff_residual"])
+                        & np.isfinite(df_flat["Yoff_residual"])
+                        & np.isfinite(df_flat["E_residual"])
+                    )
+                    df_flat = df_flat[valid_mask]
+                    n_removed = n_before_nan_filter - len(df_flat)
+                    if n_removed > 0:
+                        _logger.info(
+                            f"Removed {n_removed} events with NaN residuals "
+                            f"(fraction removed: {n_removed / n_before_nan_filter:.4f})"
+                        )
+
                 dfs.append(df_flat)
 
                 del df

src/eventdisplay_ml/evaluate.py

@@ -12,8 +12,6 @@
     mean_squared_error,
 )
 
-from eventdisplay_ml.features import target_features
-
 _logger = logging.getLogger(__name__)
 
 
@@ -69,16 +67,36 @@ def evaluate_classification_model(model, x_test, y_test, df, x_cols, name):
 
 
 def evaluate_regression_model(
-    model, x_test, y_test, df, x_cols, y_data, name, shap_per_energy=False
+    model, x_test, y_pred, y_test, df, x_cols, y_data, name, shap_per_energy=False
 ):
-    """Evaluate the trained model on the test set and log performance metrics."""
-    score = model.score(x_test, y_test)
-    _logger.info(f"XGBoost Multi-Target R^2 Score (Testing Set): {score:.4f}")
-    y_pred = model.predict(x_test)
+    """Evaluate the trained model on the test set and log performance metrics.
+
+    Parameters
+    ----------
+    model : XGBRegressor
+        Trained model.
+    x_test : pd.DataFrame
+        Test features.
+    y_pred : pd.DataFrame
+        Predicted targets (already inverse-transformed to original scale).
+    y_test : pd.DataFrame
+        True targets (in original scale).
+    df : pd.DataFrame
+        Full dataset for accessing baseline values.
+    x_cols : list
+        Feature column names.
+    y_data : pd.DataFrame
+        All target data.
+    name : str
+        Model name.
+    shap_per_energy : bool, optional
+        Whether to compute SHAP values per energy bin.
+    """
+    # Compute metrics on original-scale predictions
     mse = mean_squared_error(y_test, y_pred)
-    _logger.info(f"{name} Mean Squared Error (All targets): {mse:.4f}")
+    _logger.info(f"{name} Mean Squared Error (All targets, unscaled): {mse:.4f}")
     mae = mean_absolute_error(y_test, y_pred)
-    _logger.info(f"{name} Mean Absolute Error (All targets): {mae:.4f}")
+    _logger.info(f"{name} Mean Absolute Error (All targets, unscaled): {mae:.4f}")
 
     target_variance(y_test, y_pred, y_data.columns)
     feature_importance(model, x_cols, y_data.columns, name)
@@ -87,9 +105,8 @@ def evaluate_regression_model(
         if shap_per_energy:
             shap_feature_importance_by_energy(model, x_test, df, y_test, y_data.columns)
 
-    df_pred = pd.DataFrame(y_pred, columns=target_features("stereo_analysis"))
     calculate_resolution(
-        df_pred,
+        y_pred,
         y_test,
         df,
         percentiles=[68, 90, 95],
@@ -104,8 +121,9 @@ def target_variance(y_test, y_pred, targets):
     """Calculate and log variance explained per target."""
     y_test_np = y_test.to_numpy() if hasattr(y_test, "to_numpy") else y_test
 
-    mse_values = np.mean((y_test_np - y_pred) ** 2, axis=0)
-    variance_values = np.var(y_test_np, axis=0)
+    # Force numpy arrays so integer indexing is positional and future-proof.
+    mse_values = np.asarray(np.mean((y_test_np - y_pred) ** 2, axis=0))
+    variance_values = np.asarray(np.var(y_test_np, axis=0))
 
     _logger.info("--- Performance Per Target ---")
     for i, name in enumerate(targets):
@@ -127,14 +145,40 @@ def target_variance(y_test, y_pred, targets):
 
 def calculate_resolution(y_pred, y_test, df, percentiles, log_e_min, log_e_max, n_bins, name):
     """Compute angular and energy resolution based on predictions."""
+    # Model predicts residuals, so reconstruct full predictions and MC truth
+    # from residuals and DispBDT baseline
+    _logger.debug(
+        f"Evaluation: y_test indices min={y_test.index.min()}, max={y_test.index.max()}, len={len(y_test)}"
+    )
+    _logger.debug(
+        f"Evaluation: df shape={df.shape}, index min={df.index.min()}, max={df.index.max()}"
+    )
+
+    disp_xoff = df.loc[y_test.index, "Xoff_weighted_bdt"].values
+    disp_yoff = df.loc[y_test.index, "Yoff_weighted_bdt"].values
+
+    # Handle ErecS with proper checks for valid values
+    erec_s = df.loc[y_test.index, "ErecS"].values
+    disp_erec_log = np.where(erec_s > 0, np.log10(erec_s), np.nan)
+
+    # Reconstruct MC truth from residuals in y_test (residual = MC_true - DispBDT)
+    mc_xoff_true = y_test["Xoff_residual"].values + disp_xoff
+    mc_yoff_true = y_test["Yoff_residual"].values + disp_yoff
+    mc_e0_true = y_test["E_residual"].values + disp_erec_log
+
+    # Reconstruct predictions from residual predictions
+    mc_xoff_pred = y_pred["Xoff_residual"].values + disp_xoff
+    mc_yoff_pred = y_pred["Yoff_residual"].values + disp_yoff
+    mc_e0_pred = y_pred["E_residual"].values + disp_erec_log
+
     results_df = pd.DataFrame(
         {
-            "MCxoff_true": y_test["MCxoff"].values,
-            "MCyoff_true": y_test["MCyoff"].values,
-            "MCxoff_pred": y_pred["MCxoff"].values,
-            "MCyoff_pred": y_pred["MCyoff"].values,
-            "MCe0_pred": y_pred["MCe0"].values,
-            "MCe0": df.loc[y_test.index, "MCe0"].values,
+            "MCxoff_true": mc_xoff_true,
+            "MCyoff_true": mc_yoff_true,
+            "MCxoff_pred": mc_xoff_pred,
+            "MCyoff_pred": mc_yoff_pred,
+            "MCe0_pred": mc_e0_pred,
+            "MCe0": mc_e0_true,
         }
     )
 

src/eventdisplay_ml/features.py

@@ -16,7 +16,7 @@ def target_features(analysis_type):
         List of target feature names.
     """
     if analysis_type == "stereo_analysis":
-        return ["MCxoff", "MCyoff", "MCe0"]  # sequence matters
+        return ["Xoff_residual", "Yoff_residual", "E_residual"]  # sequence matters
     if "classification" in analysis_type:
         return []
     raise ValueError(f"Unknown analysis type: {analysis_type}")
@@ -40,6 +40,7 @@ def excluded_features(analysis_type, ntel):
     """
     if analysis_type == "stereo_analysis":
         return {
+            # Pointing corrections applied during preprocessing
             *[f"fpointing_dx_{i}" for i in range(ntel)],
             *[f"fpointing_dy_{i}" for i in range(ntel)],
         }
@@ -130,7 +131,8 @@ def _regression_features(training):
         "Ycore",
     ]
     if training:
-        return [*target_features("stereo_analysis"), *var]
+        # Load MC truth values (residuals will be computed from these)
+        return ["MCxoff", "MCyoff", "MCe0", *var]
     return var
 
 

src/eventdisplay_ml/hyper_parameters.py

@@ -10,10 +10,13 @@
     "xgboost": {
         "model": None,
         "hyper_parameters": {
-            "n_estimators": 1000,
-            "learning_rate": 0.1,  # Shrinkage
-            "max_depth": 10,
+            "n_estimators": 2000,
+            "early_stopping_rounds": 50,
+            "learning_rate": 0.05,  # Shrinkage
+            "max_depth": 6,
             "min_child_weight": 5.0,  # Equivalent to MinNodeSize=1.0% for XGBoost
+            "gamma": 0.5,
+            "reg_lambda": 2.0,
             "objective": "reg:squarederror",
             "n_jobs": 8,
             "random_state": None,