TUnuranMultiContDist.cxx

Go to the documentation of this file.

// @(#)root/unuran:$Id: TUnuranMultiContDist.cxx 35517 2010-09-21 09:43:15Z moneta $
// Authors: L. Moneta, J. Leydold Wed Feb 28 2007

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

// Implementation file for class TUnuranMultiContDist

#include "TUnuranMultiContDist.h"
#include "Math/WrappedMultiTF1.h"

#include "TF1.h"
#include <cassert>


TUnuranMultiContDist::TUnuranMultiContDist (const ROOT::Math::IMultiGenFunction & pdf, bool isLogPdf) : 
   fPdf(&pdf), 
   fIsLogPdf(isLogPdf), 
   fOwnFunc(false)
{
   //Constructor from generic function interfaces
} 


TUnuranMultiContDist::TUnuranMultiContDist (TF1 * func, unsigned int dim, bool isLogPdf) : 
   fPdf( new ROOT::Math::WrappedMultiTF1( *func, dim) ), 
   fIsLogPdf(isLogPdf), 
   fOwnFunc(true)
{
   //Constructor from a TF1 objects
} 



TUnuranMultiContDist::TUnuranMultiContDist(const TUnuranMultiContDist & rhs) : 
   TUnuranBaseDist(), 
   fPdf(0)
{
   // Implementation of copy ctor using assignment operator
   operator=(rhs);
}

TUnuranMultiContDist & TUnuranMultiContDist::operator = (const TUnuranMultiContDist &rhs) 
{
   // Implementation of assignment operator (copy only the function pointer not the function itself)
   if (this == &rhs) return *this;  // time saving self-test
   fXmin = rhs.fXmin;
   fXmax = rhs.fXmax;
   fMode = rhs.fMode;
   fIsLogPdf  = rhs.fIsLogPdf;
   fOwnFunc   = rhs.fOwnFunc;
   if (!fOwnFunc)  
      fPdf   = rhs.fPdf;
   else { 
       if (fPdf) delete fPdf;
       fPdf  = (rhs.fPdf)  ? rhs.fPdf->Clone()  : 0;  
   }
   return *this;
}

TUnuranMultiContDist::~TUnuranMultiContDist() { 
   // destructor implementation
   if (fOwnFunc && fPdf) delete fPdf; 
}


double TUnuranMultiContDist::Pdf ( const double * x) const {  
   // evaluate the distribution 
   assert(fPdf != 0);
   return (*fPdf)(x); 
}


void TUnuranMultiContDist::Gradient( const double * x, double * grad) const { 
   // do numerical derivation and return gradient in vector grad
   // grad must point to a vector of at least ndim size
   unsigned int ndim = NDim();
   for (unsigned int i = 0; i < ndim; ++i) 
      grad[i] = Derivative(x,i); 
      
   return;
}

double TUnuranMultiContDist::Derivative( const double * x, int coord) const { 
    // do numerical derivation of gradient using 5 point rule
   // use 5 point rule 

   //double eps = 0.001; 
   //const double kC1 = 8*std::numeric_limits<double>::epsilon();
   assert(fPdf != 0);

   double h = 0.001; 

   std::vector<double> xx(NDim() );
 
   xx[coord] = x[coord]+h;     double f1 = (*fPdf)(&xx.front());
   xx[coord] = x[coord]-h;     double f2 = (*fPdf)(&xx.front());

   xx[coord] = x[coord]+h/2;   double g1 = (*fPdf)(&xx.front());
   xx[coord] = x[coord]-h/2;   double g2 = (*fPdf)(&xx.front());

   //compute the central differences
   double h2    = 1/(2.*h);
   double d0    = f1 - f2;
   double d2    = 2*(g1 - g2);
   //double error  = kC1*h2*fx;  //compute the error
   double deriv = h2*(4*d2 - d0)/3.;  
   return deriv;
}

---