Merge branch 'develop' into jax_export

AMICI-dev · Nov 11, 2024 · 4a5e7d2 · 4a5e7d2
2 parents 498681a + c3bbdd9
commit 4a5e7d2
Show file tree

Hide file tree

Showing 17 changed files with 190 additions and 87 deletions.
diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml
@@ -288,7 +288,8 @@ jobs:
 
     - name: Python tests
       run: |
-        scripts/run-python-tests.sh \
+        # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved
+        scripts/run-python-tests.sh -W ignore:: \
           test_pregenerated_models.py \
           test_splines_short.py \
           test_misc.py
@@ -350,7 +351,8 @@ jobs:
 
     - name: Python tests
       run: |
-        scripts/run-python-tests.sh \
+        # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved
+        scripts/run-python-tests.sh -W ignore:: \
           --ignore=test_pregenerated_models.py \
           --ignore=test_splines_short.py \
           --ignore=test_misc.py
diff --git a/.github/workflows/test_python_ver_matrix.yml b/.github/workflows/test_python_ver_matrix.yml
@@ -33,7 +33,7 @@ jobs:
           - python-version: '3.12'
             experimental: false
           - python-version: '3.13'
-            experimental: true
+            experimental: false
 
     steps:
     - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,63 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni
 
 ## v0.X Series
 
+### v0.28.0 (2024-11-11)
+
+**Breaking changes**
+
+* Changed the default steady-state method to `integrationOnly`
+  (by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2574)
+
+  The default mode for computing steady states and sensitivities at steady
+  state was changed to `integrationOnly`
+  (from previously `integrateIfNewtonFails`).
+
+  This was done for a more robust default behavior. For example, the evaluation
+  in https://doi.org/10.1371/journal.pone.0312148 shows that - at least for
+  some models - Newton's method may easily lead to physically impossible
+  solutions.
+
+  To keep the previous behavior, use:
+  ```python
+  amici_model.setSteadyStateComputationMode(amici.SteadyStateComputationMode.integrateIfNewtonFails)
+  amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)
+  ```
+
+**Fixes**
+
+* PEtab import: **Fixed potentially incorrect sensitivities** with
+  observable/state-dependent sigmas.
+  This was fixed for all cases amici can handle, others cases will now result
+  in `ValueError`s (https://github.com/AMICI-dev/AMICI/pull/2563).
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2562
+
+* Fixed potentially incorrect disabling of Newton's method
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2576
+
+* Fixed `ModelStateDerived` copy ctor, where previously dangling pointers
+  could lead to crashes in some situations
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2583
+
+* Added missing simulation status codes
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2560
+
+* Check for unsupported observable IDs in sigma expressions
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2563
+
+
+**Features**
+
+* Optional warning in `fill_in_parameters`
+
+  by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2578
+
+**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.27.0...v0.28.0
+
 ### v0.27.0 (2024-10-21)
 
 This release comes with an **updated version of the SUNDIALS package (7.1.1)** (https://github.com/AMICI-dev/AMICI/pull/2513).

diff --git a/codecov.yml b/codecov.yml
@@ -1,8 +1,10 @@
-fixes:
-  - "build/venv/lib/python3.9/site-packages/::python/"
-  - "build/venv/lib/python3.10/site-packages/::python/"
-  - "build/venv/lib/python3.11/site-packages/::python/"
+# see https://docs.codecov.com/docs/codecovyml-reference
 
+fixes:
+  # https://docs.codecov.com/docs/fixing-paths
+  - "build/venv/lib/python[0-9.]+/site-packages/::python/"
+  - "python/sdist/build/temp_amici/CMakeFiles/amici.dir/src/::src/"
+  - "build/CMakeFiles/amici.dir/src/::src/"
 codecov:
   require_ci_to_pass: yes
 
@@ -27,6 +29,8 @@ comment:
 ignore:
   - "tests/*"
   - "tests/**/*"
+  - "build/tests/**"
+  - "amici_models/**"
 
 flags:
   python:

diff --git a/include/amici/model.h b/include/amici/model.h
@@ -2063,12 +2063,12 @@ class Model : public AbstractModel, public ModelDimensions {
 
     /** method for steady-state computation */
     SteadyStateComputationMode steadystate_computation_mode_{
-        SteadyStateComputationMode::integrateIfNewtonFails
+        SteadyStateComputationMode::integrationOnly
     };
 
     /** method for steadystate sensitivities computation */
     SteadyStateSensitivityMode steadystate_sensitivity_mode_{
-        SteadyStateSensitivityMode::integrateIfNewtonFails
+        SteadyStateSensitivityMode::integrationOnly
     };
 
     /**

diff --git a/include/amici/model_state.h b/include/amici/model_state.h
@@ -146,71 +146,37 @@ struct ModelStateDerived {
         , dwdw_(other.dwdw_)
         , dwdx_hierarchical_(other.dwdx_hierarchical_)
         , dJydy_dense_(other.dJydy_dense_) {
-        // Update the SUNContext of all matrices
-        if (J_.data()) {
-            J_.get()->sunctx = sunctx_;
+        // Update the SUNContext of all SUNDIALS objects
+        J_.set_ctx(sunctx_);
+        JB_.set_ctx(sunctx_);
+        dxdotdw_.set_ctx(sunctx_);
+        dwdx_.set_ctx(sunctx_);
+        dwdp_.set_ctx(sunctx_);
+        M_.set_ctx(sunctx_);
+        MSparse_.set_ctx(sunctx_);
+        dfdx_.set_ctx(sunctx_);
+        dxdotdp_full.set_ctx(sunctx_);
+        dxdotdp_explicit.set_ctx(sunctx_);
+        dxdotdp_implicit.set_ctx(sunctx_);
+        dxdotdx_explicit.set_ctx(sunctx_);
+        dxdotdx_implicit.set_ctx(sunctx_);
+        dx_rdatadx_solver.set_ctx(sunctx_);
+        dx_rdatadtcl.set_ctx(sunctx_);
+        dtotal_cldx_rdata.set_ctx(sunctx_);
+        dxdotdp.set_ctx(sunctx_);
+
+        for (auto& dJydy : dJydy_) {
+            dJydy.set_ctx(sunctx_);
         }
-        if (JB_.data()) {
-            JB_.get()->sunctx = sunctx_;
+        for (auto& dwdp : dwdp_hierarchical_) {
+            dwdp.set_ctx(sunctx_);
         }
-        if (dxdotdw_.data()) {
-            dxdotdw_.get()->sunctx = sunctx_;
-        }
-        if (dwdx_.data()) {
-            dwdx_.get()->sunctx = sunctx_;
-        }
-        if (dwdp_.data()) {
-            dwdp_.get()->sunctx = sunctx_;
-        }
-        if (M_.data()) {
-            M_.get()->sunctx = sunctx_;
-        }
-        if (MSparse_.data()) {
-            MSparse_.get()->sunctx = sunctx_;
-        }
-        if (dfdx_.data()) {
-            dfdx_.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_full.data()) {
-            dxdotdp_full.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_explicit.data()) {
-            dxdotdp_explicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdp_implicit.data()) {
-            dxdotdp_implicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdx_explicit.data()) {
-            dxdotdx_explicit.get()->sunctx = sunctx_;
-        }
-        if (dxdotdx_implicit.data()) {
-            dxdotdx_implicit.get()->sunctx = sunctx_;
-        }
-        if (dx_rdatadx_solver.data()) {
-            dx_rdatadx_solver.get()->sunctx = sunctx_;
-        }
-        if (dx_rdatadtcl.data()) {
-            dx_rdatadtcl.get()->sunctx = sunctx_;
-        }
-        if (dtotal_cldx_rdata.data()) {
-            dtotal_cldx_rdata.get()->sunctx = sunctx_;
-        }
-        for (auto const& dwdp : dwdp_hierarchical_) {
-            if (dwdp.data()) {
-                dwdp.get()->sunctx = sunctx_;
-            }
-        }
-        for (auto const& dwdx : dwdx_hierarchical_) {
-            if (dwdx.data()) {
-                dwdx.get()->sunctx = sunctx_;
-            }
-        }
-        if (dwdw_.data()) {
-            dwdw_.get()->sunctx = sunctx_;
-        }
-        if (dJydy_dense_.data()) {
-            dJydy_dense_.get()->sunctx = sunctx_;
+        for (auto& dwdx : dwdx_hierarchical_) {
+            dwdx.set_ctx(sunctx_);
         }
+        sspl_.set_ctx(sunctx_);
+        dwdw_.set_ctx(sunctx_);
+        dJydy_dense_.set_ctx(sunctx_);
     }
 
     /**

diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h
@@ -506,6 +506,19 @@ class SUNMatrixWrapper {
      */
     SUNContext get_ctx() const;
 
+    /**
+     * @brief Set SUNContext
+     *
+     * Update the SUNContext of the wrapped SUNMatrix.
+     *
+     * @param ctx SUNDIALS context to set
+     */
+    void set_ctx(SUNContext ctx) {
+        if (matrix_) {
+            matrix_->sunctx = ctx;
+        }
+    }
+
   private:
     /**
      * @brief SUNMatrix to which all methods are applied

diff --git a/include/amici/vector.h b/include/amici/vector.h
@@ -413,6 +413,20 @@ class AmiVectorArray {
      */
     void copy(AmiVectorArray const& other);
 
+    /**
+     * @brief Set SUNContext
+     *
+     * If any AmiVector is non-empty, this changes the current SUNContext of the
+     * associated N_Vector. If empty, do nothing.
+     *
+     * @param ctx SUNDIALS context to set
+     */
+    void set_ctx(SUNContext ctx) {
+        for (auto& vec : vec_array_) {
+            vec.set_ctx(ctx);
+        }
+    }
+
   private:
     /** main data storage */
     std::vector<AmiVector> vec_array_;

diff --git a/python/sdist/amici/petab/conditions.py b/python/sdist/amici/petab/conditions.py
@@ -41,6 +41,7 @@ def fill_in_parameters(
     scaled_parameters: bool,
     parameter_mapping: ParameterMapping,
     amici_model: AmiciModel,
+    warn_unused: bool = True,
 ) -> None:
     """Fill fixed and dynamic parameters into the edatas (in-place).
 
@@ -59,9 +60,15 @@ def fill_in_parameters(
         Parameter mapping for all conditions.
     :param amici_model:
         AMICI model.
+    :param warn_unused:
+        Whether a warning should be emitted if not all problem parameters
+        were used. I.e., if there are parameters in `problem_parameters`
+        that are not in `parameter_mapping`.
     """
-    if unused_parameters := (
-        set(problem_parameters.keys()) - parameter_mapping.free_symbols
+    if warn_unused and (
+        unused_parameters := (
+            set(problem_parameters.keys()) - parameter_mapping.free_symbols
+        )
     ):
         warnings.warn(
             "The following problem parameters were not used: "
@@ -223,6 +230,7 @@ def create_parameterized_edatas(
     scaled_parameters: bool = False,
     parameter_mapping: ParameterMapping = None,
     simulation_conditions: pd.DataFrame | dict = None,
+    warn_unused: bool = True,
 ) -> list[amici.ExpData]:
     """Create list of :class:amici.ExpData objects with parameters filled in.
 
@@ -244,6 +252,11 @@ def create_parameterized_edatas(
     :param simulation_conditions:
         Result of :func:`petab.get_simulation_conditions`. Can be provided to
         save time if this has been obtained before.
+    :param warn_unused:
+        Whether a warning should be emitted if not all problem parameters
+        were used. I.e., if there are parameters in `problem_parameters`
+        that are not in `parameter_mapping` or in the generated parameter
+        mapping.
 
     :return:
         List with one :class:`amici.amici.ExpData` per simulation condition,
@@ -282,6 +295,7 @@ def create_parameterized_edatas(
         scaled_parameters=scaled_parameters,
         parameter_mapping=parameter_mapping,
         amici_model=amici_model,
+        warn_unused=warn_unused,
     )
 
     return edatas
@@ -387,7 +401,9 @@ def create_edatas(
 
     :return:
         List with one :class:`amici.amici.ExpData` per simulation condition,
-        with filled in timepoints and data.
+        with filled in timepoints and data, but without parameter values
+        (see :func:`create_parameterized_edatas` or
+        :func:`fill_in_parameters` for that).
     """
     if simulation_conditions is None:
         simulation_conditions = (

diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py
@@ -2052,12 +2052,26 @@ def _process_log_likelihood(
                 obs["reg_symbol"] = generate_regularization_symbol(obs_id)
 
         if not event_reg:
+            sigmas = {
+                obs_id: self._sympy_from_sbml_math(sigma)
+                for obs_id, sigma in sigmas.items()
+            }
+            obs_syms = set(self.symbols[obs_symbol].keys())
+            for obs_id in self.symbols[obs_symbol]:
+                if (sigma := sigmas.get(str(obs_id))) and sigma.has(
+                    *(obs_syms - {obs_id})
+                ):
+                    raise ValueError(
+                        f"Sigma expression for {obs_id} ({sigma}) must not "
+                        f"contain any observable symbols other than {obs_id}."
+                    )
+            # check that only the corresponding observable ID is used in the
+            #  sigma formula, but not any other observable ID
+            #  https://github.com/AMICI-dev/AMICI/issues/2561
             self.symbols[sigma_symbol] = {
                 symbol_with_assumptions(f"sigma_{obs_id}"): {
                     "name": f'sigma_{obs["name"]}',
-                    "value": self._sympy_from_sbml_math(
-                        sigmas.get(str(obs_id), "1.0")
-                    ),
+                    "value": sigmas.get(str(obs_id), sp.Float(1.0)),
                 }
                 for obs_id, obs in self.symbols[obs_symbol].items()
             }

diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py
@@ -102,6 +102,7 @@ def get_results(
     sensi_order=0,
     sensi_meth=amici.SensitivityMethod.forward,
     sensi_meth_preeq=amici.SensitivityMethod.forward,
+    stst_mode=amici.SteadyStateComputationMode.integrateIfNewtonFails,
     stst_sensi_mode=amici.SteadyStateSensitivityMode.newtonOnly,
     reinitialize_states=False,
 ):
@@ -115,6 +116,7 @@ def get_results(
     solver.setSensitivityMethodPreequilibration(sensi_meth_preeq)
     solver.setSensitivityMethod(sensi_meth)
     model.setSteadyStateSensitivityMode(stst_sensi_mode)
+    model.setSteadyStateComputationMode(stst_mode)
     if edata is None:
         model.setTimepoints(np.linspace(0, 5, 101))
     else:

diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py
@@ -18,7 +18,12 @@
 def preeq_fixture(pysb_example_presimulation_module):
     model = pysb_example_presimulation_module.getModel()
     model.setReinitializeFixedParameterInitialStates(True)
-
+    model.setSteadyStateComputationMode(
+        amici.SteadyStateComputationMode.integrateIfNewtonFails
+    )
+    model.setSteadyStateSensitivityMode(
+        amici.SteadyStateSensitivityMode.integrateIfNewtonFails
+    )
     solver = model.getSolver()
     solver.setSensitivityOrder(amici.SensitivityOrder.first)
     solver.setSensitivityMethod(amici.SensitivityMethod.forward)