Changeset 7396 for trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc

Timestamp:

11/14/05 09:45:33 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc (modified) (15 diffs)

trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc

@@ -17 +17 @@
 !
 !   Author(s): Thomas Hengstebeck 2/2005 <mailto:hengsteb@physik.hu-berlin.de>
+!   Author(s): Thomas Bretz 8/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2005
@@ -31 +32 @@
 #include "MRFEnergyEst.h"
 
-#include <TFile.h>
-#include <TList.h>
-
-#include <TH1.h>
-#include <TH2.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include <TMath.h>
 #include <TVector.h>
 
@@ -45 +38 @@
 #include "MLog.h"
 #include "MLogManip.h"
+
+#include "MData.h"
+#include "MDataArray.h"
+
+#include "MRanForest.h"
+#include "MParameters.h"
 
 #include "MParList.h"
 #include "MTaskList.h"
 #include "MEvtLoop.h"
-
-#include "MRanTree.h"
-#include "MRanForest.h"
 #include "MRanForestGrow.h"
 
-#include "MData.h"
-#include "MParameters.h"
-
 ClassImp(MRFEnergyEst);
 
 using namespace std;
 
-static const TString gsDefName  = "MRFEnergyEst";
-static const TString gsDefTitle = "RF for energy estimation";
-
-// --------------------------------------------------------------------------
-//
-//  Default constructor. Set the name and title of this object
-//
+const TString MRFEnergyEst::gsDefName  = "MRFEnergyEst";
+const TString MRFEnergyEst::gsDefTitle = "RF for energy estimation";
+
 MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
     : fDebug(kFALSE), fData(0), fEnergyEst(0),
@@ -75 +64 @@
     fName  = name  ? name  : gsDefName.Data();
     fTitle = title ? title : gsDefTitle.Data();
+
+    // FIXME:
+    fNumTrees = 100; //100
+    fNumTry   = 5;   //3
+    fNdSize   = 0;   //1   0 means: in MRanForest estimated best value will be calculated
 }
 
@@ -82 +76 @@
 }
 
-// --------------------------------------------------------------------------
-//
-// Train the RF with the goven MHMatrix. The last column must contain the
-// True energy.
-//
-Int_t MRFEnergyEst::Train(const MHMatrix &matrixtrain, const TArrayD &grid)
+Int_t MRFEnergyEst::Train(const MHMatrix &matrixtrain, const TArrayD &grid, Int_t ver)
 {
     gLog.Separator("MRFEnergyEst - Train");
@@ -105 +94 @@
     }
 
-    const Int_t nbins = grid.GetSize()-1;
-    if (nbins<=0)
-    {
-        *fLog << err << "ERROR - Energy grid not vaild... abort." << endl;
-        return kFALSE;
-    }
-
     // rules (= combination of image par) to be used for energy estimation
     TFile fileRF(fFileName, "recreate");
@@ -122 +104 @@
     const Int_t nobs = 3; // Number of obsolete columns
 
+    MDataArray &dcol = *matrixtrain.GetColumns();
+
     MDataArray usedrules;
-    for(Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
-        usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());
-
-    // training of RF for each energy bin
+    for (Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
+        usedrules.AddEntry(dcol[i].GetRule());
+
+    MDataArray rules(usedrules);
+    rules.AddEntry(ver<2?"Classification":dcol[ncols-1].GetRule());
+
+    // prepare matrix for current energy bin
+    TMatrix mat(matrixtrain.GetM());
+
+    // last column must be removed (true energy col.)
+    mat.ResizeTo(nrows, ncols-nobs+1);
+
+    if (fDebug)
+        gLog.SetNullOutput(kTRUE);
+
+    const Int_t nbins = ver>0 ? 1 : grid.GetSize()-1;
     for (Int_t ie=0; ie<nbins; ie++)
     {
-        TMatrix mat1(nrows, ncols);
-        TMatrix mat0(nrows, ncols);
-
-        // prepare matrix for current energy bin
-        Int_t irow1=0;
-        Int_t irow0=0;
-
-        const TMatrix &m = matrixtrain.GetM();
-        for (Int_t j=0; j<nrows; j++)
+        switch (ver)
         {
-            const Double_t energy = m(j,ncols-1);
-
-            if (energy>grid[ie] && energy<=grid[ie+1])
-                TMatrixFRow(mat1, irow1++) = TMatrixFRow_const(m,j);
-            else
-                TMatrixFRow(mat0, irow0++) = TMatrixFRow_const(m,j);
+        case 0: // Replace Energy Grid by classification
+            {
+                Int_t irows=0;
+                for (Int_t j=0; j<nrows; j++)
+                {
+                    const Double_t energy = matrixtrain.GetM()(j,ncols-1);
+                    const Bool_t   inside = energy>grid[ie] && energy<=grid[ie+1];
+
+                    mat(j, ncols-nobs) = inside ? 1 : 0;
+
+                    if (inside)
+                        irows++;
+                }
+                if (irows==0)
+                    *fLog << warn << "WARNING - Skipping";
+                else
+                    *fLog << inf << "Training RF for";
+
+                *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irows << "/" << nrows << endl;
+
+                if (irows==0)
+                    continue;
+            }
+            break;
+
+        case 1: // Use Energy as classifier
+        case 2:
+            for (Int_t j=0; j<nrows; j++)
+                mat(j, ncols-nobs) = matrixtrain.GetM()(j,ncols-1);
+            break;
         }
 
-        const Bool_t invalid = irow1==0 || irow0==0;
-
-        if (invalid)
-            *fLog << warn << "WARNING - Skipping";
-        else
-            *fLog << inf << "Training RF for";
-
-        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irow0 << " " << irow1 << endl;
-
-        if (invalid)
-            continue;
-
-        if (fDebug)
-            gLog.SetNullOutput(kTRUE);
-
-        // last column must be removed (true energy col.)
-        mat1.ResizeTo(irow1, ncols-nobs);
-        mat0.ResizeTo(irow0, ncols-nobs);
-
-        // create MHMatrix as input for RF
-        MHMatrix matrix1(mat1, "MatrixHadrons");
-        MHMatrix matrix0(mat0, "MatrixGammas");
-
-        // training of RF
+        MHMatrix matrix(mat, &rules, "MatrixTrain");
+
+        MParList plist;
         MTaskList tlist;
-        MParList plist;
         plist.AddToList(&tlist);
-        plist.AddToList(&matrix0);
-        plist.AddToList(&matrix1);
+        plist.AddToList(&matrix);
+
+        MRanForest rf;
+        rf.SetNumTrees(fNumTrees);
+        rf.SetNumTry(fNumTry);
+        rf.SetNdSize(fNdSize);
+        rf.SetClassify(ver<2 ? 1 : 0);
+        if (ver==1)
+            rf.SetGrid(grid);
+
+        plist.AddToList(&rf);
 
         MRanForestGrow rfgrow;
-        rfgrow.SetNumTrees(fNumTrees); // number of trees
-        rfgrow.SetNumTry(fNumTry);     // number of trials in random split selection
-        rfgrow.SetNdSize(fNdSize);     // limit for nodesize
-
         tlist.AddToList(&rfgrow);
-    
+
         MEvtLoop evtloop;
-        evtloop.SetDisplay(fDisplay);
         evtloop.SetParList(&plist);
 
@@ -194 +188 @@
             gLog.SetNullOutput(kFALSE);
 
-        // Calculate bin center
-        const Double_t E = (TMath::Log10(grid[ie])+TMath::Log10(grid[ie+1]))/2;
-
-        // save whole forest
-        MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
-        forest->SetUserVal(E);
-        forest->Write(Form("%.10f", E));
+        if (ver==0)
+        {
+            // Calculate bin center
+            const Double_t E = (TMath::Log10(grid[ie])+TMath::Log10(grid[ie+1]))/2;
+
+            // save whole forest
+            rf.SetUserVal(E);
+            rf.SetName(Form("%.10f", E));
+        }
+
+        rf.Write();
     }
 
@@ -224 +222 @@
     while ((o=Next()))
     {
-        MRanForest *forest;
+        MRanForest *forest=0;
         fileRF.GetObject(o->GetName(), forest);
         if (!forest)
@@ -234 +232 @@
     }
 
+    // Maybe fEForests[0].fRules yould be used instead?
+
     if (fData->Read("rules")<=0)
     {
@@ -249 +249 @@
         return kFALSE;
 
+    cout << "MDataArray" << endl;
+
     fData = (MDataArray*)plist->FindCreateObj("MDataArray");
     if (!fData)
         return kFALSE;
 
+    cout << "ReadForests" << endl;
+
     if (!ReadForests(*plist))
     {
@@ -275 +279 @@
 }
 
-// --------------------------------------------------------------------------
-//
-//
 #include <TGraph.h>
 #include <TF1.h>
 Int_t MRFEnergyEst::Process()
 {
-    static TF1 f1("f1", "gaus");
-
     TVector event;
     if (fTestMatrix)
@@ -290 +289 @@
         *fData >> event;
 
+    // --------------- Single Tree RF -------------------
+    if (fEForests.GetEntries()==1)
+    {
+        const MRanForest *rf = (MRanForest*)fEForests[0];
+        fEnergyEst->SetVal(rf->CalcHadroness(event));
+        fEnergyEst->SetReadyToSave();
+
+        return kTRUE;
+    }
+
+    // --------------- Multi Tree RF -------------------
+    static TF1 f1("f1", "gaus");
+
     Double_t sume = 0;
     Double_t sumh = 0;
@@ -308 +320 @@
         const Double_t h = rf->CalcHadroness(event);
         const Double_t e = rf->GetUserVal();
+
         g.SetPoint(g.GetN(), e, h);
+
         sume += e*h;
         sumh += h;
+
         if (h>maxh)
         {
@@ -337 +352 @@
         fEnergyEst->SetVal(pow(10, f1.GetParameter(1)));
         break;
-
-    }
+    }
+
     fEnergyEst->SetReadyToSave();