source: branches/ichiriac/HeuristicLab.Algorithms.GDE3/GDE3.cs @ 13749

using System;
using System.Linq;
using System.Collections.Generic;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Operators;
using HeuristicLab.Optimization;
using HeuristicLab.Optimization.Operators;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.PluginInfrastructure;
using HeuristicLab.Problems.MultiObjectiveTestFunctions;
using HeuristicLab.Random;
using System.Threading;

namespace HeuristicLab.Algoritms.GDE3
{

    [Item("Generalized Differential Evolution (GDE3)", "A generalized differential evolution algorithm.")]
    [StorableClass]
    [Creatable(CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 400)]
    public class GDE3 : BasicAlgorithm 
    {
        public Func<IEnumerable<double>, double> Evaluation;
        private enum DominationResult { Dominates, IsDominated, IsNonDominated };
        IComparer<RealVector> dominance;
        public double[] ReferencePoint(int objectives)
        {
            return new double[] { 11, 11 };
        }

        public override Type ProblemType
        {
            get { return typeof(MultiObjectiveTestFunctionProblem); }
        }
         public new MultiObjectiveTestFunctionProblem Problem
        {
            get { return (MultiObjectiveTestFunctionProblem)base.Problem; }
            set { base.Problem = value; }
        }

        public ILookupParameter<DoubleMatrix> BestKnownFrontParameter
        {
            get
            {
                return (ILookupParameter<DoubleMatrix>)Parameters["BestKnownFront"];
            }
        }

        private readonly IRandom _random = new MersenneTwister();
        private int evals;

        #region ParameterNames
        private const string MaximumEvaluationsParameterName = "Maximum Evaluations";
        private const string CrossoverProbabilityParameterName = "CrossoverProbability";
        private const string PopulationSizeParameterName = "PopulationSize";
        private const string ScalingFactorParameterName = "ScalingFactor";
        #endregion

        #region ParameterProperties
        public IFixedValueParameter<IntValue> MaximumEvaluationsParameter
        {
            get { return (IFixedValueParameter<IntValue>)Parameters[MaximumEvaluationsParameterName]; }
        }
        private ValueParameter<IntValue> PopulationSizeParameter
        {
            get { return (ValueParameter<IntValue>)Parameters[PopulationSizeParameterName]; }
        }
        public ValueParameter<DoubleValue> CrossoverProbabilityParameter
        {
            get { return (ValueParameter<DoubleValue>)Parameters[CrossoverProbabilityParameterName]; }
        }
        public ValueParameter<DoubleValue> ScalingFactorParameter
        {
            get { return (ValueParameter<DoubleValue>)Parameters[ScalingFactorParameterName]; }
        }
        #endregion

        #region Properties
        public int MaximumEvaluations
        {
            get { return MaximumEvaluationsParameter.Value.Value; }
            set { MaximumEvaluationsParameter.Value.Value = value; }
        }

        public Double CrossoverProbability
        {
            get { return CrossoverProbabilityParameter.Value.Value; }
            set { CrossoverProbabilityParameter.Value.Value = value; }
        }
        public Double ScalingFactor
        {
            get { return ScalingFactorParameter.Value.Value; }
            set { ScalingFactorParameter.Value.Value = value; }
        }
        public IntValue PopulationSize
        {
            get { return PopulationSizeParameter.Value; }
            set { PopulationSizeParameter.Value = value; }
        }
        #endregion

        #region ResultsProperties
        private double ResultsBestQuality
        {
            get { return ((DoubleValue)Results["Best Quality"].Value).Value; }
            set { ((DoubleValue)Results["Best Quality"].Value).Value = value; }
        }

        private double ResultsGenerationalDistance
        {
            get { return ((DoubleValue)Results["InvertedGenerationalDistance"].Value).Value; }
            set { ((DoubleValue)Results["InvertedGenerationalDistance"].Value).Value = value; }
        }

        private double ResultsHypervolume
        {
            get { return ((DoubleValue)Results["HipervolumeValue"].Value).Value; }
            set { ((DoubleValue)Results["HipervolumeValue"].Value).Value = value; }
        }

        private DoubleMatrix ResultsBestFront
        {
            get { return (DoubleMatrix)Results["Best Front"].Value; }
            set { Results["Best Front"].Value = value; }
        }

        private int ResultsEvaluations
        {
            get { return ((IntValue)Results["Evaluations"].Value).Value; }
            set { ((IntValue)Results["Evaluations"].Value).Value = value; }
        }
        private int ResultsIterations
        {
            get { return ((IntValue)Results["Iterations"].Value).Value; }
            set { ((IntValue)Results["Iterations"].Value).Value = value; }
        }

        private DataTable ResultsQualities
        {
            get { return ((DataTable)Results["Qualities"].Value); }
        }
        private DataRow ResultsQualitiesBest
        {
            get { return ResultsQualities.Rows["Best Quality"]; }
        }

        #endregion

        [Storable]
        private RankBasedParetoFrontAnalyzer paretoFrontAnalyzer;

        [StorableConstructor]
        protected GDE3(bool deserializing) : base(deserializing) { }

        protected GDE3(GDE3 original, Cloner cloner)
          : base(original, cloner)
        {
            paretoFrontAnalyzer = (RankBasedParetoFrontAnalyzer)cloner.Clone(original.paretoFrontAnalyzer);
        }

        public override IDeepCloneable Clone(Cloner cloner)
        {
            return new GDE3(this, cloner);
        }

        public GDE3()
        {
            Parameters.Add(new FixedValueParameter<IntValue>(MaximumEvaluationsParameterName, "", new IntValue(Int32.MaxValue)));
            Parameters.Add(new ValueParameter<IntValue>(PopulationSizeParameterName, "The size of the population of solutions.", new IntValue(100)));
            Parameters.Add(new ValueParameter<DoubleValue>(CrossoverProbabilityParameterName, "The value for crossover rate", new DoubleValue(0.5)));
            Parameters.Add(new ValueParameter<DoubleValue>(ScalingFactorParameterName, "The value for scaling factor", new DoubleValue(0.5)));
            Parameters.Add(new LookupParameter<DoubleMatrix>("BestKnownFront", "The currently best known Pareto front"));
        }

        protected override void Run(CancellationToken cancellationToken)
        {
            // Set up the results display
            Results.Add(new Result("Iterations", new IntValue(0)));
            Results.Add(new Result("Evaluations", new IntValue(0)));
            Results.Add(new Result("Best Front", new DoubleMatrix()));
            Results.Add(new Result("InvertedGenerationalDistance", new DoubleValue(0)));
            Results.Add(new Result("HipervolumeValue", new DoubleValue(0)));
            Results.Add(new Result("Scatterplot", typeof(IMOFrontModel)));
            var table = new DataTable("Qualities");
            table.Rows.Add(new DataRow("Best Quality"));
            Results.Add(new Result("Qualities", table));

            //create initial population
            var population = new List<SolutionSet>(PopulationSizeParameter.Value.Value);
            var offspringPopulation = new List<SolutionSet>(PopulationSizeParameter.Value.Value);
            var matingPopulation = new List<SolutionSet>();

            for (int i = 0; i < PopulationSizeParameter.Value.Value; ++i)
            {
                var m = createIndividual();
                population.Add(m);
                var n = createEmptyIndividual();
                offspringPopulation.Add(n);
            }

            this.initProgress();
            int iterations = 1;

            while (ResultsEvaluations < MaximumEvaluations
               && !cancellationToken.IsCancellationRequested)
            {
                matingPopulation = selection(population);
                offspringPopulation = reproduction(matingPopulation, population, offspringPopulation);

                population = replacement(population, offspringPopulation);

                double[][] qualitiesFinal = new double[PopulationSizeParameter.Value.Value][];

                for (int i = 0; i < PopulationSizeParameter.Value.Value; ++i)
                {
                    qualitiesFinal[i] = new double[Problem.Objectives];
                    for (int j = 0; j < Problem.Objectives; ++j)
                    {
                        qualitiesFinal[i][j] = population[i].Quality[j];
                    }
                }

                double[][] populationFinal = new double[PopulationSizeParameter.Value.Value][];

                for (int i = 0; i < PopulationSizeParameter.Value.Value; ++i)
                {
                    populationFinal[i] = new double[Problem.ProblemSize];
                    for (int j = 0; j < Problem.ProblemSize; ++j)
                    {
                        populationFinal[i][j] = population[i].Population[j];
                    }
                }

                int objectives = Problem.Objectives;
                var optimalfront = Problem.TestFunction.OptimalParetoFront(objectives);

                IEnumerable<double[]> en = qualitiesFinal;

                double[][] opf = new double[0][];
                if (optimalfront != null)
                {
                    opf = optimalfront.Select(s => s.ToArray()).ToArray();
                }

                iterations = iterations + 1;

                IEnumerable<double[]> front = NonDominatedSelect.selectNonDominatedVectors(qualitiesFinal, Problem.TestFunction.Maximization(objectives), true);
                //update the results
                ResultsEvaluations = evals;
                ResultsIterations = iterations;
                ResultsBestFront = new DoubleMatrix(MultiObjectiveTestFunctionProblem.To2D(qualitiesFinal));
                ResultsGenerationalDistance = InvertedGenerationalDistance.Calculate(en, optimalfront, 1);
                ResultsHypervolume = Hypervolume.Calculate(front, ReferencePoint(objectives), Problem.TestFunction.Maximization(objectives));

                Results["Scatterplot"].Value = new MOSolution(qualitiesFinal, populationFinal, opf, objectives);
            }
        }
       
        protected SolutionSet createIndividual()
        {
            var dim = Problem.ProblemSize;
            var lb = Problem.Bounds[0, 0];
            var ub = Problem.Bounds[0, 1];
            var range = ub - lb;
            var v = new double[Problem.ProblemSize];
            SolutionSet solutionObject = new SolutionSet(PopulationSizeParameter.Value.Value);

            for (int i = 0; i < Problem.ProblemSize; ++i)
            {
                v[i] = _random.NextDouble() * range + lb;

            }
            var m = new RealVector(v);
            solutionObject.Population = m;
            solutionObject.Quality = Problem.Evaluate(m, _random);
            return solutionObject;
        }

        private SolutionSet createEmptyIndividual()
        {
            SolutionSet solutionObject = new SolutionSet(PopulationSizeParameter.Value.Value);
            var n = new RealVector(Problem.ProblemSize);
            solutionObject.Population = n;
            return solutionObject;
        }

        protected void initProgress()
        {
            this.evals = PopulationSizeParameter.Value.Value;
        }

        protected void updateProgres()
        {
            this.evals += PopulationSizeParameter.Value.Value;
        }

        protected List<SolutionSet> selection(List<SolutionSet> population)
        {
            List<SolutionSet> parents = new List<SolutionSet>();
            for (int i = 0; i < PopulationSizeParameter.Value.Value; i++)
            {
                int r0, r1, r2;
                //assure the selected vectors r0, r1 and r2 are different
                do
                {
                    r0 = _random.Next(0, PopulationSizeParameter.Value.Value);
                } while (r0 == i);
                do
                {
                    r1 = _random.Next(0, PopulationSizeParameter.Value.Value);
                } while (r1 == i || r1 == r0);
                do
                {
                    r2 = _random.Next(0, PopulationSizeParameter.Value.Value);
                } while (r2 == i || r2 == r0 || r2 == r1);

                parents.Add(population[r0]);
                parents.Add(population[r1]);
                parents.Add(population[r2]);
            }
            return parents;
        }

        protected List<SolutionSet> reproduction(List<SolutionSet> matingPopulation, List<SolutionSet> parentPopulation, List<SolutionSet> offspringPopulation)
        {
            for (int i = 0; i < PopulationSizeParameter.Value.Value; i++)
            {
                List<SolutionSet> parents = new List<SolutionSet>(3);
                for (int j = 0; j < 3; j++) {
                    parents.Add(matingPopulation[0]);
                    matingPopulation.RemoveAt(0);
                }

                double rnbr = _random.Next(0, Problem.ProblemSize);
                for (int m = 0; m < Problem.ProblemSize; m++)
                {
                    if (_random.NextDouble() < CrossoverProbabilityParameter.Value.Value || m == rnbr)
                    {
                        offspringPopulation[i].Population[m] = parents[2].Population[m] +
                            ScalingFactorParameter.Value.Value * (parents[1].Population[m] - parents[2].Population[0]);
                        //check the problem upper and lower bounds
                        if (offspringPopulation[i].Population[m] > Problem.Bounds[0, 1]) offspringPopulation[i].Population[m] = Problem.Bounds[0, 1];
                        if (offspringPopulation[i].Population[m] < Problem.Bounds[0, 0]) offspringPopulation[i].Population[m] = Problem.Bounds[0, 0];
                    }
                    else
                    {
                        offspringPopulation[i].Population[m] = parentPopulation[i].Population[m]; // check again
                    }
                }
                offspringPopulation[i].Quality = Problem.Evaluate(offspringPopulation[i].Population, _random);
            }
            this.updateProgres();
            return offspringPopulation;
        }

        private List<SolutionSet> replacement(List<SolutionSet> population, List<SolutionSet> offspringPopulation)
        {
            List<SolutionSet> tmpList = new List<SolutionSet>();
            for (int i = 0; i < PopulationSizeParameter.Value.Value; i++)
            {
                int result;
                result = compareDomination(population[i].Quality, offspringPopulation[i].Quality);
                if (result == -1)
                { // Solution i dominates child
                    tmpList.Add(population[i]);
                }
                else if (result == 1)
                { // child dominates
                    tmpList.Add(offspringPopulation[i]);
                }
                else
                { // the two solutions are non-dominated
                    tmpList.Add(offspringPopulation[i]);
                    tmpList.Add(population[i]);
                }
            }

            List<SolutionSet>[]ranking = computeRanking(tmpList);
            List<SolutionSet> finalPopulation = crowdingDistanceSelection(ranking);
            return finalPopulation;
        }

        private List<SolutionSet> crowdingDistanceSelection(List<SolutionSet>[] ranking)
        {
            List<SolutionSet> population = new List<SolutionSet>();
            int rankingIndex = 0;
            while (populationIsNotFull(population))
            {
                if (subFrontFillsIntoThePopulation(ranking, rankingIndex, population)){
                    addRankedSolutionToPopulation(ranking, rankingIndex, population);
                    rankingIndex++;
                } else {
                    crowdingDistanceAssignment(ranking[rankingIndex], rankingIndex);
                    addLastRankedSolutionToPopulation(ranking, rankingIndex, population);
                }
            }
            return population;
        }

        private void addLastRankedSolutionToPopulation(List<SolutionSet>[] ranking, int rankingIndex, List<SolutionSet> population)
        {
            List<SolutionSet> currentRankedFront = ranking[rankingIndex];
            currentRankedFront.Sort((x, y) => x.CrowdingDistance.CompareTo(y.CrowdingDistance));
            int i = 0;
            while (population.Count < PopulationSizeParameter.Value.Value)
            {
                population.Add(currentRankedFront[i]);
                i++;
            }
        }

        public void crowdingDistanceAssignment(List<SolutionSet> rankingSubfront, int index)
        {
            int size = rankingSubfront.Count;

            if (size == 0)
                return;

            if (size == 1)
            {
                rankingSubfront[0].CrowdingDistance = double.PositiveInfinity;
                return;
            }

            if (size == 2)
            {
                rankingSubfront[0].CrowdingDistance = double.PositiveInfinity;
                rankingSubfront[1].CrowdingDistance = double.PositiveInfinity;
                return;
            }

            //Use a new SolutionSet to evite alter original solutionSet
            List<SolutionSet> front = new List<SolutionSet>(size);
            for (int i = 0; i < size; i++)
            {
                front.Add(rankingSubfront[i]);
            }

            for (int i = 0; i < size; i++)
                rankingSubfront[i].CrowdingDistance = 0.0;

            double objetiveMaxn;
            double objetiveMinn;
            double distance;

            for (int i = 0; i < Problem.Objectives ; i++)
            {
                // Sort the population by Obj n            
                front.Sort((x, y) => x.Quality[i].CompareTo(y.Quality[i]));
                objetiveMinn = front[0].Quality[i];
                objetiveMaxn = front[front.Count - 1].Quality[i];

                //Set de crowding distance            
                front[0].CrowdingDistance = double.PositiveInfinity; 
                front[size - 1].CrowdingDistance = double.PositiveInfinity;

                for (int j = 1; j < size - 1; j++)
                {
                    distance = front[j + 1].Quality[i] - front[j - 1].Quality[i];
                    distance = distance / (objetiveMaxn - objetiveMinn);
                    distance += front[j].CrowdingDistance;
                    front[j].CrowdingDistance = distance;
                }
            }
        }

        private void addRankedSolutionToPopulation(List<SolutionSet>[] ranking, int rankingIndex, List<SolutionSet> population)
        {
            foreach(SolutionSet solution in ranking[rankingIndex])
            {
                population.Add(solution);
            }
        }

        private bool subFrontFillsIntoThePopulation(List<SolutionSet>[] ranking, int rankingIndex, List<SolutionSet> population)
        {
            return ranking[rankingIndex].Count < (PopulationSizeParameter.Value.Value - population.Count);
        }

        private bool populationIsNotFull(List<SolutionSet> population)
        {
            return population.Count < PopulationSizeParameter.Value.Value;
        }

        private List<SolutionSet>[] computeRanking(List<SolutionSet> tmpList)
        {
            // dominateMe[i] contains the number of solutions dominating i        
            int[] dominateMe = new int[tmpList.Count];

            // iDominate[k] contains the list of solutions dominated by k
            List<int>[] iDominate = new List<int>[tmpList.Count];

            // front[i] contains the list of individuals belonging to the front i
            List<int>[] front = new List<int>[tmpList.Count + 1];
            int[] solutionsRanks = new int[tmpList.Count];

            // flagDominate is an auxiliar encodings.variable
            int flagDominate;

            // Initialize the fronts 
            for (int i = 0; i < front.Length; i++)
            {
                front[i] = new List<int>();
            }

            //-> Fast non dominated sorting algorithm
            // Contribution of Guillaume Jacquenot
            for (int p = 0; p < tmpList.Count; p++)
            {
                // Initialize the list of individuals that i dominate and the number
                // of individuals that dominate me
                iDominate[p] = new List<int>();
                dominateMe[p] = 0;
            }
            for (int p = 0; p < (tmpList.Count - 1); p++)
            {
                // For all q individuals , calculate if p dominates q or vice versa
                for (int q = p + 1; q < tmpList.Count; q++)
                {
                    flagDominate = compareDomination(tmpList[p].Quality, tmpList[q].Quality);
                    if (flagDominate == -1)
                    {
                        iDominate[p].Add(q);
                        dominateMe[q]++;
                    }
                    else if (flagDominate == 1)
                    {
                        iDominate[q].Add(p);
                        dominateMe[p]++;
                    }
                }
                // If nobody dominates p, p belongs to the first front
            }
            for (int i = 0; i < tmpList.Count; i++)
            {
                if (dominateMe[i] == 0)
                {
                    front[0].Add(i);
                    tmpList[i].Rank = 0;
                }
            }

            //Obtain the rest of fronts
            int k = 0;

            while (front[k].Count != 0)
            {
                k++;
                foreach (var it1 in front[k - 1])
                {
                    foreach (var it2 in iDominate[it1])
                    {
                        int index = it2;
                        dominateMe[index]--;
                        if (dominateMe[index] == 0)
                        {
                            front[k].Add(index);
                            tmpList[index].Rank = k;
                        }
                    }
                }
            }
            //<-

            var rankedSubpopulation = new List<SolutionSet>[k];
            //0,1,2,....,i-1 are front, then i fronts
            for (int j = 0; j < k; j++)
            {
                rankedSubpopulation[j] = new List<SolutionSet>(front[j].Count);
                foreach (var it1 in front[j])
                {
                    rankedSubpopulation[j].Add(tmpList[it1]);
                }
            }
            return rankedSubpopulation;
        }

        private int compareDomination(double[] solution1, double[] solution2)
        {
            int dominate1; // dominate1 indicates if some objective of solution1 
                           // dominates the same objective in solution2. dominate2
            int dominate2; // is the complementary of dominate1.

            dominate1 = 0;
            dominate2 = 0;

            int flag; //stores the result of the comparison

            // Equal number of violated constraints. Applying a dominance Test then
            double value1, value2;
            for (int i = 0; i < Problem.Objectives; i++)
            {
                value1 = solution1[i];
                value2 = solution2[i];
                if (value1 < value2)
                {
                    flag = -1;
                }
                else if (value1 > value2)
                {
                    flag = 1;
                }
                else
                {
                    flag = 0;
                }

                if (flag == -1)
                {
                    dominate1 = 1;
                }

                if (flag == 1)
                {
                    dominate2 = 1;
                }
            }

            if (dominate1 == dominate2)
            {
                return 0; //No one dominate the other
            }
            if (dominate1 == 1)
            {
                return -1; // solution1 dominate
            }
            return 1;    // solution2 dominate   
        }
    }
}