Skip to content

abs() rework and odeint compatibility headers #67

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

Open
abgandar wants to merge 7 commits into
base: master
Choose a base branch
from
Open

abs() rework and odeint compatibility headers #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8edc8e8
Add daceAbsolute function to take DA absolute value
Apr 22, 2025
4d10f57
expose daceAbsolute in C++ interface. Can't use DA abs() to avoid con…
Apr 22, 2025
d603fc3
add double absolute() overload extension
Apr 22, 2025
b5e0c2f
The great abs() removal and AlgebraicVector cleanup (part 1)
Apr 25, 2025
6530704
boost::odeint compatibility header
Apr 25, 2025
0a922d9
Add AlgebraicVector::norm
Apr 25, 2025
8f7a843
add compatibility shims for boost::odeint vector_space_algebra
Apr 26, 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
8 changes: 5 additions & 3 deletions core/dacecompat.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,14 +63,16 @@ void dacepek(const DACEDA *ina, const unsigned int jj[], double *cjj)
}

/*! Compute absolute value of a DA object.
Same as daceNorm(ina, 0).
\param[in] ina Pointer to DA object to take absolute value of
\param[out] anorm Pointer where to store the absolute value
\deprecated Has been replaced by daceAbsoluteValue()
\sa daceAbsoluteValue()
\deprecated Has been replaced by daceNorm() and daceAbsolute()
\sa daceNorm()
\sa daceAbsolute()
*/
void daceabs(const DACEDA *ina, double *anorm)
{
*anorm = daceAbsoluteValue(ina);
*anorm = daceNorm(ina, 0);
}

/*! Compute absolute value of a DA object.
Expand Down
22 changes: 19 additions & 3 deletions core/dacemath.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -328,7 +328,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb)
continue;

const double c = ia->cc*ckon;
if(!(fabs(c) <= DACECom_t.eps))
if(LIKELY(!(fabs(c) <= DACECom_t.eps)))
{
ib->cc = c;
ib->ii = ia->ii;
Expand All @@ -345,7 +345,7 @@ void daceMultiplyDouble(const DACEDA *ina, const double ckon, DACEDA *inb)
continue;

const double c = ia->cc*ckon;
if(!(fabs(c) <= DACECom_t.eps))
if(LIKELY(!(fabs(c) <= DACECom_t.eps)))
{
if(UNLIKELY(ib >= ibmax))
{
Expand Down Expand Up @@ -653,6 +653,22 @@ void daceIntegrate(const unsigned int iint, const DACEDA *ina, DACEDA *inc)
* DACE intrinsic function routines
*********************************************************************************/

/*! Absolute value of a DA object.
Returns either the DA or the negative of the DA based on the constant part.
\param[in] ina Pointer to the DA object to operate on
\param[out] inc Pointer to the DA object to store the result in
\note This routine is aliasing safe, i.e. inc can be the same as ina.
\sa daceAbsoluteValue
\sa daceNorm
*/
void daceAbsolute(const DACEDA *ina, DACEDA *inc)
{
daceCopy(ina, inc);
const double a0 = daceGetConstant(inc);
if(a0 < 0.0)
daceMultiplyDouble(inc, -1.0, inc);
}

/*! Truncate the constant part of a DA object to an integer.
\param[in] ina Pointer to the DA object to operate on
\param[out] inc Pointer to the DA object to store the result in
Expand Down Expand Up @@ -684,7 +700,7 @@ void daceRound(const DACEDA *ina, DACEDA *inc)
void daceModulo(const DACEDA *ina, const double p, DACEDA *inc)
{
daceCopy(ina, inc);
daceSetCoefficient0(inc, 0, fmod(daceGetConstant(inc),p));
daceSetCoefficient0(inc, 0, fmod(daceGetConstant(inc), p));
}

/*! Raise a DA object to the p-th power.
Expand Down
9 changes: 0 additions & 9 deletions core/dacenorm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,15 +42,6 @@
* DACE norm and norm estimation routines
*********************************************************************************/

/*! Compute the absolute value (maximum coefficient norm) of a DA object.
\param[in] ina Pointer to the DA object to take absolute value of
\return The absolute value of ina
*/
double daceAbsoluteValue(const DACEDA *ina)
{
return daceNorm(ina, 0);
}

/*! Compute a norm of a DA object.
\param[in] ina Pointer to the DA object to take norm of
\param[in] ityp Type of norm to compute.
Expand Down
2 changes: 1 addition & 1 deletion core/include/dace/dacebase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,6 +170,7 @@ DACE_API void daceIntegrate(const unsigned int iint, const DACEDA REF(ina), DACE
/********************************************************************************
* DACE intrinsic function routines
*********************************************************************************/
DACE_API void daceAbsolute(const DACEDA REF(ina), DACEDA REF(inc));
DACE_API void daceTruncate(const DACEDA REF(ina), DACEDA REF(inc));
DACE_API void daceRound(const DACEDA REF(ina), DACEDA REF(inc));
DACE_API void daceModulo(const DACEDA REF(ina), const double p, DACEDA REF(inc));
Expand Down Expand Up @@ -213,7 +214,6 @@ DACE_API void dacePsiFunction(const DACEDA REF(ina), const unsigned int n, DACED
/********************************************************************************
* DACE norm and norm estimation routines
*********************************************************************************/
DACE_API double daceAbsoluteValue(const DACEDA REF(ina));
DACE_API double daceNorm(const DACEDA REF(ina), const unsigned int ityp);
DACE_API void daceOrderedNorm(const DACEDA REF(ina), const unsigned int ivar, const unsigned int ityp, double onorm[]);
DACE_API void daceEstimate(const DACEDA REF(ina), const unsigned int ivar, const unsigned int ityp, double c[], double err[], const unsigned int nc);
Expand Down
61 changes: 30 additions & 31 deletions interfaces/cxx/AlgebraicVector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,13 @@
#include "dace/AlgebraicVector.h"
#include "dace/AlgebraicVector_t.h"

namespace DACE{
namespace DACE {

/***********************************************************************************
* Coefficient access routines
************************************************************************************/
#ifdef WITH_ALGEBRAICMATRIX
template<> AlgebraicMatrix<double> AlgebraicVector<DA>::linear() const{
template<> AlgebraicMatrix<double> AlgebraicVector<DA>::linear() const {
/*! Return the linear part of a AlgebraicVector<T>. NOT DEFINED FOR TYPES OTHER THAN DA.
\return A AlgebraicMatrix<double> of dimension size by nvar, where size is the
size of the AlgebraicVector<T> considered and nvar is the number of variables defined
Expand All @@ -57,13 +57,13 @@ template<> AlgebraicMatrix<double> AlgebraicVector<DA>::linear() const{
const int nvar = DA::getMaxVariables();

AlgebraicMatrix<double> out(size, nvar);
for(size_t i=0; i<size; i++){
for(size_t i=0; i<size; i++) {
out.setrow( i, (*this)[i].linear() );
}
return out;
}
#else
template<> std::vector< std::vector<double> > AlgebraicVector<DA>::linear() const{
template<> std::vector< std::vector<double> > AlgebraicVector<DA>::linear() const {
/*! Return the linear part of a AlgebraicVector<T>. NOT DEFINED FOR TYPES OTHER THAN DA.
\return A std::vector< std::vector<double> >, where each std::vector<double> contains
the linear part of the corresponding DA included in the original AlgebraicVector<T>.
Expand All @@ -74,15 +74,14 @@ template<> std::vector< std::vector<double> > AlgebraicVector<DA>::linear() cons
const size_t size = this->size();

std::vector< std::vector<double> > out(size);
for(size_t i=0; i<size; i++){
for(size_t i=0; i<size; i++) {
out[i] = (*this)[i].linear();
}
return out;
}
#endif /* WITH_ALGEBRAICMATRIX */

template<> AlgebraicVector<DA> AlgebraicVector<DA>::trim(const unsigned int min, const unsigned int max) const
{
template<> AlgebraicVector<DA> AlgebraicVector<DA>::trim(const unsigned int min, const unsigned int max) const {
/*! Returns an AlgebraicVector<DA> with all monomials of order less than min and greater than max removed (trimmed). The result is copied in a new AlgebraicVector<DA>.
\param[in] min minimum order to be preserved.
\param[in] max maximum order to be preserved.
Expand All @@ -102,7 +101,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::trim(const unsigned int min,
/***********************************************************************************
* Math routines
************************************************************************************/
template<> AlgebraicVector<DA> AlgebraicVector<DA>::deriv(const unsigned int p) const{
template<> AlgebraicVector<DA> AlgebraicVector<DA>::deriv(const unsigned int p) const {
/*! Compute the derivative of a AlgebraicVector<T> with respect to variable p.
The result is copied in a new AlgebraicVector<T>. NOT DEFINED FOR TYPES OTHER THAN DA.
\param[in] p variable with respect to which the derivative is calculated.
Expand All @@ -113,13 +112,13 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::deriv(const unsigned int p)
*/
const size_t size = this->size();
AlgebraicVector<DA> temp(size);
for(size_t i=0; i<size; i++){
for(size_t i=0; i<size; i++) {
temp[i] = (*this)[i].deriv(p);}

return temp;
}

template<> AlgebraicVector<DA> AlgebraicVector<DA>::integ(const unsigned int p) const{
template<> AlgebraicVector<DA> AlgebraicVector<DA>::integ(const unsigned int p) const {
/*! Compute the integral of a AlgebraicVector<T> with respect to variable p.
The result is copied in a new AlgebraicVector<T>. NOT DEFINED FOR TYPES OTHER THAN DA.
\param[in] p variable with respect to which the integral is calculated.
Expand All @@ -130,7 +129,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::integ(const unsigned int p)
*/
const size_t size = this->size();
AlgebraicVector<DA> temp(size);
for(size_t i=0; i<size; i++){
for(size_t i=0; i<size; i++) {
temp[i] = (*this)[i].integ(p);}

return temp;
Expand All @@ -139,7 +138,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::integ(const unsigned int p)
/***********************************************************************************
* Polynomial evaluation routines
************************************************************************************/
template<> compiledDA AlgebraicVector<DA>::compile() const{
template<> compiledDA AlgebraicVector<DA>::compile() const {
/*! Compile vector of polynomials and create a compiledDA object.
\return The compiled DA object.
\note This DA specific function is only available in AlgebraicVector<DA>.
Expand All @@ -149,7 +148,7 @@ template<> compiledDA AlgebraicVector<DA>::compile() const{
return compiledDA(*this);
}

template<> AlgebraicVector<DA> AlgebraicVector<DA>::plug(const unsigned int var, const double val) const{
template<> AlgebraicVector<DA> AlgebraicVector<DA>::plug(const unsigned int var, const double val) const {
/*! Partial evaluation of vector of polynomials. In each element of the vector,
variable var is replaced by the value val. The resulting vector of DAs
is returned.
Expand All @@ -162,7 +161,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::plug(const unsigned int var,
*/
const size_t size = this->size();
AlgebraicVector<DA> temp(size);
for(size_t i=0; i<size; i++){
for(size_t i=0; i<size; i++) {
temp[i] = (*this)[i].plug(var,val);}

return temp;
Expand All @@ -177,20 +176,20 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::plug(const unsigned int var,
This is NOT intended for public use. Limited error checking is performed
in accordance with the exclusive use of this routine in map inversion.
*/
template<> void AlgebraicVector<DA>::matrix_inverse(std::vector< std::vector<double> > &A){
template<> void AlgebraicVector<DA>::matrix_inverse(std::vector< std::vector<double> > &A) {
using std::abs;

const size_t n = A.size();
std::vector<size_t> indexc(n), indexr(n), ipiv(n, 0);

for (size_t i=0; i<n; i++){
for (size_t i=0; i<n; i++) {
size_t icol = 0, irow = 0;
double big = 0.0;
for (size_t j=0; j<n; j++)
if (ipiv[j] == 0)
for (size_t k=0; k<n; k++)
if (ipiv[k] == 0)
if (abs(A[j][k]) >= big){
if (abs(A[j][k]) >= big) {
big = abs(A[j][k]);
irow = j;
icol = k;}
Expand All @@ -204,7 +203,7 @@ template<> void AlgebraicVector<DA>::matrix_inverse(std::vector< std::vector<dou
A[icol][icol] = 1.0;
for (size_t l=0; l<n; l++) A[icol][l] *= pivinv;
for (size_t ll=0; ll<n; ll++)
if (ll != icol){
if (ll != icol) {
const double temp = A[ll][icol];
A[ll][icol] = 0.0;
for (size_t l=0; l<n; l++) A[ll][l] -= A[icol][l]*temp;}
Expand All @@ -218,7 +217,7 @@ template<> void AlgebraicVector<DA>::matrix_inverse(std::vector< std::vector<dou
/*! \endcond */
#endif /* WITH_ALGEBRAICMATRIX */

template<> AlgebraicVector<DA> AlgebraicVector<DA>::invert() const{
template<> AlgebraicVector<DA> AlgebraicVector<DA>::invert() const {
/*! Invert the polynomials map given by the AlgebraicVector<DA>.
\return the inverted polynomials
\throw std::runtime_error
Expand Down Expand Up @@ -258,21 +257,21 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::invert() const{
// Compute DA representation of the inverse of the linear part of the map and its composition with non-linear part AN
AlgebraicVector<DA> Linv(nvar);
// Linv = AI*AN
for(size_t i=0; i<nvar; i++){
for(size_t i=0; i<nvar; i++) {
Linv[i] = 0.0;
for(size_t j=0; j<nvar; j++)
Linv[i] += AI[i][j]*AN[j];}
compiledDA AIoAN(Linv);
// Linv = AI*DDA
for(size_t i=0; i<nvar; i++){
for(size_t i=0; i<nvar; i++) {
Linv[i] = 0.0;
for(size_t j=0; j<nvar; j++)
Linv[i] += AI[i][j]*DDA[j];}
#endif /* WITH_ALGEBRAICMATRIX */

// Iterate to obtain the inverse map
AlgebraicVector<DA> MI = Linv;
for(unsigned int i=1; i<ord; i++){
for(unsigned int i=1; i<ord; i++) {
DA::setTO(i+1);
MI = Linv - AIoAN.eval(MI);}

Expand All @@ -282,7 +281,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::invert() const{
/********************************************************************************
* Static factory routines
*********************************************************************************/
template<> AlgebraicVector<DA> AlgebraicVector<DA>::identity(const size_t n){
template<> AlgebraicVector<DA> AlgebraicVector<DA>::identity(const size_t n) {
/*! Return the DA identity of dimension n.
\param[in] n The dimendion of the identity.
\return AlgebraicVector<DA> containing the DA identity in n dimensions
Expand All @@ -291,7 +290,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::identity(const size_t n){
error will be the result.
*/
AlgebraicVector<DA> temp(n);
for(size_t i=0; i < n; i++){
for(size_t i=0; i < n; i++) {
temp[i] = DA((int)(i+1));}

return temp;
Expand All @@ -301,7 +300,7 @@ template<> AlgebraicVector<DA> AlgebraicVector<DA>::identity(const size_t n){
* Non-member functions
************************************************************************************/
#ifdef WITH_ALGEBRAICMATRIX
template<> AlgebraicMatrix<double> linear(const AlgebraicVector<DA> &obj){
template<> AlgebraicMatrix<double> linear(const AlgebraicVector<DA> &obj) {
/*! Return the linear part of a AlgebraicVector<T>.
\param[in] obj AlgebraicVector<T> to extract linear part from
\return An AlgebraicMatrix<double> of dimensions size by nvar, where
Expand All @@ -313,7 +312,7 @@ template<> AlgebraicMatrix<double> linear(const AlgebraicVector<DA> &obj){
\sa AlgebraicVector<T>::linear
*/
#else
template<> std::vector< std::vector<double> > linear(const AlgebraicVector<DA> &obj){
template<> std::vector< std::vector<double> > linear(const AlgebraicVector<DA> &obj) {
/*! Return the linear part of a AlgebraicVector<T>. Only defined for AlgebraicVector<DA>.
\param[in] obj AlgebraicVector<T> to extract linear part from
\return A std::vector< std::vector<double> > containing the linear parts of
Expand All @@ -327,7 +326,7 @@ template<> std::vector< std::vector<double> > linear(const AlgebraicVector<DA> &
return obj.linear();
}

template<> AlgebraicVector<DA> trim(const AlgebraicVector<DA> &obj, unsigned int min, unsigned int max){
template<> AlgebraicVector<DA> trim(const AlgebraicVector<DA> &obj, unsigned int min, unsigned int max) {
/*! Returns an AlgebraicVector<DA> with all monomials of order less than min and greater than max removed (trimmed). The result is copied in a new AlgebraicVector<DA>.
\param[in] obj the AlgebraicVector<DA> to be trimmed.
\param[in] min minimum order to be preserved.
Expand All @@ -341,7 +340,7 @@ template<> AlgebraicVector<DA> trim(const AlgebraicVector<DA> &obj, unsigned int
return obj.trim(min, max);
}

template<> AlgebraicVector<DA> deriv(const AlgebraicVector<DA> &obj, const unsigned int p){
template<> AlgebraicVector<DA> deriv(const AlgebraicVector<DA> &obj, const unsigned int p) {
/*! Compute the derivative of a AlgebraicVector<T> with respect to variable p.
The result is copied in a new AlgebraicVector<T>.
\param[in] obj AlgebraicVector<T>.
Expand All @@ -355,7 +354,7 @@ template<> AlgebraicVector<DA> deriv(const AlgebraicVector<DA> &obj, const unsig
return obj.deriv(p);
}

template<> AlgebraicVector<DA> integ(const AlgebraicVector<DA> &obj, const unsigned int p){
template<> AlgebraicVector<DA> integ(const AlgebraicVector<DA> &obj, const unsigned int p) {
/*! Compute the integral of a AlgebraicVector<T> with respect to variable p.
The result is copied in a new AlgebraicVector<T>.
\param[in] obj AlgebraicVector<T>.
Expand All @@ -369,7 +368,7 @@ template<> AlgebraicVector<DA> integ(const AlgebraicVector<DA> &obj, const unsig
return obj.integ(p);
}

template<> compiledDA compile(const AlgebraicVector<DA> &obj){
template<> compiledDA compile(const AlgebraicVector<DA> &obj) {
/*! Compile vector of polynomials and create a compiledDA object.
\param[in] obj The AlgebraicVector to compile.
\return The compiled DA object.
Expand All @@ -380,7 +379,7 @@ template<> compiledDA compile(const AlgebraicVector<DA> &obj){
return obj.compile();
}

template<> AlgebraicVector<DA> plug(const AlgebraicVector<DA> &obj, const unsigned int var, const double val){
template<> AlgebraicVector<DA> plug(const AlgebraicVector<DA> &obj, const unsigned int var, const double val) {
/*! Partial evaluation of vector of polynomials. In each element of the vector,
variable var is replaced by the value val. The resulting vector of DAs
is returned.
Expand Down
Loading