45 #ifndef SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H 46 #define SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H 74 static double runtime( std::size_t noPoints, std::size_t noObjectives ) {
77 return( 0.03 * noObjectives * noObjectives * ::exp( ::log( static_cast<double>( noPoints ) ) * noObjectives * 0.5 ) );
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 );
105 template<
typename Set,
typename Extractor,
typename VectorType>
109 template<
typename VectorType>
112 template<
typename VectorType>
115 template<
typename VectorType>
118 template<
typename VectorType>
121 template<
typename VectorType>
124 template<
typename VectorType>
125 unsigned int binaryToInt(
const VectorType & bs );
127 template<
typename VectorType>
128 void intToBinary(
unsigned int i,
VectorType & result );
130 template<
typename VectorType>
133 template<
typename VectorType>
134 double getMedian(
const VectorType & bounds,
int length );
136 template<
typename Set,
typename Extractor,
typename VectorType>
140 Extractor
const& extractor,
145 std::size_t m_noObjectives;
146 std::size_t m_sqrtNoPoints;
148 template<
typename Extractor>
149 struct LastObjectiveComparator {
151 LastObjectiveComparator( Extractor
const& extractor ) : m_extractor( extractor ) {}
153 template<
typename VectorType>
155 return( m_extractor( lhs ).back() < m_extractor( rhs ).back() );
158 Extractor m_extractor;
164 template<
typename Set,
typename Extractor,
typename VectorType >
167 m_noObjectives = extractor(*constSet.begin()).
size();
168 m_sqrtNoPoints =
static_cast< std::size_t
>( ::sqrt( static_cast<double>( constSet.size() ) ) );
172 std::stable_sort(
set.begin(),
set.end(), LastObjectiveComparator<Extractor>( extractor ) );
174 if( m_noObjectives == 2 ) {
178 h = ( ::log( refPoint[0] ) - ::log( extractor(
set[0] )[0] ) ) * (::log( refPoint[1] ) - ::log( extractor(
set[0] )[1] ) );
180 h = ( refPoint[0] - extractor(
set[0] )[0] ) * ( refPoint[1] - extractor(
set[0] )[1] );
182 double diffDim1; std::size_t lastValidIndex = 0;
183 for( std::size_t i = 1; i <
set.size(); i++ ) {
185 diffDim1 = ::log( extractor(
set[lastValidIndex] )[0] ) - ::log( extractor(
set[i] )[0] );
187 diffDim1 = extractor(
set[lastValidIndex] )[0] - extractor(
set[i] )[0];
191 h += diffDim1 * ( ::log( refPoint[1] ) - ::log( extractor(
set[i] )[1] ) );
193 h += ( diffDim1 ) * ( refPoint[1] - extractor(
set[i] )[1] );
199 else if (m_noObjectives == 3)
205 std::map<double, double> front2D;
206 double prev_x2 = 0.0;
207 for (
size_t i=0; i<
set.size(); i++)
210 if (i > 0) volume += area * (std::log(x[2]) - prev_x2);
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())
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;
224 while (worse != front2D.end())
226 if (worse->second < std::log(x[1]))
break;
227 std::map<double, double>::iterator it = worse;
229 double r = (worse == front2D.end()) ? std::log(refPoint[0]): worse->first;
230 area -= (r - it->first) * (b - it->second);
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]));
239 prev_x2 = std::log(x[2]);
241 volume += area * (std::log(refPoint[2]) - prev_x2);
248 std::map<double, double> front2D;
249 double prev_x2 = 0.0;
250 for (
size_t i=0; i<
set.size(); i++)
253 if (i > 0) volume += area * (x[2] - prev_x2);
256 std::map<double, double>::iterator worse = front2D.upper_bound(x[0]);
257 double b = refPoint[1];
258 if (worse != front2D.begin())
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;
267 while (worse != front2D.end())
269 if (worse->second < x[1])
break;
270 std::map<double, double>::iterator it = worse;
272 double r = (worse == front2D.end()) ? refPoint[0]: worse->first;
273 area -= (r - it->first) * (b - it->second);
278 front2D[x[0]] = x[1];
279 double r = (worse == front2D.end()) ? refPoint[0] : worse->first;
280 area += (r - x[0]) * (b - x[1]);
284 volume += area * (refPoint[2] - prev_x2);
291 for( std::size_t i = 0; i <
set.size(); i++ ){
292 noalias(regLow) =
min(regLow,extractor(
set[i]));
294 return( stream( regLow, refPoint,
set, extractor, 0, refPoint.back() ) );
298 template<
typename VectorType>
300 for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
302 if( cuboid[i] > regionLow[i] )
308 template<
typename VectorType>
309 int HypervolumeCalculator::partCovers(
const VectorType & cuboid,
const VectorType & regionUp ) {
311 for( std::size_t i = 0; i < m_noObjectives-1; i++) {
312 if (cuboid[i] >= regionUp[i])
318 template<
typename VectorType>
319 int HypervolumeCalculator::containsBoundary(
const VectorType & cub,
const VectorType & regLow,
int split ) {
320 if( !( regLow[split] < cub[split] ) ) {
323 for (
int j = 0; j < split; j++) {
324 if (regLow[j] < cub[j]) {
332 template<
typename VectorType>
333 double HypervolumeCalculator::getMeasure(
const VectorType & regionLow,
const VectorType & regionUp ) {
336 for( std::size_t i = 0; i < m_noObjectives-1; i++) {
337 volume *= (regionUp[i] - regionLow[i]);
345 template<
typename VectorType>
347 std::size_t pile = cuboid.size();
349 for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
350 if( cuboid[i] > regionLow[i] ) {
351 if( pile != m_noObjectives ) {
362 template<
typename VectorType>
363 unsigned int HypervolumeCalculator::binaryToInt(
const VectorType & v ) {
366 for (i = 0; i < v.size(); i++) {
367 result += v[i] ? ( 1 << i ) : 0;
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++)
378 unsigned int rest = i;
382 result[idx] = (rest % 2);
389 template<
typename VectorType>
391 std::vector<int> bs( m_noObjectives-1, 1 );
395 unsigned int noSummands = binaryToInt(bs);
396 int oneCounter;
double summand;
398 for(
unsigned i = 1; i <= noSummands; i++ ) {
403 for(std::size_t j = 0; j < m_noObjectives-1; j++ ) {
405 summand *= regUp[j] - trellis[j];
408 summand *= regUp[j] - regLow[j];
411 if (oneCounter % 2 == 0)
420 template<
typename VectorType>
421 double HypervolumeCalculator::getMedian(
const VectorType & bounds,
int length) {
424 }
else if( length == 2 ) {
429 std::copy( bounds.begin(), bounds.begin() + length, v.begin() );
430 std::sort( v.begin(), v.end() );
432 return(length % 2 == 1 ? v[length/2] : (v[length/2-1] + v[(length/2)]) / 2);
435 template<
typename Set,
typename Extractor,
typename VectorType>
436 double HypervolumeCalculator::stream(
const VectorType & regionLow,
439 Extractor
const& extractor,
445 int coverIndexOld = -1;
450 double dMeasure = getMeasure(regionLow, regionUp);
451 while( cover == coverOld && coverIndex < static_cast<int>( points.size() ) ) {
452 if( coverIndexOld == coverIndex )
455 coverIndexOld = coverIndex;
457 if( covers( extractor( points[coverIndex] ), regionLow) ) {
458 cover = extractor( points[coverIndex] )[m_noObjectives-1];
459 result += dMeasure * (coverOld - cover);
467 for (c = coverIndex; c > 0; c--) {
472 if( extractor( points[c-1] )[m_noObjectives-1] == cover) {
482 bool allPiles =
true;
484 std::vector<int> piles( coverIndex );
487 for(
int i = 0; i < coverIndex; i++ ) {
488 piles[i] = isPile( extractor( points[i] ), regionLow, regionUp );
489 if (piles[i] == -1) {
498 double current = 0.0;
502 current = extractor( points[i] )[m_noObjectives-1];
504 if( extractor( points[i] )[piles[i]] < trellis[piles[i]] ) {
505 trellis[piles[i]] = extractor( points[i] )[piles[i]];
508 if (i < coverIndex) {
509 next = extractor( points[i] )[m_noObjectives-1];
517 while (next == current);
519 result += computeTrellis(regionLow, regionUp, trellis) * (next - current);
520 }
while (next != cover);
524 std::vector<double> boundaries( coverIndex );
525 std::vector<double> noBoundaries( coverIndex );
526 unsigned boundIdx = 0;
527 unsigned noBoundIdx = 0;
530 for(
int i = 0; i < coverIndex; i++ ) {
531 int contained = containsBoundary( extractor( points[i] ), regionLow, split );
532 if (contained == 1) {
534 boundaries[boundIdx] = extractor( points[i] )[split];
536 }
else if (contained == 0) {
538 noBoundaries[noBoundIdx] = extractor( points[i] )[split];
544 bound = getMedian( boundaries, boundIdx );
546 }
else if( noBoundIdx > m_sqrtNoPoints ) {
547 bound = getMedian( noBoundaries, noBoundIdx );
552 }
while (bound == -1.0);
554 Set pointsChildLow, pointsChildUp;
559 regionUpC[split] = bound;
561 regionLowC[split] = bound;
563 for(
int i = 0; i < coverIndex; i++) {
564 if( partCovers( extractor( points[i] ), regionUpC) ) {
566 pointsChildUp.push_back( points[i] );
569 if( partCovers( extractor( points[i] ), regionUp ) ) {
572 pointsChildLow.push_back( points[i] );
581 if( pointsChildUp.size() > 0 ) {
583 result += stream(regionLow, regionUpC, pointsChildUp, extractor, split, cover);
589 if (pointsChildLow.size() > 0) {
591 result += stream(regionLowC, regionUp, pointsChildLow, extractor, split, cover);