source: src/LinkedCell/LinkedCell_Controller.cpp@ 4a187d

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

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

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

#include "CodePatterns/MemDebug.hpp"

#include <set>

#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Observer/Notification.hpp"
#include "CodePatterns/Range.hpp"
#include "LinkedCell_Controller.hpp"
#include "LinkedCell_Model.hpp"
#include "LinkedCell_View.hpp"
#include "LinkedCell_View_ModelWrapper.hpp"
#include "IPointCloud.hpp"
#include "WorldTime.hpp"


using namespace LinkedCell;

double LinkedCell_Controller::lower_threshold = 1.;
double LinkedCell_Controller::upper_threshold = 20.;

/** Constructor of class LinkedCell_Controller.
 *
 */
LinkedCell_Controller::LinkedCell_Controller(const Box &_domain) :
    Observer("LinkedCell_Controller"),
    domain(_domain)
{
  // sign on to specific notifications
  domain.signOn(this, Box::MatrixChanged);
  WorldTime::getInstance().signOn(this, WorldTime::TimeChanged);

  /// Check that upper_threshold fits within half the box.
  Vector diagonal(1.,1.,1.);
  diagonal.Scale(upper_threshold);
  Vector diagonal_transformed = domain.getMinv() * diagonal;
  double max_factor = 1.;
  for (size_t i=0; i<NDIM; ++i)
    if (diagonal_transformed.at(i) > 1./max_factor)
      max_factor = 1./diagonal_transformed.at(i);
  upper_threshold *= max_factor;

  /// Check that lower_threshold is still lower, if not set to half times upper_threshold.
  if (lower_threshold > upper_threshold)
    lower_threshold = 0.5*upper_threshold;
}

/** Destructor of class LinkedCell_Controller.
 *
 * Here, we free all LinkedCell_Model instances again.
 *
 */
LinkedCell_Controller::~LinkedCell_Controller()
{
  // sign off
  domain.signOff(this, Box::MatrixChanged);
  WorldTime::getInstance().signOff(this, WorldTime::TimeChanged);

  /// we free all LinkedCell_Model instances again.
  for(MapEdgelengthModel::iterator iter = ModelsMap.begin();
      !ModelsMap.empty(); iter = ModelsMap.begin()) {
    delete iter->second;
    ModelsMap.erase(iter);
  }
}

/** Internal function to obtain the range within which an model is suitable.
 *
 * \note We use statics lower_threshold and upper_threshold as min and max
 * boundaries.
 *
 * @param distance desired egde length
 * @return range within which model edge length is acceptable
 */
const range<double> LinkedCell_Controller::getHeuristicRange(const double distance) const
{
  const double lower = 0.5*distance < lower_threshold ? lower_threshold : 0.5*distance;
  const double upper = 2.*distance > upper_threshold ? upper_threshold : 2.*distance;
  range<double> HeuristicInterval(lower, upper);
  return HeuristicInterval;
}

static size_t checkLengthsLess(
    const Vector &Lengths,
    const double min_distance)
{
  size_t counter=0;
  for (size_t i=0;i<NDIM;++i) {
    if (Lengths[i] < min_distance)
      ++counter;
  }
  return counter;
}

/** Given some given \a distance return it such that it is larger equal than some minimum
 * sensible value with respect to the current domain.
 *
 * The sensible minimum distance \f$\Delta\f$ is computed as follows:
 * \f$ \frac {l_1} \Delta \frac {l_2} \Delta \frac {l_3} \Delta \leq M \f$, where M is the
 * maximum number of linked cells. Where some additional thought must be given
 * to rounding issues. Hence, we calculate \f$ (\frac {l_i} d)^{\frac 1 d}+.5\f$,
 * where d is the assumed dimensionality of the domain: d=3 for cubic, d=2 for
 * a planar, d=1 for a bar-like system.
 *
 * \param distance distance compared with minimum value
 * \return minimum sensible value or \a distance if larger
 */
double LinkedCell_Controller::getMinimumSensibleLength(double distance) const
{
  // first check the box with respect to the given distance
  // as we must not generate more than 10^3 cells per axis due to
  // memory constraints
  const RealSpaceMatrix &M = domain.getM();
  RealSpaceMatrix UnitMatrix;
  UnitMatrix.setIdentity();
  UnitMatrix *= M;
  Vector Lengths = UnitMatrix.transformToEigenbasis();
  double max_length = 0.;
  double volume = 1.;
  double min_distance = lower_threshold;
  for (size_t dim=NDIM; dim>0; --dim) {
    const double fraction = pow(MAX_LINKEDCELLNODES+1., 1./(double)dim);
    for (size_t i=0;i<NDIM;++i) {
      max_length = std::max(Lengths[i], max_length);
      // the length dropping out for dimensionality less than 3 must appear
      // just as factor of 1 not less.
      if (Lengths[i] < min_distance)
        volume *= pow(min_distance, 1./(double)dim)+.5;
      else
        volume *= pow(Lengths[i], 1./(double)dim)+.5; //.5 for rounding issues
    }
    min_distance = volume/fraction;
    // assuming a (cubic-like) 3d system
    if (checkLengthsLess(Lengths, min_distance) == 0)
      break;
  }
  if (distance < min_distance) {
    LOG(1, "INFO: Setting LC distance from given " << distance
        << " to " << min_distance << " due to box size.");
    distance = min_distance;
  }
  return distance;
}

/** Internal function to decide whether a suitable model is present or not.
 *
 * Here, the heuristic for deciding whether a new linked cell structure has to
 * be constructed or not is implemented. The current heuristic is as follows:
 * -# the best model should have at least half the desired length (such
 *    that most we have to look two neighbor shells wide and not one).
 * -# the best model should have at most twice the desired length but
 *    no less than 1 angstroem.
 *
 * \note Dealing out a pointer is here (hopefully) safe because the function is
 * internal and we - inside this class - know what we are doing.
 *
 * @param distance edge length of the requested linked cell structure
 * @return NULL - there is no fitting LinkedCell_Model, else - pointer to instance
 */
const LinkedCell_Model *LinkedCell_Controller::getBestModel(double distance) const
{
  /// Bound distance to be within [lower_threshold, upper_threshold).
  /// Note that we need to stay away from upper boundary a bit,
  /// otherwise the distance will end up outside of the interval.
  if (distance < lower_threshold)
    distance = lower_threshold;
  if (distance > upper_threshold)
    distance = upper_threshold - std::numeric_limits<double>::round_error();

  /// Look for all models within [0.5 distance, 2. distance).
  MapEdgelengthModel::const_iterator bestmatch = ModelsMap.end();
  if (!ModelsMap.empty()) {
    for(MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
        iter != ModelsMap.end(); ++iter) {
      // check that we are truely within range
      range<double> HeuristicInterval(getHeuristicRange(iter->first));
      if (HeuristicInterval.isInRange(distance)) {
        // if it's the first match or a closer one, pick
        if ((bestmatch == ModelsMap.end())
            || (fabs(bestmatch->first - distance) > fabs(iter->first - distance)))
          bestmatch = iter;
      }
    }
  }

  /// Return best match or NULL if none found.
  if (bestmatch != ModelsMap.end())
    return bestmatch->second;
  else
    return NULL;
}

/** Internal function to insert a new model and check for valid insertion.
 *
 * @param distance edge length of new model
 * @param instance pointer to model
 */
void LinkedCell_Controller::insertNewModel(const double edgelength, const LinkedCell_Model* instance)
{
  std::pair< MapEdgelengthModel::iterator, bool> inserter =
      ModelsMap.insert( std::make_pair(edgelength, instance) );
  ASSERT(inserter.second,
      "LinkedCell_Controller::getView() - LinkedCell_Model instance with distance "
      +toString(edgelength)+" already present.");
}

/** Returns the a suitable LinkedCell_Model contained in a LinkedCell_View
 *  for the requested \a distance.
 *
 * \sa getBestModel()
 *
 * @param distance edge length of the requested linked cell structure
 * @param set of initial points to insert when new model is created (not always), should be World's
 * @return LinkedCell_View wrapping the best LinkedCell_Model
 */
LinkedCell_View LinkedCell_Controller::getView(const double distance, IPointCloud &set)
{
  // distance should be a given minimum length dependent on domain at least
  const double ConstraintDistance = getMinimumSensibleLength(distance);

  /// Look for best instance.
  const LinkedCell_Model * const LCModel_best = getBestModel(ConstraintDistance);

  /// Construct new instance if none found,
  if (LCModel_best == NULL) {
    LinkedCell_Model * const LCModel_new = new LinkedCell_Model(ConstraintDistance, domain);
    LCModel_new->insertPointCloud(set);
    insertNewModel(ConstraintDistance, LCModel_new);
    LinkedCell_View view(*LCModel_new);
    return view;
  } else {
    /// else construct interface and return.
    LinkedCell_View view(*LCModel_best);
    return view;
  }
}

/** Internal function to re-create all present and used models for the new Box.
 *
 * This is necessary in the following cases:
 * -# the Box is changed
 * -# we step on to a different time step, i.e. all atomic positions change
 *
 * The main problem are the views currently in use.
 *
 * We make use of LinkedCell:LinkedCell_View::RAIIMap as there all present are
 * listed. We go through the list, create a map with old model ref as keys to
 * just newly created ones, and finally go again through each view and exchange
 * the model against the new ones via a simple map lookup.
 *
 */
void LinkedCell_Controller::updateModels()
{
  LOG(1, "INFO: Updating all models.");

  typedef std::map<const LinkedCell_Model *, LinkedCell_Model *>  ModelLookup;
  ModelLookup models;

  // set up map, for now with NULL pointers
  for (LinkedCell_View::ModelInstanceMap::const_iterator iter = LinkedCell_View::RAIIMap.begin();
      iter != LinkedCell_View::RAIIMap.end(); ++iter) {
#ifndef NDEBUG
    std::pair< ModelLookup::iterator, bool > inserter =
#endif
        models.insert( std::pair<const LinkedCell_Model *, LinkedCell_Model *>( (*iter)->LC->getModel(), NULL) );
    LOG(2, "INFO: Added " << (*iter)->LC->getModel() << " to list of models to replace.");
    ASSERT( inserter.second,
        "LinkedCell_Controller::updateModelsForNewBoxMatrix() - failed to insert old model "
        +toString( (*iter)->LC->getModel() )+",NULL into models, is already present");
  }

  // invert MapEdgelengthModel
  LOG(2, "INFO: ModelsMap is " << ModelsMap << ".");
  typedef std::map<const LinkedCell_Model*, double > MapEdgelengthModel_inverted;
  MapEdgelengthModel_inverted ModelsMap_inverted;
  for (MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
      iter != ModelsMap.end(); ++iter) {
#ifndef NDEBUG
    MapEdgelengthModel_inverted::const_iterator assertiter = ModelsMap_inverted.find(iter->second);
    ASSERT( assertiter == ModelsMap_inverted.end(),
        "LinkedCell_Controller::updateModelsForNewBoxMatrix() - ModelsMap is not invertible, value "
        +toString(iter->second)+" is already present.");
#endif
    ModelsMap_inverted.insert( std::make_pair(iter->second, iter->first) );
  }
  LOG(2, "INFO: Inverted ModelsMap is " << ModelsMap_inverted << ".");

  // go through map and re-create models
  for (ModelLookup::iterator iter = models.begin(); iter != models.end(); ++iter) {
    // delete old model
    const LinkedCell_Model * const oldref = iter->first;
#ifndef NDEBUG
    MapEdgelengthModel_inverted::const_iterator assertiter = ModelsMap_inverted.find(oldref);
    ASSERT( assertiter != ModelsMap_inverted.end(),
        "LinkedCell_Controller::updateModelsForNewBoxMatrix() - ModelsMap_inverted does not contain old model "
        +toString(oldref)+".");
#endif
    const double distance = ModelsMap_inverted[oldref];
    // create new one, afterwards erase old model (this is for unit test to get different memory addresses)
    LinkedCell_Model * const newref = new LinkedCell_Model(distance, domain);
    MapEdgelengthModel::iterator oldmodeliter = ModelsMap.find(distance);
    delete oldmodeliter->second;
    ModelsMap.erase(oldmodeliter);
    LOG(2, "INFO: oldref is " << oldref << ", newref is " << newref << ".");
    iter->second = newref;
    // replace in ModelsMap
#ifndef NDEBUG
    std::pair< MapEdgelengthModel::iterator, bool > inserter =
#endif
        ModelsMap.insert( std::make_pair(distance, newref) );
    ASSERT( inserter.second,
        "LinkedCell_Controller::updateModelsForNewBoxMatrix() - failed to insert new model "
        +toString(distance)+","+toString(newref)+" into ModelsMap, is already present");
  }
  
  // remove all remaining active models (also those that don't have an active View on them)
  for (MapEdgelengthModel::iterator iter = ModelsMap.begin();
      !ModelsMap.empty();
      iter = ModelsMap.begin()) {
    delete iter->second;
    ModelsMap.erase(iter);
  }


  // delete inverted map for safety (values are gone)
  ModelsMap_inverted.clear();

  // go through views and exchange the models
  for (LinkedCell_View::ModelInstanceMap::const_iterator iter = LinkedCell_View::RAIIMap.begin();
      iter != LinkedCell_View::RAIIMap.end(); ++iter) {
    ModelLookup::const_iterator modeliter = models.find((*iter)->LC->getModel());
    ASSERT( modeliter != models.end(),
        "LinkedCell_Controller::updateModelsForNewBoxMatrix() - we miss a model "
        +toString((*iter)->LC->getModel())+" in ModelLookup.");
    // this is ugly but the only place where we have to set ourselves over the constness of the member variable
    if (modeliter != models.end()) {
      LOG(2, "INFO: Setting model to " << modeliter->second << " in view " << *iter << ".");
      (*iter)->LC->setModel(modeliter->second);
    }
  }
}

/** Callback function for Observer mechanism.
 *
 * @param publisher reference to the Observable that calls
 */
void LinkedCell_Controller::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_Controller::recieveNotification(Observable *publisher, Notification_ptr notification)
{
  if (publisher == &domain) {
    switch(notification->getChannelNo()) {
      case Box::MatrixChanged:
        LOG(1, "INFO: LinkedCell_Controller got update from Box.");
        updateModels();
        break;
      default:
        ASSERT(0,
            "LinkedCell_Controller::recieveNotification() - unwanted notification from Box "
            +toString(notification->getChannelNo())+" received.");
        break;
    }
  } else if (publisher == WorldTime::getPointer()) {
    switch(notification->getChannelNo()) {
      case WorldTime::TimeChanged:
        LOG(1, "INFO: LinkedCell_Controller got update from WorldTime.");
        updateModels();
        break;
      default:
        ASSERT(0,
            "LinkedCell_Controller::recieveNotification() - unwanted notification from WorldTime "
            +toString(notification->getChannelNo())+" received.");
        break;
    }
  } else {
    ELOG(1, "Notification " << notification->getChannelNo()
        << " from unknown publisher " << publisher << ".");
  }
}

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