Handle model edge cases

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.

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
         )

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,

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"))

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()

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()