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