JTable1D.cc

Go to the documentation of this file.

 
#include <iostream>
#include <iomanip>
#include <cmath>
 
#include "TRandom3.h"
 
#include "JTools/JTable1D.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/JStats.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * 2D function.
   * 
   * \param  x       x value
   * \param  y       y value
   * \return         function value
   */
  inline double f2(const double x,
                   const double y)
  {
    const double u = 1.0e2 / (x*x);
    
    return exp(-0.5*(y*y)/(u*u));
  }
}
 
 
/**
 * \file
 *
 * Example program to test interpolation between 1D tables.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  int  numberOfEvents;
  int  debug;
 
  try {
 
    JParser<> zap("Example program to test JTab1D & interpolation between 1D tables.");
 
    zap['n'] = make_field(numberOfEvents);
    zap['d'] = make_field(debug)            = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const int NX = 21;
  const int NY = 11;
 
  const JGrid<double> X(NX, -10.0, +10.0);
  const JGrid<double> Y(NY,  -5.0,  +5.0);
 
  typedef JTable1D<NY>                                             JTable1D_t;
  typedef JElement2D<double, JTable1D_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);
 
    JTable1D_t& table = buffer[x];
 
    for (int __i = 0; __i != NY; ++__i) {
 
      const double y = Y.getX(__i);
 
      table[__i] = f2(x,y);
 
      DEBUG(' ' << FIXED(5,3) << table[__i] << endl);
 
    }
  }
 
  // Test 1: basic operation on table
  const int NT = 3;
  JTable1D<NT> tab1;
  tab1[0] = 1; tab1[1] = 2; tab1[2] = 3;
  JTable1D<NT> tab2;
  tab2[0] = 4; tab2[1] = 5; tab2[2] = 6;
 
  ASSERT(tab1.getNX()         == NT, "getNX");
  ASSERT(tab2[1]              == 5,  "GetEntry");
  ASSERT(tab2.negate()[2]     == -6, "Negate");
  ASSERT((tab2.add(tab1))[0]  == -3, "Add");
  ASSERT((tab2.sub(tab1))[0]  == -4, "Subtract");
  ASSERT(tab1.mul(5)[1]       == 10, "Multiply");
  ASSERT(tab1.div(2)[1]       == 5,  "Divide");
  
  // Test 2: performance of interpolator between table cells
  JGridPolint1Function1D_t g2;
 
  g2.configure(Y);
 
  JTimer    t1("table");
  JTimer    t2("polint");
 
  JStats fa("table");
  JStats fb("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 v = f2(x,y);
 
 
    // simple lookup table
 
    t1.start();
 
    const int __i = (int) (Y.getSize() * (y - Y.getXmin()) / (Y.getXmax() - Y.getXmin()));
 
    const double w1 = buffer(x)[__i];
 
    t1.stop();
 
 
    // polynomial interpolation
 
    t2.start();
 
    const JTable1D_t& table = buffer(x);
 
    for (int __i = 0; __i != Y.getSize(); ++__i) {
      g2.getY(__i) = table[__i];
 
    }
 
    g2.compile();
 
    const double w2 = g2(y);
    t2.stop();
 
    fa.put(v - w1);
    fb.put(v - w2);
  }
 
  
  for (JStats buffer[] = { fa, fb, JStats() }, *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);
  }
}