LDA.cpp

/*
GRT MIT License
Copyright (c) <2012> <Nicholas Gillian, Media Lab, MIT>

Permission is hereby granted, free of charge, to any person obtaining a copy of this software 
and associated documentation files (the "Software"), to deal in the Software without restriction, 
including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, 
and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, 
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial 
portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT 
LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE 
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#include "LDA.h"

GRT_BEGIN_NAMESPACE

LDA::LDA(bool useScaling,bool useNullRejection,Float nullRejectionCoeff)
{
    this->useScaling = useScaling;
    this->useNullRejection = useNullRejection;
    this->nullRejectionCoeff = nullRejectionCoeff;
    classType = "LDA";
    classifierType = classType;
    classifierMode = STANDARD_CLASSIFIER_MODE;
    debugLog.setProceedingText("[DEBUG LDA]");
    errorLog.setProceedingText("[ERROR LDA]");
    trainingLog.setProceedingText("[TRAINING LDA]");
    warningLog.setProceedingText("[WARNING LDA]");
}

LDA::~LDA(void)
{
} 
    
bool LDA::train(ClassificationData trainingData){
    
    errorLog << "SORRY - this module is still under development and can't be used yet!" << std::endl;
    return false;
    
    //Reset any previous model
    numInputDimensions = 0;
    numClasses = 0;
    models.clear();
    classLabels.clear();
    trained = false;
    
    if( trainingData.getNumSamples() == 0 ){
        errorLog << "train(LabelledClassificationData trainingData) - There is no training data to train the model!" << std::endl;
        return false;
    }
    
    numInputDimensions = trainingData.getNumDimensions();
    numClasses = trainingData.getNumClasses();

 //Calculate the between scatter matrix
 MatrixFloat SB = computeBetweenClassScatterMatrix( trainingData );
 
 //Calculate the within scatter matrix
 MatrixFloat SW = computeWithinClassScatterMatrix( trainingData );


   /*

    
    //Counters and stat containers
    vector< UINT > groupLabels(numClasses);
    VectorDouble groupCounters(numClasses);
    VectorDouble priorProb(numClasses);
    MatrixFloat groupMeans(numClasses,numFeatures);
    MatrixFloat pCov(numFeatures,numFeatures);
    MatrixFloat pCovInv(numFeatures,numFeatures);
    MatrixFloat modelCoeff(numClasses,numFeatures+1);
    
    pCov.setAllValues(0);
    modelCoeff.setAllValues(0);
    
    //Set the class labels and counters
    for(UINT k=0; k<numClasses; k++){
        groupLabels[k] = trainingData.getClassTracker()[k].classLabel;
        groupCounters[k] = trainingData.getClassTracker()[k].counter;
    }
    
    //Loop over the classes to compute the group stats
    for(UINT k=0; k<numClasses; k++){
        LabelledClassificationData classData = trainingData.getClassData( groupLabels[k] );
        MatrixFloat cov(numFeatures,numFeatures);
        
        //Compute class mu
        for(UINT j=0; j<numFeatures; j++){
            groupMeans[k][j] = 0;
            for(UINT i=0; i<classData.getNumSamples(); i++){
                groupMeans[k][j] += classData[i][j];
            }
            groupMeans[k][j] /= Float(classData.getNumSamples());
        }
        
        //Compute the class covariance
        for(UINT m=0; m<numFeatures; m++){
            for(UINT n=0; n<numFeatures; n++){
                cov[m][n] = 0;
                for(UINT i=0; i<classData.getNumSamples(); i++){
                    cov[m][n] += (classData[i][m]-groupMeans[k][m]) * (classData[i][n]-groupMeans[k][n]);
                }
                cov[m][n] /= Float(classData.getNumSamples()-1);
            }
        }
        
        debugLog << "Group Cov:\n";
        for(UINT m=0; m<numFeatures; m++){
            for(UINT n=0; n<numFeatures; n++){
                debugLog << cov[m][n] << "\t";
            }debugLog << "\n";
        }debugLog << std::endl;
        
        //Set the prior probability for this class (which is just 1/numClasses)
        priorProb[k] = 1.0/Float(numClasses);
        
        //Update the main covariance matrix
        Float weight = ((classData.getNumSamples() - 1) / Float(trainingData.getNumSamples() - numClasses) );
        debugLog << "Weight: " << weight << std::endl;
        for(UINT m=0; m<numFeatures; m++){
            for(UINT n=0; n<numFeatures; n++){
                pCov[m][n] += weight * cov[m][n];
            }
        }
    }
    
    for(UINT k=0; k<numClasses; k++){
        debugLog << "GroupMu: " << groupLabels[k] << "\t";
        for(UINT j=0; j<numFeatures; j++){
            debugLog << groupMeans[k][j] << "\t";
        }debugLog << std::endl;
    }
    
    debugLog << "pCov:\n";
    for(UINT m=0; m<numFeatures; m++){
        for(UINT n=0; n<numFeatures; n++){
            debugLog << pCov[m][n] << "\t";
        }debugLog << "\n";
    }debugLog << std::endl;
    
    //Invert the pCov matrix
    LUDecomposition matrixInverter(pCov);
    if( !matrixInverter.inverse(pCovInv) ){
        errorLog << "Failed to invert pCov Matrix!" << std::endl;
        return false;
    }
    
    //Loop over classes to calculate linear discriminant coefficients
    Float sum = 0;
    vector< Float > temp(numFeatures);
    for(UINT k=0; k<numClasses; k++){
        //Compute the temporary vector
        for(UINT j=0; j<numFeatures; j++){
            temp[j] = 0;
            for(UINT m=0; m<numFeatures; m++){
                    temp[j] += groupMeans[k][m] * pCovInv[m][j];
            }
        }
        
        //Compute the model coefficients
        sum = 0;
        for(UINT j=0; j<numFeatures; j++){
            sum += temp[j]*groupMeans[k][j];
        }
        modelCoeff[k][0] = -0.5 * sum + log( priorProb[k] );
        
        for(UINT j=0; j<numFeatures; j++){
            modelCoeff[k][j+1] = temp[j];
        }
    }
    
    //Setup the models for realtime prediction
    models.resize(numClasses);
    classLabels.resize(numClasses);
    
    for(UINT k=0; k<numClasses; k++){
        classLabels[k] = groupLabels[k];
        models[k].classLabel = groupLabels[k];
        models[k].priorProb = priorProb[k];
        models[k].weights = modelCoeff.getRowVector(k);
    }
    
    //Flag that the models were successfully trained
    trained = true;
    */
    
    return true;
}

bool LDA::predict(VectorDouble inputVector){
    
    if( !trained ){
        errorLog << "predict(vector< Float > inputVector) - LDA Model Not Trained!" << std::endl;
        return false;
    }
    
    predictedClassLabel = 0;
 maxLikelihood = -10000;
    
    if( !trained ) return false;
    
 if( inputVector.getSize() != numInputDimensions ){
        errorLog << "predict(vector< Float > inputVector) - The size of the input vector (" << inputVector.getSize() << ") does not match the num features in the model (" << numInputDimensions << std::endl;
  return false;
 }
    
    //Make sure the likelihoods and distances vectors have been assigned
    if( classLikelihoods.getSize() != numClasses || classDistances.getSize() != numClasses ){
        classLikelihoods.resize(numClasses);
        classDistances.resize(numClasses);
    }
    
    //Compute the linear scores for each class
    bestDistance = 0;
    maxLikelihood = 0;
    UINT bestIndex = 0;
    Float sum = 0;
    for(UINT k=0; k<numClasses; k++){
        
        for(UINT j=0; j<numInputDimensions+1; j++){
            if( j==0 ) classDistances[k] = models[k].weights[j];
            else classDistances[k] += inputVector[j-1] * models[k].weights[j];
        }
        classLikelihoods[k] = exp( classDistances[k] );
        sum += classLikelihoods[k];
        
        if( classLikelihoods[k] > maxLikelihood ){
            bestIndex = k;
            maxLikelihood = classLikelihoods[k];
        }
    }
    
    //Normalize the likelihoods
    for(UINT k=0; k<numClasses; k++){
        classLikelihoods[k] /= sum;
    }
    
    maxLikelihood = classLikelihoods[ bestIndex ];
    
    predictedClassLabel = models[ bestIndex ].classLabel;
    
    return true;
}
    
bool LDA::saveModelToFile( std::fstream &file ) const{
    
    if(!file.is_open())
    {
        errorLog <<"saveModelToFile(fstream &file) - Could not open file to save model" << std::endl;
        return false;
    }
    
    //Write the header info
    file<<"GRT_LDA_MODEL_FILE_V1.0\n";
    file<<"NumFeatures: "<<numInputDimensions<< std::endl;
    file<<"NumClasses: "<<numClasses<< std::endl;
    file <<"UseScaling: " << useScaling << std::endl;
    file<<"UseNullRejection: " << useNullRejection << std::endl;
    
    if( useScaling ){
        file << "Ranges: \n";
        for(UINT n=0; n<ranges.size(); n++){
            file << ranges[n].minValue << "\t" << ranges[n].maxValue << std::endl;
        }
    }
    
    //Write each of the models
    for(UINT k=0; k<numClasses; k++){
        file<<"ClassLabel: "<<models[k].classLabel<< std::endl;
        file<<"PriorProbability: "<<models[k].priorProb<< std::endl;
        file<<"Weights: ";
        
        for(UINT j=0; j<models[k].getNumDimensions(); j++){
            file << "\t" << models[k].weights[j];
        }file<< std::endl;
    }
    
    return true;
}
    
bool LDA::loadModelFromFile( std::fstream &file ){
    
    trained = false;
    numInputDimensions = 0;
    numClasses = 0;
    models.clear();
    classLabels.clear();
    
    if(!file.is_open())
    {
        errorLog << "loadModelFromFile(fstream &file) - The file is not open!" << std::endl;
        return false;
    }
    
    std::string word;
    
    //Find the file type header
    file >> word;
    if(word != "GRT_LDA_MODEL_FILE_V1.0"){
        errorLog << "loadModelFromFile(fstream &file) - Could not find Model File Header" << std::endl;
        return false;
    }
    
    file >> word;
    if(word != "NumFeatures:"){
        errorLog << "loadModelFromFile(fstream &file) - Could not find NumFeatures " << std::endl;
        return false;
    }
    file >> numInputDimensions;
    
    file >> word;
    if(word != "NumClasses:"){
        errorLog << "loadModelFromFile(fstream &file) - Could not find NumClasses" << std::endl;
        return false;
    }
    file >> numClasses;
    
    file >> word;
    if(word != "UseScaling:"){
        errorLog << "loadModelFromFile(fstream &file) - Could not find UseScaling" << std::endl;
        return false;
    }
    file >> useScaling;
    
    file >> word;
    if(word != "UseNullRejection:"){
        errorLog << "loadModelFromFile(fstream &file) - Could not find UseNullRejection" << std::endl;
        return false;
    }
    file >> useNullRejection;
    
    if( useScaling ){
        //Resize the ranges buffer
        ranges.resize(numInputDimensions);
        
        file >> word;
        if(word != "Ranges:"){
            errorLog << "loadModelFromFile(fstream &file) - Could not find the Ranges" << std::endl;
            return false;
        }
        for(UINT n=0; n<ranges.size(); n++){
            file >> ranges[n].minValue;
            file >> ranges[n].maxValue;
        }
    }
    
    //Resize the buffer
    models.resize(numClasses);
    classLabels.resize(numClasses);
    
    //Load each of the K models
    for(UINT k=0; k<numClasses; k++){
        file >> word;
        if(word != "ClassLabel:"){
            errorLog << "loadModelFromFile(fstream &file) - Could not find ClassLabel for the "<<k+1<<"th model" << std::endl;
            return false;
        }
        file >> models[k].classLabel;
        classLabels[k] = models[k].classLabel;
        
        file >> word;
        if(word != "PriorProbability:"){
            errorLog << "loadModelFromFile(fstream &file) - Could not find the PriorProbability for the "<<k+1<<"th model" << std::endl;
            return false;
        }
        file >> models[k].priorProb;
        
        models[k].weights.resize(numInputDimensions+1);
        
        //Load the weights
        file >> word;
        if(word != "Weights:"){
            errorLog << "loadModelFromFile(fstream &file) - Could not find the Weights vector for the "<<k+1<<"th model" << std::endl;
            return false;
        }
        
        //Load Weights
        for(UINT j=0; j<numInputDimensions+1; j++){
            Float value;
            file >> value;
            models[k].weights[j] = value;
        }
    }
    
    //Resize the prediction results to make sure it is setup for realtime prediction
    maxLikelihood = DEFAULT_NULL_LIKELIHOOD_VALUE;
    bestDistance = DEFAULT_NULL_DISTANCE_VALUE;
    classLikelihoods.resize(numClasses,DEFAULT_NULL_LIKELIHOOD_VALUE);
    classDistances.resize(numClasses,DEFAULT_NULL_DISTANCE_VALUE);
    
    trained = true;
    
    return true;
}

MatrixFloat LDA::computeBetweenClassScatterMatrix( ClassificationData &data ){
 
 MatrixFloat sb(numInputDimensions,numInputDimensions);
 MatrixFloat classMean = data.getClassMean();
 VectorDouble totalMean = data.getMean();
 sb.setAllValues( 0 );
 
 for(UINT k=0; k<numClasses; k++){
  
  UINT numSamplesInClass = data.getClassTracker()[k].counter;
 
  for(UINT m=0; m<numInputDimensions; m++){
   for(UINT n=0; n<numInputDimensions; n++){
    sb[m][n] += (classMean[k][m]-totalMean[m]) * (classMean[k][n]-totalMean[n]) * Float(numSamplesInClass);
   }
  }
 }
 
 return sb;
}

MatrixFloat LDA::computeWithinClassScatterMatrix( ClassificationData &data ){
 
 MatrixFloat sw(numInputDimensions,numInputDimensions);
 sw.setAllValues( 0 );
 
 for(UINT k=0; k<numClasses; k++){
  
  //Compute the scatter matrix for class k
  ClassificationData classData = data.getClassData( data.getClassTracker()[k].classLabel );
  MatrixFloat scatterMatrix = classData.getCovarianceMatrix();
  
  //Add this to the main scatter matrix
  for(UINT m=0; m<numInputDimensions; m++){
   for(UINT n=0; n<numInputDimensions; n++){
    sw[m][n] += scatterMatrix[m][n];
   }
  }
 }
 
 return sw;
}

GRT_END_NAMESPACE