From 0fd102d0d8dc32e28cf3919347b09e38399a719c Mon Sep 17 00:00:00 2001 From: Phillip Schuster Date: Tue, 8 Sep 2026 14:19:00 +0200 Subject: [PATCH] fix: handle model edge cases --- docs/testing.md | 8 ++-- matilda/core.py | 8 +++- matilda/mspot_glacier.py | 9 ++-- tests/synthetic.py | 24 ++++++++++ tests/test_model_modes.py | 98 +++++++++++++++++++++++++++++++++++++++ tests/test_mspot.py | 48 ++++++++++++++++++- 6 files changed, 186 insertions(+), 9 deletions(-) diff --git a/docs/testing.md b/docs/testing.md index 6d0882a..5c877d2 100644 --- a/docs/testing.md +++ b/docs/testing.md @@ -67,7 +67,8 @@ The synthetic tests cover: - zero-glacier and fixed-glacier public simulation paths; - non-negative stores and fluxes where physically required; - calendar completeness, deterministic repetition, and input immutability; -- one model evaluation initiated through `mspot`. +- glacier loss during the first annual geometry update; +- model evaluations initiated through `mspot`, including glacier-only sampling. The runoff series used by the `mspot` test is synthetic and non-constant. Its only purpose is to exercise unit conversion, alignment, and objective-function @@ -91,9 +92,8 @@ change has been assessed and accepted. The explicit command for that step is: .venv/bin/python -m tests.generate_synthetic_references --replace ``` -Separate tests are still required for complete glacier loss, positive -cumulative mass-balance handling, and serial-versus-parallel calibration -equivalence. +Separate tests are still required for positive cumulative mass-balance +handling and serial-versus-parallel calibration equivalence. The continuous-integration workflow also builds a wheel and verifies that its package files are byte-for-byte copies of the maintained sources. diff --git a/matilda/core.py b/matilda/core.py index bfe1e19..cc69334 100644 --- a/matilda/core.py +++ b/matilda/core.py @@ -1107,6 +1107,7 @@ def updated_glacier_melt( # Setup initial variables for main loop new_area = parameter.area_glac + new_distribution = init_elev smb_cum = 0 surplus = 0 warn = True @@ -1171,6 +1172,7 @@ def updated_glacier_melt( for i in range(len(data_update.water_year.unique())): year = data_update.water_year.unique()[i] mask = data_update.water_year == year + smb_flag = False # Use updated glacier area of the previous year parameter_updated.area_glac = new_area @@ -1960,6 +1962,10 @@ def matilda_submodules( else: lookup_table = str("No lookup table generated") glacier_change = str("No glacier changes calculated") + if parameter.ele_dat is not None: + _, input_df_catchment = input_scaling(df_preproc, parameter) + else: + input_df_catchment = df_preproc.copy() else: print( @@ -2002,7 +2008,7 @@ def matilda_submodules( glacier_change = str("No glacier changes calculated") # Execute HBV module: - if glacier_profile is not None: + if glacier_profile is not None and parameter.area_glac > 0: output_HBV = hbv_simulation( input_df_catchment, parameter, glacier_area=glacier_change ) diff --git a/matilda/mspot_glacier.py b/matilda/mspot_glacier.py index f6849ff..bd392a2 100644 --- a/matilda/mspot_glacier.py +++ b/matilda/mspot_glacier.py @@ -1630,9 +1630,12 @@ def psample( **kwargs, ) - psample_setup = setup( - df, obs, target_swe, obj_func - ) # Define custom objective function using obj_func= + if glacier_only: + psample_setup = setup(df, obs, obj_func) + else: + psample_setup = setup( + df, obs, target_swe, obj_func + ) # Define custom objective function using obj_func= alg_selector = { "mc": spotpy.algorithms.mc, "sceua": spotpy.algorithms.sceua, diff --git a/tests/synthetic.py b/tests/synthetic.py index cf7b3a3..e2a97ad 100644 --- a/tests/synthetic.py +++ b/tests/synthetic.py @@ -106,6 +106,30 @@ def make_synthetic_observations() -> pd.DataFrame: return pd.DataFrame({"Date": dates, "Qobs": runoff_m3_per_second}) +def make_synthetic_glacier_profile() -> pd.DataFrame: + """Return a small glacier profile covering 20 percent of the catchment.""" + return pd.DataFrame( + { + "Elevation": [1400.0, 1600.0, 1800.0], + "Area": [0.05, 0.10, 0.05], + "WE": [4000.0, 6000.0, 5000.0], + "EleZone": [1400, 1600, 1800], + } + ) + + +def make_synthetic_mass_balance_observations() -> pd.DataFrame: + """Return annual mass-balance records for a glacier-only smoke test.""" + return pd.DataFrame( + { + "YEAR": ["2000-01-01", "2001-01-01"], + "BEGIN_PERIOD": ["2000-01-01", "2001-01-01"], + "END_PERIOD": ["2000-12-31", "2001-12-31"], + "ANNUAL_BALANCE": [-300.0, -350.0], + } + ) + + def load_synthetic_reference(name: str): """Load one maintained synthetic public-API output.""" manifest = json.loads(REFERENCE_MANIFEST_PATH.read_text(encoding="utf-8")) diff --git a/tests/test_model_modes.py b/tests/test_model_modes.py index 1a94045..a0e4d02 100644 --- a/tests/test_model_modes.py +++ b/tests/test_model_modes.py @@ -29,6 +29,7 @@ SYNTHETIC_FRAME_OUTPUTS, load_synthetic_reference, make_synthetic_forcing, + make_synthetic_glacier_profile, model_settings, ) @@ -244,3 +245,100 @@ def test_fixed_glacier_mode_is_exactly_repeatable(fixed_glacier_run): ) assert first_output[2] == second_output[2] assert first_output[4:] == second_output[4:] + + +def test_zero_glacier_evolution_matches_standard_zero_glacier_mode( + zero_glacier_run, +): + forcing = make_synthetic_forcing() + forcing_before = forcing.copy(deep=True) + profile = make_synthetic_glacier_profile() + profile_before = profile.copy(deep=True) + settings = model_settings(area_glac=0.0) + settings["elev_rescaling"] = True + + try: + with redirect_stdout(io.StringIO()): + evolving_output = matilda_simulation( + forcing, + glacier_profile=profile, + **settings, + ) + finally: + plt.close("all") + + standard_output = zero_glacier_run[0] + assert_frame_equal(forcing, forcing_before, check_exact=True) + assert_frame_equal(profile, profile_before, check_exact=True) + for position in (0, 1, 3): + assert_frame_equal( + evolving_output[position], + standard_output[position], + check_exact=True, + ) + assert evolving_output[2] == standard_output[2] + assert evolving_output[4:] == standard_output[4:] + + +def test_one_water_year_glacier_evolution_completes(): + forcing = make_synthetic_forcing() + forcing_before = forcing.copy(deep=True) + profile = make_synthetic_glacier_profile() + profile_before = profile.copy(deep=True) + settings = model_settings(area_glac=20.0) + settings.update( + { + "sim_end": "2000-09-30", + "elev_rescaling": True, + } + ) + + try: + with redirect_stdout(io.StringIO()): + output = matilda_simulation( + forcing, + glacier_profile=profile, + **settings, + ) + finally: + plt.close("all") + + assert_frame_equal(forcing, forcing_before, check_exact=True) + assert_frame_equal( + profile.loc[:, profile_before.columns], + profile_before, + check_exact=True, + ) + assert len(output[1]) == 274 + assert np.isfinite(output[1].select_dtypes(include=np.number)).all().all() + assert output[5]["glacier_area"].tolist() == [20.0] + + +def test_glacier_evolution_handles_loss_in_first_update(): + forcing = make_synthetic_forcing() + forcing_before = forcing.copy(deep=True) + profile = make_synthetic_glacier_profile() + profile["WE"] = 1.0 + profile_before = profile.copy(deep=True) + settings = model_settings(area_glac=20.0) + settings["elev_rescaling"] = True + + try: + with redirect_stdout(io.StringIO()): + output = matilda_simulation( + forcing, + glacier_profile=profile, + **settings, + ) + finally: + plt.close("all") + + assert_frame_equal(forcing, forcing_before, check_exact=True) + assert_frame_equal( + profile.loc[:, profile_before.columns], + profile_before, + check_exact=True, + ) + assert output[5]["glacier_area"].iloc[1] == 0 + assert np.isfinite(output[5]["glacier_elev"]).all() + assert np.isfinite(output[1].select_dtypes(include=np.number)).all().all() diff --git a/tests/test_mspot.py b/tests/test_mspot.py index 466789f..2f7526e 100644 --- a/tests/test_mspot.py +++ b/tests/test_mspot.py @@ -2,11 +2,13 @@ from __future__ import annotations +import random + import numpy as np from pandas.testing import assert_frame_equal, assert_index_equal import pytest -from matilda.mspot_glacier import spot_setup +from matilda.mspot_glacier import psample, spot_setup from tests.synthetic import ( AREA_CATCHMENT, SETUP_END, @@ -15,6 +17,8 @@ SIMULATION_START, calibration_parameters, make_synthetic_forcing, + make_synthetic_glacier_profile, + make_synthetic_mass_balance_observations, make_synthetic_observations, ) @@ -67,3 +71,45 @@ def test_mspot_model_evaluation_is_finite_and_has_no_file_side_effects( assert_frame_equal(forcing, forcing_before, check_exact=True) assert_frame_equal(observations, observations_before, check_exact=True) assert list(tmp_path.iterdir()) == [] + + +def test_glacier_only_psample_completes_one_evaluation(tmp_path): + numpy_random_state = np.random.get_state() + python_random_state = random.getstate() + forcing = make_synthetic_forcing() + observations = make_synthetic_mass_balance_observations() + profile = make_synthetic_glacier_profile() + + try: + np.random.seed(0) + results = psample( + forcing, + observations, + rep=1, + output=tmp_path, + dbname="glacier_only_smoke", + dbformat="ram", + set_up_start=SETUP_START, + set_up_end=SETUP_END, + sim_start=SIMULATION_START, + sim_end=SIMULATION_END, + freq="D", + lat=45.0, + area_cat=AREA_CATCHMENT, + area_glac=20.0, + ele_dat=1000.0, + ele_glac=1600.0, + glacier_profile=profile, + glacier_only=True, + obs_type="annual", + algorithm="lhs", + obj_dir="minimize", + save_sim=False, + ) + finally: + np.random.set_state(numpy_random_state) + random.setstate(python_random_state) + + assert results["best_index"] == 0 + assert np.isfinite(results["best_objf"]) + assert (tmp_path / "glacier_only_smoke_observations.csv").is_file()