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