Changeset 2297

Timestamp:

08/01/03 11:33:07 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhist/MHCamera.cc (modified) (4 diffs)
• mhist/MHCamera.h (modified) (3 diffs)
• mhist/MHStarMap.cc (modified) (1 diff)
• mhist/MHStarMap.h (modified) (2 diffs)
• mranforest/MRanForest.cc (modified) (13 diffs)
• mranforest/MRanForest.h (modified) (5 diffs)
• mranforest/MRanForestFill.cc (modified) (3 diffs)
• mranforest/MRanTree.cc (modified) (16 diffs)
• mranforest/MRanTree.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+ 2003/08/01: Thomas Bretz
+
+   * mhist/MHCamera.[h,cc]:
+     - added Fill(x, y, w)
+     - renamed GetStats to GetStatisticBox
+     
+   * mhist/MHStarMap.[h,cc]:
+     - include TH2 moved to source file
+
+   * mranforest/MRanForest.[h,cc], mranforest/MRanTree.[h,cc]:
+     - do not use all the data numbers and dimensions in thousands
+       of arguments if the data is available eg from the size of an array
+     - removed obsolete variables from header
+     - many small simplifications
+     - removed some obsolete variables from functions
+     - added many const qualifiers
+     - localized many more variables
+     
+   * mranforest/MRanForestFill.[h,cc]:
+     - default fNumTrees set to -1 tree (all trees)
+
+
+     
  2003/07/29: Thomas Bretz
  

trunk/MagicSoft/Mars/mhist/MHCamera.cc

@@ -222 +222 @@
 }
 
+// ------------------------------------------------------------------------
+//
+// Use x and y in millimeters
+//
+Int_t MHCamera::Fill(Axis_t x, Axis_t y, Stat_t w)
+{
+    for (Int_t idx=0; idx<fNcells-2; idx++)
+    {
+        MHexagon hex((*fGeomCam)[idx]);
+        if (hex.DistanceToPrimitive(x, y)>0)
+            continue;
+
+        SetUsed(idx);
+        return Fill(idx, w);
+    }
+    return -1;
+}
+
 Stat_t MHCamera::GetMean(Int_t axis) const
 {
@@ -839 +857 @@
 }
 
-TPaveStats *MHCamera::GetStats()
+TPaveStats *MHCamera::GetStatisticBox()
 {
     TObject *obj = 0;
@@ -857 +875 @@
 void MHCamera::UpdateLegend(Float_t min, Float_t max, Bool_t islog)
 {
-    TPaveStats *stats = GetStats();
+    TPaveStats *stats = GetStatisticBox();
 
     const Float_t hndc   = 0.92 - (stats ? stats->GetY1NDC() : 1);
@@ -992 +1010 @@
 // ------------------------------------------------------------------------
 //
-// Function introduced  (31-01-03)  WILL BE REMOVED IN THE FUTURE! DON'T
-// USE IT!
 //
 Int_t MHCamera::GetPixelIndex(Int_t px, Int_t py) const

trunk/MagicSoft/Mars/mhist/MHCamera.h

@@ -45 +45 @@
     void  SetRange();
 
-    TPaveStats *GetStats();
+    TPaveStats *GetStatisticBox();
 
     Int_t GetPixelIndex(Int_t px, Int_t py) const;
@@ -74 +74 @@
 
     Bool_t IsUsed(Int_t idx) const { return TESTBIT(const_cast<TArrayC&>(fUsed)[idx], kIsUsed); }
+
+    Int_t Fill(Axis_t x, Axis_t y, Stat_t w);
 
     //void     AddPixContent(Int_t idx) const { AddBinContent(idx+1); }
@@ -146 +148 @@
 
     //void SetStatusBar(TGStatusBar *bar) { fStatusBar = bar; }
+
+    const MGeomCam &GetGeomCam() const { return *fGeomCam; }
 
     ClassDef(MHCamera, 1) // Displays the magic camera

trunk/MagicSoft/Mars/mhist/MHStarMap.cc

@@ -32 +32 @@
 //
 /////////////////////////////////////////////////////////////////////////////
-
 #include "MHStarMap.h"
 
+#include <TH2.h>
 #include <TStyle.h>   // gStyle
 #include <TColor.h>   // SetRGB

trunk/MagicSoft/Mars/mhist/MHStarMap.h

@@ -6 +6 @@
 #endif
 
-#ifndef ROOT_TH2
-#include <TH2.h>
-#endif
-
+class TH2F;
 class MHillas;
 
@@ -35 +32 @@
     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
 
-    TH1 *GetHistByName(const TString name) { return fStarMap; }
+    TH1 *GetHistByName(const TString name) { return (TH1*)fStarMap; }
 
     TH2F *GetHist() { return fStarMap; }

trunk/MagicSoft/Mars/mranforest/MRanForest.cc

@@ -98 +98 @@
 }
 
+Int_t MRanForest::GetNumDim() const
+{
+    return fGammas ? fGammas->GetM().GetNcols() : 0;
+}
+
+
 Double_t MRanForest::CalcHadroness(const TVector &event)
 {
@@ -127 +133 @@
 }
 
+Int_t MRanForest::GetNumData() const
+{
+    return fHadrons && fGammas ? fHadrons->GetM().GetNrows()+fGammas->GetM().GetNrows() : 0;
+}
+
 Bool_t MRanForest::SetupGrow(MHMatrix *mhad,MHMatrix *mgam)
 {
@@ -134 +145 @@
 
     // determine data entries and dimension of Hillas-parameter space
-    fNumHad=fHadrons->GetM().GetNrows();
-    fNumGam=fGammas->GetM().GetNrows();
-    fNumDim=fHadrons->GetM().GetNcols();
-
-    if (fNumDim!=fHadrons->GetM().GetNcols()) return kFALSE;
-
-    fNumData=fNumHad+fNumGam;
+    //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(fNumData);
+    fHadTrue.Set(numdata);
     fHadTrue.Reset();
-    fHadEst.Set(fNumData);
+    fHadEst.Set(numdata);
 
     fPrior.Set(2);
     fClassPop.Set(2);
-    fWeight.Set(fNumData);
-    fNTimesOutBag.Set(fNumData);
+    fWeight.Set(numdata);
+    fNTimesOutBag.Set(numdata);
     fNTimesOutBag.Reset();
 
-    fDataSort.Set(fNumDim*fNumData);
-    fDataRang.Set(fNumDim*fNumData);
+    fDataSort.Set(dim*numdata);
+    fDataRang.Set(dim*numdata);
 
     // calculating class populations (= no. of gammas and hadrons)
     fClassPop.Reset();
-    for(Int_t n=0;n<fNumData;n++)
+    for(Int_t n=0;n<numdata;n++)
         fClassPop[fHadTrue[n]]++;
 
     // setting weights taking into account priors
     if (!fUsePriors)
-    {
         fWeight.Reset(1.);
-    }else{
+    else
+    {
         for(Int_t j=0;j<2;j++)
-            fPrior[j] *= (fClassPop[j]>=1) ?
-                Float_t(fNumData)/Float_t(fClassPop[j]):0;
-
-        for(Int_t n=0;n<fNumData;n++)
+            fPrior[j] *= (fClassPop[j]>=1) ? (Float_t)numdata/fClassPop[j]:0;
+
+        for(Int_t n=0;n<numdata;n++)
             fWeight[n]=fPrior[fHadTrue[n]];
     }
@@ -186 +198 @@
 
     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]++;
+    }
 }
 
@@ -206 +229 @@
     }
 
-    TArrayF classpopw(2);
-    TArrayI jinbag(fNumData); // Initialization includes filling with 0
-    TArrayF winbag(fNumData); // Initialization includes filling with 0
+    const Int_t numdata = GetNumData();
 
     // bootstrap aggregating (bagging) -> sampling with replacement:
@@ -215 +236 @@
     // {0,1,...,fNumData-1}, which is the set of the index numbers of
     // all events in the training sample
-    for (Int_t n=0; n<fNumData; n++)
-    {
-        const Int_t k = Int_t(gRandom->Rndm()*fNumData);
+    TArrayF classpopw(2);
+    TArrayI jinbag(numdata); // Initialization includes filling with 0
+    TArrayF winbag(numdata); // Initialization includes filling with 0
+
+    for (Int_t n=0; n<numdata; n++)
+    {
+        const Int_t k = Int_t(gRandom->Rndm()*numdata);
 
         classpopw[fHadTrue[k]]+=fWeight[k];
@@ -237 +262 @@
 
     // growing single tree
-    fRanTree->GrowTree(hadrons,gammas,fNumData,fNumDim,fHadTrue,datsortinbag,
-                       fDataRang,classpopw,jinbag,winbag,fWeight);
+    fRanTree->GrowTree(hadrons,gammas,fHadTrue,datsortinbag,
+                       fDataRang,classpopw,jinbag,winbag);
 
     // error-estimates from out-of-bag data (oob data):
@@ -249 +274 @@
     // determined from oob-data is underestimated, but can still be taken as upper limit.
 
-    for (Int_t ievt=0;ievt<fNumHad;ievt++)
+    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)
@@ -256 +285 @@
         fNTimesOutBag[ievt]++;
     }
-    for (Int_t ievt=0;ievt<fNumGam;ievt++)
-    {
-        if (jinbag[fNumHad+ievt]>0)
+
+    for (Int_t ievt=numhad;ievt<numdata;ievt++)
+    {
+        if (jinbag[ievt]>0)
             continue;
-        fHadEst[fNumHad+ievt] += fRanTree->TreeHad(gammas, ievt);
-        fNTimesOutBag[fNumHad+ievt]++;
-    }
-
+        fHadEst[ievt] += fRanTree->TreeHad(gammas, ievt-numhad);
+        fNTimesOutBag[ievt]++;
+    }
+    */
     Int_t n=0;
     Double_t ferr=0;
-    for (Int_t ievt=0;ievt<fNumData;ievt++)
+    for (Int_t ievt=0;ievt<numdata;ievt++)
         if (fNTimesOutBag[ievt]!=0)
         {
@@ -301 +331 @@
     // 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).
-
-    TArrayF v(fNumData);
-    TArrayI isort(fNumData);
+    const Int_t numdata = GetNumData();
+
+    TArrayF v(numdata);
+    TArrayI isort(numdata);
 
     const TMatrix &hadrons=fHadrons->GetM();
     const TMatrix &gammas=fGammas->GetM();
 
-    for (Int_t j=0;j<fNumHad;j++)
+    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<fNumGam;j++)
-        fHadTrue[j+fNumHad]=0;
-
-    for (Int_t mvar=0;mvar<fNumDim;mvar++)
-    {
-        for(Int_t n=0;n<fNumHad;n++)
+    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++)
         {
             v[n]=hadrons(n,mvar);
@@ -322 +357 @@
         }
 
-        for(Int_t n=0;n<fNumGam;n++)
+        for(Int_t n=0;n<numgam;n++)
         {
-            v[n+fNumHad]=gammas(n,mvar);
-            isort[n+fNumHad]=n;
+            v[n+numhad]=gammas(n,mvar);
+            isort[n+numhad]=n;
         }
 
-        TMath::Sort(fNumData,v.GetArray(),isort.GetArray(),kFALSE);
+        TMath::Sort(numdata,v.GetArray(),isort.GetArray(),kFALSE);
 
         // this sorts the v[n] in ascending order. isort[n] is the event number
@@ -334 +369 @@
         // event numbers are 0,1,...).
 
-        for(Int_t n=0;n<fNumData-1;n++)
+        for(Int_t n=0;n<numdata-1;n++)
         {
             const Int_t n1=isort[n];
             const Int_t n2=isort[n+1];
 
-            fDataSort[mvar*fNumData+n]=n1;
-            if(n==0) fDataRang[mvar*fNumData+n1]=0;
+            fDataSort[mvar*numdata+n]=n1;
+            if(n==0) fDataRang[mvar*numdata+n1]=0;
             if(v[n]<v[n+1])
             {
-                fDataRang[mvar*fNumData+n2]=fDataRang[mvar*fNumData+n1]+1;
+                fDataRang[mvar*numdata+n2]=fDataRang[mvar*numdata+n1]+1;
             }else{
-                fDataRang[mvar*fNumData+n2]=fDataRang[mvar*fNumData+n1];
+                fDataRang[mvar*numdata+n2]=fDataRang[mvar*numdata+n1];
             }
         }
-        fDataSort[(mvar+1)*fNumData-1]=isort[fNumData-1];
+        fDataSort[(mvar+1)*numdata-1]=isort[numdata-1];
     }
 }
@@ -354 +389 @@
 void MRanForest::ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag)
 {
+    const Int_t numdim=GetNumDim();
+    const Int_t numdata=GetNumData();
+
     ninbag=0;
-    for (Int_t n=0;n<fNumData;n++)
+    for (Int_t n=0;n<numdata;n++)
         if(jinbag[n]==1) ninbag++;
 
-    for(Int_t m=0;m<fNumDim;m++)
+    for(Int_t m=0;m<numdim;m++)
     {
         Int_t k=0;
         Int_t nt=0;
-        for(Int_t n=0;n<fNumData;n++)
+        for(Int_t n=0;n<numdata;n++)
         {
-            if(jinbag[datsortinbag[m*fNumData+k]]==1)
+            if(jinbag[datsortinbag[m*numdata+k]]==1)
             {
-                datsortinbag[m*fNumData+nt]=datsortinbag[m*fNumData+k];
+                datsortinbag[m*numdata+nt]=datsortinbag[m*numdata+k];
                 k++;
             }else{
-                for(Int_t j=1;j<fNumData-k;j++)
+                for(Int_t j=1;j<numdata-k;j++)
                 {
-                    if(jinbag[datsortinbag[m*fNumData+k+j]]==1)
+                    if(jinbag[datsortinbag[m*numdata+k+j]]==1)
                     {
-                        datsortinbag[m*fNumData+nt]=datsortinbag[m*fNumData+k+j];
+                        datsortinbag[m*numdata+nt]=datsortinbag[m*numdata+k+j];
                         k+=j+1;
                         break;

trunk/MagicSoft/Mars/mranforest/MRanForest.h

@@ -29 +29 @@
 class MRanTree;
 class TVector;
+class TMatrix;
 
 class MRanForest : public MParContainer
@@ -42 +43 @@
     MHMatrix *fGammas;
     MHMatrix *fHadrons;
-
-    Int_t   fNumGam;
-    Int_t   fNumHad;
-    Int_t   fNumData;
-    Int_t   fNumDim;
 
     // true  and estimated hadronness
@@ -65 +61 @@
     // estimates for classification error of growing forest
     TArrayD fTreeHad;
+
+    void InitHadEst(Int_t from, Int_t to, const TMatrix &m, TArrayI &jinbag);
 
 protected:
@@ -90 +88 @@
     MRanTree *GetCurTree() { return fRanTree; }
     Int_t      GetNumTrees() const { return fNumTrees; }
-    Int_t      GetNumData() const { return fNumData; }
-    Int_t      GetNumDim() const { return fNumDim; }
+    Int_t      GetNumData() const;
+    Int_t      GetNumDim() const;
     Double_t   GetTreeHad(Int_t i) const { return fTreeHad.At(i); }
  
@@ -99 +97 @@
     Bool_t AsciiWrite(ostream &out) const;
 
-    ClassDef(MRanForest, 1) // Storage container for tree structure
+    ClassDef(MRanForest, 0) // Storage container for tree structure
 };
 

trunk/MagicSoft/Mars/mranforest/MRanForestFill.cc

@@ -53 +53 @@
 //
 //
-MRanForestFill::MRanForestFill(const char *name, const char *title):fNumTrees(100)
+MRanForestFill::MRanForestFill(const char *name, const char *title):fNumTrees(-1)
 {
     //
@@ -98 +98 @@
 Int_t MRanForestFill::Process()
 {
-    fNum++;
-    if(!(fRanForest->AddTree(fRanTree)))
+    if (!(fRanForest->AddTree(fRanTree)))
         return kFALSE;
 
-    return fNum<fNumTrees;
+    fNum++;
+
+    return fNumTrees<0 ? kTRUE : fNum<fNumTrees;
 }
 
@@ -108 +109 @@
 {
     fRanForest->SetNumTrees(fNum);
-
     return kTRUE;
 }

trunk/MagicSoft/Mars/mranforest/MRanTree.cc

@@ -76 +76 @@
 }
 
-void MRanTree::GrowTree(const TMatrix &mhad, const TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,
-                        TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag,
-                        TArrayF &winbag,TArrayF &weight)
+void MRanTree::GrowTree(const TMatrix &mhad, const TMatrix &mgam,
+                        const TArrayI &hadtrue, TArrayI &datasort,
+                        const TArrayI &datarang, TArrayF &tclasspop, TArrayI &jinbag,
+                        const TArrayF &winbag)
 {
     // arrays have to be initialized with generous size, so number of total nodes (nrnodes)
     // is estimated for worst case
-    Int_t nrnodes=2*numdata+1;
+    const Int_t numdim =mhad.GetNcols();
+    const Int_t numdata=winbag.GetSize();
+    const Int_t nrnodes=2*numdata+1;
 
     // number of events in bootstrap sample
@@ -91 +94 @@
     TArrayI bestsplit(nrnodes);
     TArrayI bestsplitnext(nrnodes);
-    TArrayI nodepop(nrnodes);
-    TArrayI parent(nrnodes);
-    TArrayI nodex(numdata);
-    TArrayI nodestart(nrnodes);
-
-    TArrayI ncase(numdata);
-    TArrayI iv(numdim);
-    TArrayI idmove(numdata);
 
     fBestVar.Set(nrnodes);
@@ -113 +108 @@
 
     // tree growing
-    BuildTree(datasort,datarang,hadtrue,numdim,numdata,bestsplit,
-              bestsplitnext,nodepop,nodestart,tclasspop,nrnodes,
-              idmove,ncase,parent,jinbag,iv,winbag,ninbag);
+    BuildTree(datasort,datarang,hadtrue,bestsplit,
+              bestsplitnext,tclasspop,winbag,ninbag);
 
     // post processing, determine cut (or split) values fBestSplit
@@ -141 +135 @@
 }
 
-Int_t MRanTree::FindBestSplit(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim,
-                             Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop,
-                             Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase,
-                             TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t kbuild)
+Int_t MRanTree::FindBestSplit(const TArrayI &datasort,const TArrayI &datarang,
+                              const TArrayI &hadtrue,Int_t ndstart,Int_t ndend,TArrayF &tclasspop,
+                              Int_t &msplit,Float_t &decsplit,Int_t &nbest,
+                              const TArrayF &winbag)
 {
     if(!gRandom)
@@ -151 +145 @@
         return kFALSE;
     }
+
+    const Int_t nrnodes = fBestSplit.GetSize();
+    const Int_t numdata = (nrnodes-1)/2;
+    const Int_t mdim = fGiniDec.GetSize();
 
     // weighted class populations after split
@@ -161 +159 @@
     // and nsplitnext is the case number of the next larger value of msplit.
 
-    Int_t nc,nbestvar=0,k;
-    Float_t crit;
-    Float_t rrn, rrd, rln, rld, u;
+    Int_t nbestvar=0;
 
     // compute initial values of numerator and denominator of Gini-index,
     // Gini index= pno/dno
-    Float_t pno=0;
-    Float_t pdo=0;
-
+    Double_t pno=0;
+    Double_t pdo=0;
     for (Int_t j=0;j<2;j++)
     {
@@ -190 +185 @@
 
         // Gini index = rrn/rrd+rln/rld
-        rrn=pno;
-        rrd=pdo;
-        rln=0;
-        rld=0;
+        Double_t rrn=pno;
+        Double_t rrd=pdo;
+        Double_t rln=0;
+        Double_t rld=0;
 
         TArrayF wl(2); // left node
@@ -202 +197 @@
         for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++)
         {
-            nc=datasort[mvar*numdata+nsp];
-
-            u=winbag[nc];
-            k=hadtrue[nc];
-
-            rln=rln+u*(2*wl[k]+u);
-            rrn=rrn+u*(-2*wr[k]+u);
-            rld=rld+u;
-            rrd=rrd-u;
-
-            wl[k]=wl[k]+u;
-            wr[k]=wr[k]-u;
-
-            if (datarang[mvar*numdata+nc]<datarang[mvar*numdata+datasort[mvar*numdata+nsp+1]])
-            {
-                if (TMath::Min(rrd,rld)>1.0e-5)
-                {
-                    crit=(rln/rld)+(rrn/rrd);
-                    if (crit>critvar)
-                    {
-                        nbestvar=nsp;
-                        critvar=crit;
-                    }
-                }
-            }
+            const Int_t &nc=datasort[mvar*numdata+nsp];
+            const Int_t &k=hadtrue[nc];
+
+            const Float_t &u=winbag[nc];
+
+            rln+=u*(2*wl[k]+u);
+            rrn+=u*(-2*wr[k]+u);
+            rld+=u;
+            rrd-=u;
+
+            wl[k]+=u;
+            wr[k]-=u;
+
+            if (datarang[mvar*numdata+nc]>=datarang[mvar*numdata+datasort[mvar*numdata+nsp+1]])
+                continue;
+            if (TMath::Min(rrd,rld)<=1.0e-5)
+                continue;
+
+            const Double_t crit=(rln/rld)+(rrn/rrd);
+            if (crit<=critvar)
+                continue;
+
+            nbestvar=nsp;
+            critvar=crit;
         }
 
-        if (critvar>critmax) {
-            msplit=mvar;
-            nbest=nbestvar;
-            critmax=critvar;
-        }
+        if (critvar<=critmax)
+            continue;
+
+        msplit=mvar;
+        nbest=nbestvar;
+        critmax=critvar;
     }
 
     decsplit=critmax-crit0;
-    if (critmax<-1.0e10) jstat=1;
-
-    return jstat;
-}
-
-void MRanTree::MoveData(TArrayI &datasort,Int_t mdim,Int_t numdata,Int_t ndstart,
+
+    return critmax<-1.0e10 ? 1 : jstat;
+}
+
+void MRanTree::MoveData(TArrayI &datasort,Int_t ndstart,
                         Int_t ndend,TArrayI &idmove,TArrayI &ncase,Int_t msplit,
                         Int_t nbest,Int_t &ndendl)
@@ -249 +243 @@
     // the data in the part of datasort corresponding to the current node is moved to the
     // left if it belongs to the left child and right if it belongs to the right child-node.
-
-    Int_t nc,k,ih;
+    const Int_t numdata = ncase.GetSize();
+    const Int_t mdim    = fGiniDec.GetSize();
+
     TArrayI tdatasort(numdata);
 
     // compute idmove = indicator of case nos. going left
 
-    for (Int_t nsp=ndstart;nsp<=nbest;nsp++)
-    {
-        nc=datasort[msplit*numdata+nsp];
-        idmove[nc]=1;
-    }
-    for (Int_t nsp=nbest+1;nsp<=ndend;nsp++)
-    {
-        nc=datasort[msplit*numdata+nsp];
-        idmove[nc]=0;
+    for (Int_t nsp=ndstart;nsp<=ndend;nsp++)
+    {
+        const Int_t &nc=datasort[msplit*numdata+nsp];
+        idmove[nc]= nsp<=nbest?1:0;
     }
     ndendl=nbest;
@@ -271 +261 @@
     for(Int_t msh=0;msh<mdim;msh++)
     {
-        k=ndstart-1;
+        Int_t k=ndstart-1;
         for (Int_t n=ndstart;n<=ndend;n++)
         {
-            ih=datasort[msh*numdata+n];
-            if (idmove[ih]==1) {
-                k++;
-                tdatasort[k]=datasort[msh*numdata+n];
-            }
+            const Int_t &ih=datasort[msh*numdata+n];
+            if (idmove[ih]==1)
+                tdatasort[++k]=datasort[msh*numdata+n];
         }
 
         for (Int_t n=ndstart;n<=ndend;n++)
         {
-            ih=datasort[msh*numdata+n];
-            if (idmove[ih]==0){
-                k++;
-                tdatasort[k]=datasort[msh*numdata+n];
-            }
+            const Int_t &ih=datasort[msh*numdata+n];
+            if (idmove[ih]==0)
+                tdatasort[++k]=datasort[msh*numdata+n];
         }
-        for(Int_t k=ndstart;k<=ndend;k++)
-            datasort[msh*numdata+k]=tdatasort[k];
+
+        for(Int_t m=ndstart;m<=ndend;m++)
+            datasort[msh*numdata+m]=tdatasort[m];
     }
 
@@ -299 +286 @@
 }
 
-void MRanTree::BuildTree(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim,
-                         Int_t numdata,TArrayI &bestsplit,TArrayI &bestsplitnext,
-                         TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop,
-                         Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent,
-                         TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag)
+void MRanTree::BuildTree(TArrayI &datasort,const TArrayI &datarang,
+                         const TArrayI &hadtrue, TArrayI &bestsplit,
+                         TArrayI &bestsplitnext, TArrayF &tclasspop,
+                         const TArrayF &winbag, Int_t ninbag)
 {
     // Buildtree consists of repeated calls to two void functions, FindBestSplit and MoveData.
@@ -318 +304 @@
     // the next node to be split is numbered k+1.  When no more nodes can be split, buildtree
     // returns.
-
-    Int_t msplit, nbest, ndendl;
-    Float_t decsplit=0;
-    TArrayF classpop(2*nrnodes);
+    const Int_t mdim    = fGiniDec.GetSize();
+    const Int_t nrnodes = fBestSplit.GetSize();
+    const Int_t numdata = (nrnodes-1)/2;
+
+    TArrayI nodepop(nrnodes);
+    TArrayI nodestart(nrnodes);
+    TArrayI parent(nrnodes);
+
+    TArrayI ncase(numdata);
+    TArrayI idmove(numdata);
+    TArrayI iv(mdim);
+
+    TArrayF classpop(nrnodes*2);
     TArrayI nodestatus(nrnodes);
-
-    nodestart.Reset();
-    nodepop.Reset();
 
     for (Int_t j=0;j<2;j++)
@@ -331 +323 @@
 
     Int_t ncur=0;
-    nodestart[0]=0;
     nodepop[0]=ninbag;
     nodestatus[0]=2;
 
     // start main loop
-    for (Int_t kbuild=0;kbuild<nrnodes;kbuild++)
+    for (Int_t kbuild=0; kbuild<nrnodes; kbuild++)
     {
           if (kbuild>ncur) break;
@@ -346 +337 @@
           const Int_t ndend=ndstart+nodepop[kbuild]-1;
           for (Int_t j=0;j<2;j++)
-            tclasspop[j]=classpop[j*nrnodes+kbuild];
-
-          const Int_t jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata,
-                                          ndstart,ndend,tclasspop,msplit,decsplit,
-                                          nbest,ncase,jinbag,iv,winbag,kbuild);
-
-          if(jstat==1) {
+              tclasspop[j]=classpop[j*nrnodes+kbuild];
+
+          Int_t msplit, nbest;
+          Float_t decsplit=0;
+          const Int_t jstat=FindBestSplit(datasort,datarang,hadtrue,
+                                          ndstart,ndend,tclasspop,msplit,
+                                          decsplit,nbest,winbag);
+
+          if (jstat==1)
+          {
               nodestatus[kbuild]=-1;
               continue;
-          }else{
-              fBestVar[kbuild]=msplit;
-              fGiniDec[msplit]+=decsplit;
-
-              bestsplit[kbuild]=datasort[msplit*numdata+nbest];
-              bestsplitnext[kbuild]=datasort[msplit*numdata+nbest+1];
           }
 
-          MoveData(datasort,mdim,numdata,ndstart,ndend,idmove,ncase,
+          fBestVar[kbuild]=msplit;
+          fGiniDec[msplit]+=decsplit;
+
+          bestsplit[kbuild]=datasort[msplit*numdata+nbest];
+          bestsplitnext[kbuild]=datasort[msplit*numdata+nbest+1];
+
+          Int_t ndendl;
+          MoveData(datasort,ndstart,ndend,idmove,ncase,
                    msplit,nbest,ndendl);
 
@@ -377 +372 @@
           for (Int_t n=ndstart;n<=ndendl;n++)
           {
-              const Int_t nc=ncase[n];
-              const Int_t j=hadtrue[nc];
+              const Int_t &nc=ncase[n];
+              const Int_t &j=hadtrue[nc];
               classpop[j*nrnodes+ncur+1]+=winbag[nc];
           }
@@ -384 +379 @@
           for (Int_t n=ndendl+1;n<=ndend;n++)
           {
-              const Int_t nc=ncase[n];
-              const Int_t j=hadtrue[nc];
+              const Int_t &nc=ncase[n];
+              const Int_t &j=hadtrue[nc];
               classpop[j*nrnodes+ncur+2]+=winbag[nc];
           }

trunk/MagicSoft/Mars/mranforest/MRanTree.h

@@ -61 +61 @@
 
     // functions used in tree growing process
-    void GrowTree(const TMatrix &mhad, const TMatrix &mgam,Int_t numdata,
-                  Int_t numdim,TArrayI &hadtrue,
-                  TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag,
-                  TArrayF &winbag,TArrayF &weight);
+    void GrowTree(const TMatrix &mhad, const TMatrix &mgam,
+                  const TArrayI &hadtrue, TArrayI &datasort,
+                  const TArrayI &datarang,
+                  TArrayF &tclasspop, TArrayI &jinbag, const TArrayF &winbag);
 
-    Int_t FindBestSplit(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim,
-                        Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop,
-                        Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase,
-                        TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t kbuild);
+    Int_t FindBestSplit(const TArrayI &datasort, const TArrayI &datarang,
+                        const TArrayI &hadtrue,
+                        Int_t ndstart, Int_t ndend, TArrayF &tclasspop,
+                        Int_t &msplit, Float_t &decsplit, Int_t &nbest,
+                        const TArrayF &winbag);
 
-    void MoveData(TArrayI &datasort,Int_t mdim,Int_t numdata,Int_t ndstart,
-                  Int_t ndend,TArrayI &idmove,TArrayI &ncase,Int_t msplit,
-                  Int_t nbest,Int_t &ndendl);
+    void MoveData(TArrayI &datasort, Int_t ndstart, Int_t ndend,
+                  TArrayI &idmove, TArrayI &ncase, Int_t msplit,
+                  Int_t nbest, Int_t &ndendl);
 
-    void BuildTree(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim,
-                   Int_t numdata,TArrayI &bestsplit,TArrayI &bestsplitnext,
-                   TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop,
-                   Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent,
-                   TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag);
+    void BuildTree(TArrayI &datasort, const TArrayI &datarang,
+                   const TArrayI &hadtrue,
+                   TArrayI &bestsplit,TArrayI &bestsplitnext,
+                   TArrayF &tclasspop,
+                   const TArrayF &winbag,
+                   Int_t ninbag);
 
     Double_t TreeHad(const TVector &event);