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; using HeuristicLab.Algorithms.GDE3; 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 override Type ProblemType { get { return typeof(MultiObjectiveTestFunctionProblem); } } public new MultiObjectiveTestFunctionProblem Problem { get { return (MultiObjectiveTestFunctionProblem)base.Problem; } set { base.Problem = value; } } public ILookupParameter BestKnownFrontParameter { get { return (ILookupParameter)Parameters["BestKnownFront"]; } } private readonly IRandom _random = new MersenneTwister(); private int evals; #region ParameterNames private const string MaximumGenerationsParameterName = "Maximum Generations"; private const string CrossoverProbabilityParameterName = "CrossoverProbability"; private const string PopulationSizeParameterName = "PopulationSize"; private const string ScalingFactorParameterName = "ScalingFactor"; #endregion #region ParameterProperties public IFixedValueParameter MaximumGenerationsParameter { get { return (IFixedValueParameter)Parameters[MaximumGenerationsParameterName]; } } private ValueParameter PopulationSizeParameter { get { return (ValueParameter)Parameters[PopulationSizeParameterName]; } } public ValueParameter CrossoverProbabilityParameter { get { return (ValueParameter)Parameters[CrossoverProbabilityParameterName]; } } public ValueParameter ScalingFactorParameter { get { return (ValueParameter)Parameters[ScalingFactorParameterName]; } } #endregion #region Properties public int MaximumEvaluations { get { return MaximumGenerationsParameter.Value.Value; } set { MaximumGenerationsParameter.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 ResultsInvertedGenerationalDistance { get { return ((DoubleValue)Results["InvertedGenerationalDistance"].Value).Value; } set { ((DoubleValue)Results["InvertedGenerationalDistance"].Value).Value = value; } } private double ResultsHypervolume { get { return ((DoubleValue)Results["HyperVolumeValue"].Value).Value; } set { ((DoubleValue)Results["HyperVolumeValue"].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 ResultsGenerations { get { return ((IntValue)Results["Generations"].Value).Value; } set { ((IntValue)Results["Generations"].Value).Value = value; } } private double ResultsGenerationalDistance { get { return ((DoubleValue)Results["GenerationalDistance"].Value).Value; } set { ((DoubleValue)Results["GenerationalDistance"].Value).Value = value; } } private int ResultsIterations { get { return ((IntValue)Results["Iterations"].Value).Value; } set { ((IntValue)Results["Iterations"].Value).Value = value; } } private double ResultsSpacing { get { return ((DoubleValue)Results["Spacing"].Value).Value; } set { ((DoubleValue)Results["Spacing"].Value).Value = value; } } private double ResultsCrowding { get { return ((DoubleValue)Results["Crowding"].Value).Value; } set { ((DoubleValue)Results["Crowding"].Value).Value = value; } } #endregion [StorableConstructor] protected GDE3(bool deserializing) : base(deserializing) { } protected GDE3(GDE3 original, Cloner cloner) : base(original, cloner) { } public override IDeepCloneable Clone(Cloner cloner) { return new GDE3(this, cloner); } public GDE3() { Parameters.Add(new FixedValueParameter(MaximumGenerationsParameterName, "", new IntValue(1000))); Parameters.Add(new ValueParameter(PopulationSizeParameterName, "The size of the population of solutions.", new IntValue(100))); Parameters.Add(new ValueParameter(CrossoverProbabilityParameterName, "The value for crossover rate", new DoubleValue(0.5))); Parameters.Add(new ValueParameter(ScalingFactorParameterName, "The value for scaling factor", new DoubleValue(0.5))); Parameters.Add(new LookupParameter("BestKnownFront", "The currently best known Pareto front")); } protected override void Run(CancellationToken cancellationToken) { // Set up the results display Results.Add(new Result("Generations", new IntValue(0))); Results.Add(new Result("Evaluations", new IntValue(0))); Results.Add(new Result("Best Front", new DoubleMatrix())); Results.Add(new Result("Crowding", new DoubleValue(0))); Results.Add(new Result("InvertedGenerationalDistance", new DoubleValue(0))); Results.Add(new Result("GenerationalDistance", new DoubleValue(0))); Results.Add(new Result("HyperVolumeValue", new DoubleValue(0))); Results.Add(new Result("Spacing", 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)); //setup the variables List population; List offspringPopulation; SolutionSet[] parent; //initialize population population = new List(PopulationSizeParameter.Value.Value); for (int i = 0; i < PopulationSizeParameter.Value.Value; ++i) { var m = createIndividual(); population.Add(m); } this.initProgress(); int iterations = 1; while (ResultsGenerations < MaximumGenerationsParameter.Value.Value && !cancellationToken.IsCancellationRequested) { var populationSize = PopulationSizeParameter.Value.Value; // Create the offSpring solutionSet offspringPopulation = new List(PopulationSizeParameter.Value.Value * 2); for (int i = 0; i < populationSize; i++) { // Obtain parents. Two parameters are required: the population and the // index of the current individual parent = selection(population, i); SolutionSet child; // Crossover. The parameters are the current individual and the index of the array of parents child = reproduction(population[i], parent); child.Quality = Problem.Evaluate(child.Population, _random); this.updateProgres(); // Dominance test int result; result = compareDomination(population[i].Quality, child.Quality); if (result == -1) { // Solution i dominates child offspringPopulation.Add(population[i]); } else if (result == 1) { // child dominates offspringPopulation.Add(child); } else { // the two solutions are non-dominated offspringPopulation.Add(child); offspringPopulation.Add(population[i]); } } // Ranking the offspring population List[] ranking = computeRanking(offspringPopulation); int remain = populationSize; int index = 0; population.Clear(); List front = null; // Obtain the next front front = ranking[index]; while ((remain > 0) && (remain >= front.Count)) { //Assign crowding distance to individuals crowdingDistanceAssignment(front, index); //Add the individuals of this front for (int k = 0; k < front.Count; k++) { population.Add(front[k]); } // for //Decrement remain remain = remain - front.Count; //Obtain the next front index++; if (remain > 0) { front = ranking[index]; } } // remain is less than front(index).size, insert only the best one if (remain > 0) { // front contains individuals to insert while (front.Count > remain) { crowdingDistanceAssignment(front, index); int indx = getWorstIndex(front); front.RemoveAt(indx); } for (int k = 0; k < front.Count; k++) { population.Add(front[k]); } remain = 0; } iterations++; ResultsGenerations = iterations; displayResults(front); } } private void displayResults(List population) { //compute the 0 front // Return the first non-dominated front List[] rankingFinal = computeRanking(population); int objectives = Problem.Objectives; var optimalfront = Problem.TestFunction.OptimalParetoFront(objectives); double[][] opf = new double[0][]; if (optimalfront != null) { opf = optimalfront.Select(s => s.ToArray()).ToArray(); } double[][] qualitiesFinal = new double[rankingFinal[0].Count][]; for (int i = 0; i < rankingFinal[0].Count; ++i) { qualitiesFinal[i] = new double[Problem.Objectives]; for (int j = 0; j < Problem.Objectives; ++j) { qualitiesFinal[i][j] = rankingFinal[0][i].Quality[j]; } } double[][] populationFinal = new double[rankingFinal[0].Count][]; for (int i = 0; i < rankingFinal[0].Count; ++i) { populationFinal[i] = new double[Problem.ProblemSize]; for (int j = 0; j < Problem.ProblemSize; ++j) { populationFinal[i][j] = rankingFinal[0][i].Population[j]; } } IEnumerable en = qualitiesFinal; IEnumerable frontVectors = NonDominatedSelect.selectNonDominatedVectors(qualitiesFinal, Problem.TestFunction.Maximization(objectives), true); //update the results ResultsEvaluations = this.evals; ResultsBestFront = new DoubleMatrix(MultiObjectiveTestFunctionProblem.To2D(qualitiesFinal)); ResultsCrowding = Crowding.Calculate(qualitiesFinal, Problem.TestFunction.Bounds(objectives)); ResultsInvertedGenerationalDistance = InvertedGenerationalDistance.Calculate(en, optimalfront, 1); ResultsHypervolume = Hypervolume.Calculate(frontVectors, Problem.TestFunction.ReferencePoint(objectives), Problem.TestFunction.Maximization(objectives)); ResultsGenerationalDistance = GenerationalDistance.Calculate(qualitiesFinal, optimalfront, 1); Results["Scatterplot"].Value = new MOSolution(qualitiesFinal, populationFinal, opf, objectives); ResultsSpacing = Spacing.Calculate(qualitiesFinal); } private int getWorstIndex(List SolutionsList) { int result = 0; if ((SolutionsList == null) || SolutionsList.Count == 0) { result = 0; } else { SolutionSet worstKnown = SolutionsList[0], candidateSolution; int flag; for (int i = 1; i < SolutionsList.Count; i++) { candidateSolution = SolutionsList[i]; flag = compareDomination(worstKnown.Quality, candidateSolution.Quality); if (flag == -1) { result = i; worstKnown = candidateSolution; } } } return result; } 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++; } protected SolutionSet[] selection(List population, int i) { SolutionSet[] parents = new SolutionSet[3]; 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[0] = population[r0]; parents[1] = population[r1]; parents[2] = population[r2]; return parents; } protected SolutionSet reproduction(SolutionSet parent, SolutionSet[] parentsSolutions) { var individual = createEmptyIndividual(); double rnbr = _random.Next(0, Problem.ProblemSize); for (int m = 0; m < Problem.ProblemSize; m++) { if (_random.NextDouble() < CrossoverProbabilityParameter.Value.Value || m == rnbr) { double value; value = parentsSolutions[2].Population[m] + ScalingFactorParameter.Value.Value * (parentsSolutions[0].Population[m] - parentsSolutions[1].Population[m]); //check the problem upper and lower bounds if (value > Problem.Bounds[0, 1]) value = Problem.Bounds[0, 1]; if (value < Problem.Bounds[0, 0]) value = Problem.Bounds[0, 0]; individual.Population[m] = value; } else { double value; value = parent.Population[m]; individual.Population[m] = value; } } return individual; } private List replacement(List population, List offspringPopulation) { List tmpList = new List(); 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[] ranking = computeRanking(tmpList); List finalPopulation = crowdingDistanceSelection(ranking); return finalPopulation; } private List crowdingDistanceSelection(List[] ranking) { List population = new List(); 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[] ranking, int rankingIndex, List population) { List 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 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 front = new List(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[] ranking, int rankingIndex, List population) { foreach (SolutionSet solution in ranking[rankingIndex]) { population.Add(solution); } } private bool subFrontFillsIntoThePopulation(List[] ranking, int rankingIndex, List population) { return ranking[rankingIndex].Count < (PopulationSizeParameter.Value.Value - population.Count); } private bool populationIsNotFull(List population) { return population.Count < PopulationSizeParameter.Value.Value; } private List[] computeRanking(List 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[] iDominate = new List[tmpList.Count]; // front[i] contains the list of individuals belonging to the front i List[] front = new List[tmpList.Count + 1]; // flagDominate is an auxiliar encodings.variable int flagDominate; // Initialize the fronts for (int i = 0; i < front.Length; i++) { front[i] = new List(); } //-> 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(); 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[k]; //0,1,2,....,i-1 are front, then i fronts for (int j = 0; j < k; j++) { rankedSubpopulation[j] = new List(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 } } }