Feature preequilibration sensitivities via FSA (#410)

* feature(core) refactor createSteadystateSimSolver

* feature(core) enable preequilibration sensitivities via FSA

* feature(core) add option to control whether FSA is used to computed steady state sensitivities

* feature(test) add test cases for FSA steadystate sensis

* fix(test) fix solver setup routine in unittest

* fix(core) fix segfault

* feature(core) add monotonicity checks for setTimepoints

* feature(core) use std::is_sorted for monotonicity check.

* fix(core) include header for std::is_sorted

* fix(core) wrong header file for std::is_sorted

* fix(doc) update documentation

* feature(core) implement pos_pow to deal with positivity issues

* fix(core) fix pos_pow regex

* fix(steadystateproblem) use t0 not ts[0] for initializaiton of t after preequilibration

* fix(doc) fix typo in documentation

* fix(steadystateproblem) reduce sensi computation logic

* fix(sbml) enable linking of fsx0 in model header

* feature(sbml) properly process initial assignments of state variables

TODO: do this for other variables aswell.

* feature(steadystate) reset initialStates that are defined via fixedParameters after Preequilibration

* fix(python-sbml) fix identification of 0.0 values and specifiy variable name for x0_fixedParameters

* fix(steadystate) perform x0_fixedParameters update after updating the fixedParameters

* fix(python) fix some codacy issues

* fix(tests) fix codacy issues

* fix(matlab) fix compile error in CalcMD5.c

* feature(steadystate) add option for more sensitivity computation modes in the future

* fix(steadystate) remove unnecessary intialization from initial condition for non-preequilibration

* fix(matlab) fix matlab interface after update from enum classes

* fix(steadystate) fix typo and update hdf

Author: FFroehlich

include/amici/defines.h

@@ -135,14 +135,20 @@ enum class StateOrdering {
     COLAMD,
     natural
 };
+  
+/** Sensitivity computation mode in steadyStateProblem */
+enum class SteadyStateSensitivityMode {
+    newtonOnly,
+    simulationFSA
+};
 
-    /**
-     * @brief msgIdAndTxtFp
-     * @param identifier string with error message identifier
-     * @param format string with error message printf-style format
-     * @param ... arguments to be formatted
-     */
-    typedef void (*msgIdAndTxtFp)(const char *identifier, const char *format, ...);
+/**
+ * @brief msgIdAndTxtFp
+ * @param identifier string with error message identifier
+ * @param format string with error message printf-style format
+ * @param ... arguments to be formatted
+ */
+typedef void (*msgIdAndTxtFp)(const char *identifier, const char *format, ...);
 
 // clang-format on
 

include/amici/interface_matlab.h

@@ -45,13 +45,13 @@ ReturnDataMatlab *setupReturnData(mxArray *plhs[], int nlhs);
  */
 std::unique_ptr<ExpData> expDataFromMatlabCall(const mxArray *prhs[], const Model &model);
 
-void amici_dgemv(AMICI_BLAS_LAYOUT layout, AMICI_BLAS_TRANSPOSE TransA,
+void amici_dgemv(BLASLayout layout, BLASTranspose TransA,
                  const int M, const int N, const double alpha, const double *A,
                  const int lda, const double *X, const int incX,
                  const double beta, double *Y, const int incY);
 
-void amici_dgemm(AMICI_BLAS_LAYOUT layout, AMICI_BLAS_TRANSPOSE TransA,
-                 AMICI_BLAS_TRANSPOSE TransB, const int M, const int N,
+void amici_dgemm(BLASLayout layout, BLASTranspose TransA,
+                 BLASTranspose TransB, const int M, const int N,
                  const int K, const double alpha, const double *A,
                  const int lda, const double *B, const int ldb,
                  const double beta, double *C, const int ldc);

include/amici/model.h

@@ -172,8 +172,16 @@ namespace amici {
         virtual void fJv(realtype t, AmiVector *x, AmiVector *dx, AmiVector *xdot,
                              AmiVector *v, AmiVector *nJv, realtype cj) = 0;
 
+        /** Initial states
+         * @param x pointer to state variables
+         */
         void fx0(AmiVector *x);
 
+        /** Sets only those initial states that are specified via fixedParmeters
+         * @param x pointer to state variables
+         */
+        void fx0_fixedParameters(AmiVector *x);
+        
         /** Initial value for time derivative of states (only necessary for DAEs)
          * @param x0 Vector with the initial states
          * @param dx0 Vector to which the initial derivative states will be
@@ -627,6 +635,22 @@ namespace amici {
          * @return the ids
          */
         virtual std::vector<std::string> getObservableIds() const { return std::vector<std::string>(); }
+        
+        /**
+         * @brief sets the mode how sensitivities are computed in the steadystate simulation
+         * @param mode steadyStateSensitivityMode
+         */
+        void setSteadyStateSensitivityMode (const SteadyStateSensitivityMode mode) {
+            steadyStateSensitivityMode = mode;
+        }
+        
+        /**
+         * @brief gets the mode how sensitivities are computed in the steadystate simulation
+         * @return flag value
+         */
+        SteadyStateSensitivityMode getSteadyStateSensitivityMode () const {
+            return steadyStateSensitivityMode;
+        }
 
         /** number of states */
         const int nx;
@@ -708,6 +732,15 @@ namespace amici {
             throw AmiException("Requested functionality is not supported as (%s) is not implemented for this model!",__func__);
         }
 
+        /** model specific implementation of fx0
+         * @param x0 initial state
+         * @param t initial time
+         * @param p parameter vector
+         * @param k constant vector
+         **/
+        virtual void fx0_fixedParameters(realtype *x0, const realtype t, const realtype *p, const realtype *k) {
+        }
+        
         /** model specific implementation of fsx0
          * @param sx0 initial state sensitivities
          * @param t initial time
@@ -1250,6 +1283,10 @@ namespace amici {
 
         /** starting time */
         double tstart = 0.0;
+        
+        /** flag indicating whether steadystate sensivities are to be computed
+         via FSA when steadyStateSimulation is used */
+        SteadyStateSensitivityMode steadyStateSensitivityMode = SteadyStateSensitivityMode::newtonOnly;
 
     };
 

include/amici/newton_solver.h

@@ -32,7 +32,7 @@ class NewtonSolver {
 
     void getStep(int ntry, int nnewt, AmiVector *delta);
 
-    void getSensis(const int it, AmiVectorArray *sx);
+    void computeNewtonSensis(AmiVectorArray *sx);
 
     /**
      * Writes the Jacobian for the Newton iteration and passes it to the linear

include/amici/solver.h

@@ -74,7 +74,7 @@ class Solver {
      */
     virtual Solver* clone() const = 0;
 
-    void setup(ForwardProblem *fwd, Model *model);
+    void setup(AmiVector *x, AmiVector *dx, AmiVectorArray *sx, AmiVectorArray *sdx, Model *model);
 
     void setupAMIB(BackwardProblem *bwd, Model *model);
 

include/amici/steadystateproblem.h

@@ -3,6 +3,7 @@
 
 #include "amici/defines.h"
 #include "amici/vector.h"
+#include "amici/solver_cvodes.h"
 #include <amici/newton_solver.h>
 
 #include <nvector/nvector_serial.h>
@@ -33,13 +34,13 @@ class SteadystateProblem {
     void applyNewtonsMethod(ReturnData *rdata, Model *model,
                                   NewtonSolver *newtonSolver, int newton_try);
 
-    void getNewtonOutput(ReturnData *rdata, int newton_status,
-                                double run_time, int it);
+    void getNewtonOutput(ReturnData *rdata, const Model *model,
+                         int newton_status, double run_time, int it);
 
     void getSteadystateSimulation(ReturnData *rdata, Solver *solver,
                                   Model *model, int it);
 
-    std::unique_ptr<void, std::function<void(void *)>> createSteadystateSimSolver(Solver *solver, Model *model, realtype tstart);
+    std::unique_ptr<CVodeSolver> createSteadystateSimSolver(Solver *solver, Model *model, realtype tstart);
 
     /** default constructor
      * @param t pointer to time variable

include/amici/symbolic_functions.h

@@ -12,6 +12,8 @@ double Dmin(int id, double a, double b, double c);
 double max(double a, double b, double c);
 double Dmax(int id, double a, double b, double c);
 
+double pos_pow(double base, double exponent);
+    
 int isNaN(double what);
 int isInf(double what);
 double getNaN();

matlab/auxiliary/CalcMD5/CalcMD5.c

@@ -205,7 +205,7 @@ void MD5Init(MD5_CTX *context)
 void MD5Update(MD5_CTX *context, UCHAR *input, UINT inputLen)
 {
   UINT index, partLen;
-  int i, inputLenM63;
+  int inputLenM63;
 
   /*  Compute number of bytes mod 64: */
   index = (UINT)((context->count[0] >> 3) & 0x3F);
@@ -220,6 +220,7 @@ void MD5Update(MD5_CTX *context, UCHAR *input, UINT inputLen)
 
   /*  Transform as many times as possible: */
   if (inputLen >= partLen) {
+    int i;
     memcpy((POINTER)&context->buffer[index], (POINTER)input, partLen);
     MD5Transform(context->state, context->buffer);
 

matlab/auxiliary/compileAMICIDependencies.m

@@ -61,6 +61,7 @@
 function includesstr = getIncludeString(amici_root_path, sundials_path, ssparse_path)
     includesstr = '';
     includesstr = strcat(includesstr,' -I"', fullfile(sundials_path, 'include'), '"');
+    includesstr = strcat(includesstr,' -I"', fullfile(sundials_path, 'src'), '"');
     includesstr = strcat(includesstr,' -I"', fullfile(amici_root_path), '"');
     includesstr = strcat(includesstr,' -I"', fullfile(amici_root_path, 'src'), '"');
     includesstr = strcat(includesstr,' -I"', fullfile(amici_root_path, 'include'), '"');

python/amici/__init__.py

@@ -89,11 +89,11 @@ def rdataToNumPyArrays(rdata):
         
     """
     npReturnData = {'ptr': rdata}
-    fieldNames = ['t', 'x', 'x0', 'sx', 'sx0', 'y', 'sigmay', 'sy', 'ssigmay', 
-                  'z', 'rz', 'sigmaz', 'sz', 'srz', 'ssigmaz', 'sllh', 's2llh', 
+    fieldNames = ['t', 'x', 'x0', 'sx', 'sx0', 'y', 'sigmay', 'sy', 'ssigmay',
+                  'z', 'rz', 'sigmaz', 'sz', 'srz', 'ssigmaz', 'sllh', 's2llh',
                   'J', 'xdot', 'status', 'llh', 'chi2', 'res', 'sres', 'FIM',
-                  'newton_numlinsteps', 'newton_numsteps', 
-                  'numsteps', 'numrhsevals', 'numerrtestfails', 'numnonlinsolvconvfails', 
+                  'newton_numlinsteps', 'newton_numsteps',
+                  'numsteps', 'numrhsevals', 'numerrtestfails', 'numnonlinsolvconvfails',
                   'order', 'numstepsB', 'numrhsevalsB', 'numerrtestfailsB', 'numnonlinsolvconvfailsB']
 
     for field in fieldNames: