projcluster.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 _PROJCLUSTER_H
#define _PROJCLUSTER_H

#include <fstream>
#include <iostream>
// for ostream_iterator
#include <iterator>
#include <set>

using namespace std;

#include "global.h"
#include "boundeddim.h"
#include "interval.h"
#include "params.h"
#include "point.h"
#include "read.h"

// since i moved findProjectiveCluster() to MyCluster::find(),
// findProjectiveCluster() need not be a friend of
// MyClusterParams. therefore, i can move this declaration to its
// rightful place (projcluster.h) from params.h.
enum MyProjectiveClusterReturnType { ZERO_POINTS, ZERO_SEEDS,
                                     NO_RANDOM_SEEDS, NO_RANDOM_INDICES,
                                     NO_QUALITY_CLUSTERS, NO_SIGNIFICANT_CLUSTERS,
                                     CLUSTER_FOUND };


// struct storing the key params i use to control
// MyCluster::find(). these params may also be stored in
// MyClusterParams but i want to store the key subset in this
// structure and to store a vector of this subset sorted by round to
// keep a history of the parameter values.
struct MyKeyClusterParams
{
  int  numSeeds;
  int  sizeDisc;
  int  numDiscs;
  int  whichRound;
  MyNT minimumHomogeneity;
  MyProjectiveClusterReturnType found;
};



class MyPointPtr;


// what i assume is the default point set when there are multiple.
const unsigned int MY_DEFAULT_POINT_SET = 0;





class MyCluster;

// class for box with some bounded dimensions and some unbounded dimensions.
class MyBox
{
  friend class MyCluster;
public:
  MyBox()
      : numBounded(0), intervals(), realIntervals()
    {}
  

  // construct from list of bounded dimensions and bounds.
  MyBox(const vector< MyNT > &lows, const vector< MyNT > &highs,
        const vector< bool > &bounded);
  MyBox(const MyPoint &centre, const vector< MyNT > &width, const vector< bool > &bounded);

  MyBox(const MyBox& copy)
      : numBounded(copy.numBounded),
        intervals(copy.intervals), realIntervals(copy.realIntervals)
    {}
  const MyBox& operator=(const MyBox& rhs)
    {
      if (this != &rhs)
        {
          numBounded = rhs.numBounded;
          intervals.clear();
          intervals = rhs.intervals;
          realIntervals.clear();
          realIntervals = rhs.realIntervals;
        }
      return(*this);
    }
  // see p 62, item 14 of effective c++ for why i shouldn't declare a
  // virtual destructor in a non-base class.
//  
//  virtual
  ~MyBox()
    {}
  // free up memory used by MyBox.
//  virtual
  void free()
    {
      intervals.clear();
      realIntervals.clear();
    }
  
  bool operator==(const MyBox& rhs) const;
  
  unsigned int getNumBounded() const
    {
      return(numBounded);
    }
  
  
  // check if point lies inside the box.
  bool contains(const MyPoint &point) const
    {
      for (unsigned int i = 0; i < numBounded; i++)
//        if (!::contains(lows[i], highs[i], point[boundedDims[i]]))

        // should not check against realLows and realHighs since they
        // may not be defined yet! i might be calling this function
        // from MyCluster::findContainedPoints(). 
        if (!intervals[i].contains(point[intervals[i].getDimension()]))
          return(false);
      return(true);
    }

  // check whether dim is bounded in the cluster.
  bool checkBounded(int dim) const
    {
      return(realIntervals.end() != _findDimension(dim));
    }
  



  // checks if box is significant. if there is any bounded dimension
  // for which the box's side is not equal to any of the intervals for
  // that bounded dimension, returns false. otherwise, returns true. 
  // also see the comments before the implementatiom of this function
  // in projcluster.C
  bool isSignificant(vector< vector< MyInterval > >& intervals);

  // checks what fraction of the volume (should it be perimeter?) of box is contained within significant intervals. 
  MyNT isSoftContained(vector< vector< MyInterval > >& intervals);

  // return true iff every bounded interval has significance greater than maxPValue.  
  bool isSignificant(const MyPointSet& points, MyNT maxPValue);
  
  void getPValues(const MyPointSet& points, vector< MyNT >& pvalues)
    {
      for (unsigned int i = 0; i < numBounded; i++)
        pvalues.push_back(computePValue(points, i));
    }
      
  
  // set realLows and realHighs to be the bounding box of points.
  void update(const MyPointSet& points, const vector< int >& indices);
  // update realLows and realHighs to reflect addition of points[index].
  void update(const MyPointSet& points, int index);

#ifdef DEBUG  
  // don't initialise by default to NULL since the two versions of
  // MyBox::check() become ambiguous.
  bool check(const vector< MyNT >* lows,
             const vector< MyNT >* highs) const;
  bool check(const vector< MyNT >* width) const;
#endif
  
  void print(ostream& ostr, const vector< MyDimInfo >& dimensionInfo,
             MyVerbosityLevel verbosity) const;

private:
  void construct(const vector< MyNT > &lows,
                 const vector< MyNT > &highs, const vector< bool > &bounded);
  void construct(const MyPoint &centre, const vector< MyNT > &width,
                 const vector< bool > &bounded);

  // compute the p-value of the ith bounded interval.
  MyNT computePValue(const MyPointSet& points, int i);


  
  // find the entry in realIntervals that whose dimension == dim.
  vector< MyDimensionInterval >::const_iterator _findDimension(unsigned int dim) const
    {
      vector< MyDimensionInterval >::const_iterator itr;
      // find a MyDimensionInterval whose dimension == dim.
      for (itr = realIntervals.begin(); itr != realIntervals.end(); itr++)
        if (dim == (*itr).getDimension())
          break;
      return(itr);
    }
  
  // find the entry in realIntervals that whose dimension == dim.
  vector< MyDimensionInterval >::iterator _findDimension(unsigned int dim)
    {
      vector< MyDimensionInterval >::iterator itr;
      // find a MyDimensionInterval whose dimension == dim.
      for (itr = realIntervals.begin(); itr != realIntervals.end(); itr++)
        if (dim == (*itr).getDimension())
          break;
      return(itr);
    }
  
private:
  unsigned int numBounded;
  // intervals are the seed -/+ per-dim interval that the point is
  // associated with.
  vector< MyDimensionInterval > intervals;
  // realIntervals are the actual bounding box of the
  // contained points.
  vector< MyDimensionInterval > realIntervals;
};

#ifdef USE_CLUSTER_DETAILS
// detail of how the cluster is created.
struct MyClusterDetails
{
public:
  
  // how should i store a pointer to the seed? 
  //
  // (i) it is not possible to store a const MyPoint& because when i
  // construct a MyCluster, the seed i have could be allocated
  // temporarily. the seed might get overwritten or go out of scope
  // making the reference inside MyClusterDetails invalid. another
  // issue is that the seed as a const MyPoint& means that i have to
  // store a *pointer* to MyClusterDetails in MyCluster. i hate this
  // fact and i dont know if i can do anything better. the reason is
  // as follows: i want to create arrays of MyClusters (e.g., the
  // allClusters variable in cluster.C, i have to be able to call the
  // no-parameter MyCluster constructor. then when i assign a real
  // MyCluster to an "empty" MyCluster, i have to call the
  // MyClusterDetails assignment operator, which does not allow me to
  // assign to the const reference to MyPoint in MyClusterDetails. by
  // sotring a pointer to MyClusterDetails in MyCluster, i can use new
  // to assign to MyClusterDetails when i call the MyCluster
  // assignment operator.
  //
  // (ii) i could store an iterator that points to the seed inside
  // MyPointSet. however the iterator can become invalid if the vector
  // containing the points inside MyPointSet is reallocated.
  //
  // (iii) i store an index to the point. this index is invalid only
  // if i delete or insert a point with smaller index, which i believe
  // will not happen. of course, storing the index means that i have
  // to change the MyCluster constructor to accept an index.

  typedef int MyPointPtr;

  MyClusterDetails()
      : seed(), discIndices()
    {}
  MyClusterDetails(MyPointPtr s, const vector< int >& ind)
      : seed(s), discIndices(ind)
    {}
  MyClusterDetails(const MyClusterDetails& copy)
      : seed(copy.seed), discIndices(copy.discIndices)
    {}
  ~MyClusterDetails()
    {}
  // not providing an implementation.
  const MyClusterDetails& operator=(const MyClusterDetails& rhs);
  
//     const MyClusterDetails& operator=(const MyClusterDetails& rhs)
//     {
//       if (this != &rhs)
//         {
//           seed = rhs.seed;
//           discIndices.clear();
//           discIndices = rhs.discIndices;
//         }
//       return(*this);
//     }
  
  const MyPoint& getSeed(const MyPointSet& points) const
    {
      return(points[seed]);
    }
  
//   const MyPointPtr& getSeed(const MyPointSet& points) const
//     {
//       return(seed);
//     }
  
public:
  MyPointPtr seed;
//  int index;
  vector< int > discIndices;
};
#endif // USE_CLUSTER_DETAILS


// forward declarations.
// struct MyClusterLess;
// struct MyClusterMore;
class MyClusterSet;


// stores box and points (or pointers to them) that lie inside the box. 
class MyCluster
{
public:
#ifdef USE_CLUSTER_DETAILS
  typedef MyClusterDetails::MyPointPtr MyPointPtr;
#endif  
  
  MyCluster()
      : box(), numContained(0), containedPoints(), realContainedPoints(),
#ifdef USE_CLUSTER_DETAILS
        details(NULL), 
#endif     
           filledIn(false), pointsFlipped(false)
    {}

#ifdef USE_BOUNDED_DIM_TYPE
  // construct from seed and MyBoundedDim.
  MyCluster(
#ifdef USE_CLUSTER_DETAILS                
                  const MyPointPtr seedPtr, 
#endif // USE_CLUSTER_DETAILS
                  const MyBoundedDim &boundedDimInfo,
            const MyPointSet& points, const vector< int >& indices,
            const MyClusterParams &params);
#else

  // construct from list of bounded dimensions and bounds.
  MyCluster(
#ifdef USE_CLUSTER_DETAILS
                  const MyPointPtr seedPtr, 
#endif // USE_CLUSTER_DETAILS             
                  const vector< MyNT > &lows,
            const vector< MyNT > &highs, const MyPointSet& points,
            const vector< int > &indices,
            const MyClusterParams &params);
//            bool saveMemory = false);
  MyCluster(
#ifdef USE_CLUSTER_DETAILS                
                  const MyPointPtr seedPtr, 
#endif // USE_CLUSTER_DETAILS             
                  const vector< MyNT > &width,
            const MyPointSet& points,
            const vector< int > &indices,
            const MyClusterParams &params);
//            bool saveMemory = false);

  // SEEDLESS version.
  MyCluster(const vector< MyNT > &lows, const vector< MyNT > &highs,
            const MyPointSet& points, bool saveMemory = false);

  // SEEDLESS, LOWLESS, HIGHLESS version
  MyCluster(const MyPointSet& points,
            const vector< int > &indices,
            const MyClusterParams &params);
#endif
  
  MyCluster(const MyCluster& copy)
      : box(copy.box), numContained(copy.numContained),
        containedPoints(copy.containedPoints),
        realContainedPoints(copy.realContainedPoints),
#ifdef USE_CLUSTER_DETAILS
        details(copy.details), 
#endif
           filledIn(copy.filledIn),

        pointsFlipped(copy.pointsFlipped)
    {}
  
  const MyCluster& operator=(const MyCluster& rhs)
    {
      if (this != &rhs)
        {
          box = rhs.box;
          numContained = rhs.numContained;
          containedPoints.clear();
          containedPoints = rhs.containedPoints;
          realContainedPoints.clear();
          realContainedPoints = rhs.realContainedPoints;
          // don't delete details because i dont allocate details
          // using new when use operator=. if i delete, then i can
          // delete a pointer i actually need (e.g., when i assign
          // inside MY_forall).
//          delete details;
#ifdef USE_CLUSTER_DETAILS
          details = rhs.details;
#endif
          filledIn = rhs.filledIn;
          pointsFlipped = rhs.pointsFlipped;
        }
      return(*this);
    }
//  virtual
  ~MyCluster()
    {
//      delete details;
    }

  bool operator==(const MyCluster& rhs) const
    {
      if (this == &rhs)
        return(true);
      // enough to compare the boxes since the boxes store the real bounding box.
      return(box == rhs.box);
    }

    bool operator!=(const MyCluster& rhs) const
    {
      return(!(*this == rhs));
    }

  // is *this < other.
  bool MyCluster::operator<(const MyCluster &other) const
    {
      // i could use MyCluster::measure() but since it can depend on
      // shiftValue (which is yucky in its own way), i am not.
      return(getNumBounded()*getNumContained() <
             other.getNumBounded()*other.getNumContained());
    }
  
  // is *this > other.
  bool MyCluster::operator>(const MyCluster &other) const
    {
      return(!(*this < other));
    }

#ifdef USE_BOUNDED_DIM_TYPE
  void construct(const MyBoundedDim &boundedDimInfo,
                 const MyPointSet& points, const vector< int >& indices,
                 const MyClusterParams &params, bool brandNew = false);
#else  
  void construct(const vector< MyNT > &width, const MyPointSet& points,
                 bool brandNew = false);

  void construct(const MyPointSet& points,
                 const vector< int > &discIndices,
                 const MyClusterParams &params,
                 bool brandNew = false);
#endif // USE_BOUNDED_DIM_TYPE
  
  // i need this for find_greedy().
  void construct(const vector< MyNT > &lows, const vector< MyNT > &highs,
                 const MyPointSet& points, bool brandNew = false);

  // used to be findProjectiveCluster() in cluster.C. find a
  // projective cluster and set the found cluster to be *this.
  MyProjectiveClusterReturnType
  find(const MyPointSet& points, MyClusterParams& params, ofstream &ostr,
       const MyBoundedDim &boundedDimInfo,
//        const vector< MyNT >& width,
//        vector< vector< MyInterval > >* intervals = NULL,
//        vector< vector< vector< int > > >* intervalsForPoints = NULL,
       // passing MyClusterMore since i want clusters to be sorted
       // in decreasing order of measure.
       MyClusterSet *bestClusters = NULL
       )
  {

    cerr << "This version of MyCluster::find() is not implemented!" << endl;
    exit(-1);
    // call the other find() with vectors.
//     return(find(vector< MyPointSet >(1, points), params, ostr, 
//                       vector< MyBoundedDim >(1, boundedDimInfo), bestClusters));
  }

  MyProjectiveClusterReturnType
  find(vector< MyPointSet > &pointSets, 
          MyClusterParams& params, ofstream &ostr,
          //       vector< MyBoundedDim > &boundedDimInfo,
       vector< const MyBoundedDim* > &boundedDimInfo,
       MyClusterSet *bestClusters = NULL
       );

  // find a cluster using a greedy algorithm.
  MyProjectiveClusterReturnType
  find_greedy(const MyPointSet& points, MyClusterParams& params,
              ofstream &ostr,
              const vector< MyNT >& widths,
              vector< vector< MyInterval > >* intervals,
              vector< vector< vector< int > > >* intervalsForPoints);
  

  // return true if *this and other are similar, whatever that might
  // mean.
  bool MyCluster::similar(const MyCluster &other) const;
  // this functions uses an argument for the threshold for similarity.
  bool MyCluster::similar(const MyCluster &other, MyNT threshold) const;
  
  // compare this and rhs. currently, i just compute the size of the
  // symmetric difference between the two sets of contained points (in
  // numDiff) and the number of points that are in rhs but not in this
  // (in numNew).
  void compare(const MyCluster& rhs, int& numDiff, int& numNew) const;
  bool compare(vector< MyCluster >& clusters, MyVerbosityLevel verbose);
  void merge(const MyCluster& cluster);

  // compute a Shamir-like p-value based on the #points from each
  // class in the cluster.
  MyNT pValueShamir(const MyPointSet& points) const;
  // compute the p-value of the cluster based on the approximation to
  // the gaussian PDF using the central limit theorem.
  MyNT pValueGaussian(const MyPointSet &points) const;
  
  bool quality(MyNT alpha, int numPoints, int numDims) const;
  bool support(MyNT alpha, int numPoints, int numDims) const;
  MyNT homogeneity(const MyPointSet& points) const;
  // compute the type of this cluster. there is no point calling this
  // function if the class names are not provided in the command-line.
  string type(const MyPointSet &points) const;
  
  void getContainedPoints(vector< int > &cont) const
  {
    cont = containedPoints;
  }
  unsigned int getNumContained() const
    {
      return(numContained);
    }

  unsigned int getNumBounded() const
    {
      return(box.getNumBounded());
    }

  // count the number of points in cluster that are contained in this cluster.
  int countNumContained(const MyPointSet& points, const MyCluster& cluster)
    const;
  
  // count the number of points (referred to by indices) that are
  // contained in this cluster.
  int countNumContained(const MyPointSet& points, const vector< int >& indices)
    const;

  // compute the set of points contained in cluster that are contained
  // in this cluster. store the resulting set of indices in
  // pointsInside and return the number of such points.
  int findContainedPoints(const MyPointSet& points, const MyCluster& cluster,
                          vector< int >& pointsInside);
  
  // compute the set of points referred to by indices that are
  // contained in this cluster. store the resulting set of indices in
  // pointsInside and return the number of such points.
  int findContainedPoints(const MyPointSet& points,
                          const vector< int >& indices,
                          vector< int >& pointsInside) const;

  // compute the points in *this not contained in any cluster in clusters.
  int findNotContainedPoints(const MyPointSet& points,
                             const vector< MyCluster> &clusters,
                             vector< int >& pointsOutside) const;
  
  // compute the set of points contained in *this that are not
  // contained in cluster. store the resulting set of indices in
  // pointsOutside and return the number of such points.
  int findNotContainedPoints(const MyPointSet& points, 
                             const MyCluster& cluster,
                             vector< int >& pointsOutside) const;
  
  
  // compute the set of points referred to by indices that are not
  // contained in this cluster. store the resulting set of indices in
  // pointsOutside and return the number of such points.
  int findNotContainedPoints(const MyPointSet& points,
                             const vector< int >& indices,
                             vector< int >& pointsOutside) const;
  

  // compute the measure of this cluster.
  //
  // GCC3: COMMENTING OUT because gcc3 complains about ambiguous calls
  // to log10 and pow. i don't need this function anyway.
  //
//   MyNT measure(int shiftValue) const
//     {
//       // return a number equal to the concatenation of
//       // getNumContained() and getNumBounded().
//       int base = static_cast< int > (1 + log10(shiftValue));
//       return(getNumBounded()*pow(10, base) + getNumContained());
//     }

  // compute the measure of this cluster.
  MyNT measure() const
    {
      // return a number equal to the product of
      // getNumContained() and getNumBounded().
      return(getNumBounded()*getNumContained());
    }

  // see comments before definition of MyPoint::chosen to see
  // the difference between hiding and choosing a point.
  //
  // called in main() in cluster.C after a cluster is selected (when
  // CLUSTER_FOUND == found) but only when params.getPartition() is
  // true, i.e., when points are being partitioned between the
  // clusters. 
  void hideContainedPoints(MyPointSet& points) const
    {
      points.hidePoints(containedPoints);
    }

  // called in main() in cluster.C after a cluster is selected (when
  // CLUSTER_FOUND == found).
  void chooseContainedPoints(MyPointSet& points) const
    {
      points.choosePoints(containedPoints);
    }

  // hide just the cluster. i do so by replacing the coords
  // corresponding to the points and dims in the cluster by random
  // numbers.
  void hide(MyPointSet& points) const
    {
      points.hide(containedPoints, box.intervals);
    }




  // add 1 to pointCounts[i] if point i is contained in *this
  void updatePointCounts(vector< unsigned int > &pointCounts) const
    {
      const vector< int > &contPoints = realContainedPoints;
//    const vector< int > &contPoints = containedPoints;
      
      for (unsigned int i = 0; i < contPoints.size(); i++)
        pointCounts[contPoints[i]]++;
    }
  
        
  void print(ostream& ostr, MyVerbosityLevel verbosity) const;
  // see comments for definition of print() in cluster.C about whether
  // print() can be const or not.
  void print(ostream& ostr, const MyPointSet& points,
             MyVerbosityLevel verbosity, bool printItemSetFormat = false) const;
  // print info on p-values, type, quality, etc. some of them depend
  // on the existence of class names.
  void printExtra(ostream& ostr, const MyPointSet& points) const;


  // functions for adding/deleting poitns and dimensions.
  bool addPoint(const MyPointSet& points, int index);
  bool deletePoint(const MyPointSet& points, int index);
  bool addDimension(int dim, MyNT start, MyNT end, const MyPointSet& points,
                    bool check = false);
  bool addDimension(int dim, MyNT start, MyNT end, const vector< int >& contPoints);
  bool addDimension(const MyDimensionInterval& interval,
                    const vector< int >& contPoints);

  
  // TODO: DESIGN. do i need this function?  
  bool addDimension(int dim, MyNT width, const MyPointSet& points);
  bool deleteDimension(int dim, const MyPointSet& points);
  
  void anneal(const MyPointSet& points, const vector< MyNT >& widths,
                       MyClusterParams& params);
  
  bool isSignificant(vector< vector< MyInterval > >& intervals)
    {
      return(box.isSignificant(intervals));
    }
  
  MyNT isSoftContained(vector< vector< MyInterval > >& intervals)
    {
      return(box.isSoftContained(intervals));
    }

  bool isSignificant(const MyPointSet& points, MyNT maxPValue)
    {
      return(box.isSignificant(points, maxPValue));
    }
  
  void getPValues(const MyPointSet& points, vector< MyNT >& pvalues)
    {
      return(box.getPValues(points, pvalues));
    }
      
  // returns the number of containedPoints contained in [start, end]
  // along dimension dim.
  int countNumContainedForDim(int dim, MyNT start, MyNT end, const MyPointSet &points) const
    {
      int numIn = 0;
      for (unsigned int i = 0; i < numContained; i++)
        numIn += ::contains(start, end, points[containedPoints[i]][dim]);
      return(numIn);
    }

private:

  // functions that return the real index of a given sample or a given
  // gene. these functions are the only ones that actually refer to
  // gene expression data in the function name. in all of them, the
  // flipped argument indicates whether the point set has been flipped
  // or not. if not, the samples correspond to the bounded dimensions,
  // otherwise to the genes.
  unsigned int getNumSamples(bool flipped) const
    {
      if (flipped)
        return(getNumContained());
      else
        return(getNumBounded());
    }
  
  unsigned int getNumGenes(bool flipped) const
    {
      if (flipped)
        return(getNumBounded());
      else
        return(getNumContained());
    }
  
  unsigned int getSampleIndex(int i, bool flipped) const
    {
      if (flipped)
        return(containedPoints[i]);
      else
        return(box.intervals[i].getDimension());
    }
  unsigned  int getGeneIndex(int i, bool flipped) const
    {
      if (flipped)
        return(box.intervals[i].getDimension());
      else
        return(containedPoints[i]);
    }
  

  
  // free up memory used by MyCluster
//  virtual
  void free()
    {
      containedPoints.clear();
      box.free();
      filledIn = false;
    }

#ifdef DEBUG
  bool check(const MyPointSet& points) const;
#endif
  
  // find the largest class of points in the cluster.
  void _findLargestClass(const MyPointSet& points,
                         unsigned int& maxNum,
                         string& maxClass, bool real = true) const;

  // take list of points and store those that lie inside. if
  // partitionPoints is true, then i only include points in this
  // cluster that have not already been "chosen" (put in some other
  // cluster).
  void findContainedPoints(const MyPointSet &points)
    {
      for (unsigned int i = 0; i < points.size(); i++)
        // only consider points that have not already been hidden.
        // whether the hidden bit is set or not depends on the -p
        // command line parameter (params.getPartition()). a point's
        // hidden bit is set in MyCluster::hideContainedPoints().
        if (box.contains(points[i]))
          {
            realContainedPoints.push_back(i);
            if (!points[i].getHidden())
              containedPoints.push_back(i);
          }
      
    }
  
  void findContainedPoints(const MyPointSet &points, vector< int >& contPoints) const
    {
      for (unsigned int i = 0; i < points.size(); i++)
        // only consider points that have not already been
        // chosen. whether the chosen bit is set or not depend on the
        // command line parameter. a point's chosen bit is set in the
        // main function that calls findProjectiveCluster().
        if (!points[i].getHidden() && (box.contains(points[i])))
          contPoints.push_back(i);
    }

  // returns a vector of points in containedPoints that lie inside
  // [start, end] along dimension dim.
  inline int findContainedPointsForDim(int dim, MyNT start, MyNT end,
                                       const MyPointSet &points, vector< int >& newContained)
    {
      for (unsigned int i = 0; i < numContained; i++)
        if (::contains(start, end, points[containedPoints[i]][dim]))
          newContained.push_back(containedPoints[i]);
      return(newContained.size());
    }

  int findMaxNumContainedForDim(int dim, const vector< MyInterval >* intervals,
                                MyNT width, const MyPointSet& points,
                                MyNT& maxStart, MyNT& maxEnd);
  



  void mergeContainedPoints(const vector< int >& pts);

  void computeBoundedDimensions(const MyPointSet& points,
                                const MyPoint& seed, const vector< int >&discIndices,
                                const vector< MyNT > &lows, const vector< MyNT > &highs,
                                vector< bool >& bounded) const;
  void computeBoundedDimensions(const MyPointSet& points,
                                const MyPoint& seed, const vector< int >&discIndices,
                                const vector< MyNT > &width,
                                vector< bool >& bounded) const;

  void improveDimensions(const MyPointSet& points, MyClusterParams& params,
                         vector< vector< MyInterval > >* intervals);
  

private:
  MyBox box;
  unsigned int numContained;
  // the difference between these two vectors is that containedPoints
  // has points contained in box but not hidden. realContainedPoints
  // does not care when the points are hidden or not.
  vector< int > containedPoints;
  vector< int > realContainedPoints;

#ifdef USE_CLUSTER_DETAILS
  MyClusterDetails* details;
#endif
  // are all fields filled in?
  bool filledIn;

  // were the points flipped before constructing the cluster. storing
  // this bit with each cluster is a waste of space but the advantage
  // is that i can set the bit when constructing the cluster and not
  // have to worry about passing the bit to all the functions that
  // need it.
  bool pointsFlipped;
};


// define after defining MyCluster.
struct MyClusterLess :
  public binary_function< const MyCluster&, const MyCluster&, bool >
{
  bool operator()(const MyCluster &lhs, const MyCluster &rhs) const
  {
    return(lhs < rhs);
  }
};

struct MyClusterMore :
  public binary_function< const MyCluster&, const MyCluster&, bool >
{
  bool operator()(const MyCluster &lhs, const MyCluster &rhs) const
  {
    return(lhs > rhs);
  }
};



// a set of sorted clusters.
struct MyClusterSet
{
public:
  MyClusterSet(int max = 0)
    // i hope it is ok to initialise minMeasure as 0.
      :  clusters(), maxNum(max), minMeasure(0), minPtr()
    {}

  ~MyClusterSet()
    {}
  
  int size() const
    {
      return(clusters.size());
    }
  int maxSize() const
    {
      return(maxNum);
    }
  bool empty() const
    {
      return(clusters.empty());
    }
  // is the set full? if so, i cannot add any more clusters.
  bool full() const
    {
      return(maxNum == clusters.size());
    }
  bool insert(const MyCluster &cluster, const MyPointSet &points,
              const MyClusterParams &params);
  void print(ostream &ostr) const;
  void printAll(ostream &ostr, const MyPointSet &points,
                 MyClusterParams &params) const;
  void printAllTypes(ostream &ostr, const MyPointSet &points) const;
  
  
    
private:
  
private:
  // using MyClusterMore to store clusters in decreasing order of
  // measure.
  typedef set< MyCluster, MyClusterMore > _MyClusterSet;

  set< MyCluster, MyClusterMore > clusters;
  // the maximum number i want to store. if -1, i can store as many as
  // i want.
  int maxNum;
  // the smallest measure stored in clusters.
  MyNT minMeasure;
  // a pointer to the cluster with measure minMeasure.
  _MyClusterSet::iterator minPtr;
  
};



#endif // _PROJCLUSTER_H