#ifndef _HISTOGRAM_H
#define _HISTOGRAM_H

/*
   Histogram creation and manipulation routines.

   Copyright (C) 2007  David S. Smith

   This program 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.

   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
*/    

#include <iostream>
#include <fstream>
#include <valarray>
#include <vector>
#include <string>
#include <math.h>
#include <assert.h>


using std::valarray;
using std::ostream;
using std::ofstream;
using std::string;
using std::cerr;



class Histogram
{

public:

   typedef double        Real;
   typedef unsigned long UInt;
   typedef long          Int;

   Histogram (const Real lower, const Real upper, const Int  size, 
      const bool norm = false, const bool bin_log = false);

   Histogram& insert (const Real datum, const Real weight = 1.0);
   Histogram& clear  ();
   Histogram& save   (const char*  filename);
   Histogram& save   (const string filename);

   Real average    () const { return counts_.sum() / counts_.size(); }
   Real sum        () const { return counts_.sum(); }
   bool normalized () const { return norm_; }
   UInt size       () const { return size_; }

   friend ostream& operator<< (ostream& out, const Histogram& h);
   
   const Int operator[] (const Real& datum) const;
    
   Histogram& operator+= (const Real rhs);
   Histogram& operator-= (const Real rhs);
   Histogram& operator*= (const Real rhs);
   Histogram& operator/= (const Real rhs);
   
   const Histogram operator+ (const Real rhs);
   const Histogram operator- (const Real rhs);
   const Histogram operator* (const Real rhs);
   const Histogram operator/ (const Real rhs);

private:

   valarray<Real> bin_, sum_, counts_, low_part_;
   const bool     norm_, bin_log_;
   UInt           size_;

};


//-- Implementation:

Histogram::Histogram (const Real lower, const Real upper, const Int size, 
   const bool norm, const bool bin_log) 
   : norm_(norm), bin_log_(bin_log), size_(size)
{ 

   assert(bin_log_ ^ (lower <= 0.0 || upper <= 0.0));
   assert(lower < upper);
   
   bin_.resize(size+1, 0.0);
   sum_.resize(size+1, 0.0);
   low_part_.resize(size+1, 0.0);
   counts_.resize(size+1, 0.0);

   if (bin_log_)  // log spacing
      for (UInt i = 0; i < bin_.size(); ++i)
         bin_[i] = pow(10.0, (Real(i) * log10 (upper / lower) / 
            Real(bin_.size() - 1))) + log10 (lower);
   else           // linear spacing
      for (UInt i = 0; i < bin_.size(); ++i) 
         bin_[i] = Real(i) * (upper - lower) / (Real(bin_.size()) - 1.0) 
            + lower;

}


inline Histogram& Histogram::insert (const Real datum, const Real weight) 
{
   Int i = (*this)[datum];
   

#  ifdef HIST_NO_KAHAN
      sum_[i] += weight; 
#  else
      // Kahan summation algorithm to preserve low-order bits;
      // see http://en.wikipedia.org/wiki/Kahansum_mation_algorithm
      Real w       = weight - low_part_[i];
      Real t       = sum_[i] + w;
      low_part_[i] = (t - sum_[i]) - w;
      sum_[i]      = t;
#  endif

   counts_[i]  += 1.0;

   return *this;
}


Histogram& Histogram::clear()
{
   sum_      = 0.0;
   low_part_ = 0.0;
   counts_   = 0.0;
   return *this;
}


inline const long Histogram::operator[] (const Real& datum) const
{
   Int U, L, M;
   
#  ifdef HIST_DIRECT_CALC
   Real n = Real(size_), hi = bin_[n-1], lo = bin_[0];
   if (bin_log_) 
      L = Int(log10(datum/lo) * (n - 1.0) / log10(hi / lo));
   else
      L = Int((datum - lo) * (n - 1.0) / (hi - lo));
#  else
   // binary search for nearest bin (see Knuth, TAOCPv3, p. 410).
   // Directly calculating the bin would be faster but less generic.
   U = bin_.size() - 1;
   L = 0;
   do 
   {
      M = (U + L) / 2;
      if (datum < bin_[M])
         U = M;
      else 
         L = M;
   } while (U - L > 1);
#  endif
   
   return L;
}


inline Histogram& Histogram::operator+= (const Real rhs) 
   { sum_ += rhs; return *this; }

inline Histogram& Histogram::operator-= (const Real rhs) 
   { sum_ -= rhs; return *this; }
   
inline Histogram& Histogram::operator*= (const Real rhs) 
   { sum_ *= rhs; return *this; }

inline Histogram& Histogram::operator/= (const Real rhs) 
   { sum_ /= rhs; return *this; }

inline const Histogram Histogram::operator+ (const Real rhs) 
   { return Histogram(*this) += rhs; }

inline const Histogram Histogram::operator- (const Real rhs) 
   { return Histogram(*this) -= rhs; }

inline const Histogram Histogram::operator* (const Real rhs) 
   { return Histogram(*this) *= rhs; }

inline const Histogram Histogram::operator/ (const Real rhs) 
   { return Histogram(*this) /= rhs; }


inline Histogram& Histogram::save (const string filename)
{
   this->save(filename.c_str());
   return *this;
}


inline Histogram& Histogram::save (const char* filename) 
{
   ofstream os(filename);
   os << *this;
   os.close();
   return *this;
}

inline ostream& operator<< (ostream& out, const Histogram& h)
{
   out << "# bin sum counts\n";
   
   if (h.normalized())
      for (UInt i = 0; i < h.size(); ++i)
         out << h.bin_[i] << " " << h.sum_[i] / (h.bin_[i+1] - h.bin_[i])
         << " " << h.counts_[i] << "\n";
   else
      for (UInt i = 0; i < h.size(); ++i)
         out << h.bin_[i] << " " << h.sum_[i] << " " << h.counts_[i] << "\n";
   
   return out;
}


#endif /* _HISTOGRAM_H */