diff --git a/.gitignore b/.gitignore index 39f05d4405..0eab397543 100644 --- a/.gitignore +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 512d2f2b52..336615b52e 100644 --- a/CMakeLists.txt +++ b/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. diff --git a/INSTALL.md b/INSTALL.md index d19d8a2436..7299e8a08e 100755 --- a/INSTALL.md +++ b/INSTALL.md @@ -7,7 +7,9 @@ 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/) @@ -15,10 +17,6 @@ The GIT repository can currently be found at https://github.com/ICB-DCM/AMICI an 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,7 +165,9 @@ 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 @@ -175,6 +175,14 @@ 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. diff --git a/documentation/FAQ.md b/documentation/FAQ.md index f114402b26..b84f6080d5 100644 --- a/documentation/FAQ.md +++ b/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! diff --git a/include/amici/defines.h b/include/amici/defines.h index 65732c9d23..15fe280acc 100644 --- a/include/amici/defines.h +++ b/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, diff --git a/include/amici/model.h b/include/amici/model.h index 20b6e1ccbe..9e2169761c 100644 --- a/include/amici/model.h +++ b/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 ndJydy, const int nnz, const int ubw, const int lbw, amici::SecondOrderMode o2mode, const std::vector &p, std::vector k, const std::vector &plist, std::vector 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 *dJydx, const int it, + void fdJydx(std::vector &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 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 sigmay; + /** parameter derivative of data standard deviation for current timepoint * (dimension: nplist x ny, row-major) */ std::vector dsigmaydp; + /** event standard deviation for current timepoint (dimension: nz) */ std::vector sigmaz; + /** parameter derivative of event standard deviation for current timepoint * (dimension: nplist x nz, row-major) */ std::vector dsigmazdp; + /** parameter derivative of data likelihood for current timepoint * (dimension: nplist x nJ, row-major) */ std::vector dJydp; + /** parameter derivative of event likelihood for current timepoint * (dimension: nplist x nJ, row-major) */ std::vector dJzdp; /** change in x at current timepoint (dimension: nx_solver) */ std::vector deltax; + /** change in sx at current timepoint (dimension: nplist x nx_solver, * row-major) */ std::vector deltasx; + /** change in xB at current timepoint (dimension: nJ x nxtrue_cl, row-major) */ std::vector deltaxB; + /** change in qB at current timepoint (dimension: nJ x nplist, row-major) */ std::vector deltaqB; @@ -1351,23 +1387,30 @@ 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 my; + /** current event measurement (dimension: nztrue) */ std::vector mz; + + /** Sparse observable derivative of data likelihood + * (dimension nytrue, nJ x ny, ordering = ?) */ + std::vector dJydy; + /** observable derivative of data likelihood (dimension nJ x nytrue x ny, - * ordering = ?) */ - std::vector dJydy; + * ordering = ?) (only used if wasPythonGenerated()==false ) */ + std::vector dJydy_matlab; + /** observable sigma derivative of data likelihood (dimension nJ x nytrue x * ny, ordering = ?) */ std::vector dJydsigma; @@ -1375,41 +1418,51 @@ class Model : public AbstractModel { /** event ouput derivative of event likelihood (dimension nJ x nztrue x nz, * ordering = ?) */ std::vector dJzdz; + /** event sigma derivative of event likelihood (dimension nJ x nztrue x nz, * ordering = ?) */ std::vector dJzdsigma; + /** event ouput derivative of event likelihood at final timepoint (dimension * nJ x nztrue x nz, ordering = ?) */ std::vector dJrzdz; + /** event sigma derivative of event likelihood at final timepoint (dimension * nJ x nztrue x nz, ordering = ?) */ std::vector dJrzdsigma; + /** state derivative of event output (dimension: nz x nx_solver, ordering = * ?) */ std::vector dzdx; + /** parameter derivative of event output (dimension: nz x nplist, ordering = * ?) */ std::vector dzdp; + /** state derivative of event timepoint (dimension: nz x nx_solver, ordering * = ?) */ std::vector drzdx; + /** parameter derivative of event timepoint (dimension: nz x nplist, * ordering = ?) */ std::vector drzdp; + /** parameter derivative of observable (dimension: nplist x ny, row-major) */ std::vector 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 dydx; + /** tempory storage of w data across functions (dimension: nw) */ std::vector w; + /** tempory storage of sparse/dense dwdp data across functions (dimension: ndwdp) */ std::vector dwdp; - + /** tempory storage of stau data across functions (dimension: nplist) */ std::vector stau; @@ -1419,6 +1472,7 @@ class Model : public AbstractModel { /** tempory storage x_rdata (dimension: nx_rdata) */ std::vector x_rdata; + /** tempory storage sx_rdata slice (dimension: nx_rdata) */ std::vector sx_rdata; diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 74130e945c..177eab99f2 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -46,6 +46,7 @@ namespace amici { * @param ndwdp number of nonzero elements in the p derivative of the * repeating elements * @param ndxdotdw number of nonzero elements dxdotdw + * @param ndJydy number of nonzero elements dJydy * @param nnz number of nonzero elements in Jacobian * @param ubw upper matrix bandwidth in the Jacobian * @param lbw lower matrix bandwidth in the Jacobian @@ -60,14 +61,16 @@ namespace amici { 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 ndwdp, const int ndxdotdw, std::vector ndJydy, + const int nnz, const int ubw, const int lbw, const SecondOrderMode o2mode, std::vector const &p, std::vector const &k, std::vector const &plist, std::vector const &idlist, std::vector const &z2event) : Model(nx_rdata, nxtrue_rdata, nx_solver, nxtrue_solver, ny, - nytrue, nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, nnz, + nytrue, nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, + ndJydy, nnz, ubw, lbw, o2mode, p, k, plist, idlist, z2event) {} virtual void fJ(realtype t, realtype cj, AmiVector *x, AmiVector *dx, diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 34f44a2af3..798e948a27 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -46,6 +46,7 @@ namespace amici { * @param ndwdp number of nonzero elements in the p derivative of the * repeating elements * @param ndxdotdw number of nonzero elements dxdotdw + * @param ndJydy number of nonzero elements dJydy * @param nnz number of nonzero elements in Jacobian * @param ubw upper matrix bandwidth in the Jacobian * @param lbw lower matrix bandwidth in the Jacobian @@ -60,14 +61,16 @@ namespace amici { 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 ndwdp, const int ndxdotdw, std::vector ndJydy, + const int nnz, const int ubw, const int lbw, const SecondOrderMode o2mode, std::vector const &p, std::vector const &k, std::vector const &plist, std::vector const &idlist, std::vector const &z2event) : Model(nx_rdata, nxtrue_rdata, nx_solver, nxtrue_solver, ny, - nytrue, nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, nnz, + nytrue, nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, + ndJydy, nnz, ubw, lbw, o2mode, p, k, plist, idlist, z2event) {} virtual void fJ(realtype t, realtype cj, AmiVector *x, AmiVector *dx, diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h index a632f5ffe6..efc5a1d6bc 100644 --- a/include/amici/newton_solver.h +++ b/include/amici/newton_solver.h @@ -4,8 +4,7 @@ #include "amici/vector.h" #include "amici/defines.h" #include "amici/sundials_matrix_wrapper.h" - -#include "sundials/sundials_linearsolver.h" // SUNLinearSolver +#include "amici/sundials_linsol_wrapper.h" #include @@ -25,13 +24,49 @@ class AmiVector; class NewtonSolver { public: + /** + * Initializes all members with the provided objects + * + * @param t pointer to time variable + * @param x pointer to state variables + * @param model pointer to the AMICI model object + * @param rdata pointer to the return data object + */ NewtonSolver(realtype *t, AmiVector *x, Model *model, ReturnData *rdata); + /** + * Factory method to create a NewtonSolver based on linsolType + * + * @param t pointer to time variable + * @param x pointer to state variables + * @param linsolType integer indicating which linear solver to use + * @param model pointer to the AMICI model object + * @param rdata pointer to the return data object + * @param maxlinsteps maximum number of allowed linear steps per Newton step for steady state computation + * @param maxsteps maximum number of allowed Newton steps for steady state computation + * @param atol absolute tolerance + * @param rtol relative tolerance + * @return solver NewtonSolver according to the specified linsolType + */ static std::unique_ptr getSolver(realtype *t, AmiVector *x, LinearSolver linsolType, Model *model, ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol); + /** + * Computes the solution of one Newton iteration + * + * @param ntry integer newton_try integer start number of Newton solver + * (1 or 2) + * @param nnewt integer number of current Newton step + * @param delta containing the RHS of the linear system, will be + * overwritten by solution to the linear system + */ void getStep(int ntry, int nnewt, AmiVector *delta); + /** + * Computes steady state sensitivities + * + * @param sx pointer to state variable sensitivities + */ void computeNewtonSensis(AmiVectorArray *sx); /** @@ -88,16 +123,41 @@ class NewtonSolver { class NewtonSolverDense : public NewtonSolver { public: + /** + * Constructor, initializes all members with the provided objects + * and initializes temporary storage objects + * + * @param t pointer to time variable + * @param x pointer to state variables + * @param model pointer to the AMICI model object + * @param rdata pointer to the return data object + */ + NewtonSolverDense(realtype *t, AmiVector *x, Model *model, ReturnData *rdata); ~NewtonSolverDense() override; + /** + * Solves the linear system for the Newton step + * + * @param rhs containing the RHS of the linear system, will be + * overwritten by solution to the linear system + */ void solveLinearSystem(AmiVector *rhs) override; + + /** + * Writes the Jacobian for the Newton iteration and passes it to the linear + * solver + * + * @param ntry integer newton_try integer start number of Newton solver + * (1 or 2) + * @param nnewt integer number of current Newton step + */ void prepareLinearSystem(int ntry, int nnewt) override; private: /** temporary storage of Jacobian */ SUNMatrixWrapper Jtmp; - + /** dense linear solver */ SUNLinearSolver linsol = nullptr; }; @@ -110,16 +170,40 @@ class NewtonSolverDense : public NewtonSolver { class NewtonSolverSparse : public NewtonSolver { public: + /** + * Constructor, initializes all members with the provided objects, + * initializes temporary storage objects and the klu solver + * + * @param t pointer to time variable + * @param x pointer to state variables + * @param model pointer to the AMICI model object + * @param rdata pointer to the return data object + */ NewtonSolverSparse(realtype *t, AmiVector *x, Model *model, ReturnData *rdata); ~NewtonSolverSparse() override; + /** + * Solves the linear system for the Newton step + * + * @param rhs containing the RHS of the linear system,will be + * overwritten by solution to the linear system + */ void solveLinearSystem(AmiVector *rhs) override; + + /** + * Writes the Jacobian for the Newton iteration and passes it to the linear + * solver + * + * @param ntry integer newton_try integer start number of Newton solver + * (1 or 2) + * @param nnewt integer number of current Newton step + */ void prepareLinearSystem(int ntry, int nnewt) override; private: /** temporary storage of Jacobian */ SUNMatrixWrapper Jtmp; - + /** sparse linear solver */ SUNLinearSolver linsol = nullptr; }; @@ -132,11 +216,46 @@ class NewtonSolverSparse : public NewtonSolver { class NewtonSolverIterative : public NewtonSolver { public: + /** + * Constructor, initializes all members with the provided objects + * @param t pointer to time variable + * @param x pointer to state variables + * @param model pointer to the AMICI model object + * @param rdata pointer to the return data object + */ NewtonSolverIterative(realtype *t, AmiVector *x, Model *model, ReturnData *rdata); virtual ~NewtonSolverIterative() = default; + /** + * Solves the linear system for the Newton step by passing it to + * linsolveSPBCG + * + * @param rhs containing the RHS of the linear system, will be + * overwritten by solution to the linear system + */ void solveLinearSystem(AmiVector *rhs) override; + + /** + * Writes the Jacobian for the Newton iteration and passes it to the linear + * solver. + * Also wraps around getSensis for iterative linear solver. + * + * @param ntry integer newton_try integer start number of Newton solver + * (1 or 2) + * @param nnewt integer number of current Newton step + */ void prepareLinearSystem(int ntry, int nnewt) override; + + /** + * Iterative linear solver created from SPILS BiCG-Stab. + * Solves the linear system within each Newton step if iterative solver is + * chosen. + * + * @param ntry integer newton_try integer start number of Newton solver + * (1 or 2) + * @param nnewt integer number of current Newton step + * @param ns_delta Newton step + */ void linsolveSPBCG(int ntry, int nnewt, AmiVector *ns_delta); private: diff --git a/include/amici/solver.h b/include/amici/solver.h index 0d9da0d8e8..b5b9b89ac1 100644 --- a/include/amici/solver.h +++ b/include/amici/solver.h @@ -555,17 +555,23 @@ class Solver { void setInterpolationType(InterpolationType interpType); /** - * @brief sets KLU state ordering mode - * @return + * @brief Gets KLU / SuperLUMT state ordering mode + * + * @return State-ordering as integer according to + * SUNLinSolKLU::StateOrdering or SUNLinSolSuperLUMT::StateOrdering + * (which differ). */ - StateOrdering getStateOrdering() const; + int getStateOrdering() const; /** - * @brief sets KLU state ordering mode (only applies when linsol is set to - * amici.AMICI_KLU) + * @brief Sets KLU / SuperLUMT state ordering mode + * + * This only applies when linsol is set to LinearSolver::KLU or + * LinearSolver::SuperLUMT. Mind the difference between + * SUNLinSolKLU::StateOrdering and SUNLinSolSuperLUMT::StateOrdering. * @param ordering state ordering */ - void setStateOrdering(StateOrdering ordering); + void setStateOrdering(int ordering); /** * @brief returns stability limit detection mode @@ -1160,7 +1166,7 @@ class Solver { booleantype stldet = true; /** state ordering */ - StateOrdering ordering = StateOrdering::AMD; + int ordering = static_cast(SUNLinSolKLU::StateOrdering::AMD); /** maximum number of allowed Newton steps for steady state computation */ int newton_maxsteps = 0; diff --git a/include/amici/sundials_linsol_wrapper.h b/include/amici/sundials_linsol_wrapper.h index b4e12f7eac..e1351d9806 100644 --- a/include/amici/sundials_linsol_wrapper.h +++ b/include/amici/sundials_linsol_wrapper.h @@ -5,6 +5,7 @@ #include "amici/sundials_matrix_wrapper.h" #include "amici/vector.h" +#include #include #include #include @@ -13,7 +14,9 @@ #include #include #include - +#ifdef SUNDIALS_SUPERLUMT +#include +#endif #include #include @@ -187,6 +190,13 @@ class SUNLinSolDense : public SUNLinSolWrapper { */ class SUNLinSolKLU : public SUNLinSolWrapper { public: + /** KLU state reordering (different from SuperLUMT ordering!) */ + enum class StateOrdering { + AMD, + COLAMD, + natural + }; + /** * @brief Create KLU solver with given matrix * @param x A template for cloning vectors needed within the solver. @@ -228,6 +238,68 @@ class SUNLinSolKLU : public SUNLinSolWrapper { SUNMatrixWrapper A; }; +#ifdef SUNDIALS_SUPERLUMT +/** + * @brief SUNDIALS SuperLUMT sparse direct solver. + */ +class SUNLinSolSuperLUMT : public SUNLinSolWrapper { + public: + /** SuperLUMT ordering (different from KLU ordering!) */ + enum class StateOrdering { + natural, + minDegATA, + minDegATPlusA, + COLAMD, + }; + + /** + * @brief Create SuperLUMT solver with given matrix + * @param x A template for cloning vectors needed within the solver. + * @param A sparse matrix + * @param numThreads Number of threads to be used by SuperLUMT + */ + SUNLinSolSuperLUMT(N_Vector x, SUNMatrix A, int numThreads); + + /** + * @brief Create SuperLUMT solver and matrix to operate on + * + * Will set number of threads according to environment variable + * AMICI_SUPERLUMT_NUM_THREADS. Will default to 1 thread if unset. + * + * @param x A template for cloning vectors needed within the solver. + * @param nnz Number of non-zeros in matrix A + * @param sparsetype Sparse matrix type (CSC_MAT, CSR_MAT) + * @param ordering + */ + SUNLinSolSuperLUMT(AmiVector const &x, int nnz, int sparsetype, + StateOrdering ordering); + + /** + * @brief Create SuperLUMT solver and matrix to operate on + * @param x A template for cloning vectors needed within the solver. + * @param nnz Number of non-zeros in matrix A + * @param sparsetype Sparse matrix type (CSC_MAT, CSR_MAT) + * @param ordering + * @param numThreads Number of threads to be used by SuperLUMT + */ + SUNLinSolSuperLUMT(AmiVector const &x, int nnz, int sparsetype, + StateOrdering ordering, int numThreads); + + SUNMatrix getMatrix() const override; + + /** + * @brief Sets the ordering used by SuperLUMT for reducing fill in the + * linear solve. + * @param ordering + */ + void setOrdering(StateOrdering ordering); + + private: + /** Sparse matrix A for solver, only if created by here. */ + SUNMatrixWrapper A; +}; + +#endif /** * @brief SUNDIALS scaled preconditioned CG (Conjugate Gradient method) (PCG) @@ -308,7 +380,8 @@ class SUNLinSolSPBCGS : public SUNLinSolWrapper { * PREC_BOTH) * @param maxl Maximum number of solver iterations */ - SUNLinSolSPBCGS(N_Vector x, int pretype, int maxl); + SUNLinSolSPBCGS(N_Vector x, int pretype = PREC_NONE, + int maxl = SUNSPBCGS_MAXL_DEFAULT); /** * @brief SUNLinSolSPBCGS @@ -317,7 +390,8 @@ class SUNLinSolSPBCGS : public SUNLinSolWrapper { * PREC_BOTH) * @param maxl Maximum number of solver iterations */ - SUNLinSolSPBCGS(AmiVector const &x, int pretype, int maxl); + SUNLinSolSPBCGS(AmiVector const &x, int pretype = PREC_NONE, + int maxl = SUNSPBCGS_MAXL_DEFAULT); /** * @brief Sets the function pointer for ATimes @@ -449,7 +523,8 @@ class SUNLinSolSPGMR : public SUNLinSolWrapper { * PREC_BOTH) * @param maxl Maximum number of solver iterations */ - SUNLinSolSPGMR(AmiVector const &x, int pretype, int maxl); + SUNLinSolSPGMR(AmiVector const &x, int pretype = PREC_NONE, + int maxl = SUNSPGMR_MAXL_DEFAULT); /** * @brief Sets the function pointer for ATimes @@ -515,7 +590,8 @@ class SUNLinSolSPTFQMR : public SUNLinSolWrapper { * PREC_BOTH) * @param maxl Maximum number of solver iterations */ - SUNLinSolSPTFQMR(N_Vector x, int pretype, int maxl); + SUNLinSolSPTFQMR(N_Vector x, int pretype = PREC_NONE, + int maxl = SUNSPTFQMR_MAXL_DEFAULT); /** * @brief Create SPTFQMR solver @@ -524,7 +600,8 @@ class SUNLinSolSPTFQMR : public SUNLinSolWrapper { * PREC_BOTH) * @param maxl Maximum number of solver iterations */ - SUNLinSolSPTFQMR(AmiVector const &x, int pretype, int maxl); + SUNLinSolSPTFQMR(AmiVector const &x, int pretype = PREC_NONE, + int maxl = SUNSPTFQMR_MAXL_DEFAULT); /** * @brief Sets the function pointer for ATimes diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h index 18561d6286..0d773f0944 100644 --- a/include/amici/sundials_matrix_wrapper.h +++ b/include/amici/sundials_matrix_wrapper.h @@ -7,6 +7,8 @@ #include +#include + #include namespace amici { @@ -131,19 +133,12 @@ class SUNMatrixWrapper { * @return index array */ int sparsetype() const; - + /** * @brief reset data to zeroes */ void reset(); - /** - * @brief std::vector interface for multiply - * @param c output vector, may already contain values - * @param b multiplication vector - */ - void multiply(std::vector &c, const std::vector &b) const; - /** * @brief N_Vector interface for multiply * @param c output vector, may already contain values @@ -156,11 +151,16 @@ class SUNMatrixWrapper { * @param c output vector, may already contain values * @param b multiplication vector */ - void multiply(realtype *c, const realtype *b) const; + void multiply(gsl::span c, gsl::span b) const; + + /** + * @brief Set to 0.0 + */ + void zero(); private: void update_ptrs(); - + SUNMatrix matrix = nullptr; realtype *data_ptr = nullptr; sunindextype *indexptrs_ptr = nullptr; diff --git a/matlab/@amimodel/generateC.m b/matlab/@amimodel/generateC.m index 8d1d8be1ba..f40f869ac4 100644 --- a/matlab/@amimodel/generateC.m +++ b/matlab/@amimodel/generateC.m @@ -162,6 +162,7 @@ function generateC(this) fprintf(fid,[' ' num2str(this.ndwdx) ',\n']); fprintf(fid,[' ' num2str(this.ndwdp) ',\n']); fprintf(fid,[' 0,\n']); +fprintf(fid,[' {},\n']); fprintf(fid,[' ' num2str(this.nnz) ',\n']); fprintf(fid,[' ' num2str(this.ubw) ',\n']); fprintf(fid,[' ' num2str(this.lbw) ',\n']); diff --git a/matlab/examples/amiExamples.m b/matlab/examples/amiExamples.m index 9fe21d3d90..9dd4b5bb3e 100644 --- a/matlab/examples/amiExamples.m +++ b/matlab/examples/amiExamples.m @@ -13,4 +13,5 @@ example_jakstat_adjoint_hvp example_robertson example_neuron +example_dae_events path(oldpath); \ No newline at end of file diff --git a/matlab/examples/example_calvetti/example_calvetti.m b/matlab/examples/example_calvetti/example_calvetti.m new file mode 100755 index 0000000000..334f977be0 --- /dev/null +++ b/matlab/examples/example_calvetti/example_calvetti.m @@ -0,0 +1,101 @@ +function example_calvetti() +%% +% COMPILATION + +[exdir,~,~]=fileparts(which('example_calvetti.m')); +% compile the model +amiwrap('model_calvetti','model_calvetti_syms',exdir) + +% time vector +t = linspace(0,20,201); +p = []; +k = [0.29 0.74 0.44 0.08 0.27 0.18]; + +options = amioption('sensi',0,... + 'maxsteps',1e4); +D = amidata(length(t),6,0,2,4); +% load mex into memory +[~] = which('simulate_model_calvetti'); % fix for inaccessability problems +sol = simulate_model_calvetti(t,p,k,D,options); + +tic +sol = simulate_model_calvetti(t,p,k,D,options); +disp(['Time elapsed with cvodes: ' num2str(toc) ]) + +%% +% ODE15S + +y0 = [k(1); k(3); k(5); 1; 1; 1;]; +M = [1 0 0 0 0 0 + 0 1 0 0 0 0 + 0 0 1 0 0 0 + 0 0 0 0 0 0 + 0 0 0 0 0 0 + 0 0 0 0 0 0]; + +function [xdot] = dae_system(t,x,p,k,it) + if it<3 + h0 = 0; + else + h0 = 1; + end + if it<4 + h1 = 0; + else + h1 = 1; + end + xdot = [ + h0/31 - h1/31 - (100*x(1))/(899*k(1)) + (2*x(1)*(k(2)/2 - 1))/(31*k(1)*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) + 129/899; + (2*x(2)*(k(2) + k(4)/2 - 1))/(163*k(3)*(x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - (100*x(2))/(8313*k(3)) + 151/8313; + (x(3)^3*(k(2) + k(4) +k(6)/2 - 1))/(61*k(6)*k(5)^3*x(6)) - x(3)/(121999878*k(5)) + 500000/60999939; + h1/31 - h0/31 - x(4) + (100*x(1))/(899*k(1)) + (2*x(1)^2)/(k(2)*k(1)^2) - (x(1)^2*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2))/(k(2)*k(1)^2) - (2*x(1)*(k(2)/2 - 1))/(31*k(1)*(x(4)*((k(2)*k(1)^2)/x(1)^2 + (k(4)*k(3)^2)/x(2)^2) + x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - 129/899; + x(4) - x(5) + (100*x(2))/(8313*k(3)) - (2*x(2)*(k(2) + k(4)/2 - 1))/(163*k(3)*(x(5)*((k(4)*k(3)^2)/x(2)^2 + (k(6)*k(5)^2)/x(3)^2) + (k(6)*k(5)^2*x(6))/x(3)^2)) - 151/8313; + x(5) - x(6) + x(3)/(121999878*k(5)) - (x(3)^3*(k(2) + k(4) +k(6)/2 - 1))/(61*k(6)*k(5)^3*x(6)) - 500000/60999939; +]; +end + +options_ode15s = odeset('Mass',M, ... + 'RelTol',options.rtol, ... + 'AbsTol',options.atol); + +tic +X_ode15s = []; +tt = t; +tstop = [0,10,12,20]; +for it = 2:4 + [~, y] = ode15s(@(t,x) dae_system(t,x,p,k,it),t(t>=tstop(it-1) & t<=tstop(it)),y0,options_ode15s); + X_ode15s = [X_ode15s;y(1:end-1,:)]; + y0 = y(end,:); +end +X_ode15s = [X_ode15s;y(end,:)]; +disp(['Time elapsed with ode15s: ' num2str(toc) ]) + +%% +% PLOTTING +if(usejava('jvm')) + figure + c_x = get(gca,'ColorOrder'); + subplot(2,1,1) + for ix = 1:size(sol.x,2) + plot(t,sol.x(:,ix),'.-','Color',c_x(ix,:)) + hold on + plot(t,X_ode15s(:,ix),'d','Color',c_x(ix,:)) + end + legend('x1','x1_{ode15s}','x2','x2_{ode15s}', ... + 'x3','x3_{ode15s}','x4','x4_{ode15s}', ... + 'x5','x5_{ode15s}','x6','x6_{ode15s}', ... + 'Location','NorthEastOutside') + legend boxoff + xlabel('time t') + ylabel('x') + box on + subplot(2,1,2) + plot(t,abs(sol.x-X_ode15s),'--') + set(gca,'YScale','log') + legend('error x1','error x2','error x3','error x4','error x5','error x6','Location','NorthEastOutside') + legend boxoff + ylabel('x') + + set(gcf,'Position',[100 300 1200 500]) +end +end diff --git a/matlab/examples/example_calvetti/model_calvetti_syms.m b/matlab/examples/example_calvetti/model_calvetti_syms.m new file mode 100755 index 0000000000..5fc3661ebf --- /dev/null +++ b/matlab/examples/example_calvetti/model_calvetti_syms.m @@ -0,0 +1,75 @@ +function [model] = model_calvetti_syms() + +% set the parametrisation of the problem options are 'log', 'log10' and + +model.param = 'lin'; + +%% +% STATES +% create state syms +syms V1 V2 V3 f1 f2 f3 +% create state vector +model.sym.x = [V1, V2, V3, f1, f2, f3]; +%% +% PARAMETERS ( for these sensitivities will be computed ) + +% create parameter syms +% create parameter vector +model.sym.p = [ ]; +%% +% CONSTANTS ( for these no sensitivities will be computed ) +% this part is optional and can be ommited + +% create parameter syms +syms V1ss R1ss V2ss R2ss V3ss R3ss +% create parameter vector +model.sym.k = [V1ss, R1ss, V2ss, R2ss, V3ss, R3ss]; +%% +% SYSTEM EQUATIONS +% create symbolic variable for time +syms t f0 +model.sym.xdot = sym(zeros(size(model.sym.x))); +p1=1; +p2=1-R1ss; +p3=1-(R1ss+R2ss); +L1=(R1ss*(V1ss^2))^(1/3); +C1ss=V1ss/(p1-0.5*R1ss); + +C2ss=V2ss/(p2-0.5*R2ss); +L2=(R2ss*(V2ss^2))^(1/3); +C3ss=V3ss/(p3-0.5*R3ss); +L3=(R3ss*(V3ss^2))^(1/3); +R2=(L2^3)/(V2^2); +R1=(L1^3)/(V1^2); +R3=(L3^3)/(V3^2); +s=am_stepfun(t,10,1,12,0); +model.sym.xdot(1)=(1/31)*(((1.29 -(V1/V1ss))/(1.29-1))+s-2*(V1/(C1ss*((R1+R2)*f1+(R2+R3)*f2+R3*f3)))); +model.sym.xdot(2)=(1/163)*(((1.51 -(V2/V2ss))/(1.51-1))- 2*(V2/(C2ss*((R2+R3)*f2+R3*f3)))); +model.sym.xdot(3)=(1/122)*(((1000000 -(V3/V3ss))/(1000000-1))- 2*(V3/(C3ss*(R3*f3)))); +f0=(2/R1)-((R1+R2)*f1 + (R2+R3)*f2 +R3*f3)/R1; +model.sym.xdot(4)=f0-model.sym.xdot(1)-f1; +model.sym.xdot(5)=f1-model.sym.xdot(2)-f2; +model.sym.xdot(6)=f2-model.sym.xdot(3)-f3; + + +model.sym.M=[1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0;]; +%% +% INITIAL CONDITIONS +model.sym.x0(1) = V1ss; +model.sym.x0(2)=V2ss; +model.sym.x0(3)=V3ss; +model.sym.x0(4)=1; +model.sym.x0(5)=1; +model.sym.x0(6)=1; +model.sym.dx0 = sym(zeros(size(model.sym.x))); + +% OBSERVALES +model.sym.y = sym(zeros(1,6)); +model.sym.y(1) =V1; +model.sym.y(2)=V2; +model.sym.y(3)=V3; +model.sym.y(4)=f0; +model.sym.y(5)=f1; +model.sym.y(6)=f2; +end + diff --git a/models/model_calvetti/CMakeLists.txt b/models/model_calvetti/CMakeLists.txt new file mode 100644 index 0000000000..23760acd3c --- /dev/null +++ b/models/model_calvetti/CMakeLists.txt @@ -0,0 +1,94 @@ +cmake_minimum_required(VERSION 2.8) + +if(POLICY CMP0060) + cmake_policy(SET CMP0060 NEW) +endif(POLICY CMP0060) +if(POLICY CMP0065) + cmake_policy(SET CMP0065 NEW) +endif(POLICY CMP0065) + +project(model_calvetti) + +set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + +include(CheckCXXCompilerFlag) +set(MY_CXX_FLAGS -Wall -Wno-unused-function -Wno-unused-variable -Wno-unused-but-set-variable) +foreach(FLAG ${MY_CXX_FLAGS}) + unset(CUR_FLAG_SUPPORTED CACHE) + CHECK_CXX_COMPILER_FLAG(${FLAG} CUR_FLAG_SUPPORTED) + if(${CUR_FLAG_SUPPORTED}) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG}") + endif() +endforeach(FLAG) + +find_package(Amici HINTS ${CMAKE_CURRENT_LIST_DIR}/../../build) + +set(MODEL_DIR ${CMAKE_CURRENT_LIST_DIR}) + +set(SRC_LIST_LIB ${MODEL_DIR}/model_calvetti_J.cpp +${MODEL_DIR}/model_calvetti_JB.cpp +${MODEL_DIR}/model_calvetti_JDiag.cpp +${MODEL_DIR}/model_calvetti_JSparse.cpp +${MODEL_DIR}/model_calvetti_JSparseB.cpp +${MODEL_DIR}/model_calvetti_Jy.cpp +${MODEL_DIR}/model_calvetti_M.cpp +${MODEL_DIR}/model_calvetti_dJydsigma.cpp +${MODEL_DIR}/model_calvetti_dJydy.cpp +${MODEL_DIR}/model_calvetti_dwdx.cpp +${MODEL_DIR}/model_calvetti_dydx.cpp +${MODEL_DIR}/model_calvetti_root.cpp +${MODEL_DIR}/model_calvetti_sigmay.cpp +${MODEL_DIR}/model_calvetti_w.cpp +${MODEL_DIR}/model_calvetti_x0.cpp +${MODEL_DIR}/model_calvetti_xdot.cpp +${MODEL_DIR}/model_calvetti_y.cpp + +${MODEL_DIR}/wrapfunctions.cpp +) + +add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) +add_library(model ALIAS ${PROJECT_NAME}) + +target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") + +target_link_libraries(${PROJECT_NAME} + PUBLIC Upstream::amici +) + +set(SRC_LIST_EXE main.cpp) + +add_executable(simulate_${PROJECT_NAME} ${SRC_LIST_EXE}) + +target_link_libraries(simulate_${PROJECT_NAME} ${PROJECT_NAME}) + +if($ENV{ENABLE_GCOV_COVERAGE}) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage") +endif() + +## SWIG +option(ENABLE_SWIG "Build swig/python library?" ON) +if(ENABLE_SWIG) + if(NOT(${CMAKE_VERSION} VERSION_LESS 3.8)) + add_subdirectory(swig) + else() + message(WARNING "Unable to build SWIG interface, upgrade CMake to >=3.8.") + endif() +endif() + + +# +include(GNUInstallDirs) +install(TARGETS ${PROJECT_NAME} EXPORT ${PROJECT_NAME}Targets + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} +) +export(EXPORT ${PROJECT_NAME}Targets FILE ${PROJECT_NAME}Config.cmake + NAMESPACE Upstream:: + ) +# + diff --git a/models/model_calvetti/main.cpp b/models/model_calvetti/main.cpp new file mode 100644 index 0000000000..6d6caf54a8 --- /dev/null +++ b/models/model_calvetti/main.cpp @@ -0,0 +1,115 @@ +#include +#include +#include +#include +#include +#include + +#include /* AMICI base functions */ +#include /* AMICI HDF5 I/O functions */ +#include "wrapfunctions.h" /* model-provided functions */ + +/* This is a scaffold for a stand-alone AMICI simulation executable + * demonstrating + * use of the AMICI C++ API. + * + * This program reads AMICI options from an HDF5 file, prints some results + * and writes additional results to an HDF5 file. The name of the HDF5 file + * is expected as single command line argument. + * + * An initial HDF5 file with the required fields can be generated using MATLAB + * by adding the following lines + * at the end of simulate_${MODEL_NAME}.m file just before the final "end": + * + * %% Write data that is passed to AMICI to HDF5 + * hdffile = fullfile(pwd, 'mydata.h5'); + * structToHDF5Attribute(hdffile, '/options', options_ami); + * h5writeatt(hdffile, '/options', 'ts', tout); + * h5writeatt(hdffile, '/options', 'nt', numel(tout)); + * h5writeatt(hdffile, '/options', 'theta', theta); + * h5writeatt(hdffile, '/options', 'kappa', kappa); + * if(~isempty(data)) + * structToHDF5Attribute(hdffile, '/data', data); + * end + * + * ... and then running a simulation from MATLAB as usual. + * + * Default UserData settings can be written to an HDF5 file with: + * structToHDF5Attribute('test.h5', '/options', amioption()) + */ + +// Function prototypes +void processReturnData(amici::ReturnData *rdata, amici::Model *model); +void printReturnData(amici::ReturnData *rdata, amici::Model *model); + +int main(int argc, char **argv) { + // HDF5 file to read and write data (full path) + const char *hdffile; + + // Check command line arguments + if (argc != 2) { + fprintf(stderr, "Error: must provide HDF5 input file as first and only " + "argument.\n"); + return 1; + } else { + hdffile = argv[1]; + } + + auto model = getModel(); + auto solver = model->getSolver(); + + // Read AMICI settings and model parameters from HDF5 file + amici::hdf5::readModelDataFromHDF5(hdffile, *model, "/options"); + amici::hdf5::readSolverSettingsFromHDF5(hdffile, *solver, "/options"); + + // Read ExpData (experimental data for model) from HDF5 file + auto edata = amici::hdf5::readSimulationExpData(hdffile, "/data", *model); + + // Run the simulation + auto rdata = runAmiciSimulation(*solver, edata.get(), *model); + + // Do something with the simulation results + processReturnData(rdata.get(), model.get()); + + // Save simulation results to HDF5 file + amici::hdf5::writeReturnData(*rdata, hdffile, "/solution"); + + return 0; +} + +void processReturnData(amici::ReturnData *rdata, amici::Model *model) { + // show some the simulation results + printReturnData(rdata, model); +} + +void printReturnData(amici::ReturnData *rdata, amici::Model *model) { + // Print of some the simulation results + + printf("\n\nStates (xdata):\n"); + for (int i = 0; i < model->nx_rdata; ++i) { + for (int j = 0; j < model->nt(); ++j) + printf("%e\t", rdata->x[j + model->nt() * i]); + printf("\n"); + } + + printf("\nObservables (ydata):\n"); + for (int i = 0; i < model->ny; ++i) { + for (int j = 0; j < model->nt(); ++j) + printf("%e\t", rdata->y[j + model->nt() * i]); + printf("\n"); + } + + printf("\n\ndx/dt (xdotdata):\n"); + for (int i = 0; i < model->nx_solver; ++i) + printf("%e\t", rdata->xdot[i]); + + // printf("\nJacobian (jdata)\n"); + // for(int i = 0; i < nx_solver; ++i) { + // for(int j = 0; j < nx_solver; ++j) + // printf("%e\t", rdata->J[i + nx_solver * j]); + // printf("\n"); + // } + + printf("\n"); + printf("Loglikelihood (llh): %e\n", rdata->llh); +} diff --git a/models/model_calvetti/model_calvetti.h b/models/model_calvetti/model_calvetti.h new file mode 100644 index 0000000000..97d201b047 --- /dev/null +++ b/models/model_calvetti/model_calvetti.h @@ -0,0 +1,214 @@ +#ifndef _amici_model_calvetti_h +#define _amici_model_calvetti_h +/* Generated by amiwrap (R2017b) 6a604c72d192e6f8726862b089ad2c0adb456a38 */ +#include +#include +#include "amici/defines.h" +#include //SUNMatrixContent_Sparse definition +#include "amici/solver_idas.h" +#include "amici/model_dae.h" + +namespace amici { +class Solver; +} + + +extern void J_model_calvetti(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx); +extern void JB_model_calvetti(realtype *JB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx); +extern void JDiag_model_calvetti(realtype *JDiag, 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); +extern void JSparse_model_calvetti(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); +extern void JSparseB_model_calvetti(SUNMatrixContent_Sparse JSparseB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx); +extern void Jy_model_calvetti(double *nllh, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); +extern void M_model_calvetti(realtype *M, const realtype t, const realtype *x, const realtype *p, const realtype *k); +extern void dJydsigma_model_calvetti(double *dJydsigma, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); +extern void dJydy_model_calvetti(double *dJydy, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my); +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 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); +extern void y_model_calvetti(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w); + +class Model_model_calvetti : public amici::Model_DAE { +public: + Model_model_calvetti() : amici::Model_DAE(6, + 6, + 6, + 6, + 6, + 6, + 0, + 0, + 4, + 1, + 38, + 53, + 0, + 0, + {}, + 26, + 5, + 3, + amici::SecondOrderMode::none, + std::vector(0,1.0), + std::vector(6,1.0), + std::vector(), + std::vector{1, 1, 1, 0, 0, 0}, + std::vector{}) + {}; + + virtual amici::Model* clone() const override { return new Model_model_calvetti(*this); }; + + const std::string getAmiciCommit() const override { return "6a604c72d192e6f8726862b089ad2c0adb456a38"; }; + + virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { + J_model_calvetti(J, t, x, p, k, h, cj, dx, w, dwdx); + } + + virtual void fJB(realtype *JB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx) override { + JB_model_calvetti(JB, t, x, p, k, h, cj, xB, dx, dxB, w, dwdx); + } + + virtual void fJDiag(realtype *JDiag, 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 { + JDiag_model_calvetti(JDiag, t, x, p, k, h, cj, dx, w, dwdx); + } + + virtual 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); + } + + virtual void fJSparseB(SUNMatrixContent_Sparse JSparseB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx) override { + JSparseB_model_calvetti(JSparseB, t, x, p, k, h, cj, xB, dx, dxB, w, dwdx); + } + + virtual void fJrz(double *nllh, const int iz, const realtype *p, const realtype *k, const double *rz, const double *sigmaz) override { + } + + virtual void fJy(double *nllh, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) override { + Jy_model_calvetti(nllh, iy, p, k, y, sigmay, my); + } + + virtual void fJz(double *nllh, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override { + } + + virtual void fM(realtype *M, const realtype t, const realtype *x, const realtype *p, const realtype *k) override { + M_model_calvetti(M, t, x, p, k); + } + + virtual void fdJrzdsigma(double *dJrzdsigma, const int iz, const realtype *p, const realtype *k, const double *rz, const double *sigmaz) override { + } + + virtual void fdJrzdz(double *dJrzdz, const int iz, const realtype *p, const realtype *k, const double *rz, const double *sigmaz) override { + } + + virtual void fdJydsigma(double *dJydsigma, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) override { + dJydsigma_model_calvetti(dJydsigma, iy, p, k, y, sigmay, my); + } + + virtual void fdJydy(double *dJydy, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) override { + dJydy_model_calvetti(dJydy, iy, p, k, y, sigmay, my); + } + + virtual void fdJzdsigma(double *dJzdsigma, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override { + } + + virtual void fdJzdz(double *dJzdz, const int iz, const realtype *p, const realtype *k, const double *z, const double *sigmaz, const double *mz) override { + } + + virtual void fdeltaqB(double *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *xB) override { + } + + virtual void fdeltasx(double *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) override { + } + + virtual void fdeltax(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) override { + } + + virtual void fdeltaxB(double *deltaxB, 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 *xB) override { + } + + virtual void fdrzdp(double *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override { + } + + virtual void fdrzdx(double *drzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { + } + + virtual void fdsigmaydp(double *dsigmaydp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + } + + virtual void fdsigmazdp(double *dsigmazdp, const realtype t, const realtype *p, const realtype *k, const int ip) override { + } + + 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) override { + } + + 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) override { + dwdx_model_calvetti(dwdx, t, x, p, k, h, w, tcl); + } + + virtual void fdxdotdp(realtype *dxdotdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *dx, const realtype *w, const realtype *dwdp) override { + } + + virtual void fdydp(double *dydp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip, const realtype *w, const realtype *dwdp) override { + } + + virtual void fdydx(double *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dwdx) override { + dydx_model_calvetti(dydx, t, x, p, k, h, w, dwdx); + } + + virtual void fdzdp(double *dzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override { + } + + virtual void fdzdx(double *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { + } + + virtual void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx) override { + root_model_calvetti(root, t, x, p, k, h, dx); + } + + virtual void frz(double *rz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { + } + + virtual void fsigmay(double *sigmay, const realtype t, const realtype *p, const realtype *k) override { + sigmay_model_calvetti(sigmay, t, p, k); + } + + virtual void fsigmaz(double *sigmaz, const realtype t, const realtype *p, const realtype *k) override { + } + + virtual void fsrz(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) override { + } + + virtual void fstau(double *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip, const int ie) override { + } + + virtual void fsx0(realtype *sx0, const realtype t,const realtype *x0, const realtype *p, const realtype *k, const int ip) override { + } + + virtual void fsz(double *sz, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *sx, const int ip) override { + } + + virtual void fw(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override { + w_model_calvetti(w, t, x, p, k, h, tcl); + } + + virtual void fx0(realtype *x0, const realtype t, const realtype *p, const realtype *k) override { + x0_model_calvetti(x0, t, p, k); + } + + virtual void fxdot(realtype *xdot, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *w) override { + xdot_model_calvetti(xdot, t, x, p, k, h, dx, w); + } + + virtual void fy(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w) override { + y_model_calvetti(y, t, x, p, k, h, w); + } + + virtual void fz(double *z, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override { + } + +}; + +#endif /* _amici_model_calvetti_h */ diff --git a/models/model_calvetti/model_calvetti_J.cpp b/models/model_calvetti/model_calvetti_J.cpp new file mode 100644 index 0000000000..0d3a72f80d --- /dev/null +++ b/models/model_calvetti/model_calvetti_J.cpp @@ -0,0 +1,37 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void J_model_calvetti(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) { + J[0+0*6] = -cj-w[0]*(1.0E2/8.99E2)+dwdx[6]; + J[0+1*6] = dwdx[16]; + J[0+2*6] = dwdx[28]; + J[0+3*6] = dwdx[36]; + J[0+4*6] = dwdx[40]; + J[0+5*6] = dwdx[47]; + J[1+1*6] = -cj-w[6]*1.202935161794779E-2+dwdx[19]; + J[1+2*6] = dwdx[31]; + J[1+4*6] = dwdx[43]; + J[1+5*6] = dwdx[50]; + J[2+2*6] = -cj-w[31]*8.196729508204918E-9+dwdx[32]; + J[2+5*6] = dwdx[52]; + J[3+0*6] = w[0]*(1.0E2/8.99E2)-dwdx[6]+dwdx[7]*(w[23]*w[24]*2.0-w[20]*w[23]*w[24])-w[23]*w[24]*w[25]*dwdx[4]; + J[3+1*6] = -dwdx[16]-w[23]*w[24]*w[25]*dwdx[14]; + J[3+2*6] = -dwdx[28]-w[23]*w[24]*w[25]*dwdx[26]; + J[3+3*6] = -dwdx[36]-w[23]*w[24]*w[25]*dwdx[34]-1.0; + J[3+4*6] = -dwdx[40]-w[23]*w[24]*w[25]*dwdx[38]; + J[3+5*6] = -dwdx[47]-w[23]*w[24]*w[25]*dwdx[45]; + J[4+1*6] = w[6]*1.202935161794779E-2-dwdx[19]; + J[4+2*6] = -dwdx[31]; + J[4+3*6] = 1.0; + J[4+4*6] = -dwdx[43]-1.0; + J[4+5*6] = -dwdx[50]; + J[5+2*6] = w[31]*8.196729508204918E-9-dwdx[32]; + J[5+4*6] = 1.0; + J[5+5*6] = -dwdx[52]-1.0; +} + diff --git a/models/model_calvetti/model_calvetti_JB.cpp b/models/model_calvetti/model_calvetti_JB.cpp new file mode 100644 index 0000000000..5642bc792e --- /dev/null +++ b/models/model_calvetti/model_calvetti_JB.cpp @@ -0,0 +1,37 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void JB_model_calvetti(realtype *JB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx) { + JB[0+0*6] = cj+w[0]*(1.0E2/8.99E2)-dwdx[6]; + JB[0+3*6] = w[0]*(-1.0E2/8.99E2)+dwdx[6]-dwdx[7]*(w[23]*w[24]*2.0-w[20]*w[23]*w[24])+w[23]*w[24]*w[25]*dwdx[4]; + JB[1+0*6] = -dwdx[16]; + JB[1+1*6] = cj+w[6]*1.202935161794779E-2-dwdx[19]; + JB[1+3*6] = dwdx[16]+w[23]*w[24]*w[25]*dwdx[14]; + JB[1+4*6] = w[6]*(-1.202935161794779E-2)+dwdx[19]; + JB[2+0*6] = -dwdx[28]; + JB[2+1*6] = -dwdx[31]; + JB[2+2*6] = cj+w[31]*8.196729508204918E-9-dwdx[32]; + JB[2+3*6] = dwdx[28]+w[23]*w[24]*w[25]*dwdx[26]; + JB[2+4*6] = dwdx[31]; + JB[2+5*6] = w[31]*(-8.196729508204918E-9)+dwdx[32]; + JB[3+0*6] = -dwdx[36]; + JB[3+3*6] = dwdx[36]+w[23]*w[24]*w[25]*dwdx[34]+1.0; + JB[3+4*6] = -1.0; + JB[4+0*6] = -dwdx[40]; + JB[4+1*6] = -dwdx[43]; + JB[4+3*6] = dwdx[40]+w[23]*w[24]*w[25]*dwdx[38]; + JB[4+4*6] = dwdx[43]+1.0; + JB[4+5*6] = -1.0; + JB[5+0*6] = -dwdx[47]; + JB[5+1*6] = -dwdx[50]; + JB[5+2*6] = -dwdx[52]; + JB[5+3*6] = dwdx[47]+w[23]*w[24]*w[25]*dwdx[45]; + JB[5+4*6] = dwdx[50]; + JB[5+5*6] = dwdx[52]+1.0; +} + diff --git a/models/model_calvetti/model_calvetti_JDiag.cpp b/models/model_calvetti/model_calvetti_JDiag.cpp new file mode 100644 index 0000000000..6a1cf9cd56 --- /dev/null +++ b/models/model_calvetti/model_calvetti_JDiag.cpp @@ -0,0 +1,17 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void JDiag_model_calvetti(realtype *JDiag, 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) { + JDiag[0+0*6] = -cj-w[0]*(1.0E2/8.99E2)+dwdx[6]; + JDiag[1+0*6] = -cj-w[6]*1.202935161794779E-2+dwdx[19]; + JDiag[2+0*6] = -cj-w[31]*8.196729508204918E-9+dwdx[32]; + JDiag[3+0*6] = -dwdx[36]-w[23]*w[24]*w[25]*dwdx[34]-1.0; + JDiag[4+0*6] = -dwdx[43]-1.0; + JDiag[5+0*6] = -dwdx[52]-1.0; +} + diff --git a/models/model_calvetti/model_calvetti_JSparse.cpp b/models/model_calvetti/model_calvetti_JSparse.cpp new file mode 100644 index 0000000000..053841de0c --- /dev/null +++ b/models/model_calvetti/model_calvetti_JSparse.cpp @@ -0,0 +1,71 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +#include //SUNMatrixContent_Sparse definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void JSparse_model_calvetti(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) { + JSparse->indexvals[0] = 0; + JSparse->indexvals[1] = 3; + JSparse->indexvals[2] = 0; + JSparse->indexvals[3] = 1; + JSparse->indexvals[4] = 3; + JSparse->indexvals[5] = 4; + JSparse->indexvals[6] = 0; + JSparse->indexvals[7] = 1; + JSparse->indexvals[8] = 2; + JSparse->indexvals[9] = 3; + JSparse->indexvals[10] = 4; + JSparse->indexvals[11] = 5; + JSparse->indexvals[12] = 0; + JSparse->indexvals[13] = 3; + JSparse->indexvals[14] = 4; + JSparse->indexvals[15] = 0; + JSparse->indexvals[16] = 1; + JSparse->indexvals[17] = 3; + JSparse->indexvals[18] = 4; + JSparse->indexvals[19] = 5; + JSparse->indexvals[20] = 0; + JSparse->indexvals[21] = 1; + JSparse->indexvals[22] = 2; + JSparse->indexvals[23] = 3; + JSparse->indexvals[24] = 4; + JSparse->indexvals[25] = 5; + JSparse->indexptrs[0] = 0; + JSparse->indexptrs[1] = 2; + JSparse->indexptrs[2] = 6; + JSparse->indexptrs[3] = 12; + JSparse->indexptrs[4] = 15; + JSparse->indexptrs[5] = 20; + JSparse->indexptrs[6] = 26; + JSparse->data[0] = -cj-w[0]*(1.0E2/8.99E2)+dwdx[6]; + JSparse->data[1] = w[0]*(1.0E2/8.99E2)-dwdx[6]+dwdx[7]*(w[23]*w[24]*2.0-w[20]*w[23]*w[24])-w[23]*w[24]*w[25]*dwdx[4]; + JSparse->data[2] = dwdx[16]; + JSparse->data[3] = -cj-w[6]*1.202935161794779E-2+dwdx[19]; + JSparse->data[4] = -dwdx[16]-w[23]*w[24]*w[25]*dwdx[14]; + JSparse->data[5] = w[6]*1.202935161794779E-2-dwdx[19]; + JSparse->data[6] = dwdx[28]; + JSparse->data[7] = dwdx[31]; + JSparse->data[8] = -cj-w[31]*8.196729508204918E-9+dwdx[32]; + JSparse->data[9] = -dwdx[28]-w[23]*w[24]*w[25]*dwdx[26]; + JSparse->data[10] = -dwdx[31]; + JSparse->data[11] = w[31]*8.196729508204918E-9-dwdx[32]; + JSparse->data[12] = dwdx[36]; + JSparse->data[13] = -dwdx[36]-w[23]*w[24]*w[25]*dwdx[34]-1.0; + JSparse->data[14] = 1.0; + JSparse->data[15] = dwdx[40]; + JSparse->data[16] = dwdx[43]; + JSparse->data[17] = -dwdx[40]-w[23]*w[24]*w[25]*dwdx[38]; + JSparse->data[18] = -dwdx[43]-1.0; + JSparse->data[19] = 1.0; + JSparse->data[20] = dwdx[47]; + JSparse->data[21] = dwdx[50]; + JSparse->data[22] = dwdx[52]; + JSparse->data[23] = -dwdx[47]-w[23]*w[24]*w[25]*dwdx[45]; + JSparse->data[24] = -dwdx[50]; + JSparse->data[25] = -dwdx[52]-1.0; +} + diff --git a/models/model_calvetti/model_calvetti_JSparseB.cpp b/models/model_calvetti/model_calvetti_JSparseB.cpp new file mode 100644 index 0000000000..15769a394a --- /dev/null +++ b/models/model_calvetti/model_calvetti_JSparseB.cpp @@ -0,0 +1,71 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +#include //SUNMatrixContent_Sparse definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void JSparseB_model_calvetti(SUNMatrixContent_Sparse JSparseB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype cj, const realtype *xB, const realtype *dx, const realtype *dxB, const realtype *w, const realtype *dwdx) { + JSparseB->indexvals[0] = 0; + JSparseB->indexvals[1] = 1; + JSparseB->indexvals[2] = 2; + JSparseB->indexvals[3] = 3; + JSparseB->indexvals[4] = 4; + JSparseB->indexvals[5] = 5; + JSparseB->indexvals[6] = 1; + JSparseB->indexvals[7] = 2; + JSparseB->indexvals[8] = 4; + JSparseB->indexvals[9] = 5; + JSparseB->indexvals[10] = 2; + JSparseB->indexvals[11] = 5; + JSparseB->indexvals[12] = 0; + JSparseB->indexvals[13] = 1; + JSparseB->indexvals[14] = 2; + JSparseB->indexvals[15] = 3; + JSparseB->indexvals[16] = 4; + JSparseB->indexvals[17] = 5; + JSparseB->indexvals[18] = 1; + JSparseB->indexvals[19] = 2; + JSparseB->indexvals[20] = 3; + JSparseB->indexvals[21] = 4; + JSparseB->indexvals[22] = 5; + JSparseB->indexvals[23] = 2; + JSparseB->indexvals[24] = 4; + JSparseB->indexvals[25] = 5; + JSparseB->indexptrs[0] = 0; + JSparseB->indexptrs[1] = 6; + JSparseB->indexptrs[2] = 10; + JSparseB->indexptrs[3] = 12; + JSparseB->indexptrs[4] = 18; + JSparseB->indexptrs[5] = 23; + JSparseB->indexptrs[6] = 26; + JSparseB->data[0] = cj+w[0]*(1.0E2/8.99E2)-dwdx[6]; + JSparseB->data[1] = -dwdx[16]; + JSparseB->data[2] = -dwdx[28]; + JSparseB->data[3] = -dwdx[36]; + JSparseB->data[4] = -dwdx[40]; + JSparseB->data[5] = -dwdx[47]; + JSparseB->data[6] = cj+w[6]*1.202935161794779E-2-dwdx[19]; + JSparseB->data[7] = -dwdx[31]; + JSparseB->data[8] = -dwdx[43]; + JSparseB->data[9] = -dwdx[50]; + JSparseB->data[10] = cj+w[31]*8.196729508204918E-9-dwdx[32]; + JSparseB->data[11] = -dwdx[52]; + JSparseB->data[12] = w[0]*(-1.0E2/8.99E2)+dwdx[6]-dwdx[7]*(w[23]*w[24]*2.0-w[20]*w[23]*w[24])+w[23]*w[24]*w[25]*dwdx[4]; + JSparseB->data[13] = dwdx[16]+w[23]*w[24]*w[25]*dwdx[14]; + JSparseB->data[14] = dwdx[28]+w[23]*w[24]*w[25]*dwdx[26]; + JSparseB->data[15] = dwdx[36]+w[23]*w[24]*w[25]*dwdx[34]+1.0; + JSparseB->data[16] = dwdx[40]+w[23]*w[24]*w[25]*dwdx[38]; + JSparseB->data[17] = dwdx[47]+w[23]*w[24]*w[25]*dwdx[45]; + JSparseB->data[18] = w[6]*(-1.202935161794779E-2)+dwdx[19]; + JSparseB->data[19] = dwdx[31]; + JSparseB->data[20] = -1.0; + JSparseB->data[21] = dwdx[43]+1.0; + JSparseB->data[22] = dwdx[50]; + JSparseB->data[23] = w[31]*(-8.196729508204918E-9)+dwdx[32]; + JSparseB->data[24] = -1.0; + JSparseB->data[25] = dwdx[52]+1.0; +} + diff --git a/models/model_calvetti/model_calvetti_Jy.cpp b/models/model_calvetti/model_calvetti_Jy.cpp new file mode 100644 index 0000000000..7a758ba5bf --- /dev/null +++ b/models/model_calvetti/model_calvetti_Jy.cpp @@ -0,0 +1,31 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void Jy_model_calvetti(double *nllh, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) { +switch(iy){ + case 0: + nllh[0] = amici::log((sigmay[0]*sigmay[0])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[0]*sigmay[0])*pow(my[0]-y[0],2.0)*5.0E-1; + break; + case 1: + nllh[0] = amici::log((sigmay[1]*sigmay[1])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[1]*sigmay[1])*pow(my[1]-y[1],2.0)*5.0E-1; + break; + case 2: + nllh[0] = amici::log((sigmay[2]*sigmay[2])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[2]*sigmay[2])*pow(my[2]-y[2],2.0)*5.0E-1; + break; + case 3: + nllh[0] = amici::log((sigmay[3]*sigmay[3])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[3]*sigmay[3])*pow(my[3]-y[3],2.0)*5.0E-1; + break; + case 4: + nllh[0] = amici::log((sigmay[4]*sigmay[4])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[4]*sigmay[4])*pow(my[4]-y[4],2.0)*5.0E-1; + break; + case 5: + nllh[0] = amici::log((sigmay[5]*sigmay[5])*3.141592653589793*2.0)*5.0E-1+1.0/(sigmay[5]*sigmay[5])*pow(my[5]-y[5],2.0)*5.0E-1; + break; +} +} + diff --git a/models/model_calvetti/model_calvetti_M.cpp b/models/model_calvetti/model_calvetti_M.cpp new file mode 100644 index 0000000000..314fa486b9 --- /dev/null +++ b/models/model_calvetti/model_calvetti_M.cpp @@ -0,0 +1,14 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void M_model_calvetti(realtype *M, const realtype t, const realtype *x, const realtype *p, const realtype *k) { + M[0+0*6] = 1.0; + M[1+1*6] = 1.0; + M[2+2*6] = 1.0; +} + diff --git a/models/model_calvetti/model_calvetti_dJydsigma.cpp b/models/model_calvetti/model_calvetti_dJydsigma.cpp new file mode 100644 index 0000000000..f69094f7ed --- /dev/null +++ b/models/model_calvetti/model_calvetti_dJydsigma.cpp @@ -0,0 +1,31 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void dJydsigma_model_calvetti(double *dJydsigma, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) { +switch(iy){ + case 0: + dJydsigma[0+0*1] = 1.0/(sigmay[0]*sigmay[0]*sigmay[0])*pow(my[0]-y[0],2.0)*-1.0+1.0/sigmay[0]; + break; + case 1: + dJydsigma[0+1*1] = 1.0/(sigmay[1]*sigmay[1]*sigmay[1])*pow(my[1]-y[1],2.0)*-1.0+1.0/sigmay[1]; + break; + case 2: + dJydsigma[0+2*1] = 1.0/(sigmay[2]*sigmay[2]*sigmay[2])*pow(my[2]-y[2],2.0)*-1.0+1.0/sigmay[2]; + break; + case 3: + dJydsigma[0+3*1] = 1.0/(sigmay[3]*sigmay[3]*sigmay[3])*pow(my[3]-y[3],2.0)*-1.0+1.0/sigmay[3]; + break; + case 4: + dJydsigma[0+4*1] = 1.0/(sigmay[4]*sigmay[4]*sigmay[4])*pow(my[4]-y[4],2.0)*-1.0+1.0/sigmay[4]; + break; + case 5: + dJydsigma[0+5*1] = 1.0/(sigmay[5]*sigmay[5]*sigmay[5])*pow(my[5]-y[5],2.0)*-1.0+1.0/sigmay[5]; + break; +} +} + diff --git a/models/model_calvetti/model_calvetti_dJydy.cpp b/models/model_calvetti/model_calvetti_dJydy.cpp new file mode 100644 index 0000000000..d69b3d9883 --- /dev/null +++ b/models/model_calvetti/model_calvetti_dJydy.cpp @@ -0,0 +1,31 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void dJydy_model_calvetti(double *dJydy, const int iy, const realtype *p, const realtype *k, const double *y, const double *sigmay, const double *my) { +switch(iy){ + case 0: + dJydy[0+0*1] = 1.0/(sigmay[0]*sigmay[0])*(my[0]*2.0-y[0]*2.0)*-5.0E-1; + break; + case 1: + dJydy[0+1*1] = 1.0/(sigmay[1]*sigmay[1])*(my[1]*2.0-y[1]*2.0)*-5.0E-1; + break; + case 2: + dJydy[0+2*1] = 1.0/(sigmay[2]*sigmay[2])*(my[2]*2.0-y[2]*2.0)*-5.0E-1; + break; + case 3: + dJydy[0+3*1] = 1.0/(sigmay[3]*sigmay[3])*(my[3]*2.0-y[3]*2.0)*-5.0E-1; + break; + case 4: + dJydy[0+4*1] = 1.0/(sigmay[4]*sigmay[4])*(my[4]*2.0-y[4]*2.0)*-5.0E-1; + break; + case 5: + dJydy[0+5*1] = 1.0/(sigmay[5]*sigmay[5])*(my[5]*2.0-y[5]*2.0)*-5.0E-1; + break; +} +} + diff --git a/models/model_calvetti/model_calvetti_dwdx.cpp b/models/model_calvetti/model_calvetti_dwdx.cpp new file mode 100644 index 0000000000..db1e8c82f4 --- /dev/null +++ b/models/model_calvetti/model_calvetti_dwdx.cpp @@ -0,0 +1,64 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +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) { + dwdx[0] = 1.0/(x[0]*x[0]*x[0])*-2.0; + dwdx[1] = k[1]*w[15]*dwdx[0]; + dwdx[2] = dwdx[1]; + dwdx[3] = x[3]*dwdx[2]; + dwdx[4] = dwdx[3]; + dwdx[5] = -1.0/(w[20]*w[20])*dwdx[4]; + dwdx[6] = w[0]*w[14]*w[21]*(2.0/3.1E1)+w[0]*x[0]*w[14]*dwdx[5]*(2.0/3.1E1); + dwdx[7] = x[0]*2.0; + dwdx[8] = 1.0/(x[1]*x[1]*x[1])*-2.0; + dwdx[9] = k[3]*w[1]*dwdx[8]; + dwdx[10] = dwdx[9]; + dwdx[11] = x[4]*dwdx[10]; + dwdx[12] = dwdx[9]; + dwdx[13] = x[3]*dwdx[12]; + dwdx[14] = dwdx[11]+dwdx[13]; + dwdx[15] = -1.0/(w[20]*w[20])*dwdx[14]; + dwdx[16] = w[0]*x[0]*w[14]*dwdx[15]*(2.0/3.1E1); + dwdx[17] = dwdx[11]; + dwdx[18] = -1.0/(w[26]*w[26])*dwdx[17]; + dwdx[19] = w[6]*w[27]*w[29]*(2.0/1.63E2)+x[1]*w[6]*w[29]*dwdx[18]*(2.0/1.63E2); + dwdx[20] = 1.0/(x[2]*x[2]*x[2])*-2.0; + dwdx[21] = k[5]*w[4]*dwdx[20]; + dwdx[22] = dwdx[21]; + dwdx[23] = x[4]*dwdx[22]; + dwdx[24] = k[5]*w[4]*x[5]*dwdx[20]; + dwdx[25] = x[2]*2.0; + dwdx[26] = dwdx[23]+dwdx[24]; + dwdx[27] = -1.0/(w[20]*w[20])*dwdx[26]; + dwdx[28] = w[0]*x[0]*w[14]*dwdx[27]*(2.0/3.1E1); + dwdx[29] = dwdx[23]+dwdx[24]; + dwdx[30] = -1.0/(w[26]*w[26])*dwdx[29]; + dwdx[31] = x[1]*w[6]*w[29]*dwdx[30]*(2.0/1.63E2); + dwdx[32] = w[11]*w[32]*w[33]*w[34]*w[36]*(1.0/6.1E1)+x[2]*w[32]*w[33]*w[34]*w[36]*dwdx[25]*(1.0/6.1E1); + dwdx[33] = w[18]; + dwdx[34] = dwdx[33]; + dwdx[35] = -1.0/(w[20]*w[20])*dwdx[34]; + dwdx[36] = w[0]*x[0]*w[14]*dwdx[35]*(2.0/3.1E1); + dwdx[37] = w[8]; + dwdx[38] = dwdx[37]; + dwdx[39] = -1.0/(w[20]*w[20])*dwdx[38]; + dwdx[40] = w[0]*x[0]*w[14]*dwdx[39]*(2.0/3.1E1); + dwdx[41] = dwdx[37]; + dwdx[42] = -1.0/(w[26]*w[26])*dwdx[41]; + dwdx[43] = x[1]*w[6]*w[29]*dwdx[42]*(2.0/1.63E2); + dwdx[44] = k[5]*w[4]*w[5]; + dwdx[45] = dwdx[44]; + dwdx[46] = -1.0/(w[20]*w[20])*dwdx[45]; + dwdx[47] = w[0]*x[0]*w[14]*dwdx[46]*(2.0/3.1E1); + dwdx[48] = dwdx[44]; + dwdx[49] = -1.0/(w[26]*w[26])*dwdx[48]; + dwdx[50] = x[1]*w[6]*w[29]*dwdx[49]*(2.0/1.63E2); + dwdx[51] = -1.0/(x[5]*x[5]); + dwdx[52] = x[2]*w[11]*w[32]*w[33]*w[36]*dwdx[51]*(1.0/6.1E1); +} + diff --git a/models/model_calvetti/model_calvetti_dydx.cpp b/models/model_calvetti/model_calvetti_dydx.cpp new file mode 100644 index 0000000000..2918eff3dc --- /dev/null +++ b/models/model_calvetti/model_calvetti_dydx.cpp @@ -0,0 +1,22 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +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) { + dydx[0+0*6] = 1.0; + dydx[1+1*6] = 1.0; + dydx[2+2*6] = 1.0; + dydx[3+0*6] = (x[3]*2.0)/x[0]+(1.0/(k[0]*k[0])*x[0]*4.0)/k[1]-(1.0/(k[0]*k[0])*x[0]*(x[3]*((k[0]*k[0])*k[1]*1.0/(x[0]*x[0])+(k[2]*k[2])*k[3]*1.0/(x[1]*x[1]))+x[4]*((k[2]*k[2])*k[3]*1.0/(x[1]*x[1])+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2]))+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2])*x[5])*2.0)/k[1]; + dydx[3+1*6] = (1.0/(k[0]*k[0])*(x[0]*x[0])*((k[2]*k[2])*k[3]*1.0/(x[1]*x[1]*x[1])*x[3]*2.0+(k[2]*k[2])*k[3]*1.0/(x[1]*x[1]*x[1])*x[4]*2.0))/k[1]; + dydx[3+2*6] = (1.0/(k[0]*k[0])*(x[0]*x[0])*((k[4]*k[4])*k[5]*1.0/(x[2]*x[2]*x[2])*x[4]*2.0+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2]*x[2])*x[5]*2.0))/k[1]; + dydx[3+3*6] = -(1.0/(k[0]*k[0])*(x[0]*x[0])*((k[0]*k[0])*k[1]*1.0/(x[0]*x[0])+(k[2]*k[2])*k[3]*1.0/(x[1]*x[1])))/k[1]; + dydx[3+4*6] = -(1.0/(k[0]*k[0])*(x[0]*x[0])*((k[2]*k[2])*k[3]*1.0/(x[1]*x[1])+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2])))/k[1]; + dydx[3+5*6] = -(1.0/(k[0]*k[0])*(k[4]*k[4])*k[5]*(x[0]*x[0])*1.0/(x[2]*x[2]))/k[1]; + dydx[4+3*6] = 1.0; + dydx[5+4*6] = 1.0; +} + diff --git a/models/model_calvetti/model_calvetti_root.cpp b/models/model_calvetti/model_calvetti_root.cpp new file mode 100644 index 0000000000..d7db388c6f --- /dev/null +++ b/models/model_calvetti/model_calvetti_root.cpp @@ -0,0 +1,15 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void root_model_calvetti(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx) { + root[0] = -t+1.0E1; + root[1] = -t+1.2E1; + root[2] = t-1.0E1; + root[3] = t-1.2E1; +} + diff --git a/models/model_calvetti/model_calvetti_sigmay.cpp b/models/model_calvetti/model_calvetti_sigmay.cpp new file mode 100644 index 0000000000..dfd7cca3b9 --- /dev/null +++ b/models/model_calvetti/model_calvetti_sigmay.cpp @@ -0,0 +1,17 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void sigmay_model_calvetti(double *sigmay, const realtype t, const realtype *p, const realtype *k) { + sigmay[0] = 1.0; + sigmay[1] = 1.0; + sigmay[2] = 1.0; + sigmay[3] = 1.0; + sigmay[4] = 1.0; + sigmay[5] = 1.0; +} + diff --git a/models/model_calvetti/model_calvetti_w.cpp b/models/model_calvetti/model_calvetti_w.cpp new file mode 100644 index 0000000000..9ecb4c3d41 --- /dev/null +++ b/models/model_calvetti/model_calvetti_w.cpp @@ -0,0 +1,49 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void w_model_calvetti(realtype *w, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) { + w[0] = 1.0/k[0]; + w[1] = k[2]*k[2]; + w[2] = 1.0/(x[1]*x[1]); + w[3] = k[3]*w[1]*w[2]; + w[4] = k[4]*k[4]; + w[5] = 1.0/(x[2]*x[2]); + w[6] = 1.0/k[2]; + w[7] = k[5]*w[4]*w[5]; + w[8] = w[3]+w[7]; + w[9] = x[4]*w[8]; + w[10] = k[5]*w[4]*w[5]*x[5]; + w[11] = x[2]*x[2]; + w[12] = h[1]*(1.0/3.1E1); + w[13] = k[1]*(1.0/2.0); + w[14] = w[13]-1.0; + w[15] = k[0]*k[0]; + w[16] = 1.0/(x[0]*x[0]); + w[17] = k[1]*w[15]*w[16]; + w[18] = w[3]+w[17]; + w[19] = x[3]*w[18]; + w[20] = w[9]+w[10]+w[19]; + w[21] = 1.0/w[20]; + w[22] = w[0]*x[0]*w[14]*w[21]*(2.0/3.1E1); + w[23] = 1.0/(k[0]*k[0]); + w[24] = 1.0/k[1]; + w[25] = x[0]*x[0]; + w[26] = w[9]+w[10]; + w[27] = 1.0/w[26]; + w[28] = k[3]*(1.0/2.0); + w[29] = k[1]+w[28]-1.0; + w[30] = x[1]*w[6]*w[27]*w[29]*(2.0/1.63E2); + w[31] = 1.0/k[4]; + w[32] = 1.0/(k[4]*k[4]*k[4]); + w[33] = 1.0/k[5]; + w[34] = 1.0/x[5]; + w[35] = k[5]*(1.0/2.0); + w[36] = k[1]+k[3]+w[35]-1.0; + w[37] = x[2]*w[11]*w[32]*w[33]*w[34]*w[36]*(1.0/6.1E1); +} + diff --git a/models/model_calvetti/model_calvetti_x0.cpp b/models/model_calvetti/model_calvetti_x0.cpp new file mode 100644 index 0000000000..1d405146b5 --- /dev/null +++ b/models/model_calvetti/model_calvetti_x0.cpp @@ -0,0 +1,17 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void x0_model_calvetti(realtype *x0, const realtype t, const realtype *p, const realtype *k) { + x0[0] = k[0]; + x0[1] = k[2]; + x0[2] = k[4]; + x0[3] = 1.0; + x0[4] = 1.0; + x0[5] = 1.0; +} + diff --git a/models/model_calvetti/model_calvetti_xdot.cpp b/models/model_calvetti/model_calvetti_xdot.cpp new file mode 100644 index 0000000000..5ced5c757c --- /dev/null +++ b/models/model_calvetti/model_calvetti_xdot.cpp @@ -0,0 +1,17 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +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) { + xdot[0] = h[0]*(-1.0/3.1E1)+w[12]+w[22]-dx[0]-w[0]*x[0]*(1.0E2/8.99E2)+1.29E2/8.99E2; + xdot[1] = w[30]-dx[1]-x[1]*w[6]*1.202935161794779E-2+1.816432094310117E-2; + xdot[2] = w[37]-dx[2]-x[2]*w[31]*8.196729508204918E-9+8.196729508204918E-3; + xdot[3] = h[0]*(1.0/3.1E1)-x[3]-w[12]-w[22]+w[0]*x[0]*(1.0E2/8.99E2)+w[23]*w[24]*w[25]*2.0-w[20]*w[23]*w[24]*w[25]-1.29E2/8.99E2; + xdot[4] = x[3]-x[4]-w[30]+x[1]*w[6]*1.202935161794779E-2-1.816432094310117E-2; + xdot[5] = x[4]-x[5]-w[37]+x[2]*w[31]*8.196729508204918E-9-8.196729508204918E-3; +} + diff --git a/models/model_calvetti/model_calvetti_y.cpp b/models/model_calvetti/model_calvetti_y.cpp new file mode 100644 index 0000000000..b7dd0c135a --- /dev/null +++ b/models/model_calvetti/model_calvetti_y.cpp @@ -0,0 +1,17 @@ + +#include "amici/symbolic_functions.h" +#include "amici/defines.h" //realtype definition +typedef amici::realtype realtype; +#include + +using namespace amici; + +void y_model_calvetti(double *y, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w) { + y[0] = x[0]; + y[1] = x[1]; + y[2] = x[2]; + y[3] = (1.0/(k[0]*k[0])*(x[0]*x[0])*2.0)/k[1]-(1.0/(k[0]*k[0])*(x[0]*x[0])*(x[3]*((k[0]*k[0])*k[1]*1.0/(x[0]*x[0])+(k[2]*k[2])*k[3]*1.0/(x[1]*x[1]))+x[4]*((k[2]*k[2])*k[3]*1.0/(x[1]*x[1])+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2]))+(k[4]*k[4])*k[5]*1.0/(x[2]*x[2])*x[5]))/k[1]; + y[4] = x[3]; + y[5] = x[4]; +} + diff --git a/models/model_calvetti/swig/CMakeLists.txt b/models/model_calvetti/swig/CMakeLists.txt new file mode 100644 index 0000000000..616cf5e437 --- /dev/null +++ b/models/model_calvetti/swig/CMakeLists.txt @@ -0,0 +1,33 @@ +cmake_minimum_required(VERSION 3.8) # swig_add_library + +if(POLICY CMP0078) + cmake_policy(SET CMP0078 OLD) +endif(POLICY CMP0078) + +find_package(SWIG REQUIRED) +include(${SWIG_USE_FILE}) + +find_package(PythonLibs REQUIRED) +include_directories(${PYTHON_INCLUDE_DIRS}) + +set(SWIG_LIBRARY_NAME _${PROJECT_NAME}) +set(CMAKE_SWIG_FLAGS "") +set_source_files_properties(${PROJECT_NAME}.i PROPERTIES CPLUSPLUS ON) + +# swig does not use INTERFACE_INCLUDE_DIRS of linked libraries, so add manually +get_target_property(AMICI_INCLUDE_DIRS Upstream::amici INTERFACE_INCLUDE_DIRECTORIES) +include_directories(${AMICI_INCLUDE_DIRS} ..) + +swig_add_library(${PROJECT_NAME} TYPE MODULE LANGUAGE python SOURCES ${PROJECT_NAME}.i) + +swig_link_libraries(${PROJECT_NAME} + ${PYTHON_LIBRARIES} + model) + +# configure module setup script +set(SETUP_PY_IN ${Amici_DIR}/model_setup.template.py) +set(SETUP_PY_OUT ${CMAKE_CURRENT_BINARY_DIR}/setup.py) + +add_custom_target(install-python + DEPENDS ${SWIG_LIBRARY_NAME} + COMMAND python ${SETUP_PY_OUT} install) diff --git a/models/model_calvetti/swig/model_calvetti.i b/models/model_calvetti/swig/model_calvetti.i new file mode 100644 index 0000000000..29a288c129 --- /dev/null +++ b/models/model_calvetti/swig/model_calvetti.i @@ -0,0 +1,14 @@ +%module model_calvetti +%import amici.i +// Add necessary symbols to generated header + +%{ +#include "wrapfunctions.h" +#include "amici/model_ode.h" +#include "amici/model_dae.h" +using namespace amici; +%} + + +// Process symbols in header +%include "wrapfunctions.h" diff --git a/models/model_calvetti/wrapfunctions.cpp b/models/model_calvetti/wrapfunctions.cpp new file mode 100644 index 0000000000..d78ca43b7a --- /dev/null +++ b/models/model_calvetti/wrapfunctions.cpp @@ -0,0 +1,7 @@ +#include "amici/model.h" +#include "wrapfunctions.h" + +std::unique_ptr getModel() { + return std::unique_ptr(new Model_model_calvetti()); +} + diff --git a/models/model_calvetti/wrapfunctions.h b/models/model_calvetti/wrapfunctions.h new file mode 100644 index 0000000000..dbd0a7e789 --- /dev/null +++ b/models/model_calvetti/wrapfunctions.h @@ -0,0 +1,8 @@ +#ifndef _amici_wrapfunctions_h +#define _amici_wrapfunctions_h + +#include "model_calvetti.h" + +std::unique_ptr getModel(); + +#endif /* _amici_wrapfunctions_h */ diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h index 8063b12f22..c239e42acb 100644 --- a/models/model_dirac/model_dirac.h +++ b/models/model_dirac/model_dirac.h @@ -1,6 +1,6 @@ #ifndef _amici_model_dirac_h #define _amici_model_dirac_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -47,6 +47,7 @@ class Model_model_dirac : public amici::Model_ODE { 0, 0, 0, + {}, 3, 0, 1, @@ -60,7 +61,7 @@ class Model_model_dirac : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_dirac(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_dirac(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h index cbffb64161..053b359cbb 100644 --- a/models/model_events/model_events.h +++ b/models/model_events/model_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_events_h #define _amici_model_events_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -61,6 +61,7 @@ class Model_model_events : public amici::Model_ODE { 0, 0, 0, + {}, 4, 0, 1, @@ -74,7 +75,7 @@ class Model_model_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_events(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h index d9242e15df..ccc2c19954 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_h #define _amici_model_jakstat_adjoint_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -50,6 +50,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { 1, 5, 0, + {}, 18, 8, 1, @@ -63,7 +64,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index de22e5f844..4ea6082568 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_o2_h #define _amici_model_jakstat_adjoint_o2_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -50,6 +50,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { 2, 30, 0, + {}, 384, 8, 154, @@ -63,7 +64,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h index 193aa82106..76a87961e8 100644 --- a/models/model_nested_events/model_nested_events.h +++ b/models/model_nested_events/model_nested_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_nested_events_h #define _amici_model_nested_events_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -50,6 +50,7 @@ class Model_model_nested_events : public amici::Model_ODE { 0, 0, 0, + {}, 1, 0, 0, @@ -63,7 +64,7 @@ class Model_model_nested_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_nested_events(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_nested_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h index 3ef5d5f1d6..725647012c 100644 --- a/models/model_neuron/model_neuron.h +++ b/models/model_neuron/model_neuron.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_h #define _amici_model_neuron_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -64,6 +64,7 @@ class Model_model_neuron : public amici::Model_ODE { 0, 0, 0, + {}, 4, 1, 1, @@ -77,7 +78,7 @@ class Model_model_neuron : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h index 20ff8e8ab5..15428acbf7 100644 --- a/models/model_neuron_o2/model_neuron_o2.h +++ b/models/model_neuron_o2/model_neuron_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_o2_h #define _amici_model_neuron_o2_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -66,6 +66,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { 2, 0, 0, + {}, 27, 1, 8, @@ -79,7 +80,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron_o2(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h index db6cb508e2..a7203c6a8a 100644 --- a/models/model_robertson/model_robertson.h +++ b/models/model_robertson/model_robertson.h @@ -1,6 +1,6 @@ #ifndef _amici_model_robertson_h #define _amici_model_robertson_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -48,6 +48,7 @@ class Model_model_robertson : public amici::Model_DAE { 2, 1, 0, + {}, 9, 2, 2, @@ -61,7 +62,7 @@ class Model_model_robertson : public amici::Model_DAE { virtual amici::Model* clone() const override { return new Model_model_robertson(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { J_model_robertson(J, t, x, p, k, h, cj, dx, w, dwdx); diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h index 9334895607..70c1d08538 100644 --- a/models/model_steadystate/model_steadystate.h +++ b/models/model_steadystate/model_steadystate.h @@ -1,6 +1,6 @@ #ifndef _amici_model_steadystate_h #define _amici_model_steadystate_h -/* Generated by amiwrap (R2017b) b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc */ +/* Generated by amiwrap (R2017b) d43d87529b7be0983c74fe22a73d8cf2b64bc286 */ #include #include #include "amici/defines.h" @@ -47,6 +47,7 @@ class Model_model_steadystate : public amici::Model_ODE { 2, 1, 0, + {}, 9, 2, 2, @@ -60,7 +61,7 @@ class Model_model_steadystate : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_steadystate(*this); }; - const std::string getAmiciCommit() const override { return "b16b69612f19d7dd9ea78e39ea8b3c89bb4b6efc"; }; + const std::string getAmiciCommit() const override { return "d43d87529b7be0983c74fe22a73d8cf2b64bc286"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_steadystate(J, t, x, p, k, h, w, dwdx); diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000000..8917fb87d5 --- /dev/null +++ b/pytest.ini @@ -0,0 +1,4 @@ +[pytest] +# hundreds of SBML <=5.17 warnings +filterwarnings = + ignore:.*inspect.getargspec\(\) is deprecated.*:DeprecationWarning diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index 00ecfd58ce..de7e58b015 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -111,6 +111,8 @@ '(realtype *dJydy, const int iy, const realtype *p, ' 'const realtype *k, const realtype *y, ' 'const realtype *sigmay, const realtype *my)', + 'sparse': + True, }, 'dwdp': { 'signature': @@ -1290,31 +1292,34 @@ def _generateSparseSymbol(self, name): """ matrix = self.eq(name) - idx = 0 - sparseMatrix = sp.zeros(matrix.rows, matrix.cols) - symbolList = [] - sparseList = [] - symbolColPtrs = [] - symbolRowVals = [] - for col in range(0, matrix.cols): - symbolColPtrs.append(idx) - for row in range(0, matrix.rows): - if not (matrix[row, col] == 0): - symbolName = f'{name}{idx}' - sparseMatrix[row, col] = sp.Symbol(symbolName) - symbolList.append(symbolName) - sparseList.append(matrix[row, col]) - symbolRowVals.append(row) - idx += 1 - symbolColPtrs.append(idx) - sparseList = sp.Matrix(sparseList) + if name == 'dJydy': + # One entry per y-slice + self._colptrs[name] = [] + self._rowvals[name] = [] + self._sparseeqs[name] = [] + self._sparsesyms[name] = [] + self._syms[name] = [] + base_index = 0 + for iy in range(self.ny()): + symbolColPtrs, symbolRowVals, sparseList, symbolList, \ + sparseMatrix = csc_matrix(matrix[iy, :], name, + base_index=base_index) + base_index += len(symbolList) + self._colptrs[name].append(symbolColPtrs) + self._rowvals[name].append(symbolRowVals) + self._sparseeqs[name].append(sparseList) + self._sparsesyms[name].append(symbolList) + self._syms[name].append(sparseMatrix) + else: + symbolColPtrs, symbolRowVals, sparseList, symbolList, \ + sparseMatrix = csc_matrix(matrix, name) - self._colptrs[name] = symbolColPtrs - self._rowvals[name] = symbolRowVals - self._sparseeqs[name] = sparseList - self._sparsesyms[name] = symbolList - self._syms[name] = sparseMatrix + self._colptrs[name] = symbolColPtrs + self._rowvals[name] = symbolRowVals + self._sparseeqs[name] = sparseList + self._sparsesyms[name] = symbolList + self._syms[name] = sparseMatrix def _compute_equation(self, name): """computes the symbolic formula for a symbolic variable @@ -2162,9 +2167,6 @@ def _write_function_index(self, function, indextype): """ - # function signature - signature = f'(sunindextype *{indextype})' - if indextype == 'colptrs': values = self.model.colptrs(function) elif indextype == 'rowvals': @@ -2173,15 +2175,29 @@ def _write_function_index(self, function, indextype): raise ValueError('Invalid value for type, must be colptr or ' 'rowval') + # function signature + if function in multiobs_functions: + signature = f'(sunindextype *{indextype}, int index)' + else: + signature = f'(sunindextype *{indextype})' + lines = [ '#include "sundials/sundials_types.h"', '', f'void {function}_{indextype}_{self.modelName}{signature}{{', ] - lines.extend( - [' ' * 4 + f'{indextype}[{index}] = {value};' - for index, value in enumerate(values)] - ) + if function in multiobs_functions: + # list of index vectors + cases = {switch_case: [' ' * 4 + f'{indextype}[{index}] = {value};' + for index, value in enumerate(idx_vector)] + for switch_case, idx_vector in enumerate(values)} + lines.extend(getSwitchStatement('index', cases, 1)) + else: + # single index vector + lines.extend( + [' ' * 4 + f'{indextype}[{index}] = {value};' + for index, value in enumerate(values)] + ) lines.append('}') with open(os.path.join( self.modelPath, @@ -2207,7 +2223,8 @@ def _getFunctionBody(self, function, symbol): lines = [] - if min(symbol.shape) == 0: + if len(symbol) == 0 or (isinstance(symbol, sp.Matrix) + and min(symbol.shape) == 0): return lines if not self.allow_reinit_fixpar_initcond \ @@ -2235,8 +2252,12 @@ def _getFunctionBody(self, function, symbol): lines.extend(getSwitchStatement('ip', cases, 1)) elif function in multiobs_functions: - cases = {iobs : self._getSymLines(symbol[:, iobs], function, 0) - for iobs in range(self.model.ny())} + if function == 'dJydy': + cases = {iobs: self._getSymLines(symbol[iobs], function, 0) + for iobs in range(self.model.ny())} + else: + cases = {iobs : self._getSymLines(symbol[:, iobs], function, 0) + for iobs in range(self.model.ny())} lines.extend(getSwitchStatement('iy', cases, 1)) else: @@ -2306,6 +2327,9 @@ def _writeModelHeader(self): 'NDWDP': str(len(self.model.eq('dwdp'))), 'NDWDX': str(len(self.model.sparsesym('dwdx'))), 'NDXDOTDW': str(len(self.model.sparsesym('dxdotdw'))), + 'NDJYDY': 'std::vector{%s}' + % ','.join(str(len(x)) + for x in self.model.sparsesym('dJydy')), 'NNZ': str(len(self.model.sparsesym('JSparse'))), 'UBW': str(self.model.nx_solver()), 'LBW': str(self.model.nx_solver()), @@ -2339,21 +2363,21 @@ def _writeModelHeader(self): for fun in [ 'w', 'dwdp', 'dwdx', 'x_rdata', 'x_solver', 'total_cl', 'dxdotdw', - 'dxdotdp', 'JSparse', 'JSparseB', + 'dxdotdp', 'JSparse', 'JSparseB', 'dJydy' ]: tplData[f'{fun.upper()}_DEF'] = \ - get_function_definition(fun, self.modelName) + get_function_extern_declaration(fun, self.modelName) tplData[f'{fun.upper()}_IMPL'] = \ - get_function_implementation(fun, self.modelName) + get_model_override_implementation(fun, self.modelName) if fun in sparse_functions: tplData[f'{fun.upper()}_COLPTRS_DEF'] = \ - get_sunindex_definition(fun, self.modelName, 'colptrs') + get_sunindex_extern_declaration(fun, self.modelName, 'colptrs') tplData[f'{fun.upper()}_COLPTRS_IMPL'] = \ - get_sunindex_implementation(fun, self.modelName, 'colptrs') + get_sunindex_override_implementation(fun, self.modelName, 'colptrs') tplData[f'{fun.upper()}_ROWVALS_DEF'] = \ - get_sunindex_definition(fun, self.modelName, 'rowvals') + get_sunindex_extern_declaration(fun, self.modelName, 'rowvals') tplData[f'{fun.upper()}_ROWVALS_IMPL'] = \ - get_sunindex_implementation(fun, self.modelName, 'rowvals') + get_sunindex_override_implementation(fun, self.modelName, 'rowvals') if self.model.nx_solver() == self.model.nx_rdata(): tplData['X_RDATA_DEF'] = '' @@ -2639,8 +2663,8 @@ def strip_pysb(symbol): return symbol -def get_function_definition(fun, name): - """Constructs the function definition for a given function +def get_function_extern_declaration(fun, name): + """Constructs the extern function declaration for a given function Arguments: fun: function name @type str @@ -2656,8 +2680,8 @@ def get_function_definition(fun, name): f'extern void {fun}_{name}{functions[fun]["signature"]};' -def get_sunindex_definition(fun, name, indextype): - """Constructs the function definition for an index function of a given +def get_sunindex_extern_declaration(fun, name, indextype): + """Constructs the function declaration for an index function of a given function Arguments: @@ -2666,17 +2690,18 @@ def get_sunindex_definition(fun, name, indextype): indextype: index function {'colptrs', 'rowvals'} @type str Returns: - c++ function definition string + c++ function declaration string Raises: """ + index_arg = ', int index' if fun in multiobs_functions else '' return \ - f'extern void {fun}_{indextype}_{name}(sunindextype *{indextype});' + f'extern void {fun}_{indextype}_{name}(sunindextype *{indextype}{index_arg});' -def get_function_implementation(fun, name): - """Constructs the function implementation for a given function +def get_model_override_implementation(fun, name): + """Constructs amici::Model::* override implementation for a given function Arguments: fun: function name @type str @@ -2701,9 +2726,9 @@ def get_function_implementation(fun, name): ) -def get_sunindex_implementation(fun, name, indextype): - """Constructs the function implementation for an index function of a given - function +def get_sunindex_override_implementation(fun, name, indextype): + """Constructs the amici::Model:: function implementation for an index + function of a given function Arguments: fun: function name @type str @@ -2716,6 +2741,9 @@ def get_sunindex_implementation(fun, name, indextype): Raises: """ + index_arg = ', int index' if fun in multiobs_functions else '' + index_arg_eval = ', index' if fun in multiobs_functions else '' + return \ 'virtual void f{fun}_{indextype}{signature} override {{\n' \ '{ind8}{fun}_{indextype}_{name}{eval_signature};\n' \ @@ -2725,8 +2753,8 @@ def get_sunindex_implementation(fun, name, indextype): fun=fun, indextype=indextype, name=name, - signature=f'(sunindextype *{indextype})', - eval_signature=f'({indextype})', + signature=f'(sunindextype *{indextype}{index_arg})', + eval_signature=f'({indextype}{index_arg_eval})', ) @@ -2806,3 +2834,43 @@ def getSwitchStatement(condition, cases, lines.append(indentation_level * indentation_step + '}') return lines + + +def csc_matrix(matrix, name, base_index=0): + """Generates the sparse symbolic identifiers, symbolic identifiers, + sparse matrix, column pointers and row values for a symbolic + variable + + Arguments: + matrix: dense matrix to be sparsified @type sp.Matrix + name: name of the symbolic variable @type str + base_index: index for first symbol name, defaults to 0 + + Returns: + symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix + Raises: + + """ + idx = 0 + symbol_name_idx = base_index + sparseMatrix = sp.zeros(matrix.rows, matrix.cols) + symbolList = [] + sparseList = [] + symbolColPtrs = [] + symbolRowVals = [] + for col in range(0, matrix.cols): + symbolColPtrs.append(idx) + for row in range(0, matrix.rows): + if not (matrix[row, col] == 0): + symbolName = f'{name}{symbol_name_idx}' + sparseMatrix[row, col] = sp.Symbol(symbolName) + symbolList.append(symbolName) + sparseList.append(matrix[row, col]) + symbolRowVals.append(row) + idx += 1 + symbol_name_idx += 1 + + symbolColPtrs.append(idx) + sparseList = sp.Matrix(sparseList) + + return symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix diff --git a/python/amici/pandas.py b/python/amici/pandas.py index 38f77dad42..feef1a9481 100644 --- a/python/amici/pandas.py +++ b/python/amici/pandas.py @@ -417,7 +417,7 @@ def _get_specialized_fixed_parameters( cond = copy.deepcopy(condition) for field in overwrite: cond[field] = overwrite[field] - return [cond[name] for name in _get_names_or_ids( + return [float(cond[name]) for name in _get_names_or_ids( model, 'FixedParameter', by_id=by_id)] @@ -447,7 +447,7 @@ def constructEdataFromDataFrame(df, model, condition, by_id=False): # timepoints df = df.sort_values(by='time', ascending=True) - edata.setTimepoints(df['time'].values) + edata.setTimepoints(df['time'].values.astype(float)) # get fixed parameters from condition overwrite_preeq = {} @@ -488,16 +488,16 @@ def constructEdataFromDataFrame(df, model, condition, by_id=False): # fill in presimulation time if 't_presim' in condition.keys(): - edata.t_presim = condition['t_presim'] + edata.t_presim = float(condition['t_presim']) # fill in data and stds for obs_index, obs in enumerate( _get_names_or_ids(model, 'Observable', by_id=by_id)): if obs in df.keys(): - edata.setObservedData(df[obs].values, obs_index) + edata.setObservedData(df[obs].values.astype(float), obs_index) if obs + '_std' in df.keys(): edata.setObservedDataStdDev( - df[obs + '_std'].values, obs_index + df[obs + '_std'].values.astype(float), obs_index ) return edata diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index 6fcbdfe32a..6d63255733 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -6,6 +6,7 @@ import libsbml as sbml import re import math +import itertools as itt import warnings from sympy.logic.boolalg import BooleanTrue as spTrue from sympy.logic.boolalg import BooleanFalse as spFalse @@ -381,6 +382,7 @@ def getSpeciesInitial(index, conc): sbml.formulaToL3String(initial_assignment.getMath()) ) if symMath is not None: + symMath = _parse_special_functions(symMath) _check_unsupported_functions(symMath, 'InitialAssignment') speciesInitial[index] = symMath @@ -589,6 +591,7 @@ def getElementStoichiometry(element): symMath = sp.sympify(element.getStoichiometry()) else: return sp.sympify(1.0) + symMath = _parse_special_functions(symMath) _check_unsupported_functions(symMath, 'Stoichiometry') return symMath @@ -627,6 +630,7 @@ def isConstant(specie): except: raise SBMLException(f'Kinetic law "{math}" contains an ' 'unsupported expression!') + symMath = _parse_special_functions(symMath) _check_unsupported_functions(symMath, 'KineticLaw') for r in reactions: elements = list(r.getListOfReactants()) \ @@ -678,6 +682,7 @@ def processRules(self): variable = sp.sympify(rule.getVariable()) # avoid incorrect parsing of pow(x, -1) in symengine formula = sp.sympify(sbml.formulaToL3String(rule.getMath())) + formula = _parse_special_functions(formula) _check_unsupported_functions(formula, 'Rule') if variable in stoichvars: @@ -1116,7 +1121,7 @@ def _check_unsupported_functions(sym, expression_type, full_sym=None): unsupported_functions = [ sp.functions.factorial, sp.functions.ceiling, sp.functions.floor, - sp.functions.Piecewise, spTrue, spFalse, sp.function.UndefinedFunction + sp.function.UndefinedFunction ] unsupp_fun_type = next( @@ -1150,6 +1155,46 @@ def _check_unsupported_functions(sym, expression_type, full_sym=None): _check_unsupported_functions(fun, expression_type) +def _parse_special_functions(sym): + """Recursively checks the symbolic expression for functions which have be + to parsed in a special way, such as piecewise functions + + Arguments: + sym: symbolic expressions @type sympy.Basic + + Returns: + + Raises: + """ + args = tuple(_parse_special_functions(arg) for arg in sym._args) + + # Do we have piecewise expressions? + if sym.__class__.__name__ == 'piecewise': + # how many condition-expression pairs will we have? + return sp.Piecewise(*grouper(args, 2, True)) + elif isinstance(sym, (sp.Function, sp.Mul, sp.Add)): + sym._args = args + + return sym + + +def grouper(iterable, n, fillvalue=None): + """Collect data into fixed-length chunks or blocks + + E.g. grouper('ABCDEFG', 3, 'x') --> ABC DEF Gxx" + + Arguments: + iterable: any iterable + n: chunk length + fillvalue: padding for last chunk if length < n + + Returns: + itertools.zip_longest of requested chunks + """ + args = [iter(iterable)] * n + return itt.zip_longest(*args, fillvalue=fillvalue) + + def assignmentRules2observables(sbml_model, filter_function=lambda *_: True): """Turn assignment rules into observables. diff --git a/python/amici/setup.template.py b/python/amici/setup.template.py index 94c5e4d9ec..b31b3ea014 100644 --- a/python/amici/setup.template.py +++ b/python/amici/setup.template.py @@ -68,6 +68,7 @@ def getAmiciLibs(): sources=sources, include_dirs=[os.getcwd(), os.path.join(amici_path, 'include'), + os.path.join(amici_path, 'ThirdParty/gsl'), os.path.join(amici_path, 'ThirdParty/sundials/include'), os.path.join(amici_path, 'ThirdParty/SuiteSparse/include'), *h5pkgcfg['include_dirs'], diff --git a/python/bin/amici_import_petab.py b/python/bin/amici_import_petab.py index 5a7a80f185..26a5ba3dc8 100755 --- a/python/bin/amici_import_petab.py +++ b/python/bin/amici_import_petab.py @@ -91,15 +91,22 @@ def get_fixed_parameters(condition_file_name, sbml_model, into constant parameters *within* the given `sbml_model`. """ - fixed_parameters_df = petab.get_condition_df(condition_file_name) - print(f'Condition table: {fixed_parameters_df.shape}') + condition_df = petab.get_condition_df(condition_file_name) + print(f'Condition table: {condition_df.shape}') # column names are model parameter names that should be made constant - fixed_parameters = list(fixed_parameters_df.columns) + # except for any overridden parameters + # (Could potentially still be made constant, but leaving them might + # increase model reusability) + fixed_parameters = list(condition_df.columns) try: fixed_parameters.remove('conditionName') except ValueError: pass + # remove overridden parameters + fixed_parameters = [p for p in fixed_parameters + if condition_df[p].dtype != 'O'] + # must be unique assert(len(fixed_parameters) == len(set(fixed_parameters))) diff --git a/scripts/buildSundials.sh b/scripts/buildSundials.sh index b464f4ecea..3fa05920cf 100755 --- a/scripts/buildSundials.sh +++ b/scripts/buildSundials.sh @@ -10,6 +10,15 @@ AMICI_PATH=$(cd $SCRIPT_PATH/.. && pwd) SUITESPARSE_ROOT="${AMICI_PATH}/ThirdParty/SuiteSparse" SUNDIALS_BUILD_PATH="${AMICI_PATH}/ThirdParty/sundials/build/" +# enable SuperLUMT support if library exists +SuperLUMT="" +if [[ -f ${AMICI_PATH}/ThirdParty/SuperLU_MT_3.1/lib/libsuperlu_mt_PTHREAD.a ]] +then + SuperLUMT="-DSUPERLUMT_ENABLE=ON -DBLAS_ENABLE=ON \ + -DSUPERLUMT_INCLUDE_DIR=${AMICI_PATH}/ThirdParty/SuperLU_MT_3.1/SRC/ \ + -DSUPERLUMT_LIBRARY_DIR=${AMICI_PATH}/ThirdParty/SuperLU_MT_3.1/lib/" +fi + mkdir -p ${SUNDIALS_BUILD_PATH} cd ${SUNDIALS_BUILD_PATH} @@ -27,6 +36,7 @@ cmake -DCMAKE_INSTALL_PREFIX="${SUNDIALS_BUILD_PATH}" \ -DKLU_ENABLE=ON \ -DKLU_LIBRARY_DIR="${SUITESPARSE_ROOT}/lib" \ -DKLU_INCLUDE_DIR="${SUITESPARSE_ROOT}/include" \ +${SuperLUMT} \ .. make diff --git a/scripts/buildSuperLUMT.sh b/scripts/buildSuperLUMT.sh new file mode 100755 index 0000000000..0164a7af89 --- /dev/null +++ b/scripts/buildSuperLUMT.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash +# +# Build SuperLUMT +# +set -e + +SCRIPT_PATH=$(dirname $BASH_SOURCE) +AMICI_PATH=$(cd ${SCRIPT_PATH}/.. && pwd) + +cd ${AMICI_PATH}/ThirdParty + +if [[ ! -d SuperLU_MT_3.1 ]]; then + if [[ ! -f superlu_mt_3.1.tar.gz ]]; then + wget https://crd-legacy.lbl.gov/~xiaoye/SuperLU/superlu_mt_3.1.tar.gz + fi + tar -xzf superlu_mt_3.1.tar.gz SuperLU_MT_3.1/ +fi + +cd SuperLU_MT_3.1/ +# interferes with std::queue +test -f SRC/queue && mv SRC/queue SRC/queue.bak +cp MAKE_INC/make.pthread make.inc +# Add -fPIC +sed -ri 's/(CFLAGS\W*=\W*.*)(\#.*)/\1-fPIC \2/' make.inc +# Use 64bit integers +sed -ri 's/# (CFLAGS\W*+=.*-D_LONGINT.*)/\1/' make.inc +sed -ri 's/# (FFLAGS\W*+=.*-fdefault-integer-8.*)/\1/' make.inc + +make superlulib diff --git a/scripts/installAmiciSource.sh b/scripts/installAmiciSource.sh index d325901d61..8db65157a2 100755 --- a/scripts/installAmiciSource.sh +++ b/scripts/installAmiciSource.sh @@ -22,6 +22,8 @@ rm -rf ${AMICI_PATH}/python/sdist/build/ # test install from setup.py python3 -m venv ${AMICI_PATH}/build/venv --clear source ${AMICI_PATH}/build/venv/bin/activate -pip3 install --upgrade pip setuptools pkgconfig wheel scipy matplotlib pysb coverage +# install wheel separately to prevent build_wheel fail in next step +pip3 install --upgrade wheel +pip3 install --upgrade pip setuptools pkgconfig scipy matplotlib pysb coverage pip3 install --verbose -e ${AMICI_PATH}/python/sdist deactivate diff --git a/scripts/runNotebook.sh b/scripts/runNotebook.sh index 90d3952c60..437fb7217a 100755 --- a/scripts/runNotebook.sh +++ b/scripts/runNotebook.sh @@ -1,5 +1,5 @@ #!/bin/bash -# Run jupyter notebooks as given on command line, show output only on error. +# Run jupyter notebooks as given on command line, show output only on error. # If a directory is provided, run all contained notebooks non-recursively. set -x set -e @@ -26,7 +26,7 @@ if [ $# -eq 0 ]; then fi source ${AMICI_PATH}/build/venv/bin/activate -pip3 show ipython || (pip3 install --upgrade jupyter && python3 -m ipykernel install --user --name amici --display-name "Python (amici)") +pip3 show ipython || (pip3 install --upgrade jupyter jupyter_contrib_nbextensions && python3 -m ipykernel install --user --name amici --display-name "Python (amici)") for arg in "$@"; do if [ -d $arg ]; then diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index e19dd305f2..5196084de7 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -365,9 +365,6 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { if (!seflag) { solver->reInit(t, &x, &dx); - /* make time derivative consistent */ - solver->calcIC(t, &x, &dx); - if (solver->getSensitivityOrder() >= SensitivityOrder::first) { if (solver->getSensitivityMethod() == SensitivityMethod::forward) { solver->sensReInit(&sx, &sdx); @@ -541,7 +538,7 @@ void ForwardProblem::prepDataSensis(int it) { model->fdsigmaydp(it, rdata, edata); model->fdJydy(it, rdata, edata); model->fdJydsigma(it, rdata, edata); - model->fdJydx(&dJydx, it, edata); + model->fdJydx(dJydx, it, edata); model->fdJydp(it, rdata, edata); } diff --git a/src/hdf5.cpp b/src/hdf5.cpp index a706ac88c0..562816a7ea 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -590,8 +590,7 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, if(attributeExists(file, datasetPath, "ordering")) { solver.setStateOrdering( - static_cast( - getIntScalarAttribute(file, datasetPath, "ordering"))); + getIntScalarAttribute(file, datasetPath, "ordering")); } if(attributeExists(file, datasetPath, "interpType")) { diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp index dca30dd80b..f227d01064 100644 --- a/src/interface_matlab.cpp +++ b/src/interface_matlab.cpp @@ -58,6 +58,7 @@ char amici_blasCBlasTransToBlasTrans(BLASTranspose trans) { case BLASTranspose::conjTrans: return 'C'; } + throw std::invalid_argument("Invalid argument to amici_blasCBlasTransToBlasTrans"); } /*! @@ -365,7 +366,7 @@ void setSolverOptions(const mxArray *prhs[], int nrhs, Solver &solver) } if (mxGetProperty(prhs[RHS_OPTIONS], 0, "ordering")) { - solver.setStateOrdering(static_cast(dbl2int(mxGetScalar(mxGetProperty(prhs[RHS_OPTIONS], 0, "ordering"))))); + solver.setStateOrdering(dbl2int(mxGetScalar(mxGetProperty(prhs[RHS_OPTIONS], 0, "ordering")))); } if (mxGetProperty(prhs[RHS_OPTIONS], 0, "stldet")) { diff --git a/src/model.cpp b/src/model.cpp index 9c189578cf..64183ea04a 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -73,26 +73,32 @@ void Model::fsJy(const int it, const std::vector& dJydx, const AmiVect } void Model::fdJydp(const int it, ReturnData *rdata, const ExpData *edata) { - - // dJydy nJ x nytrue x ny + // dJydy nJ, nytrue x ny // dydp nplist * ny // dJydp nplist x nJ // dJydsigma - getmy(it,edata); + getmy(it, edata); std::fill(dJydp.begin(), dJydp.end(), 0.0); for (int iyt = 0; iyt < nytrue; ++iyt) { if (isNaN(my.at(iyt))) continue; - // dJydp = 1.0 * dJydp + 1.0 * dJydy * dydp - amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, - nJ, nplist(), ny, - 1.0, &dJydy.at(iyt*nJ*ny), nJ, - dydp.data(), ny, - 1.0, dJydp.data(), nJ); - + if(wasPythonGenerated()) { + // dJydp = 1.0 * dJydp + 1.0 * dJydy * dydp + for(int iplist = 0; iplist < nplist(); ++iplist) { + dJydy[iyt].multiply( + gsl::span(&dJydp.at(iplist * nJ), nJ), + gsl::span(&dydp.at(iplist * ny), ny)); + } + } else { + amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, + nJ, nplist(), ny, + 1.0, &dJydy_matlab.at(iyt*nJ*ny), nJ, + dydp.data(), ny, + 1.0, dJydp.data(), nJ); + } // dJydp = 1.0 * dJydp + 1.0 * dJydsigma * dsigmaydp amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, nplist(), ny, @@ -119,22 +125,37 @@ void Model::fdJydp(const int it, ReturnData *rdata, const ExpData *edata) { } } -void Model::fdJydx(std::vector *dJydx, const int it, const ExpData *edata) { +void Model::fdJydx(std::vector& dJydx, const int it, const ExpData *edata) { + + // dJydy: nJ, ny x nytrue + // dydx : ny x nx_solver + // dJydx: nJ x nx_solver x nt + getmy(it, edata); - // dJydy nJ x ny x nytrue - // dydx ny x nx_solver - // dJydx nJ x nx_solver x nt - getmy(it,edata); for (int iyt = 0; iyt < nytrue; ++iyt) { if (isNaN(my.at(iyt))) continue; - // dJydy A[nyt,nJ,ny] * dydx B[ny,nx_solver] = dJydx C[it,nJ,nx_solver] - // slice slice - // M K K N M N - // lda ldb ldc - amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, - nJ, nx_solver, ny, 1.0, &dJydy.at(iyt*ny*nJ), nJ, dydx.data(), ny, 1.0, - &dJydx->at(it*nx_solver*nJ), nJ); + // dJydy A[nyt,nJ,ny] * dydx B[ny,nx_solver] = dJydx C[it,nJ,nx_solver] + // slice slice + // M K K N M N + // lda ldb ldc + + if(wasPythonGenerated()) { + for(int ix = 0; ix < nx_solver; ++ix) { + dJydy[iyt].multiply( + gsl::span(&dJydx.at(it * nx_solver * nJ + + ix * nJ), nJ), + gsl::span(&dydx.at(ix * ny), ny)); + } + } else { + amici_dgemm(BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, + nJ, nx_solver, ny, 1.0, &dJydy_matlab.at(iyt*ny*nJ), nJ, dydx.data(), ny, 1.0, + &dJydx.at(it*nx_solver*nJ), nJ); + } + } + + if(alwaysCheckFinite) { + amici::checkFinite(dJydx, "dJydx"); } } @@ -837,6 +858,7 @@ Model::Model(const int nx_rdata, const int ndwdx, const int ndwdp, const int ndxdotdw, + std::vector ndJydy, const int nnz, const int ubw, const int lbw, @@ -855,6 +877,7 @@ Model::Model(const int nx_rdata, ndwdx(ndwdx), ndwdp(ndwdp), ndxdotdw(ndxdotdw), + ndJydy(std::move(ndJydy)), nnz(nnz), nJ(nJ), ubw(ubw), @@ -879,7 +902,6 @@ Model::Model(const int nx_rdata, M(nx_solver, nx_solver), my(nytrue, 0.0), mz(nztrue, 0.0), - dJydy(nJ*nytrue*ny, 0.0), dJydsigma(nJ*nytrue*ny, 0.0), dJzdz(nJ*nztrue*nz, 0.0), dJzdsigma(nJ*nztrue*nz, 0.0), @@ -908,6 +930,20 @@ Model::Model(const int nx_rdata, x_pos_tmp(nx_solver), pscale(std::vector(p.size(), ParameterScaling::none)) { + // Can't use derivedClass::wasPythonGenerated() in ctor. + // Guess we are using Python if ndJydy is not empty + if(!this->ndJydy.empty()) { + if(static_cast(nytrue) != this->ndJydy.size()) + throw std::runtime_error("Number of elements in ndJydy is not equal " + " nytrue."); + + for(int iytrue = 0; iytrue < nytrue; ++iytrue) + dJydy.push_back(SUNMatrixWrapper(nJ, ny, this->ndJydy[iytrue], + CSC_MAT)); + } else { + dJydy_matlab = std::vector(nJ*nytrue*ny, 0.0); + } + requireSensitivitiesForAllParameters(); } @@ -1193,9 +1229,10 @@ void Model::fsigmay(const int it, ReturnData *rdata, const ExpData *edata) { } } - for(int i = 0; i < nytrue; ++i) - checkSigmaPositivity(sigmay[i], "sigmay"); - + for(int i = 0; i < nytrue; ++i) { + if(edata && !std::isnan(edata->getObservedData()[it * nytrue + i])) + checkSigmaPositivity(sigmay[i], "sigmay"); + } std::copy_n(sigmay.data(), nytrue, &rdata->sigmay[it * rdata->ny]); } @@ -1323,24 +1360,49 @@ void Model::fJrz(const int nroots, ReturnData *rdata, const ExpData * /*edata*/) void Model::fdJydy(const int it, const ReturnData *rdata, const ExpData *edata) { // load measurements to my - getmy(it,edata); - std::fill(dJydy.begin(),dJydy.end(),0.0); + getmy(it, edata); + + if(wasPythonGenerated()) { + for(int iytrue = 0; iytrue < nytrue; iytrue++) { + dJydy[iytrue].zero(); + fdJydy_colptrs(dJydy[iytrue].indexptrs(), iytrue); + fdJydy_rowvals(dJydy[iytrue].indexvals(), iytrue); + + if(isNaN(my.at(iytrue))) { + continue; + } - for(int iytrue = 0; iytrue < nytrue; iytrue++){ - if(!isNaN(my.at(iytrue))){ // get dJydy slice (ny) for current timepoint and observable - fdJydy(&dJydy.at(iytrue*ny*nJ), + fdJydy(dJydy[iytrue].data(), iytrue, unscaledParameters.data(), fixedParameters.data(), gety(it,rdata), sigmay.data(), my.data()); - } - } - if(alwaysCheckFinite) { - amici::checkFinite(dJydy, "dJydy"); + if(alwaysCheckFinite) { + amici::checkFinite(ndJydy[iytrue], dJydy[iytrue].data(), "dJydy"); + } + } + } else { + std::fill(dJydy_matlab.begin(), dJydy_matlab.end(), 0.0); + for(int iytrue = 0; iytrue < nytrue; iytrue++) { + if(isNaN(my.at(iytrue))) { + continue; + } + fdJydy(&dJydy_matlab.at(iytrue*ny*nJ), + iytrue, + unscaledParameters.data(), + fixedParameters.data(), + gety(it,rdata), + sigmay.data(), + my.data()); + } + if(alwaysCheckFinite) { + // get dJydy slice (ny) for current timepoint and observable + amici::checkFinite(dJydy_matlab, "dJydy"); + } } } diff --git a/src/model_header.ODE_template.h b/src/model_header.ODE_template.h index af59b54739..bde8831f89 100644 --- a/src/model_header.ODE_template.h +++ b/src/model_header.ODE_template.h @@ -42,10 +42,9 @@ extern void dJydsigmay_TPL_MODELNAME(realtype *dJydsigmay, const int iy, const realtype *p, const realtype *k, const realtype *y, const realtype *sigmay, const realtype *my); -extern void dJydy_TPL_MODELNAME(realtype *dJydy, const int iy, - const realtype *p, const realtype *k, - const realtype *y, const realtype *sigmay, - const realtype *my); +TPL_DJYDY_DEF +TPL_DJYDY_COLPTRS_DEF +TPL_DJYDY_ROWVALS_DEF TPL_DWDP_DEF TPL_DWDX_DEF TPL_DWDX_COLPTRS_DEF @@ -116,6 +115,7 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { TPL_NDWDX, // ndwdx TPL_NDWDP, // ndwdp TPL_NDXDOTDW, // ndxdotdw + TPL_NDJYDY, // ndjydy TPL_NNZ, // nnz TPL_UBW, // ubw TPL_LBW, // lbw @@ -186,15 +186,15 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { } TPL_JSPARSE_IMPL - + TPL_JSPARSE_COLPTRS_IMPL - + TPL_JSPARSE_ROWVALS_IMPL TPL_JSPARSEB_IMPL - + TPL_JSPARSEB_COLPTRS_IMPL - + TPL_JSPARSEB_ROWVALS_IMPL /** model specific implementation of fJrz @@ -280,21 +280,6 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { dJydsigmay_TPL_MODELNAME(dJydsigma, iy, p, k, y, sigmay, my); } - /** model specific implementation of fdJydy - * @param dJydy partial derivative of time-resolved measurement negative - *log-likelihood Jy - * @param iy output index - * @param p parameter vector - * @param k constant vector - * @param y model output at timepoint - * @param sigmay measurement standard deviation at timepoint - * @param my measurement at timepoint - **/ - virtual void fdJydy(realtype *dJydy, const int iy, const realtype *p, - const realtype *k, const realtype *y, - const realtype *sigmay, const realtype *my) override { - dJydy_TPL_MODELNAME(dJydy, iy, p, k, y, sigmay, my); - } /** model specific implementation of fdJzdsigma * @param dJzdsigma Sensitivity of event measurement @@ -454,16 +439,20 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { const realtype *p, const realtype *k, const int ip) override {} + TPL_DJYDY_IMPL + TPL_DJYDY_COLPTRS_IMPL + TPL_DJYDY_ROWVALS_IMPL + TPL_DWDP_IMPL - + TPL_DWDX_IMPL - + TPL_DXDOTDW_IMPL - + TPL_DXDOTDW_COLPTRS_IMPL - + TPL_DXDOTDW_ROWVALS_IMPL - + TPL_DXDOTDP_IMPL /** model specific implementation of fdydx @@ -716,9 +705,9 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { const realtype *h) override {} TPL_X_RDATA_IMPL - + TPL_X_SOLVER_IMPL - + TPL_TOTAL_CL_IMPL /** @@ -812,7 +801,7 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { virtual const std::string getAmiciCommit() const override { return "TPL_AMICI_COMMIT_STRING"; } - + virtual bool wasPythonGenerated() const override { return true; } diff --git a/src/model_ode.cpp b/src/model_ode.cpp index b70bb6235f..9f9f93c72e 100644 --- a/src/model_ode.cpp +++ b/src/model_ode.cpp @@ -157,7 +157,8 @@ namespace amici { unscaledParameters.data(), fixedParameters.data(), h.data(), plist_[ip], w.data()); if (nw > 0) - dxdotdw.multiply(dxdotdp.data(ip), &dwdp.at(nw * ip)); + dxdotdw.multiply(gsl::span(dxdotdp.data(ip), nx_solver), + gsl::span(&dwdp.at(nw * ip), nw)); } } else { // matlab generated diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index f71a2e1959..fa62a0252f 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -1,9 +1,7 @@ #include "amici/newton_solver.h" -#include "amici/defines.h" #include "amici/model.h" #include "amici/solver.h" -#include "amici/vector.h" #include "amici/steadystateproblem.h" #include "amici/forwardproblem.h" #include "amici/rdata.h" @@ -21,14 +19,6 @@ namespace amici { NewtonSolver::NewtonSolver(realtype *t, AmiVector *x, Model *model, ReturnData *rdata) : model(model), rdata(rdata), xdot(x->getLength()), dx(x->getLength()) { - /** - * default constructor, initializes all members with the provided objects - * - * @param t pointer to time variable - * @param x pointer to state variables - * @param model pointer to the AMICI model object - * @param rdata pointer to the return data object - */ this->t = t; this->x = x; } @@ -38,23 +28,6 @@ NewtonSolver::NewtonSolver(realtype *t, AmiVector *x, Model *model, ReturnData * std::unique_ptr NewtonSolver::getSolver( realtype *t, AmiVector *x, LinearSolver linsolType, Model *model, ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol) { - /** - * Tries to determine the steady state of the ODE system by a Newton - * solver, uses forward intergration, if the Newton solver fails, - * restarts Newton solver, if integration fails. - * Computes steady state sensitivities - * - * @param t pointer to time variable - * @param x pointer to state variables - * @param linsolType integer indicating which linear solver to use - * @param model pointer to the AMICI model object - * @param rdata pointer to the return data object - * @param maxlinsteps maximum number of allowed linear steps per Newton step for steady state computation - * @param maxsteps maximum number of allowed Newton steps for steady state computation - * @param atol absolute tolerance - * @param rtol relative tolerance - * @return solver NewtonSolver according to the specified linsolType - */ std::unique_ptr solver; @@ -89,6 +62,8 @@ std::unique_ptr NewtonSolver::getSolver( throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); /* SPARSE SOLVERS */ + case LinearSolver::SuperLUMT: + throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver"); case LinearSolver::KLU: solver.reset(new NewtonSolverSparse(t, x, model, rdata)); break; @@ -107,15 +82,6 @@ std::unique_ptr NewtonSolver::getSolver( /* ------------------------------------------------------------------------- */ void NewtonSolver::getStep(int ntry, int nnewt, AmiVector *delta) { - /** - * Computes the solution of one Newton iteration - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - * @param delta containing the RHS of the linear system, will be - * overwritten by solution to the linear system - */ this->prepareLinearSystem(ntry, nnewt); @@ -126,11 +92,6 @@ void NewtonSolver::getStep(int ntry, int nnewt, AmiVector *delta) { /* ------------------------------------------------------------------------- */ void NewtonSolver::computeNewtonSensis(AmiVectorArray *sx) { - /** - * Computes steady state sensitivities - * - * @param sx pointer to state variable sensitivities - */ prepareLinearSystem(0, -1); model->fdxdotdp(*t, x, &dx); @@ -146,22 +107,11 @@ void NewtonSolver::computeNewtonSensis(AmiVectorArray *sx) { /* - Dense linear solver --------------------------------------------------- */ /* ------------------------------------------------------------------------- */ -/* Derived class for dense linear solver */ NewtonSolverDense::NewtonSolverDense(realtype *t, AmiVector *x, Model *model, ReturnData *rdata) : NewtonSolver(t, x, model, rdata), Jtmp(model->nx_solver, model->nx_solver), linsol(SUNLinSol_Dense(x->getNVector(), Jtmp.get())) { - /** - * default constructor, initializes all members with the provided objects - * and - * initializes temporary storage objects - * - * @param t pointer to time variable - * @param x pointer to state variables - * @param model pointer to the AMICI model object - * @param rdata pointer to the return data object - */ int status = SUNLinSolInitialize_Dense(linsol); if(status != AMICI_SUCCESS) throw NewtonFailure(status, "SUNLinSolInitialize_Dense"); @@ -170,15 +120,6 @@ NewtonSolverDense::NewtonSolverDense(realtype *t, AmiVector *x, Model *model, Re /* ------------------------------------------------------------------------- */ void NewtonSolverDense::prepareLinearSystem(int /*ntry*/, int /*nnewt*/) { - /** - * Writes the Jacobian for the Newton iteration and passes it to the linear - * solver - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - */ - model->fJ(*t, 0.0, x, &dx, &xdot, Jtmp.get()); int status = SUNLinSolSetup_Dense(linsol, Jtmp.get()); if(status != AMICI_SUCCESS) @@ -188,18 +129,11 @@ void NewtonSolverDense::prepareLinearSystem(int /*ntry*/, int /*nnewt*/) { /* ------------------------------------------------------------------------- */ void NewtonSolverDense::solveLinearSystem(AmiVector *rhs) { - /** - * Solves the linear system for the Newton step - * - * @param rhs containing the RHS of the linear system, will be - * overwritten by solution to the linear system - */ - int status = SUNLinSolSolve_Dense(linsol, Jtmp.get(), rhs->getNVector(), rhs->getNVector(), 0.0); // last argument is tolerance and does not have any influence on result - + if(status != AMICI_SUCCESS) throw NewtonFailure(status, "SUNLinSolSolve_Dense"); } @@ -221,15 +155,6 @@ NewtonSolverSparse::NewtonSolverSparse(realtype *t, AmiVector *x, Model *model, Jtmp(model->nx_solver, model->nx_solver, model->nnz, CSC_MAT), linsol(SUNKLU(x->getNVector(), Jtmp.get())) { - /** - * default constructor, initializes all members with the provided objects, - * initializes temporary storage objects and the klu solver - * - * @param t pointer to time variable - * @param x pointer to state variables - * @param model pointer to the AMICI model object - * @param rdata pointer to the return data object - */ int status = SUNLinSolInitialize_KLU(linsol); if(status != AMICI_SUCCESS) throw NewtonFailure(status, "SUNLinSolInitialize_KLU"); @@ -238,38 +163,22 @@ NewtonSolverSparse::NewtonSolverSparse(realtype *t, AmiVector *x, Model *model, /* ------------------------------------------------------------------------- */ void NewtonSolverSparse::prepareLinearSystem(int /*ntry*/, int /*nnewt*/) { - /** - * Writes the Jacobian for the Newton iteration and passes it to the linear - * solver - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - */ - /* Get sparse Jacobian */ model->fJSparse(*t, 0.0, x, &dx, &xdot, Jtmp.get()); int status = SUNLinSolSetup_KLU(linsol, Jtmp.get()); if(status != AMICI_SUCCESS) throw NewtonFailure(status, "SUNLinSolSetup_KLU"); -} // namespace amici +} /* ------------------------------------------------------------------------- */ void NewtonSolverSparse::solveLinearSystem(AmiVector *rhs) { - /** - * Solves the linear system for the Newton step - * - * @param rhs containing the RHS of the linear system,will be - * overwritten by solution to the linear system - */ - /* Pass pointer to the linear solver */ int status = SUNLinSolSolve_KLU(linsol, Jtmp.get(), rhs->getNVector(), rhs->getNVector(), 0.0); // last argument is tolerance and does not have any influence on result - + if(status != AMICI_SUCCESS) throw NewtonFailure(status, "SUNLinSolSolve_Dense"); } @@ -290,28 +199,11 @@ NewtonSolverIterative::NewtonSolverIterative(realtype *t, AmiVector *x, Model *m ns_t(model->nx_solver), ns_s(model->nx_solver), ns_r(model->nx_solver), ns_rt(model->nx_solver), ns_v(model->nx_solver), ns_Jv(model->nx_solver), ns_tmp(model->nx_solver), ns_Jdiag(model->nx_solver) { - /** - * default constructor, initializes all members with the provided objects - * @param t pointer to time variable - * @param x pointer to state variables - * @param model pointer to the AMICI model object - * @param rdata pointer to the return data object - */ } /* ------------------------------------------------------------------------- */ void NewtonSolverIterative::prepareLinearSystem(int ntry, int nnewt) { - /** - * Writes the Jacobian for the Newton iteration and passes it to the linear - * solver. - * Also wraps around getSensis for iterative linear solver. - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - */ - newton_try = ntry; i_newton = nnewt; if (nnewt == -1) { @@ -322,36 +214,12 @@ void NewtonSolverIterative::prepareLinearSystem(int ntry, int nnewt) { /* ------------------------------------------------------------------------- */ void NewtonSolverIterative::solveLinearSystem(AmiVector *rhs) { - /** - * Solves the linear system for the Newton step by passing it to - * linsolveSPBCG - * - * @param rhs containing the RHS of the linear system, will be - * overwritten by solution to the linear system - */ - linsolveSPBCG(newton_try, i_newton, rhs); rhs->minus(); } -void NewtonSolverIterative::linsolveSPBCG(int ntry,int nnewt, AmiVector *ns_delta) { - /** - * Iterative linear solver created from SPILS BiCG-Stab. - * Solves the linear system within each Newton step if iterative solver is - * chosen. - * - * @param ntry integer newton_try integer start number of Newton solver - * (1 or 2) - * @param nnewt integer number of current Newton step - * @param ns_delta Newton step - */ - - double rho; - double alpha; - double omega; - double res; - +void NewtonSolverIterative::linsolveSPBCG(int ntry, int nnewt, AmiVector *ns_delta) { xdot = *ns_delta; xdot.minus(); @@ -370,9 +238,9 @@ void NewtonSolverIterative::linsolveSPBCG(int ntry,int nnewt, AmiVector *ns_delt ns_v.reset(); ns_delta->reset(); ns_tmp.reset(); - rho = 1.0; - omega = 1.0; - alpha = 1.0; + double rho = 1.0; + double omega = 1.0; + double alpha = 1.0; // can be set to 0 at the moment model->fJv(*t, x, &dx, &xdot, ns_delta, &ns_Jv, 0.0); @@ -380,7 +248,7 @@ void NewtonSolverIterative::linsolveSPBCG(int ntry,int nnewt, AmiVector *ns_delt // ns_r = xdot - ns_Jv; N_VLinearSum(-1.0, ns_Jv.getNVector(), 1.0, xdot.getNVector(), ns_r.getNVector()); N_VDiv(ns_r.getNVector(), ns_Jdiag.getNVector(), ns_r.getNVector()); - res = sqrt(N_VDotProd(ns_r.getNVector(), ns_r.getNVector())); + double res = sqrt(N_VDotProd(ns_r.getNVector(), ns_r.getNVector())); ns_rt = ns_r; for (int i_linstep = 0; i_linstep < maxlinsteps; diff --git a/src/solver.cpp b/src/solver.cpp index a9fb67a2d6..72c927c95a 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -243,19 +243,19 @@ void Solver::initializeLinearSolver(const Model *model, AmiVector *x) { /* ITERATIVE SOLVERS */ case LinearSolver::SPGMR: - linearSolver = std::make_unique(*x, PREC_NONE, SUNSPGMR_MAXL_DEFAULT); + linearSolver = std::make_unique(*x); setLinearSolver(); setJacTimesVecFn(); break; case LinearSolver::SPBCG: - linearSolver = std::make_unique(*x, PREC_NONE, SUNSPBCGS_MAXL_DEFAULT); + linearSolver = std::make_unique(*x); setLinearSolver(); setJacTimesVecFn(); break; case LinearSolver::SPTFQMR: - linearSolver = std::make_unique(*x, PREC_NONE, SUNSPTFQMR_MAXL_DEFAULT); + linearSolver = std::make_unique(*x); setLinearSolver(); setJacTimesVecFn(); break; @@ -263,11 +263,24 @@ void Solver::initializeLinearSolver(const Model *model, AmiVector *x) { /* SPARSE SOLVERS */ case LinearSolver::KLU: - linearSolver = std::make_unique(*x, model->nnz, CSC_MAT, getStateOrdering()); + linearSolver = std::make_unique( + *x, model->nnz, CSC_MAT, + static_cast(getStateOrdering())); setLinearSolver(); setSparseJacFn(); break; +#ifdef SUNDIALS_SUPERLUMT + case LinearSolver::SuperLUMT: + // TODO state ordering + linearSolver = std::make_unique( + *x, model->nnz, CSC_MAT, + static_cast(getStateOrdering())); + + setLinearSolver(); + setSparseJacFn(); + break; +#endif default: throw AmiException("Invalid choice of solver: %d", static_cast(linsol)); } @@ -318,19 +331,19 @@ void Solver::initializeLinearSolverB(const Model *model, AmiVector *xB, const in /* ITERATIVE SOLVERS */ case LinearSolver::SPGMR: - linearSolverB = std::make_unique(*xB, PREC_NONE, SUNSPGMR_MAXL_DEFAULT); + linearSolverB = std::make_unique(*xB); setLinearSolverB(which); setJacTimesVecFnB(which); break; case LinearSolver::SPBCG: - linearSolverB = std::make_unique(*xB, PREC_NONE, SUNSPBCGS_MAXL_DEFAULT); + linearSolverB = std::make_unique(*xB); setLinearSolverB(which); setJacTimesVecFnB(which); break; case LinearSolver::SPTFQMR: - linearSolverB = std::make_unique(*xB, PREC_NONE, SUNSPTFQMR_MAXL_DEFAULT); + linearSolverB = std::make_unique(*xB); setLinearSolverB(which); setJacTimesVecFnB(which); break; @@ -338,11 +351,21 @@ void Solver::initializeLinearSolverB(const Model *model, AmiVector *xB, const in /* SPARSE SOLVERS */ case LinearSolver::KLU: - linearSolverB = std::make_unique(*xB, model->nnz, CSC_MAT, getStateOrdering()); + linearSolverB = std::make_unique( + *xB, model->nnz, CSC_MAT, + static_cast(getStateOrdering())); setLinearSolverB(which); setSparseJacFnB(which); break; - +#ifdef SUNDIALS_SUPERLUMT + case LinearSolver::SuperLUMT: + linearSolverB = std::make_unique( + *xB, model->nnz, CSC_MAT, + static_cast(getStateOrdering())); + setLinearSolverB(which); + setSparseJacFnB(which); + break; +#endif default: throw AmiException("Invalid choice of solver: %d", static_cast(linsol)); } @@ -742,18 +765,26 @@ void Solver::setInterpolationType(InterpolationType interpType) { this->interpType = interpType; } -StateOrdering Solver::getStateOrdering() const { +int Solver::getStateOrdering() const { return ordering; } -void Solver::setStateOrdering(StateOrdering ordering) { +void Solver::setStateOrdering(int ordering) { this->ordering = ordering; if (solverMemory && linsol == LinearSolver::KLU) { auto klu = dynamic_cast(linearSolver.get()); - klu->setOrdering(ordering); + klu->setOrdering(static_cast(ordering)); klu = dynamic_cast(linearSolverB.get()); - klu->setOrdering(ordering); + klu->setOrdering(static_cast(ordering)); + } +#ifdef SUNDIALS_SUPERLUMT + if (solverMemory && linsol == LinearSolver::SuperLUMT) { + auto klu = dynamic_cast(linearSolver.get()); + klu->setOrdering(static_cast(ordering)); + klu = dynamic_cast(linearSolverB.get()); + klu->setOrdering(static_cast(ordering)); } +#endif } int Solver::getStabilityLimitFlag() const { diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 83dcca112e..5ba9b4b0d1 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -180,10 +180,6 @@ void IDASolver::resetState(void *ami_mem, N_Vector yy0, N_Vector yp0) { /* current order */ ida_mem->ida_kk = 0; - - /* Initial setup not done yet */ - - ida_mem->ida_SetupDone = SUNFALSE; } void IDASolver::reInitPostProcessF(realtype *t, AmiVector *yout, @@ -199,7 +195,48 @@ void IDASolver::reInitPostProcessB(int which, realtype *t, AmiVector *yBout, void IDASolver::reInitPostProcess(void *ami_mem, realtype *t, AmiVector *yout, AmiVector *ypout, - realtype tout) {} + realtype tout) { + auto ida_mem = static_cast(ami_mem); + auto nst_tmp = ida_mem->ida_nst; + ida_mem->ida_nst = 0; + + auto status = IDASetStopTime(ida_mem, tout); + if(status != IDA_SUCCESS) + throw IDAException(status, "CVodeSetStopTime"); + + status = IDASolve(ami_mem, tout, t, yout->getNVector(), ypout->getNVector(), + IDA_ONE_STEP); + + if(status != IDA_SUCCESS) + throw IDAException(status, "reInitPostProcess"); + + ida_mem->ida_nst = nst_tmp+1; + if (ida_mem->ida_adjMallocDone == SUNTRUE) { + /* add new step to history array, this is copied from CVodeF */ + auto ia_mem = ida_mem->ida_adj_mem; + auto dt_mem = ia_mem->dt_mem; + + if (ida_mem->ida_nst % ia_mem->ia_nsteps == 0) { + /* currently not implemented, we should never get here as we + limit cv_mem->cv_nst < ca_mem->ca_nsteps, keeping this for + future regression */ + throw IDAException(AMICI_ERROR, "reInitPostProcess"); + } + + /* Load next point in dt_mem */ + dt_mem[ida_mem->ida_nst % ia_mem->ia_nsteps]->t = *t; + ia_mem->ia_storePnt(ida_mem, + dt_mem[ida_mem->ida_nst % ia_mem->ia_nsteps]); + + /* Set t1 field of the current ckeck point structure + for the case in which there will be no future + check points */ + ia_mem->ck_mem->ck_t1 = *t; + + /* tfinal is now set to *tret */ + ia_mem->ia_tfinal = *t; + } +} void IDASolver::reInit(realtype t0, AmiVector *yy0, AmiVector *yp0) { auto ida_mem = static_cast(solverMemory.get()); diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 44c06e2835..ae006e60d2 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -295,6 +295,7 @@ std::unique_ptr SteadystateProblem::createSteadystateSimSolver( switch(solver->getLinearSolver()) { case LinearSolver::dense: case LinearSolver::KLU: + case LinearSolver::SuperLUMT: break; default: throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "invalid solver for steadystate simulation"); diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp index e2c4c1669e..ea0b3ac8a2 100644 --- a/src/sundials_linsol_wrapper.cpp +++ b/src/sundials_linsol_wrapper.cpp @@ -1,6 +1,7 @@ -#include #include +#include + #include // bad_alloc #include @@ -367,4 +368,57 @@ int SUNNonLinSolFixedPoint::getSysFn(SUNNonlinSolSysFn *SysFn) const { return SUNNonlinSolGetSysFn_FixedPoint(solver, SysFn); } +#ifdef SUNDIALS_SUPERLUMT + +SUNLinSolSuperLUMT::SUNLinSolSuperLUMT(N_Vector x, SUNMatrix A, int numThreads) + : SUNLinSolWrapper(SUNLinSol_SuperLUMT(x, A, numThreads)) +{ + if (!solver) + throw AmiException("Failed to create solver."); +} + +SUNLinSolSuperLUMT::SUNLinSolSuperLUMT( + const AmiVector &x, int nnz, int sparsetype, + SUNLinSolSuperLUMT::StateOrdering ordering) + : A(SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype)) +{ + int numThreads = 1; + if(auto env = std::getenv("AMICI_SUPERLUMT_NUM_THREADS")) { + numThreads = std::max(1, std::stoi(env)); + } + + solver = SUNLinSol_SuperLUMT(x.getNVector(), A.get(), numThreads); + if (!solver) + throw AmiException("Failed to create solver."); + + setOrdering(ordering); +} + +SUNLinSolSuperLUMT::SUNLinSolSuperLUMT(const AmiVector &x, int nnz, + int sparsetype, StateOrdering ordering, + int numThreads) + : A(SUNMatrixWrapper(x.getLength(), x.getLength(), nnz, sparsetype)) +{ + solver = SUNLinSol_SuperLUMT(x.getNVector(), A.get(), numThreads); + if (!solver) + throw AmiException("Failed to create solver."); + + setOrdering(ordering); +} + +SUNMatrix SUNLinSolSuperLUMT::getMatrix() const +{ + return A.get(); +} + + +void SUNLinSolSuperLUMT::setOrdering(StateOrdering ordering) +{ + auto status = SUNLinSol_SuperLUMTSetOrdering(solver, static_cast(ordering)); + if (status != SUNLS_SUCCESS) + throw AmiException("SUNLinSol_SuperLUMTSetOrdering failed with %d", status); +} + +#endif + } // namespace amici diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp index 9b2275fac9..1f3667c754 100644 --- a/src/sundials_matrix_wrapper.cpp +++ b/src/sundials_matrix_wrapper.cpp @@ -10,13 +10,14 @@ namespace amici { SUNMatrixWrapper::SUNMatrixWrapper(int M, int N, int NNZ, int sparsetype) : matrix(SUNSparseMatrix(M, N, NNZ, sparsetype)) { + if (sparsetype != CSC_MAT && sparsetype != CSR_MAT) throw std::invalid_argument("Invalid sparsetype. Must be CSC_MAT or " "CSR_MAT"); if (NNZ && !matrix) throw std::bad_alloc(); - + update_ptrs(); } @@ -24,7 +25,7 @@ SUNMatrixWrapper::SUNMatrixWrapper(int M, int N) : matrix(SUNDenseMatrix(M, N)) { if (M && N && !matrix) throw std::bad_alloc(); - + update_ptrs(); } @@ -32,7 +33,7 @@ SUNMatrixWrapper::SUNMatrixWrapper(int M, int ubw, int lbw) : matrix(SUNBandMatrix(M, ubw, lbw)) { if (M && !matrix) throw std::bad_alloc(); - + update_ptrs(); } @@ -159,39 +160,31 @@ void SUNMatrixWrapper::reset() { SUNMatZero(matrix); } -void SUNMatrixWrapper::multiply(std::vector &c, - const std::vector &b) const { - if (static_cast(c.size()) != rows()) - throw std::invalid_argument("Dimension mismatch between number of rows" - "in A and elements in c"); - - if (static_cast(b.size()) != columns()) - throw std::invalid_argument("Dimension mismatch between number of cols" - "in A and elements in b"); - - multiply(c.data(), b.data()); -} - void SUNMatrixWrapper::multiply(N_Vector c, const N_Vector b) const { - if (NV_LENGTH_S(c) != rows()) - throw std::invalid_argument("Dimension mismatch between number of rows" - "in A and elements in c"); - - if (NV_LENGTH_S(b) != columns()) - throw std::invalid_argument("Dimension mismatch between number of cols" - "in A and elements in b"); - - multiply(NV_DATA_S(c), NV_DATA_S(b)); + multiply(gsl::make_span(NV_DATA_S(c), NV_LENGTH_S(c)), + gsl::make_span(NV_DATA_S(b), NV_LENGTH_S(b))); } -void SUNMatrixWrapper::multiply(realtype *c, const realtype *b) const { +void SUNMatrixWrapper::multiply(gsl::span c, gsl::span b) const { if (!matrix) return; + if (static_cast(c.size()) != rows()) + throw std::invalid_argument("Dimension mismatch between number of rows " + "in A (" + std::to_string(rows()) + ") and " + "elements in c (" + std::to_string(c.size()) + + ")"); + + if (static_cast(b.size()) != columns()) + throw std::invalid_argument("Dimension mismatch between number of cols " + "in A (" + std::to_string(columns()) + + ") and elements in b (" + + std::to_string(b.size()) + ")"); + switch (SUNMatGetID(matrix)) { case SUNMATRIX_DENSE: amici_dgemv(BLASLayout::colMajor, BLASTranspose::noTrans, rows(), - columns(), 1.0, data(), rows(), b, 1, 1.0, c, 1); + columns(), 1.0, data(), rows(), b.data(), 1, 1.0, c.data(), 1); break; case SUNMATRIX_SPARSE: @@ -222,10 +215,17 @@ void SUNMatrixWrapper::multiply(realtype *c, const realtype *b) const { } } +void SUNMatrixWrapper::zero() +{ + if(int res = SUNMatZero(matrix)) + throw std::runtime_error("SUNMatrixWrapper::zero() failed with " + + std::to_string(res)); +} + void SUNMatrixWrapper::update_ptrs() { if(!matrix) return; - + switch (SUNMatGetID(matrix)) { case SUNMATRIX_DENSE: if (columns() > 0 && rows() > 0) diff --git a/tests/cpputest/CMakeLists.txt b/tests/cpputest/CMakeLists.txt index baac740315..2553825764 100644 --- a/tests/cpputest/CMakeLists.txt +++ b/tests/cpputest/CMakeLists.txt @@ -37,6 +37,7 @@ set(TEST_MODELS events nested_events robertson + calvetti ) if(ENABLE_SWIG AND ENABLE_PYTHON) diff --git a/tests/cpputest/calvetti/CMakeLists.txt b/tests/cpputest/calvetti/CMakeLists.txt new file mode 100644 index 0000000000..388fdf9152 --- /dev/null +++ b/tests/cpputest/calvetti/CMakeLists.txt @@ -0,0 +1,19 @@ +get_filename_component(MODEL_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) +project(model_${MODEL_NAME}_test) + +set(SRC_LIST + ../main.cpp + tests1.cpp +) + +include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CppUTest_INCLUDE_DIRS}) + +add_executable(${PROJECT_NAME} ${SRC_LIST}) +add_dependencies(${PROJECT_NAME} external_model_${MODEL_NAME}) + +target_link_libraries(${PROJECT_NAME} + amici-testing + model_${MODEL_NAME} +) + +add_test(NAME ${PROJECT_NAME} COMMAND ./${PROJECT_NAME} -c) diff --git a/tests/cpputest/calvetti/tests1.cpp b/tests/cpputest/calvetti/tests1.cpp new file mode 100644 index 0000000000..ae9b28c3cf --- /dev/null +++ b/tests/cpputest/calvetti/tests1.cpp @@ -0,0 +1,26 @@ +#include "CppUTest/TestHarness.h" +#include "CppUTestExt/MockSupport.h" + +#include +#include +#include "wrapfunctions.h" + +TEST_GROUP(groupCalvetti) +{ + void setup() { + + } + + void teardown() { + + } +}; + +TEST(groupCalvetti, testSimulation) { + amici::simulateVerifyWrite("/model_calvetti/nosensi/"); +} + + + + + diff --git a/tests/cpputest/expectedResults.h5 b/tests/cpputest/expectedResults.h5 index 85031cc85b..87cbc6e87a 100644 Binary files a/tests/cpputest/expectedResults.h5 and b/tests/cpputest/expectedResults.h5 differ diff --git a/tests/cpputest/jakstat_adjoint/tests1.cpp b/tests/cpputest/jakstat_adjoint/tests1.cpp index 1948ace88d..7d3104b28d 100644 --- a/tests/cpputest/jakstat_adjoint/tests1.cpp +++ b/tests/cpputest/jakstat_adjoint/tests1.cpp @@ -41,3 +41,38 @@ TEST(groupJakstatAdjoint, testSensitivityAdjointEmptySensInd) { amici::simulateVerifyWrite("/model_jakstat_adjoint/sensiadjointemptysensind/"); } +IGNORE_TEST(groupJakstatAdjoint, testSensitivityAdjointUnusedNanOutputs) { + /* UN-IGNORE ONCE THIS MODEL HAS BEEN IMPORTED VIA PYTHON INTERFACE */ + auto model = getModel(); + auto solver = model->getSolver(); + amici::hdf5::readModelDataFromHDF5( + NEW_OPTION_FILE, *model, + "/model_jakstat_adjoint/sensiadjoint/options"); + amici::hdf5::readSolverSettingsFromHDF5( + NEW_OPTION_FILE, *solver, + "/model_jakstat_adjoint/sensiadjoint/options"); + auto edata = amici::hdf5::readSimulationExpData( + NEW_OPTION_FILE, "/model_jakstat_adjoint/sensiadjoint/data", + *model); + + // Set output parameter p[10] to NaN and remove respective measurements + // -> gradient should still be finite + + auto p = model->getParameters(); + p[10] = NAN; + model->setParameters(p); + + auto d = edata->getObservedData(); + for (int it = 0; it < edata->nt(); ++it) { + for (int iy = 0; iy < edata->nytrue(); ++iy) { + if(iy == 1) + d[it * edata->nytrue() + it] = NAN; + } + } + edata->setObservedData(d); + + auto rdata = runAmiciSimulation(*solver, edata.get(), *model); + + for(int i = 0; i < model->nplist(); ++i) + CHECK_FALSE(std::isnan(rdata->sllh[i])); +} diff --git a/tests/cpputest/robertson/tests1.cpp b/tests/cpputest/robertson/tests1.cpp index 2f170cf7d3..5e31c36a88 100644 --- a/tests/cpputest/robertson/tests1.cpp +++ b/tests/cpputest/robertson/tests1.cpp @@ -24,6 +24,16 @@ TEST(groupRobertson, testSensitivityForward) { amici::simulateVerifyWrite("/model_robertson/sensiforward/", 1e6*TEST_ATOL, 1e2*TEST_RTOL); } +TEST(groupRobertson, testSensitivityForwardDense) { + amici::simulateVerifyWrite("/model_robertson/sensiforwarddense/", 1e6*TEST_ATOL, 1e2*TEST_RTOL); +} + +TEST(groupRobertson, testSensitivityForwardSPBCG) { + amici::simulateVerifyWrite("/model_robertson/sensiforwardSPBCG/", 1e7*TEST_ATOL, 1e2*TEST_RTOL); +} + + + diff --git a/tests/cpputest/testOptions.h5 b/tests/cpputest/testOptions.h5 index 62e3b78e9e..a1e9a43c4e 100644 Binary files a/tests/cpputest/testOptions.h5 and b/tests/cpputest/testOptions.h5 differ diff --git a/tests/cpputest/testfunctions.cpp b/tests/cpputest/testfunctions.cpp index 2c6fbb639e..c2906e757d 100644 --- a/tests/cpputest/testfunctions.cpp +++ b/tests/cpputest/testfunctions.cpp @@ -89,7 +89,6 @@ void simulateVerifyWrite(const std::string& hdffileOptions, const std::string& h verifyReturnData(hdffileResults, resultPath, rdata_resimulation.get(), model.get(), atol, rtol); // verify written results verifyReturnData(hdffilewrite, writePath, rdata.get(), model.get(), atol, rtol); - //remove(hdffilewrite.c_str()); } std::unique_ptr getTestExpData(Model const& model) { @@ -156,12 +155,21 @@ void checkEqualArrayStrided(const double *expected, const double *actual, int le void verifyReturnData(std::string const& hdffile, std::string const& resultPath, const ReturnData *rdata, const Model *model, double atol, double rtol) { + CHECK_FALSE(rdata == nullptr); + + if(!hdf5::locationExists(hdffile, resultPath)) { + fprintf(stderr, "ERROR: No results available for %s!\n", + resultPath.c_str()); + return; + } // compare to saved data in hdf file H5::H5File file(hdffile, H5F_ACC_RDONLY); hsize_t m, n; + + std::vector expected; auto statusExp = hdf5::getIntScalarAttribute(file, resultPath, "status"); CHECK_EQUAL(statusExp, rdata->status); @@ -172,16 +180,28 @@ void verifyReturnData(std::string const& hdffile, std::string const& resultPath, double chi2Exp = hdf5::getDoubleScalarAttribute(file, resultPath, "chi2"); CHECK_TRUE(withinTolerance(chi2Exp, rdata->chi2, atol, rtol, 1, "chi2")); - auto expected = hdf5::getDoubleDataset2D(file, resultPath + "/x", m, n); - checkEqualArray(expected, rdata->x, atol, rtol, "x"); + if(hdf5::locationExists(file, resultPath + "/x")) { + expected = hdf5::getDoubleDataset2D(file, resultPath + "/x", m, n); + checkEqualArray(expected, rdata->x, atol, rtol, "x"); + } else { + CHECK_TRUE(rdata->x.empty()); + } // CHECK_EQUAL(AMICI_O2MODE_FULL, udata->o2mode); + + if(hdf5::locationExists(file, resultPath + "/diagnosis/J")) { + expected = hdf5::getDoubleDataset2D(file, resultPath + "/diagnosis/J", m, n); + checkEqualArray(expected, rdata->J, atol, rtol, "J"); + } else { + CHECK_TRUE(rdata->J.empty()); + } - expected = hdf5::getDoubleDataset2D(file, resultPath + "/diagnosis/J", m, n); - checkEqualArray(expected, rdata->J, atol, rtol, "J"); - - expected = hdf5::getDoubleDataset2D(file, resultPath + "/y", m, n); - checkEqualArray(expected, rdata->y, atol, rtol, "y"); + if(hdf5::locationExists(file, resultPath + "/y")) { + expected = hdf5::getDoubleDataset2D(file, resultPath + "/y", m, n); + checkEqualArray(expected, rdata->y, atol, rtol, "y"); + } else { + CHECK_TRUE(rdata->y.empty()); + } if(hdf5::locationExists(file, resultPath + "/res")) { expected = hdf5::getDoubleDataset1D(file, resultPath + "/res"); diff --git a/tests/cpputest/testfunctions.h b/tests/cpputest/testfunctions.h index fbb28b6251..32a70252b4 100644 --- a/tests/cpputest/testfunctions.h +++ b/tests/cpputest/testfunctions.h @@ -73,12 +73,12 @@ class Model_Test : public Model { const std::vector k, const std::vector plist, const std::vector idlist, const std::vector z2event) : Model(nx_rdata, nxtrue_rdata, nx_solver, nxtrue_solver, ny, nytrue, nz, - nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, nnz, ubw, lbw, o2mode, + nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, {}, nnz, ubw, lbw, o2mode, p, k, plist, idlist, z2event) {} /** default constructor */ Model_Test() - : Model(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + : Model(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, {}, 0, 0, 0, SecondOrderMode::none, std::vector(), std::vector(), std::vector(), std::vector(), std::vector()) {} diff --git a/tests/cpputest/unittests/testsSerialization.cpp b/tests/cpputest/unittests/testsSerialization.cpp index 0cfbe9328a..aae97bdce8 100644 --- a/tests/cpputest/unittests/testsSerialization.cpp +++ b/tests/cpputest/unittests/testsSerialization.cpp @@ -95,7 +95,7 @@ TEST_GROUP(dataSerialization){ solver.setNewtonMaxSteps(1e6); solver.setNewtonMaxLinearSteps(1e6); solver.setNewtonPreequilibration(true); - solver.setStateOrdering(amici::StateOrdering::COLAMD); + solver.setStateOrdering(static_cast(amici::SUNLinSolKLU::StateOrdering::COLAMD)); solver.setInterpolationType(amici::InterpolationType::polynomial); solver.setStabilityLimitFlag(0); solver.setLinearSolver(amici::LinearSolver::dense); @@ -157,11 +157,11 @@ TEST(dataSerialization, testString) { std::vector(np,0), std::vector(nx,0.0), std::vector(nz,0)); - + amici::ReturnData r(solver, &m); - + std::string serialized = amici::serializeToString(r); - + checkReturnDataEqual(r, amici::deserializeFromString(serialized)); } diff --git a/tests/cpputest/wrapTestModels.m b/tests/cpputest/wrapTestModels.m index 6417c626c2..28963e2e3a 100644 --- a/tests/cpputest/wrapTestModels.m +++ b/tests/cpputest/wrapTestModels.m @@ -92,5 +92,18 @@ function wrapTestModels() cd(fileparts(mfilename('fullpath'))); + %% EXAMPLE CALVETTI + cd([amiciPath '/examples/example_calvetti/']); + + try + [exdir,~,~]=fileparts(which('example_calvetti.m')); + amiwrap('model_calvetti', 'model_calvetti_syms', exdir); + catch err + disp(err.message) + cd(fileparts(mfilename('fullpath'))); + end + + cd(fileparts(mfilename('fullpath'))); + end diff --git a/tests/generateTestConfig/example_calvetti.py b/tests/generateTestConfig/example_calvetti.py new file mode 100755 index 0000000000..86178d684b --- /dev/null +++ b/tests/generateTestConfig/example_calvetti.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python3 +import sys +import numpy as np +from example import AmiciExample + + +class ExampleCalvetti(AmiciExample): + + def __init__(self): + AmiciExample.__init__( self ) + + self.numX = 6 + self.numP = 0 + self.numK = 6 + + self.modelOptions['theta'] = [] + self.modelOptions['kappa'] = [0.29, 0.74, 0.44, 0.08, 0.27, 0.18] + self.modelOptions['ts'] = np.linspace(0, 20, 201) + self.modelOptions['pscale'] = 0 + + self.solverOptions['atol'] = 1e-6 + self.solverOptions['rtol'] = 1e-4 + self.solverOptions['sens_ind'] = [] + self.solverOptions['sensi'] = 0 + self.solverOptions['sensi_meth'] = 1 + + +def writeNoSensi(filename): + ex = ExampleCalvetti() + + ex.writeToFile(filename, '/model_calvetti/nosensi/') + + +def main(): + if len(sys.argv) < 2: + print("Error: Must provide output file as first and only argument.") + sys.exit(1) + filename = sys.argv[1] + + writeNoSensi(filename) + + +if __name__ == "__main__": + main() diff --git a/tests/generateTestConfig/example_jakstat.py b/tests/generateTestConfig/example_jakstat.py index 381f8891cb..d914192732 100755 --- a/tests/generateTestConfig/example_jakstat.py +++ b/tests/generateTestConfig/example_jakstat.py @@ -27,7 +27,7 @@ def __init__(self): self.modelOptions['kappa'] = [1.4, 0.45] self.modelOptions['pscale'] = 2 - self.solverOptions['atol'] = 1e-12 + self.solverOptions['atol'] = 1e-16 self.solverOptions['maxsteps'] = 1e4 self.solverOptions['nmaxevent'] = 10 self.solverOptions['rtol'] = 1e-12 diff --git a/tests/generateTestConfig/example_robertson.py b/tests/generateTestConfig/example_robertson.py index a6a180454e..a9a71d4db1 100755 --- a/tests/generateTestConfig/example_robertson.py +++ b/tests/generateTestConfig/example_robertson.py @@ -21,8 +21,8 @@ def __init__(self): self.modelOptions['ts'] = [0.0, 4.0E-6, 7.030042499419172E-6, 1.2355374385909914E-5, 2.1714701757295438E-5, 3.8163819053999774E-5, 6.70733174724404E-5, 1.1788206810207238E-4, 2.071789871692485E-4, 3.641192711966091E-4, 6.39943487842423E-4, 0.0011247074791896924, 0.001976685344529533, 0.003474045495005412, 0.006105671868700933, 0.010730783181118898, 0.018859465453829577, 0.03314571091418737, 0.05825393910004978, 0.10238191690798133, 0.17993730675877775, 0.31624172843630804, 0.5557981977492555, 0.97682123781946, 1.7167737040515112, 3.017248025341849, 5.302845462360432, 9.319807242061488, 16.37966024952171, 28.787426920046055, 50.59420867421183, 88.91985930104782, 156.27759748218483, 274.6595380017199, 482.71705625573054, 848.380355168077, 1491.0374881259752, 2620.5142274381983, 4605.5815973057925, 8094.358590900622, 14225.921224892543, 25002.207701095904, 43941.64567950229, 77227.90915533014, 135728.87087581318, 238544.93266378547, 419245.253661875, 736827.9877306866, 1294983.017127056, 2275946.411607322, 4000000.0] self.modelOptions['pscale'] = 2 - self.solverOptions['atol'] = 1e-6 - self.solverOptions['rtol'] = 1e-4 + self.solverOptions['atol'] = 1e-12 + self.solverOptions['rtol'] = 1e-8 self.solverOptions['sens_ind'] = [] self.solverOptions['sensi'] = 0 self.solverOptions['sensi_meth'] = 1 @@ -43,6 +43,28 @@ def writeSensiForward(filename): ex.writeToFile(filename, '/model_robertson/sensiforward/') +def writeSensiForwardDense(filename): + ex = ExampleRobertson() + + ex.solverOptions['sens_ind'] = np.arange(0, ex.numP) + ex.solverOptions['sensi'] = 1 + ex.solverOptions['linsol'] = 1 + + ex.writeToFile(filename, '/model_robertson/sensiforwarddense/') + + +def writeSensiForwardSPBCG(filename): + ex = ExampleRobertson() + + ex.solverOptions['sens_ind'] = np.arange(0, ex.numP) + ex.solverOptions['sensi'] = 1 + ex.solverOptions['linsol'] = 7 + ex.solverOptions['atol'] = 1e-14 + ex.solverOptions['rtol'] = 1e-12 + + ex.writeToFile(filename, '/model_robertson/sensiforwardSPBCG/') + + def main(): if len(sys.argv) < 2: print("Error: Must provide output file as first and only argument.") @@ -51,6 +73,8 @@ def main(): writeNoSensi(filename) writeSensiForward(filename) + writeSensiForwardDense(filename) + writeSensiForwardSPBCG(filename) if __name__ == "__main__": diff --git a/tests/generateTestConfigurationForExamples.sh b/tests/generateTestConfigurationForExamples.sh index c30a2fdd4d..578228ee31 100755 --- a/tests/generateTestConfigurationForExamples.sh +++ b/tests/generateTestConfigurationForExamples.sh @@ -19,3 +19,4 @@ cd ${AMICI_PATH}/tests/generateTestConfig ./example_neuron.py ${TEST_FILE} ./example_robertson.py ${TEST_FILE} ./example_steadystate.py ${TEST_FILE} +./example_calvetti.py ${TEST_FILE} diff --git a/tests/testMisc.py b/tests/testMisc.py index 62ff341118..7157bae01a 100755 --- a/tests/testMisc.py +++ b/tests/testMisc.py @@ -7,7 +7,7 @@ import copy import os import unittest - +import sympy as sp class TestAmiciMisc(unittest.TestCase): """TestCase class various AMICI Python interface functions""" @@ -35,6 +35,41 @@ def test_parameterScalingFromIntVector(self): assert scale_vector[1] == amici.ParameterScaling_ln assert scale_vector[2] == amici.ParameterScaling_none + def test_csc_matrix(self): + """Test sparse CSC matrix creation""" + matrix = sp.Matrix([[1, 0], [2, 3]]) + symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix = \ + amici.ode_export.csc_matrix(matrix, 'a') + + assert symbolColPtrs == [0, 2, 3] + assert symbolRowVals == [0, 1, 1] + assert sparseList == sp.Matrix([[1], [2], [3]]) + assert symbolList == ['a0', 'a1', 'a2'] + assert str(sparseMatrix) == 'Matrix([[a0, 0], [a1, a2]])' + + def test_csc_matrix_vector(self): + """Test sparse CSC matrix creation from matrix slice""" + matrix = sp.Matrix([[1, 0], [2, 3]]) + symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix = \ + amici.ode_export.csc_matrix(matrix[:, 0], 'a') + + assert symbolColPtrs == [0, 2] + assert symbolRowVals == [0, 1] + assert sparseList == sp.Matrix([[1], [2]]) + assert symbolList == ['a0', 'a1'] + assert str(sparseMatrix) == 'Matrix([[a0], [a1]])' + + '''Test continuation of numbering of symbols''' + symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix = \ + amici.ode_export.csc_matrix(matrix[:, 1], 'a', + base_index=len(symbolList)) + + assert symbolColPtrs == [0, 1] + assert symbolRowVals == [1] + assert sparseList == sp.Matrix([[3]]) + assert symbolList == ['a2'] + assert str(sparseMatrix) == 'Matrix([[0], [a2]])' + if __name__ == '__main__': suite = unittest.TestSuite() diff --git a/tests/testModels.py b/tests/testModels.py index eb12f2e986..bf11007d46 100755 --- a/tests/testModels.py +++ b/tests/testModels.py @@ -98,6 +98,13 @@ def assert_fun(x): assert_fun, atol=1e-5, rtol=1e-3 ) + elif model_name == 'model_robertson' and \ + case == 'sensiforwardSPBCG': + verify_simulation_results( + rdata, expected_results[subTest][case]['results'], + assert_fun, + atol=1e-3, rtol=1e-3 + ) else: verify_simulation_results( rdata, expected_results[subTest][case]['results'], @@ -108,7 +115,10 @@ def assert_fun(x): case == 'sensiforwarderrorint': edata = amici.amici.ExpData(self.model.get()) - if edata and model_name != 'model_neuron_o2': + if edata and model_name != 'model_neuron_o2' and not ( + model_name == 'model_robertson' and + case == 'sensiforwardSPBCG' + ): # Test runAmiciSimulations: ensure running twice # with same ExpData yields same results if isinstance(edata, amici.amici.ExpData): @@ -159,7 +169,7 @@ def verify_simulation_results(rdata, expected_results, assert_fun, for subfield in ['J', 'xdot']: check_results(rdata, subfield, expected_results[field][subfield][()], - assert_fun, 0, 2) + assert_fun, 1e-8, 1e-8) else: if field == 's2llh': check_results(rdata, field, expected_results[field][()], diff --git a/version.txt b/version.txt index 9028ec6365..69da6ebcd0 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.10.5 +0.10.6