#include <iomanip>
#include <cassert>
#include "Riostream.h"
#include "TMath.h"
#include "TF1.h"
#include "TH1F.h"
#include "TMVA/PDF.h"
#include "TMVA/TSpline1.h"
#include "TMVA/TSpline2.h"
#include "TMVA/Version.h"
const Int_t TMVA::PDF::fgNbin_PdfHist = 10000;
const Bool_t TMVA::PDF::fgManualIntegration = kTRUE;
const Double_t TMVA::PDF::fgEpsilon = 1.0e-12;
TMVA::PDF* TMVA::PDF::fgThisPDF = 0;
ClassImp(TMVA::PDF)
#ifdef _WIN32
#pragma warning ( disable : 4355 )
#endif
TMVA::PDF::PDF()
: fUseHistogram ( kFALSE ),
fNsmooth ( 0 ),
fInterpolMethod( PDF::kSpline0 ),
fSpline ( 0 ),
fPDFHist ( 0 ),
fHist ( 0 ),
fHistOriginal ( 0 ),
fGraph ( 0 ),
fIGetVal ( 0 ),
fKDEtype ( KDEKernel::kNone ),
fKDEiter ( KDEKernel::kNonadaptiveKDE ),
fFineFactor ( 0 ),
fReadingVersion( 0 ),
fLogger ( this )
{
fgThisPDF = this;
}
TMVA::PDF::PDF( const TH1 *hist, PDF::EInterpolateMethod method, Int_t nsmooth, Bool_t checkHist )
: fUseHistogram ( kFALSE ),
fNsmooth ( nsmooth ),
fInterpolMethod( method ),
fSpline ( 0 ),
fPDFHist ( 0 ),
fHist ( 0 ),
fHistOriginal ( 0 ),
fGraph ( 0 ),
fIGetVal ( 0 ),
fKDEtype ( KDEKernel::kNone ),
fKDEiter ( KDEKernel::kNonadaptiveKDE ),
fKDEborder ( KDEKernel::kNoTreatment ),
fFineFactor ( 0. ),
fReadingVersion( 0 ),
fLogger ( this )
{
fgThisPDF = this;
if (hist == NULL) fLogger << kFATAL << "Called without valid histogram pointer!" << Endl;
if (hist->GetEntries() <= 0)
fLogger << kFATAL << "Number of entries <= 0 in histogram: " << hist->GetTitle() << Endl;
if (nsmooth<0) fLogger << kFATAL << "PDF construction called with nsmooth<0" << Endl;
fHistOriginal = (TH1F*)hist->Clone( TString(hist->GetName()) + "_original" );
fHist = (TH1F*)hist->Clone( TString(hist->GetName()) + "_smoothed" );
fHistOriginal->SetTitle( fHistOriginal->GetName() );
fHist ->SetTitle( fHist->GetName() );
fHistOriginal->SetDirectory(0);
fHist ->SetDirectory(0);
BuildPDF( checkHist );
}
TMVA::PDF::PDF( const TH1* hist, KDEKernel::EKernelType ktype, KDEKernel::EKernelIter kiter,
KDEKernel::EKernelBorder kborder, Float_t FineFactor)
: fUseHistogram ( kFALSE ),
fNsmooth (-1 ),
fInterpolMethod( PDF::kSpline0 ),
fSpline ( 0 ),
fPDFHist ( 0 ),
fHist ( 0 ),
fHistOriginal ( 0 ),
fGraph ( 0 ),
fIGetVal ( 0 ),
fKDEtype ( ktype ),
fKDEiter ( kiter ),
fKDEborder ( kborder ),
fFineFactor ( FineFactor),
fLogger ( this )
{
if (hist == NULL) fLogger << kFATAL << "Called without valid histogram pointer!" << Endl;
if (hist->GetEntries() <= 0)
fLogger << kFATAL << "Number of entries <= 0 in histogram: " << hist->GetTitle() << Endl;
fLogger << "Create "
<< ((kiter == KDEKernel::kNonadaptiveKDE) ? "nonadaptive " :
(kiter == KDEKernel::kAdaptiveKDE) ? "adaptive " : "??? ")
<< ((ktype == KDEKernel::kGauss) ? "Gauss " : "??? ")
<< "type KDE kernel for histogram: \"" << hist->GetName() << "\""
<< Endl;
fHistOriginal = (TH1F*)hist->Clone( TString(hist->GetName()) + "_original" );
fHist = (TH1F*)hist->Clone( TString(hist->GetName()) + "_smoothed" );
fHistOriginal->SetDirectory(0);
fHist ->SetDirectory(0);
FillKDEToHist();
}
TMVA::PDF::~PDF()
{
if (fSpline != NULL) delete fSpline;
if (fHist != NULL) delete fHist;
if (fPDFHist != NULL) delete fPDFHist;
if (fHistOriginal != NULL) delete fHistOriginal;
if (fIGetVal != NULL) delete fIGetVal;
if (fGraph != NULL) delete fGraph;
}
void TMVA::PDF::BuildPDF( Bool_t checkHist )
{
if (fInterpolMethod != PDF::kSpline0 && checkHist) CheckHist();
if (fNsmooth > 0) fHist->Smooth( fNsmooth );
if(fGraph) delete fGraph;
fGraph = new TGraph( fHist );
if(fSpline) { delete fSpline; fSpline=0; }
switch (fInterpolMethod) {
case kSpline0:
fUseHistogram = kTRUE;
break;
case kSpline1:
fSpline = new TMVA::TSpline1( "spline1", fGraph );
break;
case kSpline2:
fSpline = new TMVA::TSpline2( "spline2", fGraph );
break;
case kSpline3:
fSpline = new TSpline3( "spline3", fGraph );
break;
case kSpline5:
fSpline = new TSpline5( "spline5", fGraph );
break;
default:
fLogger << kWARNING << "No valid interpolation method given! Use Spline3" << Endl;
fSpline = new TMVA::TSpline2( "spline2", fGraph );
}
FillSplineToHist();
if (!UseHistogram()) {
fSpline->SetTitle( (TString)fHist->GetTitle() + fSpline->GetTitle() );
fSpline->SetName ( (TString)fHist->GetName() + fSpline->GetName() );
}
Double_t integral = GetIntegral();
if (integral < 0) fLogger << kFATAL << "Integral: " << integral << " <= 0" << Endl;
if (integral>0) fPDFHist->Scale( 1.0/integral );
}
void TMVA::PDF::FillKDEToHist()
{
fPDFHist = new TH1F( "", "", fgNbin_PdfHist, GetXmin(), GetXmax() );
fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_KDE" );
fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_KDE" );
KDEKernel *kern = new TMVA::KDEKernel(fKDEiter,
fHist,
fPDFHist->GetBinLowEdge(1),
fPDFHist->GetBinLowEdge(fPDFHist->GetNbinsX()+1),
fKDEborder,
fFineFactor);
kern->SetKernelType(fKDEtype);
Float_t histoLowEdge=fHist->GetBinLowEdge(1);
Float_t histoUpperEdge=fHist->GetBinLowEdge(fHist->GetNbinsX()+1);
for (Int_t i=1;i<fHist->GetNbinsX();i++) {
for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
fPDFHist->AddBinContent(j,fHist->GetBinContent(i)*
kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
fPDFHist->GetBinLowEdge(j+1),
fHist->GetBinCenter(i),
i)
);
}
if (fKDEborder == 3) {
if (i < fHist->GetNbinsX()/5 ) {
for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
fPDFHist->AddBinContent(j,fHist->GetBinContent(i)*
kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
fPDFHist->GetBinLowEdge(j+1),
2*histoLowEdge-fHist->GetBinCenter(i),
i)
);
}
}
if (i > 4*fHist->GetNbinsX()/5) {
for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
fPDFHist->AddBinContent(j,fHist->GetBinContent(i)*
kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
fPDFHist->GetBinLowEdge(j+1),
2*histoUpperEdge-fHist->GetBinCenter(i),
i)
);
}
}
}
}
fPDFHist->SetEntries(fHist->GetEntries());
delete kern;
Double_t integral = GetIntegral();
if (integral < 0) fLogger << kFATAL << "Integral: " << integral << " <= 0" << Endl;
if (integral>0) fPDFHist->Scale( 1.0/integral );
}
void TMVA::PDF::FillSplineToHist()
{
if (UseHistogram()) {
fPDFHist = (TH1*)fHist->Clone();
fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_spline0" );
fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_spline0" );
}
else {
fPDFHist = new TH1F( "", "", fgNbin_PdfHist, GetXmin(), GetXmax() );
fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_" + fSpline->GetTitle() );
fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_" + fSpline->GetTitle() );
for (Int_t bin=1; bin <= fgNbin_PdfHist; bin++) {
Double_t x = fPDFHist->GetBinCenter( bin );
Double_t y = fSpline->Eval( x );
if (y <= fgEpsilon) y = fHist->GetBinContent( fHist->FindBin( x ) );
fPDFHist->SetBinContent( bin, TMath::Max(y, fgEpsilon) );
}
}
fPDFHist->SetDirectory(0);
}
void TMVA::PDF::CheckHist() const
{
if (fHist == NULL) {
fLogger << kFATAL << "<CheckHist> Called without valid histogram pointer!" << Endl;
}
Int_t nbins = fHist->GetNbinsX();
Int_t emptyBins=0;
for (Int_t bin=1; bin<=nbins; bin++)
if (fHist->GetBinContent(bin) == 0) emptyBins++;
if (((Float_t)emptyBins/(Float_t)nbins) > 0.5) {
fLogger << kWARNING << "More than 50% (" << (((Float_t)emptyBins/(Float_t)nbins)*100)
<<"%) of the bins in hist '"
<< fHist->GetName() << "' are empty!" << Endl;
fLogger << kWARNING << "X_min=" << GetXmin()
<< " mean=" << fHist->GetMean() << " X_max= " << GetXmax() << Endl;
}
}
void TMVA::PDF::ValidatePDF( TH1* originalHist ) const
{
if (!originalHist) originalHist = fHistOriginal;
Int_t nbins = originalHist->GetNbinsX();
if (originalHist->GetSumw2N() == 0) originalHist->Sumw2();
Double_t chi2 = 0;
Int_t ndof = 0;
Int_t nc1 = 0;
Int_t nc2 = 0;
Int_t nc3 = 0;
Int_t nc6 = 0;
for (Int_t bin=1; bin<=nbins; bin++) {
Double_t x = originalHist->GetBinCenter( bin );
Double_t y = originalHist->GetBinContent( bin );
Double_t ey = originalHist->GetBinError( bin );
Int_t binPdfHist = fPDFHist->FindBin( x );
Double_t yref = GetVal( x );
Double_t rref = ( originalHist->GetSumOfWeights()/fPDFHist->GetSumOfWeights() *
originalHist->GetBinWidth( bin )/fPDFHist->GetBinWidth( binPdfHist ) );
if (y > 0) {
ndof++;
Double_t d = TMath::Abs( (y - yref*rref)/ey );
chi2 += d*d;
if (d > 1) { nc1++; if (d > 2) { nc2++; if (d > 3) { nc3++; if (d > 6) nc6++; } } }
}
}
fLogger << "Validation result for PDF \"" << originalHist->GetTitle() << "\"" << ": " << Endl;
fLogger << Form( " chi2/ndof(!=0) = %.1f/%i = %.2f (Prob = %.2f)",
chi2, ndof, chi2/ndof, TMath::Prob( chi2, ndof ) ) << Endl;
fLogger << Form( " #bins-found(#expected-bins) deviating > [1,2,3,6] sigmas: " \
"[%i(%i),%i(%i),%i(%i),%i(%i)]",
nc1, Int_t(TMath::Prob(1.0,1)*ndof), nc2, Int_t(TMath::Prob(4.0,1)*ndof),
nc3, Int_t(TMath::Prob(9.0,1)*ndof), nc6, Int_t(TMath::Prob(36.0,1)*ndof) ) << Endl;
}
Double_t TMVA::PDF::GetIntegral() const
{
Double_t integral = fPDFHist->GetSumOfWeights();
if (!UseHistogram()) integral *= GetPdfHistBinWidth();
return integral;
}
Double_t TMVA::PDF::IGetVal( Double_t* x, Double_t* )
{
return ThisPDF()->GetVal( x[0] );
}
Double_t TMVA::PDF::GetIntegral( Double_t xmin, Double_t xmax )
{
Double_t integral = 0;
if (fgManualIntegration) {
Int_t imin = fPDFHist->FindBin(xmin);
Int_t imax = fPDFHist->FindBin(xmax);
if (imin < 1) imin = 1;
if (imax > fPDFHist->GetNbinsX()) imax = fPDFHist->GetNbinsX();
for (Int_t bini = imin; bini <= imax; bini++) {
Double_t x = fPDFHist->GetBinCenter(bini);
Double_t dx = fPDFHist->GetBinWidth(bini);
if (bini == imin) dx = dx/2.0 + (x - xmin);
else if (bini == imax) dx = dx/2.0 + (xmax - x);
assert( dx > 0 );
integral += fPDFHist->GetBinContent(bini)*dx;
}
}
else {
if (fIGetVal == 0)
fIGetVal = new TF1( "IGetVal", PDF::IGetVal, GetXmin(), GetXmax(), 0 );
integral = fIGetVal->Integral( xmin, xmax );
}
return integral;
}
Double_t TMVA::PDF::GetVal( Double_t x ) const
{
Int_t bin = fPDFHist->FindBin(x);
bin = TMath::Max(bin,1);
bin = TMath::Min(bin,fPDFHist->GetNbinsX());
Double_t retval = 0;
if (UseHistogram()) {
retval = fPDFHist->GetBinContent( bin );
}
else {
Int_t nextbin = bin;
if ((x > fPDFHist->GetBinCenter(bin) && bin != fPDFHist->GetNbinsX()) || bin == 1)
nextbin++;
else
nextbin--;
Double_t dx = fPDFHist->GetBinCenter( bin ) - fPDFHist->GetBinCenter( nextbin );
Double_t dy = fPDFHist->GetBinContent( bin ) - fPDFHist->GetBinContent( nextbin );
retval = fPDFHist->GetBinContent( bin ) + (x - fPDFHist->GetBinCenter( bin ))*dy/dx;
}
return TMath::Max( retval, fgEpsilon );
}
ostream& TMVA::operator<< ( ostream& os, const PDF& pdf )
{
os << "NSmooth " << pdf.fNsmooth << endl;
os << "InterpolMethod " << pdf.fInterpolMethod << endl;
os << "KDE_type " << pdf.fKDEtype << endl;
os << "KDE_iter " << pdf.fKDEiter << endl;
os << "KDE_border " << pdf.fKDEborder << endl;
os << "KDE_finefactor " << pdf.fFineFactor << endl;
TH1 * histToWrite = pdf.GetOriginalHist();
const Int_t nBins = histToWrite->GetNbinsX();
os << "Histogram "
<< histToWrite->GetName()
<< " " << nBins
<< " " << setprecision(12) << histToWrite->GetXaxis()->GetXmin()
<< " " << setprecision(12) << histToWrite->GetXaxis()->GetXmax()
<< endl;
os << "Weights " << endl;
os << std::setprecision(8);
for (Int_t i=0; i<nBins; i++) {
os << std::setw(15) << std::left << histToWrite->GetBinContent(i+1) << " ";
if ((i+1)%5==0) os << endl;
}
return os;
}
istream& TMVA::operator>> ( istream& istr, PDF& pdf )
{
TString devnullS;
Int_t valI;
Int_t nbins;
Float_t xmin, xmax;
TString hname="_original";
Bool_t doneReading = kFALSE;
while (!doneReading) {
istr >> devnullS;
if (devnullS=="NSmooth") istr >> pdf.fNsmooth;
else if (devnullS == "InterpolMethod") { istr >> valI; pdf.fInterpolMethod = PDF::EInterpolateMethod(valI);}
else if (devnullS == "KDE_type") { istr >> valI; pdf.fKDEtype = KDEKernel::EKernelType(valI); }
else if (devnullS == "KDE_iter") { istr >> valI; pdf.fKDEiter = KDEKernel::EKernelIter(valI);}
else if (devnullS == "KDE_border") { istr >> valI; pdf.fKDEborder = KDEKernel::EKernelBorder(valI);}
else if (devnullS == "KDE_finefactor") {
istr >> pdf.fFineFactor;
if (pdf.GetReadingVersion() < TMVA_VERSION(3,7,3)) {
istr >> nbins >> xmin >> xmax;
doneReading = kTRUE;
}
}
else if (devnullS == "Histogram") { istr >> hname >> nbins >> xmin >> xmax;}
else if (devnullS == "Weights") { doneReading = kTRUE; }
};
TString hnameSmooth = hname;
hnameSmooth.ReplaceAll( "_original", "_smoothed" );
if (pdf.fHistOriginal != 0) delete pdf.fHistOriginal;
if (pdf.fHist != 0) delete pdf.fHist;
pdf.fHistOriginal = new TH1F( hname, hname, nbins, xmin, xmax );
pdf.fHist = new TH1F( hnameSmooth, hnameSmooth, nbins, xmin, xmax );
pdf.fHistOriginal->SetDirectory(0);
pdf.fHist->SetDirectory(0);
Float_t val;
for (Int_t i=0; i<nbins; i++) {
istr >> val;
pdf.fHistOriginal->SetBinContent(i+1,val);
pdf.fHist->SetBinContent(i+1,val);
}
if (pdf.fNsmooth>=0) pdf.BuildPDF();
else pdf.FillKDEToHist();
return istr;
}
Last change: Wed Jun 25 08:48:44 2008
Last generated: 2008-06-25 08:48
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.