added epsilon moea

eo/src/moo/eoEpsMOEA.h

@@ -0,0 +1,118 @@
+#ifndef eoEpsMOEA_h
+#define eoEpsMOEA_h
+
+#include <eoAlgo.h>
+#include <moo/eoEpsilonArchive.h>
+
+template <class EOT>
+class eoEpsMOEA : public eoAlgo<EOT> {
+ 
+    public:
+  
+    eoEpsMOEA(
+         eoContinue<EOT>& _continuator,
+         eoEvalFunc<EOT>& _eval,
+         eoGenOp<EOT>& _op,
+         const std::vector<double>& eps
+     ) : continuator(_continuator),
+	 eval (_eval),
+	 loopEval(_eval),
+	 popEval(loopEval),
+         op(_op),
+         archive(eps)
+         {}
+
+   
+    void operator()(eoPop<EOT>& pop) {
+        
+        eoPop<EOT> 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<EOT>& 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<EOT> {
+        
+        using eoPopulator< EOT >::src;
+
+        eoEpsilonArchive<EOT>& archive;
+        bool fromArchive;
+        public:
+        epsPopulator(eoEpsilonArchive<EOT>& arch, const eoPop<EOT>& pop, eoPop<EOT>& res) : eoPopulator<EOT>(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<EOT>&          continuator;
+
+  eoEvalFunc <EOT> &        eval ;
+  eoPopLoopEval<EOT>        loopEval;
+
+  eoPopEvalFunc<EOT>&       popEval;
+  
+  eoGenOp<EOT>& op;
+
+  eoEpsilonArchive<EOT> archive;
+
+};
+
+#endif
+

eo/src/moo/eoEpsilonArchive.h

@@ -0,0 +1,123 @@
+#ifndef eoEpsilonArchive_h
+#define eoEpsilonArchive_h
+
+#include <moo/eoMOFitness.h>
+#include <list>
+
+template <class EOT> 
+class eoEpsilonArchive {
+    
+    typedef typename EOT::Fitness::fitness_traits Traits;
+    
+    struct Node {
+        EOT element;
+        std::vector<double> discretized;
+
+        Node(const EOT& eo) : element(eo), discretized(Traits::nObjectives()) {}
+    };
+
+    typedef std::list<Node> archive_t;
+
+    archive_t archive;
+    std::vector<double> eps;
+    
+    public:
+    
+    static double direction(unsigned i) { return 1; }
+    static double tol() { return 1e-6; }
+
+            
+    eoEpsilonArchive(const std::vector<double>& 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<eoEpsilonArchive<EOT> >(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<EOT>& 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
+

eo/src/moo/eoFrontSorter.h

@@ -67,12 +67,7 @@ class eoFrontSorter : public eoUF< const eoPop<EOT>&, const std::vector< std::ve
     {
         fitness.resize(_pop.size());
         for (unsigned i = 0; i < _pop.size(); ++i) {
-            std::vector<double> f;
+            fitness[i] = detail::FitnessInfo(_pop[i].fitness().normalized(), i);
-            
-            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);
         }
 
         detail::front_sorter_impl(fitness, fronts, Traits::tol());
@@ -84,12 +79,7 @@ class eoFrontSorter : public eoUF< const eoPop<EOT>&, const std::vector< std::ve
     {
         fitness.resize(_pop.size());
         for (unsigned i = 0; i < _pop.size(); ++i) {
-            std::vector<double> f;
+            fitness[i] = detail::FitnessInfo(_pop[i]->fitness().normalized(), i);
-            
-            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);
         }
 
         detail::front_sorter_impl(fitness, fronts, Traits::tol());

eo/src/moo/eoMOFitness.h

@@ -30,11 +30,11 @@
 #include <vector>
 #include <stdexcept>
 #include <iostream>
-
+#include <limits>
 
 /**
  * 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 <class Traits>
+    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<double>::infinity();
+        return -std::numeric_limits<double>::infinity();
+    }
+
+    template <class Traits>
+    inline double best_possible_value(unsigned objective) {
+        return -worst_possible_value<Traits>(objective);
+    }
+
+    enum dominance_result { non_dominated_equal, first, second, non_dominated };
+
+    template <class FitnessDirectionTraits>
+    inline dominance_result check(const std::vector<double>& p1, const std::vector<double>& 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 <class FitnessTraits>
+    inline dominance_result check(const std::vector<double>& p1, const std::vector<double>& p2) {
+        return check<FitnessTraits>(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<double>(FitnessTraits::nObjectives(),def), worthValid(false) {}
+  eoMOFitness() : std::vector<double>(FitnessTraits::nObjectives()), worthValid(false) { reset(); }
+  
+  eoMOFitness(double def) : std::vector<double>(FitnessTraits::nObjectives(),def), worthValid(false) {}
 
   // Ctr from a std::vector<double>
   eoMOFitness(std::vector<double> & _v) : std::vector<double>(_v), worthValid(false) {}
 
-  /** access to the traits characteristics (so you don't have to write 
+  void reset() {
-   * a lot of typedef's around
+
-   */
+    for (unsigned i = 0; i < size(); ++i) {
-  static void setUp(unsigned _n, std::vector<bool> & _b) {FitnessTraits::setUp(_n, _b);}
+        this->operator[](i) = dominance::worst_possible_value<FitnessTraits>(i);
-  static bool maximizing(unsigned _i) { return FitnessTraits::maximizing(_i);}
+    }
+  }
+
+  // Make the objective 'feasible' by setting it to the best possible value
+  void setFeasible(unsigned objective) { this->operator[](objective) = dominance::best_possible_value<FitnessTraits>(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<FitnessTraits>& _other) const
+  {
+    dominance::dominance_result dom = dominance::check<FitnessTraits>(*this, _other);
 
-  /// Partial order based on Pareto-dominance
+    return dom == dominance::first? 1 : (dom == dominance::second? -1 : 0);
-  //bool operator<(const eoMOFitness<FitnessTraits>& _other) const
+  }
+
+  /// normalized fitness: all maximizing, removed the irrelevant ones
+  std::vector<double> normalized() const {
+        std::vector<double> 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<FitnessTraits>& _other) const
   {
-    bool dom = false;
+    return check_dominance(_other) == 1;
-
-    const std::vector<double>& performance = *this;
-    const std::vector<double>& 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;
   }
 
   /// compare *not* on dominance, but on worth 

eo/src/moo/eoNSGA_IIa_Eval.cpp

@@ -3,7 +3,7 @@
 
 namespace nsga2a {
 
-double distance(const std::vector<double>& f1, const std::vector<double>& f2, const std::vector<double>& range) {
+double calc_distance(const std::vector<double>& f1, const std::vector<double>& 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<double> mins(nDim, std::numeric_limits<double>::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<double> 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;