JPolynome5D.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "TRandom3.h"
 
#include "JTools/JMapList.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JQuantile.hh"
 
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * 5D function.
 */
inline double f5(const double x0,
                 const double x1,
                 const double x2,
                 const double x3,
                 const double x4)
{
  static const double a = 1.0;
  static const double b = 1.0;
  static const double c = 1.0;
 
  return 
    x0 * (a + x0 * (b + c*x0))  *  
    x1 * (a + x1 * (b + c*x1))  *
    x2 * (a + x2 * (b + c*x2))  *
    x3 * (a + x3 * (b + c*x3))  *
    x4 * (a + x4 * (b + c*x4));
}
 
 
/**
 * \file
 * Example program to test multi-dimensional interpolation.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
 
  int         numberOfEvents;
  int         numberOfBins;
  double      precision;
  int         debug;
 
  try {
 
    JParser<> zap("Example program to test multi-dimensional interpolation.");
 
    zap['n'] = make_field(numberOfEvents)   = 1000;
    zap['N'] = make_field(numberOfBins)     = 11;
    zap['e'] = make_field(precision)        = 1.0e-14;
    zap['d'] = make_field(debug)            = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  using namespace JPP;
 
 
  const double xmin = -1.0;
  const double xmax = +1.0;
  const double dx   = (xmax - xmin) / (numberOfBins - 1);
 
 
  typedef JGridPolint3Function1D_t                          JFunction1D_t;
 
  typedef JMAPLIST<JPolint3FunctionalGridMap,
                   JPolint3FunctionalGridMap,
                   JPolint3FunctionalGridMap,
                   JPolint3FunctionalGridMap>::maplist      JMaplist_t;
 
  typedef JMultiFunction<JFunction1D_t, JMaplist_t>         JMultiFunction_t;
 
    
  JMultiFunction_t g5;
 
   
  for (double x0 = xmin; x0 <= xmax + 0.5*dx; x0 += dx) {
    for (double x1 = xmin; x1 <= xmax + 0.5*dx; x1 += dx) {
      for (double x2 = xmin; x2 <= xmax + 0.5*dx; x2 += dx) {
        for (double x3 = xmin; x3 <= xmax + 0.5*dx; x3 += dx) {
          for (double x4 = xmin; x4 <= xmax + 0.5*dx; x4 += dx) {
            g5[x0][x1][x2][x3][x4] = f5(x0, x1, x2, x3, x4);
          }
        }
      }
    }
  }
    
  g5.compile();
  g5.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));
 
  ASSERT(numberOfEvents > 0);
 
  {
    const double x = xmin - dx;
 
    ASSERT(fabs(g5(x, 0.0, 0.0, 0.0, 0.0)) <= precision);
    ASSERT(fabs(g5(0.0, x, 0.0, 0.0, 0.0)) <= precision);
    ASSERT(fabs(g5(0.0, 0.0, x, 0.0, 0.0)) <= precision);
    ASSERT(fabs(g5(0.0, 0.0, 0.0, x, 0.0)) <= precision);
    ASSERT(fabs(g5(0.0, 0.0, 0.0, 0.0, x)) <= precision);
  }
 
  JTimer timer;
 
  JQuantile Q;
    
  for (int i = 0; i != numberOfEvents; ++i) {
    
    const double x0  = gRandom->Uniform(xmin, xmax);
    const double x1  = gRandom->Uniform(xmin, xmax);
    const double x2  = gRandom->Uniform(xmin, xmax);
    const double x3  = gRandom->Uniform(xmin, xmax);
    const double x4  = gRandom->Uniform(xmin, xmax);
      
    const double v   = f5(x0,x1,x2,x3,x4);
 
    timer.start();
 
    const double w   = g5(x0,x1,x2,x3,x4);
 
    timer.stop();
 
    Q.put(v - w);
  }
 
  if (debug >- debug_t) {
 
    Q.print(cout);
    
    timer.print(cout, 1.0 / numberOfEvents);
  }
 
  ASSERT(Q.getMean()  <= precision);
  ASSERT(Q.getSTDev() <= precision);
  
  return 0;
}