PrincipalComponentAnalysis.cpp

Go to the documentation of this file.

#define GRT_DLL_EXPORTS
#include "PrincipalComponentAnalysis.h"

GRT_BEGIN_NAMESPACE

PrincipalComponentAnalysis::PrincipalComponentAnalysis() : MLBase("PrincipalComponentAnalysis")
{
 trained = false;
 normData = false;
 numInputDimensions = 0;
 numPrincipalComponents = 0;
 maxVariance = 0;
}

PrincipalComponentAnalysis::~PrincipalComponentAnalysis(){
 
}

bool PrincipalComponentAnalysis::computeFeatureVector(const MatrixFloat &data,double maxVariance,bool normData){
 trained = false;
 this->maxVariance = maxVariance;
 this->normData = normData;
 return computeFeatureVector_(data,MAX_VARIANCE);
}

bool PrincipalComponentAnalysis::computeFeatureVector(const MatrixFloat &data,UINT numPrincipalComponents,bool normData){
 trained = false;
 if( numPrincipalComponents > data.getNumCols() ){
 errorLog << __GRT_LOG__ << " The number of principal components (";
 errorLog << numPrincipalComponents << ") is greater than the number of columns in your data (" << data.getNumCols() << ")" << std::endl;
 return false;
 }
 this->numPrincipalComponents = numPrincipalComponents;
 this->normData = normData;
 return computeFeatureVector_(data,MAX_NUM_PCS);
}

bool PrincipalComponentAnalysis::computeFeatureVector_(const MatrixFloat &data,const UINT analysisMode){
 
 trained = false;
 const UINT M = data.getNumRows();
 const UINT N = data.getNumCols();
 this->numInputDimensions = N;
 
 MatrixFloat msData( M, N );
 
 //Compute the mean and standard deviation of the input data
 mean = data.getMean();
 stdDev = data.getStdDev();
 
 if( normData ){
 //Normalize the data
 for(UINT i=0; i<M; i++)
 for(UINT j=0; j<N; j++)
 msData[i][j] = (data[i][j]-mean[j]) / stdDev[j];
 
 }else{
 //Mean Subtract Data
 for(UINT i=0; i<M; i++)
 for(UINT j=0; j<N; j++)
 msData[i][j] = data[i][j] - mean[j];
 }
 
 //Get the covariance matrix
 MatrixFloat cov = msData.getCovarianceMatrix();
 
 //Use Eigen Value Decomposition to find eigenvectors of the covariance matrix
 EigenvalueDecomposition eig;
 
 if( !eig.decompose( cov ) ){
 mean.clear();
 stdDev.clear();
 componentWeights.clear();
 sortedEigenvalues.clear();
 eigenvectors.clear();
 errorLog << __GRT_LOG__ << " Failed to decompose input matrix!" << std::endl;
 return false;
 }
 
 //Get the eigenvectors and eigenvalues
 eigenvectors = eig.getEigenvectors();
 eigenvalues = eig.getRealEigenvalues();
 
 //Any eigenvalues less than 0 are not worth anything so set to 0
 for(UINT i=0; i<eigenvalues.size(); i++){
 if( eigenvalues[i] < 0 )
 eigenvalues[i] = 0;
 }
 
 //Sort the eigenvalues and compute the component weights
 Float sum = 0;
 UINT componentIndex = 0;
 sortedEigenvalues.clear();
 componentWeights.resize(N,0);
 
 while( true ){
 Float maxValue = 0;
 UINT index = 0;
 for(UINT i=0; i<eigenvalues.size(); i++){
 if( eigenvalues[i] > maxValue ){
 maxValue = eigenvalues[i];
 index = i;
 }
 }
 if( maxValue == 0 || componentIndex >= eigenvalues.size() ){
 break;
 }
 sortedEigenvalues.push_back( IndexedDouble(index,maxValue) );
 componentWeights[ componentIndex++ ] = eigenvalues[ index ];
 sum += eigenvalues[ index ];
 eigenvalues[ index ] = 0; //Set the maxValue to zero so it won't be used again
 }
 
 Float cumulativeVariance = 0;
 switch( analysisMode ){
 case MAX_VARIANCE:
 //Normalize the component weights and workout how many components we need to use to reach the maxVariance
 numPrincipalComponents = 0;
 for(UINT k=0; k<N; k++){
 componentWeights[k] /= sum;
 cumulativeVariance += componentWeights[k];
 if( cumulativeVariance >= maxVariance && numPrincipalComponents==0 ){
 numPrincipalComponents = k+1;
 }
 }
 break;
 case MAX_NUM_PCS:
 //Normalize the component weights and compute the maxVariance
 maxVariance = 0;
 for(UINT k=0; k<N; k++){
 componentWeights[k] /= sum;
 if( k < numPrincipalComponents ){
 maxVariance += componentWeights[k];
 }
 }
 break;
 default:
 errorLog << __GRT_LOG__ << " Unknown analysis mode!" << std::endl;
 break;
 }
 
 //Get the raw eigenvalues (encase the user asks for these later)
 eigenvalues = eig.getRealEigenvalues();
 
 //Flag that the features have been computed
 trained = true;
 
 return true;
}

bool PrincipalComponentAnalysis::project(const MatrixFloat &data,MatrixFloat &prjData){
 
 if( !trained ){
 warningLog << __GRT_LOG__ << " The PrincipalComponentAnalysis module has not been trained!" << std::endl;
 return false;
 }
 
 if( data.getNumCols() != numInputDimensions ){
 warningLog << __GRT_LOG__ << " The number of columns in the input vector (" << data.getNumCols() << ") does not match the number of input dimensions (" << numInputDimensions << ")!" << std::endl;
 return false;
 }
 
 MatrixFloat msData( data );
 prjData.resize(data.getNumRows(),numPrincipalComponents);
 
 if( normData ){
 //Mean subtract the data
 for(UINT i=0; i<data.getNumRows(); i++)
 for(UINT j=0; j<numInputDimensions; j++)
 msData[i][j] = (msData[i][j]-mean[j])/stdDev[j];
 }else{
 //Mean subtract the data
 for(UINT i=0; i<data.getNumRows(); i++)
 for(UINT j=0; j<numInputDimensions; j++)
 msData[i][j] -= mean[j];
 }
 
 //Projected Data
 for(UINT row=0; row<msData.getNumRows(); row++){//For each row in the final data
 for(UINT i=0; i<numPrincipalComponents; i++){//For each PC
 prjData[row][i]=0;
 for(UINT j=0; j<data.getNumCols(); j++)//For each feature
 prjData[row][i] += msData[row][j] * eigenvectors[j][sortedEigenvalues[i].index];
 }
 }
 
 return true;
}

bool PrincipalComponentAnalysis::project(const VectorFloat &data,VectorFloat &prjData){
 
 const unsigned int N = data.getSize();
 
 if( !trained ){
 warningLog << __GRT_LOG__ << " The PrincipalComponentAnalysis module has not been trained!" << std::endl;
 return false;
 }
 
 if( N != numInputDimensions ){
 warningLog << __GRT_LOG__ << " The size of the input vector (" << N << ") does not match the number of input dimensions (" << numInputDimensions << ")!" << std::endl;
 return false;
 }
 
 VectorFloat msData = data;
 
 if( normData ){
 //Mean subtract the data
 for(UINT j=0; j<numInputDimensions; j++)
 msData[j] = (msData[j]-mean[j])/stdDev[j];
 }else{
 //Mean subtract the data
 for(UINT j=0; j<numInputDimensions; j++)
 msData[j] -= mean[j];
 }
 
 //Projected Data
 prjData.resize( numPrincipalComponents );
 for(UINT i=0; i<numPrincipalComponents; i++){//For each PC
 prjData[i]=0;
 for(UINT j=0; j<N; j++)//For each feature
 prjData[i] += msData[j] * eigenvectors[j][sortedEigenvalues[i].index];
 }
 
 return true;
}

bool PrincipalComponentAnalysis::save( std::fstream &file ) const {
 
 //Write the header info
 file << "GRT_PCA_MODEL_FILE_V1.0\n";
 
 if( !MLBase::saveBaseSettingsToFile( file ) ) return false;
 
 file << "NumPrincipalComponents: " << numPrincipalComponents << std::endl;
 file << "NormData: " << normData << std::endl;
 file << "MaxVariance: " << maxVariance << std::endl;
 
 if( trained ){
 file << "Mean: ";
 for(unsigned int i=0; i<numInputDimensions; i++){
 file << mean[i] << " ";
 }
 file << std::endl;
 
 file << "StdDev: ";
 for(unsigned int i=0; i<numInputDimensions; i++){
 file << stdDev[i] << " ";
 }
 file << std::endl;
 
 file << "ComponentWeights: ";
 for(unsigned int i=0; i<numInputDimensions; i++){
 file << componentWeights[i] << " ";
 }
 file << std::endl;
 
 file << "Eigenvalues: ";
 for(unsigned int i=0; i<numInputDimensions; i++){
 file << eigenvalues[i] << " ";
 }
 file << std::endl;
 
 file << "SortedEigenvalues: ";
 for(unsigned int i=0; i<numInputDimensions; i++){
 file << sortedEigenvalues[i].index << " ";
 file << sortedEigenvalues[i].value << " ";
 }
 file << std::endl;
 
 file << "Eigenvectors: ";
 file << eigenvectors.getNumRows() << " " << eigenvectors.getNumCols() << std::endl;
 for(unsigned int i=0; i<eigenvectors.getNumRows(); i++){
 for(unsigned int j=0; j<eigenvectors.getNumCols(); j++){
 file << eigenvectors[i][j];
 if( j+1 < eigenvectors.getNumCols() ) file << " ";
 else file << std::endl;
 }
 }
 file << std::endl;
 }
 
 return true;
}

bool PrincipalComponentAnalysis::load( std::fstream &file ) {
 
 std::string word;
 
 //Read the header info
 file >> word;
 if( word != "GRT_PCA_MODEL_FILE_V1.0" ){
 return false;
 }
 
 if( !MLBase::loadBaseSettingsFromFile( file ) ) return false;
 
 //Read the num components
 file >> word;
 if( word != "NumPrincipalComponents:" ){
 return false;
 }
 file >> numPrincipalComponents;
 
 //Read the normData
 file >> word;
 if( word != "NormData:" ){
 return false;
 }
 file >> normData;
 
 //Read the MaxVariance
 file >> word;
 if( word != "MaxVariance:" ){
 return false;
 }
 file >> maxVariance;
 
 if( trained ){
 //Read the mean vector
 file >> word;
 if( word != "Mean:" ){
 trained = false;
 return false;
 }
 mean.resize( numInputDimensions );
 
 for(unsigned int i=0; i<numInputDimensions; i++){
 file >> mean[i];
 }
 
 //Read the stddev vector
 file >> word;
 if( word != "StdDev:" ){
 trained = false;
 return false;
 }
 stdDev.resize( numInputDimensions );
 
 for(unsigned int i=0; i<numInputDimensions; i++){
 file >> stdDev[i];
 }
 
 //Read the ComponentWeights vector
 file >> word;
 if( word != "ComponentWeights:" ){
 trained = false;
 return false;
 }
 componentWeights.resize( numInputDimensions );
 
 for(unsigned int i=0; i<numInputDimensions; i++){
 file >> componentWeights[i];
 }
 
 //Read the Eigenvalues vector
 file >> word;
 if( word != "Eigenvalues:" ){
 trained = false;
 return false;
 }
 eigenvalues.resize( numInputDimensions );
 
 for(unsigned int i=0; i<numInputDimensions; i++){
 file >> eigenvalues[i];
 }
 
 //Read the SortedEigenvalues vector
 file >> word;
 if( word != "SortedEigenvalues:" ){
 trained = false;
 return false;
 }
 sortedEigenvalues.resize( numInputDimensions );
 
 for(unsigned int i=0; i<numInputDimensions; i++){
 file >> sortedEigenvalues[i].index;
 file >> sortedEigenvalues[i].value;
 }
 
 //Read the Eigenvectors vector
 file >> word;
 if( word != "Eigenvectors:" ){
 trained = false;
 return false;
 }
 UINT numRows;
 UINT numCols;
 file >> numRows;
 file >> numCols;
 eigenvectors.resize( numRows, numCols );
 
 for(unsigned int i=0; i<eigenvectors.getNumRows(); i++){
 for(unsigned int j=0; j<eigenvectors.getNumCols(); j++){
 file >> eigenvectors[i][j];
 }
 }
 }
 
 return true;
}

bool PrincipalComponentAnalysis::print( std::string title ) const{
 
 if( title != "" ){
 std::cout << title << std::endl;
 }
 if( !trained ){
 std::cout << "Not Trained!\n";
 return false;
 }
 std::cout << "NumInputDimensions: " << numInputDimensions << " NumPrincipalComponents: " << numPrincipalComponents << std::endl;
 std::cout << "ComponentWeights: ";
 for(UINT k=0; k<componentWeights.size(); k++){
 std::cout << "\t" << componentWeights[k];
 }
 std::cout << std::endl;
 std::cout << "SortedEigenValues: ";
 for(UINT k=0; k<sortedEigenvalues.size(); k++){
 std::cout << "\t" << sortedEigenvalues[k].value;
 }
 std::cout << std::endl;
 eigenvectors.print("Eigenvectors:");
 
 return true;
}

MatrixFloat PrincipalComponentAnalysis::getEigenVectors() const{
 return eigenvectors;
}

bool PrincipalComponentAnalysis::setModel( const VectorFloat &mean, const MatrixFloat &eigenvectors ){
 
 if( (UINT)mean.size() != eigenvectors.getNumCols() ){
 return false;
 }
 
 trained = true;
 numInputDimensions = eigenvectors.getNumCols();
 numPrincipalComponents = eigenvectors.getNumRows();
 this->mean = mean;
 stdDev.clear();
 componentWeights.clear();
 eigenvalues.clear();
 sortedEigenvalues.clear();
 this->eigenvectors = eigenvectors;
 
 //The eigenvectors are already sorted, so the sorted eigenvalues just holds the default index
 for(UINT i=0; i<numPrincipalComponents; i++){
 sortedEigenvalues.push_back( IndexedDouble(i,0.0) );
 }
 return true;
}

GRT_END_NAMESPACE