simplify steadystate logic

include/amici/defines.h

@@ -171,7 +171,8 @@ enum class NonlinearSolverIteration {
 /** Sensitivity computation mode in steadyStateProblem */
 enum class SteadyStateSensitivityMode {
  newtonOnly,
- simulationFSA
+ simulationFSA,
+ mixed,
 };
 
 /** State in which the steady state computation finished */
@@ -185,13 +186,6 @@ enum class SteadyStateStatus {
  success = 1
 };
 
-/** Context for which the sensitivity flag should be computed */
-enum class SteadyStateContext {
- newtonSensi = 0,
- sensiStorage = 1,
- solverCreation = 2
-};
-
 /** Damping factor flag for the Newton method */
 enum class NewtonDampingFactorMode {
  off = 0,

include/amici/model.h

@@ -1697,7 +1697,9 @@ class Model : public AbstractModel, public ModelDimensions {
  * flag indicating whether steadystate sensitivities are to be computed
  * via FSA when steadyStateSimulation is used
  */
- SteadyStateSensitivityMode steadystate_sensitivity_mode_ {SteadyStateSensitivityMode::newtonOnly};
+ SteadyStateSensitivityMode steadystate_sensitivity_mode_ {
+ SteadyStateSensitivityMode::mixed
+ };
 
  /**
  * Indicates whether the result of every call to `Model::f*` should be

include/amici/steadystateproblem.h

@@ -128,17 +128,68 @@ class SteadystateProblem {
  status) const;
 
  /**
- * @brief Checks depending on the status of the Newton solver,
- * solver settings, and the model, whether state sensitivities
- * still need to be computed via a linear system solve or stored
- * @param model pointer to the model object
- * @param solver pointer to the solver object
+ * @brief Checks whether steadystate may be computed via Newton's method
+ * @param model model instance
+ * @param solver solver instance
+ * @return true if allowed, false otherwise
+ */
+ bool checkNewtonAllowed(const Model &model, const Solver &solver) const;
+
+ /**
+ * @brief Checks whether FSA sensis provided by initialization or simulation solver can be used
+ * @param model model instance
+ * @param solver solver instance
  * @param it integer with the index of the current time step
- * @param context SteadyStateContext giving the situation for the flag
- * @return flag telling how to process state sensitivities
+ * @return true if correct, false otherwise
  */
- bool getSensitivityFlag(const Model *model, const Solver *solver, int it,
- SteadyStateContext context);
+ bool checkUseFSASensis(const Model &model,
+ const Solver &solver,
+ const int it) const;
+
+ /**
+ * @brief Checks whether settings are all compatible
+ * @param model model instance
+ * @param solver solver instance
+ * @param it integer with the index of the current time step
+ * @return true if correct, false otherwise
+ */
+ void checkSettingsCompatibility(const Model &model,
+ const Solver &solver,
+ const int it) const;
+
+ /**
+ * @brief Checks whether steadystate sensitivities need to be computed according to the implicit
+ * function theorem as solution to a linear system
+ * @param model model instance
+ * @param solver solver instance
+ * @param it integer with the index of the current time step
+ * @return true if they need to be computed, false otherwise
+ */
+ bool checkComputeImplicitSensis(const Model &model,
+ const Solver &solver,
+ const int it) const;
+
+ /**
+ * @brief Checks whether the simulator needs to compute sensi via FSA
+ * @param model model instance
+ * @param solver solver instance
+ * @param it integer with the index of the current time step
+ * @return true if computations is necessary, false otherwise
+ */
+ bool checkComputeFSASensis(const Model &model,
+ const Solver &solver,
+ const int it) const;
+
+ /**
+ * @brief Checks whether newton or FSA sensis were computed
+ * @param model model instance
+ * @param solver solver instance
+ * @param it integer with the index of the current time step
+ * @return true if sensis were computed, false otherwise
+ */
+ bool checkUseNewtonOrFSASensis(const Model &model,
+ const Solver &solver,
+ const int it) const;
 
  /**
  * @brief Computes the weighted root mean square of xdot
@@ -413,7 +464,9 @@ class SteadystateProblem {
  /** stores diagnostic information about execution success of the different
  * approaches [newton, simulation, newton] (length = 3)
  */
- std::vector<SteadyStateStatus> steady_state_status_;
+ std::vector<SteadyStateStatus> steady_state_status_ {
+ std::vector<SteadyStateStatus>(3, SteadyStateStatus::not_run)
+ };
 };
 
 } // namespace amici

python/amici/numpy.py

@@ -204,11 +204,11 @@ def __init__(self, rdata: Union[ReturnDataPtr, ReturnData]):
  'w': [rdata.nt, rdata.nw],
  'xdot': [rdata.nx_solver],
  'preeq_numlinsteps': [rdata.newton_maxsteps, 2],
- 'preeq_numsteps': [1, 3],
- 'preeq_status': [1, 3],
+ 'preeq_numsteps': [3],
+ 'preeq_status': [3],
  'posteq_numlinsteps': [rdata.newton_maxsteps, 2],
- 'posteq_numsteps': [1, 3],
- 'posteq_status': [1, 3],
+ 'posteq_numsteps': [3],
+ 'posteq_status': [3],
  'numsteps': [rdata.nt],
  'numrhsevals': [rdata.nt],
  'numerrtestfails': [rdata.nt],

python/tests/test_compare_conservation_laws_sbml.py

@@ -41,7 +41,7 @@ def generate_models():
  sbml_importer = amici.SbmlImporter(sbml_file)
 
  # Name of the model that will also be the name of the python module
- model_name =  model_output_dir ='model_constant_species'
+ model_name = model_output_dir ='model_constant_species'
  model_name_cl = model_output_dir_cl = 'model_constant_species_cl'
 
  # Define constants, observables, sigmas
@@ -101,71 +101,79 @@ def get_results(model, edata=None, sensi_order=0,
 
 
 def test_compare_conservation_laws_sbml(edata_fixture):
- # first, create the modek
+ # first, create the model
  model_with_cl, model_without_cl = generate_models()
 
  # ----- compare simulations wo edata, sensi = 0, states ------------------
  # run simulations
  rdata_cl = get_results(model_with_cl)
- assert rdata_cl['status'] == amici.AMICI_SUCCESS
+ assert rdata_cl.status == amici.AMICI_SUCCESS
  rdata = get_results(model_without_cl)
- assert rdata['status'] == amici.AMICI_SUCCESS
+ assert rdata.status == amici.AMICI_SUCCESS
 
  # compare state trajectories
- assert np.isclose(rdata['x'], rdata_cl['x']).all()
+ assert np.isclose(rdata.x, rdata_cl.x).all()
 
  # ----- compare simulations wo edata, sensi = 1, states and sensis -------
  # run simulations
  rdata_cl = get_results(model_with_cl, sensi_order=1)
- assert rdata_cl['status'] == amici.AMICI_SUCCESS
+ assert rdata_cl.status == amici.AMICI_SUCCESS
  rdata = get_results(model_without_cl, sensi_order=1)
- assert rdata['status'] == amici.AMICI_SUCCESS
+ assert rdata.status == amici.AMICI_SUCCESS
 
  # compare state trajectories
  for field in ['x', 'sx']:
  assert np.isclose(rdata[field], rdata_cl[field]).all(), field
 
- # ----- compare simulations wo edata, sensi = 0, states and sensis -------
+ # ----- compare simulations w edata, sensi = 0, states -------
  model_without_cl.setSteadyStateSensitivityMode(
  amici.SteadyStateSensitivityMode.simulationFSA
  )
 
  # run simulations
  edata, _, _ = edata_fixture
  rdata_cl = get_results(model_with_cl, edata=edata)
- assert rdata_cl['status'] == amici.AMICI_SUCCESS
+ assert rdata_cl.status == amici.AMICI_SUCCESS
+ # check that steady state computation succeeded by Newton in reduced model
+ assert np.array_equal(rdata_cl.preeq_status, [1, 0, 0])
+
  rdata = get_results(model_without_cl, edata=edata)
- assert rdata['status'] == amici.AMICI_SUCCESS
+ assert rdata.status == amici.AMICI_SUCCESS
+ # check that steady state computation succeeded by Newton in full model
+ assert np.array_equal(rdata_cl.preeq_status, [1, 0, 0])
 
  # compare preequilibrated states
  for field in ['x', 'x_ss', 'llh']:
  assert np.isclose(rdata[field], rdata_cl[field]).all(), field
 
- # ----- compare simulations wo edata, sensi = 1, states and sensis -------
+ # ----- compare simulations w edata, sensi = 1, states and sensis -------
 
  # run simulations
  rdata_cl = get_results(model_with_cl, edata=edata, sensi_order=1)
- assert rdata_cl['status'] == amici.AMICI_SUCCESS
+ assert rdata_cl.status == amici.AMICI_SUCCESS
+ # check that steady state computation succeeded by Newton in reduced
+ # model
+ assert np.array_equal(rdata_cl.preeq_status, [1, 0, 0])
+
  rdata = get_results(model_without_cl, edata=edata, sensi_order=1)
- assert rdata['status'] == amici.AMICI_SUCCESS
- # check that steady state computation succeeded only by sim in full model
- assert (rdata['preeq_status'] == np.array([-3, 1, 0])).all()
- # check that steady state computation succeeded by Newton in reduced model
- assert (rdata_cl['preeq_status'] == np.array([1, 0, 0])).all()
+ assert rdata.status == amici.AMICI_SUCCESS
+ # check that steady state computation succeeded by sim in full model
+ assert np.array_equal(rdata.preeq_status, [0, 1, 0])
 
  # compare state sensitivities with edata and preequilibration
  for field in ['x', 'x_ss', 'sx', 'llh', 'sllh']:
  assert np.isclose(rdata[field], rdata_cl[field]).all(), field
 
- # ----- check failure st.st. sensi computation if run wo CLs -------------
- # check failure of steady state senistivity computation if run wo CLs
+ # ----- check failure of sensi computation if run wo CLs -------------
+ # check failure of steady state sensitivity computation if run wo CLs
  model_without_cl.setSteadyStateSensitivityMode(
  amici.SteadyStateSensitivityMode.newtonOnly
  )
  with warnings.catch_warnings():
  warnings.filterwarnings("ignore")
  rdata = get_results(model_without_cl, edata=edata, sensi_order=1)
- assert rdata['status'] == amici.AMICI_ERROR
+ assert rdata.status == amici.AMICI_ERROR
+ assert np.array_equal(rdata.preeq_status, [-3, 1, 0])
 
 
 def test_adjoint_pre_and_post_equilibration(edata_fixture):
@@ -198,9 +206,9 @@ def test_adjoint_pre_and_post_equilibration(edata_fixture):
  reinitialize_states=reinit)
 
  # assert all are close
- assert np.isclose(rff_cl['sllh'], rfa_cl['sllh']).all()
- assert np.isclose(rfa_cl['sllh'], raa_cl['sllh']).all()
- assert np.isclose(raa_cl['sllh'], rff_cl['sllh']).all()
+ assert np.isclose(rff_cl.sllh, rfa_cl.sllh).all()
+ assert np.isclose(rfa_cl.sllh, raa_cl.sllh).all()
+ assert np.isclose(raa_cl.sllh, rff_cl.sllh).all()
 
  # --- compare fully adjoint approach to simulation with singular Jacobian ----
  raa = get_results(model, edata=edata, sensi_order=1,
@@ -209,4 +217,19 @@ def test_adjoint_pre_and_post_equilibration(edata_fixture):
  reinitialize_states=reinit)
 
  # assert gradients are close (quadrature tolerances are laxer)
- assert np.isclose(raa_cl['sllh'], raa['sllh'], 1e-5, 1e-5).all()
+ assert np.isclose(raa_cl.sllh, raa.sllh, 1e-5, 1e-5).all()
+
+ # this needs to fail due to failed RHS factorization
+ raf = get_results(model, edata=edata, sensi_order=1,
+ sensi_meth=amici.SensitivityMethod.adjoint,
+ sensi_meth_preeq=amici.SensitivityMethod.forward,
+ reinitialize_states=reinit)
+
+ assert raf.status == amici.AMICI_ERROR
+ if edata.fixedParametersPreequilibration:
+ assert np.array_equal(raf.preeq_status, [-3, 1, 0])
+ assert np.array_equal(raf.posteq_status, [0, 0, 0])
+ else:
+ assert np.array_equal(raf.preeq_status, [0, 0, 0])
+ assert np.array_equal(raf.posteq_status, [-3, 1, 0])
+