source: src/Jobs/Grid/SamplingGrid.cpp@ b123a5

/*
 * 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 "Jobs/Grid/SamplingGrid.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()
{}

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

SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
{
  // check that grids are compatible
  if (isCompatible(other)) {
    sampledvalues_t::iterator iter = sampled_grid.begin();
    sampledvalues_t::const_iterator otheriter = other.sampled_grid.begin();
    for(; iter != sampled_grid.end(); ++iter, ++otheriter )
      *iter *= (*otheriter);
  } else {
    ASSERT(0, "SamplingGrid::operator*=() - superposing incompatible 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 double SamplingGrid::getVolume() const
{
  double volume = 1.;
  for (size_t i=0;i<3;++i)
    volume *= end[i]-begin[i];
  return volume;
}

const double SamplingGrid::getWindowVolume() const
{
  double volume = 1.;
  for (size_t i=0;i<3;++i)
    volume *= end_window[i]-begin_window[i];
  return volume;
}

const size_t SamplingGrid::getTotalGridPoints() const
{
  return pow(getGridPointsPerAxis(),3);
}

const size_t SamplingGrid::getGridPointsPerAxis() const
{
  return pow(2, level);
}

const size_t SamplingGrid::getWindowGridPoints() const
{
  const double volume = getVolume();
  if (fabs(volume) < std::numeric_limits<double>::epsilon())
    return 0;
  else
    return (getWindowVolume()*getTotalGridPoints())/volume;
}

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 0.;
  } else {
    const double offset = value - begin[axis];
    const double factor = floor((double)getGridPointsPerAxis()*(offset/length));
    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 0.;
  } else {
    const double offset = value - begin[axis];
    const double factor = ceil((double)getGridPointsPerAxis()*(offset/length));
    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::setDomainSize(
    const double _begin[3],
    const double _end[3])
{
  for (size_t index=0;index<3;++index) {
    begin[index] = _begin[index];
    end[index] = _end[index];
  }
}

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::addOntoWindow(
    const double _begin_window[3],
    const double _end_window[3],
    const sampledvalues_t &_sampled_grid,
    const double prefactor)
{
#ifndef NDEBUG
  for(size_t index=0;index<3;++index) {
    ASSERT( _begin_window[index] >= begin_window[index],
        "SamplingGrid::addOntoWindow() - given window starts earlier in component "
        +toString(index)+".");
    ASSERT( _end_window[index] <= end_window[index],
        "SamplingGrid::addOntoWindow() - given window ends later 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*(_begin_window[index] - begin_window[index])+round_offset;
    length[index] = delta*(_end_window[index] - _begin_window[index])+round_offset;
    post_offset[index] = delta*(end_window[index] - _end_window[index])+round_offset;
    total[index] = delta*(end_window[index] - begin_window[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::addOntoWindow() - 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];
  ASSERT( calculated_size == _sampled_grid.size(),
      "SamplingGrid::addOntoWindow() - not enough sampled values given: "
      +toString(calculated_size)+" != "+toString(_sampled_grid.size())+".");
  ASSERT( calculated_size <=  sampled_grid.size(),
      "SamplingGrid::addOntoWindow() - not enough sampled values available: "
      +toString(calculated_size)+" <= "+toString(sampled_grid.size())+".");
  const size_t total_size = total[0]*total[1]*total[2];
  ASSERT( total_size == sampled_grid.size(),
      "SamplingGrid::addOntoWindow() - total size is not equal to number of present points: "
      +toString(total_size)+" != "+toString(sampled_grid.size())+".");
#endif
  size_t N[3];
//  size_t counter = 0;
  sampledvalues_t::iterator griditer = sampled_grid.begin();
  std::advance(griditer, pre_offset[0]*total[1]*total[2]);
  sampledvalues_t::const_iterator copyiter = _sampled_grid.begin();
  for(N[0]=0; N[0] < length[0]; ++N[0]) {
    std::advance(griditer, pre_offset[1]*total[2]);
    for(N[1]=0; N[1] < length[1]; ++N[1]) {
      std::advance(griditer, pre_offset[2]);
      for(N[2]=0; N[2] < length[2]; ++N[2]) {
        ASSERT( griditer != sampled_grid.end(),
            "SamplingGrid::addOntoWindow() - griditer is already at end of window.");
        ASSERT( copyiter != _sampled_grid.end(),
            "SamplingGrid::addOntoWindow() - griditer is already at end of window.");
        *griditer++ += prefactor*(*copyiter++);
      }
      std::advance(griditer, post_offset[2]);
    }
    std::advance(griditer, post_offset[1]*total[2]);
  }
#ifndef NDEBUG
  std::advance(griditer, post_offset[0]*total[1]*total[2]);
  ASSERT( griditer == sampled_grid.end(),
      "SamplingGrid::addOntoWindow() - griditer is not at end of window.");
  ASSERT( copyiter == _sampled_grid.end(),
      "SamplingGrid::addOntoWindow() - copyiter is not at end of window.");
#endif
  LOG(2, "DEBUG: Grid after adding other is " << 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)