1 | using System;
|
---|
2 | using System.Collections.Generic;
|
---|
3 | using System.Linq;
|
---|
4 | using HeuristicLab.Core;
|
---|
5 | using HeuristicLab.Optimization;
|
---|
6 |
|
---|
7 | namespace HeuristicLab.Algorithms.NSGA3
|
---|
8 | {
|
---|
9 | public static class NSGA3Selection
|
---|
10 | {
|
---|
11 | private const double EPSILON = 10e-6; // a tiny number that is greater than 0
|
---|
12 |
|
---|
13 | /// <summary>
|
---|
14 | /// TODO: Description If the fitnesses of all solutions in <paramref name="rt" /> do not
|
---|
15 | /// have the same number of objectives, undefined behavior happens.
|
---|
16 | /// </summary>
|
---|
17 | /// <param name="rt"></param>
|
---|
18 | /// <param name="referencePoints"></param>
|
---|
19 | /// <param name="maximization"></param>
|
---|
20 | /// <returns></returns>
|
---|
21 | public static List<Solution> SelectSolutionsForNextGeneration(List<Solution> rt, List<ReferencePoint> referencePoints, bool[] maximization, IRandom random)
|
---|
22 | {
|
---|
23 | int N = rt.Count / 2; // number of solutions in a generation
|
---|
24 | int numberOfObjectives = rt.First().Fitness.Length;
|
---|
25 | List<Solution> st = new List<Solution>();
|
---|
26 |
|
---|
27 | // Do non-dominated sort
|
---|
28 | var qualities = Utility.ToFitnessMatrix(rt);
|
---|
29 | // compute the pareto fronts using the DominationCalculator and discard the qualities
|
---|
30 | // part in the inner tuples
|
---|
31 | var fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, maximization, out int[] rank, true)
|
---|
32 | .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
|
---|
33 |
|
---|
34 | int lf = -1; // last front to be included
|
---|
35 | for (int i = 0; i < fronts.Count && st.Count < N; i++)
|
---|
36 | {
|
---|
37 | st = Utility.Concat(st, fronts[i]);
|
---|
38 | lf = i;
|
---|
39 | }
|
---|
40 | // remove useless fronts
|
---|
41 | fronts.RemoveRange(lf + 1, fronts.Count - lf - 1);
|
---|
42 |
|
---|
43 | List<Solution> nextGeneration = new List<Solution>();
|
---|
44 | if (st.Count == N) // no special selection needs to be done
|
---|
45 | nextGeneration = st;
|
---|
46 | else
|
---|
47 | {
|
---|
48 | // Add every front to the next generation list except for the one that is added only partially
|
---|
49 | nextGeneration = new List<Solution>();
|
---|
50 | for (int f = 0; f < lf; f++)
|
---|
51 | nextGeneration = Utility.Concat(nextGeneration, fronts[f]);
|
---|
52 | int k = N - nextGeneration.Count;
|
---|
53 | Normalize(st, numberOfObjectives);
|
---|
54 | Associate(fronts, referencePoints);
|
---|
55 | List<Solution> solutionsToAdd = Niching(k, referencePoints, random);
|
---|
56 | nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
|
---|
57 | }
|
---|
58 |
|
---|
59 | return nextGeneration;
|
---|
60 | }
|
---|
61 |
|
---|
62 | private static void Normalize(List<Solution> st, int numberOfObjectives)
|
---|
63 | {
|
---|
64 | // Find the ideal point
|
---|
65 | double[] idealPoint = new double[numberOfObjectives];
|
---|
66 | for (int j = 0; j < numberOfObjectives; j++)
|
---|
67 | {
|
---|
68 | // Compute ideal point
|
---|
69 | idealPoint[j] = Utility.Min(s => s.Fitness[j], st);
|
---|
70 |
|
---|
71 | // Translate objectives
|
---|
72 | foreach (var solution in st)
|
---|
73 | solution.Fitness[j] -= idealPoint[j];
|
---|
74 | }
|
---|
75 |
|
---|
76 | // Find the extreme points
|
---|
77 | Solution[] extremePoints = new Solution[numberOfObjectives];
|
---|
78 | for (int j = 0; j < numberOfObjectives; j++)
|
---|
79 | {
|
---|
80 | // Compute extreme points
|
---|
81 | double[] weights = new double[numberOfObjectives];
|
---|
82 | for (int i = 0; i < numberOfObjectives; i++) weights[i] = EPSILON;
|
---|
83 | weights[j] = 1;
|
---|
84 |
|
---|
85 | extremePoints[j] = Utility.ArgMin(s => ASF(s.Fitness, weights), st);
|
---|
86 | }
|
---|
87 |
|
---|
88 | // Compute intercepts
|
---|
89 | List<double> intercepts = GetIntercepts(extremePoints.ToList());
|
---|
90 |
|
---|
91 | // Normalize objectives
|
---|
92 | NormalizeObjectives(st, intercepts, idealPoint);
|
---|
93 | }
|
---|
94 |
|
---|
95 | private static void NormalizeObjectives(List<Solution> st, List<double> intercepts, double[] idealPoint)
|
---|
96 | {
|
---|
97 | int numberOfObjectives = st.First().Fitness.Length;
|
---|
98 |
|
---|
99 | foreach (var solution in st)
|
---|
100 | {
|
---|
101 | for (int i = 0; i < numberOfObjectives; i++)
|
---|
102 | {
|
---|
103 | if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
|
---|
104 | {
|
---|
105 | solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
|
---|
106 | }
|
---|
107 | else
|
---|
108 | {
|
---|
109 | solution.Fitness[i] = solution.Fitness[i] / EPSILON;
|
---|
110 | }
|
---|
111 | }
|
---|
112 | }
|
---|
113 | }
|
---|
114 |
|
---|
115 | private static void Associate(List<List<Solution>> fronts, List<ReferencePoint> referencePoints)
|
---|
116 | {
|
---|
117 | for (int f = 0; f < fronts.Count; f++)
|
---|
118 | {
|
---|
119 | foreach (var solution in fronts[f])
|
---|
120 | {
|
---|
121 | // find reference point for which the perpendicular distance to the current
|
---|
122 | // solution is the lowest
|
---|
123 | var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);
|
---|
124 | // associated reference point
|
---|
125 | var arp = rpAndDist.Item1;
|
---|
126 | // distance to that reference point
|
---|
127 | var dist = rpAndDist.Item2;
|
---|
128 |
|
---|
129 | if (f + 1 != fronts.Count)
|
---|
130 | {
|
---|
131 | // Todo: Add member for reference point on index min_rp
|
---|
132 | arp.NumberOfAssociatedSolutions++;
|
---|
133 | }
|
---|
134 | else
|
---|
135 | {
|
---|
136 | // Todo: Add potential member for reference point on index min_rp
|
---|
137 | arp.AddPotentialAssociatedSolution(solution, dist);
|
---|
138 | }
|
---|
139 | }
|
---|
140 | }
|
---|
141 | }
|
---|
142 |
|
---|
143 | private static List<Solution> Niching(int k, List<ReferencePoint> referencePoints, IRandom random)
|
---|
144 | {
|
---|
145 | List<Solution> solutions = new List<Solution>();
|
---|
146 | while (solutions.Count != k)
|
---|
147 | {
|
---|
148 | ReferencePoint min_rp = FindNicheReferencePoint(referencePoints, random);
|
---|
149 |
|
---|
150 | Solution chosen = SelectClusterMember(min_rp);
|
---|
151 | if (chosen == null)
|
---|
152 | {
|
---|
153 | referencePoints.Remove(min_rp);
|
---|
154 | }
|
---|
155 | else
|
---|
156 | {
|
---|
157 | min_rp.NumberOfAssociatedSolutions++;
|
---|
158 | min_rp.RemovePotentialAssociatedSolution(chosen);
|
---|
159 | solutions.Add(chosen);
|
---|
160 | }
|
---|
161 | }
|
---|
162 |
|
---|
163 | return solutions;
|
---|
164 | }
|
---|
165 |
|
---|
166 | private static ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints, IRandom random)
|
---|
167 | {
|
---|
168 | // the minimum number of associated solutions for a reference point over all reference points
|
---|
169 | int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);
|
---|
170 |
|
---|
171 | // the reference points that share the number of associated solutions where that number
|
---|
172 | // is equal to minNumber
|
---|
173 | List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();
|
---|
174 | foreach (var referencePoint in referencePoints)
|
---|
175 | if (referencePoint.NumberOfAssociatedSolutions == minNumber)
|
---|
176 | minAssociatedReferencePoints.Add(referencePoint);
|
---|
177 |
|
---|
178 | if (minAssociatedReferencePoints.Count > 1)
|
---|
179 | return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];
|
---|
180 | else
|
---|
181 | return minAssociatedReferencePoints.Single();
|
---|
182 | }
|
---|
183 |
|
---|
184 | private static Solution SelectClusterMember(ReferencePoint referencePoint)
|
---|
185 | {
|
---|
186 | Solution chosen = null;
|
---|
187 | if (referencePoint.HasPotentialMember())
|
---|
188 | {
|
---|
189 | if (referencePoint.NumberOfAssociatedSolutions == 0)
|
---|
190 | chosen = referencePoint.FindClosestMember();
|
---|
191 | else
|
---|
192 | chosen = referencePoint.RandomMember();
|
---|
193 | }
|
---|
194 | return chosen;
|
---|
195 | }
|
---|
196 |
|
---|
197 | private static double GetPerpendicularDistance(double[] values, double[] fitness)
|
---|
198 | {
|
---|
199 | double numerator = 0;
|
---|
200 | double denominator = 0;
|
---|
201 | for (int i = 0; i < values.Length; i++)
|
---|
202 | {
|
---|
203 | numerator += values[i] * fitness[i];
|
---|
204 | denominator += Math.Pow(values[i], 2);
|
---|
205 | }
|
---|
206 | double k = numerator / denominator;
|
---|
207 |
|
---|
208 | double d = 0;
|
---|
209 | for (int i = 0; i < values.Length; i++)
|
---|
210 | {
|
---|
211 | d += Math.Pow(k * values[i] - fitness[i], 2);
|
---|
212 | }
|
---|
213 | return Math.Sqrt(d);
|
---|
214 | }
|
---|
215 |
|
---|
216 | private static double ASF(double[] x, double[] weight)
|
---|
217 | {
|
---|
218 | int numberOfObjectives = x.Length;
|
---|
219 | List<int> dimensions = new List<int>();
|
---|
220 | for (int i = 0; i < numberOfObjectives; i++)
|
---|
221 | dimensions.Add(i);
|
---|
222 |
|
---|
223 | return Utility.Max((dim) => x[dim] / weight[dim], dimensions);
|
---|
224 | }
|
---|
225 |
|
---|
226 | private static List<double> GetIntercepts(List<Solution> extremePoints)
|
---|
227 | {
|
---|
228 | int numberOfObjectives = extremePoints.First().Fitness.Length;
|
---|
229 |
|
---|
230 | // Check whether there are duplicate extreme points. This might happen but the original
|
---|
231 | // paper does not mention how to deal with it.
|
---|
232 | bool duplicate = false;
|
---|
233 | for (int i = 0; !duplicate && i < extremePoints.Count; i++)
|
---|
234 | {
|
---|
235 | for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
|
---|
236 | {
|
---|
237 | // maybe todo: override Equals method of solution?
|
---|
238 | duplicate = extremePoints[i].Equals(extremePoints[j]);
|
---|
239 | }
|
---|
240 | }
|
---|
241 |
|
---|
242 | List<double> intercepts = new List<double>();
|
---|
243 |
|
---|
244 | if (duplicate)
|
---|
245 | { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
|
---|
246 | for (int f = 0; f < numberOfObjectives; f++)
|
---|
247 | {
|
---|
248 | // extreme_points[f] stands for the individual with the largest value of
|
---|
249 | // objective f
|
---|
250 | intercepts.Add(extremePoints[f].Fitness[f]);
|
---|
251 | }
|
---|
252 | }
|
---|
253 | else
|
---|
254 | {
|
---|
255 | // Find the equation of the hyperplane
|
---|
256 | List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
|
---|
257 | for (int i = 0; i < numberOfObjectives; i++)
|
---|
258 | {
|
---|
259 | b.Add(1.0);
|
---|
260 | }
|
---|
261 |
|
---|
262 | List<List<double>> a = new List<List<double>>();
|
---|
263 | foreach (Solution s in extremePoints)
|
---|
264 | {
|
---|
265 | List<double> aux = new List<double>();
|
---|
266 | for (int i = 0; i < numberOfObjectives; i++)
|
---|
267 | aux.Add(s.Fitness[i]);
|
---|
268 | a.Add(aux);
|
---|
269 | }
|
---|
270 | List<double> x = GaussianElimination(a, b);
|
---|
271 |
|
---|
272 | // Find intercepts
|
---|
273 | for (int f = 0; f < numberOfObjectives; f++)
|
---|
274 | {
|
---|
275 | intercepts.Add(1.0 / x[f]);
|
---|
276 | }
|
---|
277 | }
|
---|
278 |
|
---|
279 | return intercepts;
|
---|
280 | }
|
---|
281 |
|
---|
282 | private static List<double> GaussianElimination(List<List<double>> a, List<double> b)
|
---|
283 | {
|
---|
284 | List<double> x = new List<double>();
|
---|
285 |
|
---|
286 | int n = a.Count;
|
---|
287 | for (int i = 0; i < n; i++)
|
---|
288 | a[i].Add(b[i]);
|
---|
289 |
|
---|
290 | for (int @base = 0; @base < n - 1; @base++)
|
---|
291 | for (int target = @base + 1; target < n; target++)
|
---|
292 | {
|
---|
293 | double ratio = a[target][@base] / a[@base][@base];
|
---|
294 | for (int term = 0; term < a[@base].Count; term++)
|
---|
295 | a[target][term] = a[target][term] - a[@base][term] * ratio;
|
---|
296 | }
|
---|
297 |
|
---|
298 | for (int i = 0; i < n; i++)
|
---|
299 | x.Add(0.0);
|
---|
300 |
|
---|
301 | for (int i = n - 1; i >= 0; i--)
|
---|
302 | {
|
---|
303 | for (int known = i + 1; known < n; known++)
|
---|
304 | a[i][n] = a[i][n] - a[i][known] * x[known];
|
---|
305 | x[i] = a[i][n] / a[i][i];
|
---|
306 | }
|
---|
307 |
|
---|
308 | return x;
|
---|
309 | }
|
---|
310 | }
|
---|
311 | } |
---|