Add new class DistSamplerOptions

  add method DistSampler::Init passing the options

 add new class TFoamSampler implementing the DistSampler interface using TFoam

 fix correct number of function calls when using the Fumili-like methods


git-svn-id: http://root.cern.ch/svn/root/trunk@37232 27541ba8-7e3a-0410-8455-c3a389f83636

math/foam/inc/LinkDef.h

@@ -9,5 +9,6 @@
 #pragma link C++ class TFoamVect+;
 #pragma link C++ class TFoamCell+;
 #pragma link C++ class TFoam+;
+#pragma link C++ class TFoamSampler+;
 
 #endif

math/foam/inc/TFoamSampler.h

+// @(#)root/mathcore:$Id$
+// Author: L. Moneta Fri Sep 22 15:06:47 2006
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ *                                                                    *
+ **********************************************************************/
+// Header file for class TFoamSampler
+
+#ifndef ROOT_TFoamSampler
+#define ROOT_TFoamSampler
+
+
+#ifndef ROOT_Math_DistSampler
+#include "Math/DistSampler.h"
+#endif
+
+
+namespace ROOT { 
+   
+   namespace Fit { 
+
+      class DataRange; 
+      class BinData;
+      class UnBinData; 
+   }
+
+   namespace Math { 
+   }
+}
+
+class TFoamIntegrand; 
+
+//_______________________________________________________________________________
+/**
+   TFoamSampler class
+   class implementing  the ROOT::Math::DistSampler interface using FOAM
+   for sampling arbitrary distributions. 
+
+
+*/
+
+class TRandom; 
+class TF1; 
+class TFoam; 
+
+class TFoamSampler : public ROOT::Math::DistSampler { 
+
+public: 
+
+   /// default constructor 
+   TFoamSampler(); 
+
+
+   /// virtual destructor 
+   virtual ~TFoamSampler();
+
+
+   using DistSampler::SetFunction; 
+
+   /// set the parent function distribution to use for random sampling (one dim case)
+   void SetFunction(const ROOT::Math::IGenFunction & func)  { 
+      fFunc1D = &func; 
+      SetFunction<const ROOT::Math::IGenFunction>(func, 1);
+   }
+
+   /// set the Function using a TF1 pointer
+   void SetFunction(TF1 * pdf); 
+
+
+   /**
+      initialize the generators with the default options
+   */
+   bool Init(const char * ); 
+
+   /**
+      initialize the generators with the fiven options
+   */
+   bool Init(const ROOT::Math::DistSamplerOptions & opt ); 
+
+   /** 
+       Set the random engine to be used 
+       Needs to be called before Init to have effect
+   */ 
+   void SetRandom(TRandom * r);
+
+   /** 
+       Set the random seed for the TRandom instances used by the sampler 
+       classes
+       Needs to be called before Init to have effect      
+   */ 
+   void SetSeed(unsigned int seed);
+
+//    /* 
+//       set the mode
+//     */
+//    void SetMode(double mode) { 
+//       fMode = mode; 
+//       fHasMode = true; 
+//    }
+
+//    /*
+//      set the area 
+//     */
+//    void SetArea(double area) { 
+//       fArea = area;
+//       fHasArea = true; 
+//    }
+
+   /**
+      Get the random engine used by the sampler
+    */
+   TRandom * GetRandom();
+
+   /**
+      sample one event in one dimension
+      better implementation could be provided by the derived classes   
+   */
+//    double Sample1D();//  { 
+// //       return fFoam->Sample();
+// //    }
+
+   /**
+      sample one event in multi-dimension by filling the given array
+      return false if sampling failed
+   */
+   bool Sample(double * x);
+
+   /**
+      sample one bin given an estimated of the pdf in the bin
+      (this can be function value at the center or its integral in the bin 
+      divided by the bin width)
+      By default do not do random sample, just return the function values
+    */
+   bool SampleBin(double prob, double & value, double *error = 0);
+
+
+
+protected: 
+
+   
+private: 
+
+   // private member
+   bool                              fOneDim;      // flag to indicate if the function is 1 dimension
+//    bool                              fHasMode;     // flag to indicate if a mode is set
+//    bool                              fHasArea;     // flag to indicate if a area is set
+//    double                            fMode;        // mode of dist
+//    double                            fArea;        // area of dist
+   const ROOT::Math::IGenFunction *  fFunc1D;      // 1D function pointer
+   TFoam *                           fFoam;        // foam engine class
+   TFoamIntegrand *                  fFoamDist;    // foam distribution interface  
+
+   //ClassDef(TFoamSampler,1)  //Distribution sampler class based on UNU.RAN
+
+};
+
+
+
+#endif /* ROOT_TFoamSampler */

math/foam/src/TFoamSampler.cxx

+// @(#)root/unuran:$Id$
+// Authors: L. Moneta,  Dec 2010
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2010  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ *                                                                    *
+ **********************************************************************/
+
+// Implementation file for class TFoamSampler
+
+#include "TFoamSampler.h"
+#include "Math/DistSamplerOptions.h"
+
+#include "TFoam.h"
+#include "TFoamIntegrand.h"
+#include "Math/OneDimFunctionAdapter.h"
+#include "Math/IOptions.h"
+#include "Fit/DataRange.h"
+
+#include "TRandom.h"
+#include "TError.h"
+
+#include "TF1.h"
+#include <cassert>
+#include <cmath>
+
+class FoamDistribution : public TFoamIntegrand { 
+
+public:
+
+   FoamDistribution(const ROOT::Math::IMultiGenFunction & f, const ROOT::Fit::DataRange & range) : 
+      fFunc(f), 
+      fX(std::vector<double>(f.NDim() ) ), 
+      fMinX(std::vector<double>(f.NDim() ) ),
+      fDeltaX(std::vector<double>(f.NDim() ) )
+   {
+      assert(f.NDim() == range.NDim() );
+      std::vector<double> xmax(f.NDim() );
+      for (unsigned int i = 0; i < range.NDim(); ++i) {
+         if (range.Size(i) == 0)
+            Error("FoamDistribution","Range is not set for coordinate dim %d",i);
+         else if (range.Size(i)>1) 
+            Warning("FoamDistribution","Using only first range in coordinate dim %d",i);
+
+         std::pair<double,double> r = range(i); 
+         fMinX[i] = r.first;
+         fDeltaX[i] = r.second - r.first; 
+      }
+   } 
+   // in principle function does not need to be cloned
+
+   virtual double Density(int ndim, double * x) {
+      assert(ndim == (int) fFunc.NDim() );
+      for (int i = 0; i < ndim; ++i)
+         fX[i] = fMinX[i] + x[i] * fDeltaX[i]; 
+
+      return (fFunc)(&fX[0]);
+   }
+
+   double  MinX(unsigned int i) { return fMinX[i]; }
+   double  DeltaX(unsigned int i) { return fDeltaX[i]; }
+              
+private:
+
+   const ROOT::Math::IMultiGenFunction & fFunc;
+   std::vector<double> fX; 
+   std::vector<double> fMinX; 
+   std::vector<double> fDeltaX; 
+
+};
+
+
+
+ClassImp(TFoamSampler)
+
+TFoamSampler::TFoamSampler() : ROOT::Math::DistSampler(), 
+   fOneDim(false), 
+//    fDiscrete(false),
+//    fHasMode(false), fHasArea(false),
+//    fMode(0), fArea(0),
+   fFunc1D(0),
+   fFoam(new TFoam("FOAM")  ),
+   fFoamDist(0)
+{}
+
+TFoamSampler::~TFoamSampler() {
+   assert(fFoam != 0);
+   delete fFoam; 
+   if (fFoamDist) delete fFoamDist; 
+}
+
+bool TFoamSampler::Init(const char *) { 
+
+   // initialize using default options
+   ROOT::Math::DistSamplerOptions opt; 
+   return Init(opt);
+}
+
+bool TFoamSampler::Init(const ROOT::Math::DistSamplerOptions & opt) { 
+   // initialize foam classes using the given algorithm
+   assert (fFoam != 0 );
+   if (NDim() == 0)  {
+      Error("TFoamSampler::Init","Distribution function has not been set ! Need to call SetFunction first.");
+      return false;
+   }
+
+   // initialize the foam 
+   fFoam->SetkDim(NDim() );
+
+   // initialize random number 
+   if (!GetRandom()) SetRandom(gRandom);
+
+   // create TFoamIntegrand class 
+   if (fFoamDist) delete fFoamDist; 
+   fFoamDist = new FoamDistribution(ParentPdf(),PdfRange());
+
+   fFoam->SetRho(fFoamDist);
+
+   // get extra options 
+   ROOT::Math::IOptions * fopt = opt.ExtraOptions(); 
+   if (fopt) { 
+      int nval = 0; 
+      if (fopt->GetIntValue("nCells", nval) ) fFoam->SetnCells(nval);
+      if (fopt->GetIntValue("nCell1D", nval) ) fFoam->SetnCells(nval);
+      if (fopt->GetIntValue("nCell2D", nval) ) fFoam->SetnCells(nval);
+      if (fopt->GetIntValue("nCell3D", nval) ) fFoam->SetnCells(nval);
+      if (fopt->GetIntValue("nCellND", nval) ) fFoam->SetnCells(nval);
+
+      if (fopt->GetIntValue("nSample", nval) ) fFoam->SetnSampl(nval);
+
+      if (fopt->GetIntValue("chatLevel", nval) ) fFoam->SetChat(nval);
+   }
+   fFoam->Initialize();
+
+   return true;
+      
+}
+
+
+void TFoamSampler::SetFunction(TF1 * pdf) { 
+   // set function from a TF1 pointer 
+   SetFunction<TF1>(*pdf, pdf->GetNdim());
+} 
+
+void TFoamSampler::SetRandom(TRandom * r) { 
+   // set random generator (must be called before Init to have effect)
+   fFoam->SetPseRan(r); 
+} 
+
+void TFoamSampler::SetSeed(unsigned int seed) { 
+   // set random generator seed (must be called before Init to have effect)
+   TRandom * r = fFoam->GetPseRan();
+   if (r) r->SetSeed(seed);
+} 
+
+TRandom * TFoamSampler::GetRandom() { 
+   // get random generator used 
+   return  fFoam->GetPseRan(); 
+} 
+
+// double TFoamSampler::Sample1D() { 
+//    // sample 1D distributions
+//    return (fDiscrete) ? (double) fFoam->SampleDiscr() : fFoam->Sample(); 
+// }
+
+bool TFoamSampler::Sample(double * x) { 
+   // sample multi-dim distributions
+
+   fFoam->MakeEvent();
+   fFoam->GetMCvect(x);
+   // adjust for the range   
+   for (unsigned int i = 0; i < NDim(); ++i) 
+      x[i] = ( (FoamDistribution*)fFoamDist)->MinX(i) + ( ( (FoamDistribution*) fFoamDist)->DeltaX(i))*x[i];
+
+   return true; 
+} 
+
+
+bool TFoamSampler::SampleBin(double prob, double & value, double *error) {
+   // sample a bin according to Poisson statistics
+
+   TRandom * r = GetRandom(); 
+   if (!r) return false; 
+   value = r->Poisson(prob); 
+   if (error) *error = std::sqrt(value);
+   return true; 
+}

math/mathcore/Module.mk

@@ -61,6 +61,7 @@ MATHCOREDH2  := $(MODDIRI)/TRandom.h \
                 $(MODDIRI)/Math/BrentMinimizer1D.h \
                 $(MODDIRI)/Math/BrentRootFinder.h \
                 $(MODDIRI)/Math/DistSampler.h \
+                $(MODDIRI)/Math/DistSamplerOptions.h \
                 $(MODDIRI)/Math/GoFTest.h \
                 $(MODDIRI)/Math/SpecFuncMathCore.h \
                 $(MODDIRI)/Math/DistFuncMathCore.h

math/mathcore/inc/Fit/Chi2FCN.h

@@ -134,6 +134,7 @@ public:
 
    /// i-th chi-square residual  
    virtual double DataElement(const double * x, unsigned int i, double * g) const { 
+      if (i==0) this->UpdateNCalls();
       return FitUtil::EvaluateChi2Residual(fFunc, fData, x, i, g); 
    }
 

math/mathcore/inc/Fit/LogLikelihoodFCN.h

@@ -105,6 +105,7 @@ public:
 
    /// i-th likelihood contribution and its gradient
    virtual double DataElement(const double * x, unsigned int i, double * g) const { 
+      if (i==0) this->UpdateNCalls();
       return FitUtil::EvaluatePdf(fFunc, fData, x, i, g); 
    }
 

math/mathcore/inc/Fit/PoissonLikelihoodFCN.h

@@ -101,6 +101,7 @@ public:
 
    /// i-th likelihood element and its gradient
    virtual double DataElement(const double * x, unsigned int i, double * g) const {
+      if (i==0) this->UpdateNCalls();
       return FitUtil::EvaluatePoissonBinPdf(fFunc, fData, x, i, g);
    }
 

math/mathcore/inc/LinkDef2.h

@@ -124,6 +124,7 @@
 #pragma link C++ class ROOT::Math::BrentMinimizer1D+;
 
 #pragma link C++ class ROOT::Math::DistSampler+;
+#pragma link C++ class ROOT::Math::DistSamplerOptions+;
 #pragma link C++ class ROOT::Math::GoFTest+;
 
 #include "LinkDef_Func.h"

math/mathcore/inc/Math/DistSampler.h

@@ -37,6 +37,7 @@ namespace ROOT {
 
    namespace Math { 
 
+      class DistSamplerOptions;
 
 /**
    @defgroup Random Random number generators and generation of random number distributions
@@ -99,6 +100,16 @@ public:
    */
    virtual bool Init(const char * /* algorithm */) { return true;}
 
+   /**
+      initialize the generators with the given option
+      which my include the algorithm but also more if 
+      the method is re-impelmented by derived class
+      The default implementation calls the above method 
+      passing just the algorithm name
+   */
+   virtual bool Init(const DistSamplerOptions & opt );
+
+
    /** 
        Set the random engine to be used 
        To be implemented by the derived classes who provides 

math/mathcore/inc/Math/DistSamplerOptions.h

+// @(#)root/mathcore:$Id$ 
+// Author: L. Moneta Fri Aug 15 2008
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2008  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ *                                                                    *
+ **********************************************************************/
+
+#ifndef ROOT_Math_DistSamplerOptions
+#define ROOT_Math_DistSamplerOptions
+
+#include <string>
+
+#include <iostream>
+
+namespace ROOT { 
+   
+
+namespace Math { 
+
+
+class IOptions;      
+
+//_______________________________________________________________________________
+/** 
+    DistSampler options 
+
+    @ingroup NumAlgo
+*/
+class DistSamplerOptions {
+
+public:
+
+   // static methods for setting and retrieving the default options 
+
+   static void SetDefaultSampler(const char * type, const char * algo = 0);
+   static void SetDefaultAlgorithm(const char * algo );
+   static void SetDefaultPrintLevel(int level);
+
+   static const std::string & DefaultSampler();
+   static const std::string & DefaultAlgorithm(); 
+   static int DefaultPrintLevel(); 
+
+   /// retrieve extra options - if not existing create a IOptions 
+   static ROOT::Math::IOptions & Default(const char * name);
+
+   // find extra options - return 0 if not existing 
+   static ROOT::Math::IOptions * FindDefault(const char * name);
+
+   /// print all the default options for the name given
+   static void PrintDefault(const char * name = 0, std::ostream & os = std::cout); 
+
+public:
+
+   // constructor using the default options 
+   // pass optionally a pointer to the additional options
+   // otehrwise look if they exist for this default minimizer
+   // and in that case they are copied in the constructed instance
+   DistSamplerOptions(IOptions * extraOpts = 0);
+
+   // destructor  
+   ~DistSamplerOptions();
+
+   // copy constructor 
+   DistSamplerOptions(const DistSamplerOptions & opt);
+
+   /// assignment operators 
+   DistSamplerOptions & operator=(const DistSamplerOptions & opt);
+
+   /** non-static methods for  retrivieng options */
+
+   /// set print level
+   int PrintLevel() const { return fLevel; }
+
+   /// return extra options (NULL pointer if they are not present)
+   IOptions * ExtraOptions() const { return fExtraOptions; }
+
+   /// type of minimizer
+   const std::string & Sampler() const { return fSamplerType; }
+
+   /// type of algorithm
+   const std::string & Algorithm() const { return fAlgoType; }
+
+   /// print all the options 
+   void Print(std::ostream & os = std::cout) const; 
+
+   /** non-static methods for setting options */
+
+   /// set print level
+   void SetPrintLevel(int level) { fLevel = level; }
+
+   /// set minimizer type
+   void SetSampler(const char * type) { fSamplerType = type; }
+
+   /// set minimizer algorithm
+   void SetAlgorithm(const char *type) { fAlgoType = type; }
+
+   /// set extra options (in this case pointer is cloned)
+   void  SetExtraOptions(const IOptions & opt); 
+
+
+private:
+
+   int fLevel;               // debug print level 
+   std::string fSamplerType;   // DistSampler type (Minuit, Minuit2, etc..
+   std::string fAlgoType;    // DistSampler algorithmic specification (Migrad, Minimize, ...)
+
+   // extra options
+   ROOT::Math::IOptions *   fExtraOptions;  // extra options 
+     
+};
+
+   } // end namespace Math
+
+} // end namespace ROOT
+
+#endif

math/mathcore/src/DistSampler.cxx

@@ -11,6 +11,7 @@
 // implementation file for class DistSampler
 
 #include "Math/DistSampler.h"
+#include "Math/DistSamplerOptions.h"
 #include "Math/Error.h"
 
 #include "Math/IFunction.h"
@@ -30,6 +31,11 @@ DistSampler::~DistSampler() {
    if (fRange) delete fRange;
 }
 
+bool DistSampler::Init(const DistSamplerOptions & opt ) { 
+   // default initialization with algorithm name
+   return Init(opt.Algorithm().c_str() );
+}
+
 void DistSampler::SetRange(double xmin, double xmax, int icoord) {
    if (!fRange) { 
       MATH_ERROR_MSG("DistSampler::SetRange","Need to set function before setting the range");

math/mathcore/src/DistSamplerOptions.cxx

+// @(#)root/mathcore:$Id$ 
+// Author: L. Moneta Fri Aug 15 2008
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2008  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ *                                                                    *
+ **********************************************************************/
+
+#include "Math/DistSamplerOptions.h"
+
+#include "Math/GenAlgoOptions.h"
+
+// case of using ROOT plug-in manager
+#ifndef MATH_NO_PLUGIN_MANAGER
+#include "TEnv.h"
+#endif 
+
+#include <iomanip>
+
+namespace ROOT { 
+   
+
+namespace Math { 
+
+   namespace Sampler { 
+      static std::string gDefaultSampler = "Unuran"; // take from /etc/system.rootrc in ROOT Fitter
+      static std::string gDefaultAlgorithm = "";
+      static int  gDefaultPrintLevel  = 0; 
+   }
+
+
+void DistSamplerOptions::SetDefaultSampler(const char * type, const char * algo ) {   
+   // set the default minimizer type and algorithm
+   if (type) Sampler::gDefaultSampler = std::string(type); 
+   if (algo) Sampler::gDefaultAlgorithm = std::string(algo);
+}
+void DistSamplerOptions::SetDefaultAlgorithm(const char * algo ) {   
+   // set the default minimizer type and algorithm
+   if (algo) Sampler::gDefaultAlgorithm = std::string(algo);
+}
+void DistSamplerOptions::SetDefaultPrintLevel(int level) {
+   // set the default printing level 
+   Sampler::gDefaultPrintLevel = level; 
+}
+
+const std::string & DistSamplerOptions::DefaultAlgorithm() { return Sampler::gDefaultAlgorithm; }
+int    DistSamplerOptions::DefaultPrintLevel()       { return Sampler::gDefaultPrintLevel; }
+
+const std::string & DistSamplerOptions::DefaultSampler() 
+{ 
+   // return default minimizer
+   // if is "" (no default is set) read from etc/system.rootrc
+   // use form /etc/ ??
+
+   return Sampler::gDefaultSampler; 
+}
+
+
+DistSamplerOptions::DistSamplerOptions(IOptions * extraOpts): 
+   fLevel( Sampler::gDefaultPrintLevel),
+   fExtraOptions(extraOpts)
+{
+   // constructor using  the default options
+
+   fSamplerType = DistSamplerOptions::DefaultSampler();
+
+   fAlgoType =  DistSamplerOptions::DefaultAlgorithm();
+
+
+   // check if extra options exists (copy them if needed)
+   if (!fExtraOptions) { 
+      IOptions * gopts = FindDefault( fSamplerType.c_str() );
+      if (gopts) fExtraOptions = gopts->Clone();
+   }
+}
+
+
+DistSamplerOptions::DistSamplerOptions(const DistSamplerOptions & opt) : fExtraOptions(0) {  
+   // copy constructor 
+   (*this) = opt; 
+}
+
+DistSamplerOptions & DistSamplerOptions::operator=(const DistSamplerOptions & opt) {  
+   // assignment operator 
+   if (this == &opt) return *this; // self assignment
+   fLevel = opt.fLevel;
+   fSamplerType = opt.fSamplerType; 
+   fAlgoType = opt.fAlgoType; 
+
+   if (fExtraOptions) delete fExtraOptions; 
+   fExtraOptions = 0; 
+   if (opt.fExtraOptions)  fExtraOptions =  (opt.fExtraOptions)->Clone();
+   return *this;
+}
+
+DistSamplerOptions::~DistSamplerOptions() { 
+   if (fExtraOptions) delete fExtraOptions; 
+}
+
+void DistSamplerOptions::SetExtraOptions(const IOptions & opt) {  
+   // set extra options (clone the passed one)
+   if (fExtraOptions) delete fExtraOptions; 
+   fExtraOptions = opt.Clone(); 
+}
+
+void DistSamplerOptions::Print(std::ostream & os) const {
+   //print all the options
+   os << std::setw(25) << "DistSampler Type"        << " : " << std::setw(15) << fSamplerType << std::endl;
+   os << std::setw(25) << "DistSampler Algorithm"   << " : " << std::setw(15) << fAlgoType << std::endl;
+   os << std::setw(25) << "Print Level"            << " : " << std::setw(15) << fLevel << std::endl;
+   
+   if (ExtraOptions()) { 
+      os << fSamplerType << " specific options :"  << std::endl;
+      ExtraOptions()->Print(os);
+   }
+}
+
+IOptions & DistSamplerOptions::Default(const char * name) { 
+   // create default extra options for the given algorithm type 
+   return GenAlgoOptions::Default(name);
+}
+
+IOptions * DistSamplerOptions::FindDefault(const char * name) { 
+   // find extra options for the given algorithm type 
+   return GenAlgoOptions::FindDefault(name);
+}
+
+void DistSamplerOptions::PrintDefault(const char * name, std::ostream & os) {
+   //print default options
+   DistSamplerOptions tmp;
+   tmp.Print(os);
+   if (!tmp.ExtraOptions() ) {
+      IOptions * opt = FindDefault(name);
+      os << "Specific options for "  << name << std::endl;
+      if (opt) opt->Print(os);
+   }
+}
+
+
+
+
+} // end namespace Math
+
+} // end namespace ROOT
+