TUnuran.cxx

Go to the documentation of this file.

// @(#)root/unuran:$Id$
// 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 <cassert>
 
#include <unuran.h>
 
#include "TError.h"
 
 
TUnuran::TUnuran(TRandom * r, unsigned int debugLevel) :
   fGen(nullptr),
   fUdistr(nullptr),
   fUrng(nullptr),
   fRng(r)
{
   // constructor implementation with a ROOT random generator
   // if no generator is given the ROOT default is used
   if (fRng == nullptr) 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 != nullptr) unur_free(fGen);
   if (fUrng != nullptr) unur_urng_free(fUrng);
  // we can delete now the distribution object
   if (fUdistr != nullptr) unur_distr_free(fUdistr);
}
 
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 == nullptr) {
      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 std::unique_ptr to manage previously existing distribution objects
   TUnuranContDist * distNew = distr.Clone();
   fDist.reset(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 std::unique_ptr to manage previously existing distribution objects
   TUnuranMultiContDist * distNew = distr.Clone();
   fDist.reset(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 std::unique_ptr to manage previously existing distribution objects
   TUnuranDiscrDist * distNew = distr.Clone();
   fDist.reset(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 std::unique_ptr to manage previously existing distribution objects
   TUnuranEmpDist * distNew = distr.Clone();
   fDist.reset(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 == nullptr) return false;
   if (fGen == nullptr) return false;
 
   fUrng = unur_urng_new(&UnuranRng<TRandom>::Rndm, fRng );
   if (fUrng == nullptr) 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 continuous univariate distribution
   if (fUdistr != nullptr)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_cont_new();
   if (fUdistr == nullptr) 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 != nullptr)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_cvec_new(dist.NDim() );
   if (fUdistr == nullptr) 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 != nullptr || xmax != nullptr ) {
      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 != nullptr) {
      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 != nullptr)  unur_distr_free(fUdistr);
   if (dist.NDim() == 1)
      fUdistr = unur_distr_cemp_new();
   else
      fUdistr = unur_distr_cvemp_new(dist.NDim() );
 
   if (fUdistr == nullptr) 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 != nullptr)  unur_distr_free(fUdistr);
   fUdistr = unur_distr_discr_new();
   if (fUdistr == nullptr) return false;
   unsigned int ret = 0;
   // if a probability mesh function is provided
   if (dist.ProbVec().empty()) {
      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() {
 
   // internal function to set a method from a distribution and
   // initialize unuran with the given distribution.
   if (fUdistr == nullptr) return false;
 
   struct unur_slist *mlist = nullptr;
 
   UNUR_PAR * par = _unur_str2par(fUdistr, fMethod.c_str(), &mlist);
   if (par == nullptr) {
      Error("SetMethod","missing distribution information or syntax error");
      if (mlist != nullptr)  _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 != nullptr )  unur_free(fGen);
   fGen = unur_init(par);
   _unur_slist_free(mlist);
   if (fGen == nullptr) {
      Error("SetMethod","initializing Unuran: condition for method violated");
      return false;
   }
   return true;
 }
 
std::string TUnuran::GetInfo(bool extended)
{
   // get information string about Unuran generator
   if (!fGen) return std::string();
   return std::string(unur_gen_info(fGen, extended));
}
 
std::string TUnuran::GetGenId() const
{
   // get Unuran generator ID
   if (!fGen) return std::string();
   return std::string(unur_get_genid(fGen));
}
 
int TUnuran::GetDimension() const
{
   // get dimension of the UNURAN generator
   if (!fGen) return 0;
   return unur_get_dimension(fGen);
}
 
int TUnuran::GetDistType() const
{
   // get type of distribution
   if (!fGen) return -1;
   return unur_distr_get_type (unur_get_distr(fGen));
}
 
bool TUnuran::IsDistCont() const {
   if (!fGen) return false;
   return unur_distr_is_cont (unur_get_distr(fGen));
}
bool TUnuran::IsDistMultiCont() const {
   if (!fGen) return false;
   return unur_distr_is_cvec (unur_get_distr(fGen));
}
bool TUnuran::IsDistDiscrete() const {
   if (!fGen) return false;
   return unur_distr_is_discr (unur_get_distr(fGen));
}
bool TUnuran::IsDistEmpirical() const {
   if (!fGen) return false;
   return unur_distr_is_cemp (unur_get_distr(fGen));
}
 
int TUnuran::SampleDiscr()
{
   // sample one-dimensional distribution
   assert(fGen != nullptr);
   return unur_sample_discr(fGen);
}
 
double TUnuran::Sample()
{
   // sample one-dimensional distribution
   assert(fGen != nullptr);
   return unur_sample_cont(fGen);
}
 
bool TUnuran::SampleMulti(double * x)
{
   // sample multidimensional distribution
   if (fGen == nullptr) 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 == nullptr) 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) {
   // initialization for a Poisson
   double p[1];
   p[0] = mu;
 
   fUdistr = unur_distr_poisson(p,1);
 
   fMethod = method;
   if (fUdistr == nullptr) return false;
   if (! SetMethodAndInit() ) return false;
   if (! SetRandomGenerator() ) return false;
   return true;
}
 
bool TUnuran::InitBinomial(unsigned int ntot, double prob, const std::string & method ) {
   // initialization for a Binomial
   double par[2];
   par[0] = ntot;
   par[1] = prob;
   fUdistr = unur_distr_binomial(par,2);
 
   fMethod = method;
   if (fUdistr == nullptr) 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);
}