#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" /** * 3D function. */ inline double f3(const double x, const double y, const double z) { static const double a = 1.0; static const double b = 1.0; static const double c = 1.0; return x * (a + x * (b + c*x)) * y * (a + y * (b + c*y)) * z * (a + z * (b + c*z)); } /** * \file * * Example program to test multi-dimensional interpolation. * \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 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); } 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 g3; for (double x = xmin; x < xmax + 0.5*dx; x += dx) { for (double y = xmin; y < xmax + 0.5*dx; y += dx) { for (double z = xmin; z < xmax + 0.5*dx; z += dx) { g3[x][y][z] = f3(x,y,z); } } } g3.compile(); g3.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero)); if (numberOfEvents > 0) { JTimer timer; JQuantile Q; for (int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(xmin, xmax); const double y = gRandom->Uniform(xmin, xmax); const double z = gRandom->Uniform(xmin, xmax); const double v = f3(x,y,z); timer.start(); const double w = g3(x,y,z); timer.stop(); Q.put(v - w); } if (debug >= debug_t) { Q.print(cout); timer.print(cout, 1.0 / numberOfEvents, micro_t); } ASSERT(Q.getMean() <= precision); ASSERT(Q.getSTDev() <= precision); } else { for ( ; ; ) { cout << "> " << flush; string buffer; if (getline(cin, buffer) && !buffer.empty()) { double x, y, z; istringstream(buffer) >> x >> y >> z; try { cout << f3(x,y,z) << ' ' << g3(x,y,z) << endl; } catch(const JException& exception) { cout << exception << endl; } } else { break; } } } return 0; }