source: src/Fragmentation/Histogram/Histogram.cpp@ 3f6bbb

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder 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 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder 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 MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Histogram.cpp
 *
 *  Created on: Jul 26, 2012
 *      Author: heber
 */


// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "Histogram.hpp"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <sstream>
#include <string>

#include <boost/bind.hpp>
#include <boost/foreach.hpp>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

/** Tiny helper struct to create equally spaced bins with count of zero.
 *
 */
struct BinCreator_t {
  BinCreator_t( const double lowerend, const double _width ) :
    currentstart(lowerend),
    width(_width)
  {}
  Histogram::Bin_t operator()() {
    Histogram::Bin_t bin( make_pair( currentstart, 0.) );
    currentstart+=width;
    return bin;
  }
private:
  double currentstart;
  const double width;
};


// see http://stackoverflow.com/questions/634087/stdcopy-to-stdcout-for-stdpair
// printing a pair is not easy, especially so with ostream_iterator as it cannot find
// overload operator<<() as it is not in namespace std.
template <class T>
struct PrintPair : public std::unary_function<T, void>
{
    std::ostream& os;
    PrintPair(std::ostream& strm) : os(strm) {}

    void operator()(const T& elem) const
    {
      os << "(" << elem.first << "," << elem.second << ") ";
    }
};

std::ostream & operator<<(std::ostream &ost, const Histogram::Bin_t &elem)
{
  ost << "(" << elem.first << "," << elem.second << ") ";
  return ost;
}

Histogram::Histogram(const samples_t &samples, const int _CountBins) :
    binwidth(0.5),
    CountBins(_CountBins)
{
  // build histogram from samples
  if (!samples.empty()) {
    if (CountBins == -1) {
      CountBins = ceil(pow(samples.size(), 1./3.));
    }
    // 2. get min and max and determine width
    samples_t::const_iterator maxiter = max_element(samples.begin(), samples.end());
    samples_t::const_iterator miniter = min_element(samples.begin(), samples.end());
    ASSERT((maxiter != samples.end()) || (miniter != samples.end()),
        "Histogram::Histogram() - cannot find min/max despite non-empty range.");
      LOG(1, "DEBUG: min is " << *miniter << " and max is " << *maxiter << ".");
    binwidth = (*maxiter - *miniter)/(CountBins-1);

    // 3. create each bin
    BinCreator_t BinCreator( *miniter, binwidth );
    std::vector<Bin_t> vectorbins;
    // we need one extra bin for get...Bin()'s find to work properly
    vectorbins.resize(CountBins+1, Bin_t( make_pair(0., 0.) ) );
    std::generate( vectorbins.begin(), vectorbins.end(), BinCreator );
    bins.insert(vectorbins.begin(), vectorbins.end());

    // 4. place each sample into bin
    BOOST_FOREACH( double value, samples) {
      const Bins_t::iterator biniter = getLowerEndBin(value);
      ASSERT( biniter != bins.end(),
          "Histogram::Histogram() - cannot find bin for value from given samples.");
      biniter->second += 1.;
    }
    std::stringstream output;
    std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
    LOG(2, "DEBUG: Bins are " << output.str() << ".");
  } else {
    LOG(1, "INFO: Given vector of samples is empty.");
  }
}

Histogram& Histogram::operator+=(const Histogram &other)
{
  return *this;
}

Histogram& Histogram::operator-=(const Histogram &other)
{
  return *this;
}

bool Histogram::isEmpty() const
{
  bool status = true;
  for (Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
    status &= iter->second == 0;
  return status;
}

Histogram::Bins_t::iterator Histogram::getLowerEndBin(const double _value)
{
  // lower bound returns key that is equal or greater
  Bins_t::iterator iter = bins.lower_bound(_value);
  if (iter != bins.end()) {
    // if we are not on the boundary we always have to step back by one
    if (_value != iter->first) {
      if (iter != bins.begin()) {
          --iter;
      } else {
        iter = bins.end();
      }
    } else if (iter == --bins.end()) {
      // if we actually are on boundary of "last bin", set to end
      iter = bins.end();
    }
  }
  return iter;
}

Histogram::Bins_t::iterator Histogram::getHigherEndBin(const double _value)
{
  // upper bound return key that is strictly greater
  Bins_t::iterator iter = bins.upper_bound(_value);
  // if we are on the boundary we have to step back by one
  if (iter != bins.end())
    if (_value == iter->first)
      if (iter != bins.begin())
          --iter;

  return iter;
}