JResultHesse3D.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

#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<class JKey_t, class JValue_t, class JDistance_t>
  struct JMap_t :
    public JPolintMap<N, JKey_t, JValue_t, JGridMap,
                      JResultHesse<typename JResultType<JValue_t>::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<N,
                            JElement2D<double, double>,
                            JGridCollection,
                            JResultHesse<double> >          JFunction1D_t;
  typedef JMAPLIST<JMap_t,
                   JMap_t>::maplist                         JMaplist_t;

  typedef JMultiFunction<JFunction1D_t, JMaplist_t>         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;
}