ROOT_528-00b_version: math/mathcore/src/AdaptiveIntegratorMultiDim.cxx Source File

00001 // Implementation file for class
00002 // AdaptiveIntegratorMultiDim
00003 //
00004 #include "Math/IFunction.h"
00005 #include "Math/AdaptiveIntegratorMultiDim.h"
00006 #include "Math/Error.h"
00007 
00008 #include <cmath>
00009 
00010 namespace ROOT {
00011 namespace Math {
00012 
00013 
00014 
00015 AdaptiveIntegratorMultiDim::AdaptiveIntegratorMultiDim(double absTol, double relTol, unsigned int maxpts, unsigned int size):
00016    fDim(0), 
00017    fMinPts(0), 
00018    fMaxPts(maxpts),
00019    fSize(size), 
00020    fAbsTol(absTol),
00021    fRelTol(relTol),
00022    fResult(0), 
00023    fError(0), fRelError(0),
00024    fNEval(0),
00025    fStatus(-1),
00026    fFun(0)
00027 {
00028    // constructor - without passing a function
00029 }
00030 
00031 AdaptiveIntegratorMultiDim::AdaptiveIntegratorMultiDim( const IMultiGenFunction &f, double absTol, double relTol, unsigned int maxpts, unsigned int size):
00032    fDim(f.NDim()), 
00033    fMinPts(0), 
00034    fMaxPts(maxpts),
00035    fSize(size),
00036    fAbsTol(absTol),
00037    fRelTol(relTol),
00038    fResult(0), 
00039    fError(0), fRelError(0),
00040    fNEval(0),
00041    fStatus(-1),
00042    fFun(&f)
00043 {
00044    // constructur passing a multi-dimensional function interface
00045 }
00046 
00047 
00048 
00049 //double AdaptiveIntegratorMultiDim::Result() const { return fIntegrator->Result(); }
00050 //double AdaptiveIntegratorMultiDim::Error() const { return fIntegrator->Error(); }
00051 
00052 void AdaptiveIntegratorMultiDim::SetFunction(const IMultiGenFunction &f)
00053 {
00054    // set the integration function
00055    fFun = &f;
00056    fDim = f.NDim();
00057 }
00058 
00059 void AdaptiveIntegratorMultiDim::SetRelTolerance(double relTol){ this->fRelTol = relTol; }
00060 
00061 
00062 void AdaptiveIntegratorMultiDim::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }
00063 
00064 
00065 double AdaptiveIntegratorMultiDim::DoIntegral(const double* xmin, const double * xmax, bool absValue)
00066 {
00067    // References:
00068    //
00069    //   1.A.C. Genz and A.A. Malik, Remarks on algorithm 006:
00070    //     An adaptive algorithm for numerical integration over
00071    //     an N-dimensional rectangular region, J. Comput. Appl. Math. 6 (1980) 295-302.
00072    //   2.A. van Doren and L. de Ridder, An adaptive algorithm for numerical
00073    //     integration over an n-dimensional cube, J.Comput. Appl. Math. 2 (1976) 207-217.
00074   
00075    //to be changed later
00076    unsigned int n=fDim;
00077    bool kFALSE = false;
00078    bool kTRUE = true;
00079 
00080    double eps = fRelTol; //specified relative accuracy
00081    //output parameters
00082    fStatus = 0; //report status
00083    unsigned int nfnevl; //nr of function evaluations
00084    double relerr; //an estimation of the relative accuracy of the result
00085 
00086 
00087    double ctr[15], wth[15], wthl[15], z[15];
00088 
00089    static const double xl2 = 0.358568582800318073;//lambda_2
00090    static const double xl4 = 0.948683298050513796;//lambda_4
00091    static const double xl5 = 0.688247201611685289;//lambda_5
00092    static const double w2  = 980./6561; //weights/2^n
00093    static const double w4  = 200./19683;
00094    static const double wp2 = 245./486;//error weights/2^n
00095    static const double wp4 = 25./729;
00096 
00097    static const double wn1[14] = {     -0.193872885230909911, -0.555606360818980835,
00098                                        -0.876695625666819078, -1.15714067977442459,  -1.39694152314179743,
00099                                        -1.59609815576893754,  -1.75461057765584494,  -1.87247878880251983,
00100                                        -1.94970278920896201,  -1.98628257887517146,  -1.98221815780114818,
00101                                        -1.93750952598689219,  -1.85215668343240347,  -1.72615963013768225};
00102 
00103    static const double wn3[14] = {     0.0518213686937966768,  0.0314992633236803330,
00104                                        0.0111771579535639891,-0.00914494741655235473,-0.0294670527866686986,
00105                                        -0.0497891581567850424,-0.0701112635269013768, -0.0904333688970177241,
00106                                        -0.110755474267134071, -0.131077579637250419,  -0.151399685007366752,
00107                                        -0.171721790377483099, -0.192043895747599447,  -0.212366001117715794};
00108 
00109    static const double wn5[14] = {         0.871183254585174982e-01,  0.435591627292587508e-01,
00110                                            0.217795813646293754e-01,  0.108897906823146873e-01,  0.544489534115734364e-02,
00111                                            0.272244767057867193e-02,  0.136122383528933596e-02,  0.680611917644667955e-03,
00112                                            0.340305958822333977e-03,  0.170152979411166995e-03,  0.850764897055834977e-04,
00113                                            0.425382448527917472e-04,  0.212691224263958736e-04,  0.106345612131979372e-04};
00114 
00115    static const double wpn1[14] = {   -1.33196159122085045, -2.29218106995884763,
00116                                       -3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
00117                                       -4.76131687242798352, -5.03566529492455417, -5.17283950617283939,
00118                                       -5.17283950617283939, -5.03566529492455417, -4.76131687242798352,
00119                                       -4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
00120 
00121    static const double wpn3[14] = {     0.0445816186556927292, -0.0240054869684499309,
00122                                         -0.0925925925925925875, -0.161179698216735251,  -0.229766803840877915,
00123                                         -0.298353909465020564,  -0.366941015089163228,  -0.435528120713305891,
00124                                         -0.504115226337448555,  -0.572702331961591218,  -0.641289437585733882,
00125                                         -0.709876543209876532,  -0.778463648834019195,  -0.847050754458161859};
00126 
00127    double result = 0;
00128    double abserr = 0;
00129    fStatus  = 3;
00130    nfnevl = 0;
00131    relerr = 0;
00132    // does not work for 1D functions
00133    if (n < 2 || n > 15) { 
00134       MATH_WARN_MSGVAL("AdaptiveIntegratorMultiDim::Integral","Wrong function dimension",n); 
00135       return 0;
00136    }
00137 
00138    double twondm = std::pow(2.0,static_cast<int>(n));
00139    //unsigned int minpts = Int_t(twondm)+ 2*n*(n+1)+1;
00140 
00141    unsigned int ifncls = 0;
00142    bool  ldv   = kFALSE;
00143    unsigned int irgnst = 2*n+3;
00144    unsigned int  irlcls = (unsigned int)(twondm) +2*n*(n+1)+1;//minimal number of nodes in n dim
00145    unsigned int isbrgn = irgnst;
00146    unsigned int isbrgs = irgnst;
00147 
00148 
00149    unsigned int minpts = fMinPts; 
00150    unsigned int maxpts = std::max(fMaxPts, irlcls) ;//specified maximal number of function evaluations
00151 
00152    if (minpts < 1)      minpts = irlcls;
00153    if (maxpts < minpts) maxpts = 10*minpts;
00154 
00155    // The original agorithm expected a working space array WK of length IWK
00156    // with IWK Length ( >= (2N + 3) * (1 + MAXPTS/(2**N + 2N(N + 1) + 1))/2).
00157    // Here, this array is allocated dynamically
00158 
00159    unsigned int iwk = std::max( fSize, irgnst*(1 +maxpts/irlcls)/2 );
00160    double *wk = new double[iwk+10];
00161 
00162    unsigned int j; 
00163    for (j=0; j<n; j++) {
00164       ctr[j] = (xmax[j] + xmin[j])*0.5;//center of a hypercube
00165       wth[j] = (xmax[j] - xmin[j])*0.5;//its width
00166    }
00167 
00168    double rgnvol, sum1, sum2, sum3, sum4, sum5, difmax, f2, f3, dif, aresult;
00169    double rgncmp=0, rgnval, rgnerr;
00170 
00171    unsigned int j1, k, l, m, idvaxn=0, idvax0=0, isbtmp, isbtpp;
00172  
00173    //InitArgs(z,fParams);
00174 
00175 L20:
00176    rgnvol = twondm;//=2^n
00177    for (j=0; j<n; j++) {
00178       rgnvol *= wth[j]; //region volume
00179       z[j]    = ctr[j]; //temporary node
00180    }
00181    sum1 = (*fFun)((const double*)z);//EvalPar(z,fParams); //evaluate function
00182 
00183    difmax = 0;
00184    sum2   = 0;
00185    sum3   = 0;
00186 
00187    //loop over coordinates
00188    for (j=0; j<n; j++) {
00189       z[j]    = ctr[j] - xl2*wth[j];
00190       if (absValue) f2 = std::abs((*fFun)(z));
00191       else          f2 = (*fFun)(z);
00192       z[j]    = ctr[j] + xl2*wth[j];
00193       if (absValue) f2 += std::abs((*fFun)(z));
00194       else          f2 += (*fFun)(z);
00195       wthl[j] = xl4*wth[j];
00196       z[j]    = ctr[j] - wthl[j]; 
00197       if (absValue) f3 = std::abs((*fFun)(z));
00198       else          f3 = (*fFun)(z);
00199       z[j]    = ctr[j] + wthl[j];
00200       if (absValue) f3 += std::abs((*fFun)(z));
00201       else          f3 += (*fFun)(z);
00202       sum2   += f2;//sum func eval with different weights separately
00203       sum3   += f3;//for a given region
00204       dif     = std::abs(7*f2-f3-12*sum1);
00205       //storing dimension with biggest error/difference (?)
00206       if (dif >= difmax) {
00207          difmax=dif;
00208          idvaxn=j+1;
00209       }
00210       z[j]    = ctr[j];
00211    }
00212 
00213    sum4 = 0;
00214    for (j=1;j<n;j++) {
00215       j1 = j-1;
00216       for (k=j;k<n;k++) {
00217          for (l=0;l<2;l++) {
00218             wthl[j1] = -wthl[j1];
00219             z[j1]    = ctr[j1] + wthl[j1];
00220             for (m=0;m<2;m++) {
00221                wthl[k] = -wthl[k];
00222                z[k]    = ctr[k] + wthl[k];
00223                if (absValue) sum4 += std::abs((*fFun)(z));
00224                else            sum4 += (*fFun)(z);
00225             }
00226          }
00227          z[k] = ctr[k];
00228       }
00229       z[j1] = ctr[j1];
00230    }
00231 
00232    sum5 = 0;
00233 
00234    for (j=0;j<n;j++) {
00235       wthl[j] = -xl5*wth[j];
00236       z[j] = ctr[j] + wthl[j];
00237    }
00238 L90: //sum over end nodes ~gray codes
00239    if (absValue) sum5 += std::abs((*fFun)(z));
00240    else          sum5 += (*fFun)(z);
00241    for (j=0;j<n;j++) {
00242       wthl[j] = -wthl[j];
00243       z[j] = ctr[j] + wthl[j];
00244       if (wthl[j] > 0) goto L90;
00245    }
00246 
00247    rgncmp  = rgnvol*(wpn1[n-2]*sum1+wp2*sum2+wpn3[n-2]*sum3+wp4*sum4);
00248    rgnval  = wn1[n-2]*sum1+w2*sum2+wn3[n-2]*sum3+w4*sum4+wn5[n-2]*sum5;
00249    rgnval *= rgnvol;
00250    // avoid difference of too small numbers
00251    //rgnval = 1.0E-30;
00252    //rgnerr  = TMath::Max( std::abs(rgnval-rgncmp), TMath::Max(std::abs(rgncmp), std::abs(rgnval) )*4.0E-16 );
00253    rgnerr  = std::abs(rgnval-rgncmp);//compares estim error with expected error
00254 
00255    result += rgnval;
00256    abserr += rgnerr;
00257    ifncls += irlcls;
00258    aresult = std::abs(result);
00259    //if (result > 0 && aresult< 1e-100) {
00260    //   delete [] wk;
00261    //   fStatus = 0;  //function is probably symmetric ==> integral is null: not an error
00262    //   return result;
00263    //}
00264 
00265    //if division
00266    if (ldv) {
00267    L110:
00268       isbtmp = 2*isbrgn;
00269       if (isbtmp > isbrgs) goto L160;
00270       if (isbtmp < isbrgs) {
00271          isbtpp = isbtmp + irgnst;
00272          if (wk[isbtmp-1] < wk[isbtpp-1]) isbtmp = isbtpp;
00273       }
00274       if (rgnerr >= wk[isbtmp-1]) goto L160;
00275       for (k=0;k<irgnst;k++) {
00276          wk[isbrgn-k-1] = wk[isbtmp-k-1];
00277       }
00278       isbrgn = isbtmp;
00279       goto L110;
00280    }
00281 L140:
00282    isbtmp = (isbrgn/(2*irgnst))*irgnst;
00283    if (isbtmp >= irgnst && rgnerr > wk[isbtmp-1]) {
00284       for (k=0;k<irgnst;k++) {
00285          wk[isbrgn-k-1] = wk[isbtmp-k-1];
00286       }
00287       isbrgn = isbtmp;
00288       goto L140;
00289    }
00290 
00291 L160: //to divide or not
00292    wk[isbrgn-1] = rgnerr;//storing value & error in last 
00293    wk[isbrgn-2] = rgnval;//table records
00294    wk[isbrgn-3] = double(idvaxn);//coordinate with biggest error
00295    for (j=0;j<n;j++) {
00296       isbtmp = isbrgn-2*j-4;
00297       wk[isbtmp]   = ctr[j];
00298       wk[isbtmp-1] = wth[j];
00299    }
00300    if (ldv) {//divison along chosen coordinate
00301       ldv = kFALSE;
00302       ctr[idvax0-1] += 2*wth[idvax0-1];
00303       isbrgs += irgnst;//updating the number of nodes/regions(?)
00304       isbrgn  = isbrgs;
00305       goto L20;
00306    }
00307    //if no divisions to be made..
00308    relerr = abserr;
00309    if (aresult != 0)  relerr = abserr/aresult;
00310 
00311 
00312    if (relerr < 1e-1 && aresult < 1e-20) fStatus = 0;
00313    if (relerr < 1e-3 && aresult < 1e-10) fStatus = 0;
00314    if (relerr < 1e-5 && aresult < 1e-5)  fStatus = 0;
00315    if (isbrgs+irgnst > iwk) fStatus = 2;
00316    if (ifncls+2*irlcls > maxpts) {
00317       if (sum1==0 && sum2==0 && sum3==0 && sum4==0 && sum5==0){
00318          fStatus = 0;
00319          result = 0;
00320       }
00321       else
00322          fStatus = 1;
00323    }
00324    //..and accuracy appropriare
00325    if (relerr < eps && ifncls >= minpts) fStatus = 0;  // We do not use the absolute error.
00326    //if (relerr < eps* aresult && abserr < eps && ifncls >= minpts) fStatus = 0;
00327    if (fStatus == 3) {
00328       ldv = kTRUE;
00329       isbrgn  = irgnst;
00330       abserr -= wk[isbrgn-1];
00331       result -= wk[isbrgn-2];
00332       idvax0  = (unsigned int)(wk[isbrgn-3]);
00333       for (j=0;j<n;j++) {
00334          isbtmp = isbrgn-2*j-4;
00335          ctr[j] = wk[isbtmp];
00336          wth[j] = wk[isbtmp-1];
00337       }
00338       wth[idvax0-1]  = 0.5*wth[idvax0-1];
00339       ctr[idvax0-1] -= wth[idvax0-1];
00340       goto L20;
00341    }
00342    nfnevl = ifncls;       //number of function evaluations performed.
00343    fResult = result;
00344    fError = abserr;//wk[isbrgn-1];
00345    fRelError = relerr;
00346    fNEval = nfnevl;
00347    delete [] wk;
00348   
00349    return result;         //an approximate value of the integral
00350 }
00351 
00352 
00353   
00354 double AdaptiveIntegratorMultiDim::Integral(const IMultiGenFunction &f, const double* xmin, const double * xmax)
00355 {
00356    // calculate integral passing a function object
00357    fFun = &f;
00358    return Integral(xmin, xmax);
00359 
00360 }
00361 
00362 ROOT::Math::IntegratorMultiDimOptions  AdaptiveIntegratorMultiDim::Options() const { 
00363    // return the used options
00364    ROOT::Math::IntegratorMultiDimOptions opt; 
00365    opt.SetAbsTolerance(fAbsTol); 
00366    opt.SetRelTolerance(fRelTol); 
00367    opt.SetNCalls(fMaxPts); 
00368    opt.SetWKSize(fSize); 
00369    opt.SetIntegrator("ADAPTIVE");
00370    return opt; 
00371 }
00372 
00373 void AdaptiveIntegratorMultiDim::SetOptions(const ROOT::Math::IntegratorMultiDimOptions & opt)
00374 {
00375    //   set integration options
00376    if (opt.IntegratorType() != IntegrationMultiDim::kADAPTIVE) {
00377       MATH_ERROR_MSG("AdaptiveIntegratorMultiDim::SetOptions","Invalid options");
00378       return;
00379    }      
00380    SetAbsTolerance( opt.AbsTolerance() );
00381    SetRelTolerance( opt.RelTolerance() );
00382    SetMaxPts( opt.NCalls() );
00383    SetSize( opt.WKSize() );
00384 }
00385 
00386 } // namespace Math
00387 } // namespace ROOT
00388 
00389 
00390