source: src/Fragmentation/Summation/SetValues/SamplingGrid.cpp@ eb40e3

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. 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/Summation/SetValues/SamplingGrid.hpp"

#include <boost/assign.hpp>
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>

#include <algorithm>
#include <limits>

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

#include "LinearAlgebra/Vector.hpp"

#define MYSAMPLINGGRIDEPSILON (1e2*std::numeric_limits<double>::epsilon())

// static instances
const double SamplingGrid::zeroOffset[NDIM] = { 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[NDIM],
    const double _end[NDIM],
    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[NDIM],
    const double _end[NDIM],
    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<NDIM; ++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 double _value)
{
  std::transform(
      sampled_grid.begin(), sampled_grid.end(),
      sampled_grid.begin(),
      boost::bind(std::multiplies<double>(), _1, _value));

  return *this;
}

static void getDownsampledGrid(
    const SamplingGrid &_reference_grid,
    const SamplingGrid &_grid,
    boost::shared_ptr<SamplingGrid> &_weight_downsampled)
{
  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 int surplus_level = _reference_grid.getSurplusLevel(_grid)+round_offset;
  // need to downsample first
  _weight_downsampled.reset(new SamplingGrid); // may use default cstor as downsample does all settings
  SamplingGrid::downsample(*_weight_downsampled, _grid, _reference_grid.level-surplus_level);
}

SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
{
  // check that grids are compatible
  ASSERT(isCompatible(other),
      "SamplingGrid::operator*=() - multiplying incomatible grids is so far not in the cards.");
  const SamplingGrid *other_grid = &other;
  boost::shared_ptr<SamplingGrid> other_downsampled;
  if (!isEquivalent(other)) {
    getDownsampledGrid(*this, other, other_downsampled);
    other_grid = other_downsampled.get();
  }
  /// get minimum of window
  double min_begin_window[NDIM];
  double min_end_window[NDIM];
  bool doShrink = false;
  for (size_t index=0; index<NDIM;++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(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
  LOG(4, "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_grid->begin_window,
      other_grid->end_window,
      begin_window,
      end_window,
      sampled_grid,
      other_grid->sampled_grid,
      boost::bind(multiplyElements, _1, _2, 1.),
      sourcewindow);

  return *this;
}

void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
{
  // check that grids are compatible
  ASSERT(isCompatible(other),
      "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
  const SamplingGrid *other_grid = &other;
  boost::shared_ptr<SamplingGrid> other_downsampled;
  if (!isEquivalent(other)) {
    getDownsampledGrid(*this, other, other_downsampled);
    other_grid = other_downsampled.get();
  }
  /// get maximum of window
  double max_begin_window[NDIM];
  double max_end_window[NDIM];
  bool doExtend = false;
  for (size_t index=0; index<NDIM;++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(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
  LOG(4, "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_grid->begin_window, other_grid->end_window, other_grid->sampled_grid, prefactor);
}

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
{
  // check that grids are compatible
  ASSERT(isCompatible(weight),
      "SamplingGrid::integral() - integrating with weights from incompatible grids is so far not in the cards.");
  const SamplingGrid *weight_grid = &weight;
  boost::shared_ptr<SamplingGrid> weight_downsampled;
  if (!isEquivalent(weight)) {
    getDownsampledGrid(*this, weight, weight_downsampled);
    weight_grid = weight_downsampled.get();
  }
  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_grid->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;
}

void SamplingGrid::setWindowSize(
    const double _begin_window[NDIM],
    const double _end_window[NDIM])
{
  for (size_t index=0;index<NDIM;++index) {
    begin_window[index] = getNearestLowerGridPoint(_begin_window[index], index);
    ASSERT( (begin_window[index] - begin[index]) > -MYSAMPLINGGRIDEPSILON,
        "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]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::setWindowSize() - window ends later than domain on "
        +toString(index)+"th component.");
  }
}

void SamplingGrid::setWindow(
    const double _begin_window[NDIM],
    const double _end_window[NDIM])
{
  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[NDIM],
    const double _end[NDIM])
{
  setDomainSize(_begin, _end);
  setWindowSize(_begin, _end);
  const size_t gridpoints = getTotalGridPoints();
  sampled_grid.resize(gridpoints, 0.);
}

void SamplingGrid::extendWindow(
    const double _begin_window[NDIM],
    const double _end_window[NDIM])
{
#ifndef NDEBUG
  for(size_t index=0;index < NDIM; ++index) {
    // check that we truly have to extend the window
    ASSERT ( (begin_window[index] - _begin_window[index]) > -MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::extendWindow() - component "+toString(index)+
        " of window start is greater than old value.");
    ASSERT ( (end_window[index] - _end_window[index]) < MYSAMPLINGGRIDEPSILON,
        "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]) > -MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::extendWindow() - component "+toString(index)+
        " of window start is less than domain start.");
    ASSERT ( (_end_window[index] - end[index]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::extendWindow() - component "+toString(index)+
        " of window end is greater than domain end.");
  }
#endif
  // copy old window size and values
  double old_begin_window[NDIM];
  double old_end_window[NDIM];
  for(size_t index=0;index<NDIM;++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(6, "DEBUG: Grid after extension is " << sampled_grid << ".");
}

void SamplingGrid::shrinkWindow(
    const double _begin_window[NDIM],
    const double _end_window[NDIM])
{
#ifndef NDEBUG
  for(size_t index=0;index < NDIM; ++index) {
    // check that we truly have to shrink the window
    ASSERT ( (begin_window[index] - _begin_window[index]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::shrinkWindow() - component "+toString(index)+
        " of window start is less than old value.");
    ASSERT ( (end_window[index] - _end_window[index]) > -MYSAMPLINGGRIDEPSILON,
        "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]) > -MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::shrinkWindow() - component "+toString(index)+
        " of window start is less than domain start.");
    ASSERT ( (_end_window[index] - end[index]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::shrinkWindow() - component "+toString(index)+
        " of window end is greater than domain end.");
  }
#endif
  // copy old window size and values
  double old_begin_window[NDIM];
  double old_end_window[NDIM];
  for(size_t index=0;index<NDIM;++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(6, "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[NDIM],
    const double _end_window[NDIM],
    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[NDIM],
    const double _end_window[NDIM],
    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::getDiscreteWindowIndices(
    size_t _wbegin[NDIM],
    size_t _wlength[NDIM],
    size_t _wend[NDIM]) const
{
  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<NDIM;++index) {
    if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
      // 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 = getDeltaPerAxis(index);
      // delta is conversion factor from box length to discrete length, i.e. number of points
      _wbegin[index] = (begin_window[index] - begin[index])/delta+round_offset;
      _wlength[index] = (end_window[index] - begin_window[index])/delta+round_offset;
      _wend[index] = (end_window[index] - begin[index])/delta+round_offset;
    } else {
      _wbegin[index] = 0;
      _wlength[index] = 0;
      _wend[index] = 0;
    }
    // total is used as safe-guard against loss due to discrete conversion
    ASSERT( fabs(_wend[index] - _wbegin[index] - _wlength[index]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::getDiscreteWindowCopyIndices() - end - begin is not equal to length for "
        +toString(index)+"th component.");
  }
}

void SamplingGrid::getDiscreteWindowOffsets(
    size_t _pre_offset[NDIM],
    size_t _post_offset[NDIM],
    size_t _length[NDIM],
    size_t _total[NDIM]) const
{
  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<NDIM;++index) {
    if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
      // 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 = getDeltaPerAxis(index);
      // delta is conversion factor from box length to discrete length, i.e. number of points
      _pre_offset[index] = (begin_window[index] - begin[index])/delta+round_offset;
      _post_offset[index] = (end[index] - end_window[index])/delta+round_offset;
      _length[index] = (end_window[index] - begin_window[index])/delta+round_offset;
      _total[index] = (end[index] - begin[index])/delta+round_offset;
    } else {
      _pre_offset[index] = 0;
      _post_offset[index] = 0;
      _length[index] = 0;
      _total[index] = 0;
    }
    // total is used as safe-guard against loss due to discrete conversion
    ASSERT( (_pre_offset[index] + _post_offset[index]) + _length[index] == _total[index],
        "SamplingGrid::getDiscreteWindowCopyIndices() - pre, length, post are not equal to total for "
        +toString(index)+"th component.");
  }
}

void SamplingGrid::getDiscreteWindowCopyIndices(
    const double *larger_wbegin,
    const double *larger_wend,
    const double *smaller_wbegin,
    const double *smaller_wend,
    size_t *pre_offset,
    size_t *post_offset,
    size_t *length,
    size_t *total) const
{
  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<NDIM;++index) {
    if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
      // 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 = getDeltaPerAxis(index);
      // delta is conversion factor from box length to discrete length, i.e. number of points
      pre_offset[index] = (smaller_wbegin[index] - larger_wbegin[index])/delta+round_offset;
      length[index] = (smaller_wend[index] - smaller_wbegin[index])/delta+round_offset;
      post_offset[index] = (larger_wend[index] - smaller_wend[index])/delta+round_offset;
      total[index] = (larger_wend[index] - larger_wbegin[index])/delta+round_offset;
    } else {
      pre_offset[index] = 0;
      length[index] = 0;
      post_offset[index] = 0;
      total[index] = 0;
    }
    // total is used as safe-guard against loss due to discrete conversion
    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
        "SamplingGrid::getDiscreteWindowCopyIndices() - pre, post, and length don't sum up to total for "
        +toString(index)+"th component.");
  }
}

void SamplingGrid::addWindowOntoWindow(
    const double larger_wbegin[NDIM],
    const double larger_wend[NDIM],
    const double smaller_wbegin[NDIM],
    const double smaller_wend[NDIM],
    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<NDIM;++index) {
    ASSERT( (smaller_wbegin[index] - larger_wbegin[index]) > -MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
        +toString(index)+".");
    ASSERT( (smaller_wend[index] - larger_wend[index]) < MYSAMPLINGGRIDEPSILON,
        "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
        +toString(index)+".");
  }
#endif
  // the only issue are indices
  size_t pre_offset[NDIM];
  size_t post_offset[NDIM];
  size_t length[NDIM];
  size_t total[NDIM];
  getDiscreteWindowCopyIndices(
      larger_wbegin, larger_wend,
      smaller_wbegin, smaller_wend,
      pre_offset,
      post_offset,
      length,
      total
      );
  // 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[NDIM];
//  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() - sourceiter 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(8, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
}


bool SamplingGrid::operator==(const SamplingGrid &other) const
{
  bool status =
      static_cast<const SamplingGridProperties &>(*this)
      == static_cast<const SamplingGridProperties &>(other);
  // compare general properties
  if (status) {
    // compare windows
    for (size_t i=0; i<NDIM; ++i) {
      status &= begin_window[i] == other.begin_window[i];
      status &= end_window[i] == other.end_window[i];
    }
    // compare grids
    if (status)
      status &= sampled_grid == other.sampled_grid;
  }
  return status;
}

/** Struct contains a single point with displacements from the
 * central point and the weight in the restriction.
 */
struct PointWeight_t {
  PointWeight_t(const int &d1, const int &d2, const int &d3, const double &_weight) :
    displacement(NDIM),
    weight(_weight)
  {
    displacement[0] = d1; displacement[1] = d2; displacement[2] = d3;
  }
  typedef std::vector<int>  displacement_t;
  displacement_t displacement;
  double weight;
};

static void getLengthsOfWindow(
    int _total[NDIM],
    const SamplingGrid &_grid)
{
        const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
        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
        for (size_t index=0; index<NDIM; ++index) {
          if (fabs(_grid.end[index] - _grid.begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
            const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
            _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
          } else
            _total[index] = 0;
          // we can only assert that its atmost the maximum number of grid points
          ASSERT (_total[index] <= ::pow(2, _grid.level),
              "SamplingGrid::downsample() - total "+toString(_total[index])
              +" is not equal or less than 2^level: "+toString(_grid.level));
        }
}

//!> stencil for full weight restriction, see vmg's stencils.hpp
static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
    boost::assign::list_of
    ( PointWeight_t( 0,  0,  0, 0.125) )
    ( PointWeight_t( 1,  0,  0, 0.0625) )
    ( PointWeight_t(-1,  0,  0, 0.0625) )
    ( PointWeight_t( 0,  1,  0, 0.0625) )
    ( PointWeight_t( 0, -1,  0, 0.0625) )
    ( PointWeight_t( 0,  0,  1, 0.0625) )
    ( PointWeight_t( 0,  0, -1, 0.0625) )
    ( PointWeight_t( 1,  1,  0, 0.03125) )
    ( PointWeight_t( 1, -1,  0, 0.03125) )
    ( PointWeight_t(-1,  1,  0, 0.03125) )
    ( PointWeight_t(-1, -1,  0, 0.03125) )
    ( PointWeight_t( 0,  1,  1, 0.03125) )
    ( PointWeight_t( 0,  1, -1, 0.03125) )
    ( PointWeight_t( 0, -1,  1, 0.03125) )
    ( PointWeight_t( 0, -1, -1, 0.03125) )
    ( PointWeight_t( 1,  0,  1, 0.03125) )
    ( PointWeight_t( 1,  0, -1, 0.03125) )
    ( PointWeight_t(-1,  0,  1, 0.03125) )
    ( PointWeight_t(-1,  0, -1, 0.03125) )
    ( PointWeight_t( 1,  1,  1, 0.015625) )
    ( PointWeight_t( 1,  1, -1, 0.015625) )
    ( PointWeight_t( 1, -1,  1, 0.015625) )
    ( PointWeight_t(-1,  1,  1, 0.015625) )
    ( PointWeight_t( 1, -1, -1, 0.015625) )
    ( PointWeight_t(-1,  1, -1, 0.015625) )
    ( PointWeight_t(-1, -1,  1, 0.015625) )
    ( PointWeight_t(-1, -1, -1, 0.015625) )
;

int getValidIndex(
    const PointWeight_t::displacement_t &_disp,
    const int N[NDIM],
    const int length[NDIM])
{
  int index = 0;
  // we simply truncate in case of out of bounds access
  if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
    index += _disp[2];
  if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
    index += _disp[1]*length[2];
  if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
    index += _disp[0]*length[1]*length[2];
  return index;
}

void restrictFullWeight(
    SamplingGrid::sampledvalues_t &_coarse_level,
    const int length_c[NDIM],
    const SamplingGrid::sampledvalues_t &_fine_level,
    const int length_f[NDIM])
{
  int N_c[NDIM];
  int N_f[NDIM];
  SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
  for(N_c[0]=0, N_f[0]=0; (N_c[0] < length_c[0]) && (N_f[0] < length_f[0]); ++N_c[0], N_f[0] +=2) {
    for(N_c[1]=0, N_f[1]=0; (N_c[1] < length_c[1]) && (N_f[1] < length_f[1]); ++N_c[1], N_f[1] +=2) {
      for(N_c[2]=0, N_f[2]=0; (N_c[2] < length_c[2]) && (N_f[2] < length_f[2]); ++N_c[2], N_f[2] +=2) {
        const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
        // go through stencil and add each point relative to displacement with weight
        for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
            weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
          const PointWeight_t::displacement_t disp = weightiter->displacement;
          const int index_disp = getValidIndex(disp, N_f, length_f);
          *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
        }
        ++coarseiter;
      }
      ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
          "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
          +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
    }
    ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
        "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
        +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
  }
  ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
      "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
      +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
  ASSERT( coarseiter ==  _coarse_level.end(),
      "restrictFullWeight() - coarseiter is not at end of values.");
}

void SamplingGrid::padWithZerosForEvenNumberedSamples()
{
  size_t wbegin_index[NDIM];
  size_t wend_index[NDIM];
  size_t wlength_index[NDIM];
  getDiscreteWindowIndices(wbegin_index, wlength_index, wend_index);

  // calculate new window (always extend it such that both indices are even)
  bool changed = false;
  size_t wnewbegin_index[NDIM];
  size_t wnewend_index[NDIM];
  for(size_t i=0;i<NDIM;++i) {
    LOG(2, "INFO: window[" << i << "] starts at " << wbegin_index[i] << " and ends at "
        << wend_index[i] << " with length " << wlength_index[i]);
    // We require begin and end of window on even indices (and inside grid).
    wnewbegin_index[i] = wbegin_index[i];
    if ((wnewbegin_index[i] % (size_t)2) != 0) {
      if (wnewbegin_index[i] > 0)
        --wnewbegin_index[i];
      else
        wnewbegin_index[i] = 0;
      changed = true;
    }
    wnewend_index[i] = wend_index[i];
    if ((wnewend_index[i] % (size_t)2) != 0) {
      if (wnewend_index[i] < getGridPointsPerAxis())
        ++wnewend_index[i];
      else
        wnewend_index[i] = getGridPointsPerAxis();
      changed = true;
    }
    ASSERT( (wbegin_index[i] >= 0) && (wend_index[i] <= getGridPointsPerAxis()),
        "SamplingGrid::padWithZerosForEvenNumberedSamples() - indices "
        +toString(wbegin_index[i])+" and "+toString(wend_index[i])+" larger than grid "
        +toString(getGridPointsPerAxis())+".");
  }
  if (changed) {
    double begin_newwindow[NDIM];
    double end_newwindow[NDIM];
    for(size_t i=0;i<NDIM;++i) {
      const double delta = getDeltaPerAxis(i);
      begin_newwindow[i] = begin[i]+delta*wnewbegin_index[i];
      end_newwindow[i] = begin[i]+delta*wnewend_index[i];
    }
    LOG(2, "INFO: Original window is " << Vector(begin_window) << " <-> "  << Vector(end_window));
    LOG(2, "INFO: Padded window is " << Vector(begin_newwindow) << " <-> "  << Vector(end_newwindow));
    // extend window
    extendWindow(begin_newwindow, end_newwindow);
  }
  ASSERT( getWindowGridPoints() % (size_t)8 == 0,
      "SamplingGrid::padWithZerosForEvenNumberedSamples() - new size "
      +toString(getWindowGridPoints())+" is still not divisible by 8.");
}

void SamplingGrid::downsample(
    SamplingGrid& instance,
    const SamplingGrid& other,
    const int _level)
{
  if (&instance != &other) {
    LOG(2, "INFO: other's window is " << Vector(other.begin_window)
        << " <-> " << Vector(other.end_window));
    // take over properties
    static_cast<SamplingGridProperties &>(instance) = other;
    instance.setWindowSize(other.begin_window, other.end_window);
    LOG(2, "INFO: Set instance's window to " << Vector(instance.begin_window)
        << " <-> " << Vector(instance.end_window));
    ASSERT( _level <= other.level,
        "SamplingGrid::downsample() - desired level "+toString(_level)
        +" larger than level "+toString(other.level)+" of the given values.");
    if (_level == other.level) {
      instance.sampled_grid = other.sampled_grid;
    } else {
      // if desired level is smaller we need to downsample
      // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
      // neighbor interpolation) and always one grid level at a time till we
      // have reached the desired one

      // we need to copy the other grid because we might need to pad it with zeros anyway
      SamplingGrid FinerGrid(other);
      int length_d[NDIM];
      int length_s[NDIM];
      for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
        // pad with zeros for even indices and get length of fine grid window
        FinerGrid.padWithZerosForEvenNumberedSamples();
        getLengthsOfWindow(length_s, FinerGrid);
        ASSERT( FinerGrid.getWindowGridPoints() % 8 == 0,
            "SamplingGrid::downsample() - at level "+toString( instance.level)
            +" given grid points "+toString(FinerGrid.getWindowGridPoints())+" are not even numbered per axis anymore.");
        // re-adjust the window (length), get the corresponding window length and downsample
        instance.setWindow(FinerGrid.begin_window, FinerGrid.end_window);
        getLengthsOfWindow(length_d, instance);
        ASSERT( instance.sampled_grid.size() == FinerGrid.getWindowGridPoints()/(size_t)8,
            "SamplingGrid::downsample() - at level "+toString( instance.level)
            +" points on coarser grid "+toString(instance.sampled_grid.size())
            +" and downsampled number on finer grid "
            +toString(FinerGrid.getWindowGridPoints()/(size_t)8)+" do not match.");
        restrictFullWeight(instance.sampled_grid, length_d, FinerGrid.sampled_grid, length_s);
        // then use as new finer grid for next downsampling (if it's not the last level)
        if (instance.level > _level)
          FinerGrid = instance;
      }
      // loop stops at _level-1
      instance.level = _level;

      // and finally, renormalize downsampled grid to old value
//      instance *= other.integral()/instance.integral();
    }
  }
}

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)