Adjoint sensitivities, fix #876

include/amici/abstract_model.h

@@ -82,6 +82,22 @@ class AbstractModel {
  const AmiVector &dx, const AmiVector &xdot,
  SUNMatrix J) = 0;
 
+ /**
+ * @brief Dense Jacobian function
+ * @param t time
+ * @param cj scaling factor (inverse of timestep, DAE only)
+ * @param x state
+ * @param dx time derivative of state (DAE only)
+ * @param xB Vector with the adjoint states
+ * @param dxB Vector with the adjoint derivative states
+ * @param xBdot Vector with the adjoint right hand side (unused)
+ * @param JB dense matrix to which values of the jacobian will be written
+ */
+ virtual void fJB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB,
+ const AmiVector &dxB, const AmiVector &xBdot,
+ SUNMatrix JB) = 0;
+
  /**
  * @brief Sparse Jacobian function
  * @param t time
@@ -95,6 +111,22 @@ class AbstractModel {
  const AmiVector &x, const AmiVector &dx,
  const AmiVector &xdot, SUNMatrix J) = 0;
 
+ /**
+ * @brief Dense Jacobian function
+ * @param t time
+ * @param cj scaling factor (inverse of timestep, DAE only)
+ * @param x state
+ * @param dx time derivative of state (DAE only)
+ * @param xB Vector with the adjoint states
+ * @param dxB Vector with the adjoint derivative states
+ * @param xBdot Vector with the adjoint right hand side (unused)
+ * @param JB dense matrix to which values of the jacobian will be written
+ */
+ virtual void fJSparseB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB,
+ const AmiVector &dxB, const AmiVector &xBdot,
+ SUNMatrix JB) = 0;
+
  /**
  * @brief Diagonal Jacobian function
  * @param t time
@@ -669,13 +701,13 @@ class AbstractModel {
  const realtype *p, const realtype *k, const realtype *h,
  const realtype *w, const realtype *tcl,
  const realtype *stcl);
- 
+
  /**
  * @brief Model specific implementation for dwdp, column pointers
  * @param indexptrs column pointers
  **/
  virtual void fdwdp_colptrs(sunindextype *indexptrs);
- 
+
  /**
  * @brief Model specific implementation for dwdp, row values
  * @param indexvals row values

include/amici/backwardproblem.h

@@ -101,6 +101,13 @@ class BackwardProblem {
  */
  realtype getTnext(std::vector<realtype> const& troot, int iroot, int it);
 
+ /**
+ * @brief Computes the adjoint part of the gradient by solving a linear system.
+ *
+ * Valid only if the computed solution is a steady-state.
+ */
+ void computeIntegralForSteadyState();
+
  /**
  * @brief Compute likelihood sensitivities.
  */

include/amici/model_dae.h

@@ -92,6 +92,10 @@ class Model_DAE : public Model {
  void fJ(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector xdot,
  SUNMatrix J);
 
+ void fJB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB, const AmiVector &dxB,
+ const AmiVector &xBdot, SUNMatrix JB) override;
+
  /**
  * @brief Jacobian of xBdot with respect to adjoint state xB
  * @param t timepoint
@@ -102,7 +106,6 @@ class Model_DAE : public Model {
  * @param dxB Vector with the adjoint derivative states
  * @param JB Matrix to which the Jacobian will be written
  **/
-
  void fJB(realtype t, realtype cj, N_Vector x, N_Vector dx, N_Vector xB,
  N_Vector dxB, SUNMatrix JB);
 
@@ -121,6 +124,11 @@ class Model_DAE : public Model {
  void fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx,
  SUNMatrix J);
 
+ void fJSparseB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB,
+ const AmiVector &dxB, const AmiVector &xBdot,
+ SUNMatrix JB) override;
+
  /** JB in sparse form (for sparse solvers from the SuiteSparse Package)
  * @param t timepoint
  * @param cj scalar in Jacobian

include/amici/model_ode.h

@@ -94,6 +94,10 @@ class Model_ODE : public Model {
  **/
  void fJ(realtype t, N_Vector x, N_Vector xdot, SUNMatrix J);
 
+ void fJB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB, const AmiVector &dxB,
+ const AmiVector &xBdot, SUNMatrix JB) override;
+
  /** implementation of fJB at the N_Vector level, this function provides an
  *interface to the model specific routines for the solver implementation
  * @param t timepoint
@@ -119,6 +123,11 @@ class Model_ODE : public Model {
  */
  void fJSparse(realtype t, N_Vector x, SUNMatrix J);
 
+ void fJSparseB(const realtype t, realtype cj, const AmiVector &x,
+ const AmiVector &dx, const AmiVector &xB,
+ const AmiVector &dxB, const AmiVector &xBdot,
+ SUNMatrix JB) override;
+
  /** implementation of fJSparseB at the N_Vector level, this function
  * provides an interface to the model specific routines for the solver
  * implementation

include/amici/newton_solver.h

@@ -79,7 +79,7 @@ class NewtonSolver {
  void computeNewtonSensis(AmiVectorArray &sx);
 
  /**
- * Writes the Jacobian for the Newton iteration and passes it to the linear
+ * Writes the Jacobian (J) for the Newton iteration and passes it to the linear
  * solver
  *
  * @param ntry integer newton_try integer start number of Newton solver
@@ -88,6 +88,16 @@ class NewtonSolver {
  */
  virtual void prepareLinearSystem(int ntry, int nnewt) = 0;
 
+ /**
+ * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear
+ * solver
+ *
+ * @param ntry integer newton_try integer start number of Newton solver
+ * (1 or 2)
+ * @param nnewt integer number of current Newton step
+ */
+ virtual void prepareLinearSystemB(int ntry, int nnewt) = 0;
+
  /**
  * Solves the linear system for the Newton step
  *
@@ -125,7 +135,10 @@ class NewtonSolver {
  AmiVector *x;
  /** current state time derivative (DAE) */
  AmiVector dx;
-
+ /** current adjoint state */
+ AmiVector xB;
+ /** current adjoint state time derivative (DAE) */
+ AmiVector dxB;
 };
 
 /**
@@ -158,7 +171,7 @@ class NewtonSolverDense : public NewtonSolver {
  void solveLinearSystem(AmiVector &rhs) override;
 
  /**
- * Writes the Jacobian for the Newton iteration and passes it to the linear
+ * Writes the Jacobian (J) for the Newton iteration and passes it to the linear
  * solver
  *
  * @param ntry integer newton_try integer start number of Newton solver
@@ -167,6 +180,16 @@ class NewtonSolverDense : public NewtonSolver {
  */
  void prepareLinearSystem(int ntry, int nnewt) override;
 
+ /**
+ * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear
+ * solver
+ *
+ * @param ntry integer newton_try integer start number of Newton solver
+ * (1 or 2)
+ * @param nnewt integer number of current Newton step
+ */
+ void prepareLinearSystemB(int ntry, int nnewt) override;
+
  private:
  /** temporary storage of Jacobian */
  SUNMatrixWrapper Jtmp;
@@ -204,7 +227,7 @@ class NewtonSolverSparse : public NewtonSolver {
  void solveLinearSystem(AmiVector &rhs) override;
 
  /**
- * Writes the Jacobian for the Newton iteration and passes it to the linear
+ * Writes the Jacobian (J) for the Newton iteration and passes it to the linear
  * solver
  *
  * @param ntry integer newton_try integer start number of Newton solver
@@ -213,6 +236,16 @@ class NewtonSolverSparse : public NewtonSolver {
  */
  void prepareLinearSystem(int ntry, int nnewt) override;
 
+ /**
+ * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear
+ * solver
+ *
+ * @param ntry integer newton_try integer start number of Newton solver
+ * (1 or 2)
+ * @param nnewt integer number of current Newton step
+ */
+ void prepareLinearSystemB(int ntry, int nnewt) override;
+
  private:
  /** temporary storage of Jacobian */
  SUNMatrixWrapper Jtmp;
@@ -249,7 +282,7 @@ class NewtonSolverIterative : public NewtonSolver {
  void solveLinearSystem(AmiVector &rhs) override;
 
  /**
- * Writes the Jacobian for the Newton iteration and passes it to the linear
+ * Writes the Jacobian (J) for the Newton iteration and passes it to the linear
  * solver.
  * Also wraps around getSensis for iterative linear solver.
  *
@@ -259,6 +292,17 @@ class NewtonSolverIterative : public NewtonSolver {
  */
  void prepareLinearSystem(int ntry, int nnewt) override;
 
+ /**
+ * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear
+ * solver.
+ * Also wraps around getSensis for iterative linear solver.
+ *
+ * @param ntry integer newton_try integer start number of Newton solver
+ * (1 or 2)
+ * @param nnewt integer number of current Newton step
+ */
+ void prepareLinearSystemB(int ntry, int nnewt) override;
+
  /**
  * Iterative linear solver created from SPILS BiCG-Stab.
  * Solves the linear system within each Newton step if iterative solver is
@@ -296,6 +340,8 @@ class NewtonSolverIterative : public NewtonSolver {
  AmiVector ns_tmp;
  /** ??? */
  AmiVector ns_Jdiag;
+ /** temporary storage of Jacobian */
+ SUNMatrixWrapper ns_J;
 };