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

Enable Newton solver with simulation #1075

Merged
merged 10 commits into from
May 9, 2020
Merged

Enable Newton solver with simulation #1075

merged 10 commits into from
May 9, 2020

Conversation

paulstapor
Copy link
Contributor

This is a suggestion how to fix an existing problem in the steadystate solver (happy about alternative suggestions how to handle it):
Currently, if a model has the timepoint inf specified and forward sensis and Newton's method fails due to a poor initial guess, the steadystate solver will integrate a long time with full forward sensis.
However, this is not ideal (or even completely prohibitive) some applications:
If the Jacobian is invertible, it would be fine to simulate without sensis and then just solve the linear system of equations. This can be way faster...

One point I'm not happy about in this solution: I had to implement a function in the solver code to switch off forward sensi computation. As this one is to be used in steadystateproblem, it must be part of the public API of solver. However, this function touches stuff in the solver object, that should imho not be touchable by a public API...
Hence, I would be happy about suggestions how to solve this better...

@paulstapor paulstapor changed the title Fix newtonsolver Enable Newton solver with simulation May 6, 2020
@codecov
Copy link

codecov bot commented May 6, 2020

Codecov Report

Merging #1075 into develop will decrease coverage by 0.31%.
The diff coverage is 82.19%.

Impacted file tree graph

@@             Coverage Diff             @@
##           develop    #1075      +/-   ##
===========================================
- Coverage    72.55%   72.24%   -0.32%     
===========================================
  Files           59       59              
  Lines         9180     9216      +36     
===========================================
- Hits          6661     6658       -3     
- Misses        2519     2558      +39     
Flag Coverage Δ
#cpp 73.70% <82.19%> (-0.48%) ⬇️
#python 68.73% <ø> (+0.03%) ⬆️
Impacted Files Coverage Δ
include/amici/backwardproblem.h 100.00% <ø> (ø)
include/amici/solver.h 30.76% <ø> (ø)
include/amici/solver_cvodes.h 100.00% <ø> (ø)
include/amici/solver_idas.h 100.00% <ø> (ø)
src/backwardproblem.cpp 58.33% <40.00%> (-6.23%) ⬇️
src/solver_idas.cpp 41.82% <65.00%> (-0.07%) ⬇️
src/forwardproblem.cpp 93.02% <80.00%> (-1.03%) ⬇️
src/solver.cpp 75.28% <92.30%> (-0.04%) ⬇️
src/solver_cvodes.cpp 72.02% <94.44%> (-1.37%) ⬇️
src/steadystateproblem.cpp 79.66% <100.00%> (-2.59%) ⬇️
... and 11 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 5992f34...1e7a1e4. Read the comment docs.

@paulstapor
Copy link
Contributor Author

Okay, it is 100% reasonable, that the failing test does indeed fail. I'm just confused about the fact that it worked beforehand... 🤔 😕
Using the steady-state solver with adjoints outside preequilibration is not supposed to have a well-defined behavior (yet, will change with #917) and should trigger some raise, actually.
It's interesting that we have a test on this. It is substantially more stunning, that the test passed beforehand... 😂

I will add a jupyter notebook as documentation of the steadystate solver within this PR, which explains what the steadystate solver should do... (in my opinion). Happy to excahnge opinions on that...

@FFroehlich
Copy link
Member

This is a suggestion how to fix an existing problem in the steadystate solver (happy about alternative suggestions how to handle it):
Currently, if a model has the timepoint inf specified and forward sensis and Newton's method fails due to a poor initial guess, the steadystate solver will integrate a long time with full forward sensis.
However, this is not ideal (or even completely prohibitive) some applications:
If the Jacobian is invertible, it would be fine to simulate without sensis and then just solve the linear system of equations. This can be way faster...

One point I'm not happy about in this solution: I had to implement a function in the solver code to switch off forward sensi computation. As this one is to be used in steadystateproblem, it must be part of the public API of solver. However, this function touches stuff in the solver object, that should imho not be touchable by a public API...
Hence, I would be happy about suggestions how to solve this better...

Thats part of the public API of CVODES/IDAS, so I don't see any reason why we shouldn't expose this functionality in AMICI. The implementation absolutely fine from my point of view.

@FFroehlich
Copy link
Member

Okay, it is 100% reasonable, that the failing test does indeed fail. I'm just confused about the fact that it worked beforehand... 🤔 😕
Using the steady-state solver with adjoints outside preequilibration is not supposed to have a well-defined behavior (yet, will change with #917) and should trigger some raise, actually.
It's interesting that we have a test on this. It is substantially more stunning, that the test passed beforehand... 😂

I will add a jupyter notebook as documentation of the steadystate solver within this PR, which explains what the steadystate solver should do... (in my opinion). Happy to excahnge opinions on that...

True, we never actually check whether sllh has any meaningful content for that test, so this might have passed just by chance. Fine with me to temporarily disable the amici.SensitivityMethod_adjoint case and readd in #917.

@FFroehlich
Copy link
Member

Please have a look at my changes and if tests pass and everything looks fine, feel free to merge!

# add rdatas
rdatas[sensi_meth] = amici.runAmiciSimulation(model, solver, edata)
assert rdatas[sensi_meth]['status'] == amici.AMICI_SUCCESS

for variable in ['llh', 'sllh']:
Copy link
Contributor Author

Choose a reason for hiding this comment

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

I didn't add it in my commits, but checked it rather by hand:
In an ideal world, also sx0 and sx_ss are checked for here...

# add rdatas
rdatas[sensi_meth] = amici.runAmiciSimulation(model, solver, edata)
assert rdatas[sensi_meth]['status'] == amici.AMICI_SUCCESS

for variable in ['llh', 'sllh']:
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Suggested change
for variable in ['llh', 'sllh']:
for variable in ['llh', 'sllh', 'sx_ss', 'sx0']:


rdatas = {}
for sensi_meth in [amici.SensitivityMethod.forward,
amici.SensitivityMethod.adjoint]:
Copy link
Contributor Author

Choose a reason for hiding this comment

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

That's just a question, not necessary, but:
Do we also want to not only check FSA and ASA, but also FD?

solver.setSensitivityMethod(sensi_meth)
solver.setNewtonMaxSteps(0)
solver.SteadyStateSensitivityMethod = \
amici.SteadyStateSensitivityMode.simulationFSA
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Add another loop for Newton solver as well?
Or are we overtesting then?

@sonarcloud
Copy link

sonarcloud bot commented May 8, 2020

SonarCloud Quality Gate failed.

Bug A 0 Bugs
Vulnerability A 0 Vulnerabilities (and Security Hotspot 0 Security Hotspots to review)
Code Smell A 1 Code Smell

54.3% 54.3% Coverage
0.0% 0.0% Duplication

@paulstapor paulstapor merged commit 3ea1e29 into develop May 9, 2020
@paulstapor paulstapor deleted the fix_newtonsolver branch May 9, 2020 09:56
@FFroehlich FFroehlich mentioned this pull request May 10, 2020
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.

2 participants