Adjoint event

include/amici/abstract_model.h

@@ -74,6 +74,21 @@ class AbstractModel {
         AmiVector const& sx, AmiVector const& sdx, AmiVector& sxdot
     ) = 0;
 
+    /**
+     * @brief Residual function of adjoint state
+     * @param t time
+     * @param x state
+     * @param dx time derivative of state (DAE only)
+     * @param xB adjoint state
+     * @param dxB time derivative of adjoint state (DAE only)
+     * @param xBdot array to which values of the residual function will be
+     * written
+     */
+    virtual void fxBdot(
+        realtype t, AmiVector const& x, AmiVector const& dx,
+        AmiVector const& xB, AmiVector const& dxB, AmiVector& xBdot
+    ) = 0;
+
     /**
      * @brief Residual function backward when running in steady state mode
      * @param t time
@@ -569,11 +584,14 @@ class AbstractModel {
      * @param xdot new model right hand side
      * @param xdot_old previous model right hand side
      * @param xB current adjoint state
+     * @param xBdot right hand side of adjoint state
+     * @param tcl total abundances for conservation laws
      */
     virtual void fdeltaxB(
         realtype* deltaxB, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ie,
-        realtype const* xdot, realtype const* xdot_old, realtype const* xB
+        realtype const* xdot, realtype const* xdot_old, realtype const* xB,
+        realtype const* xBdot, realtype const* tcl
     );
 
     /**
@@ -589,11 +607,13 @@ class AbstractModel {
      * @param xdot new model right hand side
      * @param xdot_old previous model right hand side
      * @param xB adjoint state
+     * @param xBdot right hand side of adjoint state
      */
     virtual void fdeltaqB(
         realtype* deltaqB, realtype const t, realtype const* x,
         realtype const* p, realtype const* k, realtype const* h, int ip, int ie,
-        realtype const* xdot, realtype const* xdot_old, realtype const* xB
+        realtype const* xdot, realtype const* xdot_old, realtype const* xB,
+        realtype const* xBdot
     );
 
     /**

include/amici/backwardproblem.h

@@ -116,6 +116,8 @@ class BackwardProblem {
     AmiVector xQB_;
     /** array of state vectors at discontinuities*/
     std::vector<AmiVector> x_disc_;
+    /** array of old state vectors at discontinuities*/
+    std::vector<AmiVector> x_old_disc_;
     /** array of differential state vectors at discontinuities*/
     std::vector<AmiVector> xdot_disc_;
     /** array of old differential state vectors at discontinuities*/

include/amici/forwardproblem.h

@@ -89,6 +89,14 @@ class ForwardProblem {
         return x_disc_;
     }
 
+    /**
+     * @brief Accessor for x_old_disc
+     * @return x_old_disc
+     */
+    std::vector<AmiVector> const& getStatesBeforeDiscontinuities() const {
+        return x_old_disc_;
+    }
+
     /**
      * @brief Accessor for xdot_disc
      * @return xdot_disc
@@ -346,6 +354,10 @@ class ForwardProblem {
      * discontinuities, extended as needed (dimension dynamic) */
     std::vector<AmiVector> x_disc_;
 
+    /** array of old state vectors (dimension nx) for all so far encountered
+     * discontinuities, extended as needed (dimension dynamic) */
+    std::vector<AmiVector> x_old_disc_;
+
     /** array of differential state vectors (dimension nx) for all so far
      * encountered discontinuities, extended as needed (dimension dynamic) */
     std::vector<AmiVector> xdot_disc_;

include/amici/model.h

@@ -1282,10 +1282,11 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param x Current state
      * @param xdot Current residual function values
      * @param xdot_old Value of residual function before event
+     * @param xBdot right hand side of adjoint state
      */
     void addAdjointStateEventUpdate(
         AmiVector& xB, int const ie, realtype const t, AmiVector const& x,
-        AmiVector const& xdot, AmiVector const& xdot_old
+        AmiVector const& xdot, AmiVector const& xdot_old, AmiVector const& xBdot
     );
 
     /**
@@ -1297,10 +1298,12 @@ class Model : public AbstractModel, public ModelDimensions {
      * @param xB Current adjoint state
      * @param xdot Current residual function values
      * @param xdot_old Value of residual function before event
+     * @param xBdot right hand side of adjoint state
      */
     void addAdjointQuadratureEventUpdate(
-        AmiVector xQB, int const ie, realtype const t, AmiVector const& x,
-        AmiVector const& xB, AmiVector const& xdot, AmiVector const& xdot_old
+        AmiVector& xQB, int const ie, realtype const t, AmiVector const& x,
+        AmiVector const& xB, AmiVector const& xdot, AmiVector const& xdot_old,
+        AmiVector const& xBdot
     );
 
     /**

include/amici/model_dae.h

@@ -218,6 +218,11 @@ class Model_DAE : public Model {
      */
     void fxdot(realtype t, const_N_Vector x, const_N_Vector dx, N_Vector xdot);
 
+    void fxBdot(
+        realtype t, AmiVector const& x, AmiVector const& dx,
+        AmiVector const& xB, AmiVector const& dxB, AmiVector& xBdot
+    ) override;
+
     /**
      * @brief Right hand side of differential equation for adjoint state xB
      * @param t timepoint

include/amici/model_ode.h

@@ -202,14 +202,20 @@ class Model_ODE : public Model {
      */
     void fxdot(realtype t, const_N_Vector x, N_Vector xdot);
 
+    void fxBdot(
+        realtype t, AmiVector const& x, AmiVector const& dx,
+        AmiVector const& xB, AmiVector const& dxB, AmiVector& xBdot
+    ) override;
+
     /**
      * @brief Implementation of fxBdot at the N_Vector level
      * @param t timepoint
      * @param x Vector with the states
      * @param xB Vector with the adjoint states
      * @param xBdot Vector with the adjoint right hand side
      */
-    void fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot);
+    void
+    fxBdot(realtype t, const_N_Vector x, const_N_Vector xB, N_Vector xBdot);
 
     /**
      * @brief Implementation of fqBdot at the N_Vector level

matlab/@amifun/getArgs.m

@@ -72,9 +72,9 @@
         case 'deltax'
             this.argstr = '(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old)';
         case 'deltaxB'
-            this.argstr = '(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB)';
+            this.argstr = '(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot, const realtype *tcl)';
         case 'deltaqB'
-            this.argstr = '(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB)';
+            this.argstr = '(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot)';
         case 'deltasx'
             this.argstr = '(double *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl)';
         case 'dxdotdp'

models/model_calvetti/model_calvetti.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_calvetti_h
 #define _amici_model_calvetti_h
-/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */
+/* Generated by amiwrap (R2017b) 7e32da0621a3696d3bed0e2ceb0a7176767f3830 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -72,7 +72,7 @@
 
     amici::Model* clone() const override { return new Model_model_calvetti(*this); };
 
-    std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; };
+    std::string getAmiciCommit() const override { return "7e32da0621a3696d3bed0e2ceb0a7176767f3830"; };
 
     void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override {
         JSparse_model_calvetti(JSparse, t, x, p, k, h, cj, dx, w, dwdx);
@@ -112,7 +112,7 @@
     void fdJzdz(double *dJzdz, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override {
     }
 
-    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot) override {
     }
 
     void fdeltasx(double *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl) override {
@@ -121,7 +121,7 @@
     void fdeltax(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old) override {
     }
 
-    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot, const realtype *tcl) override {
     }
 
     void fdrzdp(double *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {

models/model_dirac/model_dirac.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_dirac_h
 #define _amici_model_dirac_h
-/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */
+/* Generated by amiwrap (R2017b) 7e32da0621a3696d3bed0e2ceb0a7176767f3830 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -72,7 +72,7 @@
 
     amici::Model* clone() const override { return new Model_model_dirac(*this); };
 
-    std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; };
+    std::string getAmiciCommit() const override { return "7e32da0621a3696d3bed0e2ceb0a7176767f3830"; };
 
     void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         JSparse_model_dirac(JSparse, t, x, p, k, h, w, dwdx);
@@ -108,7 +108,7 @@
     void fdJzdz(double *dJzdz, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override {
     }
 
-    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot) override {
     }
 
     void fdeltasx(double *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl) override {
@@ -119,7 +119,7 @@
         deltax_model_dirac(deltax, t, x, p, k, h, ie, xdot, xdot_old);
     }
 
-    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot, const realtype *tcl) override {
     }
 
     void fdrzdp(double *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {

models/model_events/model_events.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_events_h
 #define _amici_model_events_h
-/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */
+/* Generated by amiwrap (R2017b) 7e32da0621a3696d3bed0e2ceb0a7176767f3830 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -86,7 +86,7 @@
 
     amici::Model* clone() const override { return new Model_model_events(*this); };
 
-    std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; };
+    std::string getAmiciCommit() const override { return "7e32da0621a3696d3bed0e2ceb0a7176767f3830"; };
 
     void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         JSparse_model_events(JSparse, t, x, p, k, h, w, dwdx);
@@ -128,7 +128,7 @@
         dJzdz_model_events(dJzdz, iz, p, k, z, sigmaz, mz);
     }
 
-    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot) override {
     }
 
     void fdeltasx(double *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl) override {
@@ -138,7 +138,7 @@
     void fdeltax(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old) override {
     }
 
-    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot, const realtype *tcl) override {
     }
 
     void fdrzdp(double *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {

models/model_jakstat_adjoint/model_jakstat_adjoint.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_jakstat_adjoint_h
 #define _amici_model_jakstat_adjoint_h
-/* Generated by amiwrap (R2017b) d69c026f8f8a89a13f1fa1307548f2b8f7045fe1 */
+/* Generated by amiwrap (R2017b) 7e32da0621a3696d3bed0e2ceb0a7176767f3830 */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -75,7 +75,7 @@
 
     amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); };
 
-    std::string getAmiciCommit() const override { return "d69c026f8f8a89a13f1fa1307548f2b8f7045fe1"; };
+    std::string getAmiciCommit() const override { return "7e32da0621a3696d3bed0e2ceb0a7176767f3830"; };
 
     void fJSparse(SUNMatrixContent_Sparse JSparse, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override {
         JSparse_model_jakstat_adjoint(JSparse, t, x, p, k, h, w, dwdx);
@@ -111,7 +111,7 @@
     void fdJzdz(double *dJzdz, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override {
     }
 
-    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot) override {
     }
 
     void fdeltasx(double *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl) override {
@@ -120,7 +120,7 @@
     void fdeltax(double *deltax, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old) override {
     }
 
-    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override {
+    void fdeltaxB(double *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB, const realtype *xBdot, const realtype *tcl) override {
     }
 
     void fdrzdp(double *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {