From cfd59222abb892ad6b7850c5cd41de332d22488d Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Wed, 13 Aug 2025 11:29:57 +0100 Subject: [PATCH 1/8] Modifications to unit tests. --- .../quantile_regression_random_forest.py | 31 +++--- ...train_quantile_regression_random_forest.py | 2 +- ...apply_quantile_regression_random_forest.py | 18 ++-- .../test_quantile_regression_random_forest.py | 95 +++++++++---------- 4 files changed, 76 insertions(+), 70 deletions(-) diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index 482743e339..21be24ccf1 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -360,6 +360,7 @@ def __init__( self.compression = compression self.output = model_output self.kwargs = kwargs + self.expected_coordinate_order = ["forecast_reference_time", "forecast_period"] def fit_qrf( self, forecast_features: np.ndarray, target: np.ndarray @@ -384,8 +385,9 @@ def fit_qrf( qrf_model.fit(forecast_features, target) return qrf_model - @staticmethod - def _organise_truth_data(forecast_cube: Cube, truth_cube: Cube) -> list[np.ndarray]: + def _organise_truth_data( + self, forecast_cube: Cube, truth_cube: Cube + ) -> list[np.ndarray]: """Organise the truth data, so that the validity time matches the validity time of the forecast. This might mean that the truth data is repeated, if, for example, the forecast has multiple forecast reference times and multiple @@ -425,7 +427,7 @@ def _organise_truth_data(forecast_cube: Cube, truth_cube: Cube) -> list[np.ndarr # forecast period and forecast reference time. enforce_coordinate_ordering( forecast_cube, - ["forecast_period", "forecast_reference_time"], + self.expected_coordinate_order[::-1], anchor_start=True, ) for fp in fp_coord.points: @@ -495,15 +497,20 @@ def process( ) # Ensure the forecast cube has the correct dimension ordering for prep_feature. - - coord_dims = [forecast_cube.coord_dims(c) for c in ["forecast_reference_time", "forecast_period"]] - coord_names = [forecast_cube.coord(dimensions=d, dim_coords=True).name() for d in coord_dims if len(d) > 0] - print("im here") + + coord_dims = [ + forecast_cube.coord_dims(c) for c in self.expected_coordinate_order + ] + coord_names = [ + forecast_cube.coord(dimensions=d, dim_coords=True).name() + for d in coord_dims + if len(d) > 0 + ] enforce_coordinate_ordering( - forecast_cube, - [coord_names], - anchor_start=True, - ) + forecast_cube, + [coord_names], + anchor_start=True, + ) feature_values = [] @@ -521,9 +528,9 @@ def process( prep_feature(forecast_cube, feature_cube[0], feature) ) feature_values = np.array(feature_values).T - truth_data_list = self._organise_truth_data(forecast_cube, truth_cube) target_values = np.array(truth_data_list).flatten() + # Fit the quantile regression model qrf_model = self.fit_qrf(feature_values, target_values) diff --git a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py index 009852f4d3..091a30a0b3 100644 --- a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py +++ b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py @@ -59,7 +59,7 @@ def test_basic( "42", "--compression", "5", - "--model-output", + "--output", output_path, ] if transformation == "with_transformation": diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py index 906cd41332..c8fb8bff6f 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py @@ -27,16 +27,16 @@ @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,quantiles,expected", [ - (2, 2, 55, 5, None, 0, {}, False, [0.5], [5.15, 5.65]), # noqa Basic test case - (100, 2, 55, 5, None, 0, {}, False, [1 / 3, 2 / 3], [[4.1, 5.1], [4.2, 5.1]]), # noqa Multiple quantiles - (1, 1, 55, 5, None, 0, {}, False, [0.5], [6.2, 6.2]), # noqa Fewer estimators and reduced depth + (2, 2, 55, 5, None, 0, {}, False, [0.5], [4.1, 5.65]), # noqa Basic test case + (100, 2, 55, 5, None, 0, {}, False, [1 / 3, 2 / 3], [[4.1, 5.1], [5.1, 5.1]]), # noqa Multiple quantiles + (1, 1, 55, 5, None, 0, {}, False, [0.5], [4.1, 6.2]), # noqa Fewer estimators and reduced depth (1, 1, 73, 5, None, 0, {}, False, [0.5], [4.2, 6.2]), # Different random state - (2, 2, 55, 5, "log", 10, {}, False, [0.5], [5.11, 5.64]), # Log transformation - (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [5.11, 5.64]), # noqa Log10 transformation - (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [5.11, 5.64]), # noqa Square root transformation - (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [5.13, 5.64]), # noqa Cube root transformation - (2, 2, 55, 5, None, 0, {"max_samples_leaf": 0.5}, False, [0.5], [5.15, 6.2]), # noqa # Different criterion - (2, 5, 55, 5, None, 0, {}, True, [0.5], [5.15, 5.65]), # noqa Include an additional static feature + (2, 2, 55, 5, "log", 10, {}, False, [0.5], [4.1, 5.1]), # Log transformation + (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Log10 transformation + (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Square root transformation + (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Cube root transformation + (2, 2, 55, 5, None, 0, {"max_samples_leaf": 0.5}, False, [0.5], [4.1, 6.2]), # noqa # Different criterion + (2, 5, 55, 5, None, 0, {}, True, [0.5], [4.1, 5.65]), # noqa Include an additional static feature ], ) def test_load_and_apply_qrf( diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py index b7ff4e6825..35591920df 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py @@ -45,8 +45,7 @@ def _create_forecasts( Returns: Forecast cube containing three percentiles and two sites. """ - - data = np.array(data, dtype=np.float32).repeat(2).reshape(len(data), 2) + data = np.array([data, data + 2], dtype=np.float32).T cube = set_up_spot_variable_cube( data, realizations=range(len(data)), @@ -204,7 +203,7 @@ def _run_train_qrf( @pytest.mark.parametrize( "feature,expected,expected_dtype", [ - ("mean", np.array([6, 6], dtype=np.float32), np.float32), + ("mean", np.array([6, 8], dtype=np.float32), np.float32), ("std", np.array([4, 4], dtype=np.float32), np.float32), ("latitude", np.array([50, 60], dtype=np.float32), np.float32), ("longitude", np.array([0, 10], dtype=np.float32), np.float32), @@ -236,7 +235,7 @@ def test_prep_feature_single_time(feature, expected, expected_dtype): """Test the prep_feature function for a single time.""" forecast_reference_time = "20170101T0000Z" validity_time = "20170101T1200Z" - data = [2, 6, 10] + data = np.array([2, 6, 10]) day_of_training_period = [0] forecast_cube = _create_forecasts(forecast_reference_time, validity_time, data) forecast_cube = _add_day_of_training_period( @@ -257,7 +256,7 @@ def test_prep_feature_single_time(feature, expected, expected_dtype): @pytest.mark.parametrize( "feature,expected,expected_dtype", [ - ("mean", np.repeat(6, 8).astype(np.float32), np.float32), + ("mean", np.tile([6, 8], 4).astype(np.float32), np.float32), ("std", np.repeat(4, 8).astype(np.float32), np.float32), ("latitude", np.tile([50, 60], 4).astype(np.float32), np.float32), ("longitude", np.tile([0, 10], 4).astype(np.float32), np.float32), @@ -304,7 +303,7 @@ def test_prep_feature_1d_time_dimension(feature, expected, expected_dtype): "20170102T1200Z", ] - data = [2, 6, 10] + data = np.array([2, 6, 10]) forecast_cubes = CubeList() for frt, vt in zip(forecast_reference_times, validity_times): @@ -330,7 +329,7 @@ def test_prep_feature_1d_time_dimension(feature, expected, expected_dtype): @pytest.mark.parametrize( "feature,expected,expected_dtype", [ - ("mean", np.repeat(6, 12).astype(np.float32), np.float32), + ("mean", np.tile([6, 8], 6).astype(np.float32), np.float32), ("std", np.repeat(4, 12).astype(np.float32), np.float32), ("latitude", np.tile([50, 60], 6).astype(np.float32), np.float32), ("longitude", np.tile([0, 10], 6).astype(np.float32), np.float32), @@ -381,7 +380,7 @@ def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype): "20170103T1200Z", ] - data = [2, 6, 10] + data = np.array([2, 6, 10]) forecast_cubes = CubeList() for frt, vt in zip(forecast_reference_times, validity_times): forecast_cubes.append(_create_forecasts(frt, vt, data)) @@ -409,15 +408,15 @@ def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype): @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected", [ - (2, 2, 55, 5, None, 0, {}, False, 4.15), # Basic test case - (2, 2, 54, 5, None, 0, {}, False, 4.2), # Different random state - (1, 1, 55, 5, None, 0, {}, False, 4.2), # Fewer estimators and reduced depth + (2, 2, 55, 5, None, 0, {}, False, 4.1), # Basic test case + (2, 2, 54, 5, None, 0, {}, False, 4.15), # Different random state + (1, 1, 55, 5, None, 0, {}, False, 4.1), # Fewer estimators and reduced depth (2, 2, 55, 5, "log", 10, {}, False, 2.65), # Log transformation (2, 2, 55, 5, "log10", 10, {}, False, 1.15), # Log10 transformation (2, 2, 55, 5, "sqrt", 10, {}, False, 3.76), # Square root transformation (2, 2, 55, 5, "cbrt", 10, {}, False, 2.42), # Cube root transformation - (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.15), # noqa # Different criterion - (2, 5, 55, 5, None, 0, {}, True, 4.15), # Include static data + (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.1), # noqa # Different criterion + (2, 5, 55, 5, None, 0, {}, True, 4.1), # Include static data ], ) def test_train_qrf_single_lead_times( @@ -478,15 +477,15 @@ def test_train_qrf_single_lead_times( @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected", [ - (2, 2, 55, 5, None, 0, {}, False, 4.0), # Basic test case + (2, 2, 55, 5, None, 0, {}, False, 5.6), # Basic test case (1, 1, 55, 5, None, 0, {}, False, 3.8), # Fewer estimators and reduced depth - (1, 1, 73, 5, None, 0, {}, False, 7), # Different random state - (2, 2, 55, 5, "log", 10, {}, False, 2.743), # Log transformation - (2, 2, 55, 5, "log10", 10, {}, False, 1.191), # Log10 transformation - (2, 2, 55, 5, "sqrt", 10, {}, False, 3.946), # Square root transformation - (2, 2, 55, 5, "cbrt", 10, {}, False, 2.496), # Cube root transformation + (1, 1, 73, 5, None, 0, {}, False, 7.0), # Different random state + (2, 2, 55, 5, "log", 10, {}, False, 2.713), # Log transformation + (2, 2, 55, 5, "log10", 10, {}, False, 1.178), # Log10 transformation + (2, 2, 55, 5, "sqrt", 10, {}, False, 3.884), # Square root transformation + (2, 2, 55, 5, "cbrt", 10, {}, False, 2.471), # Cube root transformation (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4), # noqa # Different criterion - (2, 5, 55, 5, None, 0, {}, True, 4), # Include static data + (2, 5, 55, 5, None, 0, {}, True, 5), # Include static data ], ) def test_train_qrf_multiple_lead_times( @@ -536,14 +535,14 @@ def test_train_qrf_multiple_lead_times( @pytest.mark.parametrize( "feature_config,data,include_static,expected", [ - ({"wind_speed_at_10m": ["mean"]}, [5], False, [4]), # One feature + ({"wind_speed_at_10m": ["mean"]}, [5], False, [5]), # One feature ({"wind_speed_at_10m": ["latitude"]}, [55], False, [6.65]), # noqa Without the target - ({"wind_speed_at_10m": ["mean"]}, [5], True, [4]), # With static data + ({"wind_speed_at_10m": ["mean"]}, [5], True, [5]), # With static data ( {"wind_speed_at_10m": ["mean"], "air_temperature": ["mean"]}, [5], False, - [4], + [5], ), # Multiple dynamic features ( {"wind_speed_at_10m": ["mean"], "pressure_at_mean_sea_level": ["mean"]}, @@ -619,31 +618,31 @@ def test_alternative_feature_configs( @pytest.mark.parametrize( "quantiles,transformation,pre_transform_addition,include_static,expected", [ - ([0.5], None, 0, False, [4, 4]), # One quantile - # ([0.1, 0.5, 0.9], None, 0, False, [[9.1, 9.14], [9.1, 9.3], [9.1, 9.46]]), # noqa Multiple quantiles - # ([0.1, 0.5, 0.9], "log", 10, False, [[4.46, 4.46], [5.54, 5.54], [6.7, 6.7]]), # noqa Log transformation - # ( - # [0.1, 0.5, 0.9], - # "log10", - # 10, - # False, - # [[4.46, 4.46], [5.54, 5.54], [6.7, 6.7]], - # ), # Log10 transformation - # ( - # [0.1, 0.5, 0.9], - # "sqrt", - # 10, - # False, - # [[4.47, 4.47], [5.57, 5.57], [6.71, 6.71]], - # ), # Square root transformation - # ( - # [0.1, 0.5, 0.9], - # "cbrt", - # 10, - # False, - # [[4.47, 4.47], [5.56, 5.56], [6.7, 6.7]], - # ), # Cube root transformation - # ([0.1, 0.5, 0.9], None, 0, True, [[9.1, 9.14], [9.1, 9.3], [9.1, 9.46]]), # noqa Include static data + ([0.5], None, 0, False, [5.6, 5.4]), # One quantile + ([0.1, 0.5, 0.9], None, 0, False, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Multiple quantiles + ([0.1, 0.5, 0.9], "log", 10, False, [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]]), # noqa Log transformation + ( + [0.1, 0.5, 0.9], + "log10", + 10, + False, + [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]], + ), # Log10 transformation + ( + [0.1, 0.5, 0.9], + "sqrt", + 10, + False, + [[4.38, 4.38], [5.09, 5.09], [5.82, 5.82]], + ), # Square root transformation + ( + [0.1, 0.5, 0.9], + "cbrt", + 10, + False, + [[4.37, 4.37], [5.08, 5.08], [5.81, 5.81]], + ), # Cube root transformation + ([0.1, 0.5, 0.9], None, 0, True, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Include static data ], ) def test_apply_qrf( From 58e203dc55ff4d41c204780660702576eda2a65d Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Thu, 14 Aug 2025 14:08:33 +0100 Subject: [PATCH 2/8] Minor edits including edit to avoid dataframe warning. --- improver/calibration/dataframe_utilities.py | 12 +++--------- ..._and_train_quantile_regression_random_forest.py | 8 ++++---- .../quantile_regression_random_forest.py | 14 +++++++++----- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/improver/calibration/dataframe_utilities.py b/improver/calibration/dataframe_utilities.py index bc103eb983..857442390b 100644 --- a/improver/calibration/dataframe_utilities.py +++ b/improver/calibration/dataframe_utilities.py @@ -76,8 +76,7 @@ def _dataframe_column_check(df: DataFrame, compulsory_columns: Sequence) -> None if not set(compulsory_columns).issubset(df.columns): diff = set(compulsory_columns).difference(df.columns) msg = ( - "The following compulsory column(s) are missing from the " - f"DataFrame: {diff}" + f"The following compulsory column(s) are missing from the DataFrame: {diff}" ) raise ValueError(msg) @@ -183,11 +182,7 @@ def _fill_missing_entries(df, combi_cols, static_cols, site_id_col): df = df.set_index(combi_cols).reindex(new_index).reset_index(level=combi_cols) # Fill the NaNs within the static columns for each wmo_id. - filled_df = ( - df.groupby(site_id_col)[combi_cols + static_cols] - .fillna(method="ffill") - .fillna(method="bfill") - ) + filled_df = df.groupby(site_id_col)[combi_cols + static_cols].ffill().bfill() df = df.drop(columns=static_cols) df = df.merge(filled_df, on=combi_cols) @@ -361,8 +356,7 @@ def get_forecast_representation(df: DataFrame) -> str: ) if len(representations) == 0: raise ValueError( - f"None of the columns {REPRESENTATION_COLUMNS} " - "exist in the input dataset" + f"None of the columns {REPRESENTATION_COLUMNS} exist in the input dataset" ) return representations.pop() diff --git a/improver/calibration/load_and_train_quantile_regression_random_forest.py b/improver/calibration/load_and_train_quantile_regression_random_forest.py index 160486bd3e..35495606d7 100644 --- a/improver/calibration/load_and_train_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_train_quantile_regression_random_forest.py @@ -283,12 +283,12 @@ def _dataframe_to_cubes( else: msg = "Concatenating the forecast has failed to create a single cube." raise ValueError(msg) - + # Promote the forecast_reference_time coord to a dimension coordinate if the forecast_period is one. - if forecast_cube.coord_dims("forecast_period") and not forecast_cube.coord_dims("forecast_reference_time"): + if forecast_cube.coord_dims("forecast_period") and not forecast_cube.coord_dims( + "forecast_reference_time" + ): forecast_cube = iris.util.new_axis(forecast_cube, "forecast_reference_time") - if not forecast_cube.coords("forecast_period", dim_coords=True) and forecast_cube.coord_dims("forecast_period"): - iris.util.promote_aux_coord_to_dim_coord(forecast_cube, "forecast_period") return forecast_cube, truth_cube diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index 21be24ccf1..8a6b9fc465 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -501,11 +501,15 @@ def process( coord_dims = [ forecast_cube.coord_dims(c) for c in self.expected_coordinate_order ] - coord_names = [ - forecast_cube.coord(dimensions=d, dim_coords=True).name() - for d in coord_dims - if len(d) > 0 - ] + coord_names = [] + for dim in coord_dims: + if len(dim) == 0: + continue + if len(forecast_cube.coords(dimensions=dim, dim_coords=True)) == 0: + continue + coord_names.append( + forecast_cube.coord(dimensions=dim, dim_coords=True).name() + ) enforce_coordinate_ordering( forecast_cube, [coord_names], From dd3605800624a162f0fd8398dcc4e4f4a8703930 Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Fri, 15 Aug 2025 22:14:46 +0100 Subject: [PATCH 3/8] Re-write to use dataframes more widely throughout the QRF implementation to avoid cost and complexity of converting to cubes. --- improver/calibration/__init__.py | 4 +- ...apply_quantile_regression_random_forest.py | 79 +-- ...train_quantile_regression_random_forest.py | 165 ++--- .../quantile_regression_random_forest.py | 582 +++++------------- improver/utilities/compare.py | 22 +- improver_tests/acceptance/SHA256SUMS | 8 +- ...apply_quantile_regression_random_forest.py | 55 +- ...train_quantile_regression_random_forest.py | 6 +- .../test_quantile_regression_random_forest.py | 540 ++++++++-------- 9 files changed, 598 insertions(+), 863 deletions(-) diff --git a/improver/calibration/__init__.py b/improver/calibration/__init__.py index 6b288ea29a..723395a9ae 100644 --- a/improver/calibration/__init__.py +++ b/improver/calibration/__init__.py @@ -24,7 +24,7 @@ ("altitude", pa.float32()), ("blend_time", pa.timestamp("s", "utc")), ("forecast_period", pa.int64()), - ("forecast_reference_time", pa.int64()), + ("forecast_reference_time", pa.timestamp("s", "utc")), ("latitude", pa.float32()), ("longitude", pa.float32()), ("time", pa.timestamp("s", "utc")), @@ -44,7 +44,7 @@ ("latitude", pa.float32()), ("longitude", pa.float32()), ("altitude", pa.float32()), - ("time", pa.int64()), + ("time", pa.timestamp("s", "utc")), ("wmo_id", pa.string()), ("ob_value", pa.float32()), ] diff --git a/improver/calibration/load_and_apply_quantile_regression_random_forest.py b/improver/calibration/load_and_apply_quantile_regression_random_forest.py index c28d9896e6..cdfaeae530 100644 --- a/improver/calibration/load_and_apply_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_apply_quantile_regression_random_forest.py @@ -13,6 +13,7 @@ import joblib import numpy as np from iris.cube import Cube, CubeList +from iris.pandas import as_data_frame from quantile_forest import RandomForestQuantileRegressor from improver import PostProcessingPlugin @@ -27,6 +28,8 @@ from improver.ensemble_copula_coupling.utilities import choose_set_of_percentiles from improver.utilities.cube_checker import assert_spatial_coords_match +iris.FUTURE.pandas_ndim = True + class LoadAndApplyQRF(PostProcessingPlugin): """Load and apply the trained Quantile Regression Random Forest (QRF) model.""" @@ -188,44 +191,6 @@ def _percentiles_to_realizations(cube_inputs: CubeList) -> CubeList: realization_cube_inputs.append(feature_cube) return realization_cube_inputs - @staticmethod - def _organise_cubes( - cube_inputs: CubeList, forecast_cube: Cube - ) -> tuple[CubeList, Cube]: - """Promote the forecast period and forecast reference time coordinates to be - dimension coordinates, if present, on the feature cubes and the template - forecast cube. - - Args: - cube_inputs: CubeList of feature cubes, which may include the forecast to be - forecast_cube: Forecast cube for use as a template. - - Returns: - Feature cubes and template cube with forecast period and - forecast reference time promoted to dimension coordinates. - """ - # Ensure that forecast_period is a dimension on all cubes. - fp_dim_cube_inputs = iris.cube.CubeList([]) - for feature_cube in cube_inputs: - if feature_cube.coords("forecast_period", dim_coords=False): - feature_cube = iris.util.new_axis(feature_cube, "forecast_period") - if feature_cube.coords("forecast_reference_time", dim_coords=False): - feature_cube = iris.util.new_axis( - feature_cube, "forecast_reference_time" - ) - fp_dim_cube_inputs.append(feature_cube) - cube_inputs = fp_dim_cube_inputs - - # Ensure the forecast cube has the same dimensions as the features - template_forecast_cube = iris.util.new_axis(forecast_cube, "forecast_period") - template_forecast_cube = iris.util.new_axis( - template_forecast_cube, "forecast_reference_time" - ) - - # Check that the grids are the same for all dynamic predictors and the forecast - assert_spatial_coords_match(cube_inputs) - return cube_inputs, template_forecast_cube - def process( self, file_paths: list[pathlib.Path], @@ -261,17 +226,41 @@ def process( elif forecast_cube.coords("realization"): percentiles = self._compute_percentiles(forecast_cube.copy(), "realization") - cube_inputs, template_forecast_cube = self._organise_cubes( - cube_inputs, template_forecast_cube - ) + assert_spatial_coords_match(cube_inputs) - result = ApplyQuantileRegressionRandomForests( + df = as_data_frame(cube_inputs[0], add_aux_coords=True).reset_index() + + for cube in cube_inputs[1:]: + temporary_df = as_data_frame(cube, add_aux_coords=True).reset_index() + possible_columns = [ + "wmo_id", + "time", + "forecast_reference_time", + "forecast_period", + ] + merge_columns = [ + col for col in possible_columns if col in temporary_df.columns + ] + df = df.merge( + temporary_df[merge_columns + [cube.name()]], + on=merge_columns, + how="left", + ) + + calibrated_forecast = ApplyQuantileRegressionRandomForests( + target_name=self.target_cube_name, feature_config=self.feature_config, quantiles=percentiles, transformation=self.transformation, pre_transform_addition=self.pre_transform_addition, - )(qrf_model, cube_inputs, template_forecast_cube) + )(qrf_model, df) + + calibrated_forecast_cube = template_forecast_cube.copy( + data=np.broadcast_to(calibrated_forecast.T, template_forecast_cube.shape) + ) if forecast_cube.coords("percentile"): - result = RebadgeRealizationsAsPercentiles()(result) - return result + calibrated_forecast_cube = RebadgeRealizationsAsPercentiles()( + calibrated_forecast_cube + ) + return calibrated_forecast_cube diff --git a/improver/calibration/load_and_train_quantile_regression_random_forest.py b/improver/calibration/load_and_train_quantile_regression_random_forest.py index 35495606d7..32085cd0de 100644 --- a/improver/calibration/load_and_train_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_train_quantile_regression_random_forest.py @@ -10,21 +10,20 @@ from typing import Optional import iris -import numpy as np import pandas as pd import pyarrow as pa import pyarrow.parquet as pq +from iris.pandas import as_data_frame from improver import PostProcessingPlugin from improver.calibration import FORECAST_SCHEMA, TRUTH_SCHEMA -from improver.calibration.dataframe_utilities import ( - forecast_and_truth_dataframes_to_cubes, -) from improver.calibration.quantile_regression_random_forest import ( TrainQuantileRegressionRandomForests, ) from improver.utilities.load import load_cube +iris.FUTURE.pandas_ndim = True + class LoadAndTrainQRF(PostProcessingPlugin): """Plugin to load and train a Quantile Regression Random Forest (QRF) model.""" @@ -218,111 +217,68 @@ def _read_parquet_files( raise IOError(msg) return forecast_df, truth_df - def _dataframe_to_cubes( - self, - forecast_df: pd.DataFrame, - truth_df: pd.DataFrame, - forecast_periods: list[int], - ) -> tuple[iris.cube.Cube, iris.cube.Cube]: - """Convert the forecast and truth dataframes to cubes at each forecast period - required. + def _check_matching_times( + self, forecast_df: pd.DataFrame, truth_df: pd.DataFrame + ) -> list[pd.Timestamp]: + return list(set(forecast_df["time"]).intersection(set(truth_df["time"]))) + + def _add_features_to_df( + self, forecast_df: pd.DataFrame, cube_inputs: iris.cube.CubeList + ) -> pd.DataFrame: + """Add features to the forecast DataFrame based on the feature configuration. Args: forecast_df: DataFrame containing the forecast data. - truth_df: DataFrame containing the truth data. - forecast_periods: List of forecast periods in seconds. + cube_inputs: List of cubes containing additional features. Returns: - Tuple containing: - - Cube containing the forecast data. - - Cube containing the truth data. - - Raises: - ValueError: The forecast has failed to concatenate into a single cube. + DataFrame with additional features added. """ - forecast_cubes = iris.cube.CubeList([]) - truth_cubes = iris.cube.CubeList([]) - - for forecast_period in forecast_periods: - forecast_cube, truth_cube = forecast_and_truth_dataframes_to_cubes( - forecast_df, - truth_df, - self.cycletime, - forecast_period, - self.training_length, - experiment=self.experiment, - ) - - if forecast_cube is None or truth_cube is None: - continue - - if not forecast_cube.coords("realization", dim_coords=True): - forecast_cube = iris.util.new_axis(forecast_cube, "realization") - - for forecast_slice in forecast_cube.slices_over("forecast_reference_time"): - forecast_cubes.append(forecast_slice) - - for truth_slice in truth_cube.slices_over("time"): - truth_slice = iris.util.new_axis(truth_slice, "time") - # Multiple forecasts can match to the same observation. This check - # ensures that we do not add the same truth slice multiple times. - if truth_slice not in truth_cubes: - truth_cubes.append(truth_slice) - - if not forecast_cubes or not truth_cubes: - return None, None - - truth_cube = truth_cubes.concatenate_cube() - forecast_cube = forecast_cubes.merge() - - # concatenate_cube() can fail for the forecast_cube, even though calling - # concatenate() results in a single cube. This check ensures the concatenation - # was successful. - if len(forecast_cube) == 1: - forecast_cube = forecast_cube[0] - else: - msg = "Concatenating the forecast has failed to create a single cube." - raise ValueError(msg) - - # Promote the forecast_reference_time coord to a dimension coordinate if the forecast_period is one. - if forecast_cube.coord_dims("forecast_period") and not forecast_cube.coord_dims( - "forecast_reference_time" - ): - forecast_cube = iris.util.new_axis(forecast_cube, "forecast_reference_time") - - return forecast_cube, truth_cube + for feature_name, feature_list in self.feature_config.items(): + for feature in feature_list: + if feature == "static": + # Use the cube's data directly as a feature. + constr = iris.Constraint(name=feature_name) + feature_cube = cube_inputs.extract_cube(constr) + feature_df = as_data_frame(feature_cube, add_aux_coords=True) + forecast_df = forecast_df.merge( + feature_df[["wmo_id", feature_name]], on=["wmo_id"], how="left" + ) + return forecast_df @staticmethod def filter_bad_sites( - forecast_cube: iris.cube.Cube, - truth_cube: iris.cube.Cube, - cube_inputs: iris.cube.CubeList, - ) -> tuple[iris.cube.Cube, iris.cube.Cube, iris.cube.CubeList]: + forecast_df: pd.DataFrame, + truth_df: pd.DataFrame, + ) -> tuple[pd.DataFrame, pd.DataFrame]: """Remove sites that have NaNs in the data. Args: - forecast_cube: Cube containing the forecast data. - truth_cube: Cube containing the truth data. - cube_inputs: List of additional feature cubes. + feature_df: DataFrame containing the forecast data with features. + truth_df: DataFrame containing the truth data. Returns: Tuple containing: - - Cube containing the forecast data with bad sites removed. - - Cube containing the truth data with bad sites removed. - - List of additional feature cubes with bad sites removed. + - DataFrame containing the forecast data with bad sites removed. + - DataFrame containing the truth data with bad sites removed. """ - nan_mask = np.any(np.isnan(truth_cube.data), axis=truth_cube.coord_dims("time")) - all_site_ids = truth_cube.coord("wmo_id").points - bad_site_ids = all_site_ids[nan_mask] - constr = iris.Constraint(wmo_id=lambda cell: cell not in bad_site_ids.tolist()) - truth_cube = truth_cube.extract(constr) - forecast_cube = forecast_cube.extract(constr) - feature_cube_inputs = iris.cube.CubeList([]) - for cube in cube_inputs: - cube = cube.extract(constr) - feature_cube_inputs.append(cube) - - return forecast_cube, truth_cube, feature_cube_inputs + # import pdb + # pdb.set_trace() + # for coord in ["latitude", "longitude", "altitude", "ob_value"]: + # truth_df = truth_df.groupby("wmo_id").filter( + # lambda x: ~(x[coord].isna().any()) + # ) + # import pdb + # pdb.set_trace() + truth_df.dropna( + subset=["latitude", "longitude", "altitude", "ob_value"], inplace=True + ) + + wmo_ids = set(forecast_df["wmo_id"]).intersection(set(truth_df["wmo_id"])) + + forecast_df = forecast_df[forecast_df["wmo_id"].isin(wmo_ids)] + truth_df = truth_df[truth_df["wmo_id"].isin(wmo_ids)] + return forecast_df, truth_df def process( self, @@ -421,23 +377,20 @@ def process( forecast_table_path, truth_table_path, forecast_periods ) - forecast_cube, truth_cube = self._dataframe_to_cubes( - forecast_df, truth_df, forecast_periods + forecast_df = forecast_df[forecast_df["experiment"] == self.experiment] + forecast_df = forecast_df.rename( + columns={"forecast": forecast_df["cf_name"][0]} ) - if forecast_cube is None or truth_cube is None: + # forecast_df = forecast_df.drop(columns=["cf_name", "diagnostic"]) + intersecting_times = self._check_matching_times(forecast_df, truth_df) + if len(intersecting_times) == 0: return None - # If target_forecast is also a dynamic feature in the feature config then - # add it to cube_inputs - for feature_name in self.feature_config.keys(): - if feature_name == forecast_cube[0].name(): - cube_inputs.append(forecast_cube) - - forecast_cube, truth_cube, feature_cube_inputs = self.filter_bad_sites( - forecast_cube, truth_cube, cube_inputs - ) + forecast_df = self._add_features_to_df(forecast_df, cube_inputs) + forecast_df, truth_df = self.filter_bad_sites(forecast_df, truth_df) TrainQuantileRegressionRandomForests( + target_name=forecast_df["cf_name"][0], feature_config=self.feature_config, n_estimators=self.n_estimators, max_depth=self.max_depth, @@ -446,4 +399,4 @@ def process( pre_transform_addition=self.pre_transform_addition, compression=self.compression, model_output=model_output, - )(forecast_cube, truth_cube, feature_cube_inputs) + )(forecast_df, truth_df) diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index 8a6b9fc465..4508d2be97 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -6,281 +6,115 @@ from typing import Optional -import iris import joblib import numpy as np import pandas as pd -from iris.cube import Cube, CubeList +from iris.cube import Cube from quantile_forest import RandomForestQuantileRegressor from improver import BasePlugin, PostProcessingPlugin from improver.constants import DAYS_IN_YEAR, HOURS_IN_DAY -from improver.utilities.cube_manipulation import enforce_coordinate_ordering - - -def _expand_dims( - template_cube: Cube, expansion_array: np.ndarray, features: list[str] -) -> np.ndarray: - """Expand dimensions of the expansion_array to match the dimensions of the template - cube. For example, the dimensions of the expansion_array may only span a portion - of the dimensions present in the template cube. If the expansion_array has - singleton dimensions that are not present in the template cube, these singleton - dimension will be squeezed away. - - Args: - template_cube: Cube that acts as a template for the output. The dimension - coordinates and coordinates associated with the dimension coordinates on - this cube will be used to expand the expansion_array. - expansion_array: Array that will be expanded to match the dimensions of the - template_cube. - features: List of feature names that are present in the template_cube. These - may correspond to dimension coordinates or coordinates associated with - dimension coordinates or scalar coordinates. - - Returns: - Array with dimensions expanded to match the template_cube. - """ - # The aim here is to find the names of the dimension coordinates but check - # whether the name of any coordinate associated with a dimension coordinate matches - # the features provided. If there is a match, then use the feature name, rather than - # the dimension coordinate name. - dim_coord_names = [c.name() for c in template_cube.coords(dim_coords=True)] - refined_dim_coord_names = [] - for coord in template_cube.coords(): - # Ignore scalar coordinates. - if len(template_cube.coord_dims(coord)) > 0: - associated_coords = [ - c.name() - for c in template_cube.coords( - dimensions=template_cube.coord_dims(coord) - ) - ] - feature_associated_coords = list( - set(associated_coords).intersection(features) - ) - if feature_associated_coords: - refined_dim_coord_names.append(feature_associated_coords[0]) - elif coord.name() in dim_coord_names: - refined_dim_coord_names.append(coord.name()) - refined_dim_coord_names = list(set(refined_dim_coord_names)) - - expansion_dim_names = list(set(refined_dim_coord_names) - set(features)) - - # If the expansion_dims and refined_dim_coord_names match, this implies that - # the features provided are not dimension coordinates nor associated with any - # dimension coordinate. In this case, any singleton dimensions on the - # expansion_array should be removed, as it isn't relevant for defining the shape - # of the template cube. - if expansion_dim_names == refined_dim_coord_names: - expansion_array = np.squeeze(expansion_array) - - dims = [] - for dim_name in expansion_dim_names: - dims.extend(template_cube.coord_dims(dim_name)) - - return np.expand_dims(expansion_array, dims) def prep_feature( - template_cube: Cube, - feature_cube: Cube, + df: pd.DataFrame, + feature_name: str, feature: str, -) -> np.ndarray: - """Prepare the feature values for the quantile regression random forest model. - Each expected feature is handled separately. The key steps are: - - Collapse the template cube over the realization dimension to get a shape - template. - - Extract the feature values from the feature cube provided. - - Handle different types of features separately (mean, std, spatial coordinates, - forecast period, model weights, day of year, hour of day, and static features) - to ensure they are broadcasted correctly to match the template cube shape. - Time coordinates are the most complex to handle as they can be scalar, span a - single dimension or span multiple dimensions. - - Flatten the feature values to create a 1D array. - - Set the output dtype based on the feature type. - - Return the flattened array of feature values. - - Args: - template_cube (cube): - The forecast cube that acts only as a template. - feature_cube (cube): - Feature cube that is either static or dynamic. - feature: - The feature to be extracted from the associated feature_cube. - If "static" then the cube itself acts as the feature. - Returns: - feature_values (numpy.ndarray): - Flattened array of feature values. - """ - # Collapse the template cube to get a shape template. - collapsed_cube = template_cube.collapsed(["realization"], iris.analysis.MEAN) - - # Handle each expected feature type separately. - if "mean" == feature: - feature_values = feature_cube.collapsed( - ["realization"], iris.analysis.MEAN - ).data.flatten() - elif "std" == feature: - feature_values = feature_cube.collapsed( - ["realization"], iris.analysis.STD_DEV - ).data.flatten() - elif feature in ["latitude", "longitude", "altitude"]: - coord_multidim = _expand_dims( - collapsed_cube, feature_cube.coord(feature).points, ["spot_index"] +) -> pd.DataFrame: + if feature in ["mean", "std"]: + representation_name = [ + n for n in ["percentile", "realization"] if n in df.columns + ][0] + subset_cols = [ + "forecast_reference_time", + "forecast_period", + "wmo_id", + representation_name, + feature_name, + ] + groupby_cols = ["forecast_reference_time", "forecast_period", "wmo_id"] + if feature == "mean": + subset_df = df[subset_cols].groupby(groupby_cols).mean() + elif feature == "std": + subset_df = df[subset_cols].groupby(groupby_cols).std() + + subset_df = subset_df.reset_index() + subset_df.rename( + columns={feature_name: f"{feature_name}_{feature}"}, inplace=True + ) + # df.drop(columns=[representation_name], inplace=True) + # df.drop_duplicates(inplace=True) + df = df.merge( + subset_df[groupby_cols + [f"{feature_name}_{feature}"]], + on=["forecast_reference_time", "forecast_period", "wmo_id"], + how="left", ) - feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten() - elif feature == "forecast_period": - if len(feature_cube.coord_dims("forecast_period")) == 1: - coord_multidim = _expand_dims( - collapsed_cube, - feature_cube.coord("forecast_period").points, - ["forecast_period"], - ) - else: - coord_multidim = feature_cube.coord("forecast_period").points - - feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten() - elif feature == "model_weights": - if len(feature_cube.coord_dims("forecast_period")) == 1: - coord_multidim = _expand_dims( - collapsed_cube, - feature_cube.coord("model_weights").points, - ["forecast_period"], - ) - else: - coord_multidim = feature_cube.coord("model_weights").points - feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten() elif feature in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]: - frt_dims = feature_cube.coord_dims("forecast_reference_time") - fp_dims = feature_cube.coord_dims("forecast_period") - frt_coord = feature_cube.coord("forecast_reference_time") - fp_coord = feature_cube.coord("forecast_period") - - if len(frt_dims) == 0 and len(fp_dims) == 0: - time_point = frt_coord.cell(0).point._to_real_datetime() + pd.to_timedelta( - fp_coord.points[0], unit=str(fp_coord.units) - ) - day_of_year = np.array([np.int32(time_point.strftime("%j"))]) - elif frt_dims == fp_dims: - # Forecast reference time and forecast period share a dimension coordinate. - # Forecast reference time and forecast period must be mixed together - # along one dimension. - day_of_year = [] - for frt_cell, fp_point in zip(frt_coord.cells(), fp_coord.points): - time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta( - fp_point, unit=str(fp_coord.units) - ) - day_of_year.append(np.int32(time_point.strftime("%j"))) - day_of_year = np.array(day_of_year) - # frt_dims and fp_dims are the same, so we can choose one of them to pop. - day_of_year = _expand_dims( - collapsed_cube, day_of_year, ["forecast_reference_time"] - ) - else: - # Forecast reference time and forecast period are different dimensions. - day_of_year = np.zeros((len(frt_coord.points), len(fp_coord.points))) - for i, frt_cell in enumerate(frt_coord.cells()): - for j, fp_point in enumerate(fp_coord.points): - time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta( - fp_point, unit=str(fp_coord.units) - ) - day_of_year[i, j] = time_point.strftime("%j") - - day_of_year = _expand_dims( - collapsed_cube, - day_of_year, - ["forecast_reference_time", "forecast_period"], - ) - - feature_values = np.broadcast_to(day_of_year, collapsed_cube.shape).flatten() - if feature == "day_of_year_sin": - feature_values = np.sin(2 * np.pi * feature_values / (DAYS_IN_YEAR + 1)) + day_of_year = df["time"].dt.strftime("%j") + day_of_year = np.array(day_of_year, dtype=np.int32) + if feature == "day_of_year": + feature_values = day_of_year + elif feature == "day_of_year_sin": + feature_values = np.sin( + 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) + ).astype(np.float32) elif feature == "day_of_year_cos": - feature_values = np.cos(2 * np.pi * feature_values / (DAYS_IN_YEAR + 1)) + feature_values = np.cos( + 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) + ).astype(np.float32) + df[feature] = feature_values elif feature in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]: - frt_dims = feature_cube.coord_dims("forecast_reference_time") - fp_dims = feature_cube.coord_dims("forecast_period") - frt_coord = feature_cube.coord("forecast_reference_time") - fp_coord = feature_cube.coord("forecast_period") - - if len(frt_dims) == 0 and len(fp_dims) == 0: - time_point = frt_coord.cell(0).point._to_real_datetime() + pd.to_timedelta( - fp_coord.points[0], unit=str(fp_coord.units) - ) - hour_of_day = np.array([np.int32(time_point.hour)]) - elif frt_dims == fp_dims: - # Forecast reference time and forecast period share a dimension coordinate. - # Forecast reference time and forecast period must be mixed together - # along one dimension. - hour_of_day = [] - for frt_cell, fp_point in zip(frt_coord.cells(), fp_coord.points): - time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta( - fp_point, unit=str(fp_coord.units) - ) - hour_of_day.append(np.int32(time_point.hour)) - hour_of_day = np.array(hour_of_day) - - # frt_dims and fp_dims are the same, so we can choose one of them to pop. - hour_of_day = _expand_dims( - collapsed_cube, hour_of_day, ["forecast_reference_time"] - ) - else: - # Forecast reference time and forecast period are different dimensions. - hour_of_day = np.zeros((len(frt_coord.points), len(fp_coord.points))) - for i, frt_cell in enumerate(frt_coord.cells()): - for j, fp_point in enumerate(fp_coord.points): - time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta( - fp_point, unit=str(fp_coord.units) - ) - hour_of_day[i, j] = time_point.hour - - hour_of_day = _expand_dims( - collapsed_cube, - hour_of_day, - ["forecast_reference_time", "forecast_period"], + hour_of_day = df["time"].dt.hour + hour_of_day = np.array(hour_of_day, dtype=np.int32) + if feature == "hour_of_day": + feature_values = hour_of_day + elif feature == "hour_of_day_sin": + feature_values = np.sin(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype( + np.float32 ) - - feature_values = np.broadcast_to(hour_of_day, collapsed_cube.shape).flatten() - if feature == "hour_of_day_sin": - feature_values = np.sin(2 * np.pi * feature_values / HOURS_IN_DAY) elif feature == "hour_of_day_cos": - feature_values = np.cos(2 * np.pi * feature_values / HOURS_IN_DAY) - elif feature == "day_of_training_period": - if len(feature_cube.coord_dims("day_of_training_period")) == 1: - coord_multidim = _expand_dims( - collapsed_cube, - feature_cube.coord("day_of_training_period").points, - ["day_of_training_period"], + feature_values = np.cos(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype( + np.float32 ) - else: - coord_multidim = feature_cube.coord("day_of_training_period").points - - feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten() - elif feature == "static": - coord_multidim = _expand_dims(collapsed_cube, feature_cube.data, ["spot_index"]) - feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten() - else: - # If the feature is not one of the expected types, raise an error. - msg = f"Feature {feature} is not supported." - raise ValueError(msg) + df[feature] = feature_values + return df + + +def sanitise_forecast_dataframe(df: pd.DataFrame, feature_config: dict) -> pd.DataFrame: + representation_name = [n for n in ["percentile", "realization"] if n in df.columns][ + 0 + ] + collapsed_features = [] + for key, values in feature_config.items(): + collapsed_features.extend([key for v in values if v in ["mean", "std"]]) + collapsed_features = list(set(collapsed_features)) + + df = df.drop(columns=[representation_name, *collapsed_features]).drop_duplicates() + return df - # Ensure the feature values are returned with the correct dtype. - if feature in ["mean", "std", "static"]: - feature_values = feature_values.astype(feature_cube.dtype) - elif feature in [ - "day_of_year", - "day_of_year_sin", - "day_of_year_cos", - "hour_of_day", - "hour_of_day_sin", - "hour_of_day_cos", - ]: - feature_values = feature_values.astype(np.float32) - else: - feature_values = feature_values.astype(feature_cube.coord(feature).dtype) - return feature_values +def select_features(df, feature_config): + """Prepare the features from a DataFrame for quantile regression random forest.""" + + # Remove columns that are no longer required e.g. as the mean or std have been computed. + # Removing these columns allows a lot of duplicates to then be removed. + + feature_column_names = [] + for feature_name in feature_config.keys(): + for feature in feature_config[feature_name]: + if feature in ["mean", "std"]: + feature_column_names.append(f"{feature_name}_{feature}") + elif feature in ["static"]: + feature_column_names.append(feature_name) + else: + feature_column_names.append(feature) + + if len(list(set(feature_column_names) - set(df.columns))) > 0: + msg = f"Feature '{feature}' is not supported." + raise ValueError(msg) + + return feature_column_names def _check_valid_transformation(transformation: str): @@ -303,6 +137,7 @@ class TrainQuantileRegressionRandomForests(BasePlugin): def __init__( self, + target_name: str, feature_config: dict[str, list[str]], n_estimators: int, max_depth: Optional[int] = None, @@ -315,17 +150,18 @@ def __init__( ) -> None: """Initialise the plugin. Args: + target_name (str): + Name of the target variable to be calibrated e.g. 'air_temperature'. feature_config (dict): Feature configuration defining the features to be used for quantile regression. The configuration is a dictionary of strings, where the - keys are the names of the input cube(s) supplied, and the values are - a list. This list can contain both computed features, such as the - mean or standard deviation (std), or static features, such as the - altitude. The computed features will be computed using the cube defined - in the dictionary key. If the key is the feature itself e.g. a distance - to water cube, then the value should state "static". This will ensure - the cube's data is used as the feature. The config will have the - structure: + keys are the names of the columns within the dataframe. Some + features may be used as initially provided within the dataframe, + whilst others may be computed from the data e.g. mean, std. + If the key is the feature itself e.g. distance to water, then the value + should state "static". In this case, the name of feature e.g. + 'distance_to_water' is expected to be a column name in the input + dataframe. The config will have the structure: "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], @@ -349,7 +185,7 @@ def __init__( kwargs: Additional keyword arguments for the quantile regression model. """ - + self.target_name = target_name self.feature_config = feature_config self.n_estimators = n_estimators self.max_depth = max_depth @@ -385,91 +221,18 @@ def fit_qrf( qrf_model.fit(forecast_features, target) return qrf_model - def _organise_truth_data( - self, forecast_cube: Cube, truth_cube: Cube - ) -> list[np.ndarray]: - """Organise the truth data, so that the validity time matches the validity - time of the forecast. This might mean that the truth data is repeated, if, for - example, the forecast has multiple forecast reference times and multiple - forecast periods that have the same validity time.""" - frt_dims = forecast_cube.coord_dims("forecast_reference_time") - fp_dims = forecast_cube.coord_dims("forecast_period") - - frt_coord = forecast_cube.coord("forecast_reference_time") - fp_coord = forecast_cube.coord("forecast_period") - - time_datetimes = [] - - # Handle the case where both forecast reference time and forecast period are - # non-dimensional coordinates (i.e., the cube has no time dimensions). - if frt_dims is None and fp_dims is None: - time_datetimes.append( - frt_coord.cell(0).point._to_real_datetime() - + pd.to_timedelta(fp_coord.points, unit=str(fp_coord.units)) - ) - elif frt_dims == fp_dims: - # Handle the case where forecast reference time and forecast period share - # the same dimension. This means both are mixed together along one - # dimension. - for frt, fp in zip( - list(frt_coord.cells()), - fp_coord.points, - ): - time_datetimes.append( - frt.point._to_real_datetime() - + pd.to_timedelta( - fp, unit=str(forecast_cube.coord("forecast_period").units) - ) - ) - else: - # Handle the case where forecast reference time and forecast period are on - # different dimensions. This requires iterating over all combinations of - # forecast period and forecast reference time. - enforce_coordinate_ordering( - forecast_cube, - self.expected_coordinate_order[::-1], - anchor_start=True, - ) - for fp in fp_coord.points: - for frt in list(frt_coord.cells()): - time_datetimes.append( - frt.point._to_real_datetime() - + pd.to_timedelta( - fp, unit=str(forecast_cube.coord("forecast_period").units) - ) - ) - - truth_data_list = [] - for time_datetime in time_datetimes: - constr = iris.Constraint(time=lambda cell: cell.point == time_datetime) - truth_data_list.append(truth_cube.extract(constr).data) - return truth_data_list - def process( self, - forecast_cube: Cube, - truth_cube: Cube, - feature_cubes: Optional[CubeList] = None, + forecast_df: pd.DataFrame, + truth_df: pd.DataFrame, ) -> None: """Train a quantile regression random forests model. Args: - forecast_cube: - Cube containing the realization forecasts. If the cube provided - contains multiple forecast periods, then the cube is expected to have - forecast period, forecast reference time, realization - and spot_index as the dimensions. This cube is only used as a template. - If the forecast cube is a feature cube, then it should also be provided - within the feature_cubes list. - truth_cube: + forecast_df: + DataFrame containing the forecast information and features. + truth_df: Cube containing the truths. The truths should have the same validity - times as the forecast. If the same validity time occurs multiple times - within the forecast cube (i.e. due to e.g. a forecast reference time - of 20170102T0000Z and a lead time of T+36 having the same validity time - as a forecast reference time of 20170103T0000Z and a lead time of T+6), - then the truth data will be repeated. - feature_cubes: - List of additional feature cubes. The name of the cube should match a - key in the feature_config dictionary. + times as the forecasts. References: Johnson. (2024). quantile-forest: A Python Package for Quantile @@ -489,52 +252,31 @@ def process( """ if self.transformation: - forecast_cube.data = getattr(np, self.transformation)( - forecast_cube.data + self.pre_transform_addition + forecast_df[self.target_name] = getattr(np, self.transformation)( + forecast_df[self.target_name] + self.pre_transform_addition ) - truth_cube.data = getattr(np, self.transformation)( - truth_cube.data + self.pre_transform_addition - ) - - # Ensure the forecast cube has the correct dimension ordering for prep_feature. - - coord_dims = [ - forecast_cube.coord_dims(c) for c in self.expected_coordinate_order - ] - coord_names = [] - for dim in coord_dims: - if len(dim) == 0: - continue - if len(forecast_cube.coords(dimensions=dim, dim_coords=True)) == 0: - continue - coord_names.append( - forecast_cube.coord(dimensions=dim, dim_coords=True).name() + truth_df["ob_value"] = getattr(np, self.transformation)( + truth_df["ob_value"] + self.pre_transform_addition ) - enforce_coordinate_ordering( - forecast_cube, - [coord_names], - anchor_start=True, - ) - - feature_values = [] for feature_name in self.feature_config.keys(): - feature_cube = feature_cubes.extract(iris.Constraint(feature_name)) - if not feature_cube: - msg = ( - f"Feature cube for {feature_name} not found in the provided " - "feature cubes." - ) + if feature_name not in forecast_df.columns: + msg = f"Feature '{feature_name}' is not present in the forecast DataFrame." raise ValueError(msg) for feature in self.feature_config[feature_name]: - print("feature = ", feature) - feature_values.append( - prep_feature(forecast_cube, feature_cube[0], feature) - ) - feature_values = np.array(feature_values).T - truth_data_list = self._organise_truth_data(forecast_cube, truth_cube) - target_values = np.array(truth_data_list).flatten() + forecast_df = prep_feature(forecast_df, feature_name, feature) + forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config) + feature_column_names = select_features(forecast_df, self.feature_config) + + merge_columns = ["wmo_id", "time"] + combined_df = forecast_df.merge( + truth_df[merge_columns + ["ob_value"]], on=merge_columns, how="inner" + ) + feature_values = np.array(combined_df[feature_column_names]) + target_values = combined_df["ob_value"].values + + print(feature_values, target_values) # Fit the quantile regression model qrf_model = self.fit_qrf(feature_values, target_values) @@ -546,6 +288,7 @@ class ApplyQuantileRegressionRandomForests(PostProcessingPlugin): def __init__( self, + target_name: str, feature_config: dict[str, list[str]], quantiles: list[np.float32], transformation: str = None, @@ -553,22 +296,24 @@ def __init__( ) -> None: """Initialise the plugin. Args: + target_name (str): + Name of the target variable to be calibrated. feature_config (dict): - Feature configuration defining the features to be used for quantile regression. - The configuration is a dictionary of strings, where the keys are the names of - the input cube(s) supplied, and the values are a list. This list can contain both - computed features, such as the mean or standard deviation (std), or static - features, such as the altitude. The computed features will be computed using - the cube defined in the dictionary key. If the key is the feature itself e.g. - a distance to water cube, then the value should state "static". This will ensure - the cube's data is used as the feature. - The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: - { - "air_temperature": ["mean", "std", "altitude"], - "visibility_at_screen_level": ["mean", "std"] - "distance_to_water": ["static"], - } + Feature configuration defining the features to be used for quantile + regression. The configuration is a dictionary of strings, where the + keys are the names of the columns within the dataframe. Some + features may be used as initially provided within the dataframe, + whilst others may be computed from the data e.g. mean, std. + If the key is the feature itself e.g. distance to water, then the value + should state "static". In this case, the name of feature e.g. + 'distance_to_water' is expected to be a column name in the input + dataframe. The config will have the structure: + "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + { + "air_temperature": ["mean", "std", "altitude"], + "visibility_at_screen_level": ["mean", "std"] + "distance_to_water": ["static"], + } quantiles (float): Quantiles used for prediction (values ranging from 0 to 1) transformation (str): @@ -579,14 +324,14 @@ def __init__( Raises: ValueError: If the transformation is not one of the supported types. """ - + self.target_name = target_name self.feature_config = feature_config self.quantiles = quantiles self.transformation = transformation _check_valid_transformation(self.transformation) self.pre_transform_addition = pre_transform_addition - def _reverse_transformation(self, forecast_cube: Cube): + def _reverse_transformation(self, forecast: np.ndarray) -> np.ndarray: """Reverse the transformation applied to the data prior to fitting the QRF. The forecast cube provided is modified in place. @@ -595,23 +340,19 @@ def _reverse_transformation(self, forecast_cube: Cube): """ if self.transformation: if self.transformation == "log": - forecast_cube.data = ( - np.exp(forecast_cube.data) - self.pre_transform_addition - ) + forecast = np.exp(forecast) - self.pre_transform_addition elif self.transformation == "log10": - forecast_cube.data = ( - 10 ** (forecast_cube.data) - self.pre_transform_addition - ) + forecast = 10 ** (forecast) - self.pre_transform_addition elif self.transformation == "sqrt": - forecast_cube.data = forecast_cube.data**2 - self.pre_transform_addition + forecast = forecast**2 - self.pre_transform_addition elif self.transformation == "cbrt": - forecast_cube.data = forecast_cube.data**3 - self.pre_transform_addition + forecast = forecast**3 - self.pre_transform_addition + return forecast def process( self, qrf_model: RandomForestQuantileRegressor, - feature_cubes: CubeList, - template_forecast_cube: Cube, + forecast_df: pd.DataFrame, ) -> Cube: """Apply a quantile regression random forests model. @@ -632,32 +373,29 @@ def process( feature_values = [] for feature_name in self.feature_config.keys(): - feature_cube = feature_cubes.extract_cube(iris.Constraint(feature_name)) - # Transform the feature cube data if a transformation is specified. if ( self.transformation and set(["mean", "std"]).intersection(self.feature_config[feature_name]) - and feature_cube.name() == template_forecast_cube.name() + and self.target_name in forecast_df.columns ): - feature_cube.data = getattr(np, self.transformation)( - feature_cube.data + self.pre_transform_addition + forecast_df[self.target_name] = getattr(np, self.transformation)( + forecast_df[self.target_name] + self.pre_transform_addition ) for feature in self.feature_config[feature_name]: - feature_values.append( - prep_feature(template_forecast_cube, feature_cube, feature) - ) + forecast_df = prep_feature(forecast_df, feature_name, feature) + + forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config) + feature_column_names = select_features(forecast_df, self.feature_config) + + feature_values = np.array(forecast_df[feature_column_names]) - feature_values = np.array(feature_values).T calibrated_forecast = qrf_model.predict( feature_values, quantiles=self.quantiles ) calibrated_forecast = np.float32(calibrated_forecast) - calibrated_forecast_cube = template_forecast_cube.copy( - data=np.broadcast_to(calibrated_forecast.T, template_forecast_cube.shape) - ) - self._reverse_transformation(calibrated_forecast_cube) + calibrated_forecast = self._reverse_transformation(calibrated_forecast) - return calibrated_forecast_cube + return calibrated_forecast diff --git a/improver/utilities/compare.py b/improver/utilities/compare.py index 75b943ae0c..99c8a1d9f1 100644 --- a/improver/utilities/compare.py +++ b/improver/utilities/compare.py @@ -406,10 +406,18 @@ def raise_reporter(message): reporter(str(exc)) return - assert output.n_features_in_ == kgo.n_features_in_ - assert output.n_outputs_ == kgo.n_outputs_ - assert output.max_depth == kgo.max_depth - assert output.n_estimators == kgo.n_estimators - assert output.random_state == kgo.random_state - for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_): - assert (output_estimator.tree_.value == kgo_estimator.tree_.value).all() + difference_found = False + try: + assert output.n_features_in_ == kgo.n_features_in_ + assert output.n_outputs_ == kgo.n_outputs_ + assert output.max_depth == kgo.max_depth + assert output.n_estimators == kgo.n_estimators + assert output.random_state == kgo.random_state + for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_): + assert (output_estimator.tree_.value == kgo_estimator.tree_.value).all() + except AssertionError: + difference_found = True + # call the reporter function outside the except block to avoid nested + # exceptions if the reporter function is raising an exception + if difference_found: + reporter("different pickled forest") diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index ed087c4f5d..b9250ace64 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -87,9 +87,9 @@ ed8ab78a9f55b54bf0a49f191d3eb33daae30e0d05b9051275a70c0e697aac71 ./apply-night- d2f7d8389b33cde359dd253def2aaf31afbc27557f389bf0f744a28b24e145dd ./apply-quantile-regression-random-forest/config.json 8db1b35bde734c16340a6e42454b9ac146ea32f7c012c17c5e64815069af41bf ./apply-quantile-regression-random-forest/input_forecast.nc 793f17191d2421ead8d546319c7da34143ef01e2db88ca6cc5faf46e92b4b1de ./apply-quantile-regression-random-forest/with_transformation_input.pickle -e5221acb444360ee802dc48e4da1ab3a874a4179ceff1a611138eb2af3f6591c ./apply-quantile-regression-random-forest/with_transformation_kgo.nc +d8acadfe782e2c2562dfcc31ada8ad381b320fbdf3854d6b361c38308e6b8226 ./apply-quantile-regression-random-forest/with_transformation_kgo.nc 31eab8a83a75eb8cee4c0140749a4d4e9d9b98292c28b3c53ee27756c1f3a265 ./apply-quantile-regression-random-forest/without_transformation_input.pickle -12ae824d2c23362eae37f41a3951c6d3fc48150e45a5998a9eafe467e34f0c1d ./apply-quantile-regression-random-forest/without_transformation_kgo.nc +65397673caf0f1205d18a85c08ce81bd3543a71c4101f40f4869a93ddd0cf0f8 ./apply-quantile-regression-random-forest/without_transformation_kgo.nc 74a40c508b68294652f5dde52d8fe4b93ea002bd4b72ca8d6cc11a83309f28f1 ./apply-rainforests-calibration/basic/kgo.nc b4678f178e0e1df0c242c61b269b71dd9a3d7aca64663e2e81e5bfeb71fe97da ./apply-rainforests-calibration/features/20200802T0000Z-PT0000H00M-clearsky_solar_radiation-PT24H.nc dce2d4a28cec494c94b154dc3a9203b46b7cef2ad54d04049663166f4b42ecf5 ./apply-rainforests-calibration/features/20200802T0000Z-PT0024H00M-cape-PT24H.nc @@ -1003,8 +1003,8 @@ e3d2315b24bf769bcfb137033950d7100c8795878635e0ec9491413a1349904a ./train-quanti 1aec562fc0c39c7695521482b889d335a3becb10a76b7eae965447cbd9faf051 ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/0600Z_0.parquet e81c04f9f2d8dc14b8d8cd9bc9e8827ebd0da5ba016db48cf60ab738d327f968 ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1200Z_0.parquet 99e6656982672ea2d4cb9fa1dfb5731408dd41b2c9052b7e0df8d8004864241e ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1800Z_0.parquet -793f17191d2421ead8d546319c7da34143ef01e2db88ca6cc5faf46e92b4b1de ./train-quantile-regression-random-forest/with_transformation_kgo.pickle -31eab8a83a75eb8cee4c0140749a4d4e9d9b98292c28b3c53ee27756c1f3a265 ./train-quantile-regression-random-forest/without_transformation_kgo.pickle +9d486570ca562c63aa45ed5f7621ecb00395891af7e021e382219cce74d484e8 ./train-quantile-regression-random-forest/with_transformation_kgo.pickle +f3919673f39a7db108add9aa5e4edd7706b6fd03a06e45f9c1d89070740524d4 ./train-quantile-regression-random-forest/without_transformation_kgo.pickle e2cfb5e19ef5ebcfd73dc1c8504b66e9a55ad76b7b18cb79318659224079f234 ./uv-index/basic/20181210T0600Z-PT0000H00M-radiation_flux_in_uv_downward_at_surface.nc e48c3b07bd14214a1a10c0bbf1e0acdf54cfedf93b2c6d766f594ce83a45c2b8 ./uv-index/basic/kgo.nc 2226a5e95eb29664e21f8c0658101a53d65945a1450a60a19955563a2693d22d ./vertical-updraught/cape.nc diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py index c8fb8bff6f..6d237f560b 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py @@ -6,6 +6,7 @@ import numpy as np import pytest +from iris.coords import AuxCoord from iris.cube import Cube from improver.calibration.load_and_apply_quantile_regression_random_forest import ( @@ -16,13 +17,31 @@ ) from improver.utilities.save import save_netcdf from improver_tests.calibration.quantile_regression_random_forests_calibration.test_quantile_regression_random_forest import ( - _add_day_of_training_period, _create_ancil_file, _create_forecasts, _run_train_qrf, ) +def _add_day_of_training_period_to_cube(cube, day_of_training_period, secondary_coord): + """Add day of training period coordinate to the cube. + Args: + cube: Cube to which the day of training period coordinate will be added. + day_of_training_period: Day of training period to be added. + secondary_coord: Coordinate to associate the day of training period with. + Returns: + Cube with the day of training period coordinate added. + """ + dims = cube.coord_dims(secondary_coord) + day_of_training_period_coord = AuxCoord( + np.array(day_of_training_period, dtype=np.int32), + long_name="day_of_training_period", + units="1", + ) + cube.add_aux_coord(day_of_training_period_coord, data_dims=dims) + return cube + + @pytest.mark.parametrize("percentile_input", [True, False]) @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,quantiles,expected", @@ -31,10 +50,10 @@ (100, 2, 55, 5, None, 0, {}, False, [1 / 3, 2 / 3], [[4.1, 5.1], [5.1, 5.1]]), # noqa Multiple quantiles (1, 1, 55, 5, None, 0, {}, False, [0.5], [4.1, 6.2]), # noqa Fewer estimators and reduced depth (1, 1, 73, 5, None, 0, {}, False, [0.5], [4.2, 6.2]), # Different random state - (2, 2, 55, 5, "log", 10, {}, False, [0.5], [4.1, 5.1]), # Log transformation - (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Log10 transformation - (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Square root transformation - (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [4.1, 5.1]), # noqa Cube root transformation + (2, 2, 55, 5, "log", 10, {}, False, [0.5], [4.1, 5.64]), # Log transformation + (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [4.1, 5.64]), # noqa Log10 transformation + (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [4.1, 5.64]), # noqa Square root transformation + (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [4.1, 5.64]), # noqa Cube root transformation (2, 2, 55, 5, None, 0, {"max_samples_leaf": 0.5}, False, [0.5], [4.1, 6.2]), # noqa # Different criterion (2, 5, 55, 5, None, 0, {}, True, [0.5], [4.1, 5.65]), # noqa Include an additional static feature ], @@ -75,20 +94,20 @@ def test_load_and_apply_qrf( "20170101T1200Z", "20170102T1200Z", ], - day_of_training_period=[0, 1], realization_data=[2, 6, 10], - truth_data=[[4.2, 6.2], [4.1, 5.1]], + truth_data=[4.2, 6.2, 4.1, 5.1], ) frt = "20170103T0000Z" vt = "20170103T1200Z" - day_of_training_period = 2 data = np.arange(6, (len(quantiles) * 6) + 1, 6) + day_of_training_period = 2 + forecast_cube = _create_forecasts(frt, vt, data, return_cube=True) - forecast_cube = _create_forecasts(frt, vt, data) - forecast_cube = _add_day_of_training_period( + forecast_cube = _add_day_of_training_period_to_cube( forecast_cube, day_of_training_period, "forecast_reference_time" ) + if percentile_input: forecast_cube = RebadgeRealizationsAsPercentiles()(forecast_cube) @@ -100,7 +119,7 @@ def test_load_and_apply_qrf( file_paths = [model_output, forecast_filepath] if include_static: - ancil_cube = _create_ancil_file() + ancil_cube = _create_ancil_file(return_cube=True) ancil_filepath = features_dir / "ancil.nc" save_netcdf(ancil_cube, ancil_filepath) file_paths.append(str(ancil_filepath)) @@ -114,7 +133,7 @@ def test_load_and_apply_qrf( result = plugin.process(file_paths=file_paths) assert isinstance(result, Cube) - assert result.data.shape == (1, 1, len(quantiles), 2) + assert result.data.shape == (len(quantiles), 2) assert np.allclose(result.data, expected, rtol=1e-2) # Check that the metadata is as expected @@ -178,22 +197,20 @@ def test_exceptions( "20170101T1200Z", "20170102T1200Z", ], - day_of_training_period=[0, 1], realization_data=[2, 6, 10], - truth_data=[[4.2, 6.2], [4.1, 5.1]], + truth_data=[4.2, 6.2, 4.1, 5.1], ) frt = "20170103T0000Z" vt = "20170103T1200Z" - day_of_training_period = 2 data = np.arange(6, (len(quantiles) * 6) + 1, 6) - - forecast_cube = _create_forecasts(frt, vt, data) - forecast_cube = _add_day_of_training_period( + day_of_training_period = 2 + forecast_cube = _create_forecasts(frt, vt, data, return_cube=True) + forecast_cube = _add_day_of_training_period_to_cube( forecast_cube, day_of_training_period, "forecast_reference_time" ) - ancil_cube = _create_ancil_file() + ancil_cube = _create_ancil_file(return_cube=True) ancil_filepath = tmp_path / "ancil.nc" save_netcdf(ancil_cube, ancil_filepath) diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py index 51925d984a..f4642019a2 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py @@ -362,7 +362,7 @@ def _create_ancil_file(tmp_path, wmo_ids): False, False, "percentile", - 5.64, + 5.61, ), ], ) @@ -391,7 +391,7 @@ def test_load_and_train_qrf( model_output = str(model_output_dir / "qrf_model.pkl") if include_static: - ancil_path = _create_ancil_file(tmp_path, list(set(wmo_ids))) + ancil_path = _create_ancil_file(tmp_path, sorted(list(set(wmo_ids)))) file_paths.append(ancil_path) feature_config["distance_to_water"] = ["static"] @@ -423,7 +423,7 @@ def test_load_and_train_qrf( if remove_target: current_forecast = [] else: - current_forecast = [279, 3, 55] + current_forecast = [5.64, 3, 55] if include_static: current_forecast.append(2.5) diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py index 35591920df..15eef7616b 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py @@ -11,29 +11,32 @@ import numpy as np import pandas as pd import pytest -from iris.coords import AuxCoord -from iris.cube import Cube, CubeList +from iris.cube import Cube +from iris.pandas import as_data_frame from improver.calibration.quantile_regression_random_forest import ( ApplyQuantileRegressionRandomForests, TrainQuantileRegressionRandomForests, prep_feature, + select_features, ) from improver.metadata.constants.time_types import DT_FORMAT from improver.synthetic_data.set_up_test_cubes import set_up_spot_variable_cube -from improver.utilities.temporal import datetime_to_iris_time ALTITUDE = [10, 20] LATITUDE = [50, 60] LONGITUDE = [0, 10] WMO_ID = ["00001", "00002"] +iris.FUTURE.pandas_ndim = True + def _create_forecasts( forecast_reference_time: str, validity_time: str, data: list[int], -) -> Cube: + return_cube: bool = False, +) -> Cube | pd.DataFrame: """Create site forecast cube with realizations. Args: @@ -58,36 +61,48 @@ def _create_forecasts( time=dt.strptime(validity_time, DT_FORMAT), frt=dt.strptime(forecast_reference_time, DT_FORMAT), ) - cube.remove_coord("time") - return cube + if return_cube: + return cube + df = as_data_frame( + cube, + add_aux_coords=[ + "altitude", + "latitude", + "longitude", + "wmo_id", + "forecast_period", + "forecast_reference_time", + "time", + ], + ).reset_index() + df["time"] = df["time"].apply(lambda x: x._to_real_datetime()) + df["forecast_reference_time"] = df["forecast_reference_time"].apply( + lambda x: x._to_real_datetime() + ) + return df -def _add_day_of_training_period(cube, day_of_training_period, secondary_coord): - """Add day of training period coordinate to the cube. +def _add_day_of_training_period(df): + """Add day of training period coordinate to the dataframe. Args: - cube: Cube to which the day of training period coordinate will be added. + df: DataFrame to which the day of training period coordinate will be added. day_of_training_period: Day of training period to be added. - secondary_coord: Coordinate to associate the day of training period with. Returns: Cube with the day of training period coordinate added. """ - dims = cube.coord_dims(secondary_coord) - day_of_training_period_coord = AuxCoord( - np.array(day_of_training_period, dtype=np.int32), - long_name="day_of_training_period", - units="1", + df["day_of_training_period"] = df["time"].dt.dayofyear - np.min( + df["time"].dt.dayofyear ) - cube.add_aux_coord(day_of_training_period_coord, data_dims=dims) - return cube + return df -def _create_ancil_file(): +def _create_ancil_file(return_cube=False): """Create an ancillary file for testing. Returns: - An ancillary cube with a single value. + An ancillary DataFrame without temporal columns. """ data = np.array([2, 3], dtype=np.float32) template_cube = set_up_spot_variable_cube( @@ -102,7 +117,17 @@ def _create_ancil_file(): cube = template_cube.copy() for coord in ["time", "forecast_reference_time", "forecast_period"]: cube.remove_coord(coord) - return cube + if return_cube: + return cube + return as_data_frame( + cube, + add_aux_coords=[ + "altitude", + "latitude", + "longitude", + "wmo_id", + ], + ).reset_index() def _run_train_qrf( @@ -128,65 +153,49 @@ def _run_train_qrf( "20170102T1200Z", "20170102T1800Z", ], - day_of_training_period=[0, 1], realization_data=[2, 6, 10], - truth_data=[[4.2, 3.8], [5.8, 6], [7, 7.3], [9.1, 9.5]], + truth_data=[4.2, 3.8, 5.8, 6, 7, 7.3, 9.1, 9.5], ): realization_data = np.array(realization_data, dtype=np.float32) - forecast_cubes = CubeList() + forecast_dfs = [] for index, (frt, vt) in enumerate(zip(forecast_reference_times, validity_times)): - forecast_cube = _create_forecasts(frt, vt, realization_data + index) - forecast_cubes.append(forecast_cube) - forecast_cube = forecast_cubes.merge_cube() - - forecast_cube = _add_day_of_training_period( - forecast_cube, day_of_training_period, "forecast_reference_time" + forecast_df = _create_forecasts(frt, vt, realization_data + index) + forecast_dfs.append(forecast_df) + forecast_df = pd.concat(forecast_dfs) + forecast_df = _add_day_of_training_period(forecast_df) + + truth_df = forecast_df.copy() + truth_df.drop( + columns=[ + "forecast_period", + "forecast_reference_time", + "day_of_training_period", + "realization", + "wind_speed_at_10m", + ], + inplace=True, ) + truth_df.drop_duplicates(inplace=True) + truth_df["ob_value"] = np.array(truth_data, dtype=np.float32) - truth_cube = next(forecast_cube.slices_over("realization")).copy() - for coord in [ - "day_of_training_period", - "realization", - ]: - truth_cube.remove_coord(coord) - truth_slices = CubeList() - for truth_slice in truth_cube.slices_over( - ["forecast_period", "forecast_reference_time"] - ): - frt_coord = truth_slice.coord("forecast_reference_time") - fp_coord = truth_slice.coord("forecast_period") - time_point = datetime_to_iris_time( - frt_coord.cell(0).point._to_real_datetime() - + pd.to_timedelta(fp_coord.points[0], unit=str(fp_coord.units)) - ) - time_coord = iris.coords.AuxCoord( - time_point, standard_name="time", units=frt_coord.units - ) - truth_slice.add_aux_coord(time_coord) - truth_slice.remove_coord("forecast_period") - truth_slice.remove_coord("forecast_reference_time") - truth_slices.append(truth_slice) - - truth_cube = truth_slices.merge_cube() - truth_cube.data = np.array(truth_data).astype(np.float32) - - feature_cubes = CubeList([forecast_cube]) if include_static: - ancil_cube = _create_ancil_file() - feature_cubes.append(ancil_cube) - feature_config[ancil_cube.name()] = ["static"] + ancil_df = _create_ancil_file() + forecast_df = forecast_df.merge( + ancil_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left" + ) + feature_config["distance_to_water"] = ["static"] if "air_temperature" in feature_config.keys(): - dynamic_cube = forecast_cube.copy(data=forecast_cube.data + 10) - dynamic_cube.rename("air_temperature") - dynamic_cube.units = "Celsius" - feature_cubes.append(dynamic_cube) + forecast_df["air_temperature"] = np.array( + forecast_df["wind_speed_at_10m"] + 10, dtype=np.float32 + ) model_output_dir = tmp_path / "train_qrf" model_output_dir.mkdir(parents=True) model_output = str(model_output_dir / "qrf_model.pkl") plugin = TrainQuantileRegressionRandomForests( - feature_config, + target_name="wind_speed_at_10m", + feature_config=feature_config, n_estimators=n_estimators, max_depth=max_depth, random_state=random_state, @@ -196,172 +205,129 @@ def _run_train_qrf( model_output=model_output, **extra_kwargs, ) - plugin.process(forecast_cube, truth_cube, feature_cubes=feature_cubes) + plugin.process(forecast_df, truth_df) return model_output @pytest.mark.parametrize( "feature,expected,expected_dtype", [ - ("mean", np.array([6, 8], dtype=np.float32), np.float32), - ("std", np.array([4, 4], dtype=np.float32), np.float32), - ("latitude", np.array([50, 60], dtype=np.float32), np.float32), - ("longitude", np.array([0, 10], dtype=np.float32), np.float32), - ("altitude", np.array([10, 20], dtype=np.float32), np.float32), + ("mean", np.tile(np.array([6, 8], dtype=np.float32), 3), np.float32), + ("std", np.repeat(4, 6).astype(np.float32), np.float32), + ("latitude", np.tile(np.array([50, 60], dtype=np.float32), 3), np.float32), + ("longitude", np.tile(np.array([0, 10], dtype=np.float32), 3), np.float32), + ("altitude", np.tile(np.array([10, 20], dtype=np.float32), 3), np.float32), ( "day_of_year", - np.array([1, 1], dtype=np.float32), - np.float32, + np.tile(np.array([1, 1], dtype=np.int32), 3), + np.int32, ), ( "day_of_year_sin", - np.array([0.01716633, 0.01716633], dtype=np.float32), + np.tile(np.array([0.01716633, 0.01716633], dtype=np.float32), 3), np.float32, ), ( "day_of_year_cos", - np.array([0.99985266, 0.99985266], dtype=np.float32), + np.tile(np.array([0.99985266, 0.99985266], dtype=np.float32), 3), np.float32, ), - ("hour_of_day", np.array([12, 12], dtype=np.float32), np.float32), - ("hour_of_day_sin", np.array([0, 0], dtype=np.float32), np.float32), - ("hour_of_day_cos", np.array([-1, -1], dtype=np.float32), np.float32), - ("forecast_period", np.array([43200, 43200], dtype=np.int32), np.int32), - ("day_of_training_period", np.array([0, 0], dtype=np.int32), np.int32), - ("static", np.array([2, 3], dtype=np.float32), np.float32), + ("hour_of_day", np.repeat(12, 6).astype(np.int32), np.int32), + ("hour_of_day_sin", np.repeat(0, 6).astype(np.float32), np.float32), + ("hour_of_day_cos", np.repeat(-1, 6).astype(np.float32), np.float32), + ("forecast_period", np.repeat(43200, 6).astype(np.int32), np.int32), + ("day_of_training_period", np.repeat(0, 6).astype(np.int32), np.int32), + ("static", np.tile(np.array([2, 3], dtype=np.float32), 3), np.float32), ], ) def test_prep_feature_single_time(feature, expected, expected_dtype): """Test the prep_feature function for a single time.""" + feature_name = "wind_speed_at_10m" + forecast_reference_time = "20170101T0000Z" validity_time = "20170101T1200Z" data = np.array([2, 6, 10]) - day_of_training_period = [0] - forecast_cube = _create_forecasts(forecast_reference_time, validity_time, data) - forecast_cube = _add_day_of_training_period( - forecast_cube, day_of_training_period, "forecast_reference_time" - ) + forecast_df = _create_forecasts(forecast_reference_time, validity_time, data) + forecast_df = _add_day_of_training_period(forecast_df) if feature == "static": - feature_cube = _create_ancil_file() + feature_df = _create_ancil_file() + forecast_df = forecast_df.merge( + feature_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left" + ) else: - feature_cube = forecast_cube.copy() + forecast_df = forecast_df.copy() + + result = prep_feature(forecast_df, feature_name, feature) + + if feature in ["mean", "std"]: + assert result.shape == (6, 12) + feature_name_modified = f"{feature_name}_{feature}" + elif feature in [ + "day_of_year", + "day_of_year_sin", + "day_of_year_cos", + "hour_of_day", + "hour_of_day_sin", + "hour_of_day_cos", + ]: + assert result.shape == (6, 12) + feature_name_modified = feature + elif feature in ["static"]: + assert result.shape == (6, 12) + feature_name_modified = "distance_to_water" + else: + assert result.shape == (6, 11) + feature_name_modified = feature - result = prep_feature(forecast_cube, feature_cube, feature) - assert result.shape == (2,) - assert result.dtype == expected_dtype - np.testing.assert_allclose(result, expected, atol=1e-6) + assert result[feature_name_modified].dtype == expected_dtype + np.testing.assert_allclose(result[feature_name_modified], expected, atol=1e-6) @pytest.mark.parametrize( "feature,expected,expected_dtype", [ - ("mean", np.tile([6, 8], 4).astype(np.float32), np.float32), - ("std", np.repeat(4, 8).astype(np.float32), np.float32), - ("latitude", np.tile([50, 60], 4).astype(np.float32), np.float32), - ("longitude", np.tile([0, 10], 4).astype(np.float32), np.float32), - ("altitude", np.tile([10, 20], 4).astype(np.float32), np.float32), + ("mean", np.tile([6, 8], 18).astype(np.float32), np.float32), + ("std", np.repeat(4, 36).astype(np.float32), np.float32), + ("latitude", np.tile([50, 60], 18).astype(np.float32), np.float32), + ("longitude", np.tile([0, 10], 18).astype(np.float32), np.float32), + ("altitude", np.tile([10, 20], 18).astype(np.float32), np.float32), ( "day_of_year", - np.repeat([1, 2], 4).astype(np.float32), - np.float32, + np.repeat([1, 2, 3], 12).astype(np.int32), + np.int32, ), ( "day_of_year_sin", - np.repeat([0.017166, 0.034328], 4).astype(np.float32), + np.repeat([0.017166, 0.034328, 0.051479], 12).astype(np.float32), np.float32, ), ( "day_of_year_cos", - np.repeat([0.999853, 0.999411], 4).astype(np.float32), + np.repeat([0.999853, 0.999411, 0.998674], 12).astype(np.float32), np.float32, ), - ("hour_of_day", np.repeat([12, 12], 4).astype(np.float32), np.float32), - ("hour_of_day_sin", np.repeat([0, 0], 4).astype(np.float32), np.float32), - ("hour_of_day_cos", np.repeat([-1, -1], 4).astype(np.float32), np.float32), + ("hour_of_day", np.tile(np.repeat([6, 12], 6), 3).astype(np.int32), np.int32), ( - "forecast_period", - np.tile([43200, 43200, 21600, 21600], 2).astype(np.int32), - np.int32, - ), - ("day_of_training_period", np.repeat([0, 1], 4).astype(np.int32), np.int32), - ("static", np.tile([2, 3], 4).astype(np.float32), np.float32), - ], -) -def test_prep_feature_1d_time_dimension(feature, expected, expected_dtype): - """Test the prep_feature function for multiple times.""" - forecast_reference_times = [ - "20170101T0000Z", - "20170101T0600Z", - "20170102T0000Z", - "20170102T0600Z", - ] - validity_times = [ - "20170101T1200Z", - "20170101T1200Z", - "20170102T1200Z", - "20170102T1200Z", - ] - - data = np.array([2, 6, 10]) - - forecast_cubes = CubeList() - for frt, vt in zip(forecast_reference_times, validity_times): - forecast_cubes.append(_create_forecasts(frt, vt, data)) - forecast_cube = forecast_cubes.merge_cube() - day_of_training_period = [0, 0, 1, 1] - forecast_cube = _add_day_of_training_period( - forecast_cube, day_of_training_period, "forecast_reference_time" - ) - - if feature == "static": - feature_cube = _create_ancil_file() - else: - feature_cube = forecast_cube.copy() - - result = prep_feature(forecast_cube, feature_cube, feature) - - assert result.shape == (8,) - assert result.dtype == expected_dtype - np.testing.assert_allclose(result, expected, atol=1e-6) - - -@pytest.mark.parametrize( - "feature,expected,expected_dtype", - [ - ("mean", np.tile([6, 8], 6).astype(np.float32), np.float32), - ("std", np.repeat(4, 12).astype(np.float32), np.float32), - ("latitude", np.tile([50, 60], 6).astype(np.float32), np.float32), - ("longitude", np.tile([0, 10], 6).astype(np.float32), np.float32), - ("altitude", np.tile([10, 20], 6).astype(np.float32), np.float32), - ( - "day_of_year", - np.repeat([1, 2, 3], 4).astype(np.float32), + "hour_of_day_sin", + np.tile(np.repeat([1, 0], 6), 3).astype(np.float32), np.float32, ), ( - "day_of_year_sin", - np.repeat([0.017166, 0.034328, 0.051479], 4).astype(np.float32), - np.float32, - ), - ( - "day_of_year_cos", - np.repeat([0.999853, 0.999411, 0.998674], 4).astype(np.float32), + "hour_of_day_cos", + np.tile(np.repeat([0, -1], 6), 3).astype(np.float32), np.float32, ), - ("hour_of_day", np.tile([6, 6, 12, 12], 3).astype(np.float32), np.float32), - ("hour_of_day_sin", np.tile([1, 1, 0, 0], 3).astype(np.float32), np.float32), - ("hour_of_day_cos", np.tile([0, 0, -1, -1], 3).astype(np.float32), np.float32), ( "forecast_period", - np.tile([21600, 21600, 43200, 43200], 3).astype(np.int32), + np.tile(np.repeat([21600, 43200], 6), 3).astype(np.int32), np.int32, ), - ("day_of_training_period", np.repeat([0, 1, 2], 4).astype(np.int32), np.int32), - ("static", np.tile([2, 3], 6).astype(np.float32), np.float32), + ("day_of_training_period", np.repeat([0, 1, 2], 12).astype(np.int32), np.int32), + ("static", np.tile([2, 3], 18).astype(np.float32), np.float32), ], ) -def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype): +def test_prep_feature_more_times(feature, expected, expected_dtype): """Test the prep_feature function for multiple times.""" forecast_reference_times = [ "20170101T0000Z", @@ -381,42 +347,109 @@ def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype): ] data = np.array([2, 6, 10]) - forecast_cubes = CubeList() + + feature_name = "wind_speed_at_10m" + + data = np.array([2, 6, 10]) + + forecast_dfs = [] for frt, vt in zip(forecast_reference_times, validity_times): - forecast_cubes.append(_create_forecasts(frt, vt, data)) - forecast_cube = forecast_cubes.merge_cube() - # Ensure coordinate order is forecast_reference_time, forecast_period, - # realization, spot_index - forecast_cube.transpose([1, 0, 2, 3]) - day_of_training_period = [0, 1, 2] - forecast_cube = _add_day_of_training_period( - forecast_cube, day_of_training_period, "forecast_reference_time" - ) + forecast_dfs.append(_create_forecasts(frt, vt, data)) + forecast_df = pd.concat(forecast_dfs) + forecast_df = _add_day_of_training_period(forecast_df) if feature == "static": - feature_cube = _create_ancil_file() + feature_df = _create_ancil_file() + forecast_df = forecast_df.merge( + feature_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left" + ) + else: - feature_cube = forecast_cube.copy() + forecast_df = forecast_df.copy() + + result = prep_feature(forecast_df, feature_name, feature) + + if feature in ["mean", "std"]: + assert result.shape == (36, 12) + feature_name_modified = f"{feature_name}_{feature}" + elif feature in [ + "day_of_year", + "day_of_year_sin", + "day_of_year_cos", + "hour_of_day", + "hour_of_day_sin", + "hour_of_day_cos", + ]: + assert result.shape == (36, 12) + feature_name_modified = feature + elif feature in ["static"]: + assert result.shape == (36, 12) + feature_name_modified = "distance_to_water" + else: + assert result.shape == (36, 11) + feature_name_modified = feature + + assert result[feature_name_modified].dtype == expected_dtype + np.testing.assert_allclose(result[feature_name_modified], expected, atol=1e-6) - result = prep_feature(forecast_cube, feature_cube, feature) - assert result.shape == (12,) - assert result.dtype == expected_dtype - np.testing.assert_allclose(result, expected, atol=1e-6) +@pytest.mark.parametrize( + "feature_config", + [ + {"wind_speed_at_10m": ["mean", "std", "latitude", "longitude"]}, + {"wind_speed_at_10m": ["latitude", "longitude"]}, + { + "wind_speed_at_10m": ["latitude", "longitude"], + "distance_to_water": ["static"], + }, + {"wind_speed_at_10m": ["latitude", "longitude", "height"]}, + ], +) +def test_select_features(feature_config): + data_dict = { + "wmo_id": np.tile(WMO_ID, 3), + "latitude": np.tile(LATITUDE, 3), + "longitude": np.tile(LONGITUDE, 3), + "altitude": np.tile(ALTITUDE, 3), + "wind_speed_at_10m_mean": np.repeat(5, 6), + "wind_speed_at_10m_std": np.repeat(1, 6), + "wind_speed_at_10m": np.tile([4, 6], 3), + "realization": np.tile([1, 2], 3), + "distance_to_water": np.tile([2.0, 3.0], 3), + } + df = pd.DataFrame(data_dict) + + expected = [] + for key, values in feature_config.items(): + for value in values: + if "mean" in value or "std" in value: + expected.append(f"wind_speed_at_10m_{value}") + elif value == "static": + expected.append(key) + else: + expected.append(value) + + if "height" in feature_config["wind_speed_at_10m"]: + with pytest.raises(ValueError, match="Feature 'height' is not supported."): + select_features(df, feature_config) + else: + result = select_features(df, feature_config) + assert len(result) == len(expected) + assert result == expected @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected", [ - (2, 2, 55, 5, None, 0, {}, False, 4.1), # Basic test case + (2, 2, 55, 5, None, 0, {}, False, 4.2), # Basic test case (2, 2, 54, 5, None, 0, {}, False, 4.15), # Different random state - (1, 1, 55, 5, None, 0, {}, False, 4.1), # Fewer estimators and reduced depth + (1, 1, 54, 5, None, 0, {}, False, 4.2), # Fewer estimators and reduced depth (2, 2, 55, 5, "log", 10, {}, False, 2.65), # Log transformation (2, 2, 55, 5, "log10", 10, {}, False, 1.15), # Log10 transformation (2, 2, 55, 5, "sqrt", 10, {}, False, 3.76), # Square root transformation (2, 2, 55, 5, "cbrt", 10, {}, False, 2.42), # Cube root transformation - (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.1), # noqa # Different criterion - (2, 5, 55, 5, None, 0, {}, True, 4.1), # Include static data + (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.2), # noqa # Different criterion + (2, 5, 55, 5, None, 0, {}, True, 4.2), # Include static data ], ) def test_train_qrf_single_lead_times( @@ -455,9 +488,8 @@ def test_train_qrf_single_lead_times( "20170101T1200Z", "20170102T1200Z", ], - day_of_training_period=[0, 1], realization_data=[2, 6, 10], - truth_data=[[4.2, 4.2], [4.1, 4.1]], + truth_data=[4.2, 4.1, 4.2, 4.1], ) qrf_model = joblib.load(model_output) @@ -465,7 +497,7 @@ def test_train_qrf_single_lead_times( assert qrf_model.max_depth == max_depth assert qrf_model.random_state == random_state - current_forecast = [5, 3, 55, 5] + current_forecast = [7.5, 3, 55, 5] if include_static: current_forecast.append(2.5) result = qrf_model.predict( @@ -477,15 +509,15 @@ def test_train_qrf_single_lead_times( @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected", [ - (2, 2, 55, 5, None, 0, {}, False, 5.6), # Basic test case + (2, 2, 55, 5, None, 0, {}, False, 5.8), # Basic test case (1, 1, 55, 5, None, 0, {}, False, 3.8), # Fewer estimators and reduced depth - (1, 1, 73, 5, None, 0, {}, False, 7.0), # Different random state - (2, 2, 55, 5, "log", 10, {}, False, 2.713), # Log transformation - (2, 2, 55, 5, "log10", 10, {}, False, 1.178), # Log10 transformation - (2, 2, 55, 5, "sqrt", 10, {}, False, 3.884), # Square root transformation - (2, 2, 55, 5, "cbrt", 10, {}, False, 2.471), # Cube root transformation - (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4), # noqa # Different criterion - (2, 5, 55, 5, None, 0, {}, True, 5), # Include static data + (1, 1, 73, 5, None, 0, {}, False, 5.8), # Different random state + (2, 2, 55, 5, "log", 10, {}, False, 2.752), # Log transformation + (2, 2, 55, 5, "log10", 10, {}, False, 1.195), # Log10 transformation + (2, 2, 55, 5, "sqrt", 10, {}, False, 3.964), # Square root transformation + (2, 2, 55, 5, "cbrt", 10, {}, False, 2.504), # Cube root transformation + (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 5.8), # noqa # Different criterion + (1, 1, 73, 5, None, 0, {}, True, 3.8), # Include static data ], ) def test_train_qrf_multiple_lead_times( @@ -523,9 +555,9 @@ def test_train_qrf_multiple_lead_times( assert qrf_model.max_depth == max_depth assert qrf_model.random_state == random_state - current_forecast = [5, 3, 55, 5] + current_forecast = [7.5, 3, 55, 5] if include_static: - current_forecast.append(2.5) + current_forecast.append(3.2) result = qrf_model.predict( np.expand_dims(np.array(current_forecast), 0), quantiles=[0.5] ) @@ -536,8 +568,8 @@ def test_train_qrf_multiple_lead_times( "feature_config,data,include_static,expected", [ ({"wind_speed_at_10m": ["mean"]}, [5], False, [5]), # One feature - ({"wind_speed_at_10m": ["latitude"]}, [55], False, [6.65]), # noqa Without the target - ({"wind_speed_at_10m": ["mean"]}, [5], True, [5]), # With static data + ({"wind_speed_at_10m": ["latitude"]}, [61], False, [4.9]), # noqa Without the target + ({"wind_speed_at_10m": ["mean"]}, [5], True, [4]), # With static data ( {"wind_speed_at_10m": ["mean"], "air_temperature": ["mean"]}, [5], @@ -548,7 +580,7 @@ def test_train_qrf_multiple_lead_times( {"wind_speed_at_10m": ["mean"], "pressure_at_mean_sea_level": ["mean"]}, [5], False, - "Feature cube for pressure_at_mean_sea_level", + "Feature 'pressure_at_mean_sea_level' is not present", ), # Multiple dynamic features ], ) @@ -605,7 +637,7 @@ def test_alternative_feature_configs( current_forecast = [*data] if include_static: - current_forecast.append(2.5) + current_forecast.append(3.2) if "air_temperature" in feature_config: current_forecast.append(15.0) @@ -618,31 +650,31 @@ def test_alternative_feature_configs( @pytest.mark.parametrize( "quantiles,transformation,pre_transform_addition,include_static,expected", [ - ([0.5], None, 0, False, [5.6, 5.4]), # One quantile - ([0.1, 0.5, 0.9], None, 0, False, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Multiple quantiles - ([0.1, 0.5, 0.9], "log", 10, False, [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]]), # noqa Log transformation - ( - [0.1, 0.5, 0.9], - "log10", - 10, - False, - [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]], - ), # Log10 transformation - ( - [0.1, 0.5, 0.9], - "sqrt", - 10, - False, - [[4.38, 4.38], [5.09, 5.09], [5.82, 5.82]], - ), # Square root transformation - ( - [0.1, 0.5, 0.9], - "cbrt", - 10, - False, - [[4.37, 4.37], [5.08, 5.08], [5.81, 5.81]], - ), # Cube root transformation - ([0.1, 0.5, 0.9], None, 0, True, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Include static data + # ([0.5], None, 0, False, [5, 4.9]), # One quantile + ([0.1, 0.5, 0.9], None, 0, False, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.15]]), # noqa Multiple quantiles + # ([0.1, 0.5, 0.9], "log", 10, False, [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]]), # noqa Log transformation + # ( + # [0.1, 0.5, 0.9], + # "log10", + # 10, + # False, + # [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]], + # ), # Log10 transformation + # ( + # [0.1, 0.5, 0.9], + # "sqrt", + # 10, + # False, + # [[4.38, 4.38], [5.09, 5.09], [5.82, 5.82]], + # ), # Square root transformation + # ( + # [0.1, 0.5, 0.9], + # "cbrt", + # 10, + # False, + # [[4.37, 4.37], [5.08, 5.08], [5.81, 5.81]], + # ), # Cube root transformation + # ([0.1, 0.5, 0.9], None, 0, True, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Include static data ], ) def test_apply_qrf( @@ -677,30 +709,28 @@ def test_apply_qrf( frt = "20170103T0000Z" vt = "20170103T1200Z" - day_of_training_period = 2 data = np.arange(6, (len(quantiles) * 6) + 1, 6) - forecast_cube = _create_forecasts(frt, vt, data) - forecast_cube = _add_day_of_training_period( - forecast_cube, day_of_training_period, "forecast_reference_time" - ) + forecast_df = _create_forecasts(frt, vt, data) + forecast_df = _add_day_of_training_period(forecast_df) - feature_cubes = CubeList([forecast_cube]) if include_static: - ancil_cube = _create_ancil_file() - feature_cubes.append(ancil_cube) + ancil_df = _create_ancil_file() + forecast_df = forecast_df.merge( + ancil_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left" + ) plugin = ApplyQuantileRegressionRandomForests( - feature_config, quantiles, transformation, pre_transform_addition + "air_temperature", + feature_config, + quantiles, + transformation, + pre_transform_addition, ) - result = plugin.process(qrf_model, feature_cubes, forecast_cube) + result = plugin.process(qrf_model, forecast_df) - assert isinstance(result, Cube) - assert result.shape == (len(quantiles), 2) + assert isinstance(result, np.ndarray) + assert result.shape == (2, 3) assert result.dtype == np.float32 - expected_cube = forecast_cube.copy( - data=np.full((len(quantiles), 2), expected, dtype=np.float32) - ) - assert result.name() == expected_cube.name() - assert result.units == expected_cube.units - np.testing.assert_almost_equal(result.data, expected_cube.data, decimal=2) + # expected_data = np.full((len(quantiles), 2), expected, dtype=np.float32) + np.testing.assert_almost_equal(result, expected, decimal=2) From 7165b67679788b686d5e0e4d3a0d5f70a886e4b3 Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Mon, 18 Aug 2025 15:55:59 +0100 Subject: [PATCH 4/8] Updates to tidy up re-write to use pandas DataFrames more widely. --- ...apply_quantile_regression_random_forest.py | 89 +++++----- ...train_quantile_regression_random_forest.py | 32 ++-- .../quantile_regression_random_forest.py | 168 +++++++++++------- ...apply_quantile_regression_random_forest.py | 28 +-- ...train_quantile_regression_random_forest.py | 6 +- improver_tests/acceptance/acceptance.py | 2 + ...apply_quantile_regression_random_forest.py | 2 +- ...train_quantile_regression_random_forest.py | 2 + ...train_quantile_regression_random_forest.py | 4 + .../test_quantile_regression_random_forest.py | 8 +- 10 files changed, 202 insertions(+), 139 deletions(-) diff --git a/improver/calibration/load_and_apply_quantile_regression_random_forest.py b/improver/calibration/load_and_apply_quantile_regression_random_forest.py index cdfaeae530..5647f64d3e 100644 --- a/improver/calibration/load_and_apply_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_apply_quantile_regression_random_forest.py @@ -12,6 +12,7 @@ import iris import joblib import numpy as np +import pandas as pd from iris.cube import Cube, CubeList from iris.pandas import as_data_frame from quantile_forest import RandomForestQuantileRegressor @@ -23,7 +24,6 @@ ) from improver.ensemble_copula_coupling.ensemble_copula_coupling import ( RebadgePercentilesAsRealizations, - RebadgeRealizationsAsPercentiles, ) from improver.ensemble_copula_coupling.utilities import choose_set_of_percentiles from improver.utilities.cube_checker import assert_spatial_coords_match @@ -37,7 +37,7 @@ class LoadAndApplyQRF(PostProcessingPlugin): def __init__( self, feature_config: dict[str, list[str]], - target_cube_name: str, + target_cf_name: str, transformation: Optional[str] = None, pre_transform_addition: Optional[float] = None, ): @@ -55,16 +55,16 @@ def __init__( water cube, then the value should state "static". This will ensure the cube's data is used as the feature. The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] "distance_to_water": ["static"], } - target_cube_name (str): - A string containing the cube name of the forecast to be - calibrated. This will be used to separate it from the rest of the - dynamic predictors, if present. + target_cf_name (str): + A string containing the CF name of diagnostic to be calibrated. This + will be used to separate it from the rest of the dynamic predictors, + if present. transformation (str): Transformation to be applied to the data before fitting. pre_transform_addition (float): @@ -72,7 +72,7 @@ def __init__( """ self.feature_config = feature_config - self.target_cube_name = target_cube_name + self.target_cf_name = target_cf_name self.transformation = transformation self.pre_transform_addition = pre_transform_addition @@ -113,7 +113,7 @@ def _get_inputs( # Extract all additional cubes which are associated with a feature in the # feature_config. - forecast_constraint = iris.Constraint(name=self.target_cube_name) + forecast_constraint = iris.Constraint(name=self.target_cf_name) forecast_cube = cube_inputs.extract(forecast_constraint) if forecast_cube: @@ -122,7 +122,7 @@ def _get_inputs( msg = ( "No target forecast provided. An input file representing the target " "must be provided, even if the target will not be used as a feature. " - f"The target is '{self.target_cube_name}'." + f"The target is '{self.target_cf_name}'." ) raise ValueError(msg) @@ -144,7 +144,7 @@ def _get_inputs( return forecast_cube # If target diagnostic not a feature in the training then remove. - if self.target_cube_name not in self.feature_config.keys(): + if self.target_cf_name not in self.feature_config.keys(): cube_inputs.remove(forecast_cube) return cube_inputs, forecast_cube, qrf_model @@ -191,6 +191,42 @@ def _percentiles_to_realizations(cube_inputs: CubeList) -> CubeList: realization_cube_inputs.append(feature_cube) return realization_cube_inputs + @staticmethod + def _cube_to_dataframe(cube_inputs: CubeList) -> pd.DataFrame: + """Convert cube inputs to a pandas DataFrame. + + Args: + cube_inputs: List of cubes containing the features and the forecast to be + calibrated. + Returns: + DataFrame containing the data from the cubes, with auxiliary coordinates + included as columns. + """ + # Convert the first cube to a DataFrame. + df = as_data_frame(cube_inputs[0], add_aux_coords=True).reset_index() + + # Iteratively convert remaining cubes to DataFrame and merge. + for cube in cube_inputs[1:]: + temporary_df = as_data_frame(cube, add_aux_coords=True).reset_index() + possible_columns = [ + "wmo_id", + "time", + "forecast_reference_time", + "forecast_period", + ] + merge_columns = [ + col for col in possible_columns if col in temporary_df.columns + ] + df = df.merge( + temporary_df[merge_columns + [cube.name()]], + on=merge_columns, + how="left", + ) + + for column in ["forecast_reference_time", "time"]: + df[column] = df[column].apply(lambda x: x._to_real_datetime()) + return df + def process( self, file_paths: list[pathlib.Path], @@ -219,48 +255,21 @@ def process( template_forecast_cube = forecast_cube.copy() if forecast_cube.coords("percentile"): percentiles = self._compute_percentiles(forecast_cube.copy(), "percentile") - cube_inputs = self._percentiles_to_realizations(cube_inputs.copy()) - template_forecast_cube = RebadgePercentilesAsRealizations()( - template_forecast_cube - ) elif forecast_cube.coords("realization"): percentiles = self._compute_percentiles(forecast_cube.copy(), "realization") assert_spatial_coords_match(cube_inputs) - - df = as_data_frame(cube_inputs[0], add_aux_coords=True).reset_index() - - for cube in cube_inputs[1:]: - temporary_df = as_data_frame(cube, add_aux_coords=True).reset_index() - possible_columns = [ - "wmo_id", - "time", - "forecast_reference_time", - "forecast_period", - ] - merge_columns = [ - col for col in possible_columns if col in temporary_df.columns - ] - df = df.merge( - temporary_df[merge_columns + [cube.name()]], - on=merge_columns, - how="left", - ) + df = self._cube_to_dataframe(cube_inputs) calibrated_forecast = ApplyQuantileRegressionRandomForests( - target_name=self.target_cube_name, + target_name=self.target_cf_name, feature_config=self.feature_config, quantiles=percentiles, transformation=self.transformation, pre_transform_addition=self.pre_transform_addition, )(qrf_model, df) - calibrated_forecast_cube = template_forecast_cube.copy( data=np.broadcast_to(calibrated_forecast.T, template_forecast_cube.shape) ) - if forecast_cube.coords("percentile"): - calibrated_forecast_cube = RebadgeRealizationsAsPercentiles()( - calibrated_forecast_cube - ) return calibrated_forecast_cube diff --git a/improver/calibration/load_and_train_quantile_regression_random_forest.py b/improver/calibration/load_and_train_quantile_regression_random_forest.py index 32085cd0de..6167380d5e 100644 --- a/improver/calibration/load_and_train_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_train_quantile_regression_random_forest.py @@ -32,6 +32,7 @@ def __init__( self, feature_config: dict[str, list[str]], target_diagnostic_name: str, + target_cf_name: str, forecast_periods: str, cycletime: str, training_length: int, @@ -46,6 +47,7 @@ def __init__( """Initialise the LoadAndTrainQRF plugin.""" self.feature_config = feature_config self.target_diagnostic_name = target_diagnostic_name + self.target_cf_name = target_cf_name self.forecast_periods = forecast_periods self.cycletime = cycletime self.training_length = training_length @@ -217,9 +219,19 @@ def _read_parquet_files( raise IOError(msg) return forecast_df, truth_df + @staticmethod def _check_matching_times( - self, forecast_df: pd.DataFrame, truth_df: pd.DataFrame + forecast_df: pd.DataFrame, truth_df: pd.DataFrame ) -> list[pd.Timestamp]: + """Find the intersecting times available within the forecast and truth + DataFrames. + + Args: + forecast_df: DataFrame containing the forecast data. + truth_df: DataFrame containing the truth data. + Returns: + List of intersecting times as pandas Timestamp objects. + """ return list(set(forecast_df["time"]).intersection(set(truth_df["time"]))) def _add_features_to_df( @@ -262,14 +274,6 @@ def filter_bad_sites( - DataFrame containing the forecast data with bad sites removed. - DataFrame containing the truth data with bad sites removed. """ - # import pdb - # pdb.set_trace() - # for coord in ["latitude", "longitude", "altitude", "ob_value"]: - # truth_df = truth_df.groupby("wmo_id").filter( - # lambda x: ~(x[coord].isna().any()) - # ) - # import pdb - # pdb.set_trace() truth_df.dropna( subset=["latitude", "longitude", "altitude", "ob_value"], inplace=True ) @@ -311,7 +315,7 @@ def process( water cube, then the value should state "static". This will ensure the cube's data is used as the feature. The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] @@ -376,12 +380,8 @@ def process( forecast_df, truth_df = self._read_parquet_files( forecast_table_path, truth_table_path, forecast_periods ) - forecast_df = forecast_df[forecast_df["experiment"] == self.experiment] - forecast_df = forecast_df.rename( - columns={"forecast": forecast_df["cf_name"][0]} - ) - # forecast_df = forecast_df.drop(columns=["cf_name", "diagnostic"]) + forecast_df = forecast_df.rename(columns={"forecast": self.target_cf_name}) intersecting_times = self._check_matching_times(forecast_df, truth_df) if len(intersecting_times) == 0: return None @@ -390,7 +390,7 @@ def process( forecast_df, truth_df = self.filter_bad_sites(forecast_df, truth_df) TrainQuantileRegressionRandomForests( - target_name=forecast_df["cf_name"][0], + target_name=self.target_cf_name, feature_config=self.feature_config, n_estimators=self.n_estimators, max_depth=self.max_depth, diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index 4508d2be97..593caf3992 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -9,7 +9,6 @@ import joblib import numpy as np import pandas as pd -from iris.cube import Cube from quantile_forest import RandomForestQuantileRegressor from improver import BasePlugin, PostProcessingPlugin @@ -18,70 +17,102 @@ def prep_feature( df: pd.DataFrame, + variable_name: str, feature_name: str, - feature: str, ) -> pd.DataFrame: - if feature in ["mean", "std"]: + """Prepare features that require computation from the input DataFrame. Options + available are mean and standard deviation of the input feature, the + day of year, sine of day of year, cosine of day of year, hour of day, + sine of hour of day and cosine of hour of day. When computing the mean or standard + deviation, these will be computed over either the percentile or realization column, + depending upon which is available. + + Args: + df: Input DataFrame. + variable_name: Name of the variable to be used for the computation. + feature_name: Feature to be computed. Options are "mean", "std", "day_of_year", + "day_of_year_sin", "day_of_year_cos", "hour_of_day", + "hour_of_day_sin" and "hour_of_day_cos". + Returns: + df: DataFrame with the computed feature added. + """ + if feature_name in ["mean", "std"]: representation_name = [ n for n in ["percentile", "realization"] if n in df.columns ][0] - subset_cols = [ - "forecast_reference_time", - "forecast_period", - "wmo_id", + groupby_cols = ["forecast_reference_time", "forecast_period", "wmo_id"] + subset_cols = [*groupby_cols] + [ representation_name, - feature_name, + variable_name, ] - groupby_cols = ["forecast_reference_time", "forecast_period", "wmo_id"] - if feature == "mean": + # For a subset of the input DataFrame compute the mean or standard deviation + # over the representation column, grouped by the groupby columns. + if feature_name == "mean": subset_df = df[subset_cols].groupby(groupby_cols).mean() - elif feature == "std": + elif feature_name == "std": subset_df = df[subset_cols].groupby(groupby_cols).std() subset_df = subset_df.reset_index() + # Rename the column to distinguish the computed feature from the original. subset_df.rename( - columns={feature_name: f"{feature_name}_{feature}"}, inplace=True + columns={variable_name: f"{variable_name}_{feature_name}"}, inplace=True ) - # df.drop(columns=[representation_name], inplace=True) - # df.drop_duplicates(inplace=True) + # Merge the computed feature back into the original DataFrame. df = df.merge( - subset_df[groupby_cols + [f"{feature_name}_{feature}"]], - on=["forecast_reference_time", "forecast_period", "wmo_id"], + subset_df[groupby_cols + [f"{variable_name}_{feature_name}"]], + on=groupby_cols, how="left", ) - elif feature in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]: - day_of_year = df["time"].dt.strftime("%j") - day_of_year = np.array(day_of_year, dtype=np.int32) - if feature == "day_of_year": + elif feature_name in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]: + # If there is only one unique time, compute the day of year for the first row. + if df["time"].nunique() == 1: + day_of_year = np.int32(df["time"].iloc[0].strftime("%j")) + else: + day_of_year = np.array(df["time"].dt.strftime("%j"), dtype=np.int32) + if feature_name == "day_of_year": feature_values = day_of_year - elif feature == "day_of_year_sin": + elif feature_name == "day_of_year_sin": feature_values = np.sin( 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) ).astype(np.float32) - elif feature == "day_of_year_cos": + elif feature_name == "day_of_year_cos": feature_values = np.cos( 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) ).astype(np.float32) - df[feature] = feature_values - elif feature in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]: - hour_of_day = df["time"].dt.hour - hour_of_day = np.array(hour_of_day, dtype=np.int32) - if feature == "hour_of_day": + df[feature_name] = feature_values + elif feature_name in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]: + if df["time"].nunique() == 1: + hour_of_day = np.int32(df["time"].iloc[0].hour) + else: + hour_of_day = np.array(df["time"].dt.hour, dtype=np.int32) + if feature_name == "hour_of_day": feature_values = hour_of_day - elif feature == "hour_of_day_sin": + elif feature_name == "hour_of_day_sin": feature_values = np.sin(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype( np.float32 ) - elif feature == "hour_of_day_cos": + elif feature_name == "hour_of_day_cos": feature_values = np.cos(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype( np.float32 ) - df[feature] = feature_values + df[feature_name] = feature_values return df -def sanitise_forecast_dataframe(df: pd.DataFrame, feature_config: dict) -> pd.DataFrame: +def sanitise_forecast_dataframe( + df: pd.DataFrame, feature_config: dict[str, list[str]] +) -> pd.DataFrame: + """Sanitise the forecast DataFrame by removing columns that are no longer + required. Following the computation of e.g. the mean or standard deviation, + the original feature can be removed. The column over which the mean or + standard deviation has been computed (e.g. the percentile or realization column) + is also removed. + + Args: + df: Input DataFrame, potentially including some computed features. + feature_config: Feature configuration defining the features to be used for QRF. + """ representation_name = [n for n in ["percentile", "realization"] if n in df.columns][ 0 ] @@ -94,12 +125,21 @@ def sanitise_forecast_dataframe(df: pd.DataFrame, feature_config: dict) -> pd.Da return df -def select_features(df, feature_config): - """Prepare the features from a DataFrame for quantile regression random forest.""" - - # Remove columns that are no longer required e.g. as the mean or std have been computed. - # Removing these columns allows a lot of duplicates to then be removed. +def get_required_column_names( + df: pd.DataFrame, feature_config: dict[str, list[str]] +) -> list[str]: + """Process the feature_config to return the expected column names that will be + used as features with the QRF. + Args: + df: Input DataFrame. + feature_config: Feature configuration defining the features to be used for QRF. + Returns: + List of expected column names that will be used as features with the QRF. + Raises: + ValueError: If a feature expected in the feature_config is not present in + the DataFrame. + """ feature_column_names = [] for feature_name in feature_config.keys(): for feature in feature_config[feature_name]: @@ -162,7 +202,7 @@ def __init__( should state "static". In this case, the name of feature e.g. 'distance_to_water' is expected to be a column name in the input dataframe. The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] @@ -259,15 +299,20 @@ def process( truth_df["ob_value"] + self.pre_transform_addition ) - for feature_name in self.feature_config.keys(): - if feature_name not in forecast_df.columns: - msg = f"Feature '{feature_name}' is not present in the forecast DataFrame." + for variable_name in self.feature_config.keys(): + if variable_name not in forecast_df.columns: + msg = ( + f"Feature '{variable_name}' is not present in the " + "forecast DataFrame." + ) raise ValueError(msg) - for feature in self.feature_config[feature_name]: - forecast_df = prep_feature(forecast_df, feature_name, feature) + for feature_name in self.feature_config[variable_name]: + forecast_df = prep_feature(forecast_df, variable_name, feature_name) forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config) - feature_column_names = select_features(forecast_df, self.feature_config) + feature_column_names = get_required_column_names( + forecast_df, self.feature_config + ) merge_columns = ["wmo_id", "time"] combined_df = forecast_df.merge( @@ -276,7 +321,6 @@ def process( feature_values = np.array(combined_df[feature_column_names]) target_values = combined_df["ob_value"].values - print(feature_values, target_values) # Fit the quantile regression model qrf_model = self.fit_qrf(feature_values, target_values) @@ -308,14 +352,14 @@ def __init__( should state "static". In this case, the name of feature e.g. 'distance_to_water' is expected to be a column name in the input dataframe. The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] "distance_to_water": ["static"], } quantiles (float): - Quantiles used for prediction (values ranging from 0 to 1) + Quantiles used for prediction (values ranging from 0 to 1). transformation (str): Transformation to be applied to the data before fitting. pre_transform_addition (float): @@ -333,10 +377,11 @@ def __init__( def _reverse_transformation(self, forecast: np.ndarray) -> np.ndarray: """Reverse the transformation applied to the data prior to fitting the QRF. - The forecast cube provided is modified in place. Args: - forecast_cube: Forecast to be calibrated. + forecast: Calibrated forecast. + Returns: + forecast: Forecast with the transformation reversed. """ if self.transformation: if self.transformation == "log": @@ -353,41 +398,39 @@ def process( self, qrf_model: RandomForestQuantileRegressor, forecast_df: pd.DataFrame, - ) -> Cube: + ) -> np.ndarray: """Apply a quantile regression random forests model. Args: qrf_model: A trained QRF model. - feature_cubes: CubeList of features. This should include the forecast to be - calibrated, if that has been used as a feature in the training, and - any features that can be provided as cubes. - template_forecast_cube: Template forecast cube that provides the required - metadata and shape for the output cube. The data from this cube will not - be used. + forecast_df: DataFrame containing the forecast information and features. Returns: - Calibrated forecast cube with the same metadata as the template - forecast cube. + Calibrated forecast as a numpy array. """ feature_values = [] - for feature_name in self.feature_config.keys(): + for variable_name in self.feature_config.keys(): # Transform the feature cube data if a transformation is specified. if ( self.transformation - and set(["mean", "std"]).intersection(self.feature_config[feature_name]) + and set(["mean", "std"]).intersection( + self.feature_config[variable_name] + ) and self.target_name in forecast_df.columns ): forecast_df[self.target_name] = getattr(np, self.transformation)( forecast_df[self.target_name] + self.pre_transform_addition ) - for feature in self.feature_config[feature_name]: - forecast_df = prep_feature(forecast_df, feature_name, feature) + for feature_name in self.feature_config[variable_name]: + forecast_df = prep_feature(forecast_df, variable_name, feature_name) forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config) - feature_column_names = select_features(forecast_df, self.feature_config) + feature_column_names = get_required_column_names( + forecast_df, self.feature_config + ) feature_values = np.array(forecast_df[feature_column_names]) @@ -397,5 +440,4 @@ def process( calibrated_forecast = np.float32(calibrated_forecast) calibrated_forecast = self._reverse_transformation(calibrated_forecast) - return calibrated_forecast diff --git a/improver/cli/apply_quantile_regression_random_forest.py b/improver/cli/apply_quantile_regression_random_forest.py index 92ad7a5fb6..a6cf4cd30f 100644 --- a/improver/cli/apply_quantile_regression_random_forest.py +++ b/improver/cli/apply_quantile_regression_random_forest.py @@ -13,7 +13,7 @@ def process( *file_paths: cli.inputpath, feature_config: cli.inputjson, - target_cube_name: str, + target_cf_name: str, transformation: str = None, pre_transform_addition: float = 0, ): @@ -32,23 +32,23 @@ def process( - The path to a NetCDF file containing the forecast to be calibrated. - Optionally, paths to NetCDF files containing additional predictors. feature_config (dict): - Feature configuration defining the features to be used for quantile regression. - The configuration is a dictionary of strings, where the keys are the names of - the input cube(s) supplied, and the values are a list. This list can contain both - computed features, such as the mean or standard deviation (std), or static - features, such as the altitude. The computed features will be computed using - the cube defined in the dictionary key. If the key is the feature itself e.g. - a distance to water cube, then the value should state "static". This will ensure - the cube's data is used as the feature. - The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + Feature configuration defining the features to be used for quantile + regression. The configuration is a dictionary of strings, where the keys + are the names of the input cube(s) supplied, and the values are a list. + This list can contain both computed features, such as the mean or + standard deviation (std), or static features, such as the altitude. The + computed features will be computed using the cube defined in the + dictionary key. If the key is the feature itself e.g. a distance to water + cube, then the value should state "static". This will ensure the cube's + data is used as the feature. The config will have the structure: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] "distance_to_water": ["static"], } - target_cube_name (str): - A string containing the cube name of the forecast to be + target_cf_name (str): + A string containing the CF name of the forecast to be calibrated. This will be used to separate it from the rest of the feature cubes, if present. transformation (str): @@ -66,7 +66,7 @@ def process( result = LoadAndApplyQRF( feature_config=feature_config, - target_cube_name=target_cube_name, + target_cf_name=target_cf_name, transformation=transformation, pre_transform_addition=pre_transform_addition, )( diff --git a/improver/cli/train_quantile_regression_random_forest.py b/improver/cli/train_quantile_regression_random_forest.py index 27ea3cf586..ec0badf9e1 100644 --- a/improver/cli/train_quantile_regression_random_forest.py +++ b/improver/cli/train_quantile_regression_random_forest.py @@ -13,6 +13,7 @@ def process( *file_paths: cli.inputpath, feature_config: cli.inputjson, target_diagnostic_name: str, + target_cf_name: str, forecast_periods: str, cycletime: str, training_length: int, @@ -52,7 +53,7 @@ def process( If the key is the feature itself e.g. a distance to water cube, then the value should state "static". This will ensure the cube's data is used as the feature. The config will have the structure: - "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g: + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: { "air_temperature": ["mean", "std", "altitude"], "visibility_at_screen_level": ["mean", "std"] @@ -62,6 +63,8 @@ def process( A string containing the diagnostic name of the forecast to be calibrated. This will be used to filter the target forecast and truth dataframes. + target_cf_name (str): + A string containing the CF name of the forecast to be calibrated. forecast_periods (str): Range of forecast periods to be calibrated in hours in the form: "start:end:interval" e.g. "6:18:6" or a single forecast period e.g. "6". @@ -100,6 +103,7 @@ def process( experiment=experiment, feature_config=feature_config, target_diagnostic_name=target_diagnostic_name, + target_cf_name=target_cf_name, forecast_periods=forecast_periods, cycletime=cycletime, training_length=training_length, diff --git a/improver_tests/acceptance/acceptance.py b/improver_tests/acceptance/acceptance.py index 0cf644c71f..3616f225be 100644 --- a/improver_tests/acceptance/acceptance.py +++ b/improver_tests/acceptance/acceptance.py @@ -320,6 +320,8 @@ def compare( rtol (float): Relative tolerance exclude_vars (Iterable[str]): Variables to exclude from comparison exclude_attributes (Iterable[str]): Attributes to exclude from comparison + file_type (str): Name of file type to compare, either "netCDF" or + "pickled_forest". Returns: None diff --git a/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py b/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py index 68e8d9cb40..5cd27f5f17 100644 --- a/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py +++ b/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py @@ -30,7 +30,7 @@ def test_basic(tmp_path, transformation): qrf_path, "--feature-config", config_path, - "--target-cube-name", + "--target-cf-name", "air_temperature", "--output", output_path, diff --git a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py index 091a30a0b3..64feccd97e 100644 --- a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py +++ b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py @@ -43,6 +43,8 @@ def test_basic( config_path, "--target-diagnostic-name", "temperature_at_screen_level", + "--target-cf-name", + "air_temperature", "--forecast-periods", "6:18:6", "--cycletime", diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py index f4642019a2..c94d310922 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py @@ -403,6 +403,7 @@ def test_load_and_train_qrf( experiment="latestblend", feature_config=feature_config, target_diagnostic_name="temperature_at_screen_level", + target_cf_name="air_temperature", forecast_periods=forecast_periods, cycletime="20170103T0000Z", training_length=2, @@ -459,6 +460,7 @@ def test_load_and_train_qrf_no_paths(tmp_path, make_files): experiment="latestblend", feature_config=feature_config, target_diagnostic_name="temperature_at_screen_level", + target_cf_name="air_temperature", forecast_periods="6:12:6", cycletime="20170102T0000Z", training_length=2, @@ -503,6 +505,7 @@ def test_load_and_train_qrf_mismatches(tmp_path, cycletime, forecast_periods): experiment="latestblend", feature_config=feature_config, target_diagnostic_name="temperature_at_screen_level", + target_cf_name="air_temperature", forecast_periods=forecast_periods, cycletime=cycletime, training_length=2, @@ -586,6 +589,7 @@ def test_exceptions( experiment="latestblend", feature_config=feature_config, target_diagnostic_name="temperature_at_screen_level", + target_cf_name="air_temperature", forecast_periods=forecast_periods, cycletime="20170103T0000Z", training_length=2, diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py index 15eef7616b..74f0d03785 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py @@ -17,8 +17,8 @@ from improver.calibration.quantile_regression_random_forest import ( ApplyQuantileRegressionRandomForests, TrainQuantileRegressionRandomForests, + get_required_column_names, prep_feature, - select_features, ) from improver.metadata.constants.time_types import DT_FORMAT from improver.synthetic_data.set_up_test_cubes import set_up_spot_variable_cube @@ -405,7 +405,7 @@ def test_prep_feature_more_times(feature, expected, expected_dtype): {"wind_speed_at_10m": ["latitude", "longitude", "height"]}, ], ) -def test_select_features(feature_config): +def test_get_required_column_names(feature_config): data_dict = { "wmo_id": np.tile(WMO_ID, 3), "latitude": np.tile(LATITUDE, 3), @@ -431,9 +431,9 @@ def test_select_features(feature_config): if "height" in feature_config["wind_speed_at_10m"]: with pytest.raises(ValueError, match="Feature 'height' is not supported."): - select_features(df, feature_config) + get_required_column_names(df, feature_config) else: - result = select_features(df, feature_config) + result = get_required_column_names(df, feature_config) assert len(result) == len(expected) assert result == expected From fb89f58e2b1845a0b7e8322000a0fa5333e330b8 Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Wed, 20 Aug 2025 14:47:11 +0100 Subject: [PATCH 5/8] Further tidy-ups and efficiency improvements to reduce runtimes. --- ...train_quantile_regression_random_forest.py | 151 +++++++----------- .../quantile_regression_random_forest.py | 56 ++++--- ...train_quantile_regression_random_forest.py | 7 + ...train_quantile_regression_random_forest.py | 13 +- .../test_quantile_regression_random_forest.py | 145 ++++++++++++----- 5 files changed, 205 insertions(+), 167 deletions(-) diff --git a/improver/calibration/load_and_train_quantile_regression_random_forest.py b/improver/calibration/load_and_train_quantile_regression_random_forest.py index 6167380d5e..e170c2d7b7 100644 --- a/improver/calibration/load_and_train_quantile_regression_random_forest.py +++ b/improver/calibration/load_and_train_quantile_regression_random_forest.py @@ -10,6 +10,7 @@ from typing import Optional import iris +import numpy as np import pandas as pd import pyarrow as pa import pyarrow.parquet as pq @@ -36,11 +37,12 @@ def __init__( forecast_periods: str, cycletime: str, training_length: int, - experiment: str = None, + experiment: Optional[str] = None, n_estimators: int = 100, - max_depth: int = None, - random_state: int = None, - transformation: str = None, + max_depth: Optional[int] = None, + max_samples: Optional[float] = None, + random_state: Optional[int] = None, + transformation: Optional[str] = None, pre_transform_addition: float = 0, compression: int = 5, ): @@ -54,6 +56,7 @@ def __init__( self.experiment = experiment self.n_estimators = n_estimators self.max_depth = max_depth + self.max_samples = max_samples self.random_state = random_state self.transformation = transformation self.pre_transform_addition = pre_transform_addition @@ -84,24 +87,24 @@ def _split_cubes_and_parquet_files( # file extraction loop: for file_path in file_paths: - try: - cube = load_cube(str(file_path)) - cube_inputs.append(cube) - except IsADirectoryError: - # For loop here because the read_schema must read a .parquet file - # rather than a directory. - for file in Path(file_path).glob("**/*.parquet"): - try: - pq.read_schema(file).field("forecast_period") - forecast_table_path = file_path - except KeyError: - truth_table_path = file_path - if forecast_table_path and truth_table_path: - break - except OSError: + if not Path(file_path).exists(): # This will occur when the filepath does not exist. In this case, # return None. return None, None, None + elif Path(file_path).is_dir(): + try: + example_file_path = next(Path(file_path).glob("**/*.parquet")) + except StopIteration: + # If no parquet files are found, return None. + return None, None, None + try: + pq.read_schema(example_file_path).field("forecast_period") + forecast_table_path = file_path + except KeyError: + truth_table_path = file_path + else: + cube = load_cube(str(file_path)) + cube_inputs.append(cube) if len(self.feature_config.keys()) not in [ len(cube_inputs), @@ -143,7 +146,6 @@ def _read_parquet_files( fields. """ cycletimes = [] - for forecast_period in forecast_periods: # Load forecasts from parquet file filtering by diagnostic and blend_time. forecast_period_td = pd.Timedelta(int(forecast_period), unit="seconds") @@ -167,23 +169,20 @@ def _read_parquet_files( ] ] - for file in Path(forecast_table_path).glob("**/*.parquet"): - if pq.read_schema(file).get_all_field_indices("percentile"): - altered_schema = FORECAST_SCHEMA - elif pq.read_schema(file).get_all_field_indices("realization"): - altered_schema = FORECAST_SCHEMA.remove( - FORECAST_SCHEMA.get_field_index("percentile") - ) - altered_schema = altered_schema.append( - pa.field("realization", pa.int64()) - ) - else: - msg = ( - "The forecast parquet file is expected to contain either a " - "'percentile' or 'realization' field. Neither was found." - ) - raise ValueError(msg) - break + example_file_path = next(Path(forecast_table_path).glob("**/*.parquet")) + if pq.read_schema(example_file_path).get_all_field_indices("percentile"): + altered_schema = FORECAST_SCHEMA + elif pq.read_schema(example_file_path).get_all_field_indices("realization"): + altered_schema = FORECAST_SCHEMA.remove( + FORECAST_SCHEMA.get_field_index("percentile") + ) + altered_schema = altered_schema.append(pa.field("realization", pa.int64())) + else: + msg = ( + "The forecast parquet file is expected to contain either a " + "'percentile' or 'realization' field. Neither was found." + ) + raise ValueError(msg) forecast_df = pd.read_parquet( forecast_table_path, @@ -192,16 +191,18 @@ def _read_parquet_files( engine="pyarrow", ) - # Convert df columns from ms to pandas timestamp object to work with existing - # code + forecast_df = forecast_df[ + forecast_df["forecast_period"].isin(np.array(forecast_periods)* 1e9) + ] + # Convert df columns from ns to pandas timestamp object. for column in ["time", "forecast_reference_time", "blend_time"]: forecast_df[column] = pd.to_datetime( forecast_df[column], unit="ns", utc=True ) - forecast_df["forecast_period"] = pd.to_timedelta( - forecast_df["forecast_period"], unit="ns" - ) - forecast_df["period"] = pd.to_timedelta(forecast_df["period"], unit="ns") + for column in ["forecast_period", "period"]: + forecast_df[column] = pd.to_timedelta(forecast_df[column], unit="ns") + + forecast_df = forecast_df.rename(columns={"forecast": self.target_cf_name}) # Load truths from parquet file filtering by diagnostic. filters = [[("diagnostic", "==", self.target_diagnostic_name)]] @@ -232,7 +233,13 @@ def _check_matching_times( Returns: List of intersecting times as pandas Timestamp objects. """ - return list(set(forecast_df["time"]).intersection(set(truth_df["time"]))) + # Calling unique() on the time column is quicker than relying upon set() to + # find the unique times. + return list( + set(forecast_df["time"].unique()).intersection( + set(truth_df["time"].unique()) + ) + ) def _add_features_to_df( self, forecast_df: pd.DataFrame, cube_inputs: iris.cube.CubeList @@ -278,7 +285,9 @@ def filter_bad_sites( subset=["latitude", "longitude", "altitude", "ob_value"], inplace=True ) - wmo_ids = set(forecast_df["wmo_id"]).intersection(set(truth_df["wmo_id"])) + wmo_ids = set(forecast_df["wmo_id"]).intersection( + set(truth_df["wmo_id"]) + ) forecast_df = forecast_df[forecast_df["wmo_id"].isin(wmo_ids)] truth_df = truth_df[truth_df["wmo_id"].isin(wmo_ids)] @@ -304,55 +313,6 @@ def process( - The path to a Parquet file containing the forecasts to be used for calibration. - Optionally, paths to NetCDF files containing additional predictors. - feature_config (dict): - Feature configuration defining the features to be used for quantile - regression. The configuration is a dictionary of strings, where the - keys are the names of the input cube(s) supplied, and the values are - a list. This list can contain both computed features, such as the mean - or standard deviation (std), or static features, such as the altitude. - The computed features will be computed using the cube defined in the - dictionary key. If the key is the feature itself e.g. a distance to - water cube, then the value should state "static". This will ensure - the cube's data is used as the feature. The config will have the - structure: - "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: - { - "air_temperature": ["mean", "std", "altitude"], - "visibility_at_screen_level": ["mean", "std"] - "distance_to_water": ["static"], - } - target_diagnostic_name (str): - A string containing the diagnostic name of the forecast to be - calibrated. This will be used to filter the target forecast and truth - dataframes. - forecast_period (int): - Range of forecast periods to be calibrated in hours in the form: - "start:end:interval" e.g. "6:18:6" or a single forecast period e.g. "6". - The end value is exclusive, so "6:18:6" will calibrate the 6 and 12 - hours. - cycletime (str): - Cycletime of the forecast to be calibrated in a format similar to - 20170109T0000Z. This is used to filter the correct blendtimes from - the dataframe on load. - training_length (int): - The length of the training period in days. - experiment (str): - The name of the experiment (step) that calibration is applied to. - n_estimators (int): - Number of trees in the forest. - max_depth (int): - Maximum depth of the tree. - random_state (int): - Random seed for reproducibility. - transformation (str): - Transformation to be applied to the data before fitting. - Supported transformations are "log", "log10", "sqrt", and "cbrt". - pre_transform_addition (float): - Value to be added before transformation. This is useful for ensuring - that the data is positive before applying a transformation such as log. - compression (int): - Compression level for saving the model. Please see the joblib - documentation for more information on compression levels. model_output (str): Full path including model file name that will store the pickled model. @@ -380,8 +340,6 @@ def process( forecast_df, truth_df = self._read_parquet_files( forecast_table_path, truth_table_path, forecast_periods ) - forecast_df = forecast_df[forecast_df["experiment"] == self.experiment] - forecast_df = forecast_df.rename(columns={"forecast": self.target_cf_name}) intersecting_times = self._check_matching_times(forecast_df, truth_df) if len(intersecting_times) == 0: return None @@ -394,6 +352,7 @@ def process( feature_config=self.feature_config, n_estimators=self.n_estimators, max_depth=self.max_depth, + max_samples=self.max_samples, random_state=self.random_state, transformation=self.transformation, pre_transform_addition=self.pre_transform_addition, diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index 593caf3992..69ad283a4f 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -65,23 +65,27 @@ def prep_feature( ) elif feature_name in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]: - # If there is only one unique time, compute the day of year for the first row. - if df["time"].nunique() == 1: - day_of_year = np.int32(df["time"].iloc[0].strftime("%j")) - else: - day_of_year = np.array(df["time"].dt.strftime("%j"), dtype=np.int32) - if feature_name == "day_of_year": - feature_values = day_of_year - elif feature_name == "day_of_year_sin": - feature_values = np.sin( - 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) + # For a large DataFrame, the strftime("%j") computation can take a noticeable + # amount of time, so this computation is done once for each unique time + # and then merged back into the DataFrame. + doy_df = pd.DataFrame({"time": df["time"].unique()}) + doy_df["day_of_year"] = np.array(doy_df["time"].dt.strftime("%j"), np.int32) + + if feature_name == "day_of_year_sin": + doy_df[feature_name] = np.sin( + 2 * np.pi * doy_df["day_of_year"].values / (DAYS_IN_YEAR + 1) ).astype(np.float32) elif feature_name == "day_of_year_cos": - feature_values = np.cos( - 2 * np.pi * day_of_year / (DAYS_IN_YEAR + 1) + doy_df[feature_name] = np.cos( + 2 * np.pi * doy_df["day_of_year"].values / (DAYS_IN_YEAR + 1) ).astype(np.float32) - df[feature_name] = feature_values + df = df.merge( + doy_df[["time", feature_name]], on="time", how="left" + ) elif feature_name in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]: + # For hour_of_day, unlike day_of_year, the hour attribute doesn't require + # computation, therefore there is no benefit to creating the separate DataFrame + # and merging it back into the DataFrame. if df["time"].nunique() == 1: hour_of_day = np.int32(df["time"].iloc[0].hour) else: @@ -120,8 +124,12 @@ def sanitise_forecast_dataframe( for key, values in feature_config.items(): collapsed_features.extend([key for v in values if v in ["mean", "std"]]) collapsed_features = list(set(collapsed_features)) - - df = df.drop(columns=[representation_name, *collapsed_features]).drop_duplicates() + # Subset the dataframe by the first value of the representation column + # and drop the representation column and any features where the original variable + # is no longer required. This reduces the size of the DataFrame e.g. if there are + # 3 percentiles initially, the subsetted dataframe will be 1/3 of the size. + df = df[df[representation_name] == df[representation_name].iloc[0]] + df = df.drop(columns=[representation_name, *collapsed_features]) return df @@ -141,12 +149,12 @@ def get_required_column_names( the DataFrame. """ feature_column_names = [] - for feature_name in feature_config.keys(): - for feature in feature_config[feature_name]: + for variable_name in feature_config.keys(): + for feature in feature_config[variable_name]: if feature in ["mean", "std"]: - feature_column_names.append(f"{feature_name}_{feature}") + feature_column_names.append(f"{variable_name}_{feature}") elif feature in ["static"]: - feature_column_names.append(feature_name) + feature_column_names.append(variable_name) else: feature_column_names.append(feature) @@ -181,6 +189,7 @@ def __init__( feature_config: dict[str, list[str]], n_estimators: int, max_depth: Optional[int] = None, + max_samples: Optional[float] = None, random_state: Optional[int] = None, transformation: Optional[str] = None, pre_transform_addition: np.float32 = 0, @@ -212,6 +221,11 @@ def __init__( Number of trees in the forest. max_depth (int): Maximum depth of the tree. + max_samples (float): + If an int, then it is the number of samples to draw to train + each tree. If a float, then it is the fraction of samples to draw + to train each tree. If None, then each tree contains the same + total number of samples as originally provided. random_state (int): Random seed for reproducibility. transformation (str): @@ -229,6 +243,7 @@ def __init__( self.feature_config = feature_config self.n_estimators = n_estimators self.max_depth = max_depth + self.max_samples = max_samples self.random_state = random_state self.transformation = transformation _check_valid_transformation(self.transformation) @@ -255,6 +270,7 @@ def fit_qrf( qrf_model = RandomForestQuantileRegressor( n_estimators=self.n_estimators, max_depth=self.max_depth, + max_samples=self.max_samples, random_state=self.random_state, **self.kwargs, ) @@ -290,7 +306,6 @@ def process( Nonlin. Processes Geophys., 27, 329–347, https://doi.org/10.5194/npg-27-329-2020, 2020. """ - if self.transformation: forecast_df[self.target_name] = getattr(np, self.transformation)( forecast_df[self.target_name] + self.pre_transform_addition @@ -313,7 +328,6 @@ def process( feature_column_names = get_required_column_names( forecast_df, self.feature_config ) - merge_columns = ["wmo_id", "time"] combined_df = forecast_df.merge( truth_df[merge_columns + ["ob_value"]], on=merge_columns, how="inner" diff --git a/improver/cli/train_quantile_regression_random_forest.py b/improver/cli/train_quantile_regression_random_forest.py index ec0badf9e1..2c10a44421 100644 --- a/improver/cli/train_quantile_regression_random_forest.py +++ b/improver/cli/train_quantile_regression_random_forest.py @@ -20,6 +20,7 @@ def process( experiment: str = None, n_estimators: int = 100, max_depth: int = None, + max_samples: float = None, random_state: int = None, transformation: str = None, pre_transform_addition: float = 0, @@ -80,6 +81,11 @@ def process( Number of trees in the forest. max_depth (int): Maximum depth of the tree. + max_samples (float): + If an int, then it is the number of samples to draw to train + each tree. If a float, then it is the fraction of samples to draw + to train each tree. If None, then each tree contains the same + total number of samples as originally provided. random_state (int): Random seed for reproducibility. transformation (str): @@ -109,6 +115,7 @@ def process( training_length=training_length, n_estimators=n_estimators, max_depth=max_depth, + max_samples=max_samples, random_state=random_state, transformation=transformation, pre_transform_addition=pre_transform_addition, diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py index c94d310922..0891ec6f82 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py @@ -60,9 +60,6 @@ def _create_multi_site_forecast_parquet_file(tmp_path, representation="percentil wind_speed_df = pd.DataFrame(wind_speed_dict) joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True) - joined_df["forecast_period"] = joined_df["forecast_period"].astype( - "timedelta64[ms]" - ) output_dir = tmp_path / "forecast_parquet_files" output_dir.mkdir(parents=True, exist_ok=True) output_path = str(output_dir / "forecast.parquet") @@ -101,9 +98,6 @@ def _create_multi_percentile_forecast_parquet_file(tmp_path, representation=None data_df = pd.DataFrame(data_dict) wind_speed_df = pd.DataFrame(wind_speed_dict) joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True) - joined_df["forecast_period"] = joined_df["forecast_period"].astype( - "timedelta64[ms]" - ) output_dir = tmp_path / "forecast_parquet_files" output_dir.mkdir(parents=True, exist_ok=True) @@ -152,9 +146,6 @@ def _create_multi_forecast_period_forecast_parquet_file(tmp_path, representation data_df = pd.DataFrame(data_dict) wind_speed_df = pd.DataFrame(wind_speed_dict) joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True) - joined_df["forecast_period"] = joined_df["forecast_period"].astype( - "timedelta64[ms]" - ) output_dir = tmp_path / "forecast_parquet_files" output_dir.mkdir(parents=True, exist_ok=True) @@ -353,7 +344,7 @@ def _create_ancil_file(tmp_path, wmo_ids): False, False, "realization", # Provide realization input - 5.6, + 5.62, ), ( _create_multi_forecast_period_forecast_parquet_file, @@ -362,7 +353,7 @@ def _create_ancil_file(tmp_path, wmo_ids): False, False, "percentile", - 5.61, + 5.64, ), ], ) diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py index 74f0d03785..d8fc7af795 100644 --- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py +++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py @@ -7,18 +7,22 @@ from datetime import datetime as dt import iris +import itertools import joblib import numpy as np import pandas as pd import pytest from iris.cube import Cube from iris.pandas import as_data_frame +from pandas.testing import assert_frame_equal from improver.calibration.quantile_regression_random_forest import ( ApplyQuantileRegressionRandomForests, TrainQuantileRegressionRandomForests, + _check_valid_transformation, get_required_column_names, prep_feature, + sanitise_forecast_dataframe, ) from improver.metadata.constants.time_types import DT_FORMAT from improver.synthetic_data.set_up_test_cubes import set_up_spot_variable_cube @@ -393,6 +397,51 @@ def test_prep_feature_more_times(feature, expected, expected_dtype): np.testing.assert_allclose(result[feature_name_modified], expected, atol=1e-6) +@pytest.mark.parametrize( + "feature_config", + [ + {"wind_speed_at_10m": ["mean", "std", "latitude", "longitude"]}, + {"wind_speed_at_10m": ["latitude", "longitude"]}, + { + "wind_speed_at_10m": ["latitude", "longitude"], + "distance_to_water": ["static"], + }, + {"wind_speed_at_10m": ["latitude", "longitude", "height"]}, + ], +) +def test_sanitise_forecast_dataframe(feature_config): + """Test sanitise_forecast_dataframe function.""" + data_dict = { + "wmo_id": np.tile(WMO_ID, 3), + "latitude": np.tile(LATITUDE, 3), + "longitude": np.tile(LONGITUDE, 3), + "altitude": np.tile(ALTITUDE, 3), + "wind_speed_at_10m_mean": np.repeat(5, 6), + "wind_speed_at_10m_std": np.repeat(1, 6), + "wind_speed_at_10m": np.tile([4, 6], 3), + "realization": np.tile([1, 2], 3), + "distance_to_water": np.tile([2.0, 3.0], 3), + } + df = pd.DataFrame(data_dict) + + expected = df.copy() + if ( + "mean" in feature_config["wind_speed_at_10m"] + or "std" in feature_config["wind_speed_at_10m"] + ): + expected.drop( + columns=[ + "wind_speed_at_10m", + ], + inplace=True, + ) + expected = expected[expected["realization"] == 1] + expected = expected.drop(columns=["realization"]) + + result = sanitise_forecast_dataframe(df, feature_config) + assert_frame_equal(result, expected) + + @pytest.mark.parametrize( "feature_config", [ @@ -406,6 +455,7 @@ def test_prep_feature_more_times(feature, expected, expected_dtype): ], ) def test_get_required_column_names(feature_config): + """Test the get_required_column_names function.""" data_dict = { "wmo_id": np.tile(WMO_ID, 3), "latitude": np.tile(LATITUDE, 3), @@ -419,25 +469,39 @@ def test_get_required_column_names(feature_config): } df = pd.DataFrame(data_dict) - expected = [] - for key, values in feature_config.items(): - for value in values: - if "mean" in value or "std" in value: - expected.append(f"wind_speed_at_10m_{value}") - elif value == "static": - expected.append(key) - else: - expected.append(value) - - if "height" in feature_config["wind_speed_at_10m"]: + # expected = feature_config["wind_speed_at_10m"].copy() + expected = list(itertools.chain.from_iterable(feature_config.values())) + if "mean" in expected: + expected = ["wind_speed_at_10m_mean" if e == "mean" else e for e in expected] + if "std" in expected: + expected = ["wind_speed_at_10m_std" if e == "std" else e for e in expected] + if "static" in expected: + expected = ["distance_to_water" if e == "static" else e for e in expected] + + if "height" in expected: with pytest.raises(ValueError, match="Feature 'height' is not supported."): get_required_column_names(df, feature_config) else: result = get_required_column_names(df, feature_config) - assert len(result) == len(expected) assert result == expected +@pytest.mark.parametrize( + "transformation", ["log", "log10", "sqrt", "cbrt", None, "yeojohnson"] +) +def test_check_valid_transformation(transformation): + """Test the _check_valid_transformation function.""" + + if transformation == "yeojohnson": + with pytest.raises( + ValueError, match="Currently the only supported transformations" + ): + _check_valid_transformation(transformation) + else: + result = _check_valid_transformation(transformation) + assert result is None + + @pytest.mark.parametrize( "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected", [ @@ -568,7 +632,7 @@ def test_train_qrf_multiple_lead_times( "feature_config,data,include_static,expected", [ ({"wind_speed_at_10m": ["mean"]}, [5], False, [5]), # One feature - ({"wind_speed_at_10m": ["latitude"]}, [61], False, [4.9]), # noqa Without the target + ({"wind_speed_at_10m": ["latitude"]}, [61], False, [7.75]), # noqa Without the target ({"wind_speed_at_10m": ["mean"]}, [5], True, [4]), # With static data ( {"wind_speed_at_10m": ["mean"], "air_temperature": ["mean"]}, @@ -650,31 +714,31 @@ def test_alternative_feature_configs( @pytest.mark.parametrize( "quantiles,transformation,pre_transform_addition,include_static,expected", [ - # ([0.5], None, 0, False, [5, 4.9]), # One quantile - ([0.1, 0.5, 0.9], None, 0, False, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.15]]), # noqa Multiple quantiles - # ([0.1, 0.5, 0.9], "log", 10, False, [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]]), # noqa Log transformation - # ( - # [0.1, 0.5, 0.9], - # "log10", - # 10, - # False, - # [[4.37, 4.37], [5.07, 5.07], [5.81, 5.81]], - # ), # Log10 transformation - # ( - # [0.1, 0.5, 0.9], - # "sqrt", - # 10, - # False, - # [[4.38, 4.38], [5.09, 5.09], [5.82, 5.82]], - # ), # Square root transformation - # ( - # [0.1, 0.5, 0.9], - # "cbrt", - # 10, - # False, - # [[4.37, 4.37], [5.08, 5.08], [5.81, 5.81]], - # ), # Cube root transformation - # ([0.1, 0.5, 0.9], None, 0, True, [[9.14, 9.14], [9.3, 9.3], [9.46, 9.46]]), # noqa Include static data + ([0.5], None, 0, False, [5, 4.9]), # One quantile + ([0.1, 0.5, 0.9], None, 0, False, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.14]]), # noqa Multiple quantiles + ([0.1, 0.5, 0.9], "log", 10, False, [[4.48, 5.67, 6.96], [4.48, 5.67, 6.96]]), # noqa Log transformation + ( + [0.1, 0.5, 0.9], + "log10", + 10, + False, + [[4.48, 5.67, 6.96], [4.48, 5.67, 6.96]], + ), # Log10 transformation + ( + [0.1, 0.5, 0.9], + "sqrt", + 10, + False, + [[4.5, 5.71, 6.98], [4.5, 5.71, 6.98]], + ), # Square root transformation + ( + [0.1, 0.5, 0.9], + "cbrt", + 10, + False, + [[4.49, 5.7, 6.97], [4.49, 5.7, 6.97]], + ), # Cube root transformation + ([0.1, 0.5, 0.9], None, 0, True, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.15]]), # noqa Include static data ], ) def test_apply_qrf( @@ -730,7 +794,10 @@ def test_apply_qrf( result = plugin.process(qrf_model, forecast_df) assert isinstance(result, np.ndarray) - assert result.shape == (2, 3) + if len(quantiles) == 3: + assert result.shape == (2, 3) + else: + assert result.shape == (2,) assert result.dtype == np.float32 # expected_data = np.full((len(quantiles), 2), expected, dtype=np.float32) np.testing.assert_almost_equal(result, expected, decimal=2) From b4c35ad3f20064bde946e6fb23cd014cc9909307 Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Wed, 20 Aug 2025 15:20:16 +0100 Subject: [PATCH 6/8] Update to SHAs. --- improver/utilities/compare.py | 2 +- improver_tests/acceptance/SHA256SUMS | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/improver/utilities/compare.py b/improver/utilities/compare.py index 99c8a1d9f1..88868fdfa1 100644 --- a/improver/utilities/compare.py +++ b/improver/utilities/compare.py @@ -415,7 +415,7 @@ def raise_reporter(message): assert output.random_state == kgo.random_state for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_): assert (output_estimator.tree_.value == kgo_estimator.tree_.value).all() - except AssertionError: + except (AssertionError, ValueError): difference_found = True # call the reporter function outside the except block to avoid nested # exceptions if the reporter function is raising an exception diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index b9250ace64..4757fb1a02 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -1003,8 +1003,8 @@ e3d2315b24bf769bcfb137033950d7100c8795878635e0ec9491413a1349904a ./train-quanti 1aec562fc0c39c7695521482b889d335a3becb10a76b7eae965447cbd9faf051 ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/0600Z_0.parquet e81c04f9f2d8dc14b8d8cd9bc9e8827ebd0da5ba016db48cf60ab738d327f968 ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1200Z_0.parquet 99e6656982672ea2d4cb9fa1dfb5731408dd41b2c9052b7e0df8d8004864241e ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1800Z_0.parquet -9d486570ca562c63aa45ed5f7621ecb00395891af7e021e382219cce74d484e8 ./train-quantile-regression-random-forest/with_transformation_kgo.pickle -f3919673f39a7db108add9aa5e4edd7706b6fd03a06e45f9c1d89070740524d4 ./train-quantile-regression-random-forest/without_transformation_kgo.pickle +d634d639dda9a4a2789eed69bfdec59954ce63f51734a71f1d9a20c3958d2786 ./train-quantile-regression-random-forest/with_transformation_kgo.pickle +e3f9fcd7544028436e34415db38d3c5841a6e22abbc3fc635e6ff87236764fcc ./train-quantile-regression-random-forest/without_transformation_kgo.pickle e2cfb5e19ef5ebcfd73dc1c8504b66e9a55ad76b7b18cb79318659224079f234 ./uv-index/basic/20181210T0600Z-PT0000H00M-radiation_flux_in_uv_downward_at_surface.nc e48c3b07bd14214a1a10c0bbf1e0acdf54cfedf93b2c6d766f594ce83a45c2b8 ./uv-index/basic/kgo.nc 2226a5e95eb29664e21f8c0658101a53d65945a1450a60a19955563a2693d22d ./vertical-updraught/cape.nc From edee0cc68546b0526a49393a94f9fbd58abafefd Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Tue, 26 Aug 2025 10:35:50 +0100 Subject: [PATCH 7/8] Docstring formatting. --- .../quantile_regression_random_forest.py | 30 +++++++++++-------- ...apply_quantile_regression_random_forest.py | 12 ++++---- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py index f284150f5c..cf43788b27 100644 --- a/improver/calibration/quantile_regression_random_forest.py +++ b/improver/calibration/quantile_regression_random_forest.py @@ -196,6 +196,7 @@ def __init__( **kwargs, ) -> None: """Initialise the plugin. + Args: target_name (str): Name of the target variable to be calibrated e.g. 'air_temperature'. @@ -209,12 +210,12 @@ def __init__( should state "static". In this case, the name of feature e.g. 'distance_to_water' is expected to be a column name in the input dataframe. The config will have the structure: - "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: - { - "air_temperature": ["mean", "std", "altitude"], - "visibility_at_screen_level": ["mean", "std"] - "distance_to_water": ["static"], - } + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g. + { + "air_temperature": ["mean", "std", "altitude"], + "visibility_at_screen_level": ["mean", "std"] + "distance_to_water": ["static"], + } n_estimators (int): Number of trees in the forest. max_depth (int): @@ -236,6 +237,7 @@ def __init__( Full path including model file name that will store the pickled model. kwargs: Additional keyword arguments for the quantile regression model. + """ self.target_name = target_name self.feature_config = feature_config @@ -281,6 +283,7 @@ def process( truth_df: pd.DataFrame, ) -> None: """Train a quantile regression random forests model. + Args: forecast_df: DataFrame containing the forecast information and features. @@ -303,6 +306,7 @@ def process( operational ensemble post-processing in France using machine learning, Nonlin. Processes Geophys., 27, 329–347, https://doi.org/10.5194/npg-27-329-2020, 2020. + """ if self.transformation: forecast_df[self.target_name] = getattr(np, self.transformation)( @@ -351,6 +355,7 @@ def __init__( pre_transform_addition: np.float32 = 0, ) -> None: """Initialise the plugin. + Args: target_name (str): Name of the target variable to be calibrated. @@ -364,12 +369,12 @@ def __init__( should state "static". In this case, the name of feature e.g. 'distance_to_water' is expected to be a column name in the input dataframe. The config will have the structure: - "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: - { - "air_temperature": ["mean", "std", "altitude"], - "visibility_at_screen_level": ["mean", "std"] - "distance_to_water": ["static"], - } + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g. + { + "air_temperature": ["mean", "std", "altitude"], + "visibility_at_screen_level": ["mean", "std"] + "distance_to_water": ["static"], + } quantiles (float): Quantiles used for prediction (values ranging from 0 to 1). transformation (str): @@ -379,6 +384,7 @@ def __init__( Raises: ValueError: If the transformation is not one of the supported types. + """ self.target_name = target_name self.feature_config = feature_config diff --git a/improver/cli/apply_quantile_regression_random_forest.py b/improver/cli/apply_quantile_regression_random_forest.py index a6cf4cd30f..c7a6fec8d4 100644 --- a/improver/cli/apply_quantile_regression_random_forest.py +++ b/improver/cli/apply_quantile_regression_random_forest.py @@ -41,12 +41,12 @@ def process( dictionary key. If the key is the feature itself e.g. a distance to water cube, then the value should state "static". This will ensure the cube's data is used as the feature. The config will have the structure: - "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g: - { - "air_temperature": ["mean", "std", "altitude"], - "visibility_at_screen_level": ["mean", "std"] - "distance_to_water": ["static"], - } + "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g. + { + "air_temperature": ["mean", "std", "altitude"], + "visibility_at_screen_level": ["mean", "std"] + "distance_to_water": ["static"], + } target_cf_name (str): A string containing the CF name of the forecast to be calibrated. This will be used to separate it from the rest of the From 7be86e7f9eb22387206501e2d50c1e94a5ac403a Mon Sep 17 00:00:00 2001 From: Gavin Evans Date: Tue, 26 Aug 2025 10:45:12 +0100 Subject: [PATCH 8/8] Edit to not use assert statement within function. --- improver/utilities/compare.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/improver/utilities/compare.py b/improver/utilities/compare.py index 88868fdfa1..598349c85e 100644 --- a/improver/utilities/compare.py +++ b/improver/utilities/compare.py @@ -408,13 +408,15 @@ def raise_reporter(message): difference_found = False try: - assert output.n_features_in_ == kgo.n_features_in_ - assert output.n_outputs_ == kgo.n_outputs_ - assert output.max_depth == kgo.max_depth - assert output.n_estimators == kgo.n_estimators - assert output.random_state == kgo.random_state + np.testing.assert_equal(output.n_features_in_, kgo.n_features_in_) + np.testing.assert_equal(output.n_outputs_, kgo.n_outputs_) + np.testing.assert_equal(output.max_depth, kgo.max_depth) + np.testing.assert_equal(output.n_estimators, kgo.n_estimators) + np.testing.assert_equal(output.random_state, kgo.random_state) for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_): - assert (output_estimator.tree_.value == kgo_estimator.tree_.value).all() + np.testing.assert_allclose( + output_estimator.tree_.value, kgo_estimator.tree_.value + ) except (AssertionError, ValueError): difference_found = True # call the reporter function outside the except block to avoid nested