Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 7692)
+++ trunk/MagicSoft/Mars/Changelog	(revision 7693)
@@ -18,4 +18,34 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2006/05/05 Thomas Bretz
+
+   * mhflux/MHEnergyEst.cc:
+     - print the result to the all-logstream
+     - changed the fit range not to take the overflow into account
+
+   * mranforest/MRanForest.[h,cc]:
+     - Use the default Reset() instead of Reset(0)
+     - changed output while training
+     - replaced a lot of TArrays by MArrays
+     - at some points replaced gRandom->Rndm by gRandom->Integer
+     - removed some obsolete arguments from ModifyDataSort
+     - In CreateDataSort isort need not to be initialized, it is
+       done by TMath::Sort anyhow
+     - a small simplification to ModifyDataSort
+     - added some const-qualifiers in funciton calls
+
+   * mranforest/MRanTree.[h,cc]:
+     - replaced a lot of TArrays by MArrays
+     - removed some obsolete calls to Reset(0) after the instatization
+       of the array
+     - small acceleration of the averaging when calculating fBestSplit[k]
+     - at some points replaced gRandom->Rndm by gRandom->Integer
+     - directly give mean[kbuild] and square[kbuild] as an argument
+       to FindBestSplit
+     - removed the obsolste dereferencing from the call to FindBestSplit
+     - added some const-qualifiers in funciton calls
+
+
 
  2006/05/04 Thomas Bretz
Index: trunk/MagicSoft/Mars/mranforest/MRanForest.cc
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanForest.cc	(revision 7692)
+++ trunk/MagicSoft/Mars/mranforest/MRanForest.cc	(revision 7693)
@@ -51,4 +51,7 @@
 #include "MParList.h"
 
+#include "MArrayI.h"
+#include "MArrayF.h"
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -301,5 +304,5 @@
 
     fHadTrue.Set(numdata);
-    fHadTrue.Reset(0);
+    fHadTrue.Reset();
 
     for (Int_t j=0;j<numdata;j++)
@@ -312,5 +315,5 @@
     // setup labels for classification/regression
     fClass.Set(numdata);
-    fClass.Reset(0);
+    fClass.Reset();
 
     if (fClassify)
@@ -320,14 +323,14 @@
     // allocating and initializing arrays
     fHadEst.Set(numdata);
-    fHadEst.Reset(0);
+    fHadEst.Reset();
 
     fNTimesOutBag.Set(numdata);
-    fNTimesOutBag.Reset(0);
+    fNTimesOutBag.Reset();
 
     fDataSort.Set(dim*numdata);
-    fDataSort.Reset(0);
+    fDataSort.Reset();
 
     fDataRang.Set(dim*numdata);
-    fDataRang.Reset(0);
+    fDataRang.Reset();
 
     Bool_t useweights = fWeight.GetSize()==numdata;
@@ -438,7 +441,7 @@
 
         if(calcResolution)
-            *fLog << "no. of tree  no. of nodes  resolution in % (from oob-data -> overest. of error)" << endl;
+            *fLog << "TreeNum BagSize NumNodes TestSize res/% (from oob-data -> overest. of error)" << endl;
         else
-            *fLog << "no. of tree  no. of nodes  rms in % (from oob-data -> overest. of error)" << endl;
+            *fLog << "TreeNum BagSize NumNodes TestSize rms/% (from oob-data -> overest. of error)" << endl;
                      //        12345678901234567890123456789012345678901234567890
     }
@@ -450,7 +453,7 @@
     // bootstrap aggregating (bagging) -> sampling with replacement:
 
-    TArrayF classpopw(nclass);
-    TArrayI jinbag(numdata); // Initialization includes filling with 0
-    TArrayF winbag(numdata); // Initialization includes filling with 0
+    MArrayF classpopw(nclass);
+    MArrayI jinbag(numdata); // Initialization includes filling with 0
+    MArrayF winbag(numdata); // Initialization includes filling with 0
 
     float square=0;
@@ -463,5 +466,5 @@
         // all events in the training sample
   
-        const Int_t k = Int_t(gRandom->Rndm()*numdata);
+        const Int_t k = gRandom->Integer(numdata);
 
         if(fClassify)
@@ -472,8 +475,7 @@
         mean  +=fHadTrue[k]*fWeight[k];
         square+=fHadTrue[k]*fHadTrue[k]*fWeight[k];
-  
-        winbag[k]+=fWeight[k];
+
+        winbag[k]+=fWeight[k]; // Increase weight if chosen more than once
         jinbag[k]=1;
-
     }
 
@@ -483,10 +485,13 @@
     // In bagging procedure ca. 2/3 of all elements in the original
     // training sample are used to build the in-bag data
-    TArrayI datsortinbag=fDataSort;
-    Int_t ninbag=0;
-
-    ModifyDataSort(datsortinbag, ninbag, jinbag);
-
-    fRanTree->GrowTree(fMatrix,fHadTrue,fClass,datsortinbag,fDataRang,classpopw,mean,square,
+    const MArrayF hadtrue(fHadTrue.GetSize(), fHadTrue.GetArray());
+    const MArrayI fclass(fClass.GetSize(), fClass.GetArray());
+    const MArrayI datarang(fDataRang.GetSize(), fDataRang.GetArray());
+
+    MArrayI datsortinbag(fDataSort.GetSize(), fDataSort.GetArray());
+
+    ModifyDataSort(datsortinbag, jinbag);
+
+    fRanTree->GrowTree(fMatrix,hadtrue,fclass,datsortinbag,datarang,classpopw,mean,square,
                        jinbag,winbag,nclass);
 
@@ -501,12 +506,15 @@
     // determined from oob-data is underestimated, but can still be taken as upper limit.
 
+    Int_t ninbag = 0;
     for (Int_t ievt=0;ievt<numdata;ievt++)
     {
         if (jinbag[ievt]>0)
+        {
+            ninbag++;
             continue;
+        }
 
         fHadEst[ievt] +=fRanTree->TreeHad((*fMatrix), ievt);
         fNTimesOutBag[ievt]++;
-
     }
 
@@ -529,6 +537,8 @@
     //-------------------------------------------------------------------
     // give running output
-    *fLog << setw(5)  << fTreeNo;
-    *fLog << setw(18) << fRanTree->GetNumEndNodes();
+    *fLog << setw(4) << fTreeNo;
+    *fLog << Form(" %8.1f", 100.*ninbag/numdata);
+    *fLog << setw(9) << fRanTree->GetNumEndNodes();
+    *fLog << Form("  %9.1f", 100.*n/numdata);
     *fLog << Form("%18.2f", ferr*100.);
     *fLog << endl;
@@ -563,5 +573,5 @@
         {
             v[n]=(*fMatrix)(n,mvar);
-            isort[n]=n;
+            //isort[n]=n;
 
             if(TMath::IsNaN(v[n]))
@@ -607,10 +617,11 @@
 }
 
-void MRanForest::ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag)
+// Reoves all indices which are not in the bag from the datsortinbag
+void MRanForest::ModifyDataSort(MArrayI &datsortinbag, const MArrayI &jinbag)
 {
     const Int_t numdim=GetNumDim();
     const Int_t numdata=GetNumData();
 
-    ninbag=0;
+    Int_t ninbag=0;
     for (Int_t n=0;n<numdata;n++)
         if(jinbag[n]==1) ninbag++;
@@ -618,25 +629,24 @@
     for(Int_t m=0;m<numdim;m++)
     {
+        Int_t *subsort = &datsortinbag[m*numdata];
+
         Int_t k=0;
-        Int_t nt=0;
-        for(Int_t n=0;n<numdata;n++)
+        for(Int_t n=0;n<ninbag;n++)
         {
-            if(jinbag[datsortinbag[m*numdata+k]]==1)
+            if(jinbag[subsort[k]]==1)
             {
-                datsortinbag[m*numdata+nt]=datsortinbag[m*numdata+k];
+                subsort[n] = subsort[k];
                 k++;
             }else{
-                for(Int_t j=1;j<numdata-k;j++)
+                for(Int_t j=k+1;j<numdata;j++)
                 {
-                    if(jinbag[datsortinbag[m*numdata+k+j]]==1)
+                    if(jinbag[subsort[j]]==1)
                     {
-                        datsortinbag[m*numdata+nt]=datsortinbag[m*numdata+k+j];
-                        k+=j+1;
+                        subsort[n] = subsort[j];
+                        k = j+1;
                         break;
                     }
                 }
             }
-            nt++;
-            if(nt>=ninbag) break;
         }
     }
Index: trunk/MagicSoft/Mars/mranforest/MRanForest.h
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanForest.h	(revision 7692)
+++ trunk/MagicSoft/Mars/mranforest/MRanForest.h	(revision 7693)
@@ -21,4 +21,7 @@
 class TVector;
 class TObjArray;
+
+class MArrayI;
+class MArrayF;
 
 class MRanTree;
@@ -68,5 +71,5 @@
     // create and modify (->due to bagging) fDataSort
     Bool_t CreateDataSort();
-    void ModifyDataSort(TArrayI &datsortinbag, Int_t ninbag, const TArrayI &jinbag);
+    void ModifyDataSort(MArrayI &datsortinbag, const MArrayI &jinbag);
 
 public:
Index: trunk/MagicSoft/Mars/mranforest/MRanTree.cc
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanTree.cc	(revision 7692)
+++ trunk/MagicSoft/Mars/mranforest/MRanTree.cc	(revision 7693)
@@ -38,4 +38,7 @@
 #include <TRandom.h>
 
+#include "MArrayI.h"
+#include "MArrayF.h"
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -95,7 +98,7 @@
 }
 
-void MRanTree::GrowTree(TMatrix *mat, const TArrayF &hadtrue, const TArrayI &idclass,
-                        TArrayI &datasort, const TArrayI &datarang, TArrayF &tclasspop,
-                        float &mean, float &square, TArrayI &jinbag, const TArrayF &winbag,
+void MRanTree::GrowTree(TMatrix *mat, const MArrayF &hadtrue, const MArrayI &idclass,
+                        MArrayI &datasort, const MArrayI &datarang, MArrayF &tclasspop,
+                        const Float_t &mean, const Float_t &square, const MArrayI &jinbag, const MArrayF &winbag,
                         const int nclass)
 {
@@ -110,12 +113,12 @@
     for (Int_t n=0;n<numdata;n++) if(jinbag[n]==1) ninbag++;
 
-    TArrayI bestsplit(nrnodes);      bestsplit.Reset(0);
-    TArrayI bestsplitnext(nrnodes);  bestsplitnext.Reset(0);
-
-    fBestVar.Set(nrnodes);    fBestVar.Reset(0);
-    fTreeMap1.Set(nrnodes);   fTreeMap1.Reset(0);
-    fTreeMap2.Set(nrnodes);   fTreeMap2.Reset(0);
-    fBestSplit.Set(nrnodes);  fBestSplit.Reset(0);
-    fGiniDec.Set(numdim);     fGiniDec.Reset(0);
+    MArrayI bestsplit(nrnodes);
+    MArrayI bestsplitnext(nrnodes);
+
+    fBestVar.Set(nrnodes);    fBestVar.Reset();
+    fTreeMap1.Set(nrnodes);   fTreeMap1.Reset();
+    fTreeMap2.Set(nrnodes);   fTreeMap2.Reset();
+    fBestSplit.Set(nrnodes);  fBestSplit.Reset();
+    fGiniDec.Set(numdim);     fGiniDec.Reset();
 
 
@@ -139,7 +142,5 @@
         const Int_t &msp =fBestVar[k];
 
-        fBestSplit[k]  = (*mat)(bsp, msp);
-        fBestSplit[k] += (*mat)(bspn,msp);
-        fBestSplit[k] /= 2.;
+        fBestSplit[k] = ((*mat)(bsp, msp)+(*mat)(bspn,msp))/2;
     }
 
@@ -151,9 +152,9 @@
 }
 
-int MRanTree::FindBestSplitGini(const TArrayI &datasort,const TArrayI &datarang,
-                                const TArrayF &hadtrue,const TArrayI &idclass,
-                                Int_t ndstart,Int_t ndend, TArrayF &tclasspop,
-                                float &mean, float &square, Int_t &msplit,
-                                Float_t &decsplit,Int_t &nbest, const TArrayF &winbag,
+int MRanTree::FindBestSplitGini(const MArrayI &datasort,const MArrayI &datarang,
+                                const MArrayF &hadtrue,const MArrayI &idclass,
+                                Int_t ndstart,Int_t ndend, const MArrayF &tclasspop,
+                                const Float_t &mean, const Float_t &square, Int_t &msplit,
+                                Float_t &decsplit,Int_t &nbest, const MArrayF &winbag,
                                 const int nclass)
 {
@@ -161,6 +162,4 @@
     const Int_t numdata = (nrnodes-1)/2;
     const Int_t mdim = fGiniDec.GetSize();
-
-    TArrayF wr(nclass); wr.Reset(0);// right node
 
     // For the best split, msplit is the index of the variable (e.g Hillas par.,
@@ -202,6 +201,6 @@
         Double_t rld=0;
 
-        TArrayF wl(nclass); wl.Reset(0.);// left node //nclass
-        wr = tclasspop;
+        MArrayF wl(nclass); // left node //nclass
+        MArrayF wr(tclasspop);  // right node//nclass
 
         Double_t critvar=-1.0e20;
@@ -230,5 +229,4 @@
             const Double_t crit=(rln/rld)+(rrn/rrd);
 
-
             if (crit<=critvar) continue;
 
@@ -249,9 +247,9 @@
 }
 
-int MRanTree::FindBestSplitSigma(const TArrayI &datasort,const TArrayI &datarang,
-                                 const TArrayF &hadtrue, const TArrayI &idclass,
-                                 Int_t ndstart,Int_t ndend, TArrayF &tclasspop,
-                                 float &mean, float &square, Int_t &msplit,
-                                 Float_t &decsplit,Int_t &nbest, const TArrayF &winbag,
+int MRanTree::FindBestSplitSigma(const MArrayI &datasort,const MArrayI &datarang,
+                                 const MArrayF &hadtrue, const MArrayI &idclass,
+                                 Int_t ndstart,Int_t ndend, const MArrayF &tclasspop,
+                                 const Float_t &mean, const Float_t &square, Int_t &msplit,
+                                 Float_t &decsplit,Int_t &nbest, const MArrayF &winbag,
                                  const int nclass)
 {
@@ -259,6 +257,4 @@
     const Int_t numdata = (nrnodes-1)/2;
     const Int_t mdim = fGiniDec.GetSize();
-
-    float wr=0;// right node
 
     // For the best split, msplit is the index of the variable (e.g Hillas par., zenith angle ,...)
@@ -300,17 +296,16 @@
     for (Int_t mt=0; mt<fNumTry; mt++)
     {
-        const Int_t mvar=Int_t(gRandom->Rndm()*mdim);
+        const Int_t mvar= gRandom->Integer(mdim);
         const Int_t mn  = mvar*numdata;
 
         Double_t rrn=0, rrd=0, rln=0, rld=0;
-        Double_t esumr=0, esuml=0, e2sumr=0,e2suml=0;
-
-        esumr =mean;
-        e2sumr=square;
-        esuml =0;
-        e2suml=0;
+
+        Double_t esumr =mean;
+        Double_t e2sumr=square;
+        Double_t esuml =0;
+        Double_t e2suml=0;
 
         float wl=0.;// left node
-        wr = tclasspop[0];
+        float wr=tclasspop[0]; // right node
 
         Double_t critvar=critmin;
@@ -400,6 +395,6 @@
 }
 
-void MRanTree::MoveData(TArrayI &datasort,Int_t ndstart, Int_t ndend,
-                        TArrayI &idmove,TArrayI &ncase,Int_t msplit,
+void MRanTree::MoveData(MArrayI &datasort,Int_t ndstart, Int_t ndend,
+                        MArrayI &idmove,MArrayI &ncase,Int_t msplit,
                         Int_t nbest,Int_t &ndendl)
 {
@@ -410,5 +405,5 @@
     const Int_t mdim    = fGiniDec.GetSize();
 
-    TArrayI tdatasort(numdata); tdatasort.Reset(0);
+    MArrayI tdatasort(numdata);
 
     // compute idmove = indicator of case nos. going left
@@ -448,7 +443,7 @@
 }
 
-void MRanTree::BuildTree(TArrayI &datasort,const TArrayI &datarang, const TArrayF &hadtrue,
-                         const TArrayI &idclass, TArrayI &bestsplit, TArrayI &bestsplitnext,
-                         TArrayF &tclasspop, float &tmean, float &tsquare, const TArrayF &winbag,
+void MRanTree::BuildTree(MArrayI &datasort,const MArrayI &datarang, const MArrayF &hadtrue,
+                         const MArrayI &idclass, MArrayI &bestsplit, MArrayI &bestsplitnext,
+                         MArrayF &tclasspop, const Float_t &tmean, const Float_t &tsquare, const MArrayF &winbag,
                          Int_t ninbag, const int nclass)
 {
@@ -470,20 +465,20 @@
     const Int_t numdata = (nrnodes-1)/2;
 
-    TArrayI nodepop(nrnodes);     nodepop.Reset(0);
-    TArrayI nodestart(nrnodes);   nodestart.Reset(0);
-    TArrayI parent(nrnodes);      parent.Reset(0);
-
-    TArrayI ncase(numdata);       ncase.Reset(0);
-    TArrayI idmove(numdata);      idmove.Reset(0);
-    TArrayI iv(mdim);             iv.Reset(0);
-
-    TArrayF classpop(nrnodes*nclass);  classpop.Reset(0.);//nclass
-    TArrayI nodestatus(nrnodes);       nodestatus.Reset(0);
+    MArrayI nodepop(nrnodes);
+    MArrayI nodestart(nrnodes);
+    MArrayI parent(nrnodes);
+
+    MArrayI ncase(numdata);
+    MArrayI idmove(numdata);
+    MArrayI iv(mdim);
+
+    MArrayF classpop(nrnodes*nclass);//nclass
+    MArrayI nodestatus(nrnodes);
 
     for (Int_t j=0;j<nclass;j++)
         classpop[j*nrnodes+0]=tclasspop[j];
 
-    TArrayF mean(nrnodes);   mean.Reset(0.);
-    TArrayF square(nrnodes); square.Reset(0.);
+    MArrayF mean(nrnodes);
+    MArrayF square(nrnodes);
 
     mean[0]=tmean;
@@ -509,12 +504,9 @@
               tclasspop[j]=classpop[j*nrnodes+kbuild];
 
-          tmean=mean[kbuild];
-          tsquare=square[kbuild];
-
           Int_t msplit, nbest;
           Float_t decsplit=0;
 
-          if ((*this.*FindBestSplit)(datasort,datarang,hadtrue,idclass,ndstart,
-                                     ndend, tclasspop,tmean, tsquare,msplit,decsplit,
+          if ((this->*FindBestSplit)(datasort,datarang,hadtrue,idclass,ndstart,
+                                     ndend, tclasspop,mean[kbuild],square[kbuild],msplit,decsplit,
                                      nbest,winbag,nclass))
           {
Index: trunk/MagicSoft/Mars/mranforest/MRanTree.h
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanTree.h	(revision 7692)
+++ trunk/MagicSoft/Mars/mranforest/MRanTree.h	(revision 7693)
@@ -19,4 +19,7 @@
 class TRandom;
 
+class MArrayI;
+class MArrayF;
+
 class MRanTree : public MParContainer
 {
@@ -38,30 +41,30 @@
 
     int (MRanTree::*FindBestSplit)
-        (const TArrayI &, const TArrayI &, const TArrayF &, const TArrayI &,
-         Int_t, Int_t , TArrayF &, float &, float &, Int_t &, Float_t &,
-         Int_t &, const TArrayF &, const int); //!
+        (const MArrayI &, const MArrayI &, const MArrayF &, const MArrayI &,
+         Int_t, Int_t , const MArrayF &, const Float_t &, const Float_t &, Int_t &, Float_t &,
+         Int_t &, const MArrayF &, const int); //!
 
 
-    int FindBestSplitGini(const TArrayI &datasort, const TArrayI &datarang,
-                          const TArrayF &hadtrue, const TArrayI &idclass,
-                          Int_t ndstart, Int_t ndend, TArrayF &tclasspop,
-                          float &mean, float &square, Int_t &msplit,
-                          Float_t &decsplit, Int_t &nbest, const TArrayF &winbag,
+    int FindBestSplitGini(const MArrayI &datasort, const MArrayI &datarang,
+                          const MArrayF &hadtrue, const MArrayI &idclass,
+                          Int_t ndstart, Int_t ndend, const MArrayF &tclasspop,
+                          const Float_t &mean, const Float_t &square, Int_t &msplit,
+                          Float_t &decsplit, Int_t &nbest, const MArrayF &winbag,
                           const int nclass);
 
-    int FindBestSplitSigma(const TArrayI &datasort, const TArrayI &datarang,
-                           const TArrayF &hadtrue, const TArrayI &idclass,
-                           Int_t ndstart, Int_t ndend, TArrayF &tclasspop,
-                           float &mean, float &square, Int_t &msplit,
-                           Float_t &decsplit, Int_t &nbest, const TArrayF &winbag,
+    int FindBestSplitSigma(const MArrayI &datasort, const MArrayI &datarang,
+                           const MArrayF &hadtrue, const MArrayI &idclass,
+                           Int_t ndstart, Int_t ndend, const MArrayF &tclasspop,
+                           const Float_t &mean, const Float_t &square, Int_t &msplit,
+                           Float_t &decsplit, Int_t &nbest, const MArrayF &winbag,
                            const int nclass);
 
-    void MoveData(TArrayI &datasort, Int_t ndstart, Int_t ndend,
-                  TArrayI &idmove, TArrayI &ncase, Int_t msplit,
+    void MoveData(MArrayI &datasort, Int_t ndstart, Int_t ndend,
+                  MArrayI &idmove, MArrayI &ncase, Int_t msplit,
                   Int_t nbest, Int_t &ndendl);
 
-    void BuildTree(TArrayI &datasort, const TArrayI &datarang, const TArrayF &hadtrue,
-                   const TArrayI &idclass,TArrayI &bestsplit,TArrayI &bestsplitnext,
-                   TArrayF &tclasspop, float &tmean, float &tsquare, const TArrayF &winbag,
+    void BuildTree(MArrayI &datasort, const MArrayI &datarang, const MArrayF &hadtrue,
+                   const MArrayI &idclass,MArrayI &bestsplit,MArrayI &bestsplitnext,
+                   MArrayF &tclasspop, const Float_t &tmean, const Float_t &tsquare, const MArrayF &winbag,
                    Int_t ninbag, const int nclass);
 
@@ -93,7 +96,7 @@
 
     // functions used in tree growing process
-    void GrowTree(TMatrix *mat, const TArrayF &hadtrue, const TArrayI &idclass,
-                  TArrayI &datasort, const TArrayI &datarang,TArrayF &tclasspop,
-                  float &mean, float &square, TArrayI &jinbag, const TArrayF &winbag,
+    void GrowTree(TMatrix *mat, const MArrayF &hadtrue, const MArrayI &idclass,
+                  MArrayI &datasort, const MArrayI &datarang,MArrayF &tclasspop,
+                  const Float_t &mean, const Float_t &square, const MArrayI &jinbag, const MArrayF &winbag,
                   const int nclass);