productivity indicators

Author: michielkallenberg

cybench/config.py

@@ -150,11 +150,14 @@
 # Soil moisture indicators: surface moisture, root zone moisture
 SOIL_MOISTURE_INDICATORS = ["ssm"]  # , "rsm"]
 
+PRODUCTIVITY_INDICATORS = ["twso"]
+
 TIME_SERIES_INPUTS = {
-    "meteo": METEO_INDICATORS,
-    "fpar": [RS_FPAR],
-    "ndvi": [RS_NDVI],
-    "soil_moisture": SOIL_MOISTURE_INDICATORS,
+    # "meteo": METEO_INDICATORS,
+    # "fpar": [RS_FPAR],
+    # "ndvi": [RS_NDVI],
+    # "soil_moisture": SOIL_MOISTURE_INDICATORS,
+    "twso": PRODUCTIVITY_INDICATORS,
 }
 
 # Time series predictors
@@ -171,6 +174,7 @@
     RS_FPAR: "mean",
     RS_NDVI: "mean",
     "ssm": "mean",
+    "twso": "max",
 }
 
 # All predictors. Add more when available
@@ -198,10 +202,11 @@
 # Lead time for forecasting
 # Choices: "middle-of-season", "quarter-of-season",
 # "n-day(s)" where n is an integer
-FORECAST_LEAD_TIME = "middle-of-season"
+# FORECAST_LEAD_TIME = "middle-of-season"
+FORECAST_LEAD_TIME = "0-days"
 
 # Buffer period before the start of season
-SPINUP_DAYS = 90
+SPINUP_DAYS = 0
 
 # Logging
 PATH_LOGS_DIR = os.path.join(PATH_OUTPUT_DIR, "logs")

cybench/datasets/alignment.py

@@ -271,7 +271,7 @@ def align_to_crop_season_window_numpy(
     valid_start = np.zeros_like(valid_mask, dtype=bool)
     valid_end = np.zeros_like(valid_mask, dtype=bool)
 
-    tolerance = np.timedelta64(10, "D")  # indicators may not come in daily timesteps
+    tolerance = np.timedelta64(15, "D")  # indicators may not come in daily timesteps
     valid_start[valid_mask] = date_min[valid_mask] - tolerance < (
         cutoff_dates[crop_indices[valid_mask]]
         - season_window_lengths[crop_indices[valid_mask]].astype("timedelta64[D]")

cybench/models/naive_models.py

@@ -100,3 +100,65 @@ def load(cls, model_name):
             saved_model = pickle.load(f)
 
         return saved_model
+
+
+class MaxPredictorModel(BaseModel):
+    """A naive model that predicts the max value of a predefined feature column
+    for each test item. No training; prediction uses only test data.
+    """
+
+    def __init__(self, feature_column):
+        self.feature_column = feature_column
+        self._logger = logging.getLogger(__name__)
+
+    def fit(self, dataset, **fit_params):
+        """No training needed; return self."""
+        self._logger.info("MaxFromItemModel: nothing to fit.")
+        return self, {}
+
+    def predict_items(self, X: list):
+        """Predict the max value contained in each test item.
+
+        Assumes each item is a dict where item[self.feature_column]
+        is a list, array, or iterable of numeric values.
+        """
+        predictions = np.zeros(len(X))
+        for i, item in enumerate(X):
+            values = item.get(self.feature_column, None)
+
+            if values is None:
+                raise ValueError(
+                    f"Item {i} does not contain the feature column '{self.feature_column}'."
+                )
+
+            try:
+                max_val = (
+                    0.001 * np.nanmax(values) if np.any(~np.isnan(values)) else np.NAN
+                )
+                predictions[i] = max_val
+            except Exception as e:
+                raise ValueError(f"Could not compute max for item {i}: {e}")
+        return predictions, {}
+
+    def save(self, model_name):
+        """Save model, e.g. using pickle.
+
+        Args:
+          model_name: Filename that will be used to save the model.
+        """
+        with open(model_name, "wb") as f:
+            pickle.dump(self, f)
+
+    def load(cls, model_name):
+        """Deserialize a saved model.
+
+        Args:
+          model_name: Filename that was used to save the model.
+
+        Returns:
+          The deserialized model.
+        """
+        with open(model_name, "rb") as f:
+            saved_model = pickle.load(f)
+
+        return saved_model

cybench/runs/run_benchmark.py

@@ -21,7 +21,7 @@
     get_default_metrics,
     prepare_targets_preds,
 )
-from cybench.models.naive_models import AverageYieldModel
+from cybench.models.naive_models import AverageYieldModel, MaxPredictorModel
 from cybench.models.trend_models import TrendModel
 from cybench.models.sklearn_models import SklearnRidge, SklearnRandomForest
 from cybench.models.xgboost_model import XGBoostModel
@@ -43,6 +43,7 @@
 
 _BASELINE_MODEL_CONSTRUCTORS = {
     "AverageYieldModel": AverageYieldModel,
+    "MaxPredictorModel": MaxPredictorModel,
     "LinearTrend": TrendModel,
     "SklearnRidge": SklearnRidge,
     "RidgeRes": RidgeRes,
@@ -60,6 +61,8 @@
 
 BASELINE_MODELS = list(_BASELINE_MODEL_CONSTRUCTORS.keys())
 
+BASELINE_MODELS = ["AverageYieldModel", "MaxPredictorModel"]
+
 _BASELINE_MODEL_INIT_KWARGS = defaultdict(dict)
 
 NN_MODELS_EPOCHS = 50
@@ -90,6 +93,45 @@
     "epochs": NN_MODELS_EPOCHS,
     "device": "cuda" if torch.cuda.is_available() else "cpu",
 }
+_BASELINE_MODEL_INIT_KWARGS["MaxPredictorModel"] = {
+    "feature_column": "twso",
+}
+
+_BASELINE_MODEL_FIT_KWARGS["MaxPredictorModel"] = {
+    "feature_column": "twso",
+}
+
+
+def discover_datasets_from_disk(path_data_dir: str):
+    """
+    Discover datasets by checking for:
+      PATH_DATA_DIR/<crop>/<region>/twso_<crop>_<region>.csv
+    Returns list of dataset_name strings like "maize_FR".
+    """
+    found = []
+
+    # iterate crops as folders under PATH_DATA_DIR
+    if not os.path.isdir(path_data_dir):
+        raise FileNotFoundError(
+            f"PATH_DATA_DIR does not exist or is not a directory: {path_data_dir}"
+        )
+
+    for crop in sorted(os.listdir(path_data_dir)):
+        crop_dir = os.path.join(path_data_dir, crop)
+        if not os.path.isdir(crop_dir):
+            continue
+
+        # iterate regions as folders under crop
+        for region in sorted(os.listdir(crop_dir)):
+            region_dir = os.path.join(crop_dir, region)
+            if not os.path.isdir(region_dir):
+                continue
+
+            twso_file = os.path.join(region_dir, f"twso_{crop}_{region}.csv")
+            if os.path.exists(twso_file):
+                found.append(f"{crop}_{region}")
+
+    return found
 
 
 def run_benchmark(
@@ -164,6 +206,7 @@ def run_benchmark(
     else:
         sel_years = all_years
 
+    all_results = []
     for test_year in sel_years:
         train_years = [y for y in all_years if y != test_year]
         test_years = [test_year]
@@ -194,13 +237,14 @@ def run_benchmark(
         df = pd.DataFrame.from_dict(model_output)
         df[KEY_COUNTRY] = df[KEY_LOC].str[:2]
         df.set_index([KEY_COUNTRY, KEY_LOC, KEY_YEAR], inplace=True)
-        df.to_csv(os.path.join(path_results, f"{dataset_name}_year_{test_year}.csv"))
+        all_results.append(df)
 
-    df_metrics = compute_metrics(run_name, list(model_constructors.keys()))
+    df_all = pd.concat(all_results).sort_index()
 
-    return {
-        "df_metrics": df_metrics,
-    }
+    results_file = os.path.join(path_results, f"{dataset_name}.csv")
+    print(f"write results to {results_file}")
+    df_all.to_csv(results_file)
+    return
 
 
 def load_results(
@@ -228,6 +272,8 @@ def load_results(
 
     df_all = pd.DataFrame()
     for file in files:
+        if not file.lower().endswith(".csv"):
+            continue
         path = os.path.join(path_results, file)
         df = pd.read_csv(path)
         df_all = pd.concat([df_all, df], axis=0)
@@ -339,7 +385,12 @@ def run_benchmark_on_all_data():
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(prog="run_benchmark.py", description="Run cybench")
     parser.add_argument("-r", "--run-name")
-    parser.add_argument("-d", "--dataset-name")
+    parser.add_argument(
+        "-d",
+        "--dataset-name",
+        default=None,
+        help="Dataset name (e.g. maize_FR). If omitted or 'all'/'none', run all datasets.",
+    )
     parser.add_argument("-m", "--mode")
     parser.add_argument(
         "-y",
@@ -350,48 +401,45 @@ def run_benchmark_on_all_data():
         help="Test year(s)",
     )
     args = parser.parse_args()
-    dataset_name = args.dataset_name
-    assert dataset_name is not None
-
-    if args.run_name is not None:
-        run_name = args.run_name
-    else:
-        run_name = dataset_name
 
-    if args.years is None or [y.lower() for y in args.years] in (["none"], ["all"]):
-        args.years = None
+    if args.dataset_name is None or str(args.dataset_name).lower() in ("none", "all"):
+        dataset_names = discover_datasets_from_disk(PATH_DATA_DIR)
+        if not dataset_names:
+            raise FileNotFoundError(
+                f"No datasets found. Expected files like "
+                f"{PATH_DATA_DIR}/<crop>/<region>/twso_<crop>_<region>.csv"
+            )
     else:
-        args.years = [int(y) for y in args.years]
-    sel_years = args.years
-
-    if (args.mode is not None) and args.mode == "test":
-        # skipping some models
-        baseline_models = [
-            "AverageYieldModel",
-            "LinearTrend",
-            "SklearnRidge",
-            "RidgeRes",
-            "LSTM",
-            "LSTMRes",
-        ]
-        # override epochs for nn-models
-        nn_models_epochs = 5
-        results = run_benchmark(
-            run_name=run_name,
-            dataset_name=dataset_name,
-            baseline_models=baseline_models,
-            nn_models_epochs=nn_models_epochs,
-        )
-    else:
-        results = run_benchmark(
-            run_name=run_name, dataset_name=dataset_name, sel_years=sel_years
-        )
-
-    index_cols = results["df_metrics"].index.names
-    df_metrics = results["df_metrics"].reset_index()
-
-    metric_cols = [c for c in df_metrics.columns if c not in index_cols]
-
-    # Group and average all available metrics
-    agg_df = df_metrics.groupby("model")[metric_cols].mean().round(3)
-    print(agg_df)
+        dataset_names = [args.dataset_name]
+
+    for dataset_name in dataset_names:
+        run_name = args.run_name if args.run_name is not None else dataset_name
+
+        if args.years is None or [y.lower() for y in args.years] in (["none"], ["all"]):
+            args.years = None
+        else:
+            args.years = [int(y) for y in args.years]
+        sel_years = args.years
+
+        if (args.mode is not None) and args.mode == "test":
+            # skipping some models
+            baseline_models = [
+                "AverageYieldModel",
+                "LinearTrend",
+                "SklearnRidge",
+                "RidgeRes",
+                "LSTM",
+                "LSTMRes",
+            ]
+            # override epochs for nn-models
+            nn_models_epochs = 5
+            results = run_benchmark(
+                run_name=run_name,
+                dataset_name=dataset_name,
+                baseline_models=baseline_models,
+                nn_models_epochs=nn_models_epochs,
+            )
+        else:
+            run_benchmark(
+                run_name=run_name, dataset_name=dataset_name, sel_years=sel_years
+            )