JPolynome3G.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <map>

#include "JMath/JPolynome.hh"
#include "JTools/JMapList.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JToolsToolkit.hh"
#include "JTools/JFunctionalMap.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  using namespace JPP;


  /**
   * 3D polynomial function.
   */
  struct JFunction3D {

    /**
     * Result.
     */
    struct result_type {
      /**
       * Constructor.
       *
       * \param  __f         function value  f
       * \param  __fx        derivative     df/dx
       * \param  __fy        derivative     df/dy
       * \param  __fz        derivative     df/dz
       */
      result_type(const double __f,
                  const double __fx,
                  const double __fy,
                  const double __fz) :
        f (__f),
        fx(__fx),
        fy(__fy),
        fz(__fz)
      {}


      /**
       * Type conversion operator.
       *
       * \return             function value
       */
      operator double() const
      { 
        return f;
      }
    

      const double f;
      const double fx;
      const double fy;
      const double fz;
    };
  

    /**
     * Constructor.
     *
     * In each dimension, the function value and the derivative are given by: 
     * <pre>
     *      f (x) = f1.getValue(x);
     *      f'(x) = f1.getDerivative(x);
     * </pre>
     *
     * \param  f1          polynome
     */
    JFunction3D(const JPolynome& f1) :
      __f1(f1)
    {}


    /**
     * Evaluation of function value and derivatives.
     *
     * \param  x           x-value
     * \param  y           y-value
     * \param  z           z-value
     * \return             function and derivatives
     */
    inline result_type operator()(const double x,
                                  const double y,
                                  const double z) const
    {
      return result_type(f (x) * f (y) * f (z),
                         fp(x) * f (y) * f (z),
                         f (x) * fp(y) * f (z),
                         f (x) * f (y) * fp(z));
    }


    /**
     * Evaluation of integral.
     *
     * \param  xmin        minimal x-value
     * \param  xmax        maximal x-value
     * \param  ymin        minimal y-value
     * \param  ymax        maximal y-value
     * \param  zmin        minimal z-value
     * \param  zmax        maximal z-value
     * \return             integral
     */
    inline double operator()(const double xmin, const double xmax,
                             const double ymin, const double ymax,
                             const double zmin, const double zmax) const
    {
      return ((__f1.getIntegral(xmax) - __f1.getIntegral(xmin)) *
              (__f1.getIntegral(ymax) - __f1.getIntegral(ymin)) *
              (__f1.getIntegral(zmax) - __f1.getIntegral(zmin)));
    }


  protected:
    /**
     * Evaluation of function value.
     *
     * \param  x           x-value
     * \return             function value
     */
    inline double f(const double x) const
    { 
      return __f1.getValue(x);
    }


    /**
     * Evaluation of derivative.
     *
     * \param  x           x-value
     * \return             derivative
     */
    inline double fp(const double x) const
    { 
      return __f1.getDerivative(x);
    }

    JPolynome __f1;
  };


  /**
   * Test method for multi-dimensional interpolation with derivatives.
   */
  template<class JMultiFunction_t>
  int do_main(int argc, char **argv, const char* title)
  {
    using namespace std;

    int                  numberOfBins;
    JPolynome            f1;
    map<string, double>  precision;
    int                  debug;

    precision["spline"] = 1.0e-3;
    precision["polint"] = 1.0e-10;

    try {

      JParser<> zap("Example program to test multi-dimensional integration.");

      zap['N'] = make_field(numberOfBins)     = 11;
      zap['P'] = make_field(f1)               = JPolynome(3, 1.0, 1.0, 1.0);
      zap['e'] = make_field(precision)        = JPARSER::initialised();
      zap['d'] = make_field(debug)            = 3;

      zap(argc, argv);
    }
    catch(const exception &error) {
      FATAL(error.what() << endl);
    }


    if (numberOfBins <= 2) {
      FATAL("Number of bins " << numberOfBins << endl);
    }

    using namespace JPP;


    const JFunction3D f3(f1);

    const double xmin = -1.0;
    const double xmax = +1.0;
    const double dx   = (xmax - xmin) / (numberOfBins - 1);
   

    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();
    //g3.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero));

    if (debug >= debug_t) {
      cout << endl 
           << setw(8)  << left  << title        << ' ' 
           << setw(12) << right << "real"       << ' ' 
           << setw(12) << right << "calculated" << endl;
      cout << setw(8)  << left  << "integral";
      cout << ' ' << SCIENTIFIC(12,4) << f3(xmin, xmax, xmin, xmax, xmin, xmax);
      cout << ' ' << SCIENTIFIC(12,4) << getIntegral(g3);
      cout << endl;
    }

    ASSERT(fabs(f3(xmin, xmax, xmin, xmax, xmin, xmax) - getIntegral(g3)) <= precision[title] * f3(xmin, xmax, xmin, xmax, xmin, xmax));

    return 0;
  }
}


/**
 * \file
 *
 * Example program to test multi-dimensional integration.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  {
    typedef JGridPolint3Function1D_t                          JFunction1D_t;
    typedef JMAPLIST<JPolint3FunctionalGridMap,
                     JPolint3FunctionalGridMap>::maplist      JMaplist_t;
    typedef JMultiFunction<JFunction1D_t, JMaplist_t>         JMultiFunction_t;

    if (do_main<JMultiFunction_t>(argc, argv, "polint") != 0) return 1;
  }

  {
    typedef JGridSplineFunction1D_t                           JFunction1D_t;
    typedef JMAPLIST<JSplineFunctionalGridMap,
                     JSplineFunctionalGridMap>::maplist       JMaplist_t;
    typedef JMultiFunction<JFunction1D_t, JMaplist_t>         JMultiFunction_t;

    if (do_main<JMultiFunction_t>(argc, argv, "spline") != 0) return 1;
  }
}