example_NIntp.txt

The following source file shows the application of multi-dimensional interpolation calc::NIntp objects and the calc::OptIntp object which analyzes and reduces interpolation errors. The objective in this example is to show the similar uses of the three interpolation objects in optimizing and curve-fitting data.

First, a two-dimensional curve is defined by a generating function curve(). Dictated by the optimizing calc::OptIntp object, this function must accept a vector array of <xType*> which contains a evaluation data for each dimension. The generating function should return a single value of <yType>. The generating function for this example quantifies the Gauss bell-curve in two dimensions.

Parameters should be set which clarify the interpolation. The interpolation interval of the curve is from the origin (the top of the bell) out to unit length in the positive quadrant. Tolerances are set indicating the magnitude of values that will be encountered during the interpolation optimization. If the difference between two interpolation points is less than the given x-tolerance, then the points are assumed identical. The data count by dimension array, ND[NDIM], defines the structure of the interpolation table. The length of the array equals the number of dimensions NDIM in the interpolation table. Each value in the array indictes how many x-intervals the interpolation table has for the specified dimension. From the example below, the x-array of pointers to an array of data and the y-array of multi-dimensional data can be declared from values in the data count by dimension array.

The function calc::OptIntp::setSpanX() completely fills the x-array with values along equally spaced intervals. Optimization can start from any point, however, equal spacing should be a good guess.

Next, the 3 different spline objects are defined, and associated with the curve generating function. The calc::EndCube object is used with calc::NCube to better define the curve. The bell-curve has zero slope at the origin, and this fact can be passed to the cubic spline using the enumerator calc::EndCube::DY (1st derivative). The upper end of the spline is given the default "natural spline" condition using the enumerator calc::EndCube::DDY (2nd derivative) and a zero value. The 3 interpolation optimizer are given increased calc::OptIntp::GRID_SAMPLE values to improve their curve-fit. Because there are 2 dimensions, 25 error samples will be taken to determine interval error in each interval of the interpolation grid.

Once all parameters have been defined, interpolation optimization can begin. The y-arrays of the interpolation objects are initially set to zero. A call to calc::OptIntp::update() automatically fills the array with the correct data. Error data is generated from the current state of the interpolation object with calls to calc::OptIntp::gridError() and calc::OptIntp::randomError(). This data generation/print out is performed before and after the optimization to show if error has decreased.

 #include <cstring>                                    // standard library memory functions
 #include "OptIntp.h"                                  // function interpolation objects
 #include "PrtC.h"                                     // auxiliary print stream object

 using namespace std;                                  // standard C++ library namespace
 using namespace calc;                                 // CalcLib namespace

 typedef EndCube<double,double> ECUd;                  // cubic spline derivative information type
 typedef OptIntp<double,double> OPTd;                  // interpolation optimization type

 // 2D Gauss curve generating function
 double curve(const double *xVec)
    {return exp(-(xVec[0]*xVec[0]+xVec[1]*xVec[1]));}

 // results print function
 void results(ostream &ostr, const char *cDat, double va, double vb, double vc)
    {ostr<<"   "<<cDat<<"  "<<va<<"   "<<vb<<"   "<<vc<<endl;}

 // error generate function
 bool errors(OPTd& opt)
    {return (opt.update()                              // check interpolation tables
           ||opt.gridError()                           // generate grid selection errors
           ||opt.randomError());}                      // generate random selection errors

 // summary print function
 bool summary(ostream &ostr, const char *cDat, OPTd& ocu, OPTd& onu, OPTd& opo)
 {
    // generate error data
    if((errors(ocu))||(errors(onu))||(errors(opo))) return true;

    // print error data
    ostr<<"   "<<cDat<<endl;
    ostr<<"   "<<"                          NCube      NNurb      NPoly"<<endl;
    results(ostr,"Grid Error Average:   ",ocu.getError(ocu.GRID_AVG  ),
             onu.getError(onu.GRID_AVG  ),opo.getError(opo.GRID_AVG  ));
    results(ostr,"Random Error Average: ",ocu.getError(ocu.RANDOM_AVG),
             onu.getError(onu.RANDOM_AVG),opo.getError(opo.RANDOM_AVG));
    ostr<<endl;

    return false;
 }

 int main(void)
 {
    // set parameters
    const int    NDIM=2;                               // number of interpolation dimensions
    const int    NCNT0=3;                              // data count of 1st dimension
    const int    NCNT1=3;                              // data count of 2nd dimension
    const int    ND[NDIM]={NCNT0,NCNT1};               // data count by dimension array
    const double XLO =0.0;                             // independent interpolation data low value
    const double XHI =1.0;                             // independent interpolation data high value
    const double XTOL=0.001;                           // independent data tolerance
    const double YTOL=0.001;                           // dependent data tolerance

    double aX0[NCNT0];                                 // independent data array (1st dimension)
    double aX1[NCNT1];                                 // independent data array (2nd dimension)
    double *aXn[NDIM]={aX0,aX1};                       // independent data by dimension array
    double aYn[NCNT1][NCNT0];                          // dependent data array (1st dimension)

    // clear data arrays
    std::memset(aX0,0,ND[0]*sizeof(*aX0));             // standard library function from <cstring>
    std::memset(aX1,0,ND[1]*sizeof(*aX1));             // set independent data arrays to zero's
    std::memset(aYn[0],0,ND[0]*ND[1]*sizeof(**aYn));   // set dependent data array to zero's

    // spread independent data array values at even intervals
    if(OPTd::setSpanX(aXn,XLO,XHI,ND,NDIM)) return 1;

    // create interpolation objects
    const ECUd EC(0.0,ECUd::DY,0.0,ECUd::DDY);         // spline lower end point should have zero slope
    const ECUd ECN[NDIM]={EC,EC};                      // required for both dimensions
    NCube<double,double> ncu(aXn,&aYn[0][0],ND,NDIM,ECN);
    NNurb<double,double> nnu=ncu;
    NPoly<double,double> npo=nnu;

    // create optimization objects
    OPTd ocu(curve,ncu); 
    OPTd onu(curve,nnu);
    OPTd opo(curve,npo);

    ocu.setParameter(ocu.GRID_SAMPLE,5);               // set 5 test samples per dimension
    onu.setParameter(onu.GRID_SAMPLE,5);
    opo.setParameter(opo.GRID_SAMPLE,5);

    // set up print stream
    cout.precision(5);
    cout.setf(ios::showpoint,ios::showpoint);

    // print errors, optimize, print errors
    cout<<endl<<"   calc::NIntp Class Example Application"<<endl<<endl;
    if((summary(cout,"Errors Before Optimization",ocu,onu,opo))
     ||(ocu.optimize(XTOL,YTOL))
     ||(onu.optimize(XTOL,YTOL))
     ||(opo.optimize(XTOL,YTOL))
     ||(summary(cout,"Errors After Optimization",ocu,onu,opo))) return 1;

    // print out interpolation object with best results (cubic spline)
    PrtC pr;
    pr.setFormat(pr.DOUBLE,pr.DOUBLE);
    pr.stream(cout,ncu);
    cout<<flush;

    return 0;
 }

The following output file was generated by the above source code. Three different interpolation methods are used: (1) calc::NCube (cubic splines); (2) calc::NNurb (non-uniform rational Bezier splines); (3) calc::NPoly (polynomial curve-fit splines). Each interpolation object is evaluated using the same data. From the output data, the calc::NCube cubic spline is optimized to have the lowest error.

   calc::NIntp Class Example Application

   Errors Before Optimization
                             NCube      NNurb      NPoly
   Grid Error Average:     0.0074914   0.020002   0.020002
   Random Error Average:   0.0065139   0.018117   0.018117

   Errors After Optimization
                             NCube      NNurb      NPoly
   Grid Error Average:     0.0047456   0.020002   0.020002
   Random Error Average:   0.0041324   0.018117   0.018117


   // CalcLib NIntp Multi-Dimensional Tabular Interpolation Parameters

   // CalcLib EndCube Derivative Identifier Parameters
   const EndCube< double,double > endCuben     [ 2]={
      EndCube< double,double >(double(  0.00000e+00), EndCube< double,double >::DY  , double(  0.00000e+00), EndCube< double,double >::DDY ),
      EndCube< double,double >(double(  0.00000e+00), EndCube< double,double >::DY  , double(  0.00000e+00), EndCube< double,double >::DDY ),
   };

   // CalcLib Data Count by Dimension Data Array
   const int           arNn         [ 2]={ int(            3), int(            3)};

   // CalcLib Independent Data Array
   const double        arX0         [ 3]={ double(  0.00000e+00), double(  6.88000e-01), double(  1.00000e+00)};
   const double        arX1         [ 3]={ double(  0.00000e+00), double(  6.89000e-01), double(  1.00000e+00)};
   const double*       arXn         [ 2]={ arX0, arX1};

   // CalcLib Dependent Data Array
   const double        arYn         [ 3][ 3]=
      {
         { double(  1.00000e+00), double(  6.22916e-01), double(  3.67879e-01)},
         { double(  6.22059e-01), double(  3.87490e-01), double(  2.28843e-01)},
         { double(  3.67879e-01), double(  2.29158e-01), double(  1.35335e-01)},
      };

   // CalcLib 2nd Derivative Data Array
   const double        arDn         [ 3][ 3]=
      {
         { double( -3.99741e-01), double(  2.84256e-03), double(  0.00000e+00)},
         { double( -2.48662e-01), double(  1.76824e-03), double(  0.00000e+00)},
         { double( -1.47056e-01), double(  1.04572e-03), double(  0.00000e+00)},
      };

---