From 1833e78a67c4a2f0098d28c83fad9e588bbfaf61 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 16:11:46 +0200 Subject: [PATCH 1/5] remove matlab code switches --- include/amici/abstract_model.h | 123 +------- include/amici/model.h | 9 +- include/amici/model_dae.h | 1 - include/amici/model_ode.h | 37 +-- src/abstract_model.cpp | 73 ----- src/backwardproblem.cpp | 17 +- src/model.cpp | 531 +++++++++++---------------------- src/model_dae.cpp | 186 +++++------- src/model_ode.cpp | 181 ++++------- src/newton_solver.cpp | 9 - 10 files changed, 312 insertions(+), 855 deletions(-) diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index b24b9a2e43..3a9072438e 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -287,24 +287,6 @@ class AbstractModel { */ virtual void fdx0(AmiVector& x0, AmiVector& dx0); - /** - * @brief Model-specific implementation of fstau (MATLAB-only, DEPRECATED) - * @param stau total derivative of event timepoint - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param tcl total abundances for conservation laws - * @param sx current state sensitivity - * @param ip sensitivity index - * @param ie event index - */ - virtual void fstau( - realtype* stau, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, realtype const* tcl, - realtype const* sx, int ip, int ie - ); /** * @brief Model-specific implementation of fstau @@ -340,26 +322,9 @@ class AbstractModel { fy(realtype* y, realtype t, realtype const* x, realtype const* p, realtype const* k, realtype const* h, realtype const* w); - /** - * @brief Model-specific implementation of fdydp (MATLAB-only) - * @param dydp partial derivative of observables y w.r.t. model parameters p - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param ip parameter index w.r.t. which the derivative is requested - * @param w repeating elements vector - * @param dwdp Recurring terms in xdot, parameter derivative - */ - virtual void fdydp( - realtype* dydp, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, int ip, realtype const* w, - realtype const* dwdp - ); /** - * @brief Model-specific implementation of fdydp (Python) + * @brief Model-specific implementation of fdydp * @param dydp partial derivative of observables y w.r.t. model parameters p * @param t current time * @param x current state @@ -526,25 +491,7 @@ class AbstractModel { ); /** - * @brief Model-specific implementation of fdeltax (MATLAB-only) - * @param deltax state update - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param ie event index - * @param xdot new model right hand side - * @param xdot_old previous model right hand side - */ - virtual void fdeltax( - realtype* deltax, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, int ie, realtype const* xdot, - realtype const* xdot_old - ); - - /** - * @brief Model-specific implementation of fdeltax (Python-only) + * @brief Model-specific implementation of fdeltax * @param deltax state update * @param t current time * @param x current state @@ -563,31 +510,7 @@ class AbstractModel { ); /** - * @brief Model-specific implementation of fdeltasx (MATLAB-only) - * @param deltasx sensitivity update - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param w repeating elements vector - * @param ip sensitivity index - * @param ie event index - * @param xdot new model right hand side - * @param xdot_old previous model right hand side - * @param sx state sensitivity - * @param stau event-time sensitivity - * @param tcl total abundances for conservation laws - */ - virtual void fdeltasx( - realtype* deltasx, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, realtype const* w, int ip, int ie, - realtype const* xdot, realtype const* xdot_old, realtype const* sx, - realtype const* stau, realtype const* tcl - ); - - /** - * @brief Model-specific implementation of fdeltasx (Python-only) + * @brief Model-specific implementation of fdeltasx * @param deltasx sensitivity update * @param t current time * @param x current state @@ -634,24 +557,6 @@ class AbstractModel { realtype const* xB, realtype const* tcl ); - /** - * @brief Model-specific implementation of fdeltaxB (MATLAB-only!!) - * @param deltaxB adjoint state update - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param ie event index - * @param xdot new model right hand side - * @param xdot_old previous model right hand side - * @param xB current adjoint state - */ - virtual void fdeltaxB( - realtype* deltaxB, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, int ie, realtype const* xdot, - realtype const* xdot_old, realtype const* xB - ); /** * @brief Model-specific implementation of fdeltaqB @@ -676,26 +581,6 @@ class AbstractModel { realtype const* x_old, realtype const* xB ); - /** - * @brief Model-specific implementation of fdeltaqB (MATLAB-only!!) - * @param deltaqB sensitivity update - * @param t current time - * @param x current state - * @param p parameter vector - * @param k constant vector - * @param h Heaviside vector - * @param ip sensitivity index - * @param ie event index - * @param xdot new model right hand side - * @param xdot_old previous model right hand side - * @param xB adjoint state - */ - virtual void fdeltaqB( - realtype* deltaqB, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, int ip, int ie, - realtype const* xdot, realtype const* xdot_old, realtype const* xB - ); - /** * @brief Model-specific implementation of fsigmay * @param sigmay standard deviation of measurements @@ -1002,7 +887,7 @@ class AbstractModel { virtual void fdwdx_rowvals(SUNMatrixWrapper& dwdx); /** - * @brief Model-specific implementation of fdwdw, no w chainrule (Py) + * @brief Model-specific implementation of fdwdw, no w chainrule * @param dwdw partial derivative w wrt w * @param t timepoint * @param x Vector with the states diff --git a/include/amici/model.h b/include/amici/model.h index 5a032d765f..cf86ea3acf 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -72,7 +72,6 @@ enum class ModelQuantity { p, ts, dJydy, - dJydy_matlab, deltaqB, dsigmaydp, dsigmaydy, @@ -1520,13 +1519,7 @@ class Model : public AbstractModel, public ModelDimensions { [[nodiscard]] std::vector const& getReinitializationStateIdxs() const; /** - * @brief getter for dxdotdp (matlab generated) - * @return dxdotdp - */ - [[nodiscard]] AmiVectorArray const& get_dxdotdp() const; - - /** - * @brief getter for dxdotdp (python generated) + * @brief getter for dxdotdp * @return dxdotdp */ [[nodiscard]] SUNMatrixWrapper const& get_dxdotdp_full() const; diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index b871328cd3..91a5653396 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -46,7 +46,6 @@ class Model_DAE : public Model { events ) { SUNContext sunctx = derived_state_.sunctx_; - derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver, sunctx); auto const M_nnz = static_cast( std::reduce(idlist.begin(), idlist.end()) ); diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index f40014ed5e..5e970c2342 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -279,25 +279,9 @@ class Model_ODE : public Model { std::unique_ptr getSolver() override; protected: - /** - * @brief Model specific implementation for fJSparse (Matlab) - * @param JSparse Matrix to which the Jacobian will be written - * @param t timepoint - * @param x Vector with the states - * @param p parameter vector - * @param k constants vector - * @param h Heaviside vector - * @param w vector with helper variables - * @param dwdx derivative of w wrt x - **/ - virtual void fJSparse( - SUNMatrixContent_Sparse JSparse, realtype t, realtype const* x, - realtype const* p, realtype const* k, realtype const* h, - realtype const* w, realtype const* dwdx - ); /** - * @brief Model specific implementation for fJSparse, data only (Py) + * @brief Model specific implementation for fJSparse, data only * @param JSparse Matrix to which the Jacobian will be written * @param t timepoint * @param x Vector with the states @@ -355,28 +339,9 @@ class Model_ODE : public Model { realtype const* k, realtype const* h, realtype const* w ) = 0; - /** - * @brief Model specific implementation of fdxdotdp, with w chainrule - * (Matlab) - * @param dxdotdp partial derivative xdot wrt p - * @param t timepoint - * @param x Vector with the states - * @param p parameter vector - * @param k constants vector - * @param h Heaviside vector - * @param ip parameter index - * @param w vector with helper variables - * @param dwdp derivative of w wrt p - */ - virtual void fdxdotdp( - realtype* dxdotdp, realtype t, realtype const* x, realtype const* p, - realtype const* k, realtype const* h, int ip, realtype const* w, - realtype const* dwdp - ); /** * @brief Model specific implementation of fdxdotdp_explicit, no w chainrule - * (Py) * @param dxdotdp_explicit partial derivative xdot wrt p * @param t timepoint * @param x Vector with the states diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp index f7db90f2c8..f3b5e5f09a 100644 --- a/src/abstract_model.cpp +++ b/src/abstract_model.cpp @@ -54,18 +54,6 @@ void AbstractModel::fdx0(AmiVector& /*x0*/, AmiVector& /*dx0*/) { // no-op default implementation } -void AbstractModel::fstau( - realtype* /*stau*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - realtype const* /*tcl*/, realtype const* /*sx*/, int const /*ip*/, - int const /*ie*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} void AbstractModel::fstau( realtype* /*stau*/, realtype const /*t*/, realtype const* /*x*/, @@ -91,17 +79,6 @@ void AbstractModel:: ); } -void AbstractModel::fdydp( - realtype* /*dydp*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - int const /*ip*/, realtype const* /*w*/, realtype const* /*dwdp*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} void AbstractModel::fdydp( realtype* /*dydp*/, realtype const /*t*/, realtype const* /*x*/, @@ -222,17 +199,6 @@ void AbstractModel::fdrzdx( ); } -void AbstractModel::fdeltax( - realtype* /*deltax*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - int const /*ie*/, realtype const* /*xdot*/, realtype const* /*xdot_old*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} void AbstractModel::fdeltax( realtype* /*deltax*/, realtype const /*t*/, realtype const* /*x*/, @@ -247,20 +213,6 @@ void AbstractModel::fdeltax( ); } -void AbstractModel::fdeltasx( - realtype* /*deltasx*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - realtype const* /*w*/, int const /*ip*/, int const /*ie*/, - realtype const* /*xdot*/, realtype const* /*xdot_old*/, - realtype const* /*sx*/, realtype const* /*stau*/, realtype const* /*tcl*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} - void AbstractModel::fdeltasx( realtype* /*deltasx*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, @@ -290,18 +242,6 @@ void AbstractModel::fdeltaxB( ); } -void AbstractModel::fdeltaxB( - realtype* /*deltaxB*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - int const /*ie*/, realtype const* /*xdot*/, realtype const* /*xdot_old*/, - realtype const* /*xB*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} void AbstractModel::fdeltaqB( realtype* /*deltaqB*/, realtype const /*t*/, realtype const* /*x*/, @@ -317,19 +257,6 @@ void AbstractModel::fdeltaqB( ); } -void AbstractModel::fdeltaqB( - realtype* /*deltaqB*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - int const /*ip*/, int const /*ie*/, realtype const* /*xdot*/, - realtype const* /*xdot_old*/, realtype const* /*xB*/ -) { - throw AmiException( - "Requested functionality is not supported as %s is " - "not implemented for this model!", - __func__ - ); -} - void AbstractModel::fsigmay( realtype* /*sigmay*/, realtype const /*t*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*y*/ diff --git a/src/backwardproblem.cpp b/src/backwardproblem.cpp index 1c6309ef55..dbbcc17eb8 100644 --- a/src/backwardproblem.cpp +++ b/src/backwardproblem.cpp @@ -342,17 +342,12 @@ void computeQBfromQ( // set to zero first, as multiplication adds to existing value yQB.zero(); // yQB += dxdotdp * yQ - if (model.pythonGenerated) { - // fill dxdotdp with current values - auto const& plist = model.getParameterList(); - model.fdxdotdp(state.t, state.x, state.dx); - model.get_dxdotdp_full().multiply( - yQB.getNVector(), yQ.getNVector(), plist, true - ); - } else { - for (int ip = 0; ip < model.nplist(); ++ip) - yQB[ip] = dotProd(yQ, model.get_dxdotdp()[ip]); - } + // fill dxdotdp with current values + auto const& plist = model.getParameterList(); + model.fdxdotdp(state.t, state.x, state.dx); + model.get_dxdotdp_full().multiply( + yQB.getNVector(), yQ.getNVector(), plist, true + ); } SteadyStateBackwardProblem::SteadyStateBackwardProblem( diff --git a/src/model.cpp b/src/model.cpp index fc29c52ec5..3e1cc4e107 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -57,7 +57,6 @@ std::map const model_quantity_to_str{ {ModelQuantity::p, "p"}, {ModelQuantity::ts, "ts"}, {ModelQuantity::dJydy, "dJydy"}, - {ModelQuantity::dJydy_matlab, "dJydy"}, {ModelQuantity::deltaqB, "deltaqB"}, {ModelQuantity::dsigmaydp, "dsigmaydp"}, {ModelQuantity::dsigmaydy, "dsigmaydy"}, @@ -197,23 +196,7 @@ Model::Model( == gsl::narrow(simulation_parameters_.fixedParameters.size()) ); - Expects( - (events_.empty() && !pythonGenerated) - || (events_.size() == (unsigned long)ne) - ); - - if (events_.empty()) { - // for MATLAB generated models, create event objects here - for (int ie = 0; ie < ne; ie++) { - events_.emplace_back( - // The default for use_values_from_trigger_time used to be - // `true` in SBML. However, the MATLAB interface always - // implicitly used `false`, so we keep it that way until it - // will be removed completely. - std::string("event_") + std::to_string(ie), false, true, 0 - ); - } - } + Expects((events_.size() == (unsigned long)ne)); simulation_parameters_.pscale = std::vector( model_dimensions.np, ParameterScaling::none @@ -1149,26 +1132,23 @@ void Model::getObservableSigmaSensitivity( fdsigmaydp(it, edata); writeSlice(derived_state_.dsigmaydp_, ssigmay); - if (pythonGenerated) { - // ssigmay = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp - // = dsigmaydy*sy+dsigmaydp + // ssigmay = dsigmaydy*(dydx_solver*sx+dydp)+dsigmaydp + // = dsigmaydy*sy+dsigmaydp - fdsigmaydy(it, edata); + fdsigmaydy(it, edata); - // compute ssigmay = 1.0 * dsigmaydp + 1.0 * dsigmaydy * sy - // dsigmaydp C[ny,nplist] += dsigmaydy A[ny,ny] * sy B[ny,nplist] - // M N M K K N - // ldc lda ldb - setNaNtoZero(derived_state_.dsigmaydy_); - derived_state_.sy_.assign(sy.begin(), sy.end()); - setNaNtoZero(derived_state_.sy_); - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, - BLASTranspose::noTrans, ny, nplist(), ny, 1.0, - derived_state_.dsigmaydy_.data(), ny, derived_state_.sy_.data(), ny, - 1.0, ssigmay.data(), ny - ); - } + // compute ssigmay = 1.0 * dsigmaydp + 1.0 * dsigmaydy * sy + // dsigmaydp C[ny,nplist] += dsigmaydy A[ny,ny] * sy B[ny,nplist] + // M N M K K N + // ldc lda ldb + setNaNtoZero(derived_state_.dsigmaydy_); + derived_state_.sy_.assign(sy.begin(), sy.end()); + setNaNtoZero(derived_state_.sy_); + amici_dgemm( + BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, + ny, nplist(), ny, 1.0, derived_state_.dsigmaydy_.data(), ny, + derived_state_.sy_.data(), ny, 1.0, ssigmay.data(), ny + ); if (always_check_finite_) checkFinite(ssigmay, ModelQuantity::ssigmay, nplist()); @@ -1252,40 +1232,32 @@ void Model::getEventSensitivity( gsl::span sz, int const ie, realtype const t, AmiVector const& x, AmiVectorArray const& sx ) { - if (pythonGenerated) { - if (!nz) - return; + if (!nz) + return; - fdzdx(ie, t, x); - fdzdp(ie, t, x); + fdzdx(ie, t, x); + fdzdp(ie, t, x); - derived_state_.sx_.resize(nx_solver * nplist()); - sx.flatten_to_vector(derived_state_.sx_); + derived_state_.sx_.resize(nx_solver * nplist()); + sx.flatten_to_vector(derived_state_.sx_); - // compute sz = 1.0*dzdx*sx + 1.0*dzdp - // dzdx A[nz,nx_solver] * sx B[nx_solver,nplist] = sz C[nz,nplist] - // M K K N M N - // lda ldb ldc - setNaNtoZero(derived_state_.dzdx_); - setNaNtoZero(derived_state_.sx_); - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, - BLASTranspose::noTrans, nz, nplist(), nx_solver, 1.0, - derived_state_.dzdx_.data(), nz, derived_state_.sx_.data(), - nx_solver, 1.0, derived_state_.dzdp_.data(), nz - ); + // compute sz = 1.0*dzdx*sx + 1.0*dzdp + // dzdx A[nz,nx_solver] * sx B[nx_solver,nplist] = sz C[nz,nplist] + // M K K N M N + // lda ldb ldc + setNaNtoZero(derived_state_.dzdx_); + setNaNtoZero(derived_state_.sx_); + amici_dgemm( + BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, + nz, nplist(), nx_solver, 1.0, derived_state_.dzdx_.data(), nz, + derived_state_.sx_.data(), nx_solver, 1.0, derived_state_.dzdp_.data(), + nz + ); - addSlice(derived_state_.dzdp_, sz); + addSlice(derived_state_.dzdp_, sz); - if (always_check_finite_) - checkFinite(sz, ModelQuantity::sz, nplist()); - } else { - for (int ip = 0; ip < nplist(); ip++) { - fsz(&sz[ip * nz], ie, t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), sx.data(ip), plist(ip)); - } - } + if (always_check_finite_) + checkFinite(sz, ModelQuantity::sz, nplist()); } void Model::getUnobservedEventSensitivity( @@ -1310,42 +1282,33 @@ void Model::getEventRegularizationSensitivity( gsl::span srz, int const ie, realtype const t, AmiVector const& x, AmiVectorArray const& sx ) { - if (pythonGenerated) { - if (!nz) - return; - fdrzdx(ie, t, x); - fdrzdp(ie, t, x); + if (!nz) + return; - derived_state_.sx_.resize(nx_solver * nplist()); - sx.flatten_to_vector(derived_state_.sx_); + fdrzdx(ie, t, x); + fdrzdp(ie, t, x); - // compute srz = 1.0*drzdx*sx + 1.0*drzdp - // drzdx A[nz,nx_solver] * sx B[nx_solver,nplist] = srz C[nz,nplist] - // M K K N M N - // lda ldb ldc - setNaNtoZero(derived_state_.drzdx_); - setNaNtoZero(derived_state_.sx_); - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, - BLASTranspose::noTrans, nz, nplist(), nx_solver, 1.0, - derived_state_.drzdx_.data(), nz, derived_state_.sx_.data(), - nx_solver, 1.0, derived_state_.drzdp_.data(), nz - ); + derived_state_.sx_.resize(nx_solver * nplist()); + sx.flatten_to_vector(derived_state_.sx_); + + // compute srz = 1.0*drzdx*sx + 1.0*drzdp + // drzdx A[nz,nx_solver] * sx B[nx_solver,nplist] = srz C[nz,nplist] + // M K K N M N + // lda ldb ldc + setNaNtoZero(derived_state_.drzdx_); + setNaNtoZero(derived_state_.sx_); + amici_dgemm( + BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, + nz, nplist(), nx_solver, 1.0, derived_state_.drzdx_.data(), nz, + derived_state_.sx_.data(), nx_solver, 1.0, derived_state_.drzdp_.data(), + nz + ); - addSlice(derived_state_.drzdp_, srz); + addSlice(derived_state_.drzdp_, srz); - if (always_check_finite_) - checkFinite(srz, ModelQuantity::srz, nplist()); - } else { - for (int ip = 0; ip < nplist(); ip++) { - fsrz( - &srz[ip * nz], ie, t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), sx.data(ip), plist(ip) - ); - } - } + if (always_check_finite_) + checkFinite(srz, ModelQuantity::srz, nplist()); } void Model::getEventSigma( @@ -1468,24 +1431,12 @@ void Model::getEventTimeSensitivity( std::ranges::fill(stau, 0.0); - if (pythonGenerated) { - for (int ip = 0; ip < nplist(); ip++) { - fstau( - &stau.at(ip), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), dx.data(), state_.total_cl.data(), sx.data(ip), - plist(ip), ie - ); - } - } else { - for (int ip = 0; ip < nplist(); ip++) { - fstau( - &stau.at(ip), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), state_.total_cl.data(), sx.data(ip), plist(ip), - ie - ); - } + for (int ip = 0; ip < nplist(); ip++) { + fstau( + &stau.at(ip), t, computeX_pos(x), state_.unscaledParameters.data(), + state_.fixedParameters.data(), state_.h.data(), dx.data(), + state_.total_cl.data(), sx.data(ip), plist(ip), ie + ); } } @@ -1499,22 +1450,11 @@ void Model::addStateEventUpdate( std::copy_n(computeX_pos(x), nx_solver, x.data()); // compute update - if (pythonGenerated) { - fdeltax( - derived_state_.deltax_.data(), t, x.data(), - state.unscaledParameters.data(), state.fixedParameters.data(), - state.h.data(), ie, xdot.data(), xdot_old.data(), x_old.data() - ); - - } else { - // For MATLAB-imported models, use_values_from_trigger_time=true - // is not supported, and thus, x_old == x, always. - fdeltax( - derived_state_.deltax_.data(), t, x_old.data(), - state.unscaledParameters.data(), state.fixedParameters.data(), - state.h.data(), ie, xdot.data(), xdot_old.data() - ); - } + fdeltax( + derived_state_.deltax_.data(), t, x.data(), + state.unscaledParameters.data(), state.fixedParameters.data(), + state.h.data(), ie, xdot.data(), xdot_old.data(), x_old.data() + ); if (always_check_finite_) { checkFinite(derived_state_.deltax_, ModelQuantity::deltax, t); @@ -1536,25 +1476,13 @@ void Model::addStateSensitivityEventUpdate( derived_state_.deltasx_.assign(nx_solver, 0.0); // compute update - if (pythonGenerated) { - fdeltasx( - derived_state_.deltasx_.data(), t, x.data(), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), derived_state_.w_.data(), plist(ip), ie, - xdot.data(), xdot_old.data(), sx_old.data(ip), &stau.at(ip), - state_.total_cl.data(), x_old.data() - ); - } else { - // For MATLAB-imported models, use_values_from_trigger_time=true - // is not supported, and thus, x_old == x, always. - fdeltasx( - derived_state_.deltasx_.data(), t, x_old.data(), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), derived_state_.w_.data(), plist(ip), ie, - xdot.data(), xdot_old.data(), sx.data(ip), &stau.at(ip), - state_.total_cl.data() - ); - } + fdeltasx( + derived_state_.deltasx_.data(), t, x.data(), + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), derived_state_.w_.data(), plist(ip), ie, + xdot.data(), xdot_old.data(), sx_old.data(ip), &stau.at(ip), + state_.total_cl.data(), x_old.data() + ); if (always_check_finite_) { checkFinite( derived_state_.deltasx_, ModelQuantity::deltasx, nplist() @@ -1576,20 +1504,12 @@ void Model::addAdjointStateEventUpdate( derived_state_.deltaxB_.assign(nxtrue_solver * nJ, 0.0); // compute update - if (pythonGenerated) { - fdeltaxB( - derived_state_.deltaxB_.data(), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), dx.data(), ie, xdot.data(), xdot_old.data(), - x_old.data(), xB.data(), state_.total_cl.data() - ); - } else { - fdeltaxB( - derived_state_.deltaxB_.data(), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), ie, xdot.data(), xdot_old.data(), xB.data() - ); - } + fdeltaxB( + derived_state_.deltaxB_.data(), t, computeX_pos(x), + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), dx.data(), ie, xdot.data(), xdot_old.data(), + x_old.data(), xB.data(), state_.total_cl.data() + ); if (always_check_finite_) { checkFinite(derived_state_.deltaxB_, ModelQuantity::deltaxB, t); @@ -1610,21 +1530,12 @@ void Model::addAdjointQuadratureEventUpdate( for (int ip = 0; ip < nplist(); ip++) { derived_state_.deltaqB_.assign(nJ, 0.0); - if (pythonGenerated) { - fdeltaqB( - derived_state_.deltaqB_.data(), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), dx.data(), plist(ip), ie, xdot.data(), - xdot_old.data(), x_old.data(), xB.data() - ); - } else { - fdeltaqB( - derived_state_.deltaqB_.data(), t, computeX_pos(x), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), plist(ip), ie, xdot.data(), xdot_old.data(), - xB.data() - ); - } + fdeltaqB( + derived_state_.deltaqB_.data(), t, computeX_pos(x), + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), dx.data(), plist(ip), ie, xdot.data(), + xdot_old.data(), x_old.data(), xB.data() + ); for (int iJ = 0; iJ < nJ; ++iJ) xQB.at(ip * nJ + iJ) += derived_state_.deltaqB_.at(iJ); @@ -1789,7 +1700,6 @@ int Model::checkFinite( col_id += " " + getParameterIds()[plist(gsl::narrow(col))]; break; case ModelQuantity::dJydy: - case ModelQuantity::dJydy_matlab: case ModelQuantity::dJydsigma: if (hasObservableIds()) col_id += " " + getObservableIds()[col]; @@ -2117,12 +2027,7 @@ void Model::checkLLHBufferSize( ); } -void Model::initializeVectors() { - sx0data_.clear(); - if (!pythonGenerated) - derived_state_.dxdotdp - = AmiVectorArray(nx_solver, nplist(), derived_state_.sunctx_); -} +void Model::initializeVectors() { sx0data_.clear(); } void Model::fy(realtype const t, AmiVector const& x) { if (!ny) @@ -2156,22 +2061,13 @@ void Model::fdydp(realtype const t, AmiVector const& x) { /* get dydp slice (ny) for current time and parameter */ for (int ip = 0; ip < nplist(); ip++) - if (pythonGenerated) { - fdydp( - &derived_state_.dydp_.at(ip * ny), t, x_pos, - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), plist(ip), derived_state_.w_.data(), - state_.total_cl.data(), state_.stotal_cl.data(), - derived_state_.spl_.data(), derived_state_.sspl_.data() - ); - } else { - fdydp( - &derived_state_.dydp_.at(ip * ny), t, x_pos, - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), plist(ip), derived_state_.w_.data(), - derived_state_.dwdp_.data() - ); - } + fdydp( + &derived_state_.dydp_.at(ip * ny), t, x_pos, + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), plist(ip), derived_state_.w_.data(), + state_.total_cl.data(), state_.stotal_cl.data(), + derived_state_.spl_.data(), derived_state_.sspl_.data() + ); if (always_check_finite_) { checkFinite(derived_state_.dydp_, ModelQuantity::dydp, nplist()); @@ -2310,81 +2206,57 @@ void Model::fdJydy(int const it, AmiVector const& x, ExpData const& edata) { fy(edata.getTimepoint(it), x); fsigmay(it, &edata); - if (pythonGenerated) { - fdJydsigma(it, x, edata); - fdsigmaydy(it, &edata); + fdJydsigma(it, x, edata); + fdsigmaydy(it, &edata); - setNaNtoZero(derived_state_.dJydsigma_); - setNaNtoZero(derived_state_.dsigmaydy_); - for (int iyt = 0; iyt < nytrue; iyt++) { - if (!derived_state_.dJydy_.at(iyt).capacity()) - continue; - derived_state_.dJydy_.at(iyt).zero(); - fdJydy_colptrs(derived_state_.dJydy_.at(iyt), iyt); - fdJydy_rowvals(derived_state_.dJydy_.at(iyt), iyt); + setNaNtoZero(derived_state_.dJydsigma_); + setNaNtoZero(derived_state_.dsigmaydy_); + for (int iyt = 0; iyt < nytrue; iyt++) { + if (!derived_state_.dJydy_.at(iyt).capacity()) + continue; + derived_state_.dJydy_.at(iyt).zero(); + fdJydy_colptrs(derived_state_.dJydy_.at(iyt), iyt); + fdJydy_rowvals(derived_state_.dJydy_.at(iyt), iyt); - if (!edata.isSetObservedData(it, iyt)) - continue; + if (!edata.isSetObservedData(it, iyt)) + continue; - // get dJydy slice (ny) for current timepoint and observable - fdJydy( - derived_state_.dJydy_.at(iyt).data(), iyt, - state_.unscaledParameters.data(), state_.fixedParameters.data(), - derived_state_.y_.data(), derived_state_.sigmay_.data(), - edata.getObservedDataPtr(it) - ); + // get dJydy slice (ny) for current timepoint and observable + fdJydy( + derived_state_.dJydy_.at(iyt).data(), iyt, + state_.unscaledParameters.data(), state_.fixedParameters.data(), + derived_state_.y_.data(), derived_state_.sigmay_.data(), + edata.getObservedDataPtr(it) + ); - // dJydy += dJydsigma * dsigmaydy - // C(nJ,ny) A(nJ,ny) * B(ny,ny) - // sparse dense dense - 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, - derived_state_.dJydy_dense_.data(), nJ - ); + // dJydy += dJydsigma * dsigmaydy + // C(nJ,ny) A(nJ,ny) * B(ny,ny) + // sparse dense dense + 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, + derived_state_.dJydy_dense_.data(), nJ + ); - auto tmp_sparse - = SUNMatrixWrapper(derived_state_.dJydy_dense_, 0.0, CSC_MAT); - auto ret = SUNMatScaleAdd( - 1.0, derived_state_.dJydy_.at(iyt), tmp_sparse + auto tmp_sparse + = SUNMatrixWrapper(derived_state_.dJydy_dense_, 0.0, CSC_MAT); + auto ret + = SUNMatScaleAdd(1.0, derived_state_.dJydy_.at(iyt), tmp_sparse); + if (ret != SUN_SUCCESS) { + throw AmiException( + "SUNMatScaleAdd failed with status %d in %s", ret, __func__ ); - if (ret != SUN_SUCCESS) { - throw AmiException( - "SUNMatScaleAdd failed with status %d in %s", ret, __func__ - ); - } - derived_state_.dJydy_.at(iyt).refresh(); - - if (always_check_finite_) { - checkFinite( - gsl::make_span(derived_state_.dJydy_.at(iyt).get()), - ModelQuantity::dJydy, ny - ); - } } - } else { - std::ranges::fill(derived_state_.dJydy_matlab_, 0.0); - for (int iyt = 0; iyt < nytrue; iyt++) { - if (!edata.isSetObservedData(it, iyt)) - continue; - fdJydy( - &derived_state_.dJydy_matlab_.at(iyt * ny * nJ), iyt, - state_.unscaledParameters.data(), state_.fixedParameters.data(), - derived_state_.y_.data(), derived_state_.sigmay_.data(), - edata.getObservedDataPtr(it) + derived_state_.dJydy_.at(iyt).refresh(); + + if (always_check_finite_) { + checkFinite( + gsl::make_span(derived_state_.dJydy_.at(iyt).get()), + ModelQuantity::dJydy, ny ); - if (always_check_finite_) { - // get dJydy slice (ny) for current timepoint and observable - checkFinite( - gsl::span( - &derived_state_.dJydy_matlab_[iyt * ny * nJ], ny * nJ - ), - ModelQuantity::dJydy, ny - ); - } } } } @@ -2440,27 +2312,16 @@ void Model::fdJydp(int const it, AmiVector const& x, ExpData const& edata) { if (!edata.isSetObservedData(it, iyt)) continue; - if (pythonGenerated) { - // dJydp = 1.0 * dJydp + 1.0 * dJydy * dydp - for (int iplist = 0; iplist < nplist(); ++iplist) { - derived_state_.dJydy_.at(iyt).multiply( - gsl::span( - &derived_state_.dJydp_.at(iplist * nJ), nJ - ), - gsl::span( - &derived_state_.dydp_.at(iplist * ny), ny - ) - ); - } - } else { - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, - BLASTranspose::noTrans, nJ, nplist(), ny, 1.0, - &derived_state_.dJydy_matlab_.at(iyt * nJ * ny), nJ, - derived_state_.dydp_.data(), ny, 1.0, - derived_state_.dJydp_.data(), nJ + // dJydp = 1.0 * dJydp + 1.0 * dJydy * dydp + for (int iplist = 0; iplist < nplist(); ++iplist) { + derived_state_.dJydy_.at(iyt).multiply( + gsl::span(&derived_state_.dJydp_.at(iplist * nJ), nJ), + gsl::span( + &derived_state_.dydp_.at(iplist * ny), ny + ) ); } + // dJydp = 1.0 * dJydp + 1.0 * dJydsigma * dsigmaydp amici_dgemm( BLASLayout::colMajor, BLASTranspose::noTrans, @@ -2492,22 +2353,10 @@ void Model::fdJydx(int const it, AmiVector const& x, ExpData const& edata) { // M K K N M N // lda ldb ldc - if (pythonGenerated) { - for (int ix = 0; ix < nx_solver; ++ix) { - derived_state_.dJydy_.at(iyt).multiply( - gsl::span(&derived_state_.dJydx_.at(ix * nJ), nJ), - gsl::span( - &derived_state_.dydx_.at(ix * ny), ny - ) - ); - } - } else { - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, - BLASTranspose::noTrans, nJ, nx_solver, ny, 1.0, - &derived_state_.dJydy_matlab_.at(iyt * ny * nJ), nJ, - derived_state_.dydx_.data(), ny, 1.0, - derived_state_.dJydx_.data(), nJ + for (int ix = 0; ix < nx_solver; ++ix) { + derived_state_.dJydy_.at(iyt).multiply( + gsl::span(&derived_state_.dJydx_.at(ix * nJ), nJ), + gsl::span(&derived_state_.dydx_.at(ix * ny), ny) ); } } @@ -2958,45 +2807,30 @@ void Model::fdwdp(realtype const t, realtype const* x, bool include_static) { fw(t, x, include_static); derived_state_.dwdp_.zero(); - if (pythonGenerated) { - if (!ndwdp) - return; - fsspl(t); - fdwdw(t, x, include_static); - if (include_static) { - derived_state_.dwdp_hierarchical_.at(0).zero(); - fdwdp_colptrs(derived_state_.dwdp_hierarchical_.at(0)); - fdwdp_rowvals(derived_state_.dwdp_hierarchical_.at(0)); - } - fdwdp( - derived_state_.dwdp_hierarchical_.at(0).data(), t, x, - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), - state_.stotal_cl.data(), derived_state_.spl_.data(), - derived_state_.sspl_.data(), include_static - ); - - for (int irecursion = 1; irecursion <= w_recursion_depth; - irecursion++) { - derived_state_.dwdw_.sparse_multiply( - derived_state_.dwdp_hierarchical_.at(irecursion), - derived_state_.dwdp_hierarchical_.at(irecursion - 1) - ); - } - derived_state_.dwdp_.sparse_sum(derived_state_.dwdp_hierarchical_); + if (!ndwdp) + return; + fsspl(t); + fdwdw(t, x, include_static); + if (include_static) { + derived_state_.dwdp_hierarchical_.at(0).zero(); + fdwdp_colptrs(derived_state_.dwdp_hierarchical_.at(0)); + fdwdp_rowvals(derived_state_.dwdp_hierarchical_.at(0)); + } + fdwdp( + derived_state_.dwdp_hierarchical_.at(0).data(), t, x, + state_.unscaledParameters.data(), state_.fixedParameters.data(), + state_.h.data(), derived_state_.w_.data(), state_.total_cl.data(), + state_.stotal_cl.data(), derived_state_.spl_.data(), + derived_state_.sspl_.data(), include_static + ); - } else { - if (!derived_state_.dwdp_.capacity()) - return; - // matlab generated - fdwdp( - derived_state_.dwdp_.data(), t, x, state_.unscaledParameters.data(), - state_.fixedParameters.data(), state_.h.data(), - derived_state_.w_.data(), state_.total_cl.data(), - state_.stotal_cl.data(), derived_state_.spl_.data(), - derived_state_.sspl_.data() + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { + derived_state_.dwdw_.sparse_multiply( + derived_state_.dwdp_hierarchical_.at(irecursion), + derived_state_.dwdp_hierarchical_.at(irecursion - 1) ); } + derived_state_.dwdp_.sparse_sum(derived_state_.dwdp_hierarchical_); if (always_check_finite_) { checkFinite(derived_state_.dwdp_, ModelQuantity::dwdp, t); @@ -3012,7 +2846,6 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { fw(t, x, include_static); derived_state_.dwdx_.zero(); - if (pythonGenerated) { fdwdw(t, x, include_static); auto&& dwdx_hierarchical_0 = derived_state_.dwdx_hierarchical_.at(0); @@ -3040,18 +2873,6 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { } derived_state_.dwdx_.sparse_sum(derived_state_.dwdx_hierarchical_); - } else { - if (!derived_state_.dwdx_.capacity()) - return; - derived_state_.dwdx_.zero(); - 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() - ); - } - if (always_check_finite_) { checkFinite(derived_state_.dwdx_, ModelQuantity::dwdx, t); } @@ -3250,13 +3071,7 @@ std::vector const& Model::getReinitializationStateIdxs() const { return simulation_parameters_.reinitialization_state_idxs_sim; } -AmiVectorArray const& Model::get_dxdotdp() const { - assert(!pythonGenerated); - return derived_state_.dxdotdp; -} - SUNMatrixWrapper const& Model::get_dxdotdp_full() const { - assert(pythonGenerated); return derived_state_.dxdotdp_full; } diff --git a/src/model_dae.cpp b/src/model_dae.cpp index 70a4e8fcb7..cb0fcd1277 100644 --- a/src/model_dae.cpp +++ b/src/model_dae.cpp @@ -33,46 +33,34 @@ void Model_DAE::fJSparse( ) { auto const x_pos = computeX_pos(x); fdwdx(t, N_VGetArrayPointerConst(x_pos), false); - if (pythonGenerated) { - auto JSparse = SUNMatrixWrapper(J); - // python generated - derived_state_.dxdotdx_explicit.zero(); - derived_state_.dxdotdx_implicit.zero(); - if (derived_state_.dxdotdx_explicit.capacity()) { - fdxdotdx_explicit_colptrs(derived_state_.dxdotdx_explicit); - fdxdotdx_explicit_rowvals(derived_state_.dxdotdx_explicit); - fdxdotdx_explicit( - derived_state_.dxdotdx_explicit.data(), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), N_VGetArrayPointerConst(dx), - derived_state_.w_.data() - ); - } - fdxdotdw(t, x_pos, dx); - /* Sparse matrix multiplication - dxdotdx_implicit += dxdotdw * dwdx */ - derived_state_.dxdotdw_.sparse_multiply( - derived_state_.dxdotdx_implicit, derived_state_.dwdx_ - ); - derived_state_.dfdx_.sparse_add( - derived_state_.dxdotdx_explicit, 1.0, - derived_state_.dxdotdx_implicit, 1.0 - ); - fM(t, x_pos); - JSparse.sparse_add( - derived_state_.MSparse_, -cj, derived_state_.dfdx_, 1.0 - ); - } else { - fJSparse( - static_cast(SM_CONTENT_S(J)), t, + auto JSparse = SUNMatrixWrapper(J); + // python generated + derived_state_.dxdotdx_explicit.zero(); + derived_state_.dxdotdx_implicit.zero(); + if (derived_state_.dxdotdx_explicit.capacity()) { + fdxdotdx_explicit_colptrs(derived_state_.dxdotdx_explicit); + fdxdotdx_explicit_rowvals(derived_state_.dxdotdx_explicit); + fdxdotdx_explicit( + derived_state_.dxdotdx_explicit.data(), t, N_VGetArrayPointerConst(x_pos), state_.unscaledParameters.data(), - state_.fixedParameters.data(), state_.h.data(), cj, - N_VGetArrayPointerConst(dx), derived_state_.w_.data(), - derived_state_.dwdx_.data() + state_.fixedParameters.data(), state_.h.data(), + N_VGetArrayPointerConst(dx), derived_state_.w_.data() ); } + fdxdotdw(t, x_pos, dx); + /* Sparse matrix multiplication + dxdotdx_implicit += dxdotdw * dwdx */ + derived_state_.dxdotdw_.sparse_multiply( + derived_state_.dxdotdx_implicit, derived_state_.dwdx_ + ); + + derived_state_.dfdx_.sparse_add( + derived_state_.dxdotdx_explicit, 1.0, derived_state_.dxdotdx_implicit, + 1.0 + ); + fM(t, x_pos); + JSparse.sparse_add(derived_state_.MSparse_, -cj, derived_state_.dfdx_, 1.0); } void Model_DAE::fJv( @@ -169,73 +157,49 @@ void Model_DAE::fdxdotdp( auto const x_pos = computeX_pos(x); fdwdp(t, N_VGetArrayPointerConst(x_pos)); - if (pythonGenerated) { - // python generated - derived_state_.dxdotdp_explicit.zero(); - derived_state_.dxdotdp_implicit.zero(); - if (derived_state_.dxdotdp_explicit.capacity()) { - fdxdotdp_explicit_colptrs(derived_state_.dxdotdp_explicit); - fdxdotdp_explicit_rowvals(derived_state_.dxdotdp_explicit); - fdxdotdp_explicit( - derived_state_.dxdotdp_explicit.data(), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), N_VGetArrayPointerConst(dx), - derived_state_.w_.data() - ); - } - - fdxdotdw(t, x_pos, dx); - /* Sparse matrix multiplication - dxdotdp_implicit += dxdotdw * dwdp */ - derived_state_.dxdotdw_.sparse_multiply( - derived_state_.dxdotdp_implicit, derived_state_.dwdp_ - ); - - derived_state_.dxdotdp_full.sparse_add( - derived_state_.dxdotdp_explicit, 1.0, - derived_state_.dxdotdp_implicit, 1.0 + derived_state_.dxdotdp_explicit.zero(); + derived_state_.dxdotdp_implicit.zero(); + if (derived_state_.dxdotdp_explicit.capacity()) { + fdxdotdp_explicit_colptrs(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit_rowvals(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit( + derived_state_.dxdotdp_explicit.data(), t, + N_VGetArrayPointerConst(x_pos), state_.unscaledParameters.data(), + state_.fixedParameters.data(), state_.h.data(), + N_VGetArrayPointerConst(dx), derived_state_.w_.data() ); - } else { - // matlab generated - - for (int ip = 0; ip < nplist(); ip++) { - N_VConst(0.0, derived_state_.dxdotdp.getNVector(ip)); - fdxdotdp( - derived_state_.dxdotdp.data(ip), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), plist(ip), N_VGetArrayPointerConst(dx), - derived_state_.w_.data(), derived_state_.dwdp_.data() - ); - } } + + fdxdotdw(t, x_pos, dx); + /* Sparse matrix multiplication + dxdotdp_implicit += dxdotdw * dwdp */ + derived_state_.dxdotdw_.sparse_multiply( + derived_state_.dxdotdp_implicit, derived_state_.dwdp_ + ); + + derived_state_.dxdotdp_full.sparse_add( + derived_state_.dxdotdp_explicit, 1.0, derived_state_.dxdotdp_implicit, + 1.0 + ); } -void Model_DAE::fM(realtype const t, const_N_Vector x) { - if (pythonGenerated) { - /* - * non-algebraic states in python generated code always have factor - * 1 in the mass matrix, so we can easily construct this matrix here - * and avoid having to generate c++ code - */ - int ndiff = 0; - for (int ix = 0; ix < nx_solver; ix++) { - derived_state_.MSparse_.set_indexptr(ix, ndiff); - if (this->idlist.at(ix) == 1.0) { - derived_state_.MSparse_.set_data(ndiff, 1.0); - derived_state_.MSparse_.set_indexval(ndiff, ix); - ndiff++; - } +void Model_DAE::fM(realtype const /*t*/, const_N_Vector /*x*/) { + /* + * non-algebraic states in python generated code always have factor + * 1 in the mass matrix, so we can easily construct this matrix here + * and avoid having to generate c++ code + */ + int ndiff = 0; + for (int ix = 0; ix < nx_solver; ix++) { + derived_state_.MSparse_.set_indexptr(ix, ndiff); + if (this->idlist.at(ix) == 1.0) { + derived_state_.MSparse_.set_data(ndiff, 1.0); + derived_state_.MSparse_.set_indexval(ndiff, ix); + ndiff++; } - derived_state_.MSparse_.set_indexptr(nx_solver, ndiff); - assert(ndiff == derived_state_.MSparse_.capacity()); - } else { - derived_state_.M_.zero(); - auto const x_pos = computeX_pos(x); - fM(derived_state_.M_.data(), t, N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data()); } + derived_state_.MSparse_.set_indexptr(nx_solver, ndiff); + assert(ndiff == derived_state_.MSparse_.capacity()); } std::unique_ptr Model_DAE::getSolver() { @@ -521,25 +485,17 @@ void Model_DAE::fsxdot( derived_state_.J_.refresh(); } - if (pythonGenerated) { - /* copy dxdotdp and the implicit version over */ - // initialize - N_VConst(0.0, sxdot); - realtype* sxdot_tmp = N_VGetArrayPointer(sxdot); - derived_state_.dxdotdp_full.scatter( - plist(ip), 1.0, nullptr, gsl::make_span(sxdot_tmp, nx_solver), 0, - nullptr, 0 - ); - - derived_state_.J_.multiply(sxdot, sx); - derived_state_.MSparse_.multiply(sxdot, sdx, -1.0); - } else { - /* copy dxdotdp over */ - N_VScale(1.0, derived_state_.dxdotdp.getNVector(ip), sxdot); + /* copy dxdotdp and the implicit version over */ + // initialize + N_VConst(0.0, sxdot); + realtype* sxdot_tmp = N_VGetArrayPointer(sxdot); + derived_state_.dxdotdp_full.scatter( + plist(ip), 1.0, nullptr, gsl::make_span(sxdot_tmp, nx_solver), 0, + nullptr, 0 + ); - derived_state_.J_.multiply(sxdot, sx); - derived_state_.M_.multiply(sxdot, sdx, -1.0); - } + derived_state_.J_.multiply(sxdot, sx); + derived_state_.MSparse_.multiply(sxdot, sdx, -1.0); } } // namespace amici diff --git a/src/model_ode.cpp b/src/model_ode.cpp index ea76e21823..5bf3ca268c 100644 --- a/src/model_ode.cpp +++ b/src/model_ode.cpp @@ -30,40 +30,31 @@ void Model_ODE::fJSparse( void Model_ODE::fJSparse(realtype const t, const_N_Vector x, SUNMatrix J) { auto const x_pos = computeX_pos(x); fdwdx(t, N_VGetArrayPointerConst(x_pos), false); - if (pythonGenerated) { - auto JSparse = SUNMatrixWrapper(J); - // python generated - derived_state_.dxdotdx_explicit.zero(); - derived_state_.dxdotdx_implicit.zero(); - if (derived_state_.dxdotdx_explicit.capacity()) { - fdxdotdx_explicit_colptrs(derived_state_.dxdotdx_explicit); - fdxdotdx_explicit_rowvals(derived_state_.dxdotdx_explicit); - fdxdotdx_explicit( - derived_state_.dxdotdx_explicit.data(), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), derived_state_.w_.data() - ); - } - fdxdotdw(t, x_pos); - /* Sparse matrix multiplication - dxdotdx_implicit += dxdotdw * dwdx */ - derived_state_.dxdotdw_.sparse_multiply( - derived_state_.dxdotdx_implicit, derived_state_.dwdx_ - ); - - JSparse.sparse_add( - derived_state_.dxdotdx_explicit, 1.0, - derived_state_.dxdotdx_implicit, 1.0 - ); - } else { - fJSparse( - static_cast(SM_CONTENT_S(J)), t, + auto JSparse = SUNMatrixWrapper(J); + // python generated + derived_state_.dxdotdx_explicit.zero(); + derived_state_.dxdotdx_implicit.zero(); + if (derived_state_.dxdotdx_explicit.capacity()) { + fdxdotdx_explicit_colptrs(derived_state_.dxdotdx_explicit); + fdxdotdx_explicit_rowvals(derived_state_.dxdotdx_explicit); + fdxdotdx_explicit( + derived_state_.dxdotdx_explicit.data(), t, N_VGetArrayPointerConst(x_pos), state_.unscaledParameters.data(), state_.fixedParameters.data(), state_.h.data(), - derived_state_.w_.data(), derived_state_.dwdx_.data() + derived_state_.w_.data() ); } + fdxdotdw(t, x_pos); + /* Sparse matrix multiplication + dxdotdx_implicit += dxdotdw * dwdx */ + derived_state_.dxdotdw_.sparse_multiply( + derived_state_.dxdotdx_implicit, derived_state_.dwdx_ + ); + + JSparse.sparse_add( + derived_state_.dxdotdx_explicit, 1.0, derived_state_.dxdotdx_implicit, + 1.0 + ); } void Model_ODE:: @@ -148,45 +139,31 @@ void Model_ODE::fdxdotdp(realtype const t, const_N_Vector x) { auto const x_pos = computeX_pos(x); fdwdp(t, N_VGetArrayPointerConst(x_pos), false); - if (pythonGenerated) { - // python generated - derived_state_.dxdotdp_explicit.zero(); - derived_state_.dxdotdp_implicit.zero(); - if (derived_state_.dxdotdp_explicit.capacity()) { - fdxdotdp_explicit_colptrs(derived_state_.dxdotdp_explicit); - fdxdotdp_explicit_rowvals(derived_state_.dxdotdp_explicit); - fdxdotdp_explicit( - derived_state_.dxdotdp_explicit.data(), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), derived_state_.w_.data() - ); - } - - fdxdotdw(t, x_pos); - /* Sparse matrix multiplication - dxdotdp_implicit += dxdotdw * dwdp */ - derived_state_.dxdotdw_.sparse_multiply( - derived_state_.dxdotdp_implicit, derived_state_.dwdp_ - ); - - derived_state_.dxdotdp_full.sparse_add( - derived_state_.dxdotdp_explicit, 1.0, - derived_state_.dxdotdp_implicit, 1.0 + // python generated + derived_state_.dxdotdp_explicit.zero(); + derived_state_.dxdotdp_implicit.zero(); + if (derived_state_.dxdotdp_explicit.capacity()) { + fdxdotdp_explicit_colptrs(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit_rowvals(derived_state_.dxdotdp_explicit); + fdxdotdp_explicit( + derived_state_.dxdotdp_explicit.data(), t, + N_VGetArrayPointerConst(x_pos), state_.unscaledParameters.data(), + state_.fixedParameters.data(), state_.h.data(), + derived_state_.w_.data() ); - } else { - // matlab generated - for (int ip = 0; ip < nplist(); ip++) { - N_VConst(0.0, derived_state_.dxdotdp.getNVector(ip)); - fdxdotdp( - derived_state_.dxdotdp.data(ip), t, - N_VGetArrayPointerConst(x_pos), - state_.unscaledParameters.data(), state_.fixedParameters.data(), - state_.h.data(), plist(ip), derived_state_.w_.data(), - derived_state_.dwdp_.data() - ); - } } + + fdxdotdw(t, x_pos); + /* Sparse matrix multiplication + dxdotdp_implicit += dxdotdw * dwdp */ + derived_state_.dxdotdw_.sparse_multiply( + derived_state_.dxdotdp_implicit, derived_state_.dwdp_ + ); + + derived_state_.dxdotdp_full.sparse_add( + derived_state_.dxdotdp_explicit, 1.0, derived_state_.dxdotdp_implicit, + 1.0 + ); } void Model_ODE:: @@ -198,17 +175,6 @@ std::unique_ptr Model_ODE::getSolver() { return std::unique_ptr(new CVodeSolver()); } -void Model_ODE::fJSparse( - SUNMatrixContent_Sparse /*JSparse*/, realtype const /*t*/, - realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, - realtype const* /*h*/, realtype const* /*w*/, realtype const* /*dwdx*/ -) { - throw AmiException( - "Requested functionality is not supported as %s " - "is not implemented for this model!", - __func__ - ); // not implemented -} void Model_ODE::fJSparse( realtype* /*JSparse*/, realtype const /*t*/, realtype const* /*x*/, @@ -250,17 +216,6 @@ void Model_ODE::froot( ); // not implemented } -void Model_ODE::fdxdotdp( - realtype* /*dxdotdp*/, realtype const /*t*/, realtype const* /*x*/, - realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, - int const /*ip*/, realtype const* /*w*/, realtype const* /*dwdp*/ -) { - throw AmiException( - "Requested functionality is not supported as %s " - "is not implemented for this model!", - __func__ - ); // not implemented -} void Model_ODE::fdxdotdp_explicit( realtype* /*dxdotdp_explicit*/, realtype const /*t*/, realtype const* /*x*/, @@ -412,28 +367,9 @@ void Model_ODE::fqBdot( N_VConst(0.0, qBdot); fdxdotdp(t, x); - if (pythonGenerated) { - /* call multiplication */ - derived_state_.dxdotdp_full.multiply(qBdot, xB, state_.plist, true); - N_VScale(-1.0, qBdot, qBdot); - } else { - /* was matlab generated */ - for (int ip = 0; ip < nplist(); ip++) { - for (int ix = 0; ix < nxtrue_solver; ix++) - NV_Ith_S(qBdot, ip * nJ) - -= NV_Ith_S(xB, ix) * derived_state_.dxdotdp.at(ix, ip); - // second order part - for (int iJ = 1; iJ < nJ; iJ++) - for (int ix = 0; ix < nxtrue_solver; ix++) - NV_Ith_S(qBdot, ip * nJ + iJ) - -= NV_Ith_S(xB, ix) - * derived_state_.dxdotdp.at( - ix + iJ * nxtrue_solver, ip - ) - + NV_Ith_S(xB, ix + iJ * nxtrue_solver) - * derived_state_.dxdotdp.at(ix, ip); - } - } + /* call multiplication */ + derived_state_.dxdotdp_full.multiply(qBdot, xB, state_.plist, true); + N_VScale(-1.0, qBdot, qBdot); } void Model_ODE::fxBdot_ss( @@ -499,21 +435,16 @@ void Model_ODE::fsxdot( fJSparse(t, x, derived_state_.J_); derived_state_.J_.refresh(); } - if (pythonGenerated) { - /* copy dxdotdp and the implicit version over */ - // initialize - N_VConst(0.0, sxdot); - realtype* sxdot_tmp = N_VGetArrayPointer(sxdot); - - derived_state_.dxdotdp_full.scatter( - plist(ip), 1.0, nullptr, gsl::make_span(sxdot_tmp, nx_solver), 0, - nullptr, 0 - ); - } else { - /* copy dxdotdp over */ - N_VScale(1.0, derived_state_.dxdotdp.getNVector(ip), sxdot); - } + /* copy dxdotdp and the implicit version over */ + // initialize + N_VConst(0.0, sxdot); + realtype* sxdot_tmp = N_VGetArrayPointer(sxdot); + derived_state_.dxdotdp_full.scatter( + plist(ip), 1.0, nullptr, gsl::make_span(sxdot_tmp, nx_solver), 0, + nullptr, 0 + ); + derived_state_.J_.multiply(sxdot, sx); } diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index 2fa0a9adee..08b53d0249 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -64,7 +64,6 @@ void NewtonSolver::computeNewtonSensis( ); } - if (model.pythonGenerated) { for (int ip = 0; ip < model.nplist(); ip++) { N_VConst(0.0, sx.getNVector(ip)); model.get_dxdotdp_full().scatter( @@ -74,14 +73,6 @@ void NewtonSolver::computeNewtonSensis( solveLinearSystem(sx[ip]); } - } else { - for (int ip = 0; ip < model.nplist(); ip++) { - for (int ix = 0; ix < model.nx_solver; ix++) - sx.at(ix, ip) = -model.get_dxdotdp().at(ix, ip); - - solveLinearSystem(sx[ip]); - } - } } void NewtonSolver::prepareLinearSystem(Model& model, DEStateView const& state) { From 271426e92c3b762f8a10797d2285a2515a9a836e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 16:13:37 +0200 Subject: [PATCH 2/5] .. --- include/amici/abstract_model.h | 3 - include/amici/model_dimensions.h | 10 +--- include/amici/model_ode.h | 2 - src/abstract_model.cpp | 4 -- src/model.cpp | 47 ++++++++-------- src/model_ode.cpp | 2 - src/model_state.cpp | 95 +++++++++++++------------------- src/newton_solver.cpp | 16 +++--- 8 files changed, 72 insertions(+), 107 deletions(-) diff --git a/include/amici/abstract_model.h b/include/amici/abstract_model.h index 3a9072438e..c2cf6ce0f0 100644 --- a/include/amici/abstract_model.h +++ b/include/amici/abstract_model.h @@ -287,7 +287,6 @@ class AbstractModel { */ virtual void fdx0(AmiVector& x0, AmiVector& dx0); - /** * @brief Model-specific implementation of fstau * @param stau total derivative of event timepoint @@ -322,7 +321,6 @@ class AbstractModel { fy(realtype* y, realtype t, realtype const* x, realtype const* p, realtype const* k, realtype const* h, realtype const* w); - /** * @brief Model-specific implementation of fdydp * @param dydp partial derivative of observables y w.r.t. model parameters p @@ -557,7 +555,6 @@ class AbstractModel { realtype const* xB, realtype const* tcl ); - /** * @brief Model-specific implementation of fdeltaqB * @param deltaqB sensitivity update diff --git a/include/amici/model_dimensions.h b/include/amici/model_dimensions.h index 7b6ade81c5..0f883343fe 100644 --- a/include/amici/model_dimensions.h +++ b/include/amici/model_dimensions.h @@ -54,8 +54,6 @@ 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 @@ -69,8 +67,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, bool pythonGenerated = false, int ndxdotdp_explicit = 0, - int ndxdotdx_explicit = 0, int w_recursion_depth = 0 + int const lbw, int ndxdotdp_explicit = 0, int ndxdotdx_explicit = 0, + int w_recursion_depth = 0 ) : nx_rdata(nx_rdata) , nxtrue_rdata(nxtrue_rdata) @@ -99,7 +97,6 @@ struct ModelDimensions { , nJ(nJ) , ubw(ubw) , lbw(lbw) - , pythonGenerated(pythonGenerated) , ndxdotdp_explicit(ndxdotdp_explicit) , ndxdotdx_explicit(ndxdotdx_explicit) , w_recursion_depth(w_recursion_depth) { @@ -243,9 +240,6 @@ 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; diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 5e970c2342..45115aee1e 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -279,7 +279,6 @@ class Model_ODE : public Model { std::unique_ptr getSolver() override; protected: - /** * @brief Model specific implementation for fJSparse, data only * @param JSparse Matrix to which the Jacobian will be written @@ -339,7 +338,6 @@ class Model_ODE : public Model { realtype const* k, realtype const* h, realtype const* w ) = 0; - /** * @brief Model specific implementation of fdxdotdp_explicit, no w chainrule * @param dxdotdp_explicit partial derivative xdot wrt p diff --git a/src/abstract_model.cpp b/src/abstract_model.cpp index f3b5e5f09a..9edf3a9fa9 100644 --- a/src/abstract_model.cpp +++ b/src/abstract_model.cpp @@ -54,7 +54,6 @@ void AbstractModel::fdx0(AmiVector& /*x0*/, AmiVector& /*dx0*/) { // no-op default implementation } - void AbstractModel::fstau( realtype* /*stau*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, @@ -79,7 +78,6 @@ void AbstractModel:: ); } - void AbstractModel::fdydp( realtype* /*dydp*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, @@ -199,7 +197,6 @@ void AbstractModel::fdrzdx( ); } - void AbstractModel::fdeltax( realtype* /*deltax*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, @@ -242,7 +239,6 @@ void AbstractModel::fdeltaxB( ); } - void AbstractModel::fdeltaqB( realtype* /*deltaqB*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, diff --git a/src/model.cpp b/src/model.cpp index 3e1cc4e107..671c8c2184 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -2846,32 +2846,31 @@ void Model::fdwdx(realtype const t, realtype const* x, bool include_static) { fw(t, x, include_static); derived_state_.dwdx_.zero(); - fdwdw(t, x, include_static); + fdwdw(t, x, include_static); - auto&& dwdx_hierarchical_0 = derived_state_.dwdx_hierarchical_.at(0); - if (!dwdx_hierarchical_0.data() || !dwdx_hierarchical_0.capacity()) - return; + auto&& dwdx_hierarchical_0 = derived_state_.dwdx_hierarchical_.at(0); + if (!dwdx_hierarchical_0.data() || !dwdx_hierarchical_0.capacity()) + return; - if (include_static) { - dwdx_hierarchical_0.zero(); - fdwdx_colptrs(dwdx_hierarchical_0); - fdwdx_rowvals(dwdx_hierarchical_0); - } - fdwdx( - dwdx_hierarchical_0.data(), t, x, state_.unscaledParameters.data(), - state_.fixedParameters.data(), state_.h.data(), - derived_state_.w_.data(), state_.total_cl.data(), - derived_state_.spl_.data(), include_static - ); - - for (int irecursion = 1; irecursion <= w_recursion_depth; - irecursion++) { - derived_state_.dwdw_.sparse_multiply( - derived_state_.dwdx_hierarchical_.at(irecursion), - derived_state_.dwdx_hierarchical_.at(irecursion - 1) - ); - } - derived_state_.dwdx_.sparse_sum(derived_state_.dwdx_hierarchical_); + if (include_static) { + dwdx_hierarchical_0.zero(); + fdwdx_colptrs(dwdx_hierarchical_0); + fdwdx_rowvals(dwdx_hierarchical_0); + } + fdwdx( + dwdx_hierarchical_0.data(), t, x, state_.unscaledParameters.data(), + state_.fixedParameters.data(), state_.h.data(), + derived_state_.w_.data(), state_.total_cl.data(), + derived_state_.spl_.data(), include_static + ); + + for (int irecursion = 1; irecursion <= w_recursion_depth; irecursion++) { + derived_state_.dwdw_.sparse_multiply( + derived_state_.dwdx_hierarchical_.at(irecursion), + derived_state_.dwdx_hierarchical_.at(irecursion - 1) + ); + } + derived_state_.dwdx_.sparse_sum(derived_state_.dwdx_hierarchical_); if (always_check_finite_) { checkFinite(derived_state_.dwdx_, ModelQuantity::dwdx, t); diff --git a/src/model_ode.cpp b/src/model_ode.cpp index 5bf3ca268c..7cd01fc454 100644 --- a/src/model_ode.cpp +++ b/src/model_ode.cpp @@ -175,7 +175,6 @@ std::unique_ptr Model_ODE::getSolver() { return std::unique_ptr(new CVodeSolver()); } - void Model_ODE::fJSparse( realtype* /*JSparse*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, @@ -216,7 +215,6 @@ void Model_ODE::froot( ); // not implemented } - void Model_ODE::fdxdotdp_explicit( realtype* /*dxdotdp_explicit*/, realtype const /*t*/, realtype const* /*x*/, realtype const* /*p*/, realtype const* /*k*/, realtype const* /*h*/, diff --git a/src/model_state.cpp b/src/model_state.cpp index 9fb1370c04..d3caa41844 100644 --- a/src/model_state.cpp +++ b/src/model_state.cpp @@ -23,69 +23,52 @@ ModelStateDerived::ModelStateDerived(ModelDimensions const& dim) , x_rdata_(dim.nx_rdata, 0.0) , sx_rdata_(dim.nx_rdata, 0.0) , x_pos_tmp_(dim.nx_solver, sunctx_) { - // 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, sunctx_); + // size dynamically adapted for dwdx_ and dwdp_ + dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, 0, CSC_MAT, sunctx_); + dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, 0, CSC_MAT, sunctx_); - dwdw_ = SUNMatrixWrapper(dim.nw, dim.nw, dim.ndwdw, CSC_MAT, sunctx_); - // size dynamically adapted for dwdx_ and dwdp_ - dwdx_ = SUNMatrixWrapper(dim.nw, dim.nx_solver, 0, CSC_MAT, sunctx_); - dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, 0, CSC_MAT, sunctx_); - - 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( - dim.nw, dim.np, irec * dim.ndwdw + dim.ndwdp, CSC_MAT, sunctx_ - ); - dwdx_hierarchical_.emplace_back( - dim.nw, dim.nx_solver, irec * dim.ndwdw + dim.ndwdx, CSC_MAT, - sunctx_ - ); - } - 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, sunctx_ - ); - // guess size, will be dynamically reallocated - dxdotdp_implicit = SUNMatrixWrapper( - dim.nx_solver, dim.np, dim.ndwdp + dim.ndxdotdw, CSC_MAT, sunctx_ + 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( + dim.nw, dim.np, irec * dim.ndwdw + dim.ndwdp, CSC_MAT, sunctx_ ); - dxdotdx_explicit = SUNMatrixWrapper( - dim.nx_solver, dim.nx_solver, dim.ndxdotdx_explicit, CSC_MAT, + dwdx_hierarchical_.emplace_back( + dim.nw, dim.nx_solver, irec * dim.ndwdw + dim.ndwdx, CSC_MAT, sunctx_ ); - // guess size, will be dynamically reallocated - dxdotdx_implicit = SUNMatrixWrapper( - dim.nx_solver, dim.nx_solver, dim.ndwdx + dim.ndxdotdw, CSC_MAT, - sunctx_ - ); - // dynamically allocate on first call - dxdotdp_full - = SUNMatrixWrapper(dim.nx_solver, dim.np, 0, CSC_MAT, sunctx_); + } + assert( + gsl::narrow(dwdp_hierarchical_.size()) == dim.w_recursion_depth + 1 + ); + assert( + gsl::narrow(dwdx_hierarchical_.size()) == dim.w_recursion_depth + 1 + ); - for (int iytrue = 0; iytrue < dim.nytrue; ++iytrue) - dJydy_.emplace_back( - dim.nJ, dim.ny, dim.ndJydy.at(iytrue), CSC_MAT, sunctx_ - ); + dxdotdp_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndxdotdp_explicit, CSC_MAT, sunctx_ + ); + // guess size, will be dynamically reallocated + dxdotdp_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.np, dim.ndwdp + dim.ndxdotdw, CSC_MAT, sunctx_ + ); + dxdotdx_explicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndxdotdx_explicit, CSC_MAT, sunctx_ + ); + // guess size, will be dynamically reallocated + dxdotdx_implicit = SUNMatrixWrapper( + dim.nx_solver, dim.nx_solver, dim.ndwdx + dim.ndxdotdw, CSC_MAT, sunctx_ + ); + // dynamically allocate on first call + dxdotdp_full = SUNMatrixWrapper(dim.nx_solver, dim.np, 0, CSC_MAT, sunctx_); - dJydy_dense_ = SUNMatrixWrapper(dim.nJ, dim.ny, sunctx_); - } else { - dwdx_ = SUNMatrixWrapper( - dim.nw, dim.nx_solver, dim.ndwdx, CSC_MAT, sunctx_ + for (int iytrue = 0; iytrue < dim.nytrue; ++iytrue) + dJydy_.emplace_back( + dim.nJ, dim.ny, dim.ndJydy.at(iytrue), CSC_MAT, sunctx_ ); - dwdp_ = SUNMatrixWrapper(dim.nw, dim.np, dim.ndwdp, CSC_MAT, sunctx_); - dJydy_matlab_ - = std::vector(dim.nJ * dim.nytrue * dim.ny, 0.0); - } + + dJydy_dense_ = SUNMatrixWrapper(dim.nJ, dim.ny, sunctx_); } } // namespace amici diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index 08b53d0249..cca5b6cb35 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -64,15 +64,15 @@ void NewtonSolver::computeNewtonSensis( ); } - for (int ip = 0; ip < model.nplist(); ip++) { - N_VConst(0.0, sx.getNVector(ip)); - model.get_dxdotdp_full().scatter( - model.plist(ip), -1.0, nullptr, - gsl::make_span(sx.getNVector(ip)), 0, nullptr, 0 - ); + for (int ip = 0; ip < model.nplist(); ip++) { + N_VConst(0.0, sx.getNVector(ip)); + model.get_dxdotdp_full().scatter( + model.plist(ip), -1.0, nullptr, gsl::make_span(sx.getNVector(ip)), + 0, nullptr, 0 + ); - solveLinearSystem(sx[ip]); - } + solveLinearSystem(sx[ip]); + } } void NewtonSolver::prepareLinearSystem(Model& model, DEStateView const& state) { From 7b775528e1e171f2815b889634e262fb3b142dbb Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 16:14:44 +0200 Subject: [PATCH 3/5] .. --- include/amici/model_state.h | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/include/amici/model_state.h b/include/amici/model_state.h index cc857a157c..ee566d862c 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -93,7 +93,6 @@ struct ModelStateDerived { , dxdotdw_(other.dxdotdw_) , dwdx_(other.dwdx_) , dwdp_(other.dwdp_) - , M_(other.M_) , MSparse_(other.MSparse_) , dfdx_(other.dfdx_) , dxdotdp_full(other.dxdotdp_full) @@ -106,7 +105,6 @@ struct ModelStateDerived { , dtotal_cldx_rdata(other.dtotal_cldx_rdata) , dxdotdp(other.dxdotdp) , dJydy_(other.dJydy_) - , dJydy_matlab_(other.dJydy_matlab_) , dJydsigma_(other.dJydsigma_) , dJydx_(other.dJydx_) , dJydp_(other.dJydp_) @@ -152,7 +150,6 @@ struct ModelStateDerived { dxdotdw_.set_ctx(sunctx_); dwdx_.set_ctx(sunctx_); dwdp_.set_ctx(sunctx_); - M_.set_ctx(sunctx_); MSparse_.set_ctx(sunctx_); dfdx_.set_ctx(sunctx_); dxdotdp_full.set_ctx(sunctx_); @@ -210,13 +207,6 @@ struct ModelStateDerived { /** Sparse dwdp temporary storage (dimension: `nw` x `np`, nnz: `ndwdp`) */ SUNMatrixWrapper dwdp_; - /** - * Dense Mass matrix (dimension: `nx_solver` x `nx_solver`) - * - * MATLAB-generated-only, DEPRECATED. - */ - SUNMatrixWrapper M_; - /** * Sparse Mass matrix, Python-generated-only * @@ -296,12 +286,6 @@ struct ModelStateDerived { */ 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) */ From f9bfeab056c7eec5243c00319aef21f1ff9f1fc2 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 16:21:33 +0200 Subject: [PATCH 4/5] .. --- include/amici/serialization.h | 1 - models/model_calvetti_py/model_calvetti_py.h | 3 +-- models/model_dirac_py/model_dirac_py.h | 3 +-- models/model_events_py/model_events_py.h | 3 +-- models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h | 3 +-- models/model_nested_events_py/model_nested_events_py.h | 3 +-- models/model_neuron_py/model_neuron_py.h | 3 +-- models/model_robertson_py/model_robertson_py.h | 3 +-- models/model_steadystate_py/model_steadystate_py.h | 3 +-- src/model_header.template.h | 1 - 10 files changed, 8 insertions(+), 18 deletions(-) diff --git a/include/amici/serialization.h b/include/amici/serialization.h index 5d4466650c..d2240266fc 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -149,7 +149,6 @@ void serialize(Archive& ar, amici::Model& m, unsigned int const /*version*/) { ar & m.sx0data_; ar & m.nmaxevent_; ar & m.state_is_non_negative_; - ar & m.pythonGenerated; ar & m.min_sigma_; ar & m.sigma_res_; ar & m.steadystate_computation_mode_; diff --git a/models/model_calvetti_py/model_calvetti_py.h b/models/model_calvetti_py/model_calvetti_py.h index a0e20af341..846a64e0f9 100644 --- a/models/model_calvetti_py/model_calvetti_py.h +++ b/models/model_calvetti_py/model_calvetti_py.h @@ -138,7 +138,6 @@ class Model_model_calvetti_py : public amici::Model_DAE { 0, // nnz 6, // ubw 6, // lbw - true, // pythonGenerated 0, // ndxdotdp_explicit 5, // ndxdotdx_explicit 2 // w_recursion_depth @@ -558,7 +557,7 @@ class Model_model_calvetti_py : public amici::Model_DAE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_dirac_py/model_dirac_py.h b/models/model_dirac_py/model_dirac_py.h index 7a7d4a1119..3858d2625d 100644 --- a/models/model_dirac_py/model_dirac_py.h +++ b/models/model_dirac_py/model_dirac_py.h @@ -138,7 +138,6 @@ class Model_model_dirac_py : public amici::Model_ODE { 0, // nnz 2, // ubw 2, // lbw - true, // pythonGenerated 3, // ndxdotdp_explicit 3, // ndxdotdx_explicit 0 // w_recursion_depth @@ -545,7 +544,7 @@ class Model_model_dirac_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_events_py/model_events_py.h b/models/model_events_py/model_events_py.h index 114369a530..fa7415dcc1 100644 --- a/models/model_events_py/model_events_py.h +++ b/models/model_events_py/model_events_py.h @@ -138,7 +138,6 @@ class Model_model_events_py : public amici::Model_ODE { 0, // nnz 3, // ubw 3, // lbw - true, // pythonGenerated 3, // ndxdotdp_explicit 4, // ndxdotdx_explicit 0 // w_recursion_depth @@ -580,7 +579,7 @@ class Model_model_events_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h index 6dbcb93a04..89bb876f9f 100644 --- a/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h +++ b/models/model_jakstat_adjoint_py/model_jakstat_adjoint_py.h @@ -138,7 +138,6 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE { 0, // nnz 9, // ubw 9, // lbw - true, // pythonGenerated 13, // ndxdotdp_explicit 18, // ndxdotdx_explicit 0 // w_recursion_depth @@ -553,7 +552,7 @@ class Model_model_jakstat_adjoint_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_nested_events_py/model_nested_events_py.h b/models/model_nested_events_py/model_nested_events_py.h index 4f5c25497a..e56c7c6a92 100644 --- a/models/model_nested_events_py/model_nested_events_py.h +++ b/models/model_nested_events_py/model_nested_events_py.h @@ -138,7 +138,6 @@ class Model_model_nested_events_py : public amici::Model_ODE { 0, // nnz 1, // ubw 1, // lbw - true, // pythonGenerated 2, // ndxdotdp_explicit 1, // ndxdotdx_explicit 0 // w_recursion_depth @@ -553,7 +552,7 @@ class Model_model_nested_events_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_neuron_py/model_neuron_py.h b/models/model_neuron_py/model_neuron_py.h index 5b95c92f8d..26a35fcf97 100644 --- a/models/model_neuron_py/model_neuron_py.h +++ b/models/model_neuron_py/model_neuron_py.h @@ -138,7 +138,6 @@ class Model_model_neuron_py : public amici::Model_ODE { 0, // nnz 2, // ubw 2, // lbw - true, // pythonGenerated 2, // ndxdotdp_explicit 4, // ndxdotdx_explicit 0 // w_recursion_depth @@ -575,7 +574,7 @@ class Model_model_neuron_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_robertson_py/model_robertson_py.h b/models/model_robertson_py/model_robertson_py.h index 64f015a80c..6f24fe3ca3 100644 --- a/models/model_robertson_py/model_robertson_py.h +++ b/models/model_robertson_py/model_robertson_py.h @@ -138,7 +138,6 @@ class Model_model_robertson_py : public amici::Model_DAE { 0, // nnz 3, // ubw 3, // lbw - true, // pythonGenerated 5, // ndxdotdp_explicit 9, // ndxdotdx_explicit 0 // w_recursion_depth @@ -537,7 +536,7 @@ class Model_model_robertson_py : public amici::Model_DAE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/models/model_steadystate_py/model_steadystate_py.h b/models/model_steadystate_py/model_steadystate_py.h index 04ee22b23e..a67c200c1f 100644 --- a/models/model_steadystate_py/model_steadystate_py.h +++ b/models/model_steadystate_py/model_steadystate_py.h @@ -138,7 +138,6 @@ class Model_model_steadystate_py : public amici::Model_ODE { 0, // nnz 3, // ubw 3, // lbw - true, // pythonGenerated 11, // ndxdotdp_explicit 9, // ndxdotdx_explicit 0 // w_recursion_depth @@ -537,7 +536,7 @@ class Model_model_steadystate_py : public amici::Model_ODE { * @return AMICI git commit hash */ std::string getAmiciCommit() const override { - return "b12c68a7a02f1cbd33de59d47cbc0c4d77d30d6f"; + return "dcbbae0f8e9b52afd7ab8b4b7a046da36b907139"; } bool hasQuadraticLLH() const override { diff --git a/src/model_header.template.h b/src/model_header.template.h index 118b1be8bb..2cd1e1c67f 100644 --- a/src/model_header.template.h +++ b/src/model_header.template.h @@ -138,7 +138,6 @@ class Model_TPL_MODELNAME : public amici::Model_TPL_MODEL_TYPE_UPPER { 0, // nnz TPL_UBW, // ubw TPL_LBW, // lbw - true, // pythonGenerated TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit TPL_NDXDOTDX_EXPLICIT, // ndxdotdx_explicit TPL_W_RECURSION_DEPTH // w_recursion_depth From 6ff4a21ca82ceb6756e902de074312ad16b22aaa Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 17:08:26 +0200 Subject: [PATCH 5/5] fix test models, format --- cmake/clang-tools.cmake | 2 +- tests/cpp/testfunctions.h | 162 ++-- tests/cpp/unittests/testExpData.cpp | 269 +++--- tests/cpp/unittests/testMisc.cpp | 603 +++++++------ tests/cpp/unittests/testSerialization.cpp | 220 ++--- tests/cpp/unittests/testSplines.cpp | 976 +++++++++------------- 6 files changed, 1013 insertions(+), 1219 deletions(-) diff --git a/cmake/clang-tools.cmake b/cmake/clang-tools.cmake index 1fcfd760da..645f8be7fd 100644 --- a/cmake/clang-tools.cmake +++ b/cmake/clang-tools.cmake @@ -4,7 +4,7 @@ execute_process( COMMAND sh -c - "git ls-tree -r HEAD --name-only src/*.cpp include/amici/*.h | tr '\n' ' '" + "git ls-tree -r HEAD --name-only src/*.cpp include/amici/*.h tests/cpp/*.h tests/cpp/unittests/*.cpp | tr '\n' ' '" WORKING_DIRECTORY ${CMAKE_SOURCE_DIR} OUTPUT_VARIABLE ALL_CXX_SOURCE_FILES) # ########### clang-tidy ############ diff --git a/tests/cpp/testfunctions.h b/tests/cpp/testfunctions.h index ecb3ea933a..fa369fc496 100644 --- a/tests/cpp/testfunctions.h +++ b/tests/cpp/testfunctions.h @@ -33,8 +33,7 @@ std::vector getVariableNames(std::string const& name, int length); * @brief The Model_Test class is a model-unspecific implementation of model designed for unit testing. */ -class Model_Test : public Model -{ +class Model_Test : public Model { public: /** constructor with model dimensions * @param model_dimensions ModelDimensions model_dimensions, @@ -43,89 +42,96 @@ class Model_Test : public Model * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events */ - Model_Test(ModelDimensions const& model_dimensions, - SimulationParameters simulation_parameters, - const SecondOrderMode o2mode, - const std::vector idlist, - const std::vector z2event) - : Model(model_dimensions, - simulation_parameters, - o2mode, - idlist, - z2event) - {} + Model_Test( + ModelDimensions const& model_dimensions, + SimulationParameters simulation_parameters, + SecondOrderMode const o2mode, std::vector const idlist, + std::vector const z2event, std::vector const events + ) + : Model( + model_dimensions, simulation_parameters, o2mode, idlist, z2event, + events + ) {} /** default constructor */ Model_Test() : Model( - ModelDimensions(), - SimulationParameters(), - SecondOrderMode::none, - std::vector(), - std::vector()) - {} + ModelDimensions(), SimulationParameters(), SecondOrderMode::none, + std::vector(), std::vector() + ) {} - Model *clone() const override { return new Model_Test(*this); } + Model* clone() const override { return new Model_Test(*this); } std::unique_ptr getSolver() override { throw AmiException("not implemented"); } - void froot(const realtype /*t*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/, - gsl::span /*root*/) override { + void froot( + realtype const /*t*/, AmiVector const& /*x*/, AmiVector const& /*dx*/, + gsl::span /*root*/ + ) override { throw AmiException("not implemented"); } - void fxdot(const realtype /*t*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/, AmiVector & /*xdot*/) override { + void fxdot( + realtype const /*t*/, AmiVector const& /*x*/, AmiVector const& /*dx*/, + AmiVector& /*xdot*/ + ) override { throw AmiException("not implemented"); } - void fsxdot(const realtype /*t*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/, const int /*ip*/, - const AmiVector & /*sx*/, const AmiVector & /*sdx*/, - AmiVector & /*sxdot*/) override { + void fsxdot( + realtype const /*t*/, AmiVector const& /*x*/, AmiVector const& /*dx*/, + int const /*ip*/, AmiVector const& /*sx*/, AmiVector const& /*sdx*/, + AmiVector& /*sxdot*/ + ) override { throw AmiException("not implemented"); } - void fJ(const realtype /*t*/, const realtype /*cj*/, - const AmiVector & /*x*/, const AmiVector & /*dx*/, - const AmiVector & /*xdot*/, SUNMatrix /*J*/) override { + void + fJ(realtype const /*t*/, realtype const /*cj*/, AmiVector const& /*x*/, + AmiVector const& /*dx*/, AmiVector const& /*xdot*/, + SUNMatrix /*J*/) override { throw AmiException("not implemented"); } - void fJSparse(const realtype /*t*/, const realtype /*cj*/, - const AmiVector & /*x*/, const AmiVector & /*dx*/, - const AmiVector & /*xdot*/, SUNMatrix /*J*/) override { + void fJSparse( + realtype const /*t*/, realtype const /*cj*/, AmiVector const& /*x*/, + AmiVector const& /*dx*/, AmiVector const& /*xdot*/, SUNMatrix /*J*/ + ) override { throw AmiException("not implemented"); } - void fJB(const realtype /*t*/, realtype /*cj*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/, const AmiVector & /*xB*/, - const AmiVector & /*dxB*/, const AmiVector & /*xBdot*/, - SUNMatrix /*JB*/) override { + void + fJB(realtype const /*t*/, realtype /*cj*/, AmiVector const& /*x*/, + AmiVector const& /*dx*/, AmiVector const& /*xB*/, + AmiVector const& /*dxB*/, AmiVector const& /*xBdot*/, + SUNMatrix /*JB*/) override { throw AmiException("not implemented"); } - void fJSparseB(const realtype /*t*/, realtype /*cj*/, - const AmiVector & /*x*/, const AmiVector & /*dx*/, - const AmiVector & /*xB*/, const AmiVector & /*dxB*/, - const AmiVector & /*xBdot*/, SUNMatrix /*JB*/) override { + void fJSparseB( + realtype const /*t*/, realtype /*cj*/, AmiVector const& /*x*/, + AmiVector const& /*dx*/, AmiVector const& /*xB*/, + AmiVector const& /*dxB*/, AmiVector const& /*xBdot*/, SUNMatrix /*JB*/ + ) override { throw AmiException("not implemented"); } - void fJDiag(const realtype /*t*/, AmiVector & /*Jdiag*/, - const realtype /*cj*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/) override { + void fJDiag( + realtype const /*t*/, AmiVector& /*Jdiag*/, realtype const /*cj*/, + AmiVector const& /*x*/, AmiVector const& /*dx*/ + ) override { throw AmiException("not implemented"); } - void fdxdotdp(const realtype /*t*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/) override { + void fdxdotdp( + realtype const /*t*/, AmiVector const& /*x*/, AmiVector const& /*dx*/ + ) override { throw AmiException("not implemented"); } - void fJv(const realtype /*t*/, const AmiVector & /*x*/, - const AmiVector & /*dx*/, const AmiVector & /*xdot*/, - const AmiVector & /*v*/, AmiVector & /*nJv*/, - const realtype /*cj*/) override { + void + fJv(realtype const /*t*/, AmiVector const& /*x*/, AmiVector const& /*dx*/, + AmiVector const& /*xdot*/, AmiVector const& /*v*/, AmiVector& /*nJv*/, + realtype const /*cj*/) override { throw AmiException("not implemented"); } - void fxBdot_ss(const realtype /*t*/, const AmiVector & /*xB*/, - const AmiVector & /*dxB*/, AmiVector & /*xBdot*/ - ) override { + void fxBdot_ss( + realtype const /*t*/, AmiVector const& /*xB*/, AmiVector const& /*dxB*/, + AmiVector& /*xBdot*/ + ) override { throw AmiException("not implemented"); } @@ -133,10 +139,11 @@ class Model_Test : public Model throw AmiException("not implemented"); } - void writeSteadystateJB(const realtype /*t*/, realtype /*cj*/, - const AmiVector & /*x*/, const AmiVector & /*dx*/, - const AmiVector & /*xB*/, const AmiVector & /*dxB*/, - const AmiVector & /*xBdot*/) override { + void writeSteadystateJB( + realtype const /*t*/, realtype /*cj*/, AmiVector const& /*x*/, + AmiVector const& /*dx*/, AmiVector const& /*xB*/, + AmiVector const& /*dxB*/, AmiVector const& /*xBdot*/ + ) override { throw AmiException("not implemented"); } @@ -175,16 +182,17 @@ class Model_Test : public Model void simulateWithDefaultOptions(); -void simulateVerifyWrite(const std::string &path); +void simulateVerifyWrite(std::string const& path); -void simulateVerifyWrite(std::string const &path, double atol, double rtol); +void simulateVerifyWrite(std::string const& path, double atol, double rtol); -void simulateVerifyWrite(const std::string &hdffileOptions, - const std::string &hdffileResults, - const std::string &hdffilewrite, - const std::string &path, double atol, double rtol); +void simulateVerifyWrite( + std::string const& hdffileOptions, std::string const& hdffileResults, + std::string const& hdffilewrite, std::string const& path, double atol, + double rtol +); -std::unique_ptr getTestExpData(const Model &model); +std::unique_ptr getTestExpData(Model const& model); bool withinTolerance( double expected, double actual, double atol, double rtol, int index, @@ -196,18 +204,20 @@ void checkEqualArray( double rtol, std::string_view name ); -void checkEqualArray(std::vector const &expected, - std::vector const &actual, double atol, - double rtol, std::string const &name); +void checkEqualArray( + std::vector const& expected, std::vector const& actual, + double atol, double rtol, std::string const& name +); -void verifyReturnData(const std::string &hdffile, const std::string &resultPath, - const ReturnData *rdata, const Model *model, double atol, - double rtol); +void verifyReturnData( + std::string const& hdffile, std::string const& resultPath, + ReturnData const* rdata, Model const* model, double atol, double rtol +); -void verifyReturnDataSensitivities(const H5::H5File &file, - const std::string &resultPath, - const ReturnData *rdata, const Model *model, - double atol, double rtol); +void verifyReturnDataSensitivities( + const H5::H5File& file, std::string const& resultPath, + ReturnData const* rdata, Model const* model, double atol, double rtol +); void printBacktrace(int depth); diff --git a/tests/cpp/unittests/testExpData.cpp b/tests/cpp/unittests/testExpData.cpp index 649f7bb8c1..721356cdba 100644 --- a/tests/cpp/unittests/testExpData.cpp +++ b/tests/cpp/unittests/testExpData.cpp @@ -2,8 +2,8 @@ #include #include -#include #include +#include #include #include @@ -18,7 +18,7 @@ std::unique_ptr getModel(); using namespace amici; -namespace { +namespace { class ExpDataTest : public ::testing::Test { protected: @@ -30,7 +30,7 @@ class ExpDataTest : public ::testing::Test { } int nx = 1, ny = 2, nz = 3, nmaxevent = 4; - std::vector timepoints = { 1, 2, 3, 4 }; + std::vector timepoints = {1, 2, 3, 4}; std::unique_ptr model = generic_model::getModel(); @@ -41,8 +41,8 @@ class ExpDataTest : public ::testing::Test { nx, // nx_solver nx, // nxtrue_solver 0, // nx_solver_reinit - 1, // np - 3, // nk + 1, // np + 3, // nk ny, // ny ny, // nytrue nz, // nz @@ -56,136 +56,125 @@ class ExpDataTest : public ::testing::Test { 0, // ndwdp 0, // dwdw 0, // ndxdotdw - {}, // ndJydy + {0, 0}, // ndJydy 0, // ndxrdatadxsolver 0, // ndxrdatadtcl 0, // ndtotal_cldx_rdata 0, // nnz 0, // ubw 0 // lbw - ), + ), SimulationParameters( - std::vector(3, 0.0), - std::vector(1, 0.0), + std::vector(3, 0.0), std::vector(1, 0.0), std::vector(2, 1) - ), - SecondOrderMode::none, - std::vector(), - std::vector()); + ), + SecondOrderMode::none, std::vector(), std::vector(), + { + Event("e1", true, true, 1), + Event("e1", true, true, 1), + Event("e1", true, true, 1), + Event("e1", true, true, 1), + } + ); }; -TEST_F(ExpDataTest, DefaultConstructable) -{ +TEST_F(ExpDataTest, DefaultConstructable) { ExpData edata{}; ASSERT_EQ(edata.nytrue(), 0); ASSERT_EQ(edata.nztrue(), 0); ASSERT_EQ(edata.nmaxevent(), 0); } -TEST_F(ExpDataTest, ModelCtor) -{ +TEST_F(ExpDataTest, ModelCtor) { ExpData edata(model->nytrue, model->nztrue, model->nMaxEvent()); ASSERT_EQ(edata.nytrue(), model->nytrue); ASSERT_EQ(edata.nztrue(), model->nztrue); ASSERT_EQ(edata.nmaxevent(), model->nMaxEvent()); } -TEST_F(ExpDataTest, DimensionCtor) -{ +TEST_F(ExpDataTest, DimensionCtor) { ExpData edata(model->nytrue, model->nztrue, model->nMaxEvent(), timepoints); ASSERT_EQ(edata.nytrue(), model->nytrue); ASSERT_EQ(edata.nztrue(), model->nztrue); ASSERT_EQ(edata.nmaxevent(), model->nMaxEvent()); ASSERT_EQ(edata.nt(), model->nt()); checkEqualArray( - timepoints, edata.getTimepoints(), TEST_ATOL, TEST_RTOL, "ts"); + timepoints, edata.getTimepoints(), TEST_ATOL, TEST_RTOL, "ts" + ); } -TEST_F(ExpDataTest, MeasurementCtor) -{ +TEST_F(ExpDataTest, MeasurementCtor) { std::vector y(ny * timepoints.size(), 0.0); std::vector y_std(ny * timepoints.size(), 0.1); std::vector z(nz * nmaxevent, 0.0); std::vector z_std(nz * nmaxevent, 0.1); - ExpData edata(testModel.nytrue, - testModel.nztrue, - testModel.nMaxEvent(), - timepoints, - y, - y_std, - z, - z_std); + ExpData edata( + testModel.nytrue, testModel.nztrue, testModel.nMaxEvent(), timepoints, + y, y_std, z, z_std + ); ASSERT_EQ(edata.nytrue(), testModel.nytrue); ASSERT_EQ(edata.nztrue(), testModel.nztrue); ASSERT_EQ(edata.nmaxevent(), testModel.nMaxEvent()); ASSERT_EQ(edata.nt(), testModel.nt()); checkEqualArray( - timepoints, edata.getTimepoints(), TEST_ATOL, TEST_RTOL, "ts"); + timepoints, edata.getTimepoints(), TEST_ATOL, TEST_RTOL, "ts" + ); + checkEqualArray( + y, edata.getObservedData(), TEST_ATOL, TEST_RTOL, "observedData" + ); + checkEqualArray( + y_std, edata.getObservedDataStdDev(), TEST_ATOL, TEST_RTOL, + "observedDataStdDev" + ); checkEqualArray( - y, edata.getObservedData(), TEST_ATOL, TEST_RTOL, "observedData"); - checkEqualArray(y_std, - edata.getObservedDataStdDev(), - TEST_ATOL, - TEST_RTOL, - "observedDataStdDev"); + z, edata.getObservedEvents(), TEST_ATOL, TEST_RTOL, "observedEvents" + ); checkEqualArray( - z, edata.getObservedEvents(), TEST_ATOL, TEST_RTOL, "observedEvents"); - checkEqualArray(z_std, - edata.getObservedEventsStdDev(), - TEST_ATOL, - TEST_RTOL, - "observedEventsStdDev"); + z_std, edata.getObservedEventsStdDev(), TEST_ATOL, TEST_RTOL, + "observedEventsStdDev" + ); ExpData edata_copy(edata); ASSERT_EQ(edata.nytrue(), edata_copy.nytrue()); ASSERT_EQ(edata.nztrue(), edata_copy.nztrue()); ASSERT_EQ(edata.nmaxevent(), edata_copy.nmaxevent()); ASSERT_EQ(edata.nt(), edata_copy.nt()); - checkEqualArray(edata_copy.getTimepoints(), - edata.getTimepoints(), - TEST_ATOL, - TEST_RTOL, - "ts"); - checkEqualArray(edata_copy.getObservedData(), - edata.getObservedData(), - TEST_ATOL, - TEST_RTOL, - "observedData"); - checkEqualArray(edata_copy.getObservedDataStdDev(), - edata.getObservedDataStdDev(), - TEST_ATOL, - TEST_RTOL, - "observedDataStdDev"); - checkEqualArray(edata_copy.getObservedEvents(), - edata.getObservedEvents(), - TEST_ATOL, - TEST_RTOL, - "observedEvents"); - checkEqualArray(edata_copy.getObservedEventsStdDev(), - edata.getObservedEventsStdDev(), - TEST_ATOL, - TEST_RTOL, - "observedEventsStdDev"); + checkEqualArray( + edata_copy.getTimepoints(), edata.getTimepoints(), TEST_ATOL, TEST_RTOL, + "ts" + ); + checkEqualArray( + edata_copy.getObservedData(), edata.getObservedData(), TEST_ATOL, + TEST_RTOL, "observedData" + ); + checkEqualArray( + edata_copy.getObservedDataStdDev(), edata.getObservedDataStdDev(), + TEST_ATOL, TEST_RTOL, "observedDataStdDev" + ); + checkEqualArray( + edata_copy.getObservedEvents(), edata.getObservedEvents(), TEST_ATOL, + TEST_RTOL, "observedEvents" + ); + checkEqualArray( + edata_copy.getObservedEventsStdDev(), edata.getObservedEventsStdDev(), + TEST_ATOL, TEST_RTOL, "observedEventsStdDev" + ); } -TEST_F(ExpDataTest, CopyConstructable) -{ +TEST_F(ExpDataTest, CopyConstructable) { testModel.setTimepoints(timepoints); auto edata = ExpData(testModel); ASSERT_EQ(edata.nytrue(), testModel.nytrue); ASSERT_EQ(edata.nztrue(), testModel.nztrue); ASSERT_EQ(edata.nmaxevent(), testModel.nMaxEvent()); ASSERT_EQ(edata.nt(), testModel.nt()); - checkEqualArray(testModel.getTimepoints(), - edata.getTimepoints(), - TEST_ATOL, - TEST_RTOL, - "ts"); + checkEqualArray( + testModel.getTimepoints(), edata.getTimepoints(), TEST_ATOL, TEST_RTOL, + "ts" + ); } - -TEST_F(ExpDataTest, Equality) -{ +TEST_F(ExpDataTest, Equality) { auto edata = ExpData(testModel); auto edata2(edata); ASSERT_TRUE(edata == edata2); @@ -194,33 +183,28 @@ TEST_F(ExpDataTest, Equality) ASSERT_FALSE(edata == edata2); } -TEST_F(ExpDataTest, DimensionChecks) -{ +TEST_F(ExpDataTest, DimensionChecks) { std::vector bad_std(ny, -0.1); std::vector y(ny * timepoints.size(), 0.0); std::vector y_std(ny * timepoints.size(), 0.1); std::vector z(nz * nmaxevent, 0.0); std::vector z_std(nz * nmaxevent, 0.1); - ASSERT_THROW(ExpData(testModel.nytrue, - testModel.nztrue, - testModel.nMaxEvent(), - timepoints, - z, - z_std, - z, - z_std), - AmiException); - - ASSERT_THROW(ExpData(testModel.nytrue, - testModel.nztrue, - testModel.nMaxEvent(), - timepoints, - z, - bad_std, - z, - z_std), - AmiException); + ASSERT_THROW( + ExpData( + testModel.nytrue, testModel.nztrue, testModel.nMaxEvent(), + timepoints, z, z_std, z, z_std + ), + AmiException + ); + + ASSERT_THROW( + ExpData( + testModel.nytrue, testModel.nztrue, testModel.nMaxEvent(), + timepoints, z, bad_std, z, z_std + ), + AmiException + ); ExpData edata(testModel); @@ -239,21 +223,21 @@ TEST_F(ExpDataTest, DimensionChecks) std::vector bad_single_z(edata.nmaxevent() + 1, 0.0); std::vector bad_single_z_std(edata.nmaxevent() + 1, 0.1); - ASSERT_THROW(edata.setObservedData(bad_single_y, 0), - AmiException); - ASSERT_THROW(edata.setObservedDataStdDev(bad_single_y_std, 0), - AmiException); - ASSERT_THROW(edata.setObservedEvents(bad_single_z, 0), - AmiException); - ASSERT_THROW(edata.setObservedEventsStdDev(bad_single_y_std, 0), - AmiException); - - ASSERT_THROW(edata.setTimepoints(std::vector{ 0.0, 1.0, 0.5 }), - AmiException); + ASSERT_THROW(edata.setObservedData(bad_single_y, 0), AmiException); + ASSERT_THROW( + edata.setObservedDataStdDev(bad_single_y_std, 0), AmiException + ); + ASSERT_THROW(edata.setObservedEvents(bad_single_z, 0), AmiException); + ASSERT_THROW( + edata.setObservedEventsStdDev(bad_single_y_std, 0), AmiException + ); + + ASSERT_THROW( + edata.setTimepoints(std::vector{0.0, 1.0, 0.5}), AmiException + ); } -TEST_F(ExpDataTest, SettersGetters) -{ +TEST_F(ExpDataTest, SettersGetters) { ExpData edata(testModel); std::vector y(ny * timepoints.size(), 0.0); @@ -263,22 +247,22 @@ TEST_F(ExpDataTest, SettersGetters) edata.setObservedData(y); checkEqualArray( - edata.getObservedData(), y, TEST_ATOL, TEST_RTOL, "ObservedData"); + edata.getObservedData(), y, TEST_ATOL, TEST_RTOL, "ObservedData" + ); edata.setObservedDataStdDev(y_std); - checkEqualArray(edata.getObservedDataStdDev(), - y_std, - TEST_ATOL, - TEST_RTOL, - "ObservedDataStdDev"); + checkEqualArray( + edata.getObservedDataStdDev(), y_std, TEST_ATOL, TEST_RTOL, + "ObservedDataStdDev" + ); edata.setObservedEvents(z); checkEqualArray( - edata.getObservedEvents(), z, TEST_ATOL, TEST_RTOL, "ObservedEvents"); + edata.getObservedEvents(), z, TEST_ATOL, TEST_RTOL, "ObservedEvents" + ); edata.setObservedEventsStdDev(z_std); - checkEqualArray(edata.getObservedEventsStdDev(), - z_std, - TEST_ATOL, - TEST_RTOL, - "ObservedEventsStdDev"); + checkEqualArray( + edata.getObservedEventsStdDev(), z_std, TEST_ATOL, TEST_RTOL, + "ObservedEventsStdDev" + ); std::vector single_y(edata.nt(), 0.0); std::vector single_y_std(edata.nt(), 0.1); @@ -302,10 +286,12 @@ TEST_F(ExpDataTest, SettersGetters) ASSERT_THROW(edata.setObservedEvents(single_z, nz), std::exception); ASSERT_THROW(edata.setObservedEvents(single_z, -1), std::exception); - ASSERT_THROW(edata.setObservedEventsStdDev(single_z_std, nz), - std::exception); - ASSERT_THROW(edata.setObservedEventsStdDev(single_z_std, -1), - std::exception); + ASSERT_THROW( + edata.setObservedEventsStdDev(single_z_std, nz), std::exception + ); + ASSERT_THROW( + edata.setObservedEventsStdDev(single_z_std, -1), std::exception + ); ASSERT_TRUE(edata.getObservedDataPtr(0)); ASSERT_TRUE(edata.getObservedDataStdDevPtr(0)); @@ -325,23 +311,22 @@ TEST_F(ExpDataTest, SettersGetters) ASSERT_TRUE(!edata.getObservedEventsStdDevPtr(0)); checkEqualArray( - edata.getObservedData(), empty, TEST_ATOL, TEST_RTOL, "ObservedData"); - checkEqualArray(edata.getObservedDataStdDev(), - empty, - TEST_ATOL, - TEST_RTOL, - "ObservedDataStdDev"); + edata.getObservedData(), empty, TEST_ATOL, TEST_RTOL, "ObservedData" + ); + checkEqualArray( + edata.getObservedDataStdDev(), empty, TEST_ATOL, TEST_RTOL, + "ObservedDataStdDev" + ); + checkEqualArray( + edata.getObservedEvents(), empty, TEST_ATOL, TEST_RTOL, "ObservedEvents" + ); checkEqualArray( - edata.getObservedEvents(), empty, TEST_ATOL, TEST_RTOL, "ObservedEvents"); - checkEqualArray(edata.getObservedEventsStdDev(), - empty, - TEST_ATOL, - TEST_RTOL, - "ObservedEventsStdDev"); + edata.getObservedEventsStdDev(), empty, TEST_ATOL, TEST_RTOL, + "ObservedEventsStdDev" + ); } -TEST_F(ExpDataTest, RngSeed) -{ +TEST_F(ExpDataTest, RngSeed) { ReturnData rdata(CVodeSolver(), testModel); rdata.y.assign(testModel.ny * testModel.nt(), 1.0); rdata.z.assign(testModel.nz * testModel.nMaxEvent(), 1.0); diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp index 78bf872cf7..a6f18035a1 100644 --- a/tests/cpp/unittests/testMisc.cpp +++ b/tests/cpp/unittests/testMisc.cpp @@ -2,11 +2,11 @@ #include #include +#include #include #include #include #include -#include #include #include @@ -17,10 +17,7 @@ namespace amici { namespace generic_model { -std::unique_ptr getModel() -{ - return std::make_unique(); -} +std::unique_ptr getModel() { return std::make_unique(); } } // namespace generic_model } // namespace amici @@ -29,114 +26,97 @@ using namespace amici; namespace { -void -testSolverGetterSetters(CVodeSolver solver, - SensitivityMethod sensi_meth, - SensitivityOrder sensi, - InternalSensitivityMethod ism, - InterpolationType interp, - NonlinearSolverIteration iter, - LinearMultistepMethod lmm, - int steps, - int badsteps, - double tol, - double badtol); +void testSolverGetterSetters( + CVodeSolver solver, SensitivityMethod sensi_meth, SensitivityOrder sensi, + InternalSensitivityMethod ism, InterpolationType interp, + NonlinearSolverIteration iter, LinearMultistepMethod lmm, int steps, + int badsteps, double tol, double badtol +); class ModelTest : public ::testing::Test { protected: - int nx = 1, ny = 2, nz = 3, nmaxevent = 4; - std::vector p{ 1.0 }; - std::vector k{ 0.5, 0.4, 0.7 }; - std::vector plist{ 1 }; - std::vector idlist{ 0 }; - std::vector z2event{ 0, 0, 0 }; + int nx = 1, ny = 2, nz = 3, nmaxevent = 0; + std::vector p{1.0}; + std::vector k{0.5, 0.4, 0.7}; + std::vector plist{1}; + std::vector idlist{0}; + std::vector z2event{0, 0, 0}; Model_Test model = Model_Test( ModelDimensions( - nx, // nx_rdata - nx, // nxtrue_rdata - nx, // nx_solver - nx, // nxtrue_solver - 0, // nx_solver_reinit - static_cast(p.size()), // np - static_cast(k.size()), // nk - ny, // ny - ny, // nytrue - nz, // nz - nz, // nztrue - nmaxevent, // ne - 0, // ne_solver - 0, // nspl - 0, // nJ - 0, // nw - 0, // ndwdx - 0, // ndwdp - 0, // dwdw - 0, // ndxdotdw - {}, // ndJydy - 0, // ndxrdatadxsolver - 0, // ndxrdatadtcl - 0, // ndtotal_cldx_rdata - 0, // nnz - 0, // ubw - 0 // lbw - ), - SimulationParameters(k, p, plist), - SecondOrderMode::none, - idlist, - z2event); - std::vector unscaled{ NAN }; + nx, // nx_rdata + nx, // nxtrue_rdata + nx, // nx_solver + nx, // nxtrue_solver + 0, // nx_solver_reinit + static_cast(p.size()), // np + static_cast(k.size()), // nk + ny, // ny + ny, // nytrue + nz, // nz + nz, // nztrue + nmaxevent, // ne + 0, // ne_solver + 0, // nspl + 0, // nJ + 0, // nw + 0, // ndwdx + 0, // ndwdp + 0, // dwdw + 0, // ndxdotdw + {0, 0}, // ndJydy + 0, // ndxrdatadxsolver + 0, // ndxrdatadtcl + 0, // ndtotal_cldx_rdata + 0, // nnz + 0, // ubw + 0 // lbw + ), + SimulationParameters(k, p, plist), SecondOrderMode::none, idlist, + z2event, {} + ); + std::vector unscaled{NAN}; }; - -TEST_F(ModelTest, LinScaledParameterIsNotTransformed) -{ +TEST_F(ModelTest, LinScaledParameterIsNotTransformed) { model.setParameterScale(ParameterScaling::none); ASSERT_EQ(p[0], model.getParameters()[0]); } -TEST_F(ModelTest, LogScaledParameterIsTransformed) -{ +TEST_F(ModelTest, LogScaledParameterIsTransformed) { model.setParameterScale(ParameterScaling::ln); ASSERT_NEAR(std::log(p[0]), model.getParameters()[0], 1e-16); } -TEST_F(ModelTest, Log10ScaledParameterIsTransformed) -{ +TEST_F(ModelTest, Log10ScaledParameterIsTransformed) { model.setParameterScale(ParameterScaling::log10); ASSERT_NEAR(std::log10(p[0]), model.getParameters()[0], 1e-16); } -TEST_F(ModelTest, ParameterScaleTooShort) -{ - std::vector pscale(p.size() - 1, - ParameterScaling::log10); +TEST_F(ModelTest, ParameterScaleTooShort) { + std::vector pscale(p.size() - 1, ParameterScaling::log10); ASSERT_THROW(model.setParameterScale(pscale), AmiException); - } -TEST_F(ModelTest, ParameterScaleTooLong) -{ - std::vector pscale (p.size() + 1, - ParameterScaling::log10); +TEST_F(ModelTest, ParameterScaleTooLong) { + std::vector pscale(p.size() + 1, ParameterScaling::log10); ASSERT_THROW(model.setParameterScale(pscale), AmiException); } -TEST_F(ModelTest, UnsortedTimepointsThrow){ - ASSERT_THROW(model.setTimepoints(std::vector{ 0.0, 1.0, 0.5 }), - AmiException); +TEST_F(ModelTest, UnsortedTimepointsThrow) { + ASSERT_THROW( + model.setTimepoints(std::vector{0.0, 1.0, 0.5}), AmiException + ); } -TEST_F(ModelTest, ParameterNameIdGetterSetter) -{ +TEST_F(ModelTest, ParameterNameIdGetterSetter) { model.setParameterById("p0", 3.0); ASSERT_NEAR(model.getParameterById("p0"), 3.0, 1e-16); ASSERT_THROW(model.getParameterById("p1"), AmiException); - ASSERT_NEAR( - model.setParametersByIdRegex("p[\\d]+", 5.0), p.size(), 1e-16); - for (const auto& ip : model.getParameters()) + ASSERT_NEAR(model.setParametersByIdRegex("p[\\d]+", 5.0), p.size(), 1e-16); + for (auto const& ip : model.getParameters()) ASSERT_NEAR(ip, 5.0, 1e-16); ASSERT_THROW(model.setParametersByIdRegex("k[\\d]+", 5.0), AmiException); @@ -144,8 +124,9 @@ TEST_F(ModelTest, ParameterNameIdGetterSetter) ASSERT_NEAR(model.getParameterByName("p0"), 3.0, 1e-16); ASSERT_THROW(model.getParameterByName("p1"), AmiException); ASSERT_NEAR( - model.setParametersByNameRegex("p[\\d]+", 5.0), p.size(), 1e-16); - for (const auto& ip : model.getParameters()) + model.setParametersByNameRegex("p[\\d]+", 5.0), p.size(), 1e-16 + ); + for (auto const& ip : model.getParameters()) ASSERT_NEAR(ip, 5.0, 1e-16); ASSERT_THROW(model.setParametersByNameRegex("k[\\d]+", 5.0), AmiException); @@ -153,26 +134,31 @@ TEST_F(ModelTest, ParameterNameIdGetterSetter) ASSERT_NEAR(model.getFixedParameterById("k0"), 3.0, 1e-16); ASSERT_THROW(model.getFixedParameterById("k4"), AmiException); ASSERT_NEAR( - model.setFixedParametersByIdRegex("k[\\d]+", 5.0), k.size(), 1e-16); - for (const auto& ik : model.getFixedParameters()) + model.setFixedParametersByIdRegex("k[\\d]+", 5.0), k.size(), 1e-16 + ); + for (auto const& ik : model.getFixedParameters()) ASSERT_NEAR(ik, 5.0, 1e-16); - ASSERT_THROW(model.setFixedParametersByIdRegex("p[\\d]+", 5.0), AmiException); + ASSERT_THROW( + model.setFixedParametersByIdRegex("p[\\d]+", 5.0), AmiException + ); model.setFixedParameterByName("k0", 3.0); ASSERT_NEAR(model.getFixedParameterByName("k0"), 3.0, 1e-16); ASSERT_THROW(model.getFixedParameterByName("k4"), AmiException); ASSERT_NEAR( - model.setFixedParametersByNameRegex("k[\\d]+", 5.0), k.size(), 1e-16); - for (const auto& ik : model.getFixedParameters()) + model.setFixedParametersByNameRegex("k[\\d]+", 5.0), k.size(), 1e-16 + ); + for (auto const& ik : model.getFixedParameters()) ASSERT_NEAR(ik, 5.0, 1e-16); - ASSERT_THROW(model.setFixedParametersByNameRegex("p[\\d]+", 5.0), - AmiException); + ASSERT_THROW( + model.setFixedParametersByNameRegex("p[\\d]+", 5.0), AmiException + ); } -TEST_F(ModelTest, ReinitializeFixedParameterInitialStates) -{ - ASSERT_THROW(model.setReinitializeFixedParameterInitialStates(true), - AmiException); +TEST_F(ModelTest, ReinitializeFixedParameterInitialStates) { + ASSERT_THROW( + model.setReinitializeFixedParameterInitialStates(true), AmiException + ); model.setReinitializeFixedParameterInitialStates(false); ASSERT_TRUE(!model.getReinitializeFixedParameterInitialStates()); sundials::Context sunctx; @@ -180,22 +166,19 @@ TEST_F(ModelTest, ReinitializeFixedParameterInitialStates) AmiVectorArray sx(model.np(), nx, sunctx); } -TEST(SymbolicFunctionsTest, Sign) -{ +TEST(SymbolicFunctionsTest, Sign) { ASSERT_EQ(-1, sign(-2)); ASSERT_EQ(0, sign(0)); ASSERT_EQ(1, sign(2)); } -TEST(SymbolicFunctionsTest, Heaviside) -{ +TEST(SymbolicFunctionsTest, Heaviside) { ASSERT_EQ(0, heaviside(-1, 0.5)); ASSERT_EQ(0.5, heaviside(0, 0.5)); ASSERT_EQ(1, heaviside(1, 0.5)); } -TEST(SymbolicFunctionsTest, Min) -{ +TEST(SymbolicFunctionsTest, Min) { ASSERT_EQ(-1, min(-1, 2, 0)); ASSERT_EQ(-2, min(1, -2, 0)); ASSERT_TRUE(isNaN(min(getNaN(), getNaN(), 0))); @@ -203,8 +186,7 @@ TEST(SymbolicFunctionsTest, Min) ASSERT_EQ(-1, min(getNaN(), -1, 0)); } -TEST(SymbolicFunctionsTest, Max) -{ +TEST(SymbolicFunctionsTest, Max) { ASSERT_EQ(2, max(-1, 2, 0)); ASSERT_EQ(1, max(1, -2, 0)); ASSERT_TRUE(isNaN(max(getNaN(), getNaN(), 0))); @@ -212,30 +194,26 @@ TEST(SymbolicFunctionsTest, Max) ASSERT_EQ(-1, max(getNaN(), -1, 0)); } -TEST(SymbolicFunctionsTest, DMin) -{ +TEST(SymbolicFunctionsTest, DMin) { ASSERT_EQ(0, Dmin(1, -1, -2, 0)); ASSERT_EQ(1, Dmin(1, -1, 2, 0)); ASSERT_EQ(1, Dmin(2, -1, -2, 0)); ASSERT_EQ(0, Dmin(2, -1, 2, 0)); } -TEST(SymbolicFunctionsTest, DMax) -{ +TEST(SymbolicFunctionsTest, DMax) { ASSERT_EQ(1, Dmax(1, -1, -2, 0)); ASSERT_EQ(0, Dmax(1, -1, 2, 0)); ASSERT_EQ(0, Dmax(2, -1, -2, 0)); ASSERT_EQ(1, Dmax(2, -1, 2, 0)); } -TEST(SymbolicFunctionsTest, pos_pow) -{ +TEST(SymbolicFunctionsTest, pos_pow) { ASSERT_EQ(0, pos_pow(-0.1, 3)); ASSERT_EQ(pow(0.1, 3), pos_pow(0.1, 3)); } -TEST(SolverTestBasic, Equality) -{ +TEST(SolverTestBasic, Equality) { IDASolver i1, i2; CVodeSolver c1, c2; @@ -244,8 +222,7 @@ TEST(SolverTestBasic, Equality) ASSERT_FALSE(i1 == c1); } -TEST(SolverTestBasic, Clone) -{ +TEST(SolverTestBasic, Clone) { IDASolver i1; std::unique_ptr i2(i1.clone()); ASSERT_EQ(i1, *i2); @@ -256,19 +233,13 @@ TEST(SolverTestBasic, Clone) ASSERT_FALSE(*i2 == *c2); } -TEST(SolverIdasTest, DefaultConstructableAndNotLeaky) -{ - IDASolver solver; -} +TEST(SolverIdasTest, DefaultConstructableAndNotLeaky) { IDASolver solver; } -TEST(SolverIdasTest, CopyCtor) -{ +TEST(SolverIdasTest, CopyCtor) { IDASolver solver1; IDASolver solver2(solver1); } - - class SolverTest : public ::testing::Test { protected: void SetUp() override { @@ -286,7 +257,7 @@ class SolverTest : public ::testing::Test { int nx = 1, ny = 2, nz = 3, ne = 0; double tol, badtol; - std::vector timepoints = { 1, 2, 3, 4 }; + std::vector timepoints = {1, 2, 3, 4}; std::unique_ptr model = generic_model::getModel(); SensitivityMethod sensi_meth; @@ -299,127 +270,117 @@ class SolverTest : public ::testing::Test { Model_Test testModel = Model_Test( ModelDimensions( - nx, // nx_rdata - nx, // nxtrue_rdata - nx, // nx_solver - nx, // nxtrue_solver - 0, // nx_solver_reinit - 1, // np - 3, // nk - ny, // ny - ny, // nytrue - nz, // nz - nz, // nztrue - ne, // ne - 0, // ne_solver - 0, // nspl - 0, // nJ - 0, // nw - 0, // ndwdx - 0, // ndwdp - 0, // dwdw - 0, // ndxdotdw - {}, // ndJydy - 0, // ndxrdatadxsolver - 0, // ndxrdatadtcl - 0, // ndtotal_cldx_rdata - 1, // nnz - 0, // ubw - 0 // lbw - ), + nx, // nx_rdata + nx, // nxtrue_rdata + nx, // nx_solver + nx, // nxtrue_solver + 0, // nx_solver_reinit + 1, // np + 3, // nk + ny, // ny + ny, // nytrue + nz, // nz + nz, // nztrue + ne, // ne + 0, // ne_solver + 0, // nspl + 0, // nJ + 0, // nw + 0, // ndwdx + 0, // ndwdp + 0, // dwdw + 0, // ndxdotdw + {0, 0}, // ndJydy + 0, // ndxrdatadxsolver + 0, // ndxrdatadtcl + 0, // ndtotal_cldx_rdata + 1, // nnz + 0, // ubw + 0 // lbw + ), SimulationParameters( - std::vector(3, 0.0), - std::vector(1, 0.0), + std::vector(3, 0.0), std::vector(1, 0.0), std::vector(2, 1) - ), - SecondOrderMode::none, - std::vector(0, 0.0), - std::vector()); + ), + SecondOrderMode::none, std::vector(0, 0.0), + std::vector(), {} + ); CVodeSolver solver = CVodeSolver(); }; -TEST_F(SolverTest, SettersGettersNoSetup) -{ - testSolverGetterSetters(solver, - sensi_meth, - sensi, - ism, - interp, - iter, - lmm, - steps, - badsteps, - tol, - badtol); -} - -TEST_F(SolverTest, SettersGettersWithSetup) -{ +TEST_F(SolverTest, SettersGettersNoSetup) { + testSolverGetterSetters( + solver, sensi_meth, sensi, ism, interp, iter, lmm, steps, badsteps, tol, + badtol + ); +} + +TEST_F(SolverTest, SettersGettersWithSetup) { solver.setSensitivityMethod(sensi_meth); - ASSERT_EQ(static_cast(solver.getSensitivityMethod()), - static_cast(sensi_meth)); + ASSERT_EQ( + static_cast(solver.getSensitivityMethod()), + static_cast(sensi_meth) + ); auto rdata = std::make_unique(solver, testModel); AmiVector x(nx, solver.getSunContext()), dx(nx, solver.getSunContext()); - AmiVectorArray sx(nx, 1, solver.getSunContext()), sdx(nx, 1, solver.getSunContext()); + AmiVectorArray sx(nx, 1, solver.getSunContext()), + sdx(nx, 1, solver.getSunContext()); - testModel.setInitialStates(std::vector{ 0 }); + testModel.setInitialStates(std::vector{0}); solver.setup(0, &testModel, x, dx, sx, sdx); - testSolverGetterSetters(solver, - sensi_meth, - sensi, - ism, - interp, - iter, - lmm, - steps, - badsteps, - tol, - badtol); -} - -void -testSolverGetterSetters(CVodeSolver solver, - SensitivityMethod sensi_meth, - SensitivityOrder sensi, - InternalSensitivityMethod ism, - InterpolationType interp, - NonlinearSolverIteration iter, - LinearMultistepMethod lmm, - int steps, - int badsteps, - double tol, - double badtol) -{ + testSolverGetterSetters( + solver, sensi_meth, sensi, ism, interp, iter, lmm, steps, badsteps, tol, + badtol + ); +} + +void testSolverGetterSetters( + CVodeSolver solver, SensitivityMethod sensi_meth, SensitivityOrder sensi, + InternalSensitivityMethod ism, InterpolationType interp, + NonlinearSolverIteration iter, LinearMultistepMethod lmm, int steps, + int badsteps, double tol, double badtol +) { solver.setSensitivityMethod(sensi_meth); - ASSERT_EQ(static_cast(solver.getSensitivityMethod()), - static_cast(sensi_meth)); + ASSERT_EQ( + static_cast(solver.getSensitivityMethod()), + static_cast(sensi_meth) + ); solver.setSensitivityOrder(sensi); - ASSERT_EQ(static_cast(solver.getSensitivityOrder()), - static_cast(sensi)); + ASSERT_EQ( + static_cast(solver.getSensitivityOrder()), static_cast(sensi) + ); solver.setInternalSensitivityMethod(ism); - ASSERT_EQ(static_cast(solver.getInternalSensitivityMethod()), - static_cast(ism)); + ASSERT_EQ( + static_cast(solver.getInternalSensitivityMethod()), + static_cast(ism) + ); solver.setInterpolationType(interp); - ASSERT_EQ(static_cast(solver.getInterpolationType()), - static_cast(interp)); + ASSERT_EQ( + static_cast(solver.getInterpolationType()), + static_cast(interp) + ); solver.setNonlinearSolverIteration(iter); - ASSERT_EQ(static_cast(solver.getNonlinearSolverIteration()), - static_cast(iter)); + ASSERT_EQ( + static_cast(solver.getNonlinearSolverIteration()), + static_cast(iter) + ); solver.setLinearMultistepMethod(lmm); - ASSERT_EQ(static_cast(solver.getLinearMultistepMethod()), - static_cast(lmm)); + ASSERT_EQ( + static_cast(solver.getLinearMultistepMethod()), + static_cast(lmm) + ); - solver.setStabilityLimitFlag(true); + solver.setStabilityLimitFlag(true); ASSERT_EQ(solver.getStabilityLimitFlag(), true); ASSERT_THROW(solver.setNewtonMaxSteps(badsteps), AmiException); @@ -459,68 +420,81 @@ testSolverGetterSetters(CVodeSolver solver, ASSERT_EQ(solver.getAbsoluteToleranceSteadyState(), tol); } -TEST_F(SolverTest, SteadyStateToleranceFactor) -{ +TEST_F(SolverTest, SteadyStateToleranceFactor) { CVodeSolver s; // test with unset steadystate tolerances ASSERT_DOUBLE_EQ( s.getRelativeToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getRelativeTolerance()); + s.getSteadyStateToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_DOUBLE_EQ( s.getAbsoluteToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance()); + s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance() + ); ASSERT_DOUBLE_EQ( s.getRelativeToleranceSteadyStateSensi(), - s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance()); + s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_DOUBLE_EQ( s.getAbsoluteToleranceSteadyState(), - s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance()); + s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance() + ); // test with changed steadystate tolerance factor s.setSteadyStateToleranceFactor(5); ASSERT_DOUBLE_EQ( s.getRelativeToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getRelativeTolerance()); + s.getSteadyStateToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_DOUBLE_EQ( s.getAbsoluteToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance()); + s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance() + ); s.setSteadyStateSensiToleranceFactor(5); ASSERT_DOUBLE_EQ( s.getRelativeToleranceSteadyStateSensi(), - s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance()); + s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_DOUBLE_EQ( s.getAbsoluteToleranceSteadyState(), - s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance()); - + s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance() + ); // test with steadystate tolerance override tolerance factor s.setRelativeToleranceSteadyState(2); - ASSERT_NE(s.getRelativeToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getRelativeTolerance()); + ASSERT_NE( + s.getRelativeToleranceSteadyState(), + s.getSteadyStateToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_EQ(s.getRelativeToleranceSteadyState(), 2); s.setAbsoluteToleranceSteadyState(3); - ASSERT_NE(s.getAbsoluteToleranceSteadyState(), - s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance()); + ASSERT_NE( + s.getAbsoluteToleranceSteadyState(), + s.getSteadyStateToleranceFactor() * s.getAbsoluteTolerance() + ); ASSERT_EQ(s.getAbsoluteToleranceSteadyState(), 3); s.setRelativeToleranceSteadyStateSensi(4); - ASSERT_NE(s.getRelativeToleranceSteadyStateSensi(), - s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance()); + ASSERT_NE( + s.getRelativeToleranceSteadyStateSensi(), + s.getSteadyStateSensiToleranceFactor() * s.getRelativeTolerance() + ); ASSERT_EQ(s.getRelativeToleranceSteadyStateSensi(), 4); s.setAbsoluteToleranceSteadyStateSensi(5); - ASSERT_NE(s.getAbsoluteToleranceSteadyStateSensi(), - s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance()); + ASSERT_NE( + s.getAbsoluteToleranceSteadyStateSensi(), + s.getSteadyStateSensiToleranceFactor() * s.getAbsoluteTolerance() + ); ASSERT_EQ(s.getAbsoluteToleranceSteadyStateSensi(), 5); } class AmiVectorTest : public ::testing::Test { protected: - std::vector vec1{ 1, 2, 4, 3 }; - std::vector vec2{ 4, 1, 2, 3 }; - std::vector vec3{ 4, 4, 2, 1 }; + std::vector vec1{1, 2, 4, 3}; + std::vector vec2{4, 1, 2, 3}; + std::vector vec3{4, 4, 2, 1}; }; -TEST_F(AmiVectorTest, Vector) -{ +TEST_F(AmiVectorTest, Vector) { sundials::Context sunctx; AmiVector av(vec1, sunctx); N_Vector nvec = av.getNVector(); @@ -537,12 +511,11 @@ TEST_F(AmiVectorTest, Vector) ASSERT_EQ(ss.str(), "[1, 2, 4, 3]"); } -TEST_F(AmiVectorTest, VectorArray) -{ +TEST_F(AmiVectorTest, VectorArray) { sundials::Context sunctx; AmiVectorArray ava(4, 3, sunctx); AmiVector av1(vec1, sunctx), av2(vec2, sunctx), av3(vec3, sunctx); - std::vector avs{ av1, av2, av3 }; + std::vector avs{av1, av2, av3}; for (int i = 0; i < ava.getLength(); ++i) ava[i] = avs.at(i); @@ -552,7 +525,7 @@ TEST_F(AmiVectorTest, VectorArray) ASSERT_THROW(ava.flatten_to_vector(badLengthVector), AmiException); ava.flatten_to_vector(flattened); for (int i = 0; i < ava.getLength(); ++i) { - const AmiVector av = ava[i]; + AmiVector const av = ava[i]; for (int j = 0; j < av.getLength(); ++j) ASSERT_EQ(flattened.at(i * av.getLength() + j), av.at(j)); } @@ -594,7 +567,7 @@ class SunMatrixWrapperTest : public ::testing::Test { } sundials::Context sunctx; - //inputs + // inputs std::vector a{0.82, 0.91, 0.13}; std::vector b{0.77, 0.80}; SUNMatrixWrapper A = SUNMatrixWrapper(3, 2, sunctx); @@ -603,22 +576,20 @@ class SunMatrixWrapperTest : public ::testing::Test { std::vector d{1.3753, 1.5084, 1.1655}; }; -TEST_F(SunMatrixWrapperTest, SparseMultiply) -{ +TEST_F(SunMatrixWrapperTest, SparseMultiply) { auto A_sparse = SUNMatrixWrapper(A, 0.0, CSC_MAT); - auto c(a); //copy c + auto c(a); // copy c A_sparse.multiply(c, b); checkEqualArray(d, c, TEST_ATOL, TEST_RTOL, "multiply"); } -TEST_F(SunMatrixWrapperTest, SparseMultiplyEmpty) -{ +TEST_F(SunMatrixWrapperTest, SparseMultiplyEmpty) { // Ensure empty Matrix vector multiplication succeeds sundials::Context sunctx; auto A_sparse = SUNMatrixWrapper(1, 1, 0, CSR_MAT, sunctx); - std::vector b {0.1}; - std::vector c {0.1}; + std::vector b{0.1}; + std::vector c{0.1}; A_sparse.multiply(c, b); ASSERT_TRUE(c[0] == 0.1); @@ -627,31 +598,28 @@ TEST_F(SunMatrixWrapperTest, SparseMultiplyEmpty) ASSERT_TRUE(c[0] == 0.1); } -TEST_F(SunMatrixWrapperTest, DenseMultiply) -{ - auto c(a); //copy c +TEST_F(SunMatrixWrapperTest, DenseMultiply) { + auto c(a); // copy c A.multiply(c, b); checkEqualArray(d, c, TEST_ATOL, TEST_RTOL, "multiply"); } -TEST_F(SunMatrixWrapperTest, StdVectorCtor) -{ +TEST_F(SunMatrixWrapperTest, StdVectorCtor) { sundials::Context sunctx; auto b_amivector = AmiVector(b, sunctx); auto a_amivector = AmiVector(a, sunctx); } -TEST_F(SunMatrixWrapperTest, TransformThrows) -{ +TEST_F(SunMatrixWrapperTest, TransformThrows) { sundials::Context sunctx; ASSERT_THROW(SUNMatrixWrapper(A, 0.0, 13), std::invalid_argument); auto A_sparse = SUNMatrixWrapper(A, 0.0, CSR_MAT); - ASSERT_THROW(SUNMatrixWrapper(A_sparse, 0.0, CSR_MAT), - std::invalid_argument); + ASSERT_THROW( + SUNMatrixWrapper(A_sparse, 0.0, CSR_MAT), std::invalid_argument + ); } -TEST_F(SunMatrixWrapperTest, BlockTranspose) -{ +TEST_F(SunMatrixWrapperTest, BlockTranspose) { sundials::Context sunctx; SUNMatrixWrapper B_sparse(4, 4, 7, CSR_MAT, sunctx); ASSERT_THROW(B.transpose(B_sparse, 1.0, 4), std::domain_error); @@ -659,47 +627,45 @@ TEST_F(SunMatrixWrapperTest, BlockTranspose) B_sparse = SUNMatrixWrapper(4, 4, 7, CSC_MAT, sunctx); B.transpose(B_sparse, -1.0, 2); for (int idx = 0; idx < 7; idx++) { - ASSERT_EQ(SM_INDEXVALS_S(B.get())[idx], - SM_INDEXVALS_S(B_sparse.get())[idx]); + ASSERT_EQ( + SM_INDEXVALS_S(B.get())[idx], SM_INDEXVALS_S(B_sparse.get())[idx] + ); if (idx == 1) { - ASSERT_EQ(SM_DATA_S(B.get())[idx], - -SM_DATA_S(B_sparse.get())[3]); + ASSERT_EQ(SM_DATA_S(B.get())[idx], -SM_DATA_S(B_sparse.get())[3]); } else if (idx == 3) { - ASSERT_EQ(SM_DATA_S(B.get())[idx], - -SM_DATA_S(B_sparse.get())[1]); + ASSERT_EQ(SM_DATA_S(B.get())[idx], -SM_DATA_S(B_sparse.get())[1]); } else { - ASSERT_EQ(SM_DATA_S(B.get())[idx], - -SM_DATA_S(B_sparse.get())[idx]); + ASSERT_EQ(SM_DATA_S(B.get())[idx], -SM_DATA_S(B_sparse.get())[idx]); } } for (int icol = 0; icol <= 4; icol++) - ASSERT_EQ(SM_INDEXPTRS_S(B.get())[icol], - SM_INDEXPTRS_S(B_sparse.get())[icol]); + ASSERT_EQ( + SM_INDEXPTRS_S(B.get())[icol], SM_INDEXPTRS_S(B_sparse.get())[icol] + ); } -TEST_F(SunMatrixWrapperTest, TestPrint) -{ +TEST_F(SunMatrixWrapperTest, TestPrint) { std::stringstream ss; - ss << A; - EXPECT_EQ(ss.str() , "[[0.69, 0.03], [0.32, 0.44], [0.95, 0.38]]"); - ss.str(""); - ss<< B; - EXPECT_EQ(ss.str(), "[[0, 3, 1, 0], [3, 0, 0, 2], [0, 7, 0, 0], [1, 0, 0, 9]]"); -} - -TEST(UnravelIndex, UnravelIndex) -{ - EXPECT_EQ(unravel_index(0, 2), std::make_pair((size_t) 0, (size_t) 0)); - EXPECT_EQ(unravel_index(1, 2), std::make_pair((size_t) 0, (size_t) 1)); - EXPECT_EQ(unravel_index(2, 2), std::make_pair((size_t) 1, (size_t) 0)); - EXPECT_EQ(unravel_index(3, 2), std::make_pair((size_t) 1, (size_t) 1)); - EXPECT_EQ(unravel_index(4, 2), std::make_pair((size_t) 2, (size_t) 0)); - EXPECT_EQ(unravel_index(5, 2), std::make_pair((size_t) 2, (size_t) 1)); - EXPECT_EQ(unravel_index(6, 2), std::make_pair((size_t) 3, (size_t) 0)); -} - -TEST(UnravelIndex, UnravelIndexSunMatDense) -{ + ss << A; + EXPECT_EQ(ss.str(), "[[0.69, 0.03], [0.32, 0.44], [0.95, 0.38]]"); + ss.str(""); + ss << B; + EXPECT_EQ( + ss.str(), "[[0, 3, 1, 0], [3, 0, 0, 2], [0, 7, 0, 0], [1, 0, 0, 9]]" + ); +} + +TEST(UnravelIndex, UnravelIndex) { + EXPECT_EQ(unravel_index(0, 2), std::make_pair((size_t)0, (size_t)0)); + EXPECT_EQ(unravel_index(1, 2), std::make_pair((size_t)0, (size_t)1)); + EXPECT_EQ(unravel_index(2, 2), std::make_pair((size_t)1, (size_t)0)); + EXPECT_EQ(unravel_index(3, 2), std::make_pair((size_t)1, (size_t)1)); + EXPECT_EQ(unravel_index(4, 2), std::make_pair((size_t)2, (size_t)0)); + EXPECT_EQ(unravel_index(5, 2), std::make_pair((size_t)2, (size_t)1)); + EXPECT_EQ(unravel_index(6, 2), std::make_pair((size_t)3, (size_t)0)); +} + +TEST(UnravelIndex, UnravelIndexSunMatDense) { sundials::Context sunctx; SUNMatrixWrapper A = SUNMatrixWrapper(3, 2, sunctx); @@ -710,14 +676,13 @@ TEST(UnravelIndex, UnravelIndexSunMatDense) A.set_data(1, 1, 4); A.set_data(2, 1, 5); - for(int i = 0; i < 6; ++i) { + for (int i = 0; i < 6; ++i) { auto idx = unravel_index(i, A); EXPECT_EQ(A.get_data(idx.first, idx.second), i); } } -TEST(UnravelIndex, UnravelIndexSunMatSparse) -{ +TEST(UnravelIndex, UnravelIndexSunMatSparse) { sundials::Context sunctx; SUNMatrixWrapper D = SUNMatrixWrapper(4, 2, sunctx); @@ -739,43 +704,51 @@ TEST(UnravelIndex, UnravelIndexSunMatSparse) auto S = SUNSparseFromDenseMatrix(D, 1e-15, CSC_MAT); - EXPECT_EQ(unravel_index(0, S), std::make_pair((sunindextype) 2, (sunindextype) 0)); - EXPECT_EQ(unravel_index(1, S), std::make_pair((sunindextype) 3, (sunindextype) 0)); - EXPECT_EQ(unravel_index(2, S), std::make_pair((sunindextype) 0, (sunindextype) 1)); + EXPECT_EQ( + unravel_index(0, S), std::make_pair((sunindextype)2, (sunindextype)0) + ); + EXPECT_EQ( + unravel_index(1, S), std::make_pair((sunindextype)3, (sunindextype)0) + ); + EXPECT_EQ( + unravel_index(2, S), std::make_pair((sunindextype)0, (sunindextype)1) + ); SUNMatDestroy(S); } - -TEST(UnravelIndex, UnravelIndexSunMatSparseMissingIndices) -{ +TEST(UnravelIndex, UnravelIndexSunMatSparseMissingIndices) { sundials::Context sunctx; // Sparse matrix without any indices set SUNMatrixWrapper mat = SUNMatrixWrapper(2, 3, 2, CSC_MAT, sunctx); - EXPECT_EQ(unravel_index(0, mat), std::make_pair((sunindextype) -1, (sunindextype) -1)); - EXPECT_EQ(unravel_index(1, mat), std::make_pair((sunindextype) -1, (sunindextype) -1)); + EXPECT_EQ( + unravel_index(0, mat), + std::make_pair((sunindextype)-1, (sunindextype)-1) + ); + EXPECT_EQ( + unravel_index(1, mat), + std::make_pair((sunindextype)-1, (sunindextype)-1) + ); } - -TEST(ReturnCodeToStr, ReturnCodeToStr) -{ +TEST(ReturnCodeToStr, ReturnCodeToStr) { EXPECT_EQ("AMICI_SUCCESS", simulation_status_to_str(AMICI_SUCCESS)); - EXPECT_EQ("AMICI_UNRECOVERABLE_ERROR", - simulation_status_to_str(AMICI_UNRECOVERABLE_ERROR)); + EXPECT_EQ( + "AMICI_UNRECOVERABLE_ERROR", + simulation_status_to_str(AMICI_UNRECOVERABLE_ERROR) + ); } -TEST(SpanEqual, SpanEqual) -{ - std::vector a {1, 2, 3}; - std::vector b {1, 2, NAN}; +TEST(SpanEqual, SpanEqual) { + std::vector a{1, 2, 3}; + std::vector b{1, 2, NAN}; EXPECT_TRUE(is_equal(a, a)); EXPECT_TRUE(is_equal(b, b)); EXPECT_FALSE(is_equal(a, b)); } -TEST(CpuTimer, CpuTimer) -{ +TEST(CpuTimer, CpuTimer) { amici::CpuTimer timer; auto elapsed = timer.elapsed_seconds(); EXPECT_LE(0.0, elapsed); diff --git a/tests/cpp/unittests/testSerialization.cpp b/tests/cpp/unittests/testSerialization.cpp index e46f17202c..48e77cf625 100644 --- a/tests/cpp/unittests/testSerialization.cpp +++ b/tests/cpp/unittests/testSerialization.cpp @@ -1,16 +1,16 @@ +#include "testfunctions.h" +#include #include #include #include -#include "testfunctions.h" - #include #include #include -void -checkReturnDataEqual(amici::ReturnData const& r, amici::ReturnData const& s) -{ +void checkReturnDataEqual( + amici::ReturnData const& r, amici::ReturnData const& s +) { ASSERT_EQ(r.id, s.id); ASSERT_EQ(r.np, s.np); ASSERT_EQ(r.nk, s.nk); @@ -65,18 +65,26 @@ checkReturnDataEqual(amici::ReturnData const& r, amici::ReturnData const& s) ASSERT_EQ(r.cpu_timeB, s.cpu_timeB); ASSERT_EQ(r.preeq_status, s.preeq_status); - ASSERT_TRUE(r.preeq_t == s.preeq_t || - (std::isnan(r.preeq_t) && std::isnan(s.preeq_t))); - ASSERT_TRUE(r.preeq_wrms == s.preeq_wrms || - (std::isnan(r.preeq_wrms) && std::isnan(s.preeq_wrms))); + ASSERT_TRUE( + r.preeq_t == s.preeq_t + || (std::isnan(r.preeq_t) && std::isnan(s.preeq_t)) + ); + ASSERT_TRUE( + r.preeq_wrms == s.preeq_wrms + || (std::isnan(r.preeq_wrms) && std::isnan(s.preeq_wrms)) + ); ASSERT_EQ(r.preeq_numsteps, s.preeq_numsteps); EXPECT_NEAR(r.preeq_cpu_time, s.preeq_cpu_time, 1e-16); ASSERT_EQ(r.posteq_status, s.posteq_status); - ASSERT_TRUE(r.posteq_t == s.posteq_t || - (std::isnan(r.posteq_t) && std::isnan(s.posteq_t))); - ASSERT_TRUE(r.posteq_wrms == s.posteq_wrms || - (std::isnan(r.posteq_wrms) && std::isnan(s.posteq_wrms))); + ASSERT_TRUE( + r.posteq_t == s.posteq_t + || (std::isnan(r.posteq_t) && std::isnan(s.posteq_t)) + ); + ASSERT_TRUE( + r.posteq_wrms == s.posteq_wrms + || (std::isnan(r.posteq_wrms) && std::isnan(s.posteq_wrms)) + ); ASSERT_EQ(r.posteq_numsteps, s.posteq_numsteps); EXPECT_NEAR(r.posteq_cpu_time, s.posteq_cpu_time, 1e-16); @@ -106,21 +114,26 @@ class SolverSerializationTest : public ::testing::Test { solver.setMaxSteps(1e1); solver.setMaxStepsBackwardProblem(1e2); solver.setNewtonMaxSteps(1e3); - solver.setStateOrdering(static_cast(amici::SUNLinSolKLU::StateOrdering::COLAMD)); + solver.setStateOrdering( + static_cast(amici::SUNLinSolKLU::StateOrdering::COLAMD) + ); solver.setInterpolationType(amici::InterpolationType::polynomial); solver.setStabilityLimitFlag(false); solver.setLinearSolver(amici::LinearSolver::dense); solver.setLinearMultistepMethod(amici::LinearMultistepMethod::adams); - solver.setNonlinearSolverIteration(amici::NonlinearSolverIteration::newton); - solver.setInternalSensitivityMethod(amici::InternalSensitivityMethod::staggered); + solver.setNonlinearSolverIteration( + amici::NonlinearSolverIteration::newton + ); + solver.setInternalSensitivityMethod( + amici::InternalSensitivityMethod::staggered + ); solver.setReturnDataReportingMode(amici::RDataReporting::likelihood); } amici::CVodeSolver solver; }; -TEST(ModelSerializationTest, ToFile) -{ +TEST(ModelSerializationTest, ToFile) { using amici::realtype; int np = 1; int nk = 2; @@ -131,42 +144,49 @@ TEST(ModelSerializationTest, ToFile) amici::CVodeSolver solver; amici::Model_Test m = amici::Model_Test( amici::ModelDimensions( - nx, // nx_rdata - nx, // nxtrue_rdata - nx, // nx_solver - nx, // nxtrue_solver - 0, // nx_solver_reinit - np, // np - nk, // nk - ny, // ny - ny, // nytrue - nz, // nz - nz, // nztrue - ne, // ne - 0, // ne_solver - 0, // nspl - 0, // nJ - 9, // nw - 2, // ndwdx - 2, // ndwdp - 2, // dwdw - 13, // ndxdotdw - {}, // ndJydy - 9, // ndxrdatadxsolver - 0, // ndxrdatadtcl - 0, // ndtotal_cldx_rdata - 17, // nnz - 18, // ubw - 19 // lbw - ), + nx, // nx_rdata + nx, // nxtrue_rdata + nx, // nx_solver + nx, // nxtrue_solver + 0, // nx_solver_reinit + np, // np + nk, // nk + ny, // ny + ny, // nytrue + nz, // nz + nz, // nztrue + ne, // ne + 0, // ne_solver + 0, // nspl + 0, // nJ + 9, // nw + 2, // ndwdx + 2, // ndwdp + 2, // dwdw + 13, // ndxdotdw + {0, 0, 0, 0}, // ndJydy + 9, // ndxrdatadxsolver + 0, // ndxrdatadtcl + 0, // ndtotal_cldx_rdata + 17, // nnz + 18, // ubw + 19 // lbw + ), amici::SimulationParameters( - std::vector(nk, 0.0), - std::vector(np, 0.0), + std::vector(nk, 0.0), std::vector(np, 0.0), std::vector(np, 0) - ), - amici::SecondOrderMode::none, - std::vector(nx, 0.0), - std::vector(nz, 0)); + ), + amici::SecondOrderMode::none, std::vector(nx, 0.0), + std::vector(nz, 0), + { + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + } + ); { std::ofstream ofs("sstore.dat"); @@ -186,8 +206,7 @@ TEST(ModelSerializationTest, ToFile) } } -TEST(ReturnDataSerializationTest, ToString) -{ +TEST(ReturnDataSerializationTest, ToString) { using amici::realtype; int np = 1; int nk = 2; @@ -198,69 +217,76 @@ TEST(ReturnDataSerializationTest, ToString) amici::CVodeSolver solver; amici::Model_Test m = amici::Model_Test( amici::ModelDimensions( - nx, // nx_rdata - nx, // nxtrue_rdata - nx, // nx_solver - nx, // nxtrue_solver - 0, // nx_solver_reinit - np, // np - nk, // nk - ny, // ny - ny, // nytrue - nz, // nz - nz, // nztrue - ne, // ne - 0, // ne_solver - 0, // nspl - 0, // nJ - 9, // nw - 10, // ndwdx - 2, // ndwdp - 12, // dwdw - 13, // ndxdotdw - {}, // ndJydy - 9, // ndxrdatadxsolver - 0, // ndxrdatadtcl - 0, // ndtotal_cldx_rdata - 17, // nnz - 18, // ubw - 19 // lbw - ), + nx, // nx_rdata + nx, // nxtrue_rdata + nx, // nx_solver + nx, // nxtrue_solver + 0, // nx_solver_reinit + np, // np + nk, // nk + ny, // ny + ny, // nytrue + nz, // nz + nz, // nztrue + ne, // ne + 0, // ne_solver + 0, // nspl + 0, // nJ + 9, // nw + 10, // ndwdx + 2, // ndwdp + 12, // dwdw + 13, // ndxdotdw + {0, 0, 0, 0}, // ndJydy + 9, // ndxrdatadxsolver + 0, // ndxrdatadtcl + 0, // ndtotal_cldx_rdata + 17, // nnz + 18, // ubw + 19 // lbw + ), amici::SimulationParameters( - std::vector(nk, 0.0), - std::vector(np, 0.0), + std::vector(nk, 0.0), std::vector(np, 0.0), std::vector(np, 0) - ), - amici::SecondOrderMode::none, - std::vector(nx, 0.0), - std::vector(nz, 0)); + ), + amici::SecondOrderMode::none, std::vector(nx, 0.0), + std::vector(nz, 0), + { + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + amici::Event("e1", true, true, 1), + } + ); amici::ReturnData r(solver, m); r.id = "some_id"; std::string serialized = amici::serializeToString(r); checkReturnDataEqual( - r, amici::deserializeFromString(serialized)); + r, amici::deserializeFromString(serialized) + ); } -TEST_F(SolverSerializationTest, ToChar) -{ +TEST_F(SolverSerializationTest, ToChar) { int length; char* buf = amici::serializeToChar(solver, &length); - amici::CVodeSolver v = - amici::deserializeFromChar(buf, length); + amici::CVodeSolver v + = amici::deserializeFromChar(buf, length); delete[] buf; ASSERT_EQ(solver, v); } -TEST_F(SolverSerializationTest, ToStdVec) -{ +TEST_F(SolverSerializationTest, ToStdVec) { auto buf = amici::serializeToStdVec(solver); - amici::CVodeSolver v = - amici::deserializeFromChar(buf.data(), buf.size()); + amici::CVodeSolver v = amici::deserializeFromChar( + buf.data(), buf.size() + ); ASSERT_EQ(solver, v); } diff --git a/tests/cpp/unittests/testSplines.cpp b/tests/cpp/unittests/testSplines.cpp index df17cae33f..8ffb1d4b34 100644 --- a/tests/cpp/unittests/testSplines.cpp +++ b/tests/cpp/unittests/testSplines.cpp @@ -7,139 +7,128 @@ #include #include -using std::exp; +using amici::AmiException; using amici::HermiteSpline; using amici::SplineBoundaryCondition; using amici::SplineExtrapolation; -using amici::AmiException; +using std::exp; -#define ASSERT_APPROX(x, x0, rtol, atol) ASSERT_LE(std::abs((x) - (x0)), (atol) + (rtol) * std::abs(x0)) +#define ASSERT_APPROX(x, x0, rtol, atol) \ + ASSERT_LE(std::abs((x) - (x0)), (atol) + (rtol) * std::abs(x0)) void test_spline_values( HermiteSpline const& spline, - std::vector> const& expectations) -{ - for (auto const& [time, expected_value] : expectations) { - ASSERT_DOUBLE_EQ(spline.get_value(time), expected_value); - } + std::vector> const& expectations +) { + for (auto const& [time, expected_value] : expectations) { + ASSERT_DOUBLE_EQ(spline.get_value(time), expected_value); + } } void test_spline_values( HermiteSpline const& spline, std::vector> const& expectations, - const double rtol, const double atol) -{ - for (auto const& [time, expected_value] : expectations) { - ASSERT_APPROX(spline.get_value(time), expected_value, rtol, atol); - } + double const rtol, double const atol +) { + for (auto const& [time, expected_value] : expectations) { + ASSERT_APPROX(spline.get_value(time), expected_value, rtol, atol); + } } void test_spline_sensitivities( HermiteSpline const& spline, - std::vector>> const& expectations) -{ - for (auto const& [time, expected_values] : expectations) { - for (std::vector::size_type ip = 0; ip < expected_values.size(); ip++) - ASSERT_DOUBLE_EQ(spline.get_sensitivity(time, ip), expected_values[ip]); - } + std::vector>> const& expectations +) { + for (auto const& [time, expected_values] : expectations) { + for (std::vector::size_type ip = 0; ip < expected_values.size(); + ip++) + ASSERT_DOUBLE_EQ( + spline.get_sensitivity(time, ip), expected_values[ip] + ); + } } void test_spline_sensitivities( HermiteSpline const& spline, std::vector>> const& expectations, - const double rtol, const double atol) -{ - for (auto const& [time, expected_values] : expectations) { - for (std::vector::size_type ip = 0; ip < expected_values.size(); ip++) - ASSERT_APPROX(spline.get_sensitivity(time, ip), expected_values[ip], rtol, atol); - } + double const rtol, double const atol +) { + for (auto const& [time, expected_values] : expectations) { + for (std::vector::size_type ip = 0; ip < expected_values.size(); + ip++) + ASSERT_APPROX( + spline.get_sensitivity(time, ip), expected_values[ip], rtol, + atol + ); + } } -TEST(Splines, SplineUniform) -{ +TEST(Splines, SplineUniform) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {0.00, 0.0}, - {0.25, 1.74609375}, - {1.0/3, 2.0}, - {0.50, 1.3437499999999996}, - {2.0/3, 0.5}, - {0.75, 0.484375}, - {1.00, 1.0}, + {0.00, 0.0}, {0.25, 1.74609375}, + {1.0 / 3, 2.0}, {0.50, 1.3437499999999996}, + {2.0 / 3, 0.5}, {0.75, 0.484375}, + {1.00, 1.0}, }; test_spline_values(spline, expectations); ASSERT_THROW(spline.get_value(-0.05), AmiException); ASSERT_THROW(spline.get_value(1.05), AmiException); } -TEST(Splines, SplineNonUniform) -{ +TEST(Splines, SplineNonUniform) { // Non-uniform grid - HermiteSpline spline({ 0.0, 0.1, 0.5, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - false, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 0.1, 0.5, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, + SplineBoundaryCondition::given, SplineBoundaryCondition::given, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + false, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {0.00, 0.0}, - {0.05, 1.1484375}, - {0.10, 2.0}, - {0.25, 2.0498046875}, - {0.50, 0.5}, - {0.75, 0.6015625}, - {1.00, 1.0}, + {0.00, 0.0}, {0.05, 1.1484375}, {0.10, 2.0}, {0.25, 2.0498046875}, + {0.50, 0.5}, {0.75, 0.6015625}, {1.00, 1.0}, }; test_spline_values(spline, expectations); ASSERT_THROW(spline.get_value(-0.05), AmiException); ASSERT_THROW(spline.get_value(1.05), AmiException); } -TEST(Splines, SplineExplicit) -{ +TEST(Splines, SplineExplicit) { // Derivatives are given explicitly - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0, 0.75 }, - { 1.0, 0.0, 0.1, -0.1, 0.0 }, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - false, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0, 0.75}, {1.0, 0.0, 0.1, -0.1, 0.0}, + SplineBoundaryCondition::given, SplineBoundaryCondition::given, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + false, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {0.0, 0.0}, - {0.2, 1.8000000000000003}, - {0.25, 2.0}, - {0.4, 1.0243999999999998}, - {0.5, 0.5}, - {0.6, 0.6819999999999999}, - {0.75, 1.0}, - {0.8, 0.9707999999999999}, + {0.0, 0.0}, {0.2, 1.8000000000000003}, + {0.25, 2.0}, {0.4, 1.0243999999999998}, + {0.5, 0.5}, {0.6, 0.6819999999999999}, + {0.75, 1.0}, {0.8, 0.9707999999999999}, {1.0, 0.75}, }; test_spline_values(spline, expectations); @@ -147,436 +136,322 @@ TEST(Splines, SplineExplicit) ASSERT_THROW(spline.get_value(1.05), AmiException); } -TEST(Splines, SplineZeroBC) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::zeroDerivative, - SplineBoundaryCondition::zeroDerivative, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineZeroBC) { + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, + SplineBoundaryCondition::zeroDerivative, + SplineBoundaryCondition::zeroDerivative, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {0.0, 0.0}, - {0.25, 1.65234375}, - {0.5, 1.3437499999999996}, - {0.75, 0.5078125}, - {1.0, 1.0}, + {0.0, 0.0}, {0.25, 1.65234375}, {0.5, 1.3437499999999996}, + {0.75, 0.5078125}, {1.0, 1.0}, }; test_spline_values(spline, expectations); ASSERT_THROW(spline.get_value(-0.05), AmiException); ASSERT_THROW(spline.get_value(1.05), AmiException); } -TEST(Splines, SplineLogarithmic) -{ +TEST(Splines, SplineLogarithmic) { // Logarithmic parametrization - HermiteSpline spline({ 0.0, 1.0 }, - { 0.2, 2.0, 0.5, 1.0, 0.75 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - true); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.2, 2.0, 0.5, 1.0, 0.75}, {}, + SplineBoundaryCondition::given, SplineBoundaryCondition::given, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + true + ); // logarithmic_parametrization // log-space values [-1.60943791, 0.69314718, -0.69314718, 0, -0.28768207] // log-space derivatives [36, 0.3, -4, 0.5, -1.33333333] spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {0.0, 0.2}, - {0.2, 2.07939779651678}, - {0.25, 2.0}, - {0.4, 0.947459046694449}, - {0.5, 0.5}, - {0.6, 0.545987404053269}, - {0.75, 1.0}, - {0.8, 0.996753014029391}, - {1.0, 0.75}, + {0.0, 0.2}, {0.2, 2.07939779651678}, + {0.25, 2.0}, {0.4, 0.947459046694449}, + {0.5, 0.5}, {0.6, 0.545987404053269}, + {0.75, 1.0}, {0.8, 0.996753014029391}, + {1.0, 0.75}, }; test_spline_values(spline, expectations, 1e-14, 0.0); ASSERT_THROW(spline.get_value(-0.05), AmiException); ASSERT_THROW(spline.get_value(1.05), AmiException); } -TEST(Splines, SplineUniformConstantExtrapolation) -{ +TEST(Splines, SplineUniformConstantExtrapolation) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::constant, - SplineExtrapolation::constant, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::constant, + SplineExtrapolation::constant, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {-2.00, 0.0}, - {-1.00, 0.0}, - { 0.00, 0.0}, - { 0.25, 1.74609375}, - { 1.0/3, 2.0}, - { 0.50, 1.3437499999999996}, - { 2.0/3, 0.5}, - { 0.75, 0.484375}, - { 1.00, 1.0}, - { 2.00, 1.0}, - { 3.00, 1.0}, + {-2.00, 0.0}, {-1.00, 0.0}, {0.00, 0.0}, + {0.25, 1.74609375}, {1.0 / 3, 2.0}, {0.50, 1.3437499999999996}, + {2.0 / 3, 0.5}, {0.75, 0.484375}, {1.00, 1.0}, + {2.00, 1.0}, {3.00, 1.0}, }; test_spline_values(spline, expectations); } -TEST(Splines, SplineUniformLinearExtrapolation) -{ +TEST(Splines, SplineUniformLinearExtrapolation) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::linear, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::linear, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {-2.00, -12.0}, - {-1.00, -6.0}, - { 0.00, 0.0}, - { 0.25, 1.74609375}, - { 1.0/3, 2.0}, - { 0.50, 1.3437499999999996}, - { 2.0/3, 0.5}, - { 0.75, 0.484375}, - { 1.00, 1.0}, - { 2.00, 2.5}, - { 3.00, 4.0}, + {-2.00, -12.0}, {-1.00, -6.0}, {0.00, 0.0}, + {0.25, 1.74609375}, {1.0 / 3, 2.0}, {0.50, 1.3437499999999996}, + {2.0 / 3, 0.5}, {0.75, 0.484375}, {1.00, 1.0}, + {2.00, 2.5}, {3.00, 4.0}, }; test_spline_values(spline, expectations); } -TEST(Splines, SplineUniformPolynomialExtrapolation) -{ +TEST(Splines, SplineUniformPolynomialExtrapolation) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 0.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::polynomial, - SplineExtrapolation::polynomial, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {0.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::polynomial, + SplineExtrapolation::polynomial, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {-2.00, 429.0}, - {-1.00, 57.0}, - { 0.00, 0.0}, - { 0.25, 1.74609375}, - { 1.0/3, 2.0}, - { 0.50, 1.3437499999999996}, - { 2.0/3, 0.5}, - { 0.75, 0.484375}, - { 1.00, 1.0}, - { 2.00, -33.5}, - { 3.00, -248.0}, + {-2.00, 429.0}, {-1.00, 57.0}, {0.00, 0.0}, + {0.25, 1.74609375}, {1.0 / 3, 2.0}, {0.50, 1.3437499999999996}, + {2.0 / 3, 0.5}, {0.75, 0.484375}, {1.00, 1.0}, + {2.00, -33.5}, {3.00, -248.0}, }; test_spline_values(spline, expectations); } -TEST(Splines, SplineUniformPeriodicExtrapolation) -{ +TEST(Splines, SplineUniformPeriodicExtrapolation) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 1.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::periodic, - SplineBoundaryCondition::periodic, - SplineExtrapolation::periodic, - SplineExtrapolation::periodic, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {1.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::periodic, + SplineBoundaryCondition::periodic, SplineExtrapolation::periodic, + SplineExtrapolation::periodic, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {-4.0/3, 0.5}, - {-0.50, 1.2812499999999996}, - { 0.00, 1.0}, - { 0.25, 1.9140625}, - { 1.0/3, 2.0}, - { 0.50, 1.2812499999999996}, - { 2.0/3, 0.5}, - { 0.75, 0.47265625}, - { 1.00, 1.0}, - { 1.25, 1.9140625}, - { 2.75, 0.47265625}, + {-4.0 / 3, 0.5}, {-0.50, 1.2812499999999996}, + {0.00, 1.0}, {0.25, 1.9140625}, + {1.0 / 3, 2.0}, {0.50, 1.2812499999999996}, + {2.0 / 3, 0.5}, {0.75, 0.47265625}, + {1.00, 1.0}, {1.25, 1.9140625}, + {2.75, 0.47265625}, }; test_spline_values(spline, expectations); } -TEST(Splines, SplineNonUniformPeriodicExtrapolation) -{ +TEST(Splines, SplineNonUniformPeriodicExtrapolation) { // Non-uniform grid - HermiteSpline spline({ 0.0, 0.1, 0.5, 1.0 }, - { 1.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::periodic, - SplineBoundaryCondition::periodic, - SplineExtrapolation::periodic, - SplineExtrapolation::periodic, - true, // node_derivative_by_FD - false, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 0.1, 0.5, 1.0}, {1.0, 2.0, 0.5, 1.0}, {}, + SplineBoundaryCondition::periodic, SplineBoundaryCondition::periodic, + SplineExtrapolation::periodic, SplineExtrapolation::periodic, + true, // node_derivative_by_FD + false, // equidistant_spacing + false + ); // logarithmic_parametrization spline.compute_coefficients(); std::vector> expectations = { // t, expected value - {-1.90, 2.0}, - {-0.25, 0.3203125}, - { 0.00, 1.0}, - { 0.05, 1.5296875}, - { 0.10, 2.0}, - { 0.25, 1.7568359375}, - { 0.50, 0.5}, - { 0.75, 0.3203125}, - { 1.00, 1.0}, - { 1.50, 0.5}, - { 2.05, 1.5296875}, + {-1.90, 2.0}, {-0.25, 0.3203125}, {0.00, 1.0}, + {0.05, 1.5296875}, {0.10, 2.0}, {0.25, 1.7568359375}, + {0.50, 0.5}, {0.75, 0.3203125}, {1.00, 1.0}, + {1.50, 0.5}, {2.05, 1.5296875}, }; test_spline_values(spline, expectations, 1e-14, 0.0); } -TEST(Splines, SplineUniformSensitivity) -{ +TEST(Splines, SplineUniformSensitivity) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); std::vector>> expectations = { // t, expected values of sensitivities - {0.00, {3.0, 1.0, 0.0}}, - {0.25, {0.539062, 0.179688, 4.45312}}, - {1.0/3, {0.0, 0.0, 5.0}}, - {0.50, {0.1875, -0.125, 2.625}}, - {2.0/3, {0.0, 0.0, 0.0}}, - {0.75, {-1.07812, 0.179688, 0.1875}}, - {1.00, {-6.0, 1.0, 3.0}}, + {0.00, {3.0, 1.0, 0.0}}, {0.25, {0.539062, 0.179688, 4.45312}}, + {1.0 / 3, {0.0, 0.0, 5.0}}, {0.50, {0.1875, -0.125, 2.625}}, + {2.0 / 3, {0.0, 0.0, 0.0}}, {0.75, {-1.07812, 0.179688, 0.1875}}, + {1.00, {-6.0, 1.0, 3.0}}, }; test_spline_sensitivities(spline, expectations, 1e-5, 1e-6); ASSERT_THROW(spline.get_sensitivity(-0.05, 0), AmiException); - ASSERT_THROW(spline.get_sensitivity( 1.05, 1), AmiException); + ASSERT_THROW(spline.get_sensitivity(1.05, 1), AmiException); } -TEST(Splines, SplineNonUniformSensitivity) -{ - HermiteSpline spline({ 0.0, 0.1, 0.5, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - false, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineNonUniformSensitivity) { + HermiteSpline spline( + {0.0, 0.1, 0.5, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, + SplineBoundaryCondition::given, SplineBoundaryCondition::given, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + false, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); std::vector>> expectations = { // t, expected values of sensitivities - {0.00, { 3.0, 1.0, 0.0}}, - {0.05, { 1.3125, 0.4375, 2.89062}}, - {0.10, { 0.0, 0.0, 5.0}}, - {0.30, {-0.45, -0.3, 3.6}}, - {0.50, { 0.0, 0.0, 0.0}}, - {0.75, {-2.625, 0.4375, 0.921875}}, - {1.00, {-6.0, 1.0, 3.0}}, + {0.00, {3.0, 1.0, 0.0}}, {0.05, {1.3125, 0.4375, 2.89062}}, + {0.10, {0.0, 0.0, 5.0}}, {0.30, {-0.45, -0.3, 3.6}}, + {0.50, {0.0, 0.0, 0.0}}, {0.75, {-2.625, 0.4375, 0.921875}}, + {1.00, {-6.0, 1.0, 3.0}}, }; test_spline_sensitivities(spline, expectations, 1e-5, 1e-6); ASSERT_THROW(spline.get_sensitivity(-0.05, 0), AmiException); - ASSERT_THROW(spline.get_sensitivity( 1.05, 1), AmiException); + ASSERT_THROW(spline.get_sensitivity(1.05, 1), AmiException); } -TEST(Splines, SplineExplicitSensitivity) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - { 13.625, 7.5, 1.1585290151921035, 1.0 }, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - false, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineExplicitSensitivity) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, + {13.625, 7.5, 1.1585290151921035, 1.0}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + false, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; std::vector dslopesdp = { - 0.0, 0.0, 18.75, - 0.0, 1.0, 3.0, - 4.0, -0.540302, 0.0, - 0.0, 0.0, 0.0, + 0.0, 0.0, 18.75, 0.0, 1.0, 3.0, 4.0, -0.540302, 0.0, 0.0, 0.0, 0.0, }; spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); std::vector>> expectations = { // t, expected values of sensitivities - {0.00, { 3.0, 1.0, 0.0}}, - {0.25, { 0.46875, 0.109375, 4.37109}}, - {1.0/3, { 0.0, 0.0, 5.0}}, - {0.50, {-0.166667, 0.0641793, 2.625}}, - {2.0/3, { 0.0, 0.0, 0.0}}, - {0.75, {-0.75, 0.130923, 0.46875}}, - {1.00, {-6.0, 1.0, 3.0}}, + {0.00, {3.0, 1.0, 0.0}}, {0.25, {0.46875, 0.109375, 4.37109}}, + {1.0 / 3, {0.0, 0.0, 5.0}}, {0.50, {-0.166667, 0.0641793, 2.625}}, + {2.0 / 3, {0.0, 0.0, 0.0}}, {0.75, {-0.75, 0.130923, 0.46875}}, + {1.00, {-6.0, 1.0, 3.0}}, }; test_spline_sensitivities(spline, expectations, 1e-5, 0.0); ASSERT_THROW(spline.get_sensitivity(-0.05, 0), AmiException); - ASSERT_THROW(spline.get_sensitivity( 1.05, 1), AmiException); + ASSERT_THROW(spline.get_sensitivity(1.05, 1), AmiException); } -TEST(Splines, SplineZeroDerivativeSensitivity) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::zeroDerivative, - SplineBoundaryCondition::zeroDerivative, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineZeroDerivativeSensitivity) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, + SplineBoundaryCondition::zeroDerivative, + SplineBoundaryCondition::zeroDerivative, + SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); std::vector>> expectations = { // t, expected values of sensitivities - {0.00, { 3.0, 1.0, 0.0}}, - {0.25, { 0.679688, 0.226562, 4.21875}}, - {1.0/3, { 0.0, 0.0, 5.0}}, - {0.50, {0.1875, -0.125, 2.625}}, - {2.0/3, { 0.0, 0.0, 0.0}}, - {0.75, {-1.35938, 0.226562, 0.328125}}, - {1.00, {-6.0, 1.0, 3.0}}, + {0.00, {3.0, 1.0, 0.0}}, {0.25, {0.679688, 0.226562, 4.21875}}, + {1.0 / 3, {0.0, 0.0, 5.0}}, {0.50, {0.1875, -0.125, 2.625}}, + {2.0 / 3, {0.0, 0.0, 0.0}}, {0.75, {-1.35938, 0.226562, 0.328125}}, + {1.00, {-6.0, 1.0, 3.0}}, }; test_spline_sensitivities(spline, expectations, 1e-5, 0.0); ASSERT_THROW(spline.get_sensitivity(-0.05, 0), AmiException); - ASSERT_THROW(spline.get_sensitivity( 1.05, 1), AmiException); + ASSERT_THROW(spline.get_sensitivity(1.05, 1), AmiException); } -TEST(Splines, SplineLogarithmicSensitivity) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::noExtrapolation, - true, // node_derivative_by_FD - true, // equidistant_spacing - true); // logarithmic_parametrization +TEST(Splines, SplineLogarithmicSensitivity) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::noExtrapolation, + true, // node_derivative_by_FD + true, // equidistant_spacing + true + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); std::vector>> expectations = { // t, expected values of sensitivities - {0.00, { 3.0, 1.0, 0.0}}, - {0.25, { 0.585881, 0.195294, 4.38532}}, - {1.0/3, { 0.0, 0.0, 5.0}}, - {0.50, { 0.514003, -0.132395, 1.52044}}, - {2.0/3, { 0.0, 0.0, 0.0}}, - {0.75, {-0.820743, 0.13679, -0.0577988}}, - {1.00, {-6.0, 1.0, 3.0}}, + {0.00, {3.0, 1.0, 0.0}}, {0.25, {0.585881, 0.195294, 4.38532}}, + {1.0 / 3, {0.0, 0.0, 5.0}}, {0.50, {0.514003, -0.132395, 1.52044}}, + {2.0 / 3, {0.0, 0.0, 0.0}}, {0.75, {-0.820743, 0.13679, -0.0577988}}, + {1.00, {-6.0, 1.0, 3.0}}, }; test_spline_sensitivities(spline, expectations, 1e-6, 1e-6); ASSERT_THROW(spline.get_sensitivity(-0.05, 0), AmiException); - ASSERT_THROW(spline.get_sensitivity( 1.05, 1), AmiException); + ASSERT_THROW(spline.get_sensitivity(1.05, 1), AmiException); } -TEST(Splines, SplineFinalValue_ConstantExtrapolation) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::constant, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_ConstantExtrapolation) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::constant, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -588,25 +463,18 @@ TEST(Splines, SplineFinalValue_ConstantExtrapolation) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 3.0); } -TEST(Splines, SplineFinalValue_LinearExtrapolationPositiveDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LinearExtrapolationPositiveDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -618,25 +486,18 @@ TEST(Splines, SplineFinalValue_LinearExtrapolationPositiveDerivative) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 0.0); } -TEST(Splines, SplineFinalValue_LinearExtrapolationNegativeDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 0.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LinearExtrapolationNegativeDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 0.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -648,25 +509,18 @@ TEST(Splines, SplineFinalValue_LinearExtrapolationNegativeDerivative) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 0.0); } -TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 1.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 1.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -678,25 +532,18 @@ TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivative) ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(2))); } -TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivativeByBC) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 2.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::zeroDerivative, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivativeByBC) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 2.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::zeroDerivative, + SplineExtrapolation::noExtrapolation, SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -708,88 +555,69 @@ TEST(Splines, SplineFinalValue_LinearExtrapolationZeroDerivativeByBC) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 3.0); } -TEST(Splines, SplineFinalValue_PolynomialExtrapolationPositive) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { -8.0, -6.0, -1.0, -2.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::polynomial, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_PolynomialExtrapolationPositive) { + HermiteSpline spline( + {0.0, 1.0}, {-8.0, -6.0, -1.0, -2.0}, {}, + SplineBoundaryCondition::given, SplineBoundaryCondition::given, + SplineExtrapolation::noExtrapolation, SplineExtrapolation::polynomial, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); spline.compute_final_value(); spline.compute_final_sensitivity(n_params, 0, dvaluesdp, dslopesdp); ASSERT_DOUBLE_EQ(spline.get_final_value(), INFINITY); - /* NB sensitivities for this case are not implemented, since they are unlikely to be used*/ + /* NB sensitivities for this case are not implemented, since they are + * unlikely to be used*/ ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(0))); ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(1))); ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(2))); } -TEST(Splines, SplineFinalValue_PolynomialExtrapolationNegative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 2.0 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::polynomial, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_PolynomialExtrapolationNegative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 2.0}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::polynomial, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); spline.compute_final_value(); spline.compute_final_sensitivity(n_params, 0, dvaluesdp, dslopesdp); ASSERT_DOUBLE_EQ(spline.get_final_value(), -INFINITY); - /* NB sensitivities for this case are not implemented, since they are unlikely to be used*/ + /* NB sensitivities for this case are not implemented, since they are + * unlikely to be used*/ ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(0))); ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(1))); ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(2))); } -TEST(Splines, SplineFinalValue_PeriodicExtrapolation) -{ +TEST(Splines, SplineFinalValue_PeriodicExtrapolation) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 1.0, 2.0, 0.5, 1.0 }, - {}, - SplineBoundaryCondition::periodic, - SplineBoundaryCondition::periodic, - SplineExtrapolation::periodic, - SplineExtrapolation::periodic, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {1.0, 2.0, 0.5, 1.0}, {}, SplineBoundaryCondition::periodic, + SplineBoundaryCondition::periodic, SplineExtrapolation::periodic, + SplineExtrapolation::periodic, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -801,26 +629,19 @@ TEST(Splines, SplineFinalValue_PeriodicExtrapolation) ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(2))); } -TEST(Splines, SplineFinalValue_PeriodicExtrapolationConstant) -{ +TEST(Splines, SplineFinalValue_PeriodicExtrapolationConstant) { // Uniform grid - HermiteSpline spline({ 0.0, 1.0 }, - { 1.0, 1.0, 1.0, 1.0 }, - {}, - SplineBoundaryCondition::periodic, - SplineBoundaryCondition::periodic, - SplineExtrapolation::periodic, - SplineExtrapolation::periodic, - true, // node_derivative_by_FD - true, // equidistant_spacing - false); // logarithmic_parametrization + HermiteSpline spline( + {0.0, 1.0}, {1.0, 1.0, 1.0, 1.0}, {}, SplineBoundaryCondition::periodic, + SplineBoundaryCondition::periodic, SplineExtrapolation::periodic, + SplineExtrapolation::periodic, + true, // node_derivative_by_FD + true, // equidistant_spacing + false + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -832,25 +653,18 @@ TEST(Splines, SplineFinalValue_PeriodicExtrapolationConstant) ASSERT_TRUE(std::isnan(spline.get_final_sensitivity(2))); } -TEST(Splines, SplineFinalValue_LogarithmicPositiveDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 4.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - true); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LogarithmicPositiveDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 4.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + true + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -862,25 +676,18 @@ TEST(Splines, SplineFinalValue_LogarithmicPositiveDerivative) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 0.0); } -TEST(Splines, SplineFinalValue_LogarithmicNegativeDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 0.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::linear, - true, // node_derivative_by_FD - true, // equidistant_spacing - true); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LogarithmicNegativeDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 0.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::linear, + true, // node_derivative_by_FD + true, // equidistant_spacing + true + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp); @@ -892,25 +699,18 @@ TEST(Splines, SplineFinalValue_LogarithmicNegativeDerivative) ASSERT_DOUBLE_EQ(spline.get_final_sensitivity(2), 0.0); } -TEST(Splines, SplineFinalValue_LogarithmicZeroDerivative) -{ - HermiteSpline spline({ 0.0, 1.0 }, - { 2.5, 3.25, 1.0, 0.5 }, - {}, - SplineBoundaryCondition::given, - SplineBoundaryCondition::given, - SplineExtrapolation::noExtrapolation, - SplineExtrapolation::constant, - true, // node_derivative_by_FD - true, // equidistant_spacing - true); // logarithmic_parametrization +TEST(Splines, SplineFinalValue_LogarithmicZeroDerivative) { + HermiteSpline spline( + {0.0, 1.0}, {2.5, 3.25, 1.0, 0.5}, {}, SplineBoundaryCondition::given, + SplineBoundaryCondition::given, SplineExtrapolation::noExtrapolation, + SplineExtrapolation::constant, + true, // node_derivative_by_FD + true, // equidistant_spacing + true + ); // logarithmic_parametrization int n_params = 3; - std::vector dvaluesdp = { - 3.0, 1.0, 0.0, - 0.0, 0.0, 5.0, - 0.0, 0.0, 0.0, - -6.0, 1.0, 3.0 - }; + std::vector dvaluesdp + = {3.0, 1.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, -6.0, 1.0, 3.0}; auto dslopesdp = std::vector(spline.n_nodes() * n_params); spline.compute_coefficients(); spline.compute_coefficients_sensi(n_params, 0, dvaluesdp, dslopesdp);