39 #ifndef SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H 40 #define SHARK_ALGORITHMS_DIRECTSEARCH_HYPERVOLUMECALCULATOR_H 68 static double runtime( std::size_t noPoints, std::size_t noObjectives ) {
71 return( 0.03 * noObjectives * noObjectives * ::exp( ::log( static_cast<double>( noPoints ) ) * noObjectives * 0.5 ) );
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 );
99 template<
typename Set,
typename Extractor,
typename VectorType>
103 template<
typename VectorType>
106 template<
typename VectorType>
109 template<
typename VectorType>
112 template<
typename VectorType>
115 template<
typename VectorType>
118 template<
typename VectorType>
119 unsigned int binaryToInt(
const VectorType & bs );
121 template<
typename VectorType>
122 void intToBinary(
unsigned int i,
VectorType & result );
124 template<
typename VectorType>
127 template<
typename VectorType>
128 double getMedian(
const VectorType & bounds,
int length );
130 template<
typename Set,
typename Extractor,
typename VectorType>
134 Extractor
const& extractor,
139 std::size_t m_noObjectives;
140 std::size_t m_sqrtNoPoints;
142 template<
typename Extractor>
143 struct LastObjectiveComparator {
145 LastObjectiveComparator( Extractor
const& extractor ) : m_extractor( extractor ) {}
147 template<
typename VectorType>
149 return( m_extractor( lhs ).back() < m_extractor( rhs ).back() );
152 Extractor m_extractor;
158 template<
typename Set,
typename Extractor,
typename VectorType >
161 m_noObjectives = extractor(*constSet.begin()).
size();
162 m_sqrtNoPoints =
static_cast< std::size_t
>( ::sqrt( static_cast<double>( constSet.size() ) ) );
166 std::stable_sort(
set.begin(),
set.end(), LastObjectiveComparator<Extractor>( extractor ) );
168 if( m_noObjectives == 2 ) {
172 h = ( ::log( refPoint[0] ) - ::log( extractor(
set[0] )[0] ) ) * (::log( refPoint[1] ) - ::log( extractor(
set[0] )[1] ) );
174 h = ( refPoint[0] - extractor(
set[0] )[0] ) * ( refPoint[1] - extractor(
set[0] )[1] );
176 double diffDim1; std::size_t lastValidIndex = 0;
177 for( std::size_t i = 1; i <
set.size(); i++ ) {
179 diffDim1 = ::log( extractor(
set[lastValidIndex] )[0] ) - ::log( extractor(
set[i] )[0] );
181 diffDim1 = extractor(
set[lastValidIndex] )[0] - extractor(
set[i] )[0];
185 h += diffDim1 * ( ::log( refPoint[1] ) - ::log( extractor(
set[i] )[1] ) );
187 h += ( diffDim1 ) * ( refPoint[1] - extractor(
set[i] )[1] );
195 for( std::size_t i = 0; i <
set.size(); i++ ){
196 noalias(regLow) =
min(regLow,extractor(
set[i]));
198 return( stream( regLow, refPoint,
set, extractor, 0, refPoint.back() ) );
201 template<
typename VectorType>
203 for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
205 if( cuboid[i] > regionLow[i] )
211 template<
typename VectorType>
212 int HypervolumeCalculator::partCovers(
const VectorType & cuboid,
const VectorType & regionUp ) {
214 for( std::size_t i = 0; i < m_noObjectives-1; i++) {
215 if (cuboid[i] >= regionUp[i])
221 template<
typename VectorType>
222 int HypervolumeCalculator::containsBoundary(
const VectorType & cub,
const VectorType & regLow,
int split ) {
223 if( !( regLow[split] < cub[split] ) ) {
226 for (
int j = 0; j < split; j++) {
227 if (regLow[j] < cub[j]) {
235 template<
typename VectorType>
236 double HypervolumeCalculator::getMeasure(
const VectorType & regionLow,
const VectorType & regionUp ) {
239 for( std::size_t i = 0; i < m_noObjectives-1; i++) {
240 volume *= (regionUp[i] - regionLow[i]);
248 template<
typename VectorType>
250 std::size_t pile = cuboid.size();
252 for( std::size_t i = 0; i < m_noObjectives-1; i++ ) {
253 if( cuboid[i] > regionLow[i] ) {
254 if( pile != m_noObjectives ) {
265 template<
typename VectorType>
266 unsigned int HypervolumeCalculator::binaryToInt(
const VectorType & v ) {
269 for (i = 0; i < v.size(); i++) {
270 result += v[i] ? ( 1 << i ) : 0;
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++)
281 unsigned int rest = i;
285 result[idx] = (rest % 2);
292 template<
typename VectorType>
294 std::vector<int> bs( m_noObjectives-1, 1 );
298 unsigned int noSummands = binaryToInt(bs);
299 int oneCounter;
double summand;
301 for(
unsigned i = 1; i <= noSummands; i++ ) {
306 for(std::size_t j = 0; j < m_noObjectives-1; j++ ) {
308 summand *= regUp[j] - trellis[j];
311 summand *= regUp[j] - regLow[j];
314 if (oneCounter % 2 == 0)
323 template<
typename VectorType>
324 double HypervolumeCalculator::getMedian(
const VectorType & bounds,
int length) {
327 }
else if( length == 2 ) {
332 std::copy( bounds.begin(), bounds.begin() + length, v.begin() );
333 std::sort( v.begin(), v.end() );
335 return(length % 2 == 1 ? v[length/2] : (v[length/2-1] + v[(length/2)]) / 2);
338 template<
typename Set,
typename Extractor,
typename VectorType>
339 double HypervolumeCalculator::stream(
const VectorType & regionLow,
342 Extractor
const& extractor,
348 int coverIndexOld = -1;
353 double dMeasure = getMeasure(regionLow, regionUp);
354 while( cover == coverOld && coverIndex < static_cast<int>( points.size() ) ) {
355 if( coverIndexOld == coverIndex )
358 coverIndexOld = coverIndex;
360 if( covers( extractor( points[coverIndex] ), regionLow) ) {
361 cover = extractor( points[coverIndex] )[m_noObjectives-1];
362 result += dMeasure * (coverOld - cover);
370 for (c = coverIndex; c > 0; c--) {
375 if( extractor( points[c-1] )[m_noObjectives-1] == cover) {
385 bool allPiles =
true;
387 std::vector<int> piles( coverIndex );
390 for(
int i = 0; i < coverIndex; i++ ) {
391 piles[i] = isPile( extractor( points[i] ), regionLow, regionUp );
392 if (piles[i] == -1) {
401 double current = 0.0;
405 current = extractor( points[i] )[m_noObjectives-1];
407 if( extractor( points[i] )[piles[i]] < trellis[piles[i]] ) {
408 trellis[piles[i]] = extractor( points[i] )[piles[i]];
411 if (i < coverIndex) {
412 next = extractor( points[i] )[m_noObjectives-1];
420 while (next == current);
422 result += computeTrellis(regionLow, regionUp, trellis) * (next - current);
423 }
while (next != cover);
427 std::vector<double> boundaries( coverIndex );
428 std::vector<double> noBoundaries( coverIndex );
429 unsigned boundIdx = 0;
430 unsigned noBoundIdx = 0;
433 for(
int i = 0; i < coverIndex; i++ ) {
434 int contained = containsBoundary( extractor( points[i] ), regionLow, split );
435 if (contained == 1) {
437 boundaries[boundIdx] = extractor( points[i] )[split];
439 }
else if (contained == 0) {
441 noBoundaries[noBoundIdx] = extractor( points[i] )[split];
447 bound = getMedian( boundaries, boundIdx );
449 }
else if( noBoundIdx > m_sqrtNoPoints ) {
450 bound = getMedian( noBoundaries, noBoundIdx );
455 }
while (bound == -1.0);
457 Set pointsChildLow, pointsChildUp;
462 regionUpC[split] = bound;
464 regionLowC[split] = bound;
466 for(
int i = 0; i < coverIndex; i++) {
467 if( partCovers( extractor( points[i] ), regionUpC) ) {
469 pointsChildUp.push_back( points[i] );
472 if( partCovers( extractor( points[i] ), regionUp ) ) {
475 pointsChildLow.push_back( points[i] );
484 if( pointsChildUp.size() > 0 ) {
486 result += stream(regionLow, regionUpC, pointsChildUp, extractor, split, cover);
492 if (pointsChildLow.size() > 0) {
494 result += stream(regionLowC, regionUp, pointsChildLow, extractor, split, cover);