kmeans.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 * GNU General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _KMEANS_H
#define _KMEANS_H

#include "global.h"
#include "math.h"

struct MyKMeansInfo
{
public:
  MyKMeansInfo()
      : isSet(false), cost(), firstIndex()
    {}
  
  bool isSet;
  MyNT cost;
  unsigned int firstIndex;
};

inline ostream& operator<<(ostream &ostr, const MyKMeansInfo &info)
{
  ostr << info.cost << " ";
  return(ostr);
}

typedef pointer_to_binary_function< MyNT, MyNT, MyNT > MyDistanceFunction;

struct distance_Euclidean: public binary_function< MyNT, MyNT, MyNT > 
{
  MyNT operator()(MyNT __x, MyNT __y) const
    {
      return(fabs(__x - __y));
    }
};

struct distance_squared_Euclidean: public binary_function< MyNT, MyNT, MyNT > 
{
  MyNT operator()(MyNT __x, MyNT __y) const
    {
      return((__x - __y)*(__x - __y));
    }
};



// class to store info on a cluster in MyKMeansAlgo
template < class MyPointType >
class MyKMeansClusterType
{
//  typename int MyPointType;
  
public:
  void add(MyPointType &point)
    {
      _points.insert(point);
    }
  
  void add(const vector< MyPointType > &points);
  
  void remove(MyPointType &point)
    {
      _points.erase(point);
    }
  
  void remove(const vector< MyPointType > &points);

  // compute the mean.
  void mean()
    {
         _mean = 0;
         typename _MyPtsType::const_iterator itr;
      for (itr = _points.begin(); itr != _points.end(); itr++)
        _mean += *itr;

//       _mean = _points[0];
//       for (unsigned int i = 1; i < _points.size(); i++)
//         _mean += _points[i];
      _mean = _mean/_points.size();
    }

  void variance()
    {
      _variance = 0;
         typename _MyPtsType::const_iterator itr;
      for (itr = _points.begin(); itr != _points.end(); itr++)
        _variance += (*itr - _mean)*(*itr - _mean);
//       _variance = (_points[0] - _mean)*(_points[0] - _mean);
//       for (unsigned int i = 1; i < _points.size(); i++)
//         _variance += _points[i];
      _variance = _variance/_points.size();

    }

  
  MyNT distance(MyPointType point)
  {
    return(::distance_Euclidean(mean(), point))();
  }
  
private:
  MyPointType _mean;
  MyPointType _variance;
  
  typedef set< MyPointType > _MyPtsType;
  set< MyPointType > _points;
  
};



// stuff that a kmeans algorithm should support.
template < class MyKMeansPointType >
class MyKMeansAlgo
{
//  friend class MyKMeansClusterType< class MyKMeansPointType >;
  

  
  // should this typename be a template parameter?
//  typename vector< MyNT > MyKMeansPointType;
//  typename MyNT MyKMeansPointType;
  
public:

  // default or constructor for a given # clusters.
  MyKMeansAlgo(unsigned int num = 0)
      : _points(), _numClusters(num), _clusters(_numClusters)
    {}
  
  void setPoints(vector< MyKMeansPointType > &points)
    {
      _points = vector< MyKMeansPointType >(points.size());
      for (unsigned int i = 0; i < points.size(); i++)
        _points[i] = points[i];

     }
  
  void setNumClusters(unsigned int num)
    {
      _numClusters = num;
    }
  
  unsigned int getNumClusters()
    {
      return(numClusters);
    }

  // initiliase. randomly assign points to clusters.
  void init()
    {
      // static just in case i call init() tonnes of times. recreating
      // MyRand<> is very slow.
      static MyRand< int > generator;
      int index;
      for (unsigned int i = 0; i < _points.size(); i++)
        {
          index = generator(_numClusters - 1);
//          _clusters[index].add(i);
          _clusters[index].add(_points[i]);
          _whichCluster[i] = index;
        }
      for (unsigned int j = 0; j < _numClusters; j++)
        {
          _clusters[j].mean();
          _clusters[j].variance();
        }
      
    }
  
  // compute the mean of each cluster
  void estep()
    {
      for (unsigned int i = 0; i < _numClusters; i++)
        {
          _clusters[i].mean();
          _clusters[i].variance();
        }
    }
  
  // redistribute the points to the best cluster. return true iff some
  // point changes cluster.
  bool mstep()
    {
      unsigned int bestCluster;
      MyNT bestDist, dist;
      bool pointMoved = false;
      for (unsigned int i = 0; i < _points.size(); i++)
        {
          bestCluster = 0;
          bestDist = _clusters[bestCluster].distance(_points[i]);;
          for (unsigned int j = 1; j < _numClusters; j++)
            {
              dist = _clusters[j].distance(_points[i]);
              if (dist < bestDist)
                {
                  bestCluster = j;
                  bestDist = dist;;
                }
            }
          // don't do anything if bestCluster is the current cluster.
          if (bestCluster == _whichCluster[i])
            continue;
          pointMoved = true;
          // move point to the right cluster.
//          _clusters[_whichCluster[i]].remove(i);
//          _clusters[bestCluster].add(i);
          _clusters[_whichCluster[i]].remove(_points[i]);
          _clusters[bestCluster].add(_points[i]);
          _whichCluster[i] = bestCluster;
        }
      return(pointMoved);
    }
  
  
  // the actual function. returns the score.
//   MyNT operator()
//     {}
  
  

private:
  // points to cluster.
  vector< MyKMeansPointType > _points;
  
  unsigned int _numClusters;

  // how do i represent the clusters.

  // i want to store only an index in each cluster.
//  vector< MyKMeansClusterType< unsigned int > > _clusters;

    vector< MyKMeansClusterType< MyKMeansPointType > > _clusters;

  // for each point, store cluster it belongs to.
  vector< unsigned int > _whichCluster;
  
};




#endif // _KMEANS_H