#include #include #include #include #include #include "TRandom3.h" #include "JLang/JException.hh" #include "JMath/JPolynome.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JGrid.hh" #include "JTools/JQuantile.hh" #include "JTools/JToolsToolkit.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using JTOOLS::JQuantile; struct JABC { JABC(const char* title) : f (title), fp(title), v (title) {} JQuantile f; JQuantile fp; JQuantile v; }; struct JPrecision_t { JPrecision_t() : f (0.0), fp(0.0), v (0.0) {} JPrecision_t(const double f, const double fp, const double v) : f (f), fp(fp), v (v) {} friend std::istream& operator>>(std::istream& in, JPrecision_t& object) { return in >> object.f >> object.fp >> object.v; } friend std::ostream& operator<<(std::ostream& out, const JPrecision_t& object) { return out << object.f << ' ' << object.fp << ' ' << object.v; } double f; double fp; double v; }; } /** * \file * * Example program to test 1D interpolation of a polynome. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef map map_type; unsigned int numberOfEvents; unsigned int numberOfBins; JPolynome f1; map_type precision; int debug; precision["spline"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4); precision["hermite"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4); precision["grid"] = JPrecision_t(1.0e-3, 1.0e-1, 1.0e-4); precision["polint"] = JPrecision_t(1.0e-14, 1.0e-12, 1.0e-14); try { JParser<> zap("Example program to test 1D interpolation of a polynome."); zap['n'] = make_field(numberOfEvents) = 0; zap['N'] = make_field(numberOfBins) = 21; zap['P'] = make_field(f1); zap['e'] = make_field(precision) = JPARSER::initialised(); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } const int N = 3; // fixed degree of polynomial interpolation JSplineFunction1S_t spline; JHermiteSplineFunction1S_t hermite; JGridSplineFunction1S_t grid; JPolintFunction1H_t< N > polint; JPolintFunction1D_t Polint; // integrand of polint const double xmin = -1.0; const double xmax = +1.0; spline .configure(make_grid(numberOfBins, xmin, xmax), f1); hermite.configure(make_grid(numberOfBins, xmin, xmax), f1); grid .configure(make_grid(numberOfBins, xmin, xmax), f1); polint .configure(make_grid(numberOfBins, xmin, xmax), f1); spline .compile(); hermite.compile(); grid .compile(); polint .compile(); integrate(polint, Polint); Polint.compile(); if (numberOfEvents != 0) { JABC Qspline ("spline"); JABC Qhermite("hermite"); JABC Qgrid ("grid"); JABC Qpolint ("polint"); for (unsigned int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(xmin, xmax); const double y1 = f1.getValue(x); const double fp = f1.getDerivative(x); const double F1 = f1.getIntegral(x) - f1.getIntegral(xmin); Qspline .f .put(y1 - spline (x).f); Qhermite.f .put(y1 - hermite(x).f); Qgrid .f .put(y1 - grid (x).f); Qpolint .f .put(y1 - polint (x).f); Qspline .fp.put(fp - spline (x).fp); Qhermite.fp.put(fp - hermite(x).fp); Qgrid .fp.put(fp - grid (x).fp); Qpolint .fp.put(fp - polint (x).fp); Qspline .v .put(F1 - spline (x).v); Qhermite.v .put(F1 - hermite(x).v); Qgrid .v .put(F1 - grid (x).v); Qpolint .v .put(F1 - Polint (x)); } if (debug >= debug_t) { cout << "RMS: f f' F " << endl; cout << "_________________________________________" << endl; cout << "spline " << SCIENTIFIC(10,3) << Qspline .f .getSTDev() << ' ' << SCIENTIFIC(10,3) << Qspline .fp.getSTDev() << ' ' << SCIENTIFIC(10,3) << Qspline .v .getSTDev() << endl; cout << "hermite " << SCIENTIFIC(10,3) << Qhermite.f .getSTDev() << ' ' << SCIENTIFIC(10,3) << Qhermite.fp.getSTDev() << ' ' << SCIENTIFIC(10,3) << Qhermite.v .getSTDev() << endl; cout << "grid " << SCIENTIFIC(10,3) << Qgrid .f .getSTDev() << ' ' << SCIENTIFIC(10,3) << Qgrid .fp.getSTDev() << ' ' << SCIENTIFIC(10,3) << Qgrid .v .getSTDev() << endl; cout << "polint " << SCIENTIFIC(10,3) << Qpolint .f .getSTDev() << ' ' << SCIENTIFIC(10,3) << Qpolint .fp.getSTDev() << ' ' << SCIENTIFIC(10,3) << Qpolint .v .getSTDev() << endl; } ASSERT(Qspline .f .getSTDev() <= precision["spline"].f); ASSERT(Qspline .fp.getSTDev() <= precision["spline"].fp); ASSERT(Qspline .v .getSTDev() <= precision["spline"].v); ASSERT(Qhermite.f .getSTDev() <= precision["hermite"].f); ASSERT(Qhermite.fp.getSTDev() <= precision["hermite"].fp); ASSERT(Qhermite.v .getSTDev() <= precision["hermite"].v); ASSERT(Qgrid .f .getSTDev() <= precision["grid"].f); ASSERT(Qgrid .fp.getSTDev() <= precision["grid"].fp); ASSERT(Qgrid .v .getSTDev() <= precision["grid"].v); ASSERT(Qpolint .f .getSTDev() <= precision["polint"].f); ASSERT(Qpolint .fp.getSTDev() <= precision["polint"].fp); ASSERT(Qpolint .v .getSTDev() <= precision["polint"].v); } else { for ( ; ; ) { string buffer; cout << "> " << flush; getline(cin, buffer); if (buffer != "") { try { double x; istringstream(buffer) >> x; cout << "f1 " << f1 (x) << endl; cout << "spline " << spline (x).f << endl; cout << "hermite " << hermite(x).f << endl; cout << "grid spline " << grid (x).f << endl; cout << "polynomial " << polint (x).f << endl; } catch(const JException& exception) { cout << exception.what() << endl; } } else { break; } } } return 0; }