Index: src/Fragmentation/SetValues/SamplingGrid.cpp
===================================================================
--- src/Fragmentation/SetValues/SamplingGrid.cpp	(revision 7b049409753f488e726192ddd4e7610aed7e12b5)
+++ src/Fragmentation/SetValues/SamplingGrid.cpp	(revision 7b049409753f488e726192ddd4e7610aed7e12b5)
@@ -0,0 +1,599 @@
+/*
+ * 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/>.
+ */
+
+/*
+ * SamplingGrid.cpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// include headers that implement a archive in simple text format
+// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "Fragmentation/SetValues/SamplingGrid.hpp"
+
+#include <boost/bind.hpp>
+#include <limits>
+
+#include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
+
+// static instances
+const double SamplingGrid::zeroOffset[3] = { 0., 0., 0. };
+
+SamplingGrid::SamplingGrid() :
+    SamplingGridProperties()
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const double _begin[3],
+    const double _end[3],
+    const int _level) :
+  SamplingGridProperties(_begin, _end, _level)
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const double _begin[3],
+    const double _end[3],
+    const int _level,
+    const sampledvalues_t &_sampled_grid) :
+  SamplingGridProperties(_begin, _end, _level),
+  sampled_grid(_sampled_grid)
+{
+  setWindowSize(_begin, _end);
+  ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
+  SamplingGridProperties(_grid),
+  sampled_grid(_grid.sampled_grid)
+{
+  setWindowSize(_grid.begin_window, _grid.end_window);
+  ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
+  SamplingGridProperties(_props)
+{
+  setWindowSize(zeroOffset, zeroOffset);
+  ASSERT( getWindowGridPoints() == (size_t)0,
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::SamplingGrid(
+    const SamplingGridProperties &_props,
+    const sampledvalues_t &_sampled_grid) :
+  SamplingGridProperties(_props),
+  sampled_grid(_sampled_grid)
+{
+  setWindowSize(_props.begin, _props.end);
+  ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
+      "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
+}
+
+SamplingGrid::~SamplingGrid()
+{}
+
+bool SamplingGrid::isCongruent(const SamplingGrid &_props) const
+{
+  bool status = true;
+  status &= (static_cast<const SamplingGridProperties &>(*this) ==
+      static_cast<const SamplingGridProperties &>(_props));
+  for(size_t i = 0; i<3; ++i) {
+    status &= begin_window[i] == _props.begin_window[i];
+    status &= end_window[i] == _props.end_window[i];
+  }
+  return status;
+}
+
+SamplingGrid& SamplingGrid::operator=(const SamplingGrid& other)
+{
+  // check for self-assignment
+  if (this != &other) {
+    static_cast<SamplingGridProperties &>(*this) = other;
+    setWindowSize(other.begin_window, other.end_window);
+    sampled_grid = other.sampled_grid;
+  }
+  return *this;
+}
+
+static void multiplyElements(
+    double &dest,
+    const double &source,
+    const double prefactor)
+{
+  dest *= prefactor*(source);
+}
+
+
+SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
+{
+  // check that grids are compatible
+  if (isCompatible(other)) {
+    /// get minimum of window
+    double min_begin_window[3];
+    double min_end_window[3];
+    bool doShrink = false;
+    for (size_t index=0; index<3;++index) {
+      if (begin_window[index] <= other.begin_window[index]) {
+        min_begin_window[index] = other.begin_window[index];
+        doShrink = true;
+      } else {
+        min_begin_window[index] = begin_window[index];
+      }
+      if (end_window[index] <= other.end_window[index]) {
+        min_end_window[index] = end_window[index];
+      } else {
+        min_end_window[index] = other.end_window[index];
+        doShrink = true;
+      }
+    }
+    LOG(2, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
+    LOG(2, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
+    if (doShrink)
+      shrinkWindow(min_begin_window, min_end_window);
+    addWindowOntoWindow(
+        other.begin_window,
+        other.end_window,
+        begin_window,
+        end_window,
+        sampled_grid,
+        other.sampled_grid,
+        boost::bind(multiplyElements, _1, _2, 1.),
+        sourcewindow);
+  } else {
+    ASSERT(0, "SamplingGrid::operator*=() - multiplying incongruent grids is so far not in the cards.");
+  }
+  return *this;
+}
+
+void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
+{
+  /// check that grids are compatible
+  if (isCompatible(other)) {
+    /// get maximum of window
+    double max_begin_window[3];
+    double max_end_window[3];
+    bool doExtend = false;
+    for (size_t index=0; index<3;++index) {
+      if (begin_window[index] >= other.begin_window[index]) {
+        max_begin_window[index] = other.begin_window[index];
+        doExtend = true;
+      } else {
+        max_begin_window[index] = begin_window[index];
+      }
+      if (end_window[index] >= other.end_window[index]) {
+        max_end_window[index] = end_window[index];
+      } else {
+        max_end_window[index] = other.end_window[index];
+        doExtend = true;
+      }
+    }
+    LOG(2, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
+    LOG(2, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
+    if (doExtend)
+      extendWindow(max_begin_window, max_end_window);
+    /// and copy other into larger window, too
+    addOntoWindow(other.begin_window, other.end_window, other.sampled_grid, prefactor);
+  } else {
+    ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
+  }
+}
+
+const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const
+{
+  static const double round_offset(
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.); // need offset to get to round_toward_nearest behavior
+//  const double total = getTotalLengthPerAxis(axis);
+  const double delta = getDeltaPerAxis(axis);
+  if (delta == 0)
+    return 0;
+  const double length = getWindowLengthPerAxis(axis);
+  if (length == 0)
+    return 0;
+  return (size_t)(length/delta+round_offset);
+}
+
+double SamplingGrid::integral() const
+{
+  const double volume_element = getVolume()/(double)getTotalGridPoints();
+  double int_value = 0.;
+  for (sampledvalues_t::const_iterator iter = sampled_grid.begin();
+      iter != sampled_grid.end();
+      ++iter)
+    int_value += *iter;
+  int_value *= volume_element;
+  LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
+  return int_value;
+}
+
+double SamplingGrid::integral(const SamplingGrid &weight) const
+{
+  if (isCompatible(weight)) {
+    const double volume_element = getVolume()/(double)getTotalGridPoints();
+    double int_value = 0.;
+    sampledvalues_t::const_iterator iter = sampled_grid.begin();
+    sampledvalues_t::const_iterator weightiter = weight.sampled_grid.begin();
+    for (;iter != sampled_grid.end();++iter,++weightiter)
+      int_value += (*weightiter) * (*iter);
+    int_value *= volume_element;
+    //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
+    return int_value;
+  } else
+    return 0.;
+}
+
+double SamplingGrid::getNearestLowerGridPoint(
+    const double value, const size_t axis) const
+{
+  const double length = end[axis] - begin[axis];
+  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
+    return begin[axis];
+  if (value > end[axis])
+    return end[axis];
+  const double offset = value - begin[axis];
+  if (offset < 0)
+    return begin[axis];
+  // modify value a little if value actually has been on a grid point but
+  // perturbed by numerical rounding: here up, as we always go lower
+  const double factor =
+      floor((double)getGridPointsPerAxis()*(offset/length)
+          +std::numeric_limits<double>::epsilon()*1.e4);
+  return factor*(length/(double)getGridPointsPerAxis());
+}
+
+double SamplingGrid::getNearestHigherGridPoint(
+    const double value, const size_t axis) const
+{
+  const double length = end[axis] - begin[axis];
+  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
+    return end[axis];
+  if (value > end[axis])
+    return end[axis];
+  const double offset = value - begin[axis];
+  if (offset < 0)
+    return begin[axis];
+  // modify value a little if value actually has been on a grid point but
+  // perturbed by numerical rounding: here down, as we always go higher
+  const double factor =
+      ceil((double)getGridPointsPerAxis()*(offset/length)
+          -std::numeric_limits<double>::epsilon()*1.e4);
+  return factor*(length/(double)getGridPointsPerAxis());
+}
+
+void SamplingGrid::setWindowSize(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+  for (size_t index=0;index<3;++index) {
+    begin_window[index] = getNearestLowerGridPoint(_begin_window[index], index);
+    ASSERT( begin_window[index] >= begin[index],
+        "SamplingGrid::setWindowSize() - window starts earlier than domain on "
+        +toString(index)+"th component.");
+    end_window[index] = getNearestHigherGridPoint(_end_window[index], index);
+    ASSERT( end_window[index] <= end[index],
+        "SamplingGrid::setWindowSize() - window ends later than domain on "
+        +toString(index)+"th component.");
+  }
+}
+
+void SamplingGrid::setWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+  setWindowSize(_begin_window, _end_window);
+  const size_t gridpoints_window = getWindowGridPoints();
+  sampled_grid.clear();
+  sampled_grid.resize(gridpoints_window, 0.);
+}
+
+void SamplingGrid::setDomain(
+    const double _begin[3],
+    const double _end[3])
+{
+  setDomainSize(_begin, _end);
+  setWindowSize(_begin, _end);
+  const size_t gridpoints = getTotalGridPoints();
+  sampled_grid.resize(gridpoints, 0.);
+}
+
+void SamplingGrid::extendWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+#ifndef NDEBUG
+  for(size_t index=0;index < 3; ++index) {
+    // check that we truly have to extend the window
+    ASSERT ( begin_window[index] >= _begin_window[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window start is greater than old value.");
+    ASSERT ( end_window[index] <= _end_window[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window end is less than old value.");
+
+    // check that we are still less than domain
+    ASSERT ( _begin_window[index] >= begin[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window start is less than domain start.");
+    ASSERT ( _end_window[index] <= end[index],
+        "SamplingGrid::extendWindow() - component "+toString(index)+
+        " of window end is greater than domain end.");
+  }
+#endif
+  // copy old window size and values
+  double old_begin_window[3];
+  double old_end_window[3];
+  for(size_t index=0;index<3;++index) {
+    old_begin_window[index] = begin_window[index];
+    old_end_window[index] = end_window[index];
+  }
+  sampledvalues_t old_values(sampled_grid);
+  // set new window
+  setWindow(_begin_window,_end_window);
+  // now extend it ...
+  addOntoWindow(old_begin_window, old_end_window, old_values, +1.);
+  LOG(2, "DEBUG: Grid after extension is " << sampled_grid << ".");
+}
+
+void SamplingGrid::shrinkWindow(
+    const double _begin_window[3],
+    const double _end_window[3])
+{
+#ifndef NDEBUG
+  for(size_t index=0;index < 3; ++index) {
+    // check that we truly have to shrink the window
+    ASSERT ( begin_window[index] <= _begin_window[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window start is less than old value.");
+    ASSERT ( end_window[index] >= _end_window[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window end is greater than old value.");
+
+    // check that we are still less than domain
+    ASSERT ( _begin_window[index] >= begin[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window start is less than domain start.");
+    ASSERT ( _end_window[index] <= end[index],
+        "SamplingGrid::shrinkWindow() - component "+toString(index)+
+        " of window end is greater than domain end.");
+  }
+#endif
+  // copy old window size and values
+  double old_begin_window[3];
+  double old_end_window[3];
+  for(size_t index=0;index<3;++index) {
+    old_begin_window[index] = begin_window[index];
+    old_end_window[index] = end_window[index];
+  }
+  sampledvalues_t old_values(sampled_grid);
+  // set new window
+  setWindow(_begin_window,_end_window);
+  // now extend it ...
+  addIntoWindow(old_begin_window, old_end_window, old_values, +1.);
+  LOG(2, "DEBUG: Grid after extension is " << sampled_grid << ".");
+}
+
+static void addElements(
+    double &dest,
+    const double &source,
+    const double prefactor)
+{
+  dest += prefactor*(source);
+}
+
+void SamplingGrid::addOntoWindow(
+    const double _begin_window[3],
+    const double _end_window[3],
+    const sampledvalues_t &_sampled_grid,
+    const double prefactor)
+{
+  addWindowOntoWindow(
+      begin_window,
+      end_window,
+      _begin_window,
+      _end_window,
+      sampled_grid,
+      _sampled_grid,
+      boost::bind(addElements, _1, _2, boost::cref(prefactor)),
+      destwindow);
+}
+
+void SamplingGrid::addIntoWindow(
+    const double _begin_window[3],
+    const double _end_window[3],
+    const sampledvalues_t &_sampled_grid,
+    const double prefactor)
+{
+  addWindowOntoWindow(
+      _begin_window,
+      _end_window,
+      begin_window,
+      end_window,
+      sampled_grid,
+      _sampled_grid,
+      boost::bind(addElements, _1, _2, boost::cref(prefactor)),
+      sourcewindow);
+}
+
+void SamplingGrid::addWindowOntoWindow(
+    const double larger_wbegin[3],
+    const double larger_wend[3],
+    const double smaller_wbegin[3],
+    const double smaller_wend[3],
+    sampledvalues_t &dest_sampled_grid,
+    const sampledvalues_t &source_sampled_grid,
+    boost::function<void (double &, const double &)> op,
+    enum eLargerWindow larger_window)
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( smaller_wbegin[index] >= larger_wbegin[index],
+        "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
+        +toString(index)+".");
+    ASSERT( smaller_wend[index] <= larger_wend[index],
+        "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
+        +toString(index)+".");
+  }
+#endif
+  // the only issue are indices
+  const size_t gridpoints_axis = getGridPointsPerAxis();
+  size_t pre_offset[3];
+  size_t post_offset[3];
+  size_t length[3];
+  size_t total[3];
+  const double round_offset =
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.; // need offset to get to round_toward_nearest behavior
+  for(size_t index=0;index<3;++index) {
+    // we refrain from using floor/ceil as the window's starts and ends,
+    // the grids have to be compatible (equal level), should always be on
+    // discrete grid point locations.
+    const double delta = (double)gridpoints_axis/(end[index] - begin[index]);
+    // delta is conversion factor from box length to discrete length, i.e. number of points
+    pre_offset[index] = delta*(smaller_wbegin[index] - larger_wbegin[index])+round_offset;
+    length[index] = delta*(smaller_wend[index] - smaller_wbegin[index])+round_offset;
+    post_offset[index] = delta*(larger_wend[index] - smaller_wend[index])+round_offset;
+    total[index] = delta*(larger_wend[index] - larger_wbegin[index])+round_offset;
+    // total is used as safe-guard against loss due to discrete conversion
+    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
+        "SamplingGrid::addWindowOntoWindow() - pre, post, and length don't sum up to total for "
+        +toString(index)+"th component.");
+  }
+  // assert that calculated lengths match with given vector sizes
+#ifndef NDEBUG
+  const size_t calculated_size = length[0]*length[1]*length[2];
+  if (larger_window == destwindow) {
+    ASSERT( calculated_size == source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: "
+        +toString(calculated_size)+" != "+toString(source_sampled_grid.size())+".");
+    ASSERT( calculated_size <= dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough sampled values available: "
+        +toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+".");
+  } else {
+    ASSERT( calculated_size == dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: "
+        +toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+".");
+    ASSERT( calculated_size <= source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: "
+        +toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+".");
+  }
+  const size_t total_size = total[0]*total[1]*total[2];
+  if (larger_window == destwindow) {
+    ASSERT( total_size == dest_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: "
+        +toString(total_size)+" != "+toString(dest_sampled_grid.size())+".");
+  } else {
+    ASSERT( total_size == source_sampled_grid.size(),
+        "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: "
+        +toString(total_size)+" != "+toString(source_sampled_grid.size())+".");
+  }
+#endif
+  size_t N[3];
+//  size_t counter = 0;
+  sampledvalues_t::iterator destiter = dest_sampled_grid.begin();
+  sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin();
+  if (larger_window == destwindow)
+    std::advance(destiter, pre_offset[0]*total[1]*total[2]);
+  else
+    std::advance(sourceiter, pre_offset[0]*total[1]*total[2]);
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    if (larger_window == destwindow)
+      std::advance(destiter, pre_offset[1]*total[2]);
+    else
+      std::advance(sourceiter, pre_offset[1]*total[2]);
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      if (larger_window == destwindow)
+        std::advance(destiter, pre_offset[2]);
+      else
+        std::advance(sourceiter, pre_offset[2]);
+      for(N[2]=0; N[2] < length[2]; ++N[2]) {
+        ASSERT( destiter != dest_sampled_grid.end(),
+            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+        ASSERT( sourceiter != source_sampled_grid.end(),
+            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+        op(*destiter, *sourceiter);
+        ++destiter;
+        ++sourceiter;
+      }
+      if (larger_window == destwindow)
+        std::advance(destiter, post_offset[2]);
+      else
+        std::advance(sourceiter, post_offset[2]);
+    }
+    if (larger_window == destwindow)
+      std::advance(destiter, post_offset[1]*total[2]);
+    else
+      std::advance(sourceiter, post_offset[1]*total[2]);
+  }
+#ifndef NDEBUG
+  if (larger_window == destwindow)
+    std::advance(destiter, post_offset[0]*total[1]*total[2]);
+  else
+    std::advance(sourceiter, post_offset[0]*total[1]*total[2]);
+  ASSERT( destiter == dest_sampled_grid.end(),
+      "SamplingGrid::addWindowOntoWindow() - destiter is not at end of window.");
+  ASSERT( sourceiter == source_sampled_grid.end(),
+      "SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window.");
+#endif
+  LOG(2, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
+}
+
+std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
+{
+  ost << "SamplingGrid";
+  ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
+  ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
+  ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2];
+  ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2];
+  ost << ", level of " << other.level;
+  ost << " and integrated value of " << other.integral();
+  return ost;
+}
+
+template<> SamplingGrid ZeroInstance<SamplingGrid>()
+{
+  SamplingGrid returnvalue;
+  return returnvalue;
+}
+
+// we need to explicitly instantiate the serialization functions as
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)
