/** 
* @file calcKernelMatrixProperties.cpp
* @brief calc some kernel matrix properties
* @author Erik Rodner
* @date 01/15/2010

*/
#include "core/basics/Config.h"
#include "core/vector/Eigen.h"

using namespace std;
using namespace NICE;


/** 
    
    calc some kernel matrix properties 
    
*/
int main (int argc, char **argv)
{   
    std::set_terminate(__gnu_cxx::__verbose_terminate_handler);

    Config conf ( argc, argv );
    
	string in = conf.gS("main", "in");

	cout << "reading kernel matrix ..." << endl;
	ifstream ifs ( in.c_str(), ios::in );

	if ( !ifs.good() )
		fthrow(IOException, "Unable to read " << in << "." );

	Matrix kernelMatrixOrig;
	ifs >> kernelMatrixOrig;

	ifs.close();

	int norig = kernelMatrixOrig.rows();
	int morig = kernelMatrixOrig.cols();

	if ( norig != morig )
		fthrow(Exception, "Kernel matrix has to be quadratic.");

	int size = conf.gI("main", "size", 100 );
	if ( size > norig ) size = norig;
	
	Matrix kernelMatrix ( size, size );
	for ( int i = 0 ; i < kernelMatrix.rows(); i++ )
		for ( int j = 0 ; j < kernelMatrix.cols(); j++ )
			kernelMatrix(i,j) = kernelMatrixOrig(i,j);
			
	int n = kernelMatrix.rows();
	int m = kernelMatrix.cols();

	cout << "size: " << n << " x " << m << " (orig: " << norig << " x " << norig << " )" << endl;
	cout << "performing eigenvalue decomposition ..." << endl;
    Vector eigvals (n);
	eigenvalues ( kernelMatrix, &eigvals );

	double min_lambda = eigvals.Min();
	double max_lambda = eigvals.Max();

	cout << "smallest eigenvalue: " << min_lambda << endl;
	cout << "largest eigenvalue: " << max_lambda << endl;
   
    double sigma = conf.gD("main", "sigma", 0.01);
    double upper_bound = 0.5 * sigma * sigma / ( pow(min_lambda,3) );
	double lower_bound = 0.5 * sigma * sigma / ( pow(min_lambda,3) + sigma*sigma );
    cout << "inverse approximation bounds (cf. GP Multi-Task Notes): " << lower_bound << " <= |R|_2 <= " <<  upper_bound << endl;
    
    
    return 0;
}