Merge pull request #688 from ICB-DCM/develop

0.10.6 Release

.gitignore

@@ -34,6 +34,8 @@ models/model_events/build/*
 models/model_nested_events/build/*
 !models/model_robertson
 models/model_robertson/build/*
+!models/model_calvetti
+models/model_calvetti/build/*
 
 simulate_model_*_hdf.m
 simulate_model_*.m
@@ -80,6 +82,8 @@ python/.idea/*
 *.log
 *.pyc
 *.so
+*.bak
+*.sav
 
 matlab/mtoc/config/latexextras.sty-e
 matlab/mtoc/config/latexextras.sty
@@ -154,3 +158,5 @@ ThirdParty/mtocpp-master*
 ThirdParty/sundials/build/*
 ThirdParty/SuiteSparse/lib/*
 ThirdParty/SuiteSparse/share/*
+ThirdParty/SuperLU_MT_3.1/
+ThirdParty/superlu_mt_3.1.tar.gz

CMakeLists.txt

@@ -66,6 +66,17 @@ set(SUNDIALS_LIBRARIES
     )
 set(SUNDIALS_INCLUDE_DIRS "${CMAKE_SOURCE_DIR}/ThirdParty/sundials/build/include")
 
+option(SUNDIALS_SUPERLUMT_ENABLE "Enable sundials SuperLUMT?" OFF)
+if(SUNDIALS_SUPERLUMT_ENABLE)
+    set(SUNDIALS_LIBRARIES ${SUNDIALS_LIBRARIES}
+        ${SUNDIALS_LIB_DIR}/libsundials_sunlinsolsuperlumt${CMAKE_STATIC_LIBRARY_SUFFIX}
+        ${CMAKE_SOURCE_DIR}/ThirdParty/SuperLU_MT_3.1/lib/libsuperlu_mt_PTHREAD${CMAKE_STATIC_LIBRARY_SUFFIX}
+        -lblas
+        )
+    set(SUNDIALS_INCLUDE_DIRS ${SUNDIALS_INCLUDE_DIRS}
+        "${CMAKE_SOURCE_DIR}/ThirdParty/SuperLU_MT_3.1/SRC/")
+endif()
+
 set(GSL_LITE_INCLUDE_DIR "${CMAKE_SOURCE_DIR}/ThirdParty/gsl")
 
 # AMICI requires BLAS, currently Intel MKL, CBLAS or MATLAB BLAS can be used.

INSTALL.md

@@ -7,18 +7,16 @@ The sources for AMICI are accessible as
 - Source  [zip](https://github.com/ICB-DCM/AMICI/zipball/master)
 - GIT repository on  [github](https://github.com/ICB-DCM/AMICI)
 
-### Obtaining AMICI via the GIT versioning system
+If AMICI was downloaded as a zip, it needs to be unpacked in a convenient directory. If AMICI was obtained via cloning of the git repository, no further unpacking is necessary.
+
+### Obtaining AMICI via the GIT version control system
 In order to always stay up-to-date with the latest AMICI versions, simply pull it from our GIT repository and
 recompile it when a new release is available. For more information about GIT checkout their [website](http://git-scm.com/)
 
 The GIT repository can currently be found at https://github.com/ICB-DCM/AMICI and a direct clone is possible via
 
     git clone https://github.com/ICB-DCM/AMICI.git AMICI
 
-## Installation
-
-If AMICI was downloaded as a zip, it needs to be unpacked in a convenient directory. If AMICI was obtained via cloning of the git repository, no further unpacking is necessary.
-
 ### Dependencies
 
 The MATLAB interface only depends on the symbolic toolbox, which is needed for model compilation, but not simulation.
@@ -82,7 +80,7 @@ For the compilation of .mex files, MATLAB needs to be configured with a working
 For a list of supported compilers we refer to the mathworks documentation: [mathworks.com](http://mathworks.com/support/compilers/R2018b/index.html)
 Note that Microsoft Visual Studio compilers are currently not supported.
 
-###  python
+### Python
 
 To use AMICI from python, install the module and all other requirements using pip:
 
@@ -167,14 +165,24 @@ To use AMICI from C++, run the
     ./scripts/buildSuitesparse.sh
     ./scripts/buildAmici.sh
 
-script to compile amici libary. The static library file can then be linked from 
+script to compile amici library.
+
+The static library file can then be linked from
 
     ./build/libamici.a
 
 In CMake-based packages, amici can be linked via
 
     find_package(Amici)
 
+*Optional*: To build AMICI with SuperLU_MT support, run
+
+    ./scripts/buildSuperLUMT.sh
+    ./scripts/buildSundials.sh
+    cd build/
+    cmake -DSUNDIALS_SUPERLUMT_ENABLE=ON ..
+    make
+
 ## Dependencies
 
 The MATLAB interface requires the Mathworks Symbolic Toolbox for model generation via `amiwrap(...)`, but not for execution of precompiled models. Currently MATLAB R2018a or newer is not supported (see https://github.com/ICB-DCM/AMICI/issues/307)
@@ -199,3 +207,11 @@ __Algorithm 837: AMD__, an approximate minimum degree ordering algorithm, Patric
 
 __Algorithm 836: COLAMD__, a column approximate minimum degree ordering algorithm, Timothy A. Davis, John R. Gilbert, Stefan I. Larimore, Esmond G. Ng
 _ACM Transactions on Mathematical Software_, Vol 30, Issue 3, 2004, pp 377 - 380. [PDF](http://dl.acm.org/authorize?734450)
+
+
+### Optional
+
+__SuperLU_MT__, "a general purpose library for the direct solution of large,
+sparse, nonsymmetric systems of linear equations"
+(https://crd-legacy.lbl.gov/~xiaoye/SuperLU/#superlu_mt).
+SuperLU_MT is optional and is so far only available from the C++ interface.

documentation/FAQ.md

@@ -18,6 +18,12 @@ __A__: This may be due to an old compiler version. See [issue #161](https://gith
 
 ---
 
+__Q__: How are events interpreted in a DAE context?
+
+__A__: Currently we only support impulse free events. Also sensitivities have never been tested. Proceed with care and create an [issue](https://github.com/ICB-DCM/AMICI/issues) if any problems arise!
+
+---
+
 __Q__: The simulation/sensitivities I get are incorrect.
 
 __A__: There are some known issues, especially with adjoint sensitivities, events and DAEs. If your particular problem is not featured in the [issues](https://github.com/ICB-DCM/AMICI/issues) list, please add it!

include/amici/defines.h

@@ -102,7 +102,8 @@ enum class LinearSolver {
     SPGMR       = 6,
     SPBCG       = 7,
     SPTFQMR     = 8,
-    KLU         = 9
+    KLU         = 9,
+    SuperLUMT   = 10,
 };
 
 /** CVODES/IDAS forward sensitivity computation method */
@@ -131,20 +132,13 @@ enum class NonlinearSolverIteration {
     newton = 2
 };
 
-/** KLU state reordering */
-enum class StateOrdering {
-    AMD,
-    COLAMD,
-    natural
-};
-
 /** Sensitivity computation mode in steadyStateProblem */
 enum class SteadyStateSensitivityMode {
     newtonOnly,
     simulationFSA
 };
 
-/** State in which the steady state computionat finished */
+/** State in which the steady state computation finished */
 enum class NewtonStatus {
     failed=-1,
     newt=1,

include/amici/model.h

@@ -57,6 +57,8 @@ class Model : public AbstractModel {
      * @param ndwdp number of nonzero elements in the p derivative of the
      * repeating elements
      * @param ndxdotdw number of nonzero elements in the w derivative of xdot
+     * @param ndJydy number of nonzero elements in the y derivative of dJy
+     * (dimension nytrue)
      * @param nnz number of nonzero elements in Jacobian
      * @param ubw upper matrix bandwidth in the Jacobian
      * @param lbw lower matrix bandwidth in the Jacobian
@@ -70,7 +72,8 @@ class Model : public AbstractModel {
     Model(const int nx_rdata, const int nxtrue_rdata, const int nx_solver,
           const int nxtrue_solver, const int ny, const int nytrue, const int nz,
           const int nztrue, const int ne, const int nJ, const int nw,
-          const int ndwdx, const int ndwdp, const int ndxdotdw, const int nnz,
+          const int ndwdx, const int ndwdp, const int ndxdotdw,
+          std::vector<int> ndJydy, const int nnz,
           const int ubw, const int lbw, amici::SecondOrderMode o2mode,
           const std::vector<amici::realtype> &p, std::vector<amici::realtype> k,
           const std::vector<int> &plist, std::vector<amici::realtype> idlist,
@@ -132,6 +135,29 @@ class Model : public AbstractModel {
     using AbstractModel::fy;
     using AbstractModel::fz;
 
+    /**
+     * Model specific implementation of fdJydy colptrs
+     * @param indexptrs column pointers
+     * @param index ytrue index
+     */
+    virtual void fdJydy_colptrs(sunindextype *indexptrs, int index) {
+        throw AmiException("Requested functionality is not supported as %s "
+                           "is not implemented for this model!",
+                           __func__); // not implemented
+    }
+
+    /**
+     * Model specific implementation of fdxdotdw row vals
+     * @param indexptrs row val pointers
+     * @param index ytrue index
+     */
+    virtual void fdJydy_rowvals(sunindextype *indexptrs, int index) {
+        throw AmiException("Requested functionality is not supported as %s "
+                           "is not implemented for this model!",
+                           __func__); // not implemented
+    }
+
+
     /**
      * Expands conservation law for states
      * @param x_rdata pointer to state variables with conservation laws
@@ -507,11 +533,11 @@ class Model : public AbstractModel {
     /**
      * Sensitivity of time-resolved measurement negative log-likelihood Jy
      * w.r.t. state variables
-     * @param dJydx pointer to vector with values of state derivative of Jy
+     * @param dJydx vector with values of state derivative of Jy
      * @param it timepoint index
      * @param edata pointer to experimental data instance
      */
-    void fdJydx(std::vector<realtype> *dJydx, const int it,
+    void fdJydx(std::vector<realtype> &dJydx, const int it,
                 const ExpData *edata);
 
     /**
@@ -1140,7 +1166,7 @@ class Model : public AbstractModel {
      * @return that
      */
     bool getAlwaysCheckFinite() const;
-    
+
     /**
      * @brief check whether the model was generated from python
      * @return that
@@ -1176,6 +1202,8 @@ class Model : public AbstractModel {
     const int ndwdp;
     /** number of nonzero entries in dxdotdw */
     const int ndxdotdw;
+    /** number of nonzero entries in dJydy */
+    std::vector<int> ndJydy;
     /** number of nonzero entries in jacobian */
     const int nnz;
     /** dimension of the augmented objective function for 2nd order ASA */
@@ -1194,29 +1222,37 @@ class Model : public AbstractModel {
 
     /** data standard deviation for current timepoint (dimension: ny) */
     std::vector<realtype> sigmay;
+
     /** parameter derivative of data standard deviation for current timepoint
      * (dimension: nplist x ny, row-major) */
     std::vector<realtype> dsigmaydp;
+
     /** event standard deviation for current timepoint (dimension: nz) */
     std::vector<realtype> sigmaz;
+
     /** parameter derivative of event standard deviation for current timepoint
      * (dimension: nplist x nz, row-major) */
     std::vector<realtype> dsigmazdp;
+
     /** parameter derivative of data likelihood for current timepoint
      * (dimension: nplist x nJ, row-major) */
     std::vector<realtype> dJydp;
+
     /** parameter derivative of event likelihood for current timepoint
      * (dimension: nplist x nJ, row-major) */
     std::vector<realtype> dJzdp;
 
     /** change in x at current timepoint (dimension: nx_solver) */
     std::vector<realtype> deltax;
+
     /** change in sx at current timepoint (dimension: nplist x nx_solver,
      * row-major) */
     std::vector<realtype> deltasx;
+
     /** change in xB at current timepoint (dimension: nJ x nxtrue_cl, row-major)
      */
     std::vector<realtype> deltaxB;
+
     /** change in qB at current timepoint (dimension: nJ x nplist, row-major) */
     std::vector<realtype> deltaqB;
 
@@ -1351,65 +1387,82 @@ class Model : public AbstractModel {
 
     /** Sparse Jacobian (dimension: nnz)*/
     SUNMatrixWrapper J;
-    
+
     /** Sparse dxdotdw temporary storage (dimension: ndxdotdw) */
     SUNMatrixWrapper dxdotdw;
-    
+
     /** Sparse dwdx temporary storage (dimension: ndwdx) */
     SUNMatrixWrapper dwdx;
-    
+
     /** Dense Mass matrix (dimension: nx_solver x nx_solver) */
     SUNMatrixWrapper M;
 
     /** current observable (dimension: nytrue) */
     std::vector<realtype> my;
+
     /** current event measurement (dimension: nztrue) */
     std::vector<realtype> mz;
+
+    /** Sparse observable derivative of data likelihood
+     * (dimension nytrue, nJ x ny, ordering = ?) */
+    std::vector<SUNMatrixWrapper> dJydy;
+
     /** observable derivative of data likelihood (dimension nJ x nytrue x ny,
-     * ordering = ?) */
-    std::vector<realtype> dJydy;
+      * ordering = ?) (only used if wasPythonGenerated()==false ) */
+    std::vector<realtype> dJydy_matlab;
+
     /** observable sigma derivative of data likelihood (dimension nJ x nytrue x
      * ny, ordering = ?) */
     std::vector<realtype> dJydsigma;
 
     /** event ouput derivative of event likelihood (dimension nJ x nztrue x nz,
      * ordering = ?) */
     std::vector<realtype> dJzdz;
+
     /** event sigma derivative of event likelihood (dimension nJ x nztrue x nz,
      * ordering = ?) */
     std::vector<realtype> dJzdsigma;
+
     /** event ouput derivative of event likelihood at final timepoint (dimension
      * nJ x nztrue x nz, ordering = ?) */
     std::vector<realtype> dJrzdz;
+
     /** event sigma derivative of event likelihood at final timepoint (dimension
      * nJ x nztrue x nz, ordering = ?) */
     std::vector<realtype> dJrzdsigma;
+
     /** state derivative of event output (dimension: nz x nx_solver, ordering =
      * ?) */
     std::vector<realtype> dzdx;
+
     /** parameter derivative of event output (dimension: nz x nplist, ordering =
      * ?) */
     std::vector<realtype> dzdp;
+
     /** state derivative of event timepoint (dimension: nz x nx_solver, ordering
      * = ?) */
     std::vector<realtype> drzdx;
+
     /** parameter derivative of event timepoint (dimension: nz x nplist,
      * ordering = ?) */
     std::vector<realtype> drzdp;
+
     /** parameter derivative of observable (dimension: nplist x ny, row-major)
      */
     std::vector<realtype> dydp;
 
-    /** state derivative of observable (dimension: ny x nx_solver, ordering = ?)
-     */
+    /** state derivative of observable
+     * (dimension: nx_solver x ny, ordering = row-major) */
     std::vector<realtype> dydx;
+
     /** tempory storage of w data across functions (dimension: nw) */
     std::vector<realtype> w;
+
     /** tempory storage of sparse/dense dwdp data across functions
      (dimension: ndwdp)
      */
     std::vector<realtype> dwdp;
-    
+
     /** tempory storage of stau data across functions (dimension: nplist) */
     std::vector<realtype> stau;
 
@@ -1419,6 +1472,7 @@ class Model : public AbstractModel {
 
     /** tempory storage x_rdata (dimension: nx_rdata) */
     std::vector<realtype> x_rdata;
+
     /** tempory storage sx_rdata slice (dimension: nx_rdata) */
     std::vector<realtype> sx_rdata;