source: src/Fragmentation/Summation/SetValues/SamplingGrid.cpp@ 15c75f8

/*
 * 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/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(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.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(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.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 = getTotalLengthPerAxis(axis);
  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
    return begin[axis];
  if ((value - end[axis]) > -std::numeric_limits<double>::epsilon()*1.e4)
    return end[axis];
  const double offset = value - begin[axis];
  if (offset < std::numeric_limits<double>::epsilon()*1.e4)
    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 begin[axis]+factor*getDeltaPerAxis(axis);
}

double SamplingGrid::getNearestHigherGridPoint(
    const double value, const size_t axis) const
{
  const double length = getTotalLengthPerAxis(axis);
  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
    return end[axis];
  if ((value - end[axis]) > -std::numeric_limits<double>::epsilon()*1.e4)
    return end[axis];
  const double offset = value - begin[axis];
  if (offset < std::numeric_limits<double>::epsilon()*1.e4)
    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 begin[axis]+factor*getDeltaPerAxis(axis);
}

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(6, "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(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[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::getDiscreteWindowIndices(
    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 size_t gridpoints_axis = getGridPointsPerAxis();
  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) {
    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 = (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;
    } 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::getDiscreteWindowIndices() - pre, post, and length don't sum up to total for "
        +toString(index)+"th component.");
  }
}

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
  size_t pre_offset[3];
  size_t post_offset[3];
  size_t length[3];
  size_t total[3];
  getDiscreteWindowIndices(
      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[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(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<3; ++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;
}

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)