// @(#)root/tmva $Id$ // Author: Andreas Hoecker, Peter Speckmayer, Joerg Stelzer /********************************************************************************** * Project: TMVA - a Root-integrated toolkit for multivariate data analysis * * Package: TMVA * * Class : MethodFDA * * Web : http://tmva.sourceforge.net * * * * Description: * * Implementation * * * * Authors (alphabetical): * * Andreas Hoecker - CERN, Switzerland * * Peter Speckmayer - CERN, Switzerland * * Joerg Stelzer - DESY, Germany * * Maciej Kruk - IFJ PAN & AGH, Poland * * * * Copyright (c) 2005-2006: * * CERN, Switzerland * * MPI-K Heidelberg, Germany * * * * Redistribution and use in source and binary forms, with or without * * modification, are permitted according to the terms listed in LICENSE * * (http://tmva.sourceforge.net/LICENSE) * **********************************************************************************/ //_______________________________________________________________________ // // Function discriminant analysis (FDA). This simple classifier // // fits any user-defined TFormula (via option configuration string) to // // the training data by requiring a formula response of 1 (0) to signal // // (background) events. The parameter fitting is done via the abstract // // class FitterBase, featuring Monte Carlo sampling, Genetic // // Algorithm, Simulated Annealing, MINUIT and combinations of these. // // // // Can compute regression value for one dimensional output // //_______________________________________________________________________ #include "Riostream.h" #include "TList.h" #include "TFormula.h" #include "TString.h" #include "TObjString.h" #include "TRandom3.h" #include "TMath.h" #include #include #include #include #include "TMVA/ClassifierFactory.h" #include "TMVA/MethodFDA.h" #include "TMVA/Tools.h" #include "TMVA/Interval.h" #include "TMVA/Timer.h" #include "TMVA/GeneticFitter.h" #include "TMVA/SimulatedAnnealingFitter.h" #include "TMVA/MinuitFitter.h" #include "TMVA/MCFitter.h" #include "TMVA/Config.h" using std::stringstream; REGISTER_METHOD(FDA) ClassImp(TMVA::MethodFDA) //_______________________________________________________________________ TMVA::MethodFDA::MethodFDA( const TString& jobName, const TString& methodTitle, DataSetInfo& theData, const TString& theOption, TDirectory* theTargetDir ) : MethodBase( jobName, Types::kFDA, methodTitle, theData, theOption, theTargetDir ), IFitterTarget (), fFormula ( 0 ), fNPars ( 0 ), fFitter ( 0 ), fConvergerFitter( 0 ), fSumOfWeightsSig( 0 ), fSumOfWeightsBkg( 0 ), fSumOfWeights ( 0 ), fOutputDimensions( 0 ) { // standard constructor } //_______________________________________________________________________ TMVA::MethodFDA::MethodFDA( DataSetInfo& theData, const TString& theWeightFile, TDirectory* theTargetDir ) : MethodBase( Types::kFDA, theData, theWeightFile, theTargetDir ), IFitterTarget (), fFormula ( 0 ), fNPars ( 0 ), fFitter ( 0 ), fConvergerFitter( 0 ), fSumOfWeightsSig( 0 ), fSumOfWeightsBkg( 0 ), fSumOfWeights ( 0 ), fOutputDimensions( 0 ) { // constructor from weight file } //_______________________________________________________________________ void TMVA::MethodFDA::Init( void ) { // default initialisation fNPars = 0; fBestPars.clear(); fSumOfWeights = 0; fSumOfWeightsSig = 0; fSumOfWeightsBkg = 0; fFormulaStringP = ""; fParRangeStringP = ""; fFormulaStringT = ""; fParRangeStringT = ""; fFitMethod = ""; fConverger = ""; if( DoMulticlass() ) if (fMulticlassReturnVal == NULL) fMulticlassReturnVal = new std::vector(); } //_______________________________________________________________________ void TMVA::MethodFDA::DeclareOptions() { // define the options (their key words) that can be set in the option string // // format of function string: // "x0*(0)+((1)/x1)**(2)..." // where "[i]" are the parameters, and "xi" the input variables // // format of parameter string: // "(-1.2,3.4);(-2.3,4.55);..." // where the numbers in "(a,b)" correspond to the a=min, b=max parameter ranges; // each parameter defined in the function string must have a corresponding range // DeclareOptionRef( fFormulaStringP = "(0)", "Formula", "The discrimination formula" ); DeclareOptionRef( fParRangeStringP = "()", "ParRanges", "Parameter ranges" ); // fitter DeclareOptionRef( fFitMethod = "MINUIT", "FitMethod", "Optimisation Method"); AddPreDefVal(TString("MC")); AddPreDefVal(TString("GA")); AddPreDefVal(TString("SA")); AddPreDefVal(TString("MINUIT")); DeclareOptionRef( fConverger = "None", "Converger", "FitMethod uses Converger to improve result"); AddPreDefVal(TString("None")); AddPreDefVal(TString("MINUIT")); } //_______________________________________________________________________ void TMVA::MethodFDA::CreateFormula() { // translate formula string into TFormula, and parameter string into par ranges // process transient strings fFormulaStringT = fFormulaStringP; // intepret formula string // replace the parameters "(i)" by the TFormula style "[i]" for (UInt_t ipar=0; ipar Formula contains expression: \"" << Form("(%i)",ipar) << "\", " << "which cannot be attributed to a parameter; " << "it may be that the number of variable ranges given via \"ParRanges\" " << "does not match the number of parameters in the formula expression, please verify!" << Endl; } // write the variables "xi" as additional parameters "[npar+i]" for (Int_t ivar=GetNvar()-1; ivar >= 0; ivar--) { fFormulaStringT.ReplaceAll( Form("x%i",ivar), Form("[%i]",ivar+fNPars) ); } // sanity check, there should be no "xi", with 'i' a number anymore for (UInt_t ivar=GetNvar(); ivar<1000; ivar++) { if (fFormulaStringT.Contains( Form("x%i",ivar) )) Log() << kFATAL << " Formula contains expression: \"" << Form("x%i",ivar) << "\", " << "which cannot be attributed to an input variable" << Endl; } Log() << "User-defined formula string : \"" << fFormulaStringP << "\"" << Endl; Log() << "TFormula-compatible formula string: \"" << fFormulaStringT << "\"" << Endl; Log() << "Creating and compiling formula" << Endl; // create TF1 if (fFormula) delete fFormula; fFormula = new TFormula( "FDA_Formula", fFormulaStringT ); #if ROOT_VERSION_CODE >= ROOT_VERSION(5,2,0) fFormula->Optimize(); #endif // is formula correct ? if (fFormula->Compile() != 0) Log() << kFATAL << " Formula expression could not be properly compiled" << Endl; // other sanity checks if (fFormula->GetNpar() > (Int_t)(fNPars + GetNvar())) Log() << kFATAL << " Dubious number of parameters in formula expression: " << fFormula->GetNpar() << " - compared to maximum allowed: " << fNPars + GetNvar() << Endl; } //_______________________________________________________________________ void TMVA::MethodFDA::ProcessOptions() { // the option string is decoded, for availabel options see "DeclareOptions" // process transient strings fParRangeStringT = fParRangeStringP; // interpret parameter string fParRangeStringT.ReplaceAll( " ", "" ); fNPars = fParRangeStringT.CountChar( ')' ); TList* parList = gTools().ParseFormatLine( fParRangeStringT, ";" ); if ((UInt_t)parList->GetSize() != fNPars) { Log() << kFATAL << " Mismatch in parameter string: " << "the number of parameters: " << fNPars << " != ranges defined: " << parList->GetSize() << "; the format of the \"ParRanges\" string " << "must be: \"(-1.2,3.4);(-2.3,4.55);...\", " << "where the numbers in \"(a,b)\" correspond to the a=min, b=max parameter ranges; " << "each parameter defined in the function string must have a corresponding rang." << Endl; } fParRange.resize( fNPars ); for (UInt_t ipar=0; iparAt(ipar))->GetString(); Ssiz_t istr = str.First( ',' ); TString pminS(str(1,istr-1)); TString pmaxS(str(istr+1,str.Length()-2-istr)); stringstream stmin; Float_t pmin=0; stmin << pminS.Data(); stmin >> pmin; stringstream stmax; Float_t pmax=0; stmax << pmaxS.Data(); stmax >> pmax; // sanity check if (TMath::Abs(pmax-pmin) < 1.e-30) pmax = pmin; if (pmin > pmax) Log() << kFATAL << " max > min in interval for parameter: [" << ipar << "] : [" << pmin << ", " << pmax << "] " << Endl; Log() << kINFO << "Create parameter interval for parameter " << ipar << " : [" << pmin << "," << pmax << "]" << Endl; fParRange[ipar] = new Interval( pmin, pmax ); } delete parList; // create formula CreateFormula(); // copy parameter ranges for each output dimension ================== fOutputDimensions = 1; if( DoRegression() ) fOutputDimensions = DataInfo().GetNTargets(); if( DoMulticlass() ) fOutputDimensions = DataInfo().GetNClasses(); for( Int_t dim = 1; dim < fOutputDimensions; ++dim ){ for( UInt_t par = 0; par < fNPars; ++par ){ fParRange.push_back( fParRange.at(par) ); } } // ==================== // create minimiser fConvergerFitter = (IFitterTarget*)this; if (fConverger == "MINUIT") { fConvergerFitter = new MinuitFitter( *this, Form("%s_Converger_Minuit", GetName()), fParRange, GetOptions() ); SetOptions(dynamic_cast(fConvergerFitter)->GetOptions()); } if(fFitMethod == "MC") fFitter = new MCFitter( *fConvergerFitter, Form("%s_Fitter_MC", GetName()), fParRange, GetOptions() ); else if (fFitMethod == "GA") fFitter = new GeneticFitter( *fConvergerFitter, Form("%s_Fitter_GA", GetName()), fParRange, GetOptions() ); else if (fFitMethod == "SA") fFitter = new SimulatedAnnealingFitter( *fConvergerFitter, Form("%s_Fitter_SA", GetName()), fParRange, GetOptions() ); else if (fFitMethod == "MINUIT") fFitter = new MinuitFitter( *fConvergerFitter, Form("%s_Fitter_Minuit", GetName()), fParRange, GetOptions() ); else { Log() << kFATAL << " Do not understand fit method:" << fFitMethod << Endl; } fFitter->CheckForUnusedOptions(); } //_______________________________________________________________________ TMVA::MethodFDA::~MethodFDA( void ) { // destructor ClearAll(); } //_______________________________________________________________________ Bool_t TMVA::MethodFDA::HasAnalysisType( Types::EAnalysisType type, UInt_t numberClasses, UInt_t /*numberTargets*/ ) { // FDA can handle classification with 2 classes and regression with one regression-target if (type == Types::kClassification && numberClasses == 2) return kTRUE; if (type == Types::kMulticlass ) return kTRUE; if (type == Types::kRegression ) return kTRUE; return kFALSE; } //_______________________________________________________________________ void TMVA::MethodFDA::ClearAll( void ) { // delete and clear all class members // if there is more than one output dimension, the paramater ranges are the same again (object has been copied). // hence, ... erase the copied pointers to assure, that they are deleted only once. // fParRange.erase( fParRange.begin()+(fNPars), fParRange.end() ); for (UInt_t ipar=0; iparGetWeight(); if (!DoRegression()) { if (DataInfo().IsSignal(ev)) { fSumOfWeightsSig += w; } else { fSumOfWeightsBkg += w; } } fSumOfWeights += w; } // sanity check if (!DoRegression()) { if (fSumOfWeightsSig <= 0 || fSumOfWeightsBkg <= 0) { Log() << kFATAL << " Troubles in sum of weights: " << fSumOfWeightsSig << " (S) : " << fSumOfWeightsBkg << " (B)" << Endl; } } else if (fSumOfWeights <= 0) { Log() << kFATAL << " Troubles in sum of weights: " << fSumOfWeights << Endl; } // starting values (not used by all fitters) fBestPars.clear(); for (std::vector::const_iterator parIt = fParRange.begin(); parIt != fParRange.end(); parIt++) { fBestPars.push_back( (*parIt)->GetMean() ); } // execute the fit Double_t estimator = fFitter->Run( fBestPars ); // print results PrintResults( fFitMethod, fBestPars, estimator ); delete fFitter; fFitter = 0; if (fConvergerFitter!=0 && fConvergerFitter!=(IFitterTarget*)this) { delete fConvergerFitter; fConvergerFitter = 0; } } //_______________________________________________________________________ void TMVA::MethodFDA::PrintResults( const TString& fitter, std::vector& pars, const Double_t estimator ) const { // display fit parameters // check maximum length of variable name Log() << kINFO; Log() << "Results for parameter fit using \"" << fitter << "\" fitter:" << Endl; std::vector parNames; for (UInt_t ipar=0; ipar& pars ) { // compute estimator for given parameter set (to be minimised) // const Double_t sumOfWeights[] = { fSumOfWeightsSig, fSumOfWeightsBkg, fSumOfWeights }; const Double_t sumOfWeights[] = { fSumOfWeightsBkg, fSumOfWeightsSig, fSumOfWeights }; Double_t estimator[] = { 0, 0, 0 }; Double_t result, deviation; Double_t desired = 0.0; // calculate the deviation from the desired value if( DoRegression() ){ for (UInt_t ievt=0; ievtGetTarget( dim ); result = InterpretFormula( ev, pars.begin(), pars.end() ); deviation = TMath::Power(result - desired, 2); estimator[2] += deviation * ev->GetWeight(); } } estimator[2] /= sumOfWeights[2]; // return value is sum over normalised signal and background contributions return estimator[2]; }else if( DoMulticlass() ){ for (UInt_t ievt=0; ievtat(dim); Double_t t = (ev->GetClass() == static_cast(dim) ? 1.0 : 0.0 ); crossEntropy += t*log(y); } estimator[2] += ev->GetWeight()*crossEntropy; } estimator[2] /= sumOfWeights[2]; // return value is sum over normalised signal and background contributions return estimator[2]; }else{ for (UInt_t ievt=0; ievtGetWeight(); } estimator[0] /= sumOfWeights[0]; estimator[1] /= sumOfWeights[1]; // return value is sum over normalised signal and background contributions return estimator[0] + estimator[1]; } } //_______________________________________________________________________ Double_t TMVA::MethodFDA::InterpretFormula( const Event* event, std::vector::iterator parBegin, std::vector::iterator parEnd ) { // formula interpretation Int_t ipar = 0; // std::cout << "pars "; for( std::vector::iterator it = parBegin; it != parEnd; ++it ){ // std::cout << " i" << ipar << " val" << (*it); fFormula->SetParameter( ipar, (*it) ); ++ipar; } for (UInt_t ivar=0; ivarSetParameter( ivar+ipar, event->GetValue(ivar) ); Double_t result = fFormula->Eval( 0 ); // std::cout << " result " << result << std::endl; return result; } //_______________________________________________________________________ Double_t TMVA::MethodFDA::GetMvaValue( Double_t* err, Double_t* errUpper ) { // returns MVA value for given event const Event* ev = GetEvent(); // cannot determine error NoErrorCalc(err, errUpper); return InterpretFormula( ev, fBestPars.begin(), fBestPars.end() ); } //_______________________________________________________________________ const std::vector& TMVA::MethodFDA::GetRegressionValues() { if (fRegressionReturnVal == NULL) fRegressionReturnVal = new std::vector(); fRegressionReturnVal->clear(); const Event* ev = GetEvent(); Event* evT = new Event(*ev); for( Int_t dim = 0; dim < fOutputDimensions; ++dim ){ Int_t offset = dim*fNPars; evT->SetTarget(dim,InterpretFormula( ev, fBestPars.begin()+offset, fBestPars.begin()+offset+fNPars ) ); } const Event* evT2 = GetTransformationHandler().InverseTransform( evT ); fRegressionReturnVal->push_back(evT2->GetTarget(0)); delete evT; return (*fRegressionReturnVal); } //_______________________________________________________________________ const std::vector& TMVA::MethodFDA::GetMulticlassValues() { if (fMulticlassReturnVal == NULL) fMulticlassReturnVal = new std::vector(); fMulticlassReturnVal->clear(); std::vector temp; // returns MVA value for given event const TMVA::Event* evt = GetEvent(); CalculateMulticlassValues( evt, fBestPars, temp ); UInt_t nClasses = DataInfo().GetNClasses(); for(UInt_t iClass=0; iClass& parameters, std::vector& values) { // calculate the values for multiclass values.clear(); // std::copy( parameters.begin(), parameters.end(), std::ostream_iterator( std::cout, " " ) ); // std::cout << std::endl; // char inp; // std::cin >> inp; Double_t sum=0; for( Int_t dim = 0; dim < fOutputDimensions; ++dim ){ // check for all other dimensions (=classes) Int_t offset = dim*fNPars; Double_t value = InterpretFormula( evt, parameters.begin()+offset, parameters.begin()+offset+fNPars ); // std::cout << "dim : " << dim << " value " << value << " offset " << offset << std::endl; values.push_back( value ); sum += value; } // // normalize to sum of value (commented out, .. have to think of how to treat negative classifier values) // std::transform( fMulticlassReturnVal.begin(), fMulticlassReturnVal.end(), fMulticlassReturnVal.begin(), bind2nd( std::divides(), sum) ); } //_______________________________________________________________________ void TMVA::MethodFDA::ReadWeightsFromStream( std::istream& istr ) { // read back the training results from a file (stream) // retrieve best function parameters // coverity[tainted_data_argument] istr >> fNPars; fBestPars.clear(); fBestPars.resize( fNPars ); for (UInt_t ipar=0; ipar> fBestPars[ipar]; } //_______________________________________________________________________ void TMVA::MethodFDA::AddWeightsXMLTo( void* parent ) const { // create XML description for LD classification and regression // (for arbitrary number of output classes/targets) void* wght = gTools().AddChild(parent, "Weights"); gTools().AddAttr( wght, "NPars", fNPars ); gTools().AddAttr( wght, "NDim", fOutputDimensions ); for (UInt_t ipar=0; ipar= fNPars*fOutputDimensions) Log() << kFATAL << " index out of range: " << ipar << " >= " << fNPars << Endl; fBestPars[ipar] = par; ch = gTools().GetNextChild(ch); } // read formula gTools().ReadAttr( wghtnode, "Formula", fFormulaStringP ); // create the TFormula CreateFormula(); } //_______________________________________________________________________ void TMVA::MethodFDA::MakeClassSpecific( std::ostream& fout, const TString& className ) const { // write FDA-specific classifier response fout << " double fParameter[" << fNPars << "];" << std::endl; fout << "};" << std::endl; fout << "" << std::endl; fout << "inline void " << className << "::Initialize() " << std::endl; fout << "{" << std::endl; for(UInt_t ipar=0; ipar& inputValues ) const" << std::endl; fout << "{" << std::endl; fout << " // interpret the formula" << std::endl; // replace parameters TString str = fFormulaStringT; for (UInt_t ipar=0; ipar" << "http://tmva.sourceforge.net/docu/TMVAUsersGuide.pdf" << Endl; } else Log() << "http://tmva.sourceforge.net/docu/TMVAUsersGuide.pdf" << Endl; Log() << Endl; Log() << gTools().Color("bold") << "--- Performance optimisation:" << gTools().Color("reset") << Endl; Log() << Endl; Log() << "The FDA performance depends on the complexity and fidelity of the" << Endl; Log() << "user-defined discriminator function. As a general rule, it should" << Endl; Log() << "be able to reproduce the discrimination power of any linear" << Endl; Log() << "discriminant analysis. To reach into the nonlinear domain, it is" << Endl; Log() << "useful to inspect the correlation profiles of the input variables," << Endl; Log() << "and add quadratic and higher polynomial terms between variables as" << Endl; Log() << "necessary. Comparison with more involved nonlinear classifiers can" << Endl; Log() << "be used as a guide." << Endl; Log() << Endl; Log() << gTools().Color("bold") << "--- Performance tuning via configuration options:" << gTools().Color("reset") << Endl; Log() << Endl; Log() << "Depending on the function used, the choice of \"FitMethod\" is" << Endl; Log() << "crucial for getting valuable solutions with FDA. As a guideline it" << Endl; Log() << "is recommended to start with \"FitMethod=MINUIT\". When more complex" << Endl; Log() << "functions are used where MINUIT does not converge to reasonable" << Endl; Log() << "results, the user should switch to non-gradient FitMethods such" << Endl; Log() << "as GeneticAlgorithm (GA) or Monte Carlo (MC). It might prove to be" << Endl; Log() << "useful to combine GA (or MC) with MINUIT by setting the option" << Endl; Log() << "\"Converger=MINUIT\". GA (MC) will then set the starting parameters" << Endl; Log() << "for MINUIT such that the basic quality of GA (MC) of finding global" << Endl; Log() << "minima is combined with the efficacy of MINUIT of finding local" << Endl; Log() << "minima." << Endl; }