LeastContributorApproximator.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Approximately determines the individual contributing the least
5  * hypervolume.
6  *
7  *
8  *
9  * \author T.Voss
10  * \date 2010-2011
11  *
12  *
13  * \par Copyright 1995-2015 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://image.diku.dk/shark/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33 #ifndef SHARK_ALGORITHMS_DIRECT_SEARCH_INDICATORS_LEASTCONTRIBUTORAPPROXIMATOR_H
34 #define SHARK_ALGORITHMS_DIRECT_SEARCH_INDICATORS_LEASTCONTRIBUTORAPPROXIMATOR_H
35 
36 #include <shark/Core/Exception.h>
37 #include <shark/Core/Math.h>
38 
39 #include <boost/assign.hpp>
40 #include <boost/bimap.hpp>
41 #include <boost/cstdint.hpp>
42 
43 #include <algorithm>
44 #include <iostream>
45 #include <iterator>
46 #include <limits>
47 #include <vector>
48 #include <cmath>
49 
50 namespace shark {
51 
52 /**
53  * \brief Samples a random point.
54  */
55 template<typename Rng>
56 struct Sampler {
57 
58  /**
59  * \brief Samples a random point and checks whether it is dominated or not.
60  *
61  * \tparam Set The type of the set of points.
62  *
63  * \param [in] s The set of individuals to check the sampled point against.
64  * \param [in] point Iterator to the point in the supplied that shall serve as basis for the random point.
65  *
66  * \returns true if the sample was successful (i.e. non-dominated) and false otherwise.
67  *
68  */
69  template<typename Set>
70  bool operator()( const Set & s, typename Set::iterator point )
71  {
72  point->m_sample = point->m_point;
73  for( unsigned int i = 0; i < point->m_sample.size(); i++ ) {
74  point->m_sample[ i ] += Rng::uni( 0., 1. ) * ( point->m_boundingBox[ i ] - point->m_point[ i ] );
75  }
76 
77  bool successful = true;
78  for( unsigned int i = 0; i < point->m_influencingPoints.size(); i++ ) {
79  point->m_noOperations += point->m_sample.size() + 1;
80  bool inside = true;
81  for( unsigned int j = 0; j < point->m_sample.size(); j++ ) {
82  if( point->m_influencingPoints[i]->m_point[ j ] > point->m_sample[ j ] ) {
83  inside = false;
84  break;
85  }
86  }
87 
88  if( inside ) {
89  successful = false;
90  break;
91  }
92  }
93 
94  return successful;
95  }
96 };
97 
98 /**
99  * \brief Calculates bounding boxes.
100  * \tparam Set The type of the set of points.
101  */
102 template<typename Set>
104 
105  /**
106  * \brief Compares points based on their contributed volume.
107  */
109  /**
110  * \brief Compares points based on their contributed volume.
111  * \returns true if the volume contributed by i1 is larger than the volume contributed by i2, false otherwise.
112  */
113  template<typename Iterator>
114  bool operator()( Iterator i1, Iterator i2 )
115  {
116  return i1->m_overlappingVolume > i2->m_overlappingVolume;
117  }
118  };
119 
120  /**
121  * \brief Computes bounding boxes and their volume for the range of points defined by the iterators.
122  * \param [in] begin Iterator pointing to the first valid point.
123  * \param [in] end Iterator pointing right after the last valid point.
124  */
125  void operator()( typename Set::iterator begin, typename Set::iterator end )
126  {
127 
128  for( typename Set::iterator it = begin; it != end; ++it ) {
129  typename Set::value_type & p1 = *it;
130  for( typename Set::iterator itt = begin; itt != end; ++itt ) {
131 
132  if( itt == it )
133  continue;
134 
135  typename Set::value_type & p2 = *itt;
136 
137  unsigned int coordCounter = 0;
138  unsigned int coord = 0;
139  for( unsigned int o = 0; o < p1.m_point.size(); o++ ) {
140  if( p2.m_point[ o ] > p1.m_point[ o ] ) {
141  coordCounter++;
142  if( coordCounter == 2 )
143  break;
144  coord = o;
145  }
146  }
147 
148  if( coordCounter == 1 && p1.m_boundingBox[ coord ] > p2.m_point[ coord ] )
149  p1.m_boundingBox[ coord ] = p2.m_point[ coord ];
150  }
151 
152  it->m_boundingBoxVolume = 1.;
153  for( unsigned int i = 0 ; i < it->m_boundingBox.size(); i++ ) {
154 
155  it->m_boundingBoxVolume *= it->m_boundingBox[ i ] - it->m_point[ i ];
156 
157  }
158 
159  for( typename Set::iterator itt = begin; itt != end; ++itt ) {
160  if( itt == it )
161  continue;
162 
163  bool isInfluencing = true;
164  double vol = 1.;
165  for( unsigned int i = 0; i < it->m_point.size(); i++ ) {
166 
167  if( itt->m_point[ i ] >= it->m_boundingBox[ i ] ) {
168  isInfluencing = false;
169  break;
170  }
171  vol *= std::max( itt->m_point[ i ], it->m_point[ i ] ) - it->m_boundingBox[ i ];//min_dbl(_alc_front[k][l],_alc_front[i][l]) - _alc_boundBoxLower[_alc_noObjectives*i+l];
172  }
173  if( isInfluencing ) {
174  itt->m_overlappingVolume = vol;
175  it->m_influencingPoints.push_back( itt );
176  }
177  }
178 
179  std::sort( it->m_influencingPoints.begin(), it->m_influencingPoints.end(), VolumeComparator() );
180  }
181  }
182 };
183 
184 /**
185  * \brief Approximately determines the point of a set contributing the least hypervolume.
186  *
187  * See K. Bringmann, T. Friedrich. Approximating the least hypervolume contributor: NP-hard in general, but fast in practice. Proc. of the 5th International Conference on Evolutionary Multi-Criterion Optimization (EMO 2009), Vol. 5467 of LNCS, pages 6-20, Springer-Verlag, 2009.
188  *
189  * \tparam Rng The type of the Rng for sampling random points.
190  * \tparam ExactHypervolume Exact hypervolume calculator that is used to speed up the approximation scheme.
191  */
192 template<typename Rng, typename ExactHypervolume>
194 
195  /**
196  * \brief Typedefs the type including all templates.
197  */
199 
200  /**
201  * \brief Models the sampling strategy.
202  */
203  enum Strategy {
204  ENSURE_UNION_BOUND, ///< Ensures that the selected epsilon and delta are reached.
205  FAST ///< Puts a threshold on the number of samples.
206  };
207 
208  typedef boost::bimap< Strategy, std::string > registry_type;
209  static registry_type & STRATEGY_REGISTRY()
210  {
211  static registry_type registry = boost::assign::list_of< typename registry_type::relation >( ENSURE_UNION_BOUND, "EnsureUnionBound" )( FAST, "Fast" );
212  return( registry );
213  }
214 
215  friend std::ostream & operator<<( std::ostream & stream, Strategy strategy )
216  {
217  stream << STRATEGY_REGISTRY().left.find( strategy )->second;
218  return( stream );
219  }
220 
221  friend std::istream & operator>>( std::istream & stream, Strategy & strategy )
222  {
223  static std::string s;
224  stream >> s;
225  strategy = STRATEGY_REGISTRY().right.find( s )->second;
226  return( stream );
227  }
228 
229  /**
230  * \brief Default sampling strategy.
231  */
233  {
234  return( FAST );
235  }
236 
237  /**
238  * \brief Default sampling strategy.
239  */
240  static boost::uint_fast64_t DEFAULT_SAMPLE_THRESHOLD()
241  {
242  return( 1000000 );
243  }
244 
245  /**
246  * \brief Default delta value at start of algorithm.
247  */
248  static const double DEFAULT_START_DELTA()
249  {
250  return( 0.1 );
251  }
252 
253  /**
254  * \brief Default multiplier value for adjusting delta.
255  */
256  static const double DEFAULT_MULTIPLIER_DELTA()
257  {
258  return( 0.775 );
259  }
260 
261  /**
262  * \brief Default multiplier value for adjusting delta.
263  */
264  static const double DEFAULT_MINIMUM_MULTIPLIER_DELTA()
265  {
266  return( 0.2 );
267  }
268 
269  /**
270  * \brief Default threshold for sample count.
271  */
272  static const boost::uint_fast64_t DEFAULT_MAX_NUM_SAMPLES()
273  {
275  }
276 
277  /**
278  * \brief Default gamma value.
279  */
280  static const double DEFAULT_GAMMA()
281  {
282  return( 0.25 );
283  }
284 
285  /**
286  * \brief Models a point and associated information for book-keeping purposes.
287  */
288  template<typename VectorType>
289  struct Point {
290 
291  typedef typename std::list< VectorType >::iterator sample_iterator;
292  typedef typename std::list< VectorType >::const_iterator const_sample_iterator;
293 
294  Point( std::size_t noObjectives, const VectorType & point, const VectorType & refPoint ) : m_point( point ),
295  m_sample( point.size() ),
296  m_boundingBox( refPoint ),
297  m_boundingBoxVolume( 0. ),
298  m_approximatedContribution( 0. ),
299  m_noOperations( 0 ),
300  m_noSamples( 0 ),
301  m_noSuccessfulSamples( 0 )
302  {
303  }
304 
308 
309  std::list< VectorType > m_successfulSamples;
310  std::list< VectorType > m_nonSuccessfulSamples;
311 
312  std::vector< typename std::vector<Point>::const_iterator > m_influencingPoints;
313 
317 
318  std::size_t m_noOperations;
319  std::size_t m_noSamples;
321 
322  template<typename Stream>
323  void print( Stream & s ) const
324  {
325 
326  std::copy( m_point.begin(), m_point.end(), std::ostream_iterator<double>( s, " " ) );
327  // s << " || ";
328  std::copy( m_boundingBox.begin(), m_boundingBox.end(), std::ostream_iterator<double>( s, " " ) );
329  // s << " || ";
330  s << m_noOperations << " " << m_noSamples << " " << m_noSuccessfulSamples << " " << m_boundingBoxVolume << std::endl;
331 
332  }
333  };
334 
335  /**
336  * \brief Returns the supplied argument.
337  */
339  template<typename Member>
340  const Member & operator()( const Member & member ) const
341  {
342  return( member );
343  }
344  };
345 
346  double m_logFactor;
347  double m_startDelta;
350  unsigned long long m_maxNumSamples;
351  double m_gamma;
352 
353  unsigned int m_round;
354 
355  double m_errorProbability; ///<The error probability.
356  double m_errorBound; ///<The error bound
357 
359  boost::uint_fast64_t m_sampleCounter;
360  boost::uint_fast64_t m_sampleCountThreshold;
361 
362  RealVector m_reference;///< the reference calculates by updateInternals
363 
364  /**
365  * \brief C'tor
366  * \param [in] startDelta Initial delta value.
367  * \param [in] multiplierDelta Multiplier for adjusting the delta value.
368  * \param [in] minimumDeltaMultiplier Multiplier for adjusting the delta value.
369  * \param [in] maxNumSamples The maximum number of samples. If reached, the algorithm aborts.
370  * \param [in] gamma
371  * \param [in] delta the error probability of the least contributor
372  * \param [in] eps the error bound of the least contributor
373  */
374  LeastContributorApproximator( double startDelta = this_type::DEFAULT_START_DELTA(),
375  double multiplierDelta = this_type::DEFAULT_MULTIPLIER_DELTA(),
376  double minimumDeltaMultiplier = this_type::DEFAULT_MINIMUM_MULTIPLIER_DELTA(),
377  unsigned long long maxNumSamples = this_type::DEFAULT_MAX_NUM_SAMPLES(),
378  double gamma = this_type::DEFAULT_GAMMA(),
379  double delta = 1.E-2,
380  double eps = 1.E-2
381  ) : m_startDelta( startDelta ),
382  m_multiplierDelta( multiplierDelta ),
383  m_minimumMultiplierDelta( minimumDeltaMultiplier ),
384  m_maxNumSamples( maxNumSamples ),
385  m_gamma( gamma ),
386  m_round( 0 ),
387  m_errorProbability(delta),
388  m_errorBound(eps),
389  m_strategy( this_type::DEFAULT_STRATEGY() ),
390  m_sampleCounter( 0 ),
391  m_sampleCountThreshold( this_type::DEFAULT_SAMPLE_THRESHOLD() )
392  {
393 
394  }
395 
396  /**
397  * \brief Samples in the bounding box of the supplied point until a pre-defined threshold is reached.
398  * \param [in] s Set of points.
399  * \param [in] point Iterator to the point that should be sampled.
400  * \param [in] r The current round.
401  * \param [in] delta The delta that should be reached.
402  * \param [in] refPoint Reference point for hypervolume calculation/approximation.
403  */
404  template<typename Set, typename VectorType>
405  void sample( const Set & s, typename Set::iterator point, unsigned int r, double delta, const VectorType & refPoint )
406  {
407 
408  if( point->m_noSamples >= m_maxNumSamples )
409  return;
410 
411 
412  if( point->m_noOperations >= ExactHypervolume::runtime( point->m_influencingPoints.size(), refPoint.size() ) ) {
413 
414  point->m_noSamples = point->m_noSuccessfulSamples = m_maxNumSamples;
415 
416  std::vector< VectorType > neighborhood( point->m_influencingPoints.size(), refPoint );
417  for( unsigned int i = 0; i < point->m_influencingPoints.size(); i++ ) {
418  for( unsigned int j = 0; j < refPoint.size(); j++ )
419  neighborhood[i][j] = std::max( point->m_point[ j ], (point->m_influencingPoints[i])->m_point[ j ] ) + ( refPoint[ j ] - point->m_boundingBox[ j ] );
420  }
421 
423  ExactHypervolume hv;
424  point->m_approximatedContribution = point->m_boundingBoxVolume - hv( e, neighborhood, refPoint );
425  point->m_boundingBoxVolume = point->m_approximatedContribution;
426  return;
427  }
428 
429  Sampler<Rng> sampler;
430 
431  double threshold = 0.5 * ( (1. + m_gamma) * ::log( static_cast<double>( r ) ) + m_logFactor ) * sqr( point->m_boundingBoxVolume / delta );
432  for( ; point->m_noSamples == 0 || point->m_noSamples < threshold; point->m_noSamples++ ) {
433  if( m_strategy == FAST )
434  if( m_sampleCounter > m_sampleCountThreshold )
435  throw( SHARKEXCEPTION( "LeastContributorApproximator::sample(): Maximum number of total samples reached." ) );
436  if( sampler( s, point ) )
437  point->m_noSuccessfulSamples++;
438 
439  m_sampleCounter++;
440  }
441 
442  point->m_approximatedContribution = point->m_boundingBoxVolume * ( static_cast<double>( point->m_noSuccessfulSamples ) / static_cast<double>( point->m_noSamples ) );
443  }
444 
445  template<typename iterator>
446  double deltaForPoint( iterator point, unsigned int R )
447  {
448 
449  return( ::sqrt( 0.5 * ((1. + m_gamma) * ::log( static_cast<double>( R ) ) + m_logFactor ) / point->m_noSamples ) * point->m_boundingBoxVolume );
450 
451  }
452 
453  /**
454  * \brief Determines the point contributing the least hypervolume to the overall set of points.
455  * \param [in] e Extracts point information from set elements.
456  * \param [in] s Set of points.
457  * \param [in] refPoint The reference point to consider for calculating individual points' contributions.
458  */
459  template<typename Extractor, typename Set, typename VectorType>
460  std::size_t leastContributor( Extractor e, const Set & s, const VectorType & refPoint)
461  {
462  if(s.empty()) return 0;
463 
464  std::size_t noObjectives = e(s[0]).size();
465 
466  std::vector< Point<VectorType> > front;
467  front.reserve( s.size() );
468 
469  std::vector< typename std::vector< Point<VectorType> >::iterator > activePoints;
470 
471  for( typename Set::const_iterator it = s.begin(); it != s.end(); ++it ) {
472  front.push_back( Point<VectorType>( noObjectives, e( *it ), refPoint ) );
473  }
474 
475  for( typename std::vector< Point<VectorType> >::iterator it = front.begin(); it != front.end(); ++it )
476  activePoints.push_back( it );
477 
479  bbc( front.begin(), front.end() );
480 
481  double maxBoundingBoxVolume = 0.;
482 
483  for( typename std::vector< Point<VectorType> >::iterator it = front.begin(); it != front.end(); ++it )
484  maxBoundingBoxVolume = std::max( maxBoundingBoxVolume, it->m_boundingBoxVolume );
485 
486  double _delta = m_startDelta * maxBoundingBoxVolume;
487  m_logFactor = ::log( 2. * front.size() * (1. + m_gamma) / (m_errorProbability * m_gamma) );
488  double minApprox = std::numeric_limits<double>::max();
489 
490  typename std::vector<Point< VectorType > >::iterator minimalElement;
491  m_round = 0;
492  m_sampleCounter = 0;
493 
494  while( activePoints.size() > 1 ) {
495  m_round++;
496 
497  minApprox = std::numeric_limits<double>::max();
498  minimalElement = front.end();
499 
500  for( int i = 0; i < static_cast<int>( activePoints.size() ); i++ )
501  try {
502  sample( front, activePoints[i], m_round, _delta, refPoint );
503  } catch( ... ) {
504  std::size_t idx = std::distance( front.begin(), minimalElement );
505  return idx;
506  }
507 
508  for( typename std::vector< typename std::vector< Point<VectorType> >::iterator >::iterator it = activePoints.begin(); it != activePoints.end(); ++it ) {
509  if( (*it)->m_approximatedContribution < minApprox ) {
510  minApprox = (*it)->m_approximatedContribution;
511  minimalElement = *it;
512  }
513  }
514 
515  if( activePoints.size() > 2 ) {
516  try {
517  sample( front, minimalElement, m_round, m_minimumMultiplierDelta * _delta, refPoint );
518  } catch( ... ) {
519  std::size_t idx = std::distance( front.begin(), minimalElement );
520  return idx;
521  }
522  }
523 
524  minApprox = std::numeric_limits<double>::max();
525  for( typename std::vector< typename std::vector< Point<VectorType> >::iterator >::iterator it = activePoints.begin(); it != activePoints.end(); ++it ) {
526  if( (*it)->m_approximatedContribution < minApprox ) {
527  minApprox = (*it)->m_approximatedContribution;
528  minimalElement = *it;
529  }
530  }
531 
532  typename std::vector< typename std::vector< Point<VectorType> >::iterator >::iterator it = activePoints.begin();
533  while( it != activePoints.end() ) {
534  if( (*it)->m_approximatedContribution - minApprox > deltaForPoint( *it, m_round ) + deltaForPoint( minimalElement, m_round ) )
535  it = activePoints.erase( it );
536  else
537  ++it;
538  }
539 
540  if( activePoints.size() <= 1 )
541  break;
542 
543  double d = 0;
544  for( it = activePoints.begin(); it != activePoints.end(); ++it ) {
545  if( *it == minimalElement )
546  continue;
547  double nom = minApprox + deltaForPoint( minimalElement, m_round );
548  double den = (*it)->m_approximatedContribution - deltaForPoint( *it, m_round );
549 
550  if( den <= 0. )
552  else if( d < nom/den )
553  d = nom/den;
554  }
555 
556  if( d < 1. + m_errorBound )
557  break;
558 
559  _delta *= m_multiplierDelta;
560  }
561 
562  return std::distance( front.begin(), minimalElement );
563  }
564 
565  /**
566  * \brief Determines the point contributing the least hypervolume to the overall front of points.
567  *
568  * This version uses the reference point estimated by the last call to updateInternals.
569  *
570  * \param [in] extractor Extracts point information from front elements.
571  * \param [in] front pareto front of points
572  */
573  template<typename Extractor, typename ParetoFrontType>
574  std::size_t leastContributor( Extractor extractor, ParetoFrontType const& front)
575  {
576  return leastContributor(extractor,front,m_reference);
577  }
578 
579  /// \brief Updates the internal variables of the indicator using a whole population.
580  ///
581  /// Calculates the reference point of the volume from the population
582  /// using the maximum value in every dimension+1
583  /// \param extractor Extracts point information from set.
584  /// \param set The set of points.
585  template<typename Extractor, typename PointSet>
586  void updateInternals(Extractor extractor, PointSet const& set){
587  m_reference.clear();
588  if(set.empty()) return;
589 
590  //calculate reference point
591  std::size_t noObjectives = extractor(set[0]).size();
592  m_reference.resize(noObjectives);
593 
594  for( unsigned int i = 0; i < set.size(); i++ )
595  noalias(m_reference) = max(m_reference, extractor(set[i]));
596 
597  noalias(m_reference)+= 1.0;
598  }
599 };
600 }
601 
602 #endif