Let froot depend on w

WIP

This should hopefully allow for more efficient computation of `root`-derivatives
by avoiding flattening `w` into `root` which, so far, made computing `drootdt_total`
prohibitively expensive in case of large `w` dependencies in `root`.

Author: dweindl

include/amici/abstract_model.h

@@ -299,6 +299,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param tcl total abundances for conservation laws
      * @param sx current state sensitivity
@@ -307,8 +308,9 @@ class AbstractModel {
      */
     virtual void fstau(
         realtype* stau, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx,
-        realtype const* tcl, realtype const* sx, int ip, int ie
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, realtype const* tcl, realtype const* sx, int ip,
+        int ie
     );
 
     /**
@@ -542,6 +544,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ie event index
      * @param xdot new model right hand side
@@ -552,9 +555,10 @@ class AbstractModel {
      */
     virtual void fdeltaxB(
         realtype* deltaxB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ie,
-        realtype const* xdot, realtype const* xdot_old, realtype const* x_old,
-        realtype const* xB, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB,
+        realtype const* tcl
     );
 
     /**
@@ -565,6 +569,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ip sensitivity index
      * @param ie event index
@@ -575,9 +580,9 @@ class AbstractModel {
      */
     virtual void fdeltaqB(
         realtype* deltaqB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ip,
-        int ie, realtype const* xdot, realtype const* xdot_old,
-        realtype const* x_old, realtype const* xB
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ip, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB
     );
 
     /**

include/amici/model_dae.h

@@ -359,11 +359,13 @@ class Model_DAE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx Vector with the derivative states
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, double const* p,
-        double const* k, realtype const* h, realtype const* dx
+        double const* k, realtype const* h, realtype const* w,
+        realtype const* dx
     );
 
     /**

include/amici/model_ode.h

@@ -316,11 +316,13 @@ class Model_ODE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param tcl total abundances for conservation laws
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* tcl
     );
 
     /**

models/model_calvetti_py/model_calvetti_py.h

@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_calvetti_py(SUNMatrixWrapper &rowvals, int index
 
 
 
-extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -200,10 +200,10 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdeltasx(realtype *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, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -323,8 +323,8 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_calvetti_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_calvetti_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -339,7 +339,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -557,7 +557,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "b0b2684b4b67db9eadf5e47d4f87f8fe74dd9070";
+        return "unknown";
     }
 
     bool has_quadratic_llh() const override {

models/model_calvetti_py/root.cpp

@@ -5,11 +5,12 @@
 #include "x.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 
 namespace amici {
 namespace model_model_calvetti_py {
 
-void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = t - 10;
     root[1] = 10 - t;
     root[2] = 12 - t;

models/model_dirac_py/deltaqB.cpp

@@ -5,6 +5,7 @@
 #include "x.h"
 #include "p.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -13,7 +14,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 0:
             switch(ip) {

models/model_dirac_py/model_dirac_py.h

@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_dirac_py(SUNMatrixWrapper &rowvals, int index);
 
 
 
-extern void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -77,11 +77,11 @@ extern void xdot_model_dirac_py(realtype *xdot, const realtype t, const realtype
 extern void y_model_dirac_py(realtype *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 
 
-extern void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
+extern void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 extern void deltax_model_dirac_py(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, const realtype *x_old);
 extern void deltasx_model_dirac_py(realtype *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, const realtype *x_old);
 
-extern void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
+extern void deltaqB_model_dirac_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB);
 
 extern void x_solver_model_dirac_py(realtype *x_solver, const realtype *x_rdata);
 
@@ -201,11 +201,11 @@ class Model_model_dirac_py : public amici::Model_ODE {
     }
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
-        deltaqB_model_dirac_py(deltaqB, t, x, p, k, h, dx, ip, ie, xdot, xdot_old, x_old, xB);
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {
+        deltaqB_model_dirac_py(deltaqB, t, x, p, k, h, w, dx, ip, ie, xdot, xdot_old, x_old, xB);
     }
 
 
@@ -314,8 +314,8 @@ class Model_model_dirac_py : public amici::Model_ODE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_dirac_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_dirac_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -330,8 +330,8 @@ class Model_model_dirac_py : public amici::Model_ODE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
-        stau_model_dirac_py(stau, t, x, p, k, h, dx, tcl, sx, ip, ie);
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {
+        stau_model_dirac_py(stau, t, x, p, k, h, w, dx, tcl, sx, ip, ie);
     }
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
@@ -544,7 +544,7 @@ class Model_model_dirac_py : public amici::Model_ODE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "b0b2684b4b67db9eadf5e47d4f87f8fe74dd9070";
+        return "unknown";
     }
 
     bool has_quadratic_llh() const override {

models/model_dirac_py/root.cpp

@@ -9,7 +9,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_dirac_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = -p2 + t;
 }
 

models/model_dirac_py/stau.cpp

@@ -10,7 +10,7 @@
 namespace amici {
 namespace model_model_dirac_py {
 
-void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
+void stau_model_dirac_py(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie){
     switch(ie) {
         case 0:
             switch(ip) {

models/model_events_py/deltaqB.cpp

@@ -6,6 +6,7 @@
 #include "p.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 #include "xdot.h"
 #include "xdot_old.h"
 #include "x_old.h"
@@ -14,7 +15,7 @@
 namespace amici {
 namespace model_model_events_py {
 
-void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
+void deltaqB_model_events_py(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB){
     switch(ie) {
         case 2:
         case 3: