/**
 * Class with different functional needed for taking decisions.
 * based on the code from http://rosettacode.org/wiki/Verify_distribution_uniformity/Chi-squared_test
 *
 **@author I.Kulakov <I.Kulakov@gsi.de>
 */


#ifndef _TDHelper_h_
#define _TDHelper_h_

#include <iostream>
using namespace std;

#include <cmath>

template< typename T >
class TFunctionBase {
 public:
  TFunctionBase(){};
  ~TFunctionBase(){};

  virtual T operator()( T t ) const = 0;
};

template< typename T >
class TDHelper { //  TakingDecisionsHelper
 public:
  TDHelper(){};
  ~TDHelper(){};

    /// chi2 distribution for @ndf  degrees of freedom.
    /// return probability that x == @chi2 - chi2 distribution for @ndf  degrees of freedom.
  static T Chi2Probability( int ndf, T chi2 );
    /// integrated chi2 distribution for @ndf  degrees of freedom.
    /// return probability that x >= @chi2 in chi2-distribution with @ndf degrees of freedom.
  static T Chi2IProbability( int ndf, T chi2 );

    /// integrated gaussian distribution.
   /// return probability that @d >= x >= -@d in normal-distribution with @ndf degrees of freedom.
  static T GaussIProbability( T d );

  
 private:
  
    // Here are the functions required to compute the Chi-Squared probability. These are not needed if a library containing the necessary function is availabile (e.g. see Numerical Integration, Gamma function).

    /// Numerical integration method
    /// return int( f(x), x = a .. b );
  static const T PI;
  static T Simpson3_8( const TFunctionBase<T>& f, T a, T b, int N);

  
    // --  for chi2 --
  static T Gamma_Spouge( T z );

  class GICore : public TFunctionBase<T> { // base function for GammaIncomplete_Q
   public:
    GICore( T a ): fA(a){};
    ~GICore(){};

    virtual T operator()( T t ) const;
   private:
    T fA; // parameter
  };

  static T GammaIncomplete_Q( T a, T x );

      // --  for gaussian --
  class GaussCore : public TFunctionBase<T> { // base function for GammaIncomplete_Q
   public:

    virtual T operator()( T t ) const;
   private:
    static const T A; // 1/sqrt(2*PI);
  };
  

    // forbid copying
  TDHelper(const TDHelper&);
  TDHelper operator=(const TDHelper&);
}; // class TDHelper

template< typename T > const T
TDHelper<T>::PI = static_cast<float>(3.14159265358979323846);

template< typename T > const T
TDHelper<T>::GaussCore::A = static_cast<float>(0.39894228040143267794);









#include <cmath>


  /// Numerical integration method
  /// return int( f(x), x = a .. b );
template< typename T >
T TDHelper<T>::Simpson3_8( const TFunctionBase<T>& f, T a, T b, int N)
{
  int j;
  T l1;
  T h = (b-a)/N;
  T h1 = h/3.0;
  T sum = f(a) + f(b);
 
  for (j=3*N-1; j>0; j--) {
    l1 = (j%3)? 3.0 : 2.0;
    sum += l1*f(a+h1*j) ;
  }
  return h*sum*0.125;
}

  // ---------  Chi2 -----------

template< typename T >
T TDHelper<T>::Gamma_Spouge( T z )
{
  const int A = 12;

  int k;
  static T cspace[A];
  static T *coefs = NULL; // TODO rid of pointer. TODO rid of static
  T accum;
  T a = A;
 
  if (!coefs) {
    T k1_factrl = 1.0;
    coefs = cspace;
    coefs[0] = sqrt(2.0*PI);
    for(k=1; k<A; k++) {
      coefs[k] = exp(a-k) * pow(a-k,k-0.5) / k1_factrl;
      k1_factrl *= -k;
    }
  }
 
  accum = coefs[0];
  for (k=1; k<A; k++) {
    accum += coefs[k]/(z+k);
  }
  accum *= exp(-(z+a)) * pow(z+a, z+0.5);
  return accum/z;
}

template< typename T >
T TDHelper<T>::GICore::operator()( T t ) const
{
  return pow(t, fA) * exp(-t); 
}

template< typename T >
T TDHelper<T>::GammaIncomplete_Q( T a, T x ) 
{
  const T h = 1.5e-2;  // approximate integration step size

  GICore core(a-1);
  if ( a == 0.5 ) {
    const T max = 20; // tail is too small
    T y = ( x < max ) ? x : max;
    return Simpson3_8( core, y, max, static_cast<int>((max-y)/h) )/Gamma_Spouge(a);
  } else {
      // this cuts off the tail of the integration to speed things up
    T y = a-1;
    while((core(y) * (x-y) > 2.0e-8) && (y < x))   y += .4;
    if (y>x) y=x;
    return 1.0 - Simpson3_8( core, 0, y, static_cast<int>(y/h) )/Gamma_Spouge(a);
  }
}

  /// return probability that x >= @chi2 in chi2-distribution with @ndf degrees of freedom.
template< typename T >
T TDHelper<T>::Chi2IProbability( int ndf, T chi2 )
{
  return GammaIncomplete_Q( 0.5*ndf, 0.5*chi2);
}

  /// chi2 distribution for @ndf  degrees of freedom.
  /// return probability that x == @chi2 - chi2 distribution for @ndf  degrees of freedom.
template< typename T >
T TDHelper<T>::Chi2Probability( int ndf, T chi2 ){
  const T r = ndf;
  return 0.5 * pow(chi2/2, r/2-1) * exp(-chi2/2) / Gamma_Spouge(r/2) ; 
}

  // ---------  Gaussian -----------

template< typename T >
T TDHelper<T>::GaussCore::operator()( T t ) const
{
  return exp( -t*t*0.5 ) * A; 
}

   /// return probability that @d >= x >= -@d in normal-distribution with @ndf degrees of freedom.
template< typename T >
T TDHelper<T>::GaussIProbability( T d ){
  const T h = 1.5e-2;  // approximate integration step size
  return 2*Simpson3_8( GaussCore(), 0, d, static_cast<int>(d/h) );
}

#endif // _TDHelper_h_