FeatureCovar.h

Go to the documentation of this file.

/******************************************************************************
 *
 *  Copyright (C) 2008 Boris Iven, Max Hermann
 *
 *  This file is part of the "CEDRIC Event Display" application.
 *
 *****************************************************************************/

#ifndef FEATURECOVAR_H
#define FEATURECOVAR_H

#include <iostream>
#include <vector>
#include <math.h>

template<class F,int DIM_FEAT>
class FeatureCovar
{
public:
        FeatureCovar( const std::vector<F>& values_ ); 
        
        void compute();
        void print();

        unsigned int size() const { return DIM_FEAT; };
        
        double getUnitVarCovMat   (int r, int c) const;
        double getInvUnitVarCovMat(int r, int c) const;
        double getRawCovMat       (int r, int c) const;
        double getInvRawCovMat    (int r, int c) const;
        double getVarScaledCovMat (int r, int c) const;
        double getMean(int d) const;

private:
        double rawCovMat[DIM_FEAT][DIM_FEAT];
        double invRawCovMat[DIM_FEAT][DIM_FEAT];
        double unitVarCovMat[DIM_FEAT][DIM_FEAT];
        double varScaledCovMat[DIM_FEAT][DIM_FEAT];
        double invUnitVarCovMat[DIM_FEAT][DIM_FEAT];
        double means[DIM_FEAT];
        
        const std::vector<F>& values;
        std::vector<F> zeroMeanValues;
        std::vector<F> varScaledValues;
        std::vector<F> unitVarianceValues;
        
        int invert (double a[][2*DIM_FEAT], double ainv[][2*DIM_FEAT], int n);
};



//=====================================================================================================================

//              TEMPLATE IMPLEMENTATION

//=====================================================================================================================


template<class F,int DIM_FEAT>
FeatureCovar<F,DIM_FEAT>::FeatureCovar( const std::vector<F>& values_ )
: values(values_)
{       
        for(int r = 0; r < DIM_FEAT; r++)
                for(int c = 0; c < DIM_FEAT; c++)
                {       
                        means[r] = 0;
                        rawCovMat[r][c] = 0;
                        unitVarCovMat[r][c] = 0;
                        invUnitVarCovMat[r][c] = 0;
                        varScaledCovMat[r][c] = 0;
                }                       
}

template<class F,int DIM_FEAT>
void FeatureCovar<F,DIM_FEAT>::compute()
{
        using namespace std;
        
        //for matrix inversion
        double tmp[2*DIM_FEAT][2*DIM_FEAT];
        double res[2*DIM_FEAT][2*DIM_FEAT];

        cout << "Computing covariance matrix... " << endl;
        
        //vector<double> means(6,0);
        
        int num = values.size();
        
        //compute means
        
        for (int v = 0; v < num; v++)
        {
                for(int d = 0; d < DIM_FEAT; d++)
                {
                        means[d] += values[v][d] / num;
                }
        }
        
        cout << "Means: " << endl;
        
        for (int d = 0; d < DIM_FEAT; d++)
                cout << d << ": " << means[d] << endl;
        cout << endl;
        
        //compute "raw" covariance matrix
        
        for (int r = 0; r < DIM_FEAT; r++)
        {
                for (int c = 0; c < DIM_FEAT; c++)
                {
                        for (int v = 0; v < num; v++)
                        rawCovMat[r][c] += (values[v][r] - means[r]) * (values[v][c] - means[c]) / num;
                }
        }
        
        //compute inverse raw covariance matrix
        
        for (int r = 0; r < DIM_FEAT; r++)
                for (int c = 0; c < DIM_FEAT; c++)
                        tmp[r][c] = rawCovMat[r][c];
        
        invert(tmp, res, DIM_FEAT);
        
        
        
        for (int r = 0; r < DIM_FEAT; r++)
                for (int c = 0; c < DIM_FEAT; c++)
                        invRawCovMat[r][c] = res[r][c];
        
        //compute zero mean values
        zeroMeanValues = values;
        
        for (int i = 0; i < num; i++)
        {
                for (int d = 0; d < DIM_FEAT; d++)
                        zeroMeanValues[i][d] -= means[d];
        }
        
        //compute variance scaled values
        
        varScaledValues = zeroMeanValues;
        
        for (int i = 0; i < num; i++)
        {
                for (int d = 0; d < DIM_FEAT; d++)
                        varScaledValues[i][d] /= rawCovMat[d][d];
        }

        
        //compute unit variance values

        unitVarianceValues = zeroMeanValues;
        
        for (int i = 0; i < num; i++)
        {
                for (int d = 0; d < DIM_FEAT; d++)
                        unitVarianceValues[i][d] /= sqrt(rawCovMat[d][d]);
        }

        //compute unit variance covariance matrix
        
        for (int r = 0; r < DIM_FEAT; r++)
        {
                for (int c = 0; c < DIM_FEAT; c++)
                {
                        for (int v = 0; v < num; v++)
                        unitVarCovMat[r][c] += ((unitVarianceValues[v][r]) * (unitVarianceValues[v][c]) / num);
                        
                }
        }
        
        //compute mysterious variance scaled covariance matrix
        
        for (int r = 0; r < DIM_FEAT; r++)
        {
                for (int c = 0; c < DIM_FEAT; c++)
                {
                        for (int v = 0; v < num; v++)
                        varScaledCovMat[r][c] += ((varScaledValues[v][r]) * (varScaledValues[v][c]) / num);
                        
                }
        }

        
        //compute inverse covariance matrix

        
        for (int r = 0; r < DIM_FEAT; r++)
                for (int c = 0; c < DIM_FEAT; c++)
                        tmp[r][c] = unitVarCovMat[r][c];
        
        invert(tmp, res, DIM_FEAT);
        
        
        
        for (int r = 0; r < DIM_FEAT; r++)
                for (int c = 0; c < DIM_FEAT; c++)
                        invUnitVarCovMat[r][c] = res[r][c];
}

template<class F,int DIM_FEAT>
void FeatureCovar<F,DIM_FEAT>::print()
{
        using namespace std;
        
        cout << "ORIGINAL (RAW) COVARIANCE MATRIX:" << endl << endl;;
        for(int r = 0; r < DIM_FEAT; r++)
        {
                for(int c = 0; c < DIM_FEAT; c++)
                {
                        cout << rawCovMat[r][c] << "\t";
                }
                cout << endl;
        }
        cout << endl;

        cout << "INVERSE (RAW) COVARIANCE MATRIX:" << endl << endl;;
        for(int r = 0; r < DIM_FEAT; r++)
        {
                for(int c = 0; c < DIM_FEAT; c++)
                {
                        cout << invRawCovMat[r][c] << "\t";
                }
                cout << endl;
        }
        cout << endl;

        
        cout << "VAR^-1-SCALED COVARIANCE MATRIX(EXPERIMENTAL):" << endl << endl;;
        for(int r = 0; r < DIM_FEAT; r++)
        {
                for(int c = 0; c < DIM_FEAT; c++)
                {
                        cout << varScaledCovMat[r][c] << "\t";
                }
                cout << endl;
        }
        cout << endl;

        
        cout << "UNIT VARIANCE COVARIANCE MATRIX:" << endl << endl;;
        for(int r = 0; r < DIM_FEAT; r++)
        {
                for(int c = 0; c < DIM_FEAT; c++)
                {
                        cout << unitVarCovMat[r][c] << "\t";
                }
                cout << endl;
        }
        cout << endl;
        
        cout << "INVERSE UNIT VARIANCE COVARIANCE MATRIX:" << endl << endl;;
        for(int r = 0; r < DIM_FEAT; r++)
        {
                for(int c = 0; c < DIM_FEAT; c++)
                {
                        cout << invUnitVarCovMat[r][c] << "\t";
                }
                cout << endl;
        }
        cout << endl;

}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getUnitVarCovMat(int r, int c) const
{
        return unitVarCovMat[r][c];
}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getInvUnitVarCovMat(int r, int c) const
{
        return invUnitVarCovMat[r][c];
}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getRawCovMat(int r, int c) const
{
        return rawCovMat[r][c];
}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getInvRawCovMat(int r, int c) const
{
        return invRawCovMat[r][c];
}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getVarScaledCovMat(int r, int c) const
{
        return varScaledCovMat[r][c];
}

template<class F,int DIM_FEAT>
double FeatureCovar<F,DIM_FEAT>::getMean(int d) const
{
        return means[d];
}

template<class F,int DIM_FEAT>
int FeatureCovar<F,DIM_FEAT>::invert (double a[][2*DIM_FEAT], double ainv[][2*DIM_FEAT], int n) 
{
        using namespace std;
        
  int   i, j;                    // Zeile, Spalte
  int   s;                       // Elimininationsschritt
  int   pzeile;                  // Pivotzeile
  int   fehler = 0;              // Fehlerflag
  double f;                      // Multiplikationsfaktor
  const double Epsilon = 0.01;   // Genauigkeit
  double Maximum;                // Zeilenpivotisierung
  //extern FILE *fout;
  int pivot = 1;

  // ergänze die Matrix a um eine Einheitsmatrix (rechts anhängen)
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
    {
      a[i][n+j] = 0.0;
      if (i == j)
        a[i][n+j] = 1.0;
    }
  }
  
  // die einzelnen Eliminationsschritte
  s = 0;
  do {
          // Pivotisierung vermeidet unnötigen Abbruch bei einer Null in der Diagnonalen und
                // erhöht die Rechengenauigkeit
    Maximum = fabs(a[s][s]);
    if (pivot)
    {
       pzeile = s ; 
       for (i = s+1; i < n; i++)
         if (fabs(a[i][s]) > Maximum) {
           Maximum = fabs(a[i][s]) ;
           pzeile = i;
         }
    }
    fehler = (Maximum < Epsilon);

    if (fehler) break;           // nicht lösbar 

    if (pivot)
    {
      if (pzeile != s)  // falls erforderlich, Zeilen tauschen
      { double h;
        for (j = s ; j < 2*n; j++) {
          h = a[s][j];
          a[s][j] = a[pzeile][j];
          a[pzeile][j]= h;
        }
      }
    }

    // Eliminationszeile durch Pivot-Koeffizienten f = a[s][s] dividieren
    f = a[s][s];
    for (j = s; j < 2*n; j++)
      a[s][j] = a[s][j] / f;

    // Elimination --> Nullen in Spalte s oberhalb und unterhalb der Diagonalen
    // durch Addition der mit f multiplizierten Zeile s zur jeweiligen Zeile i
    for (i = 0; i < n; i++ ) {
      if (i != s) 
      {
        f = -a[i][s];                 // Multiplikationsfaktor
        for (j = s; j < 2*n ; j++)    // die einzelnen Spalten
          a[i][j] += f*a[s][j];       // Addition der Zeilen i, s
      }
    }
    s++;
  } while ( s < n ) ;

  if (fehler) 
  {
    cout <<"Inverse: Matrix ist singulär\n" << endl;
    return 0; 
  }
  // Die angehängte Einheitsmatrix Matrix hat sich jetzt in die inverse Matrix umgewandelt
  // Umkopieren auf die Zielmatrix
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      ainv[i][j] = a[i][n+j];
    }
  }
  return 1;  
}

#endif // FEATURECOVAR_H