source: Sophya/trunk/SophyaLib/NTools/integ.cc@ 3002

#include "sopnamsp.h"
#include "integ.h"
#include "generalfit.h"

// A faire :
//   1 dans GLInteg, optimiser recalcul en utilisant des poids
//     calcules sur [0,1] et chgt de variable
//   2 dans GLInteg, distinction nStep = subdivisions de l'intervalle,
//     et on applique GL d'ordre "ordre" sur chaque subdivision, en
//     utilisant le (1)

/*!
  \ingroup NTools
  \class SOPHYA::Integrator

  \brief Classe abstraite d'intégration numérique 1D.

  On fournit une fonction double f(double) au constructeur, ou
  une GeneralFunction avec des paramètres définis.
  L'objet Integrator est convertible en valeur double qui est la valeur
  de l'intégrale. Diverses méthodes permettent de choisir des options
  de calcul, et ces méthodes retournent une référence sur l'objet, pour
  permettre une notation chaînée.

  \sa TrpzInteg GLInteg
*/

//!  Constructeur par défaut. L'objet n'est pas utilisable en l'état.
Integrator::Integrator()
: mFunc(NULL), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{}

//! Constructeur à partir de la fonction double->double, et des bornes d'intégration.
Integrator::Integrator(FUNC f, double xmin, double xmax)
: mFunc(f), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(xmin), mXMax(xmax)
{}

//! Constructeur à partir de la fonction double->double, et des bornes d'intégration.
Integrator::Integrator(fun f, double xmin, double xmax)
: mFunc(new Function(f)), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(xmin), mXMax(xmax)
{}


/*!
  Constructeur sans spécifier les bornes. Elles sont positionnées
  à [0,1], et on pourra les modifier plus tard.
*/
Integrator::Integrator(FUNC f)
: mFunc(f), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{}

/*!
  Constructeur sans spécifier les bornes. Elles sont positionnées
  à [0,1], et on pourra les modifier plus tard.
*/
Integrator::Integrator(fun f)
: mFunc(new Function(f)), mGFF(NULL), mGFFParm(NULL),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{}

/*!
  Constructeur à partir d'une GeneralFunction. La fonction doit être une
  fonction de une variable, et on fournit les valeurs des paramètres ainsi
  que les bornes d'intégration.
*/
Integrator::Integrator(GeneralFunction* gff, double* par, double xmin, double xmax)
: mFunc(NULL), mGFF(gff), mGFFParm(par),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(xmin), mXMax(xmax)
{
  DBASSERT(gff->NVar() == 1);
}

/*!
  Constructeur à partir d'une GeneralFunction. La fonction doit être une
  fonction de une variable, et on fournit les valeurs des paramètres. 
  On ne spécifie pas les bornes. Elles sont positionnées
  à [0,1], et on pourra les modifier plus tard.
*/
Integrator::Integrator(GeneralFunction* gff, double* par)
: mFunc(NULL), mGFF(gff), mGFFParm(par),
  mNStep(50), mDX(-1), mReqPrec(-1),
  mXMin(0), mXMax(1)
{
  DBASSERT(gff->NVar() == 1);
}

Integrator::~Integrator()
{if(mFunc) delete mFunc;}

/*! 
  \brief Spécifie le nombre de pas pour l'intégration numérique.
  
  La signification peut dépendre de la méthode d'intégration.
*/
Integrator& 
Integrator::NStep(int n)
{
  mNStep = n;
  mDX = mReqPrec = -1;
  StepsChanged();
  return *this;
}

/*! 
  \brief Spécifie le nombre de pas pour l'intégration numérique.
  
  La signification peut dépendre de la méthode d'intégration.
*/
Integrator& 
Integrator::DX(double d)
{
  mDX = d;
  mNStep = mReqPrec = -1;
  StepsChanged();
  return *this;
}

/*!
  \brief Spécifie la précision souhaitée. 

  La signification peut dépendre de la méthode d'intégration.
  Non disponible dans toutes les méthodes d'intégration.
*/
Integrator& 
Integrator::ReqPrec(double p)
{
  DBASSERT( !"Pas encore implemente !");
  mReqPrec = p;
  mDX = mNStep = -1;
  StepsChanged();
  return *this;
}

//! Spécifie les bornes de l'intégrale.
Integrator& 
Integrator::Limits(double xmin, double xmax)
{
  mXMin = xmin;
  mXMax = xmax;
  LimitsChanged();
  return *this;
}

//! Spécifie la fonction à intégrer, sous forme double f(double).
Integrator& 
Integrator::Func(FUNC f)
{
  mFunc = f;
  mGFFParm = NULL;
  mFunc = NULL;
  FuncChanged();
  return *this;
}

/*!
  \brief Spécifie la fonction à intégrer, sous forme de GeneralFunction 
  à une variable, et les paramètres sont fournis.
*/
Integrator& 
Integrator::Func(GeneralFunction* gff, double* par)
{
  mGFF = gff;
  mGFFParm = par;
  mFunc = NULL;
  DBASSERT(gff->NVar() == 1);
  FuncChanged();
  return *this;
}

double
Integrator::FVal(double x) const
{
  DBASSERT( mFunc || (mGFF && mGFFParm) );
  return mFunc ? (*mFunc)(x) : mGFF->Value(&x, mGFFParm); 
  }

double
Integrator::ValueBetween(double xmin, double xmax)
{
  Limits(xmin,xmax);
  return Value();
}

void
Integrator::Print(int lp)
{
 cout<<"Integrator between "<<mXMin<<" and "<<mXMax
     <<" in "<<mNStep<<" steps"<<endl;
 if(lp>1) cout<<"mFunc="<<mFunc<<"mGFF="<<mGFF
              <<"mReqPrec="<<mReqPrec<<"mDX="<<mDX<<endl;
}

/*!
  \ingroup NTools
  \class SOPHYA::TrpzInteg

  \brief Implementation de Integrator par la methode des trapezes.
  
  Classe d'intégration par la méthode des trapèzes.
  Voir Integrator pour les méthodes. Le nombre de pas
  est le nombre de trapèze, le pas d'intégration est
  la largeur des trapèzez. Impossible de demander une
  précision.
  
  \sa SOPHYA::Integrator
*/


TrpzInteg::TrpzInteg()
{}

TrpzInteg::TrpzInteg(FUNC f, double xmin, double xmax)
: Integrator(f, xmin, xmax)
{}

TrpzInteg::TrpzInteg(fun f, double xmin, double xmax)
: Integrator(f, xmin, xmax)
{}

TrpzInteg::TrpzInteg(FUNC f)
: Integrator(f)
{}

TrpzInteg::TrpzInteg(fun f)
: Integrator(f)
{}

TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par, double xmin, double xmax)
: Integrator(gff, par, xmin, xmax)
{}

TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par)
: Integrator(gff, par)
{}

TrpzInteg::~TrpzInteg()
{}

//! Return the value of the integral.
double
TrpzInteg::Value()
{
  double dx = mDX;
  double nstep = mNStep;

  if (dx > 0)
    nstep = (mXMax - mXMin)/dx;

  if (nstep <= 0)
    nstep = 10;

  dx = (mXMax - mXMin) / nstep;

  double s = (FVal(mXMin) + FVal(mXMax))/2;
  double x = dx + mXMin;
  for (int i=1; i<nstep; i++, x += dx)
    s += FVal(x);
    
  return s * dx;
}



/*!
  \ingroup NTools
  \class SOPHYA::GLInteg
  
  \brief Implementation de Integrator par la methode de Gauss-Legendre.

  Classe d'intégration par la méthode de Gauss-Legendre.
  Voir Integrator pour les méthodes. 
  Pour le moment, nstep est l'ordre de la méthode.
  Il est prévu un jour de spécifier l'ordre, et que NStep
  découpe en intervalles sur chacun desquels on applique GL.
  Le principe de la méthode est de calculer les valeurs de la 
  fonction aux zéros des polynomes de Legendre. Avec les poids
  qui vont bien, GL d'ordre n est exacte pour des polynomes de 
  degré <= 2n+1 (monome le + haut x^(2*n-1).
  Impossible de demander une précision donnée.

  \sa SOPHYA::Integrator

  \warning statut EXPERIMENTAL , NON TESTE
*/



GLInteg::GLInteg()
: mXPos(NULL), mWeights(NULL)
{}


GLInteg::GLInteg(FUNC f, double xmin, double xmax)
: Integrator(f, xmin, xmax), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::GLInteg(fun f, double xmin, double xmax)
: Integrator(f, xmin, xmax), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::GLInteg(FUNC f)
: Integrator(f), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::GLInteg(fun f)
: Integrator(f), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::GLInteg(GeneralFunction* gff, double* par, double xmin, double xmax)
: Integrator(gff, par, xmin, xmax), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::GLInteg(GeneralFunction* gff, double* par)
: Integrator(gff, par), mXPos(NULL), mWeights(NULL), mOrder(8)
{
  NStep(1);
}

GLInteg::~GLInteg()
{
 if(mXPos)    delete[] mXPos;     mXPos    = NULL;
 if(mWeights) delete[] mWeights;  mWeights = NULL;
}

void
GLInteg::InvalWeights()
{
 if(mXPos)    delete[] mXPos;     mXPos    = NULL;
 if(mWeights) delete[] mWeights;  mWeights = NULL;
}


void
GLInteg::LimitsChanged()
{
 ComputeBounds();
}

void
GLInteg::StepsChanged()
{
 ComputeBounds();
}

GLInteg&
GLInteg::AddBound(double x)
{
 if (x<=mXMin || x>=mXMax) throw RangeCheckError("GLInteg::AddBound()") ; 
 // On introduira les classes d'exections apres reflexion et de maniere systematique (Rz+cmv)
 // if (x<=mXMin || x>=mXMax) throw range_error("GLInteg::AddBound  bound outside interval");
 mBounds.insert(x);
 return *this;
}


//! Definit l'ordre de la methode d'integration Gauus-Legendre
GLInteg&
GLInteg::SetOrder(int order)
{
 mOrder = order;
 InvalWeights();
 return *this;
}


//! Retourne la valeur de l'integrale
double
GLInteg::Value()
{
 if (!mXPos) ComputeWeights();
 double s = 0;
 set<double>::iterator i=mBounds.begin();
 set<double>::iterator j=i; j++;
 while (j != mBounds.end()) {
   double s1 = 0;
   double x1 = *i;
   double x2 = *j;
   for(int k=0; k<mOrder; k++)
      s1 += mWeights[k] * FVal(x1 + (x2-x1)*mXPos[k]);
   s += s1*(x2-x1);
   i++; j++;
 }
 return s;
}

void
GLInteg::ComputeBounds()
{
  mBounds.erase(mBounds.begin(), mBounds.end());
  for (int i=0; i<=mNStep; i++) 
    mBounds.insert(mXMin + (mXMax-mXMin)*i/mNStep);
}

void
GLInteg::ComputeWeights()
{
   const double EPS_gauleg = 5.0e-12;

   int m=(mOrder+1)/2;
   const double xxMin = 0;
   const double xxMax = 1;
   
   mXPos    = new double[mOrder];
   mWeights = new double[mOrder];
   
   double xm=0.5*(xxMax+xxMin);
   double xl=0.5*(xxMax-xxMin);
   for (int i=1;i<=m;i++)  {
      double z=cos(3.141592654*(i-0.25)/(mOrder+0.5));
      double z1, pp;
      do {
         double p1=1.0;
         double p2=0.0;
         for (int j=1;j<=mOrder;j++) {
            double p3=p2;
            p2=p1;
            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
         }
         pp=mOrder*(z*p1-p2)/(z*z-1.0);
         z1=z;
         z=z1-p1/pp;
      } while (fabs(z-z1) > EPS_gauleg);
      mXPos[i-1]         = xm-xl*z;
      mXPos[mOrder-i]    = xm+xl*z;
      mWeights[i-1]      = 2.0*xl/((1.0-z*z)*pp*pp);
      mWeights[mOrder-i] = mWeights[i-1];
   }

}

//! Imprime l'ordre et la valeur des poids sur cout
void
GLInteg::Print(int lp)
{
 Integrator::Print(lp);
 cout<<"GLInteg order="<<mOrder<<endl;
 if(lp>0 && mOrder>0) {
   for(int i=0;i<mOrder;i++)
     cout<<" ("<<mXPos[i]<<","<<mWeights[i]<<")";
   cout<<endl;
 }
}