33 #ifndef SHARK_ALGORITHMS_DIRECT_SEARCH_INDICATORS_LEASTCONTRIBUTORAPPROXIMATOR_H 34 #define SHARK_ALGORITHMS_DIRECT_SEARCH_INDICATORS_LEASTCONTRIBUTORAPPROXIMATOR_H 39 #include <boost/assign.hpp> 40 #include <boost/bimap.hpp> 41 #include <boost/cstdint.hpp> 55 template<
typename Rng>
69 template<
typename Set>
70 bool operator()(
const Set & s,
typename Set::iterator point )
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 ] );
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;
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 ] ) {
102 template<
typename Set>
113 template<
typename Iterator>
116 return i1->m_overlappingVolume > i2->m_overlappingVolume;
125 void operator()(
typename Set::iterator begin,
typename Set::iterator end )
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 ) {
135 typename Set::value_type & p2 = *itt;
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 ] ) {
142 if( coordCounter == 2 )
148 if( coordCounter == 1 && p1.m_boundingBox[ coord ] > p2.m_point[ coord ] )
149 p1.m_boundingBox[ coord ] = p2.m_point[ coord ];
152 it->m_boundingBoxVolume = 1.;
153 for(
unsigned int i = 0 ; i < it->m_boundingBox.size(); i++ ) {
155 it->m_boundingBoxVolume *= it->m_boundingBox[ i ] - it->m_point[ i ];
159 for(
typename Set::iterator itt = begin; itt != end; ++itt ) {
163 bool isInfluencing =
true;
165 for(
unsigned int i = 0; i < it->m_point.size(); i++ ) {
167 if( itt->m_point[ i ] >= it->m_boundingBox[ i ] ) {
168 isInfluencing =
false;
171 vol *=
std::max( itt->m_point[ i ], it->m_point[ i ] ) - it->m_boundingBox[ i ];
173 if( isInfluencing ) {
174 itt->m_overlappingVolume = vol;
175 it->m_influencingPoints.push_back( itt );
179 std::sort( it->m_influencingPoints.begin(), it->m_influencingPoints.end(),
VolumeComparator() );
192 template<
typename Rng,
typename ExactHypervolume>
208 typedef boost::bimap< Strategy, std::string > registry_type;
211 static registry_type registry = boost::assign::list_of< typename registry_type::relation >( ENSURE_UNION_BOUND,
"EnsureUnionBound" )( FAST,
"Fast" );
217 stream << STRATEGY_REGISTRY().left.find( strategy )->second;
223 static std::string s;
225 strategy = STRATEGY_REGISTRY().right.find( s )->second;
288 template<
typename VectorType>
295 m_sample( point.
size() ),
296 m_boundingBox( refPoint ),
297 m_boundingBoxVolume( 0. ),
298 m_approximatedContribution( 0. ),
301 m_noSuccessfulSamples( 0 )
322 template<
typename Stream>
326 std::copy( m_point.begin(), m_point.end(), std::ostream_iterator<double>( s,
" " ) );
328 std::copy( m_boundingBox.begin(), m_boundingBox.end(), std::ostream_iterator<double>( s,
" " ) );
330 s << m_noOperations <<
" " << m_noSamples <<
" " << m_noSuccessfulSamples <<
" " << m_boundingBoxVolume << std::endl;
339 template<
typename Member>
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,
381 ) : m_startDelta( startDelta ),
382 m_multiplierDelta( multiplierDelta ),
383 m_minimumMultiplierDelta( minimumDeltaMultiplier ),
384 m_maxNumSamples( maxNumSamples ),
387 m_errorProbability(delta),
389 m_strategy( this_type::DEFAULT_STRATEGY() ),
390 m_sampleCounter( 0 ),
391 m_sampleCountThreshold( this_type::DEFAULT_SAMPLE_THRESHOLD() )
404 template<
typename Set,
typename VectorType>
405 void sample(
const Set & s,
typename Set::iterator point,
unsigned int r,
double delta,
const VectorType & refPoint )
408 if( point->m_noSamples >= m_maxNumSamples )
412 if( point->m_noOperations >= ExactHypervolume::runtime( point->m_influencingPoints.size(), refPoint.size() ) ) {
414 point->m_noSamples = point->m_noSuccessfulSamples = m_maxNumSamples;
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 ] );
424 point->m_approximatedContribution = point->m_boundingBoxVolume - hv( e, neighborhood, refPoint );
425 point->m_boundingBoxVolume = point->m_approximatedContribution;
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++;
442 point->m_approximatedContribution = point->m_boundingBoxVolume * (
static_cast<double>( point->m_noSuccessfulSamples ) / static_cast<double>( point->m_noSamples ) );
445 template<
typename iterator>
449 return( ::sqrt( 0.5 * ((1. + m_gamma) * ::log( static_cast<double>( R ) ) + m_logFactor ) / point->m_noSamples ) * point->m_boundingBoxVolume );
459 template<
typename Extractor,
typename Set,
typename VectorType>
462 if(s.empty())
return 0;
464 std::size_t noObjectives = e(s[0]).size();
466 std::vector< Point<VectorType> > front;
467 front.reserve( s.size() );
469 std::vector< typename std::vector< Point<VectorType> >::iterator > activePoints;
471 for(
typename Set::const_iterator it = s.begin(); it != s.end(); ++it ) {
475 for(
typename std::vector<
Point<VectorType> >::iterator it = front.begin(); it != front.end(); ++it )
476 activePoints.push_back( it );
479 bbc( front.begin(), front.end() );
481 double maxBoundingBoxVolume = 0.;
483 for(
typename std::vector<
Point<VectorType> >::iterator it = front.begin(); it != front.end(); ++it )
484 maxBoundingBoxVolume =
std::max( maxBoundingBoxVolume, it->m_boundingBoxVolume );
486 double _delta = m_startDelta * maxBoundingBoxVolume;
487 m_logFactor = ::log( 2. * front.size() * (1. + m_gamma) / (m_errorProbability * m_gamma) );
490 typename std::vector<Point< VectorType > >::iterator minimalElement;
494 while( activePoints.size() > 1 ) {
498 minimalElement = front.end();
500 for(
int i = 0; i < static_cast<int>( activePoints.size() ); i++ )
502 sample( front, activePoints[i], m_round, _delta, refPoint );
504 std::size_t idx =
std::distance( front.begin(), minimalElement );
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;
515 if( activePoints.size() > 2 ) {
517 sample( front, minimalElement, m_round, m_minimumMultiplierDelta * _delta, refPoint );
519 std::size_t idx =
std::distance( front.begin(), minimalElement );
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;
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 );
540 if( activePoints.size() <= 1 )
544 for( it = activePoints.begin(); it != activePoints.end(); ++it ) {
545 if( *it == minimalElement )
547 double nom = minApprox + deltaForPoint( minimalElement, m_round );
548 double den = (*it)->m_approximatedContribution - deltaForPoint( *it, m_round );
552 else if( d < nom/den )
556 if( d < 1. + m_errorBound )
559 _delta *= m_multiplierDelta;
573 template<
typename Extractor,
typename ParetoFrontType>
576 return leastContributor(extractor,front,m_reference);
585 template<
typename Extractor,
typename Po
intSet>
588 if(
set.empty())
return;
591 std::size_t noObjectives = extractor(
set[0]).size();
592 m_reference.resize(noObjectives);
594 for(
unsigned int i = 0; i <
set.size(); i++ )
595 noalias(m_reference) =
max(m_reference, extractor(
set[i]));