#include #include #include "TRandom3.h" #include "JMath/JPolynome.hh" #include "JTools/JElement.hh" #include "JTools/JGridCollection.hh" #include "JTools/JGridMap.hh" #include "JTools/JPolint.hh" #include "JTools/JMapList.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JQuantile.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { const int N = 3; using namespace JPP; JPolynome fx; JPolynome fy; JPolynome fz; /** * 3D function. */ inline double f3(const double x, const double y, const double z) { return fx(x) * fy(y) * fz(z); } template struct JMap_t : public JPolintMap::result_type>, JDistance_t> {}; } /** * \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['X'] = make_field(fx) = JPolynome(3, 1.0, 1.0, 1.0); zap['Y'] = make_field(fy) = JPolynome(3, 1.0, 1.0, 1.0); zap['Z'] = make_field(fz) = JPolynome(3, 1.0, 1.0, 1.0); zap['e'] = make_field(precision) = 1.0e-13; 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 JPolintFunction1D, JGridCollection, JResultHesse > JFunction1D_t; typedef JMAPLIST::maplist JMaplist_t; typedef JMultiFunction JMultiFunction_t; typedef JMultiFunction_t::result_type result_type; 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(); if (numberOfEvents > 0) { JQuantile Q; JQuantile H[N][N]; 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 result_type result = g3(x,y,z); const double v = f3(x,y,z); const double w = get_value(result); Q.put(v - w); H[0][0].put(fx.getDerivative().getDerivative(x) * fy(y) * fz(z) - result.fpp.f.f); H[1][0].put(fx.getDerivative(x) * fy.getDerivative(y) * fz(z) - result.fp.fp.f); H[2][0].put(fx.getDerivative(x) * fy(y) * fz.getDerivative(z) - result.fp.f.fp); H[1][1].put(fx(x) * fy.getDerivative().getDerivative(y) * fz(z) - result.f.fpp.f); H[2][1].put(fx(x) * fy.getDerivative(y) * fz.getDerivative(z) - result.f.fp.fp); H[2][2].put(fx(x) * fy(y) * fz.getDerivative().getDerivative(z) - result.f.f.fpp); } if (debug >= debug_t) { Q.print(cout); cout << "Hessian matrix" << endl; for (int i = 0; i != N; ++i) { for (int j = 0; j <= i; ++j) { cout << ' ' << SCIENTIFIC(12,3) << H[i][j].getMean(); } cout << endl; } for (int i = 0; i != N; ++i) { for (int j = 0; j <= i; ++j) { cout << ' ' << SCIENTIFIC(12,3) << H[i][j].getSTDev(); } cout << endl; } } ASSERT(Q.getMean() <= precision); ASSERT(Q.getSTDev() <= precision); for (int i = 0; i != N; ++i) { for (int j = 0; j <= i; ++j) { ASSERT(H[i][j].getMean() <= precision); ASSERT(H[i][j].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) << ' ' << get_value(g3(x,y,z)) << endl; } catch(const JException& exception) { cout << exception << endl; } } else { break; } } } return 0; }