#ifndef __MY_PCA_H_
#define __MY_PCA_H_

#include "DataSource.h"
#include <mkl.h>
#include <cassert>
#include <iostream>
#include <vector>

/////////////////////////////////////////////////////////////////////
//my PCA class
/////////////////////////////////////////////////////////////////////
template <typename T>
class CPCA
{
public:
	//const T *const *const m_pSamples;

	std::vector<float> m_Aver;
	std::vector<float> m_Covar;					//of size m_SrcDim*m_SrcDim
	std::vector<float> m_EigenValue;			//of size m_SrcDim
	std::vector<float> m_EigenVector;			//of size m_SrcDim*m_SrcDim

private:
	int m_SrcDim;
	int m_DstDim;

public:
	CPCA(const CDataSource<T> *const pDataSource, int dstDim=0)
		:m_SrcDim(pDataSource->m_Dim),m_DstDim(dstDim)
	{
		PreProcess(pDataSource);
	}

	CPCA(T *const *const src, int srcDim, int nSamples, int dstDim = 0)
		:m_SrcDim(srcDim),m_DstDim(dstDim)
	{
		CSimDataSource<T> dataSource(src, nSamples, srcDim);
		PreProcess(&dataSource);
	}

	~CPCA(){};

	void ShowData() const
	{
		std::cout.setf (std::ios::fixed | std::ios::left);
		std::cout.precision (3);

		std::cout<<"The Average Vector:\n";
		int i,j;
		for(i=0; i<m_SrcDim; i++){
			std::cout<<m_Aver[i]<<"\t";
		}

		std::cout<<std::endl<<"\nA:The CoVaraiance Matrix:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<m_Covar[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::cout<<std::endl<<"D:The EigenValue Vector:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<m_EigenValue[i]<<"\t";
		}

		std::cout<<std::endl<<"\nV:The EigenVector Matrix:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<m_EigenVector[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::vector<float> AV(m_SrcDim*m_SrcDim, 0);
		std::vector<float> VD(m_SrcDim*m_SrcDim, 0);
		std::vector<float> D(m_SrcDim*m_SrcDim, 0);
		for(i=0; i<m_SrcDim; i++){
			D[i*m_SrcDim+i] = m_EigenValue[i];
		}

		//dsymm(side, uplo, m, n, alpha, a, lda, 
		//		b, ldb, beta, c, ldc )
		const float alpha = 1.0L;
		const float beta = 0.0L;
		dsymm("L", "U", &m_SrcDim, &m_SrcDim, &alpha, &m_Covar.front(), &m_SrcDim, 
			&m_EigenVector.front(), &m_SrcDim, &beta, &AV.front(), &m_SrcDim);

		//dgemm(transa, transb, m, n, k, alpha, a, lda, 
		//		b, ldb, beta, c, ldc)
		dgemm("N", "N", &m_SrcDim, &m_SrcDim, &m_SrcDim, &alpha, &m_EigenVector.front(), &m_SrcDim, 
			&D.front(), &m_SrcDim, &beta, &VD.front(), &m_SrcDim);

		std::cout<<"\nAV:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<AV[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::cout<<"\nVD:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<VD[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

	}

	inline void ConvertPoint(const T *const src, T *const dst)	const
	{
		int j;
		for (j=0; j<m_DstDim; j++){
			dst[j] = 0;
			for(int k=0; k<m_SrcDim; k++){
				//dst[j] += (src[k] - m_Aver[k]) * m_EigenVector[k+(m_SrcDim-j-1)*m_SrcDim];
				dst[j] += (src[k] - m_Aver[k]) * m_EigenVector[k+j*m_SrcDim];
			}
		}
	}

	void ConvertData(T** src, T** dst, int nSample) const
	{
		for (int i=0; i<nSample; i++)
			ConvertPoint(src[i], dst[i]);
	}

private:
	void PreProcess(const CDataSource<T> *const pDataSource)
	{
		assert(pDataSource->m_Dim == m_SrcDim);
		if(pDataSource->m_Dim != m_SrcDim)
			return;

		assert(m_DstDim>0 && m_DstDim<=m_SrcDim);
		if(m_DstDim<=0 || m_DstDim>m_SrcDim)
			m_DstDim = m_SrcDim;

		const int nSamples = pDataSource->m_nSamples;
		std::cout<<"\n***********PCA PreProcessing with nSamples="<<nSamples<<" srcDim="<<m_SrcDim<<" dstDim="<<m_DstDim;
		clock_t tt=clock();
		//////////////////////////////////////////////////////////////////////////
		// calculating average
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\nCalculating Average...";
		clock_t t = clock();
		std::vector<long double> sums(m_SrcDim, 0);
		int i,j;
		std::vector<T> sample(m_SrcDim);
		for(i=0; i<nSamples; i++){
			pDataSource->getData(i, &sample.front());
			for(j=0; j<m_SrcDim; j++){
				sums[j] += sample[j];
			}
		}

		m_Aver.resize(m_SrcDim);
		for(i=0; i<m_SrcDim; i++){
			m_Aver[i] = sums[i]/nSamples;
		}

		std::vector<float> matt(nSamples*m_SrcDim);
		for(i=0; i<nSamples; i++){
			pDataSource->getData(i, &sample.front());
			for(j=0; j<m_SrcDim; j++){
				matt[i*m_SrcDim+j] = (sample[j] - m_Aver[j]);
			}
		}

		std::cout<<"\t\tfinished in "<<clock()-t<<" ms";
		t=clock();

		//////////////////////////////////////////////////////////////////////////
		// calculating covariance matrix
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\nCalculating Covariance...";

		m_Covar.resize(m_SrcDim*m_SrcDim, 0);
		const float alpha = 1.0L;
 		const float beta = 0.0L;
		//cblas_ssyrk(CblasRowMajor, CblasUpper, CblasTrans, m_SrcDim, nSamples,
		//				1.0L, &matSample.front(), m_SrcDim, 0.0L, &m_Covar.front(), m_SrcDim);

		//ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
 		ssyrk("U", "N", &m_SrcDim, &nSamples, &alpha,
 					&matt.front(), &m_SrcDim, &beta, &m_Covar.front(), &m_SrcDim);

		matt.swap(std::vector<float>());			// cleaning
		std::cout<<"\tfinished in "<<clock()-t<<"ms";
 		t=clock();

		//////////////////////////////////////////////////////////////////////////
		// calculating eigenvalues and eigenvectors
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\nSolving EigenProblem...";
		m_EigenValue.resize(m_SrcDim);
		m_EigenVector.resize(m_SrcDim*m_SrcDim);
		t=clock();

		int lwork, liwork, info;
		std::vector<float> work;
		std::vector<int> iwork;

		std::vector<int> isuppz(2*m_SrcDim);
		float vl, vu, abstol=0;
		int il=m_SrcDim-m_DstDim+1, iu=m_SrcDim;
		lwork = 26*m_SrcDim;
		liwork = 10*m_SrcDim;
		work.resize(lwork);
		iwork.resize(liwork);
		matt = m_Covar;
 		//ssyevr(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, 
		//				w, z, ldz, 
		//				isuppz, work, lwork, iwork, liwork, info)

		//ssyevr("V", "A", "U", &srcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_SrcDim, 
		ssyevr("V", "I", "U", &m_SrcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_DstDim, 
			&m_EigenValue.front(), &m_EigenVector.front(), &m_SrcDim, 
			&isuppz.front(), &work.front(), &lwork, &iwork.front(), &liwork, &info);
		
		assert(info == 0);
		if(info != 0){
			std::cout<<"Failded, Abort!"<<std::endl;
			return;
		}
		//std::cout<<"resutl: "<<info<<std::endl<<m<<" eigenvalues found"<<std::endl;

/*		//alternative 1: ssyevx
		lwork = 8*m_SrcDim;
		liwork = 5*m_SrcDim;
		work.resize(lwork);
		iwork.resize(liwork);
		std::vector<int> ifail(m_SrcDim);
		matt = m_Covar;
		//ssyevx(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, 
		//		w, z, ldz
		//		work, lwork, iwork, ifail, info)
		ssyevx("V", "I", "U", &srcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_DstDim, 
			&m_EigenValue.front(), &m_EigenVector.front(), &m_SrcDim, 
			&work.front(), &lwork, &iwork.front(), &ifail.front(), &info);
		
		// alternative 2:  calculating all Eigen Values & Eigen Vectors
		j=1,i=0;
		while(j<m_SrcDim){
			j <<= 1,i++;
		}

		lwork = 3*m_SrcDim*m_SrcDim + (5 + 2*i)*m_SrcDim+1;
		liwork = 5*m_SrcDim+3;
		work.resize(lwork);
		iwork.resize(liwork);
		//dsyevd ( job , uplo , n , a , lda , w , 
		//			work , lwork , iwork , liwork , info)
		matt = m_Covar;
		ssyevd("V", "U", &m_SrcDim, &matt.front(), &m_SrcDim, &m_EigenValue.front(), 
			&work.front(), &lwork, &iwork.front(), &liwork, &info);

		//reversing;
		for(i=0; i<m_SrcDim; i++)
			std::copy(matt.begin()+i*m_SrcDim, matt.begin()+(i+1)*m_SrcDim, m_EigenVector.begin()+(m_SrcDim-i-1)*m_SrcDim);
*/

		std::cout<<"\t\tfinished in "<<clock()-t<<"ms"<<std::endl;

		//////////////////////////////////////////////////////////////////////////
		// All finished
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"***********Finished in "<<clock()-tt<<" ms"<<std::endl;
	}

};


#endif