TUnuran.cxx

Go to the documentation of this file.

// @(#)root/unuran:$Id: TUnuran.cxx 37419 2010-12-08 21:19:45Z moneta $
// Authors: L. Moneta, J. Leydold Tue Sep 26 16:25:09 2006

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// Implementation file for class TUnuran

#include "TUnuran.h"

#include "TUnuranContDist.h"
#include "TUnuranMultiContDist.h"
#include "TUnuranDiscrDist.h"
#include "TUnuranEmpDist.h"

#include "UnuranRng.h"
#include "UnuranDistrAdapter.h"

#include "TRandom.h"
#include "TSystem.h"

#include "TH1.h"

#include <cassert>


#include <unuran.h>

#include "TError.h"


TUnuran::TUnuran(TRandom * r, unsigned int debugLevel) : 
   fGen(0),
   fUdistr(0),
   fUrng(0),
   fRng(r) 
{
   // constructor implementation with a ROOT random generator
   // if no generator is given the ROOT default is used
   if (fRng == 0) fRng = gRandom; 
   // set debug level at global level 
   // (should be in a static  initialization function of the library ? )
   if ( debugLevel > 1) 
      unur_set_default_debug(UNUR_DEBUG_ALL);
   else if (debugLevel == 1) 
      unur_set_default_debug(UNUR_DEBUG_INIT);
   else
      unur_set_default_debug(UNUR_DEBUG_OFF);
      
}


TUnuran::~TUnuran() 
{
   // Destructor implementation
   if (fGen != 0) unur_free(fGen); 
   if (fUrng != 0) unur_urng_free(fUrng);
  // we can delete now the distribution object
   if (fUdistr != 0) unur_distr_free(fUdistr);
}

//private (no impl.)
TUnuran::TUnuran(const TUnuran &) 
{
   // Implementation of copy constructor.
}

TUnuran & TUnuran::operator = (const TUnuran &rhs) 
{
   // Implementation of assignment operator.
   if (this == &rhs) return *this;  // time saving self-test
   return *this;
}

bool  TUnuran::Init(const std::string & dist, const std::string & method) 
{
   // initialize with a string
   std::string s = dist + " & " + method; 
   fGen = unur_str2gen(s.c_str() ); 
   if (fGen == 0) { 
      Error("Init","Cannot create generator object"); 
      return false; 
   } 
   if (! SetRandomGenerator() ) return false; 

   return true; 
}

bool TUnuran::Init(const TUnuranContDist & distr, const std::string  & method)  
{ 
   // initialization with a distribution and and generator
   // the distribution object is copied in and managed by this class
   // use auto_ptr to manage previously existing distribution objects
   TUnuranContDist * distNew = distr.Clone(); 
   fDist = std::auto_ptr< TUnuranBaseDist>(distNew);

   fMethod = method; 
   if (! SetContDistribution(*distNew) ) return false;
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true;
}

    
bool TUnuran::Init(const TUnuranMultiContDist & distr, const std::string  & method)  
{ 
   //  initialization with a distribution and method
   // the distribution object is copied in and managed by this class
   // use auto_ptr to manage previously existing distribution objects
   TUnuranMultiContDist * distNew = distr.Clone(); 
   fDist = std::auto_ptr< TUnuranBaseDist>(distNew);

   fMethod = method;
   if (! SetMultiDistribution(*distNew) ) return false;
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true; 
}


bool TUnuran::Init(const TUnuranDiscrDist & distr, const std::string & method ) {
   //   initialization with a distribution and and generator
   // the distribution object is copied in and managed by this class
   // use auto_ptr to manage previously existing distribution objects
   TUnuranDiscrDist * distNew = distr.Clone(); 
   fDist = std::auto_ptr< TUnuranBaseDist>(distNew);

   fMethod = method; 
   if (! SetDiscreteDistribution(*distNew) ) return false;
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true;
}

bool TUnuran::Init(const TUnuranEmpDist & distr, const std::string & method ) {
   //   initialization with a distribution and and generator
   // the distribution object is copied in and managed by this class
   // use auto_ptr to manage previously existing distribution objects
   TUnuranEmpDist * distNew = distr.Clone(); 
   fDist = std::auto_ptr< TUnuranBaseDist>(distNew);
   
   fMethod = method; 
   if (distr.IsBinned()) fMethod = "hist";
   else if (distr.NDim() > 1) fMethod = "vempk";
   if (! SetEmpiricalDistribution(*distNew) ) return false;
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true;
}


bool  TUnuran::SetRandomGenerator()
{
   // set an external random generator
   if (fRng == 0) return false; 
   if (fGen == 0) return false; 

   fUrng = unur_urng_new(&UnuranRng<TRandom>::Rndm, fRng );
   if (fUrng == 0) return false; 
   unsigned int ret = 0; 
   ret |= unur_urng_set_delete(fUrng, &UnuranRng<TRandom>::Delete); 
   ret |= unur_urng_set_seed(fUrng, &UnuranRng<TRandom>::Seed);
   if (ret != 0) return false; 

   unur_chg_urng( fGen, fUrng); 
   return true;    
}

bool  TUnuran::SetContDistribution(const TUnuranContDist & dist )
{
   // internal method to set in unuran the function pointer for a continous univariate distribution 
   if (fUdistr != 0)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_cont_new(); 
   if (fUdistr == 0) return false; 
   unsigned int ret = 0; 
   ret = unur_distr_set_extobj(fUdistr, &dist);  
   if ( ! dist.IsLogPdf() ) { 
      ret |= unur_distr_cont_set_pdf(fUdistr, &ContDist::Pdf);  
      ret |= unur_distr_cont_set_dpdf(fUdistr, &ContDist::Dpdf);  
      if (dist.HasCdf() ) ret |= unur_distr_cont_set_cdf(fUdistr, &ContDist::Cdf);  
   }
   else { 
      // case user provides log of pdf 
      ret |= unur_distr_cont_set_logpdf(fUdistr, &ContDist::Pdf);  
      ret |= unur_distr_cont_set_dlogpdf(fUdistr, &ContDist::Dpdf);  
   }

   double xmin, xmax = 0; 
   if (dist.GetDomain(xmin,xmax) ) { 
      ret = unur_distr_cont_set_domain(fUdistr,xmin,xmax);  
      if (ret != 0)  { 
         Error("SetContDistribution","invalid domain xmin = %g xmax = %g ",xmin,xmax);
         return false; 
      }
   }
   if (dist.HasMode() ) { 
      ret = unur_distr_cont_set_mode(fUdistr, dist.Mode());  
      if (ret != 0)  { 
         Error("SetContDistribution","invalid mode given,  mode = %g ",dist.Mode());
         return false; 
      }
   }
   if (dist.HasPdfArea() ) { 
      ret = unur_distr_cont_set_pdfarea(fUdistr, dist.PdfArea());  
      if (ret != 0)  { 
         Error("SetContDistribution","invalid area given,  area = %g ",dist.PdfArea());
         return false; 
      }
   }

   return (ret ==0) ? true : false; 
}


bool  TUnuran::SetMultiDistribution(const TUnuranMultiContDist & dist )
{
   // internal method to set in unuran the function pointer for a multivariate distribution 
   if (fUdistr != 0)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_cvec_new(dist.NDim() ); 
   if (fUdistr == 0) return false; 
   unsigned int ret = 0; 
   ret |= unur_distr_set_extobj(fUdistr, &dist );  
   if ( ! dist.IsLogPdf() ) { 
      ret |= unur_distr_cvec_set_pdf(fUdistr, &MultiDist::Pdf);  
      ret |= unur_distr_cvec_set_dpdf(fUdistr, &MultiDist::Dpdf);  
      ret |= unur_distr_cvec_set_pdpdf(fUdistr, &MultiDist::Pdpdf);  
   }
   else { 
      ret |= unur_distr_cvec_set_logpdf(fUdistr, &MultiDist::Pdf);  
      ret |= unur_distr_cvec_set_dlogpdf(fUdistr, &MultiDist::Dpdf);  
      ret |= unur_distr_cvec_set_pdlogpdf(fUdistr, &MultiDist::Pdpdf);  
   }

   const double * xmin = dist.GetLowerDomain();
   const double * xmax = dist.GetUpperDomain();
   if ( xmin != 0 || xmax != 0 ) {
      ret = unur_distr_cvec_set_domain_rect(fUdistr,xmin,xmax);  
      if (ret != 0)  { 
         Error("SetMultiDistribution","invalid domain");
         return false; 
      }
#ifdef OLDVERS
      Error("SetMultiDistribution","domain setting not available in UNURAN 0.8.1");
#endif

   }

   const double * xmode = dist.GetMode(); 
   if (xmode != 0) { 
      ret = unur_distr_cvec_set_mode(fUdistr, xmode);  
      if (ret != 0)  { 
         Error("SetMultiDistribution","invalid mode");
         return false; 
      }
   }
   return (ret ==0) ? true : false; 
}

bool TUnuran::SetEmpiricalDistribution(const TUnuranEmpDist & dist) { 

   // internal method to set in unuran the function pointer for am empiral distribution (from histogram)
   if (fUdistr != 0)  unur_distr_free(fUdistr);
   if (dist.NDim() == 1) 
      fUdistr = unur_distr_cemp_new(); 
   else 
      fUdistr = unur_distr_cvemp_new(dist.NDim() ); 

   if (fUdistr == 0) return false; 
   unsigned int ret = 0; 


   // get info from histogram 
   if (dist.IsBinned() ) { 
      int nbins = dist.Data().size(); 
      double min = dist.LowerBin();
      double max = dist.UpperBin();
      const double * pv = &(dist.Data().front());
      ret |= unur_distr_cemp_set_hist(fUdistr, pv, nbins, min, max); 
#ifdef OLDVERS
      Error("SetEmpiricalDistribution","hist method not available in UNURAN 0.8.1");
#endif
   } 
   else { 
      const double * pv = &dist.Data().front();
      // n is number of points (size/ndim)
      int n = dist.Data().size()/dist.NDim();
      if (dist.NDim() == 1)   
         ret |= unur_distr_cemp_set_data(fUdistr, pv, n); 
      else  
         ret |= unur_distr_cvemp_set_data(fUdistr, pv, n); 
   }
   if (ret != 0) { 
      Error("SetEmpiricalDistribution","invalid distribution object");
      return false; 
   }
   return true; 
}


bool  TUnuran::SetDiscreteDistribution(const TUnuranDiscrDist & dist)
{
   // internal method to set in unuran the function pointer for a discrete univariate distribution 
   if (fUdistr != 0)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_discr_new(); 
   if (fUdistr == 0) return false; 
   unsigned int ret = 0; 
   // if a probability mesh function is provided 
   if (dist.ProbVec().size() == 0) { 
      ret = unur_distr_set_extobj(fUdistr, &dist );  
      ret |= unur_distr_discr_set_pmf(fUdistr, &DiscrDist::Pmf);  
      if (dist.HasCdf() ) ret |= unur_distr_discr_set_cdf(fUdistr, &DiscrDist::Cdf);  
 
   }
   else { 
      // case user provides vector of probabilities 
      ret |= unur_distr_discr_set_pv(fUdistr, &dist.ProbVec().front(), dist.ProbVec().size() );
   }

   int xmin, xmax = 0; 
   if (dist.GetDomain(xmin,xmax) ) { 
      ret = unur_distr_discr_set_domain(fUdistr,xmin,xmax);  
      if (ret != 0)  { 
         Error("SetDiscrDistribution","invalid domain xmin = %d xmax = %d ",xmin,xmax);
         return false; 
      }
   }
   if (dist.HasMode() ) { 
      ret = unur_distr_discr_set_mode(fUdistr, dist.Mode());  
      if (ret != 0)  { 
         Error("SetContDistribution","invalid mode given,  mode = %d ",dist.Mode());
         return false; 
      }
   }
   if (dist.HasProbSum() ) {
      ret = unur_distr_discr_set_pmfsum(fUdistr, dist.ProbSum());  
      if (ret != 0)  { 
         Error("SetContDistribution","invalid sum given,  mode = %g ",dist.ProbSum());
         return false; 
      }
   }

   return (ret ==0) ? true : false; 
}


//bool TUnuran::SetMethodAndInit(const std::string & s) { 
bool TUnuran::SetMethodAndInit() { 

   // internal function to set a method from a distribution and 
   // initialize unuran with the given distribution.
   if (fUdistr == 0) return false; 

   struct unur_slist *mlist = NULL;
   
   UNUR_PAR * par = _unur_str2par(fUdistr, fMethod.c_str(), &mlist);
   if (par == 0) {
      Error("SetMethod","missing distribution information or syntax error");
      if (mlist != 0)  _unur_slist_free(mlist);
      return false;
   }


   // set unuran to not use a private copy of the distribution object
   unur_set_use_distr_privatecopy (par, false);

   // need to free fGen if already existing ? 
   if (fGen != 0 )  unur_free(fGen);
   fGen = unur_init(par); 
   _unur_slist_free(mlist);
   if (fGen == 0) { 
      Error("SetMethod","initializing Unuran: condition for method violated");
      return false; 
   }
   return true;
 }


int TUnuran::SampleDiscr()
{
   // sample one-dimensional distribution
   assert(fGen != 0); 
   return unur_sample_discr(fGen);
}

double TUnuran::Sample()
{
   // sample one-dimensional distribution
   assert(fGen != 0); 
   return unur_sample_cont(fGen);
}

bool TUnuran::SampleMulti(double * x)
{
   // sample multidimensional distribution
   if (fGen == 0) return false;  
   unur_sample_vec(fGen,x);
   return true; 
}

void TUnuran::SetSeed(unsigned int seed) { 
   return fRng->SetSeed(seed); 
}

bool  TUnuran::SetLogLevel(unsigned int debugLevel) 
{
   if (fGen == 0) return false; 
   int ret = 0; 
   if ( debugLevel > 1) 
      ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
   else if (debugLevel == 1) 
      ret |= unur_chg_debug(fGen, UNUR_DEBUG_ALL);
   else
      ret |= unur_chg_debug(fGen, UNUR_DEBUG_OFF);

   return (ret ==0) ? true : false; 

}

bool TUnuran::InitPoisson(double mu, const std::string & method) { 
   // initializaton for a Poisson
   double p[1];
   p[0] = mu; 

   fUdistr = unur_distr_poisson(p,1);

   fMethod = method;
   if (fUdistr == 0) return false; 
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true; 
}


bool TUnuran::InitBinomial(unsigned int ntot, double prob, const std::string & method ) { 
   // initializaton for a Binomial
   double par[2];
   par[0] = ntot; 
   par[1] = prob; 
   fUdistr = unur_distr_binomial(par,2);

   fMethod = method;
   if (fUdistr == 0) return false; 
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false; 
   return true; 
}


bool TUnuran::ReInitDiscrDist(unsigned int npar, double * par) { 
   // re-initialization of UNURAN without freeing and creating a new fGen object
   // works only for pre-defined distribution by changing their parameters 
   if (!fGen ) return false;
   if (!fUdistr) return false;
   unur_distr_discr_set_pmfparams(fUdistr,par,npar);
   int iret = unur_reinit(fGen);
   if (iret) Warning("ReInitDiscrDist","re-init failed - a full initizialization must be performed");
   return (!iret); 
}

---