Add support for observable-dependent sigmas

include/amici/abstract_model.h

@@ -561,20 +561,37 @@ class AbstractModel {
      * @param t current time
      * @param p parameter vector
      * @param k constant vector
+     * @param y model output at timepoint t
      */
     virtual void fsigmay(realtype *sigmay, const realtype t, const realtype *p,
-                         const realtype *k);
+                         const realtype *k, const realtype *y);
 
     /**
-     * @brief Model-specific implementation of fsigmay
+     * @brief Model-specific implementation of fdsigmaydp
      * @param dsigmaydp partial derivative of standard deviation of measurements
      * @param t current time
      * @param p parameter vector
      * @param k constant vector
+     * @param y model output at timepoint t
      * @param ip sensitivity index
      */
     virtual void fdsigmaydp(realtype *dsigmaydp, const realtype t,
-                            const realtype *p, const realtype *k, int ip);
+                            const realtype *p, const realtype *k,
+                            const realtype *y, int ip);
+    /**
+     * @brief Model-specific implementation of fsigmay
+     * @param dsigmaydy partial derivative of standard deviation of measurements
+     * w.r.t. model outputs
+     * @param t current time
+     * @param p parameter vector
+     * @param k constant vector
+     * @param y model output at timepoint t
+     * @param ip sensitivity index
+     */
+    virtual void fdsigmaydy(realtype *dsigmaydy, const realtype t,
+                            const realtype *p, const realtype *k,
+                            const realtype *y, int ip);
+
 
     /**
      * @brief Model-specific implementation of fsigmaz

include/amici/model.h

@@ -115,6 +115,7 @@ class Model : public AbstractModel, public ModelDimensions {
     using AbstractModel::fdrzdp;
     using AbstractModel::fdrzdx;
     using AbstractModel::fdsigmaydp;
+    using AbstractModel::fdsigmaydy;
     using AbstractModel::fdsigmazdp;
     using AbstractModel::fdwdp;
     using AbstractModel::fdwdp_colptrs;
@@ -851,11 +852,13 @@ class Model : public AbstractModel, public ModelDimensions {
      * Total derivative (can be used with both adjoint and forward sensitivity).
      *
      * @param ssigmay Buffer (shape `ny` x `nplist`, row-major)
+     * @param sy Sensitivity of time-resolved observables for current timepoint
      * @param it Timepoint index
      * @param edata Pointer to experimental data instance (optional, pass
      * `nullptr` to ignore)
      */
     void getObservableSigmaSensitivity(gsl::span<realtype> ssigmay,
+                                       gsl::span<const realtype> sy,
                                        const int it, const ExpData *edata);
 
     /**
@@ -1396,6 +1399,14 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void fdsigmaydp(int it, const ExpData *edata);
 
+    /**
+     * @brief Compute partial derivative of standard deviation of measurements
+     * w.r.t. model outputs.
+     * @param it Timepoint index
+     * @param edata pointer to `amici::ExpData` data instance holding sigma values
+     */
+    void fdsigmaydy(int it, const ExpData *edata);
+
     /**
      * @brief Compute negative log-likelihood of measurements \f$ y \f$.
      *

include/amici/model_state.h

@@ -228,11 +228,16 @@ struct ModelStateDerived {
     /** data standard deviation for current timepoint (dimension: ny) */
     std::vector<realtype> sigmay_;
 
-    /** temporary storage for  parameter derivative of data standard deviation,
+    /** temporary storage for parameter derivative of data standard deviation,
      * (dimension: ny x nplist, row-major)
      */
     std::vector<realtype> dsigmaydp_;
 
+    /** temporary storage for observable derivative of data standard deviation,
+     * (dimension: ny x ny, row-major)
+     */
+    std::vector<realtype> dsigmaydy_;
+
     /** temporary storage for event-resolved observable (dimension: nz) */
     std::vector<realtype> z_;
 

matlab/@amifun/getArgs.m

@@ -80,9 +80,9 @@
         case 'dxdotdp'
             this.argstr = ['(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip' dx ', const realtype *w, const realtype *dwdp)'];
         case 'sigma_y'
-            this.argstr = '(double *sigmay, const realtype t, const realtype *p, const realtype *k)';
+            this.argstr = '(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y)';
         case 'dsigma_ydp'
-            this.argstr = '(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip)';
+            this.argstr = '(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip)';
         case 'sigma_z'
             this.argstr = '(double *sigmaz, const realtype t, const realtype *p, const realtype *k)';
         case 'dsigma_zdp'

models/model_calvetti/model_calvetti.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_calvetti_h
 #define _amici_model_calvetti_h
-/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */
+/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -22,7 +22,7 @@ extern void dJydy_model_calvetti(double *dJydy, const int iy, const realtype *p,
 extern void dwdx_model_calvetti(realtype *dwdx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 extern void dydx_model_calvetti(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx);
 extern void root_model_calvetti(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx);
-extern void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k);
+extern void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y);
 extern void w_model_calvetti(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
 extern void x0_model_calvetti(realtype *x0, const realtype t, const realtype *p, const realtype *k);
 extern void xdot_model_calvetti(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *w);
@@ -67,7 +67,7 @@ class Model_model_calvetti : public amici::Model_DAE {
 
     amici::Model* clone() const override { return new Model_model_calvetti(*this); };
 
-    std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; };
+    std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; };
 
     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);
@@ -125,7 +125,7 @@ class Model_model_calvetti : public amici::Model_DAE {
     void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {
     }
 
-    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
+    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override {
     }
 
     void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
@@ -161,8 +161,8 @@ class Model_model_calvetti : public amici::Model_DAE {
     void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {
     }
 
-    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override {
-        sigmay_model_calvetti(sigmay, t, p, k);
+    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override {
+        sigmay_model_calvetti(sigmay, t, p, k, y);
     }
 
     void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override {

models/model_calvetti/model_calvetti_sigmay.cpp

@@ -10,7 +10,7 @@ namespace amici {
 
 namespace model_model_calvetti{
 
-void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k) {
+void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) {
   sigmay[0] = 1.0;
   sigmay[1] = 1.0;
   sigmay[2] = 1.0;

models/model_dirac/model_dirac.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_dirac_h
 #define _amici_model_dirac_h
-/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */
+/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -23,7 +23,7 @@ extern void deltax_model_dirac(double *deltax, const realtype t, const realtype
 extern void dxdotdp_model_dirac(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp);
 extern void dydx_model_dirac(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx);
 extern void root_model_dirac(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
-extern void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k);
+extern void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y);
 extern void stau_model_dirac(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie);
 extern void xdot_model_dirac(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
 extern void y_model_dirac(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w);
@@ -67,7 +67,7 @@ class Model_model_dirac : public amici::Model_ODE {
 
     amici::Model* clone() const override { return new Model_model_dirac(*this); };
 
-    std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; };
+    std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; };
 
     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);
@@ -123,7 +123,7 @@ class Model_model_dirac : public amici::Model_ODE {
     void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {
     }
 
-    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
+    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override {
     }
 
     void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
@@ -159,8 +159,8 @@ class Model_model_dirac : public amici::Model_ODE {
     void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {
     }
 
-    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override {
-        sigmay_model_dirac(sigmay, t, p, k);
+    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override {
+        sigmay_model_dirac(sigmay, t, p, k, y);
     }
 
     void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override {

models/model_dirac/model_dirac_sigmay.cpp

@@ -10,7 +10,7 @@ namespace amici {
 
 namespace model_model_dirac{
 
-void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k) {
+void sigmay_model_dirac(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) {
   sigmay[0] = 1.0;
 }
 

models/model_events/model_events.h

@@ -1,6 +1,6 @@
 #ifndef _amici_model_events_h
 #define _amici_model_events_h
-/* Generated by amiwrap (R2017b) 9212bbdaf5712727ff284cc475f16d2a983f9bf2 */
+/* Generated by amiwrap (R2017b) e8e46ac8ec0c895549703eace15f2135fdb1747b */
 #include <cmath>
 #include <memory>
 #include "amici/defines.h"
@@ -32,7 +32,7 @@ extern void dydx_model_events(double *dydx, const realtype t, const realtype *x,
 extern void dzdx_model_events(double *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
 extern void root_model_events(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
 extern void rz_model_events(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h);
-extern void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k);
+extern void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y);
 extern void sigmaz_model_events(double *sigmaz, const realtype t, const realtype *p, const realtype *k);
 extern void srz_model_events(double *srz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip);
 extern void stau_model_events(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl, const realtype *sx, const int ip, const int ie);
@@ -81,7 +81,7 @@ class Model_model_events : public amici::Model_ODE {
 
     amici::Model* clone() const override { return new Model_model_events(*this); };
 
-    std::string getAmiciCommit() const override { return "9212bbdaf5712727ff284cc475f16d2a983f9bf2"; };
+    std::string getAmiciCommit() const override { return "e8e46ac8ec0c895549703eace15f2135fdb1747b"; };
 
     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);
@@ -143,7 +143,7 @@ class Model_model_events : public amici::Model_ODE {
         drzdx_model_events(drzdx, ie, t, x, p, k, h);
     }
 
-    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
+    void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const realtype *y, const int ip) override {
     }
 
     void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override {
@@ -182,8 +182,8 @@ class Model_model_events : public amici::Model_ODE {
         rz_model_events(rz, ie, t, x, p, k, h);
     }
 
-    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override {
-        sigmay_model_events(sigmay, t, p, k);
+    void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) override {
+        sigmay_model_events(sigmay, t, p, k, y);
     }
 
     void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override {

models/model_events/model_events_sigmay.cpp

@@ -10,7 +10,7 @@ namespace amici {
 
 namespace model_model_events{
 
-void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k) {
+void sigmay_model_events(double *sigmay, const realtype t, const realtype *p, const realtype *k, const realtype *y) {
   sigmay[0] = 1.0;
 }