Changeset fe8253 for src/LinkedCell


Ignore:
Timestamp:
Jan 2, 2012, 1:00:05 PM (13 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
2614e2a
Parents:
9a5649
git-author:
Frederik Heber <heber@…> (11/30/11 10:26:35)
git-committer:
Frederik Heber <heber@…> (01/02/12 13:00:05)
Message:

Added first reasonable heuristic, and added required implementation for Controller.

  • LinkedCell_Controller::getBestModel() nows contains some more advanced heuristics.
  • extended unit test on all of the added parts.
Location:
src/LinkedCell
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/LinkedCell/LinkedCell_Controller.cpp

    r9a5649 rfe8253  
    2121
    2222#include "Box.hpp"
     23#include "CodePatterns/Assert.hpp"
     24#include "CodePatterns/Range.hpp"
    2325#include "LinkedCell_Controller.hpp"
    2426#include "LinkedCell_Model.hpp"
    2527#include "LinkedCell_View.hpp"
    2628
    27 #include "CodePatterns/Assert.hpp"
    2829
    2930using namespace LinkedCell;
     31
     32double LinkedCell_Controller::lower_threshold = 1.;
     33double LinkedCell_Controller::upper_threshold = 20.;
    3034
    3135/** Constructor of class LinkedCell_Controller.
    … …  
    3438LinkedCell_Controller::LinkedCell_Controller(const Box &_domain) :
    3539    domain(_domain)
    36 {}
     40{
     41  /// Check that upper_threshold fits within half the box.
     42  Vector diagonal(1.,1.,1.);
     43  diagonal.Scale(upper_threshold);
     44  Vector diagonal_transformed = domain.getMinv() * diagonal;
     45  double max_factor = 1.;
     46  for (size_t i=0; i<NDIM; ++i)
     47    if (diagonal_transformed.at(i) > 1./max_factor)
     48      max_factor = 1./diagonal_transformed.at(i);
     49  upper_threshold *= max_factor;
     50
     51  /// Check that lower_threshold is still lower, if not set to half times upper_threshold.
     52  if (lower_threshold > upper_threshold)
     53    lower_threshold = 0.5*upper_threshold;
     54}
    3755
    3856/** Destructor of class LinkedCell_Controller.
    … …  
    4361LinkedCell_Controller::~LinkedCell_Controller()
    4462{
    45   // remove all linked cell models
     63  /// we free all LinkedCell_Model instances again.
    4664  for(MapEdgelengthModel::iterator iter = ModelsMap.begin();
    4765      !ModelsMap.empty(); iter = ModelsMap.begin()) {
    … …  
    5169}
    5270
     71/** Internal function to obtain the range within which an model is suitable.
     72 *
     73 * \note We use statics lower_threshold and upper_threshold as min and max
     74 * boundaries.
     75 *
     76 * @param distance desired egde length
     77 * @return range within which model edge length is acceptable
     78 */
     79const range<double> LinkedCell_Controller::getHeuristicRange(const double distance) const
     80{
     81  const double lower = 0.5*distance < lower_threshold ? lower_threshold : 0.5*distance;
     82  const double upper = 2.*distance > upper_threshold ? upper_threshold : 2.*distance;
     83  range<double> HeuristicInterval(lower, upper);
     84  return HeuristicInterval;
     85}
     86
    5387/** Internal function to decide whether a suitable model is present or not.
    5488 *
    5589 * Here, the heuristic for deciding whether a new linked cell structure has to
    56  * be constructed or not is implemented.
     90 * be constructed or not is implemented. The current heuristic is as follows:
     91 * -# the best model should have at least half the desired length (such
     92 *    that most we have to look two neighbor shells wide and not one).
     93 * -# the best model should have at most twice the desired length but
     94 *    no less than 1 angstroem.
    5795 *
    5896 * \note Dealing out a pointer is here (hopefully) safe because the function is
    … …  
    62100 * @return NULL - there is no fitting LinkedCell_Model, else - pointer to instance
    63101 */
    64 const LinkedCell_Model *LinkedCell_Controller::getBestModel(const double distance) const
     102const LinkedCell_Model *LinkedCell_Controller::getBestModel(double distance) const
    65103{
    66   for(MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
    67       iter != ModelsMap.end(); ++iter) {
    68     // here, we use the heuristic
    69     if (true) {
    70       return iter->second;
     104  /// Bound distance to be within [lower_threshold, upper_threshold).
     105  /// Note that we need to stay away from upper boundary a bit,
     106  /// otherwise the distance will end up outside of the interval.
     107  if (distance < lower_threshold)
     108    distance = lower_threshold;
     109  if (distance > upper_threshold)
     110    distance = upper_threshold - std::numeric_limits<double>::round_error();
     111
     112  /// Look for all models within [0.5 distance, 2. distance).
     113  MapEdgelengthModel::const_iterator bestmatch = ModelsMap.end();
     114  if (!ModelsMap.empty()) {
     115    for(MapEdgelengthModel::const_iterator iter = ModelsMap.begin();
     116        iter != ModelsMap.end(); ++iter) {
     117      // check that we are truely within range
     118      range<double> HeuristicInterval(getHeuristicRange(iter->first));
     119      if (HeuristicInterval.isInRange(distance)) {
     120        // if it's the first match or a closer one, pick
     121        if ((bestmatch == ModelsMap.end())
     122            || (fabs(bestmatch->first - distance) > fabs(iter->first - distance)))
     123          bestmatch = iter;
     124      }
    71125    }
    72126  }
    73   return NULL;
     127
     128  /// Return best match or NULL if none found.
     129  if (bestmatch != ModelsMap.end())
     130    return bestmatch->second;
     131  else
     132    return NULL;
     133}
     134
     135/** Internal function to insert a new model and check for valid insertion.
     136 *
     137 * @param distance edge length of new model
     138 * @param instance pointer to model
     139 */
     140void LinkedCell_Controller::insertNewModel(const double edgelength, LinkedCell_Model*& instance)
     141{
     142  std::pair< MapEdgelengthModel::iterator, bool> inserter =
     143      ModelsMap.insert( std::make_pair(edgelength, instance) );
     144  ASSERT(inserter.second,
     145      "LinkedCell_Controller::getView() - LinkedCell_Model instance with distance "
     146      +toString(edgelength)+" already present.");
    74147}
    75148
    … …  
    84157LinkedCell_View LinkedCell_Controller::getView(const double distance)
    85158{
    86   // look for best instance
     159  /// Look for best instance.
    87160  const LinkedCell_Model *LCModel_best = getBestModel(distance);
    88161
    89   // construct new instance if none found
     162  /// Construct new instance if none found,
    90163  if (LCModel_best == NULL) {
    91164    LinkedCell_Model *LCModel_new = new LinkedCell_Model(distance, domain);
    92     std::pair< MapEdgelengthModel::iterator, bool> inserter =
    93         ModelsMap.insert( std::make_pair(distance, LCModel_new) );
    94     ASSERT(inserter.second,
    95         "LinkedCell_Controller::getView() - LinkedCell_Model instance with distance "
    96         +toString(distance)+" already present.");
    97     LinkedCell_View interface(const_cast<const LinkedCell_Model &>(*LCModel_new));
    98     return interface;
     165    insertNewModel(distance, LCModel_new);
     166    LinkedCell_View view(const_cast<const LinkedCell_Model &>(*LCModel_new));
     167    return view;
    99168  } else {
    100     // construct interface and return
    101     LinkedCell_View interface(*LCModel_best);
    102     return interface;
     169    /// else construct interface and return.
     170    LinkedCell_View view(*LCModel_best);
     171    return view;
    103172  }
    104173}

  • src/LinkedCell/LinkedCell_Controller.hpp

    r9a5649 rfe8253  
    1616#include <map>
    1717
     18#include "CodePatterns/Range.hpp"
    1819#include "LinkedCell/LinkedCell_View.hpp"
     20
     21class LinkedCell_ControllerTest;
    1922
    2023namespace LinkedCell {
    … …  
    3033   * low computational complexity.
    3134   *
    32    * Here, it contains multiple LinkedCell_Model instances, each with different
    33    * edge lengths of the inherent mesh. You request a LinkedCell instance from
    34    * the controller, but you basically just receive a const interface reference
     35   * Here, it contains multiple LinkedCell models, each with different edge
     36   * lengths of the inherent mesh. You request a LinkedCell instance from the
     37   * controller, but you basically just receive a const interface reference
    3538   * to a present LinkedCell instance.
    3639   *
    … …  
    3841  class LinkedCell_Controller
    3942  {
     43    //!> grant unit tests access to internal parts
     44    friend class ::LinkedCell_ControllerTest;
    4045  public:
    4146    LinkedCell_Controller(const Box &_domain);
    … …  
    4651  private:
    4752    const LinkedCell_Model * getBestModel(const double distance) const;
     53    const range<double> getHeuristicRange(const double distance) const;
     54    void insertNewModel(const double edgelength, LinkedCell_Model*& instance);
    4855
    4956  private:
    … …  
    5461    MapEdgelengthModel ModelsMap;
    5562
     63    //!> lower threshold below which no new instances are generated.
     64    static double lower_threshold;
     65
     66    //!> upper threshold above which no new instances are generated.
     67    static double upper_threshold;
     68
    5669    //!> internal reference to domain required for constructing LinkedCell_Model instances.
    5770    const Box &domain;

  • src/LinkedCell/unittests/LinkedCell_ControllerUnitTest.cpp

    r9a5649 rfe8253  
    2828#include "LinearAlgebra/RealSpaceMatrix.hpp"
    2929#include "LinkedCell/LinkedCell_Controller.hpp"
     30#include "LinkedCell/LinkedCell_View.hpp"
    3031#include "LinkedCell/unittests/defs.hpp"
    3132
    … …  
    6566}
    6667
    67 
    68 /** UnitTest for LinkedCell_Controller::getViewTest().
    69  */
    70 void LinkedCell_ControllerTest::getViewTest()
    71 {
    72   CPPUNIT_ASSERT_EQUAL(true, true);
    73 
    74 }
     68/** UnitTest for LinkedCell_Controller's lower and upper thresholds.
     69 */
     70void LinkedCell_ControllerTest::thresholdTest()
     71{
     72  /// re-create instances
     73  delete controller;
     74  delete domain;
     75
     76  /// create diag(..) matrix beyond upper_threshold
     77  const double old_threshold = controller->upper_threshold;
     78  controller->lower_threshold = old_threshold*0.9;
     79  RealSpaceMatrix BoxM;
     80  BoxM.setIdentity();
     81  BoxM *= controller->upper_threshold*.5;
     82
     83  /// create Box with this matrix
     84  domain = new Box(BoxM);
     85
     86  controller = new LinkedCell::LinkedCell_Controller(*domain);
     87
     88  /// check that thresholds have been adapted
     89  CPPUNIT_ASSERT( controller->upper_threshold != old_threshold );
     90  CPPUNIT_ASSERT( controller->lower_threshold != old_threshold*0.9 );
     91}
     92
     93/** UnitTest for LinkedCell_Controller::getHeuristicRange().
     94 */
     95void LinkedCell_ControllerTest::getHeuristicRangeTest()
     96{
     97  /// re-implementing function to check is nonsense here, instead try some
     98  /// hard-coded, working values;
     99  controller->lower_threshold = 1.;
     100  controller->upper_threshold = 20.;
     101  const double inbetween = 9.5; // half and twice is definitely within both thresholds.
     102
     103  /// check distance in between
     104  range<double> interval = controller->getHeuristicRange(inbetween);
     105  CPPUNIT_ASSERT ( interval.first != controller->lower_threshold );
     106  CPPUNIT_ASSERT ( interval.last != controller->upper_threshold );
     107}
     108
     109/** UnitTest for LinkedCell_Controller::getViewTest() for getting twice the same view.
     110 */
     111void LinkedCell_ControllerTest::getView_SameViewTest()
     112{
     113  /// obtain a view
     114  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     115  LinkedCell::LinkedCell_View view = controller->getView(2.);
     116  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     117
     118  {
     119    /// get same view again and check that now new instance appears
     120    LinkedCell::LinkedCell_View view_again = controller->getView(2.);
     121    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     122  }
     123}
     124
     125/** UnitTest for LinkedCell_Controller::getViewTest() for picking two different views.
     126 */
     127void LinkedCell_ControllerTest::getView_DifferentViewTest()
     128{
     129  /// obtain a view
     130  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     131  LinkedCell::LinkedCell_View view = controller->getView(2.);
     132  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     133
     134  {
     135    /// pick another view that is not close enough
     136    LinkedCell::LinkedCell_View view_other = controller->getView(5.);
     137    CPPUNIT_ASSERT_EQUAL( (size_t)2, controller->ModelsMap.size() );
     138  }
     139}
     140
     141/** UnitTest for LinkedCell_Controller::getViewTest() for picking further views in range of present one.
     142 */
     143void LinkedCell_ControllerTest::getView_InRangeViewTest()
     144{
     145  /// obtain a view
     146  const double edgelength = 2.;
     147  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     148  LinkedCell::LinkedCell_View view = controller->getView(edgelength);
     149  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     150
     151  /// pick views that are close enough
     152  range<double> interval = controller->getHeuristicRange(edgelength);
     153  {
     154    /// ... at half lower interval half
     155    LinkedCell::LinkedCell_View view_lowerhalf = controller->getView((edgelength + interval.first)/2.);
     156    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     157  }
     158  {
     159    /// ... at half upper interval half
     160    LinkedCell::LinkedCell_View view_upperhalf = controller->getView((interval.last + edgelength)/2.);
     161    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     162  }
     163  {
     164    /// ... close to lower boundary
     165    LinkedCell::LinkedCell_View view_closelower = controller->getView(interval.first + std::numeric_limits<double>::round_error());
     166    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     167  }
     168  {
     169    /// ... close to upper boundary
     170    LinkedCell::LinkedCell_View view_closerupper = controller->getView(interval.last - std::numeric_limits<double>::round_error());
     171    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     172  }
     173  {
     174    /// on lower boundary
     175    LinkedCell::LinkedCell_View view_onlower = controller->getView(interval.first);
     176    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     177  }
     178}
     179
     180/** UnitTest for LinkedCell_Controller::getViewTest() for picking further views outside range.
     181 */
     182void LinkedCell_ControllerTest::getView_OutOfRangeViewTest()
     183{
     184  /// Here we need half of the edge length to be greater than lower_threshold
     185  const double edgelength = 2.5;
     186  CPPUNIT_ASSERT( (edgelength/2.) > controller->lower_threshold );
     187  /// obtain a view
     188  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     189  LinkedCell::LinkedCell_View view = controller->getView(edgelength);
     190  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     191
     192  /// pick views that are not close enough and check for new instance
     193  range<double> interval = controller->getHeuristicRange(edgelength);
     194  {
     195    /// ... outside lower boundary
     196    LinkedCell::LinkedCell_View view_outsidelower = controller->getView(interval.first - std::numeric_limits<double>::round_error());
     197    CPPUNIT_ASSERT_EQUAL( (size_t)2, controller->ModelsMap.size() );
     198  }
     199  {
     200    /// ... on upper boundary
     201    LinkedCell::LinkedCell_View view_onupper = controller->getView(interval.last);
     202    CPPUNIT_ASSERT_EQUAL( (size_t)3, controller->ModelsMap.size() );
     203  }
     204}
     205
     206/** UnitTest for LinkedCell_Controller::getViewTest() for picking views beneath lower threshold.
     207 */
     208void LinkedCell_ControllerTest::getView_LowerThresholdViewTest()
     209{
     210  /// obtain a view
     211  const double edgelength = 1.9*controller->lower_threshold;
     212  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     213  LinkedCell::LinkedCell_View view = controller->getView(edgelength);
     214  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     215
     216  {
     217    /// get a view at threshold and check that no new instance has been created
     218    LinkedCell::LinkedCell_View view_onlower = controller->getView(controller->lower_threshold);
     219    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     220  }
     221  {
     222    /// pick a view below 1.
     223    LinkedCell::LinkedCell_View view_beneathlower = controller->getView(0.1*controller->lower_threshold);
     224    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     225  }
     226}
     227
     228/** UnitTest for LinkedCell_Controller::getViewTest() for picking views above upper threshold.
     229 */
     230void LinkedCell_ControllerTest::getView_UpperThresholdViewTest()
     231{
     232  /// obtain a view
     233  const double edgelength = controller->upper_threshold;
     234  CPPUNIT_ASSERT_EQUAL( (size_t)0, controller->ModelsMap.size() );
     235  LinkedCell::LinkedCell_View view = controller->getView(edgelength);
     236  CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     237
     238  {
     239    /// get a view beyond threshold and check that no new instance has been created
     240    LinkedCell::LinkedCell_View view_beyondupper = controller->getView(1.1*controller->upper_threshold);
     241    CPPUNIT_ASSERT_EQUAL( (size_t)1, controller->ModelsMap.size() );
     242  }
     243
     244  {
     245    /// pick a view below threshold and check for new instance (if we make it outside acceptable range)
     246    range<double> interval = controller->getHeuristicRange(edgelength);
     247    if ( !interval.isInRange(0.1*controller->upper_threshold) ) {
     248      LinkedCell::LinkedCell_View view_beneathupper = controller->getView(0.1*controller->upper_threshold);
     249      CPPUNIT_ASSERT_EQUAL( (size_t)2, controller->ModelsMap.size() );
     250    }
     251  }
     252}

  • src/LinkedCell/unittests/LinkedCell_ControllerUnitTest.hpp

    r9a5649 rfe8253  
    2828{
    2929    CPPUNIT_TEST_SUITE( LinkedCell_ControllerTest) ;
    30     CPPUNIT_TEST ( getViewTest );
     30    CPPUNIT_TEST ( thresholdTest );
     31    CPPUNIT_TEST ( getHeuristicRangeTest );
     32    CPPUNIT_TEST ( getView_SameViewTest );
     33    CPPUNIT_TEST ( getView_DifferentViewTest );
     34    CPPUNIT_TEST ( getView_InRangeViewTest );
     35    CPPUNIT_TEST ( getView_OutOfRangeViewTest );
     36    CPPUNIT_TEST ( getView_LowerThresholdViewTest );
     37    CPPUNIT_TEST ( getView_UpperThresholdViewTest );
    3138    CPPUNIT_TEST_SUITE_END();
    3239
    … …  
    3542      void tearDown();
    3643
    37       void getViewTest();
     44      void thresholdTest();
     45      void getHeuristicRangeTest();
     46      void getView_SameViewTest();
     47      void getView_DifferentViewTest();
     48      void getView_InRangeViewTest();
     49      void getView_OutOfRangeViewTest();
     50      void getView_LowerThresholdViewTest();
     51      void getView_UpperThresholdViewTest();
    3852
    3953private:
Note: See TracChangeset for help on using the changeset viewer.