NICE_VisLearning/math/ftransform/PCA.cpp

/**
 * @file PCA.cpp
 * @brief Computation of a PCA by Eigenvalue Decomposition, Karhunen Loewe Transform
 * @author Michael Koch
 * @date 5/27/2008

 */
/** \example testPCA.cpp
 * This is an example of how to use the Principal Component Analysis (PCA) class.
 */

// #define DEBUG
#undef DEBUG

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

#include <iostream>
#include <math.h>
#include <vector>
#include <map>

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

#ifdef NICE_USELIB_TRLAN
#include "vislearning/math/algebra_trlan/EigValuesTRLAN.h"
#endif
// #include "vislearning/fourier/FourierLibrary.h"
#include "PCA.h"
#include "vislearning/baselib/Gnuplot.h"

#ifdef NICE_USELIB_ICE
#include "core/image/ImageT.h"
#include "core/vector/VectorT.h"
#include "core/vector/MatrixT.h"

#include "core/iceconversion/convertice.h"
#include "core/iceconversion/image_convertice.h"
using namespace ice;
#endif

using namespace OBJREC;

using namespace std;
using namespace NICE;

PCA::PCA(uint dim, uint maxiteration, double mindelta)
{
	init(dim, maxiteration, mindelta);
}

PCA::PCA(void)
{
	init();
}

PCA::~PCA()
{
}

void PCA::init(uint dim, uint maxiteration, double mindelta)
{
	this->targetDimension = dim;
	this->maxiteration = maxiteration;
	this->mindelta = mindelta;
}

void PCA::restore(istream & is, int format)
{

	is >> basis;
	is >> normalization;
	is >> mean;
	is >> targetDimension;
}
void PCA::store(ostream & os, int format) const
{
	os << basis << normalization << mean << targetDimension;
}
void PCA::clear()
{

}

void PCA::calculateBasis(const NICE::Matrix &features,
		const uint targetDimension, const uint mode)
{
	calculateBasis(features, targetDimension, mode, false);
}
void PCA::calculateBasis(const NICE::Matrix &features,
		const uint targetDimension, const uint mode, const bool adaptive,
		const double targetRatio)
{
	this->targetDimension = targetDimension;
	uint srcFeatureSize = features.cols();
	uint mindimension = std::min(this->targetDimension, srcFeatureSize);

#ifdef DEBUG
	cout << "calculateBasis ... stay a while and listen ..." << endl;
#endif

	NICE::Matrix eigenvectors;

	NICE::Vector eigenvalue;
#ifdef DEBUG
	cout << "EigenValue Decomposition Working..." << endl;
#endif
	//switch mode: svd, lanczos, broken lanczos, arnoldi
	//switch eigenvalue decomposition
	switch (mode)
	{
	case 0:
	{
#ifndef NICE_USELIB_ICE
		fthrow(Exception, "PCA: ICE needed for mode 0");
#else
#ifdef DEBUG
		fprintf(stderr, "PCA: warning covariance computation in progress, memory and time consuming ...!\n");
#endif
		ice::Matrix sigma;
		ice::Matrix eigenvectors_ice;
		ice::Vector eigenvalue_ice;
		ice::Matrix features_ice(NICE::makeICEMatrix(features));

		Statistics st(srcFeatureSize);
#ifdef DEBUG
		fprintf(stderr, "Transform Features ... %d %d \n", features_ice.rows(), features_ice.cols());
#endif
		Put(st, features_ice);

		ice::Vector mean_ice = Mean(st);
		mean = makeEVector(mean_ice);
#ifdef DEBUG
		cout << "MEAN: " << mean_ice << endl;
#endif
		sigma = Covariance(st);
#ifdef DEBUG
		fprintf(stderr, "Covariance Ready \n");
		fprintf(stderr, "Eigenvalue Computation ... \n");
#endif
		Eigenvalue(sigma, eigenvalue_ice, eigenvectors_ice);
		NICE::Matrix eigenvectors_tmp;
		eigenvectors_tmp = makeDoubleMatrix(eigenvectors_ice); //bad error
		eigenvectors.resize(eigenvectors_tmp.rows(), eigenvectors_tmp.cols());
		eigenvectors = eigenvectors_tmp;
		eigenvalue = makeEVector(eigenvalue_ice);

#endif
		break;
	}
	case 1:
	{

#ifdef DEBUG
		fprintf(stderr, "Calc mean \n");
#endif
		AlgebraTools::calculateMean(features, mean);
#ifdef DEBUG
		fprintf(stderr, "fastest available method \n");
		//                cerr<<"mean:"<<mean<<endl;
#endif
#ifdef NICE_USELIB_TRLAN
		EigValues *eig = new EigValuesTRLAN();//fast lanczos TRLAN
#else
		EigValues *eig = new EVArnoldi();//Arnoldi for (srcFeatureSize<n)
#endif
		NICE::Matrix features_transpose = features.transpose();
		GMCovariance C(&features_transpose);
#ifdef DEBUG
		fprintf(stderr, "Perform PCA\n");
		fprintf(stderr, "Matrix size %d %d %d %d \n", features.rows(), features.cols(), C.rows(), C.cols());
#endif
		if (adaptive)
		{
			eig->getEigenvalues(C, eigenvalue, eigenvectors, srcFeatureSize);
		}
		else
		{
			eig->getEigenvalues(C, eigenvalue, eigenvectors, mindimension);
		}

#ifdef DEBUG
		fprintf(stderr, "PCA: Computation ready \n");
		//                cerr<<eigenvectors<<endl;
#endif

#ifdef DEBUG
		fprintf(stderr, "PCA: Transpose ready \n");
		//                cerr<<eigenvectors<<endl;

#endif
		break;
	}
#ifdef USE_BROKEN_LANCZOS
		case 2:
		{
			//Experimental TODO
#ifdef DEBUG
			fprintf(stderr, "Calc mean \n");
#endif
			AlgebraTools::calculateMean(features, mean);
#ifdef DEBUG
			fprintf(stderr, "Lanczos \n");
#endif
			EigValues *eig = new EVLanczos(100, 1e-4);
			GMCovariance C(&features.transpose());
			if(adaptive)
			{
				eig->getEigenvalues(C, eigenvalue, eigenvectors, srcFeatureSize);
			}
			else
			{
				eig->getEigenvalues(C, eigenvalue, eigenvectors, mindimension); //fast lanczos
			}
#ifdef DEBUG
			fprintf(stderr, "PCA: Computation ready \n");
#endif
			break;
		}
#endif
	case 3:
	{

#ifdef DEBUG
		fprintf(stderr, "Calc mean \n");
#endif
		AlgebraTools::calculateMean(features, mean);
#ifdef DEBUG
		fprintf(stderr, "Arnoldi \n");
#endif
		EigValues *eig = new EVArnoldi(100, 1e-4);
		NICE::Matrix tmppointer = features.transpose();
		GMCovariance C(&tmppointer);
		if (adaptive)
		{
			eig->getEigenvalues(C, eigenvalue, eigenvectors, srcFeatureSize);
		}
		else
		{
			eig->getEigenvalues(C, eigenvalue, eigenvectors, mindimension); //fast arnoldi
		}
#ifdef DEBUG
		fprintf(stderr, "PCA: Computation ready \n");
#endif

#ifdef DEBUG
		fprintf(stderr, "PCA: Transpose ready \n");
#endif
		break;
	}
	default:
		fprintf(stderr, "PCA: Wrong Computation mode: %d!\n", mode);
		break;

	}
#ifdef DEBUG
	fprintf(stderr, "Eigenvalue Decomposition ready \n");
	cerr << eigenvectors << endl;
	cerr << eigenvalue << endl;

	//sort values
	fprintf(stderr, "Eigenvector-Rows:%i Eigenvector-Cols:%i\n", (int)eigenvectors.rows(), (int)eigenvectors.cols());
#endif

	multimap<double, NICE::Vector> map;
	double sumeigen = 0.0;

	NICE::Vector ratio(eigenvectors.cols());
	for (uint i = 0; i < eigenvectors.cols(); i++)//every eigenvector
	{
		NICE::Vector eigenvector(srcFeatureSize);
		for (uint k = 0; k < srcFeatureSize; k++)
		{
			eigenvector[k] = eigenvectors(k, i);
		}
		map.insert(pair<double, NICE::Vector> (eigenvalue[i], eigenvector));
		sumeigen += eigenvalue[i];
	}

	//compute basis size


	if (adaptive)
	{ //compute target dimension
		uint dimensioncount = 0;
		double addedratio = 0.0;
		multimap<double, NICE::Vector>::reverse_iterator it = map.rbegin();
		while (addedratio <= targetRatio && it != map.rend())
		{
			//calc ratio
			ratio[dimensioncount] = (*it).first / sumeigen;
			addedratio += ratio[dimensioncount];
			dimensioncount++;
			it++;
		}
		this->targetDimension = dimensioncount;
	}

	mindimension = std::min(this->targetDimension, srcFeatureSize);
	this->targetDimension = mindimension;
	basis = NICE::Matrix(srcFeatureSize, mindimension);
	//get sorted values
	uint count = 0;
	multimap<double, NICE::Vector>::reverse_iterator it = map.rbegin();
	while (count < this->targetDimension && it != map.rend())
	{
		NICE::Vector eigenvector = (*it).second;
		//put eigenvector into column
		for (uint k = 0; k < srcFeatureSize; k++)
		{
			basis(k, count) = eigenvector[k];
		}
		//calc ratio
		ratio[count] = (*it).first / sumeigen;

		count++;
		it++;
	}
	//normalization matrix / modify variance to 1 for all eigenvectors
	normalization = NICE::Matrix(mindimension, mindimension, 0);
	for (uint k = 0; k < mindimension; k++)
	{
		normalization(k, k) = 1.0 / sqrt(eigenvalue[k]);
	}

#ifdef DEBUG
	cerr << basis << endl;
	//     if (mode == 0)
	//     {
	//         cout << eigenvectors << endl;
	//         cout << basis << endl;
	//         cout << mean << endl;
	cout << "Eigenvalue-absolute:" << eigenvalue << endl;
	cout << "Eigenvalue-ratio:" << ratio << endl;
	//     }
#endif

}

NICE::Vector PCA::getFeatureVector(const NICE::Vector &data,
		const bool normalize)
{
	//free data of mean
	if (normalize)
	{

		NICE::Vector meanfree(data);
		meanfree -= mean;
		//Y=W^t * B^T
		if(normbasis.rows() == 0)
			normbasis.multiply(normalization, basis, false, true);
		NICE::Vector tmp;
		tmp.multiply(normbasis, meanfree);
		return tmp;
	}
	else
	{
		NICE::Vector meanfree(data);
		meanfree -= mean;
		//Y=W^t * B^T
		NICE::Vector tmp;
		tmp.multiply(basis, meanfree, true);
		return tmp;
	}
}

NICE::Vector PCA::getMean()
{
	return mean;
}

NICE::Matrix PCA::getBasis()
{
	return basis;
}

int PCA::getTargetDim()
{
	return targetDimension;
}