CMDS.h

Go to the documentation of this file.

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

#ifndef REDUKT_CMDS_H
#define REDUKT_CMDS_H

#define HAVE_IETL   

#include "Eigenvalue.h"
#include "Cholesky.h"
#include "ublasTools.h"

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/operation.hpp>

#ifdef HAVE_IETL
#include "Lanczos.h"
#endif

#include <assert.h>


template<class ValueType>
class CMDS
{
public:
        typedef ublas::matrix<ValueType> Matrix;
        typedef ublas::vector<ValueType> Vector;

        enum DecompositionMethod 
             { 
                         MethodQR        
                        ,MethodCholesky  
//#ifdef HAVE_IETL  
                        ,MethodLanczos   
//#endif                         
                 };

        CMDS( Matrix& D, int method = (int)MethodQR ): m_D(D), m_method(method) 
        {
                using namespace ublas;  // for column
                using namespace ublasTools;

                assert( D.size1() == D.size2() );

                // square components of D and multiply by -0.5
                for( unsigned int i=0; i < D.size1(); ++i )
                        for( unsigned int j=0; j < D.size2(); ++j )
                                D(i,j) = -0.5 * D(i,j)*D(i,j);

                // double-center D to get Torgersons scalar product matrix
                double_center( D );

                // force symmetric matrix
                if( !Eigenvalue<double>::isSymmetric( D ) )
                {
                        //cout << "Matrix not fully symmetric, assuming it should be!" << endl;
                        //cout << "Replacing upper with lower triangular matrix..." << std::flush;
                        for( unsigned int i=0; i < D.size1(); ++i )
                                for( unsigned int j=0; j < i; ++j )
                                        D(j,i) = D(i,j);
                        //cout << "done" << endl;
                }

                // compute decomposition D=XX'
                if( method == MethodCholesky )
                {               
                        Cholesky<double> cholesky( D );
                        m_D = ublas::trans( cholesky.getL() );
                
                        // no eigenvalues here ;-)              
                }
#ifdef HAVE_IETL
                if( method == MethodLanczos )
                {
                        // computation is done in project() function                    
                }
#endif
                else // MethodCholesky
                {
                        //cout << "Computing Eigendecomposition..." << std::flush;
                        Eigenvalue<double> evd( D );
                        //cout << "done" << endl;

                        // get eigenvalues
                        m_realEV = evd.getRealEigenvalues();

                        // overwrite D matrix with eigenvector matrix
                        m_D = evd.getV();
                }
                        

                // leave projection to user
        }

        const Vector& getEigenvalues() const { return m_realEV; }
        
        void project( int dim=2 )
        {       
                assert( (dim >= 1) && (dim <= (int)m_D.size1()) );
                
                if( m_method == MethodQR )
                {               
                        Matrix S( dim, dim );
                        for( int i=0; i < dim; ++i )
                                for( int j=0; j < dim; ++j )
                                        S(i,j) = (i==j ? sqrt(m_realEV[i]) : 0.0);
                
                        m_D.resize( m_D.size1(), dim );
                                
                        m_D = prod( m_D, S );
                }
#ifdef HAVE_IETL
                if( m_method == MethodLanczos )
                {
                        Lanczos<Matrix,Vector> lanczos( m_D, dim );
                        
                        m_D.resize( m_D.size1(), dim );
                        
                        m_realEV.resize( dim );
                        
                        // store eigenvectors in columns of D
                        for( int i=0; i < dim; ++i      )
                        {
                                m_realEV[i] = lanczos.getEigenvalue( i );
                                ublas::column( m_D, i ) = lanczos.getEigenvector( i );
                        }
                        
                        Matrix S( dim, dim );
                        for( int i=0; i < dim; ++i )
                                for( int j=0; j < dim; ++j )
                                        S(i,j) = (i==j ? sqrt(m_realEV[i]) : 0.0);
                
                        m_D = prod( m_D, S );                           
                }
#endif
                else // MethodCholesky
                {
                        m_D.resize( dim, m_D.size2() ); // points in columns
                        m_D = ublas::trans( m_D );      // transpose for output
                }
        }

private:
        Matrix& m_D;
        Vector  m_realEV;
        int m_method;
};

#endif // REDUKT_CMDS_H