op_kmeans_meat.hpp

Go to the documentation of this file.

// Copyright (C) 2011 the authors listed below
// http://ecocpak.sourceforge.net
// 
// Authors:
// - Dimitrios Bouzas (bouzas at ieee dot org)
// - Nikolaos Arvanitopoulos (niarvani at ieee dot org)
// - Anastasios Tefas (tefas at aiia dot csd dot auth dot gr)
// 
// This file is part of the ECOC PAK C++ library. It is 
// provided without any warranty of fitness for any purpose.
//
// You can redistribute this file and/or modify it under 
// the terms of the GNU Lesser General Public License (LGPL) 
// as published by the Free Software Foundation, either 
// version 3 of the License or (at your option) any later 
// version.
// (see http://www.opensource.org/licenses for more info)





template<typename eT>
inline
eT
op_kmeans::compute_distance
  (
  const subview_row<eT>& sample_in, 
  const subview_row<eT>& centroid_in
  )
  {
  arma_extra_debug_sigprint();

  Row<eT> temp_vec = sample_in - centroid_in;
  return accu(temp_vec % temp_vec);   
  }



template<typename eT>
inline
void
op_kmeans::compute_distance
  (
        Row<eT>& dist_out, 
            u32& winner, 
  const Mat<eT>& centroids_in, 
  const subview_row<eT>& sample_in
  )
  {
  arma_extra_debug_sigprint();
  
  const u32 k = centroids_in.n_rows;
  const u32 n_cols = centroids_in.n_cols;

  dist_out = zeros<Row<eT> >(k); 
  eT winning_dist = numeric_limits<eT>::max();
  winner = 0;

  // Iterate through centroids
  for(u32 i = 0; i < k; ++i)
    {
    Row<eT> temp_row = centroids_in.row(i) - sample_in;      
    dist_out(i) =  accu(temp_row % temp_row);
        
    if(dist_out(i) < winning_dist)
      {
      winning_dist = dist_out(i);
      winner = i;
      }       
    }

  }  



template<typename eT>
inline
void
op_kmeans::init_centroids
  (
        Mat<eT>& centroids_out, 
  const Mat<eT>& samples_in, 
  const u32 k
  )
  {
  arma_extra_debug_sigprint();
  
  const u32 n_samples = samples_in.n_rows;
  const u32 n_cols    = samples_in.n_cols;
  
  // Allocate memory for centroids matrix
  centroids_out = zeros<Mat<eT> >(k, n_cols);

  // --- Kmeans++ algorithm introduced by David Arthur & Sergei Vassilvitskii initialization --- //
  
  // Initialize current potential      
  eT currentPot = 0;
  
  // Allocate distances vector
  Col<eT> closestDistSq = zeros<Col<eT> >(n_samples);

  // Choose one random center and set the closestDistSq values
  u32 index = (eT(std::rand()) / eT(RAND_MAX)) * n_samples;
  centroids_out.row(0) = samples_in.row(index);
  
  // Compute distances and potentials
  for(u32 i = 0; i < n_samples; i++) 
    {
    closestDistSq[i] = compute_distance(samples_in.row(i), centroids_out.row(0));
    currentPot += closestDistSq[i];
    }

  // Repeat for the rest k - 1 centers to be chosen 
  for(u32 centerCount = 1; centerCount < k; centerCount++) 
    {
    // Repeat several trials
    eT bestNewPot = -1;
    u32 bestNewIndex;
    
    // Compute number of local trials
    u32 numLocalTries = 2 + std::log(eT(k));
    
    for(u32 localTrial = 0; localTrial < numLocalTries; localTrial++) 
      {
                // Choose center - have to be slightly careful to return a valid answer even accounting
          // for possible rounding errors
      eT randVal =  (eT(std::rand()) / eT(RAND_MAX)) * currentPot;
    
      for(index = 0; index < n_samples - 1; index++) 
        {
        if (randVal <= closestDistSq[index])
          {
          break;
          }
        else
          {
          randVal -= closestDistSq[index];
          }
        }

                  // Compute the new potential
      eT newPot = 0;
      
      // Iterate through samples
      for(u32 i = 0; i < n_samples; i++)
        {
        newPot += std::min(compute_distance(samples_in.row(i), samples_in.row(index)), closestDistSq[i]);
        }

      // Store the best result
      if(bestNewPot < 0 || newPot < bestNewPot) 
        {
        bestNewPot = newPot;
        bestNewIndex = index;
        }
            
      }

    // Add the appropriate center
    centroids_out.row(centerCount) = samples_in.row(bestNewIndex);
    currentPot = bestNewPot;

    for(u32 i = 0; i < n_samples; i++)
      {
      closestDistSq[i] = std::min(compute_distance(samples_in.row(i), samples_in.row(bestNewIndex)), closestDistSq[i]);
      }
  
    }
  
  }



template<typename eT>
inline
void
op_kmeans::direct_kmeans
  (
       Col<u32>& indices_out,
       Col<u32>& ranks, 
        Mat<eT>& centroids_out, 
  const Mat<eT>& samples_in,
  const u32      k, 
  const eT       lr
  )
  {
  arma_extra_debug_sigprint();
  
  #ifdef NCURSES_OUTPUT
  //cout << endl << "decoding Started on " << ctime(&tbegin);
  mvaddstr(3, 0, "Computing K-means:");
  refresh();
  #endif

  const u32 n_samples = samples_in.n_rows;
  const u32 n_cols    = samples_in.n_cols;
  
  // Check wether the user provided any samples
  arma_check(n_samples == 0, "kmeans(): Input matrix contains no samples");
  
  // Check wether the user specified the number of clusters to be greater or equal with the number of samples
  // In this case return as centroids the samples
  if(k >= n_samples)
    {
    centroids_out = samples_in;

    ranks = zeros<Col<u32> >(k);
    
    if(n_samples > 1)
      {
      indices_out = conv_to<Col<u32> >::from(linspace(0, n_samples - 1, n_samples));
      }
    else
      {
      indices_out = zeros<Col<u32> >(1);
      }
      
    return;
    }
  
  // Allocate memory for output arguments
  indices_out = zeros<Col<u32> >(n_samples);
      
  // Initialize centroids
  init_centroids(centroids_out, samples_in, k);

  // Declare distances matrix
  Mat<eT> dist_mat;

  // Allocate centintel indices  
  Col<u32> prv_inds = ones<Col<u32> >(n_samples);
  Col<u32> more_prv_inds = zeros<Col<u32> >(n_samples);
      
  // Iterations counter
  u32 iters = 1;
  
  // Repeat while there's change in centroids
  while(accu(prv_inds == indices_out) < n_samples) 
    {
    // Update centintel indices
    more_prv_inds = prv_inds;
    prv_inds = indices_out;
          
    // Initialize sum of distances and distances matrix
    dist_mat = zeros< Mat<eT> >(n_samples, k);  

    // Allocate ranks vector for each centroid -- i.e., vectors that counts the winning sample for each centroid
    ranks = zeros<Col<u32> >(k);
      
    // Sample accumulator -- accumulates winning samples for each centroid for later update
    Mat<eT> sample_accu = zeros<Mat<eT> >(k, n_cols);
    
    // Temporary vector holding the distances of the current sample from each one of the centroids
    // TODO: Can't pass dist_out.row(i) by reference to compute_distance() that's why we need this
    //       temporary row vector.
    Row<eT> dist_tmp;
     
    // Compute distances of each samples from the available centroids
    for(u32 i = 0; i < n_samples; ++i)
      {
      compute_distance(dist_tmp, indices_out(i), centroids_out, samples_in.row(i));
      dist_mat.row(i) = dist_tmp;
      sample_accu.row(indices_out(i)) += samples_in.row(i);
      ranks(indices_out(i))++;  

      // Update winning centroid by moving it towards current sample -- Actually this is the only difference with the batch variant
      // centroids_out.row(indices_out(i)) += ((eT(lr) / eT(iters)) * (samples_in.row(i) - centroids_out.row(indices_out(i))));  
       centroids_out.row(indices_out(i)) += (eT(lr) * (samples_in.row(i) - centroids_out.row(indices_out(i))));
      
      #ifdef NCURSES_OUTPUT
      mvaddch(3, 20, spin_chars[i & 3]); 
      refresh();
      #endif
      }       

    // Update centroids 
    centroids_out = sample_accu;
  
    // Normalize updated centroids according to number of winning samples
    for(u32 i = 0; i < k; ++i)
      {
      centroids_out.row(i) /= eT(ranks(i));
      
      #ifdef NCURSES_OUTPUT
      mvaddch(3, 20, spin_chars[i & 3]); 
      refresh();
      #endif
      }

    // Check for empty clusters
    if(min(ranks) == 0)
      {
      // Find the empty cluster -- centroid
      
      Col<u32> empty_indx;
      Col<u32> max_indx;
      
      // Iterate while there are empty clusters       
      empty_indx = sort_index(ranks);
      
      u32 k = 0;
      while(min(ranks) == 0)
        { 
        max_indx = sort_index(dist_mat.col(empty_indx(k)),1);             
        centroids_out.row(empty_indx(k)) = samples_in.row(max_indx(k));
        ranks(empty_indx(k)) = 1;
        indices_out(max_indx(k)) = empty_indx(k); 
        
        #ifdef NCURSES_OUTPUT
        mvaddch(3, 20, spin_chars[k & 3]); 
        refresh();
        #endif
        
        k++;             
        }
        
      // Check wether the system oscilates
      if(accu(more_prv_inds == indices_out) < n_samples)
        {
        prv_inds = indices_out;
        }
   
      }
     
    #ifdef NCURSES_OUTPUT
    mvaddch(3, 20, spin_chars[iters & 3]); 
    refresh();
    #endif
    
    // Increase number of iterations
    iters++;         
    } 
 
  }