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

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mranforest/MRanForest.cc (modified) (14 diffs)

trunk/MagicSoft/Mars/mranforest/MRanForest.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Hengstebeck 3/2003 <mailto:hengsteb@alwa02.physik.uni-siegen.de>
+!   Author(s): Thomas Hengstebeck 3/2003 <mailto:hengsteb@physik.hu-berlin.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2003
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -39 +39 @@
 // split selection (which is subject to MRanTree::GrowTree())
 //
-// Version 2:
-//  + fUserVal
-//
 /////////////////////////////////////////////////////////////////////////////
 #include "MRanForest.h"
 
-#include <TMatrix.h>
-#include <TRandom3.h>
+#include <TVector.h>
+#include <TRandom.h>
 
 #include "MHMatrix.h"
 #include "MRanTree.h"
+#include "MData.h"
+#include "MDataArray.h"
+#include "MParList.h"
 
 #include "MLog.h"
@@ -62 +62 @@
 // Default constructor.
 //
-MRanForest::MRanForest(const char *name, const char *title) : fNumTrees(100), fRanTree(NULL),fUsePriors(kFALSE), fUserVal(-1)
+MRanForest::MRanForest(const char *name, const char *title) : fClassify(1), fNumTrees(100), fNumTry(0), fNdSize(1), fRanTree(NULL), fUserVal(-1)
 {
     fName  = name  ? name  : "MRanForest";
@@ -71 +71 @@
 }
 
+MRanForest::MRanForest(const MRanForest &rf)
+{
+    // Copy constructor
+    fName  = rf.fName;
+    fTitle = rf.fTitle;
+
+    fClassify = rf.fClassify;
+    fNumTrees = rf.fNumTrees;
+    fNumTry   = rf.fNumTry;
+    fNdSize   = rf.fNdSize;
+    fTreeNo   = rf.fTreeNo;
+    fRanTree  = NULL;
+
+    fRules=new MDataArray();
+    fRules->Reset();
+
+    MDataArray *newrules=rf.fRules;
+
+    for(Int_t i=0;i<newrules->GetNumEntries();i++)
+    {
+        MData &data=(*newrules)[i];
+        fRules->AddEntry(data.GetRule());
+    }
+
+    // trees
+    fForest=new TObjArray();
+    fForest->SetOwner(kTRUE);
+
+    TObjArray *newforest=rf.fForest;
+    for(Int_t i=0;i<newforest->GetEntries();i++)
+    {
+        MRanTree *rantree=(MRanTree*)newforest->At(i);
+
+        MRanTree *newtree=new MRanTree(*rantree);
+        fForest->Add(newtree);
+    }
+
+    fHadTrue  = rf.fHadTrue;
+    fHadEst   = rf.fHadEst;
+    fDataSort = rf.fDataSort;
+    fDataRang = rf.fDataRang;
+    fClassPop = rf.fClassPop;
+    fWeight   = rf.fWeight;
+    fTreeHad  = rf.fTreeHad;
+
+    fNTimesOutBag = rf.fNTimesOutBag;
+
+}
+
 // --------------------------------------------------------------------------
-//
 // Destructor. 
-//
 MRanForest::~MRanForest()
 {
     delete fForest;
+}
+
+MRanTree *MRanForest::GetTree(Int_t i)
+{
+    return (MRanTree*)(fForest->At(i));
 }
 
@@ -88 +140 @@
 }
 
-void MRanForest::SetPriors(Float_t prior_had, Float_t prior_gam)
-{
-    const Float_t sum=prior_gam+prior_had;
-
-    prior_gam/=sum;
-    prior_had/=sum;
-
-    fPrior[0]=prior_had;
-    fPrior[1]=prior_gam;
-
-    fUsePriors=kTRUE;
+void MRanForest::SetNumTry(Int_t n)
+{
+    fNumTry=TMath::Max(n,0);
+}
+
+void MRanForest::SetNdSize(Int_t n)
+{
+    fNdSize=TMath::Max(n,1);
+}
+
+void MRanForest::SetWeights(const TArrayF &weights)
+{
+    const int n=weights.GetSize();
+    fWeight.Set(n);
+    fWeight=weights;
+
+    return;
+}
+
+void MRanForest::SetGrid(const TArrayD &grid)
+{
+    const int n=grid.GetSize();
+
+    for(int i=0;i<n-1;i++)
+        if(grid[i]>=grid[i+1])
+        {
+            *fLog<<inf<<"Grid points must be in increasing order! Ignoring grid."<<endl;
+            return;
+        }
+
+    fGrid=grid;
+
+    //*fLog<<inf<<"Following "<<n<<" grid points are used:"<<endl;
+    //for(int i=0;i<n;i++)
+    //    *fLog<<inf<<" "<<i<<") "<<fGrid[i]<<endl;
+
+    return;
 }
 
 Int_t MRanForest::GetNumDim() const
 {
-    return fGammas ? fGammas->GetM().GetNcols() : 0;
-}
-
+    return fMatrix ? fMatrix->GetNcols() : 0;
+}
+
+Int_t MRanForest::GetNumData() const
+{
+    return fMatrix ? fMatrix->GetNrows() : 0;
+}
+
+Int_t MRanForest::GetNclass() const
+{
+    int maxidx = TMath::LocMax(fClass.GetSize(),fClass.GetArray());
+
+    return int(fClass[maxidx])+1;
+}
+
+void MRanForest::PrepareClasses()
+{
+    const int numdata=fHadTrue.GetSize();
+
+    if(fGrid.GetSize()>0)
+    {
+        // classes given by grid
+        const int ngrid=fGrid.GetSize();
+
+        for(int j=0;j<numdata;j++)
+        {
+            // Array is supposed  to be sorted prior to this call.
+            // If match is found, function returns position of element.
+            // If no match found, function gives nearest element smaller
+            // than value.
+            int k=TMath::BinarySearch(ngrid, fGrid.GetArray(), fHadTrue[j]);
+
+            fClass[j]   = k;
+        }
+
+        int minidx = TMath::LocMin(fClass.GetSize(),fClass.GetArray());
+        int min = fClass[minidx];
+        if(min!=0) for(int j=0;j<numdata;j++)fClass[j]-=min;
+
+    }else{
+        // classes directly given
+        for (Int_t j=0;j<numdata;j++)
+            fClass[j] = int(fHadTrue[j]+0.5);
+    }
+
+    return;
+}
+
+/*
+Bool_t MRanForest::PreProcess(MParList *plist)
+{
+    if (!fRules)
+    {
+        *fLog << err << dbginf << "MDataArray with rules not initialized... aborting." << endl;
+        return kFALSE;
+    }
+
+    if (!fRules->PreProcess(plist))
+    {
+        *fLog << err << dbginf << "PreProcessing of MDataArray failed... aborting." << endl;
+        return kFALSE;
+    }
+
+    return kTRUE;
+}
+*/
+
+Double_t MRanForest::CalcHadroness()
+{
+    TVector event;
+    *fRules >> event;
+
+    return CalcHadroness(event);
+}
 
 Double_t MRanForest::CalcHadroness(const TVector &event)
@@ -117 +266 @@
     while ((tree=(MRanTree*)Next()))
     {
-        fTreeHad[ntree]=tree->TreeHad(event);
-        hadroness+=fTreeHad[ntree];
+        hadroness+=(fTreeHad[ntree]=tree->TreeHad(event));
         ntree++;
     }
@@ -126 +274 @@
 Bool_t MRanForest::AddTree(MRanTree *rantree=NULL)
 {
-    if (rantree)
-        fRanTree=rantree;
-    if (!fRanTree)
+    fRanTree = rantree ? rantree:fRanTree;
+
+    if (!fRanTree) return kFALSE;
+
+    MRanTree *newtree=new MRanTree(*fRanTree);
+    fForest->Add(newtree);
+
+    return kTRUE;
+}
+
+Bool_t MRanForest::SetupGrow(MHMatrix *mat,MParList *plist)
+{
+    //-------------------------------------------------------------------
+    // access matrix, copy last column (target) preliminarily
+    // into fHadTrue
+    TMatrix mat_tmp = mat->GetM();
+    int dim         = mat_tmp.GetNcols();
+    int numdata     = mat_tmp.GetNrows();
+
+    fMatrix=new TMatrix(mat_tmp);
+
+    fHadTrue.Set(numdata);
+    fHadTrue.Reset(0);
+
+    for (Int_t j=0;j<numdata;j++)
+        fHadTrue[j] = (*fMatrix)(j,dim-1);
+
+    // remove last col
+    fMatrix->ResizeTo(numdata,dim-1);
+    dim=fMatrix->GetNcols();
+
+    //-------------------------------------------------------------------
+    // setup labels for classification/regression
+    fClass.Set(numdata);
+    fClass.Reset(0);
+
+    if(fClassify) PrepareClasses();
+
+    //-------------------------------------------------------------------
+    // allocating and initializing arrays
+    fHadEst.Set(numdata);       fHadEst.Reset(0);
+    fNTimesOutBag.Set(numdata); fNTimesOutBag.Reset(0);
+    fDataSort.Set(dim*numdata); fDataSort.Reset(0);
+    fDataRang.Set(dim*numdata); fDataRang.Reset(0);
+
+    if(fWeight.GetSize()!=numdata)
+    {
+        fWeight.Set(numdata);
+        fWeight.Reset(1.);
+        *fLog << inf <<"Setting weights to 1 (no weighting)"<< endl;
+    }
+
+    //-------------------------------------------------------------------
+    // setup rules to be used for classification/regression
+    MDataArray *allrules=(MDataArray*)mat->GetColumns();
+    if(allrules==NULL)
+    {
+        *fLog << err <<"Rules of matrix == null, exiting"<< endl;
         return kFALSE;
-
-    fForest->Add((MRanTree*)fRanTree->Clone());
-
-    return kTRUE;
-}
-
-Int_t MRanForest::GetNumData() const
-{
-    return fHadrons && fGammas ? fHadrons->GetM().GetNrows()+fGammas->GetM().GetNrows() : 0;
-}
-
-Bool_t MRanForest::SetupGrow(MHMatrix *mhad,MHMatrix *mgam)
-{
-    // pointer to training data
-    fHadrons=mhad;
-    fGammas=mgam;
-
-    // determine data entries and dimension of Hillas-parameter space
-    //fNumHad=fHadrons->GetM().GetNrows();
-    //fNumGam=fGammas->GetM().GetNrows();
-
-    const Int_t dim = GetNumDim();
-
-    if (dim!=fGammas->GetM().GetNcols())
-        return kFALSE;
-
-    const Int_t numdata = GetNumData();
-
-    // allocating and initializing arrays
-    fHadTrue.Set(numdata);
-    fHadTrue.Reset();
-    fHadEst.Set(numdata);
-
-    fPrior.Set(2);
-    fClassPop.Set(2);
-    fWeight.Set(numdata);
-    fNTimesOutBag.Set(numdata);
-    fNTimesOutBag.Reset();
-
-    fDataSort.Set(dim*numdata);
-    fDataRang.Set(dim*numdata);
-
-    // calculating class populations (= no. of gammas and hadrons)
-    fClassPop.Reset();
-    for(Int_t n=0;n<numdata;n++)
-        fClassPop[fHadTrue[n]]++;
-
-    // setting weights taking into account priors
-    if (!fUsePriors)
-        fWeight.Reset(1.);
-    else
-    {
-        for(Int_t j=0;j<2;j++)
-            fPrior[j] *= (fClassPop[j]>=1) ? (Float_t)numdata/fClassPop[j]:0;
-
-        for(Int_t n=0;n<numdata;n++)
-            fWeight[n]=fPrior[fHadTrue[n]];
-    }
-
-    // create fDataSort
-    CreateDataSort();
+    }
+
+    fRules=new MDataArray(); fRules->Reset();
+    TString target_rule;
+
+    for(Int_t i=0;i<dim+1;i++)
+    {
+        MData &data=(*allrules)[i];
+        if(i<dim)
+            fRules->AddEntry(data.GetRule());
+        else
+            target_rule=data.GetRule();
+    }
+
+    *fLog << inf <<endl;
+    *fLog << inf <<"Setting up RF for training on target:"<<endl<<" "<<target_rule.Data()<<endl;
+    *fLog << inf <<"Following rules are used as input to RF:"<<endl;
+
+    for(Int_t i=0;i<dim;i++)
+    {
+        MData &data=(*fRules)[i];
+        *fLog<<inf<<" "<<i<<") "<<data.GetRule()<<endl<<flush;
+    }
+
+    *fLog << inf <<endl;
+
+    //-------------------------------------------------------------------
+    // prepare (sort) data for fast optimization algorithm
+    if(!CreateDataSort()) return kFALSE;
+
+    //-------------------------------------------------------------------
+    // access and init tree container
+    fRanTree = (MRanTree*)plist->FindCreateObj("MRanTree");
 
     if(!fRanTree)
@@ -197 +372 @@
         return kFALSE;
     }
-    fRanTree->SetRules(fGammas->GetColumns());
+
+    fRanTree->SetClassify(fClassify);
+    fRanTree->SetNdSize(fNdSize);
+
+    if(fNumTry==0)
+    {
+        double ddim = double(dim);
+
+        fNumTry=int(sqrt(ddim)+0.5);
+        *fLog<<inf<<endl;
+        *fLog<<inf<<"Set no. of trials to the recommended value of round("<<sqrt(ddim)<<") = ";
+        *fLog<<inf<<fNumTry<<endl;
+
+    }
+    fRanTree->SetNumTry(fNumTry);
+
+    *fLog<<inf<<endl;
+    *fLog<<inf<<"Following settings for the tree growing are used:"<<endl;
+    *fLog<<inf<<" Number of Trees : "<<fNumTrees<<endl;
+    *fLog<<inf<<" Number of Trials: "<<fNumTry<<endl;
+    *fLog<<inf<<" Final Node size : "<<fNdSize<<endl;
+
     fTreeNo=0;
 
     return kTRUE;
-}
-
-void MRanForest::InitHadEst(Int_t from, Int_t to, const TMatrix &m, TArrayI &jinbag)
-{
-    for (Int_t ievt=from;ievt<to;ievt++)
-    {
-        if (jinbag[ievt]>0)
-            continue;
-        fHadEst[ievt] += fRanTree->TreeHad(m, ievt-from);
-        fNTimesOutBag[ievt]++;
-    }
 }
 
@@ -224 +409 @@
     fTreeNo++;
 
+    //-------------------------------------------------------------------
     // initialize running output
+
+    float minfloat=fHadTrue[TMath::LocMin(fHadTrue.GetSize(),fHadTrue.GetArray())];
+    Bool_t calcResolution=(minfloat>0.001);
+
     if (fTreeNo==1)
     {
         *fLog << inf << endl;
-        *fLog << underline; // << "1st col        2nd col" << endl;
-        *fLog << "no. of tree    error in % (calulated using oob-data -> overestim. of error)" << endl;
+        *fLog << underline;
+
+        if(calcResolution)
+            *fLog << "no. of tree    no. of nodes    resolution in % (from oob-data -> overest. of error)" << endl;
+        else
+            *fLog << "no. of tree    no. of nodes    rms in % (from oob-data -> overest. of error)" << endl;
+                     //        12345678901234567890123456789012345678901234567890
     }
 
     const Int_t numdata = GetNumData();
-
+    const Int_t nclass  = GetNclass();
+
+    //-------------------------------------------------------------------
     // bootstrap aggregating (bagging) -> sampling with replacement:
-    //
-    // The integer k is randomly (uniformly) chosen from the set
-    // {0,1,...,fNumData-1}, which is the set of the index numbers of
-    // all events in the training sample
-    TArrayF classpopw(2);
+
+    TArrayF classpopw(nclass);
     TArrayI jinbag(numdata); // Initialization includes filling with 0
     TArrayF winbag(numdata); // Initialization includes filling with 0
 
+    float square=0; float mean=0;
+
     for (Int_t n=0; n<numdata; n++)
     {
+        // The integer k is randomly (uniformly) chosen from the set
+        // {0,1,...,numdata-1}, which is the set of the index numbers of
+        // all events in the training sample
+  
         const Int_t k = Int_t(gRandom->Rndm()*numdata);
 
-        classpopw[fHadTrue[k]]+=fWeight[k];
+        if(fClassify)
+            classpopw[fClass[k]]+=fWeight[k];
+        else
+            classpopw[0]+=fWeight[k];
+
+        mean  +=fHadTrue[k]*fWeight[k];
+        square+=fHadTrue[k]*fHadTrue[k]*fWeight[k];
+  
         winbag[k]+=fWeight[k];
         jinbag[k]=1;
-    }
-
+
+    }
+
+    //-------------------------------------------------------------------
     // modifying sorted-data array for in-bag data:
-    //
+
     // In bagging procedure ca. 2/3 of all elements in the original
     // training sample are used to build the in-bag data
@@ -261 +470 @@
     ModifyDataSort(datsortinbag, ninbag, jinbag);
 
-    const TMatrix &hadrons=fHadrons->GetM();
-    const TMatrix &gammas =fGammas->GetM();
-
-    // growing single tree
-    fRanTree->GrowTree(hadrons,gammas,fHadTrue,datsortinbag,
-                       fDataRang,classpopw,jinbag,winbag);
-
+    fRanTree->GrowTree(fMatrix,fHadTrue,fClass,datsortinbag,fDataRang,classpopw,mean,square,
+                       jinbag,winbag,nclass);
+
+    //-------------------------------------------------------------------
     // error-estimates from out-of-bag data (oob data):
     //
@@ -277 +483 @@
     // determined from oob-data is underestimated, but can still be taken as upper limit.
 
-    const Int_t numhad = hadrons.GetNrows();
-    InitHadEst(0, numhad, hadrons, jinbag);
-    InitHadEst(numhad, numdata, gammas, jinbag);
-    /*
-    for (Int_t ievt=0;ievt<numhad;ievt++)
-    {
-        if (jinbag[ievt]>0)
-            continue;
-        fHadEst[ievt] += fRanTree->TreeHad(hadrons, ievt);
+    for (Int_t ievt=0;ievt<numdata;ievt++)
+    {
+        if (jinbag[ievt]>0) continue;
+
+        fHadEst[ievt] +=fRanTree->TreeHad((*fMatrix), ievt);
         fNTimesOutBag[ievt]++;
-    }
-
-    for (Int_t ievt=numhad;ievt<numdata;ievt++)
-    {
-        if (jinbag[ievt]>0)
-            continue;
-        fHadEst[ievt] += fRanTree->TreeHad(gammas, ievt-numhad);
-        fNTimesOutBag[ievt]++;
-    }
-    */
+
+    }
+
     Int_t n=0;
-    Double_t ferr=0;
+    double ferr=0;
+
     for (Int_t ievt=0;ievt<numdata;ievt++)
-        if (fNTimesOutBag[ievt]!=0)
+    {
+        if(fNTimesOutBag[ievt]!=0)
         {
-            const Double_t val = fHadEst[ievt]/fNTimesOutBag[ievt]-fHadTrue[ievt];
+            float val = fHadEst[ievt]/float(fNTimesOutBag[ievt])-fHadTrue[ievt];
+            if(calcResolution) val/=fHadTrue[ievt];
+
             ferr += val*val;
             n++;
         }
-
+    }
     ferr = TMath::Sqrt(ferr/n);
 
+    //-------------------------------------------------------------------
     // give running output
-    *fLog << inf << setw(5) << fTreeNo << Form("%15.2f", ferr*100) << endl;
+    *fLog << inf << setw(5)  << fTreeNo;
+    *fLog << inf << setw(20) << fRanTree->GetNumEndNodes();
+    *fLog << inf << Form("%20.2f", ferr*100.);
+    *fLog << inf << endl;
 
     // adding tree to forest
@@ -318 +521 @@
 }
 
-void MRanForest::CreateDataSort()
-{
-    // The values of concatenated data arrays fHadrons and fGammas (denoted in the following by fData,
-    // which does actually not exist) are sorted from lowest to highest.
+Bool_t MRanForest::CreateDataSort()
+{
+    // fDataSort(m,n) is the event number in which fMatrix(m,n) occurs.
+    // fDataRang(m,n) is the rang of fMatrix(m,n), i.e. if rang = r:
+    //   fMatrix(m,n) is the r-th highest value of all fMatrix(m,.).
     //
-    //
-    //                   fHadrons(0,0) ... fHadrons(0,nhad-1)   fGammas(0,0) ... fGammas(0,ngam-1)
-    //                        .                 .                   .                .
-    //  fData(m,n)   =        .                 .                   .                .
-    //                        .                 .                   .                .
-    //                   fHadrons(m-1,0)...fHadrons(m-1,nhad-1) fGammas(m-1,0)...fGammas(m-1,ngam-1)
-    //
-    //
-    // Then fDataSort(m,n) is the event number in which fData(m,n) occurs.
-    // fDataRang(m,n) is the rang of fData(m,n), i.e. if rang = r, fData(m,n) is the r-th highest
-    // value of all fData(m,.). There may be more then 1 event with rang r (due to bagging).
+    // There may be more then 1 event with rang r (due to bagging).
+
     const Int_t numdata = GetNumData();
+    const Int_t dim = GetNumDim();
 
     TArrayF v(numdata);
     TArrayI isort(numdata);
 
-    const TMatrix &hadrons=fHadrons->GetM();
-    const TMatrix &gammas=fGammas->GetM();
-
-    const Int_t numgam = gammas.GetNrows();
-    const Int_t numhad = hadrons.GetNrows();
-
-    for (Int_t j=0;j<numhad;j++)
-        fHadTrue[j]=1;
-
-    for (Int_t j=0;j<numgam;j++)
-        fHadTrue[j+numhad]=0;
-
-    const Int_t dim = GetNumDim();
+
     for (Int_t mvar=0;mvar<dim;mvar++)
     {
-        for(Int_t n=0;n<numhad;n++)
+
+        for(Int_t n=0;n<numdata;n++)
         {
-            v[n]=hadrons(n,mvar);
+            v[n]=(*fMatrix)(n,mvar);
             isort[n]=n;
-        }
-
-        for(Int_t n=0;n<numgam;n++)
-        {
-            v[n+numhad]=gammas(n,mvar);
-            isort[n+numhad]=n;
+
+            if(TMath::IsNaN(v[n]))
+            {
+                *fLog << err <<"Event no. "<<n<<", matrix column no. "<<mvar;
+                *fLog << err <<" has the value NaN."<<endl;
+                return kFALSE;
+            }
         }
 
@@ -371 +557 @@
         // of that v[n], which is the n-th from the lowest (assume the original
         // event numbers are 0,1,...).
+
+        // control sorting
+        for(int n=1;n<numdata;n++)
+            if(v[isort[n-1]]>v[isort[n]])
+            {
+                *fLog << err <<"Event no. "<<n<<", matrix column no. "<<mvar;
+                *fLog << err <<" not at correct sorting position."<<endl;
+                return kFALSE;
+            }
 
         for(Int_t n=0;n<numdata-1;n++)
@@ -388 +583 @@
         fDataSort[(mvar+1)*numdata-1]=isort[numdata-1];
     }
+    return kTRUE;
 }