JTable2D.cc

Go to the documentation of this file.

 
#include <iostream>
#include <iomanip>
#include <cmath>
 
#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<double> X(NX, -10.0, +10.0);
  const JGrid<double> Y(NY,  -5.0,  +5.0);
  const JGrid<double> Z(NZ,  -5.0,  +5.0);
 
 
  typedef JTable2D<NY, NZ>                                         JTable2D_t;
  typedef JElement2D<double, JTable2D_t>                           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<JPolint1FunctionalGridMap>                      JMaplist_t;
  typedef JMultiFunction<JGridPolint1Function1D_t, JMaplist_t>     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);
  }
}