Use Model.setAlwaysCheckFinite(True) in debug builds (AMICI-dev#2262)

* Use Model.setAlwaysCheckFinite(True) in debug builds So that we automatically use the additional checks in our tests. Closes AMICI-dev#2235 * Fix sparse-`unravel_index` for the case where no indices/ptrs have been set (this is the case for matlab-generated models).
dweindl · Jan 7, 2024 · 3d01ebb · 3d01ebb
1 parent 66f4ba5
commit 3d01ebb
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 19 deletions.
diff --git a/include/amici/model.h b/include/amici/model.h
@@ -2020,7 +2020,11 @@ class Model : public AbstractModel, public ModelDimensions {
      * Indicates whether the result of every call to `Model::f*` should be
      * checked for finiteness
      */
+#ifdef NDEBUG
     bool always_check_finite_{false};
+#else
+    bool always_check_finite_{true};
+#endif
 
     /** indicates whether sigma residuals are to be added for every datapoint */
     bool sigma_res_{false};

diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py
@@ -68,10 +68,7 @@ def test_copy_constructors(pysb_example_presimulation_module):
 model_instance_settings0 = {
     # setting name: [default value, custom value]
     "AddSigmaResiduals": [False, True],
-    "AlwaysCheckFinite": [
-        False,
-        True,
-    ],
+    "AlwaysCheckFinite": [False, True],
     # Skipped due to model dependency in `'InitialStates'`.
     "FixedParameters": None,
     "InitialStates": [
@@ -132,6 +129,13 @@ def test_model_instance_settings(pysb_example_presimulation_module):
     i_getter = 0
     i_setter = 1
 
+    # the default setting for AlwaysCheckFinite depends on whether the amici
+    # extension has been built in debug mode
+    model_instance_settings0["AlwaysCheckFinite"] = [
+        model0.getAlwaysCheckFinite(),
+        not model0.getAlwaysCheckFinite(),
+    ]
+
     # All settings are tested.
     assert set(model_instance_settings0) == set(
         amici.swig_wrappers.model_instance_settings

diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp
@@ -793,20 +793,25 @@ unravel_index(sunindextype i, SUNMatrix m) {
     }
 
     if (mat_id == SUNMATRIX_SPARSE) {
-        gsl_ExpectsDebug(i < SM_NNZ_S(m));
-        sunindextype row = SM_INDEXVALS_S(m)[i];
-        sunindextype i_colptr = 0;
-        while (SM_INDEXPTRS_S(m)[i_colptr] < SM_NNZ_S(m)) {
-            if (SM_INDEXPTRS_S(m)[i_colptr + 1] > i) {
-                sunindextype col = i_colptr;
-                gsl_EnsuresDebug(row >= 0);
-                gsl_EnsuresDebug(row < SM_ROWS_S(m));
-                gsl_EnsuresDebug(col >= 0);
-                gsl_EnsuresDebug(col < SM_COLUMNS_S(m));
-                return {row, col};
-            }
-            ++i_colptr;
-        }
+        auto nnz = SM_NNZ_S(m);
+        auto ncols = SM_COLUMNS_S(m);
+        auto index_vals = SM_INDEXVALS_S(m);
+        auto index_ptrs = SM_INDEXPTRS_S(m);
+        gsl_ExpectsDebug(i < nnz);
+        sunindextype row = index_vals[i];
+        sunindextype col = 0;
+        while (col < ncols && index_ptrs[col + 1] <= i)
+            ++col;
+
+        // This can happen if indexvals / indexptrs haven't been set.
+        if(col == ncols)
+            return {-1, -1};
+
+        gsl_EnsuresDebug(row >= 0);
+        gsl_EnsuresDebug(row < SM_ROWS_S(m));
+        gsl_EnsuresDebug(col >= 0);
+        gsl_EnsuresDebug(col < ncols);
+        return {row, col};
     }
 
     throw amici::AmiException("Unimplemented SUNMatrix type for unravel_index");

diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp
@@ -689,7 +689,7 @@ TEST(UnravelIndex, UnravelIndexSunMatSparse)
     // [2, 0]
     // data [1, 2, 3]
     // colptrs [0, 2, 3]
-    // rowidxs [2, 3, 1]
+    // rowidxs [2, 3, 0]
     D.set_data(0, 0, 0);
     D.set_data(1, 0, 0);
     D.set_data(2, 0, 1);
@@ -708,6 +708,16 @@ TEST(UnravelIndex, UnravelIndexSunMatSparse)
     SUNMatDestroy(S);
 }
 
+
+TEST(UnravelIndex, UnravelIndexSunMatSparseMissingIndices)
+{
+    // Sparse matrix without any indices set
+    SUNMatrixWrapper mat = SUNMatrixWrapper(2, 3, 2, CSC_MAT);
+    EXPECT_EQ(unravel_index(0, mat.get()), std::make_pair((sunindextype) -1, (sunindextype) -1));
+    EXPECT_EQ(unravel_index(1, mat.get()), std::make_pair((sunindextype) -1, (sunindextype) -1));
+}
+
+
 TEST(ReturnCodeToStr, ReturnCodeToStr)
 {
     EXPECT_EQ("AMICI_SUCCESS", simulation_status_to_str(AMICI_SUCCESS));