/** 
 * @file EigValues.cpp
 * @brief EigValues Class
 * @author Michael Koch
 * @date 08/19/2008
*/
#ifdef NOVISUAL
#include <nice_nonvis.h>
#else
#include <nice.h>
#endif
    
#include <core/iceconversion/convertice.h>
#include <image_nonvis.h>

#include <iostream>
     
#include "EigValues.h"
#include <time.h>

#define DEBUG_ARNOLDI
     
using namespace std;
// refactor-nice.pl: check this substitution
// old: using namespace ice;
using namespace NICE;
     
void EVArnoldi::getEigenvalues(const GenericMatrix &data,Vector &eigenvalues,Matrix &eigenvectors,uint k)
{
    if(data.rows()!=data.cols())
    {
	throw("EVArnoldi: matrix has to be quadratic");
    }
          
    if(k<=0)
    {
       throw("EVArnoldi: please use k>0.");
    }
   
    uint n=data.cols();
    ice::Matrix rmatrix (k,n,0);
    ice::Matrix qmatrix (k,n,0);//saves data =eigenvectors
    ice::Vector q (n);
    ice::Vector r (n);
    
    eigenvalues.resize(n);
   
    srand48(time(NULL));
    //random initialisation
    for(uint i=0;i<rmatrix.rows();i++)
	for(uint j=0;j<rmatrix.cols();j++)
	     rmatrix[i][j]=drand48();
   
    //reduceddim
    double delta=1.0;
    uint iteration=0;
    while (delta>mindelta &&iteration<maxiterations)
    {
	ice::Vector rold (rmatrix[k-1]);
	for(uint reduceddim=0;reduceddim<k;reduceddim++)
	{
	     //get Vector r_k of Matrix
	     q=rmatrix[reduceddim];
	     q.Normalize();
	     // q=r_k/||r_k||
	     qmatrix[reduceddim]=q;
	     
	     NICE::Vector tmp;
	     data.multiply( tmp, NICE::makeEVector(q)); //r_{k+1}=A*q
	     rmatrix[reduceddim] = NICE::makeIceVectorT(tmp);
	    
	     for(uint j=0;j<reduceddim;j++)
		  rmatrix[reduceddim]=rmatrix[reduceddim]-qmatrix[j]*(qmatrix[j]*rmatrix[reduceddim]);
	}
	//convergence stuff
	ice::Vector diff=rold-rmatrix[k-1];
	delta=diff.Length();
	iteration++;
#ifdef DEBUG_ARNOLDI
	fprintf (stderr, "EVArnoldi: [%d] delta=%f\n", iteration, delta );
#endif
    }
    eigenvectors = NICE::makeDoubleMatrix ( rmatrix );
    eigenvectors.transposeInplace();
  
    for(uint i=0;i<k;i++)
    {
	// refactor-nice.pl: check this substitution
	// old: Vector tmp;
	NICE::Vector tmp;
	eigenvectors.getColumnRef(i).normalizeL2();

	data.multiply(tmp,eigenvectors.getColumn(i));
	eigenvalues[i] = tmp.scalarProduct (eigenvectors.getColumn(i));
    }
}
     
#ifdef USE_BROKEN_LANCZOS
void EVLanczos::getEigenvalues(const GenericMatrix &data,Vector &eigenvalues,Matrix &eigenvectors,uint k)
{
    if (data.rows()!=data.cols())
	throw("Input matrix is not quadratic.");
     
    uint n=data.cols();
    uint m=data.rows();  
    if(k<=0 ) k=n;
    uint kint=k;
    NICE::Matrix wmatrix(kint+1,n,0);
    NICE::Matrix vmatrix(kint+1,n,0);//saves data =eigenvectors
    NICE::Matrix t(kint,kint,0);
    eigenvalues=Vector(k);
    NICE::Matrix eigenvectorsapprox(k,n,0);
    NICE::Vector w=Vector(n);
        
    // refactor-nice.pl: check this substitution
    // old: Vector alpha=Vector(kint+1);
    NICE::Vector alpha=Vector(kint+1);
    // refactor-nice.pl: check this substitution
    // old: Vector beta=Vector(kint+1);
    NICE::Vector beta=Vector(kint+1);
          
    //random initialisation ?
    vmatrix.Set(0.0);
                    
    if(n<=0) return;

    for (uint i=0;i<vmatrix[1].size();i++)
	// refactor-nice.pl: check this substitution
	// old: vmatrix[1][i]=RandomD();
	vmatrix(1, i)=RandomD();
               
    beta[0]=0;  
    //iteration
    double delta=1.0;
    uint iteration=0;
    while(delta>mindelta && iteration<maxiterations)
    {
	// refactor-nice.pl: check this substitution
	// old: Vector check=vmatrix[kint];
	NICE::Vector check=vmatrix[kint];
	for(uint i=1;i<kint;i++)
	{
	    //iteration step
	    NICE::Vector tmp;
	    data.multiply(tmp,vmatrix[i]); 
		     
	    wmatrix[i]=tmp-(beta[i-1]*vmatrix[i-1]);//w_i=A * v[i] - beta_i * v_(i-1)
	    alpha[i-1]=wmatrix[i]*vmatrix[i];
	    wmatrix[i]=wmatrix[i]-(alpha[i-1]*vmatrix[i]);
	    beta[i]=wmatrix[i].Length();
	    vmatrix[i+1]=1.0/beta[i]*wmatrix[i];
			       
	    t(i-1, i-1)=alpha[i];
	    if(i>1)
	    {
		t(i-1, i-2)=beta[i-1];
		t(i-2, i-1)=beta[i-1];
	    }
	  
	}
	// refactor-nice.pl: check this substitution
	// old: Vector diff=check-vmatrix[kint];
	NICE::Vector diff=check-vmatrix[kint];
	delta=diff.Length();
	iteration++;
    }
    cout<<"Lanczos-Iterations:"<<iteration<<"/"<<maxiterations<<"relative Error:"<<delta<<endl;
    // refactor-nice.pl: check this substitution
    // old: Vector ice_eigval;
    NICE::Vector ice_eigval;
    // refactor-nice.pl: check this substitution
    // old: Matrix ice_eigvect;
    NICE::Matrix ice_eigvect;
    Eigenvalue ( t, ice_eigval, ice_eigvect );

    cerr << ice_eigval << endl;
	   

    for (uint l=0;l<k;l++)
	eigenvectorsapprox[l]=vmatrix[l+1];

    eigenvectors=eigenvectorsapprox;
    for (uint i=0;i<(uint)eigenvectors.rows();i++)
    {
	// refactor-nice.pl: check this substitution
	// old: Vector tmp;
	NICE::Vector tmp;
	if ( eigenvectors[i].Length() < 1e-12 )
	{
	    cerr << "Lanczos: numerical instability" << endl;
	} else {
	    eigenvectors[i].Normalize();
	}
	data.multiply(tmp,eigenvectors[i]);
	eigenvalues[i]=tmp * eigenvectors[i];
    }
}
#endif