00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "Riostream.h"
00040 #include "TMath.h"
00041 #include "TFile.h"
00042
00043 #include "TMVA/ClassifierFactory.h"
00044 #ifndef ROOT_TMVA_MethodSVM
00045 #include "TMVA/MethodSVM.h"
00046 #endif
00047 #ifndef ROOT_TMVA_Tools
00048 #include "TMVA/Tools.h"
00049 #endif
00050 #ifndef ROOT_TMVA_Timer
00051 #include "TMVA/Timer.h"
00052 #endif
00053
00054 #ifndef ROOT_TMVA_SVWorkingSet
00055 #include "TMVA/SVWorkingSet.h"
00056 #endif
00057
00058 #ifndef ROOT_TMVA_SVEvent
00059 #include "TMVA/SVEvent.h"
00060 #endif
00061
00062 #ifndef ROOT_TMVA_SVKernelFunction
00063 #include "TMVA/SVKernelFunction.h"
00064 #endif
00065
00066 #include <string>
00067
00068 const Int_t basketsize__ = 1280000;
00069 REGISTER_METHOD(SVM)
00070
00071 ClassImp(TMVA::MethodSVM)
00072
00073
00074 TMVA::MethodSVM::MethodSVM( const TString& jobName, const TString& methodTitle, DataSetInfo& theData,
00075 const TString& theOption, TDirectory* theTargetDir )
00076 : MethodBase( jobName, Types::kSVM, methodTitle, theData, theOption, theTargetDir )
00077 , fCost(0)
00078 , fTolerance(0)
00079 , fMaxIter(0)
00080 , fNSubSets(0)
00081 , fBparm(0)
00082 , fGamma(0)
00083 , fWgSet(0)
00084 , fInputData(0)
00085 , fSupportVectors(0)
00086 , fSVKernelFunction(0)
00087 , fMinVars(0)
00088 , fMaxVars(0)
00089 , fDoubleSigmaSquared(0)
00090 , fOrder(0)
00091 , fTheta(0)
00092 , fKappa(0)
00093 {
00094
00095 }
00096
00097
00098 TMVA::MethodSVM::MethodSVM( DataSetInfo& theData, const TString& theWeightFile, TDirectory* theTargetDir )
00099 : MethodBase( Types::kSVM, theData, theWeightFile, theTargetDir )
00100 , fCost(0)
00101 , fTolerance(0)
00102 , fMaxIter(0)
00103 , fNSubSets(0)
00104 , fBparm(0)
00105 , fGamma(0)
00106 , fWgSet(0)
00107 , fInputData(0)
00108 , fSupportVectors(0)
00109 , fSVKernelFunction(0)
00110 , fMinVars(0)
00111 , fMaxVars(0)
00112 , fDoubleSigmaSquared(0)
00113 , fOrder(0)
00114 , fTheta(0)
00115 , fKappa(0)
00116 {
00117
00118 }
00119
00120
00121 TMVA::MethodSVM::~MethodSVM()
00122 {
00123
00124 if (fInputData !=0) { delete fInputData; fInputData=0; }
00125 if (fSupportVectors !=0 ) { delete fSupportVectors; fSupportVectors = 0; }
00126 if (fWgSet !=0) { delete fWgSet; fWgSet=0; }
00127 if (fSVKernelFunction !=0 ) { delete fSVKernelFunction; fSVKernelFunction = 0; }
00128 }
00129
00130
00131 Bool_t TMVA::MethodSVM::HasAnalysisType( Types::EAnalysisType type, UInt_t numberClasses, UInt_t numberTargets )
00132 {
00133
00134 if (type == Types::kClassification && numberClasses == 2) return kTRUE;
00135 if (type == Types::kRegression && numberTargets == 1) return kTRUE;
00136 return kFALSE;
00137 }
00138
00139
00140 void TMVA::MethodSVM::Init()
00141 {
00142
00143
00144
00145 SetNormalised( kTRUE );
00146
00147 fInputData = new std::vector<TMVA::SVEvent*>(Data()->GetNEvents());
00148 fSupportVectors = new std::vector<TMVA::SVEvent*>(0);
00149 }
00150
00151
00152 void TMVA::MethodSVM::DeclareOptions()
00153 {
00154
00155 DeclareOptionRef( fCost, "C", "Cost parameter" );
00156 if (DoRegression()) {
00157 fCost = 0.002;
00158 }else{
00159 fCost = 1.;
00160 }
00161 DeclareOptionRef( fTolerance = 0.01, "Tol", "Tolerance parameter" );
00162 DeclareOptionRef( fMaxIter = 1000, "MaxIter", "Maximum number of training loops" );
00163 DeclareOptionRef( fNSubSets = 1, "NSubSets", "Number of training subsets" );
00164
00165
00166 DeclareOptionRef( fGamma = 1., "Gamma", "RBF kernel parameter: Gamma");
00167 }
00168
00169
00170 void TMVA::MethodSVM::DeclareCompatibilityOptions()
00171 {
00172 MethodBase::DeclareCompatibilityOptions();
00173 DeclareOptionRef( fTheKernel = "Gauss", "Kernel", "Uses kernel function");
00174
00175 DeclareOptionRef( fDoubleSigmaSquared = 2., "Sigma", "Kernel parameter: sigma");
00176
00177 DeclareOptionRef( fOrder = 3, "Order", "Polynomial Kernel parameter: polynomial order");
00178
00179 DeclareOptionRef( fTheta = 1., "Theta", "Sigmoid Kernel parameter: theta");
00180 DeclareOptionRef( fKappa = 1., "Kappa", "Sigmoid Kernel parameter: kappa");
00181 }
00182
00183
00184 void TMVA::MethodSVM::ProcessOptions()
00185 {
00186
00187 if (IgnoreEventsWithNegWeightsInTraining()) {
00188 Log() << kFATAL << "Mechanism to ignore events with negative weights in training not yet available for method: "
00189 << GetMethodTypeName()
00190 << " --> please remove \"IgnoreNegWeightsInTraining\" option from booking string."
00191 << Endl;
00192 }
00193 }
00194
00195
00196 void TMVA::MethodSVM::Train()
00197 {
00198
00199 Data()->SetCurrentType(Types::kTraining);
00200
00201 for (Int_t ievt=0; ievt<Data()->GetNEvents(); ievt++){
00202 Log() << kDEBUG << "Create event vector"<< Endl;
00203 fInputData->at(ievt) = new SVEvent(GetEvent(ievt), fCost);
00204 }
00205
00206 fSVKernelFunction = new SVKernelFunction(fGamma);
00207
00208 Log()<< kINFO << "Building SVM Working Set..."<< Endl;
00209 Timer bldwstime( GetName());
00210 fWgSet = new SVWorkingSet( fInputData, fSVKernelFunction,fTolerance, DoRegression() );
00211 Log() << kINFO <<"Elapsed time for Working Set build: "<< bldwstime.GetElapsedTime()<<Endl;
00212
00213
00214 Timer timer( GetName() );
00215 Log() << kINFO << "Sorry, no computing time forecast available for SVM, please wait ..." << Endl;
00216
00217 fWgSet->Train(fMaxIter);
00218
00219 Log() << kINFO << "Elapsed time: " << timer.GetElapsedTime()
00220 << " " << Endl;
00221
00222 fBparm = fWgSet->GetBpar();
00223 fSupportVectors = fWgSet->GetSupportVectors();
00224 }
00225
00226
00227 void TMVA::MethodSVM::AddWeightsXMLTo( void* parent ) const
00228 {
00229
00230 void* wght = gTools().AddChild(parent, "Weights");
00231 gTools().AddAttr(wght,"fBparm",fBparm);
00232 gTools().AddAttr(wght,"fGamma",fGamma);
00233 gTools().AddAttr(wght,"NSupVec",fSupportVectors->size());
00234
00235 for (std::vector<TMVA::SVEvent*>::iterator veciter=fSupportVectors->begin();
00236 veciter!=fSupportVectors->end() ; ++veciter ) {
00237 TVectorD temp(GetNvar()+4);
00238 temp[0] = (*veciter)->GetNs();
00239 temp[1] = (*veciter)->GetTypeFlag();
00240 temp[2] = (*veciter)->GetAlpha();
00241 temp[3] = (*veciter)->GetAlpha_p();
00242 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00243 temp[ivar+4] = (*(*veciter)->GetDataVector())[ivar];
00244 gTools().WriteTVectorDToXML(wght,"SupportVector",&temp);
00245 }
00246
00247 void* maxnode = gTools().AddChild(wght, "Maxima");
00248 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00249 gTools().AddAttr(maxnode, "Var"+gTools().StringFromInt(ivar), GetXmax(ivar));
00250 void* minnode = gTools().AddChild(wght, "Minima");
00251 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00252 gTools().AddAttr(minnode, "Var"+gTools().StringFromInt(ivar), GetXmin(ivar));
00253 }
00254
00255
00256 void TMVA::MethodSVM::ReadWeightsFromXML( void* wghtnode )
00257 {
00258 gTools().ReadAttr( wghtnode, "fBparm",fBparm );
00259 gTools().ReadAttr( wghtnode, "fGamma",fGamma);
00260 UInt_t fNsupv=0;
00261 gTools().ReadAttr( wghtnode, "NSupVec",fNsupv );
00262
00263 Float_t alpha=0.;
00264 Float_t alpha_p = 0.;
00265
00266 Int_t typeFlag=-1;
00267 UInt_t ns = 0;
00268 std::vector<Float_t>* svector = new std::vector<Float_t>(GetNvar());
00269
00270 if (fMaxVars!=0) delete fMaxVars;
00271 fMaxVars = new TVectorD( GetNvar() );
00272 if (fMinVars!=0) delete fMinVars;
00273 fMinVars = new TVectorD( GetNvar() );
00274 if (fSupportVectors!=0) {
00275 for (vector< SVEvent* >::iterator it = fSupportVectors->begin(); it!=fSupportVectors->end(); ++it)
00276 delete *it;
00277 delete fSupportVectors;
00278 }
00279 fSupportVectors = new std::vector<TMVA::SVEvent*>(0);
00280 void* supportvectornode = gTools().GetChild(wghtnode);
00281 for (UInt_t ievt = 0; ievt < fNsupv; ievt++) {
00282 TVectorD temp(GetNvar()+4);
00283 gTools().ReadTVectorDFromXML(supportvectornode,"SupportVector",&temp);
00284 ns=(UInt_t)temp[0];
00285 typeFlag=(int)temp[1];
00286 alpha=temp[2];
00287 alpha_p=temp[3];
00288 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00289 (*svector)[ivar]=temp[ivar+4];
00290
00291 fSupportVectors->push_back(new SVEvent(svector,alpha,alpha_p,typeFlag));
00292 supportvectornode = gTools().GetNextChild(supportvectornode);
00293 }
00294
00295 void* maxminnode = supportvectornode;
00296 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00297 gTools().ReadAttr( maxminnode,"Var"+gTools().StringFromInt(ivar),(*fMaxVars)[ivar]);
00298 maxminnode = gTools().GetNextChild(maxminnode);
00299 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00300 gTools().ReadAttr( maxminnode,"Var"+gTools().StringFromInt(ivar),(*fMinVars)[ivar]);
00301 if (fSVKernelFunction!=0) delete fSVKernelFunction;
00302 fSVKernelFunction = new SVKernelFunction(fGamma);
00303 delete svector;
00304 }
00305
00306
00307 void TMVA::MethodSVM::WriteWeightsToStream( TFile& ) const
00308 {
00309
00310
00311 }
00312
00313
00314 void TMVA::MethodSVM::ReadWeightsFromStream( istream& istr )
00315 {
00316 if (fSupportVectors !=0) { delete fSupportVectors; fSupportVectors = 0;}
00317 fSupportVectors = new std::vector<TMVA::SVEvent*>(0);
00318
00319
00320 istr >> fBparm;
00321
00322 UInt_t fNsupv;
00323
00324 istr >> fNsupv;
00325 fSupportVectors->reserve(fNsupv);
00326
00327 Float_t typeTalpha=0.;
00328 Float_t alpha=0.;
00329 Int_t typeFlag=-1;
00330 UInt_t ns = 0;
00331 std::vector<Float_t>* svector = new std::vector<Float_t>(GetNvar());
00332
00333 fMaxVars = new TVectorD( GetNvar() );
00334 fMinVars = new TVectorD( GetNvar() );
00335
00336 for (UInt_t ievt = 0; ievt < fNsupv; ievt++) {
00337 istr>>ns;
00338 istr>>typeTalpha;
00339 typeFlag = typeTalpha<0?-1:1;
00340 alpha = typeTalpha<0?-typeTalpha:typeTalpha;
00341 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00342 istr>>svector->at(ivar);
00343
00344 fSupportVectors->push_back(new SVEvent(svector,alpha,typeFlag,ns));
00345 }
00346
00347 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00348 istr >> (*fMaxVars)[ivar];
00349
00350 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++)
00351 istr >> (*fMinVars)[ivar];
00352
00353 delete fSVKernelFunction;
00354 if (fTheKernel == "Gauss" ) {
00355 fSVKernelFunction = new SVKernelFunction(1/fDoubleSigmaSquared);
00356 } else {
00357 SVKernelFunction::EKernelType k = SVKernelFunction::kLinear;
00358 if(fTheKernel == "Linear") k = SVKernelFunction::kLinear;
00359 else if (fTheKernel == "Polynomial") k = SVKernelFunction::kPolynomial;
00360 else if (fTheKernel == "Sigmoid" ) k = SVKernelFunction::kSigmoidal;
00361 else {
00362 Log() << kFATAL <<"Unknown kernel function found in weight file!" << Endl;
00363 }
00364 fSVKernelFunction = new SVKernelFunction();
00365 fSVKernelFunction->setCompatibilityParams(k, fOrder, fTheta, fKappa);
00366 }
00367 delete svector;
00368 }
00369
00370
00371 void TMVA::MethodSVM::ReadWeightsFromStream( TFile& )
00372 {
00373
00374 }
00375
00376
00377 Double_t TMVA::MethodSVM::GetMvaValue( Double_t* err, Double_t* errUpper )
00378 {
00379
00380 Double_t myMVA = 0;
00381
00382
00383 SVEvent* ev = new SVEvent( GetEvent(),0. );
00384
00385 for (UInt_t ievt = 0; ievt < fSupportVectors->size() ; ievt++) {
00386 myMVA += ( fSupportVectors->at(ievt)->GetAlpha()
00387 * fSupportVectors->at(ievt)->GetTypeFlag()
00388 * fSVKernelFunction->Evaluate( fSupportVectors->at(ievt), ev ) );
00389 }
00390
00391 delete ev;
00392
00393 myMVA -= fBparm;
00394
00395
00396 NoErrorCalc(err, errUpper);
00397
00398
00399 return 1.0/(1.0 + TMath::Exp(myMVA));
00400 }
00401
00402 const std::vector<Float_t>& TMVA::MethodSVM::GetRegressionValues()
00403 {
00404
00405 if( fRegressionReturnVal == NULL )
00406 fRegressionReturnVal = new std::vector<Float_t>();
00407 fRegressionReturnVal->clear();
00408
00409 Double_t myMVA = 0;
00410
00411 const Event *baseev = GetEvent();
00412 SVEvent* ev = new SVEvent( baseev,0. );
00413
00414 for (UInt_t ievt = 0; ievt < fSupportVectors->size() ; ievt++) {
00415 myMVA += ( fSupportVectors->at(ievt)->GetDeltaAlpha()
00416 *fSVKernelFunction->Evaluate( fSupportVectors->at(ievt), ev ) );
00417 }
00418 myMVA += fBparm;
00419 Event * evT = new Event(*baseev);
00420 evT->SetTarget(0,myMVA);
00421
00422 const Event* evT2 = GetTransformationHandler().InverseTransform( evT );
00423
00424 fRegressionReturnVal->push_back(evT2->GetTarget(0));
00425
00426 delete evT;
00427
00428 delete ev;
00429
00430 return *fRegressionReturnVal;
00431 }
00432
00433
00434 void TMVA::MethodSVM::MakeClassSpecific( std::ostream& fout, const TString& className ) const
00435 {
00436
00437 const int fNsupv = fSupportVectors->size();
00438 fout << " // not implemented for class: \"" << className << "\"" << endl;
00439 fout << " float fBparameter;" << endl;
00440 fout << " int fNOfSuppVec;" << endl;
00441 fout << " static float fAllSuppVectors[][" << fNsupv << "];" << endl;
00442 fout << " static float fAlphaTypeCoef[" << fNsupv << "];" << endl;
00443 fout << endl;
00444 fout << " // Kernel parameter(s) " << endl;
00445 fout << " float fGamma;" << endl;
00446 fout << "};" << endl;
00447 fout << "" << endl;
00448
00449
00450 fout << "inline void " << className << "::Initialize() " << endl;
00451 fout << "{" << endl;
00452 fout << " fBparameter = " << fBparm << ";" << endl;
00453 fout << " fNOfSuppVec = " << fNsupv << ";" << endl;
00454 fout << " fGamma = " << fGamma << ";" <<endl;
00455 fout << "}" << endl;
00456 fout << endl;
00457
00458
00459 fout << "inline double " << className << "::GetMvaValue__(const std::vector<double>& inputValues ) const" << endl;
00460 fout << "{" << endl;
00461 fout << " double mvaval = 0; " << endl;
00462 fout << " double temp = 0; " << endl;
00463 fout << endl;
00464 fout << " for (int ievt = 0; ievt < fNOfSuppVec; ievt++ ){" << endl;
00465 fout << " temp = 0;" << endl;
00466 fout << " for ( unsigned int ivar = 0; ivar < GetNvar(); ivar++ ) {" << endl;
00467
00468 fout << " temp += (fAllSuppVectors[ivar][ievt] - inputValues[ivar]) " << endl;
00469 fout << " * (fAllSuppVectors[ivar][ievt] - inputValues[ivar]); " << endl;
00470 fout << " }" << endl;
00471 fout << " mvaval += fAlphaTypeCoef[ievt] * exp( -fGamma * temp ); " << endl;
00472
00473 fout << " }" << endl;
00474 fout << " mvaval -= fBparameter;" << endl;
00475 fout << " return 1./(1. + exp(mvaval));" << endl;
00476 fout << "}" << endl;
00477 fout << "// Clean up" << endl;
00478 fout << "inline void " << className << "::Clear() " << endl;
00479 fout << "{" << endl;
00480 fout << " // nothing to clear " << endl;
00481 fout << "}" << endl;
00482 fout << "" << endl;
00483
00484
00485 fout << "float " << className << "::fAlphaTypeCoef[] =" << endl;
00486 fout << "{ ";
00487 for (Int_t isv = 0; isv < fNsupv; isv++) {
00488 fout << fSupportVectors->at(isv)->GetDeltaAlpha() * fSupportVectors->at(isv)->GetTypeFlag();
00489 if (isv < fNsupv-1) fout << ", ";
00490 }
00491 fout << " };" << endl << endl;
00492
00493 fout << "float " << className << "::fAllSuppVectors[][" << fNsupv << "] =" << endl;
00494 fout << "{";
00495 for (UInt_t ivar = 0; ivar < GetNvar(); ivar++) {
00496 fout << endl;
00497 fout << " { ";
00498 for (Int_t isv = 0; isv < fNsupv; isv++){
00499 fout << fSupportVectors->at(isv)->GetDataVector()->at(ivar);
00500 if (isv < fNsupv-1) fout << ", ";
00501 }
00502 fout << " }";
00503 if (ivar < GetNvar()-1) fout << ", " << endl;
00504 else fout << endl;
00505 }
00506 fout << "};" << endl<< endl;
00507 }
00508
00509
00510 void TMVA::MethodSVM::GetHelpMessage() const
00511 {
00512
00513
00514
00515
00516 Log() << Endl;
00517 Log() << gTools().Color("bold") << "--- Short description:" << gTools().Color("reset") << Endl;
00518 Log() << Endl;
00519 Log() << "The Support Vector Machine (SVM) builds a hyperplance separating" << Endl;
00520 Log() << "signal and background events (vectors) using the minimal subset of " << Endl;
00521 Log() << "all vectors used for training (support vectors). The extension to" << Endl;
00522 Log() << "the non-linear case is performed by mapping input vectors into a " << Endl;
00523 Log() << "higher-dimensional feature space in which linear separation is " << Endl;
00524 Log() << "possible. The use of the kernel functions thereby eliminates the " << Endl;
00525 Log() << "explicit transformation to the feature space. The implemented SVM " << Endl;
00526 Log() << "algorithm performs the classification tasks using linear, polynomial, " << Endl;
00527 Log() << "Gaussian and sigmoidal kernel functions. The Gaussian kernel allows " << Endl;
00528 Log() << "to apply any discriminant shape in the input space." << Endl;
00529 Log() << Endl;
00530 Log() << gTools().Color("bold") << "--- Performance optimisation:" << gTools().Color("reset") << Endl;
00531 Log() << Endl;
00532 Log() << "SVM is a general purpose non-linear classification method, which " << Endl;
00533 Log() << "does not require data preprocessing like decorrelation or Principal " << Endl;
00534 Log() << "Component Analysis. It generalises quite well and can handle analyses " << Endl;
00535 Log() << "with large numbers of input variables." << Endl;
00536 Log() << Endl;
00537 Log() << gTools().Color("bold") << "--- Performance tuning via configuration options:" << gTools().Color("reset") << Endl;
00538 Log() << Endl;
00539 Log() << "Optimal performance requires primarily a proper choice of the kernel " << Endl;
00540 Log() << "parameters (the width \"Sigma\" in case of Gaussian kernel) and the" << Endl;
00541 Log() << "cost parameter \"C\". The user must optimise them empirically by running" << Endl;
00542 Log() << "SVM several times with different parameter sets. The time needed for " << Endl;
00543 Log() << "each evaluation scales like the square of the number of training " << Endl;
00544 Log() << "events so that a coarse preliminary tuning should be performed on " << Endl;
00545 Log() << "reduced data sets." << Endl;
00546 }
