diff --git a/CMakeLists.txt b/CMakeLists.txt index 52ee00c701..3e2a7f1f95 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -135,6 +135,7 @@ set(AMICI_SRC_LIST ${CMAKE_SOURCE_DIR}/src/model.cpp ${CMAKE_SOURCE_DIR}/src/model_ode.cpp ${CMAKE_SOURCE_DIR}/src/model_dae.cpp + ${CMAKE_SOURCE_DIR}/src/model_state.cpp ${CMAKE_SOURCE_DIR}/src/newton_solver.cpp ${CMAKE_SOURCE_DIR}/src/forwardproblem.cpp ${CMAKE_SOURCE_DIR}/src/steadystateproblem.cpp diff --git a/include/amici/model.h b/include/amici/model.h index 0a8ad94d78..a450a0b030 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -6,6 +6,8 @@ #include "amici/sundials_matrix_wrapper.h" #include "amici/vector.h" #include "amici/simulation_parameters.h" +#include "amici/model_dimensions.h" +#include "amici/model_state.h" #include #include @@ -31,437 +33,6 @@ void serialize(Archive &ar, amici::Model &m, unsigned int version); } // namespace boost namespace amici { -/** - * @brief Container for model dimensions. - * - * Holds number of states, observables, etc. - */ -struct ModelDimensions { - /** Default ctor */ - ModelDimensions() = default; - - /** - * @brief Constructor with model dimensions - * @param nx_rdata Number of state variables - * @param nxtrue_rdata Number of state variables of the non-augmented model - * @param nx_solver Number of state variables with conservation laws applied - * @param nxtrue_solver Number of state variables of the non-augmented model - * with conservation laws applied - * @param nx_solver_reinit Number of state variables with conservation laws - * subject to reinitialization - * @param np Number of parameters - * @param nk Number of constants - * @param ny Number of observables - * @param nytrue Number of observables of the non-augmented model - * @param nz Number of event observables - * @param nztrue Number of event observables of the non-augmented model - * @param ne Number of events - * @param nJ Number of objective functions - * @param nw Number of repeating elements - * @param ndwdx Number of nonzero elements in the `x` derivative of the - * repeating elements - * @param ndwdp Number of nonzero elements in the `p` derivative of the - * repeating elements - * @param ndwdw Number of nonzero elements in the `w` derivative of the - * repeating elements - * @param ndxdotdw Number of nonzero elements in the \f$w\f$ derivative of - * \f$xdot\f$ - * @param ndJydy Number of nonzero elements in the \f$y\f$ derivative of - * \f$dJy\f$ (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 - */ - ModelDimensions( - const int nx_rdata, const int nxtrue_rdata, const int nx_solver, - const int nxtrue_solver, const int nx_solver_reinit, const int np, - const int nk, 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 ndwdw, const int ndxdotdw, std::vector ndJydy, - const int nnz, const int ubw, const int lbw) - : nx_rdata(nx_rdata), nxtrue_rdata(nxtrue_rdata), nx_solver(nx_solver), - nxtrue_solver(nxtrue_solver), nx_solver_reinit(nx_solver_reinit), - np(np), nk(nk), - ny(ny), nytrue(nytrue), nz(nz), nztrue(nztrue), - ne(ne), nw(nw), ndwdx(ndwdx), ndwdp(ndwdp), ndwdw(ndwdw), - ndxdotdw(ndxdotdw), ndJydy(std::move(ndJydy)), - nnz(nnz), nJ(nJ), ubw(ubw), lbw(lbw) { - Expects(nxtrue_rdata >= 0); - Expects(nxtrue_rdata <= nx_rdata); - Expects(nxtrue_solver >= 0); - Expects(nx_solver <= nx_rdata); - Expects(nxtrue_solver <= nx_solver); - Expects(nx_solver_reinit >= 0); - Expects(nx_solver_reinit <= nx_solver); - Expects(np >= 0); - Expects(nk >= 0); - Expects(nytrue <= ny); - Expects(nytrue >= 0); - Expects(nztrue >= 0); - Expects(nztrue <= nz); - Expects(ne >= 0); - Expects(nw >= 0); - Expects(ndwdx >= 0); - Expects(ndwdx <= nw * nx_solver); - Expects(ndwdp >= 0); - Expects(ndwdp <= nw * np); - Expects(ndwdw >= 0); - Expects(ndwdw <= nw * nw); - Expects(ndxdotdw >= 0); - Expects(nnz >= 0); - Expects(nJ >= 0); - Expects(ubw >= 0); - Expects(lbw >= 0); - } - - /** Number of states */ - int nx_rdata{0}; - - /** Number of states in the unaugmented system */ - int nxtrue_rdata{0}; - - /** Number of states with conservation laws applied */ - int nx_solver{0}; - - /** - * Number of states in the unaugmented system with conservation laws - * applied - */ - int nxtrue_solver{0}; - - /** Number of solver states subject to reinitialization */ - int nx_solver_reinit{0}; - - /** Number of parameters */ - int np{0}; - - /** Number of constants */ - int nk{0}; - - /** Number of observables */ - int ny{0}; - - /** Number of observables in the unaugmented system */ - int nytrue{0}; - - /** Number of event outputs */ - int nz{0}; - - /** Number of event outputs in the unaugmented system */ - int nztrue{0}; - - /** Number of events */ - int ne{0}; - - /** Number of common expressions */ - int nw{0}; - - /** - * Number of nonzero elements in the `x` derivative of the - * repeating elements - */ - int ndwdx {0}; - - /** - * Number of nonzero elements in the `p` derivative of the - * repeating elements - */ - int ndwdp {0}; - - /** - * Number of nonzero elements in the `w` derivative of the - * repeating elements - */ - int ndwdw {0}; - - /** Number of nonzero elements in the \f$w\f$ derivative of \f$xdot\f$ */ - int ndxdotdw {0}; - - /** - * Number of nonzero elements in the \f$y\f$ derivative of - * \f$dJy\f$ (dimension `nytrue`) - */ - std::vector ndJydy; - - /** Number of nonzero entries in Jacobian */ - int nnz{0}; - - /** Dimension of the augmented objective function for 2nd order ASA */ - int nJ{0}; - - /** Upper bandwidth of the Jacobian */ - int ubw{0}; - - /** Lower bandwidth of the Jacobian */ - int lbw{0}; -}; - -/** - * @brief Exchange format to store and transfer the state of the - * model at a specific timepoint. - * - * This is designed to only encompass the minimal - * number of attributes that need to be transferred. - */ -struct ModelState { - /** - * Flag indicating whether a certain Heaviside function should be active or - * not (dimension: `ne`) - */ - std::vector h; - - /** Total abundances for conservation laws - (dimension: `nx_rdata - nx_solver`) */ - std::vector total_cl; - - /** Sensitivities of total abundances for conservation laws - (dimension: `(nx_rdata-nx_solver) x np`, row-major) */ - std::vector stotal_cl; - - /** Unscaled parameters (dimension: `np`) */ - std::vector unscaledParameters; - - /** Constants (dimension: `nk`) */ - std::vector fixedParameters; - - /** - * Indexes of parameters wrt to which sensitivities are computed - * (dimension: nplist) - */ - std::vector plist; -}; - - -/** - * @brief Storage for `amici::Model` quantities computed based on - * `amici::ModelState` for a specific timepoint. - * - * Serves as workspace for a model simulation to avoid repeated reallocation. - */ -struct ModelStateDerived { - ModelStateDerived() = default; - - /** - * @brief Constructor from model dimensions. - * @param dim Model dimensions - */ - explicit ModelStateDerived(ModelDimensions const& dim) - : J_(dim.nx_solver, dim.nx_solver, dim.nnz, CSC_MAT), - JB_(dim.nx_solver, dim.nx_solver, dim.nnz, CSC_MAT), - dxdotdw_(dim.nx_solver, dim.nw, dim.ndxdotdw, CSC_MAT), - w_(dim.nw), - x_rdata_(dim.nx_rdata, 0.0), - sx_rdata_(dim.nx_rdata, 0.0), - x_pos_tmp_(dim.nx_solver) - {} - - /** Sparse Jacobian (dimension: `amici::Model::nnz`) */ - SUNMatrixWrapper J_; - - /** Sparse Backwards Jacobian (dimension: `amici::Model::nnz`) */ - SUNMatrixWrapper JB_; - - /** Sparse dxdotdw temporary storage (dimension: `ndxdotdw`) */ - SUNMatrixWrapper dxdotdw_; - - /** Sparse dwdx temporary storage (dimension: `ndwdx`) */ - SUNMatrixWrapper dwdx_; - - /** Sparse dwdp temporary storage (dimension: `ndwdp`) */ - SUNMatrixWrapper dwdp_; - - /** Dense Mass matrix (dimension: `nx_solver` x `nx_solver`) */ - SUNMatrixWrapper M_; - - /** - * Temporary storage of `dxdotdp_full` data across functions (Python only) - * (dimension: `nplist` x `nx_solver`, nnz: dynamic, - * type `CSC_MAT`) - */ - SUNMatrixWrapper dxdotdp_full; - - /** - * Temporary storage of `dxdotdp_explicit` data across functions (Python only) - * (dimension: `nplist` x `nx_solver`, nnz: `ndxdotdp_explicit`, - * type `CSC_MAT`) - */ - SUNMatrixWrapper dxdotdp_explicit; - - /** - * Temporary storage of `dxdotdp_implicit` data across functions, - * Python-only - * (dimension: `nplist` x `nx_solver`, nnz: dynamic, - * type `CSC_MAT`) - */ - SUNMatrixWrapper dxdotdp_implicit; - - /** - * Temporary storage of `dxdotdx_explicit` data across functions (Python only) - * (dimension: `nplist` x `nx_solver`, nnz: `nxdotdotdx_explicit`, - * type `CSC_MAT`) - */ - SUNMatrixWrapper dxdotdx_explicit; - - /** - * Temporary storage of `dxdotdx_implicit` data across functions, - * Python-only - * (dimension: `nplist` x `nx_solver`, nnz: dynamic, - * type `CSC_MAT`) - */ - SUNMatrixWrapper dxdotdx_implicit; - - /** - * Temporary storage of `dxdotdp` data across functions, Matlab only - * (dimension: `nplist` x `nx_solver` , row-major) - */ - AmiVectorArray dxdotdp {0, 0}; - - /** Sparse observable derivative of data likelihood, only used if - * `pythonGenerated` == `true` (dimension `nytrue`, `nJ` x `ny`, row-major) - */ - std::vector dJydy_; - - /** Observable derivative of data likelihood, only used if - * `pythonGenerated` == `false` (dimension `nJ` x `ny` x `nytrue` , - * row-major) - */ - std::vector dJydy_matlab_; - - /** Observable sigma derivative of data likelihood - * (dimension nJ x ny x nytrue, row-major) - */ - std::vector dJydsigma_; - - /** State derivative of data likelihood - * (dimension `nJ` x `nx_solver`, row-major) - */ - std::vector dJydx_; - - /** Parameter derivative of data likelihood for current timepoint - * (dimension: nJ x nplist, row-major) - */ - std::vector dJydp_; - - /** event output derivative of event likelihood - * (dimension nJ x nz x nztrue, row-major) - */ - std::vector dJzdz_; - - /** event sigma derivative of event likelihood - * (dimension nJ x nz x nztrue, row-major) - */ - std::vector dJzdsigma_; - - /** event output derivative of event likelihood at final timepoint - * (dimension nJ x nz x nztrue, row-major) - */ - std::vector dJrzdz_; - - /** event sigma derivative of event likelihood at final timepoint - * (dimension nJ x nz x nztrue, row-major) - */ - std::vector dJrzdsigma_; - - /** state derivative of event likelihood - * (dimension `nJ` x `nx_solver`, row-major) - */ - std::vector dJzdx_; - - /** parameter derivative of event likelihood for current timepoint - * (dimension: nJ x nplist x, row-major) - */ - std::vector dJzdp_; - - /** state derivative of event output - * (dimension: nz x `nx_solver`, row-major) - */ - std::vector dzdx_; - - /** parameter derivative of event output - * (dimension: nz x nplist, row-major) - */ - std::vector dzdp_; - - /** state derivative of event regularization variable - * (dimension: `nz` x `nx_solver`, row-major) - */ - std::vector drzdx_; - - /** parameter derivative of event regularization variable - * (dimension: nz x nplist, row-major) - */ - std::vector drzdp_; - - /** parameter derivative of observable - * (dimension: ny x nplist, row-major) - */ - std::vector dydp_; - - /** state derivative of time-resolved observable - * (dimension: `nx_solver` x `ny`, row-major) - */ - std::vector dydx_; - - /** temporary storage of w data across functions (dimension: nw) */ - std::vector w_; - - /** temporary storage for flattened sx, - * (dimension: `nx_solver` x `nplist`, row-major) - */ - std::vector sx_; - - /** temporary storage for `x_rdata` (dimension: `nx_rdata`) */ - std::vector x_rdata_; - - /** temporary storage for `sx_rdata` slice (dimension: `nx_rdata`) */ - std::vector sx_rdata_; - - /** temporary storage for time-resolved observable (dimension: ny) */ - std::vector y_; - - /** data standard deviation for current timepoint (dimension: ny) */ - std::vector sigmay_; - - /** temporary storage for parameter derivative of data standard deviation, - * (dimension: ny x nplist, row-major) - */ - std::vector dsigmaydp_; - - /** temporary storage for event-resolved observable (dimension: nz) */ - std::vector z_; - - /** temporary storage for event regularization (dimension: nz) */ - std::vector rz_; - - /** temporary storage for event standard deviation (dimension: nz) */ - std::vector sigmaz_; - - /** temporary storage for parameter derivative of event standard deviation, - * (dimension: nz x nplist, row-major) - */ - std::vector dsigmazdp_; - - /** temporary storage for change in x after event (dimension: `nx_solver`) */ - std::vector deltax_; - - /** temporary storage for change in sx after event - * (dimension: `nx_solver` x `nplist`, row-major) - */ - std::vector deltasx_; - - /** temporary storage for change in xB after event (dimension: `nx_solver`) */ - std::vector deltaxB_; - - /** temporary storage for change in qB after event - * (dimension: nJ x nplist, row-major) - */ - std::vector deltaqB_; - - /** temporary storage of positified state variables according to - * stateIsNonNegative (dimension: `nx_solver`) */ - AmiVector x_pos_tmp_ {0}; -}; /** * @brief The Model class represents an AMICI ODE/DAE model. diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h new file mode 100644 index 0000000000..f90317ceca --- /dev/null +++ b/include/amici/model_dimensions.h @@ -0,0 +1,177 @@ +#ifndef AMICI_MODEL_DIMENSIONS_H +#define AMICI_MODEL_DIMENSIONS_H + +#include +#include + +namespace amici { + +/** + * @brief Container for model dimensions. + * + * Holds number of states, observables, etc. + */ +struct ModelDimensions { + /** Default ctor */ + ModelDimensions() = default; + + /** + * @brief Constructor with model dimensions + * @param nx_rdata Number of state variables + * @param nxtrue_rdata Number of state variables of the non-augmented model + * @param nx_solver Number of state variables with conservation laws applied + * @param nxtrue_solver Number of state variables of the non-augmented model + * with conservation laws applied + * @param nx_solver_reinit Number of state variables with conservation laws + * subject to reinitialization + * @param np Number of parameters + * @param nk Number of constants + * @param ny Number of observables + * @param nytrue Number of observables of the non-augmented model + * @param nz Number of event observables + * @param nztrue Number of event observables of the non-augmented model + * @param ne Number of events + * @param nJ Number of objective functions + * @param nw Number of repeating elements + * @param ndwdx Number of nonzero elements in the `x` derivative of the + * repeating elements + * @param ndwdp Number of nonzero elements in the `p` derivative of the + * repeating elements + * @param ndwdw Number of nonzero elements in the `w` derivative of the + * repeating elements + * @param ndxdotdw Number of nonzero elements in the \f$w\f$ derivative of + * \f$xdot\f$ + * @param ndJydy Number of nonzero elements in the \f$y\f$ derivative of + * \f$dJy\f$ (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 + */ + ModelDimensions( + const int nx_rdata, const int nxtrue_rdata, const int nx_solver, + const int nxtrue_solver, const int nx_solver_reinit, const int np, + const int nk, 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 ndwdw, const int ndxdotdw, std::vector ndJydy, + const int nnz, const int ubw, const int lbw) + : nx_rdata(nx_rdata), nxtrue_rdata(nxtrue_rdata), nx_solver(nx_solver), + nxtrue_solver(nxtrue_solver), nx_solver_reinit(nx_solver_reinit), + np(np), nk(nk), + ny(ny), nytrue(nytrue), nz(nz), nztrue(nztrue), + ne(ne), nw(nw), ndwdx(ndwdx), ndwdp(ndwdp), ndwdw(ndwdw), + ndxdotdw(ndxdotdw), ndJydy(std::move(ndJydy)), + nnz(nnz), nJ(nJ), ubw(ubw), lbw(lbw) { + Expects(nxtrue_rdata >= 0); + Expects(nxtrue_rdata <= nx_rdata); + Expects(nxtrue_solver >= 0); + Expects(nx_solver <= nx_rdata); + Expects(nxtrue_solver <= nx_solver); + Expects(nx_solver_reinit >= 0); + Expects(nx_solver_reinit <= nx_solver); + Expects(np >= 0); + Expects(nk >= 0); + Expects(nytrue <= ny); + Expects(nytrue >= 0); + Expects(nztrue >= 0); + Expects(nztrue <= nz); + Expects(ne >= 0); + Expects(nw >= 0); + Expects(ndwdx >= 0); + Expects(ndwdx <= nw * nx_solver); + Expects(ndwdp >= 0); + Expects(ndwdp <= nw * np); + Expects(ndwdw >= 0); + Expects(ndwdw <= nw * nw); + Expects(ndxdotdw >= 0); + Expects(nnz >= 0); + Expects(nJ >= 0); + Expects(ubw >= 0); + Expects(lbw >= 0); + } + + /** Number of states */ + int nx_rdata{0}; + + /** Number of states in the unaugmented system */ + int nxtrue_rdata{0}; + + /** Number of states with conservation laws applied */ + int nx_solver{0}; + + /** + * Number of states in the unaugmented system with conservation laws + * applied + */ + int nxtrue_solver{0}; + + /** Number of solver states subject to reinitialization */ + int nx_solver_reinit{0}; + + /** Number of parameters */ + int np{0}; + + /** Number of constants */ + int nk{0}; + + /** Number of observables */ + int ny{0}; + + /** Number of observables in the unaugmented system */ + int nytrue{0}; + + /** Number of event outputs */ + int nz{0}; + + /** Number of event outputs in the unaugmented system */ + int nztrue{0}; + + /** Number of events */ + int ne{0}; + + /** Number of common expressions */ + int nw{0}; + + /** + * Number of nonzero elements in the `x` derivative of the + * repeating elements + */ + int ndwdx {0}; + + /** + * Number of nonzero elements in the `p` derivative of the + * repeating elements + */ + int ndwdp {0}; + + /** + * Number of nonzero elements in the `w` derivative of the + * repeating elements + */ + int ndwdw {0}; + + /** Number of nonzero elements in the \f$w\f$ derivative of \f$xdot\f$ */ + int ndxdotdw {0}; + + /** + * Number of nonzero elements in the \f$y\f$ derivative of + * \f$dJy\f$ (dimension `nytrue`) + */ + std::vector ndJydy; + + /** Number of nonzero entries in Jacobian */ + int nnz{0}; + + /** Dimension of the augmented objective function for 2nd order ASA */ + int nJ{0}; + + /** Upper bandwidth of the Jacobian */ + int ubw{0}; + + /** Lower bandwidth of the Jacobian */ + int lbw{0}; +}; + +} // namespace amici + +#endif // AMICI_MODEL_DIMENSIONS_H diff --git a/include/amici/model_state.h b/include/amici/model_state.h new file mode 100644 index 0000000000..578835283a --- /dev/null +++ b/include/amici/model_state.h @@ -0,0 +1,274 @@ +#ifndef AMICI_MODEL_STATE_H +#define AMICI_MODEL_STATE_H + +#include "amici/defines.h" +#include "amici/sundials_matrix_wrapper.h" +#include "amici/model_dimensions.h" + +#include + +namespace amici { + + +/** + * @brief Exchange format to store and transfer the state of the + * model at a specific timepoint. + * + * This is designed to only encompass the minimal + * number of attributes that need to be transferred. + */ +struct ModelState { + /** + * Flag indicating whether a certain Heaviside function should be active or + * not (dimension: `ne`) + */ + std::vector h; + + /** Total abundances for conservation laws + (dimension: `nx_rdata - nx_solver`) */ + std::vector total_cl; + + /** Sensitivities of total abundances for conservation laws + (dimension: `(nx_rdata-nx_solver) x np`, row-major) */ + std::vector stotal_cl; + + /** Unscaled parameters (dimension: `np`) */ + std::vector unscaledParameters; + + /** Constants (dimension: `nk`) */ + std::vector fixedParameters; + + /** + * Indexes of parameters wrt to which sensitivities are computed + * (dimension: nplist) + */ + std::vector plist; +}; + + +/** + * @brief Storage for `amici::Model` quantities computed based on + * `amici::ModelState` for a specific timepoint. + * + * Serves as workspace for a model simulation to avoid repeated reallocation. + */ +struct ModelStateDerived { + ModelStateDerived() = default; + + /** + * @brief Constructor from model dimensions. + * @param dim Model dimensions + */ + explicit ModelStateDerived(ModelDimensions const& dim); + + /** Sparse Jacobian (dimension: `amici::Model::nnz`) */ + SUNMatrixWrapper J_; + + /** Sparse Backwards Jacobian (dimension: `amici::Model::nnz`) */ + SUNMatrixWrapper JB_; + + /** Sparse dxdotdw temporary storage (dimension: `ndxdotdw`) */ + SUNMatrixWrapper dxdotdw_; + + /** Sparse dwdx temporary storage (dimension: `ndwdx`) */ + SUNMatrixWrapper dwdx_; + + /** Sparse dwdp temporary storage (dimension: `ndwdp`) */ + SUNMatrixWrapper dwdp_; + + /** Dense Mass matrix (dimension: `nx_solver` x `nx_solver`) */ + SUNMatrixWrapper M_; + + /** + * Temporary storage of `dxdotdp_full` data across functions (Python only) + * (dimension: `nplist` x `nx_solver`, nnz: dynamic, + * type `CSC_MAT`) + */ + SUNMatrixWrapper dxdotdp_full; + + /** + * Temporary storage of `dxdotdp_explicit` data across functions (Python only) + * (dimension: `nplist` x `nx_solver`, nnz: `ndxdotdp_explicit`, + * type `CSC_MAT`) + */ + SUNMatrixWrapper dxdotdp_explicit; + + /** + * Temporary storage of `dxdotdp_implicit` data across functions, + * Python-only + * (dimension: `nplist` x `nx_solver`, nnz: dynamic, + * type `CSC_MAT`) + */ + SUNMatrixWrapper dxdotdp_implicit; + + /** + * Temporary storage of `dxdotdx_explicit` data across functions (Python only) + * (dimension: `nplist` x `nx_solver`, nnz: `nxdotdotdx_explicit`, + * type `CSC_MAT`) + */ + SUNMatrixWrapper dxdotdx_explicit; + + /** + * Temporary storage of `dxdotdx_implicit` data across functions, + * Python-only + * (dimension: `nplist` x `nx_solver`, nnz: dynamic, + * type `CSC_MAT`) + */ + SUNMatrixWrapper dxdotdx_implicit; + + /** + * Temporary storage of `dxdotdp` data across functions, Matlab only + * (dimension: `nplist` x `nx_solver` , row-major) + */ + AmiVectorArray dxdotdp {0, 0}; + + /** Sparse observable derivative of data likelihood, only used if + * `pythonGenerated` == `true` (dimension `nytrue`, `nJ` x `ny`, row-major) + */ + std::vector dJydy_; + + /** Observable derivative of data likelihood, only used if + * `pythonGenerated` == `false` (dimension `nJ` x `ny` x `nytrue` , + * row-major) + */ + std::vector dJydy_matlab_; + + /** Observable sigma derivative of data likelihood + * (dimension nJ x ny x nytrue, row-major) + */ + std::vector dJydsigma_; + + /** State derivative of data likelihood + * (dimension `nJ` x `nx_solver`, row-major) + */ + std::vector dJydx_; + + /** Parameter derivative of data likelihood for current timepoint + * (dimension: nJ x nplist, row-major) + */ + std::vector dJydp_; + + /** event output derivative of event likelihood + * (dimension nJ x nz x nztrue, row-major) + */ + std::vector dJzdz_; + + /** event sigma derivative of event likelihood + * (dimension nJ x nz x nztrue, row-major) + */ + std::vector dJzdsigma_; + + /** event output derivative of event likelihood at final timepoint + * (dimension nJ x nz x nztrue, row-major) + */ + std::vector dJrzdz_; + + /** event sigma derivative of event likelihood at final timepoint + * (dimension nJ x nz x nztrue, row-major) + */ + std::vector dJrzdsigma_; + + /** state derivative of event likelihood + * (dimension `nJ` x `nx_solver`, row-major) + */ + std::vector dJzdx_; + + /** parameter derivative of event likelihood for current timepoint + * (dimension: nJ x nplist x, row-major) + */ + std::vector dJzdp_; + + /** state derivative of event output + * (dimension: nz x `nx_solver`, row-major) + */ + std::vector dzdx_; + + /** parameter derivative of event output + * (dimension: nz x nplist, row-major) + */ + std::vector dzdp_; + + /** state derivative of event regularization variable + * (dimension: `nz` x `nx_solver`, row-major) + */ + std::vector drzdx_; + + /** parameter derivative of event regularization variable + * (dimension: nz x nplist, row-major) + */ + std::vector drzdp_; + + /** parameter derivative of observable + * (dimension: ny x nplist, row-major) + */ + std::vector dydp_; + + /** state derivative of time-resolved observable + * (dimension: `nx_solver` x `ny`, row-major) + */ + std::vector dydx_; + + /** temporary storage of w data across functions (dimension: nw) */ + std::vector w_; + + /** temporary storage for flattened sx, + * (dimension: `nx_solver` x `nplist`, row-major) + */ + std::vector sx_; + + /** temporary storage for `x_rdata` (dimension: `nx_rdata`) */ + std::vector x_rdata_; + + /** temporary storage for `sx_rdata` slice (dimension: `nx_rdata`) */ + std::vector sx_rdata_; + + /** temporary storage for time-resolved observable (dimension: ny) */ + std::vector y_; + + /** data standard deviation for current timepoint (dimension: ny) */ + std::vector sigmay_; + + /** temporary storage for parameter derivative of data standard deviation, + * (dimension: ny x nplist, row-major) + */ + std::vector dsigmaydp_; + + /** temporary storage for event-resolved observable (dimension: nz) */ + std::vector z_; + + /** temporary storage for event regularization (dimension: nz) */ + std::vector rz_; + + /** temporary storage for event standard deviation (dimension: nz) */ + std::vector sigmaz_; + + /** temporary storage for parameter derivative of event standard deviation, + * (dimension: nz x nplist, row-major) + */ + std::vector dsigmazdp_; + + /** temporary storage for change in x after event (dimension: `nx_solver`) */ + std::vector deltax_; + + /** temporary storage for change in sx after event + * (dimension: `nx_solver` x `nplist`, row-major) + */ + std::vector deltasx_; + + /** temporary storage for change in xB after event (dimension: `nx_solver`) */ + std::vector deltaxB_; + + /** temporary storage for change in qB after event + * (dimension: nJ x nplist, row-major) + */ + std::vector deltaqB_; + + /** temporary storage of positified state variables according to + * stateIsNonNegative (dimension: `nx_solver`) */ + AmiVector x_pos_tmp_ {0}; +}; + + +} // namespace amici + +#endif // AMICI_MODEL_STATE_H diff --git a/matlab/@amimodel/compileAndLinkModel.m b/matlab/@amimodel/compileAndLinkModel.m index e9675a51fc..8cfd8e7431 100644 --- a/matlab/@amimodel/compileAndLinkModel.m +++ b/matlab/@amimodel/compileAndLinkModel.m @@ -196,7 +196,7 @@ function compileAndLinkModel(modelname, modelSourceFolder, coptim, debug, funs, cppsrc = {'amici', 'symbolic_functions','spline', ... 'edata','rdata', 'exception', ... 'interface_matlab', 'misc', 'simulation_parameters', ... - 'solver', 'solver_cvodes', 'solver_idas', ... + 'solver', 'solver_cvodes', 'solver_idas', 'model_state', ... 'model', 'model_ode', 'model_dae', 'returndata_matlab', ... 'forwardproblem', 'steadystateproblem', 'backwardproblem', 'newton_solver', ... 'abstract_model', 'sundials_matrix_wrapper', 'sundials_linsol_wrapper', ... diff --git a/src/model_state.cpp b/src/model_state.cpp new file mode 100644 index 0000000000..913912ae7f --- /dev/null +++ b/src/model_state.cpp @@ -0,0 +1,16 @@ +#include "amici/model_state.h" + +namespace amici { + +ModelStateDerived::ModelStateDerived(const ModelDimensions &dim) + : J_(dim.nx_solver, dim.nx_solver, dim.nnz, CSC_MAT), + JB_(dim.nx_solver, dim.nx_solver, dim.nnz, CSC_MAT), + dxdotdw_(dim.nx_solver, dim.nw, dim.ndxdotdw, CSC_MAT), + w_(dim.nw), + x_rdata_(dim.nx_rdata, 0.0), + sx_rdata_(dim.nx_rdata, 0.0), + x_pos_tmp_(dim.nx_solver) +{} + + +} // namespace amici diff --git a/swig/amici.i b/swig/amici.i index 13d904c149..fbf308b06d 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -62,6 +62,8 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData) // Include before any other header which uses enums defined there %include "amici/defines.h" +%include "amici/model_dimensions.h" +%include "amici/model_state.h" %include "amici/simulation_parameters.h" %include abstract_model.i