ENH: Mask active cell when checking if parameter is inside feasible d…

…omain

- Allow to pass any value to a non active cell and it does not return an error.
- Ruff and fprettify formatting applied
- Remove unused constant TOL_BOUNDS since F_PRECISION is now used

smash/_constant.py

@@ -34,6 +34,7 @@ def get_rr_internal_fluxes_from_structure(structure: str) -> list[str]:
     [rr_internal_fluxes.extend(MODULE_RR_INTERNAL_FLUXES[module]) for module in structure.split("-")]
     return rr_internal_fluxes
 
+
 ### FLOAT PRECISION FOR FLOAT COMPARISON ###
 F_PRECISION = 1.0e-5
 
@@ -400,7 +401,7 @@ def get_rr_internal_fluxes_from_structure(structure: str) -> list[str]:
             (1e-6, 0.999999),  # % hi
             (1e-6, 0.999999),  # % hp
             (1e-6, 0.999999),  # % ht
-            (-1e3, 0),         # % he
+            (-1e3, 0),  # % he
             (1e-6, 0.999999),  # % ha
             (1e-6, 0.999999),  # % hc
             (1e-6, 0.999999),  # % hcl
@@ -412,8 +413,6 @@ def get_rr_internal_fluxes_from_structure(structure: str) -> list[str]:
     )
 )
 
-TOL_BOUNDS = 1e-9
-
 ### OPTIMIZABLE PARAMETERS ###
 ##############################
 

smash/core/model/_standardize.py

@@ -533,17 +533,26 @@ def _standardize_rr_parameters_value(
     if not isinstance(value, (int, float, np.ndarray)):
         raise TypeError("value argument must be of Numeric type (int, float) or np.ndarray")
 
-    arr = np.array(value, ndmin=1)
-    low, upp = FEASIBLE_RR_PARAMETERS[key]
-    low_arr = np.min(arr)
-    upp_arr = np.max(arr)
+    arr = np.array(value, ndmin=1, dtype=np.float32)
 
-    if isinstance(value, np.ndarray) and value.shape != model.mesh.flwdir.shape and value.size != 1:
+    if arr.shape != model.mesh.flwdir.shape and arr.size != 1:
         raise ValueError(
             f"Invalid shape for model rr_parameter '{key}'. Could not broadcast input array from shape "
-            f"{value.shape} into shape {model.mesh.flwdir.shape}"
+            f"{arr.shape} into shape {model.mesh.flwdir.shape}"
         )
 
+    low, upp = FEASIBLE_RR_PARAMETERS[key]
+
+    if arr.shape == model.mesh.flwdir.shape:
+        # Do not check if a value is inside the feasible domain outside of active cells
+        mask = model.mesh.active_cell == 1
+    else:
+        # Check all values (size == 1)
+        mask = np.ones(arr.shape, dtype=bool)
+
+    low_arr = np.min(arr, where=mask, initial=np.inf)
+    upp_arr = np.max(arr, where=mask, initial=-np.inf)
+
     if (low_arr + F_PRECISION) <= low or (upp_arr - F_PRECISION) >= upp:
         raise ValueError(
             f"Invalid value for model rr_parameter '{key}'. rr_parameter domain [{low_arr}, {upp_arr}] is "
@@ -559,17 +568,25 @@ def _standardize_rr_states_value(
     if not isinstance(value, (int, float, np.ndarray)):
         raise TypeError("value argument must be of Numeric type (int, float) or np.ndarray")
 
-    arr = np.array(value, ndmin=1)
-    low, upp = FEASIBLE_RR_INITIAL_STATES[key]
-    low_arr = np.min(arr)
-    upp_arr = np.max(arr)
+    arr = np.array(value, ndmin=1, dtype=np.float32)
 
-    if isinstance(value, np.ndarray) and value.shape != model.mesh.flwdir.shape and value.size != 1:
+    if arr.shape != model.mesh.flwdir.shape and arr.size != 1:
         raise ValueError(
             f"Invalid shape for model {state_kind} '{key}'. Could not broadcast input array from shape "
-            f"{value.shape} into shape {model.mesh.flwdir.shape}"
+            f"{arr.shape} into shape {model.mesh.flwdir.shape}"
         )
 
+    low, upp = FEASIBLE_RR_INITIAL_STATES[key]
+    if arr.shape == model.mesh.flwdir.shape:
+        # Do not check if a value is inside the feasible domain outside of active cells
+        mask = model.mesh.active_cell == 1
+    else:
+        # Check all values (size == 1)
+        mask = np.ones(arr.shape, dtype=bool)
+
+    low_arr = np.min(arr, where=mask, initial=np.inf)
+    upp_arr = np.max(arr, where=mask, initial=-np.inf)
+
     if (low_arr + F_PRECISION) <= low or (upp_arr - F_PRECISION) >= upp:
         raise ValueError(
             f"Invalid value for model {state_kind} '{key}'. {state_kind} domain [{low_arr}, {upp_arr}] is "
@@ -585,17 +602,18 @@ def _standardize_serr_mu_parameters_value(
     if not isinstance(value, (int, float, np.ndarray)):
         raise TypeError("value argument must be of Numeric type (int, float) or np.ndarray")
 
-    arr = np.array(value, ndmin=1)
-    low, upp = FEASIBLE_SERR_MU_PARAMETERS[key]
-    low_arr = np.min(arr)
-    upp_arr = np.max(arr)
+    arr = np.array(value, ndmin=1, dtype=np.float32)
 
-    if isinstance(value, np.ndarray) and value.shape != (model.mesh.ng,) and value.size != 1:
+    if arr.shape != (model.mesh.ng,) and arr.size != 1:
         raise ValueError(
             f"Invalid shape for model serr_mu_parameter '{key}'. Could not broadcast input array from shape "
-            f"{value.shape} into shape {(model.mesh.ng,)}"
+            f"{arr.shape} into shape {(model.mesh.ng,)}"
         )
 
+    low, upp = FEASIBLE_SERR_MU_PARAMETERS[key]
+    low_arr = np.min(arr)
+    upp_arr = np.max(arr)
+
     if low_arr <= low or upp_arr >= upp:
         raise ValueError(
             f"Invalid value for model serr_mu_parameter '{key}'. serr_mu_parameter domain "
@@ -611,17 +629,18 @@ def _standardize_serr_sigma_parameters_value(
     if not isinstance(value, (int, float, np.ndarray)):
         raise TypeError("value argument must be of Numeric type (int, float) or np.ndarray")
 
-    arr = np.array(value, ndmin=1)
-    low, upp = FEASIBLE_SERR_SIGMA_PARAMETERS[key]
-    low_arr = np.min(arr)
-    upp_arr = np.max(arr)
+    arr = np.array(value, ndmin=1, dtype=np.float32)
 
-    if isinstance(value, np.ndarray) and value.shape != (model.mesh.ng,) and value.size != 1:
+    if arr.shape != (model.mesh.ng,) and arr.size != 1:
         raise ValueError(
             f"Invalid shape for model serr_sigma_parameter '{key}'. Could not broadcast input array from "
-            f"shape {value.shape} into shape {(model.mesh.ng,)}"
+            f"shape {arr.shape} into shape {(model.mesh.ng,)}"
         )
 
+    low, upp = FEASIBLE_SERR_SIGMA_PARAMETERS[key]
+    low_arr = np.min(arr)
+    upp_arr = np.max(arr)
+
     if low_arr <= low or upp_arr >= upp:
         raise ValueError(
             f"Invalid value for model serr_sigma_parameter '{key}'. serr_sigma_parameter domain "

smash/core/simulation/_doc.py

@@ -1655,8 +1655,7 @@ def _gen_docstring_from_base_doc(
 
 _forward_run_doc_appender = DocAppender(_forward_run_doc, indents=0)
 _smash_forward_run_doc_substitution = DocSubstitution(
-    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model "
-    "`smash`.",
+    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model `smash`.",
     model_return="model : `Model`\n\t It returns an updated copy of the initial Model object.",
     model_example_func="model_fwd = smash.forward_run()",
     model_example_response="model_fwd",
@@ -1672,8 +1671,7 @@ def _gen_docstring_from_base_doc(
 
 _optimize_doc_appender = DocAppender(_optimize_doc, indents=0)
 _smash_optimize_doc_substitution = DocSubstitution(
-    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model "
-    "`smash`.",
+    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model `smash`.",
     default_optimize_options_func="default_optimize_options",
     parameters_serr_mu_parameters="",
     parameters_serr_sigma_parameters="",
@@ -1699,8 +1697,7 @@ def _gen_docstring_from_base_doc(
 
 _multiset_estimate_doc_appender = DocAppender(_multiset_estimate_doc, indents=0)
 _smash_multiset_estimate_doc_substitution = DocSubstitution(
-    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model "
-    "`smash`.",
+    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model `smash`.",
     model_return="model : `Model`\n\t It returns an updated copy of the initial Model object.",
     model_example_func="model_estim = smash.multiset_estimate(model, multiset=mfr)",
     model_example_response="model_estim",
@@ -1716,8 +1713,7 @@ def _gen_docstring_from_base_doc(
 
 _bayesian_optimize_doc_appender = DocAppender(_bayesian_optimize_doc, indents=0)
 _smash_bayesian_optimize_doc_substitution = DocSubstitution(
-    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model "
-    "`smash`.",
+    model_parameter="model : `Model`\n\tPrimary data structure of the hydrological model `smash`.",
     default_optimize_options_func="default_bayesian_optimize_options",
     parameters_serr_mu_parameters="- `Model.serr_mu_parameters`",
     parameters_serr_sigma_parameters="- `Model.serr_sigma_parameters`",

smash/core/simulation/_standardize.py

@@ -225,21 +225,29 @@ def _standardize_simulation_optimize_options_bounds(
         if key in model.rr_parameters.keys:
             arr = model.get_rr_parameters(key)
             low, upp = FEASIBLE_RR_PARAMETERS[key]
+            # Do not check if a value is inside the feasible domain outside of active cells
+            mask = model.mesh.active_cell == 1
         elif key in model.rr_initial_states.keys:
             arr = model.get_rr_initial_states(key)
             low, upp = FEASIBLE_RR_INITIAL_STATES[key]
+            # Do not check if a value is inside the feasible domain outside of active cells
+            mask = model.mesh.active_cell == 1
         elif key in model.serr_sigma_parameters.keys:
             arr = model.get_serr_sigma_parameters(key)
             low, upp = FEASIBLE_SERR_SIGMA_PARAMETERS[key]
+            # Check all values
+            mask = np.ones(arr.shape, dtype=bool)
         elif key in model.serr_mu_parameters.keys:
             arr = model.get_serr_mu_parameters(key)
             low, upp = FEASIBLE_SERR_MU_PARAMETERS[key]
+            # Check all values
+            mask = np.ones(arr.shape, dtype=bool)
         # % In case we have other kind of parameters. Should be unreachable.
         else:
             pass
 
-        low_arr = np.min(arr)
-        upp_arr = np.max(arr)
+        low_arr = np.min(arr, where=mask, initial=np.inf)
+        upp_arr = np.max(arr, where=mask, initial=-np.inf)
         if (low_arr + F_PRECISION) < value[0] or (upp_arr - F_PRECISION) > value[1]:
             raise ValueError(
                 f"Invalid bounds values for parameter '{key}'. Bounds domain [{value[0]}, {value[1]}] does "
@@ -1094,11 +1102,13 @@ def _standardize_simulation_return_options(model: Model, func_name: str, return_
 
 
 def _standardize_simulation_parameters_feasibility(model: Model):
+    mask = model.mesh.active_cell == 1
     for key in model.rr_parameters.keys:
         arr = model.get_rr_parameters(key)
         low, upp = FEASIBLE_RR_PARAMETERS[key]
-        low_arr = np.min(arr)
-        upp_arr = np.max(arr)
+        # Do not check if a value is inside the feasible domain outside of active cells
+        low_arr = np.min(arr, where=mask, initial=np.inf)
+        upp_arr = np.max(arr, where=mask, initial=-np.inf)
 
         if (low_arr + F_PRECISION) <= low or (upp_arr - F_PRECISION) >= upp:
             raise ValueError(
@@ -1109,8 +1119,9 @@ def _standardize_simulation_parameters_feasibility(model: Model):
     for key in model.rr_initial_states.keys:
         arr = model.get_rr_initial_states(key)
         low, upp = FEASIBLE_RR_INITIAL_STATES[key]
-        low_arr = np.min(arr)
-        upp_arr = np.max(arr)
+        # Do not check if a value is inside the feasible domain outside of active cells
+        low_arr = np.min(arr, where=mask, initial=np.inf)
+        upp_arr = np.max(arr, where=mask, initial=-np.inf)
 
         if (low_arr + F_PRECISION) <= low or (upp_arr - F_PRECISION) >= upp:
             raise ValueError(
@@ -1132,6 +1143,7 @@ def _standardize_simulation_parameters_feasibility(model: Model):
                 f"Invalid value for model serr_mu_parameter '{key}'. serr_mu_parameter domain "
                 f"[{low_arr}, {upp_arr}] is not included in the feasible domain ]{low}, {upp}["
             )
+
     for key in model.serr_sigma_parameters.keys:
         arr = model.get_serr_sigma_parameters(key)
         # % Skip if size == 0, i.e. no gauge

smash/fcore/forward/md_simulation.f90

@@ -230,7 +230,7 @@ subroutine simulation_checkpoint(setup, mesh, input_data, parameters, output, op
 
                 rr_parameters_inc = rr_parameters_inc + 5
                 rr_states_inc = rr_states_inc + 3
-                
+
                 ! 'gr6' module
             case ("gr6")
 
@@ -266,8 +266,8 @@ subroutine simulation_checkpoint(setup, mesh, input_data, parameters, output, op
                 checkpoint_variable%ac_rr_states(:, rr_states_inc + 4) = h4
 
                 rr_parameters_inc = rr_parameters_inc + 6
-                rr_states_inc = rr_states_inc + 4                
-                
+                rr_states_inc = rr_states_inc + 4
+
                 ! 'grd' module
             case ("grd")
 

smash/fcore/operator/md_gr_operator.f90

@@ -144,21 +144,21 @@ subroutine gr_exponential_transfer(pre, be, he, qe)
         real(sp), intent(inout) :: he
         real(sp), intent(out) :: qe
         real(sp) :: he_star, ar
-        
+
         he_star = he + pre
-        ar = he_star / be
+        ar = he_star/be
         if (ar .lt. -7._sp) then
-            qe = be * exp(ar)
+            qe = be*exp(ar)
         else if (ar .gt. 7._sp) then
-            qe = he_star + be / exp(ar)
+            qe = he_star + be/exp(ar)
         else
-            qe = be * log(exp(ar) + 1._sp)
+            qe = be*log(exp(ar) + 1._sp)
         end if
-        
+
         he = he_star - qe
-        
+
     end subroutine gr_exponential_transfer
-    
+
     subroutine gr4_time_step(setup, mesh, input_data, options, returns, time_step, ac_mlt, ac_ci, ac_cp, ac_ct, &
     & ac_kexc, ac_hi, ac_hp, ac_ht, ac_qt)
 
@@ -406,14 +406,14 @@ subroutine gr6_time_step(setup, mesh, input_data, options, returns, time_step, a
 
                 end if
 
-                prr = 0.6_sp * 0.9_sp * (pr + perc) + l
-                pre = 0.4_sp * 0.9_sp * (pr + perc) + l
-                prd = 0.1_sp * (pr + perc)
+                prr = 0.6_sp*0.9_sp*(pr + perc) + l
+                pre = 0.4_sp*0.9_sp*(pr + perc) + l
+                prd = 0.1_sp*(pr + perc)
 
                 call gr_transfer(5._sp, ac_prcp(k), prr, ac_ct(k), ac_ht(k), qr)
-                
+
                 call gr_exponential_transfer(pre, ac_be(k), ac_he(k), qe)
-                
+
                 qd = max(0._sp, prd + l)
 
                 ac_qt(k) = qr + qd + qe

smash/tests/test_constant.py

@@ -347,7 +347,7 @@ def test_default_bounds_parameters():
         (1e-6, 0.999999),  # % hi
         (1e-6, 0.999999),  # % hp
         (1e-6, 0.999999),  # % ht
-        (-1e3, 0),         # % he
+        (-1e3, 0),  # % he
         (1e-6, 0.999999),  # % ha
         (1e-6, 0.999999),  # % hc
         (1e-6, 0.999999),  # % hcl