adjoints + steadystate #1099

paulstapor · 2020-05-27T06:24:50Z

Just want to see whether checks pass.
PR for this one will be done by @paszkow from his fork, after confirming his scripts still run and results are correct.

Let both methods have all parameters required by model specified implementation similarily as it is done for fJ and fJSparse methods.

The NewtonSolver contains a prepareLinearSystemB method that retrieves and passes the Jacobian JB. This allows solving linear systems JB*x = rhs.

When a linear system is prepared construct a Jacobian and use it directly in matrix-vector multiplication instead of mutliple calls to Model::fJv.

Avoid backwards integration to compute sensitivities when handled with a steady-state solution. It is enough to solve a linear system.

…ities

update testfunctions header

…dystateproblem

…ICI into paszkow-adjoint_sensitivities

paulstapor · 2020-06-18T12:45:55Z

So, from my side, this should be mergable now.
Doing some last checks on this with @paszkow , but beyond that, the code should be fine...

paulstapor · 2020-06-18T13:40:37Z

include/amici/rdata.h

@@ -225,24 +225,27 @@ class ReturnData {
    /** employed order forward problem (dimension: nt) */
    std::vector<int> order;

-    /** computation time of forward solve [s] */
+    /** computation time of forward solve [ms] */


We used to have time reporting not uniformly. At least for one of them, the unit was incorrect. Now, everything should actually be in milliseconds.
No preference, whether seconds or milliseconds are better. We just need to have it the same way for all cpu times...

preferring seconds. but either is fine, if coherently used. good to include it already in the variable/function name.

@paszkow proposed to use the chorno library (https://en.cppreference.com/w/cpp/chrono), which looks nicer and may be more precise than what we currently have...
Maybe something for the future?

FFroehlich

Looks good 👍

FFroehlich · 2020-06-18T13:56:17Z

python/tests/test_preequilibration.py

-        assert rdatas[sensi_meth]['status'] == amici.AMICI_SUCCESS
+    settings = []
+    for equil_meth in [amici.SteadyStateSensitivityMode.newtonOnly,]:
+#                       amici.SteadyStateSensitivityMode.simulationFSA,]:


Suggested change

# amici.SteadyStateSensitivityMode.simulationFSA,]:

amici.SteadyStateSensitivityMode.simulationFSA,]:

FFroehlich · 2020-06-18T13:56:47Z

python/tests/test_preequilibration.py

+    settings = []
+    for equil_meth in [amici.SteadyStateSensitivityMode.newtonOnly,]:
+#                       amici.SteadyStateSensitivityMode.simulationFSA,]:
+        for sensi_meth in [amici.SensitivityMethod.forward,


use itertools.product to reduce nesting?

FFroehlich · 2020-06-18T14:02:27Z

src/backwardproblem.cpp

+    for (int jt = model->nt() - 1; jt >= 0; jt--)
+        if (std::isinf(model->getTimepoint(it)))
+            --it;
+    if (it < model->nt() - 1){


cant we just initialize xQB to zeros for it < model->nt() - 1 and pass it to initializeB? I think that would be better readable. Also shouldn't if/else cases be switched here?

The problem ist that, if we have postequilibration (which is the case for it < model->nt() - 1), xQB already contains the quadratures from backward postequilibration. Hence, in that case, we don't want to overwrite them, but still want to reinitialize xB and dxB.
I have another idea, which should simplify this. Will commit.

FFroehlich · 2020-06-18T14:04:29Z

src/backwardproblem.cpp

-        /* reinit states */
-        solver->reInitB(which, t, xB, dxB);
-        solver->quadReInitB(which, xQB);
+        /* we still need to integrate from first datapoint to tstart */


I believe this code should be outside the loop, we only want to run this once.

True. Done.

FFroehlich · 2020-06-18T14:05:49Z

src/newton_solver.cpp

+    // Get the diagonal and ensure negativity of entries is ns_J. Note that diag(JB) = -diag(J).
+    model->fJDiag(*t, ns_Jdiag, 0.0, *x, dx);
+
+    ns_p.set(1.0);


duplicate code, move to seperate function?

yikes. This code will be removed soon anyway, current issue, see #1073.
This is really an old handcrafted reimplementation of BICGSTAB. Prefer not to do changes now, but to fully remove and correctly reimplement the remaining linear solvers (i.e., all except DENSE and KLU).

FFroehlich · 2020-06-18T14:06:35Z

src/newton_solver.cpp

@@ -271,7 +325,7 @@ void NewtonSolverIterative::linsolveSPBCG(int ntry, int nnewt,
    double alpha = 1.0;

    // can be set to 0 at the moment


I believe this comment is deprecated

dweindl · 2020-06-18T14:12:55Z

CMakeLists.txt

@@ -1,7 +1,7 @@
 #
 # Build AMICI library
 #
-cmake_minimum_required(VERSION 3.3)
+cmake_minimum_required(VERSION 3.12)


What was that for if there is no other change in CMake files?

@paszkow pointed out that using the option --parallel requires at least 3.12., this is why I just put it...
Can also remove for the moment, if this comes with issues...

Hmm... that would have to be checked in the build script. Here it is too late.

okay, I see. ... yes, okay. I'll remove for the moment. Let's do that later and not mix up too many things.

dweindl · 2020-06-18T14:21:25Z

src/backwardproblem.cpp

-    solver->setupB(&which, model->getTimepoint(it), model, xB, dxB, xQB);
+    /* If we have posteq, infty timepoints were already treated */
+    for (int jt = model->nt() - 1; jt >= 0; jt--)
+        if (std::isinf(model->getTimepoint(it)))


Suggested change

if (std::isinf(model->getTimepoint(it)))

if (std::isinf(model->getTimepoint(jt)))

?

Otherwise this loop does not make much sense, does it?

arg, true...!

dweindl · 2020-06-18T14:22:28Z

src/backwardproblem.cpp


-    --it;
+    if ((it >= 0 || discs.size() > 0) && model->getTimepoint(it) > model->t0())


Suggested change

if ((it >= 0 || discs.size() > 0) && model->getTimepoint(it) > model->t0())

if ((it >= 0 || !discs.empty()) && model->getTimepoint(it) > model->t0())

dweindl · 2020-06-18T14:24:11Z

src/newton_solver.cpp

@@ -12,11 +12,13 @@
 #include <cstring>
 #include <ctime>
 #include <cmath>
+#include <iostream>


Unused?! Leftover from debugging?

True, aye... Looks old...

FFroehlich · 2020-06-18T14:31:38Z

are JB and JSparseB generated in both MATLAB and Python interfaces?

paulstapor · 2020-06-18T18:12:02Z

are JB and JSparseB generated in both MATLAB and Python interfaces?

Yes, also Matlab generates JB and JSparseB

paulstapor · 2020-06-19T11:50:34Z

If tests pass now, would be strongly in favor of merging now...
If anyone present, for any reason knows why these changes should not be merged, let him speak now or forever hold his peace.

sonarqubecloud · 2020-06-19T12:02:14Z

SonarCloud Quality Gate failed.

0 Bugs
0 Vulnerabilities (and 0 Security Hotspots to review)
3 Code Smells

27.3% Coverage
4.7% Duplication

Łukasz Paszkowski and others added 30 commits January 19, 2020 14:12

Extract methods fJB and fJSparseB to AbstractModel.

c0f8e4f

Let both methods have all parameters required by model specified implementation similarily as it is done for fJ and fJSparse methods.

Prepare linear system for the JB matrix.

5f50a3b

The NewtonSolver contains a prepareLinearSystemB method that retrieves and passes the Jacobian JB. This allows solving linear systems JB*x = rhs.

Remove fJv calls from NewtonSolverIterative.

bb1f4b9

When a linear system is prepared construct a Jacobian and use it directly in matrix-vector multiplication instead of mutliple calls to Model::fJv.

Enable iterative Newton solver in solving JB*x=rhs.

c761883

Simplify computing adjoint sensitivities for steady-states.

8e31038

Avoid backwards integration to compute sensitivities when handled with a steady-state solution. It is enough to solve a linear system.

Merge remote-tracking branch 'upstream/develop' into adjoint_sensitiv…

58965cb

…ities

Merge branch 'develop' into adjoint_sensitivities

99b2275

update testfunctions header

24c4e24

Merge branch 'develop' into adjoint_sensitivities

994b4bf

Merge pull request #1 from ICB-DCM/paszkow-adjoint_sensitivities

bcbb501

update testfunctions header

Rename a function to enable Newton solver in addjoint problem.

3feaa55

Simplify computing the inner product.

29512c8

Allow steady-state and time itegration in adjoint problem.

17f56e8

Move Neton solver computations to SteadystateProblem.

c661beb

Merge branch 'develop' into adjoint_sensitivities

94da2fa

Fix dimension mismatch for python generated models

8a6e04e

Merge branch 'develop' into adjoint_sensitivities

0cee19b

Merge branch 'develop' into adjoint_sensitivities

777f58d

Merge branch 'develop' into adjoint_sensitivities

492a7a9

Measure computation time in backward problem using Newton solver.

f8fbf3f

try to resolve merge conflicts

06fcfa4

resolve errors

7238329

add unit test

83dbfe6

fix steadystate sensitivity computation

ddb23de

add performance test for new equilibration sensi functionality

020b53d

trying to fix

1c7d6de

fixed bugs for the moment

0ec02b8

fixing newton solver with adjoints

d217611

cleanup

4c35680

enforce minimum number for numlinsteps field in newtonsolver and stea…

1fe8132

…dystateproblem

paulstapor added 6 commits June 17, 2020 11:13

improve doc for postequilibration backward

7c0ebd8

Merge branch 'paszkow-adjoint_sensitivities' of github.com:ICB-DCM/AM…

941ab81

…ICI into paszkow-adjoint_sensitivities

fix unit in cpu time

4b4be0e

fix unit in cpu time

fc8d172

add posteq_cuptimeB

66f6b9b

resolve

55a261a

paulstapor marked this pull request as ready for review June 18, 2020 12:05

paulstapor added 2 commits June 18, 2020 14:09

Merge branch 'develop' into paszkow-adjoint_sensitivities

07aff30

Merge branch 'paszkow-adjoint_sensitivities' of github.com:ICB-DCM/AM…

1ca105a

…ICI into paszkow-adjoint_sensitivities

paulstapor requested a review from FFroehlich June 18, 2020 12:46

paulstapor commented Jun 18, 2020

View reviewed changes

Merge branch 'develop' into paszkow-adjoint_sensitivities

f1e6402

FFroehlich approved these changes Jun 18, 2020

View reviewed changes

dweindl reviewed Jun 18, 2020

View reviewed changes

paulstapor added 3 commits June 19, 2020 11:47

address code review

1eb201b

Merge branch 'develop' into paszkow-adjoint_sensitivities

17383dd

address sonarcloud

6cfa0af

paulstapor merged commit 6ee291f into develop Jun 19, 2020

This was referenced Jun 21, 2020

Documentation of Steady-state solver missing #1072

Closed

Implement adjoint steady-state sensitivities & exploit sparsity in steady-state sensitivity computation #876

Closed

paulstapor deleted the paszkow-adjoint_sensitivities branch July 6, 2020 14:39

FFroehlich mentioned this pull request Jul 16, 2020

Release 0.11.2 #1159

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

adjoints + steadystate #1099

adjoints + steadystate #1099

paulstapor commented May 27, 2020

paulstapor commented Jun 18, 2020

paulstapor Jun 18, 2020

dweindl Jun 18, 2020

paulstapor Jun 18, 2020

FFroehlich left a comment

FFroehlich Jun 18, 2020

FFroehlich Jun 18, 2020

FFroehlich Jun 18, 2020

paulstapor Jun 19, 2020

FFroehlich Jun 18, 2020

paulstapor Jun 19, 2020

FFroehlich Jun 18, 2020

paulstapor Jun 18, 2020

FFroehlich Jun 18, 2020

dweindl Jun 18, 2020

paulstapor Jun 18, 2020 •

edited

Loading

dweindl Jun 18, 2020

paulstapor Jun 18, 2020

dweindl Jun 18, 2020

paulstapor Jun 18, 2020

dweindl Jun 18, 2020

dweindl Jun 18, 2020

paulstapor Jun 18, 2020

FFroehlich commented Jun 18, 2020

paulstapor commented Jun 18, 2020

paulstapor commented Jun 19, 2020

sonarqubecloud bot commented Jun 19, 2020

	# amici.SteadyStateSensitivityMode.simulationFSA,]:
	amici.SteadyStateSensitivityMode.simulationFSA,]:

		@@ -271,7 +325,7 @@ void NewtonSolverIterative::linsolveSPBCG(int ntry, int nnewt,
		double alpha = 1.0;

		// can be set to 0 at the moment

	if (std::isinf(model->getTimepoint(it)))
	if (std::isinf(model->getTimepoint(jt)))


		--it;
		if ((it >= 0 \|\| discs.size() > 0) && model->getTimepoint(it) > model->t0())

	if ((it >= 0 \|\| discs.size() > 0) && model->getTimepoint(it) > model->t0())
	if ((it >= 0 \|\| !discs.empty()) && model->getTimepoint(it) > model->t0())

adjoints + steadystate #1099

adjoints + steadystate #1099

Conversation

paulstapor commented May 27, 2020

paulstapor commented Jun 18, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FFroehlich left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

paulstapor Jun 18, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FFroehlich commented Jun 18, 2020

paulstapor commented Jun 18, 2020

paulstapor commented Jun 19, 2020

sonarqubecloud bot commented Jun 19, 2020

paulstapor Jun 18, 2020 •

edited

Loading