Changeset 1489 for trunk/MagicSoft/Mars/mhist/MHMatrix.cc

Timestamp:

08/08/02 11:58:59 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhist/MHMatrix.cc (modified) (14 diffs)

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -45 +45 @@
 #include "MHMatrix.h"
 
+#include <fstream.h>
+
 #include <TList.h>
 #include <TArrayD.h>
@@ -54 +56 @@
 #include "MParList.h"
 
-#include "MDataChain.h"
+#include "MData.h"
+#include "MDataArray.h"
 
 ClassImp(MHMatrix);
 
-// --------------------------------------------------------------------------
-//
-// Default Constructor
+static const TString gsDefName  = "MHMatrix";
+static const TString gsDefTitle = "Multidimensional Matrix";
+
+// --------------------------------------------------------------------------
+//
+//  Default Constructor
 //
 MHMatrix::MHMatrix(const char *name, const char *title)
-    : fNumRow(0)
-{
-    fName  = name  ? name  : "MHMatrix";
-    fTitle = title ? title : "Multidimensional Matrix";
-
-    fData = new TList;
-    fData->SetOwner();
-
-    fRules = new TList;
-    fRules->SetOwner();
-}
-
-// --------------------------------------------------------------------------
-//
-// Destructor
+    : fNumRow(0), fData(NULL)
+{
+    fName  = name  ? name  : gsDefName.Data();
+    fTitle = title ? title : gsDefTitle.Data();
+}
+
+// --------------------------------------------------------------------------
+//
+//  Constructor. Initializes the columns of the matrix with the entries
+//  from a MDataArray
+//
+MHMatrix::MHMatrix(MDataArray *mat, const char *name, const char *title)
+    : fNumRow(0), fData(mat)
+{
+    fName  = name  ? name  : gsDefName.Data();
+    fTitle = title ? title : gsDefTitle.Data();
+}
+
+// --------------------------------------------------------------------------
+//
+//  Destructor. Does not deleted a user given MDataArray, except IsOwner
+//  was called.
 //
 MHMatrix::~MHMatrix()
 {
-    delete fData;
-    delete fRules;
+    if (TestBit(kIsOwner) && fData)
+        delete fData;
 }
 
@@ -99 +112 @@
     }
 
-    MDataChain &chain = *new MDataChain(rule);
-    if (!chain.IsValid())
-    {
-        *fLog << err << "Error - Rule cannot be translated... ignored." << endl;
-        delete &chain;
+    if (!fData)
+    {
+        fData = new MDataArray;
+        SetBit(kIsOwner);
+    }
+
+    fData->AddEntry(rule);
+}
+
+void MHMatrix::AddColumns(MDataArray *matrix)
+{
+    if (fM.IsValid())
+    {
+        *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl;
         return;
     }
-    fData->Add(&chain);
-
-    TNamed *name = new TNamed(rule, rule); // Fimxe, in 3.02/07 the title can't be "", why?
-    fRules->Add(name);
-}
-
-void MHMatrix::AddColumns(const MHMatrix *matrix)
-{
-    TIter Next(matrix->fRules);
-    TObject *rule=NULL;
-    while ((rule=Next()))
-        AddColumn(rule->GetName());
+
+    if (fData)
+        *fLog << warn << "Warning - columns already added... replacing." << endl;
+
+    if (fData && TestBit(kIsOwner))
+    {
+        delete fData;
+        ResetBit(kIsOwner);
+    }
+
+    fData = matrix;
 }
 
@@ -127 +148 @@
 Bool_t MHMatrix::SetupFill(const MParList *plist)
 {
-    if (fData->GetSize()==0)
-    {
-        *fLog << err << "Error - No Column specified... aborting." << endl;
+    if (!fData)
+    {
+        *fLog << err << "Error - No Columns initialized... aborting." << endl;
         return kFALSE;
     }
 
-    TIter Next(fData);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-        if (!data->PreProcess(plist))
-            return kFALSE;
-
-    return kTRUE;
+    return fData->PreProcess(plist);
 }
 
@@ -155 +170 @@
     if (!fM.IsValid())
     {
-        fM.ResizeTo(1, fData->GetSize());
+        fM.ResizeTo(1, fData->GetNumEntries());
         return;
     }
@@ -161 +176 @@
     TMatrix m(fM);
 
-    fM.ResizeTo(fM.GetNrows()*2, fData->GetSize());
+    fM.ResizeTo(fM.GetNrows()*2, fData->GetNumEntries());
 
     for (int x=0; x<m.GetNcols(); x++)
@@ -176 +191 @@
     AddRow();
 
-    int col=0;
-    TIter Next(fData);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-        fM(fNumRow-1, col++) = data->GetValue();
+    for (int col=0; col<fData->GetNumEntries(); col++)
+        fM(fNumRow-1, col) = (*fData)(col);
 
     return kTRUE;
@@ -195 +207 @@
     // It's not a fatal error so we don't need to stop PostProcessing...
     //
-    if (fData->GetSize()<1 || fNumRow<1)
+    if (fData->GetNumEntries()==0 || fNumRow<1)
         return kTRUE;
 
     TMatrix m(fM);
 
-    fM.ResizeTo(fNumRow, fData->GetSize());
+    fM.ResizeTo(fNumRow, fData->GetNumEntries());
 
     for (int x=0; x<fM.GetNcols(); x++)
@@ -287 +299 @@
 void MHMatrix::Print(Option_t *) const
 {
-    int n=0;
-
-    TIter Next(fData->GetSize() ? fData : fRules);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-    {
-        *fLog << all << " Column " << setw(3) << n++ << ": " << flush;
-        data->Print();
-        *fLog << endl;
-    }
+    if (fData)
+        fData->Print();
+    else
+        *fLog << all << "Sorry, no column information available." << endl;
 
     fM.Print();
@@ -336 +342 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Calculated the distance of vector evt from the reference sample
+// represented by the covariance metrix m.
+//  - If n<0 the kernel method is applied and
+//    sum(epx(-d)/n is returned.
+//  - For n>=0 the n nearest neighbors are summed and
+//    sqrt(sum(d))/n is returned.
+//  - if n==0 all distances are summed
+//
 Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const
 {
+    if (num==0) // may later be used for another method
+    {
+        TVector d = evt;
+        d *= m;
+        return sqrt(d*evt); 
+    }
+
     const Int_t rows = fM.GetNrows();
     const Int_t cols = fM.GetNcols();
@@ -358 +381 @@
 
         dists[i] = d2*d; // square of distance
+
+        //
+        // This corrects for numerical uncertanties in cases of very
+        // small distances...
+        //
+        if (dists[i]<0)
+            dists[i]=0;
     }
 
@@ -363 +393 @@
     TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);
 
-    const Int_t n = num<rows ? num : rows;
-
-    Double_t res = 0;
-    for (int i=0; i<n; i++)
-        res += dists[idx[i]];
-
-    return sqrt(res)/n;
-}
-
+    const Int_t n = abs(num)<rows ? abs(num) : rows;
+
+    if (num<0)
+    {
+        //
+        // Kernel function sum
+        //
+        const Double_t h = 1;
+
+        Double_t res = 0;
+        for (int i=0; i<n; i++)
+            res += exp(-dists[idx[i]]/h);
+
+        return log(res/n);
+    }
+    else
+    {
+        //
+        // Nearest Neighbor sum
+        //
+        Double_t res = 0;
+        for (int i=0; i<n; i++)
+            res += dists[idx[i]];
+
+        return sqrt(res/n);
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+// Calls calc dist. In the case of the first call the covariance matrix
+// fM2 is calculated.
+//  - If n<0 it is divided by (nrows-1)/h while h is the kernel factor.
+//
 Double_t MHMatrix::CalcDist(const TVector &evt, Int_t num)
 {
@@ -379 +434 @@
         fM2.ResizeTo(m);
         fM2 = m;
+        fM2 *= fM.GetNrows()-1;
         delete &m;
     }
@@ -391 +447 @@
     fM.ResizeTo(m);
     fM = m;
-
-    TIter Next(fRules);
-    TObject *obj = NULL;
-    while ((obj=Next()))
-        fData->Add(new MDataChain(obj->GetName()));
-}
+}
+
+void MHMatrix::StreamPrimitive(ofstream &out) const
+{
+    Bool_t data = fData && !TestBit(kIsOwner);
+
+    if (data)
+    {
+        fData->SavePrimitive(out);
+        out << endl;
+    }
+
+    out << "   MHMatrix " << GetUniqueName();
+
+    if (data || fName!=gsDefName || fTitle!=gsDefTitle)
+    {
+        out << "(";
+        if (data)
+            out << "&" << fData->GetUniqueName();
+        if (fName!=gsDefName || fTitle!=gsDefTitle)
+        {
+            if (data)
+                out << ", ";
+            out << "\"" << fName << "\"";
+            if (fTitle!=gsDefTitle)
+                out << ", \"" << fTitle << "\"";
+        }
+    }
+    out << ");" << endl;
+
+    if (fData && TestBit(kIsOwner))
+        for (int i=0; i<fData->GetNumEntries(); i++)
+            out << "   " << GetUniqueName() << ".AddColumn(\"" << (*fData)[i].GetRule() << "\");" << endl;
+}