77 TMVA::MethodHMatrix::MethodHMatrix( const
TString& jobName,
82 : TMVA::MethodBase( jobName, Types::kHMatrix, methodTitle, theData, theOption, theTargetDir )
96 : TMVA::
MethodBase(
Types::kHMatrix, theData, theWeightFile, theTargetDir )
111 fInvHMatrixS =
new TMatrixD( GetNvar(), GetNvar() );
112 fInvHMatrixB =
new TMatrixD( GetNvar(), GetNvar() );
113 fVecMeanS =
new TVectorD( GetNvar() );
114 fVecMeanB =
new TVectorD( GetNvar() );
117 SetSignalReferenceCut( 0.0 );
125 if (
NULL != fInvHMatrixS)
delete fInvHMatrixS;
126 if (
NULL != fInvHMatrixB)
delete fInvHMatrixB;
127 if (
NULL != fVecMeanS )
delete fVecMeanS;
128 if (
NULL != fVecMeanB )
delete fVecMeanB;
161 ComputeCovariance(
kTRUE, fInvHMatrixS );
162 ComputeCovariance(
kFALSE, fInvHMatrixB );
165 if (
TMath::Abs(fInvHMatrixS->Determinant()) < 10
E-24) {
166 Log() <<
kWARNING <<
"<Train> H-matrix S is almost singular with deterinant= "
168 <<
" did you use the variables that are linear combinations or highly correlated ???"
171 if (
TMath::Abs(fInvHMatrixB->Determinant()) < 10
E-24) {
172 Log() <<
kWARNING <<
"<Train> H-matrix B is almost singular with deterinant= "
174 <<
" did you use the variables that are linear combinations or highly correlated ???"
178 if (
TMath::Abs(fInvHMatrixS->Determinant()) < 10
E-120) {
179 Log() <<
kFATAL <<
"<Train> H-matrix S is singular with deterinant= "
181 <<
" did you use the variables that are linear combinations ???"
184 if (
TMath::Abs(fInvHMatrixB->Determinant()) < 10
E-120) {
185 Log() <<
kFATAL <<
"<Train> H-matrix B is singular with deterinant= "
187 <<
" did you use the variables that are linear combinations ???"
192 fInvHMatrixS->Invert();
193 fInvHMatrixB->Invert();
203 const UInt_t nvar = DataInfo().GetNVariables();
208 TMatrixD mat2(nvar, nvar); mat2 *= 0;
215 for (
Int_t i=0, iEnd=
Data()->GetNEvents(); i<iEnd; ++i) {
218 const Event* origEvt =
Data()->GetEvent(i);
222 if (IgnoreEventsWithNegWeightsInTraining() && weight <= 0)
continue;
224 if (DataInfo().IsSignal(origEvt) != isSignal)
continue;
227 GetTransformationHandler().SetTransformationReferenceClass( origEvt->
GetClass() );
228 const Event* ev = GetTransformationHandler().Transform( origEvt );
231 sumOfWeights += weight;
234 for (ivar=0; ivar<nvar; ivar++) xval[ivar] = ev->
GetValue(ivar);
237 for (ivar=0; ivar<nvar; ivar++) {
239 vec(ivar) += xval[ivar]*weight;
240 mat2(ivar, ivar) += (xval[ivar]*xval[ivar])*weight;
242 for (jvar=ivar+1; jvar<nvar; jvar++) {
243 mat2(ivar, jvar) += (xval[ivar]*xval[jvar])*weight;
244 mat2(jvar, ivar) = mat2(ivar, jvar);
250 for (ivar=0; ivar<nvar; ivar++) {
252 if (isSignal) (*fVecMeanS)(ivar) = vec(ivar)/sumOfWeights;
253 else (*fVecMeanB)(ivar) = vec(ivar)/sumOfWeights;
255 for (jvar=0; jvar<nvar; jvar++) {
256 (*mat)(ivar, jvar) = mat2(ivar, jvar)/sumOfWeights - vec(ivar)*vec(jvar)/(sumOfWeights*sumOfWeights);
271 if (s+b < 0)
Log() <<
kFATAL <<
"big trouble: s+b: " << s+b <<
Endl;
274 NoErrorCalc(err, errUpper);
276 return (b - s)/(s + b);
286 const Event* origEvt = fTmpEvent ? fTmpEvent:
Data()->GetEvent();
289 UInt_t ivar(0), jvar(0), nvar(GetNvar());
290 std::vector<Double_t> val( nvar );
294 GetTransformationHandler().SetTransformationReferenceClass( fSignalClass );
296 GetTransformationHandler().SetTransformationReferenceClass( fBackgroundClass );
298 const Event* ev = GetTransformationHandler().Transform( origEvt );
300 for (ivar=0; ivar<nvar; ivar++) val[ivar] = ev->
GetValue( ivar );
303 for (ivar=0; ivar<nvar; ivar++) {
304 for (jvar=0; jvar<nvar; jvar++) {
306 chi2 += ( (val[ivar] - (*fVecMeanS)(ivar))*(val[jvar] - (*fVecMeanS)(jvar))
307 * (*fInvHMatrixS)(ivar,jvar) );
309 chi2 += ( (val[ivar] - (*fVecMeanB)(ivar))*(val[jvar] - (*fVecMeanB)(jvar))
310 * (*fInvHMatrixB)(ivar,jvar) );
315 if (chi2 < 0)
Log() <<
kFATAL <<
"<GetChi2> negative chi2: " << chi2 <<
Endl;
360 for (ivar=0; ivar<GetNvar(); ivar++)
361 istr >> (*fVecMeanS)(ivar) >> (*fVecMeanB)(ivar);
364 for (ivar=0; ivar<GetNvar(); ivar++)
365 for (jvar=0; jvar<GetNvar(); jvar++)
366 istr >> (*fInvHMatrixS)(ivar,jvar);
369 for (ivar=0; ivar<GetNvar(); ivar++)
370 for (jvar=0; jvar<GetNvar(); jvar++)
371 istr >> (*fInvHMatrixB)(ivar,jvar);
379 fout <<
" // arrays of input evt vs. variable " << std::endl;
380 fout <<
" double fInvHMatrixS[" << GetNvar() <<
"][" << GetNvar() <<
"]; // inverse H-matrix (signal)" << std::endl;
381 fout <<
" double fInvHMatrixB[" << GetNvar() <<
"][" << GetNvar() <<
"]; // inverse H-matrix (background)" << std::endl;
382 fout <<
" double fVecMeanS[" << GetNvar() <<
"]; // vector of mean values (signal)" << std::endl;
383 fout <<
" double fVecMeanB[" << GetNvar() <<
"]; // vector of mean values (background)" << std::endl;
384 fout <<
" " << std::endl;
385 fout <<
" double GetChi2( const std::vector<double>& inputValues, int type ) const;" << std::endl;
386 fout <<
"};" << std::endl;
387 fout <<
" " << std::endl;
388 fout <<
"void " << className <<
"::Initialize() " << std::endl;
389 fout <<
"{" << std::endl;
390 fout <<
" // init vectors with mean values" << std::endl;
391 for (
UInt_t ivar=0; ivar<GetNvar(); ivar++) {
392 fout <<
" fVecMeanS[" << ivar <<
"] = " << (*fVecMeanS)(ivar) <<
";" << std::endl;
393 fout <<
" fVecMeanB[" << ivar <<
"] = " << (*fVecMeanB)(ivar) <<
";" << std::endl;
395 fout <<
" " << std::endl;
396 fout <<
" // init H-matrices" << std::endl;
397 for (
UInt_t ivar=0; ivar<GetNvar(); ivar++) {
398 for (
UInt_t jvar=0; jvar<GetNvar(); jvar++) {
399 fout <<
" fInvHMatrixS[" << ivar <<
"][" << jvar <<
"] = "
400 << (*fInvHMatrixS)(ivar,jvar) <<
";" << std::endl;
401 fout <<
" fInvHMatrixB[" << ivar <<
"][" << jvar <<
"] = "
402 << (*fInvHMatrixB)(ivar,jvar) <<
";" << std::endl;
405 fout <<
"}" << std::endl;
406 fout <<
" " << std::endl;
407 fout <<
"inline double " << className <<
"::GetMvaValue__( const std::vector<double>& inputValues ) const" << std::endl;
408 fout <<
"{" << std::endl;
409 fout <<
" // returns the H-matrix signal estimator" << std::endl;
410 fout <<
" std::vector<double> inputValuesSig = inputValues;" << std::endl;
411 fout <<
" std::vector<double> inputValuesBgd = inputValues;" << std::endl;
412 if (GetTransformationHandler().GetTransformationList().GetSize() != 0) {
414 UInt_t signalClass =DataInfo().GetClassInfo(
"Signal")->GetNumber();
415 UInt_t backgroundClass=DataInfo().GetClassInfo(
"Background")->GetNumber();
417 fout <<
" Transform(inputValuesSig," << signalClass <<
");" << std::endl;
418 fout <<
" Transform(inputValuesBgd," << backgroundClass <<
");" << std::endl;
423 fout <<
" double s = GetChi2( inputValuesSig, " <<
Types::kSignal <<
" );" << std::endl;
424 fout <<
" double b = GetChi2( inputValuesBgd, " <<
Types::kBackground <<
" );" << std::endl;
428 fout <<
" " << std::endl;
429 fout <<
" if (s+b <= 0) std::cout << \"Problem in class " << className <<
"::GetMvaValue__: s+b = \"" << std::endl;
430 fout <<
" << s+b << \" <= 0 \" << std::endl;" << std::endl;
431 fout <<
" " << std::endl;
432 fout <<
" return (b - s)/(s + b);" << std::endl;
433 fout <<
"}" << std::endl;
434 fout <<
" " << std::endl;
435 fout <<
"inline double " << className <<
"::GetChi2( const std::vector<double>& inputValues, int type ) const" << std::endl;
436 fout <<
"{" << std::endl;
437 fout <<
" // compute chi2-estimator for event according to type (signal/background)" << std::endl;
438 fout <<
" " << std::endl;
439 fout <<
" size_t ivar,jvar;" << std::endl;
440 fout <<
" double chi2 = 0;" << std::endl;
441 fout <<
" for (ivar=0; ivar<GetNvar(); ivar++) {" << std::endl;
442 fout <<
" for (jvar=0; jvar<GetNvar(); jvar++) {" << std::endl;
444 fout <<
" chi2 += ( (inputValues[ivar] - fVecMeanS[ivar])*(inputValues[jvar] - fVecMeanS[jvar])" << std::endl;
445 fout <<
" * fInvHMatrixS[ivar][jvar] );" << std::endl;
446 fout <<
" else" << std::endl;
447 fout <<
" chi2 += ( (inputValues[ivar] - fVecMeanB[ivar])*(inputValues[jvar] - fVecMeanB[jvar])" << std::endl;
448 fout <<
" * fInvHMatrixB[ivar][jvar] );" << std::endl;
449 fout <<
" }" << std::endl;
450 fout <<
" } // loop over variables " << std::endl;
451 fout <<
" " << std::endl;
452 fout <<
" // sanity check" << std::endl;
453 fout <<
" if (chi2 < 0) std::cout << \"Problem in class " << className <<
"::GetChi2: chi2 = \"" << std::endl;
454 fout <<
" << chi2 << \" < 0 \" << std::endl;" << std::endl;
455 fout <<
" " << std::endl;
456 fout <<
" return chi2;" << std::endl;
457 fout <<
"}" << std::endl;
458 fout <<
" " << std::endl;
459 fout <<
"// Clean up" << std::endl;
460 fout <<
"inline void " << className <<
"::Clear() " << std::endl;
461 fout <<
"{" << std::endl;
462 fout <<
" // nothing to clear" << std::endl;
463 fout <<
"}" << std::endl;
477 Log() <<
"The H-Matrix classifier discriminates one class (signal) of a feature" <<
Endl;
478 Log() <<
"vector from another (background). The correlated elements of the" <<
Endl;
479 Log() <<
"vector are assumed to be Gaussian distributed, and the inverse of" <<
Endl;
480 Log() <<
"the covariance matrix is the H-Matrix. A multivariate chi-squared" <<
Endl;
481 Log() <<
"estimator is built that exploits differences in the mean values of" <<
Endl;
482 Log() <<
"the vector elements between the two classes for the purpose of" <<
Endl;
483 Log() <<
"discrimination." <<
Endl;
487 Log() <<
"The TMVA implementation of the H-Matrix classifier has been shown" <<
Endl;
488 Log() <<
"to underperform in comparison with the corresponding Fisher discriminant," <<
Endl;
489 Log() <<
"when using similar assumptions and complexity. Its use is therefore" <<
Endl;
490 Log() <<
"depreciated. Only in cases where the background model is strongly" <<
Endl;
491 Log() <<
"non-Gaussian, H-Matrix may perform better than Fisher. In such" <<
Endl;
492 Log() <<
"occurrences the user is advised to employ non-linear classifiers. " <<
Endl;
void ReadWeightsFromXML(void *wghtnode)
read weights from XML file
void Init()
default initialization called by all constructors
MsgLogger & Endl(MsgLogger &ml)
void AddWeightsXMLTo(void *parent) const
create XML description for HMatrix classification
Double_t GetMvaValue(Double_t *err=0, Double_t *errUpper=0)
returns the H-matrix signal estimator
void ProcessOptions()
process user options
virtual Bool_t HasAnalysisType(Types::EAnalysisType type, UInt_t numberClasses, UInt_t numberTargets)
FDA can handle classification with 2 classes and regression with one regression-target.
Double_t GetWeight() const
return the event weight - depending on whether the flag IgnoreNegWeightsInTraining is or not...
void ReadWeightsFromStream(std::istream &istr)
read variable names and min/max NOTE: the latter values are mandatory for the normalisation in the re...
virtual ~MethodHMatrix()
destructor
void GetHelpMessage() const
get help message text
ClassImp(TMVA::MethodHMatrix) TMVA
standard constructor for the H-Matrix method
std::vector< std::vector< double > > Data
MethodHMatrix(const TString &jobName, const TString &methodTitle, DataSetInfo &theData, const TString &theOption="", TDirectory *theTargetDir=0)
TVectorT< Double_t > TVectorD
TMatrixT< Double_t > TMatrixD
void Train()
computes H-matrices for signal and background samples
Describe directory structure in memory.
void ComputeCovariance(Bool_t, TMatrixD *)
compute covariance matrix
void MakeClassSpecific(std::ostream &, const TString &) const
write Fisher-specific classifier response
static RooMathCoreReg dummy
Float_t GetValue(UInt_t ivar) const
return value of i'th variable
#define REGISTER_METHOD(CLASS)
for example
void DeclareOptions()
MethodHMatrix options: none (apart from those implemented in MethodBase)
Double_t GetChi2(Types::ESBType)
compute chi2-estimator for event according to type (signal/background)