Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Steady state sensitivity modes #1727

Merged
merged 24 commits into from
Apr 5, 2022
Merged

Steady state sensitivity modes #1727

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
4d37031
Steady State Sensitivity modes
plakrisenko Mar 15, 2022
ffec862
Merge branch 'develop' into stst_problem_settings
plakrisenko Mar 15, 2022
3e9df5b
integrationOnly and integrateIfNewtonFails in forward sensitivity case
plakrisenko Mar 21, 2022
9d31820
Merge branch 'develop' into stst_problem_settings
plakrisenko Mar 21, 2022
12eede5
turn off newton method in case integrationOnly and forward sensitivity
plakrisenko Mar 21, 2022
2709d51
python tests
plakrisenko Mar 21, 2022
e4d9222
fix
plakrisenko Mar 21, 2022
a515db3
Apply suggestions from code review
plakrisenko Mar 21, 2022
f974f7f
turn off newtin in findSteadyState
plakrisenko Mar 21, 2022
2c2a793
Merge branch 'stst_problem_settings' of github.com:AMICI-dev/AMICI in…
plakrisenko Mar 21, 2022
f5d5e08
fix test_compare_conservation_laws_sbml
plakrisenko Mar 23, 2022
df6b96b
Merge branch 'develop' into stst_problem_settings
plakrisenko Mar 25, 2022
af8e0d6
Merge branch 'develop' into stst_problem_settings
plakrisenko Mar 29, 2022
ee0ba61
fix
plakrisenko Mar 29, 2022
57f7fc4
fix
plakrisenko Mar 29, 2022
5f8dce1
take preequilibration into account
plakrisenko Mar 29, 2022
2c197a4
add extra parentheses
plakrisenko Mar 30, 2022
d32d64c
Merge branch 'develop' into stst_problem_settings
plakrisenko Apr 4, 2022
c6001d3
fix test_compare_conservation_laws_sbml
plakrisenko Apr 4, 2022
249f1a7
fix test_compare_conservation_laws_sbml
plakrisenko Apr 4, 2022
9fdcaf0
Merge branch 'develop' into stst_problem_settings
plakrisenko Apr 4, 2022
93e0431
Merge branch 'develop' into stst_problem_settings
plakrisenko Apr 4, 2022
b39b925
quick notebook fix
plakrisenko Apr 4, 2022
8a45940
Merge branch 'stst_problem_settings' of github.com:AMICI-dev/AMICI in…
plakrisenko Apr 4, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion include/amici/defines.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,8 @@ enum class NonlinearSolverIteration {
/** Sensitivity computation mode in steadyStateProblem */
enum class SteadyStateSensitivityMode {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we need some documentation of the behavior for the different values, for example in the function Model::setSteadyStateSensitivityMode(). https://github.com/AMICI-dev/AMICI/blob/master/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb also needs to be updated.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, this notebook definitely needs to be updated.

newtonOnly,
simulationFSA
integrationOnly,
integrateIfNewtonFails
};

/** State in which the steady state computation finished */
Expand Down
2 changes: 1 addition & 1 deletion python/tests/test_compare_conservation_laws_sbml.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ def test_compare_conservation_laws_sbml(models, edata_fixture):

# ----- compare simulations wo edata, sensi = 0, states and sensis -------
model_without_cl.setSteadyStateSensitivityMode(
amici.SteadyStateSensitivityMode.simulationFSA
amici.SteadyStateSensitivityMode.integrateIfNewtonFails
)

# run simulations
Expand Down
7 changes: 3 additions & 4 deletions python/tests/test_preequilibration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,8 @@ def test_equilibration_methods_with_adjoints(preeq_fixture):

rdatas = {}
equil_meths = [amici.SteadyStateSensitivityMode.newtonOnly,
amici.SteadyStateSensitivityMode.simulationFSA]
amici.SteadyStateSensitivityMode.integrationOnly,
amici.SteadyStateSensitivityMode.integrateIfNewtonFails]
sensi_meths = [amici.SensitivityMethod.forward,
amici.SensitivityMethod.adjoint]
settings = itertools.product(equil_meths, sensi_meths)
Expand Down Expand Up @@ -333,7 +334,7 @@ def test_newton_solver_equilibration(preeq_fixture):
edata.setObservedDataStdDev(np.hstack([stdy, stdy[0]]))

rdatas = {}
settings = [amici.SteadyStateSensitivityMode.simulationFSA,
settings = [amici.SteadyStateSensitivityMode.integrationOnly,
amici.SteadyStateSensitivityMode.newtonOnly]

for equil_meth in settings:
Expand All @@ -343,8 +344,6 @@ def test_newton_solver_equilibration(preeq_fixture):
model.setSteadyStateSensitivityMode(equil_meth)
if equil_meth == amici.SteadyStateSensitivityMode.newtonOnly:
solver.setNewtonMaxSteps(10)
else:
solver.setNewtonMaxSteps(0)

# add rdatas
rdatas[equil_meth] = amici.runAmiciSimulation(model, solver, edata)
Expand Down
4 changes: 2 additions & 2 deletions python/tests/test_pysb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ def get_data(model):
solver = model.getSolver()
model.setTimepoints(np.linspace(0, 60, 61))
model.setSteadyStateSensitivityMode(
amici.SteadyStateSensitivityMode.simulationFSA
amici.SteadyStateSensitivityMode.integrateIfNewtonFails
)

rdata = amici.runAmiciSimulation(model, solver)
Expand All @@ -220,7 +220,7 @@ def get_results(model, edata):
edata.reinitializeFixedParameterInitialStates = True
model.setTimepoints(np.linspace(0, 60, 61))
model.setSteadyStateSensitivityMode(
amici.SteadyStateSensitivityMode.simulationFSA
amici.SteadyStateSensitivityMode.integrateIfNewtonFails
)
return amici.runAmiciSimulation(model, solver, edata)

Expand Down
3 changes: 1 addition & 2 deletions python/tests/test_sbml_import.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,10 +205,9 @@ def test_presimulation(sbml_example_presimulation_module):
"""Test 'presimulation' test model"""
model = sbml_example_presimulation_module.getModel()
solver = model.getSolver()
solver.setNewtonMaxSteps(0)
model.setTimepoints(np.linspace(0, 60, 61))
model.setSteadyStateSensitivityMode(
amici.SteadyStateSensitivityMode.simulationFSA
amici.SteadyStateSensitivityMode.integrationOnly
)
solver.setSensitivityOrder(amici.SensitivityOrder.first)
model.setReinitializeFixedParameterInitialStates(True)
Expand Down
39 changes: 28 additions & 11 deletions src/steadystateproblem.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,16 +90,21 @@ void SteadystateProblem::workSteadyStateBackwardProblem(
void SteadystateProblem::findSteadyState(const Solver *solver, Model *model,
int it) {
steady_state_status_.resize(3, SteadyStateStatus::not_run);
bool turnOffNewton = solver->getSensitivityMethod() ==
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't we generally wan't to turn off the newton solver if model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly? Why also require forward sensitivities here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because in adjoint case steady-state computation part is independent of sensitivities computation.
If model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly in adjoint case, I think model steady-state could still be found by Newton method.

Copy link
Member

@FFroehlich FFroehlich Mar 24, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, think I am starting to understand. But that ignores solver->getSensitivityMethodPreequilibration(), right?
I think the best way is to go via SteadystateProblem::getSensitivityFlag here.

SensitivityMethod::forward && model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::integrationOnly;


/* First, try to run the Newton solver */
findSteadyStateByNewtonsMethod(model, false);
if (!turnOffNewton)
findSteadyStateByNewtonsMethod(model, false);

/* Newton solver didn't work, so try to simulate to steady state */
if (!checkSteadyStateSuccess())
findSteadyStateBySimulation(solver, model, it);

/* Simulation didn't work, retry the Newton solver from last sim state. */
if (!checkSteadyStateSuccess())
if (!turnOffNewton && !checkSteadyStateSuccess())
findSteadyStateByNewtonsMethod(model, true);

/* Nothing worked, throw an as informative error as possible */
Expand Down Expand Up @@ -237,11 +242,17 @@ void SteadystateProblem::computeSteadyStateQuadrature(const Solver *solver,
We therefore compute the integral over xB first and then do a
matrix-vector multiplication */

/* Try to compute the analytical solution for quadrature algebraically */
getQuadratureByLinSolve(model);
auto sensitivityMode = model->getSteadyStateSensitivityMode();

/* Analytical solution didn't work, perform simulation instead */
if (!hasQuadrature())
/* Try to compute the analytical solution for quadrature algebraically */
if (sensitivityMode == SteadyStateSensitivityMode::newtonOnly
|| sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails)
getQuadratureByLinSolve(model);

/* Perform simulation */
if (sensitivityMode == SteadyStateSensitivityMode::integrationOnly ||
(sensitivityMode == SteadyStateSensitivityMode::integrateIfNewtonFails
&& !hasQuadrature()))
getQuadratureBySimulation(solver, model);

/* If analytic solution and integration did not work, throw an Exception */
Expand Down Expand Up @@ -360,8 +371,10 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model,
bool forwardSensisAlreadyComputed =
solver->getSensitivityOrder() >= SensitivityOrder::first &&
steady_state_status_[1] == SteadyStateStatus::success &&
model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::simulationFSA;
(model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::integrationOnly ||
model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::integrateIfNewtonFails);

bool simulationStartedInSteadystate =
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this code requires some updates as the newton solver can't check if we started in a steadystate if it is turned off.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When is this an issue?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For models that start in steady-state (e.g., computed outside of amici or extracted from other simulations), we don't want to be forced to deactivate the newton solver. See #1422.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But thinking about it more carefully, we will only have misleading values if the newton solver is deactivated. In that case we will simply rextract sensitvities from the solver, which should be fine (but we probably perform a single, unnecessary simulaton step with sensitivities, which might be expensive).

steady_state_status_[0] == SteadyStateStatus::success &&
Expand All @@ -386,9 +399,13 @@ bool SteadystateProblem::getSensitivityFlag(const Model *model,
!simulationStartedInSteadystate;

/* When we're creating a new solver object */
bool needForwardSensiAtCreation =
needForwardSensisPreeq && model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::simulationFSA;
bool needForwardSensiAtCreation =
needForwardSensisPreeq &&
(model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::integrationOnly ||
model->getSteadyStateSensitivityMode() ==
SteadyStateSensitivityMode::integrateIfNewtonFails
);

/* Check if we need to store sensis */
switch (context) {
Expand Down
4 changes: 2 additions & 2 deletions tests/petab_test_suite/test_petab_suite.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
from amici.petab_import import import_petab_problem, PysbPetabProblem
from amici.petab_objective import (
simulate_petab, rdatas_to_measurement_df, create_parameterized_edatas)
from amici import SteadyStateSensitivityMode_simulationFSA
from amici import SteadyStateSensitivityMode

logger = get_logger(__name__, logging.DEBUG)
set_log_level(get_logger("amici.petab_import"), logging.DEBUG)
Expand Down Expand Up @@ -134,7 +134,7 @@ def check_derivatives(problem: petab.Problem, model: amici.Model) -> None:
# Required for case 9 to not fail in
# amici::NewtonSolver::computeNewtonSensis
model.setSteadyStateSensitivityMode(
SteadyStateSensitivityMode_simulationFSA)
SteadyStateSensitivityMode.integrateIfNewtonFails)

for edata in create_parameterized_edatas(
amici_model=model, petab_problem=problem,
Expand Down