COCONUT API: diffI.h Source File

00001 #ifndef __DIFFI_H__
00002 #define __DIFFI_H__
00003 
00004 #include "interval.h"
00005 #include <stdlib.h>
00006 #include <iostream>
00007 #include <coconut_vect.h>
00008 
00009 using namespace coco;
00010 
00011 const int STANDARD_MAXDEG_DIFFINTERVAL=4;
00012 
00013 /****************/
00014 /*DIFFINTERVALL*/
00015 /****************/
00016 
00017 //CLASSDEF
00018 
00019 class diffI
00020 {
00021         private:
00022                 typedef interval I;
00023 
00024         private:
00025                 I *data;
00026         public:
00027                 int MAXDEG;
00028                 inline diffI(void);
00029                 inline diffI(double a);
00030                 inline diffI(double a, double b);
00031                 inline diffI(int n);
00032                 inline diffI(double a, int n);
00033                 inline diffI(double a, double b, int n);
00034                 inline diffI(const I& a, int n);
00035                 inline diffI(const diffI &a);
00036                 virtual ~diffI();
00037         
00038                 inline diffI& operator=(const diffI &a);
00039                 inline I& operator[] (int pos) const;
00040         
00041                 inline static diffI getVar(const double &x);
00042                 inline static diffI getVar(const double &x, const double &y);
00043                 inline static diffI getVar(const double &x, int n);
00044                 inline static diffI getVar(const double &x, const double &y, int n);
00045                 inline static diffI getVar(const I &x, int n);
00046                 inline void setDeg(int n);
00047 };
00048 
00049 //OPERATORS
00050 inline diffI operator +(const diffI &a, const diffI &b);
00051 inline diffI operator +(const double &a, const diffI &b);
00052 inline diffI operator +(const diffI &a, const double &b);
00053 inline diffI operator -(const diffI &a, const diffI &b);
00054 inline diffI operator -(const double &a, const diffI &b);
00055 inline diffI operator -(const diffI &a, const double &b);
00056 inline diffI operator -(const diffI &a);
00057 inline diffI operator *(const diffI &a, const diffI &b);
00058 inline diffI operator *(const double &a, const diffI &b);
00059 inline diffI operator *(const diffI &a, const double &b);
00060 inline diffI operator /(const diffI &a, const diffI &b);
00061 inline diffI operator /(const diffI &a, const double &b);
00062 
00063 //OUTPUT
00064 std::ostream &operator << (std::ostream &s, const diffI &a);
00065 
00066 //FUNCTIONS
00067 inline diffI exp(const diffI &a);
00068 inline diffI sin(const diffI &a);
00069 inline diffI cos(const diffI &a);
00070 inline diffI sqrt(const diffI &a);
00071 
00072 inline diffI power(const diffI &a, int n);
00073 inline diffI power(const diffI &a, const diffI& b);// added by Simon Haller-> needed for ForwardHessEvaluator
00074 
00075 inline diffI log(const diffI &a);
00076 inline diffI sinh(const diffI &a);
00077 inline diffI cosh(const diffI &a);
00078 inline diffI atan2(const diffI &a, const diffI &b);
00079 inline diffI atan(const diffI &a);
00080 inline diffI square(const diffI &a);
00081 
00082 /****************/
00083 /*INLINECODE*/
00084 /****************/
00085 
00086 /*diffI*/
00087 
00088 //CONSTRUCTORS
00089 inline diffI::diffI(void)
00090 {
00091         MAXDEG=STANDARD_MAXDEG_DIFFINTERVAL;
00092         data=new I[MAXDEG+1];
00093         memset(data, 0, sizeof(I)*(MAXDEG+1));
00094 }
00095 
00096 inline diffI::diffI(int n)
00097 {
00098         MAXDEG=n;
00099         data=new I[MAXDEG+1];
00100         memset(data, 0, sizeof(I)*(MAXDEG+1));
00101 }
00102 
00103 inline diffI::diffI(double a)
00104 {
00105         MAXDEG=STANDARD_MAXDEG_DIFFINTERVAL;
00106         data= new I[MAXDEG+1];
00107         data[0]=a;
00108         memset(data+1, 0, sizeof(I)*MAXDEG);
00109 }
00110 
00111 inline diffI::diffI(double a, double b)
00112 {
00113         MAXDEG=STANDARD_MAXDEG_DIFFINTERVAL;
00114         data= new I[MAXDEG+1];
00115         data[0]=I(a,b);
00116         memset(data+1, 0, sizeof(I)*MAXDEG);
00117 }
00118 
00119 inline diffI::diffI(double a, int n)
00120 {
00121         MAXDEG=n;
00122         data= new I[MAXDEG+1];
00123         data[0]=a;
00124         memset(data+1, 0, sizeof(I)*MAXDEG);
00125 }
00126 
00127 inline diffI::diffI(double a, double b, int n)
00128 {
00129         MAXDEG=n;
00130         data= new I[MAXDEG+1];
00131         data[0]=I(a,b);
00132         memset(data+1, 0, sizeof(I)*MAXDEG);
00133 }
00134 
00135 inline diffI::diffI(const I& a, int n)
00136 {
00137         MAXDEG=n;
00138         data= new I[MAXDEG+1];
00139         data[0]=a;
00140         memset(data+1, 0, sizeof(I)*MAXDEG);
00141 }
00142 
00143 inline diffI::diffI(const diffI &d)
00144 {
00145         MAXDEG = d.MAXDEG;
00146         data = new I[MAXDEG+1];
00147         memcpy(data, d.data, (MAXDEG+1)*sizeof(I));
00148 }
00149 
00150 inline diffI& diffI::operator=(const diffI &a)
00151 {
00152         if ((MAXDEG!=-1) && (MAXDEG != a.MAXDEG)){
00153                 std::cerr<<"= degree not the same: "<< MAXDEG <<" != "<< a.MAXDEG;
00154                 abort();
00155         }
00156         
00157         if (this != &a){
00158                 if (data) delete [] data;
00159                 data = new I[a.MAXDEG+1];
00160                 for (int i=0; i<a.MAXDEG+1; i++)
00161                         data[i]=a[i];
00162         }
00163         return (*this);
00164 }
00165 
00166 
00167 inline interval &diffI::operator[] (int pos) const
00168 {
00169 //      if (pos<
00170         return data[pos];
00171 }
00172 
00173 
00174 
00175 //VARIABLES
00176 inline diffI diffI::getVar(const double &x)
00177 {
00178         diffI c;
00179         c[0]=x;
00180         c[1]=1;
00181         return c;
00182 }
00183 
00184 inline diffI diffI::getVar(const double &x, const double &y)
00185 {
00186         diffI c;
00187         c[0]=I(x,y);
00188         c[1]=1;
00189         return c;
00190 }
00191 
00192 inline diffI diffI::getVar(const double &x, int n)
00193 {
00194         diffI c(n);
00195         c[0]=x;
00196         c[1]=1;
00197         return c;
00198 }
00199 
00200 inline diffI diffI::getVar(const double &x, const double &y, int n)
00201 {
00202         diffI c(n);
00203         c[0]=I(x,y);
00204         c[1]=1;
00205         return c;
00206 }
00207 
00208 inline diffI diffI::getVar(const I &x, int n)
00209 {
00210         diffI c(n);
00211         c[0]=x;
00212         c[1]=1;
00213         return c;
00214 }
00215 
00216 inline void diffI::setDeg(int n)
00217 {
00218         __vectgen(I,tmp,n+1);
00219         if (MAXDEG>n){
00220                 for (int i=0; i<n+1; i++)
00221                         tmp[i]=data[i];
00222                 delete [] data;
00223                 data=new I[n+1];
00224                 for (int i=0; i<n+1; i++)
00225                         data[i]=tmp[i];
00226         }
00227         else {
00228                 for (int i=0; i<MAXDEG+1; i++)
00229                         tmp[i]=data[i];
00230                 delete [] data;
00231                 data=new I[n+1];
00232                 for (int i=0; i<MAXDEG+1; i++)
00233                         data[i]=tmp[i];
00234                 for (int j=MAXDEG+1; j<n+1; j++)
00235                         data[j]=0.;
00236         }
00237         MAXDEG=n;
00238 }
00239 
00240 
00241 //OPERATORS
00242 inline diffI operator +(const diffI &a, const diffI &b)
00243 {
00244         if (a.MAXDEG != b.MAXDEG){
00245                 std::cerr<<"+ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
00246                 abort();
00247         }
00248         diffI c(a.MAXDEG);
00249         for (int i=0; i<a.MAXDEG+1; i++){
00250                 c[i]=a[i] + b[i];
00251         }
00252         return c;
00253 }
00254 
00255 inline diffI operator +(const double &a, const diffI &b)
00256 {
00257         diffI c(b.MAXDEG);
00258         c=b;
00259         c[0]=b[0]+a;
00260         return c;
00261 }
00262 
00263 inline diffI operator +(const diffI &a, const double &b)
00264 {
00265         diffI c(a.MAXDEG);
00266         c=a;
00267         c[0]=a[0]+b;
00268         return c;
00269 }
00270 
00271 inline diffI operator -(const diffI &a, const diffI &b)
00272 {
00273         if (a.MAXDEG != b.MAXDEG){
00274                 std::cerr<<"- degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
00275                 abort();
00276         }
00277         diffI c(a.MAXDEG);
00278         for (int i=0; i<a.MAXDEG+1; i++){
00279                 c[i]=a[i] - b[i];
00280         }
00281         return c;
00282 }
00283 
00284 inline diffI operator -(const double &a, const diffI &b)
00285 {
00286         diffI c(b.MAXDEG);
00287         c=b;
00288         c[0]=a-b[0];
00289         return c;
00290 }
00291 
00292 inline diffI operator -(const diffI &a, const double &b)
00293 {
00294         diffI c(a.MAXDEG);
00295         c=a;
00296         c[0]=a[0]-b;
00297         return c;
00298 }
00299 
00300 inline diffI operator -(const diffI &a)
00301 {
00302         diffI c(a.MAXDEG);
00303         for (int i=0; i<a.MAXDEG+1; i++){
00304                 c[i]=-a[i];
00305         }
00306         return c;
00307 }
00308 
00309 inline diffI operator *(const diffI &a, const diffI &b)
00310 {
00311         if (a.MAXDEG != b.MAXDEG){
00312                 std::cerr<<"* degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
00313                 abort();
00314         }
00315         diffI c(a.MAXDEG);
00316         for (int i=0; i<a.MAXDEG+1; i++){
00317                 for (int j=0; j<=i; j++){
00318                         c[i]=c[i] + a[j] * b[i-j];
00319                 }
00320         }
00321         return c;
00322 }
00323 
00324 inline diffI operator *(const double &a, const diffI &b)
00325 {
00326         diffI c(b.MAXDEG);
00327         for (int i=0; i<b.MAXDEG+1; i++){
00328                 c[i]=a * b[i];
00329         }
00330         return c;
00331 }
00332 
00333 inline diffI operator *(const diffI &a, const double &b)
00334 {
00335         diffI c(a.MAXDEG);
00336         for (int i=0; i<a.MAXDEG+1; i++){
00337                 c[i]=b * a[i];
00338         }
00339         return c;
00340 }
00341 
00342 inline diffI operator /(const diffI &a, const diffI &b)
00343 {
00344         if (a.MAXDEG != b.MAXDEG){
00345                 std::cerr<<"/ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
00346                 abort();
00347         }
00348         diffI c(a.MAXDEG);
00349         c[0]=a[0] / b[0];
00350         for (int i=1; i<a.MAXDEG+1; i++){
00351                 c[i] = 0.;
00352                 for (int j=0; j<i; j++){
00353                         c[i] += c[j] * b[i-j];
00354                 }
00355                 c[i]=1. / b[0] * (a[i] - c[i]);
00356         }
00357         return c;
00358 }
00359 
00360 inline diffI operator /(const double &a, const diffI &b)
00361 {
00362         diffI c(b.MAXDEG);
00363         c[0]=a / b[0];
00364         for (int i=1; i<b.MAXDEG+1; i++){
00365                 c[i]=0.;
00366                 for (int j=0; j<i; j++){
00367                         c[i] += c[j] * b[i-j];
00368                 }
00369                 c[i]= - 1. / b[0] * c[i];
00370         }
00371         return c;
00372 }
00373 
00374 
00375 inline diffI operator /(const diffI &a, const double &b)
00376 {
00377         diffI c(a.MAXDEG);
00378         for (int i=0; i<a.MAXDEG+1; i++){
00379                 c[i]=a[i] / b;
00380         }
00381         return c;
00382 }
00383 
00384 inline diffI exp(const diffI &a)
00385 {
00386         diffI c(a.MAXDEG);
00387         c[0]=exp(a[0]);
00388         for (int i=1; i<a.MAXDEG+1; i++){
00389                 for (int j=1; j<=i; j++){
00390                         c[i]=c[i] + j* a[j] * c[i-j] / i;
00391                 }
00392         }
00393         return c;
00394 }
00395 
00396 inline diffI sin(const diffI &a)
00397 {
00398         diffI c(a.MAXDEG);
00399         diffI d(a.MAXDEG);
00400         c[0]=sin(a[0]);
00401         //cout<<a[0]<<","<<c[0]<<endl;
00402         c[0]=sin(a[0]); //wenn diese zeile geloescht wird ist c[0]=[empty]
00403         //cout<<a[0]<<","<<c[0]<<endl;
00404         //cout<<c[0];
00405         d[0]=cos(a[0]);
00406         d[0]=cos(a[0]);
00407         for (int i=1; i<a.MAXDEG+1; i++){
00408                 for (int j=1; j<=i; j++){
00409                         c[i]=c[i] + j* a[j] * d[i-j] / i;
00410                         d[i]=d[i] - j* a[j] * c[i-j] / i;
00411                 }
00412         }
00413         return c;
00414 }
00415 
00416 inline diffI cos(const diffI &a)
00417 {
00418         diffI c(a.MAXDEG);
00419         diffI d(a.MAXDEG);
00420         //I tmp;
00421         //tmp=a[0];
00422         //cout<<sin(tmp);
00423         c[0]=sin(a[0]);
00424         c[0]=sin(a[0]); //wenn diese zeile geloescht wird ist c[0]=[empty]
00425         //cout<<c[0];
00426         d[0]=cos(a[0]);
00427         d[0]=cos(a[0]);
00428         for (int i=1; i<a.MAXDEG+1; i++){
00429                 for (int j=1; j<=i; j++){
00430                         c[i]=c[i] + j* a[j] * d[i-j] / i;
00431                         d[i]=d[i] - j* a[j] * c[i-j] / i;
00432                 }
00433         }
00434         return d;
00435 }
00436 
00437 inline diffI sqrt(const diffI &a)
00438 {
00439         diffI c(a.MAXDEG);
00440         c[0]=sqrt(a[0]);
00441         for (int i=1; i<a.MAXDEG+1; i++){
00442                 c[i]=a[i];
00443                 int m = (i + (i%2))/2;
00444                 for (int j=1; j<=m-1; j++){
00445                   if(j == i-j)
00446                         c[i]=c[i] - 2.* sqr(c[j]);
00447                   else
00448                         c[i]=c[i] - 2.* c[j] * c[i-j];
00449                 }
00450                 if (i%2==0) {
00451                         int n=i/2;
00452                         c[i]=c[i]-sqr(c[n]);
00453                 }
00454                 c[i]=c[i] * 0.5 / c[0];
00455         }
00456         return c;
00457 }
00458 
00459 
00460 inline diffI power(const diffI &a, int n)
00461 {
00462         diffI c(a.MAXDEG);
00463         diffI d(a.MAXDEG);
00464         c[0]=power(a[0],n);
00465 
00466 // temporary improvement fix by Mihaly, 02.06.2009.
00467   if (a.MAXDEG>=1)
00468     c[1]=n*power(a[0],n-1)*a[1];
00469   if (a.MAXDEG>=2)
00470     c[2]=n*power(a[0],n-2)*((n-1)*sqr(a[1])+a[0]*a[2]);
00471   for (int i=3; i<a.MAXDEG+1; ++i){
00472 
00473 //      for (int i=1; i<a.MAXDEG+1; i++){
00474                 for (int j=0; j<=i-1; j++){
00475                         c[i]=c[i] + (n * (i - j) - j)* c[j] * a[i-j] / i;
00476                 }
00477                 c[i]=c[i] / a[0];
00478         }
00479         d=a;
00480         for (int i=1; i<n; i++){
00481                 d=d*a;
00482         }
00483         for (int i=0; i<a.MAXDEG+1; i++){
00484                 //std::cout<<"intersect: c"<<c[i]<<" d ";
00485                 c[i].intersectwith(d[i]);
00486                 //std::cout<<d[i]<<"\n";
00487         }
00488         return c;
00489 }
00490 
00491 inline diffI power(const diffI &a, const diffI &b)
00492 {       
00493         if (a.MAXDEG != b.MAXDEG){
00494                 std::cerr<<"+ degree not the same: "<< a.MAXDEG <<" != "<< b.MAXDEG;
00495                 abort();
00496         }
00497         diffI c(a.MAXDEG);
00498         
00499         if (a[0] == 0.0)
00500                 for (int i=0; i<a.MAXDEG+1; i++) c[i] = 0;
00501         else
00502                 c= exp(b * log(a)); // lazyman's approach of implementing the power-function...
00503         
00504         return c;
00505 } 
00506 
00507 inline diffI log(const diffI &a)
00508 {
00509         diffI c(a.MAXDEG);
00510         c[0]=log(a[0]);
00511         for (int i=1; i<a.MAXDEG+1; i++){
00512                 for (int j=1; j<=i-1; j++){
00513                         c[i]=c[i] - j* c[j] * a[i-j] / i;
00514                 }
00515                 c[i]=(c[i] + a[i]) / a[0];
00516         }
00517         return c;
00518 }
00519 
00520 inline diffI sinh(const diffI &a)
00521 {
00522         diffI c(a.MAXDEG);
00523         diffI d(a.MAXDEG);
00524         c[0]=sinh(a[0]);
00525         d[0]=cosh(a[0]);
00526         for (int i=1; i<a.MAXDEG+1; i++){
00527                 for (int j=1; j<=i; j++){
00528                         c[i]=c[i] + j* a[j] * d[i-j] / i;
00529                         d[i]=d[i] + j* a[j] * c[i-j] / i;
00530                 }
00531         }
00532         return c;
00533 }
00534 
00535 inline diffI cosh(const diffI &a)
00536 {
00537         diffI c(a.MAXDEG);
00538         diffI d(a.MAXDEG);
00539         c[0]=sinh(a[0]);
00540         d[0]=cosh(a[0]);
00541         for (int i=1; i<a.MAXDEG+1; i++){
00542                 for (int j=1; j<=i; j++){
00543                         c[i]=c[i] + j* a[j] * d[i-j] / i;
00544                         d[i]=d[i] + j* a[j] * c[i-j] / i;
00545                 }
00546         }
00547         return d;
00548 }
00549 
00550 inline diffI atan2(const diffI &a, const diffI&b)
00551 {
00552         diffI c(a.MAXDEG);
00553         diffI d(a.MAXDEG);
00554         diffI e(a.MAXDEG);
00555         c=b/a;
00556         d=1.+square(c);
00557         e[0]=atan(c[0]);
00558         for (int i=1; i<a.MAXDEG+1; i++){
00559                 e[i]=c[i]/d[0];
00560                 for (int j=1; j<=i-1; j++){
00561                         e[i]=e[i] - j* e[j] * d[i-j] / i;
00562                 }
00563         }
00564         return e;
00565 }
00566 
00567 inline diffI atan(const diffI &a)
00568 {
00569         diffI d(a.MAXDEG);
00570         diffI e(a.MAXDEG);
00571         d=1.+square(a);
00572         e[0]=atan(a[0]);
00573         for (int i=1; i<a.MAXDEG+1; i++){
00574                 e[i]=a[i]/d[0];
00575                 for (int j=1; j<=i-1; j++){
00576                         e[i]=e[i] - j* e[j] * d[i-j] / i;
00577                 }
00578         }
00579         return e;
00580 }
00581 
00582 inline diffI square(const diffI &a)
00583 {
00584   diffI c(a.MAXDEG);
00585   c[0]=sqr(a[0]);
00586   if (a.MAXDEG>=1)
00587     c[1]=2.*a[0]*a[1];
00588   if (a.MAXDEG>=2)
00589     c[2]=2.*(a[0]*a[2]+sqr(a[1]));
00590   if (a.MAXDEG>=3)  // For j+k=i, j<=k, the coefficient of a[j]*a[k] in c[i] is
00591   {                 // 2*(i over j) if j != i/2, and (i over j) if j == i/2.
00592                     // These coefficients will be computed recursively by the
00593                     // the Pascal-triangle rules, and will be stored in psc.
00594     int m;
00595     std::vector<int> psc(a.MAXDEG/2 + 1);
00596     psc[0] = 2;
00597     psc[1] = 2;
00598     for (int i=3; i<=a.MAXDEG; ++i)
00599     {
00600       m = i/2;
00601       c[i] = 2*a[0]*a[i];
00602 
00603       if (i%2 == 1)
00604       {
00605         psc[m] = psc[m-1] + 2*psc[m];
00606         c[i] = c[i] + psc[m]*a[m]*a[m+1];
00607       }
00608       else
00609       {
00610         psc[m] = psc[m-1];
00611         c[i] = c[i] + psc[m]*sqr(a[m]);
00612       }
00613 
00614       for (int j=m-1; j>=1; --j)
00615       {
00616         psc[j] = psc[j-1] + psc[j];
00617         if(j == i-j)
00618           c[i] = c[i] + psc[j]*sqr(a[j]);
00619         else
00620           c[i] = c[i] + psc[j]*a[j]*a[i-j];
00621       }
00622     }
00623   }
00624   return c;
00625 }
00626 
00627 /*
00628 // this is incorrect and is replaced by Mihaly with the one above
00629 inline diffI square(const diffI &a)
00630 {
00631         diffI c(a.MAXDEG);
00632         c[0]=sqr(a[0]);
00633         for (int i=1; i<a.MAXDEG+1; i++){
00634                 c[i]=a[i];
00635                 int m = (i + (i%2))/2;
00636                 for (int j=0; j<=m-1; j++){
00637                         c[i]=c[i] + 2.* c[j] * c[i-j];
00638                 }
00639                 if (i%2==0) {
00640                         int n=i/2;
00641                         c[i]=c[i]+sqr(c[n]);
00642                 }
00643         }
00644         return c;
00645 }
00646 */
00647 
00648 
00649 /*diffI*/
00650 inline diffI::~diffI()
00651 {
00652 //      cout<<"diffI delete";
00653         if (data) delete[] data;
00654 }
00655 
00656 inline std::ostream &operator << (std::ostream &s, const diffI &a)
00657 {
00658         diffI b(a.MAXDEG);
00659         s<<"< f0=";
00660         s<<a[0];
00661         double fak=1.;
00662         for (int i=1; i<(a.MAXDEG+1);i++)
00663         {
00664                 s<<", f";
00665                 s<<i;
00666                 fak=fak*i;
00667                 b[i]=a[i]*fak;
00668                 s<<"=" <<b[i];
00669         }
00670         s<<" >";
00671         return s;
00672 }
00673 
00674 #endif