diffNumber.h

Go to the documentation of this file.

#ifndef __DIFFNUMBER_H__
#define __DIFFNUMBER_H__

#include <stdlib.h>
#include <iostream>

#include <coconut_vect.h>

const int STANDARD_MAXDEG_DIFFNUMBER=4;

/****************/
/*DIFFNUMBER*/
/****************/

//CLASSDEF

class diffNumber
{
        private:
                double *data;
        public:
                int MAXDEG;
                inline diffNumber(void);
                inline diffNumber(double a);
                inline diffNumber(int n);
                inline diffNumber(double a, int n);
                inline diffNumber(const diffNumber &a);
                inline virtual ~diffNumber();
        
                inline diffNumber& operator=(const diffNumber &a);
                inline double& operator[] (int pos) const;
        
                inline static diffNumber getVar(const double &x);
                inline static diffNumber getVar(const double &x, int n);
                inline void setDeg(int n);
};
//OPERATORS
inline diffNumber operator +(const diffNumber &a, const diffNumber &b);
inline diffNumber operator +(const double &a, const diffNumber &b);
inline diffNumber operator +(const diffNumber &a, const double &b);
inline diffNumber operator -(const diffNumber &a, const diffNumber &b);
inline diffNumber operator -(const double &a, const diffNumber &b);
inline diffNumber operator -(const diffNumber &a, const double &b);
inline diffNumber operator -(const diffNumber &a);
inline diffNumber operator *(const diffNumber &a, const diffNumber &b);
inline diffNumber operator *(const double &a, const diffNumber &b);
inline diffNumber operator *(const diffNumber &a, const double &b);
inline diffNumber operator /(const diffNumber &a, const diffNumber &b);
inline diffNumber operator /(const diffNumber &a, const double &b);
inline diffNumber operator /(const double &a, const diffNumber &b);

//OUTPUT
inline std::ostream &operator << (std::ostream &s, const diffNumber &a);

//FUNCTIONS
inline diffNumber exp(const diffNumber &a);
inline diffNumber sin(const diffNumber &a);
inline diffNumber cos(const diffNumber &a);
inline diffNumber sqrt(const diffNumber &a);

inline diffNumber power(const diffNumber &a, int n);
inline diffNumber power(const diffNumber &a, const diffNumber &b);// added by Simon Haller-> needed for ForwardHessEvaluator

inline diffNumber log(const diffNumber &a);
inline diffNumber sinh(const diffNumber &a);
inline diffNumber cosh(const diffNumber &a);
inline diffNumber atan2(const diffNumber &a, const diffNumber &b);
inline diffNumber atan(const diffNumber &a);
inline diffNumber square(const diffNumber &a);

/****************/
/*INLINECODE*/
/****************/

/*diffNumber*/

//CONSTRUCTORS
inline diffNumber::diffNumber(void)
{
        // std::cout << "standard" << std::endl;
        MAXDEG=STANDARD_MAXDEG_DIFFNUMBER;
        data=new double[MAXDEG+1];
        for (int i=0; i<MAXDEG+1; i++)
                data[i]=0.;
}

inline diffNumber::diffNumber(int n)
{
        MAXDEG=n;
        data=new double[MAXDEG+1];
        for (int i=0; i<MAXDEG+1; i++)
                data[i]=0.;
}

inline diffNumber::diffNumber(double a)
{
        MAXDEG=STANDARD_MAXDEG_DIFFNUMBER;
        data= new double[MAXDEG+1];
        data[0]=a;
        for (int i=1; i<MAXDEG+1; i++)
                data[i]=0;
}

inline diffNumber::diffNumber(double a, int n)
{
        MAXDEG=n;
        data= new double[MAXDEG+1];
        data[0]=a;
        for (int i=1; i<MAXDEG+1; i++)
                data[i]=0;
}

inline diffNumber::diffNumber(const diffNumber &d)
{
        MAXDEG = d.MAXDEG;
        data = new double[MAXDEG+1];
        memcpy(data, d.data, (MAXDEG+1)*sizeof(double));
}


inline diffNumber& diffNumber::operator=(const diffNumber &a)
{
        if ((MAXDEG!=-1) && (MAXDEG != a.MAXDEG)){
                std::cerr<<"= degree not the same: "<< MAXDEG <<" != "<< a.MAXDEG;
                abort();
        }
        
        if (this != &a){
                
                if (data) delete [] data;
                
                data = new double[a.MAXDEG+1];
                for (int i=0; i<a.MAXDEG+1; i++)
                        data[i]=a[i];
        }
        return (*this);
}


inline double& diffNumber::operator[] (int pos) const
{
//      if (pos<
//      cout<<"diffNumber []";
        return data[pos];
}



//VARIABLES
inline diffNumber diffNumber::getVar(const double &x)
{
        diffNumber c;
        c[0]=x;
        c[1]=1;
        return c;
}

inline diffNumber diffNumber::getVar(const double &x, int n)
{
        diffNumber c(n);
        c[0]=x;
        c[1]=1;
        return c;
}

inline void diffNumber::setDeg(int n)
{
        __vectgen(double,tmp,n+1);
        if (MAXDEG>n){
                for (int i=0; i<n+1; i++)
                        tmp[i]=data[i];
                delete [] data;
                data=new double[n+1];
                for (int i=0; i<n+1; i++)
                        data[i]=tmp[i];
        }
        else {
                for (int i=0; i<MAXDEG+1; i++)
                        tmp[i]=data[i];
                delete [] data;
                data=new double[n+1];
                for (int i=0; i<MAXDEG+1; i++)
                        data[i]=tmp[i];
                for (int j=MAXDEG+1; j<n+1; j++)
                        data[j]=0.;
        }
        MAXDEG=n;
}

//OPERATORS
inline diffNumber operator +(const diffNumber &a, const diffNumber &b)
{
        if (a.MAXDEG != b.MAXDEG){
                std::cerr<<"+ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
                abort();
        }
        diffNumber c(a.MAXDEG);
        for (int i=0; i<a.MAXDEG+1; i++){
                c[i]=a[i] + b[i];
        }
        return c;
}

inline diffNumber operator +(const double &a, const diffNumber &b)
{
        diffNumber c(b.MAXDEG);
        c[0]=a+b[0];
        for (int i=1; i<b.MAXDEG+1; i++){
                c[i]=b[i];
        }
        return c;
}

inline diffNumber operator +(const diffNumber &a, const double &b)
{
        diffNumber c(a.MAXDEG);
        c[0]=a[0]+b;
        for (int i=1; i<a.MAXDEG+1; i++){
                c[i]=a[i];
        }
        return c;
}

inline diffNumber operator -(const diffNumber &a, const diffNumber &b)
{
        if (a.MAXDEG != b.MAXDEG){
                std::cerr<<"- degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
                abort();
        }
        diffNumber c(a.MAXDEG);
        for (int i=0; i<a.MAXDEG+1; i++){
                c[i]=a[i] - b[i];
        }
        return c;
}

inline diffNumber operator -(const double &a, const diffNumber &b)
{
        diffNumber c(b.MAXDEG);
        c[0]=a-b[0];
        for (int i=1; i<b.MAXDEG+1; i++){
                c[i]=-b[i];
        }
        return c;
}

inline diffNumber operator -(const diffNumber &a, const double &b)
{
        diffNumber c(a.MAXDEG);
        c[0]=a[0]-b;
        for (int i=1; i<a.MAXDEG+1; i++){
                c[i]=a[i];
        }
        return c;
}

inline diffNumber operator -(const diffNumber &a)
{
        for (int i=0; i<a.MAXDEG+1; i++){
                a[i]=-a[i];
        }
        return a;
}

inline diffNumber operator *(const diffNumber &a, const diffNumber &b)
{
        if (a.MAXDEG != b.MAXDEG){
                std::cerr<<"* degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
                abort();
        }
        diffNumber c(a.MAXDEG); //int constant;
        for (int i=0; i<a.MAXDEG+1; i++){
                //constant =1; a.MAXDEG+1
                for (int j=0; j<=i; j++){
                        c[i]=c[i] +  a[j] * b[i-j]; // error, constant missing. only valid for 0th and 1st order diffNumber! 
                }
        }
        return c;
}

inline diffNumber operator *(const double &a, const diffNumber &b)
{
        diffNumber c(b.MAXDEG);
        for (int i=0; i<b.MAXDEG+1; i++){
                c[i]=a * b[i];
        }
        return c;
}

inline diffNumber operator *(const diffNumber &a, const double &b)
{
        diffNumber c(a.MAXDEG);
        for (int i=0; i<a.MAXDEG+1; i++){
                c[i]=b * a[i];
        }
        return c;
}

inline diffNumber operator /(const diffNumber &a, const diffNumber &b)
{
        if (a.MAXDEG != b.MAXDEG){
                std::cerr<<"/ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
                abort();
        }
        diffNumber c(a.MAXDEG);
        c[0]=a[0] / b[0];
        for (int i=1; i<a.MAXDEG+1; i++){
                c[i]=0.;
                for (int j=0; j<i; j++){
                        c[i] += c[j] * b[i-j];
                }
                c[i]=1. / b[0] * (a[i] - c[i]);
        }
        return c;
}

inline diffNumber operator /(const double &a, const diffNumber &b)
{
        diffNumber c(b.MAXDEG);
        c[0]=a / b[0];
        for (int i=1; i<b.MAXDEG+1; i++){
                c[i]=0.;
                for (int j=0; j<i; j++){
                        c[i] += c[j] * b[i-j];
                }
                c[i]= - 1. / b[0] * c[i];
        }
        return c;
}

inline diffNumber operator /(const diffNumber &a, const double &b)
{
        diffNumber c(a.MAXDEG);
        for (int i=0; i<a.MAXDEG+1; i++){
                c[i]=a[i] / b;
        }
        return c;
}

inline diffNumber exp(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        c[0]=std::exp(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i; j++){
                        c[i]=c[i] + j* a[j] * c[i-j] / i;
                }
        }
        return c;
}

inline diffNumber sin(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        diffNumber d(a.MAXDEG);
        c[0]=std::sin(a[0]);
        d[0]=std::cos(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i; j++){
                        c[i]=c[i] + j* a[j] * d[i-j] / i;
                        d[i]=d[i] - j* a[j] * c[i-j] / i;
                }
        }
        return c;
}

inline diffNumber cos(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        diffNumber d(a.MAXDEG);
        c[0]=std::sin(a[0]);
        d[0]=std::cos(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i; j++){
                        c[i]=c[i] + j* a[j] * d[i-j] / i;
                        d[i]=d[i] - j* a[j] * c[i-j] / i;
                }
        }
        return d;
}

inline diffNumber sqrt(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        c[0]=std::sqrt(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                c[i]=a[i];
                int m = (i + (i%2))/2;
                for (int j=1; j<=m-1; j++){
                        c[i]=c[i] - 2.* c[j] * c[i-j];
                }
                if (i%2==0) {
                        int n=i/2;
                        c[i]=c[i]-c[n]*c[n];
                }
                c[i]=c[i] * 0.5 / c[0];
        }
        return c;
}

inline diffNumber power(const diffNumber &a, int n)
{
        diffNumber c(a.MAXDEG);
        
        if (a[0] == 0)
                for (int i=0; i<a.MAXDEG+1; i++) c[i] = 0;
        else
        {
                c[0]=std::pow(a[0],n);
        
                for (int i=1; i<a.MAXDEG+1; i++) {
                        for (int j=0; j<=i-1; j++){
                                 c[i] = c[i] + ( n * (i - j) - j)* c[j] * a[i-j] / i;
                                //c[i] = c[i] + (n * (i - j) - j)* c[j] * a[i-j] / i;
                        }
                        c[i]=c[i] / a[0];
                }
        }
        return c;       
}

inline diffNumber power(const diffNumber &a, const diffNumber &b)
{       
        if (a.MAXDEG != b.MAXDEG){
                std::cerr<<"+ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
                abort();
        }
        diffNumber c(a.MAXDEG);
        
        if (a[0] == 0.0)
                for (int i=0; i<a.MAXDEG+1; i++) c[i] = 0;
        else
                c= exp(b * log(a)); // lazyman's approach of implementing the power-function...
        
        return c;
} 


/*inline diffNumber power(const diffNumber &a, int n)
{       
        diffNumber c(a.MAXDEG); 
        c= exp(n * log(a));
        return c;
}*/

/* inline diffNumber cpower(const diffNumber &a, int n)
{
        diffNumber c(a.MAXDEG);
        c[0]=std::pow(a[0],n);
        int const1,const2;
        const1=const2=0;
        //double temp;
        //std::cout << "MAXDEG="<< a.MAXDEG << std::endl;
        switch(a.MAXDEG) {
                case 0:
                {
                        //std::cout << "h[0]=pow(" << a[0] <<", "<< n << ");" << std::endl;
                        break;
                }
                case 1:
                {
                        //std::cout << "h[0]=pow(" << a[0] <<", "<< n << ");" << std::endl;
                        c[1] = n * c[0] * a[1] / a[0]; 
                        //std::cout << "h[1]=" <<n <<" * h[0] * f[1] / f[0];" << std::endl;
                        
                        
                        break;
                }
                case 2:
                {
                        //std::cout << "h[0]=pow(" << a[0] <<", "<< n << ");" << std::endl;
                        c[1] = n * c[0] * a[1] / a[0]; 
                        //std::cout << "h[1]=" <<n <<" * h[0] * f[1] / f[0];" << std::endl;
                        c[2] = ( n * c[0] * a[2] + (n-1) * c[1] * a[1] ) / a[0]; 
                        //std::cout << "h[2]=(" <<n <<" * h[0] * f[2] + " << n-1 << "* h[1] * f[1] ) / f[0]"<<   std::endl;
                        
                        break;
                }
                case 3:
                {
                        //std::cout << "h[0]=pow(" << a[0] <<", "<< n << ");" << std::endl;
                        c[1] = n * c[0] * a[1] / a[0]; 
                        //std::cout << "h[1]=n * h[0] * f[1] / f[0];" << std::endl;
                        c[2] = ( n * c[0] * a[2]  + (n-1) * c[1] * a[1]   ) / a[0]; 
                        //std::cout << "h[2]=( n * h[0] * f[2] + n-1 * h[1] * f[1] ) / f[0]"<<   std::endl;
                        c[3] = ( n * c[0] * a[3] +  (2*n-1) * c[1] * a[2] + (n-2) * c[2] * a[1]  ) / a[0]; 
                        //std::cout << "h[3]=(n * h[0] * f[3] + (2 * n - 1) * h[1] * f[2] + (n-1) * c[2] * f[1] ) / f[0]"<<   std::endl;
                        break;
                }
                default:
                {
                        //std::cout << "h[0]=pow(" << a[0] <<", "<< n << ");" << std::endl;
                        
                        /* c[1] = n * a[1]; ///  (1, n * a[1], 0,0,0)
                        c[0] *= a[0];    ///  (a[0], n * a[1], 0,0,0)
                        
                        //c[1] *= a[0];
                        
                        //std::cout << "h[1]=" <<n <<" * h[0] * f[1] / f[0];" << std::endl;
                        c[2] = ( n * c[0] * a[2]  + (n-1) * c[1] * a[1]   ); /// (a[0], n * a[1], n * a[0] * a[2] + (n-1) * n * a[1]* a[1], 0, 0)
                        c[0] *= a[0]; /// (a[0]², n * a[1] * a[0], n * a[0] * a[2] + (n-1) * n * a[1]* a[1], 0, 0)
                        c[1] *= a[0]; /// (a[0]², n * a[1] * a[0], n * a[0] * a[2] + (n-1) * n * a[1]* a[1], 0, 0)
                        
                        //std::cout << "h[2]=( n * h[0] * f[2] + (n - 1) * h[1] * f[1] ) / f[0]"<<   std::endl;
                        c[3] = ( n * c[0] * a[3] +  (2*n-1) * c[1] * a[2] + (n-2) * c[2] * a[1]  ); 
                        c[2] *= a[0]; /// (a[0]², n * a[1] * a[0], n * a[0] * a[2] + (n-1) * n * a[1]* a[1], 0, 0)
                        c[1] *= a[0];
                        c[0] *= a[0]; /// (a[0]², n * a[1] * a[0], n * a[0] * a[2] + (n-1) * n * a[1]* a[1], 0, 0)
                
                        temp=c[0]; */
                        /*
                        int counter = a.MAXDEG;
                        
                        
                        c[0] = 1;
                        if (n < a.MAXDEG) counter=n;
                        for (int i=1; i<counter; i++) {
                        
                                        c[i]  = n   * c[0] *  a[i]; // +  ((i-1)*n-1) * c[1] * a[i-1]; 
                                        const2 = 0;
                                        const1 = 1;
                                        for (int j=1; j<i; j++){
                                                const2 = const1;
                                                const1 *=  (i-j);
                                                const1 = const1 / j;
                                                //std::cout<< "const1="<< const1 << ", const2="<< const2 << std::endl;
                                                c[i] +=  (const1 * n - const2) * c[j] * a[i-j]; 
                                        }
                                
                                        if( i<= n) {                    
                                        
                                                //c[i] +=  (n-i+1) * c[i-1] * a[1];
                                                for(int k=0; k<i; k++) c[k] *= a[0]; 
                                
                                        } 
                        }*/
                        
                        
                        /*for (int i=counter+1; i<a.MAXDEG; i++) {
                                c[i]  = n  * a[i]; // +  ((i-1)*n-1) * c[1] * a[i-1]; 
                                
                                const2 = 0;
                                const1 = 1;
                                                                
                                for (int j=1; j<i; j++){
                                        const2 = const1;
                                        const1 *=  (i-j);
                                        const1 = const1/ j;
                                
                                        c[i] +=  (const1 * n - const2) * c[j] * a[i-j]; 
                                }
                        }*/
                        
                        
                        /*
                        break;
                        
                        
                }
        }
        

        return c;
} */


inline diffNumber log(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        c[0]=std::log(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i-1; j++){
                        c[i]=c[i] - j* c[j] * a[i-j] / i;
                }
                c[i]=(c[i] + a[i]) / a[0];
        }
        return c;
}

inline diffNumber sinh(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        diffNumber d(a.MAXDEG);
        c[0]=std::sinh(a[0]);
        d[0]=std::cosh(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i; j++){
                        c[i]=c[i] + j* a[j] * d[i-j] / i;
                        d[i]=d[i] + j* a[j] * c[i-j] / i;
                }
        }
        return c;
}

inline diffNumber cosh(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        diffNumber d(a.MAXDEG);
        c[0]=std::sinh(a[0]);
        d[0]=std::cosh(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                for (int j=1; j<=i; j++){
                        c[i]=c[i] + j* a[j] * d[i-j] / i;
                        d[i]=d[i] + j* a[j] * c[i-j] / i;
                }
        }
        return d;
}

inline diffNumber atan2(const diffNumber &a, const diffNumber&b)
{
        diffNumber c(a.MAXDEG);
        diffNumber d(a.MAXDEG);
        diffNumber e(a.MAXDEG);
        c=b/a;
        d=1.+c*c;
        e[0]=std::atan(c[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                e[i]=c[i]/d[0];
                for (int j=1; j<=i-1; j++){
                        e[i]=e[i] - j* e[j] * d[i-j] / i;
                }
        }
        return e;
}

inline diffNumber atan(const diffNumber &a)
{
        diffNumber d(a.MAXDEG);
        diffNumber e(a.MAXDEG);
        d=1.+a*a;
        e[0]=std::atan(a[0]);
        for (int i=1; i<a.MAXDEG+1; i++){
                e[i]=a[i]/d[0];
                for (int j=1; j<=i-1; j++){
                        e[i]=e[i] - j* e[j] * d[i-j] / i;
                }
        }
        return e;
}

// temporary version replacing the bad one below UNTIL MAXDEG=2
inline diffNumber square(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        c[0]=a[0]*a[0];
  if (a.MAXDEG>=1)
    c[1]=2.*a[0]*a[1];
  if (a.MAXDEG>=2)
    c[2]=2.*(a[0]*a[2]+a[1]*a[1]);
        return c;
}

/*
inline diffNumber square(const diffNumber &a)
{
        diffNumber c(a.MAXDEG);
        c[0]=a[0]*a[0];
        for (int i=1; i<a.MAXDEG+1; i++){
                c[i]=a[i];
                int m = (i + (i%2))/2;
                for (int j=0; j<=m-1; j++){
                        c[i]=c[i] + 2.* c[j] * c[i-j];
                }
                if (i%2==0) {
                        int n=i/2;
                        c[i]=c[i]+c[n]*c[n];
                }
        }
        return c;
}
*/


/*DIFFNUMBER*/
inline diffNumber::~diffNumber()
{
        
        if (data) delete[] data;
}

/*ostream &operator << (ostream &s, const diffNumber &a)
{
        diffNumber b(a.MAXDEG);
        s<<"< ";
        s<<a[0];
        //double fak=1.;
        for (int i=1; i<(a.MAXDEG+1);i++)
        {
                s<<", ";
                //s<<i;
                //fak=fak*i;
                //b[i]=a[i]*fak;
                s<<a[i];
        }
        s<<" >" ;
        return s;
}*/

inline std::ostream &operator << (std::ostream &s, const diffNumber &a)
{
        diffNumber b(a.MAXDEG);
        s<<"< ";
        s<<a[0];
        double fak=1.;
        for (int i=1; i<a.MAXDEG;i++)
        {
                s<<", ";
                fak=fak*i;
                b[i]=a[i]*fak;
                s << b[i];
        }
        s<<", ";
        fak=fak*a.MAXDEG;
        b[a.MAXDEG]=a[a.MAXDEG]*fak;
        s <<b[a.MAXDEG];
        s<<" >";
        return s;
}

#endif

---