---

histogram.h

Go to the documentation of this file.

#ifndef __HISTOGRAM_H__
#define __HISTOGRAM_H__


#include <iterator>
#include <functional>
#include <vector>
#include <numeric>
#include <complex>

#include "../src/index_iterator.h"
#include "../src/genericalgs.h"
#include "../src/generic_stdlib.h"
#include "../src/bounded.h"
#include "../src/mathconstants.h"


namespace  Histo {


template<class Count=int, class Value=double>
class OneDBinList {
protected:
  std::vector<Count> counts;
  std::vector<Value> values;
  Count total_c;
  Value total_v;

  template<class C, class V>  
  friend OneDBinList<C,V> operator-(const OneDBinList<C,V> &A, const OneDBinList<C,V> &B);
  
private:
  typedef Generic::index_iterator< std::vector<Value>, Scale::Linear<Value> > index_i;
  typedef Scale::Linear<Value> ScaleType;
  
public:
  typedef typename std::vector<Value>::size_type size_type;
  typedef Value value_type;

  typedef Bounded::Magnitude Magnitude;

  typedef std::complex<Value> Vector;

  static size_type num_bins_default;

  OneDBinList(const size_type num_bins=num_bins_default)
    : counts(num_bins), values(num_bins), total_c(0), total_v(0),
      cached_vector_sum_count(0), cached_vector_sum(0,0),
      cached_sum_value_count(0),  cached_sum_value(0)
    {  assert(num_bins>0);  }


  // Could specify as h[bin]+=val instead
  OneDBinList& add(size_type bin, Value val=1){
    total_c++;
    total_v+=val;
    counts[bin]++;
    values[bin]+=val;
    return *this;
  }

  OneDBinList& keep_peak(size_type bin, Value val) {
    total_c++;
    total_v+=val;
    counts[bin]++;
    if (val > values[bin]) values[bin]=val; 
    return *this;
  }

  inline size_type num_bins() const {  return counts.size();  }

  inline Value sum_value() const {
    if (total_count() != cached_sum_value_count) {
      cached_sum_value = std::accumulate(ISEQ(values),Value(0));
      cached_sum_value_count = total_count();
    }
    return cached_sum_value;
  }
  
  inline Count sum_count()    const {  return std::accumulate(ISEQ(counts),Count(0));  }

  inline Value average_value()const {  return sum_value()/Value(num_bins());  }

  Value max_value()    const {  return *(std::max_element(ISEQ(values)));  }
  Value min_value()    const {  return *(std::min_element(ISEQ(values)));  }
  Count max_count()    const {  return *(std::max_element(ISEQ(counts)));  }
  Count min_count()    const {  return *(std::min_element(ISEQ(counts)));  }

  size_type max_value_bin() const
    {  return std::max_element(ISEQ(values))-values.begin();  }

  //Magnitude max_value_bin_mag() const {  return bin_mag(max_value_bin());  }
  
  inline Magnitude relative_selectivity() const {
    /* With only one bin, one is always fully selective; special
       case to avoid divide-by-zero errors.  */
    if (num_bins()<=1) return 1.0;
    const Magnitude proportion = safe_divide(value(max_value_bin()),sum_value());
    const Magnitude offset = 1.0/num_bins();
    const Magnitude scaled = (proportion-offset)/(1.0-offset);
    /* Since we take the peak, negative selectivities shouldn't be
       possible, but just in case... */
    return Generic::crop(Magnitude(0.0),Magnitude(1.0),scaled);
    return scaled;
  }

  Value weighted_sum() const
    {  return std::inner_product(ISEQ(values),index_i(values.begin()),Value(0));  }

  inline Magnitude weighted_average() const
    {  return safe_divide((safe_divide(weighted_sum(),sum_value())),Value(num_bins()-1));  }


  inline static size_type undefined_values() {  return undefined_vals;  }
  inline static void      undefined_values_reset(size_type val=0) {  undefined_vals=val;  }
  
  inline Value total_value()  const {  return total_v;  }
  
  inline Count total_count()  const {  return total_c;  }
  
  

  
  Vector vector_sum() const {
    if (total_count() != cached_vector_sum_count) {
      cached_vector_sum = std::inner_product
        (ISEQ(values), index_i(values.begin(),ScaleType(2*(Math::pi)/Value(num_bins()))),
         Vector(0,0),

/* Workaround for bug in at least egcs-1.1.2, possibly GCC 2.95 as well */       
#if (__GNUC__ && (__GNUC__ < 3))
         plus<Vector>(),
#else
         std::plus<Vector>(),
#endif
         
         Generic::Polar<Value>());
      cached_vector_sum_count = total_count();
    }
    
    return cached_vector_sum;
  }

  inline Vector vector_average() const {  return vector_sum()/Value(num_bins());  }

  inline Magnitude vector_direction() const {
    if (std::abs(vector_sum())==0) undefined_vals++; /* Optional check */
    return Generic::wrap(0.0,1.0,std::arg(vector_sum())/(2*(Math::pi)));
  }

  inline Magnitude vector_selectivity() const
    {  return safe_divide(std::abs(vector_sum()),sum_value());  }


  
  inline Value         value(size_type i) const {  return values[i];  }
  inline Count         count(size_type i) const {  return counts[i];  }
  inline Magnitude value_mag(size_type i) const {  return safe_divide(values[i]*1.0,total_value());  }
  inline Magnitude count_mag(size_type i) const {  return safe_divide(Magnitude(counts[i]),Magnitude(total_count()));  }


  Magnitude relative_value(size_type bin) const
    {  return safe_divide(value(bin),sum_value());  }

  
  size_type bin_number(Magnitude binmag) {
    return Generic::crop(size_type(0),size_type(num_bins())-size_type(1),
                         size_type(binmag*num_bins()+float_adjust));
  }


  size_type bin_number_inclusive(Magnitude binmag) {
    return Generic::crop(size_type(0),size_type(num_bins()),
                         size_type(binmag*(num_bins()-1)+float_adjust));
  }
  

  Magnitude bin_mag(size_type bin) {  return Magnitude(bin)/Magnitude(num_bins());  }

private:

  template <class T>
  inline T safe_divide(const T numerator, const T denominator) const 
    {  return (denominator==0? undefined_vals++, 0 : numerator/denominator);  }

  static size_type undefined_vals;

  mutable Count  cached_vector_sum_count;
  mutable Vector cached_vector_sum;
  mutable Count  cached_sum_value_count;
  mutable Value  cached_sum_value;

  static const value_type float_adjust; 
};



/* Static variables */
template<class Count, class Value> typename OneDBinList<Count,Value>::size_type
OneDBinList<Count,Value>::num_bins_default=1;  
template<class Count, class Value> typename OneDBinList<Count,Value>::size_type
OneDBinList<Count,Value>::undefined_vals=0;

template<class Count, class Value> const typename OneDBinList<Count,Value>::value_type
OneDBinList<Count,Value>::float_adjust=1E-14; 


template<class Count, class Value>
OneDBinList<Count,Value> operator-(const OneDBinList<Count,Value> &A, const OneDBinList<Count,Value> &B)
{
  typename OneDBinList<Count,Value>::size_type N = A.num_bins();
  assert(N==B.num_bins());

  OneDBinList<Count,Value> tmp(N);

  for (typename OneDBinList<Count,Value>::size_type i=0; i<N; i++) {
    tmp.values[i] = A.values[i] - B.values[i];
    tmp.counts[i] = A.counts[i] - B.counts[i];
  }
  
  tmp.total_v = A.total_v + B.total_v;
  tmp.total_c = A.total_c + B.total_c;
  
  return tmp;
}



template<class Count=int, class Value=double>
class OneDHistogram : public OneDBinList<Count,Value> {
public:
  typedef typename OneDBinList<Count,Value>::size_type size_type;
  OneDHistogram(const size_type num_bins)
    : OneDBinList<Count,Value>(num_bins) { }
};


 
template<class Count=int, class Value=double>
class OneDHistogramDifference : public OneDBinList<Count,Value> {
public:
  OneDHistogramDifference(const OneDBinList<Count,Value> &A,
                          const OneDBinList<Count,Value> &B)
    : OneDBinList<Count,Value>(A-B) { }
};
 

 
} /* namespace Histo */
#endif /* __HISTOGRAM_H__ */

---