Skip to content

Clean up symbolic_functions.{cpp,h} #2956

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dweindl merged 1 commit into from
Oct 3, 2025
Merged

Clean up symbolic_functions.{cpp,h} #2956

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions include/amici/edata.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "amici/misc.h"
#include "amici/simulation_parameters.h"

#include <string>
#include <vector>

namespace amici {
Expand Down
90 changes: 17 additions & 73 deletions include/amici/symbolic_functions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
10 changes: 6 additions & 4 deletions src/edata.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,8 @@

namespace amici {

using std::isnan;

ExpData::ExpData(int const nytrue, int const nztrue, int const nmaxevent)
: nytrue_(nytrue)
, nztrue_(nztrue)
Expand Down Expand Up @@ -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<realtype> const& ExpData::getObservedData() const {
Expand Down Expand Up @@ -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<realtype> const& ExpData::getObservedDataStdDev() const {
Expand Down Expand Up @@ -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<realtype> const& ExpData::getObservedEvents() const {
Expand Down Expand Up @@ -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;
}
Expand Down
9 changes: 5 additions & 4 deletions src/rdata.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <cmath>

namespace amici {
using std::isnan;

ReturnData::ReturnData(Solver const& solver, Model const& model)
: ReturnData(
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
);
Expand Down Expand Up @@ -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) {
Expand Down
34 changes: 17 additions & 17 deletions src/solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
#include "amici/amici.h"
#include "amici/exception.h"
#include "amici/model.h"
#include "amici/symbolic_functions.h"

#include <sundials/sundials_context.h>

Expand All @@ -12,6 +11,7 @@
#include <memory>

namespace amici {
using std::isnan;

/* Error handler passed to SUNDIALS. */
void wrapErrHandlerFn(
Expand Down Expand Up @@ -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_)
Expand Down Expand Up @@ -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));
Expand All @@ -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));
Expand Down Expand Up @@ -863,7 +863,7 @@ void Solver::setAbsoluteTolerance(double const atol) {
}

double Solver::getRelativeToleranceFSA() const {
return static_cast<double>(isNaN(rtol_fsa_) ? rtol_ : rtol_fsa_);
return static_cast<double>(isnan(rtol_fsa_) ? rtol_ : rtol_fsa_);
}

void Solver::setRelativeToleranceFSA(double const rtol) {
Expand All @@ -878,7 +878,7 @@ void Solver::setRelativeToleranceFSA(double const rtol) {
}

double Solver::getAbsoluteToleranceFSA() const {
return static_cast<double>(isNaN(atol_fsa_) ? atol_ : atol_fsa_);
return static_cast<double>(isnan(atol_fsa_) ? atol_ : atol_fsa_);
}

void Solver::setAbsoluteToleranceFSA(double const atol) {
Expand All @@ -893,7 +893,7 @@ void Solver::setAbsoluteToleranceFSA(double const atol) {
}

double Solver::getRelativeToleranceB() const {
return static_cast<double>(isNaN(rtolB_) ? rtol_ : rtolB_);
return static_cast<double>(isnan(rtolB_) ? rtol_ : rtolB_);
}

void Solver::setRelativeToleranceB(double const rtol) {
Expand All @@ -908,7 +908,7 @@ void Solver::setRelativeToleranceB(double const rtol) {
}

double Solver::getAbsoluteToleranceB() const {
return static_cast<double>(isNaN(atolB_) ? atol_ : atolB_);
return static_cast<double>(isnan(atolB_) ? atol_ : atolB_);
}

void Solver::setAbsoluteToleranceB(double const atol) {
Expand Down Expand Up @@ -971,7 +971,7 @@ void Solver::setSteadyStateToleranceFactor(double const factor) {

double Solver::getRelativeToleranceSteadyState() const {
return static_cast<double>(
isNaN(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_
isnan(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_
);
}

Expand All @@ -984,7 +984,7 @@ void Solver::setRelativeToleranceSteadyState(double const rtol) {

double Solver::getAbsoluteToleranceSteadyState() const {
return static_cast<double>(
isNaN(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_
isnan(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_
);
}

Expand All @@ -1010,7 +1010,7 @@ void Solver::setSteadyStateSensiToleranceFactor(

double Solver::getRelativeToleranceSteadyStateSensi() const {
return static_cast<double>(
isNaN(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_
isnan(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_
);
}

Expand All @@ -1023,7 +1023,7 @@ void Solver::setRelativeToleranceSteadyStateSensi(double const rtol) {

double Solver::getAbsoluteToleranceSteadyStateSensi() const {
return static_cast<double>(
isNaN(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_
isnan(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_
);
}

Expand Down
Loading
Loading