#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <cassert>
#include <limits>
#include "TMVA/Event.h"
#include "TMVA/Tools.h"
#include "TMVA/PDEFoam.h"
#include "TMVA/MsgLogger.h"
#include "TMVA/Types.h"
#ifndef ROOT_TStyle
#include "TStyle.h"
#endif
#ifndef ROOT_TObject
#include "TObject.h"
#endif
#ifndef ROOT_TH1D
#include "TH1D.h"
#endif
#ifndef ROOT_TMath
#include "TMath.h"
#endif
#ifndef ROOT_TVectorT
#include "TVectorT.h"
#endif
#ifndef ROOT_TRandom3
#include "TRandom3.h"
#endif
#ifndef ROOT_TColor
#include "TColor.h"
#endif
ClassImp(TMVA::PDEFoam)
static const Float_t kHigh= FLT_MAX;
static const Float_t kVlow=-FLT_MAX;
using namespace std;
TMVA::PDEFoam::PDEFoam() :
fName("PDEFoam"),
fDim(0),
fNCells(0),
fNBin(5),
fNSampl(2000),
fEvPerBin(0),
fMaskDiv(0),
fInhiDiv(0),
fNoAct(1),
fLastCe(-1),
fCells(0),
fHistEdg(0),
fRvec(0),
fPseRan(new TRandom3(4356)),
fAlpha(0),
fFoamType(kSeparate),
fXmin(0),
fXmax(0),
fNElements(0),
fNmin(100),
fMaxDepth(0),
fVolFrac(1.0/15.0),
fFillFoamWithOrigWeights(kFALSE),
fDTSeparation(kFoam),
fPeekMax(kTRUE),
fDistr(NULL),
fTimer(new Timer(0, "PDEFoam", kTRUE)),
fVariableNames(new TObjArray()),
fLogger(new MsgLogger("PDEFoam"))
{
if (fVariableNames)
fVariableNames->SetOwner(kTRUE);
}
TMVA::PDEFoam::PDEFoam(const TString& name) :
fName(name),
fDim(0),
fNCells(1000),
fNBin(5),
fNSampl(2000),
fEvPerBin(0),
fMaskDiv(0),
fInhiDiv(0),
fNoAct(1),
fLastCe(-1),
fCells(0),
fHistEdg(0),
fRvec(0),
fPseRan(new TRandom3(4356)),
fAlpha(0),
fFoamType(kSeparate),
fXmin(0),
fXmax(0),
fNElements(0),
fNmin(100),
fMaxDepth(0),
fVolFrac(1.0/15.0),
fFillFoamWithOrigWeights(kFALSE),
fDTSeparation(kFoam),
fPeekMax(kTRUE),
fDistr(NULL),
fTimer(new Timer(1, "PDEFoam", kTRUE)),
fVariableNames(new TObjArray()),
fLogger(new MsgLogger("PDEFoam"))
{
if(strlen(name) > 128)
Log() << kFATAL << "Name too long " << name.Data() << Endl;
if (fVariableNames)
fVariableNames->SetOwner(kTRUE);
}
TMVA::PDEFoam::~PDEFoam()
{
delete fVariableNames;
delete fTimer;
if (fDistr) delete fDistr;
if (fPseRan) delete fPseRan;
if (fXmin) delete [] fXmin; fXmin=0;
if (fXmax) delete [] fXmax; fXmax=0;
ResetCellElements();
if(fCells!= 0) {
for(Int_t i=0; i<fNCells; i++) delete fCells[i];
delete [] fCells;
}
delete [] fRvec;
delete [] fAlpha;
delete [] fMaskDiv;
delete [] fInhiDiv;
delete fLogger;
}
TMVA::PDEFoam::PDEFoam(const PDEFoam &from) :
TObject(from)
, fDim(0)
, fNCells(0)
, fNBin(0)
, fNSampl(0)
, fEvPerBin(0)
, fMaskDiv(0)
, fInhiDiv(0)
, fNoAct(0)
, fLastCe(0)
, fCells(0)
, fHistEdg(0)
, fRvec(0)
, fPseRan(0)
, fAlpha(0)
, fFoamType(kSeparate)
, fXmin(0)
, fXmax(0)
, fNElements(0)
, fNmin(0)
, fMaxDepth(0)
, fVolFrac(1.0/15.0)
, fFillFoamWithOrigWeights(kFALSE)
, fDTSeparation(kFoam)
, fPeekMax(kTRUE)
, fDistr(0)
, fTimer(0)
, fVariableNames(0)
, fLogger(new MsgLogger(*from.fLogger))
{
Log() << kFATAL << "COPY CONSTRUCTOR NOT IMPLEMENTED" << Endl;
if (fVariableNames)
fVariableNames->SetOwner(kTRUE);
}
void TMVA::PDEFoam::SetDim(Int_t kDim)
{
if (kDim < 1)
Log() << kFATAL << "<SetDim>: Dimension is zero or negative!" << Endl;
fDim = kDim;
if (fXmin) delete [] fXmin;
if (fXmax) delete [] fXmax;
fXmin = new Double_t[GetTotDim()];
fXmax = new Double_t[GetTotDim()];
}
void TMVA::PDEFoam::SetXmin(Int_t idim, Double_t wmin)
{
if (idim<0 || idim>=GetTotDim())
Log() << kFATAL << "<SetXmin>: Dimension out of bounds!" << Endl;
fXmin[idim]=wmin;
}
void TMVA::PDEFoam::SetXmax(Int_t idim, Double_t wmax)
{
if (idim<0 || idim>=GetTotDim())
Log() << kFATAL << "<SetXmax>: Dimension out of bounds!" << Endl;
fXmax[idim]=wmax;
}
void TMVA::PDEFoam::Create()
{
Bool_t addStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
if(fPseRan==0) Log() << kFATAL << "Random number generator not set" << Endl;
if(fDistr==0) Log() << kFATAL << "Distribution function not set" << Endl;
if(fDim==0) Log() << kFATAL << "Zero dimension not allowed" << Endl;
fRvec = new Double_t[fDim];
if(fRvec==0) Log() << kFATAL << "Cannot initialize buffer fRvec" << Endl;
if(fDim>0){
fAlpha = new Double_t[fDim];
if(fAlpha==0) Log() << kFATAL << "Cannot initialize buffer fAlpha" << Endl;
}
if(fInhiDiv == 0){
fInhiDiv = new Int_t[fDim];
for(Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
}
if(fMaskDiv == 0){
fMaskDiv = new Int_t[fDim];
for(Int_t i=0; i<fDim; i++) fMaskDiv[i]=1;
}
fHistEdg = new TObjArray(fDim);
for(Int_t i=0; i<fDim; i++){
TString hname, htitle;
hname = fName+TString("_HistEdge_");
hname += i;
htitle = TString("Edge Histogram No. ");
htitle += i;
(*fHistEdg)[i] = new TH1D(hname.Data(),htitle.Data(),fNBin,0.0, 1.0);
((TH1D*)(*fHistEdg)[i])->Sumw2();
}
ResetCellElements();
InitCells();
Grow();
TH1::AddDirectory(addStatus);
ResetCellElements();
}
void TMVA::PDEFoam::InitCells()
{
fLastCe =-1;
if(fCells!= 0) {
for(Int_t i=0; i<fNCells; i++) delete fCells[i];
delete [] fCells;
}
fCells = new(nothrow) PDEFoamCell*[fNCells];
if (!fCells) {
Log() << kFATAL << "not enough memory to create " << fNCells
<< " cells" << Endl;
}
for(Int_t i=0; i<fNCells; i++){
fCells[i]= new PDEFoamCell(fDim);
fCells[i]->SetSerial(i);
}
CellFill(1, 0);
for(Long_t iCell=0; iCell<=fLastCe; iCell++){
Explore( fCells[iCell] );
}
}
Int_t TMVA::PDEFoam::CellFill(Int_t status, PDEFoamCell *parent)
{
PDEFoamCell *cell;
if (fLastCe==fNCells){
Log() << kFATAL << "Too many cells" << Endl;
}
fLastCe++;
cell = fCells[fLastCe];
cell->Fill(status, parent, 0, 0);
cell->SetBest( -1);
cell->SetXdiv(0.5);
Double_t xInt2,xDri2;
if(parent!=0){
xInt2 = 0.5*parent->GetIntg();
xDri2 = 0.5*parent->GetDriv();
cell->SetIntg(xInt2);
cell->SetDriv(xDri2);
}else{
cell->SetIntg(0.0);
cell->SetDriv(0.0);
}
return fLastCe;
}
void TMVA::PDEFoam::Explore(PDEFoamCell *cell)
{
Double_t wt, dx, xBest=0, yBest;
Double_t intOld, driOld;
Long_t iev;
Double_t nevMC;
Int_t i, j, k;
Int_t nProj, kBest;
Double_t ceSum[5], xproj;
Double_t event_density = 0;
Double_t totevents = 0;
Double_t toteventsOld = 0;
PDEFoamVect cellSize(fDim);
PDEFoamVect cellPosi(fDim);
cell->GetHcub(cellPosi,cellSize);
PDEFoamCell *parent;
Double_t *xRand = new Double_t[fDim];
Double_t *volPart=0;
Double_t vol_scale = 1.0;
for (Int_t idim = 0; idim < fDim; ++idim)
vol_scale *= fXmax[idim] - fXmin[idim];
cell->CalcVolume();
dx = cell->GetVolume() * vol_scale;
intOld = cell->GetIntg();
driOld = cell->GetDriv();
toteventsOld = GetCellElement(cell, 0);
ceSum[0]=0;
ceSum[1]=0;
ceSum[2]=0;
ceSum[3]=kHigh;
ceSum[4]=kVlow;
for (i=0;i<fDim;i++) ((TH1D *)(*fHistEdg)[i])->Reset();
Double_t nevEff=0.;
for (iev=0;iev<fNSampl;iev++){
MakeAlpha();
if (fDim>0) for (j=0; j<fDim; j++) xRand[j]= cellPosi[j] +fAlpha[j]*(cellSize[j]);
wt = dx*Eval(xRand, event_density);
totevents += event_density;
nProj = 0;
if (fDim>0) {
for (k=0; k<fDim; k++) {
xproj =fAlpha[k];
((TH1D *)(*fHistEdg)[nProj])->Fill(xproj,wt);
nProj++;
}
}
ceSum[0] += wt;
ceSum[1] += wt*wt;
ceSum[2]++;
if (ceSum[3]>wt) ceSum[3]=wt;
if (ceSum[4]<wt) ceSum[4]=wt;
if (ceSum[1]>0) nevEff = ceSum[0]*ceSum[0]/ceSum[1];
else nevEff = 0;
if ( nevEff >= fNBin*fEvPerBin) break;
}
totevents *= dx;
if (fNSampl>0) totevents /= fNSampl;
if (cell==fCells[0] && ceSum[0]<=0.0){
if (ceSum[0]==0.0)
Log() << kFATAL << "No events were found during exploration of "
<< "root cell. Please check PDEFoam parameters nSampl "
<< "and VolFrac." << Endl;
else
Log() << kWARNING << "Negative number of events found during "
<< "exploration of root cell" << Endl;
}
for (k=0; k<fDim;k++){
fMaskDiv[k] =1;
if ( fInhiDiv[k]==1) fMaskDiv[k] =0;
}
kBest=-1;
nevMC = ceSum[2];
Double_t intTrue = ceSum[0]/(nevMC+0.000001);
Double_t intDriv=0.;
if (kBest == -1) Varedu(ceSum,kBest,xBest,yBest);
intDriv =sqrt(ceSum[1]/nevMC) -intTrue;
cell->SetBest(kBest);
cell->SetXdiv(xBest);
cell->SetIntg(intTrue);
cell->SetDriv(intDriv);
SetCellElement(cell, 0, totevents);
Double_t parIntg, parDriv;
for (parent = cell->GetPare(); parent!=0; parent = parent->GetPare()){
parIntg = parent->GetIntg();
parDriv = parent->GetDriv();
parent->SetIntg( parIntg +intTrue -intOld );
parent->SetDriv( parDriv +intDriv -driOld );
SetCellElement( parent, 0, GetCellElement(parent, 0) + totevents - toteventsOld);
}
delete [] volPart;
delete [] xRand;
}
void TMVA::PDEFoam::Varedu(Double_t ceSum[5], Int_t &kBest, Double_t &xBest, Double_t &yBest)
{
Double_t nent = ceSum[2];
Double_t sswAll = ceSum[1];
Double_t ssw = sqrt(sswAll)/sqrt(nent);
Double_t sswIn,sswOut,xLo,xUp;
kBest =-1;
xBest =0.5;
yBest =1.0;
Double_t maxGain=0.0;
for(Int_t kProj=0; kProj<fDim; kProj++) {
if( fMaskDiv[kProj]) {
Double_t sswtBest = kHigh;
Double_t gain =0.0;
Double_t xMin=0.0; Double_t xMax=0.0;
for(Int_t jLo=1; jLo<=fNBin; jLo++) {
Double_t aswIn=0; Double_t asswIn=0;
for(Int_t jUp=jLo; jUp<=fNBin;jUp++) {
aswIn += ((TH1D *)(*fHistEdg)[kProj])->GetBinContent(jUp);
asswIn += Sqr(((TH1D *)(*fHistEdg)[kProj])->GetBinError( jUp));
xLo=(jLo-1.0)/fNBin;
xUp=(jUp*1.0)/fNBin;
if ( (xUp-xLo) < std::numeric_limits<double>::epsilon()) sswIn=0.;
else sswIn = sqrt(asswIn) /sqrt(nent*(xUp-xLo)) *(xUp-xLo);
if ( (1.0-xUp+xLo) < std::numeric_limits<double>::epsilon()) sswOut=0.;
else if ( sswAll-asswIn < std::numeric_limits<double>::epsilon()) sswOut=0.;
else sswOut= sqrt(sswAll-asswIn)/sqrt(nent*(1.0-xUp+xLo)) *(1.0-xUp+xLo);
if( (sswIn+sswOut) < sswtBest) {
sswtBest = sswIn+sswOut;
gain = ssw-sswtBest;
xMin = xLo;
xMax = xUp;
}
}
}
Int_t iLo = (Int_t) (fNBin*xMin);
Int_t iUp = (Int_t) (fNBin*xMax);
if(gain>=maxGain) {
maxGain=gain;
kBest=kProj;
xBest=xMin;
yBest=xMax;
if(iLo == 0) xBest=yBest;
if(iUp == fNBin) yBest=xBest;
}
}
}
if( (kBest >= fDim) || (kBest<0) )
Log() << kFATAL << "Something wrong with kBest" << Endl;
}
void TMVA::PDEFoam::MakeAlpha()
{
fPseRan->RndmArray(fDim,fRvec);
for(Int_t k=0; k<fDim; k++) fAlpha[k] = fRvec[k];
}
Long_t TMVA::PDEFoam::PeekMax()
{
Long_t iCell = -1;
Long_t i;
Double_t drivMax, driv, xDiv;
Bool_t bCutNmin = kTRUE;
Bool_t bCutMaxDepth = kTRUE;
drivMax = 0;
for(i=0; i<=fLastCe; i++) {
if( fCells[i]->GetStat() == 1 ) {
driv = fCells[i]->GetDriv();
if (driv < std::numeric_limits<float>::epsilon())
continue;
xDiv = TMath::Abs(fCells[i]->GetXdiv());
if (xDiv <= std::numeric_limits<Double_t>::epsilon() ||
xDiv >= 1.0 - std::numeric_limits<Double_t>::epsilon())
continue;
if (GetMaxDepth() > 0)
bCutMaxDepth = fCells[i]->GetDepth() < GetMaxDepth();
if (GetNmin() > 0)
bCutNmin = GetCellElement(fCells[i], 0) > GetNmin();
if(driv > drivMax && bCutNmin && bCutMaxDepth) {
drivMax = driv;
iCell = i;
}
}
}
if (iCell == -1){
if (!bCutNmin)
Log() << kVERBOSE << "Warning: No cell with more than "
<< GetNmin() << " events found!" << Endl;
else if (!bCutMaxDepth)
Log() << kVERBOSE << "Warning: Maximum depth reached: "
<< GetMaxDepth() << Endl;
else
Log() << kWARNING << "<PDEFoam::PeekMax>: no more candidate cells (drivMax>0) found for further splitting." << Endl;
}
return(iCell);
}
Int_t TMVA::PDEFoam::Divide(PDEFoamCell *cell)
{
Int_t kBest;
if(fLastCe+1 >= fNCells) Log() << kFATAL << "Buffer limit is reached, fLastCe=fnBuf" << Endl;
cell->SetStat(0);
fNoAct++;
kBest = cell->GetBest();
if( kBest<0 || kBest>=fDim ) Log() << kFATAL << "Wrong kBest" << Endl;
Int_t d1 = CellFill(1, cell);
Int_t d2 = CellFill(1, cell);
cell->SetDau0((fCells[d1]));
cell->SetDau1((fCells[d2]));
Explore( (fCells[d1]) );
Explore( (fCells[d2]) );
return 1;
}
Double_t TMVA::PDEFoam::Eval(Double_t *xRand, Double_t &event_density)
{
std::vector<Double_t> xvec;
xvec.reserve(GetTotDim());
for (Int_t idim = 0; idim < GetTotDim(); ++idim)
xvec.push_back( VarTransformInvers(idim, xRand[idim]) );
return GetDistr()->Density(xvec, event_density);
}
void TMVA::PDEFoam::Grow()
{
fTimer->Init(fNCells);
Long_t iCell;
PDEFoamCell* newCell;
while ( (fLastCe+2) < fNCells ) {
iCell = PeekMax();
if ( (iCell<0) || (iCell>fLastCe) ) {
Log() << kVERBOSE << "Break: "<< fLastCe+1 << " cells created" << Endl;
for (Long_t jCell=fLastCe+1; jCell<fNCells; jCell++)
delete fCells[jCell];
fNCells = fLastCe+1;
break;
}
newCell = fCells[iCell];
OutputGrow();
if ( Divide( newCell )==0) break;
}
OutputGrow( kTRUE );
CheckAll(1);
Log() << kVERBOSE << GetNActiveCells() << " active cells created" << Endl;
}
void TMVA::PDEFoam::SetInhiDiv(Int_t iDim, Int_t inhiDiv)
{
if(fDim==0) Log() << kFATAL << "SetInhiDiv: fDim=0" << Endl;
if(fInhiDiv == 0) {
fInhiDiv = new Int_t[ fDim ];
for(Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
}
if( ( 0<=iDim) && (iDim<fDim)) {
fInhiDiv[iDim] = inhiDiv;
} else
Log() << kFATAL << "Wrong iDim" << Endl;
}
void TMVA::PDEFoam::CheckAll(Int_t level)
{
Int_t errors, warnings;
PDEFoamCell *cell;
Long_t iCell;
errors = 0; warnings = 0;
if (level==1) Log() << kVERBOSE << "Performing consistency checks for created foam" << Endl;
for(iCell=1; iCell<=fLastCe; iCell++) {
cell = fCells[iCell];
if( ((cell->GetDau0()==0) && (cell->GetDau1()!=0) ) ||
((cell->GetDau1()==0) && (cell->GetDau0()!=0) ) ) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d has only one daughter " << iCell << Endl;
}
if( (cell->GetDau0()==0) && (cell->GetDau1()==0) && (cell->GetStat()==0) ) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d has no daughter and is inactive " << iCell << Endl;
}
if( (cell->GetDau0()!=0) && (cell->GetDau1()!=0) && (cell->GetStat()==1) ) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d has two daughters and is active " << iCell << Endl;
}
if( (cell->GetPare())!=fCells[0] ) {
if ( (cell != cell->GetPare()->GetDau0()) && (cell != cell->GetPare()->GetDau1()) ) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d parent not pointing to this cell " << iCell << Endl;
}
}
if(cell->GetDau0()!=0) {
if(cell != (cell->GetDau0())->GetPare()) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d daughter 0 not pointing to this cell " << iCell << Endl;
}
}
if(cell->GetDau1()!=0) {
if(cell != (cell->GetDau1())->GetPare()) {
errors++;
if (level==1) Log() << kFATAL << "ERROR: Cell's no %d daughter 1 not pointing to this cell " << iCell << Endl;
}
}
if(cell->GetVolume()<1E-50) {
errors++;
if(level==1) Log() << kFATAL << "ERROR: Cell no. " << iCell << " has Volume of <1E-50" << Endl;
}
}
for(iCell=0; iCell<=fLastCe; iCell++) {
cell = fCells[iCell];
if( (cell->GetStat()==1) && (cell->GetVolume()<1E-11) ) {
errors++;
if(level==1) Log() << kFATAL << "ERROR: Cell no. " << iCell << " is active but Volume is 0 " << Endl;
}
}
if(level==1){
Log() << kVERBOSE << "Check has found " << errors << " errors and " << warnings << " warnings." << Endl;
}
if(errors>0){
Info("CheckAll","Check - found total %d errors \n",errors);
}
}
void TMVA::PDEFoam::PrintCell(Long_t iCell)
{
if (iCell < 0 || iCell > fLastCe) {
Log() << kWARNING << "<PrintCell(iCell=" << iCell
<< ")>: cell number " << iCell << " out of bounds!"
<< Endl;
return;
}
PDEFoamVect cellPosi(fDim), cellSize(fDim);
fCells[iCell]->GetHcub(cellPosi,cellSize);
Int_t kBest = fCells[iCell]->GetBest();
Double_t xBest = fCells[iCell]->GetXdiv();
Log() << "Cell[" << iCell << "]={ ";
Log() << " " << fCells[iCell] << " " << Endl;
Log() << " Xdiv[abs. coord.]="
<< VarTransformInvers(kBest,cellPosi[kBest] + xBest*cellSize[kBest])
<< Endl;
Log() << " Abs. coord. = (";
for (Int_t idim=0; idim<fDim; idim++) {
Log() << "dim[" << idim << "]={"
<< VarTransformInvers(idim,cellPosi[idim]) << ","
<< VarTransformInvers(idim,cellPosi[idim] + cellSize[idim])
<< "}";
if (idim < fDim-1)
Log() << ", ";
}
Log() << ")" << Endl;
fCells[iCell]->Print("1");
Log() << "Elements: [";
TVectorD *vec = (TVectorD*)fCells[iCell]->GetElement();
if (vec != NULL){
for (Int_t i=0; i<vec->GetNrows(); i++){
if (i>0) Log() << ", ";
Log() << GetCellElement(fCells[iCell], i);
}
} else
Log() << "not set";
Log() << "]" << Endl;
Log()<<"}"<<Endl;
}
void TMVA::PDEFoam::PrintCells(void)
{
for(Long_t iCell=0; iCell<=fLastCe; iCell++)
PrintCell(iCell);
}
void TMVA::PDEFoam::FillFoamCells(const Event* ev, Float_t wt)
{
std::vector<Float_t> values = ev->GetValues();
std::vector<Float_t> tvalues = VarTransform(values);
PDEFoamCell *cell = FindCell(tvalues);
SetCellElement(cell, 0, GetCellElement(cell, 0) + wt);
SetCellElement(cell, 1, GetCellElement(cell, 1) + wt*wt);
}
void TMVA::PDEFoam::ResetCellElements()
{
if (!fCells) return;
Log() << kVERBOSE << "Delete cell elements" << Endl;
for (Long_t iCell = 0; iCell < fNCells; ++iCell) {
TObject* elements = fCells[iCell]->GetElement();
if (elements) {
delete elements;
fCells[iCell]->SetElement(NULL);
}
}
}
Bool_t TMVA::PDEFoam::CellValueIsUndefined( PDEFoamCell* )
{
return kFALSE;
}
Float_t TMVA::PDEFoam::GetCellValue(const std::vector<Float_t> &xvec, ECellValue cv, PDEFoamKernelBase *kernel)
{
std::vector<Float_t> txvec(VarTransform(xvec));
if (kernel == NULL)
return GetCellValue(FindCell(txvec), cv);
else
return kernel->Estimate(this, txvec, cv);
}
std::vector<Float_t> TMVA::PDEFoam::GetCellValue( const std::map<Int_t,Float_t>& xvec, ECellValue cv )
{
std::map<Int_t,Float_t> txvec;
for (std::map<Int_t,Float_t>::const_iterator it=xvec.begin(); it!=xvec.end(); ++it)
txvec.insert(std::pair<Int_t, Float_t>(it->first, VarTransform(it->first, it->second)));
std::vector<PDEFoamCell*> cells = FindCells(txvec);
std::vector<Float_t> cell_values;
cell_values.reserve(cells.size());
for (std::vector<PDEFoamCell*>::const_iterator cell_it=cells.begin();
cell_it != cells.end(); ++cell_it)
cell_values.push_back(GetCellValue(*cell_it, cv));
return cell_values;
}
TMVA::PDEFoamCell* TMVA::PDEFoam::FindCell( const std::vector<Float_t> &xvec ) const
{
PDEFoamVect cellPosi0(GetTotDim()), cellSize0(GetTotDim());
PDEFoamCell *cell, *cell0;
cell=fCells[0];
Int_t idim=0;
while (cell->GetStat()!=1) {
idim=cell->GetBest();
cell0=cell->GetDau0();
cell0->GetHcub(cellPosi0,cellSize0);
if (xvec.at(idim)<=cellPosi0[idim]+cellSize0[idim])
cell=cell0;
else
cell=(cell->GetDau1());
}
return cell;
}
void TMVA::PDEFoam::FindCells(const std::map<Int_t, Float_t> &txvec, PDEFoamCell* cell, std::vector<PDEFoamCell*> &cells) const
{
PDEFoamVect cellPosi0(GetTotDim()), cellSize0(GetTotDim());
PDEFoamCell *cell0;
Int_t idim=0;
while (cell->GetStat()!=1) {
idim=cell->GetBest();
map<Int_t, Float_t>::const_iterator it = txvec.find(idim);
if (it != txvec.end()){
cell0=cell->GetDau0();
cell0->GetHcub(cellPosi0,cellSize0);
if (it->second <= cellPosi0[idim] + cellSize0[idim])
cell=cell0;
else
cell=cell->GetDau1();
} else {
FindCells(txvec, cell->GetDau0(), cells);
FindCells(txvec, cell->GetDau1(), cells);
return;
}
}
cells.push_back(cell);
}
std::vector<TMVA::PDEFoamCell*> TMVA::PDEFoam::FindCells(const std::vector<Float_t> &txvec) const
{
std::map<Int_t, Float_t> txvec_map;
for (UInt_t i=0; i<txvec.size(); ++i)
txvec_map.insert(std::pair<Int_t, Float_t>(i, txvec.at(i)));
std::vector<PDEFoamCell*> cells(0);
FindCells(txvec_map, fCells[0], cells);
return cells;
}
std::vector<TMVA::PDEFoamCell*> TMVA::PDEFoam::FindCells(const std::map<Int_t, Float_t> &txvec) const
{
std::vector<PDEFoamCell*> cells(0);
FindCells(txvec, fCells[0], cells);
return cells;
}
TH1D* TMVA::PDEFoam::Draw1Dim( ECellValue cell_value, Int_t nbin, PDEFoamKernelBase *kernel )
{
if ( GetTotDim()!=1 )
Log() << kFATAL << "<Draw1Dim>: function can only be used for 1-dimensional foams!"
<< Endl;
TString hname("h_1dim");
TH1D* h1=(TH1D*)gDirectory->Get(hname);
if (h1) delete h1;
h1= new TH1D(hname, "1-dimensional Foam", nbin, fXmin[0], fXmax[0]);
if (!h1) Log() << kFATAL << "ERROR: Can not create histo" << hname << Endl;
for (Int_t ibinx=1; ibinx<=h1->GetNbinsX(); ++ibinx) {
std::vector<Float_t> txvec;
txvec.push_back( VarTransform(0, h1->GetBinCenter(ibinx)) );
Float_t val = 0;
if (kernel != NULL) {
val = kernel->Estimate(this, txvec, cell_value);
} else {
val = GetCellValue(FindCell(txvec), cell_value);
}
h1->SetBinContent(ibinx, val + h1->GetBinContent(ibinx));
}
return h1;
}
TH2D* TMVA::PDEFoam::Project2( Int_t idim1, Int_t idim2, ECellValue cell_value, PDEFoamKernelBase *kernel, UInt_t nbin )
{
if ((idim1>=GetTotDim()) || (idim1<0) ||
(idim2>=GetTotDim()) || (idim2<0) ||
(idim1==idim2) )
Log() << kFATAL << "<Project2>: wrong dimensions given: "
<< idim1 << ", " << idim2 << Endl;
if (nbin>1000){
Log() << kWARNING << "Warning: number of bins too big: " << nbin
<< " Using 1000 bins for each dimension instead." << Endl;
nbin = 1000;
} else if (nbin<1) {
Log() << kWARNING << "Wrong bin number: " << nbin
<< "; set nbin=50" << Endl;
nbin = 50;
}
TString hname(Form("h_%d_vs_%d",idim1,idim2));
TH2D* h1=(TH2D*)gDirectory->Get(hname.Data());
if (h1) delete h1;
h1= new TH2D(hname.Data(), Form("var%d vs var%d",idim1,idim2), nbin, fXmin[idim1], fXmax[idim1], nbin, fXmin[idim2], fXmax[idim2]);
if (!h1) Log() << kFATAL << "ERROR: Can not create histo" << hname << Endl;
for (Int_t xbin = 1; xbin <= h1->GetNbinsX(); ++xbin) {
for (Int_t ybin = 1; ybin <= h1->GetNbinsY(); ++ybin) {
std::map<Int_t, Float_t> txvec;
txvec[idim1] = VarTransform(idim1, h1->GetXaxis()->GetBinCenter(xbin));
txvec[idim2] = VarTransform(idim2, h1->GetYaxis()->GetBinCenter(ybin));
std::vector<TMVA::PDEFoamCell*> cells = FindCells(txvec);
Float_t sum_cv = 0;
for (std::vector<TMVA::PDEFoamCell*>::const_iterator it = cells.begin();
it != cells.end(); ++it) {
PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
(*it)->GetHcub(cellPosi,cellSize);
std::vector<Float_t> tvec;
for (Int_t i=0; i<GetTotDim(); ++i) {
if ( i != idim1 && i != idim2 )
tvec.push_back(cellPosi[i] + 0.5*cellSize[i]);
else
tvec.push_back(txvec[i]);
}
if (kernel != NULL) {
sum_cv += kernel->Estimate(this, tvec, cell_value);
} else {
sum_cv += GetCellValue(FindCell(tvec), cell_value);
}
}
h1->SetBinContent(xbin, ybin, sum_cv + h1->GetBinContent(xbin, ybin));
}
}
return h1;
}
Float_t TMVA::PDEFoam::GetCellValue(const PDEFoamCell* cell, ECellValue cv)
{
switch (cv) {
case kValue:
return GetCellElement(cell, 0);
case kValueError:
return GetCellElement(cell, 1);
case kValueDensity: {
Double_t volume = cell->GetVolume();
if (volume > numeric_limits<double>::epsilon()) {
return GetCellValue(cell, kValue)/volume;
} else {
if (volume<=0){
cell->Print("1");
Log() << kWARNING << "<GetCellDensity(cell)>: ERROR: cell volume"
<< " negative or zero!"
<< " ==> return cell density 0!"
<< " cell volume=" << volume
<< " cell entries=" << GetCellValue(cell, kValue) << Endl;
} else {
Log() << kWARNING << "<GetCellDensity(cell)>: WARNING: cell volume"
<< " close to zero!"
<< " cell volume: " << volume << Endl;
}
}
}
return 0;
case kMeanValue:
return cell->GetIntg();
case kRms:
return cell->GetDriv();
case kRmsOvMean:
if (cell->GetIntg() != 0)
return cell->GetDriv()/cell->GetIntg();
else
return 0;
case kCellVolume:
return cell->GetVolume();
default:
Log() << kFATAL << "<GetCellValue>: unknown cell value" << Endl;
return 0;
}
return 0;
}
Double_t TMVA::PDEFoam::GetCellElement( const PDEFoamCell *cell, UInt_t i ) const
{
TVectorD *vec = (TVectorD*)cell->GetElement();
if (!vec || i >= (UInt_t) vec->GetNrows())
return 0;
return (*vec)(i);
}
void TMVA::PDEFoam::SetCellElement( PDEFoamCell *cell, UInt_t i, Double_t value )
{
TVectorD *vec = NULL;
if (cell->GetElement() == NULL) {
vec = new TVectorD(i+1);
vec->Zero();
(*vec)(i) = value;
cell->SetElement(vec);
} else {
vec = (TVectorD*)cell->GetElement();
if (!vec)
Log() << kFATAL << "<SetCellElement> ERROR: cell element is not a TVectorD*" << Endl;
if (i >= (UInt_t) vec->GetNrows())
vec->ResizeTo(0,i);
(*vec)(i) = value;
}
}
void TMVA::PDEFoam::OutputGrow( Bool_t finished )
{
if (finished) {
Log() << kINFO << "Elapsed time: " + fTimer->GetElapsedTime()
<< " " << Endl;
return;
}
Int_t modulo = 1;
if (fNCells >= 100) modulo = Int_t(fNCells/100);
if (fLastCe%modulo == 0) fTimer->DrawProgressBar( fLastCe );
}
void TMVA::PDEFoam::RootPlot2dim( const TString& filename, TString opt,
Bool_t createCanvas, Bool_t colors )
{
if (GetTotDim() != 2)
Log() << kFATAL << "RootPlot2dim() can only be used with "
<< "two-dimensional foams!" << Endl;
ECellValue cell_value = kValue;
Bool_t plotcellnumber = kFALSE;
Bool_t fillcells = kTRUE;
if (opt.Contains("cell_value")){
cell_value = kValue;
} else if (opt.Contains("rms_ov_mean")){
cell_value = kRmsOvMean;
} else if (opt.Contains("rms")){
cell_value = kRms;
} else {
fillcells = kFALSE;
}
if (opt.Contains("cellnumber"))
plotcellnumber = kTRUE;
std::ofstream outfile(filename, std::ios::out);
outfile<<"{" << std::endl;
if (!colors) {
outfile << "TColor *graycolors[100];" << std::endl;
outfile << "for (Int_t i=0.; i<100; i++)" << std::endl;
outfile << " graycolors[i]=new TColor(1000+i, 1-(Float_t)i/100.,1-(Float_t)i/100.,1-(Float_t)i/100.);"<< std::endl;
}
if (createCanvas)
outfile << "cMap = new TCanvas(\"" << fName << "\",\"Cell Map for "
<< fName << "\",600,600);" << std::endl;
outfile<<"TBox*a=new TBox();"<<std::endl;
outfile<<"a->SetFillStyle(0);"<<std::endl;
outfile<<"a->SetLineWidth(4);"<<std::endl;
outfile<<"TBox *b1=new TBox();"<<std::endl;
outfile<<"TText*t=new TText();"<<std::endl;
if (fillcells) {
outfile << (colors ? "gStyle->SetPalette(1, 0);" : "gStyle->SetPalette(0);")
<< std::endl;
outfile <<"b1->SetFillStyle(1001);"<<std::endl;
outfile<<"TBox *b2=new TBox();"<<std::endl;
outfile <<"b2->SetFillStyle(0);"<<std::endl;
}
else {
outfile <<"b1->SetFillStyle(0);"<<std::endl;
}
if (fillcells)
(colors ? gStyle->SetPalette(1, 0) : gStyle->SetPalette(0) );
Float_t zmin = 1E8;
Float_t zmax = -1E8;
if (fillcells) {
for (Long_t iCell=1; iCell<=fLastCe; iCell++) {
if ( fCells[iCell]->GetStat() == 1) {
Float_t value = GetCellValue(fCells[iCell], cell_value);
if (value<zmin)
zmin=value;
if (value>zmax)
zmax=value;
}
}
outfile << "// observed minimum and maximum of distribution: " << std::endl;
outfile << "// Float_t zmin = "<< zmin << ";" << std::endl;
outfile << "// Float_t zmax = "<< zmax << ";" << std::endl;
}
outfile << "// used minimum and maximum of distribution (taking into account log scale if applicable): " << std::endl;
outfile << "Float_t zmin = "<< zmin << ";" << std::endl;
outfile << "Float_t zmax = "<< zmax << ";" << std::endl;
Float_t x1,y1,x2,y2,x,y;
Float_t offs = 0.01;
Float_t lpag = 1-2*offs;
Int_t ncolors = colors ? gStyle->GetNumberOfColors() : 100;
Float_t scale = (ncolors-1)/(zmax - zmin);
PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
outfile << "// =========== Rectangular cells ==========="<< std::endl;
for (Long_t iCell=1; iCell<=fLastCe; iCell++) {
if ( fCells[iCell]->GetStat() == 1) {
fCells[iCell]->GetHcub(cellPosi,cellSize);
x1 = offs+lpag*(cellPosi[0]);
y1 = offs+lpag*(cellPosi[1]);
x2 = offs+lpag*(cellPosi[0]+cellSize[0]);
y2 = offs+lpag*(cellPosi[1]+cellSize[1]);
if (fillcells) {
Float_t value = GetCellValue(fCells[iCell], cell_value);
Int_t color;
if (colors)
color = gStyle->GetColorPalette(Int_t((value-zmin)*scale));
else
color = 1000+(Int_t((value-zmin)*scale));
outfile << "b1->SetFillColor(" << color << ");" << std::endl;
}
outfile<<"b1->DrawBox("<<x1<<","<<y1<<","<<x2<<","<<y2<<");"<<std::endl;
if (fillcells)
outfile<<"b2->DrawBox("<<x1<<","<<y1<<","<<x2<<","<<y2<<");"<<std::endl;
if (plotcellnumber) {
outfile<<"t->SetTextColor(4);"<<std::endl;
if(fLastCe<51)
outfile<<"t->SetTextSize(0.025);"<<std::endl;
else if(fLastCe<251)
outfile<<"t->SetTextSize(0.015);"<<std::endl;
else
outfile<<"t->SetTextSize(0.008);"<<std::endl;
x = offs+lpag*(cellPosi[0]+0.5*cellSize[0]);
y = offs+lpag*(cellPosi[1]+0.5*cellSize[1]);
outfile<<"t->DrawText("<<x<<","<<y<<","<<"\""<<iCell<<"\""<<");"<<std::endl;
}
}
}
outfile<<"// ============== End Rectangles ==========="<< std::endl;
outfile << "}" << std::endl;
outfile.flush();
outfile.close();
}
void TMVA::PDEFoam::FillBinarySearchTree( const Event* ev )
{
GetDistr()->FillBinarySearchTree(ev);
}
void TMVA::PDEFoam::DeleteBinarySearchTree()
{
if(fDistr) delete fDistr;
fDistr = NULL;
}