From 4f56ebe00b37a814fbc67531ee338bc591f9bc02 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 2 Oct 2025 15:48:33 +0200 Subject: [PATCH] Clean up symbolic_functions.{cpp,h} Remove MATLAB leftovers. Clean up. Related to #2727. --- include/amici/edata.h | 1 + include/amici/symbolic_functions.h | 90 ++++++------------------------ src/edata.cpp | 10 ++-- src/rdata.cpp | 9 +-- src/solver.cpp | 34 +++++------ src/symbolic_functions.cpp | 62 ++------------------ tests/cpp/unittests/testMisc.cpp | 30 ---------- 7 files changed, 51 insertions(+), 185 deletions(-) diff --git a/include/amici/edata.h b/include/amici/edata.h index 50d0545fcf..b0ee3a11a1 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -5,6 +5,7 @@ #include "amici/misc.h" #include "amici/simulation_parameters.h" +#include #include namespace amici { diff --git a/include/amici/symbolic_functions.h b/include/amici/symbolic_functions.h index ba6bf05b74..c065a74971 100644 --- a/include/amici/symbolic_functions.h +++ b/include/amici/symbolic_functions.h @@ -1,119 +1,63 @@ -#ifndef amici_symbolic_functions_h -#define amici_symbolic_functions_h +#ifndef AMICI_SYMBOLIC_FUNCTIONS_H +#define AMICI_SYMBOLIC_FUNCTIONS_H namespace amici { /** - * C implementation of log function, this prevents returning NaN values for - * negative values + * @brief A log function, that prevents returning NaN values for + * negative values. * * @param x argument * @return if(x>0) then log(x) else -Inf - * */ double log(double x); /** - * C implementation of matlab function dirac + * @brief The Dirac delta function. * * @param x argument * @return if(x==0) then INF else 0 * */ double dirac(double x); -double heaviside(double x, double x0); - -/** - * c implementation of matlab function min - * - * @param a value1 - * @param b value2 - * @param c bogus parameter to ensure correct parsing as a function - * @return if(a < b) then a else b - * - */ -double min(double a, double b, double c); - -/** - * parameter derivative of c implementation of matlab function max - * - * @param id argument index for differentiation - * @param a value1 - * @param b value2 - * @param c bogus parameter to ensure correct parsing as a function - * @return id == 1: if(a > b) then 1 else 0 - * @return id == 2: if(a > b) then 0 else 1 - * - */ -double Dmin(int id, double a, double b, double c); - -/** - * c implementation of matlab function max - * - * @param a value1 - * @param b value2 - * @return if(a > b) then a else b - * - */ -double max(double a, double b, double c); /** - * parameter derivative of c implementation of matlab function max + * The Heaviside function. * - * @param id argument index for differentiation - * @param a value1 - * @param b value2 - * @return id == 1: if(a > b) then 1 else 0 - * @return id == 2: if(a > b) then 0 else 1 + * @param x argument + * @param x0 value at x==0 + * @return if(x>0) then 1 else if (x==0) then x0 else 0 * */ -double Dmax(int id, double a, double b, double c); +double heaviside(double x, double x0); /** - * specialized pow functions that assumes positivity of the first argument + * @brief Specialized pow functions that assumes positivity of the first + * argument. * * @param base base * @param exponent exponent - * @return pow(max(base,0.0),exponen) + * @return pow(max(base,0.0), exponent) * */ double pos_pow(double base, double exponent); /** - * c++ interface to the isNaN function - * - * @param what argument - * @return isnan(what) - * - */ -int isNaN(double what); - -/** - * c++ interface to the isinf function - * - * @param what argument - * @return isnan(what) - * - */ -int isInf(double what); - -/** - * function returning nan + * @brief Get NaN value. * * @return NaN - * */ double getNaN(); /** - * c implementation of matlab function sign + * @brief The sign function. * * @param x argument - * @return 0 + * @return if(x>0) then 1 else if (x<0) then -1 else 0 * */ double sign(double x); } // namespace amici -#endif /* amici_symbolic_functions_h */ +#endif // AMICI_SYMBOLIC_FUNCTIONS_H diff --git a/src/edata.cpp b/src/edata.cpp index f1908cf1a6..d26cb6afb7 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -10,6 +10,8 @@ namespace amici { +using std::isnan; + ExpData::ExpData(int const nytrue, int const nztrue, int const nmaxevent) : nytrue_(nytrue) , nztrue_(nztrue) @@ -171,7 +173,7 @@ void ExpData::setObservedData( bool ExpData::isSetObservedData(int const it, int const iy) const { return !observed_data_.empty() - && !isNaN(observed_data_.at(it * nytrue_ + iy)); + && !isnan(observed_data_.at(it * nytrue_ + iy)); } std::vector const& ExpData::getObservedData() const { @@ -225,7 +227,7 @@ void ExpData::setObservedDataStdDev(realtype const stdDev, int const iy) { bool ExpData::isSetObservedDataStdDev(int const it, int const iy) const { return !observed_data_std_dev_.empty() - && !isNaN(observed_data_std_dev_.at(it * nytrue_ + iy)); + && !isnan(observed_data_std_dev_.at(it * nytrue_ + iy)); } std::vector const& ExpData::getObservedDataStdDev() const { @@ -265,7 +267,7 @@ void ExpData::setObservedEvents( bool ExpData::isSetObservedEvents(int const ie, int const iz) const { return !observed_events_.empty() - && !isNaN(observed_events_.at(ie * nztrue_ + iz)); + && !isnan(observed_events_.at(ie * nztrue_ + iz)); } std::vector const& ExpData::getObservedEvents() const { @@ -321,7 +323,7 @@ void ExpData::setObservedEventsStdDev(realtype const stdDev, int const iz) { bool ExpData::isSetObservedEventsStdDev(int const ie, int const iz) const { if (!observed_events_std_dev_.empty()) // avoid out of bound memory access - return !isNaN(observed_events_std_dev_.at(ie * nztrue_ + iz)); + return !isnan(observed_events_std_dev_.at(ie * nztrue_ + iz)); return false; } diff --git a/src/rdata.cpp b/src/rdata.cpp index 3e09aede03..01c6831df7 100644 --- a/src/rdata.cpp +++ b/src/rdata.cpp @@ -14,6 +14,7 @@ #include namespace amici { +using std::isnan; ReturnData::ReturnData(Solver const& solver, Model const& model) : ReturnData( @@ -346,7 +347,7 @@ void ReturnData::getDataOutput( model.getObservableSigma(slice(sigmay, it, ny), it, edata); if (edata) { - if (!isNaN(llh)) + if (!isnan(llh)) model.addObservableObjective(llh, it, sol.x, *edata); fres(it, model, sol, *edata); fchi2(it, *edata); @@ -422,14 +423,14 @@ void ReturnData::getEventOutput( slice(sigmaz, nroots_.at(ie), nz), ie, nroots_.at(ie), sol.t, edata ); - if (!isNaN(llh)) + if (!isnan(llh)) model.addEventObjective( llh, ie, nroots_.at(ie), sol.t, sol.x, *edata ); /* if called from fillEvent at last timepoint, add regularization based on rz */ - if (sol.t == model.getTimepoint(nt - 1) && !isNaN(llh)) { + if (sol.t == model.getTimepoint(nt - 1) && !isnan(llh)) { model.addEventObjectiveRegularization( llh, ie, nroots_.at(ie), sol.t, sol.x, *edata ); @@ -951,7 +952,7 @@ void ReturnData::fres( } void ReturnData::fchi2(int const it, ExpData const& edata) { - if (res.empty() || isNaN(chi2)) + if (res.empty() || isnan(chi2)) return; for (int iy = 0; iy < nytrue; ++iy) { diff --git a/src/solver.cpp b/src/solver.cpp index 7872a63b9b..0ccf8dd0a1 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -3,7 +3,6 @@ #include "amici/amici.h" #include "amici/exception.h" #include "amici/model.h" -#include "amici/symbolic_functions.h" #include @@ -12,6 +11,7 @@ #include namespace amici { +using std::isnan; /* Error handler passed to SUNDIALS. */ void wrapErrHandlerFn( @@ -604,11 +604,11 @@ bool operator==(Solver const& a, Solver const& b) { && (a.getRelativeToleranceSteadyStateSensi() == b.getRelativeToleranceSteadyStateSensi()) && (a.rtol_fsa_ == b.rtol_fsa_ - || (isNaN(a.rtol_fsa_) && isNaN(b.rtol_fsa_))) + || (isnan(a.rtol_fsa_) && isnan(b.rtol_fsa_))) && (a.atol_fsa_ == b.atol_fsa_ - || (isNaN(a.atol_fsa_) && isNaN(b.atol_fsa_))) - && (a.rtolB_ == b.rtolB_ || (isNaN(a.rtolB_) && isNaN(b.rtolB_))) - && (a.atolB_ == b.atolB_ || (isNaN(a.atolB_) && isNaN(b.atolB_))) + || (isnan(a.atol_fsa_) && isnan(b.atol_fsa_))) + && (a.rtolB_ == b.rtolB_ || (isnan(a.rtolB_) && isnan(b.rtolB_))) + && (a.atolB_ == b.atolB_ || (isnan(a.atolB_) && isnan(b.atolB_))) && (a.sensi_ == b.sensi_) && (a.sensi_meth_ == b.sensi_meth_) && (a.newton_step_steadystate_conv_ == b.newton_step_steadystate_conv_) @@ -672,8 +672,8 @@ void Solver::applyQuadTolerancesASA(int const which) const { if (sensi_ < SensitivityOrder::first) return; - realtype const quad_rtol = isNaN(quad_rtol_) ? rtol_ : quad_rtol_; - realtype const quad_atol = isNaN(quad_atol_) ? atol_ : quad_atol_; + realtype const quad_rtol = isnan(quad_rtol_) ? rtol_ : quad_rtol_; + realtype const quad_atol = isnan(quad_atol_) ? atol_ : quad_atol_; /* Enable Quadrature Error Control */ setQuadErrConB(which, !std::isinf(quad_atol) && !std::isinf(quad_rtol)); @@ -691,8 +691,8 @@ void Solver::applyQuadTolerances() const { if (sensi_ < SensitivityOrder::first) return; - realtype quad_rtolF = isNaN(quad_rtol_) ? rtol_ : quad_rtol_; - realtype quad_atolF = isNaN(quad_atol_) ? atol_ : quad_atol_; + realtype quad_rtolF = isnan(quad_rtol_) ? rtol_ : quad_rtol_; + realtype quad_atolF = isnan(quad_atol_) ? atol_ : quad_atol_; /* Enable Quadrature Error Control */ setQuadErrCon(!std::isinf(quad_atolF) && !std::isinf(quad_rtolF)); @@ -863,7 +863,7 @@ void Solver::setAbsoluteTolerance(double const atol) { } double Solver::getRelativeToleranceFSA() const { - return static_cast(isNaN(rtol_fsa_) ? rtol_ : rtol_fsa_); + return static_cast(isnan(rtol_fsa_) ? rtol_ : rtol_fsa_); } void Solver::setRelativeToleranceFSA(double const rtol) { @@ -878,7 +878,7 @@ void Solver::setRelativeToleranceFSA(double const rtol) { } double Solver::getAbsoluteToleranceFSA() const { - return static_cast(isNaN(atol_fsa_) ? atol_ : atol_fsa_); + return static_cast(isnan(atol_fsa_) ? atol_ : atol_fsa_); } void Solver::setAbsoluteToleranceFSA(double const atol) { @@ -893,7 +893,7 @@ void Solver::setAbsoluteToleranceFSA(double const atol) { } double Solver::getRelativeToleranceB() const { - return static_cast(isNaN(rtolB_) ? rtol_ : rtolB_); + return static_cast(isnan(rtolB_) ? rtol_ : rtolB_); } void Solver::setRelativeToleranceB(double const rtol) { @@ -908,7 +908,7 @@ void Solver::setRelativeToleranceB(double const rtol) { } double Solver::getAbsoluteToleranceB() const { - return static_cast(isNaN(atolB_) ? atol_ : atolB_); + return static_cast(isnan(atolB_) ? atol_ : atolB_); } void Solver::setAbsoluteToleranceB(double const atol) { @@ -971,7 +971,7 @@ void Solver::setSteadyStateToleranceFactor(double const factor) { double Solver::getRelativeToleranceSteadyState() const { return static_cast( - isNaN(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_ + isnan(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_ ); } @@ -984,7 +984,7 @@ void Solver::setRelativeToleranceSteadyState(double const rtol) { double Solver::getAbsoluteToleranceSteadyState() const { return static_cast( - isNaN(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_ + isnan(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_ ); } @@ -1010,7 +1010,7 @@ void Solver::setSteadyStateSensiToleranceFactor( double Solver::getRelativeToleranceSteadyStateSensi() const { return static_cast( - isNaN(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_ + isnan(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_ ); } @@ -1023,7 +1023,7 @@ void Solver::setRelativeToleranceSteadyStateSensi(double const rtol) { double Solver::getAbsoluteToleranceSteadyStateSensi() const { return static_cast( - isNaN(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_ + isnan(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_ ); } diff --git a/src/symbolic_functions.cpp b/src/symbolic_functions.cpp index 609d44ef70..bd6d68662b 100644 --- a/src/symbolic_functions.cpp +++ b/src/symbolic_functions.cpp @@ -8,49 +8,27 @@ #include "amici/symbolic_functions.h" #include -#include #include -#include -#if _MSC_VER && !__INTEL_COMPILER -#include -#define alloca _alloca -#elif defined(WIN32) || defined(__WIN32) || defined(__WIN32__) -// For alloca(). -#include -#else -#include -#endif +#include namespace amici { -int isNaN(double what) { return std::isnan(what); } - -int isInf(double what) { return std::isinf(what); } - -double getNaN() { return NAN; } +double getNaN() { return std::numeric_limits::quiet_NaN(); } double log(double x) { if (x <= 0) { - return -std::log(DBL_MAX); + return -std::log(std::numeric_limits::max()); } return std::log(x); } double dirac(double x) { if (x == 0.0) { - return DBL_MAX; + return std::numeric_limits::max(); } return 0.0; } -/** - * c implementation of matlab function heaviside - * - * @param x argument - * @param x0 value at x==0 - * @return if(x>0) then 1 else if (x==0) then x0 else 0 - * - */ double heaviside(double x, double x0) { if (x < 0.0) return 0.0; @@ -69,38 +47,8 @@ double sign(double x) { return 0.0; } -double max(double a, double b, double /*c*/) { - int anan = isNaN(a), bnan = isNaN(b); - if (anan || bnan) { - if (anan && !bnan) - return b; - if (!anan && bnan) - return a; - return a; - } - return (std::max(a, b)); -} - -double min(double a, double b, double c) { return (-max(-a, -b, c)); } - -double Dmax(int id, double a, double b, double /*c*/) { - if (id == 1.0) { - if (a > b) - return 1.0; - return 0.0; - } - if (a > b) { - return 0.0; - } - return 1.0; -} - -double Dmin(int id, double a, double b, double c) { - return Dmax(id, -a, -b, c); -} - double pos_pow(double base, double exponent) { - // we do NOT want to propagate NaN values here, if base is nan, so should + // we do NOT want to propagate NaN values here, if base is NaN, so should // the output be return pow(std::max(base, 0.0), exponent); } diff --git a/tests/cpp/unittests/testMisc.cpp b/tests/cpp/unittests/testMisc.cpp index a6f18035a1..ebd1475369 100644 --- a/tests/cpp/unittests/testMisc.cpp +++ b/tests/cpp/unittests/testMisc.cpp @@ -178,36 +178,6 @@ TEST(SymbolicFunctionsTest, Heaviside) { ASSERT_EQ(1, heaviside(1, 0.5)); } -TEST(SymbolicFunctionsTest, Min) { - ASSERT_EQ(-1, min(-1, 2, 0)); - ASSERT_EQ(-2, min(1, -2, 0)); - ASSERT_TRUE(isNaN(min(getNaN(), getNaN(), 0))); - ASSERT_EQ(-1, min(-1, getNaN(), 0)); - ASSERT_EQ(-1, min(getNaN(), -1, 0)); -} - -TEST(SymbolicFunctionsTest, Max) { - ASSERT_EQ(2, max(-1, 2, 0)); - ASSERT_EQ(1, max(1, -2, 0)); - ASSERT_TRUE(isNaN(max(getNaN(), getNaN(), 0))); - ASSERT_EQ(-1, max(-1, getNaN(), 0)); - ASSERT_EQ(-1, max(getNaN(), -1, 0)); -} - -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) { - 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) { ASSERT_EQ(0, pos_pow(-0.1, 3)); ASSERT_EQ(pow(0.1, 3), pos_pow(0.1, 3));