point.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 * Copyright (c) 2010-2011 Christopher L. Poirel                          *
 *                                                                        *
 * 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/>.      *
 *                                                                        *
 **************************************************************************/

// Purpose: defines point class.

#ifndef POINT_H
#define POINT_H

#include <map>
#include <set>

#include "constants.h"
// Test for GCC >= 4.3. from http://gcc.gnu.org/onlinedocs/cpp/Common-Predefined-Macros.html. See definition of GCC_VERSION macro in libutil/constants.h
#if (GCC_VERSION >= 40300)
// see http://gcc.gnu.org/gcc-4.3/changes.html
#include <tr1/unordered_set>
using namespace std::tr1;
#else
#include <ext/hash_set>
using namespace __gnu_cxx;
#endif

#include "myiterator.h"

#include "format.h"
#include "global.h"
#include "histogram.h"
#include "info.h"
#include "interval.h"
#include "params.h"
#include "random.h"

// forward declarations to avoid including params.h
// class MyClusterParams;
// enum MyVerbosityLevel;


// do not flip rows and columns in data. lots of code has to change as
// a result so i am guarding the new code with #define.
//#define DO_NOT_FLIP_DATA 1

// a base class that stores an n-dimensional point.
template <class NT>
class Point_d 
{
public:
  Point_d(int d=0)
    : dimension(0), missingDims(), coords()
    {}
  Point_d (int d, const vector< NT> &c, 
                 const set< unsigned int > *missing = NULL)
    : dimension(d), missingDims(), coords(c)
    {
         if (missing)
           missingDims = *missing;
    }
  Point_d(const Point_d& other)
      : dimension(other.dimension), missingDims(other.missingDims), coords(other.coords)
    {}
  
  
  virtual ~Point_d()
    {}
  
  unsigned int getNumDimensions() const
    {
      return(dimension);
    }
  
  NT operator[](int i) const
    {
      // TODO: BUG. check if i is valid.
      return(coords[i]);
    }

  // is the coord for dim i missing in this point?
  bool isDimMissing(unsigned int i ) const
  {
    // HACKHACKHACK. i get a segfault here, probably only when
    // missingDims has no entries. i am adding a check for the size,
    // which is quite an ugly HACKHACKHACK.
    return((0 != missingDims.size()) && (missingDims.end() != missingDims.find(i)));
//    return(missingDims.end() != missingDims.find(i));
  }

  unsigned int getNumMissingDimensions() const
    {
      return(missingDims.size());
    }
  
  void standardise()
    {
      MyNT mean = computeMean(coords, &missingDims);
      MyNT variance = computeVariance(coords, mean, &missingDims);
      standardise(mean, sqrt(variance));
    }
  

  void standardise(MyNT m, MyNT s)
    {
      for (unsigned int i = 0; i < coords.size(); i++)
        if (!isDimMissing(i))
          coords[i] = (coords[i] - m)/s;
    }
  
  

  void print(ostream& ostr) const
    {
      for (typename vector< NT>::const_iterator itr = coords.begin(); itr != coords.end(); itr++)
        ostr << "\t" << *itr;
      ostr << endl;
    }

  friend ostream& operator<<(ostream& ostr,  const Point_d &pt)
    {
      pt.print(ostr);
      return(ostr);
    }
  
  friend istream& operator>>(istream& istr, Point_d &pt)
    {
      cerr << "ERROR! istream& operator>>(istream& istr, Point_d &pt) not implemented. Exiting.\n";
      exit(-1);
      return(istr);
    }
  
  virtual void randomise()
    {
      random_shuffle(coords.begin(), coords.end());
    }
  
//protected:
    // TODO: DESIGN. should i allow user to modify any coords?
  NT& operator[](int i)
    {
      return(coords[i]);
    }
  
protected:

  void addMissingDim(unsigned int i)
  {
    missingDims.insert(i);
  }

  void deleteMissingDim(unsigned int i)
  {
    missingDims.erase(i);
  }


protected:
  unsigned int dimension;
  // list of coords with missing values.
  set< unsigned int > missingDims;
  vector< NT > coords;
};

typedef Point_d< MyNT > _MyPoint;

class MyPointInfo : public MyInfo
{
public:
  MyPointInfo(const MyString &id = "", const MyString &sid = "", const string &cname = "") 
      : MyInfo(id, sid), className(cname)
    {}
  MyPointInfo(const MyPointInfo& info)
      : MyInfo(info), className(info.className)
    {}
  MyPointInfo(const MyInfo& info)
      : MyInfo(info), className()
    {}
  virtual ~MyPointInfo()
    {}

  // return a clone of this. this function is a way of providing a
  // virtual constructor for the MyInfo class hierarchy. 
  virtual MyPointInfo* clone() const // = 0;
    {
      MyPointInfo* clone = new MyPointInfo(*this);
      clone->setClassName(className);
      return(clone);
    }

  virtual void setClassName(string name)
    {
      className = name;
    }
  virtual string getClassName() const
    {
      return(className);
    }
  
  virtual void print(ostream& ostr) const
    {
      MyInfo::print(ostr);
      ostr << "Class name = " << className << endl;
    }
  
private:
  // the name of the class the point belongs to.
  string className;
};

// define MyDimInfo to be exactly the same as MyPointInfo. this class
// is useful when i don't flip a point set.
typedef MyPointInfo MyDimInfo;
// MyDimInfo is a plain old string.
// typedef string MyDimInfo;

// the real point class with bells and whistles (e.g., info).
class MyPoint : public _MyPoint
{
  friend class MyPointSet;
  
public:
  MyPoint()
      : _MyPoint(), info(NULL), chosen(false)
    {}
  // need constructors that use both MyInfo (for readCOGData in read.C) and MyPointInfo
  MyPoint(int dim, const vector< MyNT >& coords, MyInfo *inf)
      : _MyPoint(dim, coords), info(NULL), chosen(false), hidden(false)
    {
      if (inf)
        setInfo(inf);
    }

  MyPoint(int dim, const vector< MyNT >& coords, 
                //MyPointInfo *inf, p
                MyInfo *inf, 
                set< unsigned int > *missing = NULL)
    : _MyPoint(dim, coords, missing), info(NULL), chosen(false), hidden(false)
    {
      if (inf)
        setInfo(inf);
    }

  MyPoint(const _MyPoint &pt, MyPointInfo *inf)
      : _MyPoint(pt), info(NULL), chosen(false), hidden(false)
    {
      if (inf)
        setInfo(inf);
    }
  MyPoint(const MyPoint& copy)
      : _MyPoint(copy), info(NULL), chosen(copy.chosen), hidden(copy.hidden)
    {
      if (copy.info)
        setInfo(copy.info);
    }
  const MyPoint& operator=(const MyPoint& rhs)
    {
      if (this != &rhs)
        {
          this->_MyPoint::operator=(rhs);
          if (rhs.info)
            setInfo(rhs.info);
          chosen = rhs.chosen;
          hidden = rhs.hidden;
        }
      return(*this);
    }
  
  virtual ~MyPoint()
    {
      delete info;
    }
  
  const MyInfo* getInfo() const
    {
      return(info);
    }

  MyString getName() const
    {
      if (info)
        {
          // return short name if it exists.
          if ("" != info->getShortName())
            return(info->getShortName());
          else
            return(info->getName());
        }
      else
        return(string("NONAME"));
    }

  // TODO: DESIGN. i have to set info's members when i read
  // MyAffyGene. is there a way to still keep info a const.
  //
  void setInfo(MyPointInfo* inf)
    {
#if 0      
      info = inf;
#endif
      delete info;
      
//       if (!info)
        info = inf->clone();
      // TODO: LEAK. is there a leak here? what happens to the old
      // values of name and shortName?
//       info->setName(inf->getName());
//       info->setShortName(inf->getShortName());
    }
  void setInfo(MyInfo* inf)
    {
      delete info;
      // TODO: DESIGN. why do i need to clone when a new works?
      info = new MyPointInfo(*inf);
    }

#if 0      
  void setInfo(const MyString &name, const MyString &shortName,
               const string className = "")
    {
      if (!info)
        info = new MyPointInfo;
      // TODO: LEAK. is there a leak here? what happens to the old
      // values of name and shortName?
      info->setName(name);
      info->setShortName(shortName);
      info->setClassName(className);
    }
#endif

  string getClassName() const
    {
      if (info)
        return(info->getClassName());
      else
        return("NONAME");
//        cerr << "Point " << getName() << " has no associated info to store class name with.\n";
    }
  
  void setClassName(string name)
    {
      if (info)
        info->setClassName(name);
      else
        cerr << "Point " << getName() << " has no associated info to store class name with.\n";
    }
  bool getChosen() const
    {
      return(chosen);
    }

  bool getHidden() const
    {
      return(hidden);
    }

  MyNT getMaxCoord() const
    {
      if (getNumMissingDimensions() == getNumDimensions())
        return(0);
      unsigned int i = 0;
      MyNT max;
      while (isDimMissing(i))
        i++;
      // coords[i] is not missing.
      max = coords[i];
      while (i < getNumDimensions())
        {
          if (!isDimMissing(i) && (coords[i] > max))
            max = coords[i];
          i++;
        }
      return(max);
    }
  
  MyNT getMinCoord() const
    {
      if (getNumMissingDimensions() == getNumDimensions())
        return(0);
      unsigned int i = 0;
      MyNT min;
      while (isDimMissing(i))
        i++;
      // coords[i] is not missing.
      min = coords[i];
      while (i < getNumDimensions())
        {
          if (!isDimMissing(i) && (coords[i] < min))
            min = coords[i];
          i++;
        }
      return(min);
    }

  MyNT computeNorm(unsigned int p) const;
  
  MyNT computeSparsityHoyer() const;
  

  
  int computeDiscreteCorrelation(const MyPoint &other,
                                 bool computeSign = false) const;

  MyNT computePearsonsCorrelation(const MyPoint &other,
                                  const MyNT *mean1 = NULL, const MyNT *stdDev1 = NULL,
                                  const MyNT *mean2 = NULL, const MyNT *stdDev2 = NULL)
    const;


  void computeRanks(MyPoint &rankedPoint) const;
  
private:

  // only friends and derived classes can set chosen and hidden bits.
  //
  // called by MyPointSet::choosePoints().
  void setChosen(bool c)
    {
      chosen = c;
    }
  
  // called by MyPointSet::hidePoints().
  void setHidden(bool h)
    {
      hidden = h;
    }


  // distance related functions.
  MyNT distance(const MyPoint &other)
    {
      MyNT distance = 0;
      for (unsigned int i = 0; i < getNumDimensions(); i++)
        // no need to ignore missing dimensions since they contribute 0.
        distance += coords[i]*other.coords[i];
      return(sqrt(distance));
    }

  MyNT mean() const
    {
      MyNT _mean = 0;
      for (unsigned int i = 0; i < getNumDimensions(); i++)
        {
          // ignore missing dimensions.
          if (isDimMissing(i))
            continue;
          _mean += coords[i];
        }
      _mean /= (getNumDimensions() - getNumMissingDimensions());
      return(_mean);
    }
  
  MyNT variance(MyNT mean) const
    {
      MyNT _variance = 0;
      for (unsigned int i = 0; i < getNumDimensions(); i++)
        {
          // ignore missing dimensions.
          if (isDimMissing(i))
            continue;
          _variance += (coords[i] - mean)*(coords[i] - mean);
        }
      _variance /= (getNumDimensions() - getNumMissingDimensions());
      return(_variance);
    }
  
  MyNT standardDeviation(MyNT mean) const
    {
      return(sqrt(variance(mean)));
    }
  
  
private:
  MyPointInfo* info;
  // to keep track of whether a point belongs to any cluster or not. i
  // need two variables for the following reasons: i use chosen just
  // as a statistic (to keep a count of how many points have been
  // picked). i use hidden to "hide" points and thus to control the
  // behaviour of the gropram based on the value of
  // params.getPartition(). if i use only hidden or chosen, then i have
  // to pass params.getPartition() to the MyCluster constructor, which
  // seems totally wrong. by using both chosen and hidden, i can set
  // hidden inside MyCluster::hideContainedPoints() and chosen inside
  // MyCluster::findContainedPoints().
  bool chosen;
  bool hidden;
};


typedef pair< int, MyNT > MyEntropy;

struct compareEntropies : public binary_function< MyEntropy, MyEntropy, bool> 
{
  bool operator()(const MyEntropy& __x, const MyEntropy& __y) const
    {
      return __x.second < __y.second; }
};

// inline ostream& operator<<(ostream& ostr, const MyEntropy& ent)
// {
//   ostr << "(" << ent.first << ", " << ent.second << ")";
//   return(ostr);
// }




typedef map< string, string, ltstr >  MyClassesType;
typedef set< string, ltstr > MyClassNamesType;

const unsigned int NUM_CORRELATION_HISTOGRAM_BINS = 100;
const MyNT MIN_CORRELATION = -1;
const MyNT MAX_CORRELATION = 1;
const MyNT CORRELATION_RANGE = MAX_CORRELATION - MIN_CORRELATION;


class MyPointSet
{
  friend class MyCluster;
  
  typedef map< MyNT, unsigned int > MyDimCountsType;
  typedef MyDimCountsType::const_iterator MyDimCountsConstPtr;
  typedef MyDimCountsType::iterator MyDimCountsPtr;
  
public:

  MyPointSet()
      : numPoints(0), numDimensions(0), numNotChosen(0), numNotHidden(0),
        points(),
        hiddenPoints(), hiddenPointsByClass(), 
           flipped(false),
        namesToIndices(), classNamesExist(false), 
           classNames(),
           classesToIndices(), classCounts(),
        dimInfo(), dimLows(0), dimHighs(0), dimMeans(0), dimSigmas(0), dimCounts(0),
        dimEntropies(), sortedDims(), _pointMeans(), _pointStandardDeviations()
    {}
  MyPointSet(const MyPointSet& copy)
      : numPoints(copy.numPoints), numDimensions(copy.numDimensions),
        numNotChosen(copy.numNotChosen),
        numNotHidden(copy.numNotHidden),
        points(copy.points),
        hiddenPoints(copy.hiddenPoints), 
        hiddenPointsByClass(copy.hiddenPointsByClass), 
           flipped(copy.flipped),
        namesToIndices(copy.namesToIndices),
        classNamesExist(copy.classNamesExist),
           classNames(copy.classNames),
        classesToIndices(copy.classesToIndices),
        classCounts(copy.classCounts),
        dimInfo(copy.dimInfo), dimLows(copy.dimLows), dimHighs(copy.dimHighs),
        dimMeans(copy.dimMeans), dimSigmas(copy.dimSigmas),
        dimCounts(copy.dimCounts),
        dimEntropies(copy.dimEntropies), sortedDims(copy.sortedDims),
        _pointMeans(copy._pointMeans),
        _pointStandardDeviations(copy._pointStandardDeviations)
    {}
  const MyPointSet& operator=(const MyPointSet& rhs)
    {
      if (this != &rhs)
        {
          numPoints = rhs.numPoints;
          numDimensions = rhs.numDimensions;
          numNotChosen = rhs.numNotChosen;
          numNotHidden = rhs.numNotHidden;
          points.clear();
          points = rhs.points;
          hiddenPoints.clear();
          hiddenPoints = rhs.hiddenPoints;
                hiddenPointsByClass.clear();
          hiddenPointsByClass = rhs.hiddenPointsByClass;
          flipped = rhs.flipped,
          namesToIndices.clear();
          namesToIndices = rhs.namesToIndices;
          classNamesExist = rhs.classNamesExist;
                classNames.clear();
                classNames = rhs.classNames;
                classesToIndices.clear();
          classesToIndices = rhs.classesToIndices;
                classCounts.clear();
          classCounts = rhs.classCounts;

          dimInfo.clear();
          dimInfo = rhs.dimInfo;
          dimLows.clear();
          dimLows = rhs.dimLows;
          dimHighs.clear();
          dimHighs = rhs.dimHighs;
          dimMeans.clear();
          dimMeans = rhs.dimMeans;
          dimSigmas.clear();
          dimSigmas = rhs.dimSigmas;
          
          dimCounts.clear();
          dimCounts = rhs.dimCounts;
          
          dimEntropies.clear();
          dimEntropies = rhs.dimEntropies;
          sortedDims.clear();
          sortedDims = rhs.sortedDims;
          _pointMeans = rhs._pointMeans;
          _pointStandardDeviations = rhs._pointStandardDeviations;
  }
      return(*this);
    }
  
  virtual ~MyPointSet()
    {}


  void compareTTest(string name, const MyPointSet &other, MyNT &tstat,
                    MyNT &pvalue) const;
  
  MyNT computeAverageExpressionValue(const set< string > &subset);

  void computeCoordinateDistribution(MyHistogram &histogram) const;

  MyNT computeMean(string point) const
    {
      MyPoint mean;
      computeMean(point, mean);
      return(::computeMean(mean.coords, &mean.missingDims));
    }

  bool computeMean(string point, MyPoint &mean) const;

  bool computeMeanAndVariance(string point, MyNT &mu, MyNT &var,
                              unsigned int &numCoords) const;
  

  
  void computeMeanAndStandardDeviation(MyNT &mean, MyNT &stddev);
  
  void computeSignificantValues(MyNT pvalueThreshold,
                                MyHistogram &histogram, MyNT &value);
  
  void convertToRanks();
  

  void print(string outfile);
  
  void print(ostream &ostr);
  
  // when params can have multiple input files, i need some variable
  // to say which one to read!
  virtual void read(MyClusterParams& params, unsigned int whichPointSet = 0);
  // read just with an infile. call readAffyData() directly.
  virtual void read(string infile, const MyAffyFileFormat &format = MyAffyFileFormat());

  virtual void insert(const MyPoint& point)
    {
      points.push_back(point);
      numPoints++;
    }

  virtual unsigned int size() const
    {
      return(points.size());
    }

  virtual unsigned int getNumPoints() const
    {
      return(numPoints);
    }

  virtual unsigned int getNumDimensions() const
    {
      return(numDimensions);
    }

  virtual unsigned int getNumSamples() const
    {
      if (_isFlipped())
        return(numPoints);
      return(numDimensions);
    }

  virtual unsigned int getNumGenes() const
    {
      if (_isFlipped())
        return(numDimensions);
      return(numPoints);
    }
  virtual bool isPoint(string name) const
    {
      return(namesToIndices.end() != namesToIndices.find(name));
    }
  
  virtual vector< unsigned int > getPointIndex(string name) const
    {
      map< string, vector< unsigned int >, ltstr >::const_iterator itr =
        namesToIndices.find(name);
      if (namesToIndices.end() == itr)
        {
#ifdef DEBUG
//          cerr << "Cannot find point named \"" << name << "\"" << endl;
#endif          
//          return(getNumPoints());
          // return a vector of size 1.
          return(vector<unsigned int>(1, getNumPoints()));
        }
      return(itr->second);
    }
  
  
  virtual const MyPoint& getSample(unsigned int i ) const;
  

  // is there class info stored with the points.
  bool doClassNamesExist() const
    {
      return(classNamesExist);
    }

  // get the class name of sample i. it is better to have a separate
  // function since the result depends on whether the point set is
  // flipped or not. if i didnt have this function and the point set
  // was flipped, i could call points[i].getClassName(). but i would
  // have to know in the calling function whether the point set was
  // flipped or not, which is messy.
  string getClassName(int i) const;
  // get the sample name of sample i.
  string getSampleName(int i) const;
  // get the gene name of gene i.
  string getGeneName(int i) const;
  
  // return the number of points in a particular class.
  unsigned int numInClass(string className) const
    {
      map< string, int, ltstr >::const_iterator itr;
      itr = classCounts.find(className);
      return((*itr).second);
//      return(classCounts[className]);
    }
  
  // set the class names.
  void setClasses(string className);
  void setClasses(const MyClassesType& classes);
  // read classes from a file. return the number of classes read.
  unsigned int readClasses(string fileName);

  void getClasses(MyClassNamesType &names)
    {
      names = classNames;
    }
  
  
  // return the info stored about all dimensions.
  virtual const vector< MyDimInfo >& getDimensionInfo() const
    {
      return(dimInfo);
    }

  // return info stored about dimension i.
  virtual MyDimInfo getDimensionInfo(int i) const
    {
      return(dimInfo[i]);
    }

  virtual void getDimensionNames(set< string > &names) const
    {
      for (unsigned int i = 0; i < getNumDimensions(); i++)
        names.insert(dimInfo[i].getName());
    }

  virtual void getDimCoords(int dim, vector< MyNT >& coords) const
    {
      coords = dimensions[dim];
//       for (unsigned int i = 0; i < numPoints; i++)
//         coords.push_back(points[i][dim]);
    }
  
  virtual MyNT getDimensionLow(int i) const
    {
      return(dimLows[i]);
    }
  virtual MyNT getDimensionHigh(int i) const
    {
      return(dimHighs[i]);
    }

  virtual MyNT getDimensionMean(int i) const
    {
      return(dimMeans[i]);
    }

  virtual MyNT getDimensionSigma(int i) const
    {
      return(dimSigmas[i]);
    }

  // get a vector of the distinct values for each dimension
  virtual void getDimensionValues(unsigned int i, vector< MyNT > &values)
    const;
  // return the number of coords in dimension i with coordinate = value.
  virtual unsigned int getDimensionCount(unsigned int i, MyNT value) const
    {
      return(_getDimensionCount(dimCounts[i].find(value)));
    }
  
  
  
  virtual const MyPoint& operator[](int i) const
    {
      // TODO: BUG. not checking if i is in the right range.
      return(points[i]);
    }
  
  virtual MyPoint& operator[](int i)
    {
      // TODO: BUG. not checking if i is in the right range.
      return(points[i]);
    }

  // find the min and max of the points indexed by discriminant along
  // dimension dim.
  virtual void getInterval(const vector< int > &discIndices, 
                                          unsigned int dim, MyNT &low, MyNT &high) const;


  // suggest a class (to use in MyPointSet::generate()). the current
  // algo simply iterates over all classes.
  virtual string suggestClass() const
  {
    static set< string, ltstr>::const_iterator classNameItr = 
         classNames.begin();
    ++classNameItr;
    // reset, if necessary.
    if (classNames.end() == classNameItr)
         classNameItr = classNames.begin();
    return(*classNameItr);
  }
  // generate a random index.
  virtual bool generate(int& index, bool forbidChosen = true) const;
  // generate a random set of indices.
  virtual bool generate(int num, vector< int >& indices, 
                                    // if true, don't generate indices for hidden points.
                                    bool forbidHidden = true,
                                    // only generate indices of points in this class.
                                    string classToUse = "") const;
  
  virtual int getNumNotChosen() const
    {
      return(numNotChosen);
    }
  
  virtual int getNumNotHidden() const
    {
      return(numNotHidden);
    }
  
  virtual void getPointNames(set< string > &names) const
    {
      for (unsigned int i = 0; i < getNumPoints(); i++)
        names.insert(points[i].getName());
    }

  // flipping coordinates and points.
  virtual void flip();

  void splitByClass(vector< MyPointSet > &splitPointSets,
                    map< string, unsigned int > &classToIndex);
  
  
  void randomise()
    {
      for (unsigned int i = 0; i < points.size(); i++)
        points[i].randomise();
    }
  
      

  
  // functions to set chosen and hidden bits. see comments before
  // definition of MyPoint::chosen to see the difference between
  // hiding and choosing a point.
  //
  // called by MyCluster::chooseContainedPoints().
  virtual void choosePoints(const vector< int >& containedPoints);
  
  // called by MyCluster::hideContainedPoints().
  virtual void hidePoints(const vector< int >& containedPoints);
  // called by MyCluster::hide().
  virtual void hide(const vector< int >& containedPoints,
                                const vector< MyDimensionInterval > &intervals);

  virtual MyNT computeEntropy(int dim, MyVerbosityLevel verbosity) const;
  virtual void computeEntropies(MyVerbosityLevel verbosity);
  // do some simple filtering on points.
  virtual unsigned int filter(const MyClusterParams &params);
  // a version for code that does not want to use MyClusterParams.
  unsigned int filter(MyNT maxDownRegulated, MyNT minUpRegulated);


  virtual bool computeCorrelations(string point1, string point2,
                                   vector< MyNT > &correlations,
                                   string method = "Pearsons") const;
  
  virtual void computeCorrelations(MyNT threshold, string id,
                                   bool absoluteValues = true,
                                   map< string, map< string, MyNT > >*pairs = NULL);
  virtual void printCorrelations(MyNT threshold, string outputFile,
                                 bool absoluteValues = true,
                                 map< string, map< string, MyNT > >*pairs = NULL);

  MyHistogram printCorrelationsHistogram(string outputFile) const;

  MyHistogram computeCorrelationsHistogram(map< string, map< string, MyNT > >*pairs = NULL)
    const;

  virtual MyHistogram
  computeRandomisedCorrelationsHistogram(unsigned int numRandomisations, string outputFile);


  /* Standardise each gene in the data set.

  Linearly transform each genes's coordinates to have mean 0 and standard deviation 1.

  */
  void standardise()
    {
      for (unsigned int i = 0; i < points.size(); i++)
        points[i].standardise();      
    }
  
  /* Standardise each gene in the data set to have mean m and standard deviation s.

  Linearly transform each genes's coordinates to have mean m and standard deviation s.

  */
  void standardise(MyNT m, MyNT s)
    {
      for (unsigned int i = 0; i < points.size(); i++)
        points[i].standardise(m, s);      
    }

  /* Store translations of the row names into a new name space.

     \param[in] rowNameAliases a two dimensional map, each of whose
     keys is a row name, and each of whose values is a set of aliases
     for this row name in the new name space.

     \note The method internally stores these translations, so that a
     subsequent call to MyPointSet::isPoint(name),
     MyPointSet::getPointIndex(name), or MyPointSet::computeMean(name)
     will work correctly when name belongs to the new name space.

     \warning Ensure that the keys in rowNameAliases are from the same
     namespace as the row names stored in the invocant.
   */
  void storeRowNameTranslations(const map< string, set< string > > & rowNameAliases);
  

  


  bool _isFlipped() const
    {
      return(flipped);
    }
  
private:

  // count the number of points that have not yet been chosen in a
  // cluster. see comments before definition of MyPoint::chosen to see
  // the difference between hiding and choosing a point.
  virtual void _countNumNotChosen();

  // count the number of points that have not yet been hidden in a
  // cluster. see comments before definition of MyPoint::chosen to see
  // the difference between hiding and choosing a point.
  virtual void _countNumNotHidden();

  // functions to compute and store stats on dimensions.
  void _fillIn();

  // create a map going to a point's name to its index.
  void _createNameToIndicesMap();

  // functions to return a value or its count for a dimension.
  MyNT _getCoord(MyDimCountsConstPtr itr) const
    {
      return((*itr).first);
    }
  
  void _computeDimCounts();
  
  unsigned int _getDimensionCount(MyDimCountsConstPtr itr) const
    {
      return((*itr).second);
    }

  void _computeClassCounts();

private:  

  unsigned  int numPoints;
  unsigned int numDimensions;
  unsigned int numNotChosen;
  unsigned int numNotHidden;
  vector< MyPoint > points;
  // it is useful to store the coords by dimension too. i am blowing
  // up space by a factor of two but i am saving on time since the
  // calls to MyPointSet::getDimCoords() are much faster.
  vector< vector< MyNT > > dimensions;
  // a hash set of indices corresponding to those points that have been hidden.
#if GCC_VERSION >= 40300
  unordered_set< int > hiddenPoints;
  // for each class, indices that are hidden.
  map< string, unordered_set< int > > hiddenPointsByClass;
#else
  hash_set< int > hiddenPoints;
  // for each class, indices that are hidden.
  map< string, hash_set< int > > hiddenPointsByClass;
#endif
  // have i flipped the point and dimensions?
  bool flipped;
  
  // map from a point's name to its index.
//  map< string, int, ltstr > namesToIndices;

  // for some data sets, e.g., human Affy chips, there are multiple
  // indices for each point, e.g., an entrez-gene ID, so the second
  // type in the map is a vector.
  map< string, vector< unsigned int >, ltstr > namesToIndices;

  bool classNamesExist;

  // just a collection of class names.
  set< string, ltstr> classNames;
  // map from a class to the indices of the points in that class.
  map< string, vector< int >, ltstr > classesToIndices;
  
  // map storing counts for each class.
  map< string, int, ltstr > classCounts;

  

  // something to store info about dimensions.
  // names and ids of dimensions.
  vector< MyDimInfo > dimInfo;
  // smallest coordinate for each dim.
  vector< MyNT > dimLows;
  // greatest coordinate for each dim.
  vector< MyNT > dimHighs;
  // the mean for each dim (useful in computing Gaussian p-values).
  vector< MyNT > dimMeans;
  // the sigma for each dim (useful in computing Gaussian p-values).
  vector< MyNT > dimSigmas;
  
  // for each dim, a count of how many times each value appears in a
  // point along that dim.
  vector< MyDimCountsType > dimCounts;
  
  // entropy of info for each dimension.
  vector< MyNT > dimEntropies;
  // dimension sorted in order of increasing entropy.
  vector< int > sortedDims;

  // point means
  vector< MyNT > _pointMeans;
  // point standard deviations.
  vector< MyNT > _pointStandardDeviations;
  
};


class MySetOfPointSets
{
private:
  // all the point sets.
  vector< MyPointSet > _pointSets;
  // the class/condition associated with each point set.
  vector< string > _classes;
  // the inverse of the previous map.
  map< string, unsigned int > _classToIndex;
  
public:
  // typedef vector< MyPointSet >::iterator MySetOfPointSets::iterator;
  // typedef vector< MyPointSet >::const_iterator MySetOfPointSets::const_iterator;
  
  MySetOfPointSets()
      : _pointSets(), _classToIndex()
    {}
  
  virtual ~MySetOfPointSets()
    {}

  void insert(MyPointSet &p)
    {
      _pointSets.push_back(p);
      // make sure _classes is the same length as _pointSets.
      _classes.push_back("");
    }

  void insert(MyPointSet &p, string condition)
    {
      _pointSets.push_back(p);
      _classes.push_back(condition);
      _classToIndex[condition] = _pointSets.size() - 1;
    }

  void assignClass(unsigned int index, string condition)
    {
      if (index >= _pointSets.size())
        {
          cerr << "Assinging class " << condition
               << "to a non-existant point set index " << index << "." << endl;
          return;
        }
      // forget the old info.
      _classToIndex.erase(_classes[index]);
      // record the new info.
      _classes[index] = condition;
      _classToIndex[condition] = index;
    }
  
  bool getPointSet(string condition, MyPointSet *pointSet)
    {
      map< string, unsigned int >::iterator citr = _classToIndex.find(condition);
      if (_classToIndex.end() == citr)
        {
          pointSet = NULL;
          return(false);
        }
      pointSet = &(_pointSets[citr->second]);
      return(true);
    }

  // /// \param[out] pointSet, 
  const MyPointSet * getPointSet(string condition) const
    {
      map< string, unsigned int >::const_iterator citr = _classToIndex.find(condition);
      if (_classToIndex.end() == citr)
        {
          return(NULL);
          // pointSet = NULL;
          // return(false);
        }
      return(&(_pointSets[citr->second]));
      // pointSet = &(_pointSets[citr->second]);
      // return(true);
    }
};

// next() will return a MyPointSet&.
class MySetOfPointSetsIterator: public  MyIterator< MyPointSet >
  {
  private:
    typedef vector< MyPointSet >::iterator _MySetOfPointSetsIterator;

    _MySetOfPointSetsIterator start;
    _MySetOfPointSetsIterator end;
    _MySetOfPointSetsIterator current;
    vector< MyPointSet > &pointSets;
    
  public:
    // ban this contstructor.
    MySetOfPointSetsIterator(_MySetOfPointSetsIterator _start,
                             _MySetOfPointSetsIterator _end);
    

    MySetOfPointSetsIterator(_MySetOfPointSetsIterator _start,
                             _MySetOfPointSetsIterator _end,
                             vector< MyPointSet > &_pointSets)
        : start(_start), end(_end), current(_start), pointSets(_pointSets)
      {}
    
    virtual bool hasNext()
      {
        return(current != end);
      }
    
    MyPointSet &next()
      {
        MyPointSet &_next = *current;
        current++;
        return(_next);
      }
  };



#endif //POINT_H