source: src/Fragmentation/SetValues/Histogram.cpp@ 7e19060

/*
 * 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;
}

std::ostream & operator<<(std::ostream &ost, const Histogram &histogram)
{
  for (Histogram::Bins_t::const_iterator iter = histogram.bins.begin();
      iter != histogram.bins.end(); ++iter)
    ost << *iter;
  return ost;
}

Histogram::Histogram(const samples_t &samples) :
  binwidth(0.5),
  offset(0.)
{
  if (!samples.empty()) {
    // set offset to first value
    const_cast<BinLowerEnd &>(offset) = *samples.begin();
    // set binwidth to statistical sensible value
    MinMax_t MinMax = getMinMaxFromSamples(samples);
    const_cast<double &>(binwidth) = (*(MinMax.second) - *(MinMax.first))/pow(samples.size(), 1./3.);
    // and add all samples to these histogram bins
    addSamples(samples);
  }
}

Histogram::Histogram(const samples_t &samples, const BinLowerEnd _offset, const double _binwidth) :
    binwidth(_binwidth),
    offset(_offset)
{
  addSamples(samples);
}

Histogram::MinMax_t Histogram::getMinMaxFromSamples(const samples_t &samples) const
{
  std::pair<samples_t::const_iterator, samples_t::const_iterator> returnpair;
  returnpair.first = min_element(samples.begin(), samples.end());
  returnpair.second = max_element(samples.begin(), samples.end());
  ASSERT((returnpair.second != samples.end()) || (returnpair.first != samples.end()),
      "Histogram::Histogram() - cannot find min/max despite non-empty range.");
  return returnpair;
}

void Histogram::addSamples(const samples_t &samples)
{
  // build histogram from samples
  if (!samples.empty()) {
    // 1. get min and max and determine width
    MinMax_t MinMax = getMinMaxFromSamples(samples);
//    LOG(1, "DEBUG: min is " << *(MinMax.first) << " and max is " << *(MinMax.second) << ".");

    // 2. create each bin
    {
      std::vector<Bin_t> vectorbins;
      const BinLowerEnd HistogramStart = getLowerEnd(*(MinMax.first));
      BinCreator_t BinCreator( HistogramStart, binwidth );
      // we need one extra bin for get...Bin()'s find to work properly
      const int CountBins = ceil((*(MinMax.second) - HistogramStart)/binwidth)+1;
      vectorbins.resize(CountBins+1, Bin_t( make_pair(0., 0.) ) );
      std::generate( vectorbins.begin(), vectorbins.end(), BinCreator );
      bins.insert(vectorbins.begin(), vectorbins.end());
    }

    // 3. 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.");
      // (make bin count normalized, i.e. always equal to surface area of 1)
      biniter->second += 1./binwidth;
    }
    LOG(2, "DEBUG: Bins are " << printBins() << ".");
  } else {
//    LOG(1, "INFO: Given vector of samples is empty.");
  }
}

std::string Histogram::printBins() const
{
  std::stringstream output;
  std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
  return output.str();
}

void Histogram::extendMissingBins(const BinLowerEnd LowerEnd, const BinLowerEnd NextLowerEnd)
{
  Bins_t::const_iterator loweriter = getLowerEndBin(LowerEnd);
  Bins_t::const_iterator upperiter = getHigherEndBin(NextLowerEnd);
  if (loweriter == bins.end()) {
    // we need bins beneath our first, add them
    for(BinLowerEnd offset = getLowerEnd(LowerEnd);
        offset < getLowerEnd(NextLowerEnd);
        offset += binwidth) {
      LOG(4, "DEBUG: Extending below at offset " << offset << ".");
      bins.insert( std::make_pair (offset, 0. ) );
    }
    LOG(3, "DEBUG: Bins after adding empties before start are " << printBins() << ".");
    loweriter = getLowerEndBin(LowerEnd);
    ASSERT( loweriter != bins.end(),
        "Histogram::extendMissingBins() - still no lower bound after adding empties.");
  }
  if (upperiter == bins.end()) {
    // we need bins after our last, add them
    for(BinLowerEnd offset = getLowerEnd(LowerEnd);
        offset <= getLowerEnd(NextLowerEnd)+(1.5*binwidth); /* for safety we go a little further */
        offset += binwidth) {
      LOG(4, "DEBUG: Extending above at offset " << offset << ".");
      bins.insert( std::make_pair (offset, 0. ) );
    }
    LOG(3, "DEBUG: Bins after adding empties after end are " << printBins() << ".");
    upperiter = getHigherEndBin(NextLowerEnd);
    ASSERT( upperiter != bins.end(),
        "Histogram::extendMissingBins() - still no upper bound after adding empties.");
  }
}

void Histogram::superposeOtherHistogram(const Histogram &other, const double prefactor)
{
  // go through each of the other histogram's bins
  Bins_t::const_iterator enditer = --other.bins.end(); // (except internal last one)
  for (Bins_t::const_iterator biniter = other.bins.begin();
      biniter != enditer; /* we advance ourselves in loop */) {
    const Bin_t &bin = *biniter;
    ++biniter;
    const Bin_t &nextbin = *biniter;

    LOG(4, "DEBUG: Current bin is " << bin << ", next bin is " << nextbin << ".");

    // Check first whether start or end actually fit into our histogram, if not extend.
    extendMissingBins(bin.first, nextbin.first);

    // The bin will in general not fit into one bin in this histogram, but overlap.
    // Hence, we determine the contribution of the bin to each bin in this histogram
    // its overlaps into and add this weight to the bin.
    Bins_t::const_iterator loweriter = getLowerEndBin(bin.first);
    Bins_t::const_iterator upperiter = getHigherEndBin(nextbin.first);

    ASSERT( loweriter->first < upperiter->first,
        "Histogram::superposeOtherHistogram() - the bin range is invalid.");
    LOG(4, "DEBUG: bin range here is ["
        << loweriter->first << ","  << upperiter->first << ").");

    // Next, we create a vector of offsets
    typedef std::vector< BinLowerEnd > offsets_t;
    offsets_t offsets;
    {
      offsets.push_back(bin.first);
      Bins_t::const_iterator iter = loweriter;
      for (++iter; iter != upperiter; ++iter)
        if (offsets.back() != iter->first)
          offsets.push_back(iter->first);
      if (offsets.back() != nextbin.first)
        offsets.push_back(nextbin.first);
      LOG(3, "DEBUG: Offsets are " << offsets << ".");
    }

    // then, we go through the offsets but the last one and add the respective area
    {
      offsets_t::const_iterator iter = offsets.begin();
      offsets_t::const_iterator nextiter = ++offsets.begin();
      for (; iter != --offsets.end(); ++iter, ++nextiter) {
        const double length = *nextiter - *iter;
        const double weight = bin.second * (length/binwidth);
        Bins_t::iterator filliter = getLowerEndBin(*iter);
        filliter->second += prefactor * weight;
      }
    }

    {
      std::stringstream output;
      std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
      LOG(2, "DEBUG: Bins are after summation " << output.str() << ".");
    }
  }
}

Histogram& Histogram::operator=(const Histogram &other)
{
  // check for self-assigment
  if (&other != this) {
    bins.clear();
    const_cast<BinLowerEnd &>(offset) = other.offset;
    const_cast<double &>(binwidth) = other.binwidth;
    bins.insert(other.bins.begin(), other.bins.end());
  }
  return *this;
}

Histogram& Histogram::operator+=(const Histogram &other)
{
  superposeOtherHistogram(other, +1.);
  LOG(2, "DEBUG: Bins after addition are " << printBins() << ".");
  return *this;
}

Histogram& Histogram::operator-=(const Histogram &other)
{
  superposeOtherHistogram(other, -1.);
  LOG(2, "DEBUG: Bins after subtraction are " << printBins() << ".");
  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;
}

Histogram::BinLowerEnd Histogram::getLowerEnd(const double _value) const
{
  BinLowerEnd start = _value - offset;
  // then divide by binwidth
  const int integral = floor(start/binwidth);
  //const double fraction = fmod(start,binwidth);
  start = offset + binwidth*integral;
  return start;
}

double Histogram::area() const
{
  double area = 0.;
  for(Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
    area += iter->second*binwidth;
  return area;
}

template<> Histogram ZeroInstance<Histogram>()
{
  Histogram returnvalue;
  return returnvalue;
}