Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
MethodBDT.cxx
Go to the documentation of this file.
1// Author: Andreas Hoecker, Joerg Stelzer, Helge Voss, Kai Voss, Eckhard v. Toerne, Jan Therhaag
2
3/**********************************************************************************
4 * Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
5 * Package: TMVA *
6 * Class : MethodBDT (BDT = Boosted Decision Trees) *
7 * *
8 * *
9 * Description: *
10 * Analysis of Boosted Decision Trees *
11 * *
12 * Authors (alphabetical): *
13 * Andreas Hoecker <Andreas.Hocker@cern.ch> - CERN, Switzerland *
14 * Helge Voss <Helge.Voss@cern.ch> - MPI-K Heidelberg, Germany *
15 * Kai Voss <Kai.Voss@cern.ch> - U. of Victoria, Canada *
16 * Doug Schouten <dschoute@sfu.ca> - Simon Fraser U., Canada *
17 * Jan Therhaag <jan.therhaag@cern.ch> - U. of Bonn, Germany *
18 * Eckhard v. Toerne <evt@uni-bonn.de> - U of Bonn, Germany *
19 * *
20 * Copyright (c) 2005-2011: *
21 * CERN, Switzerland *
22 * U. of Victoria, Canada *
23 * MPI-K Heidelberg, Germany *
24 * U. of Bonn, Germany *
25 * *
26 * Redistribution and use in source and binary forms, with or without *
27 * modification, are permitted according to the terms listed in LICENSE *
28 * (see tmva/doc/LICENSE) *
29 **********************************************************************************/
30
31/*! \class TMVA::MethodBDT
32\ingroup TMVA
33
34Analysis of Boosted Decision Trees
35
36Boosted decision trees have been successfully used in High Energy
37Physics analysis for example by the MiniBooNE experiment
38(Yang-Roe-Zhu, physics/0508045). In Boosted Decision Trees, the
39selection is done on a majority vote on the result of several decision
40trees, which are all derived from the same training sample by
41supplying different event weights during the training.
42
43### Decision trees:
44
45Successive decision nodes are used to categorize the
46events out of the sample as either signal or background. Each node
47uses only a single discriminating variable to decide if the event is
48signal-like ("goes right") or background-like ("goes left"). This
49forms a tree like structure with "baskets" at the end (leave nodes),
50and an event is classified as either signal or background according to
51whether the basket where it ends up has been classified signal or
52background during the training. Training of a decision tree is the
53process to define the "cut criteria" for each node. The training
54starts with the root node. Here one takes the full training event
55sample and selects the variable and corresponding cut value that gives
56the best separation between signal and background at this stage. Using
57this cut criterion, the sample is then divided into two subsamples, a
58signal-like (right) and a background-like (left) sample. Two new nodes
59are then created for each of the two sub-samples and they are
60constructed using the same mechanism as described for the root
61node. The devision is stopped once a certain node has reached either a
62minimum number of events, or a minimum or maximum signal purity. These
63leave nodes are then called "signal" or "background" if they contain
64more signal respective background events from the training sample.
65
66### Boosting:
67
68The idea behind adaptive boosting (AdaBoost) is, that signal events
69from the training sample, that end up in a background node
70(and vice versa) are given a larger weight than events that are in
71the correct leave node. This results in a re-weighed training event
72sample, with which then a new decision tree can be developed.
73The boosting can be applied several times (typically 100-500 times)
74and one ends up with a set of decision trees (a forest).
75Gradient boosting works more like a function expansion approach, where
76each tree corresponds to a summand. The parameters for each summand (tree)
77are determined by the minimization of a error function (binomial log-
78likelihood for classification and Huber loss for regression).
79A greedy algorithm is used, which means, that only one tree is modified
80at a time, while the other trees stay fixed.
81
82### Bagging:
83
84In this particular variant of the Boosted Decision Trees the boosting
85is not done on the basis of previous training results, but by a simple
86stochastic re-sampling of the initial training event sample.
87
88### Random Trees:
89
90Similar to the "Random Forests" from Leo Breiman and Adele Cutler, it
91uses the bagging algorithm together and bases the determination of the
92best node-split during the training on a random subset of variables only
93which is individually chosen for each split.
94
95### Analysis:
96
97Applying an individual decision tree to a test event results in a
98classification of the event as either signal or background. For the
99boosted decision tree selection, an event is successively subjected to
100the whole set of decision trees and depending on how often it is
101classified as signal, a "likelihood" estimator is constructed for the
102event being signal or background. The value of this estimator is the
103one which is then used to select the events from an event sample, and
104the cut value on this estimator defines the efficiency and purity of
105the selection.
106
107*/
108
109
110#include "TMVA/MethodBDT.h"
111#include "TMVA/Config.h"
112
113#include "TMVA/BDTEventWrapper.h"
116#include "TMVA/Configurable.h"
117#include "TMVA/CrossEntropy.h"
118#include "TMVA/DecisionTree.h"
119#include "TMVA/DataSet.h"
120#include "TMVA/GiniIndex.h"
122#include "TMVA/Interval.h"
123#include "TMVA/IMethod.h"
124#include "TMVA/LogInterval.h"
125#include "TMVA/MethodBase.h"
127#include "TMVA/MsgLogger.h"
129#include "TMVA/PDF.h"
130#include "TMVA/Ranking.h"
131#include "TMVA/Results.h"
133#include "TMVA/SdivSqrtSplusB.h"
134#include "TMVA/SeparationBase.h"
135#include "TMVA/Timer.h"
136#include "TMVA/Tools.h"
137#include "TMVA/Types.h"
138
139#include "TRandom3.h"
140#include "TMath.h"
141#include "TMatrixTSym.h"
142#include "TGraph.h"
143
144#include <iostream>
145#include <iomanip>
146#include <algorithm>
147#include <cmath>
148#include <numeric>
149#include <unordered_map>
150
151using std::vector;
152using std::make_pair;
153
155
156
158
159////////////////////////////////////////////////////////////////////////////////
160/// The standard constructor for the "boosted decision trees".
161
163 const TString& methodTitle,
165 const TString& theOption ) :
166 TMVA::MethodBase( jobName, Types::kBDT, methodTitle, theData, theOption)
167 , fTrainSample(0)
168 , fNTrees(0)
169 , fSigToBkgFraction(0)
170 , fAdaBoostBeta(0)
171// , fTransitionPoint(0)
172 , fShrinkage(0)
173 , fBaggedBoost(kFALSE)
174 , fBaggedGradBoost(kFALSE)
175// , fSumOfWeights(0)
176 , fMinNodeEvents(0)
177 , fMinNodeSize(5)
178 , fMinNodeSizeS("5%")
179 , fNCuts(0)
180 , fUseFisherCuts(0) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
181 , fMinLinCorrForFisher(.8) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
182 , fUseExclusiveVars(0) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
183 , fUseYesNoLeaf(kFALSE)
184 , fNodePurityLimit(0)
185 , fNNodesMax(0)
186 , fMaxDepth(0)
187 , fPruneMethod(DecisionTree::kNoPruning)
188 , fPruneStrength(0)
189 , fFValidationEvents(0)
190 , fAutomatic(kFALSE)
191 , fRandomisedTrees(kFALSE)
192 , fUseNvars(0)
193 , fUsePoissonNvars(0) // don't use this initialisation, only here to make Coverity happy. Is set in Init()
194 , fUseNTrainEvents(0)
195 , fBaggedSampleFraction(0)
196 , fNoNegWeightsInTraining(kFALSE)
197 , fInverseBoostNegWeights(kFALSE)
198 , fPairNegWeightsGlobal(kFALSE)
199 , fTrainWithNegWeights(kFALSE)
200 , fDoBoostMonitor(kFALSE)
201 , fITree(0)
202 , fBoostWeight(0)
203 , fErrorFraction(0)
204 , fCss(0)
205 , fCts_sb(0)
206 , fCtb_ss(0)
207 , fCbb(0)
208 , fDoPreselection(kFALSE)
209 , fSkipNormalization(kFALSE)
210 , fHistoricBool(kFALSE)
211{
213 fSepType = NULL;
215}
216
217////////////////////////////////////////////////////////////////////////////////
218
220 const TString& theWeightFile)
222 , fTrainSample(0)
223 , fNTrees(0)
224 , fSigToBkgFraction(0)
225 , fAdaBoostBeta(0)
226// , fTransitionPoint(0)
227 , fShrinkage(0)
228 , fBaggedBoost(kFALSE)
229 , fBaggedGradBoost(kFALSE)
230// , fSumOfWeights(0)
231 , fMinNodeEvents(0)
232 , fMinNodeSize(5)
233 , fMinNodeSizeS("5%")
234 , fNCuts(0)
235 , fUseFisherCuts(0) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
236 , fMinLinCorrForFisher(.8) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
237 , fUseExclusiveVars(0) // don't use this initialisation, only here to make Coverity happy. Is set in DeclarOptions()
238 , fUseYesNoLeaf(kFALSE)
239 , fNodePurityLimit(0)
240 , fNNodesMax(0)
241 , fMaxDepth(0)
242 , fPruneMethod(DecisionTree::kNoPruning)
243 , fPruneStrength(0)
244 , fFValidationEvents(0)
245 , fAutomatic(kFALSE)
246 , fRandomisedTrees(kFALSE)
247 , fUseNvars(0)
248 , fUsePoissonNvars(0) // don't use this initialisation, only here to make Coverity happy. Is set in Init()
249 , fUseNTrainEvents(0)
250 , fBaggedSampleFraction(0)
251 , fNoNegWeightsInTraining(kFALSE)
252 , fInverseBoostNegWeights(kFALSE)
253 , fPairNegWeightsGlobal(kFALSE)
254 , fTrainWithNegWeights(kFALSE)
255 , fDoBoostMonitor(kFALSE)
256 , fITree(0)
257 , fBoostWeight(0)
258 , fErrorFraction(0)
259 , fCss(0)
260 , fCts_sb(0)
261 , fCtb_ss(0)
262 , fCbb(0)
263 , fDoPreselection(kFALSE)
264 , fSkipNormalization(kFALSE)
265 , fHistoricBool(kFALSE)
266{
268 fSepType = NULL;
270 // constructor for calculating BDT-MVA using previously generated decision trees
271 // the result of the previous training (the decision trees) are read in via the
272 // weight file. Make sure the variables correspond to the ones used in
273 // creating the "weight"-file
274}
275
276////////////////////////////////////////////////////////////////////////////////
277/// BDT can handle classification with multiple classes and regression with one regression-target.
278
286
287////////////////////////////////////////////////////////////////////////////////
288/// Define the options (their key words). That can be set in the option string.
289///
290/// know options:
291///
292/// - nTrees number of trees in the forest to be created
293/// - BoostType the boosting type for the trees in the forest (AdaBoost e.t.c..).
294/// Known:
295/// - AdaBoost
296/// - AdaBoostR2 (Adaboost for regression)
297/// - Bagging
298/// - GradBoost
299/// - AdaBoostBeta the boosting parameter, beta, for AdaBoost
300/// - UseRandomisedTrees choose at each node splitting a random set of variables
301/// - UseNvars use UseNvars variables in randomised trees
302/// - UsePoisson Nvars use UseNvars not as fixed number but as mean of a poisson distribution
303/// - SeparationType the separation criterion applied in the node splitting.
304/// Known:
305/// - GiniIndex
306/// - MisClassificationError
307/// - CrossEntropy
308/// - SDivSqrtSPlusB
309/// - MinNodeSize: minimum percentage of training events in a leaf node (leaf criteria, stop splitting)
310/// - nCuts: the number of steps in the optimisation of the cut for a node (if < 0, then
311/// step size is determined by the events)
312/// - UseFisherCuts: use multivariate splits using the Fisher criterion
313/// - UseYesNoLeaf decide if the classification is done simply by the node type, or the S/B
314/// (from the training) in the leaf node
315/// - NodePurityLimit the minimum purity to classify a node as a signal node (used in pruning and boosting to determine
316/// misclassification error rate)
317/// - PruneMethod The Pruning method.
318/// Known:
319/// - NoPruning // switch off pruning completely
320/// - ExpectedError
321/// - CostComplexity
322/// - PruneStrength a parameter to adjust the amount of pruning. Should be large enough such that overtraining is avoided.
323/// - PruningValFraction number of events to use for optimizing pruning (only if PruneStrength < 0, i.e. automatic pruning)
324/// - NegWeightTreatment
325/// - IgnoreNegWeightsInTraining Ignore negative weight events in the training.
326/// - DecreaseBoostWeight Boost ev. with neg. weight with 1/boostweight instead of boostweight
327/// - PairNegWeightsGlobal Pair ev. with neg. and pos. weights in training sample and "annihilate" them
328/// - MaxDepth maximum depth of the decision tree allowed before further splitting is stopped
329/// - SkipNormalization Skip normalization at initialization, to keep expectation value of BDT output
330/// according to the fraction of events
331
333{
334 DeclareOptionRef(fNTrees, "NTrees", "Number of trees in the forest");
335 if (DoRegression()) {
336 DeclareOptionRef(fMaxDepth=50,"MaxDepth","Max depth of the decision tree allowed");
337 }else{
338 DeclareOptionRef(fMaxDepth=3,"MaxDepth","Max depth of the decision tree allowed");
339 }
340
341 TString tmp="5%"; if (DoRegression()) tmp="0.2%";
342 DeclareOptionRef(fMinNodeSizeS=tmp, "MinNodeSize", "Minimum percentage of training events required in a leaf node (default: Classification: 5%, Regression: 0.2%)");
343 // MinNodeSize: minimum percentage of training events in a leaf node (leaf criteria, stop splitting)
344 DeclareOptionRef(fNCuts, "nCuts", "Number of grid points in variable range used in finding optimal cut in node splitting");
345
346 DeclareOptionRef(fBoostType, "BoostType", "Boosting type for the trees in the forest (note: AdaCost is still experimental)");
347
348 AddPreDefVal(TString("AdaBoost"));
349 AddPreDefVal(TString("RealAdaBoost"));
350 AddPreDefVal(TString("AdaCost"));
351 AddPreDefVal(TString("Bagging"));
352 // AddPreDefVal(TString("RegBoost"));
353 AddPreDefVal(TString("AdaBoostR2"));
354 AddPreDefVal(TString("Grad"));
355 if (DoRegression()) {
356 fBoostType = "AdaBoostR2";
357 }else{
358 fBoostType = "AdaBoost";
359 }
360 DeclareOptionRef(fAdaBoostR2Loss="Quadratic", "AdaBoostR2Loss", "Type of Loss function in AdaBoostR2");
361 AddPreDefVal(TString("Linear"));
362 AddPreDefVal(TString("Quadratic"));
363 AddPreDefVal(TString("Exponential"));
364
365 DeclareOptionRef(fBaggedBoost=kFALSE, "UseBaggedBoost","Use only a random subsample of all events for growing the trees in each boost iteration.");
366 DeclareOptionRef(fShrinkage = 1.0, "Shrinkage", "Learning rate for BoostType=Grad algorithm");
367 DeclareOptionRef(fAdaBoostBeta=.5, "AdaBoostBeta", "Learning rate for AdaBoost algorithm");
368 DeclareOptionRef(fRandomisedTrees,"UseRandomisedTrees","Determine at each node splitting the cut variable only as the best out of a random subset of variables (like in RandomForests)");
369 DeclareOptionRef(fUseNvars,"UseNvars","Size of the subset of variables used with RandomisedTree option");
370 DeclareOptionRef(fUsePoissonNvars,"UsePoissonNvars", "Interpret \"UseNvars\" not as fixed number but as mean of a Poisson distribution in each split with RandomisedTree option");
371 DeclareOptionRef(fBaggedSampleFraction=.6,"BaggedSampleFraction","Relative size of bagged event sample to original size of the data sample (used whenever bagging is used (i.e. UseBaggedBoost, Bagging,)" );
372
373 DeclareOptionRef(fUseYesNoLeaf=kTRUE, "UseYesNoLeaf",
374 "Use Sig or Bkg categories, or the purity=S/(S+B) as classification of the leaf node -> Real-AdaBoost");
375 if (DoRegression()) {
376 fUseYesNoLeaf = kFALSE;
377 }
378
379 DeclareOptionRef(fNegWeightTreatment="InverseBoostNegWeights","NegWeightTreatment","How to treat events with negative weights in the BDT training (particular the boosting) : IgnoreInTraining; Boost With inverse boostweight; Pair events with negative and positive weights in training sample and *annihilate* them (experimental!)");
380 AddPreDefVal(TString("InverseBoostNegWeights"));
381 AddPreDefVal(TString("IgnoreNegWeightsInTraining"));
382 AddPreDefVal(TString("NoNegWeightsInTraining")); // well, let's be nice to users and keep at least this old name anyway ..
383 AddPreDefVal(TString("PairNegWeightsGlobal"));
384 AddPreDefVal(TString("Pray"));
385
386
387
388 DeclareOptionRef(fCss=1., "Css", "AdaCost: cost of true signal selected signal");
389 DeclareOptionRef(fCts_sb=1.,"Cts_sb","AdaCost: cost of true signal selected bkg");
390 DeclareOptionRef(fCtb_ss=1.,"Ctb_ss","AdaCost: cost of true bkg selected signal");
391 DeclareOptionRef(fCbb=1., "Cbb", "AdaCost: cost of true bkg selected bkg ");
392
393 DeclareOptionRef(fNodePurityLimit=0.5, "NodePurityLimit", "In boosting/pruning, nodes with purity > NodePurityLimit are signal; background otherwise.");
394
395
396 DeclareOptionRef(fSepTypeS, "SeparationType", "Separation criterion for node splitting");
397 AddPreDefVal(TString("CrossEntropy"));
398 AddPreDefVal(TString("GiniIndex"));
399 AddPreDefVal(TString("GiniIndexWithLaplace"));
400 AddPreDefVal(TString("MisClassificationError"));
401 AddPreDefVal(TString("SDivSqrtSPlusB"));
402 AddPreDefVal(TString("RegressionVariance"));
403 if (DoRegression()) {
404 fSepTypeS = "RegressionVariance";
405 }else{
406 fSepTypeS = "GiniIndex";
407 }
408
409 DeclareOptionRef(fRegressionLossFunctionBDTGS = "Huber", "RegressionLossFunctionBDTG", "Loss function for BDTG regression.");
410 AddPreDefVal(TString("Huber"));
411 AddPreDefVal(TString("AbsoluteDeviation"));
412 AddPreDefVal(TString("LeastSquares"));
413
414 DeclareOptionRef(fHuberQuantile = 0.7, "HuberQuantile", "In the Huber loss function this is the quantile that separates the core from the tails in the residuals distribution.");
415
416 DeclareOptionRef(fDoBoostMonitor=kFALSE,"DoBoostMonitor","Create control plot with ROC integral vs tree number");
417
418 DeclareOptionRef(fUseFisherCuts=kFALSE, "UseFisherCuts", "Use multivariate splits using the Fisher criterion");
419 DeclareOptionRef(fMinLinCorrForFisher=.8,"MinLinCorrForFisher", "The minimum linear correlation between two variables demanded for use in Fisher criterion in node splitting");
420 DeclareOptionRef(fUseExclusiveVars=kFALSE,"UseExclusiveVars","Variables already used in fisher criterion are not anymore analysed individually for node splitting");
421
422
423 DeclareOptionRef(fDoPreselection=kFALSE,"DoPreselection","and and apply automatic pre-selection for 100% efficient signal (bkg) cuts prior to training");
424
425
426 DeclareOptionRef(fSigToBkgFraction=1,"SigToBkgFraction","Sig to Bkg ratio used in Training (similar to NodePurityLimit, which cannot be used in real adaboost");
427
428 DeclareOptionRef(fPruneMethodS, "PruneMethod", "Note: for BDTs use small trees (e.g.MaxDepth=3) and NoPruning: Pruning: Method used for pruning (removal) of statistically insignificant branches ");
429 AddPreDefVal(TString("NoPruning"));
430 AddPreDefVal(TString("ExpectedError"));
431 AddPreDefVal(TString("CostComplexity"));
432
433 DeclareOptionRef(fPruneStrength, "PruneStrength", "Pruning strength");
434
435 DeclareOptionRef(fFValidationEvents=0.5, "PruningValFraction", "Fraction of events to use for optimizing automatic pruning.");
436
437 DeclareOptionRef(fSkipNormalization=kFALSE, "SkipNormalization", "Skip normalization at initialization, to keep expectation value of BDT output according to the fraction of events");
438
439 // deprecated options, still kept for the moment:
440 DeclareOptionRef(fMinNodeEvents=0, "nEventsMin", "deprecated: Use MinNodeSize (in % of training events) instead");
441
442 DeclareOptionRef(fBaggedGradBoost=kFALSE, "UseBaggedGrad","deprecated: Use *UseBaggedBoost* instead: Use only a random subsample of all events for growing the trees in each iteration.");
443 DeclareOptionRef(fBaggedSampleFraction, "GradBaggingFraction","deprecated: Use *BaggedSampleFraction* instead: Defines the fraction of events to be used in each iteration, e.g. when UseBaggedGrad=kTRUE. ");
444 DeclareOptionRef(fUseNTrainEvents,"UseNTrainEvents","deprecated: Use *BaggedSampleFraction* instead: Number of randomly picked training events used in randomised (and bagged) trees");
445 DeclareOptionRef(fNNodesMax,"NNodesMax","deprecated: Use MaxDepth instead to limit the tree size" );
446
447
448}
449
450////////////////////////////////////////////////////////////////////////////////
451/// Options that are used ONLY for the READER to ensure backward compatibility.
452
455
456
457 DeclareOptionRef(fHistoricBool=kTRUE, "UseWeightedTrees",
458 "Use weighted trees or simple average in classification from the forest");
459 DeclareOptionRef(fHistoricBool=kFALSE, "PruneBeforeBoost", "Flag to prune the tree before applying boosting algorithm");
460 DeclareOptionRef(fHistoricBool=kFALSE,"RenormByClass","Individually re-normalize each event class to the original size after boosting");
461
462 AddPreDefVal(TString("NegWeightTreatment"),TString("IgnoreNegWeights"));
463
464}
465
466////////////////////////////////////////////////////////////////////////////////
467/// The option string is decoded, for available options see "DeclareOptions".
468
470{
471 fSepTypeS.ToLower();
472 if (fSepTypeS == "misclassificationerror") fSepType = new MisClassificationError();
473 else if (fSepTypeS == "giniindex") fSepType = new GiniIndex();
474 else if (fSepTypeS == "giniindexwithlaplace") fSepType = new GiniIndexWithLaplace();
475 else if (fSepTypeS == "crossentropy") fSepType = new CrossEntropy();
476 else if (fSepTypeS == "sdivsqrtsplusb") fSepType = new SdivSqrtSplusB();
477 else if (fSepTypeS == "regressionvariance") fSepType = NULL;
478 else {
479 Log() << kINFO << GetOptions() << Endl;
480 Log() << kFATAL << "<ProcessOptions> unknown Separation Index option " << fSepTypeS << " called" << Endl;
481 }
482
483 if(!(fHuberQuantile >= 0.0 && fHuberQuantile <= 1.0)){
484 Log() << kINFO << GetOptions() << Endl;
485 Log() << kFATAL << "<ProcessOptions> Huber Quantile must be in range [0,1]. Value given, " << fHuberQuantile << ", does not match this criteria" << Endl;
486 }
487
488
489 fRegressionLossFunctionBDTGS.ToLower();
490 if (fRegressionLossFunctionBDTGS == "huber") fRegressionLossFunctionBDTG = new HuberLossFunctionBDT(fHuberQuantile);
491 else if (fRegressionLossFunctionBDTGS == "leastsquares") fRegressionLossFunctionBDTG = new LeastSquaresLossFunctionBDT();
492 else if (fRegressionLossFunctionBDTGS == "absolutedeviation") fRegressionLossFunctionBDTG = new AbsoluteDeviationLossFunctionBDT();
493 else {
494 Log() << kINFO << GetOptions() << Endl;
495 Log() << kFATAL << "<ProcessOptions> unknown Regression Loss Function BDT option " << fRegressionLossFunctionBDTGS << " called" << Endl;
496 }
497
498 fPruneMethodS.ToLower();
499 if (fPruneMethodS == "expectederror") fPruneMethod = DecisionTree::kExpectedErrorPruning;
500 else if (fPruneMethodS == "costcomplexity") fPruneMethod = DecisionTree::kCostComplexityPruning;
501 else if (fPruneMethodS == "nopruning") fPruneMethod = DecisionTree::kNoPruning;
502 else {
503 Log() << kINFO << GetOptions() << Endl;
504 Log() << kFATAL << "<ProcessOptions> unknown PruneMethod " << fPruneMethodS << " option called" << Endl;
505 }
506 if (fPruneStrength < 0 && (fPruneMethod != DecisionTree::kNoPruning) && fBoostType!="Grad") fAutomatic = kTRUE;
507 else fAutomatic = kFALSE;
508 if (fAutomatic && fPruneMethod==DecisionTree::kExpectedErrorPruning){
509 Log() << kFATAL
510 << "Sorry automatic pruning strength determination is not implemented yet for ExpectedErrorPruning" << Endl;
511 }
512
513
514 if (fMinNodeEvents > 0){
515 fMinNodeSize = Double_t(fMinNodeEvents*100.) / Data()->GetNTrainingEvents();
516 Log() << kWARNING << "You have explicitly set ** nEventsMin = " << fMinNodeEvents<<" ** the min absolute number \n"
517 << "of events in a leaf node. This is DEPRECATED, please use the option \n"
518 << "*MinNodeSize* giving the relative number as percentage of training \n"
519 << "events instead. \n"
520 << "nEventsMin="<<fMinNodeEvents<< "--> MinNodeSize="<<fMinNodeSize<<"%"
521 << Endl;
522 Log() << kWARNING << "Note also that explicitly setting *nEventsMin* so far OVERWRITES the option recommended \n"
523 << " *MinNodeSize* = " << fMinNodeSizeS << " option !!" << Endl ;
524 fMinNodeSizeS = TString::Format("%F3.2",fMinNodeSize);
525
526 }else{
527 SetMinNodeSize(fMinNodeSizeS);
528 }
529
530
531 fAdaBoostR2Loss.ToLower();
532
533 if (fBoostType=="Grad") {
534 fPruneMethod = DecisionTree::kNoPruning;
535 if (fNegWeightTreatment=="InverseBoostNegWeights"){
536 Log() << kINFO << "the option NegWeightTreatment=InverseBoostNegWeights does"
537 << " not exist for BoostType=Grad" << Endl;
538 Log() << kINFO << "--> change to new default NegWeightTreatment=Pray" << Endl;
539 Log() << kDEBUG << "i.e. simply keep them as if which should work fine for Grad Boost" << Endl;
540 fNegWeightTreatment="Pray";
541 fNoNegWeightsInTraining=kFALSE;
542 }
543 } else if (fBoostType=="RealAdaBoost"){
544 fBoostType = "AdaBoost";
545 fUseYesNoLeaf = kFALSE;
546 } else if (fBoostType=="AdaCost"){
547 fUseYesNoLeaf = kFALSE;
548 }
549
550 if (fFValidationEvents < 0.0) fFValidationEvents = 0.0;
551 if (fAutomatic && fFValidationEvents > 0.5) {
552 Log() << kWARNING << "You have chosen to use more than half of your training sample "
553 << "to optimize the automatic pruning algorithm. This is probably wasteful "
554 << "and your overall results will be degraded. Are you sure you want this?"
555 << Endl;
556 }
557
558
559 if (this->Data()->HasNegativeEventWeights()){
560 Log() << kINFO << " You are using a Monte Carlo that has also negative weights. "
561 << "That should in principle be fine as long as on average you end up with "
562 << "something positive. For this you have to make sure that the minimal number "
563 << "of (un-weighted) events demanded for a tree node (currently you use: MinNodeSize="
564 << fMinNodeSizeS << " ("<< fMinNodeSize << "%)"
565 <<", (or the deprecated equivalent nEventsMin) you can set this via the "
566 <<"BDT option string when booking the "
567 << "classifier) is large enough to allow for reasonable averaging!!! "
568 << " If this does not help.. maybe you want to try the option: IgnoreNegWeightsInTraining "
569 << "which ignores events with negative weight in the training. " << Endl
570 << Endl << "Note: You'll get a WARNING message during the training if that should ever happen" << Endl;
571 }
572
573 if (DoRegression()) {
574 if (fUseYesNoLeaf && !IsConstructedFromWeightFile()){
575 Log() << kWARNING << "Regression Trees do not work with fUseYesNoLeaf=TRUE --> I will set it to FALSE" << Endl;
576 fUseYesNoLeaf = kFALSE;
577 }
578
579 if (fSepType != NULL){
580 Log() << kWARNING << "Regression Trees do not work with Separation type other than <RegressionVariance> --> I will use it instead" << Endl;
581 fSepType = NULL;
582 }
583 if (fUseFisherCuts){
584 Log() << kWARNING << "Sorry, UseFisherCuts is not available for regression analysis, I will ignore it!" << Endl;
585 fUseFisherCuts = kFALSE;
586 }
587 if (fNCuts < 0) {
588 Log() << kWARNING << "Sorry, the option of nCuts<0 using a more elaborate node splitting algorithm " << Endl;
589 Log() << kWARNING << "is not implemented for regression analysis ! " << Endl;
590 Log() << kWARNING << "--> I switch do default nCuts = 20 and use standard node splitting"<<Endl;
591 fNCuts=20;
592 }
593 }
594 if (fRandomisedTrees){
595 Log() << kINFO << " Randomised trees use no pruning" << Endl;
596 fPruneMethod = DecisionTree::kNoPruning;
597 // fBoostType = "Bagging";
598 }
599
600 if (fUseFisherCuts) {
601 Log() << kWARNING << "When using the option UseFisherCuts, the other option nCuts<0 (i.e. using" << Endl;
602 Log() << " a more elaborate node splitting algorithm) is not implemented. " << Endl;
603 //I will switch o " << Endl;
604 //Log() << "--> I switch do default nCuts = 20 and use standard node splitting WITH possible Fisher criteria"<<Endl;
605 fNCuts=20;
606 }
607
608 if (fNTrees==0){
609 Log() << kERROR << " Zero Decision Trees demanded... that does not work !! "
610 << " I set it to 1 .. just so that the program does not crash"
611 << Endl;
612 fNTrees = 1;
613 }
614
615 fNegWeightTreatment.ToLower();
616 if (fNegWeightTreatment == "ignorenegweightsintraining") fNoNegWeightsInTraining = kTRUE;
617 else if (fNegWeightTreatment == "nonegweightsintraining") fNoNegWeightsInTraining = kTRUE;
618 else if (fNegWeightTreatment == "inverseboostnegweights") fInverseBoostNegWeights = kTRUE;
619 else if (fNegWeightTreatment == "pairnegweightsglobal") fPairNegWeightsGlobal = kTRUE;
620 else if (fNegWeightTreatment == "pray") Log() << kDEBUG << "Yes, good luck with praying " << Endl;
621 else {
622 Log() << kINFO << GetOptions() << Endl;
623 Log() << kFATAL << "<ProcessOptions> unknown option for treating negative event weights during training " << fNegWeightTreatment << " requested" << Endl;
624 }
625
626 if (fNegWeightTreatment == "pairnegweightsglobal")
627 Log() << kWARNING << " you specified the option NegWeightTreatment=PairNegWeightsGlobal : This option is still considered EXPERIMENTAL !! " << Endl;
628
629
630 // dealing with deprecated options !
631 if (fNNodesMax>0) {
632 UInt_t tmp=1; // depth=0 == 1 node
633 fMaxDepth=0;
634 while (tmp < fNNodesMax){
635 tmp+=2*tmp;
636 fMaxDepth++;
637 }
638 Log() << kWARNING << "You have specified a deprecated option *NNodesMax="<<fNNodesMax
639 << "* \n this has been translated to MaxDepth="<<fMaxDepth<<Endl;
640 }
641
642
643 if (fUseNTrainEvents>0){
644 fBaggedSampleFraction = (Double_t) fUseNTrainEvents/Data()->GetNTrainingEvents();
645 Log() << kWARNING << "You have specified a deprecated option *UseNTrainEvents="<<fUseNTrainEvents
646 << "* \n this has been translated to BaggedSampleFraction="<<fBaggedSampleFraction<<"(%)"<<Endl;
647 }
648
649 if (fBoostType=="Bagging") fBaggedBoost = kTRUE;
650 if (fBaggedGradBoost){
651 fBaggedBoost = kTRUE;
652 Log() << kWARNING << "You have specified a deprecated option *UseBaggedGrad* --> please use *UseBaggedBoost* instead" << Endl;
653 }
654
655}
656
657////////////////////////////////////////////////////////////////////////////////
658
660 if (sizeInPercent > 0 && sizeInPercent < 50){
661 fMinNodeSize=sizeInPercent;
662
663 } else {
664 Log() << kFATAL << "you have demanded a minimal node size of "
665 << sizeInPercent << "% of the training events.. \n"
666 << " that somehow does not make sense "<<Endl;
667 }
668
669}
670
671////////////////////////////////////////////////////////////////////////////////
672
674 sizeInPercent.ReplaceAll("%","");
675 sizeInPercent.ReplaceAll(" ","");
676 if (sizeInPercent.IsFloat()) SetMinNodeSize(sizeInPercent.Atof());
677 else {
678 Log() << kFATAL << "I had problems reading the option MinNodeEvents, which "
679 << "after removing a possible % sign now reads " << sizeInPercent << Endl;
680 }
681}
682
683////////////////////////////////////////////////////////////////////////////////
684/// Common initialisation with defaults for the BDT-Method.
685
687{
688 fNTrees = 800;
689 if (fAnalysisType == Types::kClassification || fAnalysisType == Types::kMulticlass ) {
690 fMaxDepth = 3;
691 fBoostType = "AdaBoost";
692 if(DataInfo().GetNClasses()!=0) //workaround for multiclass application
693 fMinNodeSize = 5.;
694 }else {
695 fMaxDepth = 50;
696 fBoostType = "AdaBoostR2";
697 fAdaBoostR2Loss = "Quadratic";
698 if(DataInfo().GetNClasses()!=0) //workaround for multiclass application
699 fMinNodeSize = .2;
700 }
701
702
703 fNCuts = 20;
704 fPruneMethodS = "NoPruning";
705 fPruneMethod = DecisionTree::kNoPruning;
706 fPruneStrength = 0;
707 fAutomatic = kFALSE;
708 fFValidationEvents = 0.5;
709 fRandomisedTrees = kFALSE;
710 // fUseNvars = (GetNvar()>12) ? UInt_t(GetNvar()/8) : TMath::Max(UInt_t(2),UInt_t(GetNvar()/3));
711 fUseNvars = UInt_t(TMath::Sqrt(GetNvar())+0.6);
712 fUsePoissonNvars = kTRUE;
713 fShrinkage = 1.0;
714// fSumOfWeights = 0.0;
715
716 // reference cut value to distinguish signal-like from background-like events
717 SetSignalReferenceCut( 0 );
718}
719
720
721////////////////////////////////////////////////////////////////////////////////
722/// Reset the method, as if it had just been instantiated (forget all training etc.).
723
725{
726 // I keep the BDT EventSample and its Validation sample (eventually they should all
727 // disappear and just use the DataSet samples ..
728
729 // remove all the trees
730 for (UInt_t i=0; i<fForest.size(); i++) delete fForest[i];
731 fForest.clear();
732
733 fBoostWeights.clear();
734 if (fMonitorNtuple) { fMonitorNtuple->Delete(); fMonitorNtuple=NULL; }
735 fVariableImportance.clear();
736 fResiduals.clear();
737 fLossFunctionEventInfo.clear();
738 // now done in "InitEventSample" which is called in "Train"
739 // reset all previously stored/accumulated BOOST weights in the event sample
740 //for (UInt_t iev=0; iev<fEventSample.size(); iev++) fEventSample[iev]->SetBoostWeight(1.);
741 if (Data()) Data()->DeleteResults(GetMethodName(), Types::kTraining, GetAnalysisType());
742 Log() << kDEBUG << " successfully(?) reset the method " << Endl;
743}
744
745
746////////////////////////////////////////////////////////////////////////////////
747/// Destructor.
748///
749/// - Note: fEventSample and ValidationSample are already deleted at the end of TRAIN
750/// When they are not used anymore
751
753{
754 for (UInt_t i=0; i<fForest.size(); i++) delete fForest[i];
755}
756
757////////////////////////////////////////////////////////////////////////////////
758/// Initialize the event sample (i.e. reset the boost-weights... etc).
759
761{
762 if (!HasTrainingTree()) Log() << kFATAL << "<Init> Data().TrainingTree() is zero pointer" << Endl;
763
764 if (fEventSample.size() > 0) { // do not re-initialise the event sample, just set all boostweights to 1. as if it were untouched
765 // reset all previously stored/accumulated BOOST weights in the event sample
766 for (UInt_t iev=0; iev<fEventSample.size(); iev++) fEventSample[iev]->SetBoostWeight(1.);
767 } else {
768 Data()->SetCurrentType(Types::kTraining);
769 UInt_t nevents = Data()->GetNTrainingEvents();
770
771 std::vector<const TMVA::Event*> tmpEventSample;
772 for (Long64_t ievt=0; ievt<nevents; ievt++) {
773 // const Event *event = new Event(*(GetEvent(ievt)));
774 Event* event = new Event( *GetTrainingEvent(ievt) );
775 tmpEventSample.push_back(event);
776 }
777
778 if (!DoRegression()) DeterminePreselectionCuts(tmpEventSample);
779 else fDoPreselection = kFALSE; // just to make sure...
780
781 for (UInt_t i=0; i<tmpEventSample.size(); i++) delete tmpEventSample[i];
782
783
786 for (Long64_t ievt=0; ievt<nevents; ievt++) {
787 // const Event *event = new Event(*(GetEvent(ievt)));
788 // const Event* event = new Event( *GetTrainingEvent(ievt) );
789 Event* event = new Event( *GetTrainingEvent(ievt) );
790 if (fDoPreselection){
791 if (TMath::Abs(ApplyPreselectionCuts(event)) > 0.05) {
792 delete event;
793 continue;
794 }
795 }
796
797 if (event->GetWeight() < 0 && (IgnoreEventsWithNegWeightsInTraining() || fNoNegWeightsInTraining)){
798 if (firstNegWeight) {
799 Log() << kWARNING << " Note, you have events with negative event weight in the sample, but you've chosen to ignore them" << Endl;
801 }
802 delete event;
803 }else if (event->GetWeight()==0){
804 if (firstZeroWeight) {
806 Log() << "Events with weight == 0 are going to be simply ignored " << Endl;
807 }
808 delete event;
809 }else{
810 if (event->GetWeight() < 0) {
811 fTrainWithNegWeights=kTRUE;
812 if (firstNegWeight){
814 if (fPairNegWeightsGlobal){
815 Log() << kWARNING << "Events with negative event weights are found and "
816 << " will be removed prior to the actual BDT training by global "
817 << " paring (and subsequent annihilation) with positiv weight events"
818 << Endl;
819 }else{
820 Log() << kWARNING << "Events with negative event weights are USED during "
821 << "the BDT training. This might cause problems with small node sizes "
822 << "or with the boosting. Please remove negative events from training "
823 << "using the option *IgnoreEventsWithNegWeightsInTraining* in case you "
824 << "observe problems with the boosting"
825 << Endl;
826 }
827 }
828 }
829 // if fAutomatic == true you need a validation sample to optimize pruning
830 if (fAutomatic) {
831 Double_t modulo = 1.0/(fFValidationEvents);
832 Int_t imodulo = static_cast<Int_t>( fmod(modulo,1.0) > 0.5 ? ceil(modulo) : floor(modulo) );
833 if (ievt % imodulo == 0) fValidationSample.push_back( event );
834 else fEventSample.push_back( event );
835 }
836 else {
837 fEventSample.push_back(event);
838 }
839 }
840 }
841
842 if (fAutomatic) {
843 Log() << kINFO << "<InitEventSample> Internally I use " << fEventSample.size()
844 << " for Training and " << fValidationSample.size()
845 << " for Pruning Validation (" << ((Float_t)fValidationSample.size())/((Float_t)fEventSample.size()+fValidationSample.size())*100.0
846 << "% of training used for validation)" << Endl;
847 }
848
849 // some pre-processing for events with negative weights
850 if (fPairNegWeightsGlobal) PreProcessNegativeEventWeights();
851 }
852
853 if (DoRegression()) {
854 // Regression, no reweighting to do
855 } else if (DoMulticlass()) {
856 // Multiclass, only gradboost is supported. No reweighting.
857 } else if (!fSkipNormalization) {
858 // Binary classification.
859 Log() << kDEBUG << "\t<InitEventSample> For classification trees, "<< Endl;
860 Log() << kDEBUG << " \tthe effective number of backgrounds is scaled to match "<<Endl;
861 Log() << kDEBUG << " \tthe signal. Otherwise the first boosting step would do 'just that'!"<<Endl;
862 // it does not make sense in decision trees to start with unequal number of signal/background
863 // events (weights) .. hence normalize them now (happens otherwise in first 'boosting step'
864 // anyway..
865 // Also make sure, that the sum_of_weights == sample.size() .. as this is assumed in
866 // the DecisionTree to derive a sensible number for "fMinSize" (min.#events in node)
867 // that currently is an OR between "weighted" and "unweighted number"
868 // I want:
869 // nS + nB = n
870 // a*SW + b*BW = n
871 // (a*SW)/(b*BW) = fSigToBkgFraction
872 //
873 // ==> b = n/((1+f)BW) and a = (nf/(1+f))/SW
874
875 Double_t nevents = fEventSample.size();
877 Int_t sumSig=0, sumBkg=0;
878 for (UInt_t ievt=0; ievt<fEventSample.size(); ievt++) {
879 if ((DataInfo().IsSignal(fEventSample[ievt])) ) {
880 sumSigW += fEventSample[ievt]->GetWeight();
881 sumSig++;
882 } else {
883 sumBkgW += fEventSample[ievt]->GetWeight();
884 sumBkg++;
885 }
886 }
887 if (sumSigW && sumBkgW){
888 Double_t normSig = nevents/((1+fSigToBkgFraction)*sumSigW)*fSigToBkgFraction;
889 Double_t normBkg = nevents/((1+fSigToBkgFraction)*sumBkgW); ;
890 Log() << kDEBUG << "\tre-normalise events such that Sig and Bkg have respective sum of weights = "
891 << fSigToBkgFraction << Endl;
892 Log() << kDEBUG << " \tsig->sig*"<<normSig << "ev. bkg->bkg*"<<normBkg << "ev." <<Endl;
893 Log() << kHEADER << "#events: (reweighted) sig: "<< sumSigW*normSig << " bkg: " << sumBkgW*normBkg << Endl;
894 Log() << kINFO << "#events: (unweighted) sig: "<< sumSig << " bkg: " << sumBkg << Endl;
895 for (Long64_t ievt=0; ievt<nevents; ievt++) {
896 if ((DataInfo().IsSignal(fEventSample[ievt])) ) fEventSample[ievt]->SetBoostWeight(normSig);
897 else fEventSample[ievt]->SetBoostWeight(normBkg);
898 }
899 }else{
900 Log() << kINFO << "--> could not determine scaling factors as either there are " << Endl;
901 Log() << kINFO << " no signal events (sumSigW="<<sumSigW<<") or no bkg ev. (sumBkgW="<<sumBkgW<<")"<<Endl;
902 }
903
904 }
905
906 fTrainSample = &fEventSample;
907 if (fBaggedBoost){
908 GetBaggedSubSample(fEventSample);
909 fTrainSample = &fSubSample;
910 }
911
912 //just for debug purposes..
913 /*
914 sumSigW=0;
915 sumBkgW=0;
916 for (UInt_t ievt=0; ievt<fEventSample.size(); ievt++) {
917 if ((DataInfo().IsSignal(fEventSample[ievt])) ) sumSigW += fEventSample[ievt]->GetWeight();
918 else sumBkgW += fEventSample[ievt]->GetWeight();
919 }
920 Log() << kWARNING << "sigSumW="<<sumSigW<<"bkgSumW="<<sumBkgW<< Endl;
921 */
922}
923
924////////////////////////////////////////////////////////////////////////////////
925/// O.k. you know there are events with negative event weights. This routine will remove
926/// them by pairing them with the closest event(s) of the same event class with positive
927/// weights
928/// A first attempt is "brute force", I dont' try to be clever using search trees etc,
929/// just quick and dirty to see if the result is any good
930
935 std::vector<const Event*> negEvents;
936 for (UInt_t iev = 0; iev < fEventSample.size(); iev++){
937 if (fEventSample[iev]->GetWeight() < 0) {
938 totalNegWeights += fEventSample[iev]->GetWeight();
939 negEvents.push_back(fEventSample[iev]);
940 } else {
941 totalPosWeights += fEventSample[iev]->GetWeight();
942 }
943 totalWeights += fEventSample[iev]->GetWeight();
944 }
945 if (totalNegWeights == 0 ) {
946 Log() << kINFO << "no negative event weights found .. no preprocessing necessary" << Endl;
947 return;
948 } else {
949 Log() << kINFO << "found a total of " << totalNegWeights << " of negative event weights which I am going to try to pair with positive events to annihilate them" << Endl;
950 Log() << kINFO << "found a total of " << totalPosWeights << " of events with positive weights" << Endl;
951 Log() << kINFO << "--> total sum of weights = " << totalWeights << " = " << totalNegWeights+totalPosWeights << Endl;
952 }
953
954 std::vector<TMatrixDSym*>* cov = gTools().CalcCovarianceMatrices( fEventSample, 2);
955
957
958 for (Int_t i=0; i<2; i++){
959 invCov = ((*cov)[i]);
960 if ( TMath::Abs(invCov->Determinant()) < 10E-24 ) {
961 std::cout << "<MethodBDT::PreProcessNeg...> matrix is almost singular with determinant="
962 << TMath::Abs(invCov->Determinant())
963 << " did you use the variables that are linear combinations or highly correlated?"
964 << std::endl;
965 }
966 if ( TMath::Abs(invCov->Determinant()) < 10E-120 ) {
967 std::cout << "<MethodBDT::PreProcessNeg...> matrix is singular with determinant="
968 << TMath::Abs(invCov->Determinant())
969 << " did you use the variables that are linear combinations?"
970 << std::endl;
971 }
972
973 invCov->Invert();
974 }
975
976
977
978 Log() << kINFO << "Found a total of " << totalNegWeights << " in negative weights out of " << fEventSample.size() << " training events " << Endl;
979 Timer timer(negEvents.size(),"Negative Event paired");
980 for (UInt_t nev = 0; nev < negEvents.size(); nev++){
981 timer.DrawProgressBar( nev );
982 Double_t weight = negEvents[nev]->GetWeight();
983 UInt_t iClassID = negEvents[nev]->GetClass();
984 invCov = ((*cov)[iClassID]);
985 while (weight < 0){
986 // find closest event with positive event weight and "pair" it with the negative event
987 // (add their weight) until there is no negative weight anymore
988 Int_t iMin=-1;
989 Double_t dist, minDist=10E270;
990 for (UInt_t iev = 0; iev < fEventSample.size(); iev++){
991 if (iClassID==fEventSample[iev]->GetClass() && fEventSample[iev]->GetWeight() > 0){
992 dist=0;
993 for (UInt_t ivar=0; ivar < GetNvar(); ivar++){
994 for (UInt_t jvar=0; jvar<GetNvar(); jvar++){
995 dist += (negEvents[nev]->GetValue(ivar)-fEventSample[iev]->GetValue(ivar))*
996 (*invCov)[ivar][jvar]*
997 (negEvents[nev]->GetValue(jvar)-fEventSample[iev]->GetValue(jvar));
998 }
999 }
1000 if (dist < minDist) { iMin=iev; minDist=dist;}
1001 }
1002 }
1003
1004 if (iMin > -1) {
1005 // std::cout << "Happily pairing .. weight before : " << negEvents[nev]->GetWeight() << " and " << fEventSample[iMin]->GetWeight();
1006 Double_t newWeight = (negEvents[nev]->GetWeight() + fEventSample[iMin]->GetWeight());
1007 if (newWeight > 0){
1008 negEvents[nev]->SetBoostWeight( 0 );
1009 fEventSample[iMin]->SetBoostWeight( newWeight/fEventSample[iMin]->GetOriginalWeight() ); // note the weight*boostweight should be "newWeight"
1010 } else {
1011 negEvents[nev]->SetBoostWeight( newWeight/negEvents[nev]->GetOriginalWeight() ); // note the weight*boostweight should be "newWeight"
1012 fEventSample[iMin]->SetBoostWeight( 0 );
1013 }
1014 // std::cout << " and afterwards " << negEvents[nev]->GetWeight() << " and the paired " << fEventSample[iMin]->GetWeight() << " dist="<<minDist<< std::endl;
1015 } else Log() << kFATAL << "preprocessing didn't find event to pair with the negative weight ... probably a bug" << Endl;
1016 weight = negEvents[nev]->GetWeight();
1017 }
1018 }
1019 Log() << kINFO << "<Negative Event Pairing> took: " << timer.GetElapsedTime()
1020 << " " << Endl;
1021
1022 // just check.. now there should be no negative event weight left anymore
1023 totalNegWeights = 0;
1024 totalPosWeights = 0;
1025 totalWeights = 0;
1028 Int_t nSig=0;
1029 Int_t nBkg=0;
1030
1031 std::vector<const Event*> newEventSample;
1032
1033 for (UInt_t iev = 0; iev < fEventSample.size(); iev++){
1034 if (fEventSample[iev]->GetWeight() < 0) {
1035 totalNegWeights += fEventSample[iev]->GetWeight();
1036 totalWeights += fEventSample[iev]->GetWeight();
1037 } else {
1038 totalPosWeights += fEventSample[iev]->GetWeight();
1039 totalWeights += fEventSample[iev]->GetWeight();
1040 }
1041 if (fEventSample[iev]->GetWeight() > 0) {
1042 newEventSample.push_back(new Event(*fEventSample[iev]));
1043 if (fEventSample[iev]->GetClass() == fSignalClass){
1044 sigWeight += fEventSample[iev]->GetWeight();
1045 nSig+=1;
1046 }else{
1047 bkgWeight += fEventSample[iev]->GetWeight();
1048 nBkg+=1;
1049 }
1050 }
1051 }
1052 if (totalNegWeights < 0) Log() << kFATAL << " compensation of negative event weights with positive ones did not work " << totalNegWeights << Endl;
1053
1054 for (UInt_t i=0; i<fEventSample.size(); i++) delete fEventSample[i];
1055 fEventSample = newEventSample;
1056
1057 Log() << kINFO << " after PreProcessing, the Event sample is left with " << fEventSample.size() << " events (unweighted), all with positive weights, adding up to " << totalWeights << Endl;
1058 Log() << kINFO << " nSig="<<nSig << " sigWeight="<<sigWeight << " nBkg="<<nBkg << " bkgWeight="<<bkgWeight << Endl;
1059
1060
1061}
1062
1063////////////////////////////////////////////////////////////////////////////////
1064/// Call the Optimizer with the set of parameters and ranges that
1065/// are meant to be tuned.
1066
1068{
1069 // fill all the tuning parameters that should be optimized into a map:
1070 std::map<TString,TMVA::Interval*> tuneParameters;
1071 std::map<TString,Double_t> tunedParameters;
1072
1073 // note: the 3rd parameter in the interval is the "number of bins", NOT the stepsize !!
1074 // the actual VALUES at (at least for the scan, guess also in GA) are always
1075 // read from the middle of the bins. Hence.. the choice of Intervals e.g. for the
1076 // MaxDepth, in order to make nice integer values!!!
1077
1078 // find some reasonable ranges for the optimisation of MinNodeEvents:
1079
1080 tuneParameters.insert(std::pair<TString,Interval*>("NTrees", new Interval(10,1000,5))); // stepsize 50
1081 tuneParameters.insert(std::pair<TString,Interval*>("MaxDepth", new Interval(2,4,3))); // stepsize 1
1082 tuneParameters.insert(std::pair<TString,Interval*>("MinNodeSize", new LogInterval(1,30,30))); //
1083 //tuneParameters.insert(std::pair<TString,Interval*>("NodePurityLimit",new Interval(.4,.6,3))); // stepsize .1
1084 //tuneParameters.insert(std::pair<TString,Interval*>("BaggedSampleFraction",new Interval(.4,.9,6))); // stepsize .1
1085
1086 // method-specific parameters
1087 if (fBoostType=="AdaBoost"){
1088 tuneParameters.insert(std::pair<TString,Interval*>("AdaBoostBeta", new Interval(.2,1.,5)));
1089
1090 }else if (fBoostType=="Grad"){
1091 tuneParameters.insert(std::pair<TString,Interval*>("Shrinkage", new Interval(0.05,0.50,5)));
1092
1093 }else if (fBoostType=="Bagging" && fRandomisedTrees){
1094 Int_t min_var = TMath::FloorNint( GetNvar() * .25 );
1095 Int_t max_var = TMath::CeilNint( GetNvar() * .75 );
1096 tuneParameters.insert(std::pair<TString,Interval*>("UseNvars", new Interval(min_var,max_var,4)));
1097
1098 }
1099
1100 Log()<<kINFO << " the following BDT parameters will be tuned on the respective *grid*\n"<<Endl;
1101 std::map<TString,TMVA::Interval*>::iterator it;
1102 for(it=tuneParameters.begin(); it!= tuneParameters.end(); ++it){
1103 Log() << kWARNING << it->first << Endl;
1104 std::ostringstream oss;
1105 (it->second)->Print(oss);
1106 Log()<<oss.str();
1107 Log()<<Endl;
1108 }
1109
1111 tunedParameters=optimize.optimize();
1112
1113 return tunedParameters;
1114
1115}
1116
1117////////////////////////////////////////////////////////////////////////////////
1118/// Set the tuning parameters according to the argument.
1119
1121{
1122 std::map<TString,Double_t>::iterator it;
1123 for(it=tuneParameters.begin(); it!= tuneParameters.end(); ++it){
1124 Log() << kWARNING << it->first << " = " << it->second << Endl;
1125 if (it->first == "MaxDepth" ) SetMaxDepth ((Int_t)it->second);
1126 else if (it->first == "MinNodeSize" ) SetMinNodeSize (it->second);
1127 else if (it->first == "NTrees" ) SetNTrees ((Int_t)it->second);
1128 else if (it->first == "NodePurityLimit") SetNodePurityLimit (it->second);
1129 else if (it->first == "AdaBoostBeta" ) SetAdaBoostBeta (it->second);
1130 else if (it->first == "Shrinkage" ) SetShrinkage (it->second);
1131 else if (it->first == "UseNvars" ) SetUseNvars ((Int_t)it->second);
1132 else if (it->first == "BaggedSampleFraction" ) SetBaggedSampleFraction (it->second);
1133 else Log() << kFATAL << " SetParameter for " << it->first << " not yet implemented " <<Endl;
1134 }
1135}
1136
1137////////////////////////////////////////////////////////////////////////////////
1138/// BDT training.
1139
1140
1142{
1144
1145 // fill the STL Vector with the event sample
1146 // (needs to be done here and cannot be done in "init" as the options need to be
1147 // known).
1148 InitEventSample();
1149
1150 if (fNTrees==0){
1151 Log() << kERROR << " Zero Decision Trees demanded... that does not work !! "
1152 << " I set it to 1 .. just so that the program does not crash"
1153 << Endl;
1154 fNTrees = 1;
1155 }
1156
1157 if (fInteractive && fInteractive->NotInitialized()){
1158 std::vector<TString> titles = {"Boost weight", "Error Fraction"};
1159 fInteractive->Init(titles);
1160 }
1161 fIPyMaxIter = fNTrees;
1162 fExitFromTraining = false;
1163
1164 // HHV (it's been here since looong but I really don't know why we cannot handle
1165 // normalized variables in BDTs... todo
1166 if (IsNormalised()) Log() << kFATAL << "\"Normalise\" option cannot be used with BDT; "
1167 << "please remove the option from the configuration string, or "
1168 << "use \"!Normalise\""
1169 << Endl;
1170
1171 if(DoRegression())
1172 Log() << kINFO << "Regression Loss Function: "<< fRegressionLossFunctionBDTG->Name() << Endl;
1173
1174 Log() << kINFO << "Training "<< fNTrees << " Decision Trees ... patience please" << Endl;
1175
1176 Log() << kDEBUG << "Training with maximal depth = " <<fMaxDepth
1177 << ", MinNodeEvents=" << fMinNodeEvents
1178 << ", NTrees="<<fNTrees
1179 << ", NodePurityLimit="<<fNodePurityLimit
1180 << ", AdaBoostBeta="<<fAdaBoostBeta
1181 << Endl;
1182
1183 // weights applied in boosting
1184 Int_t nBins;
1186 TString hname = "AdaBooost weight distribution";
1187
1188 nBins= 100;
1189 xMin = 0;
1190 xMax = 30;
1191
1192 if (DoRegression()) {
1193 nBins= 100;
1194 xMin = 0;
1195 xMax = 1;
1196 hname="Boost event weights distribution";
1197 }
1198
1199 // book monitoring histograms (for AdaBost only)
1200
1201 TH1* h = new TH1F(TString::Format("%s_BoostWeight",DataInfo().GetName()).Data(),hname,nBins,xMin,xMax);
1202 TH1* nodesBeforePruningVsTree = new TH1I(TString::Format("%s_NodesBeforePruning",DataInfo().GetName()).Data(),"nodes before pruning",fNTrees,0,fNTrees);
1203 TH1* nodesAfterPruningVsTree = new TH1I(TString::Format("%s_NodesAfterPruning",DataInfo().GetName()).Data(),"nodes after pruning",fNTrees,0,fNTrees);
1204
1205
1206 if(!DoMulticlass()){
1207 Results* results = Data()->GetResults(GetMethodName(), Types::kTraining, GetAnalysisType());
1208
1209 h->SetXTitle("boost weight");
1210 results->Store(h, "BoostWeights");
1211
1212
1213 // Monitor the performance (on TEST sample) versus number of trees
1214 if (fDoBoostMonitor){
1215 TH2* boostMonitor = new TH2F("BoostMonitor","ROC Integral Vs iTree",2,0,fNTrees,2,0,1.05);
1216 boostMonitor->SetXTitle("#tree");
1217 boostMonitor->SetYTitle("ROC Integral");
1218 results->Store(boostMonitor, "BoostMonitor");
1220 boostMonitorGraph->SetName("BoostMonitorGraph");
1221 boostMonitorGraph->SetTitle("ROCIntegralVsNTrees");
1222 results->Store(boostMonitorGraph, "BoostMonitorGraph");
1223 }
1224
1225 // weights applied in boosting vs tree number
1226 h = new TH1F("BoostWeightVsTree","Boost weights vs tree",fNTrees,0,fNTrees);
1227 h->SetXTitle("#tree");
1228 h->SetYTitle("boost weight");
1229 results->Store(h, "BoostWeightsVsTree");
1230
1231 // error fraction vs tree number
1232 h = new TH1F("ErrFractHist","error fraction vs tree number",fNTrees,0,fNTrees);
1233 h->SetXTitle("#tree");
1234 h->SetYTitle("error fraction");
1235 results->Store(h, "ErrorFrac");
1236
1237 // nNodesBeforePruning vs tree number
1238 nodesBeforePruningVsTree->SetXTitle("#tree");
1239 nodesBeforePruningVsTree->SetYTitle("#tree nodes");
1241
1242 // nNodesAfterPruning vs tree number
1243 nodesAfterPruningVsTree->SetXTitle("#tree");
1244 nodesAfterPruningVsTree->SetYTitle("#tree nodes");
1246
1247 }
1248
1249 fMonitorNtuple= new TTree("MonitorNtuple","BDT variables");
1250 fMonitorNtuple->Branch("iTree",&fITree,"iTree/I");
1251 fMonitorNtuple->Branch("boostWeight",&fBoostWeight,"boostWeight/D");
1252 fMonitorNtuple->Branch("errorFraction",&fErrorFraction,"errorFraction/D");
1253
1254 Timer timer( fNTrees, GetName() );
1257
1260
1261 if(fBoostType=="Grad"){
1262 InitGradBoost(fEventSample);
1263 }
1264
1265 Int_t itree=0;
1267 //for (int itree=0; itree<fNTrees; itree++) {
1268
1269 while (itree < fNTrees && continueBoost){
1270 if (fExitFromTraining) break;
1271 fIPyCurrentIter = itree;
1272 timer.DrawProgressBar( itree );
1273 // Results* results = Data()->GetResults(GetMethodName(), Types::kTraining, GetAnalysisType());
1274 // TH1 *hxx = new TH1F(TString::Format("swdist%d",itree),TString::Format("swdist%d",itree),10000,0,15);
1275 // results->Store(hxx,TString::Format("swdist%d",itree));
1276 // TH1 *hxy = new TH1F(TString::Format("bwdist%d",itree),TString::Format("bwdist%d",itree),10000,0,15);
1277 // results->Store(hxy,TString::Format("bwdist%d",itree));
1278 // for (Int_t iev=0; iev<fEventSample.size(); iev++) {
1279 // if (fEventSample[iev]->GetClass()!=0) hxy->Fill((fEventSample[iev])->GetWeight());
1280 // else hxx->Fill((fEventSample[iev])->GetWeight());
1281 // }
1282
1283 if(DoMulticlass()){
1284 if (fBoostType!="Grad"){
1285 Log() << kFATAL << "Multiclass is currently only supported by gradient boost. "
1286 << "Please change boost option accordingly (BoostType=Grad)." << Endl;
1287 }
1288
1289 UInt_t nClasses = DataInfo().GetNClasses();
1290 for (UInt_t i=0;i<nClasses;i++){
1291 // Careful: If fSepType is nullptr, the tree will be considered a regression tree and
1292 // use the correct output for gradboost (response rather than yesnoleaf) in checkEvent.
1293 // See TMVA::MethodBDT::InitGradBoost.
1294 fForest.push_back( new DecisionTree( fSepType, fMinNodeSize, fNCuts, &(DataInfo()), i,
1295 fRandomisedTrees, fUseNvars, fUsePoissonNvars, fMaxDepth,
1296 itree*nClasses+i, fNodePurityLimit, itree*nClasses+1));
1297 fForest.back()->SetNVars(GetNvar());
1298 if (fUseFisherCuts) {
1299 fForest.back()->SetUseFisherCuts();
1300 fForest.back()->SetMinLinCorrForFisher(fMinLinCorrForFisher);
1301 fForest.back()->SetUseExclusiveVars(fUseExclusiveVars);
1302 }
1303 // the minimum linear correlation between two variables demanded for use in fisher criterion in node splitting
1304
1305 nNodesBeforePruning = fForest.back()->BuildTree(*fTrainSample);
1306 Double_t bw = this->Boost(*fTrainSample, fForest.back(),i);
1307 if (bw > 0) {
1308 fBoostWeights.push_back(bw);
1309 }else{
1310 fBoostWeights.push_back(0);
1311 Log() << kWARNING << "stopped boosting at itree="<<itree << Endl;
1312 // fNTrees = itree+1; // that should stop the boosting
1314 }
1315 }
1316 }
1317 else{
1318
1319 DecisionTree* dt = new DecisionTree( fSepType, fMinNodeSize, fNCuts, &(DataInfo()), fSignalClass,
1320 fRandomisedTrees, fUseNvars, fUsePoissonNvars, fMaxDepth,
1321 itree, fNodePurityLimit, itree);
1322
1323 fForest.push_back(dt);
1324 fForest.back()->SetNVars(GetNvar());
1325 if (fUseFisherCuts) {
1326 fForest.back()->SetUseFisherCuts();
1327 fForest.back()->SetMinLinCorrForFisher(fMinLinCorrForFisher);
1328 fForest.back()->SetUseExclusiveVars(fUseExclusiveVars);
1329 }
1330
1331 nNodesBeforePruning = fForest.back()->BuildTree(*fTrainSample);
1332
1333 if (fUseYesNoLeaf && !DoRegression() && fBoostType!="Grad") { // remove leaf nodes where both daughter nodes are of same type
1334 nNodesBeforePruning = fForest.back()->CleanTree();
1335 }
1336
1339
1340 fForest.back()->SetPruneMethod(fPruneMethod); // set the pruning method for the tree
1341 fForest.back()->SetPruneStrength(fPruneStrength); // set the strength parameter
1342
1343 std::vector<const Event*> * validationSample = NULL;
1344 if(fAutomatic) validationSample = &fValidationSample;
1345 Double_t bw = this->Boost(*fTrainSample, fForest.back());
1346 if (bw > 0) {
1347 fBoostWeights.push_back(bw);
1348 }else{
1349 fBoostWeights.push_back(0);
1350 Log() << kWARNING << "stopped boosting at itree="<<itree << Endl;
1352 }
1353
1354 // if fAutomatic == true, pruneStrength will be the optimal pruning strength
1355 // determined by the pruning algorithm; otherwise, it is simply the strength parameter
1356 // set by the user
1357 if (fPruneMethod != DecisionTree::kNoPruning) fForest.back()->PruneTree(validationSample);
1358
1359 if (fUseYesNoLeaf && !DoRegression() && fBoostType!="Grad"){ // remove leaf nodes where both daughter nodes are of same type
1360 fForest.back()->CleanTree();
1361 }
1362 nNodesAfterPruning = fForest.back()->GetNNodes();
1365
1366 if (fInteractive){
1367 fInteractive->AddPoint(itree, fBoostWeight, fErrorFraction);
1368 }
1369 fITree = itree;
1370 fMonitorNtuple->Fill();
1371 if (fDoBoostMonitor){
1372 if (! DoRegression() ){
1373 if ( itree==fNTrees-1 || (!(itree%500)) ||
1374 (!(itree%250) && itree <1000)||
1375 (!(itree%100) && itree < 500)||
1376 (!(itree%50) && itree < 250)||
1377 (!(itree%25) && itree < 150)||
1378 (!(itree%10) && itree < 50)||
1379 (!(itree%5) && itree < 20)
1380 ) BoostMonitor(itree);
1381 }
1382 }
1383 }
1384 itree++;
1385 }
1386
1387 // get elapsed time
1388 Log() << kDEBUG << "\t<Train> elapsed time: " << timer.GetElapsedTime()
1389 << " " << Endl;
1390 if (fPruneMethod == DecisionTree::kNoPruning) {
1391 Log() << kDEBUG << "\t<Train> average number of nodes (w/o pruning) : "
1392 << nNodesBeforePruningCount/GetNTrees() << Endl;
1393 }
1394 else {
1395 Log() << kDEBUG << "\t<Train> average number of nodes before/after pruning : "
1396 << nNodesBeforePruningCount/GetNTrees() << " / "
1397 << nNodesAfterPruningCount/GetNTrees()
1398 << Endl;
1399 }
1401
1402
1403 // reset all previously stored/accumulated BOOST weights in the event sample
1404 // for (UInt_t iev=0; iev<fEventSample.size(); iev++) fEventSample[iev]->SetBoostWeight(1.);
1405 Log() << kDEBUG << "Now I delete the privat data sample"<< Endl;
1406 for (UInt_t i=0; i<fEventSample.size(); i++) delete fEventSample[i];
1407 for (UInt_t i=0; i<fValidationSample.size(); i++) delete fValidationSample[i];
1408 fEventSample.clear();
1409 fValidationSample.clear();
1410
1411 if (!fExitFromTraining) fIPyMaxIter = fIPyCurrentIter;
1412 ExitFromTraining();
1413}
1414
1415
1416////////////////////////////////////////////////////////////////////////////////
1417/// Returns MVA value: -1 for background, 1 for signal.
1418
1420{
1421 Double_t sum=0;
1422 for (UInt_t itree=0; itree<nTrees; itree++) {
1423 //loop over all trees in forest
1424 sum += fForest[itree]->CheckEvent(e,kFALSE);
1425
1426 }
1427 return 2.0/(1.0+exp(-2.0*sum))-1; //MVA output between -1 and 1
1428}
1429
1430////////////////////////////////////////////////////////////////////////////////
1431/// Calculate residual for all events.
1432
1433void TMVA::MethodBDT::UpdateTargets(std::vector<const TMVA::Event*>& eventSample, UInt_t cls)
1434{
1435 if (DoMulticlass()) {
1436 UInt_t nClasses = DataInfo().GetNClasses();
1437 Bool_t isLastClass = (cls == nClasses - 1);
1438
1439 #ifdef R__USE_IMT
1440 //
1441 // This is the multi-threaded multiclass version
1442 //
1443 // Note: we only need to update the predicted probabilities every
1444 // `nClasses` tree. Let's call a set of `nClasses` trees a "round". Thus
1445 // the algortihm is split in two parts `update_residuals` and
1446 // `update_residuals_last` where the latter is inteded to be run instead
1447 // of the former for the last tree in a "round".
1448 //
1449 std::map<const TMVA::Event *, std::vector<double>> & residuals = this->fResiduals;
1450 DecisionTree & lastTree = *(this->fForest.back());
1451
1452 auto update_residuals = [&residuals, &lastTree, cls](const TMVA::Event * e) {
1453 residuals[e].at(cls) += lastTree.CheckEvent(e, kFALSE);
1454 };
1455
1457 residuals[e].at(cls) += lastTree.CheckEvent(e, kFALSE);
1458
1460
1461 std::vector<Double_t> expCache(nClasses, 0.0);
1462 std::transform(residualsThisEvent.begin(),
1464 expCache.begin(), [](Double_t d) { return exp(d); });
1465
1466 Double_t exp_sum = std::accumulate(expCache.begin(),
1468 0.0);
1469
1470 for (UInt_t i = 0; i < nClasses; i++) {
1472
1473 Double_t res = (e->GetClass() == i) ? (1.0 - p_cls) : (-p_cls);
1474 const_cast<TMVA::Event *>(e)->SetTarget(i, res);
1475 }
1476 };
1477
1478 if (isLastClass) {
1479 TMVA::Config::Instance().GetThreadExecutor()
1481 } else {
1482 TMVA::Config::Instance().GetThreadExecutor()
1483 .Foreach(update_residuals, eventSample);
1484 }
1485 #else
1486 //
1487 // Single-threaded multiclass version
1488 //
1489 std::vector<Double_t> expCache;
1490 if (isLastClass) {
1491 expCache.resize(nClasses);
1492 }
1493
1494 for (auto e : eventSample) {
1495 fResiduals[e].at(cls) += fForest.back()->CheckEvent(e, kFALSE);
1496 if (isLastClass) {
1497 auto &residualsThisEvent = fResiduals[e];
1498 std::transform(residualsThisEvent.begin(),
1500 expCache.begin(), [](Double_t d) { return exp(d); });
1501
1502 Double_t exp_sum = std::accumulate(expCache.begin(),
1504 0.0);
1505
1506 for (UInt_t i = 0; i < nClasses; i++) {
1508
1509 Double_t res = (e->GetClass() == i) ? (1.0 - p_cls) : (-p_cls);
1510 const_cast<TMVA::Event *>(e)->SetTarget(i, res);
1511 }
1512 }
1513 }
1514 #endif
1515 } else {
1516 std::map<const TMVA::Event *, std::vector<double>> & residuals = this->fResiduals;
1517 DecisionTree & lastTree = *(this->fForest.back());
1518
1519 UInt_t signalClass = DataInfo().GetSignalClassIndex();
1520
1521 #ifdef R__USE_IMT
1523 double & residualAt0 = residuals[e].at(0);
1524 residualAt0 += lastTree.CheckEvent(e, kFALSE);
1525
1526 Double_t p_sig = 1.0 / (1.0 + exp(-2.0 * residualAt0));
1527 Double_t res = ((e->GetClass() == signalClass) ? (1.0 - p_sig) : (-p_sig));
1528
1529 const_cast<TMVA::Event *>(e)->SetTarget(0, res);
1530 };
1531
1532 TMVA::Config::Instance().GetThreadExecutor()
1533 .Foreach(update_residuals, eventSample);
1534 #else
1535 for (auto e : eventSample) {
1536 double & residualAt0 = residuals[e].at(0);
1537 residualAt0 += lastTree.CheckEvent(e, kFALSE);
1538
1539 Double_t p_sig = 1.0 / (1.0 + exp(-2.0 * residualAt0));
1540 Double_t res = ((e->GetClass() == signalClass) ? (1.0 - p_sig) : (-p_sig));
1541
1542 const_cast<TMVA::Event *>(e)->SetTarget(0, res);
1543 }
1544 #endif
1545 }
1546}
1547
1548////////////////////////////////////////////////////////////////////////////////
1549/// \brief Calculate residuals for all events and update targets for next iter.
1550///
1551/// \param[in] eventSample The collection of events currently under training.
1552/// \param[in] first Should be true when called before the first boosting
1553/// iteration has been run
1554///
1555void TMVA::MethodBDT::UpdateTargetsRegression(std::vector<const TMVA::Event*>& eventSample, Bool_t first)
1556{
1557 if (!first) {
1558#ifdef R__USE_IMT
1559 UInt_t nPartitions = TMVA::Config::Instance().GetThreadExecutor().GetPoolSize();
1561
1562 // need a lambda function to pass to TThreadExecutor::MapReduce
1563 auto f = [this, &nPartitions](UInt_t partition = 0) -> Int_t {
1564 Int_t start = 1.0 * partition / nPartitions * this->fEventSample.size();
1565 Int_t end = (partition + 1.0) / nPartitions * this->fEventSample.size();
1566
1567 for (Int_t i = start; i < end; ++i) {
1568 const TMVA::Event *e = fEventSample[i];
1569 LossFunctionEventInfo & lossInfo = fLossFunctionEventInfo.at(e);
1570 lossInfo.predictedValue += fForest.back()->CheckEvent(e, kFALSE);
1571 }
1572
1573 return 0;
1574 };
1575
1576 TMVA::Config::Instance().GetThreadExecutor().Map(f, seeds);
1577#else
1578 for (const TMVA::Event *e : fEventSample) {
1579 LossFunctionEventInfo & lossInfo = fLossFunctionEventInfo.at(e);
1580 lossInfo.predictedValue += fForest.back()->CheckEvent(e, kFALSE);
1581 }
1582#endif
1583 }
1584
1585 // NOTE: Set targets are also parallelised internally
1586 fRegressionLossFunctionBDTG->SetTargets(eventSample, fLossFunctionEventInfo);
1587
1588}
1589
1590////////////////////////////////////////////////////////////////////////////////
1591/// Calculate the desired response value for each region.
1592
1594{
1595 struct LeafInfo {
1597 Double_t sum2 = 0;
1598 };
1599
1600 std::unordered_map<TMVA::DecisionTreeNode*, LeafInfo> leaves;
1601 for (auto e : eventSample) {
1602 Double_t weight = e->GetWeight();
1603 TMVA::DecisionTreeNode* node = dt->GetEventNode(*e);
1604 auto &v = leaves[node];
1605 auto target = e->GetTarget(cls);
1606 v.sumWeightTarget += target * weight;
1607 v.sum2 += fabs(target) * (1.0 - fabs(target)) * weight;
1608 }
1609 for (auto &iLeave : leaves) {
1610 constexpr auto minValue = 1e-30;
1611 if (iLeave.second.sum2 < minValue) {
1612 iLeave.second.sum2 = minValue;
1613 }
1614 const Double_t K = DataInfo().GetNClasses();
1615 iLeave.first->SetResponse(fShrinkage * (K - 1) / K * iLeave.second.sumWeightTarget / iLeave.second.sum2);
1616 }
1617
1618 //call UpdateTargets before next tree is grown
1619
1620 DoMulticlass() ? UpdateTargets(fEventSample, cls) : UpdateTargets(fEventSample);
1621 return 1; //trees all have the same weight
1622}
1623
1624////////////////////////////////////////////////////////////////////////////////
1625/// Implementation of M_TreeBoost using any loss function as described by Friedman 1999.
1626
1628{
1629 // get the vector of events for each terminal so that we can calculate the constant fit value in each
1630 // terminal node
1631 // #### Not sure how many events are in each node in advance, so I can't parallelize this easily
1632 std::map<TMVA::DecisionTreeNode*,vector< TMVA::LossFunctionEventInfo > > leaves;
1633 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1634 TMVA::DecisionTreeNode* node = dt->GetEventNode(*(*e));
1635 (leaves[node]).push_back(fLossFunctionEventInfo[*e]);
1636 }
1637
1638 // calculate the constant fit for each terminal node based upon the events in the node
1639 // node (iLeave->first), vector of event information (iLeave->second)
1640 // #### could parallelize this and do the leaves at the same time, but this doesn't take very long compared
1641 // #### to the other processes
1643 iLeave!=leaves.end();++iLeave){
1644 Double_t fit = fRegressionLossFunctionBDTG->Fit(iLeave->second);
1645 (iLeave->first)->SetResponse(fShrinkage*fit);
1646 }
1647
1648 UpdateTargetsRegression(*fTrainSample);
1649
1650 return 1;
1651}
1652
1653////////////////////////////////////////////////////////////////////////////////
1654/// Initialize targets for first tree.
1655
1656void TMVA::MethodBDT::InitGradBoost( std::vector<const TMVA::Event*>& eventSample)
1657{
1658 // Should get rid of this line. It's just for debugging.
1659 //std::sort(eventSample.begin(), eventSample.end(), [](const TMVA::Event* a, const TMVA::Event* b){
1660 // return (a->GetTarget(0) < b->GetTarget(0)); });
1661 fSepType=NULL; //set fSepType to NULL (regression trees are used for both classification an regression)
1662 if(DoRegression()){
1663 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1664 fLossFunctionEventInfo[*e]= TMVA::LossFunctionEventInfo((*e)->GetTarget(0), 0, (*e)->GetWeight());
1665 }
1666
1667 fRegressionLossFunctionBDTG->Init(fLossFunctionEventInfo, fBoostWeights);
1668 UpdateTargetsRegression(*fTrainSample,kTRUE);
1669
1670 return;
1671 }
1672 else if(DoMulticlass()){
1673 UInt_t nClasses = DataInfo().GetNClasses();
1674 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1675 for (UInt_t i=0;i<nClasses;i++){
1676 //Calculate initial residua, assuming equal probability for all classes
1677 Double_t r = (*e)->GetClass()==i?(1-1.0/nClasses):(-1.0/nClasses);
1678 const_cast<TMVA::Event*>(*e)->SetTarget(i,r);
1679 fResiduals[*e].push_back(0);
1680 }
1681 }
1682 }
1683 else{
1684 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1685 Double_t r = (DataInfo().IsSignal(*e)?1:0)-0.5; //Calculate initial residua
1686 const_cast<TMVA::Event*>(*e)->SetTarget(0,r);
1687 fResiduals[*e].push_back(0);
1688 }
1689 }
1690
1691}
1692////////////////////////////////////////////////////////////////////////////////
1693/// Test the tree quality.. in terms of Misclassification.
1694
1696{
1698 for (UInt_t ievt=0; ievt<fValidationSample.size(); ievt++) {
1699 Bool_t isSignalType= (dt->CheckEvent(fValidationSample[ievt]) > fNodePurityLimit ) ? 1 : 0;
1700
1701 if (isSignalType == (DataInfo().IsSignal(fValidationSample[ievt])) ) {
1702 ncorrect += fValidationSample[ievt]->GetWeight();
1703 }
1704 else{
1705 nfalse += fValidationSample[ievt]->GetWeight();
1706 }
1707 }
1708
1709 return ncorrect / (ncorrect + nfalse);
1710}
1711
1712////////////////////////////////////////////////////////////////////////////////
1713/// Apply the boosting algorithm (the algorithm is selecte via the "option" given
1714/// in the constructor. The return value is the boosting weight.
1715
1717{
1719
1720 if (fBoostType=="AdaBoost") returnVal = this->AdaBoost (eventSample, dt);
1721 else if (fBoostType=="AdaCost") returnVal = this->AdaCost (eventSample, dt);
1722 else if (fBoostType=="Bagging") returnVal = this->Bagging ( );
1723 else if (fBoostType=="RegBoost") returnVal = this->RegBoost (eventSample, dt);
1724 else if (fBoostType=="AdaBoostR2") returnVal = this->AdaBoostR2(eventSample, dt);
1725 else if (fBoostType=="Grad"){
1726 if(DoRegression())
1727 returnVal = this->GradBoostRegression(eventSample, dt);
1728 else if(DoMulticlass())
1729 returnVal = this->GradBoost (eventSample, dt, cls);
1730 else
1731 returnVal = this->GradBoost (eventSample, dt);
1732 }
1733 else {
1734 Log() << kINFO << GetOptions() << Endl;
1735 Log() << kFATAL << "<Boost> unknown boost option " << fBoostType<< " called" << Endl;
1736 }
1737
1738 if (fBaggedBoost){
1739 GetBaggedSubSample(fEventSample);
1740 }
1741
1742
1743 return returnVal;
1744}
1745
1746////////////////////////////////////////////////////////////////////////////////
1747/// Fills the ROCIntegral vs Itree from the testSample for the monitoring plots
1748/// during the training .. but using the testing events
1749
1751{
1752 Results* results = Data()->GetResults(GetMethodName(),Types::kTraining, Types::kMaxAnalysisType);
1753
1754 TH1F *tmpS = new TH1F( "tmpS", "", 100 , -1., 1.00001 );
1755 TH1F *tmpB = new TH1F( "tmpB", "", 100 , -1., 1.00001 );
1756 TH1F *tmp;
1757
1758
1759 UInt_t signalClassNr = DataInfo().GetClassInfo("Signal")->GetNumber();
1760
1761 // const std::vector<Event*> events=Data()->GetEventCollection(Types::kTesting);
1762 // // fMethod->GetTransformationHandler().CalcTransformations(fMethod->Data()->GetEventCollection(Types::kTesting));
1763 // for (UInt_t iev=0; iev < events.size() ; iev++){
1764 // if (events[iev]->GetClass() == signalClassNr) tmp=tmpS;
1765 // else tmp=tmpB;
1766 // tmp->Fill(PrivateGetMvaValue(*(events[iev])),events[iev]->GetWeight());
1767 // }
1768
1769 UInt_t nevents = Data()->GetNTestEvents();
1770 for (UInt_t iev=0; iev < nevents; iev++){
1771 const Event* event = GetTestingEvent(iev);
1772
1773 if (event->GetClass() == signalClassNr) {tmp=tmpS;}
1774 else {tmp=tmpB;}
1775 tmp->Fill(PrivateGetMvaValue(event),event->GetWeight());
1776 }
1777 Double_t max=1;
1778
1779 std::vector<TH1F*> hS;
1780 std::vector<TH1F*> hB;
1781 for (UInt_t ivar=0; ivar<GetNvar(); ivar++){
1782 hS.push_back(new TH1F(TString::Format("SigVar%dAtTree%d",ivar,iTree).Data(),TString::Format("SigVar%dAtTree%d",ivar,iTree).Data(),100,DataInfo().GetVariableInfo(ivar).GetMin(),DataInfo().GetVariableInfo(ivar).GetMax()));
1783 hB.push_back(new TH1F(TString::Format("BkgVar%dAtTree%d",ivar,iTree).Data(),TString::Format("BkgVar%dAtTree%d",ivar,iTree).Data(),100,DataInfo().GetVariableInfo(ivar).GetMin(),DataInfo().GetVariableInfo(ivar).GetMax()));
1784 results->Store(hS.back(),hS.back()->GetTitle());
1785 results->Store(hB.back(),hB.back()->GetTitle());
1786 }
1787
1788
1789 for (UInt_t iev=0; iev < fEventSample.size(); iev++){
1790 if (fEventSample[iev]->GetBoostWeight() > max) max = 1.01*fEventSample[iev]->GetBoostWeight();
1791 }
1792 TH1F *tmpBoostWeightsS = new TH1F(TString::Format("BoostWeightsInTreeS%d",iTree).Data(),TString::Format("BoostWeightsInTreeS%d",iTree).Data(),100,0.,max);
1793 TH1F *tmpBoostWeightsB = new TH1F(TString::Format("BoostWeightsInTreeB%d",iTree).Data(),TString::Format("BoostWeightsInTreeB%d",iTree).Data(),100,0.,max);
1794 results->Store(tmpBoostWeightsS,tmpBoostWeightsS->GetTitle());
1795 results->Store(tmpBoostWeightsB,tmpBoostWeightsB->GetTitle());
1796
1798 std::vector<TH1F*> *h;
1799
1800 for (UInt_t iev=0; iev < fEventSample.size(); iev++){
1801 if (fEventSample[iev]->GetClass() == signalClassNr) {
1803 h=&hS;
1804 }else{
1806 h=&hB;
1807 }
1808 tmpBoostWeights->Fill(fEventSample[iev]->GetBoostWeight());
1809 for (UInt_t ivar=0; ivar<GetNvar(); ivar++){
1810 (*h)[ivar]->Fill(fEventSample[iev]->GetValue(ivar),fEventSample[iev]->GetWeight());
1811 }
1812 }
1813
1814
1815 TMVA::PDF *sig = new TMVA::PDF( " PDF Sig", tmpS, TMVA::PDF::kSpline3 );
1816 TMVA::PDF *bkg = new TMVA::PDF( " PDF Bkg", tmpB, TMVA::PDF::kSpline3 );
1817
1818
1819 TGraph* gr=results->GetGraph("BoostMonitorGraph");
1820 Int_t nPoints = gr->GetN();
1821 gr->Set(nPoints+1);
1822 gr->SetPoint(nPoints,(Double_t)iTree+1,GetROCIntegral(sig,bkg));
1823
1824 tmpS->Delete();
1825 tmpB->Delete();
1826
1827 delete sig;
1828 delete bkg;
1829
1830 return;
1831}
1832
1833////////////////////////////////////////////////////////////////////////////////
1834/// The AdaBoost implementation.
1835/// a new training sample is generated by weighting
1836/// events that are misclassified by the decision tree. The weight
1837/// applied is \f$ w = \frac{(1-err)}{err} \f$ or more general:
1838/// \f$ w = (\frac{(1-err)}{err})^\beta \f$
1839/// where \f$err\f$ is the fraction of misclassified events in the tree ( <0.5 assuming
1840/// demanding the that previous selection was better than random guessing)
1841/// and "beta" being a free parameter (standard: beta = 1) that modifies the
1842/// boosting.
1843
1844Double_t TMVA::MethodBDT::AdaBoost( std::vector<const TMVA::Event*>& eventSample, DecisionTree *dt )
1845{
1847
1848 std::vector<Double_t> sumw(DataInfo().GetNClasses(),0); //for individually re-scaling each class
1849
1850 Double_t maxDev=0;
1851 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1852 Double_t w = (*e)->GetWeight();
1853 sumGlobalw += w;
1854 UInt_t iclass=(*e)->GetClass();
1855 sumw[iclass] += w;
1856
1857 if ( DoRegression() ) {
1858 Double_t tmpDev = TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) );
1861 if (tmpDev > maxDev) maxDev = tmpDev;
1862 }else{
1863
1864 if (fUseYesNoLeaf){
1865 Bool_t isSignalType = (dt->CheckEvent(*e,fUseYesNoLeaf) > fNodePurityLimit );
1866 if (!(isSignalType == DataInfo().IsSignal(*e))) {
1868 }
1869 }else{
1870 Double_t dtoutput = (dt->CheckEvent(*e,fUseYesNoLeaf) - 0.5)*2.;
1872 if (DataInfo().IsSignal(*e)) trueType = 1;
1873 else trueType = -1;
1875 }
1876 }
1877 }
1878
1880 if ( DoRegression() ) {
1881 //if quadratic loss:
1882 if (fAdaBoostR2Loss=="linear"){
1884 }
1885 else if (fAdaBoostR2Loss=="quadratic"){
1887 }
1888 else if (fAdaBoostR2Loss=="exponential"){
1889 err = 0;
1890 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1891 Double_t w = (*e)->GetWeight();
1892 Double_t tmpDev = TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) );
1893 err += w * (1 - exp (-tmpDev/maxDev)) / sumGlobalw;
1894 }
1895
1896 }
1897 else {
1898 Log() << kFATAL << " you've chosen a Loss type for Adaboost other than linear, quadratic or exponential "
1899 << " namely " << fAdaBoostR2Loss << "\n"
1900 << "and this is not implemented... a typo in the options ??" <<Endl;
1901 }
1902 }
1903
1904 Log() << kDEBUG << "BDT AdaBoos wrong/all: " << sumGlobalwfalse << "/" << sumGlobalw << Endl;
1905
1906
1908 std::vector<Double_t> newSumw(sumw.size(),0);
1909
1911 if (err >= 0.5 && fUseYesNoLeaf) { // sanity check ... should never happen as otherwise there is apparently
1912 // something odd with the assignment of the leaf nodes (rem: you use the training
1913 // events for this determination of the error rate)
1914 if (dt->GetNNodes() == 1){
1915 Log() << kERROR << " YOUR tree has only 1 Node... kind of a funny *tree*. I cannot "
1916 << "boost such a thing... if after 1 step the error rate is == 0.5"
1917 << Endl
1918 << "please check why this happens, maybe too many events per node requested ?"
1919 << Endl;
1920
1921 }else{
1922 Log() << kERROR << " The error rate in the BDT boosting is > 0.5. ("<< err
1923 << ") That should not happen, please check your code (i.e... the BDT code), I "
1924 << " stop boosting here" << Endl;
1925 return -1;
1926 }
1927 err = 0.5;
1928 } else if (err < 0) {
1929 Log() << kERROR << " The error rate in the BDT boosting is < 0. That can happen"
1930 << " due to improper treatment of negative weights in a Monte Carlo.. (if you have"
1931 << " an idea on how to do it in a better way, please let me know (Helge.Voss@cern.ch)"
1932 << " for the time being I set it to its absolute value.. just to continue.." << Endl;
1933 err = TMath::Abs(err);
1934 }
1935 if (fUseYesNoLeaf)
1936 boostWeight = TMath::Log((1.-err)/err)*fAdaBoostBeta;
1937 else
1938 boostWeight = TMath::Log((1.+err)/(1-err))*fAdaBoostBeta;
1939
1940
1941 Log() << kDEBUG << "BDT AdaBoos wrong/all: " << sumGlobalwfalse << "/" << sumGlobalw << " 1-err/err="<<boostWeight<< " log.."<<TMath::Log(boostWeight)<<Endl;
1942
1943 Results* results = Data()->GetResults(GetMethodName(),Types::kTraining, Types::kMaxAnalysisType);
1944
1945
1946 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1947
1948 if (fUseYesNoLeaf||DoRegression()){
1949 if ((!( (dt->CheckEvent(*e,fUseYesNoLeaf) > fNodePurityLimit ) == DataInfo().IsSignal(*e))) || DoRegression()) {
1951
1952 if (DoRegression()) boostfactor = TMath::Power(1/boostWeight,(1.-TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) )/maxDev ) );
1953 if ( (*e)->GetWeight() > 0 ){
1954 (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
1955 // Helge change back (*e)->ScaleBoostWeight(boostfactor);
1956 if (DoRegression()) results->GetHist("BoostWeights")->Fill(boostfactor);
1957 } else {
1958 if ( fInverseBoostNegWeights )(*e)->ScaleBoostWeight( 1. / boostfactor); // if the original event weight is negative, and you want to "increase" the events "positive" influence, you'd rather make the event weight "smaller" in terms of it's absolute value while still keeping it something "negative"
1959 else (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
1960
1961 }
1962 }
1963
1964 }else{
1965 Double_t dtoutput = (dt->CheckEvent(*e,fUseYesNoLeaf) - 0.5)*2.;
1967 if (DataInfo().IsSignal(*e)) trueType = 1;
1968 else trueType = -1;
1970
1971 if ( (*e)->GetWeight() > 0 ){
1972 (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
1973 // Helge change back (*e)->ScaleBoostWeight(boostfactor);
1974 if (DoRegression()) results->GetHist("BoostWeights")->Fill(boostfactor);
1975 } else {
1976 if ( fInverseBoostNegWeights )(*e)->ScaleBoostWeight( 1. / boostfactor); // if the original event weight is negative, and you want to "increase" the events "positive" influence, you'd rather make the event weight "smaller" in terms of it's absolute value while still keeping it something "negative"
1977 else (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
1978 }
1979 }
1980 newSumGlobalw+=(*e)->GetWeight();
1981 newSumw[(*e)->GetClass()] += (*e)->GetWeight();
1982 }
1983
1984
1985 // Double_t globalNormWeight=sumGlobalw/newSumGlobalw;
1987 Log() << kDEBUG << "new Nsig="<<newSumw[0]*globalNormWeight << " new Nbkg="<<newSumw[1]*globalNormWeight << Endl;
1988
1989
1990 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
1991 // if (fRenormByClass) (*e)->ScaleBoostWeight( normWeightByClass[(*e)->GetClass()] );
1992 // else (*e)->ScaleBoostWeight( globalNormWeight );
1993 // else (*e)->ScaleBoostWeight( globalNormWeight );
1994 if (DataInfo().IsSignal(*e))(*e)->ScaleBoostWeight( globalNormWeight * fSigToBkgFraction );
1995 else (*e)->ScaleBoostWeight( globalNormWeight );
1996 }
1997
1998 if (!(DoRegression()))results->GetHist("BoostWeights")->Fill(boostWeight);
1999 results->GetHist("BoostWeightsVsTree")->SetBinContent(fForest.size(),boostWeight);
2000 results->GetHist("ErrorFrac")->SetBinContent(fForest.size(),err);
2001
2002 fBoostWeight = boostWeight;
2003 fErrorFraction = err;
2004
2005 return boostWeight;
2006}
2007
2008////////////////////////////////////////////////////////////////////////////////
2009/// The AdaCost boosting algorithm takes a simple cost Matrix (currently fixed for
2010/// all events... later could be modified to use individual cost matrices for each
2011/// events as in the original paper...
2012///
2013/// true_signal true_bkg
2014/// ----------------------------------
2015/// sel_signal | Css Ctb_ss Cxx.. in the range [0,1]
2016/// sel_bkg | Cts_sb Cbb
2017///
2018/// and takes this into account when calculating the mis class. cost (former: error fraction):
2019///
2020/// err = sum_events ( weight* y_true*y_sel * beta(event)
2021
2023{
2024 Double_t Css = fCss;
2025 Double_t Cbb = fCbb;
2026 Double_t Cts_sb = fCts_sb;
2027 Double_t Ctb_ss = fCtb_ss;
2028
2030
2031 std::vector<Double_t> sumw(DataInfo().GetNClasses(),0); //for individually re-scaling each class
2032
2033 for (vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2034 Double_t w = (*e)->GetWeight();
2036 UInt_t iclass=(*e)->GetClass();
2037
2038 sumw[iclass] += w;
2039
2040 if ( DoRegression() ) {
2041 Log() << kFATAL << " AdaCost not implemented for regression"<<Endl;
2042 }else{
2043
2044 Double_t dtoutput = (dt->CheckEvent(*e,false) - 0.5)*2.;
2046 Bool_t isTrueSignal = DataInfo().IsSignal(*e);
2048 if (isTrueSignal) trueType = 1;
2049 else trueType = -1;
2050
2051 Double_t cost=0;
2052 if (isTrueSignal && isSelectedSignal) cost=Css;
2053 else if (isTrueSignal && !isSelectedSignal) cost=Cts_sb;
2054 else if (!isTrueSignal && isSelectedSignal) cost=Ctb_ss;
2055 else if (!isTrueSignal && !isSelectedSignal) cost=Cbb;
2056 else Log() << kERROR << "something went wrong in AdaCost" << Endl;
2057
2059
2060 }
2061 }
2062
2063 if ( DoRegression() ) {
2064 Log() << kFATAL << " AdaCost not implemented for regression"<<Endl;
2065 }
2066
2067 // Log() << kDEBUG << "BDT AdaBoos wrong/all: " << sumGlobalCost << "/" << sumGlobalWeights << Endl;
2068 // Log() << kWARNING << "BDT AdaBoos wrong/all: " << sumGlobalCost << "/" << sumGlobalWeights << Endl;
2070 // Log() << kWARNING << "BDT AdaBoos wrong/all: " << sumGlobalCost << "/" << sumGlobalWeights << Endl;
2071
2072
2075
2076 Double_t boostWeight = TMath::Log((1+sumGlobalCost)/(1-sumGlobalCost)) * fAdaBoostBeta;
2077
2078 Results* results = Data()->GetResults(GetMethodName(),Types::kTraining, Types::kMaxAnalysisType);
2079
2080 for (vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2081 Double_t dtoutput = (dt->CheckEvent(*e,false) - 0.5)*2.;
2083 Bool_t isTrueSignal = DataInfo().IsSignal(*e);
2085 if (isTrueSignal) trueType = 1;
2086 else trueType = -1;
2087
2088 Double_t cost=0;
2089 if (isTrueSignal && isSelectedSignal) cost=Css;
2090 else if (isTrueSignal && !isSelectedSignal) cost=Cts_sb;
2091 else if (!isTrueSignal && isSelectedSignal) cost=Ctb_ss;
2092 else if (!isTrueSignal && !isSelectedSignal) cost=Cbb;
2093 else Log() << kERROR << "something went wrong in AdaCost" << Endl;
2094
2096 if (DoRegression())Log() << kFATAL << " AdaCost not implemented for regression"<<Endl;
2097 if ( (*e)->GetWeight() > 0 ){
2098 (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
2099 // Helge change back (*e)->ScaleBoostWeight(boostfactor);
2100 if (DoRegression())Log() << kFATAL << " AdaCost not implemented for regression"<<Endl;
2101 } else {
2102 if ( fInverseBoostNegWeights )(*e)->ScaleBoostWeight( 1. / boostfactor); // if the original event weight is negative, and you want to "increase" the events "positive" influence, you'd rather make the event weight "smaller" in terms of it's absolute value while still keeping it something "negative"
2103 }
2104
2105 newSumGlobalWeights+=(*e)->GetWeight();
2106 newSumClassWeights[(*e)->GetClass()] += (*e)->GetWeight();
2107 }
2108
2109
2110 // Double_t globalNormWeight=sumGlobalWeights/newSumGlobalWeights;
2112 Log() << kDEBUG << "new Nsig="<<newSumClassWeights[0]*globalNormWeight << " new Nbkg="<<newSumClassWeights[1]*globalNormWeight << Endl;
2113
2114
2115 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2116 // if (fRenormByClass) (*e)->ScaleBoostWeight( normWeightByClass[(*e)->GetClass()] );
2117 // else (*e)->ScaleBoostWeight( globalNormWeight );
2118 if (DataInfo().IsSignal(*e))(*e)->ScaleBoostWeight( globalNormWeight * fSigToBkgFraction );
2119 else (*e)->ScaleBoostWeight( globalNormWeight );
2120 }
2121
2122
2123 if (!(DoRegression()))results->GetHist("BoostWeights")->Fill(boostWeight);
2124 results->GetHist("BoostWeightsVsTree")->SetBinContent(fForest.size(),boostWeight);
2125 results->GetHist("ErrorFrac")->SetBinContent(fForest.size(),err);
2126
2127 fBoostWeight = boostWeight;
2128 fErrorFraction = err;
2129
2130
2131 return boostWeight;
2132}
2133
2134////////////////////////////////////////////////////////////////////////////////
2135/// Call it boot-strapping, re-sampling or whatever you like, in the end it is nothing
2136/// else but applying "random" poisson weights to each event.
2137
2139{
2140 // this is now done in "MethodBDT::Boost as it might be used by other boost methods, too
2141 // GetBaggedSample(eventSample);
2142
2143 return 1.; //here as there are random weights for each event, just return a constant==1;
2144}
2145
2146////////////////////////////////////////////////////////////////////////////////
2147/// Fills fEventSample with fBaggedSampleFraction*NEvents random training events.
2148
2149void TMVA::MethodBDT::GetBaggedSubSample(std::vector<const TMVA::Event*>& eventSample)
2150{
2151
2152 Double_t n;
2153 TRandom3 *trandom = new TRandom3(100*fForest.size()+1234);
2154
2155 if (!fSubSample.empty()) fSubSample.clear();
2156
2157 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2158 n = trandom->PoissonD(fBaggedSampleFraction);
2159 for (Int_t i=0;i<n;i++) fSubSample.push_back(*e);
2160 }
2161
2162 delete trandom;
2163 return;
2164
2165 /*
2166 UInt_t nevents = fEventSample.size();
2167
2168 if (!fSubSample.empty()) fSubSample.clear();
2169 TRandom3 *trandom = new TRandom3(fForest.size()+1);
2170
2171 for (UInt_t ievt=0; ievt<nevents; ievt++) { // recreate new random subsample
2172 if(trandom->Rndm()<fBaggedSampleFraction)
2173 fSubSample.push_back(fEventSample[ievt]);
2174 }
2175 delete trandom;
2176 */
2177
2178}
2179
2180////////////////////////////////////////////////////////////////////////////////
2181/// A special boosting only for Regression (not implemented).
2182
2183Double_t TMVA::MethodBDT::RegBoost( std::vector<const TMVA::Event*>& /* eventSample */, DecisionTree* /* dt */ )
2184{
2185 return 1;
2186}
2187
2188////////////////////////////////////////////////////////////////////////////////
2189/// Adaption of the AdaBoost to regression problems (see H.Drucker 1997).
2190
2192{
2193 if ( !DoRegression() ) Log() << kFATAL << "Somehow you chose a regression boost method for a classification job" << Endl;
2194
2195 Double_t err=0, sumw=0, sumwfalse=0, sumwfalse2=0;
2196 Double_t maxDev=0;
2197 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2198 Double_t w = (*e)->GetWeight();
2199 sumw += w;
2200
2201 Double_t tmpDev = TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) );
2202 sumwfalse += w * tmpDev;
2204 if (tmpDev > maxDev) maxDev = tmpDev;
2205 }
2206
2207 //if quadratic loss:
2208 if (fAdaBoostR2Loss=="linear"){
2209 err = sumwfalse/maxDev/sumw ;
2210 }
2211 else if (fAdaBoostR2Loss=="quadratic"){
2213 }
2214 else if (fAdaBoostR2Loss=="exponential"){
2215 err = 0;
2216 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2217 Double_t w = (*e)->GetWeight();
2218 Double_t tmpDev = TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) );
2219 err += w * (1 - exp (-tmpDev/maxDev)) / sumw;
2220 }
2221
2222 }
2223 else {
2224 Log() << kFATAL << " you've chosen a Loss type for Adaboost other than linear, quadratic or exponential "
2225 << " namely " << fAdaBoostR2Loss << "\n"
2226 << "and this is not implemented... a typo in the options ??" <<Endl;
2227 }
2228
2229
2230 if (err >= 0.5) { // sanity check ... should never happen as otherwise there is apparently
2231 // something odd with the assignment of the leaf nodes (rem: you use the training
2232 // events for this determination of the error rate)
2233 if (dt->GetNNodes() == 1){
2234 Log() << kERROR << " YOUR tree has only 1 Node... kind of a funny *tree*. I cannot "
2235 << "boost such a thing... if after 1 step the error rate is == 0.5"
2236 << Endl
2237 << "please check why this happens, maybe too many events per node requested ?"
2238 << Endl;
2239
2240 }else{
2241 Log() << kERROR << " The error rate in the BDT boosting is > 0.5. ("<< err
2242 << ") That should not happen, but is possible for regression trees, and"
2243 << " should trigger a stop for the boosting. please check your code (i.e... the BDT code), I "
2244 << " stop boosting " << Endl;
2245 return -1;
2246 }
2247 err = 0.5;
2248 } else if (err < 0) {
2249 Log() << kERROR << " The error rate in the BDT boosting is < 0. That can happen"
2250 << " due to improper treatment of negative weights in a Monte Carlo.. (if you have"
2251 << " an idea on how to do it in a better way, please let me know (Helge.Voss@cern.ch)"
2252 << " for the time being I set it to its absolute value.. just to continue.." << Endl;
2253 err = TMath::Abs(err);
2254 }
2255
2256 Double_t boostWeight = err / (1.-err);
2257 Double_t newSumw=0;
2258
2259 Results* results = Data()->GetResults(GetMethodName(), Types::kTraining, Types::kMaxAnalysisType);
2260
2261 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2262 Double_t boostfactor = TMath::Power(boostWeight,(1.-TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) )/maxDev ) );
2263 results->GetHist("BoostWeights")->Fill(boostfactor);
2264 // std::cout << "R2 " << boostfactor << " " << boostWeight << " " << (1.-TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) )/maxDev) << std::endl;
2265 if ( (*e)->GetWeight() > 0 ){
2266 Float_t newBoostWeight = (*e)->GetBoostWeight() * boostfactor;
2267 Float_t newWeight = (*e)->GetWeight() * (*e)->GetBoostWeight() * boostfactor;
2268 if (newWeight == 0) {
2269 Log() << kINFO << "Weight= " << (*e)->GetWeight() << Endl;
2270 Log() << kINFO << "BoostWeight= " << (*e)->GetBoostWeight() << Endl;
2271 Log() << kINFO << "boostweight="<<boostWeight << " err= " <<err << Endl;
2272 Log() << kINFO << "NewBoostWeight= " << newBoostWeight << Endl;
2273 Log() << kINFO << "boostfactor= " << boostfactor << Endl;
2274 Log() << kINFO << "maxDev = " << maxDev << Endl;
2275 Log() << kINFO << "tmpDev = " << TMath::Abs(dt->CheckEvent(*e,kFALSE) - (*e)->GetTarget(0) ) << Endl;
2276 Log() << kINFO << "target = " << (*e)->GetTarget(0) << Endl;
2277 Log() << kINFO << "estimate = " << dt->CheckEvent(*e,kFALSE) << Endl;
2278 }
2279 (*e)->SetBoostWeight( newBoostWeight );
2280 // (*e)->SetBoostWeight( (*e)->GetBoostWeight() * boostfactor);
2281 } else {
2282 (*e)->SetBoostWeight( (*e)->GetBoostWeight() / boostfactor);
2283 }
2284 newSumw+=(*e)->GetWeight();
2285 }
2286
2287 // re-normalise the weights
2289 for (std::vector<const TMVA::Event*>::const_iterator e=eventSample.begin(); e!=eventSample.end();++e) {
2290 //Helge (*e)->ScaleBoostWeight( sumw/newSumw);
2291 // (*e)->ScaleBoostWeight( normWeight);
2292 (*e)->SetBoostWeight( (*e)->GetBoostWeight() * normWeight );
2293 }
2294
2295
2296 results->GetHist("BoostWeightsVsTree")->SetBinContent(fForest.size(),1./boostWeight);
2297 results->GetHist("ErrorFrac")->SetBinContent(fForest.size(),err);
2298
2299 fBoostWeight = boostWeight;
2300 fErrorFraction = err;
2301
2302 return TMath::Log(1./boostWeight);
2303}
2304
2305////////////////////////////////////////////////////////////////////////////////
2306/// Write weights to XML.
2307
2308void TMVA::MethodBDT::AddWeightsXMLTo( void* parent ) const
2309{
2310 void* wght = gTools().AddChild(parent, "Weights");
2311
2312 if (fDoPreselection){
2313 for (UInt_t ivar=0; ivar<GetNvar(); ivar++){
2314 gTools().AddAttr( wght, TString::Format("PreselectionLowBkgVar%d",ivar).Data(), fIsLowBkgCut[ivar]);
2315 gTools().AddAttr( wght, TString::Format("PreselectionLowBkgVar%dValue",ivar).Data(), fLowBkgCut[ivar]);
2316 gTools().AddAttr( wght, TString::Format("PreselectionLowSigVar%d",ivar).Data(), fIsLowSigCut[ivar]);
2317 gTools().AddAttr( wght, TString::Format("PreselectionLowSigVar%dValue",ivar).Data(), fLowSigCut[ivar]);
2318 gTools().AddAttr( wght, TString::Format("PreselectionHighBkgVar%d",ivar).Data(), fIsHighBkgCut[ivar]);
2319 gTools().AddAttr( wght, TString::Format("PreselectionHighBkgVar%dValue",ivar).Data(),fHighBkgCut[ivar]);
2320 gTools().AddAttr( wght, TString::Format("PreselectionHighSigVar%d",ivar).Data(), fIsHighSigCut[ivar]);
2321 gTools().AddAttr( wght, TString::Format("PreselectionHighSigVar%dValue",ivar).Data(),fHighSigCut[ivar]);
2322 }
2323 }
2324
2325
2326 gTools().AddAttr( wght, "NTrees", fForest.size() );
2327 gTools().AddAttr( wght, "AnalysisType", fForest.back()->GetAnalysisType() );
2328
2329 for (UInt_t i=0; i< fForest.size(); i++) {
2330 void* trxml = fForest[i]->AddXMLTo(wght);
2331 gTools().AddAttr( trxml, "boostWeight", fBoostWeights[i] );
2332 gTools().AddAttr( trxml, "itree", i );
2333 }
2334}
2335
2336////////////////////////////////////////////////////////////////////////////////
2337/// Reads the BDT from the xml file.
2338
2340 UInt_t i;
2341 for (i=0; i<fForest.size(); i++) delete fForest[i];
2342 fForest.clear();
2343 fBoostWeights.clear();
2344
2345 UInt_t ntrees;
2346 UInt_t analysisType;
2348
2349
2350 if (gTools().HasAttr( parent, TString::Format("PreselectionLowBkgVar%d",0).Data())) {
2351 fIsLowBkgCut.resize(GetNvar());
2352 fLowBkgCut.resize(GetNvar());
2353 fIsLowSigCut.resize(GetNvar());
2354 fLowSigCut.resize(GetNvar());
2355 fIsHighBkgCut.resize(GetNvar());
2356 fHighBkgCut.resize(GetNvar());
2357 fIsHighSigCut.resize(GetNvar());
2358 fHighSigCut.resize(GetNvar());
2359
2362 for (UInt_t ivar=0; ivar<GetNvar(); ivar++){
2363 gTools().ReadAttr( parent, TString::Format("PreselectionLowBkgVar%d",ivar).Data(), tmpBool);
2364 fIsLowBkgCut[ivar]=tmpBool;
2365 gTools().ReadAttr( parent, TString::Format("PreselectionLowBkgVar%dValue",ivar).Data(), tmpDouble);
2366 fLowBkgCut[ivar]=tmpDouble;
2367 gTools().ReadAttr( parent, TString::Format("PreselectionLowSigVar%d",ivar).Data(), tmpBool);
2368 fIsLowSigCut[ivar]=tmpBool;
2369 gTools().ReadAttr( parent, TString::Format("PreselectionLowSigVar%dValue",ivar).Data(), tmpDouble);
2370 fLowSigCut[ivar]=tmpDouble;
2371 gTools().ReadAttr( parent, TString::Format("PreselectionHighBkgVar%d",ivar).Data(), tmpBool);
2372 fIsHighBkgCut[ivar]=tmpBool;
2373 gTools().ReadAttr( parent, TString::Format("PreselectionHighBkgVar%dValue",ivar).Data(), tmpDouble);
2374 fHighBkgCut[ivar]=tmpDouble;
2375 gTools().ReadAttr( parent, TString::Format("PreselectionHighSigVar%d",ivar).Data(),tmpBool);
2376 fIsHighSigCut[ivar]=tmpBool;
2377 gTools().ReadAttr( parent, TString::Format("PreselectionHighSigVar%dValue",ivar).Data(), tmpDouble);
2378 fHighSigCut[ivar]=tmpDouble;
2379 }
2380 }
2381
2382
2383 gTools().ReadAttr( parent, "NTrees", ntrees );
2384
2385 if(gTools().HasAttr(parent, "TreeType")) { // pre 4.1.0 version
2386 gTools().ReadAttr( parent, "TreeType", analysisType );
2387 } else { // from 4.1.0 onwards
2388 gTools().ReadAttr( parent, "AnalysisType", analysisType );
2389 }
2390
2391 void* ch = gTools().GetChild(parent);
2392 i=0;
2393 while(ch) {
2394 fForest.push_back( dynamic_cast<DecisionTree*>( DecisionTree::CreateFromXML(ch, GetTrainingTMVAVersionCode()) ) );
2395 fForest.back()->SetAnalysisType(Types::EAnalysisType(analysisType));
2396 fForest.back()->SetTreeID(i++);
2397 gTools().ReadAttr(ch,"boostWeight",boostWeight);
2398 fBoostWeights.push_back(boostWeight);
2399 ch = gTools().GetNextChild(ch);
2400 }
2401}
2402
2403////////////////////////////////////////////////////////////////////////////////
2404/// Read the weights (BDT coefficients).
2405
2407{
2408 TString dummy;
2409 // Types::EAnalysisType analysisType;
2410 Int_t analysisType(0);
2411
2412 // coverity[tainted_data_argument]
2413 istr >> dummy >> fNTrees;
2414 Log() << kINFO << "Read " << fNTrees << " Decision trees" << Endl;
2415
2416 for (UInt_t i=0;i<fForest.size();i++) delete fForest[i];
2417 fForest.clear();
2418 fBoostWeights.clear();
2419 Int_t iTree;
2421 for (int i=0;i<fNTrees;i++) {
2422 istr >> dummy >> iTree >> dummy >> boostWeight;
2423 if (iTree != i) {
2424 fForest.back()->Print( std::cout );
2425 Log() << kFATAL << "Error while reading weight file; mismatch iTree="
2426 << iTree << " i=" << i
2427 << " dummy " << dummy
2428 << " boostweight " << boostWeight
2429 << Endl;
2430 }
2431 fForest.push_back( new DecisionTree() );
2432 fForest.back()->SetAnalysisType(Types::EAnalysisType(analysisType));
2433 fForest.back()->SetTreeID(i);
2434 fForest.back()->Read(istr, GetTrainingTMVAVersionCode());
2435 fBoostWeights.push_back(boostWeight);
2436 }
2437}
2438
2439////////////////////////////////////////////////////////////////////////////////
2440
2442 return this->GetMvaValue( err, errUpper, 0 );
2443}
2444
2445////////////////////////////////////////////////////////////////////////////////
2446/// Return the MVA value (range [-1;1]) that classifies the
2447/// event according to the majority vote from the total number of
2448/// decision trees.
2449
2451{
2452 const Event* ev = GetEvent();
2453 if (fDoPreselection) {
2454 Double_t val = ApplyPreselectionCuts(ev);
2455 if (TMath::Abs(val)>0.05) return val;
2456 }
2457 return PrivateGetMvaValue(ev, err, errUpper, useNTrees);
2458
2459}
2460
2461////////////////////////////////////////////////////////////////////////////////
2462/// Return the MVA value (range [-1;1]) that classifies the
2463/// event according to the majority vote from the total number of
2464/// decision trees.
2465
2467{
2468 // cannot determine error
2469 NoErrorCalc(err, errUpper);
2470
2471 // allow for the possibility to use less trees in the actual MVA calculation
2472 // than have been originally trained.
2473 UInt_t nTrees = fForest.size();
2474
2475 if (useNTrees > 0 ) nTrees = useNTrees;
2476
2477 if (fBoostType=="Grad") return GetGradBoostMVA(ev,nTrees);
2478
2479 Double_t myMVA = 0;
2480 Double_t norm = 0;
2481 for (UInt_t itree=0; itree<nTrees; itree++) {
2482 //
2483 myMVA += fBoostWeights[itree] * fForest[itree]->CheckEvent(ev,fUseYesNoLeaf);
2484 norm += fBoostWeights[itree];
2485 }
2486 return ( norm > std::numeric_limits<double>::epsilon() ) ? myMVA /= norm : 0 ;
2487}
2488
2489
2490////////////////////////////////////////////////////////////////////////////////
2491/// Get the multiclass MVA response for the BDT classifier.
2492
2493const std::vector<Float_t>& TMVA::MethodBDT::GetMulticlassValues()
2494{
2495 const TMVA::Event *e = GetEvent();
2496 if (fMulticlassReturnVal == NULL) fMulticlassReturnVal = new std::vector<Float_t>();
2497 fMulticlassReturnVal->clear();
2498
2499 UInt_t nClasses = DataInfo().GetNClasses();
2500 std::vector<Double_t> temp(nClasses);
2501 auto forestSize = fForest.size();
2502
2503 #ifdef R__USE_IMT
2504 std::vector<TMVA::DecisionTree *> forest = fForest;
2505 auto get_output = [&e, &forest, &temp, forestSize, nClasses](UInt_t iClass) {
2507 temp[iClass] += forest[itree]->CheckEvent(e, kFALSE);
2508 }
2509 };
2510
2511 TMVA::Config::Instance().GetThreadExecutor()
2512 .Foreach(get_output, ROOT::TSeqU(nClasses));
2513 #else
2514 // trees 0, nClasses, 2*nClasses, ... belong to class 0
2515 // trees 1, nClasses+1, 2*nClasses+1, ... belong to class 1 and so forth
2516 UInt_t classOfTree = 0;
2517 for (UInt_t itree = 0; itree < forestSize; ++itree) {
2518 temp[classOfTree] += fForest[itree]->CheckEvent(e, kFALSE);
2519 if (++classOfTree == nClasses) classOfTree = 0; // cheap modulo
2520 }
2521 #endif
2522
2523 // we want to calculate sum of exp(temp[j] - temp[i]) for all i,j (i!=j)
2524 // first calculate exp(), then replace minus with division.
2525 std::transform(temp.begin(), temp.end(), temp.begin(), [](Double_t d){return exp(d);});
2526
2527 Double_t exp_sum = std::accumulate(temp.begin(), temp.end(), 0.0);
2528
2529 for (UInt_t i = 0; i < nClasses; i++) {
2530 Double_t p_cls = temp[i] / exp_sum;
2531 (*fMulticlassReturnVal).push_back(p_cls);
2532 }
2533
2534 return *fMulticlassReturnVal;
2535}
2536
2537////////////////////////////////////////////////////////////////////////////////
2538/// Get the regression value generated by the BDTs.
2539
2540const std::vector<Float_t> & TMVA::MethodBDT::GetRegressionValues()
2541{
2542
2543 if (fRegressionReturnVal == NULL) fRegressionReturnVal = new std::vector<Float_t>();
2544 fRegressionReturnVal->clear();
2545
2546 const Event * ev = GetEvent();
2547 Event * evT = new Event(*ev);
2548
2549 Double_t myMVA = 0;
2550 Double_t norm = 0;
2551 if (fBoostType=="AdaBoostR2") {
2552 // rather than using the weighted average of the tree respones in the forest
2553 // H.Decker(1997) proposed to use the "weighted median"
2554
2555 // sort all individual tree responses according to the prediction value
2556 // (keep the association to their tree weight)
2557 // the sum up all the associated weights (starting from the one whose tree
2558 // yielded the smalles response) up to the tree "t" at which you've
2559 // added enough tree weights to have more than half of the sum of all tree weights.
2560 // choose as response of the forest that one which belongs to this "t"
2561
2562 vector< Double_t > response(fForest.size());
2563 vector< Double_t > weight(fForest.size());
2565
2566 for (UInt_t itree=0; itree<fForest.size(); itree++) {
2567 response[itree] = fForest[itree]->CheckEvent(ev,kFALSE);
2568 weight[itree] = fBoostWeights[itree];
2569 totalSumOfWeights += fBoostWeights[itree];
2570 }
2571
2572 std::vector< std::vector<Double_t> > vtemp;
2573 vtemp.push_back( response ); // this is the vector that will get sorted
2574 vtemp.push_back( weight );
2576
2577 Int_t t=0;
2579 while (sumOfWeights <= totalSumOfWeights/2.) {
2580 sumOfWeights += vtemp[1][t];
2581 t++;
2582 }
2583
2584 Double_t rVal=0;
2585 Int_t count=0;
2586 for (UInt_t i= TMath::Max(UInt_t(0),UInt_t(t-(fForest.size()/6)-0.5));
2587 i< TMath::Min(UInt_t(fForest.size()),UInt_t(t+(fForest.size()/6)+0.5)); i++) {
2588 count++;
2589 rVal+=vtemp[0][i];
2590 }
2591 // fRegressionReturnVal->push_back( rVal/Double_t(count));
2592 evT->SetTarget(0, rVal/Double_t(count) );
2593 }
2594 else if(fBoostType=="Grad"){
2595 for (UInt_t itree=0; itree<fForest.size(); itree++) {
2596 myMVA += fForest[itree]->CheckEvent(ev,kFALSE);
2597 }
2598 // fRegressionReturnVal->push_back( myMVA+fBoostWeights[0]);
2599 evT->SetTarget(0, myMVA+fBoostWeights[0] );
2600 }
2601 else{
2602 for (UInt_t itree=0; itree<fForest.size(); itree++) {
2603 //
2604 myMVA += fBoostWeights[itree] * fForest[itree]->CheckEvent(ev,kFALSE);
2605 norm += fBoostWeights[itree];
2606 }
2607 // fRegressionReturnVal->push_back( ( norm > std::numeric_limits<double>::epsilon() ) ? myMVA /= norm : 0 );
2608 evT->SetTarget(0, ( norm > std::numeric_limits<double>::epsilon() ) ? myMVA /= norm : 0 );
2609 }
2610
2611
2612
2613 const Event* evT2 = GetTransformationHandler().InverseTransform( evT );
2614 fRegressionReturnVal->push_back( evT2->GetTarget(0) );
2615
2616 delete evT;
2617
2618
2619 return *fRegressionReturnVal;
2620}
2621
2622////////////////////////////////////////////////////////////////////////////////
2623/// Here we could write some histograms created during the processing
2624/// to the output file.
2625
2627{
2628 Log() << kDEBUG << "\tWrite monitoring histograms to file: " << BaseDir()->GetPath() << Endl;
2629
2630 //Results* results = Data()->GetResults(GetMethodName(), Types::kTraining, Types::kMaxAnalysisType);
2631 //results->GetStorage()->Write();
2632 fMonitorNtuple->Write();
2633}
2634
2635////////////////////////////////////////////////////////////////////////////////
2636/// Return the relative variable importance, normalized to all
2637/// variables together having the importance 1. The importance in
2638/// evaluated as the total separation-gain that this variable had in
2639/// the decision trees (weighted by the number of events)
2640
2642{
2643 fVariableImportance.resize(GetNvar());
2644 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++) {
2645 fVariableImportance[ivar]=0;
2646 }
2647 Double_t sum=0;
2648 for (UInt_t itree = 0; itree < GetNTrees(); itree++) {
2649 std::vector<Double_t> relativeImportance(fForest[itree]->GetVariableImportance());
2650 for (UInt_t i=0; i< relativeImportance.size(); i++) {
2651 fVariableImportance[i] += fBoostWeights[itree] * relativeImportance[i];
2652 }
2653 }
2654
2655 for (UInt_t ivar=0; ivar< fVariableImportance.size(); ivar++){
2656 fVariableImportance[ivar] = TMath::Sqrt(fVariableImportance[ivar]);
2657 sum += fVariableImportance[ivar];
2658 }
2659 for (UInt_t ivar=0; ivar< fVariableImportance.size(); ivar++) fVariableImportance[ivar] /= sum;
2660
2661 return fVariableImportance;
2662}
2663
2664////////////////////////////////////////////////////////////////////////////////
2665/// Returns the measure for the variable importance of variable "ivar"
2666/// which is later used in GetVariableImportance() to calculate the
2667/// relative variable importances.
2668
2670{
2671 std::vector<Double_t> relativeImportance = this->GetVariableImportance();
2672 if (ivar < (UInt_t)relativeImportance.size()) return relativeImportance[ivar];
2673 else Log() << kFATAL << "<GetVariableImportance> ivar = " << ivar << " is out of range " << Endl;
2674
2675 return -1;
2676}
2677
2678////////////////////////////////////////////////////////////////////////////////
2679/// Compute ranking of input variables
2680
2682{
2683 // create the ranking object
2684 fRanking = new Ranking( GetName(), "Variable Importance" );
2685 vector< Double_t> importance(this->GetVariableImportance());
2686
2687 for (UInt_t ivar=0; ivar<GetNvar(); ivar++) {
2688
2689 fRanking->AddRank( Rank( GetInputLabel(ivar), importance[ivar] ) );
2690 }
2691
2692 return fRanking;
2693}
2694
2695////////////////////////////////////////////////////////////////////////////////
2696/// Get help message text.
2697
2699{
2700 Log() << Endl;
2701 Log() << gTools().Color("bold") << "--- Short description:" << gTools().Color("reset") << Endl;
2702 Log() << Endl;
2703 Log() << "Boosted Decision Trees are a collection of individual decision" << Endl;
2704 Log() << "trees which form a multivariate classifier by (weighted) majority " << Endl;
2705 Log() << "vote of the individual trees. Consecutive decision trees are " << Endl;
2706 Log() << "trained using the original training data set with re-weighted " << Endl;
2707 Log() << "events. By default, the AdaBoost method is employed, which gives " << Endl;
2708 Log() << "events that were misclassified in the previous tree a larger " << Endl;
2709 Log() << "weight in the training of the following tree." << Endl;
2710 Log() << Endl;
2711 Log() << "Decision trees are a sequence of binary splits of the data sample" << Endl;
2712 Log() << "using a single discriminant variable at a time. A test event " << Endl;
2713 Log() << "ending up after the sequence of left-right splits in a final " << Endl;
2714 Log() << "(\"leaf\") node is classified as either signal or background" << Endl;
2715 Log() << "depending on the majority type of training events in that node." << Endl;
2716 Log() << Endl;
2717 Log() << gTools().Color("bold") << "--- Performance optimisation:" << gTools().Color("reset") << Endl;
2718 Log() << Endl;
2719 Log() << "By the nature of the binary splits performed on the individual" << Endl;
2720 Log() << "variables, decision trees do not deal well with linear correlations" << Endl;
2721 Log() << "between variables (they need to approximate the linear split in" << Endl;
2722 Log() << "the two dimensional space by a sequence of splits on the two " << Endl;
2723 Log() << "variables individually). Hence decorrelation could be useful " << Endl;
2724 Log() << "to optimise the BDT performance." << Endl;
2725 Log() << Endl;
2726 Log() << gTools().Color("bold") << "--- Performance tuning via configuration options:" << gTools().Color("reset") << Endl;
2727 Log() << Endl;
2728 Log() << "The two most important parameters in the configuration are the " << Endl;
2729 Log() << "minimal number of events requested by a leaf node as percentage of the " <<Endl;
2730 Log() << " number of training events (option \"MinNodeSize\" replacing the actual number " << Endl;
2731 Log() << " of events \"nEventsMin\" as given in earlier versions" << Endl;
2732 Log() << "If this number is too large, detailed features " << Endl;
2733 Log() << "in the parameter space are hard to be modelled. If it is too small, " << Endl;
2734 Log() << "the risk to overtrain rises and boosting seems to be less effective" << Endl;
2735 Log() << " typical values from our current experience for best performance " << Endl;
2736 Log() << " are between 0.5(%) and 10(%) " << Endl;
2737 Log() << Endl;
2738 Log() << "The default minimal number is currently set to " << Endl;
2739 Log() << " max(20, (N_training_events / N_variables^2 / 10)) " << Endl;
2740 Log() << "and can be changed by the user." << Endl;
2741 Log() << Endl;
2742 Log() << "The other crucial parameter, the pruning strength (\"PruneStrength\")," << Endl;
2743 Log() << "is also related to overtraining. It is a regularisation parameter " << Endl;
2744 Log() << "that is used when determining after the training which splits " << Endl;
2745 Log() << "are considered statistically insignificant and are removed. The" << Endl;
2746 Log() << "user is advised to carefully watch the BDT screen output for" << Endl;
2747 Log() << "the comparison between efficiencies obtained on the training and" << Endl;
2748 Log() << "the independent test sample. They should be equal within statistical" << Endl;
2749 Log() << "errors, in order to minimize statistical fluctuations in different samples." << Endl;
2750}
2751
2752////////////////////////////////////////////////////////////////////////////////
2753/// Make ROOT-independent C++ class for classifier response (classifier-specific implementation).
2754
2755void TMVA::MethodBDT::MakeClassSpecific( std::ostream& fout, const TString& className ) const
2756{
2757 TString nodeName = className;
2758 nodeName.ReplaceAll("Read","");
2759 nodeName.Append("Node");
2760 // write BDT-specific classifier response
2761 fout << " std::vector<"<<nodeName<<"*> fForest; // i.e. root nodes of decision trees" << std::endl;
2762 fout << " std::vector<double> fBoostWeights; // the weights applied in the individual boosts" << std::endl;
2763 fout << "};" << std::endl << std::endl;
2764
2765 if(GetAnalysisType() == Types::kMulticlass) {
2766 fout << "std::vector<double> ReadBDTG::GetMulticlassValues__( const std::vector<double>& inputValues ) const" << std::endl;
2767 fout << "{" << std::endl;
2768 fout << " uint nClasses = " << DataInfo().GetNClasses() << ";" << std::endl;
2769 fout << " std::vector<double> fMulticlassReturnVal;" << std::endl;
2770 fout << " fMulticlassReturnVal.reserve(nClasses);" << std::endl;
2771 fout << std::endl;
2772 fout << " std::vector<double> temp(nClasses);" << std::endl;
2773 fout << " auto forestSize = fForest.size();" << std::endl;
2774 fout << " // trees 0, nClasses, 2*nClasses, ... belong to class 0" << std::endl;
2775 fout << " // trees 1, nClasses+1, 2*nClasses+1, ... belong to class 1 and so forth" << std::endl;
2776 fout << " uint classOfTree = 0;" << std::endl;
2777 fout << " for (uint itree = 0; itree < forestSize; ++itree) {" << std::endl;
2778 fout << " BDTGNode *current = fForest[itree];" << std::endl;
2779 fout << " while (current->GetNodeType() == 0) { //intermediate node" << std::endl;
2780 fout << " if (current->GoesRight(inputValues)) current=(BDTGNode*)current->GetRight();" << std::endl;
2781 fout << " else current=(BDTGNode*)current->GetLeft();" << std::endl;
2782 fout << " }" << std::endl;
2783 fout << " temp[classOfTree] += current->GetResponse();" << std::endl;
2784 fout << " if (++classOfTree == nClasses) classOfTree = 0; // cheap modulo" << std::endl;
2785 fout << " }" << std::endl;
2786 fout << std::endl;
2787 fout << " // we want to calculate sum of exp(temp[j] - temp[i]) for all i,j (i!=j)" << std::endl;
2788 fout << " // first calculate exp(), then replace minus with division." << std::endl;
2789 fout << " std::transform(temp.begin(), temp.end(), temp.begin(), [](double d){return exp(d);});" << std::endl;
2790 fout << std::endl;
2791 fout << " for(uint iClass=0; iClass<nClasses; iClass++){" << std::endl;
2792 fout << " double norm = 0.0;" << std::endl;
2793 fout << " for(uint j=0;j<nClasses;j++){" << std::endl;
2794 fout << " if(iClass!=j)" << std::endl;
2795 fout << " norm += temp[j] / temp[iClass];" << std::endl;
2796 fout << " }" << std::endl;
2797 fout << " fMulticlassReturnVal.push_back(1.0/(1.0+norm));" << std::endl;
2798 fout << " }" << std::endl;
2799 fout << std::endl;
2800 fout << " return fMulticlassReturnVal;" << std::endl;
2801 fout << "}" << std::endl;
2802 } else {
2803 fout << "double " << className << "::GetMvaValue__( const std::vector<double>& inputValues ) const" << std::endl;
2804 fout << "{" << std::endl;
2805 fout << " double myMVA = 0;" << std::endl;
2806 if (fDoPreselection){
2807 for (UInt_t ivar = 0; ivar< fIsLowBkgCut.size(); ivar++){
2808 if (fIsLowBkgCut[ivar]){
2809 fout << " if (inputValues["<<ivar<<"] < " << fLowBkgCut[ivar] << ") return -1; // is background preselection cut" << std::endl;
2810 }
2811 if (fIsLowSigCut[ivar]){
2812 fout << " if (inputValues["<<ivar<<"] < "<< fLowSigCut[ivar] << ") return 1; // is signal preselection cut" << std::endl;
2813 }
2814 if (fIsHighBkgCut[ivar]){
2815 fout << " if (inputValues["<<ivar<<"] > "<<fHighBkgCut[ivar] <<") return -1; // is background preselection cut" << std::endl;
2816 }
2817 if (fIsHighSigCut[ivar]){
2818 fout << " if (inputValues["<<ivar<<"] > "<<fHighSigCut[ivar]<<") return 1; // is signal preselection cut" << std::endl;
2819 }
2820 }
2821 }
2822
2823 if (fBoostType!="Grad"){
2824 fout << " double norm = 0;" << std::endl;
2825 }
2826 fout << " for (unsigned int itree=0; itree<fForest.size(); itree++){" << std::endl;
2827 fout << " "<<nodeName<<" *current = fForest[itree];" << std::endl;
2828 fout << " while (current->GetNodeType() == 0) { //intermediate node" << std::endl;
2829 fout << " if (current->GoesRight(inputValues)) current=("<<nodeName<<"*)current->GetRight();" << std::endl;
2830 fout << " else current=("<<nodeName<<"*)current->GetLeft();" << std::endl;
2831 fout << " }" << std::endl;
2832 if (fBoostType=="Grad"){
2833 fout << " myMVA += current->GetResponse();" << std::endl;
2834 }else{
2835 if (fUseYesNoLeaf) fout << " myMVA += fBoostWeights[itree] * current->GetNodeType();" << std::endl;
2836 else fout << " myMVA += fBoostWeights[itree] * current->GetPurity();" << std::endl;
2837 fout << " norm += fBoostWeights[itree];" << std::endl;
2838 }
2839 fout << " }" << std::endl;
2840 if (fBoostType=="Grad"){
2841 fout << " return 2.0/(1.0+exp(-2.0*myMVA))-1.0;" << std::endl;
2842 }
2843 else fout << " return myMVA /= norm;" << std::endl;
2844 fout << "}" << std::endl << std::endl;
2845 }
2846
2847 fout << "void " << className << "::Initialize()" << std::endl;
2848 fout << "{" << std::endl;
2849 fout << " double inf = std::numeric_limits<double>::infinity();" << std::endl;
2850 fout << " double nan = std::numeric_limits<double>::quiet_NaN();" << std::endl;
2851 //Now for each decision tree, write directly the constructors of the nodes in the tree structure
2852 for (UInt_t itree=0; itree<GetNTrees(); itree++) {
2853 fout << " // itree = " << itree << std::endl;
2854 fout << " fBoostWeights.push_back(" << fBoostWeights[itree] << ");" << std::endl;
2855 fout << " fForest.push_back( " << std::endl;
2856 this->MakeClassInstantiateNode((DecisionTreeNode*)fForest[itree]->GetRoot(), fout, className);
2857 fout <<" );" << std::endl;
2858 }
2859 fout << " return;" << std::endl;
2860 fout << "};" << std::endl;
2861 fout << std::endl;
2862 fout << "// Clean up" << std::endl;
2863 fout << "inline void " << className << "::Clear() " << std::endl;
2864 fout << "{" << std::endl;
2865 fout << " for (unsigned int itree=0; itree<fForest.size(); itree++) { " << std::endl;
2866 fout << " delete fForest[itree]; " << std::endl;
2867 fout << " }" << std::endl;
2868 fout << "}" << std::endl;
2869 fout << std::endl;
2870}
2871
2872////////////////////////////////////////////////////////////////////////////////
2873/// Specific class header.
2874
2875void TMVA::MethodBDT::MakeClassSpecificHeader( std::ostream& fout, const TString& className) const
2876{
2877 TString nodeName = className;
2878 nodeName.ReplaceAll("Read","");
2879 nodeName.Append("Node");
2880 fout << "#include <algorithm>" << std::endl;
2881 fout << "#include <limits>" << std::endl;
2882 fout << std::endl;
2883 //fout << "#ifndef NN" << std::endl; commented out on purpose see next line
2884 fout << "#define NN new "<<nodeName << std::endl; // NN definition depends on individual methods. Important to have NO #ifndef if several BDT methods compile together
2885 //fout << "#endif" << std::endl; commented out on purpose see previous line
2886 fout << std::endl;
2887 fout << "#ifndef "<<nodeName<<"__def" << std::endl;
2888 fout << "#define "<<nodeName<<"__def" << std::endl;
2889 fout << std::endl;
2890 fout << "class "<<nodeName<<" {" << std::endl;
2891 fout << std::endl;
2892 fout << "public:" << std::endl;
2893 fout << std::endl;
2894 fout << " // constructor of an essentially \"empty\" node floating in space" << std::endl;
2895 fout << " "<<nodeName<<" ( "<<nodeName<<"* left,"<<nodeName<<"* right," << std::endl;
2896 if (fUseFisherCuts){
2897 fout << " int nFisherCoeff," << std::endl;
2898 for (UInt_t i=0;i<GetNVariables()+1;i++){
2899 fout << " double fisherCoeff"<<i<<"," << std::endl;
2900 }
2901 }
2902 fout << " int selector, double cutValue, bool cutType, " << std::endl;
2903 fout << " int nodeType, double purity, double response ) :" << std::endl;
2904 fout << " fLeft ( left )," << std::endl;
2905 fout << " fRight ( right )," << std::endl;
2906 if (fUseFisherCuts) fout << " fNFisherCoeff ( nFisherCoeff )," << std::endl;
2907 fout << " fSelector ( selector )," << std::endl;
2908 fout << " fCutValue ( cutValue )," << std::endl;
2909 fout << " fCutType ( cutType )," << std::endl;
2910 fout << " fNodeType ( nodeType )," << std::endl;
2911 fout << " fPurity ( purity )," << std::endl;
2912 fout << " fResponse ( response ){" << std::endl;
2913 if (fUseFisherCuts){
2914 for (UInt_t i=0;i<GetNVariables()+1;i++){
2915 fout << " fFisherCoeff.push_back(fisherCoeff"<<i<<");" << std::endl;
2916 }
2917 }
2918 fout << " }" << std::endl << std::endl;
2919 fout << " virtual ~"<<nodeName<<"();" << std::endl << std::endl;
2920 fout << " // test event if it descends the tree at this node to the right" << std::endl;
2921 fout << " virtual bool GoesRight( const std::vector<double>& inputValues ) const;" << std::endl;
2922 fout << " "<<nodeName<<"* GetRight( void ) {return fRight; };" << std::endl << std::endl;
2923 fout << " // test event if it descends the tree at this node to the left " << std::endl;
2924 fout << " virtual bool GoesLeft ( const std::vector<double>& inputValues ) const;" << std::endl;
2925 fout << " "<<nodeName<<"* GetLeft( void ) { return fLeft; }; " << std::endl << std::endl;
2926 fout << " // return S/(S+B) (purity) at this node (from training)" << std::endl << std::endl;
2927 fout << " double GetPurity( void ) const { return fPurity; } " << std::endl;
2928 fout << " // return the node type" << std::endl;
2929 fout << " int GetNodeType( void ) const { return fNodeType; }" << std::endl;
2930 fout << " double GetResponse(void) const {return fResponse;}" << std::endl << std::endl;
2931 fout << "private:" << std::endl << std::endl;
2932 fout << " "<<nodeName<<"* fLeft; // pointer to the left daughter node" << std::endl;
2933 fout << " "<<nodeName<<"* fRight; // pointer to the right daughter node" << std::endl;
2934 if (fUseFisherCuts){
2935 fout << " int fNFisherCoeff; // =0 if this node doesn't use fisher, else =nvar+1 " << std::endl;
2936 fout << " std::vector<double> fFisherCoeff; // the fisher coeff (offset at the last element)" << std::endl;
2937 }
2938 fout << " int fSelector; // index of variable used in node selection (decision tree) " << std::endl;
2939 fout << " double fCutValue; // cut value applied on this node to discriminate bkg against sig" << std::endl;
2940 fout << " bool fCutType; // true: if event variable > cutValue ==> signal , false otherwise" << std::endl;
2941 fout << " int fNodeType; // Type of node: -1 == Bkg-leaf, 1 == Signal-leaf, 0 = internal " << std::endl;
2942 fout << " double fPurity; // Purity of node from training"<< std::endl;
2943 fout << " double fResponse; // Regression response value of node" << std::endl;
2944 fout << "}; " << std::endl;
2945 fout << std::endl;
2946 fout << "//_______________________________________________________________________" << std::endl;
2947 fout << " "<<nodeName<<"::~"<<nodeName<<"()" << std::endl;
2948 fout << "{" << std::endl;
2949 fout << " if (fLeft != NULL) delete fLeft;" << std::endl;
2950 fout << " if (fRight != NULL) delete fRight;" << std::endl;
2951 fout << "}; " << std::endl;
2952 fout << std::endl;
2953 fout << "//_______________________________________________________________________" << std::endl;
2954 fout << "bool "<<nodeName<<"::GoesRight( const std::vector<double>& inputValues ) const" << std::endl;
2955 fout << "{" << std::endl;
2956 fout << " // test event if it descends the tree at this node to the right" << std::endl;
2957 fout << " bool result;" << std::endl;
2958 if (fUseFisherCuts){
2959 fout << " if (fNFisherCoeff == 0){" << std::endl;
2960 fout << " result = (inputValues[fSelector] >= fCutValue );" << std::endl;
2961 fout << " }else{" << std::endl;
2962 fout << " double fisher = fFisherCoeff.at(fFisherCoeff.size()-1);" << std::endl;
2963 fout << " for (unsigned int ivar=0; ivar<fFisherCoeff.size()-1; ivar++)" << std::endl;
2964 fout << " fisher += fFisherCoeff.at(ivar)*inputValues.at(ivar);" << std::endl;
2965 fout << " result = fisher > fCutValue;" << std::endl;
2966 fout << " }" << std::endl;
2967 }else{
2968 fout << " result = (inputValues[fSelector] >= fCutValue );" << std::endl;
2969 }
2970 fout << " if (fCutType == true) return result; //the cuts are selecting Signal ;" << std::endl;
2971 fout << " else return !result;" << std::endl;
2972 fout << "}" << std::endl;
2973 fout << std::endl;
2974 fout << "//_______________________________________________________________________" << std::endl;
2975 fout << "bool "<<nodeName<<"::GoesLeft( const std::vector<double>& inputValues ) const" << std::endl;
2976 fout << "{" << std::endl;
2977 fout << " // test event if it descends the tree at this node to the left" << std::endl;
2978 fout << " if (!this->GoesRight(inputValues)) return true;" << std::endl;
2979 fout << " else return false;" << std::endl;
2980 fout << "}" << std::endl;
2981 fout << std::endl;
2982 fout << "#endif" << std::endl;
2983 fout << std::endl;
2984}
2985
2986////////////////////////////////////////////////////////////////////////////////
2987/// Recursively descends a tree and writes the node instance to the output stream.
2988
2989void TMVA::MethodBDT::MakeClassInstantiateNode( DecisionTreeNode *n, std::ostream& fout, const TString& className ) const
2990{
2991 if (n == NULL) {
2992 Log() << kFATAL << "MakeClassInstantiateNode: started with undefined node" <<Endl;
2993 return ;
2994 }
2995 fout << "NN("<<std::endl;
2996 if (n->GetLeft() != NULL){
2997 this->MakeClassInstantiateNode( (DecisionTreeNode*)n->GetLeft() , fout, className);
2998 }
2999 else {
3000 fout << "0";
3001 }
3002 fout << ", " <<std::endl;
3003 if (n->GetRight() != NULL){
3004 this->MakeClassInstantiateNode( (DecisionTreeNode*)n->GetRight(), fout, className );
3005 }
3006 else {
3007 fout << "0";
3008 }
3009 fout << ", " << std::endl
3010 << std::setprecision(6);
3011 if (fUseFisherCuts){
3012 fout << n->GetNFisherCoeff() << ", ";
3013 for (UInt_t i=0; i< GetNVariables()+1; i++) {
3014 if (n->GetNFisherCoeff() == 0 ){
3015 fout << "0, ";
3016 }else{
3017 fout << n->GetFisherCoeff(i) << ", ";
3018 }
3019 }
3020 }
3021 fout << n->GetSelector() << ", "
3022 << n->GetCutValue() << ", "
3023 << n->GetCutType() << ", "
3024 << n->GetNodeType() << ", "
3025 << n->GetPurity() << ","
3026 << n->GetResponse() << ") ";
3027}
3028
3029////////////////////////////////////////////////////////////////////////////////
3030/// Find useful preselection cuts that will be applied before
3031/// and Decision Tree training.. (and of course also applied
3032/// in the GetMVA .. --> -1 for background +1 for Signal)
3033
3034void TMVA::MethodBDT::DeterminePreselectionCuts(const std::vector<const TMVA::Event*>& eventSample)
3035{
3036 Double_t nTotS = 0.0, nTotB = 0.0;
3037
3038 std::vector<TMVA::BDTEventWrapper> bdtEventSample;
3039
3040 fIsLowSigCut.assign(GetNvar(),kFALSE);
3041 fIsLowBkgCut.assign(GetNvar(),kFALSE);
3042 fIsHighSigCut.assign(GetNvar(),kFALSE);
3043 fIsHighBkgCut.assign(GetNvar(),kFALSE);
3044
3045 fLowSigCut.assign(GetNvar(),0.); // ---------------| --> in var is signal (accept all above lower cut)
3046 fLowBkgCut.assign(GetNvar(),0.); // ---------------| --> in var is bkg (accept all above lower cut)
3047 fHighSigCut.assign(GetNvar(),0.); // <-- | -------------- in var is signal (accept all blow cut)
3048 fHighBkgCut.assign(GetNvar(),0.); // <-- | -------------- in var is blg (accept all blow cut)
3049
3050
3051 // Initialize (un)weighted counters for signal & background
3052 // Construct a list of event wrappers that point to the original data
3053 for( std::vector<const TMVA::Event*>::const_iterator it = eventSample.begin(); it != eventSample.end(); ++it ) {
3054 if (DataInfo().IsSignal(*it)){
3055 nTotS += (*it)->GetWeight();
3056 }
3057 else {
3058 nTotB += (*it)->GetWeight();
3059 }
3060 bdtEventSample.push_back(TMVA::BDTEventWrapper(*it));
3061 }
3062
3063 for( UInt_t ivar = 0; ivar < GetNvar(); ivar++ ) { // loop over all discriminating variables
3064 TMVA::BDTEventWrapper::SetVarIndex(ivar); // select the variable to sort by
3065 std::sort( bdtEventSample.begin(),bdtEventSample.end() ); // sort the event data
3066
3067 Double_t bkgWeightCtr = 0.0, sigWeightCtr = 0.0;
3068 std::vector<TMVA::BDTEventWrapper>::iterator it = bdtEventSample.begin(), it_end = bdtEventSample.end();
3069 for( ; it != it_end; ++it ) {
3070 if (DataInfo().IsSignal(**it))
3071 sigWeightCtr += (**it)->GetWeight();
3072 else
3073 bkgWeightCtr += (**it)->GetWeight();
3074 // Store the accumulated signal (background) weights
3075 it->SetCumulativeWeight(false,bkgWeightCtr);
3076 it->SetCumulativeWeight(true,sigWeightCtr);
3077 }
3078
3079 //variable that determines how "exact" you cut on the preselection found in the training data. Here I chose
3080 //1% of the variable range...
3081 Double_t dVal = (DataInfo().GetVariableInfo(ivar).GetMax() - DataInfo().GetVariableInfo(ivar).GetMin())/100. ;
3082 Double_t nSelS, nSelB, effS=0.05, effB=0.05, rejS=0.05, rejB=0.05;
3084 // Locate the optimal cut for this (ivar-th) variable
3085
3086
3087
3088 for(UInt_t iev = 1; iev < bdtEventSample.size(); iev++) {
3089 //dVal = bdtEventSample[iev].GetVal() - bdtEventSample[iev-1].GetVal();
3090
3091 nSelS = bdtEventSample[iev].GetCumulativeWeight(true);
3092 nSelB = bdtEventSample[iev].GetCumulativeWeight(false);
3093 // you look for some 100% efficient pre-selection cut to remove background.. i.e. nSelS=0 && nSelB>5%nTotB or ( nSelB=0 nSelS>5%nTotS)
3094 tmpEffS=nSelS/nTotS;
3095 tmpEffB=nSelB/nTotB;
3096 tmpRejS=1-tmpEffS;
3097 tmpRejB=1-tmpEffB;
3098 if (nSelS==0 && tmpEffB>effB) {effB=tmpEffB; fLowBkgCut[ivar] = bdtEventSample[iev].GetVal() - dVal; fIsLowBkgCut[ivar]=kTRUE;}
3099 else if (nSelB==0 && tmpEffS>effS) {effS=tmpEffS; fLowSigCut[ivar] = bdtEventSample[iev].GetVal() - dVal; fIsLowSigCut[ivar]=kTRUE;}
3100 else if (nSelB==nTotB && tmpRejS>rejS) {rejS=tmpRejS; fHighSigCut[ivar] = bdtEventSample[iev].GetVal() + dVal; fIsHighSigCut[ivar]=kTRUE;}
3101 else if (nSelS==nTotS && tmpRejB>rejB) {rejB=tmpRejB; fHighBkgCut[ivar] = bdtEventSample[iev].GetVal() + dVal; fIsHighBkgCut[ivar]=kTRUE;}
3102
3103 }
3104 }
3105
3106 Log() << kDEBUG << " \tfound and suggest the following possible pre-selection cuts " << Endl;
3107 if (fDoPreselection) Log() << kDEBUG << "\tthe training will be done after these cuts... and GetMVA value returns +1, (-1) for a signal (bkg) event that passes these cuts" << Endl;
3108 else Log() << kDEBUG << "\tas option DoPreselection was not used, these cuts however will not be performed, but the training will see the full sample"<<Endl;
3109 for (UInt_t ivar=0; ivar < GetNvar(); ivar++ ) { // loop over all discriminating variables
3110 if (fIsLowBkgCut[ivar]){
3111 Log() << kDEBUG << " \tfound cut: Bkg if var " << ivar << " < " << fLowBkgCut[ivar] << Endl;
3112 }
3113 if (fIsLowSigCut[ivar]){
3114 Log() << kDEBUG << " \tfound cut: Sig if var " << ivar << " < " << fLowSigCut[ivar] << Endl;
3115 }
3116 if (fIsHighBkgCut[ivar]){
3117 Log() << kDEBUG << " \tfound cut: Bkg if var " << ivar << " > " << fHighBkgCut[ivar] << Endl;
3118 }
3119 if (fIsHighSigCut[ivar]){
3120 Log() << kDEBUG << " \tfound cut: Sig if var " << ivar << " > " << fHighSigCut[ivar] << Endl;
3121 }
3122 }
3123
3124 return;
3125}
3126
3127////////////////////////////////////////////////////////////////////////////////
3128/// Apply the preselection cuts before even bothering about any
3129/// Decision Trees in the GetMVA .. --> -1 for background +1 for Signal
3130
3132{
3133 Double_t result=0;
3134
3135 for (UInt_t ivar=0; ivar < GetNvar(); ivar++ ) { // loop over all discriminating variables
3136 if (fIsLowBkgCut[ivar]){
3137 if (ev->GetValue(ivar) < fLowBkgCut[ivar]) result = -1; // is background
3138 }
3139 if (fIsLowSigCut[ivar]){
3140 if (ev->GetValue(ivar) < fLowSigCut[ivar]) result = 1; // is signal
3141 }
3142 if (fIsHighBkgCut[ivar]){
3143 if (ev->GetValue(ivar) > fHighBkgCut[ivar]) result = -1; // is background
3144 }
3145 if (fIsHighSigCut[ivar]){
3146 if (ev->GetValue(ivar) > fHighSigCut[ivar]) result = 1; // is signal
3147 }
3148 }
3149
3150 return result;
3151}
3152
#define REGISTER_METHOD(CLASS)
for example
#define d(i)
Definition RSha256.hxx:102
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
unsigned int UInt_t
Unsigned integer 4 bytes (unsigned int)
Definition RtypesCore.h:60
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
void Print(GNN_Data &d, std::string txt="")
const_iterator begin() const
const_iterator end() const
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2290
Int_t GetN() const
Definition TGraph.h:131
virtual void Set(Int_t n)
Set number of points in the graph Existing coordinates are preserved New coordinates above fNpoints a...
Definition TGraph.cxx:2225
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:879
1-D histogram with an int per channel (see TH1 documentation)
Definition TH1.h:797
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
2-D histogram with a float per channel (see TH1 documentation)
Definition TH2.h:307
Service class for 2-D histogram classes.
Definition TH2.h:30
Absolute Deviation BDT Loss Function.
static void SetVarIndex(Int_t iVar)
static Config & Instance()
static function: returns TMVA instance
Definition Config.cxx:97
Implementation of the CrossEntropy as separation criterion.
Class that contains all the data information.
Definition DataSetInfo.h:62
static void SetIsTraining(bool on)
Implementation of a Decision Tree.
static DecisionTree * CreateFromXML(void *node, UInt_t tmva_Version_Code=262657)
re-create a new tree (decision tree or search tree) from XML
void SetTarget(UInt_t itgt, Float_t value)
set the target value (dimension itgt) to value
Definition Event.cxx:367
Implementation of the GiniIndex With Laplace correction as separation criterion.
Implementation of the GiniIndex as separation criterion.
Definition GiniIndex.h:63
Huber BDT Loss Function.
The TMVA::Interval Class.
Definition Interval.h:61
Least Squares BDT Loss Function.
The TMVA::Interval Class.
Definition LogInterval.h:83
Double_t GetMvaValue(Double_t *err=nullptr, Double_t *errUpper=nullptr) override
const std::vector< Float_t > & GetMulticlassValues() override
Get the multiclass MVA response for the BDT classifier.
void Init(void) override
Common initialisation with defaults for the BDT-Method.
void MakeClassSpecificHeader(std::ostream &, const TString &) const override
Specific class header.
void AddWeightsXMLTo(void *parent) const override
Write weights to XML.
static const Int_t fgDebugLevel
debug level determining some printout/control plots etc.
Definition MethodBDT.h:302
MethodBDT(const TString &jobName, const TString &methodTitle, DataSetInfo &theData, const TString &theOption="")
The standard constructor for the "boosted decision trees".
void BoostMonitor(Int_t iTree)
Fills the ROCIntegral vs Itree from the testSample for the monitoring plots during the training .
Double_t AdaBoostR2(std::vector< const TMVA::Event * > &, DecisionTree *dt)
Adaption of the AdaBoost to regression problems (see H.Drucker 1997).
Double_t PrivateGetMvaValue(const TMVA::Event *ev, Double_t *err=nullptr, Double_t *errUpper=nullptr, UInt_t useNTrees=0)
Return the MVA value (range [-1;1]) that classifies the event according to the majority vote from the...
void SetTuneParameters(std::map< TString, Double_t > tuneParameters) override
Set the tuning parameters according to the argument.
std::map< TString, Double_t > OptimizeTuningParameters(TString fomType="ROCIntegral", TString fitType="FitGA") override
Call the Optimizer with the set of parameters and ranges that are meant to be tuned.
Bool_t HasAnalysisType(Types::EAnalysisType type, UInt_t numberClasses, UInt_t numberTargets) override
BDT can handle classification with multiple classes and regression with one regression-target.
LossFunctionBDT * fRegressionLossFunctionBDTG
Definition MethodBDT.h:299
void DeterminePreselectionCuts(const std::vector< const TMVA::Event * > &eventSample)
Find useful preselection cuts that will be applied before and Decision Tree training.
Double_t GradBoost(std::vector< const TMVA::Event * > &, DecisionTree *dt, UInt_t cls=0)
Calculate the desired response value for each region.
Double_t AdaCost(std::vector< const TMVA::Event * > &, DecisionTree *dt)
The AdaCost boosting algorithm takes a simple cost Matrix (currently fixed for all events....
void Train(void) override
BDT training.
Double_t Boost(std::vector< const TMVA::Event * > &, DecisionTree *dt, UInt_t cls=0)
Apply the boosting algorithm (the algorithm is selecte via the "option" given in the constructor.
Double_t TestTreeQuality(DecisionTree *dt)
Test the tree quality.. in terms of Misclassification.
Double_t Bagging()
Call it boot-strapping, re-sampling or whatever you like, in the end it is nothing else but applying ...
void ReadWeightsFromStream(std::istream &istr) override
Read the weights (BDT coefficients).
void MakeClassSpecific(std::ostream &, const TString &) const override
Make ROOT-independent C++ class for classifier response (classifier-specific implementation).
void UpdateTargets(std::vector< const TMVA::Event * > &, UInt_t cls=0)
Calculate residual for all events.
void UpdateTargetsRegression(std::vector< const TMVA::Event * > &, Bool_t first=kFALSE)
Calculate residuals for all events and update targets for next iter.
Double_t GradBoostRegression(std::vector< const TMVA::Event * > &, DecisionTree *dt)
Implementation of M_TreeBoost using any loss function as described by Friedman 1999.
virtual ~MethodBDT(void)
Destructor.
Double_t GetGradBoostMVA(const TMVA::Event *e, UInt_t nTrees)
Returns MVA value: -1 for background, 1 for signal.
Double_t RegBoost(std::vector< const TMVA::Event * > &, DecisionTree *dt)
A special boosting only for Regression (not implemented).
void InitEventSample()
Initialize the event sample (i.e. reset the boost-weights... etc).
void DeclareCompatibilityOptions() override
Options that are used ONLY for the READER to ensure backward compatibility.
void WriteMonitoringHistosToFile(void) const override
Here we could write some histograms created during the processing to the output file.
void DeclareOptions() override
Define the options (their key words).
Double_t ApplyPreselectionCuts(const Event *ev)
Apply the preselection cuts before even bothering about any Decision Trees in the GetMVA .
void SetMinNodeSize(Double_t sizeInPercent)
void PreProcessNegativeEventWeights()
O.k.
void ReadWeightsFromXML(void *parent) override
Reads the BDT from the xml file.
void Reset(void) override
Reset the method, as if it had just been instantiated (forget all training etc.).
void GetHelpMessage() const override
Get help message text.
void MakeClassInstantiateNode(DecisionTreeNode *n, std::ostream &fout, const TString &className) const
Recursively descends a tree and writes the node instance to the output stream.
Double_t AdaBoost(std::vector< const TMVA::Event * > &, DecisionTree *dt)
The AdaBoost implementation.
TTree * fMonitorNtuple
monitoring ntuple
Definition MethodBDT.h:264
std::vector< Double_t > GetVariableImportance()
Return the relative variable importance, normalized to all variables together having the importance 1...
void InitGradBoost(std::vector< const TMVA::Event * > &)
Initialize targets for first tree.
const std::vector< Float_t > & GetRegressionValues() override
Get the regression value generated by the BDTs.
void GetBaggedSubSample(std::vector< const TMVA::Event * > &)
Fills fEventSample with fBaggedSampleFraction*NEvents random training events.
const Ranking * CreateRanking() override
Compute ranking of input variables.
SeparationBase * fSepType
the separation used in node splitting
Definition MethodBDT.h:229
void ProcessOptions() override
The option string is decoded, for available options see "DeclareOptions".
Virtual base Class for all MVA method.
Definition MethodBase.h:111
virtual void DeclareCompatibilityOptions()
options that are used ONLY for the READER to ensure backward compatibility they are hence without any...
Implementation of the MisClassificationError as separation criterion.
std::map< TString, Double_t > optimize()
PDF wrapper for histograms; uses user-defined spline interpolation.
Definition PDF.h:63
@ kSpline3
Definition PDF.h:70
Ranking for variables in method (implementation)
Definition Ranking.h:48
Class that is the base-class for a vector of result.
Definition Results.h:57
Implementation of the SdivSqrtSplusB as separation criterion.
Timing information for training and evaluation of MVA methods.
Definition Timer.h:58
const TString & Color(const TString &)
human readable color strings
Definition Tools.cxx:828
void ReadAttr(void *node, const char *, T &value)
read attribute from xml
Definition Tools.h:329
void * GetChild(void *parent, const char *childname=nullptr)
get child node
Definition Tools.cxx:1150
void AddAttr(void *node, const char *, const T &value, Int_t precision=16)
add attribute to xml
Definition Tools.h:347
void * AddChild(void *parent, const char *childname, const char *content=nullptr, bool isRootNode=false)
add child node
Definition Tools.cxx:1124
void UsefulSortAscending(std::vector< std::vector< Double_t > > &, std::vector< TString > *vs=nullptr)
sort 2D vector (AND in parallel a TString vector) in such a way that the "first vector is sorted" and...
Definition Tools.cxx:538
std::vector< TMatrixDSym * > * CalcCovarianceMatrices(const std::vector< Event * > &events, Int_t maxCls, VariableTransformBase *transformBase=nullptr)
compute covariance matrices
Definition Tools.cxx:1513
void * GetNextChild(void *prevchild, const char *childname=nullptr)
XML helpers.
Definition Tools.cxx:1162
Singleton class for Global types used by TMVA.
Definition Types.h:71
@ kMulticlass
Definition Types.h:129
@ kClassification
Definition Types.h:127
@ kMaxAnalysisType
Definition Types.h:131
@ kRegression
Definition Types.h:128
@ kTraining
Definition Types.h:143
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:964
virtual Int_t Read(const char *name)
Read contents of object with specified name from the current directory.
Definition TObject.cxx:667
Random number generator class based on M.
Definition TRandom3.h:27
Basic string class.
Definition TString.h:138
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2384
A TTree represents a columnar dataset.
Definition TTree.h:89
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
TSeq< unsigned int > TSeqU
Definition TSeq.hxx:204
create variable transformations
Tools & gTools()
MsgLogger & Endl(MsgLogger &ml)
Definition MsgLogger.h:148
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Definition TMath.h:720
Int_t FloorNint(Double_t x)
Returns the nearest integer of TMath::Floor(x).
Definition TMath.h:697
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:767
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Int_t CeilNint(Double_t x)
Returns the nearest integer of TMath::Ceil(x).
Definition TMath.h:685
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339