Changeset 2296 for trunk/MagicSoft/Mars/mranforest/MRanTree.cc

Timestamp:

07/29/03 13:18:57 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mranforest/MRanTree.cc (modified) (17 diffs)

trunk/MagicSoft/Mars/mranforest/MRanTree.cc

@@ -76 +76 @@
 }
 
-void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,
+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)
@@ -88 +88 @@
     for (Int_t n=0;n<numdata;n++)
         if(jinbag[n]==1) ninbag++;
-
-    // weighted class populations after split
-    TArrayF wl(2); // left node
-    TArrayF wc(2); 
-    TArrayF wr(2); // right node
-    TArrayI nc(2);
 
     TArrayI bestsplit(nrnodes);
@@ -106 +100 @@
     TArrayI idmove(numdata);
 
-    idmove.Reset();
-
     fBestVar.Set(nrnodes);
     fTreeMap1.Set(nrnodes);
@@ -123 +115 @@
     BuildTree(datasort,datarang,hadtrue,numdim,numdata,bestsplit,
               bestsplitnext,nodepop,nodestart,tclasspop,nrnodes,
-              idmove,ncase,parent,jinbag,iv,winbag,wr,wc,wl,ninbag);
+              idmove,ncase,parent,jinbag,iv,winbag,ninbag);
 
     // post processing, determine cut (or split) values fBestSplit
     Int_t nhad=mhad.GetNrows();
 
-    for(Int_t k=0;k<nrnodes;k++)
-    {
-        Int_t bsp=bestsplit[k];
-        Int_t bspn=bestsplitnext[k];
-        Int_t msp=fBestVar[k];
-
-        if (GetNodeStatus(k)!=-1)
-        {
-            fBestSplit[k] = bsp<nhad ? mhad(bsp,msp):mgam(bsp-nhad,msp);
-            fBestSplit[k] += bspn<nhad ? mhad(bspn,msp):mgam(bspn-nhad,msp);
-            fBestSplit[k] /=2.;
-        }
+    for(Int_t k=0; k<nrnodes; k++)
+    {
+        if (GetNodeStatus(k)==-1)
+            continue;
+
+        const Int_t &bsp =bestsplit[k];
+        const Int_t &bspn=bestsplitnext[k];
+        const Int_t &msp =fBestVar[k];
+
+        fBestSplit[k]  = bsp<nhad  ? mhad(bsp, msp):mgam(bsp-nhad, msp);
+        fBestSplit[k] += bspn<nhad ? mhad(bspn,msp):mgam(bspn-nhad,msp);
+        fBestSplit[k] /= 2;
     }
 
@@ -152 +144 @@
                              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,TArrayF &wr,
-                             TArrayF &wc,TArrayF &wl,Int_t kbuild)
-{
+                             TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t kbuild)
+{
+    if(!gRandom)
+    {
+        *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
+        return kFALSE;
+    }
+
+    // weighted class populations after split
+    TArrayF wc(2); 
+    TArrayF wr(2); // right node
+
     // For the best split, msplit is the index of the variable (e.g Hillas par., zenith angle ,...)
     // split on. decsplit is the decreae in impurity measured by Gini-index.
@@ -160 +161 @@
     // and nsplitnext is the case number of the next larger value of msplit.
 
-    Int_t mvar,nc,nbestvar=0,jstat,k;
-    Float_t crit,crit0,critmax,critvar=0;
+    Int_t nc,nbestvar=0,k;
+    Float_t crit;
     Float_t rrn, rrd, rln, rld, u;
 
@@ -171 +172 @@
     for (Int_t j=0;j<2;j++)
     {
-          pno+=tclasspop[j]*tclasspop[j];
-          pdo+=tclasspop[j];
-    }
-    crit0=pno/pdo;
-    jstat=0;
+        pno+=tclasspop[j]*tclasspop[j];
+        pdo+=tclasspop[j];
+    }
+
+    const Double_t crit0=pno/pdo;
+    Int_t jstat=0;
 
     // start main loop through variables to find best split,
     // (Gini-index as criterium crit)
 
-    critmax=-1.0e20;  // FIXME: Replace by a constant from limits.h
+    Double_t critmax=-1.0e20;  // FIXME: Replace by a constant from limits.h
 
     // random split selection, number of trials = fNumTry
-    if(!gRandom)
-    {
-        *fLog << err << dbginf << "gRandom not initialized... aborting." << endl;
-        return kFALSE;
-    }
     for(Int_t mt=0;mt<fNumTry;mt++)
     {
-        mvar=Int_t(mdim*gRandom->Rndm());
+        Int_t mvar=Int_t(gRandom->Rndm()*mdim);
 
         // Gini index = rrn/rrd+rln/rld
@@ -197 +194 @@
         rln=0;
         rld=0;
-        wl.Reset();
-
-        for (Int_t j=0;j<2;j++)
-        {
-            wr[j]=tclasspop[j];
-        }
-
-        critvar=-1.0e20;
+
+        TArrayF wl(2); // left node
+        wr = tclasspop;
+
+        Double_t critvar=-1.0e20;
 
         for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++)
@@ -223 +217 @@
             if (datarang[mvar*numdata+nc]<datarang[mvar*numdata+datasort[mvar*numdata+nsp+1]])
             {
-                if(TMath::Min(rrd,rld)>1.0e-5)
+                if (TMath::Min(rrd,rld)>1.0e-5)
                 {
                     crit=(rln/rld)+(rrn/rrd);
@@ -309 +303 @@
                          TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop,
                          Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent,
-                         TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc,
-                         TArrayF &wl,Int_t ninbag)
+                         TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,Int_t ninbag)
 {
     // Buildtree consists of repeated calls to two void functions, FindBestSplit and MoveData.
@@ -326 +319 @@
     // returns.
 
-    Int_t msplit,nbest,ndendl,nc,jstat,ndend,ndstart;
+    Int_t msplit, nbest, ndendl;
     Float_t decsplit=0;
-    Float_t popt1,popt2,pp;
-    TArrayF classpop;
-    TArrayI nodestatus;
-
-    nodestatus.Set(nrnodes);
-    classpop.Set(2*nrnodes);
-
-    nodestatus.Reset();
+    TArrayF classpop(2*nrnodes);
+    TArrayI nodestatus(nrnodes);
+
     nodestart.Reset();
     nodepop.Reset();
-    classpop.Reset();
-
 
     for (Int_t j=0;j<2;j++)
@@ -357 +343 @@
           // initialize for next call to FindBestSplit
 
-          ndstart=nodestart[kbuild];
-          ndend=ndstart+nodepop[kbuild]-1;
+          const Int_t ndstart=nodestart[kbuild];
+          const Int_t ndend=ndstart+nodepop[kbuild]-1;
           for (Int_t j=0;j<2;j++)
             tclasspop[j]=classpop[j*nrnodes+kbuild];
 
-          jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata,
-                              ndstart,ndend,tclasspop,msplit,decsplit,
-                              nbest,ncase,jinbag,iv,winbag,wr,wc,wl,
-                              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) {
@@ -392 +377 @@
           for (Int_t n=ndstart;n<=ndendl;n++)
           {
-              nc=ncase[n];
-              Int_t j=hadtrue[nc];
+              const Int_t nc=ncase[n];
+              const Int_t j=hadtrue[nc];
               classpop[j*nrnodes+ncur+1]+=winbag[nc];
           }
@@ -399 +384 @@
           for (Int_t n=ndendl+1;n<=ndend;n++)
           {
-              nc=ncase[n];
-              Int_t j=hadtrue[nc];
+              const Int_t nc=ncase[n];
+              const Int_t j=hadtrue[nc];
               classpop[j*nrnodes+ncur+2]+=winbag[nc];
           }
@@ -410 +395 @@
           if (nodepop[ncur+1]<=fNdSize) nodestatus[ncur+1]=-1;
           if (nodepop[ncur+2]<=fNdSize) nodestatus[ncur+2]=-1;
-          popt1=0;
-          popt2=0;
+
+          Double_t popt1=0;
+          Double_t popt2=0;
           for (Int_t j=0;j<2;j++)
           {
-            popt1+=classpop[j*nrnodes+ncur+1];
-            popt2+=classpop[j*nrnodes+ncur+2];
+              popt1+=classpop[j*nrnodes+ncur+1];
+              popt2+=classpop[j*nrnodes+ncur+2];
           }
 
           for (Int_t j=0;j<2;j++)
           {
-            if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1;
-            if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1;
+              if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1;
+              if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1;
           }
 
@@ -446 +432 @@
         {
             fNumEndNodes++;
-            pp=0;
+            Double_t pp=0;
             for (Int_t j=0;j<2;j++)
             {
@@ -482 +468 @@
 
         const Int_t m=fBestVar[kt];
-
         kt = event(m)<=fBestSplit[kt] ? fTreeMap1[kt] : fTreeMap2[kt];
     }
 
     return fBestSplit[kt];
+}
+
+Double_t MRanTree::TreeHad(const TMatrixRow &event)
+{
+    Int_t kt=0;
+    // to optimize on storage space node status and node class
+    // are coded into fBestVar:
+    // status of node kt = TMath::Sign(1,fBestVar[kt])
+    // class  of node kt = fBestVar[kt]+2 (class defined by larger
+    //  node population, actually not used)
+    // hadronness assigned to node kt = fBestSplit[kt]
+
+    for (Int_t k=0;k<fNumNodes;k++)
+    {
+        if (fBestVar[kt]<0)
+            break;
+
+        const Int_t m=fBestVar[kt];
+        kt = event(m)<=fBestSplit[kt] ? fTreeMap1[kt] : fTreeMap2[kt];
+    }
+
+    return fBestSplit[kt];
+}
+
+Double_t MRanTree::TreeHad(const TMatrix &m, Int_t ievt)
+{
+    return TreeHad(TMatrixRow(m, ievt));
 }