diff --git a/eo/src/moo/eoEpsMOEA.h b/eo/src/moo/eoEpsMOEA.h
index f847a717..9ad8157a 100644
--- a/eo/src/moo/eoEpsMOEA.h
+++ b/eo/src/moo/eoEpsMOEA.h
@@ -3,6 +3,8 @@
 
 #include <eoAlgo.h>
 #include <moo/eoEpsilonArchive.h>
+#include <utils/eoStat.h>
+
 
 template <class EOT>
 class eoEpsMOEA : public eoAlgo<EOT> {
@@ -13,14 +15,16 @@ class eoEpsMOEA : public eoAlgo<EOT> {
          eoContinue<EOT>& _continuator,
          eoEvalFunc<EOT>& _eval,
          eoGenOp<EOT>& _op,
-         const std::vector<double>& eps
+         const std::vector<double>& eps,
+         unsigned max_archive_size
      ) : continuator(_continuator),
 	 eval (_eval),
 	 loopEval(_eval),
 	 popEval(loopEval),
          op(_op),
-         archive(eps)
-         {}
+         archive(eps, max_archive_size)
+         {
+         }
 
    
     void operator()(eoPop<EOT>& pop) {
@@ -48,7 +52,11 @@ class eoEpsMOEA : public eoAlgo<EOT> {
                 }
             }
 
-        } while (continuator(pop));
+            // quite expensive to copy the entire archive time and time again, but this will make it work more seamlessly with the rest of EO
+            offspring.clear();
+            archive.appendTo(offspring);
+
+        } while (continuator(offspring)); // check archive
         
         // return archive
         pop.clear();
@@ -58,13 +66,23 @@ class eoEpsMOEA : public eoAlgo<EOT> {
     
     private :
     void update_pop(eoPop<EOT>& pop, const EOT& offspring) {
+        std::vector<unsigned> dominated;
+
         for (unsigned i = 0; i < pop.size(); ++i) {
             int dom = offspring.fitness().check_dominance(pop[i].fitness());
-            if (dom == 1) {
-                pop[i] = offspring;
-                return;
+            switch (dom) {
+                case 1 : // indy dominates this 
+                    dominated.push_back(i);
+                    break;
+                case -1 : // is dominated, do not insert
+                    return;
+                case 0: // incomparable
+                    break;
             }
-            if (dom == -1) return; //
+        }
+        
+        if (dominated.size()) {
+            pop[ dominated[ rng.random(dominated.size()) ] ] = offspring;
         }
 
         // non-dominated everywhere, overwrite random one
@@ -111,7 +129,7 @@ class eoEpsMOEA : public eoAlgo<EOT> {
   eoGenOp<EOT>& op;
 
   eoEpsilonArchive<EOT> archive;
-
+    
 };
 
 #endif
diff --git a/eo/src/moo/eoEpsilonArchive.h b/eo/src/moo/eoEpsilonArchive.h
index d20de936..00f001a1 100644
--- a/eo/src/moo/eoEpsilonArchive.h
+++ b/eo/src/moo/eoEpsilonArchive.h
@@ -11,15 +11,21 @@ class eoEpsilonArchive {
     
     struct Node {
         EOT element;
-        std::vector<double> discretized;
+        std::vector<long> discretized;
 
         Node(const EOT& eo) : element(eo), discretized(Traits::nObjectives()) {}
+    
+        dominance::dominance_result check(const Node& other) const {
+            return dominance::check_discrete(discretized, other.discretized);
+        }
     };
 
-    typedef std::list<Node> archive_t;
+    typedef std::vector<Node> archive_t;
 
     archive_t archive;
-    std::vector<double> eps;
+    std::vector<double> inv_eps;
+   
+    unsigned max_size;
     
     public:
     
@@ -27,8 +33,12 @@ class eoEpsilonArchive {
     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");    
+    eoEpsilonArchive(const std::vector<double>& eps_, unsigned max_size_ = std::numeric_limits<unsigned>::max()) : max_size(max_size_) {
+        inv_eps.resize(eps_.size());
+        for (unsigned i = 0; i < inv_eps.size(); ++i) {
+            inv_eps[i] = 1.0 / eps_[i];
+        }
+        if (inv_eps.size() != Traits::nObjectives()) throw std::logic_error("eoEpsilonArchive: need one epsilon for each objective");    
     }
 
     bool empty() { return archive.size() == 0; }
@@ -43,29 +53,31 @@ class eoEpsilonArchive {
         for (unsigned i = 0; i < eo.fitness().size(); ++i) {
             double val = Traits::direction(i) * eo.fitness()[i];
             
-            node.discretized[i] = floor(val/eps[i]);
+            node.discretized[i] = (long) floor(val*inv_eps[i]);
         }
 
         using namespace dominance;
 
-        typename archive_t::iterator boxIt = archive.end();
+        unsigned box = archive.size();
         // update archive
         
-        for (typename archive_t::iterator it = archive.begin(); it != archive.end(); ++it) {
-            dominance_result result = check<eoEpsilonArchive<EOT> >(node.discretized, it->discretized);
+
+        for (unsigned i = 0; i != archive.size(); ++i) {
+            dominance_result result = node.check(archive[i]); //check<eoEpsilonArchive<EOT> >(node.discretized, archive[i].discretized);
 
             switch (result) {
                 case first : {                          // remove dominated archive member
-                    it = archive.erase(it);
+                    std::swap( archive[i], archive.back());
+                    archive.pop_back();
+                    --i;
                     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;
+                    box = i;
                     goto exitLoop; // break            
                 }
             }
@@ -74,33 +86,60 @@ class eoEpsilonArchive {
 
 exitLoop:
         // insert 
-        if (boxIt == archive.end()) { // non-dominated, new box
+        if (box >= archive.size()) { // non-dominated, new box
             archive.push_back(node);
         } else { // fight it out
-            int dom = node.element.fitness().check_dominance( boxIt->element.fitness() );
+            int dom = node.element.fitness().check_dominance( archive[box].element.fitness() );
             
             switch (dom) {
-                case 1: *boxIt = node; break;
+                case 1: archive[box] = 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);
+                        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] * inv_eps[i];
+                            double b = Traits::direction(i) * archive[box].element.fitness()[i] * inv_eps[i];
+                            
+                            d1 += pow( a - node.discretized[i], 2.0);
+                            d2 += pow( b - node.discretized[i], 2.0);
+                        }
+                        
+                        if (d1 > d2) {
+                            archive[box] = node; 
+                        }
+                        break;
                     }
-                    
-                    if (d1 > d2) {
-                        *boxIt = node; 
-                    }
-                    break;
-                }
             }
         }
+        
+        /*
+        static unsigned counter = 0;
+        if (++counter % 500 == 0) {
+            std::vector<long> mins(archive[0].discretized.size(), std::numeric_limits<long>::max());
+            std::vector<long> maxs(archive[0].discretized.size(), std::numeric_limits<long>::min());
+            for (unsigned i = 0; i < archive.size(); ++i) {
+                for (unsigned dim = 0; dim < archive[i].discretized.size(); ++dim) {
+                    mins[dim] = std::min( mins[dim], archive[i].discretized[dim] );
+                    maxs[dim] = std::max( maxs[dim], archive[i].discretized[dim] );
+                }
+            }
+            
+            std::cout << "Range ";
+            for (unsigned dim = 0; dim < mins.size(); ++dim) {
+                std::cout << (maxs[dim] - mins[dim]) << ' ';
+            }
+            std::cout << archive.size() << std::endl;
+
+        }*/
+
+        if (archive.size() > max_size) {
+            unsigned idx = rng.random(archive.size());
+            if (idx != archive.size()-1) std::swap(archive[idx], archive.back());
+            archive.pop_back();
+        }
+
     }
 
     void appendTo(eoPop<EOT>& pop) const {
@@ -111,9 +150,7 @@ exitLoop:
 
     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;
+        return archive[i].element;
     }
 
 };
diff --git a/eo/src/moo/eoFrontSorter.cpp b/eo/src/moo/eoFrontSorter.cpp
index 3afc762a..d852b559 100644
--- a/eo/src/moo/eoFrontSorter.cpp
+++ b/eo/src/moo/eoFrontSorter.cpp
@@ -4,20 +4,6 @@ using namespace std;
 
 namespace detail {
 
-namespace {
-    struct CompareOn {
-        unsigned dim;
-        double tol;
-
-        CompareOn(unsigned d, double t) : dim(d), tol(t) {}
-
-        bool operator()(const FitnessInfo& a, const FitnessInfo& b) {
-            return  a.fitness[dim] > b.fitness[dim] && fabs(a.fitness[dim] - b.fitness[dim]) > tol;
-        }
-
-    };
-} // end anonymous namespace
-
 void one_objective(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& front, double tol)
 {
     std::sort(fitness.begin(), fitness.end(), CompareOn(0, tol));
@@ -136,6 +122,7 @@ void m_objectives(std::vector<FitnessInfo>& fitness, std::vector< std::vector<Fi
 }
 
 void front_sorter_impl(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& front_indices, double tol) {
+        
         switch (fitness[0].fitness.size())
 	{
 	    case 1:
diff --git a/eo/src/moo/eoFrontSorter.h b/eo/src/moo/eoFrontSorter.h
index f4c7c693..077da880 100644
--- a/eo/src/moo/eoFrontSorter.h
+++ b/eo/src/moo/eoFrontSorter.h
@@ -43,6 +43,18 @@ namespace detail {
           FitnessInfo(const std::vector<double>& fitness_, unsigned index_) : fitness(fitness_), index(index_) {}
     };
     
+    struct CompareOn {
+        unsigned dim;
+        double tol;
+
+        CompareOn(unsigned d, double t) : dim(d), tol(t) {}
+
+        bool operator()(const FitnessInfo& a, const FitnessInfo& b) {
+            return  a.fitness[dim] > b.fitness[dim] && fabs(a.fitness[dim] - b.fitness[dim]) > tol;
+        }
+
+    };
+    
     extern void front_sorter_impl(std::vector<FitnessInfo>& fitness, std::vector< std::vector<FitnessInfo> >& fronts, double tol);
 
 } // namespace detail
diff --git a/eo/src/moo/eoMOFitness.h b/eo/src/moo/eoMOFitness.h
index c7799763..e476e3f4 100644
--- a/eo/src/moo/eoMOFitness.h
+++ b/eo/src/moo/eoMOFitness.h
@@ -105,12 +105,44 @@ namespace dominance {
         // else b dominates a (check for non-dominance done above
         return second;
     }
+
+    template <class IntType>
+    inline dominance_result check_discrete(const std::vector<IntType>& a, const std::vector<IntType>& b) {
+
+        bool all_equal = true;
+        bool a_better_in_one = false;
+        bool b_better_in_one = false;
+
+        for (unsigned i = 0; i < a.size(); ++i) {
+        
+            if ( a[i] != b[i] ) {
+                all_equal = false;
+                if (a[i] > b[i]) {
+                    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
 
@@ -134,10 +166,22 @@ public :
     
   eoMOFitness() : std::vector<double>(FitnessTraits::nObjectives()), worthValid(false) { reset(); }
   
-  eoMOFitness(double def) : std::vector<double>(FitnessTraits::nObjectives(),def), worthValid(false) {}
+  explicit 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) {}
+  explicit eoMOFitness(std::vector<double> & _v) : std::vector<double>(_v), worthValid(false) {}
+
+  virtual ~eoMOFitness() {} // in case someone wants to subclass
+  eoMOFitness(const eoMOFitness<FitnessTraits>& org) { operator=(org); }
+
+  eoMOFitness<FitnessTraits>& operator=(const eoMOFitness<FitnessTraits>& org) {
+        
+        std::vector<double>::operator=(org);
+        worth = org.worth;
+        worthValid = org.worthValid;
+
+        return *this;
+  }
 
   void reset() {
 
@@ -177,7 +221,9 @@ public :
         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));
+            if (FitnessTraits::direction(j) != 0) {
+                f.push_back( FitnessTraits::direction(j) * this->operator[](j));
+            }
         }
         
         return f;
@@ -192,17 +238,17 @@ public :
   /// compare *not* on dominance, but on worth 
   bool operator<(const eoMOFitness<FitnessTraits>& _other) const
   {
-    return getWorth() > _other.getWorth();
+    return getWorth() < _other.getWorth();
   }
 
   bool operator>(const eoMOFitness<FitnessTraits>& _other) const
   {
-    return _other < *this;
+    return _other > *this;
   }
 
   bool operator<=(const eoMOFitness<FitnessTraits>& _other) const
   {
-    return getWorth() >= _other.getWorth();
+    return getWorth() <= _other.getWorth();
   }
 
   bool operator>=(const eoMOFitness<FitnessTraits>& _other) const
@@ -242,4 +288,11 @@ public :
 
 };
 
+namespace dominance {
+    template <class FitnessTraits>
+    inline dominance_result check(const eoMOFitness<FitnessTraits>& p1, const eoMOFitness<FitnessTraits>& p2) {
+        return check<FitnessTraits>(p1, p2, FitnessTraits::tol());   
+    }
+}
+
 #endif
diff --git a/eo/src/moo/eoNSGA_II_Eval.cpp b/eo/src/moo/eoNSGA_II_Eval.cpp
index 321bc079..72916a31 100644
--- a/eo/src/moo/eoNSGA_II_Eval.cpp
+++ b/eo/src/moo/eoNSGA_II_Eval.cpp
@@ -1,6 +1,8 @@
 #include <eoNSGA_II_Eval.h>
 
 namespace nsga2 {
+    
+  using namespace std;
 
   typedef std::pair<double, unsigned> double_index_pair;
 
@@ -15,8 +17,6 @@ namespace nsga2 {
 
   void assign_worths(const std::vector<detail::FitnessInfo>& front, unsigned rank, std::vector<double>& worths) {
   
-    unsigned i;
-
     unsigned nObjectives = front[0].fitness.size(); //traits::nObjectives(); //_pop[_cf[0]].fitness().size();
     
     std::vector<double> niche_distance(front.size());
@@ -25,7 +25,7 @@ namespace nsga2 {
     {
 
       std::vector<std::pair<double, unsigned> > performance(front.size());
-      for (i =0; i < front.size(); ++i)
+      for (unsigned i =0; i < front.size(); ++i)
       {
         performance[i].first = front[i].fitness[o];
         performance[i].second = i;
@@ -35,16 +35,18 @@ namespace nsga2 {
 
       std::vector<double> nc(front.size(), 0.0);
 
-      for (i = 1; i < front.size()-1; ++i)
+      for (unsigned i = 1; i < front.size()-1; ++i)
       { // and yet another level of indirection
+        
         nc[performance[i].second] = performance[i+1].first - performance[i-1].first;
+        
       }
 
       // set boundary at max_dist + 1 (so it will get chosen over all the others
       //nc[performance[0].second]     += 0;
-      nc[performance.back().second] += std::numeric_limits<double>::infinity(); // best on objective
+      nc[performance.back().second] = std::numeric_limits<double>::infinity(); // best on objective
       
-      for (i = 0; i < nc.size(); ++i)
+      for (unsigned i = 0; i < nc.size(); ++i)
       {
         niche_distance[i] += nc[i];
       }
@@ -53,12 +55,12 @@ namespace nsga2 {
     // now we've got niche_distances, scale them between (0, 1), making sure that infinities get maximum rank
     
     double max = 0;
-    for (unsigned i = 0; i < niche_distance[i]; ++i) {
+    for (unsigned i = 0; i < niche_distance.size(); ++i) {
         if (niche_distance[i] != std::numeric_limits<double>::infinity()) {
             max = std::max(max, niche_distance[i]);
         }
     }
-    
+     
     for (unsigned i = 0; i < front.size(); ++i) {
         double dist = niche_distance[i];
         if (dist == std::numeric_limits<double>::infinity()) {
diff --git a/eo/src/moo/eoNSGA_II_Eval.h b/eo/src/moo/eoNSGA_II_Eval.h
index de3e33f8..55b3fd60 100644
--- a/eo/src/moo/eoNSGA_II_Eval.h
+++ b/eo/src/moo/eoNSGA_II_Eval.h
@@ -34,15 +34,15 @@ class eoNSGA_II_Eval : public eoMOEval<EOT>
 
         typename eoFrontSorter<EOT>::front_t front = sorter(pop);
         
+        // calculate worths
         std::vector<double> worths(pop.size());
         for (unsigned i = 0; i < front.size(); ++i) {
             nsga2::assign_worths(front[i], front.size() - i, worths); 
         }
-
+         
+        // set worths
         for (unsigned i = 0; i < pop.size(); ++i) {
-            typename EOT::Fitness f = pop[i]->fitness();
-            f.setWorth(worths[i]);
-            pop[i]->fitness(f);
+            pop[i]->fitnessReference().setWorth( worths[i]);
         }
 
     }
diff --git a/eo/src/moo/eoNSGA_IIa_Eval.cpp b/eo/src/moo/eoNSGA_IIa_Eval.cpp
index 33b6d9ba..5218b101 100644
--- a/eo/src/moo/eoNSGA_IIa_Eval.cpp
+++ b/eo/src/moo/eoNSGA_IIa_Eval.cpp
@@ -7,109 +7,221 @@ double calc_distance(const std::vector<double>& f1, const std::vector<double>& f
     double dist = 0;
     for (unsigned i = 0; i < f1.size(); ++i) {
         double d = (f1[i] - f2[i]);
-        dist += d*d;
+        dist += pow(fabs(d), 1./2);
+        //dist += pow(fabs(d), 2.0);
     }
     return dist;
 }
 
-
 unsigned assign_worths(std::vector<detail::FitnessInfo> front, unsigned rank, std::vector<double>& worths) {
 
-unsigned nDim = front[0].fitness.size();
+    unsigned nDim = front[0].fitness.size();
 
-// find boundary points
-std::vector<unsigned> processed(nDim);
+    // find boundary points
+    std::vector<unsigned> processed(nDim);
 
-for (unsigned i = 1; i < front.size(); ++i) {
-    for (unsigned dim = 0; dim < nDim; ++dim) {
-        if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
-            processed[dim] = i;
+    for (unsigned i = 1; i < front.size(); ++i) {
+        for (unsigned dim = 0; dim < nDim; ++dim) {
+            if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
+                processed[dim] = i;
+            }
         }
     }
-}
-   
-// assign fitness to processed and store in boundaries
-std::vector<detail::FitnessInfo> boundaries;
-for (unsigned i = 0; i < processed.size(); ++i) {
-    
-    worths[ front[ processed[i] ].index] = rank;
-    
-    //cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
-
-    boundaries.push_back( front[ processed[i] ] );
-}
-rank--;
-
-// clean up processed (make unique) 
-sort(processed.begin(), processed.end()); // swap out last first
-for (unsigned i = 1; i < processed.size(); ++i) {
-    if (processed[i] == processed[i-1]) {
-        std::swap(processed[i], processed.back());
-        processed.pop_back();
-        --i;
-    }
-}
-// remove boundaries from front
-while (processed.size()) {
-    unsigned idx = processed.back();
-    //std::cout << idx << ' ' ;
-    processed.pop_back();
-    std::swap(front.back(), front[idx]);
-    front.pop_back();
-}
-//std::cout << std::endl;
-
-// calculate distances
-std::vector<double> distances(front.size(), std::numeric_limits<double>::infinity());
-
-unsigned selected = 0;
-// select based on maximum distance to nearest processed point
-for (unsigned i = 0; i < front.size(); ++i) {
-    
-    for (unsigned k = 0; k < boundaries.size(); ++k) {
-        double d = calc_distance( front[i].fitness, boundaries[k].fitness);
-        if (d < distances[i]) {
-            distances[i] = d;
-        }
-    }
-    
-    if (distances[i] > distances[selected]) {
-        selected = i;
-    }
-}
-
- 
-while (!front.empty()) {
-    
-    detail::FitnessInfo last = front[selected];
-    
-    std::swap(front[selected], front.back());
-    front.pop_back();
-
-    std::swap(distances[selected], distances.back());
-    distances.pop_back();
-
-    // set worth
-    worths[last.index] = rank--;
-
-    if (front.empty()) break;
-
-    selected = 0;
-
-    for (unsigned i = 0; i < front.size(); ++i) {
-        double d = calc_distance(front[i].fitness, last.fitness);
+       
+    // assign fitness to processed and store in boundaries
+    std::vector<detail::FitnessInfo> boundaries;
+    for (unsigned i = 0; i < processed.size(); ++i) {
         
-        if (d < distances[i]) {
-            distances[i] = d;
-        }
+        worths[ front[ processed[i] ].index] = rank;
+        
+        //cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
 
-        if (distances[i] >= distances[selected]) {
+        boundaries.push_back( front[ processed[i] ] );
+    }
+    rank--;
+
+    // clean up processed (make unique) 
+    sort(processed.begin(), processed.end()); // swap out last first
+    for (unsigned i = 1; i < processed.size(); ++i) {
+        if (processed[i] == processed[i-1]) {
+            std::swap(processed[i], processed.back());
+            processed.pop_back();
+            --i;
+        }
+    }
+    // remove boundaries from front
+    while (processed.size()) {
+        unsigned idx = processed.back();
+        //std::cout << idx << ' ' ;
+        processed.pop_back();
+        std::swap(front.back(), front[idx]);
+        front.pop_back();
+    }
+    //std::cout << std::endl;
+
+    // calculate distances
+    std::vector<double> distances(front.size(), std::numeric_limits<double>::infinity());
+
+    unsigned selected = 0;
+    // select based on maximum distance to nearest processed point
+    for (unsigned i = 0; i < front.size(); ++i) {
+        
+        for (unsigned k = 0; k < boundaries.size(); ++k) {
+            double d = calc_distance( front[i].fitness, boundaries[k].fitness);
+            if (d < distances[i]) {
+                distances[i] = d;
+            }
+        }
+        
+        if (distances[i] > distances[selected]) {
             selected = i;
         }
     }
+
+     
+    while (!front.empty()) {
+        
+        detail::FitnessInfo last = front[selected];
+        
+        std::swap(front[selected], front.back());
+        front.pop_back();
+
+        std::swap(distances[selected], distances.back());
+        distances.pop_back();
+
+        // set worth
+        worths[last.index] = rank--;
+
+        if (front.empty()) break;
+
+        selected = 0;
+
+        for (unsigned i = 0; i < front.size(); ++i) {
+            double d = calc_distance(front[i].fitness, last.fitness);
+            
+            if (d < distances[i]) {
+                distances[i] = d;
+            }
+
+            if (distances[i] >= distances[selected]) {
+                selected = i;
+            }
+        }
+    }
+
+    return rank;
 }
 
-return rank;
+unsigned assign_worths2(std::vector<detail::FitnessInfo> front, unsigned rank, std::vector<double>& worths) {
+
+    unsigned nDim = front[0].fitness.size();
+
+    // find boundary points
+    std::vector<unsigned> processed(nDim);
+
+    for (unsigned i = 1; i < front.size(); ++i) {
+        for (unsigned dim = 0; dim < nDim; ++dim) {
+            if (front[i].fitness[dim] > front[processed[dim]].fitness[dim]) {
+                processed[dim] = i;
+            }
+        }
+    }
+       
+    // assign fitness to processed and store in boundaries
+    std::vector<detail::FitnessInfo> boundaries;
+    for (unsigned i = 0; i < processed.size(); ++i) {
+        
+        worths[ front[ processed[i] ].index] = rank;
+        
+        //cout << "Boundary " << i << ' ' << front[processed[i]].index << ' ' << parents[ front[ processed[i] ].index]->fitness() << endl;
+
+        boundaries.push_back( front[ processed[i] ] );
+    }
+    rank--;
+
+    // clean up processed (make unique) 
+    sort(processed.begin(), processed.end()); // swap out last first
+    for (unsigned i = 1; i < processed.size(); ++i) {
+        if (processed[i] == processed[i-1]) {
+            std::swap(processed[i], processed.back());
+            processed.pop_back();
+            --i;
+        }
+    }
+    // remove boundaries from front
+    while (processed.size()) {
+        unsigned idx = processed.back();
+        //std::cout << idx << ' ' ;
+        processed.pop_back();
+        std::swap(front.back(), front[idx]);
+        front.pop_back();
+    }
+    //std::cout << std::endl;
+
+    // calculate distances
+    std::vector< std::vector<double> > distances(front.size(), std::vector<double>(nDim, std::numeric_limits<double>::infinity()));
+
+    double bestsum = 0;
+    unsigned selected = 0;
+    // select based on maximum distance to nearest processed point
+    for (unsigned i = 0; i < front.size(); ++i) {
+        
+        for (unsigned k = 0; k < boundaries.size(); ++k) {
+            for (unsigned dim = 0; dim < nDim; ++dim) {
+                double d = front[i].fitness[dim] - boundaries[k].fitness[dim];
+                if (d > 0 && d < distances[i][dim]) {
+                    distances[i][dim] = d;
+                }
+            }
+        }
+
+        double sum = 0;
+        for (unsigned dim = 0; dim < nDim; ++dim) sum += distances[i][dim];
+
+        if (sum > bestsum) {
+            selected = i;
+            bestsum = sum;
+        }
+    }
+
+     
+    while (!front.empty()) {
+        
+        detail::FitnessInfo last = front[selected];
+        
+        std::swap(front[selected], front.back());
+        front.pop_back();
+
+        std::swap(distances[selected], distances.back());
+        distances.pop_back();
+
+        // set worth
+        worths[last.index] = rank--;
+
+        if (front.empty()) break;
+
+        selected = 0;
+        bestsum = 0;
+
+        for (unsigned i = 0; i < front.size(); ++i) {
+            double sum = 0;
+            for (unsigned dim = 0; dim < nDim; ++dim) {
+                double d = front[i].fitness[dim] - last.fitness[dim];
+                if (d > 0 && d < distances[i][dim]) {
+                    distances[i][dim] = d;
+                }
+
+                sum += distances[i][dim];
+            }
+
+            if (sum > bestsum) {
+                selected = i;
+                bestsum = sum;
+            }
+        }
+    }
+
+    return rank;
 }
 
 }