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

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

HierarchicalClustering::HierarchicalClustering(){
    M = N = 0;
    classType = "HierarchicalClustering";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG HierarchicalClustering]");
    errorLog.setProceedingText("[ERROR HierarchicalClustering]");
    trainingLog.setProceedingText("[TRAINING HierarchicalClustering]");
    warningLog.setProceedingText("[WARNING HierarchicalClustering]");
}
    
HierarchicalClustering::HierarchicalClustering(const HierarchicalClustering &rhs){
    classType = "HierarchicalClustering";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG HierarchicalClustering]");
    errorLog.setProceedingText("[ERROR HierarchicalClustering]");
    trainingLog.setProceedingText("[TRAINING HierarchicalClustering]");
    warningLog.setProceedingText("[WARNING HierarchicalClustering]");
    *this = rhs;
}

HierarchicalClustering::~HierarchicalClustering(){

}
    
HierarchicalClustering& HierarchicalClustering::operator=(const HierarchicalClustering &rhs){
    
    if( this != &rhs ){
        
        this->M = rhs.M;
        this->N = rhs.N;
        this->clusters = rhs.clusters;
        this->distanceMatrix = rhs.distanceMatrix;
        
        //Clone the Clusterer variables
        copyBaseVariables( (Clusterer*)&rhs );
    }
    
    return *this;
}
    
bool HierarchicalClustering::deepCopyFrom(const Clusterer *clusterer){
    
    if( clusterer == NULL ) return false;
    
    if( this->getClustererType() == clusterer->getClustererType() ){
        //Clone the HierarchicalClustering values
        HierarchicalClustering *ptr = (HierarchicalClustering*)clusterer;
        
        this->M = ptr->M;
        this->N = ptr->N;
        this->clusters = ptr->clusters;
        this->distanceMatrix = ptr->distanceMatrix;
        
        //Clone the Clusterer variables
        return copyBaseVariables( clusterer );
    }
    return false;
}

bool HierarchicalClustering::reset(){
    
    Clusterer::reset();
    
    return true;
}

bool HierarchicalClustering::clear(){
    
    Clusterer::clear();
    
    M = 0;
    N = 0;
    clusters.clear();
    distanceMatrix.clear();
    
    return true;
}

bool HierarchicalClustering::train_(ClassificationData &trainingData){
 
 if( trainingData.getNumSamples() == 0 ){
  return false;
 }

    //Convert the labelled training data to a training matrix
 M = trainingData.getNumSamples();
    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 HierarchicalClustering::train_(UnlabelledData &trainingData){
    
    if( trainingData.getNumSamples() == 0 ){
  return false;
 }

    //Convert the training data into one matrix
 M = trainingData.getNumSamples();
    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 HierarchicalClustering::train_(MatrixFloat &data){
 
 trained = false;
    clusters.clear();
    distanceMatrix.clear();
    
    if( data.getNumRows() == 0 || data.getNumCols() == 0 ){
  return false;
 }
 
    //Set the rows and columns
    M = data.getNumRows();
 N = data.getNumCols();
    
    //Build the distance matrix
    distanceMatrix.resize(M,M);

    //Build the distance matrix
    for(UINT i=0; i<M; i++){
        for(UINT j=0; j<M; j++){
            if( i== j ) distanceMatrix[i][j] = grt_numeric_limits< Float >::max();
            else{
                distanceMatrix[i][j] = squaredEuclideanDistance(data[i], data[j]);
            }
        }
    }

    //Build the initial clusters, at the start each sample gets its own cluster
    UINT uniqueClusterID = 0;
    Vector< ClusterInfo > clusterData(M);
    for(UINT i=0; i<M; i++){
        clusterData[i].uniqueClusterID = uniqueClusterID++;
        clusterData[i].addSampleToCluster(i);
    }
    
    trainingLog << "Starting clustering..." << std::endl;
    
    //Create the first cluster level, each sample is it's own cluster
    UINT level = 0;
    ClusterLevel newLevel;
    newLevel.level = level;
    for(UINT i=0; i<M; i++){
        newLevel.clusters.push_back( clusterData[i] );
    }
    clusters.push_back( newLevel );
    
    //Move to level 1 and start the search
    level++;
    bool keepClustering = true;
    
    while( keepClustering ){
        
        //Find the closest two clusters within the cluster data
        Float minDist = grt_numeric_limits< Float >::max();
        Vector< Vector< UINT > > clusterPairs;
        UINT K = (UINT)clusterData.size();
        for(UINT i=0; i<K; i++){
            for(UINT j=0; j<K; j++){
                if( i != j ){
                    Float dist = computeClusterDistance( clusterData[i], clusterData[j]  );
             
                    if( dist < minDist ){
                        minDist = dist;
                        Vector< UINT > clusterPair(2);
                        clusterPair[0] = i;
                        clusterPair[1] = j;
                        clusterPairs.clear();
                        clusterPairs.push_back( clusterPair );
                    }
                    
                }
            }
        }
        
        if( minDist == grt_numeric_limits< Float >::max() ){
            keepClustering = false;
            warningLog << "train_(MatrixFloat &data) - Failed to find any cluster at level: " << level << std::endl;
            return false;
        }else{
        
            //Merge the two closest clusters together and create a new level
            ClusterLevel newLevel;
            newLevel.level = level;
            
            //Create the new cluster
            ClusterInfo newCluster;
            newCluster.uniqueClusterID = uniqueClusterID++;
            
            const UINT numClusterPairs = clusterPairs.getSize();
            
            for(UINT k=0; k<numClusterPairs; k++){
                //Add all the samples in the first cluster to the new cluster
                UINT numSamplesInClusterA = clusterData[ clusterPairs[k][0] ].getNumSamplesInCluster();
                for(UINT i=0; i<numSamplesInClusterA; i++){
                    UINT index = clusterData[ clusterPairs[k][0] ][ i ];
                    newCluster.addSampleToCluster( index );
                }
                
                //Add all the samples in the second cluster to the new cluster
                UINT numSamplesInClusterB = clusterData[ clusterPairs[k][1] ].getNumSamplesInCluster();
                for(UINT i=0; i<numSamplesInClusterB; i++){
                    UINT index = clusterData[ clusterPairs[k][1] ][ i ];
                    newCluster.addSampleToCluster( index );
                }
                
                //Compute the cluster variance
                newCluster.clusterVariance = computeClusterVariance( newCluster, data );
                
                //Remove the two cluster pairs (so they will not be used in the next search
                UINT idA = clusterData[ clusterPairs[k][0] ].getUniqueClusterID();
                UINT idB = clusterData[ clusterPairs[k][1] ].getUniqueClusterID();
                UINT numRemoved = 0;
                Vector< ClusterInfo >::iterator iter = clusterData.begin();
                while( iter != clusterData.end() ){
                    if( iter->getUniqueClusterID() == idA || iter->getUniqueClusterID() == idB ){
                        iter = clusterData.erase( iter );
                        if( ++numRemoved >= 2 ) break;
                    }else iter++;
                }
            }
            
            //Add the merged cluster to the clusterData
            clusterData.push_back( newCluster );
            
            //Add the new level and cluster data to the main cluster buffer
            newLevel.clusters.push_back( newCluster );
            
            clusters.push_back( newLevel );
            
            //Update the level
            level++;
        }
        
        //Check to see if we should stop clustering
        if( level >= M ){
            keepClustering = false;
        }
        
        if( clusterData.size() == 0 ){
            keepClustering = false;
        }
        
        trainingLog << "Cluster level: " << level << " Number of clusters: " << clusters.back().getNumClusters() << std::endl;
    }
    
    //Flag that the model is trained
    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;
}
    
bool HierarchicalClustering::printModel(){
    
    UINT K = (UINT)clusters.size();
    
    std::cout << "Hierarchical Clustering Model\n\n";
    for(UINT k=0; k<K; k++){
        UINT numClusters = clusters[k].getNumClusters();
        UINT numSamples = 0;
        for(UINT i=0; i<numClusters; i++){
            numSamples += clusters[k][i].getNumSamplesInCluster();
        }
        
        std::cout << "Level: " << clusters[k].level << "\tNumClusters: " << numClusters << "\tNumSamples: " << numSamples << std::endl;
        for(UINT i=0; i<numClusters; i++){
             std::cout << "ClusterVariance: " << clusters[k][i].clusterVariance << std::endl;
             std::cout << "Indexs: ";
            UINT numSamplesInCluster = clusters[k][i].getNumSamplesInCluster();
            for(UINT j=0; j<numSamplesInCluster; j++){
                 std::cout << clusters[k][i][j] << "\t";
            }
             std::cout << std::endl;
        }
    }

    return true;
}
    
Float HierarchicalClustering::squaredEuclideanDistance(const Float *a,const Float *b){
    Float dist = 0;
    for(UINT i=0; i<N; i++){
        dist += SQR( a[i] - b[i] );
    }
    return dist;
}
    
Float HierarchicalClustering::computeClusterDistance( const ClusterInfo &clusterA, const ClusterInfo &clusterB ){
    
    Float minDist = grt_numeric_limits< Float >::max();
    const UINT numSamplesA = clusterA.getNumSamplesInCluster();
    const UINT numSamplesB = clusterB.getNumSamplesInCluster();
    
    //Find the minimum distance between the two clusters
    for(UINT i=0; i<numSamplesA; i++){
        for(UINT j=0; j<numSamplesB; j++){
            if( distanceMatrix[ clusterA[i] ][ clusterB[j] ] < minDist ){
                minDist = distanceMatrix[ clusterA[i] ][ clusterB[j] ];
            }
        }
    }
    
    return minDist;
}
    
Float HierarchicalClustering::computeClusterVariance( const ClusterInfo &cluster, const MatrixFloat &data ){
    
    VectorFloat mean(N,0);
    VectorFloat std(N,0);
    
    //Compute the mean
    UINT numSamples = cluster.getNumSamplesInCluster();
    for(UINT j=0; j<N; j++){
        for(UINT i=0; i<numSamples; i++){
            UINT index = cluster[i];
            mean[j] += data[ index ][j];
        }
        mean[j] /= Float( numSamples );
    }
    
    //Compute the std dev
    for(UINT j=0; j<N; j++){
        for(UINT i=0; i<numSamples; i++){
            std[j] += grt_sqr( data[ cluster[i] ][j] - mean[j] );
        }
        std[j] = grt_sqrt( std[j] / Float( numSamples-1 ) );
    }
    
    Float variance = 0;
    for(UINT j=0; j<N; j++){
        variance += std[j];
    }
    return variance/N;
}

bool HierarchicalClustering::saveModelToFile( std::fstream &file ) const{
    
    if( !file.is_open() ){
        errorLog << "saveModelToFile(string filename) - Failed to open file!" << std::endl;
        return false;
    }
    
    file << "GRT_HIERARCHICAL_CLUSTERING_FILE_V1.0\n";
    
    if( !saveClustererSettingsToFile( file ) ){
        errorLog << "saveModelToFile(fstream &file) - Failed to save cluster settings to file!" << std::endl;
        return false;
    }
    
    if( trained ){
        file << "M: " << M << std::endl;
        file << "N: " << N << std::endl;
        file << "NumLevels: " << clusters.getSize() << std::endl;
        
        for(UINT i=0; i<clusters.getSize(); i++){
            file << "Level: " << clusters[i].getLevel() << std::endl;
            file << "NumClusters: " << clusters[i].getNumClusters() << std::endl;
        }
    }
    
    return true;
    
}

bool HierarchicalClustering::loadModelFromFile( std::fstream &file ){
    
    std::string word;
    
    //Clear any previous model
    clear();
    
    file >> word;
    if( word != "GRT_HIERARCHICAL_CLUSTERING_FILE_V1.0" ){
        return false;
    }
    
    if( !loadClustererSettingsFromFile( file ) ){
        errorLog << "loadModelFromFile(fstream &file) - Failed to load cluster settings from file!" << std::endl;
        return false;
    }
        
    return true;
}

GRT_END_NAMESPACE