KMeans.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 "KMeans.h"

GRT_BEGIN_NAMESPACE
    
//Register the KMeans class with the Clusterer base class
RegisterClustererModule< KMeans > KMeans::registerModule("KMeans");

//Constructor,destructor
KMeans::KMeans(const UINT numClusters,const UINT minNumEpochs,const UINT maxNumEpochs,const Float minChange,const bool computeTheta){
    
    this->numClusters = numClusters;
    this->minNumEpochs = minNumEpochs;
    this->maxNumEpochs = maxNumEpochs;
    this->minChange = minChange;
    this->computeTheta = computeTheta;
    
    numTrainingSamples = 0;
    nchg = 0;
    finalTheta = 0;
    numTrainingIterationsToConverge = 0;
    trained = false;
    
    classType = "KMeans";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG KMeans]");
    errorLog.setProceedingText("[ERROR KMeans]");
    trainingLog.setProceedingText("[TRAINING KMeans]");
    warningLog.setProceedingText("[WARNING KMeans]");
}
    
KMeans::KMeans(const KMeans &rhs){
    
    classType = "KMeans";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG KMeans]");
    errorLog.setProceedingText("[ERROR KMeans]");
    trainingLog.setProceedingText("[TRAINING KMeans]");
    warningLog.setProceedingText("[WARNING KMeans]");
    
    if( this != &rhs ){
        
        this->numTrainingSamples = rhs.numTrainingSamples;
        this->nchg = rhs.nchg;
        this->computeTheta = rhs.computeTheta;
        this->finalTheta = rhs.finalTheta;
        this->clusters = rhs.clusters;
        this->assign = rhs.assign;
        this->count = rhs.count;
        this->thetaTracker = rhs.thetaTracker;
        
        //Clone the Clusterer variables
        copyBaseVariables( (Clusterer*)&rhs );
    }
    
}

KMeans::~KMeans(){
}
    
KMeans& KMeans::operator=(const KMeans &rhs){
    
    if( this != &rhs ){
        
        this->numTrainingSamples = rhs.numTrainingSamples;
        this->nchg = rhs.nchg;
        this->computeTheta = rhs.computeTheta;
        this->finalTheta = rhs.finalTheta;
        this->clusters = rhs.clusters;
        this->assign = rhs.assign;
        this->count = rhs.count;
        this->thetaTracker = rhs.thetaTracker;
        
        //Clone the Clusterer variables
        copyBaseVariables( (Clusterer*)&rhs );
    }
    
    return *this;
}

bool KMeans::deepCopyFrom(const Clusterer *clusterer){
    
    if( clusterer == NULL ) return false;
    
    if( this->getClustererType() == clusterer->getClustererType() ){
        //Clone the KMeans values
        KMeans *ptr = (KMeans*)clusterer;
        
        this->numTrainingSamples = ptr->numTrainingSamples;
        this->nchg = ptr->nchg;
        this->computeTheta = ptr->computeTheta;
        this->finalTheta = ptr->finalTheta;
        this->clusters = ptr->clusters;
        this->assign = ptr->assign;
        this->count = ptr->count;
        this->thetaTracker = ptr->thetaTracker;
        
        //Clone the Clusterer variables
        return copyBaseVariables( clusterer );
    }
    return false;
}

bool KMeans::train_(ClassificationData &trainingData){
 
 if( trainingData.getNumSamples() == 0 ){
        errorLog << "train_(ClassificationData &trainingData) - The training data is empty!" << std::endl;
  return false;
 }
 
 //Set the numClusters as the number of classes in the training data
 numClusters = trainingData.getNumClasses();

    //Convert the labelled training data to a training matrix
 UINT M = trainingData.getNumSamples();
    UINT N = trainingData.getNumDimensions();
    MatrixFloat data(M,N);
    for(UINT i=0; i<M; i++){
        for(UINT j=0; j<N; j++){
            data[i][j] = trainingData[i][j];
        }
    }

    //Run the K-Means algorithm
    return train_( data );
}

bool KMeans::train_(UnlabelledData &trainingData){

    //Convert the training data into one matrix
 UINT M = trainingData.getNumSamples();
    UINT N = trainingData.getNumDimensions();
    MatrixFloat data(M,N);
    for(UINT i=0; i<M; i++){
        for(UINT j=0; j<N; j++){
            data[i][j] = trainingData[i][j];
        }
    }
 
 return train_(data);
}

bool KMeans::train_(MatrixFloat &data){
 
 trained = false;
 
 if( numClusters == 0 ){
        errorLog << "train_(MatrixFloat &data) - Failed to train model. NumClusters is zero!" << std::endl;
  return false;
 }
    
    if( data.getNumRows() == 0 || data.getNumCols() == 0 ){
        errorLog << "train_(MatrixFloat &data) - The number of rows or columns in the data is zero!" << std::endl;
  return false;
 }
    
 numTrainingSamples = data.getNumRows();
 numInputDimensions = data.getNumCols();

 clusters.resize(numClusters,numInputDimensions);
 assign.resize(numTrainingSamples);
 count.resize(numClusters);

 //Randomly pick k data points as the starting clusters
 Random random;
 Vector< UINT > randIndexs(numTrainingSamples);
 for(UINT i=0; i<numTrainingSamples; i++) randIndexs[i] = i;
    std::random_shuffle(randIndexs.begin(), randIndexs.end());

    //Copy the clusters
 for(UINT k=0; k<numClusters; k++){
  for(UINT j=0; j<numInputDimensions; j++){
            clusters[k][j] = data[ randIndexs[k] ][j];
  }
 }

 return trainModel( data );
}
    
bool KMeans::predict_(VectorFloat &inputVector){
    
    if( !trained ){
        return false;
 }
 
 if( inputVector.getSize() != numInputDimensions ){
  return false;
 }
    
    if( useScaling ){
        for(UINT n=0; n<numInputDimensions; n++){
            inputVector[n] = grt_scale(inputVector[n], ranges[n].minValue, ranges[n].maxValue, 0.0, 1.0);
        }
    }
 
    const Float sigma = 1.0;
    const Float gamma = 1.0 / (2.0*grt_sqr(sigma));
    Float sum = 0;
    Float dist = 0;
 UINT minIndex = 0;
 bestDistance = grt_numeric_limits< Float >::max();
 predictedClusterLabel = 0;
 maxLikelihood = 0;
 if( clusterLikelihoods.getSize() != numClusters )
        clusterLikelihoods.resize( numClusters );
    if( clusterDistances.getSize() != numClusters )
        clusterDistances.resize( numClusters );
 
 for(UINT i=0; i<numClusters; i++){
  
        //We don't need to compute the sqrt as it works without it and is faster
  dist = 0;
  for(UINT j=0; j<numInputDimensions; j++){
   dist += grt_sqr( inputVector[j]-clusters[i][j] );
  }
    
        clusterDistances[i] = dist;
        clusterLikelihoods[i] = exp( - grt_sqr(gamma * dist) ); //1.0/(1.0+dist); //This will give us a value close to 1 for a dist of 0, and a value closer to 0 when the dist is large
        
  sum += clusterLikelihoods[i];
        
  if( dist < bestDistance ){
   bestDistance = dist;
   minIndex = i;
  }
 }
 
 //Normalize the likelihood
 for(UINT i=0; i<numClusters; i++){
  clusterLikelihoods[i] /= sum;
 }
 
 predictedClusterLabel = clusterLabels[ minIndex ];
 maxLikelihood = clusterLikelihoods[ minIndex ];
    
    return true;
}

bool KMeans::trainModel(MatrixFloat &data){
    
    if( numClusters == 0 ){
        errorLog << "trainModel(MatrixFloat &data) - Failed to train model. NumClusters is zero!" << std::endl;
  return false;
 }
    
    if( clusters.getNumRows() != numClusters ){
        errorLog << "trainModel(MatrixFloat &data) - Failed to train model. The number of rows in the cluster matrix does not match the number of clusters! You should need to initalize the clusters matrix first before calling this function!" << std::endl;
  return false;
 }
    
    if( clusters.getNumCols() != numInputDimensions ){
        errorLog << "trainModel(MatrixFloat &data) - Failed to train model. The number of columns in the cluster matrix does not match the number of input dimensions! You should need to initalize the clusters matrix first before calling this function!" << std::endl;
  return false;
 }

    Timer timer;
 UINT currentIter = 0;
    UINT numChanged = 0;
 bool keepTraining = true;
    Float theta = 0;
    Float lastTheta = 0;
    Float delta = 0;
    Float startTime = 0;
    thetaTracker.clear();
    finalTheta = 0;
    numTrainingIterationsToConverge = 0;
    trained = false;
    converged = false;
    
    //Scale the data if needed
    ranges = data.getRanges();
    if( useScaling ){
        data.scale(0,1);
    }

    //Init the assign and count Vectors
    //Assign is set to K+1 so that the nChanged values in the eStep at the first iteration will be updated correctly
    for(UINT m=0; m<numTrainingSamples; m++) assign[m] = numClusters+1;
 for(UINT k=0; k<numClusters; k++) count[k] = 0;

    //Run the training loop
    timer.start();
 while( keepTraining ){
        startTime = timer.getMilliSeconds();

  //Compute the E step
  numChanged = estep( data );

        //Compute the M step
        mstep( data );

        //Update the iteration counter
  currentIter++;

  //Compute theta if needed
  if( computeTheta ){
            theta = calculateTheta(data);
            delta = lastTheta - theta;
            lastTheta = theta;
        }else theta = delta = 0;
        
        //Check convergance
  if( numChanged == 0 && currentIter > minNumEpochs ){ converged = true; keepTraining = false; }
  if( currentIter >= maxNumEpochs ){ keepTraining = false; }
  if( fabs( delta ) < minChange && computeTheta && currentIter > minNumEpochs ){ converged = true; keepTraining = false; }
        if( computeTheta )  thetaTracker.push_back( theta );
        
        trainingLog << "Epoch: " << currentIter << "/" << maxNumEpochs;
        trainingLog << " Epoch time: " << (timer.getMilliSeconds()-startTime)/1000.0 << " seconds";
        trainingLog << " Theta: " << theta << " Delta: " << delta << std::endl;
 }
    trainingLog << "Model Trained at epoch: " << currentIter << " with a theta value of: " << theta << std::endl;

    finalTheta = theta;
    numTrainingIterationsToConverge = currentIter;
 trained = true;
    
    //Setup the cluster labels
    clusterLabels.resize(numClusters);
    for(UINT i=0; i<numClusters; i++){
        clusterLabels[i] = i+1;
    }
    clusterLikelihoods.resize(numClusters,0);
    clusterDistances.resize(numClusters,0);
 
 return true;
}

UINT KMeans::estep(const MatrixFloat &data) {
  UINT k,m,n,kmin;
  Float dmin,d;
  nchg = 0;
  kmin = 0;

  //Reset Count
  for (k=0; k < numClusters; k++) count[k] = 0;

  //Search for the closest center and reasign if needed
  for (m=0; m < numTrainingSamples; m++) {
   dmin = 9.99e+99; //Set dmin to a really big value
   for (k=0; k < numClusters; k++) {
    d = 0.0;
    for (n=0; n < numInputDimensions; n++)
     d += grt_sqr( data[m][n]-clusters[k][n] );
    if (d <= dmin){ dmin = d; kmin = k; }
   }
   if ( kmin != assign[m] ){
                nchg++;
                assign[m] = kmin;
            }
   count[kmin]++;
  }
  return nchg;
}

void KMeans::mstep(const MatrixFloat &data) {
    UINT n,k,m;

    //Reset means to zero
    for (k=0; k<numClusters; k++)
        for (n=0;n<numInputDimensions;n++)
            clusters[k][n] = 0.;

    //Get new mean by adding assigned data points and dividing by the number of values in each cluster
    for(m=0; m < numTrainingSamples; m++)
        for(n=0; n < numInputDimensions; n++)
            clusters[ assign[m] ][n] += data[m][n];

    for (k=0; k < numClusters; k++) {
        if (count[k] > 0){
            Float countNorm = 1.0 / count[k];
            for (n=0; n < numInputDimensions; n++){
                clusters[k][n] *= countNorm;
            }
        }
    }
}

Float KMeans::calculateTheta(const MatrixFloat &data){

 Float theta = 0;
    Float sum = 0;
    UINT m,n,k = 0;
 for(m=0; m < numTrainingSamples; m++){
  k = assign[m];
        sum = 0;
  for(n=0; n < numInputDimensions; n++){
    sum += grt_sqr(clusters[k][n] - data[m][n]);
  }
  theta += grt_sqrt(sum);
 }
    theta /= numTrainingSamples;

 return theta;

}

bool KMeans::saveModelToFile( std::fstream &file ) const{

    if( !file.is_open() ){
        errorLog << "saveModelToFile(fstream &file) - Failed to save model, file is not open!" << std::endl;
        return false;
    }

    file << "GRT_KMEANS_MODEL_FILE_V1.0\n";
    
    if( !saveClustererSettingsToFile( file ) ){
        errorLog << "saveModelToFile(fstream &file) - Failed to save clusterer settings to file!" << std::endl;
        return false;
    }
    
    if( trained ){
        file << "Clusters:\n";

        for(UINT k=0; k<numClusters; k++){
            for(UINT n=0; n<numInputDimensions; n++){
                file << clusters[k][n] << "\t";
           }file << std::endl;
        }
    }

   return true;

}

bool KMeans::loadModelFromFile( std::fstream &file ){

    //Clear any previous model
    clear();
    
    if(!file.is_open()){
        errorLog << "loadModelFromFile(string filename) - Failed to open file!" << std::endl;
        return false;
    }

    std::string word;
    file >> word;
    if( word != "GRT_KMEANS_MODEL_FILE_V1.0" ){
    return false;
    }
    
    if( !loadClustererSettingsFromFile( file ) ){
        errorLog << "loadModelFromFile(string filename) - Failed to open file!" << std::endl;
        return false;
    }

    if( trained ){
        file >> word;
        if( word != "Clusters:" ){
            return false;
        }
        
        //Resize the buffers
        clusters.resize(numClusters,numInputDimensions);
        
        //Load the data
        for(UINT k=0; k<numClusters; k++){
            for(UINT n=0; n<numInputDimensions; n++){
                file >> clusters[k][n];
            }
        }
    }

    return true;
}
    
bool KMeans::reset(){
    Clusterer::reset();
    
    numTrainingSamples = 0;
    nchg = 0;
    finalTheta = 0;
    thetaTracker.clear();
    assign.clear();
    count.clear();
    
    return true;
}

bool KMeans::clear(){
    Clusterer::clear();
    
    numTrainingSamples = 0;
    nchg = 0;
    finalTheta = 0;
    thetaTracker.clear();
    assign.clear();
    count.clear();
    clusters.clear();
    
    return true;
}
    
bool KMeans::setComputeTheta(const bool computeTheta){
    this->computeTheta = computeTheta;
    return true;
}
    
bool KMeans::setClusters(const MatrixFloat &clusters){
    clear();
    numClusters = clusters.getNumRows();
    numInputDimensions = clusters.getNumCols();
    this->clusters = clusters;
    return true;
}

GRT_END_NAMESPACE