source: trunk/MagicSoft/Mars/mhist/MHMatrix.cc@ 1548

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MHMatrix
//
// This is a histogram container which holds a matrix with one column per
// data variable. The data variable can be a complex rule (MDataChain).
// Each event for wich Fill is called (by MFillH) is added as a new
// row to the matrix.
//
// For example:
//   MHMatrix m;
//   m.AddColumn("MHillas.fSize");
//   m.AddColumn("MMcEvt.fImpact/100");
//   m.AddColumn("HillasSource.fDist*MGeomCam.fConvMm2Deg");
//   MFillH fillm(&m);
//   taskliost.AddToList(&fillm);
//   [...]
//   m.Print();
//
/////////////////////////////////////////////////////////////////////////////
#include "MHMatrix.h"
#include <math.h> // for Alphas

#include <fstream.h>

#include <TList.h>
#include <TArrayF.h>
#include <TArrayD.h>
#include <TArrayI.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MFillH.h"
#include "MEvtLoop.h"
#include "MParList.h"
#include "MTaskList.h"

#include "MData.h"
#include "MDataArray.h"

ClassImp(MHMatrix);

const TString MHMatrix::gsDefName  = "MHMatrix";
const TString MHMatrix::gsDefTitle = "Multidimensional Matrix";

// --------------------------------------------------------------------------
//
//  Default Constructor
//
MHMatrix::MHMatrix(const char *name, const char *title)
    : 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()
{
    if (TestBit(kIsOwner) && fData)
        delete fData;
}

// --------------------------------------------------------------------------
//
// Add a new column to the matrix. This can only be done before the first
// event (row) was filled into the matrix. For the syntax of the rule
// see MDataChain.
// Returns the index of the new column, -1 in case of failure.
// (0, 1, 2, ... for the 1st, 2nd, 3rd, ...)
//
Int_t MHMatrix::AddColumn(const char *rule)
{
    if (fM.IsValid())
    {
        *fLog << warn << "Warning - matrix is already in use. Can't add a new column... skipped." << endl;
        return -1;
    }

    if (TestBit(kIsLocked))
    {
        *fLog << warn << "Warning - matrix is locked. Can't add new column... skipped." << endl;
        return -1;
    }

    if (!fData)
    {
        fData = new MDataArray;
        SetBit(kIsOwner);
    }

    fData->AddEntry(rule);
    return fData->GetNumEntries()-1;
}

void MHMatrix::AddColumns(MDataArray *matrix)
{
    if (fM.IsValid())
    {
        *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl;
        return;
    }

    if (TestBit(kIsLocked))
    {
        *fLog << warn << "Warning - matrix is locked. Can't add new columns... skipped." << endl;
        return;
    }

    if (fData)
        *fLog << warn << "Warning - columns already added... replacing." << endl;

    if (fData && TestBit(kIsOwner))
    {
        delete fData;
        ResetBit(kIsOwner);
    }

    fData = matrix;
}

// --------------------------------------------------------------------------
//
// Checks whether at least one column is available and PreProcesses all
// data chains.
//
Bool_t MHMatrix::SetupFill(const MParList *plist)
{
    if (!fData)
    {
        *fLog << err << "Error - No Columns initialized... aborting." << endl;
        return kFALSE;
    }

    return fData->PreProcess(plist);
}

// --------------------------------------------------------------------------
//
// If the matrix has not enough rows double the number of available rows.
//
void MHMatrix::AddRow()
{
    fNumRow++;

    if (fM.GetNrows() > fNumRow)
        return;

    if (!fM.IsValid())
    {
        fM.ResizeTo(1, fData->GetNumEntries());
        return;
    }

    TMatrix m(fM);

    fM.ResizeTo(fM.GetNrows()*2, fData->GetNumEntries());

    for (int x=0; x<m.GetNcols(); x++)
        for (int y=0; y<m.GetNrows(); y++)
            fM(y, x) = m(y, x);
}

// --------------------------------------------------------------------------
//
// Add the values correspoding to the columns to the new row
//
Bool_t MHMatrix::Fill(const MParContainer *par)
{
    AddRow();

    for (int col=0; col<fData->GetNumEntries(); col++)
        fM(fNumRow-1, col) = (*fData)(col);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Resize the matrix to a number of rows which corresponds to the number of
// rows which have really been filled with values.
//
Bool_t MHMatrix::Finalize()
{
    //
    // It's not a fatal error so we don't need to stop PostProcessing...
    //
    if (fData->GetNumEntries()==0 || fNumRow<1)
        return kTRUE;

    TMatrix m(fM);

    fM.ResizeTo(fNumRow, fData->GetNumEntries());

    for (int x=0; x<fM.GetNcols(); x++)
        for (int y=0; y<fM.GetNrows(); y++)
            fM(y, x) = m(y, x);

    return kTRUE;
}
/*
// --------------------------------------------------------------------------
//
// Draw clone of histogram. So that the object can be deleted
// and the histogram is still visible in the canvas.
// The cloned object are deleted together with the canvas if the canvas is
// destroyed. If you want to handle destroying the canvas you can get a
// pointer to it from this function
//
TObject *MHMatrix::DrawClone(Option_t *opt) const
{
    TCanvas &c = *MH::MakeDefCanvas(fHist);

    //
    // This is necessary to get the expected bahviour of DrawClone
    //
    gROOT->SetSelectedPad(NULL);

    fHist->DrawCopy(opt);

    TString str(opt);
    if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileX();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }
    if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileY();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }

    c.Modified();
    c.Update();

    return &c;
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the histogram into it.
// Be careful: The histogram belongs to this object and won't get deleted
// together with the canvas.
//
void MHMatrix::Draw(Option_t *opt)
{
    if (!gPad)
        MH::MakeDefCanvas(fHist);

    fHist->Draw(opt);

    TString str(opt);
    if (str.Contains("PROFX", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileX();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }
    if (str.Contains("PROFY", TString::kIgnoreCase) && fDimension==2)
    {
        TProfile *p = ((TH2*)fHist)->ProfileY();
        p->Draw("same");
        p->SetBit(kCanDelete);
    }

    gPad->Modified();
    gPad->Update();
}
*/

// --------------------------------------------------------------------------
//
// Prints the meaning of the columns and the contents of the matrix.
// Becareful, this can take a long time for matrices with many rows.
//
void MHMatrix::Print(Option_t *) const
{
    if (fData)
        fData->Print();
    else
        *fLog << all << "Sorry, no column information available." << endl;

    fM.Print();
}

const TMatrix *MHMatrix::InvertPosDef()
{
    TMatrix m(fM);

    const Int_t rows = m.GetNrows();
    const Int_t cols = m.GetNcols();

    for (int x=0; x<cols; x++)
    {
        Double_t avg = 0;
        for (int y=0; y<rows; y++)
            avg += fM(y, x);

        avg /= rows;

        for (int y=0; y<rows; y++)
            m(y, x) -= avg;
    }

    TMatrix *m2 = new TMatrix(m, TMatrix::kTransposeMult, m);

    Double_t det;
    m2->Invert(&det);
    if (det==0)
    {
        *fLog << err << "ERROR - MHMatrix::InvertPosDef failed (Matrix is sigular)." << endl;
        delete m2;
        return NULL;
    }

    // m2->Print();

    return m2;
}

// --------------------------------------------------------------------------
//
// 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();

    TArrayD dists(rows);

    //
    // Calculate:  v^T * M * v
    //
    for (int i=0; i<rows; i++)
    {
        TVector col(cols);
        col = TMatrixRow(fM, i);

        TVector d = evt;
        d -= col;

        TVector d2 = d;
        d2 *= m;

        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;
    }

    TArrayI idx(rows);
    TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);

    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)
{
    if (!fM2.IsValid())
    {
        const TMatrix &m = *InvertPosDef();
        fM2.ResizeTo(m);
        fM2 = m;
        fM2 *= fM.GetNrows()-1;
        delete &m;
    }

    return CalcDist(fM2, evt, num);
}

void MHMatrix::Reassign()
{
    TMatrix m = fM;
    fM.ResizeTo(1,1);
    fM.ResizeTo(m);
    fM = m;
}

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;
}

const TArrayI MHMatrix::GetIndexOfSortedColumn(Int_t ncol, Bool_t desc) const
{
    TMatrixColumn col(fM, ncol);

    const Int_t n = fM.GetNrows();

    TArrayF array(n);

    for (int i=0; i<n; i++)
        array[i] = col(i);

    TArrayI idx(n);
    TMath::Sort(n, array.GetArray(), idx.GetArray(), desc);

    return idx;
}

void MHMatrix::SortMatrixByColumn(Int_t ncol, Bool_t desc)
{
    TArrayI idx = GetIndexOfSortedColumn(ncol, desc);

    const Int_t n = fM.GetNrows();

    TMatrix m(n, fM.GetNcols());
    for (int i=0; i<n; i++)
    {
        TVector vold(n);
        vold = TMatrixRow(fM, idx[i]);

        TMatrixRow rownew(m, i);
        rownew = vold;
    }

    fM = m;
}

Bool_t MHMatrix::Fill(MParList *plist, MTask *read)
{
    //
    // Read data into Matrix
    //
    const Bool_t is = plist->IsOwner();
    plist->SetOwner(kFALSE);

    MTaskList tlist;
    plist->Replace(&tlist);

    MFillH fillh(this);

    tlist.AddToList(read);
    tlist.AddToList(&fillh);

    MEvtLoop evtloop;
    evtloop.SetParList(plist);

    if (!evtloop.Eventloop())
        return kFALSE;

    plist->Remove(&tlist);
    plist->SetOwner(is);

    return kTRUE;
}