source: src/LinkedCell/LinkedCell_Model.cpp@ 7d0a4e

/*
 * 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/>.
 */

/*
 * LinkedCell_Model.cpp
 *
 *  Created on: Nov 15, 2011
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "LinkedCell_Model.hpp"

#include <algorithm>
#include <boost/bind.hpp>
#include <boost/multi_array.hpp>
#include <limits>

#include "Atom/AtomObserver.hpp"
#include "Atom/TesselPoint.hpp"
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Observer/Observer.hpp"
#include "CodePatterns/Observer/Notification.hpp"
#include "CodePatterns/toString.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinkedCell/IPointCloud.hpp"
#include "LinkedCell/LinkedCell.hpp"
#include "LinkedCell/LinkedCell_Model_changeModel.hpp"
#include "LinkedCell/LinkedCell_Model_LinkedCellArrayCache.hpp"
#include "World.hpp"

// initialize static entities
LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::NearestNeighbors;


/** Constructor of LinkedCell_Model.
 *
 * @param radius desired maximum neighborhood distance
 * @param _domain Box instance with domain size and boundary conditions
 */
LinkedCell::LinkedCell_Model::LinkedCell_Model(const double radius, const Box &_domain) :
    ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
    Changes( new changeModel(radius) ),
    internal_Sizes(NULL),
    N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
    domain(_domain)
{
  // set default argument
  NearestNeighbors[0] = NearestNeighbors[1] = NearestNeighbors[2] = 1;

  // get the partition of the domain
  setPartition(radius);

  // allocate linked cell structure
  AllocateCells();

  // sign in to AtomObserver
  startListening();
}

/** Constructor of LinkedCell_Model.
 *
 * @oaram set set of points to place into the linked cell structure
 * @param radius desired maximum neighborhood distance
 * @param _domain Box instance with domain size and boundary conditions
 */
LinkedCell::LinkedCell_Model::LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain) :
    ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
    Changes( new changeModel(radius) ),
    internal_Sizes(NULL),
    N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
    domain(_domain)
{
  Info info(__func__);

  // get the partition of the domain
  setPartition(radius);

  // allocate linked cell structure
  AllocateCells();

  insertPointCloud(set);

  // sign in to AtomObserver
  startListening();
}

/** Destructor of class LinkedCell_Model.
 *
 */
LinkedCell::LinkedCell_Model::~LinkedCell_Model()
{
  // sign off from observables
  stopListening();

  // reset linked cell structure
  Reset();
  delete N;

  // delete change queue
  delete Changes;
}

/** Signs in to AtomObserver and World to known about all changes.
 *
 */
void LinkedCell::LinkedCell_Model::startListening()
{
  World::getInstance().signOn(this, World::AtomInserted);
  World::getInstance().signOn(this, World::AtomRemoved);
  AtomObserver::getInstance().signOn(this, AtomObservable::PositionChanged);
}

/** Signs off from AtomObserver and World.
 *
 */
void LinkedCell::LinkedCell_Model::stopListening()
{
  World::getInstance().signOff(this, World::AtomInserted);
  World::getInstance().signOff(this, World::AtomRemoved);
  AtomObserver::getInstance().signOff(this, AtomObservable::PositionChanged);
}


/** Allocates as much cells per axis as required by
 * LinkedCell_Model::BoxPartition.
 *
 */
void LinkedCell::LinkedCell_Model::AllocateCells()
{
  // resize array
  tripleIndex index;
  for (int i=0;i<NDIM;i++)
    index[i] = static_cast<LinkedCellArray::index>(Dimensions.at(i,i));
  N->setN().resize(index);
  ASSERT(getSize(0)*getSize(1)*getSize(2) < MAX_LINKEDCELLNODES,
      "LinkedCell_Model::AllocateCells() - Number linked of linked cell nodes exceeded hard-coded limit, use greater edge length!");
  LOG(1, "INFO: Allocating array ("
      << getSize(0) << ","
      << getSize(1) << ","
      << getSize(2) << ") for a new LinkedCell_Model.");

  // allocate LinkedCell instances
  for(index[0] = 0; index[0] != static_cast<LinkedCellArray::index>(Dimensions.at(0,0)); ++index[0]) {
    for(index[1] = 0; index[1] != static_cast<LinkedCellArray::index>(Dimensions.at(1,1)); ++index[1]) {
      for(index[2] = 0; index[2] != static_cast<LinkedCellArray::index>(Dimensions.at(2,2)); ++index[2]) {
        LOG(5, "INFO: Creating cell at " << index[0] << " " << index[1] << " " << index[2] << ".");
        N->setN()(index) = new LinkedCell(index);
      }
    }
  }
}

/** Frees all Linked Cell instances and sets array dimensions to (0,0,0).
 *
 */
void LinkedCell::LinkedCell_Model::Reset()
{
  // free all LinkedCell instances
  for(iterator3 iter3 = N->setN().begin(); iter3 != N->setN().end(); ++iter3) {
    for(iterator2 iter2 = (*iter3).begin(); iter2 != (*iter3).end(); ++iter2) {
      for(iterator1 iter1 = (*iter2).begin(); iter1 != (*iter2).end(); ++iter1) {
        delete *iter1;
      }
    }
  }
  // set dimensions to zero
  N->setN().resize(boost::extents[0][0][0]);
}

/** Inserts all points contained in \a set.
 *
 * @param set set with points to insert into linked cell structure
 */
void LinkedCell::LinkedCell_Model::insertPointCloud(IPointCloud &set)
{
  if (set.IsEmpty()) {
    ELOG(1, "set is NULL or contains no linked cell nodes!");
    return;
  }

  // put each atom into its respective cell
  set.GoToFirst();
  while (!set.IsEnd()) {
    TesselPoint *Walker = set.GetPoint();
    addNode(Walker);
    set.GoToNext();
  }
}

/** Calculates the required edge length for the given desired distance.
 *
 * We need to make some matrix transformations in order to obtain the required
 * edge lengths per axis. Goal is guarantee that whatever the shape of the
 * domain that always all points at least up to \a distance away are contained
 * in the nearest neighboring cells.
 *
 * @param distance distance of this linked cell array
 */
void LinkedCell::LinkedCell_Model::setPartition(double distance)
{
  // generate box matrix of desired edge length
  RealSpaceMatrix neighborhood;
  neighborhood.setIdentity();
  neighborhood *= distance;

  // obtain refs to both domain matrix transformations
  //const RealSpaceMatrix &M = domain.getM();
  const RealSpaceMatrix &Minv = domain.getMinv();

  RealSpaceMatrix output = Minv * neighborhood;
  std::cout << Minv << " * " << neighborhood << " = " << output << std::endl;

  Dimensions = output.invert();
  std::cout << "Dimensions are then " << Dimensions << std::endl;

  // now dimensions is floating-point, but we need it to be integer (for allocation)
  for (size_t col = 0; col < NDIM; ++col) {
    for (size_t row = 0; row < NDIM; ++row) {
      ASSERT(fabs(Dimensions.at(row,col) - Dimensions.at(col,row)) < 1.e+3*std::numeric_limits<double>::epsilon(),
          "LinkedCell_Model::setPartition() - Dimensions is not symmetric by "
          +toString(fabs(Dimensions.at(row,col) - Dimensions.at(col,row)))+ ".");
      if (col != row) {
        ASSERT(fabs(Dimensions.at(row,col)) < 1.e+3*std::numeric_limits<double>::epsilon(),
            "LinkedCell_Model::setPartition() - Dimensions is not diagonal by "
            +toString(fabs(Dimensions.at(row,col)))+".");
      } else {
        Dimensions.set(row,col, ceil(Dimensions.at(row,col)));
      }
    }
  }


  Partition = Minv*Dimensions; //

  std::cout << "Partition matrix is then " << Partition << std::endl;
}

/** Returns the number of required neighbor-shells to get all neighboring points
 * in the given \a distance.
 *
 * @param distance radius of neighborhood sphere
 * @return number of LinkedCell's per dimension to get all neighbors
 */
const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getStep(const double distance) const
{
  tripleIndex index;
  index[0] = index[1] = index[2] = 0;

  if (fabs(distance) < std::numeric_limits<double>::min())
    return index;
  // generate box matrix of desired edge length
  RealSpaceMatrix neighborhood;
  neighborhood.setIdentity();
  neighborhood *= distance;

  const RealSpaceMatrix output = Partition * neighborhood;

  //std::cout << "GetSteps: " << Partition << " * " << neighborhood << " = " << output << std::endl;

  const RealSpaceMatrix steps = output;
  for (size_t i =0; i<NDIM; ++i)
    index[i] = ceil(steps.at(i,i));
  LOG(2, "INFO: number of shells are ("+toString(index[0])+","+toString(index[1])+","+toString(index[2])+").");

  return index;
}

/** Calculates the index of the cell \a position would belong to.
 *
 * @param position position whose associated cell to calculate
 * @return index of the cell
 */
const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
{
  ASSERT(domain.isValid(position),
      "LinkedCell::LinkedCell_Model::getIndexToVector() - specified position "+toString(position)
      +"is not valid in the current domain.");
  tripleIndex index;
  Vector x(Partition*domain.enforceBoundaryConditions(position));
  LOG(2, "INFO: Transformed position is " << x << ".");
  for (int i=0;i<NDIM;i++) {
    index[i] = static_cast<LinkedCellArray::index>(floor(x[i]));
  }
  return index;
}

/** Adds an update to the list of lazy changes to add a node.
 *
 * @param Walker node to add
 */
void LinkedCell::LinkedCell_Model::addNode(const TesselPoint *Walker)
{
  LOG(2, "INFO: Requesting update to add node " << *Walker << ".");
  Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::addNode_internal, this, _1));
}

/** Adds an update to the list of lazy changes to add remove a node.
 *
 * We do nothing of Walker is not found
 *
 * @param Walker node to remove
 */
void LinkedCell::LinkedCell_Model::deleteNode(const TesselPoint *Walker)
{
  LOG(2, "INFO: Requesting update to delete node " << *Walker << ".");
  Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::deleteNode_internal, this, _1));
}

/** Adds an update to the list of lazy changes to move a node.
 *
 * @param Walker node who has moved.
 */
void LinkedCell::LinkedCell_Model::moveNode(const TesselPoint *Walker)
{
  LOG(2, "INFO: Requesting update to move node " << *Walker << " to position "
      << Walker->getPosition() << ".");
  Changes->addUpdate(Walker, 10, boost::bind(&LinkedCell_Model::moveNode_internal, this, _1));
}

/** Internal function to add a node to the linked cell structure
 *
 * @param Walker node to add
 */
void LinkedCell::LinkedCell_Model::addNode_internal(const TesselPoint *Walker)
{
  // find index
  tripleIndex index = getIndexToVector(Walker->getPosition());
  LOG(2, "INFO: " << *Walker << " goes into cell " << index[0] << ", " << index[1] << ", " << index[2] << ".");
  LOG(2, "INFO: Cell's indices are "
      << (N->getN())(index)->getIndex(0) << " "
      << (N->getN())(index)->getIndex(1) << " "
      << (N->getN())(index)->getIndex(2) << ".");
  // add to cell
  (N->setN())(index)->addPoint(Walker);
  // add to index with check for presence
  std::pair<MapPointToCell::iterator, bool> inserter = CellLookup.insert( std::make_pair(Walker, (N->setN())(index)) );
  ASSERT( inserter.second,
      "LinkedCell_Model::addNode() - Walker "
      +toString(*Walker)+" is already present in cell "
      +toString((inserter.first)->second->getIndex(0))+" "
      +toString((inserter.first)->second->getIndex(1))+" "
      +toString((inserter.first)->second->getIndex(2))+".");
}

/** Internal function to remove a node to the linked cell structure
 *
 * We do nothing of Walker is not found
 *
 * @param Walker node to remove
 */
void LinkedCell::LinkedCell_Model::deleteNode_internal(const TesselPoint *Walker)
{
  MapPointToCell::iterator iter = CellLookup.find(Walker);
  ASSERT(iter != CellLookup.end(),
      "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
  if (iter != CellLookup.end()) {
    // remove from lookup
    CellLookup.erase(iter);
    // remove from cell
    iter->second->deletePoint(Walker);
  }
}

/** Internal function to move node from current cell to another on position change.
 *
 * @param Walker node who has moved.
 */
void LinkedCell::LinkedCell_Model::moveNode_internal(const TesselPoint *Walker)
{
  MapPointToCell::iterator iter = CellLookup.find(Walker);
  ASSERT(iter != CellLookup.end(),
      "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
  if (iter != CellLookup.end()) {
    tripleIndex index = getIndexToVector(Walker->getPosition());
    if (index != iter->second->getIndices()) {
      // remove in old cell
      iter->second->deletePoint(Walker);
      // add to new cell
      N->setN()(index)->addPoint(Walker);
      // update lookup
      iter->second = N->setN()(index);
    }
  }
}

/** Checks whether cell indicated by \a relative relative to LinkedCell_Model::internal_index
 * is out of bounds.
 *
 * \note We do not check for boundary conditions of LinkedCell_Model::domain,
 * we only look at the array sizes
 *
 * @param relative index relative to LinkedCell_Model::internal_index.
 * @return true - relative index is still inside bounds, false - outside
 */
bool LinkedCell::LinkedCell_Model::checkArrayBounds(const tripleIndex &index) const
{
  bool status = true;
  for (size_t i=0;i<NDIM;++i) {
    status = status && (
        (index[i] >= 0) &&
        (index[i] < getSize(i))
        );
  }
  return status;
}

/** Corrects \a index according to boundary conditions of LinkedCell_Model::domain .
 *
 * @param index index to correct according to boundary conditions
 */
void LinkedCell::LinkedCell_Model::applyBoundaryConditions(tripleIndex &index) const
{
  for (size_t i=0;i<NDIM;++i) {
    switch (domain.getConditions()[i]) {
    case BoundaryConditions::Wrap:
      if ((index[i] < 0) || (index[i] >= getSize(i)))
        index[i] = (index[i] % getSize(i));
      break;
    case BoundaryConditions::Bounce:
      if (index[i] < 0)
        index[i] = 0;
      if (index[i] >= getSize(i))
        index[i] = getSize(i)-1;
      break;
    case BoundaryConditions::Ignore:
      if (index[i] < 0)
        index[i] = 0;
      if (index[i] >= getSize(i))
        index[i] = getSize(i)-1;
      break;
    default:
      ASSERT(0, "LinkedCell_Model::checkBounds() - unknown boundary conditions.");
      break;
    }
  }
}

/** Calculates the interval bounds of the linked cell grid.
 *
 *  The neighborhood bounds works as follows: We give the lower, left, front
 *  corner of the box and the number of boxes to go in each direction (i.e. the
 *  relative upper, right, behind corner). We assure that the first corner is
 *  within LinkedCell_Model::N, whether the relative second corner is within the
 *  domain must be assured via applying its boundary conditions, see
 *  LinkedCell_Model::applyBoundaryConditions()
 *
 * \note we check whether \a index is valid, i.e. within the range of LinkedCell_Model::N.
 *
 * \param index index to give relative bounds to
 * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
 * \return pair of tripleIndex indicating lower and upper bounds
 */
const LinkedCell::LinkedCell_Model::LinkedCellNeighborhoodBounds LinkedCell::LinkedCell_Model::getNeighborhoodBounds(
    const tripleIndex &index,
    const tripleIndex &step
    ) const
{
  LinkedCellNeighborhoodBounds neighbors;

  // calculate bounds
  for (size_t i=0;i<NDIM;++i) {
    ASSERT(index[i] >= 0,
        "LinkedCell_Model::getNeighborhoodBounds() - index "
        +toString(index)+" out of lower bounds.");
    ASSERT (index[i] < getSize(i),
        "LinkedCell_Model::getNeighborhoodBounds() - index "
        +toString(index)+" out of upper bounds.");
    switch (domain.getConditions()[i]) {
      case BoundaryConditions::Wrap:
        if ((index[i] - step[i]) < 0)
          neighbors.first[i] = getSize(i) + (index[i] - step[i]);
        else if ((index[i] - step[i]) >= getSize(i))
          neighbors.first[i] = (index[i] - step[i]) - getSize(i);
        else
          neighbors.first[i] = index[i] - step[i];
        neighbors.second[i] = 2*step[i]+1;
        break;
      case BoundaryConditions::Bounce:
        neighbors.second[i] = 2*step[i]+1;
        if (index[i] - step[i] >= 0) {
          neighbors.first[i] = index[i] - step[i];
        } else {
          neighbors.first[i] = 0;
          neighbors.second[i] = index[i] + step[i]+1;
        }
        if (index[i] + step[i] >= getSize(i)) {
          neighbors.second[i] = getSize(i) - (index[i] - step[i]);
        }
        break;
      case BoundaryConditions::Ignore:
        if (index[i] - step[i] < 0)
          neighbors.first[i] = 0;
        else
          neighbors.first[i] = index[i] - step[i];
        if ((neighbors.first[i] + 2*step[i]+1) >= getSize(i))
          neighbors.second[i] = getSize(i) - neighbors.first[i];
        else
          neighbors.second[i] = 2*step[i]+1;
        break;
      default:
        ASSERT(0, "LinkedCell_Model::getNeighborhoodBounds() - unknown boundary conditions.");
        break;
    }
    // check afterwards whether we now have correct
    ASSERT((neighbors.first[i] >= 0) && (neighbors.first[i] < getSize(i)),
        "LinkedCell_Model::getNeighborhoodBounds() - lower border "
        +toString(neighbors.first)+" is out of bounds.");
  }
  LOG(3, "INFO: Resulting neighborhood bounds are [" << neighbors.first << "]<->[" << neighbors.second << "].");

  return neighbors;
}

/** Returns a reference to the cell indicated by LinkedCell_Model::internal_index.
 *
 * \return LinkedCell ref to current cell
 */
const LinkedCell::LinkedCell& LinkedCell::LinkedCell_Model::getCell(const tripleIndex &index) const
{
  return *(N->getN()(index));
}


/** Returns size of array for given \a dim.
 *
 * @param dim desired dimension
 * @return size of array along dimension
 */
LinkedCell::LinkedCellArray::index LinkedCell::LinkedCell_Model::getSize(const size_t dim) const
{
  ASSERT((dim >= 0) && (dim < NDIM),
      "LinkedCell_Model::getSize() - dimension "
      +toString(dim)+" is out of bounds.");
  return N->getN().shape()[dim];
}

/** Callback function for Observer mechanism.
 *
 * @param publisher reference to the Observable that calls
 */
void LinkedCell::LinkedCell_Model::update(Observable *publisher)
{
  ELOG(2, "LinkedCell_Model received inconclusive general update from "
      << publisher << ".");
}

/** Callback function for the Notifications mechanism.
 *
 * @param publisher reference to the Observable that calls
 * @param notification specific notification as cause of the call
 */
void LinkedCell::LinkedCell_Model::recieveNotification(Observable *publisher, Notification_ptr notification)
{
  // either it's the World or from the atom (through relay) itself
  if (publisher == World::getPointer()) {
    switch(notification->getChannelNo()) {
      case World::AtomInserted:
        addNode(World::getInstance().lastChanged<atom>());
        break;
      case World::AtomRemoved:
        deleteNode(World::getInstance().lastChanged<atom>());
        break;
    }
  } else {
    switch(notification->getChannelNo()) {
      case AtomObservable::PositionChanged:
      {
        moveNode(dynamic_cast<const TesselPoint *>(publisher));
        break;
      }
      default:
        LOG(2, "LinkedCell_Model received unwanted notification from AtomObserver's channel "
            << notification->getChannelNo() << ".");
        break;
    }
  }
}

/** Callback function when an Observer dies.
 *
 * @param publisher reference to the Observable that calls
 */
void LinkedCell::LinkedCell_Model::subjectKilled(Observable *publisher)
{}