source: src/LinkedCell/LinkedCell_Controller.cpp@ 8c31865

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2011 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * 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 "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Range.hpp"
#include "LinkedCell_Controller.hpp"
#include "LinkedCell_Model.hpp"
#include "LinkedCell_View.hpp"
#include "IPointCloud.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) :
    domain(_domain)
{
  /// 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()
{
  /// 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;
}

/** 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, 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)
{
  /// Look for best instance.
  const LinkedCell_Model *LCModel_best = getBestModel(distance);

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