//
// C++ Implementation: TestEigenValue
//
// Description:
//
//
// Author: Michael Koch <Koch.Michael@uni-jena.de>, (C) 2009
//
// Copyright: See COPYING file that comes with this distribution
//
/**
 * @file TestEigenValue.cpp
 * @brief TestEigenValue
 * @author Michael Koch
 * @date Di Aug 4 2009
 */

#include "TestEigenValue.h"
#include <string>

#include "core/basics/cppunitex.h"
#include "core/basics/numerictools.h"
#include "core/vector/Distance.h"

#include "vislearning/math/algebra/EigValues.h"
#include "vislearning/math/algebra/AlgebraTools.h"
#include "vislearning/math/algebra/GenericMatrix.h"

#include "vislearning/nice_nonvis.h"

using namespace std;
using namespace NICE;
using namespace OBJREC;

CPPUNIT_TEST_SUITE_REGISTRATION(TestEigenValue);

void TestEigenValue::setUp()
{
}

void TestEigenValue::tearDown()
{
}

void TestEigenValue::TestEigenValueComputation()
{
    uint rows = 3;
    uint cols = rows;
    uint k = rows;
    uint maxiterations =  200;
    double mindelta = 1e-8;
    double sparse_prob = 0.0;
    int trlan_magnitude = 1;
    bool 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;


    EigValues *eig;
    for (int trlan = 0;trlan <= 1;trlan++) //this is creepy but funny
    {
        if (trlan) //this is creepy but saves lot of code
        {
#ifdef NICE_USELIB_TRLAN
            eig = new EigValuesTRLAN(trlan_magnitude);
#else
            break;
#endif
        }
        else
        {
            eig = new EVArnoldi(maxiterations, mindelta);
        }

        NICE::Vector eigvalues_dense;
        NICE::Matrix eigvect_dense;
        NICE::Vector eigvalues_sparse;
        NICE::Matrix eigvect_sparse;

        GMSlowICE Tg(T);
        eig->getEigenvalues(Tg, eigvalues_dense, eigvect_dense, k);

        GMSparse Ts(T);
        eig->getEigenvalues(Ts, eigvalues_sparse, eigvect_sparse, k);

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

//             cerr << "eigenvalue " << i << endl;
//             cerr << "lambda (dense) = " << lambda_dense << endl;
//             cerr << "lambda (sparse) = " << lambda_sparse << endl;
//             cerr << "||Av - lambda v|| (dense) = " << err_dense << endl;
//             cerr << "||Av - lambda v|| (sparse) = " << err_sparse << endl;
            CPPUNIT_ASSERT_DOUBLES_EQUAL_NOT_NAN(0.0,err_dense,1e-4);
            CPPUNIT_ASSERT_DOUBLES_EQUAL_NOT_NAN(0.0,err_sparse,1e-4);
        }
    }
}