Changeset 23c605 for src/Fragmentation


Ignore:
Timestamp:
Aug 20, 2014, 1:06:16 PM (11 years ago)
Author:
Frederik Heber <heber@…>
Children:
1cde4e8
Parents:
1ae9aa
git-author:
Frederik Heber <heber@…> (06/30/14 09:35:42)
git-committer:
Frederik Heber <heber@…> (08/20/14 13:06:16)
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.
  • TESTS: Regresssion test FragmentMolecule-cycles working again.
Location:
src/Fragmentation/Exporters
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r1ae9aa r23c605  
    4949#include <limits>
    5050#include <list>
     51#include <numeric>
    5152#include <vector>
    5253#include <map>
    … …  
    5960// static entities
    6061const double SphericalPointDistribution::warn_amplitude = 1e-2;
     62const double SphericalPointDistribution::L1THRESHOLD = 1e-2;
     63const double SphericalPointDistribution::L2THRESHOLD = 2e-1;
    6164
    6265typedef std::vector<double> DistanceArray_t;
    … …  
    6467// class generator: taken from www.cplusplus.com example std::generate
    6568struct c_unique {
    66   int current;
     69  unsigned int current;
    6770  c_unique() {current=0;}
    68   int operator()() {return current++;}
     71  unsigned int operator()() {return current++;}
    6972} UniqueNumber;
    7073
    71 inline
    72 DistanceArray_t calculatePairwiseDistances(
    73     const std::vector<Vector> &_points,
    74     const SphericalPointDistribution::IndexList_t &_indices
    75     )
    76 {
    77   DistanceArray_t result;
    78   for (SphericalPointDistribution::IndexList_t::const_iterator firstiter = _indices.begin();
    79       firstiter != _indices.end(); ++firstiter) {
    80     for (SphericalPointDistribution::IndexList_t::const_iterator seconditer = firstiter;
    81         seconditer != _indices.end(); ++seconditer) {
    82       if (firstiter == seconditer)
    83         continue;
    84       const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
    85       result.push_back(distance);
    86     }
    87   }
    88   return result;
    89 }
     74struct c_unique_list {
     75  unsigned int current;
     76  c_unique_list() {current=0;}
     77  std::list<unsigned int> operator()() {return std::list<unsigned int>(1, current++);}
     78} UniqueNumberList;
    9079
    9180/** Calculate the center of a given set of points in \a _positions but only
    … …  
    10796
    10897  return Center;
     98}
     99
     100inline
     101DistanceArray_t calculatePairwiseDistances(
     102    const SphericalPointDistribution::VectorArray_t &_points,
     103    const SphericalPointDistribution::IndexTupleList_t &_indices
     104    )
     105{
     106  DistanceArray_t result;
     107  for (SphericalPointDistribution::IndexTupleList_t::const_iterator firstiter = _indices.begin();
     108      firstiter != _indices.end(); ++firstiter) {
     109
     110    // calculate first center from possible tuple of indices
     111    Vector FirstCenter;
     112    ASSERT(!firstiter->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
     113    if (firstiter->size() == 1) {
     114      FirstCenter = _points[*firstiter->begin()];
     115    } else {
     116      FirstCenter = calculateCenter( _points, *firstiter);
     117      if (!FirstCenter.IsZero())
     118        FirstCenter.Normalize();
     119    }
     120
     121    for (SphericalPointDistribution::IndexTupleList_t::const_iterator seconditer = firstiter;
     122        seconditer != _indices.end(); ++seconditer) {
     123      if (firstiter == seconditer)
     124        continue;
     125
     126      // calculate second center from possible tuple of indices
     127      Vector SecondCenter;
     128      ASSERT(!seconditer->empty(), "calculatePairwiseDistances() - there is an empty tuple.");
     129      if (seconditer->size() == 1) {
     130        SecondCenter = _points[*seconditer->begin()];
     131      } else {
     132        SecondCenter = calculateCenter( _points, *seconditer);
     133        if (!SecondCenter.IsZero())
     134          SecondCenter.Normalize();
     135      }
     136
     137      // calculate distance between both centers
     138      double distance = 2.; // greatest distance on surface of sphere with radius 1.
     139      if ((!FirstCenter.IsZero()) && (!SecondCenter.IsZero()))
     140        distance = (FirstCenter - SecondCenter).NormSquared();
     141      result.push_back(distance);
     142    }
     143  }
     144  return result;
    109145}
    110146
    … …  
    123159  Vector dreiBein(_oldPosition);
    124160  dreiBein.VectorProduct(_RotationAxis);
     161  ASSERT( !dreiBein.IsZero(), "determineSignOfRotation() - dreiBein is zero.");
    125162  dreiBein.Normalize();
    126163  const double sign =
    127164      (dreiBein.ScalarProduct(_newPosition) < 0.) ? -1. : +1.;
    128165  LOG(6, "DEBUG: oldCenter on plane is " << _oldPosition
    129       << ", newCenter in plane is " << _newPosition
     166      << ", newCenter on plane is " << _newPosition
    130167      << ", and dreiBein is " << dreiBein);
    131168  return sign;
    … …  
    166203 */
    167204std::pair<double, double> SphericalPointDistribution::calculateErrorOfMatching(
    168     const std::vector<Vector> &_old,
    169     const std::vector<Vector> &_new,
    170     const IndexList_t &_Matching)
     205    const VectorArray_t &_old,
     206    const VectorArray_t &_new,
     207    const IndexTupleList_t &_Matching)
    171208{
    172209  std::pair<double, double> errors( std::make_pair( 0., 0. ) );
    173210
    174211  if (_Matching.size() > 1) {
    175     LOG(3, "INFO: Matching is " << _Matching);
     212    LOG(5, "INFO: Matching is " << _Matching);
    176213
    177214    // calculate all pair-wise distances
    178     IndexList_t keys(_Matching.size());
    179     std::generate (keys.begin(), keys.end(), UniqueNumber);
     215    IndexTupleList_t keys(_old.size(), IndexList_t() );
     216    std::generate (keys.begin(), keys.end(), UniqueNumberList);
     217
    180218    const DistanceArray_t firstdistances = calculatePairwiseDistances(_old, keys);
    181219    const DistanceArray_t seconddistances = calculatePairwiseDistances(_new, _Matching);
    … …  
    187225    for (;(firstiter != firstdistances.end()) && (seconditer != seconddistances.end());
    188226        ++firstiter, ++seconditer) {
    189       const double gap = *firstiter - *seconditer;
     227      const double gap = fabs(*firstiter - *seconditer);
    190228      // L1 error
    191229      if (errors.first < gap)
    … …  
    194232      errors.second += gap*gap;
    195233    }
    196   } else
    197     ELOG(3, "calculateErrorOfMatching() - Given matching's size is less than 2.");
    198   LOG(3, "INFO: Resulting errors for matching (L1, L2): "
     234  } else {
     235    // check whether we have any zero centers: Combining points on new sphere may result
     236    // in zero centers
     237    for (SphericalPointDistribution::IndexTupleList_t::const_iterator iter = _Matching.begin();
     238        iter != _Matching.end(); ++iter) {
     239      if ((iter->size() != 1) && (calculateCenter( _new, *iter).IsZero())) {
     240        errors.first = 2.;
     241        errors.second = 2.;
     242      }
     243    }
     244  }
     245  LOG(4, "INFO: Resulting errors for matching (L1, L2): "
    199246      << errors.first << "," << errors.second << ".");
    200247
    … …  
    232279 * \param _matching current matching being build up
    233280 * \param _indices contains still available indices
     281 * \param _remainingweights current weights to fill (each weight a place)
     282 * \param _remainiter iterator over the weights, indicating the current position we match
    234283 * \param _matchingsize
    235284 */
    236285void SphericalPointDistribution::recurseMatchings(
    237286    MatchingControlStructure &_MCS,
    238     IndexList_t &_matching,
     287    IndexTupleList_t &_matching,
    239288    IndexList_t _indices,
    240     unsigned int _matchingsize)
    241 {
    242   LOG(4, "DEBUG: Recursing with current matching " << _matching
     289    WeightsArray_t &_remainingweights,
     290    WeightsArray_t::iterator _remainiter,
     291    const unsigned int _matchingsize
     292    )
     293{
     294  LOG(5, "DEBUG: Recursing with current matching " << _matching
    243295      << ", remaining indices " << _indices
    244       << ", and sought size " << _matching.size()+_matchingsize);
    245   //!> threshold for L1 error below which matching is immediately acceptable
    246   const double L1THRESHOLD = 1e-2;
     296      << ", and remaining weights " << _matchingsize);
    247297  if (!_MCS.foundflag) {
    248     LOG(4, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
     298    LOG(5, "DEBUG: Current matching has size " << _matching.size() << ", places left " << _matchingsize);
    249299    if (_matchingsize > 0) {
    250300      // go through all indices
    251301      for (IndexList_t::iterator iter = _indices.begin();
    252302          (iter != _indices.end()) && (!_MCS.foundflag);) {
     303        // check whether we can stay in the current bin or have to move on to next one
     304        if (*_remainiter == 0) {
     305          // we need to move on
     306          if (_remainiter != _remainingweights.end()) {
     307            ++_remainiter;
     308          } else {
     309            // as we check _matchingsize > 0 this should be impossible
     310            ASSERT( 0, "recurseMatchings() - we must not come to this position.");
     311          }
     312        }
     313        // advance in matching to same position
     314        const size_t OldIndex = std::distance(_remainingweights.begin(), _remainiter);
     315        while (_matching.size() <= OldIndex) { // add empty lists of new bin is opened
     316          LOG(6, "DEBUG: Extending _matching.");
     317          _matching.push_back( IndexList_t() );
     318        }
     319        IndexTupleList_t::iterator filliniter = _matching.begin();
     320        std::advance(filliniter, OldIndex);
    253321        // add index to matching
    254         _matching.push_back(*iter);
    255         LOG(5, "DEBUG: Adding " << *iter << " to matching.");
     322        filliniter->push_back(*iter);
     323        --(*_remainiter);
     324        LOG(6, "DEBUG: Adding " << *iter << " to matching at " << OldIndex << ".");
    256325        // remove index but keep iterator to position (is the next to erase element)
    257326        IndexList_t::iterator backupiter = _indices.erase(iter);
    258327        // recurse with decreased _matchingsize
    259         recurseMatchings(_MCS, _matching, _indices, _matchingsize-1);
     328        recurseMatchings(_MCS, _matching, _indices, _remainingweights, _remainiter, _matchingsize-1);
    260329        // re-add chosen index and reset index to new position
    261         _indices.insert(backupiter, _matching.back());
     330        _indices.insert(backupiter, filliniter->back());
    262331        iter = backupiter;
    263332        // remove index from _matching to make space for the next one
    264         _matching.pop_back();
     333        filliniter->pop_back();
     334        ++(*_remainiter);
    265335      }
    266336      // gone through all indices then exit recursion
    … …  
    268338        _MCS.foundflag = true;
    269339    } else {
    270       LOG(3, "INFO: Found matching " << _matching);
     340      LOG(4, "INFO: Found matching " << _matching);
    271341      // calculate errors
    272342      std::pair<double, double> errors = calculateErrorOfMatching(
    … …  
    309379  MCS.foundflag = false;
    310380  MCS.bestL2 = std::numeric_limits<double>::max();
     381  // transform lists into arrays
    311382  for (WeightedPolygon_t::const_iterator iter = _polygon.begin();
    312       iter != _polygon.end(); ++iter)
     383      iter != _polygon.end(); ++iter) {
    313384    MCS.oldpoints.push_back(iter->first);
     385    MCS.weights.push_back(iter->second);
     386  }
    314387  MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
    315388
    … …  
    319392    IndexList_t indices(_newpolygon.size());
    320393    std::generate(indices.begin(), indices.end(), UniqueNumber);
    321     IndexList_t matching;
     394    IndexTupleList_t matching;
    322395
    323396    // walk through all matchings
    324     const unsigned int matchingsize = _polygon.size();
    325     ASSERT( matchingsize <= indices.size(),
    326         "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
    327     recurseMatchings(MCS, matching, indices, matchingsize);
     397    WeightsArray_t remainingweights = MCS.weights;
     398    unsigned int placesleft = std::accumulate(remainingweights.begin(), remainingweights.end(), 0);
     399    recurseMatchings(MCS, matching, indices, remainingweights, remainingweights.begin(), placesleft);
     400  }
     401  if (MCS.foundflag)
     402    LOG(3, "Found a best matching beneath L1 threshold of " << L1THRESHOLD);
     403  else {
     404    if (MCS.bestL2 < warn_amplitude)
     405      LOG(3, "Picking matching is " << MCS.bestmatching << " with best L2 error of "
     406          << MCS.bestL2);
     407    else if (MCS.bestL2 < L2THRESHOLD)
     408      ELOG(2, "Picking matching is " << MCS.bestmatching
     409          << " with rather large L2 error of " << MCS.bestL2);
     410    else
     411      ASSERT(0, "findBestMatching() - matching "+toString(MCS.bestmatching)
     412          +" has L2 error of "+toString(MCS.bestL2)+" that is too large.");
    328413  }
    329414
    330415  // combine multiple points and create simple IndexList from IndexTupleList
    331   IndexTupleList_t bestmatching;
    332   for (IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    333       iter != MCS.bestmatching.end(); ++iter)
    334     bestmatching.push_back(IndexList_t(1, *iter));
    335416  const SphericalPointDistribution::IndexList_t IndexList =
    336       joinPoints(_newpolygon, MCS.newpoints, bestmatching);
     417      joinPoints(_newpolygon, MCS.newpoints, MCS.bestmatching);
    337418
    338419  return IndexList;
    … …  
    363444      Center.Normalize();
    364445      _newpolygon.push_back(Center);
    365       LOG(5, "DEBUG: Combining " << tupleiter->size() << "points to weighted center "
     446      LOG(5, "DEBUG: Combining " << tupleiter->size() << " points to weighted center "
    366447          << Center << " with new index " << UniqueIndex);
    367448      // mark for removal
    … …  
    600681      iter != _polygon.end(); ++iter)
    601682    remainingold.push_back(iter->first);
    602   VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    603683  LOG(2, "INFO: Matching old polygon " << _polygon
    604684      << " with new polygon " << _newpolygon);
    … …  
    613693    IndexList_t bestmatching = findBestMatching(_polygon, _newpolygon);
    614694    LOG(2, "INFO: Best matching is " << bestmatching);
     695
     696    // _newpolygon has changed, so now convert to array
     697    VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    615698
    616699    // determine rotation angles to align the two point distributions with

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r1ae9aa r23c605  
    8989  //!> precalculated value for root of 3
    9090  const double SQRT_3;
     91  //!> threshold for L1 error below which matching is immediately acceptable
     92  static const double L1THRESHOLD;
     93  //!> threshold for L2 error below which matching is acceptable
     94  static const double L2THRESHOLD;
    9195
    9296  //!> typedef for a full rotation specification consisting of axis and angle.
    … …  
    103107  //!> typedef for a Vector of positions with weights
    104108  typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
     109  //!> typedef for a vector of degrees (or integral weights)
     110  typedef std::vector<unsigned int> WeightsArray_t;
    105111
    106112  //!> amplitude up to which deviations in checks of rotations are tolerated
    … …  
    112118
    113119  static std::pair<double, double> calculateErrorOfMatching(
    114       const std::vector<Vector> &_old,
    115       const std::vector<Vector> &_new,
    116       const IndexList_t &_Matching);
     120      const VectorArray_t &_old,
     121      const VectorArray_t &_new,
     122      const IndexTupleList_t &_Matching);
    117123
    118124  static Polygon_t removeMatchingPoints(
    … …  
    124130    bool foundflag;
    125131    double bestL2;
    126     IndexList_t bestmatching;
     132    IndexTupleList_t bestmatching;
    127133    VectorArray_t oldpoints;
    128134    VectorArray_t newpoints;
     135    WeightsArray_t weights;
    129136  };
    130137
    131138  static void recurseMatchings(
    132139      MatchingControlStructure &_MCS,
    133       IndexList_t &_matching,
     140      IndexTupleList_t &_matching,
    134141      IndexList_t _indices,
    135       unsigned int _matchingsize);
     142      WeightsArray_t &_remainingweights,
     143      WeightsArray_t::iterator _remainiter,
     144      const unsigned int _matchingsize
     145      );
    136146
    137147  static IndexList_t findBestMatching(

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

    r1ae9aa r23c605  
    638638}
    639639
     640/** UnitTest for matchSphericalPointDistributions() with four points and weights
     641 * not all equal to one.
     642 */
     643void SphericalPointDistributionTest::matchSphericalPointDistributionsTest_multiple()
     644{
     645  SphericalPointDistribution SPD(1.);
     646
     647  // test with four points: one point having weight of two
     648  {
     649    SphericalPointDistribution::WeightedPolygon_t polygon;
     650    polygon += std::make_pair( Vector(1.,0.,0.), 2);
     651    SphericalPointDistribution::Polygon_t newpolygon =
     652        SPD.get<4>();
     653    SphericalPointDistribution::Polygon_t expected;
     654    expected += Vector(-0.5773502691896,-5.551115123126e-17,0.8164965809277);
     655    expected += Vector(-0.5773502691896,-5.551115123126e-17,-0.8164965809277);
     656    SphericalPointDistribution::Polygon_t remaining =
     657        SphericalPointDistribution::matchSphericalPointDistributions(
     658            polygon,
     659            newpolygon);
     660//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     661    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     662  }
     663
     664  // test with five points: one point having weight of two
     665  {
     666    SphericalPointDistribution::WeightedPolygon_t polygon;
     667    polygon += std::make_pair( Vector(1.,0.,0.), 2);
     668    SphericalPointDistribution::Polygon_t newpolygon =
     669        SPD.get<5>();
     670    SphericalPointDistribution::Polygon_t expected;
     671    expected += Vector(-0.7071067811865,0.7071067811865,0);
     672    expected += Vector(-0.3535533905933,-0.3535533905933,0.8660254037844);
     673    expected += Vector(-0.3535533905933,-0.3535533905933,-0.8660254037844);
     674    SphericalPointDistribution::Polygon_t remaining =
     675        SphericalPointDistribution::matchSphericalPointDistributions(
     676            polygon,
     677            newpolygon);
     678//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     679    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     680  }
     681
     682
     683  // test with five points: one point having weight of two, one weight of one
     684  {
     685    SphericalPointDistribution::WeightedPolygon_t polygon;
     686    polygon += std::make_pair( Vector(M_SQRT1_2,M_SQRT1_2,0.), 2);
     687    polygon += std::make_pair( Vector(-1.,0.,0.), 1);
     688    SphericalPointDistribution::Polygon_t newpolygon =
     689        SPD.get<5>();
     690    SphericalPointDistribution::Polygon_t expected;
     691    expected += Vector(0.3535533786708,-0.3535533955317,-0.8660254066357);
     692    expected += Vector(0.3535534025157,-0.3535533856548,0.8660254009332);
     693    SphericalPointDistribution::Polygon_t remaining =
     694        SphericalPointDistribution::matchSphericalPointDistributions(
     695            polygon,
     696            newpolygon);
     697//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     698    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     699  }
     700
     701  // test with six points: two points each having weight of two
     702  {
     703    SphericalPointDistribution::WeightedPolygon_t polygon;
     704    polygon += std::make_pair( Vector(M_SQRT1_2,-M_SQRT1_2,0.), 2);
     705    polygon += std::make_pair( Vector(-M_SQRT1_2,M_SQRT1_2,0.), 2);
     706    SphericalPointDistribution::Polygon_t newpolygon =
     707        SPD.get<6>();
     708    SphericalPointDistribution::Polygon_t expected;
     709    expected += Vector(0.,0.,-1.);
     710    expected += Vector(0.,0.,1.);
     711    SphericalPointDistribution::Polygon_t remaining =
     712        SphericalPointDistribution::matchSphericalPointDistributions(
     713            polygon,
     714            newpolygon);
     715//    std::cout << std::setprecision(13) << "Matched polygon is " << remaining << std::endl;
     716    CPPUNIT_ASSERT_EQUAL( expected, remaining );
     717  }
     718}
     719
    640720/** UnitTest for matchSphericalPointDistributions() with five points
    641721 */

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

    r1ae9aa r23c605  
    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.