Changeset 653cea


Ignore:
Timestamp:
Sep 12, 2016, 11:48:36 PM (8 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, 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_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
0983e6
Parents:
450adf
git-author:
Frederik Heber <heber@…> (06/30/14 09:35:42)
git-committer:
Frederik Heber <heber@…> (09/12/16 23:48:36)
Message:

SphericalPointDistribution is now working with bond degree weights.

  • recurseMatching() now works on IndexTupleList_t.
  • also rewrote calculatePairwiseDistances() and calculateErrorOfMatching().
  • L1THRESHOLD in recurseMatching() moved over to class body.
  • increased verbosity level of ...Matching() functions by one, added note on eventually chosen matching and why.
  • we assert that bestL2 is not too large.
  • FIX: calculateErrorOfMatching() did not use absolute value of gap for L1 error.
  • TESTFIX: Using limited accuracy on point coordinates.
  • TESTS: Regresssion test FragmentMolecule-cycles working again.
Location:
src/Fragmentation/Exporters
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r450adf r653cea  
    4949#include <limits>
    5050#include <list>
     51#include <numeric>
    5152#include <vector>
    5253#include <map>
    … …  
    6061const double SphericalPointDistribution::SQRT_3(sqrt(3.0));
    6162const double SphericalPointDistribution::warn_amplitude = 1e-2;
     63const double SphericalPointDistribution::L1THRESHOLD = 1e-2;
     64const double SphericalPointDistribution::L2THRESHOLD = 2e-1;
    6265
    6366typedef std::vector<double> DistanceArray_t;
    … …  
    6568// class generator: taken from www.cplusplus.com example std::generate
    6669struct c_unique {
    67   int current;
     70  unsigned int current;
    6871  c_unique() {current=0;}
    69   int operator()() {return current++;}
     72  unsigned int operator()() {return current++;}
    7073} UniqueNumber;
    7174
    72 inline
    73 DistanceArray_t calculatePairwiseDistances(
    74     const std::vector<Vector> &_points,
    75     const SphericalPointDistribution::IndexList_t &_indices
    76     )
    77 {
    78   DistanceArray_t result;
    79   for (SphericalPointDistribution::IndexList_t::const_iterator firstiter = _indices.begin();
    80       firstiter != _indices.end(); ++firstiter) {
    81     for (SphericalPointDistribution::IndexList_t::const_iterator seconditer = firstiter;
    82         seconditer != _indices.end(); ++seconditer) {
    83       if (firstiter == seconditer)
    84         continue;
    85       const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
    86       result.push_back(distance);
    87     }
    88   }
    89   return result;
    90 }
     75struct c_unique_list {
     76  unsigned int current;
     77  c_unique_list() {current=0;}
     78  std::list<unsigned int> operator()() {return std::list<unsigned int>(1, current++);}
     79} UniqueNumberList;
    9180
    9281/** Calculate the center of a given set of points in \a _positions but only
    … …  
    10897
    10998  return Center;
     99}
     100
     101inline
     102DistanceArray_t calculatePairwiseDistances(
     103    const SphericalPointDistribution::VectorArray_t &_points,
     104    const SphericalPointDistribution::IndexTupleList_t &_indices
     105    )
     106{
     107  DistanceArray_t result;
     108  for (SphericalPointDistribution::IndexTupleList_t::const_iterator firstiter = _indices.begin();
     109      firstiter != _indices.end(); ++firstiter) {
     110
     111    // calculate first center from possible tuple of indices
     112    Vector FirstCenter;
     113    ASSERT(!firstiter->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
     114    if (firstiter->size() == 1) {
     115      FirstCenter = _points[*firstiter->begin()];
     116    } else {
     117      FirstCenter = calculateCenter( _points, *firstiter);
     118      if (!FirstCenter.IsZero())
     119        FirstCenter.Normalize();
     120    }
     121
     122    for (SphericalPointDistribution::IndexTupleList_t::const_iterator seconditer = firstiter;
     123        seconditer != _indices.end(); ++seconditer) {
     124      if (firstiter == seconditer)
     125        continue;
     126
     127      // calculate second center from possible tuple of indices
     128      Vector SecondCenter;
     129      ASSERT(!seconditer->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
     130      if (seconditer->size() == 1) {
     131        SecondCenter = _points[*seconditer->begin()];
     132      } else {
     133        SecondCenter = calculateCenter( _points, *seconditer);
     134        if (!SecondCenter.IsZero())
     135          SecondCenter.Normalize();
     136      }
     137
     138      // calculate distance between both centers
     139      double distance = 2.; // greatest distance on surface of sphere with radius 1.
     140      if ((!FirstCenter.IsZero()) && (!SecondCenter.IsZero()))
     141        distance = (FirstCenter - SecondCenter).NormSquared();
     142      result.push_back(distance);
     143    }
     144  }
     145  return result;
    110146}
    111147
    … …  
    124160  Vector dreiBein(_oldPosition);
    125161  dreiBein.VectorProduct(_RotationAxis);
     162  ASSERT( !dreiBein.IsZero(), "determineSignOfRotation() - dreiBein is zero.");
    126163  dreiBein.Normalize();
    127164  const double sign =
    128165      (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
    129166  LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
    130       << ", newCenter in plane is " << _newPosition
     167      << ", newCenter on plane is " << _newPosition
    131168      << ", and dreiBein is " << dreiBein);
    132169  return sign;
    … …  
    167204 */
    168205std::pair<double, double> SphericalPointDistribution::calculateErrorOfMatching(
    169     const std::vector<Vector> &_old,
    170     const std::vector<Vector> &_new,
    171     const IndexList_t &_Matching)
     206    const VectorArray_t &_old,
     207    const VectorArray_t &_new,
     208    const IndexTupleList_t &_Matching)
    172209{
    173210  std::pair<double, double> errors( std::make_pair( 0., 0. ) );
    174211
    175212  if (_Matching.size() > 1) {
    176     LOG(3, "INFO: Matching is " << _Matching);
     213    LOG(5, "INFO: Matching is " << _Matching);
    177214
    178215    // calculate all pair-wise distances
    179     IndexList_t keys(_Matching.size());
    180     std::generate (keys.begin(), keys.end(), UniqueNumber);
     216    IndexTupleList_t keys(_old.size(), IndexList_t() );
     217    std::generate (keys.begin(), keys.end(), UniqueNumberList);
     218
    181219    const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
    182220    const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
    … …  
    188226    for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
    189227        ++firstiter, ++seconditer) {
    190       const double gap = *firstiter - *seconditer;
     228      const double gap = fabs(*firstiter - *seconditer);
    191229      // L1 error
    192230      if (errors.first < gap)
    … …  
    195233      errors.second += gap*gap;
    196234    }
    197   } else
    198     ELOG(3, "calculateErrorOfMatching() - Given matching's size is less than 2.");
    199   LOG(3, "INFO: Resulting errors for matching (L1, L2): "
     235  } else {
     236    // check whether we have any zero centers: Combining points on new sphere may result
     237    // in zero centers
     238    for (SphericalPointDistribution::IndexTupleList_t::const_iterator iter = _Matching.begin();
     239        iter != _Matching.end(); ++iter) {
     240      if ((iter->size() != 1) && (calculateCenter( _new, *iter).IsZero())) {
     241        errors.first = 2.;
     242        errors.second = 2.;
     243      }
     244    }
     245  }
     246  LOG(4, "INFO: Resulting errors for matching (L1, L2): "
    200247      << errors.first << "," << errors.second << ".");
    201248
    … …  
    233280 * \param _matching current matching being build up
    234281 * \param _indices contains still available indices
     282 * \param _remainingweights current weights to fill (each weight a place)
     283 * \param _remainiter iterator over the weights, indicating the current position we match
    235284 * \param _matchingsize
    236285 */
    237286void SphericalPointDistribution::recurseMatchings(
    238287    MatchingControlStructure &_MCS,
    239     IndexList_t &_matching,
     288    IndexTupleList_t &_matching,
    240289    IndexList_t _indices,
    241     unsigned int _matchingsize)
    242 {
    243   LOG(4, "DEBUG: Recursing with current matching " << _matching
     290    WeightsArray_t &_remainingweights,
     291    WeightsArray_t::iterator _remainiter,
     292    const unsigned int _matchingsize
     293    )
     294{
     295  LOG(5, "DEBUG: Recursing with current matching " << _matching
    244296      << ", remaining indices " << _indices
    245       << ", and sought size " << _matching.size()+_matchingsize);
    246   //!> threshold for L1 error below which matching is immediately acceptable
    247   const double L1THRESHOLD = 1e-2;
     297      << ", and remaining weights " << _matchingsize);
    248298  if (!_MCS.foundflag) {
    249     LOG(4, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
     299    LOG(5, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
    250300    if (_matchingsize > 0) {
    251301      // go through all indices
    252302      for (IndexList_t::iterator iter = _indices.begin();
    253303          (iter != _indices.end()) && (!_MCS.foundflag);) {
     304        // check whether we can stay in the current bin or have to move on to next one
     305        if (*_remainiter == 0) {
     306          // we need to move on
     307          if (_remainiter != _remainingweights.end()) {
     308            ++_remainiter;
     309          } else {
     310            // as we check _matchingsize > 0 this should be impossible
     311            ASSERT( 0, "recurseMatchings() - we must not come to this position.");
     312          }
     313        }
     314        // advance in matching to same position
     315        const size_t OldIndex = std::distance(_remainingweights.begin(), _remainiter);
     316        while (_matching.size() <= OldIndex) { // add empty lists of new bin is opened
     317          LOG(6, "DEBUG: Extending _matching.");
     318          _matching.push_back( IndexList_t() );
     319        }
     320        IndexTupleList_t::iterator filliniter = _matching.begin();
     321        std::advance(filliniter, OldIndex);
    254322        // add index to matching
    255         _matching.push_back(*iter);
    256         LOG(5, "DEBUG: Adding " << *iter << " to matching.");
     323        filliniter->push_back(*iter);
     324        --(*_remainiter);
     325        LOG(6, "DEBUG: Adding " << *iter << " to matching at " << OldIndex << ".");
    257326        // remove index but keep iterator to position (is the next to erase element)
    258327        IndexList_t::iterator backupiter = _indices.erase(iter);
    259328        // recurse with decreased _matchingsize
    260         recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
     329        recurseMatchings(_MCS, _matching, _indices, _remainingweights, _remainiter, _matchingsize-1);
    261330        // re-add chosen index and reset index to new position
    262         _indices.insert(backupiter, _matching.back());
     331        _indices.insert(backupiter, filliniter->back());
    263332        iter = backupiter;
    264333        // remove index from _matching to make space for the next one
    265         _matching.pop_back();
     334        filliniter->pop_back();
     335        ++(*_remainiter);
    266336      }
    267337      // gone through all indices then exit recursion
    … …  
    269339        _MCS.foundflag = true;
    270340    } else {
    271       LOG(3, "INFO: Found matching " << _matching);
     341      LOG(4, "INFO: Found matching " << _matching);
    272342      // calculate errors
    273343      std::pair<double, double> errors = calculateErrorOfMatching(
    … …  
    310380  MCS.foundflag = false;
    311381  MCS.bestL2 = std::numeric_limits<double>::max();
     382  // transform lists into arrays
    312383  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
    313       iter != _polygon.end(); ++iter)
     384      iter != _polygon.end(); ++iter) {
    314385    MCS.oldpoints.push_back(iter->first);
     386    MCS.weights.push_back(iter->second);
     387  }
    315388  MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
    316389
    … …  
    320393    IndexList_t indices(_newpolygon.size());
    321394    std::generate(indices.begin(), indices.end(), UniqueNumber);
    322     IndexList_t matching;
     395    IndexTupleList_t matching;
    323396
    324397    // walk through all matchings
    325     const unsigned int matchingsize = _polygon.size();
    326     ASSERT( matchingsize <= indices.size(),
    327         "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
    328     recurseMatchings(MCS, matching, indices, matchingsize);
     398    WeightsArray_t remainingweights = MCS.weights;
     399    unsigned int placesleft = std::accumulate(remainingweights.begin(), remainingweights.end(), 0);
     400    recurseMatchings(MCS, matching, indices, remainingweights, remainingweights.begin(), placesleft);
     401  }
     402  if (MCS.foundflag)
     403    LOG(3, "Found a best matching beneath L1 threshold of " << L1THRESHOLD);
     404  else {
     405    if (MCS.bestL2 < warn_amplitude)
     406      LOG(3, "Picking matching is " << MCS.bestmatching << " with best L2 error of "
     407          << MCS.bestL2);
     408    else if (MCS.bestL2 < L2THRESHOLD)
     409      ELOG(2, "Picking matching is " << MCS.bestmatching
     410          << " with rather large L2 error of " << MCS.bestL2);
     411    else
     412      ASSERT(0, "findBestMatching() - matching "+toString(MCS.bestmatching)
     413          +" has L2 error of "+toString(MCS.bestL2)+" that is too large.");
    329414  }
    330415
    331416  // combine multiple points and create simple IndexList from IndexTupleList
    332   IndexTupleList_t bestmatching;
    333   for (IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    334       iter != MCS.bestmatching.end(); ++iter)
    335     bestmatching.push_back(IndexList_t(1, *iter));
    336417  const SphericalPointDistribution::IndexList_t IndexList =
    337       joinPoints(_newpolygon, MCS.newpoints, bestmatching);
     418      joinPoints(_newpolygon, MCS.newpoints, MCS.bestmatching);
    338419
    339420  return IndexList;
    … …  
    364445      Center.Normalize();
    365446      _newpolygon.push_back(Center);
    366       LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
     447      LOG(5, "DEBUG: Combining " << tupleiter->size() << " points to weighted center "
    367448          << Center << " with new index " << UniqueIndex);
    368449      // mark for removal
    … …  
    605686      iter != _polygon.end(); ++iter)
    606687    remainingold.push_back(iter->first);
    607   VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    608688  LOG(2, "INFO: Matching old polygon " << _polygon
    609689      << " with new polygon " << _newpolygon);
    … …  
    618698    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
    619699    LOG(2, "INFO: Best matching is " << bestmatching);
     700
     701    // _newpolygon has changed, so now convert to array
     702    VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    620703
    621704    // determine rotation angles to align the two point distributions with

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r450adf r653cea  
    8181  //!> precalculated value for root of 3
    8282  static const double SQRT_3;
     83  //!> threshold for L1 error below which matching is immediately acceptable
     84  static const double L1THRESHOLD;
     85  //!> threshold for L2 error below which matching is acceptable
     86  static const double L2THRESHOLD;
    8387
    8488  //!> typedef for a full rotation specification consisting of axis and angle.
    … …  
    9599  //!> typedef for a Vector of positions with weights
    96100  typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
     101  //!> typedef for a vector of degrees (or integral weights)
     102  typedef std::vector<unsigned int> WeightsArray_t;
    97103
    98104  //!> amplitude up to which deviations in checks of rotations are tolerated
    … …  
    104110
    105111  static std::pair<double, double> calculateErrorOfMatching(
    106       const std::vector<Vector> &_old,
    107       const std::vector<Vector> &_new,
    108       const IndexList_t &_Matching);
     112      const VectorArray_t &_old,
     113      const VectorArray_t &_new,
     114      const IndexTupleList_t &_Matching);
    109115
    110116  static Polygon_t removeMatchingPoints(
    … …  
    116122    bool foundflag;
    117123    double bestL2;
    118     IndexList_t bestmatching;
     124    IndexTupleList_t bestmatching;
    119125    VectorArray_t oldpoints;
    120126    VectorArray_t newpoints;
     127    WeightsArray_t weights;
    121128  };
    122129
    123130  static void recurseMatchings(
    124131      MatchingControlStructure &_MCS,
    125       IndexList_t &_matching,
     132      IndexTupleList_t &_matching,
    126133      IndexList_t _indices,
    127       unsigned int _matchingsize);
     134      WeightsArray_t &_remainingweights,
     135      WeightsArray_t::iterator _remainiter,
     136      const unsigned int _matchingsize
     137      );
    128138
    129139  static IndexList_t findBestMatching(

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.cpp

    r450adf r653cea  
    6868CPPUNIT_TEST_SUITE_REGISTRATION( SphericalPointDistributionTest );
    6969
     70/** due to root-taking in function we only have limited numerical precision,
     71 * basically half of the double range.
     72 */
     73const double CenterAccuracy = sqrt(std::numeric_limits<double>::epsilon()*1e2);
    7074
    7175void SphericalPointDistributionTest::setUp()
    … …  
    638642}
    639643
     644/** UnitTest for matchSphericalPointDistributions() with four points and weights
     645 * not all equal to one.
     646 */
     647void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_multiple()
     648{
     649  SphericalPointDistribution SPD(1.);
     650
     651  // test with four points: one point having weight of two
     652  {
     653    SphericalPointDistribution::WeightedPolygon_t polygon;
     654    polygon += std::make_pair( Vector(1.,0.,0.), 2);
     655    SphericalPointDistribution::Polygon_t newpolygon =
     656        SPD.get<4>();
     657    SphericalPointDistribution::Polygon_t expected;
     658    expected += Vector(-0.5773502691896,-5.551115123126e-17,0.8164965809277);
     659    expected += Vector(-0.5773502691896,-5.551115123126e-17,-0.8164965809277);
     660    SphericalPointDistribution::Polygon_t remaining =
     661        SphericalPointDistribution::matchSphericalPointDistributions(
     662            polygon,
     663            newpolygon);
     664//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     665//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     666    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     667  }
     668
     669  // test with five points: one point having weight of two
     670  {
     671    SphericalPointDistribution::WeightedPolygon_t polygon;
     672    polygon += std::make_pair( Vector(1.,0.,0.), 2);
     673    SphericalPointDistribution::Polygon_t newpolygon =
     674        SPD.get<5>();
     675    SphericalPointDistribution::Polygon_t expected;
     676    expected += Vector(-0.7071067811865,0.7071067811865,0);
     677    expected += Vector(-0.3535533905933,-0.3535533905933,0.8660254037844);
     678    expected += Vector(-0.3535533905933,-0.3535533905933,-0.8660254037844);
     679    SphericalPointDistribution::Polygon_t remaining =
     680        SphericalPointDistribution::matchSphericalPointDistributions(
     681            polygon,
     682            newpolygon);
     683//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     684//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     685    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     686  }
     687
     688
     689  // test with five points: one point having weight of two, one weight of one
     690  {
     691    SphericalPointDistribution::WeightedPolygon_t polygon;
     692    polygon += std::make_pair( Vector(M_SQRT1_2,M_SQRT1_2,0.), 2);
     693    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
     694    SphericalPointDistribution::Polygon_t newpolygon =
     695        SPD.get<5>();
     696    SphericalPointDistribution::Polygon_t expected;
     697    expected += Vector(0.3535533786708,-0.3535533955317,-0.8660254066357);
     698    expected += Vector(0.3535534025157,-0.3535533856548,0.8660254009332);
     699    SphericalPointDistribution::Polygon_t remaining =
     700        SphericalPointDistribution::matchSphericalPointDistributions(
     701            polygon,
     702            newpolygon);
     703//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     704//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     705    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     706  }
     707
     708  // test with six points: two points each having weight of two
     709  {
     710    SphericalPointDistribution::WeightedPolygon_t polygon;
     711    polygon += std::make_pair( Vector(M_SQRT1_2,-M_SQRT1_2,0.), 2);
     712    polygon += std::make_pair( Vector(-M_SQRT1_2,M_SQRT1_2,0.), 2);
     713    SphericalPointDistribution::Polygon_t newpolygon =
     714        SPD.get<6>();
     715    SphericalPointDistribution::Polygon_t expected;
     716    expected += Vector(0.,0.,-1.);
     717    expected += Vector(0.,0.,1.);
     718    SphericalPointDistribution::Polygon_t remaining =
     719        SphericalPointDistribution::matchSphericalPointDistributions(
     720            polygon,
     721            newpolygon);
     722//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     723//    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     724    CPPUNIT_ASSERT( areEqualToWithinBounds(expected, remaining, CenterAccuracy) );
     725  }
     726}
     727
    640728/** UnitTest for matchSphericalPointDistributions() with five points
    641729 */

  • src/Fragmentation/Exporters/unittests/SphericalPointDistributionUnitTest.hpp

    r450adf r653cea  
    3131    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_7 );
    3232    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_8 );
     33    CPPUNIT_TEST ( matchSphericalPointDistributionsTest_multiple );
    3334    CPPUNIT_TEST_SUITE_END();
    3435
    … …  
    4546      void matchSphericalPointDistributionsTest_7();
    4647      void matchSphericalPointDistributionsTest_8();
     48      void matchSphericalPointDistributionsTest_multiple();
    4749
    4850private:
Note: See TracChangeset for help on using the changeset viewer.