Change default steady-state method to integrationOnly (#2574)

**Breaking change**

Change the default mode for computing steady states and sensitivities
at steady state to `integrationOnly` (from previously `integrateIfNewtonFails`).

This is done for a more robust default behaviour. 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 behaviour, use:
```python
amici_model.setSteadyStateComputationMode(amici.SteadyStateComputationMode.integrateIfNewtonFails)
amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails)
```

Closes #2571.

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
     };
 
     /**

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:

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)

python/tests/test_pregenerated_models.py

@@ -63,6 +63,13 @@ def test_pregenerated_model(sub_test, case):
     amici.readModelDataFromHDF5(
         options_file, model.get(), f"/{sub_test}/{case}/options"
     )
+    if model_name == "model_steadystate":
+        model.setSteadyStateComputationMode(
+            amici.SteadyStateComputationMode.integrateIfNewtonFails
+        )
+        model.setSteadyStateSensitivityMode(
+            amici.SteadyStateSensitivityMode.integrateIfNewtonFails
+        )
     amici.readSolverSettingsFromHDF5(
         options_file, solver.get(), f"/{sub_test}/{case}/options"
     )

python/tests/test_swig_interface.py

@@ -105,12 +105,12 @@ def test_copy_constructors(pysb_example_presimulation_module):
     # `pysb_example_presimulation_module.getModel()`.
     "StateIsNonNegative": None,
     "SteadyStateComputationMode": [
-        2,
-        1,
+        amici.SteadyStateComputationMode.integrationOnly,
+        amici.SteadyStateComputationMode.integrateIfNewtonFails,
     ],
     "SteadyStateSensitivityMode": [
-        2,
-        1,
+        amici.SteadyStateSensitivityMode.integrationOnly,
+        amici.SteadyStateSensitivityMode.integrateIfNewtonFails,
     ],
     ("t0", "setT0"): [
         0.0,