Skip to content

Add daceIsNan and daceIsInf functions and C++ bindings #66

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
abgandar merged 4 commits into from
Apr 21, 2025
+118 −43
Merged

Add daceIsNan and daceIsInf functions and C++ bindings #66

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2922c85
Add some LIKELY keywords to basic loops, whitespace
Apr 21, 2025
72f015c
isNan and isInf functions
Apr 21, 2025
3e50d2c
white space
Apr 21, 2025
570b72e
expose isnan and isinf in C++ interface
Apr 21, 2025
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
38 changes: 37 additions & 1 deletion core/dacebasic.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -569,7 +569,7 @@ void daceTrim(const DACEDA *ina, const unsigned int imin, const unsigned int ima

/*! Copy monomials from a DA object ina to DA object inb if the same monomial
is non-zero in DA object inc, while filtering out terms below the current
cutoff
cutoff.
\param[in] ina Pointer to DA object to filter
\param[in] inb Pointer to DA object to store the filtered result in
\param[in] inc Pointer to DA object providing the filter template
Expand Down Expand Up @@ -627,4 +627,40 @@ void daceFilter(const DACEDA *ina, DACEDA *inb, const DACEDA *inc)

daceSetLength(inb, ib-ipob);
}

/*! Check each coefficient of DA object ina to see if any of them are NANs (not a number).
\param[in] ina Pointer to DA object to check
\return True (non-zero) if any of the coefficients of ina is NAN
*/
unsigned int daceIsNan(const DACEDA *ina)
{
monomial *ipoa; unsigned int ilma, illa;

daceVariableInformation(ina, &ipoa, &ilma, &illa);

for(monomial *ia = ipoa; ia < ipoa+illa; ia++)
{
if(isnan(ia->cc))
return true;
}
return false;
}

/*! Check each coefficient of DA object ina to see if any of them are INF (infinity).
\param[in] ina Pointer to DA object to check
\return True (non-zero) if any of the coefficients of ina is INF
*/
unsigned int daceIsInf(const DACEDA *ina)
{
monomial *ipoa; unsigned int ilma, illa;

daceVariableInformation(ina, &ipoa, &ilma, &illa);

for(monomial *ia = ipoa; ia < ipoa+illa; ia++)
{
if(isinf(ia->cc))
return true;
}
return false;
}
/** @}*/
8 changes: 4 additions & 4 deletions core/daceeval.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,12 +61,12 @@ double daceEvalMonomials(const DACEDA *ina, const DACEDA *inb)
monomial *const ibmax = ipob + ilmb;

double res = 0.0;
for (monomial *i = ipoa; i < ipoa + illa; i++)
for(monomial *i = ipoa; i < ipoa + illa; i++)
{
while (ib->ii < i->ii && ib < ibmax)
while(ib->ii < i->ii && ib < ibmax)
ib++;
if (ib == ibmax) break;
if (ib->ii == i->ii) res += ib->cc*i->cc;
if(ib == ibmax) break;
if(ib->ii == i->ii) res += ib->cc*i->cc;
}

return res;
Expand Down
66 changes: 33 additions & 33 deletions core/dacemath.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,14 +222,14 @@ void daceMultiplyMonomials(const DACEDA *ina, const DACEDA *inb, DACEDA *inc)
monomial *ib = ipob, *ic = ipoc;
monomial *const ibmax = ipob + ilmb, *const icmax = ipoc + ilmc;

for (monomial *i = ipoa; i < ipoa + illa; i++)
for(monomial *i = ipoa; i < ipoa + illa; i++)
{
while (ib->ii < i->ii && ib < ibmax)
while(ib->ii < i->ii && ib < ibmax)
ib++;
if (ib == ibmax) break;
if (ib->ii == i->ii)
if(ib == ibmax) break;
if(ib->ii == i->ii)
{
if (ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand Down Expand Up @@ -324,7 +324,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb)
{
for(monomial *ia = ipoa; ia < ipoa+illa; ia++)
{
if(DACECom.ieo[ia->ii] > DACECom_t.nocut)
if(UNLIKELY(DACECom.ieo[ia->ii] > DACECom_t.nocut))
continue;

const double c = ia->cc*ckon;
Expand All @@ -341,13 +341,13 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb)
monomial *const ibmax = ipob+ilmb;
for(monomial *ia = ipoa; ia < ipoa+illa; ia++)
{
if(DACECom.ieo[ia->ii] > DACECom_t.nocut)
if(UNLIKELY(DACECom.ieo[ia->ii] > DACECom_t.nocut))
continue;

const double c = ia->cc*ckon;
if(!(fabs(c) <= DACECom_t.eps))
{
if(ib >= ibmax)
if(UNLIKELY(ib >= ibmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand Down Expand Up @@ -451,13 +451,13 @@ void daceDivideByVariable(const DACEDA *ina, const unsigned int var, const unsig
const unsigned int ic1 = DACECom.ie1[i->ii];
const unsigned int ic2 = DACECom.ie2[i->ii];
const unsigned int ipow = (ic2/idiv)%ibase;
if(ipow < p)
if(UNLIKELY(ipow < p))
{
daceSetError(__func__, DACE_ERROR, 42);
daceCreateConstant(inc, 0.0);
return;
}
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand All @@ -474,13 +474,13 @@ void daceDivideByVariable(const DACEDA *ina, const unsigned int var, const unsig
const unsigned int ic1 = DACECom.ie1[i->ii];
const unsigned int ic2 = DACECom.ie2[i->ii];
const unsigned int ipow = (ic1/idiv)%ibase;
if(ipow < p)
if(UNLIKELY(ipow < p))
{
daceSetError(__func__, DACE_ERROR, 42);
daceCreateConstant(inc, 0.0);
return;
}
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand Down Expand Up @@ -534,7 +534,7 @@ void daceDifferentiate(const unsigned int idif, const DACEDA *ina, DACEDA *inc)
const unsigned int ipow = (ic2/idiv)%ibase;
if(ipow == 0 || DACECom.ieo[i->ii] > DACECom_t.nocut+1)
continue;
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand All @@ -553,7 +553,7 @@ void daceDifferentiate(const unsigned int idif, const DACEDA *ina, DACEDA *inc)
const unsigned int ipow = (ic1/idiv)%ibase;
if(ipow == 0 || DACECom.ieo[i->ii] > DACECom_t.nocut+1)
continue;
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand Down Expand Up @@ -602,15 +602,15 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc)
{
for(monomial *i = ipoa; i < ipoa+illa; i++)
{
if(DACECom.ieo[i->ii] >= DACECom_t.nocut)
if(UNLIKELY(DACECom.ieo[i->ii] >= DACECom_t.nocut))
continue;
const unsigned int ic1 = DACECom.ie1[i->ii];
const unsigned int ic2 = DACECom.ie2[i->ii];
const unsigned int ipow = (ic2/idiv)%ibase;
const double ccc = i->cc/(ipow+1);
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand All @@ -625,15 +625,15 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc)
{
for(monomial *i = ipoa; i < ipoa+illa; i++)
{
if(DACECom.ieo[i->ii] >= DACECom_t.nocut)
if(UNLIKELY(DACECom.ieo[i->ii] >= DACECom_t.nocut))
continue;
const unsigned int ic1 = DACECom.ie1[i->ii];
const unsigned int ic2 = DACECom.ie2[i->ii];
const unsigned int ipow = (ic1/idiv)%ibase;
const double ccc = i->cc/(ipow+1);
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
break;
Expand Down Expand Up @@ -2047,7 +2047,7 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co
monomial *ia = ipoa, *ib = ipob, *ic = ipoc;
monomial *const iamax = ipoa+illa, *const ibmax = ipob+illb, *const icmax = ipoc+ilmc;

if(illa > 0 && illb > 0)
if(LIKELY(illa > 0 && illb > 0))
{
// both polynomials have coefficients, merge until one runs out
unsigned int ja = ia->ii;
Expand All @@ -2057,12 +2057,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co
if(ja == jb)
{
// add the two terms
if(DACECom.ieo[ja] <= DACECom_t.nocut)
if(LIKELY(DACECom.ieo[ja] <= DACECom_t.nocut))
{
const double ccc = ia->cc*afac + ib->cc*bfac;
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
daceSetLength(inc, ilmc);
Expand All @@ -2081,12 +2081,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co
else if(ja < jb)
{
// store term a
if(DACECom.ieo[ja] <= DACECom_t.nocut)
if(LIKELY(DACECom.ieo[ja] <= DACECom_t.nocut))
{
const double ccc = ia->cc*afac;
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
daceSetLength(inc, ilmc);
Expand All @@ -2104,12 +2104,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co
else
{
// store term b
if(DACECom.ieo[jb] <= DACECom_t.nocut)
if(LIKELY(DACECom.ieo[jb] <= DACECom_t.nocut))
{
const double ccc = ib->cc*bfac;
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
daceSetLength(inc, ilmc);
Expand Down Expand Up @@ -2145,12 +2145,12 @@ void daceWeightedSum(const DACEDA *ina, const double afac, const DACEDA *inb, co

for(monomial *is = ismin; is < ismax; is++)
{
if(DACECom.ieo[is->ii] <= DACECom_t.nocut)
if(LIKELY(DACECom.ieo[is->ii] <= DACECom_t.nocut))
{
const double ccc = is->cc*fac;
if(!(fabs(ccc) <= DACECom_t.eps))
if(LIKELY(!(fabs(ccc) <= DACECom_t.eps)))
{
if(ic >= icmax)
if(UNLIKELY(ic >= icmax))
{
daceSetError(__func__, DACE_ERROR, 21);
daceSetLength(inc, ilmc);
Expand Down
2 changes: 2 additions & 0 deletions core/include/dace/dacebase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,8 @@ DACE_API void daceCopy(const DACEDA REF(ina), DACEDA REF(inb));
DACE_API void daceCopyFiltering(const DACEDA REF(ina), DACEDA REF(inb));
DACE_API void daceFilter(const DACEDA REF(ina), DACEDA REF(inb), const DACEDA REF(inc));
DACE_API void daceTrim(const DACEDA REF(ina), const unsigned int imin, const unsigned int imax, DACEDA REF(inc));
DACE_API unsigned int daceIsNan(const DACEDA REF(ina));
DACE_API unsigned int daceIsInf(const DACEDA REF(ina));

/********************************************************************************
* Basic DACE arithmetic operations
Expand Down
42 changes: 38 additions & 4 deletions interfaces/cxx/DA.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,6 +317,28 @@ DA::~DA() throw(){
/********************************************************************************
* Coefficient access routines
*********************************************************************************/
int DA::isnan() const{
/*! Check if a DA object has any NAN coefficients.
\return True is any coefficients of the DA object are NAN.
\throw DACE::DACEException
*/
const int temp = daceIsNan(m_index);
if(daceGetError()) DACEException();

return temp;
}

int DA::isinf() const{
/*! Check if a DA object has any INF coefficients.
\return True is any coefficients of the DA object are INF.
\throw DACE::DACEException
*/
const int temp = daceIsInf(m_index);
if(daceGetError()) DACEException();

return temp;
}

double DA::cons() const{
/*! Return the constant part of a DA object.
\return A double corresponding to the constant part of the DA object.
Expand Down Expand Up @@ -1826,9 +1848,6 @@ std::istream& operator>>(std::istream &in, DA &da){
// convert string vector to DA
da = DA::fromString(strs);
}



}

return in;
Expand Down Expand Up @@ -1932,6 +1951,22 @@ void DA::memdump(){
/********************************************************************************
* DACE non-member functions
*********************************************************************************/
int isnan(const DA &da) {
/*! Check if a DA object has any NAN coefficients.
\param[in] da a given DA object.
\return True if any coefficients of the DA object are NAN.
\throw DACE::DACEException
*/
return da.isnan();}

int isinf(const DA &da) {
/*! Check if a DA object has any INF coefficients.
\param[in] da a given DA object.
\return True if any coefficients of the DA object are INF.
\throw DACE::DACEException
*/
return da.isinf();}

double cons(const DA &da) {
/*! Return the constant part of a DA object.
\param[in] da a given DA object.
Expand Down Expand Up @@ -1959,7 +1994,6 @@ AlgebraicVector<DA> gradient(const DA &da) {

return da.gradient();}


DA divide(const DA &da, const unsigned int var, const unsigned int p){
/*! Divide by independent variable var raised to power p.
The result is copied in a new DA object.
Expand Down
5 changes: 4 additions & 1 deletion interfaces/cxx/include/dace/DA.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,6 @@ class DACE_API DA
private:
static bool initialized; //!< Indicates if DA::init() was called
static std::stack<unsigned int> TOstack; //!< Truncation order stack
// XXX: Mauro, MSVC is bitching around (also for other templated classes, e.g. std::string): Warning C4251 'DACE::DA::TOstack': class 'std::stack<unsigned int,std::deque<_Ty,std::allocator<_Ty>>>' needs to have dll - interface to be used by clients of class 'DACE::DA'
DACEDA m_index; //!< Index to the DA vector

public:
Expand Down Expand Up @@ -100,6 +99,8 @@ class DACE_API DA
/********************************************************************************
* Coefficient access and extraction routines
*********************************************************************************/
int isnan() const;
int isinf() const;
double cons() const; //!< Get constant part of a DA
AlgebraicVector<double> linear() const; //!< Get linear part of a DA
AlgebraicVector<DA> gradient() const; //!< Gradient vector with respect to all independent DA variables
Expand Down Expand Up @@ -259,6 +260,8 @@ class DACE_API DA
/********************************************************************************
* DACE non-member functions
*********************************************************************************/
DACE_API int isnan(const DA &da);
DACE_API int isinf(const DA &da);
DACE_API double cons(const DA &da);
DACE_API AlgebraicVector<double> linear(const DA &da);
DACE_API AlgebraicVector<DA> gradient(const DA &da);
Expand Down