#include #include #include #include "TRandom3.h" #include "JTools/JTable2D.hh" #include "JTools/JGrid.hh" #include "JTools/JElement.hh" #include "JTools/JPolint.hh" #include "JTools/JGridCollection.hh" #include "JTools/JMapList.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JQuantile.hh" #include "Jeep/JPrint.hh" #include "Jeep/JTimer.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * 3D function. * * \param x x value * \param y y value * \param z z value * \return function value */ inline double f3(const double x, const double y, const double z) { const double u = 1.0e2 / (x*x); return exp(-0.5*(y*y)/(u*u) - 0.5*(z*z)/(u*u)); } } /** * \file * * Example program to test interpolation between 2D tables. * \author mdejong */ int main(int argc, char **argv) { using namespace std; int numberOfEvents; int debug; try { JParser<> zap("Example program to test interpolation between 2D tables."); zap['n'] = make_field(numberOfEvents); zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } using namespace JPP; const int NX = 21; const int NY = 11; const int NZ = 11; const JGrid X(NX, -10.0, +10.0); const JGrid Y(NY, -5.0, +5.0); const JGrid Z(NZ, -5.0, +5.0); typedef JTable2D JTable2D_t; typedef JElement2D JElement2D_t; typedef JPolintFunction1D<2, JElement2D_t, JGridCollection> JFunction1D_t; JFunction1D_t buffer; for (int i = 0; i != X.getSize(); ++i) { DEBUG("table[" << i << "]" << endl); const double x = X.getX(i); JTable2D_t& table = buffer[x]; for (int __i = 0; __i != NY; ++__i) { for (int __j = 0; __j != NZ; ++__j) { const double y = Y.getX(__i); const double z = Z.getX(__j); table[__i][__j] = f3(x,y,z); DEBUG(' ' << FIXED(5,3) << table[__i][__j]); } DEBUG(endl); } } // interpolator between table cells typedef JMapList JMaplist_t; typedef JMultiFunction JMultiFunction_t; JMultiFunction_t g2; g2.configure(Y); for (JMultiFunction_t::iterator i = g2.begin(); i != g2.end(); ++i) { i->getY().configure(Z); } JTimer t1("table"); JTimer t2("polint"); JQuantile f1("table"); JQuantile f2("polint"); for (int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(X.getXmin(), X.getXmax()); const double y = gRandom->Uniform(Y.getXmin(), Y.getXmax()); const double z = gRandom->Uniform(Z.getXmin(), Z.getXmax()); const double v = f3(x,y,z); // simple lookup table t1.start(); const int __i = (int) (Y.getSize() * (y - X.getXmin()) / (Y.getXmax() - Y.getXmin())); const int __j = (int) (Z.getSize() * (z - Z.getXmin()) / (Z.getXmax() - Z.getXmin())); const double w1 = buffer(x)[__i][__j]; t1.stop(); // polynomial interpolation t2.start(); const JTable2D_t& table = buffer(x); for (int __i = 0; __i != Y.getSize(); ++__i) { for (int __j = 0; __j != Z.getSize(); ++__j) { g2.getY(__i).getY(__j) = table[__i][__j]; } } g2.compile(); const double w2 = g2(y,z); t2.stop(); f1.put(v - w1); f2.put(v - w2); } for (JQuantile buffer[] = { f1, f2, JQuantile() }, *i = buffer; i->getCount() != 0; ++i) { cout << i->getTitle() << endl; cout << "mean " << SCIENTIFIC(10,2) << i->getMean() << endl; cout << "RMS " << SCIENTIFIC(10,2) << i->getSTDev() << endl; } for (JTimer buffer[] = { t1, t2, JTimer() }, *i = buffer; i->getTitle() != ""; ++i) { i->print(cout, true, micro_t); } }