simplify steadystate logic

include/amici/steadystateproblem.h

@@ -128,17 +128,45 @@ 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
- * @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
+ * @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);
+
+ /**
+ * @brief Checks whether the simulation solver or initialization provided sensitivities are correct
+ * @param model model instance
+ * @param solver solver instance
+ * @return true if correct, false otherwise
+ */
+ bool checkUseFSASensis(const Model &model, const Solver &solver);
+
+ /**
+ * @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
+ * @return true if they need to be computed, false otherwise
+ */
+ bool checkComputeImplicitSensis(const Model &model, const Solver &solver);
+
+ /**
+ * @brief Checks whether the simulator needs to compute sensi via FSA
+ * @param model model instance
+ * @param solver solver instance
+ * @return true if computations is necessary, false otherwise
+ */
+ bool checkComputeFSASensis(const Model &model, const Solver &solver);
+
+ /**
+ * @brief Checks whether newton or FSA sensis were computed
+ * @param model model instance
+ * @param solver solver instance
+ * @return true if sensis were computed, false otherwise
  */
- bool getSensitivityFlag(const Model *model, const Solver *solver, int it,
- SteadyStateContext context);
+ bool checkUseNewtonOrFSASensis(const Model &model, const Solver &solver);
 
  /**
  * @brief Computes the weighted root mean square of xdot

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])
+

python/tests/test_preequilibration.py

@@ -91,10 +91,10 @@ def test_manual_preequilibration(preeq_fixture):
  assert rdata_preeq.status == amici.AMICI_SUCCESS
 
  # manual reinitialization + presimulation
- x0 = rdata_preeq['x'][0, :]
+ x0 = rdata_preeq.x[0, :]
  x0[1] = edata_presim.fixedParameters[0]
  x0[2] = edata_presim.fixedParameters[1]
- sx0 = rdata_preeq['sx'][0, :, :]
+ sx0 = rdata_preeq.sx[0, :, :]
  sx0[:, 1] = 0
  sx0[:, 2] = 0
  model.setInitialStates(x0)
@@ -103,10 +103,10 @@ def test_manual_preequilibration(preeq_fixture):
  assert rdata_presim.status == amici.AMICI_SUCCESS
 
  # manual reinitialization + simulation
- x0 = rdata_presim['x'][0, :]
+ x0 = rdata_presim.x[0, :]
  x0[1] = edata_sim.fixedParameters[0]
  x0[2] = edata_sim.fixedParameters[1]
- sx0 = rdata_presim['sx'][0, :, :]
+ sx0 = rdata_presim.sx[0, :, :]
  sx0[:, 1] = 0
  sx0[:, 2] = 0
  model.setInitialStates(x0)
@@ -136,8 +136,8 @@ def test_parameter_reordering(preeq_fixture):
 
  for ip, p_index in enumerate(plist):
  assert np.isclose(
- rdata_ordered['sx'][:, p_index, :],
- rdata_reordered['sx'][:, ip, :],
+ rdata_ordered.sx[:, p_index, :],
+ rdata_reordered.sx[:, ip, :],
  1e-6, 1e-6
  ).all(), plist
 
@@ -196,8 +196,8 @@ def test_parameter_in_expdata(preeq_fixture):
  edata.sx0 = model.getInitialStateSensitivities()
 
  # perturb model initial states
- model.setInitialStates(rdata['x_ss'] * 4)
- model.setInitialStateSensitivities(rdata['sx_ss'].flatten() / 2)
+ model.setInitialStates(rdata.x_ss * 4)
+ model.setInitialStateSensitivities(rdata.sx_ss.flatten() / 2)
 
  # set ExpData plist
  edata.plist = model.getParameterList()
@@ -241,29 +241,42 @@ def test_raise_presimulation_with_adjoints(preeq_fixture):
  model, solver, edata, edata_preeq, \
  edata_presim, edata_sim, pscales, plists = preeq_fixture
 
+ model.setSteadyStateSensitivityMode(
+ amici.SteadyStateSensitivityMode.newtonOnly
+ )
+
  # preequilibration and presimulation with adjoints:
  # this needs to fail unless we remove presimulation
  solver.setSensitivityMethod(amici.SensitivityMethod.adjoint)
 
  rdata = amici.runAmiciSimulation(model, solver, edata)
- assert rdata['status'] == amici.AMICI_ERROR
+ assert rdata.status == amici.AMICI_ERROR
+ assert np.array_equal(rdata.preeq_status, [-2, 1, 0])
+ assert np.array_equal(rdata.posteq_status, [0, 0, 0])
 
  # presimulation and postequilibration with adjoints:
  # this also needs to fail
  y = edata.getObservedData()
  stdy = edata.getObservedDataStdDev()
 
- # add infty timepoint
+ # add infty timepoint for postequilibration
  ts = np.hstack([*edata.getTimepoints(), np.inf])
  edata.setTimepoints(sorted(ts))
  edata.setObservedData(np.hstack([y, y[0]]))
  edata.setObservedDataStdDev(np.hstack([stdy, stdy[0]]))
- edata.t_presim = 0
- edata.fixedParametersPresimulation = ()
+
+ rdata = amici.runAmiciSimulation(model, solver, edata)
+ assert rdata.status == amici.AMICI_ERROR
+ assert np.array_equal(rdata.preeq_status, [-2, 1, 0])
+ assert np.array_equal(rdata.posteq_status, [0, 0, 0])
 
  # no presim any more, this should work
+ edata.t_presim = 0
+ edata.fixedParametersPresimulation = ()
  rdata = amici.runAmiciSimulation(model, solver, edata)
- assert rdata['status'] == amici.AMICI_SUCCESS
+ assert rdata.status == amici.AMICI_SUCCESS
+ assert np.array_equal(rdata.preeq_status, [-2, 1, 0])
+ assert np.array_equal(rdata.posteq_status, [-2, 1, 0])
 
 
 def test_equilibration_methods_with_adjoints(preeq_fixture):
@@ -276,7 +289,7 @@ def test_equilibration_methods_with_adjoints(preeq_fixture):
  edata.t_presim = 0.0
  edata.fixedParametersPresimulation = ()
 
- # add infty timepoint
+ # add infty timepoint for postequilibration
  y = edata.getObservedData()
  stdy = edata.getObservedDataStdDev()
  ts = np.hstack([*edata.getTimepoints(), np.inf])
@@ -296,13 +309,29 @@ def test_equilibration_methods_with_adjoints(preeq_fixture):
  equil_meth, sensi_meth = setting
  model.setSteadyStateSensitivityMode(equil_meth)
  solver.setSensitivityMethod(sensi_meth)
- solver.setNewtonMaxSteps(0)
 
  # add rdatas
  rdatas[setting] = amici.runAmiciSimulation(model, solver, edata)
  # assert successful simulation
+ if (sensi_meth, equil_meth) == (
+ amici.SensitivityMethod.adjoint,
+ amici.SteadyStateSensitivityMode.simulationFSA
+ ):
+ assert rdatas[setting].status == amici.AMICI_ERROR
+ # preeq works since we use a seperate solver
+ assert np.array_equal(rdatas[setting].preeq_status, [0, 1, 0])
+ # posteq doesn't work since we would need to activate sensis in 
+ # the solver
+ assert np.array_equal(rdatas[setting].posteq_status, [0, 0, 0])
+ else:
+ assert rdatas[setting].status == amici.AMICI_SUCCESS
+ if equil_meth == amici.SteadyStateSensitivityMode.newtonOnly:
+ target_status = [-2, 1, 0]
+ else:
+ target_status = [0, 1, 0]
+ assert np.array_equal(rdatas[setting].posteq_status, target_status)
+ assert np.array_equal(rdatas[setting].preeq_status, target_status)
 
- assert rdatas[setting]['status'] == amici.AMICI_SUCCESS
 
  for setting1, setting2 in itertools.product(settings, settings):
  # assert correctness of result
@@ -341,16 +370,14 @@ def test_newton_solver_equilibration(preeq_fixture):
  sensi_meth = amici.SensitivityMethod.forward
  solver.setSensitivityMethod(sensi_meth)
  model.setSteadyStateSensitivityMode(equil_meth)
- if equil_meth == amici.SteadyStateSensitivityMode.newtonOnly:
- solver.setNewtonMaxSteps(10)
- else:
- solver.setNewtonMaxSteps(0)
+ solver.setNewtonMaxSteps(10)
+
 
  # add rdatas
  rdatas[equil_meth] = amici.runAmiciSimulation(model, solver, edata)
 
  # assert successful simulation
- assert rdatas[equil_meth]['status'] == amici.AMICI_SUCCESS
+ assert rdatas[equil_meth].status == amici.AMICI_SUCCESS
 
  # assert correct results
  for variable in ['llh', 'sllh', 'sx0', 'sx_ss', 'x_ss']:
@@ -374,4 +401,4 @@ def test_newton_solver_equilibration(preeq_fixture):
  solver.setNewtonMaxLinearSteps(10)
  rdata_spbcg = amici.runAmiciSimulation(model, solver, edata)
 
- assert rdata_spbcg['status'] == amici.AMICI_ERROR
+ assert rdata_spbcg.status == amici.AMICI_ERROR