#include #include #include #include #include "TRandom3.h" #include "TMath.h" #include "JTools/JMapList.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JConstantFunction1D.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JQuantile.hh" #include "JTools/JToolsToolkit.hh" #include "Jeep/JTimer.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using namespace JPP; /** * 5D test function. */ inline double f5(const double x0, const double x1, const double x2, const double x3, const double x4) { return TMath::Gaus(x4, x0+x1+x2+x3, 1.0, kTRUE); } } /** * \file * * Example program to create 1D function interpolator from multi-dimensional interpolator. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; int numberOfEvents; int numberOfBins; double precision; int debug; try { JParser<> zap("Example program to create 1D function interpolator from multi-dimensional interpolator."); 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); } ASSERT(numberOfEvents > 0); const double xmin = -1.0; const double xmax = +1.0; const double dx = (xmax - xmin) / (numberOfBins - 1); typedef JSplineFunction1D, JCollection> JFunction1D_t; typedef JFunction1D_t::abscissa_type abscissa_type; typedef JFunction1D_t::value_type value_type; typedef JConstantFunction1D > JCollection1D_t; typedef JMAPLIST::maplist JMaplist_t; JMultiFunction g5; JMultiFunction h5; 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); h5[x0][x1][x2][x3][x4] = f5(x0,x1,x2,x3,x4); } } } } } g5.compile(); h5.compile(); const double x0 = +0.15; const double x1 = -0.25; const double x2 = 0.25; const double x3 = -0.15; JFunction1D_t g1; JTimer timer("4D interpolator"); timer.start(); for (int i = 0; i != 100; ++i) { copy(h5(x0,x1,x2,x3), g1); // put interpolated data into 1D functional object } timer.stop(); timer.print(cout, 1.0/100, micro_t); g1.compile(); JTimer t1("1D interpolator"); JTimer t5("5D interpolator"); JQuantile Q1("1D interpolator"); JQuantile Q5("5D interpolator"); // mini buffer to reduce overhead in timer const int N = 100; double x[N]; double y[N]; double v[N]; double w[N]; for (int i = 0; i != numberOfEvents; ++i) { for (int __i = 0; __i != N; ++__i) { const double x4 = gRandom->Uniform(xmin, xmax); x[__i] = x4; y[__i] = f5(x0,x1,x2,x3,x4); } t1.start(); for (int __i = 0; __i != N; ++__i) { v[__i] = g1(x[__i]); } t1.stop(); t5.start(); for (int __i = 0; __i != N; ++__i) { w[__i] = g5(x0,x1,x2,x3,x[__i]); } t5.stop(); for (int __i = 0; __i != N; ++__i) { Q1.put(v[__i] - y[__i]); Q5.put(w[__i] - y[__i]); } } if (debug >= notice_t) { Q1.print(cout); Q5.print(cout, false); t1.print(cout, 1.0 / (N * numberOfEvents), nano_t); t5.print(cout, 1.0 / (N * numberOfEvents), nano_t); } ASSERT(fabs(Q1.getMean() - Q5.getMean()) <= precision); ASSERT(fabs(Q1.getSTDev() - Q5.getSTDev()) <= precision); return 0; }