GSLMCIntegrator.cxx

Go to the documentation of this file.

// @(#)root/mathmore:$Id$
// Author: Magdalena Slawinska  08/2007

 /**********************************************************************
  *                                                                    *
  * Copyright (c) 2007 ROOT Foundation,  CERN/PH-SFT                   *
  *                                                                    *
  * This library is free software; you can redistribute it and/or      *
  * modify it under the terms of the GNU General Public License        *
  * as published by the Free Software Foundation; either version 2     *
  * of the License, or (at your option) any later version.             *
  *                                                                    *
  * This library is distributed in the hope that it will be useful,    *
  * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
  * General Public License for more details.                           *
  *                                                                    *
  * You should have received a copy of the GNU General Public License  *
  * along with this library (see file COPYING); if not, write          *
  * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
  * 330, Boston, MA 02111-1307 USA, or contact the author.             *
  *                                                                    *
  **********************************************************************/
//
// implementation file for class GSLMCIntegrator
// Author: Magdalena Slawinska
//
//

#include "Math/IFunction.h"
#include "Math/Error.h"

#include "GSLMonteFunctionWrapper.h"

#include "Math/GSLMCIntegrator.h"
#include "Math/GSLRndmEngines.h"
#include "GSLMCIntegrationWorkspace.h"
#include "GSLRngWrapper.h"

#include <algorithm>
#include <functional>
#include <ctype.h>   // need to use c version of tolower defined here


#include "gsl/gsl_monte_vegas.h"
#include "gsl/gsl_monte_miser.h"
#include "gsl/gsl_monte_plain.h"



namespace ROOT {
namespace Math {



// constructors

// GSLMCIntegrator::GSLMCIntegrator():
//    fResult(0),fError(0),fStatus(-1),
//    fWorkspace(0),
//    fFunction(0)
// {
//    // constructor of GSL MCIntegrator.Vegas MC is set as default integration type
//    //set random number generator
//    fRng = new GSLRngWrapper();
//    fRng->Allocate();
//    // use the default options
//    ROOT::Math::IntegratorMultiDimOptions opts;  // this create the default options
//    SetOptions(opts);
// }


GSLMCIntegrator::GSLMCIntegrator(MCIntegration::Type type, double absTol, double relTol, unsigned int calls):
   fType(type),
   fDim(0),
   fCalls((calls > 0)  ? calls : IntegratorMultiDimOptions::DefaultNCalls()),
   fAbsTol((absTol >= 0) ? absTol : IntegratorMultiDimOptions::DefaultAbsTolerance() ),
   fRelTol((relTol >= 0) ? relTol : IntegratorMultiDimOptions::DefaultRelTolerance() ),
   fResult(0),fError(0),fStatus(-1),
   fExtGen(false),
   fWorkspace(0),
   fFunction(0)
{
   // constructor of GSL MCIntegrator using enumeration as type
   SetType(type);
   //set random number generator
   fRng = new GSLRngWrapper();
   fRng->Allocate();
   // use the default options for the needed extra parameters
   // use the default options for the needed extra parameters
   if (fType == MCIntegration::kVEGAS) {
      IOptions * opts = IntegratorMultiDimOptions::FindDefault("VEGAS");
      if (opts != 0) SetParameters( VegasParameters(*opts) );
   }
   else  if (fType == MCIntegration::kMISER) {
      IOptions * opts = IntegratorMultiDimOptions::FindDefault("MISER");
      if (opts != 0)  SetParameters( MiserParameters(*opts) );
   }

}

GSLMCIntegrator::GSLMCIntegrator(const char * type, double absTol, double relTol, unsigned int calls):
   fDim(0),
   fCalls(calls),
   fAbsTol(absTol),
   fRelTol(relTol),
   fResult(0),fError(0),fStatus(-1),
   fExtGen(false),
   fWorkspace(0),
   fFunction(0)
{
   // constructor of GSL MCIntegrator. Vegas MC is set as default integration type if type == 0
   SetTypeName(type);

   //set random number generator
   fRng = new GSLRngWrapper();
   fRng->Allocate();
   // use the default options for the needed extra parameters
   if (fType == MCIntegration::kVEGAS) {
      IOptions * opts = IntegratorMultiDimOptions::FindDefault("VEGAS");
      if (opts != 0) SetParameters( VegasParameters(*opts) );
   }
   else  if (fType == MCIntegration::kMISER) {
      IOptions * opts = IntegratorMultiDimOptions::FindDefault("MISER");
      if (opts != 0)  SetParameters( MiserParameters(*opts) );
   }

}



GSLMCIntegrator::~GSLMCIntegrator()
{
   // delete workspace
   if (fWorkspace) delete fWorkspace;
   if (fRng != 0 && !fExtGen) delete fRng;
   if (fFunction != 0) delete fFunction;
   fRng = 0;

}


// disable copy ctrs


GSLMCIntegrator::GSLMCIntegrator(const GSLMCIntegrator &) :
   VirtualIntegratorMultiDim()
{}

GSLMCIntegrator & GSLMCIntegrator::operator=(const GSLMCIntegrator &) { return *this; }





void GSLMCIntegrator::SetFunction(const IMultiGenFunction &f)
{
   // method to set the a generic integration function
   if(fFunction == 0) fFunction = new  GSLMonteFunctionWrapper();
   fFunction->SetFunction(f);
   fDim = f.NDim();
}

void GSLMCIntegrator::SetFunction( GSLMonteFuncPointer f,  unsigned int dim, void * p  )
{
   // method to set the a generic integration function
   if(fFunction == 0) fFunction = new  GSLMonteFunctionWrapper();
   fFunction->SetFuncPointer( f );
   fFunction->SetParams ( p );
   fDim = dim;
}



double GSLMCIntegrator::Integral(const double* a, const double* b)
{
   // evaluate the Integral of a over the defined interval (a[],b[])
   assert(fRng != 0);
   gsl_rng* fr = fRng->Rng();
   assert(fr != 0);
   if (!CheckFunction()) return 0;

   // initialize by  creating the right WS
   // (if dimension and type are different than previous calculation)
   DoInitialize();

   if ( fType == MCIntegration::kVEGAS)
   {
      GSLVegasIntegrationWorkspace * ws = dynamic_cast<GSLVegasIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      fStatus = gsl_monte_vegas_integrate( fFunction->GetFunc(), (double *) a, (double*) b , fDim, fCalls, fr, ws->GetWS(),  &fResult, &fError);
   }
   else if (fType ==  MCIntegration::kMISER)
   {
      GSLMiserIntegrationWorkspace * ws = dynamic_cast<GSLMiserIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      fStatus = gsl_monte_miser_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(),  &fResult, &fError);
   }
   else if (fType ==  MCIntegration::kPLAIN)
   {
      GSLPlainIntegrationWorkspace * ws = dynamic_cast<GSLPlainIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      fStatus = gsl_monte_plain_integrate( fFunction->GetFunc(), (double *) a, (double *) b , fDim, fCalls, fr, ws->GetWS(),  &fResult, &fError);
   }
   /**/
   else
   {

      fResult = 0;
      fError = 0;
      fStatus = -1;
      std::cerr << "GSLIntegrator - Error: Unknown integration type" << std::endl;
      throw std::exception();
   }

   return fResult;

}


double GSLMCIntegrator::Integral(const GSLMonteFuncPointer & f, unsigned int dim, double* a, double* b, void * p )
{
   // evaluate the Integral for a function f over the defined interval (a[],b[])
   SetFunction(f,dim,p);
   return Integral(a,b);
}


/* to be added later
   double GSLMCIntegrator::Integral(GSLMonteFuncPointer f, void * p, double* a, double* b)
   {

   }

*/
//MCIntegration::Type GSLMCIntegrator::MCType() const {return fType;}

/**
   return  the Result of the last Integral calculation
*/
double GSLMCIntegrator::Result() const { return fResult; }

/**
   return the estimate of the absolute Error of the last Integral calculation
*/
double GSLMCIntegrator::Error() const { return fError; }

/**
   return the Error Status of the last Integral calculation
*/
int GSLMCIntegrator::Status() const { return fStatus; }


// setter for control Parameters  (getters are not needed so far )

/**
   set the desired relative Error
*/
void GSLMCIntegrator::SetRelTolerance(double relTol){ this->fRelTol = relTol; }

/**
   set the desired absolute Error
*/
void GSLMCIntegrator::SetAbsTolerance(double absTol){ this->fAbsTol = absTol; }

void GSLMCIntegrator::SetGenerator(GSLRandomEngine & r){
   // delete previous exist generator
   if (fRng && !fExtGen) delete fRng; 
   fRng = r.Engine();
   fExtGen = true; 
}

void GSLMCIntegrator::SetType (MCIntegration::Type type)
{
   // create workspace according to the type
   fType=type;
   if (fWorkspace != 0) {
      if (type == fWorkspace->Type() ) return;
      delete fWorkspace;  // delete because is a different type
      fWorkspace = 0;
   }
   //create Workspace according to type
   if (type == MCIntegration::kPLAIN) {
      fWorkspace = new  GSLPlainIntegrationWorkspace();
   }
   else if (type == MCIntegration::kMISER) {
      fWorkspace = new  GSLMiserIntegrationWorkspace();
   }
   else {
       // default: use  VEGAS
      if (type != MCIntegration::kVEGAS) {
         MATH_WARN_MSG("GSLMCIntegration","Invalid integration type : use Vegas as default");
         fType =  MCIntegration::kVEGAS;
      }
      fWorkspace = new  GSLVegasIntegrationWorkspace();
   }
}

void GSLMCIntegrator::SetTypeName(const char * type)
{
   // set the integration type using a string
   std::string typeName = (type!=0) ? type : "VEGAS";
   if (type == 0) MATH_INFO_MSG("GSLMCIntegration::SetTypeName","use default Vegas integrator method");
   std::transform(typeName.begin(), typeName.end(), typeName.begin(), (int(*)(int)) toupper );

   MCIntegration::Type integType =  MCIntegration::kVEGAS;  // default

   if (typeName == "PLAIN") {
      integType =  MCIntegration::kPLAIN;
   }
   else if (typeName == "MISER") {
      integType =  MCIntegration::kMISER;
   }
   else if (typeName != "VEGAS")  {
      MATH_WARN_MSG("GSLMCIntegration::SetTypeName","Invalid integration type : use Vegas as default");
   }

   // create the fWorkspace object
   if (integType != fType) SetType(integType);
}


void GSLMCIntegrator::SetMode(MCIntegration::Mode mode)
{
   //   set integration mode for VEGAS method
   if(fType ==  ROOT::Math::MCIntegration::kVEGAS)
   {
      GSLVegasIntegrationWorkspace * ws = dynamic_cast<GSLVegasIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      if(mode == MCIntegration::kIMPORTANCE) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE;
      else if(mode == MCIntegration::kSTRATIFIED) ws->GetWS()->mode = GSL_VEGAS_MODE_STRATIFIED;
      else if(mode == MCIntegration::kIMPORTANCE_ONLY) ws->GetWS()->mode = GSL_VEGAS_MODE_IMPORTANCE_ONLY;
   }

   else std::cerr << "Mode not matching integration type";
}

void GSLMCIntegrator::SetOptions(const ROOT::Math::IntegratorMultiDimOptions & opt)
{
   //   set integration options
   SetTypeName(opt.Integrator().c_str() );
   SetAbsTolerance( opt.AbsTolerance() );
   SetRelTolerance( opt.RelTolerance() );
   fCalls = opt.NCalls();

   //std::cout << fType << "   " <<  MCIntegration::kVEGAS << std::endl;

   // specific options
   ROOT::Math::IOptions * extraOpt = opt.ExtraOptions();
   if (extraOpt) {
      if (fType ==  MCIntegration::kVEGAS ) {
         VegasParameters p(*extraOpt);
         SetParameters(p);
      }
      else if (fType ==  MCIntegration::kMISER ) {
         MiserParameters p(fDim); // if possible pass dimension
         p = (*extraOpt);
         SetParameters(p);
      }
      else {
         MATH_WARN_MSG("GSLMCIntegrator::SetOptions","Invalid options set for the chosen integration type - ignore them");
      }
   }
}


void GSLMCIntegrator::SetParameters(const VegasParameters &p)
{
   // set method parameters
   if (fType ==  MCIntegration::kVEGAS)
   {
      GSLVegasIntegrationWorkspace * ws = dynamic_cast<GSLVegasIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      ws->SetParameters(p);
   }
   else
      MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
}

void GSLMCIntegrator::SetParameters(const MiserParameters &p)
{
   // set method parameters
   if (fType ==  MCIntegration::kMISER)
   {
      GSLMiserIntegrationWorkspace * ws = dynamic_cast<GSLMiserIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      ws->SetParameters(p);
   }
   else
      MATH_ERROR_MSG("GSLIntegrator::SetParameters"," Parameters not matching integration type");
}


void GSLMCIntegrator::DoInitialize ( )
{
   //    initialize by setting  integration type

   if (fWorkspace == 0) return;
   if (fDim == fWorkspace->NDim() && fType == fWorkspace->Type() )
      return; // can use previously existing ws

   // otherwise clear workspace
   fWorkspace->Clear();
   // and create a new one
   fWorkspace->Init(fDim);
}



//----------- methods specific for VEGAS

double GSLMCIntegrator::Sigma()
{
   // returns the error sigma from the last iteration of the VEGAS algorithm
   if(fType == MCIntegration::kVEGAS)
   {
      GSLVegasIntegrationWorkspace * ws = dynamic_cast<GSLVegasIntegrationWorkspace *>(fWorkspace);
      assert (ws != 0);
      return ws->GetWS()->sigma;
   }
   else
   {
      std::cerr << "Parameter not mathcing integration type";
      return 0;
   }

}


/**
*/
double GSLMCIntegrator::ChiSqr()
{
   //   returns chi-squared per degree of freedom for the estimate of the integral
   if(fType == MCIntegration::kVEGAS)
   {
      GSLVegasIntegrationWorkspace * ws = dynamic_cast<GSLVegasIntegrationWorkspace *>(fWorkspace);
      assert(ws != 0);
      return ws->GetWS()->chisq;
   }
   else
   {
      std::cerr << "Parameter not mathcing integration type";
      return 0;
   }
}



bool GSLMCIntegrator::CheckFunction()
{
   // internal method to check validity of GSL function pointer

   if (fFunction && fFunction->GetFunc() ) return true; 
   MATH_ERROR_MSG("GSLMCIntegrator::CheckFunction","Function has not been specified");
   return false;
}

const char * GSLMCIntegrator::GetTypeName() const {
   if (fType == MCIntegration::kVEGAS) return "VEGAS";
   if (fType == MCIntegration::kMISER) return "MISER";
   if (fType == MCIntegration::kPLAIN) return "PLAIN";
   return "UNDEFINED";
}

ROOT::Math::IntegratorMultiDimOptions  GSLMCIntegrator::Options() const {
   IOptions * extraOpts = ExtraOptions();
   ROOT::Math::IntegratorMultiDimOptions opt(extraOpts);
   opt.SetAbsTolerance(fAbsTol);
   opt.SetRelTolerance(fRelTol);
   opt.SetNCalls(fCalls);
   opt.SetWKSize(0);
   opt.SetIntegrator(GetTypeName() );
   return opt;
}

ROOT::Math::IOptions *  GSLMCIntegrator::ExtraOptions() const {
   if (!fWorkspace) return 0;
   return fWorkspace->Options();
}


} // namespace Math
} // namespace ROOT