GaussianMixtureModels.cpp

#include "GaussianMixtureModels.h"

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

//Constructor,destructor
GaussianMixtureModels::GaussianMixtureModels(const UINT numClusters,const UINT minNumEpochs,const UINT maxNumEpochs,const Float minChange){
    
    this->numClusters = numClusters;
    this->minNumEpochs = minNumEpochs;
    this->maxNumEpochs = maxNumEpochs;
    this->minChange = minChange;
    
    numTrainingSamples = 0;
    numTrainingIterationsToConverge = 0;
    trained = false;
    
    classType = "GaussianMixtureModels";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG GaussianMixtureModels]");
    errorLog.setProceedingText("[ERROR GaussianMixtureModels]");
    trainingLog.setProceedingText("[TRAINING GaussianMixtureModels]");
    warningLog.setProceedingText("[WARNING GaussianMixtureModels]");
}

GaussianMixtureModels::GaussianMixtureModels(const GaussianMixtureModels &rhs){
    
    classType = "GaussianMixtureModels";
    clustererType = classType;
    debugLog.setProceedingText("[DEBUG GaussianMixtureModels]");
    errorLog.setProceedingText("[ERROR GaussianMixtureModels]");
    trainingLog.setProceedingText("[TRAINING GaussianMixtureModels]");
    warningLog.setProceedingText("[WARNING GaussianMixtureModels]");
    
    if( this != &rhs ){
        
        this->numTrainingSamples = rhs.numTrainingSamples;
        this->loglike = rhs.loglike;
        this->mu = rhs.mu;
        this->resp = rhs.resp;
        this->frac = rhs.frac;
        this->lndets = rhs.lndets;
        this->det = rhs.det;
        this->sigma = rhs.sigma;
        this->invSigma = rhs.invSigma;
        
        //Clone the Clusterer variables
        copyBaseVariables( (Clusterer*)&rhs );
    }
    
}

GaussianMixtureModels::~GaussianMixtureModels(){
}

GaussianMixtureModels& GaussianMixtureModels::operator=(const GaussianMixtureModels &rhs){
    
    if( this != &rhs ){
        
        this->numTrainingSamples = rhs.numTrainingSamples;
        this->loglike = rhs.loglike;
        this->mu = rhs.mu;
        this->resp = rhs.resp;
        this->frac = rhs.frac;
        this->lndets = rhs.lndets;
        this->det = rhs.det;
        this->sigma = rhs.sigma;
        this->invSigma = rhs.invSigma;
        
        //Clone the Clusterer variables
        copyBaseVariables( (Clusterer*)&rhs );
    }
    
    return *this;
}
    
bool GaussianMixtureModels::deepCopyFrom(const Clusterer *clusterer){
    
    if( clusterer == NULL ) return false;
    
    if( this->getClustererType() == clusterer->getClustererType() ){
        //Clone the GaussianMixtureModels values
        GaussianMixtureModels *ptr = (GaussianMixtureModels*)clusterer;
        
        this->numTrainingSamples = ptr->numTrainingSamples;
        this->loglike = ptr->loglike;
        this->mu = ptr->mu;
        this->resp = ptr->resp;
        this->frac = ptr->frac;
        this->lndets = ptr->lndets;
        this->det = ptr->det;
        this->sigma = ptr->sigma;
        this->invSigma = ptr->invSigma;
        
        //Clone the Clusterer variables
        return copyBaseVariables( clusterer );
    }
    return false;
}
    
bool GaussianMixtureModels::reset(){
    
    Clusterer::reset();
    
    numTrainingSamples = 0;
    loglike = 0;
    
    return true;
}

bool GaussianMixtureModels::clear(){
    
    Clusterer::clear();
    
    numTrainingSamples = 0;
    loglike = 0;
 mu.clear();
 resp.clear();
 frac.clear();
 lndets.clear();
 det.clear();
 sigma.clear();
 invSigma.clear();
    
    return true;
}

bool GaussianMixtureModels::train_(MatrixFloat &data){
    
    trained = false;
    
    //Clear any previous training results
    det.clear();
    invSigma.clear();
    numTrainingIterationsToConverge = 0;
    
    if( data.getNumRows() == 0 ){
        errorLog << "train_(MatrixFloat &data) - Training Failed! Training data is empty!" << std::endl;
        return false;
    }
    
    //Resize the variables
    numTrainingSamples = data.getNumRows();
    numInputDimensions = data.getNumCols();
    
    //Resize mu and resp
    mu.resize(numClusters,numInputDimensions);
    resp.resize(numTrainingSamples,numClusters);
    
    //Resize sigma
    sigma.resize(numClusters);
    for(UINT k=0; k<numClusters; k++){
        sigma[k].resize(numInputDimensions,numInputDimensions);
    }
    
    //Resize frac and lndets
    frac.resize(numClusters);
    lndets.resize(numClusters);
    
    //Scale the data if needed
    ranges = data.getRanges();
    if( useScaling ){
        for(UINT i=0; i<numTrainingSamples; i++){
            for(UINT j=0; j<numInputDimensions; j++){
                data[i][j] = scale(data[i][j],ranges[j].minValue,ranges[j].maxValue,0,1);
            }
        }
    }
    
    //Pick K random starting points for the inital guesses of Mu
    Random random;
    Vector< UINT > randomIndexs(numTrainingSamples);
    for(UINT i=0; i<numTrainingSamples; i++) randomIndexs[i] = i;
    for(UINT i=0; i<numClusters; i++){
        SWAP(randomIndexs[ i ],randomIndexs[ random.getRandomNumberInt(0,numTrainingSamples) ]);
    }
    for(UINT k=0; k<numClusters; k++){
        for(UINT n=0; n<numInputDimensions; n++){
            mu[k][n] = data[ randomIndexs[k] ][n];
        }
    }
    
    //Setup sigma and the uniform prior on P(k)
    for(UINT k=0; k<numClusters; k++){
        frac[k] = 1.0/Float(numClusters);
        for(UINT i=0; i<numInputDimensions; i++){
            for(UINT j=0; j<numInputDimensions; j++) sigma[k][i][j] = 0;
            sigma[k][i][i] = 1.0e-2;   //Set the diagonal to a small number
        }
    }
    
    loglike = 0;
    bool keepGoing = true;
    Float change = 99.9e99;
    UINT numIterationsNoChange = 0;
    VectorFloat u(numInputDimensions);
 VectorFloat v(numInputDimensions);
    
    while( keepGoing ){
        
        //Run the estep
        if( estep( data, u, v, change ) ){
            
            //Run the mstep
            mstep( data );
        
            //Check for convergance
            if( fabs( change ) < minChange ){
                if( ++numIterationsNoChange >= minNumEpochs ){
                    keepGoing = false;
                }
            }else numIterationsNoChange = 0;
            if( ++numTrainingIterationsToConverge >= maxNumEpochs ) keepGoing = false;
            
        }else{
            errorLog << "train_(MatrixFloat &data) - Estep failed at iteration " << numTrainingIterationsToConverge << std::endl;
            return false;
        }
    }
    
    //Compute the inverse of sigma and the determinants for prediction
    if( !computeInvAndDet() ){
        det.clear();
        invSigma.clear();
        errorLog << "train_(MatrixFloat &data) - Failed to compute inverse and determinat!" << std::endl;
        return false;
    }
    
    //Flag that the model was 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 GaussianMixtureModels::train_(ClassificationData &trainingData){
    MatrixFloat data = trainingData.getDataAsMatrixFloat();
    return train_( data );
}

bool GaussianMixtureModels::train_(UnlabelledData &trainingData){
    MatrixFloat data = trainingData.getDataAsMatrixFloat();
    return train_( data );
}
    
bool GaussianMixtureModels::predict_(VectorFloat &x){
    
    if( !trained ){
        return false;
    }
    
    if( x.getSize() != numInputDimensions ){
        return false;
    }
    
    if( useScaling ){
        for(UINT n=0; n<numInputDimensions; n++){
            x[n] = grt_scale(x[n], ranges[n].minValue, ranges[n].maxValue, 0.0, 1.0);
        }
    }
    
    Float sum = 0;
    Float dist = 0;
    UINT minIndex = 0;
    bestDistance = 0;
    predictedClusterLabel = 0;
    maxLikelihood = 0;
    if( clusterLikelihoods.size() != numClusters )
        clusterLikelihoods.resize( numClusters );
    if( clusterDistances.size() != numClusters )
        clusterDistances.resize( numClusters );
    
    for(UINT i=0; i<numClusters; i++){
        
        dist = gauss(x,i,det,mu,invSigma);
        
        clusterDistances[i] = dist;
        clusterLikelihoods[i] = dist;
        
        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 GaussianMixtureModels::saveModelToFile( std::fstream &file ) const{
    
    if( !file.is_open() ){
        errorLog << "saveModelToFile(string filename) - Failed to open file!" << std::endl;
        return false;
    }
    
    file << "GRT_GAUSSIAN_MIXTURE_MODELS_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 << "Mu:\n";
        for(UINT k=0; k<numClusters; k++){
            for(UINT n=0; n<numInputDimensions; n++){
                file << mu[k][n] << "\t";
            }
            file << std::endl;
        }
        
        file << "Sigma:\n";
        for(UINT k=0; k<numClusters; k++){
            for(UINT i=0; i<numInputDimensions; i++){
                for(UINT j=0; j<numInputDimensions; j++){
                    file << sigma[k][i][j] << "\t";
                }
            }
            file << std::endl;
        }
        
        file << "InvSigma:\n";
        for(UINT k=0; k<numClusters; k++){
            for(UINT i=0; i<numInputDimensions; i++){
                for(UINT j=0; j<numInputDimensions; j++){
                    file << invSigma[k][i][j] << "\t";
                }
            }
            file << std::endl;
        }
        
        file << "Det:\n";
        for(UINT k=0; k<numClusters; k++){
            file << det[k] << std::endl;
        }
    }
    
    return true;
    
}

bool GaussianMixtureModels::loadModelFromFile( std::fstream &file ){
    
    //Clear any previous model
    clear();
    
    std::string word;
    file >> word;
    if( word != "GRT_GAUSSIAN_MIXTURE_MODELS_FILE_V1.0" ){
        return false;
    }
    
    if( !loadClustererSettingsFromFile( file ) ){
        errorLog << "loadModelFromFile(fstream &file) - Failed to load cluster settings from file!" << std::endl;
        return false;
    }
    
    //Load the model
    if( trained ){
        
        //Setup the memory
        mu.resize(numClusters, numInputDimensions);
        sigma.resize(numClusters);
        invSigma.resize(numClusters);
        det.resize(numClusters);
        
        //Load mu
        file >> word;
        if( word != "Mu:" ){
            clear();
            errorLog << "loadModelFromFile(fstream &file) - Failed to load Mu!" << std::endl;
            return false;
        }
        for(UINT k=0; k<numClusters; k++){
            for(UINT n=0; n<numInputDimensions; n++){
                file >> mu[k][n];
            }
        }
        
        //Load Sigma
        file >> word;
        if( word != "Sigma:" ){
            clear();
            errorLog << "loadModelFromFile(fstream &file) - Failed to load Sigma!" << std::endl;
            return false;
        }
        for(UINT k=0; k<numClusters; k++){
            sigma[k].resize(numInputDimensions, numInputDimensions);
            for(UINT i=0; i<numInputDimensions; i++){
                for(UINT j=0; j<numInputDimensions; j++){
                    file >> sigma[k][i][j];
                }
            }
        }
        
        //Load InvSigma
        file >> word;
        if( word != "InvSigma:" ){
            clear();
            errorLog << "loadModelFromFile(fstream &file) - Failed to load InvSigma!" << std::endl;
            return false;
        }
        for(UINT k=0; k<numClusters; k++){
            invSigma[k].resize(numInputDimensions, numInputDimensions);
            for(UINT i=0; i<numInputDimensions; i++){
                for(UINT j=0; j<numInputDimensions; j++){
                    file >> invSigma[k][i][j];
                }
            }
        }
        
        //Load Det
        file >> word;
        if( word != "Det:" ){
            clear();
            errorLog << "loadModelFromFile(fstream &file) - Failed to load Det!" << std::endl;
            return false;
        }
        for(UINT k=0; k<numClusters; k++){
            file >> det[k];
        }
        
        //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 GaussianMixtureModels::estep( const MatrixFloat &data, VectorFloat &u, VectorFloat &v, Float &change ){

 Float tmp,sum,max,oldloglike;
 for(UINT j=0; j<numInputDimensions; j++) u[j] = v[j] = 0;

 oldloglike = loglike;

 for(UINT k=0; k<numClusters; k++){
  Cholesky cholesky( sigma[k] );
  if( !cholesky.getSuccess() ){ return false; }
  lndets[k] = cholesky.logdet();

  for(UINT i=0; i<numTrainingSamples; i++){
   for(UINT j=0; j<numInputDimensions; j++) u[j] = data[i][j] - mu[k][j];
   if( !cholesky.elsolve(u,v) ){ return false; }
   sum=0;
   for(UINT j=0; j<numInputDimensions; j++) sum += SQR(v[j]);
   resp[i][k] = -0.5*(sum + lndets[k]) + log(frac[k]);
  }
 }

 //Compute the overall likelihood of the entire estimated paramter set
 loglike = 0;
 for(UINT i=0; i<numTrainingSamples; i++){
  sum=0;
  max = -99.9e99;
  for(UINT k=0; k<numClusters; k++) if( resp[i][k] > max ) max = resp[i][k];
  for(UINT k=0; k<numClusters; k++) sum += exp( resp[i][k]-max );
  tmp = max + log( sum );
  for(UINT k=0; k<numClusters; k++) resp[i][k] = exp( resp[i][k] - tmp );
  loglike += tmp;
 }
    
    change = (loglike - oldloglike);

 return true;
}

bool GaussianMixtureModels::mstep( const MatrixFloat &data ){

 Float wgt, sum;
 for(UINT k=0; k<numClusters; k++){
  wgt = 0.0;
  for(UINT m=0; m<numTrainingSamples; m++) wgt += resp[m][k];
  frac[k] = wgt/Float(numTrainingSamples);
  for(UINT n=0; n<numInputDimensions; n++){
   sum = 0;
   for(UINT m=0; m<numTrainingSamples; m++) sum += resp[m][k] * data[m][n];
   mu[k][n] = sum/wgt;
   for(UINT j=0; j<numInputDimensions; j++){
    sum = 0;
    for(UINT m=0; m<numTrainingSamples; m++){
     sum += resp[m][k] * (data[m][n]-mu[k][n]) * (data[m][j]-mu[k][j]);
    }
    sigma[k][n][j] = sum/wgt;
   }
  }
 }
    
    return true;

}

inline void GaussianMixtureModels::SWAP(UINT &a,UINT &b){
 UINT temp = b;
 b = a;
 a = temp;
}

bool GaussianMixtureModels::computeInvAndDet(){

 det.resize(numClusters);
 invSigma.resize(numClusters);

 for(UINT k=0; k<numClusters; k++){
  LUDecomposition lu(sigma[k]);
  if( !lu.inverse( invSigma[k] ) ){
            errorLog << "computeInvAndDet() - Matrix inversion failed for cluster " << k+1 << std::endl;
            return false;
        }
  det[k] = lu.det();
 }

    return true;

}

GRT_END_NAMESPACE