Logo ROOT   6.14/05
Reference Guide
RuleEnsemble.cxx
Go to the documentation of this file.
1 // @(#)root/tmva $Id$
2 // Author: Andreas Hoecker, Joerg Stelzer, Fredrik Tegenfeldt, Helge Voss
3 
4 /**********************************************************************************
5  * Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
6  * Package: TMVA *
7  * Class : RuleEnsemble *
8  * Web : http://tmva.sourceforge.net *
9  * *
10  * Description: *
11  * A class generating an ensemble of rules *
12  * Input: a forest of decision trees *
13  * Output: an ensemble of rules *
14  * *
15  * Authors (alphabetical): *
16  * Fredrik Tegenfeldt <Fredrik.Tegenfeldt@cern.ch> - Iowa State U., USA *
17  * Helge Voss <Helge.Voss@cern.ch> - MPI-KP Heidelberg, GER *
18  * *
19  * Copyright (c) 2005: *
20  * CERN, Switzerland *
21  * Iowa State U. *
22  * MPI-K Heidelberg, Germany *
23  * *
24  * Redistribution and use in source and binary forms, with or without *
25  * modification, are permitted according to the terms listed in LICENSE *
26  * (http://tmva.sourceforge.net/LICENSE) *
27  **********************************************************************************/
28 
29 /*! \class TMVA::RuleEnsemble
30 \ingroup TMVA
31 
32 */
33 #include "TMVA/RuleEnsemble.h"
34 
35 #include "TMVA/DataSetInfo.h"
36 #include "TMVA/DecisionTree.h"
37 #include "TMVA/Event.h"
38 #include "TMVA/MethodBase.h"
39 #include "TMVA/MethodRuleFit.h"
40 #include "TMVA/MsgLogger.h"
41 #include "TMVA/RuleFit.h"
42 #include "TMVA/Tools.h"
43 #include "TMVA/Types.h"
44 
45 #include "TRandom3.h"
46 #include "TH1F.h"
47 
48 #include <algorithm>
49 #include <cstdlib>
50 #include <iomanip>
51 #include <list>
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// constructor
55 
57  : fLearningModel ( kFull )
58  , fImportanceCut ( 0 )
59  , fLinQuantile ( 0.025 ) // default quantile for killing outliers in linear terms
60  , fOffset ( 0 )
61  , fAverageSupport ( 0.8 )
62  , fAverageRuleSigma( 0.4 ) // default value - used if only linear model is chosen
63  , fRuleFSig ( 0 )
64  , fRuleNCave ( 0 )
65  , fRuleNCsig ( 0 )
66  , fRuleMinDist ( 1e-3 ) // closest allowed 'distance' between two rules
67  , fNRulesGenerated ( 0 )
68  , fEvent ( 0 )
69  , fEventCacheOK ( true )
70  , fRuleMapOK ( true )
71  , fRuleMapInd0 ( 0 )
72  , fRuleMapInd1 ( 0 )
73  , fRuleMapEvents ( 0 )
74  , fLogger( new MsgLogger("RuleFit") )
75 {
76  Initialize( rf );
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 /// copy constructor
81 
83  : fAverageSupport ( 1 )
84  , fEvent(0)
85  , fRuleMapEvents(0)
86  , fRuleFit(0)
87  , fLogger( new MsgLogger("RuleFit") )
88 {
89  Copy( other );
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// constructor
94 
96  : fLearningModel ( kFull )
97  , fImportanceCut ( 0 )
98  , fLinQuantile ( 0.025 ) // default quantile for killing outliers in linear terms
99  , fOffset ( 0 )
100  , fImportanceRef ( 1.0 )
101  , fAverageSupport ( 0.8 )
102  , fAverageRuleSigma( 0.4 ) // default value - used if only linear model is chosen
103  , fRuleFSig ( 0 )
104  , fRuleNCave ( 0 )
105  , fRuleNCsig ( 0 )
106  , fRuleMinDist ( 1e-3 ) // closest allowed 'distance' between two rules
107  , fNRulesGenerated ( 0 )
108  , fEvent ( 0 )
109  , fEventCacheOK ( true )
110  , fRuleMapOK ( true )
111  , fRuleMapInd0 ( 0 )
112  , fRuleMapInd1 ( 0 )
113  , fRuleMapEvents ( 0 )
114  , fRuleFit ( 0 )
115  , fLogger( new MsgLogger("RuleFit") )
116 {
117 }
118 
119 ////////////////////////////////////////////////////////////////////////////////
120 /// destructor
121 
123 {
124  for ( std::vector<Rule *>::iterator itrRule = fRules.begin(); itrRule != fRules.end(); ++itrRule ) {
125  delete *itrRule;
126  }
127  // NOTE: Should not delete the histos fLinPDFB/S since they are delete elsewhere
128  delete fLogger;
129 }
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// Initializes all member variables with default values
133 
135 {
136  SetAverageRuleSigma(0.4); // default value - used if only linear model is chosen
137  fRuleFit = rf;
138  UInt_t nvars = GetMethodBase()->GetNvar();
139  fVarImportance.clear();
140  fLinPDFB.clear();
141  fLinPDFS.clear();
142  //
143  fVarImportance.resize( nvars,0.0 );
144  fLinPDFB.resize( nvars,0 );
145  fLinPDFS.resize( nvars,0 );
146  fImportanceRef = 1.0;
147  for (UInt_t i=0; i<nvars; i++) { // a priori all linear terms are equally valid
148  fLinTermOK.push_back(kTRUE);
149  }
150 }
151 
152 ////////////////////////////////////////////////////////////////////////////////
153 
154 void TMVA::RuleEnsemble::SetMsgType( EMsgType t ) {
155  fLogger->SetMinType(t);
156 }
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 ///
160 /// Get a pointer to the original MethodRuleFit.
161 
163 {
164  return ( fRuleFit==0 ? 0:fRuleFit->GetMethodRuleFit());
165 }
166 
167 ////////////////////////////////////////////////////////////////////////////////
168 ///
169 /// Get a pointer to the original MethodRuleFit.
170 
172 {
173  return ( fRuleFit==0 ? 0:fRuleFit->GetMethodBase());
174 }
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 /// create model
178 
180 {
182 
183  MakeLinearTerms();
184 
185  MakeRuleMap();
186 
187  CalcRuleSupport();
188 
189  RuleStatistics();
190 
191  PrintRuleGen();
192 }
193 
194 ////////////////////////////////////////////////////////////////////////////////
195 ///
196 /// Calculates sqrt(Sum(a_i^2)), i=1..N (NOTE do not include a0)
197 
199 {
200  Int_t ncoeffs = fRules.size();
201  if (ncoeffs<1) return 0;
202  //
203  Double_t sum2=0;
204  Double_t val;
205  for (Int_t i=0; i<ncoeffs; i++) {
206  val = fRules[i]->GetCoefficient();
207  sum2 += val*val;
208  }
209  return sum2;
210 }
211 
212 ////////////////////////////////////////////////////////////////////////////////
213 /// reset all rule coefficients
214 
216 {
217  fOffset = 0.0;
218  UInt_t nrules = fRules.size();
219  for (UInt_t i=0; i<nrules; i++) {
220  fRules[i]->SetCoefficient(0.0);
221  }
222 }
223 
224 ////////////////////////////////////////////////////////////////////////////////
225 /// set all rule coefficients
226 
227 void TMVA::RuleEnsemble::SetCoefficients( const std::vector< Double_t > & v )
228 {
229  UInt_t nrules = fRules.size();
230  if (v.size()!=nrules) {
231  Log() << kFATAL << "<SetCoefficients> - BUG TRAP - input vector wrong size! It is = " << v.size()
232  << " when it should be = " << nrules << Endl;
233  }
234  for (UInt_t i=0; i<nrules; i++) {
235  fRules[i]->SetCoefficient(v[i]);
236  }
237 }
238 
239 ////////////////////////////////////////////////////////////////////////////////
240 /// Retrieve all rule coefficients
241 
242 void TMVA::RuleEnsemble::GetCoefficients( std::vector< Double_t > & v )
243 {
244  UInt_t nrules = fRules.size();
245  v.resize(nrules);
246  if (nrules==0) return;
247  //
248  for (UInt_t i=0; i<nrules; i++) {
249  v[i] = (fRules[i]->GetCoefficient());
250  }
251 }
252 
253 ////////////////////////////////////////////////////////////////////////////////
254 /// get list of training events from the rule fitter
255 
256 const std::vector<const TMVA::Event*>* TMVA::RuleEnsemble::GetTrainingEvents() const
257 {
258  return &(fRuleFit->GetTrainingEvents());
259 }
260 
261 ////////////////////////////////////////////////////////////////////////////////
262 /// get the training event from the rule fitter
263 
265 {
266  return fRuleFit->GetTrainingEvent(i);
267 }
268 
269 ////////////////////////////////////////////////////////////////////////////////
270 /// remove rules that behave similar
271 
273 {
274  Log() << kVERBOSE << "Removing similar rules; distance = " << fRuleMinDist << Endl;
275 
276  UInt_t nrulesIn = fRules.size();
278  std::vector< Char_t > removeMe( nrulesIn,false ); // <--- stores boolean
279 
280  Int_t nrem = 0;
281  Int_t remind=-1;
282  Double_t r;
283 
284  for (UInt_t i=0; i<nrulesIn; i++) {
285  if (!removeMe[i]) {
286  first = fRules[i];
287  for (UInt_t k=i+1; k<nrulesIn; k++) {
288  if (!removeMe[k]) {
289  second = fRules[k];
290  Bool_t equal = first->Equal(*second,kTRUE,fRuleMinDist);
291  if (equal) {
292  r = gRandom->Rndm();
293  remind = (r>0.5 ? k:i); // randomly select rule
294  }
295  else {
296  remind = -1;
297  }
298 
299  if (remind>-1) {
300  if (!removeMe[remind]) {
301  removeMe[remind] = true;
302  nrem++;
303  }
304  }
305  }
306  }
307  }
308  }
309  UInt_t ind = 0;
310  Rule *theRule;
311  for (UInt_t i=0; i<nrulesIn; i++) {
312  if (removeMe[i]) {
313  theRule = fRules[ind];
314 #if _MSC_VER >= 1400
315  fRules.erase( std::vector<Rule *>::iterator(&fRules[ind], &fRules) );
316 #else
317  fRules.erase( fRules.begin() + ind );
318 #endif
319  delete theRule;
320  ind--;
321  }
322  ind++;
323  }
324  UInt_t nrulesOut = fRules.size();
325  Log() << kVERBOSE << "Removed " << nrulesIn - nrulesOut << " out of " << nrulesIn << " rules" << Endl;
326 }
327 
328 ////////////////////////////////////////////////////////////////////////////////
329 /// cleanup rules
330 
332 {
333  UInt_t nrules = fRules.size();
334  if (nrules==0) return;
335  Log() << kVERBOSE << "Removing rules with relative importance < " << fImportanceCut << Endl;
336  if (fImportanceCut<=0) return;
337  //
338  // Mark rules to be removed
339  //
340  Rule *therule;
341  Int_t ind=0;
342  for (UInt_t i=0; i<nrules; i++) {
343  if (fRules[ind]->GetRelImportance()<fImportanceCut) {
344  therule = fRules[ind];
345 #if _MSC_VER >= 1400
346  fRules.erase( std::vector<Rule *>::iterator(&fRules[ind], &fRules) );
347 #else
348  fRules.erase( fRules.begin() + ind );
349 #endif
350  delete therule;
351  ind--;
352  }
353  ind++;
354  }
355  Log() << kINFO << "Removed " << nrules-ind << " out of a total of " << nrules
356  << " rules with importance < " << fImportanceCut << Endl;
357 }
358 
359 ////////////////////////////////////////////////////////////////////////////////
360 /// cleanup linear model
361 
363 {
364  UInt_t nlin = fLinNorm.size();
365  if (nlin==0) return;
366  Log() << kVERBOSE << "Removing linear terms with relative importance < " << fImportanceCut << Endl;
367  //
368  fLinTermOK.clear();
369  for (UInt_t i=0; i<nlin; i++) {
371  }
372 }
373 
374 ////////////////////////////////////////////////////////////////////////////////
375 /// calculate the support for all rules
376 
378 {
379  Log() << kVERBOSE << "Evaluating Rule support" << Endl;
380  Double_t s,t,stot,ttot,ssb;
381  Double_t ssig, sbkg, ssum;
382  Int_t indrule=0;
383  stot = 0;
384  ttot = 0;
385  // reset to default values
386  SetAverageRuleSigma(0.4);
387  const std::vector<const Event *> *events = GetTrainingEvents();
388  Double_t nrules = static_cast<Double_t>(fRules.size());
389  Double_t ew;
390  //
391  if ((nrules>0) && (events->size()>0)) {
392  for ( std::vector< Rule * >::iterator itrRule=fRules.begin(); itrRule!=fRules.end(); ++itrRule ) {
393  s=0.0;
394  ssig=0.0;
395  sbkg=0.0;
396  for ( std::vector<const Event * >::const_iterator itrEvent=events->begin(); itrEvent!=events->end(); ++itrEvent ) {
397  if ((*itrRule)->EvalEvent( *(*itrEvent) )) {
398  ew = (*itrEvent)->GetWeight();
399  s += ew;
400  if (GetMethodRuleFit()->DataInfo().IsSignal(*itrEvent)) ssig += ew;
401  else sbkg += ew;
402  }
403  }
404  //
405  s = s/fRuleFit->GetNEveEff();
406  t = s*(1.0-s);
407  t = (t<0 ? 0:sqrt(t));
408  stot += s;
409  ttot += t;
410  ssum = ssig+sbkg;
411  ssb = (ssum>0 ? Double_t(ssig)/Double_t(ssig+sbkg) : 0.0 );
412  (*itrRule)->SetSupport(s);
413  (*itrRule)->SetNorm(t);
414  (*itrRule)->SetSSB( ssb );
415  (*itrRule)->SetSSBNeve(Double_t(ssig+sbkg));
416  indrule++;
417  }
418  fAverageSupport = stot/nrules;
420  Log() << kVERBOSE << "Standard deviation of support = " << fAverageRuleSigma << Endl;
421  Log() << kVERBOSE << "Average rule support = " << fAverageSupport << Endl;
422  }
423 }
424 
425 ////////////////////////////////////////////////////////////////////////////////
426 /// calculate the importance of each rule
427 
429 {
430  Double_t maxRuleImp = CalcRuleImportance();
431  Double_t maxLinImp = CalcLinImportance();
432  Double_t maxImp = (maxRuleImp>maxLinImp ? maxRuleImp : maxLinImp);
433  SetImportanceRef( maxImp );
434 }
435 
436 ////////////////////////////////////////////////////////////////////////////////
437 /// set reference importance
438 
440 {
441  for ( UInt_t i=0; i<fRules.size(); i++ ) {
442  fRules[i]->SetImportanceRef(impref);
443  }
444  fImportanceRef = impref;
445 }
446 ////////////////////////////////////////////////////////////////////////////////
447 /// calculate importance of each rule
448 
450 {
451  Double_t maxImp=-1.0;
452  Double_t imp;
453  Int_t nrules = fRules.size();
454  for ( int i=0; i<nrules; i++ ) {
455  fRules[i]->CalcImportance();
456  imp = fRules[i]->GetImportance();
457  if (imp>maxImp) maxImp = imp;
458  }
459  for ( Int_t i=0; i<nrules; i++ ) {
460  fRules[i]->SetImportanceRef(maxImp);
461  }
462 
463  return maxImp;
464 }
465 
466 ////////////////////////////////////////////////////////////////////////////////
467 /// calculate the linear importance for each rule
468 
470 {
471  Double_t maxImp=-1.0;
472  UInt_t nvars = fLinCoefficients.size();
473  fLinImportance.resize(nvars,0.0);
474  if (!DoLinear()) return maxImp;
475  //
476  // The linear importance is:
477  // I = |b_x|*sigma(x)
478  // Note that the coefficients in fLinCoefficients are for the normalized x
479  // => b'_x * x' = b'_x * sigma(r)*x/sigma(x)
480  // => b_x = b'_x*sigma(r)/sigma(x)
481  // => I = |b'_x|*sigma(r)
482  //
483  Double_t imp;
484  for ( UInt_t i=0; i<nvars; i++ ) {
486  fLinImportance[i] = imp;
487  if (imp>maxImp) maxImp = imp;
488  }
489  return maxImp;
490 }
491 
492 ////////////////////////////////////////////////////////////////////////////////
493 /// Calculates variable importance using eq (35) in RuleFit paper by Friedman et.al
494 
496 {
497  Log() << kVERBOSE << "Compute variable importance" << Endl;
498  Double_t rimp;
499  UInt_t nrules = fRules.size();
500  if (GetMethodBase()==0) Log() << kFATAL << "RuleEnsemble::CalcVarImportance() - should not be here!" << Endl;
501  UInt_t nvars = GetMethodBase()->GetNvar();
502  UInt_t nvarsUsed;
503  Double_t rimpN;
504  fVarImportance.resize(nvars,0);
505  // rules
506  if (DoRules()) {
507  for ( UInt_t ind=0; ind<nrules; ind++ ) {
508  rimp = fRules[ind]->GetImportance();
509  nvarsUsed = fRules[ind]->GetNumVarsUsed();
510  if (nvarsUsed<1)
511  Log() << kFATAL << "<CalcVarImportance> Variables for importance calc!!!??? A BUG!" << Endl;
512  rimpN = (nvarsUsed > 0 ? rimp/nvarsUsed:0.0);
513  for ( UInt_t iv=0; iv<nvars; iv++ ) {
514  if (fRules[ind]->ContainsVariable(iv)) {
515  fVarImportance[iv] += rimpN;
516  }
517  }
518  }
519  }
520  // linear terms
521  if (DoLinear()) {
522  for ( UInt_t iv=0; iv<fLinTermOK.size(); iv++ ) {
523  if (fLinTermOK[iv]) fVarImportance[iv] += fLinImportance[iv];
524  }
525  }
526  //
527  // Make variable importance relative the strongest variable
528  //
529  Double_t maximp = 0.0;
530  for ( UInt_t iv=0; iv<nvars; iv++ ) {
531  if ( fVarImportance[iv] > maximp ) maximp = fVarImportance[iv];
532  }
533  if (maximp>0) {
534  for ( UInt_t iv=0; iv<nvars; iv++ ) {
535  fVarImportance[iv] *= 1.0/maximp;
536  }
537  }
538 }
539 
540 ////////////////////////////////////////////////////////////////////////////////
541 /// set rules
542 ///
543 /// first clear all
544 
545 void TMVA::RuleEnsemble::SetRules( const std::vector<Rule *> & rules )
546 {
547  DeleteRules();
548  //
549  fRules.resize(rules.size());
550  for (UInt_t i=0; i<fRules.size(); i++) {
551  fRules[i] = rules[i];
552  }
554 }
555 
556 ////////////////////////////////////////////////////////////////////////////////
557 /// Makes rules from the given decision tree.
558 /// First node in all rules is ALWAYS the root node.
559 
560 void TMVA::RuleEnsemble::MakeRules( const std::vector< const DecisionTree *> & forest )
561 {
562  fRules.clear();
563  if (!DoRules()) return;
564  //
565  Int_t nrulesCheck=0;
566  Int_t nrules;
567  Int_t nendn;
568  Double_t sumnendn=0;
569  Double_t sumn2=0;
570  //
571  // UInt_t prevs;
572  UInt_t ntrees = forest.size();
573  for ( UInt_t ind=0; ind<ntrees; ind++ ) {
574  // prevs = fRules.size();
575  MakeRulesFromTree( forest[ind] );
576  nrules = CalcNRules( forest[ind] );
577  nendn = (nrules/2) + 1;
578  sumnendn += nendn;
579  sumn2 += nendn*nendn;
580  nrulesCheck += nrules;
581  }
582  Double_t nmean = (ntrees>0) ? sumnendn/ntrees : 0;
583  Double_t nsigm = TMath::Sqrt( gTools().ComputeVariance(sumn2,sumnendn,ntrees) );
584  Double_t ndev = 2.0*(nmean-2.0-nsigm)/(nmean-2.0+nsigm);
585  //
586  Log() << kVERBOSE << "Average number of end nodes per tree = " << nmean << Endl;
587  if (ntrees>1) Log() << kVERBOSE << "sigma of ditto ( ~= mean-2 ?) = "
588  << nsigm
589  << Endl;
590  Log() << kVERBOSE << "Deviation from exponential model = " << ndev << Endl;
591  Log() << kVERBOSE << "Corresponds to L (eq. 13, RuleFit ppr) = " << nmean << Endl;
592  // a BUG trap
593  if (nrulesCheck != static_cast<Int_t>(fRules.size())) {
594  Log() << kFATAL
595  << "BUG! number of generated and possible rules do not match! N(rules) = " << fRules.size()
596  << " != " << nrulesCheck << Endl;
597  }
598  Log() << kVERBOSE << "Number of generated rules: " << fRules.size() << Endl;
599 
600  // save initial number of rules
601  fNRulesGenerated = fRules.size();
602 
604 
606 
607 }
608 
609 ////////////////////////////////////////////////////////////////////////////////
610 /// Make the linear terms as in eq 25, ref 2
611 /// For this the b and (1-b) quantiles are needed
612 
614 {
615  if (!DoLinear()) return;
616 
617  const std::vector<const Event *> *events = GetTrainingEvents();
618  UInt_t neve = events->size();
619  UInt_t nvars = ((*events)[0])->GetNVariables(); // Event -> GetNVariables();
620  Double_t val,ew;
621  typedef std::pair< Double_t, Int_t> dataType;
622  typedef std::pair< Double_t, dataType > dataPoint;
623 
624  std::vector< std::vector<dataPoint> > vardata(nvars);
625  std::vector< Double_t > varsum(nvars,0.0);
626  std::vector< Double_t > varsum2(nvars,0.0);
627  // first find stats of all variables
628  // vardata[v][i].first -> value of var <v> in event <i>
629  // vardata[v][i].second.first -> the event weight
630  // vardata[v][i].second.second -> the event type
631  for (UInt_t i=0; i<neve; i++) {
632  ew = ((*events)[i])->GetWeight();
633  for (UInt_t v=0; v<nvars; v++) {
634  val = ((*events)[i])->GetValue(v);
635  vardata[v].push_back( dataPoint( val, dataType(ew,((*events)[i])->GetClass()) ) );
636  }
637  }
638  //
639  fLinDP.clear();
640  fLinDM.clear();
641  fLinCoefficients.clear();
642  fLinNorm.clear();
643  fLinDP.resize(nvars,0);
644  fLinDM.resize(nvars,0);
645  fLinCoefficients.resize(nvars,0);
646  fLinNorm.resize(nvars,0);
647 
648  Double_t averageWeight = neve ? fRuleFit->GetNEveEff()/static_cast<Double_t>(neve) : 0;
649  // sort and find limits
650  Double_t stdl;
651 
652  // find normalisation given in ref 2 after eq 26
653  Double_t lx;
654  Double_t nquant;
655  Double_t neff;
656  UInt_t indquantM;
657  UInt_t indquantP;
658 
659  MethodBase *fMB=const_cast<MethodBase *>(fRuleFit->GetMethodBase());
660 
661  for (UInt_t v=0; v<nvars; v++) {
662  varsum[v] = 0;
663  varsum2[v] = 0;
664  //
665  std::sort( vardata[v].begin(),vardata[v].end() );
666  nquant = fLinQuantile*fRuleFit->GetNEveEff(); // quantile = 0.025
667  neff=0;
668  UInt_t ie=0;
669  // first scan for lower quantile (including weights)
670  while ( (ie<neve) && (neff<nquant) ) {
671  neff += vardata[v][ie].second.first;
672  ie++;
673  }
674  indquantM = (ie==0 ? 0:ie-1);
675  // now for upper quantile
676  ie = neve;
677  neff=0;
678  while ( (ie>0) && (neff<nquant) ) {
679  ie--;
680  neff += vardata[v][ie].second.first;
681  }
682  indquantP = (ie==neve ? ie=neve-1:ie);
683  //
684  fLinDM[v] = vardata[v][indquantM].first; // delta-
685  fLinDP[v] = vardata[v][indquantP].first; // delta+
686 
687  if(!fMB->IsSilentFile())
688  {
689  if (fLinPDFB[v]) delete fLinPDFB[v];
690  if (fLinPDFS[v]) delete fLinPDFS[v];
691  fLinPDFB[v] = new TH1F(Form("bkgvar%d",v),"bkg temphist",40,fLinDM[v],fLinDP[v]);
692  fLinPDFS[v] = new TH1F(Form("sigvar%d",v),"sig temphist",40,fLinDM[v],fLinDP[v]);
693  fLinPDFB[v]->Sumw2();
694  fLinPDFS[v]->Sumw2();
695  }
696  //
697  Int_t type;
698  const Double_t w = 1.0/fRuleFit->GetNEveEff();
699  for (ie=0; ie<neve; ie++) {
700  val = vardata[v][ie].first;
701  ew = vardata[v][ie].second.first;
702  type = vardata[v][ie].second.second;
703  lx = TMath::Min( fLinDP[v], TMath::Max( fLinDM[v], val ) );
704  varsum[v] += ew*lx;
705  varsum2[v] += ew*lx*lx;
706  if(!fMB->IsSilentFile())
707  {
708  if (type==1) fLinPDFS[v]->Fill(lx,w*ew);
709  else fLinPDFB[v]->Fill(lx,w*ew);
710  }
711  }
712  //
713  // Get normalization.
714  //
715  stdl = TMath::Sqrt( (varsum2[v] - (varsum[v]*varsum[v]/fRuleFit->GetNEveEff()))/(fRuleFit->GetNEveEff()-averageWeight) );
716  fLinNorm[v] = CalcLinNorm(stdl);
717  }
718  // Save PDFs - for debugging purpose
719  if(!fMB->IsSilentFile())
720  {
721  for (UInt_t v=0; v<nvars; v++) {
722  fLinPDFS[v]->Write();
723  fLinPDFB[v]->Write();
724  }
725  }
726 }
727 
728 ////////////////////////////////////////////////////////////////////////////////
729 /// This function returns Pr( y = 1 | x ) for the linear terms.
730 
732 {
733  UInt_t nvars=fLinDP.size();
734 
735  Double_t fstot=0;
736  Double_t fbtot=0;
737  nsig = 0;
738  ntot = nvars;
739  for (UInt_t v=0; v<nvars; v++) {
740  Double_t val = fEventLinearVal[v];
741  Int_t bin = fLinPDFS[v]->FindBin(val);
742  fstot += fLinPDFS[v]->GetBinContent(bin);
743  fbtot += fLinPDFB[v]->GetBinContent(bin);
744  }
745  if (nvars<1) return 0;
746  ntot = (fstot+fbtot)/Double_t(nvars);
747  nsig = (fstot)/Double_t(nvars);
748  return fstot/(fstot+fbtot);
749 }
750 
751 ////////////////////////////////////////////////////////////////////////////////
752 /// This function returns Pr( y = 1 | x ) for rules.
753 /// The probability returned is normalized against the number of rules which are actually passed
754 
756 {
757  Double_t sump = 0;
758  Double_t sumok = 0;
759  Double_t sumz = 0;
760  Double_t ssb;
761  Double_t neve;
762  //
763  UInt_t nrules = fRules.size();
764  for (UInt_t ir=0; ir<nrules; ir++) {
765  if (fEventRuleVal[ir]>0) {
766  ssb = fEventRuleVal[ir]*GetRulesConst(ir)->GetSSB(); // S/(S+B) is evaluated in CalcRuleSupport() using ALL training events
767  neve = GetRulesConst(ir)->GetSSBNeve(); // number of events accepted by the rule
768  sump += ssb*neve; // number of signal events
769  sumok += neve; // total number of events passed
770  }
771  else sumz += 1.0; // all events
772  }
773 
774  nsig = sump;
775  ntot = sumok;
776  //
777  if (ntot>0) return nsig/ntot;
778  return 0.0;
779 }
780 
781 ////////////////////////////////////////////////////////////////////////////////
782 /// We want to estimate F* = argmin Eyx( L(y,F(x) ), min wrt F(x)
783 /// F(x) = FL(x) + FR(x) , linear and rule part
784 
786 {
787  SetEvent(e);
788  UpdateEventVal();
789  return FStar();
790 }
791 
792 ////////////////////////////////////////////////////////////////////////////////
793 /// We want to estimate F* = argmin Eyx( L(y,F(x) ), min wrt F(x)
794 /// F(x) = FL(x) + FR(x) , linear and rule part
795 
797 {
798  Double_t p=0;
799  Double_t nrs=0, nrt=0;
800  Double_t nls=0, nlt=0;
801  Double_t nt;
802  Double_t pr=0;
803  Double_t pl=0;
804 
805  // first calculate Pr(y=1|X) for rules and linear terms
806  if (DoLinear()) pl = PdfLinear(nls, nlt);
807  if (DoRules()) pr = PdfRule(nrs, nrt);
808  // nr(l)t=0 or 1
809  if ((nlt>0) && (nrt>0)) nt=2.0;
810  else nt=1.0;
811  p = (pl+pr)/nt;
812  return 2.0*p-1.0;
813 }
814 
815 ////////////////////////////////////////////////////////////////////////////////
816 /// calculate various statistics for this rule
817 
819 {
820  // TODO: NOT YET UPDATED FOR WEIGHTS
821  const std::vector<const Event *> *events = GetTrainingEvents();
822  const UInt_t neve = events->size();
823  const UInt_t nvars = GetMethodBase()->GetNvar();
824  const UInt_t nrules = fRules.size();
825  const Event *eveData;
826  // Flags
827  Bool_t sigRule;
828  Bool_t sigTag;
829  Bool_t bkgTag;
830  // Bool_t noTag;
831  Bool_t sigTrue;
832  Bool_t tagged;
833  // Counters
834  Int_t nsig=0;
835  Int_t nbkg=0;
836  Int_t ntag=0;
837  Int_t nss=0;
838  Int_t nsb=0;
839  Int_t nbb=0;
840  Int_t nbs=0;
841  std::vector<Int_t> varcnt;
842  // Clear vectors
843  fRulePSS.clear();
844  fRulePSB.clear();
845  fRulePBS.clear();
846  fRulePBB.clear();
847  fRulePTag.clear();
848  //
849  varcnt.resize(nvars,0);
850  fRuleVarFrac.clear();
851  fRuleVarFrac.resize(nvars,0);
852  //
853  for ( UInt_t i=0; i<nrules; i++ ) {
854  for ( UInt_t v=0; v<nvars; v++) {
855  if (fRules[i]->ContainsVariable(v)) varcnt[v]++; // count how often a variable occurs
856  }
857  sigRule = fRules[i]->IsSignalRule();
858  if (sigRule) { // rule is a signal rule (ie s/(s+b)>0.5)
859  nsig++;
860  }
861  else {
862  nbkg++;
863  }
864  // reset counters
865  nss=0;
866  nsb=0;
867  nbs=0;
868  nbb=0;
869  ntag=0;
870  // loop over all events
871  for (UInt_t e=0; e<neve; e++) {
872  eveData = (*events)[e];
873  tagged = fRules[i]->EvalEvent(*eveData);
874  sigTag = (tagged && sigRule); // it's tagged as a signal
875  bkgTag = (tagged && (!sigRule)); // ... as bkg
876  // noTag = !(sigTag || bkgTag); // ... not tagged
877  sigTrue = (eveData->GetClass() == 0); // true if event is true signal
878  if (tagged) {
879  ntag++;
880  if (sigTag && sigTrue) nss++;
881  if (sigTag && !sigTrue) nsb++;
882  if (bkgTag && sigTrue) nbs++;
883  if (bkgTag && !sigTrue) nbb++;
884  }
885  }
886  // Fill tagging probabilities
887  if (ntag>0 && neve > 0) { // should always be the case, but let's make sure and keep coverity quiet
888  fRulePTag.push_back(Double_t(ntag)/Double_t(neve));
889  fRulePSS.push_back(Double_t(nss)/Double_t(ntag));
890  fRulePSB.push_back(Double_t(nsb)/Double_t(ntag));
891  fRulePBS.push_back(Double_t(nbs)/Double_t(ntag));
892  fRulePBB.push_back(Double_t(nbb)/Double_t(ntag));
893  }
894  //
895  }
896  fRuleFSig = (nsig>0) ? static_cast<Double_t>(nsig)/static_cast<Double_t>(nsig+nbkg) : 0;
897  for ( UInt_t v=0; v<nvars; v++) {
898  fRuleVarFrac[v] = (nrules>0) ? Double_t(varcnt[v])/Double_t(nrules) : 0;
899  }
900 }
901 
902 ////////////////////////////////////////////////////////////////////////////////
903 /// calculate various statistics for this rule
904 
906 {
907  const UInt_t nrules = fRules.size();
908  Double_t nc;
909  Double_t sumNc =0;
910  Double_t sumNc2=0;
911  for ( UInt_t i=0; i<nrules; i++ ) {
912  nc = static_cast<Double_t>(fRules[i]->GetNcuts());
913  sumNc += nc;
914  sumNc2 += nc*nc;
915  }
916  fRuleNCave = 0.0;
917  fRuleNCsig = 0.0;
918  if (nrules>0) {
919  fRuleNCave = sumNc/nrules;
920  fRuleNCsig = TMath::Sqrt(gTools().ComputeVariance(sumNc2,sumNc,nrules));
921  }
922 }
923 
924 ////////////////////////////////////////////////////////////////////////////////
925 /// print rule generation info
926 
928 {
929  Log() << kHEADER << "-------------------RULE ENSEMBLE SUMMARY------------------------" << Endl;
930  const MethodRuleFit *mrf = GetMethodRuleFit();
931  if (mrf) Log() << kINFO << "Tree training method : " << (mrf->UseBoost() ? "AdaBoost":"Random") << Endl;
932  Log() << kINFO << "Number of events per tree : " << fRuleFit->GetNTreeSample() << Endl;
933  Log() << kINFO << "Number of trees : " << fRuleFit->GetForest().size() << Endl;
934  Log() << kINFO << "Number of generated rules : " << fNRulesGenerated << Endl;
935  Log() << kINFO << "Idem, after cleanup : " << fRules.size() << Endl;
936  Log() << kINFO << "Average number of cuts per rule : " << Form("%8.2f",fRuleNCave) << Endl;
937  Log() << kINFO << "Spread in number of cuts per rules : " << Form("%8.2f",fRuleNCsig) << Endl;
938  Log() << kVERBOSE << "Complexity : " << Form("%8.2f",fRuleNCave*fRuleNCsig) << Endl;
939  Log() << kINFO << "----------------------------------------------------------------" << Endl;
940  Log() << kINFO << Endl;
941 }
942 
943 ////////////////////////////////////////////////////////////////////////////////
944 /// print function
945 
947 {
948  const EMsgType kmtype=kINFO;
949  const Bool_t isDebug = (fLogger->GetMinType()<=kDEBUG);
950  //
951  Log() << kmtype << Endl;
952  Log() << kmtype << "================================================================" << Endl;
953  Log() << kmtype << " M o d e l " << Endl;
954  Log() << kmtype << "================================================================" << Endl;
955 
956  Int_t ind;
957  const UInt_t nvars = GetMethodBase()->GetNvar();
958  const Int_t nrules = fRules.size();
959  const Int_t printN = TMath::Min(10,nrules); //nrules+1;
960  Int_t maxL = 0;
961  for (UInt_t iv = 0; iv<fVarImportance.size(); iv++) {
962  if (GetMethodBase()->GetInputLabel(iv).Length() > maxL) maxL = GetMethodBase()->GetInputLabel(iv).Length();
963  }
964  //
965  if (isDebug) {
966  Log() << kDEBUG << "Variable importance:" << Endl;
967  for (UInt_t iv = 0; iv<fVarImportance.size(); iv++) {
968  Log() << kDEBUG << std::setw(maxL) << GetMethodBase()->GetInputLabel(iv)
969  << std::resetiosflags(std::ios::right)
970  << " : " << Form(" %3.3f",fVarImportance[iv]) << Endl;
971  }
972  }
973  //
974  Log() << kHEADER << "Offset (a0) = " << fOffset << Endl;
975  //
976  if (DoLinear()) {
977  if (fLinNorm.size() > 0) {
978  Log() << kmtype << "------------------------------------" << Endl;
979  Log() << kmtype << "Linear model (weights unnormalised)" << Endl;
980  Log() << kmtype << "------------------------------------" << Endl;
981  Log() << kmtype << std::setw(maxL) << "Variable"
982  << std::resetiosflags(std::ios::right) << " : "
983  << std::setw(11) << " Weights"
984  << std::resetiosflags(std::ios::right) << " : "
985  << "Importance"
986  << std::resetiosflags(std::ios::right)
987  << Endl;
988  Log() << kmtype << "------------------------------------" << Endl;
989  for ( UInt_t i=0; i<fLinNorm.size(); i++ ) {
990  Log() << kmtype << std::setw(std::max(maxL,8)) << GetMethodBase()->GetInputLabel(i);
991  if (fLinTermOK[i]) {
992  Log() << kmtype
993  << std::resetiosflags(std::ios::right)
994  << " : " << Form(" %10.3e",fLinCoefficients[i]*fLinNorm[i])
995  << " : " << Form(" %3.3f",fLinImportance[i]/fImportanceRef) << Endl;
996  }
997  else {
998  Log() << kmtype << "-> importance below threshold = "
999  << Form(" %3.3f",fLinImportance[i]/fImportanceRef) << Endl;
1000  }
1001  }
1002  Log() << kmtype << "------------------------------------" << Endl;
1003  }
1004  }
1005  else Log() << kmtype << "Linear terms were disabled" << Endl;
1006 
1007  if ((!DoRules()) || (nrules==0)) {
1008  if (!DoRules()) {
1009  Log() << kmtype << "Rule terms were disabled" << Endl;
1010  }
1011  else {
1012  Log() << kmtype << "Even though rules were included in the model, none passed! " << nrules << Endl;
1013  }
1014  }
1015  else {
1016  Log() << kmtype << "Number of rules = " << nrules << Endl;
1017  if (isDebug) {
1018  Log() << kmtype << "N(cuts) in rules, average = " << fRuleNCave << Endl;
1019  Log() << kmtype << " RMS = " << fRuleNCsig << Endl;
1020  Log() << kmtype << "Fraction of signal rules = " << fRuleFSig << Endl;
1021  Log() << kmtype << "Fraction of rules containing a variable (%):" << Endl;
1022  for ( UInt_t v=0; v<nvars; v++) {
1023  Log() << kmtype << " " << std::setw(maxL) << GetMethodBase()->GetInputLabel(v);
1024  Log() << kmtype << Form(" = %2.2f",fRuleVarFrac[v]*100.0) << " %" << Endl;
1025  }
1026  }
1027  //
1028  // Print out all rules sorted in importance
1029  //
1030  std::list< std::pair<double,int> > sortedImp;
1031  for (Int_t i=0; i<nrules; i++) {
1032  sortedImp.push_back( std::pair<double,int>( fRules[i]->GetImportance(),i ) );
1033  }
1034  sortedImp.sort();
1035  //
1036  Log() << kmtype << "Printing the first " << printN << " rules, ordered in importance." << Endl;
1037  int pind=0;
1038  for ( std::list< std::pair<double,int> >::reverse_iterator itpair = sortedImp.rbegin();
1039  itpair != sortedImp.rend(); ++itpair ) {
1040  ind = itpair->second;
1041  // if (pind==0) impref =
1042  // Log() << kmtype << "Rule #" <<
1043  // Log() << kmtype << *fRules[ind] << Endl;
1044  fRules[ind]->PrintLogger(Form("Rule %4d : ",pind+1));
1045  pind++;
1046  if (pind==printN) {
1047  if (nrules==printN) {
1048  Log() << kmtype << "All rules printed" << Endl;
1049  }
1050  else {
1051  Log() << kmtype << "Skipping the next " << nrules-printN << " rules" << Endl;
1052  }
1053  break;
1054  }
1055  }
1056  }
1057  Log() << kmtype << "================================================================" << Endl;
1058  Log() << kmtype << Endl;
1059 }
1060 
1061 ////////////////////////////////////////////////////////////////////////////////
1062 /// write rules to stream
1063 
1064 void TMVA::RuleEnsemble::PrintRaw( std::ostream & os ) const
1065 {
1066  Int_t dp = os.precision();
1067  UInt_t nrules = fRules.size();
1068  // std::sort(fRules.begin(),fRules.end());
1069  //
1070  os << "ImportanceCut= " << fImportanceCut << std::endl;
1071  os << "LinQuantile= " << fLinQuantile << std::endl;
1072  os << "AverageSupport= " << fAverageSupport << std::endl;
1073  os << "AverageRuleSigma= " << fAverageRuleSigma << std::endl;
1074  os << "Offset= " << fOffset << std::endl;
1075  os << "NRules= " << nrules << std::endl;
1076  for (UInt_t i=0; i<nrules; i++){
1077  os << "***Rule " << i << std::endl;
1078  (fRules[i])->PrintRaw(os);
1079  }
1080  UInt_t nlinear = fLinNorm.size();
1081  //
1082  os << "NLinear= " << fLinTermOK.size() << std::endl;
1083  for (UInt_t i=0; i<nlinear; i++) {
1084  os << "***Linear " << i << std::endl;
1085  os << std::setprecision(10) << (fLinTermOK[i] ? 1:0) << " "
1086  << fLinCoefficients[i] << " "
1087  << fLinNorm[i] << " "
1088  << fLinDM[i] << " "
1089  << fLinDP[i] << " "
1090  << fLinImportance[i] << " " << std::endl;
1091  }
1092  os << std::setprecision(dp);
1093 }
1094 
1095 ////////////////////////////////////////////////////////////////////////////////
1096 /// write rules to XML
1097 
1098 void* TMVA::RuleEnsemble::AddXMLTo(void* parent) const
1099 {
1100  void* re = gTools().AddChild( parent, "Weights" ); // this is the "RuleEnsemble"
1101 
1102  UInt_t nrules = fRules.size();
1103  UInt_t nlinear = fLinNorm.size();
1104  gTools().AddAttr( re, "NRules", nrules );
1105  gTools().AddAttr( re, "NLinear", nlinear );
1106  gTools().AddAttr( re, "LearningModel", (int)fLearningModel );
1107  gTools().AddAttr( re, "ImportanceCut", fImportanceCut );
1108  gTools().AddAttr( re, "LinQuantile", fLinQuantile );
1109  gTools().AddAttr( re, "AverageSupport", fAverageSupport );
1110  gTools().AddAttr( re, "AverageRuleSigma", fAverageRuleSigma );
1111  gTools().AddAttr( re, "Offset", fOffset );
1112  for (UInt_t i=0; i<nrules; i++) fRules[i]->AddXMLTo(re);
1113 
1114  for (UInt_t i=0; i<nlinear; i++) {
1115  void* lin = gTools().AddChild( re, "Linear" );
1116  gTools().AddAttr( lin, "OK", (fLinTermOK[i] ? 1:0) );
1117  gTools().AddAttr( lin, "Coeff", fLinCoefficients[i] );
1118  gTools().AddAttr( lin, "Norm", fLinNorm[i] );
1119  gTools().AddAttr( lin, "DM", fLinDM[i] );
1120  gTools().AddAttr( lin, "DP", fLinDP[i] );
1121  gTools().AddAttr( lin, "Importance", fLinImportance[i] );
1122  }
1123  return re;
1124 }
1125 
1126 ////////////////////////////////////////////////////////////////////////////////
1127 /// read rules from XML
1128 
1129 void TMVA::RuleEnsemble::ReadFromXML( void* wghtnode )
1130 {
1131  UInt_t nrules, nlinear;
1132  gTools().ReadAttr( wghtnode, "NRules", nrules );
1133  gTools().ReadAttr( wghtnode, "NLinear", nlinear );
1134  Int_t iLearningModel;
1135  gTools().ReadAttr( wghtnode, "LearningModel", iLearningModel );
1136  fLearningModel = (ELearningModel) iLearningModel;
1137  gTools().ReadAttr( wghtnode, "ImportanceCut", fImportanceCut );
1138  gTools().ReadAttr( wghtnode, "LinQuantile", fLinQuantile );
1139  gTools().ReadAttr( wghtnode, "AverageSupport", fAverageSupport );
1140  gTools().ReadAttr( wghtnode, "AverageRuleSigma", fAverageRuleSigma );
1141  gTools().ReadAttr( wghtnode, "Offset", fOffset );
1142 
1143  // read rules
1144  DeleteRules();
1145 
1146  UInt_t i = 0;
1147  fRules.resize( nrules );
1148  void* ch = gTools().GetChild( wghtnode );
1149  for (i=0; i<nrules; i++) {
1150  fRules[i] = new Rule();
1151  fRules[i]->SetRuleEnsemble( this );
1152  fRules[i]->ReadFromXML( ch );
1153 
1154  ch = gTools().GetNextChild(ch);
1155  }
1156 
1157  // read linear classifier (Fisher)
1158  fLinNorm .resize( nlinear );
1159  fLinTermOK .resize( nlinear );
1160  fLinCoefficients.resize( nlinear );
1161  fLinDP .resize( nlinear );
1162  fLinDM .resize( nlinear );
1163  fLinImportance .resize( nlinear );
1164 
1165  Int_t iok;
1166  i=0;
1167  while(ch) {
1168  gTools().ReadAttr( ch, "OK", iok );
1169  fLinTermOK[i] = (iok == 1);
1170  gTools().ReadAttr( ch, "Coeff", fLinCoefficients[i] );
1171  gTools().ReadAttr( ch, "Norm", fLinNorm[i] );
1172  gTools().ReadAttr( ch, "DM", fLinDM[i] );
1173  gTools().ReadAttr( ch, "DP", fLinDP[i] );
1174  gTools().ReadAttr( ch, "Importance", fLinImportance[i] );
1175 
1176  i++;
1177  ch = gTools().GetNextChild(ch);
1178  }
1179 }
1180 
1181 ////////////////////////////////////////////////////////////////////////////////
1182 /// read rule ensemble from stream
1183 
1184 void TMVA::RuleEnsemble::ReadRaw( std::istream & istr )
1185 {
1186  UInt_t nrules;
1187  //
1188  std::string dummy;
1189  Int_t idum;
1190  //
1191  // First block is general stuff
1192  //
1193  istr >> dummy >> fImportanceCut;
1194  istr >> dummy >> fLinQuantile;
1195  istr >> dummy >> fAverageSupport;
1196  istr >> dummy >> fAverageRuleSigma;
1197  istr >> dummy >> fOffset;
1198  istr >> dummy >> nrules;
1199  //
1200  // Now read in the rules
1201  //
1202  DeleteRules();
1203  //
1204  for (UInt_t i=0; i<nrules; i++){
1205  istr >> dummy >> idum; // read line "***Rule <ind>"
1206  fRules.push_back( new Rule() );
1207  (fRules.back())->SetRuleEnsemble( this );
1208  (fRules.back())->ReadRaw(istr);
1209  }
1210  //
1211  // and now the linear terms
1212  //
1213  UInt_t nlinear;
1214  //
1215  // coverity[tainted_data_argument]
1216  istr >> dummy >> nlinear;
1217  //
1218  fLinNorm .resize( nlinear );
1219  fLinTermOK .resize( nlinear );
1220  fLinCoefficients.resize( nlinear );
1221  fLinDP .resize( nlinear );
1222  fLinDM .resize( nlinear );
1223  fLinImportance .resize( nlinear );
1224  //
1225 
1226  Int_t iok;
1227  for (UInt_t i=0; i<nlinear; i++) {
1228  istr >> dummy >> idum;
1229  istr >> iok;
1230  fLinTermOK[i] = (iok==1);
1231  istr >> fLinCoefficients[i];
1232  istr >> fLinNorm[i];
1233  istr >> fLinDM[i];
1234  istr >> fLinDP[i];
1235  istr >> fLinImportance[i];
1236  }
1237 }
1238 
1239 ////////////////////////////////////////////////////////////////////////////////
1240 /// copy function
1241 
1243 {
1244  if(this != &other) {
1245  fRuleFit = other.GetRuleFit();
1246  fRuleMinDist = other.GetRuleMinDist();
1247  fOffset = other.GetOffset();
1248  fRules = other.GetRulesConst();
1249  fImportanceCut = other.GetImportanceCut();
1250  fVarImportance = other.GetVarImportance();
1251  fLearningModel = other.GetLearningModel();
1252  fLinQuantile = other.GetLinQuantile();
1253  fRuleNCsig = other.fRuleNCsig;
1255  fEventCacheOK = other.fEventCacheOK;
1258  fRuleFSig = other.fRuleFSig;
1259  fRuleMapInd0 = other.fRuleMapInd0;
1260  fRuleMapInd1 = other.fRuleMapInd1;
1261  fRuleMapOK = other.fRuleMapOK;
1262  fRuleNCave = other.fRuleNCave;
1263  }
1264 }
1265 
1266 ////////////////////////////////////////////////////////////////////////////////
1267 /// calculate the number of rules
1268 
1270 {
1271  if (dtree==0) return 0;
1272  Node *node = dtree->GetRoot();
1273  Int_t nendnodes = 0;
1274  FindNEndNodes( node, nendnodes );
1275  return 2*(nendnodes-1);
1276 }
1277 
1278 ////////////////////////////////////////////////////////////////////////////////
1279 /// find the number of leaf nodes
1280 
1281 void TMVA::RuleEnsemble::FindNEndNodes( const Node *node, Int_t & nendnodes )
1282 {
1283  if (node==0) return;
1284  if ((node->GetRight()==0) && (node->GetLeft()==0)) {
1285  ++nendnodes;
1286  return;
1287  }
1288  const Node *nodeR = node->GetRight();
1289  const Node *nodeL = node->GetLeft();
1290  FindNEndNodes( nodeR, nendnodes );
1291  FindNEndNodes( nodeL, nendnodes );
1292 }
1293 
1294 ////////////////////////////////////////////////////////////////////////////////
1295 /// create rules from the decision tree structure
1296 
1298 {
1299  Node *node = dtree->GetRoot();
1300  AddRule( node );
1301 }
1302 
1303 ////////////////////////////////////////////////////////////////////////////////
1304 /// add a new rule to the tree
1305 
1307 {
1308  if (node==0) return;
1309  if (node->GetParent()==0) { // it's a root node, don't make a rule
1310  AddRule( node->GetRight() );
1311  AddRule( node->GetLeft() );
1312  }
1313  else {
1314  Rule *rule = MakeTheRule(node);
1315  if (rule) {
1316  fRules.push_back( rule );
1317  AddRule( node->GetRight() );
1318  AddRule( node->GetLeft() );
1319  }
1320  else {
1321  Log() << kFATAL << "<AddRule> - ERROR failed in creating a rule! BUG!" << Endl;
1322  }
1323  }
1324 }
1325 
1326 ////////////////////////////////////////////////////////////////////////////////
1327 /// Make a Rule from a given Node.
1328 /// The root node (ie no parent) does not generate a Rule.
1329 /// The first node in a rule is always the root node => fNodes.size()>=2
1330 /// Each node corresponds to a cut and the cut value is given by the parent node.
1331 
1333 {
1334  if (node==0) {
1335  Log() << kFATAL << "<MakeTheRule> Input node is NULL. Should not happen. BUG!" << Endl;
1336  return 0;
1337  }
1338 
1339  if (node->GetParent()==0) { // a root node - ignore
1340  return 0;
1341  }
1342  //
1343  std::vector< const Node * > nodeVec;
1344  const Node *parent = node;
1345  //
1346  // Make list with the input node at the end:
1347  // <root node> <node1> <node2> ... <node given as argument>
1348  //
1349  nodeVec.push_back( node );
1350  while (parent!=0) {
1351  parent = parent->GetParent();
1352  if (!parent) continue;
1353  const DecisionTreeNode* dtn = dynamic_cast<const DecisionTreeNode*>(parent);
1354  if (dtn && dtn->GetSelector()>=0)
1355  nodeVec.insert( nodeVec.begin(), parent );
1356 
1357  }
1358  if (nodeVec.size()<2) {
1359  Log() << kFATAL << "<MakeTheRule> BUG! Inconsistent Rule!" << Endl;
1360  return 0;
1361  }
1362  Rule *rule = new Rule( this, nodeVec );
1363  rule->SetMsgType( Log().GetMinType() );
1364  return rule;
1365 }
1366 
1367 ////////////////////////////////////////////////////////////////////////////////
1368 /// Makes rule map for all events
1369 
1370 void TMVA::RuleEnsemble::MakeRuleMap(const std::vector<const Event *> *events, UInt_t ifirst, UInt_t ilast)
1371 {
1372  Log() << kVERBOSE << "Making Rule map for all events" << Endl;
1373  // make rule response map
1374  if (events==0) events = GetTrainingEvents();
1375  if ((ifirst==0) || (ilast==0) || (ifirst>ilast)) {
1376  ifirst = 0;
1377  ilast = events->size()-1;
1378  }
1379  // check if identical to previous call
1380  if ((events!=fRuleMapEvents) ||
1381  (ifirst!=fRuleMapInd0) ||
1382  (ilast !=fRuleMapInd1)) {
1383  fRuleMapOK = kFALSE;
1384  }
1385  //
1386  if (fRuleMapOK) {
1387  Log() << kVERBOSE << "<MakeRuleMap> Map is already valid" << Endl;
1388  return; // already cached
1389  }
1390  fRuleMapEvents = events;
1391  fRuleMapInd0 = ifirst;
1392  fRuleMapInd1 = ilast;
1393  // check number of rules
1394  UInt_t nrules = GetNRules();
1395  if (nrules==0) {
1396  Log() << kVERBOSE << "No rules found in MakeRuleMap()" << Endl;
1397  fRuleMapOK = kTRUE;
1398  return;
1399  }
1400  //
1401  // init map
1402  //
1403  std::vector<UInt_t> ruleind;
1404  fRuleMap.clear();
1405  for (UInt_t i=ifirst; i<=ilast; i++) {
1406  ruleind.clear();
1407  fRuleMap.push_back( ruleind );
1408  for (UInt_t r=0; r<nrules; r++) {
1409  if (fRules[r]->EvalEvent(*((*events)[i]))) {
1410  fRuleMap.back().push_back(r); // save only rules that are accepted
1411  }
1412  }
1413  }
1414  fRuleMapOK = kTRUE;
1415  Log() << kVERBOSE << "Made rule map for event# " << ifirst << " : " << ilast << Endl;
1416 }
1417 
1418 ////////////////////////////////////////////////////////////////////////////////
1419 /// std::ostream operator
1420 
1421 std::ostream& TMVA::operator<< ( std::ostream& os, const RuleEnsemble & rules )
1422 {
1423  os << "DON'T USE THIS - TO BE REMOVED" << std::endl;
1424  rules.Print();
1425  return os;
1426 }
Double_t fAverageSupport
Definition: RuleEnsemble.h:356
MsgLogger * fLogger
Definition: RuleEnsemble.h:385
void MakeRuleMap(const std::vector< const TMVA::Event *> *events=0, UInt_t ifirst=0, UInt_t ilast=0)
Makes rule map for all events.
J Friedman&#39;s RuleFit method.
Definition: MethodRuleFit.h:47
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:158
void SetEvent(const Event &e)
Definition: RuleEnsemble.h:145
std::vector< TH1F *> fLinPDFS
Definition: RuleEnsemble.h:352
const std::vector< const TMVA::Event *> & GetTrainingEvents() const
Definition: RuleFit.h:137
A class implementing various fits of rule ensembles.
Definition: RuleFit.h:45
RuleEnsemble()
constructor
UInt_t GetNvar() const
Definition: MethodBase.h:335
Rule * MakeTheRule(const Node *node)
Make a Rule from a given Node.
THist< 1, float, THistStatContent, THistStatUncertainty > TH1F
Definition: THist.hxx:285
Int_t CalcNRules(const TMVA::DecisionTree *dtree)
calculate the number of rules
EMsgType GetMinType() const
Definition: MsgLogger.h:71
Double_t CalcLinNorm(Double_t stdev)
Definition: RuleEnsemble.h:120
virtual ~RuleEnsemble()
destructor
const std::vector< TMVA::Rule * > & GetRulesConst() const
Definition: RuleEnsemble.h:267
Virtual base Class for all MVA method.
Definition: MethodBase.h:109
std::vector< Double_t > fLinDP
Definition: RuleEnsemble.h:347
const std::vector< Double_t > & GetVarImportance() const
Definition: RuleEnsemble.h:272
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
std::vector< Double_t > fRulePBB
Definition: RuleEnsemble.h:363
const Event * GetTrainingEvent(UInt_t i) const
Definition: RuleFit.h:132
void CleanupLinear()
cleanup linear model
void AddAttr(void *node, const char *, const T &value, Int_t precision=16)
add attribute to xml
Definition: Tools.h:353
virtual DecisionTreeNode * GetRoot() const
Definition: DecisionTree.h:93
void * AddChild(void *parent, const char *childname, const char *content=0, bool isRootNode=false)
add child node
Definition: Tools.cxx:1136
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
const TString & GetInputLabel(Int_t i) const
Definition: MethodBase.h:341
std::vector< TMVA::Rule *> fRules
Definition: RuleEnsemble.h:345
std::vector< Char_t > fLinTermOK
Definition: RuleEnsemble.h:346
Implementation of a rule.
Definition: Rule.h:48
void SetAverageRuleSigma(Double_t v)
Definition: RuleEnsemble.h:137
const std::vector< const TMVA::DecisionTree * > & GetForest() const
Definition: RuleFit.h:143
void SetMsgType(EMsgType t)
void SetImportanceRef(Double_t impref)
set reference importance
void RuleResponseStats()
calculate various statistics for this rule
double sqrt(double)
TClass * GetClass(T *)
Definition: TClass.h:577
Double_t GetRuleMinDist() const
Definition: RuleEnsemble.h:280
std::vector< Double_t > fLinNorm
Definition: RuleEnsemble.h:350
void RemoveSimilarRules()
remove rules that behave similar
void Copy(RuleEnsemble const &other)
copy function
virtual Node * GetRight() const
Definition: Node.h:88
UInt_t GetNRules() const
Definition: RuleEnsemble.h:266
Double_t PdfRule(Double_t &nsig, Double_t &ntot) const
This function returns Pr( y = 1 | x ) for rules.
virtual Node * GetLeft() const
Definition: Node.h:87
void * GetChild(void *parent, const char *childname=0)
get child node
Definition: Tools.cxx:1162
Double_t FStar() const
We want to estimate F* = argmin Eyx( L(y,F(x) ), min wrt F(x) F(x) = FL(x) + FR(x) ...
UInt_t GetClass() const
Definition: Event.h:81
void SetRules(const std::vector< TMVA::Rule *> &rules)
set rules
void Print() const
print function
std::vector< TH1F *> fLinPDFB
Definition: RuleEnsemble.h:351
void SetMinType(EMsgType minType)
Definition: MsgLogger.h:72
std::vector< Double_t > fRulePSB
Definition: RuleEnsemble.h:361
void CalcImportance()
calculate the importance of each rule
static constexpr double second
ELearningModel GetLearningModel() const
Definition: RuleEnsemble.h:262
std::vector< Double_t > fLinCoefficients
Definition: RuleEnsemble.h:349
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:533
void CleanupRules()
cleanup rules
const RuleFit * fRuleFit
Definition: RuleEnsemble.h:383
Double_t GetImportanceCut() const
Definition: RuleEnsemble.h:263
void * AddXMLTo(void *parent) const
write rules to XML
virtual Node * GetParent() const
Definition: Node.h:89
const MethodBase * GetMethodBase() const
Definition: RuleFit.h:149
void CalcVarImportance()
Calculates variable importance using eq (35) in RuleFit paper by Friedman et.al.
ROOT::R::TRInterface & r
Definition: Object.C:4
Bool_t DoLinear() const
Definition: RuleEnsemble.h:257
SVector< double, 2 > v
Definition: Dict.h:5
ELearningModel fLearningModel
Definition: RuleEnsemble.h:341
Double_t CalcRuleImportance()
calculate importance of each rule
void PrintRuleGen() const
print rule generation info
std::vector< Double_t > fLinDM
Definition: RuleEnsemble.h:348
void MakeRulesFromTree(const DecisionTree *dtree)
create rules from the decision tree structure
Double_t CoefficientRadius()
Calculates sqrt(Sum(a_i^2)), i=1..N (NOTE do not include a0)
void AddRule(const Node *node)
add a new rule to the tree
Implementation of a Decision Tree.
Definition: DecisionTree.h:64
unsigned int UInt_t
Definition: RtypesCore.h:42
std::ostream & operator<<(std::ostream &os, const BinaryTree &tree)
char * Form(const char *fmt,...)
Ssiz_t Length() const
Definition: TString.h:405
MsgLogger & Log() const
message logger
Definition: RuleEnsemble.h:386
Double_t PdfLinear(Double_t &nsig, Double_t &ntot) const
This function returns Pr( y = 1 | x ) for the linear terms.
void RuleStatistics()
calculate various statistics for this rule
void ReadAttr(void *node, const char *, T &value)
read attribute from xml
Definition: Tools.h:335
const RuleFit * GetRuleFit() const
Definition: RuleEnsemble.h:251
Tools & gTools()
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
Bool_t IsSilentFile()
Definition: MethodBase.h:370
Double_t fImportanceRef
Definition: RuleEnsemble.h:355
void ReadFromXML(void *wghtnode)
read rules from XML
const std::vector< const TMVA::Event * > * GetTrainingEvents() const
get list of training events from the rule fitter
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t GetLinQuantile() const
Definition: RuleEnsemble.h:274
void PrintRaw(std::ostream &os) const
write rules to stream
std::vector< Char_t > fEventRuleVal
Definition: RuleEnsemble.h:374
void FindNEndNodes(const TMVA::Node *node, Int_t &nendnodes)
find the number of leaf nodes
Double_t GetNEveEff() const
Definition: RuleFit.h:131
void Initialize(const RuleFit *rf)
Initializes all member variables with default values.
void SetCoefficients(const std::vector< Double_t > &v)
set all rule coefficients
double Double_t
Definition: RtypesCore.h:55
std::vector< Double_t > fVarImportance
Definition: RuleEnsemble.h:354
std::vector< Double_t > fEventLinearVal
Definition: RuleEnsemble.h:375
int type
Definition: TGX11.cxx:120
static RooMathCoreReg dummy
void * GetNextChild(void *prevchild, const char *childname=0)
XML helpers.
Definition: Tools.cxx:1174
void MakeLinearTerms()
Make the linear terms as in eq 25, ref 2 For this the b and (1-b) quantiles are needed.
std::vector< Double_t > fLinImportance
Definition: RuleEnsemble.h:353
static constexpr double s
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t fAverageRuleSigma
Definition: RuleEnsemble.h:357
void MakeRules(const std::vector< const TMVA::DecisionTree *> &forest)
Makes rules from the given decision tree.
std::vector< Double_t > fRulePBS
Definition: RuleEnsemble.h:362
void CalcRuleSupport()
calculate the support for all rules
Bool_t UseBoost() const
Definition: MethodRuleFit.h:84
ostringstream derivative to redirect and format output
Definition: MsgLogger.h:59
const std::vector< const TMVA::Event * > * fRuleMapEvents
Definition: RuleEnsemble.h:381
const MethodBase * GetMethodBase() const
Get a pointer to the original MethodRuleFit.
std::vector< Double_t > fRulePSS
Definition: RuleEnsemble.h:360
void MakeModel()
create model
Double_t GetOffset() const
Definition: RuleEnsemble.h:265
Node for the BinarySearch or Decision Trees.
Definition: Node.h:56
void SetMsgType(EMsgType t)
Definition: Rule.cxx:154
Bool_t DoRules() const
Definition: RuleEnsemble.h:258
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
const MethodRuleFit * GetMethodRuleFit() const
Get a pointer to the original MethodRuleFit.
void ResetCoefficients()
reset all rule coefficients
void GetCoefficients(std::vector< Double_t > &v)
Retrieve all rule coefficients.
Double_t EvalEvent() const
Definition: RuleEnsemble.h:416
Short_t GetSelector() const
Definition: first.py:1
const MethodRuleFit * GetMethodRuleFit() const
Definition: RuleFit.h:148
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Double_t CalcLinImportance()
calculate the linear importance for each rule
UInt_t GetNTreeSample() const
Definition: RuleFit.h:130
const Bool_t kTRUE
Definition: RtypesCore.h:87
std::vector< Double_t > fRulePTag
Definition: RuleEnsemble.h:364
Double_t fImportanceCut
Definition: RuleEnsemble.h:342
const Event * GetTrainingEvent(UInt_t i) const
get the training event from the rule fitter
Bool_t Equal(const Rule &other, Bool_t useCutValue, Double_t maxdist) const
Compare two rules.
Definition: Rule.cxx:170
std::vector< Double_t > fRuleVarFrac
Definition: RuleEnsemble.h:359
std::vector< std::vector< UInt_t > > fRuleMap
Definition: RuleEnsemble.h:378
const Event * fEvent
Definition: RuleEnsemble.h:372
void ReadRaw(std::istream &istr)
read rule ensemble from stream