// @(#)root/tmva $Id$
// Author: Alexander Voigt
/**********************************************************************************
* Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
* Package: TMVA *
* Classes: PDEFoamDecisionTree *
* Web : http://tmva.sourceforge.net *
* *
* Description: *
* Implementation of decision tree like PDEFoam *
* *
* Authors (alphabetical): *
* Tancredi Carli - CERN, Switzerland *
* Dominik Dannheim - CERN, Switzerland *
* S. Jadach - Institute of Nuclear Physics, Cracow, Poland *
* Alexander Voigt - TU Dresden, Germany *
* Peter Speckmayer - CERN, Switzerland *
* *
* Copyright (c) 2010: *
* CERN, Switzerland *
* MPI-K Heidelberg, Germany *
* *
* Redistribution and use in source and binary forms, with or without *
* modification, are permitted according to the terms listed in LICENSE *
* (http://tmva.sourceforge.net/LICENSE) *
**********************************************************************************/
//_____________________________________________________________________
//
// PDEFoamDecisionTree
//
// This PDEFoam variant acts like a decision tree and stores in every
// cell the discriminant
//
// D = #events with given class / total number of events
//
// as well as the statistical error on the discriminant. It therefore
// acts as a discriminant estimator. The decision tree-like behaviour
// is achieved by overriding PDEFoamDiscriminant::Explore() to use a
// decision tree-like cell splitting algorithm (given a separation
// type).
//
// This PDEFoam variant should be booked together with the
// PDEFoamDecisionTreeDensity density estimator, which returns the
// events in a cell without sampling.
//
//_____________________________________________________________________
#ifndef ROOT_TMVA_PDEFoamDecisionTree
#include "TMVA/PDEFoamDecisionTree.h"
#endif
#ifndef ROOT_TMVA_PDEFoamDecisionTreeDensity
#include "TMVA/PDEFoamDecisionTreeDensity.h"
#endif
ClassImp(TMVA::PDEFoamDecisionTree)
//_____________________________________________________________________
TMVA::PDEFoamDecisionTree::PDEFoamDecisionTree()
: PDEFoamDiscriminant()
, fSepType(NULL)
{
// Default constructor for streamer, user should not use it.
}
//_____________________________________________________________________
TMVA::PDEFoamDecisionTree::PDEFoamDecisionTree(const TString& name, SeparationBase *sepType, UInt_t cls)
: PDEFoamDiscriminant(name, cls)
, fSepType(sepType)
{
// Parameters:
//
// - name - name of the foam
//
// - sepType - separation type used for the cell splitting (will be
// deleted in the destructor)
//
// - cls - class to consider as signal when calcualting the purity
}
//_____________________________________________________________________
TMVA::PDEFoamDecisionTree::PDEFoamDecisionTree(const PDEFoamDecisionTree &from)
: PDEFoamDiscriminant(from)
, fSepType(from.fSepType)
{
// Copy Constructor NOT IMPLEMENTED (NEVER USED)
Log() << kFATAL << "COPY CONSTRUCTOR NOT IMPLEMENTED" << Endl;
}
//_____________________________________________________________________
TMVA::PDEFoamDecisionTree::~PDEFoamDecisionTree()
{
// Destructor
// deletes fSepType
if (fSepType)
delete fSepType;
}
//_____________________________________________________________________
void TMVA::PDEFoamDecisionTree::Explore(PDEFoamCell *cell)
{
// Internal subprogram used by Create. It explores newly defined
// cell with according to the decision tree logic. The separation
// set via the 'sepType' option in the constructor.
//
// The optimal division point for eventual future cell division is
// determined/recorded. Note that links to parents and initial
// volume = 1/2 parent has to be already defined prior to calling
// this routine.
//
// Note, that according to the decision tree logic, a cell is only
// split, if the number of (unweighted) events in each dautghter
// cell is greater than fNmin.
if (!cell)
Log() << kFATAL << "<DTExplore> Null pointer given!" << Endl;
// create edge histograms
std::vector<TH1D*> hsig, hbkg, hsig_unw, hbkg_unw;
hsig.reserve(fDim);
hbkg.reserve(fDim);
hsig_unw.reserve(fDim);
hbkg_unw.reserve(fDim);
for (Int_t idim = 0; idim < fDim; idim++) {
hsig.push_back(new TH1D(Form("hsig_%i", idim),
Form("signal[%i]", idim), fNBin, fXmin[idim], fXmax[idim]));
hbkg.push_back(new TH1D(Form("hbkg_%i", idim),
Form("background[%i]", idim), fNBin, fXmin[idim], fXmax[idim]));
hsig_unw.push_back(new TH1D(Form("hsig_unw_%i", idim),
Form("signal_unw[%i]", idim), fNBin, fXmin[idim], fXmax[idim]));
hbkg_unw.push_back(new TH1D(Form("hbkg_unw_%i", idim),
Form("background_unw[%i]", idim), fNBin, fXmin[idim], fXmax[idim]));
}
// get cell position and size
PDEFoamVect cellSize(GetTotDim()), cellPosi(GetTotDim());
cell->GetHcub(cellPosi, cellSize);
// determine lower and upper cell bound
std::vector<Double_t> lb(GetTotDim()); // lower bound
std::vector<Double_t> ub(GetTotDim()); // upper bound
for (Int_t idim = 0; idim < GetTotDim(); idim++) {
lb[idim] = VarTransformInvers(idim, cellPosi[idim] - std::numeric_limits<float>::epsilon());
ub[idim] = VarTransformInvers(idim, cellPosi[idim] + cellSize[idim] + std::numeric_limits<float>::epsilon());
}
// fDistr must be of type PDEFoamDecisionTreeDensity*
PDEFoamDecisionTreeDensity *distr = dynamic_cast<PDEFoamDecisionTreeDensity*>(fDistr);
if (distr == NULL)
Log() << kFATAL << "<PDEFoamDecisionTree::Explore>: cast failed: "
<< "PDEFoamDensityBase* --> PDEFoamDecisionTreeDensity*" << Endl;
// create TMVA::Volume object needed for searching within the BST
TMVA::Volume volume(&lb, &ub);
// fill the signal and background histograms for the given volume
distr->FillHistograms(volume, hsig, hbkg, hsig_unw, hbkg_unw);
// ------ determine the best division edge
Double_t xBest = 0.5; // best division point
Int_t kBest = -1; // best split dimension
Double_t maxGain = -1.0; // maximum gain
Double_t nTotS = hsig.at(0)->Integral(0, hsig.at(0)->GetNbinsX() + 1);
Double_t nTotB = hbkg.at(0)->Integral(0, hbkg.at(0)->GetNbinsX() + 1);
Double_t nTotS_unw = hsig_unw.at(0)->Integral(0, hsig_unw.at(0)->GetNbinsX() + 1);
Double_t nTotB_unw = hbkg_unw.at(0)->Integral(0, hbkg_unw.at(0)->GetNbinsX() + 1);
for (Int_t idim = 0; idim < fDim; ++idim) {
Double_t nSelS = hsig.at(idim)->GetBinContent(0);
Double_t nSelB = hbkg.at(idim)->GetBinContent(0);
Double_t nSelS_unw = hsig_unw.at(idim)->GetBinContent(0);
Double_t nSelB_unw = hbkg_unw.at(idim)->GetBinContent(0);
for (Int_t jLo = 1; jLo < fNBin; jLo++) {
nSelS += hsig.at(idim)->GetBinContent(jLo);
nSelB += hbkg.at(idim)->GetBinContent(jLo);
nSelS_unw += hsig_unw.at(idim)->GetBinContent(jLo);
nSelB_unw += hbkg_unw.at(idim)->GetBinContent(jLo);
// proceed if total number of events in left and right cell
// is greater than fNmin
if (!((nSelS_unw + nSelB_unw) >= GetNmin() &&
(nTotS_unw - nSelS_unw + nTotB_unw - nSelB_unw) >= GetNmin()))
continue;
Double_t xLo = 1.0 * jLo / fNBin;
// calculate separation gain
Double_t gain = fSepType->GetSeparationGain(nSelS, nSelB, nTotS, nTotB);
if (gain >= maxGain) {
maxGain = gain;
xBest = xLo;
kBest = idim;
}
} // jLo
} // idim
if (kBest >= fDim || kBest < 0) {
// No best division edge found! One must ensure, that this cell
// is not chosen for splitting in PeekMax(). But since in
// PeekMax() it is ensured that cell->GetDriv() > epsilon, one
// should set maxGain to -1.0 (or even 0.0?) here.
maxGain = -1.0;
}
// set cell properties
cell->SetBest(kBest);
cell->SetXdiv(xBest);
if (nTotB + nTotS > 0)
cell->SetIntg(nTotS / (nTotB + nTotS));
else
cell->SetIntg(0.0);
cell->SetDriv(maxGain);
cell->CalcVolume();
// set cell element 0 (total number of events in cell) during
// build-up
if (GetNmin() > 0)
SetCellElement(cell, 0, nTotS + nTotB);
// clean up
for (UInt_t ih = 0; ih < hsig.size(); ih++) delete hsig.at(ih);
for (UInt_t ih = 0; ih < hbkg.size(); ih++) delete hbkg.at(ih);
for (UInt_t ih = 0; ih < hsig_unw.size(); ih++) delete hsig_unw.at(ih);
for (UInt_t ih = 0; ih < hbkg_unw.size(); ih++) delete hbkg_unw.at(ih);
}