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.

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

merged 24 commits into from
Apr 5, 2022

Conversation

plakrisenko
Copy link
Member

@plakrisenko plakrisenko commented Mar 16, 2022

Settings for specifying sensitivities computation method at steady state have been implemented/reworked.
Three SteadyStateSensitivityMode are available for both SensitivityMethods (forward and adjoint):

  • newtonOnly: only Newton's method is used for sensitivities computation
  • integrationOnly: only integration is used for sensitivities computation, which in forward sensitivity analysis case means that numerical integration has to be used to find the steady state
  • integrateIfNewtonFails: more flexible option, where simulation is used if Newton's method fails

Note: Previously available SteadyStateSensitivityMode simulationFSA has been removed.

@codecov
Copy link

codecov bot commented Mar 21, 2022

Codecov Report

Merging #1727 (8a45940) into develop (742c464) will increase coverage by 0.06%.
The diff coverage is 100.00%.

Impacted file tree graph

@@             Coverage Diff             @@
##           develop    #1727      +/-   ##
===========================================
+ Coverage    77.82%   77.89%   +0.06%     
===========================================
  Files           74       74              
  Lines        12018    12033      +15     
===========================================
+ Hits          9353     9373      +20     
+ Misses        2665     2660       -5     
Flag Coverage Δ
cpp 75.05% <100.00%> (+0.11%) ⬆️
petab 59.18% <ø> (ø)
python 68.24% <ø> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
src/steadystateproblem.cpp 85.57% <100.00%> (+1.03%) ⬆️
src/solver_cvodes.cpp 70.43% <0.00%> (+0.16%) ⬆️
src/exception.cpp 81.08% <0.00%> (+5.40%) ⬆️

@plakrisenko plakrisenko marked this pull request as ready for review March 21, 2022 16:46
tests/petab_test_suite/test_petab_suite.py Outdated Show resolved Hide resolved
tests/petab_test_suite/test_petab_suite.py Outdated Show resolved Hide resolved
src/steadystateproblem.cpp Outdated Show resolved Hide resolved
src/steadystateproblem.cpp Outdated Show resolved Hide resolved
@FFroehlich
Copy link
Member

I would guess that the expected preeq_status and posteq_status fields in some of the python test will have to be updated.

@plakrisenko
Copy link
Member Author

I would guess that the expected preeq_status and posteq_status fields in some of the python test will have to be updated.

I should have set different setSteadyStateSensitivityMode in that test

Copy link
Member

@FFroehlich FFroehlich left a comment

Choose a reason for hiding this comment

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

great work 👍 , not sure why the last remaining test fails, don't see anything that could lead to test failure.

@@ -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

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.

Copy link
Member

@FFroehlich FFroehlich left a comment

Choose a reason for hiding this comment

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

Sorry, noticed a couple more things that need to be adressed :(

(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).

@@ -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:/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.

@@ -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

@FFroehlich FFroehlich Mar 24, 2022

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.

@plakrisenko
Copy link
Member Author

Sorry, noticed a couple more things that need to be adressed :(

You are right about the preequilibration.
I guess the logic should be rather
bool turnOffNewton = model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly && (it == -1 && solver->getSensitivityMethodPreequilibration() == SensitivityMethod::forward || solver->getSensitivityMethod() == SensitivityMethod::forward);
Right?

@FFroehlich
Copy link
Member

Sorry, noticed a couple more things that need to be adressed :(

You are right about the preequilibration. I guess the logic should be rather bool turnOffNewton = model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly && (it == -1 && solver->getSensitivityMethodPreequilibration() == SensitivityMethod::forward || solver->getSensitivityMethod() == SensitivityMethod::forward); Right?

Sorry, noticed a couple more things that need to be adressed :(

You are right about the preequilibration. I guess the logic should be rather bool turnOffNewton = model->getSteadyStateSensitivityMode() == SteadyStateSensitivityMode::integrationOnly && (it == -1 && solver->getSensitivityMethodPreequilibration() == SensitivityMethod::forward || solver->getSensitivityMethod() == SensitivityMethod::forward); Right?

Yes, that looks correct.

Copy link
Member

@FFroehlich FFroehlich left a comment

Choose a reason for hiding this comment

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

👍

@sonarcloud
Copy link

sonarcloud bot commented Apr 5, 2022

SonarCloud Quality Gate failed.    Quality Gate failed

Bug A 0 Bugs
Vulnerability A 0 Vulnerabilities
Security Hotspot A 0 Security Hotspots
Code Smell A 0 Code Smells

0.0% 0.0% Coverage
0.0% 0.0% Duplication

@plakrisenko plakrisenko merged commit 707ec6c into develop Apr 5, 2022
@FFroehlich FFroehlich deleted the stst_problem_settings branch April 5, 2022 12:08
dweindl pushed a commit that referenced this pull request Apr 5, 2022
Settings for specifying sensitivities computation method at steady state have been implemented/reworked.
Three `SteadyStateSensitivityMode` are available for both `SensitivityMethod`s (`forward` and `adjoint`):
- `newtonOnly`: only Newton's method is used for sensitivities computation
- `integrationOnly`: only integration is used for sensitivities computation, which in forward sensitivity analysis case means that numerical integration has to be used to find the steady state
- `integrateIfNewtonFails`: more flexible option, where simulation is used if Newton's method fails

Note: Previously available `SteadyStateSensitivityMode` `simulationFSA` has been removed.

Co-authored-by: Fabian Fröhlich <[email protected]>
FFroehlich pushed a commit that referenced this pull request Apr 5, 2022
* exploit stoichiometric matrix

* add solver mapping

* fixup

* fixups

* fixup

* fixup conservation laws

* Steady state sensitivity modes (#1727)

* Steady State Sensitivity modes

* integrationOnly and integrateIfNewtonFails in forward sensitivity case

* turn off newton method in case integrationOnly and forward sensitivity

* python tests

* Apply suggestions from code review

Co-authored-by: Fabian Fröhlich <[email protected]>

* turn off newtin in findSteadyState

* fix test_compare_conservation_laws_sbml

* take preequilibration into account

* fix test_compare_conservation_laws_sbml

* quick notebook fix

Co-authored-by: Fabian Fröhlich <[email protected]>

* update according to reviews

Co-authored-by: Polina Lakrisenko <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add SteadyStateSensitivityMode(s) for Adjoint sensitivities computation in steadyStateProblem
2 participants