HypervolumeCalculator.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Implementation of the exact hypervolume calculation in m dimensions.
5  *
6  * The algorithm is described in
7  *
8  * Nicola Beume und Günter Rudolph.
9  * Faster S-Metric Calculation by Considering Dominated Hypervolume as Klee's Measure Problem.
10  * In: B. Kovalerchuk (ed.): Proceedings of the Second IASTED Conference on Computational Intelligence (CI 2006),
11  * pp. 231-236. ACTA Press: Anaheim, 2006.
12  *
13  * A specialized algorithm is used for the three-objective case:
14  *
15  * M. T. M. Emmerich and C. M. Fonseca.
16  * Computing hypervolume contributions in low dimensions: Asymptotically optimal algorithm and complexity results.
17  * In: Evolutionary Multi-Criterion Optimization (EMO) 2011.
18  * Vol. 6576 of Lecture Notes in Computer Science, pp. 121--135, Berlin: Springer, 2011.
19  *
20  *
21  * \author T.Voss, O.Krause, T. Glasmachers
22  * \date 2014-2016
23  *
24  *
25  * \par Copyright 1995-2016 Shark Development Team
26  *
27  * <BR><HR>
28  * This file is part of Shark.
29  * <http://image.diku.dk/shark/>
30  *
31  * Shark is free software: you can redistribute it and/or modify
32  * it under the terms of the GNU Lesser General Public License as published
33  * by the Free Software Foundation, either version 3 of the License, or
34  * (at your option) any later version.
35  *
36  * Shark is distributed in the hope that it will be useful,
37  * but WITHOUT ANY WARRANTY; without even the implied warranty of
38  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
39  * GNU Lesser General Public License for more details.
40  *
41  * You should have received a copy of the GNU Lesser General Public License
42  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
43  *
44  */
45 #ifndef SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H
46 #define SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H
47 
48 #include <shark/LinAlg/Base.h>
49 
50 #include <algorithm>
51 #include <iostream>
52 #include <vector>
53 #include <map>
54 #include <cmath>
55 
56 namespace shark {
57  /**
58  * \brief Implementation of the exact hypervolume calculation in m dimensions.
59  *
60  * The algorithm is described in
61  *
62  * Nicola Beume und Günter Rudolph.
63  * Faster S-Metric Calculation by Considering Dominated Hypervolume as Klee's Measure Problem.
64  * In: B. Kovalerchuk (ed.): Proceedings of the Second IASTED Conference on Computational Intelligence (CI 2006),
65  * pp. 231-236. ACTA Press: Anaheim, 2006.
66  */
68 
69  /**
70  * \brief Returns an estimate on the runtime of the algorithm.
71  * \param [in] noPoints The number of points n considered in the runtime estimation.
72  * \param [in] noObjectives The number of objectives m considered in the runtime estimation.
73  */
74  static double runtime( std::size_t noPoints, std::size_t noObjectives ) {
75  if( noPoints < 10 )
76  return HypervolumeCalculator::runtime( 10, noObjectives );
77  return( 0.03 * noObjectives * noObjectives * ::exp( ::log( static_cast<double>( noPoints ) ) * noObjectives * 0.5 ) );
78  }
79 
80  /**
81  * \brief Default c'tor.
82  */
83  HypervolumeCalculator() : m_useLogHyp( false ) {
84  }
85 
86  /**
87  * \brief Serializes/Deserializes the state of the calculator to the supplied archive.
88  * \tparam Archive Archive type, needs to be a model of a boost::serialization archive.
89  * \param [in,out] archive Archive to store to/load from.
90  * \param [in] version Currently unused.
91  */
92  template<typename Archive>
93  void serialize( Archive & archive, const std::size_t version ) {
94  archive & BOOST_SERIALIZATION_NVP( m_noObjectives );
95  archive & BOOST_SERIALIZATION_NVP( m_sqrtNoPoints );
96  archive & BOOST_SERIALIZATION_NVP( m_useLogHyp );
97  }
98 
99  /**
100  * \brief Executes the algorithm.
101  * \param [in] extractor Function object \f$f\f$to "project" elements of the set to \f$\mathbb{R}^m\f$.
102  * \param [in] set The set \f$S\f$ of points for which the following assumption needs to hold: \f$\forall s \in S: \lnot \exists s' \in S: f( s' ) \preceq f( s ) \f$
103  * \param [in] refPoint The reference point \f$\vec{r} \in \mathbb{R}^m\f$ for the hypervolume calculation, needs to fulfill: \f$ \forall s \in S: s \preceq \vec{r}\f$. .
104  */
105  template<typename Set,typename Extractor, typename VectorType>
106  double operator()( Extractor const& extractor, const Set & set, const VectorType & refPoint);
107 
108  /** \cond IMPL */
109  template<typename VectorType>
110  int covers( const VectorType & cuboid, const VectorType & regionLow );
111 
112  template<typename VectorType>
113  int partCovers( const VectorType & cuboid, const VectorType & regionUp );
114 
115  template<typename VectorType>
116  int containsBoundary ( const VectorType & cub, const VectorType & regLow, int split );
117 
118  template<typename VectorType>
119  double getMeasure( const VectorType & regionLow, const VectorType & regionUp );
120 
121  template<typename VectorType>
122  int isPile( const VectorType & cuboid, const VectorType & regionLow, const VectorType & regionUp );
123 
124  template<typename VectorType>
125  unsigned int binaryToInt( const VectorType & bs );
126 
127  template<typename VectorType>
128  void intToBinary( unsigned int i, VectorType & result );
129 
130  template<typename VectorType>
131  double computeTrellis( const VectorType & regLow, const VectorType & regUp, const VectorType & trellis );
132 
133  template<typename VectorType>
134  double getMedian( const VectorType & bounds, int length );
135 
136  template<typename Set, typename Extractor, typename VectorType>
137  double stream ( const VectorType & regionLow,
138  const VectorType & regionUp,
139  const Set & points,
140  Extractor const& extractor,
141  int split,
142  double cover
143  );
144 
145  std::size_t m_noObjectives;
146  std::size_t m_sqrtNoPoints;
147  bool m_useLogHyp;
148  template<typename Extractor>
149  struct LastObjectiveComparator {
150 
151  LastObjectiveComparator( Extractor const& extractor ) : m_extractor( extractor ) {}
152 
153  template<typename VectorType>
154  bool operator()( const VectorType & lhs, const VectorType & rhs ) {
155  return( m_extractor( lhs ).back() < m_extractor( rhs ).back() );
156  }
157 
158  Extractor m_extractor;
159  };
160  /** \endcond IMPL */
161  };
162 
163  /** \cond IMPL */
164  template<typename Set,typename Extractor, typename VectorType >
165  double HypervolumeCalculator::operator()( Extractor const& extractor, const Set & constSet, const VectorType & refPoint) {
166 
167  m_noObjectives = extractor(*constSet.begin()).size();
168  m_sqrtNoPoints = static_cast< std::size_t >( ::sqrt( static_cast<double>( constSet.size() ) ) );
169 
170  Set set( constSet );
171 
172  std::stable_sort( set.begin(), set.end(), LastObjectiveComparator<Extractor>( extractor ) );
173 
174  if( m_noObjectives == 2 ) {
175 
176  double h;
177  if( m_useLogHyp )
178  h = ( ::log( refPoint[0] ) - ::log( extractor( set[0] )[0] ) ) * (::log( refPoint[1] ) - ::log( extractor( set[0] )[1] ) );
179  else
180  h = ( refPoint[0] - extractor( set[0] )[0] ) * ( refPoint[1] - extractor( set[0] )[1] );
181 
182  double diffDim1; std::size_t lastValidIndex = 0;
183  for( std::size_t i = 1; i < set.size(); i++ ) {
184  if( m_useLogHyp )
185  diffDim1 = ::log( extractor( set[lastValidIndex] )[0] ) - ::log( extractor( set[i] )[0] ); // Might be negative, if the i-th solution is dominated.
186  else
187  diffDim1 = extractor( set[lastValidIndex] )[0] - extractor( set[i] )[0];
188 
189  if( diffDim1 > 0 ) {
190  if( m_useLogHyp )
191  h += diffDim1 * ( ::log( refPoint[1] ) - ::log( extractor( set[i] )[1] ) );
192  else
193  h += ( diffDim1 ) * ( refPoint[1] - extractor( set[i] )[1] );
194  lastValidIndex = i;
195  }
196  }
197  return h;
198  }
199  else if (m_noObjectives == 3)
200  {
201  if (m_useLogHyp)
202  {
203  double volume = 0.0;
204  double area = 0.0;
205  std::map<double, double> front2D;
206  double prev_x2 = 0.0;
207  for (size_t i=0; i<set.size(); i++)
208  {
209  VectorType x = extractor(set[i]);
210  if (i > 0) volume += area * (std::log(x[2]) - prev_x2);
211 
212  // check whether x is dominated
213  std::map<double, double>::iterator worse = front2D.upper_bound(std::log(x[0]));
214  double b = std::log(refPoint[1]);
215  if (worse != front2D.begin())
216  {
217  std::map<double, double>::iterator better = worse;
218  if (better == front2D.end() || better->first > std::log(x[0])) --better;
219  if (better->second <= std::log(x[1])) continue;
220  b = better->second;
221  }
222 
223  // remove dominated points
224  while (worse != front2D.end())
225  {
226  if (worse->second < std::log(x[1])) break;
227  std::map<double, double>::iterator it = worse;
228  ++worse;
229  double r = (worse == front2D.end()) ? std::log(refPoint[0]): worse->first;
230  area -= (r - it->first) * (b - it->second);
231  front2D.erase(it);
232  }
233 
234  // insert x
235  front2D[std::log(x[0])] = std::log(x[1]);
236  double r = (worse == front2D.end()) ? std::log(refPoint[0]) : worse->first;
237  area += (r - std::log(x[0])) * (b - std::log(x[1]));
238 
239  prev_x2 = std::log(x[2]);
240  }
241  volume += area * (std::log(refPoint[2]) - prev_x2);
242  return volume;
243  }
244  else
245  {
246  double volume = 0.0;
247  double area = 0.0;
248  std::map<double, double> front2D;
249  double prev_x2 = 0.0;
250  for (size_t i=0; i<set.size(); i++)
251  {
252  VectorType x = extractor(set[i]);
253  if (i > 0) volume += area * (x[2] - prev_x2);
254 
255  // check whether x is dominated
256  std::map<double, double>::iterator worse = front2D.upper_bound(x[0]);
257  double b = refPoint[1];
258  if (worse != front2D.begin())
259  {
260  std::map<double, double>::iterator better = worse;
261  if (better == front2D.end() || better->first > x[0]) --better;
262  if (better->second <= x[1]) continue;
263  b = better->second;
264  }
265 
266  // remove dominated points
267  while (worse != front2D.end())
268  {
269  if (worse->second < x[1]) break;
270  std::map<double, double>::iterator it = worse;
271  ++worse;
272  double r = (worse == front2D.end()) ? refPoint[0]: worse->first;
273  area -= (r - it->first) * (b - it->second);
274  front2D.erase(it);
275  }
276 
277  // insert x
278  front2D[x[0]] = x[1];
279  double r = (worse == front2D.end()) ? refPoint[0] : worse->first;
280  area += (r - x[0]) * (b - x[1]);
281 
282  prev_x2 = x[2];
283  }
284  volume += area * (refPoint[2] - prev_x2);
285  return volume;
286  }
287  }
288  else
289  {
290  VectorType regLow( m_noObjectives, 1E15 );
291  for( std::size_t i = 0; i < set.size(); i++ ){
292  noalias(regLow) = min(regLow,extractor(set[i]));
293  }
294  return( stream( regLow, refPoint, set, extractor, 0, refPoint.back() ) );
295  }
296  }
297 
298  template<typename VectorType>
299  int HypervolumeCalculator::covers( const VectorType & cuboid, const VectorType & regionLow ) {
300  for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
301  // for( std::size_t i = 0; i < std::min( cuboid.size(), regionLow.size() ); i++ ) {
302  if( cuboid[i] > regionLow[i] )
303  return (0);
304  }
305  return (1);
306  }
307 
308  template<typename VectorType>
309  int HypervolumeCalculator::partCovers( const VectorType & cuboid, const VectorType & regionUp ) {
310  // for( std::size_t i = 0; i < std::min( cuboid.size(), regionUp.size() ); i++) {
311  for( std::size_t i = 0; i < m_noObjectives-1; i++) {
312  if (cuboid[i] >= regionUp[i])
313  return (0);
314  }
315  return (1);
316  }
317 
318  template<typename VectorType>
319  int HypervolumeCalculator::containsBoundary( const VectorType & cub, const VectorType & regLow, int split ) {
320  if( !( regLow[split] < cub[split] ) ) {
321  return -1;
322  } else {
323  for ( int j = 0; j < split; j++) {
324  if (regLow[j] < cub[j]) {
325  return 1;
326  }
327  }
328  }
329  return 0;
330  }
331 
332  template<typename VectorType>
333  double HypervolumeCalculator::getMeasure( const VectorType & regionLow, const VectorType & regionUp ) {
334  double volume = 1.0;
335  // for ( std::size_t i = 0; i < regionLow.size(); i++ ) {
336  for( std::size_t i = 0; i < m_noObjectives-1; i++) {
337  volume *= (regionUp[i] - regionLow[i]);
338  }
339 
340  // std::cout << "Get Measure: " << volume << std::endl;
341 
342  return( volume );
343  }
344 
345  template<typename VectorType>
346  int HypervolumeCalculator::isPile( const VectorType & cuboid, const VectorType & regionLow, const VectorType & regionUp ) {
347  std::size_t pile = cuboid.size();
348  // for( std::size_t i = 0; i < NO_OBJECTIVES - 1; i++) {
349  for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
350  if( cuboid[i] > regionLow[i] ) {
351  if( pile != m_noObjectives ) {
352  return (-1);
353  }
354 
355  pile = i;
356  }
357  }
358 
359  return (int)pile;
360  }
361 
362  template<typename VectorType>
363  unsigned int HypervolumeCalculator::binaryToInt( const VectorType & v ) {
364  int result = 0;
365  unsigned i;
366  for (i = 0; i < v.size(); i++) {
367  result += v[i] ? ( 1 << i ) : 0;
368  }
369 
370  return (result);
371  }
372 
373  template<typename VectorType>
374  void HypervolumeCalculator::intToBinary(unsigned int i, VectorType & result) {
375  for (std::size_t j = 0; j < m_noObjectives - 1; j++)
376  result[j] = 0;
377 
378  unsigned int rest = i;
379  std::size_t idx = 0;
380 
381  while (rest != 0) {
382  result[idx] = (rest % 2);
383 
384  rest = rest / 2;
385  idx++;
386  }
387  }
388 
389  template<typename VectorType>
390  double HypervolumeCalculator::computeTrellis( const VectorType & regLow, const VectorType & regUp, const VectorType & trellis ) {
391  std::vector<int> bs( m_noObjectives-1, 1 );
392 
393  double result = 0;
394 
395  unsigned int noSummands = binaryToInt(bs);
396  int oneCounter; double summand;
397 
398  for(unsigned i = 1; i <= noSummands; i++ ) {
399  summand = 1;
400  intToBinary(i, bs);
401  oneCounter = 0;
402 
403  for(std::size_t j = 0; j < m_noObjectives-1; j++ ) {
404  if (bs[j] == 1) {
405  summand *= regUp[j] - trellis[j];
406  oneCounter++;
407  } else
408  summand *= regUp[j] - regLow[j];
409  }
410 
411  if (oneCounter % 2 == 0)
412  result -= summand ;
413  else
414  result += summand;
415  }
416 
417  return(result);
418  }
419 
420  template<typename VectorType>
421  double HypervolumeCalculator::getMedian( const VectorType & bounds, int length) {
422  if( length == 1 ) {
423  return bounds[0];
424  } else if( length == 2 ) {
425  return bounds[1];
426  }
427 
428  VectorType v( length );
429  std::copy( bounds.begin(), bounds.begin() + length, v.begin() );
430  std::sort( v.begin(), v.end() );
431 
432  return(length % 2 == 1 ? v[length/2] : (v[length/2-1] + v[(length/2)]) / 2);
433  }
434 
435  template<typename Set, typename Extractor, typename VectorType>
436  double HypervolumeCalculator::stream( const VectorType & regionLow,
437  const VectorType & regionUp,
438  const Set & points,
439  Extractor const& extractor,
440  int split,
441  double cover ) {
442  double coverOld;
443  coverOld = cover;
444  int coverIndex = 0;
445  int coverIndexOld = -1;
446  int c;
447 
448  double result = 0;
449 
450  double dMeasure = getMeasure(regionLow, regionUp);
451  while( cover == coverOld && coverIndex < static_cast<int>( points.size() ) ) {
452  if( coverIndexOld == coverIndex )
453  break;
454 
455  coverIndexOld = coverIndex;
456 
457  if( covers( extractor( points[coverIndex] ), regionLow) ) {
458  cover = extractor( points[coverIndex] )[m_noObjectives-1]; //points[coverIndex * NO_OBJECTIVES + NO_OBJECTIVES - 1];
459  result += dMeasure * (coverOld - cover);
460  }
461  else
462  coverIndex++;
463  }
464 
465  // std::cout << "(II) No points: " << points.size() << std::endl;
466 
467  for (c = coverIndex; c > 0; c--) {
468  // if (points[(c - 1) * NO_OBJECTIVES + NO_OBJECTIVES - 1] == cover) {
469  /*std::cout << "\t " << extractor( points[c-1] )[m_noObjectives-1] << " vs. " << cover << std::endl;
470  if( c == 1 )
471  std::cout << "\t\t " << extractor( points[c-1] )[m_noObjectives-1] << " vs. " << cover << std::endl;*/
472  if( extractor( points[c-1] )[m_noObjectives-1] == cover) {
473  coverIndex--;
474  }
475  }
476 
477  // std::cout << "Cover index: " << coverIndex << ", split: " << split << std::endl;
478 
479  if (coverIndex == 0)
480  return (result);
481 
482  bool allPiles = true; // int i;
483 
484  std::vector<int> piles( coverIndex );
485  // int * piles = (int*)malloc(coverIndex * sizeof(int));
486 
487  for( int i = 0; i < coverIndex; i++ ) {
488  piles[i] = isPile( extractor( points[i] ), regionLow, regionUp );//isPile(points + i * NO_OBJECTIVES, regionLow, regionUp);
489  if (piles[i] == -1) {
490  allPiles = false;
491  break;
492  }
493  }
494 
495  if( allPiles ) {
496  VectorType trellis( regionUp );
497 
498  double current = 0.0;
499  double next = 0.0;
500  int i = 0;
501  do {
502  current = extractor( points[i] )[m_noObjectives-1]; // [m_noObjectives-1]; // points[i * NO_OBJECTIVES + NO_OBJECTIVES - 1];
503  do {
504  if( extractor( points[i] )[piles[i]] < trellis[piles[i]] ) {
505  trellis[piles[i]] = extractor( points[i] )[piles[i]];//points[i * NO_OBJECTIVES + piles[i]];
506  }
507  i++;
508  if (i < coverIndex) {
509  next = extractor( points[i] )[m_noObjectives-1];// points[i * NO_OBJECTIVES + NO_OBJECTIVES - 1];
510  }
511  else {
512  next = cover;
513  break;
514  }
515 
516  }
517  while (next == current);
518 
519  result += computeTrellis(regionLow, regionUp, trellis) * (next - current);
520  } while (next != cover);
521 
522  } else {
523  double bound = -1.0;
524  std::vector<double> boundaries( coverIndex );
525  std::vector<double> noBoundaries( coverIndex );
526  unsigned boundIdx = 0;
527  unsigned noBoundIdx = 0;
528 
529  do {
530  for( int i = 0; i < coverIndex; i++ ) {
531  int contained = containsBoundary( extractor( points[i] ), regionLow, split ); // containsBoundary(points + i * NO_OBJECTIVES, regionLow, split);
532  if (contained == 1) {
533  // boundaries.push_back( extractor( points[i] )[split] ); //points[i * NO_OBJECTIVES + split];
534  boundaries[boundIdx] = extractor( points[i] )[split];
535  boundIdx++;
536  } else if (contained == 0) {
537  // noBoundaries.push_back( extractor( points[i] )[split] ); // points[i * NO_OBJECTIVES + split];
538  noBoundaries[noBoundIdx] = extractor( points[i] )[split];
539  noBoundIdx++;
540  }
541  }
542 
543  if (boundIdx > 0) {
544  bound = getMedian( boundaries, boundIdx );
545  // std::cout << "Median: " << bound << std::endl;
546  } else if( noBoundIdx > m_sqrtNoPoints ) {
547  bound = getMedian( noBoundaries, noBoundIdx );
548  // std::cout << "Median: " << bound << std::endl;
549  } else {
550  split++;
551  }
552  } while (bound == -1.0);
553 
554  Set pointsChildLow, pointsChildUp;
555  // pointsChildLow.reserve( coverIndex );
556  // pointsChildUp.reserve( coverIndex );
557 
558  VectorType regionUpC( regionUp );
559  regionUpC[split] = bound;
560  VectorType regionLowC( regionLow );
561  regionLowC[split] = bound;
562 
563  for( int i = 0; i < coverIndex; i++) {
564  if( partCovers( extractor( points[i] ), regionUpC) ) {
565  // pointsChildUp.append( points[i] );
566  pointsChildUp.push_back( points[i] );
567  }
568 
569  if( partCovers( extractor( points[i] ), regionUp ) ) {
570  // pointsChildLow.append( points[i] );
571 
572  pointsChildLow.push_back( points[i] );
573  }
574  }
575 
576  // #pragma omp sections
577  {
578 
579  // #pragma omp task
580  {
581  if( pointsChildUp.size() > 0 ) {
582  // #pragma omp atomic
583  result += stream(regionLow, regionUpC, pointsChildUp, extractor, split, cover);
584  }
585  }
586 
587  // #pragma omp task
588  {
589  if (pointsChildLow.size() > 0) {
590  // #pragma omp atomic
591  result += stream(regionLowC, regionUp, pointsChildLow, extractor, split, cover);
592  }
593  }
594 
595  }
596  }
597 
598  return (result);
599  }
600  /** \endcond IMPL */
601 }
602 #endif