Index: /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/HeuristicLab.Algorithms.NSGA3-3.3.csproj
===================================================================
--- /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/HeuristicLab.Algorithms.NSGA3-3.3.csproj	(revision 17664)
+++ /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/HeuristicLab.Algorithms.NSGA3-3.3.csproj	(revision 17665)
@@ -110,4 +110,5 @@
   <ItemGroup>
     <Compile Include="NSGA3.cs" />
+    <Compile Include="NSGA3Selection.cs" />
     <Compile Include="Plugin.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
Index: /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3.cs
===================================================================
--- /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3.cs	(revision 17664)
+++ /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3.cs	(revision 17665)
@@ -26,6 +26,4 @@
     public class NSGA3 : BasicAlgorithm
     {
-        private const double EPSILON = 10e-6; // a tiny number that is greater than 0
-
         public override bool SupportsPause => false;
 
@@ -53,5 +51,5 @@
 
         [Storable]
-        private List<Solution> solutions;
+        private List<Solution> solutions; // maybe todo: rename to nextGeneration (see Run method)
 
         #endregion Storable fields
@@ -266,5 +264,11 @@
                 // create copies of generated reference points (to preserve the original ones for
                 // the next generation) maybe todo: use cloner?
-                ToNextGeneration(CreateCopyOfReferencePoints());
+                //ToNextGeneration(CreateCopyOfReferencePoints());
+
+                List<Solution> qt = Mutate(Recombine(solutions));
+                List<Solution> rt = Utility.Concat(solutions, qt);
+
+                solutions = NSGA3Selection.SelectSolutionsForNextGeneration(rt, GetCopyOfReferencePoints(), Problem.Maximization, random);
+
                 ResultsCurrentGeneration.Value++;
             }
@@ -275,5 +279,5 @@
         #region Private Methods
 
-        private List<ReferencePoint> CreateCopyOfReferencePoints()
+        private List<ReferencePoint> GetCopyOfReferencePoints()
         {
             if (ReferencePoints == null) return null;
@@ -306,289 +310,4 @@
         {
             return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random);
-        }
-
-        private void ToNextGeneration(List<ReferencePoint> referencePoints)
-        {
-            List<Solution> st = new List<Solution>();
-            List<Solution> qt = Mutate(Recombine(solutions));
-            List<Solution> rt = Utility.Concat(solutions, qt);
-            List<Solution> nextGeneration;
-
-            // Do non-dominated sort
-            var qualities = Utility.ToFitnessMatrix(rt);
-            // compute the pareto fronts using the DominationCalculator and discard the qualities
-            // part in the inner tuples
-            Fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, Problem.Maximization, out int[] rank, true)
-                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
-
-            int i = 0;
-            List<Solution> lf = null; // last front to be included
-            while (i < Fronts.Count && st.Count < PopulationSize.Value)
-            {
-                lf = Fronts[i];
-                st = Utility.Concat(st, lf);
-                i++;
-            }
-
-            if (st.Count == PopulationSize.Value) // no selection needs to be done
-                nextGeneration = st;
-            else
-            {
-                int l = i - 1;
-                nextGeneration = new List<Solution>();
-                for (int f = 0; f < l; f++)
-                    nextGeneration = Utility.Concat(nextGeneration, Fronts[f]);
-                int k = PopulationSize.Value - nextGeneration.Count;
-                Normalize(st);
-                Associate(referencePoints);
-                List<Solution> solutionsToAdd = Niching(k, referencePoints);
-                nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
-            }
-        }
-
-        private void Normalize(List<Solution> population)
-        {
-            // Find the ideal point
-            double[] idealPoint = new double[NumberOfObjectives];
-            for (int j = 0; j < NumberOfObjectives; j++)
-            {
-                // Compute ideal point
-                idealPoint[j] = Utility.Min(s => s.Fitness[j], population);
-
-                // Translate objectives
-                foreach (var solution in population)
-                    solution.Fitness[j] -= idealPoint[j];
-            }
-
-            // Find the extreme points
-            Solution[] extremePoints = new Solution[NumberOfObjectives];
-            for (int j = 0; j < NumberOfObjectives; j++)
-            {
-                // Compute extreme points
-                double[] weights = new double[NumberOfObjectives];
-                for (int i = 0; i < NumberOfObjectives; i++) weights[i] = EPSILON;
-                weights[j] = 1;
-                double func(Solution s) => ASF(s.Fitness, weights);
-                extremePoints[j] = Utility.ArgMin(func, population);
-            }
-
-            // Compute intercepts
-            List<double> intercepts = GetIntercepts(extremePoints.ToList());
-
-            // Normalize objectives
-            NormalizeObjectives(intercepts, idealPoint);
-        }
-
-        private void NormalizeObjectives(List<double> intercepts, double[] idealPoint)
-        {
-            for (int f = 0; f < Fronts.Count; f++)
-            {
-                foreach (var solution in Fronts[f])
-                {
-                    for (int i = 0; i < NumberOfObjectives; i++)
-                    {
-                        if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
-                        {
-                            solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
-                        }
-                        else
-                        {
-                            solution.Fitness[i] = solution.Fitness[i] / EPSILON;
-                        }
-                    }
-                }
-            }
-        }
-
-        private void Associate(List<ReferencePoint> referencePoints)
-        {
-            for (int f = 0; f < Fronts.Count; f++)
-            {
-                foreach (var solution in Fronts[f])
-                {
-                    // find reference point for which the perpendicular distance to the current
-                    // solution is the lowest
-                    var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);
-                    // associated reference point
-                    var arp = rpAndDist.Item1;
-                    // distance to that reference point
-                    var dist = rpAndDist.Item2;
-
-                    if (f + 1 != Fronts.Count)
-                    {
-                        // Todo: Add member for reference point on index min_rp
-                        arp.NumberOfAssociatedSolutions++;
-                    }
-                    else
-                    {
-                        // Todo: Add potential member for reference point on index min_rp
-                        arp.AddPotentialAssociatedSolution(solution, dist);
-                    }
-                }
-            }
-        }
-
-        private List<Solution> Niching(int k, List<ReferencePoint> referencePoints)
-        {
-            List<Solution> solutions = new List<Solution>();
-            while (solutions.Count != k)
-            {
-                ReferencePoint min_rp = FindNicheReferencePoint(referencePoints);
-
-                Solution chosen = SelectClusterMember(min_rp);
-                if (chosen == null)
-                {
-                    referencePoints.Remove(min_rp);
-                }
-                else
-                {
-                    min_rp.NumberOfAssociatedSolutions++;
-                    min_rp.RemovePotentialAssociatedSolution(chosen);
-                    solutions.Add(chosen);
-                }
-            }
-
-            return solutions;
-        }
-
-        private ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints)
-        {
-            // the minimum number of associated solutions for a reference point over all reference points
-            int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);
-
-            // the reference points that share the number of associated solutions where that number
-            // is equal to minNumber
-            List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();
-            foreach (var referencePoint in referencePoints)
-                if (referencePoint.NumberOfAssociatedSolutions == minNumber)
-                    minAssociatedReferencePoints.Add(referencePoint);
-
-            if (minAssociatedReferencePoints.Count > 1)
-                return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];
-            else
-                return minAssociatedReferencePoints.Single();
-        }
-
-        private Solution SelectClusterMember(ReferencePoint referencePoint)
-        {
-            Solution chosen = null;
-            if (referencePoint.HasPotentialMember())
-            {
-                if (referencePoint.NumberOfAssociatedSolutions == 0)
-                    chosen = referencePoint.FindClosestMember();
-                else
-                    chosen = referencePoint.RandomMember();
-            }
-            return chosen;
-        }
-
-        private double GetPerpendicularDistance(double[] values, double[] fitness)
-        {
-            double numerator = 0;
-            double denominator = 0;
-            for (int i = 0; i < values.Length; i++)
-            {
-                numerator += values[i] * fitness[i];
-                denominator += Math.Pow(values[i], 2);
-            }
-            double k = numerator / denominator;
-
-            double d = 0;
-            for (int i = 0; i < values.Length; i++)
-            {
-                d += Math.Pow(k * values[i] - fitness[i], 2);
-            }
-            return Math.Sqrt(d);
-        }
-
-        private double ASF(double[] x, double[] weight)
-        {
-            List<int> dimensions = new List<int>();
-            for (int i = 0; i < NumberOfObjectives; i++) dimensions.Add(i);
-            double f(int dim) => x[dim] / weight[dim];
-            return Utility.Max(f, dimensions);
-        }
-
-        private List<double> GetIntercepts(List<Solution> extremePoints)
-        {
-            // Check whether there are duplicate extreme points. This might happen but the original
-            // paper does not mention how to deal with it.
-            bool duplicate = false;
-            for (int i = 0; !duplicate && i < extremePoints.Count; i++)
-            {
-                for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
-                {
-                    // maybe todo: override Equals method of solution?
-                    duplicate = extremePoints[i].Equals(extremePoints[j]);
-                }
-            }
-
-            List<double> intercepts = new List<double>();
-
-            if (duplicate)
-            { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
-                for (int f = 0; f < NumberOfObjectives; f++)
-                {
-                    // extreme_points[f] stands for the individual with the largest value of
-                    // objective f
-                    intercepts.Add(extremePoints[f].Fitness[f]);
-                }
-            }
-            else
-            {
-                // Find the equation of the hyperplane
-                List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
-                for (int i = 0; i < NumberOfObjectives; i++)
-                {
-                    b.Add(1.0);
-                }
-
-                List<List<double>> a = new List<List<double>>();
-                foreach (Solution s in extremePoints)
-                {
-                    List<double> aux = new List<double>();
-                    for (int i = 0; i < NumberOfObjectives; i++)
-                        aux.Add(s.Fitness[i]);
-                    a.Add(aux);
-                }
-                List<double> x = GaussianElimination(a, b);
-
-                // Find intercepts
-                for (int f = 0; f < NumberOfObjectives; f++)
-                {
-                    intercepts.Add(1.0 / x[f]);
-                }
-            }
-
-            return intercepts;
-        }
-
-        private List<double> GaussianElimination(List<List<double>> a, List<double> b)
-        {
-            List<double> x = new List<double>();
-
-            int n = a.Count;
-            for (int i = 0; i < n; i++)
-                a[i].Add(b[i]);
-
-            for (int @base = 0; @base < n - 1; @base++)
-                for (int target = @base + 1; target < n; target++)
-                {
-                    double ratio = a[target][@base] / a[@base][@base];
-                    for (int term = 0; term < a[@base].Count; term++)
-                        a[target][term] = a[target][term] - a[@base][term] * ratio;
-                }
-
-            for (int i = 0; i < n; i++)
-                x.Add(0.0);
-
-            for (int i = n - 1; i >= 0; i--)
-            {
-                for (int known = i + 1; known < n; known++)
-                    a[i][n] = a[i][n] - a[i][known] * x[known];
-                x[i] = a[i][n] / a[i][i];
-            }
-
-            return x;
         }
 
Index: /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3Selection.cs
===================================================================
--- /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3Selection.cs	(revision 17665)
+++ /branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3Selection.cs	(revision 17665)
@@ -0,0 +1,327 @@
+﻿using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Core;
+using HeuristicLab.Optimization;
+
+namespace HeuristicLab.Algorithms.NSGA3
+{
+    public static class NSGA3Selection
+    {
+        private const double EPSILON = 10e-6; // a tiny number that is greater than 0
+
+        /// <summary>
+        /// TODO: Description If the fitnesses of all solutions in <paramref name="rt" /> do not
+        /// have the same number of objectives, undefined behavior happens.
+        /// </summary>
+        /// <param name="rt"></param>
+        /// <param name="referencePoints"></param>
+        /// <param name="maximization"></param>
+        /// <returns></returns>
+        public static List<Solution> SelectSolutionsForNextGeneration(List<Solution> rt, List<ReferencePoint> referencePoints, bool[] maximization, IRandom random)
+        {
+            int N = rt.Count / 2; // number of solutions in a generation
+            int numberOfObjectives = rt.First().Fitness.Length;
+            List<Solution> st = new List<Solution>();
+
+            // Do non-dominated sort
+            var qualities = Utility.ToFitnessMatrix(rt);
+            // compute the pareto fronts using the DominationCalculator and discard the qualities
+            // part in the inner tuples
+            var fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, maximization, out int[] rank, true)
+                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
+
+            int lf = -1; // last front to be included
+            for (int i = 0; i < fronts.Count && st.Count < N; i++)
+            {
+                st = Utility.Concat(st, fronts[i]);
+                lf = i;
+            }
+
+            List<Solution> nextGeneration = new List<Solution>();
+            if (st.Count == N) // no special selection needs to be done
+                nextGeneration = st;
+            else
+            {
+                // Add every front to the next generation list except for the one that is added only partially
+                nextGeneration = new List<Solution>();
+                for (int f = 0; f < lf - 1; f++)
+                    nextGeneration = Utility.Concat(nextGeneration, fronts[f]);
+                int k = N - nextGeneration.Count;
+                Normalize(st, numberOfObjectives);
+                Associate(fronts, referencePoints);
+                List<Solution> solutionsToAdd = Niching(k, referencePoints, random);
+                nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
+            }
+
+            return nextGeneration;
+        }
+
+        private static void Normalize(List<Solution> st, int numberOfObjectives)
+        {
+            // Find the ideal point
+            double[] idealPoint = new double[numberOfObjectives];
+            for (int j = 0; j < numberOfObjectives; j++)
+            {
+                // Compute ideal point
+                idealPoint[j] = Utility.Min(s => s.Fitness[j], st);
+
+                // Translate objectives
+                foreach (var solution in st)
+                    solution.Fitness[j] -= idealPoint[j];
+            }
+
+            // Find the extreme points
+            Solution[] extremePoints = new Solution[numberOfObjectives];
+            for (int j = 0; j < numberOfObjectives; j++)
+            {
+                // Compute extreme points
+                double[] weights = new double[numberOfObjectives];
+                for (int i = 0; i < numberOfObjectives; i++) weights[i] = EPSILON;
+                weights[j] = 1;
+
+                extremePoints[j] = Utility.ArgMin(s => ASF(s.Fitness, weights), st);
+            }
+
+            // Compute intercepts
+            List<double> intercepts = GetIntercepts(extremePoints.ToList());
+
+            // Normalize objectives
+            NormalizeObjectives(st, intercepts, idealPoint);
+        }
+
+        private static void NormalizeObjectives(List<Solution> st, List<double> intercepts, double[] idealPoint)
+        {
+            int numberOfObjectives = st.First().Fitness.Length;
+
+            foreach (var solution in st)
+            {
+                for (int i = 0; i < numberOfObjectives; i++)
+                {
+                    if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
+                    {
+                        solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
+                    }
+                    else
+                    {
+                        solution.Fitness[i] = solution.Fitness[i] / EPSILON;
+                    }
+                }
+            }
+            // todo: delete this
+            //for (int f = 0; f < Fronts.Count; f++)
+            //{
+            //    foreach (var solution in Fronts[f])
+            //    {
+            //        for (int i = 0; i < NumberOfObjectives; i++)
+            //        {
+            //            if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
+            //            {
+            //                solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
+            //            }
+            //            else
+            //            {
+            //                solution.Fitness[i] = solution.Fitness[i] / EPSILON;
+            //            }
+            //        }
+            //    }
+            //}
+        }
+
+        private static void Associate(List<List<Solution>> fronts, List<ReferencePoint> referencePoints)
+        {
+            for (int f = 0; f < fronts.Count; f++)
+            {
+                foreach (var solution in fronts[f])
+                {
+                    // find reference point for which the perpendicular distance to the current
+                    // solution is the lowest
+                    var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);
+                    // associated reference point
+                    var arp = rpAndDist.Item1;
+                    // distance to that reference point
+                    var dist = rpAndDist.Item2;
+
+                    if (f + 1 != fronts.Count)
+                    {
+                        // Todo: Add member for reference point on index min_rp
+                        arp.NumberOfAssociatedSolutions++;
+                    }
+                    else
+                    {
+                        // Todo: Add potential member for reference point on index min_rp
+                        arp.AddPotentialAssociatedSolution(solution, dist);
+                    }
+                }
+            }
+        }
+
+        private static List<Solution> Niching(int k, List<ReferencePoint> referencePoints, IRandom random)
+        {
+            List<Solution> solutions = new List<Solution>();
+            while (solutions.Count != k)
+            {
+                ReferencePoint min_rp = FindNicheReferencePoint(referencePoints, random);
+
+                Solution chosen = SelectClusterMember(min_rp);
+                if (chosen == null)
+                {
+                    referencePoints.Remove(min_rp);
+                }
+                else
+                {
+                    min_rp.NumberOfAssociatedSolutions++;
+                    min_rp.RemovePotentialAssociatedSolution(chosen);
+                    solutions.Add(chosen);
+                }
+            }
+
+            return solutions;
+        }
+
+        private static ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints, IRandom random)
+        {
+            // the minimum number of associated solutions for a reference point over all reference points
+            int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);
+
+            // the reference points that share the number of associated solutions where that number
+            // is equal to minNumber
+            List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();
+            foreach (var referencePoint in referencePoints)
+                if (referencePoint.NumberOfAssociatedSolutions == minNumber)
+                    minAssociatedReferencePoints.Add(referencePoint);
+
+            if (minAssociatedReferencePoints.Count > 1)
+                return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];
+            else
+                return minAssociatedReferencePoints.Single();
+        }
+
+        private static Solution SelectClusterMember(ReferencePoint referencePoint)
+        {
+            Solution chosen = null;
+            if (referencePoint.HasPotentialMember())
+            {
+                if (referencePoint.NumberOfAssociatedSolutions == 0)
+                    chosen = referencePoint.FindClosestMember();
+                else
+                    chosen = referencePoint.RandomMember();
+            }
+            return chosen;
+        }
+
+        private static double GetPerpendicularDistance(double[] values, double[] fitness)
+        {
+            double numerator = 0;
+            double denominator = 0;
+            for (int i = 0; i < values.Length; i++)
+            {
+                numerator += values[i] * fitness[i];
+                denominator += Math.Pow(values[i], 2);
+            }
+            double k = numerator / denominator;
+
+            double d = 0;
+            for (int i = 0; i < values.Length; i++)
+            {
+                d += Math.Pow(k * values[i] - fitness[i], 2);
+            }
+            return Math.Sqrt(d);
+        }
+
+        private static double ASF(double[] x, double[] weight)
+        {
+            int numberOfObjectives = x.Length;
+            List<int> dimensions = new List<int>();
+            for (int i = 0; i < numberOfObjectives; i++)
+                dimensions.Add(i);
+
+            return Utility.Max((dim) => x[dim] / weight[dim], dimensions);
+        }
+
+        private static List<double> GetIntercepts(List<Solution> extremePoints)
+        {
+            int numberOfObjectives = extremePoints.First().Fitness.Length;
+
+            // Check whether there are duplicate extreme points. This might happen but the original
+            // paper does not mention how to deal with it.
+            bool duplicate = false;
+            for (int i = 0; !duplicate && i < extremePoints.Count; i++)
+            {
+                for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
+                {
+                    // maybe todo: override Equals method of solution?
+                    duplicate = extremePoints[i].Equals(extremePoints[j]);
+                }
+            }
+
+            List<double> intercepts = new List<double>();
+
+            if (duplicate)
+            { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
+                for (int f = 0; f < numberOfObjectives; f++)
+                {
+                    // extreme_points[f] stands for the individual with the largest value of
+                    // objective f
+                    intercepts.Add(extremePoints[f].Fitness[f]);
+                }
+            }
+            else
+            {
+                // Find the equation of the hyperplane
+                List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
+                for (int i = 0; i < numberOfObjectives; i++)
+                {
+                    b.Add(1.0);
+                }
+
+                List<List<double>> a = new List<List<double>>();
+                foreach (Solution s in extremePoints)
+                {
+                    List<double> aux = new List<double>();
+                    for (int i = 0; i < numberOfObjectives; i++)
+                        aux.Add(s.Fitness[i]);
+                    a.Add(aux);
+                }
+                List<double> x = GaussianElimination(a, b);
+
+                // Find intercepts
+                for (int f = 0; f < numberOfObjectives; f++)
+                {
+                    intercepts.Add(1.0 / x[f]);
+                }
+            }
+
+            return intercepts;
+        }
+
+        private static List<double> GaussianElimination(List<List<double>> a, List<double> b)
+        {
+            List<double> x = new List<double>();
+
+            int n = a.Count;
+            for (int i = 0; i < n; i++)
+                a[i].Add(b[i]);
+
+            for (int @base = 0; @base < n - 1; @base++)
+                for (int target = @base + 1; target < n; target++)
+                {
+                    double ratio = a[target][@base] / a[@base][@base];
+                    for (int term = 0; term < a[@base].Count; term++)
+                        a[target][term] = a[target][term] - a[@base][term] * ratio;
+                }
+
+            for (int i = 0; i < n; i++)
+                x.Add(0.0);
+
+            for (int i = n - 1; i >= 0; i--)
+            {
+                for (int known = i + 1; known < n; known++)
+                    a[i][n] = a[i][n] - a[i][known] * x[known];
+                x[i] = a[i][n] / a[i][i];
+            }
+
+            return x;
+        }
+    }
+}
