NICE_VisLearning/math/progs/testEigenvalues.cpp

/** 
* @file testEigenvalues.cpp
* @brief test eigenvalues
* @author Erik Rodner
* @date 05/21/2008

*/


#include "core/vector/VectorT.h"
#include "core/vector/MatrixT.h"
#include "core/image/ImageT.h"
#include "core/imagedisplay/ImageDisplay.h"

#ifdef NICE_USELIB_ICE

#include <image_nonvis.h>
#include <core/iceconversion/convertice.h>
#include <core/vector/Distance.h>

#include "vislearning/math/algebra/EigValues.h"
#ifdef NICE_USELIB_TRLAN
#include "vislearning/math/algebra_trlan/EigValuesTRLAN.h"
#endif
#include "core/basics/Config.h"

using namespace OBJREC;

using namespace NICE;
using namespace std;


int main (int argc, char **argv)
{   
    std::set_terminate(__gnu_cxx::__verbose_terminate_handler);

    Config conf ( argc, argv );
    uint rows = conf.gI("main", "n", 2);
    uint cols = rows;
    uint k = conf.gI("main", "k", rows);
    uint maxiterations = conf.gI("main", "maxiterations", 100);
    double mindelta = conf.gD("main", "mindelta", 1e-4 );
    double sparse_prob = conf.gD("main", "sparseprob", 0.0 );
    bool trlan = conf.gB("main", "trlan", false);
    int trlan_magnitude = conf.gI("main", "trlan_magnitude", 1);
    bool init_random = conf.gB("main", "init_random", true );

    // refactor-nice.pl: check this substitution
    // old: Matrix T (rows, cols);
    NICE::Matrix T (rows, cols);
    T.set(0.0);
    if ( init_random )
	srand48(time(NULL));

    // generate random symmetric matrix
    for ( uint i = 0 ; i < rows ; i++ )
	for ( uint j = i ; j < cols ; j++ )
	{
	    if ( sparse_prob != 0.0 )
		if ( drand48() < sparse_prob )
		    continue;
	    T(i, j) = drand48();
	    T(j, i) = T(i, j);
	}

    // refactor-nice.pl: check this substitution
    // old: Vector ice_eigvalues;
    ice::Vector ice_eigvalues;
    // refactor-nice.pl: check this substitution
    // old: Matrix ice_eigvect;
    ice::Matrix ice_eigvect;
    Eigenvalue ( NICE::makeICEMatrix(T), ice_eigvalues, ice_eigvect );

    EigValues *eig; 
    
    if ( ! trlan )
    {
		fprintf (stderr, "Eigenvalue method: Arnoldi iteration\n");
		eig = new EVArnoldi( maxiterations, mindelta );
    } else {
#ifdef NICE_USELIB_TRLAN
		fprintf (stderr, "Eigenvalue method: Lanczos (TRLAN) magn=%d\n", trlan_magnitude);
		eig = new EigValuesTRLAN(trlan_magnitude);
#else 
		fprintf (stderr, "Eigenvalue method: Lanczos (TRLAN) not supported !\n");
		exit(-1);
#endif
    }
    
    NICE::Vector eigvalues_dense;
    NICE::Matrix eigvect_dense;
    NICE::Vector eigvalues_sparse;
    NICE::Matrix eigvect_sparse;

    GMSlowICE Tg (T);
    double start_dense = ice::TimeD();
    eig->getEigenvalues ( Tg, eigvalues_dense, eigvect_dense, k );
    double time_dense = ice::TimeD() - start_dense;

    GMSparse Ts (T);
    double start_sparse = ice::TimeD();
    eig->getEigenvalues ( Ts, eigvalues_sparse, eigvect_sparse, k );
    double time_sparse = ice::TimeD() - start_sparse;

    fprintf (stderr, "time: sparse %fs, dense %fs\n", time_sparse, time_dense );

    // test property
    NICE::EuclidianDistance<double> eucliddist;
    for ( uint i = 0 ; i < k ; i++ )
    {
		NICE::Vector v_dense = eigvect_dense.getColumn(i);
		double lambda_dense = eigvalues_dense[i];
		NICE::Vector Tv_dense;
		Tv_dense.multiply ( T, v_dense );
		NICE::Vector lv_dense = v_dense;
		lv_dense *= lambda_dense;
		double err_dense = eucliddist(Tv_dense, lv_dense);
		
		NICE::Vector v_sparse = eigvect_sparse.getColumn(i);
		double lambda_sparse = eigvalues_sparse[i];
			NICE::Vector Tv_sparse;
		Tv_sparse.multiply ( T, v_sparse );
		NICE::Vector lv_sparse = v_sparse;
		lv_sparse *= lambda_sparse;
		double err_sparse = eucliddist(Tv_sparse, lv_sparse);

		ice::Matrix ice_T = makeICEMatrix(T);
		cerr << T << endl;
		cerr << ice_T << endl;
		double err_ice = (ice_T*(!ice_eigvect)[i] - ice_eigvalues[i]*(!ice_eigvect)[i]).Length();
		cerr << "eigenvalue " << i << endl;
		cerr << "lambda (dense) = " << lambda_dense << endl;
		cerr << "lambda (sparse) = " << lambda_sparse << endl;
		cerr << "lambda (ICE) = " << ice_eigvalues[i] << endl;
		cerr << "||Av - lambda v|| (dense) = " << err_dense << endl;
		cerr << "||Av - lambda v|| (sparse) = " << err_sparse << endl;
		cerr << "||Av - lambda v|| (ICE) = " << err_ice << endl;
    }

    // test covariance part

    cerr << "------- Covariance Test -------" << endl;

    int n = 2*rows; // number of data elements
    NICE::Matrix data(rows, n);
    ice::Statistics stat (rows);
    for (int i = 0 ; i < n ; i++ )
    {
		for ( int j = 0 ; j < (int)rows ; j++ )
			data(j, i) = drand48() + 100.0;

		Put ( stat, NICE::makeIceVectorT( data.getColumn(i) ) );
    }
	cerr << "Mean: " << Mean(stat) << endl;
    NICE::Matrix cov = NICE::makeDoubleMatrix ( Covariance(stat) );

    GMSlowICE cov_slow ( cov );
    GMCovariance cov_cov ( &data );
    NICE::Vector x1t_slow;
    cov_slow.multiply ( x1t_slow, data.getColumn(0) );
    cerr << x1t_slow << endl;

    NICE::Vector x1t_cov;
    cov_cov.multiply ( x1t_cov, data.getColumn(0) );
    cerr << x1t_cov << endl;
    
    cerr << "Difference (simple vector multiplication): " << eucliddist( x1t_slow, x1t_cov ) << endl;

    double start_std = ice::TimeD();
    eig->getEigenvalues ( cov_slow, eigvalues_dense, eigvect_dense, k );
    double time_std = ice::TimeD() - start_std;

    NICE::Vector eigvalues_cov;
    NICE::Matrix eigvect_cov;
    
    double start_cov = ice::TimeD();
    eig->getEigenvalues ( cov_cov, eigvalues_cov, eigvect_cov, k );
    double time_cov = ice::TimeD() - start_cov;

    cerr << "Eigenvalues (std) = " << eigvalues_dense << endl;
    cerr << "Eigenvalues (covariance) = " << eigvalues_cov << endl;

    fprintf (stderr, "time: cov %fs, std %fs\n", time_std, time_cov );
    
    
    return 0;
}

#else
int main (int argc, char **argv)
{   
    throw("no ice installed\n");
    return 0;
}
#endif