001 /*
002 * Java Genetic Algorithm Library (jenetics-5.1.0).
003 * Copyright (c) 2007-2019 Franz Wilhelmstötter
004 *
005 * Licensed under the Apache License, Version 2.0 (the "License");
006 * you may not use this file except in compliance with the License.
007 * You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 *
017 * Author:
018 * Franz Wilhelmstötter (franz.wilhelmstoetter@gmail.com)
019 */
020 package io.jenetics.ext.moea;
021
022 import static java.lang.Double.POSITIVE_INFINITY;
023 import static java.lang.String.format;
024 import static java.util.Objects.requireNonNull;
025
026 import java.util.ArrayList;
027 import java.util.Arrays;
028 import java.util.Comparator;
029 import java.util.List;
030 import java.util.function.ToIntFunction;
031
032 import io.jenetics.util.ISeq;
033 import io.jenetics.util.MSeq;
034 import io.jenetics.util.ProxySorter;
035 import io.jenetics.util.Seq;
036
037 import io.jenetics.ext.internal.IntList;
038
039 /**
040 * Low-level utility methods for doing pareto-optimal calculations. This methods
041 * are mostly for users who wants to extend the existing <em>MOEA</em> classes.
042 *
043 * @author <a href="mailto:franz.wilhelmstoetter@gmail.com">Franz Wilhelmstötter</a>
044 * @version 4.1
045 * @since 4.1
046 */
047 public final class Pareto {
048
049 private Pareto() {
050 }
051
052 /* *************************************************************************
053 * Crowding distance methods.
054 * ************************************************************************/
055
056 /**
057 * The crowding distance value of a solution provides an estimate of the
058 * density of solutions surrounding that solution. The <em>crowding
059 * distance</em> value of a particular solution is the average distance of
060 * its two neighboring solutions.
061 *
062 * @apiNote
063 * Calculating the crowding distance has a time complexity of
064 * {@code O(d*n*log(n))}, where {@code d} is the number of dimensions and
065 * {@code n} the {@code set} size.
066 *
067 * @see #crowdingDistance(Seq, ElementComparator, ElementDistance, ToIntFunction)
068 *
069 * @param set the point set used for calculating the <em>crowding distance</em>
070 * @param <T> the vector type
071 * @return the crowded distances of the {@code set} points
072 * @throws NullPointerException if the input {@code set} is {@code null}
073 * @throws IllegalArgumentException if {@code set.get(0).length() < 2}
074 */
075 public static <T> double[] crowdingDistance(final Seq<? extends Vec<T>> set) {
076 return crowdingDistance(
077 set,
078 Vec::compare,
079 Vec::distance,
080 Vec::length
081 );
082 }
083
084 /**
085 * The crowding distance value of a solution provides an estimate of the
086 * density of solutions surrounding that solution. The <em>crowding
087 * distance</em> value of a particular solution is the average distance of
088 * its two neighboring solutions.
089 *
090 * @apiNote
091 * Calculating the crowding distance has a time complexity of
092 * {@code O(d*n*log(n))}, where {@code d} is the number of dimensions and
093 * {@code n} the {@code set} size.
094 *
095 * @see #crowdingDistance(Seq)
096 *
097 * @param set the point set used for calculating the <em>crowding distance</em>
098 * @param comparator the comparator which defines the (total) order of the
099 * vector elements of {@code T}
100 * @param distance the distance of two vector elements
101 * @param dimension the dimension of vector type {@code T}
102 * @param <T> the vector type
103 * @return the crowded distances of the {@code set} points
104 * @throws NullPointerException if one of the arguments is {@code null}
105 */
106 public static <T> double[] crowdingDistance(
107 final Seq<? extends T> set,
108 final ElementComparator<? super T> comparator,
109 final ElementDistance<? super T> distance,
110 final ToIntFunction<? super T> dimension
111 ) {
112 requireNonNull(set);
113 requireNonNull(comparator);
114 requireNonNull(distance);
115
116 final double[] result = new double[set.size()];
117 if (set.size() < 3) {
118 Arrays.fill(result, POSITIVE_INFINITY);
119 } else {
120 for (int m = 0, d = dimension.applyAsInt(set.get(0)); m < d; ++m) {
121 final int[] idx = ProxySorter.sort(
122 set,
123 comparator.ofIndex(m).reversed()
124 );
125
126 result[idx[0]] = POSITIVE_INFINITY;
127 result[idx[set.size() - 1]] = POSITIVE_INFINITY;
128
129 final T max = set.get(idx[0]);
130 final T min = set.get(idx[set.size() - 1]);
131 final double dm = distance.distance(max, min, m);
132
133 if (Double.compare(dm, 0) > 0) {
134 for (int i = 1, n = set.size() - 1; i < n; ++i) {
135 final double dist = distance.distance(
136 set.get(idx[i - 1]),
137 set.get(idx[i + 1]),
138 m
139 );
140
141 result[idx[i]] += dist/dm;
142 }
143 }
144 }
145 }
146
147 return result;
148 }
149
150 /* *************************************************************************
151 * Pareto ranks methods.
152 * ************************************************************************/
153
154 /**
155 * Calculates the <em>non-domination</em> rank of the given input {@code set},
156 * using the <em>natural</em> order of the elements as <em>dominance</em>
157 * measure.
158 *
159 * @apiNote
160 * Calculating the rank has a time complexity of {@code O(n^2}, where
161 * {@code n} the {@code set} size.
162 *
163 * <p>
164 * <b>Reference:</b><em>
165 * Kalyanmoy Deb, Associate Member, IEEE, Amrit Pratap,
166 * Sameer Agarwal, and T. Meyarivan.
167 * A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II,
168 * IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 6, NO. 2,
169 * APRIL 2002.</em>
170 *
171 * @param set the input set
172 * @param <T> the element type
173 * @return the <em>non-domination</em> rank of the given input {@code set}
174 */
175 public static <T> int[] rank(final Seq<? extends Vec<T>> set) {
176 return rank(set, Vec::dominance);
177 }
178
179 /**
180 * Calculates the <em>non-domination</em> rank of the given input {@code set},
181 * using the given {@code dominance} comparator.
182 *
183 * @apiNote
184 * Calculating the rank has a time and space complexity of {@code O(n^2},
185 * where {@code n} the {@code set} size.
186 *
187 * <p>
188 * <b>Reference:</b><em>
189 * Kalyanmoy Deb, Associate Member, IEEE, Amrit Pratap,
190 * Sameer Agarwal, and T. Meyarivan.
191 * A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II,
192 * IEEE TRANSACTIONS ON EVOLUTIONARY COMPUTATION, VOL. 6, NO. 2,
193 * APRIL 2002.</em>
194 *
195 * @param set the input set
196 * @param dominance the dominance comparator used
197 * @param <T> the element type
198 * @return the <em>non-domination</em> rank of the given input {@code set}
199 */
200 public static <T> int[] rank(
201 final Seq<? extends T> set,
202 final Comparator<? super T> dominance
203 ) {
204 // Pre-compute the dominance relations.
205 final int[][] d = new int[set.size()][set.size()];
206 for (int i = 0; i < set.size(); ++i) {
207 for (int j = i + 1; j < set.size(); ++j) {
208 d[i][j] = dominance.compare(set.get(i), set.get(j));
209 d[j][i] = -d[i][j];
210 }
211 }
212
213 // Compute for each element p the element q that it dominates and the
214 // number of times it is dominated. Using the names as defined in the
215 // referenced paper.
216 final int[] nq = new int[set.size()];
217 final List<IntList> fronts = new ArrayList<>();
218 IntList Fi = new IntList();
219
220 for (int p = 0; p < set.size(); ++p) {
221 final IntList Sp = new IntList();
222 int np = 0;
223
224 for (int q = 0; q < set.size(); ++q) {
225 if (p != q) {
226 // If p dominates q, add q to the set of solutions
227 // dominated by p.
228 if (d[p][q] > 0) {
229 Sp.add(q);
230
231 // Increment the domination counter of p.
232 } else if (d[q][p] > 0) {
233 np += 1;
234 }
235 }
236 }
237
238 // p belongs to the first front.
239 if (np == 0) {
240 Fi.add(p);
241 }
242
243 fronts.add(Sp);
244 nq[p] = np;
245 }
246
247 // Initialize the front counter.
248 int i = 0;
249 final int[] ranks = new int[set.size()];
250 while (!Fi.isEmpty()) {
251 // Used to store the members of the next front.
252 final IntList Q = new IntList();
253
254 for (int p = 0; p < Fi.size(); ++p) {
255 final int fi = Fi.get(p);
256 ranks[fi] = i;
257
258 // Update the dominated counts as compute next front.
259 for (int k = 0, n = fronts.get(fi).size(); k < n; ++k) {
260 final int q = fronts.get(fi).get(k);
261 nq[q] -= 1;
262
263 // q belongs to the next front.
264 if (nq[q] == 0) {
265 Q.add(q);
266 }
267 }
268 }
269
270 ++i;
271 Fi = Q;
272 }
273
274 return ranks;
275 }
276
277 /* *************************************************************************
278 * 'front'
279 * ************************************************************************/
280
281 /**
282 * Return the elements, from the given input {@code set}, which are part of
283 * the pareto front. The {@link Comparable} interface defines the dominance
284 * measure of the elements, used for calculating the pareto front.
285 * <p>
286 * <b>Reference:</b><em>
287 * E. Zitzler and L. Thiele.
288 * Multiobjective Evolutionary Algorithms: A Comparative Case Study
289 * and the Strength Pareto Approach,
290 * IEEE Transactions on Evolutionary Computation, vol. 3, no. 4,
291 * pp. 257-271, 1999.</em>
292 *
293 * @param set the input set
294 * @param <T> the element type
295 * @return the elements which are part of the pareto set
296 * @throws NullPointerException if one of the arguments is {@code null}
297 */
298 public static <T> ISeq<Vec<T>> front(final Seq<? extends Vec<T>> set) {
299 return front(set, Vec::dominance);
300 }
301
302 /**
303 * Return the elements, from the given input {@code set}, which are part of
304 * the pareto front.
305 * <p>
306 * <b>Reference:</b><em>
307 * E. Zitzler and L. Thiele.
308 * Multiobjective Evolutionary Algorithms: A Comparative Case Study
309 * and the Strength Pareto Approach,
310 * IEEE Transactions on Evolutionary Computation, vol. 3, no. 4,
311 * pp. 257-271, 1999.</em>
312 *
313 * @param set the input set
314 * @param dominance the dominance comparator used
315 * @param <T> the element type
316 * @return the elements which are part of the pareto set
317 * @throws NullPointerException if one of the arguments is {@code null}
318 */
319 public static <T> ISeq<T> front(
320 final Seq<? extends T> set,
321 final Comparator<? super T> dominance
322 ) {
323 final MSeq<T> front = MSeq.of(set);
324
325 int n = front.size();
326 int i = 0;
327 while (i < n) {
328 int j = i + 1;
329 while (j < n) {
330 if (dominance.compare(front.get(i), front.get(j)) > 0) {
331 --n;
332 front.swap(j, n);
333 } else if (dominance.compare(front.get(j), front.get(i)) > 0) {
334 --n;
335 front.swap(i, n);
336 --i;
337 break;
338 } else {
339 ++j;
340 }
341 }
342 ++i;
343 }
344
345 return front.subSeq(0, n).copy().toISeq();
346 }
347
348
349 /* *************************************************************************
350 * Common 'dominance' methods.
351 * ************************************************************************/
352
353 /**
354 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
355 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
356 *
357 * @see Vec#dominance(Comparable[], Comparable[])
358 *
359 * @param u the first vector
360 * @param v the second vector
361 * @param <C> the element type of vector <b>u</b> and <b>v</b>
362 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
363 * <b>u</b> and {@code 0} otherwise
364 * @throws NullPointerException if one of the arguments is {@code null}
365 * @throws IllegalArgumentException if {@code u.length != v.length}
366 */
367 public static <C extends Comparable<? super C>> int
368 dominance(final C[] u, final C[] v) {
369 return dominance(u, v, Comparator.naturalOrder());
370 }
371
372 /**
373 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
374 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
375 *
376 * @see Vec#dominance(Object[], Object[], Comparator)
377 *
378 * @param u the first vector
379 * @param v the second vector
380 * @param comparator the element comparator which is used for calculating
381 * the dominance
382 * @param <T> the element type of vector <b>u</b> and <b>v</b>
383 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
384 * <b>u</b> and {@code 0} otherwise
385 * @throws NullPointerException if one of the arguments is {@code null}
386 * @throws IllegalArgumentException if {@code u.length != v.length}
387 */
388 public static <T> int
389 dominance(final T[] u, final T[] v, final Comparator<? super T> comparator) {
390 requireNonNull(comparator);
391 checkLength(u.length, v.length);
392
393 return dominance(
394 u, v, u.length, (a, b, i) -> comparator.compare(a[i], b[i])
395 );
396 }
397
398 /**
399 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
400 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
401 *
402 * @see Vec#dominance(int[], int[])
403 *
404 * @param u the first vector
405 * @param v the second vector
406 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
407 * <b>u</b> and {@code 0} otherwise
408 * @throws NullPointerException if one of the arguments is {@code null}
409 * @throws IllegalArgumentException if {@code u.length != v.length}
410 */
411 public static int dominance(final int[] u, final int[] v) {
412 checkLength(u.length, v.length);
413
414 return dominance(
415 u, v, u.length, (a, b, i) -> Integer.compare(a[i], b[i])
416 );
417 }
418
419 /**
420 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
421 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
422 *
423 * @see Vec#dominance(long[], long[])
424 *
425 * @param u the first vector
426 * @param v the second vector
427 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
428 * <b>u</b> and {@code 0} otherwise
429 * @throws NullPointerException if one of the arguments is {@code null}
430 * @throws IllegalArgumentException if {@code u.length != v.length}
431 */
432 public static int dominance(final long[] u, final long[] v) {
433 checkLength(u.length, v.length);
434
435 return dominance(
436 u, v, u.length, (a, b, i) -> Long.compare(a[i], b[i])
437 );
438 }
439
440 /**
441 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
442 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
443 *
444 * @see Vec#dominance(double[], double[])
445 *
446 * @param u the first vector
447 * @param v the second vector
448 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
449 * <b>u</b> and {@code 0} otherwise
450 * @throws NullPointerException if one of the arguments is {@code null}
451 * @throws IllegalArgumentException if {@code u.length != v.length}
452 */
453 public static int dominance(final double[] u, final double[] v) {
454 checkLength(u.length, v.length);
455
456 return dominance(
457 u, v, u.length, (a, b, i) -> Double.compare(a[i], b[i])
458 );
459 }
460
461 private static void checkLength(final int i, final int j) {
462 if (i != j) {
463 throw new IllegalArgumentException(format(
464 "Length are not equals: %d != %d.", i, j
465 ));
466 }
467 }
468
469 private static <T> int dominance(
470 final T u,
471 final T v,
472 final int dimension,
473 final ElementComparator<? super T> comparator
474 ) {
475 boolean udominated = false;
476 boolean vdominated = false;
477
478 for (int i = 0; i < dimension; ++i) {
479 final int cmp = comparator.compare(u, v, i);
480
481 if (cmp > 0) {
482 udominated = true;
483 if (vdominated) {
484 return 0;
485 }
486 } else if (cmp < 0) {
487 vdominated = true;
488 if (udominated) {
489 return 0;
490 }
491 }
492 }
493
494 if (udominated == vdominated) {
495 return 0;
496 } else if (udominated) {
497 return 1;
498 } else {
499 return -1;
500 }
501 }
502
503 }
|