histogram.h

/**************************************************************************
 * Copyright (c) 2004-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 _HISTOGRAM_H
#define _HISTOGRAM_H

#include <vector>

//#include "global.h"
#include "enrichment.h"

using namespace std;

inline MyNT computeMean(const vector< MyNT > &coords, set< unsigned int > *missingIndices = NULL)
{
  MyNT numCoords = coords.size(), sum = 0;
  for (unsigned int i = 0; i < numCoords; i++)
    if (!missingIndices || (missingIndices->end() == missingIndices->find(i)))
      sum += coords[i];
  if (missingIndices)
    numCoords -= missingIndices->size();
  return(sum/numCoords);
}

inline MyNT computeVariance(const vector< MyNT > &coords, MyNT mean, set< unsigned int > *missingIndices = NULL)
{
  MyNT numCoords = coords.size(), sum = 0;
  for (unsigned int i = 0; i < numCoords; i++)
    if (!missingIndices || (missingIndices->end() == missingIndices->find(i)))
      sum += (coords[i] - mean)*(coords[i] - mean);
  if (missingIndices)
    numCoords -= missingIndices->size();
  return(sum/numCoords);
}

MyNT studentsTDistributionFunction(MyNT t, MyNT nu);

MyNT studentsTDistributionFunctionTail(MyNT t, MyNT nu);

MyNT studentsTDistributionFunctionTwoSidedTail(MyNT t, MyNT nu);




class MyHistogram
{
private:
  unsigned int _numBins;
  // each bin can store a fraction, so the type is MyNT
  vector< MyNT > _bins;
  MyNT _low;
  MyNT _high;
  // type is MyNT for the same reason as _bins.
  MyNT _numTotalValues;
//  unsigned int _numTotalValues;
public:

  MyHistogram(unsigned int numBins = 10, MyNT low = 0, MyNT high = 1)
      : _numBins(numBins), _bins(numBins, 0), _low(low), _high(high),
        _numTotalValues(0)
    {}

  virtual ~MyHistogram()
    {}

  virtual MyNT computeMean();
  
  virtual MyNT computeStandardDeviation(MyNT mean);

  virtual MyNT computeVariance(MyNT mean);

  
  virtual unsigned int getNumBins() const
    {
      return(_numBins);
    }
  
  virtual MyNT getBin(MyNT coord) const
    {
      // first find the index of the bin. then that bin's lower value.
      return(_binToValue(_valueToBin(coord)));
    }

  virtual void setBin(MyNT coord, MyNT value)
    {
//       cout << "coord is " << coord << ", value is " << value
//            << ", bin is " << _valueToBin(coord) << endl;
      _numTotalValues -= _bins[_valueToBin(coord)];
      _numTotalValues += value;
      _bins[_valueToBin(coord)] = value;
    }
  
  
  virtual MyNT getValue(MyNT coord) const
    {
      return(_bins[_valueToBin(coord)]);
    }
  

  MyNT getMin() const
    {
      return(_low);
    }
  
  MyNT getMax() const
    {
      return(_high);
    }
  
  MyNT getRange() const
    {
      return(getMax() - getMin());
    }
  
  virtual void insert(MyNT value);

  virtual void insert(const vector< MyNT > &values);
  
  // insert a vector of unsigned ints into the histogram
  // TODO eventually template MyHistogram so that this isn't necessary
  virtual void insert(const vector< unsigned int > &values);

  virtual void abs();
  virtual void accumulate(bool reverse = false);

  void constructTDistributionHistogram(unsigned int dof);

  virtual MyNT convertToPvalues(MyNT pvalueThreshold = 0.01, LIBENRICHMENT_TEST_TYPE testType = LIBENRICHMENT_NONE,
                                MyHistogram *reference = NULL);
  
  virtual void convertToPvalues(MyNT pvalueThreshold, const set< LIBENRICHMENT_TEST_TYPE > &testTypes,
                                map< LIBENRICHMENT_TEST_TYPE, MyNT > &thresholds,
                                map< LIBENRICHMENT_TEST_TYPE, MyHistogram > &histograms,
                                MyHistogram *reference = NULL);
  
  virtual const MyHistogram& operator=(const MyHistogram &other);

  virtual MyHistogram operator+(const MyHistogram &other);

  virtual MyHistogram operator-(MyNT weight);

  virtual MyHistogram operator-();

  virtual const MyHistogram &operator+=(const MyHistogram &other);

  virtual MyHistogram operator*(MyNT weight);

  virtual void operator*=(MyNT weight);

  virtual MyHistogram operator/(MyHistogram &other);
  
  virtual const MyHistogram &operator/=(MyHistogram &other);

  virtual void normalise();
  
  virtual MyNT getNumTotalValues() const
    {
      return(_numTotalValues);
    }

  virtual void print(ostream &ostr, string extra = "", vector< MyHistogram > *others = NULL);

  virtual void print(string file, string extra = "", vector< MyHistogram > *others = NULL);

  virtual void read(string infile);
  
  // compute the trapezoidal area under the curve generated by this histogram and other
  // this histogram must be monotonically increasing and other must have corresponding bins
  virtual double computeTrapezoidalAUC(const MyHistogram &other) const
  {
                double area = 0;
          for (unsigned int i = 1; i < getNumBins(); i++)
                        area += 0.5*(_bins[i] - _bins[i - 1])*(other._bins[i] + other._bins[i - 1]);
                return (area);
        }
  
  
private:
  MyNT _binToValue(unsigned int i) const
    {
      return((i*1.0/getNumBins())*getRange() + getMin());
    }
  unsigned int _valueToBin(MyNT value) const
    {
      if (value < getMin())
        {
//          cerr << "\tInserting value " << value << " into histogram whose minimum value is " << getMin() << ". Inserting " << getMin() << "." << endl;
          value = getMin();
        }
      if (value > getMax())
        {
//          cerr << "\tInserting value " << value << " into histogram whose maximum value is " << getMax() << ". Inserting " << getMax() << "." << endl;
          value = getMax();
        }
//       cout << "\tvalue is " << value
//            << ", bin is " << getNumBins()*(value - getMin())/getRange()
//            << ", int(bin) is " << int(getNumBins()*(value - getMin())/getRange())
//            << ", static cast to unsigned int is "
//            << static_cast<unsigned int>(getNumBins()*(value - getMin())/getRange())
//            << endl;
      return(static_cast<unsigned int>(
               // adding 0.5 to perform rounding. in some cases, i
               // compute the previous bin, otherwise.
               0.5
               + getNumBins()*(value - getMin())/getRange()));
    }
};



#endif // _HISTOGRAM_H