State sensitivities computation at steady state

[plakrisenko]

There are two approaches available in AMICI for state sensitivities computation at steady state: numerical integration and solving a linear system (computeNewtonSensis). There are now also three SteadyStateSensitivityModes available: newtonOnly, integrationOnly, integrateIfNewtonFails.

I guess, originally SteadyStateSensitivityMode was only meant for forward sensitivity analysis. In this pull-request #1727 its usage was extended to adjoint sensitivity analysis.

However, integrateIfNewtonFails is not quite compatible with FSA at the moment:

• in ASA case it's used here https://github.com/AMICI-dev/AMICI/blob/master/src/steadystateproblem.cpp#L281
  With this mode getQuadratureByLinSolve is attempted first and then getQuadratureBySimulation.
• for FSA the order seems to be inverted:
  https://github.com/AMICI-dev/AMICI/blob/master/src/steadystateproblem.cpp#L70-L78
  sensitivities computation by integration is part of findSteadyState if on. State sensitivities computation at steady state by solving a linear system is done in computeNewtonSensis.

One could maybe move computeNewtonSensis to here https://github.com/AMICI-dev/AMICI/blob/master/src/steadystateproblem.cpp#L121 (and after the second Newton's method attempt). Then the order will be similar to ASA case.

Or should modes for sensitivities computation at steady state be separated for FSA and ASA, is it too confusing now?

[FFroehlich]

I agree, it's all a big giant mess.

There are two things that we need to control (i) how the steady state is computed (integration/newton) and (ii) how the sensitivities are computed (integration/newton). As far as I can tell, only newton steadystate + integration sensi does not work (not sure about this for ASA), but everything else should be viable. In the current implementation (i) and (ii) can't be controlled independently by the user, but have to be set via SteadyStateSensitivityMode and SensitivityMethodPreequilibration, which adds an additional layer of complexity between pre- and post-equilibration. This leads to very complex checks in

AMICI/src/steadystateproblem.cpp

Line 112 in 81d32e7

bool turnOffNewton = solver.getNewtonMaxSteps() == 0 || (

and

AMICI/src/steadystateproblem.cpp

Line 401 in 81d32e7

bool SteadystateProblem::getSensitivityFlag(

(which, to make things worse, does different things depending on which context it is called in).

I believe the reason you think the implementation is inconsistent is because in one case you are looking at controlling (i) vs (ii), but I think the issue with the code is really more fundamental. Any effort to simplify this code is very welcome, I made some attempts at refactoring, but ultimately gave up as Paul still had intentions to fix this himself.

[dweindl]

I made some attempts at refactoring, but ultimately gave up

i think that was that: #1485
maybe something to be learned from the discussion back then...

[plakrisenko]

In the current implementation (i) and (ii) can't be controlled independently by the user, but have to be set via SteadyStateSensitivityMode and SensitivityMethodPreequilibration

But SteadyStateSensitivityMode and SensitivityMethodPreequilibration don't control (i) at all?