From 29816fc1974380edc5a30dc08702650e16a0887e Mon Sep 17 00:00:00 2001 From: maartenkeijzer Date: Tue, 4 Sep 2007 15:26:45 +0000 Subject: [PATCH] added epsilon moea --- eo/src/moo/eoEpsMOEA.h | 118 ++++++++++++++++++++++++++++ eo/src/moo/eoEpsilonArchive.h | 123 +++++++++++++++++++++++++++++ eo/src/moo/eoFrontSorter.h | 14 +--- eo/src/moo/eoMOFitness.h | 139 ++++++++++++++++++++++++--------- eo/src/moo/eoNSGA_IIa_Eval.cpp | 20 +---- 5 files changed, 349 insertions(+), 65 deletions(-) create mode 100644 eo/src/moo/eoEpsMOEA.h create mode 100644 eo/src/moo/eoEpsilonArchive.h diff --git a/eo/src/moo/eoEpsMOEA.h b/eo/src/moo/eoEpsMOEA.h new file mode 100644 index 00000000..f847a717 --- /dev/null +++ b/eo/src/moo/eoEpsMOEA.h @@ -0,0 +1,118 @@ +#ifndef eoEpsMOEA_h +#define eoEpsMOEA_h + +#include +#include + +template +class eoEpsMOEA : public eoAlgo { + + public: + + eoEpsMOEA( + eoContinue& _continuator, + eoEvalFunc& _eval, + eoGenOp& _op, + const std::vector& eps + ) : continuator(_continuator), + eval (_eval), + loopEval(_eval), + popEval(loopEval), + op(_op), + archive(eps) + {} + + + void operator()(eoPop& pop) { + + eoPop offspring; + popEval(offspring, pop); + for (unsigned i = 0; i < pop.size(); ++i) pop[i].fitnessReference().setWorth(1.0); + + do { + unsigned nProcessed = 0; + + while (nProcessed < pop.size()) { + offspring.clear(); + + epsPopulator populator(archive, pop, offspring); + op(populator); + + nProcessed += offspring.size(); + popEval(pop, offspring); + + for (unsigned i = 0; i < offspring.size(); ++i) { + offspring[i].fitnessReference().setWorth(1.0); + archive(offspring[i]); + update_pop(pop, offspring[i]); + } + } + + } while (continuator(pop)); + + // return archive + pop.clear(); + archive.appendTo(pop); + + } + + private : + void update_pop(eoPop& pop, const EOT& offspring) { + for (unsigned i = 0; i < pop.size(); ++i) { + int dom = offspring.fitness().check_dominance(pop[i].fitness()); + if (dom == 1) { + pop[i] = offspring; + return; + } + if (dom == -1) return; // + } + + // non-dominated everywhere, overwrite random one + pop[ rng.random(pop.size()) ] = offspring; + } + + class epsPopulator : public eoPopulator { + + using eoPopulator< EOT >::src; + + eoEpsilonArchive& archive; + bool fromArchive; + public: + epsPopulator(eoEpsilonArchive& arch, const eoPop& pop, eoPop& res) : eoPopulator(pop, res), archive(arch), fromArchive(true) {} + + const EOT& select() { + fromArchive = !fromArchive; + + using std::cout; + using std::endl; + + if (fromArchive && !archive.empty()) { + return archive.selectRandom(); + } + + // tournament selection on population + const EOT& eo1 = rng.choice(src); + const EOT& eo2 = rng.choice(src); + + if (eo1.fitness().dominates(eo2.fitness())) return eo1; + return eo2; // they are drawn at random, so no need to do an extra randomization step + } + + }; + + + eoContinue& continuator; + + eoEvalFunc & eval ; + eoPopLoopEval loopEval; + + eoPopEvalFunc& popEval; + + eoGenOp& op; + + eoEpsilonArchive archive; + +}; + +#endif + diff --git a/eo/src/moo/eoEpsilonArchive.h b/eo/src/moo/eoEpsilonArchive.h new file mode 100644 index 00000000..d20de936 --- /dev/null +++ b/eo/src/moo/eoEpsilonArchive.h @@ -0,0 +1,123 @@ +#ifndef eoEpsilonArchive_h +#define eoEpsilonArchive_h + +#include +#include + +template +class eoEpsilonArchive { + + typedef typename EOT::Fitness::fitness_traits Traits; + + struct Node { + EOT element; + std::vector discretized; + + Node(const EOT& eo) : element(eo), discretized(Traits::nObjectives()) {} + }; + + typedef std::list archive_t; + + archive_t archive; + std::vector eps; + + public: + + static double direction(unsigned i) { return 1; } + static double tol() { return 1e-6; } + + + eoEpsilonArchive(const std::vector& eps_) : eps(eps_) { + if (eps.size() != Traits::nObjectives()) throw std::logic_error("eoEpsilonArchive: need one epsilon for each objective"); + } + + bool empty() { return archive.size() == 0; } + + void operator()(const EOT& eo) { + + using std::cout; + using std::endl; + + // discretize + Node node(eo); + for (unsigned i = 0; i < eo.fitness().size(); ++i) { + double val = Traits::direction(i) * eo.fitness()[i]; + + node.discretized[i] = floor(val/eps[i]); + } + + using namespace dominance; + + typename archive_t::iterator boxIt = archive.end(); + // update archive + + for (typename archive_t::iterator it = archive.begin(); it != archive.end(); ++it) { + dominance_result result = check >(node.discretized, it->discretized); + + switch (result) { + case first : { // remove dominated archive member + it = archive.erase(it); + break; + } + case second : { + //cout << it->discretized[0] << ' ' << it->discretized[1] << " dominates " << node.discretized[0] << ' ' << node.discretized[1] << endl; + return; // new one does not belong in archive + } + case non_dominated : break; // both non-dominated, put in different boxes + case non_dominated_equal : { // in same box + boxIt = it; + goto exitLoop; // break + } + } + + } + +exitLoop: + // insert + if (boxIt == archive.end()) { // non-dominated, new box + archive.push_back(node); + } else { // fight it out + int dom = node.element.fitness().check_dominance( boxIt->element.fitness() ); + + switch (dom) { + case 1: *boxIt = node; break; + case -1: break; + case 0: { + double d1 = 0.0; + double d2 = 0.0; + + for (unsigned i = 0; i < node.element.fitness().size(); ++i) { + double a = Traits::direction(i) * node.element.fitness()[i]; + double b = Traits::direction(i) * boxIt->element.fitness()[i]; + + d1 += pow( (a - node.discretized[i]) / eps[i], 2.0); + d2 += pow( (b - node.discretized[i]) / eps[i], 2.0); + } + + if (d1 > d2) { + *boxIt = node; + } + break; + } + } + } + } + + void appendTo(eoPop& pop) const { + for (typename archive_t::const_iterator it = archive.begin(); it != archive.end(); ++it) { + pop.push_back( it->element ); + } + } + + const EOT& selectRandom() const { + int i = rng.random(archive.size()); + typename archive_t::const_iterator it = archive.begin(); + while (i-- > 0) it++; + return it->element; + } + +}; + + +#endif + diff --git a/eo/src/moo/eoFrontSorter.h b/eo/src/moo/eoFrontSorter.h index 9d9feb8e..f4c7c693 100644 --- a/eo/src/moo/eoFrontSorter.h +++ b/eo/src/moo/eoFrontSorter.h @@ -67,12 +67,7 @@ class eoFrontSorter : public eoUF< const eoPop&, const std::vector< std::ve { fitness.resize(_pop.size()); for (unsigned i = 0; i < _pop.size(); ++i) { - std::vector f; - - for (unsigned j = 0; j < Traits::nObjectives(); ++j) { - if (Traits::maximizing(j) != 0) f.push_back( Traits::maximizing(j) * _pop[i].fitness()[j]); - } - fitness[i] = detail::FitnessInfo(f, i); + fitness[i] = detail::FitnessInfo(_pop[i].fitness().normalized(), i); } detail::front_sorter_impl(fitness, fronts, Traits::tol()); @@ -84,12 +79,7 @@ class eoFrontSorter : public eoUF< const eoPop&, const std::vector< std::ve { fitness.resize(_pop.size()); for (unsigned i = 0; i < _pop.size(); ++i) { - std::vector f; - - for (unsigned j = 0; j < Traits::nObjectives(); ++j) { - if (Traits::maximizing(j) != 0) f.push_back( Traits::maximizing(j) * _pop[i]->fitness()[j]); - } - fitness[i] = detail::FitnessInfo(f, i); + fitness[i] = detail::FitnessInfo(_pop[i]->fitness().normalized(), i); } detail::front_sorter_impl(fitness, fronts, Traits::tol()); diff --git a/eo/src/moo/eoMOFitness.h b/eo/src/moo/eoMOFitness.h index 659a32ee..c7799763 100644 --- a/eo/src/moo/eoMOFitness.h +++ b/eo/src/moo/eoMOFitness.h @@ -30,11 +30,11 @@ #include #include #include - +#include /** * eoMOFitnessTraits: a traits class to specify - * the number of objectives and which one are maximizing or not + * the number of objectives and which one are direction or not * See test/t-eoParetoFitness for its use. * * If you define your own, make sure you make the functions static! @@ -43,11 +43,77 @@ class eoMOFitnessTraits { public : + /// Number of Objectives static unsigned nObjectives() { return 2; } - static double maximizing(int which) { return 1; } // by default: all are maximizing, zero will lead to ignored fitness, negative minimizes - static double tol() { return 1e-6; } // tolerance for distance calculations + + /// by default: all are maximizing, zero will lead to ignored fitness, negative to minimization for that objective + static double direction(unsigned which) { return 1; } + + /// tolerance for dominance check (if within this tolerance objective values are considered equal) + static double tol() { return 1e-6; } + }; +namespace dominance { + + template + inline double worst_possible_value(unsigned objective) { + double dir = Traits::direction(objective); + if (dir == 0.) return 0.0; + if (dir < 0.) return std::numeric_limits::infinity(); + return -std::numeric_limits::infinity(); + } + + template + inline double best_possible_value(unsigned objective) { + return -worst_possible_value(objective); + } + + enum dominance_result { non_dominated_equal, first, second, non_dominated }; + + template + inline dominance_result check(const std::vector& p1, const std::vector& p2, double tolerance) { + + bool all_equal = true; + bool a_better_in_one = false; + bool b_better_in_one = false; + + for (unsigned i = 0; i < p1.size(); ++i) { + + double maxim = FitnessDirectionTraits::direction(i); + double aval = maxim * p1[i]; + double bval = maxim * p2[i]; + + if ( fabs(aval-bval) > tolerance ) { + all_equal = false; + if (aval > bval) { + a_better_in_one = true; + + } else { + b_better_in_one = true; + } + // check if we need to go on + + if (a_better_in_one && b_better_in_one) return non_dominated; + } + + } + + if (all_equal) return non_dominated_equal; + + if (a_better_in_one) return first; + // else b dominates a (check for non-dominance done above + return second; + } + + template + inline dominance_result check(const std::vector& p1, const std::vector& p2) { + return check(p1, p2, FitnessTraits::tol()); + } + + +} // namespace dominance + /** eoMOFitness class: std::vector of doubles with overloaded comparison operators. Comparison is done on 'worth'. This worth needs to be set elsewhere. The template argument FitnessTraits defaults to eMOFitnessTraits, which @@ -66,16 +132,22 @@ public : typedef FitnessTraits fitness_traits; - eoMOFitness(double def = 0.0) : std::vector(FitnessTraits::nObjectives(),def), worthValid(false) {} + eoMOFitness() : std::vector(FitnessTraits::nObjectives()), worthValid(false) { reset(); } + + eoMOFitness(double def) : std::vector(FitnessTraits::nObjectives(),def), worthValid(false) {} // Ctr from a std::vector eoMOFitness(std::vector & _v) : std::vector(_v), worthValid(false) {} - /** access to the traits characteristics (so you don't have to write - * a lot of typedef's around - */ - static void setUp(unsigned _n, std::vector & _b) {FitnessTraits::setUp(_n, _b);} - static bool maximizing(unsigned _i) { return FitnessTraits::maximizing(_i);} + void reset() { + + for (unsigned i = 0; i < size(); ++i) { + this->operator[](i) = dominance::worst_possible_value(i); + } + } + + // Make the objective 'feasible' by setting it to the best possible value + void setFeasible(unsigned objective) { this->operator[](objective) = dominance::best_possible_value(objective); } void setWorth(double worth_) { worth = worth_; @@ -91,35 +163,30 @@ public : bool validWorth() const { return worthValid; } void invalidateWorth() { worthValid = false; } + + /// Check on dominance: returns 0 if non-dominated, 1 if this dominates other, -1 if other dominates this + int check_dominance(const eoMOFitness& _other) const + { + dominance::dominance_result dom = dominance::check(*this, _other); - /// Partial order based on Pareto-dominance - //bool operator<(const eoMOFitness& _other) const + return dom == dominance::first? 1 : (dom == dominance::second? -1 : 0); + } + + /// normalized fitness: all maximizing, removed the irrelevant ones + std::vector normalized() const { + std::vector f; + + for (unsigned j = 0; j < FitnessTraits::nObjectives(); ++j) { + if (FitnessTraits::direction(j) != 0) f.push_back( FitnessTraits::direction(j) * this->operator[](j)); + } + + return f; + } + + /// returns true if this dominates other bool dominates(const eoMOFitness& _other) const { - bool dom = false; - - const std::vector& performance = *this; - const std::vector& otherperformance = _other; - - for (unsigned i = 0; i < FitnessTraits::nObjectives(); ++i) - { - double maxim = FitnessTraits::maximizing(i); - double aval = maxim * performance[i]; - double bval = maxim * otherperformance[i]; - - if (fabs(aval-bval)>FitnessTraits::tol) - { - if (aval < bval) - { - return false; // cannot dominate - } - // else aval > bval - dom = true; // for the moment: goto next objective - } - //else they're equal in this objective, goto next - } - - return dom; + return check_dominance(_other) == 1; } /// compare *not* on dominance, but on worth diff --git a/eo/src/moo/eoNSGA_IIa_Eval.cpp b/eo/src/moo/eoNSGA_IIa_Eval.cpp index c7cafef9..33b6d9ba 100644 --- a/eo/src/moo/eoNSGA_IIa_Eval.cpp +++ b/eo/src/moo/eoNSGA_IIa_Eval.cpp @@ -3,7 +3,7 @@ namespace nsga2a { -double distance(const std::vector& f1, const std::vector& f2, const std::vector& range) { +double calc_distance(const std::vector& f1, const std::vector& f2) { double dist = 0; for (unsigned i = 0; i < f1.size(); ++i) { double d = (f1[i] - f2[i]); @@ -40,20 +40,6 @@ for (unsigned i = 0; i < processed.size(); ++i) { } rank--; -// calculate ranges -std::vector mins(nDim, std::numeric_limits::infinity()); -for (unsigned dim = 0; dim < boundaries.size(); ++dim) { - for (unsigned i = 0; i < boundaries.size(); ++i) { - mins[dim] = std::min( mins[dim], boundaries[i].fitness[dim] ); - } -} - -std::vector range(nDim); -for (unsigned dim = 0; dim < nDim; ++dim) { - range[dim] = boundaries[dim].fitness[dim] - mins[dim]; -} - - // clean up processed (make unique) sort(processed.begin(), processed.end()); // swap out last first for (unsigned i = 1; i < processed.size(); ++i) { @@ -81,7 +67,7 @@ unsigned selected = 0; for (unsigned i = 0; i < front.size(); ++i) { for (unsigned k = 0; k < boundaries.size(); ++k) { - double d = distance( front[i].fitness, boundaries[k].fitness, range ); + double d = calc_distance( front[i].fitness, boundaries[k].fitness); if (d < distances[i]) { distances[i] = d; } @@ -111,7 +97,7 @@ while (!front.empty()) { selected = 0; for (unsigned i = 0; i < front.size(); ++i) { - double d = distance(front[i].fitness, last.fitness, range); + double d = calc_distance(front[i].fitness, last.fitness); if (d < distances[i]) { distances[i] = d;