Fix ReturnData::{preeq_wrms,posteq_wrms} with FSA and check_sensi_steadystate_conv_=True (#2071)

Documentation says, that ReturnData::{preeq_wrms,posteq_wrms} refers to the `weighted root-mean-square of the rhs when steadystate was reached`, independent of `Solver::check_sensi_steadystate_conv_`, but this wasn't what happened. Now it is.

Closes #2066 

Furthermore, this prevents unnecessary simulation if we are already in a steadystate when entering `SteadystateProblem::runSteadystateSimulation`.

Author: dweindl

python/tests/test_preequilibration.py

@@ -2,11 +2,14 @@
 
 import itertools
 
-import amici
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
+
+import amici
 from test_pysb import get_data
 
+
 @pytest.fixture
 def preeq_fixture(pysb_example_presimulation_module):
 
@@ -115,11 +118,12 @@ def test_manual_preequilibration(preeq_fixture):
         assert rdata_sim.status == amici.AMICI_SUCCESS
 
         for variable in ['x', 'sx']:
-            assert np.isclose(
+            assert_allclose(
                 rdata_auto[variable],
                 rdata_sim[variable],
-                1e-6, 1e-6
-            ).all(), dict(pscale=pscale, plist=plist, variable=variable)
+                atol=1e-6, rtol=1e-6,
+                err_msg=str(dict(pscale=pscale, plist=plist, variable=variable))
+            )
 
 
 def test_parameter_reordering(preeq_fixture):
@@ -172,11 +176,12 @@ def test_data_replicates(preeq_fixture):
     rdata_double = amici.runAmiciSimulation(model, solver, edata)
 
     for variable in ['llh', 'sllh']:
-        assert np.isclose(
-            2*rdata_single[variable],
+        assert_allclose(
+            2 * rdata_single[variable],
             rdata_double[variable],
-            1e-6, 1e-6
-        ).all(), dict(variable=variable, sensi_meth=sensi_meth)
+            atol=1e-6, rtol=1e-6,
+            err_msg=str(dict(variable=variable, sensi_meth=sensi_meth))
+        )
 
 
 def test_parameter_in_expdata(preeq_fixture):

python/tests/test_pysb.py

@@ -19,7 +19,7 @@
 
 from amici.gradient_check import check_derivatives
 from amici.testing import skip_on_valgrind, TemporaryDirectoryWinSafe
-
+from numpy.testing import assert_allclose
 
 @skip_on_valgrind
 def test_compare_to_sbml_import(pysb_example_presimulation_module,
@@ -82,10 +82,11 @@ def test_compare_to_sbml_import(pysb_example_presimulation_module,
         elif np.isnan(rdata_pysb[field]).all():
             assert np.isnan(rdata_sbml[field]).all(), field
         else:
-            assert np.isclose(
+            assert_allclose(
                 rdata_sbml[field], rdata_pysb[field],
-                atol=1e-6, rtol=1e-6
-            ).all(), field
+                atol=1e-6, rtol=1e-6,
+                err_msg=field
+            )
 
 
 pysb_models = [

src/steadystateproblem.cpp

@@ -665,9 +665,6 @@ void SteadystateProblem::runSteadystateSimulation(
     if (backward)
         sensitivityFlag = SensitivityMethod::adjoint;
 
-    /* If run after Newton's method checks again if it converged */
-    wrms_ = getWrms(model, sensitivityFlag);
-
     int &sim_steps = backward ? numstepsB_ : numsteps_.at(1);
 
     int convergence_check_frequency = 1;
@@ -676,10 +673,30 @@ void SteadystateProblem::runSteadystateSimulation(
         convergence_check_frequency = 25;
 
     while (true) {
+        if (sim_steps % convergence_check_frequency == 0) {
+            // Check for convergence (already before simulation, since we might
+            // start in steady state)
+            wrms_ = getWrms(model, sensitivityFlag);
+            if (wrms_ < conv_thresh) {
+                if(check_sensi_conv_
+                    && sensitivityFlag == SensitivityMethod::forward) {
+                    updateSensiSimulation(solver);
+                    // getWrms needs to be called before getWrmsFSA
+                    // such that the linear system is prepared for newton-type
+                    // convergence check
+                    if (getWrmsFSA(model) < conv_thresh)
+                        break; // converged
+                } else {
+                    break; // converged
+                }
+            }
+        }
+
         /* check for maxsteps  */
         if (sim_steps >= solver.getMaxSteps()) {
             throw IntegrationFailure(AMICI_TOO_MUCH_WORK, state_.t);
         }
+
         /* increase counter */
         sim_steps++;
         /* One step of ODE integration
@@ -698,21 +715,6 @@ void SteadystateProblem::runSteadystateSimulation(
                                   xQ_);
             flagUpdatedState();
         }
-
-        /* Check for convergence */
-        if (sim_steps % convergence_check_frequency == 0) {
-            wrms_ = getWrms(model, sensitivityFlag);
-            /* getWrms needs to be called before getWrmsFSA such that the linear
-             system is prepared for newton type convergence check */
-            if (wrms_ < conv_thresh && check_sensi_conv_ &&
-                sensitivityFlag == SensitivityMethod::forward) {
-                updateSensiSimulation(solver);
-                wrms_ = getWrmsFSA(model);
-            }
-        }
-
-        if (wrms_ < conv_thresh)
-            break; // converged
     }
 
     // if check_sensi_conv_ is deactivated, we still have to update sensis