Merge b2e3519 into 0b2f47a

environment-dev.yml

@@ -8,6 +8,7 @@ dependencies:
   - numpy
   - pooch >=1.8.0
   - pydantic >=2.0,<2.5.3  # FIXME: Remove pin once our dependencies (xclim, xscen) support pydantic 2.5.3
+  - spotpy
   - statsmodels
   - xarray
   - xclim >=0.47.0  # FIXME: Remove pin once our dependencies (xclim, xscen) support pandas 2.2.0

environment.yml

@@ -8,6 +8,7 @@ dependencies:
   - numpy
   - pooch >=1.8.0
   - pydantic >=2.0,<2.5.3  # FIXME: Remove pin once our dependencies (xclim, xscen) support pydantic 2.5.3
+  - spotpy
   - statsmodels
   - xarray
   - xclim >=0.47.0  # FIXME: Remove pin once our dependencies (xclim, xscen) support pandas 2.2.0

pyproject.toml

@@ -40,6 +40,7 @@ dependencies = [
   "numpy",
   "pooch>=1.8.0",
   "pydantic>=2.0,<2.5.3",
+  "spotpy",
   "statsmodels",
   "xarray",
   "xclim>=0.47.0",

tests/test_calibration.py

@@ -0,0 +1,138 @@
+"""Test suite for the calibration algorithm in calibration.py."""
+
+# Also tests the dummy model implementation.
+import numpy as np
+import pytest
+
+from xhydro.modelling.calibration import perform_calibration
+from xhydro.modelling.hydrological_modelling import _dummy_model
+from xhydro.modelling.obj_funcs import get_objective_function, transform_flows
+
+
+def test_spotpy_calibration():
+    """Make sure the calibration works under possible test cases."""
+    bounds_low = np.array([0, 0, 0])
+    bounds_high = np.array([10, 10, 10])
+
+    model_config = {
+        "precip": np.array([10, 11, 12, 13, 14, 15]),
+        "temperature": np.array([10, 3, -5, 1, 15, 0]),
+        "qobs": np.array([120, 130, 140, 150, 160, 170]),
+        "drainage_area": np.array([10]),
+        "model_name": "Dummy",
+    }
+
+    mask = np.array([0, 0, 0, 0, 1, 1])
+
+    best_parameters, best_simulation, best_objfun = perform_calibration(
+        model_config,
+        "mae",
+        bounds_low=bounds_low,
+        bounds_high=bounds_high,
+        evaluations=1000,
+        algorithm="DDS",
+        mask=mask,
+        sampler_kwargs=dict(trials=1),
+    )
+
+    # Test that the results have the same size as expected (number of parameters)
+    assert len(best_parameters) == len(bounds_high)
+
+    # Test that the objective function is calculated correctly
+    objfun = get_objective_function(
+        model_config["qobs"],
+        best_simulation,
+        obj_func="mae",
+        mask=mask,
+    )
+
+    assert objfun == best_objfun
+
+    # Test dummy model response
+    model_config["parameters"] = [5, 5, 5]
+    qsim = _dummy_model(model_config)
+    assert qsim["qsim"].values[3] == 3500.00
+
+    # Also test to ensure SCEUA and take_minimize is required.
+    best_parameters_sceua, best_simulation, best_objfun = perform_calibration(
+        model_config,
+        "mae",
+        bounds_low=bounds_low,
+        bounds_high=bounds_high,
+        evaluations=10,
+        algorithm="SCEUA",
+    )
+
+    assert len(best_parameters_sceua) == len(bounds_high)
+
+    # Also test to ensure SCEUA and take_minimize is required.
+    best_parameters_negative, best_simulation, best_objfun = perform_calibration(
+        model_config,
+        "nse",
+        bounds_low=bounds_low,
+        bounds_high=bounds_high,
+        evaluations=10,
+        algorithm="SCEUA",
+    )
+    assert len(best_parameters_negative) == len(bounds_high)
+
+    # Test to see if transform works
+    best_parameters_transform, best_simulation, best_objfun = perform_calibration(
+        model_config,
+        "nse",
+        bounds_low=bounds_low,
+        bounds_high=bounds_high,
+        evaluations=10,
+        algorithm="SCEUA",
+        transform="inv",
+        epsilon=0.01,
+    )
+    assert len(best_parameters_transform) == len(bounds_high)
+
+
+def test_calibration_failure_mode_unknown_optimizer():
+    """Test for maximize-minimize failure mode:
+    use "OTHER" optimizer, i.e. an unknown optimizer. Should fail.
+    """
+    bounds_low = np.array([0, 0, 0])
+    bounds_high = np.array([10, 10, 10])
+    model_config = {
+        "precip": np.array([10, 11, 12, 13, 14, 15]),
+        "temperature": np.array([10, 3, -5, 1, 15, 0]),
+        "qobs": np.array([120, 130, 140, 150, 160, 170]),
+        "drainage_area": np.array([10]),
+        "model_name": "Dummy",
+    }
+    with pytest.raises(NotImplementedError) as pytest_wrapped_e:
+        best_parameters_transform, best_simulation, best_objfun = perform_calibration(
+            model_config,
+            "nse",
+            bounds_low=bounds_low,
+            bounds_high=bounds_high,
+            evaluations=10,
+            algorithm="OTHER",
+        )
+        assert pytest_wrapped_e.type == NotImplementedError
+
+
+def test_transform():
+    """Test the flow transformer"""
+    qsim = np.array([10, 10, 10])
+    qobs = np.array([5, 5, 5])
+
+    qsim_r, qobs_r = transform_flows(qsim, qobs, transform="inv", epsilon=0.01)
+    np.testing.assert_array_almost_equal(qsim_r[1], 0.0995024, 6)
+    np.testing.assert_array_almost_equal(qobs_r[1], 0.1980198, 6)
+
+    qsim_r, qobs_r = transform_flows(qsim, qobs, transform="sqrt")
+    np.testing.assert_array_almost_equal(qsim_r[1], 3.1622776, 6)
+    np.testing.assert_array_almost_equal(qobs_r[1], 2.2360679, 6)
+
+    qsim_r, qobs_r = transform_flows(qsim, qobs, transform="log", epsilon=0.01)
+    np.testing.assert_array_almost_equal(qsim_r[1], 2.3075726, 6)
+    np.testing.assert_array_almost_equal(qobs_r[1], 1.6193882, 6)
+
+    # Test Qobs different length than Qsim
+    with pytest.raises(NotImplementedError) as pytest_wrapped_e:
+        qobs_r, qobs_r = transform_flows(qsim, qobs, transform="a", epsilon=0.01)
+        assert pytest_wrapped_e.type == NotImplementedError

tests/test_hydrological_modelling.py

@@ -0,0 +1,53 @@
+"""Test suite for hydrological modelling in hydrological_modelling.py"""
+
+import numpy as np
+import pytest
+
+from xhydro.modelling.hydrological_modelling import (
+    get_hydrological_model_inputs,
+    run_hydrological_model,
+)
+
+
+def test_hydrological_modelling():
+    """Test the hydrological models as they become online"""
+    # Test the dummy model
+    model_config = {
+        "precip": np.array([10, 11, 12, 13, 14, 15]),
+        "temperature": np.array([10, 3, -5, 1, 15, 0]),
+        "qobs": np.array([120, 130, 140, 150, 160, 170]),
+        "drainage_area": np.array([10]),
+        "model_name": "Dummy",
+        "parameters": np.array([5, 5, 5]),
+    }
+    qsim = run_hydrological_model(model_config)
+    assert qsim["qsim"].values[3] == 3500.00
+
+    # Test the exceptions for new models
+    model_config.update(model_name="ADD_OTHER_HERE")
+    qsim = run_hydrological_model(model_config)
+    assert qsim == 0
+
+
+def test_import_unknown_model():
+    """Test for unknown model"""
+    with pytest.raises(NotImplementedError) as pytest_wrapped_e:
+        model_config = {"model_name": "fake_model"}
+        _ = run_hydrological_model(model_config)
+        assert pytest_wrapped_e.type == NotImplementedError
+
+
+def test_get_unknown_model_requirements():
+    """Test for required inputs for models with unknown name"""
+    with pytest.raises(NotImplementedError) as pytest_wrapped_e:
+        model_name = "fake_model"
+        _ = get_hydrological_model_inputs(model_name)
+        assert pytest_wrapped_e.type == NotImplementedError
+
+
+def test_get_model_requirements():
+    """Test for required inputs for models"""
+    model_name = "Dummy"
+    required_config = get_hydrological_model_inputs(model_name)
+    print(required_config.keys())
+    assert len(required_config.keys()) == 4

tests/test_objective_functions.py

@@ -0,0 +1,160 @@
+"""Test suite for the objective functions in obj_funcs.py."""
+
+import numpy as np
+import pytest
+
+from xhydro.modelling.obj_funcs import (
+    _get_objfun_minimize_or_maximize,
+    get_objective_function,
+)
+
+
+def test_obj_funcs():
+    """Series of tests to test all objective functions with fast test data"""
+    qobs = np.array([120, 130, 140, 150, 160, 170])
+    qsim = np.array([120, 125, 145, 140, 140, 180])
+
+    # Test that the objective function is calculated correctly
+    objfun = get_objective_function(qobs, qsim, obj_func="abs_bias")
+    np.testing.assert_array_almost_equal(objfun, 3.3333333333333335, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="abs_pbias")
+    np.testing.assert_array_almost_equal(objfun, 2.2988505747126435, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="abs_volume_error")
+    np.testing.assert_array_almost_equal(objfun, 0.022988505747126436, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="agreement_index")
+    np.testing.assert_array_almost_equal(objfun, 0.9171974522292994, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="bias")
+    np.testing.assert_array_almost_equal(objfun, -3.3333333333333335, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="correlation_coeff")
+    np.testing.assert_array_almost_equal(objfun, 0.8599102447336393, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="kge")
+    np.testing.assert_array_almost_equal(objfun, 0.8077187696552522, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="kge_mod")
+    np.testing.assert_array_almost_equal(objfun, 0.7888769531580001, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="mae")
+    np.testing.assert_array_almost_equal(objfun, 8.333333333333334, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="mare")
+    np.testing.assert_array_almost_equal(objfun, 0.05747126436781609, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="mse")
+    np.testing.assert_array_almost_equal(objfun, 108.33333333333333, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="nse")
+    np.testing.assert_array_almost_equal(objfun, 0.6285714285714286, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="pbias")
+    np.testing.assert_array_almost_equal(objfun, -2.2988505747126435, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="r2")
+    np.testing.assert_array_almost_equal(objfun, 0.7394456289978675, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="rmse")
+    np.testing.assert_array_almost_equal(objfun, 10.408329997330663, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="rrmse")
+    np.testing.assert_array_almost_equal(objfun, 0.07178158618848733, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="rsr")
+    np.testing.assert_array_almost_equal(objfun, 0.6094494002200439, 8)
+
+    objfun = get_objective_function(qobs, qsim, obj_func="volume_error")
+    np.testing.assert_array_almost_equal(objfun, -0.022988505747126436, 8)
+
+
+def test_objective_function_failure_data_length():
+    """Test for the objective function calculation failure mode:
+    qobs and qsim length are different
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = get_objective_function(
+            np.array([100, 110]),
+            np.array([100, 110, 120]),
+            obj_func="mae",
+        )
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_objective_function_failure_mask_length():
+    """Test for the objective function calculation failure mode:
+    qobs and mask length are different
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = get_objective_function(
+            np.array([100, 100, 100]),
+            np.array([100, 110, 120]),
+            obj_func="mae",
+            mask=np.array([0, 1, 0, 0]),
+        )
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_objective_function_failure_unknown_objfun():
+    """Test for the objective function calculation failure mode:
+    Objective function is unknown
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = get_objective_function(
+            np.array([100, 100, 100]),
+            np.array([100, 110, 120]),
+            obj_func="fake",
+        )
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_objective_function_failure_mask_contents():
+    """Test for the objective function calculation failure mode:
+    Mask contains other than 0 and 1
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = get_objective_function(
+            np.array([100, 100, 100]),
+            np.array([100, 110, 120]),
+            obj_func="mae",
+            mask=np.array([0, 0.5, 1]),
+        )
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_maximizer_objfun_failure_modes_bias():
+    """Test for maximize-minimize failure mode:
+    Use of bias objfun which is unbounded
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = _get_objfun_minimize_or_maximize(obj_func="bias")
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_maximizer_objfun_failure_modes_pbias():
+    """Test for maximize-minimize failure mode:
+    Use of pbias objfun which is unbounded
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = _get_objfun_minimize_or_maximize(obj_func="pbias")
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_maximizer_objfun_failure_modes_volume_error():
+    """Test for maximize-minimize failure mode:
+    Use of volume_error objfun which is unbounded
+    """
+    with pytest.raises(ValueError) as pytest_wrapped_e:
+        _ = _get_objfun_minimize_or_maximize(obj_func="volume_error")
+        assert pytest_wrapped_e.type == ValueError
+
+
+def test_maximizer_objfun_failure_modes_unknown_metric():
+    """Test for maximize-minimize failure mode:
+    Use of unknown objfun
+    """
+    with pytest.raises(NotImplementedError) as pytest_wrapped_e:
+        _ = _get_objfun_minimize_or_maximize(obj_func="unknown_of")
+        assert pytest_wrapped_e.type == NotImplementedError