diff --git a/include/amici/model.h b/include/amici/model.h index 0d197b6d6a..533f0139dc 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -112,10 +112,6 @@ class Model : public AbstractModel, public ModelDimensions { * @param o2mode Second order sensitivity mode * @param idlist Indexes indicating algebraic components (DAE only) * @param z2event Mapping of event outputs to events - * @param pythonGenerated Flag indicating matlab or python wrapping - * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit` - * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit` - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -123,9 +119,7 @@ class Model : public AbstractModel, public ModelDimensions { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, amici::SecondOrderMode o2mode, std::vector idlist, - std::vector z2event, bool pythonGenerated = false, - int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0, - int w_recursion_depth = 0, + std::vector z2event, std::map> state_independent_events = {} ); @@ -1458,9 +1452,6 @@ class Model : public AbstractModel, public ModelDimensions { */ std::vector const& getReinitializationStateIdxs() const; - /** Flag indicating Matlab- or Python-based model generation */ - bool pythonGenerated = false; - /** * @brief getter for dxdotdp (matlab generated) * @return dxdotdp @@ -2098,9 +2089,6 @@ class Model : public AbstractModel, public ModelDimensions { realtype min_sigma_{50.0}; private: - /** Recursion */ - int w_recursion_depth_{0}; - /** Simulation parameters, initial state, etc. */ SimulationParameters simulation_parameters_; diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 03ffbf8fb9..a704fcaef2 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -36,10 +36,6 @@ class Model_DAE : public Model { * @param o2mode second order sensitivity mode * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events - * @param pythonGenerated flag indicating matlab or python wrapping - * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit - * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -47,15 +43,12 @@ class Model_DAE : public Model { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, bool const pythonGenerated = false, - int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0, - int const w_recursion_depth = 0, + std::vector const& z2event, std::map> state_independent_events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit, - w_recursion_depth, state_independent_events + state_independent_events ) { derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver); auto M_nnz = static_cast( diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h index b5aa1ba21e..7b6ade81c5 100644 --- a/include/amici/model_dimensions.h +++ b/include/amici/model_dimensions.h @@ -9,7 +9,7 @@ namespace amici { /** * @brief Container for model dimensions. * - * Holds number of states, observables, etc. + * Holds number of state variables, observables, etc. */ struct ModelDimensions { /** Default ctor */ @@ -54,6 +54,11 @@ struct ModelDimensions { * @param nnz Number of nonzero elements in Jacobian * @param ubw Upper matrix bandwidth in the Jacobian * @param lbw Lower matrix bandwidth in the Jacobian + * @param pythonGenerated Flag indicating model creation from Matlab or + * Python + * @param ndxdotdp_explicit Number of nonzero elements in `dxdotdp_explicit` + * @param ndxdotdx_explicit Number of nonzero elements in `dxdotdx_explicit` + * @param w_recursion_depth Recursion depth of fw */ ModelDimensions( int const nx_rdata, int const nxtrue_rdata, int const nx_solver, @@ -64,7 +69,8 @@ struct ModelDimensions { int const ndwdw, int const ndxdotdw, std::vector ndJydy, int const ndxrdatadxsolver, int const ndxrdatadtcl, int const ndtotal_cldx_rdata, int const nnz, int const ubw, - int const lbw + int const lbw, bool pythonGenerated = false, int ndxdotdp_explicit = 0, + int ndxdotdx_explicit = 0, int w_recursion_depth = 0 ) : nx_rdata(nx_rdata) , nxtrue_rdata(nxtrue_rdata) @@ -92,7 +98,11 @@ struct ModelDimensions { , nnz(nnz) , nJ(nJ) , ubw(ubw) - , lbw(lbw) { + , lbw(lbw) + , pythonGenerated(pythonGenerated) + , ndxdotdp_explicit(ndxdotdp_explicit) + , ndxdotdx_explicit(ndxdotdx_explicit) + , w_recursion_depth(w_recursion_depth) { Expects(nxtrue_rdata >= 0); Expects(nxtrue_rdata <= nx_rdata); Expects(nxtrue_solver >= 0); @@ -128,6 +138,9 @@ struct ModelDimensions { Expects(nJ >= 0); Expects(ubw >= 0); Expects(lbw >= 0); + Expects(ndxdotdp_explicit >= 0); + Expects(ndxdotdx_explicit >= 0); + Expects(w_recursion_depth >= 0); } /** Number of states */ @@ -229,6 +242,18 @@ struct ModelDimensions { /** Lower bandwidth of the Jacobian */ int lbw{0}; + + /** Flag indicating model creation from Matlab or Python */ + bool pythonGenerated = false; + + /** Number of nonzero elements in `dxdotdx_explicit` */ + int ndxdotdp_explicit = 0; + + /** Number of nonzero elements in `dxdotdp_explicit` */ + int ndxdotdx_explicit = 0; + + /** Recursion depth of fw */ + int w_recursion_depth = 0; }; } // namespace amici diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 24d410a33c..8386eb929e 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -35,10 +35,6 @@ class Model_ODE : public Model { * @param o2mode second order sensitivity mode * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events - * @param pythonGenerated flag indicating matlab or python wrapping - * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit - * @param ndxdotdx_explicit number of nonzero elements dxdotdx_explicit - * @param w_recursion_depth Recursion depth of fw * @param state_independent_events Map of events with state-independent * triggers functions, mapping trigger timepoints to event indices. */ @@ -46,15 +42,12 @@ class Model_ODE : public Model { ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode const o2mode, std::vector const& idlist, - std::vector const& z2event, bool const pythonGenerated = false, - int const ndxdotdp_explicit = 0, int const ndxdotdx_explicit = 0, - int const w_recursion_depth = 0, + std::vector const& z2event, std::map> state_independent_events = {} ) : Model( model_dimensions, simulation_parameters, o2mode, idlist, z2event, - pythonGenerated, ndxdotdp_explicit, ndxdotdx_explicit, - w_recursion_depth, state_independent_events + state_independent_events ) {} void diff --git a/include/amici/model_state.h b/include/amici/model_state.h index 74d1966d85..6336f8735d 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -332,6 +332,9 @@ struct ModelStateDerived { /** Sparse dwdx implicit temporary storage (shape `ndwdx`) */ std::vector dwdx_hierarchical_; + + /** Temporary storage for dense dJydy (dimension: `nJ` x `ny`) */ + SUNMatrixWrapper dJydy_dense_; }; /** diff --git a/src/model.cpp b/src/model.cpp index de3f4bbb02..0c65cefe9d 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -177,19 +177,15 @@ Model::Model( ModelDimensions const& model_dimensions, SimulationParameters simulation_parameters, SecondOrderMode o2mode, std::vector idlist, std::vector z2event, - bool const pythonGenerated, int const ndxdotdp_explicit, - int const ndxdotdx_explicit, int const w_recursion_depth, std::map> state_independent_events ) : ModelDimensions(model_dimensions) - , pythonGenerated(pythonGenerated) , o2mode(o2mode) , idlist(std::move(idlist)) , state_independent_events_(std::move(state_independent_events)) , derived_state_(model_dimensions) , z2event_(std::move(z2event)) , state_is_non_negative_(nx_solver, false) - , w_recursion_depth_(w_recursion_depth) , simulation_parameters_(std::move(simulation_parameters)) { Expects( model_dimensions.np @@ -217,60 +213,6 @@ Model::Model( root_initial_values_.resize(ne, true); - /* If Matlab wrapped: dxdotdp is a full AmiVector, - if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices - */ - if (pythonGenerated) { - - derived_state_.dwdw_ = SUNMatrixWrapper(nw, nw, ndwdw, CSC_MAT); - // size dynamically adapted for dwdx_ and dwdp_ - derived_state_.dwdx_ = SUNMatrixWrapper(nw, nx_solver, 0, CSC_MAT); - derived_state_.dwdp_ = SUNMatrixWrapper(nw, np(), 0, CSC_MAT); - - for (int irec = 0; irec <= w_recursion_depth_; ++irec) { - /* for the first element we know the exact size, while for all - others we guess the size*/ - derived_state_.dwdp_hierarchical_.emplace_back( - SUNMatrixWrapper(nw, np(), irec * ndwdw + ndwdp, CSC_MAT) - ); - derived_state_.dwdx_hierarchical_.emplace_back( - SUNMatrixWrapper(nw, nx_solver, irec * ndwdw + ndwdx, CSC_MAT) - ); - } - assert( - gsl::narrow(derived_state_.dwdp_hierarchical_.size()) - == w_recursion_depth_ + 1 - ); - assert( - gsl::narrow(derived_state_.dwdx_hierarchical_.size()) - == w_recursion_depth_ + 1 - ); - - derived_state_.dxdotdp_explicit - = SUNMatrixWrapper(nx_solver, np(), ndxdotdp_explicit, CSC_MAT); - // guess size, will be dynamically reallocated - derived_state_.dxdotdp_implicit - = SUNMatrixWrapper(nx_solver, np(), ndwdp + ndxdotdw, CSC_MAT); - derived_state_.dxdotdx_explicit = SUNMatrixWrapper( - nx_solver, nx_solver, ndxdotdx_explicit, CSC_MAT - ); - // guess size, will be dynamically reallocated - derived_state_.dxdotdx_implicit - = SUNMatrixWrapper(nx_solver, nx_solver, ndwdx + ndxdotdw, CSC_MAT); - // dynamically allocate on first call - derived_state_.dxdotdp_full - = SUNMatrixWrapper(nx_solver, np(), 0, CSC_MAT); - - for (int iytrue = 0; iytrue < nytrue; ++iytrue) - derived_state_.dJydy_.emplace_back( - SUNMatrixWrapper(nJ, ny, ndJydy.at(iytrue), CSC_MAT) - ); - } else { - derived_state_.dwdx_ = SUNMatrixWrapper(nw, nx_solver, ndwdx, CSC_MAT); - derived_state_.dwdp_ = SUNMatrixWrapper(nw, np(), ndwdp, CSC_MAT); - derived_state_.dJydy_matlab_ - = std::vector(nJ * nytrue * ny, 0.0); - } requireSensitivitiesForAllParameters(); } @@ -2250,7 +2192,6 @@ void Model::fdJydy(int const it, AmiVector const& x, ExpData const& edata) { if (pythonGenerated) { fdJydsigma(it, x, edata); fdsigmaydy(it, &edata); - SUNMatrixWrapper tmp_dense(nJ, ny); setNaNtoZero(derived_state_.dJydsigma_); setNaNtoZero(derived_state_.dsigmaydy_); @@ -2275,15 +2216,15 @@ void Model::fdJydy(int const it, AmiVector const& x, ExpData const& edata) { // dJydy += dJydsigma * dsigmaydy // C(nJ,ny) A(nJ,ny) * B(ny,ny) // sparse dense dense - tmp_dense.zero(); + derived_state_.dJydy_dense_.zero(); amici_dgemm( BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, nJ, ny, ny, 1.0, &derived_state_.dJydsigma_.at(iyt * nJ * ny), nJ, - derived_state_.dsigmaydy_.data(), ny, 1.0, tmp_dense.data(), nJ + derived_state_.dsigmaydy_.data(), ny, 1.0, derived_state_.dJydy_dense_.data(), nJ ); - auto tmp_sparse = SUNMatrixWrapper(tmp_dense, 0.0, CSC_MAT); + auto tmp_sparse = SUNMatrixWrapper(derived_state_.dJydy_dense_, 0.0, CSC_MAT); auto ret = SUNMatScaleAdd( 1.0, derived_state_.dJydy_.at(iyt), tmp_sparse ); @@ -2871,8 +2812,9 @@ void Model::fsspl(realtype const t) { realtype* sspl_data = derived_state_.sspl_.data(); for (int ip = 0; ip < nplist(); ip++) { for (int ispl = 0; ispl < nspl; ispl++) - sspl_data[ispl + nspl * plist(ip)] - = splines_[ispl].get_sensitivity(t, ip, derived_state_.spl_[ispl]); + sspl_data[ispl + nspl * plist(ip)] = splines_[ispl].get_sensitivity( + t, ip, derived_state_.spl_[ispl] + ); } } @@ -2914,7 +2856,7 @@ void Model::fdwdp(realtype const t, realtype const* x, bool include_static) { derived_state_.sspl_.data(), include_static ); - for (int irecursion = 1; irecursion <= w_recursion_depth_; + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { derived_state_.dwdw_.sparse_multiply( derived_state_.dwdp_hierarchical_.at(irecursion), @@ -2966,7 +2908,7 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { derived_state_.spl_.data(), include_static ); - for (int irecursion = 1; irecursion <= w_recursion_depth_; + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { derived_state_.dwdw_.sparse_multiply( derived_state_.dwdx_hierarchical_.at(irecursion), @@ -2982,7 +2924,8 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { fdwdx( derived_state_.dwdx_.data(), t, x, state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), - derived_state_.w_.data(), state_.total_cl.data(), derived_state_.spl_.data() + derived_state_.w_.data(), state_.total_cl.data(), + derived_state_.spl_.data() ); } diff --git a/src/model_header.template.h b/src/model_header.template.h index 932fdeb1a0..bc87da5a03 100644 --- a/src/model_header.template.h +++ b/src/model_header.template.h @@ -135,7 +135,11 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_NDTOTALCLDXRDATA, // ndtotal_cldx_rdata 0, // nnz TPL_UBW, // ubw - TPL_LBW // lbw + TPL_LBW, // lbw + true, // pythonGenerated + TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit + TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit + TPL_W_RECURSION_DEPTH // w_recursion_depth ), amici::SimulationParameters( std::vector{TPL_FIXED_PARAMETERS}, // fixedParameters @@ -144,10 +148,6 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { TPL_O2MODE, // o2mode std::vector{TPL_ID}, // idlist std::vector{TPL_Z2EVENT}, // z2events - true, // pythonGenerated - TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit - TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit - TPL_W_RECURSION_DEPTH, // w_recursion_depth {TPL_STATE_INDEPENDENT_EVENTS} // state-independent events ) { root_initial_values_ = std::vector( diff --git a/src/model_state.cpp b/src/model_state.cpp index b6d1d9c850..ade1d4c847 100644 --- a/src/model_state.cpp +++ b/src/model_state.cpp @@ -19,6 +19,65 @@ ModelStateDerived::ModelStateDerived(ModelDimensions const& dim) , w_(dim.nw) , x_rdata_(dim.nx_rdata, 0.0) , sx_rdata_(dim.nx_rdata, 0.0) - , x_pos_tmp_(dim.nx_solver) {} + , x_pos_tmp_(dim.nx_solver) { + + // If Matlab wrapped: dxdotdp is a full AmiVector, + // if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices + if (dim.pythonGenerated) { + + dwdw_ = SUNMatrixWrapper(dim.nw, dim.nw, dim.ndwdw, CSC_MAT); + // size dynamically adapted for dwdx_ and dwdp_ + dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, 0, CSC_MAT); + dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, 0, CSC_MAT); + + for (int irec = 0; irec <= dim.w_recursion_depth; ++irec) { + /* for the first element we know the exact size, while for all + others we guess the size*/ + dwdp_hierarchical_.emplace_back(SUNMatrixWrapper( + dim.nw, dim.np, irec * dim.ndwdw + dim.ndwdp, CSC_MAT + )); + dwdx_hierarchical_.emplace_back(SUNMatrixWrapper( + dim.nw, dim.nx_solver, irec * dim.ndwdw + dim.ndwdx, CSC_MAT + )); + } + assert( + gsl::narrow(dwdp_hierarchical_.size()) + == dim.w_recursion_depth + 1 + ); + assert( + gsl::narrow(dwdx_hierarchical_.size()) + == dim.w_recursion_depth + 1 + ); + + dxdotdp_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndxdotdp_explicit, CSC_MAT + ); + // guess size, will be dynamically reallocated + dxdotdp_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndwdp + dim.ndxdotdw, CSC_MAT + ); + dxdotdx_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndxdotdx_explicit, CSC_MAT + ); + // guess size, will be dynamically reallocated + dxdotdx_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndwdx + dim.ndxdotdw, CSC_MAT + ); + // dynamically allocate on first call + dxdotdp_full = SUNMatrixWrapper(dim.nx_solver, dim.np, 0, CSC_MAT); + + for (int iytrue = 0; iytrue < dim.nytrue; ++iytrue) + dJydy_.emplace_back( + SUNMatrixWrapper(dim.nJ, dim.ny, dim.ndJydy.at(iytrue), CSC_MAT) + ); + + dJydy_dense_ = SUNMatrixWrapper(dim.nJ, dim.ny); + } else { + dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, dim.ndwdx, CSC_MAT); + dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, dim.ndwdp, CSC_MAT); + dJydy_matlab_ + = std::vector(dim.nJ * dim.nytrue * dim.ny, 0.0); + } +} } // namespace amici