#ifndef _utl_Integrator_h_
#define _utl_Integrator_h_

#include <vector>

#include <utl/PolynomialInterpolation.h>
#include <utl/ErrorLogger.h>


namespace utl {

  /**
    \class Integrator Integrator.h "utl/Integrator.h"
    \brief Class for integration of functions with one independent parameter.

    \author Darko Veberic
    \version $Id$
    \date 15 Mar 2007
    \ingroup math

    Example of usage:
    \code
      // define function through the functor class
      class SinFunctor {
      public:
        SinFunctor(const double frequency, const double amplitude)
          : fFrequency(frequency), fAmplitude(amplitude) { }

        double operator()(const double x) const
        { return fAmplitude * sin(fFrequency * x); }

      private:
        const double fFrequency;
        const double fAmplitude;
      };

      ...

      // one time integration
      SinFunctor sinf(3.1415, 2);
      const double intSin =  MakeIntegrator(sinf).GetIntegral(0, 1);

      ...

      // repeated use of the integrator
      double result = 0;
      SinFunctor sinf(3.1415, 10);
      const utl::Integrator<SinFunctor> integrator(sinf);
      for (double d = 0; d < 100; d += 3)
        result += integrator.GetIntegral(d, d+3);
    \endcode

    Simple functions that do not need initialization or additional parameters can be used
    directly (without the functor):
    \code
      double MyLog(const double x) { return std::log(x); }

      const double res = MakeIntegrator(MyLog).GetIntegral(1, 2);
    \endcode
    Note that some built-in function do not have a static binding (like log(), exp() etc)
    and have to be wrapped in functor or function.
  */

  template<class Functor>
  class Integrator {

  public:
    Integrator(Functor& functor, const double accuracy = 1e-5)
      : fFunctor(functor), fAccuracy(accuracy) { }

    /// final accuracy goal of the integration
    void SetAccuracy(const double accuracy) { fAccuracy = accuracy; }

    /// call this method when you do not care about the specific implementation
    double GetIntegral(const double a, const double b) const
    { return GetRombergIntegral(a, b); }


    /** \brief Romberg integration
      Setting order to 2 is equivalent to the Simpson's method.
    */
    double
    GetRombergIntegral(const double a, const double b,
                       const int order = 5, const int maxIterations = 20)
      const
    {
      const double delta = b - a;

      std::vector<double> tInt(maxIterations+1);
      std::vector<double> tStep(maxIterations+1);
      double oldInt = tInt[0] = GetBoxAverage(a, delta, 2);
      tStep[0] = delta;
      int i;
      for (i = 1; i < order; ++i) {
        oldInt = tInt[i] = GetTrapezoidalAverage(oldInt, a, delta, i);
        tStep[i] = delta/(1 << (2*i));  // h^2
      }

      const double eps = 0.5*fAccuracy;
      double extrapolation = 0;
      for (; i <= maxIterations; ++i) {
        oldInt = tInt[i] = GetTrapezoidalAverage(oldInt, a, delta, i);
        tStep[i] = delta/(1 << (2*i));  // extrapolate h^2 dependence for h->0
        const int first = i - order + 1;
        double residual;
        extrapolation =
          PolynomialInterpolation(order, &tStep[first], &tInt[first], 0, residual);
        if (fabs(residual) <= eps*fabs(extrapolation))
          return delta*extrapolation;
      }
      if (fabs(extrapolation) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*extrapolation;
    }


    /// Simpson improvement build on top of the trapezoidal integration
    double
    GetSimpsonIntegral(const double a, const double b,
                       const int minLevel = 4, const int maxLevel = 20)
      const
    {
      const double delta = b - a;
      double trapezAverage = GetBoxAverage(a, delta, 2);

      int level;
      for (level = 1; level < minLevel; ++level)
        trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level);

      double trapezOld = trapezAverage;
      trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level++);
      double simpsonOld = (4./3)*trapezAverage - (1./3)*trapezOld;

      for (; level <= maxLevel; ++level) {
        trapezOld = trapezAverage;
        trapezAverage = GetTrapezoidalAverage(trapezAverage, a, delta, level);
        const double simpsonAverage = (4./3)*trapezAverage - (1./3)*trapezOld;
        if (fabs(simpsonAverage - simpsonOld) < fAccuracy*fabs(simpsonAverage) ||
            (!simpsonAverage && !simpsonOld))
          return delta*simpsonAverage;
        simpsonOld = simpsonAverage;
      }
      if (fabs(simpsonOld) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*simpsonOld;
    }


    /// basic trapezoidal integration
    double
    GetTrapezoidalIntegral(const double a, const double b,
                           const int minLevel = 4, const int maxLevel = 20)
      const
    {
      const double delta = b - a;
      double average = GetBoxAverage(a, delta, 2);
      int level;
      for (level = 1; level <= minLevel; ++level)
        average = GetTrapezoidalAverage(average, a, delta, level);
      for (; level <= maxLevel; ++level) {
        const double old = average;
        average = GetTrapezoidalAverage(average, a, delta, level);
        if (fabs(average - old) < fAccuracy*fabs(average) || (!average && !old))
          return delta*average;
      }
      if (fabs(average) > fAccuracy)
        WARNING("Maximal level limit reached, result is not accurate.");
      return delta*average;
    }

  private:
    /**
      \brief Calculates succsessive approximations to the function average

      \param previousApproximation: result of the level-1 call
      \param a: start of the integration interval
      \param delta: width of the integration interval (b - a)
      \param level: should be greater or equal to 1 (see point scheme below)

      The evaluation points are generated according to the following
      scheme:

      \code
                         interval
      level  a                               b    n

        0    x                               x    2  <-- done by GetBoxAverage()
        1    |               x               |    1
        2    |       x               x       |    2
        3    |   x       x       x       x   |    4
        4    | x   x   x   x   x   x   x   x |    8
        5    |x x x x x x x x x x x x x x x x|   16

      x = function evaluation
      \endcode
    */

    double
    GetTrapezoidalAverage(const double previousApproximation,
                          const double a, const double delta, const int level)
      const
    {
      const int n = 1 << (level - 1);
      const double step = delta/n;
      return 0.5*(previousApproximation +
                  GetBoxAverage(a+0.5*step, step, n));
    }

    /// average of a function represented with equidistant boxes
    double
    GetBoxAverage(const double start, const double step, const int n)
      const
    {
      double sum = 0;
      for (int i = 0; i < n; ++i) {
        const double x = start + i*step;
        sum += fFunctor(x);
      }
      return sum/n;
    }

    Functor& fFunctor;
    double fAccuracy;

  };


  /// convenience factory
  template<class Functor>
  inline
  Integrator<Functor>
  MakeIntegrator(Functor& f)
  {
    return Integrator<Functor>(f);
  }
  

  template<class Functor, unsigned int dim>
  class VectorIntegrator {

  public:

    VectorIntegrator(Functor& functor, const double accuracy = 1e-5)
      : fFunctor(functor), fAccuracy(accuracy) { }

    /// final accuracy goal of the integration
    void SetAccuracy(const double accuracy) { fAccuracy = accuracy; }

    /// Romberg integration, for details see utl::Integrator
    void
    GetIntegral(double* result, const double a, const double b,
                const int order = 5, const int maxIterations = 20)
      const
    {
      for (unsigned int k = 0; k < dim; ++k) result[k] = 0;

      const double delta = b - a;

      std::vector<double> tInt[dim];
      for (unsigned int k = 0; k < dim; ++k) tInt[k].resize(maxIterations+1);
      std::vector<double> tStep(maxIterations+1);
      double oldInt[dim];
      GetBoxAverage(oldInt, a, delta, 2);
      for (unsigned int k = 0; k < dim; ++k) tInt[k][0] = oldInt[k];
      tStep[0] = delta;
      int i;
      for (i = 1; i < order; ++i) {
        GetTrapezoidalAverage(oldInt, a, delta, i);
        for (unsigned int k = 0; k < dim; ++k) tInt[k][i] = oldInt[k];
        tStep[i] = delta/(1 << (2*i));  // h^2
      }

      const double eps = 0.5*fAccuracy;
      for (unsigned int k = 0; k < dim; ++k) result[k] = 0;
      for (; i <= maxIterations; ++i) {
        GetTrapezoidalAverage(oldInt, a, delta, i);
        for (unsigned int k = 0; k < dim; ++k) tInt[k][i] = oldInt[k];
        tStep[i] = delta/(1 << (2*i));  // extrapolate h^2 dependence for h->0
        const int first = i - order + 1;

        unsigned int nAccepted = 0;
        for (unsigned int k = 0; k < dim; ++k)
        {
          double residual;
          result[k] =
            PolynomialInterpolation(order, &tStep[first], &tInt[k][first], 0, residual);
          if (fabs(residual) <= eps*fabs(result[k]))
            ++nAccepted;
        }
        if (nAccepted == dim)
        {
          for (unsigned int k = 0; k < dim; ++k) result[k] *= delta;
          return;
        }
      }

      WARNING("Maximal level limit reached, result is not accurate.");
      for (unsigned int k = 0; k < dim; ++k) result[k] *= delta;
    }

  private:

    void
    GetTrapezoidalAverage(double* previousApproximation,
                          const double a, const double delta, const int level)
      const
    {
      const int n = 1 << (level - 1);
      const double step = delta/n;
      double add[dim];
      GetBoxAverage(add, a+0.5*step, step, n);
      for (unsigned int k = 0; k < dim; ++k)
        previousApproximation[k] = 0.5*(previousApproximation[k] + add[k]);
    }

    /// average of a function represented with equidistant boxes
    void
    GetBoxAverage(double* sum, const double start, const double step, const int n)
      const
    {
      for (unsigned int k = 0; k < dim; ++k) sum[k] = 0;

      double result[dim];
      for (int i = 0; i < n; ++i) {
        const double x = start + i*step;
        fFunctor(result, x);
        for (unsigned int k = 0; k < dim; ++k) sum[k] += result[k];
      }

      for (unsigned int k = 0; k < dim; ++k) sum[k]/=n;
    }

    Functor& fFunctor;
    double fAccuracy;

  };

}


#endif