Python splines

CMakeLists.txt

@@ -121,6 +121,7 @@ add_custom_target(
 # Library source files
 set(AMICI_SRC_LIST
     ${CMAKE_SOURCE_DIR}/src/symbolic_functions.cpp
+    ${CMAKE_SOURCE_DIR}/src/splinefunctions.cpp
     ${CMAKE_SOURCE_DIR}/src/cblas.cpp
     ${CMAKE_SOURCE_DIR}/src/amici.cpp
     ${CMAKE_SOURCE_DIR}/src/misc.cpp

include/amici/abstract_model.h

@@ -4,6 +4,7 @@
 #include "amici/defines.h"
 #include "amici/sundials_matrix_wrapper.h"
 #include "amici/vector.h"
+#include "amici/splinefunctions.h"
 
 #include <sunmatrix/sunmatrix_band.h>
 #include <sunmatrix/sunmatrix_dense.h>
@@ -728,10 +729,11 @@ class AbstractModel {
      * @param k constants vector
      * @param h Heaviside vector
      * @param tcl total abundances for conservation laws
+     * @param spl spline value vector
      */
     virtual void fw(realtype *w, const realtype t, const realtype *x,
                     const realtype *p, const realtype *k, const realtype *h,
-                    const realtype *tcl);
+                    const realtype *tcl, const realtype *spl);
 
     /**
      * @brief Model-specific sparse implementation of dwdp
@@ -744,11 +746,14 @@ class AbstractModel {
      * @param w vector with helper variables
      * @param tcl total abundances for conservation laws
      * @param stcl sensitivities of total abundances for conservation laws
+     * @param spl spline value vector
+     * @param sspl sensitivities of spline values vector
      */
     virtual void fdwdp(realtype *dwdp, const realtype t, const realtype *x,
                        const realtype *p, const realtype *k, const realtype *h,
                        const realtype *w, const realtype *tcl,
-                       const realtype *stcl);
+                       const realtype *stcl, const realtype *spl, 
+                       const realtype *sspl);
 
     /**
      * @brief Model-specific implementation for dwdp, column pointers
@@ -778,7 +783,8 @@ class AbstractModel {
     virtual void fdwdp(realtype *dwdp, const realtype t, const realtype *x,
                        const realtype *p, const realtype *k, const realtype *h,
                        const realtype *w, const realtype *tcl,
-                       const realtype *stcl, int ip);
+                       const realtype *stcl, const realtype *spl, 
+                       const realtype *sspl, int ip);
 
     /**
      * @brief Model-specific implementation of dwdx, data part
@@ -790,10 +796,12 @@ class AbstractModel {
      * @param h Heaviside vector
      * @param w vector with helper variables
      * @param tcl total abundances for conservation laws
+     * @param spl spline value vector
      */
     virtual void fdwdx(realtype *dwdx, const realtype t, const realtype *x,
                        const realtype *p, const realtype *k, const realtype *h,
-                       const realtype *w, const realtype *tcl);
+                       const realtype *w, const realtype *tcl, 
+                       const realtype *spl);
 
     /**
      * @brief Model-specific implementation for dwdx, column pointers
@@ -833,6 +841,34 @@ class AbstractModel {
      * @param dwdw sparse matrix to which rowvals will be written
      */
     virtual void fdwdw_rowvals(SUNMatrixWrapper &dwdw);
+
+    /**
+     * @brief Model-specific implementation the spline constructors
+     * @param p parameter vector
+     * @param k constants vector
+     */
+    virtual std::vector<HermiteSpline> fspline_constructors(const realtype *p, 
+                                                            const realtype *k);
+
+    /**
+     * @brief Model-specific implementation the parametric derivatives
+     * of spline node values
+     * @param dspline_valuesdp vector to which derivatives will be written
+     * @param p parameter vector
+     * @param k constants vector
+     */
+    virtual void fdspline_valuesdp(realtype *dspline_valuesdp,
+                                   const realtype *p, const realtype *k);
+
+    /**
+     * @brief Model-specific implementation the parametric derivatives
+     * of slopevalues at spline nodes
+     * @param dspline_slopesdp vector to which derivatives will be written
+     * @param p parameter vector
+     * @param k constants vector
+     */
+    virtual void fdspline_slopesdp(realtype *dspline_slopesdp,
+                                   const realtype *p, const realtype *k);
 };
 
 } // namespace amici

include/amici/defines.h

@@ -213,6 +213,24 @@ enum class RDataReporting {
     likelihood,
 };
 
+/** boundary conditions for splines */
+enum class SplineBoundaryCondition {
+    given                 = -1,
+    zeroDerivative        =  0,
+    natural               =  1,
+    naturalZeroDerivative =  2,
+    periodic              =  3,
+};
+
+/** extrapolation methods for splines */
+enum class SplineExtrapolation {
+    noExtrapolation = -1,
+    constant        =  0,
+    linear          =  1,
+    polynomial      =  2,
+    periodic        =  3,
+};
+
 /**
  * Type for function to process warnings or error messages.
  */

include/amici/model.h

@@ -8,6 +8,7 @@
 #include "amici/simulation_parameters.h"
 #include "amici/model_dimensions.h"
 #include "amici/model_state.h"
+#include "amici/splinefunctions.h"
 
 #include <map>
 #include <memory>
@@ -181,6 +182,13 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     void initializeStateSensitivities(AmiVectorArray &sx, const AmiVector &x);
 
+    /**
+     * @brief Initialization of spline functions
+     */
+    void initializeSplines();
+
+    void initializeSplineSensitivities();
+
     /**
      * @brief Initialize the Heaviside variables `h` at the initial time `t0`.
      *
@@ -642,7 +650,7 @@ class Model : public AbstractModel, public ModelDimensions {
             throw AmiException("Mismatch in conservation law sensitivity size");
         state_ = state;
     };
-    
+
     /**
      * @brief Sets the estimated lower boundary for sigma_y. When :meth:`setAddSigmaResiduals` is
      * activated, this lower boundary must ensure that log(sigma) + min_sigma > 0.
@@ -651,15 +659,15 @@ class Model : public AbstractModel, public ModelDimensions {
     void setMinimumSigmaResiduals(double min_sigma) {
         min_sigma_ = min_sigma;
     }
-    
+
     /**
      * @brief Gets the specified estimated lower boundary for sigma_y.
      * @return lower boundary
      */
     realtype getMinimumSigmaResiduals() const {
         return min_sigma_;
     }
-    
+
     /**
      * @brief Specifies whether residuals should be added to account for parameter dependent sigma.
      *
@@ -673,7 +681,7 @@ class Model : public AbstractModel, public ModelDimensions {
     void setAddSigmaResiduals(bool sigma_res) {
         sigma_res_ = sigma_res;
     }
-    
+
     /**
      * @brief Checks whether residuals should be added to account for parameter dependent sigma.
      * @return sigma_res
@@ -1589,6 +1597,18 @@ class Model : public AbstractModel, public ModelDimensions {
     void fdJrzdsigma(const int ie, const int nroots, const realtype t,
                      const AmiVector &x, const ExpData &edata);
 
+    /**
+     * @brief Spline functions
+     * @param t timepoint
+     */
+    void fspl(realtype t);
+
+    /**
+     * @brief Parametric derivatives of splines functions
+     * @param t timepoint
+     */
+    void fsspl(realtype t);
+
     /**
      * @brief Compute recurring terms in xdot.
      * @param t Timepoint
@@ -1713,6 +1733,9 @@ class Model : public AbstractModel, public ModelDimensions {
      */
     ModelStateDerived derived_state_;
 
+    /** Storage for splines of the model */
+    std::vector<HermiteSpline> splines_;
+
     /** index indicating to which event an event output belongs */
     std::vector<int> z2event_;
 
@@ -1743,10 +1766,10 @@ class Model : public AbstractModel, public ModelDimensions {
      * checked for finiteness
      */
     bool always_check_finite_ {false};
-    
+
     /** indicates whether sigma residuals are to be added for every datapoint  */
     bool sigma_res_ {false};
-    
+
     /** offset to ensure positivity of sigma residuals, only has an effect when `sigma_res_` is `true`  */
     realtype min_sigma_ {50.0};
 

include/amici/model_dimensions.h

@@ -124,6 +124,9 @@ struct ModelDimensions {
     /** Number of event outputs in the unaugmented system */
     int nztrue{0};
 
+    /** numer of spline functions in the model */
+    int nspl{0};
+
     /** Number of events */
     int ne{0};
 

include/amici/model_state.h

@@ -43,6 +43,9 @@ struct ModelState {
      * (dimension: nplist)
      */
     std::vector<int> plist;
+
+    /** temporary storage for spline values */
+    std::vector<realtype> spl_;
 };
 
 
@@ -263,6 +266,9 @@ struct ModelStateDerived {
      */
     std::vector<realtype> deltaqB_;
 
+    /** temporary storage for sensitivity values of splines */
+    SUNMatrixWrapper sspl_;
+
     /** temporary storage of positified state variables according to
      * stateIsNonNegative (dimension: `nx_solver`) */
     AmiVector x_pos_tmp_ {0};