#include #include #include #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::maplist JMaplist_t; typedef JMultiFunction 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); 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; }