source: branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/SimulatedBinaryCrossover.cs

using System;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;

namespace HeuristicLab.Algorithms.NSGA3
{
    /*
     * Implements the simulated binary crossover based on the Tsung Che Chiang's NSGA3 implementation.
     * http://web.ntnu.edu.tw/~tcchiang/publications/nsga3cpp/nsga3cpp.htm
     *
     * The simulated binary crossover in HeuristicLab.Encodings.RealVectorEncoding/Crossovers/SimulatedBinaryCrossover.cs
     * is not usable, because it takes different parameters than stated in the NSGA-III paper. It takes contiguity as a
     * parameter instead of an eta value.
     */

    public static class SimulatedBinaryCrossover
    {
        private const double EPSILON = 10e-6; // a tiny number that is greater than 0

        public static Tuple<RealVector, RealVector> Apply(IRandom random, DoubleMatrix bounds, RealVector parent1, RealVector parent2, double crossoverProbability = 1.0, double eta = 20)
        {
            var child1 = new RealVector(parent1);
            var child2 = new RealVector(parent2);

            if (random.NextDouble() > crossoverProbability) return null;

            for (int i = 0; i < child1.Length; i++)
            {
                if (random.NextDouble() > 0.5) continue; // these two variables are not crossovered
                if (Math.Abs(parent1[i] - parent2[i]) < EPSILON) continue; // two values are the same

                double y1 = Math.Min(parent1[i], parent2[i]);
                double y2 = Math.Max(parent1[i], parent2[i]);

                double lb;
                double ub;
                if (bounds.Rows == 1)
                {
                    lb = bounds[0, 0];
                    ub = bounds[0, 1];
                }
                else
                {
                    lb = bounds[i, 0];
                    ub = bounds[i, 1];
                }

                double rand = random.NextDouble();

                // child 1
                double beta = 1.0 + (2.0 * (y1 - lb) / (y2 - y1));
                double alpha = 2.0 - Math.Pow(beta, -(eta + 1.0));
                double betaq = Get_Betaq(rand, alpha, eta);

                child1[i] = 0.5 * ((y1 + y2) - betaq * (y2 - y1));

                // child 2
                beta = 1.0 + (2.0 * (ub - y2) / (y2 - y1));
                alpha = 2.0 - Math.Pow(beta, -(eta + 1.0));
                betaq = Get_Betaq(rand, alpha, eta);

                child2[i] = 0.5 * ((y1 + y2) + betaq * (y2 - y1));

                // boundary checking
                child1[i] = Math.Min(ub, Math.Max(lb, child1[i]));
                child2[i] = Math.Min(ub, Math.Max(lb, child2[i]));

                if (random.NextDouble() <= 0.5)
                {
                    // swap values
                    var tmp = child1[i];
                    child1[i] = child2[i];
                    child2[i] = tmp;
                }
            }

            return Tuple.Create(child1, child2);
        }

        private static double Get_Betaq(double rand, double alpha, double eta)
        {
            double betaq;
            if (rand <= (1.0 / alpha))
                betaq = Math.Pow((rand * alpha), (1.0 / (eta + 1.0)));
            else
                betaq = Math.Pow((1.0 / (2.0 - rand * alpha)), (1.0 / (eta + 1.0)));
            return betaq;
        }
    }
}