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